基于月球借力的低能DRO 入軌策略

張晨，張皓

中國科學(xué)院 空間應(yīng)用工程與技術(shù)中心，北京 100094

隨著航天技術(shù)的發(fā)展，人類的生存空間逐漸由近地空間拓展到地月空間，2017 年美國國家航空航天局（National Aeronautics and Space Administration，NASA）提出的“ 阿爾忒彌斯（Artemis）”任務(wù)［1］將首次實(shí)現(xiàn)人類在月球及其軌道上的長期居留，并為未來復(fù)雜的載人火星任務(wù)奠定基礎(chǔ)。2021 年中國與俄羅斯就建設(shè)國際月球科研站達(dá)成了合作協(xié)議［2］?？梢灶A(yù)見在未來很長一段時(shí)間內(nèi)，地月空間都將是人類航天活動的主要內(nèi)容。在地月空間部署大型月球軌道站可以顯著降低月球開發(fā)成本并作為其他深空任務(wù)的跳板［3-5］，停泊軌道的選擇不能以單方面性能作為評價(jià)指標(biāo)，而應(yīng)權(quán)衡各種任務(wù)場景以達(dá)到綜合性能最優(yōu)。其中近直線暈軌道（Near Rectilinear Halo Orbit，NRHO）和遠(yuǎn)距離逆行軌道（Distant Retrograde Orbit，DRO）是兩類適合部署月球軌道站的三體周期軌道族，這兩類軌道具有一定的相似性，但側(cè)重點(diǎn)略有不同，NRHO 適合用于載人登月任務(wù)，而DRO 更能兼顧后續(xù)的小行星及火星任務(wù)，本文僅對基于DRO 的月球軌道站部署策略展開研究。

近年來多個(gè)任務(wù)圍繞地月三體系統(tǒng)中的DRO軌道展開，這類軌道具有長期穩(wěn)定性好和入軌代價(jià)低的特點(diǎn)。例如“小行星重定向任務(wù)（Asteroid Redirect Mission，ARM）［6］”計(jì)劃從小行星上抓取一塊巨石并拖曳至DRO 軌道上，DRO 當(dāng)中的某些軌道具有長達(dá)100 a 的穩(wěn)定性而不需要進(jìn)行頻繁的軌道維持。Artemis 的首星［7］將把獵戶座飛船送入25.5 d 的任務(wù)軌道，其中6 d 將位于DRO軌道。面向未來DRO 月球軌道站的在軌建造和貨物補(bǔ)給任務(wù)，提高DRO 入軌質(zhì)量是重要問題。

1987 年Belbruno［8］提出的弱穩(wěn)定邊界（Weak Stability Boundary，WSB）轉(zhuǎn)移技術(shù)在深空任務(wù)中獲得巨大成功，WSB 轉(zhuǎn)移依靠多個(gè)引力體的相互作用，利用混沌區(qū)域所展現(xiàn)的一些新的動力學(xué)特性降低飛行器入軌脈沖。WSB 轉(zhuǎn)移最早用于奔月軌道設(shè)計(jì)，通過選擇合適的地球發(fā)射條件，衛(wèi)星首先抵達(dá)日地引力混沌的弱穩(wěn)定邊界（大約3～5 倍地月距），衛(wèi)星在這個(gè)區(qū)域無動力飛行時(shí)，太陽引力將近地點(diǎn)從地球附近抬升至月球軌道附近，當(dāng)衛(wèi)星從外部返回月球時(shí)被低能捕獲［9］，盡管這種方法需要80～120 d 的飛行時(shí)間，但是能夠節(jié)省大約25%的月球制動脈沖。1991年WSB 轉(zhuǎn)移首次應(yīng)用于日本Hiten［10］衛(wèi)星的營救任務(wù)，隨后GRAIL（Grarity Recovery and Interior Laboratory）［11］、ESMO（the European Student Moon Orbiter）［12］、CAPSTONE（Cislunar Autonomous Positioning System Technology Operations and Navigation Experiment）［13］、KPLO（Korea Pathfinder Lunar Orbiter ）［14］等任務(wù)都采用或計(jì)劃采用WSB 轉(zhuǎn)移技術(shù)。

國內(nèi)外學(xué)者圍繞WSB 轉(zhuǎn)移軌道的設(shè)計(jì)方法展開了內(nèi)容豐富的研究，Koon 等［15］通過平動點(diǎn)附近的不變流形對WSB 的轉(zhuǎn)移機(jī)理進(jìn)行了解釋，并采用雙三體流形拼接方法構(gòu)建了低能奔月軌道。Yagasaki［16］在四體簡化模型下將地月WSB 轉(zhuǎn)移問題轉(zhuǎn)換成非線性規(guī)劃問題，通過數(shù)值迭代求解了地月低能轉(zhuǎn)移軌道。WSB 轉(zhuǎn)移不僅適用于地月低能轉(zhuǎn)移，也適用于DRO、NRHO、L4/5等地月空間典型周期軌道。Xu M 和Xu S J［17］、Tan 等［18］利 用L1/2附近不變流形設(shè)計(jì)了基于WSB 轉(zhuǎn)移的DRO 入軌策略。Scheuerle等［19］基于多步打靶和和數(shù)值延拓在四體簡化模型下計(jì)算了地球至DRO 的低能轉(zhuǎn)移軌道。Zhang 和Hou［20］在雙圓四體模型下設(shè)計(jì)了地球至L4/5的低能轉(zhuǎn)移軌道。

為了最大化飛行器入軌質(zhì)量，一些學(xué)者指出在WSB 轉(zhuǎn)移軌道的地球出發(fā)段增加月球借力（Lunar Gravity Assist，LGA）能夠有效降低火箭發(fā)射能量。為方便描述，將這種類型的轉(zhuǎn)移軌道稱為“LGA+WSB”轉(zhuǎn)移方式?！癓GA+WSB”最早出現(xiàn)于Belbruno 和Miller 在1993 年發(fā)表的文章［10］。Topputo［21］在雙圓四體模型下探索了兩脈沖地月低能轉(zhuǎn)移軌道的全局解空間，并對不同轉(zhuǎn)移方式進(jìn)行了分類，指出“LGA+WSB”轉(zhuǎn)移方式具有較低的火箭發(fā)射脈沖。Parrish 等［13］針對“CAPSTONE”任務(wù)對2 種NRHO 入軌方案進(jìn)行了對比，其中“WSB”方案火箭發(fā)射脈沖較高但是發(fā)射窗口寬松（每天都有發(fā)射窗口），“LGA+WSB”方案火箭發(fā)射脈沖較低但是發(fā)射窗口緊張（一個(gè)月一次）。上述研究沒有針對敏感的月球借力進(jìn)行單獨(dú)設(shè)計(jì)和優(yōu)化，轉(zhuǎn)移軌道中僅有很少一部分呈現(xiàn)“LGA+WSB”的轉(zhuǎn)移方式，數(shù)值計(jì)算效率較低。Scheuerle 等［19］在雙圓四體模型下從“WSB”延拓得到“LGA+WSB”轉(zhuǎn)移軌道，但沒有在星歷模型下進(jìn)行修正。Tselousova 等［22］在雙圓四體模型下基于軌道拼接的思路設(shè)計(jì)了“LGA+WSB”轉(zhuǎn)移軌道，設(shè)計(jì)過程較為繁瑣且沒有在星歷下實(shí)現(xiàn)。

為了提高星歷模型下的“LGA+WSB”轉(zhuǎn)移軌道設(shè)計(jì)效率，為了避免復(fù)雜動力學(xué)環(huán)境對轉(zhuǎn)移軌道收斂性所帶來的影響，使用“近月點(diǎn)龐加萊圖”和“v無窮匹配”獲得轉(zhuǎn)移軌道初值，再通過“多步打靶”在星歷下修正轉(zhuǎn)移軌道。改進(jìn)的方法充分優(yōu)化了月球借力參數(shù)，有效改善了“LGA+WSB”轉(zhuǎn)移方式的計(jì)算效率，所得到的轉(zhuǎn)移軌道滿足地球發(fā)射約束和高精度動力學(xué)約束。

本文結(jié)構(gòu)如下：第1 節(jié)介紹火箭發(fā)射能量和載荷質(zhì)量的關(guān)系；第2 節(jié)為動力學(xué)模型和狀態(tài)轉(zhuǎn)移矩陣；第3 節(jié)介紹DRO 軌道；第4 節(jié)介紹“LGA+WSB”轉(zhuǎn)移方式的特點(diǎn)；第5 節(jié)介紹LGA 至DRO 軌道段的設(shè)計(jì)方法；第6 節(jié)介紹近地球軌道（Lwo Earth Orbit，LEO）至LGA 軌道段的設(shè)計(jì)方法；第7 節(jié)介紹多步打靶技術(shù)；第8 節(jié)為數(shù)值仿真；第9 節(jié)為結(jié)論。

1 火箭發(fā)射能量和載荷質(zhì)量

假設(shè)衛(wèi)星初始位于近地圓形停泊軌道，沿切向施加脈沖抬升遠(yuǎn)地點(diǎn)高度。衛(wèi)星在近地點(diǎn)的速度模表示為

式中：ra和rp分別為遠(yuǎn)地點(diǎn)和近地點(diǎn)地心距；μE為地球引力常數(shù)。衛(wèi)星在近地點(diǎn)的脈沖為

發(fā)射能量C3表示為

假設(shè)衛(wèi)星近地點(diǎn)高度為200 km，圖 1 展示了LEO 發(fā)射脈沖Δνdep和發(fā)射能量C3隨遠(yuǎn)地點(diǎn)地心距ra的變化曲線。當(dāng)發(fā)射脈沖Δνdep=3.131 3 km/s時(shí)遠(yuǎn)地點(diǎn)能夠抵達(dá)月球（ra約為4×105km），發(fā)射能量C3約為-2.03（km）2/s2。發(fā)射脈沖為［3.194 3，3.206 3］ km/s 時(shí)遠(yuǎn)地點(diǎn)抵達(dá)WSB 區(qū)域ra約為［1.2，2.0］×106km，發(fā)射能量C3約為［-0.66，-0.39］（km）2/s2。

圖1 Δvdep 和C3 隨遠(yuǎn)地點(diǎn)高度變化曲線Fig 1 Δvdep and C3 according to apogee radius

采用國外公開的火箭運(yùn)力曲線［23］，圖 2 展示了獵鷹9 Block 2 火箭載荷質(zhì)量隨發(fā)射能量的變化曲線。獵鷹9 發(fā)射地月轉(zhuǎn)移軌道的載荷質(zhì)量為2 619.77 kg，直接發(fā)射至WSB 的載荷質(zhì)量為［2 500.24，2 519.34］ kg。因而如果先將衛(wèi)星發(fā)射至地月轉(zhuǎn)移軌道，再借助月球借力將遠(yuǎn)地點(diǎn)抬升至WSB，則能夠節(jié)省60～70 m/s 的速度脈沖，提升載荷質(zhì)量100.43～119.53 kg。

圖2 獵鷹9 載荷質(zhì)量隨C3 變化曲線Fig 2 Rocket capability according to characteristic energy C3 of Falcon 9 Block 2

2 星歷模型和狀態(tài)轉(zhuǎn)移矩陣

在深空任務(wù)設(shè)計(jì)時(shí)，經(jīng)常使用高精度的N體星歷模型對轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)和修正。圖 3 展示了N體星歷模型示意圖，其中k為中心天體，i為衛(wèi)星，j為其他攝動天體。定義衛(wèi)星狀態(tài)x=[rki，vki]?，衛(wèi)星的動力學(xué)方程表示為

圖3 N 體星歷模型Fig 3 N-body ephemeris model

式中：rki為衛(wèi)星i相對于中心天體k的位置矢量；vki為衛(wèi)星i相對于中心天體k的速度矢量；rkj為攝動天體j相對中心天體k的位置矢量，可以通過JPL 的DE430 星歷獲得；rji為衛(wèi)星i相對攝動天體j的位置矢量，表示為

可以發(fā)現(xiàn)攝動天體的位置矢量隨歷元變化，因而N體動力學(xué)模型是歷元的函數(shù)。

軌道修正經(jīng)常使用狀態(tài)轉(zhuǎn)移矩陣這一概念，狀態(tài)轉(zhuǎn)移矩陣在一階程度上刻畫了初始狀態(tài)改變量對終端狀態(tài)改變量所產(chǎn)生的影響［24］。狀態(tài)轉(zhuǎn)移矩陣的動力學(xué)方程表示為

式中：I為單位矩陣；Φ為狀態(tài)轉(zhuǎn)移矩陣；f(x)為衛(wèi)星動力學(xué)方程；A為動力學(xué)方程對狀態(tài)求偏導(dǎo)得到的矩陣；t0為初始時(shí)刻，t為轉(zhuǎn)移時(shí)刻。定義y為狀態(tài)變量x和狀態(tài)轉(zhuǎn)移矩陣Φ構(gòu)成的列向量，其動力學(xué)方程表示為

其中：g(y)為衛(wèi)星狀態(tài)和狀態(tài)轉(zhuǎn)移矩陣的動力學(xué)方程；“vec”是將矩陣變?yōu)榱邢蛄康乃阕印?/p>

3 DRO 軌道

DRO 是圓型限制性三體問題（Circular Restricted Three-Body Problem，CRTBP）中一類穩(wěn)定的平面軌道族［25-28］，在旋轉(zhuǎn)坐標(biāo)系下該軌道沿順時(shí)針（逆行）運(yùn)動。圖 4 展示了DRO 軌道族，軌道顏色用于區(qū)分軌道周期，L1～L5為5 個(gè)平動點(diǎn)，圖中距離單位1 LU=384 400 km。可以發(fā)現(xiàn)振幅越小的DRO 軌道周期越短。DRO 共振比是指軌道周期與月球公轉(zhuǎn)周期之比，具有典型共振比的周期軌道通?？杀ＷC衛(wèi)星與地球、月球具有周期性的幾何關(guān)系，在數(shù)值仿真中用到了共振比2∶1 的DRO 軌道。

圖4 CRTBP 下 的DRO 軌道族Fig 4 DRO family in CRTBP

以CRTBP 的DRO 軌道作為初值，使用多步打靶可獲得星歷下的DRO 軌道，定義為DRO 參考狀態(tài)，則在星歷模型下數(shù)值積分可以得到DRO 軌道上任意一點(diǎn)的狀態(tài)xdro，表示為

式中：φeph表示以為初值，從τ*到τ積分星歷模型所得到的解。

4 轉(zhuǎn)移軌道設(shè)計(jì)

基于月球借力的低能DRO 入軌策略可簡單描述為，衛(wèi)星初始位于近地圓形停泊軌道，施加第1 次脈沖后進(jìn)入地月轉(zhuǎn)移軌道（Lunar Transfer Orbit，LTO）并掠飛月球，衛(wèi)星通過月球借力抬升遠(yuǎn)地點(diǎn)高度并改變軌道傾角，之后抵達(dá)地月3～5 倍距的弱穩(wěn)定邊界，衛(wèi)星在這里持續(xù)受到太陽引力的影響并抬升近地點(diǎn)高度，當(dāng)再次返回地月空間時(shí)施加第2 次脈沖并進(jìn)入DRO 軌道，希望優(yōu)化軌道轉(zhuǎn)移策略使得兩次脈沖和最小。圖 5 展示了“LGA+WSB”的DRO 入軌示意圖，其中，Δvarr為DRO 入軌速度脈沖矢量；Δvdep為LEO 發(fā)射脈沖矢量。

圖5 LGA+WSB 的DRO 入軌示意圖Fig 5 Illustration of LGA+WSB transfer into DRO

“LGA+WSB”的轉(zhuǎn)移方式能夠有效降低任務(wù)總脈沖，但是復(fù)雜的動力學(xué)環(huán)境和月球借力給轉(zhuǎn)移軌道設(shè)計(jì)帶來解空間龐大和數(shù)值敏感等問題。本文將分別設(shè)計(jì)“LGA 至DRO 軌道段”和“LEO 至LGA 軌道段”，再使用多步打靶技術(shù)修正整條轉(zhuǎn)移軌道。

5 LGA 至DRO 軌道段

本節(jié)將構(gòu)造近月點(diǎn)龐加萊圖，使得轉(zhuǎn)移軌道從DRO 逆向積分先抵達(dá)WSB 再抵達(dá)近月點(diǎn)。圖 6 展示了DRO 入軌時(shí)刻示意圖，衛(wèi)星在DRO入軌時(shí)刻的狀態(tài)x(τf)表示為

圖6 DRO 入軌示意圖Fig 6 Illustration of DRO insertion

式中：xdro=[rdro，vdro]?為入軌時(shí)刻DRO 的狀態(tài)；rdro為入軌時(shí)刻DRO 位置矢量；vdro為入軌時(shí)刻DRO 速度矢量；Δvarr為入軌脈沖矢量；Δνarr為入軌脈沖矢量的模；α∈[0，2π)為Δvarr在e1—e2平面的投影與e1的夾角；β∈[ -π/2，π/2]為Δvarr與e1—e2平面的夾角。

圖7 展示了DRO 出發(fā)的近月點(diǎn)龐加萊圖，轉(zhuǎn)移軌道遠(yuǎn)地點(diǎn)抵達(dá)WSB 的約束表示為

圖7 DRO 出發(fā)的近月點(diǎn)龐加萊圖Fig 7 Illustration of perilune Poincare map of DRO

式中：v為衛(wèi)星速度矢量；r為衛(wèi)星位置矢量；a為衛(wèi)星加速度。轉(zhuǎn)移軌道抵達(dá)近月點(diǎn)的約束表示為

式中：rM、vM和aM分別為月球的位置、速度和加速度。受火箭滑行時(shí)間的限制，一些任務(wù)只能在月球升/降交點(diǎn)附近進(jìn)行引力輔助，轉(zhuǎn)移軌道的近月點(diǎn)還需滿足以下約束

由于衛(wèi)星以雙曲線軌道掠飛月球，衛(wèi)星離開月球影響球的雙曲線剩余速度v+∞表示為

式中：x(τlga)為衛(wèi)星在近月點(diǎn)的狀態(tài)；xM(τlga)為月球狀態(tài)為衛(wèi)星相對月球狀態(tài)。

圖8 月球借力和B 平面Fig 8 Illustration of lunar gravity assist and B-plane

雙曲線剩余速度的模表示為

式中：μM為月球引力常數(shù)；νˉ為衛(wèi)星相對月球速度矢量的模；rˉ為衛(wèi)星位置矢量的標(biāo)量。偏心率矢量e表示為

其中：e為偏心率的標(biāo)量大小。

衛(wèi)星離開影響球的雙曲線剩余速度v+∞表示為

6 LEO 至LGA 軌道段

構(gòu)造LEO 至LGA 的轉(zhuǎn)移軌道，期望在降低地球發(fā)射脈沖的同時(shí)，優(yōu)化月球借力參數(shù)，使得借力后的雙曲線剩余速度與相匹配。在二體簡化模型下構(gòu)造全局優(yōu)化問題為

式中：z∈Rn是優(yōu)化變量；J：Rn→R是目標(biāo)函數(shù)；zlr∈Rn和zur∈Rn是優(yōu)化變量邊界。優(yōu)化變量z表示為

式中：Ω和ω分別為LTO 的升交點(diǎn)赤經(jīng)和近地點(diǎn)幅角；η是LTO 的飛行時(shí)間；γ和ρ分別為月球引力輔助的相位角和月心距。表 1 列出了優(yōu)化變量和取值范圍。

表1 LEO 至LGA 軌道段優(yōu)化變量Table 1 Variables of LEO to LGA segment

假設(shè)地球圓形停泊軌道的地心距為R*，軌道傾角為i*，通過Ω和ω得到衛(wèi)星在LTO 的出發(fā)狀態(tài)x(τ1)，基于二體圓錐曲線拼接假設(shè)，求解Lambert 問題得到地球出發(fā)脈沖Δνdep和衛(wèi)星進(jìn)入月球影響球的雙曲線剩余速度圖 9 展示了B 平面和月球借力參數(shù)，其中i1，i2，i3為坐標(biāo)系的單位矢量。通過γ和ρ得到月球借力后的雙曲線剩余速度，計(jì)算過程如下：

圖9 B 平面和月球借力參數(shù)Fig 9 B-plane and lunar gravity assist parameters

為了降低地球出發(fā)脈沖Δνdep并使得盡可能接近，構(gòu)造目標(biāo)函數(shù)

需要說明的是，式（32）使得月球借力前后的2 段軌道通過v無窮進(jìn)行匹配（不是位置速度拼接）。當(dāng)全局優(yōu)化問題收斂后，可得到星歷模型下的“LGA+WSB”轉(zhuǎn)移軌道初值，計(jì)算過程不再贅述。

7 多步打靶修正

多體系統(tǒng)復(fù)雜的動力學(xué)環(huán)境使得轉(zhuǎn)移軌道對初值非常敏感，多步打靶技術(shù)可以擴(kuò)大收斂域并提高算法的魯棒性。圖 10 展示了多步打靶示意圖，整條軌道被n個(gè)離散點(diǎn)分割為n-1 條軌道段，多步打靶變量X表示為

圖10 多步打靶示意圖Fig 10 Illustration of multiple shooting

第j段軌道滿足狀態(tài)連接約束：

第j段軌道滿足時(shí)間連接約束：

第j段軌道滿足飛行時(shí)間為正的不等式約束：

引入松弛變量βj可將式（36）變?yōu)榈仁郊s束：

第1 個(gè)節(jié)點(diǎn)滿足邊界約束：

式中：R*和i*為LEO 的地心距和軌道傾角。

第n個(gè)節(jié)點(diǎn)滿足邊界約束：

至此，當(dāng)?shù)玫窖趴吮染仃嘕(X)，使用最小二乘更新方程對打靶變量Xk+1進(jìn)行修正（見式（41）），收斂后即可得到轉(zhuǎn)移軌道。在星歷修正中，動力學(xué)模型僅考慮日、地、月引力場，行星狀態(tài)通過JPL 的DE430 星歷計(jì)算。

衛(wèi)星在任務(wù)過程中共施加兩次脈沖，分別為火箭發(fā)射脈沖Δvdep和衛(wèi)星平臺施加的DRO 入軌脈沖Δνarr，任務(wù)脈沖總表示為

8 數(shù)值仿真

將在星歷模型下設(shè)計(jì)基于月球借力的低能DRO 轉(zhuǎn)移軌道，仿真參數(shù)參見表 2。其中，時(shí)間單位中的jd 為儒略日，tdb 為太陽系質(zhì)心動力學(xué)時(shí)。計(jì)算流程如下：

表2 數(shù)值仿真參數(shù)Table 2 Simulation parameters

1）網(wǎng)格化DRO 入軌時(shí)間τf和入軌速度脈沖Δvarr，獲得不同的DRO 入軌狀態(tài)x(τf)，參見式（9）和式（10）。

2）從x(τf)出發(fā)逆向積分動力學(xué)模型，參見式（8），使得衛(wèi)星先抵達(dá)WSB 再抵達(dá)近月點(diǎn)（構(gòu)造近月點(diǎn)龐加萊圖），參見式（13）～式（18），計(jì)算月球借力后的雙曲線剩余速度，參見式（20）～式（24）。

3）固定LGA 至DRO 軌道段，使用粒子群優(yōu)化LEO 至LGA 軌道段，最小化地球出發(fā)脈沖Δνdep并使得盡可 能接近（即v無窮 匹配），參見式（25）～式（32）。

4）使用多步打靶修正，得到“LGA+WSB”類型的DRO 轉(zhuǎn)移軌道，參見式（33）～式（41）。

根據(jù)以上計(jì)算流程，構(gòu)建1×107個(gè)DRO 入軌狀態(tài)，從DRO 逆向積分能夠先抵達(dá)WSB 再抵達(dá)近月點(diǎn)的僅有1 511 條轉(zhuǎn)移軌道。圖 11 展示了近月點(diǎn)龐加萊圖，圖中每個(gè)圓點(diǎn)都代表一條轉(zhuǎn)移軌道，圓點(diǎn)顏色用于區(qū)分DRO 入軌脈沖。

圖11 近月點(diǎn)龐加萊圖Fig 11 Perilune Poincare map

進(jìn)一步對比“v無窮匹配”策略對計(jì)算效率的改進(jìn)效果。首先對于圖 11 上的每一個(gè)點(diǎn)，繼續(xù)逆向積分并在近地點(diǎn)終止，得到不使用“v無窮匹配”的轉(zhuǎn)移軌道初值。圖 12 展示了這些初值的近地點(diǎn)龐加萊圖，X軸為近地點(diǎn)地心距，Y軸為近地點(diǎn)軌道傾角，紅色十字為期望的LEO 地心距和軌道傾角，可以發(fā)現(xiàn)敏感的月球借力使得紅色十字附近幾乎沒有接近的轉(zhuǎn)移軌道。如果對圖 12 直接使用多步打靶軌道修正，則僅有49 條收斂的轉(zhuǎn)移軌道。

圖12 近地點(diǎn)龐加萊圖Fig 12 Perigee Poincare map

接著采用“v無窮匹配”策略獲得轉(zhuǎn)移軌道初值，再使用多步打靶修正整條軌道，最終獲得了743 條收斂的轉(zhuǎn)移軌道。圖 13 展示了這些收斂解的散點(diǎn)圖，其中每個(gè)圓點(diǎn)都為一條收斂軌道，圓點(diǎn)顏色用于區(qū)分DRO 入軌脈沖，紅色圓圈為收斂解的帕累托前沿，紅色五角星為脈沖最低解。在數(shù)值仿真中采用Matlab 編程語言，計(jì)算設(shè)備為高性能工作站，CPU 為64 核心處理器，內(nèi)存為256 G，計(jì)算總時(shí)間約為8 h。

圖13 多步打靶修正后的收斂解Fig 13 Converged solutions with multiple shooting

需要說明的是，“v無窮匹配”策略能夠在不改變月球借力后雙曲線剩余速度的條件下，將圖12 的大多數(shù)近地點(diǎn)都移動到紅色十字附近，從而有效改進(jìn)了轉(zhuǎn)移軌道初值。就本算例而言，“v無窮匹配”策略將轉(zhuǎn)移軌道收斂率提高了約15 倍。此外，“v無窮匹配”策略可以方便拓展到NRHO、L4/5等周期軌道的入軌任務(wù)設(shè)計(jì)，計(jì)算流程并無明顯區(qū)別。

圖14 展示了圖 13 中總脈沖最低解（紅色五角星）的轉(zhuǎn)移軌道，其中紅色實(shí)線為DRO 軌道，藍(lán)色實(shí)線為地球至DRO 轉(zhuǎn)移軌道，灰色實(shí)線為月球軌道。觀察圖 14（a），衛(wèi)星從地球出發(fā)后進(jìn)入地月轉(zhuǎn)移軌道，月球借力后衛(wèi)星彈射至WSB，當(dāng)再次返回地月空間時(shí)，衛(wèi)星以切向入軌的方式進(jìn)入DRO 軌道。觀察圖 14（b），衛(wèi)星從WSB 返回地月空間時(shí)首先抵達(dá)L4附近，之后在地月引力共同影響下抵達(dá)DRO 并入軌。

圖14 總脈沖最低解的轉(zhuǎn)移軌道Fig 14 Best solution with minimum cost

表3展示了圖 13中脈沖最低解（紅色五角星）的仿真結(jié)果，LEO 出發(fā)脈沖僅需3.127 3 km/s，DRO入軌脈沖66.1 m/s，任務(wù)總脈沖3.193 4 km/s，任務(wù)總時(shí)間102.880 3 d。

表3 總脈沖最低解的仿真結(jié)果Table 3 Best solution with minimum cost

表4 展示了LEO 轉(zhuǎn)移至GEO 和DRO 的 對比，LEO 軌道高度200 km，軌道傾角28.5°，衛(wèi)星平臺采用化學(xué)推進(jìn)，比沖220 s。如果采用獵鷹9號Block 2 型運(yùn)載火箭，單次發(fā)射至多可將1 704.273 3 kg 的載荷送入GEO 軌道，可將2 553.576 8 kg 的載荷送入DRO 軌道，由于GEO 需要1.819 7 km/s 的速度脈沖改變軌道傾角和軌道圓化，使得雖然DRO 距離地球更遠(yuǎn)，但是入軌質(zhì)量卻更高，約為GEO 的1.5 倍。

表4 LEO 轉(zhuǎn)移至GEO 和DRO 的對比Table 4 Comparison between LEO to GEO and LEO to DRO

9 結(jié)論

1）改進(jìn)了星歷模型下的“LGA+WSB”轉(zhuǎn)移軌道設(shè)計(jì)方法，通過“近月點(diǎn)龐加萊圖”和“v無窮匹配”獲得較好的軌道初值，有效提高了該類型轉(zhuǎn)移軌道的計(jì)算效率。

2）通過網(wǎng)格搜索獲得了解空間的帕累托前沿，對于2：1 DRO 軌道，總脈沖最低解的飛行時(shí)間102.880 3 d，地球發(fā)射脈沖3.127 3 km/s，DRO 入軌脈沖66.1 m/s。

3）對 比LEO 至GEO 和LEO 至DRO 轉(zhuǎn) 移軌道，盡管DRO 距離地球更遠(yuǎn)，但是入軌質(zhì)量可以達(dá)到GEO 的約1.5 倍。