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

(1. 宇航動力學國家重點實驗室，陜西 西安 710043；2.西安衛(wèi)星測控中心，陜西 西安 710043)

0 引言

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

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

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

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

1 濾波計算模型

1.1 濾波算法

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

UKF算法基本方法如下：

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

觀測量預測值為：

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

1.2 狀態(tài)模型

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

1.3 觀測模型

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

1.4 起步計算

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

2 自適應濾波計算

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

此處采用2種自適應處理方法，并對其應用效果進行探討。

2.1 基于機動檢測及模型切換的自適應處理

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

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

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

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

2.2 基于加速度模型補償?shù)奈粗獧C動自適應處理

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

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

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

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

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

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

3 計算及分析

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

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

圖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

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

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

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

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

4 結束語

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

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

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

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

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