粒子滤波的基本原理粒子滤波算法广泛应用在视觉跟踪领域、通信与信号处理领域、机器人、图像处理、金融经济、以及目标定位、导航、跟踪领域,其本质是利用当前和过去的观测量来估计未知量的当前值。在粒子滤波算法中使用了大量随机样本,采用蒙特卡洛仿真来完成递推贝叶斯滤波过程,其核心是使用一组具有相应权值的随机样本(粒子)来表示状态的后验分布。该方法的基本思路是选取一个重要性概率密度函数并从中进行随机抽样,得到一些带有相应权值的随机样本后,在状态观测的基础上调节权值的大小和粒子的位置,再使用这些样本来逼近状态后验分布,最后通过这组样本的加权求和作为状态的估计值。粒子滤波不受系统模型的线性和高斯假设约束,采用样本形式而不是函数形式对状态概率密度进行描述,使其不需要对状态变量的概率分布作过多的约束,适用于任意非线性非高斯动态系统,是目前最适合于非线性、非高斯系统状态的滤波方法【ArulampalamMS,MaskellS,GordonN,etal.Atutorialonparticlefiltersforonlinenonlinear/non-GaussianBayesiantracking[J].IEEETransactionsonSignalProcessing,2002,50(2):174-188.】1动态系统的状态空间模型状态空间模型包括系统状态方程和观测方程,其通用的表示方法分别为【梁军.粒子滤波算法及其应用研究[D].哈尔滨工业大学,2009.】【黄小平,王岩,廖鹏程.粒子滤波原理及应用——MATLAB仿真[M].电子工业出版社.2017】1,kkkfXXW(1),kkkhZXV(2)其中f和h为已知函数,kW和kV是概率密度已知的随机变量,kX代表k时刻的状态量,kZ代表k时刻的观测量,kW和kV是相互独立的。关于系统的状态方程和观测方程,通常也可用1kkpXX表示状态转移模型;kkpZX表示观测似然模型;0pX表示初始状态的先验分布;0:1:kkpXZ表示系统的后验密度;1:kkpXZ表示边沿后验密度,或称为后验滤波密度。卡尔曼滤波以及粒子滤波算法的本质即是利用观测序列1:kZ对当前状态进行优化,从而得到k时刻的后验滤波密度,进而得到k时刻的状态值。【MerweRVD,DoucetA,FreitasND,etal.Theunscentedparticlefilter[C]//InternationalConferenceonNeuralInformationProcessingSystems.MITPress,2000:563-569.】2贝叶斯滤波原理贝叶斯估计是利用先验知识和实际观测值构造系统的后验概率来估计系统的状态量。若系统的先验概率密度为0:kpX,观测值为1:kZ,则系统的后验概率密度可表示为1:0:0:1:0:1:0:0:0:kkkkkkkkkpppppdZXXZXZXXX(3)从式(3)可以看出,后验概率实际上是利用观测值对先验概率进行修正。由于进行了修正,后验概率密度较先验概率密度更加接近被估计值的真实概率密度。滤波问题是为了计算后验滤波概率密度1:kkpXZ,它实际上是0:1:kkpXZ的边沿概率密度即1:0:1:01......kkkkkppddXZXZXX(4)当每一个新的观测量kZ到来时,即需要完成式(4)得到后验滤波概率密度。为清晰的给出后验滤波概率密度的逐级递推形式。首先,通过预测步骤得到不包含k时刻观测值的k时刻系统的先验滤波概率密度1:1111:11kkkkkkkpppdXZXXXZX(5)其中1kkpXX由系统状态方程决定,11:1kkpXZ是1k时刻的后验滤波概率密度。然后,根据贝叶斯公式,利用k时刻的先验滤波概率密度1:1kkpXZ推导k时刻的后验概率密度的表示形式1:11:1:1kkkkkkkkppppZXXZXZZZ(6)其中kkpZX由观测方程决定。1:11:1kkkkkkkpppdZZZXXZX(7)一般为归一化常数。因此贝叶斯滤波的递推框图如下所示0pX初始状态状态转移先验滤波密度10pXZ观测似然密度11pZX密度修正后验概率密度11pXZ状态转移先验滤波密度21pXZ观测似然密度22pZX密度修正后验概率密度22pXZ1k时刻2k时刻图1贝叶斯递推滤波框图贝叶斯递推过程中的主要的几个公式1:1111:11kkkkkkkpppdXZXXXZX(利用上一时刻的后验滤波密度,结合状态转移概率对这一时刻的先验滤波概率密度进行求解)1:11:1kkkkkkkpppdZZZXXZX(利用先验滤波概率密度以及观测似然密度得到观测概率)1:11:1:1kkkkkkkkppppZXXZXZZZ(利用这一时刻的先验滤波概率密度,结合观测似然概率密度对这一时刻的后验滤波概率密度进行求解)以上三个式子是贝叶斯滤波递推求解后验滤波概率密度的核心公式,称为贝叶斯递推滤波算法。这里给出贝叶斯滤波原理的一个简单例子。一个简单的例子:假设动态系统状态空间模型中状态方程和观测方程分别为1kkkXXW(8)kkkZXV(9)其中~0,1kNW,~0,0.5kNV,并给出初始状态00X。由此,可以得到:0X是以概率1取值为0,即00pXX(10)状态转移概率2111exp22kkkkpXXXX(11)观测似然函数密度21exp20.5kkkkpZXZX(12)该系统等效于某信号经过一个0,1N的高斯信道进行传输,然后再利用一个测量误差为0,0.5N的测量设备对该信号进行测量。当从0时刻向1时刻转移时1010001000211exp22pppdpdXZXXXXXXXXX(13)10111012211112211112112112111=expexp20.5221=exp20.5222-1113=exp-exp1331223311exp-33pppddddZZZXXZXZXXXZXXXXZZXZ(14)11101110221112121111expexp20.52211exp-33213exp112233ppppZXXZXZZZZXXZXZ(15)此时1k时刻的后验滤波概率密度11121~,33pNXZZ(16)通过式(16)可以看出,1k时刻的后验概率密度是利用观测方程对状态方程得到的先验概率密度的修正,因此较状态方程和观测方程的误差都要小。然后将1k时刻的后验概率密度作为2k时刻的先验概率密度进行递推,即从1时刻向2时刻转移过程时21211112211211212122212112221212113expexp1221223312331124expexp1862222243311exp862223pppdddXZXXXZXXZXXXXXZXZXZXXZXZ22123exp42223XZ(17)也就是说21124~,33pNXZZ(这是利用1Z对2X的估计)2122212222222121222122221122222113311expexp20.586222412112334411expexp41133111122211334exp113311pppdddZZZXXZXZXXZXZXXZZZZZZXZZ122214112333exp111126ZZZZ(18)222121:221222222121222112221213311expexp20.5862223344exp1133111141211211exp422211ppppZXXZXZZZZXXZXZZZZZXZZ(19)因此,2k时刻的后验滤波概率密度21:212414~2,11211pNXZZZ(20)这样依次递推,即可得到每一时刻2X的后验概率密度,从而实现对状态量的后验滤波。3贝叶斯重要性采样和序贯重要性采样3.1贝叶斯重要性采样在实际问题中,求出后验概率密度后,需要根据后验密度得到某个函数的估计值。然而,后验概率密度函数与其他函数乘积的积分往往很难求解。根据蒙特卡洛仿真原理,只需要按照后验概率分布对状态量进行抽样,将每一个抽样得到的状态量称为一个“粒子”,然后将每个粒子带入函数,求函数的均值。该均值可以近似看作函数的估计值。0:0:11NkkiEggiNXX(21)从上式中可以看出,当N越大该近似值越接近0:kEgX。然而,从后验密度中直接抽样是难以实现的。于是,引入重要性概率密度的参考分布:1:kokqXZ,该参考分布应该是便于抽样的,利用重要性概率密度和贝叶斯公式综合运算可得*0:0:10:*0:10:0:111NkkkikNkkiNkkkigiiNEgiNgiiXwXXwXXwX(22)其中1:0:0:*0:0:1:kkkkkkkppiqZXXwXXZ,*0:0:*0:1kkkkNkkiiiiwXwXwX(23)这样,引入重要性概率密度的参考分布后,将其转化为每个粒子的权值0:kkiwX乘入到粒子的函数值中进行加权即可。同样,注意到当参考分布与后验密度相同时,加权系数*0:1:kkkipwXZ,这时候加权系数*0:kkiwX与抽样值0:kiX没有关系是一个常数。因此0:1kkiNwX。3.2序贯重要性采样(SIS算法)上述分析虽然解决了后验概率密度难以进行采样的问题,但每当一个新的观测数据到来时,就需要重新从重要性概率密度中抽取样本