赵安近况Ｖｏｌ．１０Ｎｏ．５Ｍａｙ．２００８
第１０卷第５期
２００８年５月
０引言
粒子滤波（Ｐａｒｔｉｃｌｅｆｉｌｔｅｒ，简称ＰＦ）是一种新 的滤波算法，是递推贝叶斯估计的一种具体实现方法。它通过非参数化的蒙特卡罗模拟方法来实现递推贝叶斯滤波，其思想是利用一系列随机抽取的样本以及样本的权重来计算状态的后验概率分布（样本即粒子），从而得到状态的各种信息。当样本点数增至无穷大，蒙特卡罗模拟特性与后验概率密度函数表示等价时，其滤波精度可以逼近最优估计。该方法不受系统线性和噪声高斯分布假定的限制，能适用任何环境下任何状态的转换或测量模型。 粒子滤波算法中普遍存在的问题是退化现象。退化现象是指经过几次迭代后，除了一个粒子具有较高权值以外，其余粒子均只有非常小的权值，进而使用于近似后验概率分布的有效粒子数太少。退化现象意味着大量的计算工作都被用来更新那些对后验概率密度的估计几乎没有作用的粒子上。目前解决退化问题主要有两种方法：选取合适的重要性函数和进行重采样。其中重采样的基本思想是去除权值小的粒子
而对权值大的粒子进行复制，但这样在采样结果中就会包含许多重复点，从而丧失了粒子的多样性，故可能使采样枯竭。增加粒子数可以改善这一问题，但仿真结果表明其效果并不理想。为此，本文讨论了粒子滤波算法，并在重采样后引入一个马尔可夫
链蒙特卡罗移动步骤（ＭＣＭＣ）增加粒子的多样性。该方法可用于粒子滤波器的各种改进算法中。与增加粒子数目的方法相比，ＭＣＭＣ可以更好的提高滤波器的性能。
１
粒子滤波
１．１
贝叶斯估计
贝叶斯估计为动态系统的估计提供了一个严
谨的解决框架。对于待估计的参数，贝叶斯估计先给出该参数的先验概率分布，然后再结合样本信息得到参数的后验概率分布情况。
首先可假定系统的状态方程是：
ｘｔ＝ｆｔ（ｘｔ－１，ｕｉ）
量测方程是：
ｚｔ＝ｈｔ（ｘｔ，ｖｔ）
其中，ｆｔ和ｈｔ分别是状态函数和观测函数；ｕｔ
和ｖｔ是后验概率已知的、独立于过去和当前的噪声序列，而且它们互相独立。若已知状态的初始概率密度函数为：ｐ（ｘ０｜ｚ０）＝ｐ（ｘ０），且假设系统服从一阶马尔可夫过程，且系统观测独立，则状态预测方程和状态更新方程分别为：
ｐ（ｘｔ｜ｚ１：ｔ－１）＝!
ｐ（ｘｔ｜ｘｔ－１）ｐ（ｘｔ－１｜ｚ１：ｔ－１）ｄｘｔ－１（１）ｐ（ｘｔ｜ｚ１：ｔ）＝
ｐ（ｚｔ｜ｘｔ）ｐ（ｘｔ｜ｚ１：ｔ－１）
!ｐ（ｚｔ
｜ｘｔ
）ｐ（ｘｔ
｜ｚ
１：ｔ－１
）ｄｘｔ
（２）
在更新方程（２）中，测量值ｚｔ被用来修正后验概率密度，以获得当前状态的后验概率密度函
收稿日期：２００７－１２－２１
一种有效的粒子滤波器的改进算法
薛亚茹，李明
阎锡山日记
（西安电子科技大学雷达信号重点实验室，陕西
西安
７１００７１）
摘
要：为了解决粒子滤波算法中重采样后粒子中包含重复点过多，从而丧失了粒子的多样
性等问题，文中在重采样后引入一个马尔可夫链蒙特卡罗（ＭＣＭＣ）移动步骤来增加粒子的多样性，因而能更好地近似状态的后验概率分布。
关键词：贝叶斯估计；粒子滤波；马尔可夫链蒙特卡罗（ＭＣＭＣ）步骤；Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ算法
技术前沿
ＥｌｅｃｔｒｏｎｉｃＣｏｍｐｏｎｅｎｔ＆ＤｅｖｉｃｅＡｐｐｌｉｃａｔｉｏｎｓＶｏｌ．１０Ｎｏ．５Ｍａｙ．２００８
第１０卷第５期
２００８年５月
数。式（１）和（２）是贝叶斯估计的一般表达式，其中的积分仅可对某些线性动态系统获得解析解，而非线性非高斯系统通常是很难求解的。１．２粒子滤波
贝叶斯估计是粒子滤波的理论基础。粒子滤波是基于贝叶斯估计和蒙特卡罗法近似数值解的方法，其本质是将积分运算变成有限样本点的求和运算。它的核心思想是利用一系列随机抽取的粒子和粒子的权重来计算状态的后验概率分布，从而得到状态的估计。
在蒙特卡洛法中，后验密度函数可由下式来近似：
ｐ（ｘ
０：ｔ｜ｚ
１：ｔ
）＝１
Ｎ
n0706
Ｎ
ｉ＝１
!δ（ｘ０：ｔ－ｘ（ｉ）０：ｔ）（３）
其中，ｘ（ｉ）
０：ｔ
是采样点，δ（）是Ｄｉｒａｃｄｅｌｔａ函数。
ｘ
０：ｔ
的任意函数的期望值：
Ｅ［ｇ（ｘ
０：ｔ
）
］＝"ｇ（ｘ０：ｔ）ｐ（ｘ０：ｔ｜ｚ０：ｔ）ｄｘ０：ｔ
均可由下式近似：
Ｅ（ｇ（ｘ
０：ｔ））＝１
Ｎ
Ｎ
ｉ＝１
!ｇ（ｘ（ｉ）０：ｔ）
且当Ｎ趋于无穷大时，有：
Ｅ［ｇ（ｘ
０：ｔ）］＝Ｅ（ｇ（ｘ
０：ｔ
））
通常情况下，难以直接从ｐ（ｘ
０：ｔ｜ｚ
１：ｔ
）抽样，故可
引入一个容易抽样且已知概率密度的重要性函数
ｑ（ｘ
０：ｔ｜ｚ
１：ｔ
）来代替ｐ（ｘ
０：ｔ
｜ｚ
１：ｔ
），则：
Ｅ［ｇ（ｘ
０：ｔ
）］＝"ｇ（ｘ０：ｔ）ｐ（ｘ０：ｔ｜ｚ０：ｔ）
ｑ（ｘ
０：ｔ
｜ｚ
０：ｔ
）
ｑ（ｘ
０：ｔ
｜ｚ
０：ｔ
）ｄｘ
０：ｔ
＝"ｇ（ｘ０：ｔ）ｐ（ｚ１：ｔ｜ｘ０：ｔ）ｐ（ｘ０：ｔ）
ｐ（ｚ
１：ｔ
）ｑ（ｘ
０：ｔ
｜ｚ
１：ｔ
）
ｑ（ｘ
０：ｔ
｜ｚ
１：ｔ
）ｄｘ
０：ｔ
＝"ｇ（ｘ０：ｔ）ｗｔ（ｘ０：ｔ）
ｐ（ｚ
１：ｔ
）
ｑ（ｘ
０：ｔ
｜ｚ
１：ｔ
）ｄｘ
０：ｔ
＝
"ｇ（ｘ０：ｔ）ｗ（ｘ０：ｔ）ｑ（ｘ０：ｔ｜ｚ１：ｔ）ｄｘ０：ｔ
"ｗ（ｘ０：ｔ）ｑ（ｘ０：ｔ｜ｚ１：ｔ）ｄｘ０：ｔ
＝
Ｎ
ｉ＝１
!ｇ（ｘ（ｉ）０：ｔ）ｗ#（ｘ（ｉ）０：ｔ）
其中，
ｗ
ｔ
＝
ｐ（ｚ
１：ｔ
｜ｘ
０：ｔ
）ｐ（ｘ
０：ｔ
）
ｑ（ｘ
牵引力控制系统０：ｔ
｜ｚ
１：ｔ
）
是没有归一化的重要性权值ｗ#
（ｉ）
＝ｗ（ｉ）／
Ｎ
ｊ＝１
!ｗ（ｊ）。因此，
后验概率密度可表示为：
ｐ（ｘ
０：ｔ
｜ｚ
１：ｔ
）＝
Ｎ
ｉ＝１
!ｗ#（ｉ）δ（ｘ０：ｔ－ｘ（ｉ）０：ｔ）
事实上，ｑ（ｘ
０：ｔ
｜ｚ
１：ｔ
）的选取应尽可能的接近系统
状态的后验概率分布，选取原则是使重要性权重
的方差最小。重要性权重的递归表达式如下：
ｗ
ｔ
＝
ｐ（ｚ
１：ｔ
｜ｘ
０：ｔ
）ｐ（ｘ
０：ｔ
）
ｑ（ｘ
０：ｔ
｜ｚ
１：ｔ
）
＝
ｐ（ｚ
１：ｔ
｜ｘ
０：ｔ
）ｐ（ｘ
０：ｔ
）
ｑ（ｘ
０：ｔ－１
｜ｚ
１：ｔ－１
）ｑ（ｘ
ｔ
｜ｘ
０：ｔ－１
，ｚ
匀速直线运动
１：ｔ
）
＝ｗ
ｔ－１
ｐ（ｚ
１：ｔ
｜ｘ
０：ｔ
）ｐ（ｘ
０：ｔ
）
ｐ（ｚ
１：ｔ－１
｜ｚ
０：ｔ－１
）ｐ（ｘ
０：ｔ－１
）
１
ｑ（ｘ
ｔ
｜ｘ
０：ｔ－１
，ｚ
１：ｔ
）
＝ｗ
ｔ－１
ｐ（ｚ
ｔ
｜ｘ
ｔ
）ｐ（ｘ
ｔ
｜ｘ
ｔ－１
）
ｑ（ｘ
ｔ
应变速率
｜ｘ
０：ｔ－１
，ｚ
１：ｔ
）
（４）
粒子滤波的算法步骤如下：
（１）从ｑ（ｘ
ｔ
｜ｘ（ｉ）
０：ｔ－１
，ｚ
１：ｔ
）中随机抽取Ｎ个粒子；
（２）根据式（４）计算各粒子的重要性权值并
将其归一化；
（３）计算有效采样尺度：
Ｎ
ｅｆｆ
＝Ｎ
１＋Ｖａｒ（ｗｉ
ｔ
）
当Ｎｅｆｆ＜Ｎｔｈ（有效粒子阈值）时，可对粒子进
行重采样并将权重重新设为１／Ｎ；
（４）根据：
ｐ（ｘ
０：ｔ
｜ｚ
１：ｔ
）＝１
Ｎ
Ｎ
ｉ＝１
!δ（ｘ０：ｔ－ｘ（ｉ）０：ｔ）
即可得到状态的均值和方差等各种统计信息。
２ＭＣＭＣ步骤
减小退化的一种方法是当退化现象出现时加
入重采样步骤，重采样就是对权值大的粒子进行
复制并去除权值小的粒子。其方法是对后验密度
函数的离散近似：
ｐ（ｘ
０：ｔ
｜ｚ
１：ｔ
）＝
Ｎ
ｉ＝１
!ｗ#（ｉ）δ（ｘ０：ｔ－ｘ（ｉ）０：ｔ）
再进行一次采样后，就会重新生成一个新的
粒子集。由于新的粒子独立分布，因此，重采样
后的粒子权值均为１／Ｎ。
图２２００次以上的仿真结果
重采样在一定程度上可以减小退化问题，但因此又会引来粒子耗尽问题。即经过几次迭代后，粒子中将包含许多重复点，这在估计过程中会由于粒子失去多样性而使目标概率的离散逼近变得不精确，从而影响粒子滤波器的估计性能。在动态噪声较小时，这一问题尤为严重。为此，本文在不影响近似有效性的前提下引入了一个算法步骤来增加粒子的多样性，即ＭＣＭＣ步骤。
ＭＣＭＣ的基本思想是：如果粒子的分布服从重要性函数ｑ（ｘ!０：ｔ｜ｚ１：ｔ），则可应用一个马尔可夫链转换核!（ｘ０：ｔ｜ｘ!０：ｔ）和后验概率密度函数ｐ（ｘ０：ｔ｜ｚ１：ｔ
），来使：
!!（ｘ０：ｔ
｜ｘ
!０：ｔ
）ｑ（ｘ!０：ｔ
｜ｚ１：ｔ
）ｄｘ!０：ｔ
＝ｐ（ｘ０：ｔ
｜ｚ１：ｔ
）
由此产生的新粒子仍服从状态的后验概率分布，但新的粒子更加逼近真实目标分布。
ＭＣＭＣ转移步骤实际就是从有限混合分布１ＮＮ
ｉ＝１"!
（ｘ０：ｔ－ｘ!（ｉ）
０：ｔ
）中进行采样来增加粒子的多样性。这里不要求转换核具有各态历经性。
ＭＣＭＣ转移有很多实现算法（如Ｇｉｂｂｓ采样，ＭｅｔｒｏｐｏｌｉｓＨａｓｔｉｎｇｓ算法等）。本文采用ＭｅｔｒｏｐｏｌｉｓＨａｓｔｉｎｇｓ算法来实现ＭＣＭＣ转移，其步骤如下：
（１）产生一个在［０，１］区间上服从均匀分布
的随机数ｖ；
（２）从重要性函数ｑ（ｘ!ｔ｜ｘ（ｉ）０：ｔ－１，ｚ１：ｔ
）中采样产生备用粒子ｘ＊（ｉ）
ｔ
，ｉ＝１，……，Ｎ；
（３）如果ｖ≤ｍｉｎ｛１，α｝，则接受ｘ＊（ｉ）
ｔ
，这时ｘ
（ｉ）
０：ｔ＝（ｘ!（ｉ）
０：ｔ－１，ｘ＊（ｉ）ｔ）；否则拒绝ｘ＊（ｉ）ｔ，ｘ（ｉ）０：ｔ
＝ｘ!（ｉ）
０：ｔ。其中α称为接受比率：
α＝
ｐ（ｚｔ｜ｘ＊（ｉ）ｔ）ｐ（ｘ＊
ｔ｜ｘ!（ｉ）
ｔ－１）ｑ（ｘ!ｔ｜ｘ!（ｉ）
０：ｔ－１，ｚ１：ｔ）
ｐ（ｚｔ｜ｘ!（ｉ）ｔ）ｐ（ｘ!（ｉ）ｔ｜ｘ!（ｉ）ｔ－１）ｑ（ｘ＊（ｉ）ｔ
｜ｘ!（ｉ）
０：ｔ－１，ｚ１：ｔ）上述算法可在重采样后进行。
３仿真结果分析
本文采用１００次蒙特卡罗方式来进行仿真。
其跟踪结果如图１所示，图１中，两种算法的粒子
数均为２００，图２中的ＰＦ－ＭＣＭＣ的粒子数为２００，
ＰＦ的粒子数为３００。
由仿真结果可以看出，ＰＦ算法在增加粒子数的情况下，其性能并没有得到较大的提高，滤波效果仍然不如ＰＦ－ＭＣＭＣ算法。由于ＭＣＭＣ步骤其实只是产生一部分新的粒子来增加粒子的多样性，所以延长了算法的运行时间，其实时性要差一些。
４结束语
粒子滤波是非线性贝叶斯递推算法，它适用
于各种非线性非高斯系统模型。它本身的优势使其越来越多的应用于各个领域，同时也出现了各种改进算法。本文将ＭＣＭＣ转移步骤应用到粒子滤波中来增加粒子的多样性，从而有效克服了粒子滤波重采样后产生的粒子多样性丧失问题，同时使粒子更接近状态的后验概率分布。仿真证明，该方法可用于粒子滤波器的各种改进算法中，以提高算法性能。
图１１００次蒙特卡罗的仿真跟踪结果
（下转第７７页）
参考文献
［１］Ｎ．Ｊ．Ｇｏｒｄｏｎ，Ｄ．Ｊ．ＳａｌｍｏｎｄａｎｄＡ．Ｆ．Ｍ．Ｓｍｉｔｈ，Ｎｏｖｅｌａｐ－ｐｒｏａｃｈｔｏｎｏｎｌｉｎｅａｒ／ｎｏｎ－ＧａｕｓｓｉａｎＢａｙｅｓｉａｎｓｔａｔｅｅｓｔｉｍａ－
ｔｉｏｎ．ＩＥＥＰｒｏｃｅｅｄｉｎｇｓ－Ｆ，ｖｏｌ．１４０，ｎｏ．２，１９９３，１０７－
１３３．［２］ＡｒｕｌａｍｐａｌａｍＳ，ＭａｓｋｅｌｌＳ，ＧｏｒｄｏｎＮａｎｄＣｌａｐｐＴ，ＡＴｕｔｏｒｉａｌｏｎＰａｒｔｉｃｌｅＦｉｌｔｅｒｓｆｏｒＯｎｌｉｎｅＮｏｎ２ｌｉｎｅａｒ／
Ｎｏｎ２ｇａｕｓｓｉａｎＢａｙｅｓｉａｎＴｒａｃｋｉｎｇ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎ
ｏｎＳｉｇｎａｌＰｒｏｃｅｓｓｉｎｇ，２００２，５０（２）：１７４－１８８．［３］Ｂｅｒｇｍａｎ．ＮａｎｄＤｏｕｃｅｔ．Ａ，ＭａｒｋｏｖＣｈａｉｎＭｏｎｔｅＣａｒｌｏＤａｔａＡｓｓｏｃｉａｔｉｏｎｆｏｒＴａｒｇｅｔＴｒａｃｋｉｎｇ［Ｊ］，Ｐｒｏｃ．ＩＥＥＥＣｏｎｆ．
Ａｃｏｕｓｔ．，Ｓｐｅｅｃｈ，ＳｉｇｎａｌＰｒｏｃｅｓｓ．，２０００，ｐｐ：７３５－７４２．［４］Ｗ．Ｇｉｌｋｓ，Ｓ．Ｒｉｃｈａｒｄｓｏｎ，ａｎｄＤ．Ｓｐｉｅｇｅｌｈａｌｔｅｒ，ＭａｒｋｏｖＣｈａｉｎＭｏｎｔｅＣａｒｌｏｉｎＰｒａｃｔｉｃｅ，ｅｄｓ．ＣｈａｐｍａｎａｎｄＨａｌｌ，
１９９６．
［５］Ａｎｄｒｉｅｕ，Ｃ．ｄｅＦｒｅｉｔａｓ，Ｊ．Ｆ．Ｇ．ａｎｄＤｏｕｃｅｔ．Ａ．，ＳｅｑｕｅｎｔｉａｌＭＣＭＣｆｏｒＢａｙｅｓｉａｎｍｏｄｅｌｓｅｌｅｃｔｉｏｎ，ＩＥＥＥＨｉｇｈｅｒＯｒｄｅｒ
ＳｔａｔｉｓｔｉｃｓＷｏｒｋｓｈｏｐ，Ｃｅａｓａｒｅａ，Ｉｓｒａｅｌ，１９９９，ｐｐ．１３０－１３４．
４三种算法的比较
图７所示是三种算法在加减速时间和最高速度给定情况下的轨迹。由图７可以看出，梯形速度曲线是一个恒加速过程，它的快速性很好，但它的加速度有突变的情况，在电机驱动元件性能不能达到理想的动态响应的条件下，实际使用梯形加减速的数控系统，其起停速度轨迹并不是理想的斜线，而是还存在明显的波动，其主要表现是起动时的滞后以及到达恒速时的过冲。
指数速度曲线在加速和减速阶段结束时的平滑性是最好的，而且它的快速性也比较好，但是相对于Ｓ曲线和梯形曲线，它在加减速启动过程中的速度、加速度变化都较大，这在负载扭矩较大时，容易造成电机的堵转；
相比之下，Ｓ曲线的加速度是一个连续的变化过程，它机械运动平稳性较好。虽然从速度的快速性上
比不上梯形曲线，但是在满足一定平稳性的条件下，Ｓ曲线可以将最大加速度提高或加长匀加速阶段，以提高加速效率。而且，完整的Ｓ曲线是由多个阶段组成的，加加速段、匀加速段以及匀速段的运动时间可以取不同值，因而可以减少很多情况下的功率耗散。
５结束语
经过本文的分析可见，梯形速度曲线的加速度有突变性，实际中运用的比较少，而指数速度曲线算法比较简单，跟踪性比较强，在跟踪响应要求较高的切削加工中应用比较多。Ｓ曲线算法相对较为精确，启动和停止的平稳性都比较好，可以根据各电机性能调整速度变化率。
事实上，三种加减速曲线在软件实现上都已经完成了数字化。但是为了保证算法的精度，应用时还有必要在软件上进行修正，特别是在加速到顶峰和减速到停止的时候。
参考文献
［１］李恩林．插补原理［Ｍ］．北京：机械工业出版社，１９８４．［２］丛爽，李泽湘．实用运动控制技术［Ｍ］．北京：电子工业出版社，２００６．
图６Ｓ曲线软件流程图
１．梯形曲线２．指数曲线３．Ｓ型曲线
图７三种曲线比较图
（上接第７３页）
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议