存在非高斯重尾分布噪聲的純方位目標跟蹤算法

劉燦, 王輝, 林德福, 崔曉曦, 徐晗暉

(1.北京理工大學 宇航學院, 北京 100081; 2.中國兵器工業(yè)導(dǎo)航與控制技術(shù)研究所, 北京 100089;3.北京理工大學 設(shè)計與藝術(shù)學院, 北京 100081)

0 引言

純方位目標跟蹤算法是目標跟蹤研究中的熱點,其主要應(yīng)用于無源目標定位系統(tǒng)。相對于有源目標定位系統(tǒng),無源目標定位系統(tǒng)本身不向外發(fā)射電磁信號,而是獲取目標的輻射源來進行跟蹤定位,具有隱蔽性高、生存能力強的特點[1],在雷達被動定位、航海聲吶探測、無線傳感器網(wǎng)絡(luò)等軍民領(lǐng)域中均有著重要的應(yīng)用。

純方位目標跟蹤問題的跟蹤誤差主要來自于模型的非線性和測量噪聲的非高斯性[2],針對模型的非線性的主要挑戰(zhàn)是純方位角度測量和目標動力學之間的非線性關(guān)系。目前已有多種算法應(yīng)用于純方位目標跟蹤中,早期的研究工作使用擴展的卡爾曼濾波器(EKF)[3],EKF采用1階泰勒展開近似非線性量測方程和目標動力學方程,得到近似的線性測量方程和目標動力學方程。文獻[4-5]通過在修正的極坐標中建立純方位角目標跟蹤問題,提高了EKF的穩(wěn)定性,從而得到了修正的極坐標下的EKF(MPEKF)。除了測量方程與目標動力學的非線性影響,初始狀態(tài)值的誤差對于EKF的收斂也是至關(guān)重要的,較大的初始狀態(tài)誤差會引起發(fā)散問題[6]。此外還有更復(fù)雜的卡爾曼濾波算法應(yīng)用于純方位目標跟蹤研究,如無跡卡爾曼濾波器(UKF)[7],其是一種基于確定性樣點計算的濾波,不依賴于近似非線性狀態(tài)和量測模型。Arasaratnam等[8]提出容積卡爾曼濾波(CKF),其基于3階球面徑向容積準則,并使用一組容積點來逼近具有附加高斯噪聲的非線性系統(tǒng)的狀態(tài)均值和協(xié)方差。文獻[9]針對混合噪聲提出一種應(yīng)用于被動定位算法的自適應(yīng)容積卡爾曼濾波(ACKF)算法。而粒子濾波器(PF)是一種基于蒙特卡洛積分的最優(yōu)非線性濾波器[10],文獻[11-12]將其應(yīng)用于純方位目標跟蹤,但粒子濾波器的一個主要缺點是計算復(fù)雜度高,計算時間長,使得其在純方位目標跟蹤應(yīng)用中的實用性有限。MILLER等[13]針對純方位目標跟蹤,用構(gòu)成線性的偽線性方程代替非線性角度測量來得到線性化的遞歸貝葉斯估計的偽線性估計器(PLKF)。Nguyen[14]進一步發(fā)展了一種瞬時誤差偏置補償卡爾曼濾波(BC-PLEKF),通過偏置補償測量方程的非線性,以此提高濾波器的性能。

而在非高斯噪聲條件下,影響目標的跟蹤精度的主要因素是測量噪聲的非高斯性而不是模型的非線性[2]。在一些工程應(yīng)用中,若使用不可靠傳感器的異常測量值跟蹤敏捷目標,系統(tǒng)會存在非高斯的重尾分布過程噪聲和測量噪聲[15-16]。在實際的應(yīng)用場景中,地雜波、海雜波和偶發(fā)電磁干擾,使得雷達接收到的回波幅度突變,遠高于正常值,在參數(shù)估計中形成異常值[17]。此外,目標朝向雷達的變化可能導(dǎo)致不規(guī)則的電磁波反射亦會在角度跟蹤中產(chǎn)生異常值,此時測量噪聲不服從高斯分布,往往呈現(xiàn)為重尾分布[18]。而目標飛行器的強機動以及受到的氣流影響,會導(dǎo)致跟蹤模型的狀態(tài)值發(fā)生突變,誘導(dǎo)非高斯重尾分布系統(tǒng)噪聲[15]。傳統(tǒng)的卡爾曼濾波器的性能在這些具有非高斯重尾分布的系統(tǒng)過程和測量噪聲的工程應(yīng)用中會下降甚至失效。近年來,國內(nèi)外提出了大量的魯棒濾波器來解決具有重尾測量噪聲的線性系統(tǒng)的濾波問題,如學生t分布混合濾波器[19]、異常魯棒卡爾曼濾波器[20],以及基于學生t分布和變分貝葉斯(VB)的魯棒卡爾曼濾波器[21-23]。然而在非高斯重尾系統(tǒng)過程噪聲的情況下,這些魯棒濾波器可能會失效,因為它們都假設(shè)系統(tǒng)過程噪聲為高斯分布[19-21]。

近年來,采用變分方法來解決非高斯噪聲下的非線性估計問題成為國際研究的熱點。VB推斷是一類近似計算復(fù)雜積分的方法,主要用于復(fù)雜的統(tǒng)計模型中,可以用于近似不可觀測變量的后驗概率估計,對于給定的模型給出觀測變量的邊緣似然函數(shù)的下界,然后通過相互似然的等式進行不斷的迭代來獲得最優(yōu)解。文獻[23]提出的PEKF-VB(Proximal Extended Kalman Filter-Variational Bayes)濾波器采用KL散度作為后驗概率密度函數(shù)(PDF)的度量并推導(dǎo)出了另一種基于VB技術(shù)的非線性估計算法。文獻[24]將PEKF-VB應(yīng)用于星載雷達非線性目標跟蹤,實現(xiàn)了星載雷達的跟蹤精度的提高。文獻[25]提出一種基于學生t分布的VB卡爾曼濾波器,將一步預(yù)測的PDF和似然PDF近似為不同自由度參數(shù)的學生t分布,用VB推斷來近似重尾分布的系統(tǒng)過程噪聲和測量噪聲,提高了濾波器的估計精度。

本文基于文獻[25]提出一種新的針對非高斯重尾分布噪聲的魯棒性濾波,VB修正增益卡爾曼濾波(VBMGEKF)。本文通過構(gòu)造基于重尾分布的系統(tǒng)過程噪聲和測量噪聲,將一步預(yù)測的PDF和似然PDF分別建模為不同參數(shù)的帶異常值的高斯分布來近似重尾分布,并將未知狀態(tài)協(xié)方差矩陣的PDF模型建模為逆Wishart分布。通過引入輔助隨機變量,提出一個層次高斯狀態(tài)空間模型。在此基礎(chǔ)上,通過VB推斷得到近似后驗PDF,以及未知狀態(tài)協(xié)方差矩陣和輔助隨機變量,推導(dǎo)出系統(tǒng)具有重尾分布的過程和測量噪聲的線性狀態(tài)空間模型后驗PDF的近似閉環(huán)解析解,并結(jié)合修正增益卡爾曼濾波(MGEKF),降低純方位目標跟蹤問題中量測方程非線性的影響。仿真結(jié)果表明,該濾波器比現(xiàn)有的一些濾波算法具有較小的均方根誤差(RMSE)。

1 問題描述

考慮一個線性離散時間狀態(tài)空間模型描述的運動學系統(tǒng)為

xk=fk(xk-1)+wk

(1)

zk=hk(xk)+vk

(2)

式中:xk為目標k時刻的狀態(tài)值;zk為k時刻的測量值;fk(xk-1)為與xk-1相關(guān)的狀態(tài)預(yù)測方程;hk(xk)為系統(tǒng)狀態(tài)與測量值相關(guān)的方程;wk和vk為系統(tǒng)的過程噪聲和測量噪聲,一般假設(shè)wk和vk均為獨立的高斯分布,本文認為wk和vk為帶異常值的重尾噪聲。因此考慮將測量似然和一步預(yù)測建模為如下帶異常值的高斯分布,即

(3)

(4)

(5)

(6)

(7)

Fk-1為狀態(tài)轉(zhuǎn)移矩陣,Pk-1|k-1為狀態(tài)協(xié)方差,Qk為系統(tǒng)的過程噪聲協(xié)方差。當Qk為高斯分布時Σk=Pk|k-1,當Qk為重尾分布時Σk就不能等于Pk|k-1。本文使用VB來近似推斷未知的矩陣Σk。對于未知矩陣Σk,需要首先選擇一個共軛先驗分布,因為該共軛可以保證后驗分布具有相同的函數(shù)形式作為先驗分布。在貝葉斯統(tǒng)計中,逆Wishart分布通常用作具有已知均值的高斯分布的協(xié)方差矩陣的共軛先驗[25-26]。

(8)

因此定義Σk為逆Wishart分布,則未知矩陣的先驗分布為

Uk=(uk-n-1)Pk|k-1

(11)

式中:n為xk維數(shù)。

同時在此處使用一種啟發(fā)式算法[27],引入一個遺忘因子τ,減少前一時刻狀態(tài)協(xié)方差矩陣的影響,一般取τ∈(0,1),則

將式(12)代入式(11),可得

uk=n+τ+1

(13)

至此,一個新的基于高斯分布的狀態(tài)空間模型可由式(5)、式(6)、式(11)～式(13)組成。

(14)

(15)

圖1所示為層次高斯SSM。

圖1 層次高斯SSMFig.1 Hierarchical Gaussian State Space Model

則由式(5)～式(10)以及式(12)、式(13)構(gòu)成一個基于高斯分布的層次高斯狀態(tài)空間模型,從而將帶有重尾分布的系統(tǒng)過程和測量噪聲的線性狀態(tài)空間模型的狀態(tài)估計問題轉(zhuǎn)化為基于高斯分布的層次高斯SSM的狀態(tài)估計問題。

2 變分貝葉斯推斷

根據(jù)貝葉斯定理,本文可以得到由式(5)～式(13)構(gòu)成的層次高斯狀態(tài)空間模型的聯(lián)合后驗分布p(Φk|z1:k,Ψk),即

(16)

式中:Φk={xk,Rk,Σk,λk,ξk},Ψk={uk,Uk,ak,bk,ek,fk}分別為模型的未知變量和先驗值的參數(shù)。由于聯(lián)合后驗分布PDFp(Φk|z1:k,Ψk)難以求得解析解,本文采取VB方法來獲取聯(lián)合后驗分布的次優(yōu)近似解[21]。VB使用一個分布q(Φk)來近似后驗分布,對數(shù)似然邊際logp(z1:k|Φk)可表示為

logp(z1:k|Φk)=
F(q(Φk))+KL(q(Φk)‖p(Φ|z1:k,Ψk))

(17)

式中:F(·)為logp(z1:k|Φk)的下界函數(shù);q(Φk)為近似后驗分布;p(z1:k|Φk)為真實后驗分布;KL(·)為q(Φk)和p(z1:k|Φk)之間的Kullback-liebler散度,用來衡量兩個概率分布值之間的差異,KL越小,表示兩個分布間的差異越小,且KL非負。

(18)

(19)

VB方法的目標是通過最小化KL來近似真實的后驗值。由于KL是非負的,由式(17)可知這一目標可以通過最大化F(q(Φk))來實現(xiàn)。q(Φk)的最優(yōu)解可由式(20)[28]得到:

logq(Δ)=EΦk(-Δ)[logp(z1:k|Φk)p(Φk|Ψk)]

(20)

將式(21)代入式(20),得

(22)

將式(22)代入式(20)并令Δ=ξk,可得

(23)

由式(23)可知q(ξk)是Gamma分布,

(24)

(25)

(26)

(27)

當選擇i次迭代先驗密度的參數(shù)作為下一次迭代后驗密度的參數(shù)值時,可能會發(fā)生欠擬合[21],在此引入退火算法,設(shè)定退火系數(shù)κ∈[0,1],則

(28)

將式(22)代入式(20)并令Δ=Σk,可得

(29)

由式(29)可知q(Σk)為逆Wishart分布,

(30)

(31)

(32)

(33)

(34)

將式(22)代入式(20)并令Δ=λk,可得

(35)

由式(35)可知q(λk)為Gamma分布,即

(36)

(37)

(38)

(39)

(40)

則

(41)

將式(22)代入式(20)并令Δ=xk,可得

(42)

此時一步預(yù)測PDFp(i+1)(xk|z1:k-1)和測量似然PDFp(i+1)(zk|xk)為

(43)

(44)

(45)

(46)

3 VBMGEKF三維純方位目標跟蹤算法

3.1 修正增益卡爾曼濾波MGEKF

圖2 三維空間目標定位示意圖Fig.2 Schematic diagram of 3D spatial target positioning

(47)

(48)

此時狀態(tài)協(xié)方差矩陣為

(49)

式中:I為單位矩陣;Kk為卡爾曼增益;

(50)

式(50)中Gk的推導(dǎo)詳見附錄B,其中

(51)

(52)

(53)

(54)

(55)

(56)

(57)

3.2 VBMGEKF算法步驟

本文所提出的具有重尾分布噪聲的線性狀態(tài)空間模型的變分貝葉斯修正增益卡爾曼濾波器的一個步長的運算步驟如下:

步驟1輸入k時刻的初始參數(shù):k-1|k-1,Hk,ak,bk,ek,fk,Pk-1|k-1,Fk-1,zk,Qk,Rk,n,τ,L,κ。

步驟2預(yù)測步驟

k|k-1=Fk-1k-1|k-1

(58)

設(shè)定跟蹤模型為勻速運動模型,T為步長,I3為維度為3的單位矩陣。

(59)

(60)

步驟3變分貝葉斯近似更新

fori=1:L

1)初始化:

uk=n+τ+1,Uk=τPk|k-1,k|k=k|k-1,Pk|k=Pk|k-1

2)VB迭代

代入式(24)更新q(i+1)(ξk);

代入式(30)更新q(i+1)(Σk);

代入式(36)更新q(i+1)(λk);

end

3)修正增益卡爾曼濾波MGEKF

(61)

(62)

(63)

(64)

(65)

4 仿真分析

對本文提出的VBMGEKF算法針對純方位目標跟蹤問題進行仿真驗證,仿真模擬單站固定無源系統(tǒng)對機動目標的定位跟蹤,以驗證算法的跟蹤性能和魯棒性;另外將本文提出的VBMGEKF算法與EKF、UKF、PEKF-VB、VBEKF進行比較。

4.1 重尾噪聲建模

構(gòu)建重尾過程噪聲wk與測量噪聲vk,按照文獻[30]設(shè)定過程噪聲協(xié)方差Qk和測量噪聲協(xié)方差矩陣Rk為緩慢時變矩陣,分別為

(66)

(67)

式中:N為仿真時間,此處設(shè)置為150 s;Δt=T;I2為維度為2的單位矩陣;q和g用于調(diào)節(jié)時變噪聲協(xié)方差Qk和Rk的大小,此處設(shè)置q=1,g=100。

設(shè)定非高斯重尾分布的系統(tǒng)過程噪聲以及測量噪聲為

(68)

式中:a為異常值大小的調(diào)節(jié)參數(shù);w.p表示取該項分布的概率。

狀態(tài)值為xk=[rxkrvkrzkvxkvvkvzk]T,設(shè)定x0=[4 000 2 000 10 500 650 150 1]T為初始狀態(tài)值,P=diag[1 0002500250025025021]為初始狀態(tài)協(xié)方差。為比較現(xiàn)有濾波器和本文提出的濾波器的性能,選擇位置和速度誤差的均方根為性能標準,其定義如下:

(69)

(70)

圖3 跟蹤目標軌跡場景示意圖Fig.3 Diagram of target tracking trajectory scene

4.2 仿真校驗

4.2.1a為固定值

分別取a=300、a=500和a=700,此時系統(tǒng)的過程噪聲wk以及測量噪聲vk均為帶異常值的非高斯重尾分布噪聲。仿真結(jié)果如圖4～圖15及表1所示。本文的算法參數(shù)設(shè)置為τ=10-5,L=4,ak=5,bk=5,ek=5,fk=5。其中圖4、圖8、圖12為目標跟蹤距離誤差圖,圖5、圖9、圖13為相對跟蹤誤差,即目標跟蹤誤差與實際距離的比值,圖6、圖10、圖14為跟蹤目標x軸方向的速度誤差圖,圖7、圖11、圖15為跟蹤目標x軸方向的加速度誤差圖,表1為不同噪聲下的各算法的目標跟蹤距離誤差的RMSE。

圖4 目標跟蹤距離誤差RMSE(a=300)Fig.4 RMSE in distance error of target tracking (a=300)

圖5 目標跟蹤相對誤差RMSE (a=300)Fig.5 RMSE inrelative distance error of target tracking (a=300)

圖6 目標跟蹤速度誤差RMSE(a=300)Fig.6 RMSE in velocity error of target tracking (a=300)

圖7 目標跟蹤加速度誤差RMSE (a=300)Fig.7 RMSE in acceleration error of target tracking (a=300)

圖8 目標跟蹤距離誤差RMSE(a=500)Fig.8 RMSE in distance error of target tracking (a=500)

圖9 目標跟蹤相對誤差RMSE (a=500)Fig.9 RMSE in distance relative error of target tracking (a=500)

圖10 目標跟蹤速度誤差RMSE(a=500)Fig.10 RMSE in velocity error of target tracking (a=500)

圖12 目標跟蹤距離誤差RMSE(a=700)Fig.12 RMSE in distance error of target tracking (a=700)

圖13 目標跟蹤相對誤差RMSE (a=700)Fig.13 RMSE in relative distance error of target tracking (a=700)

圖14 目標跟蹤速度誤差RMSE(a=700)Fig.14 RMSE in velocity error of target tracking (a=700)

圖15 目標跟蹤加速度誤差RMSE(a=700)Fig.15 RMSE in acceleration error of target tracking (a=700)

表1 不同噪聲分布下的平均RMSE

從圖4、圖8、圖12中均可以看到,所有算法的跟蹤誤差都隨著時間增大,這是因為目標實際運動軌跡為由近到遠,被測角的角速度減小,目標仍然維持高機動,跟蹤難度上升,因此如圖4、圖8、圖12所示,所有濾波算法的跟蹤距離誤差的變化趨勢均為由小變大,其中本文算法的RMSE相對于其余算法最低。觀察目標跟蹤相對誤差RMSE,各個算法在異常值處于不同水平時的趨勢相同,以圖5為例,可以看到,在40 s之前所有算法的相對誤差均呈現(xiàn)下降趨勢,但是由于EKF算法與PEKF-VB算法將噪聲的先驗分布建模為高斯分布,其在非高斯重尾噪聲下濾波器性能下降嚴重,兩個算法均在40 s后逐漸發(fā)散,且在105 s后PEKF-VB算法發(fā)散更嚴重;UKF算法由于采用確定性樣點來逼近狀態(tài)的后驗

概率密度,具有一定的魯棒性,因此具有比EKF算法、PEKF-VB算法更好的跟蹤性能,但是仍然在40 s之后發(fā)散,但相對誤差相對于EKF算法與PEKF-VB算法更低;本文所提的VBMGEKF算法與文獻[25]提出的VBEKF算法,將系統(tǒng)過程噪聲與測量噪聲均建模為重尾噪聲,本文算法在35 s處逐漸收斂到15%,平均較VBEKF算法低2%左右,可以看出本文的算法性能更佳。另外從圖6、圖10、圖14及其局部放大圖中可以看到,本文算法在一開始對目標速度的估計相對于其他算法就具有更小的誤差。在圖7、圖11、圖15所示的加速度誤差圖中,可以看到本文所提算法與VBEKF算法對加速度的估計相對于其余算法具有更好的效果。

從表1中可以看出,本文算法與VBEKF均具有較小的RMSE,在a=500、a=300和a=700時,即存在不同程度異常值的重尾分布噪聲下具有較好性能,這是因為本文算法與VBEKF算法均將一步預(yù)測的PDF和似然PDF建模為重尾分布,通過VB推斷得到狀態(tài)協(xié)方差與測量噪聲的后驗分布信息,因此比其他算法有更佳的性能,但同時本文還結(jié)合MGEKF,降低量測方程非線性的影響,因此如圖4～圖15所示,本文算法相比其他算法在跟蹤距離誤差、相對跟蹤誤差、目標速度誤差、加速度誤差中均具有較小的RMSE。

圖16 目標跟蹤距離誤差RMSE Fig.16 RMSE in distance error of target tracking

圖17 目標跟蹤相對誤差RMSEFig.17 RMSE in relative distance error of target tracking

圖18 目標跟蹤速度誤差RMSEFig.18 RMSE in velocity error of target tracking

圖19 目標跟蹤加速度誤差RMSEFig.19 RMSE in acceleration error of target tracking

7 結(jié)論

1)本文針對純方位目標跟蹤中存在的非高斯重尾分布噪聲問題提出了一種新的基于高斯分布的變分貝葉斯自適應(yīng)魯棒濾波算法VBMGEKF,通過構(gòu)造基于高斯分布的層次高斯模型來擬合重尾分布的系統(tǒng)過程噪聲與測量噪聲,使用變分貝葉斯近似,迭代得到模型的未知變量和先驗值的參數(shù),推斷出系統(tǒng)狀態(tài)協(xié)方差與測量噪聲協(xié)方差,并且針對純方位目標跟蹤的觀測方程采用MGEKF,降低觀測方程非線性帶來的誤差。

2)經(jīng)過數(shù)值仿真,本文驗證了在純方位目標跟蹤問題中VBMGEKF相對其他非線性濾波EKF、UKF、PEKF-VB、VBEKF在兩種重尾分布噪聲下均具有更高的跟蹤精度與魯棒性。

3)在后續(xù)的研究工作中,將深入探索本文算法在無源多站協(xié)同探測中的應(yīng)用,進一步提高目標跟蹤精度。