查巍巍 張 勇 應(yīng)曉慧 周 源 李知杰

（武漢船舶通信研究所 武漢 430064）

1 引言

粒子濾波是一種基于蒙特卡羅方法和遞推貝葉斯估計(jì)的統(tǒng)計(jì)方法［1］，適用于任何用到非線性、非高斯遞推貝葉斯估計(jì)的系統(tǒng)，精度可以逼近最優(yōu)估計(jì)。然而由于其計(jì)算量龐大，若采用集中式結(jié)構(gòu)，粒子信息計(jì)算和通信所產(chǎn)生的巨大能耗會(huì)導(dǎo)致中心節(jié)點(diǎn)的癱瘓，因此有必要探索分布式的處理方法。

現(xiàn)有的分布式粒子濾波（Distributed Particle Fifter，DPF）主要包括兩個(gè)方面：一種基于信息融合技術(shù)［2］即每個(gè)節(jié)點(diǎn)獨(dú)立地運(yùn)行完整的粒子濾波形成一次狀態(tài)估計(jì)，然后交由融合中心進(jìn)行融合；另一種基于并行計(jì)算，即將粒子濾波算法分解成可并行計(jì)算的進(jìn)程，多個(gè)節(jié)點(diǎn)協(xié)作共同完成一次狀態(tài)估計(jì)。由于粒子濾波的計(jì)算量巨大，單個(gè)節(jié)點(diǎn)上獨(dú)立運(yùn)行完整的粒子濾波算法將產(chǎn)生大量能耗，因此在硬件實(shí)現(xiàn)上多偏向于并行計(jì)算。圖1給出了分布式粒子濾波并行計(jì)算模型原理圖。

Bashi等設(shè)計(jì)了三種分布式粒子濾波算法［3］：全局DPF（Global DPF）、局部 DPF（Local DPF）和壓縮 DPF（Compressed DPF）。三種方法都是基于并行計(jì)算，不同主要在于重采樣，全局DPF和壓縮DPF采用的是全局重采樣，局部DPF采用的是局部重采樣，然而若采用全局重采樣則需要上傳大量粒子的信息，增加了通信損耗，若采用局部重采樣，會(huì)出現(xiàn)部分子節(jié)點(diǎn)由于粒子退化對(duì)全局狀態(tài)估計(jì)作用甚小，導(dǎo)致最后全局估計(jì)誤差偏大。

圖1 分布式粒子濾波并行計(jì)算原理圖

為了克服全局重采樣帶來的大通信量弊端，有人提出了基于GMM模型的DPF［4］，主要改進(jìn)在于將傳輸?shù)牧Ｗ佑肎MM分布近似，然而每次狀態(tài)估計(jì)都需要計(jì)算GMM模型參數(shù)，會(huì)產(chǎn)生大量的計(jì)算能耗。因此，本文給出了一種精度能達(dá)到全局重采樣，而且不需要傳輸具體的粒子信息的分布式粒子濾波并行計(jì)算方法。其中重采樣算法來源于比例分配重采樣［5］，只是將該算法在粒子集上并行實(shí)現(xiàn)，并采用粒子群優(yōu)化算法解決部分子節(jié)點(diǎn)由于粒子退化帶來的計(jì)算量失衡的問題。仿真結(jié)果表明，本文的算法在保持負(fù)載均衡的同時(shí)減少了通信損耗，并且可以取得較好的估計(jì)精度。

2 粒子濾波算法

粒子濾波器的本質(zhì)是遞歸貝葉斯濾波，根據(jù)貝葉斯估計(jì)理論，一階馬 爾科夫序列｛xk｝k＝1，2，…的估 計(jì)問題是通過觀測序列｛zk｝k＝1，2，…和狀態(tài)轉(zhuǎn)移概率密度函數(shù)p（xk｜xk－1）得到后驗(yàn)概率密度函數(shù)p（xk｜zk），即：

遞推過程可分為預(yù)測和更新兩步

預(yù)測：

更新：

根據(jù)蒙特卡羅原理，后驗(yàn)概率分布可以用一系列帶權(quán)重的隨機(jī)采樣點(diǎn)逼近，即：

由于從后驗(yàn)概率中很難直接采樣，所以引入容易采樣的已知分布的重要密度函數(shù)：q（x0：k｜z1：k）。粒子從重要密度函數(shù)中采樣得到，則權(quán)值為

一般選擇狀態(tài)轉(zhuǎn)移概率密度函數(shù)p（xk｜xk－1）作為重要密度函數(shù)。

粒子濾波算法由四步組成，描述如下：

1）采樣。從p（xk｜xk－1）中采樣N 個(gè)粒子

2）權(quán)值計(jì)算。

并且歸一化

則可得k時(shí)刻的狀態(tài)估計(jì)為

3）重采樣。得到新的粒子集合

4）時(shí)刻k＝k＋1，轉(zhuǎn)到第1）步。

3 分布式處理

分布式粒子濾波算法就是將式（6）～（8）在粒子集上并行計(jì)算，這三個(gè)式子可以很好的在粒子集上進(jìn)行拆分計(jì)算，但是難點(diǎn)在于重采樣的并行處理，不同的處理思想直接影響到分布式粒子濾波的執(zhí)行效率和估計(jì)精度。

常見的重采樣算法有多項(xiàng)式重采樣算法［1］、殘差重采樣算法［6～7，9］、分層重采樣算法［6～8］、系統(tǒng)重采樣算法［6～9］。然而這些重采樣算法都建立在粒子權(quán)重的累積分布函數(shù)，需要對(duì)全部的權(quán)重信息集中排序，若將這些算法并行計(jì)算，則不可避免進(jìn)行粒子權(quán)重信息的交換，將帶來額外的通信能耗。

Bolic提出比例分配重采樣算法，該算法按粒子歸一化的權(quán)重進(jìn)行按比例采樣，建立在權(quán)重和之上，不需要交換每個(gè)粒子的權(quán)重信息，因此可以很好地進(jìn)行并行計(jì)算。

接下來可以進(jìn)行并行重采樣計(jì)算，偽代碼如下：

Step1：中心節(jié)點(diǎn)根據(jù)子節(jié)點(diǎn)的權(quán)重和進(jìn)行第一次一次比例分配重采樣。

Step2：中心節(jié)點(diǎn)將 ｛u′ ，C｝變量傳給子節(jié)點(diǎn)，子節(jié)點(diǎn)對(duì)本地粒子進(jìn)行第二次比例分配重采樣。

其中，中心節(jié)點(diǎn)產(chǎn)生的N＿resample（j）表示第j個(gè)節(jié)點(diǎn)上應(yīng)該重采樣N＿resample（j）個(gè)粒子，子節(jié)點(diǎn)上n＿resam－表示第j個(gè)子節(jié)點(diǎn)上第i個(gè)粒子應(yīng)該被重采樣n＿resa－次。而變量u′的計(jì)算剛好保證第j個(gè)子節(jié)點(diǎn)上通過step2得到的粒子數(shù)之和與中心節(jié)點(diǎn)計(jì)算的相同，即：

也就是說局部獨(dú)立重采樣的結(jié)果和全局重采樣的結(jié)果是一樣的，圖2為并行比例分配重采樣的計(jì)算原理圖。

圖2 并行比例分配重采樣原理圖

其中節(jié)點(diǎn)1和節(jié)點(diǎn)2的重采樣并行進(jìn)行，U和U′是并行重采樣用到的隨機(jī)數(shù)，u為集中式重采樣用到的隨機(jī)數(shù)，n（i）為第i個(gè)粒子重采樣后重復(fù)的個(gè)數(shù)，若

則并行重采樣和集中式重采樣能得到相同的結(jié)果，即相同的n值。

根據(jù)Step1得

由集中式比例分配重采樣算法得

即U′（1）＝u（6），以此類推，式（11）得證。也就是說比例分配重采樣并行計(jì)算的關(guān)鍵在于每個(gè)子節(jié)點(diǎn)上第一個(gè)隨機(jī)數(shù)的產(chǎn)生，即變量u′的計(jì)算。

4 分布式處理帶來的負(fù)載失衡問題

通過第3節(jié)可以很好地將粒子濾波算法并行計(jì)算，但是在若干次估計(jì)之后，一些子節(jié)點(diǎn)由于粒子退化使得權(quán)重值和很小，中心節(jié)點(diǎn)比例分配重采樣之后使得這些節(jié)點(diǎn)只需要產(chǎn)生少量的粒子數(shù)，即N＿resample很小，極端情況下，N＿resample只有一個(gè)非零，也就是說粒子全部在一個(gè)節(jié)點(diǎn)上計(jì)算，分布式粒子濾波退化成集中式粒子濾波。

上述問題可以理解為并行計(jì)算系統(tǒng)中的負(fù)載失衡問題，即有些節(jié)點(diǎn)超載，有些節(jié)點(diǎn)輕載，必須采用一種有效的算法均衡各個(gè)負(fù)載之間的計(jì)算任務(wù)。本文采用粒子群優(yōu)化算法（PSO）對(duì)出現(xiàn)輕載的節(jié)點(diǎn)上的粒子進(jìn)行優(yōu)化，使其總體權(quán)重和變大，那么在下一次狀態(tài)估計(jì)中，N＿resample變大，在一定程度上減輕負(fù)載失衡的問題。

粒子群優(yōu)化算法（PSO）采用兩個(gè)極值來更新粒子的速度和位置，使得整個(gè)粒子群向最優(yōu)位置靠近，一個(gè)為個(gè)體極值Pi，即粒子本身從初始態(tài)到當(dāng)前態(tài)經(jīng)歷的最優(yōu)位置，另一個(gè)為種群極值Pg，即整個(gè)粒子群經(jīng)歷過的最優(yōu)位置。利用粒子群優(yōu)化算法對(duì)采樣過程進(jìn)行優(yōu)化，使得粒子向高似然概率區(qū)域運(yùn)動(dòng)。本文采用一種改進(jìn)的PSO算法［10］，該方法基于一個(gè)高斯分布來不斷更新粒子的速度，其收斂性好于經(jīng)典的PSO算法。

一般狀態(tài)空間模型分為狀態(tài)轉(zhuǎn)移模型和量測模型：

其中wk和Vk分別為過程噪聲和量測噪聲，并且是相互獨(dú)立，協(xié)方差分別為Qk和Rk的零均值加性噪聲序列，因此真實(shí)值應(yīng)該在觀測值周圍服從均值為0，協(xié)方差Rk的噪聲分布，定義適應(yīng)度函數(shù)為

PSO更新過程為

經(jīng)過M次迭代后，粒子的最優(yōu)值符合某閾值ε時(shí)，說明粒子群已經(jīng)分布在真實(shí)狀態(tài)附近，則停止優(yōu)化。

圖3是一個(gè)用PSO算法解決負(fù)載失衡問題的示意，其中總粒子數(shù)為100，子節(jié)點(diǎn)數(shù)為5。其中大括號(hào)內(nèi)的值表示每個(gè)節(jié)點(diǎn)重采樣之后的粒子數(shù)。

圖3 利用PSO算法解決負(fù)載失衡的例子

引入PSO算法的目的是為了優(yōu)化節(jié)點(diǎn)上的粒子，在一定程度上解決負(fù)載失衡的問題，但同時(shí)帶來了額外的計(jì)算量，并且PSO優(yōu)化打破了比例分配重采樣的全局并行實(shí)現(xiàn)，而是用局部重采樣代替，這樣將會(huì)帶來估計(jì)誤差，在實(shí)際中需要權(quán)衡PSO優(yōu)化的次數(shù)。

5 完整的分布式粒子濾波算法

用于粒子濾波的總粒子數(shù)目為N，用于計(jì)算的子節(jié)點(diǎn)數(shù)目為S，下標(biāo)k表示時(shí)間點(diǎn)，i表示第i個(gè)粒子，j表示第j個(gè)子節(jié)點(diǎn)。則分布式粒子濾波算法描述如下：

1）子節(jié)點(diǎn)采樣。

2）子節(jié)點(diǎn)權(quán)值計(jì)算。

聚合數(shù)據(jù)如下：

4）中心節(jié)點(diǎn)進(jìn)行狀態(tài)估計(jì)。

5）中心節(jié)點(diǎn)進(jìn)行重采樣。判斷N＿resample，若N＿resample（j）／N＜α，則置位PSO優(yōu)化標(biāo)志，將數(shù)據(jù)｛u′，Ck，F(xiàn)｝發(fā)送到子節(jié)點(diǎn)，其中α為判斷負(fù)載是否失衡的閾值，F(xiàn)為PSO優(yōu)化的標(biāo)志位。

6）子節(jié)點(diǎn)進(jìn)行重采樣。若F相應(yīng)位被置位，那么相應(yīng)子節(jié)點(diǎn)先調(diào)用PSO算法，然后采用局部DPF中的局部重采樣；若F沒有被置位，則子節(jié)點(diǎn)進(jìn)行比例分配重采樣。然后回到第1）步。

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

本文在Matlab平臺(tái)下進(jìn)行仿真，假設(shè)目標(biāo)在二維平面內(nèi)做勻速直線運(yùn)動(dòng)。

系統(tǒng)狀態(tài)轉(zhuǎn)移方程為

觀測方程為

其中，F(xiàn)＝［1T00 0100 001T 0001］；G＝［T2／20 T0 0T2／2 0T］；Wk為過程噪聲，其協(xié)方差矩陣Q＝diag（0.12，0.12）；Vk為觀測噪聲，這里采用閃爍噪聲，其密度函數(shù)可表示為

協(xié)方差P＝（1－eta）＊P1＋eta＊P2，其中

取σr1＝σr2＝50；σb1＝1°，σb2＝5°；eta＝0.1

初始狀態(tài)選為x0＝［50000 300 50000 －100］T，仿真點(diǎn)為50，采用不同的粒子數(shù)和用于計(jì)算的子節(jié)點(diǎn)數(shù)，分別進(jìn)行100次蒙特卡羅實(shí)驗(yàn)，本文的算法和其它兩種分布式粒子濾波算法的誤差比較如表1所示。

表1 不同粒子數(shù)、子節(jié)點(diǎn)數(shù)時(shí)的誤差

表1中的結(jié)果可以表示成圖4～圖6的形式。分別對(duì)應(yīng)節(jié)點(diǎn)數(shù)為1、5和10時(shí)三種算法的誤差曲線。

圖4 1個(gè)節(jié)點(diǎn)的誤差曲線

圖5 5個(gè)節(jié)點(diǎn)的誤差曲線

圖6 10個(gè)節(jié)點(diǎn)的誤差曲線

對(duì)于GDPF誤差在粒子數(shù)增加到2000時(shí)可以降到很小，但是該算法需要傳輸大量的粒子，增加通信損耗；對(duì)于LDPF算法，由于局部重采樣導(dǎo)致誤差比同粒子數(shù)時(shí)GDPF的要大，但是不需要傳送任何粒子信息，實(shí)現(xiàn)方便；本文所采用的算法也不需要傳送任何粒子信息，并且由于引入了PSO算法，優(yōu)化了用于估計(jì)的粒子，因此誤差降低，但是PSO運(yùn)行的次數(shù)越多會(huì)帶來額外的計(jì)算開銷?？梢园l(fā)現(xiàn)在粒子數(shù)為100時(shí)，容易出現(xiàn)負(fù)載失衡，也就是調(diào)用PSO算法的次數(shù)會(huì)很多，誤差減小的同時(shí)帶來很大的PSO計(jì)算開銷，當(dāng)粒子數(shù)上升為2000時(shí)，調(diào)用PSO算法的次數(shù)為0，也就是沒有出現(xiàn)負(fù)載失衡的情況，此時(shí)誤差值和節(jié)點(diǎn)數(shù)為1時(shí)相近，即比例分配重采樣的并行計(jì)算可以達(dá)到與集中式計(jì)算同樣的精度，證明了并行處理過程是可行的。

［1］Gordon N，Salmond D.Novel appr oach to non－linear and non－Gaussian Bayesian state estimation［J］.Procof I nstitute E lectr ic E ngineering，1993，140（2）：107－113.

［2］Coates M.Distributed particle fiflters for sensor networks.PIPSN，2004：99－107.

［3］Bashi A S，Jilkov V P，Li X R，et al.Distributed implementations of particlefilters［C］／／Proceeding of the Sixth International Conference of Information Fusion，2003，2：1164－1171.

［4］Sheng X H，Hu Y H.Parameswaran Ramanathan.Distributed particle filter with GMM approximation for multiple targets［C］／／4th International Symposium on.Information Pocessing in Sensor Networks，2005：181－188.

［5］Bolic M，Djuric P M，Hong S J.Resamping alogorithms and architectures for distributed partical filters［J］.IEEE Trans.on Signal Processin，2005，53（7）：2442－2450.

［6］R V Merwe，A Doucet，Nando De Freitas，et al.The Unsented Particle Filter［R］／／Technical Report，CUED／F－INPENG／TR 380.UK：Engineering Department，Cambridge University，2000.

［7］J D Hol.Resampling in Particle Filters［R］／／Intership report，LiTH－ISY－EX－ET－0283－2004.Sweden：Link ping University，2004.

［8］Miodrag Bolic'，Peter M Djuric'，Sangjin Hong.New Resampling Algorithms for Particle Filters［J］.ICASSP2003（S1845－5921），2003，2：589－592.

［9］KROHLING R A.Gaussian swarm：a novel particle swarm optimization algorithm［C］／／Proceedings of the IEEE Conference on Cybernetics and Intelligent Systems（CIS）.Singapore：IEEE Press，2004：372－376.

［10］Peng zhang，Ming Li，Yan Wu.An improved particle filter algorithm based on Markov Random Field modeling in stationary wavelet domain for SAR image despeckling［J］.Pattern Recognition Letters，2012，33（10）：1316－1328.