CN104462015A - 处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法 - Google Patents

Abstract

本发明公开了一种处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法。首先，给出状态预测初始值和预测误差协方差初始值。接着，利用近似方法分解得到非高斯Lévy噪声的分解值，推导状态近似值和量测输出近似值，并由此计算相应系统噪声协方差和量测噪声协方差。然后，利用当前状态估计值计算下一时刻状态预测值，利用当前时刻的估计误差协方差和系统噪声协方差计算该下一时刻预测误差协方差。最后，结合所得到的状态预测值更新状态估计值，利用预测误差协方差更新估计误差协方差。本发明由于解决了非高斯Lévy噪声下分数阶线性离散系统的状态估计问题，从而拓展了分数阶理论的应用范围，而且本发明易于与已有的状态估计软件相结合。

Description

处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法

技术领域

本发明涉及一种处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法，属于系统分析与处理技术领域。

背景技术

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

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

发明内容

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

由于非高斯噪声普遍存在于实际系统中，那么利用传统的卡尔曼滤波方法对这些系统进行状态估计时不能得出理想结果。处理非高斯Lévy噪声的状态更新方法通过近似处理得到状态量和量测量的近似值，并由此推导系统噪声和量测噪声协方差，进而对得到状态预测值和预测协方差值，最后对状态估计值和估计误差协方差进行更新。该方法与常规动态状态估计程序相结合，能很好地处理非高斯噪声情况下的状态估计问题。

技术方案：一种处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法，所述方法是在计算机中依次按以下步骤实现的：

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

(2)、对非高斯Lévy噪声进行近似处理，推导得出状态近似值和量测输出近似值，计算步骤为：

${\stackrel{\&OverBar;}{x}}_{k+1}=\left\{\begin{array}{cc}{\mathrm{\&Omega;}}_1+{\mathrm{\&delta;}}_1\&CenterDot;\mathrm{sign}\left({\mathrm{\&Omega;}}_2\right)& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ x_{k+1}& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

${\stackrel{\&OverBar;}{y}}_{k+1}=\left\{\begin{array}{cc}C{\tilde{x}}_{k+1}+{\mathrm{\&delta;}}_2\&CenterDot;\mathrm{sign}(y_{k+1}-C{\tilde{x}}_{k+1})& \mathrm{abs}(y_{k+1}-C{\tilde{x}}_{k+1})\&GreaterEqual;{\mathrm{\&delta;}}_2\\ y_{k+1}& \mathrm{abs}(y_{k+1}-C{\tilde{x}}_{k+1})<{\mathrm{\&delta;}}_2\end{array}\right.$

式中

${\mathrm{\&Omega;}}_1=A{\tilde{x}}_k+{\mathrm{Bu}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}$

${\mathrm{\&Omega;}}_2=x_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}$

(3)利用已得到的状态近似值和量测输出近似值，计算系统噪声协方差和量测噪声协方差，计算步骤为：

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k\\ {+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})}^T+(A+{\mathrm{\&gamma;}}_1){\tilde{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\tilde{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

${\stackrel{\&OverBar;}{R}}_{k+1}=({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1}){({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})}^T+C{\tilde{p}}_{k+1}{C}^T$

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

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\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}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}\end{array}\right.$

(5)、利用当前时刻估计误差协方差和系统噪声协方差，计算下一时刻预测误差协方差，计算步骤为：

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

(6)、利用预测误差协方差和量测噪声协方差计算滤波增益矩阵，计算步骤为：

$K_{k+1}={\tilde{p}}_{k+1}{C}^T{(C{\tilde{p}}_{k+1}{C}^T+{\stackrel{\&OverBar;}{R}}_{k+1})}^{-1}$

(7)、利用状态预测值，已得到的滤波增益矩阵和量测输出近似值，更新状态估计值，计算步骤为：

${\hat{x}}_{k+1}={\tilde{x}}_{k+1}+K_{k+1}({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})$

(8)、利用滤波增益矩阵和预测误差协方差更新估计误差协方差，计算步骤为：

${\hat{p}}_{k+1}=(I-K_{k+1}C){\tilde{p}}_{k+1}$

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

传统的卡尔曼滤波(Kalman Filter，KF)算法处理的是整数阶系统，并且要求系统噪声和量测噪声均为高斯白噪声，这意味着系统的阶次必须是整数，噪声的概率分布是正态分布，而且满足它的二阶矩不相关。这些要求只是理想化的定义，在现实生活中很难得以满足，比如一些粘弹性结构、有耗网络和扩散波等，这些本身就有着分数阶的性质，用整数阶模型很难对这些系统进行建模。另一方面，实际系统中的噪声大多为非高斯白噪声，这样的系统处理起来更具有一定的复杂性和挑战性。

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

附图说明

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

图2为本发明提出的状态估计方法得到的仿真值，其中：图2(a)是变量1的真实值与估计值，图2(b)是变量2的真实值与估计值。

具体实施方式

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

如图1所示，处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法，包括如下步骤：

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

(2)、对非高斯Lévy噪声进行近似处理，推导得出状态近似值和量测输出近似值。

(3)利用已得到的状态近似值和量测输出近似值，计算系统噪声协方差和量测噪声协方差。

(4)、利用当前时刻状态估计值，计算下一时刻状态预测值。

(5)、利用当前时刻估计误差协方差和系统噪声协方差，计算下一时刻预测误差协方差。

(6)、利用预测误差协方差和量测噪声协方差计算滤波增益矩阵。

(7)、利用状态预测值，已得到的滤波增益矩阵和量测输出近似值，更新状态估计值。

(8)、利用滤波增益矩阵和预测误差协方差更新估计误差协方差。

(9)、判断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\left[{\hat{e}}_{k+1}{\hat{e}}_{k+1}^T\right]$

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

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

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

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

${\mathrm{\&gamma;}}_j=\mathrm{diag}\left(\right[\left(\begin{array}{c}{\mathrm{\&alpha;}}_1\\ j\end{array}\right)\left(\begin{array}{c}{\mathrm{\&alpha;}}_2\\ j\end{array}\right)...\left(\begin{array}{c}{\mathrm{\&alpha;}}_N\\ j\end{array}\right)\left]\right),$ 其中 $\left(\begin{array}{c}\mathrm{\&alpha;}\\ j\end{array}\right)=\left\{\begin{array}{cc}1& j=0\\ \frac{\mathrm{\&alpha;}(\mathrm{\&alpha;}-1)...(\mathrm{\&alpha;}-j+1)}{j!}& j>0\end{array}\right.$

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

在以上所示分数阶线性离散系统模型中，系统噪声和量测噪声均为非高斯Lévy噪声，由于其具有无限方差特性，因此不能用来直接计算协方差矩阵，必须用近似处理的方法得出具有高斯白噪声性质的近似值才可以用来计算协方差矩阵。近似处理后状态向量和量测输出向量可以表示为：

${\stackrel{\&OverBar;}{x}}_{k+1}=\left\{\begin{array}{cc}{\mathrm{\&Omega;}}_1+{\mathrm{\&delta;}}_1\&CenterDot;\mathrm{sign}\left({\mathrm{\&Omega;}}_2\right)& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ x_{k+1}& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

${\stackrel{\&OverBar;}{y}}_{k+1}=\left\{\begin{array}{cc}C{\tilde{x}}_{k+1}+{\mathrm{\&delta;}}_2\&CenterDot;\mathrm{sign}(y_{k+1}-C{\tilde{x}}_{k+1})& \mathrm{abs}(y_{k+1}-C{\tilde{x}}_{k+1})\&GreaterEqual;{\mathrm{\&delta;}}_2\\ y_{k+1}& \mathrm{abs}(y_{k+1}-C{\tilde{x}}_{k+1})<{\mathrm{\&delta;}}_2\end{array}\right.$

式中

${\mathrm{\&Omega;}}_1=A{\tilde{x}}_k+{\mathrm{Bu}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}$

${\mathrm{\&Omega;}}_2=x_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}$

δ₁表示所取的第一个域值，δ₂表示所取的第二个域值，表示k+1时刻状态预测值，sign(x)为符号函数，如果x>0，返回1；如果x<0，返回-1。abs表示取绝对值函数。

系统噪声协方差阵和量测噪声协方差阵可分别表示为：

${\stackrel{\&OverBar;}{Q}}_k=E\left[{\stackrel{\&OverBar;}{w}}_k{\stackrel{\&OverBar;}{w}}_k^T\right]$

${\stackrel{\&OverBar;}{R}}_{k+1}=E\left[{\stackrel{\&OverBar;}{v}}_{k+1}{\stackrel{\&OverBar;}{v}}_{k+1}^T\right]$

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

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

将上述两式分别代入系统和量测噪声协方差阵可得：

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k\\ {+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})}^T+(A+{\mathrm{\&gamma;}}_1){\tilde{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\tilde{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

${\stackrel{\&OverBar;}{R}}_{k+1}=({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1}){({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})}^T+C{\tilde{p}}_{k+1}{C}^T$

式中表示第k时刻状态误差协方差。

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

${\hat{x}}_k=E\left[x_k\right|y_k^{*}]$

${\hat{P}}_k=\mathrm{cov}\left[x_k\right|y_k^{*}]$

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

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

${\tilde{x}}_{k+1}=E\left[x_{k+1}\right|y_k^{*}]$

${\tilde{P}}_{k+1}=\mathrm{cov}\left[x_{k+1}\right|y_k^{*}]$

根据以上公式得到k时刻到k+1时刻的状态估计值及估计误差协方差分别为：

${\hat{x}}_{k+1}=E\left[x_{k+1}\right|y_{k+1}^{*}]$

${\hat{P}}_{k+1}=\mathrm{cov}\left[x_{k+1}\right|y_{k+1}^{*}]$

状态量估计步骤总结如下：

步骤1：在当前控制量和量测累计值基础上，利用当前k时刻的状态估计量计算k+1时刻的状态预测量；

步骤2：利用当前k时刻的估计误差协方差计算k+1时刻的预测误差协方差；

步骤3：利用k+1时刻的状态预测量更新状态估计量；

步骤4：利用k+1时刻的预测误差协方差更新估计误差协方差。

根据以上步骤，k+1时刻状态量的预测值为：

$\begin{array}{c}{\tilde{x}}_{k+1}=E\left[x_{k+1}\right|y_k^{*}]\\ =E\left[({\mathrm{Ax}}_k+{\mathrm{Bu}}_k+{\stackrel{\&OverBar;}{w}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j})\right|y_k^{*}]\\ =\mathrm{AE}\left[x_k\right|y_k^{*}]+{\mathrm{Bu}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_{j}E[x_{k+1-j}\left|y_k^{*}\right]\end{array}$

上式中最后一项可近似处理为：

$E\left[x_{k+1-j}\right|y_k^{*}]\&cong;[x_{k+1-j}\left|y_{k+1-j}^{*}\right]$

则k+1时刻状态量的预测值为：

${\tilde{x}}_{k+1}=A{\hat{x}}_k+{\mathrm{Bu}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}$

k+1时刻的预测误差协方差为：

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

式中表示第k时刻状态估计误差协方差。

接下来利用k+1时刻的状态预测量更新状态估计量，计算步骤为：

$\begin{array}{c}{\tilde{x}}_{k+1}=E\left[x_{k+1}\right|y_{k+1}^{*}]\\ =E\left[({\mathrm{Ax}}_k+{\mathrm{Bu}}_k+{\stackrel{\&OverBar;}{w}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j})\right|y_{k+1}^{*}]\\ ={\tilde{x}}_{k+1}+K_{k+1}({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})\end{array}$

其中，K_k+1为卡尔曼滤波增益，可表示为：

$K_{k+1}={\tilde{p}}_{k+1}{C}^T{(C{\tilde{p}}_{k+1}{C}^T+{\stackrel{\&OverBar;}{R}}_{k+1})}^{-1}$

最后，得到估计误差协方差更新值为：

$\begin{array}{c}{\tilde{P}}_{k+1}=\mathrm{cov}\left[x_{k+1}\right|y_{k+1}^{*}]\\ =\mathrm{cov}\left[({\mathrm{Ax}}_k+{\mathrm{Bu}}_k+{\stackrel{\&OverBar;}{w}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{x}_{k+1-j})\right|y_{k+1}^{*}]\\ =(I-K_{k+1}C){\tilde{p}}_{k+1}\end{array}$

处理非高斯Lévy噪声的分数阶卡尔曼滤波算法的完整的状态估计计算公式如下：

a.预测步：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\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}{\stackrel{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\hat{x}}_{k+1-j}\\ {\tilde{p}}_{k+1}=(A+{\mathrm{\&gamma;}}_1){\hat{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+{\stackrel{\&OverBar;}{Q}}_k+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\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+{\stackrel{\&OverBar;}{R}}_{k+1})}^{-1}\\ {\hat{x}}_{k+1}={\tilde{x}}_{k+1}+K_{k+1}({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})\\ {\hat{p}}_{k+1}=(I-K_{k+1}C){\tilde{p}}_{k+1}\end{array}\right.$

式中：为系统噪声协方差阵，为量测噪声协方差阵，并且有

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-{\mathrm{Bu}}_k\\ {+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})}^T+(A+{\mathrm{\&gamma;}}_1){\tilde{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\tilde{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}.$

${\stackrel{\&OverBar;}{R}}_{k+1}=({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1}){({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})}^T+C{\tilde{p}}_{k+1}{C}^T$

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

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

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}=A{x}_k+{\mathrm{Bu}}_k+{\stackrel{\&OverBar;}{w}}_k\\ {\stackrel{\&OverBar;}{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{x}_{k+1}-\underset{j=1}{\stackrel{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{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}0\\ 1\end{array}\right],$ C＝[0.1 0.3]， $\mathrm{\&alpha;}=\left[\begin{array}{c}0.7\\ 1.2\end{array}\right],\mathrm{\&delta;}=\left[\begin{array}{c}20\\ 20\end{array}\right]$

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

Claims (1)

1.一种处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法，其特征在于，包括以下步骤：

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

(2)、对非高斯Lévy噪声进行近似处理，推导得出状态近似值和量测输出近似值，计算步骤为：

${\stackrel{\&OverBar;}{x}}_{k+1}=\left\{\begin{array}{cc}{\mathrm{\&Omega;}}_1+{\mathrm{\&delta;}}_1\&CenterDot;\mathrm{sign}\left({\mathrm{\&Omega;}}_2\right)& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ x_{k+1}& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

${\stackrel{\&OverBar;}{y}}_{k+1}=\left\{\begin{array}{cc}C{\tilde{x}}_{k+1}+{\mathrm{\&delta;}}_2\&CenterDot;\mathrm{sign}(y_{k+1}-C{\tilde{x}}_{k+1})& \mathrm{abs}(y_{k+1}-C{\tilde{x}}_{k+1})\&GreaterEqual;{\mathrm{\&delta;}}_2\\ y_{k+1}& \mathrm{abs}(y_{k+1}-C{\tilde{x}}_{k+1})<{\mathrm{\&delta;}}_2\end{array}\right.$

式中

${\mathrm{\&Omega;}}_1=A{\tilde{x}}_k+B{u}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}$

${\mathrm{\&Omega;}}_2=x_{k+1}-A{\tilde{x}}_k+B{u}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}$

(3)利用已得到的状态近似值和量测输出近似值，计算系统噪声协方差和量测噪声协方差，计算步骤为：

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-B{u}_k+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j})({\stackrel{\&OverBar;}{x}}_{k+1}-A{\tilde{x}}_k-B{u}_k\\ +\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\tilde{x}}_{k+1-j}{)}^T+(A+{\mathrm{\&gamma;}}_1){\tilde{p}}_k{(A+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\tilde{p}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

${\stackrel{\&OverBar;}{R}}_{k+1}=({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1}){({\stackrel{\&OverBar;}{y}}_{k+1}-C{\stackrel{\&OverBar;}{x}}_{k+1})}^T+C{\tilde{p}}_{k+1}{C}^T$

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

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

(5)、利用当前时刻估计误差协方差和系统噪声协方差，计算下一时刻预测误差协方差，计算步骤为：

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

(6)、利用预测误差协方差和量测噪声协方差计算滤波增益矩阵，计算步骤为：

$K_{k+1}={\tilde{p}}_{k+1}{C}^T{(C{\tilde{p}}_{k+1}{C}^T+{\stackrel{\&OverBar;}{R}}_{k+1})}^{-1}$

(7)、利用状态预测值，已得到的滤波增益矩阵和量测输出近似值，更新状态估计值，计算步骤为：

${\hat{x}}_{k+1}={\tilde{x}}_{k+1}+K_{k+1}({\stackrel{\&OverBar;}{y}}_{k+1}-C{\tilde{x}}_{k+1})$

(8)、利用滤波增益矩阵和预测误差协方差更新估计误差协方差，计算步骤为：

${\hat{p}}_{k+1}=(I-K_{k+1}C){\tilde{p}}_{k+1}$

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