多高超声速飞行器静态协同再入制导方法

乔 浩，李师尧，李新国

(1. 西北工业大学航天学院， 西安 710072； 2. 陕西省空天飞行器设计重点实验室， 西安 710072)

0 引 言

多飞行器的协同能够大幅提高系统整体的作战效能。按照各成员在飞行过程中是否进行通信，可将协同分为动态协同与静态协同两种。由于高超声速飞行器再入过程中高度、速度变化剧烈，且通常要经历“黑障”区域，因此仅再入初期进行信息交互的静态协同更适用于此类飞行器的协同问题。多高超声速飞行器的协同再入是现有单成员再入任务的一个扩展，不仅具有单个飞行器再入的所有特点，同时也具有其独特之处，即：在满足常规到达目标区域的同时，对各成员到达时间提出要求。

现有导弹协同方面的研究工作主要集中在常规战术导弹层面。2006年，Jeon等[1]首先提出了针对反舰导弹命中时间可控的制导方法；2010年，Shin等[2]研究了针对海洋目标与空中目标的协同防御问题，将协同问题由单纯攻击扩展到主动防御场景； 2015年，Zhang等[3]采用基于偏置比例导引法的分布式协同制导方法研究了多枚反舰导弹的同时攻击问题，并考虑了视场约束。2016年，吴云洁等[4]研究了主从编队反舰导弹在末制导阶段位置与姿态的一体化控制问题，采用动态面控制方法设计了串级型模型，设计了闭环的制导算法。与上述方法不同，Zhao等[5]提出一种轨迹曲线矫正的方法，将寻的过程的轨迹定义为特定形式的多项式，通过调节系数控制导弹满足终端落角及时间约束，具有较好的几何意义。2018年，吕腾等[6]研究了有向拓扑下无径向速度测量的多导弹协同制导问题，实现了平面内的多弹攻击问题；同年，张帅等[7]研究了带有引诱角色的有限时间协同制导问题，能够适应不同的拦截角以及初始发射条件。张聪等[8]针对反舰弹道导弹研究了制导与控制的一体化协同飞行方法，采用领-从弹的结构实现了有攻击角度约束的多弹协同攻击。2019年，赵恩娇等[9]研究了通信拓扑切换下的多飞行器协同拦截方法，解决了机动目标的协同攻击问题。

在上述协同制导问题中，大多引入了恒定速度或可控速度假设，将原动力学问题转化为简单的运动学问题进行研究。但这一简化对于高超声速飞行器并不适用。在已公开文献中，国内外针对高超声速再入飞行器的协同问题研究较少，涉及高超协同的大部分研究对象为高超声速巡航导弹，研究过程与常规导弹相类似，不能用于协同再入问题的研究。

2015年，赵启伦等[10]研究了高超声速武器与常规导弹协同攻击的可行性，形成多对一的打击态势。同年，李聪颖等[11]将再入导弹应用于反舰场景，利用SQP优化方法研究了多弹协同饱和攻击的最优弹道，证明了同一攻击过程有一定的时间可调节范围，并研究了序列再入且同时命中目标的协同攻击问题。此外，王少平等[12]也针对这类导弹的协同攻击进行了可行性层面的讨论。

2017年，Liang等[13]对双高超声速滑翔导弹的协同制导问题进行了研究，是该类型导弹协同方法研究在国内外公开文献中的首次。同年，Chu等[14]将高超声速滑翔导弹的协同再入问题转化为优化问题，采用改进的模型预测静态规划方法实现多枚导弹同时命中目标。2018年，Wang等[15]将航向偏差角定义为速度前置角，并以之作为控制量，在视线坐标系内研究了多枚再入滑翔导弹的协同制导与控制问题。2019年，Yu等[16]提出一种解析的多弹协同再入制导方法，降低了总飞行时间的预测与校正难度。但由于该解析方法建立在简化的气动力模型之上，要求飞行器具有特定形式的升阻比(文中采用常值与二次型升阻比进行研究)，对于以无动力滑翔为主要飞行方式的高超声速再入飞行器并不适用。同年，Li等[17]采用阻力加速度剖面与预测校正方法相结合的思路，研究了多高超声速滑翔导弹时间协同再入制导问题，取得了较好的效果。但由于使用了大圆弧假设，需要在每个制导周期对参考阻力加速度剖面进行更新，计算量较大，增加了算法的复杂性。

基于参考轨迹的制导方法是再入问题中形式简单、求解迅速的一类方法。Saraf等[18]提出的衍化阻力加速度方法具有较好的轨迹长度处理能力。沈振等[19]在再入走廊内采用两段折线与走廊边界设计参考阻力加速度剖面，提高了标准轨迹制导方法的实时性。王涛等[20]在阻力加速度方法中引入预测校正思想，由走廊上、下边界插值获得标准轨迹，并在三维空间对航向进行校正，提升了再入飞行能力。

本文提出一种基于参考轨迹的多高超声速飞行器静态协同再入制导方法。采用真实轨迹长度进行设计，避免了扰动环境下的参考轨迹不足；依据轨迹长度与飞行时间的对应关系，以公共轨迹长度作为参考轨迹设计时的协调变量；之后采用时域信息转换的方法，将协同逻辑转变为容易实现的状态收敛问题。最后给出完整的协同再入制导方法。

1 再入动力学模型

以总能量E为自变量，假设地球为均质圆球，并忽略地球自转，对应的动力学模型为：

(1)

式中：下标i为成员编号；V为速度大小；r为地心距；R0为地球半径；θ,φ分别为经、纬度；γ、ψ分别为航迹倾角与航向角；Rflight为轨迹长度；t为飞行时间；g为地球引力加速度；σ为倾侧角；L与D分别为升力、阻力加速度。总能量E定义为：

(2)

式中：μ为地球引力常量。在实际求解中将其归一化，定义初始能量为0，终端能量为1，任意时刻能量表示为(E-E0)/(Ef-E0)。路径约束主要有：

(3)

对于升阻比较大的再入飞行器，通常为改善轨迹质量，还需要考虑平衡滑翔约束，取平衡滑翔约束中的倾侧角σ=0°，有平衡滑翔条件为：

L-g+V2/r=0

(4)

此外，还要考虑终端高度hf,速度Vf,剩余航程Stogo,f以及航向偏差Δψf等的约束。

2 轨迹长度与飞行时间

现有的轨迹设计与制导方法在解决轨迹长度问题时往往基于大圆弧假设，使得轨迹长度需要在制导环节进行更新迭代。考虑到再入任务一旦确定，飞行器与目标之间对应的初始纵程Rl随即确定，因此将真实轨迹长度定义为：

Rflight=PmRl

(5)

式中：Pm定义为机动系数，其最大、最小值对应了最长、最短实际轨迹长度问题。可采用伪谱法进行优化确定[21]。以CAV-H为研究对象，采用表1数据作为标准再入的初始与终端状态。

表1 初始和终端要求Table 1 Initial and terminal requirements

对应有Pm∈[1,1.1373]。本文在进行参考轨迹设计时采用与衍化阻力加速度方法类似的设计思路。由于轨迹长度根据飞行任务直接给定，因此不需要进行轨迹长度的迭代。参考阻力加速度剖面的跟踪器采用二阶误差动力学模型进行设计[20]。

选取不同的机动系数进行轨迹设计，统计各机动系数对应的时间，获得二者对应关系如图1所示。可见，再入飞行时间在一定范围内与机动系数近似成线性关系。即机动系数越大，再入的总飞行时间越长。由二者之间的拟合曲线可以大致预测某一机动系数条件下的总飞行时间，为针对舰船等时间敏感目标的再入问题提供了设计依据。

图1 机动系数与总飞行时间Fig.1 Maneuver coefficient-total flight time

需要指出的是，在设计参考轨迹时同一个机动系数对应的阻力加速度剖面形式并不唯一，即同一个轨迹长度可以对应不同的参考阻力加速度剖面。在仿真中发现，即使轨迹长度相同，不同形式的阻力加速度剖面也会对应不同的飞行时间。轨迹长度与阻力加速度剖面不一一对应这一问题从侧面反映了再入过程的复杂性。因此需要研究阻力加速度剖面形式对总飞行时间的影响。令Pm=1.08，分别取靠近再入走廊上边界与下边界两种形式，对应阻力加速度剖面与飞行时间如图2和图3所示。

图2 指定机动系数两种不同阻力加速度剖面Fig.2 Different drag profiles with determined Pm

图3 不同剖面对应的飞行时间Fig.3 Time profiles with different drag profiles

可见，在给定机动系数之后，满足轨迹长度的参考阻力加速度剖面不唯一。其中慢速轨迹对应的阻力加速度剖面在初始时相对较高，但在能量约0.4之后，对应的阻力加速度明显更小。时间的演化曲线表明，在飞行的大部分能量区域内，接近准平衡滑翔边界的阻力加速度剖面对应的时间变化更慢，仅在能量约0.8的位置时间变化开始加快。对应的高度变化以及三维曲线如图4和图5所示。

图4 高度与能量关系Fig.4 Altitude vs. normalized energy

图5 快速与慢速到达三维轨迹Fig.5 Trajectories with different arrival time

图4中能量0.4对应的高度约46 km，在此之前慢速轨迹对应的下降速率较快，飞行器快速下降至大气较为稠密的高度，随之以更接近准平衡滑翔的轨迹飞行，轨迹高度更高，对应的阻力更小，滑翔时间也更长。由三维轨迹可以看出，末段慢速轨迹对应的轨迹更陡，相对俯冲更快，与能量0.8之后的时间演化曲线一致。

综上可知，机动系数与轨迹长度之间的对应关系是唯一的，轨迹长度与总飞行时间之间对应关系的唯一性是有条件的，当参考阻力加速度剖面的设计规则给定后，轨迹长度与飞行时间一一对应，且近似成正比例关系。当参考阻力加速度剖面的形式不限制时，即使给定轨迹长度，总飞行时间也具有一定的可调范围，其中慢速飞行对应的弹道形式为：初始快速下降进入稠密大气层，并在该区间维持较长时间飞行，末段弹道快速下压。

3 基于参考轨迹的再入协同框架

考虑一种协同攻击任务：多枚导弹自相距数百公里的阵位发射，要求同时到达指定区域。假设各成员型号相同，助推段末端或天基释放对应的初始能量状态相同，即初始速度、高度相同。令各成员初始航迹倾角一致，初始航向瞄准目标。利用机动系数与轨迹长度一一对应，而轨迹长度在剖面形式给定的情况下与到达时间一一对应的特点，对多阵位发射问题设计如图6所示的双层再入协同框架。

图6 基于参考轨迹的双层协同框架Fig.6 Cooperative framework based on reference trajectory

其中协调层用于在任务之前收集各成员的阵位信息，根据各成员的机动性能，确定各自可选轨迹长度的范围，通过协调单元进行比对，确定公共轨迹长度。由于不同成员具有的初始速度、高度一致，终端约束一致，当公共轨迹长度确定后即可进行公共阻力加速度剖面的设计，之后将公共阻力加速度剖面作为协调信息分发给各成员。执行层在获得分发信息之后，各成员求解本地轨迹曲率子问题，获得各自参考轨迹。结合不同的轨迹跟踪方法即构成基于参考轨迹的协同再入制导方法。此外，该协同框架对初始阵位同样有一定要求，即各成员飞行器轨迹长度区间的交集需要非空。

首先确定各飞行器对应的初始纵程Rl,i，之后利用该信息采用优化方法确定各阵位对应的机动系数区间。由于机动系数与初始纵程有关，当初始纵程不同时无法直接进行比对(相当于度量标准不同)，因此实际中直接利用轨迹长度区间进行比对，即：Pm,i×Rl,i。结合轨迹长度与飞行时间在有限条件下一一对应的特点，协同轨迹设计转化为依据各成员阵位信息，求取符合各成员机动性能的公共轨迹长度。为方便表述，此处以Rf代表成员飞行器的轨迹长度，各阵位轨迹长度区间为[Rfi,min,Rfi,max]，对应的公共轨迹长度Rflight须满足：

(6)

可见，公共轨迹长度受各成员最长与最短轨迹的限制。在确定公共轨迹长度之后，即以该轨迹长度作为公共参考阻力加速度剖面的设计依据，由协调单元完成轨迹长度子问题，并进行公共参考阻力加速度剖面的分发；轨迹曲率子问题由各成员本地解算。相比独立再入各成员的计算量大幅减少。

轨迹曲率问题采用三段式动态航向偏差走廊的方式进行求解[22]。航向偏差走廊的迭代更新方法为：

(7)

4 用于协同再入的参考轨迹跟踪方法

在完成多成员协同再入的标准轨迹设计之后，各成员还需要对标准轨迹进行跟踪以完成整个再入协同制导。由于各成员具有公共的纵向参考轨迹，因此可采用二阶跟踪误差模型设计统一的轨迹跟踪器，由各成员飞行器各自在线根据与参考状态的偏差进行纵向轨迹跟踪。轨迹曲率子问题采用动态航向偏差走廊在轨迹跟踪过程中进行航向控制。

本质上讲，阻力加速度跟踪器是一种通过跟踪参考阻力加速度剖面，实现高度、速度匹配的设计方式，而速度与高度有多种组合可以满足阻力加速度的需求，因此实际阻力加速度关于参考阻力加速度剖面最终的收敛，并不能说明轨迹长度及期望飞行时间与设计值一致。对Pm=1.08在表4所示的扰动环境下进行500次蒙特卡洛打靶，对应的阻力加速度剖面跟踪效果如图7所示。

图7 扰动环境下阻力加速度跟踪曲线Fig.7 Tracking curves of drag profile with disturbance

可以看出，不同扰动条件下阻力加速度跟踪曲线的过渡方式、过渡快慢存在差异，这就导致了实际轨迹长度、总飞行时间必然产生较大的偏差。直接跟踪参考阻力加速度剖面对轨迹长度、总飞行时间没有直接的跟踪效果的这一特点，限制了跟踪阻力加速度剖面方法在协同任务中的使用，因此本文采用另一种轨迹跟踪思路。在协同任务中，各成员依据协调信息(公共轨迹长度)确定公共参考阻力加速度剖面之后，可以结合本地轨迹曲率问题完成参考轨迹的全状态生成。这些状态包含了丰富的轨迹信息，其中高度、速度、轨迹长度(与纵程不同)、航迹倾角为纵向参数，航向偏差、经度、纬度为侧向参数，飞行时间依据公式(1)确定。由于各成员对应的纵向参考轨迹是相同的，可以考虑使用状态跟踪器对纵向状态进行跟踪，轨迹的曲率问题由动态航向偏差走廊进行控制。在由参考轨迹获得的丰富状态信息中，时间是一个确定的序列，参考轨迹一旦确定，飞行过程中的时间变化是相同的，即由协调变量确定的参考轨迹给出了一个协调一致的期望终端时间。因此可以对原能量域中的状态信息关于飞行时间进行插值，获得时域内的全状态信息，此时协同轨迹跟踪问题转化为：各成员采用相应的轨迹跟踪器，使最终到达期望末端时刻的各飞行状态符合误差要求，协同的效能评估可采用到达指定时间各成员的状态偏差进行评定。

作为状态跟踪器中较为成熟的方法，线性二次调节器(LQR)最早被应用于线性系统的最优控制问题中，其增益的求解形式简单，对不同任务有较好的适应性，因此被很多学者应用于标准轨迹制导方法的设计[18]。需要注意的是，LQR方法是一种状态跟踪器，需要以状态信息的形式给出参考轨迹。由参考阻力加速度剖面解算出状态参数既可通过对阻力加速度剖面求导获得，也可通过沿参考阻力加速度剖面进行标称跟踪获得。前者需要采用平均高度等假设，相对而言后者获得的状态信息更为精确，求解过程也更为简便。因此本文在使用LQR方法进行轨迹跟踪时采用沿公共参考阻力加速度剖面进行标称跟踪的方法获得参考轨迹的全状态信息。基于参考轨迹的多阵位协同再入制导流程如图8所示。

图8 协同再入制导流程Fig.8 Flow chart of cooperative reentry guidance

5 仿真校验

5.1 多阵位信息协同

假设协同再入任务八个备选阵位如表2所示。

表2 各成员初阵位信息Table 2 Initial states of different members

采用优化的方法求解各成员初始阵位对应的最长、最短轨迹以及轨迹长度区间如图9,图10所示。可见，最短轨迹接近于沿初始弹目连线，但在一些阵位要稍大，这主要是由于在一些距离目标位置较近的阵位，最短轨迹受到过程约束以及倾侧角演化速率的限制，需要进行一定的侧向飞行以满足热流、过载与动压的约束。

图9 不同初始阵位对应的最长、最短轨迹Fig.9 The longest and shortest trajectories of each member

图10 不同阵位最大、最小轨迹长度统计Fig.10 Statistics of maximum and minimum trajectories length

可见，依据式(7)该阵型对应的公共轨迹长度区间为空。主要原因是M3距离目标过远，需要相对其余成员提前发射。因此除M3外7个成员参与此次协同，最长、最短轨迹分别为7325.8 km、 6707.5 km。考虑到轨迹长度在下边界附近会造成倾侧角反转次数过多，在无特殊要求时应与下边界保持一定范围，此处取Rflight=6850 km为仿真校验的公共轨迹长度。

5.2 公共参考轨迹设计与标称跟踪

确定公共轨迹长度之后，即可确定满足不同阵位需求的公共参考阻力加速度剖面，如图11所示。

图11 扰动环境下阻力加速度跟踪曲线Fig.11 Tracking curves of drag profile with disturbance

各成员相应的速度-高度、航迹倾角变化曲线如图12和图13所示。

图12 参考高度-速度变化曲线Fig.12 Reference altitude-velocity

图13 参考航迹倾角变化曲线Fig.13 Reference flight path angle

由协调信息确定的各成员参考总飞行时间为1410 s。当纵向参数与终端高度、速度约束确定后，再入飞行的总时间也是确定的。因此当各成员遵守相同的协调信息时，各成员在标称环境下的时间演化曲线也是一致的，如图14所示。

图14 时间演化曲线Fig.14 Flight time

图15 各成员倾侧角指令Fig.15 Bank angle of each member

各成员对应的实际倾侧角指令如图15所示。可以看出，各成员在相应能量点处的倾侧角幅值一致，仅倾侧角的符号不同，这一不同是由不同阵位对航向的不同需求造成的，动态的航向偏差走廊通过调整倾侧角符号的反转位置以及反转次数，使自不同阵位出发的各个成员最终均到达目标区域。仿真结果如图16至图18所示。

对各成员的终端状态进行统计，如表3所示。可见，标称环境下的跟踪结果符合终端状态偏差的要求。

表3 标称状态下各成员终端偏差Table 3 Terminal deviation of different members in nominal state

图16 各成员地面轨迹投影Fig.16 Ground track of each member

图17 各成员航向偏差Fig.17 Heading error of each member

图18 各成员三维轨迹Fig.18 Three dimensional trajectory of each member

5.3 模拟扰动环境下参考轨迹跟踪

对于侧向运动不大的飞行任务，一般只需采用固定的航向偏差走廊即可，但对于本文的多飞行器协同任务，不同阵位对应的侧向运动范围很大，因此各成员在实际飞行过程中应考虑对飞行方向的控制，以满足末端剩余航程与航向偏差的要求。以根据M2成员参考轨迹获得的增益矩阵作为公共增益，各成员飞行中依据本地状态与参考状态的偏差结合预置的增益进行制导指令更新。在模拟扰动环境下考虑侧向制导逻辑，对参与协同任务的7个成员各自进行200次三维蒙特卡洛打靶。模拟扰动环境的偏差分布类型与参数如表4所示。

各成员终端状态以及终端状态偏差分布的统计打靶结果如图19所示。

由于最终的协同效果由所有成员到达截止时间的状态偏差进行评价，此处将所有成员的统计打靶结果视作同一个待统计过程，相应的偏差及相关参数的统计结果如表5所示。

表4 模拟扰动环境Table 4 Disturbance environment for simulation

图19 多成员蒙特卡洛仿真结果Fig.19 Multi-member Monte Carlo simulation results

表5 多成员蒙特卡洛仿真结果统计Table 5 Multi-member Monte Carlo simulation results statistics

在纵向轨迹确定的情况下，各成员终端的剩余航程与航向偏差主要受侧向制导逻辑的影响。动态的航向偏差走廊本质上与直接搜索最佳的倾侧角符号反转点等效，都是通过调节倾侧角的反转时机进行飞行方向的调整。各成员终端剩余航程偏差与航向偏差的统计结果说明了侧向制导逻辑的有效性。

各成员三维蒙特卡洛打靶对应的地面轨迹与三维轨迹分别如图20和图21所示。

图20 不同阵位地面轨迹Fig.20 Ground tracks for different initial positions

图21 三维轨迹Fig.21 Three dimensional trajectories

不同阵位对应的终端位置如图22所示。可以看出，各成员最终都落在了距目标一定距离的指定区域内，轨迹的终端位置较为集中。

图22 不同初始阵位对应的终端位置分布Fig.22 Terminal positions for different initial positions

由于再入过程中的过程约束主要与纵向的轨迹参数有关，而多阵位协同问题中各成员具有相同的纵向轨迹参数，因此在相同的扰动环境下，各成员的过程约束演化过程一致。以M2成员为例进行说明，500次蒙特卡洛仿真对应的过程约束如图23所示。可见，在给定的模拟扰动环境下，对应的热流、过载与动压满足约束条件，确保了各成员在协同再入过程中，均处在满足再入约束的飞行状态内。

图23 扰动环境下的热流、动压与过载约束Fig.23 Path constraints in disturbed environments

6 结 论

针对多再入飞行器的多阵位协同问题，本文给出一种基于公共参考轨迹的静态协同再入制导方法。研究表明：采用实际轨迹长度替代传统的大圆弧假设，能够更好描述飞行过程中的真实轨迹；机动系数在量化真实轨迹长度的同时，能够建立实际轨迹长度与总飞行时间之间的联系，且二者表现出一种有条件的一一对应关系，当阻力加速度剖面的形式确定后，确定的轨迹长度对应唯一的飞行时间序列，这一点为多成员的协同再入提供了切入点；基于参考轨迹的双层再入协同框架能够有效完成不同阵位多成员的协同再入问题，且具有形式简单，流程清晰的特点；提出的时域信息提取与协同逻辑转换策略可以将原协同再入问题转换为简单的截止时间状态收敛问题，对静态协同问题具有较好的适应性；以LQR方法为基础的状态跟踪器能够较好适应多阵位成员的参考轨迹跟踪问题。统计打靶的结果表明，本文所提方法能够在保证不同阵位飞行器同时到达目标区域的同时，使终端状态偏差符合需求，且具有一定的鲁棒性。