粒子滤波：非线性和非高斯系统的状态估计
粒子滤波：非线性和非高斯系统的状态估计
什么是粒子滤波？
在工程、金融、机器人等领域，我们经常需要从一系列带有噪声的观测数据中，推断出系统内部真实的、不可直接测量的状态。例如，通过雷达信号估计飞机的位置与速度，或根据股票价格波动推断市场波动率。
经典的方法如卡尔曼滤波，在处理线性且噪声服从高斯分布的系统时表现完美。但现实世界充满了非线性与非高斯性：雷达测量角度与距离的关系是非线性的，金融收益率常常是尖峰厚尾的。粒子滤波正是为打破这些限制而生，它通过一群随机样本（粒子）来近似任意分布，灵活应对非线性和非高斯场景。
为何需要粒子滤波？
贝叶斯滤波框架
所有的状态估计问题都可以纳入贝叶斯滤波框架。我们有两个核心方程：
- 状态方程：$\mathbf{x}k = f(\mathbf{x}{k-1}, \mathbf{u}_k) + \mathbf{v}_k$
描述状态如何随时间演变，$f$ 可以是任意非线性函数，$\mathbf{v}_k$ 是任意分布的过程噪声。 - 观测方程：$\mathbf{z}_k = h(\mathbf{x}_k) + \mathbf{w}_k$
描述我们如何观测到状态，$h$ 可以是任意非线性函数，$\mathbf{w}_k$ 是任意分布的观测噪声。
目标是在获得所有历史观测 $\mathbf{z}_{1:k}$ 后，计算状态的后验概率密度 $p(\mathbf{x}k | \mathbf{z}{1:k})$。贝叶斯滤波给出了这个后验的递推公式，但对于一般的非线性和非高斯模型，递推中的积分通常没有解析解。粒子滤波使用蒙特卡洛方法突破了这一瓶颈。
核心思想：用粒子表示信念
粒子滤波将后验分布表示为一组带权重的随机样本（粒子）：
$$ p(\mathbf{x}k | \mathbf{z}{1:k}) \approx \sum_{i=1}^{N} w_k^{(i)} \delta(\mathbf{x}_k - \mathbf{x}_k^{(i)}) $$
其中每个粒子 $\mathbf{x}_k^{(i)}$ 是状态空间中的一个可能取值，$w_k^{(i)}$ 是其归一化权重（$\sum w_k^{(i)} = 1$），代表该粒子能解释观测数据的程度。粒子数 $N$ 越多，近似越精确。
序贯重要性采样
如何递推地更新粒子集？序贯重要性采样（SIS）提供了基本框架。
重要性采样原理
如果我们想从目标分布 $p(\mathbf{x})$ 中采样但难以直接做到，可以找一个容易采样的提议分布 $q(\mathbf{x})$，采样后赋予样本权重 $w \propto p(\mathbf{x})/q(\mathbf{x})$，这样加权样本即可表示 $p(\mathbf{x})$。
递推权重更新
在滤波场景下，我们需要递推采样。设当前时刻的提议分布可以分解为 $q(\mathbf{x}{0:k}|\mathbf{z}{1:k}) = q(\mathbf{x}k|\mathbf{x}{0:k-1},\mathbf{z}{1:k}) q(\mathbf{x}{0:k-1}|\mathbf{z}{1:k-1})$。从上一代的粒子 $\mathbf{x}{0:k-1}^{(i)}$ 出发，按 $q(\mathbf{x}k|\mathbf{x}{0:k-1}^{(i)},\mathbf{z}_{1:k})$ 延伸出新粒子 $\mathbf{x}_k^{(i)}$，权重更新遵循：
$$ w_k^{(i)} \propto w_{k-1}^{(i)} \frac{p(\mathbf{z}_k|\mathbf{x}k^{(i)}) p(\mathbf{x}k^{(i)}|\mathbf{x}{k-1}^{(i)})} {q(\mathbf{x}k^{(i)}|\mathbf{x}{0:k-1}^{(i)},\mathbf{z}{1:k})} $$
最简且最常用的提议分布是状态转移先验：$q(\mathbf{x}k|\mathbf{x}{0:k-1},\mathbf{z}{1:k}) = p(\mathbf{x}k|\mathbf{x}{k-1})$。此时权重更新简化为 $w_k^{(i)} \propto w{k-1}^{(i)} p(\mathbf{z}_k|\mathbf{x}_k^{(i)})$，即粒子只需要通过状态方程进行预测，然后根据新观测的似然调整权重。
序贯重要性重采样
单纯的SIS面临权重退化问题：经过几次迭代，绝大多数粒子的权重会趋近于零，只有极少数粒子有效，造成大量计算浪费。
有效样本量
衡量退化程度的指标是有效样本量（Effective Sample Size, ESS），近似为：
$$ \hat{N}{\text{eff}} = \frac{1}{\sum{i=1}^N (w_k^{(i)})^2} $$
当 $\hat{N}_{\text{eff}}$ 低于某一阈值（如 $N/2$）时，执行重采样。
重采样操作
重采样从当前粒子集合中，依据权重有放回地抽取 $N$ 次，权重大的粒子被多次复制，权重小的粒子被淘汰，并将所有权重重置为 $1/N$。常用重采样方法：
- 多项式重采样：直接按权重抽样，简单但方差大。
- 系统重采样：将 $[0,1)$ 分成 $N$ 个等宽区间，每个区间随机偏移后联合抽样，方差更小。
- 残差重采样：保留粒子的确定性复制部分，其余随机，兼具效率与低方差。
重采样会牺牲粒子多样性，但保证了粒子集持续聚焦在高概率区域。
标准粒子滤波算法
结合SIS与重采样，得到最通用的采样重要性重采样粒子滤波：
- 初始化：$k=0$，从先验 $p(\mathbf{x}_0)$ 采样 $N$ 个粒子 $\mathbf{x}_0^{(i)}$，权重 $w_0^{(i)} = 1/N$。
- 对于每个时刻 $k=1,2,\dots$：
- 预测：对每个粒子 $i$，从提议分布采样新粒子
$\mathbf{x}k^{(i)} \sim q(\mathbf{x}k | \mathbf{x}{0:k-1}^{(i)}, \mathbf{z}{1:k})$
（通常使用状态转移先验：$\mathbf{x}_k^{(i)} \sim p(\mathbf{x}k | \mathbf{x}{k-1}^{(i)})$） - 更新权重：按公式计算非归一化权重
$\tilde{w}k^{(i)} = w{k-1}^{(i)} \frac{p(\mathbf{z}_k|\mathbf{x}k^{(i)}) p(\mathbf{x}k^{(i)}|\mathbf{x}{k-1}^{(i)})} {q(\mathbf{x}k^{(i)}|\mathbf{x}{0:k-1}^{(i)},\mathbf{z}{1:k})}$
若用先验提议，简化为 $\tilde{w}k^{(i)} = w{k-1}^{(i)} p(\mathbf{z}_k|\mathbf{x}_k^{(i)})$。
归一化：$w_k^{(i)} = \tilde{w}k^{(i)} / \sum{j=1}^N \tilde{w}_k^{(j)}$。 - 计算有效样本量：若 $\hat{N}{\text{eff}} < N{\text{th}}$，执行重采样：
从 ${\mathbf{x}_k^{(i)}}$ 中以权重 ${w_k^{(i)}}$ 重采样 $N$ 次，得到新粒子集，置所有 $w_k^{(i)} = 1/N$。
- 预测：对每个粒子 $i$，从提议分布采样新粒子
- 状态估计：输出后验均值 $\hat{\mathbf{x}}k = \sum{i=1}^N w_k^{(i)} \mathbf{x}_k^{(i)}$ 或最大后验等。
进阶变体
为提升性能，出现了多种改进算法：
辅助粒子滤波
在重采样之前，先通过一个辅助变量评估“当前粒子在未来时步的表现”，使高似然的粒子更容易存活，进一步减缓退化。尤其适用于观测噪声较小、似然函数尖锐的场景。
正则化粒子滤波
重采样后从连续分布中平滑粒子，避免完全相同的粒子副本。例如用核密度估计的带宽对重采样粒子添加微小噪声，维持多样性并改善连续状态空间的估计精度。
应用实例
目标跟踪
雷达测量目标距离 $r$ 和方位角 $\theta$，测量方程为： $$ z_k = \begin{bmatrix} \sqrt{x_k^2 + y_k^2} \ \arctan(y_k/x_k) \end{bmatrix} + \mathbf{w}_k $$ 非线性严重，噪声可能非高斯。粒子滤波直接基于此模型计算粒子似然，稳健估计目标位置。
机器人定位
移动机器人在已知地图中，通过激光或视觉特征更新位姿。观测函数体现了高度非线性与多峰分布（走廊环境下的位置歧义）。粒子滤波能同时维持多个位姿假设，是蒙特卡洛定位的基石。
金融波动率建模
资产收益率的波动率是隐变量，观测方程通常是线性的，但状态方程中的波动率动态常为非线性且噪声服从厚尾分布。粒子滤波能灵活捕捉波动率的突变性和聚集性。
粒子滤波的优势与挑战
优势：
- 处理任意非线性、非高斯状态空间模型。
- 全局逼近，可表示多峰分布。
- 实现直观，易于扩展。
挑战与对策：
- 粒子退化与样本贫化：重采样可缓解退化，但过度重采样降低多样性。引入马尔可夫链蒙特卡洛移动步骤或正则化可改善。
- 高维状态下的维度灾难：所需粒子数随维度指数增长。此时需结合边缘化、Rao-Blackwellization等技术降维。
- 提议分布选择：使用最优提议分布 $p(\mathbf{x}k|\mathbf{x}{k-1},\mathbf{z}_k)$ 能显著提高效率，但在非线性下常需近似。
结语
粒子滤波是贝叶斯状态估计的强大工具，它以纯粹的蒙特卡洛方法突破了传统滤波的线性高斯牢笼。理解其“预测-更新-重采样”的核心流程，便能搭建出解决众多实际非线性估计问题的基础框架。当你下次面对难以线性化或噪声非高斯的动态系统时，粒子滤波或许就是最灵活的答案。