CN103684352A - 基于差分演化的粒子滤波方法 - Google Patents

Abstract

本发明公开了基于差分演化的粒子滤波方法，属于信号处理、非线性滤波领域。该发明是对粒子滤波的改进，主要创新集中在：一、为了充分利用最新的观测信息，采用无迹卡尔曼滤波来产生粒子滤波的重要性分布；二、对重要性分布产生的采样粒子不再做传统重采样操作，而是对采样粒子做差分迭代优化处理，从而得到最优的粒子点集。本发明一改传统重采样算法对粒子的简单复制与丢弃的做法，而是采用差分演化算法对粒子集进行进化寻优，有效缓解了粒子滤波中的粒子退化和贫化问题，提高了粒子的利用率，具有更高的估计精度，在信号处理、非线性滤波领域具有较好的应用前景。

Description

基于差分演化的粒子滤波方法

技术领域

本发明属于信号处理、非线性滤波领域，涉及粒子滤波技术，具体说是基于差分演化的粒子滤波方法。 

背景技术

非线性系统状态估计问题广泛存在于目标跟踪、信号处理、自动控制、人工智能、无线通信和金融分析等诸多领域。 

扩展卡尔曼滤波器（Extented Kalman Filter,EKF）是最经典的非线性滤波方法之一，其基本思想是通过对非线性系统进行线性化操作，实现对非线性系统状态的估计。由于其线性化误差较大，且许多实际问题中很难得到非线性函数的Jacobian矩阵，使其估计精度和使用范围受到了很大限制。之后，出现了无需计算Jacobian矩阵的无迹卡尔曼滤波器（Unscented Kalman Filter,UKF），但是对于非线性、非高斯性较强的问题，UKF的收敛性会急剧下降，甚至出现发散现象。 

为了解决非线性、非高斯系统的估计问题，1993年Gordon等人提出一种基于SIS的Bootstrap滤波方法，即粒子滤波（PF或SIR）算法，其思想是利用一组带有相应权重的样本（或粒子）来描述系统随机变量的后验概率密度，以样本均值代替积分运算，从而获得状态最小方差估计的过程。粒子滤波可完整的反映状态的后验分布，理论上适用于任何的非线性系统。 

尽管如此，粒子滤波存在粒子退化现象，即经过若干次迭代后，除了一个粒子外，几乎所有粒子都只具有非常小的权值，这就意味着大量的计算花在了那些对后验概率估计几乎没有影响的粒子上。 

选择好的重要性分布，是改善粒子退化的方法之一。无迹粒子滤波（UPF）正是利用了这一点，它利用UKF来产生重要性分布函数， 同时利用了当前时刻新的观测值，从而可获得较高的估计精度，但存在计算量大的问题。 

重采样是减弱粒子退化的另一种方法，虽然重采样后的粒子权值不为零，但由于传统重采样将权值大的粒子多次复制，同时抛弃权值较小的粒子，如此下去，有效粒子样本被逐渐耗尽，从而导致新的问题，即粒子样本贫化现象。UPF采用了传统重采样方法，显然无法避免粒子贫化问题的出现。 

发明内容

本发明的目的在于避免上述背景技术中的不足，提供一种能在一定程度上解决样本贫化，降低计算量，提高粒子利用率的基于差分演化的粒子滤波方法。 

本发明的目的是这样实现的：基于差分演化的粒子滤波方法，其特征在于包括步骤如下： 

a.在时间t＝0时，从先验概率分布p(x₀)中随机采样，获得初始化粒子

其中i＝1,2,...,N代表粒子个数，N为大于1的整数，同时设每个初始化粒子的权值为

并计算粒子的均值与方差； 

b.根据当前时刻的观测，利用UKF算法对前一时刻的粒子进行滤波处理，得到重要性函数，并从重要性函数采样得到当前时刻的重要性采样粒子； 

c.计算当前时刻重要性采样粒子的权值并对其归一化； 

d.利用差分演化重采样方法对当前时刻重要性采样粒子进行优化； 

e.对差分演化重采样优化后的带权粒子群加权后就得到了当前时刻的滤波结果； 

f.继续步骤b～e，直至完成所有时刻的滤波处理。 

其中，步骤b中所述的重要性函数为

所述的从重要性函数采样得到当前时刻的重要性采样粒子

过程为： 

其中

表示由UKF算法滤波得到的粒子估计值的均值，

表示由UKF算法估计得到的粒子估计值的方差，

表示均值为方差为

的高斯分布函数或正态分布函数。 

其中，步骤c中所述的对粒子权值

归一化具体为： 

为归一化后的粒子权值。 

其中，步骤d具体包括以下步骤： 

d1.把步骤b中产生的采样粒子作为差分演化算法中的种群样本，构成当前时刻的原始种群，把相应的粒子权重作为当前时刻原始种群样本的适应函数； 

d2.对当前时刻原始种群中的样本做差分变异，得到当前时刻变异种群； 

d3.对当前时刻的原始种群和变异种群中的样本做差分交叉，得到当前时刻试验种群，并计算其适应函数； 

d4.对当前时刻的原始种群和试验种群中的样本做差分选择，得到下一代的原始种群； 

d5.继续步骤d2～d4，直到种群可以满足设定的终止条件或者种群繁殖代数达到设定的最大值为止。 

其中，步骤d2中所述的对当前时刻原始种群P_X,m＝(X_i,m)中的样本X_i,m＝(x_d,i,m)_{d＝1,2,...,D}做差分变异，得到当前时刻变异种群，具体实现如下： 

从初始种群中随机选择两个样本做比例差分，再与种群中的随机样本相加得到变异种群P_V,m＝(V_i,m)中的变异样本V_i,m＝(v_d,i,m)_{d＝1,2,...,D}： 

V_i,m＝F₁·X_r1,m+F₂·(X_r2,m-X_r3,m) 

其中，i=1,2,...,NP为种群样本个数，种群样本个数与粒子数相同，即NP＝N，每个样本由D维样本数值构成；比例因子F₁,F₂∈(0,1+)用来控制样本发生变异的速度；X_r2,m-X_r3,m为差异向量，反映群体的差异度，r1，r2，r3是每次种群变异时随机产生的序号，且i≠r1≠r2≠r3；m表示种群的进化代数。 

其中，步骤d3中所述的对当前时刻的原始种群和变异种群中的样本做差分交叉，得到当前时刻试验种群，具体实现如下： 

交叉是对原种群P_X,m和变异种群P_V,m中样本数值的融合交叉，以产生试验种群P_U,m＝(U_i,m)： 

其中，交叉概率Cr∈[0,1]，用来控制变异样本数值对试验样本数值产生贡献的大小，若随机数rand_d(0,1)小于Cr，试验样本数值由变异样本数值贡献，否则由原样本数值贡献；r(D)是[1,D]间的随机数。 

其中，步骤d4具体实现如下： 

若U_i,m的适应函数值大于X_i,m的适应函数值，在下一代种群中就用U_i,m来取代X_i,m，否则保留X_i,m： 

其中，ψ(U_i,m)为U_i,m的适应函数；ψ(X_i,m)为X_i,m的适应函数。 

本发明与现有粒子滤波方法相比，有益效果主要体现在： 

本发明充分利用最新观测值，采用UKF算法来产生重要性分布，减小采样粒子与真实后验分布产生样本的偏差；对重要性分布产生的采样粒子不再做传统重采样操作，而是采用差分演化重采样方法对采 样粒子进行优化选择，有效缓解了粒子滤波中的粒子退化和贫化问题，提高了粒子的利用率，具有更高的估计精度。 

附图说明

图1是本发明所述基于差分演化的粒子滤波方法的流程图； 

图2是本发明与已有粒子滤波方法的状态估计性能对比图； 

图3是本发明与已有粒子滤波方法在不同粒子数时状态估计的均方根误差均值比较图。 

具体实施方式

以下结合附图1-3举例对本发明做进一步的说明。 

本发明基于差分演化的粒子滤波方法的基本步骤如下： 

a.在时间t＝0时，从先验概率分布p(x₀)中随机采样，获得初始化粒子

其中i＝1,2,...,N代表粒子个数，N为大于1的整数，同时设每个初始化粒子的权值为

并计算粒子的均值

与方差 

${\stackrel{\‾}{x}}_0^{\left(i\right)}=E\left[x_0^{\left(i\right)}\right]$

$P_0^{\left(i\right)}=E\left[(x_0^{\left(i\right)}-{\stackrel{\‾}{x}}_0^{\left(i\right)}){(x_0^{\left(i\right)}-{\stackrel{\‾}{x}}_0^{\left(i\right)})}^T\right]$

b.根据当前时刻的观测，利用UKF算法对前一时刻的粒子进行滤波处理，得到重要性函数，并从重要性函数采样得到当前时刻的重要性采样粒子； 

该步骤中所述的重要性函数为

所述的从重要性函数采样得到当前时刻的重要性采样粒子过程为：

其中

表示由UKF算法滤波得到的粒子估计值的均值，

表示由UKF算法估计得到的粒子估计值的方差，

表示均值为方差为

的高斯分布函数或正态分布函数。 

该步骤具体实现为： 

b1.UKF算法初始化： 

将粒子

的维数扩展，得到UKF算法中的初始粒子均值

及方差

${\stackrel{\‾}{x}}_0^{\left(i\right)a}=E\left[x_0^{\left(i\right)a}\right]=E\left[\begin{array}{ccc}{\stackrel{\‾}{x}}_0^{\left(i\right)}& 0& 0\end{array}\right]$

$P_0^{\left(i\right)a}=E\left[(x_0^{\left(i\right)a}-{\stackrel{\‾}{x}}_0^{\left(i\right)a}){(x_0^{\left(i\right)a}-{\stackrel{\‾}{x}}_0^{\left(i\right)a})}^T\right]=\left[\begin{array}{ccc}{P}_0^{\left(i\right)}& 0& 0\\ 0& Q& 0\\ 0& 0& R\end{array}\right]$

其中，Q为系统过程噪声误差的协方差阵，R为测量噪声误差的协方差阵； 

b2.计算每个粒子对应的2n_a+1个sigma点

${\mathrm{\χ}}_{t-1}^{\left(i\right)a}={\left[\begin{array}{ccc}{\stackrel{\‾}{x}}_{t-1}^{\left(i\right)a}& {\stackrel{\‾}{x}}_{t-1}^{\left(i\right)a}+\sqrt{(n_a+\mathrm{\λ})P_{t-1}^{\left(i\right)}}& {\stackrel{\‾}{x}}_{t-1}^{\left(i\right)a}-\sqrt{(n_a+\mathrm{\λ})P_{t-1}^{\left(i\right)}}\end{array}\right]}_{({2n}_a+1)\×1}$

sigma点的权重系数：

计算方差的权重系数： $w_0^{\left(c\right)}=\mathrm{\λ}/(n_a+\mathrm{\λ})+(1+{\mathrm{\α}}^2+\mathrm{\β})$

$w_j^{\left(m\right)}=w_j^{\left(c\right)}=1/\left\{2(n_a+\mathrm{\λ})\right\},j=\mathrm{1,2},...,{2n}_a$

其中，n_a为粒子的维数，λ为比例参数，α为正比例参数，β为影响方差计算的参数； 

b3.计算时间更新后的粒子

及其方差

并计算粒子

对应的测量预测

${\mathrm{\χ}}_{t|t-1}^{\left(i\right)x}=f({\mathrm{\χ}}_{t-1}^{\left(i\right)x},{\mathrm{\χ}}_{t-1}^{\left(i\right)v}),{\stackrel{\‾}{x}}_{t|t-1}^{\left(i\right)}=\underset{j=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{w}_j^{\left(m\right)}{\mathrm{\χ}}_{j,t|t-1}^{\left(i\right)x}$

$P_{t|t-1}^{\left(i\right)}=\underset{j=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{w}_j^{\left(c\right)}{[{\mathrm{\χ}}_{j,t|t-1}^{\left(i\right)x}-{\stackrel{\‾}{x}}_{t|t-1}^{\left(i\right)}][{\mathrm{\χ}}_{j,t|t-1}^{\left(i\right)x}-{\stackrel{\‾}{x}}_{t|t-1}^{\left(i\right)}]}^T$

$y_{t|t-1}^{\left(i\right)}=g({\mathrm{\χ}}_{t-1}^{\left(i\right)x},{\mathrm{\χ}}_{t-1}^{\left(i\right)n}),{\stackrel{\‾}{y}}_{t|t-1}^{\left(i\right)}=\underset{j=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{w}_j^{\left(m\right)}{y}_{j,t|t-1}^{\left(i\right)}$

其中，f(·)为非线性系统的状态转移矩阵函数，g(·)为非线性系统的观测矩阵函数。 

b4.测量更新，计算新粒子及其方差

$P_{y_t{y}_t}=\underset{j=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{w}_j^{\left(c\right)}{[y_{j,t|t-1}^{\left(i\right)}-{\stackrel{\‾}{y}}_{t|t-1}^{\left(i\right)}][y_{j,t|t-1}^{\left(i\right)}-{\stackrel{\‾}{y}}_{t|t-1}^{\left(i\right)}]}^T$

$P_{x_t{y}_t}=\underset{i=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{w}_j^{\left(c\right)}{[{\mathrm{\χ}}_{i,t|t-1}^{\left(i\right)x}-{\stackrel{\‾}{x}}_{t|t-1}^{\left(i\right)}][y_{j,t|t-1}^{\left(i\right)}-{\stackrel{\‾}{y}}_{t|t-1}^{\left(i\right)}]}^T$

$K_t=P_{x_t{y}_t}{P}_{y_t{y}_t}^{-1}$

${\stackrel{\‾}{x}}_t^{\left(i\right)}={\stackrel{\‾}{x}}_{t|t-1}^{\left(i\right)}+K_t(y_t-{\stackrel{\‾}{y}}_{t|t-1}^{\left(i\right)}),$ y_t是观测量 

b5.构造重要性函数

并从中采样： 

其中，

表示均值为

方差为

的高斯分布函数或正态分布函数； 

c.计算当前时刻重要性采样粒子的权值并对其归一化。 

计算每个粒子的权值

${\mathrm{\ω}}_t^{\left(i\right)}\∝{\mathrm{\ω}}_{t-1}^{\left(i\right)}p\left(y_t\right|{\hat{x}}_t^{\left(i\right)})$

步骤c中所述的对粒子权值归一化具体为：

其中，

为归一化后的粒子权值。 

d.利用差分演化重采样方法对当前时刻重要性采样粒子进行优化，具体步骤为： 

d1.把步骤b中产生的采样粒子

作为DE算法中的样本X_i,m＝(x_d,i,m)_{d＝1,2,...,D}，每个样本由D维样本数值构成；所有样本构成初始种群P₁＝(X_i,1),i=1,2,...,NP，此时N=NP；把相应的粒子权重

作为种群的适应函数ψ(X_i,m)，其中m表示种群的进化代数。 

d2.对种群中的样本做差分变异： 

从初始种群中随机选择两个样本做比例差分，再与种群中的随机 样本相加得到变异种群P_V,m＝(V_i,m)中的变异样本V_i,m＝(v_d,i,m)_{d＝1,2,...,D}： 

V_i,m＝F₁·X_r1,m+F₂·(X_r2,m-X_r3,m) 

其中，比例因子F₁,F₂∈(0,1+)用来控制样本发生变异的速度，X_r2,m-X_r3,m为差异向量，它反映群体的差异度，r1，r2，r3是每次种群变异时随机产生的序号，且i≠r1≠r2≠r3。 

d3.对种群中的样本做差分交叉： 

交叉是对原种群P_X,m和变异种群P_V,m中样本数值的融合交叉，以产生试验种群P_U,m＝(U_i,m)： 

其中，交叉概率Cr∈[0,1]，用来控制变异样本数值对试验样本数值产生贡献的大小，若随机数rand_d(0,1)小于Cr，试验样本数值由变异样本数值贡献，否则由原样本数值贡献。r(D)是[1,D]间的随机数，用来防止试验种群对变异种群进行复制，也即为了保证试验种群中至少有一个样本数值来自于变异种群。 

d4.对种群中的样本做差分选择： 

若U_i,m的目标函数值大于X_i,m的目标函数值，在下一代种群中就用U_i,m来取代X_i,m，否则保留X_i,m： 

d5.继续步骤d2～d4，直到种群可以满足终止条件或者种群繁殖代数达到设定的最大值为止。 

e.输出系统的状态估计值：

其中是对非线性系统状态变量的估计值，

是采样粒子经差分演化算法优化后的粒子，

是优化粒子

的权值。 

f.继续步骤b～e，直至完成所有时刻的滤波处理。 

实施实例： 

利用惯用的非线性模型对滤波性能进行验证，其状态方程和测量方程为： 

$\left\{\begin{array}{c}{x}_k=1+\sin \left(0.04\mathrm{\&pi;k}\right)+{0.5x}_{k-1}+u_{k-1}\\ z_k=\left\{\begin{array}{cc}{0.2x}_k^2+v_k& 1\&le;k\&le;30\\ {0.5x}_k-2+v_k& 30<k\&le;T\end{array}\right.\end{array}\right.$

其中过程噪声u_k服从Gamma(3,2)分布，测量噪声v_k服从

的高斯分布，观测时间T＝50s。 

试验中，分别采用SIR，正则粒子滤波（RPF），UPF与本文方法进行比较。本发明方法的参数设置为F＝0.9,Cr＝0.6，种群最大繁殖代数为G＝20，UKF中的参数设置与参考文献1：Merwe,R.,Doucet,A.,Freitas,N.,and Wan,E.:‘The unscented particle filter’,Technical report CUED/F-INFENG/TR380,Cambridge University Engineering Department,2000,pp.1-40中的设置相同。 

图2是几种滤波方法经500次蒙特卡洛试验得到的估计均方根误差（RMSE）图，粒子数N=10。其中DSIR是将SIR中的传统重采样替换为差分演化后的滤波方法，可以看出只替换重采样部分，算法的估计精度高于SIR方法的精度，说明利用差分演化进行重采样是有效的；DPF1，DPF2，DPF3，DPF4都属于本发明基于差分演化的粒子滤波方法的范畴，区别在于所采用的差分变异方式不同，它们分别采用的差分变异为：V_i,m＝F₁·X_r1,m+F₂·(X_r2,m-X_r3,m)，V_l,g＝X_r1,g+F×(X_r2,g-X_r3,g+X_r4,g-X_r5,g)，V_l,g＝X_best,g+F×(X_r2,g-X_r3,g)，V_l,g＝X_best,g+F×(X_r2,g-X_r3,g+X_r4,g-X_r5,g)，其中X_best,g是当前种群中适应函数值最大的样本。 

从图2可以看出，本发明方法的估计均方根误差明显小于SIR、 RPF、UPF的估计均方根误差。在DPF1，DPF2，DPF3，DPF4中，DPF3的估计精度最高。 

图3是在不同粒子数目（N=10，100，200）情况下，几种滤波算法的性能比较，其中蒙特卡洛试验次数为200。可以看出，DPF3只使用10粒子的估计精度高于SIR、RPF使用200个粒子的估计精度，DPF3使用100粒子的估计精度好于UPF使用200个粒子的估计精度。相比其他方法，DPF3使用较少的粒子就能实现对非线性系统的精确估计，可见本发明方法的粒子利用率更高。 

正是由于差分演化的引入，避免了传统重采样对粒子的简单复制与抛弃，而是采用差分演化算法对粒子集进行进化寻优，有效缓解了粒子滤波中的粒子退化和贫化问题，提高了粒子的利用率，具有更高的估计精度。 

Claims (7)

1.基于差分演化的粒子滤波方法，其特征在于包括步骤如下：

a.在时间t＝0时，从先验概率分布p(x₀)中随机采样，获得初始化粒子

其中i＝1,2,...,N代表粒子个数，N为大于1的整数，同时设每个初始化粒子的权值为

并计算粒子的均值与方差；

b.根据当前时刻的观测，利用UKF算法对前一时刻的粒子进行滤波处理，得到重要性函数，并从重要性函数采样得到当前时刻的重要性采样粒子；

c.计算当前时刻重要性采样粒子的权值并对其归一化；

d.利用差分演化重采样方法对当前时刻重要性采样粒子进行优化；

e.对差分演化重采样优化后的带权粒子群加权后就得到了当前时刻的滤波结果；

f.继续步骤b～e，直至完成所有时刻的滤波处理。

2.根据权利要求1所述的基于差分演化的粒子滤波方法，其特征在于，步骤b中所述的重要性函数为

所述的从重要性函数采样得到当前时刻的重要性采样粒子过程为：

其中

表示由UKF算法滤波得到的粒子估计值的均值，

表示由UKF算法估计得到的粒子估计值的方差，

表示均值为

方差为

的高斯分布函数或正态分布函数。

3.根据权利要求1所述的基于差分演化的粒子滤波方法，其特征在于，步骤c中所述的对粒子权值

归一化具体为：

其中，

为归一化后的粒子权值。

4.根据权利要求1所述的基于差分演化的粒子滤波方法，其特征在于：步骤d具体包括以下步骤：

d1.把步骤b中产生的采样粒子作为差分演化算法中的种群样本，构成当前时刻的原始种群，把相应的粒子权重作为当前时刻原始种群样本的适应函数；

d2.对当前时刻原始种群中的样本做差分变异，得到当前时刻变异种群；

d3.对当前时刻的原始种群和变异种群中的样本做差分交叉，得到当前时刻试验种群，并计算其适应函数；

d4.对当前时刻的原始种群和试验种群中的样本做差分选择，得到下一代的原始种群；

d5.继续步骤d2～d4，直到种群可以满足设定的终止条件或者种群繁殖代数达到设定的最大值为止。

5.根据权利要求4所述的基于差分演化的粒子滤波方法，其特征在于，步骤d2中所述的对当前时刻原始种群P_X,m＝(X_i,m)中的样本X_i,m＝(x_d,i,m)_{d＝1,2,...,D}做差分变异，得到当前时刻变异种群，具体实现如下：

从初始种群中随机选择两个样本做比例差分，再与种群中的随机样本相加得到变异种群P_V,m＝(V_i,m)中的变异样本V_i,m＝(v_d,i,m)_{d＝1,2,...,D}：

V_i,m＝F₁·X_r1,m+F₂·(X_r2,m-X_r3,m)

其中，i=1,2,...,NP为种群样本个数，种群样本个数与粒子数相同，即NP＝N，每个样本由D维样本数值构成；比例因子F₁,F₂∈(0,1+)用来控制样本发生变异的速度；X_r2,m-X_r3,m为差异向量，反映群体的差异度，r1，r2，r3是每次种群变异时随机产生的序号，且i≠r1≠r2≠r3；m表示种群的进化代数。

6.根据权利要求4所述的基于差分演化的粒子滤波方法，其特征在于：步骤d3中所述的对当前时刻的原始种群和变异种群中的样本做差分交叉，得到当前时刻试验种群，具体实现如下：

交叉是对原种群P_X,m和变异种群P_V,m中样本数值的融合交叉，以产生试验种群P_U,m＝(U_i,m)：

其中，交叉概率Cr∈[0,1]，用来控制变异样本数值对试验样本数值产生贡献的大小，若随机数rand_d(0,1)小于Cr，试验样本数值由变异样本数值贡献，否则由原样本数值贡献；r(D)是[1,D]间的随机数。

7.根据权利要求6所述的基于差分演化的粒子滤波方法，其特征在于，步骤d4具体实现如下：

若U_i,m的适应函数值大于X_i,m的适应函数值，在下一代种群中就用U_i,m来取代X_i,m，否则保留X_i,m：

其中，ψ(U_i,m)为U_i,m的适应函数；ψ(X_i,m)为X_i,m的适应函数。

Images