魏二虎,楊洪洲,2,張帥,劉經(jīng)南,易慧

(1.武漢大學(xué)測(cè)繪學(xué)院,武漢430079;2.卡爾加里大學(xué)測(cè)繪工程系,加拿大卡爾加里T2N 1N4; 3.武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,武漢430079;4.北京測(cè)繪科學(xué)研究院,北京100830)

脈沖星非實(shí)時(shí)平差的火星探測(cè)自主導(dǎo)航模型

魏二虎1,楊洪洲1,2,張帥1,劉經(jīng)南3,易慧4

(1.武漢大學(xué)測(cè)繪學(xué)院,武漢430079;2.卡爾加里大學(xué)測(cè)繪工程系,加拿大卡爾加里T2N 1N4; 3.武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,武漢430079;4.北京測(cè)繪科學(xué)研究院,北京100830)

針對(duì)基于X射線脈沖星觀測(cè)的火星探測(cè)器自主導(dǎo)航,研究了幾種不同的實(shí)時(shí)自適應(yīng)方法,包括:擴(kuò)展卡爾曼濾波(EKF),自適應(yīng)擴(kuò)展卡爾曼濾波(AEKF)和魯棒自適應(yīng)擴(kuò)展濾波(RAEKF)。首先根據(jù)脈沖星導(dǎo)航原理,模擬了觀測(cè)值:脈沖到達(dá)時(shí)刻;接著,分別利用擴(kuò)展卡爾曼濾波,自適應(yīng)擴(kuò)展卡爾曼濾波和魯棒自適應(yīng)擴(kuò)展濾波方法估算出探測(cè)器的位置和速度;最后,上述幾種濾波軌道與STK模擬的標(biāo)稱軌道較差,然后比較它們的濾波精度發(fā)現(xiàn):AEKF和RAEKF的精度相對(duì)較高,AEKF的三個(gè)軸向?yàn)V波位置精度達(dá)到:X軸優(yōu)于100 m、Y和Z軸優(yōu)于30 m,優(yōu)于VLBI技術(shù)的km量級(jí),和Doppler技術(shù)的精度相當(dāng)。

X射線脈沖星;自主導(dǎo)航;自適應(yīng)擴(kuò)展卡爾曼濾波;魯棒擴(kuò)展卡爾曼濾波

0 引言

火星探測(cè)是繼月球探測(cè)之后又一個(gè)深空探測(cè)的熱門領(lǐng)域,探測(cè)器的導(dǎo)航和定位決定著火星探測(cè)活動(dòng)的成敗與否。俄羅斯的“福布斯號(hào)”火星探測(cè)器的發(fā)射失敗再一次說(shuō)明了自主導(dǎo)航在火星探測(cè)活動(dòng)中的重要性。應(yīng)用于月球探測(cè)的導(dǎo)航技術(shù)如:VLBI和多普勒技術(shù)需要地面站的跟蹤、觀測(cè)和通信,它們不能滿足火星探測(cè)的精度要求[1-2]。此外,火星探測(cè)器與地球的遙遠(yuǎn)距離也制約著這些技術(shù)的應(yīng)用。因此,利用自主導(dǎo)航技術(shù)的火星探測(cè)器的發(fā)展可以彌補(bǔ)傳統(tǒng)導(dǎo)航定位技術(shù)的缺陷,降低探測(cè)器對(duì)地面跟蹤站的依賴,從而最大限度地確?；鹦翘綔y(cè)活動(dòng)的成功實(shí)施?；赬射線脈沖星的自主導(dǎo)航是一項(xiàng)新的極具潛力的自主導(dǎo)航技術(shù),給深空探測(cè)器實(shí)現(xiàn)自主導(dǎo)航帶來(lái)了希望。

本文正是基于這一背景,對(duì)基于X射線脈沖星的火星探測(cè)器自主導(dǎo)航技術(shù)進(jìn)行了研究。首先,介紹了相關(guān)的坐標(biāo)系統(tǒng)及其之間的轉(zhuǎn)換關(guān)系;其次,選擇了具有高觀測(cè)精度的脈沖星用于自主導(dǎo)航實(shí)驗(yàn);最后,結(jié)合火星探測(cè)器的動(dòng)力學(xué)模型和脈沖星觀測(cè)模型設(shè)計(jì)了一款自適應(yīng)擴(kuò)展卡爾曼濾波器,并開展了相關(guān)的數(shù)值實(shí)驗(yàn)和結(jié)果分析。

1 相關(guān)坐標(biāo)系統(tǒng)

在基于X射線脈沖星的火星探測(cè)器自主導(dǎo)航中涉及的坐標(biāo)系包括:太陽(yáng)系質(zhì)心坐標(biāo)系(SSB),火星赤道慣性坐標(biāo)系(SMCI)和火星中心慣性坐標(biāo)系(SECI)。文獻(xiàn)[3-4]對(duì)上述幾種坐標(biāo)系的定義做了相關(guān)介紹。這幾種坐標(biāo)系之間的轉(zhuǎn)換關(guān)系可以用下面的公式表示

式中:RMCI2ECI是由SMCI到SECI的轉(zhuǎn)換矩陣;RSB是探測(cè)器在SSB中的位置向量;RMCI是探測(cè)器在SMCI中的位置向量;RM是火星質(zhì)心在SSB中的位置向量,可以從JPL的DE421星歷表中查詢得到;F、I和E是由SMCI向SECI轉(zhuǎn)換的歐拉角,也可以由DE 421獲得。

因?yàn)榛鹦翘綔y(cè)器動(dòng)力學(xué)模型是在SMCI中描述的,而脈沖觀測(cè)模型是在SECI中建立的。需要將脈沖星時(shí)間預(yù)報(bào)模型從太陽(yáng)系質(zhì)心轉(zhuǎn)換到火星質(zhì)心,轉(zhuǎn)換精度達(dá)到0.1μs5。

2 自適應(yīng)擴(kuò)展卡爾曼濾波

2.1 軌道動(dòng)力學(xué)模型

環(huán)火探測(cè)器的軌道動(dòng)力學(xué)模型是非線性的,在慣性坐標(biāo)系SMCI中可以由下述公式來(lái)描述。

式中:r=(x,y,z)是探測(cè)器的位置向量;μM代表火星引力常數(shù);W代表狀態(tài)噪聲矩陣,是方差為Q的白噪聲;ar表示模型攝動(dòng)加速度。研究表明火星非球形攝動(dòng)的J2項(xiàng)是主要的攝動(dòng)因素[5],在軌道動(dòng)力學(xué)模型中僅考慮J2項(xiàng)攝動(dòng)因素,由此可得

因此探測(cè)器的狀態(tài)向量可以定義為

式中:v是火星探測(cè)器的速度向量,火星探測(cè)器的動(dòng)力學(xué)模型可以表述為

已知某一歷元探測(cè)器的初始狀態(tài)值,可以利用式(6)通過(guò)數(shù)值積分獲得下一歷元的狀態(tài)值。

式(6)描述的是一個(gè)非線性系統(tǒng),而傳統(tǒng)的卡爾曼濾波主要針對(duì)線性系統(tǒng)問(wèn)題。擴(kuò)展卡爾曼濾波卻可以應(yīng)用于非線性系統(tǒng)的狀態(tài)估計(jì)。

經(jīng)過(guò)線性化后,環(huán)火探測(cè)器的系統(tǒng)擾動(dòng)方程可表示為[6]

經(jīng)離散化,式(7)可寫作

式中:Φk|k-1是系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣。

2.2 狀態(tài)轉(zhuǎn)移矩陣的計(jì)算

文章采用Markley方法計(jì)算狀態(tài)轉(zhuǎn)移矩陣。它的計(jì)算需要利用tk和tk-1的探測(cè)器狀態(tài)值,μM、J2、Δt和火星半徑等。該方法在很短的時(shí)間間隔Δt內(nèi),用泰勒級(jí)數(shù)展開近似計(jì)算狀態(tài)轉(zhuǎn)移矩陣。

根據(jù)Markley方法的原理,狀態(tài)轉(zhuǎn)移矩陣可以近似表示為式(9),其中:Φ(t0,t0)=I。

由于加速度對(duì)位置的偏導(dǎo)數(shù)矩陣G是狀態(tài)值的函數(shù),可以通過(guò)數(shù)值積分獲得[78],因此,狀態(tài)轉(zhuǎn)移矩陣是可以通過(guò)積分計(jì)算的。如果僅考慮中心天體的J2項(xiàng)攝動(dòng)加速度[9]影響,矩陣G可表示為

2.3 測(cè)量模型

用于卡爾曼濾波算法的測(cè)量模型可以表達(dá)為

式中:Z表示脈沖觀測(cè)值;v為白噪聲向量,代表隨機(jī)測(cè)量誤差,具有下列統(tǒng)計(jì)特性:

式中:δkj是克羅內(nèi)克符號(hào),當(dāng)時(shí)其值為1,否則為0;v可以通過(guò)分析脈沖星TOA觀測(cè)精度得到,R表示觀測(cè)噪聲的協(xié)方差矩陣。

Hk是觀測(cè)方程矩陣,根據(jù)基于X射線脈沖星的環(huán)火探測(cè)器自主導(dǎo)航原理可知

式中:nT表示在SECI的脈沖星方向矢量,因?yàn)槊}沖星與火星距離遙遠(yuǎn),所以可以認(rèn)為該方向矢量與在SSB相同。

3 數(shù)值實(shí)驗(yàn)

為了驗(yàn)證基于X射線脈沖星的自主導(dǎo)航技術(shù)的有效性,需要開展有效的數(shù)值實(shí)驗(yàn),文章設(shè)計(jì)了幾款擴(kuò)展卡爾曼濾波器用于數(shù)值計(jì)算。

3.1 脈沖星測(cè)距精度

在利用脈沖星導(dǎo)航過(guò)程中,需要選擇最合適的脈沖星。除了要考慮脈沖星的位置、基本屬性等,還要顧及脈沖星測(cè)量精度。脈沖到達(dá)時(shí)刻的測(cè)量精度可以用下述公式描述[8]

式中:σr表示利用脈沖星測(cè)距精度;σTOA代表脈沖到達(dá)時(shí)刻TOA的測(cè)量精度;SNR是脈沖星脈沖星號(hào)的信噪比;BX表示X射線背景輻射流量;FX為X射線光子流量;A為探測(cè)器有效探測(cè)面積;T是觀測(cè)時(shí)刻;pf為脈沖流量比;d為脈沖寬度比;W50%表示脈沖輪廓寬度的50%;c是真空中光速值。

根據(jù)上述性質(zhì),實(shí)驗(yàn)選取了7顆脈沖星列于表1中。它們的測(cè)距精度信息也列于表中[1011]。

表1 優(yōu)選脈沖星的測(cè)距精度Table 1 Ranging accuracy with observation time

3.2 觀測(cè)值仿真

仿真實(shí)驗(yàn)中,環(huán)火探測(cè)器的標(biāo)稱軌道是利用Satellite Tool Kit(STK)軟件仿真得到的,歷元間隔為500 s。仿真軌道的初始軌道參數(shù)[12]見表2。

表2 火星探測(cè)器初始軌道參數(shù)Table 2 Mars parameters For observations simulation

在觀測(cè)值模擬過(guò)程中,需要顧及探測(cè)器上的原子鐘誤差和將時(shí)間相位模型由太陽(yáng)系質(zhì)心轉(zhuǎn)移到火星質(zhì)心而引起的轉(zhuǎn)換誤差。實(shí)驗(yàn)中,在觀測(cè)值中分別加入了10-9s的原子鐘誤差和10-7s的模型轉(zhuǎn)換誤差[8]。此外,還在觀測(cè)數(shù)據(jù)中隨機(jī)加入了表1中的測(cè)距誤差,使得模擬觀測(cè)值更加趨近于真實(shí)情況。

3.3 濾波

數(shù)值實(shí)驗(yàn)應(yīng)用三種數(shù)值處理方法,包括EKF、AEKF和RAEKF。其中,AEKF又有兩種不同的自適應(yīng)因子選取方法。兩種自適應(yīng)因子的選取依據(jù)下述公式

在RAEKF中,等價(jià)權(quán)的計(jì)算公式為式(19)和式(20),魯棒最小方差估值的表達(dá)式為式(21)。

式中:k0和k1設(shè)置為1.5和2.5。濾波的初始參數(shù)設(shè)置如下:

3.4 結(jié)果和討論

濾波計(jì)算中的一步預(yù)測(cè)值的計(jì)算是利用RKF7(8)數(shù)值積分器計(jì)算得到的。表3是具體的殘差信息。從中可以看出,AEKF1、AEKF2和RAEKF濾波器的濾波結(jié)果明顯優(yōu)于EKF的結(jié)果;前三者的RMS值非常接近,其中又以AEKF1的精度最好,X、Y和Z三個(gè)軸向的位置精度分別達(dá)到了94.1 m、19.5 m和22.3 m。而X方向的精度較之Y和Z偏差可能與導(dǎo)航脈沖星的分布有關(guān)。

AEKF2要達(dá)到或接近AEKF1的精度有一個(gè)先決條件就是,它只能解決類似式(24)給出的觀測(cè)方程,即必須具有多余觀測(cè)值才能夠高精度地解算狀態(tài)參數(shù)值。因此AEKF1濾波器具有更好的實(shí)用性。

RAEKF與AEKF2類似,都要求多余觀測(cè)。但是RAEKF的濾波精度略低于前兩者,其原因可能是閾值設(shè)置而導(dǎo)致的自適應(yīng)因子不恰當(dāng)。所以需要繼續(xù)改善自適應(yīng)因子的選擇方法,以便使得濾波精度進(jìn)一步提高。

表3 EKF、AEKF1、AEKF2、RAEKF濾波結(jié)果比較Table 3 Comparison of the results from EKF、AEKF1、AEKF2、RAEKF m

從表3可以看出,3個(gè)軸向的速度誤差基本上都分布在-0.1～1 m/s之間,其原因可能與式(21)中觀測(cè)方程矩陣Hk的構(gòu)造有關(guān),Hk的最后三列的值均為0。在數(shù)值計(jì)算后,誤差被吸收進(jìn)位置參數(shù)中,從而確保了速度參數(shù)的高精度。在接下來(lái)的試驗(yàn)中,在初始狀態(tài)值中分別加入了1 km的位置誤差和0.5 m/s的速度誤差。實(shí)驗(yàn)發(fā)現(xiàn),上述4種濾波器具有類似的結(jié)果,其中尤以AEKF1的精度最好。圖1和圖2是AEKF1的濾波值與標(biāo)稱值的偏差示意圖。

圖1 AEKF1濾波值與標(biāo)稱狀態(tài)的位置偏差Fig.1 Position bias between AEKF1 orbit and nominal orbit

可以看到,AEKF1在100歷元后具有較高的收斂精度。因此,除了具有較高的精度外,AEKF1還能夠抵抗初始誤差的影響,獲得精度優(yōu)于1 km的位置濾波值和1 m/s的速度濾波值。

圖2 AEKF1濾波值與標(biāo)稱狀態(tài)的速度偏差Fig.2 Velocity bias between AEKF1 orbit and nominal orbit

4 結(jié) 論

上面的理論分析和數(shù)值仿真實(shí)驗(yàn)表明:基于X射線脈沖星的自主導(dǎo)航技術(shù)能夠很好地勝任環(huán)火探測(cè)器的自主導(dǎo)航任務(wù)。實(shí)驗(yàn)數(shù)據(jù)顯示,自主導(dǎo)航系統(tǒng)狀態(tài)估計(jì)值在三個(gè)軸向X、Y和Z的位置精度分別可以達(dá)到:94.1 m、19.5 m和22.3 m。但是在實(shí)驗(yàn)中,若要獲得上述精度的濾波值,需要觀測(cè)7顆具有較高測(cè)距精度的脈沖星,這或許會(huì)限制系統(tǒng)的適用性。不過(guò)未來(lái)可以利用多頻段天線或者多臺(tái)獨(dú)立的天線同時(shí)觀測(cè)來(lái)解決這一問(wèn)題。

在上面的實(shí)驗(yàn)中,AEKF1濾波器表現(xiàn)良好而且適用范圍較廣、受條件限制較小,因此該方法可以應(yīng)用于火星探測(cè)器的自主導(dǎo)航中,從而彌補(bǔ)地面跟蹤系統(tǒng)(包括VLBI、ΔDOR和ΔDOD)和測(cè)距技術(shù)的不足,尤其是當(dāng)?shù)孛嬲靖櫜坏教綔y(cè)器的時(shí)候。基于脈沖星的自主導(dǎo)航技術(shù)與地面跟蹤系統(tǒng)的優(yōu)勢(shì)互補(bǔ)為未來(lái)的深空探測(cè)器導(dǎo)航指明了方向。

在本文的實(shí)驗(yàn)中,火星大氣阻力、太陽(yáng)光壓和火星非球形引力的高階項(xiàng)等攝動(dòng)力因素并沒(méi)有考慮進(jìn)探測(cè)器的動(dòng)力學(xué)模型中。若建立顧及上述諸多攝動(dòng)因素的動(dòng)力學(xué)模型,可以進(jìn)一步提高探測(cè)器的導(dǎo)航定位精度。因此,如何完善火星探測(cè)器的軌道動(dòng)力學(xué)模型需要進(jìn)一步的研究。

[1]Cao J F,Huang Y,Hu X,et al.Mars express tracking and orbit determination trails with Chinese VLBI network[J]. Chinese Sci Bull,2010(55):3654-3660.

[2]李金玲,劉鵬,馬茂莉,等.中國(guó)VLBI網(wǎng)火星快車衛(wèi)星跟蹤資料的定位歸算[J].宇航學(xué)報(bào),2010,31(7):1718-1723.[Li J L,Liu P,Ma M L,et al.Positioning reduction of the mars express tracking data by the Chinese VLBI network[J].J. Astronautics,2010,31(7):1718-1723.]

[3]Wei J,Xu S.Analysis on the influence of natural convection during melting of phase change materia[J].Chinese Space Sci.Tech,2009,29(2):65-70.

[4]Brukhart P D.MSL update to Mars coordinate frame definitions[EB/OL].[2014-07-29].http://www. coursehero.comfile2541877 HW4.2006-08-15:1-20.

[5]劉林,湯靖師.火星軌道器的運(yùn)動(dòng)變化特征[J].宇航學(xué)報(bào), 2008,3(2):462-466.[Liu L,Tang J S.Orbit variation characteristics of the Mars'orbiters[J].J.Astronautics, 2008,3(2):462-466.]

[6]帥平,李明,陳紹龍,等.X射線脈沖星導(dǎo)航原理與方法[M].北京:中國(guó)宇航出版社,2009:466-470.[Shui P,Li M,Chen S L,et al.Principles and methods of X-ray pulsars navigation [M].Beijing:China Astronautic Publishing House,2009: 466-470.]

[7]Markley J.Approximate cartesian state transition matrix[J]. Journal of the Astronautical Science,1986,34(2):161-169.

[8]Chiaradia A P M,Kuga H K.Comparision between two methods to calculate the transition matrix of orbit motion[J]. Mathematical Problems in Engineering,2012,2012 (768973):1-13.

[9]Moraes R V,Silva A A,Kuga H K.Simple orbit determination using GPS based on a least-squares algorithm employing sequential given rotations[J].Mathematical Problems in Engineering,2007(49781):1-8.

[10]Sheikh S I.The use of variable celestial X-ray sources for spacecraft for navigation[D].Maryland:Department of Aerospace Enginerring,University of Maryland,2005.

[11]Kargaltsev O,et al.X-ray pulsars database[EB/OL]. 2012.http://www2.astro.psu.edu/users/green/ psrdatabase.

[12]Lemoine F G,Smith D E,Rowlands D D,et al.An improved solution of the gravity field of Mars from Mars global surveyer[J].J Geophys.Res,2001,106(E10): 23359-23376.

[13]Wei E,Jin S G,Zhang A,et al,Autonomous navigation of Mars probe using X-ray pulsars:modeling and results[J]. Adv.Space Res,2013,51(5):849-857.

[14]Wei E,Jin S G,Yang H,et al.Simulation and results on realtime positioning of Chang'E-3 rover with the same-beam VLBI observations[J].Planetary and Space Science,2013,84:20-27.

[15]Yang Y,He H,Xu G.A new adaptively robust filtering for kinematic geodetic positioning[J].Journal of Geodesy,2001, 75(2):109-116.

通信地址:湖北省武漢市珞喻路129號(hào)武漢大學(xué)測(cè)繪學(xué)院空間定位與導(dǎo)航工程研究所(430079)

電話:(027)68758505

E-mail:ehwei@sgg.whu.edu.cn

[責(zé)任編輯:高莎]

Modeling on Autonomous Navigation of Mars Probe with Pulsars and Nonreal-Time Adjustment Methods

WEI Erhu1,YANG Hongzhou1,2,ZHANG Shuai1,LIU Jingnan3,YI Hui4
(1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China;2.Department of Geomatics Engineering, University of Calgary,Calgary T2N 1N4,Canada;3.GNSS Center,Wuhan University,Wuhan 430079,China; 4.Chinese Academy of Surveying and Mapping,Beijing 100830,China)

To investigate the autonomous navigation of probe orbiting the Mars orbit X-ray pulsars,different real-time adjustment methods are applied in this paper,including the extended Kalman filter(EKF),adaptive extended Kalman filter(AEKF)and robust adaptive extended Kalman filter(RAEKF).Firstly,the observation data,i.e.,time of arrival,is simulated according to the principles of pulsars navigation.Secondly,extended Kalman filter(EKF),adaptive extended Kalman filter(AEKF)and robust adaptive extended Kalman filter (RAEKF)are used to get the position and velocity of the probe.Thirdly,the adjustment results are compared with the orbit generated by STK,the results show AEKF and RAEKF performs better than EKF,the result of AEKF is better than the km level of VLBI and can reach the precision level of triple pass Doppler,leveling at 94.1 m in X-axis,19.5 m in Y-axis and 22.3 m in Z-axis.

X-ray pulsars;autonomous navigation;adaptive extended Kalman filter(AEKF);robust adaptive extended Kalman filter(RAEKF)

P1

:A

:2095-7777(2014)04-0298-05

10.15982/j.issn.2095-7777.2014.04.009

魏二虎(1961—),男,教授,博士生導(dǎo)師,主要研究方向?yàn)榭臻g大地測(cè)量。

2014-07-29

2014-08-09

國(guó)家自然科學(xué)基金資助項(xiàng)目(41374012)