一種重尾量測噪聲下的高超聲速飛行器跟蹤算法

惲鵬 李星秀 吳盤龍 何山

摘要：為了解決在實際系統(tǒng)中因野值干擾帶來的高超聲速飛行器跟蹤精度下降的問題，提出了一種交互多模型變分貝葉斯濾波算法（IMM-VB），該算法通過子模型權重與馬爾可夫轉移矩陣獲取子模型的狀態(tài)預測值。隨后采用具有重尾特性的學生t分布取代高斯分布來描述量測模型，并利用VB算法實現(xiàn)子模型的量測協(xié)方差與狀態(tài)的聯(lián)合估計。最后在交互式多模型（IMM）框架下更新子模型權重與目標狀態(tài)的融合輸出。仿真結果表明，在野值觀測條件下該算法比IMM算法具有更高的跟蹤精度。

關鍵詞：高超聲速飛行器;野值;學生七分布;重尾分布;IMM-VB

中圖分類號：TP24 文獻標志碼：A

高超聲速飛行器具有飛行速度快、機動靈活、突防能力強、有效載荷大等優(yōu)點，因此受到了世界各國的關注。目前而言，交互式多模型（IMM）算法是解決高超聲速飛行器跟蹤的有效途徑之一。然而該算法是在廣義偽貝葉斯理論基礎上形成的，其前提是量測噪聲與過程噪聲均為高斯分布。在實際的高超聲速飛行器跟蹤過程中，由于大氣、云層等因素的干擾，噪聲的統(tǒng)計特性往往不符合高斯分布[1～3]，導致IMM算法的跟蹤精度較差，無法滿足實際系統(tǒng)的需求。

為了提高實際工程中的目標跟蹤精度，許多研究學者對當前的跟蹤算法進行了改進。為了解決重尾噪聲線性系統(tǒng)的濾波問題，參考文獻[4]提出了一種學生t濾波，該濾波基于噪聲的重尾特性，將噪聲建模為學生t分布，有效地克服了野值的干擾，提高了算法的魯棒性。然而當過程噪聲與量測噪聲的自由度參數(shù)不同時，該濾波無法實現(xiàn)狀態(tài)的更新。為了處理野值干擾導致的非高斯重尾噪聲下的目標跟蹤問題，參考文獻[5]提出了一種最大熵卡爾曼濾波器，然而在誤差協(xié)方差估計中缺少足夠的理論基礎，導致該濾波器的跟蹤精度不佳。參考文獻[6]利用變分貝葉斯法聯(lián)合估計目標的狀態(tài)分布與噪聲分布，有效地提高了目標的跟蹤精度，然而該算法是基于單模型設計的，因此不適合跟蹤具有強機動性的目標，如高超聲速飛行器等。

IMM因其計算量小、跟蹤精度高等優(yōu)點被廣泛應用于各類飛行器的跟蹤。參考文獻[7]將常加速度模型（CA）/勻速模型（CV）/SINGER模型相結合，并基于徑向速度估計補償技術，有效地提高了臨近空間高超聲速飛行器的跟蹤精度。參考文獻[8]將IMM與容積卡爾曼濾波[9]（CKF）相結合，提出了一種IMM-CKF算法，用以提高高速滑翔飛行器的跟蹤精度。參考文獻[10]將標簽伯努利[11]（ LMB）與無跡卡爾曼濾波[12]（UKF）相結合，提出了一種廣義LMB-UKF算法，實現(xiàn)了多個高超聲速飛行器的檢測與跟蹤。然而，上述算法將噪聲均建模為高斯分布，當噪聲為非高斯分布或者噪聲協(xié)方差與建模時給定數(shù)值差距過大時[13～15]，都將嚴重降低算法的跟蹤精度。因此不適用于野值干擾下高超聲速目標的跟蹤。

本文提出了一種IMM-VB算法用來解決野值觀測條件下的高超聲速目標跟蹤問題。首先將CA模型與當前統(tǒng)計模型（CS）[16]進行交互，獲取相應的預測值與權重，隨后基于噪聲的重尾特性，將量測噪聲建模為學生t分布，利用VB算法，實現(xiàn)狀態(tài)與量測噪聲的聯(lián)合估計。最后利用IMM-VB算法和IMM算法對高超聲速飛行器進行跟蹤。仿真結果表明，IMM-VB算法比IMM算法具有更高的跟蹤精度。

1 基礎算法

1.1 IMM算法

IMM算法是一種基于馬爾可夫轉移矩陣的廣義偽貝葉斯算法[17，18]。該算法根據(jù)目標的先驗信息，首先給予相應子模型的狀態(tài)分量與權重，隨后，對子模型進行濾波，同時根據(jù)量測信息與基于各個模型更新后的狀態(tài)對不同模型賦予相應的權重，并通過加權求和更新目標狀態(tài)。針對機動性較為復雜的目標，IMM算法相較于單模型算法具有更好的跟蹤精度。下面介紹一下IMM算法的4個步驟：狀態(tài)交互、狀態(tài)濾波、模型概率更新和狀態(tài)融合輸出。

1.1.1 狀態(tài)交互

假設k時刻各個子模型的目標狀態(tài)、協(xié)方差、權重分別為X_j（k|k）、P_j（k|k）、μ（k|k），j為子模型的編號，對目標各個子模型進行交互運算，獲取各個子模型的狀態(tài)X_j⁰（k|k）與協(xié)方差P_j⁰（k|k）：式中：N為模型總數(shù)，X_j（k|k）為模型i到模型j的轉移權重，π_ij為模型i到模型j的馬爾科夫轉移概率。

1.1.2 狀態(tài)濾波

狀態(tài)濾波過程與單模型濾波過程相同，依據(jù)卡爾曼濾波即可實現(xiàn)子模型下的狀態(tài)更新。將X_j⁰（k|k），P_j⁰（k|k）作為子模型的輸入進行狀態(tài)預測與更新：式中：X_j（k+1|k）為狀態(tài)預測值P_j（k+1|k）為預測協(xié)方差，K_j（k+1）為濾波增益，X_j（k+1|k+1）t狀態(tài)更新值，P_j（k+1|k+1）為協(xié)方差更新值，F(xiàn)_j（k）為狀態(tài)轉移矩陣，H_j（k+1）為量測矩陣Q_j（k）為過程噪聲協(xié)方差，R_j（k+1）為量測噪聲協(xié)方差。

1.1.3 模型概率更新

通過濾波過程中獲取的子模型信息與信息協(xié)方差，構建量測與模型匹配的似然函數(shù)∧_j（k+1）：式中：m為量測信息中的元素數(shù)目，子模型信息v_j（k+1）與信息協(xié)方差S_j（k+1）分別為：

通過狀態(tài)交互后的子模型權重以及相應的似然函數(shù)更新模型概率μ_j（k+1|k+1）：

1.1.4 狀態(tài)融合輸出

對各子模型的模型概率與狀態(tài)信息進行加權求和，獲取融合后的狀態(tài)X（k+1|k+1）與協(xié)方差P（k+1|k+1）：

1.2 VB算法

當實際的噪聲不滿足高斯分布時，貝葉斯公式很難推斷出目標的后驗分布。而VB算法通過迭代逼近的手段可以獲取近似的后驗分布[19]，有效地解決了非高斯條件下狀態(tài)后驗分布求解的問題。為了衡量近似分布的逼近程度，需要利用Kullback-Leibler散度（KLD）函數(shù)：式中：KL（q（Ψ））≥0且僅q（Ψ）=p（Ψ|X），p（Ψ|X）為參數(shù)集的后驗分布，q（Ψ）為參數(shù)集的近似分布，且KL（q（Ψ））越小，表示q（Ψ）的逼近程度越高。令下界函數(shù)L（q（Ψ））=

P（Ψ，X）表示Ψ和X的聯(lián)合概率密度函數(shù)，由式（16）可以看出，若想使KL（q（Ψ））最小化，則需要使L（q（Ψ））越大越好。基于平均場理論可以求得滿足L（q（Ψ））最大化的參數(shù)近似分布：式中：E_Ψi-[1np（Ψ，X）]表示除參數(shù)Ψ_i外1np（Ψ，X）的數(shù)學期望，Ψ_i為參數(shù)集Ψ的第i個參數(shù)。在對式（17）進行求解的過程中，需要使用到上一次迭代步驟中的參數(shù)分布，因此，求解式（17）是一個不斷迭代更新的過程，此時可以設置相應的迭代步數(shù)或者對相鄰迭代的下界函數(shù)差值設置閾值來終止迭代過程。最后迭代過程獲取的參數(shù)分布即為VB求解出的近似分布。

2 IMM-VB算法

由于量測信息存在野值的干擾，導致IMM算法的跟蹤精度下降，本文將IMM與VB相結合提出一種IMM-VB算法，可以有效地克服跟蹤過程中的野值干擾，提高高超聲速飛行器的跟蹤精度，其算法流程如圖1所示。

該算法本質上還是基于IMM思想進行，其狀態(tài)交互過程與IMM算法相同，這里就不再闡述。在濾波過程中由于量測信息存在野值干擾，導致其統(tǒng)計特性呈重尾分布，因此將量測噪聲建模為學生t分布，可得子模型的似然函數(shù)式中： St（;μ，Σ，v）表示均值為μ、協(xié)方差為Σ、自由度參數(shù)為v的學生t分布，而學生t分布是高斯分布與伽馬分布乘積的積分和，因此可以將式（18）改寫為：式中：N（;μ，Σ）表示均值為μ，協(xié)方差為Σ的高斯分布，G（;α，β）表示形狀參數(shù)α、尺度參數(shù)為β的伽馬分布。引入輔助隨機變量λ（k+1），并將其建模為伽馬分布，則可以從式（19）獲取概率密度函數(shù)p（λ_j（k+1））和p（z（k+1）|X_j（k+1），λ_j（k+1））：

利用式（20）、式（21）可以獲取聯(lián)合概率密度函數(shù)

結合式（4）、式（5）、式（20）、式（21）可以計算出聯(lián)合概率密度函數(shù)p（λ_j（k+1），X_j（k+1），z_1：k+1）：

利用式（23）可求取1np（λ_j（k+1），X_j（k+1），z_1：k+1）式中：n為狀態(tài)信息中的元素數(shù)目，利用式（17）、式（24）迭代獲取相應參數(shù)的近似分布：式中：C為保證qⁱ⁺¹（λ_j（k+1》積分為1的常量，后面的C均指代這樣的常量，但是對于每一個公式而言，其值可能不同。由式（27）可得：同理，可獲取狀態(tài)的近似分布：

由式（32）可推導出：

由高斯乘積定理可得：

進行完子模型濾波后，在模型概率更新過程中，由于輔助隨機變量的存在，導致子模型的新息協(xié)方差發(fā)生改變，此時似然函數(shù)∧_j（k + 1）為：

最后，隨后通過式（12）～式（14）、式（35）、式（36）、式（39）可以獲取融合后的目標狀態(tài)信息。

3 仿真試驗

本文對高超聲速飛行器[20]跟蹤進行了仿真試驗。以CA、CS為交互子模型，分別設計了相應的IMM算法和IMM-VB算法對該目標進行跟蹤。其中目標的狀態(tài)由二維坐標中的位置、速度、加速度組成：

馬爾可夫轉移矩陣為：

量測噪聲滿足：

式（42）表明野值干擾概率為0.1，且給出了有野值和無野值下的量測噪聲分布。過程噪聲滿足：

當子模型為CA時：當子模型為CS時：式中：采樣時間T=1s，a_max為目標加速度最大絕對值，a（k）為上一時刻目標加速度估計值，目標機動頻率α取值1/60。輔助隨機變量的自由度參數(shù)v取值50，初始時刻子模型的權重均為0.5，目標初始狀態(tài)和協(xié)方差分別為：

根據(jù)所給的條件設計相應的算法，獲取試驗信息。圖2展示了高超聲速飛行器的真實軌跡和兩種算法的估計軌跡。

圖3、圖4展示了兩種算法在不同坐標軸上的位置跟蹤均方根誤差，圖5、圖6展示了兩種算法不同坐標軸上的速度跟蹤均方根誤差?？梢钥闯?，IMM-VB相較于IMM算法具有更高位置和速度跟蹤精度。

4 結束語

針對野值干擾下的高超聲速目標跟蹤問題，本文提出了一種IMM-VB算法，并通過仿真試驗驗證了該算法相較于IMM算法更適用于跟蹤野值干擾下的高超聲速目標。然而由于VB的存在，導致算法的計算時間過長，因此后續(xù)的研究過程中，應考慮如何降低算法的復雜度。另外，該算法只針對于單目標，因此如何設計針對于多目標跟蹤的算法，也是后續(xù)工作的研究重點。

參考文獻

[1]Agamennoni G，Nieto J I，Nebot E M.Approximate inferencein state-space models with heavy tailed noise[J].IEEE Transac-tions on Signal Processing，2012，60（10）：5024-5037.

[2]Huang Y L，Zhang Y G，Li N，et al.Robust Students t basednonlinear filter and smoother[J].IEEE Transactions on Aero-space and Electronic Systems，2016，52（5）：2586-2596.

[3]Karlgaard C D，Schaub H.Huber-based divided difference fil-tering[J].Journal of Guidance，Control and Dynamics，2007，30（3）：885-891.

[4]Roth M，Ozkan E，Gustafsson F.A students t filter for heavytailed process and measurement noise[C]//（ICASSP），2013IEEE International Conference on IEEE，2013：5770-5774.

[5]Chen B D，Wang J J，Zhao H Q，et al.Convergence of a fixed-point algorithm under maximum correntropy criterion[J].IEEESignal Processing Letters，2015，22（10）：1723-1727.

[6]Huang Y L，Zhang Y G，Li N，et al.A novel robust studentst-based Kalman filter[J].IEEE Transactions on Aerospace andElectronic Systems，2017，53（3）：1545-1554.

[7]Zhang X Y，Wang G H，Song Z Y，et al.Hypersonic sliding targettracking in near space[J].Defence Technology，2015，11（4）：370-381.

[8]Fan Y，Lu F，Zhu W X，et al.A hybrid model algorithm forhypersonic glide vehicle maneuver tracking based on the aerody-namic model[J].Applied Sciences，2017，7（2）：159.

[9]Arasaratnam I，Haykin S.Cubature kalman filters[J].IEEETransactions on Automatic Control，2009，54（6）：1254-1269.

[10]Cui H，Liu W，Zhu S，et al.A hypersonic vehicle tracking algo-rithm based on the UKF generalized labeled multi-Bemoulli fil-ter[C]//Control Conference（CCC），2016 35th Chinese，D”？E，2016：4911-4916.

[11]Reuter S，VoBT，Vo B N，et al.The labeled multi Bernoullifilter[J].IEEE Transactions on Signal Processing，2014，62（12）：3246-3260.

[12]Julier S J，Uhlmann J K，Durrant-Whyte H F.A new approachfor filtering nonlinear systems[C] //American Control Confer-ence，Proceedings ofthe 1995，IEEE，1995（3）：1628-1632.

[13]Wang C C，Wu P L，Bo Y M，et al.Maneuvering target trackingwith adaptive debiased converted measurement filter[J].Journalof Chinese Inertial Technology，2017，25（1）：136-140.

[14]Deng Y H，WuPL，LiXX，et al.Decoupled unbiased con-verted measurements adaptive gating CPHD filter for airbornemulti-target tracking[J].Journal of Aeronautics，Astronauticsand Aviation，2017，49（1）：123-131.

[15]Yu M J.INS/GPS integration system using adaptive filter forestimating measurement noise variance[J].IEEE Transac-tions on Aerospace and Electronic Systems，2012，48（2）：1786-1792.

[16]Jiang B，Sheng W，Zhang R，et al.Range tracking methodbased on adaptive“current”statistical model with velocity pre-diction[J].Signal Processing，2017，131：261-270.

[17]Cao H J，Zhang X Y，Wu P L.Combining theIMM with secondorder debiased converted measurement for detonation control[J].International Journal of Computational Science and Engineering，2017，14（1）：84-91.

[18]Wu P L，LiXX，Zhang L Z，et al.Tracking algorithm withradar and infrared sensors using a novel adaptive grid interactingmultiple model[J].IET Science，Measurement&Technology，2014，8（5）：270-276.

[19]Zhu H，Leung H， He Z S.A variational Bayesian approach torobust sensor fusion based on student-t distribution[J].Informa-tion Sciences，2013，221：201-214.

[20]何廣軍，李檳檳，何其芳，等.基于IMM模型的周期性跳躍運動高超聲速飛行器跟蹤算法[J].彈箭與制導學報，2015，35（01）：37-40.