一種適用于大幅值Quasi周期軌道的數(shù)值構(gòu)造方法

李瑞龍，朱戰(zhàn)霞

(西北工業(yè)大學(xué)航天學(xué)院，西安 710072)

0 引 言

近年來，各航天大國都在加速布局地月空間開發(fā)與深空探測[1-2]，平動點(diǎn)任務(wù)是其中的重要組成部分，各航天大國均制定了未來數(shù)年的平動點(diǎn)任務(wù)規(guī)劃。平動點(diǎn)及周期軌道是圓型限制性三體問題(CRTBP)的特殊動力學(xué)結(jié)構(gòu)，可以滿足特定的空間科學(xué)任務(wù)的需求。以美國的“阿爾忒彌斯”計(jì)劃為例，提出把Gateway空間站部署在地月系統(tǒng)的近直線暈軌道(NRHO)上，為載人登月實(shí)現(xiàn)中轉(zhuǎn)[3-4]。

對CRTBP來說，quasi周期軌道相比周期軌道是更豐富且普遍的存在[5]。在一條周期軌道的附近存在大量的quasi周期軌道，它們不僅具有周期軌道的眾多優(yōu)點(diǎn)，而且增廣的自由度擴(kuò)展了軌道的設(shè)計(jì)空間，quasi周期運(yùn)動也是對高精度動力學(xué)模型更好的近似[6]。此外，由于分布在高維的不變環(huán)面上，使得其在軌道轉(zhuǎn)移及編隊(duì)飛行等方面更具有先天的動力學(xué)優(yōu)勢[7]?；谝陨咸匦?，以往的平動點(diǎn)任務(wù)也多次采用了quasi周期軌道作為任務(wù)軌道，比如NASA的SOHO[8]和“阿爾忒彌斯”[9]等。隨著未來平動點(diǎn)飛行任務(wù)的科學(xué)需求與約束日益復(fù)雜化，quasi周期軌道能夠提供更加靈活的標(biāo)稱軌道參考，具有廣闊的應(yīng)用價(jià)值。

然而，quasi周期軌道更加復(fù)雜的動力學(xué)特性給計(jì)算帶來了更大的難度，尤其是對大幅值quasi周期軌道而言。已有不少學(xué)者研究了quasi周期軌道的計(jì)算問題，提出了半解析法和全數(shù)值法兩大類方法。Farquhar等[10]最先采用Lindstadt-Poincaré攝動法給出了quasi-Halo軌道的三階近似解析解，這種方法的精度較為有限，后續(xù)需要配合數(shù)值方法進(jìn)一步修正。Gómez等[11]采用攝動法計(jì)算了Halo軌道附近的quasi-Halo軌道，他們構(gòu)造了軌道的25階近似解析解，并提出多步打靶的思想進(jìn)一步對解析解進(jìn)行修正，該方法可以得到較為準(zhǔn)確的quasi-Halo軌道，但是解析解的構(gòu)造較為繁瑣，且通用性較低。Jorba等[12]提出用中心流形約化的方法，以Fourier級數(shù)的截?cái)嘈问接?jì)算quasi周期軌道。然而，無論是攝動法還是Fourier級數(shù)，都存在一定程度的截?cái)嗯c近似，后續(xù)以此為初值進(jìn)一步用數(shù)值方法修正時(shí)，收斂域較小，一般只能給出中心周期軌道附近很小幅值范圍內(nèi)的quasi周期軌道。

隨著計(jì)算機(jī)的發(fā)展，全數(shù)值方法逐漸取代半解析法，成為quasi周期軌道的主流計(jì)算方法。Howell等[13]提出全數(shù)值的微分修正算法用于計(jì)算Lissajous軌道，但是得到的并不是不變環(huán)面意義下的嚴(yán)格quasi周期軌道。Mondelo[14]提出采用Fourier級數(shù)對2維不變環(huán)面的不變曲線進(jìn)行參數(shù)化描述，以此實(shí)現(xiàn)數(shù)值計(jì)算，但是該方法計(jì)算效率較低，他們使用并行計(jì)算進(jìn)行加速。Kolemen等[15]提出多Poincaré截面法，通過截?cái)郌ourier級數(shù)把截面上閉合的軌道參數(shù)化，通過Newton迭代使得閉合的軌道滿足特定的約束條件，以此計(jì)算L2附近的quasi周期軌道，該方法具有較低的計(jì)算成本，在靠近混沌的區(qū)域依然具有較高的收斂性，不過他們并未對其他類型軌道的適用性進(jìn)行分析。Schilder等[16]提出用偏微分方程來計(jì)算quasi周期軌道的不變環(huán)面，但是求解偏微分方程較為耗費(fèi)計(jì)算資源。Olikara等[17-18]基于頻閃映射法，使用Fourier級數(shù)計(jì)算高維環(huán)面上的不變曲線，成功計(jì)算了圓型和橢圓型限制性三體系統(tǒng)的quasi-Halo軌道。Baresi等[19]綜合比較了偏微分方程法、多Poincaré截面法和頻閃映射法的特點(diǎn)。McCarthy[7]在Olikara等工作的基礎(chǔ)上，計(jì)算了等能量、等頻率比、等映射時(shí)間三種quasi周期軌道族，該方法在計(jì)算Jacobi矩陣時(shí)較為復(fù)雜。Lujan等[3]基于頻閃映射法，對地月系統(tǒng)L2點(diǎn)quasi-Halo軌道進(jìn)行單參數(shù)延拓，在頻域研究了全軌道族的特征，并得到了倍周期分岔的幾種quasi周期軌道。已有的研究工作表明，數(shù)值方法是構(gòu)造quasi周期軌道的有效手段，然而過去的研究對象一般集中于Halo軌道、Vertical軌道附近小范圍的quasi周期軌道，而且較少涉及三角平動點(diǎn)等情況，對于具有多quasi運(yùn)動方向的軌道也少有研究。

為此，本文旨在研究一種構(gòu)造大幅值quasi周期軌道的全數(shù)值方法，在Olikara等頻閃映射法[17-18]的基礎(chǔ)上，重點(diǎn)提高方法對不同類型大幅值軌道的通用性和收斂性。該方法采用多Poincaré截面把不變曲面劃分為若干段，設(shè)定若干約束，再通過微分修正的思想計(jì)算不變環(huán)面，最后利用數(shù)值延拓獲得大幅值quasi周期軌道及等映射間隔的軌道族。仿真表明，該方法不僅在計(jì)算大幅值的quasi-Vertical等共線平動點(diǎn)軌道時(shí)收斂性好，對于三角平動點(diǎn)及多對復(fù)特征值的quasi周期軌道同樣具有良好的適用性。

1 動力學(xué)模型與研究對象

1.1 圓型限制性三體模型

本文基于圓型限制性三體問題進(jìn)行建模。在會合坐標(biāo)系，航天器的動力學(xué)方程為：

(1)

其中，

(2)

1.2 quasi周期軌道

周期軌道的單值矩陣M的特征值總是成對出現(xiàn)[20]，且其中有一對為實(shí)數(shù)1，即λ2=1/λ2=1，剩下的兩對λi和1/λi(i=1,3)為復(fù)數(shù)(或?qū)崝?shù))。當(dāng)單值矩陣M具有在單位圓上的復(fù)特征值λc時(shí)，在周期軌道的附近存在quasi周期運(yùn)動，這些quasi周期軌道構(gòu)成了2維不變環(huán)面S。2維不變環(huán)面S可用兩個(gè)獨(dú)立的角度來描述，分別稱為經(jīng)度θ1和緯度θ2，相應(yīng)的頻率為ω1和ω2[18]，如圖1所示：

圖1 二維不變環(huán)面與quasi周期軌道Fig.1 Two-dimensional invariant torus and the quasi periodic orbit

Poincaré截面∑(θ1=θ1p)與2維不變環(huán)面S的交集構(gòu)成一條不變曲線Γ，不變曲線Γ上的狀態(tài)具有相同的經(jīng)度θ1p。對于一條不變曲線，quasi周期軌道從該不變曲線出發(fā)，積分時(shí)間T1=2π/ω1后，會返回至該不變曲線(如圖1所示)，T1稱為映射時(shí)間。

出發(fā)點(diǎn)與返回點(diǎn)之間存在角度差ρ，有如下關(guān)系：

(3)

圖2 中心周期軌道與quasi周期軌道的幅值Fig.2 Amplitude of central periodic orbit and the quasi periodic orbit

本文用quasi周期軌道與中心周期軌道的幅值之比來描述軌道幅值的大小。由于相關(guān)研究較為有限，沒有統(tǒng)一的標(biāo)準(zhǔn)，本文認(rèn)為Bx/Ax,By/Ay或Bz/Az三者任意一個(gè)大于1.2的quasi周期軌道屬于“大幅值”的范疇。

2 環(huán)面多步微分修正

由于大幅值quasi周期軌道所在的2維不變環(huán)面維度高、非線性強(qiáng)，對數(shù)值構(gòu)造算法的收斂性有較高要求，針對這一問題，本文在Olikara等頻閃映射法[17-18]的基礎(chǔ)上，提出映射間隔約束的環(huán)面多步微分修正算法，該算法對于大幅值的quasi周期軌道具有良好的收斂性和通用性。

在不變曲線Γ(θ1)上等間隔地選取N個(gè)緯度θ2i=2πi/N所對應(yīng)的特征點(diǎn)，其中i=1,2,…,N，這N個(gè)特征點(diǎn)的狀態(tài)稱為特征狀態(tài)，可表示為：

xi=x(t,θ1,θ2i)

(4)

圖3以N=5為例展示了特征點(diǎn)的選取原則。

圖3 等緯度地在不變曲線上選取5個(gè)特征點(diǎn)Fig.3 Five feature points are selected on the invariant curve with equal latitude

為了增強(qiáng)微分修正算法的收斂性，本文選取多個(gè)Poincaré截面，把2維不變環(huán)面劃分為若干段，相應(yīng)的形成若干條不變曲線，如圖4所示。其中，quasi周期軌道從不變曲線1出發(fā)，在其他不變曲線滿足拼接的約束條件，接著返回不變曲線1。

圖4 4個(gè)Poincaré截面和4條不變曲線Fig.4 Four Poincaré sections and four invariant curves

為了方便起見，不妨選取等間隔的M個(gè)經(jīng)度θ1j=2πj/M對應(yīng)的截面，其中j=1,2,…,M，此時(shí)每兩個(gè)不變環(huán)面的映射時(shí)間間隔為Tf=T1/M。在每個(gè)經(jīng)度θ1j所對應(yīng)的不變環(huán)面上，再以圖3的方式等間隔地選取N個(gè)緯度，則第j條不變曲線上的第i個(gè)特征狀態(tài)可表示為：

(5)

據(jù)此，本文用各不變曲線上所有特征狀態(tài)的Jacobi常數(shù)的平均值，代表quasi周期軌道的Jacobi常數(shù)Jq，即：

(6)

(7)

環(huán)面的多步微分修正需要滿足特定的約束條件：

1)拼接約束。第j=1,…,M-1條不變曲線Γ(θ1j)上的特征狀態(tài)，積分Tf后，需要與第j+1條不變曲線Γ(θ1,j+1)上的特征狀態(tài)重合。

2)返回約束。第j=M條不變曲線Γ(θ1M)上的特征狀態(tài)，積分Tf后，再經(jīng)過角度變換-ρ，需要與第1條不變曲線Γ(θ1,1)上的特征狀態(tài)重合。

3)映射間隔約束。每兩個(gè)不變環(huán)面間的映射時(shí)間間隔Tf需要等于設(shè)定的需求值Td。

映射間隔約束是附加的約束條件，所以本文稱之為“映射間隔約束的環(huán)面多步微分修正”算法。綜合以上三個(gè)約束條件，約束變量F可以表示為：

(8)

(9)

R表示角度變換-ρ的算子，即：

R·x(t,θ1j,θ2i+ρ)=x(t,θ1j,θ2i)

(10)

約束變量F(X)對自由變量X的雅可比矩陣為：

(11)

該雅可比矩陣解析形式的分量不易獲得，可以采用前向差分或中心差分等數(shù)值方法計(jì)算近似的雅可比矩陣。

由于自由變量X比約束變量F(X)的維度多1，多步微分修正的迭代過程可使用最小二乘的形式：

Xk+1=Xk-DFT(Xk)[DF(Xk)DFT(Xk)]-1F(Xk)

(12)

可通過如下的“特征值方法”[7]構(gòu)造自由變量的初始猜想X0。記λc對應(yīng)的復(fù)特征向量為vc，對于周期軌道上的特征狀態(tài)x*(t,θ1)，沿vc的方向進(jìn)行攝動δx0=εvc，則不變曲線Γ(θ1)上的任意緯度θ2對應(yīng)的特征狀態(tài)x可以近似為：

x(t,θ1,θ2)≈x*(t,θ1)+εvceθ2=

(13)

示意圖如圖5所示。

圖5 特征值方法構(gòu)造初始猜想的原理示意圖Fig.5 Schematic diagram of theprinciple of constructing the initial guess by the eigenvalue method

對于Tf，可以利用周期軌道的周期T來近似，即：

(14)

對于角度差ρ，取沿特征向量方向的攝動量δx0，根據(jù)單值矩陣M的性質(zhì)，積分一個(gè)周期后攝動量為：

(15)

通過與式(13)對比，角度差可以近似為：

(16)

3 環(huán)面的數(shù)值延拓

當(dāng)攝動量ε中等偏小時(shí)，上述特征值方法能夠給出適當(dāng)?shù)淖杂勺兞康某跏疾孪?，環(huán)面的多步微分修正算法收斂性較好。但是，當(dāng)攝動量ε超過一定的程度，自由變量的初始猜想不再準(zhǔn)確，受其制約，環(huán)面的多步微分修正算法無法保證收斂，此時(shí)可以使用延拓法計(jì)算quasi周期軌道族，實(shí)現(xiàn)計(jì)算大幅值軌道的目的。延拓法利用已有的quasi周期軌道的信息構(gòu)造新的初始猜想，故相比特征值方法更加準(zhǔn)確。

3.1 自然參數(shù)延拓法

(17)

圖6 NPC和PAC兩種算法計(jì)算quasi周期軌道族的流程圖Fig.6 Flow charts for calculating the quasi periodic orbital families of NPC and PAC

本文選擇環(huán)面的角度差ρ作為延拓變量，由于目標(biāo)是構(gòu)造等映射間隔的quasi周期軌道族，因而把映射間隔的需求量Td設(shè)為固定值。(下文所述的NPC延拓法都是以環(huán)面角度差ρ為延拓變量，并且固定映射間隔需求量Td)

3.2 偽弧長延拓法

NPC延拓法在特征曲線回折處會出現(xiàn)失效的情況，采用偽弧長延拓法(Pseudo-arclength continu-ation,PAC)進(jìn)行延拓可以解決這一問題[7,21]。偽弧長延拓法給約束變量增加了一個(gè)維度，增廣的約束變量為：

(18)

增廣的雅可比矩陣為：

(19)

自由變量初始猜想的構(gòu)造形式為：

(20)

此時(shí)，增廣約束變量與自由變量的維度相同，微分修正的迭代過程為：

Xk+1=Xk-DG-1(Xk)G(Xk)

(21)

為了實(shí)現(xiàn)計(jì)算等映射間隔的quasi周期軌道族，偽弧長延拓法的計(jì)算過程中，Td始終是固定值。使用偽弧長延拓法計(jì)算quasi周期軌道族完整的流程如圖6(b)所示。

NPC和PAC兩種延拓法均固定映射時(shí)間間隔的需求值Td，因此，族內(nèi)每條軌道都具有相同的映射時(shí)間T1。NPC的優(yōu)點(diǎn)是延拓步長s有明確的物理意義，代表族參數(shù)ρ的改變量，而PAC的延拓步長s沒有明確的物理意義。NPC的缺點(diǎn)是將在族參數(shù)ρ達(dá)到極值時(shí)延拓失敗，而PAC的延拓進(jìn)程則不會受到極值的制約。綜上所述，兩種延拓法都可以計(jì)算quasi周期軌道族，也都有各自的優(yōu)缺點(diǎn)，可以互補(bǔ)使用。

4 數(shù)值仿真與分析

本節(jié)將以地月CRTBP系統(tǒng)的L1點(diǎn)Vertical軌道、L4點(diǎn)Vertical軌道、L4點(diǎn)Axial軌道、遠(yuǎn)距離逆行軌道(Distant retrograde orbit,DRO)等為中心周期軌道，利用映射間隔約束的環(huán)面多步微分修正算法，考慮NPC與PAC兩種延拓法，設(shè)計(jì)各中心周期軌道對應(yīng)的quasi周期軌道族。地月CRTBP系統(tǒng)的參數(shù)及中心周期軌道的初值見表1。

表1 地月CRTBP系統(tǒng)的參數(shù)及中心周期軌道的初值Table 1 Parameters of the Earth-Moon CRTBP and initial states of the central periodic orbit

4.1 quasi-Vertical軌道族算例

NPC與PAC具有不同的延拓機(jī)理，因此計(jì)算出的軌道族將存在差異。本部分以L1點(diǎn)的Vertical軌道為中心周期軌道，在相同仿真參數(shù)的條件下，對比NPC與PAC的計(jì)算結(jié)果，相關(guān)仿真參數(shù)見表2。

表2 quasi-Vertical軌道族的仿真參數(shù)Table 2 Parameters of the quasi-Vertical orbital family

表2中，NPC的延拓步數(shù)80是能保證算法收斂的最大步數(shù)，即步數(shù)再增加，算法將發(fā)散；PAC的延拓步數(shù)38再增加，計(jì)算出的軌道與上一個(gè)步長無明顯差別，因此人為地不再增加步數(shù)，結(jié)束延拓進(jìn)程。圖7給出了NPC與PAC兩種方法延拓過程的Poincaré截面，每個(gè)黑色實(shí)心圓點(diǎn)代表quasi-Vertical軌道穿過xy平面一次?？梢钥闯?，隨著延拓的進(jìn)行，截面的y向幅值逐漸增大。最后一個(gè)延拓步長所對應(yīng)環(huán)面稱為族末端環(huán)面，PAC族末端環(huán)面的By/Ay為47.4，而NPC族末端環(huán)面這一值為39.2，即PAC軌道族末端環(huán)面的y向幅值比NPC軌道族更大。

圖7 兩種方法延拓過程的Poincaré截面Fig.7 Poincaré sections of the two continuation processes

兩種延拓方法給出的都是等映射間隔的quasi-Vertical軌道族，然而Poincaré截面與族末端環(huán)面呈現(xiàn)差異的原因是，NPC在族參數(shù)ρ減小到一定程度時(shí)延拓失敗，而PAC不受族參數(shù)ρ的制約，故給出的是完整的軌道族。圖8給出了兩種方法的軌道族的Jacobi常數(shù)、映射時(shí)間間隔Tf、ρ等特征參數(shù)之間的差值與幅值Bx之間的關(guān)系，當(dāng)幅值Bx增大到0.0475時(shí)，NPC算法終止。Tf的差值始終為0，說明兩個(gè)軌道族具有相同的映射時(shí)間間隔，其他兩個(gè)特征參數(shù)差值的量級非常小，綜合起來說明兩個(gè)軌道族的自由變量基本一致，兩種方法給出的是同一個(gè)軌道族，區(qū)別在于PAC的是完整的軌道族，而NPC的是部分的軌道族。

圖8 兩種軌道族特征參數(shù)之間的差值與幅值的關(guān)系Fig.8 Relation between difference of characteristic parameters and amplitude for the two orbital families

4.2 三角平動點(diǎn)quasi周期軌道族算例

前面分析了兩種延拓法對于L1點(diǎn)的Vertical軌道的適用性，本部分以L4點(diǎn)的Vertical軌道和L4點(diǎn)的Axial軌道為中心周期軌道，采用能給出完整軌道族的PAC算法進(jìn)行延拓，驗(yàn)證算法對三角平動點(diǎn)quasi周期軌道族的適用性，相關(guān)仿真參數(shù)見表3。

表3 三角平動點(diǎn)軌道族的仿真參數(shù)Table 3 Parameters ofthe triangular libration point orbital family

圖9(a)是L4點(diǎn)quasi-Vertical軌道族末端環(huán)面，Jacobi常數(shù)為2.6183，x向幅值比為8.34。圖9(b)是L4點(diǎn)quasi-Axial軌道族末端環(huán)面，Jacobi常數(shù)為2.1744，x向幅值比為1.37。

圖9 quasi-Vertical和quasi-Axial軌道族末端環(huán)面Fig.9 Terminal invariant torus of quasi-Vertical and quasi-Axial orbital families

可以看出，對于L4點(diǎn)quasi-Vertical和quasi-Axial軌道族，族末端軌道與中心周期軌道的x向幅值比都超過了1.2，說明算法對計(jì)算三角平動點(diǎn)大幅值quasi周期軌道具有良好的適用性。

4.3 多對復(fù)特征值的quasi周期軌道族算例

以上算例所選取的中心周期軌道，只有一對在單位圓上的復(fù)特征值，所以只有一個(gè)quasi周期運(yùn)動的方向，單對復(fù)特征值也是以往的研究所較多涉及的。但是，CRTBP還存在具有兩對在單位圓上的復(fù)特征值的軌道，比如DRO軌道等，這些軌道相應(yīng)地具有兩個(gè)quasi周期運(yùn)動的方向。

近年來，DRO軌道因其良好的穩(wěn)定性，在月球空間站、小行星防御等空間任務(wù)中得到了廣泛的研究[22-23]。本算例以目前研究較多的2∶1共振DRO軌道為中心周期軌道，采用PAC算法進(jìn)行quasi周期軌道族的延拓，分析其兩個(gè)quasi周期運(yùn)動方向的特性。由于PAC算法生成的是等映射時(shí)間間隔的軌道族，因此族內(nèi)每條軌道的映射時(shí)間T1都等于中心DRO軌道的周期，所以生成的quasi-DRO軌道依然具有2∶1共振特性。相關(guān)仿真參數(shù)及2∶1共振DRO軌道的6個(gè)特征值見表4。

表4 2∶1共振quasi-DRO軌道族的仿真參數(shù)Table 4 Parameters of the 2∶1 resonance quasi-DRO families

對于特征值λ3和1/λ3，quasi周期運(yùn)動的方向在xy平面內(nèi)，相應(yīng)地，quasi周期軌道也在xy平面內(nèi)。圖10(a)是xy平面內(nèi)quasi-DRO軌道族的末端環(huán)面和軌道，y方向的幅值比為1.28。

對于特征值λ1和1/λ1，quasi周期運(yùn)動的方向在xy平面外，圖10(b)是xy平面外quasi-DRO軌道族末端環(huán)面和軌道，它是三維的軌道，即z方向具有quasi周期運(yùn)動的分量。因?yàn)橹行腄RO軌道不具有z向幅值，不妨用Bz/Ax來描述幅值的大小，該條軌道的這一參數(shù)為0.35。由于具有z方向的運(yùn)動分量，相比傳統(tǒng)的平面DRO軌道，它作為任務(wù)軌道更具有天然的優(yōu)勢，其一，軌道仍然具有2∶1的共振特性；其二，通過選取適當(dāng)?shù)膠向幅值等參數(shù)，可以規(guī)避日食及運(yùn)行到月球背面時(shí)的通信遮擋問題[24]。

圖10 面內(nèi)和面外quasi周期運(yùn)動方向的軌道族末端環(huán)面Fig.10 Terminal invariant torus of planar and out-of-plane quasi-periodic motions

綜合比較各種quasi周期軌道族末端環(huán)面3個(gè)方向的幅值比，見表5。幾乎每種周期軌道族末端環(huán)面都有1個(gè)方向的幅值比大于1.2，說明本文的數(shù)值構(gòu)造算法能夠適應(yīng)較多類型的大幅值quasi周期軌道。

表5 quasi周期軌道族末端環(huán)面3個(gè)方向的幅值比Table 5 Amplitude ratio of 3 directions for terminal invariant torus of quasi-periodic orbital family

5 結(jié) 論

本文提出了一種利用映射間隔約束的環(huán)面多重微分修正算法，結(jié)合NPC或PAC計(jì)算大幅值quasi周期軌道及軌道族的策略。本文方法能夠有效地解決共線平動點(diǎn)、三角平動點(diǎn)、多對復(fù)特征值等情況下大幅值quasi周期軌道及軌道族的數(shù)值計(jì)算問題。