CN105930640A - 一种处理Lévy噪声的分数阶卡尔曼滤波方法 - Google Patents

Abstract

本发明公开了一种处理Lévy噪声的分数阶卡尔曼滤波方法。首先，给出状态预测初值和预测误差协方差初值。接着，对Lévy噪声进行处理，计算当前时刻量测噪声协方差和当前时刻最优滤波增益。然后，利用当前时刻最优滤波增益和状态预测值计算当前时刻状态估计值，利用预测误差协方差计算当前时刻估计误差协方差，再利用当前时刻状态估计值和估计误差协方差计算当前时刻系统噪声协方差。最后，利用当前时刻状态估计值对下一时刻状态预测值进行更新，利用当前时刻估计误差协方差对下一时刻预测误差协方差进行更新。本发明能解决分数阶线性离散系统在非高斯噪声下的状态估计问题，而且易于与已有的状态估计软件相结合。

Description

一种处理Lévy噪声的分数阶卡尔曼滤波方法

技术领域

本发明涉及一种处理Lévy噪声的分数阶卡尔曼滤波方法，属于系统分析与处理技术领域。

背景技术

系统分析与处理，旨在研究特定系统结构中各部分(各子系统)的相互作用，系统的对外接口与界面，以及系统整体的行为、功能和局限，从而为系统未来的变迁与有关决策提供参考和依据，其目标之一在于改善决策过程及系统性能，以达到系统的整体最优。在系统分析与处理领域，状态估计起着至关重要的作用。依观测数据与被估计状态在时间上的相对关系，状态估计又可分为平滑、滤波和预测3种情形。为了估计t时刻的状态x(t)，如果可用的信息包括t以前的观测值，就是平滑问题。如果可用的信息是时刻t的观测值，估计可实时地进行，称为滤波问题。如果必须时刻(t-Δ)以前的观测来估计经历了Δ时间之后的状态x(t)，则是预测问题。传统的状态估计方法只考虑到整数阶次的系统，大多用来描述平滑问题和滤波问题，而对预测问题则不能进行描述；分数阶次具有全局相关性，能较好地体现系统函数发展的历史依赖过程，因此能处理预测问题。

状态估计是卡尔曼滤波的重要组成部分，其对于了解和控制一个系统具有重要意义。传统的卡尔曼滤波方法要求系统为整数阶次，并且系统噪声和量测噪声均为高斯白噪声，这些过于理想化的要求在实际系统中很难得到满足。相反，实际存在的系统往往受到非高斯噪声的影响，而且并非全部系统都能以整数阶来建模。

发明内容

发明目的：基于以上分析，本发明采用分数阶理论，处理非高斯Lévy噪声下的状态估计问题，以期提高状态估计的估计精度，进而提高整个数据系统的质量和可靠性。

由于非高斯噪声普遍存在于实际系统中，那么利用传统的卡尔曼滤波方法对这些系统进行状态估计时不能得出理想结果。一种处理Lévy噪声的分数阶卡尔曼滤波方法首先对Lévy噪声进行处理，计算当前时刻量测噪声协方差和当前时刻最优滤波增益。然后，利用当前时刻最优滤波增益和状态预测值计算当前时刻状态估计值，利用预测误差协方差计算当前时刻估计误差协方差，再利用当前时刻状态估计值和估计误差协方差计算当前时刻系统噪声协方差。最后，利用当前时刻状态估计值对下一时刻状态预测值进行更新，利用当前时刻估计误差协方差对下一时刻预测误差协方差进行更新。该方法与常规动态状态估计程序相结合，能很好地处理非高斯噪声情况下的状态估计问题。

技术方案：一种处理Lévy噪声的分数阶卡尔曼滤波方法，所述方法是在计算机中依次按以下步骤实现的：

(1)、初始化，包括：设定状态预测量的初值和预测误差协方差的初值；

(2)、对当前时刻系统噪声和量测噪声进行处理，处理步骤为：

$`w_k=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_1\&CenterDot;sign\left({\stackrel{\&OverBar;}{w}}_k\right)& if& \left|{\stackrel{\&OverBar;}{w}}_k\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ {\stackrel{\&OverBar;}{w}}_k& if& \left|{\stackrel{\&OverBar;}{w}}_k\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

$`v_k=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_2\&CenterDot;sign\left({\stackrel{\&OverBar;}{v}}_k\right)& if& \left|{\stackrel{\&OverBar;}{v}}_k\right|\&GreaterEqual;{\mathrm{\&delta;}}_2\\ {\stackrel{\&OverBar;}{v}}_k& if& \left|{\stackrel{\&OverBar;}{v}}_k\right|<{\mathrm{\&delta;}}_2\end{array}\right.$

式中，`w<sub>k</sub>表示系统噪声的近似值，下标k表示第k时刻，δ<sub>1</sub>表示选取的第一个阈值，sign(x)表示符号函数，如果x＞0，返回1；如果x＜0，返回-1。表示含Lévy序列的系统噪声，|x|表示取x的绝对值，`v_k表示量测噪声的近似值，δ₂表示选取的第二个阈值，表示含Lévy序列的量测噪声；

(3)、利用系统噪声和量测噪声的近似值，带入原系统得到新的分数阶线性系统模型，新模型可表示为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}={\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w_k\\ x_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j}\\ y_k={\mathrm{Cx}}_k+`v_k\end{array}\right.$

式中，Δ表示分数阶算子，上标α表示分数阶阶次向量，x_k+1表示第k+1时刻状态向量，A、B和C表示适当维数的已知矩阵，u_k表示第k时刻控制输入向量，`w<sub>k</sub>表示第k时刻系统噪声的近似值，y<sub>k</sub>表示第k时刻量测向量，`v_k表示第k时刻量测噪声的近似值，γ_j表示

其中

α_N表示第N个分数阶阶次值，j＝1，2，…，k+1；

(4)、计算当前时刻量测噪声协方差，计算步骤为：

$`R_k=\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;{\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;}^T+C{\tilde{p}}_k{C}^T$

式中，`R_k表示第k时刻量测噪声协方差，y_k表示第k时刻量测向量，C表示适当维数的已知矩阵，上标T表示转置，表示第k时刻状态的预测量，表示第k时刻状态预测误差协方差；

(5)、计算当前时刻最优滤波增益，计算步骤为：

$K_k={\tilde{p}}_k{C}^T{(C{\tilde{p}}_k{C}^T+`R_k)}^{-1}$

式中，K_k表示第k时刻最优滤波增益；

(6)、计算当前时刻状态估计值，计算步骤为：

${\hat{x}}_k={\tilde{x}}_k+K_k\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;$

(7)、计算当前时刻估计误差协方差，计算步骤为：

${\hat{p}}_k=(I-K_{k}C){\hat{p}}_k$

(8)、计算当前时刻系统噪声协方差，计算步骤为：

$\begin{array}{c}`Q_k=(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\left(-1\right)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}){(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\left(-1\right)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j})}^T\\ +(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

式中，`Q_k表示第k时刻系统噪声协方差，表示第k时刻状态估计向量，u_k表示控制输入向量，x_k+1表示第k+1时刻状态向量，表示第k+1-j时刻状态估计误差协方差；

(9)、利用当前时刻状态估计值更新下一时刻状态预测值，计算步骤为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\tilde{x}}_{k+1}=A{\tilde{x}}_k+{\mathrm{Bu}}_k\\ {\tilde{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\tilde{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}\end{array}\right.$

式中Δ表示分数阶算子，上标α表示分数阶阶次向量；

(10)、利用当前时刻估计误差协方差更新下一时刻状态预测协方差，计算步骤为：

${\tilde{p}}_{k+1}=(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+`Q_k+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T$

(11)、判断k+1是否大于等于步长L(可自己设定)，如果是，结束计算，否则返回步骤(2)进行下一次估计。

本发明提出的一种处理Lévy噪声的分数阶卡尔曼滤波方法，在传统卡尔曼滤波的基础上运用分数阶理论，并且通过非高斯Lévy噪声替代传统的高斯白噪声，利用近似代替方法分解得到Lévy噪声的分解值，重新计算系统噪声和量测噪声的协方差，进而得到当前时刻状态估计值和估计误差协方差，最后对下一时刻状态预测值和预测误差协方差进行更新。本发明由于结合分数阶理论和非高斯Lévy噪声，能有效地解决整数阶系统的局限性问题和非高斯噪声下的状态估计问题。

附图说明

图1为本发明实施例的方法流程图；

图2为本发明提出的方法得到的仿真值。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，一种处理Lévy噪声的分数阶卡尔曼滤波方法，包括如下步骤：

(1)、初始化，包括：设定状态预测量的初值和预测误差协方差的初值。

(2)、对当前时刻系统噪声和量测噪声进行处理。

(3)、利用系统噪声和量测噪声的近似值，带入原系统得到新的分数阶线性系统模型。

(4)、计算当前时刻量测噪声协方差和当前时刻最优滤波增益。

(5)、利用当前时刻最优滤波增益和状态预测值计算当前时刻状态估计值。

(6)、利用当前预测误差协方差和最优滤波增益计算当前时刻估计误差协方差。

(7)、利用当前时刻状态估计值和估计误差协方差计算当前时刻系统噪声协方差。

(8)、利用当前时刻状态估计值更新下一时刻状态预测值。

(9)、利用当前时刻估计误差协方差更新下一时刻状态预测协方差

(10)、判断k+1是否大于等于步长L，如果是，结束计算，否则返回步骤(2)进行下一次估计。

目前分数阶线性系统状态估计模型主要采用的是20世纪70年代初由R.E.Larson和Debs提出的卡尔曼滤波(Kalman Filter，KF)算法发展而来的分数阶卡尔曼滤波(Fractional Kalman Filter，FKF)算法。卡尔曼滤波是一种递推回归方法，它的基本思想是根据新量测的数据和由前面一步量测数据计算而得的估计值，再来推算出新的估计值，即

新估计值＝旧估计值十修正值

分数阶卡尔曼滤波的本质是求出k+1时刻状态量真值x_k+1的最优估计值估计的准则是以状态量的估计误差协方差阵最小为目标函数，即：

$\min {\hat{P}}_{k+1|}=\min E\&lsqb;{\hat{e}}_{k+1}{\hat{e}}_{k+1}^T\&rsqb;$

式中E为求期望函数，x_k+1为k+1时刻状态量的真值，为k+1时刻状态量的估计值。

由于实际中的一些系统具有分数阶性质，并且存在非高斯白噪声，因此传统卡尔曼滤波算法在应用上受到了很大的限制，分数阶卡尔曼滤波算法具有全局相关性，能较好地体现系统函数发展的历史依赖过程，因而能处理具有分数阶性质的系统。对于具有非高斯Lévy噪声的分数阶线性离散系统，模型描述可以表示为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}=A{\stackrel{\&OverBar;}{x}}_k+{\mathrm{Bu}}_k+{\stackrel{\&OverBar;}{w}}_k\\ {\stackrel{\&OverBar;}{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}\\ {\stackrel{\&OverBar;}{y}}_k={\mathrm{Cx}}_k+{\stackrel{\&OverBar;}{v}}_k\end{array}\right.$

式中Δ表示分数阶算子，α表示分数阶阶次向量，表示第k+1时刻含有Lévy噪声的状态向量，A、B和C表示适当维数的已知矩阵，u_k表示第k时刻控制输入矩阵，表示第k时刻量测输出向量，和均表示非高斯Lévy噪声，两者协方差矩阵分别为和γ_j表示

其中

式中α_N表示第N个分数阶阶次值，j＝1，2，…，k+1。

在以上所示分数阶线性离散系统模型中，系统噪声和量测噪声均为非高斯Lévy噪声，由于其具有无限方差特性，因此不能用来直接计算协方差矩阵，必须对它们进行近似处理。处理方法为：

$`w_k=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_1\&CenterDot;sign\left({\stackrel{\&OverBar;}{w}}_k\right)& if& \left|{\stackrel{\&OverBar;}{w}}_k\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ {\stackrel{\&OverBar;}{w}}_k& if& \left|{\stackrel{\&OverBar;}{w}}_k\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

$`v_k=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_2\&CenterDot;sign\left({\stackrel{\&OverBar;}{v}}_k\right)& if& \left|{\stackrel{\&OverBar;}{v}}_k\right|\&GreaterEqual;{\mathrm{\&delta;}}_2\\ {\stackrel{\&OverBar;}{v}}_k& if& \left|{\stackrel{\&OverBar;}{v}}_k\right|<{\mathrm{\&delta;}}_2\end{array}\right.$

式中，`w<sub>k</sub>表示系统噪声的近似值，下标k表示第k时刻，δ<sub>1</sub>表示选取的第一个阈值，sign(x)表示符号函数，如果x＞0，返回1；如果x＜0，返回-1。表示含Lévy序列的系统噪声，|x|表示取x的绝对值，`v_k表示量测噪声的近似值，δ₂表示选取的第二个阈值，表示含Lévy序列的量测噪声。

对系统噪声和量测噪声近似处理后，得到新分数阶非线性系统模型，可表示为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}={\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w_k\\ x_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j}\\ y_k={\mathrm{Cx}}_k+`v_k\end{array}\right.$

式中，Δ表示分数阶算子，上标α表示分数阶阶次向量，x_k+1表示第k+1时刻状态向量，A、B和C表示适当维数的已知矩阵，u_k表示第k时刻控制输入向量，`w<sub>k</sub>表示第k时刻系统噪声的近似值，y<sub>k</sub>表示第k时刻量测向量，`v_k表示第k时刻量测噪声的近似值。

第k时刻量测噪声协方差可分别表示为：

$`R_k=E\&lsqb;`v_k`v_k^T\&rsqb;$

利用新模型中的表达式进行代换，可得到：

`v_k＝y_k-Cx_k

将上式代入量测噪声协方差阵可得：

$`R_k=\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;{\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;}^T+C{\tilde{p}}_k{C}^T$

式中，表示第k时刻状态的预测量，上标T表示转置，表示第k时刻状态预测误差协方差。

由已得到的第k时刻量测噪声协方差，计算第k时刻最优滤波增益，计算公式为：

$K_k={\tilde{p}}_k{C}^T{(C{\tilde{p}}_k{C}^T+`R_k)}^{-1}$

系统在第k时刻状态估计量和估计误差协方差可分别表示为：

${\hat{x}}_k=E\&lsqb;x_k|y_k^{*}\&rsqb;$

${\hat{P}}_k=\mathrm{cov}\&lsqb;x_k|y_k^{*}\&rsqb;$

式中，表示k-1时刻之前所有控制量和量测量的累计值，E表示取期望函数，cov表示取协方差函数。

根据上述两式的定义公式，可得到具体表示式为：

$\begin{array}{c}{\hat{x}}_k=E\&lsqb;x_k|y_k^{*}\&rsqb;\\ =E\&lsqb;{\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w-\underset{j=1}{\overset{k}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k-j}|y_k^{*}\&rsqb;\\ ={\tilde{x}}_k+K_k\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;\end{array}$

$\begin{array}{c}{\hat{P}}_k=\mathrm{cov}\&lsqb;x_k|y_k^{*}\&rsqb;\\ =\mathrm{cov}\&lsqb;{\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w-\underset{j=1}{\overset{k}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k-j}|y_k^{*}\&rsqb;\\ =(I-K_{k}C){\tilde{p}}_k\end{array}$

第k时刻系统噪声协方差可分别表示为：

$`Q_k=E\&lsqb;`w_k`w_k^T\&rsqb;$

利用新模型中的表达式进行代换，可得到：

$`w=x_{k+1}-{\mathrm{Ax}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j}$

将上式代入系统噪声协方差阵可得：

$\begin{array}{c}`Q_k=(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\left(-1\right)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}){(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\left(-1\right)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j})}^T\\ +(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

式中，表示第k时刻状态估计向量，u_k表示控制输入向量，表示第k时刻状态估计误差协方差。

由k时候变到k+1时刻时，k+1时刻状态预测量和预测误差协方差分别为：

${\tilde{x}}_{k+1}=E\&lsqb;x_{k+1}|y_k^{*}\&rsqb;$

${\tilde{P}}_{k+1}=\mathrm{cov}\&lsqb;x_{k+1}|y_k^{*}\&rsqb;$

具体计算公式如下：

$\begin{array}{c}{\tilde{x}}_{k+1}=E\&lsqb;x_{k+1}|y_k^{*}\&rsqb;\\ =E\&lsqb;{\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j}|y_k^{*}\&rsqb;\\ =A{\hat{x}}_k+{\mathrm{Bu}}_k-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}\end{array}$

$\begin{array}{c}{\tilde{P}}_{k+1}=\mathrm{cov}\&lsqb;x_{k+1}|y_k^{*}\&rsqb;\\ =\mathrm{cov}\&lsqb;{\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j}|y_k^{*}\&rsqb;\\ =(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+`Q_k+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

一种处理Lévy噪声的分数阶卡尔曼滤波方法的完整计算公式如下：

a.预测步：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\tilde{x}}_{k+1}=A{\tilde{x}}_k+{\mathrm{Bu}}_k\\ {\tilde{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\tilde{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}\\ {\tilde{p}}_{k+1}(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+`Q_k+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}\right.$

b.估计步：

$\left\{\begin{array}{c}{K}_{k+1}={\tilde{p}}_{k+1}{C}^T{(C{\tilde{p}}_{k+1}{C}^T+`R_{k+1})}^{-1}\\ {\hat{x}}_{k+1}={\tilde{x}}_{k+1}+K_{k+1}(y_{k+1}-C{\tilde{x}}_{k+1})\\ {\hat{p}}_{k+1}=(I-K_{k+1}C){\tilde{p}}_{k+1}\end{array}\right.$

式中：`Q<sub>k</sub>表示第k时刻系统噪声协方差，`R_k表示第k时刻量测噪声协方差，并且有

$\begin{array}{c}`Q_k=(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\left(-1\right)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}){(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\left(-1\right)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j})}^T\\ +(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\\ `R_k=\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;{\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;}^T+C{\tilde{p}}_k{C}^T\end{array}$

下面介绍本发明的一个实施例：

考虑非高斯Lévy噪声下分数阶线性离散系统模型

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}=A{\stackrel{\&OverBar;}{x}}_k+{\mathrm{Bu}}_k+{\stackrel{\&OverBar;}{w}}_k\\ {\stackrel{\&OverBar;}{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}\\ {\stackrel{\&OverBar;}{y}}_k={\mathrm{Cx}}_k+{\stackrel{\&OverBar;}{v}}_k\end{array}\right.$

式中

$A=\left[\begin{array}{cc}0& 1\\ -0.1& -0.2\end{array}\right],B=\left[\begin{array}{c}1\\ 1\end{array}\right],C=\left[\begin{array}{cc}0.1& 1\\ 0.6& 1\end{array}\right],\mathrm{\&alpha;}=\left[\begin{array}{c}0.7\\ 1.2\end{array}\right],\mathrm{\&delta;}=\left[\begin{array}{c}4\\ 5\end{array}\right]$

运用本发明所提到的方法进行仿真，结果如图2所示。

Claims (1)

1.一种处理Lévy噪声的分数阶卡尔曼滤波方法，其特征在于，包括以下步骤：

(1)、初始化，包括：设定状态预测量的初值和预测误差协方差的初值；

(2)、对当前时刻系统噪声和量测噪声进行处理，处理步骤为：

$`w_k=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_1\&CenterDot;sign\left({\stackrel{\&OverBar;}{w}}_k\right)& if& \left|{\stackrel{\&OverBar;}{w}}_k\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ {\stackrel{\&OverBar;}{w}}_k& if& \left|{\stackrel{\&OverBar;}{w}}_k\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

$`v_k=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_2\&CenterDot;sign\left({\stackrel{\&OverBar;}{v}}_k\right)& if& \left|{\stackrel{\&OverBar;}{v}}_k\right|\&GreaterEqual;{\mathrm{\&delta;}}_2\\ {\stackrel{\&OverBar;}{v}}_k& if& \left|{\stackrel{\&OverBar;}{v}}_k\right|<{\mathrm{\&delta;}}_2\end{array}\right.$

式中，`w<sub>k</sub>表示系统噪声的近似值，下标k表示第k时刻，δ<sub>1</sub>表示选取的第一个阈值，sign(x)表示符号函数，如果x＞0，返回1；如果x＜0，返回-1。表示含Lévy序列的系统噪声，|x|表示取x的绝对值，`v_k表示量测噪声的近似值，δ₂表示选取的第二个阈值，表示含Lévy序列的量测噪声；

(3)、利用系统噪声和量测噪声的近似值，带入原系统得到新的分数阶线性系统模型，新模型可表示为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^a{x}_{k+1}={\mathrm{Ax}}_k+{\mathrm{Bu}}_k+`w_k\\ x_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j}\\ y_k={\mathrm{Cx}}_k+`v_k\end{array}\right.$

式中，Δ表示分数阶算子，上标α表示分数阶阶次向量，x_k+1表示第k+1时刻状态向量，A、B和C表示适当维数的已知矩阵，u_k表示第k时刻控制输入向量，`w<sub>k</sub>表示第k时刻系统噪声的近似值，y<sub>k</sub>表示第k时刻量测向量，`v_k表示第k时刻量测噪声的近似值，γ_j表示

其中

α_N表示第N个分数阶阶次值，j＝1，2，…，k+1；

(4)、计算当前时刻量测噪声协方差，计算步骤为：

$`R_k=\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;{\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;}^T+C{\tilde{p}}_k{C}^T$

式中，`R_k表示第k时刻量测噪声协方差，y_k表示第k时刻量测向量，C表示适当维数的已知矩阵，上标T表示转置，表示第k时刻状态的预测量，表示第k时刻状态预测误差协方差；

(5)、计算当前时刻最优滤波增益，计算步骤为：

$K_k={\tilde{p}}_k{C}^T{(C{\tilde{p}}_k{C}^T+`R_k)}^{-1}$

式中，K_k表示第k时刻最优滤波增益；

(6)、计算当前时刻状态估计值，计算步骤为：

${\hat{x}}_k={\tilde{x}}_k+K_k\&lsqb;y_k-C{\tilde{x}}_k\&rsqb;$

(7)、计算当前时刻估计误差协方差，计算步骤为：

${\hat{p}}_k=(I-K_{k}C){\hat{p}}_k$

(8)、计算当前时刻系统噪声协方差，计算步骤为：

$\begin{array}{c}`Q_k=\left(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}\right){\left(x_{k+1}-A{\hat{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}\right)}^T\\ +(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=1}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

式中，`Q_k表示第k时刻系统噪声协方差，表示第k时刻状态估计向量，u_k表示控制输入向量，x_k+1表示第k+1时刻状态向量，表示第k+1-j时刻状态估计误差协方差；

(9)、利用当前时刻状态估计值更新下一时刻状态预测值，计算步骤为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^a{\tilde{x}}_{k+1}=A{\hat{x}}_k+{\mathrm{Bu}}_k\\ {\tilde{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\tilde{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\&Sigma;}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}\end{array}\right.$

式中Δ表示分数阶算子，上标α表示分数阶阶次向量；

(10)、利用当前时刻估计误差协方差更新下一时刻状态预测协方差，计算步骤为：

${\tilde{p}}_{k+1}=(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+`Q_k+\underset{j=2}{\overset{k+1}{\&Sigma;}}{\mathrm{\&gamma;}}_j{\hat{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T$

(11)、判断k+1是否大于等于步长L(可自己设定)，如果是，结束计算，否则返回步骤(2)进行下一次估计。