地月L1點(diǎn)低能轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化

摘 要：針對(duì)地月空間平動(dòng)點(diǎn)周期軌道與近地軌道之間的低能轉(zhuǎn)移問(wèn)題，提出一種地月L1（Earth-Moon L EML1）點(diǎn)Halo軌道到地球靜止軌道（geostationary Earth orbit， GEO）的四脈沖低能轉(zhuǎn)移軌道的設(shè)計(jì)方法。所提方法在擾動(dòng)流形和Lambert弧段拼接的三脈沖轉(zhuǎn)移軌道設(shè)計(jì)基礎(chǔ)上，從分析軌道雅可比常數(shù)變化與速度增量關(guān)系的角度出發(fā)設(shè)計(jì)四脈沖低能轉(zhuǎn)移軌道。數(shù)值仿真結(jié)果表明，四脈沖優(yōu)化模型比三脈沖模型效率更高，可以得到更優(yōu)的轉(zhuǎn)移方案，有效解決了優(yōu)化過(guò)程中由于搜索空間大、極值數(shù)量多而導(dǎo)致的優(yōu)化結(jié)果不佳的問(wèn)題。所提設(shè)計(jì)方法可以用于EML1其他周期軌道族與各類近地軌道的相互轉(zhuǎn)移問(wèn)題研究。

關(guān)鍵詞： 地月空間; 平動(dòng)點(diǎn); 周期軌道; 不變流形; 低能轉(zhuǎn)移軌道

中圖分類號(hào)： V 412.4+1 文獻(xiàn)標(biāo)志碼： A""" DOI：10.12305/j.issn.1001-506X.2024.10.28

Design and optimization of Earth-Moon L1 low-energy transfer orbit

QIAO Chenyuan， YANG Leping^*

（College of Aerospace Science and Engineering， National University of Defense Technology，

Changsha 410073， China）

Abstract： Aiming at the low-energy transfer problem between the periodic orbit of the cislunar space libration point and near-Earth orbit， a design method of 4-impulses low-energy transfer orbit from Earth-Moon L1 （EML1） point Halo orbit to geostationary Earth orbit （GEO） is proposed. Based on the design of 3-impulses transfer orbit which applies disturbed manifold and Lambert arc， the 4-impulses low-energy transfer orbit is designed by analyzing the relationship between the change of Jacobi constant and the velocity increment. Numerical simulation results show that the 4-impulse optimization model is more efficient than 3-impulse model and a better transfer scheme can be obtained， which effectively solves the problem of obtaining poor optimization results due to large search space and numerous extreme values in the optimization process. The proposed design method can be used to study the transfer between other periodic orbit families of EML1 and various near-Earth orbits as well.

Keywords： cislunar space; libration point; periodic orbit; invariant manifold; low-energy transfer orbit

0 引 言

時(shí)至今日，人類的太空活動(dòng)大多數(shù)集中于近地空間，受制于技術(shù)和成本要求，地球同步軌道至月球軌道的空間以及月球軌道外的空間則相對(duì)平靜。地球靜止軌道（geostationary Earth orbit， GEO）和地月平動(dòng)點(diǎn)因其特殊的空間位置備受關(guān)注。對(duì)于搶占地月空間資源、構(gòu)建地月空間有利態(tài)勢(shì)而言，研究GEO軌道與地月平動(dòng)點(diǎn)軌道之間的機(jī)動(dòng)轉(zhuǎn)移具有重要價(jià)值。一方面，利用地月平動(dòng)點(diǎn)軌道到GEO軌道的高效轉(zhuǎn)移，未來(lái)可以用發(fā)射地月平動(dòng)點(diǎn)載荷的運(yùn)載器“拼車”搭載GEO載荷，實(shí)現(xiàn)運(yùn)載能力最大化，并提高GEO載荷發(fā)射靈活性、隱蔽性;另一方面，利用地月空間動(dòng)力學(xué)特性，可以將平動(dòng)點(diǎn)作為近地空間到地月空間的轉(zhuǎn)移中樞，實(shí)現(xiàn)近地軌道到地月空間軌道、月球軌道的低能往返，未來(lái)可以基于此構(gòu)建地月空間中的“轉(zhuǎn)移運(yùn)輸網(wǎng)絡(luò)”。

針對(duì)平動(dòng)點(diǎn)轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題，一般的轉(zhuǎn)移軌道設(shè)計(jì)思路為直接轉(zhuǎn)移^［13］。Parker^［2］等深入研究地月系統(tǒng)下的直接轉(zhuǎn)移軌道，指出直接轉(zhuǎn)移的時(shí)間約在3天到2個(gè)月，并計(jì)算各次速度增量的范圍和總共耗費(fèi)能量的范圍，其中在Halo軌道上施加的速度脈沖一般在500～600 m/s。隨著對(duì)平動(dòng)點(diǎn)軌道認(rèn)知的深入，不變流形逐漸被應(yīng)用到低能量轉(zhuǎn)移軌道的設(shè)計(jì)中^［4］。研究發(fā)現(xiàn)，利用不變流形設(shè)計(jì)的轉(zhuǎn)移軌道所需的速度增量小于直接轉(zhuǎn)移方式，在Halo軌道流形上運(yùn)動(dòng)的航天器能夠以很小的速度增量進(jìn)入或離開(kāi)Halo軌道，無(wú)需在Halo軌道上施加較大的速度脈沖。利用不變流形設(shè)計(jì)的轉(zhuǎn)移軌道主要針對(duì)日地L1（Sun-Earth L SEL1）點(diǎn)周期軌道與近地軌道之間的轉(zhuǎn)移^［57］、地月平動(dòng)點(diǎn)與近月軌道之間的轉(zhuǎn)移^［89］和不同平動(dòng)點(diǎn)周期軌道之間的轉(zhuǎn)移^［1011］。對(duì)于地月平動(dòng)點(diǎn)與近地軌道的轉(zhuǎn)移問(wèn)題，連一君^［12］采用EML1周期軌道不穩(wěn)定流形上轉(zhuǎn)移點(diǎn)與近地軌道構(gòu)建Lambert轉(zhuǎn)移的方式設(shè)計(jì)二脈沖低能轉(zhuǎn)移軌道;張漢清等^［13］利用擾動(dòng)流形研究EML1到地球的低能轉(zhuǎn)移問(wèn)題，對(duì)象為地球近地軌道與Lyapunov軌道的二脈沖共面轉(zhuǎn)移;Rosales等^［14］在雙圓模型下研究地球到EML2-Halo軌道的轉(zhuǎn)移;Neelakantan等^［15］在橢圓限制性三體問(wèn)題下借助不變流形設(shè)計(jì)二脈沖轉(zhuǎn)移方案;美國(guó)國(guó)家航天局（National Aeronautics and Space Administration， NASA）^［16］利用不變流形設(shè)計(jì)多種從近地軌道到EML1-Halo軌道的三脈沖轉(zhuǎn)移軌道，并分別給出轉(zhuǎn)移方案和控制策略。近年來(lái)的研究方向逐漸轉(zhuǎn)向?qū)⑿⊥屏εc不變流形相結(jié)合實(shí)現(xiàn)近地軌道到平動(dòng)點(diǎn)的轉(zhuǎn)移^［1720］，以及轉(zhuǎn)向利用不變流形實(shí)現(xiàn)平動(dòng)點(diǎn)周期軌道到遠(yuǎn)距離逆行軌道（distant retrograde orbit， DRO）^［2123］和近直線暈軌道（near-rectilinear halo orbit， NRHO）^［2426］的轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題。

以上研究主要針對(duì)三體系統(tǒng)中平動(dòng)點(diǎn)周期軌道與小質(zhì)量天體之間的轉(zhuǎn)移以及平動(dòng)點(diǎn)附近軌道之間的轉(zhuǎn)移，如地球到SEL1點(diǎn)，地月L1（Earth-Moon L EML1）點(diǎn)到月球以及EML1和EML2、DRO、NRHO之間的轉(zhuǎn)移^［27］，因?yàn)檫@些平動(dòng)點(diǎn)周期軌道的不變流形可以與小質(zhì)量天體或者其他流形相交，進(jìn)而可以通過(guò)拼接流形來(lái)實(shí)現(xiàn)低能轉(zhuǎn)移軌道設(shè)計(jì)。而對(duì)于大質(zhì)量天體的轉(zhuǎn)移，例如EML1到地球，由于流形并不會(huì)到達(dá)地球附近，這類轉(zhuǎn)移軌道設(shè)計(jì)相對(duì)繁瑣，相關(guān)的研究也較少。文獻(xiàn)［12］的方法簡(jiǎn)便，能夠快速求解出優(yōu)化結(jié)果，但是沒(méi)有充分利用地月空間的動(dòng)力學(xué)特征，使得求解結(jié)果還有進(jìn)一步優(yōu)化的空間。文獻(xiàn)［13］的方法僅適用于平面轉(zhuǎn)移，且優(yōu)化參數(shù)多，對(duì)優(yōu)化算法提出了較高的要求，并且該模型要求龐加萊截面圖相交，這一條件屬于過(guò)約束條件，進(jìn)一步阻礙了模型求解。

因此，針對(duì)EML1點(diǎn)與近地軌道的轉(zhuǎn)移問(wèn)題，本文將文獻(xiàn)［13］擾動(dòng)流形的設(shè)計(jì)思路由二維擴(kuò)展到三維，采用Lambert弧段拼接的方法將優(yōu)化模型的過(guò)約束條件進(jìn)行松弛，建立由EML1點(diǎn)到GEO的三脈沖低能轉(zhuǎn)移軌道設(shè)計(jì)模型。針對(duì)優(yōu)化過(guò)程極值數(shù)目多、搜索空間大的特點(diǎn)，從分析軌道雅可比常數(shù)變化與速度增量關(guān)系的角度出發(fā)，提出了四脈沖低能轉(zhuǎn)移軌道設(shè)計(jì)方法。最后，通過(guò)數(shù)值仿真比較了不同模型的優(yōu)化結(jié)果和計(jì)算效率。

1 軌道動(dòng)力學(xué)描述

1.1 圓形限制性三體問(wèn)題

針對(duì)地月空間轉(zhuǎn)移軌道設(shè)計(jì)問(wèn)題，基于圓形限制性三體問(wèn)題（circular restricted three-body problem， CRTBP）進(jìn)行分析。在地月空間的CRTBP中，通常在質(zhì)心旋轉(zhuǎn)坐標(biāo)系O-XYZ中描述航天器的運(yùn)動(dòng)，如圖1所示。

其中，O-XYZ為慣性系，原點(diǎn)O為地月系統(tǒng)質(zhì)心，z軸與慣性系下的Z軸指向一致，x軸與地月連線重合并由第一主天體（地球）指向第二主天體（月球）方向，y軸根據(jù)右手系確定。

在質(zhì)心旋轉(zhuǎn)系中，航天器的動(dòng)力學(xué)方程可以表示為

d²xdτ²－2dydτ=Ωx

d²ydτ²+2dxdτ=Ωy

d²zdτ²=Ωz（1）

式中：

Ω=12（x²+y²）+1－μr₁+μr₂

r₁=（x+μ）²+y²+z²

r₂=（x－1+μ）²+y²+z²（2）

為了簡(jiǎn)化動(dòng)力學(xué)方程，并在數(shù)值積分中獲得更好的數(shù)值精度，對(duì)式（1）中的各個(gè)物理量全部進(jìn)行了無(wú)量綱化處理，表1羅列了一些典型的單位換算關(guān)系，下文若無(wú)特殊說(shuō)明均使用無(wú)量綱單位。

式（1）中存在一個(gè)積分常數(shù)C：

C=2Ω－（x·²+y·²+z·²）（3）

稱作雅可比積分，由于其具有能量量綱，又被稱為雅可比能量。雅可比積分與航天器機(jī)械能的關(guān)系為C=－2E。

1.2 平動(dòng)點(diǎn)周期軌道與狀態(tài)轉(zhuǎn)移矩陣

航天器在旋轉(zhuǎn)坐標(biāo)系下的運(yùn)動(dòng)狀態(tài)可以表示為X=［x，y，z，x·，y·，z·］T，式（1）可以寫(xiě)為

X·（t）=f（X（t），t）（4）

EML1點(diǎn)是位于地球和月球之間的共線平動(dòng)點(diǎn)，在平動(dòng)點(diǎn)上的航天器將與兩個(gè)主天體保持相對(duì)靜止的位置關(guān)系，其附近存在的周期軌道能夠?yàn)殚L(zhǎng)期的觀測(cè)以及航天器各類對(duì)月任務(wù)部署提供有利的條件。周期軌道和擬周期軌道根據(jù)其空間特性可以分為不同的軌道族^［2829］，本文針對(duì)與月球軌道共面的Lyapunov軌道與三維的Halo軌道進(jìn)行研究，如圖2所示。

狀態(tài)轉(zhuǎn)移矩陣是動(dòng)力學(xué)中狀態(tài)方程對(duì)于初始狀態(tài)的導(dǎo)數(shù)。在CRTBP中，狀態(tài)轉(zhuǎn)移矩陣反映了參考軌道對(duì)于初始狀態(tài)微小擾動(dòng)δX（t₀）的線性化特征，即

δX（t）=Φ（t，t₀）δX（t₀）（5）

狀態(tài)轉(zhuǎn)移矩陣滿足

Φ·（t，t₀）=M（X）Φ（t，t₀）

Φ（t₀，t₀）=I₆（6）

式中：I₆為6階單位矩陣;M（X）為雅可比矩陣，滿足

M（X）=O_3×3I₃

Ω_XXC（7）

式中：C=020

－200

000; Ω_XX=Ω_xxΩ_xyΩ_xz

Ω_yxΩ_yyΩ_yz

Ω_zxΩ_zyΩ_zz為Ω關(guān)于r=［x，y，z］T的二階偏導(dǎo)數(shù)。根據(jù)式（6）可以計(jì)算出任意時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣。

1.3 不變流形與龐加萊截面

不變流形是與周期軌道光滑連接的一簇軌道構(gòu)成的流管，分為遠(yuǎn)離周期軌道的不穩(wěn)定流形和進(jìn)入周期軌道的穩(wěn)定流形，如圖3所示。航天器在流形上的演化不需要耗費(fèi)能量。不變流形的計(jì)算需要利用單值矩陣的特性。單值矩陣Φ（T，0）反映了周期軌道中航天器初始的偏差經(jīng)過(guò)一個(gè)周期的映射關(guān)系，其特征值滿足

λ₁gt;1

λ₂=1/λ₁

λ₃=λ₄=1

λ₅=λ-₆， |λ₅|=1（8）

λ₁的模大于1，說(shuō)明其對(duì)應(yīng)特征向量方向的偏差會(huì)逐漸放大，則該特征向量指向了不穩(wěn)定的方向，記為Vu;相應(yīng)地，λ₂的模小于1，說(shuō)明其對(duì)應(yīng)特征向量方向的偏差會(huì)逐漸減小，對(duì)應(yīng)的方向指向穩(wěn)定方向，把該特征向量記為Vs。假設(shè)周期軌道上航天器的狀態(tài)為X₀，則不穩(wěn)定流形和穩(wěn)定流形的初值Xu（X₀）、Xs（X₀）分別為

Xu（X₀）=X₀±ε（Vu/Vu）

Xs（X₀）=X₀±ε（Vs/Vs）

（9）

式中：·為L(zhǎng)₂范數(shù);ε為小量;±表示不變流形向周期軌道兩側(cè)的兩個(gè)分支。對(duì)于不穩(wěn)定流形而言，流形代表的是質(zhì)點(diǎn)離開(kāi)周期軌道后的運(yùn)動(dòng)狀態(tài)，需要對(duì)Xu（X₀）進(jìn)行正向積分;對(duì)于穩(wěn)定流形而言，流形代表的是質(zhì)點(diǎn)趨向周期軌道之前的運(yùn)動(dòng)狀態(tài)，需要對(duì)Xs（X₀）進(jìn)行逆向積分。

對(duì)于CRTBP中不變流形的研究，希望在保留其空間結(jié)構(gòu)的基礎(chǔ)上盡可能降低相空間的維數(shù)，龐加萊截面是實(shí)現(xiàn)上述目的非常有效的方式。傳統(tǒng)的龐加萊截面法選取x軸平面或者y軸平面為龐加萊截面，選?。▁，x·）或者（y，y·）作為狀態(tài)點(diǎn)在龐加萊截面上的坐標(biāo)，但這種方法主要針對(duì)周期軌道的研究，這里采用一種以角度和距離作為截面坐標(biāo)的方法，稱為θ-r截面法^［30］。

首先，針對(duì)平面CRTBP來(lái)分析，即忽略z軸方向的運(yùn)動(dòng)。設(shè)龐加萊截面與z軸正方向夾角為α，截面上的交點(diǎn)狀態(tài)為［x，y，v_x，v_y］T，令θ為速度方向與截面法向的夾角，r為交點(diǎn)距離地球的距離，如圖4所示，則

r=（x－μ）²+y²

θ=α+π2－arctanv_yv_x（10）

當(dāng)問(wèn)題上升至三維空間，龐加萊截面的結(jié)構(gòu)也會(huì)隨之變得更加復(fù)雜，具體而言就是θ-r截面法中的每一個(gè)參量都要上升一維。這里取龐加萊截面為從地球出發(fā)、與xOy平面垂直的截面，即平面CRTBP下的龐加萊截面增加z軸方向的自由度。由于空間維度上升，龐加萊截面上的截面坐標(biāo)也要擴(kuò)充到4個(gè)，即［r，θ，β，γ］，其中各分量的含義如圖5所示。令龐加萊截面P的法向量為n_P，地球指向航天器的矢量為r，由n_P和r確定的平面為P₁，P₁的法向量為n，由n和n_P確定的平面為P₂，則θ代表速度矢量在P₁的投影與n_P間的夾角，β代表速度矢量在P₂的投影與n_P的夾角，γ代表r與xOy平面的夾角。

龐加萊截面的核心思想是通過(guò)降低維數(shù)來(lái)簡(jiǎn)化不變流形的研究。上述龐加萊截面無(wú)論是在平面CRTBP下或是在三維空間下都省去了一個(gè)位置的維度和一個(gè)速度的維度，這種方法可以在后續(xù)的軌道設(shè)計(jì)過(guò)程中忽略無(wú)關(guān)因素的干擾，利用龐加萊截面圖迅速找到合適的軌道拼接點(diǎn)。圖6展示了三維空間中Halo軌道不穩(wěn)定流形的龐加萊截面。

在旋轉(zhuǎn)坐標(biāo)系下，GEO軌道會(huì)隨著月球公轉(zhuǎn)而不斷向西移動(dòng)，最終形成類似手鐲的形狀。同樣，可以用龐加萊截面的方法來(lái)對(duì)GEO軌道進(jìn)行研究：圖7展示了GEO軌道在0°龐加萊截面上的截面圖，可以看到GEO軌道的r坐標(biāo)和θ坐標(biāo)變?yōu)榱艘粋€(gè)固定值，β坐標(biāo)和γ坐標(biāo)則在白赤交角范圍之內(nèi)隨軌道西進(jìn)不斷往復(fù)。因?yàn)閞坐標(biāo)始終等于GEO軌道的軌道半徑，θ坐標(biāo)表示速度方向與周向方向的夾角，GEO軌道是圓軌道，速度方向始終與周向方向重合，因此夾角保持0°不變。根據(jù)對(duì)稱性，GEO軌道在任意角度的龐加萊截面內(nèi)都有如圖7所示的截面圖。

2 地月低能轉(zhuǎn)移軌道設(shè)計(jì)

2.1 三脈沖低能轉(zhuǎn)移軌道設(shè)計(jì)

利用平動(dòng)點(diǎn)實(shí)現(xiàn)的空間飛行任務(wù)中，大多數(shù)是在日地系統(tǒng)中開(kāi)展的，因?yàn)镾EL1的流形能夠直接到達(dá)地球附近，航天器可以直接進(jìn)入流形，使得轉(zhuǎn)移軌道的設(shè)計(jì)相對(duì)簡(jiǎn)單。但是在地月系統(tǒng)中，EML1點(diǎn)的流形無(wú)法直接到達(dá)地球附近，因此需要采用其他方法來(lái)設(shè)計(jì)軌道。參考文獻(xiàn)［13］，這里采用一種擾動(dòng)流形的方式實(shí)現(xiàn)軌道拼接，目標(biāo)軌道是白道面內(nèi)的共面地球同步軌道（geosynchronous orbit， GSO）。EML1點(diǎn)的不穩(wěn)定流形無(wú)法到達(dá)共面GSO的高度，如果設(shè)法在流形上的某處進(jìn)行軌道機(jī)動(dòng)，原不穩(wěn)定流形就會(huì)在機(jī)動(dòng)的作用下發(fā)生變形，稱為擾動(dòng)流形。通過(guò)控制機(jī)動(dòng)的大小和方向，就可以控制擾動(dòng)流形與GSO相交，進(jìn)而設(shè)計(jì)出Lyapunov軌道到共面GSO的轉(zhuǎn)移軌道。

如圖8所示，可以看到該擾動(dòng)流形剛好與GSO相切，其相切點(diǎn)（也就是龐加萊截面圖中的交點(diǎn)）就是流形拼接的點(diǎn)，由此也就可以確定轉(zhuǎn)移軌道在不穩(wěn)定流形中的位置，進(jìn)一步可確定轉(zhuǎn)移軌道在Lyapunov周期軌道起點(diǎn)的位置。選擇相切點(diǎn)作為拼接點(diǎn)的原因是希望避免航天器在調(diào)整速度方向上額外耗費(fèi)速度增量。同理，在三維的轉(zhuǎn)移軌道設(shè)計(jì)中，也希望找到擾動(dòng)流形與GEO龐加萊截面的交點(diǎn)作為軌道的拼接點(diǎn)，但是地月空間的動(dòng)力學(xué)特性復(fù)雜，GEO的龐加萊截面與擾動(dòng)流形的截面是四維空間中的閉合曲線，令其相交非常困難，這導(dǎo)致龐加萊截面相交的約束條件過(guò)于嚴(yán)格，難以展開(kāi)優(yōu)化計(jì)算。這里通過(guò)Lambert轉(zhuǎn)移實(shí)現(xiàn)擾動(dòng)流形末端與GEO的拼接。

假設(shè)擾動(dòng)流形在截面Σ上的坐標(biāo)集合為P，GEO的坐標(biāo)集合為Q，P和Q上的點(diǎn)為

X_P=［r_P，θ_P，β_P，γ_P］， X_P∈P

X_Q=［r_Q，θ_Q，β_Q，γ_Q］， X_Q∈Q（11）

選擇P和Q在龐加萊截面空間中距離最短的點(diǎn)來(lái)確定Lambert弧段連接的起點(diǎn)和終點(diǎn)，即

minxX_P－X_Q

s.t.

x=［X_P，X_Q］T

X_P∈P

X_Q∈Q

（12）

這里選擇距離最短的點(diǎn)是因?yàn)辇嫾尤R截面中的4個(gè)坐標(biāo)代表了位置的偏差和速度方向的偏差，這樣選擇能夠盡可能減少位置和速度調(diào)整所需耗費(fèi)的額外的速度增量。令轉(zhuǎn)移時(shí)間為T_l，根據(jù)X_P，X_Q計(jì)算出對(duì)應(yīng)的旋轉(zhuǎn)坐標(biāo)系坐標(biāo)，再將坐標(biāo)轉(zhuǎn)換至地心慣性系，即可轉(zhuǎn)化為L(zhǎng)ambert問(wèn)題并進(jìn)行計(jì)算。

這里需要強(qiáng)調(diào)的是，盡管Lambert問(wèn)題是二體運(yùn)動(dòng)下的邊值問(wèn)題，但是該問(wèn)題仍可以適用于CRTBP下的邊值問(wèn)題求解。通過(guò)前期選取龐加萊截面上最小距離的點(diǎn)，待計(jì)算的Lambert弧段的起點(diǎn)與終點(diǎn)都位于GEO附近，該弧段可以看作是月球引力攝動(dòng)下的二體Lambert轉(zhuǎn)移，使用二體Lambert轉(zhuǎn)移計(jì)算出的軌道與實(shí)際的拼接弧段相差不大。為了使得計(jì)算的弧段更精確，可以先按照二體Lambert轉(zhuǎn)移計(jì)算弧段的初值，再根據(jù)CRTBP模型設(shè)計(jì)微分校正法對(duì)弧段起點(diǎn)處施加的機(jī)動(dòng)進(jìn)行校正，以獲得更符合三體運(yùn)動(dòng)規(guī)律的Lambert拼接弧段。將上述計(jì)算Lambert弧段的計(jì)算過(guò)程表示為

［Δv₂，Δv₃］=Lambert（X_P，X_Q，T_l）（13）

Lambert弧段各參量含義如圖9所示。這種方法并不需要擾動(dòng)流形的截面與GEO截面相交，相當(dāng)于將優(yōu)化模型的強(qiáng)約束條件轉(zhuǎn)換為選取X_P和X_Q，并通過(guò)式（13）轉(zhuǎn)換為L(zhǎng)ambert轉(zhuǎn)移的速度增量，添加到優(yōu)化函數(shù)中，這種松弛處理可以解決限制條件過(guò)于嚴(yán)格導(dǎo)致難以開(kāi)展優(yōu)化的問(wèn)題，對(duì)優(yōu)化算法運(yùn)行的效率有非常大的提升作用。

下面建立三維空間中地月轉(zhuǎn)移軌道優(yōu)化模型。假設(shè)航天器在雅可比能量為C=3.138 4的Halo軌道上運(yùn)動(dòng)，目標(biāo)是進(jìn)入GEO。對(duì)于給定的ε值，可以計(jì)算出該周期軌道的不穩(wěn)定流形W_u。在進(jìn)入不穩(wěn)定流形Δt時(shí)間后，施加大小為Δv₁的速度增量，方向?yàn)樗俣确较蚶@z軸順時(shí)針旋轉(zhuǎn)角度α_dis，航天器進(jìn)入擾動(dòng)流形W_dis，運(yùn)行至角度為α的龐加萊截面時(shí)停止積分，積分終止條件為

crit（X）=arctanyx+μ－α（14）

令W_dis在龐加萊截面上形成的截面圖為集合P，GEO軌道的龐加萊截面為Q。選取P和Q上距離最近的點(diǎn)X_P和X_Q，對(duì)于給定的時(shí)間T_l計(jì)算出Lambert轉(zhuǎn)移軌道，由此可以計(jì)算出Lambert弧段始末的速度增量Δv₂和Δv₃，之后航天器即可切入GEO。設(shè)全過(guò)程的飛行時(shí)間為τ，則需要被優(yōu)化的參數(shù)為

x=［ε，Δt，Δv₁，α_dis，α，T_l］T（15）

能量最優(yōu)的地月轉(zhuǎn)移軌道優(yōu)化模型為

minx f₁（x）=Δv₁+Δv₂+Δv₃

s.t. x=［ε，Δt，Δv₁，α_dis，α，T_l］T（16）

能量時(shí)間最優(yōu)的地月轉(zhuǎn)移軌道優(yōu)化模型為

minx f₂（x）=Δv₁+Δv₂+Δv₃+kτ

s.t. x=［ε，Δt，Δv₁，α_dis，α，T_l］T（17）

2.2 四脈沖低能轉(zhuǎn)移軌道設(shè)計(jì)

上述擾動(dòng)流形拼接的方法在優(yōu)化過(guò)程中基本是在盲目搜索，對(duì)于如何施加速度增量、在什么位置施加速度增量等問(wèn)題模型并未給出明確的結(jié)果，這也就產(chǎn)生了搜索空間龐大、極值點(diǎn)數(shù)目多、優(yōu)化算法容易陷入極小值等問(wèn)題。下面提出一種通過(guò)分析軌道雅可比常數(shù)變化與速度增量之間的解析關(guān)系來(lái)設(shè)計(jì)四脈沖低能轉(zhuǎn)移軌道的方法。

選擇一條按照三脈沖設(shè)計(jì)方法得到的低能轉(zhuǎn)移軌道，計(jì)算其雅可比常數(shù)，結(jié)果如圖10所示?？梢钥吹?，Halo軌道的C值較小，能量較高，而GEO的C值較大，能量較低。從Halo軌道到GEO的轉(zhuǎn)移實(shí)質(zhì)上需要提高航天器的C值，降低能量。而設(shè)計(jì)低能轉(zhuǎn)移軌道，就是希望能夠以最少的速度增量Δv來(lái)將能量降低到對(duì)應(yīng)水平。下面從軌道的能量變化入手來(lái)進(jìn)行分析。

根據(jù)式（3）可知雅可比常數(shù)C的計(jì)算公式為

C=2Ω－v²（18）

當(dāng)航天器位置確定時(shí)，對(duì)式（18）兩邊同時(shí)計(jì)算變分，有

δC=－2vδv（19）

整理后可得

δv=－δC2v（20）

式（20）提供了很多信息：首先，v代表速度的模，有v≥0，因此δC與δv的符號(hào)相反，這就說(shuō)明如果希望提高C值，就需要對(duì)航天器進(jìn)行減速機(jī)動(dòng);其次，對(duì)于特定大小的δv，可以看到v越大，則δC越大，說(shuō)明在速度較大的地方施加機(jī)動(dòng)，相同的δv可以使航天器的C值提高更多一些;最后，δC的大小取決于速度大小的變化δv，這就說(shuō)明速度增量Δv要盡可能使得速度大小改變，而不是使得速度方向改變，因此進(jìn)行減速機(jī)動(dòng)時(shí)Δv的方向要盡可能與速度方向共線。同理，如果需要調(diào)整速度方向，則應(yīng)該在速度較小的地方施加機(jī)動(dòng)，盡可能減小不使速度大小發(fā)生變化的Δv分量。

基于上述分析提出低能轉(zhuǎn)移的機(jī)動(dòng)法則：一是用于航天器減速的機(jī)動(dòng)要在速度較大處施加，且速度增量方向要與速度方向共線;二是用于航天器速度調(diào)整的機(jī)動(dòng)要在速度較小處施加。基于以上法則，可以得到如下優(yōu)化模型。

假設(shè)航天器在雅可比能量為C=3.138 4的Halo軌道上運(yùn)動(dòng)，目標(biāo)是進(jìn)入GEO。對(duì)于給定的ε值，可以計(jì)算出該周期軌道的不穩(wěn)定流形W_u。在進(jìn)入不穩(wěn)定流形后，在最小速度處施加機(jī)動(dòng)Δv₁=［Δv_1x，Δv_1y，Δv_1z］T，調(diào)整速度方向進(jìn)入擾動(dòng)流形?？梢酝ㄟ^(guò)下面方法判別是否到達(dá)最小速度：令v和a為速度和加速度矢量，判別函數(shù)為crit（X）=a·v，當(dāng)判別函數(shù)正向過(guò)零時(shí)，說(shuō)明此刻速度達(dá)到極小值。同理，當(dāng)判別函數(shù)負(fù)向過(guò)零時(shí)，速度達(dá)到極大值。在調(diào)整速度后，在到達(dá)速度最大處進(jìn)行減速機(jī)動(dòng)Δv₂，方向與速度方向相反。運(yùn)動(dòng)T時(shí)間后機(jī)動(dòng)進(jìn)入Lambert轉(zhuǎn)移弧段。因此，

x=［ε，Δv_1x，Δv_1y，Δv_1z，Δv₂，T，T_l］T（21）

可以看到，雖然待優(yōu)化的參量上升到7個(gè)，但是該設(shè)計(jì)方案從三脈沖機(jī)動(dòng)改變?yōu)樗拿}沖機(jī)動(dòng)，相較之前可以探索更多可能的低能轉(zhuǎn)移軌道。這種方法的缺陷是設(shè)計(jì)模型時(shí)完全從節(jié)省能量的角度出發(fā)，軌道的結(jié)構(gòu)被高度約束，因此采用能量時(shí)間最優(yōu)目標(biāo)的意義并不大，這里只用這種方法設(shè)計(jì)優(yōu)化目標(biāo)為能量最優(yōu)的模型，優(yōu)化模型為

minx f₃（x）=Δv₁+Δv₂+Δv₃+Δv₄

s.t. x=［ε，Δv_1x，Δv_1y，Δv_1z，Δv₂，T，T_l］T

Δv₁=［Δv_1x，Δv_1y，Δv_1z］T

［Δv₃，Δv₄］=Lambert（X_P，X_Q，T_l）（22）

3 地月低能轉(zhuǎn)移軌道設(shè)計(jì)

對(duì)于上述優(yōu)化模型，采用復(fù)合粒子群算法進(jìn)行求解^［31］。該算法在前期強(qiáng)調(diào)局部最優(yōu)粒子的作用，以提高算法的全局搜索能力;后期則強(qiáng)調(diào)全局最優(yōu)粒子的作用，提高算法的收斂性和精度。

3.1 三脈沖轉(zhuǎn)移優(yōu)化結(jié)果

3.1.1 能量最優(yōu)

在能量最優(yōu)的模型下，航天器在進(jìn)入不變流形后一段時(shí)間內(nèi)施加機(jī)動(dòng)進(jìn)入擾動(dòng)流形，到達(dá)合適的窗口后機(jī)動(dòng)進(jìn)行Lambert轉(zhuǎn)移并進(jìn)入GEO。在能量最優(yōu)的目標(biāo)下，航天器同樣要在地月空間附近多次盤旋以尋找最佳的轉(zhuǎn)移窗口，因此需要花費(fèi)較長(zhǎng)的時(shí)間。能量最優(yōu)優(yōu)化目標(biāo)下計(jì)算的轉(zhuǎn)移軌道如圖11所示，各次機(jī)動(dòng)的速度增量及時(shí)間如表2所示。

3.1.2 能量時(shí)間最優(yōu)

在能量時(shí)間最優(yōu)的模型下，航天器在進(jìn)入擾動(dòng)流形階段需要額外花費(fèi)更多能量，但是轉(zhuǎn)移的總時(shí)間下降到能量最優(yōu)的50%以下。同時(shí)，通過(guò)與直接轉(zhuǎn)移軌道對(duì)比，可以發(fā)現(xiàn)能量時(shí)間最優(yōu)的軌跡路徑基本相同，在之后的優(yōu)化中可以根據(jù)先驗(yàn)信息縮小對(duì)參數(shù)的搜索范圍以獲得能量更優(yōu)或者時(shí)間更短的解。能量時(shí)間最優(yōu)優(yōu)化目標(biāo)下計(jì)算的轉(zhuǎn)移軌道如圖12所示，各次機(jī)動(dòng)的速度增量及時(shí)間如表3所示。

3.2 四脈沖轉(zhuǎn)移優(yōu)化結(jié)果

結(jié)果表明，四脈沖轉(zhuǎn)移模型得到的最優(yōu)解在能量方面優(yōu)于單純使用擾動(dòng)流形拼接法得到的解，最終得到總速度增量只需1.469 5 km/s的轉(zhuǎn)移方案，有效解決了低能轉(zhuǎn)移軌道的優(yōu)化問(wèn)題。在處理復(fù)雜的軌道轉(zhuǎn)移問(wèn)題時(shí)，不妨從機(jī)理分析出發(fā)，提前推測(cè)最優(yōu)解可能出現(xiàn)的位置，并據(jù)此重新設(shè)計(jì)優(yōu)化模型，這樣可以有效簡(jiǎn)化或壓縮搜索空間，更快得到更優(yōu)的解。

在四脈沖優(yōu)化模型中，根據(jù)低能轉(zhuǎn)移的機(jī)動(dòng)法則，可以在最大速度機(jī)動(dòng)Δv₂后再次選擇速度極大/極小位置進(jìn)行機(jī)動(dòng)，構(gòu)建五脈沖、六脈沖或者更多脈沖機(jī)動(dòng)的轉(zhuǎn)移模型，但是這里不提倡采納更多脈沖的轉(zhuǎn)移方案。原因如下：第一，通過(guò)計(jì)算圖13中各個(gè)位置的速度大小發(fā)現(xiàn)，在Δv₂機(jī)動(dòng)后，后續(xù)不存在比Δv₁速度更小的位置，調(diào)整速度方向的機(jī)動(dòng)代價(jià)較大;第二，等待后續(xù)速度極大點(diǎn)增加脈沖機(jī)動(dòng)次數(shù)可能會(huì)進(jìn)一步節(jié)省燃料，但是需要更長(zhǎng)的轉(zhuǎn)移時(shí)間以及更大的軌道控制成本，節(jié)省的速度增量難以彌補(bǔ)更長(zhǎng)轉(zhuǎn)移時(shí)間帶來(lái)的軌道修正所需的速度增量。因此，這里暫時(shí)不考慮更多脈沖的優(yōu)化模型。

表4所示為四脈沖能量最優(yōu)模型機(jī)動(dòng)參數(shù)的統(tǒng)計(jì)情況。

四脈沖模型與三脈沖模型的本質(zhì)區(qū)別在于四脈沖模型通過(guò)式（20）將各次機(jī)動(dòng)的位置確定為最大/最小位置處，符合軌道力學(xué)原理，仿真證明了該方法的有效性。通過(guò)跟蹤兩次優(yōu)化過(guò)程中每代的最優(yōu)速度增量（見(jiàn)圖14）可以看出，四脈沖轉(zhuǎn)移模型在算法迭代到一半時(shí)已經(jīng)達(dá)到三脈沖模型計(jì)算出的最優(yōu)值，通過(guò)后續(xù)的優(yōu)化可以計(jì)算出更優(yōu)的結(jié)果，更適于算法進(jìn)行尋優(yōu)。從計(jì)算結(jié)果來(lái)看，四脈沖優(yōu)化模型所需的時(shí)間也要優(yōu)于三脈沖模型，因?yàn)樗拿}沖模型在遇到速度極大/極小點(diǎn)處會(huì)提前終止數(shù)值積分，相較三脈沖模型計(jì)算量更小，計(jì)算效率更高。兩種模型的優(yōu)化參數(shù)對(duì)比如表5所示。

4 結(jié)束語(yǔ)

本文針對(duì)EML1點(diǎn)Halo軌道到GEO的低能轉(zhuǎn)移需求，在擾動(dòng)流形和Lambert弧段拼接的三脈沖轉(zhuǎn)移軌道設(shè)計(jì)基礎(chǔ)上，從分析軌道雅可比常數(shù)變化與速度增量關(guān)系的角度出發(fā)設(shè)計(jì)四脈沖低能轉(zhuǎn)移軌道。仿真結(jié)果表明，該方法計(jì)算效率更高，可以得到更優(yōu)的轉(zhuǎn)移方案，有效解決了優(yōu)化過(guò)程中由于搜索空間大、極值數(shù)量多而導(dǎo)致優(yōu)化結(jié)果不佳的問(wèn)題，可以用于地月低能轉(zhuǎn)移軌道的高效計(jì)算。最終得到總速度增量只需1.469 5 km/s的低能轉(zhuǎn)移方案，有效解決了低能轉(zhuǎn)移軌道的設(shè)計(jì)與優(yōu)化問(wèn)題。通過(guò)修改GEO參數(shù)，上述模型還可以擴(kuò)展至設(shè)計(jì)往返Halo軌道和不同軌道半徑和傾角的近地軌道的轉(zhuǎn)移方案，對(duì)于地月平動(dòng)點(diǎn)轉(zhuǎn)移軌道的設(shè)計(jì)與應(yīng)用、相關(guān)參數(shù)選擇等具有重要的借鑒意義。

參考文獻(xiàn)

［1］ RAUSCH R R. Earth to Halo orbit transfer trajectories［D］. West Lafayette： Indiana： Purdue University， 2005.

［2］ PARKER J S， ANDERSON R L. Low-energy lunar trajectory design［M］. Hoboken： Wiley， 2014.

［3］ TAN M H， ZHANG K， WANG J Y. Optimization of bi-impulsive Earth-Moon transfers using periodic orbits［J］. Astrophysics and Space Science， 202 366： 19.

［4］ ZHANG R Y. A review of periodic orbits in the circular restricted three-body problem［J］. Journal of Systems Engineering and Electronics， 202 33（3）： 612646.

［5］ GMEZ G， KOON W S， LO M W， et al. Connecting orbits and invariant manifolds in the spatial restricted three-body problem［J］. Nonlinearity， 2004， 17（5）： 15711606.

［6］ BARDEN B T. Using stable manifolds to generate transfers in the circular restricted problem of three bodies［D］. West Lafayette： Purdue University， 1994.

［7］ HOWELL K C， MAINS D L， BARDEN B T. Transfer trajectories from Earth parking orbits to Sun-Earth Halo orbits［C］∥Proc.of the AAS/AIAA Astrodynamics Specialist Conference， 1994： 94160.

［8］ JIN Y， XU B. Three-maneuver transfers from the cislunar L2 Halo orbits to low lunar orbits［J］. Advances in Space Research， 202 69（2）： 989999.

［9］ SANTOS L B T， SOUSA-SILVA P A， TERRA M O， et al. Optimal transfers from Moon to L2 Halo orbit of the Earth-Moon system［J］. Advances in Space Research， 202 70（11）： 33623372.

［10］ 徐明， 徐世杰. 地月系平動(dòng)點(diǎn)及Halo軌道的應(yīng)用研究［J］. 宇航學(xué)報(bào)， 2006， 27（4）： 695699.

XU M， XU S J. The application of libration points and Halo orbits in the earth-moon system to space mission design［J］. Journal of Astronautics， 2006， 27（4）： 695699.

［11］ 于錫崢， 鄭建華， 高懷寶， 等. 地月系L1和L2點(diǎn)間轉(zhuǎn)移軌道設(shè)計(jì)［J］. 吉林大學(xué)學(xué)報(bào)（工學(xué)版）， 2008， 38（3）： 741745.

YU X Z， ZHENG J H， GAO H B， et al. Transfer trajectory design between L1 and L2 in the Earth-Moon system［J］. Journal of Jilin University （Engineering and Technology Edition）， 2008， 38（3）： 741745.

［12］ 連一君. 基于三體平動(dòng)點(diǎn)的低能轉(zhuǎn)移軌道設(shè)計(jì)研究［D］. 長(zhǎng)沙： 國(guó)防科學(xué)技術(shù)大學(xué)， 2008.

LIAN Y J. Research on low-cost transfer trajectory design based on three-body libration points［D］. Changsha： National University of Defense Technology， 2008.

［13］ 張漢清， 李言俊. 地月三體問(wèn)題下L1-地球低能轉(zhuǎn)移軌道設(shè)計(jì)［J］. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)， 201 43（5）： 8488.

ZHANG H Q， LI Y J. Design of L1-Earth low energy transfer trajectory in Earth-Moon three-body problem［J］. Journal of Harbin Institute of Technology， 201 43（5）： 8488.

［14］ ROSALES J J， JORBA ， JORBA-CUSC M. Transfers from the Earth to L2 Halo orbits in the Earth-Moon bicircular problem［J］. Celestial Mechanics and Dynamical Astronomy， 202 133： 55.

［15］ NEELAKANTAN R， RAMANAN R V. Two-impulse transfer to multi-revolution Halo orbits in the Earth-Moon elliptic restricted three body problem framework［J］. Journal of Astrophysics and Astronomy， 202 43（2）： 5066.

［16］ FOLTA D， BOSANAC N， ELLIOTT I， et al. Astrodynamics convention and modeling reference for Lunar， Cislunar， and libration point orbits［R］. Hampton： NASA Langley Research Center， 2022： 127145.

［17］ SINGH S K， ANDERSON B D， TAHERI E， et al. Exploiting manifolds of L1 Halo orbits for end-to-end Earth-Moon low-thrust trajectory design［J］. Acta Astronautica， 202 183（1）： 255272.

［18］ OSHIMA K. Optimization-blemaided， low-energy transfers via vertical instability around Earth-Moon L1［J］. Journal of Guidance， Control， and Dynamics， 202 44（2）： 389398.

［19］ DU C， STARINOVA O L， LIU Y. Transfer between the planar Lyapunov orbits around the Earth-Moon L2 point using low-thrust engine［J］. Acta Astronautica， 202 201（1）： 513525.

［20］ DU C， STARINOVA O L， LIU Y. Low-thrust transfer dynamics and control between Halo orbits in the Earth-Moon system by means of invariant manifold［J］. IEEE Trans.on Aerospace and Electronic Systems， 2023， 59（4）： 34523462.

［21］ WANG Y， ZHANG R K， ZHANG C， et al. Transfers between NRHOs and DROs in the Earth-Moon system［J］. Acta Astronautica， 202 186（1）： 6073.

［22］ CAPDEVILA L， GUZZETTI D， HOWELL K. Various transfer options from Earth into distant retrograde orbits in the vicinity of the Moon［J］. Advances in the Astronautical Sciences， 2014， 152： 36593678.

［23］ CAPDEVILA L R， HOWELL K C. A transfer network linking Earth， Moon， and the triangular libration point regions in the Earth-Moon system［J］. Advances in Space Research， 2018， 62（7）： 18261852.

［24］ ZIMOVAN E M， HOWELL K C， DAVIS D C. Near rectilinear Halo orbits and their application in cislunar space［C］∥Proc.of the 3rd AIAA Conference on Dynamics and Control of Space Systems， 2017： 2040.

［25］ BOUDAD K K， HOWELL K C， DAVIS D C. Dynamics of synodic resonant near rectilinear Halo orbits in the bicircular four-body problem［J］. Advances in Space Research， 2020， 66（9）： 21942214.

［26］ TROFIMOV S， SHIROBOKOV M， TSELOUSOVA A， et al. Transfers from near-rectilinear Halo orbits to low-perilune orbits and the Moon’s surface［J］. Acta Astronautica， 2020， 167（1）： 260271.

［27］ 李翔宇， 喬棟， 程潏. 三體軌道動(dòng)力學(xué)研究進(jìn)展［J］. 力學(xué)學(xué)報(bào)， 202 53（5）： 12231245.

LI X Y， QIAO D， CHENG Y. Progress of three-body orbital dynamics study［J］. Chinese Journal of Theoretical and Applied Mechanics， 202 53（5）： 12231245.

［28］ CHOW C C， WETTERER C J， HILL K， et al. Cislunar periodic orbit families and expected observational features［C］∥Proc.of the Advanced Maui Optical and Space Surveillance Technologies Conference， 2021.

［29］ 錢霙婧， 荊武興， 劉玥， 等. 地月平動(dòng)點(diǎn)擬周期軌道設(shè)計(jì)方法［J］. 系統(tǒng)工程與電子技術(shù)， 2014， 36（8）： 15861594.

QIAN Y J， JING W X， LIU Y， et al. Design of quasi periodic orbit about the translunar libration point［J］. Systems Engineering and Electronics， 2014， 36（8）： 15861594.

［30］ 李言俊， 張科， 呂梅柏， 等. 利用拉格朗日點(diǎn)的深空探測(cè)技術(shù)［M］. 西安： 西北工業(yè)大學(xué)出版社， 2015： 4449.

LI Y J， ZHANG K， LYU M B， et al. Deep space exploration techniques using Lagrange points［M］. Xi’an： Northwestern Polytechnical University Press， 2015： 4449.

［31］ WU L H， WANG Y N， YUAN X F， et al. Research and application of compound optimum model particle swarm optimization［J］. Journal of Systems Engineering and Electronics， 2006， 28（7）： 10871090.

作者簡(jiǎn)介

喬琛遠(yuǎn)（1999—），男，碩士研究生，主要研究方向?yàn)榈卦驴臻g軌道動(dòng)力學(xué)與控制。

楊樂(lè)平（1964—），男，教授，博士，主要研究方向?yàn)楹教烊蝿?wù)規(guī)劃、空間電磁操控。