粒子滤波算法原理和仿真
1 引言
粒子滤波(Particle Filter, PF)是一种基于蒙特卡洛(Monte Carlo, MC)方法的递推贝叶斯滤波算法。其核心思想是通过从状态空间寻找的一系列随机样本来近似系统变量的概率密度函数，以样本均值代替积分运算，从而获得状态的最小方差估计。其中从状态空间中抽取的样本称为“粒子”。一般地，随着粒子数目的增加，粒子的概率密度函数就逐渐逼近状态的概率密度函数，从而达到最优贝叶斯估计的效果。
2 粒子滤波原理 2.1 系统的动态空间
对于被观测对象的状态，可以通过以下非线性离散系统来描述：
xt?f(xt?1,wt?1) (1)
zt?h(xt,vt) (2)
以上为系统的状态方程和观测方程。其中，f(?)为状态函数，h(?)为观测函数，xt是系统在时间t的状态变量，wt为对应的过程噪声，zt是系统在时间t的观测值，vt为对应的观测噪声。
从贝叶斯估计角度来看，状态估计问题就是根据观测信息z0:t构造状态的概率密度函数p(x0:t|z0:t)，从而估计在系统在任何状态下的滤波值。设系统状态序列函数为gt，则有：
E?gt(x0:t)???gt(x0:t)p(x0:tz0:t)dx0:x (3)
根据蒙特卡洛方法，后验概率分布可以用有限的离散样本来近似，由大数定律，当系统粒子数N→∞时，期望E[gt(x0:t)]可近似为：
1N(i)E?gt(x0:t)???gt(x0:t) (4)
Ni?1(i) 式中{x0:t: i=1,2,...N}为状态空间中按p(x0:t|z0:t)得到的采样点。
2.2 重要性采样
在粒子采集过程中，p(x0:t|z0:t)往往是未知且多变的，因此可先从一个已知且容易采样的参考分布q(x0:t|z0:t)中抽样，再通过对抽样粒子集进行加权求和来估计系统的状态值，即：
E?gt(x0:t)???gt(x0:t)p(x0:tz0:t)dx0:x??gt(x0:t)p(x0:tz0:t)q(x0:tz0:t)q(x0:tz0:t)dx0:t (5)
令ωt(x0:t) = p(z0:t|x0:t)p(x0:t) ∕ q(x0:t|z0:t)，则式(5)可表示为：
E?gt(x0:tg(x)w(x)q(x?)???w(x)q(xzt0:ttt0:t0:t0:t0:t0:tz0:t)dx0:x)dx0:x (6)
按照式(4)，式(6)可近似为：
E?gt(x0:t)??1N?g(xti?1N(i)0:t(i))?t(x0:t)(i)0:t1NNi?1??(xti?1N) (7)
(i)(i)??gt(x0:t)?t(x0:t)(i)(i) 其中?t(x0:t)为?t(x0:的归一化权值，x)0:t是由q( x0:t ?z0:t )采样获得的粒子。 t(i)2.3 序贯重要性抽样
贝叶斯估计是一个序列估计问题，通过序贯重要性抽样(Sequential Importance Sampling, SIS)建立采样粒子的序列关系。在t+1时刻采样时，不改变状态序列过去的样本集，而采用递归的形式计算重要性权值。因此参考分布可表示为：
q(x0:tz0:t)?q(x0:t?1z0:t?1)q(xtx0:t?1,z0:t) (8)
同时，假设系统状态是一阶马尔可夫过程，即xt只与xt-1相关，zt只与xt相关，因此：
p(x0:t)?p(x0:t?1)p(xtx0:t?1)?p(x0:t?1)p(xtxt?1) (9)
p(z0:tx0:t)?p(ztz0:t?1,x0:t)p(z0:t?1x0:t)?p(ztxt)p(z0:t?1x0:t?1) 再将式(8) ~ (10)代入ωt，可推知：
(10)
?t??t?1p(ztxt)p(xtxt?1)q(xtx0:t?1,z0:t) (11)
一般q( xt|x0:t-1 , z0:t ) =p(xt|xt-1, zt)为参考分布的最优分布，其概率密度仅依赖于xt-1和zt。但这种选择在实际中往往难以实现，因此多选用q( xt|x0:t-1 , z0:t ) = p(xt|xt-1)来近似。此时系统仅记录当前状态粒子xt(i)和其权重?t(i)，则由式(11)可知：
i)i?t(i)??t(?1p(ztxt) (12)
2.4 粒子滤波重采样
粒子退化是序贯重要性采样算法的主要问题，即经过若干次递推后，权值方差会逐渐增大，这样大量权值都集中在少数粒子上，而多数粒子则对系统状态的估计影响甚微，由此严重影响了滤波器的性能。重采样是解决粒子退化的重要方法，其核心思想是减少或剔除权值较小的粒子，而复制权值较大的粒子，最后将所有粒子的权值设为1/N。粒子滤波重采样方法的如下图所示：
图1 粒子滤波重采样示意
(i)(i) 图中xt(i)为采样得到的粒子，对应的权值为?t；xt为重采样所得的粒子，对应权值为1/N。目前常用的重采样算法包括随机重采样(random resampling)、多项式重采样(multinomial resampling)、系统重采样(systematic resampling)和残差重采样(residual resampling)，几种算法的计算步骤分别如下： 1. 随机重采样
设u~U(0,1]，{ui : i=1, 2, .., N }独立同分布。定义函数D(?)，若：
ui?(??,??t(j)](j)tj?1j?1(m)则D(ui) = m，且xt?xt。
(i)m?1mm?N (13)
2. 多项式重采样
设u~U(0,1]，{ui : i=1, 2, .., N }独立同分布，定义：
1??ui?ui?1uii?1?u?uNN?N(i)i?1,...,N?1 (14)
(m)对函数D(?)，由式(13)可知：D(ui) = m，且xt?xt。
3. 系统重采样 4. 残差重采样
综上所述，粒子滤波算法的一般流程为：
(i) （1）粒子集初始化，由先验概率密度p(x0)产生粒子集{x0}，并使所有粒子权值为1/N； (i)(i)(i) （2）按照xt~q(xtxt?1,zt)对t时刻的系统进行重要性采样，并通过当前观测值zt
根据式(12)计算粒子权值?t(i)；
（3）计算归一化权值?t，并按式(15)计算有效粒子数Neff；
(i)Neff?1?(?i?1N(i)2t (15)
) （4）若Neff小于设定阈值Nth，则对粒子进行重采样，得到新的粒子集{xt,?t}，其中?(i)t=1/N，否则
(i)(i){xt,?t}?{xt(i),?t}；
(i)(i)(i) （5）由式(7)推知系统的状态估计值xt为：
xt??xt?t (16)
i?1N(i)(i)再转到步骤（2），对t+1时刻的系统进行重要性采样。
3 粒子滤波仿真
首先构造系统的状态模型和观测模型：
xt?xk?125xt?1??8cos(1.2(t?1))?wt (17) 21?xt2?1xt2zk??vt (18)
20 式中，过程噪声wk ~ N(0, 5)，观测噪声vk ~ N(0,5)。设系统的状态初始值x0 = 0.1，初始分布p(x0) ~ N(0,2)。观测时长T = 200，每个时刻的采样粒子数N = 500。观测时间内，系统
状态值和观测值如图2所示：
图2 系统状态值和观测值
图中的理论状态值是系统在无过程噪声(wt = 0)的情况下，通过严格时序递推得到的理想
值，可以看出，它与系统的真实状态之间存在差异；理论观测值反映系统真实状态的理论观测情况(vt = 0)，而真实观测值是系统在观测噪声的干扰下，监测所得的实际数值。粒子滤波将通过观的真实值估计系统的状态，并期望滤波结果接近于真实值。对仿真系统的滤波结果如图3所示：
图3 系统的滤波效果
由图3知，粒子滤波很好地实现了对系统状态的估计。其中，为了进一步评估滤波的效果，将状态的估计值代入式(18)，并令vt = 0，得到滤波后的观测值。通过式(19)可以计算当前观测信号Z的信噪比SNR：
zzSNR?10?lg (19) T(z?Z)(z?Z) 式中z为观测信号的理论值，代表Z中真实有用的成分。表1给出几组不同滤波器参数下的观测信号的信噪比：
表1 不同参数滤波效果比较
观测信号 采样方差 SNR(dB) z — Tz1 1
z2 2 z3 5 z4 10 -1.239 3.992 4.052 4.747 3.303 表中z为滤波前的观测信号，z1~z4为按不同方差的分布进行粒子抽样，然后滤波所得的观测信号。由表1知，粒子滤波使信号的信噪比有了明显提升，而滤波效果与采样粒子的方差相关。根据上表及前文内容可知，当粒子的分布接近过程噪声时，往往可以得到较好的滤波效果。