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

摘要

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

主权项

基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法，实施该制导方法前需要进行求解解析，具体步骤如下：步骤1：建立飞行器再入问题的运动方程；首先定义两个坐标系，地心赤道旋转即GER坐标系：坐标系原点在地心，ze轴指向北极；xe和ye轴是赤道平面内相互垂直的两条轴线；GER坐标系以ze为中心，以与地球自转角速度相同的角速度旋转，地球自转角速度ωe为7.292116×10‑5rad/s；当地北东下即NED坐标系：坐标系原点o为飞行器质心M和地心连线与地球表面的交点，即飞行器在地面的垂直投影点；x轴指向当地的北向，y轴指向当地的东向，z轴铅垂向下；由于大气认为是相对地球静止的，选择在GER坐标系中处理再入问题较合适；将再入飞行器看作质点，在球形地球模型下，飞行器运动方程为dλdt=Vcos(γ)sin(ψ)(Re+H)cos(φ)---(1)]]>dφdt=Vcos(γ)cos(ψ)(Re+H)---(2)]]>dHdt=Vsin(γ)---(3)]]>dVdt=-Dm-gsin(γ)+ωe2(Re+H)cos2(φ)sin(γ)-ωe2(Re+H)sin(φ)cos(φ)cos(γ)cos(ψ)---(4)]]>dγdt=1V[Lcos(σ)m-gcos(γ)+V2cos(γ)Re+H+2Vωecos(φ)sin(ψ)+ωe2(Re+H)cos2(φ)cos(γ)+ωe2(Re+H)sin(φ)cos(φ)sin(γ)cos(ψ)]---(5)]]>dψdt=1V[Lsin(σ)mcos(γ)+V2cos(γ)sin(ψ)tan(φ)Re+H+2Vωesin(φ)ωe2(Re+H)sin(φ)cos(φ)sin(ψ)cos(γ)-2Vωecos(φ)tan(γ)cos(ψ)]---(6)]]>在上述运动方程中，t是时间，λ是经度，φ是纬度，H是高度，V是对旋转地球的相对速度，γ是飞行的弹道倾角，Ψ是航向角，m是飞行器质量，σ是倾侧角，Re是地球的平均半径，其值为6356.766km，L和D分别为升力和阻力，g是当地的重力加速度，GER坐标系到NED坐标系的坐标转换矩阵为TGERNED=-cos(λ)sin(φ)-sin(λ)sin(φ)cos(φ)-sin(λ)cos(λ)0-cos(λ)cos(φ)-sin(λ)cos(φ)-sin(φ)---(7)]]>步骤2：建立再入问题的过程约束和终端约束；通常，飞行器再入飞行过程中需要满足热流率动压q和过载n的约束，约束条件为q=12ρV2≤qmax---(8)]]>Q·=kQρV3.15≤Q·max---(9)]]>n=Lmg0≤nmax---(10)]]>对于CAV模型，kQ＝1.5×10‑8，qmax＝150kpa，nmax＝2，其中g0为海平面重力加速度，ρ为大气密度；当CAV飞行到距离目标水平距离STAEM＝50km时，再入滑翔段结束，转入末制导段；设计的滑翔段末端条件为：航向角误差|ΔψTAEM|≤5deg，对旋转地球的相对速度为VTAEM＝2000m/s，高度为HTAEM＝25km，倾侧角|σTAEM|≤30deg；其中，STAEM，ΔΨTAEM和VTAEM是判断制导方法优劣的主要要素；下标“TAEM”代表可重复使用飞行器RLV中的末端能量管理段TAEM；步骤3：引入广义赤道和广义子午线；为了求解出解析解，先建立辅助地心惯性即AGI坐标系，并定义广义赤道和子午线；在惯性空间中，目标T的位置是随着旋转地球一起旋转的，飞行器M飞行到预测命中点P，并不是飞行器再入飞行开始阶段的目标位置，预测命中点P大致预测为λP=λT+ωetgoφP=φTHP=HT---(11)]]>在此，λP、φP和HP分别为碰撞点P的经纬度和高度；λT、φT和HT分别为目标的经纬度和高度，tgo是剩余飞行时间，由下式近似计算tgo=2sgoV+VTAEM---(12)]]>式中，sgo是剩余射程，假设飞行器亚轨道再入点和地心的连线与地球表面交点为M，交点M与目标位置T之间的地球表面大圆弧长度即为剩余射程；由于单位向量的点乘即是两向量之间夹角的余弦值，得sgo=Recos-1(x^EMGER·x^ETGER)---(13)]]>上式中，和分别为和的单位向量，为地心E指向飞行器位置M的向量，为地心E指向目标位置T的向量；上标“GER”代表其在GER坐标系中的分量，在ze轴的分量是sin(φ)，在赤道面的分量是cos(φ)，可得在xe和ye轴的分量为cos(λ)cos(φ)和sin(λ)cos(φ)；故为x^EMGER=cos(λ)cos(φ)sin(λ)cos(φ)sin(φ)T---(14)]]>上式中，上标“T”代表矩阵或者向量的转置，类似的，可以得出为x^ETGER=cos(λT)cos(φT)sin(λT)cos(φT)sin(φT)T---(15)]]>为地心E指向P的向量，的单位向量为x^EPGER=cos(λp)cos(φP)sin(λp)cos(φP)sin(φP)T---(16)]]>定义AGI坐标系如下：坐标系原点在地心；轴指向飞行器位置M；轴在由M、E和P组成的平面内，并垂直于轴；轴与上述两轴组成右手坐标；AGI坐标系中轴的单位向量为x~GER=x^EMGER---(17)]]>由于轴与轴垂直，运用格莱姆‑施密特正交化方法，得到AGI坐标系的轴单位向量为y~GER=x^EPGER-(x^EPGER·x^EMGER)x^EMGER||x^EPGER-(x^EPGER·x^EMGER)x^EMGER||---(18)]]>最后可以得到轴的单位向量为z~GER=x~GER×y~GER---(19)]]>在此，地心E、飞行器当前位置M和预测命中点P组成的平面为MEP平面，定义MEP平面与地球表面的交线为广义赤道平面，与广义赤道正交的大圆弧的一半定义为广义子午线，可知广义子午线端点在轴上；通过此定义，再入飞行器CAV的位置用广义经度广义纬度和广义高度表示；为分析方便，假设广义的0°子午线通过CAV初始再入点，故CAV初始位置表示为λ~0=0φ~0=0H~0=H---(20)]]>尽管由于AGI坐标系随着CAV的运动而旋转，在此为方便分析，可认为是静止的；CAV相对于AGI坐标系的速度，用表示，是相对地球速度和牵连速度的矢量和，初始速度为V~0NED=Vcos(γ)cos(ψ)Vcos(γ)sin(ψ)+ωe(Re+H)cos(φ)-Vsin(γ)---(21)]]>在此，上标“NED”代表向量在NED坐标系中的分量形式；定义的大小为广义速度定义和当地水平面夹角为广义弹道倾角定义的水平分量于当地广义子午线之间夹角为广义航向角由式(21)得到，初始的和为V~0=[V2+2Vωe(Re+H)cos(φ)cos(γ)sin(ψ)+ωe2(Re+H)2cos2(φ)]12---(22)]]>γ~0=sin-1(Vsin(γ)/V~0)---(23)]]>由于亚轨道再入初始点M点沿广义子午线上的切线与平行，因此与当地水平分量之间夹角即为初始航向角为ψ~0=cos-1(z~GER·V~H0GER/||V~H0GER||);y~GER·V~H0GER≥0-cos-1(z~GER·V~H0GER/||V~H0GER||);y~GER·V~H0GER<0---(24)]]>上式中的和是由公式(18)和(19)得到的，由如下公式得到V~H0GER=(TGERNED)TV~H0NED---(25)]]>V~H0NED=Vcos(γ)cos(ψ)Vcos(γ)sin(ψ)+ωe(Re+H)cos(φ)0---(26)]]>上式中，是由公式(7)给出P点的广义经纬度和高度为λ~P=cos-1(x^EPGER·x~GER)φ~P=0H~P=HT---(27)]]>前面说明了对旋转地球模型下的再入制导的末端条件，这些条件需要转移到AGI坐标系中，定义为x~1GER=x^EPGER-x~GER---(28)]]>通过运用格莱姆‑施密特正交化方法，得到与正交的向量为x2GER=x1GER-(x1GER·x^EPGER)x^EPGER---(29)]]>定义x3为与过P点与水平面、广义赤道交线都平行的单位向量，x3在NEDP坐标系中表示为x3NEDP=TGERNEDPx2GER/||x2GER||---(30)]]>上式中，缩写词“NEDP”代表在位置P处的NED坐标系，将λ＝λP和φ＝φP带入到公式(7)中，即得到从GER坐标系到NEDP坐标系的坐标转换矩阵定义相对于AGI坐标系的最终速度向量为在AGI坐标系中，再入制导方法导引飞行器近似沿着广义赤道到达P，并且使得广义弹道倾角接近0，在此为方便分析，假设与x3平行，制导律也能够通过最后一次倾侧角翻转来消除其引起的误差；定义相对于旋转地球的飞行器最终速度为VTAEM，除此之外，为得到地球旋转引起的牵连速度假设滑翔段结束距离P的距离STAEM远小于地球半径，指向当地东向，大小为VeP=ωe(Re+HTAEM)cos(φP)---(31)]]>速度向量的大小为在以上前提情况下，可得V~f=V~fx3---(32)]]>(V~f-VeP)2=(VTAEM)2---(33)]]>将上式展开得V~f2-2V~fVePcos(θP)-2(VeP)2=VTAEM2---(34)]]>cos(θP)=x3NEDP|y---(35)]]>上式中，是在y轴上的分量，进而得到的估计值为V~f=VePcos(θP)+(VeP)2(cos2(θP)-1)+VTAEM2---(36)]]>定义相对于惯性空间AGI坐标系的单位质量机械能为表示为E~=12V~2-μRe+H~---(37)]]>式中，μ是重力常数，则最终的末端约束条件转化成相对于惯性空间的能量为E~f=12V~f2-μRe+HTAEM---(38)]]>步骤4：运动方程线性化；在分析时，AGI坐标系认为是静止的，得到运动方程为dλ~dt=V~cos(γ~)sin(ψ~)(Re+H~)cos(φ~)---(39)]]>dφ~dt=V~cos(γ~)cos(ψ~)(Re+H~)---(40)]]>dH~dt=V~sin(γ~)---(41)]]>dV~dt=-Dm-gsin(γ~)---(42)]]>dγ~dt=1V~[Lcos(σ)m-gcos(γ~)+V~2cos(γ~)Re+H~]---(43)]]>dψ~dt=1V~[Lsin(σ)mcos(γ~)+V~2cos(γ~)sin(ψ~)tan(φ~)Re+H~]---(44)]]>其中能量对时间的导数为dE~dt=-DV~m---(45)]]>定义射程xD和横向射程xC为和将公式(39)(40)(44)与公式(45)联立得到dxDdE~=-mcos(γ~)sin(ψ~)Dcos(φ~)Re(Re+H~)---(46)]]>dxCdE~=-mcos(γ~)cos(ψ~)DRe(Re+H~)---(47)]]>dψ~dE~=-m~DV~2[Lsin(σ)mcos(γ~)+V~2cos(γ~)sin(ψ~)Re+H~tan(xCRe)]---(48)]]>令L1＝Lcos(σ)，L2＝Lsin(σ)；假设飞行器以平稳滑翔条件即Steady Glide Condition,SGC滑翔飞行，飞行器的广义弹道倾角变化率为0，即由公式(43)得L1m-gcos(γ~)+V~2cos(γ~)Re+H~=0---(49)]]>经过转换可得cos(γ~)=L1mg-mV~2Re+H~---(50)]]>由式(37)可得V~2=2E~+2μRe+H~---(51)]]>在AGI坐标系中，由于制导律导引飞行器近似沿着广义赤道飞向P，可知因此假设和tan(xC/Re)＝xC/Re，将式(50)、(51)带入到式(46)、(47)和(48)，运用上述假设，得dxDdE~=-L1DRe-μ(Re+H~)-2E~---(52)]]>dxCdE~=L1DReψ~-μRe+H~-2E~-π2L1DRe-μRe+H~-2E~---(53)]]>dψ~dE~=L1D1μRe+H~+2E~xCRe-L2D12E~+2μRe+H~---(54)]]>为减小微分计算的复杂性，式(48)中的分母假设为1；步骤5：求解解析解；由于对式(52)(53)和(54)的解影响很小，故在分析时，设为一常量令R*＝Re+H*；定义纵向升阻比为L1/D＝Lcos(σ)/D，将L1/D写成如下形式L1D=a2E~2+a1E~+a0---(55)]]>定义水平升阻比为L2/D＝Lsin(σ)/D，将L2/D写成如下形式L2D=b2E~2+b1E~+b0---(56)]]>在此，设计L1/D和L2/D曲线，通过同时调节攻角和倾侧角来跟踪参考升阻比曲线；①射程表达式将式(55)带入式(52)，积分可得射程的表达式为xD(E~,E~0)=14a2Re(E~2-E~02)+12(a1-a2μ2R*)Re(E~-E~0)+4a0(R*)2-2μR*a1+μ2a28(R*)2Reln(2R*E~+μ2R*E~0+μ)---(57)]]>②横向射程和方位角表达式将式(53)和(54)结合，得到dxCdE~dψ~dE~=L1D1μR*+2E~0-Re1Re0xCψ~+π2L1DReμR*+2E~-L2D12E~+2μR*T---(58)]]>令f1(E~)=L1D1μR*+2E~---(59)]]>f2(E~)=π2L1DReμR*+2E~---(60)]]>f3(E~)=-L2D12E~+2μR*---(61)]]>f4(E~,E~0)=∫E0E~f1(x~E)dx~E=xD(E~,E~0)Re---(62)]]>A=0-Re1Re0---(63)]]>式(58)写成如下形式dxCdE~dψ~dE~=f1(E~)AxCψ~+f2(E~)f3(E~)---(64)]]>引入谱分解方法来求解上述特殊类型的变系数线性系统，定义为Q(E~,E~0)=exp(-A∫E~0E~f1(x~E)dx~E)=exp(-Af4(E~,E~0))---(65)]]>在式(64)两边乘以得到exp(-Af4(E~,E~0))dxCdE~dψ~dE~-exp(-Af4(E~,E~0))f1(E~)AxCψ~=exp(-Af4(E~,E~0))f2(E~)f3(E~)---(66)]]>上式(66)写成exp(-Af4(E~,E~0))dxCdE~dψ~dE~+ddE~[-Af4(E~,E~0)]xCψ~=exp(-Af4(E~,E~0))f2(E~)f3(E~)---(67)]]>为求取其微分，反向利用点乘法则，可得ddE~[-Af4(E~,E~0)]xCψ~=exp(-Af4(E~,E~0))f2(E~)f3(E~)---(68)]]>对式(68)两边积分，得exp(-Af4(E~,E~0))xCψ~-exp(-Af4(E~,E~0))xC0ψ~0=∫E~0E~exp(-Af4(E~,E~0))f2(E~)f3(E~)dE~---(69)]]>注意到其中02×2是2×2阶零矩阵，I2×2为2×2单位阵，的逆为Φ(E~,E~0)=[Q(E~,E~0)]-1=exp(Af4(E~,E~0))---(70)]]>将式(70)左乘以式(69)得xCψ~=Φ(E~,E~0)xC0ψ0+∫E~0E~Φ(E~,x~E)f2(x~E)f3(x~E)dx~E---(71)]]>上式中，为状态转移矩阵，通过对矩阵A谱分解，由矩阵论知识得Φ(E~,E~0)=exp(Af4(E~,E~0))=exp(λ1f4(E~,E~0))G1+exp(λ2f4(E~,E~0))G2---(72)]]>上式中，λ1＝i和λ2＝‑i是矩阵A的特征值，其中G1和G2是矩阵A的谱投影，为G1=A-λ2Iλ1-λ2=1/2Rei/2-i/(2Re)1/2---(73)]]>G2=A-λ1Iλ2-λ1=1/2-Rei/2i/(2Re)1/2---(74)]]>进而可得Φ(E~,E~0)=cos(f4(E~,E~0))-Resin(f4(E~,E~0))sin(f4(E~,E~0))/Recos(f4(E~,E~0))---(75)]]>和∫E~0E~Φ(E~,x~E)f2(x~E)f3(x~E)dx~E=∫E~0E~cos(f4(E~,x~E))f2(x~E)-Resin(f4(E~,x~E))f3(x~E)sin(f4(E~,x~E))f2(x~E)/Re+cos(f4(E~,x~E))f3(x~E)dx~E---(76)]]>由式(60)和(62)，有∫E~0E~cos(f4(E~,x~E))f2(x~E)dx~E=-πRe2∫E~0E~cos(f4(E~,x~E))df4(E~,x~E)=πRe2sin(f4(E~,E~0))---(77)]]>1Re∫E~0E~sin(f4(E~,x~E))f2(x~E)dx~E=-π2∫E~0E~sin(f4(E~,x~E))df4(E,x~E)=π2-π2cos(f4(E~,E~0))---(78)]]>这样就求得，横向射程和航向角的解析公式为xC(E~,E~0)=xC0cos(f4(E~,E~0))-Reψ~0sin(f4(E~,E~0))+πRe2sin(f4(E~,E~0))-Re∫E~0E~sin(f4(E~,x~E))f3(x~E)dx~E---(79)]]>ψ~(E~,E~0)=xC0Resin(f4(E~,E~0))ψ~0cos(f4(E~,E~0))+π2-π2cos(f4(E~,E~0))+∫E~0E~cos(f4(E~,x~E))f3(x~E)dx~E---(80)]]>在此，需要提到的是在求解再入制导方法时，飞行器飞行的攻角为方案攻角，通过调节倾侧角来跟踪规划的L1/D曲线，这意味着L2/D曲线不能任意规划，也不能像式(56)一样用有限次数多项式表示；然而由于L2/D的表达式只包含并没有在上述微分中体现，因此上述两个解析解仍然适用；步骤6：验证上述解析解可行性；该基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法，具体步骤如下：步骤A：规划快速下降段的制导方案；CAV在与运载器分离之后进入无动力滑翔状态，开始进入快速下降段，直到开始满足平稳滑翔条件即SGC时结束；在这一段，由于大气密度ρ非常小，飞行器快速下降，高度快速降低；随着高度下降，大气密度升高，飞行器的热流率急剧升高，而最大热流率大致出现在这段的结束；为了能满足热流率约束，设计此段飞行器飞行时保持最大攻角(αmax)并保持倾侧角为0，这样使得在下降段飞行的更高，从而降低热流率；当时，说明此时开始升力已经足够大，满足平稳滑翔条件；为保证飞行器能从DP平滑转移到SGP，下式(81)给出攻角方案，指令攻角和倾侧角为αcmd=αmax,γ·<0deg/sΔγΔγ1αmax+Δγ1-ΔγΔγ1[αplan+kγΔγ],γ·≥0deg/s---(81)]]>σcmd＝σplan＝0deg                   (82)Δγ＝γplan‑γ                  (83)上式中，αplan是SGP的规划的攻角，在式(84)中给出；γplan是平稳滑翔的弹道倾角，γ为当前飞行器的弹道倾角，Δγ1是当时的Δγ，kγ是值为5的常数；由式(81)可知，Δγ从Δγ1逐渐趋近于0的同时，指令攻角αcmd从αmax平滑趋向于αplan；当Δγ＝0，时，下降段DP结束，因为之后飞行器将向上飞，在此处飞行器相对于惯性空间的单位质量接机械能表示为步骤B：规划平稳滑翔段的制导方案；在平稳滑翔段，飞行器沿着规划的攻角曲线飞行；制导律利用射程方向的解析表达式实时更新L1/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=12V2-μRe+H---(84)]]>为了能够飞出更远的射程，飞行器在平稳滑翔段大部分时间内以α1＝10deg飞行，这样近似保持最大的升阻比(L/D)max；当飞行器飞行到轨迹末段时，平滑的将αplan减小到α2＝6deg，将高度H降低到HTAEM进入末制导段，这样保证飞行器有足够的动压满足末端机动飞行，在此设计αplan为αplan={α1;E2≤E≤E0(E2-EE2-Ef)2(α2-α1)+α1;Ef≤E<E2---(85)]]>上式中E0是初始的能量，平稳滑翔飞行结束时的E2可设为‑5.6×107J/kg，Ef是由终端约束所决定的能量，通过将VTAEM和HTAEM带入式(84)中得到；当αplan曲线确定后，就确定了高度走廊，高度走廊的下边界(Hmin)由前面所述的过程约束决定；除此之外忽略地球旋转的影响，假设倾侧角σ＝0，这样通过求解式(5)得到高度走廊的上边界(Hmax)；下面给出相对能量E的升阻比(L/D)plan曲线，由于αplan曲线已经确定，(L/D)plan是马赫数Ma的函数，而Ma是速度V和高度H的函数；所以，确定了高度曲线，就得到相应的(L/D)plan曲线；实际上，当飞行器高速滑翔时，升阻比(L/D)很轻微的随Ma变化，但是由于高度曲线被限定在很小的走廊中，H对Ma影响很小，所以高度对(L/D)plan曲线影响很小；在此，假设飞行器沿着高度走廊的中间飞行，高度为H(E)=Hmax(E)+Hmin(E)2---(86)]]>当给定E和H时，得到V和Ma，然后利用αplan和Ma，得到相对于E的(L/D)plan曲线，最大的升阻比可达3；在利用解析解快速生成参考弹道的过程中，为了确定第一次翻转节点需要知道在惯性空间中的基准升阻比曲线这就需要建立E与之间的转换关系，由式22、37和91可知，E与之间的关系为E~=E+Vωe(Re+H)cos(φ)cos(γ)sin(ψ)+0.5ωe2(Re+H)2cos2(φ)---(87)]]>由于在实际飞行中，飞行任务剩下的弹道曲线并不明确，故上式中的关系并不实用，通常E0=-3.5×104kJ/kgEf=-6×104kJ/kgV0ωe(Re+H0)≈3.3×103kJ/kgVfωe(Re+Hf)≈0.93×103kJ/kgωe2(Re+H0)2≈0.22×103kJ/kg---(88)]]>从中看出式(87)中的非线性项远小于线性项，即线性项占主要作用，因此在实际仿真结果中与E看起来近似直线；这里与E之间近似采用如下线性转化关系xE=Ef-EE~f-E~(x~E-E~)+E---(89)]]>上式中，E与为当前相对旋转地球和相对于惯性空间的能量值，Ef与即式(37)为相对旋转地球和相对于惯性空间的终端能量值；xE和任意时刻相对旋转地球和相对于惯性空间的能量值，由于(L/D)plan曲线一大部分是常值，由式(89)确立的误差较小，并且误差通过第二次倾侧角翻转来修正，在此，令与相关的(L/D)plan曲线为由下式得到在下式(91)中，当指令攻角αplan曲线和高度曲线确立之后，将其转化为最大倾侧角σmax约束；当飞行器以平稳滑翔条件飞行时，增大σ会使得高度H降低，当将H降低到Hmin时，即得到最大的σmax；如式(92)所示，在平稳滑翔时，纵向升力L1与重力和向心力近似平衡，由式(91)中右边第一项，求得σmax；若飞行轨迹有振荡，式(91)右边第二项可以消除振荡，因此，若高度H快速下降，到达(H)E‑(Hmin)E>0，式(91)右边第二项会降低σmax，使得飞行器能够拉起，保持在Hmin之上；在此，函数f对E求导，写成fE；σmax=arccos(L1Lmax)+kσ(dHdE-dHmindE)---(91)]]>L1=m(g-V~02Re+Hmin(E))---(92)]]>Lmax＝CLplan[0.5ρ(Hmin)V2]Sref            (93)上式中，kσ设为‑50，CLplan为αplan的升力系数，由式(3)(4)和(84)可知，E的变化率为dEdt=VV·+gH·=-DVm+Vωe2(Re+H)cos2(φ)sin(γ)-Vωe2(Re+H)sin(φ)cos(φ)cos(γ)cos(ψ)---(94)]]>忽略地球旋转影响，由式(3)和式(94)，得dHdE≈-msin(γ)D---(95)]]>(2)规划纵向升阻比L1/D曲线再入制导方法利用射程上的解析解实时规划纵向升阻比曲线，然后通过调节倾侧角跟踪规划的升阻比曲线，以满足终端速度约束；当E>E2时，(L/D)plan认为是常量，故时，设计纵向升阻比L1/D曲线为水平段；当时，(L/D)plan随着E下降而下降，导致L1/D不能认为是常数；但是，在此可认为(L1/D)1到(L1/D)2的转变是线性的，L1/D曲线设计为L1/D(x~E)={(L1/D)1;E~2≤x~E≤E~1x~E-E~fE~2-E~f(L1/D)1+E~2-x~EE~2-E~f(L1/D)2;E~f≤x~E<E~2---(96)]]>上式中，是将xE＝E2带入到式(89)中得到的，(L1/D)1和(L1/D)2是需要实时修正的参量，由于倾侧角σ的余弦为纵向升阻比与总升阻比之比，而且最终的倾侧角约束要求接近0，(L1/D)2则由下式实时计算(L1/D)2=(L/D)plan(Ef)(L/D)est(E)(L/D)plan(E)---(97)]]>其中，(L/D)est(E)是由惯性导航系统即Inertial Navigation System,INS得到的，(L/D)est(E)是在当前能量下的规划升阻比，(L/D)plan(Ef)是在终点能量下的规划升阻比；在AGI坐标系中，剩余射程为xDf=Reλ~P-STAEM---(98)]]>上式中，由式(26)求取，下面计算(L1/D)1：a.当将式(96)带入式(52)，然后积分，得xD(E~2,E~)+xD(E~f,E~2)=xDf---(99)]]>其中当时有xD(E~2,E~)=(L1/D)1Re2ln(2R*E~2+μ2R*E~+μ)---(100)]]>当时有xD(E~f,E~)=Re[(L1/D)1-(L1/D)2](E~f-E~)2(E~2-E~f)-Re(L1/D)12(E~2-E~f)(E~f+μ2R*)ln(2R*E~f+μ2R*E~+μ)+Re(L1/D)22(E~2-E~f)(E~2+μ2R*)ln(2R*E~f+μ2R*E~+μ)---(101)]]>其中，将带入式(101)得通过求解式(99)，得(L1/D)1=c1c2---(102)]]>其中c1=xDf-12Re(L1/D)2-Re(L1/D)22(E~2-E~f)(E~2+μ2R*)ln(2R*E~f+μ2R*E~2+μ)---(103)]]>c2=Re2ln(2R*E~2+μ2R*E~+μ)-12Re-12ReE~2-E~f(E~f+μ2R*)ln(2R*E~f+μ2R*E~2+μ)---(104)]]>b.当在这种情形下，由于大部分时间飞行器处于AAP阶段，不需要跟踪纵向升阻比L1/D曲线，因此其不需要再更新；(3)确定第一次倾侧角翻转节点再入制导方法利用横向射程的解析解表达式(79)，来决定第一次倾侧角翻转的能量节点倾侧角翻转是规划的，规划的横向升阻比L2/D曲线为L2/D(x~E){sgn·|(L2/D)(x~E)|x~E>E~3orE~f≤x~E≤E~4-sgn·|(L2/D)(x~E)|E~4<x~E≤E~3---(105)]]>上式中，是第二次倾侧角翻转的能量节点，在此假设剩余飞行中的L/D微分由INS给出，预测的为sgn的符号由再入初始条件决定sgn={1;ψ~0∈[0,π/2)-1;ψ~0∈[π/2,π]---(107)]]>从上面的结果看出，这里安排两次倾侧角翻转，第一次在第二次在并且第二次倾侧角翻转离终点状态较近，这样选取接近末端的有下列好处：的调节修正之前积累的误差，并且使得后面剩余较短射程内的误差较小，近似可以忽略；为方便利用横向射程解，构造如下积分函数F(x~E2,x~E1)=∫x~E1x~E2sin(f4(E~f,x~E))f5(x~E)dx~E---(108)]]>f5(x~E)=-|(L2/D)(x~E)|2x~E+2μR*---(109)]]>上述积分函数与横向射程解析解积分形式类似，由于纵向升阻比L1/D曲线是分段函数，写成如下分段函数形式f4(E~f,x~E)={[xD(E~f,E~2)+xD(E~2,x~E)]/Re;x~E≥E~2xD(E~f,x~E)/Re;x~E<E~2---(110)]]>其中，和由式(101)求得，由式(100)求得；由于假设最终的横向射程解析式仅是的函数，将初始条件式(20)、(21)、(22)、(23)和L2/D曲线式(105)、(106)、(107)带入横向射程解析式(79)得到最终形式为xCf(E~3)=(π2-ψ~0)Resin(xDfRe)-sgnReF(E~f,E~)+2sgnReF(E~4,E~3)---(111)]]>其导数为xCf′(E~3)=-2sgnResin[f4(E~f,E~3)]f5(E~3)---(112)]]>在此，运用牛顿法求解的解E~3(k+1)=E~3(k)-xCf(E~3(k))xCf′(E~3(k))---(113)]]>为了提升运算效率，做三次运算；第一次运算是在平稳滑翔段的开始，得到当T1＝200s，时，第二次计算，得到当T2＝30s，时，第三次计算得到在这里，是因为倾侧角变化率的限制所加的延迟，在式(115)中给出，aD是由惯导系统所测得的当前的阻力加速度；(4)确定基准倾侧角由于L1＝Lcosσ，L1/D＝L/Dcosσ，为跟踪纵向升阻比L1/D曲线，基准的倾侧角为σplan=sgn·cos-1(L1/D(L/D)est);E~>E~3+ΔE~-sgn·cos-1(L1/D(L/D)est);E~4+ΔE~<E~≤E~3+ΔE~---(114)]]>其中，是考虑到飞行器的滚转速率的限制，提前时间Δt翻转基准倾侧角，其值为ΔE~=aDV~0Δt---(115)]]>上式中，Δt完成倾侧角翻转所需时间的一半，用下式估计Δt=|σσ·max|---(116)]]>其中，σ是当前倾侧角，是最大滚转速率；(5)攻角指令αcmd和倾侧角指令σcmd若将基准攻角αplan和倾侧角σplan直接用作攻角和倾侧角指令，再入飞行轨迹将会有微弱的长周期振动，在此，下面给出一种消除长周期振动的方法；升力系数的垂直分量CL1和水平分量CL2，分别由基准攻角αplan和倾侧角σplan决定；为消除振动，在垂直方向上附加一个与垂直速度分量相反的升力，即在CL1上加上ΔCL1，设计ΔCL1为ΔCL1=CLαkγ(γplan-γ)---(117)]]>其中，是升力线斜率，γplan为平稳滑翔所需要的近似弹道倾角，在后面式(126)中给出；kγ为增益系数，值为5，若飞行器上升较快，γplan‑γ<0，ΔCL1<0，造成升力L1降低，从而抑制飞行器上升过快，同样的可以抑制下降过快；c.指令攻角αcmd推导对于给定的马赫数Ma，攻角是升力系数CL的函数，α＝fα(CL)，其一阶Taylor展开为αcmd=fα((CL1+ΔCL1)2+CL22≈fα((CL1)2+CL22)+fα′((CL1)2+CL22)CL1(CL1)2+CL22ΔCL1---(118)]]>注意到有fα′(CL12+CL22)=dαdCL=1CLα---(119)]]>将式(117)和(119)代入到式(118)中，得αcmd＝αplan+cos(σplan)kγ(γplan‑γ)               (120)d.指令倾侧角σcmd推导σcmd=arctan(CL2CL1+ΔCL1)---(121)]]>其一阶Taylor近似为σcmd≈arctan(CL2CL1)-CL2CL12+CL22ΔCL1=σplan-CL2CLCLαCLkγ(γplan-γ)≈σplan-sin(σplan)kγ(γplan-γ)αplan1---(122)]]>其中，对于CAV，令αplan＝10deg；为了保证飞行中满足各个过程约束，倾侧角有最大界限为σcmd={σmax,σcmd>σmax-σmax,σcmd<-σmax---(123)]]>其中，σmax为最大倾侧角约束，由式(91)得到；e.基准弹道倾角γplan推导γplan为平稳滑翔的近似弹道倾角，由式(5)，假设忽略地球自转影响，得Lplancos(σplan)-mgcos(γplan)+mV2cos(γplan)Re+H=0---(124)]]>将上式对时间求导得1mdCLplandEdEdt12ρV2Srefcos(σplan)+1mCLplan12dρdHdHdtV2Srefcos(σplan)+1mCLplanρVdVdtSrefcos(σplan)+1mLplanddE[cos(σplan)]dEdt-dgdHdHdtcos(γplan)+gdγplandtsin(γplan)+1R0+H[2VdVdtcos(γplan)-V2sin(γplan)dγplandt]-1(R0+H)2V2dHdtcos(γplan)=0---(125)]]>其中，由式(94)得dEdt=-DVm---(126)]]>为简化式(125)，假设γplan和g的导数为0，除此之外，γplan≈0，令cos(γplan)＝1，sin(γplan)＝γplan，将式(3)和(4)带入式(125)中，忽略地球自转影响，得到γplan=-Dplanmgd1d2---(127)]]>上式中，Dplan＝CDplanqSref，CDplan是攻角为αplan时的阻力系数，ρ是大气密度，Sref为参考面积，d1和d2分别表示如下d1=ρV2Srefcos(σplan)2mdCLplandE+2R0+H+CLplanρSrefcos(σplan)m+LplanmddE[cos(σplan)]---(128)]]>d2=-CLplanV2Srefcos(σplan)2mgdρdH+2R0+H+CLplanρSrefcos(σplan)m+V2(R0+H)2g---(129)]]>前面提到，将函数f对E的导数定为fE；因为CLplan(E)已经规划，(CLplan)E是已知的；在平稳滑翔段用式(130)求得[cos(σplan)]E；然而，在弹道下压段，由于[cos(σplan)]E的计算公式较复杂，用相邻两时刻的cos(σplan)的差分计算；ddE[cos(σplan)]=ddE[L1/D(L/D)est]=1(L/D)estdL1/DdE-L1/D(L/D)est(L/D)pland(L/D)plandE---(130)]]>其中，(L/D)plan(E)是规划曲线，[(L/D)plan]E是已知的，(L1/D)E可由式(89)和(96)得到；步骤C：规划弹道下压段的制导方案；飞行器经过第二次倾侧角翻转之后，进入弹道下压段即AAP；通过降低攻角，飞行器快速进入稠密大气层，这样能够使得飞行器在末制导段有较大的动压进行机动；这就使得在这段中，不能近似接近0，所以使得上述解析解难以满足终端约束条件；在此用一种与前面不同的制导方案：在最后一次倾侧角翻转之前，为得到需求的最终速度，在线的弹道仿真修正能量节点在最后一次倾侧角翻转之后，采用比例导引法即Proportional Navigation,PN导引，下面首先介绍弹道下压段的制导方案，再给出修正的方法；①基准倾侧角在弹道下压段，αplan曲线与平稳滑翔段相同，而倾侧角由比例导引律确定；当飞行器接近目标时，忽略地球曲率的影响，MP视线的方位角的变化率为ψ·MP=V~0cos(γ~0)sin(π/2-ψ~0)xDP---(131)]]>其中，由比例导引律得到的需用横向机动加速度为aL2=kPNψ·MPV~0cos(γ~0)---(132)]]>其中，为了防止初始倾侧角饱和，令kPN从2逐渐变化到4，为kPN=2xDfxDfE4+4(1-xDfxDfE4)---(133)]]>其中，xDf是第二次倾侧角翻转的能量节点；另外，升力的垂直分量需要与重力和向心力平衡，可得aL1≈g-V~02Re+H---(134)]]>进而可以得到基准倾侧角为σplan=sgn·arctan(aL2aL1)E~≤E~4+ΔE~---(135)]]>其中，由式(115)求得；②指令攻角和指令倾侧角在飞行器进入弹道下压段之前，需要获得射程‑能量参考曲线并且需要修正在弹道下压段，通过调节攻角跟踪曲线，如式(136)所示，消除由于干扰引起的最终速度误差；注意到在弹道下压段中，大部分时间升阻比L/D对攻角导数大于0，意味着增加攻角会增加升阻比，造成最终速度升高，为了消除长周期振动，σcmd与式(122)一样；αcmd=αplan+kα[sgo(E)-sgoref(E)]---(136)]]>σcmd=σplan-sin(σplan)kγ(γplan-γ)αplan1---(137)]]>其中，kα值为(5π/18)×10‑6，也就是说0.5deg的攻角修正，会引起10km的射程散布；在弹道下压段中，在计算式(127)中的γplan时，由于[cos(σplan)]E表达式较复杂，其对E的导数由差分求得；同时，为了满足过程约束，仍然需要最大倾侧角约束σcmd={σmax,σcmd>σmax-σmax,σcmd<-σmax---(138)]]>③第二次倾侧角翻转在弹道下压段，不需要跟踪L1/D曲线，需要调整翻转能量节点来满足最终速度要求；下面分析与最终速度Vf之间的关系：若降低则翻转时间推迟，时刻的航向角误差增加，为消除航向角误差，需要更大的横向机动加速度，因而倾侧角的大小会增大，造成纵向升阻比L1/D减小，导致最终速度降低；所以，降低造成最终速度Vf降低；在第一次倾侧角翻转之后，在此采用预测校正法修正通过弹道仿真预测最终速度，之后运用secant法修正当前的状态作为仿真的初始状态，由于INS能够实时测量当前气动加速度，通过对比实际和理想气动加速度，得到当前气动系数的偏差；弹道预测仿真假设当前估计的偏差即为剩余飞行段的气动系数偏差，除此之外，仿真中忽略式(136)中右边的第二项，当Sgo＝STAEM时，仿真结束，得到预测的最终速度修正的secant法式为E~4(k+1)=E~4(k)-(Vf(E~4(k))-VTAEM)(E~4(k)-E~4(k-1))(Vf(E~4(k))-Vf(E~4(k-1)))---(139)]]>其中，为第k次仿真的为减小运算次数，的修正有两次：第一次预测校正过程在第一次倾侧角翻转之后，得到T1＝60s时，此时进行第二次预测校正过程，得到之后可令当CAV飞行到Sgo＝STAEM时，仿真结束。