从贝叶斯滤波到粒子滤波

本文转载于知乎：从贝叶斯滤波到粒子滤波^[1] (文末直达) ，作者：石鹏^[2]

注意：由于公式序号较多，在公众号会出现排版问题，所以小编将公序号进行了剔除已经删除最后应用实例部分。如果对本文公式推导等感兴趣，可以参见原文。该系列还有几篇高质量文章，感兴趣可以关注这位博主的推文。

小编有话说

目录

0 回顾与引言

在此前的文章：从概率到贝叶斯滤波^[3]中，我们曾经提到，贝叶斯滤波框架下求解预测步中的先验概率密度函数 、更新步中的归一化常数 、状态量最优估计 时均涉及到无穷积分，大多数情况无法得到解析解，使得贝叶斯滤波算法的直接应用十分困难。为了解决贝叶斯滤波积分难的问题，通常从两个角度出发：

(1) 作理想假设

• 线性高斯假设状态转移函数 和观测函数 均为线性函数，过程噪声随机变量 和观测噪声随机变量 均服从均值为 0 的正态分布。线性高斯问题可通过卡尔曼滤波（Kalman Filter）进行解决，详细内容请参考此前文章

• 非线性高斯假设状态转移函数 和（或）观测函数 为非线性函数，过程噪声随机变量 和观测噪声随机变量 均服从均值为 0 的正态分布。非线性高斯问题可通过扩展卡尔曼滤波（Extended Kalman Filter）和无迹卡尔曼滤波（Unscented Kalman Filter）进行解决。扩展卡尔曼滤波通过一阶泰勒级数展开将非线性高斯系统的状态转移函数 和（或）观测函数 线性化，然后采用标准卡尔曼滤波框架实现状态量的滤波过程，详细内容请参考此前文章

无迹卡尔曼滤波基于无迹变换，无迹变换研究的是如何通过确定的采样点捕获经非线性变换的高斯随机变量的后验分布的问题，通过无迹变换得到相应的统计特性后，再结合标准卡尔曼滤波框架，便得到无迹卡尔曼滤波，详细内容请参考此前文章

(2) 化连续为离散

对于强非线性非高斯问题，需要将无穷积分转化为离散的数值积分——即化整为零，由此可以引申出本文的主题——粒子滤波（Particle Filter，PF）。

粒子滤波是贝叶斯滤波的一种非参数实现，所谓非参数，即不对滤波状态量的后验概率密度作任何假设。粒子滤波的主要思想是用一系列从后验得到的带权重的随机采样表示后验。从采样的角度考虑，粒子滤波与无迹卡尔曼滤波相似，区别在于，无迹卡尔曼滤波使用 sigma 确定性采样，通过无迹变换计算 sigma 样本点的位置与权重；而粒子滤波使用蒙特卡罗随机采样从建议分布中得到样本（粒子），并通过观测值更新粒子权重，针对粒子的权值退化问题，还涉及粒子的重采样步骤。粒子滤波算法广泛用于解决无人车的定位问题。

关键词： 贝叶斯，大数定律，蒙特卡罗，重要性采样，序贯重要性采样，粒子退化，重采样，采样重要性重采样，无人车定位

1 蒙特卡罗方法

假设存在某一连续型随机变量 ，其概率密度函数为 ，则 X 的数学期望为：

若存在另一连续型随机变量 ，满足 ，则 的数学期望为：

正如我们所知道的，上式中的 和 可能很难得到解析解。蒙特卡罗（Monte Carlo）方法告诉我们，可以通过对随机变量的概率密度进行随机采样，并对样本进行加权求和来近似随机变量的期望，如此一来，积分问题便转化为有限样本点的求和问题。

假设 是从 中进行随机采样得到的样本集合，样本数量为 N，每一个样本（粒子） 代表 的一种可能状态，即

根据辛钦大数定律， 的样本均值依概率 1 收敛于数学期望，即对 ，有

式 (1.4) 和 (1.5) 意味着，当样本数量 N 足够大时，有

2 重要性采样

2.1 什么是重要性采样

如上述，如果我们能够对概率密度函数 进行随机采样，便可通过对样本进行加权求和来近似随机变量 和 的数学期望，然而现实情况是， 可能很难直接采样，甚至根本无法采样。此时，可以选择从一个更加容易采样的概率密度函数 中进行随机采样得到样本集合 ，并通过 近似估计

其中

上述的采样方式称为重要性采样（Importance Sampling），其中的更易采样的概率密度函数 称为建议分布（Proposal Distribution）, 也称重要性密度（Importance Density）或重要性函数 (Importance Function) 。

显然

故, 式 (2.1) 可改写为

我们记 , 并记

我们称之为粒子 的非归一化的重要性权重 (Unnormalized Importance Weight) , 我们称之为归一化的重要性权重（Normalized Importance Weight）, 则式 (2.5) 等价于

从式 (2.7) 中我们不难发现, 归一化的重要性权重 即对应粒子 的离散概率值 (概率质 量), 只要得到了粒子 及其对应的归一化的重要性权重 , 便可通过式 (2.7) 近似估计期 望 。简言之, 重要性采样要解决的是原分布难以采样甚至无法采样的问题。

2.2 序贯重要性采样

上述的重要性采样过程针对的是单一随机变量的估计, 而对于贝叶斯估计而言, 我们需要处理的是 从 0 时刻到 时刻的随机过程

意即, 我们需要估计的是条件期望

类似 节中的推导过程, 我们从一个更加容易采样的概率密度函数 中进行随机采样 得到样本集合 , 并通过 近似估计条件期望

其中

是由每一时刻的第 个样本粒子依次增广构成的粒子簇

令

显然

故, 式 (2.10) 可改写为

我们记 , 并记

我们称之为粒子簇 的非归一化的重要性权重, 我们称之为归一化的重要性权重, 则 式 (2.15) 等价于

从式 (2.17) 中我们不难发现, 归一化的重要性权重 即对应粒子簇 的离散概率值（概率 质量) , 只要得到了粒子簇 及其对应的归一化的重要性权重 , 便可通过式 (2.17) 近似估 计条件期望 。 式 (2.16) 中的重要性权重是批量形式, 存在这样一个问题, 每个新的时刻的观测数据到来时都需 要重新批量计算整个状态序列的重要性权重, 随着时间的推移, 这无疑将带来巨大的资源消耗。为 了递推地更新相邻时刻间粒子的重要性权重, 我们引入序贯重要性采样 (Sequential Importance Sampling, SIS）算法。 根据贝叶斯公式、马尔可夫状态独立性假设和马尔可夫观测独立性假设, 我们很容易得到原分布的递推形式

假定建议分布有如下的递推形式

将式 (2.18) 和式 (2.19) 代入式 (2.13)

根据式 (2.19) 对建议分布递推形式的假设, 式 (2.11) 粒子簇的批量采样可以很自然地转化为单个粒 子的序贯采样形式:

将式 (2.21) 代入式 (2.20) 便得到重要性权重的递推更新形式:

若我们只关心当前 时刻状态的估计结果, 则根据系统的马尔可夫特性, 有

此时的建议分布只与 时刻的状态 和 时刻的观测结果 有关, 这意味着, 在实际应 用时, 我们无需存储更早时刻的历史状态序列和观测序列。根据式 (2.23), 式 (2.21) 演变为

将式 (2.23) 和式 (2.24) 代入式 (2.22), 可以得到重要性权重最终的递推更新形式

注意，式 (2.25) 中的 。按照式 (2.25) 递推地得到 时刻每个粒子的重要性权重后，再对权重作归一化处理，最终便可根据式 (2.17) 近似估计条件期望。式 (2.24)、(2.25)、(2.16) 和 (2.17) 构成了 SIS 算法的核心流程，下面我们总结下 SIS 的整体框架伪代码。

简言之，序贯重要性采样是为了序贯地进行粒子采样，并递推地更新粒子的重要性权重。

3 粒子退化问题

3.1 概述

3.1.1 什么是粒子退化

SIS 算法在经历次多次迭代后，粒子重要性权重的方差可能将变得很大，从而引发粒子退化问题（Particle Degeneracy Problem）。所谓粒子退化，指的是大量粒子中只有少数粒子具有较高权重，而绝大多数粒子的权重都很小甚至接近于 0，导致计算加权均值时大量的运算资源被浪费在了小权重粒子上。粒子退化问题发生的根本原因是建议分布与真实分布的不匹配。

3.1.2 粒子退化程度的度量

为了度量粒子退化问题的严重程度，我们引入 有效样本数 ，并将其定义如下

的直接计算并不容易, 我们通常代为使用其估计量

越小, 则意味着粒子退化越严重。抑制粒子退化问题最直观的方法是增加采样粒子数量 , 但这无疑会增加计算负担, 降低算法实时性, 以下两个手段则更具实操性:

• 选择合适的建议分布
• 在 SIS 算法流程结束后, 实施粒子重采样 (Resampling) 步骤

3.2 建议分布的选择

最优的建议分布无疑就是真实的后验分布本身

但这个结论对于问题的求解并没有实际意义, 因为正是由于后验分布末知或难于积分我们才引入了 建议分布的概念。建议分布的设计已经发展出很多形式, 并由此衍生出种类繁多的粒子滤波变种:

• 以 UKF 生成建议分布的无迹粒子滤波 (Unscented Particle Filter, UPF)
• 以 EKF 生成建议分布的扩展卡尔曼粒子滤波（Extended Kalman Particle Filter, EKPF）
• 通过构建辅助变量, 提升和观测更为匹配的粒子被采样的概率的辅助粒子滤波（Auxiliary Particle Filter, APF)

UPF、EKPF、APF 等并不是本文的主题, 这里我们介绍一种精度尚可, 但实现更加简单的建议分 布形式一一先验状态转移概率分布。

根据式 (2.19) 和式 (2.23) 我们不难发现, 只要条件建议分布 一经确定, 建议分布 便自然确定, 假设我们选取系统的先验状态转移概率分布 作为条件建议 分布

根据式 (3.4), 式 (2.24) 演变为

从式 (3.5) 我们可以知道, 只要确定了系统先验的状态转移概率分布, 便可以根据系统的状态转移 方程进行粒子序贯采样。根据式 (3.4) 和式 (3.5), 式 (2.25) 演变为

从式 (3.6) 我们可以知道, 只要完成了粒子的序贯采样, 并确定了系统的似然概率分布, 便可以根 据系统的观测方程以及上一时刻粒子的重要性权重递推地更新当前时刻粒子的重要性权重（非归一化的）。

3.3 重采样

重采样也可有效抑制粒子退化问题。所谓重采样, 指的是在得到当前时刻的粒子集及每个粒子归一 化的重要性权重 后, 根据每个粒子的权重需要进行重新采样, 粒子权重越高, 被重新采样到的概率也越高, 这意味着, 某些粒子在重采样后可能会被复制多份, 而某些粒子在重 采样后可能直接不存在了, 如下图所示 (参考 23, P263)。 Resampling 重采样示意图 重采样过程等价于对下面这个离散分布进行采样

我们记重采样得到的粒子集及每个粒子各自的权重为 , 因为是直接对离散分布 进行采样, 故

针对重采样过程, 目前已发展出多种方法, 如多项式重采样（Multinomial Resampling）、分层重 采样（Stratified Resampling）、系统重采样（Systematic Resampling）、残差重采样（Residual Resampling）等, 后面我们会对这几种方法作简单介绍并给出详尽的伪代码实现框架。

重采样也可有效抑制粒子退化问题。所谓重采样，指的是在得到当前时刻的粒子集及每个粒子归一化的重要性权重 后, 根据每个粒子的权重需要进行重新采样, 粒子权重越高, 被重新采样到的概率也越高, 这意味着, 某些粒子在重采样后可能会被复制多份, 而某些粒子在重 采样后可能直接不存在了, 如下图所示 (参考 23, P263)。

重采样过程等价于对下面这个离散分布进行采样

我们记重采样得到的粒子集及每个粒子各自的权重为 , 因为是直接对离散分布 进行采样, 故

针对重采样过程, 目前已发展出多种方法, 如多项式重采样（Multinomial Resampling）、分层重 采样（Stratified Resampling）、系统重采样（Systematic Resampling）、残差重采样（Residual Resampling）等, 后面我们会对这几种方法作简单介绍并给出详尽的伪代码实现框架。

重采样也可有效抑制粒子退化问题。所谓重采样，指的是在得到当前时刻的粒子集及每个粒子归一化的重要性权重 后, 根据每个粒子的权重需要进行重新采样, 粒子权重越高, 被重新采样到的概率也越高, 这意味着, 某些粒子在重采样后可能会被复制多份, 而某些粒子在重 采样后可能直接不存在了, 如下图所示 (参考 23, P263)。

3.3.1 几种经典的重采样算法

3.3.1.1 多项式重采样

又称简单随机重采样（Simple Random Resampling）。其核心思想很简单, 首先, 计算粒子归一 化权重的累积分布

然后, 生成服从均匀分布 内的随机数

使用二分搜索查找该随机数在累积分布中所处的位置, 得到重采样粒子对应的原粒子的索引 i, 多项式重采样算法的时间复杂度为 , 其中的对数项源于二分搜索的使用。下面我们 给出多项式重采样的伪代码。

其中, 表示粒子总数, 表示重采样的粒子数 。在单独使用多项式重采样时, ; 在用于残差重采样的步骤二时, $\left.N_{r}<n\right),\lfloor\cdot\rfloor$ 表示向下取整。在计算粒子权重的累="" 积分布时,="" 为避免圆整误差,="" 需要令="" $c_{n-1}="1$" 。此外,="" 使用二分搜索时需要注意搜索区间的设="" 置。<="" p="">

3.3.1.2 分层重采样

分层重采样将粒子归一化重要性权重的累积分布等分为 层子区间，对每层子区间，使用一个介于子区间左右边界值的随机数来查找该层子区间重采样粒子所对应的原粒子索引，分层重采样算法的时间复杂度为 。下面我们给出分层重采样的伪代码。

其中， 表示粒子总数， 表示重采样的粒子数（。在单独使用分层重采样时，；在用于残差重采样的步骤二时，）。

3.3.1.3 系统重采样

系统重采样与分层重采样类似，唯一的区别在于，系统重采样中的每层子区间在查找各自重采样粒子对应的原粒子索引时使用的是同一个随机数，意即，系统重采样只涉及一次随机数生成操作，系统重采样算法的时间复杂度同样为 。下面我们给出系统重采样的伪代码。

其中， 表示粒子总数， 表示重采样的粒子数（。在单独使用系统重采样时，；在用于残差重采样的步骤二时，）。

3.3.1.4 残差重采样

残差重采样分两步进行：

步骤一：确定性拷贝采样

对原粒子集的每个粒子 根据 作为复制次数进行拷贝采样, 其中

例如, 原粒子集共有 个粒子, 第 个粒子 的归一化的重要性权重为 , 则在残差重采样的步骤一中, 将被拷贝 次。在步骤一中, 总的重采样粒子数为

每个拷贝后的重采样粒子的权重置为 。

步骤二：残差随机采样

显然, 由于忽略了 的小数部分, 经过步骤一后重采样粒子集中还缺少 个粒子, 对于剩余的待重采样的粒子, 我们结合上面的重采样方法进行随机采样。首先, 计算原粒 子集中每个粒子经过确定性拷贝采样后新的归一化权重 (这里我们称之为归一化的残差权 重, 上角标 代表 residual)

然后, 将 和 代入多项式重采样、分层重采样或系统重采样中的任意一种, 完成剩余重 采样粒子的采样。 下面我们给出残差重采样的伪代码, 其中步骤二使用的是多项式重采样。

3.3.2 重采样的副作用

如果粒子的权值退化问题非常严重，那么重采样后的粒子将会是极少数个别粒子的大量副本，这意味粒子的多样性严重丧失，此时的粒子集已经无法很好地刻画原本的概率密度函数，我们称这种现象为粒子的样本贫化（Sample Impoverishment）。样本贫化极有可能导致滤波器发散，为了处理样本贫化问题，已经发展出很多方法，例如正则粒子滤波（Regularized Particle Filter，RPF），本文不对此作展开，有兴趣可以查阅相关文献。

3.3.3 基本的粒子滤波器

将 SIS 与重采样进行结合便构成了基本的粒子滤波器，下面我们总结下整体框架的伪代码（以 SIS + SystematicResampling 为例）。

其中， 是用于决定是否执行重采样步骤的阈值，通常取 。

4 SIR 滤波器

这里我们介绍一种常用的经典粒子滤波算法实现——采样重要性重采样（Sampling Importance Resampling，SIR）滤波。

在 3.2 节的结论中我们已经知道，假设我们取系统的先验状态转移概率分布 作为条件建议分布，便可以根据系统的状态转移方程进行粒子的序贯采样；同时，只要完成了粒子的序贯采样，并确定了系统的似然概率分布，便可以根据系统的观测方程以及上一时刻粒子的重要性权重递推地更新当前时刻粒子的重要性权重（非归一化的）。据此，我们对 3.3.3 节中的基本粒子滤波算法框架中的建议分布稍加修改，并在每次计算完归一化的重要性权重后都实施一次重采样，便可得到 SIR 滤波算法框架：

从框架中我们可以看出，SIR 的粒子序贯采样依赖上一时刻的粒子状态及系统的状态转移方程，我们称之为预测步；SIR 的粒子重要性权重序贯更新依赖上一时刻的粒子重要性权重及预测步粒子对应的似然概率，我们称之为更新步。预测步只改变粒子的位置，不改变粒子的权重；更新步只改变粒子的权重，不改变粒子的位置。SIR 算法的流程示意图如下图所示（出自参考 24，P681）。

6 总结

首先，我们回顾了用于解决线性高斯问题的卡尔曼滤波，用于解决非线性高斯问题的扩展卡尔曼滤波与无迹卡尔曼滤波，为解决强非线性非高斯问题我们引入了粒子滤波的概念。

粒子滤波基于随机采样实现。我们从蒙特卡罗方法讲起，为解决原分布难以采样甚至无法采样的问题引入了重要性采样。为序贯地进行粒子采样，并递推地更新粒子的重要性权重，又引入了序贯重要性采样。序贯重要性采样存在粒子退化问题，为应对该问题，可以选择恰当的建议分布并实施重采样步骤，文中介绍了四种重采样方法。通过将系统的状态转移概率分布选择为建议分布，并辅以重采样方法，我们由序贯重要性采样得到了采样重要性重采样滤波器。

参考资料