近地小行星探測任務月球借力逃逸轉移混合優化方法

井 泉，劉靖怡，李明濤

(1.中國科學院 國家空間科學中心，北京 101499； 2.中國科學院大學，北京 100049)

小行星作為太陽系形成之初遺留的原始物質，保存著太陽系形成與早期演化信息，對研究太陽系起源與早期演化、地球生命和水的來源等重大科學問題具有重要意義。對小行星的探索在一定程度上還能牽引深空探測技術的發展。

最早的小行星探測是由美國的Galileo木星探測器完成的，探測器借力飛行序列為金星-地球-地球，第一次地球借力后于1991年12月順訪了小行星Gaspre，在第二次地球借力后于1993年8月順訪了小行星Ida[1]。此后，美國的Dawn探測器是第一個探測小行星帶的探測器，采用電推進和一次火星借力，分別于2012年2月、2015年9月抵達Vesta、Ceres進行探測[2]；OSIRIS-Rex探測器采用深空機動和一次地球借力，于2018年11月抵達小行星Bennu，計劃完成采樣并返回[3]。日本Hayabusa探測器利用地球借力與小推力轉移，于2005年9月13日飛抵Itokawa進行觀測[4]；Hayabusa 2探測器采用小推力與地球借力的方式實現軌道轉移，完成了對小行星1999 JU3采樣返回[5]。

對于深空探測，通過月球借力(lunar gravity assist， LGA)逃逸來轉移，可以減小發射能量，提升送入深空的有效載荷質量[6]。月球借力最早可以追溯到美國宇航局的ISEE-3任務，1982年完成任務后飛行器利用月球借力從日地L1點轉移至L2附近探測地球磁尾，并于1983年12月22日再次利用月球借力逃逸前往與彗星交會[7]。日本航天局1990年的Hiten任務中，飛行器完成10次月球借力[8]；之后的Planet-A與Planet-B計劃中，均利用月球引力，前者用于研究考慮太陽引力影響下月球借力的效果，后者兩次月球借力后前往火星[9]。Gil-Fernandez等[10]以美歐探測器ExoMars為背景，研究了單次、二次、三次月球借力發射探測器的情況，并且增加脈沖輔助，發現ExoMars探測器在脈沖輔助的作用下，入軌質量比單次借力和兩次借力分別提升150 kg和270 kg。Casalino等[11]以探測近地小行星為背景，研究了兩次月球借力的情形。對于多次借力，一般有Backflip借力、Coplanar借力以及共振借力。Backflip借力指航天器在一次借力后轉過180°的地心角，然后進行第二次借力。Coplanar借力指第一次借力后的軌道平面與月球軌道共面[11]。共振借力指一次借力后的軌道周期與借力天體軌道周期比為p∶q(p、q均為正整數)。He等[12]研究了小推力情況下多次月球共振借力的小行星采樣返回軌道設計，推導了p∶q共振借力計算方法；楊洪偉等[13]解析推導了n∶1(n為正整數)共振借力軌道的充要條件。

本文考慮化學推進，采用圓錐曲線拼接法[14]。任務軌道分為兩段，第一段為探測器從地球停泊軌道出發，經月球借力轉移至地球影響球邊界；第二段為日心段，從影響球邊界至小行星。單次月球借力深空逃逸軌道設計方法有兩種：正向法，設計順序與航天器運行順序相同[15],缺點是由于優化變量較多，使用全局優化算法優化耗時長；反向法，首先求解最優日心轉移軌道，然后根據飛行器在地球影響球邊界處的狀態求解地心段軌道[16],缺點是依賴于日心段解，設計時沒有考慮月球借力窗口約束，并不一定是全局最優，并且由于模型采取了較多假設，解的精度較差。

本文結合正向法精度高與反向法收斂速度快的優點，提出了混合優化方法。并以交會近地小行星為例開展數值仿真，驗證了本文方法的有效性。

1 模型介紹

1.1 正向法模型

正向法的月球借力轉移軌道各時刻如下：t0時刻從地球停泊軌道出發；t0+Δt1到達月球附近，借力；t0+Δt1+Δt2時刻到達地球影響球邊界，逃逸；t0+Δt1+Δt2+αtf時刻深空機動；t0+Δt1+Δt2+tf時刻到達小行星[17]。任務時序見圖1。

圖1 正向法月球借力逃逸轉移示意Fig.1 Lunar gravity assist transfer trajectory in forward patched method

飛行器以雙曲線軌道逃逸，軌道方程為

R=|a|(ecoshF-1)

(1)

式中：R為探測器到地心的距離，a為軌道半長軸，F為偏近點角。偏心率e為

(2)

式中：Rp為逃逸時地心距；V0為逃逸速度。F與真近點角f的關系為

(3)

在計算中，取R為t0時刻地月距離，由式(1)～(3)可以求得F、f，進而求得近地點到月球的轉移時間Δt。

借力時月球的赤經λm和赤緯Bm，需滿足ω0=U-f，Ω0=λm-ΔΩ0,見圖2。

圖2 航天器軌道根數與月球交會點關系Fig.2 Relationship between orbital elements and Moon intersection

根據球面三角關系，便得到借力前的其余軌道根數。給出一組借力高度hm和角度參數ψm，即可求出借力后探測器的位置和速度，借力后至地球影響球邊界的飛行時間為Δt2。再進行日心段的遞推，采取先滑行、后機動的策略，轉移總時間tf，滑行時長占轉移時間的占比為α。求解以施加機動處為起點、到達時刻小行星位置為終點的Lambert問題，即可求出兩處所需速度脈沖Δv1、Δv2。以上過程形成如下優化問題：

優化變量X=(t0,V0,hm,ψm,tf,α)；目標函數J=Δv1+Δv2。

1.2 反向法模型

圖3 反向法月球借力逃逸轉移示意Fig.3 Lunar gravity assist transfer trajectory in backward patched method

日心段轉移模型與正向法相同。優化變量X=(t0,v∞,A,δ,tf,α)，v∞、A、δ分別為雙曲超速的模長、赤經、赤緯。t0為飛行器在地球影響球邊界處逃逸時刻，此時飛行器位置速度視作[rE,vE+v∞]，rE、vE為地球的位置和速度。求解對應的Lambert問題，可得到機動處與到達小行星時所需速度脈沖Δv1、Δv2。目標函數J=Δv1+Δv2。

1)月球至逃逸的軌道設計

為了求解這一段的軌道，定義新坐標系——月球軌道面坐標系，并且將各量無量綱化，詳細過程參見文獻[16]。

日心段優化后可得逃逸地球時的雙曲超速v∞。將v∞坐標變換后到月球軌道面坐標系，記為vesc。對于正負vesc,z，逃逸和升交點赤經在同側、異側記作iz=±1，見圖4。

圖4 逃逸雙曲超速地心示意Fig.4 Geometry of escape hyperbola

軌道的半長軸由vesc的大小確定。指向升交點的單位矢量un=(cosΩ2,sinΩ2,0)，可求出vesc與un間的夾角α。

由圖4可得到如下的角度關系

(4)

式中Φ=arccos(1/e2)。使用ia=±1分別標記在升交點處借力與降交點處，真近點角對應ν2=-ω2與ν2=-ω2+π。整理(4)式可得關于e2的方程

(5)

至此，給定初始猜測值Ω2，可依次得到軌道根數、轉移時間Δt2、月球真實位置，根據位置可求出下一步迭代Ω2的初值。經過若干次迭代后，可快速收斂到真實的Ω2，得到月球至地球影響球邊界的轉移軌道。

2)地月轉移軌道設計

這段軌道的升交點赤經Ω1與前一段的升交點赤經Ω2關系為Ω1=Ω2或Ω1=Ω2+π。假設月球沿圓軌道運行，u,v,w分別是航天器與月球相遇時的徑向、切向和法向速度，VM是月球的速度，則滿足

(6)

式中u、v、w可以用下式表示：

(7)

(8)

(9)

式中ν1是和月球相遇時轉移軌道的真近點角,rp為逃逸時地心距,整理方程(6)～(9)可得

(10)

式(10)為關于e1的一元二次方程，可直接解出e1。進而解出地月轉移軌道根數，可求得借力前的速度矢量和相對月球的雙曲超速。借力前后的雙曲超速已知，便可求出轉角，進而解出借力半徑。為了盡量充分利用月球引力并保證飛越時飛行器的安全，本文設置最低借力高度為100 km。

2 混合法模型

首先利用反向法計算效率高的特點，采用多目標優化算法，快速求解出N組分布在不同參數空間的參考軌道；然后將參考軌道參數代入反向法的優化模型，使用局部優化算法在參考軌道附近搜索；最后選取速度增量較小的M組解(M

圖5 混合法計算流程圖Fig.5 Algorithm of hybrid method

本文多目標優化算法采用NSGA-Ⅱ算法[18]，該算法主要使用Pareto前沿(Pareto front)機制：對于一個決策變量，如果不存在其他決策變量可以支配它，則稱該決策變量為非支配解。非支配解意味著無法在改進任何目標函數的同時不削弱至少一個其他的目標函數。所以每一個非支配解都是多目標優化問題的一個最優解，這些解在目標空間可視作在一個超曲面上，即Pareto前沿。并且算法中引入擁擠度的概念，使得到的解盡量均勻地分布在目標空間內。

通過多目標優化，可以得到分散的若干個解。這樣可以為下一步序列二次優化提供較好的初值。

序列二次優化將非線性問題轉化為一系列二次規劃子問題來求解[19]。利用多目標優化得到的初值在初值附近優化，優化速度快，并可以確保優化出的解都達到各自的局部最優(即從圖6所示三角形代表的解優化至五角星解)。

圖6 模型差異示意Fig.6 Illustration of difference between two models

最后將上一步中得到的解帶入正向法模型。由于正向法的模型比反向法的更精確，兩種模型下的解空間可能有所差異，見圖6。因局部優化受初值影響較大，全局搜索算法[20]會在反向法給出的初始解附近生成若干起始點，然后基于正向法模型進行局部優化(即從圖6中五角星代表的解優化至星號解)。在保證解滿足正向法精度的同時，也提高了收斂到全局最優解的概率。

混合優化方法的好處是：充分利用多目標優化解的多樣性、局部優化算法快速收斂性和全局搜索算法對初值的魯棒性，可以加強解的穩定性，并提升計算速度。

3 仿真結果

3.1 算例1

以探測1989 ML小行星為例，開展月球借力逃逸轉移軌道優化設計。1989 ML的軌道根數見表1。

表1 1989 ML歷元及軌道根數(日心黃道坐標系)Tab.1 Epoch and orbital elements of 1989 ML (referenced in the heliocentric ecliptic coordinate frame)

正向法設置發射窗口為2024～2027年，轉移時間100～600 d；反向法設置發射窗口為2025年5月～2025年11月，轉移時間200～600 d。

設置遺傳算法種群數1 000，最大代數500。仿真計算機CPU為i7-7700，主頻3.60 GHz，四核并行。正向法、反向法、混合法的仿真結果見表2。

表2 1989 ML各方法仿真結果Tab.2 Results of 1989 ML obtained by different methods

從表2可以看出：直接利用遺傳算法開展正向法優化，解的穩定性較差，10次仿真只有2次收斂到最優解；并且單次優化用時688 s，效率在幾種方法中最低。直接利用遺傳算法開展反向法優化，雖然單次優化用時只有100 s，但從1.480 km·s-1的速度增量和其他量的值，可以看出并沒有優化至最優解；并且10次優化只有3次收斂到該解，穩定性較差。混合法中，多目標優化兩個優化目標分別是速度增量與發射C3；設置保留6組多目標優化個體(N=6)。反向法多目標優化時Pareto前沿見圖7，當前前沿包含的解對應的軌道速度增量分布在1.5～5.5 km·s-1，C3分布在0～4.5 km2·s-2，可見，通過多目標優化能夠獲得具備多樣性的解。局部優化后選擇速度增量小的前3組解(M=3)，再使用全局搜索算法優化。仿真結果顯示，混合法對局部優化后的3組解進行了有效優化，10次仿真中有8次收斂到最優解，混合優化法解的穩定性相比正向法提升4倍。單次優化用時相比正向法節省355 s，用時縮短50%以上。

圖7 1989 ML速度增量-C3 ParetoFig.7 1989 ML Pareto front of Δv and C3

交會小行星1989 ML混合法優化結果對應的軌道見圖8、9。

圖8 1989 ML月球借力轉移示意Fig.8 Illustration of LGA transfer to 1989 ML

圖9 1989 ML日心轉移軌道示意Fig.9 Illustration of heliocentric transfer to 1989 ML

3.2 算例2

以探測2003 SM84小行星為例，開展月球借力逃逸轉移軌道優化設計。2003 SM84的軌道根數見表3。

表3 2003 SM84歷元及軌道根數(日心黃道坐標系)Tab.3 Epoch and orbital elements of 2003 SM84 (referenced in the heliocentric ecliptic coordinate frame)

正向法設置發射窗口為2027～2030年,轉移時間100～600 d;反向法設置發射窗口為2027年5月～2027年11月，轉移時間200～600 d。

仿真結果見表4。從表4可以看出：利用遺傳算法開展正向法優化，10次仿真只有2次收斂到最優解，解的穩定性較差；平均單次優化用時623 s，效率在幾種方法中最低。利用遺傳算法開展反向法優化，10次仿真有4次收斂到最優解，雖然單次優化用時少，但由于反向法的模型簡單，得出的軌道并不是最優解。混合法中，多目標優化時Pareto前沿見圖10。當前前沿包含的解對應的軌道速度增量分布在1.5～3 km·s-1，C3分布在0～5 km2·s-2，可見，通過多目標優化獲得了具備多樣性的解。仿真結果顯示，10次仿真中有6次能收斂到最優解，穩定性相比正向法提升3倍。單次優化用時相比正向法節省312 s，用時縮短約50%。有效地提升了解的穩定性和計算效率。

表4 2003 SM84各方法仿真結果Tab.4 Results of 2003 SM84 obtained by different methods

圖10 2003 SM84速度增量-C3 ParetoFig.10 2003 SM84 Pareto front of Δv and C3

交會小行星2003 SM84混合法優化結果對應的軌道見圖11、12。

圖11 2003 SM84月球借力轉移示意Fig.11 Illustration of LGA transfer to 2003 SM84

圖12 2003 SM84日心轉移軌道示意Fig.12 Illustration of heliocentric transfer to 2003 SM84

4 結 論

針對正向法全局優化效率低、收斂穩定性差、反向法模型精度不高的缺點，為進一步提升算法收斂穩定性和收斂速度，本文提出了結合兩種方法優點的混合優化方法：以多目標優化反向法快速求解出的Pareto解集為初始參考軌道，局部優化后再利用正向法精確模型下在參考軌道附近開展全局搜索，獲得最優轉移軌道。最后以1989 ML和2003 SM84小行星為探測目標進行仿真，結果顯示解的收斂穩定性相比正向法提升3倍以上，計算用時縮短50%以上。驗證了該方法的有效性。