基于解析剖面的时间协同再入制导

2019-03-29王肖郭杰唐胜景祁帅

航空学报 2019年3期

王肖，郭杰，唐胜景，祁帅

北京理工大学宇航学院，北京 100081

高超声速滑翔飞行器作为作战武器使用时，具有速度快、射程远、精度高、机动和突防能力强等特点，近年来受到世界各国的广泛关注[1-2]。与此同时，针对高超声速目标威胁，各国相继研发了THAAD、宙斯盾、C-400等防空反导武器系统，使得单个高超声速滑翔飞行器的突防能力和作战效能大大降低[3]。在此背景下，发展多高超声速滑翔飞行器协同打击技术，可大大提高对防空反导武器系统的突防概率，实现对目标的饱和攻击[4]。

以文献[5-7]为代表的基于协调变量的双层协同制导架构，以包含协调变量的协调策略为上层协调控制，以满足飞行器特点的带约束制导律为底层控制。以文献[8-10]为代表的“领弹-从弹”协同制导架构，通过使从弹跟踪领弹运动状态以实现多导弹协同制导。文献[11-12]提出了一种双阶段协同制导架构，第1阶段通过一致性算法使过渡状态一致，然后切换至比例导引实现终端协同，无需剩余时间估计。现有的协同制导大部分是基于以上3种协同架构。

然而，上述协同制导均研究的是末制导段。对于高超声速滑翔飞行器而言，其再入飞行段占整个返回段飞行时间的90%以上[13]。若不考虑再入段的时间协同，则末制导段初始状态可能相差很大，而末制导段的时间调节范围有限，难以保证最终协同效果。因此，有必要研究时间协同的多飞行器再入制导，使多个飞行器在同一时刻到达指定位置，从而为协同末制导提供良好的交班条件[14-15]。

上述针对末制导段的协同方法也难以直接运用到协同再入制导中。这些方法一般研究的是定常速度下的飞行器运动学。然而对于再入飞行，由于整个过程中空气密度变化很大，速度变化范围大，飞行器动力学变化剧烈，热流、过载、动压等过程约束严苛，仅仅考虑定常速度下的飞行器运动学难以满足再入飞行需求。此外，再入过程中“黑障”区的存在以及各飞行器间距离较大超过通讯距离，使得基于通讯的各类协同制导方法也不再适用[16]。因此，必须考虑再入飞行的特殊性，设计区别于传统协同末制导的时间协同再入制导律。

目前关于时间协同再入制导的研究较少。文献[17]基于模型预测静态规划设计了协同再入制导律，但该方法需预知各飞行器的飞行时间，从而相应地改变发射时间以实现终端时间一致，难以称为协同制导。文献[13]首次提出了一种时间可控再入制导律，并基于多飞行器时间协调信息设计了协同再入策略。该方法通过BP神经网络在线预测剩余飞行时间，但需针对特定飞行器进行大量离线训练。通过改变航向角走廊宽度增加侧向机动以改变飞行时间，但由于再入飞行器的侧向机动能力较弱，在不改变纵向动力学的情况下该方法调节时间能力较弱。该文献显示其时间调整范围仅为整个再入时间4%～5%，这就大大降低了该方法的实际意义，且较宽的航向角走廊可能导致飞行器最终错过目标。此外，其协调时间的计算还依赖于弹间通讯，不适用于再入飞行。

再入制导方法一般分为标准轨迹制导和预测校正制导两类[18]。由于飞行时间是在线状态量，可优先考虑使用预测校正制导，通过某种方法在线预测剩余飞行时间，校正飞行轨迹，从而实现时间约束以及协同飞行。

基于以上分析，本文提出一种基于高度-速度剖面的时间协同制导律。首先在高度-速度剖面内设计了由两个轨迹参数确定的参考轨迹，则剩余航程和飞行时间可表示为两个轨迹参数的函数。通过在线预测剩余飞行航程和时间并校正两个轨迹参数，实现了时间约束再入制导。在此基础上针对多飞行器协同再入任务设计了协同策略。该策略无需离线训练和弹间通讯且时间可控范围更大，更加适用于实际的再入过程。仿真结果说明了本文方法的有效性。

1 时间协同再入问题描述

1.1 运动方程

假设地球为均匀圆球，且不考虑自转，则多个高超声速滑翔飞行器的三自由度动力学方程描述为[19]

(1)

式中：i代表第i个飞行器;r为地心距;V为飞行速度;g为重力加速度;λ和φ分别为地球经度和纬度;θ为航迹角;ψ为航向角;σ为倾侧角；L和D分别为升力加速度和阻力加速度，其计算公式可表示为

(2)

式中：m为质量;S为特征面积;ρ=ρ0e-h/7 110为大气密度，ρ0=1.225 kg/m3为海平面大气密度，h=r-R0为高,R0为地球半径;CL和CD分别为升力和阻力系数，是迎角α和速度V的函数。文献[20]给出了基于非线性最小二乘辨识的再入飞行器解析气动系数模型，具有较高精度:

(3)

式中：CLi、CDi(i=0,1,2,3)为辨识出的气动系数，具体数值参见文献[20]。

1.2 再入约束

再入过程约束包括热流约束、过载约束、动压约束和准平衡滑翔条件[19]

(4)

(5)

(6)

(7)

终端约束包括终端高度约束、速度约束及经、纬度约束：

(8)

式中：tf,i为第i个飞行器的终端飞行时间；rf、Vf、λf和φf分别为相应的终端约束值。对于时间协同再入问题，还有终端飞行时间约束：

tf,1=tf,2=…=tf,n=tf

(9)

式中：tf为预设的协同飞行时间。

2 时间约束再入制导律

2.1 飞行剖面设计

再入过程中考虑到初始段热保护要求和后续航程要求，通常采用分段函数形式的迎角方案：

(10)

式中：α1为一个较大的迎角；α2为最大升阻比迎角；Va和Vb为给定的速度值。

在给定迎角方案后，即可将热流、过载、动压约束和准平衡滑翔条件转换为高度-速度平面内的再入飞行走廊：

(11)

初始下降段飞行器高度较高，气动力作用很弱，难以对轨迹进行有效控制，通常采用常值倾侧角开环制导，并当满足一定条件后转入滑翔段[21-22]。

针对滑翔段纵向制导，为充分利用走廊宽度和飞行器的再入能力，本文设计了如下的3段多项式形式的解析剖面：

(12)

式中：kij(i=1,2;j=0,1,2,3)分别为多项式系数；V1,V2∈(Vf,Vtran)为滑翔段内选定的两个速度值；(Vtran、Htran)为初始下降段与滑翔段间的过渡点；H1、H2分别为V1、V2处的选定高度且满足：

Hi=Hmax(Vi)-k1(Hmax(Vi)-Hmin(Vi))

i=1,2

(13)

其中：Hmax(Vi)、Hmin(Vi)分别为飞行走廊上对应速度Vi处的高度上边界和下边界；k1∈(0,0.98)为待设计的高度系数。当k1增大时，中段直线会更靠近走廊下界，整个H-V剖面也随之下降。

k1确定后，中段直线段即确定。另外两段均为三次多项式，确定三次多项式参考轨迹需要4组约束条件。当V1≤V≤Vtran时，考虑H-V剖面内滑翔段初始点(Vtran,Htran)及中段直线段连接点(V1,H1)处连续性和光滑性要求即可确定。

当Vf≤V

H3=Hmax(V3)-k2(Hmax(V3)-Hmin(V3))

(14)

式中：k2∈(0,0.98)为第2个待设计高度系数。

于是确定了整个H-V剖面，如图1所示。该剖面分为3段。其中第1段多项式用于将飞行器快速拉升，避免超出走廊下边界。中段直线段通过调整高度系数k1可充分利用走廊宽度，以调节飞行器航程和飞行时间，同时保证轨迹在走廊以内。第3段多项式一方面满足终端高度、速度约束，另一方面通过高度系数k2调整飞行宽度，从而调节飞行器航程和飞行时间。

将再入段的H-V剖面设计为解析多项式是较为传统的方法。但传统方法均只有一个轨迹参数以满足航程约束，相比之下本文设计的剖面有两个待设计的轨迹参数，能够更加充分地利用走廊高度和改变飞行器纵向动力学，从而同时满足射程和飞行时间约束。

图1 H-V剖面内的纵向轨迹Fig.1 Longitudinal trajectory in H-V profile

2.2 剩余飞行航程和时间预测

本文在每次轨迹预测中，同时预测剩余航程和剩余飞行时间。对于剩余航程sgo，在假设飞行轨迹近似于大圆弧下有

(15)

根据运动方程式(1)可得

(16)

则剩余航程对速度的导数为

(17)

考虑到dt=-dtgo，则由式(16)，剩余飞行时间tgo对速度的导数为

(18)

由式(2)知，阻力加速度D仅与高度、速度有关，在H-V剖面内仅与V有关。

由运动方程式(1)取高度对速度的微分有

(19)

(20)

(21)

进一步考虑在初始、终端速度确定的情况下，所设计的H-V剖面仅与两个高度系数k1、k2有关，于是式(20)、式(21)可写为

sgo=f(k1,k2)

(22)

tgo=g(k1,k2)

(23)

式中：f(·)、g(·)分别为剩余航程和飞行时间，为两个高度系数的函数。

图2和图3分别展示了以CAV-H飞行器为例，在Vf=2 000 m/s、Hf=25 km约束和迎角方案式(10)下，数值仿真得到的航程、飞行时间与高度系数k1、k2间的关系。从图中可知，随着高度系数减小，航程、飞行时间均单调增大。由式(20)、式(21)分析知，高度系数减小，则飞行轨迹越接近走廊上边界，飞行高度增大，阻力加速度减小，航程与飞行时间均增大。理论分析与仿真结果一致。

图2 航程与高度系数的关系Fig.2 Relationship between range and altitude coefficients

图3 飞行时间与高度系数的关系Fig.3 Relationship between flight time and altitude coefficients

2.3 校正算法

传统预测校正算法往往采用牛顿迭代法或割线法以单航程约束校正倾侧角。本文在每个制导周期内，以航程、时间双约束，基于H-V剖面同时对两个高度系数进行校正：

(24)

式中：s1为理想剩余航程，可由当前位置与终端位置求出；t1为理想剩余时间，可由终端时间减去已飞时间得到。上述两个变量在每个制导周期内均为已知量，则方程组式(24)为关于高度系数k1、k2的二元非线性方程组

(25)

式中：F(k1,k2)=sgo-s1,G(k1,k2)=tgo-t1。利用二元非线性方程组求根的牛顿迭代法可快速求解方程组式(25)[23]:

(26)

(27)

式中：Δ为一正小量。利用式(26)计算得到下一次迭代的k1、k2值，重复迭代，直到F<ε1且G<ε2时停止迭代，得到满足精度的k1、k2的解。ε1、ε2为允许的航程和时间误差。

实际仿真中发现，当飞行过程中存在参数扰动且其影响较大时，方程组式(25)可能得不到满足精度的解。此时，取目标函数

(28)

式中：s0为航程误差的归一化参数；t0为时间误差的归一化参数。于是将航程、时间双约束的预测校正问题，转化为使目标函数式(28)最小的参数优化问题。利用参数优化的牛顿迭代法可快速求解出使目标函数最小的高度系数k1、k2[24]:

(29)

高度系数确定之后则参考H-V剖面确定。根据运动学方程式(1)，选取控制输入为cosσ，对高度求二阶微分，得到动力学系统为

(30)

式中：a=-Dsinθ-g+V2cos2θ/r;b=Lcosθ。设计控制律跟踪参考高度Hd:

(31)

式中：ξ、ω分别为阻尼比和自然频率。则控制输入cosσ为

(32)

于是得到了倾侧角的幅值。

2.4 侧向制导

侧向制导采用经典的航向角走廊方法确定倾侧角符号。当航向角误差超过预设的误差走廊时，改变倾侧角符号使其重新回到走廊内；当航向角误差未超过误差走廊时，保持倾侧角符号不变。航向角走廊方法相比于一次或两次翻转方法，倾侧角符号翻转次数较多，但对在线参数扰动鲁棒性更好。

本文中的时间约束再入制导律的核心思路是将传统再入制导律中航程约束的单轨迹参数搜索问题扩展为航程、时间双约束的双轨迹参数搜索问题。这种思路既适用于离线的轨迹规划，也适用于在线预测校正制导；既可基于高度-速度剖面，也可基于阻力加速度-能量剖面乃至倾侧角剖面。轨迹参数既可选为两个高度系数，也可选择任意两个可同时影响轨迹高度的参数。

本文沿用了传统预测校正制导的一般方法，通过在线数值积分计算式(20)、式(21)，但必须考虑计算实时性与精度之间的矛盾。

1) 注意到式(20)、式(21)以速度为积分变量，积分终止条件为终端速度约束，于是积分步长可取为(V-Vf)/N，V为当前速度，N为积分步数，可取为100。即采用变步长积分，而积分总步数一定，当距离目标较远时，速度相差较大，对制导精度要求较低，积分步长较大；当距离目标较近时，积分步长较小，以提高预测精度。

2) 对于制导周期，即预测校正周期，在再入开始时刻，距离目标较远，对制导精度要求较低，可将制导周期设置为多个仿真步长。随着飞行器距离目标越来越近，每次预测校正的耗时逐渐减小，可逐渐减小制导周期至每一步预测校正以保证终端精度。

3) 传统数值积分多采用四阶龙格库塔法，每次积分需要4次右端函数计算。本文采用Adams预估校正法，每次积分仅需两次右端函数计算，在精度相当的情况下计算量为四阶龙格库塔法的一半，从而大大减少了计算时间。

3 多飞行器协同策略

3.1 飞行时间可调范围

在多飞行器协同再入之前，需要确定协同飞行时间，为此先确定各飞行器的飞行时间可调范围。图3通过数值仿真得到了飞行时间与高度系数k1、k2间的关系曲面，但此曲面是在无航程约束下的。对于某个确定的再入任务，其航程约束一定，此时飞行时间与高度系数间的关系退化为三维空间中的一条曲线。

图4为在Vf=2 000 m/s、Hf=25 km、航程sf=8 700 km约束下数值仿真得到的飞行时间与高度系数间的关系。从图中可知，随着高度系数k2减小、k1增大，飞行时间单调增大。

图4 航程约束下的飞行时间与高度系数关系Fig.4 Relationship between flight time and altitude coefficients under a certain range constraint

将H-V剖面的第3段多项式轨迹称为滑翔后半段，前两段多项式轨迹称为滑翔前半段。分析剩余航程与飞行时间的表达式(20)与式(21)知，两者被积函数分母相同，但剩余航程的被积函数分子为Vcosθ，而剩余时间的分子为1。这意味着在相同分母下，剩余航程受速度的影响比剩余时间更大。由于滑翔前半段速度显著大于滑翔后半段速度，可推测剩余航程受滑翔前半段影响较大，剩余时间则主要由滑翔后半段决定。因此，随着滑翔后半段高度增大，高度系数k2减小，阻力加速度减小，飞行时间增大。同时，为使总航程保持一定，前半段高度系数k1增大。理论分析与图4中的仿真结果一致。

于是，对于航程确定的再入任务，其最大飞行时间对应高度系数k2最小，即在滑翔后半段轨迹与走廊上边界相切时得到；最小飞行时间对应高度系数k2最大，即在滑翔后半段轨迹与走廊下边界相切时得到。据此思路，可在飞行之前求出两条相切轨迹，求出对应飞行时间，作为飞行时间可调范围。

图5为在Vf=2 000 m/s、Hf=25 km、sf=8 700 km约束下的对应于最大、最小飞行时间的两条轨迹，最大飞行时间与最小飞行时间相差227 s，时间可调范围约占总飞行时间的15%，大大高于文献[13]中的4%～5%。由于充分利用了纵向剖面，本文方法的时间调节范围较大。

图5 对应于最大、最小飞行时间的两条轨迹Fig.5 Two trajectories corresponding to maximum and minimum flight times

3.2 协同策略

协同再入的目的是使多飞行器在同一时刻到达指定的末制导交班点，为此需要先确定协同飞行时间。为保证协同飞行时间存在，即各飞行器时间可调范围有交集，应使再入段各飞行器的初始待飞航程大致相当。若各飞行器的初始待飞航程相差太大，无法实现协同再入。而高超声速滑翔飞行器整个飞行过程分为助推段、再入段和末制导段[13]，再入段的初始条件是助推段的终端条件。于是在发射之前，通过设计助推段程序方案以保证再入段待飞航程大致相当，从而确保协同飞行时间存在以及协同再入的可行性。于是可求取如式(33)的协同飞行时间:

(33)

式中：tmax,i为第i个飞行器的最大飞行时间；tmin,i为第i个飞行器的最小飞行时间。

再入任务开始后，初始下降段各飞行器气动力和机动能力较弱，采用常值倾侧角开环制导，仅需将运动方向对准目标，不需考虑时间约束。在所有飞行器进入滑翔段后，气动力和机动能力增强，于是各飞行器独立地以tf为时间约束执行本地时间约束再入制导律，通过在线预测剩余航程和时间，实时校正飞行轨迹，最终实现时间协同再入任务。

考虑到再入过程中“黑障”区的存在以及各飞行器间距离较大，大大超过通讯距离，该协同策略避免了集中式或分布式通讯结构的使用，更加适用于实际的再入过程。

4 仿真分析

仿真以美国通用航空器CAV-H为对象[25]。该飞行器质量为907 kg，参考面积为0.483 9 m2。迎角方案参数:α1=20°,α2=10°,Va=6 500 m/s,Vb=5 000 m/s。高度-速度剖面参数:V1=7 000 m/s,V2=5 400 m/s,V3=3 000 m/s。归一化参数s0=5 km,t0=5 s。控制律参数ξ=0.7,ω=0.1。航向角误差门限选取为10°。过程约束选取为

本文分别在标称条件和扰动条件下针对时间约束再入制导律进行仿真验证，然后进行多飞行器协同再入仿真以验证协同策略的有效性。

4.1 标称条件下多任务仿真

首先在标称条件下针对不同航程、时间约束约束下的再入任务验证时间约束再入制导律的有效性。4个不同任务的初始、终端条件设置见表1。以终端速度为仿真截止条件计算终端误差，终端误差见表2。

图6(a)为标称条件下的高度-速度曲线。4个任务下的轨迹均比较平滑，无明显振荡，基本满足准平衡滑翔条件，并能到达指定的终端速度高度。并且随着航程的增大，H-V曲线逐渐抬高，始终保持在再入走廊内，满足过程约束。图6(b) 为标称条件下的地面轨迹。表2为4个任务下的终端误差。基于解析H-V剖面的设计方法使得终端速度、高度、航程误差均较小，终端时间误差也保持在0.4 s以内，从而验证了时间约束再入制导律对再入时间的可控性，为实现多飞行器协同再入奠定了基础。

表1 再入任务Table 1 Entry mission

表2 终端误差(标称条件)Table 2 Terminal errors (standard conditions)

图6 标称条件下的高度-速度曲线和地面轨迹Fig.6 Height-velocity profile and ground tracks in standard conditions

4.2 扰动条件下时间约束仿真

为了验证时间约束再入制导律在扰动条件下的精度和鲁棒性，针对表2中任务2进行200次蒙特卡罗仿真。仿真中各初始状态偏差和参数偏差假设符合正态分布，其偏差限如表3所示。

图7(a)为扰动条件下的高度-速度曲线。扰动条件下的轨迹比较平滑，基本满足准平衡滑翔条件，且均分布在再入走廊内，保证满足过程约束。图7(b)中的终端高度-速度误差散布显示终端高度误差基本在400 m以内，速度误差在0.3 m/s 以内，具有一定的精度。图7(c)为扰动条件下的地面轨迹。图7(d)为扰动条件下的落点散布。落点沿飞行方向散布，且散布误差基本在10 km以内，满足再入制导要求。图8为终端飞行时间误差直方图。飞行时间误差分布在3.5 s以内，且分布范围比文献[13]中小，从而验证了时间约束再入制导律在扰动条件下对飞行时间的可控性。