CN104751012A - 沿飞行弹道的扰动引力快速逼近方法 - Google Patents

Abstract

本发明以弹道导弹为研究对象，针对导弹飞行过程中的扰动引力快速赋值问题，首次提出了一种沿飞行弹道的全程扰动引力快速赋值方法。首先，根据发射任务生成标准弹道；其次，生成以标准导弹为基准的飞行管道，进行空域剖分；之后，通过点质量法或高阶球谐函数法对节点扰动引力进行赋值；最终，在计算单元内部基于网函数逼近理论计算实际弹道上任意一点的扰动引力值，完成沿飞行弹道的扰动引力全程快速赋值计算。该方法能够实现沿任意飞行弹道的扰动引力快速赋值，其赋值精度满足弹道计算要求，赋值速度和数据存储量等指标远优于现有方法，解决了扰动引力弹上实时赋值从无到有的问题。

Description

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

技术领域

本发明属于飞行器动力学建模领域，特别涉及沿飞行弹道的扰动引力快速逼近方法。

背景技术

持续提升射前机动快速发射性能及落点精度是我国新一代弹道导弹发展的必然趋势。地球扰动引力模型的构建是射前准备中的必要环节，其构建速度直接影响发射快速性。在影响弹道导弹落点精度的因素中，飞行过程中的扰动引力是除地球扁率和垂线偏差以外产生制导方法误差的主要来源。因此，建立扰动引力场模型快速重构方法及飞行过程中的高精度扰动引力赋值方法具有重大的军事意义和工程价值。

当前，地球物理和大地测量领域的学者提出了若干地球外部空间扰动引力的求解方法，但由于这些方法的落脚点并不在于其在弹道计算中的应用，因此通常未对计算规模和计算速度有所要求，其算法的复杂性使之无法满足模型快速构建及飞行过程中实时赋值的要求。目前唯一有资料表明用于弹道导弹扰动引力计算的是点质量法，该方法欲求解空间任意一点的扰动引力，需要基于该区域点质量数据以及通过对大量扰动质量的引力进行求和，存储量大且计算规模大，因此仍然无法满足新形势下的快速性和高精度要求。

因此，亟待建立一种面向快速机动发射应用的扰动引力计算方法，该方法应具有“计算速度快、适应区域广、修正模型精”的特征。该问题的难点在于，在保证必要赋值精度的前提下满足弹载计算机存储量及计算能力的要求。其核心是寻求一种易于实现的模型重构策略，建立一种形式简单、计算效率高和逼近程度高的数值函数，使之既能精确反映扰动引力的空间特征，又能满足弹道计算要求。其间涉及的问题主要包括空间剖分建模、重构函数确定、效能损失控制以及重构精度与速度优化抉择等。

针对上述问题，提出一种沿飞行弹道的扰动引力快速逼近方法，该方法首先需在飞行任务要求下生成一条标准弹道，其次需生成沿标准弹道的飞行管道，在管道内部实现空域剖分，最后基于网函数逼近理论实现实际弹道上任意点的扰动引力快速赋值。建立的方法力争填补面向弹道导弹应用的扰动引力实时快速赋值方法的空白，实现我国扰动引力场模型从“事先构建”向“射前动态构建”的跨越，为新一代弹道导弹实现既定区域内快速机动发射和高命中精度提供技术基础和方法支撑。

发明内容

本发明针对沿飞行弹道的扰动引力快速逼近问题，提出一种基于沿标准弹道构建“漏斗形”飞行管道和网函数逼近理论的快速赋值方法。

该方法综合利用标准弹道生成、飞行管道建立、计算点位置确定和网函数逼近理论实现沿飞行弹道的扰动引力快速逼近。方法借鉴有限元的基本思想，首先基于标准弹道生成求解区域，之后将求解区域离散，根据一定的准则建立起单元的泛函表达式，由此得到以单元节点泛函为未知数的线性方程组，通过解方程组得到问题的数值解。基于该方法的扰动引力快速赋值基本思路如下：首先，根据发射任务要求确定一条不考虑扰动引力的标准弹道；其次，基于标准弹道生成“漏斗形”飞行管道，在飞行管道内完成空域剖分并确定节点位置；再次，采用点质量法、球谐函数法等方法计算得到的扰动引力对节点赋值；最后，在诸元或弹上导航计算时，根据导弹位置判断其所在的单元，并根据当前位置与所在单元各节点相对位置关系以及各节点的扰动引力值，以一定的逼近准则快速计算当前位置对应扰动引力值。

本发明的技术方案主要包括以下步骤：

第一步，标准弹道生成

根据发射任务要求，生成不考虑扰动引力的标准弹道。基于地球圆球假设，认为惯性空间中被动段弹道为椭圆弹道的一部分。

第二步，飞行管道生成

(1)主动段飞行管道生成

在发射坐标系内按如下方式生成主动段飞行管道：

①在标准弹道上依次选定基准点d₀,d₁,d₂,d₃...，相邻两点间沿y轴方向的距离为δy₁,δy₂,δy₃...；

②以d_i为几何中心，在垂直于y轴的平面内生成边长分别为δx_i和δz_i的四边形；

③令δx_i＜δx_i+1，δy_i＜δy_i+1，δz_i＜δz_i+1，依次连接各长方形顶点k_i，构成主动段“漏斗形”飞行管道。

(2)被动段飞行管道生成

在局部轨道坐标系内按如下方式生成被动段飞行管道：

①在标准弹道上选定基准点d₀,d₁,d₂,...，各点对应的真近点角分别为f₀,f₁,f₂,...，设δf_j＝f_j-f_j-1(j＝1,2,...)；

②以d_j为原点，建立局部轨道坐标系d_j-r_jβ_jξ_j，r_j轴沿地心矢方向，ξ_j轴与动量矩方向一致，β_j轴与r_j、ξ_j轴构成右手坐标系；

③在轨道坐标系内垂直于β_j轴的平面内生成边长分别为δr_j和δξ_j的四边形；

④令δr_j＜δr_j+1，δξ_j＜δξ_j+1，δf_j＜δf_j+1，依次连接各长方形顶点k_j，构成被动段“漏斗形”飞行管道。

第三步，节点坐标位置确定

(1)主动段节点坐标位置确定

根据沿飞行弹道的空域剖分方法，确定基准点d_i对应的四个节点在发射系中的坐标分别为k₁(x_i0+δx_i/2,y_i0,z_i0+δz_i/2)，k₂(x_i0+δx_i/2,y_i0,z_i0-δz_i/2)，k₃(x_i0-δx_i/2,y_i0,z_i0-δz_i/2)，k₄(x_i0-δx_i/2,y_i0,z_i0+δz_i/2)。

(2)被动段节点坐标位置确定

根据沿飞行弹道的空域剖分方法，确定基准点d_j对应的四个节点在局部轨道坐标系中的坐标分别为k₁(δr_j/2,0,δξ_j/2)，k₂(δr_j/2,0,-δξ_j/2)，k₃(-δr_j/2,0,-δξ_j/2)，k₄(-δr_j/2,0,δξ_j/2)。

设被动段实际弹道上某点A在地心惯性坐标系中坐标为(x′,y′,z′)，在局部轨道坐标系中的坐标为(x,y,z)，根据如下坐标转换关系，求得节点在地心惯性系中的坐标：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=M_3(\mathrm{\ω}+f)\·M_1\left(i\right)\·M_3\left(\mathrm{\Ω}\right)\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]-\left[\begin{array}{c}r\\ 0\\ 0\end{array}\right]---\left(1\right)$

其中，r为地心距r，f为真近点角，ω为近地点角距，Ω为升交点角距，i为轨道倾角。M_i(i＝1,2,3)为坐标转换矩阵，具体形式如下：

$\begin{array}{c}{M}_1\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]\\ M_2\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& 0& -\sin \mathrm{\α}\\ 0& 1& 0\\ \sin \mathrm{\α}& 0& \cos \mathrm{\α}\end{array}\right]\\ M_3\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]\end{array}---\left(2\right)$

第四步，节点扰动引力赋值

(1)主动段点质量赋值法

通过多层点质量求解算法求得N个质点的质量M_j、深度D_j和球坐标(j＝1,2,…,N)，可得扰动引力三分量的计算式：

其中，ρ_P＝R+H_P为计算点的地心距离，其它量的计算式为：R_j＝R-D_j。令为P点的球坐标，则

$r_{\mathrm{Pj}}={({\mathrm{\ρ}}_P^2+R_j^2-2{\mathrm{\ρ}}_P{R}_j\cos {\mathrm{\ψ}}_{\mathrm{Pj}})}^{1/2}---\left(4\right)$

(2)被动段球谐函数法赋值法

地球外部空间任意点P相对于旋转地球的地心距r、地心纬度经度λ为已知时，由球谐函数级数形式表示的地球引力位V为：

$V=\frac{\mathrm{\μ}}{r}[1+\underset{n=2}{\overset{s}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a_e}{r}\right)}^n\·({\stackrel{\‾}{C}}_{\mathrm{nm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin \mathrm{m\λ})\·{\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\sin \mathrm{\φ}\right)]---\left(8\right)$

令除主要位系数外的位系数都为零，得地球正常引力位

$\tilde{U}=\frac{\mathrm{\μ}}{r}[1+{\stackrel{\‾}{C}}_{20}{\left(\frac{a_e}{r}\right)}^2{\stackrel{\‾}{P}}_{20}\left(\sin \mathrm{\φ}\right)]---\left(9\right)$

真实引力位与正常引力位之差，即为扰动引力位T：

$T=V-\tilde{U}=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{s}{\mathrm{\Σ}}}{\left(\frac{a_e}{r}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}\·({\stackrel{\‾}{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\sin \mathrm{\φ}\right)]---\left(10\right)$

其中：

与T相应的引力即为扰动引力加速度即

$\mathrm{\δ}\stackrel{\→}{g}=\mathrm{grad}T---\left(12\right)$

则扰动引力加速度在天东北坐标系O_E-REN中的三个分量δg_R,δg_E,δg_N为：

第五步，当前计算点单元判断

(1)主动段单元判断方法

①假设第一个单元序号为1；

②确定实际弹道上某点A在发射坐标系中坐标为A(x^*,y^*,z^*)；

③比较y^*与基准点d_i的坐标y_i,0(i＝0,1,2,…)的大小。若y^*≥y_i,0且y^*≤y_i+1,0，则可确定A所在的单元序号为i+1。

(2)被动段单元判断算法

①假设第一个单元序号为1；

②确定实际弹道上某点A在发射坐标系中坐标为A(x^*,y^*,z^*)；

③求出A点对应的真近点角f^*；

④比较f^*与基准点d_j对应的真近点角y_j,0(j＝0,1,2,…)的大小。若f^*≥f_j,0且f^*≤f_j+1,0，则可确定A所在的单元序号为j+1。

第六步，单元内部逼近计算

基于网函数逼近理论进行单元内部的逼近计算。以主动段某计算单元为例，记8个顶点对应的扰动引力值分别为g_i,1、g_i,2、g_i,3、g_i,4、g_i+1,1、g_i+1,2、g_i+1,3、g_i+1,4；称12条棱为计算单元上的1-网，定义在其上的函数为1-网函数，记为f_i(x,y,z)(i＝0,1,…,11)。令L(x)、L(y)、L(z)分别为关于x、y、z的一次Lagrange插值算符，其插值基函数为

令为3维1-网函数插值算符，则

将作用于1-网函数，即可求得单元内任意一点A(x,y,z)的值F(x,y,z)，

垂直于飞行弹道的1-网函数由其对应的两节点插值求得，以f₆(x,y,z)的求取为例：

沿飞行弹道的1-网函数由其对应的两节点及两侧节点通过加权三点插值求得，以f₉(x,y,z)的求取为例：

$\begin{array}{c}{f}_9(x,y,z)=\frac{1}{2}L\left(y\right)\{g_{i-\mathrm{1,4}},g_{i,4},g_{i+\mathrm{1,4}}\}+\frac{1}{2}L\left(y\right)\{g_{i,4},g_{i+\mathrm{1,4}},g_{i+\mathrm{2,4}}\}\\ =\frac{1}{2}{\left[\begin{array}{c}\frac{(y-y_i)(y-y_{i+1})}{\mathrm{\δ}{y}_{i-1}(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i-1})}\\ -\frac{(y-y_{i-1})(y-y_{i+1})}{\mathrm{\δ}{y}_{i-1}\mathrm{\δ}{y}_i}\\ \frac{(y-y_i)(y-y_{i-1})}{\mathrm{\δ}{y}_i(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i-1})}\end{array}\right]}^T\left[\begin{array}{c}{g}_{i-\mathrm{1,4}}\\ g_{i,4}\\ g_{i+\mathrm{1,4}}\end{array}\right]+=\frac{1}{2}{\left[\begin{array}{c}\frac{(y-y_{i+1})(y-y_{i+2})}{\mathrm{\δ}{y}_i(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i+1})}\\ -\frac{(y-y_i)(y-y_{i+2})}{\mathrm{\δ}{y}_i\mathrm{\δ}{y}_{i+1}}\\ \frac{(y-y_i)(y-y_{i+1})}{\mathrm{\δ}{y}_{i+1}(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i+1})}\end{array}\right]}^T\end{array}\left[\begin{array}{c}{g}_{i,4}\\ g_{i+\mathrm{1,4}}\\ g_{i+\mathrm{2,4}}\end{array}\right]---\left(18\right)$

至此，经过上述六步可最终建立沿飞行弹道的扰动引力快速逼近方法，该方法具有“计算速度快、适应区域广、修正模型精”的特征，能满足快速机动发射及弹上实时计算的要求。与现有扰动引力赋值方法相比，本发明提出的方法具有以下优点：

1)首次提出一种面向弹道导弹应用的扰动引力快速赋值方法，其模型构建的快速性可使我国扰动引力场模型实现从“事先构建”向“射前动态构建”的跨越，存储量和计算规模的轻量化可突破飞行过程中扰动引力实时计算面临的重大工程挑战；

2)首次提出一种以标准弹道为基准的“漏斗形”飞行管道构建策略，在被动段引入局部轨道坐标系以简化计算，从而极大地缩小求解区域和减少节点数量，并能有效防止实际飞行弹道超出管道；

3)在计算单元内部采用网函数逼近理论进行任意点的扰动引力赋值，其插值函数类具有Coons型结构以及明显的统计特征，便于计算机实现；利用引力异常变化趋势构造最佳的1-网逼近曲线并以此为界调控整个曲面的形状和趋向，构建分块表征模型，因此具有较高的逼近精度；

4)将“延拓”思想引入1-网函数的求解算法，利用当前单元节点信息和相邻单元节点信息共同确定沿飞行弹道的1-网函数，能有效保证网格边界数据的连续性和平滑性；由于采用的数据都是已有节点的数据，因此无需增加额外的存储量和计算自由度；

5)该方法具有计算速度快、适应区域广、修正模型精、弹上存储量小的特征，可将扰动引力场模型重构时间控制在1分钟以内，将二次逼近误差在扰动引力赋值总误差中占比控制在20％以内，将逼近误差造成的导弹落点偏差控制在8m以内；重构方法适应于任意弹道，不出现奇点。

附图说明

图1为本发明主动段飞行管道示意图；

图2为本发明被动段飞行管道示意图；

图3为本发明网函数计算单元示意图；

图4为本发明沿飞行弹道的1-网函数计算单元示意图；

图5为本发明不同发射点弹道扰动引力快速赋值结果；

图6为本发明不同发射方位角弹道扰动引力快速赋值结果；

图7为本发明不同射程弹道扰动引力快速赋值结果；

图8为本发明二次逼近误差导致的落点偏差。

具体实施方式

下面结合具体实施例，对本发明作进一步的说明：

以某型号弹道导弹为仿真对象；主动段起始单元大小为δx₀＝δy₀＝δz₀＝0.5km，单元边长增幅为δx_i+1-δx_i＝δy_i+1-δy_i＝δz_i+1-δz_i＝2km；被动段起始单元大小为δr₀＝δξ₀＝200km，δf₀＝2°，单元边长增幅为δr_j+1-δr_j＝δξ_j+1-δξ_j＝20km，δf_j+1-δf_j＝0.5°。

仿真初始条件设置为：(1)射程为12000km，发射方位角为α＝90°，发射点分别选择为①平原地区(经度λ＝115°E，纬度高度H＝0km)；②丘陵地区(经度λ＝110°E，纬度高度H＝1km)；③特大山区(经度λ＝80°E，纬度高度H＝3km)；(2)射程为12000km，发射点位于特大山区(经度λ＝80°E，纬度高度H＝3km)，发射方位角分别选择为①正北(α＝0°)；②正东(α＝90°)；③正南(α＝180°)；④正西(α＝270°)；(3)发射点位于特大山区(经度λ＝80°E，纬度高度H＝3km)，发射方位角为α＝90°，射程分别选择为①12000km；②8000km；③5000km。

仿真计算机配置为：Intel(R)Core(TM)i5-3470CPU3.20GHz，内存为3.46GB。软件环境为Window XP操作系统，计算程序基于VC++6.0开发。

其具体步骤如下：

第一步，标准弹道生成

根据发射任务要求，生成不考虑扰动引力的标准弹道。基于地球圆球假设，认为惯性空间中被动段弹道为椭圆弹道的一部分。

第二步，飞行管道生成

(1)主动段飞行管道生成

在发射坐标系内按如下方式生成主动段飞行管道：

①在标准弹道上依次选定基准点d₀,d₁,d₂,d₃...，相邻两点间沿y轴方向的距离为δy₁,δy₂,δy₃...；

②以d_i为几何中心，在垂直于y轴的平面内生成边长分别为δx_i和δz_i的四边形；

③令δx_i＜δx_i+1，δy_i＜δy_i+1，δz_i＜δz_i+1，依次连接各长方形顶点k_i，构成主动段“漏斗形”飞行管道，见图1。

(2)被动段飞行管道生成

在局部轨道坐标系内按如下方式生成被动段飞行管道：

①在标准弹道上选定基准点d₀,d₁,d₂,...，各点对应的真近点角分别为f₀,f₁,f₂,...，设δf_j＝f_j-f_j-1(j＝1,2,...)；

②以d_j为原点，建立局部轨道坐标系d_j-r_jβ_jξ_j，r_j轴沿地心矢方向，ξ_j轴与动量矩方向一致，β_j轴与r_j、ξ_j轴构成右手坐标系；

③在轨道坐标系内垂直于β_j轴的平面内生成边长分别为δr_j和δξ_j的四边形；

④令δr_j＜δr_j+1，δξ_j＜δξ_j+1，δf_j＜δf_j+1，依次连接各长方形顶点k_j，构成被动段“漏斗形”飞行管道，见图2。

第三步，节点坐标位置确定

(1)主动段节点坐标位置确定

根据沿飞行弹道的空域剖分方法，确定基准点d_i对应的四个节点在发射系中的坐标分别为k₁(x_i0+δx_i/2,y_i0,z_i0+δz_i/2)，k₂(x_i0+δx_i/2,y_i0,z_i0-δz_i/2)，k₃(x_i0-δx_i/2,y_i0,z_i0-δz_i/2)，k₄(x_i0-δx_i/2,y_i0,z_i0+δz_i/2)。

(2)被动段节点坐标位置确定

根据沿飞行弹道的空域剖分方法，确定基准点d_j对应的四个节点在局部轨道坐标系中的坐标分别为k₁(δr_j/2,0,δξ_j/2)，k₂(δr_j/2,0,-δξ_j/2)，k₃(-δr_j/2,0,-δξ_j/2)，k₄(-δr_j/2,0,δξ_j/2)。

设被动段实际弹道上某点A在地心惯性坐标系中坐标为(x′,y′,z′)，在局部轨道坐标系中的坐标为(x,y,z)，根据如下坐标转换关系，求得节点在地心惯性系中的坐标：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=M_3(\mathrm{\ω}+f)\·M_1\left(i\right)\·M_3\left(\mathrm{\Ω}\right)\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]-\left[\begin{array}{c}r\\ 0\\ 0\end{array}\right]---\left(19\right)$

其中，r为地心距r，f为真近点角，ω为近地点角距，Ω为升交点角距，i为轨道倾角。M_i(i＝1,2,3)为坐标转换矩阵，具体形式如下：

$\begin{array}{c}{M}_1\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]\\ M_2\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& 0& -\sin \mathrm{\α}\\ 0& 1& 0\\ \sin \mathrm{\α}& 0& \cos \mathrm{\α}\end{array}\right]\\ M_3\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]\end{array}---\left(20\right)$

第四步，节点扰动引力赋值

(1)主动段点质量赋值法

通过多层点质量求解算法求得N个质点的质量M_j、深度D_j和球坐标(j＝1,2,…,N)，可得扰动引力三分量的计算式：

其中，ρ_P＝R+H_P为计算点的地心距离，其它量的计算式为：R_j＝R-D_j。令为P点的球坐标，则

$r_{\mathrm{Pj}}={({\mathrm{\ρ}}_P^2+R_j^2-2{\mathrm{\ρ}}_P{R}_j\cos {\mathrm{\ψ}}_{\mathrm{Pj}})}^{1/2}---\left(22\right)$

(2)被动段球谐函数法赋值法

地球外部空间任意点P相对于旋转地球的地心距r、地心纬度经度λ为已知时，由球谐函数级数形式表示的地球引力位V为：

$V=\frac{\mathrm{\μ}}{r}[1+\underset{n=2}{\overset{s}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a_e}{r}\right)}^n\·({\stackrel{\‾}{C}}_{\mathrm{nm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin \mathrm{m\λ})\·{\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\sin \mathrm{\φ}\right)]---\left(26\right)$

令除主要位系数外的位系数都为零，得地球正常引力位

$\tilde{U}=\frac{\mathrm{\μ}}{r}[1+{\stackrel{\‾}{C}}_{20}{\left(\frac{a_e}{r}\right)}^2{\stackrel{\‾}{P}}_{20}\left(\sin \mathrm{\φ}\right)]---\left(27\right)$

真实引力位与正常引力位之差，即为扰动引力位T：

$T=V-\tilde{U}=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{s}{\mathrm{\Σ}}}{\left(\frac{a_e}{r}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}\·({\stackrel{\‾}{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\sin \mathrm{\φ}\right)]---\left(28\right)$

其中：

与T相应的引力即为扰动引力加速度即

$\mathrm{\δ}\stackrel{\→}{g}=\mathrm{grad}T---\left(30\right)$

则扰动引力加速度在天东北坐标系O_E-REN中的三个分量δg_R,δg_E,δg_N为：

第五步，当前计算点单元判断

(1)主动段单元判断方法

①假设第一个单元序号为1；

②确定实际弹道上某点A在发射坐标系中坐标为A(x^*,y^*,z^*)；

③比较y^*与基准点d_i的坐标y_i,0(i＝0,1,2,…)的大小。若y^*≥y_i,0且y^*≤y_i+1,0，则可确定A所在的单元序号为i+1。

(2)被动段单元判断算法

①假设第一个单元序号为1；

②确定实际弹道上某点A在发射坐标系中坐标为A(x^*,y^*,z^*)；

③求出A点对应的真近点角f^*；

④比较f^*与基准点d_j对应的真近点角y_j,0(j＝0,1,2,…)的大小。若f^*≥f_j0且f^*≤f_j+1,0，则可确定A所在的单元序号为j+1。

第六步，单元内部逼近计算

基于网函数逼近理论进行单元内部的逼近计算。以主动段某计算单元为例，记8个顶点对应的扰动引力值分别为g_i,1、g_i,2、g_i,3、g_i,4、g_i+1,1、g_i+1,2、g_i+1,3、g_i+1,4；称12条棱为计算单元上的1-网，定义在其上的函数为1-网函数，记为f_i(x,y,z)(i＝0,1,…,11)。令L(x)、L(y)、L(z)分别为关于x、y、z的一次Lagrange插值算符，其插值基函数为

令为3维1-网函数插值算符，则

将作用于1-网函数，即可求得单元内任意一点A(x,y,z)的值F(x,y,z)，

垂直于飞行弹道的1-网函数由其对应的两节点插值求得，以f₆(x,y,z)的求取为例：

沿飞行弹道的1-网函数由其对应的两节点及两侧节点通过加权三点插值求得，以f₉(x,y,z)的求取为例：

$\begin{array}{c}{f}_9(x,y,z)=\frac{1}{2}L\left(y\right)\{g_{i-\mathrm{1,4}},g_{i,4},g_{i+\mathrm{1,4}}\}+\frac{1}{2}L\left(y\right)\{g_{i,4},g_{i+\mathrm{1,4}},g_{i+\mathrm{2,4}}\}\\ =\frac{1}{2}{\left[\begin{array}{c}\frac{(y-y_i)(y-y_{i+1})}{\mathrm{\δ}{y}_{i-1}(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i-1})}\\ -\frac{(y-y_{i-1})(y-y_{i+1})}{\mathrm{\δ}{y}_{i-1}\mathrm{\δ}{y}_i}\\ \frac{(y-y_i)(y-y_{i-1})}{\mathrm{\δ}{y}_i(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i-1})}\end{array}\right]}^T\left[\begin{array}{c}{g}_{i-\mathrm{1,4}}\\ g_{i,4}\\ g_{i+\mathrm{1,4}}\end{array}\right]+=\frac{1}{2}{\left[\begin{array}{c}\frac{(y-y_{i+1})(y-y_{i+2})}{\mathrm{\δ}{y}_i(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i+1})}\\ -\frac{(y-y_i)(y-y_{i+2})}{\mathrm{\δ}{y}_i\mathrm{\δ}{y}_{i+1}}\\ \frac{(y-y_i)(y-y_{i+1})}{\mathrm{\δ}{y}_{i+1}(\mathrm{\δ}{y}_i+\mathrm{\δ}{y}_{i+1})}\end{array}\right]}^T\end{array}\left[\begin{array}{c}{g}_{i,4}\\ g_{i+\mathrm{1,4}}\\ g_{i+\mathrm{2,4}}\end{array}\right]---\left(36\right)$

图5-图7给出了将本发明及1080阶球谐函数法应用于不同发射点、不同发射方位角及不同射程弹道计算得到的扰动引力赋值结果。表1-表3为以1080阶球谐函数法计算结果为近似真值，对本发明赋值误差的统计结果。

表1不同发射点弹道扰动引力赋值误差统计结果(单位：mgal)

表2不同发射方位角弹道扰动引力赋值误差统计结果(单位：mgal)

表3不同射程弹道扰动引力赋值误差统计结果(单位：mgal)

由表1-表3可见，本发明的平均赋值误差控制在0.01mgal量级，可见本发明具有较高的赋值精度

由图5-图7可见，本发明计算得到的沿弹道的全程扰动引力赋值结果与1080阶球谐函数法计算得到的沿弹道的全程扰动引力赋值结果能较高程度的吻合。

图8给出了将本发明应用于不同射程、不同发射方位角弹道计算得到的由赋值误差导致的落点偏差。图中将该方法应用于弹道计算中得到的弹道落点修正量与将1080阶球谐函数应用于弹道计算中得到的弹道落点修正量进行比较，二者之差即为由逼近误差导致的落点偏差。结果表明，各种情况下由逼近误差导致的落点偏差均不超过8m，说明该方法具有极小的方法误差，能满足弹道计算的要求。

表4扰动引力模型重构速度和存储量分析结果

表4给出了本发明扰动引力模型重构速度和存储量分析结果，在算例所示的网格划分下，对于某射程为12000km的弹道，主动段只需23个网格、600个存储数据，被动段只需19个网格、480个存储数据，全程只需17.531s的时间即可完成扰动引力模型重构。重构速度和数据存储量满足快速发射及弹上计算的要求。

表5不同方法应用于弹道计算的时间对比分析结果

表5为将不同扰动引力赋值方法应用于弹道计算的时间对比分析结果。结果表明，本发明的计算时间相比1080阶球谐函数法及Stokes方法提高了30多倍，因此在计算速度上具有极大的优势。

综合上述仿真结果可获得以下结论：

1)本发明的逼近方法误差小，由此导致的落点偏差小于8m，满足赋值精度要求；

2)本发明以牺牲可容忍的小量精度为代价，极大地提高了计算速度，使扰动引力射前模型快速构建和弹上计算变为可能，具有突破性的工程价值；

3)本发明提出的方法对弹载计算机存储量要求小，能够满足弹上实时计算的要求；

4)本发明能够适应各种情况下的沿弹道扰动引力赋值计算，弹道积分无奇点。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种沿飞行弹道的扰动引力快速逼近方法，包括以下步骤：

首先，根据发射任务要求确定一条不考虑扰动引力的标准弹道；其次，基于标准弹道生成“漏斗形”飞行管道，在飞行管道内完成空域剖分并确定节点位置；再次，采用点质量法、球谐函数法计算得到的扰动引力对节点赋值；最后，在诸元或弹上导航计算时，根据导弹位置判断其所在的单元，并根据当前位置与所在单元各节点相对位置关系以及各节点的扰动引力值，快速计算当前位置对应扰动引力值。

2.根据权利要求1所述的一种沿飞行弹道的扰动引力快速逼近方法，其特征在于，包括以下步骤：

A、生成标准弹道：根据发射任务要求，生成不考虑扰动引力的标准弹道；

B、生成飞行管道：

1)主动段飞行管道生成

在发射坐标系内按如下方式生成主动段飞行管道：

①在标准弹道上依次选定基准点d₀,d₁,d₂,d₃...，相邻两点间沿y轴方向的距离为δy₁,δy₂,δy₃...；

②以d_i为几何中心，在垂直于y轴的平面内生成边长分别为δx_i和δz_i的四边形；

③令δx_i＜δx_i+1，δy_i＜δy_i+1，δz_i＜δz_i+1，依次连接各长方形顶点k_i，构成主动段“漏斗形”飞行管道；

2)被动段飞行管道生成

在局部轨道坐标系内按如下方式生成被动段飞行管道：

①在标准弹道上选定基准点d₀,d₁,d₂,...，各点对应的真近点角分别为f₀,f₁,f₂,...，设δf_j＝f_j-f_j-1(j＝1,2,...)；

②以d_j为原点，建立局部轨道坐标系d_j-r_jβ_jξ_j，r_j轴沿地心矢方向，ξ_j轴与动量矩方向一致，β_j轴与r_j、ξ_j轴构成右手坐标系；

③在轨道坐标系内垂直于β_j轴的平面内生成边长分别为δr_j和δξ_j的四边形；

④令δr_j＜δr_j+1，δξ_j＜δξ_j+1，δf_j＜δf_j+1，依次连接各长方形顶点k_j，构成被动段“漏斗形”飞行管道；

C、节点坐标位置确定：

1)主动段节点坐标位置确定

根据沿飞行弹道的空域剖分方法，确定基准点d_i对应的四个节点在发射系中的坐标分别为k₁(x_i0+δx_i/2,y_i0,z_i0+δz_i/2)，k₂(x_i0+δx_i/2,y_i0,z_i0-δz_i/2)，k₃(x_i0-δx_i/2,y_i0,z_i0-δz_i/2)，k₄(x_i0-δx_i/2,y_i0,z_i0+δz_i/2)；

2)被动段节点坐标位置确定

根据沿飞行弹道的空域剖分方法，确定基准点d_j对应的四个节点在局部轨道坐标系中的坐标分别为k₁(δr_j/2,0,δξ_j/2)，k₂(δr_j/2,0,-δξ_j/2)，k₃(-δr_j/2,0,-δξ_j/2)，k₄(-δr_j/2,0,δξ_j/2)；

设被动段实际弹道上某点A在地心惯性坐标系中坐标为(x′,y′,z′)，在局部轨道坐标系中的坐标为(x,y,z)，根据如下坐标转换关系，获得节点在地心惯性系中的坐标：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=M_3(\mathrm{\ω}+f)\·M_1\left(i\right)\·M_3\left(\mathrm{\Ω}\right)\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]-\left[\begin{array}{c}r\\ 0\\ 0\end{array}\right]---I$

其中，r为地心距r，f为真近点角，ω为近地点角距，Ω为升交点角距，i为轨道倾角；M_i(i＝1,2,3)为坐标转换矩阵，具体形式如下：

$M_1\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]$

$M_2\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& 0& -\sin \mathrm{\α}\\ 0& 1& 0\\ \sin \mathrm{\α}& 0& \cos \mathrm{\α}\end{array}\right]---\mathrm{II}$

$M_3\left(\mathrm{\α}\right)=\left[\begin{array}{ccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]$

D、节点扰动引力赋值：

1)、主动段点质量赋值法：

通过多层点质量求解算法求得N个质点的质量M_j、深度D_j和球坐标(j＝1,2,…,N)，可得扰动引力三分量的计算式：

其中，ρ_P＝R+H_P为计算点的地心距离，其它量的计算式为：R_j＝R-D_j；令为P点的球坐标，则

$r_{\mathrm{Pj}}={({\mathrm{\ρ}}_P^2+R_j^2-2{\mathrm{\ρ}}_P{R}_j\cos {\mathrm{\ψ}}_{\mathrm{Pj}})}^{1/2}---\mathrm{IV}$

2)、被动段球谐函数法赋值法

地球外部空间任意点P相对于旋转地球的地心距r、地心纬度经度λ为已知时，由球谐函数级数形式表示的地球引力位V为：

$V=\frac{\mathrm{\μ}}{r}[1+\underset{n=2}{\overset{s}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a_e}{r}\right)}^n\·({\stackrel{\‾}{C}}_{\mathrm{nm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin \mathrm{m\λ})\·{\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\sin \mathrm{\φ}\right)]---\mathrm{VIII}$

令除主要位系数外的位系数都为零，得地球正常引力位

$\tilde{U}=\frac{\mathrm{\μ}}{r}[1+{\stackrel{\‾}{C}}_{20}{\left(\frac{a_e}{r}\right)}^2{\stackrel{\‾}{P}}_{20}\left(\sin \mathrm{\φ}\right)]---\mathrm{IX}$

真实引力位与正常引力位之差，即为扰动引力位T：

$T=V-\tilde{U}=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{s}{\mathrm{\Σ}}}{\left(\frac{a_e}{r}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\sin \mathrm{\φ}\right)---X$

其中：

与T相应的引力即为扰动引力加速度即

$\mathrm{\δ}\stackrel{\→}{g}=\mathrm{grad}T---\mathrm{XII}$

则扰动引力加速度在天东北坐标系O_E-REN中的三个分量δg_R,δg_E,δg_N为：

E、当前计算点单元判断：

1)主动段单元判断方法的过程如下：

①假设第一个单元序号为1；

②确定实际弹道上某点A在发射坐标系中坐标为A(x^*,y^*,z^*)；

③比较y^*与基准点d_i的坐标y_i,0(i＝0,1,2,…)的大小；若y^*≥y_i,0且y^*≤y_i+1,0，则可确定A所在的单元序号为i+1；

2)被动段单元判断算法的过程如下：

①假设第一个单元序号为1；

②确定实际弹道上某点A在发射坐标系中坐标为A(x^*,y^*,z^*)；

③求出A点对应的真近点角f^*；

④比较f^*与基准点d_j对应的真近点角y_j,0(j＝0,1,2,…)的大小；若f^*≥f_j,0且f^*≤f_j+1,0，则可确定A所在的单元序号为j+1；

F、单元内部逼近计算；

基于网函数逼近理论进行单元内部的逼近计算；以主动段某计算单元为例，记8个顶点对应的扰动引力值分别为g_i,1、g_i,2、g_i,3、g_i,4、g_i+1,1、g_i+1,2、g_i+1,3、g_i+1,4；称12条棱为计算单元上的1-网，定义在其上的函数为1-网函数，记为f_i(x,y,z)(i＝0,1,…,11)；令L(x)、L(y)、L(z)分别为关于x、y、z的一次Lagrange插值算符，其插值基函数为：

令l⁽³⁾为3维1-网函数插值算符，则

l⁽³⁾＝L(x)L(y)+L(y)L(z)+L(z)L(x)-2L(x)L(y)L(z)   ⅩⅤ

将l⁽³⁾作用于1-网函数，即可求得单元内任意一点A(x,y,z)的值F(x,y,z)，

垂直于飞行弹道的1-网函数由其对应的两节点插值获得，沿飞行弹道的1-网函数由其对应的两节点及两侧节点通过加权三点插值获得。

3.根据权利要求2所述的一种沿飞行弹道的扰动引力快速逼近方法，其特征在于：

步骤F中，垂直于飞行弹道的1-网函数由其对应的两节点插值获得过程以f₆(x,y,z)的求取为例：

4.根据权利要求2所述的一种沿飞行弹道的扰动引力快速逼近方法，其特征在于：

步骤F中，沿飞行弹道的1-网函数由其对应的两节点及两侧节点通过加权三点插值获得过程以f₉(x,y,z)的求取为例：

$\begin{array}{c}{f}_9(x,y,z)=\frac{1}{2}L\left(y\right)\{g_{i-\mathrm{1,4}},g_{i,4},g_{i+\mathrm{1,4}}\}+\frac{1}{2}L\left(y\right)\{g_{i,4},g_{i+\mathrm{1,4}},g_{i+\mathrm{2,4}}\}\\ =\frac{1}{2}{\left[\begin{array}{c}\frac{(y-y_i)(y-y_{i+1})}{{\mathrm{\δy}}_{i-1}({\mathrm{\δy}}_i+{\mathrm{\δy}}_{i-1})}\\ -\frac{(y-y_{i-1})(y-y_{i+1})}{{\mathrm{\δy}}_{i-1}{\mathrm{\δy}}_i}\\ \frac{(y-y_i)(y-y_{i-1})}{{\mathrm{\δy}}_i({\mathrm{\δy}}_i+{\mathrm{\δy}}_{i-1})}\end{array}\right]}^T\left[\begin{array}{c}{g}_{i-\mathrm{1,4}}\\ g_{i,4}\\ g_{i+\mathrm{1,4}}\end{array}\right]+=\frac{1}{2}{\left[\begin{array}{c}\frac{(y-y_{i+1})(y-y_{i+2})}{{\mathrm{\δy}}_i({\mathrm{\δy}}_i+{\mathrm{\δy}}_{i+1})}\\ -\frac{(y-y_i)(y-y_{i+2})}{{\mathrm{\δy}}_i{\mathrm{\δy}}_{i+1}}\\ \frac{(y-y_i)(y-y_{i+1})}{{\mathrm{\δy}}_{i+1}({\mathrm{\δy}}_i+{\mathrm{\δy}}_{i+1})}\end{array}\right]}^T\left[\begin{array}{c}{g}_{i,4}\\ g_{i+\mathrm{1,4}}\\ g_{i+\mathrm{2,4}}\end{array}\right]\end{array}----\mathrm{XVIII}$