張俊根

(北方民族大學(xué)電氣信息工程學(xué)院，寧夏銀川 750021)

1 引言

在純方位多目標(biāo)跟蹤問(wèn)題中，由于目標(biāo)信息的不確定性，目標(biāo)數(shù)目每個(gè)時(shí)刻都會(huì)發(fā)生變化.此外，目標(biāo)運(yùn)動(dòng)建模與量測(cè)建模坐標(biāo)系轉(zhuǎn)換過(guò)程中以及被動(dòng)傳感器自身物理特性造成的模型非線性問(wèn)題，以及量測(cè)信息的不完備性，均給目標(biāo)跟蹤帶來(lái)了很大困難.如何根據(jù)純方位量測(cè)信息對(duì)多個(gè)目標(biāo)進(jìn)行實(shí)時(shí)、有效地跟蹤，一直是學(xué)術(shù)界和工程應(yīng)用領(lǐng)域的研究熱點(diǎn)和難點(diǎn)問(wèn)題[1-2].

相對(duì)于其他傳統(tǒng)的多目標(biāo)跟蹤算法，基于隨機(jī)集理論的概率假設(shè)密度(probability hypothesis density，PHD)濾波算法，可以將復(fù)雜的多目標(biāo)狀態(tài)空間的運(yùn)算轉(zhuǎn)換為單目標(biāo)狀態(tài)空間內(nèi)的運(yùn)算，有效避免了多目標(biāo)跟蹤中復(fù)雜的數(shù)據(jù)關(guān)聯(lián)組合問(wèn)題[3-6].而勢(shì)概率假設(shè)密度(cardinalized probability hypothesis density，CPHD)濾波算法，能充分利用多目標(biāo)密度的信息，并且不需要限定目標(biāo)數(shù)目服從泊松分布，受到了更廣泛的關(guān)注，許多學(xué)者開展了相關(guān)的研究[7-12].文獻(xiàn)[13]基于著名的Faàdi Bruno行列式，推導(dǎo)了CPHD模型中的衍生目標(biāo)廣義勢(shì)分布預(yù)測(cè)方程，提出了一種通用勢(shì)分布預(yù)測(cè)方程的可處理遞推計(jì)算方法.文獻(xiàn)[14]研究了基于最適高斯近似和強(qiáng)跟蹤的多模型伽馬高斯逆威夏特-概率假設(shè)密度濾波的群跟蹤算法，能有效處理群目標(biāo)的合并和衍生情況.文獻(xiàn)[15]針對(duì)擴(kuò)展目標(biāo)跟蹤問(wèn)題，提出了一種擴(kuò)展的標(biāo)記箱粒子CPHD濾波器，可以改善目標(biāo)數(shù)目估計(jì)的精度，同時(shí)能獲得目標(biāo)的軌跡.文獻(xiàn)[16]針對(duì)多傳感器多目標(biāo)跟蹤中的分布式傳感器控制問(wèn)題，利用隨機(jī)有限集進(jìn)行建模，提出一種基于多目標(biāo)戰(zhàn)術(shù)重要性評(píng)估的多傳感器控制策略.

本文針對(duì)混合線性/非線性目標(biāo)狀態(tài)模型，將CPHD和Rao-Blackwellised的思想[17]有效結(jié)合起來(lái)，提出一種新的混合粒子濾波CPHD實(shí)現(xiàn)算法，將目標(biāo)狀態(tài)劃分為線性和非線性這兩種結(jié)構(gòu)信息，線性狀態(tài)信息利用卡爾曼濾波(Kalman filter，KF)，非線性狀態(tài)信息利用粒子濾波(particle filter，PF)，對(duì)多目標(biāo)的概率假設(shè)密度和勢(shì)分布進(jìn)行預(yù)測(cè)與估計(jì)，再運(yùn)用Mean-Shift方法提取目標(biāo)狀態(tài)估計(jì)值，并對(duì)算法的復(fù)雜度進(jìn)行了推導(dǎo)和分析，最后，將算法應(yīng)用到多傳感器純方位目標(biāo)跟蹤場(chǎng)景中，仿真實(shí)驗(yàn)結(jié)果表明，所提出的算法能夠?qū)崿F(xiàn)對(duì)目標(biāo)數(shù)和目標(biāo)狀態(tài)的準(zhǔn)確估計(jì)，改善目標(biāo)跟蹤精度.

2 問(wèn)題描述

2.1 隨機(jī)集純方位多目標(biāo)跟蹤問(wèn)題

考慮下面的被動(dòng)多傳感器純方位多目標(biāo)跟蹤系統(tǒng):

在隨機(jī)集多目標(biāo)跟蹤問(wèn)題中，多目標(biāo)的狀態(tài)和傳感器量測(cè)可以建模為隨機(jī)有限集的形式[3]，即Xk={xk，1，…，xk，Nk}∈F(X)為目標(biāo)狀態(tài)集，Zk={zk，1，…，zk，Mk}∈F(Z)為量測(cè)集，其中:F(X)和F(Z)分別是X和Z上的所有有限子集的集合;Nk及Mk分別表示k時(shí)刻目標(biāo)個(gè)數(shù)及量測(cè)個(gè)數(shù)，其中某些量測(cè)可能源于雜波.

若k-1時(shí)刻狀態(tài)隨機(jī)集為Xk-1，則k時(shí)刻的狀態(tài)隨機(jī)集Xk可表示為

其中:Γk表示k時(shí)刻新生目標(biāo)狀態(tài)隨機(jī)集;Bk|k-1(x)表示k時(shí)刻從目標(biāo)x衍生分裂出來(lái)的目標(biāo)狀態(tài)隨機(jī)集;Sk|k-1(x)表示從k-1時(shí)刻到k時(shí)刻仍然存活的目標(biāo)狀態(tài)隨機(jī)集.

目標(biāo)的量測(cè)隨機(jī)集Zk可表示為

其中:Kk表示虛警或雜波量測(cè)隨機(jī)集;Θk(x)表示源于真實(shí)目標(biāo)的量測(cè)隨機(jī)集.

2.2 CPHD濾波

相比于PHD濾波，CPHD濾波聯(lián)合估計(jì)目標(biāo)強(qiáng)度和勢(shì)分布(即目標(biāo)數(shù)n的概率分布)，可以改善狀態(tài)估計(jì)性能[7].

令Dk|k-1和pk|k-1分別表示k-1時(shí)刻預(yù)測(cè)多目標(biāo)強(qiáng)度和勢(shì)分布，Dk和pk表示k時(shí)刻后驗(yàn)多目標(biāo)強(qiáng)度和勢(shì)分布，則CPHD迭代公式為

其中:N(x;m，P)表示均值為m、協(xié)方差為P的高斯密度，Qk-1是過(guò)程噪聲協(xié)方差，Rk是量測(cè)噪聲協(xié)方差.

3 多目標(biāo)跟蹤粒子CPHD濾波

在很多應(yīng)用中，目標(biāo)狀態(tài)空間包含線性和非線性兩部分，針對(duì)這類混合線性/非線性狀態(tài)模型，利用貝葉斯原理將線性狀態(tài)從系統(tǒng)模型中分離出來(lái)，并結(jié)合卡爾曼濾波器和粒子濾波器對(duì)目標(biāo)的線性狀態(tài)和非線性狀態(tài)分別進(jìn)行估計(jì)[17]，本文結(jié)合RBPF的思想，估計(jì)出多目標(biāo)PHD和勢(shì)分布，由于多目標(biāo)跟蹤的復(fù)雜性，使得迭代估計(jì)所得的PHD形狀很不規(guī)則，運(yùn)用Mean-Shift方法，可以提取出精確的密度函數(shù)峰值位置[18]，可以提高目標(biāo)跟蹤精度，減小估計(jì)方差.

3.1 算法描述

假定目標(biāo)不發(fā)生衍生，各個(gè)目標(biāo)的動(dòng)態(tài)系統(tǒng)同時(shí)存在線性和非線性的狀態(tài)，且噪聲是加性的，那么目標(biāo)濾波模型可以表示為線性和非線性的形式:

且wk與vk相互獨(dú)立.

算法由預(yù)測(cè)和更新兩部分構(gòu)成:

預(yù)測(cè):假定已知k-1時(shí)刻的多目標(biāo)后驗(yàn)概率假設(shè)密度Dk-1和勢(shì)分布pk-1，并且Dk-1可以用一組帶權(quán)值的粒子來(lái)表達(dá):

為了得到Dk|k-1的粒子近似，對(duì)式(22)中的每一項(xiàng)都運(yùn)用Rao-Blackwellised方法進(jìn)行混合濾波.首先，對(duì)于存活目標(biāo)，由式(19)預(yù)測(cè)非線性狀態(tài)粒子

線性狀態(tài)粒子可由KF方程獲得

其中

計(jì)算粒子的權(quán)值

勢(shì)分布pk|k-1(n)通過(guò)式(5)來(lái)計(jì)算.

更新:已知k-1時(shí)刻的多目標(biāo)預(yù)測(cè)概率假設(shè)密度Dk|k-1和勢(shì)分布pk|k-1，并且Dk|k-1可由粒子來(lái)描述.運(yùn)用KDE理論[19]可得多目標(biāo)隨機(jī)集的更新PHD:

勢(shì)分布pk(n)利用式(7)計(jì)算.

目標(biāo)數(shù)的估計(jì)為

利用Mean-Shift算法，提取密度函數(shù)所有的峰值位置.對(duì)于粒子其Mean-Shift向量為

本文所提算法:最密集的位置，即

1) 給定Dk-1和pk-1

預(yù)測(cè):

2) 對(duì)于i=1，…，Lk-1，利用式(23)計(jì)算

7) 根據(jù)式(5)預(yù)測(cè)勢(shì)分布pk|k-1.

更新:

9) 根據(jù)式(28)估計(jì)目標(biāo)強(qiáng)度函數(shù)Dk;

10) 利用式(7)更新勢(shì)分布pk，并利用式(31)估計(jì)目標(biāo)數(shù)nk;

3.2 算法復(fù)雜度分析

這一小節(jié)引入等效浮點(diǎn)運(yùn)算(equivalent flop，EF)測(cè)度來(lái)分析算法的復(fù)雜度[21]，一個(gè)運(yùn)算的EF復(fù)雜度定義為浮點(diǎn)運(yùn)算(flops)的數(shù)量，這里一次浮點(diǎn)運(yùn)算表示兩個(gè)浮點(diǎn)數(shù)的一個(gè)加、減、乘或除.從算法描述看出，算法的計(jì)算量主要體現(xiàn)在更新步驟，其中包括漏檢目標(biāo)和檢測(cè)目標(biāo)的粒子似然比的計(jì)算[8]、粒子Mean-Shift向量的計(jì)算以及重采樣的計(jì)算，對(duì)于每次時(shí)間迭代，算法更新步驟總的加法或乘法次數(shù)可以表示為量測(cè)數(shù)Mk的一個(gè)多項(xiàng)式

其中:Lk表示k時(shí)刻的粒子數(shù)，Nmax為最大目標(biāo)數(shù)，系數(shù)c1，c2，c3分別用來(lái)等效高斯似然、Mean-Shift向量以及重采樣的計(jì)算復(fù)雜度.

從式(33)可以看出，在給定最大目標(biāo)數(shù)Nmax的情況下:1)隨著量測(cè)數(shù)的增大，算法復(fù)雜度呈三次方增長(zhǎng);2)隨著粒子數(shù)的增大，算法復(fù)雜度呈二次方增長(zhǎng).利用門選技術(shù)，刪除一些不相關(guān)的量測(cè)，可以減小算法的計(jì)算量.

4 仿真實(shí)驗(yàn)與分析

對(duì)所提算法進(jìn)行性能仿真實(shí)驗(yàn)，利用目標(biāo)數(shù)估計(jì)的統(tǒng)計(jì)值和OSPA距離評(píng)價(jià)算法的性能，并與PFCPHD，RBPF-PHD這兩種濾波算法作對(duì)比.

系統(tǒng)模型如式(1)-(2)，其中目標(biāo)數(shù)隨時(shí)間變化，假設(shè)有3個(gè)被動(dòng)傳感器對(duì)目標(biāo)進(jìn)行純方位跟蹤，假定各傳感器的量測(cè)數(shù)據(jù)已完成配準(zhǔn)及關(guān)聯(lián).系統(tǒng)具體參數(shù)設(shè)置如下:

采樣周期Ts=1 s;目標(biāo)系統(tǒng)矩陣

過(guò)程噪聲Q=diag{[0.01 0.01]};3個(gè)傳感器位置分別為(-8，-10)km，(8，-10)km和(8，13.86)km.

被動(dòng)傳感器對(duì)目標(biāo)只能測(cè)得方位角，量測(cè)函數(shù)

其中:(xo，yo)是傳感器位置，量測(cè)噪聲標(biāo)準(zhǔn)差都設(shè)為σβ=0.01，目標(biāo)存活概率Ps=0.99，檢測(cè)概率Pd=0.98.

新生目標(biāo)隨機(jī)集的強(qiáng)度函數(shù)為

其中:

雜波均勻分布于量測(cè)空間，數(shù)目服從參數(shù)為λ=3的泊松分布;OSPA距離參數(shù)p=2，c=100;粒子數(shù)500，仿真次數(shù)為100.

圖1為多目標(biāo)在x和y方向上的運(yùn)動(dòng)軌跡，圖2為各算法目標(biāo)數(shù)估計(jì)均值，圖3給出了各算法目標(biāo)數(shù)估計(jì)標(biāo)準(zhǔn)差.可以看到，PF-CPHD，RBPF-PHD和本文算法的目標(biāo)數(shù)估計(jì)都能收斂到真實(shí)值，RBPF-PHD 濾波算法對(duì)目標(biāo)數(shù)的估計(jì)能較快收斂;本文算法對(duì)目標(biāo)數(shù)估計(jì)的方差最小，表明所提算法的目標(biāo)數(shù)估計(jì)要更加可靠.另外，從圖2可知，當(dāng)目標(biāo)數(shù)發(fā)生變化時(shí)，PFCPHD和本文算法對(duì)目標(biāo)數(shù)的估計(jì)具有較大的延時(shí)，而RBPF-PHD濾波算法的目標(biāo)數(shù)估計(jì)延時(shí)較小，原因主要在于RBPF-PHD濾波算法的目標(biāo)數(shù)估計(jì)方差較大，估計(jì)的結(jié)果不可靠，導(dǎo)致濾波器更容易受新來(lái)量測(cè)的影響，響應(yīng)較快.

圖1 實(shí)際的目標(biāo)運(yùn)動(dòng)軌跡Fig.1 The true target tracks in x and y-coordinates

圖2 目標(biāo)數(shù)估計(jì)均值Fig.2 The estimated mean of cardinality

圖3 目標(biāo)數(shù)估計(jì)標(biāo)準(zhǔn)差Fig.3 The estimated STD of cardinality

圖4給出了3種算法目標(biāo)跟蹤的OSPA距離.可以看出，在多數(shù)時(shí)刻，本文算法的跟蹤誤差比PF-CPHD和RBPF-PHD濾波算法都要小;在目標(biāo)數(shù)發(fā)生變化的時(shí)候，RBPF-PHD算法的跟蹤精度要好于另外兩種算法，主要原因是PF-CPHD和本文算法在目標(biāo)數(shù)發(fā)生變化時(shí)，對(duì)目標(biāo)數(shù)的估計(jì)存在較大的延時(shí)，OSPA距離懲罰的程度比RBPF-PHD濾波算法要大;而在目標(biāo)數(shù)穩(wěn)定以后，本文算法的OSPA距離相對(duì)較小.

圖4 目標(biāo)跟蹤的OSPA距離Fig.4 100 MC average OSPA distance

圖5給出了這3種算法的平均OSPA距離隨粒子數(shù)變化的對(duì)比.可以看出，本文算法的平均OSPA距離比PF-CPHD 和RBPF-PHD 算法都要小;當(dāng)粒子數(shù)由300增大到500時(shí)，RBPF-PHD濾波算法的平均OSPA距離有較大的改善，當(dāng)粒子數(shù)繼續(xù)增大時(shí)，性能改善不明顯;PF-CPHD和本文算法隨著粒子數(shù)的增大，平均OSPA距離都將減小.當(dāng)然，隨著粒子數(shù)的增多，各算法的運(yùn)算時(shí)間也將增大.圖6給出了3種算法平均運(yùn)算時(shí)間對(duì)比，本文算法由于要聯(lián)合估計(jì)目標(biāo)強(qiáng)度及目標(biāo)數(shù)的概率分布，并且運(yùn)用了Mean-Shift算法進(jìn)行狀態(tài)提取，運(yùn)算量較大.

圖5 各算法的平均OSPA距離對(duì)比Fig.5 Comparison of the mean of MC average OSPA distance

圖6 各算法的平均運(yùn)算時(shí)間對(duì)比Fig.6 Comparison of average operation time

為了驗(yàn)證本文算法的復(fù)雜度與量測(cè)數(shù)的關(guān)系，對(duì)仿真實(shí)驗(yàn)重新設(shè)計(jì)，假定目標(biāo)數(shù)不隨時(shí)間變化，通過(guò)增大目標(biāo)數(shù)來(lái)增大量測(cè)數(shù)，雜波強(qiáng)度設(shè)為3，粒子數(shù)為500，其他仿真參數(shù)不變，算法復(fù)雜度用運(yùn)算時(shí)間來(lái)代表.圖7給出了本文算法的平均運(yùn)算時(shí)間與量測(cè)數(shù)的關(guān)系，為便于對(duì)比分析，圖中畫出了一個(gè)關(guān)于量測(cè)數(shù)的三次方曲線方程.可以看出，本文算法復(fù)雜度由于與多個(gè)參量有關(guān)系，與量測(cè)數(shù)近似滿足隨著量測(cè)數(shù)的增大復(fù)雜度呈三次方增長(zhǎng)的關(guān)系.

圖7 本文算法平均運(yùn)算時(shí)間與量測(cè)數(shù)的關(guān)系Fig.7 The relationship between the average operation time and the number of measurements

5 結(jié)語(yǔ)

本文針對(duì)混合線性/非線性目標(biāo)狀態(tài)空間跟蹤模型，提出了一種多傳感器純方位多目標(biāo)跟蹤粒子CPHD濾波算法，通過(guò)挖掘混合線性/非線性狀態(tài)模型的結(jié)構(gòu)信息，結(jié)合PF與KF對(duì)各個(gè)目標(biāo)的狀態(tài)進(jìn)行預(yù)測(cè)與估計(jì)，以便更好地估計(jì)目標(biāo)強(qiáng)度函數(shù)及勢(shì)分布，運(yùn)用Mean-Shift算法提取多目標(biāo)密度函數(shù)的峰值作為目標(biāo)狀態(tài)估計(jì)值，并分析了算法的復(fù)雜度，仿真結(jié)果驗(yàn)證了算法的有效性.算法雖然提升了跟蹤精度，但帶來(lái)的計(jì)算復(fù)雜度較高，如何利用被動(dòng)傳感器的分布式計(jì)算來(lái)減小算法的運(yùn)行時(shí)間，將是后續(xù)研究的方向.