
2.4 一般系统噪声密度函数下状态估计的粒子滤波器设计
粒子滤波(Particle Filtering,PF)[1]主要使用样本点来逼近动态系统的非高斯性质,同时描述噪声的高阶矩信息。它不同于KF理论,对系统的扰动噪声没有高斯噪声这一要求,也没有限定系统中非线性函数的具体形式。随着计算机技术的发展,PF在经济、军事和雷达等领域得到了广泛的应用。
2.4.1 非线性系统描述
为简单起见,这里仅讨论下列情况下的非线性模型:

(2.4.1)

(2.4.2)
其中,w(k) 和v(k) 都是建模误差序列。状态过程随机变量具有概率密度分布
p(x)
(2.4.3)
在式(2.4.1)中

(2.4.4)
进一步,在式(2.4.2)中

(2.4.5)
另外,已知初始条件,即系统状态初始值x(0) 的统计特性如下:

(2.4.6)
在粒子滤波中,需要对

(2.4.7)
进行积分运算。其中,x∈R1p(x) 是待估计随机状态变量x的概率密度函数;f(x) 是关于p(x) 的可积函数。若随机采样N个服从p(x) 分布的样本粒子点

则概率密度函数就可以近似表示为p(x)≈pN(x),N→∞,即

(2.4.8)
其中,N是一个相当大的正整数;δ(x−xi) 是狄拉克函数。当粒子数目N趋近于无穷大时,式(2.4.7)可以近似为如下形式:

(2.4.9)
数学期望是每次实验的结果与该结果概率的乘积和,本质上是对函数与函数概率乘积的积分,即

(2.4.10)
若存在N个粒子点,则使用概率密度函数q(x) 来取点就可以将求取数学期望的过程进行简化。

(2.4.11)
使用蒙特卡罗方法,求取关于f(x0→k) 的条件期望。

(2.4.12)
设积分中的中间项为权重值,表示为wk(x0→k) ,那么当粒子点的数目为N时,条件期望就可以近似地以如下形式表示。

(2.4.13)
由上述推论,可以得到PF的滤波步骤如下。
第一步:基于条件概率密度的粒子点集采样。

(2.4.14)
第二步:采样时刻间PF粒子点权重递归转移更新。

(2.4.15)
第三步:当前时刻粒子点权重归一化计算。

(2.4.16)
第四步:计算当前时刻状态xk的估计值。

(2.4.17)
计算每个粒子点的权重乘积和,从而得到滤波器的最终状态估计值。
PF方法是一种非高斯滤波方法,适用于非高斯系统。PF方法也存在一些问题,如一些重要函数的选取方案、粒子多样性的缺乏、粒子退化、计算复杂度较高等问题。
对于x∈Rn,也有类似表述。
2.4.2 PF基本原理
PF方法采用采样逼近,其基本思想是构造一个基于采样样本的后验概率密度函数。利用N个粒子及相应的权重构成的集合,表示系统的后验概率密度函数p(x0:k|z1:k) 。其中,
是从后验概率分布状态空间中采样的集合。各个样本点的粒子权重取值为
,且粒子的权重和为1。基于粒子采样集合,k时刻的后验分布概率密度可以近似表示为[2-3]

(2.4.18)
通过这种逼近方法,能够将较为复杂的积分运算转化为求和运算。例如,若需要求解的函数的期望为

(2.4.19)
则基于采样点的逼近求解公式为

(2.4.20)
在实际应用系统中,直接从后验概率分布中获取有效的样本点是比较困难的。有效获取后验分布的样本,能够降低统计估计方差,这是提高PF方法滤波性能的关键。重要性采样方法(Importance Sampling Method)的引入能够有效提高采样效率。在此方法中,通过运用重要性采样密度q(x0:k|z1:k) 进行采样,可以有效地避免直接从后验概率密度中采样的困难。引入采样密度q(x0:k|z1:k) 后,式(2.4.19)可以改写为如下形式:

(2.4.21)
其中,。
根据式(2.4.21),从重要性采样密度中独立采样N个样本点,由式(2.4.20)可得

(2.4.22)
其中,为归一化后的权重。
可由下式计算:

(2.4.23)
假设在k−1 时刻,已获取后验概率密度p(x0:k−1|z1:k−1) 逼近的粒子集合,则下一步可以用一组新的样本集来逼近k时刻的后验概率密度p(x0:k|z1:k) 。
为方便进行滤波过程中的递归计算,将重要性密度函数分解如下:

(2.4.24)
然后,将从重要性采样密度中新获得的粒子,加入到已获取的粒子集
中,最后得到新的粒子集
。利用马尔可夫假设,上式可以进一步写成如下形式:

(2.4.25)
将式(2.4.25)代入式(2.4.23)可得如下形式:

(2.4.26)
其中,和
分别为似然函数和概率转移密度函数;
为重要性采样密度。由式(2.4.26)可知,若能够有效地选择q(⋅) ,则可以有效地递归更新粒子的权重值。后验滤波密度函数
可近似解析为如下形式:

(2.4.27)
根据以上分析,m阶PF方法可以总结如下。
第一步:从q(⋅) 中采样获取N个粒子:

(2.4.28)
然后得到新的粒子集:。
第二步:根据式(2.4.26),计算更新后的粒子对应的权重值。
第三步:权重归一化处理:

(2.4.29)
第四步:返回第一步。
对x∈Rn,也有类似描述。
2.4.3 PF仿真验证
考虑如下非线性系统:

(2.4.30)

(2.4.31)
并假设系统过程噪声w(k),k=0,1,2,⋯ 服从均值为0、协方差为1的均匀分布,观测噪声v(k),k=0,1,2,⋯ 服从均值为0、协方差为1的均匀分布,且w(k) 和v(k) 相互独立;初始状态x(0)=0.1 ,且系统初始状态和过程噪声及观测噪声相互独立,粒子数设定为N=100 。以下仿真结果基于50次蒙特卡罗仿真实验。
仿真结果如图2.4.1和图2.4.2所示。其中,图2.4.1显示的是状态值和估计值,图2.4.2显示的是状态估计误差。

图2.4.1 状态值和估计值

图2.4.2 状态估计误差
从图2.4.1可以看出,针对非线性较强的系统,采用PF能够有效地对目标状态进行估计,并且随着时间的推移,滤波效果并没有出现明显的发散现象。从图2.4.2能够看出状态估计误差随时间的变化情况,虽然在某些时刻误差较大,但这些时刻在整个滤波过程中所占的比例是比较小的,滤波的整体效果还是比较理想的。
2.4.4 滤波方法分析
PF方法是一种基于贝叶斯估计的滤波方法,其针对的是非线性高斯和非高斯系统,在解决非线性非高斯系统的参数估计和状态滤波问题方面有其自身的优势,因此广受国内外学者的青睐。但是相对于高斯系统中以KF为代表的滤波方法而言,PF方法还不够成熟,仍有大量的问题有待解决,主要表现在以下几个方面。
(1)重要性函数的选择会直接影响PF方法的性能。虽然一些学者针对实际滤波问题提出了改进的方法,但对通用系统缺乏统一的借鉴。一般根据特定系统的先验信息进行数学建模,在获得系统更多的信息后,再选择重要性函数。这已经成为对重要性函数进行选择和研究的主要方法。
(2)针对粒子匮乏现象的经典重采样虽然具有良好的效果,但随着粒子数的增加,计算复杂度将会呈几何级数增加,对于要求实时性的目标系统,PF方法在该系统中的实现将会成为关键的问题。在采样过程中,引入一些成熟的寻优方法,用于更快地获取能够反映系统概率特征的典型样本,将成为研究重采样方法的重点。
(3)从PF方法所涉及的数学基础来看,PF方法的收敛性问题目前还未彻底解决。若能在应用过程中有效解决其收敛性问题,将对粒子的匮乏现象起到抑制的作用。如何评估PF方法的估计性能也是一个需要认真考虑和解决的问题。
(4)将多种非线性滤波方法进行结合。虽然PF方法能够解决非线性高斯及非高斯系统中的滤波问题,但由于该方法存在运算实时性问题,以及初始状态概率选取问题,因此与工程应用尚有一定距离。根据系统在不同阶段、不同时刻表现出的不同统计特性,可将PF方法与其他非线性滤波方法进行结合,从而有效避免非线性系统的线性逼近。对于减小非线性系统模型线性化后的高阶截断误差,以及非高斯噪声带来的影响,这也是一种值得研究的方法。
(5)关于PF方法的硬件实现。PF方法的硬件实现对于提高该方法的运算速度和鲁棒性具有重要的作用。PF方法硬件实现的核心思想是:将PF方法划分为初始采样、重采样、状态更新等过程,利用流水线方法实现这些过程的并行处理。但能够应用于实际系统中的粒子滤波器尚未研制成功。
(6)将PF方法推广到新的应用领域。PF方法出现较晚,目前仅应用于有限的几个领域中。对于非线性系统中的滤波问题,在卡尔曼滤波方法应用受限的情况下,可尝试采用PF方法。
目前,PF方法在国内外广受关注且发展很快,在统计学、模式识别、人工智能、计算机视觉、经济学、通信、信号处理等领域取得了一些研究成果。若能有效解决上述问题,必将有效促进PF方法理论和应用的发展。