基于B平面算法的火星探測(cè)器中途軌道修正

唐 衛(wèi)

(上海航天控制技術(shù)研究所，上海200233)

基于B平面算法的火星探測(cè)器中途軌道修正

唐 衛(wèi)

(上海航天控制技術(shù)研究所，上海200233)

探測(cè)器在日心過(guò)渡段軌道的軌道修正問題是火星探測(cè)工程中需要解決的關(guān)鍵問題.首先，引入B平面坐標(biāo)系，推導(dǎo)了B平面參數(shù)的計(jì)算公式，驗(yàn)證了速度修正量與目標(biāo)參數(shù)間存在的線性化關(guān)系，給出了算法流程圖.其次，在引入合理的速度執(zhí)行誤差與測(cè)軌誤差的情況下，提出了一個(gè)包含三次主修正，一次機(jī)動(dòng)修正，一次備份修正的控制策略，仿真驗(yàn)證了控制策略的可行性.

火星探測(cè);軌道修正;B平面

逃離地球引力場(chǎng)后，星際探測(cè)器進(jìn)入日心過(guò)渡軌道.由于在逃逸引力場(chǎng)時(shí)不可避免存在變軌誤差和定軌誤差等，因此需要在日心過(guò)渡軌道段進(jìn)行軌道修正以滿足設(shè)定的要求[1].20世紀(jì)60年代，中途制導(dǎo)計(jì)算方法和制導(dǎo)策略的研究得到了許多學(xué)者的重視.數(shù)值搜索的方法計(jì)算量大，不利于軌道控制的分析和實(shí)際在軌操作，所以一般采用線性修正的方法.在小偏差情況下，線性方法具有精度足夠高，計(jì)算耗時(shí)少.Kizner發(fā)現(xiàn)建立在目標(biāo)星的B平面參數(shù)與軌道狀態(tài)量的偏差量之間存在著線性關(guān)系[2]，而B平面參數(shù)很容易轉(zhuǎn)換為目標(biāo)量，由此可計(jì)算得到修正誤差所需要的速度沖量，Kizner的工作大大促進(jìn)了深空軌道制導(dǎo)的研究，行星際飛行軌道修正的目標(biāo)參數(shù)(如軌道傾角，軌道高度)都能夠轉(zhuǎn)換為B平面參數(shù).Hintz和Chadwick基于B平面參數(shù)給出了一種中途制導(dǎo)軌道修正方法[3]，該方法用參數(shù)化數(shù)據(jù)來(lái)滿足目標(biāo)參數(shù)約束、定軌誤差等因素.本文采用B平面算法進(jìn)行中途軌道修正，并給出中途軌道控制策略.

1 B平面坐標(biāo)系

所謂B平面就是包含二體假設(shè)下的雙曲線軌跡的焦點(diǎn)(對(duì)于火星段雙曲線軌道，焦點(diǎn)就是火心)，且垂直于雙曲線的進(jìn)入漸近線的平面.B平面和探測(cè)器軌道平面相交，在空間得到一條交線.B矢量沿著這條線，以焦點(diǎn)為起點(diǎn)，終點(diǎn)在進(jìn)入漸近線穿過(guò)B平面處.設(shè)進(jìn)入漸近線單位矢量為S，引入行星赤道面法線或黃道面法線單位矢量N，二者叉乘得到垂直于S并穿過(guò)行星中心的矢量T:T=S×N，將S和T叉乘可得到矢量R:R=S×T，B平面內(nèi)B矢量在T和R方向投影，構(gòu)成了瞄準(zhǔn)參數(shù)B·T和B·R.其空間幾何關(guān)系如圖1所示.

圖1 B平面空間幾何Fig.1 Spatial geometry of B-plane

2 軌道控制算法

2.1 由位置矢量r和速度矢量ν計(jì)算B矢量

設(shè)其軌道面法線單位矢量n，則n=r×ν/|r× ν|，偏心率矢量為:e=ν×h/μ-r/|r|，其中μ為中心星體的引力位常數(shù)，h為探測(cè)器相對(duì)于火心的動(dòng)量矩，h=r×ν.由于S在探測(cè)器軌道平面內(nèi)，選取軌道平面內(nèi)兩個(gè)標(biāo)準(zhǔn)正交的矢量組成一個(gè)坐標(biāo)系，就能唯一確定它的位置，取e為X軸，e和n叉乘為Y軸.設(shè)S與e的夾角為θ，飛入漸近線矢量為:S= cosθe+sinθ(n×e)，則:B=b(S×n)，其中b為雙曲線短半軸.進(jìn)一步求得BT和BR:BT=B·T，BR= B·R.

2.2 由軌道要素計(jì)算B矢量

下面取B平面的參考平面為火星赤道平面，定義火星近焦點(diǎn)坐標(biāo)系如圖2所示.

圖2 火星近焦點(diǎn)坐標(biāo)系Fig.2 Mars near focus coordinate system

取偏心率單位矢量e為X軸，它以焦點(diǎn)為起點(diǎn)，背向近星體點(diǎn).n和e叉乘構(gòu)成Y軸，B矢量和S矢量分別可以表示為:B=b[sinθ -cosθ 0]T，S= [cosθ sinθ 0]T，其中θ=arcos(1/e)，e=1+點(diǎn)赤經(jīng)為Ω，近火點(diǎn)幅角為ω的目標(biāo)火星軌道，赤道面法向量N在近焦點(diǎn)坐標(biāo)系中為:N=[sinωsin i cosωsin i cos i]T，則T和R方向矢量可表示為:T= S×N，R=S×T，同理于2.1節(jié)有BT=B·T，BR= B·R，設(shè)下標(biāo)NORM表示相應(yīng)于標(biāo)準(zhǔn)軌道的值，則B平面的誤差散布為

2.3 B平面參數(shù)偏導(dǎo)矩陣的求法

已經(jīng)得到了B平面參數(shù)及B平面參數(shù)誤差的表達(dá)式，但是仍然需要建立B平面參數(shù)與需要修正的速度增量之間的關(guān)系.凱茲納發(fā)現(xiàn)目標(biāo)星B平面參數(shù)偏差與軌道狀態(tài)量偏差量之間存在很好的線性關(guān)系.下面給出這種線性關(guān)系的說(shuō)明.

首先，定義參數(shù)TOF表示從修正點(diǎn)到達(dá)B平面的飛行時(shí)間，這樣B參數(shù)實(shí)際上就由BT、BR和TOF構(gòu)成.定義B平面參數(shù)對(duì)于速度改變量的偏導(dǎo)數(shù).在火星赤道慣性坐標(biāo)系下，記V沿X軸方向的速度改變量為△VX，那么定義 BT對(duì) VX的偏導(dǎo)數(shù)如式(2)所示:

同理，BR、TOF對(duì)VX的偏導(dǎo)數(shù)以及BT、BR和TOF對(duì)VY和VZ的偏導(dǎo)數(shù)，由此得到B平面參數(shù)對(duì)V的偏導(dǎo)數(shù)矩陣Bp如式(3)所示:

考察偏導(dǎo)數(shù)矩陣Bp對(duì)△V的敏感度，將V沿三個(gè)坐標(biāo)軸改變量的取值范圍定為10-18km/s～101km/s，計(jì)算偏導(dǎo)數(shù)矩陣Bp中每個(gè)元素的值，以X軸方向?yàn)槔?，如圖3所示.

圖3 B平面參數(shù)對(duì)VX的偏導(dǎo)Fig.3 Partial derivative of B-plane parameters with respect to VX

從圖3可以看出，當(dāng)△V的X方向分量大致在10-12km/s～10-1km/s的量級(jí)范圍內(nèi)變化時(shí)，B參數(shù)對(duì)VX的偏導(dǎo)數(shù)基本為常數(shù)，Y、Z兩軸類似，因此Bp是常數(shù)矩陣，所以當(dāng)速度的改變量在上述范圍內(nèi)時(shí)，△B與△V有良好的線性關(guān)系，即有式(4):

速度修正量計(jì)算流程圖如下所示:

圖4 速度修正量計(jì)算流程圖Fig.4 Flow chart of velocity correction algorithm

由于 B平面需要以火星赤道慣性系作為參考平面，因此在計(jì)算實(shí)際軌道 B平面參數(shù)時(shí)，需要將當(dāng)前的位置矢量和速度矢量經(jīng)過(guò)坐標(biāo)系轉(zhuǎn)換，由日心黃道坐標(biāo)系轉(zhuǎn)到火星赤道慣性系下.圖5給出了各坐標(biāo)系之間的幾何關(guān)系:

(1)火星赤道坐標(biāo)系O-NYZ，N是X軸指向，XY坐標(biāo)面即火星赤道

(2)火星軌道坐標(biāo)系O-N'YZ，N'是X軸指向，XY坐標(biāo)面即火星軌道面

圖5 火星中心坐標(biāo)系示意圖Fig.5 Sketch of Mars center coordinate system

(3)日心黃道坐標(biāo)系O-γYZ，γ是X軸指向，XY坐標(biāo)面即黃道面

上述坐標(biāo)系中的O是火星的質(zhì)量中心，N是火星軌道與火星赤道的交點(diǎn)，相當(dāng)于地球的軌道(即黃道)與地球赤道的交點(diǎn)γ.N'是火星軌道與黃道的交點(diǎn).γ是黃道與地球赤道的交點(diǎn)，即春分點(diǎn).由此，探測(cè)器在日心黃道系到火星赤道系的坐標(biāo)旋轉(zhuǎn)矩陣為:RX(-ε)RZ(ψ)RX(I)RZ(Ω).4個(gè)旋轉(zhuǎn)角的計(jì)算公式如式(4)和式(5)所示:

其中d的單位是地球日，JD(t)是t時(shí)刻對(duì)應(yīng)的儒略日，JD(2000.0)是歷元時(shí)刻對(duì)應(yīng)的儒略日，其值為2451545.0，T即世紀(jì)數(shù).

3 控制策略與仿真分析

控制策略需要回答的問題是何時(shí)修正，每一次修正多少.關(guān)于修正時(shí)機(jī)的選擇，本文做簡(jiǎn)單的討論，并給出一個(gè)易于工程實(shí)現(xiàn)的策略.火星探測(cè)的軌道飛行時(shí)間很長(zhǎng)，對(duì)于一次修正后，其殘差沿軌道傳播仍有可能超出任務(wù)要求，所以一般要進(jìn)行多次.修正過(guò)早可能由于定軌精度不夠而使修正失去意義，修正過(guò)晚有可能因?yàn)檎`差傳播過(guò)大而超出修正能力[1].20世紀(jì)60年代星際制導(dǎo)研究快速發(fā)展，如何確定最優(yōu)修正次數(shù)和修正時(shí)刻以使修正燃料最少是研究的熱點(diǎn)，Pfeiffer研究了動(dòng)態(tài)規(guī)劃方法[4]，Breakwell提出的間距比(Space Ratio)方法[5]，該方法的主要思想是盡快進(jìn)行第一次修正，其后各次修正時(shí)間間隔與前一時(shí)間間隔之比為常數(shù)，直至滿足誤差要求，不再進(jìn)行修正.美國(guó)勇氣號(hào)(Spirit)和機(jī)遇號(hào)(Opportunity)火星探測(cè)器[6]的軌道修正策略中也充分重視了盡早修正這一原則.這里可以借鑒等間距比的思想，首先盡早修正.

本節(jié)結(jié)合本文涉及的火星探測(cè)器的實(shí)際情況，給出了三種修正方案，分別為修正一次，修正二次，修正三次，最終通過(guò)仿真說(shuō)明了，修正三次的方案合理且能夠?qū)崿F(xiàn)軌道修正的目的，最終給出了在三次修正基礎(chǔ)上余留一次機(jī)動(dòng)修正和一次備份修正的控制策略.結(jié)合火星探測(cè)器軌控發(fā)動(dòng)機(jī)性能、測(cè)軌精度、和姿控能力，估計(jì)中途修正速度大小執(zhí)行誤差± 0.5%，速度方向執(zhí)行誤差±0.1°，測(cè)軌精度上限為20～30km.目標(biāo)軌道軌道傾角 95°，近火點(diǎn)高度300km.

3.1 修正一次方案分析

選擇發(fā)射機(jī)會(huì)為2011年11月6日發(fā)射，考慮只修正一次的情況，表1給列出了在不同的修正時(shí)間進(jìn)行修正需要的速度增量以及相關(guān)參數(shù).

表1 不同時(shí)刻修正第一次的情況Tab.1 First correction at different times

表1中第一列對(duì)應(yīng)的是探測(cè)器出地球影響球的修正時(shí)刻.從表1可以看出隨著修正時(shí)刻的延后，所需要的速度修正量越來(lái)越大.通過(guò)仿真發(fā)現(xiàn)，用 B平面法試圖在第32天后任選一個(gè)時(shí)刻進(jìn)行修正時(shí)，算法無(wú)法收斂，這是由于計(jì)算出的|△V|的量級(jí)超過(guò)了10-1km/s，超出了偏導(dǎo)數(shù)矩陣Bp對(duì)△V的線性化范圍.因此，需要盡早進(jìn)行修正，本文將第一次修正的時(shí)間點(diǎn)放在探測(cè)器剛出地球影響球時(shí).相對(duì)于探測(cè)器實(shí)際的速度而言，第一次修正的速度改變量是比較小的，后續(xù)軌道修正就是針對(duì)其殘差的進(jìn)行的，量級(jí)會(huì)更小.實(shí)際情況中，需要考慮執(zhí)行誤差、測(cè)軌誤差等因素的影響，因此給出在考慮中途修正速度大小誤差、速度方向誤差和測(cè)軌誤差的情況下，修正后最終到達(dá)火星時(shí) B平面參數(shù)誤差的分布情況圖6所示.

圖6 第一次修正后B平面參數(shù)誤差分布Fig.6 B-plane parameters error distribution after first correction

由上述仿真實(shí)例可以分析得出:一次修正方案，在考慮了誤差后，最終探測(cè)器到達(dá)火星時(shí)距離目標(biāo)軌道偏移量較大，超出了誤差要求，這是由于一次修正時(shí)，修正時(shí)刻較早，微小的誤差在未來(lái)長(zhǎng)時(shí)間的軌道運(yùn)行中不斷積累，導(dǎo)致最終無(wú)法到達(dá)預(yù)定目標(biāo).因此，一次修正方案無(wú)法滿足任務(wù)要求，需要增加中途修正次數(shù).

3.2 修正兩次方案分析

根據(jù)控制策略，第一次修正時(shí)刻為探測(cè)器出地球影響球時(shí)，即距發(fā)射日期3.12天.第一次修正量:

△V=[7.30627 19.88104 -26.28166]T(m/s)，|△V|=33.75445(m/s)，考慮存在測(cè)軌誤差、速度執(zhí)行誤差(包括方向和大小)，取隨機(jī)誤差下的某一次實(shí)際速度改變量:△V=[7.32264 19.83375 -26.29379]T(m/s)，選擇不同的時(shí)刻進(jìn)行第二次修正，求得速度修正量，仿真結(jié)果如表2所示.

由表2可以得出結(jié)論:第二次修正的速度量的量級(jí)小于第一次修正的速度量級(jí)，且第二次修正的速度量值隨著時(shí)間間隔的增加而增加，但是增加的比較緩慢.如表2所示，第二次修正的速度修正量隨時(shí)間的增幅很有限，因此，第二次修正時(shí)間點(diǎn)的選擇限制可以放寬要求.本文設(shè)計(jì)在距第一次修正后的4個(gè)月進(jìn)行第二次修正，計(jì)算的修正量如表2所示，同理于第一次修正，同時(shí)考慮測(cè)軌誤差，速度執(zhí)行誤差，最終到達(dá)火星時(shí)B平面參數(shù)誤差的分布情況的如圖7所示.

表2 不同時(shí)刻第二次修正的情況Tab.2 Second correction at different times

圖7 第二次修正后B平面參數(shù)誤差分布圖Fig.7 B-plane parameters error distribution after second correction

從圖7中可以看出，第二次修正后，探測(cè)器最終到達(dá)火星時(shí)，B平面參數(shù)誤差大幅度減小，但是在BT方向上仍有250km的誤差，因此需要繼續(xù)增加修正次數(shù).

3.3 修正三次方案分析

當(dāng)探測(cè)器接近火星時(shí)，其入軌精度需要盡可能的高.因此在修正兩次方案的基礎(chǔ)上，增加一次探測(cè)器進(jìn)入火星影響球之前的修正，形成三次修正的方案.本文設(shè)計(jì)第三次修正距離第二次修正間隔175天，這樣可以保證第三次修正點(diǎn)距離探測(cè)器到達(dá)火星影響球還有3～5天的時(shí)間余量.

第二次修正量:△V=[-3.49298 3.44780 0.71025]T(m/s)，取隨機(jī)誤差影響下的某一次實(shí)際速度改變量:△V=[-3.4773 3.42703 0.70552]T(m/s)，計(jì)算得到第三次修正需要的速度增量△V=[-0.01431 -0.12062 0.02096]T(m/s)，最終到達(dá)火星時(shí)B平面參數(shù)誤差的分布情況的如圖8所示.

圖8 第三次修正后B平面參數(shù)誤差分布圖Fig.8 B-plane parameter error distribution after third correction

從圖8中可以看出，經(jīng)過(guò)三次修正，探測(cè)器到達(dá)火星影響球時(shí)，B平面參數(shù)誤差可以達(dá)到25km.如果需要增加最終入射精度，即滿足軌道傾角和近火點(diǎn)高度要求，可以增加一次機(jī)動(dòng)修正.

以探測(cè)器到達(dá)火星影響球時(shí)進(jìn)行一次機(jī)動(dòng)修正為例，經(jīng)過(guò)三次修正后，考慮誤差影響，隨機(jī)抽取探測(cè)器的一組狀態(tài)值，r=[0.433077 5.336899 2.1767830]T×105km，υ=[-0.20324 -2.55276

-1.00312]T(km/s)，對(duì)應(yīng)的軌道傾角為95.86°，近火點(diǎn)高度300.002km，可見由于考慮了實(shí)際的誤差，導(dǎo)致探測(cè)器進(jìn)入火星時(shí)不在預(yù)定軌道平面內(nèi)，可以進(jìn)行一次機(jī)動(dòng)修正，計(jì)算得所需速度修正量為:

△V=[-0.5724 0.0271 0.0556]T(m/s)，修正后軌道傾角為95.00°，近火點(diǎn)高度為300.000km.在此基礎(chǔ)上，可以再增加一次備份修正，目的是補(bǔ)救第三、四次修正可能存在的意外造成的修正失敗.綜上所述，控制策略如表3所示.

表3 控制策略Tab.1 Control strategy

4 結(jié) 論

1)當(dāng)速度的變化量在10-12km/s～10-1km/s的量級(jí)范圍內(nèi)變化時(shí)，B平面參數(shù)對(duì)速度的偏導(dǎo)數(shù)矩陣是常數(shù)陣，可以直接建立B參數(shù)變化量與速度變化量之間線性關(guān)系.

2)給出了一個(gè)以三次修正為主，包含一次機(jī)動(dòng)修正，一次備份修正的控制策略，通過(guò)仿真，驗(yàn)證了控制策略的可行性.

[1] 楊嘉墀.航天軌道動(dòng)力學(xué)與控制(下)[M].北京:宇航出版社，1995 Yang JC.Spacecraft orbital dynamics and control(volume two)[M].Beijing:Aerospace Press，1995

[2] Kizner W.A method of describingmiss distances for lunar and interplanetary trajectories[R].Jet Propulsion Laboratory，1959，NASA-CR-63483;REPT.-674 65N85813

[3] Hintz C R，Chadwich C.A design technique for trajectory correction maneuvers[J].The Journal of the Astronautical Sciences，1985，33(4):429-443

[4] Pfeiffer C G.A dynam ic programming analysis ofmultiple of guidance correction of trajectory[J].AIAA Journal，1965，3(9):1674-1681

[5] Breakwell J V.Fuel requirements for crude interplanetary guidance[J].Advances in Astronautical Science，1960，(5):53-65

[6] Christopher L P，Behzad R，Julie A K.Mars exp loration rovers propulsivemaneuver design[C].AIAA/AAS Astrodynamics Specialist Conference and Exhibit，Providence，Rhode Island，16-19 August，2004

B-Plane M ethod-Based M idcourse Trajectory Correction M aneuver for M ars Probe

TANGWei
(Shanghai Institute of Spaceflight Control Technology，Shanghai 200233，China)

It is trajectory correction maneuver(TCM)that is absolutely important and has to be solved as the key problem during the Earth-to-Mars trajectory in the Mars exploration.Firstly，the B plane reference frame based on the target planet and target parameters are introduced.Then formulas of B plane parameters are deduced.In succession，it is validated that linear relationship exits between target parameters and the velocity correction value.Meanwhile the corresponding arithmetic flow chart is proposed.In the next place，taking the rational velocity maneuver error and the orbitmeasurement error into account，a control strategy containing threemain TCMs，one contingent TCM and one backup TCM is proposed and is proved possible and feasible by the simulation eventually.

Mars exploration;trajectory correction maneuver;B-plane

V423

A

1674-1579(2012)06-0050-06

唐 衛(wèi)(1986—)，男，助理工程師，研究方向?yàn)楹教炱髯藨B(tài)、軌道控制.

2011-11-03

DO I:10.3969/j.issn.1674-1579.2012.06.011