沿飞行弹道的扰动引力逼近方法

周欢， 丁智坚， 郑伟

(1.中国工程物理研究院 总体工程研究所， 四川 绵阳 621999； 2.中国空气动力研究与发展中心， 四川 绵阳 621000；3.国防科技大学 航天科学与工程学院， 湖南 长沙 410073)

0 引言

地球外部空间扰动引力作用于飞行过程中的导弹，影响其导航解算和制导计算，进而影响落点精度，故需要在发射前构建引力模型。现有点质量法、球谐函数、Stokes方法等空间扰动引力计算方法难以满足诸元计算及飞行过程中扰动引力赋值的存储量和计算速度要求，由此需要建立快速逼近理论。快速逼近理论的核心是建立形式简单、计算效率高、存储量小和适应范围广的数值函数，既能精确地反映扰动引力空间特征，又能满足弹道计算要求。其中涉及到空间结构剖分、重构函数确定、效能损失控制以及重构精度与速度优化抉择等问题。

国外，点质量法[1]曾应用于民兵Ⅲ导弹的重力场计算中。为提高解算速度，一些学者探讨了谱方法[2-5]和数值逼近方法在扰动引力求解中的应用。Junkins[6]首次提出了重力位的有限元表达，并与Engels等[7]探讨了其在惯性导航计算中的应用。国内，广义延拓逼近[8]、内插[9-10]、三次等距B样条函数[11]等方法被用于求解主动段扰动引力，球谐函数换极法[12]被用于求解被动段扰动引力。然而，上述方法或难以满足弹上存储量和计算效率的要求，或因无法兼顾数据的空间趋势信息而提高逼近精度。

针对以上问题，本文提出了一种基于网函数逼近理论的沿飞行弹道扰动引力快速重构方法，推导了该逼近方法的余项误差，分析了影响赋值精度的因素，验证了该方法在各种情况下的适应性。

1 沿飞行弹道的空域剖分

基于有限元思想，按如下思路构建扰动引力快速计算模型：1)确定不考虑扰动引力的参考弹道；2)形成以参考弹道为中心的飞行管道，进行空域剖分并确定各节点位置；3)通过高精度引力计算方法进行节点扰动引力赋值；4)判断计算点所在单元及其与该单元各节点的相对位置关系；5)根据一定的逼近方法，由节点扰动引力计算当前点扰动引力。扰动引力对落点精度的影响随飞行时间累积，拟采用“漏斗形”管道以减少单元数和防止实际弹道超出管道。

以主动段飞行管道的构建为例。利用主动段弹道相对平直的特性，在发射坐标系内生成管道，进行空域剖分：1)在标准弹道上依次选定基准点d0，d1，d2，d3，…，di，相邻两点间沿y轴方向的距离为δy1,δy2,δy3,…，δyi；2)以di为几何中心，在平行于Oxz的平面内生成边长分别为δxi和δzi的长方形；3)令δxi<δxi+1，δyi<δyi+1，δzi<δzi+1，依次连接各长方形顶点ki，构成主动段“漏斗形”飞行管道。采用类似方法在轨道坐标系内构建被动段飞行管道。

2 扰动引力快速逼近方法

快速逼近方法是高精度扰动引力快速赋值的主要方法，扰动引力的快速逼近主要涉及两个问题：一是判断计算点所在单元；二是单元内部的逼近计算。前者通过比较计算点与di的相对位置关系获得，后者基于网函数逼近理论实现。

由Cook[13]提出的网函数逼近理论是一种多方向拟合一次误差调整的多变元函数逼近方法，其本质上是单变元函数的Lagrange插值算法的极自然推广。其插值函数类是简单的高阶偏微分方程边值问题的解，具有Coons型结构以及明显的统计特征，算法简单。该方法的基本思想是通过n维欧式空间中网格边界上值来近似计算网格中任意点的值，与单元内部的扰动引力逼近问题契合。

令L(x)、L(y)、L(z)分别为关于x、y、z的一次Lagrange插值算符，其插值基函数为

(1)

记主动段某六面体单元的8个顶点为ki,1、ki,2、ki,3、ki,4、ki+1,1、ki+1,2、ki+1,3、ki+1,4，其扰动引力值分别为gi,1、gi,2、gi,3、gi,4、gi+1,1、gi+1,2、gi+1,3、gi+1,4；记12条棱为L0～L11，称其为计算单元上的1-网，定义在1-网上的扰动引力值为fi(x,y,z)(i=0,1,…,11)。根据网函数逼近理论，已知1-网函数，可由三维1-网插值方法求取单元内任意一点的值。令l(3)为三维1-网插值算符，则

l(3)=L(x)L(y)+L(y)L(z)+L(z)L(x)-
2L(x)L(y)L(z).

(2)

据此，可求得单元内任意一点A(x,y,z)的值F(x,y,z)，

F(x,y,z)=F1(x,y,z)+F2(x,y,z)+
F3(x,y,z)-F4(x,y,z)，

(3)

式中：

(4)

节点扰动引力可由传统扰动引力计算方法获得，而fi(x,y,z)(i=0,1,…,7)则通过对每条棱所对应的两节点插值获得。对于沿弹道方向的棱fi(x,y,z)(i=8,9,10,11)，为保证飞行管道的平滑，纳入相邻网格数据的空间趋势信息；同时为避免因多点插值引起的龙格现象，采用加权3点插值的方式求取。以f9(x,y,z)的求解为例，

(5)

通过相同方法求取其他3条沿弹道方向的棱上的扰动引力值。至此，已完成了主动段单元内部任意一点扰动引力的求解，通过类似方法可完成被动段扰动引力的快速逼近。

3 逼近精度分析

3.1 插值余项及误差估计

称R(3)[F]=F(x,y,z)-l(3)[fi(x,y,z)]为函数F(x,y,z)的l(3)插值余项，记γ(x,y,z)=R(3)[F]。当F(x,y,z)有直到6阶的连续偏导数，则可推导γ(x,y,z)的估计式，

(6)

式中：

(7)

(8)

η1、ζ1为Oyz平面上的给定点坐标，ζ2、ξ2为Oxz平面上的给定点坐标，ξ3、η3为Oxy平面上的给定点坐标，ξ0、η0、ζ0为Oxyz空间坐标系下的给定点坐标。

由(6)式可知，γ(x,y,z)取决于计算单元大小以及表征F(x,y,z)变化率的高阶偏导数。因此，单元越小，数据变化越平缓，计算点离节点越近，则逼近程度越高。

3.2 逼近精度影响因素分析

本节通过数值仿真，对3.1节的结论进行验证。以分层点质量法进行节点扰动引力赋值。针对某一固定空域进行等大小的网格剖分，以网函数方法和分层点质量法分别计算空域内1 000个随机点的扰动引力，二者之差即为逼近误差。

3.2.1 改变网格大小

计算原点坐标(λ,φ,H)=(110°,40°,0 km)，其中λ为经度，φ为纬度，H为高程；计算范围(Δλ,Δφ,ΔH)=(±1.5°,±1.5°,300 km)。选取3种大小的网格：

1)尺寸1网格(dλ,dφ,dH)=(0.25°,0.25°,5 km)；

2)尺寸2网格(dλ,dφ,dH)=(0.5°,0.5°,10 km)；

3)尺寸3网格(dλ,dφ,dH)=(0.75°,0.75°,15 km)。

扰动引力三分量逼近误差的统计结果见表1.

表1 不同网格大小的逼近误差统计

3.2.2 改变网格形状

计算原点坐标与3.2.1节相同，选取体积大致相等、形状不同的4种网格：

1)形状1网格(dλ,dφ,dH)=(0.1°,0.1°,10 km)；

2)形状2网格(dλ,dφ,dH)=(0.2°,0.2°,2.5 km)；

3)形状3网格(dλ,dφ,dH)=(0.05°,0.05°,40 km)；

4)形状4网格(dλ,dφ,dH)=(0.2°,0.1°,20 km)。

4种网格形状的示意图见图1. 图1中，R为沿地球矢径方向，L为沿经度方向，B为沿纬度方向。扰动引力三分量逼近误差的统计结果见表2.

表2 不同网格形状的赋值误差统计

3.2.3 改变数据空间变化趋势

基于同一地区扰动引力的变化随高程增加而趋于平缓的结论，设计数值仿真算例。计算范围及计算原点经度和纬度同3.2.2节，计算原点高程分别为：1)H=0 km；2)H=20 km；3)H=40 km；4)H=60 km. 网格大小为(dλ,dφ,dH)=(0.5°,0.5°,10 km)。扰动引力三分量逼近误差的统计结果见表3.

表3 数据变化趋势对赋值误差统计结果的影响

由表2的计算结果可以看出，网格形状越趋近正方体，逼近误差越小。由表3的计算结果可以看出，随高程增加，扰动引力数据趋于平缓，逼近误差逐渐减小。对比表1和表2的结果可知，形状2～形状4的网格远小于尺寸1～尺寸3的网格，但逼近精度并没有得到明显提高，甚至低于尺寸1的赋值精度。由此可知，网格形状比网格大小对逼近精度具有更大的影响。由表1～表3数值仿真得到的结论与3.1节中理论推导结论一致。

4 快速逼近方法的适应性分析

由于地球外部空间扰动引力场受地球表面形状及地球内部物质分布的影响，通过将快速逼近方法应用于不同发射区域、不同射向和不同射程的弹道导弹全程弹道计算中检验方法的适应性。适应性判断准则如下：

1)扰动引力计算模型适应于任意弹道，全程弹道积分不出现奇异；

2)实际飞行弹道不超出预先构建的“漏斗形”弹道管道逼近空间；

3)根据弹道导弹射前引力模型构建的一般要求，扰动引力单向分量的逼近误差均值应不大于1 mgal；

4)根据弹道导弹对落点精度的一般要求，由扰动引力二次逼近误差导致的落点偏差应不超过100 m.

4.1 单弹头全程扰动引力快速赋值结果分析

设置主动段起始单元边长为0.5 km，单元边长增幅为2 km. 被动段起始单元边长为δr1=δξ1=200 km，δf1=2°，单元边长增幅为20 km和0.5°. 节点扰动引力及近似真值均采用1 080阶球谐函数进行计算。

4.1.1 不同发射区域

针对射程为12 000 km，发射方位角为90°的弹道开展研究。发射点分别为：1)平原λ=115°，φ=35°，H=0 km；2)丘陵λ=110°，φ=27°，H=1 km；3)特大山区λ=80°，φ=42°，H=3 km. 赋值结果见表4.

表4 不同发射区域的赋值误差统计结果

4.1.2 不同发射方位角

针对射程为12 000 km，发射点位于特大山区的弹道开展研究。发射方位角α分别为：1)正北α=0°；2)正东α=90°；3)正南α=180°；4)正西α=270°. 赋值结果见表5.

表5 不同发射方位角的赋值误差统计结果

4.1.3 不同射程

针对发射点位于特大山区，发射方位角α=90°，射程分别为5 000 km、8 000 km和12 000 km的弹道开展研究。赋值结果见表6.

表6 不同射程弹道的赋值误差统计结果

4.1节的仿真结果表明，本文建立的扰动引力快速逼近方法能够适应各种条件下的弹道计算，全程弹道积分无奇异且不存在实际飞行弹道超出预设逼近空间的情况。表4～表6的统计结果表明，在不同发射区域、不同射向和不同射程的弹道导弹全程弹道计算中，扰动引力三向分量的逼近误差均值控制在10-2mgal量级，满足单向分量逼近误差均值不超过1 mgal的要求。

4.2 双弹头全程扰动引力快速赋值结果分析

本节将扰动引力快速逼近方法应用于多头分导弹道导弹，以考察方法的适应性。假定导弹射向为正北，发射点的经度、纬度和高程分别为111°、39°和1 500 m. 假设分导弹头的射程小于主弹头，且主弹头落点和分导弹头落点纵程相差约1 000 km，横程相差约160 km. 节点扰动引力采用1 080阶球谐函数法赋值，并以球谐函数法计算结果作为计算参考值。考虑到分导过程发生在主动段结束以后，分导弹头与主弹头的弹道仅在被动段不同，因此只针对被动段扰动引力赋值情况进行分析。被动段空域剖分参数设置如下：1)起始网格大小δz0=δr0=200 km，δf0=2.0°；2)相邻网格边长增幅δzi+1-δzi=δri+1-δri=20 km，δfi+1-δfi=0.5°. 经仿真可以得到，主弹头和分导弹头的被动段弹道及其弹道管道，如图2所示。

上述空域剖分方案对应的网格总数为25个，节点数为88个，总存储量为528个数据，可以满足弹上数据存储量要求。根据上述单元划分，将快速赋值方法的计算结果与球谐函数方法计算结果进行比较，可以得到主弹头及分导弹头被动段扰动引力逼近误差，如图3和图4所示。

由图3和图4可知，在被动段飞行的大部分时段，主弹头和分导弹头扰动引力逼近误差均可控制在1 mgal以内，具有较高的赋值精度。仅当飞行至被动段末端，由于漏斗形弹道管道形成的逼近空间逐渐增大，导致赋值精度下降，但逼近误差仍控制在10 mgal以内。随着飞行时间增加，扰动引力赋值误差对落点精度的影响逐渐降低，上述赋值精度满足应用需求。

在主弹头和分导弹头的落点偏差计算中，分别基于球谐函数方法和快速赋值方法求解扰动引力项。表7给出了两种扰动计算方法对应的弹头落点纵向偏差、横向偏差和总偏差。以球谐函数的计算结果为参考值，视两种方法计算结果之差为快速赋值方法的逼近误差对应的落点偏差。结果表明，快速赋值方法逼近误差对应的落点偏差可控制在米级，具有较高精度。

表7 扰动引力计算误差对应的落点偏差

4.3 逼近误差对落点精度的影响分析

针对发射点为λ=110°，φ=40°，H=1.5 km，射程分别为5 000 km、8 000 km和12 000 km，发射方位角为αi=-180°+i×5°(i=0,1,2,…,72)的弹道开展研究。弹道积分中的扰动引力项分别采用快速逼近方法和1 080阶球谐函数方法计算得到，视两种弹道计算结果对应的落点偏差之差为赋值误差引起的落点偏差，计算结果如图5所示。结果表明，各种情况下赋值误差引起的落点偏差均不超过8 m，满足实际应用要求。

综合4.1节～4.3节的仿真结果，认为本文建立的扰动引力快速逼近方法满足逼近算法适应性判断准则，能够适应于不同情况下的弹道计算。

5 讨论

5.1 计算时间及存储量

本文提出的扰动引力快速逼近方法主要包括引力模型构建及飞行过程中的赋值计算，其中前者主要包括空域剖分、节点坐标确定及节点扰动引力赋值，后者主要包括当前计算点所在单元判断及单元内部的逼近计算。存储数据主要包括各节点位置三分量及各节点扰动引力三分量。根据第4.1节的网格划分条件，表8列出了3种射程弹道的模型构建速度及存储量分析结果，表9列出了不同扰动引力计算方法下射程为12 000 km弹道的计算时间。仿真计算所使用的计算机基本配置为：Celeron(R)CPU 2.53 GHz，内存大小为512 MB. 软件环境为Window XP操作系统，计算程序基于VC++6.0开发。结果表明，本文提出的方法在计算时间和存储量方面具有很大的优势。

表8 重构速度及存储量分析

表9 不同计算方法对应的计算时间

5.2 赋值精度及适应性

针对某射程为12 000 km，射向为正北方向的弹道，在不考虑节点扰动引力赋值误差情况下，基于网函数方法和型函数方法分别计算沿不同高度弹道的扰动引力。与第4节相同，仍以1 080阶球谐函数方法的扰动引力计算结果作为近似真值，将网函数与型函数的计算结果分别与之比较，可以得到各自的赋值误差均方差，如表10所示。

表10 网函数和型函数赋值误差统计结果

由表10可以看出，网函数的赋值精度略高。从理论上讲，网函数与型函数都是Lagrange插值的推广，因此具有量级相当的逼近精度。但与型函数等内插计算方法相比，基于网函数逼近理论的赋值方法除利用当前单元节点的离散数据外，还考虑了相邻单元节点数据的空间趋势信息，根据引力异常变化趋势构造最佳的1-网逼近曲线，进而以逼近曲线为边界曲线来调控整个曲面的形状和趋向，构建分块表征模型。这一算法相当于在离散的节点数据基础上追加了一定的边界条件，因此提高了精度。由于采用的数据都是已有节点的数据，因此无需增加额外的存储量。

该方法在不同发射区域、不同射向和不同射程弹道中的应用结果表明，其赋值误差在10-2mgal量级，由此引起的落点偏差在8 m以内，具有精度高、计算速度快、存储量小和适应范围广等特点。

6 结论

本文基于标准弹道构建飞行管道，实现对沿飞行弹道附近空域的有限元剖分，首次将网函数逼近理论应用于扰动引力快速赋值计算，提出了一种沿飞行弹道的扰动引力快速计算方法。在理论推导该方法的插值余项误差基础上，基于某一区域的扰动引力测量数据，数值模拟了不同因素对赋值精度的影响。以1 080阶球谐函数构建的地球引力模型为近似真值，分析了提出的快速计算方法在不同情况下弹道导弹全程弹道计算中的应用情况。结果表明，该方法产生的赋值误差约为10-2mgal量级，由此引起的落点偏差在8 m以内，全程弹道计算时间远小于其他方法。因此，本文提出的方法能够实现沿飞行弹道任意点扰动引力三分量的快速赋值，其赋值精度、计算速度和存储量均满足弹道计算要求，具有广泛的适应性。