基于外測(cè)觀測(cè)的再入飛行轉(zhuǎn)彎方向辨識(shí)方法

淡 鵬， 王 丹, 李恒年

(1. 宇航動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 陜西西安 710043； 2. 西安衛(wèi)星測(cè)控中心， 陜西西安 710043)

0 引 言

航天器再入返回[1-3]地球過程中，為使其著陸于特定的位置或區(qū)域，再入飛行時(shí)需要不斷調(diào)整飛行方向，使其按照事先設(shè)計(jì)好的制導(dǎo)方法向目標(biāo)點(diǎn)或目標(biāo)區(qū)域飛行。在一些機(jī)動(dòng)能力較強(qiáng)的航天器再入飛行中，有時(shí)需要對(duì)其飛行過程的轉(zhuǎn)彎方向進(jìn)行實(shí)時(shí)辨識(shí)，以便更好地判斷及預(yù)測(cè)其實(shí)際的飛行路徑情況(有時(shí)飛行器因?yàn)殡x軌制動(dòng)控制差異、大氣環(huán)境變化影響等而未能按照事先設(shè)定的路徑飛行)。

由于黑障[4-6]的影響，航天器再入飛行過程部分時(shí)段會(huì)出現(xiàn)遙測(cè)下傳中斷的現(xiàn)象，此時(shí)就無法從遙測(cè)數(shù)據(jù)中獲取有效的飛行姿態(tài)及氣動(dòng)參數(shù)信息，故而辨識(shí)不出其轉(zhuǎn)彎方向。這種情況下，來自于地面外測(cè)[1]設(shè)備的觀測(cè)數(shù)據(jù)就成為黑障中或姿態(tài)不穩(wěn)定過程中航天器飛行狀態(tài)估計(jì)的主要可依賴數(shù)據(jù)。另外在一些非合作目標(biāo)的再入觀測(cè)中，目前也主要依賴?yán)走_(dá)[7]等外測(cè)測(cè)量數(shù)據(jù)。

因缺乏姿態(tài)信息，航天器再入飛行過程中僅外測(cè)觀測(cè)下的轉(zhuǎn)彎方向辨識(shí)只能從外測(cè)數(shù)據(jù)中的測(cè)距、方位角、仰角、測(cè)速(測(cè)距變化率)等觀測(cè)量上入手，而這幾類數(shù)據(jù)運(yùn)動(dòng)特征不明顯，不能直接轉(zhuǎn)換出速度信息，從而給轉(zhuǎn)彎方向的辨識(shí)帶來一定困難。

關(guān)于再入飛行轉(zhuǎn)彎方向的實(shí)時(shí)辨識(shí)問題，經(jīng)查詢很難找到類似的解決案例或文獻(xiàn)作為參考。為此，針對(duì)部分航天器再入返回工程實(shí)施過程中遇到的這一實(shí)際問題，在對(duì)再入飛行過程受力分析的基礎(chǔ)上給出了一種轉(zhuǎn)彎方向的實(shí)時(shí)在線辨識(shí)方法。

1 轉(zhuǎn)彎方向辨識(shí)方法

1.1 再入過程受力情況

航天器再入段飛行過程主要受到地球引力、大氣阻力、氣動(dòng)升力的作用。而其他攝動(dòng)力的量級(jí)較小，在再入飛行過程中可忽略不計(jì)[8]。這幾個(gè)主要的作用力如圖1所示。

圖1 再入過程主要受力示意圖

僅考慮J2項(xiàng)時(shí)，地球引力加速度在地心矢徑r和地球自轉(zhuǎn)軸矢量we上的分量gr,gw分別為

(1)

式中：J2=1.082 63×10-3為二階帶諧項(xiàng)系數(shù)；μ=3.986 005×1014m3/s2為地球引力常數(shù)；Re為地球赤道半徑。

設(shè)ρ為大氣密度(kg/m3)，v是航天器相對(duì)大氣的速度，其矢量模為v(m/s)，Sref為飛行器參考面積(m2)，m為飛行器質(zhì)量(kg)，CL和CD為升力系數(shù)和阻力系數(shù)(無量綱)。

大氣阻力在地球固連系(與J2000慣性系可相互轉(zhuǎn)換)下的加速度為

(2)

氣動(dòng)升力在地球固連系下的加速度為

(3)

式中，C0為航天器所受總升力再入縱平面的側(cè)向分量在地固系下的單位矢量，L0為航天器總升力再入縱平面內(nèi)分量在地固系下的單位矢量,γ為傾側(cè)角。

從這幾個(gè)主要作用力的代數(shù)式可見，再入飛行過程中軌道面的法線方向主要受到氣動(dòng)阻力法向分量、氣動(dòng)升力法向分量、地球非球形攝動(dòng)力法向分量的影響。而實(shí)際計(jì)算中發(fā)現(xiàn)，地球非球形攝動(dòng)力的法線分量相對(duì)其他幾個(gè)量明顯較小[8]。

由此，航天器實(shí)際再入飛行過程中，通過對(duì)姿態(tài)角的控制，特別是在其發(fā)生滾轉(zhuǎn)以后，升力在當(dāng)?shù)劂U垂面和水平面上的分量會(huì)發(fā)生改變，從而使航天器具有了一定的縱向機(jī)動(dòng)和側(cè)向機(jī)動(dòng)能力，進(jìn)而可實(shí)現(xiàn)方向轉(zhuǎn)彎。

也就是說，通過實(shí)時(shí)計(jì)算當(dāng)前軌道平面法向的加速度或速度變化即可對(duì)滾轉(zhuǎn)過程發(fā)生的轉(zhuǎn)彎方向進(jìn)行辨識(shí)。

基于以上分析，本文對(duì)滾轉(zhuǎn)方向辨識(shí)的思路是：當(dāng)只有外測(cè)測(cè)量數(shù)據(jù)時(shí)，先通過濾波算法使用外測(cè)觀測(cè)量計(jì)算出飛行彈道的速度矢量，然后計(jì)算當(dāng)前速度與前一點(diǎn)的速度差(速度增量)在前一點(diǎn)速度法線方向的分量；進(jìn)而通過對(duì)一系列時(shí)間點(diǎn)處的法線方向速度分量的滑窗統(tǒng)計(jì)來獲取對(duì)航天器轉(zhuǎn)彎方向的判斷。

1.2 再入飛行外彈道計(jì)算模型

再入飛行過程外彈道測(cè)量數(shù)據(jù)(主要為測(cè)距、測(cè)速、測(cè)角)無法直接計(jì)算出速度參數(shù)，除幾何方法、差分方法外，可采用濾波算法來實(shí)現(xiàn)飛行彈道中速度參數(shù)的估計(jì)?？紤]到再入飛行時(shí)受力情況較為復(fù)雜，無姿態(tài)及其他氣動(dòng)參數(shù)遙測(cè)信息時(shí)(僅外測(cè)觀測(cè)下)，較難進(jìn)行飛行狀態(tài)的動(dòng)力學(xué)建模[1]，故此處考慮采用非動(dòng)力學(xué)建模算法進(jìn)行濾波計(jì)算。

在非動(dòng)力學(xué)建模方法中，“當(dāng)前”統(tǒng)計(jì)模型[9-11]是機(jī)動(dòng)目標(biāo)建模的一種較常用算法?！爱?dāng)前”統(tǒng)計(jì)模型把機(jī)動(dòng)目標(biāo)加速度的一步預(yù)測(cè)看作是“當(dāng)前”機(jī)動(dòng)加速度，并采用該加速度作為修正瑞利分布的均值來實(shí)現(xiàn)目標(biāo)自適應(yīng)跟蹤。由于采用非零均值的機(jī)動(dòng)加速度和相應(yīng)機(jī)動(dòng)加速度方差的自適應(yīng)調(diào)整，因此具有較好的機(jī)動(dòng)目標(biāo)跟蹤能力。

無跡卡爾曼濾波[12-14](UKF)是由牛津大學(xué)機(jī)器人技術(shù)研究所Julier和Uhlman等人于1995年提出的，后來被推廣到非線性動(dòng)力學(xué)系統(tǒng)參數(shù)估計(jì)中。UKF算法的核心是UT變換，它選擇2n+1個(gè)Sigma采樣點(diǎn)來逼近樣本非線性變換參量的矩，能達(dá)到真實(shí)值的三階精度。同時(shí)UKF算法不需要求導(dǎo)計(jì)算Jacobian矩陣，不需要清楚地了解非線性函數(shù)的具體形式，便于進(jìn)行模塊化開發(fā)，因而在工程中有較大的應(yīng)用價(jià)值。

基于以上考慮，本文將采用當(dāng)前統(tǒng)計(jì)模型及UKF濾波算法來實(shí)現(xiàn)僅外測(cè)觀測(cè)下飛行彈道中速度參數(shù)的估計(jì)。

“當(dāng)前”統(tǒng)計(jì)模型下的狀態(tài)外推模型為

X(k+1|k)=Φ(k+1|k)X(k|k)+

(4)

因x、y、z各方向相互正交，則有Φ=diag{Φx,Φy,Φz(mì)},U=diag{Ux,Uy,Uz}。

狀態(tài)噪聲協(xié)方差陣

Q=E{W(k)·WT(k)}=diag{Qx,Qy,Qz}

(5)

對(duì)x、y、z各方向，有

(6)

x、y、z三個(gè)方向上噪聲協(xié)方差陣Qi(i=x,y,z)的計(jì)算公式為

(7)

式中，

(8)

(9)

這樣就建立了一種外測(cè)觀測(cè)下的基于當(dāng)前統(tǒng)計(jì)模型的UKF濾波再入飛行彈道估計(jì)算法。

1.3 濾波觀測(cè)模型建立

由于外彈道測(cè)量數(shù)據(jù)中各類觀測(cè)數(shù)據(jù)的質(zhì)量不一，一般情況下測(cè)距數(shù)據(jù)質(zhì)量較高、而測(cè)角數(shù)據(jù)質(zhì)量較差，為此，濾波時(shí)將觀測(cè)模型建立在測(cè)站地平坐標(biāo)系下，以便充分利用不同精度的觀測(cè)量。則觀測(cè)方程為

(10)

1.4 濾波初始化方法

濾波計(jì)算時(shí)的一個(gè)重要內(nèi)容是起步初值的計(jì)算，考慮到再入過程中外測(cè)跟蹤數(shù)據(jù)建立在測(cè)站地平坐標(biāo)系下，無法直接轉(zhuǎn)換為J2000.0坐標(biāo)系的位置速度，但可以由測(cè)距和測(cè)角值計(jì)算出位置矢量，因此可通過多個(gè)點(diǎn)的曲線擬合完成起步計(jì)算。

取多個(gè)連續(xù)的觀測(cè)點(diǎn)，設(shè)第j個(gè)點(diǎn)處測(cè)量值時(shí)間、測(cè)距、方位角、仰角分別為tj，ρj,Aj,Ej。由測(cè)站坐標(biāo)可計(jì)算出測(cè)站東北天坐標(biāo)系到J2000.0地心慣性系的轉(zhuǎn)換矩陣為Mtj，則可計(jì)算出各點(diǎn)在測(cè)站東北天坐標(biāo)系下的位置矢量為

Rdj=[ρjcosEjsinAjρjcosEjcosAjρjsinEj]

j=0,1,2,…

(11)

則各點(diǎn)在J2000.0坐標(biāo)系下的位置矢量為

(12)

進(jìn)而對(duì)該位置數(shù)列使用最小二乘曲線擬合后求導(dǎo)的方法可計(jì)算出速度初值。

狀態(tài)方程中加速度等其他參數(shù)的起步初值可設(shè)置為0。這樣就得到了濾波的起步初值。

1.5 單點(diǎn)轉(zhuǎn)彎方向辨識(shí)

則前一時(shí)刻軌道面的法線方向單位矢量為

Nk=rk×vk/‖rk‖/‖vk‖

(13)

當(dāng)前時(shí)刻與前一時(shí)刻的速度差為

(14)

給定一個(gè)小值(定義為σ，σ≥0)，用于判斷轉(zhuǎn)彎方向是否為0。

則可定義轉(zhuǎn)彎方向的辨別函數(shù)為

(15)

當(dāng)F取值為1時(shí)，表示左轉(zhuǎn)(順著飛行速度方向觀察時(shí))；取值-1時(shí)，表示右轉(zhuǎn)；取值為0時(shí)，表示未轉(zhuǎn)彎(或值太小無法辨識(shí))。

這樣就得到了當(dāng)前時(shí)刻的轉(zhuǎn)彎方向的辨識(shí)值。

此方法是利用速度增量在軌道面法向上的分量進(jìn)行轉(zhuǎn)彎方向的辨識(shí)計(jì)算，考慮到轉(zhuǎn)彎時(shí)軌道面法線方向上的速度增量將直接在下一時(shí)刻的速度矢量上體現(xiàn)出來，因而上述計(jì)算過程也可直接簡(jiǎn)化為計(jì)算k+1時(shí)刻的速度vk+1在k時(shí)刻軌道面法線方向Nk上的投影大小的判斷。即用λk+1=vk+1·Nk代替ηk+1，并運(yùn)用該判別函數(shù)進(jìn)行轉(zhuǎn)彎方向的判斷。

1.6 多點(diǎn)滑窗統(tǒng)計(jì)

由于濾波計(jì)算時(shí)存在“震蕩”收斂問題，以及某點(diǎn)數(shù)據(jù)質(zhì)量較差[1]時(shí)可能影響單點(diǎn)判斷結(jié)果的準(zhǔn)確性。因而僅僅使用上面的單點(diǎn)轉(zhuǎn)彎方向辨識(shí)方法可能會(huì)出現(xiàn)錯(cuò)判的情況。為此，實(shí)際工程使用中，還不能直接使用上面的單點(diǎn)判別結(jié)果，而需要進(jìn)行多點(diǎn)數(shù)據(jù)的綜合判斷。

為此，在飛行過程的實(shí)時(shí)監(jiān)視計(jì)算中，可采用多點(diǎn)滑窗統(tǒng)計(jì)的方法來進(jìn)行轉(zhuǎn)彎方向的辨識(shí)。

設(shè)固定使用當(dāng)前時(shí)刻及其之前n-1個(gè)時(shí)間點(diǎn)上的單點(diǎn)辨別結(jié)果進(jìn)行統(tǒng)計(jì)計(jì)算，這n個(gè)時(shí)間點(diǎn)上的辨識(shí)結(jié)果分別為f1,f2,…,fn，則可統(tǒng)計(jì)出值為1的點(diǎn)數(shù)為n1，值為-1的點(diǎn)數(shù)為n-1。設(shè)綜合判斷閾值為ε(0.5<ε≤1，值越大越可靠，但判斷不出的可能性也將變大)，當(dāng)n1/n>ε時(shí)可認(rèn)為發(fā)生左向轉(zhuǎn)彎，當(dāng)n-1/n>ε時(shí)認(rèn)為發(fā)生右向轉(zhuǎn)彎，否則認(rèn)為無法辨識(shí)出轉(zhuǎn)彎方向(或未發(fā)生轉(zhuǎn)彎)。

2 計(jì)算及分析

為了驗(yàn)證算法的正確性，采用某航天器再入返回過程的理論彈道數(shù)據(jù)進(jìn)行仿真計(jì)算(該航天器再入飛行中的側(cè)滑角為小量，可忽略)，按多站接力方式生成仿真數(shù)據(jù)，并對(duì)各站的測(cè)距、測(cè)角、測(cè)速分別設(shè)置50 m、0.02°、0.05 m/s的隨機(jī)誤差。

使用本文算法實(shí)時(shí)計(jì)算轉(zhuǎn)彎方向，并與其在軌道系下的理論滾動(dòng)姿態(tài)角給出的轉(zhuǎn)彎方向進(jìn)行對(duì)比。

對(duì)軌道系下的滾動(dòng)姿態(tài)角(記為φ)，給定一個(gè)比較閾值σ2(σ2>0，此算例中取為0.01°)，結(jié)合滾動(dòng)角的定義，當(dāng)φ>σ2時(shí)，認(rèn)為右轉(zhuǎn)；當(dāng)φ<-σ2時(shí)，認(rèn)為左轉(zhuǎn)。

為了對(duì)比明顯，將兩種方法表示的左轉(zhuǎn)及右轉(zhuǎn)值進(jìn)行區(qū)分。

方法1為采用本文濾波及多點(diǎn)滑窗算法的判別結(jié)果，左轉(zhuǎn)辨識(shí)結(jié)果取值為1，右轉(zhuǎn)辨識(shí)結(jié)果取值為-1；方法2為滾動(dòng)姿態(tài)角判別出的轉(zhuǎn)彎方向，左轉(zhuǎn)取值為0.9，右轉(zhuǎn)取值為-0.9；無法判斷時(shí)兩種方法均取值為0。

使用這兩種方法計(jì)算的轉(zhuǎn)彎方向值對(duì)比如圖2所示。

圖2 兩種方法轉(zhuǎn)彎方向?qū)Ρ?/p>

需說明的是，因?yàn)樗抡娴耐鉁y(cè)測(cè)量量有些弧段沒有觀測(cè)值(設(shè)備跟蹤不上)，故圖2中部分時(shí)段只有理論姿態(tài)角的計(jì)算值，而無外測(cè)濾波計(jì)算值的判別結(jié)果(方法1)。

從圖中可見，兩種方法判別的轉(zhuǎn)彎方向基本上是一致的，這說明本文算法是可行的。由于濾波震蕩問題使其部分點(diǎn)存在誤判的情況，但大多數(shù)點(diǎn)均判別正確。另外，濾波輸出的判斷結(jié)果存在稍微滯后的現(xiàn)象(1～5 s，各處存在差異)。但總體上來看，本文算法是有效的，對(duì)轉(zhuǎn)彎方向判別的正確率較高。

造成濾波判斷稍滯后的原因在于：1)濾波計(jì)算結(jié)果相對(duì)于測(cè)量數(shù)據(jù)時(shí)標(biāo)的滯后；2)多點(diǎn)滑窗統(tǒng)計(jì)造成的統(tǒng)計(jì)滯后。第一種因素改善稍顯困難，需從提高數(shù)據(jù)的質(zhì)量，以及尋找新的濾波算法模型上下功夫；第二種可從減少滑窗時(shí)長(zhǎng)，或增加每秒時(shí)長(zhǎng)內(nèi)數(shù)據(jù)采樣個(gè)數(shù)等方面入手，在采樣頻率不變的情況下，減少時(shí)長(zhǎng)可能會(huì)使得誤判的幾率放大，所以需要權(quán)衡處理。因?yàn)楸疚牡幕八惴ㄊ菚r(shí)間上逐點(diǎn)滑動(dòng)的，所以算法本身不會(huì)帶來處理延遲。

需要說明的是，本文方法在轉(zhuǎn)彎方向辨識(shí)時(shí)未考慮側(cè)滑角的影響(當(dāng)側(cè)滑角很小時(shí))，這樣才可直接與滾動(dòng)角判斷出的轉(zhuǎn)彎方向比較。否則，在實(shí)際計(jì)算和比較過程中還需要考慮側(cè)滑角的影響。

3 結(jié)束語

針對(duì)僅外測(cè)觀測(cè)下的再入飛行過程中轉(zhuǎn)彎方向的實(shí)時(shí)辨識(shí)問題，提出了一種基于UKF濾波及速度增量(或速度)法向分解的解決方法。從仿真計(jì)算結(jié)果來看，本文算法是可行的。

在使用該方法時(shí)需要注意的是：對(duì)于飛行中轉(zhuǎn)彎方向的實(shí)時(shí)辨識(shí)所用的判斷閾值不應(yīng)太小，否則可能將極小作用力的加速度考慮進(jìn)去。也就是說，只有法線分量絕對(duì)值大于某一值后，才能認(rèn)為轉(zhuǎn)彎方向可判斷(可通過事先對(duì)理論飛行彈道分析的方法給定閾值)。

考慮到一些航天器外形比較復(fù)雜，致使返回過程的受力飛行分析更加困難，下一步將在考慮側(cè)滑角的轉(zhuǎn)彎方向辨識(shí)方面進(jìn)行研究。