基于深度强化学习与高度速率反馈的再入制导方法

2022-09-27武天才王宏伦刘一恒

无人系统技术 2022年4期

武天才，王宏伦，刘一恒，任斌，余跃

（1.北京航空航天大学自动化科学与电气工程学院，北京 100191；2.北京航空航天大学高等理工学院，北京 100191；3.北京航空航天大学飞行器控制一体化技术重点实验室，北京 100191；4.北京航天自动控制研究所，北京 100854）

1 引言

高超声速飞行器由于其飞行速度快、飞行空域广、机动性能好、远程打击精确等优势，在军事和民用领域备受关注[1-2]。然而，在再入飞行过程中，由于受到复杂环境、模型不确定性和严苛路径约束等因素的制约，再入制导技术成为高超声速飞行技术领域的核心和关键技术之一，也是当前各国空天领域研究的重点与难点[3]。

高超声速飞行器的再入制导算法一般分为标称轨迹制导与预测校正制导两种[4]。虽然标称轨迹制导方法具有计算量小、实时性高和对机载计算机性能要求较低等优点，但是其鲁棒性较差，初始散布偏差和环境不确定性会对其制导性能产生严重影响。近年来，随着机载计算机性能的不断提高，具有在线校正制导指令能力的预测校正制导方法逐渐成为主流[5]。文献[6]基于基础的预测校正制导方法，给出了适用于不同种类高超声速飞行器再入制导方案的统一设计框架。为解决热流密度、动压和过载等路径约束对飞行器的限制，文献[7]提出了一种受限预测校正制导方法，该方法借助准平衡滑翔条件（Quasi-Equilibrium Glide Condition,QEGC），巧妙地将再入路径约束转换为对倾侧角幅值的约束。

对于具有高升阻比的高超声速飞行器而言，另一个需要关注的现象是：在再入滑翔过程中，飞行器高度会呈现一种周期性振荡的现象。这种振荡会使飞行器逾越路径约束，甚至使准平衡滑翔条件失效，严重影响系统的稳定性。文献[8]在预测校正制导方法的基础上，设计了一种参数化的反馈控制律，可以根据飞行器高度变化速率进行倾侧角幅值修正来抑制振荡现象。文献[9]则是以高度变化率及空速为输入，设计模糊控制器对倾侧角幅值进行调节来抑制振荡。但是这两种方法一定程度上依赖于参数或模糊规则设计，需要人工不断迭代校正才能提升效果。

近年来，以深度学习、深度强化学习为代表的新一代人工智能方法广泛应用于各类复杂系统以提高传统方法的性能[10-11]。为解决无人机在复杂环境下的航路规划问题，文献[12]将深度神经网络与传统流体扰动方法结合，提升了算法的实时性和自适应能力。文献[13]将导弹目标拦截中的制导控制一体化问题转化为强化学习优化问题，通过设计启发式的奖励函数使导弹可以低能耗且快速稳定地拦截目标。文献[14]利用深度强化学习Q-learning 算法在线决策飞行器倾侧角的符号指令，消除了传统方法在侧向制导逻辑设计时增加的一些不必要约束，并有望应用于规避多禁飞区的轨迹规划研究。文献[15]利用深度强化学习算法在线决策飞行器倾侧角指令。与传统方法相比，该方法充分发挥飞行器的宽域飞行优势，进一步拓展了飞行剖面。

基于上述分析，为解决具有高升阻比的高超声速飞行器在再入制导中高度存在的周期性振荡现象，本文在传统预测校正制导方法的基础上引入高度速率反馈，并以此为基础构建深度强化学习训练环境，采用深度确定性策略梯度算法（Deep Deterministic Policy Gradient,DDPG）训练动作网络，以根据飞行器当前状态输出合适的反馈增益，然后，将训练好的动作网络提取至上述传统制导方法框架中，实现反馈增益的在线实时调整，接着，进一步结合飞行器高度速率反馈进行倾侧角指令补偿量的在线实时求取，来抑制高度的周期性振荡，有效避免了传统平稳滑翔制导算法设计中需要依据人工经验对算法参数进行不断迭代校正的这一冗余过程。

2 再入制导问题

2.1 三自由度运动学模型

为方便高超声速飞行器再入制导问题的研究，在考虑地球自转的基础上，建立无量纲化的高超声速飞行器三自由度质点模型[16]：

式中，r表示无量纲化的飞行器地心距，V表示无量纲化的飞行器速度，λ和φ分别表示飞行器所处位置的经度和纬度，θ和ψs分别表示飞行器的航迹倾角和航迹偏角，γs表示飞行器的倾侧角，L=QSC L/mg0和D=QSC D/mg0分别表示飞行器的升力加速度和阻力加速度，Q表示飞行器动压，CL和CD分别表示升力和阻力系数，一般为与飞行器攻角α相关的非线性函数，m表示飞行器质量，g0表示海平面重力加速度，CV，Cθ和Cψs分别表示由地球自转造成的作用项，定义如下：

高超声速飞行器时间、地心距、速度和地球自转角速度对应的无量纲化过程定义如下：

式中，R0表示地球半径，分别表示真实的时间、地心距、飞行速度和地球自转角速度。

2.2 再入过程约束

高超声速飞行器再入飞行过程考虑的约束包含：包括热流密度约束Q、动压约束q和过载约束n在内的“硬”约束条件和包括准平衡滑翔条件在内的“软”约束条件，具体计算公式如下[8]：

式中，KQ表示与飞行器结构和材料相关的常值参数；ρ表示大气密度，本项目中采用指数大气模型ρ=ρ0exp(-βh)；β=1.3785e-4,qmax,nmax分别表示热流密度约束、动压约束和过载约束的最大幅值。

2.3 再入终端约束

由于高超声速飞行器再入终端时间不固定，引入如下能量形式的e作为自变量[5]：

式中，能量e为单调递增函数，可作为飞行器是否抵达再入终端的判断条件。

再入终端约束主要包括高度约束、速度约束和位置（经纬度）约束，用公式表示如下：

式中，ef=1/r f-Vf2/2表示终端能量，不失一般性，高度约束在本文中通过终端地心距r(ef)约束来表示，rf,V f,λ f,φf分别表示期望的终端地心距、速度和经纬度。

3 再入制导算法设计

本文所提的基于深度强化学习和高度速率反馈的再入制导方法以传统预测校正制导方法为基础，进行初始纵向和侧向制导指令的求取；接着，引入高度速率反馈进行倾侧角幅值补偿量的求取，来抑制高度存在的周期性振荡现象。针对反馈增益的设计问题，本文采用深度强化学习方法进行解决：首先，构建以传统预测校正制导和高度速率反馈为基础的深度强化学习训练环境；然后，采用深度确定性策略梯度算法进行离线训练；最后，提取出训练好的动作网络，与上述基础制导算法结合进行在线使用，实现反馈增益的实时调整与倾侧角幅值补偿量的实时求取，所提制导算法结构示意图如图1所示。

图1 再入制导算法结构图Fig.1 Schematic diagram of the reentry guidance algorithm

3.1 纵向制导算法设计

纵向制导算法设计包括攻角剖面和倾侧角幅值剖面的设计。

首先给出攻角剖面设计方法，考虑到飞行器初始下降段的热防护要求，本文中攻角剖面采用升力式飞行器常用的分段标称攻角剖面，计算公式如下[9]：

式中，αmax,αmin为再入攻角边界值；V1,V2为飞行器攻角剖面临界速度。因此，倾侧角γs成为改变飞行器再入轨迹的唯一控制变量。

根据设计好的攻角剖面，再入过程约束式(11)可以转换为高度-速度（h-）剖面飞行再入走廊的边界（如图2所示）[15]：

图2 h -平面内再入走廊示意图Fig.2 Schematic diagram of reentry corridor inh - plane

为进一步方便在制导算法设计中考虑上述约束的限制，借助准平衡滑翔条件，再入走廊的下边界约束可进一步转换为倾侧角幅值的约束：

纵向制导算法设计的另一个目标是进行倾侧角幅值剖面的设计，为了提升制导性能，倾侧角幅值剖面通常设计为如下与能量相关的线性分段函数[6]：

式中，能量e定义与式(12)相同，γs,0为每个制导周期求得的倾侧角幅值，γs,f是一个预先设计的常数，例如，可令γs,f=γs,0。因此，只要确定了γs,0的取值，就可以得到整个倾侧角幅值剖面。

式中，λpt和φpt为预测的终端落点的经纬度。

当前状态到期望终端状态的待飞航程为：

待飞航程偏差可定义为：f(γs,0)=Spt-Sf。然后，采用割线法进行迭代求取使待飞航程偏差f(γs,0)=0的解[5]：

3.2 侧向制导算法设计

侧向制导算法的目标是确定倾侧角指令的符号，本文通过航向角误差走廊来进行倾侧角翻转逻辑的设计。

首先，给出飞行器当前位置到终端目标点的视线方位角Φ的定义[5]：

可以得到航向角误差：Δψs=ψs-Φ。为满足终端精度要求，采用如图3所示的航向角误差走廊。

图3 航向角误差走廊示意图Fig.3 Schematic diagram of heading angle error corridor

倾侧角翻转逻辑设计为：当航向角误差 Δψs超出误差走廊的上边界 Δψs,max(e)时，倾侧角指令符号为负；当航向角误差 Δψs超出走廊的下边界 Δψs,min(e)时，倾侧角指令为正；当航向角误差Δψs位于误差走廊内时，倾侧角符号保持不变。倾侧角翻转逻辑可用公式表示为：

3.3 基于深度强化学习与高度速率反馈的平稳滑翔制导律设计

本文平稳滑翔制导律的设计思想是：在传统预测校正制导算法基础上引入高度速率反馈，并以此为基础构建深度强化学习训练环境，然后通过深度确定性策略梯度算法训练动作网络，使该网络可以根据飞行器当前状态输出反馈增益。然后，将训练好的动作网络提取至上述制导算法框架中，实现反馈增益的在线实时调整，进一步结合高度速率反馈计算得出倾侧角指令补偿量，来抑制飞行器高度振荡。

深度强化学习解决的是马尔可夫决策问题（Markov Decision Process,MDP）。MDP 包括状态、动作、策略和奖励等要素，其中策略一般由深度神经网络建模，由强化学习算法训练，所采用的状态输入、动作输出以及奖励函数需要进行设计。

3.3.1 深度强化学习状态建模

所采用的状态需要能完整反映高超当前的状态，本项目采用的深度强化学习状态设计为

式中，r,V,θ,ψs表示飞行器当前状态，表示飞行器当前攻角指令、倾侧角指令幅值和倾侧角指令符号，Δ,λ,Δφ表示飞行器当前位置与期望目标位置之差，表示飞行器当前高度速率与参考高度速率之差。

对于具有高升阻比的高超声速飞行器而言，在再入制导中高度存在周期性振荡的现象，如图4所示。通常，高超声速飞行器再入大气层后会在首次拉起过程中产生最明显的轨迹振荡，之后，振荡幅值逐渐衰弱。参考高度速率求取采用文献[8]提出的方法，在轨迹振荡的多个波峰和波谷分别进行插值计算得到高度速率的上边界和下边界，然后计算算数平均值即可得到。

图4 参考高度变化率示意图Fig.4 Schematic diagram of reference value of height gradient

3.3.2 深度强化学习动作建模

本文设计的深度强化学习动作为

式中，kDRL表示动作网络给出的反馈增益。

3.3.3 深度强化学习奖励函数设计

本文设计的深度强化学习奖励函数为

式中，rdis表示与飞行器和期望目标位置距离相关的奖励函数，表示与飞行器高度速率和参考高度速率偏差相关的奖励函数，rV和rh表示飞行器终端高度、速度和期望终端高度、速度相关的奖励函数，rdis，，rV和rh分别设计如下

式中，dis=arccos(sinφsin(φf)+cosφcos(φf)cos(λf-λ))表示飞行器当前位置与期望目标位置的距离，Rf表示期望目标区域半径，“500”为强奖励项，表示飞行器达到目标区域后给一个极大的正奖励；“-dis”为弱激励项，只有飞行器离终点的距离越来越近，负奖励才会越来越小，从而引导飞行器向终点飞行。

式中，为飞行器高度速率与参考高度速率之差的奖励增益，本文取=200，表示当飞行器高度速率接近参考高度速率时，所给负奖励较小；当飞行器高度速率远离参考高度速率时，给较大的负奖励，来驱使飞行器尽可能按照参考高度速率飞行，抑制飞行器高度振荡现象。

式中，kh为飞行器终端高度与期望终端高度之差的奖励增益，本文取kh=1，IsDone表示当前训练回合的终止条件，由于倾侧角补偿量的引入，会影响飞行器终端高度的制导精度。因此，引入飞行器终端高度与期望终端高度之差的奖励项，来增强飞行器终端高度的制导精度。

式中，kV为飞行器终端速度与期望终端速度之差的奖励增益，本文取kV=1。同样，由于倾侧角补偿量的引入，会影响飞行器终端速度的制导精度。因此，引入飞行器终端速度与期望终端速度之差的奖励项，来增强飞行器终端速度的制导精度。

当前训练回合的终止条件IsDone设置为：

式中，e＞ef表示飞行器当前能量大于期望终端能量，达到了终止条件；(h＞hmax)∪ (h＜hmin)和(V＞Vmax)∪ (V＜Vmin)表示飞行器当前状态超过合理范围。

3.3.4 基于DDPG 的深度强化学习算法

本文DDPG[17]算法进行深度强化学习动作网络训练，深度强化学习DDPG 算法结构示意图如图5所示。深度强化学习DDPG 算法流程如下述伪代码所示。

图5 深度强化学习DDPG 算法结构示意图Fig.5 Structure diagram of DDPG algorithm

对上述过程进行充分的迭代，可以得到训练好的动作现实网络。

将动作现实网络输出的反馈增益kDRL与飞行器当前高度速率与参考高度速率之差相乘得到倾侧角指令补偿量γs,DRL，再将γs,DRL线性叠加到预测校正制导律的输出上：

4 仿真结果及分析

本文以高超声速飞行器CAV-H 为仿真对象来验证所提算法的有效性。飞行器初始状态为：高度h0=80 km，速度=7100 m/s，经度λ0=10°，纬度φ0=-20°，航迹倾角θ0=-1°，航迹偏角ψs,0=45°。飞行器终端约束为：高度hf=20 km，速度=1800 m/s，经度λf=90°，纬度φf=30°，终端位置误差Rf≤10 km。飞行器过程约束为：热流密度约束=1.5 MW/m2，动压约束qmax=200 kPa，过载约束nmax=4.5。

4.1 深度强化学习训练结果

根据第3 节设计的深度强化学习训练环境开展所提制导方法的训练，训练选取的飞行器初始状态散布偏差和参数摄动情况如表1所示。

表1 初始状态散布偏差和参数摄动Table 1 Initial state dispersion and parameter perturbation

图6给出了深度强化学习训练过程的奖励函数平均值曲线，可以发现DDPG 算法可在训练过程的第100 个回合左右使奖励函数进入收敛状态，验证了所提制导方法通过深度强化学习方法训练动作网络以实现反馈增益在线自适应调整的可行性。

图6 DDPG 训练过程中的奖励函数曲线Fig.6 Reward function curve in DDPG training process

4.2 对比仿真

在标称条件下进行对比仿真，来验证所提制导算法的有效性。采用传统预测校正制导方法（为表述简洁，以下称为“传统方法”）和文献[8]所提的基于倾侧角反馈控制的预测校正再入制导方法（以下称为“传统反馈方法”）作为对比项进行仿真对比，仿真结果如图7～15 所示。

由图7给出的飞行器高度曲线可得，在传统方法作用下，飞行器高度存在明显的振荡现象；在传统反馈方法作用下，飞行器高度振荡现象虽有些许改善，但是首次拉起时飞行器弹起高度较高，极有可能在再次下落时导致热流密度超越限幅，影响飞行器的热防护安全；相比之下，在所提方法的作用下，动作网络可以根据飞行器当前状态给出一个更优的反馈增益系数，进而计算得到一个更加合适的倾侧角补偿量，使得飞行器在再入过程中可以平稳飞行。

图7 对比仿真中的高度曲线Fig.7 Altitude curve in comparative simulation

由图8给出的飞行器经纬度曲线可以发现，在传统方法、传统反馈方法和所提方法作用下，飞行器终端位置误差分别为：5.25 km，4.93 km，4.25 km，均符合制导精度要求，所提方法增加了对高度振荡现象的抑制，但并未影响最终的制导精度，验证了所提方法的有效性。对图7和图9的终端误差进行统计，三种方法的终端高度误差的绝对值分别为：0.85 km，0.96 km，0.06 km、三种方法的终端速度误差的绝对值分别为：3.95 m/s，2.18 m/s，0.15 m/s，可以发现，所提方法在终端高度和终端速度方面均取得了最优的制导精度，这是由于在深度强化学习奖励函数设置中加入了与终端高度和终端速度的奖励项，动作网络在输出反馈增益时，选择了使奖励更高的输出。

图8 对比仿真中的经纬度曲线Fig.8 Latitude and longitude curve in comparative simulation

图9 对比仿真中的速度曲线Fig.9 Velocity curve in comparative simulation

图10给出了三种方法的倾侧角曲线，可以发现所提方法仅在倾侧角的幅值上做出了修正，没有增加制导指令反转的次数，体现了所提方法的可靠性。图11给出了所提方法的倾侧角补偿量曲线，由图可得，所提方法从飞行器首次拉起时开始进行倾侧角补偿，首先补偿一个较大倾侧角，分散飞行器升力作用，来避免飞行器弹起更高的高度；然后，根据飞行器当前高度与期望高度速率关系进行倾侧角指令微调，实现飞行器平稳再入滑翔。

图10 对比仿真中的倾侧角曲线Fig.10 Bank angle curve in comparative simulation

图11 对比仿真中的倾侧角补偿量曲线Fig.11 Compensated bank angle curve in comparative simulation

图12～14 给出了在三种方法作用下飞行器热流密度、动压和过载曲线，由热流密度曲线可以明显发现在传统方法和传统反馈方法的作用下，由于未对飞行器高度振荡现象进行抑制或进行有效抑制，在飞行器高度振荡过程中，极有可能逾越热流密度的约束，给飞行器结构带来损害；而在所提方法的作用下，高度振荡现象被有效补偿，热流密度远小于限幅值，取得了更好的热防护效果。由动压曲线可以明显发现，虽然在三种方法作用下，飞行器均满足动压约束限幅，但所提方法的动压更加平稳，保证了飞行的稳定性。由过载曲线可以明显发现，由于传统方法和传统反馈方法未能对飞行器高度振荡进行有效抑制，在飞行末段存在过载超过限幅的情况，体现出飞行器再入过程中需要对高度振荡现象抑制的必要性，也再次验证所提方法的有效性。

图12 对比仿真中的热流密度曲线Fig.12 Heating rate curve in comparative simulation

图13 对比仿真中的动压曲线Fig.13 Dynamic pressure curve in comparative simulation

图14 对比仿真中的过载曲线Fig.14 Overload curve in comparative simulation

图15给出了在每个制导周期所提方法的制导指令解算时间统计图，可以发现制导指令的解算时间不超过0.25 s，制导指令解算时间平均值小于0.1s，满足制导周期（T=1s）的实时性要求。可以发现，将训练好的动作网络引入制导算法框架中并不会严重影响制导指令的解算时间。这是因为在在线使用过程中，仅需要根据当前状态进行动作网络前项传播过程的计算，计算较为简单，不会占用大量计算资源和消耗大量的计算时间。

4.3 蒙特卡洛仿真

为进一步验证所提算法的自适应能力和鲁棒性，在表1所述的初始状态散布偏差和参数摄动范围内进行50 次蒙特卡洛仿真，仿真结果如图16～20 所示。

图16～17 给出了蒙特卡洛仿真中的经纬度曲线及终端精度统计图。可以发现，在飞行器初始状态存在散布偏差和气动参数存在摄动时，在所提制导算法作用下，飞行器终端位置误差均保持在10 km 范围内，验证了所提算法的有效性和鲁棒性。

图16 蒙特卡洛仿真中的经纬度曲线Fig.16 Latitude and longitude curve in Monte Carlo simulation

图17 蒙特卡洛仿真中的终端精度统计图Fig.17 Terminal accuracy statistics in Monte Carlo simulation

图18～19 给出了蒙特卡洛仿真中的高度和速度曲线，终端高度和速度误差的绝对值不超过2 km 和100 m/s。图20～21 给出了蒙特卡洛仿真中的倾侧角和倾侧角补偿量曲线，在所提制导方法作用下，动作网络可根据飞行器当前状态实时选取合适的反馈增益，计算倾侧角反馈补偿量，有效抑制高度存在的周期性振荡现象，使整个飞行过程趋于平稳，验证了本文所提方法的自适应能力。