适用于非线性系统状态的粒子估计方法

摘要

适用于非线性系统状态的粒子估计方法。本发明属于工程领域，具体涉及到非线性系统状态估计的技术领域。它是为解决对非线性系统状态进行估计时，因系统状态观测似然函数具有双峰特性使状态估计精度较差的问题。它在标准粒子滤波方法的基础上，在粒子权值更新步骤中将粒子观测向量和系统状态观测向量的夹角作为计算粒子权值的参数，并在系统状态估计步骤中加入了平滑操作。本发明适用于一般非线性系统的状态估计问题，当系统状态观测似然函数具有双峰特性时，该方法的状态估计精度较高。

说明书

技术领域

本发明属于工程领域，具体涉及的是非线性系统状态估计的技术领域。

背景技术

非线性系统状态估计问题广泛存在于信号处理及其相关领域中，例如视觉跟踪、语音识别、机器人定位、通信和故障检测等。尽管经过了近半个世纪的发展，但非线性系统状态估计方法还较不成熟。一方面，这是由于非线性系统的复杂性、多样性所造成的；另一方面，由于非线性系统中各种随机噪声的存在，导致状态估计更加困难，影响了状态估计的准确率，为非线性系统状态估计带来了较大的困难。20世纪80年代以来，非线性系统状态估计方法受到越来越多的关注，其中最常用的是扩展卡尔曼滤波器。但扩展卡尔曼滤波仅对某些特定的非线性系统有效，对于一般的非线性系统扩展卡尔曼滤波不能保证其收敛且滤波误差较大。近年来一种在非线性系统状态估计方面具有独到优势的新滤波方法一粒子滤波得到了迅速的发展，它有效克服了扩展卡尔曼滤波的缺点，并已被成功应用于众多非线性系统状态估计问题中，例如，美国国家航空和宇宙航行局研制的“火星漫游者”机器人、舰船涡轮发动机、涡轮喷气飞机发动机等。粒子滤波是递归贝叶斯状态估计的一种模拟实现方法。它利用大量随机样本(或称为粒子)代表待估计状态的后验概率密度函数。这些粒子通过系统模型按时间顺序向前传播，从而得到大量系统状态的路径样本。各时刻系统状态的后验概率密度函数用这些粒子来代表。目前比较有实用价值的粒子滤波方法有标准粒子滤波方法(SIR)、辅助粒子滤波方法(APF)、正则化粒子滤波方法(RPF)和高斯粒子滤波方法(GPF)。其中，辅助粒子滤波方法、正则化粒子滤波方法和高斯粒子滤波方法均是标准粒子滤波方法的改进方法。

如果用x_k表示k时刻的系统状态，用y_k表示x_k的观测值，用{x(i)：i＝1，2，...，N}代表k时刻的所有粒子，用p(y|x)表示观测似然函数，则标准粒子滤波方法的步骤可描述为：

步骤一、粒子预测步骤--将k时刻系统状态先验分布作为重要性密度函数并从中抽取随机粒子{x(i)：i＝1，2，...N}；

步骤二、粒子权值更新步骤-得到k时刻系统状态观测值y_k后，利用观测似然函数计算k时刻上述各粒子的似然函数值并将其作为k时刻各粒子的权值即{q_k(i)＝p(y_k|x_k(i))：i＝1，2，...，N}；

步骤三、系统状态估计步骤--利用这些粒子及其权值组成的离散概率分布作为k时刻系统状态后验概率分布，去近似k时刻系统状态的真实概率分布，并用于k时刻系统状态估计；

步骤四、重采样步骤--在上述离散概率分布中重新采样得到新的粒子群，各粒子的权值等于1/N，然后返回步骤一。

辅助粒子滤波方法在标准粒子滤波方法的基础上，采用辅助方法预先从前一时刻的粒子样本中选出较好的预测粒子，以期其粒子组成的离散概率分布接近系统状态真实概率分布的可能性更高。实际应用中该方法的估计精度在系统噪声较小时好于标准粒子滤波方法，但当系统噪声较大时精度不如标准粒子滤波方法。正则化粒子滤波方法改变了标准粒子滤波方法步骤四中的重采样。不像标准粒子滤波方法在离散概率分布上重采样，正则化粒子滤波方法利用正则化方法获得一个连续概率分布并从中重采样。该方法改善了标准粒子滤波方法中存在的粒子多样性匮乏问题，但它有一个理论上的缺陷，即所得的粒子离散概率分布不再能保证是渐进地收敛到系统状态真实概率分布。实际应用中当粒子多样性匮乏问题较严重时该方法好于标准粒子滤波方法。高斯粒子滤波方法在标准粒子滤波方法的基础上，采用一个高斯分布来近似系统状态真实概率分布，并且没有重采样步骤。当系统状态真实概率分布服从高斯分布时该方法的状态估计精度优于标准粒子滤波方法。

综上所述，以上四种方法都具有各自的优点，在其所适合的应用范围内都是比较有效的。然而在实际非线性系统状态估计中，系统状态观测似然函数常常具有双峰特性，这通常导致上述方法的状态估计精度较差。

发明内容

本发明是为解决现有粒子滤波方法在对非线性系统状态进行估计时，因系统状态观测似然函数具有双峰特性而导致估计精度较差的问题，进而提供的一种适用于非线性系统状态的粒子估计方法。

本发明的步骤为：

步骤一、粒子预测步骤：将k时刻系统状态先验分布作为重要性密度函数并从中抽取随机粒子{x(i)：i＝1，2，...，N}；

步骤二、粒子权值更新步骤：

(1)得到k时刻系统状态观测值y_k后，利用观测似然函数计算k时刻各粒子的似然函数值为{p(y_k|x_k(i))}，其中i＝1，2，...，N；

(2)计算各粒子观测向量V(i)和系统状态观测向量V之间的夹角θ(i)，其中i＝1，2，...，N；

(3)计算k时刻各粒子的权值q_k(i)如下：

q_k(i)＝e^α/θ(i)×p(y_k|x_k(i))，i＝1，2，..，N

计算k-L+l时刻各粒子的权值q_k-L+l(i)如下：

q_k-L+l(i)＝e^α/θ(i)×p(y_k-L+l|x_k-L+l(i))，i＝1，2，..，N

其中，α是一个预先指定的比例因子；

步骤三、系统状态估计步骤：

(1)k时刻系统状态的滤波估计--利用k时刻粒子及其权值组成的离散概率分布作为k时刻系统状态后验概率分布，去近似k时刻系统状态的真实概率分布，并用于k时刻系统状态估计；

(2)k-L+l时刻系统状态的平滑估计--利用k-L+l时刻粒子及其权值组成的离散概率分布作为k-L+l时刻系统状态后验概率分布，去近似k-L+l时刻系统状态的真实概率分布，并用于k-L+l时刻系统状态估计；

(3)返回步骤一；

用x_k表示k时刻的系统状态，用y_k表示x_k的观测值，用{x(i)：i＝1，2，...，N}代表k时刻的所有粒子，用p(y|x)表示观测似然函数，用{x_j(i)：j＝k-L+1，...，k}代表k-L+1时刻到k时刻的所有粒子，用V(i)＝{y_k-L+1(i)，y_k-L+2(i)，...，y_k(i)}表示粒子观测向量，用V＝{y_k-L+1，y_k-L+2，...，y_k}表示系统状态观测向量，用θ(i)表示这两种向量的夹角。这里，L和l是两个预先指定的常量且1≤l＜L，α是一个预先指定的比例因子。

本发明的粒子权值更新步骤和系统状态估计步骤不同于现有的粒子滤波方法，本发明方法在粒子权值更新步骤中将粒子观测向量和系统状态观测向量的夹角作为计算粒子权值的参数，并在系统状态估计步骤中加入了平滑操作。这些特征都是现有的粒子滤波方法所不具备的。本发明适用于一般非线性系统状态估计问题的求解，在系统状态观测似然函数具有双峰特性时该方法优于现有粒子滤波方法，能有效地提高非线性系统状态的估计精度。

附图说明

图1是本发明的方法与SIR、APF、RPF和GPF算法的性能对比实验效果图。

具体实施方式

具体实施方式一：

本发明基于下述理论基础：如果一个粒子的观测路径靠近系统状态的观测路径，则该粒子的路径应靠近系统状态的路径。如果将一个粒子的一系列观测值和系统状态的一系列观测值分别看作一个粒子观测向量和系统状态观测向量，则可用这两个向量的夹角来度量它们之间的相似程度。当该夹角较小时就表示该粒子靠近系统状态。

用x_k表示k时刻的系统状态，用y_k表示x_k的观测值，用{x(i)：i＝1，2，...，N}代表k时刻的所有粒子，用p(y|x)表示观测似然函数，用{x_j(i)：j＝k-L+1，...，k}代表k-L+1时刻到k时刻的所有粒子，用V(i)＝{y_k-L+1(i)，y_k-L+2(i)，...，y_k(i)}表示粒子观测向量，用V＝{y_k-L+1，y_k-L+2，...，y_k}表示系统状态观测向量，用θ(i)表示这两种向量的夹角。这里，L和l是两个预先指定的常量且1≤l＜L，α是一个预先指定的比例因子。

本实施方式的非线性系统状态粒子估计方法步骤为：

步骤一、粒子预测步骤：将k时刻系统状态先验分布作为重要性密度函数并从中抽取随机粒子{x(i)：i＝1，2，...，N}；

步骤二、粒子权值更新步骤：

(1)得到k时刻系统状态观测值y_k后，利用观测似然函数计算k时刻各粒子的似然函数值为{p(y_k|x_k(i))}，其中i＝1，2，...，N；

(2)计算各粒子观测向量V(i)和系统状态观测向量V之间的夹角θ(i)，其中i＝1，2，...，N；

(3)计算k时刻各粒子的权值q_k(i)如下：

q_k(i)＝e^α/θ(i)×p(y_k|x_k(i))，i＝1，2，..，N

计算k-L+l时刻各粒子的权值q_k-L+l(i)如下：

q_k-L+l(i)＝e^α/θ(i)×p(y_k-L+l|x_k-L+l(i))，i＝1，2，..，N

其中，α是一个预先指定的比例因子；

步骤三、系统状态估计步骤：

(1)k时刻系统状态的滤波估计--利用k时刻粒子及其权值组成的离散概率分布作为k时刻系统状态后验概率分布，去近似k时刻系统状态的真实概率分布，并用于k时刻系统状态估计；

(2)k-L+l时刻系统状态的平滑估计--利用k-L+l时刻粒子及其权值组成的离散概率分布作为k-L+l时刻系统状态后验概率分布，去近似k-L+l时刻系统状态的真实概率分布，并用于k-L+l时刻系统状态估计；

(3)返回步骤一。

本具体实施方式的非线性系统状态粒子估计方法步骤二中第(2)子步骤所述的向量V(i)和V之间的夹角θ(i)的计算步骤为：

步骤1、计算V(i)和V的内积如下

(V(i)，V)＝y_k-L+1(i)·y_k-L+1+y_k-L+2(i)·y_k-L+2+…+y_k(i)·y_k

步骤2、分别计算V(i)和V的范数如下

$\left|\right|V\left(i\right)\left|\right|=\sqrt{{\left(y_{k-L+1}\left(i\right)\right)}^2+{\left(y_{k-L+2}\left(i\right)\right)}^2+...+{\left(y_k\left(i\right)\right)}^2}$

$\left|\right|V\left|\right|=\sqrt{{\left(y_{k-L+1}\right)}^2+{\left(y_{k-L+2}\right)}^2+...+{\left(y_k\right)}^2}$

步骤3、计算V(i)和V之间的夹角θ(i)如下

$\theta \left(i\right)=\arccos \left(\frac{(V\left(i\right),V)}{\left|\right|V\left(i\right)\left|\right|\times \left|\right|V\left|\right|}\right)$

其中，arccos(.)是反余弦函数。

本方法在被应用于如下具有代表性的观测似然函数有双峰特性的非线性系统中时，本方法的误差远小于SIR、APF、RPF和GPF方法。

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

$y_k=\frac{x_k^2}{20}+v_k$

使用粒子总数N＝100、u_k-1～N(0，1)、v_k～N(0，1)，x₀＝0.1；一次实验运行时间为5000个时间点，共进行100次实验运行；取L＝5、I＝2和α＝1.8；采用误差均方根值(RMSE)评价算法性能。实验结果显示在图1中，其中PE代表本方法结果，SIR代表标准粒子滤波方法，APF代表辅助粒子滤波方法，RPF代表正则化粒子滤波方法，GPF代表高斯粒子滤波方法。