基于虚拟再入角的快速离轨制动制导方法

2023-04-02权申明王竹晁涛杨明

兵工学报 2023年3期

权申明 , 王竹晁涛杨明

(1. 哈尔滨工业大学控制与仿真中心，黑龙江哈尔滨150080; 2. 上海机电工程研究所，上海 201109)

0 引言

飞行器在再入点的状态如再入速度倾角、再入位置、再入速度等，对飞行器在大气层内的飞行状态和飞行能力有至关重要的影响[1-2]。良好的再入状态，对飞行器在大气层内的驻点热流密度、总加热量等物理约束、滑翔飞行的调节能力产生积极的作用，同时能够增加飞行器的安全性和稳定性，减轻制导控制系统的压力。再入点各状态约束值由飞行器设计的总体单位根据任务需求及飞行器性能计算得到，设计离轨段制导控制算法时可认为是已知的。

飞行器离轨制动段制导算法可以规划出一条满足再入点状态的转移轨道。为方便研究，在初步设计轨道时，通常假设飞行器在轨道转移过程中可瞬间获得所需的速度冲量且仅受地球中心引力场作用。

文献[3]指出，良好的离轨制动飞行状态是确保天地往返飞行器安全、准确返回既定着陆场区的前提。传统钝头体再入返回舱，受到横向机动能力和搜救能力限制，必须保证再入点有足够的控制精度。离轨制动可以视为轨道转移的一种形式，现有的轨道转移轨迹规划算法设计时，发动机主要有脉冲推力(视为瞬时脉冲)、小推力和有限推力三种形式。其中，基于脉冲推力的计算较为简单，相对成熟的Lambert 制导算法可在较短时间内计算出变轨点的速度增量，然而实际飞行过程中，发动机难以提供任意大小的发动机推力；小推力轨道转移通常用在深空探测等场景，推力较小，因此需要数日才能完成某任务；有限推力离轨制动形式在航天器返回、载荷释放等场景下应用广泛，对该问题的研究具有工程意义。具体地，有限推力离轨规划又可分为逼近式策略以及优化式策略。

在逼近式策略方面，文献[4]提出两次“推-滑”的离轨制导策略，分析了推力沿不同方向施加时能量和动量矩的变化规律，给出临界高度的计算公式及两次“推-滑”的判断条件。文献[5]将离轨制动过程分成能量耗散段和闭路导引控制段，探讨了在燃料随机耗尽情况下推力方向对再入点参数的影响，分析了推力方向切换与能量窗口的关系。文献[6]利用轨道瞬时根数进行离轨制动算法设计，将离轨制动问题转化为制动起点的纬度幅角的单参数搜索问题，并考虑了J2 项摄动、初始状态偏差、质量与推力误差的影响。文献[7]借鉴运载火箭能量管理的思路，为实现耗尽关机，将离轨制动过程分为能量耗散与闭环制导控制两个阶段，在运算效率与制导精度方面均有一定的提升。逼近式策略方面的研究基于机理分析，从能量角度设计制动方案，通常仅考虑再入点速度大小和方向约束，由于没有考虑时间代价，传统的逼近式策略难以满足快速离轨需求。

在优化式策略方面，文献[8]提出的方法能够实现指定再入点位置的飞行器离轨任务，但在内层计算过程中需要反复调用求解非线性规划求解器，因此难以实现制动速度的在线计算。文献[9]考虑再入返回飞行器对再入点高度和再入角的要求，推算了离轨制动的轨道参数及所需速度增量，同时指出了惯性系速度倾角和地固系下的再入角差异较小。文献[10]针对建立了有限推力转移轨迹优化的凸优化模型，并提出一种迭代逼近算法，但每步迭代计算中的凸优化求解耗时在1.24~3.23 s 之间，算法难以实现在线优化。文献[11-12]针对燃料最优问题进行研究，在飞行器总体优化设计阶段具有重要意义，能够节约大量成本。文献[13]基于最优控制理论推导了离轨机动的最优性条件，给出了协态变量和状态变量的解析解，进而转化为两点边值问题，但初值的选取对求解结果有较大影响。文献[14]针对转移时间最优和燃料最优两种性能指标，将原问题转化为非线性规划问题，将Bezier 曲线法得到的结果作为初值，采用内点法与序列二次规划算法进行求解，虽然对初值的计算进行了改进与优化，但是仍难以做到在线计算。文献[15]分别针对共面、异面轨道转移问题，研究了基于凸优化方法的轨迹规划，并选用不同指标进行了轨道机动能力分析，但是算法实时性较差，难以工程实现。文献[16]提出一种组合制导方法，前期采用燃料最优的开环策略，采用基于标称轨迹的闭环制导，受开环控制的影响，该方法的调节能力有待进一步研究。文献[17]针对动能弹垂直再入问题利用自适应伪谱法对有限推力制动与连续小推力制动两种方式进行了研究。上述研究成果往往需要借助工具箱，如SNOPT、SQP、CVX 等[18]。

由以上分析可以看出，基于逼近式和优化式策略的离轨制动算法均存在一定的不足：逼近式策略算法简单，计算效率高，易于工程实现，然而少有考虑时间最优的离轨需求；优化式策略可以求解同时考虑多个再入状态约束下的快速离轨制动问题，然而优化计算效率较低，难以实时在线计算。随着空间飞行器智能化需求的不断提高，再入飞行任务已经从到达地面固定目标点，逐步向着覆盖全球任意位置发展，因此传统的固定点离轨制动、固定的再入难以适应这种需求。对于已经部署就位的航天器，其上所携带燃料固定，为了应对一定范围内可变的任务要求，飞行器需要具备一定的再入点状态调节能力。

本文研究在离轨制动任务确定之后，燃料有限情况下的多终端约束的时间最优离轨制动算法。本文的主要创新点如下：

1)结合优化算法的结果，通过其最优轨道形式，提出燃料约束情况下的有限推力快速离轨制动制导算法；

2)提出虚拟再入角的概念，将时间最优快速离轨制动问题参数化为虚拟再入角的求解，简化优化问题；

3)为满足离轨制导算法的实时性要求，本文所提出的离轨制动算法计算简单、快速且具有较强的收敛性。

1 问题描述

飞行器轨道转移示意图如图1所示。图1中，M为飞行器接到任务时的位置，K为离轨制动初始位置，E为再人点位置，C为地球表面目标点，h为轨道高度，hf为离轨段终端高度。选取轨道转移时间固定或能量最小等为指标后，可唯一确定出转移轨道。

图1 轨道转移示意图Fig. 1 Schematic diagram of orbital transfer

由于推力有限且连续作用，飞行状态不能根据自由段解析弹道的特征解析计算，而是需要通过积分弹道方程来获得。飞行器在大气层外的运动模型可以描述为

式中：r为飞行器位置矢量，r为飞行器质心距地心的径向距离，r=‖r‖；v为飞行器速度矢量；µ为地球引力常数；pe为发动机推力，pe表示制动推力大小；m为飞行器质量；Isp为发动机比冲；g0为地球表面重力加速度大小。

为了更加清晰地描述飞行器相对于地球的运动情况，给出惯性系下飞行器运动模型的另一种表示：

式中：v为飞行器相对地速；θ为飞行器速度倾角，在上为正；λ和φ分别为飞行器星下点的经度和纬度；χ为弹道偏角，顺时针为正；α表示制动角，即推力方向与速度方向的夹角(类似于攻角的定义)；g为重力加速度，g=g0(R0/r)2，R0为地球平均半径。

2 离轨制动算法设计

优化式策略通常基于最优控制相关理论，求解复杂、效率低，往往需要调用优化工具箱得到最优结果；逼近式策略基于能量守恒等原理，结合卫星轨道理论，实时计算可行的控制量，逐渐调整轨迹，可逐步逼近期望终端状态。

2.1 冲量假设的离轨制动算法

冲量假设简化了飞行器离轨制动的运动过程，认为飞行器可瞬间获得需要速度，给离轨制动段的初步设计带来了方便。

针对更加一般性离轨制动问题，飞行器一旦进入离轨制动阶段，为实现高度的降低，飞行轨迹为椭圆轨道。

设v1为当前位置需用速度，r1表示由地心指向飞行器实时位置的向量，r2表示由地心指向期望再入点位置的向量，h0为飞行器离轨点高度，θ1为根据能量守恒和动量矩不变原理反求出r1处的期望速度倾角。冲量假设下的离轨制动问题描述为：已知r1处的初始速度v0和初始速度倾角θ0，r2处的终端速度v2和终端速度倾角θf，求解冲量假设下的离轨段制动参数Δθ和Δv，v1和制动速度的夹角为φ。矢量r1和r2的大小分别为r1和r2，v1和v2的大小分别为v1和v2。图2 为椭圆轨道离轨段航迹示意图。

图2 椭圆轨道离轨段航迹示意图Fig. 2 Deorbit trajectory of an elliptical orbit

下面针对图2 所示的任意椭圆轨道进行算法设计。根据能量守恒和动量矩不变原理，有

通过式(3)可以求解出θ1和v1，代入下式中：

由此可以得到离轨段制动参数Δθ和Δv，若能够瞬时沿着φ的方向提供一个Δv，则飞行器将沿着新轨道运动，最终到达再入点。

2.2 有限推力离轨制动算法

通常情况下，天基再入滑翔飞行器的离轨制动发动机推力有限，因此2.1 节的冲量假设离轨制动制导方法难以适用。此外，即使发动机可瞬时提供任意方向及大小的推力，由于空间中摄动因素的存在，其离轨段终端将产生较大偏差，影响再入飞行。因此，采用有限推力方式的离轨制动制导方法研究更有实际意义，脉冲假设下的计算结果可作为参考。

由于燃料的消耗，飞行器质量不断减少，根据给定的发动机参数，可以利用推力和比冲进行计算。

在有限推力模型中，采用冲量假设下的离轨制动制导方法，可以求解得到冲量假设下的离轨段离轨参数Δθ和Δv。在求解有限推力下的离轨段轨迹时，给定天基再入滑翔飞行器的制动点状态及再入点状态约束，计算推力大小与方向变化曲线，控制飞行器按照期望再入角、再入速度到达再入点高度。

结合2.1 节瞬时脉冲仿真结果，设计一种采用逼近式策略的考虑有限推力的制导步骤如下：

步骤1根据飞行器实时状态、再入角度及速度，由式(3)按照瞬时脉冲方式确定出当前时刻所需的速度和角度信息。

步骤2结合飞行器当前角度和速度信息，计算制动速度Δv及制动角度。

步骤3沿着制动角度方向，改变飞行器速度增量为Δv，发动机采用最大推力。

步骤4若所需速度小于其设定阈值，即，则关闭发动机，飞行器沿椭圆轨道自由运动，否则循环步骤1~步骤3 的计算。

2.3 快速离轨制动算法设计

通过分析不难发现，以上两种方法没有考虑时间最优性指标，仅考虑了再入点高度、速度大小与速度倾角，因此求得的离轨制动轨迹难以满足快速离轨的需求。

通过最优控制理论中的时间最优结论可知，时间最优情况下控制量通常是Bang-Bang 控制，类似于开关控制量，离轨制动阶段飞行器所携带的燃料有限，因此假设时间最优的离轨制动问题中，制动发动机推力大小与制动角度形式如图3 所示。图3中，O点为原点，α轴表示制动角，t轴为时间，α1(t)为0~t1时间内的制动角变化曲线，α2(t)为t2~t3时间内的制动角变化曲线，pemax为最大制动推力，t1为制动发动机首次关机时刻，t2为第2 次开机时刻，t3为第2 次关机时刻，t4为到达再入点的时刻。

图3 时间最优制动发动机工作状态示意图Fig. 3 Schematic diagram of time-optimal deorbit

求解α(t)，使得终端时间tf最小，同时满足初始各状态与终端各状态，可描述为

式中：u为制动角变化曲线α(t)；Δm为飞行器总质量变化量；mfuel为燃料质量。

可将式(5)描述的问题转化为非线性规划问题，进而基于伪谱法进行求解，但是该类方法计算效率较差，工程实现性不强。还可将该非线性规划问题转化为凸优化问题，凸优化在求解最优控制问题的快速性与最优性上均有显著提升[19]，然而现阶段凸优化算法仍难以满足在线实时计算的速度需求。

由于传统的有限推力离轨制动算法并未考虑时间最优性，存在一定不足。设计改进的有限推力离轨制动算法如图 4 虚线曲线所示，引入虚拟再入角的概念，在前期采用虚拟再入角，后期采用实际期望再入角作为约束。

图4 改进方法与传统方法弹道倾角变化曲线示意图Fig. 4 Curve of the flight path angle with the improved method and the traditional method

定义1虚拟再入角：

在再入返回的初期，设计不同于实际期望终端速度倾角θf的某固定或者时变再入角作为制导输入指令，称角度为虚拟再入角(在后期仍使用θf作为期望再入角)。

在离轨制动前期主要控制再入速度倾角与再入速度大小，后期根据文献[4]中得到的垂直速度方向施加力不改变能量大小的结论，通过预测剩余燃料消耗时间，进行正推与反推的能量耗散策略。后期使用剩余燃料朝着目标点飞行，由于后期时间相对短，终端速度不会偏差太大。本文设计的快速离轨制动制导仿真流程如图 5 所示。

图5 快速离轨制动制导仿真流程图Fig. 5 Flow chart of simulation of fast deorbit guidance

因此，可将问题式(5)转化为如下问题：

并使用凸优化方法对关键参数 Δt1、Δt2，以及第一阶段的虚拟期望终端倾角求解。

基于冲量假设的有限推力下离轨制动算法存在的问题在于没有考虑快速离轨的需求，但在再入速度倾角和速度大小的控制上，计算量小且精度较高。基于凸优化的离轨制动可以考虑时间最优约束，仍难以实时计算，但是通过与有限推力离轨制动算法相结合，将最优结果作为有限推力制动方法虚拟终端约束，从而给出快速离轨制动制导算法。

2.4 基于凸优化方法的子问题求解

文献[19]给出了一种可以处理非线性等式的凸优化问题的方法。在时间最优的离轨制动问题中，给定起始时刻的高度、速度，要求飞行器以最短时间到达某一高度，并满足终端的速度大小及倾角要求。

使用凸优化对问题进行求解需要分为三个步骤。

2.4.1 状态方程归一化

式中：σ为松弛变量，代表推力-加速度。设状态量为x=(r,v,z)T，控制量为u=(τ,σ)T，将状态方程写为如下矩阵形式：

当使用迭代方法求解最优解时，式(10)可看作线性方程。

2.4.2 求解初始轨迹

使用迭代方法求解时间最优离轨制动问题时，首先需要一条初始轨迹，在此基础进行优化。

求解时间最优离轨问题，需要将离轨过程分为N段，每一段时间间隔为t0/N，求解N+ 1个状态变量及控制变量。飞行器初始位置、速度与终端时刻位置、速度约束可表示为

离轨过程约束为式(9)~式(13)，由于地心距r为归一化后的变量，在离轨过程中近似为1，在对结果精确度要求不高的情况下，式(9)中的r可近似为1，变为线性方程。

由于所有约束及目标函数均为凸的，此问题为一标准的凸优化问题。

2.4.3 迭代求解时间最优离轨制动问题

求解时间最优离轨制动问题不指定固定的再入点以及离轨时间，只要求到达固定的再入高度时满足指定的速度大小及倾角，因此需要引入时间变量。此时式(9)与式(13)均为非线性等式约束。同时，为保证结果的精确，需要在初始解的基础上进行多次迭代求解。

当凸优化问题中含有非线性等式约束hi(y) = 0,i= 1，… ,q时，假设当前迭代结果为y[k]，y= （f x,u,tf）。求解y[k+1]的过程如下：

在y[k]处对函数进行泰勒展开，hi(y) =0可以写为

同时，又有展开式

凸优化问题的约束为

优化结果为yp，将在y[k]处1阶展开：

式(18)两侧同乘(yp-y[k])T/2，得到

因此便将非凸约束hi(y) =0转换为凸约束。

时间最优的离轨制动问题中，除了考虑式(9)~式(11)的约束之外，其余约束为

以离轨段终端时间tf最短为指标。式(9)的离散形式以及式(22)为非线性约束，将其按前述可得到凸函数的约束条件，并附加约束

即可求得y[k+1]以及u[k+1]，多次迭代直至收敛，得到时间最优离轨制动结果。

至此，将离轨制动问题转化为凸优化问题，使用2阶锥规划方法进行求解，求解工具详见文献[20]。

2.5 基于最优轨迹的虚拟指令计算

根据凸优化求解得到的优轨迹中发动机推力大小曲线，可以看出时间最优的离轨制动问题中，控制量呈现“推-滑-推”的特点，从而根据其时间变化规律可得 Δt1、 Δt2，根据速度、倾角变化曲线，根据式(4)可以反求出虚拟再入角θ˜f大小。

3 数值仿真及分析

考虑有两台推力发动机，每台推力为500 N，比冲为300 s，离轨到再入阶段单独计算燃料质量，燃料为250 kg。

根据2.1~2.3 节给出的不同方法，考虑从某圆轨道进行离轨制动，初始高度h0=500 km，初始速度为v0=7 616.7 m/s，初始倾角θ0= 0°，期望终端高度hf=120 km、终端速度vf=7 950 m/s、终端倾角θf=-2 .5°为例，分别给出不同方法的仿真结果。需要指出的是，本文所提出的快速离轨制动制导方法采用类似的策略在横向进行虚拟指令设计后，可适用于三维空间中的异面离轨制动，仿真中仅以纵向运动为例。

3.1 基于凸优化仿真结果

将快速离轨制动问题转化为凸优化问题，进而使用2阶锥规划进行求解。仿真结果如图6所示。通过多约束下时间最优离轨制动的仿真计算，可以看出飞行轨迹能够实现时间最优的指标；从推力大小曲线可以看出，发动机工作状态类似Bang-Bang控制，与最优控制的时间最优控制问题结果一致；由图6(f)可以看出，飞行器终端时刻燃料用尽，终端各状态满足要求。同时，为了验证方法有效性，放宽燃料限制，当燃料足够时，推力可以长时间最大推力工作，通过改变制动角来实现时间最优指标。

图6 多约束下时间最优离轨制动状态量变化曲线Fig. 6 Curves of time-optimal deorbit with multiple constraints

综合以上仿真曲线可以看出，时间最优下的仿真结果虽比有限推力离轨制动的制动时间指标更优，但是从编程实现及运算效率上来看，代价较大。即使在初值合理的情况下，单次优化耗时在 1 s 以上，因此难以实现在线制导的目标。

3.2 快速离轨制动算法仿真

为验证本文算法的有效性，选取不同终端状态θf=-1 .5° 、θf=-2 .5° 、θf=-3 .5° 、θf=- 4.5° 进行仿真，基于文献[6]中在线制导算法的仿真结果如图 7 所示，基于本文快速离轨制动制导算法仿真结果如图 8 所示。

对比图7、图8仿真结果曲线可以看出，两种方法均能满足终端高度、速度、倾角的终端约束。在传统算法中，未考虑快速离轨制动需求时，燃料有剩余，且随着终端倾角(幅值)的增大，离轨制动段飞行时间减小，燃料消耗增加，这是由于空间轨道中，速度方向的改变需要较大的速度增量。对比以上图7(f)、图8(f)子图燃料质量变化曲线可知，为满足离轨制动快速性的需求，本文算法能够将燃料耗尽，离轨段节约时间最大值约为500 s。在计算效率方面，单步计算耗时均在5 ms以内，能够满足在线制导的需求。

图7 传统离轨制动制导算法仿真结果Fig. 7 Simulation results of the traditional deorbit guidance algorithm

图8 快速离轨制动制导算法仿真结果Fig. 8 Simulation results of the fast deorbit guidance algorithm

3.3 蒙特卡洛仿真

在实际离轨制动过程中，飞行器发动机受到各种偏差因素的影响，本文主要考虑推力大小摄动Δpe、推力方向摄动Δφ、燃料质量摄动Δm，各偏差均为均匀分布如表1 所示。

表1 参数摄动表Table 1 Parameter perturbations

进行1 000 组蒙特卡洛仿真实验，仿真结果如图 9 所示(Rf表示离轨制动段射程)。

图9 蒙特卡洛仿真参数散布图Fig. 9 Distribution of Monte Carlo simulation parameters

根据终端状态散布结果可以看出，该算法对终端速度与终端倾角控制精度较高，其中终端速度偏差为±0.2 m/s，终端倾角偏差为±0.002°。由于该算法未对射程进行精确控制，该状态的绝对偏差较大为±50 km，然而该值相对于整个再入阶段相对误差较小，约为0.5%，同时考虑到后续的再入滑翔阶段，对射程调节能力较强，因此离轨制动阶段射程的影响可忽略。