[发明专利]一种处理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|\&lt;{\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|\&lt;{\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<sub>k</sub>表示第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<sub>k</sub>表示第k时刻系统噪声协方差，表示第k时刻状态估计向量，u<sub>k</sub>表示控制输入向量，x<sub>k+1</sub>表示第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)进行下一次估计。