多禁飞区高超声速滑翔飞行器再入机动制导

高 杨,蔡光斌,2,张胜修,徐 慧

(1.火箭军工程大学 导弹工程学院, 西安 710025; 2.西北工业大学 航天学院, 西安 710072)

近年来，高超声速滑翔飞行器的再入机动制导研究逐渐成为世界各大国航天工业发展的一个热点[1-5]。高超声速滑翔飞行器再入飞行必须满足复杂的约束条件，其中包括传统的热流率、动压、过载约束和终端约束，也包括自然、军事等因素造成的复杂禁飞区约束[4]。飞行器的禁飞区是一种路径约束，在飞行器最优轨迹的解算中，如果禁飞区的数量越多或者模型比较复杂，那么飞行器在轨迹优化的路径约束就越多，优化轨迹的解算难度就越大。现今研究禁飞区突防的成果较少，且主要集中在离线求取最优解，其中，伪谱法是数值求解方法的代表，采用不同的配点和多项式插值方法达到全局优化的目的，这些方法无法应用到在线的机动制导中，尤其是在禁飞区信息较复杂或者信息缺失的情况下[6]。因而，预测反馈的制导机制在飞行中作用明显[7]，一些基于误差校正的神经网络与系统辨识方法被广泛应用[8-11]，然而这些研究中对于未知形状的禁飞区约束条件研究并不成熟，当预测数据因为禁飞区等因素急剧变化时，会导致反馈误差变大而导致制导错误。部分学者在研究昆虫生物学触角反馈机制的基础上，使用较少的“触角”探测前方未知区域并进行反馈，从而制定制导方案[12-14]。文献[12]在移动机器人路径规划研究中，首先提出了基于触角的方法，其中触角代表了机器人可能的运动方向。文献[14]使飞行器连续产生左右双触角对前方进行探测，通过跟踪标准轨迹并使用侧滑角正负反转策略规避不同形状的禁飞区。

鉴于现有研究结果中存在问题，本文改进了轨迹预测时触角的设计与侧滑角的瞬变策略，引入了航向角误差限制条件和延时计数开关，计算了规避策略的优先级，轨迹能够顺利规避禁飞区并到达终点。通过鲁棒性仿真验证了该策略预测过程计算时间足够小，满足机动制导策略的要求。

1 再入模型构建

高超声速滑翔飞行器的动力学模型以及再入飞行过程中的约束条件模型是高超声速滑翔飞行器机动制导过程的基础，下面对这些模型与条件分别进行定义。

1.1 高超声速滑翔飞行器再入运动学建模

针对高超声速滑翔飞行器在再入飞行阶段的机动制导过程，如图1所示，在地心赤道旋转坐标系内，忽略地球自转角的影响，建立高超声速滑翔飞行器的动力学方程[20]：

其中，r是高超声速滑翔飞行器地心距，V是飞行器的地球相对速度，ψ与γ是飞行器的航向角与航迹角，θ和φ分别是飞行器所处的经度纬度，m是飞行器的质量，g是重力加速度，σ是飞行器的侧滑角。D=ρV2SrefCD/2与L=ρV2SrefCL/2是飞行器在飞行过程中气动阻力与升力，其中ρ是空气密度，Sref是飞行器的参考横截面积，CL与CD分别是与飞行器攻角有关的空气动力学参数。不同的飞行器由于设计的不同，导致升力阻力的空气动力学参数不同以及升阻比的不同。

图1 地心赤道旋转坐标系下飞行器动力学方程建模

1.2 常规约束条件建模

在临近空间的长距离、无动力飞行过程中，高超声速滑翔飞行器会受到多种约束的影响，这些因素都会影响制导的结果。本文横向制导策略在预测的过程中，每条预测“触角”上预测时的每一个积分步长端点为一个飞行器约束检测点，预测计算时每到一个检测点都要检测是否遵循常规约束条件：

一是起始状态与终止状态的端点约束，它是对飞行器起始与终止位置、速度和姿态的约束，也是飞行器不同状态转换的交接班条件，起始条件通常为给定的任务起始状态[19]，令X=[r,θ,φ,V,γ,ψ]为飞行器状态矩阵，εX为较小的常值矩阵，终端约束模型建立如下：

|X(tf)-Xf|≤εX

二是飞行器的常规路径约束，分别是飞行器的热流率约束、过载约束和动压约束[10]，建立模型分别为

Q=KQρ0.5V3.15≤Qmax

三是航向角约束。在高超声速滑翔飞行器因禁飞区等限制原因机动飞行的过程中，为了保证飞行器转弯规避禁飞区后仍能到达目的地，需要对飞行过程中的航向角进行约束。这种约束使得飞行器能够有足够的距离通过对侧滑角和攻角的控制，把飞行器制导到终点位置。首先将航向角约束定义为在轨迹检测点上的飞行器实时航向与目标地点在经纬度平面的夹角。图2中先求出要研究的飞行器的最大转弯轨迹，在该轨迹上求出各检测点Pi(i=1,2,…,m)的速度方向与终点连线在经纬度平面上的夹角，为速度航向角误差限制Ψi(i=1,2,…,m)，即

Ψi=ψi-arctan((θi-θf)/(φi-ff))

式中，[θi,φi,ψi]为检测点Pi位置的经纬度坐标和航向角，[θf,φf]是终点的经纬度坐标。

图3中是两种飞行器Shuttle与Common Aero Vehicle-HPMARV(CAV_H)根据最大转弯轨迹计算的航向角约束以及在文献[4，14]中按照经验设定的对CAV-H飞行器的航向角约束。

图2 最大转弯轨迹与检测点上的航向角

图3 速度-航向角约束变化

1.3 多禁飞区约束建模

在飞行器实际飞行中禁飞区的形状是多样的，用解析式来表示很难，而用多个圆柱形禁飞区组合表示是比较方便的做法。假设这些禁飞区都是如图4所示的圆柱形的区域，在俯视视角下如图5所示，高超声速滑翔飞行器面对的多个小型禁飞区的复杂约束情形以及各禁飞区与飞行器的距离。

图4 圆柱体禁飞区示意图

图5 高超声速滑翔飞行器多禁飞区约束情形

将圆柱形禁飞区约束简化为在经纬度上的二维约束，针对n个圆柱形禁飞区中的第i(i=1,2,…,n)个禁飞区，令禁飞区的中心经纬度坐标为(θi,φi)，半径为Ri，其在经纬度上的禁飞区边界表示为

飞行器与各禁飞区中心之间的距离定义为{d1,d2,…,dn}，当di≤Ri(i=1,2,…,n)时，飞行器被判定进入第i个禁飞区，飞行任务失败。在未知禁飞区具体形状的情况下，飞行器只能探测到禁飞区的边界。

2 基于“触角”预测的侧滑角瞬变机动制导策略

2.1 基于线性二次调节器的纵向制导方法

在机动制导过程中，制导算法分为纵向和横向两个制导方法，纵向平面跟踪方法通过线性二次调节器(Linear Quadratic Regulator,LQR)实现[17-19]，即

其中，Δα和Δσ为攻角和侧滑角的调整量，[Δr,ΔV,Δγ]为实际轨迹与标准轨迹的地心距、速度和航迹角差值，因为是纵向平面的轨迹跟踪，所以状态量上只选择了这3个关系密切的变量。

横向制导方法是一种侧滑角符号转换方法，它与纵向制导策略同时进行，横向制导策略的周期一般大于纵向制导策略周期，使侧滑角在飞行器纵向制导的跟踪过程中保持不变，即

ΔTLat=kΔT

其中，ΔTLat与ΔT分别是横向和纵向制导策略的计算周期，k取整数，k值越大，横向机动策略计算的次数就越少，总机动策略计算的平均时间就越小，k值越小，对禁飞区的探测效果越好。本文中两个计算周期都取1 s。

2.2 基于动力学方程的三“触角”探测模型

在对高超声速滑翔飞行器机动制导的过程中，只有最开始解算出的标准轨迹是已知的，而前方飞行通道的实时禁飞状况由于不能提前探知，需要飞行器采取预测校正的制导方法进行探测。本文采取的基于“触手”的方法需要对前方发出3个触手来预测轨迹的飞行状态。这种预测轨迹是由计算机根据飞行器的动力学方程积分产生的，其中积分的步长越小，积分产生的预测路径越符合实际轨迹。而高超声速滑翔飞行器动力学方程较为复杂，预测轨迹的积分速度较慢，因而产生的预测轨迹越少越好。在文献[14]中，每一个预测周期产生两条“触角”轨迹线，极大地缩减了计算量，但是缺少了直线前进的“触角”预测轨迹线，导致在可以无机动地直线飞行时，采取了左右摆动的复杂制导策略，侧滑角正负瞬变对控制系统造成的压力也较大。

本文采取的是前向“触角”与左右两侧“触角”结合的三“触角”预测方式，如图6，当飞行器在前向触角遇到禁飞区时，即飞行段OA段，采取三触角与双触角的反馈结果是基本相同的，飞行器只能在两侧的触角中选择其中的较优路线进行机动。当飞行器在前向触角没有遇到禁飞区而停止时，即飞行段AB段，飞行器需要判断前向触角和另外两侧触角的反馈信息，文献[14]中提到的方法，在类似A点的飞行状态时，由于左右触角给出的反馈信息的实时变化，飞行器会进行左右摇摆的机动，直到因为航向角约束的逐渐减小而最后使轨迹逼近终点；而采取类似于B点的三触角探测策略，飞行器的前向触角在AB段没有禁飞区约束，因而触角的终点更贴近终点，飞行器在经过短暂的左右机动调整后可以沿直线飞行到终点。左右两侧触角预测轨迹的侧滑角符号一正一负，数值上使用飞行器的最大侧滑角度数。如果使用度数较小的侧滑角进行两侧轨迹的预测，大多时候会因为机动角度过小而无法使飞行器飞跃禁飞区；如果在原有三条“触角”的基础上，再加两条小角度的预测线，会增加计算的时间，为横向制导策略的选择增加困难，同时因为高超声速滑翔飞行器高灵敏的特性而增加控制的难度。

图6 三“触角”预测示意图

飞行器每一个预测周期的三个“触角”的起点是相同的，它们都以飞行器当前时刻仿真积分一个时间步长的状态为预测起点，当预测结束并且获得相应的控制指令后，飞行器也刚刚飞至下一时刻，由于积分时间步长固定且较小，预测起点与实际状态误差较小。三个“触角”的终止条件如下：

C1：积分时间超过了理论上飞行的最大任务时间，本文采取的最大任务时间是3 000 s；

C2：“触角”的末端触碰到了禁飞区的边缘或者因为积分步长的原因进入了禁飞区；

C3：预测轨迹线上的航向角超出了该速度所对应的最大航向角约束；

C4：预测轨迹顺利到达终点位置。

其中，前3个条件属于禁忌类条件，飞行器应避免飞行到这三类“触角”的终点，第4个终止条件属于允许类条件，飞行器应优先选择这样终止的“触手”制导飞行。其中的禁飞区终止条件是禁忌条件中禁忌程度最高的，而时间终止条件的成立几率较小，由此得到终止条件所对应的“触角”选择优先级Ka，如表1所示。

表1 终止条件优先级

3条轨迹预测只能以这4种情况结束，并将终止时的状态和终止原因反馈给横向制导方法，由横向制导方法选择三“触角”的其中一条进行机动制导。

2.3 基于侧滑角瞬变的横向制导方法

文献[4]和文献[14]等提到的侧滑角逆变策略是一种提供侧滑角符号的方法，侧滑角的绝对值是基本不变的，其中文献[14]在这一基础上根据LQR方法其侧滑角的数值进行了微调，但都属于侧滑角逆变策略，即在机动时将飞行器的侧滑角变号，改变飞行器的空气动力方向。本文将侧滑角归零这一选择加入到侧滑角的瞬变策略中，使得飞行器能够在不需要机动的时候选择前向飞行，或者在侧滑角进行正负转换时，作为稳定控制系统的过渡状态，因此飞行器的侧滑角度数有3种可能，即

σ=sgn(σ)·σmax,sgn(σ)∈{-1,0,1}

其中，σmax是侧滑角最大值，左右触角分别对应符号函数sgn(σ)的-1与1。符号函数sgn(σ)的赋值由触角反馈的总优先级K决定。总优先级K由终止条件优先级Ka和触角终止点与终点的距离优先级Kb共同决定。因此，每一个触角的总优先级如下式所示：

Ki=-Ka·Kb,i=Left,Mid,Right

其中，Kb根据触角终止位置的经纬度坐标与飞行路径终点经纬度坐标的距离排序，距离越小，该触角的终止位置越靠近终点，该触角越应该被选择，即根据这一距离的远中近3个等级，Kb分别赋值为{1.1,1.0,0.9}。式中Ki的优先级越大，代表该触角更应该被选择作为下一个时刻的制导指令，即：

为了避免飞行器侧滑角瞬变角度过大或者因为优先级的振荡变化造成的侧滑角振荡变化而控制系统不能完成这样的变化，加入一种延时锁定开关。如图 7所示，由侧滑角符号函数计算出来的结果进入延时锁定开关后，如果开关关闭，延时计数开关开始计数，控制系统中的侧滑角不变；如果开关打开，则重新判断。

图7 侧滑角计算流程框图

当前时刻为j时刻，预测周期为ΔTLat，则控制系统获得的侧滑角σ′由下式计算得到：

其中，σj、σj-ΔTLat与σj+ΔTLat代表j、j-ΔTLat与j+ΔTLat三个时刻的侧滑角符号函数输入值，延时开关的延时时间最小值根据飞行器自身特性来决定，本文中的延时时间统一取10 s。在延时过程中，统计式(17)中出现的侧滑角同类振荡变化次数，振荡次数超过一定限制，则打开延时开关，并将σj+ΔTLat输入控制系统。

3 仿真算例与分析

本节采用CAV-H和Shuttle的高超声速滑翔飞行器模型在多禁飞区条件下验证基于三触手的侧滑角瞬变机动制导策略。飞行器Shuttle和CAV-H的模型数据分别来源于文献[19]和文献[20]，表2是两个飞行器标准轨迹的飞行数据[14,19]：

根据表2初始数据获得两飞行器的初始标准轨迹高度曲线，如图8所示。本文两飞行器的标准轨迹中侧滑角设为0，攻角的变化幅度较小，导致飞行时间较长。

3.1 未知多禁飞区情况下的机动制导仿真

为检测机动制导算法的有效性，设置了11个半径为1°的小型圆柱禁飞区，它们的中心经纬度坐标如表3所示，组成了复杂的多禁飞区模型。

表2 飞行器飞行数据

图8 标准轨迹高度曲线

序号经度纬度序号经度纬度130-3765-1.523038651.5350-3980345001080055031180-36400

根据表3的禁飞区信息，按照图9的算法流程采取两种方法求解并比较。一种是利用GPOPS II工具包求解，该求解方法是在离线状态下进行的，因此使用已知禁飞区具体信息求解，同时侧滑角在最小值与最大值构成的区间内自由取值。另一方面采取本文的机动制导策略进行在线求解，求解过程只能通过触角的探测，反馈给制导策略禁飞区部分轮廓信息，侧滑角的取值只有最大值、最小值和零。如图10和图11所示，基于“触角”的方法与基于GPOPS II的方法都可以使两种飞行器顺利规避禁飞区到达终点，同时，如图12和图13所示，由于基于“触角”的方法的侧滑角的选取区间相对简单，路径图中的轨迹相对平滑。

图9 仿真算法流程框图

图10 CAV-H飞行器规避禁飞区路径

图11 Shuttle飞行器规避禁飞区路径

图12 CAV-H控制变量变化

图14和图15显示，基于“触角”的机动制导策略在规避禁飞区的同时，由于速度和高度严格跟踪标准轨迹，虽然在算法中触角终止条件没有判别是否符合3个经典路径约束条件，机动制导的结果依旧满足约束条件。

图14 CAV-H热流率、动压、过载变化

图15 Shuttle热流率、动压、过载变化

由图16中可以看出：飞行过程中的单次预测周期计算时间是呈下降趋势的，其中的大小变化主要是由预测时飞行器位置决定的，飞行器距离禁飞区、航向角约束或终点越近，计算的时间就越短。

图16 飞行过程中单次预测周期计算时间

3.2 模型的鲁棒性与计算时间

为了验证本文基于三“触角”的侧滑角瞬变机动制导模型，对两个飞行器进行鲁棒性仿真实验。采用蒙特卡洛方法，对初始状态中的纬度、航迹角、航向角以及飞行器的质量、面积实施3°、3°、10°和5%、5%的离散，经仿真获得了每个飞行器的1 000次结果，从图17和图18中可以看出虽然有一些仿真路径的机动幅度较大，但是所有仿真轨迹均能顺利规避本文所设计的复杂禁飞区到达各自的标准轨迹终点。

图17 CAV-H蒙特卡洛仿真结果

图18 Shuttle蒙特卡洛仿真结果

本文进行实验时使用的计算机CPU主频为2.6 GHz，软件为Matlab。在仿真实验中，统计每一个仿真实验的单次预测周期计算时间的最大值、最小值、平均值，得到图19和图20，可以看出，本文提出的机动制导策略在预测时最大预测时间不超过0.25 s，而且平均预测时间集中在0.01 s，完全符合机动制导的时间要求。同时，图21与图22展示了仿真中计算得到的热流率、动压和过载最大值均满足表2的约束条件。

图20 Shuttle单预测周期计算时间统计图

图21 CAV-H热流率、动压、过载最大值点集

图22 Shuttle热流率、动压、过载最大值点集

4 结论

本文通过对基于“触角”的预测校正制导方法的改进，综合航向角约束和侧滑角延时计数开关设计，提出了基于三“触角”的侧滑角瞬变机动制导策略，解决了高超声速滑翔飞行器在多禁飞区信息未知情况下再入机动制导的问题。通过仿真验证了该方法能够有效解算出多禁飞区情况下高超声速滑翔飞行器机动制导的在线飞行方案，通过蒙特卡洛仿真实验验证，突出体现了该方法计算时间上能够满足机动制导的要求。