查巍巍 張 勇 應(yīng)曉慧 周 源 李知杰
(武漢船舶通信研究所 武漢 430064)
粒子濾波是一種基于蒙特卡羅方法和遞推貝葉斯估計(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ù)牧W佑肎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ì)精度。
粒子濾波器的本質(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)步。
分布式粒子濾波算法就是將式(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ì)算。
通過第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ù)。
用于粒子濾波的總粒子數(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)步。
本文在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.