CN101819682A - 基于马尔科夫链蒙特卡洛粒子滤波的目标跟踪方法 - Google Patents

Abstract

本发明提供的是一种基于马尔科夫链蒙特卡洛粒子滤波的目标跟踪方法：1、初始时刻，由初始分布中得到一组初始粒子，并设置其初始的均值和方差；2、重要性采样；3、权值更新；4、得到归一化的权值；5、再采样；6、引入MCMC移动步骤；7、状态更新。本发明通过MCMC移动步骤，将粒子推向先验分布和后验分布都较大的区域，改善粒子多样性，在一定程度上抑制样本贫化问题。样本贫化问题的解决使得算法再采样的效果得到保证，进而可以提高滤波的精度。MCMC移动步骤比较容易实现，因此可以与其它改进步骤结合运用对粒子滤波进行优化。MCMC步骤的加入增加了滤波方法的运算量，但另一方面降低了精确估计所需要的粒子数，提高了滤波的效率。

Description

基于马尔科夫链蒙特卡洛粒子滤波的目标跟踪方法

技术领域

本发明提供的是一种函数跟踪方法，具体涉及融合了马尔科夫链蒙特卡洛移动步骤的粒子滤波方法(PF-MCMC)。

背景技术

非线性滤波方法目前广泛应用于信号处理以及导航制导、目标跟踪、金融分析、人工智能等相关领域。最早出现的非线性滤波算法是扩展卡尔曼滤波(EKF)，EKF的核心思想是对随机系统的非线性模型作线性近似，而且其噪声是基于高斯假设的，因此对于强非线性非高斯模型的滤波效果不佳，在工程应用中，EKF仅对某些特定的模型有效，对于很多非线性系统不能保证收敛和滤波精度。无迹卡尔曼滤波(UKF)也是较为常用的一种非线性滤波算法，UKF无须将模型进行线性化，直接使用非线性模型，避免了局部线性化引入的误差，避免在强非线性系统中出现发散。但由于UKF也是基于高斯假设，所以在工程应用中的很多非高斯模型上很难取得理想的滤波效果。以非参数化的蒙特卡罗模拟方法为基础的粒子滤波(PF)可以有效的弥补EKF和UKF的不足，PF方法的核心是利用些随机样本(粒子)来表示系统随机变量的后验概率密度，能得到基于物理模型的近似最优数值解，而不是对近似模型进行最优滤波。

PF最常见的问题是粒子退化现象，即经过几次迭代，除一个粒子外，所有粒子都只具有微小的权值，这意味着大量的计算工作都被用来更新那些对后验概率密度的估计几乎没有影响的粒子上。解决该问题的一个关键技术是再采样策略，其基本思想是通过后验概率密度函数再采样N次，产生新的粒子集，由于再采样是独立同分布，粒子权值被重新设置为1/N，再采样带来的负面问题是样本贫化现象，即高权值的粒子被多次复制，采样结果中包含了很多重复点，不能有效的反应状态变量的概率分布，从而丧失了粒子的多样性，甚至会导致滤波发散。因此，一个好的再采样算法应该能够在增加粒子多样性和减少小权值粒子数目之间进行有效平衡。

马尔可夫链蒙特卡罗(MCMC)算法是一种利用马尔科夫链无记忆性的随机模拟方法。其基本思想是：构造一条非周期、不可约的马尔科夫链(X₀，X₁，X₂…)，并且它有不变分布π，当这一马尔科夫链迭代次数接近无穷大时，得到的X_n的分布可以无限接近π。利用这种方法我们可以找到一条以P为不变分布的马尔科夫链，进而可以得到近似服从P的随机样本。MCMC方法主要应用在下面几个方面：用于生成比较复杂的随机数，可以实现对高维分布的采样，也是重要性采样的一种方法；实现高维积分的数值计算；用被估参数的贝叶斯分布的采样来估计参数；求复杂样本空间上函数的极值(模拟退火)。

发明内容

本发明的目的在于提供一种能提高滤波精度，改善粒子多样性的基于马尔科夫链蒙特卡洛粒子滤波的目标跟踪方法。

本发明的目的是这样实现的：

主要包括如下步骤：

第一步，初始化粒子，并设置其初始的均值和方差。

第二步，使用UKF/SSUKF或混合建议分布更新粒子：

通过标准无迹粒子滤波算法(UKF)或基于超球面采样的无迹粒子滤波(SSUKF)进行重要性采样即通过UKF/SSUKF获得粒子集的均值

和方差P_k ⁱ，从UKF/SSUKF结果

中采样产生N个粒子。

混合建议分布即通过UKF/SSUKF得到的分布产生一部分粒子，剩余粒子由先验分布产生。

第三步，权值更新 $w_k^i=w_{k-1}^i\frac{p\left(y_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{\mathrm{\π}\left(x_k^i\right|x_{k-1}^i,y_{1:k})}$

上式中的重要性概率密度函数 $\mathrm{\π}\left(x_k^i\right|x_{k-1}^i,y_{1:k})=N({\stackrel{\‾}{x}}_k^i,P_k^i),$ 引入了最新的观测值，因此改善了滤波器的性能。

第四步，得到归一化的权值

第五步，再采样。

第六步，引入MCMC移动步骤，使用马特普利斯-黑斯廷(Metropolis-Hasting)算法的具体步骤。

第七步，状态更新。

至此，算法一个时间步的运算结束，进入下一时间步的运算，转入第二步运行。

本发明针对再采样步骤引起的样本贫化现象，在这一步骤之后引入MCMC移动步骤，使得粒子群有序移动，得到更好的分布，改善了多样性，能够在定程度上减轻样本的贫化现象。

MCMC移动步骤本来用于不确定参数的系统辨识，Metropolis-Hasting算法是一种典型的MCMC方法，在再采样步骤中引入MCMC移动步骤可以提高滤波精度和对粒子个数进行调整。

在k时刻，利用马尔可夫过渡核函数

来生成一组新的样本用以代替重采样后的样本集

此时核函数

满足下列恒定性条件

$\∫\mathrm{\κ}\left(x_{0:k}\right|{\tilde{x}}_{0:k})p\left({\tilde{x}}_{0:k}\right|y_{1:k})d{\tilde{x}}_{0:k}=p\left(x_{0:k}\right|y_{1:k})$

因此，新的粒子集

仍然服从p(x_0：k|y_1：k)的分布。

Metropolis-Hasting算法的基本步骤如下：

第一步，任意选取马尔科夫链的初始状态X₀＝x。

第二步，从重要性分布中产生备选择点

其中k为算法迭代步数。

第三步，由均匀分布U(0，1)中产生一个参考分布，即v～U(0，1)

第四步，判断是否接受移动

如果

则接受移动，令否则拒绝移动，丢弃保留再采样的粒子，即

从上述步骤可以得知，MCMC方法主要是把粒子推向先验密度和后验密度都较大的区域。可以动态调整系统的性能和控制粒子数目的增减。

在粒子滤波的步骤中，最为重要的是重要性采样和再采样步骤。在再采样步骤之中加入MCMC算法。通过利用某一门限值判断是否接受移动，如果移动则将粒子推向后验密度较大的区域，如果不移动，则说明再采样的粒子得到了应有的效果无需移动。因此，这一步骤从一定程度上保证了再采样步骤的效果，相应也保证了算法的滤波精度。但是，加入了MCMC步骤之后，也相应加重了算法的计算负担，使得滤波的快速性有所下降。因此在需要保证快速性的场合，可以考虑使用相对简单的重要性采样方式或者对粒子数进行一定的控制，以减小算法整体的计算量。

本发明具有以下优点：

第一，通过MCMC步骤，粒子推向在先验密度和后验密度都较大的区域，从而在一定程度上解决了样本贫化问题。

第二，保证了再采样步骤的效果，提高了滤波算法的精度。

第三，在传统的SIR中或者改进后的EKPF、UPF等粒子滤波算法中都可以加入MCMC算法，应用范围较广，均可以取得比较理想的效果。

本发明针对粒子滤波的粒子退化以及再采样引起的样本贫化问题，在粒子滤波中引入统计了统计理论中的马尔科夫链-蒙特卡洛方法，在每个粒子上增加一个稳定分布作为后验概率密度的马尔科夫链蒙特卡洛移动步骤，由此提高滤波精度，改善粒子多样性。

附图说明

图1为PF-MCMC滤波方法流程图；

图2为UPF与UPF-MCMC在一次独立实验中的函数跟踪(N＝200)；

图3为SSUPF与SSUPF-MCMC在一次独立实验中的函数跟踪(N＝200)；

图4为UPF与混合分布UPF-MCMC在一次独立实验中的函数跟踪(N＝200)。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

本发明PF-MCMC滤波方法的基本步骤如下：

第一步，初始时刻，由初始分布p(x₀)中得到一组初始粒子，并设置其初始的均值和方差。

第二步，更新粒子：

使用UKF/SSUKF或混合建议分布更新粒子：

通过UKF或SSUKF进行重要性采样即通过UKF/SSUKF获得粒子集的均值和方差P_k ⁱ，从UKF/SSUKF结果

中采样产生N个粒子。

混合建议分布即通过UKF/SSUKF得到的分布产生一部分粒子，剩余粒子由先验分布产生，步骤如下：

(1)对每一个采样点x_k-1 ⁱ，应用UKF/SSUKF得到粒子集的均值

和方差P_k ⁱ。

(2)引入一个比例参数a(0＜a＜1)，从UKF结果

中产生a×N个粒子，从先验概率分布p(x_k|x_k-1 ⁱ)中产生(1-a)×N个粒子。

第三步，权值更新：

$w_k^i=w_{k-1}^i\frac{p\left(y_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{\mathrm{\π}\left(x_k^i\right|x_{k-1}^i,y_{1:k})}$

上式中的重要性概率密度函数引入了最新的观测值，因此改善了滤波器的性能。

第四步，得到归一化的权值

第五步，再采样：

定义N_eff来衡量有效粒子数量。如N_eff＜N_thr，即对

重新采样，产生新的集合{x_k ⁱ，i＝1，...，N}，重新设定粒子的权值为：

第六步，引入MCMC移动步骤：

即运用Metropolis-Hasting算法的具体步骤

第七步，状态更新：

至此，算法一个时间步的运算结束，进入下一时间步的运算，转入第二步运行。

实施实例一

单变量非稳定增长模型(UNGM)是在经济领域应用较为广泛的一种非线性模型，利用这一模型对PF-MCMC滤波方法的跟踪性能和一般粒子滤波性能进行比较，该模型是研究比较各种粒子滤波的算法性能的标准验证程序之一，其状态方程和观测方程如下：

$x_t=\frac{x_{k-1}}{2}+\frac{{25x}_{k-1}}{1+x_{k-1}^2}+8\cos \left(1.2t\right)+u_t$

$y_t=\frac{x_t^2}{20}+v_t$

其中观测噪声v_t～N(0，1)，系统噪声为如下的高斯和分布：

p(u_t)＝0.8N(u；0，1)+0.2N(u；0，10)

给定系统状态条件下的测量似然函数在y_t＞0时是双模态的，在y_t＜0是单模态的。由于状态初始值和滤波迭代步数与粒子滤波算法无关，因此可以不考虑这两个参数。另外，状态初始方差、测量噪声方差、过程噪声方差与重采样算法关系不大，而且在实际应用中，这三个参数都是由传感器等测量仪给出的。仿真实验中的滤波参数如下：状态初值为0.1，状态初始方差2，采样粒子数100，迭代时间步数100，再采样阈值0.9。混合建议分布MCMC滤波选用SSUKF和先验分布共同产生建议分布，比例参数a＝0.4。

仿真数据在表1中显示，实验表明，MCMC方法把粒子推向高似然函数区域，在过程噪声比较小的情况下，UPF-MCMC和SSUPF-MCMC和混合建议分布MCMC算法完全没有采样枯竭问题。有效样本个数可以衡量粒子滤波算法的退化程度，对于SIR滤波N_eff＝25.947；而UPF-MCMC滤波达到了N_eff＝72.214；SSUPF-MCMC滤波N_eff＝73.526；混合建议分布MCMC滤波N_eff＝73.029。

表1

	UPF	UPF-MCMC	SSUPF	SSUPF-MCMC	混合建议分布MCMC
						RMSE均值	4.1206	3.2093	4.1238	3.2253	3.3762
RMSE方差	1.1107	1.0961	1.1149	1.0980	1.0913

实施实例二

利用下面的非线性模型对滤波跟踪性能进行验证，其状态方程和观测方程如下：

$\left\{\begin{array}{c}{x}_k=1+\sin \left((4e-2)\mathrm{\π}(k-1)\right)+{0.5x}_{k-1}+v_{k-1}\\ y_k=\left\{\begin{array}{cc}{0.2x}_k^2+n_k& k\≤30\\ {0.5x}_k-2+n_k& k>30\end{array}\right.\end{array}\right.$

其中过程噪声v_t服从Gamma(3，2)分布，量测噪声n_k服从高斯分布N(0，10^-5)。系统状态的估计为

一次独立实验的均方误差为：

在实验中，分别采用UPF以及SSUPF与本文提出的UPF-MCMC、SSUPF-MCMC、混合建议分布MCMC方法进行比较。采用的粒子数N为200个，量测时间为T＝60，进行100次独立的实验，UPF的UT变换参数为α＝1，β＝0，κ＝2。混合建议分布MCMC滤波选用SSUKF和先验分布共同产生建议分布，比例参数a＝0.4。

图2给出了给出了UPF、UPF-MCMC在一次独立实验中的状态估计结果。。图3给出了SSUPF、SSUPF-MCMC在一次独立实验中的状态估计误差。图4给出了SSUPF、混合建议分布MCMC在一次独立实验中的状态估计误差。从以上图中可以看出，对于非线性、非高斯问题，PF算法得到的状态估计在部分时刻会严重偏离真实值，UPF与SSUPF算法得到的状态估计可以较好的吻合真实状态，在加入了MCMC步骤的三种算法的跟踪性能进一步提高。从图2到图4可以看出，在非线性较强的前30个时间步，加入了MCMC的三种算法的跟踪性能明显好于其余三种，而在后30时间步，各种算法的跟踪性能很相近。

表2

	UPF	UPF-MCMC	SSUPF	SSUPF-MCMC	混合建议分布MCMC
						RMSE均值	0.05201	0.05032	0.05227	0.05068	0.05104

	UPF	UPF-MCMC	SSUPF	SSUPF-MCMC	混合建议分布MCMC
						RMSE方差	0.003706	0.001695	0.003683	0.001702	0.001828
运算时间(S)	2.7291	6.6283	1.9883	4.3977	3.2165

通过上述模型将加入了MCMC步骤的粒子滤波跟踪方法与UPF和SSUPF三种方法的跟踪性能进行对比。主要对比数据为经过50次运行6种粒子滤波算法的均方根误差(RMSE)、运算时间，如表2所示。通过对这些性能指标的比较可以得到，在加入了MCMC步骤的三种算法中，UPF-MCMC、SSUPF-MCMC和混合建议分布MCMC算法的跟踪性能比较好，算法精度较高。但所用的计算时间大大超过前两种算法。由于在sigma采样时节省了一定的计算成本，SSUPF-MCMC的计算效率好于UPF-MCMC。由于采用了计算效率较高的重要性采样方式，所以混合建议分布MCMC算法的计算效率进一步有所改善。

Claims (1)

1.一种基于马尔可夫链蒙特卡罗粒子滤波的目标跟踪方法，其特征是主要包括如下步骤：

第一步，初始时刻，由初始分布p(x₀)中得到一组初始粒子，并设置其初始的均值和方差；

第二步，重要性采样

(1)根据标准无迹粒子滤波算法以及基于超球面采样的无迹粒子滤波对粒子的状态进行更新；或采取混合建议分布，即一部分粒子通过标准无迹粒子滤波或算法以及基于超球面采样的无迹粒子滤得到，剩余部分由先验分布产生；

(2)求粒子集的均值

和方差P_k ⁱ；

(3)从重要性密度函数中抽取粒子；

第三步，权值更新 $w_k^i=w_{k-1}^i\frac{p\left(y_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{\mathrm{\π}\left(x_k^i\right|x_{k-1}^i,y_{1:k})}$

上式中，w_k ⁱ为粒子对应的权值，p(·|·)为概率密度函数，重要性概率密度函数

为由均值和方差得到的建议分布；

第四步，得到归一化的权值 ${\tilde{w}}_k^i=\frac{{\stackrel{i}{w}}_k}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^j};$

第五步，再采样：定义阈值N_eff来衡量有效粒子数量，N_thr为有效粒子数，如N_eff＜N_thr，即对粒子集

重新采样，产生新的集合{x_k ⁱ，i＝1，...，N}，重新设定粒子的权值为：w_k ⁱ＝1/N，N为粒子数；

第六步，引入马尔科夫链蒙特卡洛移动步骤；

第七步，状态更新： $x\left(k\right|k)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^i{w}_k^i.$

Images