火星探測(cè)最優(yōu)離軌制導(dǎo)方法研究

滕 銳，韓宏偉，喬 棟

（1. 中國(guó)運(yùn)載火箭技術(shù)研究院 空間物理重點(diǎn)實(shí)驗(yàn)室，北京 100076；2. 北京理工大學(xué) 宇航學(xué)院，北京 100081；3. 深空探測(cè)自主導(dǎo)航與控制工信部重點(diǎn)實(shí)驗(yàn)室，北京 100081）

引 言

在火星著陸探測(cè)過程中，EDL（Entry，Descent，and Landing）過程的成功與否直接決定了整個(gè)任務(wù)的成敗，而EDL過程前的離軌制動(dòng)又是決定EDL過程順利執(zhí)行的關(guān)鍵環(huán)節(jié)，因此其在探測(cè)器著陸過程中起到?jīng)Q定性的作用。離軌機(jī)動(dòng)是飛行器從初始環(huán)繞軌道至行星大氣邊緣的制動(dòng)轉(zhuǎn)移過程，是飛行器在進(jìn)入行星大氣前執(zhí)行的重要機(jī)動(dòng)過程。由于大氣進(jìn)入始于大氣邊緣，并且進(jìn)入過程的初始狀態(tài)對(duì)于整個(gè)進(jìn)入乃至著陸過程均會(huì)產(chǎn)生較大的影響，因此，行星大氣進(jìn)入前需要保證理想的進(jìn)入初始狀態(tài)[1]。大氣進(jìn)入的初始狀態(tài)主要通過離軌機(jī)動(dòng)來實(shí)現(xiàn)，離軌機(jī)動(dòng)過程的性能直接決定了進(jìn)入初始狀態(tài)的精度，因此，高效且精確的離軌機(jī)動(dòng)策略是保證行星大氣進(jìn)入過程順利進(jìn)行的前提[2-4]。由于離軌制動(dòng)問題在整個(gè)火星探測(cè)過程中舉足輕重，因此，離軌制動(dòng)相關(guān)的軌道設(shè)計(jì)和制導(dǎo)控制等關(guān)鍵技術(shù)已被美國(guó)等航天強(qiáng)國(guó)列為航天技術(shù)中重要的一類技術(shù)研究[5]。

對(duì)于離軌制動(dòng)問題，國(guó)內(nèi)外很多學(xué)者對(duì)離軌策略進(jìn)行了諸多的研究。早期的研究主要將離軌過程看成一種沖量式制動(dòng)問題，也即在給定終端大氣進(jìn)入接口狀態(tài)約束和最優(yōu)性能指標(biāo)的前提下，給出單脈沖和多脈沖離軌制動(dòng)問題的最優(yōu)解[6-9]。為保證解的存在性和求解過程的收斂性，此類問題一般將終端大氣進(jìn)入狀態(tài)約束設(shè)定為范圍約束，因此不同目標(biāo)函數(shù)下最優(yōu)解對(duì)應(yīng)的大氣進(jìn)入狀態(tài)存在差異[10]?？紤]到實(shí)際任務(wù)執(zhí)行中，飛行器發(fā)動(dòng)機(jī)需執(zhí)行多弧段有限推力機(jī)動(dòng)，從而實(shí)現(xiàn)離軌制動(dòng)過程，因此，后期的研究以實(shí)際機(jī)動(dòng)形式對(duì)應(yīng)的問題為研究對(duì)象。對(duì)于此類問題，目前的研究均以固定的大氣進(jìn)入速度和飛行路徑角為終端約束，通過設(shè)定離軌機(jī)動(dòng)弧段和滑行弧段的組合形式，在推導(dǎo)得到各個(gè)時(shí)間節(jié)點(diǎn)解析的狀態(tài)變量和協(xié)態(tài)變量的基礎(chǔ)上，采用多重打靶得到離軌制導(dǎo)軌跡[11]。該制導(dǎo)方法的優(yōu)勢(shì)在于不僅保證了方法的實(shí)際可操作性，又保證了推力方向的快速生成，因此該方法目前已成為最具優(yōu)勢(shì)的離軌制導(dǎo)策略[12-13]。

本文以火星進(jìn)入離軌機(jī)動(dòng)問題為研究對(duì)象，以大氣入口狀態(tài)為離軌機(jī)動(dòng)終端約束，給出了一種 “推–滑”兩弧段離軌機(jī)動(dòng)模式的最優(yōu)離軌制導(dǎo)方法，此外通過蒙特卡羅仿真驗(yàn)證了方法魯棒性，并分析了離軌機(jī)動(dòng)過程的隨機(jī)偏差對(duì)火星大氣入口狀態(tài)精度的影響。

1 最優(yōu)有限推力離軌制導(dǎo)問題描述

在火星探測(cè)過程中，當(dāng)飛行器從環(huán)繞停泊軌道至進(jìn)入大氣的過程之間，需執(zhí)行連續(xù)推力機(jī)動(dòng)實(shí)現(xiàn)離軌過程，因此，整個(gè)離軌過程從問題本質(zhì)上屬于多重弧段（即包括多段推力弧和滑行弧）有限推力軌道機(jī)動(dòng)問題。本文考慮到實(shí)際任務(wù)時(shí)間限定和有利實(shí)際工程應(yīng)用等客觀因素，以一次“推–滑”兩弧段轉(zhuǎn)移形式為離軌過程飛行器的機(jī)動(dòng)模式，并在此基礎(chǔ)上給出最優(yōu)制導(dǎo)方法。

盡管離軌機(jī)動(dòng)本質(zhì)屬于有限推力機(jī)動(dòng)過程，但部分約束形式等方面又具有其獨(dú)特的一面。由于一般進(jìn)入大氣前的停泊軌道高度較低，因此，為保證制導(dǎo)算法在進(jìn)行離軌軌跡求解時(shí)的實(shí)時(shí)性要求，其動(dòng)力學(xué)方程的形式、狀態(tài)變量以及協(xié)狀態(tài)變量的求解等均會(huì)在一定的精度要求和假設(shè)前提下進(jìn)行簡(jiǎn)化，所以動(dòng)力學(xué)模型以及問題的描述形式也不同于一般的有限推力軌道機(jī)動(dòng)問題。

1.1 動(dòng)力學(xué)模型

通過火星慣性系下的動(dòng)力學(xué)模型來描述離軌機(jī)動(dòng)問題。其對(duì)應(yīng)的無量綱動(dòng)力學(xué)方程為

其中：r為火心慣性系下的無量綱位置矢量；V為火心慣性系下的無量綱速度矢量；m為有量綱的飛行器質(zhì)量；T為飛行器發(fā)動(dòng)機(jī)推力，Isp為飛行器發(fā)動(dòng)機(jī)比沖；g0為火星表面引力加速度；υ 為推力方向在火心慣性系下的單位矢量。在無量綱化過程中，距離的無量綱參數(shù)為火星半徑R；速度的無量綱參數(shù)為；時(shí)間的無量綱參數(shù)為。

由于大氣進(jìn)入前的停泊軌道高度較低，且位置矢徑大小對(duì)應(yīng)的高度介于停泊軌道高度與大氣邊緣高度之間，因此，整個(gè)離軌過程位置矢徑大小變化微小?；谠撎匦?，在每次迭代制導(dǎo)的過程中，位置矢徑可近似為[1]

其中：r0為每次制導(dǎo)迭代初始時(shí)刻的無量綱位置矢徑大小。因此，式（1）中動(dòng)力學(xué)方程中速度矢量的微分項(xiàng)可重新表示為[1]

1.2 離軌機(jī)動(dòng)最優(yōu)性條件

基于式（1）和式（3）的動(dòng)力學(xué)模型，在給出最優(yōu)性條件之前，首先給出協(xié)狀態(tài)方程的表達(dá)式。假設(shè)速度矢量、位置矢量以及質(zhì)量對(duì)應(yīng)的協(xié)狀態(tài)量分別為pr，pV和pm，基于最優(yōu)控制理論[14]，系統(tǒng)的哈密頓函數(shù)為

其中，協(xié)狀態(tài)量pr和pV對(duì)應(yīng)的微分方程為

在離軌過程中，由于“推–滑”離軌機(jī)動(dòng)模式下，飛行器發(fā)動(dòng)機(jī)機(jī)動(dòng)模式固定，所以此處控制變量只有“推”弧段的推力方向單位矢量。根據(jù)極大值原理，最優(yōu)離軌過程對(duì)應(yīng)的發(fā)動(dòng)機(jī)推力方向單位矢量υ 為

上述式（6）中給出的推力方向單位矢量 υ?是最優(yōu)離軌機(jī)動(dòng)過程的必要條件。

1.3 協(xié)態(tài)變量和狀態(tài)變量的解析解

在如式（2）的假設(shè)前提下，離軌機(jī)動(dòng)的動(dòng)力學(xué)方程中部分項(xiàng)被線性化，而式（5）中的協(xié)狀態(tài)微分方程直接具有線性微分方程的基本形式，這種特性給微分方程的解析求解提供了一定的便利。因此，本文采用一種將協(xié)態(tài)變量和狀態(tài)變量變換的方法，求解得到狀態(tài)方程和協(xié)態(tài)方程的解析解[1-2,12]。事實(shí)上，這種解析求解方法目前已經(jīng)成熟應(yīng)用于離軌制導(dǎo)和入軌制導(dǎo)等問題中。下面將給出該方法求解狀態(tài)和協(xié)態(tài)變量解析解的具體過程。

首先，對(duì)協(xié)態(tài)變量進(jìn)行變化，即引入一個(gè)新的量為

對(duì)于新變量，此處省略推導(dǎo)過程，直接給出對(duì)應(yīng)的從式（5）變形得到的微分方程的解析解，即

其中：t和t0分別是當(dāng)前時(shí)間和初始時(shí)間，I3為3階單位矩陣，λ (t0)為初始時(shí)刻對(duì)應(yīng)的值。

其次，對(duì)于狀態(tài)變量，此處同樣需進(jìn)行變換從而使得狀態(tài)變量的微分方程可以解析得到，引入一個(gè)新的量為

根據(jù)參考文獻(xiàn)[12]中的方法，在給出如式（9）的變換之后，從狀態(tài)變量轉(zhuǎn)換得到的新變量對(duì)應(yīng)的微分方程具有如式（10）的解析解形式，即

其中： Θ (t0)為初始時(shí)刻對(duì)應(yīng)的值，為

此外，式（10）中的向量 Γ (t,t0)的表達(dá)式為兩個(gè)子向量的組合，其具體表達(dá)式為

其中，子向量的表達(dá)式分別為

對(duì)于式（13）和式（14）中子向量的值，在每一時(shí)刻均采用數(shù)值積分求解，一般情況下，只要數(shù)據(jù)節(jié)點(diǎn)之間的間隔足夠小，其積分精度能夠得到保證。為了保證求解速度的要求，兩個(gè)時(shí)間節(jié)點(diǎn)間隔之間只需進(jìn)行一次積分即可，數(shù)值積分方法諸如4階龍格庫塔法以及米勒法等均能夠以較高精度實(shí)現(xiàn)對(duì)上述積分函數(shù)的求解。具體求解公式此處不再贅述。

1.4 離軌過程的終端約束及兩點(diǎn)邊值問題

1.3小節(jié)給出了狀態(tài)變量和協(xié)態(tài)變量的解析計(jì)算方法，對(duì)于離軌制導(dǎo)問題，從問題本質(zhì)上屬于典型的兩點(diǎn)邊值問題，因此，下面需要分別確定兩點(diǎn)邊值問題對(duì)應(yīng)的變量及終端約束。對(duì)于離軌機(jī)動(dòng)過程，一般初始狀態(tài)是確定的，初始協(xié)態(tài)變量和總的離軌時(shí)間是可變的，且終端狀態(tài)受約束。在離軌制導(dǎo)過程中，制導(dǎo)參數(shù)包括初始協(xié)態(tài)變量和總離軌時(shí)間，但是由于終端狀態(tài)變量受約束，因此直接根據(jù)變分法得到的制導(dǎo)終端約束會(huì)包含拉格朗日乘子，因此下面將推導(dǎo)給出只包含狀態(tài)變量和協(xié)態(tài)變量的制導(dǎo)終端約束。

對(duì)于火星離軌機(jī)動(dòng)過程，一般包括如下3個(gè)終端狀態(tài)約束，即

其中：rf和Vf分別是終端位置矢量和速度矢量，而rEI、VEI和 γEI分別為基于上述3個(gè)終端狀態(tài)約束。對(duì)于最優(yōu)離軌過程，其橫截條件為

其中：tf是終端時(shí)間；μ1、μ2和μ3是拉格朗日乘子。

為了消去式（18）中的拉格朗日乘子，下面將結(jié)合式（15）～（17）與式（18）給出只包含狀態(tài)和協(xié)狀態(tài)的制導(dǎo)終端約束。首先對(duì)式（18）等式左右分別轉(zhuǎn)置之后右乘位置矢量rf和速度矢量Vf，從而引入式（15）～（17）的等式關(guān)系，也即

通過式（19）中的3個(gè)公式，可以直接求解得到3個(gè)拉格朗日乘子 μ1、μ2和μ3，然后將其帶入式（18），可得到不含拉格朗日乘子的橫截條件，其最終終端約束方程為

其中，中間矩陣A為

式（20）即為最優(yōu)“推–滑”離軌制導(dǎo)對(duì)應(yīng)的由狀態(tài)變量和協(xié)態(tài)變量構(gòu)成的終端約束。此外，由于終端時(shí)間自由，因此，需要額外增加終端哈密頓函數(shù)為0的約束，即

至此，制導(dǎo)過程對(duì)應(yīng)的變量和約束方程已給出。

在火星離軌制導(dǎo)過程中，為保證收斂性和求解速度每個(gè)制導(dǎo)周期內(nèi)問題的求解均采用萊文貝格–馬夸特（Levenberg-Marquardt）方法[15]。一般情況下，制導(dǎo)變量初值選取的優(yōu)劣對(duì)首次制導(dǎo)方程求解的收斂速度有明顯影響，因此，此處基于實(shí)際數(shù)值仿真結(jié)果，給出較為一般性得制導(dǎo)變量初值，即協(xié)態(tài)變量初值此外，離軌總時(shí)間t的初f值選取對(duì)收斂速度影響不大，一般選取為初始軌道周期的一般較為合理。下文將給出離軌制導(dǎo)迭代流程。

2 仿真結(jié)果與分析

為驗(yàn)證上述最優(yōu)“推–滑”離軌制導(dǎo)方法的有效性和魯棒性，本節(jié)以火星探測(cè)任務(wù)為背景，給出相關(guān)的仿真結(jié)果。首先，給定火星探測(cè)器在進(jìn)入大氣前所在的停泊軌道根數(shù)，其數(shù)值通過轉(zhuǎn)換即可得到離軌機(jī)動(dòng)的初始狀態(tài)變量，如表1所示。

表1 離軌制導(dǎo)的初始軌道根數(shù)Table 1 Initial orbit elements of deorbit guidance

此外，給定火星半徑R= 3 396 km，因此，初始軌道實(shí)際上為高度為350 km的圓軌道。在離軌機(jī)動(dòng)過程中，飛行器的發(fā)動(dòng)機(jī)推力T= 1 000 N，發(fā)動(dòng)機(jī)比沖Isp=311 s，飛行器初始質(zhì)量假設(shè)為4 000 kg。由于整個(gè)離軌過程以飛行器到達(dá)大氣進(jìn)入邊緣為終止位置，因此，此處給定大氣邊緣高度為128 km。此外，所有仿真均假定目標(biāo)大氣入口處速度約束為VEI= 3 516.2 m/s，飛行路徑角γEI= –2.358°，而進(jìn)入點(diǎn)的位置不進(jìn)行約束。

2.1 火星離軌機(jī)動(dòng)最優(yōu)解

基于上述離軌過程狀態(tài)變量初值和發(fā)動(dòng)機(jī)參數(shù)，本節(jié)給出了不考慮隨機(jī)偏差影響的火星“推–滑”離軌機(jī)動(dòng)最優(yōu)解。圖1給出了最優(yōu)“推–滑”離軌機(jī)動(dòng)對(duì)應(yīng)的轉(zhuǎn)移軌跡。從圖中可以看出，整個(gè)過程機(jī)動(dòng)弧段較短，最終實(shí)現(xiàn)的是類似于一個(gè)脈沖機(jī)動(dòng)的效果，而滑行弧段較長(zhǎng)。此外，從轉(zhuǎn)移軌跡也能看出，在終端速度和飛行路徑角約束下，整個(gè)機(jī)動(dòng)幾乎在同一軌道面內(nèi)。

圖1 最優(yōu)離軌過程軌跡Fig. 1 Optimal deorbit trajectory

通過仿真可知，整個(gè)“推–滑”離軌機(jī)動(dòng)最優(yōu)解對(duì)應(yīng)的離軌總時(shí)間為2 241.46 s，其中機(jī)動(dòng)弧段總時(shí)間為295.31 s，而滑行弧段歷時(shí)1 946.15 s。因此通過發(fā)動(dòng)機(jī)參數(shù)推算，整個(gè)離軌機(jī)動(dòng)所消耗的燃料為255.696 kg，對(duì)應(yīng)的速度增量為76.29 m/s。對(duì)于機(jī)動(dòng)弧段，圖2給出了該階段對(duì)應(yīng)的發(fā)動(dòng)機(jī)推力單位方向矢量在慣性系下的分布。為了更為直觀展示推力方向與速度的夾角分布，圖2還給出了推力方向角，此處推力方向角定義為推力方向與飛行器速度方向的夾角，從結(jié)果可以看出，整個(gè)離軌過程，推力幾乎近似沿速度的反方向，也即減速制動(dòng)從而使得飛行器降軌進(jìn)入大氣。此外，通過仿真計(jì)算，表2給出了最優(yōu)解對(duì)應(yīng)的初始協(xié)態(tài)變量的值。

圖2 推力單位方向矢量分量與推力方向角變化曲線Fig. 2 Components of Thrust Unit Vector and Thrust Angle

表2 最優(yōu)離軌過程對(duì)應(yīng)的協(xié)態(tài)變量初值Table 2 Initial costate value of optimal deorbit process

2.2 蒙特卡羅仿真分析

為了驗(yàn)證上述最優(yōu)“推–滑”離軌制導(dǎo)方法的魯棒性和精度，本節(jié)將在不確定性隨機(jī)偏差存在的基礎(chǔ)上進(jìn)行蒙特卡羅仿真，并分析不確定性因素的存在對(duì)“推–滑”離軌制導(dǎo)過程的影響。在制導(dǎo)過程中，給定機(jī)動(dòng)弧段制導(dǎo)周期為5 s，在每一次制導(dǎo)起始增加隨機(jī)偏差，根據(jù)工程經(jīng)驗(yàn)，此處給定位置偏差的3σ值為10 m，速度偏差的3σ值為0.1 m/s，蒙特卡羅打靶次數(shù)為1 000。

圖3 最優(yōu)離軌過程軌跡蒙特卡羅仿真分布Fig. 3 Monte Carlo simulation of optimal deorbit trajectories

圖4 離軌終端狀態(tài)及總離軌時(shí)間蒙特卡羅仿真散布Fig. 4 Dispersions of deorbit terminal state and total deorbit time

圖3給出了蒙特卡羅仿真得到的最優(yōu)離軌過程軌跡分布，從軌跡圖可以看出盡管隨機(jī)偏差存在，但整個(gè)軌跡散布范圍很小。圖4分別給出了終端狀態(tài)和離軌總時(shí)間在隨機(jī)偏差存在下的散布情況。從圖中可以看出，進(jìn)入點(diǎn)路徑角的散布與標(biāo)稱值 γEI相差最大約為0.046°，而進(jìn)入點(diǎn)速度的散布與標(biāo)稱值VEI相差最大約為0.5 m/s，因此，整個(gè)制導(dǎo)結(jié)果精度滿足進(jìn)入點(diǎn)狀態(tài)偏差的精度的要求。從離軌總時(shí)間散布來看，其受隨機(jī)偏差的影響較為明顯，但其與最優(yōu)解的值2 241.46 s相比最大偏差僅在8 s左右。此外，為反映機(jī)動(dòng)過程推力方向的分布，圖5給出了推力單位方向矢量分量與推力方向角在隨機(jī)偏差影響下的散布情況，可以看出盡管隨機(jī)偏差存在，但推力方向仍處于標(biāo)稱推力方向附近。

圖5 推力單位方向矢量分量與推力方向角的蒙特卡羅仿真分布Fig. 5 Distributions of Thrust Unit Vector components and Thrust Angle

3 結(jié)束語

針對(duì)火星探測(cè)中的離軌機(jī)動(dòng)問題，本文給出了一種基于解析狀態(tài)變量和協(xié)態(tài)變量的最優(yōu)“推–滑”離軌制導(dǎo)方法。該方法在位置矢徑局部線性化的假設(shè)下分別給出了最優(yōu)“推–滑”離軌機(jī)動(dòng)的必要條件和只由狀態(tài)變量和協(xié)態(tài)變量構(gòu)成的制導(dǎo)終端約束。通過仿真分析發(fā)現(xiàn)，該方法所得到的離軌過程近似為機(jī)動(dòng)，且機(jī)動(dòng)弧段的推力方向幾乎沿速度反方向。蒙特卡羅仿真表明該方法在滿足制導(dǎo)精度的前提下具有較強(qiáng)的魯棒性。上述研究結(jié)果可為我國(guó)未來火星探測(cè)離軌機(jī)動(dòng)方案設(shè)計(jì)提供技術(shù)支持和理論依據(jù)。