圓型限制性三體問題相對運動解析研究

周 敬，胡 軍，張 斌

(1. 北京控制工程研究所， 北京 100190; 2. 空間智能控制技術國防科技重點實驗室， 北京 100190)

0 引 言

近年來，隨著航天技術的快速發展，深空探測受到世界各航天大國與組織越來越多的重視。按照我國航天事業發展的指南，未來航天的三個重點方向之一就是深空探測[1]。因此，與深空探測密切相關的各項基礎和應用研究便成為了我國航天事業中的一個熱點研究方向。

不同于經典二體問題下的近地空間航天活動，深空探測所涉及的運動模型更多的是三體問題。由于動力學本質上的區別，三體問題相對二體問題更加復雜。在三體問題的研究中，最簡單、應用最廣泛的模型是圓型限制性三體問題(Circular restricted three-body problem，CRTBP)。在CRTBP中存在5個動平衡點，也稱為平動點(Libration points)或拉格朗日點(Lagrangian points)。這5個平動點相對主天體的位置是固定的，其中L1,L2和L3與兩主天體共線(L1位于主天體之間的連線上，L2和L3位于延長線上)，稱為共線平動點；L4和L5均與兩個主天體構成等邊三角形，稱為三角平動點。平動點(尤其是共線平動點)及其附近的周期/擬周期軌道所具有的獨特位置優勢使其在深空探測任務設計中占據著重要的地位，已經成為宇宙觀測、天文研究的理想場所和星際高速公路(Interplanetary superhighway，IPS)的門戶站。

為了實現重返月球，美國航天局和其他太空研究機構已經考慮在地-月L1點建立空間站，以支持載人登月、載人火星探測甚至宇宙航行。與近地軌道空間站一樣，平動點空間站的建設和運行同樣離不開交會對接(Rendezvous and docking，RVD)技術的支持，而RVD技術的基礎則是相對運動。此外，與近地軌道航天任務發展趨勢一樣，深空探測任務發展趨勢也必將由單顆航天器獨立工作向多顆航天器協同工作轉變，如TPF、Darwin、MAXIM計劃等，而多顆航天器協同工作所采用的編隊飛行技術更需要相對運動研究作為基石。因此，對三體問題平動點附近相對運動的研究兼具理論研究和工程應用價值。

目前，涉及三體問題相對運動的研究工作主要有：侯錫云等[2]系統地闡述了限制性三體問題中共線平動點的動力學特征，給出了這類平動點附近的周期軌道和擬周期軌道的計算方法；徐明[3]詳細介紹了平動點軌道的發展歷史，并深入剖析了平動點附近的相空間結構；文獻[4]通過建立改進型龐加萊截面圖，同時結合狀態轉移矩陣和打靶法，建立了一種計算三體問題周期軌道的新方法；Richardson等[5]將推導Halo軌道三階解析解的思想應用到編隊飛行，得到了非線性相對運動方程的三階展開式；Gómez等[6]將其關于Lissajous軌道的Fourier算法等研究成果推廣到相對運動動力學，得到了平動點相對運動的高階展開式；文獻[7-9]基于Richardson三階解析解研究了Halo軌道上編隊飛行航天器的相對運動構型特性；Howell等[10]運用瞄準法獲得了期望的日-地平動點附近航天器的非自然編隊構型；文獻[11]基于單值矩陣和Fourier級數，將編隊設計問題從微分方程轉化為代數方程，分別設計了沿Halo軌道的長期和短期自然松散編隊構型；Luo等[12]提出了一種新型的、基于遍歷Poincaré映射的數值搜索方法，獲得了CRTBP擬Halo軌道附近的自然編隊飛行構型；Luquett等[13]以全非線性形式給出了限制性三體問題下的自然相對運動動力學模型，然后將其線性化并轉化為狀態空間的形式，獲得了相對運動的近似解析解；文獻[14]應用Lindstedt-Poincaré方法和多項式展開的方法，構造了限制性三體問題平動點附近周期相對運動構型的解析解；Cont等[15]通過線性化的方法對限制性三體問題平動點附近的相對運動和交會問題進行解析分析，并與數值仿真結果進行了對比驗證；Ferrari等[16]通過分析初始構型和關鍵參數的方法，研究了圓型和橢圓型限制性三體問題下共線平動點周期軌道附近無控編隊飛行的演化過程；黃海濱等[17]運用Legendre偽譜法和協同進化粒子群(CPSO)方法研究了深空環境下衛星編隊飛行隊形重構實時重規劃問題；Arnot等[18]在CRTBP模型下，應用最小位置反饋方法構造了新的穩定的、調頻的相對運動軌道；Case[19]運用線性化相對運動動力學和非線性微分修正方法實現了CRTBP下平動點軌道的交會；Prioroc等[20]利用極值搜索方法對時延反饋機構和PD控制器的增益進行優化，實現了小推力下的日-地CRTBP平動點附近的編隊飛行；Jung等[21]運用開關-Hamiltonian保結構控制方法實現了地-月L2點Halo軌道上的編隊飛行控制；文獻[22]設計了基于輸出調節理論的位置保持控制器，實現了CRTBP平動點附近航天器的位置保持和編隊飛行控制；Prioroc等[23]運用非線性控制的方法研究了Halo軌道編隊飛行的相對運動控制問題；張楷田等[24]基于CRTBP模型推導了航天器編隊飛行的非線性動力學方程，然后應用線性自抗擾控制技術實現了編隊飛行控制；姜春生等[25]結合擴張觀測器和自抗擾控制方法，實現了對日-地平動點附近的航天器編隊飛行的控制；Yin等[26]運用自抗擾控制方法解決了考慮不確定性和攝動情況下的日-地CRTBP下L2平動點附近的多航天器編隊飛行控制問題；Wang等[27]結合反饋控制和分布式自適應同步方法，實現了速度不可測和質量不確定情況下平動點軌道附近的多航天器編隊飛行；王峰等[28]建立日-地L2平動點編隊飛行相對運動QLPV動力學模型，并提出一種改進的多項式特征結構配置方法實現日-地L2平動點編隊飛行高精度相對位置保持；文獻[29]使用多重打靶修正方法改進了中心航天器的擬Halo參考軌道，并使用線性化的變分方程來進行位置保持控制，最終運用Hamiltonian保結構控制方法實現了L1點擬周期軌道附近編隊飛行的穩定控制。

從上述研究中可以發現，有關三體問題下的相對運動研究多以編隊飛行為研究背景，通過將軌道絕對運動的近似解析解直接作差，得到相對運動的近似解析解；或者將相對運動動力學方程線性化，然后運用線性系統理論對其求解，獲得相對運動的近似解析解；亦或者將相對運動問題具化為實現某一具體任務而引申出的控制問題加以解決。這三種方式在多數情況下可以滿足任務需求，但不便于進一步獲得三體問題下相對運動的本質規律。基于此，本文借鑒Floquet理論的研究思想，在對單值矩陣進行分析的基礎上，從更本質的角度出發，對三體問題中的共線平動點附近周期/擬周期軌道下的相對運動進行研究。

1 圓型限制性三體問題

限制性三體問題研究的是一個小天體或航天器在兩個大天體(主天體)引力場中的運動情況，其中小天體或航天器的質量相比于兩主天體可忽略不計，不會對兩主天體之間的相互運動產生影響。當主天體圍繞公共質心作圓軌道運動時，即為圓型限制性三體問題(CRTBP)。為方便起見，本文以地-月系統CRTBP中的L1點作為研究對象開展相對運動研究。本文研究方法同樣適用于其他三體系統共線平動點附近的相對運動。

1.1 動力學模型

為方便描述航天器運動和簡化計算，在CRTBP中通常需要對相關物理量進行無量綱化處理，并以時間作為獨立變量。相應的質量[M]、長度[L]和時間[T]的歸一化單位取為

(1)

式中:m1，m2為兩主天體質量；L12為兩主天體之間的距離；G為萬有引力常數。

為了更直觀、更清晰地描述航天器在CRTBP中的運動，通常選擇在會合坐標系(也稱為質心旋轉坐標系)下進行研究。如圖1所示，會合坐標系的原點位于兩主天體的公共質心，x軸由較大主天體指向較小主天體，z軸與主天體系統角動量方向平行，y軸滿足右手坐標系定則。

圖1 地-月系統CRTBP下的會合坐標系和L1會合坐標系

(2)

式中：Ω表示偽勢能函數。

(3)

式中：X1表示在會合坐標系下L1點距離主天體公共質心的距離，γ1為L1點與最近主天體之間的距離。

因此，L1合坐標系下航天器運動的動力學方程為

(4)

1.2 近似解析解

在會合坐標系下，假設[δX,δY,δZ]T表示航天器相對L1平動點的位置矢量。根據文獻[25]，共線平動點附近航天器運動的近似解析解為

(5)

式中：Ai(i=1,…,6)是由初始條件確定的積分常數，并在本文中作為Lissajous軌道和Lyapunov軌道的軌道參數，參數k1,k2以及λ1,λ3,λ5的詳細計算公式可以參考文獻[30]。

由于Halo軌道是共線平動點附近一類特殊的三維周期軌道，且常作為深空探測航天器的工作軌道，具有更加重要的研究和應用價值。因此，Richardson[31]基于Lindstedt-Poincaré方法獲得了Halo軌道的三階近似解析解，如式(6)所示，相關參數的含義及計算過程可以參考文獻[31]。

(6)

此外，關于Halo軌道三階近似解析解的改進可以參考文獻[32]。

1.3 微分修正

由于CRTBP具有強非線性、不穩定性和“混沌”性，第1.2節獲得的近似解析解在真實動力學模型下積分一段時間后會很快發散。為了獲得精確、收斂的周期軌道，可以采用微分修正方法。有關微分修正方法的詳細介紹可以參考文獻[1]。微分修正方法主要利用周期/擬周期軌道的對稱性來修正近似解析解提供的軌道初值，使修正后的初值進行數值積分后生成的周期軌道可以自然維持多圈而不發散。下面分別介紹微分修正方法在Halo軌道、Lyapunov軌道和Lissajous軌道上的具體應用過程。

1.3.1Halo軌道

由于Halo軌道關于XZ平面對稱，根據參考文獻[1]和近似解析解(6)，微分修正的初始狀態x0設置為

(7)

(8)

式中：Φ(i,j)表示矩陣Φ中第i行，第j列元素。式(8)即為Halo軌道初值的微分修正公式。詳細修正過程可以參考文獻[1]。

1.3.2Lyapunov軌道

由于Lyapunov軌道是一種XY平面軌道且關于x軸對稱，根據參考文獻[1]和近似解析解(5)，微分修正的初始狀態x0設置為

(9)

(10)

式(10)即為Lyapunov軌道初值的微分修正公式。詳細修正過程可以參考文獻[1]。

1.3.3Lissajous軌道

Lissajous軌道是一種擬周期軌道，而非嚴格意義上的周期軌道，因此其微分修正過程較Halo軌道和Lyapunov軌道更為復雜。Howell等[33]提出使用兩層微分修正方法來求解CRTBP下的Lissajous軌道，考慮到兩層微分修正方法的復雜性，本文提出一種較簡單的修正方法。鑒于Lissajous軌道在XY平面內的運動與Z軸的運動是解耦的，且兩部分運動的周期也不相同，因此可以在XY平面和Z軸上分別、同時進行微分修正，最終獲得精確的Lissajous軌道初值。

首先考慮XY平面內的運動，Lissajous軌道在XY平面內的運動軌道實際上為Lyapunov軌道，因此在XY平面上的微分修正可以沿用Lyapunov軌道的修正方法。

根據近似解析解(5)，微分修正的初始狀態x0設置為

(11)

對應的第一個期望狀態設置為

(12)

經過時間t1后的真實狀態x1為

(13)

則真實狀態x1和期望狀態xd1之差δx1為

(14)

(15)

然后，考慮Z軸上的運動，對應的第二個期望狀態設置為

(16)

經過時間t2后的真實狀態x2為

(17)

則真實狀態x2和期望狀態xd2之差δx2為

(18)

對式(18)進行反解，得到修正量δz0和δt2的計算公式

(19)

式(15)和式(19)即為Lissajous軌道初值的微分修正公式。上述兩部分微分修正同時進行，便可以獲得精確的Lissajous軌道初值。

2 相對運動

相對運動通常涉及兩航天器，即目標航天器和追蹤航天器。假設目標航天器的運行軌道為地-月L1點附近的一條周期/擬周期軌道，追蹤航天器運行在目標航天器附近。故本文相對運動特指，以某個Halo/Lyapunov/Lissajous軌道為參考軌道，研究航天器在此參考軌道附近相對此參考軌道的運動。

圖2 相對運動示意圖

本文研究相對運動的目的是獲得與近地軌道CW方程類似的CRTBP下相對運動的近似解析解。下面以Halo軌道作為目標航天器的運行軌道為例，介紹相對運動的研究方法。

在CRTBP中，單值矩陣(Monodromy Matrix)M是指一個軌道周期T對應的周期軌道的狀態轉移矩陣(State transition matrix,STM)Φ(t0+t,x0)，即

M=Φ(t0+T,x0)

(20)

雖然單值矩陣沒有精確的解析表達式，但可以通過數值積分的方式獲得精確的數值解，詳細的計算方法可以參考文獻[25]。

單值矩陣共存在6個特征值，表示為(λ1,λ2,λ3,λ4,λ5,λ6)，可以劃分為三對，即

(21)

第一對特征值(λ1,λ2)是互為倒數的正實數，與參考軌道的雙曲特性有關。λ1>1是絕對值最大的特征值，其對應的特征向量代表著參考軌道的發散方向，因此λ1及其對應的特征向量在一定程度上代表著參考軌道的不穩定特性。與之相反，由于λ2=1/λ1<1，故λ2及其對應的特征向量在一定程度上表示參考軌道的穩定特性。

第二對特征值(λ3,λ4)的數值均為1，在一定程度上代表著參考軌道的既不發散也不收斂特性，即參考軌道的中性特性。

第三對特征值(λ5,λ6)為模值為1的共軛復數，在一定程度上代表著參考軌道的中心流形，即參考軌道的周期特性。

綜上所述，單值矩陣的六個特征值和對應的特征向量提供了一種具有幾何意義的表示，即不同的特征值和對應的特征向量可以代表相對運動中具有不同特性的分量模態，如不穩定模態、穩定模態、中性模態以及周期模態，這也是Floquet理論[34]的核心思想。因此，通過對單值矩陣的分析，可以建立由一組獨立基向量所形成的一個局部坐標系統，并以此局部坐標系統來描述相對運動。借鑒Floquet理論的思想，這組基向量可以由單值矩陣的特征向量形成，而這組基向量也就是本文提出的相對運動模態。

假設相對運動模態矩陣E(t)由6個相對運動模態組成

E(t)=[e1(t),e2(t),e3(t),e4(t),e5(t),e6(t)]

(22)

式中：ei(t),i=1,2,3,…,6即為6個相對運動模態，分別與單值矩陣的6個特征值和特征向量相對應。

通過對單值矩陣的分析，可以令Halo軌道上當前狀態對應的單值矩陣的特征向量組成相對運動模態矩陣E(t)的初始值E0，即

E0=[v1,v2,v3,v4,real(v5),imag(v6)]

(23)

式中：vi(i=1,…,6)即為單值矩陣的6個特征向量，可以通過求解單值矩陣獲得，real(v5)表示特征向量的實部，imag(v6)表示特征向量的虛部。不難看出rank(E0)=6，即這組基向量是線性無關的，可以用來描述6維空間中的任意向量，即可以用來描述任意初始狀態的相對運動。

根據Halo軌道的性質[30]以及狀態轉移矩陣性質，相對運動模態矩陣滿足如下關系

E(t)=Φ(t)E0

(24)

在獲得相對運動模態變化規律E(t)的基礎上，便可以將相對運動看做上述6個相對運動模態的線性組合。假設追蹤航天器和目標航天器的運動狀態分別記為x′和x，則兩者之間的差值即為航天器間的相對運動δx，

(25)

將相對運動δx表示成相對運動模態的線性組合形式，即

(26)

式中：C(t)為各個模態在相對運動中所占的權值矩陣，不難發現其為一個對角矩陣，即

C(t)=diag(c1(t),c2(t),c3(t),c4(t),

c5(t),c6(t))

(27)

而相對運動的初始狀態δx0可以表示為

δx0=C0E0

(28)

根據文獻[35]，由于相對運動的尺度相對參考軌道尺度來說為一小量，因此可以將追蹤航天器相對于目標航天器的相對運動δx看做目標航天器狀態x的一種攝動，則根據狀態轉移矩陣的性質，有

δx(t)=Φ(t)δx0

(29)

將式(29)與式(24)和式(28)聯立，有

(30)

從式(30)可以得到

C(t)=C0

(31)

即權值矩陣C(t)為一常值矩陣。

這樣，便獲得了Halo軌道為參考軌道的相對運動近似解析解

δx(t)=C0Φ(t)E0

(32)

由于任意時刻的狀態轉移矩陣Φ(t)可以事先獲得，即作為已知量，依據此近似解析解，便可以獲得任意時刻的相對運動狀態。

由于單值矩陣M是非線性系統周期解的一種共有的特性，而Lyapunov軌道同Halo軌道一樣均為非線性CRTBP動力學方程的周期解；此外，雖然Lissajous軌道僅為擬周期軌道，但實際上其XY平面內的運動周期與Z軸上的運動周期相差并不大，在一定精度條件下也可以視為周期軌道。因此，上述Halo軌道下的相對運動分析方法同樣適用于Lyapunov軌道和Lissajous軌道下的相對運動。

根據上述分析，本文提出的CRTBP共線平動點附近相對運動近似解析解的適用條件如下：

1)目標航天器運行的參考軌道為Halo軌道、Lyapunov軌道時的相對運動，以及精度要求不高的情況下，也適用于Lissajous軌道為參考軌道的相對運動。

2)相對運動相比參考軌道的尺度較小，即相對運動可以作為參考軌道的攝動來處理。

對比參考文獻[11]，本文提出的相對運動研究方法不僅適用于Halo軌道，也適用于Lyapunov軌道和Lissajous軌道，是一種在Floquet理論基礎上改進的、通用的、基于單值矩陣的CRTBP共線平動點附近周期/擬周期軌道下相對運動的解析研究方法。

3 仿真校驗

為校驗本文所提出的基于單值矩陣的CRTBP共線平動點附近周期/擬周期軌道下的相對運動解析研究方法的正確性，本節分別以地-月系統L1點附近的Halo軌道、Lissajous軌道和Lyapunov軌道作為目標航天器的參考軌道，研究追蹤航天器與目標航天器之間的相對運動。

3.1 Halo軌道

假設目標航天器和追蹤航天器均運行在各自的Halo軌道上，相應的Z軸振幅分別記Az1和Az2，相對運動時間記為tTOF，xT0為目標航天器的初始狀態，xC0為追蹤航天器的初始狀態，以上各參數在L1會合坐標系下的無量綱單位(Dimensionless unit, DU)取值如表1所示。

Halo軌道相對運動的6個模態的初值以及一個軌道周期內的變化趨勢分別如表2和圖3所示。從仿真結果可以看出，模態1呈發散趨勢，證明了其為不穩定模態；模態2呈收斂趨勢，證明了其為穩定模態；模態3和模態4呈既不發散也不收斂趨勢，證明了其為中性模態；模態5和模態6呈周期變化趨勢，證明了其為周期模態。

表1 Halo軌道相對運動相關參數

表2 Halo軌道相對運動6個模態初始值

圖3 Halo軌道相對運動6個模態位置空間的變化規律

圖4表示真實Halo軌道相對運動與反演Halo軌道相對運動的對比。真實相對運動為追蹤航天器的軌道數據與目標航天器的軌道數據之差；反演相對運動是由相對運動近似解析解式(32)生成。從圖4可以看出，反演相對運動與真實相對運動相一致，雖然由于線性化的原因會不可避免地出現一定的誤差，但依然可以說明本方法的有效性。

圖4 Halo軌道相對運動仿真結果

3.2 Lissajous軌道

假設目標航天器和追蹤航天器均運行在各自的Lissajous軌道上，相應的軌道參數[A1,A2,A3,A4,A5,A6]分別記A1和A2，相對運動時間記為tTOF，xT0為目標航天器的初始狀態，xC0為追蹤航天器的初始狀態，以上各參數在L1會合坐標系下的無量綱單位取值如表3所示。

表3 Lissajous軌道相對運動相關參數

Lissajous軌道相對運動的6個模態的初值以及一個周期內的變化趨勢分別如表4和圖6所示。從仿真結果可以看出，6個相對運動模態同樣均與其自身的性質相符合。

圖5表示真實Lissajous軌道相對運動與反演Lissajous軌道相對運動的對比。從圖5可以看出，反演相對運動與真實相對運動相一致，誤差很小，說明了本方法的有效性。

圖5 Lissajous軌道相對運動仿真結果

表4 Lissajous軌道相對運動6個模態初始值

3.3 Lyapunov軌道

假設目標航天器和追蹤航天器均運行在各自的Lyapunov軌道上，相應的軌道參數[A1,A2,A3,A4,A5,A6]分別記A1和A2，相對運動時間記為tTOF，xT0為目標航天器的初始狀態，xC0為追蹤航天器的初始狀態，以上各參數在L1會合坐標系下的無量綱單位取值如表5所示。

Lyapunov軌道相對運動的6個模態的初值以及一個周期內的變化趨勢分別如表6和圖7所示。從仿真結果可以看出，6個相對運動模態同樣均與其自身的性質相符合。

圖8表示真實Lyapunov軌道相對運動與反演Lyapunov軌道相對運動的對比。從圖8可以看出，反演相對運動與真實相對運動相一致，誤差很小，說明了本方法的有效性。

表5 Lyapunov軌道相對運動相關參數

表6 Lyapunov軌道相對運動6個模態初始值

圖7 Lyapunov軌道相對運動6個模態位置子空間的變化規律

圖8 Lyapunov軌道相對運動仿真結果

4 結 論

本文針對圓型限制性三體問題共線平動點附近周期/擬周期軌道下的相對運動問題，給出了一種通用的、基于單值矩陣分析的相對運動近似解析解。該方法通過對三體問題下周期/擬周期軌道單值矩陣的分析，依據單值矩陣的特征值和對應的特征向量建立了6個相對運動模態，具體分為不穩定模態、穩定模態、中性模態以及周期模態，然后將航天器間的相對運動表示為上述模態的線性組合，最終獲得了相對運動的近似解析解。最后，以地-月系統L1點為研究對象，對參考軌道分別為Halo軌道、Lissajous軌道和Lyapunov軌道時的相對運動模態和相對運動進行了仿真分析，說明了方法的有效性。