高斯粒子流濾波器

張宏欣，周穗華，馮士民

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

高斯粒子流濾波器

張宏欣，周穗華，馮士民

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

粒子流濾波器以粒子流速度場描述隨機(jī)樣本從先驗分布到后驗分布的演化，實現(xiàn)對系統(tǒng)狀態(tài)的貝葉斯估計.針對其一般解計算復(fù)雜、難于濾波求解的問題，導(dǎo)出一種高斯假設(shè)條件下的粒子流濾波器.在線性高斯條件下推導(dǎo)了速度場的解析解；證明了當(dāng)演化步長趨近于0時，該解析解與Kalman-Bucy濾波器的解具有一致的形式；基于該解導(dǎo)出了非線性高斯系統(tǒng)速度場的表達(dá)式，并進(jìn)一步利用Unscented變換近似求解.通過若干仿真算例表明，高斯粒子流濾波器放寬了系統(tǒng)噪聲為高斯型的限制，其精度優(yōu)于經(jīng)典非線性高斯濾波器，計算復(fù)雜度低于一般粒子濾波器，且具有良好的穩(wěn)定性.

非線性濾波；貝葉斯估計；粒子流濾波器；速度場；Unscented變換

1 引言

貝葉斯非線性濾波基于非線性動態(tài)系統(tǒng)模型，通過帶噪聲的觀測采樣序列近似計算出目標(biāo)狀態(tài)的后驗分布統(tǒng)計量，來實現(xiàn)對系統(tǒng)狀態(tài)的估計［1］，近幾十年來已經(jīng)在目標(biāo)跟蹤、導(dǎo)航等工程領(lǐng)域得到了廣泛應(yīng)用.根據(jù)近似方法的不同，可分為解析近似方法和仿真近似方法.解析法基于卡爾曼濾波器（Kalman Filter，KF）框架，假設(shè)后驗分布為高斯型，近似計算非線性變換后的狀態(tài)均值和誤差矩陣，利用卡爾曼濾波框架進(jìn)行求解，其中，擴(kuò)展卡爾曼濾波（Extended Kalman Filter，EKF）基于一階偽線性化，對于高階模型精度不高；以無跡卡爾曼濾波（Unscented Kalman Filter，UKF）為代表的σ點卡爾曼濾波器［2，3］采用一組確定的樣本點更精確的捕捉均值和方差信息，對高斯分布可達(dá)到三階精度；但對于某些高階模型（如指數(shù)函數(shù)，三角函數(shù)等），在高維情況下仍可能引入較大誤差，而捕捉更高階矩信息則需要更龐大的Sigma點集和克服數(shù)值問題［4］.

仿真近似方法采用大量隨機(jī)樣本經(jīng)過特定“變換”來近似后驗分布.粒子濾波［5］（Particle Filter，PF）從建議分布中抽取滿足后驗分布樣本（重采樣），但并未改變樣本本身，因此建議分布樣本需足量覆蓋后驗分布樣本，否則將導(dǎo)致濾波性能不穩(wěn)定甚至發(fā)散［6］，該問題在狀態(tài)維數(shù)較高時尤為明顯，通常需要大量隨機(jī)樣本（隨問題維數(shù)呈指數(shù)增長［7］）.針對此問題，反饋型的粒子濾波算法通過對每個先驗分布樣本疊加反饋量來構(gòu)造建議分布樣本，普遍做法是對每個樣本進(jìn)行非線性卡爾曼濾波，從而衍生出EKF-PF及UKF-PF［8］濾波器，這類算法的時間復(fù)雜度較高，且在重要性權(quán)值計算中很容易出現(xiàn)數(shù)值不穩(wěn)定情況.

粒子流濾波器［9～12］（Particle Flow Filter，PFF）是一種無需構(gòu)造建議分布和重采樣過程的反饋式粒子濾波器，通過構(gòu)造同倫（Homotopy）函數(shù)，將先驗分布粒子到后驗分布粒子的變化過程描述為粒子的“流動”，結(jié)合描述概率分布隨時間變化的Fokker-Planck方程（FPE），得到粒子流速度場（以下簡稱速度場）的偏微分方程，求出速度場即可利用數(shù)值方法得到后驗分布樣本；其優(yōu)勢在于無需構(gòu)造建議分布和重采樣過程，且適于求解高維問題［6］.文獻(xiàn)［9］中通過直接求解泊松方程給出了粒子流速度場理論表達(dá)式的Monte-Carlo （MC）積分解，其計算復(fù)雜度過高，難于實現(xiàn).

與文獻(xiàn)［9］不同，本文從粒子流濾波器的基本理論出發(fā)，推導(dǎo)出一種高斯型粒子流濾波器（Gaussian Particle Flow Filter，GPFF）.假設(shè)系統(tǒng)噪聲為高斯分布的基礎(chǔ)上，根據(jù)速度場所滿足的偏微分方程，在線性高斯條件下推導(dǎo)了速度場的解析解，證明了該解與連續(xù)時間條件下的Kalman-Bucy濾波解具有一致的形式；并將該解推廣到非線性系統(tǒng)，推導(dǎo)了速度場的期望形式解，結(jié)合無跡變換（Unscented Transform，UT）進(jìn)行近似求解.由于GPFF采用粒子進(jìn)行非線性傳播，放寬了對系統(tǒng)噪聲的高斯型限制［13］，相比EKF和UKF而言具有較高的精度；保留了PFF無需構(gòu)造建議分布及重采樣的優(yōu)點，相對PF來講降低了所需的粒子數(shù)量，減小了計算量，且穩(wěn)定性得到提高.通過對若干問題模型進(jìn)行仿真驗證了GPFF的有效性.

2 反饋粒子濾波器

考慮如下濾波問題模型：

若將狀態(tài)向量x視為R Rd空間標(biāo)準(zhǔn)基下的某一點的位置坐標(biāo)，則f（x，λ）可以看做λ時刻該點上的速度，若f對所有粒子都有定義，則它表示粒子從先驗分布到后驗分布“流動”的速度場.因此，若能求得f，即可通過數(shù)值積分方法計算后驗粒子.

假設(shè)變化過程不存在噪聲，根據(jù)連續(xù)時間隨機(jī)濾波理論，q（x，λ）滿足零擴(kuò)散項（diffusion term）的Fokker-Planck方程：

其中▽·為散度算子.對式（4）取自然對數(shù)得，

對式（7）求關(guān)于λ的偏導(dǎo)數(shù)，并結(jié)合式（6）可得：

式（9）是難處理的，這是因為：（a）獲得合適的β值是困難的［9］；（b）q（x，λ）求值需要采用分布擬合算法，計算代價過高；（c）格林公式僅適用于d＞2情況；（d）已經(jīng)證明［14］，這種形式的解在數(shù)值上是不穩(wěn)定的.因此需要尋求其他方法.事實上，文獻(xiàn)［9］中僅給出了使用式（9）計算速度場的結(jié)果，并未給出基于該式的濾波結(jié)果.

3 粒子流速度場在線性高斯條件下的解

3.1解析解推導(dǎo)

對于線性系統(tǒng)，假設(shè)先驗概率密度 g（x）滿足高斯分布；似然函數(shù)h（x）為高斯型，即

分別為λ條件下的后驗均值及協(xié)方差矩陣.在線性高斯條件下，設(shè)粒子流速度場f（x，λ）具有如下形式

其中A（λ）∈R Rd×d，b（λ）∈R Rd.記f=f（x，λ），注意到式（8）右端的形式，則可將其重寫為

將式（11）、（12）和（14）代入式（8）可得：

可得

根據(jù)式（18）右端的形式，可設(shè)

將式（19）和（20）代入式（15），可以得到如下關(guān)系

為減小計算量，利用如下關(guān)系

則A（λ）可僅由先驗協(xié)方差矩陣表示如下

將式（25）和式（26）代入式（20），則b（λ）可僅由先驗信息表示為

式（13）、（24）和（27）給出了粒子流速度場在線性高斯條件下的解析解.

3.2速度場與Kalman-Bucy濾波器的一致性

由考慮如下的連續(xù)時間線性隨機(jī)系統(tǒng)

其中xt為連續(xù)時間狀態(tài)向量，F(xiàn)t為漂移項；Ht為觀測矩陣，wt，vt分別為滿足零均值的高斯過程的狀態(tài)噪聲和觀測噪聲，且

對于連續(xù)時間線性高斯系統(tǒng)，Kalman-Bucy濾波器［16］以線性微分方程形式給出了最小均方無偏估計解及其協(xié)方差矩陣Pt，即

利用式（26）求Pλ關(guān)于λ的導(dǎo)數(shù)，可得

ε為一任意小的正數(shù).將式（37）代入式（36），化簡得

利用式（23），得到Pλ關(guān)于λ的導(dǎo)數(shù)為

對式（22）取關(guān)于x的期望可得，

綜上可知線性高斯條件下λ變化過程中，速度場與零漂移項和零過程噪聲系統(tǒng)的Kalman-Bucy濾波器解是一致的.

4 非線性濾波

設(shè)非線性條件下速度場的解具有如下形式：

將h（x）在狀態(tài)預(yù)測值^xk|k-1附近做一階泰勒展開，即

將式（41）和式（42）代入式（14），利用2.1節(jié)的推導(dǎo)過程可得

為將式（43）表示為統(tǒng)計量形式，利用式（23），并注意到

其中

式（41）和（46）給出了速度場在非線性高斯條件下的近似解.

采用無跡變換（UT）近似求解式（44）及式（45），UT通過一組確定的采樣點來近似計算隨機(jī)變量經(jīng)過非線性變換后的均值和方差，首先生成Sigma點以及相應(yīng)用于計算均值的權(quán)值和用于計算方差的權(quán)值即：

關(guān)于狀態(tài)的統(tǒng)計量直接通過粒子實現(xiàn)，則先驗統(tǒng)計量為

其中x（i）k|k-1為第i個先驗粒子，N為粒子個數(shù).利用導(dǎo)出的速度場計算公式，對狀態(tài)在λ∈（0，1］內(nèi)進(jìn)行數(shù)值積分即可求出后驗粒子，相應(yīng)的后驗統(tǒng)計量可根據(jù)式（51）和（52）類似地得到.

5 仿真結(jié)果及分析

為驗證 GPFF的有效性，對三個典型算例，即連續(xù)Wiener過程加速度（CWPA）模型、一維非穩(wěn)定增長（UNG）模型以及只測角跟蹤（BOT）模型進(jìn)行仿真，CWPA模型為線性模型，用于在高斯噪聲和非高斯噪聲條件下將GPFF與KF進(jìn)行比較；UNG模型用于在精度上將GPFF與EKF、UKF及PF等經(jīng)典非線性濾波器進(jìn)行比較.BOT模型考察多維和非理想初值條件下各濾波器性能.GPFF的粒子數(shù)量均為100個，步長Δλ=0.05.

5.1平面CWPA模型

平面CPWA模型為目標(biāo)加速度受隨機(jī)擾動的連續(xù)時間運動模型，以隨機(jī)微分方程描述如下：

其中，

k=1，2，…，T，T為觀測點數(shù)，

噪聲協(xié)方差矩陣

Δt為離散化時間間隔.

以平面位置坐標(biāo)作為觀測量，觀測方程為

觀測矩陣H=（I2×202×4）.設(shè)Δt=0.1，T=80.初始狀態(tài)為x0=06×1，濾波器初始分布為N（x0，P0|0），初始誤差矩陣P0|0=diag（0.1，0.1，0.1，0.1，0.5，0.5）.分別在高斯噪聲和混合高斯噪聲條件下進(jìn)行仿真，即

為直觀比較，采用由下式定義的維度無關(guān)均方誤差（MSE）

圖1（a）～（b）分別給出了在不同類型觀測噪聲下，KF與GPFF的100次MC計算MSE對比.由圖可知，在高斯噪聲條件下，GPFF與KF誤差近似相同，這與2.2中證明的結(jié)論是一致的；而在非高斯噪聲條件下，KF誤差增幅較大，且變得不穩(wěn)定，而GPFF仍具有較小的誤差，這表明GPFF保持了粒子型算法對系統(tǒng)噪聲的適應(yīng)性.

圖2給出了非高斯噪聲條件下GPFF與KF的跟蹤結(jié)果對比，可知在非高斯觀測噪聲下，KF跟蹤軌跡與真實軌跡產(chǎn)生了較大偏差，而GPFF仍可準(zhǔn)確地跟蹤目標(biāo).

5.2UNG模型

該模型是非線性濾波中的基準(zhǔn)問題［5］，其狀態(tài)空間模型為

圖3給出了各個算法 MSE的100次 MC仿真計算結(jié)果，其中PF粒子數(shù)N=100.由于UNG是一個不含采樣時間的完全離散模型，且包含高階非線性項，導(dǎo)致EKF的一階線性化誤差較大，穩(wěn)定性較差；GPFF在MSE和穩(wěn)定性上均顯著地優(yōu)于EKF以及同樣采用UT變換的UKF.

表1進(jìn)一步給出了各個算法100次MC計算的MSE均值和方差以及單次平均計算時間.其中，對于PF 和GPFF分別給出了 N=100，N=500和N=1000時的計算結(jié)果（如括號內(nèi)所示），執(zhí)行時間在MATLAB R2013（8.1.0.604）版本和2.5GHz雙核計算機(jī)條件下測得.由表可知，作為高斯型濾波器，GPFF精度優(yōu)于EKF和UKF；而作為粒子型算法，GPFF的計算量隨粒子數(shù)量的增長變化是較慢的；PF采用重采樣方法，其計算時間隨粒子數(shù)量增長較快，因此GPFF的計算效率高于PF.綜合圖表結(jié)果看，由于PF對樣本分布不基于任何假設(shè)，對于本例的情況（一維模型和理想初始條件），其精度最高.而在算例3中將會看到，對于多維問題，PF的性能將會大幅下降.

表1100　次MC計算結(jié)果：MSE均值，方差，計算時間

5.3BOT模型

BOT模型［15］是只以方位角測量值作為觀測量的跟蹤模型，其狀態(tài)模型為離散Wiener速度模型，

觀測噪聲vk～N（02×1，R），噪聲協(xié)方差矩陣R= 0.052I2×2.觀測點坐標(biāo)，設(shè)Δt=0.1s，時刻k=1，2，…，T，T= 200.狀態(tài)初值x0=（0，36，0.2，-1）T，濾波器初始分布為N（^x0，P0|0），初始估計值^x0=（-6，17，0，0）T，初始誤差矩陣P0|0=（25，25，1，1）T.PF粒子個數(shù)N=500.

采用各算法進(jìn)行100次Monte-Carlo計算，對于k時刻的某一個狀態(tài)分量εk，按下式計算各個狀態(tài)分量的MSE

其中Nm為Monte-Carlo仿真次數(shù)，^εk為k時刻狀態(tài)估計值.由上式可知，MSEkε是100次Monte-Carlo計算的平方誤差在特定時刻上的均值，可衡量算法隨觀測值信息增加的收斂性能.其計算結(jié)果如圖4（a）～（d）所示，由圖可知，EKF、UKF和GPFF均能夠收斂到較小的誤差范圍，而PF則出現(xiàn)了顯著偏差，且隨時間變化的趨勢是起伏的，可見在多維模型和不理想初始值條件下，經(jīng)典PF的性能較差；進(jìn)一步地可以看出，相較EKF和 UKF，GPFF能夠更快的收斂，且誤差較小.

圖5（a）～（d）給出了各個算法MSE的100次MC計算結(jié)果，可以看到，PF對于目標(biāo)縱軸坐標(biāo)的估計產(chǎn)生了很大誤差，而對于目標(biāo)橫軸速度估計的MSE則較小，且PF穩(wěn)定性較差，這種現(xiàn)象是由于y的估計初值與真實初值偏差最大，而.x的估計初值與真實初始偏差很小造成的，這從圖4中也可得出.

圖4和圖5表明一定粒子數(shù)量條件下，對于多維估計問題，PF的性能對于初始條件是十分敏感的，若初始條件不夠理想的，其平均性能甚至?xí)佑贓KF算法.而GPFF對于每個分量估計的MSE均明顯地低于其他算法，具有較高的精度，且在所有仿真中都保持了穩(wěn)定的性能.

表2中給出了各個算法條件下，100次MC計算所得的狀態(tài)分量估計MSE均值，方差和執(zhí)行時間，其中對于PF分別給出了粒子數(shù)時的結(jié)果.可知各算法中GPFF 的MSE均值為最小，且穩(wěn)定性好；UKF精度高于 EKF；PF算法在粒子數(shù)量為4000時給出了較好的估計結(jié)果，但此時其精度仍低于GPFF，表明在不理想初值和/或多維模型條件下，PF若無足夠粒子來代表多維樣本空間，或建議分布本身就與后驗分布相去甚遠(yuǎn)（不理想的初值），則重采樣方法此時是低效的.

表2　MSE均值、標(biāo)準(zhǔn)差和執(zhí)行時間

6 總結(jié)

針對粒子流濾波器中速度場難于計算的問題，基于粒子流濾波器基本理論，導(dǎo)出一種高斯粒子流濾波器.在線性高斯條件下推導(dǎo)了速度場的解析解；證明了當(dāng)粒子流演化步長趨近于0時速度場解析解與Kalman-Bucy濾波器解的一致性；并將算法推廣到非線性系統(tǒng)，利用Unscented變換近似求解.仿真算例表明，算法放寬了系統(tǒng)噪聲為高斯型的限制，其精度優(yōu)于EKF 和UKF等經(jīng)典非線性高斯濾波器，計算復(fù)雜度低于一般粒子濾波器，且性能穩(wěn)定，具有一定的應(yīng)用價值.

［1］James C.Bayesian Signals Processing［M］.Hoboken：John Wiley&Sons Press，2009.3-10.

［2］Simon J，Jeffery U，Hugh D.A new method for the nonlinear transformation of means and covariance in filters and estimators［J］.IEEE Transactions on Auto Control，2000，45（3）：477-481.

［3］Arasaratnam I，Haykin S.Cubature Kalman filters［J］. IEEE Transactions on Auto Control，2009，56（6）：1254 -1269.

［4］Rudolph M.Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models［D］.Oregon：Oregon Health&Science University，2004.251-256.

［5］Gordon N，Salmond D，Smith A.Novel approach to nonlinear/non-Gaussian Bayesian state estimation［J］.IEE Radar and Signal Processing，1993，140（2）：107-113.

［6］Fred D.Nonlinear filters：Beyond Kalman filter［J］.IEEE Aerospace&Electronics Systems Magazine，2005，43（8）：57-69.

［7］Synder C，Bengtsson T，Bickel P，et al.Obstacles to highdimensional particle filtering［J］.Monthly Weather Review，2008，136（1）：4629-4640.

［8］Rudolph M.The Unscented Particle Filter［R］.Cambridge：Cambridge University，2000.

［9］Daum F.Coulomb's law particle flow for nonlinear filter［A］.Proceedings of SPIE on Signal Processing and Sensor Fusion［C］.USA：SPIE Press，2011.3351-3362.

［10］Daum F.Numerical experiments with Coulomb’s law for particle flow in nonlinear filters［A］.Proceedings of SPIE Conference［C］.San Diego：SPIE Press，2011.2137 -2146.

［11］Fred D，Jim H.Particle flow with non-zero diffusion for nonlinear filters［A］.Proceedings of SPIE on Signal Processing and Sensor Fusion［C］.USA：SPIE Press，2013. 1745-1755.

［12］Tao D，Mark C.Implementation of the daum-huang exactflow particle filter［A］.IEEE Statistical Signal Processing Workshop［C］.Ann Arbor：IEEE Press，2012.257-261.

［13］Kotecha J，Djuric P.Gaussian particle filtering［J］.IEEE Transactions on Signal Processing，2003，51（10）：2592 -2601.

［14］Tao Y，Mehta P G，Meyn S P.Feedback particle filter ［J］.IEEE Transactions on Automatic Control，2013，58 （10）：2037-2043.

［15］Bar-Shalom Y，Xiaorong Li.Estimation with Applications to Tracking and Navigation［M］.New York：John Wiley &Sons Press，2005.50-54.

［16］Simon D.Optimal State Estimation［M］.New York：John Wiley&Sons Press，2006.22-27.

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

E-mail：mylifeforthebattle@hotmail.com

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

Gaussian Particle Flow Filter

ZHANG Hong-xin，ZHOU Sui-hua，F(xiàn)ENG Shi-min

（Department of Weapon Engineering，Naval University of Engineering，Wuhan，Hubei 430033，China）

Particle flow filter formulated the dynamics from prior samples with posterior samples with particle flow velocity field to perform Bayesian estimation of system state.To address difficulties of particle velocity field computation in present particle flow filter，a novel particle flow filter based on Gaussian assumption was proposed.The analytical solution of velocity field under linear Gaussian condition was derived.The consistency of this analytical solution with Kalman-Bucy filter for continuous system，when discrete dynamic step goes to zero，was proved.The solution was finally extended to obtain the nonlinear Gaussian velocity field expression which can be approximated by using unscented transformation.Several simulations revealed the effectiveness over classic nonlinear Gaussian assumption on accuracy and particle filter on efficiency and stability.

nonlinear filtering；Bayesian estimation；particle flow filter；velocity field；Unscented transformation

TP202

A

0372-2112（2016）04-0795-09

電子學(xué)報URL：http：//www.ejournal.org.cn 10.3969/j.issn.0372-2112.2016.04.007

2014-10-28；

2015-03-23；責(zé)任編輯：孫瑤