步士超，王宇凱,2，李文丹，李爽，牛升達(dá)

（1. 南京航空航天大學(xué) 航天新技術(shù)實(shí)驗(yàn)室，南京 210016；2. 上海微小衛(wèi)星工程中心，上海 201210；3. 上海衛(wèi)星工程研究所，上海 201109）

0 引 言

小行星因其包含了太陽(yáng)系起源、行星演化等方面的重要信息，以及其可能蘊(yùn)藏的豐富礦藏資源和潛在的對(duì)地球撞擊的威脅，已經(jīng)成為21世紀(jì)國(guó)際深空探測(cè)活動(dòng)的重點(diǎn)目標(biāo)。隨著我國(guó)深空探測(cè)持續(xù)不斷向前推進(jìn)，開(kāi)展對(duì)具有獨(dú)特科學(xué)價(jià)值的小行星的探測(cè)已經(jīng)成為深空探測(cè)技術(shù)發(fā)展的必然趨勢(shì)。雙小行星系統(tǒng)是眾多小行星類群中的一種特殊形式，是指兩顆小行星圍繞其共同的質(zhì)心旋轉(zhuǎn)，同時(shí)兩者的質(zhì)心圍繞太陽(yáng)做公轉(zhuǎn)運(yùn)動(dòng)。研究表明雙小行星系統(tǒng)是太陽(yáng)系中十分普遍的物理現(xiàn)象，大約有16%的近地小行星和主帶小行星屬于雙小行星系統(tǒng)[1]。對(duì)兩顆小行星構(gòu)成的雙小行星系統(tǒng)進(jìn)行探測(cè)具有多種探測(cè)單顆小行星所不具備地獨(dú)特科學(xué)價(jià)值，例如，許多雙小行星系統(tǒng)是由單顆小行星分裂而演化導(dǎo)致，于是，原本隱藏于小行星內(nèi)部的深層結(jié)構(gòu)便由于分裂而直觀地展現(xiàn)在分裂后的雙小行星系統(tǒng)的表面上，通過(guò)對(duì)這種結(jié)構(gòu)的直觀分析，可以獲得對(duì)小行星結(jié)構(gòu)構(gòu)造、動(dòng)力學(xué)特性等方面的更深層洞見(jiàn)。

雙小行星系統(tǒng)引力場(chǎng)與動(dòng)力學(xué)建模是進(jìn)行弱引力雙小行星系統(tǒng)探測(cè)任務(wù)設(shè)計(jì)時(shí)首先需要解決的問(wèn)題，更是研究雙小行星系統(tǒng)附近的飛行器軌道動(dòng)力學(xué)特性和相關(guān)任務(wù)軌道設(shè)計(jì)的基礎(chǔ)。

因此，本文擬采用不同的引力場(chǎng)簡(jiǎn)化模型，針對(duì)弱引力雙小行星系統(tǒng)的引力場(chǎng)進(jìn)行建模，對(duì)雙小行星系統(tǒng)動(dòng)力學(xué)特性進(jìn)行較全面的研究分析，為我國(guó)正在論證中的小行星采樣返回探測(cè)任務(wù)提供一定的技術(shù)支持。

1 限制性全三體問(wèn)題

1996年以前，針對(duì)形狀不規(guī)則的小行星的引力場(chǎng)建模問(wèn)題，學(xué)者們一般將其近似為三軸橢球體；1996年以后，多面體模型法在小行星的引力場(chǎng)建模中得到了廣泛的應(yīng)用，從而使人們能更精確地研究小行星附近的航天器軌道動(dòng)力學(xué)問(wèn)題。例如，Scheeres和Hudson等分別以形狀不規(guī)則的小行星4769 Castalia以及4179 Toutatis為例研究了其附近的軌道動(dòng)力學(xué)問(wèn)題[2-3]，詳細(xì)地分析了其重力場(chǎng)的特征情況、近小行星軌道動(dòng)力學(xué)方程等方面的問(wèn)題。在2005年后，D. J. Scheeres等學(xué)者針對(duì)雙小行星系統(tǒng)附近的航天器運(yùn)動(dòng)的軌道動(dòng)力學(xué)問(wèn)題等開(kāi)展了研究與分析[4-6]，深入分析了航天器在雙小行星系統(tǒng)中的軌道特點(diǎn)，并且研究了小行星系統(tǒng)三體問(wèn)題中小天體的非球形攝動(dòng)對(duì)軌道的影響。

但是，雙小行星系統(tǒng)附近的探測(cè)器動(dòng)力學(xué)建模問(wèn)題比較復(fù)雜，如果把兩個(gè)小行星的不規(guī)則引力場(chǎng)和自旋轉(zhuǎn)等因素都考慮進(jìn)去就是非常復(fù)雜的全三體問(wèn)題，如圖1所示，全三體問(wèn)題到目前為止基本上處于無(wú)解狀態(tài)。北京理工大學(xué)的李翔宇和喬棟等針對(duì)1996FG3雙小行星系統(tǒng)，采用雙球體模型，對(duì)逃逸軌道展開(kāi)了設(shè)計(jì)[7]，如圖2所示。但是，雙球體模型相對(duì)于實(shí)際的雙小行星系統(tǒng)的誤差比較大。Bellerose和Scheeres將雙小行星系統(tǒng)其中一個(gè)天體的質(zhì)量分布引入考慮，以均勻三軸橢球體和球體分別來(lái)近似兩個(gè)天體，構(gòu)成限制性全三體問(wèn)題（橢球體–球體模型），如圖3所示，并分析了雙小行星系統(tǒng)中航天器的運(yùn)動(dòng)規(guī)律，如平動(dòng)點(diǎn)、雅可比常數(shù)、周期軌道穩(wěn)定區(qū)域、雙星間轉(zhuǎn)移軌道等等[4,8-9]。

圖2 球體–球體模型示意圖Fig. 2 Schematic diagram of the sphere-sphere model

圖3 橢球體–球體模型示意圖Fig. 3 Schematic diagram of the ellipsoid-sphere model

1.1 采用橢圓積分計(jì)算的橢球體-球體模型

首先，我們以Bellerose和Scheeres提出的一個(gè)簡(jiǎn)化的雙小行星系統(tǒng)橢球體–球體模型為例（圖3）[8-9]?？紤]一個(gè)質(zhì)點(diǎn)或航天器在雙小行星系統(tǒng)的引力場(chǎng)范圍內(nèi)。M1為球體的質(zhì)量，M2為橢球體的質(zhì)量，rb為球體相對(duì)于橢球體的位置矢量，re為橢球體相對(duì)于系統(tǒng)質(zhì)心CM的位置矢量，rs為球體相對(duì)于系統(tǒng)質(zhì)心CM的位置矢量，為航天器相對(duì)于質(zhì)心CM的位置矢量。

假設(shè)航天器（或質(zhì)點(diǎn)）對(duì)雙小行星系統(tǒng)的運(yùn)動(dòng)沒(méi)有影響，那么航天器的運(yùn)動(dòng)方程為

設(shè)橢球體最長(zhǎng)軸的半長(zhǎng)軸為α，那么在此半徑下，系統(tǒng)平均角速度為在這里，我們以α和n為標(biāo)準(zhǔn)作為歸一化的長(zhǎng)度和角速度單位，那因此方程（1）歸一化后，變?yōu)?/p>

其中

式（3）是在如圖3所示的長(zhǎng)軸穩(wěn)定結(jié)構(gòu)下ω的取值情況所得。

1.2 采用橢圓積分計(jì)算的橢球體–橢球體模型

在前文的橢球體–球體模型的基礎(chǔ)上，本節(jié)提出了一種改進(jìn)的限制性橢球體–橢球體模型，如圖4所示。限制性橢球體–橢球體模型本質(zhì)上是將“橢球體–球體模型”中的“球體”替換為一個(gè)“限制性橢球體”，即將圖3中的球體M1改進(jìn)為一個(gè)限制性橢球體，將其三軸改為并且使用橢圓積分來(lái)計(jì)算其引力勢(shì)能，以便使更新后的雙小行星系統(tǒng)引力場(chǎng)模型能更準(zhǔn)確地描述探測(cè)器近雙小行星系統(tǒng)動(dòng)力學(xué)問(wèn)題。

圖4 限制性橢球體–橢球體模型示意圖Fig. 4 Schematic diagram of the restrictive ellipsoid-ellipsoid model

在此限制性橢球體–橢球體模型下

2 采用二階二次球諧函數(shù)模型的雙小行星系統(tǒng)引力場(chǎng)建模

在1996年，Kaula等學(xué)者建立了球諧函數(shù)模型理論[10]，其使用一系列的球諧函數(shù)展開(kāi)式來(lái)逼近小行星的引力場(chǎng)。球諧函數(shù)模型具有計(jì)算效率高、有明確解析表達(dá)式和各階攝動(dòng)大小明顯等優(yōu)點(diǎn)，因此其廣泛地應(yīng)用于天體力學(xué)和人造航天器軌道動(dòng)力學(xué)等方面。在前文中，都是采用橢圓積分來(lái)計(jì)算各種模型中的引力勢(shì)。但是，在計(jì)算過(guò)程中過(guò)多地引入積分計(jì)算環(huán)節(jié)，會(huì)大大地增大計(jì)算量和計(jì)算難度，消耗更多的計(jì)算時(shí)間。因此，在本節(jié)中，采用計(jì)算效率高，無(wú)積分環(huán)節(jié)的二階二次球諧函數(shù)模型來(lái)進(jìn)行雙小行星系統(tǒng)的引力場(chǎng)建模[11]。

小行星的引力勢(shì)函數(shù)可由球諧函數(shù)展開(kāi)式來(lái)逼近，如下式所示

胡維多等學(xué)者[11-12]研究發(fā)現(xiàn)采用二階二次球諧函數(shù)模型來(lái)計(jì)算小行星引力場(chǎng)既能得到比較精確的結(jié)果，還能兼顧計(jì)算的簡(jiǎn)便程度。在本節(jié)中，忽略球諧函數(shù)模型的高階項(xiàng)，取二階二次球諧函數(shù)模型，因此小行星引力勢(shì)為

其中：δ為緯度；λ為經(jīng)度。

當(dāng)以中心引力體的質(zhì)心作為坐標(biāo)系的原點(diǎn)時(shí)，那么二階二次球諧函數(shù)模型中的參數(shù)取值為而當(dāng)該中心引力體固連坐標(biāo)系的三軸沿著中心體的主慣量軸方向時(shí)，那么參數(shù)因此，在這種坐標(biāo)系下，利用二階二次球諧函數(shù)模型計(jì)算的天體引力場(chǎng)可以簡(jiǎn)化為

其中，兩個(gè)系數(shù)C20和C22與中心引力體的3個(gè)主慣量軸的轉(zhuǎn)動(dòng)慣量有如下關(guān)系[11-12]

假設(shè)中心天體為橢球體，其3個(gè)主軸分別為α、β和γ，那么各個(gè)坐標(biāo)軸的轉(zhuǎn)動(dòng)慣量為

2.1 球體–球體模型

在球體–球體模型中，天體M1和M2的引力勢(shì)均用球體模型來(lái)計(jì)算。則天體M1和M2的勢(shì)函數(shù)為

2.2 橢球體–球體模型

Bellerose在橢球體–球體模型中采用橢圓積分來(lái)計(jì)算橢球體的引力勢(shì)。在本節(jié)中，改進(jìn)采用計(jì)算效率高，無(wú)積分環(huán)節(jié)的二階二次球諧函數(shù)模型計(jì)算天體的引力勢(shì)。則天體M1和M2的勢(shì)函數(shù)為

2.3 橢球體–橢球體模型

在第2.2節(jié)的橢球體–球體模型的基礎(chǔ)上，本節(jié)提出了一種改進(jìn)的限制性橢球體–橢球體模型。限制性橢球體–橢球體模型本質(zhì)上是將“橢球體–球體模型”中的“球體”替換為一個(gè)“限制性橢球體”，即將球體M1改進(jìn)為一個(gè)限制性橢球體，將其三軸改為并且使用二階二次球諧函數(shù)模型來(lái)計(jì)算其引力勢(shì)能，以便使更新后的雙小行星系統(tǒng)引力場(chǎng)模型能更準(zhǔn)確地描述探測(cè)器近雙小行星系統(tǒng)動(dòng)力學(xué)問(wèn)題。因此，在此模型下，天體M1和M2的勢(shì)函數(shù)為

3 仿真結(jié)果

在前文中，針對(duì)弱引力雙小行星系統(tǒng)引力場(chǎng)建模問(wèn)題，分別采用復(fù)雜度和精度依次遞增的球體–球體模型、橢球體–球體模型和限制性橢球體–橢球體模型進(jìn)行建模，并分別采用橢圓積分和二階二次球諧函數(shù)來(lái)計(jì)算其引力勢(shì)能。本節(jié)針對(duì)平面圓形限制性三體問(wèn)題（Circular Restricted Three-Body Problem，CRTBP），以雙小行星系統(tǒng)1999KW4為例，分別對(duì)前文給出的不同弱引力雙小行星系統(tǒng)引力場(chǎng)模型進(jìn)行勢(shì)能的仿真分析。雙小行星系統(tǒng)1999KW4限制性橢球體–橢球體模型歸一化后的參數(shù)如表1所示。

表1 雙小行星系統(tǒng)1999KW4歸一化后參數(shù)Table 1 The normalized parameters of the restricted ellipsoidellipsoid model of binary system 1999KW4

雙小行星系統(tǒng)1999KW4中天體M1和M2的勢(shì)函數(shù)分別為等效勢(shì)能函數(shù)為

則不同引力場(chǎng)模型下的等效勢(shì)能函數(shù)曲面分別如圖5和圖6所示。

由圖5和圖6可以看出，等效勢(shì)能曲面有2個(gè)極大值點(diǎn)和3個(gè)鞍點(diǎn)（對(duì)應(yīng)5個(gè)平動(dòng)點(diǎn)），圖中也畫(huà)出了在不同高度（對(duì)應(yīng)不同的雅克比能量）的等高線，這些等高線就是對(duì)應(yīng)雅克比能量下的零速度曲線，圖7以等高線圖的形式更清楚地描述了零速度曲線隨雅克比能量變化的情況。由圖7可以看出，零速度曲線關(guān)于x軸對(duì)稱，隨著雅克比能量的增大，對(duì)應(yīng)的可運(yùn)動(dòng)區(qū)域增大，并且該區(qū)域的拓?fù)浣Y(jié)構(gòu)也發(fā)生了變化。

圖6 平面圓形限制性三體問(wèn)題等效勢(shì)能函數(shù)曲面（二階二次球諧函數(shù)計(jì)算引力勢(shì)能）Fig. 6 The effective potential of CRTBP （ellipsoid potential energy is written in terms of the second degree second order spherical harmonic function）

在CRTBP中，存在著5個(gè)平動(dòng)點(diǎn)（拉格朗日點(diǎn)），并且拉格朗日點(diǎn)是梯度矢量場(chǎng)的零點(diǎn)，見(jiàn)式（31），該矢量場(chǎng)也是等效勢(shì)能所產(chǎn)生的加速度場(chǎng)，不同引力勢(shì)模型下的等效勢(shì)能的梯度矢量場(chǎng)及5個(gè)拉格朗日點(diǎn)如圖8所示。

圖7 平面圓形限制性三體問(wèn)題等效勢(shì)能函數(shù)的等高線圖Fig. 7 The contours of effective potential in CRTBP

接下來(lái)以1999KW4為例分別比較不同模型下的平動(dòng)點(diǎn)位置偏差。首先，計(jì)算由橢圓積分計(jì)算引力勢(shì)的雙橢球體模型下的平動(dòng)點(diǎn)位置坐標(biāo)，并與文獻(xiàn)[9]中橢球體-球體模型的結(jié)果進(jìn)行對(duì)比，結(jié)果如表2所示。然后，計(jì)算不同模型下分別由橢圓積分和由二階二次球諧函數(shù)計(jì)算引力勢(shì)的平動(dòng)點(diǎn)位置坐標(biāo)，如表3、表4所示。

圖8 等效勢(shì)能的梯度矢量場(chǎng)及5個(gè)拉格朗日點(diǎn)（二階二次球諧函數(shù)計(jì)算引力勢(shì)能Fig. 8 The gradient vector field of and five Lagrangian points

觀察表2可知，雙橢球模型與文獻(xiàn)[9]中橢球-球模型的平動(dòng)點(diǎn)位置最大偏差為0.178 6 %，考慮到雙橢球模型更接近于1999KW4的實(shí)際形狀，因此用雙橢球體模型是較為精確的。觀察表3和表4，可以發(fā)現(xiàn)，在利用橢球體-球體模型和橢球體-橢球體模型進(jìn)行雙小行星系統(tǒng)的引力場(chǎng)建模時(shí)，采用橢圓積分進(jìn)行計(jì)算的結(jié)果和采用二階二次球諧函數(shù)模型進(jìn)行計(jì)算的結(jié)果是十分相近的，最大偏差為0.049 7 %，因此采用二階二次球諧函數(shù)模型進(jìn)行引力勢(shì)的計(jì)算也是能得到很精確的結(jié)果的。而比較計(jì)算時(shí)間，從表3中，可以看出，在橢球體-球體模型中，采用橢圓積分進(jìn)行計(jì)算需要消耗0.53 s的時(shí)間，而采用二階二次球諧函數(shù)模型進(jìn)行計(jì)算消耗時(shí)間僅僅是0.015 s，遠(yuǎn)遠(yuǎn)小于0.53 s。而在橢球體-橢球體模型中，采用橢圓積分進(jìn)行計(jì)算需要計(jì)算更多的積分環(huán)節(jié)，從表4可以看出，其所消耗的時(shí)間也從0.53 s增加到了1.389 s，而采用二階二次球諧函數(shù)模型進(jìn)行計(jì)算所消耗的時(shí)間并沒(méi)有明顯的增加，只有0.016 s，遠(yuǎn)遠(yuǎn)小于1.389 s。因此，綜上所述，可以證明本文提出的二階二次球諧函數(shù)計(jì)算引力勢(shì)的橢球體-橢球體模型計(jì)算精度高，復(fù)雜程度更低，計(jì)算量更少，計(jì)算速度更快，是十分有效的。

表2 橢圓積分計(jì)算引力勢(shì)的平動(dòng)點(diǎn)位置坐標(biāo)對(duì)比Table 2 Comparison of coordinates of libration points calculated by elliptic integrals

表3 橢球體-球體模型平動(dòng)點(diǎn)位置坐標(biāo)對(duì)比Table 3 Comparison of coordinates of libration points in ellipsoid - sphere model

表4 橢球體-橢球體模型平動(dòng)點(diǎn)位置坐標(biāo)對(duì)比Table 4 Comparison of coordinates of libration points in ellipsoid - ellipsoid model

4 結(jié) 論

由于雙小行星系統(tǒng)附近的探測(cè)器動(dòng)力學(xué)建模問(wèn)題比較復(fù)雜，而且如果把兩個(gè)小行星的不規(guī)則引力場(chǎng)和自旋轉(zhuǎn)等因素都考慮進(jìn)去就是非常復(fù)雜的全三體問(wèn)題，其目前為止基本上處于無(wú)解狀態(tài)。因此，針對(duì)弱引力雙小行星系統(tǒng)的引力場(chǎng)建模問(wèn)題，本文采用復(fù)雜度和精度依次遞增的球體–球體模型、橢球體–球體模型和改進(jìn)的限制性橢球體–橢球體模型來(lái)進(jìn)行引力場(chǎng)建模，從而比較精確地刻畫(huà)雙小行星系統(tǒng)和探測(cè)器構(gòu)成的限制性全三體問(wèn)題的動(dòng)力學(xué)模型。并針對(duì)不同的引力場(chǎng)建模方法進(jìn)行了仿真，以雙小行星系統(tǒng)1999KW4為例，給出了不同引力場(chǎng)模型下的等效勢(shì)能曲面及等高線圖，最后還給出了等效勢(shì)能的梯度矢量場(chǎng)及5個(gè)拉格朗日點(diǎn)。最后針對(duì)不同模型下的平動(dòng)點(diǎn)位置坐標(biāo)偏差進(jìn)行了比較。本文提出的二階二次球諧函數(shù)計(jì)算引力勢(shì)的橢球體–橢球體模型計(jì)算精度高，復(fù)雜程度更低，計(jì)算量更少，計(jì)算速度更快，能夠較精確的對(duì)雙小行星系統(tǒng)進(jìn)行引力場(chǎng)建模。未來(lái)更進(jìn)一步，還可以結(jié)合多面體模型法和球諧函數(shù)模型法來(lái)更精確地描述雙小行星系統(tǒng)中不規(guī)則天體的引力場(chǎng)模型。

[1]Pravec P，Harris A W. Binary asteroid population：1. angular momentum content [J]. Icarus，2007，190（1）：250-259.

[2]Scheeres D J，Ostro S J，Hudson R S，et al. Orbits close to asteroid 4769 Castalia [J]. Icarus，1996，121：67-87.

[3]Scheeres D J，Ostro S J，Hudson R S，et al. Dynamics of orbits close to asteroid 4179 Toutatis [J]. Icarus，1998，132：53-79.

[4]Bellerose J，Scheeres D J. Stability of equilibrium points in the restricted full three body problem[J]. Acta Astronautica，2007，60：141-152.

[5]Fahnestock E G，Scheeres D J. Simulation of the full two rigid body problem using polyhedral mutual potential and potential derivatives approach [J]. Celestial Mechanics and Dynamical Astronomy，2006，96：317-339.

[6]Fahnestock E G，Lee T，Leok M，et al. Polyhedral potential and variational integrator computation of the full two body problem[C]//The 2006 AIAA/AAS Astrodynamics Specialist Meeting.[S.l.]：AIAA，2006.

[7]李翔宇，喬棟，崔平遠(yuǎn). 雙體小行星1996FG3捕獲與逃逸軌道設(shè)計(jì)[C]//上海航天技術(shù)研究院第二屆小行星探測(cè)學(xué)術(shù)研討論文集.上海：上海航天技術(shù)研究院，2014.Li X Y，Qiao D，Cui P Y. The capture and escape trajectory design of the binary asteroid 1996FG3 [C]// The 2nd Shanghai Academy of Spaceflight Technology Asteroid Exploration Conference. Shanghai：Shanghai Academy of spaceflight Technology，2014.

[8]Bellerose J，Scheeres D J. Restricted full three-body problem：application to binary system 1999 KW4 [J]. Journal of Guidance，Control，and Dynamics，2008，31（1）：162-171.

[9]Bellerose J. The restricted full three body problem：applications to binary asteroid exploration[D]. Ann Arbor：University OF Michigan，2008.

[10]Kaula W M. Theory of satellite geodesy：applications of satellites to geodesy [M]. Mineola：Dover publications，2000：4-8.

[11]Hu W. Orbital motion in uniformly rotating second degree and order gravity fields [D]. Ann Arbor：the University of Michigan，2002.

[12]Hu W，Scheeres D J. Numerical determination of stability regions for orbital motion in uniformly rotating second degree and order gravity fields [J]. Planetary and Space Science，2004，52（8）：685-692.