一種離散軌道數(shù)據(jù)約束下的地月三體軌道脈沖轉(zhuǎn)移算法

王義宇，羅宇航，徐田來，包為民，袁 帥，張澤旭，李宸碩，胡志杰

（1.哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001；2.中國航天科技集團(tuán)有限公司，北京 100048）

引言

近年來，不同太空機(jī)構(gòu)使用多種類型的航天器重新訪問了月球，各航天強(qiáng)國均提出更富挑戰(zhàn)性的探月計劃[1-4]。不管是針對月球的長期探測和建設(shè)月球空間站等在軌任務(wù)，還是針對載人登月和建造月球基地等著陸任務(wù)，借助地月空間內(nèi)的多種特殊軌道從而以更低的成本完成任務(wù)成為關(guān)鍵問題，而其中設(shè)計各個特殊軌道之間的轉(zhuǎn)移軌道更是所有任務(wù)的重要基礎(chǔ)[5-6]。

目前，地月空間內(nèi)的局部周期軌道得到了廣泛關(guān)注?？紤]到軌道能量與穩(wěn)定性，研究對象主要以繞月的大幅值逆行軌道（Distant Retrograde Orbit，DRO）和繞平動點(diǎn)的Lissajous軌道、Halo軌道族以及其中較為特殊的近直線暈軌道（Near-Rectilinear Halo Orbit，NRHO）為主[7-8]。

各學(xué)者根據(jù)不同軌道的不同特性設(shè)計了多種脈沖轉(zhuǎn)移軌道的建模和計算方法：孫俞等[9]針對平動點(diǎn)軌道間兩脈沖轉(zhuǎn)移問題設(shè)計了三體Lambert求解算法、潘迅等[10]基于主矢量理論提出一種快速計算雙脈沖轉(zhuǎn)移軌道的方法、Tan等[11]利用周期軌道研究了低能量雙脈沖轉(zhuǎn)移問題、 Zimovan-Spreen等[12]對NRHO和DRO進(jìn)行了動力學(xué)分析并根據(jù)流形設(shè)計了轉(zhuǎn)移軌道、曾豪等[13]采用配點(diǎn)法與局部優(yōu)化算法構(gòu)造出滿足近月點(diǎn)約束的環(huán)月軌道與平動點(diǎn)軌道之間燃料最優(yōu)往返轉(zhuǎn)移、Oshima[14]以平面李雅普諾夫軌道的垂直穩(wěn)定流形為基礎(chǔ)設(shè)計了NRHO與DRO之間的轉(zhuǎn)移軌道，也有學(xué)者提出采用深度強(qiáng)化學(xué)習(xí)如PPO、MCTS等普適方法設(shè)計地月空間轉(zhuǎn)移軌道[15-16]。

從研究對象看，已有研究主要集中于地月停泊軌道與特殊軌道之間的轉(zhuǎn)移，而對于特殊軌道之間的轉(zhuǎn)移問題研究較少；從優(yōu)化問題模型看，對于兩脈沖轉(zhuǎn)移，大多建模為三體模型下的兩點(diǎn)邊值問題[17-19]，多脈沖轉(zhuǎn)移?？紤]脈沖點(diǎn)位置等特殊約束并建模為非線性優(yōu)化問題[11,20]；從求解算法看，一般采用數(shù)值迭代方法或成熟的非線性求解器進(jìn)行求解[20]，但解的收斂性和最優(yōu)性對初值依賴度較高，求解難度較大，一般會采用智能算法（如粒子群算法、遺傳算法、差分進(jìn)化算法等[21-22]）或根據(jù)起始和目標(biāo)軌道特性[13,20]進(jìn)行初值的簡單搜索，普適性較低。

基于以上問題，本文以地月空間中特殊軌道之間轉(zhuǎn)移為研究對象，引入離散軌道數(shù)據(jù)作為轉(zhuǎn)移軌道的起點(diǎn)和終點(diǎn)約束，提高對不同軌道的適應(yīng)能力，并將最小脈沖轉(zhuǎn)移問題進(jìn)行拆分重構(gòu)，轉(zhuǎn)化為雙層優(yōu)化問題，降低問題求解難度，同時兼顧全局最優(yōu)性和計算效率，在雙層優(yōu)化模型的基礎(chǔ)上提出了雙層優(yōu)化求解算法，避免了初值敏感問題，實(shí)現(xiàn)了多種類型軌道間的最小脈沖轉(zhuǎn)移軌道設(shè)計。

1 圓形限制性三體問題模型

1.1 動力學(xué)模型

圓型限制性三體問題（Circular Restricted ThreeBody Problem，CRTBP）是對實(shí)際三體問題的一種理想近似，是研究三體軌道動力學(xué)行為的基礎(chǔ)。在CRTBP下，可定義會合坐標(biāo)系B-xyz，其坐標(biāo)原點(diǎn)在三體系統(tǒng)的質(zhì)心B上，x軸由系統(tǒng)質(zhì)心B指向第二主天體的質(zhì)心，xy坐標(biāo)面即兩個主天體相對運(yùn)動平面，x軸與y、z軸構(gòu)成右手坐標(biāo)系，如圖1所示。

圖1 圓型限制性三體動力學(xué)模型Fig.1 Circular restricted three-body dynamic model

在三體會合坐標(biāo)系中，航天器的動力學(xué)模型可描述為[9]

其中：x，y，z表示航天器在三體會合坐標(biāo)系中的位置；μ表示兩主天體的質(zhì)量比參數(shù)；r1、r2分別表示航天器到兩主天體的距離。

另外在CRTBP模型下，存在雅可比常數(shù)可表征三體系下航天器所在軌道的能量

1.2 CRTBP下特殊軌道

圓型限制性三體問題動力學(xué)模型具有5個平衡點(diǎn)，也稱為平動點(diǎn)。圍繞平動點(diǎn)與主天體也存在著很多特殊的周期解，構(gòu)成了三體模型中特殊的周期軌道和擬周期軌道。主要包含繞平動點(diǎn)的周期軌道如Halo、Lissajous軌道和繞天體的周期軌道如大幅值逆行軌道等，根據(jù)不同要求，往往可得到同類型的一族軌道。圖2中繪出了典型的局部周期軌道族（Halo、Lissajous、DRO），一般對會合坐標(biāo)系下速度、距離等進(jìn)行歸一化處理，以下仿真圖中距離的單位取為歸一化單位，即1LEM=384 400 km。

圖2 地月三體系統(tǒng)下的局部周期軌道示意圖Fig.2 Local periodic orbit under Earth-Moon CRTBP

Halo軌道族及其特殊的一類近直線暈軌道在未來具有較為突出的應(yīng)用價值，對月球探測以及月球南北極著陸具有重要意義。而大幅值逆行周期軌道DRO具有較高的穩(wěn)定性，因此適宜于月球空間站或探測器的長期在軌任務(wù)需求[19]。

此外不變流形也是基于CRTBP模型的主要理論之一。在CRTBP模型下，存在穩(wěn)定流形流入平動點(diǎn)軌道，還存在受擾動后流出平動點(diǎn)軌道的不穩(wěn)定流形。對于與其它特殊軌道相接的不變流形可設(shè)計低能轉(zhuǎn)移軌道，對于不與目標(biāo)軌道相接的流形，也可通過較小的速度脈沖實(shí)現(xiàn)軌道轉(zhuǎn)移[23-24]。Halo軌道的穩(wěn)定流形與不穩(wěn)定流形示意圖如圖3所示。

圖3 Halo軌道的穩(wěn)定流形與不穩(wěn)定流形Fig.3 Stable and unstable manifolds of Halo orbits

然而即使CRTBP是三體問題的簡化模型，其動力學(xué)模型仍較為復(fù)雜。對于前文提及的特殊周期軌道以及不變流形，難以給出精確的解析解。目前主流的求解方式仍是在近似解析解的基礎(chǔ)上通過微分修正或多點(diǎn)打靶等方法進(jìn)行積分修正，進(jìn)而得到較為精確的數(shù)值解，并一般以一定時間步長離散的軌道數(shù)據(jù)形式給出。

2 最小脈沖轉(zhuǎn)移問題建模

地月系統(tǒng)下特殊軌道的最小脈沖轉(zhuǎn)移問題可定義為：在多次脈沖消耗總?cè)剂匣蚩偰芰孔钚〉那闆r下，實(shí)現(xiàn)起始軌道到目標(biāo)軌道的轉(zhuǎn)移。本文討論最一般的情況：轉(zhuǎn)移軌道起點(diǎn)與終點(diǎn)不固定、轉(zhuǎn)移時間不固定、起始軌道與目標(biāo)軌道類型不固定?？紤]理想假設(shè)：在脈沖時刻，航天器在地月三體系統(tǒng)空間中的位置連續(xù)、速度產(chǎn)生突變。

為了盡可能地保證軌道的精確，本文在進(jìn)行轉(zhuǎn)移軌道的設(shè)計時，將起始軌道與目標(biāo)軌道以離散數(shù)據(jù)構(gòu)成的數(shù)據(jù)表的形式給出，從而約束轉(zhuǎn)移軌道起點(diǎn)與終點(diǎn)的位置與速度。這種方法不僅保證了起點(diǎn)與終點(diǎn)的精度，同時也極大拓展了算法的普適性。起始與目標(biāo)軌道不僅可以是特定類型的某一條軌道、某一類型的軌道族，同時也適用于借助不變流形的轉(zhuǎn)移軌道設(shè)計，即不穩(wěn)定流形作為起始軌道或穩(wěn)定流形作為目標(biāo)軌道。

考慮最一般的n+1次脈沖情況，脈沖時刻分別為t=[t0,t1,t2,···,tn]1×(n+1)，分別對應(yīng)著速度改變量為Δv=[Δv0,Δv1,Δv2,···,Δvn]3×(n+1)?？啥x指標(biāo)為燃料最優(yōu)或能量最優(yōu)[25]

約束條件主要包括動力學(xué)約束、始末狀態(tài)約束、脈沖約束、脈沖點(diǎn)約束以及轉(zhuǎn)移時間約束等。

1）動力學(xué)約束

對圓型限制性三體問題下的動力學(xué)方程式（1）進(jìn)行變換，可寫為

定義各個脈沖時刻對應(yīng)的航天器狀態(tài)為

X0,X1,X2,···,Xn

在沒有脈沖作用時，航天器應(yīng)始終滿足動力學(xué)方程，可表示為積分形式

2）始末狀態(tài)約束

始末狀態(tài)約束要求轉(zhuǎn)移軌道的起點(diǎn)位于起始軌道上、終點(diǎn)位于目標(biāo)軌道上。由前文提出本文起始軌道與目標(biāo)軌道基于離散軌道數(shù)據(jù)構(gòu)成的數(shù)據(jù)表給出，可分別定義起點(diǎn)可行數(shù)集與終點(diǎn)可行數(shù)集

其中：p、q為軌道數(shù)據(jù)的離散點(diǎn)數(shù)目；Xcj與Xtl分別對應(yīng)離散的軌道數(shù)據(jù)，且每點(diǎn)應(yīng)包含時刻、位置與速度，共7維，即

轉(zhuǎn)移軌道起點(diǎn)和軌道終點(diǎn)為X0和Xn，即第一次脈沖時刻以及最后一次脈沖時刻對應(yīng)的航天器狀態(tài)。始末狀態(tài)約束即要求航天器在起始時刻未施加脈沖時的狀態(tài)是起始軌道數(shù)據(jù)中的一點(diǎn)，且在終端時刻施加脈沖后的狀態(tài)是目標(biāo)軌道數(shù)據(jù)中的一點(diǎn)，則約束可表述為

3）其它約束

對于各個脈沖時刻要求位置連續(xù)，脈沖前后可表示為

此外還可能對脈沖點(diǎn)存在其它限制，如要求脈沖點(diǎn)的位置等?？山y(tǒng)一寫為如下形式

考慮脈沖上限，可表示為如下約束

考慮轉(zhuǎn)移時間上限，同時時間序列應(yīng)為遞增順序，可表示為如下約束

綜上，可建立最小脈沖轉(zhuǎn)移問題如問題1所示。

問題1：

3 脈沖轉(zhuǎn)移問題重構(gòu)與算法設(shè)計

3.1 雙層優(yōu)化問題建模

將所有優(yōu)化變量拆分為兩部分：起點(diǎn)與終點(diǎn)的狀態(tài)、脈沖時間序列為上層變量，脈沖序列為下層變量。相應(yīng)的，問題1的約束條件也可拆分為兩部分，從而重構(gòu)為雙層優(yōu)化問題如問題2所示。

問題2：

雙層優(yōu)化實(shí)質(zhì)上是將最小脈沖轉(zhuǎn)移問題拆分為兩步進(jìn)行：第一步在起始軌道和目標(biāo)軌道上選點(diǎn)，并給出這兩點(diǎn)之間的脈沖時間序列；第二步建立最優(yōu)控制問題模型，在確定起點(diǎn)和終點(diǎn)以及轉(zhuǎn)移時間的情況下僅優(yōu)化控制量即優(yōu)化脈沖序列，優(yōu)化指標(biāo)為燃料最優(yōu)或能量最優(yōu)。上層問題求解方法選用具有自治能力的粒子群算法（Adaptive Genetic Particle Swarm Optimization，AGPSO）的智能算法，下層問題求解方法選用序列二次規(guī)劃算法（Sequential Quadratic Programming，SQP）。

3.2 上層優(yōu)化問題求解算法設(shè)計

由于行數(shù)/列數(shù)與此行/列的數(shù)據(jù)存在固定的映射關(guān)系，不同行/列的映射關(guān)系并不相同，每個行數(shù)/列數(shù)與其對應(yīng)的數(shù)據(jù)均可描述為某個特定的未知函數(shù)

對于解析的迭代方法來說，不能顯式表達(dá)的函數(shù)往往不利于求解。但得益于智能算法的普適性，可通過尋找行數(shù)或者列數(shù)實(shí)現(xiàn)降維：原本要優(yōu)化的起點(diǎn)和終點(diǎn)狀態(tài)共12維的變量，通過表格數(shù)據(jù)的自然約束降維為2維，即起始軌道數(shù)據(jù)表中的某行或某列以及目標(biāo)軌道數(shù)據(jù)表中的某行或某列（對于一族軌道數(shù)據(jù)，可通過三維表格進(jìn)行約束，此時優(yōu)化變量降維為4維），可建立行數(shù)/列數(shù)的約束如下

因此可重新構(gòu)建上層優(yōu)化問題如下。

問題3：

由于上層優(yōu)化問題優(yōu)化變量與優(yōu)化指標(biāo)關(guān)系不明確，且優(yōu)化變量的初值難以設(shè)計，因此選擇智能算法進(jìn)行求解。綜合考慮計算效率與收斂性，選擇AGPSO進(jìn)行優(yōu)化。在傳統(tǒng)的PSO算法中，所有粒子在局部搜索和全局搜索的效果方面都是相同的，因此可以認(rèn)為所有粒子是一個具有同一種策略的群。各個粒子的位置和速度更新可表示為

然而在任何基于種群的優(yōu)化算法中，使用具有共同目標(biāo)的不同自治群體，理論上可以使算法具備更多的隨機(jī)或定向的搜索。AGPSO的核心思路為利用不同的策略來更新c1和c2，對自治群進(jìn)行數(shù)學(xué)建模[26]。換言之，這些群體在關(guān)注個體和全局表現(xiàn)的程度上表現(xiàn)得不同。本文設(shè)計4個群體，即設(shè)計4種c1和c2的變化策略，更新函數(shù)如表1所示。

總之，這種傳統(tǒng)教學(xué)與PBL教學(xué)相結(jié)合的教學(xué)方式，理論成績與實(shí)踐考核等方面均有提高，學(xué)生滿意度高，既提升了師資隊伍水平，又培養(yǎng)出敏而好學(xué)的學(xué)生，加強(qiáng)了學(xué)生人文素養(yǎng)的培養(yǎng)，激發(fā)了學(xué)習(xí)的熱情。教學(xué)模式有機(jī)的結(jié)合能夠發(fā)揮兩方面的優(yōu)勢，更好地促進(jìn)健康管理的理論與教學(xué)實(shí)踐。

表1 4種 c1和 c2的更新策略Table 1 Four updating formula for c1 and c2

后文仿真中，上層優(yōu)化問題均基于AGPSO算法進(jìn)行求解，其它參數(shù)如表2所示。

表2 AGPSO算法相關(guān)參數(shù)表Table 2 Parameters of AGPSO algorithm

3.3 下層優(yōu)化問題求解算法設(shè)計

經(jīng)過拆分處理后，下層優(yōu)化問題的自由變量僅剩3(n+1)維的脈沖序列 Δv。根據(jù)上層優(yōu)化所得到的起點(diǎn)與終點(diǎn)狀態(tài)以及各個脈沖時刻，可建立最優(yōu)控制問題如下。

問題4：

問題4仍為含約束的非線性優(yōu)化問題，但相比于問題1，變量個數(shù)大大減少、求解難度大大降低。當(dāng)轉(zhuǎn)移過程僅為兩脈沖時，問題可直接建模為三體模型下的Lambert問題，可采用如牛頓–同倫法[27]、擬線性化–局部變分迭代法[28]等算法求解。但涉及多脈沖或包含其它約束限制時，固定轉(zhuǎn)移時間的兩點(diǎn)邊值約束問題求解方法不再適用。針對此類非線性問題，采用SQP算法。求解過程本文不再贅述，相關(guān)參數(shù)選取如表3所示。

表3 SQP算法相關(guān)參數(shù)表Table 3 Parameters of SQP algorithm

雖然SQP算法本身較為成熟，但其仍存在局部收斂、初值敏感等問題[6]。對于本文所構(gòu)建的非線性規(guī)劃問題4來說，由于其優(yōu)化指標(biāo)為優(yōu)化變量的2–范數(shù)之和或2–范數(shù)平方和的形式，所期望脈沖盡可能小，即問題4的解盡可能靠近 0。考慮到SQP算法的局部收斂性，選取初值為Δv=0，認(rèn)為問題4的期望解在其鄰域內(nèi)①下層優(yōu)化問題4的求解受上層優(yōu)化問題3的結(jié)果影響。當(dāng)SQP算法初值選為0時，上層結(jié)果可能導(dǎo)致下層優(yōu)化不收斂，表明此時求解的脈沖較大，不在0的鄰域內(nèi)。同時也表明了當(dāng)次求解上層結(jié)果的不合理，因此可對此結(jié)果賦予大指標(biāo)作為懲罰，進(jìn)入下一次迭代求解。。在一定程度上避免了SQP對初值敏感、難以設(shè)計迭代初值的問題。

通過結(jié)合AGPSO與SQP對問題3與問題4雙層優(yōu)化模型進(jìn)行迭代求解，可得到設(shè)置最大迭代次數(shù)內(nèi)的最小脈沖，整體求解框架如圖4所示。

圖4 基于AGPSO和SQP的雙層優(yōu)化求解框架Fig.4 Two-layer optimization solution framework based on AGPSO and SQP

4 仿真分析

4.1 指定軌道間的轉(zhuǎn)移軌道

對于指定軌道間的脈沖轉(zhuǎn)移問題，本文選取3種軌道轉(zhuǎn)移場景進(jìn)行仿真：①繞L1點(diǎn)NRHO軌道到繞L2點(diǎn)Halo軌道的轉(zhuǎn)移；②繞L1點(diǎn)Lissajous軌道到DRO軌道的轉(zhuǎn)移；③繞L1點(diǎn)NRHO軌道到DRO軌道的轉(zhuǎn)移。

1）繞L1點(diǎn)NRHO到繞L2點(diǎn)Halo的轉(zhuǎn)移

選擇振幅為125 000 km的L1點(diǎn)NRHO和振幅為100 000 km的L2點(diǎn)Halo作為起始軌道和目標(biāo)軌道，離散點(diǎn)均為2 000個。采用兩脈沖與三脈沖兩種方式進(jìn)行軌道轉(zhuǎn)移，指標(biāo)選擇為能量最優(yōu)。通過算法可求得轉(zhuǎn)移軌道分別如圖5（a）、（b）所示，指標(biāo)隨迭代的變化曲線如圖5（c）所示。

圖5 L1點(diǎn)NRHO–L2點(diǎn)Halo的兩脈沖與三脈沖轉(zhuǎn)移Fig.5 Transfer from L1 NRHO to L2 Halo using two impulses/three impulses

表4 L1點(diǎn)NRHO–L2點(diǎn)Halo的脈沖轉(zhuǎn)移優(yōu)化結(jié)果Table 4 Transfer Impulse from L1 NRHO to L2 Halo

2）繞L1點(diǎn)Lissajous軌道到DRO軌道的轉(zhuǎn)移

選擇振幅為125 000 km的L1點(diǎn)Lissajous軌道和振幅為496 600 km的DRO作為起始軌道和目標(biāo)軌道，離散點(diǎn)均為2 000個。同樣采用兩脈沖與三脈沖兩種方式進(jìn)行軌道轉(zhuǎn)移，指標(biāo)選擇為能量最優(yōu)。通過算法可求得轉(zhuǎn)移軌道分別如圖6（a）、（b）所示，指標(biāo)隨迭代的變化曲線如圖6（c）所示。

圖6 L1點(diǎn)Lissajous軌道–DRO的兩脈沖與三脈沖轉(zhuǎn)移Fig.6 Transfer from L1 Lissajous to DRO using two impulses/three impulses

由圖6（a）、（b）所示，轉(zhuǎn)移軌道通過兩脈沖和三脈沖完成L1點(diǎn)Lissajous軌道至DRO的轉(zhuǎn)移，同時實(shí)現(xiàn)起始點(diǎn)、目標(biāo)點(diǎn)、中間脈沖點(diǎn)以及轉(zhuǎn)移時間的選取。根據(jù)圖6（c）可得，隨著算法迭代計算，能量指標(biāo)逐漸收斂至較低的值，從而實(shí)現(xiàn)設(shè)定代數(shù)內(nèi)的最優(yōu)能量轉(zhuǎn)移，可得到計算結(jié)果如表5所示。

表5 L1點(diǎn)Lissajous軌道–DRO的脈沖轉(zhuǎn)移優(yōu)化結(jié)果Table 5 Transfer Impulse from L1 Lissajous orbit to DRO

3）繞L1點(diǎn)NRHO軌道到DRO軌道的轉(zhuǎn)移

選擇振幅為125 000 km的L1點(diǎn)NRHO軌道和振幅為285 600 km的DRO作為起始軌道和目標(biāo)軌道，離散點(diǎn)均為2 000個。采用兩脈沖、三脈沖以及月球借力3種方式進(jìn)行軌道轉(zhuǎn)移，指標(biāo)選擇為能量最優(yōu)。

其中，月球借力轉(zhuǎn)移策略為：在三脈沖轉(zhuǎn)移的基礎(chǔ)上再限制中間脈沖點(diǎn)位置作為問題4的約束條件。要求中間脈沖點(diǎn)為近月點(diǎn)，即脈沖點(diǎn)相對月球中心的位置矢量與速度矢量垂直，同時要求中間脈沖點(diǎn)距月球表面高度小于500 km，可表示為

通過算法可求得轉(zhuǎn)移軌道分別如圖7（a）、（b）、（c）所示，指標(biāo)隨迭代的變化曲線如圖7（d）所示。

圖7 L1點(diǎn)NRHO–DRO的兩脈沖、三脈沖與月球借力轉(zhuǎn)移Fig.7 Transfer from L1 NRHO to DRO using two impulses/three impulses/gravity assist

由圖7（a）、（b）、（c）所示，轉(zhuǎn)移軌道通過兩脈沖、三脈沖以及借月飛行3種方式完成L1點(diǎn)NRHO軌道至DRO的轉(zhuǎn)移，同時實(shí)現(xiàn)起始點(diǎn)、目標(biāo)點(diǎn)、中間脈沖點(diǎn)以及轉(zhuǎn)移時間的選取。根據(jù)圖7（d）可得，隨著算法迭代計算，3種不同轉(zhuǎn)移方式的能量指標(biāo)均能逐漸收斂，可得到計算結(jié)果如表6所示。

表6 L1點(diǎn)NRHO軌道–DRO的脈沖轉(zhuǎn)移優(yōu)化結(jié)果Table 6 Transfer Impulse from L1 NRHO to DRO

另外由圖7（c）可見，月球借力飛行的中間脈沖點(diǎn)也滿足所設(shè)要求，能夠?qū)崿F(xiàn)借助月球引力進(jìn)行轉(zhuǎn)移。但添加近月脈沖的限制條件后，算法在50個粒子、100次迭代的情況下，月球借力飛行的結(jié)果并不如兩脈沖與三脈沖轉(zhuǎn)移。然而根據(jù)圖7（a）、（b）兩圖可直觀看出，兩脈沖與三脈沖轉(zhuǎn)移的軌道實(shí)際上也借助了月球引力從而改變了速度方向，而三脈沖不限制中間點(diǎn)脈沖，可能存在更優(yōu)的轉(zhuǎn)移軌道，可見本文提出的算法具有較強(qiáng)的搜索能力。

4.2 軌道族間的轉(zhuǎn)移軌道

對于軌道族間的脈沖轉(zhuǎn)移問題，本文選取兩種軌道轉(zhuǎn)移場景進(jìn)行仿真：①繞L1點(diǎn)Halo軌道族到繞L2點(diǎn)Halo軌道族的轉(zhuǎn)移；②繞L1點(diǎn)Halo軌道族到DRO軌道族的轉(zhuǎn)移。

1）L1點(diǎn)Halo軌道族到L2點(diǎn)Halo軌道族的轉(zhuǎn)移

選擇振幅范圍為67 850～96 600 km的繞L1點(diǎn)Halo軌道族和振幅范圍為64 350～99 250 km的繞L2點(diǎn)Halo軌道族作為備選的起始軌道和目標(biāo)軌道，離散的軌道數(shù)據(jù)構(gòu)成三維矩陣，維度分別為：6 × 2000 × 80和6 × 2 000 × 80。采用兩脈沖、三脈沖兩種方式進(jìn)行軌道轉(zhuǎn)移，指標(biāo)選擇為能量最優(yōu)。通過算法可求得轉(zhuǎn)移軌道分別如圖8（a）、（b）所示，指標(biāo)隨迭代的變化曲線如圖8（c）所示。

圖8 L1點(diǎn)Halo軌道族–L2點(diǎn)Halo軌道族的兩脈沖與三脈沖轉(zhuǎn)移Fig.8 Transfer from L1 Halo family to L2 Halo family using two impulses/three impulses

由圖8（a）、（b）所示，算法選擇了兩個軌道族中距離最接近的兩條軌道，即能量相差最少的軌道為最優(yōu)轉(zhuǎn)移方案。通過兩脈沖和三脈沖完成L1點(diǎn)Halo軌道至L2點(diǎn)Halo軌道的轉(zhuǎn)移，同時實(shí)現(xiàn)起始點(diǎn)、目標(biāo)點(diǎn)、中間脈沖點(diǎn)以及轉(zhuǎn)移時間的選取。根據(jù)圖8（c）可得，隨著算法迭代計算，能量指標(biāo)逐漸收斂至較低的值，從而實(shí)現(xiàn)設(shè)定代數(shù)內(nèi)的最優(yōu)能量轉(zhuǎn)移，可得到計算結(jié)果如表7所示。

表7 L1點(diǎn)Halo軌道族–L2點(diǎn)Halo軌道族的脈沖轉(zhuǎn)移結(jié)果Table 7 Transfer Impulse from L1 Halo family to L2 Halo family

2）L1點(diǎn)Halo軌道族到DRO軌道族的轉(zhuǎn)移

選擇振幅范圍為67 850～96 600 km繞L1點(diǎn)Halo軌道族和振幅范圍為86 450～285 600 km的DRO軌道族作為備選的起始軌道和目標(biāo)軌道，離散的軌道數(shù)據(jù)構(gòu)成三維矩陣，維度分別為：6 × 2 000 × 80和6 × 2 000 × 80。采用兩脈沖、三脈沖兩種方式進(jìn)行軌道轉(zhuǎn)移，指標(biāo)選擇為能量最優(yōu)。通過算法可求得轉(zhuǎn)移軌道分別如圖9（a）、（b）所示，指標(biāo)隨迭代的變化曲線如圖9（c）所示。

圖9 L1 Halo軌道族–DRO軌道族兩脈沖與三脈沖轉(zhuǎn)移Fig.9 Transfer from L1Halo family to DRO family using two impulses/three impulses

由圖9（a）、（b）所示，算法通過兩脈沖和三脈沖完成L1點(diǎn)Halo軌道至DRO軌道的轉(zhuǎn)移，同時實(shí)現(xiàn)起始點(diǎn)、目標(biāo)點(diǎn)、中間脈沖點(diǎn)以及轉(zhuǎn)移時間的選取。使用兩脈沖和三脈沖的差別很小，轉(zhuǎn)移軌道重合度較高。根據(jù)圖9（c）可得，隨著算法迭代計算，能量指標(biāo)逐漸收斂，可得到計算結(jié)果如表8所示。

表8 L1點(diǎn)Halo軌道族–DRO軌道族的脈沖轉(zhuǎn)移結(jié)果Table 8 Transfer impulses from L1 Halo family to DRO family

4.3 借助流形的轉(zhuǎn)移軌道

對于借助平動點(diǎn)軌道流形的脈沖轉(zhuǎn)移問題，仍然可將流形數(shù)據(jù)視為三體模型下的一類特殊軌道族。本文選取兩種軌道轉(zhuǎn)移場景進(jìn)行仿真：①繞L1點(diǎn)Halo軌道流形到繞L2點(diǎn)Halo軌道的轉(zhuǎn)移；② 繞L1點(diǎn)Halo軌道流形到DRO軌道的轉(zhuǎn)移。

1）L1點(diǎn)Halo軌道流形到L2點(diǎn)Halo軌道的轉(zhuǎn)移

選擇振幅為48 320 km繞L1點(diǎn)Halo軌道的不穩(wěn)定流形和振幅為52 900 km的繞L2點(diǎn)Halo軌道作為起始軌道和目標(biāo)軌道，離散的軌道數(shù)據(jù)構(gòu)成三維矩陣和二維矩陣，維度分別為：6 × 2 000 × 50和6 × 2 000。通過算法可求得兩脈沖轉(zhuǎn)移軌道分別如圖10（a）所示，指標(biāo)隨迭代的變化曲線如圖10（b）所示。

圖10 L1 Halo流形–L2 Halo兩脈沖轉(zhuǎn)移Fig.10 Transfer from manifolds of L1 Halo to L2 Halo using two impulses

2）L1點(diǎn)Halo軌道流形到DRO軌道的轉(zhuǎn)移

選擇振幅為48 320 km繞L1點(diǎn)Halo軌道的不穩(wěn)定流形和振幅為285 600 km的DRO作為起始軌道和目標(biāo)軌道，離散的軌道數(shù)據(jù)構(gòu)成三維矩陣和二維矩陣，維度分別為：6 × 2 000 × 50和6 × 2 000。采用兩脈沖進(jìn)行軌道轉(zhuǎn)移，指標(biāo)選擇為能量最優(yōu)。通過算法可求得轉(zhuǎn)移軌道分別如圖11（a）所示，指標(biāo)隨迭代的變化曲線如圖11（b）所示。

由于擾動施加在平動點(diǎn)軌道的不同位置可生成不同的流形，其離散數(shù)據(jù)可認(rèn)為與特殊軌道族類似。由圖10（a）、圖11（a）所示，算法通過兩脈沖完成L1點(diǎn)Halo軌道的流形至L2點(diǎn)Halo軌道和DRO軌道的轉(zhuǎn)移，同時實(shí)現(xiàn)起始點(diǎn)、目標(biāo)點(diǎn)、中間脈沖點(diǎn)以及轉(zhuǎn)移時間的選取。根據(jù)圖10（b）、圖11（b）可得，兩轉(zhuǎn)移軌道的指標(biāo)隨著算法迭代計算逐漸收斂，可得此算法對借助流形的轉(zhuǎn)移軌道設(shè)計仍然適用。脈沖優(yōu)化結(jié)果如表9所示。

表9 借助流形的脈沖轉(zhuǎn)移結(jié)果Table 9 Impulse transfer by means of manifolds

根據(jù)4.1～4.3節(jié)仿真結(jié)果可得本文所提出的分層重構(gòu)的建模方法能夠適應(yīng)各種不同特殊軌道之間的多脈沖轉(zhuǎn)移問題，也具有添加其它約束條件的潛力（如月球借力中的高度約束和近月點(diǎn)約束）。另外，如3.3節(jié)中對算法初值的分析，本文仿真中選擇所有脈沖矢量為 0矢量作為迭代初值，說明本算法無需對不同轉(zhuǎn)移任務(wù)設(shè)計不同的脈沖初值，算法普適性較高。

4.4 對比仿真

文獻(xiàn)[20]詳細(xì)分析了地月空間中NRHO與DRO之間的兩脈沖轉(zhuǎn)移問題，并根據(jù)軌道特性設(shè)計了轉(zhuǎn)移軌道的初值，文中關(guān)于初值的討論可總結(jié)如下。

1）遍歷搜索NRHO中的500個離散點(diǎn)作為起始點(diǎn)，DRO同樣被離散為數(shù)百個點(diǎn)，從軌跡上的點(diǎn)到DRO離散點(diǎn)的最小距離被視為到DRO的近似距離；

2）選取出發(fā)脈沖前后速度的比值：假設(shè)出發(fā)脈沖與出發(fā)點(diǎn)速度平行，且比值范圍在0.5～1.5；

3）兩參數(shù)確定的軌跡將在最大傳輸時間內(nèi)進(jìn)行積分，使用無量綱單位為20，大約相當(dāng)于87 d，分為200組進(jìn)行遍歷；

4）在搜索階段，軌跡可能多次接近DRO，從每組到DRO的最小距離中選取所有的最小距離，并將傳輸軌跡和DRO上的相應(yīng)點(diǎn)視為潛在的插入點(diǎn)。提取搜索參數(shù)、DRO上的潛在插入點(diǎn)和飛行時間作為初始猜測。然后，在每組中選擇估計傳輸成本最低的前100個初始猜測，以供后續(xù)優(yōu)化使用。

與文獻(xiàn)[20]中選擇相同振幅的NRHO與DRO分別作為起始軌道與目標(biāo)軌道，同時轉(zhuǎn)移時間也選擇相同條件，參數(shù)如表10所示。將參數(shù)代入本文所提出的雙層優(yōu)化算法框架中，可得到轉(zhuǎn)移軌道如圖12所示，所需脈沖與文獻(xiàn)數(shù)據(jù)的對比如表11所示。

表10 對比仿真相關(guān)參數(shù)表Table 10 Parameters of comparative simulation

表11 脈沖結(jié)果對比Table 11 Comparison of transfer impulses

可見，在軌道與轉(zhuǎn)移時間完全相同的仿真情況下，通過本文所提方法進(jìn)行計算得到的轉(zhuǎn)移脈沖可與文獻(xiàn)數(shù)據(jù)達(dá)到更好的最優(yōu)程度，且節(jié)省了約27 m/s的速度脈沖。這主要由于文獻(xiàn)中所提初值設(shè)計方法有一定的局限性，在進(jìn)行多個參數(shù)的選擇時容易陷入局部最優(yōu)情況，而本文通過對問題的分層建模求解，實(shí)現(xiàn)了更好的全局最優(yōu)性。同時無需經(jīng)過初值設(shè)計，也同樣適用于其它軌道；文獻(xiàn)中僅研究了兩脈沖轉(zhuǎn)移問題，所用初值設(shè)計方法也局限于兩脈沖，本文所提方法同樣適用于多脈沖轉(zhuǎn)移，實(shí)現(xiàn)了更好的普適性。

5 結(jié)論

本文基于離散的軌道數(shù)據(jù)，在地月圓型限制性三體模型下提出了一種多脈沖轉(zhuǎn)移軌道的設(shè)計方法?？紤]起點(diǎn)與終點(diǎn)不固定、轉(zhuǎn)移時間不固定、起始軌道與目標(biāo)軌道可變化等多種不確定因素，通過拆分變量將最小脈沖問題重構(gòu)為雙層優(yōu)化問題，設(shè)計了一種雙層優(yōu)化算法：上層兼顧了全局最優(yōu)性和計算效率，并通過引入離散軌道數(shù)據(jù)約束實(shí)現(xiàn)了對多種類型軌道的普遍適用；下層考慮始末狀態(tài)約束、時間約束、脈沖約束與脈沖點(diǎn)約束等實(shí)現(xiàn)最小脈沖問題的局部最優(yōu)，同時避免算法對初值的敏感性。

仿真結(jié)果表明，本文所提出的對離散軌道數(shù)據(jù)進(jìn)行始末狀態(tài)約束能夠適用于多種確定軌道、不同軌道族以及平動點(diǎn)軌道流形之間的相互轉(zhuǎn)移，同時提出的雙層優(yōu)化算法也能夠滿足多種約束并實(shí)現(xiàn)多脈沖下的近似最小脈沖轉(zhuǎn)移。對算法進(jìn)一步拓展，也可對更高精度星歷模型下的轉(zhuǎn)移軌道進(jìn)行計算。該研究成果為未來長期繞月探測、載人登月以及建設(shè)月球空間站等提供了有力的技術(shù)支持，通過優(yōu)化的軌道設(shè)計，將為這些任務(wù)的執(zhí)行提供更高的效率和可行性。未來的研究可以在本文基礎(chǔ)上進(jìn)一步探索，優(yōu)化和改進(jìn)該算法，以滿足更復(fù)雜的轉(zhuǎn)移需求，并拓展到其它三體問題和空間探索領(lǐng)域。