羅 燕，秦 同，喬 棟

（1. 北京理工大學宇航學院，北京100081；2. 清華大學精儀系，北京100084）

1 引 言

小天體是太陽系中較古老的一類天體，探測小天體，獲取其物理參數(shù)，是研究太陽系形成與演化、生命起源與進化，以及抵御外來天體撞擊的重要技術途徑，是人類探索和開發(fā)宇宙資源的重要技術手段。同時，也為試驗驗證空間技術，尤其是深空探測技術提供了平臺。目前，小天體物理參數(shù)的獲取方式主要包括近距離光學成像、多譜段光譜測量和表面采樣探測。此類方式可以感知小天體的尺寸、形狀、自旋等物理參數(shù)，以及小天體表面地形、地表物質組成等。但光學觀測的方式精度較低，且無法獲得小天體內(nèi)部質量與密度分布等信息。小天體的引力場以及自旋參數(shù)是小天體探測的基本內(nèi)容之一。感知高精度的小天體引力場以及自旋，特別是不規(guī)則形狀小天體引力場與自旋，有助于探索小天體的內(nèi)部結構、質量與密度分布等信息，揭示小天體起源，獲得新的科學發(fā)現(xiàn)。

小天體引力場以及自旋信息的獲取通常有兩種方式。其一，從傳統(tǒng)的光學觀測中提取光度或光學特征點等有效信息，以此為觀測量反演自旋角速度以及外部形貌［1-2］，在均勻密度的假設下，獲得質量分布，從而建立引力場模型［3-5］。由于這類方法的信息源來自小天體表面光學觀測，因此獲得的自旋角速度信息以及引力場模型較為粗糙，且無法體現(xiàn)小天體內(nèi)部結構與質量信息。另一種方法是通過探測器軌道數(shù)據(jù)進行反演。未來的探測器需具有全面的感知能力，包括環(huán)境感知及自身狀態(tài)感知［6-7］。由于探測器的運動由其軌道動力學環(huán)境，即探測器所處引力場決定，且小天體自旋會影響引力場的空間分布，因此借助精確的軌道數(shù)據(jù)，可對小天體引力場與自旋進行精確感知［8］。NEAR 探測器在對Eros 433 小行星的近距離飛越過程中，通過地面的無線電測量與星載光學測量數(shù)據(jù)不僅實現(xiàn)了高精度定軌［9］，同時還對Eros 433 的引力場進行了估計，并在美國國家航空航天局官網(wǎng)公布了引力場球諧系數(shù)模型的前8 階球諧系數(shù)，是通過軌道數(shù)據(jù)感知小天體引力場的先例。

高精度測量與動力學估計是感知小天體自旋角速度與引力場模型的前提。現(xiàn)階段，對小天體探測器的測量主要為地面無線電測量，包括測距、測速、甚長基線干涉測量，以及星載光學相機測視線與激光測距儀測距［10-12］。綜合各類測量信息，結合較粗糙的先驗引力場模型，可以估計探測器在小天體坐標系的軌道。然而，地面測控定軌自主性差、精度隨距離增加而降低、且測量幾何構型變化弱，無法實現(xiàn)小天體坐標系下的精密定軌。此外，小天體自旋與引力場感知需要長時間連續(xù)測量估計，依靠地面測控成本高昂。光學測量雖自主性強，但受光學導航原理（圖像特征提取與匹配）的限制，參數(shù)估計精度較低。因此，傳統(tǒng)的測量定軌方式在定軌精度上只能滿足工程需求，無法滿足高精度引力場與自旋參數(shù)感知的需求。

本文針對無人系統(tǒng)在小天體科學探測方面的應用，研究基于軌道動力學特性的小天體引力場與自旋角速度感知問題，即如何通過小天體探測器之間的器間測距信息對小天體引力場模型參數(shù)與自旋角速度進行修正。該問題本質上為參數(shù)估計問題，可用的測量信息為器間測距，需估計的參數(shù)為引力場模型參數(shù)、自旋角速度以及探測器軌道。在二體動力學問題中，由于引力場的對稱性，器間相對測量僅能確定探測器軌道的大小、形狀與尺寸，但星座的旋轉方位信息，即各個軌道升交點赤經(jīng)、軌道傾角以及近地點幅角均不可估計［13-14］。原因在于，星座整體旋轉時，器間的測距、測速信息不會發(fā)生變化，導致器間測量導航系統(tǒng)不完全可觀。但在小天體的不規(guī)則引力場中，引力場各向差異顯著，不再呈現(xiàn)對稱性，因此星座旋轉后，各軌道的大小、形狀尺寸也會隨之改變，從而引起器間測量的變化。如此一來，器間測量對軌道各個元素均有一定的敏感性，整個星座定軌系統(tǒng)完全可觀。因此，小天體引力場的不規(guī)則性雖然增加了小天體附近軌道設計難度，但也由于其復雜性，改善了參數(shù)感知系統(tǒng)的可觀性，是參數(shù)感知的有利因素。

小天體的不規(guī)則引力場是實現(xiàn)參數(shù)感知的關鍵因素，同時，由于先驗引力場模型以及自旋角速度存在誤差，在進行定軌估計時，除了各軌道狀態(tài)，還需要將引力場模型系數(shù)以及自旋角速度作為待估狀態(tài)變量進行估計。因此，基于動力學的小天體參數(shù)感知實質上為動力學參數(shù)估計。以小天體球諧系數(shù)模型為例，通過三維建模與均勻密度假設得到的引力場模型中各階球諧系數(shù)誤差較大，將各階球諧系數(shù)以及自旋角速度作為狀態(tài)變量，擴充至動力學估計系統(tǒng)的狀態(tài)向量中，依靠器間觀測量，通過非線性濾波算法同時估計各階球諧系數(shù)、自旋角速度以及探測器軌道參數(shù)。

2 小天體附近探測器軌道動力學

2.1 小天體引力場模型

小行星形狀不規(guī)則，內(nèi)部密度的分布也不均勻，這使得中心天體對環(huán)繞飛行器的引力作用也不嚴格等于等效的質心引力，而是等于中心天體各部分質量引力作用之和，即引力位函數(shù)：

其中，r為飛行器的位置，s為小天體內(nèi)某一點的位置，上述積分應該對整個中心天體體積進行積分，中心天體對環(huán)繞飛行器的引力可表示為：

由于我們無法得到精確的中心天體的質量分布結果，所以上式并不能用來計算中心天體的引力勢。常用的中心天體引力的計算模型有球諧系數(shù)模型［15］、多面體模型［16］、質點群模型等［17］。本文用球諧函數(shù)的級數(shù)展開式表示小行星的引力勢能，能夠更直觀地描述小天體模型精度階次。球諧系數(shù)模型如式（3）所示，小行星引力向量可以通過對勢能求梯度得到。

其中：

其中，φ和λ分別為飛行器在小天體固連坐標系下的緯度和經(jīng)度，φ′和λ′分別為積分單元的緯度和經(jīng)度，Cnm和Snm為非歸一化的引力場模型諧函數(shù)和為歸一化的引力場模型諧函數(shù)。Pnm為勒讓德多項式，δnm為克羅內(nèi)克符號（即n，m相同為1，不同為0），GM為引力常數(shù)，R為天體平均半徑。

由于Cnm和Snm的計算同樣需要知道質量分布，故一般也無法計算。目前而言，諧函數(shù)和一般是發(fā)射衛(wèi)星并通過定軌來獲得，從而通過換算得到Cnm和Snm。

m=n= 0的項為中心天體的質心引力部分，當n，m不為0 時，代表中心天體引力位偏離中心質點的程度。

2.2 小天體附近探測器軌道動力學模型

在小天體固連坐標系下，探測器軌道動力學模型如式（9）所示：

式中，xi，yi，zi為探測器i的三軸位置，ω為小天體的自旋角速度，nx，ny為未建模的干擾加速度，V為小天體的引力位勢函數(shù)球諧系數(shù)模型。

3 引力場與自旋參數(shù)感知系統(tǒng)

3.1 系統(tǒng)狀態(tài)模型

在小天體固連坐標系下建立引力場與自旋感知系統(tǒng)，待估系統(tǒng)變量為兩環(huán)繞探測器的位置、速度、引力場球諧系數(shù)、模型系數(shù)，以及自旋角速度，系統(tǒng)狀態(tài)變量如式（10）所示：

小天體各級引力場參數(shù)以及自旋角速度視為常數(shù)，其狀態(tài)模型如式（11）所示：

由式（10）～（11）可得到引力場反演的系統(tǒng)狀態(tài)模型如式（12）所示：

f(X)的具體表達式參照式（10）～（11）。

在兩探測器繞飛小天體的過程中，通過器間特高頻頻段無線電通信，實現(xiàn)相對距離測量，測量模型如式（13）所示：

式中，wρ為測量噪聲，在此假設為高斯白噪聲。

3.2 系統(tǒng)狀態(tài)估計方法

針對小天體引力場反演這一非線性估計問題，采用無跡卡爾曼濾波算法實現(xiàn)狀態(tài)最優(yōu)估計，給出tk時刻的系統(tǒng)先驗狀態(tài)Xk和狀態(tài)誤差協(xié)方差矩陣Pk，后續(xù)的估計流程如下：

（1）計算2n+ 1個sigma點及其權值。

式中，λ=α2(n+k)-n，α決定sigma 點的散布程度，通常取一個小的正值，k通常取0，β用來描述系統(tǒng)狀態(tài)的分析信息，高斯分布情況下，β最優(yōu)取值為2。

（2）計算sigma 點通過非線性函數(shù)f(X)的傳播結果。

進而可得：

（3）觀測更新。

4 仿真結果及分析

對目標小行星2016HO3 進行引力場與自旋參數(shù)感知仿真分析，假設目標天體等效半徑50m，密度1.26×103kg/m3（標準C 類小行星密度，2016HO3富含鐵質，實際密度要大），套用Bennu 的不規(guī)則引力場模型，前四階引力場參考模型參數(shù)如表1 所示。假設先驗引力場各階模型誤差為50%，兩探測器在本體系下的初始位置誤差標準差為10m，器間測距精度為10cm/1cm，探測器繞飛高度1km/500m/100m。

表1 前八階引力場參考模型參數(shù)Table 1 The first 8th order gravity model parameters

當測量精度為10cm 時，不同軌道高度的引力場反演仿真結果如圖1 所示。軌道高度1km 時，通過240h（包含不可見弧段時長）的測量，可以反演出二階引力場，軌道500m時，通過120h觀測可以反演出三階引力場，軌道高度100m時，通過240h觀測可以反演出四階引力場。不同軌道高度下探測器位置收斂精度近似，三軸位置精度約為20cm，如圖2所示。

圖1 測量精度10cm時不同軌道高度下感知的引力場模型參數(shù)誤差Fig.1 The gravity parameter errors with 10cm measurement accuracy in situations with different orbits

圖2 測量精度10cm時探測器軌道確定精度Fig.2 The orbit determination errors with 10cm measurement accuracy

考慮到2016HO3 的半徑約為50m，因此100m軌道高度已接近工程可實現(xiàn)的極限最低高度。為反演更高階次的引力場模型，需進一步提高測量精度。當測量精度達到1cm 時，100m 高度的探測器器間測距反演引力場結果如圖3 所示，對應的探測器軌道精度如圖4 所示。從仿真結果可以看出，當測量精度提升后，引力場反演精度明顯提高，在100m 軌道高度進行器間測量反演，可修正六階引力場參數(shù)。同時，探測器的定軌精度也從20cm 提高至5cm。通過引力場參數(shù)計算得到的真實引力加速度與感知加速度仿真對比如圖5所示。圖6給出了自旋角速度的誤差。假設光學觀測獲取的先驗自旋角速度誤差為1%，通過動力學感知方法，經(jīng)過170h 連續(xù)觀測后，自旋角速度誤差降低至0.01%。

圖3 測量精度1cm時100m軌道高度120h觀測感知的引力場模型參數(shù)誤差Fig.3 The gravity parameter errors with 1cm measurement accuracy and 100m-height orbit within 120h

5 結 論

圖4 測量精度1cm時探測器軌道確定精度Fig. 4 The orbit determination errors with 1cm measurement accuracy

圖5 小天體表面真實引力加速度與感知加速度對比Fig. 5 The comparion of real and perceived acceleration errors

圖6 自旋角速度感知誤差Fig. 6 The self-rotation speed error

本文研究了基于軌道動力學特性的小天體引力場與自旋角速度感知方法，利用兩個小天體環(huán)繞探測器之間的無線電測距信息，可同時估計探測器在小天體固連坐標系下的軌道、各階球諧系數(shù)模型參數(shù)以及小天體自旋角速度。影響引力場與自旋感知精度的因素主要包括探測器軌道高度和測量精度兩方面。仿真結果驗證了軌道高度越低，測量精度越高，可以感知的引力場階次越高，自旋角速度精度越高。精確的引力場信息是實施小天體探測任務的重要保障，本文的研究對空間環(huán)境感知、小天體探測中的參數(shù)估計有重要參考價值，創(chuàng)新點在于以器間測距信息為輸入對小天體的引力場進行了較為精確的估計，同時提高了探測器的定軌精度。但利用球諧系數(shù)計算小行星非球星引力存在著一定的缺點，小行星形狀越不規(guī)則時誤差越大，本文僅套用Bennu 的引力場模型進行了模擬仿真，后期可利用仿真或真實拍攝的圖片初步估計小行星的形狀，并采用更高階的非線性估計方法，進一步驗證該方法和提高對小行星引力場參數(shù)估計的精度。