CN104392047B - 一种基于平稳滑翔弹道解析解的快速弹道规划方法 - Google Patents

Abstract

本发明公开了一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：步骤1：滑翔段弹道规划问题建模；步骤2：滑翔段弹道规划变量设计；步骤3：滑翔段弹道解析解求解；步骤4：滑翔段弹道终端速度控制方案；步骤5：滑翔段弹道再入走廊调整方案；步骤6：滑翔段弹道规划初值生成；步骤7：滑翔段弹道规划流程设计。本发明提出以纵向机动加速度比例系数和横向机动加速度比例系数为滑翔段弹道规划变量，使得运动方程中弹道倾角、弹道偏角、高度、经度与纬度的微分方程中不含速度项；本发明获得了固定纵向机动加速度比例系数和横向机动加速度比例系数对应的滑翔段弹道解析解。

Description

一种基于平稳滑翔弹道解析解的快速弹道规划方法

技术领域

本发明涉及一种基于平稳滑行弹道解析解的快速弹道规划方法，属于航天技术、武器技术领域。

背景技术

随着高超声速技术的快速发展，再入弹道规划成为了研究热点。滑翔段是再入弹道的主要构成部分，决定了再入飞行的航程和机动范围。滑翔段弹道快速规划不仅可以分析高超声速飞行器的性能，还可用于在线弹道规划和预测制导，具有较高的研究价值。

由于滑翔段弹道具有飞行时间长，对规划变量高度敏感，且再入过程非线性约束强，使得滑翔弹道的可行域较窄，采用传统优化算法求解时往往需要很长的计算时间。尽管伪谱法、SQP等优化算法在计算效率上有较大的提高，其计算耗时依然在秒级以上，并且当过程约束碰边界时，其计算效率会大大降低，甚至影响到结果的收敛性。为了提高滑翔段弹道规划的速度，一些特殊的约束被引入滑翔段弹道规划，如阻力曲线、等热流曲线、等动压曲线、平衡滑翔条件等。这些约束的共同特点是将滑翔段弹道限定在某种特殊的平稳滑翔状态，从而降低了弹道规划问题的敏度。但由于在弹道规划时还需要进行数值积分，限制了规划速度的提高。为了摆脱对弹道的积分的依赖，动态逆被引入滑翔段弹道规划。该方法直接规划几何弹道形状，并通过调整弹道形状来满足约束要求。但由于滑翔段前后的环境变化剧烈，使得滑翔段全程弹道难以用低阶曲线拟合。

发明内容

本发明的目的是通过求解平稳滑翔弹道高精度解析解，获得不依赖弹道积分的滑翔段弹道快速规划方法，为分析高超声速飞行器性能和预测制导提供技术支持。

一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：

步骤1：滑翔段弹道规划问题建模；

步骤2：滑翔段弹道规划变量设计；

步骤3：滑翔段弹道解析解求解；

步骤4：滑翔段弹道终端速度控制方案；

步骤5：滑翔段弹道再入走廊调整方案；

步骤6：滑翔段弹道规划初值生成；

步骤7：滑翔段弹道规划流程设计。

本发明的优点在于：

(1)提出以纵向机动加速度比例系数和横向机动加速度比例系数为滑翔段弹道规划变量，使得运动方程中弹道倾角、弹道偏角、高度、经度与纬度的微分方程中不含速度项；

(2)获得了固定纵向机动加速度比例系数和横向机动加速度比例系数对应的滑翔段弹道解析解，包括高度解析解、飞行距离解析解、弹道偏角解析解、经度解析解和纬度解析解；

(3)给出了纵向机动加速度比例系数和横向机动加速度比例系数的初值生成方法；

(4)提出以纵向机动加速度比例系数、横向机动加速度比例系数、横向反转弹道倾角和初始弹道倾角增量为滑翔段弹道规划变量，其中：纵向机动加速度比例系数对应滑翔段飞行距离；横向机动加速度比例系数对应终端速度；横向反转弹道倾角对应终端经纬度；初始弹道倾角增量对应过程约束，上述四者之间相互独立。

(5)本发明采用滑翔弹道解析解规确定弹道，而积分仅仅用于终端速度校正，使得弹道规划具有非常快的速度。

附图说明

图1是平稳滑翔弹道规划建模流程图；

图2是大圆坐标系；

图3是与V_f的关系；

图4是Δγ对过程约束的影响；

图5是球面飞行距离估计；

图6是k_ε与s_f的关系；

图7是滑翔段弹道快速规划方法流程；

图8是纵向弹道解析解精度验证；

图9是横向弹道解析解精度验证；

图10是Bell解析解；

图11是可达区域边界弹道攻角曲线族；

图12是可达区域边界弹道倾侧角曲线族；

图13是可达区域边界横向弹道曲线族；

图14是可达区域边界纵向弹道曲线族；

图15是可达区域边界弹道速度曲线族；

图16是可达区域边界弹道热流密度曲线族；

图17是可达区域边界弹道动压曲线族；

图18是可达区域边界弹道过载曲线族；

图19是滑翔段弹道规划可达区域；

图20是平面规划弹道；

图21是规划弹道攻角曲线；

图22是规划弹道倾侧角曲线。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明在分析运动方程的耦合条件的基础上，提出了一种以固定纵向机动加速度比例系数和横向机动加速度比例系数为控制变量的弹道模式。该弹道模式下，速度与其他状态变量之间解耦，从而获得了解耦的弹道形状微分方程，并利用该微分方程求解得到了高精度的三维弹道形状解析解。此基础上，构造了滑翔段弹道的快速规划方法：首先利用弹道形状解析解规划出满足位置约束的弹道；然后通过调整横向机动大小来满足终端速度约束；最后通过调整初始弹道倾角来满足过程约束要求。

本发明是一种基于平稳滑行弹道解析解的快速弹道规划方法，规划方法建模流程如图1所示，规划流程如图7所示，包括以下几个步骤：

步骤1：滑翔段弹道规划问题建模

为了便于理论分析，假设地球为圆球体，不考虑地球自转的影响，则半速度坐标系下的三自由度质点运动方程如下所示，

式中，h为飞行器的高度，r为从地心到飞行器质心的距离，两者之间的关系为h＝r-R₀，其中R₀为地球半径；θ、φ分别为经度和纬度；s为飞行距离；V为相对地球的速度；γ为当地弹道倾角；ψ为弹道偏角。和分别为高度、飞行距离、经度、纬度、弹道倾角、弹道偏角和速度对时间的导数。σ为倾侧角；g＝μ/r²为当地重力加速度，μ为地球重力加速度相关的常数；L和D分别为升力加速度和阻力加速度，其表达式为，

L＝ρV²S_refC_L/(2m) (2)

D＝ρV²S_refC_D/(2m) (3)

式中，m为飞行器质量；S_ref为气动参考面积；C_L和C_D分别为升力系数和阻力系数，与马赫数Ma和攻角α相关；ρ为大气密度，通常采用指数大气模型，如下所示，

ρ＝ρ_seae^-βh (4)

式中，ρ_sea为海平面大气密度；β为指数大气模型常数，通常取β＝1/7200m^-1。

在滑翔段的弹道规划中还需考虑热流密度、动压、过载等过程约束，如下所示

式中，为热流密度，为最大热流密度，k为常系数；q为动压，q_max为最大动压；n为总过载，n_max为最大总过载，g₀＝9.81m/s²。

步骤2：滑翔段弹道规划变量设计

引入纵向机动加速度a_ε和横向机动加速度a_β，表达式如下，

a_ε＝Lcosσ+(V²/r-g)cosγ (6)

a_β＝Lsinσ/cosγ+(V²/r)cosγsinψtanφ (7)

式中，a_ε和a_β分别为纵向机动加速度和横向机动加速度。将式(6)和式(7)分别带入式(1)，则弹道倾角微分方程和弹道偏角微分方程可以化为，

假设a_ε＜0，则滑翔段的γ单调递减，因此可将运动方程中自变量替换为γ。将式(1)除以式(8)中的弹道倾角微分方程，得

式中，和分别为高度、飞行距离、弹道偏角、经度、纬度和速度对γ的微分。其中，第一行为滑翔段纵向运动微分方程，第二行为滑翔段横向运动微分方程，第三行为速度微分方程。

由式(9)可以看出，若a_ε和a_β均与V²成正比，则可在纵向运动微分方程和横向运动微分方程中均不显含速度项；同时由于θ和φ的微分方程分别与sinψ和cosψ相关，若要进行解析积分，则ψ的表达式越简单越好，因此a_ε和a_β之间也成正比。另外，由式(7)可以看出，为了使规划弹道的倾侧角变化幅度较小，a_β还需与ρ成正比。总的来说，可设滑翔段的a_ε与a_β形式如下，

式中，k_ε和k_β分别为纵向机动加速度比例系数和横向机动加速度比例系数。将式(10)带入式(9)，并取sinγ＝γ、cosγ＝1(在滑翔段γ为小量)，则有

由式(11)至式(12)可知，只要给定k_ε和k_β，就可对上述微分方程进行解析求解。因此以k_ε和k_β作为规划变量，则可以将滑翔段弹道规划问题转化为求解非线性方程的问题。

步骤3：滑翔段弹道解析解求解；

设滑翔段的初始弹道倾角为γ₀，初始高度为h₀，初始弹道偏角为ψ₀、初始经度和纬度为分别为θ₀和φ₀；滑翔段的终点高度为h_f，对于给定的k_ε和k_β，则可对式(11)至(15)式进行解析求解。

第一、高度解析解

将式(4)带入式(11)可得，

对式(16)积分可得，

式中，γ₀为初始弹道倾角；ρ₀为滑翔段初始大气密度。对式(17)进一步整理可得，

式中，C₁为一常数，取值为

第二、飞行距离解析解

将式(18)带入式(12)可得，

设初始飞行距离s₀＝0，则式(19)积分可得，

s＝f_s(γ)-f_s(γ₀) (20)

式中，f_s(γ)为由式(19)不定积分获得的函数，其表达式与C₁的正负相关，如下所示，

式中，x为自变量。

第三、道偏角解析解

进一步，由式(13)积分可得，

ψ＝ψ₀+(k_β/k_ε)(γ-γ₀) (21)

式中，ψ₀为初始弹道偏角。当滑翔段弹道位于赤道附近时，ψ的取值在π/2附近；若弹道不在赤道附近，则可建立滑翔段起点和终点决定的大圆坐标系，则在大圆坐标系下的弹道偏角ψ还在π/2附近。(如图2所示，坐标原点位于地心，x_e1轴由地心指向滑翔起点，y_e1位于由地心、滑翔起点、滑翔终点决定的大圆内，且与初始速度方向的夹角小于90deg，z_e1由右手定则决定。)

第四、纬度解析解

将cosψ在π/2附近5阶泰勒展开，可得，

取y＝π/2-ψ，其中y为积分中间变量，并将式(18)、式(21)和式(22)带入式(15)可得，

式中，C₂为一常量，满足C₂＝γ₀+k_ε(π/2-ψ₀)/k_β。取，

式中，C_a0、C_a1和C_a2均为常量。将式(24)带入式(23)可得，

进一步取，

式中，C_b1、C_b2、C_b3和C_b4均为常量。利用式(26)可将式(25)化为，

对式(27)积分可得，

φ＝φ₀+f_φ(y)-f_φ(y₀) (28)

式中，φ₀为初始纬度；y₀与初始弹道偏角相关，取值为y₀＝π/2-ψ₀；f_φ(y)为与y相关的函数，如下所示

式中，g_φ(y)为与y相关的函数，如下所示

第五、经度解析解

式(18)、式(21)和式(29)分别给出了滑翔段弹道的高度、弹道偏角和纬度随弹道倾角变化规律的解析解，将它们分别代入式(14)可得，

可以看出，式(30)的右边仅与γ相关，但复杂度较高，难以解析求解，因此可利用Guass-Legendre求积公式直接获得结果，具体表达式如下，

式中，γ_i为[γ,γ₀]中的第i个高斯节点；A_i为求积系数；n为高斯节点总数。

第六、解析解精度分析

在上面的解析解中，弹道偏角解析解(式21)无任何假设，属于准确解；高度解析解(式18)的误差来源于sinγ一阶泰勒展开，通常γ＜2deg，因此高度解析解的相对误差大约为10^-4；飞行距离解析解(式20)并未引入新的误差项，因此相对误差与高度相同；纬度解析解(式29)中引入了cosψ的五阶泰勒展开，若ψ的变化范围为[0,π]，则展开误差最大为0.0047，因此纬度解析解的相对误差介于0.1％至1％；而求解经度解析解时使用了纬度解析解，并且还存在求积误差，因此其相对误差最大，在1％左右。

步骤4：滑翔段弹道终端速度控制方案；

在步骤3中，给出了给定k_ε和k_β条件下的滑翔段弹道解析解，但没有给出速度的解析解。为了获得终端速度的大小，需要利用上述解析解对速度进行数值积分。

设负比能量e＝μ/r-V²/2，则由式(1)可得，

式中，为负比能量的导数。由式(8)、式(10)和式(32)可得，

由式(33)可以看出，C_D是速度求解的关键。通过情况下，C_D可写成C_L和e的函数，如下所示，

C_D＝f_CD(C_L,e) (34)

式中，f_CD为升力系数与阻力系数的关系函数。求解式(33)积分的关键为求解C_L，由式(6)、(7)式和式(10)可得，

另外，由式(4)和式(18)可得，

由(35)、式(36)和式(37)可求得C_L，将C_L带入式(34)可得C_D的大小。

最终，将r和C_D的表达式带入式(33)，并采用龙格-库塔法进行积分，可获得终端负比能量，从而解得终端的速度大小。需要指出的是，由式(33)易知C_D越大则终端速度越小，而C_L越大则C_D越大，因此调整C_L可以控制终端速度。由于k_ε通常为一小量，还决定了射程的大小，因此不能用于速度调节；因而只能采用k_β来调整终端速度。图3给出了与V_f的关系，可以看出两者之间的线性化程度非常高，有利于弹道规划时的数值求解。

当规划的弹道包含倾侧反转点时，则反转点前后的k_β需要满足如下连接条件，

式中，k_β-和k_β+分别为反转点前后的横向机动加速度比例系数；ρ_c、ψ_c和φ_c分别为γ_c处的大气密度、弹道偏角和纬度。

步骤5：滑翔段弹道再入走廊调整方案；

采用初始弹道倾角调整再入弹道在走廊中的位置。由图4可知，初始弹道倾角γ₀越大，则再入弹道的高度越高，距离由最大热流密度约束、最大动压约束和最大过载约束确定的再入走廊下边界越远。但考虑到滑翔段攻角的平滑性，初始弹道倾角取值如下所示，

γ₀＝γ^*+Δγ (39)

式中，Δγ为调整弹道再入走廊位置的初始弹道倾角分量；γ*为满足攻角平滑性要求的初始弹道倾角分量，满足，

γ^*＝-2g/(K_NβV₀ ²) (40)

式中K_N＝Lcosσ/D，为纵向升阻比。需要指出的是，由于本发明规划的滑翔段弹道满足平稳滑翔条件，因而热流密度、动压和过载变化较为平缓，峰值也较小，因此Δγ通常为零。

步骤6：滑翔段弹道规划初值生成；

在进行滑翔段弹道规划前，需要对s_f、γ^*、k_ε和k_β的初值进行估算。

第一、s_f估算

球面坐标系下，(θ₀,φ₀)和(θ_f,φ_f)的距离为(如图5所示)，

式中，s_s为(θ₀,φ₀)至(θ_f,φ_f)的球面距离；和分别满足，

考虑到初始弹道偏角偏离射面的情况，终端飞行距离为，

式中，其中为滑翔段起点子午线与射面所在大圆的夹角。

第二、γ^*估算

在滑翔段，纵向升阻比与终端飞行距离存在如下关系，

将式(42)带入式(40)，可以求得γ^*的估值为，

第三、k_ε估算

设Δγ＝0，则γ₀＝γ^*。由式(17)和式(20)可知，γ_f、k_ε和s_f满足如下关系，

式中，γ_f和ρ_f分别为终端弹道倾角和终端大气密度。由式(44)可知，给定s_f，则存在唯一的k_ε与之对应(如图6所示)。将式(41)得到的s_f带入式(44)求解即可获得k_ε的估值。

第四、k_β估算

由图2可知，在k_β较小时，与V_f的线性关系较好。因此在估算k_β初值时的原则为宁小勿大，这样有利于之后进行迭代修正。由式(21)知，在已知(γ-γ₀)/k_ε的条件下，k_β仅与ψ_f-ψ₀相关。假设ψ_f-ψ₀仅用于修正ψ₁₀的影响，则根据平面圆弧弹道假设有，

ψ_f-ψ₀＝2ψ₁₀

从而可得k_β的初始估值为，

k_β＝2ψ₁₀C_N1/(γ-γ₀) (45)

步骤7：滑翔段弹道规划流程设计；

滑翔段弹道规划任务为寻找合适的k_ε、k_β、Δγ和倾侧角反转点γ_c，使得再入飞行器从滑翔起点(θ₀,φ₀,r₀,V₀,ψ₀)至滑翔终点(θ_f,φ_f,r_f,V_f)。上述四个规划变量分别对应着不同的约束要求，彼此之间几乎不耦合，其中：k_ε用于调整滑翔段的飞行距离；k_β用于调整滑翔段终点的速度；γ_c用于调整滑翔段的终点位置；Δγ将用于调整再入弹道高度，以满足过程约束要求。利用滑翔段弹道形状解析解，可将滑翔段弹道规划分拆成弹道形状规划和终端速度校正两部分。其中，弹道形状规划不需要进行弹道积分，解决终端位置约束问题；终端端速度校正则需要对速度微分方程进行积分，解决终端速度约束和过程约束问题。完整的滑翔段弹道流程如图7所示，具体说明如下：

第一、根据步骤6，获取s_f、γ^*、k_ε和k_β的规划初值，并假设Δγ＝0和γ_c＝γ₀。

第二、根据给定的Δγ、γ^*、γ_c、k_ε和k_β，利用步骤3计算θ(γ_f)和φ(γ_f)。其中，γ_c为反转点弹道倾角；θ(γ_f)和φ(γ_f)为规划弹道终端经度和纬度。

第三、修正γ_c，使得终端纬度满足要求，修正的方法为牛顿迭代法，迭代的终止条件为|φ(γ_f)-φ_f|＜φ_limit。式中φ_limit为纬度规划精度，考虑步骤3中的解析弹道计算精度，可取φ_limit＝0.01deg。若|φ(γ_f)-φ_f|≥φ_limit则转向步骤7中的第二步。

第四、修正s_f(或者k_ε)使得终端经度满足要求，修正方法为牛顿迭代法，迭代终止条件为|θ(γ_f)-θ_f|＜θ_limit。式中θ_limit为经度规划精度，考虑到步骤3中的解析弹道计算精度，可取θ_limit＝0.1deg。若|θ(γ_f)-θ_f|≥θ_limit则转向步骤7中的第二步。

第五、根据Δγ、γ^*、γ_c、k_ε和k_β，利用步骤4积分得到滑翔段弹道终端速度V(γ_f)，并采用牛顿迭代法修正k_β，使得|V(γ_f)-V_f|＜V_limit。式中V_limit为速度规划精度，考虑到步骤3解析弹道计算精度，可取V_limit＝10m/s。若|V(γ_f)-V_f|≥V_limit则转向步骤7中的第二步。

第六，根据上一步的积分弹道判断是否满足过程约束，并对Δγ进行如下修正，

式中，h为实际飞行高度；h_limit为由最大热流密度约束、最大动压约束和最大过载约束确定的再入走廊高度下边界；min(h-h_limit)＜0表示滑翔段弹道超出再入走廊，需要增大Δγ，并转向步骤7中的第二步；min(h-h_limit)≥0表示滑翔段弹道满足再入走廊要求；γ_lim为min(h-h_limit)取极值时的弹道倾角；Δγ^(old)和Δγ^(new)分别为修正前和修正后的初始弹道倾角修正量；h(γ_lim)、s(γ_lim)和h_limit(γ_lim)分别为弹道倾角取γ_lim时的高度、飞行距离和高度下边界。

实施案例：

为了检验解析求解算法的精度，选用CAV作为计算模型，进行数值仿真效验。仿真采用无转动圆球模型和指数大气模型，参数设置如表格1。仿真计算机CPU为Core i3-2120，内存4GB，仿真环境为Matlab。

表格1仿真参数设置

第一，解析解精度效验

在进行弹道规划前，首先对弹道形状解析解精度进行验证，并与Bell解析解进行对比。图8和图9分别给出了纵向弹道解析解与横向弹道解析解与弹道积分结果的对比，可以看出，弹道形状解析解与弹道积分结果几乎完全重合，具有极高的精度。具体来说，采用解析解直接预测滑翔段终点位置，其高度误差小于1m；弹道偏角误差小于0.01deg；大圆坐标系下的纬度误差小于0.1deg，经度误差小于0.3deg；验证了文中的误差分析。图10给出了Bell给出的横向弹道解析解，可以看出Bell解析解只在较小的范围内具有较高的精度，当飞行距离较远时，横程偏差较大。这是由于在Bell解析解中忽略了地球曲率对横向弹道的影响，而采用平面坐标系进行求解，带来了较大的求解误差。

第二，弹道规划可达区域分析

在验证了弹道形状解析解的精度之后，还需分析弹道规划的可达区域。在求解可达区域时，假定弹道偏角为固定值，因而主要通过倾侧反转来调节终端位置，因此无反转的弹道即为弹道规划可达区域的边界。图11至图18给出了边界弹道曲线族。由图11和图12可知，规划弹道的攻角曲线和倾侧角曲线非常光滑且变化幅度较小；由图13可知滑翔段弹道规划的可达区域纵程覆盖范围约为11000千米，横程覆盖范围约为6700千米；由图14和图15可知规划弹道的终端高度和终端速度满足要求；由图16至图18可知规划弹道满足过程约束要求。对于目标在可达区域内的情况，由于仅仅增加了倾侧角反转点来调节终端经纬度，因此其攻角曲线、倾侧曲线、纵向弹道曲线和过程约束曲线的趋势将与图11至图18中给出的结果类似。

在图13的基础上，图19进一步获得了滑翔段规划弹道的可达区域。与平稳滑翔最大可达区域相比，本发明获得的可达区域的要小一些，最大射程损失约为81.7km，最小射程损失约为302.8km，最大横程损失约为923.9km。这是由于本发明选用固定的k_β作为横向机动规划变量，使得在大机动情况下倾侧角存在较大的变化(如图11所示)，这会增加横程机动的能量损耗，从而降低横向机动能力。当射程较小或者横程较大时，需要较大的横向机动，因而出现图18所示的可达区域损失。

第三，滑翔段弹道快速规划

弹道规划的初值和终值如前所示，采用图5给出的滑翔段弹道规划流程，对图19给出的可达区域内的目标点进行弹道规划，得到的结果如表格1和图20至图22所示。由表格2可以看出，基于反馈线性化解析解的滑翔段弹道规划方法能够在0.03s左右的时间内规划出一条满足约束要求的再入弹道，并且弹道终端具有较高的位置精度。不足在于规划弹道的终端速度误差较大，这是由于终端速度是根据几何弹道、采用动态逆方法积分获得的，使得弹道形状上的微小偏差通过求解速度的积分得以放大。图20给出了规划弹道的形状；图21给出出了规划弹道的攻角曲线，可以看出攻角变化较为光滑，变化幅度小于4deg；图22给出了倾侧角曲线，同样也较为光滑。

表格2滑翔段弹道规划结果

Claims (1)

1.一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：

步骤1：滑翔段弹道规划问题建模

设地球为圆球体，不考虑地球自转，则半速度坐标系下的三自由度质点运动方程如下所示，

$\begin{array}{c}\begin{array}{cc}\stackrel{\·}{h}=V\sin \mathrm{\γ}& \stackrel{\·}{\mathrm{\θ}}=V\cos \mathrm{\γ}\sin \mathrm{\ψ}/\left(r\cos \mathrm{\φ}\right)\\ \stackrel{\·}{s}=V{\mathrm{cos\γR}}_0/r& \stackrel{\·}{\mathrm{\φ}}=V\cos \mathrm{\γ}\cos \mathrm{\ψ}/r\\ \stackrel{\·}{\mathrm{\γ}}=\[L\cos \mathrm{\σ}+(V^2/r-g)\cos \mathrm{\γ}\]/V& \stackrel{\·}{\mathrm{\ψ}}=\[L\sin \mathrm{\σ}/\cos \mathrm{\γ}+(V^2/r)\cos \mathrm{\γ}\sin \mathrm{\ψ}\tan \mathrm{\φ}\]/V\end{array}\\ \stackrel{\·}{V}=-D-g\sin \mathrm{\γ}\end{array}---\left(1\right)$

式中，h为飞行器的高度，r为从地心到飞行器质心的距离，h＝r-R₀，其中R₀为地球半径；θ、φ分别为经度和纬度；s为飞行距离；V为相对地球的速度；γ为当地弹道倾角；ψ为弹道偏角；和分别为高度、飞行距离、经度、纬度、弹道倾角、弹道偏角和速度对时间的导数；σ为倾侧角；g＝μ/r²为当地重力加速度，μ为地球重力加速度相关的常数；L和D分别为升力加速度和阻力加速度，其表达式为，

L＝ρV²S_refC_L/(2m) (2)

D＝ρV²S_refC_D/(2m) (3)

式中，m为飞行器质量；S_ref为气动参考面积；C_L和C_D分别为升力系数和阻力系数，与马赫数Ma和攻角α相关；ρ为大气密度，如下所示，

ρ＝ρ_seae^-βh (4)

式中，ρ_sea为海平面大气密度；β为指数大气模型常数；

过程约束如下所示

$\begin{array}{ccc}\stackrel{\·}{Q}=k\sqrt{\mathrm{\ρ}}{V}^{3.15}\≤{\stackrel{\·}{Q}}_{max}& q=\frac{1}{2}{\mathrm{\ρV}}^2\≤q_{\max }& n=\sqrt{L^2+D^2}/g_0\≤n_{\max }\end{array}---\left(5\right)$

式中，为热流密度，为最大热流密度，k为常系数；q为动压，q_max为最大动压；n为总过载，n_max为最大总过载，g₀＝9.81m/s²；

步骤2：滑翔段弹道规划变量设计

引入纵向机动加速度a_ε和横向机动加速度a_β，表达式如下，

a_ε＝Lcosσ+(V²/r-g)cosγ (6)

a_β＝Lsinσ/cosγ+(V²/r)cosγsinψtanφ (7)

式中，a_ε和a_β分别为纵向机动加速度和横向机动加速度；将式(6)和式(7)分别带入式(1)，则弹道倾角微分方程和弹道偏角微分方程为，

$\begin{array}{cc}\stackrel{\·}{\mathrm{\γ}}=a_{\mathrm{\ϵ}}/V& \stackrel{\·}{\mathrm{\ψ}}=a_{\mathrm{\β}}/V\end{array}---\left(8\right)$

假设a_ε＜0，则滑翔段的γ单调递减，将运动方程中自变量替换为γ；将式(1)除以式(8)中的弹道倾角微分方程，得

$\begin{array}{c}\begin{array}{cc}\frac{dh}{d\mathrm{\γ}}=\frac{V^2\sin \mathrm{\γ}}{a_{\mathrm{\ϵ}}}& \frac{ds}{d\mathrm{\γ}}=\frac{R_0{V}^2\cos \mathrm{\γ}}{{\mathrm{ra}}_{\mathrm{\ϵ}}}\end{array}\\ \begin{array}{ccc}\frac{d\mathrm{\ψ}}{d\mathrm{\γ}}=\frac{a_{\mathrm{\β}}}{a_{\mathrm{\ϵ}}}& \frac{d\mathrm{\θ}}{d\mathrm{\γ}}=\frac{V^2\cos \mathrm{\γ}\sin \mathrm{\ψ}}{{\mathrm{ra}}_{\mathrm{\ϵ}}\cos \mathrm{\φ}}& \frac{d\mathrm{\φ}}{d\mathrm{\γ}}=\frac{V^2\cos \mathrm{\γ}\cos \mathrm{\ψ}}{{\mathrm{ra}}_{\mathrm{\ϵ}}}\end{array}\\ \frac{dV}{d\mathrm{\γ}}=-\frac{(D+g\sin \mathrm{\γ})V}{a_{\mathrm{\ϵ}}}\end{array}---\left(9\right)$

式中，和分别为高度、飞行距离、弹道偏角、经度、纬度和速度对γ的微分；其中，第一行为滑翔段纵向运动微分方程，第二行为滑翔段横向运动微分方程，第三行为速度微分方程；

设滑翔段的a_ε与a_β形式如下，

$\begin{array}{cc}{a}_{\mathrm{\ϵ}}=\frac{{\mathrm{\ρV}}^2{S}_{ref}{k}_{\mathrm{\ϵ}}}{2m}& a_{\mathrm{\β}}=\frac{{\mathrm{\ρV}}^2{S}_{ref}{k}_{\mathrm{\β}}}{2m}\end{array}---\left(10\right)$

式中，k_ε和k_β分别为纵向机动加速度比例系数和横向机动加速度比例系数；将式(10)带入式(9)，并取sinγ＝γ、cosγ＝1，则有

$\frac{dh}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\mathrm{\γ}---\left(11\right)$

$\frac{ds}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\frac{R_0}{r}---\left(12\right)$

$\frac{d\mathrm{\ψ}}{d\mathrm{\γ}}=\frac{k_{\mathrm{\β}}}{k_{\mathrm{\ϵ}}}---\left(13\right)$

$\frac{d\mathrm{\θ}}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\frac{sin\mathrm{\ψ}}{rcos\mathrm{\φ}}---\left(14\right)$

$\frac{d\mathrm{\φ}}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\frac{cos\mathrm{\ψ}}{r}---\left(15\right)$

以k_ε和k_β作为规划变量，将滑翔段弹道规划问题转化为求解非线性方程的问题；

步骤3：滑翔段弹道解析解求解；

设滑翔段的初始弹道倾角为γ₀，初始高度为h₀，初始弹道偏角为ψ₀、初始经度和纬度为分别为θ₀和φ₀；滑翔段的终点高度为h_f，对于给定的k_ε和k_β，对式(11)至(15)式进行解析求解；

第一、高度解析解

将式(4)带入式(11)可得，

$e^{-\mathrm{\β}h}dh=\frac{2m}{{\mathrm{\ρ}}_{sea}{S}_{ref}{k}_{\mathrm{\ϵ}}}\mathrm{\γ}d\mathrm{\γ}---\left(16\right)$

对式(16)积分可得，

${\mathrm{\γ}}^2-{\mathrm{\γ}}_0^2=\frac{S_{ref}{k}_{\mathrm{\ϵ}}}{m\mathrm{\β}}({\mathrm{\ρ}}_0-\mathrm{\ρ})---\left(17\right)$

式中，γ₀为初始弹道倾角；ρ₀为滑翔段初始大气密度；对式(17)进一步整理可得，

$\frac{m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}=\frac{1}{\mathrm{\β}(C_1-{\mathrm{\γ}}^2)}---\left(18\right)$

式中，C₁为一常数，取值为

第二、飞行距离解析解

将式(18)带入式(12)可得，

$\frac{ds}{d\mathrm{\γ}}=\frac{-2R_0}{\mathrm{\β}r({\mathrm{\γ}}^2-C_1)}---\left(19\right)$

设初始飞行距离s₀＝0，则式(19)积分可得，

s＝f_s(γ)-f_s(γ₀) (20)

式中，f_s(γ)为由式(19)不定积分获得的函数，其表达式与C₁的正负相关，如下所示，

$f_s\left(x\right)=\left\{\begin{array}{cc}-\frac{2R_0}{\mathrm{\&beta;}r\sqrt{-C_1}}arctan\left(\frac{x}{\sqrt{-C_1}}\right)& C_1<0\\ 2R_0/\left(\mathrm{\&beta;}rx\right)& C_1=0\\ -\frac{R_0}{\mathrm{\&beta;}r\sqrt{C_1}}ln\left(\frac{x-\sqrt{C_1}}{x+\sqrt{C_1}}\right)& C_1>0\end{array}\right.$

式中，x为自变量；

第三、道偏角解析解

进一步，由式(13)积分可得，

ψ＝ψ₀+(k_β/k_ε)(γ-γ₀) (21)

式中，ψ₀为初始弹道偏角；

第四、纬度解析解

将cosψ在π/2附近5阶泰勒展开，可得，

$cos\mathrm{\&psi;}=(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&psi;})-\frac{1}{6}{(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&psi;})}^3+\frac{1}{120}{(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&psi;})}^5---\left(22\right)$

取y＝π/2-ψ，其中y为积分中间变量，并将式(18)、式(21)和式(22)带入式(15)可得，

$\frac{d\mathrm{\&phi;}}{dy}=\left(\frac{k_{\mathrm{\&epsiv;}}}{{\mathrm{r\&beta;k}}_{\mathrm{\&beta;}}}\right)\frac{2y-y^3/3+y^5/60}{{\&lsqb;C_2-(k_{\mathrm{\&epsiv;}}/k_{\mathrm{\&beta;}})y\&rsqb;}^2-C_1}---\left(23\right)$

式中，C₂为一常量，满足C₂＝γ₀+k_ε(π/2-ψ₀)/k_β；取，

$\begin{array}{c}\begin{array}{cc}{C}_{a0}=k_{\mathrm{\&beta;}}/\left(60{\mathrm{r\&beta;k}}_{\mathrm{\&epsiv;}}\right)& C_{a1}=2C_2{k}_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}}\end{array}\\ C_{a2}=(C_2^2-C_1){(k_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}})}^2\end{array}---\left(24\right)$

式中，C_a0、C_a1和C_a2均为常量；将式(24)带入式(23)可得，

$\frac{d\mathrm{\&phi;}}{dy}=\frac{C_{a0}(y^5-20y^3+120y)}{y^2-C_{a1}y+C_{a2}}---\left(25\right)$

进一步取，

$\begin{array}{c}{C}_{b1}=C_{a1}^2-C_{a2}-20\\ C_{b2}=C_{a1}^3-2C_{a1}{C}_{a2}-2C_{a1}\\ C_{b3}=C_{a1}^4-3C_{a1}^2{C}_{a2}-20C_{a1}^2+C_{a2}^2+20C_{a2}+120\\ C_{b4}=-C_{a1}^3{C}_{a2}+2C_{a1}{C}_{a2}^2+20C_{a1}{C}_{a2}\end{array}---\left(26\right)$

式中，C_b1、C_b2、C_b3和C_b4均为常量；利用式(26)可将式(25)化为，

$\begin{array}{c}\frac{d\mathrm{\&phi;}}{dy}=C_{a0}(y^3+C_{a1}{y}^{2}y+C_{b1}y+C_{b2})\\ +\frac{(C_{a0}{C}_{b3}/2)(2y-C_{a1})}{y^2-C_{a1}y+C_{a2}}+\frac{C_{a0}(C_{b3}{C}_{a1}/2+C_{b4})}{y^2-C_{a1}y+C_{a2}}\end{array}---\left(27\right)$

对式(27)积分可得，

φ＝φ₀+f_φ(y)-f_φ(y₀) (28)

式中，φ₀为初始纬度；y₀与初始弹道偏角相关，取值为y₀＝π/2-ψ₀；f_φ(y)为与y相关的函数，如下所示

$\begin{array}{c}{f}_{\mathrm{\&phi;}}\left(y\right)=C_{a0}(\frac{y^4}{4}+\frac{C_{a1}{y}^3}{3}+\frac{C_{b1}{y}^2}{2}+C_{b2}y)\\ +\frac{C_{a0}{C}_{b3}}{2}\ln (y^2-C_{a1}y+C_{a2})\\ +C_{a0}(C_{b3}{C}_{a1}/2+C_{b4})g_{\mathrm{\&phi;}}\left(y\right)\end{array}---\left(29\right)$

式中，g_φ(y)为与y相关的函数，如下所示

$g_{\mathrm{\&phi;}}\left(y\right)=\left\{\begin{array}{cc}\frac{2}{\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}arctan\left(\frac{2y-C_{a1}}{\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}\right)& C_{a1}^2<4C_{a2}\\ -\frac{2}{2y-C_{a1}}& C_{a1}^2=4C_{a2}\\ \frac{1}{\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}ln\left(\frac{2y-C_{a1}-\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}{2y-C_{a1}+\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}\right)& C_{a1}^2>4C_{a2}\end{array}\right.$

第五、经度解析解

式(18)、式(21)和式(29)分别给出了滑翔段弹道的高度、弹道偏角和纬度随弹道倾角变化规律的解析解，将它们分别代入式(14)可得，

$\frac{d\mathrm{\&theta;}}{d\mathrm{\&gamma;}}=\frac{2sin\&lsqb;{\mathrm{\&psi;}}_0+(k_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}})(\mathrm{\&gamma;}-{\mathrm{\&gamma;}}_0)\&rsqb;}{\mathrm{\&beta;}r(C_1-{\mathrm{\&gamma;}}^2)cos\&lsqb;{\mathrm{\&phi;}}_0+f_{\mathrm{\&phi;}}\left(x\right)-f_{\mathrm{\&phi;}}\left(x_0\right)\&rsqb;}---\left(30\right)$

利用Guass-Legendre求积公式直接获得结果，具体表达式如下，

$\mathrm{\&theta;}={\mathrm{\&theta;}}_0+\frac{\mathrm{\&gamma;}-{\mathrm{\&gamma;}}_0}{2}\underset{i=1}{\overset{n}{\&Sigma;}}{A}_i{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_1\left({\mathrm{\&gamma;}}_i\right)---\left(31\right)$

式中，γ_i为[γ,γ0]中的第i个高斯节点；A_i为求积系数；n为高斯节点总数；

步骤4：滑翔段弹道终端速度控制方案；

设负比能量e＝μ/r-V²/2，则由式(1)可得，

$\stackrel{\&CenterDot;}{e}=VD---\left(32\right)$

式中，为负比能量的导数；由式(8)、式(10)和式(32)可得，

$\frac{de}{d\mathrm{\&gamma;}}=\frac{2(\mathrm{\&mu;}/r-e)C_D}{k_{\mathrm{\&epsiv;}}}---\left(33\right)$

由式(33)看出，C_D是速度求解的关键；通过情况下，C_D可写成C_L和e的函数，如下所示，

C_D＝f_CD(C_L,e) (34)

式中，f_CD为升力系数与阻力系数的关系函数；求解式(33)积分的关键为求解C_L，由式(6)、(7)式和式(10)可得，

$C_{L}cos\mathrm{\&sigma;}=\frac{m\&lsqb;g/(\mathrm{\&mu;}/r-e)-2/r\&rsqb;cos\mathrm{\&gamma;}}{{\mathrm{\&rho;S}}_{ref}}+k_{\mathrm{\&epsiv;}}---\left(35\right)$

$C_{L}sin\mathrm{\&sigma;}=k_{\mathrm{\&beta;}}-\frac{2m{\left(cos\mathrm{\&gamma;}\right)}^{2}sin\mathrm{\&psi;}tan\mathrm{\&phi;}}{{\mathrm{\&rho;rS}}_{ref}}---\left(36\right)$

另外，由式(4)和式(18)可得，

$r=R_0-\frac{1}{\mathrm{\&beta;}}ln\&lsqb;\frac{m\mathrm{\&beta;}(C_1-{\mathrm{\&gamma;}}^2)}{{\mathrm{\&rho;}}_0{S}_{ref}{k}_{\mathrm{\&epsiv;}}}\&rsqb;---\left(37\right)$

由(35)、式(36)和式(37)可求得C_L，将C_L带入式(34)可得C_D的大小；

最终，将r和C_D的表达式带入式(33)，并采用龙格-库塔法进行积分，获得终端负比能量，从而解得终端的速度大小；

当规划的弹道包含倾侧反转点时，则反转点前后的k_β需要满足如下连接条件，

$k_{\mathrm{\&beta;}-}+k_{\mathrm{\&beta;}+}=\frac{4m{\left({\mathrm{cos\&gamma;}}_c\right)}^2{\mathrm{sin\&psi;}}_c{\mathrm{tan\&phi;}}_c}{{\mathrm{\&rho;}}_c{\mathrm{rS}}_{ref}}---\left(38\right)$

式中，k_β-和k_β+分别为反转点前后的横向机动加速度比例系数；ρ_c、ψ_c和φ_c分别为γ_c处的大气密度、弹道偏角和纬度；

步骤5：滑翔段弹道再入走廊调整方案；

初始弹道倾角取值如下所示，

γ₀＝γ^*+Δγ (39)

式中，Δγ为调整弹道再入走廊位置的初始弹道倾角分量；γ^*为满足攻角平滑性要求的初始弹道倾角分量，满足，

${\mathrm{\&gamma;}}^{*}=-2g/\left(K_N{\mathrm{\&beta;V}}_0^2\right)---\left(40\right)$

式中K_N＝L cosσ/D，为纵向升阻比；

步骤6：滑翔段弹道规划初值生成；

在进行滑翔段弹道规划前，需要对s_f、γ^*、k_ε和k_β的初值进行估算；

第一、s_f估算

球面坐标系下，(θ₀,φ₀)和(θ_f,φ_f)的距离为，

$s_s=R_{0}arccos\&lsqb;{\left({\stackrel{\&RightArrow;}{r}}_0\right)}^T{\stackrel{\&RightArrow;}{r}}_f\&rsqb;$

式中，s_s为(θ₀,φ₀)至(θ_f,φ_f)的球面距离；和分别满足，

$\begin{array}{cc}{\stackrel{\&RightArrow;}{r}}_0=\left[\begin{array}{c}{\mathrm{cos\&phi;}}_0{\mathrm{cos\&theta;}}_0\\ {\mathrm{cos\&phi;}}_0{\mathrm{sin\&theta;}}_0\\ {\mathrm{sin\&phi;}}_0\end{array}\right]& {\stackrel{\&RightArrow;}{r}}_f=\left[\begin{array}{c}{\mathrm{cos\&phi;}}_f{\mathrm{cos\&theta;}}_f\\ {\mathrm{cos\&phi;}}_f{\mathrm{sin\&theta;}}_f\\ {\mathrm{sin\&phi;}}_f\end{array}\right]\end{array}$

则终端飞行距离为，

$s_f=\frac{{\mathrm{\&psi;}}_{10}{s}_s}{{\mathrm{sin\&psi;}}_{10}}---\left(41\right)$

式中，其中为滑翔段起点子午线与射面所在大圆的夹角；

第二、γ^*估算

在滑翔段，纵向升阻比与终端飞行距离存在如下关系，

$K_N=\frac{2s_f}{rln\&lsqb;(gr-V_f^2)/(gr-V_0^2)\&rsqb;}---\left(42\right)$

将式(42)带入式(40)，求得γ^*的估值为，

${\mathrm{\&gamma;}}^{*}=-\frac{grln\&lsqb;(gr-V_f^2)/(gr-V_0^2)\&rsqb;}{s_f{\mathrm{\&beta;V}}_0^2}---\left(43\right)$

第三、k_ε估算

设Δγ＝0，则γ₀＝γ^*；由式(17)和式(20)可知，γ_f、k_ε和s_f满足如下关系，

$\begin{array}{c}{\mathrm{\&gamma;}}_f^2-{\mathrm{\&gamma;}}_0^2=\frac{S_{ref}{k}_{\mathrm{\&epsiv;}}}{m\mathrm{\&beta;}}({\mathrm{\&rho;}}_0-{\mathrm{\&rho;}}_f)\\ s_f=f_s\left({\mathrm{\&gamma;}}_f\right)-f_s\left({\mathrm{\&gamma;}}_0\right)\end{array}---\left(44\right)$

式中，γ_f和ρ_f分别为终端弹道倾角和终端大气密度；由式(44)可知，给定s_f，则存在唯一的k_ε与之对应；将式(41)得到的s_f带入式(44)求解即可获得k_ε的估值；

第四、k_β估算

假设ψ_f-ψ₀仅用于修正ψ₁₀的影响，则根据平面圆弧弹道假设有，

ψ_f-ψ₀＝2ψ₁₀

k_β的初始估值为，

k_β＝2ψ₁₀C_N1/(γ-γ₀) (45)

步骤7：滑翔段弹道规划流程设计；

具体为：

第一、根据步骤6，获取s_f、γ^*、k_ε和k_β的规划初值，并假设Δγ＝0和γ_c＝γ₀；

第二、根据给定的Δγ、γ^*、γ_c、k_ε和k_β，利用步骤3计算θ(γ_f)和φ(γ_f)；其中，γ_c为反转点弹道倾角；θ(γ_f)和φ(γ_f)为规划弹道终端经度和纬度；

第三、修正γ_c，使得终端纬度满足要求，修正的方法为牛顿迭代法，迭代的终止条件为|φ(γ_f)-φ_f|＜φ_limit；式中φ_limit为纬度规划精度；若|φ(γ_f)-φ_f|≥φ_limit则转向步骤7中的第二步；

第四、修正s_f或者k_ε使得终端经度满足要求，修正方法为牛顿迭代法，迭代终止条件为|θ(γ_f)-θ_f|＜θ_limit；式中θ_limit为经度规划精度；若|θ(γ_f)-θ_f|≥θ_limit则转向步骤7中的第二步；

第五、根据Δγ、γ^*、γ_c、k_ε和k_β，利用步骤4积分得到滑翔段弹道终端速度V(γ_f)，并采用牛顿迭代法修正k_β，使得式中为速度规划精度；若|V(γ_f)-V_f|≥V_limit则转向步骤7中的第二步；

第六，根据上一步的积分弹道判断是否满足过程约束，并对Δγ进行如下修正，

${\mathrm{\&Delta;\&gamma;}}^{\left(new\right)}=\left\{\begin{array}{cc}{\mathrm{\&Delta;\&gamma;}}^{\left(old\right)}-\frac{h\left({\mathrm{\&gamma;}}_{\lim }\right)-h_{\lim it}\left({\mathrm{\&gamma;}}_{\lim }\right)}{s\left({\mathrm{\&gamma;}}_{\lim }\right)}& \min (h-h_{\lim it})<0\\ {\mathrm{\&Delta;\&gamma;}}^{\left(old\right)}& \min (h-h_{\lim it})\&GreaterEqual;0\end{array}\right.$

式中，h为实际飞行高度；h_limit为由最大热流密度约束、最大动压约束和最大过载约束确定的再入走廊高度下边界；min(h-h_limit)＜0表示滑翔段弹道超出再入走廊，需要增大Δγ，并转向步骤7中的第二步；min(h-h_limit)≥0表示滑翔段弹道满足再入走廊要求；γ_lim为min(h-h_limit)取极值时的弹道倾角；Δγ^(old)和Δγ^(new)分别为修正前和修正后的初始弹道倾角修正量；h(γ_lim)、s(γ_lim)和h_limit(γ_lim)分别为弹道倾角取γ_lim时的高度、飞行距离和高度下边界。