張 揚,馬天力,高 嵩,陳超波

(西安工業(yè)大學電子信息工程學院,陜西 西安 710021)

1 引言

近年來,隨著航空、航天、航海事業(yè)的不斷發(fā)展,目標跟蹤技術(shù)在對海底、海面、陸地和空中等領(lǐng)域取得了極大的應(yīng)用。在目標跟蹤過程中,傳感器測量信息不僅包含目標信號,同時也包含隨機量測噪聲。目標跟蹤要解決的基本問題就是從被噪聲污染的量測信息中盡可能充分地濾除干擾噪聲,使系統(tǒng)的狀態(tài)估計更加準確。傳統(tǒng)的濾波方法多采用卡爾曼濾波器[1],其基于最小均方誤差準則,對于準確建模且噪聲服從高斯分布的動態(tài)系統(tǒng),KF能夠遞歸得到狀態(tài)的一致最小方差線性無偏估計,但要求量測噪聲滿足高斯分布的假設(shè)。通常,量測噪聲的統(tǒng)計特性可以事先由傳感器的物理特性得到。但對于跟蹤系統(tǒng)而言,由于外界干擾、加速度的物理特性和環(huán)境等因素影響,很難用一個準確的統(tǒng)計特性來表示量測噪聲,所以其統(tǒng)計特性往往是未知且時變的。針對量測噪聲時變的目標跟蹤問題,一般采用自適應(yīng)卡爾曼濾波器(AKF),例如Sage-Husa自適應(yīng)濾波器(SHAKF)[2-4],基于新息的自適應(yīng)濾波器(IAKF)[5-8]和交互式多模型自適應(yīng)濾波器(IMMAKF)[9,10]。SHAKF可以用最大后驗準則估計量測噪聲的統(tǒng)計量,但是它不能保證量測噪聲協(xié)方差的收斂性,易造成濾波發(fā)散[11]。IAKF可以通過新息協(xié)方差對量測噪聲進行估計和調(diào)整,使濾波算法能更好地適應(yīng)變化的噪聲統(tǒng)計特性[12]。但是,其只適用于緩慢變化的量測噪聲協(xié)方差,因為實現(xiàn)量測噪聲協(xié)方差的可靠估計需要很大的數(shù)據(jù)窗口[13]。IMMAKF可以使用多個模型來描述目標的運動狀態(tài),通過有效的加權(quán)融合進行目標運動狀態(tài)估計,很好地克服了單模型濾波算法估計誤差較大的問題,但其計算比較復雜[14]。

在實際的工程應(yīng)用中,跟蹤目標容易受到來自外界的干擾,量測噪聲存在野值,導致噪聲分布出現(xiàn)重尾特征。例如在雷達跟蹤系統(tǒng)中,由于目標位置的隨機擺動,使得對目標位置的量測伴隨著閃爍噪聲的出現(xiàn),這種閃爍噪聲導致系統(tǒng)的量測呈現(xiàn)非高斯特性[15]。直接使用上述方法會導致算法性能降低,濾波精度下降。為解決重尾量測噪聲的濾波問題,Guo等[16]提出了一種混合高斯估計算法(GM-CDMKF)解決非高斯噪聲的目標跟蹤問題,在不增加過多計算的情況下,改善了非高斯噪聲的濾波性能。文獻[17]研究了用混合拉普拉斯模型處理重尾噪聲問題,該方法使用梯度下降法求解負似然對數(shù)最小值來估計參數(shù)。Roth等[18]將量測噪聲的統(tǒng)計特性近似為T分布,推導出了一種魯棒T分布濾波算法,該方法借助T分布對量測異常值建模,最終通過T分布來近似系統(tǒng)的后驗概率密度函數(shù)。Zhu等[19]提出基于T分布量測噪聲的VBT濾波方法,利用變分貝葉斯學習算法進行模型的參數(shù)估計,可以解決重尾噪聲模型的濾波問題,提高了估計精度。但是基于T分布的濾波算法性能可能會因自由度參數(shù)的選擇不當而降低。文獻[20]將粒子濾波算法用于T分布的非高斯噪聲模型,該算法基于蒙特卡洛思想,將粒子進行重采樣以獲得足夠的統(tǒng)計數(shù)據(jù),但其在改善跟蹤性能的同時也需要大量的粒子,使計算量大大增加[21]。在目標跟蹤的實際應(yīng)用中,為了保持最優(yōu)濾波系統(tǒng)的有效性,需要比較精確的表征傳感數(shù)據(jù)的統(tǒng)計特性。而可靠性較差的傳感器誘導出的異常量測值,傳感器發(fā)生故障或受到環(huán)境中的電磁干擾,量測模型往往伴隨著不可避免的重尾偏斜噪聲,這會嚴重降低以高斯分布或T分布為前提假設(shè)的目標跟蹤系統(tǒng)的濾波性能。

針對量測噪聲重尾不對稱條件下的線性目標跟蹤系統(tǒng)的狀態(tài)估計問題,本文提出了一種解決重尾不對稱量測噪聲的VBST算法。每次時間更新后,用變分貝葉斯學習估計得到近似后驗分布,再根據(jù)量測噪聲對目標的速度和位置狀態(tài)進行更新,從而達到自適應(yīng)濾波的目的。通過與傳統(tǒng)的卡爾曼濾波算法、VBT算法以及粒子濾波算法對比,驗證本文所提算法的濾波性能和魯棒性。

2 問題描述

建立如(1)式的偏斜T分布量測噪聲的系統(tǒng)模型

xk=Axk|k-1+wk

zk=Cxk+ek

(1)

其中,xk∈n是狀態(tài)向量。A∈n×n是狀態(tài)轉(zhuǎn)移矩陣,wk∈n是過程噪聲,假設(shè)其服從均值為零,方差為Q的高斯白噪聲。zk∈m是在k時刻的量測值,C∈m×n是量測模型矩陣。ek∈m是獨立于wk的重尾不對稱的非高斯量測噪聲

(2)

當i=1時,式(2)為一維偏斜T分布,其概率密度函數(shù)可以表示為[22]

(3)

其中

(4)

(5)

T(·;0,1,v)是T分布的累積分布函數(shù)。Γ(·)表示伽馬函數(shù)。圖1給出了不同偏斜系數(shù)下偏斜T分布的概率密度函數(shù)圖[23],其中,偏斜T分布的其它參數(shù)設(shè)置為:μ=0,σ=1,v=4。

圖1 不同偏斜系數(shù)下的偏斜T分布概率密度函數(shù)

如圖1所示,若偏斜系數(shù)大于0,偏斜T的概率密度函數(shù)左偏;若偏斜系數(shù)小于0,偏斜T分布的概率密度函數(shù)右偏。

當δ趨于0時,偏斜T分布近似為學生T分布;當δ趨于0且v趨于∞時,偏斜T分布近似為正態(tài)分布[20]。

根據(jù)狀態(tài)空間模型,p(xk|z1:k-1)的概率密度函數(shù)為

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

(6)

z1,…zN是來自ST(zk;Cxk,R,Δ,v)的N個獨立的觀測值,其似然函數(shù)為[24]:

(7)

(8)

(9)

其中,N(·)表示高斯分布,HN(·)表示半正態(tài)分布,G(·)表示伽馬分布。

基于貝葉斯規(guī)則,其聯(lián)合概率分布為

(10)

圖2表示偏斜T分布量測噪聲系統(tǒng)的圖模型,圖中方框表示先驗的超參數(shù),深色的圓形表示系統(tǒng)的量測。

圖2 偏斜T分布的圖模型

3 基于偏斜T噪聲的魯棒狀態(tài)估計

3.1 變分貝葉斯原理分析

基于貝葉斯定理

(11)

其中,z表示觀測數(shù)據(jù)的集合,模型的未知參數(shù)為φ={X,u,Λ},隱變量為X={x1,x2,…,xN},先驗的超參數(shù)ψ={Δ,v}。則

p(z|φ,ψ)=p(z|φ)

(12)

根據(jù)式(12)可以將式(11)寫為

(13)

邊緣似然函數(shù)的對數(shù)logp(z|φ)表示為

logp(z|φ)=F(q(φ))+KL(q(φ)‖p(φ|z,ψ))

(14)

其中,F(q(φ))為自由能量函數(shù)[25-26]。

(15)

KL(q(φ)‖p(φ|z,ψ))是近似后驗q(φ)與真實后驗p(φ|z,ψ)之間的Kullback-Liebler散度[27]。

(16)

變分的思想是尋找一個簡單易求的概率分布q(φ)去逼近真實后驗分布p(φ|z,ψ),即讓KL散度值達到最小。由于KL散度是非負的,當且僅當q(φ)=p(φ|z,ψ)時取值為零,最小化KL散度相當于最大化自由能量函數(shù),因此可以通過使自由能量達到最大來實現(xiàn)目標。

根據(jù)均值域逼近理論[28],q(φ)可以分解為

(17)

在VB步驟中,分解的真實后驗分布與近似后驗分布的最小KL散度見下式

(18)

其中,Q=qx(xk)qu(uk)qΛ(Λk)。

根據(jù)式(17),近似后驗分布表示為

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

(19)

則z,φ聯(lián)合分布的對數(shù)表示如下式所示

log(p(z|φ)p(φ|ψ))

(20)

3.2 遞歸VB卡爾曼濾波的偏斜T算法

通過求解logqx(xk),logqu(uk)和logqΛ(Λk)可以得到其估計的后驗概率密度函數(shù)。

對于目標狀態(tài)

logqx(xk)

(21)

其后驗概率密度函數(shù)的形式如式(22)所示

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

(22)

其中,狀態(tài)xk|k,協(xié)方差矩陣Pk|k和卡爾曼濾波增益Kx的更新方程如下式所示。

(23)

Pk|k=(I-KkC)Pk|k-1

(24)

(25)

同理,對于模型的未知參數(shù)uk

logqu(uk)

(26)

其后驗概率密度函數(shù)的形式如式(27)所示。

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

(27)

其中,TN(·)表示截斷正態(tài)分布。qu(uk)的后驗概率密度參數(shù)的更新方程[29]如下式所示

uk|k=Ku(zk-Cxk|k)

(28)

(29)

=Δ(Δ2+R)-1

(30)

對于模型的未知參數(shù)Λk

logqΛ(Λk)

(31)

其后驗概率密度函數(shù)的形式如下式所示

(32)

其中

φk=R-1((zk-Cxk|k)(zk-Cxk|k)T+CPk|kCT)

(33)

其中,式(21)(26)(31)中的cx、cu和cΛ是和變量xk、uk和Λk相關(guān)的常數(shù)。

4 仿真研究

為了驗證所提算法的有效性,分別采用卡爾曼濾波算法(KF),VBT算法,粒子濾波算法(PF)以及VBST算法對目標跟蹤系統(tǒng)進行仿真。

假設(shè)目標的初始位置在x方向500m處,速度為1.5m/s,在y方向500m處,速度為1m/s。觀測站的位置在原點(0,0),粒子數(shù)設(shè)為1000。協(xié)方差矩陣P0=diag([2,10-2,2,10-2]),采樣時間為0.1s,仿真時間200s,過程噪聲為零均值的高斯白噪聲, 其協(xié)方差矩陣Q=diag([10-2,10-3,10-2,10-3])。

4.1 實驗一:非時變偏斜T噪聲模型

假設(shè)量測噪聲的偏斜系數(shù)為0.9,自由度為10,R=2×I2×2。本文主要通過四種算法的均方根誤差(RMSE)來說明算法的濾波性能。

圖3為真實軌跡和KF算法、VBT算法、粒子濾波算法以及VBST算法的濾波軌跡。從圖中可以直觀看出VBST算法的濾波性能要明顯優(yōu)于另外三種算法,粒子濾波算法次之,KF和VBT的算法相對而言較差。

圖3 真實軌跡及四種算法濾波后軌跡

圖4和圖5是四種算法的位置均方根誤差和速度均方根誤差,從圖中可以看出,VBST算法具有更小的位置和速度均方根誤差,其濾波性能更好。主要原因是VBST算法不斷的更新隱變量和參數(shù),直到算法收斂,使得到的近似后驗分布逼近真實后驗分布。雖然VBT算法也用到了變分貝葉斯方法,但是T分布不能很好地擬合重尾非對稱噪聲,所以在本系統(tǒng)的量測噪聲模型中,其估計的后驗分布與真實后驗分布相差太大,造成濾波性能下降。

圖4 四種算法的位置均方根誤差

圖5 四種算法的速度均方根誤差

4.2 實驗二:時變偏斜T噪聲模型

為了驗證所提算法的魯棒性,考慮時變的偏斜T分布噪聲模型,實驗參數(shù)設(shè)置如下:t<150s設(shè)置偏斜系數(shù)為0.9,自由度為10,R=2×I2×2。t∈[150,200]s改變噪聲分布,增大噪聲協(xié)方差,設(shè)置偏斜系數(shù)為-0.6,自由度為40,R=5×I2×2,其仿真結(jié)果如圖6、7和8所示。

圖6 時變噪聲下的四種算法濾波軌跡

從圖6可以看出,改變噪聲分布后,系統(tǒng)的量測噪聲在x方向525m處突然增大,KF和VBT算法的濾波效果明顯變差,粒子濾波算法也有輕微的波動,但是VBST的濾波效果并未受到很大影響。

圖7和圖8為四種算法在不同量測噪聲下的位置和速度均方根誤差,從圖中可以看出,在150s改變噪聲分布后,KF和VBT算法的均方根誤差明顯增大,粒子濾波算法相比較KF和VBT算法而言,濾波效果較好,位置和速度的RMSE值也小,而VBST算法的位置均方根誤差和速度均方根誤差都比較小,說明在時變噪聲條件下本文所提算法和粒子濾波算法較于其它兩種算法而言魯棒性更強。

圖7 時變噪聲下的四種算法位置均方根誤差

圖8 時變噪聲下的四種算法速度均方根誤差

以實驗一的仿真模型為例,主要通過求解變分下界函數(shù)L,也表示自由能量函數(shù)F,分析VBST算法的收斂性。設(shè)置迭代次數(shù)為30,進行100次蒙特卡羅實驗,仿真結(jié)果如圖9所示,由下圖可以看出,L最終穩(wěn)定在-0.63±0.03,說明算法收斂。

圖9 變分自由能量變化過程

表1為單次蒙特卡羅實驗四種算法的相對運行時間對比,假設(shè)卡爾曼濾波算法運算時間為單位1。

表1 不同算法的運行時間

從表1可以看出,粒子濾波算法的計算量較大,用時較長。KF和VBT算法的運算時間短,但濾波效果差。VBST算法濾波性能較好,運行時間較粒子濾波而言小的多。從圖8圖9可以看出粒子濾波和VBST算法的濾波性能相似,但從算法的復雜度來看,粒子濾波的時間復雜度為O(n+Nn+2Nn2),VBST算法的時間復雜度為O(Mn),其中,n=200,N=1000,M=30,所以粒子濾波的計算要比VBST算法復雜得多。綜上分析,本文所提的VBST算法濾波精度更高,性能更加優(yōu)越。

5 結(jié)論

本文考慮量測噪聲為偏斜T分布時的目標跟蹤問題,提出了一種基于變分貝葉斯理論的VBST濾波算法。實驗結(jié)果表明,在量測噪聲為偏斜T分布時,本文所提算法具有較高的估計精度和較好的濾波效果。在時變噪聲條件下,VBST算法也有很好的魯棒性。鑒于VBST算法的優(yōu)越性能,未來考慮非線性系統(tǒng)模型和混合不確定量測噪聲模型。