高斯混合采样粒子滤波算法

高斯混合采样粒子滤波算法
贝叶斯方法为动态系统的估计问题提供了一类严谨的解决框架。

它利用已知的信息建立系统的概率密度函数可以得到对系统状态估计的最优解。

对于线性高斯的估计问题，期望的概率密度函数仍是高斯分布，它的分布特性可用均值和方差来描述。

卡尔曼滤波器很好地解决了这类估计问题[1]。

对于非线性系统的估计问题，最经典并得到广泛应用的方法以扩展的卡尔曼滤波为代表，这类方法需要对模型进行线性化，同时要求期望的概率密度函数满足高斯分布，然而在对实际系统建模时，模型往往是非线性非高斯的。

此时，最优估计很难实现。

粒子（particle）滤波器——序列重要性采样粒子滤波器，是一种适用于强非线性、无高斯约束的基于模拟的统计滤波器。

它利用一定数量的粒子来表示随机变量的后验概率分布，从而可以近似得到任意函数的数学期望，并且能应用于任意非线性随机系统。

本文介绍一种估计性能更好的粒子滤波算法——高斯混合采样粒子滤波器(GMSPPF)，相比通常意义上的粒子滤波算法（SIR-PF），GMSPPF粒子滤波器具有更小的系统状态估计的均方误差和均值。

贝叶斯滤波问题
贝叶斯滤波用概率统计的方法从已观察到的数据中获得动态状态空间（DSS）模型参数。

在DSS模型中，包含状态和观测两个方程。

其中状态转移方程（State Equation）通常写作
（1）这里，是已知，且是白噪声独立的随机序列，而且分布是已知的。

观测方程表达式写为
（2）这里：是白噪声序列，独立且分布已知。

并且满足。

图1描述了DSS模型中状态转移和似然函数的关系。

假设初始时刻系统的状态分布已知，k时刻的已知信息序列表示。

图1 动态状态空间模型（DSSM）
这样，贝叶斯估计的问题理解为：利用观测到的信息Yk，求解系统状态的概率分布。

若系统状态的变化是隐马尔柯夫过程，即当前系统的状态信息只与上一个时刻的状态有关，可以通过预测和更新的途径求解。

（3）这里：
（4）假设xk，wk是相互独立的随机变量，满足。

于是，参考（1）式可以把（4）式写为
（5）其中，是采样函数。

当是已知时，xk可以通过确定性方程（1）得到。

依据贝叶斯准则，系统状态估计量
（6）其中，
（7）另外，在给定xk，vk，分布的条件下，yk的条件概率依据测量方程（2）可以表示为如下形式
（8）由（6）式可以看出，后验概率密度包含3个部分。

先验概率似然函数和证据。

如何获得这三项的近似是贝叶斯滤波的核心问题。

更新方程（5）中观测值用来对的先验预测值修正，从而获得状态的后验概率。

方程（3）和（6）的递归关系构成了求解贝叶斯估计问题
的两个步骤：预测与更新。

如果(1)，(2)中的hk，fk是线性的，且噪声wk，vk满足高斯白噪声，可以把贝叶斯估计问题简化为卡尔曼分析解。

但这类问题仅仅是实际问题中很小的一个部分。

对于更多的问题，很难得到分析解。

只有通过对问题的近似线性处理(扩展卡尔曼滤波)或其它途径（蒙特卡洛方法）实现非线性、非高斯问题的解。

依据后面分析问题需要，这里重点对蒙特卡洛方法积分进行说明。

蒙特卡洛方法
在过去的二十多年，蒙特卡洛方法得到了很大的发展。

其优点就是用系列满足条件的采样点及其权重来表示后验概率密度。

蒙特卡洛方法采用统计抽样和估计对数学问题进行求解。

按照其用途，可以把蒙特卡洛方法分为三类：蒙特卡洛抽样、计算、优化。

其中，蒙特卡洛抽样是寻找有效的、方差很小的、用于估计的抽样方法。

蒙特卡洛计算则是设计产生满足特定要求随机数的随机发生器的问题。

而蒙特卡洛优化是采用蒙特卡洛思想对实际中的非凸非差分函数优化求解。

对于，可以由概率空间p(x)中抽取N个样本，用近似值作为的解。

大数定理证明：收敛于，并且满足条件。

这里，是的方差。

不同于确定性的数字计算，蒙特卡洛近似的一个重要特点就是估计的精度独立于状态空间的维数。

而且，积分估计的方差与采样点的个数成反比。

显然，蒙特卡洛近似方法的关键点有两个：首先如何由一个样本空间中抽取N个采样点，用来表征后验概率密度。

其次就是计算。

重要性抽样（Important Sampling）解决了如何借助于已知分布来对实现有效采样的问题，由Marshall 1965年提出。

当数据空间十分巨大时，重要性抽样只对其中“重要”区域进行采样，节省了计算量。

对于高维采样空间模型，如统计物理学、贝叶斯统计量，这一点尤为重要。

重要性抽样的中心思想是选择一个覆盖真实分布p（x）的建议分布q（x）。

这样，
(9) 对q（x）作蒙特卡洛抽样，假设粒子数目为N，有
(10) 其中，称为重要性权重，再作归一处理，
(11) 是归一化权重。

为了减小估计的方差，选择的建议性分布q（x）与p（x）尽可能匹配。

通常，建议分布q（x）需要一个长的拖尾，这样可以解决区间之外的干扰。

确切的说，匹配的q（x）必须与p（x）f（x）成正比。

当q（x）与p（x）不匹配时，w（x(i)）是不均匀分布的，在整个递归迭代的过程中，存在大量的权值极小的样本，而这些样本对估计的贡献很小。

事实上，权值较大的少数样本决定蒙特卡洛采样的估计精度。

大量时间损耗在这些“无关紧要”的粒子计算上，即所谓的粒子退化现象（Degeneracy Problem）。

目前，标准的粒子滤波器选择先验概率（Prior）作为建议分布。

对于粒子退化现象，采样—重要性重采样方法给出了很好的解决途径。

其基本思想就是通过在两次重要性采样之间增加重采样步骤，消除权值较小的样本，并对权值较大的样本复制，降低了计算的复杂度。

在o(N)时间复杂度范围内可以已排序的均匀分布序列作重采样处理。

对重采样（Resampling）处理，新的采样结果放在数组，具体的算法用伪码语言写为如下的形式：
步骤1：令这里必须注意是随机变量的累计概率密度序列。

步骤2：初始假设，当，产生一组序列分布。

对一个固定的j，分别用逐一比较，一旦，就可以得到一组新的样本集合。

如此循环直到。

需要说明的是，重采样方法在消除粒子退化问题的
同时，也带来了其它两个问题：首先，降低了粒子运算并行执行的可能性；其次，由于权值较大的粒子多次被选择，粒子的多样性减少。

这种情况尤其在小过程噪声条件下表现更为明
显[11]。

图2 SIR-PF重要性采样与重采样示意图
GMSPPF滤波算法
如前所述，利用序列重要性采样和重采样的方法，粒子滤波可以有效的递归更新后验概率的分布。

但是，由于对粒子未加假设，大量的粒子在处理非线性、非高斯问题时出现了计算的高复杂性问题。

另外，由于少数权值较大的粒子反复被选择，粒子坍塌明显。

文献提出了在重要性采样步骤的建议分布的生成阶段“搬运”粒子到似然较高区域，可以缓解坍塌，同时提高估计的性能。

但是不可避免的是对每一个粒子的后验概率处理，使得计算的复杂性进一步加剧。

鉴于此种情况，这里介绍一种新颖的高斯混合采样粒子滤波器（Gaussian Mixture Sigma Point Particle Filter，GMSPPF）。

GMSPPF算法利用有限高斯混合模型表征后验概率分布情况，可以通过基于重要性采样的加权的后验粒子，借助于加权的期望最大化算法（Weighted Expection Maximization）替换标准重采样步骤，降低粒子坍塌效应。

基于高斯混合近似的采样卡尔曼滤波器
根据最优滤波理论，一个概率密度p(x)都可以写作高斯混合模型（Gaussian Mixture Model）。

即，这里，G是高斯分量的个数，是高斯分量的权重，是以向量为均值，以p(g)为协方差矩阵的随机向量x的高斯分布。

考虑DSS状态转移方程和观测方程，假设先验概率及噪声密度服从高斯混合模型（GMM）。

这样，预测的先验概率密度满足，更新后。

这里，。

在此基础之上，预测的先验概率和后验概率对应的均值和方差可以通过采样卡尔曼滤波器（Sigma Point KF）计算。

基于观测更新的重要性采样（Important Sampling）
前已叙及重要性抽样是一种蒙特卡洛方法，即用一组带有权值的样本数据来表征随机变量的概率密度。

利用DSS模型的一阶马尔柯夫本质和给定状态的观测值依赖性，可以推导递归的权值更新方程，这里仅对于给定的粒子而言。

在GMSPPF算法中，用GMM近似来。

作为建议分布。

由于包含了最新的样本数据，使得粒子聚集在高似然区域，一定程度减少了粒子坍塌效应。

另外，使用预测的先验概率平滑权值更新方程中的，这是因为GMSPPF算法用GMM表示后验概率，本次后验同时又是下一个时间步的先验概率，GMM模型中高斯核对后验概率做了平滑处理。

基于观测更新步骤的重要性采样方法中对粒子不作任何假设，对非线性、非高斯问题具有很强的鲁棒性。

采用加权的EM算法做重采样和GMM还原
基于观测更新步骤的重要性采样输出是一组加权的粒子，在标准的粒子滤波器中，这些粒子必须作重采样处理丢弃小权值粒子，同时对权值较大的粒子做放大处理。

通过这种处理，可以有效的防止粒子集合的方差增加太快。

不幸的是，重采样步骤只对当观测似然微弱、大量粒子聚集极少数粒子副本情况有效。

在GMSPPF算法中，采用加权的期望最大（Weighted Expection Maximization）直接得到GMM模型，实现对加权粒子的最大似然拟合，这就相当于对粒子的后验概率做了平滑，避免了粒子坍塌问题，同时，GMM 模型中的高斯核的个数减少到G，防止其呈指数级增长，降低了算法复杂度。

为了比较算法的性能，系统状态估计的条件均值，均方误差（Error Convariane）可以通过两个方法计算，即在加权的EM算法平滑之前，用下面公式
求解，描述了系统的均值与均方误差性能。