王文強(qiáng)

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

1 引言

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

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

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

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

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

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

圖1 MEADS測(cè)量單元安裝位置Fig.1 MEADS port geometry

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

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

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

3.1 表面壓力模型

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

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

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

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

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

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

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

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

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

圖4 測(cè)壓管安裝位置Fig.4 Pressure port geometry

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

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

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

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

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

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

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

3.2 卡爾曼濾波

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

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

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

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

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

式中Ψ,,ΘΦ分別為本體坐標(biāo)系相對(duì)于當(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)系就可以得到一個(gè)測(cè)量點(diǎn)量測(cè)方程，將所有測(cè)量點(diǎn)的量測(cè)方程聯(lián)立為量測(cè)方程組，將該量測(cè)方程組離散化，即將k=1,2,3,…時(shí)刻的量測(cè)值表示為每個(gè)時(shí)刻狀態(tài)向量和量測(cè)噪聲的函數(shù)hk:

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

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

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

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

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

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

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

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

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

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

對(duì)于每個(gè)k，執(zhí)行如下遞推計(jì)算:

式中為k時(shí)刻狀態(tài)量的先驗(yàn)估計(jì);為k時(shí)刻狀態(tài)量的后驗(yàn)估計(jì)。

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

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

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

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

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

4 精度分析

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

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

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

5 結(jié)束語(yǔ)

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

References)

[1] 歐陽(yáng)自遠(yuǎn), 李春來(lái). 深空探測(cè)的進(jìn)展與我國(guó)深空探測(cè)的發(fā)展戰(zhàn)略[C]. 2002年國(guó)際深空探測(cè)技術(shù)和應(yīng)用會(huì)議論文集，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] 榮偉, 陳旭, 陳國(guó)良. 火星探測(cè)著陸系統(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.