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

惲鵬 李星秀 吳盤龍 何山

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

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

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

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

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

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

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

1 基礎算法

1.1 IMM算法

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

1.1.1 狀態交互

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

1.1.2 狀態濾波

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

1.1.3 模型概率更新

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

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

1.1.4 狀態融合輸出

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

1.2 VB算法

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

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

2 IMM-VB算法

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

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

利用式（20）、式（21）可以獲取聯合概率密度函數

結合式（4）、式（5）、式（20）、式（21）可以計算出聯合概率密度函數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為狀態信息中的元素數目，利用式（17）、式（24）迭代獲取相應參數的近似分布：式中：C為保證qⁱ⁺¹（λ_j（k+1》積分為1的常量，后面的C均指代這樣的常量，但是對于每一個公式而言，其值可能不同。由式（27）可得：同理，可獲取狀態的近似分布：

由式（32）可推導出：

由高斯乘積定理可得：

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

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

3 仿真試驗

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

馬爾可夫轉移矩陣為：

量測噪聲滿足：

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

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

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

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

4 結束語

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

參考文獻

[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.