基于核密度的动态初始化重置粒子滤波
开篇：润墨网以专业的文秘视角，为您筛选了一篇基于核密度的动态初始化重置粒子滤波范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！
文章编号:1001-9081(2012)01-0295-04 doi:10.3724/SP.J.1087.2012.00295
摘 要:针对粒子滤波过程中，长时间的重采样造成的粒子多样性枯竭，由此导致目标跟踪中出现的精度下降及跟踪轨迹大幅振荡的现象，通过对采样粒子分布规律的研究，根据粒子枯竭的程度设置重置门限，在滤波过程中实时地检测粒子枯竭参数，当粒子的枯竭超过设置门限时，采用重置初始化粒子的方法来缓解采样粒子的枯竭趋势，有效地增加了长时间大量重采样后粒子的多样性，避免了粒子所含信息过多的丢失，显著地提高了粒子滤波的精度，在二维目标跟踪模型中应用所提算法并进行仿真实验，仿真结果证明了算法的可行有效。
关键词:
核密度；粒子滤波；采样粒子；机动目标跟踪；重采样
中图分类号: TP391.413; TP183 文献标志码:A
Abstract: It has been found that the accuracy of particle filtering is much lower when the maneuvering target tracking process has been executed for a long time. The reason for this problem is that the diversity of the sampled particles is rapidly lost because of the excessive resampling. Therefore, the track of the maneuvering target estimated by the particle filtering will be widely wiggly from the true one. Through the research of the distribution of the sampled particles, a new algorithm was proposed. And a detected threshold was set to detect if the particle was dried up badly. When the particle was dried up badly, the particles of the state-space would be reset to relax the degree, so the new particles could contain more distribution information. The new algorithm has a high capability in the simulation of the 2-D maneuvering target tracking.
Key words: kernel density; particle filtering; sampled particle; maneuvering target tracking; resampling
0 引言
在目标跟踪中，应用粒子滤波算法进行滤波跟踪时，当跟踪时间不是很长时，跟踪精度较高，但当跟踪时间较长时，出现跟踪轨迹在真实轨迹附近大幅振荡的现象，严重地影响了跟踪精度，给目标跟踪带来了不良的影响。
通过对粒子滤波算法的研究发现，引起这种现象的原因是由于粒子滤波长时间大量重采样导致采样粒子的多样性枯竭，致使粒子所包含的信息大量丢失，进而造成滤波精度大幅下降。
针对上述目标跟踪中所出现的现象，通过对采样粒子分布规律的研究，提出新型的动态初始化重置算法，设置粒子枯竭程度的重置门限，在滤波过程中实时地检测粒子枯竭统计量，当粒子的枯竭超过设置门限时，采用重置初始粒子的方式来解决粒子多样性枯竭问题，提高粒子滤波的跟踪精度。
1 粒子滤波
1.1 粒子滤波原理
粒子滤波算法［1-3］以贝叶斯滤波为理论基础，运用蒙特卡洛方法，通过重要密度函数进行粒子采样来近似后验概率密度p(xk|zk)，在最优估计理论的指导下，通过递推加权的方法来代替积分运算，并可在线计算状态向量的最优估计值。由于其非参数化的性质，使其从根本上摆脱了非线性滤波过程中系统、观测噪声必须为高斯分布的限制条件，因此其应用范围更加广泛。
在k时刻的状态的后验概率密度分布可离散地加权为：
p(xk|z1:k)≈∑Nsi=1w(i)kδx(i)k(xk－x(i)k)（1）
其中的权值通过重要采样法实现，粒子集{x(i)0:k,w(i)k}Nsi=1可由重要密度函数采样获得。
若粒子集{x(i)0:k}Nsi=1可由重要密度函数q(x0:k|z1:k)得到，则权值为：
w(i)k∝p(x(i)0:k|z1:k)q(x(i)0:k|z1:k)(2)
若重要密度能分解为：
q(x0:k|z1:k)=q(xk|x0:k－1,z1:k)q(x0:k－1|z1:k－1)（3）
则通过由q(xk|x0:k－1,z1:k)得到的粒子{x(i)k}Nsi=1与由q(x0:k－1|z1:k－1)得到的粒子集{x(i)0:k－1}Nsi=1可得到新的粒子集{x(i)0:k}Nsi=1。
由于后验概率密度函数［4］可表示为：
p(x0:k|z1:k－1)=p(x0:k|zk,z0:k－1)=p(zk|x0:k,z1:k－1)p(x0:k|z1:k－1)p(zk|z1:k－1)=p(zk|x0:k,z1:k－1)p(xk|x0:k－1,z1:k－1)p(zk|z1:k－1)p(x0:k－1|z1:k－1)=
p(zk|xk)p(xk|xk－1)p(zk|z1:k－1)p(x0:k－1|z1:k－1)∝p(zk|xk)p(xk|xk－1)p(x0:k－1|z1:k－1)（4）
将式（3）、（4）代入式（2），即可得到重要性权值更新公式(5)：
w(i)k∝p(zk|x(i)k)p(x(i)k|x(i)k－1)p(x(i)0:k－1|z1:k－1)q(x(i)k|x(i)0:k－1,z1:k)q(x(i)0:k－1|z1:k－1)=w(i)k－1p(zk|x(i)k)p(x(i)k|x(i)k－1))q(x(i)k|x(i)0:k－1,z1:k)（5）
如果q(xk|x0:k－1,z1:k)=q(xk|xk－1,zk)，则重要密度函数仅依赖于xk－1和zk，在计算时，仅需要存储{x(i)k}Nsi=1，而不必关心粒子集{x(i)k－1}Nsi=1和过去测量值z1:k－1。修正后的权值为：
w(i)k∝w(i)k－1p(zk|x(i)k)p(x(i)k|x(i)k－1)q(x(i)k|x(i)k－1,zk)（6）
1.2 重采样
在上述粒子滤波过程中，经过几次迭代,可能只有一个粒子有非零权值，其他粒子只具有微小权值的退化现象(degeneracy)［3,5-7］。退化现象意味着大量的计算工作都被用来更新那些对p(xk|z1:k)的估计几乎没有影响的粒子上。对退化现象的一个恰当的测度是有效采样尺度Neff,有效采样尺度定义为：
Neff=Ns1+Var(w(i)k)
一般情况下无法确切计算Neff的数值，但可以通过式(7)近似估计：
N^eff=1/∑Nsi=1(w(i)k)2（7）
当Neff≤Ns时，Neff越小，意味着退化现象越严重。一般情况下，通过在序贯重要采样(Sequential Importance Sampling, SIS)请补充SIS的中文名称和英文全称。之后进行重要性样本的重采
最常用的重采样方法是随机采样方法，其过程是：首先产生Ns个在［0,1］均匀分布的随机数{ui:i=1,2,…,Ns}，然后通过搜索算法找到满足以下条件的整数m使得∑m－1j=0w(j)k
记录样本x(m)k并将其作为新样本集中的采样，将区间［0,1］以λi=∑Nsj=0w(j)k(i=1,2,…,Ns)分成Ns个小区间，当随机数ui落在第m个区间Im=(λm－1,λ］“=”是否应该用“∈”来替换，因为是区间的取值，应该是用“∈”吧？时，对应样本x(m)k进行复制。
设粒子和其权值组成粒子组{x(i)0:k,w(i)k}Nsi=1，重采样的目的就是去掉权重较小的粒子，将较大权重对应的粒子组{x(i)0:k,w(i)k}Nsi=1复制成若干个权值为1Ns的粒子组x(i)0:k,1NsNsi=1，但过多地复制大权值的粒子会造成样本枯竭但过多的复制大权值的粒子的负作用“负作用”是专用术语，还是应该改为“副作用”？另外，此句不通顺，请核实并调整。就是造成了样本枯竭，较大权值的粒子被多次选取复制，采样结果中包含了多个重复的粒子点，从而损失了粒子的多样性。如图1所示。
标准粒子滤波算法［3-4,8］选择最易于实现的先验转移概率密度作为重要密度函数，即：
q(x(i)k|x(i)k－1,zk)=p(x(i)k|x(i)k－1)（8）
将式（8）代入式（6），可将重要性权值化简为：
w(i)k∝w(i)k－1p(zk|x(i)k)（9）
将权值w(i)k归一化，即w(i)k=w(i)k/∑Nsi=1w(i)k。
而后验概率密度p(xk|z1:k)可表示为：
p(xk|z1:k)≈∑Nsi=1w(i)kδ(xk－x(i)k)
2 跟踪中的振荡现象及解决方法
2.1 振荡现象
在标准粒子滤波过程中，每个迭代点都要进行重采样处理并且把重采样后的粒子权值设定为等权值的1/N，所以粒子多样性快速枯竭，当滤波时间过长时，导致采样粒子多样性枯竭严重，滤波精度下降，在二维目标跟踪中则表现为跟踪轨迹振荡现象。
由图2可以看出：在滤波前期，跟踪效果较好，当目标机动后，粒子滤波通过小幅振荡来适应目标机动，快速收敛，但由于过量的重采样，导致采样粒子多样性枯竭信息大量丢失，造成后面的跟踪大幅度地振荡，偏离目标轨迹，跟踪精度急剧下降。
产生跟踪轨迹大幅振荡的深层次根本原因是，采样粒子的性能下降，造成跟踪精度下降，长时间的重采样导致了采样粒子多样性的枯竭，性能下降。
图3是粒子滤波执行100步的过程中每次采样粒子的分布密度图，颜色越深的地方代表粒子的密度越高。从图中可以清楚地看出：前期粒子的分布比较广泛；而滤波后期粒子的分布逐渐集中在几个中心点周围，采样粒子的多样性枯竭，导致滤波效率下降。
由图4可以看出当粒子滤波执行到第25步时，采样粒子的范围在［-30,10］区间内，当执行到第50步时为［0,8］，在第75步时为［10,18］，在第100步时为［-10,20］，可以看出采样粒子范围在不断减小的情况下越来越集中，最后在第100步时，粒子主要集中在0和10这两个点附近，形成两个尖锐的波峰，导致采样粒子丢失过多的采样信息，造成滤波性能不断下降。
2.2 动态初始化重置算法
针对上述滤波过程中存在的问题，通过分析采样粒子的分布规律，提出了动态初始化重置粒子滤波算法，设置粒子枯竭门限，当采样粒子枯竭程度超过门限时，在当前估计状态的基础上对系统状态向量进行初始化重置来增加采样粒子的多样性。
通过大量的仿真实验，发现基于“核密度(kernel density)”的粒子统计量可以较好地描述粒子的分布情况，定量地给出粒子的枯竭程度。
其中“核(kernel)”用于称呼任意的光滑函数K，它满足K(x)≥0并且∫K(x)dx=1, ∫xK(x)dx=0, σ2K∫x2K(x)dx>1。采样粒子x(i)k的分布可以近似地用核函数描述：
p(x(i)k)≈∑Nsi=1ω(i)kKh(xk－x(i)k)（10）
其中：Kh(x)=1hnxK(x/h)是重新调整过的核密度K(•)，h>0是核的带宽，nx是采样粒子的x(i)k的维数，ω(i)k是归一化权值。
选择适当的核函数K(•)和带宽h使得采样粒子的统计积分均方误差均值最小，达到最佳的统计效果：
MISE()=E∫［(x(i)k)－p(x(i)k)］2 dx(i)k
其中(x(i)k)为采样粒子概率分布p(x(i)k)的近似。选择最优Epanechnikov核：
Kopt=nx+22cnx(1－x2), x
其中cnx是Rnx的单位球的体积。
当密度函数是单位协方差矩阵的高斯分布时，带宽的最佳选择为：
hopt=AN1/(nx+4)s（11）
其中A=［8c－1nx(nx+4)(2π)nx］1/(nx+4)。
式(11)成立的条件为采样粒子的分布满足高斯分布，分析发现一般在粒子滤波的后期采样粒子的分布出现双峰分布，但该最佳带宽对描述双峰的采样粒子分布仍有指导作用。
设采样粒子序列{x(i)k,w(i)k}Nsi=1，根据核密度的估计理论，确定带宽h，生成统计序列{m(i)k}Nmi=1,其此处语句不通顺，请作相应调整。为各个带宽内含有的粒子个数，其中Nm=round(max(x(j)k)－min(x(j)k)h)为采样粒子的分布范围按照带宽划分的区域个数。
m(i)k=f^Ns(x)•Ns=∑Nsj=11hK(x－x(j)kh)x=i•h+min(x(j)k)
实际上，m(i)k就是落入某一带宽内的粒子个数，即集合{x(j)k|min(x(i)k)+(j－1)•h≤x(j)k
计算上述序列{m(i)k}Nmi=1的统计量σ2k：
σ2k=∑Nmi=1(m(i)k－k)2（12）
其中k为序列{m(i)k}Nmi=1的均值：
k=1Nm•∑Nmi=1m(i)k（13）
上述方法适用于单纯的单峰函数分布，当粒子分布为双峰甚至更复杂时，只比较统计数量σ2k已不能反映粒子的好坏了。
从图4中此处指代不清晰，指代哪个图，请表明。可以看出，当分布为双峰时，大量的粒子集中在两个中心点周围，而大部分区域没有粒子分布，导致粒子分布的过于集中、多样性枯竭，因此在求解统计量σ2k时，应该先从统计序列{m(i)k}Nmi=1中剔除没有粒子存在的零值点，经实验分析求得的统计量2k很好地定量地反映了粒子分布状态的好坏，既适用于粒子分布的单峰也适用于双峰的复杂情况。
通过剔除粒子统计序列{m(i)k}Nmi=1中的零值点得到新的粒子统计序列{(i)k}mi=1，计算统计序列的统计量2k：
2k=∑Nmi=1((i)k－k)2（14）
其中k为粒子统计序列{(i)k}mi=1的均值：
k=1m•∑Nmi=1(i)k（15）
通过上述分析，为简化算法，在SIR算法的基础上，动态重置粒子算法归纳如下。
1)设置粒子枯竭门限σ2threshold。
2)标准粒子滤波算法。
3)在每次重采样后，根据式(14)计算粒子分布方差2k。
程序前
For i=1:Ns
根据当前状态向量，重新初始化粒子
END For
程序后
5)标准粒子滤波算法。
上述过程中，σ2threshold值的选择可以依据经验值，也可以取粒子滤波前端的几个粒子序列统计量2k的均值。
3 仿真验证
仿真实验中选取系统状态方程为：
xk=xk－12+25•xk－11+x2k－1+8•cos(1.2•(k－1))+vk－1yk=x2k20+wk
其中vk-1、wk为服从N(0,1)的高斯白噪声。应用本文新算法对上述系统方程进行滤波计算并与标准PF算法进行对比。
通过仿真实验，计算评判总体跟踪效果的指标，得均方根误差(Root Mean Square Error, RMSE)如下。
PF算法的RMSE为：
RMSE=1N∑Nk=1(xk－k)2=3.2163
动态初始化重置PF算法的RMSE为：
RMSE=1N∑Nk=1(xk－k)2=3.1589
通过上述仿真实验结果可知在一维情况下，该新型算法的跟踪效率较标准粒子滤波算法得到了有效地提高。
下面将该新型算法应用到二维机动目标跟踪中。
目标状态方程采用文献［8］中运动模型，其状态方程描述如下：
Xk+1=f(Xk)+BvVk
zk=arctan(ηk/ξk)+wk
其中：Xk=(ξk ξ •k ηk k)T，ξk和ηk目标在二维坐标系中的坐标位置，ξ •k和k分别为对应方向上的速度分量。系统方程中：
f(X)=AX=1T000100001T0001X
Bv=0.5T0T000.5T0T
Vk=(vx,vy)Tk为系统白噪声，满足零均值方差为diag(102,102)的高斯分布；wk为观测噪声，满足方差为(3/180×π)2的高斯分布；目标初始状态为X0=(400,40,4800,－30)T。
在上述系统状态方程的基础上，运用标准粒子滤波和新型动态初始化重置粒子滤波算法对机动目标轨迹进行滤波跟踪，粒子采样数为300，仿真时间步长为100,采样时间T为1s。
从图5中可以看出标准粒子滤波算法在跟踪的后期效率急剧下降，对目标的跟踪严重偏离目标真实轨迹，已不能有效跟踪目标轨迹。而动态初始化重置粒子滤波算法因为实时监控粒子的枯竭程度，所以保证了粒子的多样性，使采样粒子保留了更多的信息，从图5中可以清楚地看到，该新型算法对机动目标跟踪更准确，性能更好，尤其是在目标做大机动运动时，其跟踪性能明显好于标准粒子滤波，其振荡幅度较标准粒子滤波更小，收敛速度更快，精度更高。
计算评判总体跟踪效果的指标，均方根误差RMSE得：PF算法RMSE为RMSE1=233.851207；动态初始化重置PF算法RMSE为RMSE2=194.152953。
4 结语
在研究机动目标跟踪的过程中，发现随着跟踪时间的加长，粒子滤波的跟踪效果越来越差，并出现大幅度的振荡现象，通过分析采样粒子的分布规律发现，由于过度地重采样造成粒子多样性枯竭，导致粒子中所含信息的丢失，造成滤波性能下降。在核密度统计理论的基础上，提出了实时检测粒子枯竭程度，采取粒子适时初始化重置的方法来缓解粒子枯竭的趋势，通过二维机动目标跟踪的仿真实验可知，新型动态初始化重置粒子滤波算法的性能优越，对长时间的机动目标跟踪有着较好的适应能力。
参考文献:
[1]
JOHANSEN A M, DOUCET A. A note on auxiliary particle filters ［J］. Statistics and Probability Letters, 2008, 78(12): 1498-1504.
[2]
KWORK C, FOX D, MEILA M. Real-time particle filters ［J］. Proceedings of the IEEE, 2004, 92(3): 469-484.
[3]
DOUCET A, JOHANSEN A M. A tutorial on particle filtering and smoothing: Fifteen years later ［R］. Vancouver, BC: University of British Columbia, Department of Statistics, 2008: 12-14.
[4]
DOUCET A, BRIERS M, SENECAL S. Efficient block sampling strategies for sequential Monte Carlo methods ［J］. Journal of Computational and Graphical Statistics, 2006, 15(3): 693-711.
[5]
VO B-N, VO B-T, SINGH S. Sequential Monte Carlo methods for static parameter estimation in random set models ［C］// Proceedings of the 2004 Intelligent Sensors, Sensor Networks and Information Processing Conference. Washington, DC: IEEE Computer Society, 2004: 313-318.
[6]
CLARK D E, BELL J. Convergence results for the particle PHD filter ［J］. IEEE Transactions on Signal Processing, 2006, 54(7): 2652-2661.
[7]
冯驰,王萌,汲清波,等.粒子滤波器重采样算法的分析与比较［J］.系统仿真学报,2009,21(4):1101-1110.
[8]
KARLSSON R, BERGMAN N. Auxiliary particle filters for tracking a maneuvering target ［C］// Proceedings of the 39th IEEE Conference on Decision and Control. Washington, DC: IEEE Computer Society, 2000: 3891-3895.
[9]
FEARNHEAD P. Computational methods for complex stochastic systems: A review of some alternatives to MCMC ［J］. Statistics and Computing, 2008, 18(2): 151-171.
[10]
杜正聪,唐斌,李可.混合退火粒子滤波器［J］.物理学报,2006,55(3):999 -1004.
[11]
BRASNETT P, MIHAYLOVA L, BULL D, et al. Sequential Monte Carlo tracking by fusing multiple cues in video sequences ［J］. Image and Vision Computing, 2007, 25(8): 1217-1227.
[12]
KIM J S, SERPEDIN E, SHIN D R. Improved particle filtering-based estimation of the number of competing stations in IEEE 802.11 networks ［J］. IEEE Signal Processing Letters, 2008, 15(1): 87-90.
[13]
莫以为,萧德云.基于进化粒子滤波器的混合系统故障诊断［J］.控制与决策,2004,19(6):611-615.
[14]
胡士强,敬忠良.粒子滤波算法综述［J］.控制与决策,2005,20(4):361-365.
[15]
GIREMUS A, TOURNERET J Y, CALMETTES V. A particle filtering approach for joint detection/estimation of multipath effects on GPS measurements ［J］. IEEE Transactions on Signal Processing, 2007, 55(4): 1275-1285.
[16]
DUAN ZHUOHUA, CAI ZIXING, YU JINXIA. An adaptive particle filter for soft fault compensation of mobile robots ［J］. Science in China Series F: Information Sciences, 2008, 51(12): 2033-2046.
收稿日期:2011-06-27;修回日期:2011-09-05。
作者简介:
白剑锋（1985-），男，河北沧州人，硕士研究生，主要研究方向：目标跟踪、信息融合； 南建国（1969-），男，陕西西安人，副教授，博士，主要研究方向：目标跟踪、信息融合； 邬蒙（1979-），男，陕西西安人，讲师，硕士，主要研究方向：目标跟踪、信息融合； 查翔(1988-)，男，山东泰安人，硕士研究生，主要研究方向：目标跟踪、信息融合。