粒子滤波理论

粒子滤波理论模拟方法来实现递推贝叶斯滤波，适用、目标跟踪、本章首先概述用于求解目标状态后验概率的贝叶多样性丧失问题，提出了一种量子进化粒子滤波算法。贝叶斯滤波动态系统的目标跟踪问题可以通过图2.1所示的状态空间模型来描述。本节在贝叶斯滤波框架下讨论目标跟踪问题。图2.1 状态空间模型Fig.2.1 Statespacemodel在目标跟踪问题中，动态系统的状态空间模型可描述为xkf(xk1)uk1ykh(xk)vk

(2.1)f(h(xkykuk为过程vkXkx0:kx0x1,⋯xk与Yky1:ky1,⋯,yk}xkxk-1有关。另外ykkxk有关。贝式求解后验概率密度pXk|Yk或滤波概率密度pxk|Yk)。新过程则利用最新的测量值对先验概率密度进行修正，得到后验概率密度。xk1|Yk1，贝叶斯滤波的具体过程如下：p(xk,xk1|Yk1)p(xk|xk1,Yk1)p(xk1|Yk1)

(2.2)xk1xk与Yk1相互独立，因此

(2.3)xk1Chapman-Komolgorov方程(2.4)pxk|Yk1pxk|Yk：kyk后，利用贝叶斯公式对先验概率密度进行更新，得到后验概率p(x

xk,Yk1)p(xk|Yk1)k k p(y

k|Yk1)\*MERGEFORMAT(2.5)ykxk决定，即)p(yk|xk)因此p(x

|Yk1)k k p(y

k|Yk1)\*MERGEFORMAT(2.7)yk|Yk1为归一化常数

|Yk1)p(yk|xk)p(xk|Yk1)dxk(MAP)准则或最小均方误差(MMSE)准则，将具有极大后验概率密度的状态或条件均值作为系统状态的估计值，即ˆAP

argminp(x

|Y)ˆˆ[f(x|Y]

f(f(x)p(x

k k kxk\*MERGEFORMAT(2.9)阶或二阶Taylor级数展开的扩展Kalman滤波得到广泛应用[119]滤中方案便是基于蒙特卡洛模拟的粒子滤波器。粒子滤波2050年代，Hammersley便采用基于序贯重要性采样(Sequentialimportance2060年代后期，HandschinMayne122]2070Hachnkahi以及Zaritskii等学者的一系列研究工作使得序贯蒙特卡洛方法得到进一步发展[123][126]。限卡洛方法成功解决了一系列高维积分问题[127-130]。SmithGelfand提出的采样-重采样思想Bayesian推理提供了一种易于实现的计算策略[131]。随后，SmithGordon2090年代初将重采样(Resampling)步骤引入到粒子滤波中，在一定程度上解决了序贯重要性采样的权值退化问题，并由此产生了第一个可实现的SIR(SamplingimportanceNodestar21世纪，粒子滤波器成为[133-135]，IEEE出版的论文集“SequentialMonteCarloMethodsinPractice”对粒子滤波器进行了详细介绍[136]。贝叶斯重要性采样蒙特卡洛模拟是一诺特卡洛模拟方法利用所求状态空间中大量的样本点来近似逼近待估计变量的后验概率分布，px|YNx(ii1,⋯N，则有k k k1N (i)N(xkxk)Ni1

(2.11)xk为连续变量(x-xk为单位冲激函数(狄拉克函数)，即x-xk)0,xxk，且(x)dx1。xkP(xk|Yk可近似逼近为P(x

1

x(i))k k N

i1

k k(2.12)

x(i))1,x

x(i);

x(i))0,x

x(i)。k k k k k k k k1图2.2 经验概率分布函数Fig.2.2 Empiricalprobabilitydistributionfunctionkkx(i)为从后验概率密度函数pxkk估计可以用求和方式逼近，即

f(xk的期望1N (i)Nf(xk)Ni1

(2.13)构造概率模型。对于本身具有随机性质的问题，主要工作是正确地描述和模拟这个概用蒙特卡洛方法求解需要事先构造一个人为的概率过程，将它的某些参量视为问题的解。从指定概率分布中采样。产生服从己知概率分布的随机变量是实现蒙特卡洛方法模拟试验的关键步骤。建立各种估计量的估计。一般说来，构造出概率模型并能从中抽样后，便可进行现模拟试验。随后，就要确定一个随机变量，将其作为待求解问题的解进行估计。qxk|Yk，从中生成采样粒子，利用这些随机样本的加权和来逼近后验滤波概率密度px|Y)，如图2.3所示。令{x(iw(ii1,⋯.N}x(i)为k k k k kkk时刻第iw(i，则后验滤波概率密度可以表示为kp(x

|Y)w(i)(x

x(i))Nk k k k N(2.14)其中，wwk

p(x(i)|Y) k kk q(x(i)|Yk

(2.15)图2.3 重要性采样Fig.2.3 Importancesampling当采样粒子的数目很大时，式(2.14)便可近似逼近真实的后验概率密度函数。任意函数的期望估计为1N (i)

p(x(i)|Y) 1N

(i) (i)

NN

) k kk q(x(i)|Yk

NN

)wk(2.16)序贯重要性采样算法在基于重要性采样的蒙特卡洛模拟方法中，估计后验滤波概率需要利用所有的观测数x0:k|y1:k)q(x0:k1|y1:k1)q(xk|x0:k1,y1:k)

k(2.18)k(2.19)后验概率密度函数的递归形式可以表示为p(x

0:k k

yk|Yk1)p(yk|x0:k,Yk1)p(xk|x0:k1,Yk1)p(x0:k1|Yk1)xk)p(xk|xk1)p(x0:k1|Yk1)p(yk|Yk1)(2.20)k粒子权值w(i)的递归形式可以表示为kwwk

p(x(i)|Y) 0:k k0:k q(x(i)|Y0:k p(y|x(i))p(x(i)|x(i))p(x(i) |Y ) k k k 1 1q(x(i)|x(i)

,Y)q(x(i) |Y )k k

k1(i)

|x(i))p(x(i)|x(i))k k k 1

(2.21)q(x(i)|x(i)

,Y)k k通常，需要对粒子权值进行归一化处理，即(i)

w(i)̃k

k Nkw(i)Nki1

(2.22)计。序贯重要性采样算法的流程可以用如下伪代码描述：[{x(i),w(i)}N

]SIS({x(i),w(i)}N,Y)k k Fori=1:N

1 i1 kq(x(i)|x(i),Y)x(i)；k 0:k1 k kkw(i；kEndFor几乎不起作用的粒子更新上，使得估计性能下降。通常采用有Neff来衡量粒子权值的退化程度，即eff N N/(1var(w*(i)eff

(2.23)*(i)

p(x(i)|y)wk k 1:k

(2.24)q(x |q(x |x )k k1 1:k有效粒子数越小，表明权值退化越严重。Neff可以近似为1ˆefN1

(2.25)k(w(i))2ki1在进行序贯重要性采样时，若ˆef小于事先设定的某一阈值，则应当采取一些措施加(1)选择合适的重要性概率密度函数；(2)在序贯重要性采样之后，采用重采样方法。重要密度函数的选择pxk|xk1作为重要性概率密度函数。此时，粒子的权值为k w(i)w(i)k

|x(i))

k 虑最新观测数据所提供的信息，从中抽取的样本与真实后验分布产生的样本存在一定的偏k k {w(k

Doucet等给出的最优重要性概率密度函数为q(x(i)|x(i),y)p(x(i)|x(i),y)k 1 k k 1 kp(y

|x(i),x(i))p(x(i)|x(i)) k k 1 k 1|x p|x

(i)k p(y

|x(i))p(x(i)|x(i)) k k k k1p(y

|x(i))k k1

(2.27)此时，粒子的权值为k w(i)w(i)k

(i)|x k |x (2.28)以px(i)|x(i),y作为重要性概率密度函数需要对其直接采样。x为有限离k 1 k k散状态或px(i)|x(i),ypy|x(i))才存在解析解。在实际情况中，构造k 1 k k k1从2.4所示。图2.4 移动粒子至高似然区域Fig.2.4 Movethesamplesinthepriortoregionsofhighlikelihoodkk1时刻最有前途(预测似然度大)的粒子扩展k时刻[137]，从而生成采样粒子。与SIR滤波器相比，当粒子的似然函数位于先验分布的引入辅助变量mk1时刻的粒子列表，应用贝叶斯定理，联合概率密度函数p(y

|xm)p(x|xm)wmk k k 1 1生成{x(im(i)}N

的重要性概率密度函数q(x,m|x

(2.29)k i1

k 0:k1 1:kq(x,m|x ,y)p(y|m)p(x|xm)wm

(2.30)k 0:k1 1:k k k k k1 k1mp(x

或预测均k kmE{x|xm。

k k k k k1qx|m,y)p(x|xm

，由于k 1:k k k1

(2.31)则有q(m|y)p(y|m)wm

(2.32)1:k k k 1kw(i)为kw(i)wm(i)

p(y|x(i)

)p(xi|xm(i))k k

p(y|x(i))k k k q(x,m|xk

|m(i))k k k

k k k(2.33)采用局部线性化的方法来逼近pxk|xk1,yk)是另一种提高粒子采样效率的有效方法。扩展Kalman粒子滤波与UncentedUKF或者EKF，并将它作为重要性概率密度函数[138]。另外，利用似然函数的梯度信息，采用牛顿迭代[139]或均值漂移等方法移动粒子至高似然区域，也是一种可行的方案，如图2.5所示。以上这些方法的共同特点是将最新的观测数据融入到系统状态的转移过程中，引导粒子到高似然区密度函数所需的粒子数。,N1}ki),}kk ̃(i),(i)k k i),̃i)k kk图2.5 结合均值漂移的粒子滤波算法Fig.2.5 Particlefiltercombinedwithmeanshift重采样方法Bootstrap的p(̃(i)x(i))̃(i){̃(i),̃(i)N更新为{x(i),1/NNk k k k k 1 k 1式(Multinomialresampling)重采样、残差重采样(Residualresampling)、分层重采样(Stratifiedresampling)与系统重采样(Systematicresampling)等。残余重采样法具有效率高、实现方便的

NiÑ(i)，其中

k k k k̃(i)N1Ñ(i)Ni选择余下的Nk k k的主要过程为

NNi2.6所示。残余重采样Ni1Nj,j1,⋯Nk。生成N[，1]lˆk；k 对于每个i，寻找归一化权值累计量大于或等于i的最小标号mm1

l

i落在区间[

xm2.6所示。mm1kkx(i经重采样后的个数为步骤(3)Nimm1kk图2.6 残差重采样Fig.2.6 ResidualResampling重采样并没有从根本上解决权值退化问题。重采样后的粒子之间不再是统计独立关系，给估计结果带来额外的方差。重采样破坏了序贯重要性采样算法的并行性，不利于VLSI硬目之间进行有效折衷。k1N个权值为1/Nx(ikk̃(i)kk̃(i)重新设置为1/N。k图2.7 SIR算法示意图Fig.2.7 SIRalgorithm标准的粒子滤波算法流程为：k0：对于i12,⋯Np(x)生成采样粒子{x(i)}N0k1,2⋯，循环执行以下步骤：

0 k ①重要性采样：对于i,⋯N{̃(i)k k̃(i)k②{̃(i)̃(i)}{x(i)1/N}；k k kk

Ñ̃̃ k k k退化现象，但重采样方法也会带来一些其它的问题。重采样需要综合所有的粒子才能实现，法是增加马尔可夫蒙特卡洛(Markovchainmontecarlo,MCMC)移动步骤[141]卡洛方法(GibbsMetropolis-Hastings采样等)利用不可约马尔可夫过程可逆平稳分布p̃k|kkk|̃k的马尔科夫链变换之后，若满足

k(k|̃k)p(̃k|k)d̃kp(k|k)(2.34)更为有利的位置。qx中生成粒子的具体步骤为：从[0,1]vU[0,1；k k x*px*|x(ik k q(x*)p(x(i)|x*) * *k k ，否则，