空天飞行器返回滑翔段在线制导方法

2021-03-13周宏宇王小刚赵亚丽崔乃刚

宇航学报 2021年2期

周宏宇，王小刚，赵亚丽，崔乃刚

(1.哈尔滨工业大学航天学院，哈尔滨 150001；2.北京航天晨信科技有限公司，北京 102308)

0 引言

随着航天技术的快速发展和航天活动的多元化与频繁化，航天发射的经济性和灵活性愈发重要；在此背景下，可完全重复使用的空天往返飞行器应运而生。不同于运载火箭和航天飞机，空天飞行器采用多种动力模式入轨，上下两级均有自主返回能力，在灵活性和经济性上更具优势，已成为各航天强国的研究热点[1-4]。美国近年来完成了多次空天飞行器相关试验[5-6]，印度于2016年完成了可重复使用验证机首飞[7]，日本也于2019年制定了新一代空天飞行器研制计划。

返回滑翔段是空天飞行器实现可重复使用的关键阶段。由于再入速度高、滑翔距离远，为实现安全返回，必须满足热环境、过载和拟平衡滑翔以及攻角和倾侧角等控制约束[8]；同时还要满足高度、速度、倾角、位置和航向等终端约束，从而顺利过渡至末端能量管理段。此外，考虑初始返回状态的不确定性和诸多干扰因素的影响，必须对滑翔段轨迹进行在线制导修正。

目前可用于空天飞行器返回滑翔段的制导方法包括标称轨迹法[9-11]、预测校正法[12-14]以及在线优化法[15-16]等。其中文献[9]给出了一种较为典型的标称轨迹法，该方法事先设定攻角剖面并在飞行走廊中迭代倾侧角，然后使用PID控制器跟踪标称攻角及倾侧角指令。预测校正法通过在线求解终端状态修正飞行指令，无需轨迹跟踪器，抗干扰性更强。在线优化法通过反复计算最优轨迹修正飞行指令，能够满足多种飞行约束并保证飞行指标的最优性，具有很好的应用前景。

上述方法技术成熟、应用广泛，但在空天飞行器应用中还需进一步研究与改进。首先，标称轨迹法一般需要满足小扰动假设并预设攻角形式，制约了任务的灵活性。预测校正法需在线计算终端偏差，但常见的基于数值解或解析解的预测方法难以兼顾计算量和计算精度。由于滑翔段航程远、飞行环境复杂、非线性耦合强[17]，快速求解最优轨迹并非易事，故在线优化法需同时考虑最优性和收敛性的问题。为解决上述问题，诸多学者开始探索新的理论方法和技术途径。基于滑模控制理论，文献[18]设计了一种能够克服大偏差的自适应轨迹跟踪器。文献[19]设计了一种改进型扰动观测器，提高了标称高度和速度剖面的跟踪精度。文献[20]采用反馈线性化理论，基于高低角和方位角信息实现了终端速度的快速预测。文献[21]提出了一种快速优化算法，保证了最优轨迹的在线求解效率。同时一些混合方法正不断涌现，如将标称轨迹法和在线优化法结合，从而实现制导精度和计算量的最佳平衡。

针对滑翔轨迹在线重构需求，考虑气动、约束、轨迹和指标间的复杂耦合关系[22]，从提高在线求解效率出发，设计了一种新的滑翔段飞行剖面。推导了滑翔段运动状态、过程约束和性能指标的解析表达式，降低了优化问题的求解难度。在此基础上提出了一种双层在线制导方法：内层实时重构飞行剖面并借助解析解控制终端速度；外层使用混合优化算法优化飞行剖面，满足过程约束的同时降低过载需求。该方法无需事先计算攻角与倾侧角，无需设计轨迹跟踪器，无需在线积分预测终端状态，能够自动满足多种约束条件并保证轨迹的最优性。数学仿真验证了方法的正确性和有效性；在非标称状态下，终端状态相对偏差小于3%。

1 返回滑翔段运动数学模型

考虑地球自转，空天飞行器在返回滑翔段中的运动方程如下：

(1)

式中：V为速度，γ为飞行路径角，r为质心到地心的距离，ψ为航向角，θ为经度，φ为地心纬度；m为质量，D为阻力，L为升力，σ为倾侧角，g为地球引力加速度，ωe为地球自转角速度。

阻力和升力的计算方式为：

(2)

式中：cD和cL分别为阻力系数和升力系数，Sref为飞行器特征面积，q=ρV2/2为飞行动压，其中ρ为大气密度，可采用指数函数计算：

ρ=ρ0exp(-βh)

(3)

式中：β=1.41×10-4m-1与ρ0=1.225 kg/m3为常系数；h=r-Re为飞行高度，其中Re为地球平均半径。

设气动系数按如下规律变化：

(4)

式中：cD0,cL0,β1,β2,β3,β4为给定参数，α为攻角。

定义当前位置与终端位置间的航程为剩余航程，记为RL，则有：

(5)

式中：μL为剩余航程对应的地心航程角，θT和φT表示终端位置；ζ=ψw-ψ为航向偏差，其中ψw为期望航向角，计算方式如下：

(6)

空天飞行器在返回滑翔段中必须满足诸多约束条件。为保证顺利交班，设终端约束如下：

Φf=[hf,γf,Vf,RLf,ζf]T-[hc,γc,Vc,0,0]T=0

(7)

式中：下标“f”表示实际终端状态，下标“c”表示期望终端状态。

同时，设置过程约束如下：

(8)

式中：下标“max”表示该过程变量允许的最大值，NL=L/(mg)为速度坐标系下的法向过载，qr为驻点热流。对于空天飞行器，热流的计算方式如下：

(9)

式中：RN为驻点曲率半径，取为0.08 m。

此外，滑翔段还应满足拟平衡滑翔条件，具体表示为如下形式：

(10)

2 返回滑翔段飞行剖面设计

2.1 纵向剖面设计

设计飞行高度为剩余航程的函数：

(11)

式中：ai(i=0,1,…,n)为未知系数。取n=5。

1)攻角确定

(12)

(13)

其中

(14)

根据式(1)和(2)，式(12)和(13)可表示为：

(15)

(16)

(17)

定义：

(18)

(19)

(20)

因此飞行剖面上各处的升力为：

(21)

注1.式中令ψw=arccos(Kw)，结合式可得：

(22)

2)剖面确定

由于n=5，需建立六个方程来求解F(RL)中的系数a0～a5。由于初始和终端高度已知，由式(11)可得：

(23)

式中：下标“0”表示滑翔段初始运动状态。同时根据初始和终端飞行路径角，由式(12)和(15)可得：

(24)

(25)

假设cosζ≈1，式(24)和(25)可改写为：

(26)

(27)

定义：

(28)

类似于式(23)，可以得到另外两个方程：

(29)

式中：h1和h2分别为剩余航程2/3和1/3处的高度值。定义剖面设计变量为uh=[h1,h2]T，因此只需两个参数便能够由式(23)～(29)求出a0～a5。

最后，给出剖面设计的性能指标函数。由式(21)可知，升力由f″(RL)决定；因此在满足约束条件的前提下，若能使|f″(RL)|尽量小，则需用过载和攻角会更小。因此设计性能指标如下：

minJglide=max(|f″(RL)|)

(30)

2.2 侧向剖面设计

设计倾侧角为航向角偏差的函数：

(31)

式中：为ζm为给定正数，kσ为待定正数；设t0为倾侧角翻转时刻，则Δt=t-t0且σ0=σ(t0)。当t=t0时，σ=σ0；当t足够大时，σ由航向偏差ζ决定。

(32)

定义Kσ=ζσmax/ζm-σ0，考虑到e-kσΔt≤1：

(33)

2.3 滑翔段解析解推导

2.3.1运动状态的解析解

由于高度已定义为剩余航程的函数，以剩余航程为自变量，滑翔段高度的解析解即式(11)。

设cosζ≈1，cosγ≈1和sinγ≈γ；根据式(1)、(5)和(12)可得飞行路径角的解析解为：

(34)

为求得速度的解析解，由动力学方程可得：

(35)

滑翔段中阻力通常远大于引力沿速度方向的分量，因此式(35)可改写为：

(36)

定义τD=cD0Srefρ0/(2m)，则有：

(37)

由式(11)、(34)和(37)可得：

(38)

对上式两端积分得：

(39)

式中：CV为积分常数。定义：

(40)

因此滑翔段速度的解析解为：

V=[τDy(RL)+CV]-1/β1

(41)

由于RL0对应V0，可得CV=V0-β1-τDy(RL0)。同时由于RLf=0，故y(RLf)=0。因此终端速度为：

(42)

2.3.2过程约束和性能指标的解析解

根据式(41)，动压的解析表达式为：

(43)

同时，驻点热流的解析解为：

(44)

式中：ks=9.98×10-7，此时qr的单位为W/m2。

设cosζ=cosγ=1，则升力的解析解为：

(45)

因此法向过载的解析表达式为：

(46)

欲求性能指标即|f″(RL)|的极值，令f‴(RL)=0；由式(13)可得：

(47)

其中

(48)

求解方程F2(RL)G1-2G2=0，便可得到f(RL)二阶导数的极值，进而求出性能指标函数Jglide。

3 返回滑翔段在线制导算法

3.1 在线制导方法

通过设计纵向剖面和侧向过载，终端高度、飞行路径角、剩余航程和滑翔偏差得以自动满足；而式表明终端速度是剩余航程的函数，因此可通过控制横向机动改变剩余航程的变化率，使剩余航程匹配期望的终端速度。

通过设计路径点来改变剩余航程。如图1所示；其中p1为当前位置，p2为路径点，p3为目标点。R12为p1p2间的航程，R23为p2p3间的航程，R13为p1p3间的航程。为满足终端速度约束Vfc，由式(40)和(42)可得剩余航程应满足：

图1 剩余航程控制方法

(49)

记方程(49)的解为RLw，则路径点p2的选择必须满足R12+R23=RLw。参考式可得：

(50)

式中：θ2和φ2为路径点p2的经纬度，ψ12为指向p2的期望航向角。

令R12=R23=RLw/2，由球面三角形定理可得：

(51)

式中：ψ13表示经过p2后指向目标点的期望航向角。因此，路径点p2的位置为：

(52)

设从起点到当前位置飞过的航程为Rcov，制导周期为ΔRcov；在线制导方法具体如下所示：

1)路径点更新

每五个周期(5ΔRcov)更新一次期望剩余航程RLw和路径点p2；

2)剖面优化计算

每四个周期(4ΔRcov)计算一次uh=[h1,h2]T的最优值，使指标Jguide达到最小；否则uh的值不变；

3)剖面在线重构

每一个周期(ΔRcov)根据uh和当前运动状态，由式(23)～(29)更新多项式F(RL)的系数a0～a5；

4)分别由式(21)和(31)计算攻角和倾侧角；

5)更新参数τD；

6)积分动力学方程更新运动状态；

7)返回步骤1)。

步骤1)是为了修正终端速度。步骤2)是为了更新剖面重构基准uh=[h1,h2]T，修正速度解析值与真实值间的偏差，满足过程约束并保证轨迹的最优性。步骤3)是为了重构飞行剖面，从而满足终端高度、位置和飞行路径角约束。

3.2 混合优化方法

上述在线滑翔剖面优化问题实际上是一个含有五个过程(不等式)约束的双参数寻优问题。只含不等式约束的优化问题可通过罚函数法转化为无约束优化问题：

(53)

式中：rΦ为罚因子，NΦ为不等式约束的数量。上述无约束优化问题可采用CGM求解。基本CGM如下：

xl+1=xl+αldl

(54)

(55)

式中：αl,dl和βl分别为搜索步长、搜索方向和搜索因子；x=Uh=[h1,h2]T为优化变量，l为迭代次数；Gl为函数Jnew关于xl的梯度，可采用差分法求得。

标准CGM中的搜索步长为单调递减，在迭代后期计算效率较低；因此按如下方式选取步长αl：

(56)

式中：δ1∈(0,1/2]和δ2∈(δ1,1)为预设常数。

为避免计算变量βl，将式(55)第二项改写为：

dl=-QlGl,l>0

(57)

式中：Ql为二维方阵。定义sl=dl-dl-1,zl=Gl-Gl-1以及二维单位方阵I2，则Ql可按如下方式计算：

(58)

改进后的CGM如下：

1)令迭代次数l=0，初始化x0；

2)由式(55)第一项和式(57)确定搜索方向dl；

3)由式(56)确定搜索步长αl；

(1)给定αl的值，求出xl+1以及相应的性能指标梯度G(xl+1)；

(2)根据式(56)调整αl值，直至式(56)成立；

4)由式(54)计算xl+1；

5)根据xl+1判断算法是否收敛：

(59)

式中：εj和εx为收敛精度。若满足式(59)则终止迭代并输出xl+1，否则l加1并返回步骤2)。

CGM在初值适合时收敛很快，而文献[23]指出在惯性权重的值满足一定条件时，PSO算法将快速收敛于近优解附近；因此，PSO[24]可作为CGM的初值求解器，由于初步计算h1和h2。

4 仿真分析

为验证方法的正确性，建立表1所示仿真场景，其中初始航向角由式(6)求出，故初始航向偏差为零。设飞行器质量分别为6000 kg,5500 kg和5000 kg，特征面积分别为4.5 m2,4.0 m2和3.5 m2。

表1 仿真场景设定

首先不考虑干扰和偏差，得到仿真结果如图2～图7所示。图2表明高度变化平缓，图4表明飞行路径角的绝对值始终小于2°(大多数时间小于1°)，满足拟平衡滑翔条件和解析解中的前提假设。

图2 滑翔段飞行剖面

图3 飞行速度随时间变化情况

图4 飞行路径角随时间变化情况

图5 过程约束随时间变化情况

图6 飞行指令随时间变化情况

图7 三维飞行轨迹

由于在线搜索剖面设计变量时使用了大量解析表达式，即使需要计算约20×15=300次性能指标，PSO也能很快得到初值。以场景1中首次(零时刻)优化计算为例，在东芝L001笔记本电脑(双核4 GB内存、英特尔酷睿i3 M330/2.13 GHz处理器)上基于MATLAB 2016a执行混合优化算法，具体计算过程如图8所示。

图8 混合优化算法计算过程

由于上述滑翔段剖面设计方法能够自动满足终端高度、飞行路径角和位置约束，在此仅关注终端速度偏差；如图3所示，三种场景下，终端速度偏差分别为23.25 m/s、20.98 m/s和14.31 m/s；相对误差分别为1.94%、1.95%和1.19%。

图5表明各项过程约束均得到满足。图6表明攻角指令十分平滑。由于攻角源于函数F(RL)的二阶导，因此攻角变化率是连续的，同时倾侧角指令也是连续光滑，这对控制系统是很有利的。在路径点处，期望航向的变化导致了侧向过载方向的改变，因此倾侧角会出现突变，同时攻角也会出现突变；但在倾侧角计算过程中已经考虑了精度变化率约束，因此路径点处的飞行指令变化幅度满足控制系统的能力约束。同时，在滑翔段的绝大多数时间里，攻角和倾侧角指令均是平滑连续的。

图8中，PSO用时0.2 s(优化结果为h1=36.5 km，h2=35.3 km)，CGM用时0.06 s(迭代4次后收敛，结果为h1=38.4 km，h2=35.6 km)，说明基于解析解的在线优化算法满足计算效率要求。另外，图8表明PSO算法在进化第11代时已经收敛，说明可以用更少的进化代数或粒子数完成初值计算。但PSO具有随机性，故不宜采用过少的进化代数或粒子数。

其次考虑干扰因素，设置服从均匀分布的滑翔段干扰因素如表2所示。以表1中的场景1为例进行500次蒙特卡洛打靶仿真，结果见图9～图12。

图12 过程约束满足情况

表2 滑翔段干扰因素

如图9所示，终端速度偏差的绝对值在500次仿真中的最大值为34.98 m/s(相对误差2.91%)，平均值为20.34 m/s(相对误差1.69%)。

图9 终端速度偏差直方图

图10和图11表明，终端约束在非标称条件下依然能够得到满足，同时飞行路径角值保持在0°附近，满足拟平衡滑翔条件和解析解中的前提假设。图12给出了单次仿真中各过程变量的最大值，可以看出在500次仿真中，过程约束均得到满足。

图10 滑翔段飞行剖面

图11 飞行路径角随时间变化情况

最后指出，仿真中取更新周期ΔRcov=50 km；由表1中的初始速度可知，ΔRcov对应的飞行时间大于10 s。因此，路径点更新周期(5ΔRcov)大于50 s，而飞行剖面优化周期(4ΔRcov)大于40 s。基于上述剖面设计方法和解析解，该更新周期在仿真中足够完成在线剖面优化以及三维滑翔轨迹重构。