王文強

(北京空間機電研究所，北京 100094)

1 引言

作為深空探測的主要內(nèi)容之一，火星探測是一個國家航天科技實力的集中體現(xiàn)[1]。從美國和俄羅斯/蘇聯(lián)幾十年的火星探測任務(wù)中可以發(fā)現(xiàn)，實現(xiàn)軟著陸探測是獲取火星相關(guān)資料的主要的探測方式和目標(biāo)。而要開展火星表面軟著陸探測，首先需要獲取任務(wù)設(shè)計所必須的大氣密度、風(fēng)速隨高度的變化曲線，建立初步的火星表面大氣模型[2]。

目前，我國對火星大氣的了解幾乎完全基于國外公開發(fā)表的相關(guān)文獻資料?？紤]到文獻來源的可靠性以及火星大氣參數(shù)的不穩(wěn)定性，通過火星探測獲取數(shù)據(jù)建立火星大氣模型是我國開展進一步火星探測的必要步驟。

另一方面，由于火星著陸任務(wù)中探測器飛行時間較長，同時再考慮到火星大氣的不確定性等因素，實際著陸彈道將嚴重偏離標(biāo)準彈道。因此記錄著陸過程的氣動參數(shù)并重建進入彈道是對火星探測器進行任務(wù)分析和性能評估的重要措施，同時也為后續(xù)探測器任務(wù)設(shè)計提供參考依據(jù)。

本文以“火星科學(xué)實驗室”(Mars Science Laboratory, MSL)為研究對象，討論了火星探測中大氣模型以及進入彈道的重建方法，為我國未來將要開展的火星探測以及建立火星大氣模型提供依據(jù)。

2 MSL數(shù)據(jù)獲取系統(tǒng)介紹

為了測量進入艙在進入火星大氣時的大氣環(huán)境和熱環(huán)境，MSL進入艙安裝了1套火星進入減速著陸儀器系統(tǒng)(Mars Entry, Descent, and Landing Instrumentation, MEDLI)[3-4]，如圖1所示，圖中不同顏色代表預(yù)估的防熱罩表面壓力大小，紅色處壓力最大，綠色最小。

圖1 MEADS測量單元安裝位置Fig.1 MEADS port geometry

MEDLI由安裝在MSL防熱罩上的1組傳感器組成，依據(jù)傳感器類型和主要功能可以分為3個部分，分別為傳感器集成插頭(Mars Integrated Sensor Plug, MISP)，數(shù)據(jù)采集系統(tǒng)(Mars Entry Atmospheric Data System, MEADS)以及傳感器電子支持設(shè)備 (Sensor Support Electronics, SSE)。MISP系統(tǒng)由7個安裝在防熱罩里的集成感應(yīng)塞組成，每個感應(yīng)塞包括1個衰退傳感器和4個熱電偶傳感器，用于測量進入大氣時的熱環(huán)境，得到防熱罩表面溫度分布情況。MEADS由 7個安裝在防熱罩上的壓力測量單元組成(圖 1中P1～P7)，測量飛行時防熱罩表面的壓力分布，在測得進入艙速度的情況下，MEADS的測量數(shù)據(jù)還可以用于估計當(dāng)?shù)氐拇髿饷芏纫约帮L(fēng)速。MEADS的工作頻率為 8Hz，相應(yīng)的采樣周期為 0.125s，從進入大氣時開始啟動，持續(xù)工作至防熱罩分離。SSE的任務(wù)是為MISP和MEADS提供輸入信號，并保存記錄相應(yīng)的測量數(shù)據(jù)。本文主要介紹MEADS的組成以及相應(yīng)的估計模型和算法。

根據(jù)MSL任務(wù)介紹，MEADS的主要任務(wù)是單獨利用壓力的測量值估計進入彈道的氣動參數(shù)，并要求對攻角和側(cè)滑角的估計偏差不超過 0.5°(3σ)，動壓估計偏差不超過 2%(3σ)，馬赫數(shù)估計偏差不超過0.1(3σ)[5-6]。另一方面，利用機載慣性測量單元(Inertial Measurement Unit, IMU)測量得到艙體角速度和加速度數(shù)據(jù)后，還可以結(jié)合估計得到的氣動參數(shù)實現(xiàn)對來流密度以及風(fēng)速的估計，從而建立火星大氣模型。

3 數(shù)學(xué)模型

3.1 表面壓力模型

傳統(tǒng)的大氣參數(shù)重建都是基于速管大氣數(shù)據(jù)系統(tǒng)模型(Flush Atmospheric Data System，F(xiàn)ADS)進行參數(shù)估計，一般采用牛頓流體近似理論建立表面壓力分布模型。其中表面壓力是測量單元位置、攻角、側(cè)滑角、馬赫數(shù)等參數(shù)的函數(shù)。圖2即為某個測量單元處速度關(guān)于測壓管軸線的入射角示意圖。

圖中，θ為入射角;v為航天器速度。測壓管測得的壓力包括當(dāng)?shù)氐撵o壓以及沿著測壓管軸線方向的動壓，即

圖2 入射角示意圖Fig.2 Incidence angle

式中p為某個測量單元的測量值;p∞為當(dāng)?shù)仂o壓;ρ為當(dāng)?shù)卮髿饷芏?q∞為動壓。整理后可得

式中pt為當(dāng)?shù)氐目倝?R為靜壓與總壓之比，是馬赫數(shù)的函數(shù)，具體表達式為:

式中γ為比熱比;Ma為馬赫數(shù)。

對于火星大氣，取γ=1.335，馬赫數(shù)與R的函數(shù)關(guān)系如圖3所示:

圖3 γ=1.335時馬赫數(shù)與R的函數(shù)關(guān)系Fig.3 Mach vs pressure ratio R forγ=1.335

入射角θ是壓力測量單元位置以及攻角和側(cè)滑角的函數(shù)。圖4為探測器軸向和側(cè)向視圖，顯示了其中1個測量單元在防熱罩上的安裝位置:

圖4 測壓管安裝位置Fig.4 Pressure port geometry

圖中，壓力測量單元的安裝位置由安裝錐角λ和指針角φ確定。對于每個測量單元，錐角和指針角均為確定的常數(shù)。當(dāng)進入艙的攻角為α，側(cè)滑角為β時，對圖4中的測壓單元進行分析，如圖5所示:

圖5中O-XYZ為進入艙本體坐標(biāo)系，X軸沿著進入艙縱軸指向前，Z軸在縱向?qū)ΨQ面內(nèi)，垂直于X軸指向上，Y軸由右手坐標(biāo)系確定。OA為某個壓力測量單元的軸線方向，與進入艙的夾角即為錐角λ，OD為進入艙的速度v，OC為速度在OXZ平面的投影，α和β分別為攻角和側(cè)滑角，速度與壓力測量單元的軸線夾角θ即為入射角，面AOB與面BOC之間的夾角即為上文定義的指針角φ。

圖5 α,β,λ和θ之間的角度關(guān)系Fig.5 Relationship betweenα,β,λandθ

在四面角O-ACD、O-ACB中分別利用四面角的余弦和正弦定理可以得到入射角與測量單元位置、攻角以及側(cè)滑角的函數(shù)關(guān)系為:

同理，對于其他測量單元，也可以得到與上式類似的關(guān)系式(只是安裝錐角λ和指針角φ數(shù)值不同)。

將公式(3)和公式(4)代入公式(2)就可以得到某個測量單元測得的壓力值關(guān)于測量點位置、總壓、馬赫數(shù)、攻角以及側(cè)滑角的函數(shù)關(guān)系，將7個測量單元的關(guān)系式組成1個函數(shù)方程組，即為表面壓力分布模型。

綜上，基于傳統(tǒng)的速管大氣數(shù)據(jù)系統(tǒng)模型，通過推導(dǎo)測壓單元測量值與安裝位置之間的函數(shù)關(guān)系，建立了MSL防熱罩表面的壓力分布模型，用于下文的狀態(tài)估計。

3.2 卡爾曼濾波

上文已經(jīng)建立了表面壓力分布模型，在此基礎(chǔ)上利用卡爾曼濾波估計氣動參數(shù)和大氣參數(shù)。卡爾曼濾波器可以將MEADS和IMU測量結(jié)果通過一種緊耦合的方式直接融合，從而利用MEADS測量的壓力值和IMU測量得到的加速度及角速度值直接估計大氣數(shù)據(jù)和氣動參數(shù)，其優(yōu)越性主要體現(xiàn)在可以通過反復(fù)迭代基本消除測量噪聲的影響，從而利用有限的量測值估計得到精度較高的狀態(tài)量。實現(xiàn)的方法是將壓力模型改寫為來流氣動特性(壓力、密度、風(fēng)速)和空速的函數(shù)，構(gòu)建量測方程，計算得到這些中間變量之后，攻角、側(cè)滑角、馬赫數(shù)以及動壓也可以進一步計算得到[7]。

基于這種思路，濾波器選擇的狀態(tài)變量x為:

式中r為進入艙質(zhì)心到火星質(zhì)心距離;φ是赤維;ei為水平坐標(biāo)系到本體坐標(biāo)系的變換四元數(shù)[8];u,v,w分別為進入艙地速在當(dāng)?shù)厮阶鴺?biāo)系中的投影;uw,vw,ww為風(fēng)速在當(dāng)?shù)厮阶鴺?biāo)系中的投影。

將圖5中進入艙空速投影到本體坐標(biāo)系上，沿3個軸的投影值分別為ub,vb,wb。則氣動參數(shù)(攻角、側(cè)滑角)可以表示為:

ub,vb,wb同時也是風(fēng)速和進入軌道的函數(shù):

式中Ψ,,ΘΦ分別為本體坐標(biāo)系相對于當(dāng)?shù)厮阶鴺?biāo)系(北東下坐標(biāo)系)的偏航角、俯仰角、滾轉(zhuǎn)角;G為相應(yīng)的轉(zhuǎn)移矩陣。

聯(lián)立公式(1)和公式(4)可得:

上式中再利用式(5)和式(6)中攻角、側(cè)滑角與當(dāng)?shù)厮阶鴺?biāo)系中u,v,w,uw,vw,ww以及歐拉角之間的關(guān)系就可以得到一個測量點量測方程，將所有測量點的量測方程聯(lián)立為量測方程組，將該量測方程組離散化，即將k=1,2,3,…時刻的量測值表示為每個時刻狀態(tài)向量和量測噪聲的函數(shù)hk:

式中pk為k時刻的所有測量點測量值矩陣;xk為k時刻狀態(tài)向量;vk為k時刻的量測噪聲。

濾波器動力學(xué)方程(狀態(tài)方程)如下所述:

式中μ為火星的引力常量;φ為赤緯; 矩陣G為當(dāng)?shù)厮阶鴺?biāo)系到本體坐標(biāo)系的轉(zhuǎn)移矩陣;xω，yω和ωz為IMU測量得到的本體坐標(biāo)系中3個坐標(biāo)軸方向的轉(zhuǎn)動角速度;au，av和aw是進入艙質(zhì)心加速度在本體系中的投影，與IMU測量的IMU單元安裝支架處的加速度之間的關(guān)系為:

式中ax，ay和az為IMU測量得到的加速度分量;xm，ym和zm為IMU在本體坐標(biāo)系中的安裝位置;Gr為IMU坐標(biāo)系到本體坐標(biāo)系的轉(zhuǎn)移矩陣。

于是濾波器狀態(tài)方程和量測方程分別可以表示為:

式中w是分別是系統(tǒng)過程噪聲;Rk為量測噪聲協(xié)方差;Q為過程噪聲協(xié)方差;f為狀態(tài)量微分關(guān)系。

由狀態(tài)方程可以看出，系統(tǒng)過程噪聲主要指的是IMU的測量噪聲，由IMU性能決定，屬于零均值高斯白噪聲[9]。量測方程中，量測噪聲主要來自于測量單元的安裝誤差和測量時間誤差。其中測量單元的安裝誤差為15mm(3σ)，可以換算成安裝錐角λ和指針角φ的偏差，均為零均值高斯白噪聲; 測量時間誤差指的是不同測量單元測量得到同一組壓力值時的時間差，包括系統(tǒng)固有誤差以及隨機偏差。其中，固有誤差可以通過地面試驗進行矯正，隨機誤差約為±0.025s(3σ)，同樣為零均值白噪聲，符合卡爾曼濾波器的要求[10]。

濾波器初值設(shè)定為:

式中為零時刻狀態(tài)量的后驗估計，等于狀態(tài)量初值的均值E(x0)。

式中為零時刻狀態(tài)量估計誤差協(xié)方差;表示對x0的估計。

對于每個k，執(zhí)行如下遞推計算:

式中為k時刻狀態(tài)量的先驗估計;為k時刻狀態(tài)量的后驗估計。

式中=P為k時刻先驗狀態(tài)估計協(xié)方差;為k時刻后驗狀態(tài)估計協(xié)方差;為狀態(tài)轉(zhuǎn)移矩陣;為噪聲轉(zhuǎn)移矩陣。

式中Hk,i為k時刻第i次迭代測量值關(guān)于狀態(tài)量的偏微分矩陣。

式中Mk,i為k時刻第i次迭代的測量噪聲轉(zhuǎn)移矩陣。

式中Κk,i為k時刻第i次迭代的增益矩陣。

圖6 卡爾曼濾波流程Fig. 6 Rocess of the Kalman Filter

4 精度分析

上文中，基于卡爾曼濾波建立了用于估計氣動參數(shù)和大氣參數(shù)的拓展卡爾曼濾波器。由于目前無法獲得量測值以及部分由IMU測得的狀態(tài)量，為了研究卡爾曼濾波器的估計精度，選擇部分狀態(tài)量(X、Y、Z位置矢量)，利用上述濾波器對進入艙進行軌道位置估計，與彈道仿真的軌道信息對比計算位置誤差，并與最小二乘的估計誤差進行比較，其結(jié)果如圖7所示:

圖7 卡爾曼濾波誤差Fig. 7 Errors of the Kalman Filter

由上圖可以看出，相比于普通的最小二乘法狀態(tài)估計，卡爾曼濾波器的估計誤差更小，經(jīng)過約 10s時間誤差基本收斂為零，最大誤差不超過10m，具有一定的優(yōu)越性。因此，可以將卡爾曼濾波方法應(yīng)用于火星探測中進入彈道和大氣模型重建。

5 結(jié)束語

本文首先對“火星科學(xué)實驗室”的MEADS系統(tǒng)的組成以進行了介紹，以此為基礎(chǔ)詳討論了表面壓力分布模型的建立過程。通過在防熱罩上安裝的7個壓力測量單元，可以測得進入火星大氣過程中進入艙表面的壓力分布情況。若探測器安裝了用于姿態(tài)軌道控制的機載IMU設(shè)備，就可以利用表面壓力分布模型和IMU測量值，通過迭代卡爾曼濾波比較精確的計算每個時刻火星探測器的氣動參數(shù)和周圍大氣的密度和靜壓。這一方面可以重建較精確的進入彈道用于任務(wù)分析，也可以估計得到進入段的大氣密度以及靜壓隨高度的變化情況，建立火星大氣模型。本文的研究可以為將來我國開展火星探測以及建立火星大氣模型提供一定的參考依據(jù)。

References)

[1] 歐陽自遠, 李春來. 深空探測的進展與我國深空探測的發(fā)展戰(zhàn)略[C]. 2002年國際深空探測技術(shù)和應(yīng)用會議論文集，2002.OU Yangziyuan, LI Chunlai. Progress of Deep Space Exploration and Development Strategy in China[C]. The Deep Space Exploration Technology and Application Conference Proceedings, 2002.(in Chinese)

[2] 榮偉, 陳旭, 陳國良. 火星探測著陸系統(tǒng)開傘控制方法研究[J]. 航天返回與遙感, 2007, 28(4): 6–11.RONG Wei, CHEN Xu, CHEN Guoliang. The Control Method Study on the Parachute Deployment for the Mars Exploration Landing System[J]. Spacecraft Recovery & Remote Sensing, 2007, 28(4): 6–11. (in Chinese)

[3] Christopher D K, Roger E B. Mars Entry Atmospheric Data System Modeling and Algorithm Development[R].NASA-0025419, 2009.

[4] Gazarik M J, Wright M J. Overview of the MEDLI Project[C]. IEEE 2008 Aerospace Conference, 2008.

[5] Pruett C D, Wolf H, Heck M L. Innovative Air Data System for the Space Shuttle Orbiter[J]. Journal of Spacecraft and Rockets,1983, 20(1): 61–69.

[6] Siemers P M, Henry M W, Eades J B. Shuttle Entry Air Data System(SEADS)-advanced Air Data System Results: Air Data Across the Entry Speed Range[R]. NASA CP-3248, 1995.

[7] Kasich D C, Cheng P Y. Flush Port/Inertially Blended Air Data Estimator[C]. AIAA-910670, 1991.

[8] John A C, Amanda M V, Robert D B. Tatistical Reconstruction of Mars Entry, Descent, and Landing Trajectories of Atmospheric Profiles[C]. AIAA-6192, 2007.

[9] Karlgaard C D, Tartabini P V, Blanchard R C. Hyper-X Post-flight Trajectory Reconstruction[J]. Journal of Spacecraft and Rockets, 2006, 43(1): 105–115.

[10] Dan S. Optimal State Estimation[M]. Cleveland State University: John Wiley & Sons, Inc, 2006.