王 瑞

(西安電子科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，陜西西安 710071)

非線性非高斯?fàn)顟B(tài)空間模型［1－3］的最優(yōu)估計(jì)［4］在信號(hào)處理、金融、無(wú)線通訊等領(lǐng)域具有重要的應(yīng)用。針對(duì)實(shí)際問(wèn)題中出現(xiàn)的非線性非高斯模型，出現(xiàn)了一種被稱(chēng)為“序貫重要性采樣(SIS)”的Monte Carlo方法，其通過(guò)離散的隨機(jī)測(cè)度逼近概率分布。然而由于高度的計(jì)算復(fù)雜性和退化問(wèn)題，1993年 Gordon提出了重采樣概念，克服了早期算法的退化問(wèn)題，出現(xiàn)了第一個(gè)可操作的Monte Carlo濾波器?，F(xiàn)代計(jì)算技術(shù)使Monte Carlo濾波方法得到了迅速發(fā)展，現(xiàn)在通稱(chēng)為粒子濾波器［5］。本文是在非線性非高斯?fàn)顟B(tài)空間模型中假設(shè)某些參數(shù)未知，采用粒子濾波算法及改進(jìn)后的輔助粒子濾波算法分別對(duì)狀態(tài)進(jìn)行估計(jì)。粒子濾波算法［6－8］源于 Monte Carlo［9－10］的思想，其核心思想是通過(guò)從后驗(yàn)概率中抽取的隨機(jī)狀態(tài)粒子來(lái)表達(dá)其分布，是一種順序重要性采樣法(Sequential Importance Sampling)。

本文介紹了非線性混合效應(yīng)狀態(tài)空間模型，粒子濾波算法相關(guān)的序貫重要性采樣(SIS)及采樣重要性重抽樣(SIR)算法，并引入一種新的改進(jìn)粒子濾波算法—輔助粒子濾波算法(APF)［11－12］，通過(guò)實(shí)例模型進(jìn)行仿真，進(jìn)一步說(shuō)明了APF算法是可行的，且相比其他算法更加有效。

1 非線性混合效應(yīng)狀態(tài)空間模型

非線性狀態(tài)空間模型是一種常用的動(dòng)態(tài)系統(tǒng)模型，其主要有兩部分

其中，θ是固定效應(yīng);bj是隨機(jī)效應(yīng)，服從均值為0，協(xié)方差為D的正態(tài)分布。文中 θ被當(dāng)作未知參數(shù)，式(1)～式(3)被稱(chēng)為混合效應(yīng)狀態(tài)空間模型(NLMESSM)。這里研究的是單個(gè)目標(biāo)，即當(dāng)參數(shù)未知得到一個(gè)目標(biāo)的觀測(cè)信號(hào)值時(shí)，用序貫重要性采樣和改進(jìn)的輔助粒子濾波算法方法進(jìn)行在線估計(jì)狀態(tài)。

2 粒子濾波

序貫重要性采樣算法是一種通過(guò)蒙特卡洛模擬實(shí)現(xiàn)貝葉斯濾波器的技術(shù)，其核心思想是利用一系列隨機(jī)樣本的加權(quán)和來(lái)表示所需后驗(yàn)概率密度，從而得到狀態(tài)的估計(jì)值。

假設(shè)觀測(cè)值y1∶t已知，所求的狀態(tài)設(shè)為，狀態(tài)的后驗(yàn)密度為，則狀態(tài)估計(jì)的推導(dǎo)過(guò)程如下

則

在初始時(shí)刻假定隨機(jī)樣本{θj}mj=1是從參數(shù)的先驗(yàn)分布p(θ)中抽取的，則在下一時(shí)刻，從樣本中按照每個(gè)粒子對(duì)應(yīng)的權(quán)重進(jìn)行重采樣得到新樣本。以此類(lèi)推，假設(shè)在t時(shí)刻，得到一組參數(shù)樣本，這些樣本近似服從參數(shù)的后驗(yàn)分布p(θ具體更新步驟如下:

(1)初始時(shí)刻從先驗(yàn)分布 p(θ)中抽取樣本{θj}mj=1。

(3)按照權(quán)重wjt對(duì)樣本{θj}mj=1進(jìn)行重抽樣，得到t+1時(shí)刻樣本{θ*j}mj=1

因此，由式(5)可得

通過(guò)粒子濾波算法，用樣本均值來(lái)近似這一期望。定義ωt為粒子的權(quán)重，則

假設(shè)在t時(shí)刻有狀態(tài)變量的隨機(jī)樣本{xi}N，當(dāng)ti=1參數(shù)θj和觀測(cè)值yt+1為已知時(shí)，具體的更新步驟如下:(1)初始時(shí)從中得到樣本 x*it+1，j=1，2，…，N。(2)計(jì)算每個(gè)樣本對(duì)應(yīng)的權(quán)重 ωit+1 ∝。標(biāo)準(zhǔn)化權(quán)重。(3)重采樣則有

將式(7)代入式(6)可得

即得到了狀態(tài)的估計(jì)值x^t。

3 輔助粒子濾波

以上采樣重抽樣(SIS)算法雖簡(jiǎn)單易求，便于抽樣，但由于其僅從粒子運(yùn)動(dòng)和之前的一些狀態(tài)盲目的進(jìn)行抽樣，因而可能會(huì)使大量的低權(quán)值粒子丟失，最終導(dǎo)致更高的蒙特卡洛且濾波性能更差。一個(gè)可行的提高性能方法就是使用輔助粒子濾波(APF)。其以SIS為基礎(chǔ)，假設(shè)在t時(shí)刻，得到一組樣本{xj，θj∶j=1，…，t t N}，及相對(duì)應(yīng)的權(quán)重{ωjt∶j=1，…，N}。需注意的是，θ的下標(biāo)t表示參數(shù)來(lái)自t時(shí)刻的后驗(yàn)密度，而并非隨時(shí)間變化。有了t+1時(shí)刻的觀測(cè)值yt+1，需從p(xt+1，y1∶t+1)中產(chǎn)生樣本，則由貝葉斯理論可知

為了推導(dǎo)出參數(shù)和狀態(tài)的聯(lián)合濾波，明顯要處理關(guān)于未知的密度函數(shù)p(θ的問(wèn)題。

在t時(shí)刻，設(shè)已知后驗(yàn)參數(shù)樣本θjt和權(quán)重 ωjt，j=1，…，N，用一種離散的蒙特卡洛來(lái)近似p()，記和Vt分別為蒙特卡洛后驗(yàn)密度p(θ)的均值和方差矩陣。則平滑核密度形式為N

步驟1 用{μjt+1，mjt}來(lái)確定{xt，θ}的先驗(yàn)點(diǎn)估計(jì)，j=1，…，N，其中可由狀態(tài)方程計(jì)算得到，且

步驟2 按照概率權(quán)重gjt+1從集合{1，2，…，N}中抽取輔助整數(shù)變量，稱(chēng)為采樣指數(shù)k，其中g(shù)jt+1∝ωjtp

步驟3 從第k個(gè)核密度的正態(tài)分量中抽取新的參數(shù)向量 θkt+1，即 θkt+1～N(·mkt，h2Vt)。

步驟6 多次重復(fù)步驟2～步驟5產(chǎn)生最終后驗(yàn)近似所需的{xkt+1，θkt+1}及權(quán)重 ωkt+1。

通過(guò)以上得到的一組狀態(tài)樣本及其權(quán)重，對(duì)這些樣本進(jìn)行加權(quán)平均就是此時(shí)刻所估計(jì)的狀態(tài)，即

4 仿真

4.1 一維非線性追蹤模型

這兩個(gè)模型被稱(chēng)為狀態(tài)空間模型，其來(lái)自Andrade(1987)，是一類(lèi)非線性追蹤模型。這是一個(gè)關(guān)于目標(biāo)跟蹤的模型，假設(shè)研究對(duì)象為單個(gè)目標(biāo)，要對(duì)此目標(biāo)的狀態(tài)進(jìn)行估計(jì)。其中Vt～N(0，δ2v)，Wt～ N(0，δ2w)，定義時(shí)間長(zhǎng)度 T=100，α =0.5，β =25，γ =8，δ2v=10，δ2w=1，這里需固定狀態(tài)誤差δ2v和觀測(cè)誤差δ2w，取狀態(tài)方程中的參數(shù)α，β，γ作為未知參數(shù)，其中α，β，γ的分布為a ～N(0.5，0.03)，β ～ N(25，3)，γ ～ N(8，0.8)。令向量 θ=(α，β，γ)，則有混合效應(yīng)模型 θj=α +bj，j=1，2，…，m。其中，固定效應(yīng)θ=(0.5 25 8)，bj是隨機(jī)效應(yīng)，服從均值為0，協(xié)方差為D的正態(tài)分布。這里從初始狀態(tài)分布N(0.1，2)中分別抽取N個(gè)粒子，從參數(shù)θ的先驗(yàn)分布中抽取M個(gè)參數(shù)。

在表1中，針對(duì)不同的狀態(tài)和參數(shù)粒子數(shù)50，100，150，200，分別列出了SIR和APF在不同粒子個(gè)數(shù)時(shí)的均方根誤差(RMS)，其中N為粒子個(gè)數(shù)。估計(jì)值與真實(shí)值的均方根誤差為

表1 SIR與APF在粒子個(gè)數(shù)不同時(shí)各自的均方根誤差

4.2 仿真結(jié)果分析

由表1可看出，通過(guò)SIR濾波計(jì)算得到的RMS比APF所得到的RMS要大，其中在N=50時(shí)，兩種算法的RMS相差最大，隨著粒子數(shù)N的增多，差距逐漸減小，這說(shuō)明粒子數(shù)目的多少是影響仿真結(jié)果的一個(gè)因素，粒子數(shù)目越多，仿真效果越好。

在一維非線性追蹤模型中，兩種算法的時(shí)間長(zhǎng)度tf=100，粒子個(gè)數(shù)N=200，參數(shù)個(gè)數(shù)M=200。SIR濾波仿真結(jié)果如圖1所示，圖中紅線代表真實(shí)狀態(tài)，綠線代表參數(shù)未知時(shí)的狀態(tài)估計(jì)?？煽闯?，兩條曲線存在著差距，狀態(tài)估計(jì)的波動(dòng)范圍遠(yuǎn)小于真實(shí)值，這可能是由于參數(shù)的重要性密度函數(shù)選擇不合適，從而導(dǎo)致參數(shù)的樣本遠(yuǎn)離真實(shí)值。APF濾波仿真結(jié)果如圖2所示，從圖中可看出，APF算法逼近真值程度明顯更優(yōu)，且狀態(tài)估計(jì)結(jié)果與真實(shí)值之間相差較小，這是由于APF算法中對(duì)參數(shù)樣本按照各自的權(quán)重重新選擇參數(shù)的核密度，使得從核密度中產(chǎn)生的樣本更加接近真實(shí)值，從而使權(quán)值變化幅度更小，計(jì)算得到的RMS則也較小。

圖1 真實(shí)狀態(tài)及SIR估計(jì)狀態(tài)與時(shí)間t的關(guān)系

圖2 真實(shí)狀態(tài)及APF估計(jì)狀態(tài)與時(shí)間t的關(guān)系

5 結(jié)束語(yǔ)

基于混合效應(yīng)狀態(tài)空間模型，針對(duì)狀態(tài)估計(jì)提出了SIR算法和APF算法。這兩個(gè)算法各有利弊，如SIR濾波雖其先驗(yàn)分布函數(shù)簡(jiǎn)單易求，且易于抽樣，但由于與外部觀測(cè)值無(wú)關(guān)，因而產(chǎn)生的粒子大部分位于后驗(yàn)概率分布尾部，使得權(quán)值變化較大，仿真效果并不理想。而APF算法與外部觀測(cè)信息有關(guān)，使得采樣粒子更接近真實(shí)狀態(tài)，但由于過(guò)程相對(duì)復(fù)雜，因此實(shí)時(shí)性比SIR稍微差。整體而言，APF濾波仍比SIR濾波有大幅的改進(jìn)，故在處理此類(lèi)問(wèn)題中占有較大優(yōu)勢(shì)。

［1］趙昕東，耿鵬.基于Bayesian Gibbs Sampler的狀態(tài)空間模型估計(jì)方法研究及其在中國(guó)潛在產(chǎn)出估計(jì)上的應(yīng)用［J］.統(tǒng)計(jì)研究，2009，26(9):55 －63.

［2］錢(qián)偉民，柴根象.縱向數(shù)據(jù)混合效應(yīng)模型的統(tǒng)計(jì)分析［J］.數(shù)學(xué)年刊，2004，25(4):445 －456.

［3］包云霞，魯法明，劉福升.非線性狀態(tài)空間模型的EHMM抽樣法［J］.山東科技大學(xué)學(xué)報(bào)，2005，24(3):104 －106.

［4］侯代文，殷福亮.非線性系統(tǒng)中狀態(tài)和參數(shù)聯(lián)合估計(jì)的雙重粒子濾波方法［J］.電子與信息學(xué)報(bào)，2008，30(9):2128－2133.

［5］楊小軍，潘泉，王睿，等.粒子濾波進(jìn)展與展望［J］.控制理論與應(yīng)用，2006，23(2):261 －267

［6］劉維亭，戴曉強(qiáng)，朱志宇.基于重要性重采樣粒子濾波器的機(jī)動(dòng)目標(biāo)跟蹤方法［J］.江蘇科技大學(xué)學(xué)報(bào)，2007，21(1):37－41.

［7］江寶安，盧煥章.粒子濾波器及其在目標(biāo)跟蹤中的應(yīng)用［J］.雷達(dá)科學(xué)與技術(shù)，2003，1(3):170 －174.

［8］Sanjeev M，Arulam Palam，Simon Maskell，et al.A tutorial on particle filters for online nonlinear/non－gaussian bayesian tracking［J］.IEEE Transactions on Signal Processing，2002，50(2):174－188.

［9］李軍科，張串，吳建軍.基于蒙特卡洛方法的粒子濾波算法研究［J］.電腦與信息技術(shù)，2008，16(1):49 －53.

［10］李東風(fēng)，謝衷潔.Monte Carlo濾波新進(jìn)展及其應(yīng)用［J］.數(shù)學(xué)進(jìn)展，2004，33(4):415 －424.

［11］趙梅，張三同，朱剛.輔助粒子濾波算法及仿真舉例［J］.北京交通大學(xué)學(xué)報(bào)，2006，30(2):24 －28.

［12］王洪有.基于輔助粒子濾波算法的紅外目標(biāo)跟蹤［J］.應(yīng)用光學(xué)，2010，31(1):132 －135.