粒子滤波理论

本文由最基础的贝叶斯估计开始介绍，再引出蒙特卡罗采样，重要性采样，SIS粒子滤波，重采样，基本粒子滤波Generic Particle Filter，SIR粒子滤波，这些概念的引进，都是为了解决上一个概念中出现的问题而环环相扣的。
粒子滤波通过非参数化的蒙特卡洛(Monte Carlo)模拟方法来实现递推贝叶斯滤波，适用于任何能用状态空间模型描述的非线性系统，精度可以逼近最优估计。粒子滤波器具有简单、易于实现等特点，它为分析非线性动态系统提供了一种有效的解决方法，从而引起目标跟踪、信号处理以及自动控制等领域的广泛关注。

粒子滤波(PF: Particle Filter)的思想基于蒙特卡洛方法(Monte Carlo methods)，它是利用粒子集来表示概率，可以用在任何形式的状态空间模型上。其核心思想是通过从后验概率中抽取的随机状态粒子来表达其分布，是一种顺序重要性采样法(Sequential Importance Sampling)。

一、贝叶斯滤波

动态系统的目标跟踪问题可以通过图所示的状态空间模型来描述。本节在贝叶斯滤波框架下讨论目标跟踪问题。

假设有一个系统，我们知道它的状态方程，和测量方程如下：
xk=fk(xk−1,vk−1)如: xk=xk−12+25xk−11+xk−12+8cos⁡(1.2(k−1))+vk(1)x_{k}=f_{k}\left(x_{k-1}, v_{k-1}\right) \text { 如: } x_{k}=\frac{x_{k-1}}{2}+\frac{25 x_{k-1}}{1+x_{k-1}^{2}}+8 \cos (1.2(k-1))+v_{k} (1) xk=fk(xk−1,vk−1) 如: xk=2xk−1+1+xk−1225xk−1+8cos(1.2(k−1))+vk(1)
yk=hk(xk,nk)如: yk=xk220+nk(2)y_{k}=h_{k}\left(x_{k}, n_{k}\right) \quad \text { 如: } y_{k}=\frac{x_{k}^{2}}{20}+n_{k}(2) yk=hk(xk,nk) 如: yk=20xk2+nk(2)

其中x为系统状态，y为测量到的数据，f,h是状态转移函数和测量函数，v,n为过程噪声和测量噪声，噪声都是独立同分布的。上面对应的那个例子将会出现在程序中。

从贝叶斯理论的观点来看，状态估计问题（目标跟踪、信号滤波）就是根据之前一系列的已有数据y1:ky_{1: k}y1:k（后验知识）递推的计算出当前状态xkx_{k}xk的可信度。这个可信度就是概率公式p(xk∣y1:k)p\left(x_{k} \mid y_{1:k}\right)p(xk∣y1:k)，它需要通过预测和更新两个步奏来递推的计算。

1.1两个过程

预测过程是利用系统模型(状态方程1)预测状态的先验概率密度，也就是通过已有的先验知识对未来的状态进行猜测，即p(xk∣xk−1)p\left(x_{k} \mid x_{k-1}\right)p(xk∣xk−1)。

更新过程则利用最新的测量值对先验概率密度进行修正，得到后验概率密度，也就是对之前的猜测进行修正。

在处理这些问题时，一般都先假设系统的状态转移服从一阶马尔科夫模型，即当前时刻的状态x(k)只与上一个时刻的状态x(k-1)有关。这是很自然的一种假设，就像小时候玩飞行棋，下一时刻的飞机跳到的位置只由当前时刻的位置和骰子决定。同时，假设k时刻测量到的数据yky_{k}yk只与当前的状态xkx_{k}xk有关，如上面的状态方程2。

为了进行递推，不妨假设已知k-1时刻的概率密度函数p(xk−1∣y1:k−1)p\left(x_{k-1} \mid y_{1:k-1}\right)p(xk−1∣y1:k−1);

1.1.1预测过程

预测：由上一时刻的概率密度p(xk−1∣y1:k−1)p\left(x_{k-1} \mid y_{1:k-1}\right)p(xk−1∣y1:k−1)得到p(xk∣y1:k−1)p\left(x_{k} \mid y_{1:k-1}\right)p(xk∣y1:k−1)，这个公式的含义是既然有了前面（1：k-1）时刻的测量数据，那就可以预测一下状态xkx_{k}xk出现的概率。

计算推导如下：
p(xk∣y1:k−1)=∫p(xk,xk−1∣y1:k−1)dxk−1=∫p(xk∣xk−1,y1:k−1)p(xk−1∣y1:k−1)dxk−1=∫p(xk∣xk−1)p(xk−1∣y1:k−1)dxk−1\begin{aligned} p\left(x_{k} \mid y_{1: k-1}\right) &=\int p\left(x_{k}, x_{k-1} \mid y_{1: k-1}\right) d x_{k-1} \\ &=\int p\left(x_{k} \mid x_{k-1}, y_{1: k-1}\right) p\left(x_{k-1} \mid y_{1: k-1}\right) d x_{k-1} \\ &=\int p\left(x_{k} \mid x_{k-1}\right) p\left(x_{k-1} \mid y_{1: k-1}\right) d x_{k-1} \end{aligned} p(xk∣y1:k−1)=∫p(xk,xk−1∣y1:k−1)dxk−1=∫p(xk∣xk−1,y1:k−1)p(xk−1∣y1:k−1)dxk−1=∫p(xk∣xk−1)p(xk−1∣y1:k−1)dxk−1

等式的第一行到第二行纯粹是贝叶斯公式的应用。第二行得到第三行是由于一阶马尔科夫过程的假设，状态xkx_{k}xk只由xk−1x_{k-1}xk−1决定。

两个问题：

第一：既然p(xk−1∣xk−1,y1:k−1)=p(xk∣xk−1)p\left(x_{k-1} \mid x_{k-1}, y_{1: k-1}\right)=p\left(x_{k} \mid x_{k-1}\right)p(xk−1∣xk−1,y1:k−1)=p(xk∣xk−1)，xkx_{k}xk都只由xk−1x_{k-1}xk−1决定了，即p(xk∣xk−1)p\left(x_{k} \mid x_{k-1}\right)p(xk∣xk−1)，在这里弄一个p(xk∣y1:k−1)p\left(x_{k} \mid y_{1:k-1}\right)p(xk∣y1:k−1)是什么意思？

这两个概率公式含义是不一样的，前一个是纯粹根据模型进行预测，x(k)x\left(k\right)x(k)实实在在的由x(k−1)x\left(k-1\right)x(k−1)决定，后一个是既然测到的数据和状态是有关系的，现在已经有了很多测量数据y了，那么我可以根据已有的经验对你进行预测，只是猜测xkx_{k}xk，而不能决定xkx_{k}xk;

第二：上面公式的最后一行p(xk−1∣y1:k−1)p\left(x_{k-1} \mid y_{1: k-1}\right)p(xk−1∣y1:k−1)是假设已知的，但是p(xk∣xk−1)p\left(x_{k} \mid x_{k-1}\right)p(xk∣xk−1)怎么得到呢？

其实p(xk∣xk−1)p\left(x_{k} \mid x_{k-1}\right)p(xk∣xk−1)它由系统状态方程决定，它的概率分布形状和系统的过程噪声vk−1v_{k-1}vk−1形状一模一样。如何理解呢?观察状态方程(1)式我们知道，x(k)=Constant⁡(x(k−1))+v(k−1)x(k)=\operatorname{Constant}(x(k-1))+v(k-1)x(k)=Constant(x(k−1))+v(k−1) 也就是(xk)\left(x_{k}\right)(xk)由一个通过(xk−1)\left(x_{k-1}\right)(xk−1)计算出的常数叠加一个噪声得到。看下面的图：

如果没有噪声，xkx_{k}xk完全由xk−1x_{k-1}xk−1计算得到，也就没由概率分布这个概念了，由于出现了噪声，所以xkx_{k}xk不好确定，他的分布就如同图中的阴影部分，实际上形状和噪声是一样的，只是进行了一些平移。理解了这一点，对粒子滤波程序中，状态xkx_{k}xk的采样的计算很有帮助，要采样xkx_{k}xk，直接采样一个过程噪声，再叠加上f(xk−1)f(x_{k-1})f(xk−1) 这个常数就行了。

1.1.2更新过程

更新：由p(xk∣y1:k−1)p\left(x_{k} \mid y_{1: k-1}\right)p(xk∣y1:k−1)(上式所得)得到p(xk∣y1:k)p\left(x_{k} \mid y_{1: k}\right)p(xk∣y1:k)后验概率。这个后验概率才是真正有用的东西，上一步还只是预测，这里又多了k时刻的测量，对上面的预测再进行修正，就是滤波了。这里的后验概率也将是代入到下次的预测，形成递推。
推导：
p(xk∣y1:k)=p(yk∣xk,y1:k−1)p(xk∣y1:k−1)p(yk∣y1:k−1)=p(yk∣xk)p(xk∣y1:k−1)p(yk∣y1:k−1)\begin{aligned} p\left(x_{k} \mid y_{1: k}\right) &=\frac{p\left(y_{k} \mid x_{k}, y_{1: k-1}\right) p\left(x_{k} \mid y_{1: k-1}\right)}{p\left(y_{k} \mid y_{1: k-1}\right)} \\ &=\frac{p\left(y_{k} \mid x_{k}\right) p\left(x_{k} \mid y_{1: k-1}\right)}{p\left(y_{k} \mid y_{1: k-1}\right)} \end{aligned} p(xk∣y1:k)=p(yk∣y1:k−1)p(yk∣xk,y1:k−1)p(xk∣y1:k−1)=p(yk∣y1:k−1)p(yk∣xk)p(xk∣y1:k−1)
其中归一化常数：
p(yk∣y1:k−1)=∫p(yk∣xk)p(xk∣y1:k−1)dxkp\left(y_{k} \mid y_{1: k-1}\right)=\int p\left(y_{k} \mid x_{k}\right) p\left(x_{k} \mid y_{1: k-1}\right) d x_{k} p(yk∣y1:k−1)=∫p(yk∣xk)p(xk∣y1:k−1)dxk
等式第一行到第二行是因为测量方程知道, y(k)只与x(k)有关，p(yk∣xk)p\left(y_{k} \mid x_{k}\right)p(yk∣xk)也称之为似然函数，由测量方程决定。也和上面的推理一样，yk=h(xk)+nky_{k}=h\left(x_{k}\right)+n_{k}yk=h(xk)+nk, x(k)部分是常数，p(yk∣xk)p\left(y_{k} \mid x_{k}\right)p(yk∣xk)也是只和测量噪声n(k)的概率分布有关，注意这个也将为SIR粒子滤波里权重的采样提供编程依据。
但是，请注意上面的推导过程中需要用到积分，这对于一般的非线性，非高斯系统，很难得到后验概率的解析解。为了解决这个问题，就得引进蒙特卡洛采样。

二、蒙特卡洛采样

假设我们能从一个目标概率分布p(x)p\left(x\right)p(x)中采样到一系列的样本（粒子）x1,...,xnx_{1},...,x_{n}x1,...,xn；（至于怎么生成服从p(x)分布的样本，这个问题先放一放），那么就能利用这些样本去估计这个分布的某些函数的期望值。譬如：
E(f(x))=∫abf(x)p(x)dxE(f(x))=\int_{a}^{b} f(x) p(x) d x E(f(x))=∫abf(x)p(x)dx
Var⁡(f(x))=E(f(x)−E(f(x)))2=∫ab(f(x)−E(f(x)))2p(x)dx\operatorname{Var}(f(x))=E(f(x)-E(f(x)))^{2}=\int_{a}^{b}(f(x)-E(f(x)))^{2} p(x) d x Var(f(x))=E(f(x)−E(f(x)))2=∫ab(f(x)−E(f(x)))2p(x)dx
上面的式子其实都是计算期望的问题，只是被积分的函数不同。

蒙特卡洛采样的思想就是用平均值来代替积分，求期望：
E(f(x))≈f(x1)+…+f(xn)nE(f(x)) \approx \frac{f\left(x_{1}\right)+\ldots+f\left(x_{n}\right)}{n} E(f(x))≈nf(x1)+…+f(xn)
这可以从大数定理的角度去理解它。我们用这种思想去指定不同的f(x)f(x)f(x)以便达到估计不同东西的目的。比如：要估计一批同龄人的体重，不分男女，在大样本中男的有100个，女的有20个，为了少做事，我们按比例抽取10个男的，2个女的，测算这12个人的体重求平均就完事了。注意这里的按比例抽取，就可以看成从概率分布p(x)中进行抽样。

由上面我们知道，蒙特卡洛采样可以用来估计概率，而在上一节中，贝叶斯后验概率的计算里要用到积分，为了解决这个积分难的问题，可以用蒙特卡洛采样来代替计算后验概率。

假设可以从后验概率中采样到N个样本，那么后验概率的计算可表示为：
p^(xn∣y1:k)=1N∑i=1Nδ(xn−xn(i))≈p(xn∣y1:k)\hat{p}\left(x_{n} \mid y_{1: k}\right)=\frac{1}{N} \sum_{i=1}^{N} \delta\left(x_{n}-x_{n}^{(i)}\right) \approx p\left(x_{n} \mid y_{1: k}\right) p^(xn∣y1:k)=N1i=1∑Nδ(xn−xn(i))≈p(xn∣y1:k)

其中，在这个蒙特卡洛方法中，我们定义f(x)=δ(xn−xn(i))f(x)=\delta\left(x_{n}-x_{n}^{(i)}\right)f(x)=δ(xn−xn(i)),是狄拉克函数(dirac delta function)，跟上面的指示函数意思差不多。
看到这里，既然用蒙特卡洛方法能够用来直接估计后验概率，现在估计出了后验概率，那到底怎么用来做图像跟踪或者滤波呢？要做图像跟踪或者滤波，其实就是想知道当前状态的期望值：
E[f(xn)]≈∫f(xn)p^(xn∣y1:k)dxn=1N∑i=1N∫f(xn)δ(xn−xn(i))dxn=1N∑i=1Nf(xn(i))\begin{aligned} E\left[f\left(x_{n}\right)\right] \approx & \int f\left(x_{n}\right) \hat{p}\left(x_{n} \mid y_{1: k}\right) d x_{n} \\ &=\frac{1}{N} \sum_{i=1}^{N} \int f\left(x_{n}\right) \delta\left(x_{n}-x_{n}^{(i)}\right) d x_{n} \\ &=\frac{1}{N} \sum_{i=1}^{N} f\left(x_{n}^{(i)}\right) \end{aligned} E[f(xn)]≈∫f(xn)p^(xn∣y1:k)dxn=N1i=1∑N∫f(xn)δ(xn−xn(i))dxn=N1i=1∑Nf(xn(i))
也就是用这些采样的粒子的状态值直接平均就得到了期望值，也就是滤波后的值，这里的f(x)f(x)f(x) 就是每个粒子的状态函数。这就是粒子滤波了，只要从后验概率中采样很多粒子，用它们的状态求平均就得到了滤波结果。
思路看似简单，但是要命的是，后验概率不知道啊，怎么从后验概率分布中采样！所以这样直接去应用是行不通的，这时候得引入重要性采样这个方法来解决这个问题。

三、重要性采样

无法从目标分布中采样，就从一个已知的可以采样的分布里去采样如 q(x|y)，这样上面的求期望问题就变成了：
E[f(xk)]=∫f(xk)p(xk∣y1:k)q(xk∣y1:k)q(xk∣y1:k)dxk=∫f(xk)p(y1:k∣xk)p(xk)p(y1:k)q(xk∣y1:k)q(xk∣y1:k)dxk=∫f(xk)Wk(xk)p(y1:k)q(xk∣y1:k)dxk\begin{aligned} E\left[f\left(x_{k}\right)\right] &=\int f\left(x_{k}\right) \frac{p\left(x_{k} \mid y_{1: k}\right)}{q\left(x_{k} \mid y_{1: k}\right)} q\left(x_{k} \mid y_{1: k}\right) d x_{k} \\ &=\int f\left(x_{k}\right) \frac{p\left(y_{1: k} \mid x_{k}\right) p\left(x_{k}\right)}{p\left(y_{1: k}\right) q\left(x_{k} \mid y_{1: k}\right)} q\left(x_{k} \mid y_{1: k}\right) d x_{k} \\ &=\int f\left(x_{k}\right) \frac{W_{k}\left(x_{k}\right)}{p\left(y_{1: k}\right)} q\left(x_{k} \mid y_{1: k}\right) d x_{k} \end{aligned} E[f(xk)]=∫f(xk)q(xk∣y1:k)p(xk∣y1:k)q(xk∣y1:k)dxk=∫f(xk)p(y1:k)q(xk∣y1:k)p(y1:k∣xk)p(xk)q(xk∣y1:k)dxk=∫f(xk)p(y1:k)Wk(xk)q(xk∣y1:k)dxk
其中
Wk(xk)=p(y1:k∣xk)p(xk)q(xk∣y1:k)∝p(xk∣y1:k)q(xk∣y1:k)W_{k}\left(x_{k}\right)=\frac{p\left(y_{1: k} \mid x_{k}\right) p\left(x_{k}\right)}{q\left(x_{k} \mid y_{1: k}\right)} \propto \frac{p\left(x_{k} \mid y_{1: k}\right)}{q\left(x_{k} \mid y_{1: k}\right)} Wk(xk)=q(xk∣y1:k)p(y1:k∣xk)p(xk)∝q(xk∣y1:k)p(xk∣y1:k)
由于：
p(y1:k)=∫p(y1:k∣xk)p(xk)dxkp\left(y_{1: k}\right)=\int p\left(y_{1: k} \mid x_{k}\right) p\left(x_{k}\right) d x_{k} p(y1:k)=∫p(y1:k∣xk)p(xk)dxk
所以上式可以进一步写成（3）式：
E[f(xk)]=1p(y1:k)∫f(xk)Wk(xk)q(xk∣y1:k)dxk=∫f(xk)Wk(xk)q(xk∣y1:k)dxk∫p(y1:k∣xk)p(xk)dxk=∫f(xk)Wk(xk)q(xk∣y1:k)dxk∫Wk(xk)q(xk∣y1:k)dxk=Eq(xk∣y1:k)[Wk(xk)f(xk)]Eq(xk∣y1:k)[Wk(xk)]\begin{aligned} E\left[f\left(x_{k}\right)\right]=& \frac{1}{p\left(y_{1: k}\right)} \int f\left(x_{k}\right) W_{k}\left(x_{k}\right) q\left(x_{k} \mid y_{1: k}\right) d x_{k} \\ &=\frac{\int f\left(x_{k}\right) W_{k}\left(x_{k}\right) q\left(x_{k} \mid y_{1: k}\right) d x_{k}}{\int p\left(y_{1: k} \mid x_{k}\right) p\left(x_{k}\right) d x_{k}} \\ &=\frac{\int f\left(x_{k}\right) W_{k}\left(x_{k}\right) q\left(x_{k} \mid y_{1: k}\right) d x_{k}}{\int W_{k}\left(x_{k}\right) q\left(x_{k} \mid y_{1: k}\right) d x_{k}} \\ =& \frac{E_{q\left(x_{k} \mid y_{1: k}\right)}\left[W_{k}\left(x_{k}\right) f\left(x_{k}\right)\right]}{E_{q\left(x_{k} \mid y_{1: k}\right)}\left[W_{k}\left(x_{k}\right)\right]} \end{aligned} E[f(xk)]==p(y1:k)1∫f(xk)Wk(xk)q(xk∣y1:k)dxk=∫p(y1:k∣xk)p(xk)dxk∫f(xk)Wk(xk)q(xk∣y1:k)dxk=∫Wk(xk)q(xk∣y1:k)dxk∫f(xk)Wk(xk)q(xk∣y1:k)dxkEq(xk∣y1:k)[Wk(xk)]Eq(xk∣y1:k)[Wk(xk)f(xk)]
上面的期望计算都可以通过蒙特卡洛方法来解决它，也就是说，通过采样N个样本{xk(i)}∼q(xk∣y1:k)\left\{x_{k}^{(i)}\right\} \sim q\left(x_{k} \mid y_{1: k}\right){xk(i)}∼q(xk∣y1:k),用样本的平均来求它们的期望，所以上面的（3）式可以近似为（4）：
E[f(xk)]≈1N∑i=1NWk(xk(i))f(xk(i))1N∑i=1NWk(xk(i))=∑i=1NW~k(xk(i))f(xk(i))\begin{aligned} E\left[f\left(x_{k}\right)\right] \approx & \frac{\frac{1}{N} \sum_{i=1}^{N} W_{k}\left(x_{k}^{(i)}\right) f\left(x_{k}^{(i)}\right)}{\frac{1}{N} \sum_{i=1}^{N} W_{k}\left(x_{k}^{(i)}\right)} \\ &=\sum_{i=1}^{N} \tilde{W}_{k}\left(x_{k}^{(i)}\right) f\left(x_{k}^{(i)}\right) \end{aligned} E[f(xk)]≈N1∑i=1NWk(xk(i))N1∑i=1NWk(xk(i))f(xk(i))=i=1∑NW~k(xk(i))f(xk(i))
其中:
W~k(xk(i))=Wk(xk(i))∑i=1NWk(xk(i))\tilde{W}_{k}\left(x_{k}^{(i)}\right)=\frac{W_{k}\left(x_{k}^{(i)}\right)}{\sum_{i=1}^{N} W_{k}\left(x_{k}^{(i)}\right)} W~k(xk(i))=∑i=1NWk(xk(i))Wk(xk(i))
这就是归一化以后的权重，而之前在(2)式中的那个权重是没有归一化的。

注意上面的(4)式，它不再是（1）式中所有的粒子状态直接相加求平均了，而是一种加权和的形式。不同的粒子都有它们相应的权重，如果粒子权重大，说明信任该粒子比较多。

到这里已经解决了不能从后验概率直接采样的问题，但是上面这种每个粒子的权重都直接计算的方法，效率低，因为每增加一个采样，p(x(k)∣y(1:k))p(x(k) \mid y(1: k))p(x(k)∣y(1:k))都得重新计算，并且还不好计算这个式子。所以求权重时能否避开计算p(x(k)∣y(1:k))p(x(k) \mid y(1: k))p(x(k)∣y(1:k))？而最佳的形式是能够以递推的方式去计算权重，这就是所谓的序贯重要性采样（SIS），粒子滤波的原型。

下面开始权重w递推形式的推导：
假设重要性概率密度函数q(x0:k∣y1:k)q\left(x_{0:k} \mid y_{1: k}\right)q(x0:k∣y1:k)（这里x的下标是0:k，也就是说粒子滤波是估计过去所有时刻的状态的后验）。假设它可以分解为：
q(x0:k∣y1:k)=q(x0:k−1∣y1:k−1)q(xk∣x0:k−1,y1:k)q\left(x_{0:k} \mid y_{1:k}\right)=q\left(x_{0:k-1} \mid y_{1: k-1}\right) q\left(x_{k} \mid x_{0: k-1}, y_{1: k}\right) q(x0:k∣y1:k)=q(x0:k−1∣y1:k−1)q(xk∣x0:k−1,y1:k)
后验概率密度函数的递归形式可以表示为：
p(x0:k∣Yk)=p(yk∣x0:k,Yk−1)p(x0:k∣Yk−1)p(yk∣Yk−1)=p(yk∣x0.k,Yk−1)p(xk∣x0.k−1,Yk−1)p(x0.k−1∣Yk−1)p(yk∣Yk−1)=p(yk∣xk)p(xk∣xk−1)p(x0:k−1∣Yk−1)p(yk∣Yk−1)∝p(yk∣xk)p(xk∣xk−1)p(x0:k−1∣Yk−1)\begin{array}{c} p\left(x_{0: k} \mid Y_{k}\right)=\frac{p\left(y_{k} \mid x_{0: k}, Y_{k-1}\right) p\left(x_{0: k} \mid Y_{k-1}\right)}{p\left(y_{k} \mid Y_{k-1}\right)} \\ =\frac{p\left(y_{k} \mid x_{0 . k}, Y_{k-1}\right) p\left(x_{k} \mid x_{0 . k-1}, Y_{k-1}\right) p\left(x_{0 . k-1} \mid Y_{k-1}\right)}{p\left(y_{k} \mid Y_{k-1}\right)} \\ =\frac{p\left(y_{k} \mid x_{k}\right) p\left(x_{k} \mid x_{k-1}\right) p\left(x_{0: k-1} \mid Y_{k-1}\right)}{p\left(y_{k} \mid Y_{k-1}\right)} \\ \propto p\left(y_{k} \mid x_{k}\right) p\left(x_{k} \mid x_{k-1}\right) p\left(x_{0: k-1} \mid Y_{k-1}\right) \end{array} p(x0:k∣Yk)=p(yk∣Yk−1)p(yk∣x0:k,Yk−1)p(x0:k∣Yk−1)=p(yk∣Yk−1)p(yk∣x0.k,Yk−1)p(xk∣x0.k−1,Yk−1)p(x0.k−1∣Yk−1)=p(yk∣Yk−1)p(yk∣xk)p(xk∣xk−1)p(x0:k−1∣Yk−1)∝p(yk∣xk)p(xk∣xk−1)p(x0:k−1∣Yk−1)
其中，为了表示方便，将y(1:k)y(1:k)y(1:k) 用 Y(k)Y(k)Y(k)来表示，注意 YYY 与 yyy 的区别。同时，上面这个式子和上一节贝叶斯滤波中后验概率的推导是一样的，只是之前的x(k)变成了这里的x(0:k)，就是这个不同，导致贝叶斯估计里需要积分，而这里后验概率的分解形式却不用积分。
粒子权值的递归形式可以表示为(5):
wk(i)∞p(x0:k(i)∣Yk)q(x0.k(i)∣Yk)=p(yk∣xk(i))p(xk(i)∣xk−1(i))p(x0:k−1(i)∣Yk−1)q(xk(i)∣x0.k−1(i),Yk)q(x0:k−1(i)∣Yk−1)=wk−1(i)p(yk∣xk(i))p(xk(i)∣xk−1(i))q(xk(i)∣x0.k−1(i),Yk)\begin{array}{l} w_{k}^{(i)} \infty \frac{p\left(x_{0: k}^{(i)} \mid Y_{k}\right)}{q\left(x_{0 . k}^{(i)} \mid Y_{k}\right)} \\ =\frac{p\left(y_{k} \mid x_{k}^{(i)}\right) p\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}\right) p\left(x_{0: k-1}^{(i)} \mid Y_{k-1}\right)}{q\left(x_{k}^{(i)} \mid x_{0 . k-1}^{(i)}, Y_{k}\right) q\left(x_{0: k-1}^{(i)} \mid Y_{k-1}\right)} \\ =w_{k-1}^{(i)} \frac{p\left(y_{k} \mid x_{k}^{(i)}\right) p\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}\right)}{q\left(x_{k}^{(i)} \mid x_{0 . k-1}^{(i)}, Y_{k}\right)} \end{array} wk(i)∞q(x0.k(i)∣Yk)p(x0:k(i)∣Yk)=q(xk(i)∣x0.k−1(i),Yk)q(x0:k−1(i)∣Yk−1)p(yk∣xk(i))p(xk(i)∣xk−1(i))p(x0:k−1(i)∣Yk−1)=wk−1(i)q(xk(i)∣x0.k−1(i),Yk)p(yk∣xk(i))p(xk(i)∣xk−1(i))
注意，这种权重递推形式的推导是在前面（2）式的形式下进行推导的，也就是没有归一化。而在进行状态估计的公式为∑i=1NW~k(xk(i))f(xk(i))\sum_{i=1}^{N} \tilde{W}_{k}\left(x_{k}^{(i)}\right) f\left(x_{k}^{(i)}\right)∑i=1NW~k(xk(i))f(xk(i))这个公式中的的权重是归一化以后的，所以在实际应用中，递推计算出w(k)后,要进行归一化，才能够代入(4)式中去计算期望。同时，上面(5)式中的分子是不是很熟悉，在上一节贝叶斯滤波中我们都已经做了介绍，p( y|x ),p( x(k)|x(k-1) )的形状实际上和状态方程中噪声的概率分布形状是一样的，只是均值不同了。因此这个递推的(5)式和前面的非递推形式相比，公式里的概率都是已知的，权重的计算可以说没有编程方面的难度了。权重也有了以后，只要进行稍微的总结就可以得到SIS Filter。

四、Sequential Importance Sampling (SIS) Filter

在实际应用中我们可以假设重要性分布q()满足：
q(xk∣x0:k−1,y1:k)=q(xk∣xk−1,yk)q\left(x_{k} \mid x_{0: k-1}, y_{1: k}\right)=q\left(x_{k} \mid x_{k-1}, y_{k}\right)q(xk∣x0:k−1,y1:k)=q(xk∣xk−1,yk)
这个假设说明重要性分布只和前一时刻的状态x(k-1)以及测量y(k)有关了，那么(5)式就可以转化为：
wk(i)∝wk−1(i)p(yk∣xk(i))p(xk(i)∣xk−1(i))q(xk(i)∣xk−1(i),yk)w_{k}^{(i)} \propto w_{k-1}^{(i)} \frac{p\left(y_{k} \mid x_{k}^{(i)}\right) p\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}\right)}{q\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}, y_{k}\right)} wk(i)∝wk−1(i)q(xk(i)∣xk−1(i),yk)p(yk∣xk(i))p(xk(i)∣xk−1(i))
在做了这么多假设和为了解决一个个问题以后，终于有了一个像样的粒子滤波算法了，他就是序贯重要性采样滤波。

这个算法就是粒子滤波的前身了。只是在实际应用中，又发现了很多问题，如粒子权重退化的问题，因此就有了重采样( resample )，就有了基本的粒子滤波算法,还有就是重要性概率密度q()的选择问题。

在这一章中，我们是用的重要性采样这种方法去解决的后验概率无法采样的问题。实际上，关于如何从后验概率采样，也就是如何生成特定概率密度的样本，有很多经典的方法（如拒绝采样，Markov Chain Monte Carlo，Metropolis-Hastings算法，Gibbs采样）。

五、重采样

在应用SIS 滤波的过程中，存在一个退化的问题。就是经过几次迭代以后，很多粒子的权重都变得很小，可以忽略了，只有少数粒子的权重比较大。并且粒子权值的方差随着时间增大，状态空间中的有效粒子数较少。随着无效采样粒子数目的增加，使得大量的计算浪费在对估计后验滤波概率分布几乎不起作用的粒子上，使得估计性能下降，如图所示。

通常采用有效粒子数来衡量粒子权值的退化程度，即
Neff=N/(1+var⁡(wk∗(i)))N_{e f f}=N /\left(1+\operatorname{var}\left(w_{k}^{*\left(i\right)}\right)\right) Neff=N/(1+var(wk∗(i)))
wk∗(i)=p(xk(i)∣y1:k)q(xk(i)∣xk−1(i),y1:k)w_{k}^{*_{(i)}}=\frac{p\left(x_{k}^{(i)} \mid y_{1: k}\right)}{q\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}, y_{1: k}\right)} wk∗(i)=q(xk(i)∣xk−1(i),y1:k)p(xk(i)∣y1:k)
这个式子的含义是，有效粒子数越小，即权重的方差越大，也就是说权重大的和权重小的之间差距大，表明权值退化越严重。在实际计算中，有效粒子数可以近似为：
N^eff≈1∑i=1N(wk(i))2\hat{N}_{e f f} \approx \frac{1}{\sum_{i=1}^{N}\left(w_{k}^{(i)}\right)^{2}} N^eff≈∑i=1N(wk(i))21
在进行序贯重要性采样时，若上式小于事先设定的某一阈值，则应当采取一些措施加以控制。克服序贯重要性采样算法权值退化现象最直接的方法是增加粒子数，而这会造成计算量的相应增加，影响计算的实时性。因此，一般采用以下两种途径：(1)选择合适的重要性概率密度函数；(2)在序贯重要性采样之后，采用重采样方法。
对于第一种方法：选取重要性概率密度函数的一个标准就是使得粒子权值的方差最小。关于这部分内容，还是推荐百度文库的那篇文章《粒子滤波理论》，他在这里也引申出来了几种不同的粒子滤波方法。

这里着重讲第二种方法，重采样。

重采样的思路是：既然那些权重小的不起作用了，那就不要了。要保持粒子数目不变，得用一些新的粒子来取代它们。找新粒子最简单的方法就是将权重大的粒子多复制几个出来，至于复制几个？那就在权重大的粒子里面让它们根据自己权重所占的比例去分配，也就是老大分身分得最多，老二分得次多，以此类推。下面以数学的形式来进行说明。
前面已经说明了求某种期望问题变成了这种加权和的形式：
p(xk∣y1:k)=∑i=1Nwk(i)δ(xk−xk(i))p\left(x_{k} \mid y_{1: k}\right)=\sum_{i=1}^{N} w_{k}^{(i)} \delta\left(x_{k}-x_{k}^{(i)}\right) p(xk∣y1:k)=i=1∑Nwk(i)δ(xk−xk(i))
通过重采样以后，希望表示成：
p~(xk∣y1:k)=∑j=1N1Nδ(xk−xk(j))=∑i=1NniNδ(xk−xk(i))\tilde{p}\left(x_{k} \mid y_{1: k}\right)=\sum_{j=1}^{N} \frac{1}{N} \delta\left(x_{k}-x_{k}^{(j)}\right)=\sum_{i=1}^{N} \frac{n_{i}}{N} \delta\left(x_{k}-x_{k}^{(i)}\right) p~(xk∣y1:k)=j=1∑NN1δ(xk−xk(j))=i=1∑NNniδ(xk−xk(i))
注意对比(1)和(2)。xk(i)x_{k}^{(i)}xk(i)是第k时刻的粒子。xk(j)x_{k}^{(j)}xk(j)是k时刻重采样以后的粒子。xk(i)x_{k}^{(i)}xk(i)其中n(i)是指粒子在产生新的粒子集xk(j)x_{k}^{(j)}xk(j)时被复制的次数。（2）式中第一个等号说明重采样以后，所有的粒子权重一样，都是1/N，只是有的粒子多出现了n(i)次。
思路有了，就看具体的操作方法了。在《On resampling algorithms for particle fi lters》这篇paper里讲了四种重采样的方法。这四种方法大同小异。如果你接触过遗传算法的话，理解起来就很容易，就是遗传算法中那种轮盘赌的思想。

在这里用个简单的例子来说明：

假设有3个粒子，在第k时刻的时候，他们的权重分别是0.1, 0.1 ,0.8, 然后计算他们的概率累计和(matlab 中为cumsum() )得到： [0.1, 0.2, 1]。接着，我们用服从[0,1]之间的均匀分布随机采样3个值，假设为0.15 , 0.38 和 0.54。也就是说，第二个粒子复制一次，第三个粒子复制两次。
对于上面的过程，还可以对着下面的图加深理解：

将重采样的方法放入之前的SIS算法中，便形成了基本粒子滤波算法。

标准的粒子滤波算法流程为∶
（1）粒子集初始化，k =0∶
对于i=1，2…，N，由先验p(x0)p\left(x_{0}\right)p(x0)生成采样粒子{x0(i)}i=1N\left\{x_{0}^{(i)}\right\}_{i=1}^{N}{x0(i)}i=1N点
（2）对于k= 1，2…，N，循环执行以下步骤∶

① 重要性采样∶对于i=1，2，…，N，从重要性概率密度中生成采样粒子{x~k(i)}i=1N\left\{\tilde{x}_{k}^{(i)}\right\}_{i=1}^{N}{x~k(i)}i=1N，计算粒子权值w~k(i)\tilde{w}_{k}^{(i)}w~k(i)，并进行归一化;
② 重采样∶对粒子集{x~k(i),w~k(i)}\left\{\tilde{x}_{k}^{(i)}, \tilde{w}_{k}^{(i)}\right\}{x~k(i),w~k(i)}进行重采样，重采样后的粒子集为{xk(i),1/N}\left\{x_{k}^{(i)}, 1 / N\right\}{xk(i),1/N};
③ 输出∶ 计算k时刻的状态估计值∶x^k=∑i=1Nx~k(i)w~k(i)\hat{x}_{k}=\sum_{i=1}^{N} \tilde{x}_{k}^{(i)} \tilde{w}_{k}^{(i)}x^k=∑i=1Nx~k(i)w~k(i)。
重采样的思路很简单，但是当你仔细分析权重的计算公式时：
Wk(xk)=p(y1:k∣xk)p(xk)q(xk∣y1:k)∝p(xk∣y1:k)q(xk∣y1:k)W_{k}\left(x_{k}\right)=\frac{p\left(y_{1: k} \mid x_{k}\right) p\left(x_{k}\right)}{q\left(x_{k} \mid y_{1: k}\right)} \propto \frac{p\left(x_{k} \mid y_{1: k}\right)}{q\left(x_{k} \mid y_{1: k}\right)} Wk(xk)=q(xk∣y1:k)p(y1:k∣xk)p(xk)∝q(xk∣y1:k)p(xk∣y1:k)
会有疑问，权重大的就多复制几次，这一定是正确的吗？权重大，如果是分子大，即后验概率大，那说明确实应该在后验概率大的地方多放几个粒子。但权重大也有可能是分母小造成的，这时候的分子也可能小，也就是实际的后验概率也可能小，这时候的大权重可能就没那么优秀了。何况，这种简单的重采样会使得粒子的多样性丢失，到最后可能都变成了只剩一种粒子的分身。在遗传算法中好歹还引入了变异来解决多样性的问题。当然，粒子滤波里也有专门的方法：正则粒子滤波，有兴趣的可以查阅相关资料。

至此，整个粒子滤波的流程已经清晰明朗了，在实际应用中还有一些不确定的就是重要性概率密度的选择。在下一章中，首先引出 SIR 粒子滤波，接着用SIR滤波来进行实践应用。

六、Sampling Importance Resampling Filter (SIR)

SIR滤波器很容易由前面的基本粒子滤波推导出来，只要对粒子的重要性概率密度函数做出特定的选择即可。在SIR中，选取：
q(xk(i)∣xk−1(i),yk)=p(xk(i)∣xk−1(i))q\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}, y_{k}\right)=p\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}\right) q(xk(i)∣xk−1(i),yk)=p(xk(i)∣xk−1(i))
p( x(k)|x(k-1) )这是先验概率，在第一章贝叶斯滤波预测部分已经说过怎么用状态方程来得到它。将这个式子代入到第二章SIS推导出的权重公式中：
wk(i)∝wk−1(i)p(yk∣xk(i))p(xk(i)∣xk−1(i))q(xk(i)∣xk−1(i),yk)w_{k}^{(i)} \propto w_{k-1}^{(i)} \frac{p\left(y_{k} \mid x_{k}^{(i)}\right) p\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}\right)}{q\left(x_{k}^{(i)} \mid x_{k-1}^{(i)}, y_{k}\right)} wk(i)∝wk−1(i)q(xk(i)∣xk−1(i),yk)p(yk∣xk(i))p(xk(i)∣xk−1(i))
得到:
wk(i)∝wk−1(i)p(yk∣xk(i))w_{k}^{(i)} \propto w_{k-1}^{(i)} p\left(y_{k} \mid x_{k}^{(i)}\right) \quad wk(i)∝wk−1(i)p(yk∣xk(i))
由之前的重采样我们知道，实际上每次重采样以后，有wk−1(i)=1Nw_{k-1}^{(i)}=\frac{1}{N}wk−1(i)=N1。

所以(1)式可以进一步简化成：
wk(i)∝p(yk∣xk(i))w_{k}^{(i)} \propto p\left(y_{k} \mid x_{k}^{(i)}\right) wk(i)∝p(yk∣xk(i))
这里又出来一个概率采样p(yk∣xki)p\left(y_{k} \mid x_{k}^{i}\right)p(yk∣xki) ，实际怎么得到这个概率，程序里面又怎么去采样呢？
先搞清这个概率p(yk∣xki)p\left(y_{k} \mid x_{k}^{i}\right)p(yk∣xki)的含义，它表示在状态x出现的条件下，测量y出现的概率。在机器人定位里面就是，在机器人处于位姿x时，此时传感器数据y出现的概率。更简单的例子是，我要找到一个年龄是14岁的男孩（状态x），身高为170（测量y）的概率。要知道y出现的概率，需要知道此时y的分布。这里以第一篇文章的状态方程为例，由系统状态方程可知，测量是在真实值附近添加了一个高斯噪声。因此，y的分布就是以真实测量值为均值，以噪声方差为方差的一个高斯分布。因此，权重的采样过程就是：当粒子处于x状态时，能够得到该粒子的测量y。要知道这个测量y出现的概率，就只要把它放到以真实值为均值，噪声方差为方差的高斯分布里去计算就行了：
w=η(2πΣ)−1/2exp⁡{−12(ytrue −y)Σ−1(ytrue −y)}w=\eta(2 \pi \Sigma)^{-1 / 2} \exp \left\{-\frac{1}{2}\left(y_{\text {true }}-y\right) \Sigma^{-1}\left(y_{\text {true }}-y\right)\right\} w=η(2πΣ)−1/2exp{−21(ytrue −y)Σ−1(ytrue −y)}
到这里，就可以看成SIR只和系统状态方程有关了，不用自己另外去设计概率密度函数，所以在很多程序中都是用的这种方法。

下面以伪代码的形式给出SIR滤波器：

在上面算法中，每进行一次，都必须重采样一次，这是由于权重的计算方式决定的。

分析上面算法中的采样，发现它并没有加入测量y(k)。只是凭先验知识p(x(k)∣x(k−1))p( x(k)|x(k-1) )p(x(k)∣x(k−1))进行的采样，而不是用的修正了的后验概率。所以这种算法存在效率不高和对奇异点(outliers)敏感的问题。但不管怎样，SIR确实简单易用。

七、粒子滤波的应用

在这里主要以第一章的状态方程作为例子进行演示。
xk=xk−12+25xk−11+xk−12+8cos⁡(1.2(k−1))+vkx_{k}=\frac{x_{k-1}}{2}+\frac{25 x_{k-1}}{1+x_{k-1}^{2}}+8 \cos (1.2(k-1))+v_{k} xk=2xk−1+1+xk−1225xk−1+8cos(1.2(k−1))+vk
yk=xk220+nky_{k}=\frac{x_{k}^{2}}{20}+n_{k} yk=20xk2+nk
在这个存在过程噪声和量测噪声的系统中，估计状态x(k)。

%% SIR粒子滤波的应用，算法流程参见博客http://blog.csdn.net/heyijia0327/article/details/40899819
clear all
close all
clc
%% initialize the variables
x = 0.1; % initial actual state
x_N = 1; % 系统过程噪声的协方差  (由于是一维的，这里就是方差)
x_R = 1; % 测量的协方差
T = 75;  % 共进行75次
N = 100; % 粒子数，越大效果越好，计算量也越大    %initilize our initial, prior particle distribution as a gaussian around
%the true initial value    V = 2; %初始分布的方差
x_P = []; % 粒子
% 用一个高斯分布随机的产生初始的粒子
for i = 1:N    x_P(i) = x + sqrt(V) * randn;
end    z_out = [x^2 / 20 + sqrt(x_R) * randn];  %实际测量值
x_out = [x];  %the actual output vector for measurement values.
x_est = [x]; % time by time output of the particle filters estimate
x_est_out = [x_est]; % the vector of particle filter estimates.    for t = 1:T    x = 0.5*x + 25*x/(1 + x^2) + 8*cos(1.2*(t-1)) +  sqrt(x_N)*randn;    z = x^2/20 + sqrt(x_R)*randn;    for i = 1:N    %从先验p(x(k)|x(k-1))中采样    x_P_update(i) = 0.5*x_P(i) + 25*x_P(i)/(1 + x_P(i)^2) + 8*cos(1.2*(t-1)) + sqrt(x_N)*randn;    %计算采样粒子的值，为后面根据似然去计算权重做铺垫    z_update(i) = x_P_update(i)^2/20;    %对每个粒子计算其权重，这里假设量测噪声是高斯分布。所以 w = p(y|x)对应下面的计算公式    P_w(i) = (1/sqrt(2*pi*x_R)) * exp(-(z - z_update(i))^2/(2*x_R));    end    % 归一化.    P_w = P_w./sum(P_w);    %% Resampling这里没有用博客里之前说的histc函数，不过目的和效果是一样的    for i = 1 : N    x_P(i) = x_P_update(find(rand <= cumsum(P_w),1));   % 粒子权重大的将多得到后代    end                                                     % find( ,1) 返回第一个 符合前面条件的数的 下标    %状态估计，重采样以后，每个粒子的权重都变成了1/N    x_est = mean(x_P);    % Save data in arrays for later plotting    x_out = [x_out x];    z_out = [z_out z];    x_est_out = [x_est_out x_est];    end    t = 0:T;
figure(1);
clf
plot(t, x_out, '.-b', t, x_est_out, '-.r','linewidth',3);
set(gca,'FontSize',12); set(gcf,'Color','White');
xlabel('time step'); ylabel('flight position');
legend('True flight position', 'Particle filter estimate');

滤波后的结果如下：

这是粒子滤波的一个应用，还有一个目标跟踪(matlab)，机器人定位(python)的例子，我一并放入压缩文件，供大家下载，下载请点击。请原谅博主的这一点点自私。两个程序得效果如下：

粒子滤波从推导到应用这个系列到这里就结束了。结合前面几章的问题起来看，基本的粒子滤波里可改进的地方很多，正由于此才诞生了很多优化了的算法，而这篇博客只理顺了基本算法的思路，希望有帮到大家。
另外，个人感觉粒子滤波和遗传算法真是像极了。同时，如果你觉得这种用很多粒子来计算的方式效率低，在工程应用中不好接受，推荐看看无味卡尔曼滤波（UKF）,他是有选择的产生粒子，而不是盲目的随机产生。

参考：
1.particle filtering—粒子滤波（讲的很通俗易懂）

2. Taylan Cemgil 《A Tutorial Introduction to Monte Carlo methods, Markov Chain Monte Carlo and Particle Filtering》

3. M. Sanjeev Arulampalam 《A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking》

4. ZHE CHEN 《Bayesian Filtering: From Kalman Filters to Particle Filters, and Beyond》

5.百度文库《粒子滤波理论》

6. Haykin 《Neural Networks and learning Machines 》Chapter 14

7. 统计之都 <LDA-math-MCMC 和 Gibbs Sampling>

8. Gabriel A. Terejanu 《Tutorial on Monte Carlo Techniques》

Particle Filtering粒子滤波相关推荐

Particle Filter 粒子滤波
Particle Filter 粒子滤波贝叶斯滤波重要性采样 Sequential Importance Sampling (SIS) Filter 重采样 Sampling Importance ...
particle filter粒子滤波
粒子滤波可以用在非高斯分布下的线性和非线性系统.卡尔曼只能用在高斯线性系统,扩展卡尔曼和无迹卡尔曼可以用在高斯非线性系统.适应的场景:KF<-EKF<-UKF<-PF 基本步骤: 1 ...
particle filtering---粒子滤波
particle filtering---粒子滤波一.贝叶斯滤波二.蒙特卡洛采样三.重要性采样四.Sequential Importance Sampling (SIS) Filter 五.重 ...
粒子滤波（particle filtering）梳理
本文根据MLaPP第23章,从importance sampling开始梳理,直至导出粒子滤波,并给出相应的范例程序.rejection sampling太简单,就不介绍了. 1. Importanc ...
[转]粒子滤波（particle filtering）的思路发展过程及应用（详细深度好文）
粒子滤波作为视觉SLAM中后端进行状态估计的主要算法之一,很好的完成了扩展卡尔曼滤波无法有效处理的复杂状态方程下的状态估计任务.这篇文章详细地描述了粒子滤波的思想历程,即如何一步步从简单的状态估计.采 ...
Particle Filter Tutorial 粒子滤波：从推导到应用（二）
二.蒙特卡洛采样假设我们能从一个目标概率分布p(x)中采样到一系列的样本(粒子),(至于怎么生成服从p(x)分布的样本,这个问题先放一放),那么就能利用这些样本去估计这个分布的某些函数的期望值.譬如 ...
Particle Filter Tutorial 粒子滤波：从推导到应用（一）
前言: 早几年前写博客的时候,编辑公式都是在线网站编辑的,有时候这个网站不稳定,贴的公式看不到.如果有需要,pdf版下载地址点击打开链接,给您阅读带来的不便表示抱歉,祝大家好运. 博主在自主学习粒子滤 ...
高斯粒子滤波matlab,粒子滤波(Particle filter)matlab实现 | 学步园
粒子滤波是以贝叶斯推理和重要性采样为基本框架的.因此,想要掌握粒子滤波,对于上述两个基本内容必须有一个初步的了解.贝叶斯公式非常perfect,但是在实际问题中,由于变量维数很高,被积函数很难积分,常 ...
Particle Filter Tutorial 粒子滤波：从推导到应用（四）
六.Sampling Importance Resampling Filter (SIR) SIR滤波器很容易由前面的基本粒子滤波推导出来,只要对粒子的重要性概率密度函数做出特定的选择即可.在SIR中 ...
Particle Filter Tutorial 粒子滤波：从推导到应用（三）
五.重采样在应用SIS 滤波的过程中,存在一个退化的问题.就是经过几次迭代以后,很多粒子的权重都变得很小,可以忽略了,只有少数粒子的权重比较大.并且粒子权值的方差随着时间增大,状态空间中的有效粒子数 ...

Particle Filtering粒子滤波