第2章 粒子滤波技术及其收敛性分析证明
第2章 粒子滤波技术及其收敛性分析证明
粒子滤波(particle filter)，又称序贯蒙特卡罗方法，是一种基于蒙特卡罗方法和递推贝叶斯估计的统计滤波方法，它依据大数定理采用蒙特卡罗方法来求解贝叶斯估计中的积分运算。其基本思想是：首先依据系统状态向量的经验条件分布在状态空间产生一组随机样本的集合，这些样本称为粒子；然后根据量测量不断地调整粒子的权重和位置，通过调整后的粒子的信息，修正最初的经验条件分布。其核心思想是：用由粒子及其权重组成的离散随机测度近似相关的概率分布，并且根据算法递推更新离散随机测度。当样本容量很大时，这种蒙特卡罗描述就近似于状态变量真实的后验概率密度函数。这种技术适用于任何能用状态空间模型以及传统的卡尔曼滤波表示的非高斯背景的非线性随机系统，精度可以逼近最优估计，是一种很有效的非线性滤波技术。
2.1滤波问题常用框架
2.1.1 滤波常用框架
滤波的概念其实是结合在线得到的量测值对系统状态(参数或隐变量)进行序贯估计的问题。
序贯估计问题介绍：
1、后验概率(posterior density)：p?Xk|Zk?。其中Xk??x1,x2,,xk?为系统状态变量，
Zk??z1,z2,,zk?为系统量测值。
2、递推计算得到边缘后验概率，即滤波密度，p?xk|Zk?，这里就省却了保留存储所有的历史状态。
3、在获得滤波密度的前提下，系统状态的一系列估计值可以得到：
1)均值(系统状态的最小均方差估计)：
?k?E?xk|Zk???xkp?xk|Zk?dxkx
2)模式、中值、置信度、峰态等。
zk?2zk?1zkp?zk|xk?observedunobservedxk?2p?xk|xk?1?xk?1xk
图 2.1 滤波的常规框架
Fig 2.1 the conventional framework of filtering
7
江苏科技大学工学硕士学位论文
2.1.2 动态空间模型
时变问题分析的一般动态状态空间模型分为状态转移模型p?xk|xk?1?和量测模型；zk?Ry为系p?zk|xk?，其中xk?Rnx代表系统在时间k的状态变量（隐含变量或参数）
n统在时间k的量测值。对于非线性非高斯过程，其模型可表示为
xk?f?xk?1,vk?1? (2.1)
zk?h?x? (2.2) k,nk其中zk?Ry表示输出量测值，xk?Rnx为系统的状态变量；vk?Rnv和nk?Rnn分别为过程噪声和量测噪声，vk?1和nk是相互独立、协方差分别为Q?和R?k的零均值加性噪声序
kn列；f:R?R，h:R?R为有界非线性映射。
为了使状态方程和量测方程能描述更为一般的动态系统，这里作以下假设： 1、状态符合一阶马尔可夫过程。
nxnxnnnnp?xk|xk?1,xk?2,2、量测相对于状态是独立的。
x0??p?xk|xk?1?
p?zk|xk,A??p?zk|xk?
2.1.3 递推贝叶斯估计
根据给定的状态空间模型，可以递推得到如下滤波概率密度：
p?xk|Zk????p?Zk|xk?p?xk?p?Yk?p?yk,Zk?1|xk?p?xk?p?zk,Yk?1?p?yk|Zk?1,xk?p?Zk?1|xk?p?xk?p?zk|Yk?1?p?Zk?1?p?zk|Zk?1,xk?p?xk|Zk?1?p?Zk?1?p?xk??p?zk|Zk?1?p?Zk?1?p?xk??p?zk|xk?p?xk|Zk?1?p?zk|Zk?1? (2.3)
即
p?xk|Zk??p?zk|xk?p?xk|Zk?1?likehood?priorlikehood?prior??p?zk|Zk?1?evidencenormalizing constant
其中，p?zk|xk?为似然函数，p?xk|Zk?1?为先验概率，转移概率为p?xk|xk?1?，似然函数由量测模型得出，先验概率由转移概率和上一相邻时刻的滤波密度p?x率可由状态方程得到。p?zkk?1|Zk?1?推得。而转移概
|Zk?1?为归一化常数，可由量测模型和似然函数决定，即
p?zk|Zk?1???p?zk|xk?p?xk|Zk?1?dxk。
8
第2章 粒子滤波技术及其收敛性分析证明
2.2 粒子滤波理论
粒子滤波提供了一种方便而有效的从非高斯、非线性、高维的量测数据中计算后验概率的方法，非常灵活，容易实现，可以并行化，实用性强，在工程学、应用统计学等领域获得了最新发展。下面介绍粒子滤波基本原理[1-2,21,25-26]。 2.2.1 标准粒子滤波算法
本章采用(2.1)和(2.2)中描述的状态空间模型来描述粒子滤波算法。 假设动态系统的状态空间模型如下： 系统状态模型为xk?fk?xk?1,vk?1?，如(2.1)，系统量测模型为zk?hk?xk,uk?，如(2.2)。
为了计算在k时刻的后验概率密度p?x0:k|z1:k?，有两个步骤：预测和更新。
1、预测：假设k?1时刻p?x0:k?1|z1:k?1?已知。先验概率密度p?xk|xk?1?由(2.1)得到。预测公式如下
p?x0:k|z1:k?1?=?p?xk|xk?1?p?x0:k?1|z1:k?1?dxk?1 (2.4)
2、更新：由k时刻的量测值z，根据贝叶斯规则对预测进行更新。更新公式如下
kp?x0:k|z1:k?=p(zk|xk)p?xk|xk?1??p(zk|xk)p?xk|xk?1?dxk (2.5)
式中的积分通常难以获得解析解，通常基于蒙特卡罗方法转化为有限样本点的求和运
iiNi,N}是支持样本算。假设{x0:k,wk}i?1表示后验概率为p?x0:k|z1:k?的随机粒子集合，其中{x0:k,i?1,集，相应的权值为{wki,i?1,,N}，将权值归一化处理后，后验概率可以近似为
iip?x0:k|z1:k???wk??x0:k?x0:k? (2.6)
i?1N其中
ip(x0:k|z1:k)w?iq(x0:k|z1:k)ikiiiip(zk|xk)p(xk|xk?1)p(x0:k?1|z1:k?1)?iiiq(xk|x0:k?1,z1:k)q(x0:k?1,z1:k?1)i?wk?1p(zk|x)p(x|x)iiq(xk|x0:k?1,z1:k)ikikik?1 (2.7)
q??为重要密度函数，通常，重要性密度函数选取为转移先验，且假设(2.1)符合一阶
9
江苏科技大学工学硕士学位论文
马尔可夫过程，则有
iiq?xk|x0:k?1,z1:k??q?xk|xk?1,zk?
iiq?xk|xk?1,zk??p?xk|xk?1?
代入(2.7)得
iiiwk?wkpz|x??1kk? (2.8)
若在对量测精度要求低的场合，这种选取方法能够获得较好的结果。但是由于这种选取方法没有考虑当前的量测值，从重要性概率密度中取样得到的样本与从真实后验概率密度采样得到的样本有很大的偏差。尤其当似然函数位于系统状态转移概率密度的尾部或似然函数呈尖峰状态时，这种偏差就更加明显。似然函数与先验分布关系的具体情况如下图所示：
p(zk|xk)p?xk|xk?1?
图2.2 似然函数呈尖峰状态
Fig 2.2 Likelihood function in peak condition
p(zk|xk)p?xk|xk?1?
图2.3 似然函数位于先验分布的尾部
Fig 2.3 Likelihood function at the tail of the prior distribution
2.2.2 粒子集的退化和重采样
SIS存在一个很严重的缺陷，总要性权重的方差随着时间随机递增，使得粒子的权重
iii集中到少数粒子上，这是因为wk?p?x0:k|z1:k?/q?x0:k|z1:k?，当用重要性函数替代后验概率分布作为采样函数，理想情况是重要性函数非常接近后验概率分布，也就是希望重要性函数的方差基本为零。即
iii?varq(?|z1:k)?px|z/qx|z?varw?????0:k1:k?q(?|z1:k)k??0 ?0:k1:k因此，重要性权重方差的增长给采样的准确性带来很坏的影响，它经常使得粒子的权
10