劉磊 劉勇 曹建峰 唐歌實(shí) 胡松杰

（1航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京100094）（2北京航天飛行控制中心，北京100094）

1 引言

限制性三體問(wèn)題的5個(gè)平動(dòng)點(diǎn)及其附近的周期軌道由于獨(dú)特的空間位置和動(dòng)力學(xué)特性，非常適合于空間觀測(cè)和中繼通信等應(yīng)用，或者作為深空探測(cè)的轉(zhuǎn)移中樞[1]。目前為止，國(guó)際上已經(jīng)成功發(fā)射了10個(gè)平動(dòng)點(diǎn)探測(cè)器；2011年我國(guó) “嫦娥二號(hào)”衛(wèi)星在探月任務(wù)結(jié)束后，也進(jìn)入了日-地月系L2點(diǎn)擬周期軌道開(kāi)展探測(cè)[2]；此外國(guó)內(nèi)外還有多項(xiàng)平動(dòng)點(diǎn)探測(cè)任務(wù)處于計(jì)劃之中。在這些平動(dòng)點(diǎn)任務(wù)的設(shè)計(jì)階段，為了滿足航天測(cè)控能力約束，尤其是探測(cè)任務(wù)目標(biāo)的需求，必須首先了解目標(biāo)平動(dòng)點(diǎn)周期軌道的周期和運(yùn)動(dòng)范圍等特性。

平動(dòng)點(diǎn)周期軌道的理論基礎(chǔ)為天體力學(xué)中的三體問(wèn)題，尤其是圓型限制性三體問(wèn)題，目前國(guó)內(nèi)外對(duì)其動(dòng)力學(xué)特性分析和相關(guān)解的研究較為充分[3-7]，在Halo軌道族計(jì)算尤其是族參數(shù)選擇、族周期和幅值的變化特性等方面的研究則相對(duì)有限。為了全面了解平動(dòng)點(diǎn)大范圍Halo軌道的周期和運(yùn)動(dòng)范圍等特性，可以研究平動(dòng)點(diǎn)附近的Halo軌道族。不過(guò)，由于Halo軌道族計(jì)算受族參數(shù)影響較大，在某些特定族參數(shù)下大幅值Halo軌道迭代不收斂，或者迭代結(jié)果與初值相差較大，不能得到完整連續(xù)的Halo軌道族。

因此，本文基于微分修正和延拓法研究了Halo軌道族計(jì)算及其運(yùn)動(dòng)特性，克服了傳統(tǒng)Halo軌道設(shè)計(jì)方法在大幅值Halo軌道和完整Halo軌道族數(shù)值計(jì)算上的缺陷，得到了日-地月系，尤其是地月系共線平動(dòng)點(diǎn)連續(xù)完整的大范圍南北Halo軌道族，可用于我國(guó)未來(lái)深空探測(cè)的平動(dòng)點(diǎn)任務(wù)設(shè)計(jì)。

2 數(shù)學(xué)模型

2.1 Halo軌道設(shè)計(jì)

限制性三體問(wèn)題中，第三體的質(zhì)量遠(yuǎn)小于兩個(gè)大天體的質(zhì)量。首先將質(zhì)量、長(zhǎng)度和時(shí)間單位歸一化[6]，記μ為質(zhì)量稍小的大天體相對(duì)于兩個(gè)大天體質(zhì)量之和的比值；其次以兩大天體的質(zhì)心為原點(diǎn)O，兩大天體的連線為x軸，相對(duì)運(yùn)動(dòng)平面為xy平面，建立旋轉(zhuǎn)坐標(biāo)系或稱會(huì)合坐標(biāo)系。在圓型限制性三體問(wèn)題下，建立小天體在坐標(biāo)系Oxyz中的動(dòng)力學(xué)方程：

式中Ω為等效勢(shì)[6]。

在平動(dòng)點(diǎn)附近的運(yùn)動(dòng)相對(duì)于平動(dòng)點(diǎn)的偏差Δx、Δy和Δz為小量的情況下，將式（1）在平動(dòng)點(diǎn)附近作一階展開(kāi)，得到平動(dòng)點(diǎn)附近周期運(yùn)動(dòng)的解析形式為

式中Ax和Az分別為x和z方向的運(yùn)動(dòng)幅值；η和ζ分別為x和z方向的運(yùn)動(dòng)頻率，當(dāng)運(yùn)動(dòng)幅值充分大時(shí)，非線性項(xiàng)的影響可能使得二者相等，從而形成Halo軌道；φ和ψ分別為x和z方向的初始相位。

由式（2）可見(jiàn)，Halo軌道關(guān)于xz平面對(duì)稱，且在xz平面附近運(yùn)動(dòng)方向垂直于xz平面，這一性質(zhì)對(duì)于Halo軌道的數(shù)值設(shè)計(jì)至關(guān)重要。Halo軌道數(shù)值設(shè)計(jì)的迭代初值可采用式（2）。如果對(duì)于初值要求較高，則可以基于式（2），由Lindstedt-Poincaré法得到共線平動(dòng)點(diǎn)Halo軌道的三階近似解[6]。

微分修正中，利用Halo軌道的對(duì)稱性，可選擇初始狀態(tài)x0位于xz平面，在動(dòng)力學(xué)模型式（1）下由x0積分t時(shí)刻至末狀態(tài)xd，通過(guò)修正x0和t使得xd位于xz平面且滿足Halo軌道特性，即x和z方向速度為0，設(shè)xd的偏差量為δxd，修正量分別為δx0和δt，則迭代過(guò)程為

式中Φij（i＝2，4，6；j＝1，3，5）表示狀態(tài)轉(zhuǎn)移陣Φ（t1；t0）的第i行第j列個(gè)分量[3]。Φ（t1；t0）建立了由初始狀態(tài)到末狀態(tài)的單步線性估計(jì)，也體現(xiàn)了末狀態(tài)對(duì)初始狀態(tài)的敏感性，如果考慮其余天體引力和光壓等攝動(dòng)，則修改式（1）右端即可求解出對(duì)應(yīng)的狀態(tài)轉(zhuǎn)移矩陣，即得完整力模型下的Halo軌道（族）。

2.2 Halo軌道族參數(shù)與延拓計(jì)算

式（3）為欠定方程組，一般固定一個(gè)待求參數(shù)后有確定解，這個(gè)固定參數(shù)即定義為Halo軌道族的族參數(shù)，由此可選擇x0、z0、0或周期T作為族參數(shù)，不過(guò)有些情況下會(huì)出現(xiàn)迭代不收斂情況。同時(shí)需指出，除族參數(shù)外的其余參數(shù)應(yīng)同步求解，文獻(xiàn)[3]中雖然也是固定z0，但是迭代時(shí)先積分至xz平面，然后改正x0和z0使得末狀態(tài)和為0，該方法對(duì)于大幅值Halo軌道會(huì)出現(xiàn)迭代發(fā)散現(xiàn)象，其原因在于預(yù)先施加約束而構(gòu)造的迭代序列偏離了真實(shí)的軌道解。

式（3）還可以采用最小二乘法求解，其優(yōu)點(diǎn)在于收斂性較好，大幅值Halo軌道也可以收斂。不過(guò)，由于所有參數(shù)均參與迭代改正，沒(méi)有反映族特征的明顯族參數(shù)，此時(shí)雖然可以考慮選擇相關(guān)的Halo軌道幅值A(chǔ)x或Az作為族參數(shù)，然而數(shù)值計(jì)算表明，軌道幅值較大時(shí)的最終結(jié)果與迭代初值相差較大，往往收斂于平面Lyapunov軌道，且得到的Halo軌道族連續(xù)性較差。

選定族參數(shù)后，基于式（3）計(jì)算Halo軌道族。鑒于微分修正法極大地依賴于x0與真實(shí)解的逼近程度，大幅值Halo軌道計(jì)算時(shí)可能出現(xiàn)迭代發(fā)散，考慮采用延拓法計(jì)算Halo軌道族初值。延拓法是一種求解非線性方程組的大范圍收斂算法[8]，對(duì)x0與真實(shí)解的逼近程度要求較低。對(duì)于如下非線性方程組

延拓法的思想就是構(gòu)造一族映象H∶D×[0，1]?Rn＋1→Rn，且有

其中F0X（）＝0的解X0已知，H（X，1）＝0的解即式（4）的解，于是式（4）的求解轉(zhuǎn)化為求如下同倫方程組的解Xν（）

式（3）對(duì)應(yīng)于延拓法中的待求問(wèn)題F（X），除族參數(shù)以外的變量即待求參數(shù)X，族參數(shù)即同倫方程組變量ν。利用延拓法計(jì)算Halo軌道族，首先在平動(dòng)點(diǎn)附近利用Halo軌道高階解和微分修正法，計(jì)算初始標(biāo)稱軌道，即式（5）中F0X（）＝0的解X0；其次，選擇同倫方程組變量ν，選取步長(zhǎng)Δν，利用式（3）計(jì)算族參數(shù)增加步長(zhǎng)Δν后X0的其余狀態(tài)參數(shù)的改變量，從而得到新的X0；最后，利用微分修正法迭代求解族參數(shù)增加步長(zhǎng)Δν后的Halo軌道。

為了分析Halo軌道族的運(yùn)動(dòng)區(qū)域，定義微分修正后的Halo軌道在x和z方向的幅值為

3 數(shù)值仿真

由于地月系的μ值遠(yuǎn)大于日-地月系，因而二者Halo軌道族特性差別較大，同時(shí)L1和L2點(diǎn)附近的動(dòng)力學(xué)特性差別使得各自附近的運(yùn)動(dòng)也不相同。因此，本節(jié)數(shù)值計(jì)算地月系和日-地月系L1和L2點(diǎn)的南北Halo軌道族，給出軌道族的周期和運(yùn)動(dòng)幅值變化，通過(guò)計(jì)算給出最大的固定延拓步長(zhǎng)，同時(shí)研究族參數(shù)選擇對(duì)Halo軌道族計(jì)算的影響。

3.1 日-地月系Halo軌道族

選擇x0作為族參數(shù)計(jì)算L1點(diǎn)北族Halo軌道，如圖1所示。圖1中○表示各Halo軌道的初始位置狀態(tài)，●和■分別表示L1點(diǎn)和地月系質(zhì)心，曲線為軌道族中的Halo軌道。

由圖1可見(jiàn)，Halo軌道族在x方向由L1點(diǎn)拓展至地月系質(zhì)心，但是y和z方向上運(yùn)動(dòng)范圍有限。L1點(diǎn)附近的Halo軌道相對(duì)于yz面傾角最小，在地月系質(zhì)心附近傾角最大。另外，仿真中發(fā)現(xiàn)延拓步長(zhǎng)最大可達(dá)10-5，如果進(jìn)一步增大步長(zhǎng)，則對(duì)初始軌道要求苛刻，迭代一般收斂于平面Lyapunov軌道。

L1點(diǎn)北族Halo軌道族的軌道周期和運(yùn)動(dòng)幅值隨族參數(shù)x0的變化如圖2中的T1和Sx1及Sz1所示。為直觀起見(jiàn)，圖2中將x0轉(zhuǎn)換為相對(duì)于地月系質(zhì)心的距離dx。以x0作為族參數(shù)計(jì)算L1點(diǎn)南族Halo軌道族，其周期T和運(yùn)動(dòng)幅值的變化與圖2中的T1和Sx1及Sz1完全相同，充分反映了南北Halo軌道族的對(duì)稱性。此外，以x0作為族參數(shù)計(jì)算L2點(diǎn)南北族Halo軌道族，其周期T約為110～180天，運(yùn)動(dòng)幅值Sx和Sz與L1點(diǎn)Halo軌道族基本相同。

選擇0作為族參數(shù)無(wú)法計(jì)算得到L1點(diǎn)南北族整個(gè)Halo軌道族。不過(guò)，可以獲得L2點(diǎn)南北族整個(gè)Halo軌道族，尤其是可以將軌道族延拓到地月系質(zhì)心附近很小范圍內(nèi)，軌道周期和幅值的變化如圖2中的T2和運(yùn)動(dòng)幅值Sx2及Sz2所示。

圖1 日-地月系L1點(diǎn)北族Halo軌道族Fig.1 North halo family of the libration point L1in the Sun-Earth／Moon system

圖2 日-地月系L1點(diǎn)和L2點(diǎn)Halo軌道族的周期與幅值Fig.2 Periods and scopes of the halo family of the libration point L1and L2in the Sun-Earth／Moon system

由圖2可見(jiàn)，周期T2的變化范圍大于T1，二者最大值較為接近，但是T2最小可達(dá)約48天。Halo軌道族的運(yùn)動(dòng)幅值Sz大于Sx，Sz最大可達(dá)1.2×106km，Sx最大僅為5×105km。另外，雖然L1和L2點(diǎn)Halo軌道族在x和z方向的運(yùn)動(dòng)幅值的變化趨勢(shì)有所差異，但是它們的變化范圍大致相同，且均在平動(dòng)點(diǎn)與地月系質(zhì)心之間的某個(gè)位置出現(xiàn)最大值。

如果選擇z0作為族參數(shù)，則無(wú)法得到L1點(diǎn)和L2點(diǎn)的整個(gè)南北Halo軌道族。經(jīng)過(guò)分析可知，Halo軌道族在xz平面上沿z軸并非單調(diào)變化，固定步長(zhǎng)延拓失效，需要采用較為復(fù)雜的延拓步長(zhǎng)調(diào)整策略。如果選擇T作為族參數(shù)，可以得到L1點(diǎn)附近絕大部分南北Halo軌道族。至于L2點(diǎn)附近，選擇T作為族參數(shù)的結(jié)果優(yōu)于選擇x0，其原因在于固定延拓步長(zhǎng)下x0失效時(shí)仍然可以用T延拓。

本文計(jì)算了日-地月系L1和L2點(diǎn)的南北Halo軌道族，每個(gè)軌道族均用4個(gè)族參數(shù)進(jìn)行計(jì)算，共計(jì)16種情況，限于篇幅，且南北軌道族關(guān)于xy平面對(duì)稱，僅給出L1和L2點(diǎn)北族Halo軌道族仿真結(jié)果，如表1所示，其中Spx和Spz分別為Halo軌道族在x和z方向上距地月系質(zhì)心的變化范圍，P為軌道周期范圍，τ為最大延拓步長(zhǎng)。除τ外，其余各參數(shù)均給出變化范圍。

表1 日-地月系Halo軌道族特性Tab.1 Characteristics of the halo family in the Sun-Earth／Moon system

由表1可見(jiàn)，Sx的最小范圍在以0為族參數(shù)時(shí)最小可達(dá)到約700km，選擇其余族參數(shù)時(shí)最大不超過(guò)2.2×105km，文獻(xiàn)[6]中由三階解析解給出的最小值約為2×105km，由此可見(jiàn)數(shù)值計(jì)算得到的Halo軌道族可達(dá)范圍遠(yuǎn)大于三階解析解結(jié)果，同時(shí)，日-地月系L1點(diǎn)Halo軌道族應(yīng)選擇T或x作為族參數(shù)，L2點(diǎn)軌道族應(yīng)選擇或T作為族參數(shù)。表1中的各項(xiàng)參數(shù)可為平動(dòng)點(diǎn)任務(wù)設(shè)計(jì)提供非常有益的參考。

3.2 地月系Halo軌道族

地月L1點(diǎn)南北軌道族，以z0作為族參數(shù)的范圍最大，如圖3所示，注意此處AU為地月距，■表示月球。

由圖3可見(jiàn)，軌道族從L1點(diǎn)向月球拓展，到達(dá)月球后仍可進(jìn)一步拓展，z0最大可以達(dá)到1個(gè)地月距左右，將之與圖2對(duì)比，充分說(shuō)明了不同三體系統(tǒng)的周期軌道運(yùn)動(dòng)范圍并不完全相同。軌道周期T由13.5天降低到最低約8天后，又快速增加到約12天，Sx和Sz的變化也較為復(fù)雜，且南族Halo軌道族的變化與北族類似，限于篇幅，此處略去變化圖。

以其余族參數(shù)得到的L1軌道族中，僅以x0為族參數(shù)的結(jié)果與以z0為族參數(shù)的結(jié)果接近，以0和T為族參數(shù)的結(jié)果較差。對(duì)于L2點(diǎn)南北軌道族，以0作為族參數(shù)的結(jié)果范圍最大，其周期和幅值變化趨勢(shì)與以0為族參數(shù)的日-地月系L2點(diǎn)結(jié)果相同。以x0和T為族參數(shù)的結(jié)果略差于以0為族參數(shù)的結(jié)果，以z0為族參數(shù)的結(jié)果最差。因此，地月系L1點(diǎn)Halo軌道族應(yīng)選擇z或x作為族參數(shù)，L2點(diǎn)軌道族應(yīng)選擇y方向速度或T作為族參數(shù)。

地月系L1和L2點(diǎn)的北Halo軌道族特性如表2所示，其中僅給出較好的族參數(shù)結(jié)果。

圖3 地月系L1點(diǎn)北族Halo軌道族Fig.3 North halo family of the libration point L1in the Earth-Moon system

表2 地月系Halo軌道族特性Tab.2 Characteristics of the halo family in the Earth-Moon system

由表2可見(jiàn)，地月系L1點(diǎn)Halo軌道族在x和z方向的可達(dá)范圍和運(yùn)動(dòng)幅值遠(yuǎn)大于L2點(diǎn)Halo軌道族，不過(guò)軌道周期的變化范圍小于L2點(diǎn)Halo軌道族。另外，L2點(diǎn)Halo軌道族的延拓步長(zhǎng)最大可以達(dá)到0.1，L1點(diǎn)也可以達(dá)到10-3，遠(yuǎn)遠(yuǎn)大于表1中日-地月系情況。

4 結(jié)束語(yǔ)

在圓型限制性三體問(wèn)題下，基于延拓法研究了固定步長(zhǎng)下族參數(shù)選擇對(duì)Halo軌道族計(jì)算的影響，給出了各Halo軌道族的軌道周期和空間位置變化特性。研究表明，日-地月系L1點(diǎn)Halo軌道族應(yīng)選擇軌道周期T或會(huì)合坐標(biāo)系x坐標(biāo)作為族參數(shù)，L2點(diǎn)軌道族應(yīng)選擇y方向速度或T作為族參數(shù)；地月系L1點(diǎn)Halo軌道族應(yīng)選擇z或x作為族參數(shù)，L2點(diǎn)軌道族應(yīng)選擇y方向速度或T作為族參數(shù)，且數(shù)值計(jì)算得到的Halo軌道族運(yùn)動(dòng)范圍大于解析解。

文中僅用固定延拓步長(zhǎng)即可得到大幅值周期軌道和連續(xù)完整的周期軌道族，可用于任意三體系統(tǒng)平動(dòng)點(diǎn)周期軌道族計(jì)算，不過(guò)研究中發(fā)現(xiàn)固定延拓步長(zhǎng)情況下族參數(shù)的選擇影響很大，如果不是單調(diào)變化，則無(wú)法得到整個(gè)Halo軌道族。因此，后續(xù)研究計(jì)劃將考慮族參數(shù)的變化特性，制定變步長(zhǎng)計(jì)算策略。

[1]LO M W.The interplanetary superhighway and the origins program，IEEE Space 2002Conference[C].IEEE，Big Sky，MT，US，March 9-16，2002.

[2]WU WEIREN，LIU YONG，LIU LEI，et al.Pre-LOI trajectory maneuvers of CHANG′E-2libration point mission [J].SCIENCE CHINA，2012，55（6）：1249-1258.

[3]CHOW NAOMI，GRALLA ERICA.Low earth orbit constellation design using the Earth-Moon L1point[D].New Jersey：Princeton University，2004.

[4]HOWELL K C.Families of orbits in the vicinity of the collinear libration points[J].Journal of the Astronautic Sciences，2001，49（1）：107-125.

[5]CORS J M，PINYOL C，SOLER J.Periodic solutions in the spatial elliptic restricted three-body problem [J].Physica D，2001，154：195-206.

[6]RICHARDSON D L.Analytic construction of periodic orbits about the collinear points [J].Celestial Mechanics and Dynamical Astronomy，1980，22（3）：241-253.

[7]侯錫云.平動(dòng)點(diǎn)的動(dòng)力學(xué)特征及其應(yīng)用 [D].南京：南京大學(xué)，2008.HOU XIYUN.Dynamics and their applications of libration points[D].Nanjing：Nanjing University，2008.

[8]李慶揚(yáng)，莫孜中，祁力群.非線性方程組的數(shù)值解法 [M].北京：科學(xué)出版社，1997.LI QINGYANG，MO ZIZHONG，QI LIQUN.Numerical methods for non-linear equation group [M].Beijing：Science Press，1997.