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

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

0 引 言

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

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

1 軟著陸動力學模型

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

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

(1)

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

(2)

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

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

(3)

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

(4)

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

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

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

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

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

(5)

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

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

(6)

尋求最優(yōu)控制變量u*=[θ*,φ*]T，使得探測器在短時間由式(3)的初始狀態(tài)完成動力下降過程，并滿足式(4)的終端約束。與之對應的狀態(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 改進的制導律設計

結合月球探測器軟著陸的動力下降段特點：經(jīng)過霍夫曼轉移，探測器在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)的探測器動力學方程可以寫作：

(13)

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

(14)

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

(15)

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

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

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

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

本文將采用基于標稱軌跡的離散化方法解決分段最優(yōu)控制問題，月球軟著陸最優(yōu)制導問題轉化為含有一般約束的數(shù)學規(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]。為此，引入時間尺度變換法，假設新的時間變量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)

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

(25)

由此，轉換狀態(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，由改進的多項式制導方法得到的著陸時間為607.956 s，利用CVP方法得到的最終著陸時間為605.85 s,減少了2.1 s。

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

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

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

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

5 著陸精度仿真分析

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

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

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

圖8給出了700次打靶計算中動力下降段結束時高度誤差的散布情況。圖8(a)所示是起始點誤差在±500 m情況下動力下降段結束時高度誤差和著陸速度的散布情況，可以看出在絕大多數(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情況下動力下降段結束時軌跡的散布情況，可以看到動力下降段結束時，著陸軌跡大致分布在同一平面上，證明算法滿足任務要求的著陸高度精度要求。

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

6 結 論

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