基于月球借力的低能DRO 入軌策略

張晨，張皓

中國科學院 空間應用工程與技術中心，北京 100094

隨著航天技術的發展，人類的生存空間逐漸由近地空間拓展到地月空間，2017 年美國國家航空航天局（National Aeronautics and Space Administration，NASA）提出的“ 阿爾忒彌斯（Artemis）”任務［1］將首次實現人類在月球及其軌道上的長期居留，并為未來復雜的載人火星任務奠定基礎。2021 年中國與俄羅斯就建設國際月球科研站達成了合作協議［2］。可以預見在未來很長一段時間內，地月空間都將是人類航天活動的主要內容。在地月空間部署大型月球軌道站可以顯著降低月球開發成本并作為其他深空任務的跳板［3-5］，停泊軌道的選擇不能以單方面性能作為評價指標，而應權衡各種任務場景以達到綜合性能最優。其中近直線暈軌道（Near Rectilinear Halo Orbit，NRHO）和遠距離逆行軌道（Distant Retrograde Orbit，DRO）是兩類適合部署月球軌道站的三體周期軌道族，這兩類軌道具有一定的相似性，但側重點略有不同，NRHO 適合用于載人登月任務，而DRO 更能兼顧后續的小行星及火星任務，本文僅對基于DRO 的月球軌道站部署策略展開研究。

近年來多個任務圍繞地月三體系統中的DRO軌道展開，這類軌道具有長期穩定性好和入軌代價低的特點。例如“小行星重定向任務（Asteroid Redirect Mission，ARM）［6］”計劃從小行星上抓取一塊巨石并拖曳至DRO 軌道上，DRO 當中的某些軌道具有長達100 a 的穩定性而不需要進行頻繁的軌道維持。Artemis 的首星［7］將把獵戶座飛船送入25.5 d 的任務軌道，其中6 d 將位于DRO軌道。面向未來DRO 月球軌道站的在軌建造和貨物補給任務，提高DRO 入軌質量是重要問題。

1987 年Belbruno［8］提出的弱穩定邊界（Weak Stability Boundary，WSB）轉移技術在深空任務中獲得巨大成功，WSB 轉移依靠多個引力體的相互作用，利用混沌區域所展現的一些新的動力學特性降低飛行器入軌脈沖。WSB 轉移最早用于奔月軌道設計，通過選擇合適的地球發射條件，衛星首先抵達日地引力混沌的弱穩定邊界（大約3～5 倍地月距），衛星在這個區域無動力飛行時，太陽引力將近地點從地球附近抬升至月球軌道附近，當衛星從外部返回月球時被低能捕獲［9］，盡管這種方法需要80～120 d 的飛行時間，但是能夠節省大約25%的月球制動脈沖。1991年WSB 轉移首次應用于日本Hiten［10］衛星的營救任務，隨后GRAIL（Grarity Recovery and Interior Laboratory）［11］、ESMO（the European Student Moon Orbiter）［12］、CAPSTONE（Cislunar Autonomous Positioning System Technology Operations and Navigation Experiment）［13］、KPLO（Korea Pathfinder Lunar Orbiter ）［14］等任務都采用或計劃采用WSB 轉移技術。

國內外學者圍繞WSB 轉移軌道的設計方法展開了內容豐富的研究，Koon 等［15］通過平動點附近的不變流形對WSB 的轉移機理進行了解釋，并采用雙三體流形拼接方法構建了低能奔月軌道。Yagasaki［16］在四體簡化模型下將地月WSB 轉移問題轉換成非線性規劃問題，通過數值迭代求解了地月低能轉移軌道。WSB 轉移不僅適用于地月低能轉移，也適用于DRO、NRHO、L4/5等地月空間典型周期軌道。Xu M 和Xu S J［17］、Tan 等［18］利 用L1/2附近不變流形設計了基于WSB 轉移的DRO 入軌策略。Scheuerle等［19］基于多步打靶和和數值延拓在四體簡化模型下計算了地球至DRO 的低能轉移軌道。Zhang 和Hou［20］在雙圓四體模型下設計了地球至L4/5的低能轉移軌道。

為了最大化飛行器入軌質量，一些學者指出在WSB 轉移軌道的地球出發段增加月球借力（Lunar Gravity Assist，LGA）能夠有效降低火箭發射能量。為方便描述，將這種類型的轉移軌道稱為“LGA+WSB”轉移方式。“LGA+WSB”最早出現于Belbruno 和Miller 在1993 年發表的文章［10］。Topputo［21］在雙圓四體模型下探索了兩脈沖地月低能轉移軌道的全局解空間，并對不同轉移方式進行了分類，指出“LGA+WSB”轉移方式具有較低的火箭發射脈沖。Parrish 等［13］針對“CAPSTONE”任務對2 種NRHO 入軌方案進行了對比，其中“WSB”方案火箭發射脈沖較高但是發射窗口寬松（每天都有發射窗口），“LGA+WSB”方案火箭發射脈沖較低但是發射窗口緊張（一個月一次）。上述研究沒有針對敏感的月球借力進行單獨設計和優化，轉移軌道中僅有很少一部分呈現“LGA+WSB”的轉移方式，數值計算效率較低。Scheuerle 等［19］在雙圓四體模型下從“WSB”延拓得到“LGA+WSB”轉移軌道，但沒有在星歷模型下進行修正。Tselousova 等［22］在雙圓四體模型下基于軌道拼接的思路設計了“LGA+WSB”轉移軌道，設計過程較為繁瑣且沒有在星歷下實現。

為了提高星歷模型下的“LGA+WSB”轉移軌道設計效率，為了避免復雜動力學環境對轉移軌道收斂性所帶來的影響，使用“近月點龐加萊圖”和“v無窮匹配”獲得轉移軌道初值，再通過“多步打靶”在星歷下修正轉移軌道。改進的方法充分優化了月球借力參數，有效改善了“LGA+WSB”轉移方式的計算效率，所得到的轉移軌道滿足地球發射約束和高精度動力學約束。

本文結構如下：第1 節介紹火箭發射能量和載荷質量的關系；第2 節為動力學模型和狀態轉移矩陣；第3 節介紹DRO 軌道；第4 節介紹“LGA+WSB”轉移方式的特點；第5 節介紹LGA 至DRO 軌道段的設計方法；第6 節介紹近地球軌道（Lwo Earth Orbit，LEO）至LGA 軌道段的設計方法；第7 節介紹多步打靶技術；第8 節為數值仿真；第9 節為結論。

1 火箭發射能量和載荷質量

假設衛星初始位于近地圓形停泊軌道，沿切向施加脈沖抬升遠地點高度。衛星在近地點的速度模表示為

式中：ra和rp分別為遠地點和近地點地心距；μE為地球引力常數。衛星在近地點的脈沖為

發射能量C3表示為

假設衛星近地點高度為200 km，圖 1 展示了LEO 發射脈沖Δνdep和發射能量C3隨遠地點地心距ra的變化曲線。當發射脈沖Δνdep=3.131 3 km/s時遠地點能夠抵達月球（ra約為4×105km），發射能量C3約為-2.03（km）2/s2。發射脈沖為［3.194 3，3.206 3］ km/s 時遠地點抵達WSB 區域ra約為［1.2，2.0］×106km，發射能量C3約為［-0.66，-0.39］（km）2/s2。

圖1 Δvdep 和C3 隨遠地點高度變化曲線Fig 1 Δvdep and C3 according to apogee radius

采用國外公開的火箭運力曲線［23］，圖 2 展示了獵鷹9 Block 2 火箭載荷質量隨發射能量的變化曲線。獵鷹9 發射地月轉移軌道的載荷質量為2 619.77 kg，直接發射至WSB 的載荷質量為［2 500.24，2 519.34］ kg。因而如果先將衛星發射至地月轉移軌道，再借助月球借力將遠地點抬升至WSB，則能夠節省60～70 m/s 的速度脈沖，提升載荷質量100.43～119.53 kg。

圖2 獵鷹9 載荷質量隨C3 變化曲線Fig 2 Rocket capability according to characteristic energy C3 of Falcon 9 Block 2

2 星歷模型和狀態轉移矩陣

在深空任務設計時，經常使用高精度的N體星歷模型對轉移軌道進行設計和修正。圖 3 展示了N體星歷模型示意圖，其中k為中心天體，i為衛星，j為其他攝動天體。定義衛星狀態x=[rki，vki]?，衛星的動力學方程表示為

圖3 N 體星歷模型Fig 3 N-body ephemeris model

式中：rki為衛星i相對于中心天體k的位置矢量；vki為衛星i相對于中心天體k的速度矢量；rkj為攝動天體j相對中心天體k的位置矢量，可以通過JPL 的DE430 星歷獲得；rji為衛星i相對攝動天體j的位置矢量，表示為

可以發現攝動天體的位置矢量隨歷元變化，因而N體動力學模型是歷元的函數。

軌道修正經常使用狀態轉移矩陣這一概念，狀態轉移矩陣在一階程度上刻畫了初始狀態改變量對終端狀態改變量所產生的影響［24］。狀態轉移矩陣的動力學方程表示為

式中：I為單位矩陣；Φ為狀態轉移矩陣；f(x)為衛星動力學方程；A為動力學方程對狀態求偏導得到的矩陣；t0為初始時刻，t為轉移時刻。定義y為狀態變量x和狀態轉移矩陣Φ構成的列向量，其動力學方程表示為

其中：g(y)為衛星狀態和狀態轉移矩陣的動力學方程；“vec”是將矩陣變為列向量的算子。

3 DRO 軌道

DRO 是圓型限制性三體問題（Circular Restricted Three-Body Problem，CRTBP）中一類穩定的平面軌道族［25-28］，在旋轉坐標系下該軌道沿順時針（逆行）運動。圖 4 展示了DRO 軌道族，軌道顏色用于區分軌道周期，L1～L5為5 個平動點，圖中距離單位1 LU=384 400 km。可以發現振幅越小的DRO 軌道周期越短。DRO 共振比是指軌道周期與月球公轉周期之比，具有典型共振比的周期軌道通常可保證衛星與地球、月球具有周期性的幾何關系，在數值仿真中用到了共振比2∶1 的DRO 軌道。

圖4 CRTBP 下 的DRO 軌道族Fig 4 DRO family in CRTBP

以CRTBP 的DRO 軌道作為初值，使用多步打靶可獲得星歷下的DRO 軌道，定義為DRO 參考狀態，則在星歷模型下數值積分可以得到DRO 軌道上任意一點的狀態xdro，表示為

式中：φeph表示以為初值，從τ*到τ積分星歷模型所得到的解。

4 轉移軌道設計

基于月球借力的低能DRO 入軌策略可簡單描述為，衛星初始位于近地圓形停泊軌道，施加第1 次脈沖后進入地月轉移軌道（Lunar Transfer Orbit，LTO）并掠飛月球，衛星通過月球借力抬升遠地點高度并改變軌道傾角，之后抵達地月3～5 倍距的弱穩定邊界，衛星在這里持續受到太陽引力的影響并抬升近地點高度，當再次返回地月空間時施加第2 次脈沖并進入DRO 軌道，希望優化軌道轉移策略使得兩次脈沖和最小。圖 5 展示了“LGA+WSB”的DRO 入軌示意圖，其中，Δvarr為DRO 入軌速度脈沖矢量；Δvdep為LEO 發射脈沖矢量。

圖5 LGA+WSB 的DRO 入軌示意圖Fig 5 Illustration of LGA+WSB transfer into DRO

“LGA+WSB”的轉移方式能夠有效降低任務總脈沖，但是復雜的動力學環境和月球借力給轉移軌道設計帶來解空間龐大和數值敏感等問題。本文將分別設計“LGA 至DRO 軌道段”和“LEO 至LGA 軌道段”，再使用多步打靶技術修正整條轉移軌道。

5 LGA 至DRO 軌道段

本節將構造近月點龐加萊圖，使得轉移軌道從DRO 逆向積分先抵達WSB 再抵達近月點。圖 6 展示了DRO 入軌時刻示意圖，衛星在DRO入軌時刻的狀態x(τf)表示為

圖6 DRO 入軌示意圖Fig 6 Illustration of DRO insertion

式中：xdro=[rdro，vdro]?為入軌時刻DRO 的狀態；rdro為入軌時刻DRO 位置矢量；vdro為入軌時刻DRO 速度矢量；Δvarr為入軌脈沖矢量；Δνarr為入軌脈沖矢量的模；α∈[0，2π)為Δvarr在e1—e2平面的投影與e1的夾角；β∈[ -π/2，π/2]為Δvarr與e1—e2平面的夾角。

圖7 展示了DRO 出發的近月點龐加萊圖，轉移軌道遠地點抵達WSB 的約束表示為

圖7 DRO 出發的近月點龐加萊圖Fig 7 Illustration of perilune Poincare map of DRO

式中：v為衛星速度矢量；r為衛星位置矢量；a為衛星加速度。轉移軌道抵達近月點的約束表示為

式中：rM、vM和aM分別為月球的位置、速度和加速度。受火箭滑行時間的限制，一些任務只能在月球升/降交點附近進行引力輔助，轉移軌道的近月點還需滿足以下約束

由于衛星以雙曲線軌道掠飛月球，衛星離開月球影響球的雙曲線剩余速度v+∞表示為

式中：x(τlga)為衛星在近月點的狀態；xM(τlga)為月球狀態為衛星相對月球狀態。

圖8 月球借力和B 平面Fig 8 Illustration of lunar gravity assist and B-plane

雙曲線剩余速度的模表示為

式中：μM為月球引力常數；νˉ為衛星相對月球速度矢量的模；rˉ為衛星位置矢量的標量。偏心率矢量e表示為

其中：e為偏心率的標量大小。

衛星離開影響球的雙曲線剩余速度v+∞表示為

6 LEO 至LGA 軌道段

構造LEO 至LGA 的轉移軌道，期望在降低地球發射脈沖的同時，優化月球借力參數，使得借力后的雙曲線剩余速度與相匹配。在二體簡化模型下構造全局優化問題為

式中：z∈Rn是優化變量；J：Rn→R是目標函數；zlr∈Rn和zur∈Rn是優化變量邊界。優化變量z表示為

式中：Ω和ω分別為LTO 的升交點赤經和近地點幅角；η是LTO 的飛行時間；γ和ρ分別為月球引力輔助的相位角和月心距。表 1 列出了優化變量和取值范圍。

表1 LEO 至LGA 軌道段優化變量Table 1 Variables of LEO to LGA segment

假設地球圓形停泊軌道的地心距為R*，軌道傾角為i*，通過Ω和ω得到衛星在LTO 的出發狀態x(τ1)，基于二體圓錐曲線拼接假設，求解Lambert 問題得到地球出發脈沖Δνdep和衛星進入月球影響球的雙曲線剩余速度圖 9 展示了B 平面和月球借力參數，其中i1，i2，i3為坐標系的單位矢量。通過γ和ρ得到月球借力后的雙曲線剩余速度，計算過程如下：

圖9 B 平面和月球借力參數Fig 9 B-plane and lunar gravity assist parameters

為了降低地球出發脈沖Δνdep并使得盡可能接近，構造目標函數

需要說明的是，式（32）使得月球借力前后的2 段軌道通過v無窮進行匹配（不是位置速度拼接）。當全局優化問題收斂后，可得到星歷模型下的“LGA+WSB”轉移軌道初值，計算過程不再贅述。

7 多步打靶修正

多體系統復雜的動力學環境使得轉移軌道對初值非常敏感，多步打靶技術可以擴大收斂域并提高算法的魯棒性。圖 10 展示了多步打靶示意圖，整條軌道被n個離散點分割為n-1 條軌道段，多步打靶變量X表示為

圖10 多步打靶示意圖Fig 10 Illustration of multiple shooting

第j段軌道滿足狀態連接約束：

第j段軌道滿足時間連接約束：

第j段軌道滿足飛行時間為正的不等式約束：

引入松弛變量βj可將式（36）變為等式約束：

第1 個節點滿足邊界約束：

式中：R*和i*為LEO 的地心距和軌道傾角。

第n個節點滿足邊界約束：

至此，當得到雅克比矩陣J(X)，使用最小二乘更新方程對打靶變量Xk+1進行修正（見式（41）），收斂后即可得到轉移軌道。在星歷修正中，動力學模型僅考慮日、地、月引力場，行星狀態通過JPL 的DE430 星歷計算。

衛星在任務過程中共施加兩次脈沖，分別為火箭發射脈沖Δvdep和衛星平臺施加的DRO 入軌脈沖Δνarr，任務脈沖總表示為

8 數值仿真

將在星歷模型下設計基于月球借力的低能DRO 轉移軌道，仿真參數參見表 2。其中，時間單位中的jd 為儒略日，tdb 為太陽系質心動力學時。計算流程如下：

表2 數值仿真參數Table 2 Simulation parameters

1）網格化DRO 入軌時間τf和入軌速度脈沖Δvarr，獲得不同的DRO 入軌狀態x(τf)，參見式（9）和式（10）。

2）從x(τf)出發逆向積分動力學模型，參見式（8），使得衛星先抵達WSB 再抵達近月點（構造近月點龐加萊圖），參見式（13）～式（18），計算月球借力后的雙曲線剩余速度，參見式（20）～式（24）。

3）固定LGA 至DRO 軌道段，使用粒子群優化LEO 至LGA 軌道段，最小化地球出發脈沖Δνdep并使得盡可 能接近（即v無窮 匹配），參見式（25）～式（32）。

4）使用多步打靶修正，得到“LGA+WSB”類型的DRO 轉移軌道，參見式（33）～式（41）。

根據以上計算流程，構建1×107個DRO 入軌狀態，從DRO 逆向積分能夠先抵達WSB 再抵達近月點的僅有1 511 條轉移軌道。圖 11 展示了近月點龐加萊圖，圖中每個圓點都代表一條轉移軌道，圓點顏色用于區分DRO 入軌脈沖。

圖11 近月點龐加萊圖Fig 11 Perilune Poincare map

進一步對比“v無窮匹配”策略對計算效率的改進效果。首先對于圖 11 上的每一個點，繼續逆向積分并在近地點終止，得到不使用“v無窮匹配”的轉移軌道初值。圖 12 展示了這些初值的近地點龐加萊圖，X軸為近地點地心距，Y軸為近地點軌道傾角，紅色十字為期望的LEO 地心距和軌道傾角，可以發現敏感的月球借力使得紅色十字附近幾乎沒有接近的轉移軌道。如果對圖 12 直接使用多步打靶軌道修正，則僅有49 條收斂的轉移軌道。

圖12 近地點龐加萊圖Fig 12 Perigee Poincare map

接著采用“v無窮匹配”策略獲得轉移軌道初值，再使用多步打靶修正整條軌道，最終獲得了743 條收斂的轉移軌道。圖 13 展示了這些收斂解的散點圖，其中每個圓點都為一條收斂軌道，圓點顏色用于區分DRO 入軌脈沖，紅色圓圈為收斂解的帕累托前沿，紅色五角星為脈沖最低解。在數值仿真中采用Matlab 編程語言，計算設備為高性能工作站，CPU 為64 核心處理器，內存為256 G，計算總時間約為8 h。

圖13 多步打靶修正后的收斂解Fig 13 Converged solutions with multiple shooting

需要說明的是，“v無窮匹配”策略能夠在不改變月球借力后雙曲線剩余速度的條件下，將圖12 的大多數近地點都移動到紅色十字附近，從而有效改進了轉移軌道初值。就本算例而言，“v無窮匹配”策略將轉移軌道收斂率提高了約15 倍。此外，“v無窮匹配”策略可以方便拓展到NRHO、L4/5等周期軌道的入軌任務設計，計算流程并無明顯區別。

圖14 展示了圖 13 中總脈沖最低解（紅色五角星）的轉移軌道，其中紅色實線為DRO 軌道，藍色實線為地球至DRO 轉移軌道，灰色實線為月球軌道。觀察圖 14（a），衛星從地球出發后進入地月轉移軌道，月球借力后衛星彈射至WSB，當再次返回地月空間時，衛星以切向入軌的方式進入DRO 軌道。觀察圖 14（b），衛星從WSB 返回地月空間時首先抵達L4附近，之后在地月引力共同影響下抵達DRO 并入軌。

圖14 總脈沖最低解的轉移軌道Fig 14 Best solution with minimum cost

表3展示了圖 13中脈沖最低解（紅色五角星）的仿真結果，LEO 出發脈沖僅需3.127 3 km/s，DRO入軌脈沖66.1 m/s，任務總脈沖3.193 4 km/s，任務總時間102.880 3 d。

表3 總脈沖最低解的仿真結果Table 3 Best solution with minimum cost

表4 展示了LEO 轉移至GEO 和DRO 的 對比，LEO 軌道高度200 km，軌道傾角28.5°，衛星平臺采用化學推進，比沖220 s。如果采用獵鷹9號Block 2 型運載火箭，單次發射至多可將1 704.273 3 kg 的載荷送入GEO 軌道，可將2 553.576 8 kg 的載荷送入DRO 軌道，由于GEO 需要1.819 7 km/s 的速度脈沖改變軌道傾角和軌道圓化，使得雖然DRO 距離地球更遠，但是入軌質量卻更高，約為GEO 的1.5 倍。

表4 LEO 轉移至GEO 和DRO 的對比Table 4 Comparison between LEO to GEO and LEO to DRO

9 結論

1）改進了星歷模型下的“LGA+WSB”轉移軌道設計方法，通過“近月點龐加萊圖”和“v無窮匹配”獲得較好的軌道初值，有效提高了該類型轉移軌道的計算效率。

2）通過網格搜索獲得了解空間的帕累托前沿，對于2：1 DRO 軌道，總脈沖最低解的飛行時間102.880 3 d，地球發射脈沖3.127 3 km/s，DRO 入軌脈沖66.1 m/s。

3）對 比LEO 至GEO 和LEO 至DRO 轉 移軌道，盡管DRO 距離地球更遠，但是入軌質量可以達到GEO 的約1.5 倍。