制導(dǎo)率未知時(shí)返回彈道自適應(yīng)抗差濾波計(jì)算

淡 鵬，李恒年，麻蔚然，王 丹

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

0 引言

航天器返回過程彈道計(jì)算[1-2]是其落點(diǎn)預(yù)報(bào)的基礎(chǔ)，為落點(diǎn)預(yù)報(bào)提供初始位置速度，因而在返回任務(wù)實(shí)施過程中有著重要意義。外測跟蹤數(shù)據(jù)是返回過程的一類重要跟蹤數(shù)據(jù)，但其觀測量一般為測站地平坐標(biāo)系下的測距、方位角、仰角及測距變化率數(shù)據(jù)，要將其轉(zhuǎn)換出含有x，y，z三個(gè)分量的速度數(shù)據(jù)還需要采用濾波等其他較復(fù)雜的計(jì)算方法。

在航天器升力式返回地球過程中，為使其著陸于事先設(shè)定的區(qū)域，常常需要進(jìn)行制導(dǎo)控制[3-5]，不同的制導(dǎo)率設(shè)計(jì)對飛行彈道及落點(diǎn)位置的影響巨大，且一些航天器返回過程制導(dǎo)方法設(shè)計(jì)較為復(fù)雜，使其飛行過程受力情況變化明顯，給制導(dǎo)率未知時(shí)的返回彈道估計(jì)帶來困難。

針對此類不能準(zhǔn)確建立外推模型的機(jī)動[6]過程濾波問題，常用的計(jì)算方法是采用當(dāng)前統(tǒng)計(jì)模型、多項(xiàng)式模型等，文獻(xiàn)[7-9]對該2類模型進(jìn)行探討，但當(dāng)前統(tǒng)計(jì)及多項(xiàng)式模型對觀測數(shù)據(jù)的變化較為敏感，抗差性稍差。文獻(xiàn)[10]在橫向和縱向上建立再入目標(biāo)模型，并引入一種改進(jìn)的自適應(yīng)無跡卡爾曼濾波(Unscented Kalman Filter,UKF)算法，但該方法不適用于需要計(jì)算三維位置速度分量的情況；文獻(xiàn)[11]將重點(diǎn)放在了各種濾波方法使用上，對模型建立方法考慮較少；文獻(xiàn)[12]采用3站測距、測速及UKF濾波算法估計(jì)再入彈道，但不適用于少于3站測量的情況。

針對這些情況，在李恒年等人[13-14]提出的變質(zhì)量動力學(xué)模型基礎(chǔ)上，對其進(jìn)行擴(kuò)充和自適應(yīng)改進(jìn)處理，并采用UKF方法，實(shí)現(xiàn)對制導(dǎo)率未知狀態(tài)下的返回彈道位置速度三維分量的自適應(yīng)抗差估計(jì)計(jì)算。

1 濾波計(jì)算模型

1.1 濾波算法

目前常用的濾波計(jì)算框架有擴(kuò)展卡爾曼粒子濾波(Extended Kalman Particle Filter，EKPF)、UKF、容積卡爾曼濾波(Cubature Kalman Filter，CKF)和粒子濾波(PF)等，EKPF易于實(shí)現(xiàn)，但其線性化過程會引入模型誤差，且需計(jì)算非線性函數(shù)的Jacob矩陣；CKF采用等權(quán)值的一組采樣點(diǎn)進(jìn)行濾波計(jì)算，雖然計(jì)算量較UKF小，但實(shí)際計(jì)算表明，因各采樣點(diǎn)等權(quán)值，使得某些狀態(tài)下可能使非線性濾波的穩(wěn)定性減弱；PF實(shí)現(xiàn)稍顯復(fù)雜，故本文采用UKF算法。

UKF算法基本方法如下：

狀態(tài)估計(jì)協(xié)方差為(Q為狀態(tài)噪聲協(xié)方差陣)：

觀測量預(yù)測值為：

狀態(tài)向量與觀測向量之間相關(guān)協(xié)方差矩陣為：

1.2 狀態(tài)模型

以文獻(xiàn)[13]中的機(jī)動推力下的濾波狀態(tài)模型為基礎(chǔ)，建立系統(tǒng)狀態(tài)模型，由于原模型在計(jì)算加速度時(shí)，需要使用姿態(tài)信息，增加了算法的局限性。為此，此處將該模型中的加速度項(xiàng)進(jìn)行三維擴(kuò)展，在J2000慣性系下建立濾波系統(tǒng)狀態(tài)向量為：

1.3 觀測模型

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

1.4 起步計(jì)算

再入過程外測跟蹤數(shù)據(jù)建立在測站地平坐標(biāo)系下，無法直接轉(zhuǎn)換為J2000坐標(biāo)系下的位置速度，但可以由測距、測角值計(jì)算出位置矢量。據(jù)此，可用多個(gè)點(diǎn)進(jìn)行起步，并使用二次或三次多項(xiàng)式在x，y，z各分量上進(jìn)行多項(xiàng)式擬合，進(jìn)而由最小二乘算法計(jì)算出擬合點(diǎn)的位置分量變化(速度)值。狀態(tài)向量起步值中后4項(xiàng)可給定為0。

2 自適應(yīng)濾波計(jì)算

由于返回過程制導(dǎo)率未知，即存在未知的機(jī)動過程，導(dǎo)致動力學(xué)模型不能反映實(shí)際飛行狀態(tài)，為此，需要對濾波過程進(jìn)行自適應(yīng)處理。

此處采用2種自適應(yīng)處理方法，并對其應(yīng)用效果進(jìn)行探討。

2.1 基于機(jī)動檢測及模型切換的自適應(yīng)處理

機(jī)動檢測[15]的基本思想是，機(jī)動發(fā)生時(shí)目標(biāo)的狀態(tài)估計(jì)將出現(xiàn)偏差，導(dǎo)致濾波的新息(殘差)序列統(tǒng)計(jì)特性發(fā)生變化。為此，可根據(jù)濾波新息構(gòu)造二階統(tǒng)計(jì)特性，來實(shí)現(xiàn)機(jī)動的檢測。

對于新息序列，D(k)服從自由度為m的X2分布。當(dāng)機(jī)動發(fā)生時(shí)，新息序列的統(tǒng)計(jì)特性發(fā)生變化，不再是均值為零的高斯白噪聲，D(k)值可能出現(xiàn)較大的變化。為此，可根據(jù)返回彈道特性，對D(k)值設(shè)置一定的檢測門限，當(dāng)多個(gè)點(diǎn)的D(k)值連續(xù)超過該檢測門限時(shí)，認(rèn)為發(fā)生機(jī)動不采用單點(diǎn)檢測的原因在于觀測數(shù)據(jù)中可能存在野值，會影響判斷的準(zhǔn)確性。實(shí)現(xiàn)時(shí)，可采用多點(diǎn)滑窗方法對D(k)值序列進(jìn)行檢測。

當(dāng)檢測出機(jī)動發(fā)生時(shí)，通過調(diào)整濾波參數(shù)P和Q來完成不同模型的切換，使其濾波狀態(tài)快速適應(yīng)機(jī)動變化。

在連續(xù)多點(diǎn)D(k)值變化較大，且持續(xù)時(shí)間較長時(shí)，可通過重啟濾波器的方法來實(shí)現(xiàn)更快的收斂及對發(fā)散的抑制。

2.2 基于加速度模型補(bǔ)償?shù)奈粗獧C(jī)動自適應(yīng)處理

機(jī)動發(fā)生時(shí)，UKF的狀態(tài)模型將與實(shí)際情況不匹配，則在濾波計(jì)算中，若將Q矩陣設(shè)置為固定值，則可能與實(shí)際飛行變化出現(xiàn)較大偏差，導(dǎo)致濾波不能適應(yīng)機(jī)動變化，可能需要較長時(shí)間才能收斂，甚至?xí)霈F(xiàn)濾波發(fā)散情況。為此，考慮對濾波過程的狀態(tài)噪聲協(xié)方差矩陣Q進(jìn)行實(shí)時(shí)補(bǔ)償，來自適應(yīng)未知機(jī)動的發(fā)生。

考慮到機(jī)動發(fā)生時(shí)，加速度發(fā)生較大變化，狀態(tài)模型中的加速度項(xiàng)將與實(shí)際情況不符，因而對Q的實(shí)時(shí)補(bǔ)償將重點(diǎn)放在加速度分量的補(bǔ)償上。

加速度模型誤差近似值為：

進(jìn)而，對狀態(tài)噪聲協(xié)方差矩陣更新為：

式中，λ為彈道機(jī)動頻率系數(shù)(0<λ<1)，可取固定值(如0.001)。當(dāng)機(jī)動加速度小時(shí)，或需要提高抗差效果時(shí)，λ應(yīng)取小量。

這樣，即可由采樣點(diǎn)狀態(tài)量預(yù)測均值和濾波狀態(tài)量改進(jìn)值實(shí)現(xiàn)對狀態(tài)協(xié)方差矩陣的自動更新。

3 計(jì)算及分析

為了測試制導(dǎo)率未知狀態(tài)下返回彈道自適應(yīng)濾波算法的適應(yīng)性及應(yīng)用效果，分別采用當(dāng)前統(tǒng)計(jì)模型EKF算法、基于機(jī)動檢測及模型切換的自適應(yīng)濾波算法(此處記為MOUKF)，基于加速度模型補(bǔ)償?shù)奈粗獧C(jī)動自適應(yīng)濾波算法(記作COUKF)等3種算法進(jìn)行返回過程彈道濾波計(jì)算。以某次衛(wèi)星理論返回彈道為基礎(chǔ)，仿真多個(gè)連續(xù)接力的測站的外測跟蹤數(shù)據(jù)，并為測距、方位角、仰角和測距變化率分別加上10 m,0.01°，0.01°，0.01 m/s的隨機(jī)噪聲(1σ)。

3種算法計(jì)算的彈道近地點(diǎn)高度及理論曲線局部圖如圖1和圖2所示(該曲線相對于高度曲線，對濾波穩(wěn)定性及適應(yīng)性反應(yīng)更加明顯)。

圖1 Hp曲線第一段局部放大圖Fig.1 Partial enlarged drawing of the first section of Hp curve

圖2 Hp曲線第二段局部放大圖Fig.2 Partial enlarged drawing of the second section of Hp curve

由計(jì)算可知，MOUKF和COUKF均能夠自適應(yīng)返回段飛行彈道的變化，曲線震蕩幅度均小于EKF。但MOUKF對機(jī)動適應(yīng)能力沒有COUKF強(qiáng)，出現(xiàn)了多次曲線分段或重起步的現(xiàn)象。而COUKF在自適應(yīng)性、抗差性等方面均表現(xiàn)出了較好的適應(yīng)能力，且實(shí)現(xiàn)過程比MOUKF簡化。

為測試COUKF中加速度機(jī)動系數(shù)λ的影響，對λ分別取0.000 01和0.001，所得到的近地點(diǎn)高度曲線變化如圖3所示。

圖3 機(jī)動系數(shù)影響對比Fig.3 Comparison of the influence of maneuver coefficient

由圖3可以看出，λ取值小時(shí)，濾波震蕩幅度更小，抗差能力更強(qiáng)；但取值大時(shí)，對機(jī)動的反應(yīng)更加迅速。使用時(shí)可根據(jù)彈道特點(diǎn)合理取值。

4 結(jié)束語

針對外測跟蹤下的制導(dǎo)率未知狀態(tài)下返回彈道自適應(yīng)濾波問題，建立了包含10維數(shù)據(jù)的狀態(tài)模型，并分別采用基于機(jī)動檢測及模型切換的濾波算法、基于加速度模型補(bǔ)償?shù)奈粗獧C(jī)動濾波算法進(jìn)行了自適應(yīng)計(jì)算。從計(jì)算結(jié)果可以得出：

① 文中給出的濾波計(jì)算狀態(tài)模型及觀測模型是可行的，濾波算法雖然未采用制導(dǎo)率建立外推模型，但仍然能夠有效適應(yīng)返回彈道的飛行特點(diǎn)；

② 文中給出的基于機(jī)動檢測及模型切換的自適應(yīng)處理方法、基于加速度模型補(bǔ)償?shù)奈粗獧C(jī)動自適應(yīng)處理方法是可行的；

③ 基于加速度模型補(bǔ)償?shù)奈粗獧C(jī)動自適應(yīng)處理方法在對機(jī)動的適應(yīng)能力上優(yōu)于基于機(jī)動檢測及模型切換的算法，且比后者更易于實(shí)現(xiàn)；

應(yīng)該看到，制導(dǎo)率未知狀態(tài)下返回彈道的計(jì)算是一個(gè)較復(fù)雜的工程問題，實(shí)際飛行狀態(tài)及觀測數(shù)據(jù)質(zhì)量可能出現(xiàn)異常或較大偏離等情況。下一步將重點(diǎn)放在算法的健壯性及優(yōu)化上開展研究。