王鵬，武小宇，張立華

（1.航天東方紅衛(wèi)星有限公司，北京 100094；2.北京理工大學 宇航學院，北京 100081；3.飛行器動力學與控制教育部重點實驗室，北京 100081）

引 言

通過對小天體探測可以了解太陽系演化、生命起源等基本規(guī)律，并有可能在未來實現(xiàn)小天體資源利用。隨著航天技術(shù)的發(fā)展，小天體探測逐漸成為人類深空探測任務的熱點目標之一[1-2]。1997年, 第一個專門探測小行星的太空計劃NEAR[3]成功交會了433 Eros 小行星[3-4]，并在途中飛越了253 Mathilde小行星?！吧羁?號”探測器于1999年造訪了9969 Braille 小行星[5]。2005年，日本“隼鳥號”完成了“絲川”小行星的采樣任務并返回地球，使人類第一次從小行星上獲取了樣本[6]。美國“黎明號”探測器于2007年飛往 Vesta和 Ceres 小行星并進行探測[7]。美國國家航空航天局（NASA）宣布了近地小行星采樣返回任務，計劃于2022年發(fā)射探測器，2028年到達101955 1999RQ36 小行星采集土壤并于2032年送回地球[8]。

迄今為止，通過天文觀測和探測任務，人類已經(jīng)基本確定了太陽系小天體的化學組成和內(nèi)部結(jié)構(gòu)，小天體中蘊藏有豐富資源已成為共識[7]。2015年美國通過的《2015外空資源探索和利用法》標志著小天體商業(yè)開發(fā)熱潮的開始。小天體抵近資源勘查是滿足小天體資源商業(yè)開發(fā)風險控制的可行途徑，但傳統(tǒng)抵近探測模式又很難滿足商業(yè)開發(fā)對成本的要求。以歐洲最為成功的小天體抵近任務ROSETTA為例，花費了7億歐元，10年時間，利用甩擺軌道技術(shù)[9]實現(xiàn)了一顆彗星的附著和2顆小天體的飛越[10]。其采用的一顆探測器串行探測多顆小天體的方式，使得多個小天體探測任務緊密耦合，任務復雜度高、風險大，并不完全適合商業(yè)勘查的要求。

隨著微納飛行器相關(guān)技術(shù)的不斷進步，50 kg級小天體交會探測任務正成為現(xiàn)實。微納飛行器能夠在相對短的周期（12～18個月）實現(xiàn)從任務概念設(shè)計到任務實施，其尺寸、重量和功耗都相對較小，成本低廉。如果多個微納小天體探測器能夠以一箭多星的方式發(fā)射至日地拉格朗日點進行停泊與軌道保持，而后對多個小天體進行一對一交會探測，則能夠有效縮短探測任務的時間并降低經(jīng)濟成本、提高探測效能。利用一顆CZ3A或CZ4C火箭可以實現(xiàn)6～10顆50 kg級的微納飛行器小天體交會勘查任務，經(jīng)濟成本有可能降到30萬人民幣/千克，任務全壽命周期有可能縮短到5年以內(nèi)，是實現(xiàn)“一探多”的一種可行途徑。

以飛行時間為約束，降低微納飛行器的星際轉(zhuǎn)移成本是多小天體并行探測的關(guān)鍵技術(shù)之一。本文在圓型限制性三體模型下構(gòu)建了結(jié)合行星借力與不變流形機制的低能量轉(zhuǎn)移方案。引入了擾動流形思想，實現(xiàn)與目標小天體的快速交會。進一步，提出了一種基于多小天體探測的全局搜索方法，該方法基于在日地軌道上停泊的微納飛行器集群，在線確定探測目標，逐次離軌進行探測，從而達到發(fā)射一次、探測多個目標的“并行”探測方式。本文數(shù)值仿真采用不同類型的近地小天體（2009BD，2008EA9，99942 Apophis）作為示例進行數(shù)值仿真分析，結(jié)果表明，該方案能夠大幅度降低轉(zhuǎn)移任務的燃料消耗，并縮短轉(zhuǎn)移時間。

1 動力學原理

1.1 動力學模型

在近地小天體探測過程中，微納飛行器主要受到太陽、出發(fā)行星兩個天體的引力作用，為便于分析，本文采用圓型限制性三體模型對探測器的軌道動力學進行建模，忽略探測器對太陽和行星運動的影響，并且假設(shè)太陽和行星繞兩者的共同質(zhì)心作圓周運動，具體如圖1所示。

圖1 圓型限制性三體問題Fig.1 The sketch of circular restricted three-body problem

定義旋轉(zhuǎn)坐標系oxyz：原點為太陽與地球的共同質(zhì)心，x軸為原點指向地球方向，y軸在太陽—行星運動平面內(nèi)指向地球運動方向，z軸構(gòu)成右手坐標系。質(zhì)量比 μ為行星質(zhì)量與M?之 比。則在旋轉(zhuǎn)坐標系oxyz中，探測器的動力學方程可以描述為

其中：? 為探測器的偽勢能函數(shù)，具有如下表達式

在圓型限制性三體模型中，探測器的運動只存在一個獨立積分，即雅可比積分

圓型限制性三體模型中存在5個動平衡點，如圖1所示。其中 L1、L2和L3位于太陽與行星的連線上，稱為共線平衡點。L4和L5分別與太陽和行星構(gòu)成等邊三角形，稱為三角平衡點。平衡點及其附近周期和擬周期軌道存在與之相連的穩(wěn)定和不穩(wěn)定流形。不變流形的漸近特性為設(shè)計低能量轉(zhuǎn)移軌道提供了可能，探測器可以沿著穩(wěn)定或不穩(wěn)定流形逐漸地趨近或者脫離（擬）周期軌道。

1.2 三體問題下的流形結(jié)構(gòu)

為獲得halo軌道的穩(wěn)定和不穩(wěn)定流形，首先在halo軌道上選取一參考點X0=X(t0)（本文選取halo軌道正向穿越xy平面的點為參考點），計算halo軌道對應的單值矩陣M(τ=1,t0)，τ∈[0,1]為halo軌道周期歸一化后的時間參數(shù)，則τ對應的時間為

求取單值矩陣M(τ=1,t0)的特征值及對應的特征向量，得到穩(wěn)定特征向量ηS(X0)和不穩(wěn)定特征向量ηU(X0)。在[0,1]區(qū) 間內(nèi)對τ進行等分離散化，則各離散點對應的穩(wěn)定和不穩(wěn)定特征向量可以用下式計算得到

其中：M(τ,t0)為 τ時刻相對參考點的狀態(tài)轉(zhuǎn)移矩陣。

進一步，采用攝動法可以求得各離散點對應的穩(wěn)定和不穩(wěn)定不變流形初值如下

其中：ε 為微小的攝動，一般選取ε 使位置攝動處于20～50 km范圍內(nèi)。

以式（6）為初值，分別通過逆向和正向遞推軌道可以得到halo軌道的穩(wěn)定和不穩(wěn)定流形。圖2給出了太陽—地球系統(tǒng)中 L1和 L2點附近的halo軌道及與其相連的不變流形。由圖中可以看出，halo軌道周圍存在4個不同的不變流形分支。在局部區(qū)域內(nèi)，halo軌道的左端和右端各存在兩個分支，一支為穩(wěn)定流形，另外一支為不穩(wěn)定流形。

1.3 小天體探測交會模型

微納飛行器從地球停泊軌道出發(fā)轉(zhuǎn)移至目標小天體的過程（如圖3所示）可以描述為：地球停泊微納飛行器集群中的一飛行器施加機動 ?vD到達日地/地月halo軌道，進而繼續(xù)施加速度機動 ?vM到達停泊halo軌道（只針對地月halo軌道出發(fā)情況）。當探測器進入目標影響范圍時，在目標小天體附近施加機動?vA使探測器進入目標附近，實現(xiàn)目標觀測軌道的入軌。從微納飛行器的飛行歷程來看，整條轉(zhuǎn)移軌道可以分為halo軌道出發(fā)、出發(fā)行星逃逸、星際轉(zhuǎn)移、目標行星捕獲和觀測軌道入軌5個階段。

圖2 太陽–地球系統(tǒng)halo軌道不變流形結(jié)構(gòu)Fig.2 The invariant manifolds of halo orbits in the Sun-Earth system

圖3 微納飛行器從地球停泊軌道出發(fā)轉(zhuǎn)移至目標小天體的過程Fig.3 The process of a cubesat transfering from the Earth parking orbit to the target asteroid

2 多小天體并行探測關(guān)鍵技術(shù)

2.1 擾動流形技術(shù)

設(shè)計小天體行星際探測任務轉(zhuǎn)移軌道時，燃料消耗和轉(zhuǎn)移時間是通常考慮的兩個標準。合理利用halo軌道的不變流形可以降低轉(zhuǎn)移過程中的燃料消耗，同時不變流形引入也必將增加探測器的轉(zhuǎn)移時間。例如探測月球的低能量轉(zhuǎn)移軌道，其本質(zhì)是利用太陽的長時間引力攝動增大軌道的近地點直至月球軌道。影響探測器在不變流形上飛行時間的是halo軌道的穩(wěn)定系數(shù)，若減小穩(wěn)定系數(shù)則可使探測器能夠較快地脫離halo軌道。計算傳統(tǒng)不變流形時本質(zhì)上就是通過施加微小攝動改變其穩(wěn)定系數(shù)。如果適當?shù)卦龃髷z動，進一步減小其穩(wěn)定系數(shù)，減少流形在halo軌道附近的運動時間，則既可以利用不變流形軌道的低能量特性，又能夠使探測器快速脫離或抵達halo軌道。為此，本文引入擾動流形進行探測器轉(zhuǎn)移軌道的設(shè)計。本文定義的擾動流形具有以下3個特點：①僅對速度矢量施加攝動；②速度攝動方向仍沿單值矩陣穩(wěn)定和不穩(wěn)定特征矢量方向；③速度攝動范圍為1 .0×10?6～3.0×10?1km/s。擾動流形的速度初值可以寫成

其中：?V為 施加的速度攝動，χS和 χU分別為穩(wěn)定和不穩(wěn)定特征向量的速度分量。

圖4給出了太陽—地球系統(tǒng) L1點 附近z方向幅值為1.6×105km北向halo軌道的擾動流形結(jié)構(gòu)，速度攝動大小為 ?V=100m/s。 由圖中可以看出，L1點附近halo軌道的擾動流形結(jié)構(gòu)仍然由4個分支構(gòu)成，并且擾動流形與傳統(tǒng)不變流形的演化趨勢基本一致。

圖4 日地系統(tǒng)L1點Halo軌道擾動流形結(jié)構(gòu)Fig.4 The perturbation invariant manifolds of the Halo orbits near the L1 point in the Sun-Earth system

與傳統(tǒng)不變流形不同的是，擾動流形是在較大的速度攝動下產(chǎn)生的，其可以更快地脫離halo軌道區(qū)域，從而節(jié)省探測器的飛行時間。為考察擾動流形對飛行時間的影響，設(shè)定如下龐加萊截面x=1?μ，計算L1點附近halo軌道右側(cè)不穩(wěn)定流形第一次到達該截面的時間，如圖5所示，圖中橫軸為一個halo軌道周期內(nèi)halo軌道上不同流形分支對應的出發(fā)時間，縱軸為流形從halo軌道到龐加萊截面的飛行時間。

由圖5可以看到，對于 L1點 附近z方向幅值為1.6×105km北向halo軌道，其右端的傳統(tǒng)不穩(wěn)定流形到達龐加萊截面的飛行時間為200天左右，而擾動流形到達龐加萊截面的時間則大大減少，例如?V=50m/s擾動流形最小飛行時間為73.79天， ?V=300m/s擾動流形最小飛行時間僅為37.98天。由此可見，擾動流形能夠以較小的燃料消耗代價大大縮短探測器在流形上的飛行時間。

圖5 擾動流形對飛行時間的影響Fig.5 The influence of the perturbation invariant manifolds on the fly time

2.2 停泊halo軌道期間小天體發(fā)射機會搜索技術(shù)

對給定的出發(fā)和目標halo軌道，利用行星借力和擾動流形的低能量轉(zhuǎn)移軌道關(guān)鍵參數(shù)包括

其中：t0為探測器從初始halo軌道出發(fā)的時間；τD表征了出發(fā)時探測器在halo軌道上的位置；?VD為出發(fā)不穩(wěn)定擾動流形的速度攝動大??；?vD為在出發(fā)行星近拱點處施加的軌道機動脈沖；tS為由出發(fā)行星近拱點到目標行星近拱點的飛行時間；?vA為在目標行星近拱點處施加的軌道機動脈沖； ?VA為目標穩(wěn)定擾動流形的速度攝動大?。?τA表征了探測器進入目標halo軌道時在halo軌道上的位置。

性能指標可綜合考慮燃料消耗和轉(zhuǎn)移時間，選取為

其中：tD為出發(fā)halo軌道到行星近拱點的飛行時間；tD為目標行星近拱點到halo軌道的飛行時間；β>0為罰因子。

轉(zhuǎn)移軌道必須滿足出發(fā)和目標halo軌道的軌線約束。

其中：ρ為探測器的軌道狀態(tài)；XD和XA分別為出發(fā)和目標halo軌道上的狀態(tài)；tf=t0+tD+tS+tA。

另外，由于需要在行星近拱點施加軌道機動，轉(zhuǎn)移軌道還應該滿足如下內(nèi)點約束

其中：ξD和ξA分別為探測器相對于出發(fā)和目標行星的位置矢量；tPD=t0+tD；tPA=tPD+tS。

上述模型構(gòu)成了一個復雜的多點邊值問題，理論上可采用非線性規(guī)劃等技術(shù)直接求解。然而，由于此類型轉(zhuǎn)移軌道途經(jīng)多個引力場，動力學非線性強，軌道約束對設(shè)計參數(shù)異常敏感，導致數(shù)值優(yōu)化算法很難收斂，提供合理的設(shè)計初值是解決該問題的關(guān)鍵。針對該問題，本文提出一種基于模型拼接的初始轉(zhuǎn)移軌道構(gòu)建方法。由前面的討論可知，整條轉(zhuǎn)移軌道由5個軌道段構(gòu)成，5個分段的軌道特性各不相同。綜合探測器考慮飛行時間和受力特點，可分別在不同動力學模型下對各分段軌道進行討論。對于halo軌道出發(fā)和入軌段，在圓型限制性三體模型中討論。對于出發(fā)行星逃逸和目標行星捕獲段，在以行星為中心的二體模型中討論，且忽略其飛行時間。對于日心轉(zhuǎn)移段，則只考慮太陽的引力影響。需要指出的是，以上討論均考慮星歷約束下的行星相位。

根據(jù)以上假設(shè)，本文提出的初始轉(zhuǎn)移軌道搜索方法可以描述為

1）在日心二體模型中，對于給定的出發(fā)行星逃逸時間段和日心轉(zhuǎn)移時間段，利用二體Lambert算法對日心轉(zhuǎn)移軌道進行遍歷搜索，求得各遍歷點對應的探測器逃逸出發(fā)行星時的時間tD、逃逸雙曲線超速、日心轉(zhuǎn)移段飛行時間tS以及到達目標行星時的雙曲線超速。

2）在圓型限制性三體模型中，分別計算出發(fā)halo軌道不穩(wěn)定流形（傳統(tǒng)不變流形）和目標halo軌道穩(wěn)定流形對應的第一次行星近拱點龐加萊映射，并將其轉(zhuǎn)化到以行星為原點的慣性坐標系中。

3）在以行星為中心的二體模型中，分別完成不變流形近拱點軌道狀態(tài)與步驟1）各遍歷點對應和的匹配，計算得到在行星近拱點需要施加的速度脈沖?vD和?vA，確定對應最小的日心轉(zhuǎn)移軌道參數(shù)tD、、tS和。

4）在圓型限制性三體模型下，分別計算出發(fā)和目標halo軌道對應的不同速度攝動?V和不同出發(fā)位置τ擾動流形的第一次近拱點龐加萊映射，并將其轉(zhuǎn)化到以行星為原點的慣性坐標系中。

5）在以行星為中心的二體模型中，分別完成擾動流形近拱點軌道狀態(tài)與步驟3）計算的和的匹配，并根據(jù)性能指標（9）確定最佳的機動脈沖?vD和?vA，進一步確定halo出發(fā)和入軌段軌道參數(shù)t0、τD、?VD、?VA和τA。

該搜索方法通過對完整轉(zhuǎn)移軌道進行分段處理，分別在不同的動力學模型下對關(guān)鍵參數(shù)進行搜索，有效降低了轉(zhuǎn)移軌道對關(guān)鍵參數(shù)的敏感度。該搜索算法本質(zhì)上分為兩個層次，第一層是確定日心轉(zhuǎn)移軌道參數(shù)，第二層是確定其它4個軌道段的參數(shù)。通過第一層的搜索，出發(fā)行星和目標行星附近的軌道實現(xiàn)了解耦，可獨立進行關(guān)鍵參數(shù)搜索，降低了搜索難度。

3 數(shù)值仿真分析

本文分別選取了具有較高探測價值的近地小行星2009BD（Aollo型，S類）、2008EA9（Amor，C類）以及99942 Apophis（Aten型，S類），上述目標同時考慮到探測的可行性與科學價值的多樣性，可以有效地覆蓋當今小天體探測所關(guān)注的熱點問題，以2009BD為例，小行星基本軌道參數(shù)如表1。

表1 目標小行星軌道根數(shù)Table 1 The orbital elements of the target asteroid

選取日地 L1點z方 向幅值為1 .6×105km的北向halo軌道作為出發(fā)軌道，微納飛行器逃逸地球影響球的時間區(qū)間為2018年1月1日至2018年12月31日，日心轉(zhuǎn)移段最大飛行時間為500天。速度攝動 ?V搜索范圍為10?6～300m/s。采用本文所提搜索方法，可以獲得采用傳統(tǒng)不變流形時在地球和小行星交會施加機動脈沖之和的等高線圖（圖6）。

圖6 地球–目標小行星（2009BD為例）halo軌道間轉(zhuǎn)移近拱點速度增量等高線圖Fig.6 The contour map of the periapsis velocity increment to the target 2009 BD asteroid based on the halo orbit

由圖6可以看出，等高線圖將解空間分成兩個區(qū)域，并且上半?yún)^(qū)的總速度增量更小，可以尋找到總速度增量在 3 km/s以下的轉(zhuǎn)移機會，但飛行時間較長，約為250天。利用該等高線圖，可以容易得到總速度增量最小解，約為 2 .72km/s。與之對應的地球逃逸時間為2018年6月7日，日心轉(zhuǎn)移段飛行時間263天，總速度增量為0.40 km/s。根據(jù)該轉(zhuǎn)移機會，采用軌道匹配算法對擾動流形參數(shù) τ和 ?V進行搜索，分別繪制了在地球和小行星附近需要施加的速度脈沖和探測器在擾動流形上飛行時間的能量等高線圖，如圖6～7所示。

圖7 地球近拱點速度增量等高線圖Fig.7 The contour map of the periapsis velocity increment near the Earth

由圖7和圖8可以看出，對于日地halo軌道出發(fā)段，在近拱點需要施加的速度增量存在多個局部極小和極大值，而由出發(fā)halo軌道到地球近拱點的飛行時間則隨著攝動速度?V的增大而增大。選擇合適的轉(zhuǎn)移軌道需要綜合考慮機動速度增量和飛行時間，對于圖7中的3個局部極小值，解1的飛行時間過長，解2具有最小的速度增量，飛行時間約為140天，解3雖然速度增量有所增加，但飛行時間則大大減少，約為74天。因此，可以選取解2和解3作為可行的出發(fā)halo軌道離軌機會。通過以上分析，可以獲得可行的出發(fā)halo軌道離軌機會，如表2所示。

圖8 出發(fā)halo軌道到目標小行星的飛行時間Fig.8 The flight time from the halo orbit to the target asteroid

表2 出發(fā)halo軌道離軌機會Table 2 The departure chances from halo orbits

表2中IM為傳統(tǒng)不變流形解，PM為擾動流形解，tD和tA分別為探測器在出發(fā)擾動流形和捕獲擾動流形上的飛行時間。可以看出，EPM1不僅比EIM大大減少了微納飛行器在流形上的飛行時間，并且在地球近拱點需要施加的速度增量更小，這主要是由于速度攝動改變了halo軌道流形的近拱點分布，增加了解的多樣性。同樣，MPM1和MPM2相比EIM也同時節(jié)省了飛行時間和速度增量，這充分證明了擾動流形在構(gòu)建快速低能量轉(zhuǎn)移軌道方面的優(yōu)勢。由表2可以看到，由于引入了行星借力機制，使得轉(zhuǎn)移軌道的燃料消耗大幅降低。另一方面，速度擾動有效地減少了探測器在流形上的飛行時間，獲得的轉(zhuǎn)移軌道兼具快速低能的特性。

通過上述獲取的離軌機會參數(shù)，微納飛行器轉(zhuǎn)移至目標小行星（以2009BD為例）軌跡如圖9～10所示。

圖9 從日地halo軌道轉(zhuǎn)移至目標小行星軌跡Fig.9 The trajectory from the halo orbit in the Sun-Earth system to the target asteroid

圖10 從地月halo軌道轉(zhuǎn)移至目標小行星軌跡前段Fig.10 The trajectory from the Earth-Moon system to the target asteroid

通過對上述3顆目標小行星進行轉(zhuǎn)移軌跡設(shè)計，可以發(fā)現(xiàn)總速度增量明顯小于前人軌跡設(shè)計方案[11]，如表3所示。

表3 轉(zhuǎn)移至目標小行星總速度增量Table 3 The total velocity increment during the transfer process to the target asteroids

4 結(jié) 論

本文首次提出了基于微納飛行器的多小天體并行探測軌跡設(shè)計方案，所做的研究工作具有如下幾個特點：①針對傳統(tǒng)多小天體串行探測任務設(shè)計復雜、周期長，成本高的特點，本文利用搭載微納飛行器的母艦集群停泊在日地/地月halo軌道的方式，提出通過在線選取探測目標，逐個釋放微納飛行器，與目標小行星進行交會探測的方式，有效解耦多小天體探測任務、降低任務成本并且達到探測目標的多樣化，實現(xiàn)一次發(fā)射，多目標探測的“并行”探測方式；②不變流形和行星借力是實現(xiàn)低能量深空探測任務的重要天體力學機理，本文基于不變流形的演化規(guī)律將兩種機理結(jié)合起來，構(gòu)建了大幅度降低燃料消耗的轉(zhuǎn)移軌道方案；③針對傳統(tǒng)不變流形轉(zhuǎn)移時間長的問題，引入擾動流形思想，通過施加較大的速度攝動減小周期軌道的穩(wěn)定系數(shù)，縮短了探測器在流形上的飛行時間，并增加了解空間的多樣性；④提出了一種多目標小行星轉(zhuǎn)移軌跡設(shè)計搜索方法，可對星歷模型下結(jié)合擾動流形與行星借力機制的轉(zhuǎn)移軌道進行大范圍全局搜索，為探測任務設(shè)計提供可行的快速低能量轉(zhuǎn)移方案初始解集。相比傳統(tǒng)方案，本文獲得的轉(zhuǎn)移軌道在燃料消耗和飛行時間上都具有較大的改進，可以為未來的小天體商業(yè)勘查提供有益參考。