彭 坤，孫國江，楊 雷，徐 明，奉振球

(1.中國空間技術(shù)研究院載人航天總體部,北京100094；2.北京航空航天大學宇航學院,北京100191)

1 引言

相比于冒險進入遙遠的太陽系深處,增加人類在地月空間的探索活動更加實際可行。同時,月球及其鄰近空間可作為深空探測技術(shù)的天然試驗場和載人深空探測的中轉(zhuǎn)站［1］,故月球探測活動的實施可為人類深空探測打下堅實的技術(shù)基礎(chǔ)。地月轉(zhuǎn)移軌道設(shè)計是實施月球探測任務(wù)的關(guān)鍵。不同于中心引力體環(huán)繞軌道,地月轉(zhuǎn)移軌道涉及兩個中心引力體,非線性強,設(shè)計難度較大。

目前,已有大量學者對地月轉(zhuǎn)移軌道設(shè)計方法進行研究。 張漢清等［2］和何巍等［3］采用引入太陽引力攝動的弱穩(wěn)定邊界理論,在四體問題模型下設(shè)計了地月低能轉(zhuǎn)移軌道。徐明等［4］則在四體模型基礎(chǔ)上設(shè)計了小推力變軌方式的地月轉(zhuǎn)移軌道。不管是地月低能轉(zhuǎn)移軌道還是小推力地月轉(zhuǎn)移軌道,其轉(zhuǎn)移時間較長,通常為幾十到幾百天?？紤]到轉(zhuǎn)移時間約束,月球探測任務(wù)往往采用兩脈沖變軌形式的地月轉(zhuǎn)移軌道,在近地軌道進行加速,到達近月點后實施制動以進入目標環(huán)月軌道,其轉(zhuǎn)移時間僅為3～5天［5-6］。對于兩脈沖地月轉(zhuǎn)移軌道,一般設(shè)計步驟為先利用低精度模型設(shè)計軌道初值,再在高精度模型中迭代搜索精確軌道。常用的低精度模型有雙二體模型、圓型限制性三體模型和偽狀態(tài)模型。雙二體模型是將地月轉(zhuǎn)移軌道分為地心段和月心段,在月球影響球處進行拼接,并以拼接點處軌道參數(shù)作為設(shè)計變量［7-8］。該方法中設(shè)計變量物理意義不明顯,迭代搜索復雜,且地月轉(zhuǎn)移過程動力學模型不連續(xù),精度稍低。目前圓型限制性三體模型主要用于求解自由返回型的地月轉(zhuǎn)移軌道,根據(jù)自由返回軌道對稱性以近月點參數(shù)作為設(shè)計變量［9-10］。張科等［11］利用圓型限制性三體模型求解了雙脈沖地月轉(zhuǎn)移軌道,但其求得的地月轉(zhuǎn)移軌道轉(zhuǎn)移時間較大,不適合工程應(yīng)用。偽狀態(tài)模型本質(zhì)上是限制性三體問題的封閉近似解法［12］,可較好近似地球月球雙引力場作用下的軌道推演,可用于一般地月轉(zhuǎn)移軌道和自由返回軌道的初值求解［13-14］。 高玉東等［15］則采用分層搜索思想,利用B平面參數(shù)直接在高精度模型中搜索地月轉(zhuǎn)移軌道,但僅針對極月目標軌道進行求解。楊維廉等［16］以地月轉(zhuǎn)移軌道出發(fā)速度、升交點赤經(jīng)和近地點幅角為設(shè)計變量,詳細推導了地月轉(zhuǎn)移軌道的轉(zhuǎn)移矩陣和修正方法,分析了地月轉(zhuǎn)移軌道的特性,但未給出具體的初值估計方法。賀波勇等［17］給出了地月轉(zhuǎn)移軌道半長軸、升交點赤經(jīng)和近地點幅角的初值估計方法,且算例中估計值非常接近精確值,其初值依賴近月點月球位置,沒有反映出地月轉(zhuǎn)移軌道特性。

由文獻［18］知,地月轉(zhuǎn)移過程中航天器所受的地球引力和月球引力影響最大,因此忽略地球和月球的非球形影響,以及太陽引力攝動影響,在三體模型下設(shè)計地月轉(zhuǎn)移軌道不影響其本質(zhì)特性。同時,根據(jù)文獻［16］的研究,與月球軌道面共面的二維地月轉(zhuǎn)移軌道和三維地月轉(zhuǎn)移軌道特性類似,故可通過共面轉(zhuǎn)移情況分析地月轉(zhuǎn)移軌道基本特性,為三維地月轉(zhuǎn)移軌道設(shè)計提供依據(jù)。本文建立三體模型下地心旋轉(zhuǎn)系動力學模型,并以地月轉(zhuǎn)移加速速度增量和地心旋轉(zhuǎn)系對準角為控制變量,提出其二體模型下初值估計解析計算公式,針對不同運行方向地月轉(zhuǎn)移軌道提出相應(yīng)的搜索方法,采用序列二次規(guī)劃算法對地月轉(zhuǎn)移軌道進行求解,分析地月轉(zhuǎn)移軌道的全局特性,以及地月轉(zhuǎn)移軌道與設(shè)計變量之間的關(guān)系,其結(jié)論可應(yīng)用于高保真動力學模型中的精確初值猜測,提高高保真模型下地月轉(zhuǎn)移軌道的設(shè)計效率。

2 系統(tǒng)模型

在三體模型下,假設(shè)地球、月球和航天器均被看作質(zhì)點,月球繞地球做勻速圓周運動,航天器只受到地球和月球的萬有引力影響。該三體模型假設(shè)與圓型限制性三體模型假設(shè)相同。不同于圓型限制性三體模型中將坐標原點選在地月系統(tǒng)中心,為揭示雙脈沖地月轉(zhuǎn)移軌道控制變量與轉(zhuǎn)移時間的關(guān)系,本文以地心為中心建立動力學方程。航天器地月轉(zhuǎn)移加速時刻的地月連線為x軸,由地心指向月心為正,z軸為月球公轉(zhuǎn)運動角動量方向,y軸由右手定則確定,建立地心慣性系,則該坐標系下航天器動力學方程如式(1)所示：

式中,R為航天器相對地球的位置矢量,Rm為航天器相對月球的位置矢量,RME為月球相對地心的位置矢量,μE為地心引力,μM為月心引力。與地心慣性系類似,以地心為中心,地月連線為x軸,可建立地心旋轉(zhuǎn)系。將式(1)從地心慣性系轉(zhuǎn)移到地心旋轉(zhuǎn)系,只考慮xy面內(nèi)運動,并以極坐標系代替直角坐標系可得式(2)：

式中,［r,θ,vr,vθ］ 為地心旋轉(zhuǎn)極坐標系統(tǒng)(圖1)下航天器的極徑、極角、徑向速度和橫向速度,ω為月球繞地球公轉(zhuǎn)角速度,d為地月距離,Rm為航天器與月球的距離。

地月轉(zhuǎn)移軌道初始狀態(tài)為近地停泊圓軌道,滿足軌道高度約束,施加切向脈沖Δv后進入地月轉(zhuǎn)移軌道；終端狀態(tài)為近月點,其月心航跡角γm為0,且需滿足環(huán)月軌道高度約束。設(shè)地月轉(zhuǎn)移加速時刻為0,終端時刻為tf,則其初始條件和目標約束可分別表示為式(3)和(4)：

式中,RE和RM分別為地球和月球的半徑,he和hm分別為近地停泊軌道高度和環(huán)月軌道高度,rm( tf)和vm( tf)分別為地心旋轉(zhuǎn)系下終端時刻航天器相對月球的位置和速度矢量,其表達式如式(5)所示：

3 地月轉(zhuǎn)移軌道求解方法

3.1 求解模型

由系統(tǒng)模型可知,在初始條件中 θ0()和vθ0()未知。本文設(shè)置初始時刻地月轉(zhuǎn)移加速切向脈沖Δv,以及極徑r與地心轉(zhuǎn)移系－x軸的夾角α(簡稱地心旋轉(zhuǎn)系對準角,逆時針方向為正)為控制變量,其與初始條件關(guān)系如式(6)所示：

為提高搜索速度,將航天器相對月心的航跡角作為終止條件,即式(4)的第二分式。同時增加對航跡角變化率的判斷,進一步確定終端點為近月點,如式(7)所示：

轉(zhuǎn)移時間作為目標約束,則相應(yīng)的目標約束變?yōu)槭?8)：

式中,T為目標轉(zhuǎn)移時間。通過以上轉(zhuǎn)換,可將地月轉(zhuǎn)移軌道轉(zhuǎn)化為非線性規(guī)劃問題,控制變量為Δv,α[],目標約束如式(8)所示,狀態(tài)方程如式(2)所示,積分終止條件為式(7)。

3.2 初值估計

地月轉(zhuǎn)移軌道大部分受地球引力影響,在初值估計時可將其假設(shè)為二體模型下的霍曼轉(zhuǎn)移軌道,即一個近地距為RE+he,遠地距為d+RM+hm的地心橢圓軌道。根據(jù)二體模型近地點速度計算公式可算得Δv的值。同時二體模型假設(shè)下,地月轉(zhuǎn)移軌道拱線應(yīng)對準近月點時刻的月球位置,則α應(yīng)為航天器地月轉(zhuǎn)移過程中月球在白道面內(nèi)的掃角。由以上分析可給出控制變量的初值估計解析公式,如式(9)所示?？紤]月球引力的影響,初值估計中適當增加地月轉(zhuǎn)移加速速度增量值。

3.3 搜索策略

通過式(9)可計算出控制變量Δv,α[]的初值,代入求解模型中計算出目標約束式(8)的誤差,通過修正算法不斷修正控制變量,從而得到滿足目標約束的地月轉(zhuǎn)移軌道。本文采用成熟的序列二次規(guī)劃算法SNOPT［19］對此非線性規(guī)劃問題進行求解。

對于二維地月轉(zhuǎn)移軌道,滿足相同近地軌道高度、環(huán)月軌道高度和轉(zhuǎn)移時間約束的軌道共有四條,分別對應(yīng)地心順行月心順行、地心順行月心逆行、地心逆行月心順行和地心逆行月心逆行,如圖2所示(與地心旋轉(zhuǎn)系的旋轉(zhuǎn)方向相同為順行,反之為逆行)。圖2(a)對應(yīng)地心均為順行方向,月心分別為順行和逆行方向的兩條地月轉(zhuǎn)移軌道；圖2(b)對應(yīng)月心均為順行方向,地心分別為順行和逆行方向的兩條地月轉(zhuǎn)移軌道。

對于地心順行和逆行的情況,可通過設(shè)置初始時刻橫向速度vθ方向的方法進行區(qū)分,如式(10)所示：

對于月心順行和逆行的情況,兩者在大部分階段軌跡近似相同,只是在靠近月球段略有不同,若直接采用SNOPT對地月轉(zhuǎn)移軌道進行求解時無法進行區(qū)分。為明確地月轉(zhuǎn)移軌道月心運行方向,本文設(shè)置月心順行和逆行搜索策略來控制其搜索方向。比較月心順行和逆行的地月轉(zhuǎn)移軌道可知,月心逆行軌道的近月點在月球外側(cè),因此其地月轉(zhuǎn)移加速Δv和地心旋轉(zhuǎn)系對準角α的值相對較大。根據(jù)這一特性設(shè)計如下月心順逆行搜索策略：先按標稱程序進行搜索,若滿足月心順逆行要求,則直接輸出結(jié)果；若不滿足月心順逆行要求,則對當前已搜到的控制變量進行修正,再采用SNOPT進行搜索,其修正公式如式(11)所示：

4 仿真驗證

4.1 仿真參數(shù)

設(shè)地球和月球引力常數(shù)分別為 μE＝398 600.4415km3/s和μM＝ 4902．8000 km3/s,地球和月球的半徑分別為 RE＝6378．140 km和RM＝1738．000 km,地月間距離 d ＝384 400 km,月球公轉(zhuǎn)角速度 ω ＝ 2．649 07 ×10－6rad／s。 為驗證地月轉(zhuǎn)移軌道求解方法,設(shè)近地軌道高度he＝200 km,環(huán)月軌道高度hm＝100 km［6］。

4.2 單個地月轉(zhuǎn)移軌道算例

將地月轉(zhuǎn)移軌道設(shè)置為地心順行月心逆行軌道,令轉(zhuǎn)移時間為3天。采用本文提出的搜索算法進行求解,仿真結(jié)果如表1所示。由表1可知,搜索算法僅迭代3步就搜索到地月轉(zhuǎn)移軌道,耗時僅為5.68 s,驗證了算法的可行性和快速性。同時從表1可看出,地月轉(zhuǎn)移速度增量Δv的初值估計和真實值幾乎相同；對準角α的初值估計和真實值僅相差約9°,在0°～360°范圍內(nèi)看相差也較小,說明了初值估計方法的正確性。

表1 單個地月轉(zhuǎn)移軌道搜索結(jié)果Table 1 Search results of single trans-lunar trajectory

通過坐標變換可繪制地月轉(zhuǎn)移軌道分別在地心旋轉(zhuǎn)系和地心慣性系的軌跡,如圖3所示。由圖3(a)可知,地月轉(zhuǎn)移軌道大部分軌跡都在月球影響球之外；由圖3(b)可知,地月轉(zhuǎn)移軌道拱線近似對準近月點時刻的月球位置,驗證了初值估計方法中二體模型假設(shè)的合理性,同時解釋了為何初值估計與真實值誤差較小。

4.3 四種類型軌道對比算例

令轉(zhuǎn)移時間為3天,分別搜索地心順逆行和月心順逆行4種類型組合的地月轉(zhuǎn)移軌道,其搜索結(jié)果如表2所示。由表2可知,對于地心段運行方向,地心順行軌道的地月轉(zhuǎn)移加速Δv比相應(yīng)的地心逆行軌道略小,與初值估計近似相同；地心順行軌道的對準角α比相應(yīng)的地心逆行軌道大,且大于初值估計,而地心逆行軌道對準角α小于初值估計。該結(jié)果說明了地心順行軌道要盡量向近月點時刻月球位置前方瞄準；而地心逆行軌道要盡量向近月點時刻月球位置后方瞄準。對于月心段運行方向,月心順行軌道地月轉(zhuǎn)移加速Δv比相應(yīng)月心逆行軌道小,其對準角α也比相應(yīng)月心逆行軌道略小。該結(jié)果說明了與月心順心軌道相比,月心逆行軌道的對準角要略偏向前方,速度增量略大,驗證了不同類型軌道搜索策略的正確性。

表2 4種類型地月轉(zhuǎn)移軌道搜索結(jié)果Table 2 Search results of four types of trans-lunar trajectories

為了驗證本文搜索算法的魯棒性,對轉(zhuǎn)移時間為1～7天的4種類型地月轉(zhuǎn)移軌道進行搜索,每隔0.5天搜索一次,其搜索結(jié)果如圖4所示。在搜索過程中,除個別軌道搜索時間略長,大部分軌道均在10 s以內(nèi)收斂,說明該算法能適應(yīng)大范圍轉(zhuǎn)移時間變化；同時搜索方向準確,能找到指定類型的地月轉(zhuǎn)移軌道。

5 軌道特性分析

在4種類型地月轉(zhuǎn)移軌道中,地心逆行軌道不能充分利用地球自轉(zhuǎn)速度從而降低火箭運載能力,一般不予采用。對于載人月球探測任務(wù),其地月轉(zhuǎn)移軌道要求設(shè)計為自由返回軌道,而自由返回軌道的近月點為月心逆行軌道。以下主要對地心順行月心逆行的地月轉(zhuǎn)移軌道特性進行分析,為載人月球探測地月轉(zhuǎn)移軌道選取和設(shè)計提供依據(jù)。

5.1 轉(zhuǎn)移時間影響分析

轉(zhuǎn)移時間是地月轉(zhuǎn)移軌道的關(guān)鍵設(shè)計因素［20］,其決定了地月轉(zhuǎn)移加速速度增量和近月制動速度增量的大小。取T＝［1,7］d,近地軌道高度he＝200 km,環(huán)月軌道高度hm＝100 km,分析此轉(zhuǎn)移時間范圍內(nèi)地月轉(zhuǎn)移軌道兩端點速度增量大小,如圖5所示。隨著轉(zhuǎn)移時間的增大,地月轉(zhuǎn)移加速和近月制動速度增量逐漸減小,到3天后兩速度增量變化趨于平緩,最小速度增量出現(xiàn)在5天左右。故對于無人地月轉(zhuǎn)移軌道,要求變軌總速度增量最小,其飛行時間可取為5天轉(zhuǎn)移時間；對于載人地月轉(zhuǎn)移軌道,要綜合考慮速度增量和轉(zhuǎn)移時間,可取為3天轉(zhuǎn)移時間。

在初值估計過程中,地月轉(zhuǎn)移加速速度增量和對準角均基于二體模型進行估值,未考慮月球引力影響。以下比較初值估計與真實值之間的誤差,分析不同轉(zhuǎn)移時間下月球引力對二體模型初值估計的修正量,如圖6所示。對于地月轉(zhuǎn)移加速速度增量,轉(zhuǎn)移時間大于3天時,初值估計與真實值相差不大,且均可近似為常值；轉(zhuǎn)移時間小于3天時,轉(zhuǎn)移時間越小,地月轉(zhuǎn)移加速速度增量真實值與初值估計相差越大。故當轉(zhuǎn)移時間小于3天時,月球引力對地月轉(zhuǎn)移加速速度增量影響大,初值估計時需要進行大幅修正。對于對準角,轉(zhuǎn)移時間小于2天時,對準角隨轉(zhuǎn)移時間遞減,真實值與初值估計相差較大；轉(zhuǎn)移時間大于2天時,對準角隨轉(zhuǎn)移時間遞增,真實值與初值估計誤差逐漸減小,到轉(zhuǎn)移時間為3天時,誤差在10°之內(nèi)。因此,對于轉(zhuǎn)移時間3～7天的地月轉(zhuǎn)移軌道,地月轉(zhuǎn)移加速速度增量近似維持不變,對準角與轉(zhuǎn)移時間近似成線性關(guān)系本文的初值估計方法能對控制變量進行高精度估計。

5.2 近地點高度和近月點高度影響分析

近地點高度和近月點高度決定了月球探測任務(wù)近地停泊軌道和環(huán)月軌道的高度,其選取除考慮火箭運載約束和環(huán)月軌道穩(wěn)定性外,還應(yīng)考慮變軌速度增量大小。以T＝3 d為例,分析不同近地點高度和近月點高度下地月轉(zhuǎn)移軌道地月轉(zhuǎn)移加速和近月制動速度增量變化情況。為提高收斂速度,根據(jù)轉(zhuǎn)移時間影響分析結(jié)果,對對準角α初值估計進行修正,其修正公式為α＝ωT+8°()。

設(shè)近地點高度he＝［200,1000］km,近月點高度hm＝100 km,地月轉(zhuǎn)移軌道變軌速度增量變化曲線如圖7所示。隨著近地點高度增加,地月轉(zhuǎn)移加速速度增量逐漸減小,變化范圍約為200 m/s；近月制動速度增量也逐漸減小,但變化范圍僅為3 m/s,可近似認為不變。

設(shè)近月軌道高度hm＝［100,5000］km,近地軌道高度he＝200 km,地月轉(zhuǎn)移軌道變軌速度增量的變化曲線如圖8所示。隨著近月點軌道高度增加,地月轉(zhuǎn)移加速速度增量逐漸增大,但變化范圍僅為4 m/s,可近似認為不變；而近月制動速度增量逐漸減小,變化范圍約為150 m/s。

6 結(jié)論

本文在三體模型下建立了地心旋轉(zhuǎn)系地月轉(zhuǎn)移軌道求解模型,給出了控制變量初值估計解析式,提出針對地心順逆行和月心順逆行4種類型軌道的搜索策略,采用序列二次規(guī)劃算法進行了求解。仿真結(jié)果表明該方法搜索速度快,魯棒性好,能適應(yīng)1～7天轉(zhuǎn)移時間內(nèi)不同類型地月轉(zhuǎn)移軌道的搜索。當轉(zhuǎn)移時間范圍為3～7天時,初值估計方法能夠?qū)φ鎸嵵颠M行準確估計。同時,根據(jù)初值估計與真實值的誤差,能夠得到月球引力對二體模型解析解的修正量,為后續(xù)高保真模型中地月轉(zhuǎn)移軌道設(shè)計提供更精確的初值。

轉(zhuǎn)移時間增加會導致地月轉(zhuǎn)移軌道地月轉(zhuǎn)移加速和近月制動速度增量的減少,當轉(zhuǎn)移時間為3天時趨于平緩,可將3天～5天作為月球探測任務(wù)中地月轉(zhuǎn)移軌道轉(zhuǎn)移時間設(shè)計約束。近地點高度增加會導致地月轉(zhuǎn)移軌道地月轉(zhuǎn)移加速和近月制動速度增量的減少,主要影響地月轉(zhuǎn)移加速速度增量,近月制動速度增量幾乎不變。近月點高度增加會導致地月轉(zhuǎn)移軌道地月轉(zhuǎn)移加速速度增量的增加和近月制動速度增量的減少,主要影響近月制動速度增量,地月轉(zhuǎn)移加速速度增量幾乎不變。