基于軌道動(dòng)力學(xué)特性的小天體物理參數(shù)感知方法

羅 燕，秦 同，喬 棟

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

1 引 言

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

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

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

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

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

2 小天體附近探測器軌道動(dòng)力學(xué)

2.1 小天體引力場模型

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

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

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

其中：

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

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

m=n= 0的項(xiàng)為中心天體的質(zhì)心引力部分，當(dāng)n，m不為0 時(shí)，代表中心天體引力位偏離中心質(zhì)點(diǎn)的程度。

2.2 小天體附近探測器軌道動(dòng)力學(xué)模型

在小天體固連坐標(biāo)系下，探測器軌道動(dòng)力學(xué)模型如式（9）所示：

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

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

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

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

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

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

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

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

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

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

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

（1）計(jì)算2n+ 1個(gè)sigma點(diǎn)及其權(quán)值。

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

（2）計(jì)算sigma 點(diǎn)通過非線性函數(shù)f(X)的傳播結(jié)果。

進(jìn)而可得：

（3）觀測更新。

4 仿真結(jié)果及分析

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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