CN101436251A - 基于拟蒙特卡罗采样的并行高斯粒子滤波方法 - Google Patents

Abstract

本发明公开了一种基于拟蒙特卡罗采样的并行高斯粒子滤波方法，属于信号处理技术领域，主要解决目前基于蒙特卡罗采样的滤波方法在样本数较少时对非线性动态系统状态估计性能降低的问题。其滤波步骤包括并行拟蒙特卡罗序列产生步骤和并行高斯粒子滤波步骤，即首先采用跳跃拟蒙特卡罗序列产生方法，并行产生后续滤波所需的P条拟蒙特卡罗随机样本序列；然后分别通过P个执行单元，将每条拟蒙特卡罗随机样本序列转换为服从指定高斯分布的拟高斯样本序列；同时对非线性动态系统的状态进行时间更新和状态更新；最后对所有执行单元进行综合处理，完成对非线性动态系统的滤波。本发明具有滤波性能高和实时性好的优点，可用于信号处理、自动化和人工智能领域。

Description

基于拟蒙特卡罗采样的并行高斯粒子滤波方法

技术领域

本发明属于信号处理技术领域，涉及非线性滤波，特别是一种并行高斯粒子滤波方法，用于非线性动态系统的状态估计。

背景技术

非线性滤波技术广泛应用于信号处理、自动化、计算机视觉和人工智能等许多领域。扩展卡尔曼滤波EKF是非线性滤波的经典方法。但由于这种扩展卡尔曼滤波是把非线性函数用它的二阶泰勒展开式来代替，会产生近似误差，精度较低，在应用上容易造成滤波发散的情况。无迹卡尔曼滤波UKF是非线性滤波的另一种方法，无迹卡尔曼滤波不用线性化方法逼近非线性函数，而是用有限个能完全描述状态的均值和方差的样本点来表示状态的概率密度，通过直接应用状态和量测方程来映射这些样本点，加权求和得到更新的均值和方差。虽然无迹卡尔曼滤波比扩展卡尔曼滤波有更高的精度，但是由于无迹卡尔曼滤波和扩展卡尔曼滤波都是线性卡尔曼滤波的变形或改进形式，因此应用上也受线性卡尔曼滤波的制约，在对非高斯系统进行状态估计时性能会变差。

成熟有效的粒子滤波是英国学者Gordon等人于1993年提出的，其基本思想是用一组带权的随机样本表示状态的后验概率密度。粒子滤波没有卡尔曼的限制，可以处理非线性和非高斯的系统，并且容易实现并行化。粒子滤波自提出以来，广泛应用于信号处理、自动化、计算机视觉和人工智能等领域对非线性非高斯系统进行状态估计，成为研究的热点。该基本粒子滤波方法的不足之处是存在粒子退化，一般采用称为“重采样”的算法对其进行改进，但是这种重采样不仅占用了很多计算机资源并且不易于并行实现。

Jayesh H.Kotecha等人于2003年提出高斯粒子滤波GPF，该高斯粒子滤波方法在高斯滤波的基础上用重要性采样方法计算更新高斯模型的均值和协方差矩阵两个参数。GPF要优于传统的高斯滤波器，较粒子滤波而言，GPF不需要重采样步骤，更容易实现并行化。在系统状态可以用单峰高斯模型描述的情况下，精度不在粒子滤波之下。

由于基本粒子滤波和高斯粒子滤波都是基于蒙特卡罗采样，采样粒子容易在状态空间形成“空隙和簇”，导致算法估计性能下降。虽然通过增加粒子数可以减缓这种现象，但同时增加了算法的复杂度。高斯粒子滤波虽然避免了基本粒子滤波必须的重采样，提高了计算效率，但采用蒙特卡罗采样所带来的估计性能下降问题仍然存在。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于拟蒙特卡罗QMC采样的并行高斯粒子滤波方法，以在保证估计性能的条件下提高运算速度。

实现本发明技术目的的方案，包括以下步骤：

并行拟蒙特卡罗序列产生步骤：采用跳跃拟蒙特卡罗序列产生方法，并行产生后续滤波所需的P条拟随机样本序列；

并行高斯粒子滤波步骤：通过P个执行单元，分别将每条拟随机样本序列转换为服从指定高斯分布的拟高斯样本序列，然后同时对非线性动态系统的状态进行时间更新和状态更新，最后对所有执行单元进行综合处理，完成对非线性动态系统的滤波。

所述的并行产生后续滤波所需的P条拟随机样本序列，包括以下步骤：

(1)产生伪随机数x；

(2)以x为初始值产生一个长度为P的串行拟蒙特卡罗序列sq[p]，p＝1，2，...，P；

(3)分别以sq[1]为初始值产生跳跃拟蒙特卡罗序列Q₁，以sq[2]为初始值产生跳跃拟蒙特卡罗序列Q₂，以sq[p]为初始值产生跳跃拟蒙特卡罗序列Q_P；

(4)将产生的P条拟随机样本序列分别输入到P个执行单元。

所述的将每条拟随机样本序列转换为服从指定高斯分布的拟高斯样本序列，其步骤是：首先采用逆采样将所述的拟随机样本序列转换为服从单位高斯分布的拟高斯样本序列；再通过线性变换，把服从单位高斯分布的拟高斯样本序列转换为服从指定均值和协方差矩阵的拟高斯样本序列。

所述的对非线性动态系统的状态进行时间更新和状态更新，其步骤是：首先，每个执行单元将所述的服从指定高斯分布的拟高斯样本序列代入非线性动态系统状态方程，得到时间更新后的样本；然后根据观测数据，计算每个样本的权值，分别计算每个执行单元内样本的加权均值、加权协方差矩阵与权值和，得到更新后的状态。

所述的对所有执行单元进行综合处理，其步骤是：首先计算所有执行单元中样本的权值和，均值和与协方差矩阵和；然后将均值和与协方差矩阵和除以权值和得到非线性动态系统状态估计结果。

本发明具有以下优点：

1.由于本发明采用了拟蒙特卡罗QMC方法，产生了比蒙特卡罗MC随机样本更均匀的低偏差样本点，因而在使用较少样本点的情况下，提高了高斯粒子滤波的精度。

2.由于本发明采用了并行拟蒙特卡罗QMC序列产生方法，保证了各并行单元样本的统计关系，提高了滤波性能的稳定性。

仿真结果表明，在样本数相对较少时，本发明可以更均匀的逼近待估计状态的后验概率密度，从而比基本粒子滤波和高斯粒子滤波取得更优的滤波性能；而本发明采用并行处理很大的提高了滤波的实时性。

附图说明

图1是本发明的整体流程图；

图2是含有目标的模拟红外灰度图像；

图3是模拟红外图像中目标的运动轨迹；

图4是样本数目为50时本发明与高斯粒子滤波、粒子滤波仿真误差比较图；

图5是样本数目为100时本发明与高斯粒子滤波、粒子滤波仿真误差比较图；

图6是样本数目为400时本发明与高斯粒子滤波、粒子滤波仿真误差比较图；

图7是本发明执行单元个数与并行加速比的关系图。

具体实施方式

通过建立非线性系统动态模型来实施本发明提出的滤波方法，具体模型如下：

状态方程：X_t＝f(X_t-1)+w_t    (1)

观测方程：Z_t＝h(X_t)+e_t    (2)

其中，f(·)，h(·)为有界的非线性函数，X_t为系统在t时刻的状态，Z_t为系统在t时刻的观测值；w_t为过程噪声，e_t为观测噪声。

本发明滤波方法所涉及的参数包括：

N：为采样数目，P＝2^k：为并行单元数目，p＝1，2，...，P：  为并行单元的序号。

参照图1，本发明方法包括并行拟蒙特卡罗序列产生步骤和并行高斯粒子滤波步骤。

一.并行拟蒙特卡罗序列产生步骤

采用跳跃拟蒙特卡罗序列产生方法，并行产生后续滤波所需的P条拟随机样本序列，具体步骤如下：

1)产生伪随机数x；

采用线性同余法，产生一个伪随机数为：

x_t+1＝ax_tmod M，t＝0，1，...    (3)

其中，a和M分别是乘子和模，t表示时间。

2)以x为初始值产生一个长度为P的串行拟蒙特卡罗序列sq[p]，p＝1，2，...，P；

本发明产生的串行拟蒙特卡罗序列是Sobol序列，该序列的产生需要一组“直接数”。在本发明实施之前，该“直接数”已计算完毕并存入固定寄存器，便于实时调用。该“直接数”的计算包括如下步骤：

2a.选择简单多项式：

P_s＝x^d+a₁x^d-1+...+a_d-1x+1，a_i∈{0，1}    (4)

2b.按照给定初始值迭代产生“直接数”v_i：

$v_i=a_1{v}_{i-1}\⊕a_2{v}_{i-2}\⊕...\⊕a_{d-1}{v}_{i-d+1}\⊕v_{i-d}\⊕[v_{i-d}/2^d],i>d---\left(5\right)$

本发明中使用的Sobol序列可以通过当前的Sobol数sq[p]与选定的“直接数”v按照式(6)运算，迭代P次后得到串行拟蒙特卡罗序列sq[p]，p＝1，2，...，P：

$\mathrm{sq}[p+1]=\mathrm{sq}\left[p\right]\⊕v_{z\left[p\right]}---\left(6\right)$

其中“直接数”v_z[p]的下标z[p]是p的二进制形式的最右端0的位数，

表示按位异或运算。

3)分别以sq[1]为初始值产生跳跃拟蒙特卡罗序列Q₁，以sq[2]为初始值产生跳跃拟蒙特卡罗序列Q₂，以sq[p]为初始值产生跳跃拟蒙特卡罗序列Q_P；

以sq[1]为例，跳跃拟蒙特卡罗序列Q₁的产生，包括以下步骤：

3a.利用存储器中的“直接数”生成“跳跃直接数”：

$V_{z[\mathrm{pP}-1]}=v_{z[\mathrm{pP}-1]}\⊕v_{z[\mathrm{pP}-2]}\⊕...\⊕v_{z[(p-1)P-1]}---\left(7\right)$

其中，“直接数”和“跳跃直接数”的下标z[i]是i的二进制形式的最右端0的位数，i＝pP-1，pP-2，…，(p-1)P-1，

表示按位异或运算。

3b.将“跳跃直接数”存储于固定的寄存器中，便于实时调用；

3c.以sq[1]作为初始值，与选定的“跳跃直接数”进行迭代N/P次运算，得到跳跃拟蒙特卡罗序列Q₁：

$Q_{1,i+1}=V_{z[\mathrm{iP}-1]}\⊕Q_{1,i}---\left(8\right)$

其中，V_z[iP-1]是事先计算好并存放在存储器的跳跃直接数，z[iP-1]是iP-1的二进制形式的最右端0的位数，i＝1，2，…，N/P。

按照同样的步骤，分别以sq[2]为初始值产生跳跃拟蒙特卡罗序列Q₂，以sq[p]为初始值产生跳跃拟蒙特卡罗序列Q_P。

4)将产生的P条拟随机样本序列分别输入到P个执行单元。

本发明中并行拟蒙特卡罗序列产生完毕，将得到的P条拟随机样本序列分别提供给后续操作的P个执行单元完成滤波。

二.并行高斯粒子滤波步骤

1)假设t时刻，第p个执行单元将拟随机样本序列转换为服从均值为高斯分布的拟高斯样本序列，具体如下：

1a.将所述的第p条拟随机样本Q_p，t-1通过逆采样转换为服从单位高斯分布的拟高斯样本序列：

Y_p，t-1＝Gau^-1(Q_p，t-1)    (9)

其中，Gau^-1为单位高斯分布函数的逆函数；

1b.将拟高斯样本序列Y_p，t-1线性变换为服从前一时刻均值向量μ_t-1和协方差矩阵∑_t-1的样本序列X_p，t-1：

X_p，t-1＝R_t-1Y_p，t-1+μ_t-1    (10)

其中， ${\mathrm{\Σ}}_{t-1}=R_{t-1}^T{R}_{t-1}.$

2)将服从指定高斯分布的拟蒙特卡罗样本序列X_p，t-1代入非线性系统模型的状态方程式(1)，得到更新后的样本序列：

X_p，t＝f(X_p，t-1)+w_t。    (11)

3)根据当前时刻得到的观测值Z_t，计算每个更新后的样本权值：

W_p，t＝q(Z_t|X_p，t)    (12)

其中q(·)为似然函数，由观测方程式(2)得到。

4)根据得到的新样本序列X_p，t和样本权值W_p，t，分别计算第p个执行单元中N/P个加权样本的和

协方差矩阵

以及权值的和为：

${\mathrm{\μ}}_t^p=\underset{i=1}{\overset{N/P}{\mathrm{\Σ}}}{W}_{p,t}^i{X}_{p,t}^i---\left(13a\right)$

${\mathrm{\Σ}}_t^p=\underset{i=1}{\overset{N/P}{\mathrm{\Σ}}}{W}_{p,t}^i{X}_{p,t}^i{\left(X_{p,t}^i\right)}^T---\left(13b\right)$

${\mathrm{SW}}_t^p=\underset{i=1}{\overset{N/P}{\mathrm{\Σ}}}{W}_{p,t}^i---\left(13c\right)$

该加权样本的和

为更新后的状态。

5)根据得到的每个执行单元中加权样本的和

协方差矩阵

以及权值的和

p＝1，2，…，P，计算系统样本权值总和SW、系统状态均值向量μ_t与协方差矩阵∑_n为：

$\mathrm{SW}=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\mathrm{SW}}_t^p---\left(14a\right)$

${\mathrm{\μ}}_t=\frac{1}{\mathrm{SW}}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\μ}}_t^p---\left(14b\right)$

${\mathrm{\Σ}}_t=\frac{1}{\mathrm{SW}}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\left({\mathrm{\Σ}}_t^p\right)-{\mathrm{\μ}}_t{\left({\mathrm{\μ}}_t\right)}^T---\left(14c\right)$

均值向量μ_t与协方差矩阵∑_t即为t时刻滤波结果。

6)将协方差矩阵∑_t进行Cholesky分解，得到一个上三角矩阵Rt和下三角矩阵

即 ${\mathrm{\Σ}}_t=R_t{R}_t^T,$ R_t将用于t+1时刻拟高斯样本的产生，完成t+1时刻的滤波。

非线性动态系统经过多次滤波后，得到的输出结果将满足性能要求，在新的观测值不断得到的情况下，滤波操作将一直持续。

本发明的效果通过仿真实例进一步说明：

1、仿真条件

选择一个非线性运动目标跟踪的例子，模拟产生了一组包含一个弱小目标的红外灰度图像数据，共有80帧，大小为40×40，并加入仿真的弱小目标，如图2所示。其中，图2(a)是序列中第5帧模拟图像，图2(b)是序列中第25帧模拟图像，目标肉眼无法识别。弱小目标从图像序列的第1帧出现一直到第80帧，在二维平面内作非线性运动，其运动轨迹如图3所示。目标动态方程为：

x_t+1＝Fx_t+w_t    (15)

其中，t是图像帧数， $x_t={\left[\begin{array}{ccccc}{x}_t& {\stackrel{\·}{x}}_t& y_t& {\stackrel{\·}{y}}_t& I_t\end{array}\right]}^T$ 是状态向量；(x_t，y_t)，

和I_t分别表示第t帧目标的位置，速度和强度； $F=\left[\begin{array}{ccc}{F}_s& 0& 0\\ 0& F_s& 0\\ 0& 0& 1\end{array}\right],$ $F_s=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right];$ w_t为协方差矩阵为w的零均值白噪声， $W=\left[\begin{array}{ccc}{Q}_s& 0& 0\\ 0& Q_s& 0\\ 0& 0& q_{2}T\end{array}\right],$ $Q_s=q_1\left[\begin{array}{cc}{T}^3/3& T^3/3\\ T^3/2& T\end{array}\right],$ T为采样周期，q₁和q₂分别表示目标运动噪声和目标强度噪声的功率谱密度。本仿真实例取T＝1，q₁＝0.005，q₂＝0.01，∑＝0.7，目标初始状态给定为x₀＝[70.470.318]^T。

观测模型：

其中表示目标信号对观测的贡献，

为噪声：

$h_t^{(i,j)}\left(x_t\right)\≈\frac{{\mathrm{\Δ}}_x{\mathrm{\Δ}}_y{I}_t}{2\mathrm{\π}{\mathrm{\Σ}}^2}\exp \{-\frac{{({\mathrm{i\Δ}}_x-x_t)}^2+{({\mathrm{i\Δ}}_y-y_t)}^2}{2{\mathrm{\Σ}}^2}\}---\left(17\right)$

$e_t^{(i,j)}~N(0,{\mathrm{\σ}}^2).---\left(18\right)$

2.仿真内容及结果

1)根据仿真条件，以给定初始状态作为本发明滤波方法的初始值，根据事先存储的“直接数”生成滤波所需要的拟蒙特卡罗样本；按照发明中所述的并行高斯滤波方法进行时间更新和观测更新，最后计算得到目标的估计状态。同时，为了验证本发明的优点，将本发明与基本粒子滤波PF和高斯粒子滤波GPF进行性能对比，粒子数目分别为N＝50，100，400。进行100次蒙特卡罗仿真得到三种滤波方法在信噪比为2时的均方根误差比较，如图4、图5和图6所示，其中均方根误差RMSE如下式定义：

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

表1对比了在不同粒子数目情况下三种滤波方法的失跟率，实验分别在粒子数为50，100和400的情况下进行了100次蒙特卡罗仿真。

表1  三种滤波方法的失跟率对比

 

粒子数	PF	GPF	本发明方法
				50	54％	22％	10％
100	47％	10％	4％
				400	8％	3％	1％

从表1可以看出，在粒子数较少时，PF和GPF的失跟率都非常大，而本发明方法的失跟率则很小，而粒子数增加到400时，三者的失跟率都减小，但本发明方法依然优于PF和GPF。

从图4、图5和图6可以看出，在粒子数目较少时，本发明进行估计的RMSE明显小于PF和GPF，估计性能优于PF和GPF，而随着粒子数的增加，三种滤波器的估计性能差别逐渐减小。

通过理论分析和实验结果，可以看出在粒子数相对较少时，粒子的最优分布很重要，在同样数目的粒子情况下，本发明方法可以更均匀的逼近待估计状态的后验概率密度，从而比PF和GPF取得更优的滤波性能。尽管三种滤波算法在粒子数很大时估计精度相近，但本发明方法用较少的粒子可以达到比较好的估计性能，因此在实际应用中会非常有用。

3)为了分析本发明并行处理的优势，根据图1所示方法流程，按照Amdahl定律对本发明的处理速度作了仿真分析。

Amdalh定律如下式：

$S_P=\frac{1}{(s+m/P)}---\left(20\right)$

其中，S_P表示并行处理与串行处理的加速比，s为串行处理工作量，m表示并行处理工作量，P为执行单元个数。

仿真中，串行处理的计算时间≤0.1％，根据Amdalh定律，随着执行单元数目的增加，并行处理的加速比变化为图7虚线所示。实际仿真中，在处理器个数分别为2、4、8、16、32、64的情况下得到加速比如图7的实线所示，比串行处理占0.1％情况下的理论估计值还要好，可见发明中串行处理的计算时间比估计的0.1％小，因此发明采用并行系统处理能获得很好的效率。

Claims (5)

1.一种基于拟蒙特卡罗采样的并行高斯粒子滤波方法，包括：

并行拟蒙特卡罗序列产生步骤：采用跳跃拟蒙特卡罗序列产生方法，并行产生后续滤波所需的P条拟随机样本序列；

并行高斯粒子滤波步骤：通过P个执行单元，分别将每条拟随机样本序列转换为服从指定高斯分布的拟高斯样本序列，然后同时对非线性动态系统的状态进行时间更新和状态更新，最后对所有执行单元进行综合处理，完成对非线性动态系统的滤波。

2.根据权利1所述的并行高斯粒子滤波方法，其中所述的并行产生后续滤波所需的P条拟随机样本序列，包括以下步骤；

(1)产生伪随机数x；

(2)以x为初始值产生一个长度为P的串行拟蒙特卡罗序列sq[p]，p＝1，2，...，P；

(3)分别以sq[1]为初始值产生跳跃拟蒙特卡罗序列Q1，以sq[2]为初始值产生跳跃拟蒙特卡罗序列Q₂，以sq[p]为初始值产生跳跃拟蒙特卡罗序列Q_P；

(4)将产生的P条拟随机样本序列分别输入到P个执行单元。

3.根据权利1所述的并行高斯粒子滤波方法，其中所述的将每条拟随机样本序列转换为服从指定高斯分布的拟高斯样本序列，首先采用逆采样将所述的拟随机样本序列转换为服从单位高斯分布的拟高斯样本序列；再通过线性变换，把服从单位高斯分布的拟高斯样本序列转换为服从指定均值和协方差矩阵的拟高斯样本序列。

4.根据权利1或3所述的并行高斯粒子滤波方法，其中所述的对非线性动态系统的状态进行时间更新和状态更新，包括以下步骤：

(1)每个执行单元将所述的服从指定高斯分布的拟高斯样本序列代入非线性动态系统状态方程，得到时间更新后的样本；

(2)根据观测数据，计算每个样本的权值，分别计算每个执行单元内样本的加权均值、加权协方差矩阵与权值和，得到更新后的状态。

5.根据权利1所述的并行高斯粒子滤波方法，其中所述的对所有执行单元进行综合处理，包括计算所有执行单元中样本的权值和，均值和与协方差矩阵和，并将均值和与协方差矩阵和除以权值和得到非线性动态系统状态估计结果。

Images