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

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

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

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

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

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

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

1 問題描述

建立擴展目標狀態空間模型:

(1)

Ck=λXk⊕Uk

(2)

式中:λ為散射因子,表示目標擴展狀態對量測值的影響程度,λ∈(0,1];Xk為k時刻目標形態矩陣;Uk為傳感器量測誤差邊界;⊕表示Minkowski和。

對于上述擴展目標跟蹤系統,在UBB噪聲條件下如何對目標運動和擴展狀態進行估計將是本文解決的核心問題。

2 基于集員濾波的擴展目標跟蹤方法

2.1 集員濾波相關基礎理論

對UBB噪聲條件下的擴展目標進行跟蹤,其主要思路是利用基于區間數學理論的集員估計方法[20],下面首先介紹集合相關定義及其運算性質。

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

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

(3)

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

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

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

(4)

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

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

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

(5)

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

圖1 外定界橢球

因此,需要根據定理2和定理3計算最優準則下的外定界橢球E(x,S)。

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

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

(6)

證明:利用支持函數可將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),根據橢球定義:

i=1,2

(10)

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

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

(11)

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

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

(12)

根據定義2,上式可變換為:

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

αTx2+(αTS2α)1/2

(13)

則外定界橢球中心x為:

x=x1+x2

(14)

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

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

(15)

(16)

當(γ-1)(ρ-1)≥1時,式(16)成立,可行標量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關于參數p的函數。根據以下定理計算最優橢球參數p。

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

(19)

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

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

(20)

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

(21)

2.2 基于集員濾波的擴展目標跟蹤方法

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

2.2.1 運動狀態預測

在卡爾曼濾波基礎上,狀態集合一步預測為:

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

(22)

預測協方差誤差陣為:

(23)

根據定理2,結合狀態空間模型,可得橢球集合Φkχk-1的支持函數如下:

(24)

由定理1可得狀態集合預測的支持函數:

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

(25)

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

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

(26)

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

(27)

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

(28)

(29)

2.2.2 擴展形態預測

對于擴展形態Xk的預測,假設其在時間間隔內目標形態大小不發生改變[9],則k時刻擴展形態的預測Xk,k-1為k-1時刻擴展目標狀態Xk-1的更新結果,即:

Xk,k-1=Xk-1

(30)

2.2.3 運動狀態更新

(31)

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

(32)

(33)

根據定理1和定理2,結合式(1),可得狀態更新集合χk的支持函數:

(34)

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

(35)

同樣,根據定理4可得關于可行標量pSk的函數表示為:

(36)

則最小跡橢球參數pSk為:

(37)

2.2.4 擴展形態更新

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

(38)

(39)

因量測橢球集合邊界由目標擴展狀態以及傳感器量測誤差邊界信息構成。因此可將其作為目標擴展形態誤差邊界,即包含k時刻所有量測橢球集合的最小邊界為目標擴展形態最大誤差邊界。由于利用橢球集合Minkowski和無法求得包含所有量測橢球的最小外接橢球,因此利用凸包計算幾何中Graham scan算法計算包含k時刻所有量測橢球集合的最小邊界集合A。

(40)

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

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

(a)

(41)

圖3 E1與E2橢球Minkowski差

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

(42)

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

(43)

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

(44)

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

(45)

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

(46)

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

(47)

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

(48)

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

3 仿真實驗

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

3.1 機動場景仿真實驗

(49)

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

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

圖4 4種算法擴展目標跟蹤結果

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

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

(50)

圖6為4種算法100次Monte-Carlo實驗下的GW距離。從圖中可以看出,相比于MEM-EKF、MEM-SOEKF和RMF算法,本文所提算法具有更小的GW距離,且RMF算法估計性能較差。主要原因在于RMF算法將擴展形狀建模為隨機矩陣,導致迭代更新過程中無法處理橢圓長短半軸與角度之間的不確定性導致對目標擴展狀態估計精度較低。雖然MEM-EKF和MEM-SOEKF比RMF算法對目標擴展狀態估計精度較高,但在UBB噪聲條件下受噪聲適應性的影響其估計精度低于本文所提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 參數影響分析

本文所提算法假設量測噪聲為散射源與誤差量測噪聲集合的Minkowski和,初始參數散射因子λ以及量測誤差噪聲邊界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為不同噪聲影響參數λ下SMF-EOT算法位置RMSE和GW距離對比。從圖中可以看出,影響參數λ值越小,SMF-EOT算法收斂性越好。

(a)位置

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

(a)位置

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

4 結語

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