平動(dòng)點(diǎn)雙脈沖轉(zhuǎn)移軌道的快速計(jì)算方法

潘 迅，楊 瑞，泮斌峰，唐 碩

(1. 西北工業(yè)大學(xué)航天學(xué)院，西安710072；2. 航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安710072)

?

平動(dòng)點(diǎn)雙脈沖轉(zhuǎn)移軌道的快速計(jì)算方法

潘 迅1,2，楊 瑞1,2，泮斌峰1,2，唐 碩1

(1. 西北工業(yè)大學(xué)航天學(xué)院，西安710072；2. 航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安710072)

針對(duì)限制性三體問(wèn)題中的平動(dòng)點(diǎn)雙脈沖轉(zhuǎn)移，提出一種高效的計(jì)算方法。通過(guò)利用基于二維插值的數(shù)值流形近似方法對(duì)流形進(jìn)行近似計(jì)算，同時(shí)利用二體模型下的圓錐曲線近似流形拼接段，根據(jù)經(jīng)典軌道要素推導(dǎo)得到完成拼接所需的速度增量，避免在優(yōu)化過(guò)程中對(duì)流形的重復(fù)積分計(jì)算，以及在三體模型下對(duì)拼接段的迭代計(jì)算，從而顯著提高計(jì)算效率。然后推導(dǎo)得到三體問(wèn)題下的主矢量理論，可將其用于對(duì)優(yōu)化所得的雙脈沖轉(zhuǎn)移軌道進(jìn)行燃料最優(yōu)性的驗(yàn)證。最后，以航天器從近地圓軌道到地月系L1點(diǎn)的halo軌道的雙脈沖轉(zhuǎn)移為例進(jìn)行數(shù)值仿真，驗(yàn)證數(shù)值流形近似算法和二體模型近似脈沖的有效性，并表明該方法在優(yōu)化過(guò)程中具有高效性。

圓限制性三體問(wèn)題；平動(dòng)點(diǎn)；不變流形近似；主矢量理論；軌道優(yōu)化；地月轉(zhuǎn)移

0 引 言

隨著嫦娥計(jì)劃的進(jìn)行，我國(guó)的月球探測(cè)正在有條不紊的推進(jìn)過(guò)程中。相比于傳統(tǒng)的二體模型，圓限制性三體問(wèn)題模型能更精確地描述地月空間，且存在平動(dòng)點(diǎn)、周期軌道等眾多動(dòng)力學(xué)特性。月球附近存在兩個(gè)平動(dòng)點(diǎn)L1點(diǎn)和L2點(diǎn)，前者位于地月之間，可作為地月轉(zhuǎn)移的中轉(zhuǎn)站，后者位于地月連線的延長(zhǎng)線上，可對(duì)地球和月球背面的連續(xù)不間斷通信提供支持[1]，也可作為深空探測(cè)的基站，對(duì)這兩個(gè)平動(dòng)點(diǎn)的合理利用有利于后續(xù)深空探測(cè)的開(kāi)展。文獻(xiàn)[2-4]對(duì)平動(dòng)點(diǎn)附近的空間結(jié)構(gòu)，動(dòng)力學(xué)特性及其在深空探測(cè)中應(yīng)用進(jìn)行了詳細(xì)的介紹。 2010年ARTEMIS任務(wù)的P1和P2航天器分別進(jìn)入了地月L2點(diǎn)和地月L1點(diǎn)的擬周期軌道，成為世界上第一次地月平動(dòng)點(diǎn)任務(wù)[5]。我國(guó)的嫦娥五號(hào)試驗(yàn)器也于2014年對(duì)地月系L2點(diǎn)進(jìn)行了探測(cè)。

脈沖推進(jìn)作為目前的主要推進(jìn)方式，短期內(nèi)在實(shí)際任務(wù)的開(kāi)展和應(yīng)用中的地位不可動(dòng)搖。對(duì)于航天器從地球到平動(dòng)點(diǎn)的脈沖轉(zhuǎn)移軌道的優(yōu)化設(shè)計(jì)，目前有較為豐富的研究成果。文獻(xiàn)[6]利用三體Lambert問(wèn)題的迭代算法，對(duì)地球到平動(dòng)點(diǎn)周期軌道的雙脈沖轉(zhuǎn)移進(jìn)行了相關(guān)研究；文獻(xiàn)[7]結(jié)合網(wǎng)格搜索算法和遺傳算法對(duì)地月系平動(dòng)點(diǎn)轉(zhuǎn)移軌道的雙脈沖轉(zhuǎn)移進(jìn)行優(yōu)化；文獻(xiàn)[8]利用同構(gòu)映射將平面圓限制性三體問(wèn)題的維度從4維降為3維，然后運(yùn)用幾何方法對(duì)地球到平動(dòng)點(diǎn)的轉(zhuǎn)移軌道進(jìn)行優(yōu)化；文獻(xiàn)[9]基于擾動(dòng)流形和軌道角動(dòng)量對(duì)地月系平面限制性三體問(wèn)題的地球停泊軌道和地月L1點(diǎn)之間的多脈沖轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)；文獻(xiàn)[10]利用限制性三體問(wèn)題的動(dòng)力學(xué)系統(tǒng)理論研究了從地球到地月系L3點(diǎn)halo軌道的雙脈沖轉(zhuǎn)移；文獻(xiàn)[11]對(duì)地月系下雙脈沖轉(zhuǎn)移進(jìn)行了詳細(xì)研究，并對(duì)不同流形拼接點(diǎn)和不同halo軌道的轉(zhuǎn)移進(jìn)行了對(duì)比分析；此外，文獻(xiàn)[12]針對(duì)在日地系和地月系耦合雙三體模型下的脈沖轉(zhuǎn)移進(jìn)行了研究。

在三體模型下對(duì)平動(dòng)點(diǎn)轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)時(shí)，為節(jié)省能量，均利用了穩(wěn)定不變流形能無(wú)動(dòng)力趨向于周期軌道的特性。不變流形可通過(guò)對(duì)流形初始點(diǎn)進(jìn)行數(shù)值積分得到。在優(yōu)化過(guò)程中，需要對(duì)從近地停泊軌道進(jìn)入流形的位置進(jìn)行優(yōu)化，即選取合適的流形拼接點(diǎn)。在優(yōu)化過(guò)程中需要多次積分得到不同的流形狀態(tài)點(diǎn)作為拼接點(diǎn)，計(jì)算量較大，增加了優(yōu)化難度。文獻(xiàn)[13]提出了數(shù)值流形近似計(jì)算方法，先利用二維三次卷積插值得到流形近似值，然后根據(jù)能量對(duì)其進(jìn)行修正，得到更為精確的流形狀態(tài)點(diǎn)，避免了重復(fù)積分過(guò)程；文獻(xiàn)[14]也對(duì)該流形近似計(jì)算方法進(jìn)行了研究，并將其應(yīng)用于基于形狀法的低能轉(zhuǎn)移軌道。

在得到脈沖轉(zhuǎn)移軌道后，需要對(duì)該軌道的最優(yōu)性進(jìn)行驗(yàn)證。Lawden[15]在1963年提出二體動(dòng)力學(xué)的主矢量理論，可用于確定是否為最優(yōu)飛行軌道。文獻(xiàn)[16]對(duì)三體問(wèn)題下的月地返回軌道進(jìn)行優(yōu)化，并利用主矢量理論驗(yàn)證了其優(yōu)化結(jié)果的最優(yōu)性。

本文對(duì)地月系圓限制性三體問(wèn)題下的平動(dòng)點(diǎn)雙脈沖轉(zhuǎn)移軌道進(jìn)行研究。對(duì)于近地停泊軌道與不變流形的拼接段，以二體模型下的圓錐曲線進(jìn)行近似，得到優(yōu)化結(jié)果后再進(jìn)行修正，減小對(duì)拼接段脈沖的優(yōu)化難度及優(yōu)化時(shí)間，同時(shí)結(jié)合數(shù)值流形近似計(jì)算方法，縮短了優(yōu)化時(shí)間。與文獻(xiàn)中的研究相比，采用本文的方法能在極短時(shí)間內(nèi)得到優(yōu)化結(jié)果，顯著提高優(yōu)化效率，并經(jīng)主矢量理論驗(yàn)證，優(yōu)化結(jié)果為最優(yōu)脈沖轉(zhuǎn)移軌道。

1 限制性三體問(wèn)題

對(duì)于地月系圓限制性三體問(wèn)題，航天器在轉(zhuǎn)移過(guò)程中同時(shí)受到兩個(gè)主天體地球和月球的影響，其在旋轉(zhuǎn)坐標(biāo)系中的動(dòng)力學(xué)模型為

(1)

式中：μ表示主天體的質(zhì)量比，是三體系統(tǒng)的唯一參數(shù)，在地月系中取值為0.01215；r1、r2表示航天器與地球和月球之間的距離，分別為

(2)

平動(dòng)點(diǎn)、周期軌道和不變流形是限制性三體問(wèn)題的主要特征。不變流形分為趨近于平動(dòng)點(diǎn)周期軌道的穩(wěn)定流形和遠(yuǎn)離周期軌道的不穩(wěn)定流形，利用不變流形可進(jìn)行主天體與平動(dòng)點(diǎn)周期軌道之間的轉(zhuǎn)移軌道設(shè)計(jì)。在計(jì)算過(guò)程中，通過(guò)施加微小擾動(dòng)，使航天器偏離原周期軌道，得到流形積分初始點(diǎn)，對(duì)其進(jìn)行積分則可得到相應(yīng)不變流形，其初始點(diǎn)選取方法為

(3)

式中：x0為周期軌道上的點(diǎn)，Y為對(duì)應(yīng)x0的單位特征向量，上角標(biāo)s和u分別表示穩(wěn)定流形不穩(wěn)定流形，Xs(x0)和Xu(x0)分別為穩(wěn)定流形和不穩(wěn)定流形的積分初始點(diǎn)，ε表示小擾動(dòng)量，在本文中取值為2×10-6。對(duì)流形積分初始點(diǎn)進(jìn)行積分，即可得到相應(yīng)不變流形。L1點(diǎn)halo軌道的穩(wěn)定流形如圖1所示。

圖1 L1點(diǎn)halo軌道的穩(wěn)定流形Fig.1 The stable invariant manifolds of L1 halo orbit

2 不變流形數(shù)值近似計(jì)算

由于不變流形不能到達(dá)地球附近區(qū)域，在對(duì)轉(zhuǎn)移軌道進(jìn)行設(shè)計(jì)時(shí)，航天器從近地停泊軌道出發(fā)后需要選擇合適的位置進(jìn)入流形，即對(duì)流形拼接點(diǎn)進(jìn)行優(yōu)化。傳統(tǒng)的計(jì)算方法為通過(guò)積分得到流形點(diǎn)，由于在優(yōu)化過(guò)程中，對(duì)選取不同拼接點(diǎn)時(shí)的性能進(jìn)行分析對(duì)比，需要大量計(jì)算流形狀態(tài)點(diǎn)，而基于數(shù)值積分的計(jì)算方法需要對(duì)整條流形進(jìn)行積分才能得到某個(gè)狀態(tài)點(diǎn)，計(jì)算量較大，計(jì)算效率較為低下。Topputo針對(duì)該問(wèn)題提出了數(shù)值近似快速計(jì)算方法，通過(guò)二維插值，對(duì)流形狀態(tài)點(diǎn)進(jìn)行插值近似，在優(yōu)化過(guò)程中不需要對(duì)流形重復(fù)積分，從而在滿足計(jì)算精度的前提下，大大提高了計(jì)算效率[13]。

在積分得到流形狀態(tài)點(diǎn)的過(guò)程中，只與參數(shù)τ和t有關(guān)，其中τ為對(duì)周期軌道初始點(diǎn)進(jìn)行積分，并得到式中x0的積分時(shí)間，t為對(duì)穩(wěn)定流形初始點(diǎn)Xs(x0)的積分時(shí)間，因此通過(guò)二維插值則可得到流形狀態(tài)點(diǎn)。首先通過(guò)積分得到流形插值的原始數(shù)據(jù)，對(duì)積分流形的網(wǎng)格劃分為

(4)

式中：T1為平動(dòng)點(diǎn)周期軌道的周期，T2為流形的最大積分時(shí)間，對(duì)于穩(wěn)定流形有T2<0，對(duì)不穩(wěn)定流形則有T2>0，N和M為離散點(diǎn)數(shù)量。根據(jù)τi和tj則可積分得到相應(yīng)的穩(wěn)定流形點(diǎn)xs,ij(τ,t)。

得到插值的原始數(shù)據(jù)xs,ij后，對(duì)任意給定的(τ,t)進(jìn)行計(jì)算的插值公式為

(5)

式中：ci+n,j+m為插值樣點(diǎn)，h1=T1/(N-1)和h2=T2/(M-1)為離散步長(zhǎng)，u為插值內(nèi)核。

對(duì)于插值內(nèi)核，設(shè)u(0)=1，u(n)=0，對(duì)其他值的計(jì)算式為

(6)

插值樣點(diǎn)c為六維向量，分別對(duì)應(yīng)流形狀態(tài)點(diǎn)的三維位置矢量和三維速度矢量，計(jì)算公式為

(7)

(8)

(9)

3 二體模型脈沖近似

在得到流形拼接點(diǎn)后，需要對(duì)近地停泊軌道和流形進(jìn)行拼接。在限制性三體問(wèn)題中，三體Lambert轉(zhuǎn)移需要進(jìn)行多次迭代計(jì)算，計(jì)算過(guò)程較為復(fù)雜。由于月球引力遠(yuǎn)小于地球引力，因此在計(jì)算拼接段時(shí)，可先利用地球-航天器二體模型進(jìn)行近似計(jì)算，然后在三體模型下進(jìn)行修正，從而降低計(jì)算難度和減小計(jì)算量。

在二體模型下，根據(jù)經(jīng)典軌道六要素，對(duì)于拼接段圓錐曲線，有

(10)

(11)

rp=a(1-e)

(12)

(13)

然后可計(jì)算得到航天器在拼接點(diǎn)所需施加的速度脈沖為

ΔvMI=(1-k)vt

(14)

到達(dá)近地停泊軌道后需要施加的速度脈沖為

(15)

可得到二體模型近似下所需的速度增量為

Δvapp=ΔvE+ΔvMI

(16)

Φ(tp+Δt,xpr+Δxr)=Φ(tp,xpr)+

(17)

(18)

兩個(gè)修正量對(duì)應(yīng)兩個(gè)約束方程，因此可求解得到Δt和p，從而計(jì)算得到三體模型下的拼接段以及完成拼接所需的速度增量ΔvEm和ΔvMIm。

4 主矢量理論

在1963年提出的主矢量理論中，Lawden給出的最優(yōu)脈沖的必要條件為[15]：

4)對(duì)于所有內(nèi)點(diǎn)脈沖(初始時(shí)刻脈沖和終端時(shí)候脈沖除外)成立。

本文對(duì)限制性三體問(wèn)題下的主矢量進(jìn)行推導(dǎo)。在三體模型中，推力加速度近似為脈沖表示

(19)

(20)

式中：ti為施加脈沖時(shí)刻。

根據(jù)最優(yōu)控制理論，構(gòu)造哈密爾頓方程為

(21)

(22)

(24)

得到脈沖優(yōu)化結(jié)果后，根據(jù)所施加的脈沖計(jì)算得到λv(t0)和λv(tf)，再以λr(t0)為迭代變量，根據(jù)邊界條件λv(tf)求解出λr(t0)，然后積分得到主矢量λv在時(shí)間區(qū)間[t0,tf]上的值，判斷是否滿足所有Lawden必要條件，從而可判斷優(yōu)化得到雙脈沖轉(zhuǎn)移軌道是否為最優(yōu)脈沖轉(zhuǎn)移軌道。

5 數(shù)值仿真

考慮航天器從近地軌道出發(fā)，通過(guò)施加第一次脈沖進(jìn)入拼接段，在流形拼接點(diǎn)處施加第二次脈沖進(jìn)行L1點(diǎn)halo軌道的穩(wěn)定流形，然后沿流形無(wú)動(dòng)力的趨近于halo軌道。地球平均半徑為6378.137km，近地停泊軌道為高度300km的圓軌道，L1點(diǎn)halo軌道的幅值A(chǔ)z為5000km。本算例中的計(jì)算在CPU為i7-4700MQ，主頻為2.40GHz，內(nèi)存為8G的計(jì)算機(jī)上進(jìn)行，計(jì)算平臺(tái)為MATLAB2015a。

考慮到實(shí)際任務(wù)中的飛行時(shí)間限制，流形積分時(shí)間上限為50天，同時(shí)為保證數(shù)值近似流形的精度，在halo軌道上均勻選取200個(gè)點(diǎn)，在每條流形上選取2001個(gè)點(diǎn)，作為流形插值數(shù)據(jù)庫(kù)，即式中N=200，M=2001。為驗(yàn)證數(shù)值近似流形的精度，將其與積分得到的流形點(diǎn)進(jìn)行誤差對(duì)比分析。定義流形誤差為

(25)

誤差結(jié)果如表1所示，誤差的對(duì)數(shù)lg(ξ)的分布情況如圖2所示，可知平均精度為1.6529×10-6，相當(dāng)于是0.6353km，在地月轉(zhuǎn)移軌道設(shè)計(jì)過(guò)程中，該精度已屬于較高精度，因此認(rèn)為數(shù)值近似流形的精度足夠滿足轉(zhuǎn)移軌道初步設(shè)計(jì)的要求。

表1 數(shù)值近似流形的誤差情況Table 1 The error of the approximation invariant manifolds

選取τ=0時(shí)的流形，在該流形積分時(shí)間為[-50, -10]天的范圍內(nèi)，根據(jù)時(shí)間均勻選擇160個(gè)點(diǎn)，對(duì)二體模型近似脈沖和三體模型修正后脈沖進(jìn)行對(duì)比，如圖3所示，其中ΔvEm和ΔvMIm為三體模型修正后所需的脈沖，對(duì)比分析可知，二體模型下的近似脈沖能很好地估計(jì)完成拼接段轉(zhuǎn)移所需的速度增量。

圖2 數(shù)值近似流形的誤差分析lg(ξ)Fig.2 The error pattern lg(ξ) of approximation invariant manifolds over the evaluation grid

對(duì)于拼接點(diǎn)的確定，利用MATLAB自帶函數(shù)fmincon進(jìn)行優(yōu)化，優(yōu)化變量為[τ,t]。在優(yōu)化過(guò)程中，只利用二體模型脈沖近似，在得到優(yōu)化結(jié)果后再在三體模型下進(jìn)行修正，修正精度要求為1km。由于fmincon是局部?jī)?yōu)化函數(shù)，選取不同的初始值會(huì)得到不同的優(yōu)化結(jié)果。多次選取不同的初始值，可得到4組優(yōu)化結(jié)果，與圖3中的單條流形的4個(gè)極小值點(diǎn)相對(duì)應(yīng)。優(yōu)化結(jié)果數(shù)據(jù)如表2所示，其中優(yōu)化時(shí)間為只考慮fmincon函數(shù)優(yōu)化的時(shí)間，4組結(jié)果所需優(yōu)化時(shí)間均小于0.1s。這4組結(jié)果對(duì)應(yīng)的轉(zhuǎn)移軌道分別如圖4～7所示，其中左圖為旋轉(zhuǎn)坐標(biāo)系下的轉(zhuǎn)移軌道，右圖為地心慣性系下的轉(zhuǎn)移軌道。從圖4～7的(a)圖可以看出，優(yōu)化得到的4個(gè)拼接點(diǎn)均為流形的遠(yuǎn)地點(diǎn)，這是由于在遠(yuǎn)地點(diǎn)受到的地球引力較小，克服引力所消耗的能量較小，脈沖的利用率較大，因此完成轉(zhuǎn)移所需的總速度增加較??；從右圖可以看出，拼接點(diǎn)在流形上的位置相差約一個(gè)周期，結(jié)合表2中航天器沿流形的運(yùn)動(dòng)時(shí)間可知，流形繞地球運(yùn)動(dòng)一周所需時(shí)間約為10天。

圖3 二體模型近似脈沖與三體模型修正后脈沖對(duì)比Fig.3 Comparison of the two-body approximate impulse and the impulse modified in CRTBP

優(yōu)化結(jié)果ΔvMImΔvLEOmΔvtot轉(zhuǎn)移段時(shí)間/d流形段時(shí)間/d優(yōu)化時(shí)間/sA0．61633．06793．68423．858417．60340．0599B0．56103．06333．62443．447028．00000．0614C0．56793．06333．63123．416638．00450．0495D0．59073．06363．65443．433248．03220．0835Larsen等[7]0．663．063．723．018．8———Pontani等[8]0．5613．0163．57720．183．21062．8?

圖4 雙脈沖轉(zhuǎn)移優(yōu)化結(jié)果AFig.4 The results of two-impulse transfers: A

Larsen等[7]和Pontani等[8]均對(duì)航天器從地球到L1點(diǎn)周期軌道的雙脈沖轉(zhuǎn)移進(jìn)行了研究，但僅考慮xy平面內(nèi)的運(yùn)動(dòng)，目標(biāo)軌道為L(zhǎng)1點(diǎn)Lyapunov軌道。前者的近地停泊軌道高度為300km，后者的近地停泊軌道高度為500km，優(yōu)化結(jié)果在表2中描述。文獻(xiàn)[7]所需的速度增量小于本文優(yōu)化結(jié)果主要是因?yàn)槠浣赝２窜壍赖母叨容^高，克服地球引力所需能量較??；但對(duì)于相同高度的近地停泊軌道，本文的B組結(jié)果不論從轉(zhuǎn)移時(shí)間還是所需速度增量，都要優(yōu)于文獻(xiàn)[8]中的結(jié)果。文獻(xiàn)[11]對(duì)地月系下的雙脈沖轉(zhuǎn)移進(jìn)行了詳細(xì)研究，當(dāng)目標(biāo)halo軌道的橫坐標(biāo)為319052km時(shí)，對(duì)拼接點(diǎn)進(jìn)行優(yōu)化時(shí)，得到多個(gè)極小值點(diǎn)，該拼接點(diǎn)均為流形遠(yuǎn)地點(diǎn)，所需最小速度增量約為3.62km/s，與本文的優(yōu)化結(jié)果類似。更為重要的是，利用文獻(xiàn)[7]中的網(wǎng)格法進(jìn)行優(yōu)化時(shí)需耗費(fèi)較長(zhǎng)時(shí)間，文本作者復(fù)現(xiàn)文獻(xiàn)[8]中的結(jié)果需用時(shí)1062.8s，利用[11]中計(jì)算單條轉(zhuǎn)移軌道時(shí)也需要對(duì)施加脈沖進(jìn)行打靶求解，而利用數(shù)值流形近似和二體模型脈沖近似方法，積分得到流形插值數(shù)據(jù)需耗時(shí)26.7s，然后利用fmincon函數(shù)求得上文中一個(gè)極小值解所需優(yōu)化時(shí)間不到0.1s，因此本文的快速計(jì)算方法顯著提高了計(jì)算效率。

圖5 雙脈沖轉(zhuǎn)移優(yōu)化結(jié)果BFig.5 The results of two-impulse transfers: B

圖6 雙脈沖轉(zhuǎn)移優(yōu)化結(jié)果CFig.6 The results of two-impulse transfers: C

圖7 雙脈沖轉(zhuǎn)移優(yōu)化結(jié)果DFig.7 The results of two-impulse transfers: D

圖8 B組優(yōu)化結(jié)果的主矢量曲線Fig.8 The primer vector curve of result B

6 結(jié) 論

本文對(duì)限制性三體問(wèn)題的平動(dòng)點(diǎn)雙脈沖轉(zhuǎn)移軌道的快速計(jì)算方法進(jìn)行研究。利用數(shù)值流形近似方法對(duì)流形進(jìn)行插值并修正，從而避免了在優(yōu)化過(guò)程中重復(fù)積分計(jì)算流形。對(duì)于確定的流形拼接點(diǎn)，以二體模型下的圓錐曲線作為拼接段的近似，根據(jù)經(jīng)典軌道要素推導(dǎo)得到完成拼接所需的速度增量，再在三體模型下進(jìn)行修正。最后利用推導(dǎo)得到的三體模型下的主矢量理論對(duì)優(yōu)化結(jié)果是否滿足燃料最優(yōu)進(jìn)行驗(yàn)證。通過(guò)數(shù)值仿真表明，結(jié)合這兩種方法能在較短時(shí)間內(nèi)得到優(yōu)化結(jié)果，計(jì)算效率得到顯著提升；并選取其中一條轉(zhuǎn)移軌道，利用主矢量理論驗(yàn)證該脈沖轉(zhuǎn)移軌道的燃料最優(yōu)性，顯示利用該方法在一定程度上能得到燃料最優(yōu)轉(zhuǎn)移軌道，從而表明該方法的合理性和高效性。

[1] Farquhar R, Kamel A. Quasi-periodic orbits about the translunar libration point [J]. Celestial Mechanics and Dynamical Astronomy, 1973, 7(4): 458-473.

[2] 徐明. 平動(dòng)點(diǎn)軌道的動(dòng)力學(xué)與控制研究綜述[J]. 宇航學(xué)報(bào), 2009, 30(4):1299-1313. [Xu Ming. Overview of orbital dynamics and control for libration point Orbits [J]. Journal of Astronautics, 2009, 30(4):1299 -1313.]

[3] 徐明, 徐世杰. 地-月系平動(dòng)點(diǎn)及halo軌道的應(yīng)用研究[J]. 宇航學(xué)報(bào), 2006, 27(4):695-699. [Xu Ming, Xu Shi-jie. The application of libration points and halo orbits in the Earth-Moon system to space mission design [J]. Journal of Astronautics, 2006, 27(4):695 -699.]

[4] 侯錫云, 劉林. 共線平動(dòng)點(diǎn)的動(dòng)力學(xué)特征及其在深空探測(cè)中的應(yīng)用[J]. 宇航學(xué)報(bào), 2008, 29(3):736-747. [Hou Xi-yun, Liu Lin. The dynamics and applications of the collinear libration points in deep space exploration [J]. Journal of Astronautics, 2008, 29(3):736 -747.]

[5] Cosgrove D, Frey S, Bester M, et al. Navigating THEMIS to the ARTEMIS low-energy lunar transfer trajectory[C]. SpaceOps 2010 Conference, Alabama, USA. April 25-30, 2010.

[6] 連一君. 基于三體平動(dòng)點(diǎn)的低能轉(zhuǎn)移軌道設(shè)計(jì)研究[D]. 長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué), 2008. [Lian Yi-jun. Research on low-cost transfer trajectories design based on three-body libration points [D]. Changsha: National University of Defense Technology, 2008. ]

[7] Larsen A, Anthony W, Critz T, et al. Optimal transfers with guidance to the Earth-Moon L1and L3libration points using invariant manifolds: a preliminary study [C]. AIAA/AAS Astrodynamics Specialist Conference, Minnesota, USA. August 13-16, 2012.

[8] Pontani M, Teofilatto P. Low-energy Earth-Moon transfers involving manifolds through isomorphic mapping [J]. Acta Astronautica, 2013, 91(10):96-106.

[9] 張漢清, 李言俊. 地月三體問(wèn)題下L1-地球低能轉(zhuǎn)移軌道設(shè)計(jì)[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2011, 43(5):84-88. [Zhang Han-qing, Li Yan-jun. Design of L1-Earth low energy transfer trajectory in Earth-Moon three-body problem [J]. Journal of Harbin Institute of Technology, 2011, 43(5):84-88.]

[10] Davis K, Born G, Butcher E. Transfers to Earth-Moon L3 halo orbits [J]. Acta Astronautica, 2013, 88(3):116 -128.

[11] Parker J S, Born G H. Direct lunar halo orbit transfers [J]. The Journal of the Astronautical Sciences, 2008, 56(4):441-476.

[12] Parker J S, Anderson R L, Peterson A. Surveying ballistic transfers to low lunar orbit [J]. Journal of Guidance Control & Dynamics, 2013, 36(5):1501 -1511.

[13] Topputo F. Fast numerical approximation of invariant manifolds in the circular restricted three-body problem [J]. Communications in Nonlinear Science & Numerical Simulation, 2016, 32:89-98.

[14] 張仁勇. 深空探測(cè)小推力低能轉(zhuǎn)移軌道設(shè)計(jì)方法研究[D]. 西安：西北工業(yè)大學(xué), 2015. [Zhang Ren-yong. Research on low-thrust low-energy transfer trajectory design for deep space exploration [D]. Xi’an: Northwestern Polytechnical University, 2015. ]

[15] Lawden D F. Optimal trajectories for space navigation [J]. Mathematical Gazette, 1964, 48(366):478.

[16] Shen H, Casalino L. Indirect optimization of three-dimensional multiple-impulse Moon-to-Earth transfers [J]. The Journal of the Astronautical Sciences, 2014, 61(3):1-20.

通信地址：陜西省西安市碑林區(qū)友誼西路127號(hào)西北工業(yè)大學(xué)251信箱(710072)

電話：15191433469

E-mail：xpan2012@gmail.com

An Efficient Calculation Method for Two-Impulse Transfers to Libration Point

PAN Xun1,2, YANG Rui1,2, PAN Bin-feng1,2,TANG Shuo1

(1. School of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China;2. National Key Laboratory of Aerospace Flight Dynamics, Xi’an 710072, China)

An efficient calculation method is proposed for the optimization of the two-impulse transfers to the librations in circular restricted three-body problem. With the invariant manifolds approximation based on a two-dimensional interpolation, the computation of the great number of manifold insertions is much more efficient. While the conic curve in a two-body problem used as an approximation of the transfer segment connecting the low Earth orbit and the stable manifold, the maneuvers can be deduced with the classical orbital elements, thus the iterative computation of the maneuvers is avoided in the optimization process. Then the primer vector of CRTBP is deduced to verify the optimality of the two-impulse transfer. At last, a numerical simulation of the two-impulse transfer from a low Earth orbit to aL1point halo orbit is presented, verifying the validity of the two-body impulse approximation and the manifold approximation, and demonstrating the efficiency of the method.

Circular restricted three-body problem; Libration point; Invariant manifolds approximation; Primer vector theory; Trajectory optimization; Earth-Moon transfers

2017-01-09;

2017-04-21

國(guó)家自然科學(xué)基金(11672234)

V448.2

A

1000-1328(2017)06-0574-09

10.3873/j.issn.1000-1328.2017.06.003

潘 迅(1990-)，男，博士生，主要研究方向?yàn)轱w行動(dòng)力學(xué)與控制，深空探測(cè)軌道優(yōu)化設(shè)計(jì)。