自由返回軌道的P平面參數(shù)快速設(shè)計方法

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

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

0 引 言

2020年中國完成了月球采樣返回任務(wù)[1]，載人登月或?qū)⑹俏磥碓虑蛱綔y的重要選項。作為載人登月任務(wù)關(guān)鍵技術(shù)的載人登月軌道設(shè)計，首要考慮的因素即航天員的安全，必須具備任務(wù)故障時將航天員安全送回地球的能力[2]，繞月自由返回軌道由此成為軌道設(shè)計的首選項。借助繞月自由返回軌道，飛船由地球飛抵月球后，可以經(jīng)月球借力而無動力返回地球，從而實現(xiàn)應(yīng)急救生目的。美國阿波羅13號飛船曾利用該類型軌道，在飛船發(fā)生故障后將3名航天員安全送回地球[3]。

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

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

1 動力學(xué)模型

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

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

(1)

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

2 P平面參數(shù)

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

2.1 P平面參數(shù)定義

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

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

(2)

(3)

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

P矢量可表示為

(4)

P平面參數(shù)為

(5)

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

2.2 目標(biāo)P平面參數(shù)計算

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

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

(6)

根據(jù)半長軸a與目標(biāo)近心距rp，可計算出半通徑p，即

(7)

當(dāng)目標(biāo)參數(shù)為再入點地心距re和再入角γ時，

(8)

(9)

式中：φ的兩個解分別對應(yīng)升軌和降軌，當(dāng)cosi>sinα?xí)r式(9)無解，由此限定了目標(biāo)軌道傾角的取值。

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

(10)

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

2.3 P平面參數(shù)連續(xù)性分析

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

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

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

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

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

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

3 自由返回軌道快速設(shè)計

本節(jié)將在上一節(jié)的基礎(chǔ)上，以近月點參數(shù)為出發(fā)點，研究基于P平面參數(shù)的自由返回軌道快速設(shè)計方法，主要包括白道面內(nèi)對稱自由返回軌道設(shè)計和空間繞月自由返回軌道設(shè)計。

3.1 白道面內(nèi)自由返回軌道設(shè)計

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

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

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

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

(11)

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

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

(12)

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

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

(13)

(14)

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

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

(15)

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

(16)

(17)

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

(18)

(19)

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

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

圖6 rpe與x0和vy的關(guān)系Fig.6 Curve of relationship between rpe and x0 and vy

當(dāng)x0為正值或負(fù)值時，自由返回軌道相對月心分別為逆行或順行軌道，結(jié)合返回軌道相對于地心的運動方向，白道面內(nèi)自由返回軌道共有4種構(gòu)型。

3.2 空間繞月自由返回軌道設(shè)計

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

若目標(biāo)軌道參數(shù)包括入軌點經(jīng)度和返回軌道經(jīng)度，可以將T0和x0與前述4個參數(shù)一起作為設(shè)計變量，構(gòu)建一個6對6的計算模型，可對T0和x0做如下調(diào)整

(20)

發(fā)射和返回的時間差ΔTi和ΔTc根據(jù)軌道與發(fā)射點和落點的經(jīng)度差計算。

3.3 數(shù)值計算討論

由于繞月自由返回軌道繞飛地球和月球，動力學(xué)環(huán)境變化復(fù)雜，數(shù)值積分時采用自適應(yīng)變步長積分器較好。

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

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

由于地心和月心附近軌道分別存在對應(yīng)升軌和降軌的P平面參數(shù)，再結(jié)合地心和月心處的順行和逆行方向，空間繞月自由返回軌道共有16種構(gòu)型。

4 數(shù)值仿真

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

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

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

表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語言編寫的相關(guān)代碼，求解一條空間雙二體模型繞月自由返回軌道僅約耗時0.2 ms，至于考慮日地月引力和地球J2項的高精度力模型，采用KSG積分器，求解一條繞月自由返回軌道僅約耗時0.5 s。在相似計算條件下，應(yīng)用文獻[15]的方法計算同樣的轉(zhuǎn)移軌道則需要至少1個小時，因此基于P平面參數(shù)的自由返回軌道設(shè)計方法具有非常高的計算效率。

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

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

圖9 轉(zhuǎn)移時間與|x0|關(guān)系曲線Fig.9 Relation curves between transfer time and |x0|

5 結(jié) 論

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