CN103152014A - Metropolis-Hastings变异粒子群重采样粒子滤波器实现方法 - Google Patents

Abstract

本发明涉及基于Metropolis-Hastings变异的粒子群重采样粒子滤波器的实现方法。为了解决粒子滤波在粒子数量较少时估计精度不高的问题，本发明提供了一种基于Metropolis-Hastings(MH)变异的粒子群重采样及粒子滤波器的实现方法。该方法将Metropolis-Hastings(MH)移动作为粒子群优化的变异算子，通过将MH变异规则与粒子群的速度-位置搜索过程相结合，使得重采样后的粒子群更接近真实的后验概率密度分布，有效解决了一般的变异粒子群算法容易发散的问题，加快了粒子滤波在序贯估计过程中的收敛速度，提高了其估计精度。仿真试验证明，基于MH变异的粒子群优化粒子滤波器可以有效地克服粒子贫化现象，改善对非线性系统的跟踪估计效果。

Description

Metropolis-Hastings变异粒子群重采样粒子滤波器实现方法

技术领域

本发明属于数字信号处理领域，更具体的涉及非线性滤波领域，提供一种粒子滤波重采样以及粒子滤波器的实现方法。

背景技术

粒子滤波器在非高斯、非线性系统状态估计方面的应用非常广泛，特别在导航制导、目标跟踪、金融分析、人工智能、盲信号处理等领域已表现出较好的跟踪估计性能。但粒子滤波器存在的一个不可避免的问题就是粒子退化（ParticleDegeneracy）现象，即粒子集经过若干次迭代之后，除少数粒子外，大部分粒子仅具有微小的权值（粒子权值小意味着对后验概率密度的贡献小），并且小权值的粒子还要参加后续的迭代计算，增加无用的计算量。

粒子重采样是解决粒子退化的有效手段，常规的粒子重采样方法（如残差重采样、分层重采样、系统重采样等）虽然能一定程度上解决退化问题，但却减少了粒子集合的多样性，即样本贫化（Particle Impoverishment）。样本贫化产生的原因是，每次重采样后将改善粒子的权值，使其不再为零，但由于重采样将高权值的粒子过多复制，有效粒子数目不断减少，经过若干次递推采样后，有效粒子将被重采样复制耗尽，直至仅剩最后一个权值为1的样本，此时，样本的分布已退化为单点分布，不能反映状态的真实后验分布。

为解决上述问题，研究者们提出将生物进化思想引入到粒子重采样中，并且取得了较好的研究结果。Park.S（2007）、叶龙（2007）、胡振涛（2009）等提出将遗传算法和粒子重采样相结合，虽然改善了重采样后粒子集合的多样性，但却急剧增加了计算量。粒子群优化算法（Particle Swarm Optimized,PSO）是一种较新的模拟生物进化算法，其模拟鸟群的觅食行为，通过种群中个体之间的集体协作使群体达到最优，沈艳（2005）的研究结果表明粒子群优化算法的计算效率优于遗传算法，因此，方正（2006）、佟国锋（2007）将标准粒子群优化规则引入到粒子重采样中，在减少计算量的前提下保证了滤波器的估计精度。然而，Kennedy.J（1995）的研究表明，标准粒子群优化算法存在早熟收敛现象，容易陷入局部最优解。针对这一问题，吕振肃（2004）、王海峰（2009）、陈建超（2009）提出了普通变异的粒子群优化改进，增强了粒子群跳出局部最优解的能力，但其在变异移动过程中未考虑状态的概率密度分布，进而增加了变异后发散的风险。相比普通变异方法而言，马尔可夫链蒙特卡罗(Markov chain Monte Carlo，MCMC)Metropolis-Hastings(MH)移动变异方法提供了一种从后验分布中采样和变异优质样本的机制，即根据MCMC方法建立马尔可夫链，使它的平稳分布与后验分布相同。当马尔可夫链收敛时，变异后重采样的粒子可以看作由后验分布中抽取的样本，而MH移动选择可以从建立的马尔可夫链中选择优质样本。

发明内容

本发明提供一种基于Metropolis-Hastings变异的粒子群重采样及粒子滤波器实现方法（Metropolis-Hastings Mutation Particle Swarm OptimizedParticle Filter,MHMPSOPF），其目的重在提供一种利用较少粒子数量保证粒子多样性并达到较高估计精度的粒子滤波器实现方法。

为达到上述目的，本发明采用了如下技术方案：

粒子滤波的理论基础来源于贝叶斯滤波原理，其实质是利用已知信息来构造系统状态度量的后验概率密度。即对于一个非线型离散系统：

$X_k=F_k(X_{k-1},V_k)$

$Y_k=H_k(X_k,U_k)$

其中，X_k与Y_k分别表示目标状态值与观测值，F_k(·)表示系统的非线性状态转移函数，V_k为系统的状态转移噪声；H_k(·)表示系统的非线性观测函数，U_k为观测噪声。粒子滤波的目的就是通过观测值Y_k估计状态X_k。根据贝叶斯滤波的思想，状态X_k的信息可以由后验概率密度p(X_k|Y_k)估计得到，粒子滤波则通过寻找一组在状态空间中传播的随机样本（粒子）对后验概率分布进行近似，当样本点数足够大时，便可获得状态的最小方差估计。为了获得有效粒子集合，粒子滤波采用序贯重要性采样(Sequential Importance Sampling)的方法，通过对粒子集合进行重要性递归修正，进而获得后验概率的近似分布。为了解决粒子序贯估计中的粒子贫化问题，人们提出了残差重采样、系统重采样等策略，但上述重采样策略都需要大量粒子逼近系统的状态值，当粒子数较少时将丧失粒子集的多样性，从而引起较大的估计误差。

[0008]作为一种有效的仿生优化方法，粒子群优化机制为解决粒子贫化问题提供了重要的指导思想。其实现原理如下：

第一：初始化粒子集合

其中，1≤i≤N，代表粒子集合中的第i个粒子，N为粒子总数量；并随机设定粒子的初始位置

和初始速度

其中，1≤k≤K，代表第k个采样点，K为信号的总采样点数；

第二：根据系统的状态转移函数F_k(·)，进行粒子的状态预测，即

第三：利用系统的观测方程计算粒子的似然分布值：

其中Y_k为观测值，Y_pred为各个粒子预测的观测值，R_k为观测噪声的方差，p(Y_k|X_k)为似然函数；

第四：更新粒子状态权值，即 ${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^{i}p\left(Y_k\right|X_k),$ 并根据公式 ${\hat{N}}_{\mathrm{eff}}=1/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left({\mathrm{\ω}}_k^i\right)}^2$ 判断有效粒子数；

第五：若有效粒子数低于阈值，则将似然函数计算公式作为粒子群优化的适应度函数，似然分布值作为粒子的适应值执行第六；否则返回执行第二；

第六：根据粒子群优化规则进行重采样，即：对于单个粒子，将其当前适应值与其所经历过的最佳位置P_t的适应值进行比较，若大于粒子最佳位置P_g的适应值，则将当前适应值作为最佳位置；对于全部粒子，将每个粒子所经历过最佳位置P_t的适应值与粒子群的全局最佳位置P_g的适应值进行比较，若大于全局最佳位置P_g的适应值，则将其作为群体当前的全局最佳位置；

[0015]根据以下公式更新每个粒子的速度和位置：

$v_{k+1}^i=\mathrm{\λ}*v_k^i+c_1*\mathrm{rand}*(p_t-x_k^i)+c_2*\mathrm{rand}*(p_g-x_k^i)$

$x_{k+1}^i=x_k^i+v_{k+1}^i$

其中，λ、c₁、c₂为模型参数，rand是正态分布的随机数。

由于标准的粒子群优化算法易出现收敛早熟问题，目前提出的变异粒子群优化方法，虽然增强了粒子群跳出局部最优解的能力，但在移动变异过程中未考虑粒子群的后验概率密度分布，进而增加了变异后发散的风险。为了克服普通变异方法的心声粒子不易收敛的缺点，提出一种基于MCMC-MH变异的粒子群优化方法，以依据后验概率密度分布产生有效的粒子，提高粒子群的多样性提高跟踪精度。

MCMC方法提供了一种从后验分布中抽取样本的机制，即根据MCMC方法建立马尔可夫链，使它的平稳分布与后验分布相同。当马尔可夫链收敛时，重采样的粒子可以看作由后验分布中抽取的样本。MH移动变异是MCMC的一种实现途径，其通过构造一个平稳分布为建议分布的马尔可夫链来获得新粒子集，其实现原理为：

根据给定的转移函数（本文中为粒子的建议分布）进行扰动，产生一个候选粒子x'，如

$q(x,x^{\′})=q\left(\right|x-x^{\′}\left|\right)\∝\exp \left[\frac{{(x-x^{\′})}^2}{{2\mathrm{\δ}}^2}\right],$ δ为正态分布的方差；

计算接受率

$q(x,x^{\′})=\min [\frac{\mathrm{\π}\left(x^{\′}x\right)q(x^{\′},x)}{\mathrm{\π}\left(x^{\′}x\right)q(x,x^{\′})},1]$

按照u～U(0,1)均匀采样；如果u≤q(x,x')，则x^k+1=x'；否则，x^k+1=x；

根据上述的MH移动变异规则，改进标准粒子群优化粒子滤波算法，即在粒子群的速度-位置更新之后，对最佳位置粒子进行MH移动变异。其实现原理为在前述的标准粒子群优化粒子滤波的第六之后增加第七；

第七：将当前最佳粒子作为变异公式中的x，依公式产生一个新的x'，并依据选择公式计算变异后粒子的接受率，如果u≤q(x,x')，则选择变异后粒子为最佳粒子；否则保持不变。

增加了MH变异步骤后，基于MH变异的粒子群优化粒子滤波的实现流程见附图1所示。

附图说明

图1为本发明的实施过程流程图；

图2为标准粒子滤波(Particle filter,SRPF)、粒子群优化粒子滤波(Particle swarm optimized particle filter,PSOPF)、变异粒子群优化粒子滤波(Mutation particle swarm optimized particle filter,MPSOPF)和本文算法(MHMPSOPF)对UNGM系统模型（一维变量）的跟踪估计结果；

图3为标准粒子滤波(Particle filter,SRPF)、粒子群优化粒子滤波(Particle swarm optimized particle filter,PSOPF)、变异粒子群优化粒子滤波(Mutation particle swarm optimized particle filter,MPSOPF)和本文算法(MHMPSOPF)对Bearings-only系统模型（多维变量）的跟踪估计结果；

具体实施方式

下面详细叙述本发明的具体实现过程，并结合非线性系统模型进行举例。

本发明的MH变异粒子群优化重采样粒子滤波的执行步骤如下：

步骤1：初始化粒子集合

其中，1≤i≤N，代表粒子集合中的第i个粒子，N为粒子总数量；并随机设定粒子的初始位置和初始速度其中，1≤k≤K，代表第k个采样点，K为信号的总采样点数；

步骤2：根据系统的状态转移函数F_k(·)，进行粒子的状态预测，即

步骤3：利用系统的观测方程计算粒子的似然分布值：其中Y_k为观测值，Y_pred为各个粒子预测的观测值，R_k为观测噪声的方差，p(Y_k|X_k)为似然函数；

步骤4：更新粒子状态权值，即 ${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^{i}p\left(Y_k\right|X_k),$ 并根据公式 ${\hat{N}}_{\mathrm{eff}}=1/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left({\mathrm{\ω}}_k^i\right)}^2$ 判断有效粒子数；

步骤5：若有效粒子数低于阈值，则将似然函数计算公式作为粒子群优化的适应度函数，似然分布值作为粒子的适应值执行步骤6；否则返回执行步骤2；

步骤6：根据粒子群优化规则进行重采样，即：

（1）对于单个粒子，将其当前适应值与其所经历过的最佳位置P_t的适应值进行比较，若大于粒子最佳位置P_g的适应值，则将当前适应值作为最佳位置；

（2）对于全部粒子，将每个粒子所经历过最佳位置P_t的适应值与粒子群的全局最佳位置P_g的适应值进行比较，若大于全局最佳位置P_g的适应值，则将其作为群体当前的全局最佳位置；

（3）根据以下公式更新每个粒子的速度和位置：

$v_{k+1}^i=\mathrm{\λ}*v_k^i+c_1*\mathrm{rand}*(p_t-x_k^i)+c_2*\mathrm{rand}*(p_g-x_k^i)$

$x_{k+1}^i=x_k^i+v_{k+1}^i$

其中，λ、c₁、c₂为模型参数，rand是正态分布的随机数。

步骤7：MH移动变异的实现流程为：

（1）根据给定的转移函数（本发明中为粒子的建议分布）进行扰动，产生一个候选x'，如

$q(x,x^{\′})=q\left(\right|x-x^{\′}\left|\right)\∝\exp \left[\frac{{(x-x^{\′})}^2}{{2\mathrm{\δ}}^2}\right],$ δ为正态分布的方差；

（2）计算接受率

$q(x,x^{\′})=\min [\frac{\mathrm{\π}\left(x^{\′}x\right)q(x^{\′},x)}{\mathrm{\π}\left(x^{\′}x\right)q(x,x^{\′})},1]$

（3）按照u～U(0,1)均匀采样；

（4）如果u≤q(x,x')，则x^k+1=x'；否则，x^k+1=x；

（5）将当前最佳粒子作为变异公式中的x，依公式产生一个新的x'，并计算变异后粒子的接受率，如果u≤q(x,x')，则选择变异后粒子为最佳粒子；否则保持不变。

利用单变量非静态增长模型(Univariate nonstationary growth model,UNGM)和纯角度跟踪模型(Bearings-only tracking model,BOTM)进行了性能仿真分析。由于本文是针对普通变异粒子群优化粒子滤波算法做的改进，并且基于系统重采样（System resampling,SR）的标准粒子滤波作为基础方法被广泛采用。因此，为了充分比较，本文对标准粒子滤波(Particle filter,SRPF)、粒子群优化粒子滤波(Particle swarm optimized particle filter,PSOPF)、变异粒子群优化粒子滤波(Mutation particle swarm optimized particle filter,MPSOPF)和本文算法(MHMPSOPF)进行了分析比较。以下主要比较几种方法得到的置信区间和RMSE。

置信区间为：

$({\hat{x}}_k-\frac{{\mathrm{\σ}}_k}{\sqrt{N}}{U}_{\mathrm{\α}/2},{\hat{x}}_k+\frac{{\mathrm{\σ}}_k}{\sqrt{N}}{U}_{\mathrm{\α}/2})$

其中，N为粒子数量，

为k时刻的状态估计值，σ_k为状态估计的方差，U_α/2为标准正态分布α/2上的侧分位数。

RMSE定义为：

$\mathrm{RMSE}={\left[\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{({\hat{x}}_k-x_k^i)}^2\right]}^{1/2}$

仿真分析中粒子群优化的参数设置为，λ=0.65，c₁=2.1，c₂=2.1。每个数据模型进行100次蒙特卡罗仿真实验。实验环境为Intel Core Duo2.66GHz，Matlab2008B。

式中，u_k～N(0,10)、v_k～N(0,1)为相互独立的高斯白噪声。模型参数h₁=0.5，h₂=25，h₃=8，K=50。置信区间设置为α=0.05。

图2所示的结果是在取10个粒子的单次跟踪结果的置信区间分析比较，4种方法对应的RMSE见表1。图2(2-1)中SRPF得到的置信区间最窄，是4种方法中的最优置信区间，但跟踪估计误差最大；图2(2-2)中粒子群优化重采样得到的置信区间大，有效地增加了粒子多样性，但跟踪估计精度提高较小；图2(2-4)中MH变异的粒子群优化重采样跟踪估计差最小，且置信区间较窄，表明增加了多样性，而且保持了粒子的精度。

表1UNGM模型四种算法的性能比较（100次实验平均）

从表1可以看出，MHMPSOPF算法仅10个粒子的滤波精度就优于SRPF算法30个粒子的滤波精度。在粒子数较少时，MHMPSOPF算法比SRPF算法在跟踪精度上提升很多，而在粒子数较多时性能上的提高就不明显，这是因为粒子数的增加缓解了粒子贫乏和退化问题。通过对算法的估计精度和运行时间进行比较可以看出：虽然粒子数量相同时MHMPSOPF耗时最多，但其用较少的粒子数量即可达到较高的估计精度，综合估计精度与粒子数量而言，提高了算法的总体运行效率。因此，在实际算法选择时，MHMPSOPF粒子滤波是一个实用有效的算法。

Bearings-only模型实验：

表2给出的是四种算法100次仿真实验的RMSE均值。图3所示的结果是取50个粒子的单次跟踪效果及置信区间分析比较，4种方法对应的RMSE见表2。可以看出，表2与图3所示的跟踪效果及置信区间分布与表1和图2所示结果相同，可以看出，本文算法也能较好地解决高维模型的状态跟踪和估计问题。

通过上述实验可以发现，本文在MH变异粒子群优化规则的基础上，利用MH变异算子来改进传统的全变异粒子群优化规则，进而提出了一种基于MH变异粒子群优化的粒子滤波重采样方法，本发明由于按照状态的后验概率密度产生新粒子，保证了重采样粒子的有效性，提高了滤波器的跟踪估计精度。仿真结果表明，在粒子数量较少的情况下，MHMPSOPF算法的滤波性能明显优于PF、PSOPF、MPSOPF算法，且少量增加了运算时间，为减少粒子滤波所需的粒子数量和提高粒子滤波的跟踪估计精度提供了一条有效途径。

表2Bearings-only模型四种算法的性能比较（100次实验平均）

Claims (4)

1.一种基于Metropolis-Hastings变异的粒子群重采样粒子滤波器的实现方法，其特征主要包括如下步骤：

步骤1，初始化粒子集合

其中，1≤i≤N，代表粒子集合中的第i个粒子，N为粒子总数量；并随机设定粒子的初始位置和初始速度

其中，1≤k≤K，代表第k个采样点，K为信号的总采样点数；

步骤2，根据系统的状态转移函数F_k(·)，进行粒子的状态预测，即

步骤3，利用系统的观测方程计算粒子的似然分布值,

其中X_k为目标状态值，Y_k为观测值，Y_pred为各个粒子预测的观测值，R_k为观测噪声的方差，p(Y_k|X_k)为似然函数；

步骤4，更新粒子状态权值，即 ${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^{i}p\left(Y_k\right|X_k),$ 并根据公式 ${\hat{N}}_{\mathrm{eff}}=1/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left({\mathrm{\ω}}_k^i\right)}^2$ 判断有效粒子数；

步骤5，若有效粒子数低于阈值，则将似然函数计算公式作为粒子群优化的适应度函数，似然分布值作为粒子的适应值执行步骤6；否则返回执行步骤2。

步骤6，根据粒子群优化规则进行重采样，更新每个粒子的速度和位置；

步骤7，在粒子群的速度、位置更新之后，对最佳位置粒子进行MH移动变异。

2.如权利要求1所述的基于Metropolis-Hastings变异的粒子群重采样粒子滤波器的实现方法，其特征在于，所述步骤6根据粒子群优化规则进行重采样，更新每个粒子的速度和位置具体为：

1)对于单个粒子，将其当前适应值与其所经历过的最佳位置P_t的适应值进行比较，若大于粒子最佳位置P_g的适应值，则将当前适应值作为最佳位置；

2)对于全部粒子，将每个粒子所经历过最佳位置P_t的适应值与粒子群的全局最佳位置P_g的适应值进行比较，若大于全局最佳位置P_g的适应值，则将_Pt作为群体当前的全局最佳位置；

3)根据以下公式更新每个粒子的速度和位置，

$v_{k+1}^i=\mathrm{\λ}*v_k^i+c_1*\mathrm{rand}*(p_t-x_k^i)+c_2*\mathrm{rand}*(p_g-x_k^i),$

$x_{k+1}^i=x_k^i+v_{k+1}^i$

其中，λ、c₁、c₂为模型参数，rand是正态分布的随机数。

3.如权利要求2所述的基于Metropolis-Hastings变异的粒子群重采样粒子滤波器的实现方法，其特征在于，所述步骤7在粒子群的速度、位置更新之后，对最佳位置粒子进行MH移动变异具体为：

4)根据给定的转移函数进行扰动产生一个候选粒子x'如

$q(x,x^{\′})=q\left(\right|x-x^{\′}\left|\right)\∝\exp \left[\frac{{(x-x^{\′})}^2}{{2\mathrm{\δ}}^2}\right],$

其中δ为正态分布的方差；

5)计算接受率，

$q(x,x^{\′})=\min [\frac{\mathrm{\π}\left(x^{\′}x\right)q(x^{\′},x)}{\mathrm{\π}\left(x^{\′}x\right)q(x,x^{\′})},1];$

6)按照u～U(0,1)均匀采样；

7)如果u≤q(x,x')，则x^k+1=x'；否则，x^k+1=x。

8)将当前最佳粒子作为变异公式中的x，依公式产生一个新的x'，并计算变异后粒子的接受率，如果u≤q(x,x')，则选择变异后粒子为最佳粒子；否则保持不变。

4.如权利要求3所述的基于Metropolis-Hastings变异的粒子群重采样粒子滤波器的实现方法，其特征在于，所述给定的转移函数为粒子的建议分布。

Images