專利名稱：基于多建議分布的粒子濾波方法
技術(shù)領(lǐng)域：
本發(fā)明涉及基于多建議分布的粒子濾波方法，要求保護(hù)的技術(shù)方案屬于信號(hào)處理、人工智能和計(jì)算機(jī)視覺(jué)領(lǐng)域。
背景技術(shù)：
非線性濾波問(wèn)題在很多領(lǐng)域的問(wèn)題中都涉及到，其中包括信號(hào)處理、金融、人工智能和計(jì)算機(jī)視覺(jué)等等。解決非線性濾波問(wèn)題最為普遍的一種方法是使用擴(kuò)展卡爾曼濾波器(EKF)，但是EKF只適用于弱非線性系統(tǒng)，對(duì)于強(qiáng)非線性系統(tǒng)很容易導(dǎo)致發(fā)散。還有一種解決非線性濾波問(wèn)題的辦法是無(wú)跡卡爾曼濾波器(UKF)。UKF和EKF不同，它不使用局部線性化技術(shù)，而直接使用系統(tǒng)的非線性方程進(jìn)行運(yùn)算，從而能夠避免局部線性化所引入的誤差，避免在強(qiáng)非線性系統(tǒng)中出現(xiàn)發(fā)散。但是，UKF的使用不適用于一般的非高斯模型。
另一種解決非線性濾波問(wèn)題的方法就是粒子濾波器(Particle Filter)。它的基本思想就是使用一組帶有權(quán)值的粒子(樣本)集合來(lái)近似表示解決問(wèn)題時(shí)需要的系統(tǒng)后驗(yàn)概率密度。粒子濾波器自提出以來(lái)已經(jīng)廣泛應(yīng)用到信號(hào)處理、金融、人工智能、計(jì)算機(jī)視覺(jué)以及機(jī)器人等領(lǐng)域中，成為各國(guó)研究者關(guān)注的焦點(diǎn)。
粒子濾波算法中的一項(xiàng)關(guān)鍵技術(shù)是建議分布的選擇。選擇好的建議分布能夠有效的避免退化問(wèn)題的影響，提高系統(tǒng)的狀態(tài)估計(jì)精度。目前，常見(jiàn)的建議分布有先驗(yàn)概率分布(Transition Prior，TP)、EKF、UKF等。以這三種建議分布為基礎(chǔ)的粒子濾波器，分別被稱為基本粒子濾波器(Generic PF)、擴(kuò)展卡爾曼粒子濾波器(EKPF)和無(wú)跡粒子濾波器(UPF)。但是先驗(yàn)概率分布沒(méi)有考慮當(dāng)前時(shí)刻新的觀測(cè)值的影響，從而影響了粒子濾波器估計(jì)精度；而EKF采用局部線性化技術(shù)，引入了過(guò)多的截?cái)嗾`差，對(duì)估計(jì)精度產(chǎn)生影響；以UKF為建議分布的粒子濾波器具有較高的時(shí)間耗費(fèi)。

發(fā)明內(nèi)容
針對(duì)上述問(wèn)題，本發(fā)明提出一種基于多建議分布的粒子濾波算法。
該算法包括以下基本步驟一.初始狀態(tài)下，首先從初始分布p(x0)中提取一組(N個(gè))樣本粒子，并設(shè)置其初始的均值和方差。
二.對(duì)于c*N個(gè)粒子(1)首先使用UKF進(jìn)行粒子傳遞和更新，獲得關(guān)于該時(shí)刻狀態(tài)的估計(jì)xk ukf(i).
(2)使用EKF重復(fù)粒子更新過(guò)程，并以第(1)步中得到的狀態(tài)估計(jì)作為輸入。得到該時(shí)刻最終的估計(jì)量xk(i)和
(3)根據(jù)(2)中得到的估計(jì)量，構(gòu)成建議分布N(xk(i)， )，從中抽取新的樣本粒子，并賦予權(quán)值wki。
三.對(duì)于剩余的(1-c)*N個(gè)粒子，使用非線性系統(tǒng)方程進(jìn)行更新，得到狀態(tài)估計(jì)預(yù)測(cè)量xk|k-1(i)，并從先驗(yàn)概率分布p(xk|k-1(i)|xk-1(i))中抽取粒子，賦予權(quán)值wki。
四.對(duì)所有粒子的權(quán)值進(jìn)行歸一化處理。
五.執(zhí)行再采樣過(guò)程，對(duì)得到的新的粒子集合中的每一個(gè)粒子賦予相同的權(quán)值1/N。
六.輸出七.如繼續(xù)迭代，轉(zhuǎn)到二，否則結(jié)束。
本發(fā)明采用分治采樣策略，將系統(tǒng)所需粒子分成兩部分，分別由不同的建議分布中提取。
一部分(假設(shè)占百分比為c)由混合的建議分布提取，即由UKF和EKF組成的混合建議分布。其余的一部分(1-c)，由先驗(yàn)概率分布p(xk|xk-1)中提取。以此為基礎(chǔ)的粒子濾波算法不僅能夠產(chǎn)生較高的估計(jì)精度，而且能夠降低粒子濾波器的運(yùn)行時(shí)間。用戶可以根據(jù)自身對(duì)實(shí)際要求的不同，調(diào)整參數(shù)c，以滿足他們對(duì)運(yùn)行效率和精度的不同要求。
假設(shè)非線性動(dòng)態(tài)系統(tǒng)的狀態(tài)空間模型為xk＝fk(xk-1，vk-1)yk＝hk(xk，uk)其中xk表示系統(tǒng)在k時(shí)刻所處的狀態(tài)，yk表示k時(shí)刻的觀測(cè)值。函數(shù)f和h表示系統(tǒng)的狀態(tài)轉(zhuǎn)移和測(cè)量模型。vk和uk分別表示系統(tǒng)噪聲與測(cè)量噪聲。
本發(fā)明提供的多建議分布粒子濾波器采用了EKF、UKF和TP作為建議分布，將總的粒子數(shù)分成兩部分，分別由不同的建議分布中抽取，具體實(shí)現(xiàn)過(guò)程包括以下步驟第一步.k＝0時(shí)，初始化一組粒子作為初始狀態(tài)下的粒子集合，并設(shè)置初始的均值和協(xié)方差矩陣，設(shè)置方法同UPF算法中的初始設(shè)置。
第二步.對(duì)于c*N個(gè)粒子，1.首先使用UKF更新粒子(1)選擇一組sigma點(diǎn)，構(gòu)成sigma矩陣。
χk-1(i)a=x‾t-1(i)ax‾t-1(i)a±(na+λ)Pt-1(i)a]]>(2)根據(jù)動(dòng)態(tài)系統(tǒng)的非線性方程，將這些sigma點(diǎn)向前傳遞，得到系統(tǒng)狀態(tài)的均值的預(yù)測(cè)量xk|k-1ukf(i)和協(xié)方差的預(yù)測(cè)量Pk|k-1ukf(i)，以及與測(cè)量值有關(guān)的預(yù)測(cè)量yk|k-1ukf(i)(預(yù)測(cè))
(3)獲取當(dāng)前時(shí)刻的測(cè)量值yk之后，并根據(jù)該測(cè)量值對(duì)(2)中得到的一步預(yù)測(cè)值(更新)x‾k(i)ukf=x‾k|k-1(i)ukf+Kk(yk-y‾k|k-1(i)ukf)]]>其中Kk=PxkykPy~ky~k-1]]>(卡爾曼增益)Py~ky~k=Σj=0na2Wj(c)[Yj,k|k-1(i)-y‾k|k-1(i)][Yj,k|k-1(i)-y‾k|k-1(i)]T]]>Px~ky~k=Σj=0na2Wj(c)[χj,k|k-1(i)-x‾k|k-1(i)][Yj,k|k-1(i)-y‾k|k-1(i)]T]]>這樣就獲得了估計(jì)量xk ukf(i)。
其中，xa＝[xTvTuT]T為擴(kuò)張(augmented)狀態(tài)變量，由系統(tǒng)狀態(tài)x，系統(tǒng)噪聲v和測(cè)量噪聲u組成。χa＝[(χx)T(χv)T(χu)T]T為sigma矩陣。na＝nx+nv+nu為xa的維數(shù)。χk-1(i)a表示k-1時(shí)刻由sigma點(diǎn)構(gòu)成的sigma矩陣，χj，k|k-1(i)為χk|k-1(i)的第j個(gè)元素。Wi(m)和Wi(c)表示sigma向量的權(quán)值，分別用來(lái)求均值和協(xié)方差估計(jì)，其計(jì)算方法在文獻(xiàn)(Rudolph van der Merwe等,The Unscented Particle Filter)中。λ＝α2(nx+k)-nx是一個(gè)尺度調(diào)節(jié)因子，α決定了選擇的sigma點(diǎn)在其均值x附近的分布情況，通常將α設(shè)置為一個(gè)很小的正值(如0.001)。k是次級(jí)尺度調(diào)節(jié)因子，通常設(shè)置為0，β是用來(lái)結(jié)合關(guān)于x的分布的先驗(yàn)知識(shí)(對(duì)于高斯分布，β的最佳取值為2)。 是矩陣(nx+λ)Px平方根的第i行，2.使用EKF更新粒子，并以xk ukf(i)作為輸入，即令xk-1(i)=x‾k(i)ukf.]]>(1)傳遞粒子，并得到相關(guān)的一步預(yù)測(cè)量
x‾k|k-1(i)ekff(x‾k-1(i))=f(x‾k(i)ukf)]]>Pk|k-1(i)ekf=Fk(i)Pk-1(i)FkT(i)+Gk(i)QkGkT(i)]]>(2)計(jì)算系統(tǒng)狀態(tài)轉(zhuǎn)移模型的雅可比矩陣F和G，以及測(cè)量模型的雅可比矩陣H和U.并計(jì)算出卡爾曼增益。
Kk=Pk|k-1(i)ekfHkT(i)[Uk(i)RkUkT(i)+Hk(i)Pk|k-1(i)ekfHkT(i)]-1]]>(3)根據(jù)當(dāng)前的時(shí)刻的觀測(cè)值yk，修正(1)中的一步預(yù)測(cè)量，得到修正的的估計(jì)量P^k(i)ekf=Pk|k-1(i)ekf-KkHk(i)Pk|k-1(i)ekf]]>x‾k(i)ekf=x‾k|k-1(i)ekf+P^k(i)ekfHkT(i)Rk-1(yk-h(x‾k|k-1(i)ekf))]]>(4)令x‾k(i)=x‾k(i)ekf,P^k(i)=P^k(i)ekf]]>3.從它們所組成的建議分布中抽取新的粒子x^k(i)~q(xk(i)|x0:k-1(i),y1:k)=N(x‾k(i),P^k(i))]]>并賦予其權(quán)值wki。
其中，Q和R分別表示系統(tǒng)噪聲和測(cè)量噪聲的方差第三步.對(duì)于剩余的(1-c)*N個(gè)粒子，使用先驗(yàn)概率分布進(jìn)行更新。
1.直接根據(jù)系統(tǒng)的非線性方程計(jì)算出狀態(tài)估計(jì)量x‾k|k-1(i)=f(xk-1(i))]]>2.從先驗(yàn)分布密度函數(shù)p(xk|k-1(i)|xk-1(i))中抽取粒子，即 并為其賦予權(quán)值wki。
第四步.對(duì)所有的粒子的權(quán)值進(jìn)行歸一化處理wki=wki/Σwkj]]>第五步.執(zhí)行再采樣過(guò)程去除那些權(quán)值較小的粒子，復(fù)制權(quán)值較大的粒子，得到一組近似服從后驗(yàn)概率分布p(x0:k|z1:k)的粒子集合，并為其中的每個(gè)粒子賦予相同的權(quán)值1/N。
第六步.輸出階段。
第七步.k＝k+1，如繼續(xù)迭代，轉(zhuǎn)到第二步循環(huán)，否則結(jié)束。
本發(fā)明具有以下優(yōu)點(diǎn)1.采用多建議分布方案，融合了各種獨(dú)立建議分布在粒子濾波器框架內(nèi)的優(yōu)點(diǎn)。
2.采用分治采樣策略，在粒子濾波器的運(yùn)行時(shí)間和估計(jì)精度等方面作折衷，用戶可以根據(jù)自己的需求適當(dāng)調(diào)整參數(shù)，獲得對(duì)運(yùn)行時(shí)間或精度的要求。


圖1-多建議分布的粒子濾波流程圖；圖2-在不同的c取值下，多建議分布粒子濾波器算法與其它粒子濾波器算法的估計(jì)均方誤差對(duì)比曲線；圖3-在不同的c取值下，多建議分布粒子濾波器算法與其它粒子濾波器算法的運(yùn)行時(shí)間對(duì)比曲線；圖4-原始信號(hào)曲線；圖5-使用本發(fā)明進(jìn)行濾波后的曲線；具體實(shí)施方式
下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明作進(jìn)一步說(shuō)明。
本發(fā)明的算法包括以下基本步驟一.初始狀態(tài)下，首先從初始分布p(x0)中提取一組(N個(gè))樣本粒子，并設(shè)置其初始的均值和方差。
二.對(duì)于c*N個(gè)粒子(1)首先使用UKF進(jìn)行粒子傳遞和更新，獲得關(guān)于該時(shí)刻狀態(tài)的估計(jì)xk ukf(i)。
(2)使用EKF重復(fù)粒子更新過(guò)程，并以第(1)步中得到的狀態(tài)估計(jì)作為輸入。得到該時(shí)刻最終的估計(jì)量xk(i)和 (3)根據(jù)(2)中得到的估計(jì)量，構(gòu)成建議分布N(xk(i)， )，從中抽取新的樣本粒子，并賦予權(quán)值wki。
三.對(duì)于剩余的(1-c)*N個(gè)粒子，使用非線性系統(tǒng)方程進(jìn)行更新，得到狀態(tài)估計(jì)預(yù)測(cè)量xk|k-1(i)，并從先驗(yàn)概率分布p(xk|k-1(i)|xk-1(i))中抽取粒子，賦予權(quán)值wki。
四.對(duì)所有粒子的權(quán)值進(jìn)行歸一化處理。
五.執(zhí)行再采樣過(guò)程，對(duì)得到的新的粒子集合中的每一個(gè)粒子賦予相同的權(quán)值1/N。
六.輸出七.如繼續(xù)迭代，轉(zhuǎn)到二，否則結(jié)束。
實(shí)施例1.通過(guò)一個(gè)非線性動(dòng)態(tài)系統(tǒng)的狀態(tài)空間模型如下xk＝1+sin(0.04π(k-1))+0.5xk-1+vk-1(狀態(tài)轉(zhuǎn)移模型)
zk=0.2xk2+ukk≤300.5xk-2+ukk>30]]>(測(cè)量模型)對(duì)本發(fā)明與其它幾種粒子濾波器算法進(jìn)行對(duì)比。使用粒子數(shù)目N＝200，觀測(cè)時(shí)間T＝60。程序在CPU為2.67GHz，內(nèi)存為1GB的微型計(jì)算機(jī)上運(yùn)行100次。
每次運(yùn)行程序的輸出為粒子集合的均值x^k=1NΣj=1Nxkj.]]>每一次運(yùn)行中，均方誤差的計(jì)算公式為MSE=(1TΣk=1T(x^k-xk)2)1/2]]>取c＝0.3，在時(shí)間耗費(fèi)和精度方面進(jìn)行對(duì)比時(shí)，本發(fā)明算法產(chǎn)生的估計(jì)均方誤差比UPF低80％，運(yùn)行時(shí)間比UPF低約50％。在c＝0.4時(shí)，本發(fā)明算法的均方誤差比UPF低72％，運(yùn)行時(shí)間比UPF低44％。
圖2所示，c在取不同的值時(shí)，本發(fā)明與其它粒子濾波器算法均方誤差的對(duì)比，圖3為相應(yīng)的運(yùn)行時(shí)間對(duì)比。隨著c值的增加，本發(fā)明提出的算法在運(yùn)行時(shí)間上也隨之增加，但是本發(fā)明所產(chǎn)生的均方誤差總體上是呈現(xiàn)遞減的趨勢(shì)。
實(shí)施例2.將本發(fā)明實(shí)施于鉆井系統(tǒng)中的隨鉆測(cè)井儀。隨鉆測(cè)井儀器分為井上部分和井下部分。井下部分的測(cè)量信號(hào)經(jīng)泥漿傳送到地面。井上部分首先要對(duì)接收到的井下信號(hào)進(jìn)行濾波，然后是解碼并算出物理參數(shù)值。本發(fā)明的作用就是將帶噪聲的編碼信號(hào)進(jìn)行濾波，得到去噪的信號(hào)。將本發(fā)明實(shí)施于儀器的井上部分。本發(fā)明能夠得到理想的濾波效果(圖4，圖5)，為之后數(shù)據(jù)的處理提供了基礎(chǔ)。
權(quán)利要求
1.基于多建議分布的粒子濾波方法，其特征在于包括以下步驟1)首先初始化一組(N個(gè))樣本粒子，并設(shè)置其初始的均值和方差；2)對(duì)于c*N個(gè)粒子(1)首先使用UKF進(jìn)行粒子傳遞和更新，獲得關(guān)于該時(shí)刻狀態(tài)的估計(jì)xk(i)ukf；(2)使用EKF重復(fù)粒子更新過(guò)程，并以第(1)步中得到的狀態(tài)估計(jì)作為輸入。得到該時(shí)刻最終的估計(jì)量xk(i)和 (3)根據(jù)(2)中得到的估計(jì)量，構(gòu)成建議分布N(xk(i)， )，并從中抽取新的樣本粒子，并賦予權(quán)值wki；3)對(duì)于剩余的(1-c)*N個(gè)粒子，使用非線性系統(tǒng)方程進(jìn)行更新，得到狀態(tài)估計(jì)預(yù)測(cè)量xk|k-1(i)，并從先驗(yàn)概率分布p(xk|k-1(i)|xk-1(i))中抽取粒子，賦予權(quán)值wki；4)對(duì)所有粒子的權(quán)值進(jìn)行歸一化處理；5)執(zhí)行再采樣過(guò)程，對(duì)得到的新的粒子集合中的每一個(gè)粒子賦予相同的權(quán)值1/N；6)輸出；7)如繼續(xù)迭代，轉(zhuǎn)到二，否則結(jié)束。
全文摘要
本發(fā)明涉及基于多建議分布的粒子濾波方法，屬于信號(hào)處理、人工智能和計(jì)算機(jī)視覺(jué)領(lǐng)域。本發(fā)明在粒子濾波器的框架內(nèi)使用多種建議分布先驗(yàn)概率分布、擴(kuò)展卡爾曼濾波器、unscented卡爾曼濾波器等，采用分治采樣策略來(lái)處理樣本粒子，把總的粒子數(shù)分成多個(gè)部分，分別從不同的建議分布中提取，這樣可以降低粒子濾波器的運(yùn)行時(shí)間，提高運(yùn)行效率，而不損失粒子濾波器的估計(jì)精度。用戶可以根據(jù)自己對(duì)時(shí)間和精度的需求進(jìn)行參數(shù)設(shè)置。在涉及非線性濾波問(wèn)題的領(lǐng)域具有廣泛的應(yīng)用前景。
文檔編號(hào)G06F17/00GK101055563SQ200710099440
公開(kāi)日2007年10月17日 申請(qǐng)日期2007年5月21日 優(yōu)先權(quán)日2007年5月21日
發(fā)明者趙清杰, 王法勝, 王巍 申請(qǐng)人:北京理工大學(xué)