張宏欣，周穗華，馮士民

(海軍工程大學(xué)兵器工程系，湖北武漢430033)

?

基于徑向基-Galerkin解的反饋粒子濾波器

張宏欣，周穗華，馮士民

(海軍工程大學(xué)兵器工程系，湖北武漢430033)

摘要：反饋粒子濾波器(FPF)是一種連續(xù)時間最優(yōu)貝葉斯估計器.針對FPF對系統(tǒng)采樣率要求較高的問題，提出一種基于徑向基函數(shù)-Galerkin法求解的反饋粒子濾波器.該算法推導(dǎo)了反饋增益勢函數(shù)所滿足偏微分方程的弱形式，結(jié)合徑向基函數(shù)對勢函數(shù)進行近似，利用Galerkin法和蒙特卡羅積分得到了反饋增益近似解，并給出了一種徑向基參數(shù)選取方法，從數(shù)值上分析了徑向基函數(shù)參數(shù)選取對于濾波精度的影響.仿真算例表明反饋粒子濾波器在低系統(tǒng)采樣率下會嚴重發(fā)散，而本文算法能夠避免這一問題，且提高了FPF在低系統(tǒng)采樣率下的濾波精度和穩(wěn)定性.

關(guān)鍵詞：非線性濾波;貝葉斯濾波;反饋粒子濾波器; Galerkin法;徑向基函數(shù)

1　引言

粒子濾波是一類基于遞推貝葉斯法則，利用樣本統(tǒng)計量對系統(tǒng)狀態(tài)進行估計的方法.在處理非高斯非線性系統(tǒng)的參數(shù)/狀態(tài)濾波問題方面具有獨特的優(yōu)勢，已被廣泛應(yīng)用于目標(biāo)跟蹤、通信信號處理以及計算機視覺等領(lǐng)域［1，2］.

一般針對粒子濾波算法研究主要基于文獻［3］中提出的自舉濾波(Bootstrap Filter)框架，即序貫建議分布采樣(SIS)或建議分布采樣-重采樣(SIR)，若建議分布產(chǎn)生的樣本中不能足量覆蓋后驗分布樣本，則容易導(dǎo)致粒子貧化.針對貧化問題已有大量改進算法［4～8］，主要針對重采樣改進或(和)建議分布改進.重采樣改進，如馬爾科夫蒙特卡洛(MCMC)采樣、輔助粒子濾波(APF)、正則化粒子濾波(RPF)等，目的在于從建議樣本中更準(zhǔn)確的抽取后驗樣本，并未從根本上克服貧化問題;建議分布改進，如無跡卡爾曼-粒子濾波(UKFPF)、容積卡爾曼-粒子濾波(CKF-PF)及混合粒子濾波(MPF)等，目的是構(gòu)造更接近后驗概率密度的建議分布，這類改進方法對于不同系統(tǒng)模型的性能差異較大，難以達到貝葉斯意義下的最優(yōu)估計.

根據(jù)以上研究發(fā)現(xiàn)，傳統(tǒng)粒子濾波及其改進算法的主要問題在于建議分布難于構(gòu)造，以及重采樣方法在多維空間難于抽樣［5］.反饋粒子濾波器［9～13］(Feedback Particle Filter，F(xiàn)PF)是一種連續(xù)時間最優(yōu)貝葉斯估計器.通過對先驗樣本構(gòu)造反饋，采用Kullback-Leibler散度(KLD)來度量反饋后概率分布與真實后驗分布之間的差異，將最小化KLD等效為歐拉-拉格朗日邊值題求解反饋增益，使得由先驗樣本能夠直接得出后驗樣本，從而不需要再構(gòu)造建議分布和進行重采樣過程.

由于FPF在連續(xù)時間條件下導(dǎo)出，本身需要較小的采樣間隔(≤0.05)來保證其精度，而現(xiàn)有算法在求解過程中采用了常數(shù)反饋增益近似，這種方法在不合適的采樣間隔條件下將產(chǎn)生誤差累積，導(dǎo)致精度下降甚至濾波發(fā)散，針對此問題，本文提出一種能適用于較大采樣間隔的FPF算法.利用徑向基函數(shù)［14］(Radial Basis Function，RBF)良好的多維非線性函數(shù)擬合能力近似反饋增益的勢函數(shù)，結(jié)合穩(wěn)定性較好的Galerkin法［15～19］，且Galerkin法基于積分形式，因此利用Monte Carlo積分即可以較小的計算代價將問題轉(zhuǎn)化為矩陣方程.相比常數(shù)反饋近似，RBF-Galerkin解提高了反饋增益近似精度，避免了常數(shù)反饋近似在較大采樣間隔下的誤差累積，放寬了FPF算法對采樣間隔的要求，具有較高的濾波精度和較好的穩(wěn)定性.

2　反饋粒子濾波器

考慮如下連續(xù)時間非線性濾波問題:

其中，xt∈d為d維狀態(tài)矢量，yt∈m為m維觀測矢量;函數(shù)f:d→d，h:d→m; wt∈d，vt∈m分別為系統(tǒng)噪聲和觀測噪聲，均服從互不相關(guān)的零均值高斯分布，協(xié)方差矩陣分別為函數(shù)h的第j個分量為函數(shù)hj，hj:d→，j =1，2，…，m.

非線性濾波的目的是根據(jù)歷史觀測估計出狀態(tài)xt的后驗分布.連續(xù)時間條件下的后驗分布p*定義為，對于任意d空間上的可測集A，有

其中Zt為歷史觀測的σ域，即Zt: =σ(Zs: s≤t).

令dt≈Δt，此時后驗分布p*可寫為離散時間的遞推貝葉斯形式:

其中py|x為似然函數(shù)，py為與狀態(tài)無關(guān)的歸一化常數(shù).px為狀態(tài)先驗分布.在離散時間條件下，對先驗粒子疊加反饋，t時刻加入反饋后的粒子為:

其中Kt:d→d×m為反饋增益函數(shù)，K = { K1，K2，…，Km}，Kj:d→d; It為更新項(Innovation Term)，i = 1，2，…，N，N為粒子個數(shù)，

其中，D為p與真實后驗分布p*的K-L信息散度:

式(8)是一個變分問題，可以導(dǎo)出［11］，當(dāng)Δt→0，若p(x，0)= p*(x，0)，則p(x，t)= p*(x，t)，t＞0.函數(shù)Cj(x，t):d→d是滿足如下偏微分方程(PDE)的解，

則最優(yōu)反饋函數(shù)的每個元素為:

其中Clj(x，t)為Cj(x，t)的第l個元素，函數(shù)矩陣C(x，t)= { C1(x，t)，C2(x，t)，…，Cm(x，t)} .

疊加反饋后的粒子的連續(xù)時間系統(tǒng)方程為:

式(9～11)稱為反饋式粒子濾波器(FPF).

3　反饋增益的RBF-Galerkin解

由前述可知，F(xiàn)PF算法的關(guān)鍵在于求解PDE(9).RBF求解PDE常用的是配點法［14］(Collocation Methods)，即選取分別滿足邊界條件和PDE解的一組點直接代入PDE求解;而對于本文中的二階時變PDE，配點法會使得(1)系數(shù)矩陣條件數(shù)過大;(2)需要擬合出概率分布p(x，t).這均會導(dǎo)致求解過程中的數(shù)值不穩(wěn)定.Galerkin法基于雙線性形式，其系數(shù)矩陣是數(shù)值穩(wěn)定的.綜合考慮RBF的擬合能力和Galerkin法的穩(wěn)定性，本節(jié)導(dǎo)出反饋增益的RBF-Galerkin解.

3.1Galerkin弱形式的導(dǎo)出

記p = p(x，t)，Cj= Cj(x，t)，hj= hj(x).假設(shè)函數(shù)Cj代表的場是無旋的，即▽×Cj= 0，則存在勢函數(shù)φj(x，t)使得Cj=▽φj，結(jié)合式(9)可得:

記L2(d; p)為d中對概率密度p(·，t)平方可積的函數(shù)構(gòu)成的希爾伯特空間; Hk(d; p)為d上k階導(dǎo)數(shù)在L2(d; p)中的函數(shù)構(gòu)成的希爾伯特空間.若存在φj∈H1(d; p)，使得對于任意H1(d; p)中的函數(shù)ψ(x，t):d→，有

式(13)為式(9)的弱形式，ψ=ψ(x，t)稱為試函數(shù).式(13)左端可展開為:

對式(14)右端第一項利用高維分部積分可得:

其中Γ為積分邊界，假設(shè)邊界上的積分為0，將式(15)代入式(14)得:

則(13)可重寫為:

式(17)稱為PDE(12)在H1(d; p)上的Galerkin弱形式，Galerkin法的目的即是求解φj∈H1(d; p)，使得對任意ψ ∈H1(d; p)滿足式(17).將式(17)寫為期望形式如下:

文獻［13］中給出的常數(shù)近似增益定義為:

3.2增益解的RBF-Galerkin近似

RBF為一組在d中徑向?qū)ΨQ的標(biāo)量值函數(shù)，為了在由RBF張成的子空間Vn?H1(d; p)中逼近目標(biāo)函數(shù)，將增益函數(shù)的勢近似為RBF的線性組合，略去與求解無關(guān)的變量下標(biāo)，則:

其中，θε(·)為RBF，為RBF的中心點，λj，j =1，2，…，L為RBF下待求解的坐標(biāo)，ε為形參數(shù)，由于Galerkin法求解要求θε(·)∈H1(d; p).因此選擇高斯型徑向函數(shù)，

若要在Vn中近似φj，則需對任意試函數(shù)ψ∈Vn滿足式(18)，為此試函數(shù)選擇為RBF的任意線性組合:

其中aj，j = 1，2，…，L為任一組不為零的常數(shù)，記，并將式(22)和(20)代入(18)，得到

消去ai，式(23)可轉(zhuǎn)化為關(guān)于λ的線性方程組Aλ= b，

此時系數(shù)矩陣A是對稱的.進一步利用Monte-Carlo積分，式(24)可以由粒子閉合地得出，

由上式可解出λ，并注意到Cj=▽φj，可對每個粒子求得RBF-Galerkin反饋增益近似解:

3.3RBF參數(shù)選擇

RBF-Galerkin方法需要選擇的參數(shù)為基函數(shù)個數(shù)L，形參數(shù)ε以及中心點c = { c1，c2，…，cL} .目前關(guān)于ε 和cj的選擇已有廣泛研究，文獻［18］中采用貪心算法給出了一種最優(yōu)的cj選擇方法，但計算量較大，不適合濾波求解.而對于基中心點的選擇，一方面需使其位于增益解較大的區(qū)域［19］，此時形參數(shù)ε的選擇可以放寬;另一方面由于濾波算法執(zhí)行必須基于離散化模型，在高維和低采樣率條件下，F(xiàn)PF算法必須在較小的區(qū)域進行求解，初始狀態(tài)誤差需取得較小，否則會出現(xiàn)嚴重發(fā)散，此時解是足夠平滑的，這要求基中心點靠近積分邊界，在此條件下，RBF應(yīng)足夠“平”，即參數(shù)ε需取很小的值以保證求解穩(wěn)定性.對于基個數(shù)，較大的L會使得系數(shù)矩陣條件數(shù)過大，數(shù)值上不穩(wěn)定;過小則會使精度下降甚至濾波發(fā)散.

濾波問題通常希望參數(shù)能夠在每個時刻通過樣本統(tǒng)計特性確定，綜上，本文選擇基中心點為，

對于形參數(shù)，采用與文獻［19］中類似的選取方法，引入無量綱的Mahalanobis距離，選取形參數(shù)

其中ravg為狀態(tài)粒子的平均Mahalanobis距離，α為尺度因子，該參數(shù)決定了RBF的平滑程度.

RBF相關(guān)參數(shù)選擇將在數(shù)值仿真部分中進一步分析.需要指出，RBF參數(shù)與具體問題模型有關(guān)，目前仍未有確定的通用方法.

4　結(jié)果及分析

本節(jié)通過對一個非線性目標(biāo)運動模型［13］進行仿真實驗，對比分析FPF算法和本文算法的濾波性能.該模型描述目標(biāo)以恒定角速度做平面螺旋運動，狀態(tài)矢量xt=［xt，yt］T∈2為t時刻目標(biāo)的平面直角坐標(biāo)，目標(biāo)狀態(tài)滿足如下隨機微分方程，

其中vt，1，vt，2為相互獨立的狀態(tài)噪聲分量，且vt，i～為狀態(tài)噪聲標(biāo)準(zhǔn)差;

其中，λ，Θ和ρ為實值參量; 1(ρ，∞)(·)為指示函數(shù)

仿真初始條件為x0=［0.5，-0.5］T;模型參數(shù)λ= 2，Θ= 50，ρ= 9.狀態(tài)初始分布為p0= N(［0，0］T，42I2×2)，σv=0.1，σw= 0.27，粒子數(shù)為100.仿真時間T 為15s.RBF參數(shù)為α=0.0006，κ=20.

圖1為Δt =0.06時，目標(biāo)軌跡與兩種算法一次跟蹤軌跡對比，可見，F(xiàn)PF算法在該采樣間隔下出現(xiàn)了偏差;而本文算法則能夠準(zhǔn)確地跟蹤目標(biāo)軌跡.

觀測模型為

為分析FPF算法和本文算法在不同采樣間隔下的濾波性能，在相同初始條件下Δt分別取0.04，0.06，0.08，0.1時對兩種算法進行100次Monte

對比結(jié)果如圖2(a)～(d)所示，圖中橫坐標(biāo)為Monte Carlo仿真序號，縱坐標(biāo)為RMSE.由圖中可以看出，隨著Δt增大，F(xiàn)PF算法由于誤差累積產(chǎn)生了較大的估計誤差，逐漸出現(xiàn)了濾波發(fā)散現(xiàn)象，Δt =0.1時，F(xiàn)PF算法在大部分實驗中嚴重發(fā)散，在第63次和76次仿真實驗中出現(xiàn)了超過20m的RMSE值(出于整體直觀考慮未在圖中示出)，而本文算法在不同采樣間隔條件下的RMSE是穩(wěn)定的，且在所有實驗中都保持了較小的估計誤差.

圖3(a)～(b)分別給出了Δt = 0.05條件下，兩種算法疊加反饋后粒子的歸一化似然函數(shù)值.從圖中可看出，約在10.5s之后，F(xiàn)PF反饋后的粒子的似然函數(shù)值均接近于0，這表明在一定觀測噪聲方差條件下，由于常數(shù)反饋增益解的不精確性，使得反饋后的部分粒子仍位于低值似然區(qū)域，在逐步迭代的過程中引起了誤差累積，導(dǎo)致在某一時刻全部粒子均位于似然區(qū)域之外;而本文算法反饋后的粒子在全仿真時間內(nèi)都具有較高似然值，表明RBF-Galerkin反饋增益解精度高于常數(shù)增益解.

表1中給出了不同采樣間隔條件下兩種算法100次Monte Carlo計算的RMSE均值和方差.從表中數(shù)據(jù)可知，采樣間隔越大，F(xiàn)PF的濾波精度和穩(wěn)定性越差，在較大采樣間隔時均出現(xiàn)了嚴重發(fā)散;而本文算法在不同采樣間隔條件下都保持了較高的濾波精度和穩(wěn)定性，且兩者均較大幅度優(yōu)于FPF算法.同時可以看到，采樣間隔越小，兩種算法的濾波精度差異越小，且均能取得較高的濾波精度.由此可見:(1)反饋增益函數(shù)所滿足的PDE是在連續(xù)時間條件下導(dǎo)出的，當(dāng)Δt→0時，PDE(9)本身是精確的，在此條件下，常數(shù)增益近似產(chǎn)生的反饋誤差并不會導(dǎo)致誤差累積;(2)作為濾波問題“最優(yōu)”準(zhǔn)則的最小均方誤差(MMSE)等同于貝葉斯方法中的后驗均值估計，若反饋增益能使得粒子“移動”到顯著似然區(qū)域，則此時反饋增益解的不精確對濾波解精度的影響是較小的［21］.(3)在采樣間隔較小時，精確求解PDE(9)能夠避免誤差累積.

表1中同時給出了各個采樣間隔條件下兩種算法的平均計算時間，該時間在MATLAB R2011b(7.13.0.564)仿真軟件和2.5GHz雙核計算機環(huán)境下測得.可知本文算法對采樣時間的放寬和精度的提高是以增加計算量為代價的.與經(jīng)典FPF算法相比，需要計算出Galerkin系數(shù)矩陣，因此其時間復(fù)雜度要一定程度地高于經(jīng)典FPF算法.Carlo計算，按下式統(tǒng)計狀態(tài)估計的均方根誤差(RMSE)，

表1　100次Monte Carlo計算結(jié)果：RMSE均值，方差及計算時間

為分析不同RBF參數(shù)設(shè)置對濾波精度的影響，取Δt =0.06，取α和κ為在一定范圍內(nèi)的變化值，對于每個參數(shù)取值進行100次Monte Carlo計算，得到每個參數(shù)取值條件下的RMSE均值如圖4所示，從圖中可看出，當(dāng)κ值較小(20)時，在0.002～0.0004的α值范圍內(nèi)均能取得較高的濾波精度，但在α值進一步減小時，較小的κ值會導(dǎo)致濾波結(jié)果不穩(wěn)定;而隨著κ值增大和α減小，可取得較平穩(wěn)且精度較高的結(jié)果;圖中同時給出了FPF算法獨立于參數(shù)取值進行100次Monte Carlo計算結(jié)果，可以看出在較大參數(shù)范圍內(nèi)，本文算法均能取得高于FPF算法的濾波精度.

5　總結(jié)

針對反饋粒子濾波器中常數(shù)反饋增益近似方法在較大采樣率下濾波發(fā)散的問題，本文通過徑向基函數(shù)擬合反饋增益的勢函數(shù)，結(jié)合Galerkin法和Monte Carlo積分推導(dǎo)了反饋增益的RBF-Galerkin解，給出了一種徑向基函數(shù)參數(shù)選取方法，得出了一種適用于較大采樣間隔的FPF算法，仿真結(jié)果表明算法精度較高，穩(wěn)定性好，復(fù)雜度適中，同時也表明:(1)本文算法的反饋增益精度高于常數(shù)近似方法，能顯著提高FPF算法在低系統(tǒng)采樣率下的濾波性能，放寬了FPF算法對于采樣間隔要求;(2)系統(tǒng)采樣率較低時，精確求解反饋增益函數(shù)能夠避免誤差累積導(dǎo)致濾波發(fā)散;(3)對于濾波問題，RBF參數(shù)與問題模型、狀態(tài)噪聲、觀測噪聲以及初始條件都有關(guān)，如何完全自適應(yīng)的選取最優(yōu)RBF參數(shù)是一個十分復(fù)雜的問題，這是需要進一步研究的問題.

參考文獻

［1］Djuric P M，Kotecha J H，Zhang J Q.Particle filtering［J］.IEEE Signal Processing Magazine，2003，20(5): 19 -38.

［2］Micheal I，Andrew B.Conditional density propagation for visual tracking［J］.International Journal of Computer Vision，1998，29(1): 5 -28.

［3］Gordon J.Novel approach to nonlinear/non-Gaussian Bayesian state estimation［A］.IEEE Proceedings on Sensor Fusion［C］.IEEE，1993.107 -113.

［4］Zhe Chen.Bayesian filtering: from Kalman filters to particle filters［J］.Statistics，2003，182(1): 1 -69.

［5］Arulampalam M S，Simon M，Neil G.A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking［J］.IEEE Transactions on Signal Processing，2002，50(2): 174 -188.

［6］Fred D.Nonlinear filters: beyond Kalman filter［J］.IEEE Aerospace and Electronic System Magazine，2005，43(8): 57 -69.

［7］Rudolph M.Sigma-point Kalman filters for probabilistic inference in dynamic state-space models［D］.Oregon: Oregon Health＆Science University，2004.251 -256.

［8］Yardim C，Michalopoulou Z H，Gerstoft P.An overview of sequential Bayesian filtering in ocean acoustics［J］.IEEE Journal of Oceanic Engineering，2011，36(1): 71 -89.

［9］Tao Y，Mehta P G，Meyn S P.Feedback particle filter with me an-field coupling［A］.IEEE 50th Conference on Decision and Control［C］.Florida: IEEE，2011.7909 -7916.

［10］Tao Y，Mehta P G，Meyn S P.A mean-field control-oriented approach to particle filtering［A］.American Control Conference［C］.San Francisco，2011.2037 -2043.

［11］Tao Y，Richard S，Sean M，et al.Multivariable feedback particle filter［A］.IEEE 51st Annual Conference on Decision and Control［C］.Florida: IEEE，2012.4063 -4070.

［12］Tao Y，Mehta P G，Meyn S P.Feedback particle filter［J］.IEEE Transactions on Automatic Control，2013，58(10):2037 -2043.

［13］Tilton A.A comparative study of nonlinear filtering techniques［A］.IEEE 16th International Conference on Information Fusion［C］.Istanbul: IEEE，2013.1827 -1834.

［14］Larsson E.A numerical study of some radial basis function based solution methods for elliptic PDEs［J］.Computers＆Mathematics with Applications，2003，46(1): 891 -902.

［15］Zienkiewicz C.The Finite Element Method for Solid and Structure Mechanics［M］.Oxford: Elsevier Press，2005.17 -25.

［16］Chen J S，et al.A stabilized conforming nodal integration for Galerkin meshfree methods［J］.International Journal for Numerical Methods in Engineering，2001，50(2):435 -466.

［17］Wenland H.Meshless Galerkin methods using RBFs［J］.Mathematics of Computation，1999，68(228 ): 1521 -1531.

［18］Ling L，Schaback R.An improved subspace selection algorithm for meshless collocation methods［J］.International Journal of Numerical Methods Engineering，2009，80: 1623 -1639.

［19］Katharina K，Elisabeth L.A Galerkin radial basis function method for the Schrodinger equation［J］.SIAM Journal of Scientific Computations，2013，35(6): 2832 -2855.

［20］張賢達.矩陣分析與應(yīng)用［M］.北京: Springer-清華大學(xué)出版社，2005.225 -227.Zhang XianDa.Matrix Analysis and Application［M］.Beijing: Springer-Tsinghua University Press，2005.225 -227.(in Chinese)

［21］Daum F.Coulomb＇s law particle flow for nonlinear filter ［A］.Proceedings of SPIE on Signal Processing and Sensor Fusion［C］.San Diego: SPIE，2011.3351 -3362.

張宏欣男，1987年12月出生，陜西漢中人.2010年畢業(yè)于西安理工大學(xué)，現(xiàn)為海軍工程大學(xué)博士生，從事統(tǒng)計信號處理及目標(biāo)跟蹤相關(guān)研究.

周穗華男，1962年10月出生，廣東五華人，1984年畢業(yè)于海軍工程學(xué)院，1990年在海軍工程學(xué)院獲得博士學(xué)位.現(xiàn)為海軍工程大學(xué)教授，從事軍用目標(biāo)特性信息處理及武器系統(tǒng)總體設(shè)計方面研究.

Feedback Particle Filter Using Radial Basis Functions Based Galerkin Method

ZHANG Hong-xin，ZHOU Sui-hua，F(xiàn)ENG Shi-min
(Department of Weapon Engineering，Naval University of Engineering，Wuhan，Hubei 430033，China)

Abstract:A radial basis function(RBF)Galerkin solution based feedback particle filter is proposed to resolve the divergency problem existing in present particle filter when the continuity of system model is violated.A weak formulation of the PDE regarding to the potential of feedback gain is firstly derived，then the RBFs are employed to approximate the potential function.Finally the feedback gain solution is obtained using Galerkin method and Monte Carlo integral，also the method for choosing RBF parameters is provided and analyzed numerically.It is demonstrated that the present FPF diverges under low system sample rate，whereas our proposed feedback particle filter is nevertheless effective，with preferable tracking accuracy and stability under low system sample rate.

Key words:nonlinear filtering; Bayesian filtering; feedback particle filer(FPF); Galerkin method; radial basis functions(RBF)

作者簡介

收稿日期:2014-04-14;修回日期: 2014-08-03;責(zé)任編輯:李勇鋒

DOI:電子學(xué)報URL：http: / /www.ejournal.org.cn10.3969/j.issn.0372-2112.2016.01.014

中圖分類號：TP202

文獻標(biāo)識碼：A

文章編號：0372-2112(2016)01-0095-06