初彥峰，穆榮軍，梁浩，崔乃剛

（1.哈爾濱工業(yè)大學航天學院，哈爾濱 150001；2.北京宇航系統(tǒng)工程研究所，北京 100076）

0 引言

作為太陽系起源的原始物質，小行星蘊含了行星形成及動力學過程的信息，其為研究太陽系演化、深空資源開發(fā)、行星防御等空間科學問題提供了重要途徑［1］。小行星探測已成為當代太空爭奪的新前沿，也是一個國家綜合國力的體現(xiàn)。最初人類只能利用天文望遠鏡對小行星進行觀測，而隨著科學技術的發(fā)展，自20 世紀80 年代以來，人類已經(jīng)實現(xiàn)了對小行星的環(huán)繞飛行及采樣返回等探測任務［2］。

NASA 發(fā)射的伽利略木星探測器于1991 年對加斯普拉小行星進行了環(huán)繞飛行，這是人類第一次實現(xiàn)對小行星近距離探測，并于1993年飛越了編號為243的Ida小行星［3］；2018年黎明號探測器完成對小行星灶神星（Vesta）和谷神星（Ceres）的探測［4］；2020 年奧西里斯-雷克斯探測器實現(xiàn)了在Bennu 小行星上著陸及采樣，之后其攜帶約250 g 樣本于2023 年在美國猶他州安全著陸，這是NASA 首個小行星采樣返回任務［5］。日本宇宙航空研究開發(fā)機構（JAXA）于2014年發(fā)射了國際上首顆采樣返回式小行星探測器——隼鳥2號，并于2020 年實現(xiàn)了對小行星龍宮（1999JU3）采集樣本及返回的探測任務［6］。ESA 和NASA 聯(lián)合推出了小行星撞擊與偏轉評估計劃，于2022 年利用DART 航天器撞擊Didymos 小行星，以迫使其偏離運行軌道［7］。我國嫦娥二號月球探測器于2012年對Toutatis小行星進行了交會探測，拍攝了多張高清的小行星光學圖像［8］。

小行星引力場建模是探測器實現(xiàn)精準著陸的基礎。由于形狀不規(guī)則，小行星附近的引力場呈現(xiàn)出不規(guī)則性。球諧系數(shù)引力場模型使用球形勢的諧波級數(shù)展開，所以該方法對于非球形天體引力場的描述有局限性；內(nèi)球諧引力場模型的收斂域可以擴展到小天體表面任意一點，但該方法仍無法使用統(tǒng)一表達式實現(xiàn)引力場的全局表征［9］。質點群法采用有限個質點填充小天體內(nèi)部空間的方式進行建模，當觀察點靠近小天體表面時引力場的精度會迅速退化，因此該方法并不適用于小天體表面的動力學研究［10］。多面體引力場模型可以精準地描述任意均質多面體的引力勢、引力加速度和引力梯度張量，被譽為小天體引力場的優(yōu)美解［10-11］。多面體引力場被視作高精準的模型，NASA 已將該模型用于433Eros小行星探測任務的規(guī)劃中［12］。

著陸軌跡對探測器能否安全著陸至關重要。當探測器接近小行星表面時將處于潛在的碰撞危險中，并且探測器將著陸到地形復雜區(qū)域以采集到高價值的樣本。這都要求探測器的著陸軌跡滿足一定的約束條件，同時還需考慮如燃料消耗和著陸時間等性能指標。小行星著陸軌跡優(yōu)化方法已經(jīng)成為一個熱門的研究問題。文獻［13］采用內(nèi)球諧引力場模型對小行星引力場進行建模，通過變量松弛和離散化處理將非凸的3-DOF 著陸問題轉化為更容易求解的等效凸優(yōu)化問題。文獻［14］針對分段狀態(tài)約束的軌跡優(yōu)化問題，采用歸一化方法，并引入了時間膨脹系數(shù)將其轉化成為固定時間的軌跡優(yōu)化問題。文獻［15］提出了基于凸優(yōu)化的軌跡優(yōu)化方法，通過建立最小著陸誤差目標函數(shù)，實現(xiàn)了不同飛行時間下太陽輻射壓力-控制的軌跡優(yōu)化，然后采用信賴域約束得到一系列SOCP 問題。文獻［16］在Lambert 求解器的基礎上提出了一種變飛行時間軌跡優(yōu)化算法，用于解決小行星著陸過程中的燃料消耗問題。

上述方法都是基于3-DOF 的制導算法，即只涉及探測器質心的平移運動，而不考慮探測器姿態(tài)變化。3-DOF動力著陸制導起源于二十世紀的阿波羅登月計劃［17］，其實現(xiàn)相對容易。對于平動和轉動緊密耦合的未來小行星探測任務，3-DOF 制導方案可能無法滿足探測器的任務需求，實現(xiàn)位置和姿態(tài)同步控制的6-DOF制導算法更具優(yōu)勢。

軌跡優(yōu)化方法包括直接法和間接法。間接法采用Pontryagin 原理求解，并高度依賴于初始估計的準確性，所以使用間接法求解比較困難［18］。直接法將軌跡優(yōu)化問題轉化為參數(shù)優(yōu)化問題，進而使用非線性規(guī)劃對參數(shù)優(yōu)化問題進行求解。直接法主要包括偽譜法、凸優(yōu)化方法和打靶法，其中凸優(yōu)化方法具有計算速度快和實時性高的特點，被廣泛地應用于飛行器軌跡優(yōu)化問題［19］。

本文提出一種6-DOF 小行星著陸軌跡序列凸優(yōu)化方法。首先使用多面體引力場模型建立不規(guī)則小行星的引力場；通過對6-DOF 動力學模型及其約束進行線性化和離散處理，將原著陸軌跡優(yōu)化問題轉化為二階錐規(guī)劃問題，然后通過虛擬控制和信賴域增強算法的魯棒性；最后通過模擬著陸433Eros小行星驗證了算法的可行性和有效性。

1 小行星著陸動力學模型

1.1 多面體引力場模型

多面體引力場模型適合于描述形狀不規(guī)則小天體的引力場，通過一系列三維實體來逼近小行星形狀。Werner等［20］利用Gauss 散度定理推導出了多面體的引力場模型，由面引力和邊引力兩部分組成。不規(guī)則小行星附近空間任意一點的引力加速度可表示為：

式中：?U(r)是引力加速度；G為引力常數(shù)；ρ為小行星密度；nf是多面體的面數(shù)；ne是多面體的邊數(shù)；rf和re分別是探測器到第f個面（三角形）和第e條邊的向量；Ff和Ee是和多面體形狀相關的并矢對稱矩陣；ωf和Le是描述計算點與多面體面和邊相對關系的無量綱因子［10，16］。

1.2 小行星著陸動力學模型

本文涉及2 種與探測器運動相關的坐標系，如圖1所示。

圖1 坐標系的定義Fig.1 Definition of coordinate systems

小行星固聯(lián)坐標系OI-XIYIZI（簡寫為OI系）：坐標原點位于小行星的質心，OIZI軸為小行星自轉軸，OIXI軸為小行星最大慣性主軸，OIYI與OIXI，OIZI構成右手坐標系。

探測器本體坐標系OB-XBYBZB（簡寫為OB系）：坐標原點位于探測器的質心，OBXB，OBYB和OBZB分別為探測器的慣性主軸。

探測器的質量消耗計算如下：

式中：m(t)是探測器質量；TB(t)是探測器在OB系下的推力；α=-1 ()Ispg為質量流量系數(shù)，Isp和g分別為推力器比沖和地球標準重力常數(shù)。

探測器的動力學方程可以表示為：

式中：ω=[ω1ω2ω3]T。

四元數(shù)和歐拉角（偏航角ψ、俯仰角θ和滾轉角γ）的轉換關系定義如下：

1.3 狀態(tài)約束和控制約束

首先探測器的質量m(t)應始終大于它的干重mdry，質量約束如下：

輸出推力T(t)應介于最大推力Tmax和最小推力Tmin之間，推力約束如下：

輸出力矩M(t)小于最大力矩Mmax即可，約束如下：

考慮到小行星形狀不規(guī)則，并且存在著自轉，探測器下降時可能會發(fā)生硬著陸現(xiàn)象。為了防止探測器撞擊到小行星，需要對探測器的飛行路徑進行約束。從距離小行星較遠的位置到著陸點上方，探測器主要在小行星的外部空間飛行，因此形成包含小行星的橢球，也就是說橢球體積應大于小行星的實際尺寸。如圖2 所示，最外層黃色的網(wǎng)格代表橢球，最里面灰色的是小行星（本文繪制的小行星是433Eros，其半長軸為16 km × 6.5 km × 6.5 km）。橢球以外的點滿足≥1，整理后得到的橢球約束數(shù)學表達式如下：

圖2 橢球約束Fig.2 The ellipsoid constraint

式中：a，b，c是橢球的半長軸，由其構成的橢球要想將433Eros 小行星包圍，需要滿足a＞16 km，b＞6.5 km，c＞6.5 km；rI=[xyz]T是探測器位置；R是橢球半長軸形成的對角矩陣；tc為動力下降段到最終著陸段的時間。

當探測器進入到最終著陸段時，探測器與著陸點的距離足夠近，則橢球約束不再適用。要實現(xiàn)高精度著陸，在最終著陸段需要確保探測器配備的相機能始終對著陸點進行監(jiān)測。相機視野范圍應該包括著陸點，也就是說相機視野范圍大于相機光軸和視線向量之間的夾角，所形成的視場角約束如圖3所示，定義如下：

圖3 視場角約束Fig.3 The field-of-view constraint

式中：β是相機視場角；ρB是相機在探測器上的安裝位置；dB是相機光軸向量；tf為最終著陸時間；rI和rf分別為OI系下探測器當前位置和最終著陸位置。

2 小行星著陸序列凸優(yōu)化算法

2.1 狀態(tài)方程的凸化

系統(tǒng)狀態(tài)方程式（3）及式（10）、式（11）等都是非凸的，可對其進行線性化，進而獲得凸化的表達式。定義探測器的狀態(tài)向量為x=[rTvTωTqT]T∈R13，控制變量為u=[TTMT]T∈R6。系統(tǒng)狀態(tài)方程可表示為：

式（12）的線性化形式為：

式中相關的系數(shù)矩陣如下：

將探測器的狀態(tài)向量分為平動狀態(tài)xp=[r v]T和轉動狀態(tài)xr=[ω q]T，則式（13）可以分解為2 部分，如下：

其中相關的系數(shù)矩陣如式（16）～（17）所示

式中：E3代表三階單位對角矩陣。

式中：Cω=[03×6，E3，03×4]，Cq=[04×9，E4]。

根據(jù)式（15）～（17）可得式（13）中相關矩陣的具體形式如下：

2.2 約束的凸化

引入松弛因子Γ=[ΓxΓyΓz]T代替控制力T=[ΤxΤyΤz]T，兩者關系可以定義為Tj，j=x，y，z（代表估計值），則有：

對橢球約束（10）進行一階泰勒展開，則有：

同理，對視場角約束（11）進行一階泰勒展開，則有：

式中：h(x)計算如下：

2.3 離散化

將飛行時間[0，tf]離散化為N個子區(qū)間，設每個區(qū)間長度為Δt，則有：

系統(tǒng)的狀態(tài)方程離散化后可以表示為：

式中：?k∈K，K={0，1，2，…，N-1}，i代表迭代次數(shù)，相關系數(shù)矩陣定義如下：

當?k∈Kc，Kc={0，1，2，…，Nc}，Nc=tcΔt，橢球約束離散化后：

當?k∈Kf，Kf={Nc+1，…，N}，視場角約束離散化后：

控制力約束離散化后可以表示為：

2.4 虛擬控制和信賴域

由于對系統(tǒng)方程和約束進行了線性化，這樣的做法將存在人為不可行性，對于飛行時間較短的軌跡優(yōu)化問題可能存在不可行的解［21］。通過在線性化方程中引入虛擬控制項Vx可以解決該類問題，則式（24）可表示為：

橢球約束和視場角約束同樣引入Ve和Vf，則式（26）和式（27）改寫為：

將Vx，Ve和Vf加入到最小燃料消耗的目標函數(shù)中，則有：

式中：Wx＞0，We＞0和Wf＞0 為相應的權重系數(shù)；在求解過程中，虛擬控制項可視作強制懲罰項以保證算法的收斂性。

為確保線性化后能夠捕獲原始問題的非線性，可以引入信賴域約束使前后兩步迭代之間沒有明顯的偏離［22］。在燃料最優(yōu)的目標函數(shù)中，狀態(tài)方程和約束的非線性主要是和xi(k)相關的，因此只引入狀態(tài)信賴域約束即可。設第i次和i-1次的狀態(tài)偏差為：

與狀態(tài)相關的信賴域約束可以表示為：

上述約束引入到目標函數(shù)后，式（32）可進一步改寫為：

式中：WΔ是信賴域約束的權重系數(shù)。

2.5 序列凸優(yōu)化算法

經(jīng)過上述凸化、離散化以及引入虛擬控制和信賴域后，小行星著陸軌跡凸優(yōu)化子問題，即二階錐規(guī)劃問題（SOCP）如下所示：

式中：x0和xf分別表示初始狀態(tài)和預計著陸狀態(tài)；ef∈R13表示允許著陸誤差；mwet和mdry分別表示探測器初始質量和探測器的凈質量。

設最小狀態(tài)偏差為ε，最大迭代次數(shù)為Niter，初始參考著陸參數(shù)表示為{x0(k)，u0(k)，m0(k)}，則小行星著陸軌跡序列凸優(yōu)化算法流程如圖4所示。

圖4 小行星著陸軌跡序列凸優(yōu)化算法流程Fig.4 Flow chart of sequence convex optimization algorithm for asteroid landing trajectory

首先利用初始狀態(tài)x0和預計著陸狀態(tài)xf計算初始參考著陸參數(shù)，然后利用其計算系數(shù)矩陣。在i小于最大迭代次數(shù)時，利用第i-1 次迭代的結果進行凸優(yōu)化求解，該過程需滿足狀態(tài)約束、控制約束及邊界條件，并且根據(jù)探測器位置矢量ri-1(k)實時更新引力加速度?U(ri-1(k))；優(yōu)化后可得到第i次迭代結果{xi(k)，ui(k)，mi(k)}，將第i-1 次結果和第i次結果進行比較，兩者小于一定范圍（最小狀態(tài)偏差ε），可認為已經(jīng)收斂并輸出優(yōu)化軌跡；否則以此次計算結果重新執(zhí)行凸優(yōu)化。以此類推，直到狀態(tài)偏差滿足最小狀態(tài)偏差或者達到最大迭代次數(shù)。

3 仿真驗證

為了快速進行序列求解，初始軌跡可利用初始狀態(tài)x0和預計著陸狀態(tài)xf均勻插值獲取，如下：

燃料消耗會導致慣性張量發(fā)生變化，定義慣性張量更新方式，如下：

式中：Ji(k)代表第i次迭代k時刻的慣性張量；mi(k)是第i次迭代k時刻探測器質量。

為了校驗所提出算法的有效性，本文以433Eros小行星為典型分析對象，它的自轉角速度和密度分別為3.31×10-4rad/s和2.67×103kg/m3。仿真采用了856 個頂點和1 708 個面的小行星多面體模型，相關數(shù)據(jù)可以在PDS 網(wǎng)站上獲?。?3］。時間參數(shù)設定為動力下降段至最終著陸段歷時tc=900 s 及最終著陸時間tf=1 200 s，允許狀態(tài)偏差為ε=10-3，其他參數(shù)如表1和表2所示。

表1 仿真參數(shù)Table 1 Simulation parameters

表2 初始狀態(tài)、預計著陸狀態(tài)及容許著陸誤差Table 2 Initial state，final landing state and tolerating landing errors

本文數(shù)值模擬使用了MATLAB 軟件下的CVX工具箱進行解算。根據(jù)上述模型和算法，小行星著陸軌跡序列凸優(yōu)化仿真結果如圖5～12 所示；為了直觀地表示探測器姿態(tài)變化，本文將四元數(shù)轉換為歐拉角（圖8）。

圖5 位置變化曲線Fig.5 Curves of position

圖6 速度變化曲線Fig.6 Curves of velocity

圖7 角速度變化曲線Fig.7 Curves of angular velocity

圖8 姿態(tài)變化曲線Fig.8 Curves of attitude

如圖5～8所示，所提出的算法可以引導探測器到達目標著陸點，探測器的最終著陸狀態(tài)誤差小于允許著陸誤差ef，利用其優(yōu)化出來的軌跡是有效的。

從圖9～10 可知，探測器的推力和力矩均滿足控制約束。探測器的質量變化如圖11所示，最終剩余質量為1 394.8 kg，也就是說消耗燃料5.2 kg，滿足質量約束。

圖9 推力變化曲線Fig.9 Curves of thrust

圖10 力矩變化曲線Fig.10 Curves of moment

圖11 質量變化曲線Fig.11 Curve of mass

式（10）表明橢球約束與探測器的位置相關，進而影響探測器著陸過程中的速度等參數(shù)；經(jīng)計算，探測器的位置在動力下降段[0，tc]滿足橢球約束。

式（11）表明視場角約束的作用是限制探測器姿態(tài)運動，即探測器要持續(xù)調(diào)整姿態(tài)才可以追蹤著陸點；如圖12 所示，探測器的視場角β在最終著陸段(tc，tf]滿足視場角約束。

圖12 視場角變化曲線Fig.12 Curve of angle of field of view

為了進一步驗證軌跡優(yōu)化算法對初始狀態(tài)偏差的魯棒性，本文進行了500次蒙特卡洛打靶仿真。初始位置和速度參數(shù)如表3 所示，其他未列舉的參數(shù)與表1～2中相同。

表3 蒙特卡洛打靶仿真初始位置和速度Table 3 Initial positions and velocities for Monte Carlo simulation

著陸誤差統(tǒng)計如表4 所示，各個參數(shù)的單位與表2一致，其中最終著陸誤差都達到了邊界值（允許著陸誤差ef）。由仿真結果可知，即使存在初始狀態(tài)偏差，在制導算法作用下，探測器仍能以較小的狀態(tài)誤差到達著陸點，表明所提出的制導算法有很強的魯棒性。

表4 著陸狀態(tài)誤差統(tǒng)計Table 4 Error statistics for landing state

4 結論

本文提出了一種6-DOF 小行星著陸軌跡序列凸優(yōu)化算法。通過線性化和離散化將小行星著陸軌跡優(yōu)化問題轉化為一個可以進行序列凸優(yōu)化求解的SOCP 問題；并引入虛擬控制項和信賴域，以保證線性化后采用凸優(yōu)化方法是可行的。以形狀不規(guī)則的433Eros 小行星為著陸目標，采用多面體引力場模型對該小行星進行引力場建模；仿真結果表明探測器能在節(jié)省燃料的同時進行高精度著陸，并滿足各項約束。所提出小行星著陸軌跡優(yōu)化算法可以擴展到6 自由度耦合的其他行星著陸問題，如火星或月球探測；同時，該方法精度高、魯棒性強等優(yōu)勢對機載應用非常有吸引力。