parfor 模式在粒子概率假設(shè)密度濾波中的應(yīng)用

陳金廣，馬全海

(西安工程大學(xué)計(jì)算機(jī)科學(xué)學(xué)院，西安 710048)

針對(duì)多目標(biāo)跟蹤，Mahler 在貝葉斯框架的基礎(chǔ)上提出了概率假設(shè)密度濾波(probability hypothesis density filter，PHD)，通過(guò)遞推多目標(biāo)后驗(yàn)的一階統(tǒng)計(jì)矩，從一階統(tǒng)計(jì)矩中提取目標(biāo)的狀態(tài)［1］。目前，PHD 濾波主要有2 種實(shí)現(xiàn)方式:一種是基于粒子概率假設(shè)密度濾波(Particle PHD，P-PHD)(也叫序貫蒙特卡羅(Sequential Monte Carlo，SMC)的實(shí)現(xiàn)方式)［2］;另外一種方式是高斯混合概率假設(shè)密度濾波器(Gaussian Mixture PHD，GM-PHD)［3］。這2 種方式雖然都是由Vo 提出的，但是適用的情況卻不同。前者用諸多的粒子來(lái)近似PHD 濾波的一階統(tǒng)計(jì)矩，能適用于非線性非高斯系統(tǒng);而后者僅適用于線性高斯系統(tǒng)。

基于粒子的PHD 濾波雖然有其自身的優(yōu)點(diǎn)，但是隨著粒子數(shù)目的增多會(huì)使得計(jì)算量相應(yīng)增大，從而導(dǎo)致時(shí)間復(fù)雜度也隨之增加。文獻(xiàn)［4］中提出一種非序貫重采樣算法，能夠?qū)崿F(xiàn)P-PHD 濾波算法中權(quán)值更新步驟和重釆樣步驟的并行操作，在同樣的濾波效果下，降低了時(shí)間復(fù)雜度。文獻(xiàn)［5］中根據(jù)相對(duì)熵來(lái)動(dòng)態(tài)計(jì)算每個(gè)目標(biāo)所分配的粒子數(shù)目，而不是固定的給每個(gè)目標(biāo)分配固定的粒子數(shù)，從而提高了濾波效率，減少了計(jì)算量。文獻(xiàn)［6］中則采用擬蒙特卡羅方法并行處理多目標(biāo)跟蹤的數(shù)據(jù)關(guān)聯(lián)，使得運(yùn)算效率得到提高。文獻(xiàn)［7］中間接地使用并行算法，將目標(biāo)分為已有目標(biāo)和新生目標(biāo)兩類，利用極大似然門限來(lái)分配兩類目標(biāo)對(duì)應(yīng)的最優(yōu)觀測(cè)數(shù)據(jù)，然后分別進(jìn)行目標(biāo)狀態(tài)更新，從而提高了濾波性能。文獻(xiàn)［8］中則是通過(guò)自適應(yīng)門限濾掉雜波，然后僅利用落入門限內(nèi)的量測(cè)進(jìn)行更新，大大降低了時(shí)間復(fù)雜度。文獻(xiàn)［9］中闡述了粒子濾波具有可并行化的特點(diǎn)，并用OpenMP 實(shí)現(xiàn)了粒子濾波目標(biāo)跟蹤的并行化處理，提高了粒子濾波的計(jì)算效率。文獻(xiàn)［10］中提出了獨(dú)立性分布式粒子濾波并行計(jì)算。并驗(yàn)證了在處理時(shí)間上比粒子濾波更有優(yōu)勢(shì)。文獻(xiàn)［11］中將粒子濾波和GPU 結(jié)合起來(lái)，充分發(fā)揮GPU(Graphics Processing Unit)的并行計(jì)算功能，在濾波效果相當(dāng)?shù)那闆r下，顯著減少了處理時(shí)延。文獻(xiàn)［12］中把并行算法引入到重采樣中，降低了重采樣的處理時(shí)延。

上述這些算法表明了并行計(jì)算在粒子濾波中具有可行性。文獻(xiàn)［4］中只是在理論上闡述，沒(méi)有實(shí)現(xiàn)粒子濾波的并行處理。文獻(xiàn)［11］中用OpenMP 語(yǔ)言實(shí)現(xiàn)了基于粒子濾波的圖像并行處理跟蹤，減少了處理時(shí)延。本研究把Matlab 中的parfor 并行計(jì)算結(jié)構(gòu)引入到基于粒子的概率假設(shè)密度濾波的更新步驟中，實(shí)現(xiàn)了并行實(shí)驗(yàn)。并行算法和傳統(tǒng)的PPHD 濾波算法相比在濾波精度上基本相同，隨著并行核數(shù)的增加，更新步驟的時(shí)間復(fù)雜度明顯減少。

1 背景知識(shí)

1.1 多目標(biāo)隨機(jī)集模型

在基于隨機(jī)有限集(RFS)多目標(biāo)濾波中，分為多目標(biāo)狀態(tài)集和多目標(biāo)量測(cè)集。k 時(shí)刻的狀態(tài)集和量測(cè)集可表示為:Xk={xk，1，…，xk，N}，Zk={ zk，1，…，zk，M}。其中M，N 對(duì)應(yīng)表示量測(cè)目標(biāo)個(gè)數(shù)和目標(biāo)個(gè)數(shù)。也可表示為

式中:Kk表示所有雜波的量測(cè)集;表示所有真實(shí)目標(biāo)所產(chǎn)生的量測(cè)集。

1.2 基于粒子的概率假設(shè)密度濾波(P-PHDF)算法

對(duì)應(yīng)的，PHD 的粒子實(shí)現(xiàn)也對(duì)應(yīng)為預(yù)測(cè)步驟和更新步驟。預(yù)測(cè)分為對(duì)新生目標(biāo)的預(yù)測(cè)以及對(duì)存活目標(biāo)的預(yù)測(cè)兩部分，需要預(yù)測(cè)的是對(duì)應(yīng)粒子的狀態(tài)和權(quán)重。對(duì)于新生目標(biāo)，其目標(biāo)狀態(tài)和權(quán)值的預(yù)測(cè)方程為

對(duì)于存活目標(biāo)，目標(biāo)狀態(tài)和權(quán)重的預(yù)測(cè)由方程式(4)得到

式中:Lk-1、Jk分別為存活目標(biāo)和新生目標(biāo)對(duì)應(yīng)粒子的數(shù)目;對(duì)應(yīng)表示k -1 時(shí)刻和k 時(shí)刻第i 個(gè)粒子的狀態(tài);為k 時(shí)刻的第i 粒子對(duì)應(yīng)的權(quán)重;p(·|Zk)和q(·|Zk)分別對(duì)應(yīng)于新生目標(biāo)和存活目標(biāo)的重要性函數(shù)為新生目標(biāo)的概率假設(shè)密度函數(shù)。其中:ek|k-1(·|·)表示k-1 時(shí)刻的目標(biāo)在k 時(shí)刻仍然存活的概率，fk|k-1(·|·)表示單個(gè)目標(biāo)的狀態(tài)轉(zhuǎn)移概率密度，βk|k-1(·|·)表示k -1 時(shí)刻目標(biāo)衍生的概率假設(shè)密度。

在基于粒子的PHD 濾波的更新步驟中，只需更新粒子的權(quán)重，粒子的狀態(tài)無(wú)需更新。權(quán)重更新方程為:

2 parfor 在P-PHD 中的應(yīng)用

并行計(jì)算是一種許多指令同時(shí)進(jìn)行的計(jì)算模式。它的思想是在同時(shí)進(jìn)行的前提下，將計(jì)算的過(guò)程分解成小部分，之后以并發(fā)方式來(lái)加以解決。并行計(jì)算將一個(gè)復(fù)雜的任務(wù)分配給多個(gè)處理器去處理，而不是由一個(gè)處理器去順序執(zhí)行，從而能夠減少處理的時(shí)延，大大提高運(yùn)行的效率。Matlab自帶的工具箱可以在多核計(jì)算機(jī)上解決計(jì)算問(wèn)題，其并行處理結(jié)構(gòu)為任務(wù)并行(parfor 循環(huán))和數(shù)據(jù)并行。①使用parfor循環(huán)執(zhí)行任務(wù)并行時(shí)，Matlab 會(huì)將獨(dú)立的循環(huán)迭代自動(dòng)分配給多個(gè)Matlab worker。工具箱中的parfor 循環(huán)提供了一種在多個(gè)Matlab worker 間分配任務(wù)的方式。使用該循環(huán)，可以將獨(dú)立的循環(huán)迭代自動(dòng)分配給多個(gè)Matlab worker。parfor(并行for 循環(huán))結(jié)構(gòu)管理著Matlab 客戶端會(huì)話與worker 之間的數(shù)據(jù)和代碼傳輸。它會(huì)自動(dòng)檢測(cè)是否有worker，如果沒(méi)有，則會(huì)還原為串行方式。②執(zhí)行數(shù)據(jù)并行時(shí)，對(duì)于需要大型數(shù)據(jù)集處理的Matlab 算法，Parallel Computing Toolbox 提供了分布式數(shù)組、并行函數(shù)，以及使用spmd (單程序多數(shù)據(jù))關(guān)鍵詞注釋代碼區(qū)段的功能，可用于在數(shù)個(gè)worker 上并行執(zhí)行。這些并行結(jié)構(gòu)可處理worker 間通信，并協(xié)調(diào)后臺(tái)的并行計(jì)算。

Matlab 的并行計(jì)算實(shí)質(zhì)還是主從結(jié)構(gòu)的分布式計(jì)算。當(dāng)初始化Matlab 并行計(jì)算環(huán)境時(shí)，最初的Matlab 進(jìn)程自動(dòng)成為主節(jié)點(diǎn)，同時(shí)初始化多個(gè)(具體個(gè)數(shù)手動(dòng)設(shè)定)Matlab 計(jì)算子節(jié)點(diǎn)。Parfor 的作用就是讓這些子節(jié)點(diǎn)同時(shí)運(yùn)行Parfor語(yǔ)句段中的代碼。Parfor 運(yùn)行之初，主節(jié)點(diǎn)會(huì)將Parfor 循環(huán)程序之外變量傳遞給計(jì)算子節(jié)點(diǎn)。子節(jié)點(diǎn)運(yùn)算過(guò)程時(shí)互不干擾，運(yùn)算完畢，則應(yīng)該有相應(yīng)代碼將各子節(jié)點(diǎn)得到的結(jié)果組合到同一個(gè)數(shù)組變量中，并返回到Matlab 主節(jié)點(diǎn)。在Matlab 并行工具箱中，已經(jīng)提到了parfor 并行循環(huán)對(duì)于蒙特卡羅效果較好，本研究將parfor 并行循環(huán)引入到基于粒子的概率假設(shè)密度濾波中，利用其并行計(jì)算來(lái)減少運(yùn)行時(shí)延。

在P-PHD 濾波框架中，其主要步驟分為初始化、預(yù)測(cè)、更新、估計(jì)目標(biāo)數(shù)、重采樣和聚類等。在時(shí)間復(fù)雜度方面，重采樣和聚類較大，更新步驟次之，預(yù)測(cè)和其他步驟的時(shí)延基本可以忽略。傳統(tǒng)的P-PHD 濾波中，除更新步驟外，其他步驟暫不具備并行計(jì)算的條件，或者采用并行計(jì)算效果不明顯。理由如下:預(yù)測(cè)步驟中，分為存活目標(biāo)的預(yù)測(cè)和新生目標(biāo)的預(yù)測(cè)，對(duì)于存活目標(biāo)可以引入parfor 并行計(jì)算，但是因?yàn)槠湔加玫臅r(shí)延很小，故加速效果不明顯。重采樣的目的是為了去除權(quán)重小的粒子，復(fù)制權(quán)重大的粒子，在傳統(tǒng)的PPHD 的框架下不能用parfor 并行循環(huán)計(jì)算。聚類步驟是為了提取目標(biāo)狀態(tài)，其方法本身限制了parfor 循環(huán)的使用。

由式(5)～式(7)得出，在P-PHD 濾波的串行更新步驟中，需要計(jì)算所有量測(cè)對(duì)于每個(gè)粒子權(quán)重的影響，由于每個(gè)粒子都是相互獨(dú)立的，這就滿足了并行計(jì)算的條件。在預(yù)測(cè)步驟中，已經(jīng)得到了Lk=Lk-1+Jk個(gè)粒子的狀態(tài)和權(quán)重，即接下來(lái)是要對(duì)每個(gè)粒子權(quán)重的所有量測(cè)進(jìn)行并行更新。對(duì)于Lk=Lk-1+Jk個(gè)粒子，使用任務(wù)并行計(jì)算parfor 循環(huán)并行更新粒子權(quán)重，步驟如下:①開(kāi)啟并行parfor循環(huán)執(zhí)行步驟②～⑨;②對(duì)于n 個(gè)量測(cè)值Zm×n，循環(huán)③～⑦來(lái)計(jì)算其對(duì)粒子權(quán)重的影響; ③根據(jù)公式θm= arctan計(jì)算粒子對(duì)應(yīng)3 個(gè)傳感器的集中式融合量測(cè)值; ④根據(jù)③得出的融合量測(cè)值，計(jì)算似然g =N(θ1，Z(1，m)，σw(1))N(θ2，Z(2，m)，σw(2))N(θ3，Z(3，m)，σw(3));⑤由雜波平均數(shù)計(jì)算雜波強(qiáng)度kz=r/(2π)3;⑥ 由②中的似然 g 計(jì)算權(quán)值 w1= pdg/(kz+⑦累加所有量測(cè)對(duì)粒子權(quán)值的影響w2=w2+w1;⑧由公式得出更新后粒子權(quán)值;⑨最后得出

3 仿真結(jié)果與分析

假設(shè)在一個(gè)二維的場(chǎng)景中，監(jiān)視區(qū)域范圍為［-40，60］×［-40，60］。采樣周期選為T =1，共仿真40 步。目標(biāo)運(yùn)動(dòng)滿足線性模型，其運(yùn)動(dòng)方程為

式中，xk表示目標(biāo)運(yùn)動(dòng)狀態(tài)，其表達(dá)式為xk=［x，vx，y，vy］T，其中(x，y)表示目標(biāo)的位置信息，(vx，vy)表示目標(biāo)的速度信息;過(guò)程噪聲并且雜波平均數(shù)為λ=3，服從泊松分布，并假設(shè)其均勻地分布在整個(gè)觀測(cè)空間。本實(shí)驗(yàn)中設(shè)置3 個(gè)傳感器來(lái)跟蹤目標(biāo)，其位置坐標(biāo)分別為(-40，-40)，(10，60)，(60，-40)。假設(shè)傳感器只能獲得目標(biāo)的方位角θk，并且對(duì)方位角進(jìn)行集中式融合處理，觀測(cè)方程為

其中:x，y 分別對(duì)應(yīng)目標(biāo)的真實(shí)運(yùn)動(dòng)狀態(tài)的x 軸和y 軸坐標(biāo)，si表示第i 個(gè)傳感器。量測(cè)噪聲diag(［0.000 1 0.000 1 0.000 1］)，vk，wk相互獨(dú)立。雜波服從均值為r 的泊松分布。設(shè)置目標(biāo)檢測(cè)概率pD=1，無(wú)虛警，目標(biāo)的存活概率pS=0.99。假設(shè)沒(méi)有目標(biāo)衍生。表1 為4個(gè)目標(biāo)的初始狀態(tài)和存活時(shí)刻。

表1 目標(biāo)運(yùn)動(dòng)狀態(tài)

針對(duì)粒子濾波的更新步驟引入并行計(jì)算功能，并行更新每個(gè)粒子的權(quán)重，粒子數(shù)從500 ～5 000 分別開(kāi)啟1 ～4 核做并行計(jì)算，其中1 核即為串行執(zhí)行方式。各種情況的運(yùn)行時(shí)間統(tǒng)計(jì)如表2 所示。綜合表2 可以看出，并行模式下，核數(shù)越多，實(shí)時(shí)性越好。這是因?yàn)樵诟虏襟E中，對(duì)于每個(gè)粒子，都需要獨(dú)立計(jì)算其對(duì)于所有量測(cè)值更新之后的權(quán)重。并行計(jì)算可以并發(fā)的計(jì)算其權(quán)重，而不是順序執(zhí)行，進(jìn)而提高了實(shí)時(shí)性。需要指出的是，并行模式所用的核數(shù)與時(shí)間性能的提高并不是呈倍數(shù)關(guān)系，這是因?yàn)椴⑿心Ｊ街校⑿杏?jì)算的同步和結(jié)果的匯總也需要消耗時(shí)間。

表2 并行更新運(yùn)行時(shí)間

仿真環(huán)境為操作系統(tǒng)為win 7 旗艦版; 處理器為Intel(R)Core i5 -2400 CPU@3.10GHz; 內(nèi)存為2G，32 位操作系統(tǒng)。Matlab 版本為R2013a。本算法采用Wasserstein 距離作為性能指標(biāo)。以2 000 個(gè)粒子，開(kāi)啟雙核并行，對(duì)并行算法和串行粒子概率假設(shè)密度濾波進(jìn)行對(duì)比，仿真結(jié)果如圖1 ～圖4 所示。

圖1 給出了目標(biāo)做線性運(yùn)動(dòng)時(shí)的真實(shí)運(yùn)動(dòng)軌跡，其中位于坐標(biāo)軸上的3 個(gè)三角形表示傳感器的位置，標(biāo)明了其視場(chǎng)范圍;圖2 則給出了真實(shí)坐標(biāo)與所跟蹤的X 和Y 坐標(biāo)，在第10、11、28 時(shí)刻略差于串行P-PHD 濾波，第7 時(shí)刻略優(yōu)于串行P-PHD 濾波，其他時(shí)刻2 種算法所跟蹤的目標(biāo)位置基本相同。圖3 給出了整個(gè)觀測(cè)時(shí)間內(nèi)SMC-PHD 濾波算法以及并行算法對(duì)目標(biāo)數(shù)的估計(jì)和實(shí)際目標(biāo)數(shù)的對(duì)比，第7 時(shí)刻串行的P-PHD 濾波誤估計(jì)了目標(biāo)數(shù)。從圖4 可以看出，在第7、29 時(shí)刻并行算法的估計(jì)誤差比傳統(tǒng)的P-PHD 濾波算法大，在第6 時(shí)刻較串行P-PHD 稍小，其他時(shí)刻基本相同。實(shí)驗(yàn)結(jié)果可以表明并行算法和串行算法運(yùn)算結(jié)果稍有差異，這是因?yàn)閷?shí)驗(yàn)過(guò)程中每次仿真所產(chǎn)生的數(shù)據(jù)具有隨機(jī)波動(dòng)性，因而盡管并行算法滿足運(yùn)算結(jié)果的封閉性，但是不同次的仿真結(jié)果是會(huì)出現(xiàn)差異的。

圖1 目標(biāo)運(yùn)動(dòng)軌跡

圖2 估計(jì)目標(biāo)X 和Y 坐標(biāo)

圖3 目標(biāo)數(shù)

圖4 Wasserstein 距離

4 結(jié)束語(yǔ)

針對(duì)傳統(tǒng)P-PHD 更新步驟中采用串行更新粒子的權(quán)重所存在的不足，通過(guò)引入并行計(jì)算更新粒子的權(quán)重。大量的實(shí)驗(yàn)表明，在不失濾波性能的情況下，隨著核數(shù)的增加，更新步驟的運(yùn)算時(shí)間得以減少，運(yùn)算效率較串行計(jì)算有明顯提高。但在P-PHD 濾波中引入并行計(jì)算，要求數(shù)據(jù)必須是相互獨(dú)立的，否則不能實(shí)現(xiàn)并行計(jì)算。

［1］Mahler R. Multitarget Bayes filtering via first-order multitarget moments［J］.IEEE Transactions on Aerospace and Electronic Systems，2003，39(4):1152-1178.

［2］Vo B N，Singh S，Doucet A.Sequential Monte Carlo implementation of the PHD filter for multi-target tracking［C］//In Proceedings of the International Conference on Information Fusion.Cairns，Australia，2003:792-799.

［3］Vo B N，Ma W K.The Gaussian mixture probability hypothesis density filter［J］.IEEE Transactions on Signal Processing，2006，54(11):4091-4104.

［4］鄭云美.實(shí)時(shí)概率假設(shè)密度粒子濾波的算法研究［D］.浙江:浙江大學(xué)，2013.

［5］李威，韓崇昭，閆小喜.基于相對(duì)熵的概率假設(shè)密度濾波器序貫蒙特卡羅實(shí)現(xiàn)方式［J］. 控制與決策，2014，29(6):997-1002.

［6］張俊根，姬紅兵，蔡紹曉.并行高斯粒子濾波被動(dòng)多目標(biāo)跟蹤新算法［J］.西安電子科技大學(xué)學(xué)報(bào):自然科學(xué)版，2011，38(1):117-122.

［7］章濤，來(lái)燃，吳仁彪.觀測(cè)最優(yōu)分配的GM-PHD 多目標(biāo)跟蹤算法［J］.信號(hào)處理，2014，30(12):1419-1426.

［8］章濤，吳仁彪. 自適應(yīng)門限GM-PHD 多目標(biāo)跟蹤算法［J］.數(shù)據(jù)采集與處理，2014，29(4):549-554.

［9］王愛(ài)俠，李晶皎，王青，等.基于多核的并行粒子濾波運(yùn)動(dòng)目標(biāo)跟蹤［J］.計(jì)算機(jī)科學(xué)，2012，38(8):296-299.

［10］王美.粒子濾波并行處理的算法研究［D］.西安:西安電子科技大學(xué)，2014.

［11］孫偉平，向杰，陳加忠，等.基于GPU 的粒子濾波并行算法［J］.華中科技大學(xué)學(xué)報(bào):自然科學(xué)版，2011，39(5):63-66.

［12］Gong P，Basciftci Y O，Ozguner F. A Parallel Resampling Algorithm for Particle Filtering on Shared-Memory Architectures［C］//Parallel and Distributed Processing Symposium Workshops ＆ PhD Forum (IPDPSW). Shanghai，China，2012:1477-1483.