有界噪聲條件下基于集員濾波的擴展目標(biāo)跟蹤方法

劉 盼,馬天力,張 榮,李頤果

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

目標(biāo)跟蹤基于最優(yōu)估計原理,采用相關(guān)濾波算法對受到噪聲干擾的量測信息進行處理,從而準(zhǔn)確估計目標(biāo)特征,其廣泛應(yīng)用于智能駕駛、視頻監(jiān)控、防空反導(dǎo)等領(lǐng)域[1-5]。傳統(tǒng)目標(biāo)跟蹤算法主要以“點目標(biāo)”假設(shè)為前提,假定一個目標(biāo)最多產(chǎn)生一個量測。隨著傳感器分辨率的提高,一個目標(biāo)往往可占據(jù)多個分辨單元,即可產(chǎn)生多個量測,該目標(biāo)被稱為擴展目標(biāo)。不同于傳統(tǒng)的點目標(biāo)跟蹤,擴展目標(biāo)跟蹤可以從量測集合中提取目標(biāo)擴展?fàn)顟B(tài)信息(例如目標(biāo)大小、形狀和朝向等[6-7]。因其更貼近實際過程,近幾年吸引了國內(nèi)外學(xué)者的廣泛關(guān)注。

Koch[8]于2008年首次提出基于隨機矩陣模型(random matrix model,RMM)的擴展目標(biāo)跟蹤框架,其將目標(biāo)橢圓輪廓描述為二維對稱正定(symmetric and positive definite,SPD)隨機矩陣,并假設(shè)目標(biāo)散射源均勻散布于目標(biāo)輪廓表面,并利用貝葉斯濾波對擴展目標(biāo)運動狀態(tài)與擴展形態(tài)進行估計。該模型僅考慮目標(biāo)散射源統(tǒng)計特性,未考慮實際過程中存在的傳感器測量誤差對系統(tǒng)的影響。為此,Fledman等[9]提出混合加性量測噪聲表示模型,其主要思想是量測的散布程度受到擴展目標(biāo)形態(tài)與實際觀測噪聲共同作用,具體表現(xiàn)為兩者線性組合的形式,但是其難以采用基于貝葉斯理論的濾波方法對其狀態(tài)與擴展形態(tài)的后驗概率密度進行求解。針對這一問題,蘭劍等[10-11]提出量測噪聲近似策略,通過構(gòu)建形態(tài)觀測矩陣對混合加性量測噪聲進行近似。Liu[12]則構(gòu)建乘性誤差模型(multiplicative error model,MEM),該模型通過引入高斯乘性噪聲項對目標(biāo)散射源分布進行描述,進而利用擴展卡爾曼濾波(extended Kalman filter,EKF)和二階擴展卡爾曼濾波(second-order EKF,SO-EKF)等[13-15]方法對擴展目標(biāo)狀態(tài)和形態(tài)進行估計。

在實際應(yīng)用過程中,由于傳感器系統(tǒng)受到內(nèi)部熱噪聲和外部信號擾動影響,量測噪聲可能存在非高斯性,若采用上述方法對擴展目標(biāo)進行估計,極易導(dǎo)致目標(biāo)跟蹤精度降低。鑒于此,Zhang等[16]將量測噪聲通過偏正態(tài)分布進行表示,利用變分貝葉斯推理計算系統(tǒng)后驗概率密度函數(shù)。LI Yawen、Gao Lei和陳輝等提出基于Student’s t分布的量測噪聲統(tǒng)計模型,并分別利用變分推理、魯棒Student’s t濾波算法對擴展目標(biāo)運動狀態(tài)與形態(tài)進行計算[17-19]。

現(xiàn)有的擴展目標(biāo)跟蹤方法均是基于概率框架,其假設(shè)擴展目標(biāo)狀態(tài)以及量測服從某一特定分布,并利用貝葉斯濾波方法對其目標(biāo)狀態(tài)與形態(tài)進行估計。事實上,因目標(biāo)加速度物理特性,外界未知環(huán)境不確定干擾,使得量測噪聲具有未知但有界特性(unknown but bounded,UBB),難以運用概率框架下的相關(guān)濾波算法進行求解。因此,針對UBB噪聲條件下擴展目標(biāo)狀態(tài)估計問題,本文提出一種基于集員濾波的擴展目標(biāo)跟蹤方法(extended object tracking based on set membership filter,SMF-EOT),通過UBB橢球集合對量測噪聲進行表示,利用集員濾波方法對狀態(tài)集合參數(shù)進行計算,在形態(tài)更新過程中,結(jié)合Graham scan策略,求得包含目標(biāo)形態(tài)最大誤差的最小邊界矩陣,通過仿射變換和偏移超曲面對擴展目標(biāo)形態(tài)進行更新,最終獲得擴展目標(biāo)運動狀態(tài)與擴展形態(tài)。

1 問題描述

建立擴展目標(biāo)狀態(tài)空間模型:

(1)

Ck=λXk⊕Uk

(2)

式中:λ為散射因子,表示目標(biāo)擴展?fàn)顟B(tài)對量測值的影響程度,λ∈(0,1];Xk為k時刻目標(biāo)形態(tài)矩陣;Uk為傳感器量測誤差邊界;⊕表示Minkowski和。

對于上述擴展目標(biāo)跟蹤系統(tǒng),在UBB噪聲條件下如何對目標(biāo)運動和擴展?fàn)顟B(tài)進行估計將是本文解決的核心問題。

2 基于集員濾波的擴展目標(biāo)跟蹤方法

2.1 集員濾波相關(guān)基礎(chǔ)理論

對UBB噪聲條件下的擴展目標(biāo)進行跟蹤,其主要思路是利用基于區(qū)間數(shù)學(xué)理論的集員估計方法[20],下面首先介紹集合相關(guān)定義及其運算性質(zhì)。

定義1橢球集合E(x,S)表示為:

E(x,S)={y∈Rn|(y-x)TS-1(y-x)≤1}

(3)

式中:x為橢球中心;S為SPD矩陣,表示橢球形狀大小。

定義2橢球集合E(x,S)的支持函數(shù)[20]表示為:

η(E(x,S))=αTx+(αTSα)1/2,α∈Rn

(4)

定理1橢球集合E(x1,S1)與E(x2,S2)支持函數(shù)的Minkowski和[20]為:

η(E(x1,S1)⊕E(x2,S2))=

η(E(x1,S1))+η(E(x2,S2))

(5)

雖然定理1給出了E(x1,S1)和E(x2,S2)支持函數(shù)的Minkowski和,但其并不是一個確定大小的橢球。因此,需找到包含Minkowski和的外定界橢球E(x,S),如圖1所示。

圖1 外定界橢球

因此,需要根據(jù)定理2和定理3計算最優(yōu)準(zhǔn)則下的外定界橢球E(x,S)。

定理2E(x,S)為橢球集合,Σ為已知n維方陣,則:

ΣE(x,S)=E(Σx,ΣSΣT)

(6)

證明:利用支持函數(shù)可將E(x,S)表示為:

(7)

定理3已知橢球集合E(x1,S1)和E(x2,S2),包含2個橢球的Minkowski和的外定界橢球E(x,S)可表示為:

E(x,S)=E(x1,S1)⊕E(x2,S2)=

E(x1+x2,S(p))

(8)

其中:

(9)

證明考慮橢球集合E(x1,S1)和E(x2,S2),根據(jù)橢球定義:

i=1,2

(10)

假定外接橢球集合為E(x,S)。

E(x,S)={y∈Rn|(y-x)TS-1(y-x)≤1}

(11)

若外接定界橢球集合E(x,S)能夠包含2個橢球集合,則其支持函數(shù)必須滿足下列不等式:

η(E(x,S))≥η(E(x1,S1))+η(E(x2,S2))

(12)

根據(jù)定義2,上式可變換為:

αTx+(αTSα)1/2≥αTx1+(αTS1α)1/2+

αTx2+(αTS2α)1/2

(13)

則外定界橢球中心x為:

x=x1+x2

(14)

根據(jù)式(13),將外界橢球邊界S通過S1和S2的線性組合表示為:

γαTS1α+ραTS2α≥αTS1α+2(αTS1α)1/2·(αTS2α)1/2+αTS2α

(15)

(16)

當(dāng)(γ-1)(ρ-1)≥1時,式(16)成立,可行標(biāo)量p需滿足γ-1=p-1,ρ-1≥p,則:

γ=1+p-1,ρ=1+p

(17)

即S(p)可以表示為:

S(p)=(1+p-1)S1+(1+p)S2,p>0

(18)

通過定理3可知,2個橢球集合Minkowski和的外定界橢球是形狀矩陣S關(guān)于參數(shù)p的函數(shù)。根據(jù)以下定理計算最優(yōu)橢球參數(shù)p。

定理4已知橢球集合E(x1,S1)和E(x2,S2),則半軸平方和最小意義下最小跡橢球參數(shù)p。

(19)

證明最小化橢球E(x,S(p))的最小跡等價于求取如下函數(shù)的最小值。

f(p)=tr(S(p))=tr((1+p-1)S1+(1+p)S2)

(20)

將式(20)對變量p求導(dǎo)并令其導(dǎo)數(shù)等于零,即可求得函數(shù)f(p)極值時p的值:

(21)

2.2 基于集員濾波的擴展目標(biāo)跟蹤方法

令擴展目標(biāo)狀態(tài)集合為χk-1?E(xk-1,Sk-1),xk-1為狀態(tài)集合中心,Pk-1為協(xié)方差矩陣,Sk-1表示狀態(tài)集合邊界。由于擴展目標(biāo)狀態(tài)向量為集合形式,則目標(biāo)運動狀態(tài)估計由點估計變?yōu)闋顟B(tài)可行集的估計?；诩瘑T濾波的擴展目標(biāo)跟蹤方法對目標(biāo)狀態(tài)與形態(tài)的估計分別包括預(yù)測步驟和更新步驟。

2.2.1 運動狀態(tài)預(yù)測

在卡爾曼濾波基礎(chǔ)上,狀態(tài)集合一步預(yù)測為:

χk,k-1=Φkχk-1⊕wk

(22)

預(yù)測協(xié)方差誤差陣為:

(23)

根據(jù)定理2,結(jié)合狀態(tài)空間模型,可得橢球集合Φkχk-1的支持函數(shù)如下:

(24)

由定理1可得狀態(tài)集合預(yù)測的支持函數(shù):

η(E(xk,k-1,Sk,k-1))=

(25)

要使E(xk,k-1,Sk,k-1)能夠包含式中2個橢球的Minkowski和,則其支持函數(shù)必須滿足:

η(E(xk,k-1,Sk,k-1))≥

(26)

根據(jù)定理2和定理3,包含狀態(tài)集合χk,k-1的外定界橢球E(xk,k-1,Sk,k-1)的中心值xk,k-1及橢球大小的矩陣Sk,k-1為:

(27)

式(27)中,需計算最小跡橢球參數(shù)pSk,k-1使橢球E(xk,k-1,Sk,k-1)為包含橢球ΦkE(xk-1,Sk-1)和E(0,Qk)的Minkowski和最小外定界橢球。根據(jù)定理4可計算得到關(guān)于可行標(biāo)量pSk,k-1的最小跡函數(shù)和半軸平方和最小意義下最小跡橢球參數(shù)pSk,k-1。

(28)

(29)

2.2.2 擴展形態(tài)預(yù)測

對于擴展形態(tài)Xk的預(yù)測,假設(shè)其在時間間隔內(nèi)目標(biāo)形態(tài)大小不發(fā)生改變[9],則k時刻擴展形態(tài)的預(yù)測Xk,k-1為k-1時刻擴展目標(biāo)狀態(tài)Xk-1的更新結(jié)果,即:

Xk,k-1=Xk-1

(30)

2.2.3 運動狀態(tài)更新

(31)

式中:濾波增益Kk及更新誤差協(xié)方差陣Pk分別為:

(32)

(33)

根據(jù)定理1和定理2,結(jié)合式(1),可得狀態(tài)更新集合χk的支持函數(shù):

(34)

類比于時間更新步驟,根據(jù)定義2和定理3可得包含狀態(tài)更新集合χk的橢球E(xk,Sk)的中心值xk和橢球形狀矩陣Sk。

(35)

同樣,根據(jù)定理4可得關(guān)于可行標(biāo)量pSk的函數(shù)表示為:

(36)

則最小跡橢球參數(shù)pSk為:

(37)

2.2.4 擴展形態(tài)更新

根據(jù)定理3和定理4,則第r個量測橢球集合邊界為:

(38)

(39)

因量測橢球集合邊界由目標(biāo)擴展?fàn)顟B(tài)以及傳感器量測誤差邊界信息構(gòu)成。因此可將其作為目標(biāo)擴展形態(tài)誤差邊界,即包含k時刻所有量測橢球集合的最小邊界為目標(biāo)擴展形態(tài)最大誤差邊界。由于利用橢球集合Minkowski和無法求得包含所有量測橢球的最小外接橢球,因此利用凸包計算幾何中Graham scan算法計算包含k時刻所有量測橢球集合的最小邊界集合A。

(40)

式中:dot(·)為向量點積運算符。

根據(jù)極角τi大小順序?qū)螶進行排序得到Js。將Js中的每一個點ai(i≥3)與其余各點aj在二維坐標(biāo)平面中求向量叉積σi,j=ai×aj。若σi,j≤0,表明ai是最小邊界上的點;若σi,j>0,表明ai不是最小邊界上的點。經(jīng)過上述計算可得到包含所有量測橢圓的最小邊界點的集合A=[a1,a2,…,am-1,am]。以a3與a2,a1為例,如圖2所示。

(a)

(41)

圖3 E1與E2橢球Minkowski差

為求得E1?E2,首先利用Givens Rotation[23]對E2進行分解,計算橢圓長短半軸aE2,bE2與旋轉(zhuǎn)矩陣RE2,其次利用線性變換對E1進行參數(shù)化。

(42)

式中:Di為橢球邊界任意一點;ΛE1=diag(a1,b1)。

(43)

式中:ΛE2=diag(aE2/bE2,1)。將式(43)代入式(44)中,可得到經(jīng)仿射變換后E1橢球邊界的參數(shù)表達式:

(44)

為了獲得內(nèi)切于E1的各圓E2的中心,利用超曲面gofs(φ)對仿射變換后E1?E2進行計算。

(45)

所求的偏移超曲面gofs(φ)為仿射變換后E1與E2Minkowski差的邊界。如圖3所示,因仿射變換后其形狀大小和方向等狀態(tài)發(fā)生變化,需對gofs(φ)進行仿射逆變換,使E1回到初始狀態(tài)。則E1?E2邊界的封閉解為:

(46)

根據(jù)中心xk到邊界geb(φ)的距離計算E1與E2的Minkowski差所表示橢圓的長短半軸和旋轉(zhuǎn)矩陣,進而可得表示E1與E2的Minkowski差橢圓的矩陣Xd,由式可知Xd可近似為集合λXk,k和Uk的Minkowski和,Xd=λXk,k⊕Uk。為了求解擴展目標(biāo)形態(tài)更新矩陣Xk,k,則需要計算Xd與誤差邊界Uk的Minkowski差。即計算Xd?Uk邊界的封閉解geb,Xd(φ)。

(47)

式中:RUk為Uk的旋轉(zhuǎn)矩陣;gofs,Xd表示仿射變換后Xd?Uk的邊界表達式。

(48)

式中:Dofs(φ)為經(jīng)仿射變換后橢圓Xd關(guān)于角度φ的點集,通過中心xk到geb,Xd(φ)的距離對擴展目標(biāo)形態(tài)更新橢圓長短軸和角度進行計算,獲得更新后的擴展目標(biāo)形態(tài)Xk,k。

3 仿真實驗

為了驗證所提基于集員濾波的擴展目標(biāo)跟蹤方法(extended object tracking based on set membership filter,SMF-EOT)的有效性,首先采用RMF[8]、MEM-EKF[14]、MEM-SOEKF[15]在UBB噪聲條件下對擴展目標(biāo)跟蹤性能進行對比驗證。其次,采用不同有界噪聲參數(shù)對本文所提算法跟蹤性能進行仿真實驗分析。

3.1 機動場景仿真實驗

(49)

式中:Θ為系統(tǒng)機動時間常數(shù),Θ=40 s;過程噪聲集合邊界Qk=diag[0.52;0.52]。觀測矩陣Hk=[1 0 0],量測噪聲集合Ck=λXk⊕Uk,其中λ=0.25為噪聲參數(shù)散射因子,量測誤差邊界Uk=diag{202,102}。采樣間隔T=10 s。擴展目標(biāo)在時間t∈[30,70]和t∈[90,130]分別做角速率為4.5 rad/s和-4.5 rad/s的轉(zhuǎn)彎運動。

假設(shè)擴展目標(biāo)初始狀態(tài)橢球χ0=E(x0,S0),目標(biāo)初始位置x0=[800,-200,9.82,-9.82,0,0]T,S0為初始橢球形狀大小矩陣,初始狀態(tài)協(xié)方差矩陣P0=diag[702,702,102,102,52,52],系統(tǒng)仿真時間為160 s。圖4為4種算法一次實驗的擴展目標(biāo)跟蹤結(jié)果,可以看出在UBB噪聲條件下,本文所提SMF-EOT算法性能優(yōu)于RMF、MEM-EKF、MEM-SOEKF算法。

圖4 4種算法擴展目標(biāo)跟蹤結(jié)果

圖5為4種算法的位置和速度RMSE的100次Monte-Carlo實驗仿真結(jié)果?？梢钥闯?與MEM-EKF、MEM-SOEKF和RMF算法相比,本文所提出的SMF-EOT算法具有更小的位置RMSE。當(dāng)目標(biāo)機動時,MEM-EKF和RMF的位置和速度RMSE增大,主要是由于MEM-EKF和RMF算法假設(shè)的噪聲統(tǒng)計特性與實際噪聲不匹配,從而導(dǎo)致其位置和速度估計精度下降。雖然MEM-SOEKF算法通過二階泰勒級數(shù)展開對非線性進行近似,在目標(biāo)機動時估計精度并未下降,但該方法同樣僅適用于噪聲服從高斯分布的擴展目標(biāo)跟蹤系統(tǒng),在UBB噪聲條件下其位置和速度RMSE仍高于SMF-EOT算法。

(a)位置

表1為4種算法的平均均方根誤差(average root mean square error,ARMSE)。可以看出,所提算法的位置ARMSE相比與RMF、MEM-EKF和MEM-SOEKF算法分別提高了69.2%、77.05%、和65.74%;其速度ARMSE分別提高了49.69%、49.92%和48.7%。

表1 4種算法的位置和速度ARMSE對比

為對四種算法目標(biāo)擴展形態(tài)估計的準(zhǔn)確度進行評價,采用高斯威斯頓距離(gaussian wasserstein distance,GWD)作為評價指標(biāo)[23],該距離通過目標(biāo)位置和形狀誤差計算目標(biāo)擴展?fàn)顟B(tài)誤差,其計算公式如下:

(50)

圖6為4種算法100次Monte-Carlo實驗下的GW距離。從圖中可以看出,相比于MEM-EKF、MEM-SOEKF和RMF算法,本文所提算法具有更小的GW距離,且RMF算法估計性能較差。主要原因在于RMF算法將擴展形狀建模為隨機矩陣,導(dǎo)致迭代更新過程中無法處理橢圓長短半軸與角度之間的不確定性導(dǎo)致對目標(biāo)擴展?fàn)顟B(tài)估計精度較低。雖然MEM-EKF和MEM-SOEKF比RMF算法對目標(biāo)擴展?fàn)顟B(tài)估計精度較高,但在UBB噪聲條件下受噪聲適應(yīng)性的影響其估計精度低于本文所提SMF-EOT算法。

圖6 4種算法的GW距離

表2為4種算法的(average Gaussian wasserstein distance,AGWD)對比,所提算法的AGWD相對于RMF、MEM-EKF和MEM-SOEKF分別提高了49.59%、54.43%和36.73%。

表2 4種算法的AGWD對比

3.2 參數(shù)影響分析

本文所提算法假設(shè)量測噪聲為散射源與誤差量測噪聲集合的Minkowski和,初始參數(shù)散射因子λ以及量測誤差噪聲邊界U可能會影響算法性能。

考慮量測誤差噪聲邊界U=diag{202,102}時,散射因子λ分別為0.25、0.5和1對本文所提算法進行100次Monte-Carlo實驗;散射因子λ=0.25時,量測誤差噪聲邊界U分別為diag{202,102}、diag{1002,502}和diag{1502,1002}對本文所提算法進行100次Monte-Carlo實驗。

圖7為不同噪聲影響參數(shù)λ下SMF-EOT算法位置RMSE和GW距離對比。從圖中可以看出,影響參數(shù)λ值越小,SMF-EOT算法收斂性越好。

(a)位置

圖8為量測誤差噪聲邊界U下SMF-EOT算法位置RMSE和GW距離對比?？梢钥闯隽繙y誤差噪聲邊界U越小,量測誤差邊界越小,SMF-EOT算法對目標(biāo)運動和形態(tài)估計精度越高。

(a)位置

綜上所述,在UBB噪聲條件下,與MEM-EKF、MEM-SOEKF和RMF算法RMSE相比,本文所提SMF-EOT算法具有更小的位置、速度和橢圓GW距離,估計性能更好。主要原因在于SMF-EOT算法考慮噪聲邊界已知但其統(tǒng)計特性未知的擴展目標(biāo)跟蹤系統(tǒng),假設(shè)系統(tǒng)噪聲為UBB噪聲并將其建模為橢球集合,利用集員估計理論以及橢球Minkowski差對目標(biāo)運動和擴展?fàn)顟B(tài)進行估計。然而MEM-EKF、MEM-SOEKF和RMF算法使用貝葉斯規(guī)則迭代估計目標(biāo)狀態(tài),噪聲假設(shè)局限于高斯等隨機分布,沒有考慮到有界噪聲的擴展目標(biāo)跟蹤問題導(dǎo)致其跟蹤性能降低。

4 結(jié)語

本文針對有界噪聲條件下的橢圓擴展目標(biāo)跟蹤問題,提出基于集員濾波的橢圓擴展目標(biāo)跟蹤方法。其將系統(tǒng)噪聲建模為橢球集合噪聲,采用集員濾波方法對目標(biāo)運動狀態(tài)進行估計,在對目標(biāo)擴展?fàn)顟B(tài)估計時利用Minkowski和理論獲取量測橢圓并用Graham scan算法對其進行融合,同時結(jié)合橢圓封閉形式的Minkowski差求解得到目標(biāo)擴展?fàn)顟B(tài)。數(shù)值模擬仿真實驗結(jié)果表明,在基于有界噪聲假設(shè)的擴展目標(biāo)跟蹤系統(tǒng)中,本文所提算法對擴展目標(biāo)運動狀態(tài)和形態(tài)具有較高的估計精度。未來的研究方向可以考慮解決有界噪聲條件下非凸形狀擴展目標(biāo)的跟蹤問題。