杜燕茹，李翔宇，韓宏偉，喬 棟

（1.北京理工大學(xué)宇航學(xué)院，北京 100081；2.深空自主導(dǎo)航與控制工信部重點實驗室，北京 100081）

引 言

作為小行星的重要存在形式，雙體小行星在近地小行星、主帶小行星與特洛伊小行星中都已被觀測到[1]。在近地小行星中，雙體小行星的數(shù)量約占到其總數(shù)的15±4%[2-3]。由于雙體小行星系統(tǒng)中天體間所具有的獨特運動特性可為研究行星系統(tǒng)的演化提供重要線索，因而成為目前小行星探測的熱點目標(biāo)。雙體小行星模型的不確定性及其附近復(fù)雜多變的動力學(xué)環(huán)境對雙體小行星探測任務(wù)的設(shè)計提出了挑戰(zhàn)。

近年來雙體小行星的穩(wěn)定特性及其附近復(fù)雜的軌道動力學(xué)問題引起了學(xué)者們的廣泛關(guān)注。Scheeres[4]基于全二體模型在雙體小行星能量與角動量守恒條件下，研究了雙體小行星系統(tǒng)旋轉(zhuǎn)運動與平移運動中能量與角動量的轉(zhuǎn)換關(guān)系，并給出了雙體小行星系統(tǒng)Hill穩(wěn)定性準(zhǔn)則和發(fā)生碰撞的條件?；谇蚝蜋E球構(gòu)成的雙體小行星系統(tǒng)，Scheeres提出了評估平衡狀態(tài)下能量穩(wěn)定性與譜穩(wěn)定性的方法[5]。Bellerose 等基于橢球-球構(gòu)成的全二體問題[6-7]，研究了相對平衡狀態(tài)與穩(wěn)定性問題，推導(dǎo)了穩(wěn)定性的解析判據(jù)，并針對不同構(gòu)型下的雙體小行星系統(tǒng)求解出了穩(wěn)定區(qū)域，研究發(fā)現(xiàn)：與限制性三體問題中的穩(wěn)定區(qū)域相比，橢球體-球體系統(tǒng)的穩(wěn)定范圍更小[8]。

本文在前人研究的基礎(chǔ)上，通過將橢球體與球體模型拓展至雙橢球體模型構(gòu)成的全二體問題，研究了雙體小行星系統(tǒng)的平衡態(tài)及穩(wěn)定性，并分析了不同物理參數(shù)對雙體小行星平衡態(tài)穩(wěn)定性的影響，將為未來小行星探測提供理論基礎(chǔ)。

1 雙體小行星動力學(xué)與運動方程

1.1 全二體問題

全二體問題研究兩個主天體在考慮形狀因素以及質(zhì)量分布的情況下，受相互引力作用而進(jìn)行的運動，且運動被限制在一個平面內(nèi)，如圖1所示。其中，O點為系統(tǒng)的質(zhì)心，P1和P2為構(gòu)成雙體小行星系統(tǒng)的兩個形狀任意、質(zhì)量分別為M1與M2的主天體。在該雙體小行星系統(tǒng)中，分別引入慣性坐標(biāo)系OXYZ與旋轉(zhuǎn)坐標(biāo)系Oxyz。慣性坐標(biāo)系OXYZ的坐標(biāo)原點為雙體小行星系統(tǒng)質(zhì)心O，X軸指向慣性空間的固定方向，Z軸指向主天體軌道運動的角動量方向，X、Y、Z軸構(gòu)成右手坐標(biāo)系。旋轉(zhuǎn)坐標(biāo)系Oxyz的坐標(biāo)原點為系統(tǒng)質(zhì)心O，x軸由主天體P1質(zhì)心指向主天體P2質(zhì)心，z軸指向主天體軌道運動的角動量方向，x、y、z構(gòu)成右手坐標(biāo)系。此外，分別以小行星的質(zhì)心為原點建立小行星本體坐標(biāo)系，選取本體系的xi、yi和zi(i=1,2)軸與小行星主軸方向一致，其中zi軸垂直于軌道面[9]。假設(shè)雙體小行星系統(tǒng)質(zhì)心O的運動相對于小行星間的相對運動可忽略，且小行星僅圍繞慣性主軸進(jìn)行自旋，因此，雙體小行星間的相對運動可簡化為用參數(shù)r，?i(i=1,2)，θ表示的四自由度問題，其中r為兩小行星質(zhì)心間的距離，?i(i=1,2)為小行星本體系xi(i=1,2)軸與小行星質(zhì)心連線的夾角，θ為質(zhì)心連線相對于慣性系的旋轉(zhuǎn)角度。

圖1 平面全二體問題示意圖Fig.1 Planar full two-body problem

1.2 雙體小行星系統(tǒng)勢能與動能

對于空間全二體問題，慣性系中雙體小行星系統(tǒng)的勢能為[5]

方程(1)給出了將慣性矩陣展開至二階形式的系統(tǒng)勢能。其中：Tr 表示矩陣的跡；為兩個小天體質(zhì)心在慣性系中的相對位置矢量；1,2)為單位化后的慣性矩陣，Mi(i=1,2)為兩個小天體的質(zhì)量；Ai(i=1,2)為由慣性坐標(biāo)系轉(zhuǎn)換至小天體本體系的旋轉(zhuǎn)矩陣。對于橢球體模型，假設(shè)質(zhì)量均勻分布且橢球體的半長軸分別為αi>βi>γi(i=1,2)，則繞主軸的慣性矩可以表示為

根據(jù)系統(tǒng)動能的一般性表達(dá)形式，在慣性系OXYZ中可以推導(dǎo)得到雙體小行星系統(tǒng)的動能，在假設(shè)主天體的旋轉(zhuǎn)僅圍繞主慣性矩發(fā)生的前提下，系統(tǒng)動能在旋轉(zhuǎn)坐標(biāo)系下的表達(dá)式為

系統(tǒng)的能量可以表示為E=T+V。

1.3 雙體小行星系統(tǒng)運動方程

對于雙體小行星系統(tǒng)等經(jīng)典物理系統(tǒng)，定義拉格朗日量L為系統(tǒng)動能與系統(tǒng)勢能之差，即

選取兩主天體質(zhì)心間的相對距離r、雙體小行星系統(tǒng)角速度以及主天體相對于旋轉(zhuǎn)坐標(biāo)系的旋轉(zhuǎn)角度?i用于表示平面全二體問題中的運動方程。

對于q=[r,θ,?1,?2]，根據(jù)拉格朗日方程?L/?qi=可得雙體小行星系統(tǒng)的運動方程為

2 雙橢球體系統(tǒng)平衡態(tài)與穩(wěn)定性分析

2.1 雙橢球體系統(tǒng)平衡態(tài)的狀態(tài)約束

在求解雙體小行星系統(tǒng)平衡態(tài)對應(yīng)的約束方程時，角動量是其中的一個關(guān)鍵參數(shù)。分析式(1)與式(3)中的雙體小行星系統(tǒng)動能表達(dá)式與勢能表達(dá)式可發(fā)現(xiàn)，系統(tǒng)勢能與動能的值與雙體小行星系統(tǒng)旋轉(zhuǎn)角度θ無關(guān)，即雙體小行星系統(tǒng)角動量守恒，角動量K的表達(dá)式為

求解雙體小行星系統(tǒng)的平衡狀態(tài)通常將全二體問題運動學(xué)方程中關(guān)于時間求導(dǎo)數(shù)的變量設(shè)為零，可求解平衡狀態(tài)所滿足的約束。因此雙體小行星系統(tǒng)處于平衡狀態(tài)的條件為，系統(tǒng)能量E=T+V需滿足對各變量求偏導(dǎo)的值均為零。

由于式(1)的勢能表達(dá)式僅適用于平面全二體問題，因此系統(tǒng)的能量表達(dá)式僅適用于平面全二體問題。平面全二體問題的平衡態(tài)需滿足的約束為

通過式(6)中的4個約束可以確定雙體小行星系統(tǒng)平衡狀態(tài)應(yīng)滿足的3個約束方程為

條件③：Er=0。

本文僅考慮?1=?2=0 的條件下，雙體小行星系統(tǒng)的相對運動情況。此時相對平衡態(tài)下對應(yīng)的雙體小行星系統(tǒng)角動量平方的值K02與系統(tǒng)能量的大小E0可以分別從式(6)計算得出。

2.2 雙橢球體系統(tǒng)平衡態(tài)穩(wěn)定性的分析

在進(jìn)行雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性分析之前，首先需要引入系統(tǒng)零速度狀態(tài)的概念，也即兩天體相對運動為零時狀態(tài)變量構(gòu)成的等式約束。將相對平衡狀態(tài)r0處的雙體小行星系統(tǒng)能量E0與角動量的平方K02帶入零速度狀態(tài)公式，且將平衡狀態(tài)對應(yīng)的距離r0作為參數(shù)，可以獲得雙體小行星系統(tǒng)運動可達(dá)范圍的控制方程。由零速度狀態(tài)公式推導(dǎo)得到的控制方程為

其中，A～G為該五次方程的各階系數(shù)，受小行星質(zhì)量分布以及平衡狀態(tài)時兩小行星的距離r0的影響。計算驗證式(7)的五次方程有二重根r=r0，故可因式分解得到簡化的控制方程為

利用勞斯-赫爾維茨判據(jù)可得結(jié)論：當(dāng)式(8)中A> 0 時，該方程的根均分布在復(fù)平面的左半平面；當(dāng)A<0 時該方程在復(fù)平面的右半平面有一個根，即存在一個正實根。對該三次不等式進(jìn)行討論研究，可獲得雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的3種情況：①運動可達(dá)范圍無上界的不穩(wěn)定平衡態(tài)；②運動可達(dá)范圍有上界的不穩(wěn)定平衡態(tài)；③穩(wěn)定的平衡態(tài)。

由于平衡態(tài)的穩(wěn)定性本質(zhì)上關(guān)注的是能量穩(wěn)定性與Hill 穩(wěn)定性。能量穩(wěn)定性為一類基礎(chǔ)的穩(wěn)定性類型，決定了系統(tǒng)能否通過能量耗散的方式達(dá)到能量較低的狀態(tài)。雙體小行星系統(tǒng)具有能量穩(wěn)定性意味著該雙體小行星系統(tǒng)已達(dá)到在給定角動量下能量最低的構(gòu)型狀態(tài)。Hill穩(wěn)定性描述雙體小行星系統(tǒng)在相互引力作用下是否具有發(fā)生逃逸能力的穩(wěn)定性，與能量穩(wěn)定性相比，Hill穩(wěn)定性研究對象不是平衡狀態(tài)附近的運動，而是關(guān)注雙體小行星系統(tǒng)運動的絕對界限。所以以上3 種平衡態(tài)穩(wěn)定性根據(jù)其參數(shù)特性又可歸結(jié)為：①Hill 不穩(wěn)定平衡態(tài)；②Hill 穩(wěn)定、能量不穩(wěn)定平衡態(tài)；③能量穩(wěn)定平衡態(tài)。下面對這3類穩(wěn)定平衡態(tài)分別進(jìn)行分析。

1）Hill不穩(wěn)定平衡態(tài)

A>0 時，不等式的根均分布于復(fù)平面的左半平面，則式(8)中控制方程始終成立，即r的取值范圍沒有約束，r的取值范圍無上界，對應(yīng)的雙體小行星系統(tǒng)在受擾動的情況下可發(fā)生逃逸，且r0處對應(yīng)的自由能始終為正。

A=0時，若2A+B≠0，通過推導(dǎo)分析可知簡化后控制不等式始終成立。若2A+B=0，則三次不等式可以簡化為2Gr/r0+G≥0，顯然控制不等式同樣恒成立。

圖2 運動可達(dá)范圍無上界的不穩(wěn)定平衡態(tài)Fig.2 Unstable equilibrium，the motion is allowable on all parts of the curve beneath the equilibrium energy

由以上分析可知，當(dāng)A>0與A=0時相應(yīng)的控制方程始終成立，即對應(yīng)的雙體小行星系統(tǒng)在受擾動的情況下可發(fā)生逃逸，運動范圍r無上界，該情況為運動可達(dá)范圍無上界的不穩(wěn)定平衡態(tài)。圖2所示為系統(tǒng)能量與相對距離的關(guān)系，平衡狀態(tài)位于r/r0=1，運動可達(dá)區(qū)域為位于平衡狀態(tài)能量E0之下且零能量曲面之上對應(yīng)的r的取值范圍。

2）Hill穩(wěn)定、能量不穩(wěn)定平衡態(tài)

在A<0 的情況下，控制方程有一實數(shù)根均分布于復(fù)平面的右半平面，其余兩根分布在復(fù)平面的左半平面。假設(shè)控制方程的正實根為X，若X> 1 即r >r0，此時r有上界M，且上界M=Xr0應(yīng)滿足M>r0。如果對平衡態(tài)r/r0=1處的雙體小行星系統(tǒng)施加擾動，則運動可達(dá)區(qū)域可能達(dá)到上界M，即雙體小行星系統(tǒng)內(nèi)兩個主天體的相對距離的值可能在r0和M之間變動，運動可達(dá)區(qū)域為位于平衡狀態(tài)能量之下的曲線部分對應(yīng)的距離。在該情況下，r/r0=1處的雙體小行星系統(tǒng)的系統(tǒng)總能量為負(fù)值，但為不穩(wěn)定狀態(tài)。雙體小行星系統(tǒng)能量與相對距離的關(guān)系如圖3所示。

圖3 雙橢球體系統(tǒng)下運動受限的平衡態(tài)Fig.3 Equilibrium with limited movement in the double ellipsoid system

3）能量穩(wěn)定平衡態(tài)

在A< 0 且X≤ 1 即r

圖4 雙橢球體系統(tǒng)能量最小的平衡態(tài)Fig.4 Minimal energy equilibrium in the double ellipsoid system

3 雙體小行星系統(tǒng)物理參數(shù)對平衡態(tài)穩(wěn)定性影響分析

基于上述關(guān)于雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的結(jié)論，研究雙體小行星系統(tǒng)平衡態(tài)受小行星形狀參數(shù)、質(zhì)量分?jǐn)?shù)以及雙體小行星系統(tǒng)平衡態(tài)時小行星間相對距離等物理參數(shù)的影響。

假設(shè)小行星Pi為密度均勻一致的雙橢球體模型，Pi的橢球體三軸滿足αi>βi>γi。根據(jù)2.2小節(jié)中對雙體小行星系統(tǒng)平衡狀態(tài)穩(wěn)定性的分類，分別選取形狀參數(shù)一致的雙體小行星系統(tǒng)，并通過調(diào)整雙體小行星主天體的構(gòu)型與質(zhì)心間距離，分別獲得3種不同穩(wěn)定性情況的系統(tǒng)平衡狀態(tài)。各系統(tǒng)主天體形狀參數(shù)如表1所示。

表1 3種不同穩(wěn)定性情況參數(shù)選取Table 1 Parameters for three different kinds of equilibrium m

隨后對主天體P2質(zhì)心沿徑向施加5 N·s 的沖量，并以施加沖量后的雙體小行星系統(tǒng)狀態(tài)參數(shù)為狀態(tài)初值，依照第2節(jié)中提供的雙體小行星系統(tǒng)運動方程沿時間正向遞推，從而研究雙體小行星受到小擾動后的運動變化。遞推時長選取為25 h。

圖5為分別針對具有能量穩(wěn)定性、Hill 穩(wěn)定且能量不穩(wěn)定性、Hill不穩(wěn)定性的平衡狀態(tài)雙體小行星系統(tǒng)，施加小擾動力后主天體質(zhì)心間的運動變化規(guī)律。由仿真結(jié)果可知，對于處于平衡狀態(tài)的雙體小行星系統(tǒng)受到擾動后主天體質(zhì)心間的距離變化可依據(jù)不同的穩(wěn)定特性分為3 類。針對能量穩(wěn)定系統(tǒng)施加小擾動時，雙體小行星系統(tǒng)主天體質(zhì)心間距離具有保持原平衡狀態(tài)不改變的特性。在平衡狀態(tài)的雙體小行星系統(tǒng)具有Hill穩(wěn)定的情況中，施加小擾動使得兩個主天體質(zhì)心間的相對距離在一定范圍內(nèi)改變，系統(tǒng)呈現(xiàn)有界的相對運動。對于平衡狀態(tài)下Hill不穩(wěn)定的雙體小行星系統(tǒng)，施加極其微小的擾動即可使系統(tǒng)主天體質(zhì)心間的相對距離達(dá)到無窮遠(yuǎn)，即發(fā)生逃逸。

圖5 能量穩(wěn)定、Hill穩(wěn)定但能量不穩(wěn)定、Hill不穩(wěn)定雙體小行星系統(tǒng)質(zhì)心距離變化Fig.5 Centroid distance of energy stable，Hill stable but energy unstable，Hill unstable binary asteroid systems

由2.2節(jié)中關(guān)于能量穩(wěn)定性與Hill穩(wěn)定性的結(jié)論，繪制質(zhì)量分?jǐn)?shù)的值在[0,1]之間變化的不同橢球體形狀下曲線。圖6～8 中的距離長度已相對于兩橢球體的半長軸之和α1+α2單位化。將雙球體模型的曲線與不同形狀參數(shù)下雙橢球體的曲線結(jié)果對比，顯示出不同形狀對應(yīng)結(jié)果的差異。

圖6～8 展示了考慮兩天體形狀為橢球時，因形狀參數(shù)的改變而產(chǎn)生的能量穩(wěn)定極限、Hill穩(wěn)定極限以及角動量相同時對應(yīng)的共軛解的變化。由圖6～8可知，在改變形狀參數(shù)的過程中穩(wěn)定性界限的特性未發(fā)生質(zhì)的改變，比如對于Hill穩(wěn)定性的極限，存在某一部分使得對于任意天體質(zhì)心間距離，雙體小行星系統(tǒng)均具有Hill穩(wěn)定性，對于能量穩(wěn)定性則不存在這種情況；Hill穩(wěn)定極限對應(yīng)的兩天體間距離始終小于能量穩(wěn)定極限對應(yīng)的距離。微小的橢球率可一定程度改變系統(tǒng)的穩(wěn)定性極限，且兩主天體三軸比例不同的情況會使得對應(yīng)的能量穩(wěn)定性極限、Hill穩(wěn)定性極限以及共軛解的曲線不再具有對稱性。下面分析小行星形狀參數(shù)變化對雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的影響。假設(shè)主天體形狀參數(shù)滿足βi=γi，選取小行星歸一化后的的長半軸為α1=6，α2=4，兩個小行星質(zhì)心間的初始平衡態(tài)距離為r*=10。

圖6 不同形狀參數(shù)天體的能量穩(wěn)定質(zhì)心間距離極限隨質(zhì)量分?jǐn)?shù)的變化Fig.6 Centroid distance limit of energy stable system for a range of different shape bodies

圖7 不同形狀參數(shù)天體Hill穩(wěn)定質(zhì)心間距離極限隨質(zhì)量分?jǐn)?shù)的變化Fig.7 Centroid distance limit of Hill stable system for a range of different shape bodies

圖8 與接觸系統(tǒng)具有相同角動量的不同形狀天體的共軛解隨質(zhì)量分?jǐn)?shù)的變化Fig.8 Conjugate stable solutions for a range of different shape bodies.The conjugate solution has the same angular momentum as the solution for a touching system

圖9給出了控制方程的三次項系數(shù)A與控制方程實數(shù)解Xr的變化曲線。由圖9可知，每一條曲線均具有相似的變化趨勢，即當(dāng)主天體P1短半軸與長半軸的比值β1/α1極小時，控制方程三次項系數(shù)A>0，此時雙體小行星系統(tǒng)處于可能發(fā)生逃逸的不穩(wěn)定狀態(tài)；隨著β1/α1略微增大，存在小部分區(qū)間可滿足A<0且Xr>1，該情況下雙體小行星系統(tǒng)運動可達(dá)范圍具有上界M=Xr r*，為能量不穩(wěn)定的相對平衡狀態(tài)；主星形狀參數(shù)β1/α1繼續(xù)增大，雙體小行星系統(tǒng)的結(jié)果為A<0且Xr<1，即平衡狀態(tài)能量穩(wěn)定。而針對不同取值的β2/α2可發(fā)現(xiàn)，主天體P2短半軸與長半軸的比值的增大將使得雙體小行星系統(tǒng)達(dá)到穩(wěn)定狀態(tài)所需的小行星P1的形狀參數(shù)比例β1/α1更小。

下面分析不同小行星質(zhì)量分?jǐn)?shù)對雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的影響。令主天體的三軸比例滿足1:0.75:0.5，選取小行星P2的長半軸為α2=1，兩個小行星質(zhì)心間的距離為r*=10，給出控制方程的三次項系數(shù)A與控制方程實數(shù)解Xr的變化曲線如圖10所示。圖中變化曲線與圖9(b)中變化趨勢相近，出現(xiàn)該現(xiàn)象的原因為在假設(shè)密度相同的條件下，β1/α1比值的增加將導(dǎo)致質(zhì)量分?jǐn)?shù)的增大。并且在質(zhì)量分?jǐn)?shù)極大時，即小行星P1的質(zhì)量M1遠(yuǎn)大于M2時，雙體小行星平衡狀態(tài)對應(yīng)的A與Xr增大，雙體小行星系統(tǒng)開始由能量穩(wěn)定的相對平衡狀態(tài)向Hill 穩(wěn)定且能量不穩(wěn)定狀態(tài)轉(zhuǎn)變。

圖9 β1/α1與β2/α2對控制方程的三次項系數(shù)A與控制方程實數(shù)解Xr的影響Fig.9 Effect of β1/α1 and β2/α2 on oefficient A and real solution Xr in the control equation

圖10 質(zhì)量分?jǐn)?shù)對雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的影響Fig.10 Influence of the mass ratio on the stability of equilibrium

最后分析小行星間相對距離變化對雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的影響。假設(shè)兩個小行星具有相同的形狀參數(shù)，選取小行星的長半軸α1=α2=1，為保證兩個小行星初始狀態(tài)不為碰撞狀態(tài)，質(zhì)心間的距離取值需r*>α1+α2，即在[2,10]的區(qū)間內(nèi)改變雙體小行星系統(tǒng)的質(zhì)心距離。

控制方程的三次項系數(shù)A與控制方程實數(shù)解Xr的隨r*變化的曲線如圖11所示。由圖11可知，隨著雙星間距離的增加，A為負(fù)值且持續(xù)減小，Xr始終為小于1的正數(shù)，即雙體小行星系統(tǒng)的相對平衡狀態(tài)將在r*增大的過程中始終保持穩(wěn)定狀態(tài)。

圖11 小行星間距離r0對雙體小行星系統(tǒng)平衡態(tài)穩(wěn)定性的影響Fig.11 Influence of the centroid distance on the stability of equilibrium

4 結(jié) 論

本文研究了雙體小行星系統(tǒng)在考慮形狀參數(shù)與質(zhì)量分布的情況下，在相互引力作用運動的平衡態(tài)與穩(wěn)定性問題。文中基于雙橢球構(gòu)成的全二體模型從分析雙體小行星系統(tǒng)的運動可達(dá)范圍的角度給出了確定雙體小行星系統(tǒng)平衡狀態(tài)的穩(wěn)定性的方法。研究發(fā)現(xiàn)小行星形狀參數(shù)β1/α1的增大可使得雙體小行星系統(tǒng)由具有逃逸潛能的不穩(wěn)定平衡態(tài)，經(jīng)歷運動范圍存在上界的平衡狀態(tài)，最終轉(zhuǎn)變?yōu)榉€(wěn)定的平衡狀態(tài)；β2/α2增大將使得雙體小行星系統(tǒng)最終達(dá)到穩(wěn)定狀態(tài)對應(yīng)的β1/α1臨界值更小。雙體小行星系統(tǒng)質(zhì)量分?jǐn)?shù)M1/(M1+M2)的增大對穩(wěn)定性的影響與形狀參數(shù)產(chǎn)生的影響類似，但在M1>>M2時系統(tǒng)由穩(wěn)定平衡狀態(tài)趨近于有上界的不穩(wěn)定狀態(tài)。系統(tǒng)平衡態(tài)時小行星質(zhì)心間的距離r0的改變不會對系統(tǒng)的穩(wěn)定性產(chǎn)生質(zhì)的影響。這些研究結(jié)論對于未來雙體小行星附近探測器的運動研究具有重要的理論指導(dǎo)意義。