CN105740506A - 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 - Google Patents

Abstract

本发明公开了一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法，包括以下步骤：建立弹道包络侧向宽度与纵程的数学模型；建立换极坐标系；建立换极坐标系中飞行器动力学模型；计算换极坐标系弹道三维包络；换极坐标系全域空域剖分；建立一般坐标系空间包络扰动引力重构模型；当前弹道位置所在网格判断；建立计算网格内局域坐标系；网格内部扰动引力网函数逼近计算。本发明涉及的弹道空间包络计算方法可覆盖任意情况下的临近空间大范围机动弹道，根据本发明方法进行沿弹道的扰动引力计算，具有计算规模小、计算速度快、计算精度高的特征，能满足弹道规划、制导计算等对计算量和计算速度的要求。

Description

沿临近空间大范围机动弹道空间包络的扰动引力逼近方法

技术领域

本发明涉及一种飞行器动力学建模领域的沿临近空间大范围机动弹道空间包络的扰动引力快速计算方法，尤其涉及一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法。

背景技术

临近空间大范围机动弹道特指高超声速滑翔飞行器包含初始下降段弹道和滑翔段弹道在内的弹道，其起点为再入点，终点为滑翔段弹道终点。该类弹道长时间处于临近空间，具有侧向大范围机动的特性，受地球外部空间扰动引力影响显著。扰动引力通过影响导航解算而影响飞行器当前状态量的确定，进而对制导计算和弹道精度产生影响。有计算表明，由扰动引力引起的临近空间大范围机动弹道终端位置偏差可达几十公里，飞行器可能因无法进入末制导有效半径而导致飞行任务失败，因此有必要在导航解算和制导计算中考虑扰动引力的影响。

对于采用纯惯性导航方案的飞行器，弹载设备无法测量扰动引力，引力项精度只能通过射前构建扰动引力模型并进行飞行过程中实时计算来保证。射前模型的重构速度受飞行器发射快速性指标的约束，实时计算方法的计算量、计算速度和存储量等受弹载计算机数据存储能力和计算能力的约束。现有地球物理和大地测量领域的传统计算方法，通常具有数据存储量大、计算规模大、计算速度慢的特征，难以满足射前引力模型快速构建及飞行过程中快速计算的要求。与沿弹道式弹道的扰动引力快速计算问题相比，沿临近空间大范围机动弹道的扰动引力快速计算具有更加严格的要求。首先，弹道侧向机动范围可达几千公里，弹道形态因飞行任务不同而迥异，因此建立的方法须能适应射前飞行任务临时变更的情况。其次，由于弹道具有大范围侧向机动的特性，需合理确定逼近空域范围以防止实际弹道超出逼近区间。最后，由于弹道处于扰动引力变化复杂的临近空间，对逼近方法的精度具有更高要求。

目前尚没有一种成熟的方法，能够达到沿任意临近空间大范围机动弹道空间包络的扰动引力快速计算的要求。

发明内容

本发明的目的就在于为了解决上述问题而提供一种逼近精度高、计算速度快、数据存储量小的沿临近空间大范围机动弹道空间包络的扰动引力逼近方法。

本发明通过以下技术方案来实现上述目的：

一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法，包括以下步骤：

(1)建立弹道包络侧向宽度与纵程的数学模型；

(2)建立换极坐标系；

(3)建立换极坐标系中飞行器动力学模型；

(4)计算换极坐标系弹道三维包络；

(5)换极坐标系全域空域剖分；

(6)建立一般坐标系空间包络扰动引力重构模型；

(7)当前弹道位置所在网格判断；

(8)建立计算网格内局域坐标系；

(9)网格内部扰动引力网函数逼近计算。

作为优选，所述步骤(1)的方法为：

令弹道再入点为I，其经度为λ_I、地心纬度为φ_I，弹道落点为T，其经度为λ_T、地心纬度为φ_T；以由再入点和落点确定的再入大圆弧平面为对称面，称沿射向方向对称面以左弹道为左侧机动弹道，沿射向方向对称面以右弹道为右侧机动弹道；记最大左侧机动弹道距对称面的最大距离为左边界B_l，记最大右侧机动弹道距对称面的最大距离为右边界B_r；

根据给定的飞行器再入点飞行状态参数、弹道终端飞行状态参数、飞行过程约束条件及终端约束条件，针对再入点为(λ_I,φ_I)＝(0°,0°)、落点为(λ_ki,0°)(i＝1,2,3,4,…)的低空大范围机动弹道进行弹道计算，记各落点对应的弹道纵程为L_i，由下式(1)计算L_i：

L_i＝R_earccos(sinφ_Isinφ_T+cosφ_Icosφ_Tcos(λ_T-λ_I))(1)

其中，R_e为地球半径；

通过改变再入点初始速度方位角获得最大左侧机动弹道和最大右侧向机动弹道，记各落点对应的弹道侧向包络宽度为W_io，由下式(2)计算W_io：

W_io＝B_l-B_r(2)

取C_w倍W_io为实际侧向包络宽度W_i，对L_i(i＝1,2,3,4,…)和W_i(i＝1,2,3,4,…)进行拟合，可建立如下式(3)所示的弹道包络侧向宽度与纵程的数学模型，确定模型系数为a_i：

$W=\underset{i=0}{\overset{term-1}{\Σ}}{a}_i{L}^i---\left(3\right);$

所述步骤(2)的方法为：

引入一个换极坐标系：用表示换极坐标系中各物理量，用X表示一般坐标系中各物理量，按如下方式建立换极坐标系：

首先，定义一个再入大圆弧平面作为换极赤道平面：对目标点确定的情况，将再入点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面；对于目标点未确定的情况，根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面；

然后，基于换极赤道平面定义换极坐标系O_E为地心，轴沿再入点地心矢径方向，轴在换极赤道面内垂直于轴指向目标点方向，轴与轴、轴构成右手系；

所述步骤(3)的方法为：

在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程，其飞行状态量为换极后的经度地心纬度航迹偏航角速度速度倾角和地心距

$\left\{\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \widehat{\mathrm{\θ}}\\ \frac{d\widehat{\mathrm{\λ}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\σ}}}{\hat{r}\cos \widehat{\mathrm{\φ}}}\\ \frac{d\widehat{\mathrm{\φ}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\σ}}}{\hat{r}}\\ \frac{d\widehat{\mathrm{\σ}}}{dt}=\frac{L\sin \mathrm{\υ}}{\hat{V}\cos \widehat{\mathrm{\θ}}}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\θ}}\tan \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\σ}}-{\hat{g}}_{{\mathrm{\ω}}_e}\frac{\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\σ}}}{\hat{V}\cos \widehat{\mathrm{\θ}}}+C_{\mathrm{\σ}}+{\tilde{C}}_{\mathrm{\σ}}\\ \frac{d\widehat{\mathrm{\θ}}}{dt}=\frac{L}{\hat{V}}\cos \mathrm{\υ}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\θ}}+\frac{{\hat{g}}_r\cos \widehat{\mathrm{\θ}}}{\hat{V}}+\frac{{\hat{g}}_{{\mathrm{\ω}}_e}}{\hat{V}}(\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\φ}}-\cos \widehat{\mathrm{\σ}}\sin \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\φ}})+C_{\mathrm{\θ}}+{\tilde{C}}_{\mathrm{\θ}}\\ \frac{d\hat{V}}{dt}=-D+{\hat{g}}_r\sin \widehat{\mathrm{\θ}}+{\hat{g}}_{{\mathrm{\ω}}_e}(\cos \widehat{\mathrm{\σ}}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\φ}}+\sin \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\φ}})+{\tilde{C}}_V\end{array}\right.---\left(4\right)$

其中，C_σ、C_θ为哥氏加速度项，和为牵连加速度项，

$\left\{\begin{array}{c}{C}_{\mathrm{\σ}}=(2{\mathrm{\ω}}_{ex}-2\tan \widehat{\mathrm{\θ}}({\mathrm{\ω}}_{ey}\sin \widehat{\mathrm{\σ}}+{\mathrm{\ω}}_{ez}\cos \widehat{\mathrm{\σ}}\left)\right)\\ {\tilde{C}}_{\mathrm{\σ}}=-\frac{\hat{r}}{\hat{V}\cos \widehat{\mathrm{\θ}}}({\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ey}\cos \widehat{\mathrm{\σ}}-{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ez}\sin \widehat{\mathrm{\σ}})\\ C_{\mathrm{\θ}}=2({\mathrm{\ω}}_{ez}\sin \widehat{\mathrm{\σ}}-{\mathrm{\ω}}_{ey}\cos \widehat{\mathrm{\σ}})\\ {\tilde{C}}_{\mathrm{\θ}}=\frac{\hat{r}}{\hat{V}}\[{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ey}\sin \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\σ}}+{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ez}\sin \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\σ}}+({\mathrm{\ω}}_{ey}^2+{\mathrm{\ω}}_{ez}^2)\cos \widehat{\mathrm{\θ}}\]\\ {\tilde{C}}_V=\hat{r}\[-{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ey}\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\σ}}-{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ez}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\σ}}+({\mathrm{\ω}}_{ey}^2+{\mathrm{\ω}}_{ez}^2)\sin \widehat{\mathrm{\θ}}\]\end{array}\right.---\left(5\right)$

其中，

$\{\begin{array}{c}{\mathrm{\ω}}_{ex}={\mathrm{\ω}}_e(\cos \widehat{\mathrm{\λ}}\cos \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p{\mathrm{cosA}}_p+\sin \widehat{\mathrm{\λ}}\cos \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p{\mathrm{sinA}}_p+\sin \widehat{\mathrm{\φ}}{\mathrm{sin\φ}}_p)\\ {\mathrm{\ω}}_{ey}={\mathrm{\ω}}_e(-\sin \widehat{\mathrm{\λ}}{\mathrm{cos\φ}}_p{\mathrm{cosA}}_p+\cos \widehat{\mathrm{\λ}}{\mathrm{cos\φ}}_p{\mathrm{sinA}}_p)\\ {\mathrm{\ω}}_{ez}={\mathrm{\ω}}_e(-\cos \widehat{\mathrm{\λ}}\sin \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p{\mathrm{cosA}}_p-\sin \widehat{\mathrm{\λ}}\sin \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p{\mathrm{sinA}}_p+\cos \widehat{\mathrm{\φ}}{\mathrm{sin\φ}}_p)\end{array}---\left(6\right)$

其中，ω_e为地球旋转加速度矢量，λ_p和φ_p为换极后极点P的经度和地心纬度，A_P为P的方位角；

根据换极坐标系定义，一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定义一致，

$r=\hat{r},\mathrm{\θ}=\widehat{\mathrm{\θ}},V=\hat{V}---\left(7\right)$

定义

$\begin{array}{c}\left[\begin{array}{ccc}{\mathrm{cos\φ}}_f{\mathrm{cos\λ}}_f& {\mathrm{cos\φ}}_f{\mathrm{sin\λ}}_f& {\mathrm{sin\φ}}_f\\ -{\mathrm{sin\ψ}}_f{\mathrm{sin\λ}}_f-{\mathrm{cos\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{cos\λ}}_f& {\mathrm{sin\ψ}}_f{\mathrm{cos\λ}}_f-{\mathrm{cos\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{sin\λ}}_f& {\mathrm{cos\ψ}}_f{\mathrm{cos\φ}}_f\\ {\mathrm{cos\ψ}}_f{\mathrm{sin\λ}}_f-{\mathrm{sin\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{cos\λ}}_f& -{\mathrm{cos\ψ}}_f{\mathrm{cos\λ}}_f-{\mathrm{sin\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{sin\λ}}_f& {\mathrm{sin\ψ}}_f{\mathrm{cos\φ}}_f\end{array}\right]\\ \stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right]\end{array}---\left(8\right)$

其中，ψ_f为点F的方位角，

$\{\begin{array}{c}{G}_{11}\cos \mathrm{\φ}\cos \mathrm{\λ}+G_{12}\cos \mathrm{\φ}\sin \mathrm{\λ}+G_{13}\sin \mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{k}_1\\ G_{21}\cos \mathrm{\φ}\cos \mathrm{\λ}+G_{22}\cos \mathrm{\φ}\sin \mathrm{\λ}+G_{23}\sin \mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{k}_2\\ G_{31}\cos \mathrm{\φ}\cos \mathrm{\λ}+G_{32}\cos \mathrm{\φ}\sin \mathrm{\λ}+G_{33}\sin \mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{k}_3\end{array}---\left(9\right)$

$\{\begin{array}{c}{G}_{11}\cos \widehat{\mathrm{\φ}}\cos \widehat{\mathrm{\λ}}+G_{21}\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\λ}}+G_{31}\sin \widehat{\mathrm{\φ}}\stackrel{\mathrm{\Δ}}{=}{\tilde{k}}_1\\ G_{12}\cos \widehat{\mathrm{\φ}}\cos \widehat{\mathrm{\λ}}+G_{22}\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\λ}}+G_{32}\sin \widehat{\mathrm{\φ}}\stackrel{\mathrm{\Δ}}{=}{\tilde{k}}_2\\ G_{13}\cos \widehat{\mathrm{\φ}}\cos \widehat{\mathrm{\λ}}+G_{23}\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\λ}}+G_{33}\sin \widehat{\mathrm{\φ}}\stackrel{\mathrm{\Δ}}{=}{\tilde{k}}_3\end{array}---\left(10\right)$

由一般坐标系中λ和φ确定换极坐标系中和的表达式为：

$\begin{array}{cc}\left\{\begin{array}{c}\cos \widehat{\mathrm{\λ}}=k_1/\sqrt{{k_1}^2+{k_2}^2}\\ \sin \widehat{\mathrm{\λ}}=k_2/\sqrt{{k_1}^2+{k_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \widehat{\mathrm{\φ}}=k_3\\ \cos \widehat{\mathrm{\φ}}=\sqrt{{k_1}^2+{k_2}^2}\end{array}\right.---\left(11\right)\end{array}$

由换极坐标系中λ和φ确定一般坐标系中λ和φ的表达式为：

$\begin{array}{cc}\left\{\begin{array}{c}\cos \mathrm{\λ}={\tilde{k}}_1/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\\ \sin \mathrm{\λ}={\tilde{k}}_2/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \mathrm{\φ}={\tilde{k}}_3\\ \cos \mathrm{\φ}=\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.\end{array}---\left(12\right)$

由一般坐标系中σ确定换极坐标系中的表达式为：

$\widehat{\mathrm{\σ}}=\mathrm{\σ}+\mathrm{\η}---\left(13\right)$

其中，

$\{\begin{array}{c}\sin \mathrm{\η}=\frac{\sin (\mathrm{\λ}-{\mathrm{\λ}}_P){\mathrm{cos\φ}}_P}{\cos \widehat{\mathrm{\φ}}}\\ \cos \mathrm{\η}=-\cos (A_P-\widehat{\mathrm{\λ}})\cos (\mathrm{\λ}-{\mathrm{\λ}}_P)+\sin (A_P-\widehat{\mathrm{\λ}})\sin (\mathrm{\λ}-{\mathrm{\λ}}_P){\mathrm{sin\φ}}_P\end{array}---\left(14\right);$

所述步骤(4)的方法为：

根据步骤(3)所述坐标转换关系，根据弹道再入点经度λ_I、再入点地心纬度φ_I、落点经度λ_T、落点地心纬度φ_T计算得到换极弹道参数，由计算得 记换极弹道落点经度为

由式(15)计算换极弹道纵程

$\hat{L}=R_{e}arccos(sin{\widehat{\mathrm{\φ}}}_{I}sin{\widehat{\mathrm{\φ}}}_T+cos{\widehat{\mathrm{\φ}}}_{I}cos{\widehat{\mathrm{\φ}}}_{T}cos({\widehat{\mathrm{\λ}}}_T-{\widehat{\mathrm{\λ}}}_I\left)\right)---\left(15\right)$

其中，R_e为地球半径；

将代入式(16)，可计算得到换极弹道侧向包络宽度

$\hat{W}=\underset{i=0}{\overset{term-1}{\Σ}}{a}_i{\hat{L}}^i---\left(16\right)$

其中，模型系数a_i由步骤(1)确定；

将换极弹道侧向包络描述为长度为宽度为的矩形，其边界由式(17)确定：

$\begin{array}{cc}\{\begin{array}{c}{\hat{E}}_{\mathrm{\λ}1}={\mathrm{\λ}}_I\\ {\hat{E}}_{\mathrm{\λ}2}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\end{array},& \left\{\begin{array}{c}{\hat{E}}_{\mathrm{\φ}1}=\frac{90\hat{W}}{{\mathrm{\πR}}_e}\\ {\hat{E}}_{\mathrm{\φ}2}=\frac{90\hat{W}}{{\mathrm{\πR}}_e}\end{array}\right.\end{array}---\left(17\right)$

其中，为换极系侧向包络东向下边界，为换极系侧向包络东向上边界；为换极系侧向包络北向下边界，为换极系侧向包络北向上边界；

侧向包络对应的经度范围Δλ及地心纬度范围Δφ由式(18)确定：

$\left\{\begin{array}{c}\mathrm{\Δ}\widehat{\mathrm{\λ}}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\\ \mathrm{\Δ}\widehat{\mathrm{\φ}}=\frac{180\hat{W}}{{\mathrm{\πR}}_e}\end{array}\right.---\left(18\right)$

以弹道高度范围的C_r倍作为弹道包络的垂向范围

所述步骤(5)的方法为：

记换极系空间包络确定的空域为Ω，按式(19)计算Ω边界：

$\begin{array}{ccc}\{\begin{array}{c}{\hat{F}}_{r1}={\hat{r}}_I-\mathrm{\Δ}\hat{r}-2\·dr\\ {\hat{F}}_{r2}={\hat{r}}_I+2\·dr\end{array},& \{\begin{array}{c}{\hat{F}}_{\mathrm{\λ}1}={\widehat{\mathrm{\λ}}}_I-2\·d\mathrm{\λ}\\ {\hat{F}}_{\mathrm{\λ}2}={\widehat{\mathrm{\λ}}}_I+\mathrm{\Δ}\widehat{\mathrm{\λ}}+2\·d\mathrm{\λ}\end{array},& \left\{\begin{array}{c}{\hat{F}}_{\mathrm{\φ}1}={\widehat{\mathrm{\φ}}}_I-\mathrm{\Δ}\widehat{\mathrm{\φ}}/2\\ {\hat{F}}_{\mathrm{\φ}2}={\widehat{\mathrm{\φ}}}_I+\mathrm{\Δ}\widehat{\mathrm{\φ}}/2\end{array}\right.---\left(19\right)\end{array}$

其中，为换极系Ω天向下边界，为换极系Ω天向上边界，为换极系Ω东向下边界，为换极系Ω东向上边界，为换极系Ω北向下边界，为换极系Ω北向上边界，dr、dλ和dφ分别为空域剖分间隔；

按照天向东向北向的间隔将由式(19)确定的空间包络空域Ω均匀剖分为q个互不重叠的子域Ω_e(e＝1,2,…,q)，记网格节点坐标为

$\mathrm{\Ω}=\underset{e=1}{\overset{q}{\∪}}{\mathrm{\Ω}}_e---\left(20\right);$

所述步骤(6)的方法为：

根据换极系网格节点坐标由步骤(3)所述坐标转换关系计算一般坐标系网格节点坐标

由球谐函数方法计算一般坐标系网格节点扰动引力位T_G，

$T_G=\frac{\mathrm{\μ}}{r_{gi}}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_{gj}+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_{gj}){\stackrel{\‾}{P}}_{mn}\left({\mathrm{sin\φ}}_{gk}\right)---\left(21\right)$

其中，

${\stackrel{\‾}{C}}_{nm}^{*}=\left\{\begin{array}{ccc}0& n=2& m=0\\ {\stackrel{\‾}{C}}_{nm}& else& \end{array}\right.---\left(22\right)$

与T_G相应的引力即为扰动引力加速度δg，即

δg＝gradT_G(23)

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

$\left\{\begin{array}{c}{\mathrm{\δg}}_R=\frac{\∂T}{\∂r_{gi}}=-\frac{\mathrm{\μ}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\Σ}}(n+1){\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_{gj}+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_{gj}){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_{gk}\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_{gi}{\mathrm{cos\φ}}_{gk}}\frac{\∂T}{\∂{\mathrm{\λ}}_{gj}}=-\frac{1}{{\mathrm{cos\φ}}_{gk}}\frac{\mathrm{\μ}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\Σ}}m({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{sinm\λ}}_{gj}+{\stackrel{\‾}{S}}_{nm}{\mathrm{cosm\λ}}_{gj}){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_{gk}\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_{gi}}\frac{\∂T}{\∂{\mathrm{\φ}}_{gk}}=\frac{\mathrm{\μ}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_{gj}+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_{gj})\frac{d{\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_{gk}\right)}{{\mathrm{d\β}}_{gk}}\end{array}\right.---\left(24\right)$

存储一般坐标系节点位置三分量和节点扰动引力三分量，完成空间包络扰动引力重构模型构建；

所述步骤(7)的方法为：

令当前弹道位置坐标为A(r,λ,φ)，当满足如式(25)所示的条件时，

$\{\begin{array}{c}{r}_{gi}\&le;r<r_{gi+1}\\ {\mathrm{\&lambda;}}_{gj}\&le;\mathrm{\&lambda;}<{\mathrm{\&lambda;}}_{gj+1}\\ {\mathrm{\&phi;}}_{gk}\&le;\mathrm{\&phi;}<{\mathrm{\&phi;}}_{gk+1}\end{array}---\left(25\right)$

判断A(r,λ,φ)位于由节点N(r_gi,λ_gj,φ_gk)、N(r_gi,λ_gj+1,φ_gk)、N(r_gi,λ_gj,φ_gk+1)、N(r_gi,λ_gj+1,φ_gk+1)、N(r_gi+1,λ_gj,φ_gk)、N(r_gi+1,λ_gj+1,φ_gk)、N(r_gi+1,λ_gj,φ_gk+1)、N(r_gi+1,λ_gj+1,φ_gk+1)所确定的网格中；

所述步骤(8)的方法为：

局部坐标系由半径r_l＝r_gi+dr/2的球面、经度λ_l＝λ_gj+dλ/2的子午面、纬度φ_l＝φ_gk+dφ/2的纬圈的交线组成，原点l(r_l,λ_l,φ_l)为三交线的交点，原点局部坐标为l(0,0,0)，ξ、η、ζ分别沿原点l的天向、北向和东向，则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为：

$\left\{\begin{array}{c}\mathrm{\&xi;}=r-r_l\\ \mathrm{\&eta;}=rcos\mathrm{\&phi;}(\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_l)\\ \mathrm{\&zeta;}=r(\mathrm{\&phi;}-{\mathrm{\&phi;}}_l)\end{array}\right.---\left(26\right)$

单元顶点A_i(r_i,λ_i,φ_i)的局部坐标A_i(ξ_i,η_i,ζ_i)为：

$\left\{\begin{array}{c}{\mathrm{\&xi;}}_i=r_i-r_i\\ {\mathrm{\&eta;}}_i=r_i{\mathrm{cos\&phi;}}_i({\mathrm{\&lambda;}}_i-{\mathrm{\&lambda;}}_l)\\ {\mathrm{\&zeta;}}_i=r_i({\mathrm{\&phi;}}_i-{\mathrm{\&phi;}}_l)\end{array}\right.---\left(27\right);$

所述步骤(9)的方法为：

记某计算单元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(ξ,η,ζ)(i＝0,1,…,11)；令L(ξ)、L(η)、L(ζ)分别为关于ξ、η、ζ的一次Lagrange插值算符，其插值基函数为：

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

l⁽³⁾＝L(ξ)L(η)+L(η)L(ζ)+L(ζ)L(ξ)-2L(ξ)L(η)L(ζ)

将l⁽³⁾作用于1-网函数，即可求得单元内任意一点A(ξ,η,ζ)的值δg(ξ,η,ζ)，

至此，经过上述九步可最终建立一种沿临近空间大范围机动弹道空间包络的扰动引力网函数逼近方法，方法可适应沿空间包络内任意弹道的扰动引力快速计算要求。

本发明的有益效果在于：

1、本发明首次针对临近空间大范围机动弹道的扰动引力快速计算问题，提出一种基于弹道空间包络和网函数逼近理论的扰动引力计算方法，方法的计算规模、计算速度明显优于传统扰动引力计算方法，可适应沿弹道空间包络内任意弹道的扰动引力快速计算，计算精度满足临近空间大范围机动弹道的精度要求；

2、本方法基于弹道空间包络计算构建扰动引力重构模型，可有效避免实际弹道因侧向大范围机动而超出逼近区间的问题，能覆盖考虑禁飞区及最大侧向机动等极端情况，对各种情况下的弹道均具有良好的适应性；

3、本方法采用的弹道空间包络计算方法能够根据再入点和目标点快速确定弹道可达域，计算得到的包络可覆盖具各种情况下的机动弹道；

4、本方法引入极点变换的思想以简化弹道空间包络边界的确定以及空域剖分建模过程，引入局部坐标系以简化逼近方法的建立；

5、本发明首次将网函数逼近理论应用到针对临近空间大范围机动弹道的扰动引力快速计算中，网函数的插值函数类具有Coons型结构以及明显的统计特征，便于计算机实现，且具有较高的逼近精度。

附图说明

图1为不同射向弹道侧向包络示意图；

图2为东向不同纵程弹道侧向包络示意图；

图3为滑翔弹道三维包络示意图；

图4为换极坐标系及飞行状态量定义示意图；

图5为弹道再入点I和换极后极点P的位置关系示意图；

图6为换极坐标系航迹偏航角的位置关系示意图；

图7-1为换极坐标系空间包络空域剖分节点示意图之俯视图；

图7-2为换极坐标系空间包络空域剖分节点示意图之侧视图；

图8-1为一般坐标系空间包络空域剖分节点示意图之俯视图；

图8-2为一般坐标系空间包络空域剖分节点示意图之侧视图；

图9为空间包络节点扰动引力天向分量示意图；

图10为空间包络节点扰动引力北向分量示意图；

图11为空间包络节点扰动引力东向分量示意图；

图12为网函数计算单元示意图；

图13为沿不同区域弹道扰动引力逼近结果对比示意图；

图14为沿不同区域弹道扰动引力逼近误差对比示意图；

图15为沿不同射向弹道扰动引力逼近结果对比示意图；

图16为沿不同射向弹道扰动引力逼近误差对比示意图；

图17为沿不同射程弹道扰动引力逼近结果对比示意图；

图18为沿不同射程弹道扰动引力逼近误差对比示意图；

图19为扰动引力引起的弹道终端状态偏差(初始方位角0°～360°)示意图；

图20为逼近误差引起的弹道终端状态偏差(初始方位角0°～360°)示意图。

具体实施方式

下面结合实施例和附图对本发明作进一步说明：

实施例：

基于CAV-H飞行器模型进行仿真，基本仿真条件设置为：①质量m＝908kg，参考面积Sm＝0.48375m2；②再入点状态参数：速度V_I＝6500m/s，速度倾角θ_I＝0°，高度H_I＝80km；③滑翔段终端状态参数：速度V_k＝2500m/s，速度倾角θ_k＝0°，高度H_k＝30km，终端航迹偏航角偏差不大于±5°；④飞行约束条件：最大热流密度最大动压q_max＝100kPa，最大过载n_max＝3g；⑤终端结束时距目标点的待飞航程：S_togo＝100km。

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

本发明所述沿临近空间大范围机动弹道空间包络的扰动引力逼近方法，包括以下步骤：

(1)建立弹道包络侧向宽度与纵程的数学模型，具体方法为：

令弹道再入点为I，其经度为λ_I、地心纬度为φ_I，弹道落点为T，其经度为λ_T、地心纬度为φ_T；以由再入点和落点确定的再入大圆弧平面为对称面，称沿射向方向对称面以左弹道为左侧机动弹道，沿射向方向对称面以右弹道为右侧机动弹道；记最大左侧机动弹道距对称面的最大距离为左边界B_l，记最大右侧机动弹道距对称面的最大距离为右边界B_r；

针对不同射向的临近空间大范围机动弹道进行仿真，其再入点经度为λ_I＝0°、再入点地心纬度为φ_I＝0°，落点经度λ_T和落点地心纬度φ_T见表1：

表1不同射向弹道落点坐标

射向	正北	正东	正南	正西
					经度(°)	0	60	0	-60
地心纬度(°)	60	0	-60	0

计算得到不同射向弹道侧向包络示意图见图1，不同射向弹道侧向包络边界见表2。由图1和表2可见，由于地球自转和非球形因素的影响，东向弹道包络最大、南北向弹道包络次之，西向弹道包络最小；南北向弹道包络向东有小范围偏移。

表2不同射向弹道侧向包络边界

射向	正东	正南	正西	正北
					左边界Bl(km)	1597	1072	1072	1072
右边界Br(km)	1597	1513	1072	1528

针对射向为正东、纵程不同的临近空间大范围机动弹道进行仿真。弹道再入点经度为λ_I＝0°，再入点地心纬度为φ_I＝0°；落点地心纬度为0°，落点经度分别为50°、60°、70°和80°。

根据给定的飞行器再入点飞行状态参数、弹道终端飞行状态参数、飞行过程约束条件及终端约束条件，针对再入点为(λ_I,φ_I)＝(0°,0°)、落点为(λ_ki,0°)(i＝1,2,3,4,…)的低空大范围机动弹道进行弹道计算，记各落点对应的弹道纵程为L_i，由下式(1)计算L_i：

L_i＝R_earccos(sinφ_Isinφ_T+cosφ_Icosφ_Tcos(λ_T-λ_I))(1)

其中，R_e为地球半径；

计算得到不同落点对应的纵程L_i见表3。

表3不同落点对应的纵程

落点经度(°)	50	60	70	80
					纵程Li(km)	5560	6671	7784	8896

通过改变再入点初始速度方位角获得最大左侧机动弹道和最大右侧向机动弹道，不同纵程弹道侧向包络示意图见图2。记各落点对应的弹道侧向包络宽度为W_io，由下式(2)计算W_io：

W_io＝B_l-B_r(2)

考虑建模误差，取C_w倍W_io为实际侧向包络宽度W_i，这里取1.5倍W_io为实际侧向包络宽度W_i。计算得到不同纵程弹道侧向包络宽度见表4。

表4不同纵程弹道侧向包络宽度

对L_i(i＝1,2,3,4,…)和W_i(i＝1,2,3,4,…)进行拟合，可建立如下式(3)所示的弹道包络侧向宽度与纵程的数学模型，确定模型系数为确定模型系数为a_i：

$W=\underset{i=0}{\overset{term-1}{\&Sigma;}}{a}_i{L}^i---\left(3\right)$

具体为：W＝-1.064813L+0.000879L²-1.225047e^-7L³+4.607571e^-12L⁴，临近空间大范围机动弹道三维包络示意图见图3。

(2)建立换极坐标系，具体方法为：

引入一个换极坐标系：用表示换极坐标系中各物理量，用X表示一般坐标系中各物理量，按如下方式建立换极坐标系：

首先，定义一个再入大圆弧平面作为换极赤道平面：对目标点确定的情况，将再入点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面；对于目标点未确定的情况，根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面；

然后，基于换极赤道平面定义换极坐标系O_E为地心，轴沿再入点地心矢径方向，轴在换极赤道面内垂直于轴指向目标点方向，轴与轴、轴构成右手系，见图4。

(3)建立换极坐标系中飞行器动力学模型，具体方法为：

在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程，其飞行状态量为换极后的经度地心纬度航迹偏航角速度速度倾角和地心距

$\left\{\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \widehat{\mathrm{\&theta;}}\\ \frac{d\widehat{\mathrm{\&lambda;}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&sigma;}}}{\hat{r}\cos \widehat{\mathrm{\&phi;}}}\\ \frac{d\widehat{\mathrm{\&phi;}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&sigma;}}}{\hat{r}}\\ \frac{d\widehat{\mathrm{\&sigma;}}}{dt}=\frac{L\sin \mathrm{\&upsi;}}{\hat{V}\cos \widehat{\mathrm{\&theta;}}}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\&theta;}}\tan \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&sigma;}}-{\hat{g}}_{{\mathrm{\&omega;}}_e}\frac{\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&sigma;}}}{\hat{V}\cos \widehat{\mathrm{\&theta;}}}+C_{\mathrm{\&sigma;}}+{\tilde{C}}_{\mathrm{\&sigma;}}\\ \frac{d\widehat{\mathrm{\&theta;}}}{dt}=\frac{L}{\hat{V}}\cos \mathrm{\&upsi;}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\&theta;}}+\frac{{\hat{g}}_r\cos \widehat{\mathrm{\&theta;}}}{\hat{V}}+\frac{{\hat{g}}_{{\mathrm{\&omega;}}_e}}{\hat{V}}(\cos \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&phi;}}-\cos \widehat{\mathrm{\&sigma;}}\sin \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&phi;}})+C_{\mathrm{\&theta;}}+{\tilde{C}}_{\mathrm{\&theta;}}\\ \frac{d\hat{V}}{dt}=-D+{\hat{g}}_r\sin \widehat{\mathrm{\&theta;}}+{\hat{g}}_{{\mathrm{\&omega;}}_e}(\cos \widehat{\mathrm{\&sigma;}}\cos \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&phi;}}+\sin \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&phi;}})+{\tilde{C}}_V\end{array}\right.---\left(4\right)$

其中，C_σ、C_θ为哥氏加速度项，和为牵连加速度项，

$\left\{\begin{array}{c}{C}_{\mathrm{\&sigma;}}=(2{\mathrm{\&omega;}}_{ex}-2\tan \widehat{\mathrm{\&theta;}}({\mathrm{\&omega;}}_{ey}\sin \widehat{\mathrm{\&sigma;}}+{\mathrm{\&omega;}}_{ez}\cos \widehat{\mathrm{\&sigma;}}\left)\right)\\ {\tilde{C}}_{\mathrm{\&sigma;}}=-\frac{\hat{r}}{\hat{V}\cos \widehat{\mathrm{\&theta;}}}({\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ey}\cos \widehat{\mathrm{\&sigma;}}-{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ez}\sin \widehat{\mathrm{\&sigma;}})\\ C_{\mathrm{\&theta;}}=2({\mathrm{\&omega;}}_{ez}\sin \widehat{\mathrm{\&sigma;}}-{\mathrm{\&omega;}}_{ey}\cos \widehat{\mathrm{\&sigma;}})\\ {\tilde{C}}_{\mathrm{\&theta;}}=\frac{\hat{r}}{\hat{V}}\&lsqb;{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ey}\sin \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&sigma;}}+{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ez}\sin \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&sigma;}}+({\mathrm{\&omega;}}_{ey}^2+{\mathrm{\&omega;}}_{ez}^2)\cos \widehat{\mathrm{\&theta;}}\&rsqb;\\ {\tilde{C}}_V=\hat{r}\&lsqb;-{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ey}\cos \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&sigma;}}-{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ez}\cos \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&sigma;}}+({\mathrm{\&omega;}}_{ey}^2+{\mathrm{\&omega;}}_{ez}^2)\sin \widehat{\mathrm{\&theta;}}\&rsqb;\end{array}\right.---\left(5\right)$

其中，

$\{\begin{array}{c}{\mathrm{\&omega;}}_{ex}={\mathrm{\&omega;}}_e(\cos \widehat{\mathrm{\&lambda;}}\cos \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{cosA}}_p+\sin \widehat{\mathrm{\&lambda;}}\cos \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{sinA}}_p+\sin \widehat{\mathrm{\&phi;}}{\mathrm{sin\&phi;}}_p)\\ {\mathrm{\&omega;}}_{ey}={\mathrm{\&omega;}}_e(-\sin \widehat{\mathrm{\&lambda;}}{\mathrm{cos\&phi;}}_p{\mathrm{cosA}}_p+\cos \widehat{\mathrm{\&lambda;}}{\mathrm{cos\&phi;}}_p{\mathrm{sinA}}_p)\\ {\mathrm{\&omega;}}_{ez}={\mathrm{\&omega;}}_e(-\cos \widehat{\mathrm{\&lambda;}}\sin \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{cosA}}_p-\sin \widehat{\mathrm{\&lambda;}}\sin \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{sinA}}_p+\cos \widehat{\mathrm{\&phi;}}{\mathrm{sin\&phi;}}_p)\end{array}---\left(6\right)$

其中，ω_e为地球旋转加速度矢量，λ_p和φ_p为换极后极点P的经度和地心纬度，A_P为P的方位角，见图5。

根据换极坐标系定义，一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定义一致，

$r=\hat{r},\mathrm{\&theta;}=\widehat{\mathrm{\&theta;}},V=\hat{V}---\left(7\right)$

定义

$\begin{array}{c}\left[\begin{array}{ccc}{\mathrm{cos\&phi;}}_f{\mathrm{cos\&lambda;}}_f& {\mathrm{cos\&phi;}}_f{\mathrm{sin\&lambda;}}_f& {\mathrm{sin\&phi;}}_f\\ -{\mathrm{sin\&psi;}}_f{\mathrm{sin\&lambda;}}_f-{\mathrm{cos\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{cos\&lambda;}}_f& {\mathrm{sin\&psi;}}_f{\mathrm{cos\&lambda;}}_f-{\mathrm{cos\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{sin\&lambda;}}_f& {\mathrm{cos\&psi;}}_f{\mathrm{cos\&phi;}}_f\\ {\mathrm{cos\&psi;}}_f{\mathrm{sin\&lambda;}}_f-{\mathrm{sin\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{cos\&lambda;}}_f& -{\mathrm{cos\&psi;}}_f{\mathrm{cos\&lambda;}}_f-{\mathrm{sin\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{sin\&lambda;}}_f& {\mathrm{sin\&psi;}}_f{\mathrm{cos\&phi;}}_f\end{array}\right]\\ \stackrel{\mathrm{\&Delta;}}{=}\left[\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right]\end{array}---\left(8\right)$

其中，ψ_f为点F的方位角，

$\{\begin{array}{c}{G}_{11}\cos \mathrm{\&phi;}\cos \mathrm{\&lambda;}+G_{12}\cos \mathrm{\&phi;}\sin \mathrm{\&lambda;}+G_{13}\sin \mathrm{\&phi;}\stackrel{\mathrm{\&Delta;}}{=}{k}_1\\ G_{21}\cos \mathrm{\&phi;}\cos \mathrm{\&lambda;}+G_{22}\cos \mathrm{\&phi;}\sin \mathrm{\&lambda;}+G_{23}\sin \mathrm{\&phi;}\stackrel{\mathrm{\&Delta;}}{=}{k}_2\\ G_{31}\cos \mathrm{\&phi;}\cos \mathrm{\&lambda;}+G_{32}\cos \mathrm{\&phi;}\sin \mathrm{\&lambda;}+G_{33}\sin \mathrm{\&phi;}\stackrel{\mathrm{\&Delta;}}{=}{k}_3\end{array}---\left(9\right)$

$\{\begin{array}{c}{G}_{11}\cos \widehat{\mathrm{\&phi;}}\cos \widehat{\mathrm{\&lambda;}}+G_{21}\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&lambda;}}+G_{31}\sin \widehat{\mathrm{\&phi;}}\stackrel{\mathrm{\&Delta;}}{=}{\tilde{k}}_1\\ G_{12}\cos \widehat{\mathrm{\&phi;}}\cos \widehat{\mathrm{\&lambda;}}+G_{22}\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&lambda;}}+G_{32}\sin \widehat{\mathrm{\&phi;}}\stackrel{\mathrm{\&Delta;}}{=}{\tilde{k}}_2\\ G_{13}\cos \widehat{\mathrm{\&phi;}}\cos \widehat{\mathrm{\&lambda;}}+G_{23}\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&lambda;}}+G_{33}\sin \widehat{\mathrm{\&phi;}}\stackrel{\mathrm{\&Delta;}}{=}{\tilde{k}}_3\end{array}---\left(10\right)$

由一般坐标系中λ和φ确定换极坐标系中和的表达式为：

$\begin{array}{cc}\left\{\begin{array}{c}\cos \widehat{\mathrm{\&lambda;}}=k_1/\sqrt{{k_1}^2+{k_2}^2}\\ \sin \widehat{\mathrm{\&lambda;}}=k_2/\sqrt{{k_1}^2+{k_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \widehat{\mathrm{\&phi;}}=k_3\\ \cos \widehat{\mathrm{\&phi;}}=\sqrt{{k_1}^2+{k_2}^2}\end{array}\right.---\left(11\right)\end{array}$

由换极坐标系中λ和φ确定一般坐标系中λ和φ的表达式为：

$\begin{array}{cc}\left\{\begin{array}{c}\cos \mathrm{\&lambda;}={\tilde{k}}_1/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\\ \sin \mathrm{\&lambda;}={\tilde{k}}_2/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \mathrm{\&phi;}={\tilde{k}}_3\\ \cos \mathrm{\&phi;}=\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.\end{array}---\left(12\right)$

由一般坐标系中σ确定换极坐标系中的表达式为：

$\widehat{\mathrm{\&sigma;}}=\mathrm{\&sigma;}+\mathrm{\&eta;}---\left(13\right)$

其中，

$\{\begin{array}{c}\sin \mathrm{\&eta;}=\frac{\sin (\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_P){\mathrm{cos\&phi;}}_P}{\cos \widehat{\mathrm{\&phi;}}}\\ \cos \mathrm{\&eta;}=-\cos (A_P-\widehat{\mathrm{\&lambda;}})\cos (\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_P)+\sin (A_P-\widehat{\mathrm{\&lambda;}})\sin (\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_P){\mathrm{sin\&phi;}}_P\end{array}---\left(14\right);$

式(13)的位置关系见图6。

(4)计算换极坐标系弹道三维包络，具体方法为：

针对再入点经度λ_I＝100°、再入点地心纬度φ_I＝20°、落点经度λ_T＝165°、落点地心纬度φ_T＝50°的临近空间大范围机动弹道进行计算。根据步骤(3)所述坐标转换关系，计算得到换极弹道参数，由计算得

由式(15)计算换极弹道纵程

$\hat{L}=R_{e}arccos(sin{\widehat{\mathrm{\&phi;}}}_{I}sin{\widehat{\mathrm{\&phi;}}}_T+cos{\widehat{\mathrm{\&phi;}}}_{I}cos{\widehat{\mathrm{\&phi;}}}_{T}cos(\widehat{\mathrm{\&lambda;}}-{\widehat{\mathrm{\&lambda;}}}_I\left)\right)---\left(15\right)$

其中，R_e为地球半径；

由计算得将代入式(16)，可计算得到换极弹道侧向包络宽度

$\hat{W}=\underset{i=0}{\overset{term-1}{\&Sigma;}}{a}_i{\hat{L}}^i---\left(16\right)$

其中，模型系数a_i由步骤(1)确定；

由计算得将换极弹道侧向包络描述为长度为宽度为的矩形，其边界由式(17)确定：

$\begin{array}{cc}\{\begin{array}{c}{\hat{E}}_{\mathrm{\&lambda;}1}={\mathrm{\&lambda;}}_I\\ {\hat{E}}_{\mathrm{\&lambda;}2}=\frac{180\hat{L}}{{\mathrm{\&pi;R}}_e}\end{array},& \left\{\begin{array}{c}{\hat{E}}_{\mathrm{\&phi;}1}=\frac{90\hat{W}}{{\mathrm{\&pi;R}}_e}\\ {\hat{E}}_{\mathrm{\&phi;}2}=\frac{90\hat{W}}{{\mathrm{\&pi;R}}_e}\end{array}\right.\end{array}---\left(17\right)$

其中，为换极系侧向包络东向下边界，为换极系侧向包络东向上边界；为换极系侧向包络北向下边界，为换极系侧向包络北向上边界

侧向包络对应的经度范围Δλ及地心纬度范围Δφ由式(18)确定：

$\left\{\begin{array}{c}\mathrm{\&Delta;}\widehat{\mathrm{\&lambda;}}=\frac{180\hat{L}}{{\mathrm{\&pi;R}}_e}\\ \mathrm{\&Delta;}\widehat{\mathrm{\&phi;}}=\frac{180\hat{W}}{{\mathrm{\&pi;R}}_e}\end{array}\right.---\left(18\right)$

以弹道高度范围的C_r倍作为弹道包络的垂向范围

由计算可得，考虑机动和偏差，取 $\mathrm{\&Delta;}\hat{r}=50km.$

(5)换极坐标系全域空域剖分，具体方法为：

记换极系空间包络确定的空域为Ω，按式(19)计算Ω边界：

$\begin{array}{ccc}\{\begin{array}{c}{\hat{F}}_{r1}={\hat{r}}_I-\mathrm{\&Delta;}\hat{r}-2\&CenterDot;dr\\ {\hat{F}}_{r2}={\hat{r}}_I+2\&CenterDot;dr\end{array},& \{\begin{array}{c}{\hat{F}}_{\mathrm{\&lambda;}1}={\widehat{\mathrm{\&lambda;}}}_I-2\&CenterDot;d\mathrm{\&lambda;}\\ {\hat{F}}_{\mathrm{\&lambda;}2}={\widehat{\mathrm{\&lambda;}}}_I+\mathrm{\&Delta;}\widehat{\mathrm{\&lambda;}}+2\&CenterDot;d\mathrm{\&lambda;}\end{array},& \left\{\begin{array}{c}{\hat{F}}_{\mathrm{\&phi;}1}={\widehat{\mathrm{\&phi;}}}_I-\mathrm{\&Delta;}\widehat{\mathrm{\&phi;}}/2\\ {\hat{F}}_{\mathrm{\&phi;}2}={\widehat{\mathrm{\&phi;}}}_I+\mathrm{\&Delta;}\widehat{\mathrm{\&phi;}}/2\end{array}\right.---\left(19\right)\end{array}$

其中，为换极系Ω天向下边界，为换极系Ω天向上边界，为换极系Ω东向下边界，为换极系Ω东向上边界，为换极系Ω北向下边界，为换极系Ω北向上边界，dr、dλ和dφ分别为空域剖分间隔。

取计算得到的换极系空间包络空域范围见表5。

表5换极系空间包络空域范围

按照天向东向北向的间隔将由式(19)确定的空间包络空域Ω均匀剖分为q个互不重叠的子域Ω_e(e＝1,2,…,q)，记网格节点坐标为

$\mathrm{\&Omega;}=\underset{e=1}{\overset{q}{\&cup;}}{\mathrm{\&Omega;}}_e---\left(20\right)$

换极坐标系空间包络空域剖分节点示意图见图7。

(6)建立一般坐标系空间包络扰动引力重构模型，具体方法为：

根据换极系网格节点坐标由步骤(3)所述坐标转换关系计算一般坐标系网格节点坐标N(r_gi,λ_gj,φ_gk)(gi,gj,gk＝0,1,2,…)。一般坐标系空间包络空域剖分节点示意图见图8。

由球谐函数方法计算一般坐标系网格节点扰动引力位T_G，

$T_G=\frac{\mathrm{\&mu;}}{r_{gi}}\underset{n=2}{\overset{s}{\&Sigma;}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\&Sigma;}}({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{cosm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{sinm\&lambda;}}_{gj}){\stackrel{\&OverBar;}{P}}_{mn}\left({\mathrm{sin\&phi;}}_{gk}\right)---\left(21\right)$

其中，

${\stackrel{\&OverBar;}{C}}_{nm}^{*}=\left\{\begin{array}{ccc}0& n=2& m=0\\ {\stackrel{\&OverBar;}{C}}_{nm}& else& \end{array}\right.---\left(22\right)$

与T_G相应的引力即为扰动引力加速度δg，即

δg＝gradT_G(23)

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

$\left\{\begin{array}{c}{\mathrm{\&delta;g}}_R=\frac{\&part;T}{\&part;r_{gi}}=-\frac{\mathrm{\&mu;}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\&Sigma;}}(n+1){\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\&Sigma;}}({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{cosm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{sinm\&lambda;}}_{gj}){\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)\\ {\mathrm{\&delta;g}}_G=\frac{1}{r_{gi}{\mathrm{cos\&phi;}}_{gk}}\frac{\&part;T}{\&part;{\mathrm{\&lambda;}}_{gj}}=-\frac{1}{{\mathrm{cos\&phi;}}_{gk}}\frac{\mathrm{\&mu;}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\&Sigma;}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\&Sigma;}}m({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{sinm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{cosm\&lambda;}}_{gj}){\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)\\ {\mathrm{\&delta;g}}_G=\frac{1}{r_{gi}}\frac{\&part;T}{\&part;{\mathrm{\&phi;}}_{gk}}=\frac{\mathrm{\&mu;}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\&Sigma;}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\&Sigma;}}({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{cosm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{sinm\&lambda;}}_{gj})\frac{d{\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)}{{\mathrm{d\&beta;}}_{gk}}\end{array}\right.---\left(24\right)$

存储一般坐标系节点位置三分量和节点扰动引力三分量，完成空间包络扰动引力重构模型构建。空间包络节点扰动引力示意图见图9、图10和图11。

(7)当前弹道位置所在网格判断，具体方法为：

令当前弹道位置坐标为A(r,λ,φ)，当满足如式(25)所示的条件时，

$\{\begin{array}{c}{r}_{gi}\&le;r<r_{gi+1}\\ {\mathrm{\&lambda;}}_{gj}\&le;\mathrm{\&lambda;}<{\mathrm{\&lambda;}}_{gj+1}\\ {\mathrm{\&phi;}}_{gk}\&le;\mathrm{\&phi;}<{\mathrm{\&phi;}}_{gk+1}\end{array}---\left(25\right)$

判断A(r,λ,φ)位于由节点N(r_gi,λ_gj,φ_gk)、N(r_gi,λ_gj+1,φ_gk)、N(r_gi,λ_gj,φ_gk+1)、N(r_gi,λ_gj+1,φ_gk+1)、N(r_gi+1,λ_gj,φ_gk)、N(r_gi+1,λ_gj+1,φ_gk)、N(r_gi+1,λ_gj,φ_gk+1)、N(r_gi+1,λ_gj+1,φ_gk+1)所确定的网格中。

(8)建立计算网格内局域坐标系，具体方法为：

局部坐标系由半径r_l＝r_gi+dr/2的球面、经度λ_l＝λ_gj+dλ/2的子午面、纬度φ_l＝φ_gk+dφ/2的纬圈的交线组成，原点l(r_l,λ_l,φ_l)为三交线的交点，原点局部坐标为l(0,0,0)，ξ、η、ζ分别沿原点l的天向、北向和东向，则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为：

$\left\{\begin{array}{c}\mathrm{\&xi;}=r-r_l\\ \mathrm{\&eta;}=rcos\mathrm{\&phi;}(\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_l)\\ \mathrm{\&zeta;}=r(\mathrm{\&phi;}-{\mathrm{\&phi;}}_l)\end{array}\right.---\left(26\right)$

单元顶点A_i(r_i,λ_i,φ_i)的局部坐标A_i(ξ_i,η_i,ζ_i)为：

$\left\{\begin{array}{c}{\mathrm{\&xi;}}_i=r_i-r_i\\ {\mathrm{\&eta;}}_i=r_i{\mathrm{cos\&phi;}}_i({\mathrm{\&lambda;}}_i-{\mathrm{\&lambda;}}_l)\\ {\mathrm{\&zeta;}}_i=r_i({\mathrm{\&phi;}}_i-{\mathrm{\&phi;}}_l)\end{array}\right.---\left(27\right).$

(9)网格内部扰动引力网函数逼近计算，具体方法为：

基于网函数逼近理论进行单元内部的逼近计算。记某计算单元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(ξ,η,ζ)(i＝0,1,…,11)。网函数计算单元示意图见图12。

令L(ξ)、L(η)、L(ζ)分别为关于ξ、η、ζ的一次Lagrange插值算符，其插值基函数为：

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

l⁽³⁾＝L(ξ)L(η)+L(η)L(ζ)+L(ζ)L(ξ)-2L(ξ)L(η)L(ζ)

将l⁽³⁾作用于1-网函数，即可求得单元内任意一点A(ξ,η,ζ)的值δg(ξ,η,ζ)，

取节点扰动引力采用1080阶球谐函数进行赋值。以1080阶球谐函数计算得到的扰动引力值作为近似真值，对比分析逼近方法在不同情况下的赋值情况。

针对纵程为6700km，初始方位角为90°的滑翔弹道进行计算。弹道再入点和落点分别为：①平原地区λ_{I_A1}＝115°，φ_{I_A1}＝35°，λ_{T_A1}＝180°，φ_{T_A1}＝16.5°；②丘陵地区λ_{I_A2}＝110°，φ_{I_A2}＝27°，λ_{T_A2}＝173°，φ_{T_A2}＝13°；③特大山区λ_{I_A3}＝80°，φ_{I_A3}＝42°，λ_{T_A3}＝147°，φ_{T_A3}＝19.3°。沿不同区域弹道扰动引力逼近结果示意图见图13，逼近误差见图14和表6。

表6不同区域弹道扰动引力逼近精度统计

针对纵程为7000km，弹道起点位于特大山区(λ_I＝80°，φ_I＝42°)，初始方位角分别为0°、180°和270°的滑翔弹道进行计算。滑翔弹道终点分别为：①正北λ_{T_σ1}＝80°，φ_{T_σ1}＝102.3°；②正南λ_{T_σ2}＝80°，φ_{T_σ2}＝-18.3°；③正西λ_{T_σ3}＝12.5°，φ_{T_σ3}＝20°。沿不同射向弹道扰动引力逼近结果示意图见图15，逼近误差见图16和表7。

表7不同射向弹道扰动引力逼近精度统计

针对弹道起点位于平原地区(λ_I＝115°，φ_I＝35°)，初始方位角为90°，射程分别为6000km、7000km和8000km的滑翔弹道进行计算。滑翔弹道终点分别为：①射程6000km弹道λ_{T_r1}＝174.3°，φ_{T_r1}＝19.8°；②射程7000km弹道λ_{T_r2}＝182.2°，φ_{T_r2}＝15°；③射程8000km弹道λ_{T_r3}＝189.9°，φ_{T_r3}＝10°。沿不同射程弹道扰动引力逼近结果示意图见图17，逼近误差见图18和表8。

表8不同射程弹道扰动引力逼近精度统计

针对弹道起点为λ_I＝0°，φ_I＝0°，初始方位角为i×10°(i＝0,1,2,…,35)的7000km射程滑翔弹道进行计算。由扰动引力引起的弹道终端状态偏差见图19。分别采用逼近方法和1080阶球谐函数方法求解弹道计算中的扰动引力，视两者落点偏差之差为逼近误差引起的落点偏差。计算得到逼近误差引起的弹道终端状态偏差见图20。

由图13～图20及表6～表8可见，本发明提出的方法能够适用于不同滑翔弹道的计算。在的网格划分下，沿滑翔弹道的平均逼近误差可控制在0.1mgal～1mgal量级。由逼近误差导致的滑翔弹道终端位置偏差在10m量级，使弹道精度提高约1e-3个量级。

表9采用不同扰动引力计算方法所需的弹道计算时间

计算方法	计算时间(s)
		网函数逼近方法	61.257
1080阶球谐函数方法	125970.109
		Stokes方法	104653.258

采用不同扰动引力计算方法所需的弹道计算时间见表9。由表9可见，采用本发明方法进行弹道计算中的扰动引力计算，单条弹道的计算时间远小于采用传统扰动引力计算方法所需的弹道计算时间。

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

1、由本发明确定的方法建立沿临近空间大范围机动弹道空间包络的扰动引力重构模型，对位于不同区域，具有不同射程和不同射向的弹道均有较好的适应性，弹道积分无奇点。

2、通过本发明计算得到的扰动引力具有较高的逼近精度。在10km/1°/1°的网格划分下，可使沿弹道的平均扰动引力逼近误差控制在0.1mgal～1mgal量级，由逼近误差导致的滑翔弹道终端位置偏差控制在10m量级，弹道精度提高约1e-3个量级。

3、本方法具有计算速度快的特征。采用本发明方法进行弹道计算中的扰动引力计算，单条弹道的计算时间远小于采用传统扰动引力计算方法所需的弹道计算时间。

上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。

Claims (2)

1.一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法，其特征在于：包括以下步骤：

(1)建立弹道包络侧向宽度与纵程的数学模型；

(2)建立换极坐标系；

(3)建立换极坐标系中飞行器动力学模型；

(4)计算换极坐标系弹道三维包络；

(5)换极坐标系全域空域剖分；

(6)建立一般坐标系空间包络扰动引力重构模型；

(7)当前弹道位置所在网格判断；

(8)建立计算网格内局域坐标系；

(9)网格内部扰动引力网函数逼近计算。

2.根据权利要求1所述的沿临近空间大范围机动弹道空间包络的扰动引力逼近方法，其特征在于：

所述步骤(1)的方法为：

令弹道再入点为I，其经度为λ_I、地心纬度为φ_I，弹道落点为T，其经度为λ_T、地心纬度为φ_T；以由再入点和落点确定的再入大圆弧平面为对称面，称沿射向方向对称面以左弹道为左侧机动弹道，沿射向方向对称面以右弹道为右侧机动弹道；记最大左侧机动弹道距对称面的最大距离为左边界B_l，记最大右侧机动弹道距对称面的最大距离为右边界B_r；

根据给定的飞行器再入点飞行状态参数、弹道终端飞行状态参数、飞行过程约束条件及终端约束条件，针对再入点为(λ_I,φ_I)＝(0°,0°)、落点为(λ_ki,0°)(i＝1,2,3,4,…)的低空大范围机动弹道进行弹道计算，记各落点对应的弹道纵程为L_i，由下式(1)计算L_i：

L_i＝R_earccos(sinφ_Isinφ_T+cosφ_Icosφ_Tcos(λ_T-λ_I))(1)

其中，R_e为地球半径；

通过改变再入点初始速度方位角获得最大左侧机动弹道和最大右侧向机动弹道，记各落点对应的弹道侧向包络宽度为W_io，由下式(2)计算W_io：

W_io＝B_l-B_r(2)

取C_w倍W_io为实际侧向包络宽度W_i，对L_i(i＝1,2,3,4,…)和W_i(i＝1,2,3,4,…)进行拟合，可建立如下式(3)所示的弹道包络侧向宽度与纵程的数学模型，确定模型系数为a_i：

$W=\underset{i=0}{\overset{term-1}{\&Sigma;}}{a}_i{L}^i---\left(3\right);$

所述步骤(2)的方法为：

引入一个换极坐标系：用表示换极坐标系中各物理量，用X表示一般坐标系中各物理量，按如下方式建立换极坐标系：

首先，定义一个再入大圆弧平面作为换极赤道平面：对目标点确定的情况，将再入点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面；对于目标点未确定的情况，根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面；

然后，基于换极赤道平面定义换极坐标系O_E为地心，轴沿再入点地心矢径方向，轴在换极赤道面内垂直于轴指向目标点方向，轴与轴、轴构成右手系；

所述步骤(3)的方法为：

在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程，其飞行状态量为换极后的经度地心纬度航迹偏航角速度速度倾角和地心距

$\{\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \widehat{\mathrm{\&theta;}}\\ \frac{d\widehat{\mathrm{\&lambda;}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&sigma;}}}{\hat{r}\cos \widehat{\mathrm{\&phi;}}}\\ \frac{d\widehat{\mathrm{\&phi;}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&sigma;}}}{\hat{r}}\\ \frac{d\widehat{\mathrm{\&sigma;}}}{dt}=\frac{L\sin \mathrm{\&upsi;}}{\hat{V}\cos \widehat{\mathrm{\&theta;}}}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\&theta;}}\tan \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&sigma;}}-{\hat{g}}_{{\mathrm{\&omega;}}_e}\frac{\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&sigma;}}}{\hat{V}\cos \widehat{\mathrm{\&theta;}}}+C_{\mathrm{\&sigma;}}+{\tilde{C}}_{\mathrm{\&sigma;}}\\ \frac{d\widehat{\mathrm{\&theta;}}}{dt}=\frac{L}{\hat{V}}\cos \mathrm{\&upsi;}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\&theta;}}+\frac{{\hat{g}}_r\cos \widehat{\mathrm{\&theta;}}}{\hat{V}}+\frac{{\hat{g}}_{{\mathrm{\&omega;}}_e}}{\hat{V}}(\cos \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&phi;}}-\cos \widehat{\mathrm{\&sigma;}}\sin \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&phi;}})+C_{\mathrm{\&theta;}}+{\tilde{C}}_{\mathrm{\&theta;}}\\ \frac{d\hat{V}}{dt}=-D+{\hat{g}}_r\sin \widehat{\mathrm{\&theta;}}+{\widehat{\mathrm{\&theta;}}}_{{\mathrm{\&omega;}}_e}(\cos \widehat{\mathrm{\&sigma;}}\cos \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&phi;}}+\sin \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&phi;}})+{\tilde{C}}_V\end{array}---\left(4\right)$

其中，C_σ、C_θ为哥氏加速度项，和为牵连加速度项，

$\left\{\begin{array}{c}{C}_{\mathrm{\&sigma;}}=(2{\mathrm{\&omega;}}_{ex}-2\tan \widehat{\mathrm{\&theta;}}({\mathrm{\&omega;}}_{ey}\sin \widehat{\mathrm{\&sigma;}}+{\mathrm{\&omega;}}_{ez}\cos \widehat{\mathrm{\&sigma;}}))\\ {\tilde{C}}_{\mathrm{\&sigma;}}=-\frac{\hat{r}}{\hat{V}\cos \widehat{\mathrm{\&theta;}}}({\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ey}\cos \widehat{\mathrm{\&sigma;}}-{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ez}\sin \widehat{\mathrm{\&sigma;}})\\ C_{\mathrm{\&theta;}}=2({\mathrm{\&omega;}}_{ez}\sin \widehat{\mathrm{\&sigma;}}-{\mathrm{\&omega;}}_{ey}\cos \widehat{\mathrm{\&sigma;}})\\ {\tilde{C}}_{\mathrm{\&theta;}}=\frac{\hat{r}}{\hat{V}}\&lsqb;{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ey}\sin \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&sigma;}}+{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ez}\sin \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&sigma;}}+({\mathrm{\&omega;}}_{ey}^2+{\mathrm{\&omega;}}_{ez}^2)\cos \widehat{\mathrm{\&theta;}}\&rsqb;\\ {\tilde{C}}_V=\hat{r}\&lsqb;-{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ey}\cos \widehat{\mathrm{\&theta;}}\sin \widehat{\mathrm{\&sigma;}}-{\mathrm{\&omega;}}_{ex}{\mathrm{\&omega;}}_{ez}\cos \widehat{\mathrm{\&theta;}}\cos \widehat{\mathrm{\&sigma;}}+({\mathrm{\&omega;}}_{ey}^2+{\mathrm{\&omega;}}_{ez}^2)\sin \widehat{\mathrm{\&theta;}}\&rsqb;\end{array}\right.---\left(5\right)$

其中，

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_{ex}={\mathrm{\&omega;}}_e(\cos \widehat{\mathrm{\&lambda;}}\cos \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{cosA}}_p+\sin \widehat{\mathrm{\&lambda;}}\cos \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{sinA}}_p+\sin \widehat{\mathrm{\&phi;}}{\mathrm{sin\&phi;}}_p)\\ {\mathrm{\&omega;}}_{ey}={\mathrm{\&omega;}}_e(-\sin \widehat{\mathrm{\&lambda;}}{\mathrm{cos\&phi;}}_p{\mathrm{cosA}}_p+\cos \widehat{\mathrm{\&lambda;}}{\mathrm{cos\&phi;}}_p{\mathrm{sinA}}_p)\\ {\mathrm{\&omega;}}_{ez}={\mathrm{\&omega;}}_e(-\cos \widehat{\mathrm{\&lambda;}}\sin \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{cosA}}_p-\sin \widehat{\mathrm{\&lambda;}}\sin \widehat{\mathrm{\&phi;}}{\mathrm{cos\&phi;}}_p{\mathrm{sinA}}_p+\cos \widehat{\mathrm{\&phi;}}{\mathrm{sin\&phi;}}_p)\end{array}\right.---\left(6\right)$

其中，ω_e为地球旋转加速度矢量，λ_p和φ_p为换极后极点P的经度和地心纬度，A_P为P的方位角；

根据换极坐标系定义，一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定义一致，

$r=\hat{r},\mathrm{\&theta;}=\widehat{\mathrm{\&theta;}},V=\hat{V}---\left(7\right)$

定义

$\begin{array}{c}\left[\begin{array}{ccc}{\mathrm{cos\&phi;}}_f{\mathrm{cos\&lambda;}}_f& {\mathrm{cos\&phi;}}_f{\mathrm{sin\&lambda;}}_f& {\mathrm{sin\&phi;}}_f\\ -{\mathrm{sin\&psi;}}_f{\mathrm{sin\&lambda;}}_f-{\mathrm{cos\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{cos\&lambda;}}_f& {\mathrm{sin\&psi;}}_f{\mathrm{cos\&lambda;}}_f-{\mathrm{cos\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{sin\&lambda;}}_f& {\mathrm{cos\&psi;}}_f{\mathrm{cos\&phi;}}_f\\ {\mathrm{cos\&psi;}}_f{\mathrm{sin\&lambda;}}_f-{\mathrm{sin\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{cos\&lambda;}}_f& -{\mathrm{cos\&psi;}}_f{\mathrm{cos\&lambda;}}_f-{\mathrm{sin\&psi;}}_f{\mathrm{sin\&phi;}}_f{\mathrm{sin\&lambda;}}_f& {\mathrm{sin\&psi;}}_f{\mathrm{cos\&phi;}}_f\end{array}\right]\\ \stackrel{\mathrm{\&Delta;}}{=}\left[\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{\mathrm{21}}& G_{\mathrm{22}}& G_{23}\\ G_{\mathrm{31}}& G_{\mathrm{32}}& G_{33}\end{array}\right]\end{array}---\left(8\right)$

其中，ψ_f为点F的方位角，

$\left\{\begin{array}{c}{G}_{11}\cos \mathrm{\&phi;}\cos \mathrm{\&lambda;}+G_{12}\cos \mathrm{\&phi;}\sin \mathrm{\&lambda;}+G_{13}\sin \mathrm{\&phi;}\stackrel{\mathrm{\&Delta;}}{=}{k}_1\\ G_{21}\cos \mathrm{\&phi;}\cos \mathrm{\&lambda;}+G_{22}\cos \mathrm{\&phi;}\sin \mathrm{\&lambda;}+G_{23}\sin \mathrm{\&phi;}\stackrel{\mathrm{\&Delta;}}{=}{k}_2\\ G_{31}\cos \mathrm{\&phi;}\cos \mathrm{\&lambda;}+G_{32}\cos \mathrm{\&phi;}\sin \mathrm{\&lambda;}+G_{33}\sin \mathrm{\&phi;}\stackrel{\mathrm{\&Delta;}}{=}{k}_3\end{array}\right.---\left(9\right)$

$\left\{\begin{array}{c}{G}_{11}\cos \widehat{\mathrm{\&phi;}}\cos \widehat{\mathrm{\&lambda;}}+G_{21}\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&lambda;}}+G_{31}\sin \widehat{\mathrm{\&phi;}}\stackrel{\mathrm{\&Delta;}}{=}{\tilde{k}}_1\\ G_{12}\cos \widehat{\mathrm{\&phi;}}\cos \widehat{\mathrm{\&lambda;}}+G_{22}\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&lambda;}}+G_{32}\sin \widehat{\mathrm{\&phi;}}\stackrel{\mathrm{\&Delta;}}{=}{\tilde{k}}_2\\ G_{13}\cos \widehat{\mathrm{\&phi;}}\cos \widehat{\mathrm{\&lambda;}}+G_{23}\cos \widehat{\mathrm{\&phi;}}\sin \widehat{\mathrm{\&lambda;}}+G_{33}\sin \widehat{\mathrm{\&phi;}}\stackrel{\mathrm{\&Delta;}}{=}{\tilde{k}}_3\end{array}\right.---\left(10\right)$

由一般坐标系中λ和φ确定换极坐标系中和的表达式为：

$\begin{array}{cc}\left\{\begin{array}{c}\cos \widehat{\mathrm{\&lambda;}}=k_1/\sqrt{{k_1}^2+{k_1}^2}\\ sin\widehat{\mathrm{\&lambda;}}=k_2/\sqrt{{k_1}^2+{k_1}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \widehat{\mathrm{\&phi;}}=k_3\\ \cos \widehat{\mathrm{\&phi;}}=\sqrt{{k_1}^2+{k_1}^2}\end{array}\right.\end{array}---\left(11\right)$

由换极坐标系中λ和φ确定一般坐标系中λ和φ的表达式为：

$\begin{array}{cc}\left\{\begin{array}{c}\cos \mathrm{\&lambda;}={\tilde{k}}_1/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\\ sin\mathrm{\&lambda;}={\tilde{k}}_2/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \mathrm{\&phi;}={\tilde{k}}_3\\ \cos \mathrm{\&phi;}=\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.\end{array}---\left(12\right)$

由一般坐标系中σ确定换极坐标系中的表达式为：

$\widehat{\mathrm{\&sigma;}}=\mathrm{\&sigma;}+\mathrm{\&eta;}---\left(13\right)$

其中，

$\left\{\begin{array}{c}\sin \mathrm{\&eta;}=\frac{\sin (\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_P){\mathrm{cos\&phi;}}_P}{\cos \widehat{\mathrm{\&phi;}}}\\ \cos \mathrm{\&eta;}=-\cos (A_P-\widehat{\mathrm{\&lambda;}})\cos (\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_P)+\sin (A_P-\widehat{\mathrm{\&lambda;}})\sin (\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_P){\mathrm{sin\&phi;}}_P\end{array}\right.---\left(14\right);$

所述步骤(4)的方法为：

根据步骤(3)所述坐标转换关系，根据弹道再入点经度λ_I、再入点地心纬度φ_I、落点经度λ_T、落点地心纬度φ_T计算得到换极弹道参数，由计算得 记换极弹道落点经度为

由式(15)计算换极弹道纵程

$\hat{L}=R_{e}arccos(sin{\widehat{\mathrm{\&phi;}}}_I\sin {\widehat{\mathrm{\&phi;}}}_T+cos{\widehat{\mathrm{\&phi;}}}_{I}cos{\widehat{\mathrm{\&phi;}}}_{T}cos({\widehat{\mathrm{\&lambda;}}}_T-{\widehat{\mathrm{\&lambda;}}}_I\left)\right)---\left(15\right)$

其中，R_e为地球半径；

将代入式(16)，可计算得到换极弹道侧向包络宽度

$\hat{W}=\underset{i=0}{\overset{term-1}{\&Sigma;}}{a}_i{\hat{L}}^i---\left(16\right)$

其中，模型系数a_i由步骤(1)确定；

将换极弹道侧向包络描述为长度为宽度为的矩形，其边界由式(17)确定：

$\{\begin{array}{c}{\hat{E}}_{\mathrm{\&lambda;}1}={\mathrm{\&lambda;}}_I\\ {\hat{E}}_{\mathrm{\&lambda;}2}=\frac{180\hat{L}}{{\mathrm{\&pi;R}}_e}\end{array},\left\{\begin{array}{c}{\hat{E}}_{\mathrm{\&phi;}1}=-\frac{90\hat{W}}{{\mathrm{\&pi;R}}_e}\\ {\hat{E}}_{\mathrm{\&phi;}2}=\frac{90\hat{W}}{{\mathrm{\&pi;R}}_e}\end{array}\right.---\left(17\right)$

其中，为换极系侧向包络东向下边界，为换极系侧向包络东向上边界；为换极系侧向包络北向下边界，为换极系侧向包络北向上边界；

侧向包络对应的经度范围Δλ及地心纬度范围Δφ由式(18)确定：

$\left\{\begin{array}{c}\mathrm{\&Delta;}\widehat{\mathrm{\&lambda;}}=\frac{180\hat{L}}{{\mathrm{\&pi;R}}_e}\\ \mathrm{\&Delta;}\widehat{\mathrm{\&phi;}}=\frac{180\hat{W}}{{\mathrm{\&pi;R}}_e}\end{array}\right.---\left(18\right)$

以弹道高度范围的C_r倍作为弹道包络的垂向范围

所述步骤(5)的方法为：

记换极系空间包络确定的空域为Ω，按式(19)计算Ω边界：

$\left\{\begin{array}{c}{\hat{F}}_{r1}={\hat{r}}_I-\mathrm{\&Delta;}\hat{r}-2\&CenterDot;dr\\ {\hat{F}}_{r2}={\hat{r}}_I+2\&CenterDot;dr\end{array}\right.,\{\begin{array}{c}{\hat{F}}_{\mathrm{\&lambda;}1}={\widehat{\mathrm{\&lambda;}}}_I-2\&CenterDot;d\mathrm{\&lambda;}\\ {\hat{F}}_{\mathrm{\&lambda;}2}={\widehat{\mathrm{\&lambda;}}}_I+\mathrm{\&Delta;}\widehat{\mathrm{\&lambda;}}+2\&CenterDot;d\mathrm{\&lambda;}\end{array},\left\{\begin{array}{c}{\hat{F}}_{\mathrm{\&phi;}1}={\widehat{\mathrm{\&phi;}}}_I-\mathrm{\&Delta;}\widehat{\mathrm{\&phi;}}\mathrm{/2}\\ {\hat{F}}_{\mathrm{\&phi;}2}={\widehat{\mathrm{\&phi;}}}_I+\mathrm{\&Delta;}\widehat{\mathrm{\&phi;}}\mathrm{/2}\end{array}\right.---\left(19\right)$

其中，为换极系Ω天向下边界，为换极系Ω天向上边界，为换极系Ω东向下边界，为换极系Ω东向上边界，为换极系Ω北向下边界，为换极系Ω北向上边界，dr、dλ和dφ分别为空域剖分间隔；

按照天向东向北向的间隔将由式(19)确定的空间包络空域Ω均匀剖分为q个互不重叠的子域Ω_e(e＝1,2,…,q)，记网格节点坐标为 $\hat{N}({\hat{r}}_{gi},{\widehat{\mathrm{\&lambda;}}}_{gj},{\widehat{\mathrm{\&phi;}}}_{gk})(gi,gj,gk=0,1,2,...),$

$\mathrm{\&Omega;}=\underset{e=1}{\overset{q}{\&cup;}}{\mathrm{\&Omega;}}_e---\left(20\right);$

所述步骤(6)的方法为：

根据换极系网格节点坐标由步骤(3)所述坐标转换关系计算一般坐标系网格节点坐标 $N(r_{gi},{\mathrm{\&lambda;}}_{gj},{\mathrm{\&phi;}}_{gk})(gi,gj,gk=0,1,2,...);$

由球谐函数方法计算一般坐标系网格节点扰动引力位T_G，

$T_G=\frac{\mathrm{\&mu;}}{r_{gi}}\underset{n=2}{\overset{s}{\mathrm{\&Sigma;}}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{cosm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{sinm\&lambda;}}_{gj}){\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)---\left(21\right)$

其中，

${\stackrel{\&OverBar;}{C}}_{nm}^{*}=\left\{\begin{array}{ccc}0& n=2& m=0\\ {\stackrel{\&OverBar;}{C}}_{nm}& else& \end{array}\right.---\left(22\right)$

与T_G相应的引力即为扰动引力加速度δg，即

δg＝gradT_G(23)

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

$\left\{\begin{array}{c}{\mathrm{\&delta;g}}_R=\frac{\&part;T}{\&part;r_{gi}}=-\frac{\mathrm{\&mu;}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\mathrm{\&Sigma;}}}(n+1){\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{cosm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{sinm\&lambda;}}_{gj}){\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)\\ {\mathrm{\&delta;g}}_G=\frac{1}{r_{gi}{\mathrm{cos\&phi;}}_{gk}}\frac{\&part;T}{\&part;{\mathrm{\&lambda;}}_{gj}}=-\frac{1}{{\mathrm{cos\&phi;}}_{gk}}\frac{\mathrm{\&mu;}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\mathrm{\&Sigma;}}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}m({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{sinm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{cosm\&lambda;}}_{gj}){\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)\\ {\mathrm{\&delta;g}}_G=\frac{1}{r_{gi}}\frac{\&part;T}{\&part;{\mathrm{\&phi;}}_{gk}}=\frac{\mathrm{\&mu;}}{{r_{gi}}^2}\underset{n=2}{\overset{s}{\mathrm{\&Sigma;}}}{\left(\frac{a_e}{r_{gi}}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}m({\stackrel{\&OverBar;}{C}}_{nm}^{*}{\mathrm{cosm\&lambda;}}_{gj}+{\stackrel{\&OverBar;}{S}}_{nm}{\mathrm{sinm\&lambda;}}_{gj})\frac{d{\stackrel{\&OverBar;}{P}}_{nm}\left({\mathrm{sin\&phi;}}_{gk}\right)}{{\mathrm{d\&phi;}}_{gk}}\end{array}\right.---\left(24\right)$

存储一般坐标系节点位置三分量和节点扰动引力三分量，完成空间包络扰动引力重构模型构建；

所述步骤(7)的方法为：

令当前弹道位置坐标为A(r,λ,φ)，当满足如式(25)所示的条件时，

$\left\{\begin{array}{c}{r}_{gi}\&le;r<r_{gi+1}\\ {\mathrm{\&lambda;}}_{gj}\&le;\mathrm{\&lambda;}<{\mathrm{\&lambda;}}_{gj+1}\\ {\mathrm{\&phi;}}_{gk}\&le;\mathrm{\&phi;}<{\mathrm{\&phi;}}_{gk+1}\end{array}\right.---\left(25\right)$

判断A(r,λ,φ)位于由节点N(r_gi,λ_gj,φ_gk)、N(r_gi,λ_gj+1,φ_gk)、N(r_gi,λ_gj,φ_gk+1)、N(r_gi,λ_gj+1,φ_gk+1)、N(r_gi+1,λ_gj,φ_gk)、N(r_gi+1,λ_gj+1,φ_gk)、N(r_gi+1,λ_gj,φ_gk+1)、N(r_gi+1,λ_gj+1,φ_gk+1)所确定的网格中；

所述步骤(8)的方法为：

局部坐标系由半径r_l＝r_gi+dr/2的球面、经度λ_l＝λ_gj+dλ/2的子午面、纬度φ_l＝φ_gk+dφ/2的纬圈的交线组成，原点l(r_l,λ_l,φ_l)为三交线的交点，原点局部坐标为l(0,0,0)，ξ、η、ζ分别沿原点l的天向、北向和东向，则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为：

$\left\{\begin{array}{c}\mathrm{\&xi;}=r-r_l\\ \mathrm{\&eta;}=rcos\mathrm{\&phi;}(\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_l)\\ \mathrm{\&zeta;}=r(\mathrm{\&phi;}-{\mathrm{\&phi;}}_l)\end{array}\right.---\left(26\right)$

单元顶点A_i(r_i,λ_i,φ_i)的局部坐标A_i(ξ_i,η_i,ζ_i)为：

$\left\{\begin{array}{c}{\mathrm{\&xi;}}_i=r_i-r_l\\ {\mathrm{\&eta;}}_i=r_i{\mathrm{cos\&phi;}}_i({\mathrm{\&lambda;}}_i-{\mathrm{\&lambda;}}_l)\\ {\mathrm{\&zeta;}}_i=r_i({\mathrm{\&phi;}}_i-{\mathrm{\&phi;}}_l)\end{array}\right.---\left(27\right);$

所述步骤(9)的方法为：

记某计算单元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(ξ,η,ζ)(i＝0,1,…,11)；令L(ξ)、L(η)、L(ζ)分别为关于ξ、η、ζ的一次Lagrange插值算符，其插值基函数为：

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

l⁽³⁾＝L(ξ)L(η)+L(η)L(ζ)+L(ζ)L(ξ)-2L(ξ)L(η)L(ζ)

将l⁽³⁾作用于1-网函数，即可求得单元内任意一点A(ξ,η,ζ)的值δg(ξ,η,ζ)，