粒子滤波器

粒子滤波器（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.