基于彈道逃逸和小推力捕獲的地月轉移軌道設計

彭坤，黃震，楊宏，張柏楠

中國空間技術研究院 載人航天總體部，北京 100094

隨著美俄聯合開展“深空之門”地月空間設施方案設計，開發地月空間已成為世界各航天機構的共識，同時往返地月空間也有助于獲取月球以遠載人深空探測所需的技術和經驗[1-2]。而地月軌道轉移方式及其轉移軌道設計是開發地月空間的基礎。目前已實施的月球探測工程主要有兩種地月軌道轉移方式：雙脈沖地月軌道轉移和全小推力地月軌道轉移。美國的阿波羅載人登月工程[3]和中國的嫦娥月球探測工程[4]都采用了前者，即在近地軌道施加脈沖進入地月轉移軌道以實現彈道逃逸，在近月點施加脈沖進入環月軌道以實現彈道捕獲。其優點是軌道轉移時間短，僅為3～5天，但由于采用化學推進方式，比沖低，推進劑消耗大。歐空局的SMAT-1衛星則采用全小推力變軌方式實施地球逃逸和月球捕獲[5]。其優點是小推力發動機比沖高，節省推進劑，但由于推力小，飛行時間長約兩百天，同時在地球逃逸過程中需要反復穿越范艾倫輻射帶(高度范圍為內帶1 500～5 000 km，外帶13 000～20 000 km)，對航天器的電子設備是嚴峻考驗。為此，本文提出一種彈道逃逸和小推力捕獲相結合的地月軌道轉移方式，先利用運載火箭在近地軌道上脈沖加速將航天器送入地月轉移軌道；進入地月空間后航天器利用自身的高比沖小推力發動機擇機開機進行月球捕獲，以取代近月點的脈沖制動，從而減少航天器燃料消耗。

基于彈道逃逸和小推力捕獲的地月轉移軌道設計主要涉及脈沖地月轉移軌道設計、小推力地月轉移軌道設計以及兩段軌道拼接點選取等問題。脈沖地月轉移軌道不同于單中心引力體脈沖變軌軌道，其非線性強，一般先在低保真模型中計算初值，再在高保真模型中修正初值以得到精確值。目前常用的低保真動力學模型有雙二體模型[6]、三體模型[7,8]，以及偽狀態模型[9,10]。也有學者直接在高保真模型中采用分層搜索流程[11,12]或基于高精度初值估計及改進修正方法[13]直接求解脈沖地月轉移軌道。小推力轉移軌道設計方法分為三類：間接法、直接法和混合法。Ranieri和Ocampo[14]基于間接法思想，用曲線擬合技術估計地球逃逸段和月球捕獲段的伴隨變量初值，并將月心伴隨變量轉換到地心系進行求解。Lee和Bang[15]則基于間接法理論采用最優螺旋曲線初值估計方法進行伴隨變量初值估計并通過多層搜索逐步確定地月轉移軌道。Betts和Erb[5]采用直接法將地月轉移軌道離散化，再利用序列二次規劃算法進行求解。由于地月小推力轉移軌道飛行時間長，為保證軌道精度，離散化的狀態變量和控制變量多達20萬個，對于一般計算機而言屬于海量計算。Pierson和Kluever[16]、Gao[17]等采用混合法思想，將地心段和月心段的伴隨變量初值作為優化變量，采用非線性規劃算法進行優化求解。徐明和徐世杰[18]則基于地月系統L1點Halo軌道構造了小推力地月轉移自由飛行段，并設計合適控制律完成地球加速段和月球減速段軌道設計。對于基于彈道逃逸和小推力捕獲的地月轉移軌道，逃逸段和捕獲段軌道的拼接也是其設計難點，若小推力發動機開機過早會過多消耗推進劑；若小推力發動機開機過晚則有可能不足以抵消脈沖地月轉移軌道的對月速度，從而導致航天器月球捕獲失敗而逃離月球。目前還未有針對基于彈道逃逸和小推力捕獲的地月轉移軌道的綜合設計方法研究。

考慮地月轉移軌道主要受地球和月球引力影響[19]，本文在三體模型下建立地月轉移軌道的動力學模型。對于彈道逃逸段，建立二維地心旋轉系模型，采用地心旋轉系對準角和地月轉移加速速度增量作為控制變量，給出其初值估計的解析式，應用序列二次規劃算法求解出滿足目標約束的脈沖地月轉移軌道。對于小推力捕獲段，建立二維月心旋轉系模型，進行逆向軌道仿真以提高軌道設計收斂性；以近月點伴隨變量為優化變量，以月心距為基準進行逃逸段和捕獲段軌道拼接以及目標約束設置；采用能量匹配觀點預估轉移時間，同時應用最優螺旋曲線伴隨變量解析式預估初值范圍；最后利用引導型人工免疫算法(Guiding Artificial Immune Algorithm，GAIA)對小推力月球捕獲軌道進行優化求解，從而得到整條基于彈道逃逸和小推力捕獲的地月轉移軌道。

1 系統模型

基于彈道逃逸和小推力捕獲的地月轉移軌道飛行過程為：航天器在近地圓軌道(Low Earth Oribt，LEO)上運行，在A點施加切向脈沖Δv1，進入脈沖逃逸地月轉移軌道(Ballistic Escape Trajectory, BET)，并在BET上無動力飛行；設到達B點時航天器小推力發動機開機進行減速制動，進入小推力月球捕獲軌道(Low-Thrust Capture Trajectory, LTCT)，直至航天器到達近月點C進入目標環月軌道(CircumLunar Orbit，CLO)，其飛行軌跡示意圖如圖1所示。

忽略影響較小的攝動力，在圓型限制性三體模型下建立地心慣性系下航天器的動力學方程，如式(1)所示。其中月球繞地球以ω做勻速圓周運動，航天器只受到地球和月球引力作用。

(1)

式中：μE和μM分別為地球和月球的引力常數；Re為航天器相對地心的位置矢量，其標量長度為Re；Rm為航天器相對月心的位置矢量，其標量長度為Rm；RME為月球相對地心的位置矢量，其標量長度為RME。

設航天器在A點出發時刻為t0=0，以t0時刻地月連線為X軸，由地心指向月心，Z軸為月球軌道角動量方向，建立地心旋轉系。將式(1)從地心慣性系轉換到地心旋轉系中，并進行極坐標變換x=rcosθ，y=rsinθ可得

圖1 基于彈道逃逸和小推力捕獲的 地月轉移軌道軌跡Fig.1 Trans-lunar trajectory based on ballistic escape and low-thrust capture

(2)

對于小推力月球捕獲軌道段，可建立月心旋轉系下航天器極坐標運動方程，并加入小推力發動機攝動加速度，如式(3)所示。此時，地球作為次天體，其月心旋轉系坐標為(-d,0)。

(3)

圖2 二維地月極坐標系統Fig.2 Two-dimensional Earth-Moon polar coordinate system

2 彈道逃逸地月轉移軌道設計

圖3 地心旋轉系下彈道逃逸軌道與 地月的幾何關系Fig.3 Geometrical relationship between BET, Earth and Moon in geocentric rotating frame

(4)

為減少搜索過程中的軌道仿真時間，設置月心旋轉系航跡角γ2的余弦值Γγ2滿足式(5)的時刻tf為軌道積分終止條件，以保證航天器終端時刻處于近月點。

(5)

(6)

為滿足近月距和轉移時間約束，設置目標約束條件為

(7)

為增加軌道搜索的收斂性，需提供精確的設計變量初值。雙脈沖地月轉移軌道大部分軌道段主要受地球引力影響，本文將其近似為二體模型下的霍曼轉移軌道，此時月球可看作一個普通航天器，忽略其引力影響。通過霍曼轉移公式可給出計算Δv1的解析式；同時由霍曼轉移軌道拱線應對準目標航天器交會位置可知，α等于轉移時間ttransfer內月球公轉轉過的角度(如圖3所示)。根據以上分析可給出設計變量初值估計解析式：

(8)

根據設計變量的初值估計值，本文采用常用的序列二次規劃算法SNOPT[20]進行迭代求解。以下給出he=200 km，hm=2 500 km，ttransfer=5 d的地心順行月心逆行的彈道逃逸地月轉移軌道設計結果，如表1和圖4所示。從表1可看出，對于5天地月轉移軌道，初值估計值非常接近精確值；同時本文設計方法搜索速度快，迭代時間不到4 s，迭代4次就收斂到精確地月轉移軌道。圖4給出了地心旋轉系下彈道逃逸地月轉移軌道的軌跡。從圖4可看出，該軌道大部分軌道段位于月球影響球外，主要受地球引力影響，驗證了設計變量初值估計二體模型假設的合理性。

表1 彈道逃逸地月轉移軌道搜索結果

圖4 地心旋轉系下彈道逃逸地月轉移 軌道軌跡Fig.4 Ballistic escape trans-lunar trajectory in geocentric rotating frame

圖5 月心旋轉系下彈道逃逸地月轉移軌道狀態變量Fig.5 State variables of ballistic escape trans-lunar trajectory in selenocentric rotating frame

3 小推力月球捕獲軌道設計

不同于脈沖轉移軌道，小推力轉移軌道是一條漸變螺旋軌道，需要設計轉移過程中每一時刻的推力矢量。對于小推力月球捕獲這類多圈軌道，轉移時間長且月心距變化范圍大。若采用間接法，由于多圈軌道的傳遞性，進一步增大了伴隨變量初值的敏感性，初值猜測困難[17]。若采用直接法將軌道進行離散，為保證軌道設計精度，涉及海量設計變量搜索，計算過程復雜[5]。因此，本文基于混合法思想，采用GAIA[21]對伴隨變量進行搜索求解。同時，小推力月球捕獲軌道起始點B處于地月之間引力敏感區域(如圖1所示)，動力學非線性強，而終端點C處于月球引力場中，軌道特性穩定。若直接從B點軌道順推設計小推力軌道，初值敏感性進一步加強。為增加軌道搜索的收斂性，本文采用軌道逆推的方法設計小推力月球捕獲軌道，以終端點C為軌道仿真起始點，而將起始點B作為軌道仿真終端點。

(9)

根據極大值原理，小推力最優控制角為

(10)

設C點時刻為tC=0，考慮目標環月軌道為逆行軌道，則軌道逆推的小推力月球捕獲軌道初始條件為

(11)

式中：由于近月點C的極角θ2(0)自由，相應的伴隨變量λθ2(0)=0。

(12)

一般拼接點的月心距較大，較小的θ2誤差可導致較大的位置誤差，因此θ2誤差無法真實表征位置誤差。這里將θ2約束轉換為相對距離誤差，則終端約束可進一步變為

(13)

小推力月球捕獲軌道要求燃料消耗最小，同時考慮終端約束，將其作為罰函數引入評價函數。由于GAIA搜索極大值，故將評價函數設為如下形式，作為GAIA的親和度值[21]。

{r2(-tLT)·

σ2|vr2(-tLT)-fvr2(r2(-tLT))|-

σ3|vθ2(-tLT)-fvθ2(r2(-tLT))|

(14)

式中：σ1、σ2和σ3為罰函數系數，采用模擬退火算法自動調整系數取值，σi=max{1,10/β(k/kmax-1)},i=1,2,3；k和kmax分別為當前尋優代數和最大尋優代數。

(15)

根據文獻[15]的研究結論，τ可設為-0.012，κ為伴隨變量等比例縮放系數。本文中κ的取值原則是保證3個伴隨變量初值均大于1，使得3個伴隨變量初值的搜索間隔足夠大，從而使GAIA能進行充分的全局搜索。對于小推力捕獲軌道的近月點C無角度要求，故θ2(0)∈[0°,360°]。

圖6 彈道逃逸地月轉移軌道月心能量Fig.6 Energy of ballistic escape trans-lunar trajectory in selenocentric frame

模型下小推力捕獲軌道的能量和極徑變化曲線，直至與彈道逃逸軌道的能量和極徑完全匹配，此時的轉移時間tLTguess可作為初值。

4 仿真算例及分析

設LEO為he=200 km的順行軌道，CLO為hm=2 500 km的逆行軌道。設航天器到達近月點質量為m(0)=1 000 kg。考慮發動機推質比為4×10-3m/s2[15]，小推力發動機參數設為T=4 N。參考文獻[17]，比沖取為w=3 500×9.8 m/s。彈道逃逸軌道選擇ttransfer=5 d的脈沖地月轉移軌道，其搜索結果見表1和圖4。

取κ=50 000，根據初值估計公式可計算出λr2guess=12.55 723，λvr2guess=600，λvθ2guess=50 000。考慮一定搜索區域和各伴隨變量的敏感性，設置伴隨變量初值的搜索范圍為λr2(0)∈[0.6λr2guess,1.4λr2guess]，λvr2(0)∈[-10λvr2guess,10λvr2guess]，λvθ2(0)=[1.2λvθ2guess,0.8λvθ2guess]。采用切向推力模式計算出脈沖逃逸軌道和小推力捕獲軌道的能量-月心距曲線及其匹配點，如圖7所示。匹配結果為tLTguess=4.425 d，此時r2=121 832.801 km，ε2=0.277 km2/s2。設左右偏差1.5 d，則小推力捕獲軌道轉移時間搜索范圍為tLT∈[3,6] d。近月點C處極角搜索范圍為θ2(0)∈[0°,360°] 。

根據以上搜索范圍采用GAIA進行尋優搜索，種群個數為N=100，最大尋優代數為kmax=60，搜索到的彈道逃逸和小推力捕獲的地月轉移軌道結果如表2和圖8～圖10所示。由表2可知，小推力捕獲軌道終端位置速度誤差非常小，滿足與彈道逃逸軌道的拼接點終端約束。小推力捕獲軌道燃料消耗為39.670 kg，轉移時間為3.940 d。

由圖8可知，根據初值預估的搜索范圍，GAIA在第10代就收斂到最優值附近，搜索速度快。圖9繪制了小推力捕獲軌道的軌跡，由圖可知航天器經過約4圈軌道進入環月軌道。圖10給出了小推力發動機的控制角變化曲線，由圖可知越靠近月球，控制角振蕩幅值越小，近似于切向制動。

圖7 彈道逃逸軌道和小推力捕獲軌道的月心能量Fig.7 Energies of ballistic escape trajectory and low- thrust capture trajectory in selenocentric frame

表2 小推力月球捕獲軌道搜索結果Table 2 Search results of low-thrust capture trajectory

圖8 最大親和度隨迭代代數變化曲線Fig.8 Maximum affinity value of every generation

圖9 月心旋轉系下小推力捕獲軌道軌跡Fig.9 Low-thrust capture trajectory in selenocentric rotating frame

圖10 小推力捕獲軌道的控制角Fig.10 Steering angle of low-thrust capture trajectory

圖11繪制了整條彈道逃逸和小推力捕獲地月轉移軌道在地心旋轉系和地心慣性系下的軌跡，其中虛線軌跡為不進行小推力捕獲而沿著原脈沖逃逸軌道自由飛行的虛擬彈道捕獲軌道(Virtual Ballistic Capture Trajectory，VBCT)。航天器在彈道逃逸軌道上飛行3.617 d后，小推力發動機開機進入小推力捕獲軌道；在小推力捕獲軌道上飛行3.940 d后到達目標環月軌道。其中，航天器在彈道逃逸軌道上飛行約1.2 h后地心高度超過20 000 km，即穿過范艾倫輻射帶。

圖11 基于彈道逃逸和小推力捕獲的 整條地月轉移軌道Fig.11 Whole trans-lunar trajectory based on ballistic escape and low-thrust capture

圖12給出了整條地月轉移軌道月心旋轉系下的狀態變量變化曲線，并與虛擬彈道捕獲軌道進行對比。其中，兩者的終端月心距和月心徑向速度相同，而極角和月心橫向速度不同。兩者月心距相同說明均滿足目標環月軌道高度約束，月心徑向速度相同且為0說明了終端點均為近月點。極角不同說明兩者的目標環月軌道入軌點不同，月心橫向速度不同是因為虛擬彈道捕獲軌道還未實施制動機動。若采用脈沖制動，需要速度增量647.507 m/s，設化學推進發動機比沖為310×9.8 m/s，則燃料消耗為199.568 kg，是小推力捕獲軌道燃料消耗的5倍，說明了采用小推力捕獲可大大節省燃料消耗。同時，采用小推力捕獲方式總飛行時間為7.557 d，僅比5 d雙脈沖地月轉移軌道多2.557 d。因此，相對于全小推力地月轉移軌道200 d左右的轉移時間，以及雙脈沖地月轉移軌道3～5 d的轉移時間，基于彈道逃逸和小推力捕獲的地月轉移軌道的轉移時間適中。

圖12 月心旋轉系下整條地月轉移軌道狀態變量Fig.12 State variables of whole trans-lunar trajectory in selenocentric rotating frame

5 結 論

1)對于運載火箭發射進入地月轉移軌道的電推進月球探測器以及月球軌道空間站任務，提出了一種基于彈道逃逸和小推力捕獲的新型地月軌道轉移方式，并建立了二維極坐標系統模型。

2)將對準角和速度增量作為彈道逃逸軌道的設計變量，提出了初值預估解析式，仿真結果表明初值估計接近真實值，求解模型易于收斂。

3)利用月心距解決了小推力捕獲軌道和彈道逃逸軌道的拼接問題。采用能量匹配和最優螺旋軌道伴隨變量關系猜測初值范圍，并以軌道逆推方式提高收斂性，實現了最優小推力捕獲軌道的快速搜索。

4)與雙脈沖地月轉移軌道相比，基于彈道逃逸和小推力捕獲的地月轉移軌道近月制動燃料消耗約減少了80%，提高了航天器有效載荷比重或使航天器保留更多推進劑以延長其環月飛行壽命，同時其飛行時間僅增加幾天時間。

5)與全小推力地月轉移軌道相比，基于彈道逃逸和小推力捕獲的地月轉移軌道飛行時間降低了一個數量級，同時其利用彈道逃逸軌道不到2小時即可快速穿越地球輻射帶。