基于高精度纵、横程解析预测方法的平稳滑翔再入制导律

摘要

基于高精度纵、横程解析预测方法的平稳滑翔再入制导律，该方法有三大步骤：步骤1：规划快速下降段的制导方案；步骤2：规划平稳滑翔段的制导方案；步骤3：规划弹道下压段的制导方案；基于能量的纵向射程、横向射程和航向角的解析解，求解解析解的具体步骤如下：步骤1：建立飞行器再入问题的运动方程；步骤2：建立再入问题的过程约束和终端约束；步骤3：引入广义赤道和广义子午线；步骤4：运动方程线性化；步骤5：求解解析解；步骤6：验证上述解析解可行性。此制导方法能够导引飞行器飞行到距离目标一定距离的指定高度，并能够满足相应的过程约束，末速度约束和方位角误差约束，具有很好的鲁棒性。

说明书

技术领域

本发明属于航天技术、高超声速飞行器再入制导技术领域，涉及一种基于高精度纵、横 程解析预测方法的平稳滑翔再入制导律，尤其涉及一种新的基于能量-纵程及能量-横程高精 度解析预测方法的通用飞行器(CAV)的平稳滑翔再入制导方法。具体给出能满足过程约束， 终端的速度、高度、航向角误差和射程约束要求的再入制导方法。

背景技术

通用飞行器(CAV)是一种高超声速飞行器，它通过运载火箭或者弹道导弹助推到亚轨 道，之后再入到大气层中无动力飞行。CAV的飞行过程可以大致分成两个阶段：再入滑翔段 和末端制导段。CAV与运载器分离之后，开始再入滑翔段，CAV在再入过程中需要在满足热 流率约束、动压约束和过载约束等过程约束的前提下，进行横向机动，一直飞行到距离目标 一定距离和高度时，进入末制导段。

飞行器再入制导有两种形式：跟踪参考轨迹和预测校正制导。航天飞机再入制导是典型 的参考轨迹跟踪制导，航天飞机再入飞行前先设计一个阻力-速度参考剖面来满足射程要求， 其在实际飞行中通过调节倾侧角来跟踪阻力-速度参考剖面。当航向角误差超过预先设计的误 差门限曲线时，通过倾侧角翻转来实现横向射程控制。另一种形式的再入制导方法是预测校 正制导，是以消除实际轨道的预报落点和预定落点位置之间的偏差为目的的制导方法。关于 高超声速飞行器滑翔再入问题，还有很多研究给出解析形式的结果，这些结果大多假定升阻 比和弹道倾角为常数，或为某些系数的函数形式，并忽略地球曲率的影响。

发明内容

1、目的：本发明给出了一种新的基于能量-纵程及能量-横程高精度解析预测方法的通用 飞行器(CAV)的平稳滑翔再入制导方法，是一种高超声速飞行器的再入制导方法。在推导 过程中，通过将升阻比L/D和倾侧角设计成与能量有关的函数，并利用线性化方法得到一种 特殊类型的变系数线性系统。为求解此系统，提出一种新的基于谱分解的求解方法，然后根 据此求解方法得到飞行器再入飞行射程、横向射程和航向角的解析解。随后基于这些解析解， 给出新的再入制导方法，此制导方法能够导引飞行器飞行到距离目标一定距离的指定高度， 并能够满足相应的过程约束，末速度约束和方位角误差约束，具有很好的鲁棒性。

2、技术方案：在给出本发明的制导方法前，对于高超声速飞行器再入飞行问题，本发明 给出一种新的基于能量的纵向射程、横向射程和航向角的解析解，求解解析解的具体步骤如 下：

步骤1：建立飞行器再入问题的运动方程；

如附图1所示，在推导飞行器运动方程之前，首先定义两个坐标系。

地心赤道旋转(GER)坐标系：坐标系原点在地心，z_e轴指向北极；x_e和y_e轴是赤道平 面内相互垂直的两条轴线。GER坐标系以z_e为中心，以与地球自转角速度相同的角速度旋转， 地球自转角速度ω_e为7.292116×10^-5rad/s。

当地北东下(NED)坐标系：坐标系原点o为飞行器质心M和地心连线与地球表面的交 点，即飞行器在地面的垂直投影点；x轴指向当地的北向，y轴指向当地的东向，z轴铅垂向 下。

由于大气可以认为是相对地球静止的，选择在GER坐标系中处理再入问题较合适。将再 入飞行器看作质点，在球形地球模型下，飞行器运动方程为

$\frac{\mathrm{d\lambda }}{\mathrm{dt}}=\frac{V\cos \left(\gamma \right)\sin \left(\psi \right)}{(R_e+H)\cos \left(\text{φ}\right)}---\left(1\right)$

$\frac{\mathrm{d\phi }}{\mathrm{dt}}=\frac{V\cos \left(\gamma \right)\cos \left(\psi \right)}{(R_e+H)}---\left(2\right)$

$\frac{\mathrm{dH}}{\mathrm{dt}}=V\sin \left(\gamma \right)---\left(3\right)$

$\left(\begin{array}{c}\frac{\mathrm{dV}}{\mathrm{dt}}=-\frac{D}{m}\text{-gsin}\left(\gamma \right)+{{\omega }_e}^2(R_e+H){\cos }^2\left(\phi \right)\sin \left(\gamma \right)\\ -{{\omega }_e}^2(R_e+H)\sin \left(\phi \right)\cos \left(\phi \right)\cos \left(\gamma \right)\cos \left(\psi \right)\end{array}\right)---\left(4\right)$

$\left(\begin{array}{c}\frac{\mathrm{d\gamma }}{\mathrm{dt}}=\frac{1}{V}[\frac{L\cos \left(\sigma \right)}{m}-g\cos \left(\gamma \right)+\frac{V^2\cos \left(\gamma \right)}{R_e+H}+2V{\omega }_e\cos \left(\phi \right)\sin \left(\psi \right)\\ +{{\omega }_e}^2(R_e+H){\cos }^2\left(\phi \right)\cos \left(\gamma \right)+{{\omega }_e}^2(R_e+H)\sin \left(\phi \right)\cos \left(\phi \right)\sin \left(\gamma \right)\cos \left(\psi \right)]\end{array}\right)---\left(5\right)$

$\left(\begin{array}{c}\frac{\mathrm{d\psi }}{\mathrm{dt}}=\frac{1}{V}[\frac{L\sin \left(\sigma \right)}{m\cos \left(\gamma \right)}+\frac{V^2\cos \left(\gamma \right)\sin \left(\psi \right)\tan \left(\phi \right)}{R_e+H}+2V{\omega }_e\sin \left(\phi \right)\\ -\frac{{{\omega }_e}^2(R_e+H)\sin \left(\phi \right)\cos \left(\phi \right)\sin \left(\psi \right)}{\cos \left(\gamma \right)}-2V{\omega }_e\cos \left(\phi \right)\tan \left(\gamma \right)\cos \left(\psi \right)]\end{array}\right)---\left(6\right)$

在上述运动方程中，如图1所示，t是时间，λ是经度，φ是纬度，H是高度，V是对旋转地 球的相对速度，γ是飞行的弹道倾角，ψ是航向角，m是飞行器质量，σ是倾侧角，R_e是地 球的平均半径，其值为6356.766km，L和D分别为升力和阻力，g是当地的重力加速度。GER 坐标系到NED坐标系的坐标转换矩阵为

$T_{\mathrm{GER}}^{\mathrm{NED}}=\left(\begin{array}{ccc}-\cos \left(\lambda \right)\sin \left(\phi \right)& -\sin \left(\lambda \right)\sin \left(\phi \right)& \cos \left(\phi \right)\\ -\sin \left(\lambda \right)& \cos \left(\lambda \right)& 0\\ -\cos \left(\lambda \right)\cos \left(\phi \right)& -\sin \left(\lambda \right)\cos \left(\phi \right)& -\sin \left(\phi \right)\end{array}\right)---\left(7\right)$

步骤2：建立再入问题的过程约束和终端约束；

通常，飞行器再入飞行过程中需要满足热流率动压q和过载n的约束，约束条件为

$q=\frac{1}{2}\rho V^2\le q_{\max }---\left(8\right)$

$\stackrel{.}{Q}=k_Q\sqrt{\rho }{V}^{3.15}\le {\stackrel{.}{Q}}_{\max }---\left(9\right)$

$n=\frac{L}{{\mathrm{mg}}_0}{\le n}_{\max }---\left(10\right)$

对于本文所要研究的CAV模型，k_Q＝1.5×10^-8，q_max＝150kpa，n_max＝2， 其中g₀为海平面重力加速度，ρ为大气密度。

当CAV飞行到距离目标水平距离S_TAEM＝50km时，再入滑翔段结束，转入末制导段。设 计的滑翔段末端条件为：航向角误差|Δψ_TAEM|≤5deg，对旋转地球的相对速度为 V_TAEM＝2000m/s，高度为H_TAEM＝25km，倾侧角|σ_TAEM|≤30deg。其中S_TAEM，Δψ_TAEM和V_TAEM是判断制导方法优劣的主要要素。下标“TAEM”代表可重复使用飞行器(RLV)中的末端 能量管理段(TAEM)，在此，为了方便，仍如此沿用。

步骤3：引入广义赤道和广义子午线；

为了求解出解析解，先建立辅助地心惯性(AGI)坐标系，并定义广义赤道和子午线。

在惯性空间中，目标T的位置是随着旋转地球一起旋转的，所以，如附图2所示，飞行 器M飞行到预测命中点P，并不是飞行器再入飞行开始阶段的目标位置。预测命中点P可以 大致预测为

$\left(\begin{array}{c}{\lambda }_P={\lambda }_T+{\omega }_e{t}_{\mathrm{go}}\\ {\phi }_P={\phi }_T\\ H_P=H_T\end{array}\right)---\left(11\right)$

在此，λ_P、φ_P和H_P分别为碰撞点P的经纬度和高度。λ_T、φ_T和H_T分别为目标的经纬度和 高度。t_go是剩余飞行时间，可以由下式近似计算

$t_{\mathrm{go}}=\frac{2s_{\mathrm{go}}}{V+V_{\mathrm{TAEM}}}---\left(12\right)$

式中，s_go是剩余射程，假设飞行器亚轨道再入点和地心的连线与地球表面交点为M，交点M 与目标位置T之间的地球表面大圆弧长度即为剩余射程。由于单位向量的点乘即是两向量之 间夹角的余弦值，可得

$s_{\mathrm{go}}=R_e{\cos }^{-1}({\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}·{\hat{x}}_{\mathrm{ET}}^{\mathrm{GER}})---\left(13\right)$

上式中，和分别为和的单位向量，为地心E指向飞行器位置M的向量， 为地心E指向目标位置T的向量。上标“GER”代表其在GER坐标系中的分量。如图1 所示，在z_e轴的分量是sin(φ)，在赤道面的分量是cos(φ)，可得在x_e和y_e轴的分量为 cos(λ)cos(φ)和sin(λ)cos(φ)。故为

${\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}={\left(\begin{array}{ccc}\cos \left(\lambda \right)\cos \left(\phi \right)& \sin \left(\lambda \right)\cos \left(\phi \right)& \sin \left(\phi \right)\end{array}\right)}^T---\left(14\right)$

上式中，上标“T”代表矩阵或者向量的转置。类似的，可以得出为

${\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}={\left(\begin{array}{ccc}\cos \left({\lambda }_T\right)\cos \left({\phi }_T\right)& \sin \left({\lambda }_T\right)\cos \left({\phi }_T\right)& \sin \left({\phi }_T\right)\end{array}\right)}^T---\left(15\right)$

为地心E指向P的向量，的单位向量为

${\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}={\left(\begin{array}{ccc}\cos \left({\lambda }_P\right)\cos \left({\phi }_P\right)& \sin \left({\lambda }_P\right)\cos \left({\phi }_P\right)& \sin \left({\phi }_P\right)\end{array}\right)}^T---\left(16\right)$

如图2所示，定义AGI坐标系如下：坐标系原点在地心；轴指向飞行器位置M；轴 在由M、E和P组成的平面内，并垂直于轴；轴与上述两轴组成右手坐标系。AGI坐标 系中轴的单位向量为

${\hat{x}}^{\mathrm{GER}}={\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}---\left(17\right)$

由于轴与轴垂直，运用格莱姆-施密特正交化方法，可以得到AGI坐标系的轴单位向量 为

${\tilde{y}}^{\mathrm{GER}}=\frac{{\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}-({\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}·{\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}){\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}}{\left|\right|{\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}-({\hat{x}}_{\mathrm{EP}}^{\mathrm{GMR}}·{\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}){\hat{x}}_{\mathrm{EM}}^{\mathrm{GER}}\left|\right|}---\left(18\right)$

最后可以得到轴的单位向量为

${\tilde{z}}^{\mathrm{GER}}={\tilde{x}}^{\mathrm{GER}}\times {\tilde{y}}^{\mathrm{GER}}---\left(19\right)$

在此，地心E、飞行器当前位置M和预测命中点P组成的平面为MEP平面，定义MEP 平面与地球表面的交线为广义赤道平面，与广义赤道正交的大圆弧的一半定义为广义子午线， 可知广义子午线端点在轴上。通过此定义，再入飞行器(CAV)的位置可以用广义经度 广义纬度和广义高度表示。为分析方便，假设广义的0°子午线通过CAV初始再入点， 故CAV初始位置表示为

$\left(\begin{array}{c}{\tilde{\lambda }}_0=0\\ {\tilde{\phi }}_0=0\\ {\tilde{H}}_0=H\end{array}\right)---\left(20\right)$

尽管由于AGI坐标系随着CAV的运动而旋转，在此为方便分析，可认为是静止的。CAV 相对于AGI坐标系的速度，用表示，是相对地球速度和牵连速度的矢量和，初始速度为

${\tilde{v}}_0^{\mathrm{NED}}=\left(\begin{array}{c}V\cos \left(\gamma \right)\cos \left(\psi \right)\\ V\cos \left(\gamma \right)\sin \left(\psi \right)+{\omega }_e(R_e+H)\cos \left(\phi \right)\\ -V\sin \left(\gamma \right)\end{array}\right)---\left(21\right)$

在此，上标“NED”代表向量在NED坐标系中的分量形式。定义的大小为广义速度定义和当地水平面夹角为广义弹道倾角定义的水平分量于当地广义子午线之间夹角 为广义航向角由式(21)可得到，初始的和为

${\tilde{V}}_0={[V^2+2V{\omega }_e(R_e+H)\cos \left(\phi \right)\cos \left(\gamma \right)\sin \left(\psi \right)+{{\omega }_e}^2{(R_e+H)}^2{\cos }^2\left(\phi \right)]}^{\frac{1}{2}}---\left(22\right)$

${\tilde{\gamma }}_0={\sin }^{-1}(V\sin \left(\gamma \right)/{\tilde{V}}_0)---\left(23\right)$

由于亚轨道再入初始点M点沿广义子午线上的切线与平行，因此与当地水平分量 之间夹角即为初始航向角为

上式中的和是由公式(18)和(19)得到的，由如下公式得到

${\tilde{v}}_{H0}^{\mathrm{GER}}={\left(T_{\mathrm{GER}}^{\mathrm{NED}}\right)}^T{\tilde{v}}_{H0}^{\mathrm{NED}}---\left(25\right)$

${\tilde{v}}_{H0}^{\mathrm{NED}}=\left(\begin{array}{c}V\cos \left(\gamma \right)\cos \left(\psi \right)\\ V\cos \left(\gamma \right)\sin \left(\psi \right)+{\omega }_e(R_e+H)\cos \left(\phi \right)\\ 0\end{array}\right)---\left(26\right)$

上式中，是由公式(7)给出。由图2可知，P点的广义经纬度和高度为

$\left(\begin{array}{c}{\tilde{\lambda }}_P={\cos }^{-1}({\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}·{\tilde{x}}^{\mathrm{GER}})\\ {\tilde{\phi }}_p=0\\ {\tilde{H}}_p=H_T\end{array}\right)---\left(27\right)$

前面说明了对旋转地球模型下的再入制导的末端条件，这些条件需要转移到AGI坐标系 中，定义为

${\tilde{x}}_1^{\mathrm{GER}}={\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}-{\tilde{x}}^{\mathrm{GER}}---\left(28\right)$

通过运用格莱姆-施密特正交化方法，得到与正交的向量为

$x_2^{\mathrm{GER}}=x_1^{\mathrm{GER}}-(x_1^{\mathrm{GER}}·{\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}){\hat{x}}_{\mathrm{EP}}^{\mathrm{GER}}---\left(29\right)$

定义x₃为与过P点与水平面、广义赤道交线都平行的单位向量。x₃在NEDP坐标系中可表示为

$x_3^{\mathrm{NEDP}}=T_{\mathrm{GER}}^{\mathrm{NEDP}}{x}_2^{\mathrm{GER}}/\left|\right|x_2^{\mathrm{GER}}\left|\right|---\left(30\right)$

上式中，缩写词“NEDP”代表在位置P处的NED坐标系，将λ＝λ_P和φ＝φ_P带入到公式(7) 中，即得到从GER坐标系到NEDP坐标系的坐标转换矩阵

定义相对于AGI坐标系的最终速度向量为在AGI坐标系中，再入制导律导引飞行 器近似沿着广义赤道到达P，并且使得广义弹道倾角接近0，在此为方便分析，如图3所示， 假设与x₃平行。尽管这样假设不准确，会造成一定的误差，但是制导律也能够通过最后一 次倾侧角翻转来消除其引起的误差。定义相对于旋转地球的飞行器最终速度为V_TAEM。除此之 外，为得到地球旋转引起的牵连速度假设滑翔段结束距离P的距离S_TAEM远小于地球半 径。指向当地东向，大小为

$v_e^P={\omega }_e(R_e+H_{\mathrm{TAEM}})\cos \left({\phi }_P\right)---\left(31\right)$

速度向量的大小为在以上前提情况下，可得

${\tilde{v}}_f={\tilde{V}}_f{x}_3---\left(32\right)$

由图3可得

${({\tilde{v}}_f-v_e^P)}^2={\left(v_{\mathrm{TAEM}}\right)}^2---\left(33\right)$

将上式展开得

${{\tilde{V}}_f}^2-2{\tilde{V}}_f{V}_e^P\cos \left({\theta }_P\right)-2{\left(V_e^P\right)}^2=V_{\mathrm{TAEM}}^2---\left(34\right)$

$\cos \left({\theta }_P\right)=x_3^{\mathrm{NEDP}}{|}_y---\left(35\right)$

上式中，是在y轴上的分量。进而可得到的估计值为

${\tilde{V}}_f=V_e^P\cos \left({\theta }_P\right)+\sqrt{{\left(V_e^P\right)}^2({\cos }^2\left({\theta }_P\right)-1)+V_{\mathrm{TAEM}}^2}---\left(36\right)$

定义相对于惯性空间AGI坐标系的单位质量机械能为可表示为

$\tilde{E}=\frac{1}{2}{\tilde{V}}^2-\frac{\mu }{R_e+\tilde{H}}---\left(37\right)$

式中，μ是重力常数，则最终的末端约束条件转化成相对于惯性空间的能量为

${\tilde{E}}_f=\frac{1}{2}{{\tilde{V}}_f}^2-\frac{\mu }{R_e+H_{\mathrm{TAEM}}}---\left(38\right)$

步骤4：运动方程线性化；

在分析时，AGI坐标系可以认为是静止的，得到运动方程为

$\frac{d\tilde{\lambda }}{\mathrm{dt}}=\frac{\tilde{V}\cos \left(\tilde{\gamma }\right)\sin \left(\tilde{\psi }\right)}{(R_e+\tilde{H})\cos \left(\tilde{\phi }\right)}---\left(39\right)$

$\frac{d\tilde{\phi }}{\mathrm{dt}}=\frac{\tilde{V}\cos \left(\tilde{\gamma }\right)\cos \left(\tilde{\psi }\right)}{(R_e+\tilde{H})}---\left(40\right)$

$\frac{d\tilde{H}}{\mathrm{dt}}=\tilde{V}\sin \left(\tilde{\gamma }\right)---\left(41\right)$

$\frac{d\tilde{V}}{\mathrm{dt}}=-\frac{D}{m}-g\sin \left(\tilde{\gamma }\right)---\left(42\right)$

$\frac{d\tilde{\gamma }}{\mathrm{dt}}=\frac{1}{\tilde{V}}[\frac{L\cos \left(\sigma \right)}{m}-g\cos \left(\tilde{\gamma }\right)+\frac{{\tilde{V}}^2\cos \left(\tilde{\gamma }\right)}{R_e+\tilde{H}}]---\left(43\right)$

$\frac{d\tilde{\psi }}{\mathrm{dt}}=\frac{1}{\tilde{V}}[\frac{L\sin \left(\sigma \right)}{m\cos \left(\tilde{\gamma }\right)}+\frac{{\tilde{V}}^2\cos \left(\tilde{\gamma }\right)\sin \left(\tilde{\psi }\right)\tan \left(\tilde{\phi }\right)}{R_e+\tilde{H}}---\left(44\right)$

其中能量对时间的导数为

$\frac{d\tilde{E}}{\mathrm{dt}}=-\frac{D\tilde{V}}{m}---\left(45\right)$

定义射程x_D和横向射程x_C为和将公式(39，40，44)与公式(45)联立得到

$\frac{{\mathrm{dx}}_D}{d\tilde{E}}=-\frac{m\cos \left(\tilde{\gamma }\right)\sin \left(\tilde{\psi }\right)}{D\cos \left(\tilde{\phi }\right)}\frac{R_e}{(R_e+\tilde{H})}---\left(46\right)$

$\frac{{\mathrm{dx}}_C}{d\tilde{E}}=-\frac{m\cos \left(\tilde{\gamma }\right)\cos \left(\tilde{\psi }\right)}{D}\frac{R_e}{(R_e+\tilde{H})}---\left(47\right)$

$\frac{d\tilde{\psi }}{d\tilde{E}}=-\frac{\tilde{m}}{D{\tilde{V}}^2}[\frac{L\sin \left(\sigma \right)}{m\cos \left(\tilde{\gamma }\right)}+\frac{{\tilde{V}}^2\cos \left(\tilde{\gamma }\right)\sin \left(\tilde{\psi }\right)}{R_e+\tilde{H}}\tan \left(\frac{x_C}{R_e}\right)]---\left(48\right)$

令L₁＝Lcos(σ)，L₂＝Lsin(σ)。假设飞行器以平稳滑翔条件(Steady Glide Condition，SGC) 滑翔飞行，飞行器的广义弹道倾角变化率为0，即由公式(43)可得

$\frac{L_1}{m}-g\cos \left(\tilde{\gamma }\right)+\frac{{\tilde{V}}^2\cos \left(\tilde{\gamma }\right)}{R_e+\tilde{H}}=0---\left(49\right)$

经过转换可得

$\cos \left(\tilde{\gamma }\right)=\frac{L_1}{\mathrm{mg}-\frac{m{\tilde{V}}^2}{R_e+\tilde{H}}}---\left(50\right)$

由式(37)可得

${\tilde{V}}^2=2\tilde{E}+\frac{2\mu }{R_e+\tilde{H}}---\left(51\right)$

在AGI坐标系中，由于制导律导引飞行器近似沿着广义赤道飞向P，可知$\tilde{\phi }=x_C/R_e\approx 0.$因此可以假设$\cos \left(\tilde{\psi }\right)\approx \pi /2-\tilde{\psi },$$\sin \left(\tilde{\psi }\right)=1,$$\cos \left(\tilde{\phi }\right)=1,$和tan(x_C/R_e)＝x_C/R_e， 将式(50)(51)带入到使(46-48)，运用上述假设，可得

$\frac{{\mathrm{dx}}_D}{d\tilde{E}}=-\frac{L_1}{D}\frac{R_e}{-\frac{\mu }{(R_e+\tilde{H})}-2\tilde{E}}---\left(52\right)$

$\frac{{\mathrm{dx}}_C}{d\tilde{E}}=\frac{L_1}{D}\frac{R_e\tilde{\psi }}{-\frac{\mu }{R_e+\tilde{H}}-2\tilde{E}}-\frac{\pi }{2}\frac{L_1}{D}\frac{R_e}{-\frac{\mu }{R_e+\tilde{H}}-2\tilde{E}}---\left(53\right)$

$\frac{d\tilde{\psi }}{d\tilde{E}}=\frac{L_1}{D}\frac{1}{\frac{\mu }{R_e+\tilde{H}}+2\tilde{E}}\frac{x_C}{R_e}-\frac{L_2}{D}\frac{1}{2\tilde{E}+\frac{2\mu }{R_e+\tilde{H}}}---\left(54\right)$

为减小微分计算的复杂性，式(48)中的分母可假设为1。

步骤5：求解解析解；

由于对式(52-54)的解影响很小。故在分析时，可以设为一常量 令R^*＝R_e+H^*。定义纵向升阻比为L₁/D＝Lcos(σ)/D，将L₁/D写成 如下形式

$\frac{L_1}{D}=a_2{\tilde{E}}^2+a_1\tilde{E}+a_0---\left(55\right)$

定义水平升阻比为L₂/D＝Lsin(σ)/D，将L₂/D写成如下形式

$\frac{L_2}{D}=b_2{\tilde{E}}^2+b_1\tilde{E}+b_0---\left(56\right)$

在此，设计L₁/D和L₂/D曲线，可以通过同时调节攻角和倾侧角来跟踪参考升阻比曲线。

1.射程表达式

将式(55)带入式(52)，积分可得射程的表达式为

$\left(\begin{array}{c}{x}_D(\tilde{E},{\tilde{E}}_0)=\frac{1}{4}{a}_2{R}_e({\tilde{E}}^2-{\tilde{E}}_0^2)+\frac{1}{2}(a_1-\frac{a_2\mu }{{2R}^{*}})R_e\left(\widetilde{E-}{\tilde{E}}_0\right)\\ +\frac{{4a}_0{\left(R^{*}\right)}^2-{2\mathrm{\mu R}}^{*}{a}_1+{\mu }^2{a}_2}{8{\left(R^{*}\right)}^2}{R}_e\ln \left(\frac{{2R}^{*}\tilde{E}+\mu }{{2R}^{*}{\tilde{E}}_0+\mu }\right)\end{array}\right)---\left(57\right)$

2.横向射程和方位角表达式

将式(53)和(54)结合，得到

$\left(\begin{array}{c}\frac{{\mathrm{dx}}_C}{d\tilde{E}}\\ \frac{d\tilde{\psi }}{d\tilde{E}}\end{array}\right)=\frac{L_1}{D}\frac{1}{\frac{\mu }{R^{*}}+2\tilde{E}}\left(\begin{array}{cc}0& -R_e\\ \frac{1}{R_e}& 0\end{array}\right)\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)+{\left(\begin{array}{cc}\frac{\pi }{2}\frac{L_1}{D}\frac{R_e}{\frac{\mu }{R^{*}}+2\tilde{E}}& -\frac{L_2}{D}\frac{1}{2\tilde{E}+\frac{2\mu }{R^{*}}}\end{array}\right)}^T---\left(58\right)$

令

$f_1\left(\tilde{E}\right)=\frac{L_1}{D}\frac{1}{\frac{\mu }{R^{*}}+2\tilde{E}}---\left(59\right)$

$f_2\left(\tilde{E}\right)=\frac{\pi }{2}\frac{L_1}{D}\frac{R_e}{\frac{\mu }{R^{*}}+2\tilde{E}}---\left(60\right)$

$f_3\left(\tilde{E}\right)=-\frac{L_2}{D}\frac{1}{2\tilde{E}+\frac{2\mu }{R^{*}}}---\left(61\right)$

$f_4(\tilde{E},{\tilde{E}}_0)=\underset{E_0}{\overset{\tilde{E}}{\int }}{f}_1\left({\tilde{x}}_E\right)d{\tilde{x}}_E=\frac{x_D(\tilde{E},{\tilde{E}}_0)}{R_e}---\left(62\right)$

$A=\left(\begin{array}{cc}0& -R_e\\ \frac{1}{R_e}& 0\end{array}\right)---\left(63\right)$

式(58)可写成如下形式

$\left(\begin{array}{c}\frac{{\mathrm{dx}}_C}{d\tilde{E}}\\ \frac{d\tilde{\psi }}{d\tilde{E}}\end{array}\right)=f_1\left(\tilde{E}\right)A\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)+\left(\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right)---\left(64\right)$

下面我们引入谱分解方法来求解上述特殊类型的变系数线性系统。定义为

$Q(\tilde{E},{\tilde{E}}_0)=\exp (-A\underset{E_0}{\overset{\tilde{E}}{\int }}{f}_1\left({\tilde{x}}_E\right)d{\tilde{x}}_E)=\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))---\left(65\right)$

在式(64)两边乘以得到

$\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}\frac{{\mathrm{dx}}_C}{d\tilde{E}}\\ \frac{d\tilde{\psi }}{d\tilde{E}}\end{array}\right)-\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))f_1\left(\tilde{E}\right)A\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)=\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right)---\left(66\right)$

上式(66)可以写成

$\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}\frac{d{x}_C}{d\tilde{E}}\\ \frac{d\tilde{\psi }}{d\tilde{E}}\end{array}\right)+\frac{d}{d\tilde{E}}[-A{f}_4(\tilde{E},{\tilde{E}}_0)]\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)=\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right)---\left(67\right)$

为求取其微分，反向利用点乘法则，可得

$\frac{d}{d\tilde{E}}[-A{f}_4(\tilde{E},{\tilde{E}}_0)]\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)=\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right)---\left(68\right)$

对式(68)两边积分，得

$\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)-\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}{x}_{C0}\\ {\tilde{\psi }}_0\end{array}\right)=\underset{E_0}{\overset{\tilde{E}}{\int }}\exp (-A{f}_4(\tilde{E},{\tilde{E}}_0))\left(\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right)d\tilde{E}---\left(69\right)$

注意到其中0_2×2是2×2阶零矩阵，I_2×2为2×2单位阵。 的逆为

$\Phi (\tilde{E},{\tilde{E}}_0)={\left[Q(\tilde{E},{\tilde{E}}_0)\right]}^{-1}=\exp \left(A{f}_4(\tilde{E},{\tilde{E}}_0)\right)---\left(70\right)$

将式(70)左乘以式(69)得

$\left(\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right)=\Phi (\tilde{E},{\tilde{E}}_0)\left(\begin{array}{c}{x}_{C0}\\ {\psi }_0\end{array}\right)+\underset{E_0}{\overset{\tilde{E}}{\int }}\Phi (\tilde{E},{\tilde{x}}_E)\left(\begin{array}{c}{f}_2\left({\tilde{x}}_E\right)\\ f_3\left({\tilde{x}}_E\right)\end{array}\right)d{\tilde{x}}_E---\left(71\right)$

上式中，为状态转移矩阵。通过对矩阵A谱分解，由矩阵论知识可得

$\Phi (\tilde{E},{\tilde{E}}_0)=\exp \left(A{f}_4(\tilde{E},{\tilde{E}}_0)\right)=\exp \left({\lambda }_1{f}_4(\tilde{E},{\tilde{E}}_0)\right)G_1+\exp \left({\lambda }_2{f}_4(\tilde{E},{\tilde{E}}_0)\right)G_2---\left(72\right)$

上式中，λ₁＝i和λ₂＝-i是矩阵A的特征值，其中G₁和G₂是矩阵A的谱投影，为

$G_1=\frac{A-{\lambda }_{2}I}{{\lambda }_1-{\lambda }_2}=\left(\begin{array}{cc}1/2& R_{e}i/2\\ -i/\left(2R_e\right)& 1/2\end{array}\right)---\left(73\right)$

$G_2=\frac{A-{\lambda }_{1}I}{{\lambda }_2-{\lambda }_1}=\left(\begin{array}{cc}1/2& -R_{e}i/2\\ i/\left(2R_e\right)& 1/2\end{array}\right)---\left(74\right)$

进而可得

$\Phi (\tilde{E},{\tilde{E}}_0)=\left(\begin{array}{cc}\cos \left(f_4(\tilde{E},{\tilde{E}}_0)\right)& -R_e\sin \left(f_4(\tilde{E},{\tilde{E}}_0)\right)\\ \sin \left(f_4(\tilde{E},{\tilde{E}}_0)\right)/R_e& \cos \left(f_4(\tilde{E},{\tilde{E}}_0)\right)\end{array}\right)---\left(75\right)$

和

$\underset{E_0}{\overset{\tilde{E}}{\int }}\Phi (\tilde{E},{\tilde{x}}_E)\left(\begin{array}{c}{f}_2\left({\tilde{x}}_E\right)\\ f_3\left({\tilde{x}}_E\right)\end{array}\right)d{\tilde{x}}_E=\underset{E_0}{\overset{\tilde{E}}{\int }}\left(\begin{array}{c}\cos \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_2\left({\tilde{x}}_E\right)-R_e\sin \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_3\left({\tilde{x}}_E\right)\\ \sin \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_2\left({\tilde{x}}_E\right)/R_e+\cos \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_3\left({\tilde{x}}_E\right)\end{array}\right)d{\tilde{x}}_E---\left(76\right)$

由式(60)和(62)，有

$\underset{E_0}{\overset{\tilde{E}}{\int }}\cos (f_4\tilde{E},{\tilde{x}}_E)f_2\left({\tilde{x}}_E\right)d{\tilde{x}}_E=-\frac{\pi R_e}{2}\underset{E_0}{\overset{\tilde{E}}{\int }}\cos \left(f_4(\tilde{E},{\tilde{x}}_E)\right)d{f}_4(\tilde{E},{\tilde{x}}_E)=\frac{\pi R_e}{2}\sin \left(f_4(\tilde{E},{\tilde{E}}_0)\right)---\left(77\right)$

$\frac{1}{R_e}\underset{E_0}{\overset{\tilde{E}}{\int }}\sin \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_2\left({\tilde{x}}_E\right)d{\tilde{x}}_E=-\frac{\pi }{2}\underset{E_0}{\overset{\tilde{E}}{\int }}\sin \left(f_4(\tilde{E},{\tilde{x}}_E)\right)d{f}_4(E,{\tilde{x}}_E)=\frac{\pi }{2}-\frac{\pi }{2}\cos \left(f_4(\tilde{E},{\tilde{E}}_0)\right)---\left(78\right)$

这样就可以求得，横向射程和航向角的解析公式为

$\left(\begin{array}{c}{x}_C(\tilde{E},{\tilde{E}}_0)=x_{C0}\cos \left(f_4(\tilde{E},{\tilde{E}}_0)\right)-R_e{\tilde{\psi }}_0\sin \left(f_4(\tilde{E},{\tilde{E}}_0)\right)\\ +\frac{\pi R_e}{2}\sin \left(f_4(\tilde{E},{\tilde{E}}_0)\right)-R_e\underset{E_0}{\overset{\tilde{E}}{\int }}\sin \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_3\left({\tilde{x}}_E\right)d{\tilde{x}}_E\end{array}\right)---\left(79\right)$

$\left(\begin{array}{c}\tilde{\psi }(\tilde{E},{\tilde{E}}_0)=\frac{x_{C0}}{R_e}\sin \left(f_4(\tilde{E},{\tilde{E}}_0)\right)+{\tilde{\psi }}_0\cos \left(f_4(\tilde{E},{\tilde{E}}_0)\right)+\frac{\pi }{2}\\ -\frac{\pi }{2}\cos \left(f_4(\tilde{E},{\tilde{E}}_0)\right)+\underset{E_0}{\overset{\tilde{E}}{\int }}\cos \left(f_4(\tilde{E},{\tilde{x}}_E)\right)f_3\left({\tilde{x}}_E\right)d{\tilde{x}}_E\end{array}\right)---\left(80\right)$

在此，需要提到的是在求解再入制导律时，飞行器飞行的攻角为方案攻角，通过调节倾侧角 来跟踪规划的L₁/D曲线，这意味着L₂/D曲线不能任意规划，也不能像式(56)一样用有限次数 多项式表示。然而由于L₂/D的表达式只包含并没有在上述微分中体现，因此上述两 个解析解仍然适用。

步骤6：验证上述解析解可行性；

下面给出一个解析解和弹道仿真的对比实例。仿真中使用的是CAV模型，初始条件为 ${\tilde{V}}_0=7200m/s,{\tilde{H}}_0=60\mathrm{km},{\tilde{\gamma }}_0=0\mathrm{deg},{\tilde{\psi }}_0=90\mathrm{deg}.$终端条件为${\tilde{E}}_f=-5.5\times {10}^4\mathrm{kj}/\mathrm{kg}.$仿真 中考虑两种情形，在第一种情形中，CAV再入滑翔段距离较长，设计的参考曲线为

L₁/D＝2.5                 (81)

在第二种情形中，CAV再入滑翔段较短，设计的参考曲线为

L₁/D＝1.5          (83)

文章“A Closed-Form Solution to Lifting Reentry”中，Bell在假设L/D和倾侧角σ为常数 的前提下，求解出射程，横向射程和航向角的解析解形式，其方法忽略了地球曲率对航向角 产生的影响，即式(43)中二次项。在此为作对比，Bell方法的解析解分别表示为和 为了便于区别，将本文中给出的解析解用Yu-Chen法(YCGF)表示，Bell的解析解用 Bell法(BGF)表示。根据上述参考剖面，可将BGF修改为如下形式

$x_D^{\mathrm{Bell}}(\tilde{V},{\tilde{V}}_0)=L_1/D\frac{{\left(R^{*}\right)}^2}{\mu }{\int }_{\tilde{V}}^{{\tilde{V}}_0}\frac{\mu v\cos ({\tilde{\psi }}^{\mathrm{Bell}}(v,{\tilde{V}}_0)-\frac{\pi }{2})}{\mu -v^2{R}^{*}}\mathrm{dv}---\left(85\right)$

$x_D^{\mathrm{Bell}}(\tilde{V},{\tilde{V}}_0)=L_1/D\frac{{\left(R^{*}\right)}^2}{\mu }{\int }_{\tilde{V}}^{{\tilde{V}}_0}\frac{\mu v\sin ({\tilde{\psi }}^{\mathrm{Bell}}(v,{\tilde{V}}_0)-\frac{\pi }{2})}{\mu -v^2{R}^{*}}\mathrm{dv}---\left(86\right)$

其中，

上式中，通过弹道仿真可算得在时，大约为5802.29m/s。弹道仿真最终 的末速度为在第一种情形中，k_sgn＝1，在第二种情形中，k_sgn＝-1。在计 算解析解积分时采用20个点的高斯-勒让得积分法。在弹道仿真中，采用了消除长周期振动 的弹道阻尼技术。除此之外，为了能够更好的跟踪参考曲线，仿真中忽略了对倾侧角变化速 率的限制。仿真程序用C语言编写，在装有因特尔酷睿-2核处理器的计算机上运行得到下 面的结果。

仿真结果如附图4-6所示，由下表1中可以看出，YCGF具有很好的精度，而其运算速 度比弹道仿真小两个数量级。

首先，我们分析射程方向结果。由于在做微分运算时假设这就意味着忽略了 航向角误差，由图中看出，YGCF法的射程结果比弹道仿真的结果要大。而在实际运用时， 再入制导律可以通过多次倾侧角翻转来减小航向角误差，提高YCGF法射程的准确性。在这 里也给出了BGF法的仿真结果作为对比。由于BGF法在推导射程解析式时考虑了航向角误 差，BGF法在射程上的精度要比YCGF法要高，在第二种情形中比较明显。然而，相对于 YCGF法，在在线实时规划L₁/D剖面时BGF法需要更大的运算载荷，因为在线规划曲线是 需要考虑射程和航向角的交互影响，而且需要数值方法计算积分。

接下来分析横向射程和航向角的结果。从图表的结果中可以看出来BGF法在横向射程和 航向角上的精度要比YCGF方法低，尤其是在射程较大时显得更加明显。这是由于BGF法在 做微分计算时，省略了式(44)括弧中二次项。除此之外，在YCGF法中射程表达式中含有一 个积分项，由于YCGF法使得倾侧角翻转只需做几次更新，而不是实时规划，所以运算速度 更快。

表1弹道仿真与解析解的结果对比

本发明在给出上述新的高超声速飞行器再入飞行问题射程、横向射程和航向角的解析解 后，根据这些解析解本发明给出一种新的再入制导方法。首先，根据飞行器再入飞行的弹道 特点，可以将整个再入飞行过程分成三个阶段：快速下降段(Descent Phase，DP)，平稳滑翔段 (Steady Glide Phase，SGP)，弹道下压段(Altitude Adjustment Phase，AAP)。制导律分别给出了每 段的相应的制导方法。具体阐述如下：

步骤1：规划快速下降段的制导方案；

CAV在与运载器分离之后进入无动力滑翔状态，开始进入快速下降段，直到开始满足平 稳滑翔条件(SGC)时结束。在这一段，由于大气密度ρ非常小，飞行器快速下降，高度快 速降低。随着高度下降，大气密度升高，飞行器的热流率急剧升高，而最大热流率大致出现 在这段的结束。为了能满足热流率约束，设计此段飞行器飞行时保持最大攻角(α_max)并保持 倾侧角为0，这样可以使得在下降段飞行的更高，从而降低热流率。当时，说明此时开 始升力已经足够大，可以满足平稳滑翔条件。为保证飞行器能从DP平滑转移到SGP，下式 (88)给出攻角方案。指令攻角和倾侧角为

σ_cmd＝σ_plan＝0deg         (89)

Δγ＝γ_plan-γ         (90)

上式中，α_plan是SGP的规划的攻角，在式(91)中给出。γ_plan是平稳滑翔的弹道倾角，将在后 面给出说明。γ为当前飞行器的弹道倾角，Δγ₁是当时的Δγ，k_γ是值为5的常数。由式 (88)可知，Δγ从Δγ₁逐渐趋近于0的同时，指令攻角α_cmd从α_max平滑趋向于α_plan。当Δγ＝0， 时，下降段DP结束，因为之后飞行器将向上飞。在此处飞行器相对于惯性空间的单位 质量接机械能表示为

步骤2：规划平稳滑翔段的制导方案；

在平稳滑翔段，飞行器沿着规划的攻角曲线飞行。制导律利用射程方向的解析表达式实 时更新L₁/D曲线，并通过调节倾侧角跟踪纵向升阻比曲线。制导律利用横向射程的表达式来 确定倾侧角的第一次翻转节点。下面给出详细结果。

1.规划的攻角(α_plan)和升阻比(L/D)_plan

(L/D)_plan是方案攻角(α_plan)对应的升阻比曲线。L/D通常随着攻角和马赫数(Ma)变化。 Ma受大气运动影响，但是在此可以认为大气是相对旋转地球静止的。定义相对于旋转地球 的单位质量机械能E，如式(91)所示，在此，将方案攻角α_plan和(L/D)_plan表示为能量E的函 数形式。这样做的优势在于不同飞行任务的α_plan和(L/D)_plan不需要再进行调整修改。

$E=\frac{1}{2}{V}^2-\frac{\mu }{R_e+H}---\left(91\right)$

为了能够飞出更远的射程，飞行器在平稳滑翔段大部分时间内以α₁＝10deg飞行，这样 可以近似保持最大的升阻比(L/D)_max。当飞行器飞行到轨迹末段时，平滑的将α_plan减小到 α₂＝6deg，可以将高度H降低到H_TAEM进入末制导段，这样可以保证飞行器有足够的动压满 足末端机动飞行。在此设计α_plan为

${\alpha }_{\mathrm{plan}}=\left(\begin{array}{cc}{\alpha }_1;& E_2\le E\le E_0\\ {\left(\frac{E_2-E}{E_2-E_f}\right)}^2({\alpha }_2-{\alpha }_1)+{\alpha }_1;& E_f\le E<E_2\end{array}\right)---\left(92\right)$

上式中E₀是初始的能量，平稳滑翔飞行结束时的E₂可设为-5.6×10⁷J/kg，E_f是由终端约束 所决定的能量，可以通过将V_TAEM和H_TAEM带入式(91)中得到。

当α_plan曲线确定后，就确定了高度走廊，如图7所示。高度走廊的下边界(H_min)由前面 所述的过程约束决定。除此之外忽略地球旋转的影响，假设倾侧角σ＝0，，这样可以 通过求解式(5)得到高度走廊的上边界(H_max)。

下面给出相对能量E的升阻比(L/D)_plan曲线。由于α_plan曲线已经确定，(L/D)_plan是马赫 数Ma的函数，而Ma是速度V和高度H的函数。所以，如果确定了高度曲线，就可以得到相 应的(L/D)_plan曲线。实际上，当飞行器高速滑翔时，升阻比L/D很轻微的随Ma变化，但是 由于高度曲线被限定在很小的走廊中，H对Ma影响很小，所以高度对(L/D)_plan曲线影响很 小。在此，可以假设飞行器沿着高度走廊的中间飞行，高度为

$H\left(E\right)=\frac{H_{\max }\left(E\right)+H_{\min }\left(E\right)}{2}---\left(93\right)$

当给定E和H时，可以得到V和Ma，然后利用α_plan和Ma，可以得到相对于E的(L/D)_plan曲线，如图8所示。从图中可以看出，最大的升阻比可达3。

在利用解析解快速生成参考弹道的过程中，为了确定第一次翻转节点需要知道在惯 性空间中的基准升阻比曲线这就需要建立E与之间的转换关系。由式 (22,37,91)中可知，E与之间的关系为

$\tilde{E}=E+V{\omega }_e(R_e+H)\cos \left(\phi \right)\cos \left(\gamma \right)\sin \left(\psi \right)+0.5{{\omega }_e}^2{(R_e+H)}^2{\cos }^2\left(\phi \right)---\left(94\right)$

由于在实际飞行中，飞行任务剩下的弹道曲线并不明确，故上式中的关系并不实用。通常 E₀＝-3.5×10⁴kJ/kg，E_f＝-6×10⁴kJ/kg，V₀ω_e(R_e+H₀)≈3.3×10³kJ/kg， V_fω_e(R_e+H_f)≈0.93×10³kJ/kg和ω_e²(R₄+H₀)²≈0.22×10³kJ/kg。从中可以看出式(94)中的非 线性项远小于线性项，即线性项占主要作用，因此在实际仿真结果中与E看起来近似直线。 这里与E之间近似采用如下线性转化关系

$x_E=\frac{E_f-E}{{\tilde{E}}_f-\tilde{E}}({\tilde{x}}_E-\tilde{E})+E---\left(95\right)$

上式中，E与为当前相对旋转地球和相对于惯性空间的能量值，E_f与(式(37))为相对 旋转地球和相对于惯性空间的终端能量值。x_E和任意时刻相对旋转地球和相对于惯性空间 的能量值。由于(L/D)_plan曲线一大部分可认为是常值，由式(95)确立的误差较小，并且误 差可以通过第二次倾侧角翻转来修正。在此，令与相关的(L/D)_plan曲线为可以 由下式得到

在下式(97)中，当指令攻角α_plan曲线和高度曲线确立之后，可以将其转化为最大倾侧角 σ_max约束。当飞行器以平稳滑翔条件飞行时，增大σ会使得高度H降低。当将H降低到H_min时，即可得到最大的σ_max。如式(98)所示，在平稳滑翔时，纵向升力L₁与重力和向心力近似 平衡，由式(97)中右边第一项，可以求得σ_max。若飞行轨迹有振荡，式(97)右边第二项可以消 除振荡。因此，若高度H快速下降，到达(H)_E-(H_min)_E＞0，式(97)右边第二项会降低σ_max， 使得飞行器能够拉起，保持在H_min之上。在此，为书写方便，函数f对E求导，写成(f)_E。

${\sigma }_{\max }=\arccos \left(\frac{L_1}{L_{\max }}\right)+k_{\sigma }(\frac{\mathrm{dH}}{\mathrm{dE}}-\frac{{\mathrm{dH}}_{\min }}{\mathrm{dE}})---\left(97\right)$

$L_1=m(g-\frac{{{\tilde{V}}_0}^2}{R_e+H_{\min }\left(E\right)})---\left(98\right)$

L_max＝C_Lplan[0.5ρ(H_min)V²]S_ref      (99)

上式中，k_σ可设为-50。C_Lplan为α_plan的升力系数，由式(3，4，91)可知，E的变化率为

$\left(\begin{array}{c}\frac{\mathrm{dE}}{\mathrm{dt}}=V\stackrel{·}{V}+g\stackrel{·}{H}\\ =-\frac{\mathrm{DV}}{m}+V{{\omega }_e}^2(R_e+H){\cos }^2\left(\phi \right)\sin \left(\gamma \right)-V{{\omega }_e}^2(R_e+H)\sin \left(\phi \right)\cos \left(\phi \right)\cos \left(\gamma \right)\cos \left(\psi \right)\end{array}\right)---\left(100\right)$

忽略地球旋转影响，由式(3)和式(100)，可得

$\frac{\mathrm{dH}}{\mathrm{dE}}\approx -\frac{m\sin \left(\gamma \right)}{D}---\left(101\right)$

2.规划纵向升阻比L₁/D曲线

再入制导律利用射程上的解析解实时规划纵向升阻比曲线，然后通过调节倾侧角跟踪规 划的升阻比曲线，以满足终端速度约束。如图8所示，当E＞E₂时，(L/D)_plan可以认为是常 量。故时，设计纵向升阻比L₁/D曲线为水平段。当时，(L/D)_plan随着E下降而 下降，导致L₁/D不能认为是常数。但是，在此可认为(L₁/D)₁到(L₁/D)₂的转变是线性的，L₁/D 曲线设计为

$L_1/D\left({\tilde{x}}_E\right)=\left(\begin{array}{c}{(L_1/D)}_1;{\tilde{E}}_2\le {\tilde{x}}_E\le {\tilde{E}}_1\\ \frac{{\tilde{x}}_E-{\tilde{E}}_f}{{\tilde{E}}_2-{\tilde{E}}_f}{(L_1/D)}_1+\frac{{\tilde{E}}_2-{\tilde{x}}_E}{{\tilde{E}}_2-{\tilde{E}}_f},{(L_1/D)}_2;{\tilde{E}}_f\le {\tilde{x}}_E<{\tilde{E}}_2\end{array}\right)---\left(102\right)$

上式中，是将x_E＝E₂带入到式(95)中得到的。(L₁/D)₁和(L₁/D)₂是需要实时修正的参量。 由于倾侧角σ的余弦为纵向升阻比与总升阻比之比，而且最终的倾侧角约束要求接近0。 (L₁/D)₂则由下式实时计算

${(L_1/D)}_2={(L/D)}_{\mathrm{plan}}\left(E_f\right)\frac{{(L/D)}_{\mathrm{est}}\left(E\right)}{{(L/D)}_{\mathrm{plan}}\left(E\right)}---\left(103\right)$

其中，(L/D)_est(E)是由惯性导航系统(Inertial Navigation System，INS)得到的，(L/D)_est(E) 是在当前能量下的规划升阻比，(L/D)_plan(E_f)是在终点能量下的规划升阻比。在AGI坐标系 中，剩余射程为

$x_{\mathrm{Df}}=R_e{\tilde{\lambda }}_P-S_{\mathrm{TAEM}}---\left(104\right)$

上式中，由式(26)求取，下面计算(L₁/D)₁：

a.当$\tilde{E}>{\tilde{E}}_2$

将式(102)带入式(52)，然后积分，得

$x_D({\tilde{E}}_2,\tilde{E})+x_D({\tilde{E}}_f,{\tilde{E}}_2)=x_{\mathrm{Df}}---\left(105\right)$

其中当时有

$x_D({\tilde{E}}_2,\tilde{E})=\frac{{(L_1/D)}_1{R}_e}{2}\ln \left(\frac{{2R}^{*}{\tilde{E}}_2+\mu }{{2R}^{*}\tilde{E}+\mu }\right)---\left(106\right)$

当$\tilde{E}{<\tilde{E}}_2$时有

$x_D({\tilde{E}}_f,\tilde{E})=\left(\begin{array}{c}\frac{R_e[{(L_1/D)}_1-{(L_1/D)}_2]({\tilde{E}}_f-\tilde{E})}{2({\tilde{E}}_2-{\tilde{E}}_f)}-\frac{R_e{(L_1/D)}_1}{2({\tilde{E}}_2-{\tilde{E}}_f)}({\tilde{E}}_f+\frac{\mu }{{2R}^{*}})\ln \left(\frac{{2R}^{*}{\tilde{E}}_f+\mu }{{2R}^{*}\tilde{E}+\mu }\right)\\ +\frac{R_e{\left(L_{1}D\right)}_2}{2({\tilde{E}}_2-{\tilde{E}}_f)}({\tilde{E}}_2+\frac{\mu }{{2R}^{*}})\ln \left(\frac{{2R}^{*}{\tilde{E}}_f+\mu }{{2R}^{*}\tilde{E}+\mu }\right)\end{array}\right)---\left(107\right)$

其中，将带入式(107)可得通过求解式(105)，得

${(L_1/D)}_1=\frac{c_1}{c_2}---\left(108\right)$

其中

$c_1=x_{\mathrm{Df}}-\frac{1}{2}{R}_e{(L_1/D)}_2-\frac{R_e{(L_1/D)}_2}{2({\tilde{E}}_2-{\tilde{E}}_f)}({\tilde{E}}_2+\frac{\mu }{{2R}^{*}})\ln \left(\frac{{2R}^{*}{\tilde{E}}_f+\mu }{{2R}^{*}{\tilde{E}}_2+\mu }\right)---\left(109\right)$

$c_2=\frac{R_e}{2}\ln \left(\frac{{2R}^{*}{\tilde{E}}_2+\mu }{{2R}^{*}\tilde{E}+\mu }\right)-\frac{1}{2}{R}_e-\frac{1}{2}\frac{R_e}{{\tilde{E}}_2-{\tilde{E}}_f}({\tilde{E}}_f+\frac{\mu }{{2R}^{*}})\ln \left(\frac{{2R}^{*}{\tilde{E}}_f+\mu }{{2R}^{*}{\tilde{E}}_2+\mu }\right)---\left(110\right)$

b.当$\tilde{E}{<\tilde{E}}_2$

在这种情形下，由于大部分时间飞行器处于AAP阶段，不需要跟踪纵向升阻比L₁/D曲 线，因此其不需要再更新。

3.确定第一次倾侧角翻转节点

再入制导律利用横向射程的解析解表达式(式79)，来决定第一次倾侧角翻转的能量节点 倾侧角翻转是规划的，规划的横向升阻比L₂/D曲线为

上式中，是第二次倾侧角翻转的能量节点。在此假设剩余飞行中的L/D微分由INS给出， 预测的为

sgn的符号由再入初始条件决定

$\mathrm{sgn}=\left(\begin{array}{c}1;{\tilde{\psi }}_0\in [0,\pi /2)\\ -1;{\tilde{\psi }}_0\in [\pi /2,\pi ]\end{array}\right)---\left(113\right)$

从上面的结果可以看出，这里安排两次倾侧角翻转，第一次在第二次在并且 第二次倾侧角翻转离终点状态较近，这样选取接近末端的有下列好处：的调 节可以修正之前积累的误差，并且使得后面剩余较短射程内的误差较小，近似可以忽略。为 方便利用横向射程解，构造如下积分函数

$F({\tilde{x}}_{E2},{\tilde{x}}_{E1})={\int }_{{\tilde{x}}_{E1}}^{{\tilde{x}}_{E2}}\sin \left(f_4({\tilde{E}}_f,{\tilde{x}}_E)\right)f_5\left({\tilde{x}}_E\right)d{\tilde{x}}_E---\left(114\right)$

$f_5\left({\tilde{x}}_E\right)=-\frac{\left|(L_2/D)\left({\tilde{x}}_E\right)\right|}{2{\tilde{x}}_E+\frac{2\mu }{R^{*}}}---\left(115\right)$

上述积分函数与横向射程解析解积分形式类似。由于纵向升阻比L₁/D曲线是分段函数， 也可以写成如下分段函数形式

$f_4({\tilde{E}}_f,{\tilde{x}}_E)=\left(\begin{array}{c}[x_D({\tilde{E}}_f,{\tilde{E}}_2)+x_D({\tilde{E}}_2,{\tilde{x}}_E)]/R_e;{\tilde{x}}_E\ge {\tilde{E}}_2\\ x_D({\tilde{E}}_f,{\tilde{x}}_E)/R_e;{\tilde{x}}_E<{\tilde{E}}_2\end{array}\right)---\left(116\right)$

其中，和由式(107)求得，可由式(106)求得。由于假设最终的横向射程解析式仅是的函数。将初始条件式(20-23)和L₂/D曲线式(111-113)带入横向 射程解析式(79)得到最终形式为

$x_{\mathrm{Cf}}\left({\tilde{E}}_3\right)=(\frac{\pi }{2}-{\tilde{\psi }}_0)R_e\sin \left(\frac{x_{\mathrm{Df}}}{R_e}\right)-\mathrm{sgn}{R}_{e}F({\tilde{E}}_f,\tilde{E})+2\mathrm{sgn}{R}_{e}F({\tilde{E}}_4,{\tilde{E}}_3)---\left(117\right)$

其导数为

$x_{\mathrm{Cf}}^{\prime }\left({\tilde{E}}_3\right)=-2\mathrm{sgn}{R}_e\sin \left[f_4({\tilde{E}}_f,{\tilde{E}}_3)\right]f_5\left({\tilde{E}}_3\right)---\left(118\right)$

在此，运用牛顿法求解的解

${\tilde{E}}_3^{(k+1)}={\tilde{E}}_3^{\left(k\right)}-\frac{x_{\mathrm{Cf}}\left({\tilde{E}}_3^{\left(k\right)}\right)}{x_{\mathrm{Cf}}^{\prime }\left({\tilde{E}}_3^{\left(k\right)}\right)}---\left(119\right)$

为了提升运算效率，做三次运算。第一次运算是在平稳滑翔段的开始，得到当 T₁＝200s，时，第二次计算，得到当T₂＝30s， 时，第三次计算得到在这里，是因为倾侧角变化率的限 制所加的延迟，在式(121)中给出。a_D是由惯导系统所测得的当前的阻力加速度。

4.确定基准倾侧角

由于L₁＝L cos(σ)，L₁/D＝L/D cos(σ)，为跟踪纵向升阻比L₁/D曲线，基准的倾侧角为

${\sigma }_{\mathrm{plan}}=\left(\begin{array}{c}{\mathrm{sgn}·\cos }^{-1}\left(\frac{L_1/D}{{(L/D)}_{\mathrm{est}}}\right);\tilde{E}>{\tilde{E}}_3+\Delta \tilde{E}\\ {-\mathrm{sgn}·\cos }^{-1}\left(\frac{L_1/D}{{(L/D)}_{\mathrm{est}}}\right);{\tilde{E}}_4+\Delta \tilde{E}<\tilde{E}\le {\tilde{E}}_3+\Delta \tilde{E}\end{array}\right)---\left(120\right)$

其中，是考虑到飞行器的滚转速率的限制，提前时间Δt翻转基准倾侧角，其值为

$\Delta \tilde{E}=a_D{\tilde{V}}_0\mathrm{\Delta t}---\left(121\right)$

上式中，Δt完成倾侧角翻转所需时间的一半，可以用下式估计

$\mathrm{\Delta t}=\left|\frac{\sigma }{{\stackrel{·}{\sigma }}_{\max }}\right|---\left(122\right)$

其中，σ是当前倾侧角，是最大滚转速率。

5.攻角指令α_cmd和倾侧角指令σ_cmd

若将基准攻角α_plan和倾侧角σ_plan直接用作攻角和倾侧角指令，再入飞行轨迹将会有微弱 的长周期振动(Phugoid Oscillation)。在此，下面给出一种消除长周期振动的方法。如图9所 示，升力系数的垂直分量C_L1和水平分量C_L2，分别由基准攻角α_plan和倾侧角σ_plam决定。为消 除振动，在垂直方向上附加一个与垂直速度分量相反的升力，即在C_L1上加上ΔC_L1，设计ΔC_L1为

${\mathrm{\Delta C}}_{L1}=C_L^{\alpha }{k}_{\gamma }({\gamma }_{\mathrm{plan}}-\gamma )---\left(123\right)$

其中，是升力线斜率。γ_plan为平稳滑翔所需要的近似弹道倾角，在后面式(132)中给出。k_γ为增益系数，值为5。若飞行器上升较快，γ_plan-γ＜0，ΔC_L1＜0，造成升力L₁降低，从而抑 制飞行器上升过快，同样的可以抑制下降过快。

a.指令攻角α_cmd推导

对于给定的马赫数Ma，攻角是升力系数C_L的函数，α＝f_α(C_L)，其一阶Taylor展开为

${\alpha }_{\mathrm{cmd}}=f_{\alpha }(\sqrt{{(C_{L1}+{\mathrm{\Delta C}}_{L1})}^2+C_{L2}^2}\approx f_{\alpha }\left(\sqrt{{\left(C_{L1}\right)}^2+C_{L2}^2}\right)+\frac{f_{\alpha }^{\prime }\left(\sqrt{{\left(C_{L1}\right)}^2+C_{L2}^2}\right)C_{L1}}{\sqrt{{\left(C_{L1}\right)}^2+C_{L2}^2}}{\mathrm{\Delta C}}_{L1}---\left(124\right)$

注意到有

$f_{\alpha }^{\prime }\left(\sqrt{C_{L1}^2+C_{L2}^2}\right)=\frac{\mathrm{d\alpha }}{{\mathrm{dC}}_L}=\frac{1}{C_L^{\alpha }}---\left(125\right)$

将式(123，125)代入到式(124)中，得

α_cmd＝α_plan+cos(σ_plan)k_γ(γ_plan-γ)   (126)

b.指令倾侧角σ_cmd推导

由图9可知

${\sigma }_{\mathrm{cmd}}=\arctan \left(\frac{C_{L2}}{C_{L1}+{\mathrm{\Delta C}}_{L1}}\right)---\left(127\right)$

其一阶Taylor近似为

$\left(\begin{array}{c}{\sigma }_{\mathrm{cmd}}\approx \arctan \left(\frac{C_{L2}}{C_{L!}}\right)-\frac{C_{L2}}{C_{L1}^2+C_{L2}^2}\Delta C_{L1}\\ ={\sigma }_{\mathrm{plan}}-\frac{C_{L2}}{C_L}\frac{C_L^{\alpha }}{C_L}{k}_{\gamma }({\gamma }_{\mathrm{plan}}-\gamma )\\ \approx {\sigma }_{\mathrm{plan}}-\sin \left({\sigma }_{\mathrm{plan}}\right)\frac{k_{\gamma }({\gamma }_{\mathrm{plan}}-\gamma )}{{\alpha }_{\mathrm{plan}1}}\end{array}\right)---\left(128\right)$

其中，对于CAV，令α_plan＝10deg。为了保证飞行中满足各个过程约束，倾侧角有最大界限 为

其中，σ_max为最大倾侧角约束，由式(97)得到。

c.基准弹道倾角γ_plan推导

γ_plan为平稳滑翔的近似弹道倾角。由式(5)，假设忽略地球自转影响，可得

$L_{\mathrm{plan}}\cos \left({\sigma }_{\mathrm{plan}}\right)-\mathrm{mg}\cos \left({\gamma }_{\mathrm{plan}}\right)+\frac{m{V}^2\cos \left({\gamma }_{\mathrm{plan}}\right)}{R_e+H}=0---\left(130\right)$

将上式对时间求导得

$\left(\begin{array}{c}\frac{1}{m}\frac{d{C}_{\mathrm{Lplan}}}{\mathrm{dE}}\frac{\mathrm{dE}}{\mathrm{dt}}\frac{1}{2}\rho V^2{S}_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)+\frac{1}{m}{C}_{\mathrm{Lplan}}\frac{1}{2}\frac{\mathrm{d\rho }}{\mathrm{dH}}\frac{\mathrm{dH}}{\mathrm{dt}}{V}^2{S}_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)\\ +\frac{1}{m}{C}_{\mathrm{Lplan}}\mathrm{\rho V}\frac{\mathrm{dV}}{\mathrm{dt}}{S}_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)+\frac{1}{m}{L}_{\mathrm{plan}}\frac{d}{\mathrm{dE}}\left[\cos \left({\sigma }_{\mathrm{plan}}\right)\right]\frac{\mathrm{dE}}{\mathrm{dt}}\\ -\frac{\mathrm{dg}}{\mathrm{dH}}\frac{\mathrm{dH}}{\mathrm{dt}}\cos \left({\gamma }_{\mathrm{plan}}\right)+g\frac{d{\gamma }_{\mathrm{plan}}}{\mathrm{dt}}\sin \left({\gamma }_{\mathrm{plan}}\right)+\frac{1}{R_0+H}[2V\frac{\mathrm{dV}}{\mathrm{dt}}\cos \left({\gamma }_{\mathrm{plan}}\right)-V^2\sin \left({\gamma }_{\mathrm{plan}}\right)\frac{d{\gamma }_{\mathrm{plan}}}{\mathrm{dt}}]\\ -\frac{1}{{(R_0+H)}^2}{V}^2\frac{\mathrm{dH}}{\mathrm{dt}}\cos \left({\gamma }_{\mathrm{plan}}\right)=0\end{array}\right)---\left(131\right)$

其中，由式(100)可得

$\frac{\mathrm{dE}}{\mathrm{dt}}=-\frac{\mathrm{DV}}{m}---\left(132\right)$

为简化式(131)，假设γ_plan和g的导数为0。除此之外，γ_plan≈0，令cos(γ_plan)＝1， sin(γ_plan)＝γ_plan。将式(3-4)带入式(131)中，忽略地球自转影响，可以得到

${\gamma }_{\mathrm{plan}}=-\frac{D_{\mathrm{plan}}}{\mathrm{mg}}\frac{d_1}{d_2}---\left(133\right)$

上式中，D_plan＝C_DplanqS_ref，C_Dplan是攻角为α_plan时的阻力系数，ρ是大气密度，S_ref为参考面 积，d₁和d₂分别表示如下

$d_1=\frac{\rho V^2{S}_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)}{2m}\frac{d{C}_{\mathrm{Lplan}}}{\mathrm{dE}}+\frac{2}{R_0+H}+\frac{C_{\mathrm{Lplan}}\rho S_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)}{m}+\frac{L_{\mathrm{plan}}}{m}\frac{d}{\mathrm{dE}}\left[\cos \left({\sigma }_{\mathrm{plan}}\right)\right]---\left(134\right)$

$d_2=-\frac{C_{\mathrm{Lplan}}{V}^2{S}_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)}{2\mathrm{mg}}\frac{\mathrm{d\rho }}{\mathrm{dH}}+\frac{2}{R_0+H}+\frac{C_{\mathrm{Lplan}}\rho S_{\mathrm{ref}}\cos \left({\sigma }_{\mathrm{plan}}\right)}{m}+\frac{V^2}{{(R_0+H)}^{2}g}---\left(135\right)$

前面提到，将函数f对E的导数定为(f)_E。因为C_Lplan(E)已经规划，(C_Lplan)_E是已知的。在 平稳滑翔段可用式(136)求得[cos(σ_plan)]_E。然而，在弹道下压段，由于[cos(σ_plan)]_E的计算 公式较复杂，可用相邻两时刻的cos(σ_plan)的差分计算。

$\frac{d}{\mathrm{dE}}\left[\cos \left({\sigma }_{\mathrm{plan}}\right)\right]=\frac{d}{\mathrm{dE}}\left[\frac{L_1/D}{{(L/D)}_{\mathrm{est}}}\right]=\frac{1}{{(L/D)}_{\mathrm{est}}}\frac{d{L}_1/D}{\mathrm{dE}}-\frac{L_1/D}{{(L/D)}_{\mathrm{est}}{(L/D)}_{\mathrm{plan}}}\frac{d{(L/D)}_{\mathrm{plan}}}{\mathrm{dE}}---\left(136\right)$

其中，(L/D)_plan(E)是规划曲线，((L/D)_plan)_E是已知的，(L₁/D)_E可由式(95)和(102)得到。

步骤3：规划弹道下压段的制导方案；

飞行器经过第二次倾侧角翻转之后，进入弹道下压段(AAP)。通过降低攻角，飞行器快 速进入稠密大气层，这样能够使得飞行器在末制导段有较大的动压进行机动。这就使得在这 段中，不能近似接近0，所以使得上述解析解难以满足终端约束条件。在此用一种与前面不 同的制导方案：在最后一次倾侧角翻转之前，为得到需求的最终速度，在线的弹道仿真修正 能量节点在最后一次倾侧角翻转之后，采用比例导引法(Proportional Navigation，PN) 导引，下面首先介绍弹道下压段的制导方案，再给出修正的方法。

1.基准倾侧角

在弹道下压段，α_plan曲线与平稳滑翔段相同，而倾侧角由比例导引律确定。当飞行器接 近目标时，可以忽略地球曲率的影响。如图2所示，MP视线的方位角的变化率为

${\stackrel{·}{\psi }}_{\mathrm{MP}}=\frac{{\tilde{V}}_0\cos \left({\tilde{\gamma }}_0\right)\sin (\pi /2-{\tilde{\psi }}_0)}{x_{\mathrm{DP}}}---\left(137\right)$

其中，由比例导引律得到的需用横向机动加速度为

$a_{L2}=k_{\mathrm{PN}}{\stackrel{·}{\psi }}_{\mathrm{MP}}{\tilde{V}}_0\cos \left({\tilde{\gamma }}_0\right)---\left(138\right)$

其中，为了防止初始倾侧角饱和，令k_PN从2逐渐变化到4，为

$k_{\mathrm{PN}}=2\frac{x_{\mathrm{Df}}}{x_{\mathrm{Df}}^{E4}}+4(1-\frac{x_{\mathrm{Df}}}{x_{\mathrm{Df}}^{E4}})---\left(139\right)$

其中，x_Df是第二次倾侧角翻转的能量节点。另外，升力的垂直分量需要与重力和向心力平衡， 可得

$a_{L1}\approx g-\frac{{\tilde{V}}_0^2}{R_e+H}---\left(140\right)$

进而可以得到基准倾侧角为

${\sigma }_{\mathrm{plan}}=\mathrm{sgn}·\arctan \left(\frac{a_{L2}}{a_{L1}}\right),\tilde{E}\le {\tilde{E}}_4+\Delta \tilde{E}---\left(141\right)$

其中，由式(121)求得。

2.指令攻角和指令倾侧角

在飞行器进入弹道下压段之前，需要获得射程-能量参考曲线并且需要修正在弹道下压段，通过调节攻角跟踪曲线，如式(142)所示，可以消除由于干扰引起的最 终速度误差。注意到在弹道下压段中，大部分时间升阻比L/D对攻角导数大于0，意味着增 加攻角会增加升阻比，造成最终速度升高。为了消除长周期振动，σ_cmd与式(128)一样。

${\alpha }_{\mathrm{cmd}}={\alpha }_{\mathrm{plan}}+k_{\alpha }[s_{\mathrm{go}}\left(E\right)-s_{\mathrm{go}}^{\mathrm{ref}}\left(E\right)]---\left(142\right)$

${\sigma }_{\mathrm{cmd}}={\sigma }_{\mathrm{plan}}-\sin \left({\sigma }_{\mathrm{plan}}\right)\frac{k_{\gamma }({\gamma }_{\mathrm{plan}}-\gamma )}{{\alpha }_{\mathrm{plan}1}}---\left(143\right)$

其中，k_α值为(5π/18)×10^-6，也就是说0.5deg的攻角修正，会引起10km的射程散布。在弹 道下压段中，在计算式(133)中的γ_plan时，由于[cos(σ_plan)]_E表达式较复杂，其对E的导数由 差分求得。同时，为了满足过程约束，仍然需要最大倾侧角约束

3.第二次倾侧角翻转

在弹道下压段，不需要跟踪L₁/D曲线，需要调整翻转能量节点来满足最终速度要求。 下面分析与最终速度V_f之间的关系：若降低则翻转时间推迟，时刻的航向角误差 增加，为消除航向角误差，需要更大的横向机动加速度，因而倾侧角的大小会增大，造成纵 向升阻比L₁/D减小，导致最终速度降低。所以，降低造成最终速度V_f降低。

在第一次倾侧角翻转之后，在此采用预测校正法修正通过弹道仿真预测最终速度， 之后运用secant法修正当前的状态作为仿真的初始状态，由于INS能够实时测量当前气 动加速度，通过对比实际和理想气动加速度，可以得到当前气动系数的偏差。弹道预测仿真 假设当前估计的偏差即为剩余飞行段的气动系数偏差。除此之外，仿真中忽略式(142)中右边 的第二项。当S_go＝S_TAEM时，仿真结束，可以得到预测的最终速度修正的secant 法式为

${\tilde{E}}_4^{(k+1)}={\tilde{E}}_4^{\left(k\right)}-\frac{(V_f\left({\tilde{E}}_4^{\left(k\right)}\right)-V_{\mathrm{TAEM}})({\tilde{E}}_4^{\left(k\right)}-{\tilde{E}}_4^{(k-1)})}{(V_f\left({\tilde{E}}_4^{\left(k\right)}\right)-V_f\left({\tilde{E}}_4^{(k-1)}\right))}---\left(145\right)$

其中，为第k次仿真的为减小运算次数，的修正有两次：第一次预测校正过程在 第一次倾侧角翻转之后，得到T₁＝60s时，此时进行第二次预 测校正过程，得到之后可令

当CAV飞行到S_go＝S_TAEM时，仿真结束。

3、优点及功效：本发明的优点在于：

(1)本发明给出新的关于能量的射程、横向射程和航向角的解析解，求取解析解时考虑 了地球曲率的影响，精度更高。

(2)相较于传统的跟踪参考轨迹和预测校正制导，本发明通过新的解析解规划再入制导 方法，能很好的避免对整个飞行过程中的重复积分，降低运算载荷，也便于工程实现。

(3)本发明将再入飞行过程分为三段，分别规划每段的再入制导方法。为了减小飞行控 制系统(FCS)的负荷，再入飞行设计成只有两侧倾侧角翻转，相比与传统方法的多次翻转， 有很大的提升。

(4)由于典型的再入飞行轨迹有一定的长周期振动，本发明给出三维的抑制长周期振动 的方法，使得飞行轨迹更加平稳。

(5)由于参考剖面和倾侧角翻转是在线规划的，本发明给出的制导方法能够很好的在一 系列气动和大气参数散布情况下，保持很好的鲁棒性。

附图说明

图1是GER坐标系：E-x_ey_ez_e和NED坐标系o-xyz的示意图

图2是辅助地心惯性(AGI)坐标系的示意图

图3是在AF坐标系下的终点速度近似计算示意图

图4是解析解和弹道仿真结果中的射程对比结果图

图5是解析解和弹道仿真结果中的横向射程对比结果图

图6是解析解和弹道仿真结果中的航向角对比结果图

图7是规划的高度走廊图

图8是基准升阻比曲线图

图9是附加升力系数抑制弹道振动的示意图

图10是制导方案流程框图

图11是三种理想情形仿真结果的高度-速度曲线图

图12是三种理想情形仿真结果的地面投影曲线图

图13是三种理想情形仿真结果的弹道倾角曲线图

图14是三种理想情形仿真结果的攻角曲线图

图15是三种理想情形仿真结果的倾侧角曲线图

图16是气动拉偏情形下的升阻比拉偏百分比曲线图

图17是气动拉偏情形下仿真结果的高度-速度曲线图

图18是气动拉偏情形下仿真结果的地面投影曲线图

图19是气动拉偏情形下仿真结果的攻角曲线图

图20是气动拉偏情形下仿真结果的倾侧角曲线图

图21是YCEG法Monto Carlo仿真的高度-速度曲线图

图22是LGG法Monto Carlo仿真的高度-速度曲线图

图23是YCEG法Monto Carlo仿真的地面投影曲线图

图24是LGG法Monto Carlo仿真的地面投影曲线图

图25是YCEG法Monto Carlo仿真的攻角曲线图

图26是LGG法Monto Carlo仿真的攻角曲线图

图27是YCEG法Monto Carlo仿真的倾侧角曲线图

图28是LGG法Monto Carlo仿真的倾侧角曲线图

图29是YCEG法Monto Carlo仿真的热流率曲线图

图30是LGG法Monto Carlo仿真的热流率曲线图

图31是YCEG法Monto Carlo仿真的最终状态图

图32是LGG法Monto Carlo仿真的最终状态图

图中符号说明如下：

图1中：x_e、y_e和z_e轴分别为地心赤道旋转(GER)坐标系的三个坐标轴，x、y和z轴 分别为当地北东下(NED)坐标系的三个坐标轴，E为地心，M是飞行器质心，o为飞行器 在地面投影，λ是经度，φ是纬度，H是高度，V是对旋转地球的相对速度，V_H为V的水平 分量，γ是飞行的弹道倾角，ψ是航向角，ω_e为地球自转角速度。

图2中：E为地心，M是飞行器质心，T为飞行器再入时目标位置，P为预测命中点。 和轴分别为辅助地心惯性(AGI)坐标系三个坐标轴。

图3中：P为预测命中点，x_P和y_P轴分别为P点NED坐标系的坐标轴，为最终速度 向量，V_TAEM为相对于旋转地球模型的飞行器最终速度，为地球旋转引起的牵连速度，θ_P为与y_P之间夹角，x₃为与过P点与水平面和广义赤道交线平行的单位向量。

图4-6中：横坐标的为相对惯性空间的单位质量机械能。

图7-8中：横坐标的E为相对旋转地球的单位质量机械能，图7中的H_min为高度走廊的 下边界，H_max为高度走廊的上边界，图8中的L/D_plan为参考的升阻比。

图9中：M为飞行器质心，C_L1和C_L2分别为升力系数的垂直分量和水平分量。ΔC_L1为 抑制长周期振动附加的升力系数垂直分量。

图13中：γ为弹道倾角，γ_plan为平稳滑翔所需要的近似弹道倾角。

图11、17、21和22中的H_min为高度走廊的下边界。

具体实施方式

下面具体结合附图以及仿真实例对本发明做进一步详细的说明。

图1是GER坐标系：E-x_ey_ez_e和NED坐标系o-xyz的示意图，形象显示飞行器再入模 型的相关位置和角度关系。图2是辅助地心惯性(AGI)坐标系的示意图，显示AGI 系的定义。图3是在AF坐标系下的终点速度近似计算示意图，帮助求解旋转地球模型下的 相关速度关系。图4-6分别是解析解和弹道仿真结果中的纵向射程对比、横向射程以及航向 角对比结果图，可以很清楚显示本文解析解的可行性。图7是规划的高度走廊图，是由相关 的过程约束所确定的。图8是基准升阻比曲线图，对于不同的飞行任务，基准升阻比曲线近 似不变。图9是附加升力系数抑制弹道振动的示意图。为形象说明本发明的流程，图10给出 了本发明的制导方案流程框图。

下面具体结合几个仿真案例来展示该制导方法的可行性及优势。

本发明用于再入制导仿真的模型是CAV，其最大升阻比可达3。由于飞行控制系统的限 制，约束攻角的变化率为：和倾侧角约束为：|σ|≤80deg、和飞行器的初始条件为：H₀＝80km，λ₀＝0rad，φ₀＝0rad，V₀＝6500m/s， γ₀＝0rad，

案例1：无干扰理想情形

首先观察三种理想的情形，不考虑干扰的影响，目标之间距离较远。情形1，目标位置 为λ_T1R_e＝7500km，φ_T1R_e＝1000km；情形2，λ_T2R_e＝10000km，φ_T2R_e＝0；情形3， λ_T3R_e＝12000km，φ_T3R_e＝2000km。仿真结果如图11-15所示，较小射程情形下，为跟踪纵 向升阻比L₁/D曲线，会降低倾侧角σ的大小，从而降低高度H。图11是三种理想情形仿真 结果的高度-速度曲线图，可以看出高度速度曲线始终保持在下边界之上，能满足相关的过程 约束。图12是三种理想情形仿真结果的地面投影曲线图，显示飞行器能够准确飞行到距离目 标一定距离上空，并近似指向目标，航向角误差很小。图13是三种理想情形仿真结果的弹道 倾角曲线图，显示γ_plan是对γ的很好的近似。图14是三种理想情形仿真结果的攻角曲线图， 为抑制长周期振动，攻角并没有保持规划的值，而是有很小的反馈修正。图15是三种理想情 形仿真结果的倾侧角曲线图。仿真结果还显示，消除长周期振动的方法能够很好的消除由于 倾侧角翻转造成的轨迹振荡。

表2三种理想情形的仿真结果

案例2：气动拉偏情形

在实际飞行中，存在气动参数散布。一般情况下，气动拉偏系数是攻角和马赫数的函数。 为简化考虑，下面给出线性的气动力拉偏系数模型

${\delta }_{\mathrm{CL}}={\delta }_{\mathrm{CL}0}+k_{\mathrm{\delta CL}}^{\mathrm{Ma}}\frac{\mathrm{Ma}-15}{17}+k_{\mathrm{\delta CL}}^{\alpha }\frac{18}{\pi }(\alpha -\frac{\pi }{18})---\left(146\right)$

${\delta }_{\mathrm{CD}}={\delta }_{\mathrm{CD}0}+k_{\mathrm{\delta CD}}^{\mathrm{Ma}}\frac{\mathrm{Ma}-15}{17}+k_{\mathrm{\delta CD}}^{\alpha }\frac{18}{\pi }(\alpha -\frac{\pi }{18})---\left(147\right)$

上式中，δ_CL是升力系数C_L拉偏系数，δ_CD是阻力系数C_D拉偏系数。由于L/D对射程有 重要影响，下面给出升阻比L/D较大偏差两种情形。这里仅考虑阻力系数拉偏的情况，即 δ_CL＝0。情形1，令δ_CD0＝-20％，和情形2，令δ_CD0＝25％，和目标位置为λ_T1R_e＝9000km，φ_T1R_e＝1000km。气动拉偏结果如附图16所示， 仿真结果如附图17-20所示。图17是气动拉偏情形下仿真结果的高度-速度曲线图，高度速 度曲线始终保持在下边界之上。图18是气动拉偏情形下仿真结果的地面投影曲线图，显示飞 行器能够准确飞行到距离目标一定距离上空，并近似指向目标，航向角误差很小。图19是气 动拉偏情形下仿真结果的攻角曲线图。图20是气动拉偏情形下仿真结果的倾侧角曲线图。仿 真结果显示若C_D增加，导致L/D降低，倾侧角σ大小降低，从而降低横向机动能力，为满足 横向射程要求，第一次翻转时间推迟。经过第二次倾侧角翻转之后，为跟踪剩余射程-能量曲 线，需要轻微的调整α。

表3气动拉偏结果

案例3：Monte Carle仿真情形

表4正态分布的随机变量数值特性

Monte Carlo仿真可以检验再入制导律的鲁棒性。如表4所示，给出了再入初始条件，气 动模型和大气模型的散布。式(146-147)给出了气动拉偏模型。δ_ρ代表大气密度的百分比拉偏， 代表东西向的风速，代表南北向的风速。在每次仿真中，δ_ρ、和为常数。 尽管散布模型较为简单，并不与实际情形相符合，但是这些情形比实际情况更为严格，能够 很好的验证制导律的鲁棒性。作为对比，同时给出文献“Constrained Predictor-Corrector Entry  Guidance”中的制导律Monte Carlo仿真结果。为区分两种方法，本文中的制导方法称为 Yu-Chen's Entry Guidance(YCEG)，文献中的制导方法称为Lu's Gliding Guidance(LGG)。为满 足热流率约束，对LGG做轻微修改，在飞行器快速下降段，攻角最大，倾侧角为0。除此之 外，LGG的攻角曲线也是通过式(92)给出的，由于LGG最终的倾侧角趋近于30deg，为达到 H_TAEM，将α₂设为7deg。多次仿真绘制仿真结果如图21-32所示。图21和22分别是YCEG 法和LGG法Monto Carlo仿真的高度-速度曲线图。图23和24分别是YCEG法和LGG法 Monto Carlo仿真的地面投影曲线图，可以看出两种方法都能使得飞行器准确飞行到距离目标 一定距离上空，并近似指向目标，航向角误差较小。图25和26分别是YCEG法和LGG法 Monto Carlo仿真的攻角曲线图，YCEG法中攻角的轻微的反馈修正能够很好的为抑制飞行轨 迹的长周期振动。图27和28分别是YCEG法和LGG法Monto Carlo仿真的倾侧角曲线图。 从仿真结果可以看出，YCEG能够很有效的消除由倾侧角翻转造成的振荡。在LGG中，由于 消除振荡的增益随时间下降，在飞行末端，振荡不能够很好的抑制。在YCEG中，整个飞行 过程只有两次倾侧角翻转，而且位置比较集中。但是在LGG中，需要由3-5次倾侧角翻转， 而且由于L/D的散布，次数和位置不确定，更大的L/D会增加倾侧角翻转的次数。图29和 30分别为YCEG法和LGG法Monto Carlo仿真的热流率曲线图，热流率都没有超过约束值。 图31和32分别是YCEG法和LGG法Monto Carlo仿真的最终状态图。从图31-32中，可以 看出来，YCEG的最终航向角误差远小于LGG。尽管有几个仿真结果，可能是由于大气密度 ρ的下降，造成出现H稍小于H_min的情况，但是过程约束仍然能满足。