面向非線性MTT的多模型泊松多伯努利混合濾波算法

陳嵩杰,李 波,張 露

(遼寧工業大學 電子與信息工程學院,遼寧 錦州 121001)

0 引 言

多目標跟蹤[1](Multi Target Tracking,MTT)是一個重要的研究方向,適用于航空航天、空中交管和自動駕駛等領域[2,3].考慮MTT共軛性所反映的目標先驗與后驗分布概率密度具有相同形式的特點[4],基于多伯努利共軛先驗理論的帶標簽廣義標簽多伯努利(Generalized Label Multi-Bernoulli,GLMB)濾波[5]和不帶標簽的泊松多伯努利混合(Poisson Multi-Bernoulli Mixture Filter,PMBM)濾波[6]已有研究.其中,PMBM濾波可為目標概率密度提供精確閉式解,計算代價和精度較GLMB濾波器更優.文獻[7]提出了一種多模型多伯努利濾波算法,給出了該算法在線性高斯模型的實現;文獻[8]提出了一種多伯努利混合濾波的高斯實現;文獻[9]提出了一種基于多模型泊松多伯努利混合濾波的高斯逆威沙特實現,利用高斯逆威沙特分布估計目標質心,并引入強跟蹤濾波修正狀態模型協方差矩陣.

針對非線性模型的MTT問題,現有文獻主要采用序貫蒙特卡洛(Sequential Monte Carlo,SMC)和高斯混合(Gaussian Mixture,GM)兩種方式實現.其中,SMC實現[10,11]存在粒子退化問題,濾波計算繁瑣,且跟蹤速度慢;GM實現[12]的計算復雜程度相對較低,適用于線性高斯模型的MTT.基于此,文獻[13]提出了一種非線性模型的擴展卡爾曼濾波(Extended Kalman Filter,EKF)的PMBM算法,但EKF在強非線性系統中存在較低的魯棒性,其雅克比矩陣計算也較為困難;文獻[14]提出了一種非線性模型的(Cubature Kalman Filter-Multi-Bernoulli,CKF-MB)和(Square-rooted CKF-MB,SCKF-MB)濾波算法,但在低檢測概率時,目標勢估計存在一定偏差;文獻[15]提出了一種用于非線性模型δ-GLMB濾波的積分卡爾曼高斯混合實現,但有效性不高.針對非線性高雜波環境的MTT,文獻[16]提出了一種均方根容積卡爾曼濾波與δ-GLMB結合的GM實現,但仍存在平衡濾波運算時間和跟蹤精度難以平衡的問題.

針對MTT的非線性與低檢測概率問題,本文在非線性模型下融合SCKF與高斯混合多模型泊松多伯努利混合濾波(GM Multiple Model PMBM,GM-MM-PMBM),提出一種基于SCKF的GM-MM-PMBM濾波算法.首先,推導出MM-PMBM預測和更新的GM實現.然后,采用誤差協方差的平方根方式進行遞推,應用等權值容積點集計算多目標密度函數的均值與協方差矩陣.最后,結合SCKF對GM-MM-PMBM濾波高斯項進行預測和更新,采用最優子模式[17,18](Optimal Sub-Pattern Assignment,OSPA)距離評價該算法的性能.

1 相關工作

1.1 隨機有限集

在隨機有限集(Random Finite Set,RFS)理論框架中[19],假設單目標狀態x∈nx,多目標狀態X∈F(nx),X為單目標狀態向量的集合,F(nx)為nx空間中所有有限子集集合,目標狀態可由量測集合Z∈F(nz)中的量測觀察.

假設k時刻的狀態集為X={x1,…,xn};測量集為Z={Z1,…,Zn}∪κ(z),κ(z)為雜波量測集合且滿足κ(z)=λc(z),λ為泊松率,c(z)為雜波空間概率分布函數,Zi為目標i產生的量測集合,κ(z)與Z1…Zn相互獨立.

1.2 PMBM共軛先驗

PMBM濾波可由泊松和多伯努利混合兩部分表示[20].其中,泊松部分用于表示未檢目標;多伯努利混合部分則表示已檢目標[21].假設X=Xu+∪Xd為目標集合,Xu和Xd分別為未檢目標和已檢目標,則共軛先驗可定義為:

(1)

(2)

(3)

式中,wi,a為權重,+∪表示不相交并集;fp(·)為泊松部分概率密度函數;μ(·)為泊松強度;fmbm(·)為多伯努利混合部分概率密度函數;A為全局假設索引集合,a為所有全局假設索引,n為所有伯努利分量的數目.

于是,第a個全局假設中第i個伯努利量為:

(4)

式中,ri,a(·)和fi,a(·)分別為第a個全局假設中第i個的目標存在概率和伯努利分量.

1.3 標準PMBM算法

1.3.1 預測過程

由式(1)～式(3)可知,PMBM濾波器的預測過程包含泊松部分和多伯努利混合部分.假設k-1時刻泊松強度為μk-1(xk-1),則k時刻泊松強度預測密度函數為:

μk|k-1(xk)=λb(xk)+〈μk-1(xk-1),pSf(xk|xk-1)〉

(5)

式中,λb(·)為新生強度,pS為目標存活概率,f(xk|xk-1)為狀態轉移函數,[·,·]表示函數的內積.

同時,多伯努利部分預測參數可表示為:

(6)

1.3.2 更新過程

更新過程包含3個部分:未檢目標的泊松強度更新、首次檢測目標的泊松強度更新、已檢目標的伯努利更新[22].

未檢目標的泊松強度更新為:

(7)

首次檢測目標的伯努利更新為:

(8)

(9)

已檢目標(漏檢目標和量測目標)的伯努利更新為:

(10)

(11)

可以看出,泊松強度更新和伯努利更新過程是分別進行的,首次檢測目標的伯努利更新和已檢目標的伯努利更新差異是前者中的伯努利分量含泊松分量.

2 多伯努利混合濾波非線性實現

由于標準PMBM算法在遞歸步驟中存在難以求解RFS積分的問題[23],因此基于高斯線性模型推導出該算法的GM實現,將積分形式傳遞的伯努利參數轉化為高斯項加權求和形式.

為表述簡潔,下文中的各函數主要以變量ξ定義,且假設檢測概率pd(·)和存活概率pS(·)與模型相關,滿足pd(xk,ξ)=pd(ξ),pS(xk,ξ)=pS(ξ).于是,PMBM算法的GM實現的預測和更新遞推如下.

2.1 GM-MM-PMBM算法

2.1.1 預測過程

基于高斯分布N(·),新生目標泊松強度的GM形式可記為:

(12)

假設k-1時刻泊松強度的GM形式為:

(13)

則k時刻預測強度的GM形式為:

(14)

式中,Fk(ξ)為模型ξ的狀態轉移矩陣,Qk(ξ)為模型ξ的過程噪聲協方差矩陣,fk|k-1(ξ|ξ′)為模型ξ′～模型ξ的轉移概率.

假設k-1時刻伯努利密度函數為:

(15)

則預測過程伯努利分量的權重保持不變,其存在概率滿足:

(16)

于是,多伯努利部分預測密度函數為:

(17)

2.1.2 更新過程

未檢目標的泊松強度更新:假設泊松預測強度GM形式由式(19)給出.

(18)

則更新后的后驗強度為:

μk|k(xk,ξ)=(1-pd)μk|k-1(xk,ξ)

(19)

首次檢測到目標的泊松強度更新:

(20)

(21)

式中,相關參數定義如下:

已檢目標的伯努利更新:對于漏檢目標定義下述公式.

(22)

(23)

(24)

對于量測目標則有:

(25)

(26)

(27)

式中,相關參數定義如下:

2.2 基于SCKF的GM-MM-PMBM算法

考慮運動模型ξ,假設目標狀態方程和量測方程分別為:

(28)

式中,xk和zk分別為狀態向量和量測向量;fk(·)和hk(·)分別為狀態函數和量測函數;uk-1是均值為0、協方差矩陣為Qk-1的過程噪聲;vk是均值為0、協方差矩陣為Rk的量測噪聲.

假設k-1時刻誤差協方差的平方根滿足下三角陣Γk-1|k-1=Tria(A),Tria(·)為矩陣的三角化運算,于是基于SCKF的非線性濾波過程如下:

2.2.1 預測過程

由Cholesky分解計算預測高斯項的協方差矩陣:

(29)

利用三階球面-徑向容積準則生成容積點χq,k-1及權重wq=1/2n,q=1,2,…,2n,這里系統狀態維數設為2n.于是,計算非線性狀態函數的傳遞容積點:

(30)

計算狀態預測均值:

(31)

計算預測誤差協方差矩陣的平方根:

(32)

2.2.2 更新過程

產生容積點χq,+并計算非線性量測函數傳遞容積點:

(33)

計算量測預測均值:

(34)

計算新息協方差平方根:

(35)

計算互協方差矩陣:

(36)

計算卡爾曼增益:

(37)

計算狀態均值更新:

(38)

計算協方差矩陣的平方根更新:

(39)

3 仿真驗證與結果分析

3.1 仿真情形及參數設置

目標運動模型ξ設為勻速(Constant Velocity,CV)運動和勻速轉彎(Constant Turn,CT)運動,其狀態轉移方程分別為:

式中,wk-1是均值為0、協方差矩陣為Qk-1的過程噪聲.

目標非線性量測模型方程為:

觀測矩陣和量測噪聲協方差矩陣分別為:

式中,標準差σε=10,I2為二維單位方陣.此外,目標存活概率為pS=0.99,由泊松分布對雜波建模.模型1定義為CV運動,設過程噪聲為5m/s2;模型2定義為CT運動,設過程噪聲為5 m/s2,逆時針轉動角度為3°/s;模型3定義為CT運動,設過程噪聲為5m/s2,順時針轉動角度為3°/s.于是,定義CV和CT運動的線性狀態轉移矩陣如下:

假設目標新生模型為泊松形式,新生目標泊松強度GM為:

所有存活目標在前20s進行CV運動,在第21～ 40s內按逆時針方向由CV運動轉換為CT運動,在第60～80s內按順時針方向轉換為CT運動,模型間轉移概率矩陣定義為:

圖1給出了檢測區域和3個目標的非線性運動軌跡.

圖1 目標真實軌跡Fig.1 Target real trajectory

3.2 實驗結果分析

為比較本文SCKF-GM-MM-PMBM算法與EKF-PMBM算法[13]、SCKF-GM-MM-MB算法[14]、SMC-MM-PMBM算法的跟蹤性能,本實驗進行100次SMC仿真.

3.2.1 固定雜波密度及高檢測概率情形下

假設檢測概率為pd=0.96,雜波密度為λC=0.6×10-3/m2,目標勢估計及誤差估計仿真結果如圖2所示.

圖2 高檢測概率下勢估計Fig.2 Cardinality estimation under high detection probability

可以看出,隨著時間的變化,各算法都能較好完成勢估計.EKF-PMBM算法、SCKF-GM-MM-MB算法和本文算法都能接近目標真實數目,而經典SMC-MM-PMBM算法在粒子數目較少時,會出現低于目標真實值的估計情形;并在第15s和第30s時,SMC-MM-PMBM算法估計時延較長且對目標勢估計有明顯偏差,EKF-PMBM算法和SCKF-GM-MM-MB算法的目標數目估計時延較低,而本文算法的估計時延和消耗較最低,能有效地估計目標勢變化.

本文采用OSPA距離作為MTT性能的評價指標,該距離越大表示算法的總體精度越差,設定參誤差調節參數c為100,階次參數p為1[24].各濾波算法的OSPA距離對比如圖3所示.可以看出,各算法具有相似的OSPA誤差,在第5s和第18s時,SCKF-GM-MM-MB算法給出的OSPA距離略優于SCKF-GM-MM-PMBM算法.若目標在第20s出現,在第60s和第80s時目標消失,本文算法優于SCKF-GM-MM-MB算法,其原因為PMBM算法中泊松部分保留了漏檢分量和剛消失的目標分量,能及時檢測新生目標并延遲對消亡目標的響應.EKF-PMBM算法和SMC-MM-PMBM算法的OSPA距離相比于另兩種算法偏大,其原因為強非線性下EKF需先進行線性處理,忽略了高階項,出現發散問題,導致跟蹤目標丟失;而SMC的采樣粒子較多易發生粒子退化現象,若當采樣粒子數目較少時,目標跟蹤性能和運算速度降低.

圖3 OSPA距離對比Fig.3 OSPA distance comparison

3.2.2 固定雜波密度及低檢測概率情形下

假設目標檢測概率為pd=0.6,仿真結果如圖4所示.

圖4 低檢測概率下勢估計Fig.4 Cardinality estimation under low detection probability

對比高檢測概率環境,可以看出SCKF-GM-MM-PMBM算法的勢估計與SCKF-GM-MM-MB算法、EKF-PMBM算法和SMC-MM-PMBM算法有著明顯的差異,SCKF-GM-MM-PMBM算法具有更好的勢估計誤差.

各算法的OSPA距離對比如圖5所示.可以看出,當第62s目標消亡時,本文算法的OSPA距離誤差較SCKF-GM-MM-MB算法大,其原因為較低的Pd使本文算法在更新過程中的泊松分量權重增大,無法及時刪除剛剛消失的目標狀態.

圖5 OSPA距離對比Fig.5 OSPA distance comparison

表1綜合比較了OSPA距離誤差和單步平均運行時間.可以看出,本文算法在低檢測概率時下具有更好的跟蹤性能.

表1 不同檢測概率的OSPA距離和單步平均運行時間Table 1 OSPA distance and single-step average running time under different detection probabilities

3.2.3 不同雜波密度情形下

為進一步驗證所提算法的跟蹤性能,給出了各算法在不同雜波密度λC=(0.6,1.2,1.8,2.4,3.0)×10-3/m2時的勢估計誤差與OSPA距離對比.

圖6和圖7分別為勢估計誤差對比圖與OSPA距離對比圖.可以看出,當雜波密度增大時,各算法的勢估計誤差與OSPA距離誤差增大,EKF-PMBM算法僅次于SCKF-GM-MM-MB算法,其原因為弱非線性環境下EKF-PMBM算法可優先考慮,在強非線性情形下EKF算法與線性假設原則沖突,導致濾波發散并產生估計誤差,但本文算法在不同雜波密度情況下,其勢估計誤差和OSPA距離誤差都要優于其他3種算法.

圖6 勢估計誤差對比Fig.6 Cardinality estimation and error comparison

圖7 OSPA距離對比Fig.7 OSPA distance comparison

表2對比了當檢測概率pd=0.96時不同雜波密度環境下的OSPA距離與單步平均運行時間.可以看出,當同一檢測概率時,各算法的OSPA距離和運行時間都隨雜波密度的增加而增大.當不同雜波時,本文算法具有較高跟蹤精度,運行時間低于其他算法.當雜波密度較高時,由于EKF-PMBM算法進行雅克比矩陣計算,導致運行時間增加;而SMC-MM-PMBM算法所需的采樣粒子數目不斷增加,導致運行時間增加.

表2 不同雜波密度下的OSPA距離和單步平均運行時間Table 2 OSPA distance and single-step average running time under different clutter densities

綜上,本文算法在不同檢測概率情形下能以較高的精度估計目標的勢和運動狀態.尤其是在低檢測概率情形下,其跟蹤精度和運行時間優于其他3種算法;在高檢測概率不同雜波密度情形下,各算法的OSPA距離與運行時間具有一定差距,本文算法的OSPA距離優于其他3種算法.

4 總 結

本文針對MTT存在的非線性與低檢測概率問題,提出了MM-PMBM算法的GM實現,在低檢測概率環境下為目標概率密度提供精確閉式解,其跟蹤精度和運行時間上相比MM-MB算法有所提高;在非線性系統模型下,融合SCKF濾波框架,避免了計算協方差矩陣的平方根,其數值穩定性和濾波性能得以保證;在不同雜波密度情形下,該算法可準確估計目標的勢與運動狀態,提高了MTT的跟蹤精度.由不同情形下的仿真實驗結果可知,本文算法在非線性情形中OSPA距離誤差與運行時間均最小.

針對MTT的運行時間,接下來將深入研究本文算法用于解決機動目標、擴展目標和群目標等比較復雜的問題,進一步節省運行時間.