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

(1.北京航天飛行控制中心，北京100094；2.西安衛(wèi)星測控中心宇航動力學國家重點實驗室，西安710043；3.中國航天員科研訓練中心，北京100094)

1 引言

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

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

2 動力學模型

2.1 圓型限制性三體問題

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

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

在歸一化單位下，第三體P在會合坐標系中的動力學方程[4]如式(2)：

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

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

2.2 Halo軌道及不變流形

Halo軌道是共線平動點附近的三維周期軌道，在平動點任務中被廣泛應用。CR3BP會合坐標系下，Halo軌道關于xz面對稱，與xz面交于兩點，通常取其中距離x軸較遠的點與x軸的距離作為表征Halo軌道大小的參數(shù)，稱之為法向幅值Az，計算過程可參考文獻[4]。

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

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

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

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

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

本文針對LEO出發(fā)達到LL2點空間站的貨物補給運輸任務，綜合不變流形與月球借力兩種節(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)移軌道段和不變流形軌道段，二者在月球附近實現(xiàn)拼接；整個轉(zhuǎn)移過程共需要2次脈沖機動，第一次機動的位置為近地LEO，作用是使航天器逃離LEO進入地月轉(zhuǎn)移軌道，第二次機動的位置為月球附近，作用是使航天器插入穩(wěn)定流形。

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

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

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

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

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

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

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

3.2 求解策略

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

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

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

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

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

圖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個離散點分布圖Fig.3 Distribution map of the 360 discrete points on Halo orbit

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

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

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

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

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

式中，uE為地球引力常數(shù)；r1為LEO軌道半長軸；r2為φ＝90°插入點的地心距。約束條件見式(9)，當約束條件滿足時，通過活力公式即可解析計算出ΔVLEO。 優(yōu)化目標為總?cè)剂舷淖钌?，目標函?shù)見式(10)。優(yōu)化得到的φ＝90°插入點對應的最優(yōu)插入脈沖記為對應的入軌點切向逃逸脈沖大小記為。

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

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

3.3 算例驗證

LEO-LL2 Halo間接轉(zhuǎn)移軌道設計實例參數(shù)及約束配置如下：LEO高度hLEO＝185 km；目標流形軌道近月點高度hS＝200 km；地月轉(zhuǎn)移段飛行時間TES≤ 7.4 d；Halo軌道幅值Az＝10000 km，為LL2點南向Halo軌道。

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

圖6 hS＝200 km的目標流形軌道對應的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)移軌道特性分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

圖11 不同幅值Halo軌道對應的目標流形的三維圖Fig.11 The 3D maps of target manifold orbits corresponding to the Halo orbits with different Az