沿滑翔弹道的扰动引力重构模型优化方法.pdf

本发明公开了一种沿滑翔弹道的扰动引力重构模型优化方法，包括以下步骤：首先，计算换极坐标系弹道三维包络；其次，进行换极坐标系全域空域剖分；再次，建立一般坐标系扰动引力全域重构模型；然后，建立沿弹道的扰动引力局域重构模型；基于此，建立飞行过程中扰动引力快速逼近算法；之后，开展试验设计并构建重构模型的代理模型；最后，针对代理模型建立优化方法。根据本发明方法进行沿滑翔弹道的扰动引力重构模型优化，可以获得给定精度要求下的存储量最小模型和给定存储量要求下的精度最高模型，能够满足弹上实时计算对数据存储量、计算速度和重构精度的要求，具有较好的工程应用前景。

1.一种沿滑翔弹道的扰动引力重构模型优化方法，其特征在于：包括以下步骤：
第一步，建立弹道包络侧向宽度与纵程的数学模型；
第二步，建立换极坐标系；
第三步，建立换极坐标系中飞行器动力学模型；
第四步，计算换极坐标系弹道三维包络；
第五步，换极坐标系全域空域剖分；
第六步，建立一般坐标系扰动引力全域重构模型；
第七步，建立沿弹道的扰动引力局域重构模型；
第八步，建立飞行过程中扰动引力实时逼近算法；
第九步，获取代理模型训练样本计算条件；
第十步，建立重构模型的代理模型；
第十一步，建立重构模型的优化算法。
2.根据权利要求1所述的沿滑翔弹道的扰动引力重构模型优化方法，其特征在于：所述
滑翔弹道特指高超声速滑翔飞行器全程弹道的一部分，弹道起点为再入点，弹道终点为中
末交班点。
3.根据权利要求2所述的沿滑翔弹道的扰动引力重构模型优化方法，其特征在于：
第一步，建立弹道包络侧向宽度与纵程的数学模型
令弹道再入点为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_e arccos(sinφ_I sinφ_T+cosφ_I cosφ_T cos(λ_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)$
第二步，建立换极坐标系
引入一个换极坐标系，为表述方便，用表示换极坐标系中各物理量，用X表示一般坐标
系中各物理量，按如下方式建立换极坐标系：
①定义一个再入大圆弧平面作为换极赤道平面：1)对目标点确定的情况，将再入点和
目标点地心矢径构成的再入大圆弧平面作为换极赤道平面；2)对于目标点未确定的情况，
根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面；
②基于换极赤道平面定义换极坐标系O_E为地心，轴沿再入点地心矢径方
向，轴在换极赤道面内垂直于轴指向目标点方向，轴与轴、轴构成右手系；
第三步，建立换极坐标系中飞行器动力学模型
在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程，其飞行状态量为换极
后的经度地心纬度航迹偏航角速度速度倾角和地心距
$\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)$
其中，
$\left\{\begin{array}{c}{\mathrm{\ω}}_{ex}={\mathrm{\ω}}_e(\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)\\ {\mathrm{\ω}}_{ey}={\mathrm{\ω}}_e(-\sin \widehat{\mathrm{\λ}}{\mathrm{cos\φ}}_p\cos A_p+\cos \widehat{\mathrm{\λ}}{\mathrm{cos\φ}}_p\sin A_p)\\ {\mathrm{\ω}}_{ez}={\mathrm{\ω}}_e(-\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)\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)$
定义
$\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的方位角；
$\left\{\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}\right.---\left(9\right)$
$\left\{\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}\right.---\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.\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)$
第四步，计算换极坐标系弹道三维包络
根据第三步所述坐标转换关系，可根据弹道再入点经度λ_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倍作为弹道包络的垂向范围
第五步，换极坐标系全域空域剖分
记扰动引力全域重构模型为Field模型，其确定的空域为Ω；为保证首个和末个延拓点
有意义，按式(19)计算换极系Ω边界，
$\left\{\begin{array}{c}{\hat{F}}_{r1}={\hat{r}}_I-\mathrm{\Δ}\hat{r}-2\·dr\\ {\hat{F}}_{r2}={\hat{r}}_I+2\·dr\end{array}\right.,\left\{\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}\right.,\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)$
其中，为换极系Ω天向下边界，为换极系Ω天向上边界；为换极系Ω东向下边
界，为换极系Ω东向上边界；为换极系Ω北向下边界，为换极系Ω北向上边界；dr、
dλ和dφ分别为空域剖分间隔；
按照天向dr、东向dλ、北向dφ的间隔将由式(19)确定的换极系全域空域Ω均匀剖分为
q个互不重叠的子域Ω_e(e＝1,2,…,q)，记网格节点坐标为
$\mathrm{\Ω}=\underset{e=1}{\overset{q}{\∪}}{\mathrm{\Ω}}_e---\left(20\right)$
第六步，建立一般坐标系扰动引力全域重构模型
根据换极系网格节点坐标由第三步所述坐标转换关系计算一般坐标系网
格节点坐标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}{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}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_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}^{*}\sin {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\cos {\mathrm{m\λ}}_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}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G)\frac{d{\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)}{{\mathrm{d\φ}}_G}\end{array}\right.---\left(24\right)$
存储一般坐标系节点位置三分量和节点扰动引力三分量，完成扰动引力全域重构模型
(Field模型)构建；
第七步，建立沿弹道的扰动引力局域重构模型
记沿弹道的扰动引力局域重构模型为Channel模型；
根据飞行任务进行弹道规划，计算参考弹道；
判断参考弹道历经的三层延拓网格，从Field模型数据中读取历经延拓网格的节点位
置三分量及节点扰动引力三分量，完成扰动引力局域重构模型(Channel模型)构建；
第八步，建立飞行过程中扰动引力实时逼近算法
令子域Ω_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{\η}=rcos\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)$
令延拓域Ω′_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}{ccc}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}minI\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)第五步所述的重构模型网格划分间隔dr、dλ、dφ；
(2)第八步所述的延拓域节点数s；
(3)第八步所述的三元多项式类项数t；
待设计因素的设计范围为设定值，选取依据为重构精度要求；
经实验设计，可获得N_experi个设计因素组合Xⁱ＝(dr,dλ,dφ,t,s)，其中i＝1,2,…,
N_experi，N_experi为设定值；
第十步，建立重构模型的代理模型
根据第九步计算得到的不同水平下的设计因素组合Xⁱ＝(dr,dλ,dφ,t,s)开展考虑扰
动引力的滑翔弹道仿真；弹道仿真条件由飞行任务确定；扰动引力根据第一步至第八步所
述的扰动引力重构方法计算得到；
经弹道仿真，可以获取不同设计因素组合Xⁱ＝(dr,dλ,dφ,t,s)对应的Channel模型节
点数N_ch-node以及由扰动引力引起的滑翔弹道终端位置偏差dD；
以Xⁱ＝(dr,dλ,dφ,t,s)作为训练代理模型的输入对，以Yⁱ＝(N,dD)作为训练代理的输
出对，通过遗传算法或其他拟合方法建立描述输入对与输出对之间函数关系的代理模型；
第十一步，建立重构模型的优化算法
考虑如下两类优化问题：
(1)优化问题一：给定弹道精度要求，建立存储量最小模型；
目标函数：minN；
待设计变量：dr、dλ、dφ、t、s；
约束条件：dD∈[dD_min,dD_max]；
待设计变量取值范围：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,dφ_max]、t∈
[t_min,t_max]、s∈[s_min,s_max]；
(2)优化问题二：给定存储量要求，建立重构精度最高模型；
目标函数：mindD；
待设计变量：dr、dλ、dφ、t、s；
约束条件：N∈(0,N_max]；
待设计变量取值范围：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,dφ_max]、t∈
[t_min,t_max]、s∈[s_min,s_max]；
选取一种全局优化方法求解上述优化问题，最终可建立给定精度要求下的存储量最小
模型和给定存储量要求下的精度最高模型。

沿滑翔弹道的扰动引力重构模型优化方法

技术领域

本发明属于飞行器动力学建模领域，特别涉及沿滑翔弹道的扰动引力重构模型优
化方法。

背景技术

本发明涉及的滑翔弹道特指高超声速滑翔飞行器包含初始下降段弹道和滑翔段
弹道在内的弹道，其起点为再入点，终点为滑翔段中末交班点。滑翔弹道长时间处于临近空
间，受地球扰动引力影响显著。计算表明，由扰动引力引起的滑翔弹道终端位置偏差可达几
十公里，因此有必要在弹道计算和制导解算中考虑扰动引力的影响。现有扰动引力计算方
法通常具有计算规模大、数据存储量大、计算时间长的特点，无法满足弹上实时计算对轻量
化和计算效率的要求。为寻求最优重构模型参数，提出一种沿滑翔弹道的扰动引力重构模
型优化方法。方法首先计算换极坐标系弹道三维包络；其次进行换极坐标系全域空域剖分；
再次建立一般坐标系扰动引力全域重构模型；然后建立沿弹道的扰动引力局域重构模型；
基于此建立飞行过程中扰动引力快速逼近算法；之后开展试验设计并构建重构模型的代理
模型；最后针对代理模型建立优化方法。方法可获得给定精度要求下的存储量最小模型以
及给定存储量要求下的精度最高模型，具有较好的工程应用前景。

发明内容

本发明针对沿滑翔弹道的扰动引力重构模型优化问题，提出一种基于代理模型的
重构模型的沿滑翔弹道的扰动引力重构模型优化方法，可获得给定精度要求下的存储量最
小模型以及给定存储量要求下的精度最高模型，且具有精度高、速度快、存储量小的特征。

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

一种沿滑翔弹道的扰动引力重构模型优化方法，包括以下步骤：

第一步，建立弹道包络侧向宽度与纵程的数学模型；

第二步，建立换极坐标系；

第三步，建立换极坐标系中飞行器动力学模型；

第四步，计算换极坐标系弹道三维包络；

第五步，换极坐标系全域空域剖分；

第六步，建立一般坐标系扰动引力全域重构模型；

第七步，建立沿弹道的扰动引力局域重构模型；

第八步，建立飞行过程中扰动引力实时逼近算法；

第九步，获取代理模型训练样本计算条件；

第十步，建立重构模型的代理模型；

第十一步，建立重构模型的优化算法。

具体地，所述滑翔弹道特指高超声速滑翔飞行器全程弹道的一部分，弹道起点为
再入点，弹道终点为中末交班点。

对上述步骤做进一步说明：

第一步，建立弹道包络侧向宽度与纵程的数学模型

令弹道再入点为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，对和W_i(i＝1,
2,3,4,…)进行拟合，可建立如式(3)所示的弹道包络侧向宽度与纵程的数学模型，确定模
型系数a_i，

第二步，建立换极坐标系

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

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

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

第三步，建立换极坐标系中飞行器动力学模型

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

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

其中，

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

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

定义

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

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

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

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

其中，

第四步，计算换极坐标系弹道三维包络

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

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

其中，R_e为地球半径；

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

其中，模型系数a_i由第一步确定；

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

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

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

考虑机动和偏差，以弹道高度范围的C_r倍作为弹道包络的垂向范围

第五步，换极坐标系全域空域剖分

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

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

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

第六步，建立一般坐标系扰动引力全域重构模型

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

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

其中，

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

δg＝gradT_G (23)

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

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

第七步，建立沿弹道的扰动引力局域重构模型

记沿弹道的扰动引力局域重构模型为Channel模型；

根据飞行任务进行弹道规划，计算参考弹道；

判断参考弹道历经的三层延拓网格，从Field模型数据中读取历经延拓网格的节
点位置三分量及节点扰动引力三分量，完成扰动引力局域重构模型(Channel模型)构建；

第八步，建立飞行过程中扰动引力实时逼近算法

令子域Ω_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(ξ,η,ζ)为，

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

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

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

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

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

令

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

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

第九步，获取代理模型训练样本计算条件

基于全因子设计、正交设计或拉丁超立方设计开展实验设计，获取不同待设计因
素在设计范围内不同水平下的量值；

待设计因素包括：

(1)第五步所述的重构模型网格划分间隔dr、dλ、dφ；

(2)第八步所述的延拓域节点数s；

(3)第八步所述的三元多项式类项数t；

待设计因素的设计范围为设定值，选取依据为重构精度要求；

经实验设计，可获得N_experi个设计因素组合Xⁱ＝(dr,dλ,dφ,t,s)，其中i＝1,
2,…,N_experi，N_experi为设定值；

第十步，建立重构模型的代理模型

根据第九步计算得到的不同水平下的设计因素组合Xⁱ＝(dr,dλ,dφ,t,s)开展考
虑扰动引力的滑翔弹道仿真；弹道仿真条件由飞行任务确定；扰动引力根据第一步至第八
步所述的扰动引力重构方法计算得到；

经弹道仿真，可以获取不同设计因素组合Xⁱ＝(dr,dλ,dφ,t,s)对应的Channel模
型节点数N_ch-node以及由扰动引力引起的滑翔弹道终端位置偏差dD；

以Xⁱ＝(dr,dλ,dφ,t,s)作为训练代理模型的输入对，以Yⁱ＝(N,dD)作为训练代理
的输出对，通过遗传算法或其他拟合方法建立描述输入对与输出对之间函数关系的代理模
型；

第十一步，建立重构模型的优化算法

考虑如下两类优化问题：

(1)优化问题一：给定弹道精度要求，建立存储量最小模型。

目标函数：minN；

待设计变量：dr、dλ、dφ、t、s；

约束条件：dD∈[dD_min,dD_max]；

待设计变量取值范围：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,d
φ_max]、t∈[t_min,t_max]、s∈[s_min,s_max]；

(2)优化问题二：给定存储量要求，建立重构精度最高模型。

目标函数：mindD；

待设计变量：dr、dλ、dφ、t、s；

约束条件：N∈(0,N_max]；

待设计变量取值范围：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,d
φ_max]、t∈[t_min,t_max]、s∈[s_min,s_max]；

选取一种全局优化方法求解上述问题，最终可建立给定精度要求下的存储量最小
模型和给定存储量要求下的精度最高模型。

至此，经过上述十一步可最终建立一种沿滑翔弹道的扰动引力重构模型优化方
法。方法可获得给定精度要求下的存储量最小模型以及给定存储量要求下的精度最高模
型，且具有精度高、速度快、存储量小的特征。与现有研究基础相比，本发明提出的方法具有
以下优点：

1)首次提出一种针对沿滑翔弹道扰动引力重构模型的优化方法，方法能快速完成
沿任意滑翔弹道的扰动引力计算，具有逼近精度高、计算速度快、数据存储量小的特征。

2)提出一种基于代理模型开展重构模型优化设计的思路，代理模型基于试验设计
方法和高精度拟合方法构建，具有较高的拟合精度，终端位置偏差的拟合误差能控制在
0.5％以内，节点数基本无误差。

3)建立的优化方法可以获得任意给定精度要求下的存储量最小模型，当终端位置
偏差限定为5m时，最优模型对应的弹上数据存储量仅为377，可以满足弹上数据存储量的要
求。

4)建立的优化方法可以获得任意给定存储量要求下的精度最高模型，当节点数限
定400时，最优模型对应的弹道终端位置偏差仅为4.030m，可满足重构精度要求。

附图说明

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

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

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

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

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

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

图7为侧向全域网格与局域网格示意图；

图8为纵向全域网格与局域网格示意图；

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

图10为RBF神经网络结构示意图；

具体实施方式

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

基于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-3470CPU@3.20GHz，内存为3.46GB。软件
环境为Window XP操作系统，计算程序基于VC++6.0开发。

其具体步骤如下：

第一步，建立弹道包络侧向宽度与纵程的数学模型

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

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

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

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

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

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

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

其中，R_e为地球半径。

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

表3不同落点对应的纵程

通过改变再入点初始速度方位角获得由飞行器能力决定的最大左侧机动弹道和
最大右侧向机动弹道。不同纵程弹道侧向包络示意图见图2。

记各落点对应的弹道侧向包络宽度为W_io。由式(33)计算W_io，

W_io＝B_l-B_r (33)

考虑建模误差，取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,…)进行拟合，可建立如式(34)所示的弹
道包络侧向宽度与纵程的数学模型，

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

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

第二步，建立换极坐标系

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

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

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

第三步，建立换极坐标系中飞行器动力学模型

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

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

其中，

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

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

定义

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

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

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

由一般坐标系中σ确定换极坐标系中的表达式为，其位置关系图见图6，

其中，

第四步，计算换极坐标系弹道三维包络

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

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

其中，R_e为地球半径。

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

其中，模型系数a_i由第一步确定。

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

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

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

由计算可得，考虑机动和偏差，取

第五步，换极坐标系全域空域剖分

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

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

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

表5换极系全域空域范围

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

第六步，建立一般坐标系扰动引力全域重构模型

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

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

其中，

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

δg＝gradT_G (54)

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

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

第七步，建立沿弹道的扰动引力局域重构模型

记沿弹道的扰动引力局域重构模型为Channel模型。

根据飞行任务进行弹道规划，计算参考弹道。

判断参考弹道历经的三层延拓网格，从Field模型数据中读取历经延拓网格的节
点位置三分量及节点扰动引力三分量，完成扰动引力局域重构模型(Channel模型)构建。

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

第八步，建立飞行过程中扰动引力实时逼近算法

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

建立子域局部曲线坐标系，在局部坐标系中描述单元节点和计算点坐标。子域Ω_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(ξ,η,ζ)为，

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

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

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

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

则有，

上述问题可描述为，

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

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

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

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

写成矩阵形式，即

其中，

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

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

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

第九步，获取代理模型训练样本计算条件

以t＝20，s＝32，dλ＝dφ的简化情况为例进行说明。待设计变量为dr、dφ，其取值
范围为5km≤dr≤15km，1°≤dφ≤20°。选取样本数为50，基于最优拉丁超立方实验设计方
法进行实验设计，设计结果如表6所示。

表6实验设计结果

第十步，建立重构模型的代理模型

根据第九步计算得到的不同水平下的设计因素组合Xⁱ＝(dr,dφ)开展考虑扰动
引力的滑翔弹道仿真。弹道仿真条件为：纵程6700km，初始方位角90°，再入点和落点位于特
大山区(λ₀＝80°，φ₀＝42°，λ_k＝147°，φ_k＝19.3°)。扰动引力根据第一步至第八步所述的
扰动引力重构方法计算得到。

经弹道仿真，可以获取不同设计因素组合Xⁱ＝(dr,dφ)对应的Channel模型节点
数N_ch-node以及由扰动引力引起的滑翔弹道终端位置偏差dD，见表7。

表7设计因素组合对应的节点数和弹道终端位置偏差

以Xⁱ＝(dr,dλ,dφ,t,s)作为训练代理模型的输入对，以Yⁱ＝(N,dD)作为训练代理
的输出对，通过RBF神经网络建立描述输入对与输出对之间函数关系的代理模型，见图10。

输入4组随机选择的设计变量，通过将代理模型与原始模型的计算结果进行比较，
计算代理模型的拟合精度。代理模型对弹道终端位置偏差的拟合误差见表8，代理模型对节
点数的拟合误差见表9。可以看出，代理模型具有较高的拟合精度高，终端位置偏差的拟合
误差可以控制在0.5％以内，节点数基本无误差。

表8代理模型对弹道终端位置偏差的拟合误差

表9代理模型对节点数的拟合误差

第十一步，建立重构模型的优化算法

基于多岛遗传算法求解优化问题，其主要参数设置见表10。

表10多岛遗传算法参数设置

优化问题一的数学描述为：

目标函数：minN；

待设计变量：dr、dφ；

约束条件：0m≤dD≤5m；

待设计变量取值范围：5km≤dr≤15km，1°≤dφ≤20°。

经875次计算寻找到满足约束条件的最优网格划分。最优重构模型对应的dr为
10.040km，dφ为3.349°，终端位置偏差为4.570m，节点数为377。由此，建立了给定精度要求
下的存储量最小模型。

优化问题二的数学描述为：

目标函数：mindD；

待设计变量：dr、dφ；

约束条件：0＜N≤400；

待设计变量取值范围：5km≤dr≤15km，1°≤dφ≤20°。

经846次计算寻找到满足约束条件的最优网格划分。最优重构模型对应的dr为
10.189km，dφ为3.111°，终端位置偏差为4.030m，节点数为399。由此，建立了给定存储量要
求下的精度最高模型。

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

1)重构方法能快速完成沿滑翔弹道的扰动引力计算，具有计算精度高且弹上数据
存储量小的特征。

2)代理模型具有较高的拟合精度，终端位置偏差的拟合误差能控制在0.5％以内，
节点数基本无误差；

3)通过本文建立的优化方法，可以获得任意给定精度要求下的存储量最小模型和
任意给定存储量要求下的精度最高模型。当终端位置偏差限定为5m时，最优模型对应的弹
上数据存储量仅为377，可以满足弹上数据存储量的要求。当节点数限定400时，最优模型对
应的弹道终端位置偏差仅为4.030m，可满足重构精度要求。

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

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