CN105760573A - 沿临近空间大范围机动弹道的扰动引力延拓逼近方法 - Google Patents

Abstract

本发明公开了一种沿临近空间大范围机动弹道的扰动引力延拓逼近方法，包括以下步骤：第一步，建立弹道包络侧向宽度与纵程的数学模型；第二步，建立换极坐标系；第三步，建立换极坐标系中飞行器动力学模型；第四步，计算换极坐标系弹道三维包络；第五步，换极坐标系全域空域剖分；第六步，建立一般坐标系扰动引力全域重构模型；第七步，建立沿弹道的扰动引力局域重构模型；第八步，建立子域局部曲线坐标系；第九步，飞行过程中扰动引力延拓逼近计算。本发明具有逼近精度高、计算速度快、数据存储量小的特征，能适应沿任意临近空间大范围机动弹道扰动引力快速计算的要求，且具备快速响应飞行任务临时变更的能力。

Description

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

技术领域

本发明属于飞行器动力学建模领域，特别涉及沿临近空间大范围机动弹道的扰动引力延拓逼近方法。

背景技术

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

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

发明内容

本发明针对沿临近空间大范围机动弹道的扰动引力计算问题，首次提出一种基于延拓逼近理论的高精度快速计算方法。提出的方法具有逼近精度高、计算速度快、数据存储量小的特征，能适应沿任意临近空间大范围机动弹道扰动引力快速计算的要求，且具备快速响应飞行任务临时变更的能力。

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

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

A1：建立弹道包络侧向宽度与纵程的数学模型

令弹道再入点为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)$

A2：建立换极坐标系

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

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

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

A3：建立换极坐标系中飞行器动力学模型

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

$\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 \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)$

其中，

$\left\{\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}\right.---\left(6\right)$

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

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

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

定义

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

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

$\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.\end{array}---\left(11\right)$

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

$\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)$

其中，

$\left\{\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}\right.---\left(14\right)$

A4：计算换极坐标系弹道三维包络

根据第三步所述坐标转换关系，可根据弹道再入点经度λ_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由第一步确定，

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

$\left\{\begin{array}{c}{\hat{E}}_{\mathrm{\λ}1}={\mathrm{\λ}}_I\\ {\hat{E}}_{\mathrm{\λ}2}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\end{array}\right.,\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.---\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倍作为弹道包络的垂向范围

A5：换极坐标系全域空域剖分

记扰动引力全域重构模型为Field模型，其确定的空域为Ω，为保证首个和末个延拓点有意义，按式(19)计算换极系Ω边界，

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

按照天向dr、东向dλ、北向dφ的间隔将由式(19)确定的换极系全域空域Ω均匀剖分为q个互不重叠的子域Ω_e(e＝1,2,…,q)，记网格节点坐标为 $\hat{N}({\hat{r}}_G,{\widehat{\mathrm{\λ}}}_G,{\widehat{\mathrm{\φ}}}_G),$

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

A6：建立一般坐标系扰动引力全域重构模型

根据换极系网格节点坐标由第三步所述坐标转换关系计算一般坐标系网格节点坐标N(r_G,λ_G,φ_G)，

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

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

其中，

${\stackrel{\‾}{C}}_{nm}^{*}=\left\{\begin{array}{cc}0& \begin{array}{cc}n=2& m=0\end{array}\\ {\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_G}=-\frac{\mathrm{\μ}}{{r_G}^2}{\underset{n=2}{\overset{s}{\Σ}}(n+1)\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_G+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G{\mathrm{cos\φ}}_G}\frac{\∂T}{\∂{\mathrm{\λ}}_G}=-\frac{1}{{\mathrm{cos\φ}}_G}\frac{\mathrm{\μ}}{{r_G}^2}{\underset{n=2}{\overset{s}{\Σ}}\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}m({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{sinm\λ}}_G+{\stackrel{\‾}{S}}_{nm}{\mathrm{cosm\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G}\frac{\∂T}{\∂{\mathrm{\φ}}_G}=\frac{\mathrm{\μ}}{{r_G}^2}{\underset{n=2}{\overset{s}{\Σ}}\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_G+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_G)\frac{d{\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)}{{\mathrm{d\φ}}_G}\end{array}\right.---\left(24\right)$

存储一般坐标系节点位置三分量和节点扰动引力三分量，完成扰动引力全域重构模型(Field模型)构建；

A7：建立沿弹道的扰动引力局域重构模型

记沿弹道的扰动引力局域重构模型为Channel模型，根据飞行任务进行弹道规划，计算参考弹道，判断参考弹道历经的三层延拓网格，从Field模型数据中读取历经延拓网格的节点位置三分量及节点扰动引力三分量，完成扰动引力局域重构模型(Channel模型)构建；

A8：建立子域局部曲线坐标系

令子域Ω_e及其相邻网格确定的空域为延拓域Ω′_e，有记Ω_e包含的节点数为r，Ω′_e包含的节点数为s，有s＞r，

建立子域局部曲线坐标系，在局部坐标系中描述单元节点和计算点坐标，子域Ω_e的局部曲线坐标系由半径r_l＝r_G+dr/2的球面、经度λ_l＝λ_G+dλ/2的子午面、纬度φ_l＝φ_G+dφ/2的纬圈的交线组成，原点l(r_l,λ_l,φ_l)为三交线的交点，原点局部坐标为l(0,0,0)，ξ、η、ζ分别沿原点l的天向、北向和东向，则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为，

$\left\{\begin{array}{c}\mathrm{\ξ}=r-r_l\\ \mathrm{\η}=r\cos \mathrm{\φ}(\mathrm{\λ}-{\mathrm{\λ}}_l)\\ \mathrm{\ζ}=r(\mathrm{\φ}-{\mathrm{\φ}}_l)\end{array}\right.---\left(25\right)$

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

$\left\{\begin{array}{c}{\mathrm{\ξ}}_i=r_i-r_l\\ {\mathrm{\η}}_i=r_i{\mathrm{cos\φ}}_i({\mathrm{\λ}}_i-{\mathrm{\λ}}_l)\\ {\mathrm{\ζ}}_i=r_i({\mathrm{\φ}}_i-{\mathrm{\φ}}_l)\end{array}\right.---\left(26\right)$

A9：飞行过程中扰动引力延拓逼近计算

令延拓域Ω′_e节点扰动引力值为u_i，在子域局部曲线坐标系中，可将节点扰动引力描述为关于局部坐标的函数，

u(ξ,η,ζ):R³→R(27)

在Ω_e上构造u的一个近似函数U:Ω_e→R，满足U(x_i)＝u_i(i＝1,2,…,n)，在子域Ω_e上取三元多项式类

$\begin{array}{c}\left\{g_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\right\}=\{1,\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ},{\mathrm{\ξ}}^2,{\mathrm{\η}}^2,{\mathrm{\ζ}}^2,\mathrm{\ξ}\mathrm{\η},\mathrm{\ξ}\mathrm{\ζ},\mathrm{\μ}\mathrm{\ζ},{\mathrm{\ξ}}^3,{\mathrm{\η}}^3,{\mathrm{\ζ}}^3,\\ {\mathrm{\ξ}}^2\mathrm{\η},{\mathrm{\ξ}}^2\mathrm{\ζ},{\mathrm{\ξ\η}}^2,{\mathrm{\ξ\ζ}}^2,{\mathrm{\η\ζ}}^2,{\mathrm{\η}}^2\mathrm{\ζ},\mathrm{\ξ}\mathrm{\η}\mathrm{\ζ},\mathrm{...}\}\end{array}---\left(28\right)$

的前t项为插值基函数，且r＜t＜s，即令

$\begin{array}{cc}U(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\underset{j=1}{\overset{t}{\Σ}}{a}_j{g}_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})& (\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\∈{\mathrm{\Ω}}_e\end{array}---\left(29\right)$

令

$\left\{\begin{array}{c}G={\left\{g_j({\mathrm{\ξ}}_i,{\mathrm{\η}}_i,{\mathrm{\ζ}}_i)\right\}}_{ij},i=r+1,r+2,...,s;j=1,2,...t\\ G_I={\left\{g_j({\mathrm{\ξ}}_i,{\mathrm{\η}}_i,{\mathrm{\ζ}}_i)\right\}}_{ij},i=1,2,...,r;j=1,2,...t\\ u={\[u_{r+1},u_{r+2},...,u_s\]}^T\\ u_I={\[u_1,u_2,...,u_r\]}^T\\ a={\[a_1,a_2,...,a_t\]}^T\end{array}\right.---\left(30\right)$

由式(31)求解待定系数a，

$\left\{\begin{array}{c}\min I\left(a\right)={(Ga-u)}^T(Ga-u)\\ \begin{array}{cc}s.t.& G_{I}a-u_I=0\end{array}\end{array}\right.---\left(31\right)$

将求解得到的a代入式(29)，即可根据节点扰动引力值u_i求解延拓域内任意一点的扰动引力值。

本发明的有益效果在于：

具有逼近精度高、计算速度快、数据存储量小的特征，能够满足在飞行过程中沿任意弹道快速计算扰动引力的要求，与现有研究基础相比，本发明提出的方法具有以下优点：

1)首次提出一种针对临近空间大范围机动弹道的扰动引力实时计算方法。该方法具有逼近精度高、计算速度快、数据存储量小的特征，可适应沿任意临近空间大范围机动弹道的扰动引力快速计算，且具备快速响应飞行任务临时变更的能力，将该方法应用于弹道计算中的扰动引力计算，单条弹道计算时间仅为1080阶球谐函数方法的6‰；

2)首次提出一种在发射准备阶段建立覆盖弹道三维包络的全域扰动引力重构模型策略，将全域重构模型数据作为“知识库”离线存储在地面设备，既减轻了射前模型构建的负担，又可以适应实际应用中可能存在的射前临时变更飞行任务的情况；

3)首次提出一种基于射前弹道规划和离线数据库快速构建沿参考弹道的扰动引力局域重构模型的策略，将射前模型构建简化为有效数据判读，避免多次调用高精度扰动引力计算方法求解节点扰动引力，既缩短了射前模型构建时间，又降低了弹载计算机存储量；

4)采用多层网格构建沿弹道的局域重构模型，能够有效避免弹道因大范围侧向机动超出逼近区间而导致逼近计算的发散，建立延拓逼近方法时，通过引入局部曲面坐标系使逼近函数的推导大为简化；

5)首次将延拓逼近方法应用到沿弹道的扰动引力计算中，延拓逼近使用了外节点信息，充分考虑了网格之间的联系，由于没有增加全域的节点数，计算规模没有增加，该方法集合了拟合方法和插值方法的优势，在计算量相同的情况下具有更高的赋值精度，在相同的逼近精度要求下可减小计算规模；

6)在延拓逼近方法的系数求解中，可以利用矩阵分块技术降低求解问题阶数，由于主要计算量集中到型函数矩阵的求解上，而该矩阵的值只与当前点所在延拓域有关，只有当起算点超出子域后才需要进行更新，因此具有计算量小的特征。

附图说明

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

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

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

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

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

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

图7为侧向无禁飞区全域网格与局域网格示意图；

图8为纵向无禁飞区全域网格与局域网格示意图；

图9为侧向有禁飞区全域网格与局域网格示意图；

图10为纵向有禁飞区全域网格与局域网格示意图；

图11为三维子域及延拓域示意图；

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

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

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

图15为沿不同机动弹道扰动引力逼近结果示意图；

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

具体实施方式

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

基于CAV-H飞行器模型进行仿真，基本仿真条件设置为：①质量m＝908kg，参考面积Sm＝0.48375m²；②再入点状态参数：速度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开发。

其具体步骤如下：

A1：建立弹道包络侧向宽度与纵程的数学模型

令弹道再入点为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°。记各落点对应的弹道纵程为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)

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

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

对L_i(i＝1,2,3,4,…)和W_i(i＝1,2,3,4,…)进行拟合，建立如式(3)所示的弹道包络侧向宽度与纵程的数学模型，根据得出

W＝-1.064813L+0.000879L²-1.225047e^-7L³+4.607571e^-12L⁴(3)

临近空间大范围机动弹道三维包络示意图见图3。

A2：建立换极坐标系

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

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

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

A3：建立换极坐标系中飞行器动力学模型

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

$\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 \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)$

其中，

$\left\{\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}\right.---\left(6\right)$

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

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

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

定义

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

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

$\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.\end{array}---\left(11\right)$

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

$\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)$

其位置关系图见图6，

其中，

$\left\{\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}\right.---\left(14\right)$

A4：计算换极坐标系弹道三维包络

根据步骤A3中的坐标转换关系，针对再入点经度λ_I＝0°、再入点地心纬度φ_I＝0°、落点经度λ_T＝30°、落点地心纬度φ_T＝40°的临近空间大范围机动弹道进行计算。根据第三步所述坐标转换关系，可计算得到换极坐标系弹道参数

由式(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由第一步确定，

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

$\left\{\begin{array}{c}{\hat{E}}_{\mathrm{\λ}1}={\mathrm{\λ}}_I\\ {\hat{E}}_{\mathrm{\λ}2}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\end{array}\right.,\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.---\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)$

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

A5：换极坐标系全域空域剖分

记扰动引力全域重构模型为Field模型，其确定的空域为Ω，为保证首个和末个延拓点有意义，按式(19)计算换极系Ω边界，

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

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

表5换极系全域空域范围

按照天向dr、东向dλ、北向dφ的间隔将由式(19)确定的换极系全域空域Ω均匀剖分为q个互不重叠的子域Ω_e(e＝1,2,…,q)，记网格节点坐标为 $\hat{N}({\hat{r}}_G,{\widehat{\mathrm{\λ}}}_G,{\widehat{\mathrm{\φ}}}_G),$

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

A6：建立一般坐标系扰动引力全域重构模型

根据换极系网格节点坐标由第三步所述坐标转换关系计算一般坐标系网格节点坐标N(r_G,λ_G,φ_G)，

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

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

其中，

${\stackrel{\‾}{C}}_{nm}^{*}=\left\{\begin{array}{cc}0& \begin{array}{cc}n=2& m=0\end{array}\\ {\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_G}=-\frac{\mathrm{\μ}}{{r_G}^2}{\underset{n=2}{\overset{s}{\Σ}}(n+1)\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_G+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G{\mathrm{cos\φ}}_G}\frac{\∂T}{\∂{\mathrm{\λ}}_G}=-\frac{1}{{\mathrm{cos\φ}}_G}\frac{\mathrm{\μ}}{{r_G}^2}{\underset{n=2}{\overset{s}{\Σ}}\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}m({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{sinm\λ}}_G+{\stackrel{\‾}{S}}_{nm}{\mathrm{cosm\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G}\frac{\∂T}{\∂{\mathrm{\φ}}_G}=\frac{\mathrm{\μ}}{{r_G}^2}{\underset{n=2}{\overset{s}{\Σ}}\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}{\mathrm{cosm\λ}}_G+{\stackrel{\‾}{S}}_{nm}{\mathrm{sinm\λ}}_G)\frac{d{\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)}{{\mathrm{d\φ}}_G}\end{array}\right.---\left(24\right)$

存储一般坐标系节点位置三分量和节点扰动引力三分量，完成扰动引力全域重构模型(Field模型)构建；

A7：建立沿弹道的扰动引力局域重构模型

记沿弹道的扰动引力局域重构模型为Channel模型，根据飞行任务进行弹道规划，计算参考弹道，判断参考弹道历经的三层延拓网格，从Field模型数据中读取历经延拓网格的节点位置三分量及节点扰动引力三分量，完成Channel模型构建；

针对第四步所述的不考虑禁飞区的临近空间大范围机动弹道，其侧向无禁飞区全域网格与局域网格示意图见图7，纵向无禁飞区全域网格与局域网格示意图见图8。由图7和图8可见，局域网格包含在全局网格之中，且能覆盖可行弹道，说明本方法建立的全域重构模型和局域重构模型可行。

针对再入点经度λ_I＝100°、再入点地心纬度φ_I＝20°；落点经度λ_T＝165°、落点地心纬度φ_T＝50°；禁飞区中心经度λ_nf＝138°，中心地心纬度φ_nf＝40°，高度H_nf＝100km，半径R_nf＝1000km的临近空间大范围机动弹道进行计算。计算得到侧向有禁飞区全域网格与局域网格示意图见图9，纵向有禁飞区全域网格与局域网格示意图见图10。由图9和图10可见，可行弹道向左侧绕飞规避了禁飞区，且能被局域网格覆盖，说明本方法建立的全域重构模型和局域重构模型适用于考虑禁飞区的情况。

A8：建立子域局部曲线坐标系

令子域Ω_e及其相邻网格确定的空域为延拓域Ω′_e，有记Ω_e包含的节点数为r，Ω′_e包含的节点数为s，有s＞r，三维子域及延拓域示意图见图11。

建立子域局部曲线坐标系，在局部坐标系中描述单元节点和计算点坐标，子域Ω_e的局部曲线坐标系由半径r_l＝r_G+dr/2的球面、经度λ_l＝λ_G+dλ/2的子午面、纬度φ_l＝φ_G+dφ/2的纬圈的交线组成，原点l(r_l,λ_l,φ_l)为三交线的交点，原点局部坐标为l(0,0,0)，ξ、η、ζ分别沿原点l的天向、北向和东向，则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为，

$\left\{\begin{array}{c}\mathrm{\ξ}=r-r_l\\ \mathrm{\η}=r\cos \mathrm{\φ}(\mathrm{\λ}-{\mathrm{\λ}}_l)\\ \mathrm{\ζ}=r(\mathrm{\φ}-{\mathrm{\φ}}_l)\end{array}\right.---\left(25\right)$

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

$\left\{\begin{array}{c}{\mathrm{\ξ}}_i=r_i-r_l\\ {\mathrm{\η}}_i=r_i{\mathrm{cos\φ}}_i({\mathrm{\λ}}_i-{\mathrm{\λ}}_l)\\ {\mathrm{\ζ}}_i=r_i({\mathrm{\φ}}_i-{\mathrm{\φ}}_l)\end{array}\right.---\left(26\right)$

A9：飞行过程中扰动引力延拓逼近计算

令延拓域Ω′_e节点扰动引力值为u_i，在子域局部曲线坐标系中，可将节点扰动引力描述为关于局部坐标的函数，

u(ξ,η,ζ):R³→R(27)

在Ω_e上构造u的一个近似函数U:Ω_e→R，满足U(x_i)＝u_i(i＝1,2,…,n)，在子域Ω_e上取三元多项式类

$\begin{array}{c}\left\{g_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\right\}=\{1,\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ},{\mathrm{\ξ}}^2,{\mathrm{\η}}^2,{\mathrm{\ζ}}^2,\mathrm{\ξ}\mathrm{\η},\mathrm{\ξ}\mathrm{\ζ},\mathrm{\μ}\mathrm{\ζ},{\mathrm{\ξ}}^3,{\mathrm{\η}}^3,{\mathrm{\ζ}}^3,\\ {\mathrm{\ξ}}^2\mathrm{\η},{\mathrm{\ξ}}^2\mathrm{\ζ},{\mathrm{\ξ\η}}^2,{\mathrm{\ξ\ζ}}^2,{\mathrm{\η\ζ}}^2,{\mathrm{\η}}^2\mathrm{\ζ},\mathrm{\ξ}\mathrm{\η}\mathrm{\ζ},\mathrm{...}\}\end{array}---\left(28\right)$

的前t项为插值基函数，且r＜t＜s，即令

$\begin{array}{cc}U(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\underset{j=1}{\overset{t}{\Σ}}{a}_j{g}_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})& (\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\∈{\mathrm{\Ω}}_e\end{array}---\left(29\right)$

取r＝8、s＝32、t＝20。令子域Ω_e节点编号依次为0～7，延拓域Ω′_e中除Ω_e节点外的其他节点编号为8～31，见图11。则有，

$\left\{\begin{array}{c}G={\left\{g_j(x_i,y_i,z_i)\right\}}_{ij},i=9,10,...,32;j=1,2,...t\\ G_I={\left\{g_j(x_i,y_i,z_i)\right\}}_{ij},i=1,2,...,8;j=1,2,...t\\ u={\[u_9,u_{10},...,u_{32}\]}^T\\ u_I={\[u_1,u_2,...,u_8\]}^T\\ a={\[a_1,a_2,...,a_{20}\]}^T\end{array}\right.---\left(30\right)$

其中，待定系数a可由式(31)求解，

$\left\{\begin{array}{c}\min I\left(a\right)={(Ga-u)}^T(Ga-u)\\ \begin{array}{cc}s.t.& G_{I}a-u_I=0\end{array}\end{array}\right.---\left(31\right)$

引入拉格朗日乘子l＝[l₁,l₂,…,l₈]^T，构造式(32)所示的拉格朗日函数，则式(32)的最小值minL即为式(31)所描述问题的解。

L(a,l)＝(Ga-u)^T(Ga-u)+2(G_Ia-u_I)l(32)

根据优化原理，minL可由式(33)确定，

$\left\{\begin{array}{c}\frac{\∂L}{\∂a}=2{(Ga-u)}^{T}G+2l^T{G}_I=0\\ \frac{\∂L}{\∂l}=2{(G_{I}a-u_I)}^T=0\end{array}\right.---\left(33\right)$

写成矩阵形式，即

$\left[\begin{array}{cc}A& C^T\\ C& 0\end{array}\right]\left[\begin{array}{c}a\\ l\end{array}\right]=\left[\begin{array}{c}{F}_0\\ F_1\end{array}\right]---\left(34\right)$

其中，

$\left\{\begin{array}{c}A=G^{T}G\\ C=G_1\end{array}\right.,\left\{\begin{array}{c}{F}_0=G^{T}u\\ F_1=u_I\end{array}\right.---\left(35\right)$

为提高式(34)的求解效率，引入型函数。注意到A为方阵且可逆，令

$\left\{\begin{array}{c}B={\left(G^{T}G\right)}^{-1}\\ D_{22}={(-G_I{\mathrm{BG}}_I^T)}^{-1}\\ D_{12}=-{\mathrm{BG}}_I^T{D}_{22}\\ D_{11}=B-D_{12}{G}_I^{T}B\end{array}\right.---\left(36\right)$

则可将28阶矩阵求逆问题转化为20阶和8阶矩阵的求逆问题，减少了计算量。求解式(34)可得，

$a=\[\begin{array}{cc}{D}_{11}& D_{12}\end{array}\]\left[\begin{array}{c}{F}_0\\ F_1\end{array}\right]=\[\begin{array}{cc}{D}_{11}{G}^T& D_{12}\end{array}\]\left[\begin{array}{c}u\\ u_I\end{array}\right]---\left(37\right)$

主要的计算量集中在矩阵[D₁₁G^TD₁₂]的求解上，而该矩阵的值只与当前点所在子域及其延拓域有关，因此只有当起算点超出子域后才需要进行更新。将求解得到的a代入式(29)，即可根据节点扰动引力值u_i求解延拓域内任意一点的扰动引力值。

取节点扰动引力采用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°。沿不同区域弹道扰动引力逼近结果示意图见图12，逼近误差见表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°。沿不同射向弹道扰动引力逼近结果示意图见图13，逼近误差见表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°。沿不同射程弹道扰动引力逼近结果示意图见图14，逼近误差见表8。

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

针对起点为λ_I＝0°，φ_I＝0°；终点为λ_T＝50°，φ_T＝0°的射程约5560km的东向发射弹道进行计算。分别考虑：①右侧最大范围机动弹道，初始方位角为σ_max＝126.2°；②左侧最大范围机动弹道，初始方位角为σ_min＝53.8°；③考虑禁飞区的中程机动弹道，禁飞区中心经度为λ_nf＝38°，中心地心纬度为φ_nf＝4°，高度为H_nf＝100km，半径R_nf＝1000km。沿不同机动弹道扰动引力逼近结果示意图见图15，逼近误差见表9。

表9不同机动弹道扰动引力逼近精度统计

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

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

不同网格划分对应的弹上存储量及弹道终端状态偏差见表10，全域模型存储量及重构时间见表11，局域模型存储量及重构时间见表12，采用不同扰动引力计算方法所需的弹道计算时间见表13。计算表明，在10km/20°/20°的网格划分下，弹上数据存储量为816，射前模型重构时间约为0.768s，弹道计算时间约为76.781s，由逼近误差引起的弹道终端位置偏差小于100m。相较于1080阶球谐函数法，采用延拓逼近方法进行弹道计算中的扰动引力计算，单条弹道计算时间仅为前者的6‰。

表10不同网格划分对应的弹上存储量及弹道终端状态偏差

表11全域模型存储量及重构时间

表12局域模型存储量及数据获取时间

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

计算方法	计算时间(s)
		延拓逼近方法(10km/1.2°/1.2°)	82.484
延拓逼近方法(10km/20°/20°)	76.781
		50阶球谐函数方法	281.094
1080阶球谐函数方法	125970.109
		Stokes方法	104653.258

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

1)由本发明确定的方法进行沿临近空间大范围机动弹道的扰动引力计算，对位于不同区域，具有不同射程、不同射向和不同机动模式的弹道均有较好的适应性。

2)相较于非延拓逼近方法，通过本发明计算得到的扰动引力具有更高逼近精度，在10km/1.2°/1.2°和10km/20°/20°的网格划分下，可使由逼近误差导致的滑翔弹道终端位置偏差分别小于1m和100m，使弹道精度分别提高约1e4个量级和1e2个量级。

3)方法具有存储量小、计算速度快的特征。在10km/20°/20°的网格划分下，弹上数据存储量仅为816，射前模型重构时间约为0.768s，弹道计算时间约为76.781s。相较于1080阶球谐函数法，采用延拓逼近方法进行弹道计算中的扰动引力计算，单条弹道计算时间仅为前者的6‰。

4)提出的方法具有逼近精度高、计算速度快、适应范围广的特征，且具备快速响应飞行任务临时变更的能力。

以上仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围内。

Claims (1)

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

A1：建立弹道包络侧向宽度与纵程的数学模型

令弹道再入点为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)$

A2：建立换极坐标系

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

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

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

A3：建立换极坐标系中飞行器动力学模型

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

$\left\{\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \widehat{\mathrm{\θ}}\\ \frac{d\hat{r}}{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}}\left(\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\φ}}-\cos \widehat{\mathrm{\σ}}\sin \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\φ}}\right)+C_{\mathrm{\θ}}+{\tilde{C}}_{\mathrm{\θ}}\\ \frac{d\tilde{V}}{dt}=-D+{\hat{g}}_r\sin \widehat{\mathrm{\θ}}+{\hat{g}}_{{\mathrm{\ω}}_e}\left(\cos \widehat{\mathrm{\σ}}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\φ}}+\sin \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\φ}}\right)+{\tilde{C}}_V\end{array}\right.---\left(4\right)$

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

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

其中，

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

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

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

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

定义

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

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

$\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.\end{array}---\left(11\right)$

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

$\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)$

其中，

$\left\{\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}\right.---\left(14\right)$

A4：计算换极坐标系弹道三维包络

根据第三步所述坐标转换关系，可根据弹道再入点经度λ_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由第一步确定，

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

$\left\{\begin{array}{c}{\hat{E}}_{\mathrm{\λ}1}={\mathrm{\λ}}_I\\ {\hat{E}}_{\mathrm{\λ}2}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\end{array}\right.,\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.---\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倍作为弹道包络的垂向范围

A5：换极坐标系全域空域剖分

记扰动引力全域重构模型为Field模型，其确定的空域为Ω，为保证首个和末个延拓点有意义，按式(19)计算换极系Ω边界，

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

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

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

A6：建立一般坐标系扰动引力全域重构模型

根据换极系网格节点坐标由第三步所述坐标转换关系计算一般坐标系网格节点坐标N(r_G,λ_G,φ_G)，

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

$T_G=\frac{\mathrm{\μ}}{r_G}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\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_G}=-\frac{\mathrm{\μ}}{{r_G}^2}\underset{n=2}{\overset{s}{\mathrm{\Σ}}}\left(n+1\right){\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}\left({\stackrel{\‾}{C}}_{nm}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G\right){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G{\mathrm{cos\φ}}_G}\frac{\∂T}{\∂r_G}=-\frac{1}{{\mathrm{cos\φ}}_G}\frac{\mathrm{\μ}}{{r_G}^2}\underset{n=2}{\overset{s}{\mathrm{\Σ}}}{\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}m\left({\stackrel{\‾}{C}}_{nm}^{*}\sin {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\cos {\mathrm{m\λ}}_G\right){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G}\frac{\∂T}{\∂{\mathrm{\φ}}_G}=\frac{\mathrm{\μ}}{{r_G}^2}\underset{n=2}{\overset{s}{\mathrm{\Σ}}}{\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}\left({\stackrel{\‾}{C}}_{nm}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G\right)\frac{d{\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)}{{\mathrm{d\φ}}_G}\end{array}\right.---\left(24\right)$

存储一般坐标系节点位置三分量和节点扰动引力三分量，完成扰动引力全域重构模型(Field模型)构建；

A7：建立沿弹道的扰动引力局域重构模型

记沿弹道的扰动引力局域重构模型为Channel模型，根据飞行任务进行弹道规划，计算参考弹道，判断参考弹道历经的三层延拓网格，从Field模型数据中读取历经延拓网格的节点位置三分量及节点扰动引力三分量，完成扰动引力局域重构模型(Channel模型)构建；

A8：建立子域局部曲线坐标系

令子域Ω_e及其相邻网格确定的空域为延拓域Ω′_e，有记Ω_e包含的节点数为r，Ω′_e包含的节点数为s，有s＞r，

建立子域局部曲线坐标系，在局部坐标系中描述单元节点和计算点坐标，子域Ω_e的局部曲线坐标系由半径r_l＝r_G+dr/2的球面、经度λ_l＝λ_G+dλ/2的子午面、纬度φ_l＝φ_G+dφ/2的纬圈的交线组成，原点l(r_l,λ_l,φ_l)为三交线的交点，原点局部坐标为l(0,0,0)，ξ、η、ζ分别沿原点l的天向、北向和东向，则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为，

$\left\{\begin{array}{c}\mathrm{\ξ}=r-r_l\\ \mathrm{\η}=r\cos \mathrm{\φ}(\mathrm{\λ}-{\mathrm{\λ}}_l)\\ \mathrm{\ζ}=r(\mathrm{\φ}-{\mathrm{\φ}}_l)\end{array}\right.---\left(25\right)$

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

$\left\{\begin{array}{c}{\mathrm{\ξ}}_i=r_i-r_l\\ {\mathrm{\η}}_i=r_i{\mathrm{cos\φ}}_i({\mathrm{\λ}}_i-{\mathrm{\λ}}_l)\\ {\mathrm{\ζ}}_i=r_i({\mathrm{\φ}}_i-{\mathrm{\φ}}_l)\end{array}\right.---\left(26\right)$

A9：飞行过程中扰动引力延拓逼近计算

令延拓域Ω′_e节点扰动引力值为u_i，在子域局部曲线坐标系中，可将节点扰动引力描述为关于局部坐标的函数，

u(ξ,η,ζ):R³→R(27)

在Ω_e上构造u的一个近似函数U:Ω_e→R，满足U(x_i)＝u_i(i＝1,2,…,n)，在子域Ω_e上取三元多项式类

$\begin{array}{c}\left\{g_j\left(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ}\right)\right\}=\{1,\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ},{\mathrm{\ξ}}^2,{\mathrm{\η}}^2,{\mathrm{\ζ}}^2,\mathrm{\ξ}\mathrm{\η},\mathrm{\ξ}\mathrm{\ζ},\mathrm{\μ}\mathrm{\ζ},{\mathrm{\ξ}}^3,{\mathrm{\η}}^3,{\mathrm{\ζ}}^3,\\ {\mathrm{\ξ}}^2\mathrm{\η},{\mathrm{\ξ}}^2\mathrm{\ζ},{\mathrm{\ξ\η}}^2,{\mathrm{\ξ\ζ}}^2,{\mathrm{\η\ζ}}^2,{\mathrm{\η}}^2\mathrm{\ζ},\mathrm{\ξ}\mathrm{\η}\mathrm{\ζ},\mathrm{...}\}\end{array}---\left(28\right)$

的前t项为插值基函数，且r＜t＜s，即令

$U(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\underset{j=1}{\overset{t}{\Σ}}{a}_j{g}_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ}),(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\∈{\mathrm{\Ω}}_e---\left(29\right)$

令

$\left\{\begin{array}{c}G={\left\{g_j\left({\mathrm{\ξ}}_i,{\mathrm{\η}}_i,{\mathrm{\ζ}}_i\right)\right\}}_{ij},i=r+1,r+2,\mathrm{...},s;j=1,2,\mathrm{...}t\\ G_I={\left\{g_j\left({\mathrm{\ξ}}_i,{\mathrm{\η}}_i,{\mathrm{\ζ}}_i\right)\right\}}_{ij},i=1,2,\mathrm{...},r;j=1,2,\mathrm{...}t\\ u={\[u_{r+1},u_{r+2},\mathrm{...},u_s\]}^T\\ u_I={\[u_1,u_2,\mathrm{...}{u}_r\]}^T\\ a={\[a_1,a_2,\mathrm{...},a_t\]}^T\end{array}\right.---\left(30\right)$

由式(31)求解待定系数a，

$\left\{\begin{array}{cc}min& I\left(a\right)={(Ga-u)}^T(Ga-u)\\ s.t.& G_{I}a-u_I=0\end{array}\right.---\left(31\right)$

将求解得到的a代入式(29)，即可根据节点扰动引力值u_i求解延拓域内任意一点的扰动引力值。