基于普適變量法的火星探測器軌道初步設(shè)計及仿真

周 杰，張樹瑜，劉付成

（1.上海航天控制工程研究所，上海 200233；2.上海航天技術(shù)研究院，上海 201109）

0 引言

火星探測是深空探測中的一個熱點(diǎn)領(lǐng)域，歷來受到世界航天大國的重視。嫦娥探月工程拉開了我國深空探測的序幕，火星探測成為下一個目標(biāo)?；鹦翘綔y器軌道初步設(shè)計方法中，重點(diǎn)是求解蘭伯特問題。GAUSS最早獲得了巨大成功，但其算法遠(yuǎn)不能滿足當(dāng)代航天動力學(xué)的要求［1］。BATTIN，VAUGHAN分別改進(jìn)了Guass算法，使它對所有可能實(shí)現(xiàn)的軌道均收斂，同時滿足較高的數(shù)值精度，但計算方法較復(fù)雜［2－3］。BATE，BOND 提出并改進(jìn)的普適變量法，其優(yōu)點(diǎn)是使用同一組公式就可計算各種圓錐曲線的軌道，且計算過程較簡單，便于計算機(jī)編程，同時其精度和收斂性也較好［4－5］。本文基于該方法對基于普適變量法的火星探測器軌道初步設(shè)計及仿真進(jìn)行了研究。

1 火星探測器軌道初步設(shè)計

根據(jù)能量獲取方式，火星探測器軌道可分為大推力變軌、小推力變軌和利用天體引力輔助變軌三種基本類型。本文討論典型大推力變軌方式。

探測器在從地球逃逸軌道向環(huán)繞火星的目標(biāo)軌道飛行過程中，始終受到地球、太陽、火星和其他行星引力的影響。這是一個N體問題，無法獲得解析解。本文引入影響球概念，影響球模型是將N體問題轉(zhuǎn)為一系列的二體問題，可描述為：每個大質(zhì)量的天體都有一定的引力支配區(qū)域，探測器在一個天體影響球內(nèi)只受該天體單獨(dú)的引力作用，當(dāng)探測器遠(yuǎn)離該天體時，飛行到該天體引力影響球邊界，則開始進(jìn)入另一個天體的引力支配區(qū)域。

根據(jù)影響球概念，完整的火星探測器軌道需經(jīng)過地球發(fā)射和逃逸段的地球影響球、地火轉(zhuǎn)移段的太陽影響球和被火星捕獲段的火星影響球三個不同的區(qū)域。采用的求解影響球模型方法為圓錐曲線拼接法，它基于假設(shè)：

a）在探測器飛行的任意時刻，探測器僅受到一個天體引力作用；

b）假設(shè)地球和火星的運(yùn)行軌道為同面的圓軌道。

在上述假設(shè)下，地火轉(zhuǎn)移能量最省軌道是霍曼轉(zhuǎn)移軌道，其近日點(diǎn)位于地球的日心軌道軌道上，遠(yuǎn)日點(diǎn)位于火星的日心軌道上，如圖1所示。

圖1 火星探測器軌道Fig.1 Trajectory of Mars probe

1.1 地火轉(zhuǎn)移軌道段

地火轉(zhuǎn)移軌道的初步設(shè)計步驟如下。

步驟1：確定從地球出發(fā)日期和到達(dá)火星日期。

步驟2：查星歷表，確定出發(fā)時在日心慣性坐標(biāo)系中地球的位置R1和速度v1，以及到達(dá)時火星的位置R2和速度v2。

步驟3：求解蘭伯特問題，確定出發(fā)和到達(dá)時探測器的雙曲剩余速度。

步驟1中，理論上出發(fā)日期可任選，但實(shí)際受能量的限制，只能在特定的一段時間內(nèi)選取，稱之為發(fā)射窗口。對不同的出發(fā)和到達(dá)日期，探測器實(shí)現(xiàn)地火轉(zhuǎn)移所需的能量均不同。實(shí)際上，能量最省發(fā)射機(jī)會出現(xiàn)的頻率取決于地球和火星的會合周期。

設(shè)地球和火星環(huán)繞太陽運(yùn)動的角速度分別為ωE，ωM，則角速度差

地球和火星相對位置重復(fù)出現(xiàn)兩次的時間間隔

若令PE，PM分別為地球和火星環(huán)繞太陽運(yùn)動的軌道周期，則地球和火星的會合周期還可表示為

由式（2）或式（3）可得，地球和火星的會合周期為777.9d，即能量最省的發(fā)射機(jī)會約每26個月會出現(xiàn)1次。發(fā)射日期可選在能量最省發(fā)射日期附近，此時地球、火星與太陽連線間的夾角約44°。

步驟2中，查星歷表，得到出發(fā)和到達(dá)日期地球和火星在日心黃道慣性坐標(biāo)系（X軸指向J2000春分點(diǎn)方向，Z軸沿黃道面法向，Y軸由右手準(zhǔn)則確定）中的狀態(tài)向量。

步驟3中，確定出發(fā)和到達(dá)時探測器的雙曲剩余速度v∞。該速度是指探測器在地球影響球邊界處相對地球的速度?？蓪⒉襟E2中求出的出發(fā)時地球的日心位置矢量和到達(dá)時火星的日心位置矢量近似視作飛行器在離開地球影響球和到達(dá)火星影響球時的位置矢量，通過求解蘭伯特問題可得探測器在地球影響球邊界處（或火星影響球邊界處）的日心位置矢量和日心速度矢量，這樣就可解得探測器的逃逸雙曲剩余速度（或到達(dá)雙曲剩余速度）。

求解蘭伯特問題是地火轉(zhuǎn)移軌道設(shè)計中的關(guān)鍵。蘭伯特問題可描述為：橢圓弧上兩點(diǎn)間的飛行時間只取決于橢圓的半長軸、弧上兩點(diǎn)至焦點(diǎn)的距離之和，以及連接弧上兩點(diǎn)的弧長。其數(shù)學(xué)表達(dá)式為

式中：t為橢圓弧上兩點(diǎn)間的飛行時間；a為橢圓軌道的半長軸；r1，r2為弧上兩點(diǎn)到焦點(diǎn)的距離；c為連接弧上兩點(diǎn)的弧長。對蘭伯特問題，一般解法為：給定r1，r2以及從r1至r2的飛行時間t，求航天器在r1，r2處的速度矢量v1，v2［6］。

由共面矢量的基本定理：若A，B，C為共面矢量，且A，B不共線，則C可由A，B線性組合表示。因開普勒運(yùn)動限制在一平面內(nèi)，有

式中：f，g為拉格朗日系數(shù)?？赏茖?dǎo)出拉格朗日系數(shù)與真近點(diǎn)角變化值Δθ間的關(guān)系以及v1，v2的表達(dá)式為

式中：μ為中心天體引力常數(shù)；L為角動量，且L＝｜r1×v2｜［6］。

拉格朗日系數(shù)確定后，蘭伯特問題即可獲解。本文用BATE的普適變量法求解蘭伯特問題。

設(shè)普適變量χ與近點(diǎn)角間的關(guān)系為（近地點(diǎn)處t0＝0）

式中：E，F(xiàn)分別為橢圓和雙曲線軌道的偏近點(diǎn)角。

設(shè)t0為χ等于零的時刻，則時刻t0＋Δt的χ值可由普適開普勒方程

迭代解出。此處：r0，vr0分別為時刻t＝t0的半徑和徑向速度；α為長半軸a的倒數(shù)，α＜0，α＝0，α＞0分別對應(yīng)為雙曲線、拋物線和橢圓；C（z），S（z）為Stumpff函數(shù)，且

其中：z＝αχ2，z＜0，z＝0，z＞0分別對應(yīng)為雙曲線、拋物線和橢圓。

用χ，C（z），S（z）表示的拉格朗日系數(shù)為

上述各式中未知量為L，χ，z，已知量為Δθ，Δt，r1，r2。用普適變量法求解蘭伯特問題時需用逐次逼近法，求解過程如下：

a）計算r1，r2，Δθ。

c）令

e）計算拉格朗日系數(shù)

f）計算v1，v2。

至此，采用普適變量法求解蘭伯特問題過程結(jié)束。

在上述第二步中已求得地球和火星在日心慣性系中的狀態(tài)向量Rearth，vearth，RMars，vMars。在第三步中，假設(shè)探測器在地球影響球邊界處的日心位置矢量即為地球的位置矢量R1＝Rearth，在火星影響球邊界處的日心位置矢量即為火星的位置矢量R2＝RMars，求解蘭伯特問題得到探測器在地球影響球邊界處的日心速度vdeparture和在火星影響球邊界處的日心速度varrive。則逃逸雙曲剩余速度v∞，departure（探測器在地球影響球邊界處相對地球的速度）和到達(dá)雙曲剩余速度v∞，arrive（探測器在火星影響球邊界處相對火星的速度）分別為

引入深空探測器軌道設(shè)計中的重要參數(shù)C3，定義C3＝（v∞）2。設(shè)探測器在半徑為r的地球停泊圓軌道上加速進(jìn)入一逃逸雙曲軌道，其速度增量可表示為

可見，發(fā)射C3越大，探測器逃逸地球引力場所需的速度增量就越大，其大小表征了轉(zhuǎn)移軌道對發(fā)射能量需求。同樣，到達(dá)C3與探測器到達(dá)火星近火點(diǎn)時制動所需的Δv直接相關(guān)，表征了轉(zhuǎn)移軌道對探測器機(jī)動能力需求。地火轉(zhuǎn)移軌道設(shè)計中，節(jié)省能量是考慮的首要問題，故希望的發(fā)射C3和到達(dá)C3較小。

不同出發(fā)日期和到達(dá)日期對應(yīng)的發(fā)射C3和到達(dá)C3不同。本文選定一段范圍的出發(fā)日期和到達(dá)日期（如出發(fā)日期為2013年10月1日至2014年3月1日，到達(dá)日期為2014年5月20日至2015年1月20日），用上述解法求得每一對出發(fā)日期和到達(dá)日期相對應(yīng)的發(fā)射C3和到達(dá)C3，結(jié)果分別如圖2、3所示。

圖2 地球發(fā)射C3等高圖Fig.2 Lnuch C3contour map

圖3 到達(dá)火星C3等高圖Fig.3 Arrival C3contour map

由圖2、3可知：發(fā)射能量呈山脊分布，共分為上下兩個區(qū)域，每個區(qū)域有一局部極小點(diǎn)，其中下半?yún)^(qū)對應(yīng)的是Ⅰ類轉(zhuǎn)移軌道，其日心轉(zhuǎn)角小于180°，上半?yún)^(qū)對應(yīng)的是Ⅱ類轉(zhuǎn)移軌道，其日心轉(zhuǎn)角大于180°。相對而言，Ⅱ類轉(zhuǎn)移軌道的最小能量略小于Ⅰ類轉(zhuǎn)移軌道，但所花時間更長，兩者各有優(yōu)缺點(diǎn)。

同理可得到達(dá)火星C3等高圖。對某些探測任務(wù)，常希望用較小的速度脈沖制動就能形成環(huán)繞火星的目標(biāo)軌道，這樣就可通過到達(dá)火星C3等高圖確定合適的發(fā)射機(jī)會。

另外，從發(fā)射和到達(dá)C3等高圖還可知：在2013年12月1日至2013年12月20日期間發(fā)射所需的C3能量較小，且到達(dá)火星制動所需的能量也較小，相對應(yīng)的到達(dá)火星的日期為2014年10月1日至2014年10月30日。對于發(fā)射質(zhì)量約1 000kg的探測器，我國CZ－3A運(yùn)載火箭所能提供的最大C3能量大小約9.5km2／s2，即在發(fā)射C3等高圖中C3值小于9.5的區(qū)域均可行。如以2013年12月15日為出發(fā)日期，2014年10月20日為到達(dá)日期，由蘭伯特問題求解得到發(fā)射C3≈8.955 49km2／s2，到達(dá)C3≈12.052 8km2／s2，飛行時間約309d。

1.2 地球逃逸軌道段

作過地心與v∞平行的直線，繞此直線旋轉(zhuǎn)逃逸雙曲線可得一旋轉(zhuǎn)曲面，其近地點(diǎn)在空間形成一圓形軌跡，該曲面上任一雙曲線均可作為逃逸雙曲線。

為保證發(fā)射面與地球停泊軌道面及逃逸雙曲線在同一平面內(nèi)，設(shè)A為過地心和v∞反向延長線上停泊軌道高度上的點(diǎn)（如圖4所示），則發(fā)射場、地心和點(diǎn)A確定的平面即為發(fā)射軌道面，探測器應(yīng)在點(diǎn)A從發(fā)射軌道進(jìn)入停泊軌道，滑行至前述圓形軌跡上任一點(diǎn)，即在逃逸雙曲線近地點(diǎn)上加速進(jìn)入逃逸軌道。

圖4 發(fā)射方位角和發(fā)射窗口Fig.4 Lunch azimuth and Lunch window

對某一發(fā)射點(diǎn)，因安全考慮對發(fā)射方位角有限制。1d中，待發(fā)射點(diǎn)隨地球自轉(zhuǎn)至點(diǎn)1時，對應(yīng)發(fā)射的最小方位角，此時發(fā)射窗口打開；轉(zhuǎn)至點(diǎn)2時，對應(yīng)最大方位角，此時發(fā)射窗口關(guān)閉。

1.3 中途修正原理

根據(jù)前述圓錐曲線拼接法設(shè)計了初步的地火轉(zhuǎn)移軌道，但由于發(fā)射產(chǎn)生的誤差及探測器在飛行過程中受到的各種攝動力的影響，實(shí)際軌道會產(chǎn)生偏差，偏離預(yù)定的軌道，若不進(jìn)行軌道修正，則其到達(dá)火星時會由于非線性效應(yīng)產(chǎn)生巨大偏差。計算表明：在300km的地球停泊軌道，近地點(diǎn)發(fā)動機(jī)關(guān)機(jī)時速度變化0.01%（1.1m／s）時，日心霍曼轉(zhuǎn)移橢圓目標(biāo)軌道半徑將變化153 000km；發(fā)動機(jī)關(guān)機(jī)時近地點(diǎn)半徑變化0.01%（0.67km）時也將產(chǎn)生偏差70 000km。因此，為到達(dá)預(yù)定的火星探測目標(biāo)軌道，須在地火巡航段完成對探測器進(jìn)行多次軌道修正，一般地火轉(zhuǎn)移段需作2～3次軌道修正，首次修正一般在探測器脫離地球影響球后初期進(jìn)行，以后根據(jù)實(shí)際再進(jìn)行數(shù)次微小修正。

軌道修正一般是改變探測器的速度。如當(dāng)探測器飛行到火星的影響球范圍時，其相對火星的速度可認(rèn)為是v∞，則通常對其進(jìn)行調(diào)整，設(shè)其改變量為Δv∞，以使調(diào)整后的探測器按預(yù)定的要求到達(dá)火星附近（如軌道傾角為90°，近火點(diǎn)高度為250km），但該關(guān)系為非線性，利用微分修正法方法尋找合適的Δv∞較難，有時甚至不能收斂。

研究發(fā)現(xiàn)用迭代算法修正B平面上的誤差時收斂性非常好，易計算獲得修正誤差所需的速度增量，因此被廣泛用于行星際軌道的制導(dǎo)和精確設(shè)計。B平面是通過目標(biāo)星中心且垂直于雙曲軌道入射漸近線矢量S＾的平面（如圖5所示），在B平面上建立一平面坐標(biāo)系，其原點(diǎn)O為目標(biāo)星質(zhì)心；N＾為行星赤道面法向或黃道面發(fā)向，則B平面上兩坐標(biāo)軸單位矢量為

式中：b＝‖B‖；θ為雙曲軌道入射漸近線與火星近焦點(diǎn)坐標(biāo)系Xp間的夾角；iM，ωM分別為探測器軌道相對于火星赤道的軌道傾角和近火點(diǎn)角距。

圖5 B平面Fig.5 B－plane

則由式（25）、（26），B平面參數(shù)可表示為

當(dāng)iM＝90°（目標(biāo)軌道為火星極地軌道）時，可得特殊形式

由上建立了B平面參數(shù)和軌道狀態(tài)參數(shù)間的關(guān)系，在地火轉(zhuǎn)移軌道的終端約束條件（一般是軌道傾角和近火點(diǎn)高度）給定后，可將其轉(zhuǎn)換為B平面參數(shù)，再用不同搜索算法進(jìn)行計算，直至B平面參數(shù)得到滿足，此時的轉(zhuǎn)移軌道便為期望軌道。

2 仿真

基于本文方法求得的轉(zhuǎn)移軌道參數(shù)，用STK軟件的軌道機(jī)動模塊astrogator工具設(shè)計一完整的火星探測器軌道：探測器在任意時刻的受力模型為真實(shí)力模型，即考慮各種星體的引力攝動，并在轉(zhuǎn)移軌道段增加兩次中途修正。采用的星歷為JPL的DE405星表，中途修正方法為B平面法，目標(biāo)軌道為近火點(diǎn)高度250km，軌道周期50h的大橢圓極地軌道。

選擇發(fā)射日期為2013年12月15日，到達(dá)日期為2014年10月20日，用蘭伯特定理可得地火轉(zhuǎn)移軌道 的 初 始 參 數(shù) 為：C3＝8.955 49km2／s2，α＝203.84°，δ＝25.006 2°。假設(shè)探測器從西昌衛(wèi)星發(fā)射中心（東經(jīng)102.0°，北緯28.3°）發(fā)射，先發(fā)射至高度300km的圓停泊軌道，然后滑行至合適的地點(diǎn)點(diǎn)火加速進(jìn)入地球逃逸雙曲軌道。雙曲軌道的參數(shù)滿足地火轉(zhuǎn)移軌道的初始參數(shù)，發(fā)射方位角98°。地火轉(zhuǎn)移階段進(jìn)行2次軌道修正。第一次在探測器剛離開地球影響球時進(jìn)行，第二次在探測器剛進(jìn)入火星影響球時進(jìn)行。在火星軌道捕獲階段，當(dāng)探測器運(yùn)行到近火點(diǎn)為250km時，給其一速度增量進(jìn)行制動，使探測器進(jìn)入一周期為50h的大橢圓極火軌道。仿真所得參數(shù)見表1，地球逃逸軌道、地火轉(zhuǎn)移軌道和火星捕獲軌道如圖6所示。

3 結(jié)束語

本文對火星探測器軌道的設(shè)計進(jìn)行了研究。先求解了基于二體假設(shè)的蘭伯特問題，得到了能量較優(yōu)的發(fā)射窗口，再討論了基于B平面法的中途軌道修正，并用STK軟件模擬了火星探測器從我國西昌衛(wèi)星發(fā)射中心發(fā)射、中途軌道修正并被火星捕獲的全過程。仿真結(jié)果可知：圓錐曲線拼接法是行星際軌道初步設(shè)計的基本方法，該法獲得的參數(shù)滿足行星際軌道初步設(shè)計的精度要求，可作為行星際軌道精確設(shè)計的良好初值；用B平面法作為中途軌道修正方法，可使搜索算法有良好的收斂性；組合使用基于普適變量法的蘭伯特問題數(shù)值求解和STK仿真軟件，可有效獲得真實(shí)力模型下的行星際軌道設(shè)計。本文的火星探測器軌道設(shè)計未考慮中途軌道修正策略，尚有待深入研究。

表1 仿真所得參數(shù)Tab.1 Parameter gained by simulation

圖6 地球逃逸軌道、地火轉(zhuǎn)移軌道和火星捕獲軌道仿真結(jié)果Fig.6 Simulation results of earth departure orbit，earth－Mars transfer trajectory and mars capture orbit

［1］ GAUSS D F.Theory of the motion of the heavenly bodies moving about the sun in conic sections［M］.New York：Dover Publications，1963.

［2］ BATTIN R H.An introduction to the mathematics and method of astrodynamics［M］.New York：AIAA Education Serious，1987.

［3］ BATTIN R H，VAUGHAN R M.An elegant lambert algorithm［J］.Journal of Guidance，Control and Dynamics，1984，7：662－670.

［4］ BATE R R，MUELLER D，WHITE J E.Fundamentals of astrodynamics［M］.New York：Dover Publications，1971.

［5］ BOND V，MODERN A.Astrodynamics：fundmentals and perturbation methods［M］.Princeton University Press，1996.

［6］ CURTIS H D.Orbital mechanics for engineering students［M］. Elsevier： Butterworth－Heinemann，2005.