自由返回軌道的P平面參數快速設計方法

劉 勇，劉 磊，曹鵬飛，張 堯

(1.北京航天飛行控制中心，北京 100094；2.航天飛行動力學技術重點實驗室，北京 100094)

0 引 言

2020年中國完成了月球采樣返回任務[1]，載人登月或將是未來月球探測的重要選項。作為載人登月任務關鍵技術的載人登月軌道設計，首要考慮的因素即航天員的安全，必須具備任務故障時將航天員安全送回地球的能力[2]，繞月自由返回軌道由此成為軌道設計的首選項。借助繞月自由返回軌道，飛船由地球飛抵月球后，可以經月球借力而無動力返回地球，從而實現應急救生目的。美國阿波羅13號飛船曾利用該類型軌道，在飛船發生故障后將3名航天員安全送回地球[3]。

在繞月自由返回軌道設計方面，雙二體模型、圓型限制性三體問題模型(Circular restricted three-body problem,CR3BP)等計算量較小，適用于軌道特性分析和初步軌道設計。Penzo等[4]在雙二體模型下研究了自由返回軌道特性。Miele等[5]基于CR3BP模型驗證了地月空間鏡像軌道的存在，并利用鏡像軌道理論對自由返回軌道進行優化設計。郗曉寧等[6]提出了一種基于雙二體模型圓錐曲線拼接的初步設計方法和全攝動模型精確修正的分層軌道設計方法。白玉鑄等[7]將月地轉移段作為應急返回軌道，只對近地點高度進行約束，在高精度模型下采用微分修正法求解了自由返回軌道。黃文德等[8]在雙二體模型下設計了載人登月自由返回軌道。張磊等[9]提出了橢圓軌道B平面參數，采用三級修正策略完成高精度模型下軌道求解，不過在軌道偏心率接近1時可能出現軌道計算結果的偏心率大于1，由此導致B矢量失效的問題。Peng等[10-11]采用粒子群算法和序列二次規劃(Sequence quadratic program,SQP)對自由返回軌道進行分步設計，同時研究了自由返回軌道在會合坐標系下的偏差傳播特性。劉林等[12]基于CR3BP模型構造了適用于4種類型自由返回軌道的設計方法，然后在真實引力模型進行了修正設計，計算速度快但約束考慮不全，只考慮高度約束沒有考慮傾角約束。Li等[13]綜合自由返回軌道和Hybrid軌道提出了多段自由返回軌道。賀波勇等[14]為降低月球引力甩擺段強非線性的影響，提出了由近月點參數出發逆/正求解自由返回軌道的思路。曹鵬飛等[15]進一步利用差分進化算法與SQP構造了一種自由返回軌道混合-分層優化策略，有效提升了計算精度，不過計算時間仍在小時量級。路毅等[16]以月球背面載人登月為背景，基于微分修正方法設計了自由返回軌道，研究了從自由返回軌道轉移至地月L2點Halo軌道的最優轉移問題。彭坤等[17]利用自由返回軌道對稱性建立了求解模型，研究了軌道設計方法和特性。周晚萌等[18]在文獻[14]基礎上，定義一組新的近月點狀態偽參數，提出一種基于混合多圓錐方法的自由返回軌道設計方法，提高了初值精度。He等[19]以基于近地軌道(LEO)交會的載人月球任務為背景，提出了一種基于高精度動力學模型的自適應LEO相位自由返回軌道設計方法。

本文針對B平面參數在偏心率等于1時不連續導致軌道求解收斂困難等問題，基于圓錐曲線半通徑，提出了P平面參數的概念并將其應用到自由返回軌道的求解中。首先，對P平面參數的定義和計算方法進行闡述，并將其與B平面參數進行了連續性對比分析。其次，研究了基于P平面參數的自由返回軌道求解方法，并在簡化模型和高精度力模型下設計了多種類型的自由返回軌道，驗證了方法的普適性與有效性。最后，通過大量仿真算例分析了自由返回軌道特性。研究結果可直接應用于中國后續月球探測任務軌道設計。

1 動力學模型

在自由返回軌道設計中，通常的思路是首先使用二體模型、雙二體模型和限制性三體模型等簡化力模型設計初步軌道，進而使用考慮攝動因素的高精度力模型修正初步軌道[20]。

在地心或月心慣性坐標系中，考慮日地月質點引力、地月非球形攝動和大氣阻力等攝動力，探測器的動力學方程為

(1)

式中：r為地心或日心位置矢量；μ為中心體引力常數；AN為第三體引力攝動加速度，可通過JPL DE421星歷計算得到第三體位置；ANSE為地球非球形攝動，取WGS84引力場模型6×6階計算；ANSM為月球非球形攝動，取LP165P引力場模型6×6階計算；AR為太陽光壓攝動；AD為大氣阻力攝動加速度。

2 P平面參數

文獻[9]對B平面參數優勢與局限性進行了充分了論證，這里不再贅述。本節提出P平面參數，并給出定義、計算方法以和線性度分析。相比B平面參數，P平面參數應用范圍更加廣泛，可以適用于雙曲線、橢圓和拋物線軌道設計。

2.1 P平面參數定義

半通徑是圓錐曲線的通用參數，常用字母p表示。對于橢圓軌道、拋物線軌道和雙曲線軌道，p均是有限值，且在偏心率跨越1時保持連續，與短半軸b相比，具有較好的連續性。

如圖1所示，給出了橢圓軌道P平面示意圖。由圖可知，P平面參數選用通用的近心點e矢量代替雙曲漸近線的S矢量，選用P矢量代替B矢量，P平面即為過天體中心且垂直于e矢量的平面。對于拋物線軌道和雙曲線軌道，上述定義方法同樣適用，不再贅述。下面給出P平面參數的求解方法。

(2)

(3)

圖1 橢圓軌道的P平面Fig.1 P plane for the elliptical orbit

P矢量可表示為

(4)

P平面參數為

(5)

可見，P平面參數在幾何意義上與B平面參數類似，但是由于半通徑對于圓錐曲線軌道普遍適用，因而應用范圍更廣。

2.2 目標P平面參數計算

上文給出了根據探測器當前狀態計算實際P平面參數的方法。在實際應用中，需要根據軌道傾角和近心距，或再入點高度和再入角等目標參數計算目標P平面參數，以此作為微分修正目標量。

引入假設：假設目標軌道e矢量與半長軸a不變，給出目標P平面參數求解方法，如圖2所示。圖中α為P平面與參考平面的夾角，為

(6)

根據半長軸a與目標近心距rp，可計算出半通徑p，即

(7)

當目標參數為再入點地心距re和再入角γ時，

(8)

(9)

式中：φ的兩個解分別對應升軌和降軌，當cosi>sinα時式(9)無解，由此限定了目標軌道傾角的取值。

根據式(8)和式(9)，可得P平面參數如下

(10)

圖2 P平面參數求解Fig.2 Solution of the P-plane parameters

2.3 P平面參數連續性分析

雙曲線軌道B平面參數使用時要求偏心率大于1.0。擴展的橢圓B平面參數[9]，適用于橢圓軌道和雙曲線軌道，但其在偏心率為1時無效。設置月心J2000系近月點位置速度：x=1940.9 km，vy=1.124 km/s，vz=1.947 km/s，其他位置和速度分量為0。改變Y方向速度分量，使軌道偏心率在0.4～1.8之間連續變化，并計算對應的B平面參數和P平面參數的值。得到偏心率與B平面參數的關系和偏心率與P平面參數的關系，分別如圖3和圖4所示。

由圖3可知，在e=1處，B平面參數存在明顯的大值。因此，當使用擴展橢圓B平面參數作為目標參數進行微分修正時，迭代過程中易存在發散的情況，導致不能獲得期望的結果，因此需要設計新的目標參數。

圖3 擴展橢圓B平面參數隨偏心率的變化Fig.3 Variation of the extended elliptic B-plane parameters with eccentricity

由圖4可知，當圓錐曲線在橢圓、拋物線、雙曲線之間變換時，P矢量是連續的，不存在奇點或者無定義的情況。

圖4 P平面參數隨偏心率的變化Fig.4 Variation of the P-plane parameters with eccentricity

由上述分析可知，與擴展橢圓B平面參數相比，P平面參數可以適應偏心率小于1、等于1和大于1的情況，具有良好的連續性。另外，對于行星際引力輔助[22]或地火轉移至火星捕獲時，也可能遇到剩余速度接近于0(即偏心率接近于1)的情況，在進行軌道設計時，該方法同樣適用。

3 自由返回軌道快速設計

本節將在上一節的基礎上，以近月點參數為出發點，研究基于P平面參數的自由返回軌道快速設計方法，主要包括白道面內對稱自由返回軌道設計和空間繞月自由返回軌道設計。

3.1 白道面內自由返回軌道設計

作為最簡單的自由返回軌道類型，白道面內的對稱自由返回軌道可以為地月自由返回軌道設計提供良好初值，因此給出基于P平面參數的白道面內對稱自由返回軌道設計方法。

定義瞬時地月慣性系，與探測器近月點時刻的白道坐標系平行，如圖5所示。選擇月心瞬時地月慣性系X軸上的近月點作為初始點，位置速度記為(x0,0,0,vy)，逆向積分設計地月轉移軌道，正向積分設計月地返回軌道，在x0給定時，僅vy為設計變量。

圖5 白道面內對稱自由返回軌道狀態Fig.5 Status of the symmetrical free-return orbit in the plane of lunar orbit

根據x0和vy，易求得白道面內的月心軌道參數，參考文獻[21]，以及出影響球邊界的徑向速度vr和垂直徑向速度vu為

(11)

式中：θ為影響球邊界點的真近點角。

根據θ計算出飛行時間Δt，從而得到月球運行的角度β=ωLΔt，以及地心瞬時地月慣性系下月球的位置矢量RL和速度矢量VL

(12)

式中：L是地月距；V是月球繞地球的飛行速度。

于是可得探測器在月心瞬時地月慣性系中的位置矢量Rm和速度矢量Vm為

(13)

(14)

式(13)和式(14)分別為初始點在月球外側和內側的情況，Rs為影響球半徑。

于是，探測器的地心位置矢量Re和速度矢量Ve

(15)

進而得到地心軌道的半通徑pe參數和近地點地心距rpe

(16)

(17)

令Re的分量為xe和ye，Ve的分量為vex和vey，則式(16)和(17)中

(18)

(19)

以上分析的是出影響球邊界情況，進入影響球邊界情況與之相同，且二者軌道關于地月軸對稱。

數值分析rpe與x0和vy的關系，如圖6所示，各曲線起點的vy均為月心距x0時的月球逃逸速度。

圖6 rpe與x0和vy的關系Fig.6 Curve of relationship between rpe and x0 and vy

當x0為正值或負值時，自由返回軌道相對月心分別為逆行或順行軌道，結合返回軌道相對于地心的運動方向，白道面內自由返回軌道共有4種構型。

3.2 空間繞月自由返回軌道設計

實際任務中，無論是根據任務需要還是受攝動因素影響，繞月自由軌道都不可能僅位于白道面內，因此基于P平面參數設計空間繞月自由返回軌道，軌道力模型可以為雙二體模型或高精度力模型。

在瞬時地月慣性系下，空間繞月自由返回軌道初始點選取在xoz平面內，此時y0=0，在給定坐標x0后設計變量為[z0,vx,vy,vz]T。軌道計算方法與白道面內繞月自由返回軌道設計相同，即由初始點逆向積分計算得到入軌參數，正向積分計算得到返回參數，然后根據二者與設計目標偏差修正設計變量，直到偏差滿足設計要求。

空間繞月自由返回軌道計算步驟如下：

1)給定初始時刻T0和x0；

2)根據T0計算月球位置速度，建立瞬時地月慣性系，并得到瞬時地月慣性系到地心慣性系的轉換矩陣Lia；

3)采用3.1小節方法計算白道面內對稱繞月自由返回軌道的vy 0，得到設計變量初值[0,0,vy 0,0]T；

4)修正[z0,vx,vy,vz]T使得軌道參數滿足設計要求，得到包含入軌點、近月點和返回點的繞月自由返回軌道。

從計算效率角度，步驟(4)可以直接采用牛頓迭代法或Broyden法等微分改正方法，計算過程如下：

1)將初始點位置速度[x0,0,z0,vx,vy,vz]T轉換成軌道根數；

2)求解進入影響球時的真近點角θi；

3)根據θi和T0計算進影響球的時刻Ti與月心瞬時地月慣性系得位置速度RA和VA；

4)計算Ti時刻月球位置速度RL和VL，以及對應地心慣性系的位置速度Re和Ve；

5)根據Re和Ve計算入軌時刻的P平面參數和偏差；

6)根據真近點角θo=-θi和T0計算進影響球的時刻Tc，然后重復步驟3)～5)計算返回或再入時刻的P平面參數與偏差;

7)根據P平面參數偏差，采用微分改正算法修正設計參數。

若目標軌道參數包括入軌點經度和返回軌道經度，可以將T0和x0與前述4個參數一起作為設計變量，構建一個6對6的計算模型，可對T0和x0做如下調整

(20)

發射和返回的時間差ΔTi和ΔTc根據軌道與發射點和落點的經度差計算。

3.3 數值計算討論

由于繞月自由返回軌道繞飛地球和月球，動力學環境變化復雜，數值積分時采用自適應變步長積分器較好。

對于前往地球的轉移軌道，可以采用積分到近地點方式，同時為了減小計算消耗，可以根據瞬時軌道用二體模型估算到達近地點的飛行時間，然后限制積分最大時間低于估算時間，重復估算和積分，直到近地點為止。在上述積分計算之前，最好先固定積分1天以飛出月球的引力范圍。與用近地點時刻作為設計變量相比，積分到近地點的方法可以減少兩個變量，由此減小計算量。

基于P平面參數的繞月自由返回軌道設計具有較好的收斂性，高精度模型下的軌道設計可以直接使用白道面內雙二體解作為初值，甚至可以直接使用模值較大的vy 0作為初值。

由于地心和月心附近軌道分別存在對應升軌和降軌的P平面參數，再結合地心和月心處的順行和逆行方向，空間繞月自由返回軌道共有16種構型。

4 數值仿真

入軌點傾角設置為28.5°即地心順行軌道，近地點高度設置為200 km，返回點設置為傾角43.0°，再入點高度設置為120 km，由此存在8種構型的地月自由返回軌道。以北京時2022年8月9日07∶58∶50的月球位置速度建立瞬時地月慣性系，x0分別取±6000 km。設計得到雙二體模型和高精度力模型下的自由返回軌道分別如表1和表2所示。

表1和表2中的軌道都為近地點的軌道參數，前4種構型的x0為6000 km，即月心逆行軌道，其總飛行時間約6天，后4種構型的x0為-6000 km，即月心順行軌道，其總飛行時間約14天。對比表1和表2可見，雙二體模型的偏差并不大，因此可以用雙二體模型進行初步窗口搜索和分析。同時可見，如果發射載人登月軌道，應該采樣月心逆行軌道。對應飛行軌跡如圖7和圖8所示。

進一步，分析各種軌道類型對應的x0與轉移時間關系，如圖9所示。月心順行軌道的轉移時間隨|x0|單調下降，月心逆行軌道的轉移時間則隨|x0|單調上升。x0對總轉移時間影響很大，特別是近月距較低時的月心順行軌道轉移時間變化劇烈，而月心逆行軌道轉移時間變化較緩。當x0在26000 km左右時，月心順行和逆行軌道的轉移時間接近。

表1 雙二體模型下的自由返回軌道Table 1 Free return trajectories with the dual two-body model

表2 高精度力模型下的自由返回軌道Table 2 Free return trajectories with the high-precision dynamics model

為了驗證計算效率，采用速度為2.3 GHz的單核CPU運行C語言編寫的相關代碼，求解一條空間雙二體模型繞月自由返回軌道僅約耗時0.2 ms，至于考慮日地月引力和地球J2項的高精度力模型，采用KSG積分器，求解一條繞月自由返回軌道僅約耗時0.5 s。在相似計算條件下，應用文獻[15]的方法計算同樣的轉移軌道則需要至少1個小時，因此基于P平面參數的自由返回軌道設計方法具有非常高的計算效率。

圖7 月心逆行的4種構型自由返回軌道三維視圖Fig.7 Three-dimensional view of free return trajectories in 4 lunar centric retrograde configurations

圖8 月心順行的4種構型自由返回軌道三維視圖Fig.8 Three-dimensional view of the free return orbit in 4 lunar centric anterograde configurations

圖9 轉移時間與|x0|關系曲線Fig.9 Relation curves between transfer time and |x0|

5 結 論

本文提出了一種基于P平面參數的自由返回軌道設計方法，研究結果表明，P平面參數不僅具有類似于B平面參數線性特征，而且適用于各種類型的圓錐曲線軌道，因而適用性更廣。基于P平面參數的繞月自由返回軌道設計方法，具有軌道構型清晰直觀，求解計算高效，尤其適合繞月自由返回軌道設計，以及其他天體接近軌道設計。此外，從數值仿真結果可見，月心逆行自由返回軌道總飛行時間明顯低于月心順行軌道，可作為未來載人登月任務首選軌道。