基于蒙特卡洛方法的高斯混合采樣粒子濾波算法研究

1、基于蒙特卡洛方法的高斯混合采樣粒子濾波算法研究摘要本文提出了一種標(biāo)準(zhǔn)粒子濾波器的改進(jìn)算法高斯混合采樣粒子濾波算法GSPPF。仿真結(jié)果說(shuō)明，新算法在大幅降低計(jì)算復(fù)雜度的前提下，具有比標(biāo)準(zhǔn)粒子濾波算法(SIR-PPF)更好估計(jì)性能.關(guān)鍵詞卡爾曼濾波；粒子濾波；序列蒙特卡洛；貝葉斯濾波；高斯混合采樣貝葉斯方法為動(dòng)態(tài)系統(tǒng)的估計(jì)問(wèn)題提供了一類嚴(yán)謹(jǐn)?shù)慕鉀Q框架。它利用的信息建立系統(tǒng)的概率密度函數(shù)可以得到對(duì)系統(tǒng)狀態(tài)估計(jì)的最優(yōu)解。對(duì)于線性高斯的估計(jì)問(wèn)題，期望的概率密度函數(shù)仍是高斯分布，它的分布特性可用均值和方差來(lái)描繪。卡爾曼濾波器很好地解決了這類估計(jì)問(wèn)題1。對(duì)于非線性系統(tǒng)的估計(jì)問(wèn)題，最經(jīng)典并得到廣泛應(yīng)用的方法以

2、擴(kuò)展的卡爾曼濾波為代表，這類方法需要對(duì)模型進(jìn)展線性化，同時(shí)要求期望的概率密度函數(shù)滿足高斯分布，然而在對(duì)實(shí)際系統(tǒng)建模時(shí)，模型往往是非線性非高斯的。此時(shí)，最優(yōu)估計(jì)很難實(shí)現(xiàn)。粒子partile濾波器序列重要性采樣粒子濾波器，是一種適用于強(qiáng)非線性、無(wú)高斯約束的基于模擬的統(tǒng)計(jì)濾波器2。它利用一定數(shù)量的粒子來(lái)表示隨機(jī)變量的后驗(yàn)概率分布，從而可以近似得到任意函數(shù)的數(shù)學(xué)期望，并且能應(yīng)用于任意非線性隨機(jī)系統(tǒng)。本文介紹一種估計(jì)性能更好的粒子濾波算法高斯混合采樣粒子濾波器(GSPPF)，相比通常意義上的粒子濾波算法SIR-PF，GSPPF粒子濾波器具有更小的系統(tǒng)狀態(tài)估計(jì)的均方誤差和均值。貝葉斯濾波用概率統(tǒng)計(jì)的方法從

3、已觀察到的數(shù)據(jù)中獲得動(dòng)態(tài)狀態(tài)空間DSS模型參數(shù)。在DSS模型中，包含狀態(tài)和觀測(cè)兩個(gè)方程34。其中狀態(tài)轉(zhuǎn)移方程StateEquatin通常寫作1這里，是，且是白噪聲獨(dú)立的隨機(jī)序列，而且分布是的。觀測(cè)方程表達(dá)式寫為2這里：是白噪聲序列，獨(dú)立且分布。并且滿足。圖1描繪了DSS模型中狀態(tài)轉(zhuǎn)移和似然函數(shù)的關(guān)系。假設(shè)初始時(shí)刻系統(tǒng)的狀態(tài)分布，k時(shí)刻的信息序列表示。圖1動(dòng)態(tài)狀態(tài)空間模型DSS這樣，貝葉斯估計(jì)的問(wèn)題理解為：利用觀測(cè)到的信息Yk，求解系統(tǒng)狀態(tài)的概率分布。假設(shè)系統(tǒng)狀態(tài)的變化是隱馬爾柯夫過(guò)程，即當(dāng)前系統(tǒng)的狀態(tài)信息只與上一個(gè)時(shí)刻的狀態(tài)有關(guān)，可以通過(guò)預(yù)測(cè)和更新的途徑求解。3這里：4假設(shè)xk，k是互相獨(dú)立的

4、隨機(jī)變量，滿足。于是，參考1式可以把4式寫為5其中，是采樣函數(shù)。當(dāng)是時(shí)，xk可以通過(guò)確定性方程1得到。根據(jù)貝葉斯準(zhǔn)那么，系統(tǒng)狀態(tài)估計(jì)量6其中，7另外，在給定xk，vk，分布的條件下，yk的條件概率根據(jù)測(cè)量方程2可以表示為如下形式8由6式可以看出，后驗(yàn)概率密度包含3個(gè)部分。先驗(yàn)概率似然函數(shù)和證據(jù)。如何獲得這三項(xiàng)的近似是貝葉斯濾波的核心問(wèn)題。更新方程5中觀測(cè)值用來(lái)對(duì)的先驗(yàn)預(yù)測(cè)值修正，從而獲得狀態(tài)的后驗(yàn)概率。方程3和6的遞歸關(guān)系構(gòu)成了求解貝葉斯估計(jì)問(wèn)題的兩個(gè)步驟：預(yù)測(cè)與更新。假設(shè)(1)，(2)中的hk，fk是線性的，且噪聲k，vk滿足高斯白噪聲，可以把貝葉斯估計(jì)問(wèn)題簡(jiǎn)化為卡爾曼分析解。但這類問(wèn)題僅僅

5、是實(shí)際問(wèn)題中很小的一個(gè)部分。對(duì)于更多的問(wèn)題，很難得到分析解。只有通過(guò)對(duì)問(wèn)題的近似線性處理(擴(kuò)展卡爾曼濾波)或其它途徑蒙特卡洛方法實(shí)現(xiàn)非線性、非高斯問(wèn)題的解。根據(jù)后面分析問(wèn)題需要，這里重點(diǎn)對(duì)蒙特卡洛方法積分進(jìn)展說(shuō)明。在過(guò)去的二十多年，蒙特卡洛方法得到了很大的開展。其優(yōu)點(diǎn)就是用系列滿足條件的采樣點(diǎn)及其權(quán)重來(lái)表示后驗(yàn)概率密度。蒙特卡洛方法采用統(tǒng)計(jì)抽樣和估計(jì)對(duì)數(shù)學(xué)問(wèn)題進(jìn)展求解。按照其用途，可以把蒙特卡洛方法分為三類5：蒙特卡洛抽樣、計(jì)算、優(yōu)化。其中，蒙特卡洛抽樣是尋找有效的、方差很小的、用于估計(jì)的抽樣方法。蒙特卡洛計(jì)算那么是設(shè)計(jì)產(chǎn)生滿足特定要求隨機(jī)數(shù)的隨機(jī)發(fā)生器的問(wèn)題。而蒙特卡洛優(yōu)化是采用蒙特卡洛思想

6、對(duì)實(shí)際中的非凸非差分函數(shù)優(yōu)化求解。對(duì)于，可以由概率空間p(x)中抽取N個(gè)樣本，用近似值作為的解。大數(shù)定理證明：收斂于，并且滿足條件。這里，是的方差。不同于確定性的數(shù)字計(jì)算，蒙特卡洛近似的一個(gè)重要特點(diǎn)就是估計(jì)的精度獨(dú)立于狀態(tài)空間的維數(shù)。而且，積分估計(jì)的方差與采樣點(diǎn)的個(gè)數(shù)成反比。顯然，蒙特卡洛近似方法的關(guān)鍵點(diǎn)有兩個(gè)：首先如何由一個(gè)樣本空間中抽取N個(gè)采樣點(diǎn)，用來(lái)表征后驗(yàn)概率密度。其次就是計(jì)算。重要性抽樣IprtantSapling解決了如何借助于分布來(lái)對(duì)實(shí)現(xiàn)有效采樣的問(wèn)題，由arshall1965年提出。當(dāng)數(shù)據(jù)空間非常宏大時(shí)，重要性抽樣只對(duì)其中“重要區(qū)域進(jìn)展采樣，節(jié)省了計(jì)算量。對(duì)于高維采樣空間模型，

7、如統(tǒng)計(jì)物理學(xué)、貝葉斯統(tǒng)計(jì)量，這一點(diǎn)尤為重要。重要性抽樣的中心思想是選擇一個(gè)覆蓋真實(shí)分布px的建議分布qx8。這樣，(9)對(duì)qx作蒙特卡洛抽樣，假設(shè)粒子數(shù)目為N，有(10)其中，稱為重要性權(quán)重，再作歸一處理，(11)是歸一化權(quán)重。為了減小估計(jì)的方差，選擇的建議性分布qx與px盡可能匹配。通常，建議分布qx需要一個(gè)長(zhǎng)的拖尾，這樣可以解決區(qū)間之外的干擾。確切的說(shuō)，匹配的qx必須與pxfx成正比9。當(dāng)qx與px不匹配時(shí)，x(i)是不均勻分布的，在整個(gè)遞歸迭代的過(guò)程中，存在大量的權(quán)值極小的樣本，而這些樣本對(duì)估計(jì)的奉獻(xiàn)很校事實(shí)上，權(quán)值較大的少數(shù)樣本決定蒙特卡洛采樣的估計(jì)精度。大量時(shí)間損耗在這些“無(wú)關(guān)緊要的

8、粒子計(jì)算上，即所謂的粒子退化現(xiàn)象DegenerayPrble。目前，標(biāo)準(zhǔn)的粒子濾波器選擇先驗(yàn)概率Prir作為建議分布。對(duì)于粒子退化現(xiàn)象，采樣重要性重采樣方法給出了很好的解決途徑。其根本思想就是通過(guò)在兩次重要性采樣之間增加重采樣步驟，消除權(quán)值較小的樣本，并對(duì)權(quán)值較大的樣本復(fù)制，降低了計(jì)算的復(fù)雜度。在(N)時(shí)間復(fù)雜度范圍內(nèi)可以已排序的均勻分布序列作重采樣處理。對(duì)重采樣Resapling處理，新的采樣結(jié)果放在數(shù)組，詳細(xì)的算法用偽碼語(yǔ)言寫為如下的形式：步驟1：令這里必須注意是隨機(jī)變量的累計(jì)概率密度序列。步驟2：初始假設(shè)，當(dāng)，產(chǎn)生一組序列分布。對(duì)一個(gè)固定的j，分別用逐一比較，一旦，就可以得到一組新的樣本

9、集合。如此循環(huán)直到。需要說(shuō)明的是，重采樣方法在消除粒子退化問(wèn)題的同時(shí)，也帶來(lái)了其它兩個(gè)問(wèn)題：首先，降低了粒子運(yùn)算并行執(zhí)行的可能性；其次，由于權(quán)值較大的粒子屢次被選擇，粒子的多樣性減少。這種情況尤其在小過(guò)程噪聲條件下表現(xiàn)更為明顯11。圖2SIR-PF重要性采樣與重采樣示意圖如前所述，利用序列重要性采樣和重采樣的方法，粒子濾波可以有效的遞歸更新后驗(yàn)概率的分布。但是，由于對(duì)粒子未加假設(shè)，大量的粒子在處理非線性、非高斯問(wèn)題時(shí)出現(xiàn)了計(jì)算的高復(fù)雜性問(wèn)題。另外，由于少數(shù)權(quán)值較大的粒子反復(fù)被選擇，粒子坍塌明顯。文獻(xiàn)4提出了在重要性采樣步驟的建議分布的生成階段“搬運(yùn)粒子到似然較高區(qū)域，可以緩解坍塌，同時(shí)進(jìn)步估計(jì)

10、的性能。但是不可防止的是對(duì)每一個(gè)粒子的后驗(yàn)概率處理，使得計(jì)算的復(fù)雜性進(jìn)一步加劇。鑒于此種情況，這里介紹一種新穎的高斯混合采樣粒子濾波器GaussianixtureSigaPintPartileFilter，GSPPF。GSPPF算法利用有限高斯混合模型表征后驗(yàn)概率分布情況，可以通過(guò)基于重要性采樣的加權(quán)的后驗(yàn)粒子，借助于加權(quán)的期望最大化算法EightedExpetinaxiizatin交換標(biāo)準(zhǔn)重采樣步驟，降低粒子坍塌效應(yīng)。4.1基于高斯混合近似的采樣卡爾曼濾波器根據(jù)最優(yōu)濾波理論，一個(gè)概率密度p(x)都可以寫作高斯混合模型Gaussianixturedel。即，這里，G是高斯分量的個(gè)數(shù)，是高斯分量

11、的權(quán)重，是以向量為均值，以p(g)為協(xié)方差矩陣的隨機(jī)向量x的高斯分布。考慮DSS狀態(tài)轉(zhuǎn)移方程和觀測(cè)方程，假設(shè)先驗(yàn)概率及噪聲密度服從高斯混合模型G。這樣，預(yù)測(cè)的先驗(yàn)概率密度滿足，更新后。這里，。在此根底之上，預(yù)測(cè)的先驗(yàn)概率和后驗(yàn)概率對(duì)應(yīng)的均值和方差可以通過(guò)采樣卡爾曼濾波器SigaPintKF計(jì)算。4.2基于觀測(cè)更新的重要性采樣IprtantSapling前已敘及重要性抽樣是一種蒙特卡洛方法，即用一組帶有權(quán)值的樣本數(shù)據(jù)來(lái)表征隨機(jī)變量的概率密度。利用DSS模型的一階馬爾柯夫本質(zhì)和給定狀態(tài)的觀測(cè)值依賴性，可以推導(dǎo)遞歸的權(quán)值更新方程，這里僅對(duì)于給定的粒子而言。在GSPPF算法中，用G近似來(lái)。作為建議分布

12、。由于包含了最新的樣本數(shù)據(jù)，使得粒子聚集在高似然區(qū)域，一定程度減少了粒子坍塌效應(yīng)。另外，使用預(yù)測(cè)的先驗(yàn)概率平滑權(quán)值更新方程中的，這是因?yàn)镚SPPF算法用G表示后驗(yàn)概率，本次后驗(yàn)同時(shí)又是下一個(gè)時(shí)間步的先驗(yàn)概率，G模型中高斯核對(duì)后驗(yàn)概率做了平滑處理。基于觀測(cè)更新步驟的重要性采樣方法中對(duì)粒子不作任何假設(shè)，對(duì)非線性、非高斯問(wèn)題具有很強(qiáng)的魯棒性。4.3采用加權(quán)的E算法做重采樣和G復(fù)原基于觀測(cè)更新步驟的重要性采樣輸出是一組加權(quán)的粒子，在標(biāo)準(zhǔn)的粒子濾波器中，這些粒子必須作重采樣處理丟棄小權(quán)值粒子，同時(shí)對(duì)權(quán)值較大的粒子做放大處理。通過(guò)這種處理，可以有效的防止粒子集合的方差增加太快。不幸的是，重采樣步驟只對(duì)當(dāng)觀

13、測(cè)似然微弱、大量粒子聚集極少數(shù)粒子副本情況有效。在GSPPF算法中，采用加權(quán)的期望最大EightedExpetinaxiizatin直接得到G模型，實(shí)現(xiàn)對(duì)加權(quán)粒子的最大似然擬合，這就相當(dāng)于對(duì)粒子的后驗(yàn)概率做了平滑，防止了粒子坍塌問(wèn)題，同時(shí)，G模型中的高斯核的個(gè)數(shù)減少到G，防止其呈指數(shù)級(jí)增長(zhǎng)，降低了算法復(fù)雜度。為了比較算法的性能，系統(tǒng)狀態(tài)估計(jì)的條件均值，均方誤差Errrnvariane可以通過(guò)兩個(gè)方法計(jì)算，即在加權(quán)的E算法平滑之前，用下面公式求解，描繪了系統(tǒng)的均值與均方誤差性能。這里，給定系統(tǒng)狀態(tài)估計(jì)問(wèn)題的算法評(píng)估模型(12)是噪聲，。另外，非平穩(wěn)觀測(cè)模型(13)，其中，觀測(cè)噪聲服從高斯分布。假

14、設(shè)給定含噪的系統(tǒng)狀態(tài)觀測(cè)值yk，采用兩種不同的算法：標(biāo)準(zhǔn)的粒子濾波算法SIR-PF以及GSPPF算法對(duì)系統(tǒng)的狀態(tài)xk估計(jì)。每次實(shí)驗(yàn)共做150次，每次的觀察樣本重新產(chǎn)生，SIR-PF算法中粒子的個(gè)數(shù)是250個(gè)。GSPPF算法中采用兩種方案：第一種方案用5個(gè)高斯核擬合狀態(tài)后驗(yàn)概率。狀態(tài)噪聲vk，觀測(cè)噪聲nk各用一個(gè)高斯核擬合。第二種方案那么用3個(gè)高斯核擬合Gaa(3，2)分布的拖尾狀態(tài)噪聲，這里擬合方法采用E算法。圖3、圖4描繪了系統(tǒng)的隱狀態(tài)和觀測(cè)值及SIR-PF，GSPPF算法系統(tǒng)狀態(tài)的估計(jì)值。圖3SIR-PF粒子濾波器狀態(tài)估計(jì)圖4GSPPF粒子濾波器狀態(tài)估計(jì)采用4.3部分的均方誤差和均值計(jì)算公

15、式對(duì)不同算法對(duì)系統(tǒng)狀態(tài)估計(jì)性能作了比對(duì)。圖3、圖4曲線說(shuō)明，在系統(tǒng)的觀測(cè)噪聲nk均方誤差很小，而過(guò)程噪聲服從具有長(zhǎng)的拖尾分布時(shí)，采用轉(zhuǎn)移概率作為建議分布的標(biāo)準(zhǔn)粒子濾波器性能很差。這是因?yàn)橛^測(cè)方程中峰值似然函數(shù)和系統(tǒng)狀態(tài)急劇的跳躍變化產(chǎn)生的結(jié)果。盡管可以通過(guò)采樣卡爾曼(Siga-Pint)濾波器將粒子向似然峰值區(qū)域搬動(dòng)解決這一問(wèn)題，但是也使得計(jì)算量加大。GSPPF算法兩種不同方案都具有比SIR-PF更好的系統(tǒng)狀態(tài)估計(jì)性能，均方誤差比后者數(shù)量級(jí)降低了1/103-1/104。與1個(gè)高斯核擬合過(guò)程噪聲的GSPPF算法比較，3個(gè)高斯核擬合算法性能更好，但時(shí)間復(fù)雜度同樣有所進(jìn)步。由于GSPPF算法在大幅度

16、降低了算法的計(jì)算復(fù)雜度同時(shí)，可以獲得準(zhǔn)確的系統(tǒng)估計(jì)性能。所以說(shuō)，GSPPF算法為粒子濾波理論實(shí)時(shí)應(yīng)用，如目的定位單目的與多目的、時(shí)變信道估計(jì)、圖像增強(qiáng)、機(jī)器故障診斷以及語(yǔ)音信號(hào)處理等提供了一個(gè)新的方案。1Y.HandR.K.Lee，ABayesianapprahtprblesinsthastiestiatinandntrlIEEETrans.Autat.ntr.vl.A-9.pp.333-3392A.Duet，N.Freitas，N.Grdn.SequentialntearlethdsinPratie.Springer3B.D.AndersnandJ.B.re.ptialfiltering.P

17、rentieHallEngldliff，NJ.19794N.J.Grdn，D.J.Salnd，A.F.Sith，Nvelapprahtnnlinear/nn-GaussianBayesianstateestiatin，IEEpreedingsvl140，N2，April19935ULLER，“ntearlintegratiningeneraldynaidels“ntep.ath.1991，115，pp，145-1636FredriGustafssn，NilasBergan，PartilefiltersfrPsitin，NavigatinandTraking“，F(xiàn)inalversinfrIEEE

18、TransatinsnSignalPressingSpeialissuenntearlethdsfrstatistialsignal7ARNAUDDUET，SINGDSILL，nsequentialntearlsaplingethdfrBayesianfiltering“statistisandputing(2000)，10，197-208，reivedJuly1998andaeptedAugust19998JayeshH.KtehaandPetar.Djurri，GaussianPartileFilteringInPr.rkshpStatistialSignalPress.Singapre，Aug.20019J.S.LiuR.hen.SequentialntearlethdsfrDynaialSystes.JurnalftheAerianStatistialAssiatin，1998，Vlue93.pp.1032-104410ZHEHEN，Bayesianfiltering：FrKalanfiltertpartilefil