蔡如華，楊 標(biāo)，吳孫勇,2

（1.桂林電子科技大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院，廣西 桂林 541004；2.廣西密碼學(xué)與信息安全重點實驗室，廣西 桂林 541004）

0 引言

多目標(biāo)跟蹤（multi-target tracking,MTT）算法的有效性和實時性一直以來是國內(nèi)外學(xué)者關(guān)注的熱點。針對多目標(biāo)跟蹤問題，Mahler 提出了隨機(jī)有限集[1]（random finite set,RFS）理論，將多目標(biāo)狀態(tài)集和多目標(biāo)觀測集建模為隨機(jī)有限集的形式，有效地解決了基于多假設(shè)跟蹤[2]和聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)[3]在多目標(biāo)跟蹤過程中產(chǎn)生的組合爆炸問題，并提出了基于隨機(jī)有限集的概率假設(shè)密度[4]（probability hypothesis density filter,PHD）濾波器，PHD濾波器能夠?qū)δ繕?biāo)狀態(tài)和數(shù)目進(jìn)行有效的估計。

傳統(tǒng)的粒子濾波器是以點粒子的形式來模擬目標(biāo)狀態(tài)。在進(jìn)行蒙特卡洛實驗時，需要大量的點粒子來對目標(biāo)進(jìn)行捕獲和持續(xù)跟蹤，所以會產(chǎn)生較大的計算量，為減輕計算量，文獻(xiàn)[5]提出了基于區(qū)間分析[6-7]的箱粒子濾波算法，之后在文獻(xiàn)[8]中結(jié)合PHD濾波器實現(xiàn)了箱粒子PHD濾波器，從而減輕了計算量。2015年國內(nèi)學(xué)者宋驪平[9]將箱粒子PHD濾波用于擴(kuò)展目標(biāo)跟蹤，之后在文獻(xiàn)[10]中使用箱粒子PHD濾波對群目標(biāo)進(jìn)行了有效的跟蹤。文獻(xiàn)[11]提出了未知雜波狀態(tài)下基于箱粒子濾波的PHD算法。

但對于低可觀測目標(biāo)的跟蹤，單傳感器PHD濾波器已不能滿足需求，而多傳感器技術(shù)的應(yīng)用則能夠很好地解決這一問題。早在2009年，Mahler[12-13]便給出了多傳感器的實現(xiàn)。文獻(xiàn)[14-15]分別給出了多傳感器PHD濾波的序列蒙特卡洛實現(xiàn)和高斯實現(xiàn)，文獻(xiàn)[16]利用先驗信息對下一時刻可能存在目標(biāo)的區(qū)域進(jìn)行壓縮，并利用多個被動式傳感器同時對可能存在目標(biāo)的區(qū)域進(jìn)行觀測，以多個被動式傳感器的觀測角度相交區(qū)域以及傳感器位置計算得到目標(biāo)位置區(qū)間，作為目標(biāo)區(qū)間量測，根據(jù)區(qū)間量測值產(chǎn)生箱粒子，再利用箱粒子PHD濾波對目標(biāo)進(jìn)行跟蹤，旨在去除虛假量測降低運算量，提升計算效率。

本文則是針對低可觀測目標(biāo)難以跟蹤以及點粒子濾波產(chǎn)生較大計算量的問題，利用多傳感器技術(shù)在低可觀測目標(biāo)檢測上的優(yōu)勢以及箱粒子濾波能夠提升計算效率的優(yōu)勢，提出MS-BOX-PHD算法。并通過仿真實驗驗證了MS-BOX-PHD算法對多目標(biāo)狀態(tài)以及數(shù)目估計的穩(wěn)定有效性，以及MS-BOX-PHD算法相對于區(qū)間量測下多傳感器標(biāo)準(zhǔn)PHD 粒子濾波（multi-sensor standard probability hypothesis density particle filter with interval measurement,IM-PHD-PF）算法在計算效率上的優(yōu)勢。

1 系統(tǒng)模型

在多目標(biāo)跟蹤系統(tǒng)中，會有多個目標(biāo)隨機(jī)的出現(xiàn)和消失在觀測區(qū)域中。全文單目標(biāo)狀態(tài)用小寫字母x表示，其中(x,y)包含目標(biāo)位置信息，包含目標(biāo)速度信息；多目標(biāo)狀態(tài)用大寫字母X表示；單目標(biāo)量測用小寫字母z，z= [θ;r]，其中θ表示角度，r表示距離；多目標(biāo)量測集用大寫字母Z表示。

1.1 多目標(biāo)狀態(tài)轉(zhuǎn)移模型

考慮系統(tǒng)：

式中：k為時間指標(biāo)；f為非線性的狀態(tài)轉(zhuǎn)移函數(shù)；wk表示k時刻的過程噪聲。對于k時刻多目標(biāo)狀態(tài)可以描述為隨機(jī)有限集：

式中：N表示k時刻的目標(biāo)個數(shù)。xk,i表示k時刻的第i個目標(biāo)，可能為新生目標(biāo)，也可能為k－1時刻存活并轉(zhuǎn)移后的目標(biāo)。

1.2 多目標(biāo)多傳感器量測模型

對于k時刻的目標(biāo)xk,i，其可能被傳感器檢測到并產(chǎn)生量測值zk,i，也可能未被傳感器檢測到，沒有產(chǎn)生量測值，同樣在k時刻傳感器也可能產(chǎn)生虛假量測。

由第s∈{1,2,…,S}個傳感器產(chǎn)生的量測集合Zks可表示為：

本文采用非線性量測模型，假設(shè)第s個傳感器位置為(dxs,dys)，則量測方程為：

2 收縮算法

收縮算法以區(qū)間分析[6-7]為基礎(chǔ)，通過區(qū)間收縮，保證更新之后的箱粒子擁有一個合適的區(qū)間大小，并利用收縮前后的箱粒子區(qū)間大小來求解似然函數(shù)，從而對箱粒子的權(quán)重進(jìn)行更新。

假設(shè)原箱粒子為[x]，區(qū)間量測值記為[z]，h為量測方程，若則[x′]＝[x]，但是似然度否則收縮后的箱粒子滿足：

收縮之后得到的箱粒子[x′]的似然度可以表示為：

式中：|[?](i)|表示箱粒子[?]的第i維的長度。

3 MS-BOX-PHD濾波

假設(shè)在k－1時刻的后驗PHD由一組服從均勻分布的箱粒子集描述為：

式中：[xk－1,i]為箱粒子支撐集；Lk－1為k－1時刻箱粒子的個數(shù)。

3.1 MS-BOX-PHD濾波預(yù)測

若γk(x)表示新生目標(biāo)貢獻(xiàn)強(qiáng)度，βSur,k(xk)表示k－1時刻到k時刻存活下來的目標(biāo)貢獻(xiàn)強(qiáng)度，表示k－1時刻目標(biāo)x空間分布概率密度，表示目標(biāo)x由k－1時刻到k時刻的轉(zhuǎn)移密度，則k時刻預(yù)測PHD可以由k時刻的預(yù)測箱粒子參數(shù)集表示為：

新生箱粒子根據(jù)k－1時刻量測產(chǎn)生[5]，在量測上服從均勻分布，為新生箱粒子集合，新生目標(biāo)貢獻(xiàn)強(qiáng)度γk(x)可由其表示為：

取[wk－1]為k－1 時刻的過程噪聲區(qū)間，[fk－1]為k－1 時刻的狀態(tài)轉(zhuǎn)移函數(shù)的包含函數(shù)；pSur([?])表示箱粒子[?]的存活概率，則可由狀態(tài)轉(zhuǎn)移方程得到： 區(qū)1

3.2 MS-BOX-PHD濾波更新

3.2.1 量測融合

量測融合分兩步進(jìn)行：

Step 1：輸入k時刻S個傳感器產(chǎn)生的量測集將S個傳感器產(chǎn)生的區(qū)間量測轉(zhuǎn)換到同一坐標(biāo)系下，步驟如下：

輸入k時刻區(qū)間量測和傳感器坐標(biāo)：

對S個傳感器循環(huán)：

對k時刻第s個傳感器產(chǎn)生所有量測循環(huán)：

Step 2：對S個傳感器產(chǎn)生的量測集進(jìn)行融合。若各傳感器區(qū)間量測間有交集，則保留相交的區(qū)間量測的交集以及其他互不相交的區(qū)間量測值作為當(dāng)前時刻區(qū)間量測集。

以兩個傳感器為例，對k時刻S個傳感器產(chǎn)生的區(qū)間量測進(jìn)行分析，共分為3種情況（如圖1(a)、1(b)、1(c)）。

圖1 傳感器量測分析Fig.1 Sensor measurements analysis

假設(shè)k時刻傳感器1 產(chǎn)生量測集經(jīng)過坐標(biāo)轉(zhuǎn)換后記為傳感器2 產(chǎn)生量測集經(jīng)過坐標(biāo)轉(zhuǎn)換后為記為

顯然因為量測融合同時考慮到多個傳感器的區(qū)間量測，包含目標(biāo)真實信息的可能性也會更大。若其中一個傳感器無法對目標(biāo)進(jìn)行有效檢測時，其他傳感器也能夠?qū)ζ錂z測并產(chǎn)生量測。若多個傳感器同時檢測到目標(biāo)，則通過情形3可知，通過區(qū)間量測相交運算則能夠?qū)δ繕?biāo)量測進(jìn)行壓縮，保證較小的量測區(qū)間即可包含目標(biāo)信息。

通過量測融合后得到的k時刻量測融合集Zk仍然是一個隨機(jī)有限集，包括真實目標(biāo)產(chǎn)生量測和雜波。將其看作是由一個虛擬傳感器產(chǎn)生的量測集，若第s個傳感器的雜波空間分布率和雜波空間密度分別用λs和pdfs表示，則虛擬傳感器的量測集中雜波概率假設(shè)密度為：

若第s個傳感器的檢測概率表示為pDs，則虛擬傳感器的檢測概率為：

3.2.2 濾波更新

輸入k時刻S個傳感器融合后的區(qū)間量測集Zk,對預(yù)測得到的箱粒子集進(jìn)行更新。

4 數(shù)值實驗

實驗計算機(jī)系統(tǒng)為Windows 10 專業(yè)版，64位操作系統(tǒng)。處理器為Intel(R) Core(TM) i5-6500 CPU @3.20 GHz，運行內(nèi)存為8.00 GB。

實驗中設(shè)置所有目標(biāo)在觀測區(qū)域中服從勻速直線（constant velocity,CV）運動模型：

式中：過程噪聲wk為服從協(xié)方差為Q的高斯分布：

本文綜合考慮位置估計誤差和勢估計誤差，使用最優(yōu)子模式分配距離[17]（optimal subpattern assignment,OSPA）來對MS-BOX-PHD濾波器的跟蹤性能進(jìn)行評估。對于任意的兩個多目標(biāo)集合X和Y，若X中包含 |m X=個狀態(tài)，Y中包含 |n Y=個狀態(tài)，則X和Y之間的OSPA 距離定義為：

1）當(dāng)m≤n時：

2）當(dāng)m≥n時：

3）當(dāng)m＝n＝0時：

式中：c為目標(biāo)個數(shù)誤差懲罰參數(shù)；p為階數(shù)，用來懲罰目標(biāo)位置估計誤差。本文實驗中c＝100，p＝1。

4.1 場景1

設(shè)置傳感器1位于坐標(biāo)原點[0,0]的觀測角度范圍為θ1∈[0,π/2]，距離范圍為r1∈[0 m,1500 m]，雜波率為λ1＝2；傳感器2位于[400,0]，觀測角度范圍為θ2∈[π/2,π]，距離范圍為r2∈[0 m,1500 m]，雜波率為λ2＝2；傳感器3位于[700,600]，觀測角度范圍為θ2∈[－π/2,－π]，距離范圍為r3∈[0 m,1500 m]，雜波率為λ3＝2。場景1 中共有10個目標(biāo)隨機(jī)出現(xiàn)和消失在所有傳感器的觀測區(qū)域內(nèi)。試驗中所有目標(biāo)存活并進(jìn)行狀態(tài)轉(zhuǎn)移的概率為pSur＝0.98。

圖2為單次蒙特卡洛實驗下，傳感器檢測概率為0.65時的目標(biāo)真實位置以及傳感器產(chǎn)生的區(qū)間量測。表1為10個目標(biāo)的信息。包括目標(biāo)的初始位置信息、速度信息、出生時刻和死亡時刻。

圖2 目標(biāo)真實位置和區(qū)間量測Fig.2 True targets location and interval measurements

表1 目標(biāo)信息Table1 Information of targets

圖3為單次蒙特卡洛實驗下，MS-BOX-PHD濾波器在傳感器檢測概率都為0.65時對目標(biāo)的跟蹤效果圖，由圖3可知所提MS-BOX-PHD濾波器在檢測概率較低時也能夠有效地對多目標(biāo)狀態(tài)進(jìn)行跟蹤。

圖4為100次蒙特卡洛實驗下，MS-BOX-PHD濾波器和Single-BOX-PHD濾波器在不同檢測概率下對目標(biāo)的勢估計圖，由圖 4可以看出 Single- BOX-PHD濾波在檢測概率逐漸減小時，對目標(biāo)個數(shù)的估計也隨之減少。相對于Single-BOX-PHD濾波器，本文所提MS-BOX-PHD濾波器雖然對目標(biāo)的個數(shù)估計也會受到檢測概率的影響，但受影響幅度明顯較小，在檢測概率為0.95、0.85、0.75、0.65時都能有效地對多目標(biāo)的個數(shù)進(jìn)行估計。

圖3 在x軸和y軸目標(biāo)狀態(tài)估計Fig.3 Targets state estimation in x-axis and y-axis

圖5為100次蒙特卡洛實驗下，MS-BOX-PHD濾波和Single-BOX-PHD濾波在不同檢測概率下對目標(biāo)狀態(tài)估計的OSPA 距離。

圖5(b)為位置OSPA 距離，由圖5(b)可以看出對于位置估計 MS-BOX-PHD濾波器與 Single- BOX-PHD濾波器相差不大，Single-BOX-PHD濾波器位置OSPA 距離普遍稍小是因為隨著檢測概率的減小，目標(biāo)個數(shù)的估計會越來越少，當(dāng)少估計目標(biāo)個數(shù)時候，OSPA 距離只會統(tǒng)計估計到的目標(biāo)在位置上的誤差，所以隨著檢測概率的減小，在位置上的OSPA距離卻越來越小。對于MS-BOX-PHD濾波器，能夠在檢測概率較小時對目標(biāo)個數(shù)進(jìn)行正確的估計，但在位置估計上會統(tǒng)計所有估計到的目標(biāo)的位置OSPA距離，所以普遍稍大。且檢測概率越大，對位置的估計會越精確；圖5(c)為勢估計OSPA 距離，由圖5(c)可以看出在對目標(biāo)個數(shù)估計上Single-BOX-PHD濾波器明顯會受到檢測概率的影響，檢測概率越低，對目標(biāo)個數(shù)的估計誤差也越大；圖5(a)為綜合考慮目標(biāo)位置估計和個數(shù)估計的OSPA 距離，由圖5(a)可以看出MS-BOX-PHD濾波器對多目標(biāo)估計的性能會受到檢測概率的影響，但所受影響相對于Single-BOX-PHD濾波器受檢測概率的影響要小得多，且只有在目標(biāo)新生或者死亡時才會出現(xiàn)較大的估計誤差，如第1、13、15、24、28、41、45時刻都有目標(biāo)新生。說明MS-BOX-PHD濾波器在檢測概率較低時，能夠有效地估計多目標(biāo)的狀態(tài)和個數(shù)。

圖4 不同檢測概率下勢估計（100MC）Fig.4 Cardinality estimation under different detection probabilities (100MC)

4.2 場景2

設(shè)置傳感器1和傳感器2的觀測角度范圍都為θ∈[0 rad,π rad]，距離觀測范圍都為r∈[37.5 m,1000 m]，雜波率為λ1＝λ2＝2；場景2 中共有2個目標(biāo)，兩個目標(biāo)的航跡都在傳感器1和傳感器2的觀測區(qū)域內(nèi)。目標(biāo)1的初始狀態(tài)為[－800 m;16 m/s;100 m;7 m/s]，存活時間為1 s→100 s，目標(biāo)2的初始狀態(tài)為[800 m;－16 m/s;100 m;7 m/s]，存活時間為30 s→90 s。試驗中所有目標(biāo)存活并進(jìn)行狀態(tài)轉(zhuǎn)移的概率為pSur＝0.98，每個目標(biāo)被傳感器檢測到的概率都是pD＝0.65。

圖6為100次蒙特卡洛實驗下，MS-BOX-PHD濾波器和粒子數(shù)分別為3000個、4000個、4500個、5000個、8000個時，IM-PHD-PF 濾波器對目標(biāo)勢估計對比圖。

由圖6可以看出IM-PHD-PF 濾波器隨著粒子數(shù)的增多，對勢的估計也會越來越準(zhǔn)確。但在粒子數(shù)為3000個、4000個、4500個、5000個時都沒有箱粒子數(shù)為8個時MS-BOX-PHD濾波器對目標(biāo)勢估計準(zhǔn)確，而在IM-PHD-PF 濾波器粒子數(shù)提高到8000個時，IM-PHD-PF 濾波器對目標(biāo)的估計個數(shù)已經(jīng)接近于真實目標(biāo)個數(shù)。

圖7為100次蒙特卡洛實驗下，MS-BOX-PHD濾波器和粒子數(shù)分別為3000個、4000個、4500個、5000個、8000個時，IM-PHD-PF 濾波器對多目標(biāo)狀態(tài)估計的OSPA 距離。

圖5 不同檢測概率下OSPA 距離估計（100MC）Fig.5 OSPA distance estimation under different detection probabilities (100MC)

圖6 MS-BOX-PHD濾波器和不同粒子數(shù)下IM-PHD-PF 濾波器對目標(biāo)勢估計（100MC）Fig.6 Cardinality estimation by utilizing MS-BOX-PHD filter and IM-PHD-PF filter with different number of particles,respectively (100MC)

圖7 MS-BOX-PHD濾波器和不同粒子數(shù)下IM-PHD-PF 濾波器對目標(biāo)狀態(tài)估計的OSPA 距離（100MC）Fig.7 OSPA distance estimation by utilizing MS-BOX-PHD filter and IM-PHD-PF filter with different number of particles,respectively (100MC)

其中圖7(c)為勢估計OSPA 距離，圖7(b)為位置估計OSPA 距離，圖7(b)中，IM-PHD-PF 濾波器在位置估計OSPA 距離要比MS-BOX-PHD濾波器位置估計OSPA 距離小，因為箱粒子是一個區(qū)間，最終是使用箱粒子中心位置來對目標(biāo)位置進(jìn)行估計的，所以偏差要大于點粒子下IM-PHD-PF 濾波器對目標(biāo)位置的估計；圖7(c)則說明，粒子數(shù)越少IM-PHD-PF濾波器在對目標(biāo)進(jìn)行跟蹤時，越容易丟失目標(biāo)；圖7(a)綜合考慮目標(biāo)的位置估計誤差和勢估計誤差，圖7(a)表明MS-BOX-PHD濾波器僅使用8個箱粒子對多目標(biāo)狀態(tài)估計的性能，要比IM-PHD-PF 濾波器使用3000、4000、4500、5000個粒子的情況下要穩(wěn)定得多。

圖7(a)中MS-BOX-PHD濾波器的OSPA 距離均值為28.4835，運行時間為14.97 s。

表2為100次蒙特卡洛實驗下IM-PHD-PF 濾波器使用不同粒子數(shù)下的運行時間以及OSPA 距離均值，由表2可以看出隨著粒子數(shù)目的增多，區(qū)間量測下多傳感器標(biāo)準(zhǔn)PHD 粒子濾波器對多目標(biāo)狀態(tài)估計性能越來越好，但運行時間花費也越來越大。

綜合表2以及圖6和圖7可知，在IM-PHD-PF濾波器使用4500個粒子時運行時間為24.37 s，OSPA 距離均值為30.14；MS-BOX-PHD濾波器使用8個箱粒子，運行時間為14.97 s，OSPA距離均值為28.4835。說明IM-PHD-PF濾波器和MS-BOX-PHD濾波器在對多目標(biāo)狀態(tài)估計性能相近時候，MS-BOX-PHD濾波器要比IM-PHD-PF 濾波器在計算效率上提升了38.57%。

表2 IM-PHD-PF 濾波器運行時間（100MC）Table2 Time cost for IM-PHD-PF filter

5 結(jié)論

本文利用多傳感器技術(shù)，結(jié)合BOX PF 以及PHD濾波器給出了MS-BOX-PHD算法的實現(xiàn)。并通過數(shù)值實驗，驗證了對于低可觀測目標(biāo)MS-BOX-PHD算法的優(yōu)越跟蹤性能，并且驗證了與區(qū)間量測下多傳感器標(biāo)準(zhǔn)PHD 粒子濾波器相比較，當(dāng)跟蹤性能接近時，MS-BOX-PHD濾波器在計算效率上提升了38.57%.