“再入制导”是指高超声速飞行器按照一定的导引规律将飞行器从初始投放点引导到目标状态点，期间满足多种过程约束、控制约束和终端约束。从本质上讲，轨迹优化设计问题是一个复杂的非线性最优控制问题，它的实时性规划一直是难题。传统的标称轨迹制导方法采用离线轨迹设计与在线标称轨迹跟踪相结合的制导方案，此举规避了轨迹在线设计的难题，使得再入制导工程化成为可能^[1]。然而，标称轨迹制导方法的离线轨迹设计的方案决定了它只能应用于特定场景下的特定任务，从根本上限制了制导的鲁棒性以及任务的灵活性，随着近年来任务需求的不断增强，它的弊端不断突显，包括耗费人力物力、鲁棒性有限、任务灵活性低等。

20世纪90年代末，在新一代可重复使用运载飞行器发展背景下，美国NASA牵头启动了先进制导与控制计划(AG & C)^[2], 旨在为下一代飞行器提供更加自主、可靠与省钱的再入制导，掀起了新型再入制导方法研究的热潮。新型再入制导技术的发展主要围绕2个方向，包括标称轨迹在线重构^[3-7]以及预测校正制导方法^[8-12]。2种方法中，可行轨迹在线快速设计都是技术关键。1997年，Lu等^[13-14]引用伪谱法概念以阻力为自变量，以总吸热量最小为优化目标，采用成熟的序列二次规划方法实现快速阻力加速度生成。2003年，Shen等^[10]提出了拟平衡滑翔条件，基于约束生成的倾侧角走廊实现特定射程的倾侧角剖面快速生成。2014年，Lu^[15]在总结近年来轨迹设计结果的基础上，进一步提出了一套通用的轨迹快速规划方案，并率先提出采用时标分离理论对轨迹进行修正，可以得到比较平缓的再入轨迹。在文[16]的工作基础上，雍恩米等^[17]进一步研究考虑地球旋转影响和大横程的三维参考轨迹快速生成。以上方法取得不同技术突破，但是由于设计参数与设计目标之间没有显式的函数关系，算法大多采取周期内多次弹道积分，并依据轨迹预测信息提取出来的梯度信息迭代到目标参数值。这种周期内多次轨迹预测和迭代的方法给实际制导工程化带来严峻的挑战。与直接打靶法相比，伪谱法同时离散状态和控制变量，不需要对状态积分，具有求解速度更快的显著优势，在轨迹快速规划方面具有广泛应用前景，受到研究者的持续关注^[18-22]。

结合以上研究现状分析和问题描述，本文提出一种基于在线模型辨识的可行轨迹实时规划算法被提出。该算法将不断产生于周期间的轨迹预测信息收集起来，并采用递推最小二乘估计算法对求解模型在线辨识。辨识得到的模型可以在周期间传递，进而应用于可行轨迹的参数在线修正。不断修正的可行轨迹与经典的标称轨迹跟踪制导方法结合起来，可以构成一种新型复合制导方案。制导多任务仿真验证了整个方案的有效性和可行性。

1 再入制导问题 1.1 再入运动模型

本论文以高超声速滑翔式再入飞行器为研究对象，其考虑地球自转的3自由度再入运动模型表示如下^[23]：

$ \left\{ \begin{array}{l} \dot r = v\sin \theta ,\\ \dot \lambda = \frac{v}{{r\cos \phi }}\cos \theta \sin \psi ,\\ \dot \phi = \frac{v}{r}\cos \theta \cos \psi ,\\ {{\dot S}_e} = \frac{v}{r}\cos \theta ,\\ v = - \frac{D}{{{m_0}}} - g\sin \theta + \omega _e^2r\cos \phi \left( {\sin \theta \cos \phi - \cos \theta \sin \phi \cos \psi } \right),\\ \dot \theta = \frac{1}{v}\left[ {\frac{{L\cos \sigma }}{{{m_0}}} + \left( {\frac{{{v^2}}}{r} - g} \right)\cos \theta + 2{\omega _e}v\cos \phi \sin \psi + } \right.\\ \;\;\;\;\;\left. {\omega _e^2r\cos \phi \left( {\cos \theta \cos \phi + \sin \theta \sin \phi \cos \psi } \right)} \right],\\ \dot \psi = \frac{1}{v}\left[ {\frac{{L\sin \sigma }}{{{m_0}\cos \theta }} + \frac{{{v^2}}}{r}\cos \theta \sin \psi \tan \phi + \frac{{\omega _e^2r}}{{\cos \theta }}\cos \psi \sin \phi \cos \phi - } \right.\\ \;\;\;\;\;\;\left. {2{\omega _e}v\left( {\cos \phi \cos \psi \tan \theta - \sin \phi } \right)} \right]. \end{array} \right. $

(1)

其中：r为地心距；R₀为地球半径；λ和ϕ分别为经度和纬度；v为飞行器相对于地球大气的速度；θ为速度矢量与水平面的夹角，向上为正；ψ为速度矢量与当地正北方向的夹角，顺时针为正；σ为飞行器倾侧角，右偏为正；ω_e为地球自转角速度；S_e是单位为弧度的射程角，衡量飞行器累积射程；m₀为飞行器质量，保持恒定；g表示飞行器受到地球的引力加速度值；L=0.5ρv²C_LS_ref和D=0.5ρv²C_DS_ref为飞行器升力和阻力大小；C_L和C_D为升力和阻力系数，与攻角和Mach有关；S_ref为气动参考面积。攻角为速度的函数，本论文采用三段式：

$ {\alpha _{{\rm{ref}}}} = \left\{ \begin{array}{l} {\alpha _{\max }},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;v > {v_1};\\ \frac{{{\alpha _{\max }} - {\alpha _{\min }}}}{{{v_1} - {v_2}}}\left( {v - {v_2}} \right),\;\;\;\;\;{v_2} < v < {v_1};\\ {\alpha _{\min }},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;v < {v_2}. \end{array} \right. $

(2)

其中：最大攻角α_max=20°, 最小攻角α_min=8.5°, 2个速度节点v₁=4 700 m/s，v₂=3 100 m/s。

1.2 再入轨迹约束条件

再入过程约束主要包括热峰值约束、过载约束和动压约束^[16]。

$ \dot Q\left( \tau \right) = \frac{{{C_1}}}{{\sqrt {{R_{\rm{d}}}} }}{\left( {\frac{\rho }{{{\rho _0}}}} \right)^{0.5}}{\left( {\frac{v}{{{V_C}}}} \right)^{3.15}} \le {\dot Q_{\max }}, $

(3)

$ n\left( \tau \right) = \sqrt {{{\bar L}^2} + {{\bar D}^2}} = q\sqrt {C_L^2 + C_D^2} \frac{{{S_{{\rm{ref}}}}}}{{{m_0}{g_0}}} \le {n_{\max }}, $

(4)

$ q\left( \tau \right) = \frac{1}{2}\rho {v^2} \le {q_{\max }}. $

(5)

其中：C₁=11 030，为与飞行器特性相关的常数；R_d=0.1 m，为飞行器头部曲率半径；V_C为第一宇宙速度，大小为$\sqrt {{R_0}{g_0}} $；L和D为无量纲加速度；m₀g₀为无量纲化因子。

考虑到与终端能量管理端的交接，再入段轨迹设计的终端需满足给定高度、速度条件以及位置条件，即

$ \begin{array}{*{20}{c}} {h\left( {{\tau _f}} \right) = {h_f},v\left( {{\tau _f}} \right) = {v_f},}\\ {\lambda \left( {{\tau _f}} \right) = {\lambda _f},\phi \left( {{\tau _f}} \right) = {\phi _f}.} \end{array} $

(6)

本文倾侧角为唯一控制变量，倾侧角可用区间为0°~90°，倾侧角最大调节速度为8°/s。制导周期采用0.1 s。

1.3 仿真对象

为验证本论文提出算法的有效性，引用美国波音公司1998年设计的带控制翼的锥形体再入机动飞行器CAV-L作为仿真对象，总体参数如表 1所示，其他总体气动参数详见文[17]和[23]。

2 飞行器控制参数化

在初始投放阶段，大气密度稀薄，飞行器机动能力有限，因此飞行器不作轨迹调整，沿固定倾侧角飞行。整个再入制导过程分为初始下滑段和拟平衡滑翔段，初始下滑段固定倾侧角设计方法和截止条件详见文[9]，在此不再赘述。

飞行器进入拟平衡滑翔段，将有足够的气动能力进行轨迹调整，所以本文的可行轨迹预测校正算法针对的是平衡滑翔段。

2.1 复合HV走廊

飞行器高度-速度(height velocity, HV)走廊是轨迹设计中最常用的一种飞行走廊，它具有求解简单，物理含义明确的优点，且在实际轨迹设计中应用方便。飞行器复合走廊由常规HV走廊H_CUB和H_CLB、考虑末端状态的末段飞行包络H_TUB和H_TLB复合而成，具体设计过程详见文[23]，如图 1所示。复合走廊上界记为H_Ub，下界记为H_Lb。复合走廊为飞行器提供了可行的飞行空间，原则上，在走廊内部，飞行器满足滑翔条件、过程约束和终端高度、速度约束。复合走廊的构建为飞行器可行轨迹实时解算奠定了基础。

2.2 射程与高度的关系

根据式(1)，并忽略地球自转项。可以得到射程与速度的简化导数关系为

$ \frac{{{\rm{d}}{S_{\rm{e}}}}}{{{\rm{d}}v}} = \frac{{v\cos \theta /r}}{{ - D/{m_0} - g\sin \theta }}. $

(7)

在平衡滑翔假设下，可以粗略认定r≈R₀，θ≈0，将D具体表达式代入式(7)中可以得到

$ \frac{{{\rm{d}}{S_{\rm{e}}}}}{{{\rm{d}}v}} = \frac{{{m_0}}}{{ - \frac{1}{2}\rho v{C_L}{S_{{\rm{ref}}}}{R_0}}}. $

(8)

其中：C_L是关于α的函数，也就是速度的函数，对应关系是明确的；ρ=ρ₀e^(-h/h_s)，ρ是关于h的递减函数，所以dS_e/dv是关于高度h的递增函数。根据以上简化公式，可以认定，飞行高度越高，飞行器远程滑翔距离S_go越大，S_go为飞行器从当前状态起始到终端速度v_f的实际飞行射程。

2.3 控制剖面参数化

可行轨迹快速规划的本质是如何在有限时间内规划一条满足多种约束条件的可行剖面。本文在复合走廊构建的基础上，提出将复合走廊上下界加权得到高度-速度剖面的可行轨迹生成方法，表示为

$ {h_{{\rm{ref}}}} = \omega \cdot {H_{{\rm{UB}}}} + \left( {1 - \omega } \right) \cdot {H_{{\rm{LB}}}}. $

(9)

其中ω为走廊上下界的加权系数。由复合走廊的推导过程可知，设计的h_ref是关于速度的函数，且h_ref自动满足飞行滑翔能力、过程约束、终端约束。在纵向轨迹中，控制约束暂不考虑下，射程成唯一等式约束。由于ω∈[0, 1]的调整可实现高度剖面的升降，结合前文论述的高度与射程的关系，确立通过在线调整ω实现射程S_go精确满足的控制参数化方案。

3 轨迹在线迭代预测校正算法 3.1 基于Gauss积分法的射程快速预测

鉴于权重系数ω值大小可实现设计剖面h_ref的高低调整，从而影响射程，所以ω与射程之间可以建立一定的量化函数关系。

已知高度剖面，射程的推导公式为

$ \begin{array}{*{20}{c}} {{S_{{\rm{pre}}}} = \int_{{v_{{\rm{now}}}}}^{{v_f}} {\frac{{v\cos \theta /\left( {{R_0} + {h_{{\rm{ref}}}}} \right)}}{{ - D - g\sin \theta }}{\rm{d}}v} = }\\ {\int_{{v_{{\rm{now}}}}}^{{v_f}} {f\left( {{h_{{\rm{ref}}}},v} \right){\rm{d}}v} .} \end{array} $

(10)

其中：S_pre为预测射程，r=R₀+h_ref，h_Now和v_Now为当前高度、速度。结合式(1)中第1和4项可以反估算出弹道倾角：

$ \sin \theta = \frac{{ - D \cdot \left( {{\rm{d}}{h_{{\rm{ref}}}}/{\rm{d}}v} \right)}}{{{m_0}\left( {g \cdot \left( {{\rm{d}}{h_{{\rm{ref}}}}/{\rm{d}}v} \right) + v} \right)}}. $

(11)

其中dh_ref/dv≈(h_Now－h_f)/(v_Now－v_f)。

由式(10)可知，射程的估算需要射程积分，而飞行器远程滑翔的特性决定了射程积分的代价是很高的，因此本文引进Lagrange多项式对被积函数进行拟合，配点选作LGL点，它们是多项式$\left( {1 - {\tau ^2}} \right){\dot P_{K - 1}}\left( \tau \right)$的根，其中K阶Legendre多项式P_K(τ)的表达式为

$ {P_K}\left( \tau \right) = \frac{1}{{{2^K}K!}}\frac{{{\rm{d}}k}}{{{\rm{d}}{\tau ^K}}}\left[ {{{\left( {{\tau ^2} - 1} \right)}^K}} \right],K = 0,1, \cdots ,n. $

(12)

因此改进射程被积函数为

$ f\left( {{h_{{\rm{ref}}}},v} \right) \approx \sum\limits_{i = 1}^K {{{\tilde L}_i}\left( \tau \right)f\left( {{H_{{\rm{ref}}}}\left( {{\tau _i}} \right),{v_i}} \right)} . $

(13)

其中${\tilde L_i}\left( \tau \right) = \prod\limits_{j = 0, {\rm{ }}j \ne i}^K {\frac{{\tau - {\tau _j}}}{{{\tau _i} - {\tau _j}}}} $。

最终，射程预测公式可以改进为

$ {S_{{\rm{pre}}}} = \int_{{v_{{\rm{Now}}}}}^{{v_f}} {\sum\limits_{i = 1}^K {{{\tilde L}_i}\left( \tau \right)f\left( {{H_{{\rm{ref}}}}\left( {{\tau _i}} \right),{v_i}} \right){\rm{d}}v} } . $

(14)

进而，

$ {S_{{\rm{pre}}}} = \frac{{{v_f} - {v_{{\rm{Now}}}}}}{2}\sum\limits_{i = 1}^K {{w_k}f\left( {{H_{{\rm{ref}}}}\left( {{\tau _i}} \right),{v_i}} \right)} . $

(15)

其中：${w_k} = \smallint _{{\rm{ - }}1}^1{\tilde L_i}\left( \tau \right){\rm{d}}\tau $为Gauss权重，由离线设计完成。

对比式(10)和(15)可以看出，通过对积分函数f(h_ref, v)采用Lagrange多项式进行拟合，最终将积分问题转换为两个向量内积问题，从而大大提高算法运行速度。

3.2 基于递推最小二乘估计的模型在线辨识

传统的轨迹校正方法是单个周期内采用优化手段确定满足约束的控制参数，经典的如Lu等^[15]采用Gauss-Newton求解满足射程的初值倾侧角。传统的轨迹校正方法存在两点不足。首先，周期与周期之间虽然参数可以传递，但是模型信息是分割的，也就是说之前求解的信息没有充分利用起来。其次，也正是每个周期迭代开始时，对模型了解不充分，参数求解需要多次迭代，才能结合梯度或者模型信息，求解出高精度的解，这是由问题的属性决定的。

结合式(9)和(15)，可以建立走廊权重系数w与预测射程的函数关系为

$ {S_{{\rm{pre}}}} = f\left( {\omega ,{x_{{\rm{Now}}}}} \right). $

(16)

其中x_Now代表当前状态。分析对应函数可以得出以下结论：首先，函数f(ω, x_Now)随ω单调连续变化；其次，函数f(ω, x_Now)随着状态的转移在缓慢变化；再次，S_pre由于模型简化、拟合误差以及仿真截止误差等原因，实际ω预测射程S_pre值会存在噪声。

基于以上理由，需要采用一套可靠、在线运行的模型辨识算法对模型进行持续辨识与传递，并加以利用。考虑S_pre与ω为非线性关系，在当前状态下按照小扰动理论展开为

$ {S_{{\rm{pre}}}} = f\left( {\omega ,{x_{{\rm{Now}}}}} \right) \approx a \cdot \omega + b. $

(17)

即，将f(ω, x_Now)在当前状态下一阶Taylor展开，b表示当前状态下理论ω对应的射程，a表示当前状态下单位w对应的射程变化量。鉴于在实际飞行中，ω只进行小幅度变化，S_pre≈a·ω+b实际上给出了非线性函数关系S_pre=f(ω, x_Now)在当前ω值附近切线，而不是整个ω可行域的线性拟合。值得注意的是，随着状态的转移以及ω的变化，线性化模型是不断变化的，因此需要不断根据周期内产生的轨迹预测信息对线性模型进行修正和更新，以满足当前状态下的线性拟合。

考虑线性离散观测模型为

$ \mathit{\boldsymbol{h}}_i^{\rm{T}}\mathit{\boldsymbol{x}} + {v_i} = {z_i}. $

(18)

其中：增益矩阵${\mathit{\pmb{h}}_i} = \left[ \begin{array}{l} {\omega _i}\\ \;1 \end{array} \right]$和观测值为z_i=S_pre(ω_i)都是已知的；$\mathit{\pmb{x}} = \left[ \begin{array}{l} a\\ b \end{array} \right] \in {\mathbf{R} ^2}$是估计参数向量，${v_i} \in {{{\mathbf{R} }}^2}$是量测噪声，不可测量。将经过k次观测得到的数据集合起来，得到式(18)的矩阵形式，即

$ \mathit{\boldsymbol{H}}{ = _k} + {\mathit{\boldsymbol{V}}_k} = {\mathit{\boldsymbol{Z}}_k}. $

(19)

其中

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_k} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{h}}_1^{\rm{T}}}\\ \vdots \\ {\mathit{\boldsymbol{h}}_k^{\rm{T}}} \end{array}} \right] \in {R^{k \times 2}},{\mathit{\boldsymbol{V}}_k} = \left[ {\begin{array}{*{20}{c}} {{v_1}}\\ \vdots \\ {{v_k}} \end{array}} \right] \in {R^k},}\\ {{\mathit{\boldsymbol{Z}}_k} = \left[ {\begin{array}{*{20}{c}} {{z_1}}\\ \vdots \\ {{z_k}} \end{array}} \right] \in {R^k}.} \end{array} $

(20)

使得二次型指标最小的最优观测值${\hat {\mathit{\pmb{x}}}_k}$为

$ {\mathit{\boldsymbol{\hat x}}_k} = {\left( {\mathit{\boldsymbol{H}}_k^{\rm{T}}\mathit{\boldsymbol{Q}}_k^{ - 1}{\mathit{\boldsymbol{H}}_k}} \right)^{ - 1}}\mathit{\boldsymbol{H}}_k^{\rm{T}}\mathit{\boldsymbol{Q}}_k^{ - 1}{\mathit{\boldsymbol{Z}}_k}. $

(21)

其中Q_k^－1为正定矩阵，用来决定误差向量$\left( {{\mathit{\pmb{Z}}_k} - {\mathit{\pmb{H}}_k}{{\hat {\mathit{\pmb{x}}}}_k}} \right)$的分量权重。

在线参数估计中，新的信息不断得到，本文引进递推的最小二乘估计来实现参数的在线修正。假设在时刻得到新的观测值z_k+1，即

$ \mathit{\boldsymbol{h}}_{k + 1}^{\rm{T}}\mathit{\boldsymbol{x}} + {v_{k + 1}} = {z_{k + 1}}, $

(22)

记

$ {\mathit{\boldsymbol{P}}_k} = {\left( {\mathit{\boldsymbol{H}}_k^{\rm{T}}\mathit{\boldsymbol{Q}}_k^{ - 1}{\mathit{\boldsymbol{H}}_k}} \right)^{ - 1}}. $

(23)

则由式(21)可以得到

$ {\mathit{\boldsymbol{\hat x}}_k} = {\mathit{\boldsymbol{P}}_k}\mathit{\boldsymbol{H}}_k^{\rm{T}}\mathit{\boldsymbol{Q}}_k^{ - 1}{\mathit{\boldsymbol{Z}}_k}. $

(24)

设前后信息的加权矩阵为

$ \mathit{\boldsymbol{Q}}_{k + 1}^{ - 1} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}_k^{ - 1}}&{}\\ {}&{R_{k + 1}^{ - 1}} \end{array}} \right]. $

(25)

其中：Q_k^－1为已有信息向量权重，R_k+1^－1为新信息的权重，鉴于P_k与P_k+1存在前后的递推关系为:

$ {\mathit{\boldsymbol{K}}_{k + 1}} = {\mathit{\boldsymbol{P}}_k}{\mathit{\boldsymbol{h}}_{k + 1}}{\left( {{R_{k + 1}} + \mathit{\boldsymbol{h}}_{k + 1}^{\rm{T}}{\mathit{\boldsymbol{P}}_k}{\mathit{\boldsymbol{h}}_{k + 1}}} \right)^{ - 1}}, $

(26)

$ {\mathit{\boldsymbol{P}}_{k + 1}} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_{k + 1}}\mathit{\boldsymbol{h}}_{k + 1}^{\rm{T}}} \right){\mathit{\boldsymbol{P}}_k}. $

(27)

则最终可以得到${\hat {\mathit{\pmb{x}}}_{k + 1}}$的递推关系为

$ {\mathit{\boldsymbol{\hat x}}_{k + 1}} = {\mathit{\boldsymbol{\hat x}}_k} + {\mathit{\boldsymbol{K}}_{k + 1}}\left( {{z_{k + 1}} - \mathit{\boldsymbol{h}}_{k + 1}^{\rm{T}}{{\mathit{\boldsymbol{\hat x}}}_k}} \right). $

(28)

3.3 跨周期迭代算法

由节3.2可知，采用递推的最小二乘估计方法可以将飞行器累积的轨迹预测信息收集起来，并将辨识得到的模型信息传递给下一个制导周期。鉴于在每个周期算法运行时，有比较可靠的模型信息，参数迭代算法可以大大简化。

根据参数的最优估计值为${\hat {\mathit{\pmb{x}}}_{k + 1}}$，可以确定飞行器当前状态下走廊权重系数w与预测射程的函数关系为:

$ \left[ {\begin{array}{*{20}{c}} {\hat a}\\ {\hat b} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat x}}}_{k + 1}}\left( 1 \right)}\\ {{{\mathit{\boldsymbol{\hat x}}}_{k + 1}}\left( 2 \right)} \end{array}} \right], $

(29)

$ {S_{{\rm{pre}}}} \approx \hat a \cdot \omega + \hat b. $

(30)

其中$\hat a$等于${\hat {\mathit{\pmb{x}}}_{k + 1}}$的第1个元素值，代表了单位ω引起的射程的变化。采用New-Raphson方法实现ω的迭代修正，修正公式为

$ {\omega _{k + 1}} = {\omega _k} + \left( {{S_{{\rm{pre}}}}\left( {{\omega _k}} \right) - {S_{{\rm{go}}}}} \right)/\hat a. $

(31)

修正后的标准高度-速度剖面为

$ {h_{{\rm{ref}}}}\left( {k + 1} \right) = {\omega _{k + 1}} \cdot {\mathit{\boldsymbol{H}}_{{\rm{UB}}}} + \left( {1 - {\omega _{k + 1}}} \right) \cdot {\mathit{\boldsymbol{H}}_{{\rm{LB}}}}. $

(32)

最终，一套完整的基于飞行器复合HV走廊的可行轨迹跨周期快速预测校正算法设计完成，如图 2所示。新算法在以下几个方面区别于传统的预测校正方法：首先在常规HV走廊的基础上，提出了复合HV走廊的概念，基于复合HV走廊完成可行轨迹控制参数化建模，将复杂多约束问题转化为单射程等式约束求根问题；其次，为提高射程预测的速度，引进Lagrange多项式对积分函数进行拟合，将积分问题转化为显式函数计算问题；然后，基于递推的最小二乘估计方法将历史预测信息收集起来并进行挖掘，得到可靠的权重与射程对应函数。最后，采用Newton-Raphson方法完成ω的修正。在此必须补充说明一点的是本文设计的算法是跨周期运行的，也就是每个周期只需要一次射程预测即可实现轨迹的调整。

4 跟踪控制器设计

本文设计的可行剖面是高度剖面，考虑到高度剖面变化太快会引起倾侧角跟踪指令的大幅跳跃，因此在轨迹跟踪器之前对高度进行滤波，滤波算法为

$ \begin{array}{*{20}{c}} {h_{{\rm{cmd}}}^k = \left( {1 - {k_{{\rm{Pf}}}}} \right) \cdot h_{{\rm{cmd}}}^{k - 1} + {k_{{\rm{Pf}}}} \cdot h_{{\rm{ref}}}^k + }\\ {{k_{{\rm{If}}}}\int {\left( {h_{{\rm{cmd}}}^{k - 1} - h_{{\rm{ref}}}^{k - 1}} \right){\rm{d}}t} .} \end{array} $

(33)

其中：h_ref^k和h_ref^k－1为当前采样点和前一采样点对应的参考高度，h_cmd^k和h_cmd^k－1为当前采样点和前一采样点滤波后的高度值，k_Pf为滤波器比例参数，k_If为滤波器积分参数。

高度的跟踪是通过调整倾侧角实现的。首先对高度进行连续二次求导，得到

$ \ddot h = \left[ {\frac{{L\cos \sigma }}{m} + \left( {\frac{{{v^2}}}{r} - g} \right)} \right]\cos \theta - \left( {\frac{D}{m} + g\sin \theta } \right)\sin \theta . $

(34)

令Δh=h－h_cmd，配置误差Δh收敛律为

$ \Delta \ddot h + 2 \cdot {\lambda _h}\Delta \dot h + \lambda _h^2\Delta h = 0. $

(35)

其中λ_h为高度收敛系数。最终，综合式(34)和(35)可以反算出倾侧角值，实现对纵向高度和高度下降速度的精确跟踪。通过飞行过程中高度剖面的不断变化，可实现末端速度的精确控制，因此本文制导方案不需要设计专门的速度跟踪控制器。本论文中横向制导采用参考文[16]中的航向角误差走廊方案。

5 数字仿真测试

为验证本文提出算法的有效性，开展两组试验，分别测试算法的实时性、可靠性和标准环境下的任务自主能力。

5.1 算法实时性测试

此部分试验主要对比3种参数迭代方法的求解时间，3种方法分别为基于积分射程预测+Gauss-Newton方法求解(I & GN)，基于Gauss积分法+Gauss-Newton方法求解(G & GN)，以及基于Gauss积分公式+跨周期New-Raphson方法求解(G & PCNR)。

从求解时间对比来看，I & GN与G & GN求解时间的对比主要是显示Lagrange插值多项式对积分对象拟合的效果，对比显示求解速度提高约46倍。G & GN与G & PCNR的时间对比主要显示Gauss-Newton方法与本文提出的基于模型在线辨识的跨周期New-Raphson方法的对比，前后信息的综合利用可以提高算法约6.2倍。3种方法的实时性在表 2中更加直观。

5.2 标准环境下多任务仿真试验

本文随机选取4个典型再入过程，再入点分别为case1(160, 5)、case2(180, 5)、case3(175, 10)、case4(175, 0)。

图 3—6给出了4个单独任务的仿真图，从图中可以看出：首先，从制导效果来看，飞行器能实现精确地再入飞行，横侧向控制良好，且图 4给出的高度-速度剖面可以看出，飞行器平稳再入滑翔，虽然存在轨迹超过上界，短时间内不满足平衡滑翔的状况，但是软约束不会影响飞行的安全，且所有飞行轨迹高于走廊下界，说明任务再入过程都严格满足过程约束限制，此外终端高度和速度约束也较好地得到满足。其次，从算法迭代稳定性来看，算法能够迅速收敛到目标值，且随着射程需求的改变，ω能够稳定调整，从高度-速度剖面上看，不同射程需求对应不同的飞行高度剖面。再次，为消除横向机动带来的误差损失，图 5中射程系数ω随着飞行慢慢变大，这也验证了制导算法在线自主调整能力。图 6给出了算法的射程调整精度和调节时间统计。可以看出，整个再入过程，飞行器通过对过去信息的挖掘及建模，可实现非常准确的跨周期参数迭代，迭代精度和实时性都非常高。

6 结论

高超声速再入制导轨迹设计是一个复杂的多约束优化问题，本文提出一套跨周期运行的轨迹预测校正算法，并结合轨迹跟踪控制器设计，实现飞行器自主和高精度再入复合制导。

本文基于复合HV走廊将复杂轨迹设计问题转化为一元函数求根问题，在此基础上，引入Lagrange多项式拟合被积函数，实现轨迹射程的快速预测，并采用递推最小二乘估计方法实现在线信息挖掘与模型辨识，最终实现可靠求解模型的传递。基于已知求解模型，提出跨周期迭代算法，实现求解参数的跨周期、高精度迭代求解。综合以上技术手段，最终将再入制导轨迹预测校正算法控制在严格制导周期以内(MATLAB仿真下小于0.1 s)。结合轨迹跟踪和多任务仿真验证了算法的实时性和自主性。