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

摘要

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

主权项

基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法，实施该制导方法前需要进行求解解析，具体步骤如下：步骤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{d\lambda }{dt}=\frac{Vcos\left(\gamma \right)sin\left(\psi \right)}{(R_e+H)cos\left(\phi \right)}---\left(1\right)$$\frac{d\phi }{dt}=\frac{Vcos\left(\gamma \right)cos\left(\psi \right)}{(R_e+H)}---\left(2\right)$$\frac{dH}{dt}=Vsin\left(\gamma \right)---\left(3\right)$$\begin{array}{c}\frac{dV}{dt}=-\frac{D}{m}-g\sin \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}---\left(4\right)$$\begin{array}{c}\frac{d\gamma }{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}+2{\mathrm{V\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}---\left(5\right)$$\begin{array}{c}\frac{d\psi }{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}+2{\mathrm{V\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)}-2{\mathrm{V\omega }}_e\cos \left(\phi \right)\tan \left(\gamma \right)\cos \left(\psi \right)]\end{array}---\left(6\right)$在上述运动方程中，t是时间，λ是经度，φ是纬度，H是高度，V是对旋转地球的相对速度，γ是飞行的弹道倾角，Ψ是航向角，m是飞行器质量，σ是倾侧角，R_e是地球的平均半径，其值为6356.766km，L和D分别为升力和阻力，g是当地的重力加速度，GER坐标系到NED坐标系的坐标转换矩阵为$T_{GER}^{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}{\mathrm{\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的位置是随着旋转地球一起旋转的，飞行器M飞行到预测命中点P，并不是飞行器再入飞行开始阶段的目标位置，预测命中点P大致预测为$\left\{\begin{array}{c}{\lambda }_P={\lambda }_T+{\omega }_e{t}_{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_{go}=\frac{2s_{go}}{V+V_{TAEM}}---\left(12\right)$式中，s_go是剩余射程，假设飞行器亚轨道再入点和地心的连线与地球表面交点为M，交点M与目标位置T之间的地球表面大圆弧长度即为剩余射程；由于单位向量的点乘即是两向量之间夹角的余弦值，得$s_{go}=R_e{\cos }^{-1}({\hat{x}}_{EM}^{GER}·{\hat{x}}_{ET}^{GER})---\left(13\right)$上式中，和分别为和的单位向量，为地心E指向飞行器位置M的向量，为地心E指向目标位置T的向量；上标“GER”代表其在GER坐标系中的分量，在z_e轴的分量是sin(φ)，在赤道面的分量是cos(φ)，可得在x_e和y_e轴的分量为cos(λ)cos(φ)和sin(λ)cos(φ)；故为${\hat{x}}_{EM}^{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}}_{ET}^{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}}_{EP}^{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)$定义AGI坐标系如下：坐标系原点在地心；轴指向飞行器位置M；轴在由M、E和P组成的平面内，并垂直于轴；轴与上述两轴组成右手坐标；AGI坐标系中轴的单位向量为${\tilde{x}}^{GER}={\hat{x}}_{EM}^{GER}---\left(17\right)$由于轴与轴垂直，运用格莱姆‑施密特正交化方法，得到AGI坐标系的轴单位向量为${\tilde{y}}^{GER}=\frac{{\hat{x}}_{EP}^{GER}-({\hat{x}}_{EP}^{GER}·{\hat{x}}_{EM}^{GER}){\hat{x}}_{EM}^{GER}}{||{\hat{x}}_{EP}^{GER}-({\hat{x}}_{EP}^{GER}·{\hat{x}}_{EM}^{GER}){\hat{x}}_{EM}^{GER}||}---\left(18\right)$最后可以得到轴的单位向量为${\tilde{z}}^{GER}={\tilde{x}}^{GER}\times {\tilde{y}}^{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^{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+2{\mathrm{V\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}\left(Vsin\right(\gamma )/{\tilde{V}}_0)---\left(23\right)$由于亚轨道再入初始点M点沿广义子午线上的切线与平行，因此与当地水平分量之间夹角即为初始航向角为${\tilde{\psi }}_0=\left\{\begin{array}{cc}{\cos }^{-1}({\tilde{z}}^{GER}·{\tilde{V}}_{H0}^{GER}/||{\tilde{V}}_{H0}^{GER}||);& {\tilde{y}}^{GER}·{\tilde{V}}_{H0}^{GER}\ge 0\\ -{\cos }^{-1}({\tilde{z}}^{GER}·{\tilde{V}}_{H0}^{GER}/||{\tilde{V}}_{H0}^{GER}||);& {\tilde{y}}^{GER}·{\tilde{V}}_{H0}^{GER}<0\end{array}\right.---\left(24\right)$上式中的和是由公式(18)和(19)得到的，由如下公式得到${\tilde{V}}_{H0}^{GER}={\left(T_{GER}^{NED}\right)}^T{\tilde{V}}_{H0}^{NED}---\left(25\right)$${\tilde{V}}_{H0}^{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)给出P点的广义经纬度和高度为$\left\{\begin{array}{c}{\tilde{\lambda }}_P={\cos }^{-1}({\hat{x}}_{EP}^{GER}·{\tilde{x}}^{GER})\\ {\tilde{\phi }}_P=0\\ {\tilde{H}}_P=H_T\end{array}\right.---\left(27\right)$前面说明了对旋转地球模型下的再入制导的末端条件，这些条件需要转移到AGI坐标系中，定义为${\tilde{x}}_1^{GER}={\hat{x}}_{EP}^{GER}-{\tilde{x}}^{GER}---\left(28\right)$通过运用格莱姆‑施密特正交化方法，得到与正交的向量为$x_2^{GER}=x_1^{GER}-(x_1^{GER}·{\hat{x}}_{EP}^{GER}){\hat{x}}_{EP}^{GER}---\left(29\right)$定义x₃为与过P点与水平面、广义赤道交线都平行的单位向量，x₃在NEDP坐标系中表示为$x_3^{NEDP}=T_{GER}^{NEDP}{x}_2^{GER}/||x_2^{GER}||---\left(30\right)$上式中，缩写词“NEDP”代表在位置P处的NED坐标系，将λ＝λ_P和φ＝φ_P带入到公式(7)中，即得到从GER坐标系到NEDP坐标系的坐标转换矩阵定义相对于AGI坐标系的最终速度向量为在AGI坐标系中，再入制导方法导引飞行器近似沿着广义赤道到达P，并且使得广义弹道倾角接近0，在此为方便分析，假设与x₃平行，制导律也能够通过最后一次倾侧角翻转来消除其引起的误差；定义相对于旋转地球的飞行器最终速度为V_TAEM，除此之外，为得到地球旋转引起的牵连速度假设滑翔段结束距离P的距离S_TAEM远小于地球半径，指向当地东向，大小为$V_e^P={\omega }_e(R_e+H_{TAEM})cos\left({\phi }_P\right)---\left(31\right)$速度向量的大小为在以上前提情况下，可得${\tilde{V}}_f={\tilde{V}}_f{x}_3---\left(32\right)$${({\tilde{V}}_f-V_e^P)}^2={\left(V_{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_{TAEM}^2---\left(34\right)$$cos\left({\theta }_P\right)=x_3^{NEDP}{|}_y---\left(35\right)$上式中，是在y轴上的分量，进而得到的估计值为${\tilde{V}}_f=V_e^{P}cos\left({\theta }_P\right)+\sqrt{{\left(V_e^P\right)}^2\left({\cos }^2\right({\theta }_P)-1)+V_{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_{TAEM}}---\left(38\right)$步骤4：运动方程线性化；在分析时，AGI坐标系认为是静止的，得到运动方程为$\frac{d\tilde{\lambda }}{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 }}{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}}{dt}=\tilde{V}sin\left(\tilde{\gamma }\right)---\left(41\right)$$\frac{d\tilde{V}}{dt}=-\frac{D}{m}-g\sin \left(\tilde{\gamma }\right)---\left(42\right)$$\frac{d\tilde{\gamma }}{dt}=\frac{1}{\tilde{V}}[\frac{Lcos\left(\sigma \right)}{m}-gcos\left(\tilde{\gamma }\right)+\frac{{\tilde{V}}^{2}cos\left(\tilde{\gamma }\right)}{R_e+\tilde{H}}]---\left(43\right)$$\frac{d\tilde{\psi }}{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}}{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{mcos\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{mcos\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}-gcos\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}{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，可知因此假设和tan(x_C/R_e)＝x_C/R_e，将式(50)、(51)带入到式(46)、(47)和(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)(53)和(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曲线，通过同时调节攻角和倾侧角来跟踪参考升阻比曲线；①射程表达式将式(55)带入式(52)，积分可得射程的表达式为$\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(\tilde{E}-{\tilde{E}}_0)\\ +\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}---\left(57\right)$②横向射程和方位角表达式将式(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{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}{f}_1({\tilde{x}}_E\left)d{\tilde{x}}_E\right)=\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)---\left(65\right)$在式(64)两边乘以得到$\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\left[\begin{array}{c}\frac{{\mathrm{dx}}_C}{d\tilde{E}}\\ \frac{d\tilde{\psi }}{d\tilde{E}}\end{array}\right]-\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)f_1\left(\tilde{E}\right)A\left[\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right]=\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\left[\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right]---\left(66\right)$上式(66)写成$\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\left[\begin{array}{c}\frac{{\mathrm{dx}}_C}{d\tilde{E}}\\ \frac{d\tilde{\psi }}{d\tilde{E}}\end{array}\right]+\frac{d}{d\tilde{E}}[-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0)]\left[\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right]=\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\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}}[-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0)]\left[\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right]=\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\left[\begin{array}{c}{f}_2\left(\tilde{E}\right)\\ f_3\left(\tilde{E}\right)\end{array}\right]---\left(68\right)$对式(68)两边积分，得$\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\left[\begin{array}{c}{x}_C\\ \tilde{\psi }\end{array}\right]-\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\left[\begin{array}{c}{x}_{C0}\\ {\tilde{\psi }}_0\end{array}\right]=\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\exp (-{\mathrm{Af}}_4(\tilde{E},{\tilde{E}}_0\left)\right)\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({\mathrm{Af}}_4\right(\tilde{E},{\tilde{E}}_0\left)\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{{\tilde{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({\mathrm{Af}}_4\right(\tilde{E},{\tilde{E}}_0\left)\right)=\exp \left({\lambda }_1{f}_4\right(\tilde{E},{\tilde{E}}_0\left)\right)G_1+\exp \left({\lambda }_2{f}_4\right(\tilde{E},{\tilde{E}}_0\left)\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\right(\tilde{E},{\tilde{E}}_0\left)\right)& -R_e\sin \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)\\ \sin \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)/R_e& \cos \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)\end{array}\right]---\left(75\right)$和$\underset{{\tilde{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{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\left[\begin{array}{c}\cos \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_2\left({\tilde{x}}_E\right)-R_e\sin \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_3\left({\tilde{x}}_E\right)\\ \sin \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_2\left({\tilde{x}}_E\right)/R_e+\cos \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_3\left({\tilde{x}}_E\right)\end{array}\right]d{\tilde{x}}_E---\left(76\right)$由式(60)和(62)，有$\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\cos \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_2\left({\tilde{x}}_E\right)d{\tilde{x}}_E=-\frac{{\mathrm{\pi R}}_e}{2}\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\cos \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right){\mathrm{df}}_4(\tilde{E},{\tilde{x}}_E)=\frac{{\mathrm{\pi R}}_e}{2}\sin \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)---\left(77\right)$$\frac{1}{R_e}\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\sin \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_2\left({\tilde{x}}_E\right)d{\tilde{x}}_E=-\frac{\pi }{2}\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\sin \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)d{f}_4(E,{\tilde{x}}_E)=\frac{\pi }{2}-\frac{\pi }{2}\cos \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)---\left(78\right)$这样就求得，横向射程和航向角的解析公式为$\begin{array}{c}{x}_C(\tilde{E},{\tilde{E}}_0)=x_{C0}\cos \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)-R_e{\tilde{\psi }}_0\sin \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)\\ +\frac{{\mathrm{\pi R}}_e}{2}\sin \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)-R_e\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\sin \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_3\left({\tilde{x}}_E\right)d{\tilde{x}}_E\end{array}---\left(79\right)$$\begin{array}{c}\tilde{\psi }(\tilde{E},{\tilde{E}}_0)=\frac{x_{C0}}{R_e}\sin \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right){\tilde{\psi }}_0\cos \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)+\frac{\pi }{2}\\ -\frac{\pi }{2}\cos \left(f_4\right(\tilde{E},{\tilde{E}}_0\left)\right)+\underset{{\tilde{E}}_0}{\overset{\tilde{E}}{\int }}\cos \left(f_4\right(\tilde{E},{\tilde{x}}_E\left)\right)f_3\left({\tilde{x}}_E\right)d{\tilde{x}}_E\end{array}---\left(80\right)$在此，需要提到的是在求解再入制导方法时，飞行器飞行的攻角为方案攻角，通过调节倾侧角来跟踪规划的L₁/D曲线，这意味着L₂/D曲线不能任意规划，也不能像式(56)一样用有限次数多项式表示；然而由于L₂/D的表达式只包含并没有在上述微分中体现，因此上述两个解析解仍然适用；步骤6：验证上述解析解可行性；该基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法，具体步骤如下：步骤A：规划快速下降段的制导方案；CAV在与运载器分离之后进入无动力滑翔状态，开始进入快速下降段，直到开始满足平稳滑翔条件即SGC时结束；在这一段，由于大气密度ρ非常小，飞行器快速下降，高度快速降低；随着高度下降，大气密度升高，飞行器的热流率急剧升高，而最大热流率大致出现在这段的结束；为了能满足热流率约束，设计此段飞行器飞行时保持最大攻角(α_max)并保持倾侧角为0，这样使得在下降段飞行的更高，从而降低热流率；当时，说明此时开始升力已经足够大，满足平稳滑翔条件；为保证飞行器能从DP平滑转移到SGP，下式(81)给出攻角方案，指令攻角和倾侧角为${\alpha }_{cmd}=\left\{\begin{array}{cc}{\alpha }_{max},& \stackrel{·}{\gamma }<0\mathrm{deg}/s\\ \frac{\Delta \gamma }{{\mathrm{\Delta \gamma }}_1}{\alpha }_{max}+\frac{{\mathrm{\Delta \gamma }}_1-\Delta \gamma }{{\mathrm{\Delta \gamma }}_1}[{\alpha }_{plan}+k_{\gamma }\Delta \gamma ],& \stackrel{·}{\gamma }\ge 0\mathrm{deg}/s\end{array}\right.---\left(81\right)$σ_cmd＝σ_plan＝0deg                   (82)Δγ＝γ_plan‑γ                  (83)上式中，α_plan是SGP的规划的攻角，在式(84)中给出；γ_plan是平稳滑翔的弹道倾角，γ为当前飞行器的弹道倾角，Δγ₁是当时的Δγ，k_γ是值为5的常数；由式(81)可知，Δγ从Δγ₁逐渐趋近于0的同时，指令攻角α_cmd从α_max平滑趋向于α_plan；当Δγ＝0，时，下降段DP结束，因为之后飞行器将向上飞，在此处飞行器相对于惯性空间的单位质量接机械能表示为步骤B：规划平稳滑翔段的制导方案；在平稳滑翔段，飞行器沿着规划的攻角曲线飞行；制导律利用射程方向的解析表达式实时更新L₁/D曲线，并通过调节倾侧角跟踪纵向升阻比曲线，制导律利用横向射程的表达式来确定倾侧角的第一次翻转节点，下面给出详细结果；(1)规划的攻角(α_plan)和升阻比(L/D)_plan(L/D)_plan是方案攻角(α_plan)对应的升阻比曲线，(L/D)通常随着攻角和马赫数(Ma)变化，Ma受大气运动影响，在此认为大气是相对旋转地球静止的；定义相对于旋转地球的单位质量机械能E，如式(84)所示，在此，将方案攻角α_plan和(L/D)_plan表示为能量E的函数形式；这样做的优势在于不同飞行任务的α_plan和(L/D)_plan不需要再进行调整修改；$E=\frac{1}{2}{V}^2-\frac{\mu }{R_e+H}---\left(84\right)$为了能够飞出更远的射程，飞行器在平稳滑翔段大部分时间内以α₁＝10deg飞行，这样近似保持最大的升阻比(L/D)_max；当飞行器飞行到轨迹末段时，平滑的将α_plan减小到α₂＝6deg，将高度H降低到H_TAEM进入末制导段，这样保证飞行器有足够的动压满足末端机动飞行，在此设计α_plan为${\alpha }_{plan}=\{\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}---\left(85\right)$上式中E₀是初始的能量，平稳滑翔飞行结束时的E₂可设为‑5.6×10⁷J/kg，E_f是由终端约束所决定的能量，通过将V_TAEM和H_TAEM带入式(84)中得到；当α_plan曲线确定后，就确定了高度走廊，高度走廊的下边界(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(86\right)$当给定E和H时，得到V和Ma，然后利用α_plan和Ma，得到相对于E的(L/D)_plan曲线，最大的升阻比可达3；在利用解析解快速生成参考弹道的过程中，为了确定第一次翻转节点需要知道在惯性空间中的基准升阻比曲线这就需要建立E与之间的转换关系，由式22、37和91可知，E与之间的关系为$\tilde{E}=E+{\mathrm{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(87\right)$由于在实际飞行中，飞行任务剩下的弹道曲线并不明确，故上式中的关系并不实用，通常$\left\{\begin{array}{c}{E}_0=-3.5\times {10}^{4}kJ/kg\\ E_f=-6\times {10}^{4}kJ/kg\\ V_0{\omega }_e(R_e+H_0)\approx 3.3\times {10}^{3}kJ/kg\\ V_f{\omega }_e(R_e+H_f)\approx 0.93\times {10}^{3}kJ/kg\\ {{\omega }_e}^2{(R_e+H_0)}^2\approx 0.22\times {10}^{3}kJ/kg\end{array}\right.---\left(88\right)$从中看出式(87)中的非线性项远小于线性项，即线性项占主要作用，因此在实际仿真结果中与E看起来近似直线；这里与E之间近似采用如下线性转化关系$x_E=\frac{E_f-E}{{\tilde{E}}_f-\tilde{E}}({\tilde{x}}_E-\tilde{E})+E---\left(89\right)$上式中，E与为当前相对旋转地球和相对于惯性空间的能量值，E_f与即式(37)为相对旋转地球和相对于惯性空间的终端能量值；x_E和任意时刻相对旋转地球和相对于惯性空间的能量值，由于(L/D)_plan曲线一大部分是常值，由式(89)确立的误差较小，并且误差通过第二次倾侧角翻转来修正，在此，令与相关的(L/D)_plan曲线为由下式得到在下式(91)中，当指令攻角α_plan曲线和高度曲线确立之后，将其转化为最大倾侧角σ_max约束；当飞行器以平稳滑翔条件飞行时，增大σ会使得高度H降低，当将H降低到H_min时，即得到最大的σ_max；如式(92)所示，在平稳滑翔时，纵向升力L₁与重力和向心力近似平衡，由式(91)中右边第一项，求得σ_max；若飞行轨迹有振荡，式(91)右边第二项可以消除振荡，因此，若高度H快速下降，到达(H)_E‑(H_min)_E>0，式(91)右边第二项会降低σ_max，使得飞行器能够拉起，保持在H_min之上；在此，函数f对E求导，写成f_E；${\sigma }_{max}=arccos\left(\frac{L_1}{L_{\max }}\right)+k_{\sigma }(\frac{dH}{dE}-\frac{{\mathrm{dH}}_{\min }}{dE})---\left(91\right)$$L_1=m(g-\frac{{\tilde{V}}_0^2}{R_e+H_{\min }\left(E\right)})---\left(92\right)$L_max＝C_Lplan[0.5ρ(H_min)V²]S_ref            (93)上式中，k_σ设为‑50，C_Lplan为α_plan的升力系数，由式(3)(4)和(84)可知，E的变化率为$\begin{array}{c}\frac{dE}{dt}=V\stackrel{·}{V}+g\stackrel{·}{H}\\ =-\frac{DV}{m}+{{\mathrm{V\omega }}_e}^2(R_e+H){\cos }^2\left(\phi \right)\sin \left(\gamma \right)-{{\mathrm{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}---\left(94\right)$忽略地球旋转影响，由式(3)和式(94)，得$\frac{dH}{dE}\approx -\frac{msin\left(\gamma \right)}{D}---\left(95\right)$(2)规划纵向升阻比L₁/D曲线再入制导方法利用射程上的解析解实时规划纵向升阻比曲线，然后通过调节倾侧角跟踪规划的升阻比曲线，以满足终端速度约束；当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)=\{\begin{array}{cc}{(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}---\left(96\right)$上式中，是将x_E＝E₂带入到式(89)中得到的，(L₁/D)₁和(L₁/D)₂是需要实时修正的参量，由于倾侧角σ的余弦为纵向升阻比与总升阻比之比，而且最终的倾侧角约束要求接近0，(L₁/D)₂则由下式实时计算${(L_1/D)}_2={(L/D)}_{plan}\left(E_f\right)\frac{{(L/D)}_{est}\left(E\right)}{{(L/D)}_{plan}\left(E\right)}---\left(97\right)$其中，(L/D)_est(E)是由惯性导航系统即Inertial Navigation System,INS得到的，(L/D)_est(E)是在当前能量下的规划升阻比，(L/D)_plan(E_f)是在终点能量下的规划升阻比；在AGI坐标系中，剩余射程为$x_{Df}=R_e{\tilde{\lambda }}_P-S_{TAEM}---\left(98\right)$上式中，由式(26)求取，下面计算(L₁/D)₁：a.当将式(96)带入式(52)，然后积分，得$x_D({\tilde{E}}_2,\tilde{E})+x_D({\tilde{E}}_f,{\tilde{E}}_2)=x_{Df}---\left(99\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(100\right)$当时有$\begin{array}{c}{x}_D({\tilde{E}}_f,\tilde{E})=\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{(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}+\mu }\right)\end{array}---\left(101\right)$其中，将带入式(101)得通过求解式(99)，得${(L_1/D)}_1=\frac{c_1}{c_2}---\left(102\right)$其中$c_1=x_{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(103\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(104\right)$b.当在这种情形下，由于大部分时间飞行器处于AAP阶段，不需要跟踪纵向升阻比L₁/D曲线，因此其不需要再更新；(3)确定第一次倾侧角翻转节点再入制导方法利用横向射程的解析解表达式(79)，来决定第一次倾侧角翻转的能量节点倾侧角翻转是规划的，规划的横向升阻比L₂/D曲线为$L_2/D\left({\tilde{x}}_E\right)\{\begin{array}{cc}\mathrm{sgn}·\left|(L_2/D)\left({\tilde{x}}_E\right)\right|& {\tilde{x}}_E>{\tilde{E}}_{3}or{\tilde{E}}_f\le {\tilde{x}}_E\le {\tilde{E}}_4\\ -\mathrm{sgn}·\left|(L_2/D)\left({\tilde{x}}_E\right)\right|& {\tilde{E}}_4<{\tilde{x}}_E\le {\tilde{E}}_3\end{array}---\left(105\right)$上式中，是第二次倾侧角翻转的能量节点，在此假设剩余飞行中的L/D微分由INS给出，预测的为sgn的符号由再入初始条件决定$\mathrm{sgn}=\{\begin{array}{cc}1;& {\tilde{\psi }}_0\in [0,\pi /2)\\ -1;& {\tilde{\psi }}_0\in [\pi /2,\pi ]\end{array}---\left(107\right)$从上面的结果看出，这里安排两次倾侧角翻转，第一次在第二次在并且第二次倾侧角翻转离终点状态较近，这样选取接近末端的有下列好处：的调节修正之前积累的误差，并且使得后面剩余较短射程内的误差较小，近似可以忽略；为方便利用横向射程解，构造如下积分函数$F({\tilde{x}}_{E2},{\tilde{x}}_{E1})={\int }_{{\tilde{x}}_{E1}}^{{\tilde{x}}_{E2}}sin\left(f_4\right({\tilde{E}}_f,{\tilde{x}}_E\left)\right)f_5\left({\tilde{x}}_E\right)d{\tilde{x}}_E---\left(108\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(109\right)$上述积分函数与横向射程解析解积分形式类似，由于纵向升阻比L₁/D曲线是分段函数，写成如下分段函数形式$f_4({\tilde{E}}_f,{\tilde{x}}_E)=\{\begin{array}{cc}[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}---\left(110\right)$其中，和由式(101)求得，由式(100)求得；由于假设最终的横向射程解析式仅是的函数，将初始条件式(20)、(21)、(22)、(23)和L₂/D曲线式(105)、(106)、(107)带入横向射程解析式(79)得到最终形式为$x_{Cf}\left({\tilde{E}}_3\right)=(\frac{\pi }{2}-{\tilde{\psi }}_0)R_{e}sin\left(\frac{x_{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(111\right)$其导数为$x_{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(112\right)$在此，运用牛顿法求解的解${\tilde{E}}_3^{(k+1)}={\tilde{E}}_3^{\left(k\right)}-\frac{x_{Cf}\left({\tilde{E}}_3^{\left(k\right)}\right)}{x_{Cf}^{\prime }\left({\tilde{E}}_3^{\left(k\right)}\right)}---\left(113\right)$为了提升运算效率，做三次运算；第一次运算是在平稳滑翔段的开始，得到当T₁＝200s，时，第二次计算，得到当T₂＝30s，时，第三次计算得到在这里，是因为倾侧角变化率的限制所加的延迟，在式(115)中给出，a_D是由惯导系统所测得的当前的阻力加速度；(4)确定基准倾侧角由于L₁＝Lcosσ，L₁/D＝L/Dcosσ，为跟踪纵向升阻比L₁/D曲线，基准的倾侧角为${\sigma }_{plan}=\left\{\begin{array}{cc}\mathrm{sgn}·{\cos }^{-1}\left(\frac{L_1/D}{{(L/D)}_{est}}\right);& \tilde{E}>{\tilde{E}}_3+\Delta \tilde{E}\\ -\mathrm{sgn}·{\cos }^{-1}\left(\frac{L_1/D}{{(L/D)}_{est}}\right);& {\tilde{E}}_4+\Delta \tilde{E}<\tilde{E}\le {\tilde{E}}_3+\Delta \tilde{E}\end{array}\right.---\left(114\right)$其中，是考虑到飞行器的滚转速率的限制，提前时间Δt翻转基准倾侧角，其值为$\Delta \tilde{E}=a_D{\tilde{V}}_0\Delta t---\left(115\right)$上式中，Δt完成倾侧角翻转所需时间的一半，用下式估计$\Delta t=\left|\frac{\sigma }{{\stackrel{·}{\sigma }}_{max}}\right|---\left(116\right)$其中，σ是当前倾侧角，是最大滚转速率；(5)攻角指令α_cmd和倾侧角指令σ_cmd若将基准攻角α_plan和倾侧角σ_plan直接用作攻角和倾侧角指令，再入飞行轨迹将会有微弱的长周期振动，在此，下面给出一种消除长周期振动的方法；升力系数的垂直分量C_L1和水平分量C_L2，分别由基准攻角α_plan和倾侧角σ_plan决定；为消除振动，在垂直方向上附加一个与垂直速度分量相反的升力，即在C_L1上加上ΔC_L1，设计ΔC_L1为${\mathrm{\Delta C}}_{L1}=C_L^{\alpha }{k}_{\gamma }({\gamma }_{plan}-\gamma )---\left(117\right)$其中，是升力线斜率，γ_plan为平稳滑翔所需要的近似弹道倾角，在后面式(126)中给出；k_γ为增益系数，值为5，若飞行器上升较快，γ_plan‑γ<0，ΔCL₁<0，造成升力L₁降低，从而抑制飞行器上升过快，同样的可以抑制下降过快；c.指令攻角α_cmd推导对于给定的马赫数Ma，攻角是升力系数C_L的函数，α＝f_α(C_L)，其一阶Taylor展开为${\alpha }_{cmd}=f_{\alpha }(\sqrt{{(C_{L1}+{\mathrm{\Delta C}}_{L1})}^2+C_{L2}^2}\approx f_{\alpha }(\sqrt{{\left(C_{L1}\right)}^2+C_{L2}^2})+\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(118\right)$注意到有$f_{\alpha }^{\prime }\left(\sqrt{C_{L1}^2+C_{L2}^2}\right)=\frac{d\alpha }{{\mathrm{dC}}_L}=\frac{1}{C_L^{\alpha }}---\left(119\right)$将式(117)和(119)代入到式(118)中，得α_cmd＝α_plan+cos(σ_plan)k_γ(γ_plan‑γ)               (120)d.指令倾侧角σ_cmd推导${\sigma }_{cmd}=arctan\left(\frac{C_{L2}}{C_{L1}+{\mathrm{\Delta C}}_{L1}}\right)---\left(121\right)$其一阶Taylor近似为$\begin{array}{c}{\sigma }_{cmd}\approx \arctan \left(\frac{C_{L2}}{C_{L1}}\right)-\frac{C_{L2}}{C_{L1}^2+C_{L2}^2}{\mathrm{\Delta C}}_{L1}\\ ={\sigma }_{plan}-\frac{C_{L2}}{C_L}\frac{C_L^{\alpha }}{C_L}{k}_{\gamma }({\gamma }_{plan}-\gamma )\\ \approx {\sigma }_{plan}-\sin \left({\sigma }_{plan}\right)\frac{k_{\gamma }({\gamma }_{plan}-\gamma )}{{\alpha }_{plan1}}\end{array}---\left(122\right)$其中，对于CAV，令α_plan＝10deg；为了保证飞行中满足各个过程约束，倾侧角有最大界限为${\sigma }_{cmd}=\{\begin{array}{cc}{\sigma }_{\max },& {\sigma }_{cmd}>{\sigma }_{\max }\\ -{\sigma }_{\max },& {\sigma }_{cmd}<-{\sigma }_{\max }\end{array}---\left(123\right)$其中，σ_max为最大倾侧角约束，由式(91)得到；e.基准弹道倾角γ_plan推导γ_plan为平稳滑翔的近似弹道倾角，由式(5)，假设忽略地球自转影响，得$L_{plan}cos\left({\sigma }_{plan}\right)-mgcos\left({\gamma }_{plan}\right)+\frac{{\mathrm{mV}}^{2}cos\left({\gamma }_{plan}\right)}{R_e+H}=0---\left(124\right)$将上式对时间求导得$\begin{array}{c}\frac{1}{m}\frac{{\mathrm{dC}}_{Lplan}}{dE}\frac{dE}{dt}\frac{1}{2}{\mathrm{\rho V}}^2{S}_{ref}\cos \left({\sigma }_{plan}\right)+\frac{1}{m}{C}_{Lplan}\frac{1}{2}\frac{d\rho }{dH}\frac{dH}{dt}{V}^2{S}_{ref}\cos \left({\sigma }_{plan}\right)\\ +\frac{1}{m}{C}_{Lplan}\rho V\frac{dV}{dt}{S}_{ref}\cos \left({\sigma }_{plan}\right)+\frac{1}{m}{L}_{plan}\frac{d}{dE}\left[\cos \left({\sigma }_{plan}\right)\right]\frac{dE}{dt}\\ -\frac{dg}{dH}\frac{dH}{dt}\cos \left({\gamma }_{plan}\right)+g\frac{{\mathrm{d\gamma }}_{plan}}{dt}\sin \left({\gamma }_{plan}\right)+\frac{1}{R_0+H}[2V\frac{dV}{dt}\cos \left({\gamma }_{plan}\right)-V^2\sin \left({\gamma }_{plan}\right)\frac{{\mathrm{d\gamma }}_{plan}}{dt}]\\ -\frac{1}{{(R_0+H)}^2}{V}^2\frac{dH}{dt}\cos \left({\gamma }_{plan}\right)=0\end{array}---\left(125\right)$其中，由式(94)得$\frac{dE}{dt}=-\frac{DV}{m}---\left(126\right)$为简化式(125)，假设γ_plan和g的导数为0，除此之外，γ_plan≈0，令cos(γ_plan)＝1，sin(γ_plan)＝γ_plan，将式(3)和(4)带入式(125)中，忽略地球自转影响，得到${\gamma }_{plan}=-\frac{D_{plan}}{mg}\frac{d_1}{d_2}---\left(127\right)$上式中，D_plan＝C_DplanqS_ref，C_Dplan是攻角为α_plan时的阻力系数，ρ是大气密度，S_ref为参考面积，d₁和d₂分别表示如下$d_1=\frac{{\mathrm{\rho V}}^2{S}_{ref}\cos \left({\sigma }_{plan}\right)}{2m}\frac{{\mathrm{dC}}_{Lplan}}{dE}+\frac{2}{R_0+H}+\frac{C_{Lplan}{\mathrm{\rho S}}_{ref}\cos \left({\sigma }_{plan}\right)}{m}+\frac{L_{plan}}{m}\frac{d}{dE}\left[\cos \left({\sigma }_{plan}\right)\right]---\left(128\right)$$d_2=-\frac{C_{Lplan}{V}^2{S}_{ref}\cos \left({\sigma }_{plan}\right)}{2mg}\frac{d\rho }{dH}+\frac{2}{R_0+H}+\frac{C_{Lplan}{\mathrm{\rho S}}_{ref}\cos \left({\sigma }_{plan}\right)}{m}+\frac{V^2}{{(R_0+H)}^{2}g}---\left(129\right)$前面提到，将函数f对E的导数定为f_E；因为C_Lplan(E)已经规划，(C_Lplan)_E是已知的；在平稳滑翔段用式(130)求得[cos(σ_plan)]_E；然而，在弹道下压段，由于[cos(σ_plan)]_E的计算公式较复杂，用相邻两时刻的cos(σ_plan)的差分计算；$\frac{d}{dE}\left[cos\left({\sigma }_{plan}\right)\right]=\frac{d}{dE}\left[\frac{L_1/D}{{(L/D)}_{est}}\right]=\frac{1}{{(L/D)}_{est}}\frac{{\mathrm{dL}}_1/D}{dE}-\frac{L_1/D}{{(L/D)}_{est}{(L/D)}_{plan}}\frac{d{(L/D)}_{plan}}{dE}---\left(130\right)$其中，(L/D)_plan(E)是规划曲线，[(L/D)_plan]_E是已知的，(L₁/D)_E可由式(89)和(96)得到；步骤C：规划弹道下压段的制导方案；飞行器经过第二次倾侧角翻转之后，进入弹道下压段即AAP；通过降低攻角，飞行器快速进入稠密大气层，这样能够使得飞行器在末制导段有较大的动压进行机动；这就使得在这段中，不能近似接近0，所以使得上述解析解难以满足终端约束条件；在此用一种与前面不同的制导方案：在最后一次倾侧角翻转之前，为得到需求的最终速度，在线的弹道仿真修正能量节点在最后一次倾侧角翻转之后，采用比例导引法即Proportional Navigation,PN导引，下面首先介绍弹道下压段的制导方案，再给出修正的方法；①基准倾侧角在弹道下压段，α_plan曲线与平稳滑翔段相同，而倾侧角由比例导引律确定；当飞行器接近目标时，忽略地球曲率的影响，MP视线的方位角的变化率为${\stackrel{·}{\psi }}_{MP}=\frac{{\tilde{V}}_{0}cos\left({\tilde{\gamma }}_0\right)sin(\pi /2-{\tilde{\psi }}_0)}{x_{DP}}---\left(131\right)$其中，由比例导引律得到的需用横向机动加速度为$a_{L2}=k_{PN}{\stackrel{·}{\psi }}_{MP}{\tilde{V}}_{0}cos\left({\tilde{\gamma }}_0\right)---\left(132\right)$其中，为了防止初始倾侧角饱和，令k_PN从2逐渐变化到4，为$k_{PN}=2\frac{x_{Df}}{x_{Df}^{E4}}+4(1-\frac{x_{Df}}{x_{Df}^{E4}})---\left(133\right)$其中，x_Df是第二次倾侧角翻转的能量节点；另外，升力的垂直分量需要与重力和向心力平衡，可得$a_{L1}\approx g-\frac{{\tilde{V}}_0^2}{R_e+H}---\left(134\right)$进而可以得到基准倾侧角为$\begin{array}{cc}{\sigma }_{plan}=\mathrm{sgn}·arctan\left(\frac{a_{L2}}{a_{L1}}\right)& \tilde{E}\le {\tilde{E}}_4+\Delta \tilde{E}\end{array}---\left(135\right)$其中，由式(115)求得；②指令攻角和指令倾侧角在飞行器进入弹道下压段之前，需要获得射程‑能量参考曲线并且需要修正在弹道下压段，通过调节攻角跟踪曲线，如式(136)所示，消除由于干扰引起的最终速度误差；注意到在弹道下压段中，大部分时间升阻比L/D对攻角导数大于0，意味着增加攻角会增加升阻比，造成最终速度升高，为了消除长周期振动，σ_cmd与式(122)一样；${\alpha }_{cmd}={\alpha }_{plan}+k_{\alpha }[s_{go}\left(E\right)-s_{go}^{ref}\left(E\right)]---\left(136\right)$${\sigma }_{cmd}={\sigma }_{plan}-\sin \left({\sigma }_{plan}\right)\frac{k_{\gamma }({\gamma }_{plan}-\gamma )}{{\alpha }_{plan1}}---\left(137\right)$其中，k_α值为(5π/18)×10^‑6，也就是说0.5deg的攻角修正，会引起10km的射程散布；在弹道下压段中，在计算式(127)中的γ_plan时，由于[cos(σ_plan)]_E表达式较复杂，其对E的导数由差分求得；同时，为了满足过程约束，仍然需要最大倾侧角约束${\sigma }_{cmd}=\{\begin{array}{cc}{\sigma }_{\max },& {\sigma }_{cmd}>{\sigma }_{\max }\\ -{\sigma }_{\max },& {\sigma }_{cmd}<-{\sigma }_{\max }\end{array}---\left(138\right)$③第二次倾侧角翻转在弹道下压段，不需要跟踪L₁/D曲线，需要调整翻转能量节点来满足最终速度要求；下面分析与最终速度V_f之间的关系：若降低则翻转时间推迟，时刻的航向角误差增加，为消除航向角误差，需要更大的横向机动加速度，因而倾侧角的大小会增大，造成纵向升阻比L₁/D减小，导致最终速度降低；所以，降低造成最终速度V_f降低；在第一次倾侧角翻转之后，在此采用预测校正法修正通过弹道仿真预测最终速度，之后运用secant法修正当前的状态作为仿真的初始状态，由于INS能够实时测量当前气动加速度，通过对比实际和理想气动加速度，得到当前气动系数的偏差；弹道预测仿真假设当前估计的偏差即为剩余飞行段的气动系数偏差，除此之外，仿真中忽略式(136)中右边的第二项，当S_go＝S_TAEM时，仿真结束，得到预测的最终速度修正的secant法式为${\tilde{E}}_4^{(k+1)}={\tilde{E}}_4^{\left(k\right)}-\frac{\left(V_f\right({\tilde{E}}_4^{\left(k\right)})-V_{TAEM})({\tilde{E}}_4^{\left(k\right)}-{\tilde{E}}_4^{(k-1)})}{\left(V_f\right({\tilde{E}}_4^{\left(k\right)})-V_f({\tilde{E}}_4^{(k-1)}\left)\right)}---\left(139\right)$其中，为第k次仿真的为减小运算次数，的修正有两次：第一次预测校正过程在第一次倾侧角翻转之后，得到T₁＝60s时，此时进行第二次预测校正过程，得到之后可令当CAV飞行到S_go＝S_TAEM时，仿真结束。