基于偶極子模型的雙小行星系統(tǒng)平衡點特性研究

李 旭，錢霙婧，楊曉東，張 偉

(1. 北京工業(yè)大學(xué)材料與制造學(xué)部，北京 100124；2. 北京市非線性振動與機械結(jié)構(gòu)強度重點實驗室，北京 100124)

小行星是太陽系內(nèi)類似行星環(huán)繞太陽運動，但體積和質(zhì)量比行星小得多的天體。19世紀初，意大利天文學(xué)家Piazzi發(fā)現(xiàn)了第一顆小行星1 Ceres。截至2014年已經(jīng)編目的小行星已近50萬顆且在逐年增加，其中包括大量雙小行星系統(tǒng)或多小行星系統(tǒng)[1]。1993年伽利略飛船在前往木星的路途中發(fā)現(xiàn)了Ida和它的衛(wèi)星Dactyl，這是人類第一次發(fā)現(xiàn)的雙小行星系統(tǒng)[2]。隨著天文觀測技術(shù)的不斷提高和航天技術(shù)不斷發(fā)展，比如太陽帆航天器被廣泛應(yīng)用于深空探測任務(wù)[3?4]，被發(fā)現(xiàn)的雙小行星數(shù)量不斷增加，據(jù)估計雙小行星系統(tǒng)約占近地小行星總數(shù)的15%±4%[5]。近年來，雙小行星探測已發(fā)展為未來探測任務(wù)的首選目標之一。歐洲航天局ESA(European Space Agency)提出的MarcoPolo-R任務(wù)中，將雙小行星系統(tǒng)(175706) 1996 FG3的第二主星作為目標小行星進行著陸[6]。美國國家宇航局NASA(National Aeronautics and Space Administration)和歐洲ESA聯(lián)合提出的小行星撞擊和偏轉(zhuǎn)評估任務(wù)(asteroid impact and deflection assessment,AIDA)對存在潛在危險的近地雙小行星(65803)Didymos進行探測器撞擊的測試，以評估驗證近地小行星防御的可能性[7]。美國噴氣推進實驗室(JPL)提出雙小行星(65803) Didymos定位探測(Binary asteroid in-situ explorer mission, BASIX)任務(wù)，擬計劃對雙小行星系統(tǒng)進行定量測量[8]。對雙小行星系統(tǒng)的研究有助于了解太陽系中天體的形成和演化機理，具有獨特科學(xué)價值及經(jīng)濟價值。而深入了解雙小行星系統(tǒng)附近的動力學(xué)環(huán)境，則對軌道設(shè)計和控制策略至關(guān)重要，也是相關(guān)探測任務(wù)成功的關(guān)鍵和基礎(chǔ)。

根據(jù)雙小行星系統(tǒng)的旋轉(zhuǎn)周期和形態(tài)，當系統(tǒng)的旋轉(zhuǎn)周期和兩顆小行星自轉(zhuǎn)周期相同時，稱為雙同步雙小行星系統(tǒng)(doubly synchronous system)。當兩顆小行星的自轉(zhuǎn)周期和系統(tǒng)旋轉(zhuǎn)周期均不同時，稱為異步雙小行星系統(tǒng)(asynchronous system)，此時系統(tǒng)最為復(fù)雜。當一顆小行星自轉(zhuǎn)和系統(tǒng)旋轉(zhuǎn)周期一致時，稱為單同步雙小行星系統(tǒng) (singly synchronous system)[9]，這種情況一般發(fā)生在大尺度比的雙天體和多天體系統(tǒng)中。單同步和雙同步雙星是人類探測任務(wù)最有可能的目標，因為它們有穩(wěn)定的狀態(tài)、形態(tài)特征和豐富的動力學(xué)行為。

由于小行星的不規(guī)則形狀，在其附近運行的探測器所受到的引力場非常復(fù)雜。建立合理的引力場模型是研究這些不規(guī)則小行星附近動力學(xué)的第一步。對于單個小行星的引力場建模問題，經(jīng)歷了從簡單到復(fù)雜、從低維到高維的研究階段。經(jīng)典的引力場建模方法是級數(shù)展開法(series expansion method)[10]。對于近球形的小行星，球諧函數(shù)攝動展開方法具有良好的效果[11]。然而，大部分小行星為強不規(guī)則天體，表面附近Legendre級數(shù)難以收斂[12]。為了改進球諧函數(shù)法，Bierly[13]和MacMillan[14]提出橢球諧函數(shù)法。Romain等[15]對Wirtanen彗星的引力場運用橢球諧函數(shù)法進行建模并對著陸過程進行了模擬。為了進一步精確描述小天體不規(guī)則引力場，1993年Werner等[16]獨立推導(dǎo)了均質(zhì)多面體引力場建模方法。之后，Werner與Scheeres[17]合作改進了上述方法，并成功應(yīng)用到小行星4769 Castalia引力場的建模。此外Geissler等[18]提出質(zhì)點群法來近似小行星整體的引力場。然而，多面體法和質(zhì)點群法的計算時間會隨引入面元或點質(zhì)量的數(shù)量而增加，而且很難分析動力學(xué)特性和模型參數(shù)(例如，形狀尺寸參數(shù)、質(zhì)量參數(shù)、自旋參數(shù)等)之間的關(guān)系，即很難通過參數(shù)變化探究得到某一類具有相似外形小行星引力場動力學(xué)特性。而簡化模型法則能做到這一點，如細直棒模型[19?24]、偶極子模型[25?26]、雙直棒模型[27]、圓環(huán)模型[28?29]、圓餅?zāi)Ｐ蚚30?32]、三角盤和正方形盤[33]、三粒子粒桿模型[34?35]、偶極子段模型[36]、立方體模型[37?39]、啞鈴體模型[40 ? 41]等。

雙小行星系統(tǒng)被認為是由一顆小行星俘獲附近的另一顆小行星或分裂演化而形成的[42?43]。為了簡化雙小行星系統(tǒng)，Wang等[44]采用圓形限制性三體模型(i.e. CRTBP)來描述雙小行星附近的運動，提出一種捕獲策略并將其運用到90 Antiope雙小行星系統(tǒng)和216 Kleopatra小行星系統(tǒng)中。然而，當雙小行星系統(tǒng)中兩天體距離比較近，質(zhì)量分布和形狀強不規(guī)則時，兩天體不能簡單地視為兩質(zhì)量點，CRTBP模型也不再適用。為了準確地研究雙星系統(tǒng)附近的軌道運動，研究者們努力發(fā)展適合描述航天器在雙星系統(tǒng)重力場中的軌道運動的模型和方法。與單個小行星相似，雙小行星系統(tǒng)中的天體也可以被合適地視為球體、橢球體或者不規(guī)則體。一些簡化模型法被成功應(yīng)用到雙小行星系統(tǒng)中。Bellerose和Scheeres[45?46]研究了球—橢球雙星系統(tǒng)的周期軌道，特別是平動點周期軌道。Chappaz和Howell[47]研究了球—橢球雙星系統(tǒng)和橢球—橢球雙星系統(tǒng)。Linder等[21]用細直棒和一個點質(zhì)量來模擬雙小行星系統(tǒng)中的主星和次星，并發(fā)現(xiàn)細直棒附近存在穩(wěn)定的同步軌道、一般的混沌軌道、不穩(wěn)定周期軌道、旋轉(zhuǎn)穩(wěn)定軌道等。Blaikie等[48]研究了2個細直棒組成的雙小行星系統(tǒng)在相互影響下各自的軌道與姿態(tài)響應(yīng)。Jain和Sinha[49?50]考慮兩個細直棒在同步旋轉(zhuǎn)狀態(tài)下系統(tǒng)平衡點附近的動力學(xué)行為。Ferrari等[51]使用3個粒子模擬雙小行星，用2個不同的三體系統(tǒng)構(gòu)建模型，他們稱該模型為拼接三體問題(patched three-body problem)。類似地，Torres Dos Santos等[52]將一個主小行星視為球體、將另一個不規(guī)則小行星視為旋轉(zhuǎn)的質(zhì)量偶極子去構(gòu)建同步雙小行星系統(tǒng)模型，并對系統(tǒng)的平衡點進行了分析。除了研究平動點周期軌道，雙小行星系統(tǒng)的共振軌道也被研究和分析。楊雅迪等[53]基于雙橢球模型，給出同步雙小行星系統(tǒng)中共振軌道的設(shè)計方法并以雙小行星系統(tǒng)1999 KW4為例設(shè)計了1∶1、1∶2、1∶3、1∶4、2∶3共振軌道族。Nadoushan和Assadian[54]假設(shè)2個小行星為一個球體和一個橢球體，研究了自旋-軌道共振，同時還計算了一個假定的雙小行星系統(tǒng)的龐加萊截面。Shang等[55]基于軌道的對稱性搜索了橢球—橢球雙星系統(tǒng)的全局周期軌道，得到了不同的共振軌道族。

小行星平衡點的位置為進行長期科學(xué)觀測提供了難得的有利條件。利用平衡點處的動力學(xué)特性，可以設(shè)計很多具有實際應(yīng)用價值的軌道，對擴展探測任務(wù)的科學(xué)回報具有重要意義。一般國內(nèi)外學(xué)者都是重點關(guān)注小行星自身平衡點的存在性、個數(shù)與穩(wěn)定性問題。Elipe和Riaguas[56]計算了旋轉(zhuǎn)對數(shù)函數(shù)并找到有限長細直棒引力場中的4個外部平衡點，同時討論了這4個平衡點的線性穩(wěn)定性。Scheeres等[57]發(fā)現(xiàn)了小行星25143 Itokawa引力場中的4個外部平衡點及相應(yīng)的位置坐標。Mondelo等[58]計算了小行星4 Vesta引力場中的4個外部平衡點的坐標及穩(wěn)定性；他們的研究表明，其中2個平衡點是穩(wěn)定的，另外2個是不穩(wěn)定的。Liu等[37]給出了當引力加速度與離心加速度之比為1時的旋轉(zhuǎn)均質(zhì)立方體引力場中的平衡點位置。Yu和Baoyin[59]基于多面體模型分析了小行星216 Kleopatra引力場中的4個外部平衡點及其局部流形，并給出了這些平衡點的位置坐標、特征值與穩(wěn)定性。Scheeres[60]發(fā)現(xiàn)了小行星1580 Betulia引力場中的6個外部平衡點，以及彗星67P/CG引力場中的4個外部平衡點。Jiang等[61]建立了一般旋轉(zhuǎn)不規(guī)則小天體引力場中平衡點附近的運動理論，根據(jù)特征方程的特征值的不同情況對小行星的平衡點進行了拓撲分類。之后Jiang[62]還討論了小行星平衡點附近的周期軌道和平衡點的拓撲類型關(guān)系。Wang、Jiang和Gong[63]計算了23個強不規(guī)則小天體的平衡點的個數(shù)、分布情況及穩(wěn)定性。Jiang等[64]還研究了22 Kalliope小行星的平衡點在自轉(zhuǎn)速率和密度的影響下的運動。在動力學(xué)研究領(lǐng)域廣泛應(yīng)用的參數(shù)化分析方法[65?66]也適用于小行星平衡點相關(guān)特性的研究，比如Zeng等[36]建立偶極子段模型，參數(shù)化研究了系統(tǒng)平衡點的位置和穩(wěn)定性變化趨勢。這些工作對本文的平衡點動力學(xué)特性的研究具有一定的指導(dǎo)意義。

由文獻可知，簡化模型適合前期的定性分析，幫助理解或認識實際不規(guī)則的雙小行星系統(tǒng)特性，且具有易于推導(dǎo)的引力勢表達式，容易展開動力學(xué)特性分析。偶極子模型作為一種簡化模型，其能夠有效地近似細長類小行星引力場。文章針對主星是細長型的小行星，而次星是小而規(guī)則的天體的雙小行星系統(tǒng)，建立了普適性的引力場模型。針對同步旋轉(zhuǎn)情況，對系統(tǒng)的平動點特性進行研究，細致分析系統(tǒng)參數(shù)對同步雙小行星系統(tǒng)平衡點的影響。在非同步情況下，給出了系統(tǒng)等效平衡點的運動規(guī)律。最后，結(jié)合路徑搜索修正法與偽弧長延拓方法有效地計算得到共線平衡點附近的周期軌道族。本文的研究可以為未來這一類雙小行星系統(tǒng)的探測提供參考。

1 雙小行星系統(tǒng)動力學(xué)模型

本文關(guān)注主星是細長型的小行星，而次星是小而規(guī)則的天體的雙小行星系統(tǒng)，并采用偶極子—粒子模型對系統(tǒng)進行描述，如圖1所示。系統(tǒng)中，主小行星P12為細長型強不規(guī)則小行星，由2個偶極子P1和P2(質(zhì)量分別是m1和m2)和一根長度固定為ld的無質(zhì)量桿細桿構(gòu)成[33]。主小行星P12繞其質(zhì)心B勻速自轉(zhuǎn)運動。次小行星P4為質(zhì)點，其質(zhì)量為m4，且m4＜ m1+ m2。質(zhì)量點P1、P2和P4在同一平面內(nèi)運動。主小行星P12質(zhì)心和次小行星P4之間距離為h，且h＞ld。P3為探測器，其質(zhì)量為m3，探測器運動對其他引力體不產(chǎn)生影響。在系統(tǒng)中，P1和P2到質(zhì)心B的距離分別為d1和d2，且有l(wèi)d= d1+ d2。P4到質(zhì)心O的距離為s1，而主小行星質(zhì)心B到系統(tǒng)質(zhì)心O之間的距離為s2，h = s1+ s2。

圖1 雙小行星系統(tǒng)簡化模型示意圖Fig. 1 The simplified model of binary asteroid

為方便描述，文中定義如下3個坐標系統(tǒng)：

1) 雙小行星系統(tǒng)慣性坐標系(O-XYZ)：原點O是雙小行星系統(tǒng)的質(zhì)心；X-Y坐標面即兩個小行星相對運動平面，X軸方向的選擇，對應(yīng)初始時刻t=t0時，兩個小行星處于該坐標軸上，且正方向由系統(tǒng)質(zhì)心O指向主小行星質(zhì)心B，Z軸垂直于雙小行星相對運動平面且和系統(tǒng)角動量方向一致；OY軸與另外兩軸形成右手坐標系。

2) 系統(tǒng)旋轉(zhuǎn)坐標系(B-XrYrZr)：原點B位于主小行星質(zhì)心，該坐標系相對于O-XYZ坐標系的角速度為Ω1=φ˙Z，即兩個小行星圍繞系統(tǒng)質(zhì)心系的運動角速度。Xr軸正方向由P4指向主小行星質(zhì)心B。BZr軸和OZ軸平行且正方向一致，BYr軸與另外兩軸形成右手坐標系。

3) 主小行星體坐標系(B-xyz)：原點B位于主小行星質(zhì)心，該坐標系相對于O-XYZ坐標系的絕對自轉(zhuǎn)角速度為Ω2=γ˙Z，Bx軸與偶極子段重合沿P1指向P2方向；Bz軸為自轉(zhuǎn)軸且其正方向為主小行星自轉(zhuǎn)角動量方向；By軸與另外兩軸形成右手坐標系。此外，主小行星體坐標系(B-xyz)以角速度Ω=θ˙Zr相對于坐標系(B-XrYrZr)做勻速旋轉(zhuǎn)運動。由于Zr和Z方向一致，Ω亦可表示為Ω=Z。根據(jù)模型描述易知，Ω2=Ω1+Ω。

在慣性系中，則探測器P3的加速度可表示為：

式中：Ri=‖Ri‖(i=1,2,3)表示探測器P3相對Pi的距離；ρ為探測器P3的位置矢量，可表達為：

式中，(x,y,z)是探測器P3在主小行星體坐標系中位置矢量R的坐標。對式(2)進行求導(dǎo)可得：

式中：上標“I”“rel”分別表示相對雙小行星系統(tǒng)慣性坐標系和主小行星體坐標系進行求導(dǎo)。將對應(yīng)的角速度Ω1、Ω2代入式(3)并求導(dǎo)，可得探測器P3的速度和加速度：

接下來，對于式(1)右邊部分進行計算，將Ri(i=1,2,3)表示為：

將式(6)代入式(1)右邊部分，式(5)代入式(1)左邊部分，兩邊對比，取x前面的系數(shù)，可以得到：

為了分析問題和計算上的方便，在求動力學(xué)方程時，往往采用無量綱形式，本文中將各物理量相應(yīng)的長度單位[L]、質(zhì)量單位[M]和時間單位[T]分別取為：

由此，在無量綱系統(tǒng)中，主小行星相對于慣性系的旋轉(zhuǎn)角速度為1，且 P1、P2和P4的質(zhì)量有：

通過主小行星質(zhì)心和雙小行星系統(tǒng)質(zhì)心可知，在無量綱系統(tǒng)中：

并得到：

式中，上標“^”表示對應(yīng)的無量綱量。用ν表示在無量綱系統(tǒng)中雙小行星系統(tǒng)相對慣性系的旋轉(zhuǎn)角速度，則：

式中，k為“受力比”，其反映了小行星的自身特性，具有如下的形式：

式中，G是引力常數(shù)。

由此，對式(7)進行無量綱化，得到：

同理，可以得到另外兩個方向的方程。此外，由雙小行星系統(tǒng)和主小行星均為勻速自轉(zhuǎn)可知：

代入后得到探測器P3的動力學(xué)方程。為方便后續(xù)研究工作，用不帶“^”的符號表示無量綱物理量，則各分量形式如下：

式中：θ=(1?v)×t+θ0，t為無量綱意義下的時間，θ0為主小行星的初始相位，即在初始時刻，正向Xr軸與正向x軸之間的夾角。

式(17)即為本文關(guān)注的雙小行星系統(tǒng)普適性引力場模型。根據(jù)v的定義，當v=1時，式(17)表示同步雙小行星系統(tǒng)的動力學(xué)方程；當0＜v＜1時，式(17)表示非同步雙小行星系統(tǒng)的動力學(xué)方程。

2 同步雙小行星系統(tǒng)平衡點及其穩(wěn)定性

當雙小行星系統(tǒng)同步運轉(zhuǎn)時，即雙小行星系統(tǒng)旋轉(zhuǎn)的角速度和主小行星自轉(zhuǎn)角速度一致時，ν=1。此時動力學(xué)方程為：

式中，V為有效引力勢，其表達式如下：

通過式(18)可以得到一個能量積分：

式中，C為雅克比常數(shù)，決定了探測器被允許的運動區(qū)域。

應(yīng)該注意的是，雖然式(18)是建立在B-xyz坐標系下的動力學(xué)方程，但是對于同步雙小行星系統(tǒng)，由式(18)計算得到的平衡點仍然是整個系統(tǒng)的平衡點。此外，由上述分析可知受力比k、初始相位θ0、質(zhì)量參數(shù)μ1和μ4對平衡點的分布及穩(wěn)定性都有著重要影響。本章將研究系統(tǒng)參數(shù)對于同步雙小行星系統(tǒng)平衡點分布及穩(wěn)定性的影響。

2.1 平衡點分布特性研究

對于同步雙小行星系統(tǒng)，給定系統(tǒng)的初始相位后，系統(tǒng)的平衡點是固定不變的。由于偶極子—粒子模型關(guān)于B-xy對稱，所以系統(tǒng)平衡點都位于該平面內(nèi)，非共線平動解滿足：

而共線平衡點則滿足：

進而通過牛頓迭代法可以得到平衡點的精確位置。

圖2是系統(tǒng)參數(shù)[μ1，μ4，k]T=[0.7，0.05，9]T且初始相位為0時，系統(tǒng)零速度面等高線以及平衡點的相對位置情況(未考慮小行星內(nèi)部的平衡點)，其中橫軸和縱軸是無量綱意義下的x軸和y軸。該參數(shù)條件下平衡點的具體坐標列在表1中。在本文中，坐標系統(tǒng)中的x和 y都代表無量綱坐標。

圖2 μ1=0.7，μ4=0.05，k = 9時的平衡點以及零速度面Fig. 2 Equilibrium points and zero-velocity curves whenμ1=0.7, μ4=0.05 and k = 9

表1 圖2中平衡點位置Table1 The positions of the equilibrium points in Figure 2

表2進一步給出了平衡點位置隨系統(tǒng)參數(shù)[μ1，μ4，k]T變化時的移動規(guī)律。表中：斜箭頭符號“↗”表示參數(shù)增大；符號“≡”表示參數(shù)值固定不變；4個水平方向箭頭表示平衡點坐標的移動方向。例如，當受力比k和μ1固定不變時，共線平衡點E2隨著μ4的增加向左移動遠離主小行星，E1向右移動靠近主小行星，E3向左移動靠近主小行星。非共線平衡點E4的x分量向左移動，y分量向上移動遠離雙小行星系統(tǒng)，E5與E4始終保持對稱。當μ1和μ4固定不變時，共線平衡點E1和E2隨著受力比k的增加向左移動，而共線平衡點E3向右移動，非共線平衡點E4的y分量向上移動，總之，平衡點隨著受力比k的增加而全部遠離雙小行星系統(tǒng)質(zhì)心。當k和μ4固定不變時，所有平衡點隨著μ1的增加變化趨勢都會出現(xiàn)轉(zhuǎn)折，比如，共線平衡點E3隨著μ1的增加先向右移動再向左移動，即共線平衡點E3先遠離系統(tǒng)質(zhì)心再靠近系統(tǒng)質(zhì)心。

表2 系統(tǒng)平衡點隨系統(tǒng)參數(shù)變化的移動規(guī)律Table2 Variation trends of coordinates for the equilibriumpoints with the change of system parameters

圖3給出μ4=0.04時，不同μ1下，非共線平衡點E4位置隨k值從9.025減小到0.125的變化情況。由圖3可知，對于不同的μ1，平衡點的變化趨勢都是隨著k減小，x坐標逐漸增加，y坐標逐漸減小，而且μ1值越小，x變化量越大。由于雙小行星系統(tǒng)中主次小行星之間的距離應(yīng)該大于主星質(zhì)心到兩個偶極子之間的距離，否則小行星之間會產(chǎn)生碰撞。圖中“×”表示不符合文中模型假設(shè)的平衡點，即使得h＜μ1或h＜1?μ1的平衡點。

圖3 μ4=0.04，不同μ1下，非共線平衡點E4隨k從9.025減小到0.125的位置變化情況Fig. 3 When μ4=0.04, different locations of E4 with variations of μ1 and k

圖4中紅點表示k=0.1時不同μ1下，非共線平衡點E4隨著μ4從0.5減小到0 的位置變化情況。由圖4可知，對于不同的μ1，平衡點的變化趨勢都是隨著μ4減小，x坐標逐漸增加，y坐標逐漸減小，而且μ1值越小，x變化量越大。當μ1≤0.5且μ4≠0時，由于系統(tǒng)中次星的存在時，不同于單偶極子模型描述的系統(tǒng)[67]，雙小行星系統(tǒng)的非共線平衡點E4并沒有消失，依然存在。而當μ4=0時，由于系統(tǒng)退化為單偶極子模型，系統(tǒng)的拓撲構(gòu)型發(fā)生變化，無論μ1如何取值，非共線平衡點E4(或E5)均會消失，而產(chǎn)生新的共線平衡點。為了討論一般性，兩個極限情況k→0和k→∞沒有在本文進行探究。此外，當μ1＞0.5且μ4≠0時，計算所得的平衡點均不符合文中模型的描述，即h＜μ1或h＜1?μ1，用“×”表示。

圖4 k=0.1，不同μ1下，非共線平衡點E4隨μ4從0.5減小到0的位置變化情況Fig. 4 When k=0.1, different locations of E4 with variations of μ1 and μ4

2.2 平衡點穩(wěn)定性

同步雙小行星系統(tǒng)式(18)線性部分的特征方程為：

特征方程(23)的根決定了平衡點周圍的運動特性。如果特征根具有純虛根或具有負實部的復(fù)根，則平衡點是穩(wěn)定的。具體的判斷條件如下：

基于以上穩(wěn)定條件，通過牛頓迭代計算發(fā)現(xiàn)，共線平衡點是不穩(wěn)定的，而非共線平衡點存在滿足線性穩(wěn)定性的參數(shù)域。由于多數(shù)雙小行星系統(tǒng)中，次小行星的質(zhì)量比主小行星的質(zhì)量小很多，即μ4較小。因此仿真中，μ4取值為0.005～0.040，μ1取值為0～1，k則取值1～100，穩(wěn)定邊界曲線結(jié)果如圖5(a)所示。當系統(tǒng)參數(shù)取值位于該曲面上方區(qū)域時，系統(tǒng)具有穩(wěn)定的非共線平衡點；反之則存在不穩(wěn)定的非共線平衡點。圖5(b)和圖5(c)分別為圖5(a)在μ1-k平面和μ4-k平面的投影。

從圖5(b)中可以看到，當μ4比較大時，穩(wěn)定邊界曲線峰值較高，對應(yīng)k值較大，且曲線向右偏斜；而當μ4較小時，穩(wěn)定邊界曲線峰值下降，穩(wěn)定域變大，對應(yīng)k值較小。當μ4=0.005時，其曲線基本關(guān)于μ1=0.5對稱，非常接近單偶極子模型的穩(wěn)定域情況[67]。由此可知，次小行星質(zhì)量越小，系統(tǒng)平衡點的穩(wěn)定邊界所對應(yīng)的k值越小。對于特例μ1=0或者μ1=1(兩者所對應(yīng)的系統(tǒng)一致)，此時雙小行星系統(tǒng)的主小行星已經(jīng)簡化為一球體，系統(tǒng)退化為限制性三體問題系統(tǒng)，穩(wěn)定域內(nèi)k的最小值為1，與限制性三體問題結(jié)論一致。

當μ1固定不變時，非共線平衡點E4穩(wěn)定性隨參數(shù)k和μ4變化規(guī)律如圖5(c)所示。各曲線下方為不穩(wěn)定區(qū)域，上方為穩(wěn)定區(qū)域。由圖5(c)可知，當μ1的值保持不變并不斷增大μ4時，受力比k的臨界值會不斷增大。特別的，當μ1≤0.2或μ1=1時，k的臨界值幾乎為一常值，使得平衡點的不穩(wěn)定區(qū)域近似為一個矩形。這表明，同步雙小行星系統(tǒng)在μ1≤0.2或μ1=1后，次星對系統(tǒng)非共線平衡點穩(wěn)定性的影響已非常小，幾乎可以忽略不記。

圖5 非共線平衡點的線性穩(wěn)定域Fig. 5 Linear stable area for non-collinear equilibrium points

另外發(fā)現(xiàn)隨著μ4逐漸增大，平衡點穩(wěn)定性的規(guī)律會發(fā)生變化。當μ4=0.05時，非共線平衡點不再具有以上規(guī)律，只是μ1和k在很小的時候平衡點具有穩(wěn)定性，如圖6所示。圖中陰影區(qū)域是非共線平衡點的穩(wěn)定域。

圖6 μ4=0.05時非共線平衡點的線性穩(wěn)定域Fig. 6 Linear stable area for non-collinear equilibrium points when μ4=0.05

綜上所述，偶極子—粒子模型在系統(tǒng)質(zhì)量固定且次星質(zhì)量μ4較小的情況下，受力比k越大，非共線平衡點越穩(wěn)定。即，主小行星系統(tǒng)質(zhì)量固定且桿長度不變的情況下，主小行星自轉(zhuǎn)角速度越小，平衡點越趨于穩(wěn)定。

2.3 初始相位對平衡點的影響

由式可知，有效勢能V與初始相位θ0相關(guān)。不同的初始相位會造成系統(tǒng)的對稱性發(fā)生變化，使得系統(tǒng)平衡點的個數(shù)及位置出現(xiàn)變化。

在 本 小 節(jié) 以 系 統(tǒng) 參 數(shù) [μ1，μ4，k]T=[0.7，0.05，9]T為例來說明初始相位對平衡點的影響。對于不同的初始相位，系統(tǒng)平衡點的個數(shù)和位置可以通過零速度曲線進行判斷，然后再用牛頓迭代法計算出平衡點的精確值，見附錄A。由于系統(tǒng)關(guān)于原點對稱，所以文中考慮0°～180°相位變化如圖7(a)～圖7(f)所示。

從圖7(a)～圖7(f)可知，隨著初始相位變化，所有平衡點的位置均發(fā)生變化。在初始相位從0°～90°變化過程中，次小行星的位置發(fā)生變化，它會“擠壓”非共線平衡點E4，同時會“吸引”另一個非共線平衡點E5。而當次小行星到達初始相位為0°時平衡點E4的周圍時，在“吸引”過程中將“拉扯”出2個新的平衡點E6和E7，平衡點的個數(shù)由5變?yōu)?。在初始相位從90°～180°變化過程中，次小行星會將繼續(xù)“擠壓”平衡點E4，最終使得平衡點E3和E4“湮滅”，平衡點的個數(shù)減少為5。

圖7 平衡點和初始相位從 0°～180°變化雙小行星系統(tǒng)在x-y平面的零速度曲線Fig. 7 Equilibrium points and the zero-velocity curves in the equatorial plane when the initial phase changes from 0°～180°

圖7給出了平衡點隨著相位θ0的變化變化情況，此時θ0表示二元系統(tǒng)的相對構(gòu)型。而自然系4 2 0統(tǒng)在能量上傾向于繞其最大慣性軸旋轉(zhuǎn)，這一構(gòu)型有能量穩(wěn)定的解[68?69]，所以圖7(a)和圖7(g)(初始相位為0°或180°)代表的雙小行星系統(tǒng)是穩(wěn)定的。

3 非同步雙小行星系統(tǒng)的等效平衡點運動軌跡

當0＜v＜1時，式(17)表示非同步雙小行星系統(tǒng)的動力學(xué)方程。令方程右側(cè)等于0能夠得到非同步情況下的等效平衡點。假設(shè)系統(tǒng)參數(shù) [μ1，μ4，k，v]T=[0.7，0.05，9，0.1]T，并設(shè)θ0=0。初始時刻t=0時，各等效平衡點位置情況如圖8所示。值得注意的是，為了更清楚地表示等高線和平衡點位置的關(guān)系，次小行星外側(cè)存在一等效平衡點E5并未在圖中表示。

圖8 初始時刻非同步雙小行星系統(tǒng)的零速度曲線Fig. 8 Equilibrium points in the non-synchronous state and the zero-velocity curves at the initial moment

在一個θ的周期內(nèi)，圖9、圖10、圖11～圖13分別展示了這5個等效平衡點在x-y平面的軌跡，其中黑點表示在初始時刻等效平衡點的位置?？梢钥闯?，在該模型的假設(shè)下，非同步雙小行星系統(tǒng)的等效平衡點的軌跡是周期性的。其中等效平衡點E1和E2、等效平衡點E3和E4的軌跡變化規(guī)律相似，等效平衡點E1和E2在x方向上的變化小，等效平衡點E3和E4在y方向上的變化小。等效平衡點E5的軌跡呈近圓形軌跡。

圖9 等效平衡點E1的軌跡Fig. 9 Trajectory of equivalent equilibrium point E1

圖10 等效平衡點E2的軌跡Fig. 10 Trajectory of equivalent equilibrium point E2

圖11 等效平衡點E3的軌跡Fig. 11 Trajectory of equivalent equilibrium point E3

圖13 等效平衡點E5的軌跡Fig. 13 Trajectory of equivalent equilibrium point E5

4 同步雙小行星系統(tǒng)共線平衡點附近的軌道族

同步雙小行星系統(tǒng)的共線平衡點具有特殊的空間位置和動力學(xué)特征。利用共線平衡點在旋轉(zhuǎn)坐標系中位置固定的特性，可以實現(xiàn)雙小行星系統(tǒng)的觀測任務(wù)。例如，在圍繞平衡點E1的Lyapunov軌道上，探測器能夠同時實現(xiàn)觀測兩個小行星的任務(wù)。然而，文中偶極子—粒子模型不同于限制性三體模型，無法采用解析法求解系統(tǒng)的三階解析解。因此本節(jié)將基于平動點周期軌道垂直穿越龐加萊截面的幾何對稱性，利用路徑搜索修正法探索周期軌道[70?71]。再使用偽弧長延拓法尋找同步雙小行星系統(tǒng)共線平衡點附近的周期軌道族。

圖12 等效平衡點E4的軌跡Fig. 12 Trajectory of equivalent equilibrium point E4

4.1 同步雙小行星系統(tǒng)的周期軌道

在x軸上取一點[x0, 0]，作為平動點周期軌道搜索起點，軌道的初值速度設(shè)置為[0,]。雅可比常量前文已給出，那么對于初始速度，有：

式中，V是起點[x0, 0]處的有效勢能。

通過改變C的大小，能夠得到一系列起點在[x0, 0]的軌道。記錄穿越x軸時的速度和對應(yīng)的雅可比常數(shù)[C]。在-C曲線中，0對應(yīng)的初始狀態(tài)是周期軌道初值的理想猜測值。

假定系統(tǒng)參數(shù)為[μ1, μ4, k]T=[0.7, 0.05, 9]T，初始相位為180°。本節(jié)將利用上述方法來搜索E1、E2和E3附近的周期軌道。在共線平衡點E1附近，設(shè)置搜索起點x0=1.5。圖14是雅克比常數(shù)C與軌道首次穿越x軸時速度分量的關(guān)系。從圖中可以看出，曲線與=0交于A、B兩點，運用單步打靶法可求得A和B處的軌道精確初始狀態(tài)為[1.500 954 738 784 62;0;0;0; 1.981 369 681 868 90;0]T和[1.500 496 445 354 81;0;0;0;1.683 561 820 199 68;0]T，對應(yīng)周期分別是3.700 567 127 127 2和3.856 824 801 326 6。

圖14 起點為x0=1.5時雅克比常數(shù)C與速度的關(guān)系曲線Fig. 14 The relationship between Jacobian constant C and velocity component when x0=1.5

采用圖14中A、B處的初值，通過積分可得到圖15中的軌道。圖中實線代表B處軌道，虛線代表A處軌道。其中，軌道A是圍繞次小行星的逆行周期軌道，而軌道B是圍繞共線平衡點E1的Lyapunov軌道。

圖15 A、B處的周期軌道Fig. 15 The periodic orbits in A and B

類似的，在共線平衡點E2附近，利用上述方法可得到雅克比常數(shù)C與從x0=2.5出發(fā)的軌道首次穿越x軸時速度分量的關(guān)系曲線，如圖16所示。圖中，D、E處均可能存在Lyapunov軌道。利用單步打靶法可求得精確初始狀態(tài)分別為[2.500 303 091 710 3;0;0;0;0.904 595 587 122 23;0]T和[2.502 000 755 613 91;0;0;0;0.692 197 762 539 07;0]T，且相應(yīng)周期分別是3.750 056 994 207 15和2.488 268 157 419 96。

圖16 起點為x0=2.5時雅克比常數(shù)C與速度的關(guān)系曲線Fig. 16 The relationship between Jacobian constant C and velocity component when x0=2.5

采用圖16中D、E處的初值，通過積分可得到圖17中的軌道。圖中，實線代表D處軌道，虛線代表E處軌道。其中，軌道D是圍繞平衡點E2的平面Lyapunov軌道，軌道E則是圍繞次小行星P4的順行軌道。

圖17 D、E處的周期軌道Fig. 17 The periodic orbits in D and E

在共線平衡點E3附近，利用上述方法可得到雅克比常數(shù)C與從x0=?2.5出發(fā)的軌道首次穿越x軸時速度分量的關(guān)系曲線，如圖18所示。圖中，F(xiàn)、G處均可能存在Lyapunov軌道。由此，求得周期軌道的精確初始狀態(tài)分別為[?2.501 175 834 544 529;0;0;0;2.300 538 550 025 692;0]T和[?2.498 584 107 628 348;0;0;0;0.757 759 3649 790 15;0]T，且相應(yīng)周期分別是2.960 621 259 436 86和5.373 663 202 852 81。

圖18 起點為x0=?2.5時雅克比常數(shù)C與速度的關(guān)系曲線Fig. 18 The relationship between Jacobian constant C and velocity component when x0=?2.5

采用圖18中F、G處的初值，通過積分可得到圖19中的軌道。圖中，實線代表G處軌道，虛線代表F處軌道。其中，軌道G是圍繞平衡點E3的平面Lyapunov軌道。值得注意的是，因為軌道F與主小行星會發(fā)生碰撞，所以軌道F并不存在。

圖19 F、G處的周期軌道Fig. 19 The periodic orbits in F and G

4.2 共線平衡點附近的軌道族

通過路徑搜索法發(fā)現(xiàn)了雙小行星系統(tǒng)在平衡點E1、E2和E3周圍存在平面周期軌道，再采用偽弧長延拓法則可以找到周期軌道所對應(yīng)的軌道族[72]。具體步驟如下：

1) 從路徑搜索方法中可以得到周期軌道的一個初值猜測，修正得到第一個初值狀態(tài)。略微改變搜索起始點，能得到與第一個初始狀態(tài)非常接近的第二個初值狀態(tài)。

2) 定義狀態(tài)參數(shù)Z。

對于平面問題，這5個參數(shù)確定了周期軌道的解。根據(jù)第1)步可以得到在第(j?2)條軌道和第(j?1)條軌道的準確周期解(j ≥ 3)，并表示為：

由此，選取在第(j?1)條軌道處的演化方向為：

3) 根據(jù)在第(j?1)條軌道處的演化方向和準確初值，得到第j條軌道的預(yù)測初值：

式中，Δs是很小的預(yù)測步長。

4) 通過單步打靶法修正得到第j條軌道準確初值狀態(tài)。再重復(fù)以上步驟。

根據(jù)上述數(shù)值延拓方法對4.1節(jié)中B、D、G處的Lyapunov軌道進行延拓，如圖20所示，其中：實線軌道代表圍繞E1的軌道，稱為軌道族1；虛線軌道代表圍繞E2的軌道，稱為軌道族2；點劃線軌道代表圍繞E3的軌道，稱為軌道族3。軌道族1、軌道族2和軌道族3的初始狀態(tài)，周期以及雅克比常數(shù)的情況如圖21所示。其中x軸、y軸和z軸分別代表軌道初始時刻的位置、速度和雅克比常數(shù)，顏色代表軌道的周期。平衡點E1、E2和E3的雅克比常數(shù)分別是C1= 15.954 356 339 033 744，C2=15.120 681 652 208 280，C3=14.130 225 232 500 724。

圖21 平面Lyapunov軌道族初始位置、周期以及雅克比常數(shù)關(guān)系圖Fig. 21 Initial state, period and Jacobi constant of the planar Lyapunov families

不難發(fā)現(xiàn)，這些軌道都是與系統(tǒng)自轉(zhuǎn)角速度φ˙ 、 主星自轉(zhuǎn)角速度接近1∶1共振的軌道族。共振軌道越接近平衡點，軌道周期越小，軌道的雅克比常數(shù)越接近平衡點處的雅克比常數(shù)。

5 結(jié)論

本文采用偶極子—粒子模型研究了主星為細長類小行星的雙小行星系統(tǒng)平衡點附近的動力學(xué)特性和周期軌道。

文章采用偶極子模型和質(zhì)量點分別描述細長類主小行星和次小行星，建立了雙小行星系統(tǒng)普適性引力場模型。當系統(tǒng)旋轉(zhuǎn)參數(shù)v=1時，表示同步雙小行星系統(tǒng)；當0＜v＜1時，表示非同步雙小行星系統(tǒng)。其次，文章重點討論了同步雙小行星系統(tǒng)的平衡點的分布特征、存在性及穩(wěn)定性。通過分析可知，在系統(tǒng)質(zhì)量固定且次星質(zhì)量μ4較小的情況下，受力比k越大，非共線平衡點越穩(wěn)定。即，主小行星系統(tǒng)質(zhì)量固定且桿長度不變的情況下，主小行星自轉(zhuǎn)角速度越小，平衡點越趨于穩(wěn)定。此外，發(fā)現(xiàn)隨著初始相位的變化，系統(tǒng)平衡點將隨之出現(xiàn)或湮滅。對于非同步情況，給出一個示例并得到了系統(tǒng)等效平衡點的運動軌跡。最后，采用路徑搜索法和偽弧長延拓法得到了同步雙小行星系統(tǒng)共線平衡點附近的3類1∶1共振軌道族。

本文給出的普適性引力場模型、平衡點的分布和穩(wěn)定性規(guī)律分析結(jié)論以及共振軌道族設(shè)計方法為雙小行星系統(tǒng)探測任務(wù)提供了理論基礎(chǔ)。