月球探測器動力下降段最優(yōu)軌跡參數(shù)化方法

喬衍迪,張澤旭,2,鄧涵之,徐田來,2

(1.哈爾濱工業(yè)大學(xué)深空探測基礎(chǔ)研究中心，哈爾濱 150001；2. 陜西省組合與智能導(dǎo)航重點實驗室，西安 710000；3. 濟南大學(xué)自動化與電氣工程學(xué)院，濟南 250022)

0 引 言

近年來，包括中美在內(nèi)的多個國家組織均推出了自己的探月計劃，月球探測重新成為深空探測任務(wù)的焦點。月球探測器軟著陸過程主要分為動力下降段、姿態(tài)調(diào)整段和垂直下降段三個階段，其中動力下降段要求探測器距離從距月面15 km下降到3 km，探測器的速度從1.7 km/s降至0附近，為保證后續(xù)姿態(tài)調(diào)整段具有充分的高度與速度余量，要求動力下降段的末制導(dǎo)具有高精度的速度和高度。但是，從月球軌道離軌進行霍曼轉(zhuǎn)移，進入動力下降段容易出現(xiàn)高度控制偏差，本文充分考慮動力下降段的起始高度波動影響，提出一種高精度參數(shù)化制導(dǎo)方法，在燃料更省的基礎(chǔ)上實現(xiàn)高精度動力下降段末制導(dǎo)。

隨著航天技術(shù)的發(fā)展，對月球探測器實現(xiàn)自主高精度制導(dǎo)的需求更加迫切。月球軟著陸制導(dǎo)方法主要分為三類：重力轉(zhuǎn)彎制導(dǎo)法、標稱軌道制導(dǎo)法和顯式制導(dǎo)法，國內(nèi)外對各種制導(dǎo)方法已經(jīng)開展了多年深入的研究工作。重力轉(zhuǎn)彎制導(dǎo)法原理簡單，多出現(xiàn)在20世紀低成本探月任務(wù)設(shè)計中。標稱軌跡法即設(shè)計一條著陸軌跡，探測器根據(jù)位置和速度信息跟蹤理想軌跡實現(xiàn)軟著陸。文獻[1-2]將月球軟著陸軌跡進行離散化，再利用序列二次規(guī)劃方法(Sequential quadratic programming,SQP)進行數(shù)值求解。文獻[3]針對軟著陸與采樣返回軌道進行B樣條數(shù)值逼近，在復(fù)雜約束下得到了高精度著陸軌跡，但是標稱軌跡法對于初始偏差敏感，不一定能夠?qū)崿F(xiàn)完美跟蹤。顯式制導(dǎo)法根據(jù)探測器的實時位置和速度信息，實時計算制導(dǎo)指令。文獻[4-5]利用高斯-偽譜法實現(xiàn)了探測器二維定點軟著陸軌道優(yōu)化，求解了月面安全的著陸區(qū)范圍。文獻[6-8]對閉環(huán)軟著陸軌跡跟蹤算法進行了研究，有效地減小了制導(dǎo)誤差，提高了制導(dǎo)系統(tǒng)的穩(wěn)定性。文獻[9]考慮到月球軟著陸任務(wù)的工程實際，設(shè)計了可以實現(xiàn)機載計算機自主規(guī)劃的月球軟著陸顯式制導(dǎo)律，該方法計算簡單，可以快速規(guī)劃出著陸軌跡。文獻[10]提出了混合制導(dǎo)方法(Hybrid guidance law,HGL)，但是該制導(dǎo)律只適用于軟著陸姿態(tài)調(diào)整段，以實現(xiàn)目標著陸點的重決策；近年來隨著最優(yōu)問題的研究深入，更多的數(shù)學(xué)規(guī)劃方法被應(yīng)用到時深空探測領(lǐng)域。文獻[11-12]運用凸優(yōu)化理論解決月球軟著陸問題常出現(xiàn)的局部無解問題，得到的最終著陸高度誤差約200 m。文獻[13]將狀態(tài)和控制變量進行參數(shù)化，結(jié)合迭代的方法，解決了非線性系統(tǒng)最優(yōu)控制問題。文獻[14-16]結(jié)合了控制變量參數(shù)化方法和強化技術(shù)，并成功應(yīng)用于月球軟著陸軌跡的優(yōu)化問題，但構(gòu)建的動力學(xué)參數(shù)復(fù)雜，控制變量較多，不適合運用探測器的機載計算機解算，且需要用到專業(yè)的數(shù)學(xué)軟件，實用性較差。本文綜合考慮月球探測器軟著陸動力下降段的快速性和工程實用性，簡化控制變量，提出了一種改進的燃料最優(yōu)的顯式制導(dǎo)律，并利用控制變量參數(shù)化方法，快速規(guī)劃著陸軌跡，引入時間尺度變換方法，求解得到了軟著陸動力下降段最優(yōu)著陸軌跡。

1 軟著陸動力學(xué)模型

月球探測器軟著陸動力下降段中，探測器從15 km的橢圓軌道下降至月球表面3 km處，且橫向行程較長，但是歷時較短，因此可以將月球視為均勻球體，忽略月球自轉(zhuǎn)，建立三維軟著陸模型[17]。

用U，V，W表示探測器在探測器著陸坐標系中的速度分量，有：

(1)

式中:α和β為著陸坐標系的方位角和仰角,R為探測器與月心的距離。在不考慮攝動影響且忽略月球自轉(zhuǎn)情況下，得到月球探測器動力下降段的動力學(xué)模型，有：

(2)

式中:μ為月心引力常數(shù),F為探測器制動發(fā)動機推力且為常數(shù),Φ為探測器推力仰角,Θ為推力的推力傾角,Isp為制動發(fā)動機的比沖,gE為地球重力加速度,m為探測器的總質(zhì)量。

探測器軟著陸過程動力下降段的初始位置和初始速度可寫為：

(3)

動力下降階段的終端約束條件除了速度分量約束外，還包括高度約束：

(4)

2 月球探測器軟著陸燃料最優(yōu)軌道求解

2.1 最優(yōu)控制問題

利用最優(yōu)控制方法直接求解月球探測器軟著陸問題比較復(fù)雜，動力學(xué)模型具有較強的非線性，且難以獲得解析形式的表達式。通常需要利用給定初值進行迭代的方法求解，不利于在飛行器上實現(xiàn)自主控制。因此，需要對軟著陸模型進行簡化。

引入兩個假設(shè)：1)假設(shè)月球表面為平面，引力場均勻，且月球引力加速度為常數(shù)，則月心慣性坐標系和探測器著陸坐標系Oxoyozo的瞬時坐標軸一致，用u，v，w表示探測器的速度分量；2)假設(shè)探測器采用變推力發(fā)動機，即制動加速度不隨著探測器的質(zhì)量發(fā)生變化，是常值。

基于假設(shè)，式(1)所表示的動力學(xué)模型可簡化為：

(5)

式中:aN表示探測器推力F與質(zhì)量m的比值,gm為月球的引力加速度,θ為簡化后探測器推力仰角,φ為推力傾角。

探測器在動力下降段需要消耗大量的燃料抵消初始速度，因此應(yīng)將燃料作為最優(yōu)控制問題的性能指標。對于動力下降段的連續(xù)制動控制，燃料最優(yōu)目標函數(shù)：

(6)

尋求最優(yōu)控制變量u*=[θ*,φ*]T，使得探測器在短時間由式(3)的初始狀態(tài)完成動力下降過程，并滿足式(4)的終端約束。與之對應(yīng)的狀態(tài)變量[x*,y*,z*]T即為動力下降段的最優(yōu)著陸軌跡。

對于月球探測器軟著陸動力下降段的最優(yōu)控制問題，通常根據(jù)狀態(tài)變量和控制變量的初值，猜測協(xié)狀態(tài)變量或中間變量，利用迭代的方法對當前狀態(tài)進行優(yōu)化，即將剩余時間tgo劃分出多個時間節(jié)點，分段進行控制變量u*=[θ*,φ*]T的優(yōu)化，狀態(tài)變量不斷更迭，直到滿足終值約束。

2.2 改進的制導(dǎo)律設(shè)計

結(jié)合月球探測器軟著陸的動力下降段特點：經(jīng)過霍夫曼轉(zhuǎn)移，探測器在xoOozo平面內(nèi)水平初始速度遠大于另外兩個方向的初始速度。因此，制動發(fā)動機的主要工作是抵消水平初始速度u0，則對于控制變量u=[θ,φ]T，推力傾角φ≈180°，并將推力仰角θ視為小量。

利用龐得里亞金極大值原理，可以得到Hamiltonian方程：

H=1+λuaNcosθcosφ+λvaNcosθsinφ+

λw(aNsinθ-gm)+λxu+λyv+λzw+

(7)

式中:λi(i=u,v,w,x,y,z)為共軛變量，且滿足：

(8)

則最優(yōu)解為:

(9)

(10)

簡化為：

(11)

tanθ≈sinθ=κ1+κ2t

(12)

式(5)的探測器動力學(xué)方程可以寫作：

(13)

求解兩點邊值問題，得到ψ0、κ1、κ2和tgo。

(14)

則月球探測器軟著陸的動力下降段的歷時tgo為：

(15)

改進的顯式制導(dǎo)律在探測器軟著陸過程中無需解算協(xié)狀態(tài)變量，根據(jù)當前的速度解算得到推力仰角，探測器高度的表達式為三次多項式，通過多步迭代可以獲得探測器的最優(yōu)燃料制導(dǎo)律的著陸軌跡。

3 控制變量參數(shù)化方法求解著陸軌跡

3.1 控制變量參數(shù)化

上節(jié)中，根據(jù)龐得里亞金極大值原理將月球軟著陸轉(zhuǎn)化為兩點邊值問題求解，但由于協(xié)狀態(tài)變量的初值沒有物理意義，迭代計算引入了積分問題，計算量較大。因此，本節(jié)采用控制變量參數(shù)化(CVP)方法，將月球軟著陸軌道離散化，分段求解最優(yōu)控制問題。

本文將采用基于標稱軌跡的離散化方法解決分段最優(yōu)控制問題，月球軟著陸最優(yōu)制導(dǎo)問題轉(zhuǎn)化為含有一般約束的數(shù)學(xué)規(guī)劃問題，通過內(nèi)點法求解[18]。

(16)

對于式(16)，尋找控制變量u，使目標函數(shù):

(17)

最小,Φ0(t,x(t),u(t))其中為終端性能指標,L0(t,x(t),u(t))為積分性能指標。

(18)

3.2 時間尺度變換

若時間序列不固定那么求取參數(shù)的梯度信息變得很困難[20]。為此，引入時間尺度變換法，假設(shè)新的時間變量s∈[0,1]，將不確定的下降時間t1重定義到s∈[0,1]中：

(19)

(20)

(21)

t1(1)=tf
t1(0)=0

定義:

wp(s)=up(t1(s))

則

(22)

將時間尺度變換式(16)代入非線性系統(tǒng)中：

δp,σp)υp(s)

(23)

及

(24)

由于軟著陸任務(wù)動力下降段的主要任務(wù)是抵消探測器較大的水平速度，使其在姿態(tài)調(diào)整段的初始速度接近于0，且下降高度不低于3 km，將速度u，w與月心距r作為路徑約束。最優(yōu)控制問題的目標函數(shù)寫作：

(25)

由此，轉(zhuǎn)換狀態(tài)約束進入目標函數(shù)中，采用變分法給出系統(tǒng)(16)梯度信息，定義變分系統(tǒng)如下：

(26)

則一階梯度信息可以通過下式得到：

(27)

由此采用變分法，可以求解狀態(tài)方程組，具有較好的數(shù)值穩(wěn)定性和較高的精度，所需求解的微分方程不多，不會給求解帶來負擔。先給出梯度信息與終端約束之差ε的值，檢驗是否符合約束條件，若不符合，重新規(guī)劃控制變量參數(shù)的取值范圍。顯然，對于每條軌跡一定能夠找到足夠小的取值范圍使約束條件滿足，當減小到給定值時這個計算過程就停止。很明顯，當誤差足夠小時，能夠達到滿足要求的近似。誤差定為ε≤10-3。

下面給出基于標稱軌跡的控制變量參數(shù)化算法具體的實現(xiàn)流程圖，如圖1所示。

圖1 基于標稱軌跡的控制變量參數(shù)化算法流程圖Fig.1 Flow chart of trajectory generation algorithm based on control variable parameterization method

4 仿真校驗

月球探測器的初始質(zhì)量m0=750 kg；制動發(fā)動機的最大推力為Fmax=1750 N，比沖為Isp=316 s；初始水平速度為U0=1692 m/s，其它兩個方向初始速度分量為W0=V0=0。月球引力常數(shù)μ=4.903×1012m3/s2，近月點與月心距離R0=1753 km，最終著陸點與月球表面的距離為RF=1741.03 km，最終著陸速度約束分別為|Uf|≤1 m/s,|Wf|≤1.5 m/s,|Vf|≤3.8 m/s，由改進的多項式制導(dǎo)方法得到的著陸時間為607.956 s，利用CVP方法得到的最終著陸時間為605.85 s,減少了2.1 s。

如圖2所示，動力下降段結(jié)束時，由改進的多項式制導(dǎo)律得到的探測器的最終質(zhì)量為396.94 kg，由CVP算法得到的最終質(zhì)量為398.69 kg，燃料少消耗了1.75 kg。

圖2 質(zhì)量變化曲線Fig.2 The mass curve

如圖3所示，月球探測器在動力下降段結(jié)束時三軸速度分別為Uf=0.309 m/s,Wf=-1.094 m/s,Vf=2.981 m/s，滿足著陸速度約束，且豎直速度分量較小，經(jīng)過姿態(tài)調(diào)整段，可以實現(xiàn)豎直安全著陸。

圖3 速度變化曲線Fig.3 The velocity curve of soft landing

如圖4和圖5所示，月球軟著陸的最終著陸高度為1741.03 km，水平行程約為638.7 km。

圖4 高度變化曲線Fig.4 The height curve

圖5 著陸軌跡曲線Fig.5 The trajectory of soft landing

圖6所示為月球探測器的推力仰角的變化曲線，可以看出采用改進的多項式制導(dǎo)律在動力下降段末端需要進行固定參數(shù)，否則會產(chǎn)生跳變，而CVP方法得到的控制曲線更加平緩，在450 s后基本不再需要再規(guī)劃，更具有實現(xiàn)性。

圖6 控制變量θ變化曲線Fig.6 The control variable θ curve

5 著陸精度仿真分析

下面給出測量誤差情況下采用變推力方案的著陸參數(shù)和著陸誤差分布情況。以下關(guān)于誤差的分析均采用蒙特卡羅打靶,打靶次數(shù)為700次，假設(shè)各誤差均符合正態(tài)分布。

圖7給出了月球探測器在月心慣性坐標系3個方向上初始著陸點在1753±1 km的散布示意圖。本文中給出的動力下降段起始點誤差散布分別在500 m、1 km和3 km范圍的球體內(nèi)。

圖7 著陸起始點分布示意圖Fig.7 Initial landing site dispersion

圖8給出了700次打靶計算中動力下降段結(jié)束時高度誤差的散布情況。圖8(a)所示是起始點誤差在±500 m情況下動力下降段結(jié)束時高度誤差和著陸速度的散布情況，可以看出在絕大多數(shù)情況下，高度誤差不大于30 m，符合著陸安全要求，且90%以上的豎直速度小于1 m/s，水平速度不大于1.5 m/s；圖8(b)所示為起始點誤差±1000 m時的誤差分布情況，可以看出軌跡的高度誤差依然保持在50 m以內(nèi)，大致滿足σ分布準則，且豎直方向的速度在±1 m/s以內(nèi)，水平方向的速度絕大多數(shù)在1～2 m/s之間，符合仿真給了的精度要求；圖8(c)所示為高度誤差±3000 m時的誤差分布情況，多數(shù)誤差分布在±100 m以內(nèi)，滿足σ分布準則，半數(shù)以上分布在±20 m內(nèi)，豎直方向的著陸速度不大于1 m/s，水平方向的著陸速度不大于4 m/s，為姿態(tài)調(diào)整段安全著陸提供保障。

圖8 高度誤差分布情況Fig.8 Dispersion in final altitude

圖9給出了700次打靶計算中著陸軌跡的散布情況。圖9(a)、(b)和(c)分別給出了起始點誤差±500 m、±1000 m和±3000 m情況下動力下降段結(jié)束時軌跡的散布情況，可以看到動力下降段結(jié)束時，著陸軌跡大致分布在同一平面上，證明算法滿足任務(wù)要求的著陸高度精度要求。

圖9 軌跡分布圖Fig.9 3D landing trajectories MC simulation

6 結(jié) 論

本文提出了一種改進的顯式制導(dǎo)方法，簡化控制變量，避免了制導(dǎo)末段剩余時間發(fā)散問題，并利用控制變量參數(shù)化方法將著陸軌跡分段，引入時間尺度變換方法，運用經(jīng)典的非線性規(guī)劃方法求解每段控制變量最優(yōu)數(shù)值解，最終得到月球探測器動力下降段燃料最優(yōu)著陸軌跡。仿真結(jié)果表明，相比于傳統(tǒng)的顯式制導(dǎo)律，本文提出的方法燃料消耗更少。蒙特卡羅打靶實驗證明，本文提出的算法在初始著陸點存在±20%誤差的情況下，依然能夠保證足夠的著陸高度和速度裕量，保障姿態(tài)調(diào)整段和垂直著陸段探測器安全著陸。