火星探测最优离轨制导方法研究

滕 锐，韩宏伟，乔 栋

（1. 中国运载火箭技术研究院 空间物理重点实验室，北京 100076；2. 北京理工大学 宇航学院，北京 100081；3. 深空探测自主导航与控制工信部重点实验室，北京 100081）

引 言

在火星着陆探测过程中，EDL（Entry，Descent，and Landing）过程的成功与否直接决定了整个任务的成败，而EDL过程前的离轨制动又是决定EDL过程顺利执行的关键环节，因此其在探测器着陆过程中起到决定性的作用。离轨机动是飞行器从初始环绕轨道至行星大气边缘的制动转移过程，是飞行器在进入行星大气前执行的重要机动过程。由于大气进入始于大气边缘，并且进入过程的初始状态对于整个进入乃至着陆过程均会产生较大的影响，因此，行星大气进入前需要保证理想的进入初始状态[1]。大气进入的初始状态主要通过离轨机动来实现，离轨机动过程的性能直接决定了进入初始状态的精度，因此，高效且精确的离轨机动策略是保证行星大气进入过程顺利进行的前提[2-4]。由于离轨制动问题在整个火星探测过程中举足轻重，因此，离轨制动相关的轨道设计和制导控制等关键技术已被美国等航天强国列为航天技术中重要的一类技术研究[5]。

对于离轨制动问题，国内外很多学者对离轨策略进行了诸多的研究。早期的研究主要将离轨过程看成一种冲量式制动问题，也即在给定终端大气进入接口状态约束和最优性能指标的前提下，给出单脉冲和多脉冲离轨制动问题的最优解[6-9]。为保证解的存在性和求解过程的收敛性，此类问题一般将终端大气进入状态约束设定为范围约束，因此不同目标函数下最优解对应的大气进入状态存在差异[10]。考虑到实际任务执行中，飞行器发动机需执行多弧段有限推力机动，从而实现离轨制动过程，因此，后期的研究以实际机动形式对应的问题为研究对象。对于此类问题，目前的研究均以固定的大气进入速度和飞行路径角为终端约束，通过设定离轨机动弧段和滑行弧段的组合形式，在推导得到各个时间节点解析的状态变量和协态变量的基础上，采用多重打靶得到离轨制导轨迹[11]。该制导方法的优势在于不仅保证了方法的实际可操作性，又保证了推力方向的快速生成，因此该方法目前已成为最具优势的离轨制导策略[12-13]。

本文以火星进入离轨机动问题为研究对象，以大气入口状态为离轨机动终端约束，给出了一种 “推–滑”两弧段离轨机动模式的最优离轨制导方法，此外通过蒙特卡罗仿真验证了方法鲁棒性，并分析了离轨机动过程的随机偏差对火星大气入口状态精度的影响。

1 最优有限推力离轨制导问题描述

在火星探测过程中，当飞行器从环绕停泊轨道至进入大气的过程之间，需执行连续推力机动实现离轨过程，因此，整个离轨过程从问题本质上属于多重弧段（即包括多段推力弧和滑行弧）有限推力轨道机动问题。本文考虑到实际任务时间限定和有利实际工程应用等客观因素，以一次“推–滑”两弧段转移形式为离轨过程飞行器的机动模式，并在此基础上给出最优制导方法。

尽管离轨机动本质属于有限推力机动过程，但部分约束形式等方面又具有其独特的一面。由于一般进入大气前的停泊轨道高度较低，因此，为保证制导算法在进行离轨轨迹求解时的实时性要求，其动力学方程的形式、状态变量以及协状态变量的求解等均会在一定的精度要求和假设前提下进行简化，所以动力学模型以及问题的描述形式也不同于一般的有限推力轨道机动问题。

1.1 动力学模型

通过火星惯性系下的动力学模型来描述离轨机动问题。其对应的无量纲动力学方程为

其中：r为火心惯性系下的无量纲位置矢量；V为火心惯性系下的无量纲速度矢量；m为有量纲的飞行器质量；T为飞行器发动机推力，Isp为飞行器发动机比冲；g0为火星表面引力加速度；υ 为推力方向在火心惯性系下的单位矢量。在无量纲化过程中，距离的无量纲参数为火星半径R；速度的无量纲参数为；时间的无量纲参数为。

由于大气进入前的停泊轨道高度较低，且位置矢径大小对应的高度介于停泊轨道高度与大气边缘高度之间，因此，整个离轨过程位置矢径大小变化微小。基于该特性，在每次迭代制导的过程中，位置矢径可近似为[1]

其中：r0为每次制导迭代初始时刻的无量纲位置矢径大小。因此，式（1）中动力学方程中速度矢量的微分项可重新表示为[1]

1.2 离轨机动最优性条件

基于式（1）和式（3）的动力学模型，在给出最优性条件之前，首先给出协状态方程的表达式。假设速度矢量、位置矢量以及质量对应的协状态量分别为pr，pV和pm，基于最优控制理论[14]，系统的哈密顿函数为

其中，协状态量pr和pV对应的微分方程为

在离轨过程中，由于“推–滑”离轨机动模式下，飞行器发动机机动模式固定，所以此处控制变量只有“推”弧段的推力方向单位矢量。根据极大值原理，最优离轨过程对应的发动机推力方向单位矢量υ 为

上述式（6）中给出的推力方向单位矢量 υ∗是最优离轨机动过程的必要条件。

1.3 协态变量和状态变量的解析解

在如式（2）的假设前提下，离轨机动的动力学方程中部分项被线性化，而式（5）中的协状态微分方程直接具有线性微分方程的基本形式，这种特性给微分方程的解析求解提供了一定的便利。因此，本文采用一种将协态变量和状态变量变换的方法，求解得到状态方程和协态方程的解析解[1-2,12]。事实上，这种解析求解方法目前已经成熟应用于离轨制导和入轨制导等问题中。下面将给出该方法求解状态和协态变量解析解的具体过程。

首先，对协态变量进行变化，即引入一个新的量为

对于新变量，此处省略推导过程，直接给出对应的从式（5）变形得到的微分方程的解析解，即

其中：t和t0分别是当前时间和初始时间，I3为3阶单位矩阵，λ (t0)为初始时刻对应的值。

其次，对于状态变量，此处同样需进行变换从而使得状态变量的微分方程可以解析得到，引入一个新的量为

根据参考文献[12]中的方法，在给出如式（9）的变换之后，从状态变量转换得到的新变量对应的微分方程具有如式（10）的解析解形式，即

其中： Θ (t0)为初始时刻对应的值，为

此外，式（10）中的向量 Γ (t,t0)的表达式为两个子向量的组合，其具体表达式为

其中，子向量的表达式分别为

对于式（13）和式（14）中子向量的值，在每一时刻均采用数值积分求解，一般情况下，只要数据节点之间的间隔足够小，其积分精度能够得到保证。为了保证求解速度的要求，两个时间节点间隔之间只需进行一次积分即可，数值积分方法诸如4阶龙格库塔法以及米勒法等均能够以较高精度实现对上述积分函数的求解。具体求解公式此处不再赘述。

1.4 离轨过程的终端约束及两点边值问题

1.3小节给出了状态变量和协态变量的解析计算方法，对于离轨制导问题，从问题本质上属于典型的两点边值问题，因此，下面需要分别确定两点边值问题对应的变量及终端约束。对于离轨机动过程，一般初始状态是确定的，初始协态变量和总的离轨时间是可变的，且终端状态受约束。在离轨制导过程中，制导参数包括初始协态变量和总离轨时间，但是由于终端状态变量受约束，因此直接根据变分法得到的制导终端约束会包含拉格朗日乘子，因此下面将推导给出只包含状态变量和协态变量的制导终端约束。

对于火星离轨机动过程，一般包括如下3个终端状态约束，即

其中：rf和Vf分别是终端位置矢量和速度矢量，而rEI、VEI和 γEI分别为基于上述3个终端状态约束。对于最优离轨过程，其横截条件为

其中：tf是终端时间；µ1、µ2和µ3是拉格朗日乘子。

为了消去式（18）中的拉格朗日乘子，下面将结合式（15）～（17）与式（18）给出只包含状态和协状态的制导终端约束。首先对式（18）等式左右分别转置之后右乘位置矢量rf和速度矢量Vf，从而引入式（15）～（17）的等式关系，也即

通过式（19）中的3个公式，可以直接求解得到3个拉格朗日乘子 µ1、µ2和µ3，然后将其带入式（18），可得到不含拉格朗日乘子的横截条件，其最终终端约束方程为

其中，中间矩阵A为

式（20）即为最优“推–滑”离轨制导对应的由状态变量和协态变量构成的终端约束。此外，由于终端时间自由，因此，需要额外增加终端哈密顿函数为0的约束，即

至此，制导过程对应的变量和约束方程已给出。

在火星离轨制导过程中，为保证收敛性和求解速度每个制导周期内问题的求解均采用莱文贝格–马夸特（Levenberg-Marquardt）方法[15]。一般情况下，制导变量初值选取的优劣对首次制导方程求解的收敛速度有明显影响，因此，此处基于实际数值仿真结果，给出较为一般性得制导变量初值，即协态变量初值此外，离轨总时间t的初f值选取对收敛速度影响不大，一般选取为初始轨道周期的一般较为合理。下文将给出离轨制导迭代流程。

2 仿真结果与分析

为验证上述最优“推–滑”离轨制导方法的有效性和鲁棒性，本节以火星探测任务为背景，给出相关的仿真结果。首先，给定火星探测器在进入大气前所在的停泊轨道根数，其数值通过转换即可得到离轨机动的初始状态变量，如表1所示。

表1 离轨制导的初始轨道根数Table 1 Initial orbit elements of deorbit guidance

此外，给定火星半径R= 3 396 km，因此，初始轨道实际上为高度为350 km的圆轨道。在离轨机动过程中，飞行器的发动机推力T= 1 000 N，发动机比冲Isp=311 s，飞行器初始质量假设为4 000 kg。由于整个离轨过程以飞行器到达大气进入边缘为终止位置，因此，此处给定大气边缘高度为128 km。此外，所有仿真均假定目标大气入口处速度约束为VEI= 3 516.2 m/s，飞行路径角γEI= –2.358°，而进入点的位置不进行约束。

2.1 火星离轨机动最优解

基于上述离轨过程状态变量初值和发动机参数，本节给出了不考虑随机偏差影响的火星“推–滑”离轨机动最优解。图1给出了最优“推–滑”离轨机动对应的转移轨迹。从图中可以看出，整个过程机动弧段较短，最终实现的是类似于一个脉冲机动的效果，而滑行弧段较长。此外，从转移轨迹也能看出，在终端速度和飞行路径角约束下，整个机动几乎在同一轨道面内。

图1 最优离轨过程轨迹Fig. 1 Optimal deorbit trajectory

通过仿真可知，整个“推–滑”离轨机动最优解对应的离轨总时间为2 241.46 s，其中机动弧段总时间为295.31 s，而滑行弧段历时1 946.15 s。因此通过发动机参数推算，整个离轨机动所消耗的燃料为255.696 kg，对应的速度增量为76.29 m/s。对于机动弧段，图2给出了该阶段对应的发动机推力单位方向矢量在惯性系下的分布。为了更为直观展示推力方向与速度的夹角分布，图2还给出了推力方向角，此处推力方向角定义为推力方向与飞行器速度方向的夹角，从结果可以看出，整个离轨过程，推力几乎近似沿速度的反方向，也即减速制动从而使得飞行器降轨进入大气。此外，通过仿真计算，表2给出了最优解对应的初始协态变量的值。

图2 推力单位方向矢量分量与推力方向角变化曲线Fig. 2 Components of Thrust Unit Vector and Thrust Angle

表2 最优离轨过程对应的协态变量初值Table 2 Initial costate value of optimal deorbit process

2.2 蒙特卡罗仿真分析

为了验证上述最优“推–滑”离轨制导方法的鲁棒性和精度，本节将在不确定性随机偏差存在的基础上进行蒙特卡罗仿真，并分析不确定性因素的存在对“推–滑”离轨制导过程的影响。在制导过程中，给定机动弧段制导周期为5 s，在每一次制导起始增加随机偏差，根据工程经验，此处给定位置偏差的3σ值为10 m，速度偏差的3σ值为0.1 m/s，蒙特卡罗打靶次数为1 000。

图3 最优离轨过程轨迹蒙特卡罗仿真分布Fig. 3 Monte Carlo simulation of optimal deorbit trajectories

图4 离轨终端状态及总离轨时间蒙特卡罗仿真散布Fig. 4 Dispersions of deorbit terminal state and total deorbit time

图3给出了蒙特卡罗仿真得到的最优离轨过程轨迹分布，从轨迹图可以看出尽管随机偏差存在，但整个轨迹散布范围很小。图4分别给出了终端状态和离轨总时间在随机偏差存在下的散布情况。从图中可以看出，进入点路径角的散布与标称值 γEI相差最大约为0.046°，而进入点速度的散布与标称值VEI相差最大约为0.5 m/s，因此，整个制导结果精度满足进入点状态偏差的精度的要求。从离轨总时间散布来看，其受随机偏差的影响较为明显，但其与最优解的值2 241.46 s相比最大偏差仅在8 s左右。此外，为反映机动过程推力方向的分布，图5给出了推力单位方向矢量分量与推力方向角在随机偏差影响下的散布情况，可以看出尽管随机偏差存在，但推力方向仍处于标称推力方向附近。

图5 推力单位方向矢量分量与推力方向角的蒙特卡罗仿真分布Fig. 5 Distributions of Thrust Unit Vector components and Thrust Angle

3 结束语

针对火星探测中的离轨机动问题，本文给出了一种基于解析状态变量和协态变量的最优“推–滑”离轨制导方法。该方法在位置矢径局部线性化的假设下分别给出了最优“推–滑”离轨机动的必要条件和只由状态变量和协态变量构成的制导终端约束。通过仿真分析发现，该方法所得到的离轨过程近似为机动，且机动弧段的推力方向几乎沿速度反方向。蒙特卡罗仿真表明该方法在满足制导精度的前提下具有较强的鲁棒性。上述研究结果可为我国未来火星探测离轨机动方案设计提供技术支持和理论依据。