地月L2點(diǎn)空間站轉(zhuǎn)移軌道設(shè)計(jì)與特性分析

曹鵬飛，賀波勇，劉景勇，彭祺擘

(1.北京航天飛行控制中心，北京100094；2.西安衛(wèi)星測(cè)控中心宇航動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安710043；3.中國(guó)航天員科研訓(xùn)練中心，北京100094)

1 引言

地月L2點(diǎn)(以下簡(jiǎn)稱LL2點(diǎn))位于月球背面約6.5×104km處，在月球探測(cè)中具有重要作用。早在Apollo時(shí)代，F(xiàn)arquhar就發(fā)現(xiàn)LL2點(diǎn)附近存在大量三維周期軌道，即Halo軌道[1]，并首次提出了Halo軌道空間站的概念[2-3]，指出與環(huán)月軌道空間站相比，LL2點(diǎn)Halo軌道空間站在月球探測(cè)方面更有優(yōu)勢(shì)，包括與月面站有更長(zhǎng)的通信時(shí)長(zhǎng)、與地球通信無(wú)遮擋、兼顧月球背面通信和運(yùn)營(yíng)方式相對(duì)自主等[4]。但由于資金等諸多原因，該方案未被實(shí)施。1978年，美國(guó)成功發(fā)射了世界首顆平動(dòng)點(diǎn)探測(cè)器ISEE-3[5]，在國(guó)際上掀起了一輪平動(dòng)點(diǎn)轉(zhuǎn)移軌道研究熱潮。Gomez等[6]通過(guò)引入動(dòng)力系統(tǒng)理論，發(fā)現(xiàn)共線平動(dòng)點(diǎn)附近存在大量無(wú)動(dòng)力飛行的不變流形。Parker等[7]基于不變流形，通過(guò)大量仿真數(shù)據(jù)研究了從近地停泊軌道(Low Earth Orbit，LEO)出發(fā)到達(dá)LL2點(diǎn)Halo軌道的轉(zhuǎn)移軌道特性，發(fā)現(xiàn)在月球附近插入不變流形，可同時(shí)節(jié)省一定的燃料和時(shí)間成本。Gordon[8]基于月球借力與不變流形，通過(guò)采用改進(jìn)的微分修正方法設(shè)計(jì)了LEO出發(fā)到達(dá)LL2點(diǎn)Halo軌道的轉(zhuǎn)移軌道。Li等[9]在Gordon的研究基礎(chǔ)上，提出了基于月球借力的LL2點(diǎn)Halo軌道三脈沖轉(zhuǎn)移策略，進(jìn)一步降低了燃料消耗。張景瑞等[10]基于遺傳算法與微分修正算法結(jié)合的混合優(yōu)化技術(shù)，研究了不同月球借力約束下的LL2點(diǎn)Halo軌道轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題。

上述學(xué)者對(duì)基于不變流形與月球借力的LL2點(diǎn)間接轉(zhuǎn)移軌道設(shè)計(jì)方法做了比較深入的研究，但對(duì)該類轉(zhuǎn)移軌道的一般特性規(guī)律研究尚不充分。本文以未來(lái)LL2點(diǎn)Halo軌道空間站的貨物補(bǔ)給任務(wù)為背景，結(jié)合解析初值猜想和局部梯度優(yōu)化，提出一種串行間接轉(zhuǎn)移軌道優(yōu)化方法，并在此基礎(chǔ)上，通過(guò)大量仿真算例研究不變流形插入相位、月球借力高度以及Halo軌道幅值等參數(shù)對(duì)轉(zhuǎn)移軌道燃料消耗的影響特性。以期為未來(lái)LL2點(diǎn)Halo軌道低能轉(zhuǎn)移任務(wù)軌道方案制定提供參考。

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

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

本文所選用的動(dòng)力學(xué)模型為圓型限制性三體問(wèn)題(Circular Restricted Three Body Problem，CR3BP)模型，即在一個(gè)三體系統(tǒng)中，假設(shè)質(zhì)量分別為M1和M2的兩大天體P1和P2，繞其公共質(zhì)心做勻速圓周運(yùn)動(dòng)，研究質(zhì)量為m的第三體P在P1和P2引力下的運(yùn)動(dòng)規(guī)律。三個(gè)天體的質(zhì)量關(guān)系為M1＞M2?m，即P對(duì)P1和P2的運(yùn)動(dòng)影響可忽略不計(jì)[4]。

CR3BP中常用的坐標(biāo)系為質(zhì)心會(huì)合坐標(biāo)系[4]，其原點(diǎn)位于P1和P2的公共質(zhì)心O，x軸由P1指向P2，z軸指向系統(tǒng)角動(dòng)量矢量方向，y軸與其他兩軸構(gòu)成右手坐標(biāo)系。引入歸一化單位，長(zhǎng)度單位為P1和P2質(zhì)心間的距離，質(zhì)量單位為M1、M2之和，P1和P2的相對(duì)運(yùn)動(dòng)周期為2π。 引入質(zhì)量比μ如式(1)：

在歸一化單位下，第三體P在會(huì)合坐標(biāo)系中的動(dòng)力學(xué)方程[4]如式(2)：

其中，U為與位置相關(guān)的偽勢(shì)能函數(shù)，表達(dá)式為式(3)：

式中，r1和r2分別為P到P1和P2的距離，表達(dá)式如式(4)、(5)：

2.2 Halo軌道及不變流形

Halo軌道是共線平動(dòng)點(diǎn)附近的三維周期軌道，在平動(dòng)點(diǎn)任務(wù)中被廣泛應(yīng)用。CR3BP會(huì)合坐標(biāo)系下，Halo軌道關(guān)于xz面對(duì)稱，與xz面交于兩點(diǎn)，通常取其中距離x軸較遠(yuǎn)的點(diǎn)與x軸的距離作為表征Halo軌道大小的參數(shù)，稱之為法向幅值A(chǔ)z，計(jì)算過(guò)程可參考文獻(xiàn)[4]。

不變流形是與Halo軌道自然相銜接的空間軌道，分為穩(wěn)定流形與不穩(wěn)定流形，在低成本深空軌道設(shè)計(jì)中被廣泛應(yīng)用。將Halo軌道離散為不動(dòng)點(diǎn)，并求解不動(dòng)點(diǎn)處的單值矩陣Φ(T，0)，即一個(gè)Halo軌道周期后的狀態(tài)轉(zhuǎn)移矩陣。每一不動(dòng)點(diǎn)對(duì)應(yīng)的Φ(T，0)矩陣均存在三對(duì)乘積為1的特征值：λu＞1為不穩(wěn)定特征值，對(duì)應(yīng)不穩(wěn)定特征向量Vu；λs＜1為穩(wěn)定特征值，對(duì)應(yīng)穩(wěn)定特征向量Vs；一對(duì)互為共軛的復(fù)數(shù)特征值和一對(duì)值為1的特征值。在不動(dòng)點(diǎn)處施加一個(gè)小擾動(dòng)即可得到不變流形積分初值。

對(duì)于不穩(wěn)定流形，其積分初值Xu0如式(6)：

對(duì)于穩(wěn)定流形，其積分初值Xs0如式(7)：

式中：X0為不動(dòng)點(diǎn)狀態(tài)量；d為擾動(dòng)變量。對(duì)于地月系統(tǒng)，為了滿足方程線性，d一般取30～70 km。

2.3 LEO-LL2 Halo間接轉(zhuǎn)移軌道方案

本文針對(duì)LEO出發(fā)達(dá)到LL2點(diǎn)空間站的貨物補(bǔ)給運(yùn)輸任務(wù)，綜合不變流形與月球借力兩種節(jié)省燃料方式，提出LEO-LL2 Halo間接轉(zhuǎn)移軌道方案，如圖1所示。

圖1 LEO-LL2 Halo間接轉(zhuǎn)移軌道示意圖Fig.1 Schematic diagram of LEO-LL2 Halo indirect transfer trajectory

由圖1可知：LEO-LL2 Halo間接轉(zhuǎn)移軌道共包括2段，分別是地月轉(zhuǎn)移軌道段和不變流形軌道段，二者在月球附近實(shí)現(xiàn)拼接；整個(gè)轉(zhuǎn)移過(guò)程共需要2次脈沖機(jī)動(dòng)，第一次機(jī)動(dòng)的位置為近地LEO，作用是使航天器逃離LEO進(jìn)入地月轉(zhuǎn)移軌道，第二次機(jī)動(dòng)的位置為月球附近，作用是使航天器插入穩(wěn)定流形。

3 LEO-LL2 Halo間接轉(zhuǎn)移軌道設(shè)計(jì)

本節(jié)將CR3BP下LEO-LL2 Halo間接轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題，轉(zhuǎn)換為有約束非線性規(guī)劃求解問(wèn)題。建立優(yōu)化模型，并在此基礎(chǔ)上，提出一種解析初值猜想與局部梯度優(yōu)化相結(jié)合的求解策略。

3.1 軌道優(yōu)化模型

下面從設(shè)計(jì)變量、約束條件和目標(biāo)函數(shù)三方面介紹其優(yōu)化模型。

設(shè)計(jì)變量為逃逸點(diǎn)高度hE、逃逸點(diǎn)切向逃逸脈沖ΔVLEO、穩(wěn)定流形(可由近月點(diǎn)高度hS確定)、地月轉(zhuǎn)移段轉(zhuǎn)移時(shí)間TES、穩(wěn)定流形插入位置XIS，穩(wěn)定流形插入脈沖ΔVIS以及Halo軌道幅值A(chǔ)z。設(shè)計(jì)變量均已知時(shí)，在CR3BP下可唯一確定一條LEO-LL2 Halo間接轉(zhuǎn)移軌道。其中，hE＝hLEO，由發(fā)射系統(tǒng)確定；hS根據(jù)測(cè)控系統(tǒng)精度確定；Az由任務(wù)目標(biāo)給定。因此，一條LEO-LL2 Halo間接轉(zhuǎn)移軌道可由式(8)所示的4元素確定：

約束條件主要考慮轉(zhuǎn)移時(shí)間和近地點(diǎn)高度約束。參照嫦娥四號(hào)中繼星任務(wù)[11]，本文將工程與軌道約束簡(jiǎn)化為式(9)：

由于空間站貨物補(bǔ)給任務(wù)對(duì)轉(zhuǎn)移時(shí)間的要求并不高，本節(jié)在軌道設(shè)計(jì)時(shí)弱化了時(shí)間約束，以總?cè)剂舷臑閮?yōu)化目標(biāo)，即目標(biāo)函數(shù)為式(10)：

3.2 求解策略

針對(duì)上述優(yōu)化模型求解問(wèn)題，提出一種解析初值猜想與局部梯度優(yōu)化相結(jié)合的求解策略，具體流程如圖2所示。

1)根據(jù)火箭運(yùn)載能力與任務(wù)目標(biāo)，確定LEO軌道高度hLEO、穩(wěn)定流形近月點(diǎn)高度hS和目標(biāo)Halo軌道。

2)計(jì)算目標(biāo)流形軌道。首先，將Halo軌道按時(shí)間等間距離散成360份，并按時(shí)間順序?qū)﹄x散點(diǎn)進(jìn)行編號(hào)n＝1～360，如圖3所示；其次，計(jì)算360個(gè)離散點(diǎn)對(duì)應(yīng)的穩(wěn)定流形，并選出近月點(diǎn)高度與hS一致的穩(wěn)定流形，作為目標(biāo)流形軌道。

圖4給出了滿足不同近月點(diǎn)高度約束的穩(wěn)定流形，由圖可知，當(dāng)近月點(diǎn)高度給定時(shí)，將存在兩條穩(wěn)定流形與之對(duì)應(yīng)，分別稱之為“A型”和“B型”穩(wěn)定流形。其中，B型穩(wěn)定流形在繞月后徑直奔向Halo軌道，因此不適用于從LEO出發(fā)到達(dá)Halo軌道的轉(zhuǎn)移任務(wù)。本文選用A型穩(wěn)定流形軌道作為目標(biāo)流形軌道。

定義目標(biāo)流形軌道上位于月球附近的積分點(diǎn)為插入點(diǎn)，并將其相對(duì)月球的相位記為φ。如圖5所示，φ與插入點(diǎn)位置XIS構(gòu)成一一對(duì)應(yīng)關(guān)系。

圖2 基于CR3BP的LEO-LL2 Halo轉(zhuǎn)移軌道求解流程圖Fig.2 Flow chart of solution of LEO-LL2 Halo transfer trajectory based on CR3BP

圖3 Halo軌道上360個(gè)離散點(diǎn)分布圖Fig.3 Distribution map of the 360 discrete points on Halo orbit

圖4 滿足近月點(diǎn)高度約束的穩(wěn)定流形Fig.4 Stable manifold orbits satisfying the attitude constraint at perilune

圖5 不變流形插入點(diǎn)相位定義Fig.5 Definition of phase of the insertion point on invariant manifold orbit

3)基于時(shí)間逆向積分求解思路，通過(guò)采用SQP算法優(yōu)化一條從LEO出發(fā)到達(dá)φ＝90°插入點(diǎn)的地月轉(zhuǎn)移軌道。優(yōu)化變量如式(11)所示：

迭代初值如式(12)所示：

其中，ΔVIS0由二體模型霍曼轉(zhuǎn)移理論[12]給出，計(jì)算公式如式(13)：

式中，uE為地球引力常數(shù)；r1為L(zhǎng)EO軌道半長(zhǎng)軸；r2為φ＝90°插入點(diǎn)的地心距。約束條件見(jiàn)式(9)，當(dāng)約束條件滿足時(shí)，通過(guò)活力公式即可解析計(jì)算出ΔVLEO。 優(yōu)化目標(biāo)為總?cè)剂舷淖钌伲繕?biāo)函數(shù)見(jiàn)式(10)。優(yōu)化得到的φ＝90°插入點(diǎn)對(duì)應(yīng)的最優(yōu)插入脈沖記為對(duì)應(yīng)的入軌點(diǎn)切向逃逸脈沖大小記為。

4)通過(guò)采用數(shù)值延拓策略，計(jì)算從LEO出發(fā)到達(dá)其他插入點(diǎn)的地月轉(zhuǎn)移軌道。由于步驟3已經(jīng)得到從LEO出發(fā)到達(dá)φ＝90°插入點(diǎn)的地月轉(zhuǎn)移軌道，在求解LEO出發(fā)到達(dá)φ＝91°或φ＝89°插入點(diǎn)的地月轉(zhuǎn)移軌道時(shí)，將步驟3中優(yōu)化變量初值設(shè)為ΔVφ＝90LEO即可。依次遞推，可計(jì)算出LEO出發(fā)到達(dá)任意插入點(diǎn)的地月轉(zhuǎn)移軌道。

5)分析插入點(diǎn)相位與總速度增量之間的關(guān)系，即可得到燃料最優(yōu)的LEO-LL2 Halo間接轉(zhuǎn)移軌道。

3.3 算例驗(yàn)證

LEO-LL2 Halo間接轉(zhuǎn)移軌道設(shè)計(jì)實(shí)例參數(shù)及約束配置如下：LEO高度hLEO＝185 km；目標(biāo)流形軌道近月點(diǎn)高度hS＝200 km；地月轉(zhuǎn)移段飛行時(shí)間TES≤ 7.4 d；Halo軌道幅值A(chǔ)z＝10000 km，為L(zhǎng)L2點(diǎn)南向Halo軌道。

通過(guò)采用3.2節(jié)的求解策略，計(jì)算滿足上述約束且燃料最優(yōu)的LEO-LL2 Halo間接轉(zhuǎn)移軌道，計(jì)算得到的軌道飛行軌跡如圖6所示，相應(yīng)的軌道參數(shù)見(jiàn)表1。計(jì)算在計(jì)算機(jī)CPU為2.67 GHz的MATLAB環(huán)境下進(jìn)行，單條軌道的計(jì)算時(shí)長(zhǎng)不超過(guò)30 s。因此，利用該方法進(jìn)行LEO-LL2 Halo間接轉(zhuǎn)移軌道設(shè)計(jì)時(shí)，收斂速度很快，適合用于大規(guī)模的軌道特性分析。

圖6 hS＝200 km的目標(biāo)流形軌道對(duì)應(yīng)的LEO-LL2 Halo間接轉(zhuǎn)移軌道Fig.6 LEO-LL2 Halo indirect transfer trajectory corresponding to the target manifold orbit with 200 km perilune attitude

表1 燃料消耗最優(yōu)的LEO-LL2 Halo間接轉(zhuǎn)移軌道參數(shù)Table 1 Parameters of LEO-(LL2)Halo indirect transfer trajectory with optimal fuel consumption

4 LEO-LL2 Halo間接轉(zhuǎn)移軌道特性分析

在任務(wù)設(shè)計(jì)階段，工程人員往往對(duì)單條軌道參數(shù)并不關(guān)心，而是更加注重這一類軌道的一般規(guī)律。本節(jié)在3.3節(jié)算例的基礎(chǔ)上，通過(guò)仿真實(shí)驗(yàn)對(duì)滿足約束的LEO-LL2 Halo間接轉(zhuǎn)移軌道特性進(jìn)行分析。

4.1 穩(wěn)定流形插入點(diǎn)相位影響分析

針對(duì)插入點(diǎn)相位和轉(zhuǎn)移燃料消耗之間的關(guān)系，以3.3節(jié)算例為基礎(chǔ)，目標(biāo)流形軌道插入點(diǎn)相位φ設(shè)為0～180°，其余參數(shù)設(shè)置保持不變。采用3.2節(jié)求解策略，逐條計(jì)算出從LEO出發(fā)達(dá)到φ＝0～180°插入點(diǎn)的地月轉(zhuǎn)移軌道，得到插入點(diǎn)相位φ與總速度增量ΔVtotal的關(guān)系，如圖7所示。由圖可知：①插入點(diǎn)相位φ對(duì)轉(zhuǎn)移速度增量ΔVtotal的影響較大，當(dāng)φ≈120°時(shí)，ΔVtotal取極小值，約為3.4818 km/s。 ②當(dāng)φ＜120°時(shí)，φ越小，ΔVtotal越大。產(chǎn)生上述現(xiàn)象的原因：當(dāng)φ減小時(shí)，插入點(diǎn)與目標(biāo)流形近月點(diǎn)之間的間距將增大，繼而將導(dǎo)致月球借力效果變差，ΔVtotal增加(圖5)。③當(dāng)φ＞120°時(shí)，φ增大，ΔVtotal將略增大。由圖4可知，這一點(diǎn)是由目標(biāo)流形的“走向”引起的。

4.2 穩(wěn)定流形近月點(diǎn)高度影響分析

圖7 插入點(diǎn)相位與轉(zhuǎn)移總速度增量的關(guān)系Fig.7 Relationship between the phase of insertion point and the total velocity increment of transfer consumption

針對(duì)目標(biāo)流形近月點(diǎn)高度和轉(zhuǎn)移所需的最小燃料消耗之間的關(guān)系，以3.3節(jié)算例為基礎(chǔ)，目標(biāo)流形軌道近月點(diǎn)高度hS設(shè)為100～5000 km，其余參數(shù)設(shè)置保持不變。選出近月點(diǎn)高度hS為0～5000 km的A型目標(biāo)穩(wěn)定流形族如圖8所示。采用3.2節(jié)求解策略，逐條計(jì)算出從LEO出發(fā)達(dá)到所有目標(biāo)穩(wěn)定流形的燃料最優(yōu)的地月轉(zhuǎn)移軌道。得到近月點(diǎn)高度與總速度增量的關(guān)系如圖9所示。

由圖9可知：目標(biāo)流形近月點(diǎn)高度hS對(duì)總速度增量ΔVtotal有一定影響，但整體影響不大；當(dāng)hS取0～1300 km時(shí)，ΔVtotal相對(duì)較小且波動(dòng)變化不大；當(dāng)hS＞1300 km時(shí)，hS越大，轉(zhuǎn)移所需的速度增量越多，該結(jié)論與文獻(xiàn)[9]中借力飛行理論相吻合。

4.3 Halo軌道法向幅值影響分析

圖8 滿足近月點(diǎn)高度約束的穩(wěn)定流形Fig.8 Stable manifold orbits satisfying the attitude constraint at perilune

圖9 目標(biāo)流形軌道近月點(diǎn)高度與轉(zhuǎn)移總脈沖關(guān)系Fig.9 Relationship between attitude of perilune of the target manifold orbit and the total velocity increment of transfer consumption

針對(duì)Halo軌道法向幅值和轉(zhuǎn)移所需的最小燃料消耗之間的關(guān)系，仍以3.3節(jié)算例為基礎(chǔ)，Halo軌道幅值A(chǔ)z＝5000、10000、……、30000 km，目標(biāo)流形軌道插入點(diǎn)相位φ＝120°，近月點(diǎn)高度hS＝200 km，其余參數(shù)設(shè)置保持不變。采用3.2節(jié)求解策略，逐條計(jì)算出從LEO出發(fā)到達(dá)不同Halo軌道的燃料最優(yōu)LEO-LL2 Halo間接轉(zhuǎn)移軌道，得到Halo軌道幅值A(chǔ)z與總速度增量ΔVtotal的關(guān)系如圖10所示。

由圖10可知，Halo軌道的幅值A(chǔ)z對(duì)轉(zhuǎn)移所需總速度增量的影響較為顯著，Az越大，越大，反之亦然。該現(xiàn)象可通過(guò)目標(biāo)流形軌道近月點(diǎn)時(shí)刻的速度矢量與白道面的夾角β解釋。不同Az的Halo軌道對(duì)應(yīng)的目標(biāo)流形軌道如圖11所示。

圖10 Halo軌道法向幅值與總脈沖消耗關(guān)系Fig.10 Relationship between the normal amplitude of Halo orbit and the total velocity increment of transfer consumption

由圖11可知：對(duì)于小幅值(如Az＝5000 km)的Halo軌道，β比較小，月球借力效果較為顯著；而對(duì)于較大幅值(如Az＝30 000 km)的Halo軌道，β則比較大，月球借力效果不佳，導(dǎo)致轉(zhuǎn)移所需的燃料消耗較多。因此，LEO-LL2 Halo間接轉(zhuǎn)移軌道方案更適用于從LEO出發(fā)到達(dá)小幅值Halo軌道的轉(zhuǎn)移任務(wù)。

5 結(jié)論

本文研究了不變流形和月球借力相結(jié)合的間接轉(zhuǎn)移軌道設(shè)計(jì)方法及軌道參數(shù)特性。建立的LEO-LL2 Halo間接轉(zhuǎn)移軌道優(yōu)化模型，提出的解析初值搜索和局部梯度優(yōu)化相結(jié)合的串行優(yōu)化求解方法，通過(guò)算例得到了驗(yàn)證。仿真實(shí)驗(yàn)表明，不變流形插入相位和Halo軌道幅值對(duì)間接轉(zhuǎn)移軌道燃料消耗的影響比較大，而月球借力高度的影響則比較小。研究結(jié)論可為未來(lái)LL2點(diǎn)低能轉(zhuǎn)移任務(wù)的軌道方案制定提供參考。

圖11 不同幅值Halo軌道對(duì)應(yīng)的目標(biāo)流形的三維圖Fig.11 The 3D maps of target manifold orbits corresponding to the Halo orbits with different Az