惲鵬 李星秀 吳盤龍 何山

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

關(guān)鍵詞：高超聲速飛行器;野值;學(xué)生七分布;重尾分布;IMM-VB

中圖分類號(hào)：TP24 文獻(xiàn)標(biāo)志碼：A

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

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

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

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

1 基礎(chǔ)算法

1.1 IMM算法

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

1.1.1 狀態(tài)交互

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

1.1.2 狀態(tài)濾波

狀態(tài)濾波過程與單模型濾波過程相同，依據(jù)卡爾曼濾波即可實(shí)現(xiàn)子模型下的狀態(tài)更新。將X_j⁰（k|k），P_j⁰（k|k）作為子模型的輸入進(jìn)行狀態(tài)預(yù)測(cè)與更新：式中：X_j（k+1|k）為狀態(tài)預(yù)測(cè)值P_j（k+1|k）為預(yù)測(cè)協(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)轉(zhuǎn)移矩陣，H_j（k+1）為量測(cè)矩陣Q_j（k）為過程噪聲協(xié)方差，R_j（k+1）為量測(cè)噪聲協(xié)方差。

1.1.3 模型概率更新

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

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

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

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

1.2 VB算法

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

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

2 IMM-VB算法

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

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

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

結(jié)合式（4）、式（5）、式（20）、式（21）可以計(jì)算出聯(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)信息中的元素?cái)?shù)目，利用式（17）、式（24）迭代獲取相應(yīng)參數(shù)的近似分布：式中：C為保證qⁱ⁺¹（λ_j（k+1》積分為1的常量，后面的C均指代這樣的常量，但是對(duì)于每一個(gè)公式而言，其值可能不同。由式（27）可得：同理，可獲取狀態(tài)的近似分布：

由式（32）可推導(dǎo)出：

由高斯乘積定理可得：

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

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

3 仿真試驗(yàn)

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

馬爾可夫轉(zhuǎn)移矩陣為：

量測(cè)噪聲滿足：

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

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

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

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

4 結(jié)束語(yǔ)

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

參考文獻(xiàn)

[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模型的周期性跳躍運(yùn)動(dòng)高超聲速飛行器跟蹤算法[J].彈箭與制導(dǎo)學(xué)報(bào)，2015，35（01）：37-40.