王志文，張洪華

(北京控制工程研究所，北京100190)

0 引 言

月球上升段由3個階段構(gòu)成，如圖1所示.第一階段是100m的垂直上升，該階段所用時間很短，主要目的是使探測器離開月面，達到動力上升所需要的條件.第二階段是單軸旋轉(zhuǎn)，該階段要求在盡量短的時間里調(diào)整上升器的姿態(tài)，使推力方向的指向從垂直方向轉(zhuǎn)為動力上升段開始時所要求的最優(yōu)指向.第一和第二階段，一般采用開環(huán)制導(dǎo)律.第三階段是動力上升段，即主動段.以美國Altair月球著陸器為例，該階段要求上升器到達近月點約15km、遠月點約75km的橢圓軌道.該階段飛行時間較長，一般有300多秒，故需要設(shè)計閉環(huán)最優(yōu)控制使得該段所消耗的燃料最少[1-3].

圖1 月球上升器3個階段[1]Fig.1 Three stages for lunar ascent

本文使用線性協(xié)方差的方法研究第三階段制導(dǎo)，即主動段制導(dǎo)，在存在初始狀態(tài)偏差和系統(tǒng)參數(shù)偏差的情況下，終端高度、速度大小以及飛行路徑角的偏差.其中，主動段制導(dǎo)采用的制導(dǎo)律為PEG+保持制導(dǎo)[4].在用線性協(xié)方差分析PEG制導(dǎo)性能之前，首先需要對PEG進行線性化.文獻[4]使用了一種數(shù)值方法實現(xiàn)了對PEG的線性化處理，但這種方法運算量大，需要時間較長，且面臨著步長選擇的問題.本文提出一種解析的方法，可以得到PEG3 個制導(dǎo)參數(shù)的線性化顯式表達，相比而言可以大幅減少運算量和仿真時間.通過對PEG線性化處理，可以考察主動段制導(dǎo)存在初始狀態(tài)偏差(位置偏差和速度偏差)、系統(tǒng)參數(shù)偏差(比沖偏差、質(zhì)量偏差和質(zhì)量流量偏差等)以及執(zhí)行機構(gòu)偏差的情況下，終端時刻的制導(dǎo)偏差.

1 制導(dǎo)律簡介

如引言所述，本文采用PEG+保持制導(dǎo).PEG是一種預(yù)測校正的迭代制導(dǎo)律.

考慮運動方程

式中，F(xiàn)為推力，m為當(dāng)前質(zhì)量，g為當(dāng)前位置月球重力加速度，uf為推力方向單位矢量，且具有如下形式:

式中，λv為速度增量方向上單位矢量垂直于表示推力方向的變化.tλ為制導(dǎo)時間，選擇合適的數(shù)值，可使得推力引起的速度變化沿著λv方向.uf，三者之間的關(guān)系如圖2所示.PEG制導(dǎo)的主要問題是得到這3個制導(dǎo)參數(shù)，然后按照式(2)就可以計算任意時刻推力方向的單位矢量.

圖2 制導(dǎo)參數(shù)之間的關(guān)系Fig.2 Relations of guidance parameters

關(guān)于PEG制導(dǎo)的詳細(xì)推導(dǎo)過程這里不再贅述，讀者可以參考相關(guān)文獻[4-6].需要說明的是，PEG以速度增量vgo作為其獨立變量，使用預(yù)測校正的迭代方法產(chǎn)生滿足終端要求的vgo，在此基礎(chǔ)上得到PEG制導(dǎo)律3個制導(dǎo)參數(shù)λv，，tλ的顯式表達.

2 線性協(xié)方差分析

線性協(xié)方差分析作為一種誤差分析方法已經(jīng)被成功應(yīng)用于交會對接、地月軌跡、行星任務(wù)的GNC系統(tǒng)和導(dǎo)航分析.和蒙特卡洛方法相比，線性協(xié)方差方法速度快，一次運行就可以得到誤差的分析結(jié)果，并且可以隨意選擇誤差組合，分析不同誤差帶來的影響.

將線性協(xié)方差用于主動段制導(dǎo)誤差分析，首先面臨的一個困難是將PEG制導(dǎo)律線性化.如前所述，PEG制導(dǎo)律復(fù)雜，高度非線性，并且包含預(yù)測校正的迭代求解過程，所以將其線性化比較困難.本文給出一種近似的PEG線性化解析方法，得到了PEG3 個制導(dǎo)參數(shù)的線性化顯式表達.

制導(dǎo)在tgo=tf時轉(zhuǎn)入保持模式，從蒙特卡洛仿真的角度來看，不同樣本從PEG轉(zhuǎn)入保持模式的時間都是不一樣的，所以針對這種基于事件觸發(fā)而不是時間觸發(fā)的制導(dǎo)模式轉(zhuǎn)換，線性協(xié)方差分析需要進行相應(yīng)的修正，這也就是文獻[4]中提到的事件觸發(fā)的問題.

最后，當(dāng)制導(dǎo)轉(zhuǎn)入保持模式后，推力單位矢量保持不變，但這個時候推力單位矢量的偏導(dǎo)并不為零，也并不等于從PEG轉(zhuǎn)入保持模式時刻的偏導(dǎo)值，即線性協(xié)方差分析還需要解決保持模式中推力方向偏導(dǎo)問題[4].

下面將根據(jù)這3個問題依次進行說明.

2.1 系統(tǒng)狀態(tài)定義

影響月球上升段主動段制導(dǎo)精度的有:初始位置偏差δr，初始速度偏差δv，比沖偏差δIsp，初始質(zhì)量偏差δm，初始質(zhì)量流量偏差δdm，推力大小偏差δF，以及執(zhí)行機構(gòu)偏差.其中執(zhí)行機構(gòu)參數(shù)偏差包括零偏誤差b，標(biāo)度因數(shù)誤差f，失準(zhǔn)誤差κ.此外對制導(dǎo)精度有影響的還有重力模型偏差εgrav，以及環(huán)境噪聲等.

設(shè)系統(tǒng)狀態(tài)為

2.2 PEG線性化

PEG制導(dǎo)律線性化中所有的偏差都是相對于標(biāo)稱軌跡來說的.標(biāo)稱軌跡是由理想初始條件并忽略所有中間過程誤差生成的.規(guī)定下面符號均表示x的標(biāo)稱值.

PEG制導(dǎo)律線性化指的是線性化PEG制導(dǎo)律的3個制導(dǎo)參數(shù)由PEG的制導(dǎo)方程可知

線性化式(3)和(4)有

由上述J，L的定義可知，其偏差δJ，δL可以寫成關(guān)于的線性表達，如下所示:

綜上可以看出要得到δλv，δtλ，只需要得到δvgo就可以了.

線性化式(5)有

同理δS，δQ可寫為如下的線性表達:

式中，vd為期望的終端速度，為當(dāng)前的速度.

線性化有

經(jīng)過推導(dǎo)，δvgo可以寫為關(guān)于的線性表達

近似認(rèn)為

結(jié)合式(8)、(11)、(12)可知，δvgo可以進一步寫為如下的線性表達:

至此，可以得到δvgo的最簡線性化表達

對于rgo的線性化δrgo，思路和上面類似.需要注意的是在PEG制導(dǎo)方程中，rgo的計算分兩步:第一步按著計算，然后利用約束修正rgo.所以在線性化rgo的時候也應(yīng)該按著這樣的兩步線性化.最終，δrgo也可以寫為如下的最簡線性化表達:

得到了vgo，rgo的線性化表達，代入就可以得到PEG3個制導(dǎo)參數(shù)的線性化表達

2.3 事件觸發(fā)

如前所述，設(shè)定PEG轉(zhuǎn)出的時間閥值tf，當(dāng)tgo=tf時制導(dǎo)轉(zhuǎn)入保持模式.由于這種模式轉(zhuǎn)換不是基于固定的時間，而是基于離散的事件(tgo=tf)，所以線性協(xié)方差需要進行相應(yīng)的修正[4].

關(guān)于包含事件觸發(fā)的線性協(xié)方差分析，文獻[7]中給出了較為詳細(xì)的推導(dǎo).本文將其原理應(yīng)用于主動段模式轉(zhuǎn)換時線性協(xié)方差分析.對一般情況，事件觸發(fā)可以用數(shù)學(xué)表達式表示成如下關(guān)于狀態(tài)的函數(shù)[7]:

式中，x表示系統(tǒng)狀態(tài)，當(dāng)?shù)仁綕M足時意味著該事件在te時刻發(fā)生了.在本文中，該事件可以具體表示為

這里省略具體的推導(dǎo)過程，只給出最后的結(jié)論[7]

2.4 保持模式中單位推力的偏導(dǎo)

當(dāng)tgo=tf時制導(dǎo)轉(zhuǎn)入保持模式，在這個階段，系統(tǒng)將保持PEG制導(dǎo)最后一步的推力方向不變，直到tf時間后主發(fā)動機關(guān)機，制導(dǎo)結(jié)束.

雖然在保持模式中推力方向保持不變，但推力方向的偏導(dǎo)并不是零.這是因為，從蒙特卡洛仿真的角度分析，不同的樣本在進入保持模式時推力方向不同，而線性協(xié)方差中推力方向的偏導(dǎo)統(tǒng)計的就是不同樣本相對于標(biāo)稱軌跡的偏差.而且保持模式中推力方向的偏導(dǎo)也不等于PEG最后時刻(進入保持模式之前)的偏導(dǎo)，否則就意味著保持模式中推力方向的偏差是由該階段系統(tǒng)狀態(tài)的變化所導(dǎo)致的[4].

文獻[4]給出了一種解決方法，即引入一個虛擬狀態(tài)uf，uf表示推力單位矢量.這時系統(tǒng)的狀態(tài)重寫為

uf由PEG制導(dǎo)得到，每調(diào)用一次PEG，就更新一次uf，所以可以認(rèn)為uf的動力學(xué)為[4]

uf的更新方程如下[4]:

于是可以得到更新后系統(tǒng)協(xié)方差矩陣

式(20)表示每調(diào)用一次PEG，系統(tǒng)協(xié)方差矩陣就更新一次，直到保持階段，不再調(diào)用PEG.

2.5 系統(tǒng)誤差源建模

本文將執(zhí)行機構(gòu)偏差中包含的f，κ，b以及重力模型偏差εgrav均視為一階馬爾科夫過程[8].

對每個執(zhí)行機構(gòu)參數(shù)pj，其狀態(tài)協(xié)方差為時間常數(shù)為τj，白噪聲方差如下:

其動力學(xué)方程可寫為

對重力模型偏差εgrav，其狀態(tài)協(xié)方差為相關(guān)距離為dk，白噪聲方差如下:

其動力學(xué)方程可寫成如下的分量形式:

系統(tǒng)參數(shù)認(rèn)為保持不變，故

3 仿真及分析

上升器主動段終端關(guān)注高度，速度大小，以及飛行路徑角的大小.本文按著PEG+保持的制導(dǎo)模式(tf=8s)，用上述線性協(xié)方差方法分析在存在初始狀態(tài)偏差、系統(tǒng)參數(shù)偏差以及執(zhí)行機構(gòu)偏差的情況下，終端高度、速度大小、飛行路徑角的偏差.為了驗證上述線性協(xié)方差的有效性，引入了相同初始條件下的300次蒙特卡洛仿真.仿真中用到的執(zhí)行機構(gòu)參數(shù)偏差和噪聲如表1[9]所示.

表1 仿真參數(shù)Tab.1 Simulation parameters

按著初始仿真條件，分為如下兩種情況:

其中情況2相對于情況1，增加了參數(shù)vex，m，dm的偏差，以考察系統(tǒng)參數(shù)偏差對PEG制導(dǎo)性能的影響.

圖3和圖4給出了兩種不同初始條件下線性協(xié)方差(LC)和蒙特卡洛(MC)關(guān)于高度偏差δh，速度大小偏差δv，以及飛行路徑角偏差δγ的對比曲線.從圖中可以看出，在305s之前有一個突變，這是由蒙特卡洛在處理從PEG轉(zhuǎn)入保持模式的數(shù)據(jù)時重置不同樣本的仿真時間為一個共同的仿真時間e(e為標(biāo)稱情況下從PEG轉(zhuǎn)入保持模式的時間)所造成的.這樣處理的結(jié)果方便MC與LC進行直接比較分析.同時從圖中還可以看出，經(jīng)過8s的保持模式，在制導(dǎo)結(jié)束時刻，LC和MC的結(jié)果重合較好，且都趨向于零.說明在存在上述偏差的情況下，PEG+保持模式的制導(dǎo)模式能夠保證終端精度要求.

從圖中可以看出，LC和MC的仿真曲線在過程中間時刻有一些偏差，這主要是由于本文提出的線性協(xié)方差分析是一種近似的分析，本身PEG線性化存在誤差，但實際中更關(guān)注的是終端時刻的狀態(tài)偏差，從圖中兩種情況的仿真分析可以看出，終端時刻LC和MC兩條仿真曲線基本重合，表明上述的線性協(xié)方差方法可以較好地估計終端時刻的狀態(tài)偏差.

圖3 情況1仿真圖Fig.3 Simulation results in Condition one

在存在初始狀態(tài)偏差，系統(tǒng)參數(shù)偏差以及執(zhí)行機構(gòu)偏差的情況下，制導(dǎo)終端狀態(tài)偏差都趨于零，可以看出PEG具有較好的魯棒性.然而這種分析其實暗含了一個前提，就是系統(tǒng)的這些初始參數(shù)偏差，狀態(tài)偏差都必須是已知的，也就是“可觀的”.換句話說，制導(dǎo)計算機中存儲的初始狀態(tài)以及系統(tǒng)參數(shù)必須和實際上升器初始時刻的狀態(tài)和參數(shù)一致.在實際應(yīng)用中，上述兩部分對應(yīng)的參數(shù)通常不能保證一致，所以需要在線估計實際上升器的參數(shù)，并把辨識得到的參數(shù)傳給制導(dǎo)計算機來產(chǎn)生制導(dǎo)指令，這樣就可以保證制導(dǎo)計算機中的上升器參數(shù)和實際上升器的參數(shù)保持一致.

圖4 情況2仿真圖Fig.4 Simulation results in condition 2

4 結(jié) 論

本文提出了一種解析的PEG線性化方法，應(yīng)用于月球上升器主動段誤差分析.研究了主動段在存在初始狀態(tài)偏差和系統(tǒng)參數(shù)偏差以及執(zhí)行機構(gòu)偏差的情況下，終端時刻的高度偏差、速度大小偏差、以及飛行路徑角偏差.最后通過和蒙特卡洛的仿真對比，驗證了本文提出的PEG線性化方法的有效性，可以看出在系統(tǒng)參數(shù)“可觀”的情況下，PEG具有較好的魯棒性.同時需要說明的是，本文提出的PEG線性化分析沒有考慮導(dǎo)航誤差以及姿態(tài)控制方面的影響，這也是后續(xù)研究中需要著重考慮的問題.

[1]KOS L D，POLSGROVE T T，SOSTARIC R R，et al.Altair descent and ascent reference trajectory design and initial dispersion analyes[R].20100035768，2010.

[2]SOSTARIC R R，MERRIAM R S.Lunar ascent and rendezvous trajectory design[R].20080009584，2008.

[3]BENNETT F V.Apollo experience reportmission planning for lunar module descent and ascent[R].19720018205，1972.

[4]ROSE M B，GELLER D.Linear covariance techniques for powered ascent[R].AIAA20108175，2010.

[5]MCHENRY R L， BRAND T J， LONG A D， et al， Space shuttle ascent guidance， navigation， and control[J].Journal of Astronautical Sciences， 1979，28(1):1-38.

[6]SPRINGMANN P，PROULX R，F(xiàn)ILL T.Lunar descent using sequential engine shutdown[R].2006-6678，2006.

[7]GELLER D K，ROSE M B，WOFFINDEN D C.Event trigger in linear covariance analysis with application to orbital rendezvous[J].Journal of Guidance， Control，and Dynamics， 2009，32(1):102-111.

[8]GELLER D K，CHRISTENSEN D P.Linear covariance analysis for powered lunar descent and landing[J].Journal of Spacecraft and Rockets，2009，46(6):1231-1235.

[9]MOESSER T J.Guidance and navigation linear covariance analysis for lunar powered descent[D].Logan Utah State:Utah State University，2010.