CN101055563A - 基于多建议分布的粒子滤波方法 - Google Patents

Abstract

本发明涉及基于多建议分布的粒子滤波方法，属于信号处理、人工智能和计算机视觉领域。本发明在粒子滤波器的框架内使用多种建议分布：先验概率分布、扩展卡尔曼滤波器、unscented卡尔曼滤波器等，采用分治采样策略来处理样本粒子，把总的粒子数分成多个部分，分别从不同的建议分布中提取，这样可以降低粒子滤波器的运行时间，提高运行效率，而不损失粒子滤波器的估计精度。用户可以根据自己对时间和精度的需求进行参数设置。在涉及非线性滤波问题的领域具有广泛的应用前景。

Description

基于多建议分布的粒子滤波方法

技术领域

本发明涉及基于多建议分布的粒子滤波方法，要求保护的技术方案属于信号处理、人工智能和计算机视觉领域。

背景技术

非线性滤波问题在很多领域的问题中都涉及到，其中包括信号处理、金融、人工智能和计算机视觉等等。解决非线性滤波问题最为普遍的一种方法是使用扩展卡尔曼滤波器(EKF)，但是EKF只适用于弱非线性系统，对于强非线性系统很容易导致发散。还有一种解决非线性滤波问题的办法是无迹卡尔曼滤波器(UKF)。UKF和EKF不同，它不使用局部线性化技术，而直接使用系统的非线性方程进行运算，从而能够避免局部线性化所引入的误差，避免在强非线性系统中出现发散。但是，UKF的使用不适用于一般的非高斯模型。

另一种解决非线性滤波问题的方法就是粒子滤波器(Particle Filter)。它的基本思想就是使用一组带有权值的粒子(样本)集合来近似表示解决问题时需要的系统后验概率密度。粒子滤波器自提出以来已经广泛应用到信号处理、金融、人工智能、计算机视觉以及机器人等领域中，成为各国研究者关注的焦点。

粒子滤波算法中的一项关键技术是建议分布的选择。选择好的建议分布能够有效的避免退化问题的影响，提高系统的状态估计精度。目前，常见的建议分布有先验概率分布(Transition Prior，TP)、EKF、UKF等。以这三种建议分布为基础的粒子滤波器，分别被称为：基本粒子滤波器(Generic PF)、扩展卡尔曼粒子滤波器(EKPF)和无迹粒子滤波器(UPF)。但是先验概率分布没有考虑当前时刻新的观测值的影响，从而影响了粒子滤波器估计精度；而EKF采用局部线性化技术，引入了过多的截断误差，对估计精度产生影响；以UKF为建议分布的粒子滤波器具有较高的时间耗费。

发明内容

针对上述问题，本发明提出一种基于多建议分布的粒子滤波算法。

该算法包括以下基本步骤：

一.初始状态下，首先从初始分布p(x₀)中提取一组(N个)样本粒子，并设置其初始的均值和方差。

二.对于c*N个粒子：

(1)首先使用UKF进行粒子传递和更新，获得关于该时刻状态的估计x_k ukf ⁽ⁱ⁾.

(2)使用EKF重复粒子更新过程，并以第(1)步中得到的状态估计作为输入。得到该时刻最终的估计量 x_k ⁽ⁱ⁾和

(3)根据(2)中得到的估计量，构成建议分布N( x_k ⁽ⁱ⁾，

)，从中抽取新的样本粒子，并赋予权值w_k ⁱ。

三.对于剩余的(1-c)*N个粒子，使用非线性系统方程进行更新，得到状态估计预测量 x_k|k-1 ⁽ⁱ⁾，并从先验概率分布p(x_k|k-1 ⁽ⁱ⁾|x_k-1 ⁽ⁱ⁾)中抽取粒子，赋予权值w_k ⁱ。

四.对所有粒子的权值进行归一化处理。

五.执行再采样过程，对得到的新的粒子集合中的每一个粒子赋予相同的权值1/N。

六.输出

七.如继续迭代，转到二，否则结束。

本发明采用分治采样策略，将系统所需粒子分成两部分，分别由不同的建议分布中提取。

一部分(假设占百分比为c)由混合的建议分布提取，即：由UKF和EKF组成的混合建议分布。其余的一部分(1-c)，由先验概率分布p(x_k|x_k-1)中提取。以此为基础的粒子滤波算法不仅能够产生较高的估计精度，而且能够降低粒子滤波器的运行时间。用户可以根据自身对实际要求的不同，调整参数c，以满足他们对运行效率和精度的不同要求。

假设非线性动态系统的状态空间模型为：

                       x_k＝f_k(x_k-1，v_k-1)

                       y_k＝h_k(x_k，u_k)

其中x_k表示系统在k时刻所处的状态，y_k表示k时刻的观测值。函数f和h表示系统的状态转移和测量模型。v_k和u_k分别表示系统噪声与测量噪声。

本发明提供的多建议分布粒子滤波器采用了EKF、UKF和TP作为建议分布，将总的粒子数分成两部分，分别由不同的建议分布中抽取，具体实现过程包括以下步骤：

第一步.k＝0时，初始化一组粒子作为初始状态下的粒子集合，并设置初始的均值和协方差矩阵，设置方法同UPF算法中的初始设置。

第二步.对于c*N个粒子，

1.首先使用UKF更新粒子：

(1)选择一组sigma点，构成sigma矩阵。

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

(2)根据动态系统的非线性方程，将这些sigma点向前传递，得到系统状态的均值的预测量 x_k|k-1ukf ⁽ⁱ⁾和协方差的预测量P_k|k-1ukf ⁽ⁱ⁾，以及与测量值有关的预测量y_k|k-1ukf ⁽ⁱ⁾(预测)：

(3)获取当前时刻的测量值y_k之后，并根据该测量值对(2)中得到的一步预测值(更新)：

${{\stackrel{\‾}{x}}_k^{\left(i\right)}}_{\mathrm{ukf}}={{\stackrel{\‾}{x}}_{k|k-1}^{\left(i\right)}}_{\mathrm{ukf}}{+K}_k(y_k-{{\stackrel{\‾}{y}}_{k|k-1}^{\left(i\right)}}_{\mathrm{ukf}})$

其中 $K_k=P_{x_k{y}_k}{P}_{{\tilde{y}}_k{\tilde{y}}_k}^{-1}$ (卡尔曼增益)

$P_{{\tilde{y}}_k{\tilde{y}}_k}=\underset{j=0}{\overset{n_a^2}{\mathrm{\Σ}}}{W}_j^{\left(c\right)}[Y_{j,k|k-1}^{\left(i\right)}-{\stackrel{\‾}{y}}_{k|k-1}^{\left(i\right)}]{[Y_{j,k|k-1}^{\left(i\right)}-{\stackrel{\‾}{y}}_{k|k-1}^{\left(i\right)}]}^T$

$P_{{\tilde{x}}_k{\tilde{y}}_k}=\underset{j=0}{\overset{n_a^2}{\mathrm{\Σ}}}{W}_j^{\left(c\right)}[{\mathrm{\χ}}_{j,k|k-1}^{\left(i\right)}-{\stackrel{\‾}{x}}_{k|k-1}^{\left(i\right)}]{[Y_{j,k|k-1}^{\left(i\right)}-{\stackrel{\‾}{y}}_{k|k-1}^{\left(i\right)}]}^T$

这样就获得了估计量 x_k ukf ⁽ⁱ⁾。

其中，x^a＝[x^T v^T u^T]^T为扩张(augmented)状态变量，由系统状态x，系统噪声v和测量噪声u组成。χ^a＝[(χ^x)^T (χ^v)^T (χ^u)^T]^T为sigma矩阵。n_a＝n_x+n_v+n_u为x^a的维数。χ_k-1 ^(i)a表示k-1时刻由sigma点构成的sigma矩阵，χ_j，k|k-1 ⁽ⁱ⁾为χ_k|k-1 ⁽ⁱ⁾的第j个元素。W_i ^(m)和W_i ^(c)表示sigma向量的权值，分别用来求均值和协方差估计，其计算方法在文献(Rudolph van der Merwe等,The Unscented Particle Filter)中。λ＝α²(n_x+k)-n_x是一个尺度调节因子，α决定了选择的sigma点在其均值 x附近的分布情况，通常将α设置为一个很小的正值(如0.001)。k是次级尺度调节因子，通常设置为0，β是用来结合关于x的分布的先验知识(对于高斯分布，β的最佳取值为2)。

是矩阵(n_x+λ)P_x平方根的第i行，

2.使用EKF更新粒子，并以 x_k ukf ⁽ⁱ⁾作为输入，即：令 $x_{k-1}^{\left(i\right)}={{\stackrel{\‾}{x}}_k^{\left(i\right)}}_{\mathrm{ukf}}.$

(1)传递粒子，并得到相关的一步预测量：

${{\stackrel{\‾}{x}}_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}f\left({\stackrel{\‾}{x}}_{k-1}^{\left(i\right)}\right)=f\left({{\stackrel{\‾}{x}}_k^{\left(i\right)}}_{\mathrm{ukf}}\right)$

${P_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}=F_k^{\left(i\right)}{P}_{k-1}^{\left(i\right)}{F}_k^{T\left(i\right)}+G_k^{\left(i\right)}{Q}_k{G}_k^{T\left(i\right)}$

(2)计算系统状态转移模型的雅可比矩阵F和G，以及测量模型的雅可比矩阵H和U.并计算出卡尔曼增益。

$K_k={P_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}{H}_k^{T\left(i\right)}{[U_k^{\left(i\right)}{R}_k{U}_k^{T\left(i\right)}+H_k^{\left(i\right)}{P_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}{H}_k^{T\left(i\right)}]}^{-1}$

(3)根据当前的时刻的观测值y_k，修正(1)中的一步预测量，得到修正的的估计量：

${{\hat{P}}_k^{\left(i\right)}}_{\mathrm{ekf}}={P_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}-K_k{H}_k^{\left(i\right)}{P_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}$

${{\stackrel{\‾}{x}}_k^{\left(i\right)}}_{\mathrm{ekf}}={{\stackrel{\‾}{x}}_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}+{{\hat{P}}_k^{\left(i\right)}}_{\mathrm{ekf}}{H}_k^{T\left(i\right)}{R}_k^{-1}(y_k-h\left({{\stackrel{\‾}{x}}_{k|k-1}^{\left(i\right)}}_{\mathrm{ekf}}\right))$

(4)令 ${\stackrel{\‾}{x}}_k^{\left(i\right)}={{\stackrel{\‾}{x}}_k^{\left(i\right)}}_{\mathrm{ekf}},{\hat{P}}_k^{\left(i\right)}={{\hat{P}}_k^{\left(i\right)}}_{\mathrm{ekf}}$

3.从它们所组成的建议分布中抽取新的粒子：

${\hat{x}}_k^{\left(i\right)}~q\left(x_k^{\left(i\right)}\right|x_{0:k-1}^{\left(i\right)},y_{1:k})=N({\stackrel{\‾}{x}}_k^{\left(i\right)},{\hat{P}}_k^{\left(i\right)})$

并赋予其权值w_k ⁱ。

其中，Q和R分别表示系统噪声和测量噪声的方差

第三步.对于剩余的(1-c)*N个粒子，使用先验概率分布进行更新。

1.直接根据系统的非线性方程计算出状态估计量：

${\stackrel{\‾}{x}}_{k|k-1}^{\left(i\right)}=f\left(x_{k-1}^{\left(i\right)}\right)$

2.从先验分布密度函数p(x_k|k-1 ⁽ⁱ⁾|x_k-1 ⁽ⁱ⁾)中抽取粒子，即：

并为其赋予权值w_k ⁱ。

第四步.对所有的粒子的权值进行归一化处理：

$w_k^i=w_k^i/\mathrm{\Σ}{w}_k^j$

第五步.执行再采样过程：

去除那些权值较小的粒子，复制权值较大的粒子，得到一组近似服从后验概率分布p(x_0:k|z_1:k)的粒子集合，并为其中的每个粒子赋予相同的权值1/N。

第六步.输出阶段。

第七步.k＝k+1，如继续迭代，转到第二步循环，否则结束。

本发明具有以下优点：

1.采用多建议分布方案，融合了各种独立建议分布在粒子滤波器框架内的优点。

2.采用分治采样策略，在粒子滤波器的运行时间和估计精度等方面作折衷，用户可以根据自己的需求适当调整参数，获得对运行时间或精度的要求。

附图说明

图1-多建议分布的粒子滤波流程图；

图2-在不同的c取值下，多建议分布粒子滤波器算法与其它粒子滤波器算法的估计均方误差对比曲线；

图3-在不同的c取值下，多建议分布粒子滤波器算法与其它粒子滤波器算法的运行时间对比曲线；

图4-原始信号曲线；

图5-使用本发明进行滤波后的曲线；

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

本发明的算法包括以下基本步骤：

一.初始状态下，首先从初始分布p(x₀)中提取一组(N个)样本粒子，并设置其初始的均值和方差。

二.对于c*N个粒子：

(1)首先使用UKF进行粒子传递和更新，获得关于该时刻状态的估计x_k ukf ⁽ⁱ⁾。

(2)使用EKF重复粒子更新过程，并以第(1)步中得到的状态估计作为输入。得到该时刻最终的估计量 x_k ⁽ⁱ⁾和

(3)根据(2)中得到的估计量，构成建议分布N( x_k ⁽ⁱ⁾， )，从中抽取新的样本粒子，并赋予权值w_k ⁱ。

三.对于剩余的(1-c)*N个粒子，使用非线性系统方程进行更新，得到状态估计预测量 x_k|k-1 ⁽ⁱ⁾，并从先验概率分布p(x_k|k-1 ⁽ⁱ⁾|x_k-1 ⁽ⁱ⁾)中抽取粒子，赋予权值w_k ⁱ。

四.对所有粒子的权值进行归一化处理。

五.执行再采样过程，对得到的新的粒子集合中的每一个粒子赋予相同的权值1/N。

六.输出

七.如继续迭代，转到二，否则结束。

实施例1.通过一个非线性动态系统的状态空间模型如下：

x_k＝1+sin(0.04π(k-1))+0.5x_k-1+v_k-1      (状态转移模型)

$z_k=\left\{\begin{array}{cc}0.2x_k^2+u_k& k\≤30\\ 0.5x_k-2+u_k& k>30\end{array}\right.$ (测量模型)

对本发明与其它几种粒子滤波器算法进行对比。使用粒子数目N＝200，观测时间T＝60。程序在CPU为2.67GHz，内存为1GB的微型计算机上运行100次。

每次运行程序的输出为粒子集合的均值： ${\hat{x}}_k=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^j.$ 每一次运行中，均方误差的计算公式为：

$\mathrm{MSE}={\left(\frac{1}{T}\underset{k=1}{\overset{T}{\mathrm{\Σ}}}{({\hat{x}}_k-x_k)}^2\right)}^{1/2}$

取c＝0.3，在时间耗费和精度方面进行对比时，本发明算法产生的估计均方误差比UPF低80％，运行时间比UPF低约50％。在c＝0.4时，本发明算法的均方误差比UPF低72％，运行时间比UPF低44％。

图2所示，c在取不同的值时，本发明与其它粒子滤波器算法均方误差的对比，图3为相应的运行时间对比。随着c值的增加，本发明提出的算法在运行时间上也随之增加，但是本发明所产生的均方误差总体上是呈现递减的趋势。

实施例2.将本发明实施于钻井系统中的随钻测井仪。随钻测井仪器分为井上部分和井下部分。井下部分的测量信号经泥浆传送到地面。井上部分首先要对接收到的井下信号进行滤波，然后是解码并算出物理参数值。本发明的作用就是将带噪声的编码信号进行滤波，得到去噪的信号。将本发明实施于仪器的井上部分。本发明能够得到理想的滤波效果(图4，图5)，为之后数据的处理提供了基础。

Claims (1)

1、基于多建议分布的粒子滤波方法，其特征在于：包括以下步骤：

1)首先初始化一组(N个)样本粒子，并设置其初始的均值和方差；

2)对于c*N个粒子：

(1)首先使用UKF进行粒子传递和更新，获得关于该时刻状态的估计x_k ⁽ⁱ⁾ _ukf；

(2)使用EKF重复粒子更新过程，并以第(1)步中得到的状态估计作为输入。得到该时刻最终的估计量 x_k ⁽ⁱ⁾和

(3)根据(2)中得到的估计量，构成建议分布N( x_k ⁽ⁱ⁾，

)，并从中抽取新的样本粒子，并赋予权值w_k ⁱ；

3)对于剩余的(1-c)*N个粒子，使用非线性系统方程进行更新，得到状态估计预测量 x_k|k-1 ⁽ⁱ⁾，并从先验概率分布p(x_k|k-1 ⁽ⁱ⁾|x_k-1 ⁽ⁱ⁾)中抽取粒子，赋予权值w_k ⁱ；

4)对所有粒子的权值进行归一化处理；

5)执行再采样过程，对得到的新的粒子集合中的每一个粒子赋予相同的权值1/N；

6)输出；

7)如继续迭代，转到二，否则结束。

Images