陳培 韓錦飛 賴玉敏 孫秀聰 譚龍玉

(1 北京航空航天大學(xué)宇航學(xué)院，北京 100191) (2 上海航天控制研究所，上海 201109) (3 上?？臻g智能控制技術(shù)重點(diǎn)實(shí)驗(yàn)室，上海 201109)

航天器的自主導(dǎo)航在航天任務(wù)中有著舉足輕重的地位，對(duì)地球附近(地球中低軌道高度)和近地行星附近(近地行星中低軌道高度)的航天器導(dǎo)航都是如此：地球附近的航天器任務(wù)中，目前大多使用衛(wèi)星定位(如GPS、北斗等)，其優(yōu)點(diǎn)是定位精度高，且可以實(shí)現(xiàn)實(shí)時(shí)定位導(dǎo)航，但其也存在著電磁波信號(hào)易受到干擾和欺騙的問題[1]，尤其是軍用航天器，為保證安全性與穩(wěn)定性，其對(duì)自主導(dǎo)航的需求更為迫切；另外，近年來對(duì)于近地行星如火星等的探測任務(wù)逐漸增多，而目前GPS等衛(wèi)星導(dǎo)航在其他行星附近處于完全不可用的狀態(tài)，近地行星附近的探測器也需要一種更加有效且自主的導(dǎo)航方式。對(duì)于近地行星附近的航天器導(dǎo)航，目前存在的導(dǎo)航方式有限。常見的航天器自主導(dǎo)航方式有磁場導(dǎo)航、星光折射導(dǎo)航、X脈沖星導(dǎo)航、慣性導(dǎo)航等[2]。磁場導(dǎo)航利用地球附近每個(gè)點(diǎn)與地磁矢量存在的一一對(duì)應(yīng)關(guān)系進(jìn)行匹配導(dǎo)航[3]。其優(yōu)點(diǎn)是利用了地球本身的磁場信息，可以實(shí)現(xiàn)實(shí)時(shí)測量解算。但地磁導(dǎo)航也存在著地磁場模型精度不夠高的缺點(diǎn)，而且由于地球本身的磁場信號(hào)比較弱，測量過程中很容易受到周圍磁場干擾。而且在一些行星，如金星、火星，據(jù)目前的資料顯示，其上是沒有磁場的，因此無法利用磁場對(duì)其附近地航天器進(jìn)行導(dǎo)航。星光折射屬于天文導(dǎo)航，容易受到天體可見性，以及大氣模型的影響，而大氣模型隨著季節(jié)、維度、太陽活動(dòng)而不斷變化，因此利用星光折射法對(duì)地球附近的航天器進(jìn)行導(dǎo)航精度并不高[4]；另外，由于近地行星周圍沒有像地球一樣厚厚的大氣層，因此，星光折射無法應(yīng)用。X射線脈沖星導(dǎo)航能夠提供一種獨(dú)立的絕對(duì)時(shí)空基準(zhǔn)，解決星座長時(shí)間自主運(yùn)行的問題，為航天器自主導(dǎo)航提供了一種新的思路和可行途徑。而其缺點(diǎn)也很明顯，系統(tǒng)非常復(fù)雜，且只能實(shí)現(xiàn)相對(duì)慣性系的導(dǎo)航，相對(duì)于行星本體系的導(dǎo)航精度沒有保證[5]；慣性導(dǎo)航的優(yōu)點(diǎn)是在短時(shí)間內(nèi)的導(dǎo)航精度很高，但由于眾所周知的誤差累積問題，長時(shí)間導(dǎo)航時(shí)會(huì)產(chǎn)生很大的導(dǎo)航誤差，因此無法勝任長時(shí)間的自主導(dǎo)航任務(wù)[6]。以上導(dǎo)航方式都沒辦法對(duì)實(shí)現(xiàn)近地行星附近的導(dǎo)航，因此需要尋找一種能夠完全自主且能夠長時(shí)間在行星附近進(jìn)行導(dǎo)航的方法。

近年來，世界航天強(qiáng)國也各自制定了相應(yīng)的航天計(jì)劃：美國的“空間探索計(jì)劃”、俄羅斯的未來航天發(fā)展規(guī)劃、歐空局制定的“曙光計(jì)劃”與“歐洲阿波羅計(jì)劃”等超大規(guī)模衛(wèi)星探測計(jì)劃等[7]。與此同時(shí)，我國在深空探測領(lǐng)域也取得了顯著成就，對(duì)于月球以及火星的探索都展開了研究工作。但是，目前我國的航天探測任務(wù)主要以地面測控為主，有測控精度低、通訊時(shí)延長等缺點(diǎn)[8]。其對(duì)于地面設(shè)備的依賴使得其自主生存能力大大減弱。而且地面遙測遙控的方式無法滿足交會(huì)飛越、下降著陸等任務(wù)階段的實(shí)時(shí)性要求[9]，在此情況下，航天器的自主導(dǎo)航顯得尤為重要。本文在上述研究基礎(chǔ)上，針對(duì)目前近地行星附近的自主導(dǎo)航問題，提出了近地行星附近的引力梯度導(dǎo)航方法，并以火星為例，初步探索該方法應(yīng)用于火星附近自主導(dǎo)航的可行性。

1 引力梯度匹配定位研究

隨著大地測量技術(shù)的飛速發(fā)展，近地行星的引力場模型已經(jīng)具有很高的精度，而且重力梯度儀也有很高的測量精度，歐空局在2009年已經(jīng)將高精度的靜電重力梯度儀配置于GOCE衛(wèi)星上進(jìn)行引力梯度測量。另外，引力梯度導(dǎo)航技術(shù)日趨成熟，并且已經(jīng)成功應(yīng)用于水下潛器和低空飛行器的慣性輔助導(dǎo)航系統(tǒng)中，并成功提高了慣性導(dǎo)航的精度[10-11]。但是在航天器導(dǎo)航上的應(yīng)用卻很少有人研究。筆者已經(jīng)將全張量重力梯度應(yīng)用到繞地飛行衛(wèi)星導(dǎo)航系統(tǒng)。對(duì)于中心引力模型或只考慮J2項(xiàng)影響的模型，利用特征值分解可以直接獲得估計(jì)位置，同時(shí)采用概率統(tǒng)計(jì)的方式解決了解的符號(hào)模糊問題[12]；對(duì)于考慮高階項(xiàng)的EGM2008重力場模型，利用卡爾曼濾波可以獲得位置的較為精確估計(jì)值[13]。以上定位方法均利用GOCE衛(wèi)星實(shí)測數(shù)據(jù)進(jìn)行了可行性驗(yàn)證，取得了成功。

目前對(duì)于火星引力場的研究，如文獻(xiàn)[14]研究了目前兩種最新的火星引力場模型JGMRO-120D和GMM-3，引力異常和星歷積分兩個(gè)方面的分析結(jié)果表明，對(duì)于10 m量級(jí)的定軌精度，這兩個(gè)模型均可滿足導(dǎo)航需求。由于JGMRO-120D模型可以直接從NASA官網(wǎng)下載，模型文件中含有火星引力場模型的時(shí)變修正公式，且此修正公式相對(duì)于GMM-3模型更為簡單，使用也更加方便，因此本文采用JGMRO-120D模型，并在使用的過程中，考慮了火星引力場的時(shí)變特性，加入了時(shí)變修正項(xiàng)。通過仿真算例，實(shí)現(xiàn)利用引力梯度對(duì)航天器進(jìn)行自主導(dǎo)航，為近地行星附近的自主導(dǎo)航方式的研究開辟一條新的道路。下面首先對(duì)引力梯度的基本概念加以闡釋。引力梯度是引力加速度對(duì)位置向量的導(dǎo)數(shù)，其單位是E(E?tv?s, 1E=10-9s-2)，也是引力勢對(duì)位置向量的二階導(dǎo)數(shù)。它是一個(gè)二維矩陣，稱為引力梯度張量(gravity gradient tensor，GGT)矩陣。

(1)

式中：g是引力加速度向量；r是位置向量；(gx,gy,gz)和(x,y,z)分別是g和r在行星固連坐標(biāo)系下對(duì)應(yīng)方向上的分量。

火星周圍空間各處引力場強(qiáng)度分布不同，它是由火星本身決定的。和火星相固連的每一點(diǎn)都有唯一的GGT，即位置和引力梯度張量是一一對(duì)應(yīng)的。重力梯度儀并不能直接測量在火星固連坐標(biāo)系下的引力梯度Vij，而只能測量在重力梯度儀本體坐標(biāo)系(Gradiometer Reference Frame ，GRF)下的引力梯度。兩個(gè)坐標(biāo)系的定義分別如圖1所示：O-XYZ表示火星固連坐標(biāo)系，O為火星質(zhì)心，OZ由地心指向協(xié)議北極點(diǎn)，OX由火星質(zhì)心指向火星歷元平赤道面與地球J2000歷元平赤道面的交點(diǎn)，OY由右手空間直角坐標(biāo)系可得。而GRF則與重力梯度儀安裝角度有關(guān)，這里可認(rèn)為其定義見圖1，O′表示重力梯度儀中心，O′Z′沿著OO＇方向并遠(yuǎn)離火星質(zhì)心，O′X′表示沿著軌道切線方向并與速度方向相同，而O′Y′則垂直于軌道平面[14]。以VM和VG分別表示火星固連坐標(biāo)系和重力梯度儀本體坐標(biāo)系下的GGT，則VM可以用下面的轉(zhuǎn)換方式轉(zhuǎn)換到VG。

(2)

圖1 火星固連坐標(biāo)系(MCF)和重力梯度儀坐標(biāo)系(GRF)Fig.1 Mars-centered Mars-fixed reference frame and gravity gradiometer reference frame

由于GGT是對(duì)稱的，所以重力梯度儀只需輸出VG中的6個(gè)元素，因此在構(gòu)造觀測方程時(shí)，引入下面兩個(gè)列向量

(3)

則觀測方程可以寫成

(4)

(5)

2 引力梯度匹配定位方法

與地球引力場模型相比，火星引力場模型的確定更為復(fù)雜，其中涉及到火星引力場季節(jié)性變化。火星非球形引力場模型與地球不同，其非球形引力位中的球諧項(xiàng)系數(shù)基本都要比地球相應(yīng)值大一個(gè)量級(jí)，尤其是赤道橢率項(xiàng)，其大小接近于火星動(dòng)力學(xué)扁率項(xiàng) 。由于火星引力場的這些特性，使得基于火星引力場的定位更加困難，因此，本文在建模時(shí)加入了對(duì)火星重力場模型時(shí)變量的修正。下面介紹一種將特征值分解和最小二乘迭代結(jié)合起來進(jìn)行位置估計(jì)的方法。由于引力梯度是位置向量的非線性函數(shù)，故所要研究的問題是非線性問題。這里用最小二乘迭代的方法對(duì)位置進(jìn)行求解。在非線性最小二乘問題中，初值的選取是比較關(guān)鍵的。在本方法中，初值用特征值分解法給出。在中心引力場模型中，將引力加速度對(duì)位置向量求一階導(dǎo)數(shù)，得到中心引力場下GGT在行星固連坐標(biāo)系下的表達(dá)式

(6)

式中：μ表示火星引力常數(shù)；r表示重力梯度儀的火心距，r=|r|；對(duì)求得的GGT進(jìn)行特征值分解得到如下

(7)

Γ中對(duì)角線元素即為GGT的特征值，Φ中的3個(gè)列向量分別對(duì)應(yīng)其特征向量。從Γ中可以看出，特征值只取決于r，無論其在坐標(biāo)系中的任何方位，只要距離球心的距離相等，GGT分解出的特征值就相等；對(duì)應(yīng)的特征向量只取決于φ和λ，二者分別代表了行星上的經(jīng)度和緯度,見圖1，無論其距離球心的距離如何變化，只要它所對(duì)應(yīng)的經(jīng)緯度沒有發(fā)生變化，特征向量就不變。事實(shí)上，由球坐標(biāo)系與笛卡爾坐標(biāo)系之間的關(guān)系可知，Φ中的第三個(gè)特征向量對(duì)應(yīng)著它們特征值中r的方向，是位置向量的單位向量，表示指向。即

r=rΦ3

(8)

式中：Φ3表示特征向量矩陣Φ的第3列。

用ξ表示特征值中的最大值，η為對(duì)應(yīng)的特征向量。由式(7)可知

(9)

求得位置向量為

(10)

這里需要注意的是，特征向量存在方向上的不確定性，也就是說通過特征值分解法我們解算出的應(yīng)該是兩個(gè)關(guān)于球心對(duì)稱的位置向量。這樣在給定初值的時(shí)候就會(huì)出現(xiàn)兩個(gè)位置初值，從而造成結(jié)果的歧義性。要消除這一問題，如果知道大致的初始位置，則可以采用比較法排除差距較大的值。

用特征值分解法給出位置初值后，采用非線性最小二乘法迭代求解位置，進(jìn)一步提高位置解的精度。假設(shè)Q為重力梯度儀測量噪聲方差矩陣。用非線性最小二乘迭代方法求解的迭代方程

(11)

(12)

高階次引力場模型的Hk表達(dá)式比較復(fù)雜，這里不再詳細(xì)列出。對(duì)于初值，Hk可由式(12)計(jì)算，在迭代過程中，Hk由高階次的公式得到。位置估計(jì)的協(xié)方差矩陣可以由協(xié)方差公式得到

(13)

由于估計(jì)的協(xié)方差表示估計(jì)結(jié)果在各個(gè)方向上對(duì)真值的偏離程度，可以用其來描述最終定位的精度，同樣，我們也可以用其平方根來表示精度。這里我們用其平方根，即使用標(biāo)準(zhǔn)差來表示位置估計(jì)的精度。方法流程如圖2所示。

圖2 行星附近引力梯度導(dǎo)航方法流程圖

3 仿真分析

基于以上原理，在MATLAB軟件中，以火星為例進(jìn)行仿真驗(yàn)證。衛(wèi)星軌道高度設(shè)定為距離火星表面300 km，軌道偏心率為0.05。動(dòng)力學(xué)模型為NASA發(fā)布的最新的火星引力場JGMRO-120D模型，其最高階數(shù)為120?；鹦侵亓瞿Ｐ途哂袝r(shí)變特性，在JGMRO-120D模型文件中，給出了模型的時(shí)變修正項(xiàng)

dJ3=2.958 705 874 08×10-9×sin [w×(t-t0)]

(14)

式中：dJ3表示火星重力場模型的J3項(xiàng)修正系數(shù);w為火星自轉(zhuǎn)的平均角速度，w=191.39 ()/a；(t-t0)為距離1999年1月1日的時(shí)間長度，單位為a。在仿真時(shí)已經(jīng)加入了火星重力場的時(shí)變特性。考慮到計(jì)算速度，仿真采用前20階球諧項(xiàng)系數(shù)，假定重力梯度儀固連坐標(biāo)系與軌道坐標(biāo)系重合。設(shè)定仿真采樣頻率為0.1 Hz，即采樣間隔為10 s，一共采集6000個(gè)歷元的數(shù)據(jù)。仿真分兩次分別添加兩種強(qiáng)度不同的噪聲，一種與當(dāng)前重力梯度儀的測量噪聲相當(dāng)，每個(gè)方向都設(shè)為正態(tài)分布的高斯白噪聲，其標(biāo)準(zhǔn)差為0.1 E(100 mE)；另一種是標(biāo)準(zhǔn)差為10 mE的高斯白噪聲。

當(dāng)測量噪聲為100 mE量級(jí)時(shí)，測量噪聲協(xié)方差陣為

(15)

由特征值分解法可以得到位置解，由于實(shí)際解算時(shí)采用的是高階引力場模型，因此需要將此位置解作為初始值利用最小二乘法進(jìn)行迭代求解。將估計(jì)值與仿真的真值比較，可以得到在火星固連坐標(biāo)系下位置估算的3個(gè)方向誤差。對(duì)于估計(jì)方差，采用單歷元的測量數(shù)據(jù)進(jìn)行分析，每個(gè)歷元有6個(gè)測量數(shù)據(jù)，迭代得到對(duì)應(yīng)的方差和標(biāo)準(zhǔn)差，如圖3所示。

圖3 噪聲100 mE時(shí)位置估算偏差及相應(yīng)的標(biāo)準(zhǔn)差Fig.3 Position estimation error and corresponding standard deviation in noise 100mE

仿真結(jié)果表明，在300 km高度處，測量噪聲在100 mE量級(jí)時(shí)，x方向估計(jì)誤差約為154.3 m，y方向誤差約為170.3 m，z方向估計(jì)誤差約為139.3 m，三維誤差約為268.9 m。

隨著重力梯度儀測量精度的提高，測量噪聲越來越小，未來有可能將測量噪聲降低到10 mE以下。這里仿真得到測量噪聲為10 mE時(shí)的結(jié)果：當(dāng)測量噪聲設(shè)置為10 mE時(shí)，用同樣的方法得到火星固連坐標(biāo)系下仿真定位結(jié)果，如圖4所示。

圖4 噪聲10 mE時(shí)位置估算偏差及相應(yīng)的標(biāo)準(zhǔn)差Fig.4 Position estimation error and corresponding standard deviation in noise 10mE

仿真結(jié)果表明，在測量噪聲為10 mE量級(jí)時(shí)，x方向估計(jì)誤差約為15.5 m，y方向誤差約為16.7 m，z方向估計(jì)誤差約為14.1 m，三維估計(jì)誤差約為26.8 m，不同測量噪聲水平下解算結(jié)果如表1所示。

表1 不同測量噪聲水平下系統(tǒng)定位誤差

從表1中分析可知，當(dāng)重力梯度儀測量噪聲水平在10 mE時(shí)，由此方法可以得到30 m以內(nèi)的定位結(jié)果，當(dāng)重力梯度儀測量噪聲水平在100 mE時(shí)，定位誤差也能控制在300 m以內(nèi)。由于目前重力梯度儀的測量水平越來越高，美國馬里蘭大學(xué)的Paik團(tuán)隊(duì)甚至研制出了測量水平在1 mE量級(jí)的超導(dǎo)重力梯度儀[15]，由此可以看到，此種方法在已知高精度引力場的行星附近進(jìn)行導(dǎo)航是完全可行的。

4 結(jié)束語

在近地行星附近，引力梯度導(dǎo)航是一種實(shí)用性非常強(qiáng)的導(dǎo)航方法。本文通過理論分析了引力梯度進(jìn)行近地行星導(dǎo)航的可行性，給出了一種基于引力梯度的特征值分解導(dǎo)航算法，最后以火星為例進(jìn)行了數(shù)值仿真。結(jié)果表明，在當(dāng)前100 mE測量噪聲水平，定位誤差達(dá)到268 m左右，具有較高的精度；如果未來10 mE水平的重力梯度儀能夠得到應(yīng)用，則可以得到更加精確的導(dǎo)航結(jié)果。與其它方法相比，本方法是首次利用引力梯度對(duì)火星附近的航天器進(jìn)行導(dǎo)航驗(yàn)證。方法整體思路簡單，易于實(shí)現(xiàn)，且不用依賴任何地基設(shè)備即可對(duì)行星附近航天器進(jìn)行長時(shí)間、全天候的實(shí)時(shí)自主導(dǎo)航，大大增加了此種導(dǎo)航方法在行星附近航天器上的應(yīng)用價(jià)值。另外，文中采取的模型僅為較低階的火星引力場模型，利用仿真數(shù)據(jù)驗(yàn)證方法的可行性與可能達(dá)到的定位精度。在實(shí)際航天任務(wù)中，可采用更高階的火星引力場模型，以進(jìn)一步提高定位精度。雖然本文以火星為例，但是此種方法適用于地球、類地行星及其它所有引力場接近有能力場且比較穩(wěn)定的各類較大質(zhì)量星體的中低軌道衛(wèi)星的自主導(dǎo)航，應(yīng)用前景十分廣闊。

參考文獻(xiàn)(References)

[1] 趙金磊. GPS定位及欺騙干擾技術(shù)[D]. 西安: 西安電子科技大學(xué), 2014

Zhao Jinlei. GPS positioning and deception jamming technology[D]. Xi’an: Xi’an Electronic and Science University, 2014 (in Chinese)

[2] 姜竹青. 自主導(dǎo)航中濾波算法的研究及應(yīng)用[D]. 北京：北京郵電大學(xué)，2014

Jiang Zuqing. Research and application of filtering algorithm in autonomous navigation[D]. Beijing: Beijing University of Posts and Telecommunications, 2014

[3] 寇義民. 地磁導(dǎo)航關(guān)鍵技術(shù)研究[D]. 哈爾濱：哈爾濱工業(yè)大學(xué), 2010

Kou Yimin. Research on the key technology of geomagnetic navigation[D]. Harbin: Harbin Institute of Technology, 2010. (in Chinese)

[4] 寧曉琳, 王龍華, 白鑫貝, 等. 一種星光折射衛(wèi)星自主導(dǎo)航系統(tǒng)方案設(shè)計(jì)[J]. 宇航學(xué)報(bào), 2012, 33(11): 1601-1610

Ning Xiaolin, Wang Longhua, Bai Xinbei, et al. A scheme design of satellite autonomous navigation system based on stellar refraction[J]. Journal of Astronautics, 2012, 33(11): 1601-1610 (in Chinese)

[5] 孫景榮. X射線脈沖星導(dǎo)航及其增強(qiáng)方法研究[D]. 西安：西安電子科技大學(xué), 2014

Sun Jingrong. X ray pulsar navigation and its enhancement method [D]. Xi’an: Xi’an Electronic and Science University, 2014 (in Chinese)

[6] 吳俊偉, 曾啟明, 聶莉娟. 慣性導(dǎo)航系統(tǒng)的誤差估計(jì)[J]. 中國慣性技術(shù)學(xué)報(bào), 2002, 10(6):1-5

Wu Junwei， Tsang K M, Nie Lijuan. Estimation of the INS’s errors[J]. Journal of Chinese Inertial Technology, 2002, 10(6):1-5 (in Chinese)

[7] 韓鴻碩. 21世紀(jì)國外深空探測發(fā)展計(jì)劃及進(jìn)展[J]. 航天器工程, 2008, 17(3): 1-22

Han Hongsuo. 21stcentury foreign deep space exploration development plans and their progress [J]. Spacecraft Engineering, 2008, 17(3): 1-22 (in Chinese)

[8] 王大軼, 黃翔宇. 深空探測自主導(dǎo)航與控制技術(shù)綜述[J]. 空間控制技術(shù)與應(yīng)用, 2009, 35(03):6-12

Wang Dayi, Huang Xiangyu. Review of the autonomous navigation and control technology of deep space exploration[J]. Space Control and Technology and Application, 2009, 35(3): 6-12 (in Chinese)

[9] 吳偉仁, 于登云. 深空探測發(fā)展與未來關(guān)鍵技術(shù)[J]. 深空探測學(xué)報(bào), 2014, 1(1): 5-17

Wu Weiren, Yu Dengyun. Development of deep space exploration and its future key technologies[J]. Journal of Deep Space Exploration, 2014, 1(1): 5-17 (in Chinese)

[10] 李蘭玉. 基于全張量重力梯度的水下導(dǎo)航技術(shù)研究[D]. 武漢：武漢科技大學(xué), 2014

Li Lanyu. Underwater navigation technology based on full tensor gravity gradient[D]. Wuhan:Wuhan University of Science and Technology, 2014 (in Chinese)

[11] 王文晶. 基于重力和環(huán)境特征的水下導(dǎo)航定位方法研究[D]. 哈爾濱工程大學(xué), 2009

Wang Wenjing. Underwater navigation methods based on gravity and environmental features[D]. Harbin: Harbin Engineering University, 2009 (in Chinese)

[12] Chen P, Sun X, Han C. Gravity gradient tensor eigende composition for spacecraft positioning[J]. Journal of Guidance Control & Dynamics, 2015, 38: 2200-2206

[13] Sun X, Chen P, Macabiau C, et al. Autonomous orbit determination via Kalman filtering of gravity gradients[J]. IEEE Transactions on Aerospace and Electronic Systems, 2016, 52(5)

[14] 曹建峰, 劉磊, 黃勇, 等. 火星指向模型與重力場模型的發(fā)展回顧與使用[J]. 天文學(xué)進(jìn)展, 2017, 35(1):127-139

Cao Jianfeng, Liu Lei, Huang Yong, et al. Review and the utilization of martian orientation model and gravity field model[J]. Progress in Astronomy, 2017, 35(1):127-139 (in Chinese)

[15] Moody M V, Paik H J, Canavan E R. Three-axis superconducting gravity gradiometer for sensitive gravity experiments[J]. Review of Scientific Instruments, 2002, 73(11): 3957-3974 (in Chinese)