楊雅迪,陳 奇,李翔宇,喬 棟(1. 北京理工大學(xué)宇航學(xué)院,北京 100081;2. 深空探測自主導(dǎo)航與控制工信部重點實驗室,北京 100081)
自20世紀90年代以來,小行星探測成為深空探測的熱點。在小行星探測任務(wù)中,具有獨特動力學(xué)特性的雙小行星系統(tǒng)引起了科學(xué)家的關(guān)注[1]。通過對雙小行星系統(tǒng)展開實地探測與采樣,可獲取小行星內(nèi)部構(gòu)造、物質(zhì)組成和動力學(xué)演化等方面的重要信息,能得到對單顆小行星探測時無法獲得的獨特科學(xué)回報。
目前各國提出了多項針對雙小行星系統(tǒng)的探測計劃。美國NASA與歐洲ESA聯(lián)合提出了小行星撞擊和偏轉(zhuǎn)評估(AIDA)任務(wù)[2],計劃對雙小行星65803 Didymos進行撞擊,評估驗證近地小行星防御的可行性;美國JPL同時提出雙小行星定位探測(BASiX)任務(wù)[3],以雙小行星65803 Didymos為目標,研究微重力環(huán)境中瓦礫堆小行星的起源及演化。在雙小行星系統(tǒng)的探測任務(wù)中,復(fù)雜探測任務(wù)的開展首先需要獲得系統(tǒng)動力學(xué)環(huán)境相關(guān)的先驗信息,而系統(tǒng)完整的引力場模型必須在探測器靠近探測目標,并展開長時間觀測與繞飛后才能獲得[4-5]。同時通過環(huán)繞探測,也將為著陸點選擇和著陸軌跡設(shè)計提供參考[6-7]。因此選擇合適的環(huán)繞探測軌道對探測任務(wù)的成功開展具有舉足輕重的作用。為了滿足小行星附近尤其是雙小行星系統(tǒng)中的探測任務(wù)需求,需要多樣性豐富、燃料消耗少且可保持長期穩(wěn)定的探測軌道。共振軌道可以較好地滿足這幾大特性,因此研究小行星附近的共振軌道對未來探測任務(wù)具有較大的參考價值。
在對共振軌道的研究中,2013和2014年,Vaquero和Howell[8-9]對地月系統(tǒng)共振軌道的動力學(xué)結(jié)構(gòu)進行了深入細致的分析,借助龐加萊截面設(shè)計了二維及三維共振軌道,并提出利用共振軌道實現(xiàn)低能量轉(zhuǎn)移的軌道設(shè)計方法。2014年,Antoniadou和Voyatzis[10]研究了三體系統(tǒng)下地外行星系統(tǒng)平面及空間中共振比為4∶3,3∶2,5∶2,3∶1及4∶1的共振軌道,研究發(fā)現(xiàn)穩(wěn)定共振軌道普遍具有較大的偏心率。同年,于洋[11]研究了單顆小行星赤道面附近的1∶1共振軌道,從能量角度揭示了赤道面附近的1∶1共振的動力學(xué)本質(zhì),確定了形成1∶1共振拋射軌道的臨界條件。2016年,Anderson等[12]基于圓型限制性三體模型,利用網(wǎng)格搜索方法對木衛(wèi)系統(tǒng)中不穩(wěn)定共振軌道進行了搜索,并對不同共振比軌道進行了特性分析。研究發(fā)現(xiàn)在穩(wěn)定共振軌道的基礎(chǔ)上,可利用連續(xù)法計算得到不穩(wěn)定的共振軌道。此后,F(xiàn)eng等[13]對接觸式雙小行星附近動力學(xué)環(huán)境進行了分析,發(fā)現(xiàn)了三類不同共振比的共振軌道。
此前對共振軌道的研究,主要集中在地月和木衛(wèi)系統(tǒng)[14]等質(zhì)量比較大的系統(tǒng)。雙小行星系統(tǒng)具有質(zhì)量比較小,形狀不規(guī)則和相對距離接近等特點,導(dǎo)致其附近動力學(xué)特性復(fù)雜。而針對雙小行星系統(tǒng)內(nèi)的共振軌道尚未得到系統(tǒng)性研究。因此本文將以雙小行星系統(tǒng)探測任務(wù)為背景,給出一種具有普適性的雙小行星系統(tǒng)共振軌道設(shè)計方法,對不同共振比的平面及三維共振軌道進行設(shè)計,并分析共振軌道隨軌道周期及能量的演化規(guī)律。本文的研究可以為未來雙小行星系統(tǒng)探測軌道設(shè)計提供參考。
本文采用雙橢球體模型表征同步雙小行星系統(tǒng),建立質(zhì)心旋轉(zhuǎn)坐標系O-xyz如圖1所示。其中A1為雙體小行星系統(tǒng)中質(zhì)量較大的小行星(主星),A2為質(zhì)量較小的小行星(從星),B為航天器。x軸沿A1,A2之間的連線由A1指向A2,z軸平行于軌道角動量矢量方向,y軸與x軸、z軸構(gòu)成右手坐標系。
圖1 雙小行星系統(tǒng)橢球體-橢球體模型Fig.1 The double ellipsoid model of binary asteroid
假設(shè)A1與A2同步自旋,且慣性主軸均位于x軸上。它們的半主軸分別為a1,b1,c1及a2,b2,c2,且滿足0≤c1≤b1≤a1,0≤c2≤b2≤a2。記A1和A2的質(zhì)量分別為M1,M2,m(m?M1,M2),則系統(tǒng)的質(zhì)量比為:
(1)
(2)
Rb=[x,y,z]T為航天器相對于系統(tǒng)質(zhì)心的位置矢量。
(3)
則探測器在質(zhì)心旋轉(zhuǎn)系下的運動方程可表示為:
(4)
Us=(1-μ)U1(ρ-R1)+μU2(ρ-R2)
(5)
其中,U1,U2分別代表雙小行星系統(tǒng)中兩顆均質(zhì)三軸橢球體的外部勢能,
(6)
(7)
式中:λ1,λ2分別為φ(ρ-R1,λ1)與φ(ρ-R2,λ2)的根,v為積分變量。
(8)
(9)
(10)
(11)
式(4)表示成分量形式為:
(12)
式中:Ea,Eb,Ec為代表雙橢球質(zhì)量分布的橢圓積分,具體求解過程可參考文獻[15]。
首先忽略A2的質(zhì)量,基于二體模型求解共振軌道,作為雙小行星系統(tǒng)內(nèi)共振軌道的初值,假設(shè)探測器B繞A1飛行p圈所需的時間與A2繞A1飛行q圈所需的時間相同,p和q均為正整數(shù),B與A2的軌道周期比為:
(13)
其中,np,nq分別為B與A2的平均角速度,
(14)
Tp,Tq分別為B與A2的周期。μ1=GM1為A1的引力常數(shù)。根據(jù)此定義,在慣性系中,探測器初始狀態(tài)可由式(15)~(18)確定:
(15)
p1=a(1-e2)
(16)
(17)
(18)
其中,p1為半正焦弦,r0為初始距離大小,v0為初始速度。θ為真近角,在近心點處,θ=0°,在遠心點處,θ=180°。探測器B的周期Tp可以由共振比p:q,A2周期Tq和式(13)計算得到:
(19)
共振軌道的確定還需要知道偏心率e的大小,選擇近心點或遠心點作為探測器的初始狀態(tài),可更方便地得到平面共振軌道。因此本文將通過選取近心點位置rp(或遠心點位置ra)和已求得的半長軸a來求解偏心率e。
若初始狀態(tài)完全確定,即可解析求得共振軌道在慣性坐標系中的解。假設(shè)雙小行星中主從星的半徑分別為r1=12.38 km和r2=3.60 km,相對距離R=33 km,密度ρ=2.91 g/cm3,同步旋轉(zhuǎn)周期為2.917 h。初始時刻為A1,A2和B三者共線的時刻,可求得僅考慮主星質(zhì)量下的3∶4共振軌道,選擇雙星相對距離作為單位長度,在旋轉(zhuǎn)坐標系下共振軌道如圖2所示。
圖2 二體系統(tǒng)中3∶4共振軌道Fig.2 3∶4 resonant orbit in binary asteroid system based on two-body model
將第1節(jié)二體模型中求得的共振軌道作為初值,利用改進并行打靶法修正可得雙小行星系統(tǒng)下的共振軌道。修正分兩步進行,第一次修正得到忽略小行星形狀,僅考慮從星質(zhì)量的共振軌道,并將它的初始狀態(tài)作為第二次修正的初值,從而得到雙橢球模型下具有周期性的共振軌道。然后,將該共振軌道作為初值,利用改進的連續(xù)法延拓得到一組共振比相同的共振軌道族。最后,利用動力學(xué)分岔理論,沿共振軌道的特征方向延拓生成新的空間共振軌道族。
由于雙小行星系統(tǒng)內(nèi)動力學(xué)環(huán)境的復(fù)雜性,共振軌道對初值敏感,需要對二體模型下的軌道初值進行修正。本文采用并行打靶法進行軌道設(shè)計,該方法相對傳統(tǒng)的打靶法而言,有更好的穩(wěn)定性和更大的收斂域;同時相對于常用的差分方法,要求解的非線性方程組的規(guī)模要小得多、計算速度快。
圖3 改進的并行打靶法示意圖Fig.3 Improved parallel shooting method
(20)
引入松弛變量β,則約束方程可表示為
(21)
終端點的約束向量為
(22)
為了得到周期性的軌道,打靶過程中第一個節(jié)點一定在超平面上,因此有y1=yh=0。
將積分時間記作T,則改進的多次打靶法中的自由變量X為
(23)
X由6n+2個元素構(gòu)成。約束函數(shù)F(X)為
F(X)由6n+1個分量組成。則Jacobi矩陣DF(X)為(6n+1)×(6n+2)維矩陣,具體表達式為:
DF(X)=
(24)
其中Φ(t2,t1)表示從時刻t1至t2的狀態(tài)轉(zhuǎn)移矩陣,
A=diag(-1, -1, -1, -1, 0, 1)
(25)
(26)
(27)
利用改進的并行打靶法,將初始軌道分成若干段,首先考慮雙小行星系統(tǒng)中的子星的質(zhì)量,將系統(tǒng)的質(zhì)量比由0變?yōu)檎鎸嵵颠M行第一次修正,圖4給出不同質(zhì)量比下雙小行星系統(tǒng)中的3∶4共振軌道的演化情況。
圖4 μ變化下共振軌道演化Fig.4 Evolution of resonance orbits under varying μ
系統(tǒng)質(zhì)量比μ從0逐漸增大為真實值0.0245,圖4(a)給出了不同系統(tǒng)質(zhì)量比下的共振軌道演化情況,發(fā)現(xiàn)隨著質(zhì)量比μ的增大,共振軌道所對應(yīng)的Jacobi常數(shù)增大,同時共振軌道與主星間距離逐漸增大,與子星的距離減小。μ=0.0245時球體-球體模型中的最終修正軌道如圖4(b)所示。
進一步考慮小行星形狀變化,可得雙橢球模型下的共振軌道。假設(shè)雙小行星的形狀參數(shù)為
a1=7.25 km,b1=5.90 km,c1=5.55 km
a2=1.90 km,b2=1.75 km,c2=1.75 km
圖5(a)給出形狀變化下雙小行星系統(tǒng)中的3∶4共振軌道的演化情況,隨著主從星形狀連續(xù)變化,3∶4共振軌道所對應(yīng)的Jacobi常數(shù)值不斷增大。軌道同時遠離主星和子星,但偏移距離小于質(zhì)量比變化造成的改變。利用改進并行打靶法和兩步修正策略,可以實現(xiàn)雙橢球模型下精確共振軌道的設(shè)計,最終的共振軌道如圖5(b)所示。
圖5 形狀變化下共振軌道演化Fig.5 Evolution of resonance orbits under varying shape
采用連續(xù)法可以從已知共振軌道中延拓出相同共振比的軌道,降低了軌道設(shè)計難度,實現(xiàn)共振軌道族的快速生成。本文采用偽弧長連續(xù)法對軌道族進行求解[16-17]。偽弧長連續(xù)法在與軌道族相切的方向上取連續(xù)變量,保證所求解處在同一平面內(nèi),避免解的空間不確定性。
(28)
其中,Δk為選取的與軌道族相切方向的步長,偽弧長延拓法中的約束矩陣表達式如下:
(29)
(30)
(31)
偽弧長連續(xù)法具有魯棒性強的優(yōu)點,可保證每一次迭代所求解的唯一性。根據(jù)第2.1節(jié)確定的3∶4共振軌道,通過偽弧長連續(xù)法可求得共振軌道族如圖6所示。
圖6 雙小行星3∶4共振軌道族Fig.6 3∶4 resonance orbit family
軌道的穩(wěn)定性是判斷一條軌道是否適用于探測任務(wù)的重要依據(jù)之一,穩(wěn)定共振軌道的軌道維持所需的燃料較省。此外,共振軌道族的穩(wěn)定性改變對應(yīng)于軌道特征結(jié)構(gòu)的變化,可用于尋找新的共振軌道族,從而為任務(wù)軌道的設(shè)計提供更多思路。
利用單值矩陣可對軌道穩(wěn)定性進行分析,定義單值矩陣M為一個軌道周期T處的狀態(tài)轉(zhuǎn)移矩陣值,即
(32)
F(λ)=λ4-αλ3+βλ2-αλ+1
(33)
則剩余4個特征值可令F(λ)=0求得,當|λi|≤1(i=3,…,6)時,初始狀態(tài)加上其對應(yīng)的特征向量方向上的微小變化,周期軌道穩(wěn)定;當|λi|>1(i=3,…,6)時,周期軌道不穩(wěn)定。引入穩(wěn)定性指數(shù)ν,定義穩(wěn)定性指數(shù)ν為互為倒數(shù)的一對特征根的平均值,即
(34)
穩(wěn)定性指數(shù)ν可用來分析給定解的穩(wěn)定性。對于給定的單值矩陣來說,除了均為1的一對特征根,由另外兩對特征根可得到兩個不同的穩(wěn)定性指數(shù)νh,νv。其中,“水平”穩(wěn)定性指數(shù)ν=νh反映了軌道平面中的穩(wěn)定性,“垂直”穩(wěn)定性指數(shù)ν=νv反映了軌道平面外的穩(wěn)定性。同理,如果穩(wěn)定性指數(shù)ν大于1,則軌道穩(wěn)定,反之則軌道不穩(wěn)定。通過引入穩(wěn)定性指數(shù),可以較為全面地考慮單值矩陣特征根對軌道穩(wěn)定性帶來的影響,特別地,考慮引入最大穩(wěn)定性指數(shù)νmax=max{|νh|,|νv|},則若νmax≤1,軌道穩(wěn)定,反之則軌道不穩(wěn)定。
對于一組軌道族的特征根來說,若其結(jié)構(gòu)發(fā)生顯著改變,軌道的穩(wěn)定性突然發(fā)生改變,表明軌道可能出現(xiàn)分岔[18],利用分岔理論可以從分岔軌道延拓得到新的共振軌道族。在共振軌道中常見的分岔類型包括正切分岔及倍增分岔。
圖7為正切分岔與倍增分岔的示意圖,當一對原本位于單位圓上(即|λ|=1)的共軛復(fù)數(shù)特征根在復(fù)平面上的(-1,0)處(即ν=-1處)交匯,之后沿著實數(shù)軸Re(λ)變化時,出現(xiàn)正切分岔。一般在正切分岔中,軌道穩(wěn)定階數(shù)將發(fā)生改變。
圖7 不同類型的分岔形式Fig.7 Different types of bifurcation
當一對原本位于單位圓上的共軛復(fù)數(shù)特征根在復(fù)平面上的點(1,0)處(即ν=1處)交匯,之后沿著實數(shù)軸Re(λ)變化時,出現(xiàn)倍增分岔。在倍增分岔前后,軌道穩(wěn)定階數(shù)不變,新的軌道族的周期為分岔點處初始軌道周期的兩倍。
利用分岔理論,可在二維平面共振軌道的基礎(chǔ)上對雙小行星系統(tǒng)中的三維空間中的共振軌道進行設(shè)計。選取出現(xiàn)分岔的共振軌道,計算其狀態(tài)轉(zhuǎn)移矩陣將這條分岔軌道初始狀態(tài)的z分量加一個小的攝動量。固定攝動分量的值,利用改進的多次打靶法可找到一條三維共振軌道,利用連續(xù)法即可得到三維的共振軌道。
根據(jù)第2節(jié)提出的共振軌道計算與延拓方法,本文選取的雙體小行星系統(tǒng)1999 KW4[19],求解其附近的平面和空間共振軌道族,并討論共振軌道特性隨軌道周期及能量的演化規(guī)律。雙體小行星系統(tǒng)1999 KW4的三軸橢球體形狀模型及相關(guān)參數(shù)如表1所示,系統(tǒng)質(zhì)量比為0.0541。
圖8給出橢球體-橢球體模型中的平面共振軌道,選取的共振比均小于等于1,因為大于1的共振比的軌道大部分與主星相交,即發(fā)生碰撞,且周期小、軌道速度快,對工程探測而言意義不大。
表1 1999 KW4三軸橢球體模型Table 1 Tri-axial ellipsoid model of 1999 KW4
圖8 雙體小行星不同共振比的平面軌道族Fig.8 Resonant orbit families with different resonant ratios
從圖8所示的共振軌道圖可以看出,環(huán)的個數(shù)總是和共振比p∶q中p的值一致,即共振比為1∶1,1∶2,1∶3,1∶4的共振軌道中均只有一個環(huán),而共振比為2∶3的共振軌道有2個環(huán),共振比為3∶4的共振軌道有3個環(huán)。
由于同一族軌道的各參數(shù)連續(xù)變化,所以選取各圖中軌道延拓的內(nèi)外邊界軌道(非軌道族的邊界)進行分析,表2列出了它們的非零初始狀態(tài)、周期和雅克比常數(shù)C的值。
表2 平面共振軌道的初始狀態(tài)、軌道周期和雅克比常數(shù)值Table 2 Initial state, period and Jacobi constant of the planar resonant orbits
由圖9可知,三維共振軌道族與對應(yīng)的平面共振軌道族的形狀保持一定的相似性;共振比為1∶1的分岔軌道的初始位置與主星較近,共振比為1∶2,1∶3,1∶4的分岔軌道的初始位置則較靠近從星,共振比為2∶3的分岔軌道有兩族,其中一族的初始位置靠近主星,另一族的初始位置靠近從星。
表3中列出了三維空間共振軌道的非零初始狀態(tài)、周期和Jacobi常數(shù)C的值。從表3可以看出,隨著三維共振軌道初始狀態(tài)z0分量絕對值增大,軌道周期增大、Jacobi常數(shù)減小。從軌道周期的數(shù)據(jù)可以看出,探測器周期與雙小行星系統(tǒng)的公轉(zhuǎn)周期之比均不為整數(shù)比,隨著軌道族的延拓(往z0分量的絕對值增大的方向)更接近于整數(shù)比。
本文主要研究了雙體小行星中共振軌道設(shè)計方法。利用改進并行打靶法和分岔理論,提出雙小行星系統(tǒng)平面共振軌道修正方法和三維共振軌道生成方法,并利用連續(xù)法實現(xiàn)共振軌道族的快速延拓。以雙小行星系統(tǒng)1999 KW4為例,設(shè)計了共振比為1∶1,1∶2,1∶3,1∶4,2∶3的平面和空間共振軌道族,并分析了共振軌道的特性及軌道周期和軌道能量的變化規(guī)律。研究發(fā)現(xiàn)雙小行星系統(tǒng)中共振軌道“軌道環(huán)”的個數(shù)和共振比p∶q中p的值總是保持一致。由于雙小行星系統(tǒng)復(fù)雜的動力學(xué)環(huán)境與特性,探測器的軌道周期與較小天體的公轉(zhuǎn)周期比不為整數(shù),但隨著軌道能量的增大,比值將趨近于整數(shù)。本文給出的共振軌道設(shè)計方法和軌道運動規(guī)律研究,可為未來雙小行星系統(tǒng)探測的任務(wù)軌道設(shè)計提供一定的借鑒。
表3 共振軌道的初始狀態(tài)、軌道周期和雅克比常數(shù)值Table 3 Initial state, period and Jacobi constant of three dimensional resonant orbits