基于粒子群優(yōu)化的月球捕獲制動(dòng)方法研究

郭繼峰，白成超，賀國(guó)平

（哈爾濱工業(yè)大學(xué)航天工程系，哈爾濱150001）

1 引言

月球探測(cè)是航天技術(shù)發(fā)展的主要領(lǐng)域之一。目前，我國(guó)已經(jīng)成功實(shí)施了月球探測(cè)一期工程（繞）和二期工程（落），三期工程（回）也已立項(xiàng)實(shí)施，而軌道設(shè)計(jì)是探月任務(wù)的關(guān)鍵技術(shù)。月球軌道分為停泊、地月轉(zhuǎn)移、捕獲制動(dòng)、環(huán)繞和著陸五類，其中地月轉(zhuǎn)移軌道相對(duì)于月球?yàn)殡p曲線軌道，在探月衛(wèi)星飛至近月點(diǎn)時(shí)必須進(jìn)行減速機(jī)動(dòng)以被月球引力捕獲，否則衛(wèi)星將飛離地月系統(tǒng)且不會(huì)返回［1］。因此，捕獲軌道是軌道設(shè)計(jì)中的核心之一，直接關(guān)系到軌道器任務(wù)實(shí)施的成敗。我國(guó)嫦娥一號(hào)和嫦娥二號(hào)探月衛(wèi)星采用近月點(diǎn)三次減速機(jī)動(dòng)的方案實(shí)現(xiàn)月球捕獲，最終進(jìn)入平均高度分別為200 km和100 km的極月圓形目標(biāo)軌道［2?3］。 嫦娥三號(hào)探測(cè)器變推力發(fā)動(dòng)機(jī)經(jīng)過(guò)361 s制動(dòng)后，順利進(jìn)入距月面平均高度約100 km的環(huán)月軌道［3］。

探月衛(wèi)星軌道的捕獲制動(dòng)策略可分為固定推力方向制動(dòng)和變推力方向制動(dòng)，它對(duì)捕獲時(shí)燃料消耗、捕獲后目標(biāo)軌道確定等任務(wù)規(guī)劃和探月總體設(shè)計(jì)有直接的影響［4］。另外，由于探測(cè)器可攜帶燃料有限，對(duì)軌道捕獲策略進(jìn)行優(yōu)化設(shè)計(jì)，使完成捕獲任務(wù)的燃料消耗最小也十分重要，常用的優(yōu)化算法主要包括粒子群算法、蟻群算法、模擬退火算法、遺傳算法、差分進(jìn)化算法等。并且，Myatt證明了相比于遺傳算法和粒子群算法，差分進(jìn)化算法在深空引力輔助機(jī)動(dòng)軌道設(shè)計(jì)中的優(yōu)勢(shì)［5］；而B(niǎo)essette比較了粒子群算法和差分進(jìn)化算法在軌道設(shè)計(jì)問(wèn)題中的適用性［6］；Vasile進(jìn)一步提出了一種結(jié)合差分進(jìn)化算法和單純形算法中啟發(fā)式的混合算法，并仿真驗(yàn)證了該混合算法在深空軌道設(shè)計(jì)中的優(yōu)越性［7］。

目前，國(guó)內(nèi)很多學(xué)者也致力于捕獲制動(dòng)段探月軌道設(shè)計(jì)與優(yōu)化策略，并且取得了一定得成果。蔣小勇提出用動(dòng)力學(xué)分解方法來(lái)快速獲得小推力軌道的近似解析解［8］。黃文德以載人登月軌道和中止軌道一體化設(shè)計(jì)為指導(dǎo)，提出了以快速設(shè)計(jì)結(jié)果為初值的載人登月攝動(dòng)軌道設(shè)計(jì)方法，同時(shí)針對(duì)設(shè)計(jì)參數(shù)多、對(duì)初值比較敏感等特點(diǎn)，采取參數(shù)群體優(yōu)化的方法進(jìn)行軌道設(shè)計(jì)［9］。李軍鋒研究了捕獲制動(dòng)結(jié)束后衛(wèi)星繞月球軌道半長(zhǎng)軸和偏心率給定情況下的燃料最優(yōu)問(wèn)題，通過(guò)求解最優(yōu)控制對(duì)應(yīng)的兩點(diǎn)邊值問(wèn)題，得到了燃料最優(yōu)情況下的初始條件、捕獲策略和燃料消耗，驗(yàn)證了最優(yōu)捕獲策略為捕獲發(fā)動(dòng)機(jī)開(kāi)機(jī)后以最大幅值工作至捕獲結(jié)束，推力的方向近似沿著速度相反方向［4］；他還與其合作者基于有限推力模型，以繞月軌道的半長(zhǎng)軸和偏心率為目標(biāo)約束，分析了燃料最優(yōu)捕獲、姿態(tài)勻速轉(zhuǎn)動(dòng)和姿態(tài)慣性定向捕獲策略，并分別利用最優(yōu)控制和粒子群優(yōu)化算法針對(duì)不同捕獲策略進(jìn)行了求解［10］。王云財(cái)?shù)葘?duì)雙曲線軌道的近心點(diǎn)高度進(jìn)行了優(yōu)化處理，并在此基礎(chǔ)上對(duì)整個(gè)捕獲過(guò)程中的推力方向和發(fā)動(dòng)機(jī)的開(kāi)關(guān)機(jī)時(shí)間進(jìn)行優(yōu)化，最終得到整個(gè)捕獲過(guò)程中的控制參數(shù)，進(jìn)一步降低了捕獲過(guò)程中的燃料消耗［11］。

對(duì)于月球捕獲制動(dòng)段軌道設(shè)計(jì)優(yōu)化問(wèn)題，以上提到的方法存在求解分析過(guò)程復(fù)雜的不足，本文在考慮地球攝動(dòng)的前提下，基于固定推力方向制動(dòng)策略進(jìn)行月球捕獲制動(dòng)研究，以燃料消耗以及最終飛行器近月點(diǎn)高度偏差為指標(biāo)，采用簡(jiǎn)單有效的粒子群優(yōu)化算法，進(jìn)行捕獲制動(dòng)軌跡優(yōu)化。

2 月球捕獲制動(dòng)問(wèn)題描述

2.1 動(dòng)力學(xué)模型

在月球J2000慣性坐標(biāo)系下，考慮月球引力和地球的攝動(dòng)力，探測(cè)器軌道動(dòng)力學(xué)方程如式（1）［3］：

其中，m為探測(cè)器的質(zhì)量；Isp為推力器比沖；g為地球重力加速度；r和v分別為探測(cè)器在月球J2000慣性坐標(biāo)系下的位置矢量和速度矢量，且r＝‖r‖；rE為地球在月球J2000慣性坐標(biāo)系下的位置矢量，rE＝‖rE‖，考慮地球?qū)μ綔y(cè)器的攝動(dòng)時(shí)，不考慮地球扁率的影響；rpE為探測(cè)器相對(duì)地球的位置矢量，且rpE＝ r－ rE，μL、μE分別為月球和地球的引力常數(shù)；F為捕獲制動(dòng)時(shí)推力器產(chǎn)生的推力。忽略太陽(yáng)光壓微小攝動(dòng)力。

2.2 固定推力方向制動(dòng)策略

固定推力方向制動(dòng)策略工程上實(shí)現(xiàn)簡(jiǎn)單，是一種很可取的捕獲制動(dòng)方式。定義探測(cè)器近月球點(diǎn)的位置和速度單位矢量分別為＾r(nóng)p和＾vp，常推力方向制動(dòng)中，令推力矢量保持在＾r(nóng)p和＾vp確定的平面內(nèi)，推力方向角定義為α，從＾vp反方向向上測(cè)量為正，從＾r(nóng)p正方向向下測(cè)量為負(fù)，如圖1所示，則存在推力方向平衡關(guān)系如式（2）：

開(kāi)始點(diǎn)火的位置矢量與預(yù)測(cè)近月點(diǎn)的位置矢量的夾角為點(diǎn)火角，如圖2所示。

圖1 推力方向及推力方向角示意圖Fig.1 Schematic diagram of thrust direction and thrust angle

圖2 捕獲制動(dòng)角度關(guān)系與點(diǎn)火角定義Fig.2 Angle relations of braking process and igni?tion angle

假設(shè)開(kāi)始進(jìn)行連續(xù)推力捕獲制動(dòng)時(shí)，飛行器的軌道傾角滿足指標(biāo)要求，僅需進(jìn)行平面內(nèi)制動(dòng)。由此可以確定需要計(jì)算的參數(shù)包括：初始點(diǎn)火角θ、推力方向角α以及點(diǎn)火時(shí)長(zhǎng)Tp。

3 粒子群優(yōu)化

3.1 粒子群優(yōu)化算法

在PSO算法中，優(yōu)化問(wèn)題在所有搜索范圍內(nèi)的解向量都對(duì)應(yīng)著多維度的有限空間里的一個(gè)點(diǎn)，我們稱之為“位置”，而有一定數(shù)量的占據(jù)隨機(jī)位置的初始解向量，被稱為“粒子”，每一個(gè)粒子都占據(jù)一個(gè)位置，每一個(gè)的粒子所在的位置都對(duì)應(yīng)一個(gè)適應(yīng)值，稱為“適應(yīng)度”，每一個(gè)粒子還都有一個(gè)決定下一次迭代時(shí)位置移動(dòng)量的值，稱為“速度”，在迭代過(guò)程中，粒子們就根據(jù)上一次迭代時(shí)的速度、自身經(jīng)歷過(guò)的適應(yīng)度最高的位置和所有粒子經(jīng)歷過(guò)的適應(yīng)度最高的位置決定的速度對(duì)位置進(jìn)行更新［12］。

在一個(gè)D維的搜索空間中，由n個(gè)粒子組成的種群X可以表示為式（3）：

其中第i個(gè)粒子為一個(gè)如式（4）所示的D維向量：

該空間里，第i個(gè)粒子的速度可以表示為式（5）：

在迭代搜索過(guò)程中，定義如下兩個(gè)目標(biāo)位置作為迭代的參考位置：

1）單個(gè)粒子自身在迭代過(guò)程經(jīng)歷過(guò)的適應(yīng)度最高的位置pi（局部最優(yōu)），如式（6）所示：

2）全部粒子在迭代過(guò)程中經(jīng)歷過(guò)的適應(yīng)度最高的位置pg（全局最優(yōu)），如式（7）所示，pg在每次迭代中都是唯一的：

所示，速度和位置更新公式就如式（8）～（9）所示［12］：

其中 i ＝ 1，2，…，n，d ＝ 1，2，…，D；k 為當(dāng)前迭代次數(shù)；ω為慣性權(quán)重；c1，c2為加速常數(shù)；ξ，η為 0，1[ ]區(qū)間的隨機(jī)數(shù)為了防止粒子的盲目搜索，一般建議將位置或者速度限制在一定的搜索區(qū)間。PSO算法流程如圖3所示。

圖3 PSO算法流程Fig.3 Flow chart of PSO algorithm

如前所述，對(duì)于捕獲制動(dòng)軌道優(yōu)化問(wèn)題，需要計(jì)算的參數(shù)包括初始點(diǎn)火角θ、推力方向角α以及點(diǎn)火時(shí)長(zhǎng)Tp三個(gè)。由于制動(dòng)過(guò)程中軌道偏心率單調(diào)下降，所以優(yōu)化過(guò)程中，將軌道偏心率的值作為制動(dòng)過(guò)程結(jié)束的判斷條件，故點(diǎn)火時(shí)長(zhǎng)不作為粒子群優(yōu)化的變量，優(yōu)化變量只有點(diǎn)火角和推力方向角，而約束只剩下近月點(diǎn)高度一個(gè)。因子此處D ＝ 2，Xi＝ [ Xi1Xi2]，Xi1表示推力方向角，Xi2表示點(diǎn)火角；vi＝ [ vi1vi2]，對(duì)應(yīng)為兩個(gè)優(yōu)化變量的變化速度。

3.2 慣性權(quán)重的選擇

慣性權(quán)重ω體現(xiàn)的是粒子繼承先前的速度的能力，一個(gè)較大的慣性權(quán)值有利于全局搜索，而一個(gè)較小的慣性權(quán)值更有利于局部搜索。為了更好的平衡算法的全局搜索與局部搜索能力，采用式（10）所示慣性權(quán)重選擇方法［12］：

其中ωstart為初始慣性權(quán)重，ωend為迭代至最大次數(shù)時(shí)的慣性權(quán)重，k為當(dāng)前迭代次數(shù)，Tmax為最大迭代次數(shù)。一般來(lái)說(shuō)，慣性權(quán)值ωstart＝0.9，ωend＝0.4左右算法性能最好。這樣，隨著迭代的進(jìn)行，慣性權(quán)重由0.9線性遞減至0.4，迭代初期較大的慣性權(quán)重使算法保持了較強(qiáng)的全局搜索能力，而迭代后期較小的慣性權(quán)重有利于算法進(jìn)行更精確的局部搜索。

4 基于粒子群優(yōu)化的固定推力制動(dòng)方法驗(yàn)證

4.1 參數(shù)設(shè)置

在Visual Studio 2015下，編寫C＋＋仿真主程序，進(jìn)行仿真。根據(jù)月球半徑及其影響球半徑，假設(shè)給定以下仿真參數(shù)［13］：

末端近月點(diǎn)距離約束：rpf＝1938 km；

最終繞月橢圓軌道離心率約束：e＝0.6845 ；

優(yōu)化指標(biāo)：J＝min( Tp)。

固定推力制動(dòng)策略中，主要參數(shù)有點(diǎn)火角，推力方向角與點(diǎn)火時(shí)間（制動(dòng)時(shí)間）三個(gè)，因此仿真設(shè)置了3個(gè)優(yōu)化變量，2個(gè)等式約束，1個(gè)優(yōu)化指標(biāo)。由于制動(dòng)過(guò)程中軌道偏心率單調(diào)下降，所以優(yōu)化過(guò)程中，將軌道偏心率的值作為制動(dòng)過(guò)程結(jié)束的判斷條件，故點(diǎn)火時(shí)長(zhǎng)不作為粒子群優(yōu)化的變量，所以優(yōu)化變量只有點(diǎn)火角和推力方向角，而約束只剩下近月點(diǎn)高度一個(gè)。使用PSO尋優(yōu)算法優(yōu)化，通過(guò)罰函數(shù)將等式約束去掉。粒子群優(yōu)化基本參數(shù)設(shè)定如表1所示。

表1 參數(shù)搜索區(qū)間Table 1 Search interval of parameters

從而得到新的粒子群優(yōu)化指標(biāo)如式（11）：

假定初始雙曲線軌道參數(shù)如表2所示，在理想姿態(tài)無(wú)誤差情況下，探測(cè)器飛入月球主引力區(qū)并制動(dòng)，制動(dòng)前探測(cè)器質(zhì)量為120 t，主發(fā)動(dòng)機(jī)推力大小為200 kN，比沖為900 s，優(yōu)化軌道制動(dòng)參數(shù)，考察制動(dòng)結(jié)果和過(guò)程。

表2 初始雙曲線軌道近月點(diǎn)參數(shù)Table 2 Initial hyperbolic orbit near point parameters

計(jì)算地球引力時(shí)，使用JPL DE421星歷表查詢地月間的相對(duì)運(yùn)動(dòng)參數(shù)，并假定起始時(shí)刻為2017年4月21日0時(shí)0分0秒。

4.2 仿真結(jié)果

優(yōu)化得到的制動(dòng)變軌三要素如表3所示：

表3 固定推力方向制動(dòng)策略參數(shù)Table 3 Parameters of braking strategy in fixed thrust direction

制動(dòng)后軌道六根數(shù)如表4所示：

表4 制動(dòng)后軌道六根數(shù)Table 4 Six orbital elements after braking

制動(dòng)后指標(biāo)變量如表5所示。

表5 制動(dòng)結(jié)束后指標(biāo)變量Table 5 Indicator variables after braking

近月點(diǎn)處的距離、速度、離心率、軌道傾角的設(shè)定值、真實(shí)值及設(shè)定值和真實(shí)值之間的偏差如表6所示。

表6 多種參數(shù)的偏差Table 6 Deviation of parameters

由表6可看出，制動(dòng)結(jié)束后，近月點(diǎn)高度誤差為0.000 27 km，小于100 m；近月點(diǎn)速度偏差為 0.000 513 km／s，小于 1 m／s；離心率的偏差為0.001 058，誤差百分率約為1%；軌道傾角偏差為0.000 11°，小于 0.05°。 以上結(jié)果說(shuō)明，固定推力方向制動(dòng)策略可以通過(guò)優(yōu)化，在消耗推進(jìn)劑盡可能少的同時(shí)滿足指定的終端狀態(tài)約束。

從飛行器進(jìn)入月球影響球開(kāi)始到捕獲制動(dòng)結(jié)束后飛行器繞月球飛行近一圈為止飛行器的運(yùn)動(dòng)軌跡仿真圖如圖4所示。其中紅色曲線為制動(dòng)前的軌跡，藍(lán)色曲線為制動(dòng)過(guò)程中的軌跡，綠色曲線為捕獲制動(dòng)結(jié)束后的繞月球飛行的橢圓軌跡。

圖4 從進(jìn)入月球影響球開(kāi)始的整體捕獲制動(dòng)軌跡Fig.4 Capture trajectory after entering the influence ball of the moon

圖5 為捕獲制動(dòng)過(guò)程軌跡的放大圖，同樣，紅色曲線為制動(dòng)前的軌跡，藍(lán)色曲線為制動(dòng)過(guò)程中的軌跡，綠色曲線為捕獲制動(dòng)結(jié)束后的繞月球飛行的橢圓軌跡的一部分。

制導(dǎo)過(guò)程中以離心率減小到一定值為終止制動(dòng)條件，離心率隨時(shí)間的變化如圖6所示。

圖5 捕獲制動(dòng)過(guò)程軌跡放大圖Fig.5 Enlarged view of capture process trajectory

圖6 捕獲制動(dòng)過(guò)程中離心率隨時(shí)間的變化Fig.6 Change of eccentricity over time during cap?ture process

優(yōu)化指標(biāo)是使燃耗最少，因此，捕獲制動(dòng)過(guò)程中，對(duì)飛行器質(zhì)量變化的追蹤也至關(guān)重要。飛行器質(zhì)量隨時(shí)間變化的過(guò)程如圖7所示。

圖7 飛行器質(zhì)量隨時(shí)間變化Fig.7 Mass change of the spacecraft over time

捕獲制動(dòng)過(guò)程中，軌道能量和軌道偏心率是單調(diào)下降的，由于是平面內(nèi)制動(dòng)，軌道傾角只由攝動(dòng)力影響，所以軌道傾角變化很小，由軌道近月點(diǎn)高度在固定推力方向制動(dòng)完成時(shí)正好達(dá)到200 km的目標(biāo)近月點(diǎn)高度。

5 結(jié)論

本文在考慮地球攝動(dòng)的情況下，基于固定推力方向制動(dòng)策略，以燃料消耗以及最終飛行器近月點(diǎn)高度偏差為指標(biāo)，以最終繞月軌道偏心率以及最終繞月軌道近心點(diǎn)高度為約束，采用粒子群算法進(jìn)行捕獲制動(dòng)段的軌道優(yōu)化。仿真結(jié)果表明，本文方法可有效完成月球?qū)μ綔y(cè)器的捕獲，探月衛(wèi)星初始以雙曲線軌道進(jìn)入月球影響球，通過(guò)捕獲制動(dòng)，最終到達(dá)繞月大橢圓軌道上。最終狀態(tài)約束誤差表明，本方法能在滿足所有約束的情況下最小化燃耗，這驗(yàn)證了基于粒子群優(yōu)化的固定推力方向制動(dòng)方法的有效性。此外，本文沒(méi)有復(fù)雜的方程分析，這說(shuō)明了本文方法的簡(jiǎn)單與易實(shí)現(xiàn)性。

目前在考慮地球?qū)μ綔y(cè)器的攝動(dòng)影響時(shí)，沒(méi)有考慮地球的扁率，認(rèn)為地球是一個(gè)規(guī)則的球體，此外沒(méi)有考慮太陽(yáng)及光壓等攝動(dòng)項(xiàng)的影響，在今后的工作中，將進(jìn)一步考慮這些攝動(dòng)影響，建立更為精確的動(dòng)力學(xué)模型。同時(shí)，將對(duì)多種捕獲制動(dòng)策略展開(kāi)研究，并進(jìn)行對(duì)比分析。

［1］ 李革非，韓潮.月球捕獲控制分析［C］／／全國(guó)第十二屆空間及運(yùn)動(dòng)體控制技術(shù)學(xué)術(shù)會(huì)議，2006：240?244.

Li Gefei， Han Chao.Analysis on Lunar capture control［C］ ／／Twelfth National Conference on Space and Movement Control Technology， 2006： 240?244.（in Chinese）

［2］ 王宏，董光亮，胡小工，等.月球探測(cè)器捕獲分析及評(píng)估［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)， 2007，24（6）： 406?409.

Wang Hong， Dong Guangliang， Hu Xiaogong， et al.Analysis and assessment of the Lunar probe captured by Moon［J］.Journal of Zhengzhou Institute of Surveying and Mapping，2007， 24（6）： 406?409.（in Chinese）

［3］ 盛偉強(qiáng).火星近火點(diǎn)捕獲制動(dòng)姿軌一體化控制［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2015.

Sheng Weiqiang.Study on Coupling Control of Attitude and Orbit for Mars Orbit Insertion at Periareon［D］.Harbin： Har?bin Institute of Technology， 2015.（in Chinese）

［4］ 李軍鋒.基于精確動(dòng)力學(xué)模型的火星探測(cè)軌道設(shè)計(jì)研究［D］.北京：清華大學(xué)，2012.

Li Junfeng.Research on Orbit Design for Mars Explorer by the Accurate Dynamic［D］.Beijing： Tsinghua University，2012.（in Chinese）

［5］ Myatt D R，Becerra V M，Nasuto SJ，et al.Advanced global optimization tools for mission analysis and design［R］.Final Report of ESA Ariadna ITT AO4532 ／18138 ／04 ／NL／MV，2004.

［6］ Bessette C R，Spencer D B.Performance comparison of sto?chastic search algorithms on the interplanetary gravity?assist trajectory problem［J］.Journal of Spacecraft and Rockets，2015， 44（3）： 722?724.

［7］ Vasile M，Minisci E.Analysis of some global optimization al?gorithms for space trajectory design［J］.Journal of Spacecraft and Rockets， 2010， 47（2）： 334?344.

［8］ 蔣小勇.深空探測(cè)器小推力軌道優(yōu)化設(shè)計(jì)方法及其應(yīng)用研究［D］.長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué)，2014.

Jiang Xiaoyong.Study on Low?Thrust Trajectory Design and Optimization Method for Deep Space Explorer and Application［D］.Changsha： National University of Defense Technology，2014.（in Chinese）

［9］ 黃文德.載人登月中止軌道的特性分析與優(yōu)化設(shè)計(jì)［D］.長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué)，2011.

Huang Wende.Characteristic Analysis and Optimal Design of Abort Trajectories for Manned Lunar Landing Mission［D］.Changsha： National University of Defense Technology， 2011.（in Chinese）

［10］ 李軍鋒，龔勝平.有限推力模型火星探測(cè)捕獲策略分析［J］.中國(guó)科學(xué)： 物理學(xué) 力學(xué) 天文學(xué)，2013（06）：781?786.

Li Junfeng， Gong Shengping.Analysis of capture strategies for Mars explorer with finite?thrust［J］.Sci Sin?Phys Mech As?tron， 2013（06）： 781?786.（in Chinese）

［11］ 王云財(cái)，劉輝，胡曉賽.一種改進(jìn)的火星探測(cè)捕獲段軌道設(shè)計(jì)方法［C］／／全國(guó)航天動(dòng)力學(xué)與控制學(xué)術(shù)會(huì)議.2015.

Wang Yuncai， Liu Hui， Hu Xiaosai.An improved Martian exploration capture trajectory design method［C］／／National Conference on Space Dynamics and Control， 2015.（in Chi?nese）

［12］ 史峰，王輝，郁磊，等.MATLAB智能算法30個(gè)案例分析［M］.北京：北京航空航天大學(xué)出版社，2011：102?107.

Shi Feng， Wang Hui， Yu Lei， et al.30 Cases Analysis of MATLAB Intelligent Algorithm［M］.Beijing： Beihang Uni?versity Press， 2011： 102?107.（in Chinese）

［13］ 郗曉寧，曾國(guó)強(qiáng)，任萱，等.月球探測(cè)器軌道設(shè)計(jì)［M］.北京：國(guó)防工業(yè)出版社，2001：240?248.

Xi Xiaoning， Zeng Guoqiang， Ren Xuan， et al.Lunar Explo?ration Orbit Design［M］.Beijing： National Defence Industry Press， 2001： 240?248.（in Chinese）