觸地關機模式下的著陸器軟著陸穩定性研究

董洋， 王春潔,2,*， 吳宏宇， 丁宗茂， 滿劍鋒

(1. 北京航空航天大學 機械工程及自動化學院, 北京 100083；2. 北京航空航天大學 虛擬現實技術與系統國家重點實驗室, 北京 100083)

腿式著陸器因其結構簡單、可靠性較高等特點，在各國地外星球探測任務中普遍被選擇。目前，腿式著陸器性能的研究和驗證方法主要有物理樣機實驗和虛擬樣機仿真。前者由于成本較高，且其他星球的重力環境在實驗過程中很難模擬，故難以研究著陸器在多種工況下的性能。隨著計算機運算能力的快速增長，后者的高效率和低成本等優點不斷凸顯，逐漸成為著陸器研制過程中不可或缺的環節。

文獻[1-3]提出了地球重力下的著陸沖擊實驗方法，研究了腿式著陸器在典型工況下的軟著陸性能。文獻[4-7]建立了著陸器的動力學仿真模型，研究了典型工況下的著陸器軟著陸性能，其中文獻[6]通過對比仿真與實驗結果，證明了建模的準確性。文獻[8-9]基于理論分析給出著陸參數穩定性邊界，其中文獻[9]通過實驗對邊界進行驗證。文獻[10-11]基于腿式著陸器非0高度關機軟著陸模式下的動力學仿真模型，結合離散化搜尋的仿真策略，給出了主要工況參數的二維穩定性邊界，其中文獻[11]在給出穩定性邊界的過程中引入了響應面方法，提高了計算效率。文獻[12-14]在腿式著陸器的動力學仿真模型中添加了主發動機與姿態控制發動機的推力等效模型，研究了著陸器的觸地關機軟著陸穩定性，但未給出直觀的穩定性邊界。

有鑒于此，本文以某型著陸器為研究對象，利用ADAMS軟件建立其觸地關機軟著陸模式下的動力學仿真模型，結合仿真模型與優化方法確定了多組極惡劣地形工況，并結合徑向基函數(RBF)神經網絡代理模型，給出了速度的三維穩定性邊界。

1 著陸器動力學建模

1.1 著陸器整機模型

本文研究的某型著陸器裝有4套中心對稱分布于主體的緩沖機構，如圖1所示，著陸器是由多個構件通過運動副連接且受到復雜的外力作用的多體系統，分析著陸器軟著陸性能，需基于多體動力學理論，建立著陸器的多體動力學模型。其基本動力學方程如下[15]:

(1)

(2)

(3)

式中：M為廣義質量矩陣，由著陸器各構件質量和轉動慣量構成；Ф為運動副約束方程矢量；q為廣義坐標矢量；Фq為約束方程的雅克比矩陣；Q為廣義力，主要考慮發動機推力和接觸力；λ為拉格朗日乘子矢量；t為時間。結合所建立的著陸器的動力學模型，利用ADAMS軟件進行仿真分析。

緩沖機構主要由主支柱、緩沖桿、左右輔助支柱和足墊組成[12-14]。緩沖桿在受到沖擊時會發生塑性變形，其建模方法參考文獻[12]的剛體等效法。如圖2所示，αe、Me和Se分別為緩沖桿的等效轉角、等效彎矩和等效軸向位移，先在有限元軟件中建立緩沖桿的模型，進行仿真分析，可以得到圖3所示的αe-Me和αe-Se曲線，然后將曲線導入剛體模型之中即可完成緩沖桿的等效建模。

輔助支柱內外筒之間的緩沖元件力學特性通過有限元分析獲得，進而確定輔助支柱的壓縮位移DS和緩沖力FS之間關系如圖4所示。將DS-FS曲線導入剛體模型之中完成其等效建模。

足墊和著陸面的接觸采用庫倫摩擦模型模擬切向力，采用非線性彈簧阻尼模型模擬法向碰撞力[12]。

圖2 緩沖桿等效模型Fig.2 Equivalent model of buffer rod

圖3 緩沖桿力學特性Fig.3 Mechanical behavior of buffer rod

圖4 緩沖元件力學特性Fig.4 Mechanical behavior of cushion element

1.2 發動機動力系統等效建模

著陸器主體上裝有一個用于制動的主發動機和多臺用于姿態調節的姿控發動機。觸地關機軟著陸模式下，發動機作用對著陸器的軟著陸穩定性不能忽視，因此需要建立考慮發動機作用效果的著陸器動力學仿真模型。

1.2.1 姿控發動機控制策略

除主發動機外，觸地關機軟著陸模式下還有10臺姿控發動機處于開啟狀態，各發動機的位置如圖5所示。其中，2臺位于主體上方，編號A1、A2，8臺位于著陸器主體下方，編號B1～B8。

圖5 坐標系定義及發動機位置Fig.5 Coordinate system definition and engine position

此外，圖5還包括了與姿態控制相關的坐標系。其中，O-XYZ為地面固定坐標系(右手系)，X軸沿重力方向指向天空，Z軸水平指向下坡方向。Oc-XcYcZc為隨體控制坐標系(右手系)，原點在著陸平臺中心，Xc軸垂直著陸平臺表面且指向天空，Zc軸在著陸平臺內且在俯視圖中指向相鄰兩足墊的中點。定義O-XYZ按照Z-Y-X的順序旋轉到與Oc-XcYcZc平行位置時繞三軸轉動的角度分別為θx、θy、θz；同時定義Oc-XcYcZc相對O-XYZ的角速度矢量在Oc-XcYcZc三軸上的投影分別為ωx、ωy、ωz。將上述6個參數作為控制輸入。

為提升著陸器姿態的穩定性，A1、A2發動機的推力如式(4)和式(5)所示；B系列發動機的推力公式形式相同，均參考文獻[12]建立，以B1、B2為例，其推力如式(6)所示：

FA1=

(4)

FA2=

(5)

FBi=

i=1,2

(6)

式中：ai、bi、ki為控制策略中的相關參數；FM1與FM2分別為2 個系列發動機推力最大值。

1.2.2 主發動機控制策略

主發動機位于主體下表面中心且噴口方向垂直主體下表面。正常工作時，其推力方程為[12]

(7)

式中：m為著陸器整機質量；vx為豎直速度值。

觸地關機軟著陸模式下，若足墊與星球表面接觸，則其內部的傳感器將發出主發動機關機信號，但主發動機由于關機后效其推力不能馬上減小到零，該情況會對著陸穩定性產生影響[16]，因此，為了提高著陸穩定性，姿控發動機會在主發動機收到關機指令之后延遲關機。關機信號發出后，由于信號傳輸的延遲，主發動機推力下降開始時間會滯后。其推力變化過程如圖6所示。圖中：t0為觸地時刻，t1時刻發出關機信號，t2時刻推力開始下降，t3時刻姿控發動機關閉。

圖6 主發動機推力變化過程Fig.6 Thrust changing process of main engine

1.2.3 發動機推力等效建模

本文通過在發動機安裝位置添加與發動機噴口方向相反的力來模擬發動機的推力作用效果，并利用測量函數實時反饋θx、θy、θz、ωx、ωy、ωz與vx的值。在真實情況下，由于傳感器按給定周期采集和處理信號，故推力值不隨位姿和速度參數連續變化[17]，因此，本文在根據式(7)計算推力值的同時，也按照給定周期給模型中的推力賦值，以模擬發動機的真實工作情況。

2 著陸工況參數與軟著陸穩定性判據

2.1 著陸工況參數

本文考慮的著陸工況參數有2類：

1) 地形參數。包括著陸面坡度(α)、足墊滑移狀態(K)和著陸偏航角(ψ)。由于著陸器安裝有地貌識別系統，正常情況下會降落在平原地區，坡度一般不超過12°，本文考慮著陸面坡度為0°～12°。正常情況下足墊可以在著陸面上滑移，但如果足墊遇到石塊或者小坑，其將處于卡死狀態，本文用K=0表示足墊可滑移，K=1表示足墊卡死。考慮著陸器的幾何對稱性，ψ取0°～45°即可以表征所有用于評判著陸器軟著陸穩定性的姿態。

2) 著陸速度。指著陸器相對地面固定坐標系的速度，為空間矢量。為方便后文的研究，著陸速度采用柱坐標系(vh，β，vx)描述，vh(0～2.5 m/s)為速度在ZOY平面投影的模長，即水平速度大小；β(0°～360°)為從+X方向看，Z軸繞到水平速度方向轉過的角度；vx為豎直速度大小(0～5 m/s)，其方向沿-X方向。

部分著陸工況參數如圖7所示。

圖7 部分著陸工況參數示意圖Fig.7 Schematic diagram of partial landing conditions parameters

2.2 軟著陸穩定性判據

本文重點研究著陸器的軟著陸穩定性，根據著陸器構型，確定其穩定性評判參數如下：

1) 緩沖桿最大等效轉角αem(單位為(°))。緩沖桿通過塑性變形吸收沖擊能量，若等效轉角過大其極易發生損壞，故αem應小于許用值α*。

2) 輔助支柱最大壓縮行程DSM(單位為mm)。一方面應避免DSM超過設計壓縮范圍，另一方面，DSM過大易導致著陸平臺和坡面夾角過大，不利于探測車的轉移，故DSM應小于許用值D*。

3) 著陸器質心到4個翻倒平面距離的最小值LD(單位為mm)。定義通過相鄰兩個足墊中心的豎直平面為翻倒平面[11]，為保證著陸器在軟著陸過程中不發生翻倒，LD應大于許用值L*。

4) 主發動機尾噴管底面中心到著陸面的最小距離HM(單位為mm)。觸地關機軟著陸過程中，主發動機大部分時間處于開啟狀態，考慮HM過小容易導致粉塵堵塞尾噴管或凸起撞擊尾噴管等情況，故HM應大于許用值H*。

3 極惡劣地形的確定

著陸工況參數中，α、K、ψ依賴于星球的地形，難以主動預測或者控制，而vx、vh、β作為著陸器自身的參數，具有一定的可控性，為了研究可控參數對著陸器軟著陸穩定性的具體影響方式，有必要事先確定α、K、ψ的最壞取值。基于動力學模型，本文采用優化方法尋找考慮著陸器水平速度、使穩定性指標αem、DSM、LD、HM最差的α、K、ψ取值。為了降低優化維度、提高計算效率，設定vx=2.5 m/s、vh=1.0 m/s，考慮水平速度大小固定的情況下，β取值的不同可能帶來α、K、ψ最差取值的不同，故將β作為尋找最差工況的變量之一。綜上，以尋找LD最差的工況為例，建立優化模型如下：

minLD(X)X=[α,K,ψ,β]T∈Ω

(8)

式中：Ω為參數的取值范圍。

本文使用多島遺傳算法(MIGA)、結合動力學模型實現優化計算，算法參數配置如表1所示。

經優化計算，得到4組極惡劣工況如表2所示。

表1 優化參數設置Table 1 Setting of optimal parameters

表2 極惡劣工況Table 2 Worst landing conditions

4 可控工況參數分析

4.1 單穩定性指標的速度穩定性邊界

研究vx、vh、β對著陸器軟著陸穩定性的具體影響方式的過程中，需進行大規模仿真計算，為了提高效率，本文引入RBF神經網絡代理模型建立極惡劣工況下的vx、vh、β與4項軟著陸穩定性指標值之間的映射關系[18]。通過隨機抽點仿真，將仿真結果導入神經網絡模型進行訓練，用400個獨立的樣本測試模型，訓練終止準則為均方根相對誤差值RMSE(見式(9))小于0.05。

(9)

得到神經網絡模型后，本文采用控制變量法，令vx、vh、β在取值范圍內等間距抽取100個點，共得到106組樣本點，將其代入神經網絡模型中，計算得到106組軟著陸穩定性指標值。綜合考慮著陸器的軟著陸可靠性與模型誤差等因素，設定本文的各項軟著陸穩定性指標許用值分別為α*=30°、D*=70 mm、L*=50 mm、H*=250 mm，進而給出了不同vx取值下反映vh、β與αem、DSM、LD、HM關系的云圖如圖8所示，vh單位為m/s。

此外，本文還在柱坐標系下給出了vx、vh、β三維穩定性邊界(許用值等值面)，如圖9所示。

1) 緩沖桿最大等效轉角分析定性指標的速度三維穩定性邊界如圖9(a)所示，圖中曲面為αem=30°的等值面。

由圖8(a)可知，在相同vx與vh下，β對αem影響非常顯著，β越接近180°、αem越大。因此，圖8(a)所示的αem許用值等值線接近直線，并且集中在120°<β<240°范圍內，即偏坡上方向；另外隨著vx增大，不穩定區域面積增大。以αem為穩

2) 輔助支柱最大壓縮行程分析

由圖8(b)可知，在相同vx下，DSM的許用值等值線接近環形，可知β對DSM影響比較小，DSM的主要影響因素為vh。此外，隨著vx增大，不穩定區域面積增大。以DSM為穩定性指標的速度三維穩定性邊界如圖9(b)所示，圖中曲面為DSM=70mm的等值面。

3) 質心距翻倒平面最小距離分析

由圖8(c)可知，LD的許用值等值線集中在270°<β <360°方向，即偏坡下方向，表明水平速度方向接近坡下方向時著陸器更容易翻倒。以LD為穩定性指標的速度三維穩定性邊界如圖9(c)所示，圖中曲面為LD=50mm的等值面。

4) 尾噴管底面中心到著陸面最小距離分析

由圖8(d)可知，在相同vx下，HM的許用值等值線非常光滑且接近圓，說明β對HM基本沒有影響；以HM為穩定性指標的速度三維穩定性邊界如圖9(d)所示，圖中曲面為HM=250mm的等值面，圖9(d)直觀顯示邊界集中在vx=4～5m/s區間，一方面表明該指標主要受vx影響，同時vx小于4m/s時，HM不會超過安全閾值。

4.2 綜合多項穩定性指標的速度穩定性邊界

將圖9所示的全部三維穩定性邊界向同一柱坐標系內整合，得到圖10所示的綜合三維穩定性邊界，同時得到不同豎直速度下的切面圖如圖11所示，圖中穩定區域為各穩定區域的交集,vh單位為m/s。

將圖10所示的三維穩定性邊界向同一豎直面內做旋轉投影，得到消除β影響的4項軟著陸穩定性指標的保守速度穩定性邊界，綜合4條保守邊界，可最終確定綜合多項穩定性指標的速度穩定性邊界，如圖12所示。由圖12可以直觀地確定保證著陸器絕對穩定著陸的速度取值范圍。

圖10 綜合指標下著陸速度三維穩定性邊界Fig.10 3D stability boundary determined by landing velocity considering comprehensive indicators

圖11 不同vx下綜合穩定性邊界Fig.11 Comprehensive stability boundary under different vx

圖12 vx-vh綜合穩定性邊界Fig.12 vx-vh comprehensive stability boundary

5 結 論

本文建立某型著陸器觸地關機模式下的動力學仿真模型，著重分析著陸速度對穩定性的影響，得到：

1) 結合動力學仿真模型及多島遺傳優化算法，確定了考慮著陸器水平速度下的極惡劣地形參數。

2) 基于仿真模型計算結果，利用RBF神經網絡建立了反映極惡劣地形下著陸器速度參數與穩定性指標值之間映射關系的代理模型，進而計算得到了軟著陸穩定性指標的速度穩定性邊界，確定了保證著陸器安全軟著陸的速度取值范圍。

3) 本文研究方法具有低成本、高效率等特點，且適用于其他型號著陸器的研制過程。