一种基于直方图估计粒子滤波算法的非线性动态系统信号处理方法

摘要

一种基于直方图估计粒子滤波算法的非线性动态系统信号处理方法，设非线性动态系统的状态空间模型为：x

说明书

技术领域

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

背景技术

动态系统的状态估计问题涉及很多领域，尤其是信号处理、人工智能和计算机视觉领域。传统的卡尔曼滤波只适用于线性高斯系统，而扩展卡尔曼滤波也只能应对系统的弱非线性。因此，适用于非线性、非高斯系统的粒子滤波备受关注。

粒子滤波是一种基于蒙特卡洛模拟和递推贝叶斯估计的滤波方法。它采用粒子描述状态空间，使用一组带权重的粒子近似表示系统的后验概率密度，并通过模型方程和观测信息实现递推的估计过程。常见的粒子滤波算法包括SIR粒子滤波、辅助粒子滤波(APF)、正则粒子滤波(RPF)、高斯粒子滤波(GPF)和无迹粒子滤波(UPF)。

粒子滤波算法中的两个关键技术是建议分布的选取和重采样算法。常用的建议分布虽然容易获取，但易造成粒子权重的急剧退化，使得有效粒子数大大减少，尤其是当似然分布位于先验分布的尾部时。同时，虽然重采样算法解决了粒子权重的退化问题，但常见的重采样算法只是对粒子的简单复制和剔除，从而又导致了粒子匮乏问题(粒子多样性的丧失)。

为了克服现有基于粒子滤波算法的动态系统信号处理方法中的由于粒子权重易急剧退化和粒子匮乏，不能有效代表后验密度，因而无法有效处理非线性和非高斯信号的不足，本发明提出一种能更好地估计系统的后验概率密度、粒子多样性好、能有效处理非线性和非高斯信号的基于直方图估计粒子滤波算法的非线性动态系统信号处理方法。

为了解决上述技术问题提出的技术方案为：

一种基于直方图估计粒子滤波算法的非线性动态系统信号处理方法，设非线性动态系统的状态空间模型为：

x_k＝f(x_k-1)+v_k-1

z_k＝h(x_k)+n_k

其中，x_k和z_k分别表示系统在k时刻的状态和观测值，f(x_k-1)和h(x_k)分别表示系统的状态转移方程和观测方程，x_k-1表示系统在k时刻的状态，x_k-1和n_k分别表示系统噪声和观测噪声；

所述非线性动态系统信号处理方法包括以下步骤：

第一步，根据k-1时刻的N个粒子i＝1，2，...，N，通过状态转移方程得到k时刻的N个预测粒子i＝1，2，...，N；

第二步，确定预测粒子所占据的状态空间，并将其均分为箱，箱的个数记为M；

第三步，根据划分的箱估计系统的先验概率密度，各箱对应的密度d_i为：

$d^i=\frac{N^i}{N\times h},i=\mathrm{1,2},...,M$

其中，Nⁱ表示落入对应箱的预测粒子数，h表示一个箱所占据的超空间的体积。

第四步，估计各箱对应的后验概率密度wⁱ：

wⁱ＝dⁱ×p(z_k|xⁱ)，i＝1，2，...，M

其中，xⁱ表示第i个箱的中心；

第五步，对各箱的后验概率密度进行归一化：

$w^i=\frac{w^i}{{\Sigma }_{j=1}^M{w}^j},i=\mathrm{1,2}...,M$

第六步，从M个箱中重采样N个粒子作为k时刻系统的后验密度分布，i＝1，2，...，N，其中，同一个箱中重采样的粒子在对应的状态空间中进行均匀采样；

第七步，输出系统状态的估计值

$x_k^u\approx {\Sigma }_{i=1}^M{w}^i\times x^i.$

本发明的技术构思为：该算法的一个递推过程包括以下基本步骤：

1)、根据k-1时刻的N个粒子得到k时刻的N个预测粒子。

2)、确定预测粒子所占据的状态空间，并将其均分为箱。

3)、根据划分的箱估计系统的先验概率密度。

4)、估计各箱的后验概率密度。

5)、对各箱的后验概率密度归一化。

6)、重采样。

7)、输出。

本发明具有以下优点：

1、采用直方图估计中的箱来估计系统的先验和后验概率密度，避免了粒子权重的计算，更有效地代表了系统的后验概率密度，且减少了计算量。

2、在箱中进行重采样，避免了粒子匮乏问题。

附图说明

图1为直方图估计粒子滤波算法的流程图。

图2为直方图估计粒子滤波算法与其他粒子滤波算法在不同粒子数下的均方根误差的均值对比曲线示意图。

图3为直方图估计粒子滤波算法与其他粒子滤波算法在不同粒子数下的运行时间对比曲线示意图。

具体实施方式

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

参照图1～图3，一种基于直方图估计粒子滤波算法的非线性动态系统信号处理方法，采用粒子描述动态系统的状态空间，设非线性动态系统的状态空间模型为：

x_k＝f(x_k-1)+v_k-1

z_k＝h(x_k)+n_k

其中，x_k和z_k分别表示系统在k时刻的状态和观测值，f(x_k-1)和h(x_k)分别表示系统的状态转移方程和观测方程，x_k-1表示系统在k时刻的状态，v_k-1和n_k分别表示系统噪声和观测噪声；

非线性动态系统信号处理方法包括以下步骤：

第一步，根据k-1时刻的N个粒子i＝1，2，...，N，通过状态转移方程得到k时刻的N个预测粒子i＝1，2，...，N；

第二步，确定预测粒子所占据的状态空间，并将其均分为箱，箱的个数记为M；

第三步，根据划分的箱估计系统的先验概率密度，各箱对应的密度dⁱ为：

$d^i=\frac{N^i}{N\times h},i=\mathrm{1,2},...,M$

其中，Nⁱ表示落入对应箱的预测粒子数，h表示一个箱所占据的超空间的体积。

第四步，估计各箱对应的后验概率密度wⁱ：

wⁱ＝dⁱ×p(z_k|xⁱ)，i＝1，2，...，M

其中，xⁱ表示第i个箱的中心；

第五步，对各箱的后验概率密度进行归一化：

$w^i=\frac{w^i}{{\Sigma }_{j=1}^M{w}^j},i=\mathrm{1,2}...,M$

第六步，从M个箱中重采样N个粒子作为k时刻系统的后验密度分布，i＝1，2，...，N，其中，同一个箱中重采样的粒子在对应的状态空间中进行均匀采样；

第七步，输出系统状态的估计值

$x_k^u\approx {\Sigma }_{i=1}^M{w}^i\times x^i.$

本实施例通过一个非线性动态系统的状态估计对本发明和其它几种粒子滤波算法进行比较。系统的状态空间模型如下：

$x_k=\frac{x_{k-1}}{2}+\frac{{25x}_{k-1}}{1+x_{k-1}^2}+8\cos \left(1.2k\right)+v_{k-1}$

$z_k=\frac{x_k^2}{20}+n_k$

其中，系统噪声方差和观测噪声方差分别取为10和1。设定观测时间为100，运行次数为100，在粒子数N分别取100、200、300、400、500、600时，本发明提出的算法与其他粒子滤波算法所产生的均方根误差(RMSE)均值和所运行的时间分别如图2和图3所示。

从图2和图3中可以看出，本发明算法在不同粒子数目下的RMSE均值明显优于除UPF外的其它算法，而运行时间也显著降低。当N＝600时，本发明算法的运行时间只有SIR的约57％，UPF的约10％。