不確定重尾量測噪聲干擾下的魯棒目標跟蹤算法

馬天力， 張 揚， 劉 盼， 高 嵩

(西安工業大學電子信息工程學院，西安，710021)

目標跟蹤是在已知受噪聲干擾的量測信息條件下，利用相關濾波算法對目標狀態(位置，速度等)進行估計的過程，其在彈道導彈防御、空中目標監控和水下目標探測等領域應用廣泛[1-2]。然而系統的量測數據與目標動態參數間的關系大多為非線性，一般采用點估計和概率密度估計[3]的方法對此類系統狀態進行計算。點估計方法是通過非線性逼近策略輸出狀態均值和協方差矩陣獲得系統狀態。根據逼近策略可分為函數逼近、統計量逼近和隨機模型逼近；其代表算法分別為擴展卡爾曼濾波[4]、無跡卡爾曼濾波[5]和容積卡爾曼濾波(cubature kalman filter, CKF)[6]。概率密度估計則是通過構建狀態的概率密度表達，計算狀態后驗概率密度函數，并根據后驗信息獲得系統狀態[7]，該方法需要大量狀態對后驗概率密度函數進行描述，計算量較大，因此，實際系統大多采用點估計器進行狀態估計。

一般來說，跟蹤系統中量測噪聲的統計特性往往是未知且時變的。針對該類問題，文獻[8]提出改進Sage-Husa自適應Kalman濾波算法，運用協方差匹配判據來判斷濾波發散趨勢，并引入自適應衰減因子修正預測誤差協方差，從而調整濾波增益陣，抑制濾波發散。文獻[9]提出基于最大后驗概率的自適應無跡卡爾曼濾波算法，利用輸出的量測信息在線更新噪聲均值和協方差。Chang等[10]通過Huber代價函數重新表述量測噪聲，無需對非線性方程進行線性化處理，提高了濾波算法的估計精度。彭美康等[11]提出基于魯棒M估計的自適應容積卡爾曼濾波算法，通過修正觀測向量的誤差協方差，將損失函數最小化，在觀測噪聲污染率較高的情況下能夠自適應抑制野值對系統的干擾。Sarkka等[12]提出了變分貝葉斯自適應卡爾曼濾波算法，將變分貝葉斯方法與卡爾曼濾波結合，對模型狀態和時變量測噪聲參數進行聯合遞歸估計，一定程度上降低時變量測噪聲方差對系統的影響。雖然上述方法可以有效解決時變量測噪聲條件下非線性系統狀態估計問題，但均需要滿足量測噪聲服從高斯分布的假設，具有一定局限性[13]。

解決時變非高斯噪聲條件下的非線性濾波問題主要核心是對非高斯噪聲進行描述。Fan等[14]運用混合高斯模型對非高斯量測噪聲進行表示，利用集成卡爾曼濾波方法估計系統狀態；董祥祥等[15]運用Pearson Type VII分布對重尾噪聲進行建模，結合魯棒容積卡爾曼濾波器，實現狀態向量、模型參數的聯合估計；Zhu等[16]利用Student-T分布描述非高斯噪聲模型，采用變分學習算法對系統模型和噪聲參數進行聯合估計。在此基礎上，Huang等[17]將Student-T 分布加權積分分解為球面積分和徑向積分，并利用球面徑向變換規則進行積分數值化求解，提出基于Student-T分布的隨機容積濾波器。然而在實際目標跟蹤過程中，受傳感器設備老化、故障或電磁干擾等因素影響，導致量測噪聲出現不對稱重尾性，利用混合高斯模型，Student-T分布模型難以對量測噪聲進行準確表述，不準確的模型描述將嚴重降低非線性目標跟蹤系統的濾波性能。Wang等[18]用Skew-T模型對非對稱量測噪聲進行建模，提出基于Skew-T分布擴展卡爾曼濾波算法(Skew-T extended kalman filter, SkewEKF)，假設噪聲參數已知，一旦先驗噪聲信息與實際值存在較大偏差，難以保證準確的濾波效果。

因此，本文針對不確定重尾量測噪聲條件下非線性系統狀態估計問題，提出基于變分推理的魯棒容積卡爾曼濾波(variational inference robust cubature kalman filter, VIRCKF)算法，該方法運用Skew-T分布對不對稱重尾噪聲進行表示，然后基于三階球面徑向規則計算狀態和量測的容積點，并通過采樣所得的容積點逼近系統狀態和協方差，同時結合變分推理學習估計得到近似后驗分布，并對系統狀態進行更新，從而獲得準確的目標狀態。

1 問題描述

考慮不對稱重尾量測噪聲下的目標跟蹤問題, 建立系統模型：

xk=Axk-1+wk-1

(1)

zk=h(xk)+ek

(2)

式中：xk為k時刻系統狀態向量；A為狀態轉移矩陣；wk-1為過程噪聲，假設其服從零均值、協方差為Qk-1的高斯分布；h(·)為量測方程；zk是k時刻的傳感器接收到的量測值；ek為不對稱重尾量測噪聲，且Cov(wk,ek)=0。

(3)

式中：ST(·)表示Skew-T分布；[ek]i表示k時刻噪聲的第i個元素；μ、σ和δ分別是位置參數，擴展參數和模型偏斜系數，v表示自由度，σ、δ、v>0。

以一維Skew-T分布為例，當i=1時，式(3)可表示為[19]：

(4)

且

(5)

(6)

式中：t(·;μ,σ2,v)為Student-T分布的概率密度函數；T(·;0,1,v)是Student-T分布的累積分布函數；Γ(·)表示伽馬函數。

Skew-T分布、Student-T分布和高斯分布的概率密度函數如圖1所示。當δ趨于0時，Skew-T分布接近Student-T分布；當δ趨于0且v趨于∞時，Skew-T分布接近正態分布。由圖1可以看出，Skew-T分布偏斜更明顯，更適合對不對稱重尾噪聲進行模擬。

圖1 3種分布的概率密度函數

2 基于變分推理的魯棒容積卡爾曼濾波器

2.1 變分推斷

假設系統狀態向量先驗服從高斯分布，則式(1)表示為p(xk|xk|k-1)=N(xk;Axk|k-1,Qk-1)，系統一步預測的概率密度函數p(xk|z1:k-1)為：

p(xk|z1:k-1)=N(xk;xk|k-1,Pk|k-1)

(7)

結合式(2)和(3)得到量測似然函數：

p(zk|xk)∶ST(zk;h(xk),R,Δ,v)

(8)

由于式(8)中似然函數為Skew-T分布，其不滿足貝葉斯推斷中參數屬于共軛指數分布簇的假設條件[20]，造成參數難以進行求解。因此，將式(8)的似然函數表示為：

(9)

(10)

(11)

式中：NH(·)為半正態分布，G(·)為伽馬分布。

根據貝葉斯規則[21]，參數集合的后驗概率密度函數計算為：

(12)

式(12)的難點在于邊緣似然函數p(z1:k)的計算，故引入變分貝葉斯理論，考慮尋找一個簡單近似分布q(xk,uk,Λk)去逼近真實后驗分布p(xk,uk,Λk|z1:k)，取p(z1:k)的對數形式：

KL(q(xk,uk,Λk)‖p(xk,uk,Λk|z1:k))

(13)

KL(q(xk,uk,Λk)‖p(xk,uk,Λk|z1:k)) 為KL散度，其用來衡量近似分布q(xk,uk,Λk)與真實后驗分布p(xk,uk,Λk|z1:k)之間的差異[22]，由于KL散度非負，只需令變分下界F(q(xk,uk,Λk))達到最大，即可獲得后驗概率最優近似分布。

系統聯合概率密度函數為：

(14)

根據均值域逼近理論[23]，后驗概率密度函數的近似因子分解形式為：

p(xk,uk,Λk|z1:k)≈qx(xk)qu(uk)qΛ(Λk)

(15)

2.2 變分貝葉斯容積卡爾曼濾波

非線性濾波問題可以歸結為一個積分計算問題，被積函數為非線性函數和密度乘積的形式。容積卡爾曼濾波算法是將積分轉換為球面徑向形式，采用球面徑向規則進行數值積分。由于包含Skew-T分布的后驗概率密度函數形式復雜，考慮將變分推斷引入CKF算法。在狀態預測階段，因狀態向量服從高斯分布，可用三階球面徑向規則計算標準的高斯加權積分，即通過一組容積點數值化計算系統的狀態均值和協方差。

(16)

?i,k-1|k-1=Sk-1|k-1ξi+xk-1|k-1

(17)

依據式(1)和式(7)，狀態向量的先驗均值xk|k-1和協方差矩陣Pk|k-1可表示為：

(18)

(19)

對容積點?i,k-1|k-1進行更新，

?i,k|k-1=A?i,k-1|k-1

(20)

則狀態向量先驗均值和協方差矩陣可表示為：

(21)

(22)

在量測更新階段，量測服從Skew-T分布，通過變分推斷計算其近似形式，用更新的容積點計算系統量測均值和協方差。其量測容積點計算如下：

(23)

?i,k|k-1=Sk|k-1ξi+xk|k-1

(24)

結合式(7)和(8)，量測向量的先驗概率密度p(zk|z1:k-1)為：

(25)

量測向量先驗均值zk|k-1和協方差矩陣Pzz,k|k-1表示為：

(26)

(27)

狀態向量和量測的互協方差矩陣Pxz,k|k-1為:

(28)

其中：

p(xk,zk|z1:k-1)=p(zk|xk)p(xk|z1:k-1)=

ST(zk;h(xk),R,Δ,v)N(xk;xk|k-1,Pk|k-1)

(29)

考慮Skew-T分布用其近似分布式(9)代替，式(27)～(28)可以表示為式(32)～(33) ：

zi,k|k-1=h(?i,k|k-1)

(30)

(31)

(32)

(33)

系統狀態的后驗概率密度函數為：

p(xk,uk,Λk|z1:k)=

p(xk|z1:k-1)p(Λk)·p(uk|Λk)p(zk|xk,uk,Λk)

(34)

由于后驗概率密度函數形式復雜，根據式(15)，通過求解logqx(xk)、logqu(uk)和logqΛ(Λk)可得到系統的狀態和協方差，具體計算過程如下：

logqx(xk)=

(35)

由于變分貝葉斯近似適用于所有共軛指數域中的分布函數，那么后驗可由其先驗唯一確定[24]，所以其后驗概率密度函數的形式為：

qx(xk)=N(xk;xk|k,Pk|k)

(36)

式中：更新后的卡爾曼濾波增益Kk、系統狀態xk|k和協方差矩陣Pk|k分別為：

(37)

xk|k=xk|k-1+Kk(zk-k|k-1-Δūk)

(38)

(39)

同理，對于隱變量uk，有:

logqu(uk)=

(40)

其后驗概率密度函數的形式為：

qu(uk)=NT(uk;uk|k,Uk|k)

(41)

式中：NT(·)表示截斷正態分布。qu(uk)的后驗概率密度參數的更新方程為：

uk|k=Ku(zk-h(xk|k))

(42)

(43)

(44)

對于似然函數參數Λk，有：

logqΛ(Λk)

=

(45)

其后驗概率密度函數的形式為：

(46)

其中：

φk=R-1((zk-h(xk|k))(zk-h(xk|k))T+Pzz,k|k)+

(47)

式中：式(35)、(40)、(45)中的cx、cu和cΛ是與變量xk、uk和Λk相關的常數。

3 仿真實驗

為了驗證所提算法的有效性，考慮時變Skew-T分布噪聲模型，分別采用CKF、SkewEKF以及VIRCKF算法對目標跟蹤系統進行仿真實驗。本文主要通過3種算法的均方根誤差(RMSE)來說明算法的濾波性能。令目標運動模型為勻速模型，目標運動方程為：

xk=Axk-1+wk-1

(48)

(49)

式中：采樣時間T=1 s，過程噪聲協方差矩陣Q=diag([5×10-2,5×10-3,5×10-2,5×10-3])。

量測方程為：

(50)

式中：zk=[rkθk]T，rk和θk分別表示k時刻目標徑向距離與方位角。量測噪聲參數設置如下：

假設量測噪聲協方差初始值R0=I2，令目標的初始位置為x0|0=[500 m,25 m/s,500 m,-20 m/s]T，狀態協方差矩陣初始值P0=diag([2,10-2,2,10-2])。仿真時間為200 s。對3種算法分別進行100次Monte Carlo實驗。

圖2、圖3分別為3種算法在X、Y方向的位置和速度均方根誤差。從圖2可以看出，當t>100 s時，由于量測噪聲的突然增大，CKF算法的均方根誤差出現波動；當t>150 s時，量測噪聲再次增大，CKF算法的RMSE也隨之增大。SkewEKF算法的RMSE明顯小于CKF，因為SkewEKF算法的本質是利用泰勒級數展開將非線性系統進行線性化，存在一定誤差。而本文所提算法通過構建近似分布，迭代計算，直至逼近系統真實后驗概率密度函數，具有自適應性。從圖3可以看出，在150 s噪聲增大時，3種算法速度RMSE均增大，但VIRCKF算法RMSE更小。

圖2 不同方向的位置均方根誤差

圖3 不同方向的速度均方根誤差

通過計算VIRCKF算法的變分下界函數，對算法的收斂性進行分析。圖4為實驗過程中系統狀態的變分下界函數。可以看出，該算法在X方向位置變分下界穩定在[-3.40,-3.00]，在Y方向最終穩定在[-4.41,-3.65]，表明算法穩定收斂。

圖4 變分下界函數

圖5為量測噪聲的近似分布擬合真實分布的過程，可以看出，隨著參數的迭代優化，近似分布會不斷逼近真實分布。

圖5 量測噪聲的真實分布與近似分布的擬合過程

圖6 變分下界函數隨值的變化

圖7 變分下界函數隨值的變化

分別對3種算法的計算復雜度進行分析，CKF算法的時間復雜度為O(6Nn)，SkewEKF的時間復雜度為O(Nm)，VIRCKF算法的時間復雜度為O(8Nm+6Nnm+10N)。其中，N為總運行時間，2n表示容積點數，m為變分迭代次數。雖然VIRCKF的時間復雜度更高，但綜合3種算法的濾波效果和魯棒性分析，本文所提的VIRCKF算法濾波精度更高，性能更加優越。

4 結語

本文針對不確定重尾量測噪聲條件下非線性系統目標跟蹤問題，提出基于變分推理的魯棒容積卡爾曼濾波算法。將Skew-T分布的概率密度函數用近似形式代替，然后利用球面徑向規則對函數進行數值積分求解，同時結合變分推理學習估計得到系統的近似聯合后驗分布。實驗結果表明，在量測噪聲為重尾分布且參數不確定時，本文算法具有較高的估計精度和較好的濾波效果。未來可以探索更有效的估計算法，解決不確定過程噪聲與量測噪聲下的機動目標跟蹤問題。