火星著陸任務落點誤差快速分析方法

任高峰，朱圣英

(1.哈爾濱工業大學航天學院，150080 哈爾濱;2.北京理工大學宇航學院，100081 北京)

火星著陸任務落點誤差快速分析方法

任高峰1，朱圣英2

(1.哈爾濱工業大學航天學院，150080 哈爾濱;2.北京理工大學宇航學院，100081 北京)

在火星著陸任務設計中，為高效地評估進入點導航誤差等動力學方程中的不確定因素對著陸精度的影響，引入Galerkin投影法則，提出一種火星著陸落點不確定度快速預報與分析方法.該方法采用Wiener-Askey正交多項式逼近著陸器狀態，將落點誤差的預報問題轉換為1個確定性微分方程的求解問題.最后以我國未來實施火星著陸探測為背景，分別用本文算法和蒙特卡洛方法進行數值仿真分析，結果表明，本文算法能夠對落點偏差進行快速分析，準確預測其均值和方差，且計算效率明顯提高.

落點誤差;正交多項式;火星著陸;Galerkin投影

對于未來的火星探測任務，為了提高科學回報和降低任務成本，選擇在科學價值高、地形復雜的區域進行著陸成為1個必然的要求和發展趨勢［1］.然而，探測器在火星大氣進入點處的導航控制誤差，探測器的氣動參數以及火星大氣模型的不確定性，都會嚴重影響著陸器最終的著陸精度，甚至關乎任務的成敗［2］.因此，分析這些偏差以及不確定性對著陸點的影響，是1項必不可少的工作;另外，隨著著陸精度要求的提高，在設計著陸器進入大氣的標稱軌跡時，同時考慮了不確定性因素影響的魯棒標稱軌跡設計，成為未來標稱軌跡設計發展的1個重要方向，而不確定性輸入在系統中如何演化，是魯棒軌跡規劃的重要研究內容［3］.綜上可見，針對火星著陸任務，發展一種快速的落點不確定度分析方法，對降低未來火星著陸設計周期和成本，提高設計效率很有意義.

對于動力學系統中不確定性的傳播問題，目前存在的方法主要有蒙特卡洛仿真方法、線性化方法、無味變換方法以及Galerkin投影法等.文獻［4］對蒙特卡洛方法進行了詳細的闡述，這種方法的優點在于計算實施簡單，可以精確得到狀態量精確的統計特性，不足之處在于它的收斂速度比較慢;針對特定的問題，文獻［5］提出了準蒙特卡洛方法，文獻［6］提出了馬爾科夫鏈蒙特卡洛方法，從一定程度上提高了計算效率;針對不確定性輸入為高斯白噪聲，動力學系統非線性強度在可接受范圍內的情況，文獻［7］利用線性化的方法將非線性系統泰勒展開，以線性理論為基礎，以誤差方差陣為表現形式，實現不確定性的遞推.文獻［8］提出一種確定的取樣方法，利用一系列sigma點來逼近系統輸入的概率分布密度，然后將非線性系統作用于sigma點，利用轉換后的點來近似系統輸出的統計特性，所取sigma點個數會隨著系統維數及精度要求的提高而迅速增加.Galerkin投影法以混沌多項式展開為基礎，因其具有嚴密的數學理論體系，成為目前研究的1個熱點，在隨機流體力學［9］，隨機有限元分析［10］，固體力學［11］，空間技術［12］等領域得到了應用，但其在動力學控制與動態估計中應用很少.

本文基于Galerkin投影法，對火星著陸系統落點偏差分析問題進行了研究.首先對火星著陸器動力學系統及不確定參數進行描述，然后利用Wiener-Askey正交多項式表示著陸器的狀態和不確定參數，根據Galerkin投影法則，將隨機動力學系統轉化為1個高維的確定性系統，最后采用牛頓積分求解正交多項式的時變系數，進一步得到著陸器狀態的均值和方差.

1 問題描述

忽略火星自轉，假設火星表面大氣相對于火星靜止，則火星著陸縱向動力學可以描述為

其中:h表示著陸器距離火星表面的距離;v表示著陸器速度的大小;γ表示航跡角;μ表示火星引力系數;Rm表示火星半徑;B表示著陸器的彈道系數;k表示著陸器的升阻比;φ表示傾側角;λ表示大氣模型不確定性因子;ρ表示火星大氣密度，其與著陸器距離火星表面高度的關系如式(2)所示，它是根據NASA開發的火星大氣模型Mars-Gram所生成的數據進行最小二乘擬合得到的［3］.

由于測控能力的限制，著陸器在進入點處的位置、速度等狀態不能精確已知，并且著陸器的氣動參數依靠地面試驗也不能夠精確建模.對應的，本文假設上述動力學方程中的狀態初值和系統參數均存在不確定性，表示為

其中 δh0、δv0、δγ0、δB、δk、δλ 為符合一定分布的隨機量.

根據以上描述，火星著陸系統落點偏差的分析問題可以抽象成為在初始狀態和系統參數均存在不確定性的隨機微分方程的求解問題.下面將利用隨機投影法對上述問題進行求解.

2 Galerkin投影法

2.1 Galerkin投影法概述

Galerkin 投影法首先由 Wiener［13］提出，用于對符合高斯分布的隨機過程進行建模分析;學者修東斌和 Karniadakis［14］利用 Askey正交多項式體系以及Cameron Martin定理將其推廣到符合多種分布的隨機過程，并對這一推廣的收斂性進行了證明;隨后這種方法在眾多領域得到了廣泛的應用［9－12］.本文利用上述成果，將帶有不確定性的著陸器狀態用一組正交多項式逼近，此多項式的自變量為符合一定分布的隨機變量;隨后將多項式表示的狀態帶入動力學方程中，形成1個顯式表達的隨機微分方程，最后利用Galerkin投影法則，將其轉換為1個高維的確定性微分方程，進而通過對此確定性微分方程求解實現對著陸器狀態不確定性進行分析.

2.2 Askey正交多項式體系

根據Cameron Martin定理及其推廣，任何二階隨機過程X(t)可以用式(3)在L2意義上進行近似逼近:

其中In是以隨機變量ξ為自變量的n階混沌多項式，為表示及后續計算方便，上式可以重新表示為以下形式:

其中 Φi(ξ)與 In(ξi1(θ)，…，ξi3(θ))存在一一映射關系.

Wiener Askey混沌多項式存在以下正交關系:

其中δij為克羅內克δ函數，并且〈·，·〉表示整體內積，假設W(ξ)為對應正交多項式的權重函數，對應于連續隨機變量和離散隨機變量，整體內積可以分別定義為

或者

需要特別指出的是，當權重函數對應于某一隨機變量的概率密度函數時，則此隨機變量可以用對應于相應權重函數的一階正交多項式精確表示.并且，對應于常用的隨機變量分布類型，如高斯分布等，與一些正交多項式存在著一一對應關系，如表1所示.

表1 正交多項式與隨機變量類型之間對應關系

2.3 高維確定性微分方程

利用Askey混沌多項式將式(1)～(2)描述的隨機微分方程轉換為確定性微分方程，首先要解決的問題是將系統的不確定性輸入用混沌多項式表示，由于系統的不確定性參數(彈道系數B，升阻比k，以及火星大氣不確定性因子λ)只與任務實施前對著陸器氣動特性及火星大氣的認知程度有關，在整個著陸過程中統計特性并不改變，本文假設這些不確定參數符合高斯分布，則可以用一階Hermite正交多項式表示為

對于系統的狀態(高度h、速度v以及航跡角γ)，在初始時刻通常以高斯分布的形式給出，但在著陸過程中，這些狀態量之間相互影響，并且受不確定參數的影響，在著陸器著陸過程中，狀態量的統計特性隨時間發生變化，對于每一時刻，系統狀態對應的隨機變量可以用混沌多項式表示，因此對于整個時間序列中的系統狀態，可以利用帶有時變系數的P階正交多項式表示:

將式(4)～(5)帶入到方程(1)～(2)描述的隨機微分方程中，可以得到

為了使上式更加簡便，做以下假設:

并利用Galerkin投影法則將式(6)投影到正交多項式的基Hm(ζ)上，則可以得到

其中，m=0，1，…，P.方程(10)為與隨機微分方程(1)等價的高維確定性微分方程，可以通過龍哥－庫塔等算法對此微分方程進行求解，并通過方程(5)可以得到系統狀態的統計特性.

由此，xi可以通過以下方程組求解:

其中M的元素為

2.3 狀態統計特性的求取

在利用式(10)求得混沌多項式的系數之后，著陸器的各狀態量可以用式(5)表示，然后根據定義，可求得各個狀態量的均值和方差，由于結構類似，本文僅給出高度的均值和方差計算公式，其它狀態量的方差和均值可以類比得到，高度的均值為

另外，需要指出兩點:1)本文方法不僅能得到隨機變量的方差和均值，而且可以根據定義通過式(5)求得隨機變量的任意階矩;2)本文僅針對不確定性輸入符合高斯分布的情況進行了分析，如果不確定性輸入服從其他標準分布，可以根據表1選擇對應的正交多項式來代替上述公式中的Hermite多項式，最后所得統計特性計算公式形式與式(10)～(12)完全一樣.

3 數學仿真與分析

為驗證所提算法的可行性，同時分析各種不確定性因素對落點偏差的影響，本文針對不同的初始誤差，分別利用蒙特卡洛仿真方法和本文方法進行數學仿真，著陸器氣動參數與標稱軌跡初始條件如表2所示.

表2 著陸器相關參數與標稱軌跡初始條件

假設系統中的彈道系數、升阻比以及大氣密度建模均存在10%的不確定性，而初始狀態高度存在5 km的不確定性，初始速度存在5 m/s的不確定性，初始航跡角存在0.03°的不確定性，同樣假設符合高斯分布.

仿真采用二階Hermite多項式，蒙特卡洛仿真進行10 000次，仿真計算機CPU為Q8400，主頻2.66 GHz，內存3.49 G，仿真軟件用Mathworks公司的Matlab2009a.

從圖1～圖6可以看出，在各種仿真條件下，本文基于混沌多項式展開的方法(PCE)與蒙特卡洛仿真方法(MC)所得到的結果都能夠很好吻合，為了更清楚地評價本文方法結果的精確程度，定義相對誤差:

圖1 考慮高度不確定性

圖2 考慮速度不確定性

圖3 考慮航跡角不確定性

圖4 考慮大氣密度不確定性

圖5 考慮升阻比不確定性

圖6 考慮彈道系數不確定性

表3 均值－標準差估計誤差及相對誤差

從表3可以看出，所有狀態的均值估計相對誤差均小于7%，標準差估計誤差小于12%，但蒙特卡洛仿真所消耗的時間為2 360.7 s，本文算法消耗的時間為15.09 s，計算效率明顯提高.另外，從圖4和圖6可以看出，僅考慮10%的大氣密度不確定性和僅考慮10%的彈道系數不確定性時，系統狀態量統計特性變化趨勢一致，原因分析如下:

根據式(1)可知，大氣密度不確定因子和彈道系數在動力學方程中總是同時出現，只是所處的位置不同.僅考慮大氣密度不確定性時，相當于在對應項處乘以隨機變量1+0.1ζ，僅考慮彈道系數不確定性時，相當于在對應項處除以隨機變量1+0.1ζ，其中ζ是符合標準正態分布的隨機變量，對于小量0.1ζ有

又因為對于零均值的標準正態分布ζ來說，1+0.1ζ和1-0.1ζ符合相同的分布，所以僅考慮大氣密度不確定性和僅考慮彈道系數不確定性兩種情況結果類似.

4 結論

本文針對火星著陸落點誤差評估這一問題做了相關研究，主要工作和研究結論如下:

1)提出了一種基于Galerkin投影的火星著陸落點誤差分析方法，該方法以Askey正交多項式和Galerkin投影法則為理論基礎，將描述著陸器動力學方程的隨機微分方程轉換成為1個等效的高維確定性微分方程，通過求解此確定性微分方程，實現對著陸器狀態不確定性進行快速分析.

2)以系統不確定性符合高斯分布為例，分別針對6種情況對算法的有效性進行了仿真驗證，通過和蒙特卡羅仿真結果對比可以看到:利用2階Hermite多項式對系統狀態進行逼近，則本文算法對狀態的均值估計誤差小于7%，標準差估計誤差小于12%，并且計算效率大為提高;

3)通過仿真發現大氣密度不確定性和彈道系數不確定性所引起的著陸器狀態統計特性變化趨勢大體一致，并對產生這種現象的機理在理論上進行分析;

4)這種算法可以推廣到系統不確定性為非高斯白噪聲的情況，只需相應的調整概率密度函數，選擇合適的正交多項式代替仿真中的Hermite多項式，算法流程及推導的等效高維確定性微分方程形式不變.

［1］BRAUN R D，MANNING R M.Mars exploration entry，descent，and landing challenges［J］.Journal of Spacecraft and Rockets，2007，44(2):310 －313.

［2］WOLF A，GRAVES C，POWELL R，et al.Systems for pinpoint landing at Mars［C］//14th AIAA/AAS Space Flight Mechanics Meeting.Pasadena:Jet Propulsion Laboratory，National Aeronautics and Space Administration，2004:1－20.

［3］SHEN Hai-jun，SEYWALD H，POWELL R.Desensitizing the pin-point landing trajectory on Mars［C］//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Pasadena:Jet Propulsion Laboratory，National Aeronautics and Space Administration，2008:1130－1143.

［4］FISHMAN G S.Monte carlo:concepts，algorithms and applications［M］.New York:Springer-Verlag，1996:66 －70.

［5］FOX B.Strategies for quasi-monte carlo［M］.Boston:Kluwer Academic Publishers，1999:305 －323.

［6］GAMERMAN D，LOPES H.Markov chain monte carlo:stochastic simulation for bayesian inference［M］.New York:Chapman＆ Hall/CRC，2006:121－136.

［7］BROWN R G，HWANG P Y C，Introduction to random signals and applied Kalman filtering:with MATLAB exercises and solutions［M］.Wiley:Chichester，1992:335 －377.

［8］WAN，E，VAN DER MERWE R.The unscented Kalman filter for nonlinear estimation［C］//Adaptive Systems for Signal Processing，Communications，and Control Symposium.Piscataway:IEEE Inc，2000:153－158.

［9］XIU Dong-bin，KARNIADAKIS G E.Modeling uncertainty in flow simulations via generalized polynomial chaos［J］.Journal of Computational Physics，2003，187(1):137－167.

［10］SPANOS P，GHANEM R.Stochastic finite element expansion for random media［J］.Journal of Engineering Mechanics，1988，115(5):1035 －1053.

［11］GHANEM R.Ingredients for a general purpose stochastic finite elements implementation［J］.Computer Methods in Applied Mechanics and Engineering，1999，168(1):19－34.

［12］PRABHAKAR A，FISHER J，BHATTACHARYA R.Polynomial chaos-based analysis of probabilistic uncertainty in hypersonic flight dynamics［J］.Journal of Guidance，Control，and Dynamics，2010，33(1):222 －234.

［13］WIENER N.The homogeneous chaos［J］.American Journal of Mathematics，1938，60(4):897－936.

［14］XIU Dong-bin，KARNIADAKIS G E.The wiener-askey polynomial chaos for stochastic differential equations［J］.SIAM Journal on Scientific Computing，2002，24(2):619－644.

A new landing site uncertainty analysis method for mars entry mission

REN Gao-feng1，ZHU Sheng-ying2

(1.School of Astronautics，Harbin Institute of Technology，150080 Harbin ，China;2.School of Aerospace Engineering，Beijing Institute of Technology，100081 Beijing，China)

To evaluate how the uncertain factors in the entry dynamics such as the navigation errors at the entry point affect the landing precision effectively，by introducing the Galerkin projection，a rapid landing error prediction method for Mars entry is proposed.In this method，the Polynomial Chaos are used to approximate the vehicle＇s states，and the problem is converted to a deterministic dynamical systems.Applying the method to a representative entry scenario of our country＇s future Mars entry mission，simulation results indicate that，compared with Monte Carlo method，the proposed method is able to predict the mean and covariance with little error and more computational efficiency.

landing error;orthogonal polynomial;mars landing;galerkin projection

V448.15

A

0367－6234(2012)07－0014－07

2011－06－23.

國家高技術研究發展計劃重點資助項目(2010122206).

任高峰(1983—)，男，博士研究生.

任高峰，stone－breaker@126.com.

(編輯 張 宏)