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

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

中圖分類號： V 412.4+1 文獻標志碼： 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ù)集中于近地空間，受制于技術(shù)和成本要求，地球同步軌道至月球軌道的空間以及月球軌道外的空間則相對平靜。地球靜止軌道（geostationary Earth orbit， GEO）和地月平動點因其特殊的空間位置備受關(guān)注。對于搶占地月空間資源、構(gòu)建地月空間有利態(tài)勢而言，研究GEO軌道與地月平動點軌道之間的機動轉(zhuǎn)移具有重要價值。一方面，利用地月平動點軌道到GEO軌道的高效轉(zhuǎn)移，未來可以用發(fā)射地月平動點載荷的運載器“拼車”搭載GEO載荷，實現(xiàn)運載能力最大化，并提高GEO載荷發(fā)射靈活性、隱蔽性;另一方面，利用地月空間動力學特性，可以將平動點作為近地空間到地月空間的轉(zhuǎn)移中樞，實現(xiàn)近地軌道到地月空間軌道、月球軌道的低能往返，未來可以基于此構(gòu)建地月空間中的“轉(zhuǎn)移運輸網(wǎng)絡(luò)”。

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

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

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

1 軌道動力學描述

1.1 圓形限制性三體問題

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

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

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

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）

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

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

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

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

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

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

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

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

狀態(tài)轉(zhuǎn)移矩陣是動力學中狀態(tài)方程對于初始狀態(tài)的導數(shù)。在CRTBP中，狀態(tài)轉(zhuǎn)移矩陣反映了參考軌道對于初始狀態(tài)微小擾動δ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的二階偏導數(shù)。根據(jù)式（6）可以計算出任意時刻的狀態(tài)轉(zhuǎn)移矩陣。

1.3 不變流形與龐加萊截面

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

λ₁gt;1

λ₂=1/λ₁

λ₃=λ₄=1

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

λ₁的模大于1，說明其對應(yīng)特征向量方向的偏差會逐漸放大，則該特征向量指向了不穩(wěn)定的方向，記為Vu;相應(yīng)地，λ₂的模小于1，說明其對應(yīng)特征向量方向的偏差會逐漸減小，對應(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₂范數(shù);ε為小量;±表示不變流形向周期軌道兩側(cè)的兩個分支。對于不穩(wěn)定流形而言，流形代表的是質(zhì)點離開周期軌道后的運動狀態(tài)，需要對Xu（X₀）進行正向積分;對于穩(wěn)定流形而言，流形代表的是質(zhì)點趨向周期軌道之前的運動狀態(tài)，需要對Xs（X₀）進行逆向積分。

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

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

r=（x－μ）²+y²

θ=α+π2－arctanv_yv_x（10）

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

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

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

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

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

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

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

假設(shè)擾動流形在截面Σ上的坐標集合為P，GEO的坐標集合為Q，P和Q上的點為

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

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

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

minxX_P－X_Q

s.t.

x=［X_P，X_Q］T

X_P∈P

X_Q∈Q

（12）

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

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

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

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

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

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

令W_dis在龐加萊截面上形成的截面圖為集合P，GEO軌道的龐加萊截面為Q。選取P和Q上距離最近的點X_P和X_Q，對于給定的時間T_l計算出Lambert轉(zhuǎn)移軌道，由此可以計算出Lambert弧段始末的速度增量Δv₂和Δv₃，之后航天器即可切入GEO。設(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）

能量時間最優(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è)計

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

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

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

C=2Ω－v²（18）

當航天器位置確定時，對式（18）兩邊同時計算變分，有

δC=－2vδv（19）

整理后可得

δv=－δC2v（20）

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

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

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

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

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

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

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

3.1.1 能量最優(yōu)

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

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

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

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

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

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

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

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

4 結(jié)束語

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

參考文獻

［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］ 徐明， 徐世杰. 地月系平動點及Halo軌道的應(yīng)用研究［J］. 宇航學報， 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點間轉(zhuǎn)移軌道設(shè)計［J］. 吉林大學學報（工學版）， 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］ 連一君. 基于三體平動點的低能轉(zhuǎn)移軌道設(shè)計研究［D］. 長沙： 國防科學技術(shù)大學， 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］ 張漢清， 李言俊. 地月三體問題下L1-地球低能轉(zhuǎn)移軌道設(shè)計［J］. 哈爾濱工業(yè)大學學報， 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］ 李翔宇， 喬棟， 程潏. 三體軌道動力學研究進展［J］. 力學學報， 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］ 錢霙婧， 荊武興， 劉玥， 等. 地月平動點擬周期軌道設(shè)計方法［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］ 李言俊， 張科， 呂梅柏， 等. 利用拉格朗日點的深空探測技術(shù)［M］. 西安： 西北工業(yè)大學出版社， 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.

作者簡介

喬琛遠（1999—），男，碩士研究生，主要研究方向為地月空間軌道動力學與控制。

楊樂平（1964—），男，教授，博士，主要研究方向為航天任務(wù)規(guī)劃、空間電磁操控。