周 敬，胡 軍，張 斌

(1. 北京控制工程研究所， 北京 100190; 2. 空間智能控制技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室， 北京 100190)

0 引 言

近年來(lái)，隨著航天技術(shù)的快速發(fā)展，深空探測(cè)受到世界各航天大國(guó)與組織越來(lái)越多的重視。按照我國(guó)航天事業(yè)發(fā)展的指南，未來(lái)航天的三個(gè)重點(diǎn)方向之一就是深空探測(cè)[1]。因此，與深空探測(cè)密切相關(guān)的各項(xiàng)基礎(chǔ)和應(yīng)用研究便成為了我國(guó)航天事業(yè)中的一個(gè)熱點(diǎn)研究方向。

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

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

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

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

1 圓型限制性三體問(wèn)題

限制性三體問(wèn)題研究的是一個(gè)小天體或航天器在兩個(gè)大天體(主天體)引力場(chǎng)中的運(yùn)動(dòng)情況，其中小天體或航天器的質(zhì)量相比于兩主天體可忽略不計(jì)，不會(huì)對(duì)兩主天體之間的相互運(yùn)動(dòng)產(chǎn)生影響。當(dāng)主天體圍繞公共質(zhì)心作圓軌道運(yùn)動(dòng)時(shí)，即為圓型限制性三體問(wèn)題(CRTBP)。為方便起見(jiàn)，本文以地-月系統(tǒng)CRTBP中的L1點(diǎn)作為研究對(duì)象開(kāi)展相對(duì)運(yùn)動(dòng)研究。本文研究方法同樣適用于其他三體系統(tǒng)共線平動(dòng)點(diǎn)附近的相對(duì)運(yùn)動(dòng)。

1.1 動(dòng)力學(xué)模型

為方便描述航天器運(yùn)動(dòng)和簡(jiǎn)化計(jì)算，在CRTBP中通常需要對(duì)相關(guān)物理量進(jìn)行無(wú)量綱化處理，并以時(shí)間作為獨(dú)立變量。相應(yīng)的質(zhì)量[M]、長(zhǎng)度[L]和時(shí)間[T]的歸一化單位取為

(1)

式中:m1，m2為兩主天體質(zhì)量；L12為兩主天體之間的距離；G為萬(wàn)有引力常數(shù)。

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

圖1 地-月系統(tǒng)CRTBP下的會(huì)合坐標(biāo)系和L1會(huì)合坐標(biāo)系

(2)

式中：Ω表示偽勢(shì)能函數(shù)。

(3)

式中：X1表示在會(huì)合坐標(biāo)系下L1點(diǎn)距離主天體公共質(zhì)心的距離，γ1為L(zhǎng)1點(diǎn)與最近主天體之間的距離。

因此，L1合坐標(biāo)系下航天器運(yùn)動(dòng)的動(dòng)力學(xué)方程為

(4)

1.2 近似解析解

在會(huì)合坐標(biāo)系下，假設(shè)[δX,δY,δZ]T表示航天器相對(duì)L1平動(dòng)點(diǎn)的位置矢量。根據(jù)文獻(xiàn)[25]，共線平動(dòng)點(diǎn)附近航天器運(yùn)動(dòng)的近似解析解為

(5)

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

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

(6)

此外，關(guān)于Halo軌道三階近似解析解的改進(jìn)可以參考文獻(xiàn)[32]。

1.3 微分修正

由于CRTBP具有強(qiáng)非線性、不穩(wěn)定性和“混沌”性，第1.2節(jié)獲得的近似解析解在真實(shí)動(dòng)力學(xué)模型下積分一段時(shí)間后會(huì)很快發(fā)散。為了獲得精確、收斂的周期軌道，可以采用微分修正方法。有關(guān)微分修正方法的詳細(xì)介紹可以參考文獻(xiàn)[1]。微分修正方法主要利用周期/擬周期軌道的對(duì)稱性來(lái)修正近似解析解提供的軌道初值，使修正后的初值進(jìn)行數(shù)值積分后生成的周期軌道可以自然維持多圈而不發(fā)散。下面分別介紹微分修正方法在Halo軌道、Lyapunov軌道和Lissajous軌道上的具體應(yīng)用過(guò)程。

1.3.1Halo軌道

由于Halo軌道關(guān)于XZ平面對(duì)稱，根據(jù)參考文獻(xiàn)[1]和近似解析解(6)，微分修正的初始狀態(tài)x0設(shè)置為

(7)

(8)

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

1.3.2Lyapunov軌道

由于Lyapunov軌道是一種XY平面軌道且關(guān)于x軸對(duì)稱，根據(jù)參考文獻(xiàn)[1]和近似解析解(5)，微分修正的初始狀態(tài)x0設(shè)置為

(9)

(10)

式(10)即為L(zhǎng)yapunov軌道初值的微分修正公式。詳細(xì)修正過(guò)程可以參考文獻(xiàn)[1]。

1.3.3Lissajous軌道

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

首先考慮XY平面內(nèi)的運(yùn)動(dòng)，Lissajous軌道在XY平面內(nèi)的運(yùn)動(dòng)軌道實(shí)際上為L(zhǎng)yapunov軌道，因此在XY平面上的微分修正可以沿用Lyapunov軌道的修正方法。

根據(jù)近似解析解(5)，微分修正的初始狀態(tài)x0設(shè)置為

(11)

對(duì)應(yīng)的第一個(gè)期望狀態(tài)設(shè)置為

(12)

經(jīng)過(guò)時(shí)間t1后的真實(shí)狀態(tài)x1為

(13)

則真實(shí)狀態(tài)x1和期望狀態(tài)xd1之差δx1為

(14)

(15)

然后，考慮Z軸上的運(yùn)動(dòng)，對(duì)應(yīng)的第二個(gè)期望狀態(tài)設(shè)置為

(16)

經(jīng)過(guò)時(shí)間t2后的真實(shí)狀態(tài)x2為

(17)

則真實(shí)狀態(tài)x2和期望狀態(tài)xd2之差δx2為

(18)

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

(19)

式(15)和式(19)即為L(zhǎng)issajous軌道初值的微分修正公式。上述兩部分微分修正同時(shí)進(jìn)行，便可以獲得精確的Lissajous軌道初值。

2 相對(duì)運(yùn)動(dòng)

相對(duì)運(yùn)動(dòng)通常涉及兩航天器，即目標(biāo)航天器和追蹤航天器。假設(shè)目標(biāo)航天器的運(yùn)行軌道為地-月L1點(diǎn)附近的一條周期/擬周期軌道，追蹤航天器運(yùn)行在目標(biāo)航天器附近。故本文相對(duì)運(yùn)動(dòng)特指，以某個(gè)Halo/Lyapunov/Lissajous軌道為參考軌道，研究航天器在此參考軌道附近相對(duì)此參考軌道的運(yùn)動(dòng)。

圖2 相對(duì)運(yùn)動(dòng)示意圖

本文研究相對(duì)運(yùn)動(dòng)的目的是獲得與近地軌道CW方程類似的CRTBP下相對(duì)運(yùn)動(dòng)的近似解析解。下面以Halo軌道作為目標(biāo)航天器的運(yùn)行軌道為例，介紹相對(duì)運(yùn)動(dòng)的研究方法。

在CRTBP中，單值矩陣(Monodromy Matrix)M是指一個(gè)軌道周期T對(duì)應(yīng)的周期軌道的狀態(tài)轉(zhuǎn)移矩陣(State transition matrix,STM)Φ(t0+t,x0)，即

M=Φ(t0+T,x0)

(20)

雖然單值矩陣沒(méi)有精確的解析表達(dá)式，但可以通過(guò)數(shù)值積分的方式獲得精確的數(shù)值解，詳細(xì)的計(jì)算方法可以參考文獻(xiàn)[25]。

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

(21)

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

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

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

綜上所述，單值矩陣的六個(gè)特征值和對(duì)應(yīng)的特征向量提供了一種具有幾何意義的表示，即不同的特征值和對(duì)應(yīng)的特征向量可以代表相對(duì)運(yùn)動(dòng)中具有不同特性的分量模態(tài)，如不穩(wěn)定模態(tài)、穩(wěn)定模態(tài)、中性模態(tài)以及周期模態(tài)，這也是Floquet理論[34]的核心思想。因此，通過(guò)對(duì)單值矩陣的分析，可以建立由一組獨(dú)立基向量所形成的一個(gè)局部坐標(biāo)系統(tǒng)，并以此局部坐標(biāo)系統(tǒng)來(lái)描述相對(duì)運(yùn)動(dòng)。借鑒Floquet理論的思想，這組基向量可以由單值矩陣的特征向量形成，而這組基向量也就是本文提出的相對(duì)運(yùn)動(dòng)模態(tài)。

假設(shè)相對(duì)運(yùn)動(dòng)模態(tài)矩陣E(t)由6個(gè)相對(duì)運(yùn)動(dòng)模態(tài)組成

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

(22)

式中：ei(t),i=1,2,3,…,6即為6個(gè)相對(duì)運(yùn)動(dòng)模態(tài)，分別與單值矩陣的6個(gè)特征值和特征向量相對(duì)應(yīng)。

通過(guò)對(duì)單值矩陣的分析，可以令Halo軌道上當(dāng)前狀態(tài)對(duì)應(yīng)的單值矩陣的特征向量組成相對(duì)運(yùn)動(dòng)模態(tài)矩陣E(t)的初始值E0，即

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

(23)

式中：vi(i=1,…,6)即為單值矩陣的6個(gè)特征向量，可以通過(guò)求解單值矩陣獲得，real(v5)表示特征向量的實(shí)部，imag(v6)表示特征向量的虛部。不難看出rank(E0)=6，即這組基向量是線性無(wú)關(guān)的，可以用來(lái)描述6維空間中的任意向量，即可以用來(lái)描述任意初始狀態(tài)的相對(duì)運(yùn)動(dòng)。

根據(jù)Halo軌道的性質(zhì)[30]以及狀態(tài)轉(zhuǎn)移矩陣性質(zhì)，相對(duì)運(yùn)動(dòng)模態(tài)矩陣滿足如下關(guān)系

E(t)=Φ(t)E0

(24)

在獲得相對(duì)運(yùn)動(dòng)模態(tài)變化規(guī)律E(t)的基礎(chǔ)上，便可以將相對(duì)運(yùn)動(dòng)看做上述6個(gè)相對(duì)運(yùn)動(dòng)模態(tài)的線性組合。假設(shè)追蹤航天器和目標(biāo)航天器的運(yùn)動(dòng)狀態(tài)分別記為x′和x，則兩者之間的差值即為航天器間的相對(duì)運(yùn)動(dòng)δx，

(25)

將相對(duì)運(yùn)動(dòng)δx表示成相對(duì)運(yùn)動(dòng)模態(tài)的線性組合形式，即

(26)

式中：C(t)為各個(gè)模態(tài)在相對(duì)運(yùn)動(dòng)中所占的權(quán)值矩陣，不難發(fā)現(xiàn)其為一個(gè)對(duì)角矩陣，即

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

c5(t),c6(t))

(27)

而相對(duì)運(yùn)動(dòng)的初始狀態(tài)δx0可以表示為

δx0=C0E0

(28)

根據(jù)文獻(xiàn)[35]，由于相對(duì)運(yùn)動(dòng)的尺度相對(duì)參考軌道尺度來(lái)說(shuō)為一小量，因此可以將追蹤航天器相對(duì)于目標(biāo)航天器的相對(duì)運(yùn)動(dòng)δx看做目標(biāo)航天器狀態(tài)x的一種攝動(dòng)，則根據(jù)狀態(tài)轉(zhuǎn)移矩陣的性質(zhì)，有

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

(29)

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

(30)

從式(30)可以得到

C(t)=C0

(31)

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

這樣，便獲得了Halo軌道為參考軌道的相對(duì)運(yùn)動(dòng)近似解析解

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

(32)

由于任意時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣Φ(t)可以事先獲得，即作為已知量，依據(jù)此近似解析解，便可以獲得任意時(shí)刻的相對(duì)運(yùn)動(dòng)狀態(tài)。

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

根據(jù)上述分析，本文提出的CRTBP共線平動(dòng)點(diǎn)附近相對(duì)運(yùn)動(dòng)近似解析解的適用條件如下：

1)目標(biāo)航天器運(yùn)行的參考軌道為Halo軌道、Lyapunov軌道時(shí)的相對(duì)運(yùn)動(dòng)，以及精度要求不高的情況下，也適用于Lissajous軌道為參考軌道的相對(duì)運(yùn)動(dòng)。

2)相對(duì)運(yùn)動(dòng)相比參考軌道的尺度較小，即相對(duì)運(yùn)動(dòng)可以作為參考軌道的攝動(dòng)來(lái)處理。

對(duì)比參考文獻(xiàn)[11]，本文提出的相對(duì)運(yùn)動(dòng)研究方法不僅適用于Halo軌道，也適用于Lyapunov軌道和Lissajous軌道，是一種在Floquet理論基礎(chǔ)上改進(jìn)的、通用的、基于單值矩陣的CRTBP共線平動(dòng)點(diǎn)附近周期/擬周期軌道下相對(duì)運(yùn)動(dòng)的解析研究方法。

3 仿真校驗(yàn)

為校驗(yàn)本文所提出的基于單值矩陣的CRTBP共線平動(dòng)點(diǎn)附近周期/擬周期軌道下的相對(duì)運(yùn)動(dòng)解析研究方法的正確性，本節(jié)分別以地-月系統(tǒng)L1點(diǎn)附近的Halo軌道、Lissajous軌道和Lyapunov軌道作為目標(biāo)航天器的參考軌道，研究追蹤航天器與目標(biāo)航天器之間的相對(duì)運(yùn)動(dòng)。

3.1 Halo軌道

假設(shè)目標(biāo)航天器和追蹤航天器均運(yùn)行在各自的Halo軌道上，相應(yīng)的Z軸振幅分別記Az1和Az2，相對(duì)運(yùn)動(dòng)時(shí)間記為tTOF，xT0為目標(biāo)航天器的初始狀態(tài)，xC0為追蹤航天器的初始狀態(tài)，以上各參數(shù)在L1會(huì)合坐標(biāo)系下的無(wú)量綱單位(Dimensionless unit, DU)取值如表1所示。

Halo軌道相對(duì)運(yùn)動(dòng)的6個(gè)模態(tài)的初值以及一個(gè)軌道周期內(nèi)的變化趨勢(shì)分別如表2和圖3所示。從仿真結(jié)果可以看出，模態(tài)1呈發(fā)散趨勢(shì)，證明了其為不穩(wěn)定模態(tài)；模態(tài)2呈收斂趨勢(shì)，證明了其為穩(wěn)定模態(tài)；模態(tài)3和模態(tài)4呈既不發(fā)散也不收斂趨勢(shì)，證明了其為中性模態(tài)；模態(tài)5和模態(tài)6呈周期變化趨勢(shì)，證明了其為周期模態(tài)。

表1 Halo軌道相對(duì)運(yùn)動(dòng)相關(guān)參數(shù)

表2 Halo軌道相對(duì)運(yùn)動(dòng)6個(gè)模態(tài)初始值

圖3 Halo軌道相對(duì)運(yùn)動(dòng)6個(gè)模態(tài)位置空間的變化規(guī)律

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

圖4 Halo軌道相對(duì)運(yùn)動(dòng)仿真結(jié)果

3.2 Lissajous軌道

假設(shè)目標(biāo)航天器和追蹤航天器均運(yùn)行在各自的Lissajous軌道上，相應(yīng)的軌道參數(shù)[A1,A2,A3,A4,A5,A6]分別記A1和A2，相對(duì)運(yùn)動(dòng)時(shí)間記為tTOF，xT0為目標(biāo)航天器的初始狀態(tài)，xC0為追蹤航天器的初始狀態(tài)，以上各參數(shù)在L1會(huì)合坐標(biāo)系下的無(wú)量綱單位取值如表3所示。

表3 Lissajous軌道相對(duì)運(yùn)動(dòng)相關(guān)參數(shù)

Lissajous軌道相對(duì)運(yùn)動(dòng)的6個(gè)模態(tài)的初值以及一個(gè)周期內(nèi)的變化趨勢(shì)分別如表4和圖6所示。從仿真結(jié)果可以看出，6個(gè)相對(duì)運(yùn)動(dòng)模態(tài)同樣均與其自身的性質(zhì)相符合。

圖5表示真實(shí)Lissajous軌道相對(duì)運(yùn)動(dòng)與反演Lissajous軌道相對(duì)運(yùn)動(dòng)的對(duì)比。從圖5可以看出，反演相對(duì)運(yùn)動(dòng)與真實(shí)相對(duì)運(yùn)動(dòng)相一致，誤差很小，說(shuō)明了本方法的有效性。

圖5 Lissajous軌道相對(duì)運(yùn)動(dòng)仿真結(jié)果

表4 Lissajous軌道相對(duì)運(yùn)動(dòng)6個(gè)模態(tài)初始值

3.3 Lyapunov軌道

假設(shè)目標(biāo)航天器和追蹤航天器均運(yùn)行在各自的Lyapunov軌道上，相應(yīng)的軌道參數(shù)[A1,A2,A3,A4,A5,A6]分別記A1和A2，相對(duì)運(yùn)動(dòng)時(shí)間記為tTOF，xT0為目標(biāo)航天器的初始狀態(tài)，xC0為追蹤航天器的初始狀態(tài)，以上各參數(shù)在L1會(huì)合坐標(biāo)系下的無(wú)量綱單位取值如表5所示。

Lyapunov軌道相對(duì)運(yùn)動(dòng)的6個(gè)模態(tài)的初值以及一個(gè)周期內(nèi)的變化趨勢(shì)分別如表6和圖7所示。從仿真結(jié)果可以看出，6個(gè)相對(duì)運(yùn)動(dòng)模態(tài)同樣均與其自身的性質(zhì)相符合。

圖8表示真實(shí)Lyapunov軌道相對(duì)運(yùn)動(dòng)與反演Lyapunov軌道相對(duì)運(yùn)動(dòng)的對(duì)比。從圖8可以看出，反演相對(duì)運(yùn)動(dòng)與真實(shí)相對(duì)運(yùn)動(dòng)相一致，誤差很小，說(shuō)明了本方法的有效性。

表5 Lyapunov軌道相對(duì)運(yùn)動(dòng)相關(guān)參數(shù)

表6 Lyapunov軌道相對(duì)運(yùn)動(dòng)6個(gè)模態(tài)初始值

圖7 Lyapunov軌道相對(duì)運(yùn)動(dòng)6個(gè)模態(tài)位置子空間的變化規(guī)律

圖8 Lyapunov軌道相對(duì)運(yùn)動(dòng)仿真結(jié)果

4 結(jié) 論

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