賈建華 王 琪

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

引力輔助技術(shù) (gravity-assist)，又稱(chēng)借力飛行或引力甩擺，是一種航天任務(wù)中經(jīng)常用來(lái)減少燃料消耗的方法.目前已經(jīng)有多個(gè)任務(wù)利用引力輔助技術(shù)設(shè)計(jì)了行星際探測(cè)軌道，完成了太陽(yáng)系八大行星的探測(cè)任務(wù)［1］.

引力輔助技術(shù)起源于18世紀(jì)，Laplace、leverrier、Tisserand用此來(lái)解釋行星對(duì)彗星軌道的攝動(dòng)，Lawden首先提出可將引力輔助技術(shù)應(yīng)用于飛行器［2］.目前對(duì)引力輔助技術(shù)的研究主要集中在應(yīng)用方面，如引力輔助序列的選?。?－4］、小推力引力輔助的研究［5－6］、氣動(dòng)引力輔助的研究［7－8］，而對(duì)引力輔助機(jī)理本身的探討和研究較少.1988年文獻(xiàn)［9］對(duì)二維引力輔助機(jī)理進(jìn)行了分析，給出了前側(cè)飛越使飛行器能量減少、后側(cè)飛越使飛行器能量增加等重要結(jié)論，揭開(kāi)了人類(lèi)對(duì)引力輔助機(jī)理認(rèn)識(shí)的新篇章，但討論僅限于平面圓型限制性三體問(wèn)題.1997年文獻(xiàn) ［10］采用平面橢圓型限制性三體模型對(duì)引力輔助機(jī)理進(jìn)行了深入研究，結(jié)果表明，偏心率和真近點(diǎn)角對(duì)引力輔助軌道具有重要影響.文獻(xiàn) ［11－12］也對(duì)基于平面限制性三體問(wèn)題模型的引力輔助機(jī)理進(jìn)行了一些討論.以上的研究都只限于二維引力輔助，飛行器運(yùn)動(dòng)被限制在借力天體的運(yùn)動(dòng)平面內(nèi)，這顯然不能滿(mǎn)足實(shí)際任務(wù)的需要.目前，對(duì)三維引力輔助機(jī)理的研究尚不多見(jiàn)，文獻(xiàn)［13］初步討論了面外角對(duì)引力輔助軌道的影響，但是，對(duì)面外角的研究較為簡(jiǎn)略，而且對(duì)飛行器近心點(diǎn)速度有一定限制，因此需要對(duì)三維引力輔助機(jī)理進(jìn)行進(jìn)一步的研究.

為此，本文在前人研究工作的基礎(chǔ)上，通過(guò)引入兩個(gè)參數(shù)，將對(duì)引力輔助機(jī)理的探討從平面橢圓型限制性三體模型拓展到三維橢圓型限制性三體模型，重點(diǎn)討論了這兩個(gè)參數(shù)對(duì)引力輔助前后的飛行器軌道類(lèi)型、軌道能量和軌道傾角的影響.最后，以地月系統(tǒng)為例，針對(duì)不同任務(wù)要求，給出了引力輔助參數(shù)的選擇區(qū)域，并對(duì)引力輔助參數(shù)進(jìn)行了優(yōu)化.

1 問(wèn)題描述

在三維橢圓型限制性三體模型下，三維引力輔助可通過(guò)6個(gè)自由參數(shù)來(lái)描述 (如圖1所示):①Vp:飛越借力天體時(shí)近心點(diǎn)處的速度;②rp:飛越借力天體時(shí)近心點(diǎn)距離;③ν:借力天體的真近點(diǎn)角;④ψ:rp在x-y平面的投影與M1，M2連線的夾角;⑤α:rp與x-y平面的夾角;⑥β:VP與x-y平面的夾角.

圖1 三維橢圓型限制性三體模型引力輔助示意圖

2 三維橢圓型限制性三體模型

假設(shè)兩個(gè)主天體分別為M1和M2，繞其公共質(zhì)心做橢圓運(yùn)動(dòng)，第三體M3(本文中為飛行器)的質(zhì)量遠(yuǎn)遠(yuǎn)小于兩個(gè)主天體的質(zhì)量，在兩個(gè)主天體的作用下運(yùn)動(dòng)，而不影響兩個(gè)主天體的運(yùn)動(dòng).M3的運(yùn)動(dòng)軌跡不限于兩個(gè)主天體的運(yùn)動(dòng)平面內(nèi)，可以是三維的.

為了討論和計(jì)算上的方便，采用無(wú)量綱形式，定義:以M1和M2相對(duì)軌道的半長(zhǎng)軸為單位長(zhǎng)度1;平均角速度為1;單位質(zhì)量為m1+m2，其中，m1，m2分別是 M1和 M2的真實(shí)質(zhì)量，這樣，M2的質(zhì)量可以用μ=m2/(m1+m2)表示，M1的質(zhì)量為 (1－μ);定義單位時(shí)間，使得主天體的運(yùn)動(dòng)周期為2π;M1-M2系統(tǒng)的引力常數(shù) (萬(wàn)有引力常數(shù)G=6.67×10－11m3·kg－1·s－2與系統(tǒng)總質(zhì)量的乘積)定義為1.

用于分析此模型的參照系有很多，本文采用慣性坐標(biāo)系.在慣性系中，坐標(biāo)原點(diǎn)為M1和M2的公共質(zhì)心;橫軸x是初始時(shí)刻連接M1和M2的連線;縱軸y過(guò)原點(diǎn)垂直于x軸并且在兩個(gè)主天體的運(yùn)動(dòng)平面內(nèi);z軸符合右手定則.在此參考系中，M1和M2的位置分別為

其中，r=(1－e2)/(1+ecosν)，為M1和M2間的距離;ν為M2的真近點(diǎn)角;e為橢圓軌道的偏心率.

M3的動(dòng)力學(xué)方程可表示為

其中，r1，r2分別為M3到M1和M2的距離.

此外，時(shí)間和真近點(diǎn)角的關(guān)系為

其中，p=a(1－e2)，為橢圓的半正交弦.

可以求解出此慣性系下，M3的能量EI和角動(dòng)量CI分別為

其中，RI是M3的矢徑;VI是M3的速度向量.則M3的軌道傾角i=arccos(CIz/|CI|)，其中，CIz為CI的z軸分量.

3 數(shù)值算法

由于橢圓型限制性三體問(wèn)題目前無(wú)法得到軌道變化的解析表達(dá)式，所以只能通過(guò)數(shù)值積分進(jìn)行分析.數(shù)值算法流程如下:

1)給參數(shù)Vp，rp，ν，ψ，α，β，e，μ賦初值.

2)計(jì)算M3的初始位置和初始速度:

3)在初始條件下，對(duì)動(dòng)力學(xué)方程進(jìn)行向前積分，直到M3與M2的距離大于指定距離d時(shí)停止.計(jì)算出引力輔助之后飛行器的能量E+和角動(dòng)量C+.

4)在初始條件下，對(duì)動(dòng)力學(xué)方程進(jìn)行向后積分，直到M3與M2的距離大于指定距離d時(shí)停止.計(jì)算出引力輔助之前飛行器的能量E－和角動(dòng)量C－.

本文所有數(shù)值模擬，指定距離d都為0.5，通過(guò)變化參數(shù)的初值，可以研究參數(shù)對(duì)引力輔助軌道的影響.

4 數(shù)值仿真

根據(jù)軌道動(dòng)力學(xué)的基本定義:軌道能量為“正”表示雙曲線軌道， “負(fù)”表示橢圓軌道;角動(dòng)量z軸分量為“正”表示順行軌道，“負(fù)”表示逆行軌道.根據(jù)得到的E－，C－，E+，C+，可以將引力輔助軌道劃分為16種類(lèi)型，分別用前16個(gè)英文字母表示，具體見(jiàn)表1.

表1 英文字母與軌道類(lèi)型的對(duì)應(yīng)關(guān)系

關(guān)于引力輔助參數(shù)Vp，rp，ν，ψ對(duì)飛行器軌道的影響，文獻(xiàn) ［9－13］已經(jīng)進(jìn)行了討論，在此不再贅述.本文將重點(diǎn)討論參數(shù)α和β對(duì)引力輔助前后的飛行器軌道類(lèi)型、軌道能量和軌道傾角的影響.

4.1 參數(shù)α的分析

4.1.1 α對(duì)引力輔助前后軌道類(lèi)型的影響

本文對(duì)大量不同的α取值進(jìn)行了仿真，特選取了其中較有代表性的情況，見(jiàn)圖2.圖2所示是 α 分別為0°，60°，90°，120°時(shí)的 ψ-Vp圖，圖中字母代表對(duì)應(yīng)的軌道類(lèi)型.其他參數(shù)取值分別為:μ=0.012 1，rp=0.004 95，e=0，ν=0，β=0.此組參數(shù)下，0＜ψ＜π圖形與π＜ψ＜2π圖形關(guān)于ψ=π鏡像 (鏡像規(guī)律:B和E，C和I，D和M，G和J，H和N，L和O互為鏡像，A，F(xiàn)，K，P不變)，因此只給出ψ為π～2π時(shí)情況;α與360°－α的圖形完全相同，因此只給出α為0°～180°情況.

圖2 α對(duì)引力輔助軌道類(lèi)型的影響

當(dāng)α=β=0°時(shí)，三維橢圓型限制性三體問(wèn)題就退化為平面橢圓型限制性三體問(wèn)題，三維引力輔助退化為二維引力輔助，本文結(jié)果與文獻(xiàn)［10］的結(jié)果完全相同，說(shuō)明本文對(duì)三維引力輔助機(jī)理的研究可以很好地退化到二維情況.

當(dāng)α≠0°時(shí)，引力輔助為三維情況.從圖2可以看到:隨著α的變化，各種軌道類(lèi)型區(qū)域均發(fā)生變化.現(xiàn)將其變化規(guī)律總結(jié)如下:

1)橢圓軌道區(qū)域 (引力輔助前后都為橢圓軌道，包含類(lèi)型A，B，E，F(xiàn))總是出現(xiàn)在左下方.α從0°～90°，整體區(qū)域增大，A，F(xiàn)增大，B減小;α=90°時(shí)B消失;α從90°～180°，整體區(qū)域減小，A，F(xiàn)減小，E(B的鏡像類(lèi)型)增大.

2)雙曲軌道區(qū)域 (引力輔助前后都為雙曲軌道，包含類(lèi)型K，L，O，P)總是出現(xiàn)在右上方.α從0°～90°，整體區(qū)域增大，K，P增大，L減小;α=90°時(shí)L消失;α從90°～180°，整體區(qū)域減小，K、P減小，O(L的鏡像類(lèi)型)增大.

3)逃逸軌道區(qū)域 (引力輔助前后分別為橢圓軌道和雙曲軌道，包含類(lèi)型I，J，M，N)總是夾在橢圓軌道區(qū)域和雙曲軌道區(qū)域之間.α從0°～90°，整體區(qū)域縮小，I，J，N均減小;α=90°時(shí)均消失;α為90°～180°時(shí)，不出現(xiàn).

4)俘獲軌道區(qū)域 (引力輔助前后分別為雙曲軌道和橢圓軌道，包含類(lèi)型C，D，G，H)總是夾在橢圓軌道區(qū)域和雙曲軌道區(qū)域之間.α為0°～90°時(shí)，不出現(xiàn);α從90°～180°，整體區(qū)域增大，C，G，H均增大.

5)逆行軌道區(qū)域 (引力輔助前后都為逆行軌道，包含類(lèi)型F，H，N，P)總是出現(xiàn)在左上方.α從0°～90°，整體區(qū)域略有增大，F(xiàn)，P增大，N減小;α=90°時(shí)N消失;α從90°～180°，整體區(qū)域略有減小，F(xiàn)，P減小，H(N的鏡像類(lèi)型)增大.

6)順行軌道區(qū)域 (引力輔助前后都為順行軌道，包含類(lèi)型A，C，I，K)總是出現(xiàn)在右半部和底部 ，一直占據(jù)大部分.α從0°～90°，整體區(qū)域增大，A，K增大，I減小;α=90°時(shí)，I消失;α從90°～180°，整體區(qū)域減小，A，K減小，C(I的鏡像類(lèi)型)增大.

7)逆行變順行軌道區(qū)域 (引力輔助前后分別為逆行軌道和順行軌道，包含類(lèi)型B，D，J，L)總是夾在逆行軌道區(qū)域和順行軌道區(qū)域之間.α從0°～90°，整體區(qū)域縮小，B，J，L均減小;α=90°時(shí)均消失;α為90～180°時(shí)，不出現(xiàn).

8)順行變逆行軌道區(qū)域 (引力輔助前后分別為順行軌道和逆行軌道，包含類(lèi)型:E，G，M，O)總是夾在逆行軌道區(qū)域和順行軌道區(qū)域之間.α為 0°～90°時(shí)，不出現(xiàn);α從 90°～180°，整體區(qū)域增大，E，G，O均增大.

4.1.2 α對(duì)飛行器能量的影響

圖3所示為α對(duì)飛行器能量變化的影響.其中參數(shù)設(shè)置如下:μ=0.012 1，rp=0.004 95，e=0，ν=0，β=0，ψ =20°，Vp=2.6.

圖3 α對(duì)軌道能量變化的影響

可以看到，曲線呈鐘形，關(guān)于 α=180°對(duì)稱(chēng).α=0°或180°時(shí)，引力輔助前后飛行器軌道能量改變最大.當(dāng)α=90°或270°時(shí)，Δe=0.本文對(duì)其他rp值進(jìn)行了模擬，所得曲線仍呈鐘形.

圖4 α對(duì)軌道傾角變化的影響

4.1.3 α對(duì)飛行器軌道傾角的影響

圖4所示為α對(duì)飛行器軌道傾角變化的影響.參數(shù)設(shè)置與圖3相同.

可以看到，圖形關(guān)于α=180°對(duì)稱(chēng).α=0°～180°曲線形狀類(lèi)似于正弦曲線.α =0°，90°，180°，270°時(shí)，Δi=0.本文對(duì)其他 rp值進(jìn)行了模擬，所得曲線形狀類(lèi)似.

4.2 參數(shù)β的分析

4.2.1 β對(duì)引力輔助前后軌道類(lèi)型的影響

本文對(duì)大量的不同β取值進(jìn)行了仿真，特選取了其中較有代表性的情況，見(jiàn)圖5.圖5所示是β分別為 0°，67.5°，75°，90°，105°，112.5°，180°時(shí)的ψ-Vp圖，其他參數(shù)取值分別為:μ=0.0121，rp=0.00495，e=0，ν=0，α=0.此組參數(shù)下，0＜ψ＜π圖形與π＜ψ＜2π圖形關(guān)于ψ=π鏡像 (鏡像規(guī)則與4.1.1節(jié)中相同)，因此只給出ψ為π～2π時(shí)情況;β與360°－β的圖形完全相同，因此只給出β為0°～180°情況.

從圖5可以看到:隨著β的變化，各種軌道類(lèi)型區(qū)域均發(fā)生變化.現(xiàn)將其變化規(guī)律總結(jié)如下:

1)橢圓軌道區(qū)域總是出現(xiàn)在左下方、右下方.β從0°～180°，左下方區(qū)域減小，右下方區(qū)域增大;B，F(xiàn)初始在左下角，不斷減小至消失，后又從右下角出現(xiàn)，并不斷增大.

2)雙曲軌道區(qū)域總是出現(xiàn)在上部.β從0°～90°，整體區(qū)域增大，P，L先增大后減小;β=75°時(shí) P，L消失，只剩 K;β從 90°～180°，整體區(qū)域減小;P，L又從右上角出現(xiàn).

3)逃逸軌道區(qū)域總是出現(xiàn)在橢圓軌道區(qū)域和雙曲軌道區(qū)域之間.β從0°～67.5°，J，N均減小;β=67.5°時(shí)J，N消失，只剩下I;β從112.5°～180°，J，N再次出現(xiàn)并不斷增大.

4)俘獲軌道區(qū)域沒(méi)有出現(xiàn)在圖中，在ψ為0°～180°中會(huì)出現(xiàn).

5)逆行軌道區(qū)域，β為0°～67.5時(shí)，出現(xiàn)在左上方;β為0°～67.5°時(shí)，出現(xiàn)在右上方.β從0°～67.5°，整體區(qū)域減小;β=67.5°時(shí)消失;β從112.5～180°，F(xiàn)，N，P都再次出現(xiàn)并不斷增大.

圖5 β對(duì)引力輔助軌道類(lèi)型的影響

6)順行軌道區(qū)域一直占據(jù)大部分.β從0°～75°，整體區(qū)域增大，從右側(cè)不斷向左側(cè)擴(kuò)張(主要是K在擴(kuò)張);β為75°～105°時(shí)獨(dú)占整個(gè)平面;β從105°～180°，整體區(qū)域減小，向左側(cè)縮小 (主要是K在縮小).

7)逆行變順行軌道區(qū)域總是出現(xiàn)在逆行軌道區(qū)域和順行軌道區(qū)域之間.β從0°～75°，向上部移動(dòng)，β=75°時(shí)均消失;β為75°～105°時(shí)，不出現(xiàn);β為105°～180°時(shí)，B，J，L 再次出現(xiàn).

8)順行變逆行軌道區(qū)域沒(méi)有出現(xiàn)在圖中，在ψ為0°～180°中會(huì)出現(xiàn).

4.2.2 β對(duì)飛行器能量的影響

圖6所示為β對(duì)飛行器能量變化的影響.參數(shù)設(shè)置如下:μ=0.0121，rp=0.004 95，e=0，ν=0，α =0，ψ =20°，Vp=2.6.

可以看到，曲線呈鐘形，關(guān)于β=180°對(duì)稱(chēng);曲線變化幅度很小，說(shuō)明不同的β引起的軌道能量差別很小.本文對(duì)其他rp值進(jìn)行了模擬，所得曲線仍呈鐘形.

4.2.3 β對(duì)飛行器軌道傾角的影響

圖6 β對(duì)軌道能量變化的影響

圖7和圖8所示為不同rp下β對(duì)飛行器軌道傾角變化的影響.其余參數(shù)設(shè)置與圖6相同.

圖7 rp=0.008時(shí)β對(duì)軌道傾角變化的影響

圖8 rp=0.00495時(shí)β對(duì)軌道傾角變化的影響

從圖7可以看出:曲線關(guān)于β=180°對(duì)稱(chēng).軌道傾角變化的兩個(gè)峰值約在β=150°和β=210°處取得，而當(dāng)β=0°或180°時(shí)，軌道傾角變化為0.從圖8可以看出:此參數(shù)下，軌道傾角變化較大，β=180°時(shí)，軌道傾角變化為180°，說(shuō)明軌道發(fā)生順逆轉(zhuǎn)換.本文對(duì)其他rp值進(jìn)行了模擬，發(fā)現(xiàn)β=180°時(shí)，使軌道發(fā)生順逆轉(zhuǎn)換的臨界rp值約為0.00536.

綜合圖7和圖8，可以看到:①圖形總是關(guān)于β=180°對(duì)稱(chēng);②β=180°時(shí)，可能的Δi只有兩個(gè)值:0°或 180°.

5 以地月系統(tǒng)為例進(jìn)行參數(shù)優(yōu)化

上述研究對(duì)航天器軌道設(shè)計(jì)具有重要指導(dǎo)意義.在特定任務(wù)中，往往需要獲得特定軌道類(lèi)型，如俘獲與逃逸、順行與逆行等，而且往往有一定參數(shù)限制，特別對(duì)Vp，rp會(huì)有一定限制.如果任務(wù)目標(biāo)是對(duì)M2進(jìn)行數(shù)據(jù)采集，則往往需要設(shè)計(jì)使Vp和rp最小;如果需要躲避M2的大氣層、射線等，則需要使Vp和rp最大.

本文以地月系統(tǒng)為例進(jìn)行說(shuō)明.地月系統(tǒng)中，e=0.005 49，μ=0.012 15.限制 rp=0.00504(距月球表200 km)，ν=90°，α =60°，β=30°，限制 ψ 范圍為 180°～360°，分別求逃逸、俘獲的最大Vp.

從圖9可以看出:想要獲得逃逸軌道 (I，J，M，N)，最大 Vp為2.75，軌道類(lèi)型為 I.想要獲得俘獲軌道 (C，D，G，H)，最大Vp為3.17，軌道類(lèi)型為H.

圖9 ψ-Vp圖

6 結(jié)論

為了研究三維引力輔助技術(shù)，本文通過(guò)引入?yún)?shù)α和β，將平面橢圓型限制性三體模型拓展到三維橢圓型限制性三體模型.深入討論了α和β對(duì)飛行器軌道類(lèi)型、軌道能量和軌道傾角的影響.研究結(jié)果表明:

1)α從0°～90°，橢圓軌道、雙曲軌道、逆行軌道、順行軌道整體區(qū)域都增大，逃逸軌道、俘獲軌道、逆行變順行軌道、順行變逆行軌道整體區(qū)域都減小;α從90～180°，變化趨勢(shì)相反.

2)α=90°時(shí)，只存在A，F(xiàn)，K，P四種軌道類(lèi)型，意味著引力輔助前后軌道類(lèi)型保持不變.α=90°或270°時(shí)，引力輔助前后飛行器軌道能量、軌道傾角均保持不變.α=0°或180°時(shí)，引力輔助前后飛行器軌道能量改變最大.

3)限定π ＜ψ ＜2π，則只有當(dāng)0°＜α＜90°時(shí)，才存在逃逸軌道和逆行變順行軌道;只有當(dāng)90°＜α＜180°時(shí)，才存在俘獲軌道和順行變逆行軌道.0＜ψ＜π時(shí)的規(guī)律，可根據(jù)鏡像規(guī)則得到.

4)β從0°～75°，逆行軌道、逆行變順行軌道、順行變逆行軌道整體區(qū)域都減小，順行軌道整體區(qū)域增大;β從150°～180°，變化趨勢(shì)相反.

5)β為75°～105°時(shí)，只存在順行軌道，說(shuō)明引力輔助前后的軌道只可能是順行軌道.

6)不論β取何值，只有當(dāng)π＜ψ＜2π時(shí)，才存在逃逸軌道和逆行變順行軌道;只有當(dāng)0＜ψ＜π時(shí)，才存在俘獲軌道和順行變逆行軌道.

以上的研究和分析揭示了參數(shù)α和β對(duì)飛行器軌道特性的影響.通過(guò)選取適當(dāng)?shù)膮?shù)α和β，可以獲得特定任務(wù)所需要的軌道類(lèi)型、軌道能量和軌道傾角.這些研究可以為引力輔助軌道的設(shè)計(jì)提供依據(jù).

References)

［1］Curtis H D.軌道力學(xué) ［M］.周建華，等譯.北京:科學(xué)出版社，2009:308－309

Curtis H D.Orbital mechanics for engineering students［M］.Translated by Zhou Jianhua，et al.Beijing:Science Press，2009:308－309(in Chinese)

［2］Flandro G A.From instrumented comets to grand tours:on the history of gravity assist［R］.AIAA 2011-0176，2001

［3］Strange N J，Longuski J M.A graphical method for gravity-assist trajectory design ［J］.Journal of Spacecraft and Rockets，2002，39(1):9－16

［4］Chen K J，Kloster K W，Longuski J M.A graphical method for preliminary design of low-thrust gravity-assist Trajectories［R］.AIAA 2008-6952，2008

［5］Carnelli I，Dachwald B，Vasile M.Evolutionary neurocontrol:a novel method for low-thrust gravity-assist trajectory optimization［J］.Journal of Guidance，Control and Dynamics，2009，32(2):615－624

［6］Pascale P，De Vasile M.Preliminary design of low-thrust multiple gravity-assist trajectories［J］.Journal of Spacecraft and Rockets，2006，43(5):1065－1076

［7］Sims J A，Longuski J M，Patel M R.Aerogravity-assist trajectories to the outer planets and the effect of drag［J］.Journal of Spacecraft and Rockets，2000，37(1):49 －55

［8］Armellin R，le Lavagna M，Ercoli-Finzi A.Aero-gravity assist maneuvers:coupled trajectory and vehicle shape optimization［J］.Celestial Mechanics and Dynamical Astronomy，2006，95:391－405

［9］Broucke R A.The celestial mechanics of gravity assist［C］//AIAA/AAS Astrodynamics Conference.Washington D C:American Institute of Aeronautics and Astronautics，1988:69－78

［10］Prado A F B A.Close-approach trajectories in the elliptic restricted problem ［J］.Journal of Guidance，Control and Dynamics，1997，20(4):797 － 802

［11］喬棟，崔平遠(yuǎn)，崔祜濤.基于圓型限制性三體模型的借力飛行機(jī)理研究［J］.宇航學(xué)報(bào)，2009，30(1):82－87

Qiao Dong，Cui Pingyuan，Cui Hutao.Research on gravity-assist mechanism in circular restricted three-body model［J］.Journal of Astronautics，2009，30(1):82－87(in Chinese)

［12］喬棟，崔平遠(yuǎn)，尚海濱.基于橢圓型限制性三體模型的借力飛行機(jī)理研究［J］.宇航學(xué)報(bào)，2010，31(1):36－43

Qiao Dong，Cui Pingyuan，Shang Haibin.Research on gravity-assist mechanism in elliptic restricted three-body model［J］.Journal of Astronautics，2010，31(1):36－43(in Chinese)

［13］Felipe G，Prado A F B A.Classification of out-of-plane swing-by trajectories［J］.Journal of Guidance，Control and Dynamics，1999，22(5):643 －649