基于Sigma点<italic>H</italic>∞滤波的拟蒙特卡罗粒子滤波算法

引用本文

孔云波, 冯新喜, 鹿传国, 刘振涛. 基于Sigma点H_∞滤波的拟蒙特卡罗粒子滤波算法 . 44(6): 1831-1837
KONG Yun-bo, FENG Xin-xi, LU Chuan-guo, LIU Zhen-tao. Quasi-Monte Carlo particle filter algorithm based on sigma pointH_∞ . Journal of Jilin University Engineering and Technology Edition, 44(6): 1831-1837 复制到剪切板

Permissions

基于Sigma点H_∞滤波的拟蒙特卡罗粒子滤波算法

孔云波¹, 冯新喜¹, 鹿传国¹, 刘振涛²

1.空军工程大学信息与导航学院,西安 710077

2.中国人民解放军 95072部队,广西南宁 530000

作者简介:孔云波(1987-),男,博士研究生.研究方向:多传感器信息融合.E-mail:kongyunbo123@163.com

基金:陕西省自然科学基金项目(2012JM8023);

摘要

在滤波算法中,用Sigma点H_∞滤波来产生重要性概率密度函数,由于Sigma点H_∞滤波对不确定观测噪声具有较强的鲁棒性,而且在滤波过程中考虑了最新的观测值,因此由其产生的重要性函数更逼近于真实的后验概率分布。同时在重采样阶段,利用拟蒙特卡罗重采样算法进行重采样,有效地克服了粒子退化现象并提高了状态估计精度。仿真结果表明了所提算法的可行性和有效性。

关键词: 通信技术; Sigma点转换; H∞滤波; 准蒙特卡罗; 粒子滤波; 非线性系统

中图分类号:TN953 文献标志码:A 文章编号:1671-5497(2014)06-1831-07

Quasi-Monte Carlo particle filter algorithm based on sigma pointH_∞

KONG Yun-bo¹, FENG Xin-xi¹, LU Chuan-guo¹, LIU Zhen-tao²

1. Information and Navigation Institute, Airforce Engineering University, Xi'an 710077,China

2. Unit of 95072 of the PLA, Nanning 530000,China

Abstract

A new filtering algorithm is proposed. In this algorithm, the probability density function is generated by the sigma pointH_∞ filter. The sigma pointH_∞ filter has very high accuracy and strong robustness to uncertain observation noise, and the filter algorithm integrates the new observation. So the generated probability density function can reasonably approximate the real posterior probability distribution of the system state. At the resampling step, the degeneration problem is effectively overcome and the accuracy of state estimation is improved by using the quasi-Monte Carlo-based resampling algorithm. Simulation results demonstrated the feasibility and effectiveness of the proposed algorithm.

Keyword: communication; Sigma point transformation; H∞ filtering; quasi-Monte Carlo; particle filtering; nonlinear system

Show Figures

0 引言

粒子滤波是一种基于贝叶斯估计的非线性滤波算法,目前已被广泛应用于实现非线性、非高斯环境下的参数估计和状态估计。粒子滤波借助于一个容易采样的重要性密度函数得到满足后验概率分布的样本粒子,重要性密度函数选择的优劣直接决定了粒子滤波的性能。文献[1]最早提出了使用先验转移密度作为提议分布的Bootstrap粒子滤波算法。该算法简单、易于实现,但由于没有考虑最新的观测信息常导致粒子集偏离真实的后验概率。文献[2]提出了利用扩展卡尔曼滤波(Extended kalman filter,EKF)来产生重要密度函数。文献[3]提出了用无味卡尔曼滤波(Unscented kalman filter,UKF)来产生重要性函数,使滤波器的性能又进一步提高。文献[4]提出了一种基于扩展 H_∞粒子滤波算法,该算法利用扩展 H_∞滤波(Extended H_∞ filter,EHF)产生重要性函数,使得在进一步提高滤波器性能的同时降低了算法的复杂度。此外,常用的建议分布有转移先验建议分布、混合建议分布、退火先验建议分布、似然函数建议分布、基于梯度的转移密度概率建议分布等^{[ 5]}。这些方法在一定程度上避免了粒子退化现象,提高了粒子滤波算法的估计精度。

上述基于重要性采样的粒子滤波都存在一个主要问题是粒子退化现象。抑制退化现象的最主要方法是增加粒子数和重采样。重采样的目的是去掉权值较小的粒子数目,保留权值大的粒子。但经过重采样,权值大的粒子被多次复制,丧失了粒子的多样性。对此,文献[6]提出了在重采样后增加马尔科夫蒙特卡罗移动步骤的方法,该方法可以减弱粒子间的相关性,抑制样本的贫化,但不能保证其收敛性。文献[7]提出了使用高斯分布近似状态后验概率密度的高斯粒子滤波。文献[8]提出了一种结合准蒙特卡罗方法的粒子滤波算法,在重要性采样后,将生成的随机化蒙特卡罗序列分别映射到以权值大的粒子为核心的独立空间上,避免了直接对采样空间进行预测,同时又保证了样本的多样性。

本文提出了一种新的基于Sigma点 H_∞滤波的拟蒙特卡罗粒子滤波算法,在重要性采样阶段,利用Sigma点 H_∞滤波来产生重要性密度函数,然后利用拟蒙特卡罗重采样算法进行重采样。仿真试验表明,该算法在提高估计精度的同时,可以有效地抑制样本枯竭现象。

1 Sigma点 H_∞滤波算法

1.1 H_∞滤波方法

H_∞滤波方法是将鲁棒控制中的 H_∞范数应用于滤波,以解决系统噪声的统计特性不确定问题。该方法采用极大极小准则,无需知道噪声的统计特性等先验信息,在最大噪声作用下寻求干扰引起的误差范数最小。因此,该方法具有较强的鲁棒性。

考虑离散的线性系统:

$\begin{matrix} \begin{matrix} X_{k + 1} = F_{k} X_{k} + B_{k} W_{k} (1) \\ Z_{k} = H_{k} X_{k} + V_{k} (2) \end{matrix} \end{matrix}$

式中: X_k为系统的状态向量; F_k为状态转移矩阵, B_k为系统噪声矩阵; H_k为量测矩阵; W_k和 V_k为互不相关的白噪声,其方差为 Q_k和 R_k。

若系统过程噪声和测量噪声 W_k和 V_k满足 $\begin{matrix} \overset{N}{\sum_{k = 0}} \end{matrix} \begin{matrix} ‖ W_{k} ‖^{2} \end{matrix}$ <∞, $\begin{matrix} \overset{N}{\sum_{k = 0}} \end{matrix} \begin{matrix} ‖ V_{k} ‖^{2} \end{matrix}$ <∞。如果定义代价函数:

$\begin{matrix} J = \frac{\overset{N}{\sum_{k = 0}} ‖ X_{k} - {\hat{X}}_{k} ‖^{2}}{‖ X_{0} - {\hat{X}}_{0} ‖_{P {(0)}^{- 1}}^{2} + \overset{N}{\sum_{k = 0}} (‖ W_{k} ‖^{2} + ‖ V_{k} ‖^{2})} (3) \end{matrix}$

则 H_∞滤波问题使系统的输出 X_k在各种可能的条件下满足:

$\begin{matrix} X_{k} = \underset{x \in R^{n_{x}}}{argmin} ‖ J ‖_{x} (4) \end{matrix}$

式中:

$\begin{matrix} ‖ X_{0} - {\hat{X}}_{0} ‖_{P {(0)}^{- 1}}^{2} \end{matrix}$ 􀰛 $\begin{matrix} (X_{0} - {\hat{X}}_{0})^{T} \end{matrix}$ P(0) ^-¹( X₀ - $\begin{matrix} {\hat{X}}_{0} \end{matrix}$ )。

由于最优 H_∞滤波器在一般情况下难以求得解析形式的解。因此,最优 H_∞滤波问题实质是一种次优滤波算法,设定一门限,使得代价函数满足:

$\begin{matrix} ‖ J ‖_{x} = \sup_{x \in R^{n_{x}}} J \leq γ^{2} (5) \end{matrix}$

在式(5)所示的约束条件下,可得 H_∞滤波方程,即:

$\begin{matrix} \{\begin{matrix} X_{k | k - 1} = F_{k} X_{k - 1} \\ S_{k} = (1 - P_{k - 1} / γ^{2} + {H^{T}}_{k} H_{k} P_{k - 1})^{- 1} \\ K_{k} = F_{k} P_{k - 1} S_{k} {F^{T}}_{k} \\ X_{k} = X_{k | k - 1} + K_{k} (Z_{k} - H_{k} X_{k | k - 1}) \\ P_{k} = F_{k} P_{k - 1} S_{k} {F^{T}}_{k} + W_{k + 1} \end{matrix} (6) \end{matrix}$

从递推表达式可以看出, H_∞滤波器的滤波性能与 γ的取值具有较大的关系,当 γ取最小值时, H_∞滤波器具有最好的鲁棒性,但此时方差估计不一定达到最小;而当 γ→ ¥时,方差估计达到最小,但滤波器的鲁棒性较差。所以, H_∞滤波器可以理解为二范数和无穷范数间的折中,通过选取适当的 γ,可以同时具有鲁棒性和估计方差最小的特性。

1.2 Sigma点 H_∞滤波方法

类似于UKF滤波器,文献[9]将 H_∞滤波应用到非线性估计领域,提出了Sigma点 H_∞滤波方法。Sigma点转换根据加权统计线性回归(WSLR)的思想,通过真正的非线性函数,计算随机变量经过非线性变换后的均值和方差。由于加权统计线性回归考虑了随机变量的先验统计特性,因此能获得比截断泰勒级数更小的线性化误差。

设 n维随机向量 x,其统计特性为( $\begin{matrix} \bar{x} \end{matrix}$ , P_x), x通过非线性函数 h映射后产生随机变量 z,即 z=h( x)。设 z的统计特性为( $\begin{matrix} \bar{x} \end{matrix}$ , P_x),协方差为 P_xz= $\begin{matrix} E [(x - \bar{x}) {(z - \bar{z})}^{T} \end{matrix}$ ]。则基于随机向量 x,利用Sigma点转换可近似求得 $\begin{matrix} \bar{z} \end{matrix}$ 、 P_z和 P_xz的值。

为了将Sigma点转换方法应用于 H_∞滤波器,采用式(7)进行状态估计协方差阵递推^{[ 9]}:

$\begin{matrix} \begin{matrix} P_{k} = P_{k | k - 1} - PP (7) \\ 式中 : \\ PP = [\begin{matrix} P_{xz} & P_{k | k - 1} \end{matrix}] [\begin{matrix} P_{zz} + I & {P^{T}}_{xz} \\ P_{xz} & P_{k | k - 1} - γ^{2} I \end{matrix}] [\begin{matrix} {P^{T}}_{xz} \\ {P^{T}}_{k | k - 1} \end{matrix}] \end{matrix} \end{matrix}$

将Sigma点转换方法应用于 H_∞滤波器,就可以得到基于Sigma点的 H_∞非线性滤波方法。

2 本文算法

2.1 拟蒙特卡罗重采样算法

拟蒙特卡罗^{[ 10]}(Quasi-Montor Carlo,QMC)的主要思想是用更加规则分布的确定性超均匀分布序列来构造蒙特卡罗(Montor Carlo,MC)方法中用随机序列所形成的近似。这种确定性方法可以选择采样空间最优的可能性分布点,避免发生在MC随机采样过程中粒子之间的间隙过大、粒子层叠等现象,提高了采样效率。为克服粒子退化问题,本文采用拟蒙特卡罗重采样(Quasi-Montor Carlo resample,QMCR)算法,具体参见文献[8]。

2.2 Sigma点 H_∞拟蒙特卡罗粒子滤波算法

利用Sigma点 H_∞滤波产生重要性函数,利用拟蒙特卡罗算法进行重采样,就得到了Sigma点 H_∞拟蒙特卡罗粒子滤波算法(Quasi-Monte Carlo particle filter based on Sigma point H_∞,SHQPF)。由于Sigma点 H_∞滤波对不确定观测噪声具有较强的鲁棒性,而且在滤波过程中考虑了最新的观测值,所以根据该重要性函数产生的样本更接近于真实的采样样本。基于Sigma点 H_∞拟蒙特卡罗粒子滤波算法的步骤如下:

步骤1 初始化。

k=0,从先验分布中产生样本:

$\begin{matrix} \{x_{0}^{i}, i = 1,2, \dots, N_{s}\} 。 \end{matrix}$

步骤2 重要性采样。

i=1,2,…, N_s,在每一时刻用Sigma点 H_∞滤波器算法更新粒子。

$\begin{matrix} \begin{matrix} 选取粒子 : x_{k - 1}^{ia} = \\ [\begin{matrix} {\bar{x}}_{k - 1}^{a} & {\bar{x}}_{k - 1}^{a} + ζ \sqrt[]{P_{k - 1}^{ia}} & {\bar{x}}_{k - 1}^{a} - ζ \sqrt[]{P_{k - 1}^{ia}} \end{matrix}] \\ 时间更新 : x_{k | k - 1}^{ix} = f (x_{k - 1}^{ix}) (8) \\ {\bar{x}}_{k | k - 1}^{i} = \overset{2 n_{a}}{\sum_{j = 0}} W_{j}^{(m)} x_{j, k | k - 1}^{ix} (9) \\ P_{k | k - 1}^{i} = \overset{2 n_{a}}{\sum_{j = 0}} W_{j}^{(c)} [x_{j, k | k - 1}^{ix} - x_{k | k - 1}^{i}] \\ [x_{j, k | k - 1}^{ix} - x_{k | k - 1}^{i}]^{T} + Q_{k - 1}^{i} (10) \\ z_{k | k - 1}^{i} = h (x_{k | k - 1}^{ix}) (11) \\ {\bar{z}}_{k | k - 1}^{i} = \overset{2 n_{a}}{\sum_{j = 0}} W_{j}^{(m)} z_{j, k | k - 1}^{i} (12) \\ 预测更新 : P_{zz} = \\ \overset{2 n_{a}}{\sum_{j = 0}} W_{j}^{(c)} [z_{j, k | k - 1}^{i} - {\bar{z}}_{k | k - 1}^{i}] [z_{j, k | k - 1}^{i} - {\bar{z}}_{k | k - 1}^{i}]^{T} (13) \\ P_{xz} = \\ \overset{2 n_{a}}{\sum_{j = 0}} W_{j}^{(c)} [x_{j, k | k - 1}^{ix} - {\bar{x}}_{k | k - 1}^{i}] [z_{j, k | k - 1}^{i} - {\bar{z}}_{k | k - 1}^{i}]^{T} (14) \\ K_{k} = P_{xz} P_{zz}^{- 1} (15) \\ {\bar{x}}_{k}^{i} = {\bar{x}}_{k | k - 1}^{i} + K_{k} (z_{k} - {\bar{z}}_{k | k - 1}^{i}) (16) \\ P_{k} = P_{k | k - 1} - PP (17) \\ PP = \\ [\begin{matrix} P_{xz} & P_{k | k - 1} \end{matrix}] [\begin{matrix} P_{zz} + I & {P^{T}}_{xz} \\ P_{xz} & P_{k | k - 1} - γ^{2} I \end{matrix}] [\begin{matrix} {P^{T}}_{xz} \\ {P^{T}}_{k | k - 1} \end{matrix}] (18) \end{matrix} \end{matrix}$

采样粒子:

$\begin{matrix} {\hat{x}}_{k}^{i} ~ q ({\bar{x}}_{k}^{i} | {\bar{x}}_{0 : k - 1}^{i}, z_{1 : k}) = N ({\bar{x}}_{k}^{i}, {\hat{P}}_{k}^{i}) (19) \end{matrix}$

步骤3 计算各粒子权重。

$\begin{matrix} w_{k}^{i} = w_{k - 1}^{i} \frac{p (z_{k} | {\hat{x}}_{k}^{i}) p ({\hat{x}}_{k}^{i} | x_{k - 1}^{i})}{q ({\hat{x}}_{k}^{i} | x_{0 : k - 1}^{i}, z_{1 : k})} (20) \end{matrix}$

步骤4 重采样。

计算 N_eff,并判断 N_eff <N_T是否成立,其中 N_T表示门限值。若成立,则调用QMCR算法进行重采样,得到新的粒子集合{ $\begin{matrix} x_{0 : k}^{i *} \end{matrix}$ , i=1,2,…, N_s}。

步骤5 状态估计。

$\begin{matrix} \begin{matrix} X_{k} = \overset{N_{s}}{\sum_{i = 1}} w_{k}^{i} {\bar{x}}_{k}^{i} (21) \\ P_{k} = \overset{N_{s}}{\sum_{i = 1}} w_{k}^{i} [{\bar{x}}_{k}^{i} - X_{k}] [{\bar{x}}_{k}^{i} - X_{k}]^{T} (22) \end{matrix} \end{matrix}$

3 仿真试验及结果分析

为了验证本文算法的性能,分4个试验讨论算法在不同仿真条件下的参数变化情况和滤波性能。试验一讨论了参数 γ对算法滤波性能的影响;试验二将QMCR算法和残差重采样算法、多项式重采样算法、系统重采样算法、分层重采样算法进行比较;试验三对所提算法与扩展卡尔曼粒子滤波(Extended kalman particle filter,EKPF)算法、无迹卡尔曼粒子滤波(Unscented particle filter,UPF)算法的估计误差、收敛情况和跟踪稳定性等性能进行比较;试验四对所提算法的有界收敛界进行了讨论。

3.1 试验一

过程模型:

$\begin{matrix} x (t) = 0.5 x (t - 1) + \sin (0.04 πt) + 1 + w (t) (23) \end{matrix}$

量测模型:

$\begin{matrix} z (t) = \{\begin{matrix} \frac{x {(t)}^{2}}{50} + v (t), t \leq 30 \\ \frac{x {(t)}^{2}}{20} - 2 + v (t), t > 30 \end{matrix} (24) \end{matrix}$

式中: w( t) ~Γ(3,2); v( t) ~N(0,1)。

初始状态 x( t) =1,仿真时长 T=50,粒子数为300。分别计算算法在参数 γ取不同值时估计值与真值之间的均方根误差( RMSE),仿真次数为100次。

$\begin{matrix} RMSE (\hat{x}) = {\{\frac{1}{T} \overset{T}{\sum_{t = 1}} {[x (t) - \hat{x} (t)]}^{2}\}}^{\frac{1}{2}} (25) \end{matrix}$

由1.1节理论分析可知,在本文算法中,参数 γ的取值对滤波性能影响较大。随着 γ的逐渐增大, H_∞滤波器的状态估计会逼近最优估计值,但此时的鲁棒性变差。随着 γ的逐渐减小, H_∞滤波器鲁棒性逐渐增加,但状态估计误差逐渐增大。一般情况下,通过使给定 γ的不断减小,使其逐步逼近Sigma点 H_∞拟蒙特卡罗粒子滤波估计问题对应的 γ₀时,就可以得到Sigma点 H_∞拟蒙特卡罗粒子滤波的最优估计,其中前提条件是 γ>γ₀,否则Sigma点 H_∞拟蒙特卡罗粒子滤波将会无解。图1给出了100次蒙特卡罗仿真下,算法滤波误差随参数 γ的变化情况。从图1中可以看出,当参数 γ值从10逐渐减小到2时,算法滤波的均方根误差逐渐减小,当参数的值从2减小到0.1时,均方根误差逐渐增大,因此,可以推断最优解大概在 γ=2 .0附近。

	Figure Option View Download New Window
	图1 SHQPF算法在 γ取不同值时均方根误差曲线图Fig.1 RMSE of SHQPF in the different γ

3.2 试验二

选取具有典型的非线性特征的单变量非静态增长模型(UNGM),其过程模型和量测模型如下:

过程模型:

$\begin{matrix} \begin{matrix} x (t) = 0.5 x (t - 1) + \frac{25 x (t - 1)}{1 + [x {(t - 1)]}^{2}} + \\ 8 \cos [1.2 (t - 1)] + w (t) (26) \end{matrix} \end{matrix}$

量测模型:

$\begin{matrix} z (t) = \frac{x {(t)}^{2}}{20} + v (t) (27) \end{matrix}$

式中: w( t)和 v( t)为零均值高斯噪声。

该系统的特征是高度非线性,似然函数呈双峰。

重要性重采样的目的是去掉权值较小的粒子数目,保留权值大的粒子。重要性重采样不依赖于系统模型和具体的应用,仅与重采样算法本身有关^{[ 11]}。本试验将拟蒙特卡罗重采样算法与系统重采样算法、残差重采样算法、多项式重采样算法和分层重采样算法进行比较,重点研究了在采样粒子取不同值时,5种重采样算法的优劣。

图2给出了5种算法的状态估计曲线。由图2可以看出,采用5种重采样算法得到的均方根数据较为接近,为了更清楚地显示他们的差别,采用文献[11]的方法,对均方根数据作如下标准化处理:

$\begin{matrix} d_{i}^{j} = \frac{d_{i}^{j} - Mi n_{i}}{Ma x_{i} - Mi n_{i}} (28) \end{matrix}$

式中: i=1,…, M; j=1,…, N; M、 N分别表示行数和列数; $\begin{matrix} d_{i}^{j} \end{matrix}$ 表示第 i行,第 j列数据; Min_i、 Max_i分别表示第 i行数据的最小值和最大值。

	Figure Option View Download New Window
	图2 状态估计曲线Fig.2 Line of state estimation

图3和表1给出了采样粒子数分别取200、300、400、500时的标准化均方根误差图和平均有效样本数。从图3可以看出,多项式重采样算法在粒子数 N为200、300、400和500时出现峰值较多,效果极差;残差重采样算法在粒子数 N为200和300时峰值较多,在 N为400和500时存在峰值,但已明显减少,因此,残差重采样算法综合表现一般;分层重采样算法在粒子数 N为200、300、400和500时只出现少数峰值,结果较好,优于残差重采样;系统重采样与分层重采样效果较为相似,但整体优于分层重采样;拟蒙特卡罗重采样算法随着采样粒子数的增多,峰值数逐渐减少,在粒子数 N为400时峰值为1个,在 N为500时峰值为0个。另外,从图3中可以看出,拟蒙特卡罗算法的平均值远小于其他4种重采样算法,从整体上看,拟蒙特卡罗重采样算法得到的结果是最优的。

在粒子滤波方法中,有效采样粒子数(Effective sampling size,ESS)是衡量粒子退化程度的一种参数。表1给出了5种算法在不同粒子数下的平均有效样本数。从表1可以看出,随着粒子数 N的增加,各重采样算法的平均有效样本数逐渐增加,但当粒子数相同时,拟蒙特卡罗算法平均有效样本数大于其他4种算法。有效样本数越少,说明粒子退化现象越严重。由此可见,与残差重采样算法、多项式重采样算法、系统重采样算法、分层重采样算法相比,拟蒙特卡罗重采样算法更能抑制粒子退化现象。

	Figure Option View Download New Window
	图3 标准化均方根误差曲线Fig.3 Line of standard REMS

表1 平均有效样本数 Table 1 Mean number of effective sample

3.3 试验三

纯方位跟踪模型(Bearing-only tracking,BOT)是一种具有强非线性的目标跟踪模型。本试验考虑二维平面上的一个点目标的运动情况,假设传感器仅能获得目标的方位信息。传感器的扫描周期 T=1 s,方位角标准差为0.03 rad,跟踪时间为300 s。目标做匀速直线运动,其初始位置为(0 m,600 m),速度为2.57 m/s,航向为5π/6;为了保证系统的可观测性,传感器必须满足一定的机动条件。设置传感器系统初始位置为(0 m,0 m),初始速度为1.542 m/s,初始航向为π/6,在1~300 s做匀速转弯运动,角加速度为1/300π rad/s²。目标与传感器的运动轨迹如图4所示。

	Figure Option View Download New Window
	图4 目标和传感器的运动轨迹Fig.4 Trace of target and sensor

表2给出了各滤波算法的滤波精度和滤波时间随采样粒子数的变化情况。从表2可以看出,所有滤波方法的估计误差都会随着粒子数目的增加而逐渐减少。当粒子数相同时,EKPF、UPF、SHQPF的滤波误差明显减小,但运算所需时间依次增加,即精度的提高是以牺牲运算时间为代价的。与EKPF和UPF相比,SHQPF均方根误差明显较低,说明SHQPF明显提高了目标跟踪的精度,且运行时间与UPF相当。同时在100次仿真试验中,UPF和SHQPF具有较好的收敛性,没有出现跟踪发散现象,而EKPF出现了7次跟踪发散现象。

表2 各算法跟踪性能比较 Table 2 Estimation performance comparison of all kinds filter

3.4 试验四

在粒子滤波算法中,由于各粒子间的交互作用使得统计独立的假设不再成立,所以粒子滤波算法收敛性的分析非常复杂。文献[12]在大量工程实践的基础上提出了粒子滤波算法的有限收敛界(Limited convergence bound,LCB)的概念。因此本文使用有限收敛界进一步分析算法的收敛问题。

选取试验二给出的单变量非静态增长模型(UNGM),该系统是一个复杂的非线性递推Bayesian滤波问题,其特征是高度非线性和其最优解不存在。采用SHQPF对非线性系统进行状态估计,并求出该算法的有限收敛界 N_LCB。该试验以粒子数50为起点,粒子增量为50,波动误差为0.01,蒙特卡洛仿真次数为100。使用均方根误差( RMSE)的时间平均作为衡量SHQPF估计误差的指标,令 N_T为用于时间平均的样点数; $\begin{matrix} {\tilde{x}}_{k} \end{matrix}$ =x_k- $\begin{matrix} {\hat{x}}_{k} \end{matrix}$ 为估计残差,则其估计式为:

$\begin{matrix} \bar{RMSE} ({\hat{x}}_{k}) = \frac{1}{N_{T}} \overset{N_{T}}{\sum_{k = 1}} {[\frac{1}{N_{MC}} \overset{N_{MC}}{\sum_{i = 1}} ({\tilde{x}}_{k}^{i})^{2}]}^{\frac{1}{2}} (29) \end{matrix}$

依据有限收敛界的定义可以计算得 N_LCB=200,即所提算法有界收敛。另外,从图5可以看出,随着粒子数的增大,平均 RMSE值逐渐降低,当粒子数小于200时,平均 RMSE值降低明显,粒子数的增加使得滤波精度明显提高,超过200后,滤波精度增加趋于平缓,说明SHQPF随着粒子数的增大而逐渐收敛,但这种变化趋势并非是均匀的。在实际情况中,当粒子数超过某个门限后,在粒子数增加有限的情况下,对滤波估计性能的影响已经很弱,这与一般粒子滤波是一致的。

	Figure Option View Download New Window
	图5 平均 RMSE值和粒子数Fig.5 Mean RMSE and particle number

4 结束语

提出了一种新的基于Sigma点 H_∞滤波的拟蒙特卡罗粒子滤波算法。算法利用Sigma点 H_∞滤波产生重要性密度函数,由于Sigma点 H_∞滤波对不确定观测噪声具有较强的鲁棒性,而且在滤波过程中考虑了最新的观测值,因此由其产生的重要性函数更逼近于真实的后验概率分布。同时,该算法无需知道噪声的统计特性等先验信息,摆脱了应用环境的束缚。另外,在重采样阶段引入拟蒙特卡罗重采样又克服了传统粒子滤波算法粒子退化现象,提高了滤波器的估计能力。仿真结果表明,本文算法是有效且可行的。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Gordon N J, Salmond D J. Novel approach to nonlinear/non-Gaussian Bayesian state estimation[J]. IEE Proceedings Part F, Radar, Sonar and Navigation, 1993, 140(2): 107-113. [本文引用:1]
[2]	Doucet A, Godsill S J, Andrieu C. On sequential Monte Carlo sampling methods for Bayesian filtering[J]. Statistics and Computing, 2000, 10(3): 197-208. [本文引用:1] [JCR: 1.977]
[3]	Julier S J, Uhlmann J K. Unscented filtering and nonlinear estimation[J]. IEEE Transactions on Signal Processing, 2004, 92(3): 401-422. [本文引用:1] [JCR: 2.813]
[4]	万洋, 王首勇, 于兴伟. 一种扩展H_∞粒子滤波方法[J]. 信号处理, 2010, 26(6): 869-874. Wan Yang, Wang Shou-yong, Yu Xing-wei. A extended H_∞ particle filter algorithm[J]. Signal Processing, 2010, 26(6): 869-874. [本文引用:1] [JCR: 1.851]
[5]	朱志宇. 粒子滤波算法及其应用[M]. 北京: 科学出版社, 2010. [本文引用:1]
[6]	Andrieu C, Freitas J F G, Doucet A. Sequential MCMC for Bayesian model selection[C]∥IEEE Higher Order Statistics Workshop, Ceasarea, 1999: 130-134. [本文引用:1]
[7]	Kotecha J H, Djuric P M. Gaussian particle filtering[J]. IEEE Transactions on Signal Processing, 2003, 51(10): 2592-2601. [本文引用:1] [JCR: 2.813]
[8]	赵玲玲, 马培军, 苏小红. 一种快速准蒙特卡罗粒子滤波算法[J]. 自动化学报, 2010, 36(9): 1351-1356. Zhao Ling-ling, Ma Pei-jun, Su Xiao-hong. A fast quasi-Monte Carlo-based particle filter algorithm[J]. Acta Automatica Sinica, 2010, 36(9): 1351-1356. [本文引用:1] [CJCR: 0.572]
[9]	侯代文, 殷福亮, 陈喆. 基于sigma点H_∞滤波的说话人跟踪方法[J]. 信号处理, 2009, 25(3): 374-378. Hou Dai-wen, Yin Fu-liang, Chen Zhe. Sigma point H_∞ filtering method for speaker tracking[J]. Signal Processing, 2009, 25(3): 374-378. [本文引用:1] [JCR: 1.851]
[10]	Guo D, Wang X D. Quasi-Monte Carlo filtering in nonlinear dynamic systems[J]. IEEE Transactions on Signal Processing, 2006, 54(6): 2087-2098. [本文引用:1] [JCR: 2.813]
[11]	吴宝成. 粒子滤波重采样算法研究及其应用[D]. 哈尔滨: 哈尔滨工业大学计算机科学与技术学院, 2006. Wu Bao-cheng. Research and application of particle filter resampling algorithms[D]. Harbin: School of Computer Science and Technology, Harbin Institute of Technology, 2006. [本文引用:1]
[12]	程水英, 张剑云. 裂变自举例子滤波[J]. 电子学报, 2008, 36(3): 500-504. Cheng Shui-ying, Zhang Jian-yun. Fission bootstrap particle filtering[J]. Acta Electronica Sinice, 2008, 36(3): 500-504. [本文引用:1]

1993

0.0

2000

1.977

0.0

. 2000, 10(3):197-208 DOI:10.1023/A:1008935410038

On sequential Monte Carlo sampling methods for Bayesian filtering

1.University of Cambridge 2.University of Cambridge 3.University of Cambridge

In this article, we present an overview of methods for sequential simulation from posterior distributions. These methods are of particular interest in Bayesian filtering for discrete time dynamic models that are typically nonlinear and non-Gaussian. A general importance sampling framework is developed that unifies many of the methods which have been proposed over the last few decades in several different scientific disciplines. Novel extensions to the existing methods are also proposed. We show in particular how to incorporate local linearisation methods similar to those which have previously been employed in the deterministic filtering literature; these lead to very effective importance distributions. Furthermore we describe a method which uses Rao-Blackwellisation in order to take advantage of the analytic structure present in some important classes of state-space models. In a final section we develop algorithms for prediction, smoothing and evaluation of the likelihood in dynamic models.

2004

2.813

0.0

2010

1.851

0.0

2010

0.0

... 此外,常用的建议分布有转移先验建议分布、混合建议分布、退火先验建议分布、似然函数建议分布、基于梯度的转移密度概率建议分布等^[5] ...

1999

0.0

2003

2.813

0.0

2010

0.0

0.572

. 2010, 36(9):1351-1356 DOI:10.3724/SP.J.1004.2010.01351

A fast quasi-Monte Carlo-based particle filter algorithm

1. School of Computer Science and Technology, Harbin Institute of Technology, Harbin 150001

Particle filters have a high computational complexity when using quasi-Monte Carlo (QMC) methods, and are subject to the sample impoverishment caused by resampling step. To solve these problems, a new particle filter algorithm based on QMC method is proposed. It generates the randomized QMC points after sample importance, and then transforms them into some independent sub-spaces, whose kernels are the particles with heavy weights, to avoid predicting the sampling space and preserve the diversity of samples. The simulation results suggest that the algorithm can escape successfully from the sample impoverishment, provide more accurate estimators than the Monte Carlo (MC) method, meanwhile has a computational cost similar to the general particle filter.

针对准蒙特卡罗(Quasi-Monte Carlo, QMC)方法应用于粒子滤波采样时计算复杂度高, 以及粒子滤波中重采样步骤引起样本枯竭的问题, 提出一种结合准蒙特卡罗方法的粒子滤波算法, 在重要性采样后, 将生成的随机化QMC序列分别映射到以大权重粒子为核心的独立子空间上, 避免了直接对采样空间进行预测, 同时又保持了样本多样性. 实验结果表明该方法可以有效抑制样本枯竭现象, 获得了高于蒙特卡罗(Monte Carlo, MC)方法的估计精度, 而计算效率与粒子滤波相近.

2009

1.851

0.0

... 滤波器,采用式(7)进行状态估计协方差阵递推^[9]: ...

2006

2.813

0.0

... 1 拟蒙特卡罗重采样算法拟蒙特卡罗^[10](Quasi-Montor Carlo,QMC)的主要思想是用更加规则分布的确定性超均匀分布序列来构造蒙特卡罗(Montor Carlo,MC)方法中用随机序列所形成的近似 ...

2006

0.0

... 重要性重采样不依赖于系统模型和具体的应用,仅与重采样算法本身有关^[11] ...

2008

0.0

. 2008, 36(3):500-504

Fission bootstrap particle filtering

Electronic and Engineering Institute,Hefei,Anhui 230037,China

Among the various particle filtering (PF) evolutions,the bootstrap PF (BPF) is the most classical and popular algorithm.However,it is subject to severe sample impoverishment after resampling.To overcome the above problem,the fission BPF (FBPF) algorithm is proposed,in which the preprocess including weights sorting,particle reproducing by fission,and weights normalizing is inserted before the original resampling step as soon as the importance weights degenerate severely.The results of Monte Carlo simulations about a typical severely nonlinear filtering problem with bimodal posterior density have demonstrated that compared to the BPF,the more robust FBPF can escape successfully from the sample impoverishment problem as well as preserve the estimation accuracy and computation burden.

自举粒子滤波(BPF:Bootstrap Particle Filtering)是一种经典而应用广泛的粒子滤波算法,但其重采样后常会引起严重的样本枯竭问题.本文提出在权值蜕化较为严重时,在原先的重采样前增加SFN预处理,即权值排序、裂变繁殖(fission)和权值归一,得到裂变BPF(FBPF)算法.针对一个典型的后验密度为双峰的强非线性滤波估计问题,通过Monte Carlo仿真表明,FBPF算法在保持与BPF算法相当的估计精度和运算时间的条件下,克服了样本枯竭问题,算法的鲁棒性更强.