專利名稱：基于擬蒙特卡羅采樣的并行高斯粒子濾波方法
技術領域：
本發明屬于信號處理技術領域，涉及非線性濾波，特別是一種并行高斯粒子濾波方法，用于非線性動態系統的狀態估計。
背景技術：
非線性濾波技術廣泛應用于信號處理、自動化、計算機視覺和人工智能等許多領域。擴展卡爾曼濾波EKF是非線性濾波的經典方法。但由于這種擴展卡爾曼濾波是把非線性函數用它的二階泰勒展開式來代替，會產生近似誤差，精度較低，在應用上容易造成濾波發散的情況。無跡卡爾曼濾波UKF是非線性濾波的另一種方法，無跡卡爾曼濾波不用線性化方法逼近非線性函數，而是用有限個能完全描述狀態的均值和方差的樣本點來表示狀態的概率密度，通過直接應用狀態和量測方程來映射這些樣本點，加權求和得到更新的均值和方差。雖然無跡卡爾曼濾波比擴展卡爾曼濾波有更高的精度，但是由于無跡卡爾曼濾波和擴展卡爾曼濾波都是線性卡爾曼濾波的變形或改進形式，因此應用上也受線性卡爾曼濾波的制約，在對非高斯系統進行狀態估計時性能會變差。
成熟有效的粒子濾波是英國學者Gordon等人于1993年提出的，其基本思想是用一組帶權的隨機樣本表示狀態的后驗概率密度。粒子濾波沒有卡爾曼的限制，可以處理非線性和非高斯的系統，并且容易實現并行化。粒子濾波自提出以來，廣泛應用于信號處理、自動化、計算機視覺和人工智能等領域對非線性非高斯系統進行狀態估計，成為研究的熱點。該基本粒子濾波方法的不足之處是存在粒子退化，一般采用稱為“重采樣”的算法對其進行改進，但是這種重采樣不僅占用了很多計算機資源并且不易于并行實現。
Jayesh H.Kotecha等人于2003年提出高斯粒子濾波GPF，該高斯粒子濾波方法在高斯濾波的基礎上用重要性采樣方法計算更新高斯模型的均值和協方差矩陣兩個參數。GPF要優于傳統的高斯濾波器，較粒子濾波而言，GPF不需要重采樣步驟，更容易實現并行化。在系統狀態可以用單峰高斯模型描述的情況下，精度不在粒子濾波之下。
由于基本粒子濾波和高斯粒子濾波都是基于蒙特卡羅采樣，采樣粒子容易在狀態空間形成“空隙和簇”，導致算法估計性能下降。雖然通過增加粒子數可以減緩這種現象，但同時增加了算法的復雜度。高斯粒子濾波雖然避免了基本粒子濾波必須的重采樣，提高了計算效率，但采用蒙特卡羅采樣所帶來的估計性能下降問題仍然存在。
發明內容
本發明的目的在于克服上述已有技術的不足，提出一種基于擬蒙特卡羅QMC采樣的并行高斯粒子濾波方法，以在保證估計性能的條件下提高運算速度。
實現本發明技術目的的方案，包括以下步驟 并行擬蒙特卡羅序列產生步驟采用跳躍擬蒙特卡羅序列產生方法，并行產生后續濾波所需的P條擬隨機樣本序列； 并行高斯粒子濾波步驟通過P個執行單元，分別將每條擬隨機樣本序列轉換為服從指定高斯分布的擬高斯樣本序列，然后同時對非線性動態系統的狀態進行時間更新和狀態更新，最后對所有執行單元進行綜合處理，完成對非線性動態系統的濾波。
所述的并行產生后續濾波所需的P條擬隨機樣本序列，包括以下步驟 (1)產生偽隨機數x； (2)以x為初始值產生一個長度為P的串行擬蒙特卡羅序列sq[p]，p＝1，2，...，P； (3)分別以sq[1]為初始值產生跳躍擬蒙特卡羅序列Q1，以sq[2]為初始值產生跳躍擬蒙特卡羅序列Q2，以sq[p]為初始值產生跳躍擬蒙特卡羅序列QP； (4)將產生的P條擬隨機樣本序列分別輸入到P個執行單元。
所述的將每條擬隨機樣本序列轉換為服從指定高斯分布的擬高斯樣本序列，其步驟是首先采用逆采樣將所述的擬隨機樣本序列轉換為服從單位高斯分布的擬高斯樣本序列；再通過線性變換，把服從單位高斯分布的擬高斯樣本序列轉換為服從指定均值和協方差矩陣的擬高斯樣本序列。
所述的對非線性動態系統的狀態進行時間更新和狀態更新，其步驟是首先，每個執行單元將所述的服從指定高斯分布的擬高斯樣本序列代入非線性動態系統狀態方程，得到時間更新后的樣本；然后根據觀測數據，計算每個樣本的權值，分別計算每個執行單元內樣本的加權均值、加權協方差矩陣與權值和，得到更新后的狀態。
所述的對所有執行單元進行綜合處理，其步驟是首先計算所有執行單元中樣本的權值和，均值和與協方差矩陣和；然后將均值和與協方差矩陣和除以權值和得到非線性動態系統狀態估計結果。
本發明具有以下優點 1.由于本發明采用了擬蒙特卡羅QMC方法，產生了比蒙特卡羅MC隨機樣本更均勻的低偏差樣本點，因而在使用較少樣本點的情況下，提高了高斯粒子濾波的精度。
2.由于本發明采用了并行擬蒙特卡羅QMC序列產生方法，保證了各并行單元樣本的統計關系，提高了濾波性能的穩定性。
仿真結果表明，在樣本數相對較少時，本發明可以更均勻的逼近待估計狀態的后驗概率密度，從而比基本粒子濾波和高斯粒子濾波取得更優的濾波性能；而本發明采用并行處理很大的提高了濾波的實時性。
圖1是本發明的整體流程圖； 圖2是含有目標的模擬紅外灰度圖像； 圖3是模擬紅外圖像中目標的運動軌跡； 圖4是樣本數目為50時本發明與高斯粒子濾波、粒子濾波仿真誤差比較圖； 圖5是樣本數目為100時本發明與高斯粒子濾波、粒子濾波仿真誤差比較圖； 圖6是樣本數目為400時本發明與高斯粒子濾波、粒子濾波仿真誤差比較圖； 圖7是本發明執行單元個數與并行加速比的關系圖。
具體實施例方式 通過建立非線性系統動態模型來實施本發明提出的濾波方法，具體模型如下 狀態方程Xt＝f(Xt-1)+wt(1) 觀測方程Zt＝h(Xt)+et(2) 其中，f(·)，h(·)為有界的非線性函數，Xt為系統在t時刻的狀態，Zt為系統在t時刻的觀測值；wt為過程噪聲，et為觀測噪聲。
本發明濾波方法所涉及的參數包括 N為采樣數目，P＝2k為并行單元數目，p＝1，2，...，P 為并行單元的序號。
參照圖1，本發明方法包括并行擬蒙特卡羅序列產生步驟和并行高斯粒子濾波步驟。
一.并行擬蒙特卡羅序列產生步驟 采用跳躍擬蒙特卡羅序列產生方法，并行產生后續濾波所需的P條擬隨機樣本序列，具體步驟如下 1)產生偽隨機數x； 采用線性同余法，產生一個偽隨機數為 xt+1＝axtmod M，t＝0，1，...(3) 其中，a和M分別是乘子和模，t表示時間。
2)以x為初始值產生一個長度為P的串行擬蒙特卡羅序列sq[p]，p＝1，2，...，P； 本發明產生的串行擬蒙特卡羅序列是Sobol序列，該序列的產生需要一組“直接數”。在本發明實施之前，該“直接數”已計算完畢并存入固定寄存器，便于實時調用。該“直接數”的計算包括如下步驟 2a.選擇簡單多項式 Ps＝xd+a1xd-1+...+ad-1x+1，ai∈{0，1}(4) 2b.按照給定初始值迭代產生“直接數”vi 本發明中使用的Sobol序列可以通過當前的Sobol數sq[p]與選定的“直接數”v按照式(6)運算，迭代P次后得到串行擬蒙特卡羅序列sq[p]，p＝1，2，...，P 其中“直接數”vz[p]的下標z[p]是p的二進制形式的最右端0的位數，
表示按位異或運算。
3)分別以sq[1]為初始值產生跳躍擬蒙特卡羅序列Q1，以sq[2]為初始值產生跳躍擬蒙特卡羅序列Q2，以sq[p]為初始值產生跳躍擬蒙特卡羅序列QP； 以sq[1]為例，跳躍擬蒙特卡羅序列Q1的產生，包括以下步驟 3a.利用存儲器中的“直接數”生成“跳躍直接數” 其中，“直接數”和“跳躍直接數”的下標z[i]是i的二進制形式的最右端0的位數，i＝pP-1，pP-2，…，(p-1)P-1，
表示按位異或運算。
3b.將“跳躍直接數”存儲于固定的寄存器中，便于實時調用； 3c.以sq[1]作為初始值，與選定的“跳躍直接數”進行迭代N/P次運算，得到跳躍擬蒙特卡羅序列Q1 其中，Vz[iP-1]是事先計算好并存放在存儲器的跳躍直接數，z[iP-1]是iP-1的二進制形式的最右端0的位數，i＝1，2，…，N/P。
按照同樣的步驟，分別以sq[2]為初始值產生跳躍擬蒙特卡羅序列Q2，以sq[p]為初始值產生跳躍擬蒙特卡羅序列QP。
4)將產生的P條擬隨機樣本序列分別輸入到P個執行單元。
本發明中并行擬蒙特卡羅序列產生完畢，將得到的P條擬隨機樣本序列分別提供給后續操作的P個執行單元完成濾波。
二.并行高斯粒子濾波步驟 1)假設t時刻，第p個執行單元將擬隨機樣本序列轉換為服從均值為高斯分布的擬高斯樣本序列，具體如下 1a.將所述的第p條擬隨機樣本Qp，t-1通過逆采樣轉換為服從單位高斯分布的擬高斯樣本序列 Yp，t-1＝Gau-1(Qp，t-1)(9) 其中，Gau-1為單位高斯分布函數的逆函數； 1b.將擬高斯樣本序列Yp，t-1線性變換為服從前一時刻均值向量μt-1和協方差矩陣∑t-1的樣本序列Xp，t-1 Xp，t-1＝Rt-1Yp，t-1+μt-1(10) 其中， 2)將服從指定高斯分布的擬蒙特卡羅樣本序列Xp，t-1代入非線性系統模型的狀態方程式(1)，得到更新后的樣本序列 Xp，t＝f(Xp，t-1)+wt。(11) 3)根據當前時刻得到的觀測值Zt，計算每個更新后的樣本權值 Wp，t＝q(Zt|Xp，t)(12) 其中q(·)為似然函數，由觀測方程式(2)得到。
4)根據得到的新樣本序列Xp，t和樣本權值Wp，t，分別計算第p個執行單元中N/P個加權樣本的和
協方差矩陣
以及權值的和
為 該加權樣本的和
為更新后的狀態。
5)根據得到的每個執行單元中加權樣本的和
協方差矩陣
以及權值的和
p＝1，2，…，P，計算系統樣本權值總和SW、系統狀態均值向量μt與協方差矩陣∑n為 均值向量μt與協方差矩陣∑t即為t時刻濾波結果。
6)將協方差矩陣∑t進行Cholesky分解，得到一個上三角矩陣Rt和下三角矩陣
即Rt將用于t+1時刻擬高斯樣本的產生，完成t+1時刻的濾波。
非線性動態系統經過多次濾波后，得到的輸出結果將滿足性能要求，在新的觀測值不斷得到的情況下，濾波操作將一直持續。
本發明的效果通過仿真實例進一步說明 1、仿真條件 選擇一個非線性運動目標跟蹤的例子，模擬產生了一組包含一個弱小目標的紅外灰度圖像數據，共有80幀，大小為40×40，并加入仿真的弱小目標，如圖2所示。其中，圖2(a)是序列中第5幀模擬圖像，圖2(b)是序列中第25幀模擬圖像，目標肉眼無法識別。弱小目標從圖像序列的第1幀出現一直到第80幀，在二維平面內作非線性運動，其運動軌跡如圖3所示。目標動態方程為 xt+1＝Fxt+wt(15) 其中，t是圖像幀數，是狀態向量；(xt，yt)，
和It分別表示第t幀目標的位置，速度和強度；wt為協方差矩陣為w的零均值白噪聲，T為采樣周期，q1和q2分別表示目標運動噪聲和目標強度噪聲的功率譜密度。本仿真實例取T＝1，q1＝0.005，q2＝0.01，∑＝0.7，目標初始狀態給定為x0＝[70.470.318]T。
觀測模型
其中
表示目標信號對觀測的貢獻，
為噪聲 2.仿真內容及結果 1)根據仿真條件，以給定初始狀態作為本發明濾波方法的初始值，根據事先存儲的“直接數”生成濾波所需要的擬蒙特卡羅樣本；按照發明中所述的并行高斯濾波方法進行時間更新和觀測更新，最后計算得到目標的估計狀態。同時，為了驗證本發明的優點，將本發明與基本粒子濾波PF和高斯粒子濾波GPF進行性能對比，粒子數目分別為N＝50，100，400。進行100次蒙特卡羅仿真得到三種濾波方法在信噪比為2時的均方根誤差比較，如圖4、圖5和圖6所示，其中均方根誤差RMSE如下式定義 表1對比了在不同粒子數目情況下三種濾波方法的失跟率，實驗分別在粒子數為50，100和400的情況下進行了100次蒙特卡羅仿真。
表1 三種濾波方法的失跟率對比 從表1可以看出，在粒子數較少時，PF和GPF的失跟率都非常大，而本發明方法的失跟率則很小，而粒子數增加到400時，三者的失跟率都減小，但本發明方法依然優于PF和GPF。
從圖4、圖5和圖6可以看出，在粒子數目較少時，本發明進行估計的RMSE明顯小于PF和GPF，估計性能優于PF和GPF，而隨著粒子數的增加，三種濾波器的估計性能差別逐漸減小。
通過理論分析和實驗結果，可以看出在粒子數相對較少時，粒子的最優分布很重要，在同樣數目的粒子情況下，本發明方法可以更均勻的逼近待估計狀態的后驗概率密度，從而比PF和GPF取得更優的濾波性能。盡管三種濾波算法在粒子數很大時估計精度相近，但本發明方法用較少的粒子可以達到比較好的估計性能，因此在實際應用中會非常有用。
3)為了分析本發明并行處理的優勢，根據圖1所示方法流程，按照Amdahl定律對本發明的處理速度作了仿真分析。
Amdalh定律如下式 其中，SP表示并行處理與串行處理的加速比，s為串行處理工作量，m表示并行處理工作量，P為執行單元個數。
仿真中，串行處理的計算時間≤0.1％，根據Amdalh定律，隨著執行單元數目的增加，并行處理的加速比變化為圖7虛線所示。實際仿真中，在處理器個數分別為2、4、8、16、32、64的情況下得到加速比如圖7的實線所示，比串行處理占0.1％情況下的理論估計值還要好，可見發明中串行處理的計算時間比估計的0.1％小，因此發明采用并行系統處理能獲得很好的效率。
權利要求
1.一種基于擬蒙特卡羅采樣的并行高斯粒子濾波方法，包括
并行擬蒙特卡羅序列產生步驟采用跳躍擬蒙特卡羅序列產生方法，并行產生后續濾波所需的P條擬隨機樣本序列；
并行高斯粒子濾波步驟通過P個執行單元，分別將每條擬隨機樣本序列轉換為服從指定高斯分布的擬高斯樣本序列，然后同時對非線性動態系統的狀態進行時間更新和狀態更新，最后對所有執行單元進行綜合處理，完成對非線性動態系統的濾波。
2.根據權利1所述的并行高斯粒子濾波方法，其中所述的并行產生后續濾波所需的P條擬隨機樣本序列，包括以下步驟；
(1)產生偽隨機數x；
(2)以x為初始值產生一個長度為P的串行擬蒙特卡羅序列sq[p]，p＝1，2，...，P；
(3)分別以sq[1]為初始值產生跳躍擬蒙特卡羅序列Q1，以sq[2]為初始值產生跳躍擬蒙特卡羅序列Q2，以sq[p]為初始值產生跳躍擬蒙特卡羅序列QP；
(4)將產生的P條擬隨機樣本序列分別輸入到P個執行單元。
3.根據權利1所述的并行高斯粒子濾波方法，其中所述的將每條擬隨機樣本序列轉換為服從指定高斯分布的擬高斯樣本序列，首先采用逆采樣將所述的擬隨機樣本序列轉換為服從單位高斯分布的擬高斯樣本序列；再通過線性變換，把服從單位高斯分布的擬高斯樣本序列轉換為服從指定均值和協方差矩陣的擬高斯樣本序列。
4.根據權利1或3所述的并行高斯粒子濾波方法，其中所述的對非線性動態系統的狀態進行時間更新和狀態更新，包括以下步驟
(1)每個執行單元將所述的服從指定高斯分布的擬高斯樣本序列代入非線性動態系統狀態方程，得到時間更新后的樣本；
(2)根據觀測數據，計算每個樣本的權值，分別計算每個執行單元內樣本的加權均值、加權協方差矩陣與權值和，得到更新后的狀態。
5.根據權利1所述的并行高斯粒子濾波方法，其中所述的對所有執行單元進行綜合處理，包括計算所有執行單元中樣本的權值和，均值和與協方差矩陣和，并將均值和與協方差矩陣和除以權值和得到非線性動態系統狀態估計結果。
全文摘要
本發明公開了一種基于擬蒙特卡羅采樣的并行高斯粒子濾波方法，屬于信號處理技術領域，主要解決目前基于蒙特卡羅采樣的濾波方法在樣本數較少時對非線性動態系統狀態估計性能降低的問題。其濾波步驟包括并行擬蒙特卡羅序列產生步驟和并行高斯粒子濾波步驟，即首先采用跳躍擬蒙特卡羅序列產生方法，并行產生后續濾波所需的P條擬蒙特卡羅隨機樣本序列；然后分別通過P個執行單元，將每條擬蒙特卡羅隨機樣本序列轉換為服從指定高斯分布的擬高斯樣本序列；同時對非線性動態系統的狀態進行時間更新和狀態更新；最后對所有執行單元進行綜合處理，完成對非線性動態系統的濾波。本發明具有濾波性能高和實時性好的優點，可用于信號處理、自動化和人工智能領域。
文檔編號G06K9/00GK101436251SQ20081023276
公開日2009年5月20日 申請日期2008年12月22日 優先權日2008年12月22日
發明者姬紅兵, 斌 武, 曦 陳 申請人:西安電子科技大學