粒子濾波器

粒子濾波器（particle filter）是一種使用蒙特卡羅方法的遞歸濾波器。它透過一組具有權重的隨機樣本（粒子）表示隨機事件的後驗機率，並從含有雜訊或不完整的觀測序列中，估計動態系統的狀態。粒子濾波器可以運用在任何狀態空間的模型上。

同時，粒子濾波器是卡爾曼濾波器的一般化方法。卡爾曼濾波器建立在線性的狀態空間和高斯分布的雜訊上；而粒子濾波器的狀態空間模型可以是非線性，且雜訊分布可以是任何型式。

## 歷史

### 啟發式演演算法

從統計和機率的角度來看，粒子濾波器屬於遺傳算法分支過程的類別，並且是指平均場類型相互作用的粒子方法。這些粒子方法的解釋取決於科學專業。在進化計算中，平均場遺傳類型粒子方法通常用作啟發式和自然搜索算法（也稱為元啟發算法）。在計算物理學和分子化學中，它們用於解決費曼-卡茲(Feynman-Kac)路徑積分問題，或者用於計算玻爾茲曼-吉布斯(Boltzmann-Gibbs)度量，薛丁格算子的最高特徵值和基態。在生物學與遺傳學中，它們還代表了某些環境中一群個體或基因的進化。

平均場類型進化計算技術的起源於1950年艾倫·圖靈在遺傳類型突變選擇學習機上的開創性工作，以及1954年紐澤西州普林斯頓高級研究所的尼爾斯·艾爾·巴里切利(Nils Aall Barricelli)的文章。統計學中的粒子濾波器最早可追溯到1950年代中期。哈默斯利(Hammersley)等人在1954年提出的「窮人的蒙特卡洛」法暗示了當今使用的遺傳類型粒子過濾方法。1963年，尼爾斯·艾爾·巴里切利(Nils Aall Barricelli)模擬了一種遺傳類型演算法，以模仿個人玩簡單遊戲的能力。在進化計算文獻中，遺傳類型突變選擇算法通過1970年代初期約翰·霍蘭德(John Holland)的開創性工作而流行，特別是1975年出版的著作。

在生物學與遺傳學中，澳大利亞遺傳學家亞歷克斯·弗雷澤（Alex Fraser）在1957年還發表了一系列有關人工選擇生物的遺傳類型模擬的論文。由生物學家進行的進化計算機模擬在1960年代初變得更加普遍，並且該方法在弗雷澤和伯內爾（Burnell）（1970）和克羅斯比（Crosby）（1973）的書中進行了描述。弗雷澤的模擬包括了現代突變選擇遺傳粒子演算法的所有基本要素。

從數學觀點來看，給定一些部分和雜訊觀測結果的信號隨機狀態的條件分佈，是由信號的隨機軌跡上的費曼-卡茲（Feynman-Kac）機率描述的，該機率由一系列似然勢函數加權。量子蒙特卡洛法，更具體地說是擴散蒙特卡洛法，也可以解釋為費曼-卡茲路徑積分的平均場遺傳類型粒子近似。量子蒙特卡羅方法的起源通常歸因於恩里科·費米(Enrico Fermi)和羅伯特·里克特邁耶（Robert Richtmyer），他們於1948年開發了一種中子鏈反應的均值場粒子解釋方法，但它是第一個類似啟發式和遺傳類型的粒子演算法（又名「重採樣」或「重新配置」蒙特卡洛方法）是在1984年由Jack H. Hetherington提出的用於估計量子系統（在簡化的矩陣模型中）的基態能量的方法。還可以引用1951年出版的西奧多·愛德華·哈里斯（Theodore E. Harris）和赫爾曼·康恩（Herman Kahn）的開創性著作，該著作使用平均場，但類似於啟發式遺傳方法，用於估計粒子傳輸能量。在分子化學中，遺傳啟發式粒子方法（又名修剪和富集策略）的使用可以追溯到1955年馬歇爾·N·羅森布魯斯（Marshall N. Rosenbluth）和阿里安娜·W·羅森布盧斯（Arianna W. Rosenbluth）的開創性工作。

遺傳粒子演算法在高級信號處理和貝葉斯推斷中的使用是最近的。1993年1月，北川源四郎開發了「蒙特卡羅濾波器」，該文在1996年出現了輕微修改。1993年4月，戈登等人在他們的開創性著作中發表了遺傳類型演算法在貝葉斯統計推斷中的應用。作者將其演算法命名為「自舉濾波器」，並證明與其他濾波方法相比，其自舉演算法不需要對該狀態空間或系統噪聲進行任何假設。獨立的是皮埃爾·德爾·莫拉爾(Pierre Del Moral)和希米爾孔·卡瓦略(Himilcon Carvalho)，皮埃爾·德爾·莫拉爾(Pierre Del Moral)，安德烈·莫寧(André Monin)和熱拉爾·薩魯特(Gérard Salut)在1990年代中期發表的有關粒子過濾器的論文。1989-1992年初，皮埃爾·德爾·莫拉爾(P.Del Moral)，J.C.諾耶(J.C. Noyer)，G.Rigal和G.Salut在LAAS-CNRS中也開發了信號處理技術，並通過STCAN（服務技術）進行了一系列受限和分類研究des Constructions and Armes Navales），IT公司DIGILOG和LAAS-CNRS（系統分析和體系結構實驗室）來解決RADAR/SONAR和GPS信號處理問題。

### 數學基礎

從1950年到1996年，所有關於粒子過濾器、遺傳算法的出版物，包括在計算物理和分子化學中引入的修剪和重採樣蒙特卡洛方法，都提出了適用於不同情況下的類似於自然和啟發式的演算法，而沒有任何證據證明它們的一致性，也沒有討論估計的偏差以及基於族譜和祖先樹的演算法。

這些粒子演算法的數學基礎和首次嚴格的分析是由皮埃爾·德爾·莫拉爾（Pierre Del Moral）在1996年提出的。他的著作還包含了似然函數和未標準化的條件機率測度的粒子近似的無偏性質的證明。文章介紹的似然函數的無偏粒子估計器如今已用於貝式統計推理中。

直到1990年代末，Dan Crisan、Jessica Gaines、Terry Lyons、Pierre Del Moral等人提出了具有不同種群大小的分支類型粒子方法。P. Del Moral，A。Guionnet和L. Miclo在2000年開發了該領域的進一步發展。第一個中心極限定理是由Pierre Del Moral和Alice Guionnet於1999年以及Pierre Del Moral和Laurent Miclo於2000年提出的。Pierre於1990年代末提出了關於粒子過濾器時間參數的第一個統一收斂結果。Del Moral和Alice Guionnet。2001年，P。Del Moral和L. Miclo首次對基於族譜樹的粒子濾波器平滑器進行了嚴格的分析。

關於Feynman-Kac粒子方法論的理論和相關的粒子過濾器演算法已在2000年和2004年的書中得到發展。這些抽象的機率模型封裝了遺傳類型演算法，粒子和自舉濾波器，交互的卡爾曼濾波器（又稱Rao–Blackwellized粒子濾波器），重要性採樣和重採樣樣式的粒子濾波器技術，包括基於族譜樹的方法和用於解決濾波和平滑問題的粒子後向方法。其他類別的粒子濾波方法包括基於族譜樹的模型、後向馬爾可夫粒子模型、自適應平均場粒子模型、島型粒子模型和粒子馬爾可夫鏈蒙特卡羅方法。

## 濾波問題

### 目標

粒子濾波器的目標是通過觀測值估計狀態變量的後驗概率。粒子濾波器是針對隱馬爾可夫模型設計的，系統中既包含隱藏的變量，也包含可觀測的變量。可觀測變量（觀測過程）與隱藏變量（狀態過程）通過某種已知的函數關係相關聯。類似地，定義狀態變量演化的動態系統的機率描述是已知的。

一個通用粒子濾波器通過觀察測量過程估計隱藏狀態的後驗分布。考察下圖所示的狀態空間：

$$\begin{array}{cccccccccc} X_0&\to &X_1&\to &X_2&\to&X_3&\to &\cdots&\text{signal}\\ \downarrow&&\downarrow&&\downarrow&&\downarrow&&\cdots&\\ Y_0&&Y_1&&Y_2&&Y_3&&\cdots&\text{observation} \end{array}$$

濾波問題是要通過給定的觀測過程 $Y_0,\cdots,Y_k,$ 在任意時刻 k 的值，依次估計隱藏狀態 $X_k$ 的值。

所有的對 $X_k$ 的貝葉斯估計服從後驗分布 $p(x_k\mid y_0,y_1,\dots,y_k)$。粒子濾波器方法使用經驗測量和通用粒子演算法，提供了這些條件概率的近似值。與之對比，馬爾可夫鏈蒙特卡洛（MCMC）或重要性採樣方式會對整個後驗概率$p(x_0,x_1,\dots,x_k\mid y_0,y_1,\dots,y_k)$進行建模。

### 信號-觀測模型

粒子濾波器能夠從一系列含有雜訊或不完整的觀測值中，估計出動態系統的內部狀態。在動態系統的分析中，需要兩個模型，一個用來描述狀態隨時間的變化（系統模型），另一個用來描述每個狀態下觀測到的雜訊（觀測模型），將這兩個模型都用機率來表示。

在許多情況下，每得到一個新的觀測值時，都必須對系統做出一次估計，利用遞迴濾波器，能夠有效地達到這樣的目的。遞迴濾波器會對得到的資料做連續處理，而非分批處理，因此不需要將完整的資料儲存起來，也不需要在得到新的觀測值時，將現有的資料重新做處理。遞迴濾波器包含兩個步驟：

預測：利用系統模型，由前一個狀態的資訊預測下一個狀態的機率密度函數。

更新：利用最新的觀測值，修改預測出的機率密度函數。

藉由貝葉斯推斷（Baysian inference），我們可以描述出狀態空間的機率，並在得到新的觀測值時，對系統做出更新，因而達成上述目的。

### 近似貝式計算模型

在某些問題中，在信號的隨機狀態下，觀測值的條件分佈可能不具有密度，或者不可能或太複雜而無法計算。在這種狀況下，我們需要求近似值。其中一個策略是藉由馬爾可夫鏈$\mathcal X_k=\left(X_k,Y_k\right)$來取代信號 $X_k$，採用虛擬觀測的形式

$$\mathcal Y_k=Y_k+\epsilon \mathcal V_k\quad\mbox{for some parameter}\quad\epsilon\in [0,1]$$

對於具有已知機率密度函數的獨立序列中的某些序列，核心的概念是觀測下式

$$\text{Law}\left(X_k|\mathcal Y_0=y_0,\cdots, \mathcal Y_k=y_k\right)\approx_{\epsilon\downarrow 0} \text{Law}\left(X_k|Y_0=y_0,\cdots, Y_k=y_k\right)$$

在部分的觀測 $\mathcal Y_0=y_0,\cdots, \mathcal Y_k=y_k,$ 之下，粒子濾波器是和馬爾可夫鏈 $\mathcal X_k=\left(X_k,Y_k\right)$ 有相關的，定義為根據演化出的 $\mathbb R^{d_x+d_y}$，在一些明顯的濫用性符號 $p(\mathcal Y_k|\mathcal X_k)$ 下，帶有似然函數的粒子。這些機率性的技術跟近似貝式計算非常相關。在量子濾波器中，這些近似貝式量子濾波器的技術在1998年被皮埃爾·德爾·莫拉爾（P. Del Moral）, 讓·雅各德(J. Jacod)和菲力普·普羅特(P. Protter)所採用。且被皮埃爾·德爾·莫拉爾，阿諾·杜斯(A. Doucet)和阿傑·賈斯拉(A. Jasra)更進一步發展。

### 非線性濾波方程

藉由貝式定理在條件機率中可以得出

$$p(x_0, \cdots, x_k|y_0,\cdots,y_k) =\frac{p(y_0,\cdots,y_k|x_0, \cdots, x_k) p(x_0,\cdots,x_k)}{p(y_0,\cdots,y_k)}$$

當

$$\begin{align} p(y_0,\cdots,y_k) &=\int p(y_0,\cdots,y_k|x_0,\cdots, x_k) p(x_0,\cdots,x_k) dx_0\cdots dx_k \\ p(y_0,\cdots, y_k|x_0,\cdots ,x_k) &=\prod_{l=0}^{k} p(y_l|x_l) \\ p(x_0,\cdots, x_k) &=p_0(x_0)\prod_{l=1}^{k} p(x_l|x_{l-1}) \end{align}$$

粒子濾波器也是一個近似值，當有足夠的粒子時，結果會更加的準確。非線性濾波方程可以藉由遞迴得出

$$\begin{align} p(x_k|y_0,\cdots,y_{k-1}) &\stackrel{\text{updating}}{\longrightarrow} p(x_k|y_0,\cdots,y_k)=\frac{p(y_k|x_k)p(x_k|y_0,\cdots,y_{k-1})}{\int p(y_k|x'_k)p(x'_k|y_0,\cdots,y_{k-1})dx'_k} \\ &\stackrel{\text{prediction}}{\longrightarrow}p(x_{k+1}|y_0,\cdots,y_k)=\int p(x_{k+1}|x_k) p(x_k|y_0,\cdots,y_k) dx_k \end{align}$$

依照慣例 $p(x_0|y_0,\cdots,y_{k-1})=p(x_0)$ 在k=0時。非線性濾波方程問題的關鍵在於依序計算這些條件機率分布。

### 費曼-卡茲公式

固定一個時間範圍和一個序列的觀測 $Y_0=y_0,\cdots,Y_n=y_n$，在k=0....n，我們定義:

$$G_k(x_k)=p(y_k|x_k).$$

在這種表示法中，對於從原點k = 0到時間k = n的 軌跡集合上的任何有界函數F，我們都可以用費曼-卡茲公式來表示

$$\begin{align} \int F(x_0,\cdots,x_n) p(x_0,\cdots,x_n|y_0,\cdots,y_n) dx_0\cdots dx_n &= \frac{\int F(x_0,\cdots,x_n) \left\{\prod\limits_{k=0}^{n} p(y_k|x_k)\right\}p(x_0,\cdots,x_n) dx_0\cdots dx_n}{\int \left\{\prod\limits_{k=0}^{n} p(y_k|x_k)\right\}p(x_0,\cdots,x_n) dx_0\cdots dx_n}\\ &=\frac{E\left(F(X_0,\cdots,X_n)\prod\limits_{k=0}^{n} G_k(X_k)\right)}{E\left(\prod\limits_{k=0}^{n} G_k(X_k)\right)} \end{align}$$

這些費曼-卡茲路徑積分模型出現在各種科學學科中，包括計算物理、生物學、資訊理論和電腦科學。他們的解釋取決於應用領域。舉例來說，如果我們選擇狀態空間某些子集的指標函數$G_n(x_n)=1_A(x_n)$，它們代表馬爾可夫鏈在給定管中的條件分佈。我們會得到

$$E\left(F(X_0,\cdots,X_n) | X_0\in A, \cdots, X_n\in A\right) =\frac{E\left(F(X_0,\cdots,X_n)\prod\limits_{k=0}^{n} G_k(X_k)\right)}{E\left(\prod\limits_{k=0}^{n} G_k(X_k)\right)}$$

和

$$P\left(X_0\in A,\cdots, X_n\in A\right)=E\left(\prod\limits_{k=0}^{n} G_k(X_k)\right)$$

只要標準化常數嚴格為正。

## 非線性貝葉斯追蹤

在動態系統中，我們常常需要對物體做追蹤。假設有一動態系統，已知其狀態序列$\{x_k, k\in N\}$定義為

$$x_k=f_k(x_{k-1},v_{k-1})$$

其中$f_k:R^{n_x}\times R^{n_v}\rightarrow R^{n_x}$是狀態轉移函數，可以是非線性的函數，$\{v_{k-1},k\in N\}$是一個獨立且相同分布（i.i.d.）的過程雜訊序列，$n_x$和$n_v$分別代表狀態和過程雜訊向量的維度，$N$為自然數的集合。追蹤的目標是要遞迴地從觀測值$y_k$估計出$x_k$，而觀測值$y_k$定義為

$$y_k=h_k(x_k,n_k)$$

其中$h_k:R^{n_x}\times R^{n_n}\rightarrow R^{n_y}$是觀測函數，可以是非線性的函數，$\{n_k,k\in N\}$是一個獨立且相同分布（i.i.d.）的觀測雜訊序列，$n_y$和$n_n$分別代表觀測值和觀測雜訊向量的維度。

從貝葉斯理論的觀點來看，追蹤問題是要根據到時間$k$為止的已知資訊$y_{1:k}$，遞迴地計算出時間$k$的狀態$x_k$的可信度，這個可信度用機率密度函數$p(x_k\mid y_{1:k})$來表示。假設初始機率密度函數$p(x_0\mid y_0)\equiv p(x_0)$(稱為先驗機率)為已知，則利用「預測」和「更新」兩個步驟就能遞迴地計算出機率密度函數$p(x_k\mid y_{1:k})$。

在此假設在時間$k-1$的機率密度函數$p(x_{k-1}\mid y_{1:k-1})$為已知。

### 預測

利用查普曼-科爾莫戈羅夫等式(Chapman–Kolmogorov equation)，可以由狀態轉移函數與時間$k-1$的機率密度函數$p(x_{k-1}\mid y_{1:k-1})$，計算出時間$k$的先驗機率$p(x_k\mid y_{1:k-1})$

$$\begin{align} p(x_k\mid y_{1:k-1}) & = \int p(x_k,x_{k-1}\mid y_{1:k-1})d{x_{k-1}} \\ & = \int p(x_k\mid x_{k-1}, y_{1:k-1})p(x_{k-1}\mid y_{1:k-1})d{x_{k-1}} \\ & = \int p(x_k\mid x_{k-1})p(x_{k-1}\mid y_{1:k-1})d{x_{k-1}} \\ \end{align}$$

其中，由於狀態轉移模型被假設為一階馬可夫過程，時間$k-1$的狀態只由時間$k$決定，因此$p(x_k\mid x_{k-1}, y_{1:k-1})=p(x_k\mid x_{k-1})$。機率模型$p(x_k\mid x_{k-1})$由狀態轉移函數$x_k=f_k(x_{k-1},v_{k-1})$和$v_{k-1}$的統計值得到。

### 更新

在時間$k$，我們得到觀測值$y_k$，因此可以利用貝氏定理，由先驗機率$p(x_k\mid y_{1:k-1})$得到後驗機率$p(x_k\mid y_{1:k})$，也就是考慮觀測值後得到的機率。

$$p(x_k\mid y_{1:k})={p(y_k\mid x_k)p(x_k\mid y_{1:k-1})\over p(y_k\mid y_{1:k-1})}$$

其中的歸一化常數為

$$p(y_k\mid y_{1:k-1})=\int p(y_k\mid x_k)p(x_k\mid y_{1:k-1})dx_k$$

其中的似然函數$p(y_k\mid x_k)$由觀測函數$y_k=h_k(x_k,n_k)$和$n_k$的統計值得到。

上述「預測」與「更新」的遞迴關係，是貝葉斯最佳解的基本概念，然而公式中運用到的積分，對於一般非線性、非高斯的系統，難以得到解析解，因此需要利用蒙地卡羅方法來近似貝葉斯最佳解。

## 序列重要性重採樣（Sequential Importance Resampling, SIR）

### 序列重要性採樣（Sequential Importance Sampling, SIS）

SIR是由SIS加上重採樣（resampling）改良而來，在SIS中，我們將後驗機率$p(x_k\mid y_{1:k})$用$N$個隨機採樣的樣本(稱為粒子)與各自的權重表示為

$${\{x_k^{(i)},w_k^{(i)}\}}_{i=1}^N$$

其中的權重是根據重要性採樣的規則產生，且必須符合

$$\sum_{i=1}^N w_k^{(i)}=1$$

SIS是將重要性採樣遞迴執行的一種方法，根據重要性採樣的理論，一個函數$f$的期望值可以用加權平均來近似

$$\int f(x_k)p(x_k\mid y_{1:k})dx_k\approx \sum_{i=1}^N w_k^{(i)}f(x_k^{(i)})=1$$

在每一次遞迴過程中，我們希望由前一次採樣的權重，計算出下一次採樣的權重。假設採樣的樣本分布可以表示為

$x^{(i)}\sim q(x),$ $i=1,...,N$

其中$q(x)$稱為重要性密度（importance density）。若樣本$x_k^{(i)}$是由重要性密度$q(x_k\mid y_{1:k})$抽取出來，則權重可表示為

$$w_k^{(i)}\propto{p(x_k^{(i)}\mid y_{1:k})\over q(x_k^{(i)}\mid y_{1:k})}$$

我們將重要性密度分解為

$$q(x_k\mid y_{1:k})=q(x_k\mid x_{k-1},y_{1:k})q(x_{k-1}\mid y_{1:k-1})$$

再將後驗機率表示為

$$\begin{align} p(x_k\mid y_{1:k}) & = {p(y_k\mid x_k,y_{1:k-1})p(x_k\mid y_{1:k-1})\over p(y_k\mid y_{1:k-1})} \\ & = {p(y_k\mid x_k,y_{1:k-1})p(x_k\mid x_{k-1},y_{1:k-1})p(x_{k-1}\mid y_{1:k-1})\over p(y_k\mid y_{1:k-1})} \\ & = {p(y_k\mid x_k)p(x_k\mid x_{k-1})p(x_{k-1}\mid y_{1:k-1})\over p(y_k\mid y_{1:k-1})} \\ & \propto p(y_k\mid x_k)p(x_k\mid x_{k-1})p(x_{k-1}\mid y_{1:k-1}) \\ \end{align}$$

則權重的遞迴式可以表示為

$$\begin{align} w_k^{(i)} & \propto{p(y_k\mid x_k^{(i)})p(x_k^{(i)}\mid x_{k-1}^{(i)})p(x_{k-1}^{(i)}\mid y_{1:k-1})\over q(x_k^{(i)}\mid x_{k-1}^{(i)},y_{1:k})q(x_{k-1}^{(i)}\mid y_{1:k-1})} \\ & = w_{k-1}^{(i)}{p(y_k\mid x_k^{(i)})p(x_k^{(i)}\mid x_{k-1}^{(i)})\over q(x_k^{(i)}\mid x_{k-1}^{(i)},y_{1:k})} \\ \end{align}$$

### 重採樣（resampling）

SIS可能會造成退化問題（degeneracy problem），也就是在經過幾次遞迴後，很多粒子的權重都變小到可以忽略不計，只剩少數粒子的權重較大，如此會浪費大量的計算量在幾乎沒有作用的粒子上，而使估計性能下降。由於在遞迴過程中，權重的變異數只會愈來愈大，因此退化問題無法被避免。

為了評估退化問題，我們定義有效粒子數為

$$\widehat{N_{eff}}={1\over\sum_{i=1}^N (w_k^{(i)})^2}$$

在進行SIS時，若有效粒子數小於某一閾值，則對粒子做重採樣，即可減緩退化問題。重採樣的概念是去除權重過小的粒子，專注於權重較大的粒子。進行重採樣時，要由現有的粒子分布取樣，產生一組新的粒子${\{x_k^{(i)*}\}}_{i=1}^N$，由於產生出的新樣本為獨立且相同分布（i.i.d.），因此將權重重新設定為$w_k^{(i)}={1\over N}$。

## 參考文獻

- M. S. Arulampalam, S. Maskell, N. Gordon and T. Clapp, "A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking." In IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174-188, Feb 2002.
- A. Doucet, N. de Freitas, N. Gordon, "An Introduction to Sequential Monte Carlo Methods." In A. Doucet, N. de Freitas, N. Gordon (eds.) Sequential Monte Carlo Methods in Practice. Statistics for Engineering and Information Science. Springer, New York, NY, 2001
- A. Doucet, "On sequential Monte Carlo methods for Bayesian filtering." Dept. Eng., Univ. Cambridge, UK, Tech. Rep., 1998.