CN103825576B - 非线性系统的多项式滤波故障检测方法 - Google Patents

Abstract

本发明公开了一种非线性系统的多项式滤波故障检测方法，包括：多项式逼近步骤，将非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别用多项式表示，获得由μ阶多项式和μ+1阶多项式余项构成的多项式逼近模型；在无故障情况下，将所述多项式逼近模型中的μ+1阶多项式余项表示为系数为比例矩阵和不确定矩阵的积的低阶多项式；基于所述比例矩阵和不确定矩阵获得增广状态演化方程；滤波器设计步骤，针对所述增广状态演化方程构建滤波器函数获得第k+1时刻的系统状态变量的估计值；确定所述滤波增益系数以使得增广状态估计误差的方差的上界在第k+1时刻最小；故障检测步骤，确定第k时刻的检测阈值，根据所述检测阈值确定故障检测策略。

Description

非线性系统的多项式滤波故障检测方法

技术领域

本发明涉及故障检测领域，具体地说，涉及一种非线性系统的多项式滤波故障检测方法。

背景技术

由于非线性在真实系统中的广泛存在，非线性系统的滤波与控制问题得到大量的研究关注。如果系统中的非线性不能得到很好地处理，则可能导致震荡甚至发散。扩展Kalman滤波器（EKF）是最小方差准则下处理非线性估计问题的一种经典方法，目前也已经有部分研究希望可以对扩展Kalman滤波器（EKF）方法进行改进，以提升其处理非线性和随机性的能力。

多项式Kalman滤波器（PEKF）基于非线性函数的多项式逼近对扩展Kalman滤波器（EKF）进行了扩展。EKF方法只考虑非线性的线性化部分，而KalmanPEKF利用阶次为μ的多项式来估计非线性。PEKF处理由系统原始状态的Kronecker幂组成的扩展状态，目前已经得到了大量研究。但是由于高阶多项式余项的表达式非常复杂，目前的PEKF研究都忽略了高阶多项式余项，即忽略了多项式逼近误差，因此可能导致无法得到对非线性系统的精确估计结果。

随着现代系统运行安全性要求的提高，故障检测（FD）问题在实际应用中得到越来越多的重视。由于设计合适的滤波器可以有效的检测到系统中异常情况的发生，因此目前出现了大量基于估计的故障检测方法。在现有的多项式Kalman滤波器（PEKF）方法中，由于高阶多项式余项被忽略导致设计的滤波器不够合理，因此不能准确检测故障。

基于上述情况，亟需一种非线性系统的多项式滤波故障检测方法在考虑高阶多项式余项的情况下准确检测故障。

发明内容

本发明针对现有技术中存在的上述问题，提供了一种非线性系统的多项式滤波故障检测方法，包括：

多项式逼近步骤，将非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别用多项式表示，获得由μ阶多项式和μ+1阶多项式余项构成的多项式逼近模型；在无故障情况下，将所述多项式逼近模型中的μ+1阶多项式余项表示为系数为比例矩阵和不确定矩阵的积的低阶多项式，所述不确定矩阵的范数有界；基于所述比例矩阵和不确定矩阵获得增广状态演化方程；其中，1≤μ,1≤m≤μ，μ和m均为大于等于1的正整数，μ是人为设定的逼近阶次；

滤波器设计步骤，针对所述增广状态演化方程构建滤波器函数，所述滤波器函数根据增广状态下第k时刻的系统状态变量，残差信号，以及滤波增益系数获得第k+1时刻的系统状态变量的估计值；确定所述滤波增益系数以使得增广状态估计误差的方差的上界在第k+1时刻最小；

故障检测步骤，根据所述增广状态演化方程确定第k时刻的检测阈值，根据所述检测阈值确定故障检测策略。

根据本发明的一实施例，所述多项式逼近步骤中，基于所述比例矩阵和不确定矩阵获得增广状态演化方程包括：

在有故障的情况下，依据所述比例矩阵和不确定矩阵针对所述非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别构建系统演化方程；

定义由非线性系统的系统状态变量和测量输出变量的kronecker幂构成的增广状态，基于所述系统状态转移方程获得增广状态演化方程。

根据本发明的一实施例，所述多项式逼近步骤中，所述多项式逼近模型为：

$\left\{\begin{array}{c}{x}_{k+1}^{\left[m\right]}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{G}_{m,i}({\tilde{x}}_k,{\mathrm{\μ}}_k,v_k,f_k){(x_k-{\tilde{x}}_k)}^{\left[i\right]}+G_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,u_k,v_k,f_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}\\ y_k^{\left[m\right]}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{H}_{m,i}({\tilde{x}}_k,w_k){(x_k-{\tilde{x}}_k)}^{\left[i\right]}+H_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,w_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}\end{array}\right.$

其中，为第k+1时刻系统状态变量m阶kronecker幂，为第k时刻测量输出变量的m阶kronecker幂，θ∈[0,1]，x_k为系统状态变量，为系统状态变量x_k的估计值，u_k为控制输入，y_k为测量输出变量，f_k为故障信号，w_k为测量噪声，v_k为系统过程噪声，w_k和v_k为零均值且互相独立；

$G_{m,i}({\tilde{x}}_k,u_k,v_k,f_k)=\frac{1}{i!}({\▿}_x^{\left[i\right]}\⊗{(g+v+f)}^{\left[m\right]}){|}_{(x={\tilde{x}}_k,u=u_k,f=f_k)},$

$H_{m,i}({\tilde{x}}_k,w_k)=\frac{1}{i!}({\▿}_x^{\left[i\right]}\⊗{(h+w)}^{\left[m\right]}){|}_{(x={\tilde{x}}_k,w=w_k)},$

运算符定义为

${\▿}_x^{\left[0\right]}\⊗x=x,$

${\▿}_x^{[i+1]}\⊗x={\▿}_x\⊗{\▿}_x^{\left[i\right]}\⊗x,i>0,$

其中g表示系统函数，h表示测量函数。

根据本发明的一实施例，所述多项式逼近步骤中，在无故障情况下，将所述多项式逼近模型中的μ+1阶多项式余项表示为

$G_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,{\mathrm{\μ}}_k,v_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}=L_{g,m,k}{\mathrm{\Δ}}_{g,m,k}{(x_k-{\tilde{x}}_k)}^{\left[\mathrm{\μ}\right]}$

$H_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,u_k,v_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}=L_{h,m,k}{\mathrm{\Δ}}_{h,m,k}{(x_k-{\tilde{x}}_k)}^{\left[\mathrm{\μ}\right]}$

其中，L_g,m,k和L_h,m,k为比例矩阵，△_g,m,k和△_h,m,k为不确定矩阵，且有‖△_g,m,k‖≤1和‖△_h,m,k‖≤1。

根据本发明的一实施例，所述在有故障的情况下，依据所述比例矩阵和不确定矩阵针对所述非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别构建故障状态时的系统演化方程为：

其中，

$A_{m,j,k}=\underset{i=j}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(({\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]})\⊗{\mathrm{\ζ}}_{v,k,m-p})M_i^j(I_{n^j}\⊗{(-{\tilde{x}}_k)}^{[i-j]}),$

$C_{m,j,k}=\underset{i=j}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(({\▿}_x^{\left[i\right]}\⊗h^{\left[p\right]})\⊗{\mathrm{\ζ}}_{w,k,m-p})M_i^j(I_{n^j}\⊗{(-{\tilde{x}}_k)}^{[i-j]}),$

$v_{m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(\left(({\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]}){(x_k-{\tilde{x}}_k)}^{\left[i\right]}\right)\⊗(v_k^{[m-p]}-{\mathrm{\ζ}}_{v,k,m-p})),$

$w_{m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(\left(({\▿}_x^{\left[i\right]}\⊗p^{\left[p\right]}){(x_k-{\tilde{x}}_k)}^{\left[i\right]}\right)\⊗(w_k^{[m-p]}-{\mathrm{\ζ}}_{w,k,m-p})),$

其中，为矩阵因数，是维数为n^j的单位矩阵，

根据本发明的一实施例，所述增广状态定义为：

$x_k^{\left(\mathrm{\μ}\right)}=\left[\begin{array}{c}1\\ x_k\\ x_k^{\left[2\right]}\\ .\\ .\\ .\\ x_k^{\left[\mathrm{\μ}\right]}\end{array}\right],y_k^{\left(\mathrm{\μ}\right)}=\left[\begin{array}{c}1\\ y_k\\ y_k^{\left[2\right]}\\ .\\ .\\ .\\ y_k^{\left[\mathrm{\μ}\right]}\end{array}\right],$

所述增广状态演化方程由下式表示：

其中，为增广状态下第k+1时刻的系统状态变量，为增广状态下第k时刻的测量输出变量，为增广状态下的过程噪声，为增广状态下的测量噪声，

$v_k^{\left(\mathrm{\μ}\right)}={(v_{0,k}^T,v_{1,k}^T,...,v_{\mathrm{\μ},k}^T)}^T,w_k^{\left(\mathrm{\μ}\right)}={(w_{0,k}^T,w_{1,k}^T,...,w_{\mathrm{\μ},k}^T)}^T,$ Υ_k＝[Υ_0，k，Υ_1，k，...，Υ_μ，k]

L_g,k=diag{L_g,0,k,L_g,1,k,…,L_g,μ,k},L_h,k=diag{L_h,0,k,L_h，1,k,…,L_h,μ,k},

△_g,k=diag{△_g,0,k,△_g,1,k,…,△_g,μ,k},△_h,k=diag{△_h,0,k,△_h,1,k,…,△_h,μ,k}。

根据本发明的一实施例，所述滤波器设计步骤中，所述滤波器函数表示为：

$\begin{array}{c}{\tilde{x}}_{k+1}^{\left(\mathrm{\μ}\right)}=A_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)}+K_k{r}_k\\ =A_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)}+K_k(y_k^{\left(\mathrm{\μ}\right)}-C_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)})\end{array}$

其中，为第k时刻的残差信号，为增广状态下第k时刻的系统状态变量，为第k+1时刻的系统状态变量的估计值，K_k为滤波增益系数；

所述增广状态估计误差为

根据本发明的一实施例，所述故障检测步骤中，所述检测阈值为：

其中，

${\mathrm{\δ}}_{v,m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}\left|\right|M_m^p\left|\right|\left|\right|{\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]}\left|\right|{\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}(s_{v,k,m-p}+|\left|{\mathrm{\ζ}}_{v,k,m-p}\right),$

${\mathrm{\δ}}_{w,m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}\left|\right|M_m^p\left|\right|\left|\right|{\▿}_x^{\left[i\right]}\⊗h^{\left[p\right]}\left|\right|{\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}(s_{w,k,m-p}+|\left|{\mathrm{\ζ}}_{w,k,m-p}\right|\left|\right),$

${\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}=\underset{q=0}{\overset{i}{\mathrm{\Σ}}}\left|\right|M_i^q(I_{n^q}\⊗{(-{\tilde{x}}_k)}^{[i-q]})\left|\right|\left(\right|\left|{\tilde{x}}_{q,k}^{\left(\mathrm{\μ}\right)}\right||+{\mathrm{\δ}}_{e,k}),$

${\mathrm{\δ}}_{v,k}=\sqrt{\underset{m=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{\left({\mathrm{\δ}}_{v,m,k}\right)}^2},$

${\mathrm{\δ}}_{w,k}=\sqrt{\underset{m=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{\left({\mathrm{\δ}}_{w,m,k}\right)}^2},$

δ_e,k的初始值： ${\mathrm{\δ}}_{e,i,0}=s_{0,i}+\left|\right|{\mathrm{\ζ}}_{0,i}\left|\right|,{\mathrm{\δ}}_{e,0}=\sqrt{\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{\mathrm{\δ}}_{e,i,0}^2}.$

根据本发明的一实施例，所述故障检测步骤中的故障检测策略为：对于第k时刻，当所述残差信号的范数小于所述检测阈值时，系统正常；当所述残差信号的范数大于所述检测阈值时，系统故障。

本发明针对一类非线性系统，给出了一种新颖的状态估计与故障诊断的多项式方法。对非线性系统进行阶次为μ的多项式逼近，其中高阶多项式余项则利用系数为范数有界不确定矩阵的低阶多项式来刻画。在考虑高阶多项式余项的情况下，通过状态增广得到由原系统状态Kronecker幂组成的增广状态演化规律，针对增广之后的状态设计估计及故障检测器。滤波器参数设计问题转化为最小二乘滤波问题，通过求解一系列线性矩阵递推方程，可得到滤波器的参数设计方案，使得在无故障情况下，增广状态估计误差方差上界最小。通过分析估计状态和多项式系数，计算故障检测阈值，并设计对应的检测策略。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细的说明。

图1是根据本发明实施例一的非线性系统的多项式滤波故障检测方法流程图；

图2是本发明实施例二的非线性系统真实状态与估计状态的曲线图；

图3是本发明实施例二的非线性系统残差信号范数与故障检测阈值的曲线图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。

需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。另外，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

实施例一

以下参考图1来详细说明非线性系统的多项式滤波故障检测方法。包括：

步骤S101，多项式逼近步骤，将非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别用多项式表示，获得由μ阶多项式和μ+1阶多项式余项构成的多项式逼近模型；其中，1≤μ,1≤m≤μ，μ和m均为大于等于1的正整数，μ是人为设定的逼近阶次；

具体的，将离散非线性系统表示为：

$\left\{\begin{array}{c}{x}_{k+1}=g(x_k,u_k)+v_k+f_k,\\ y_k=h\left(x_k\right)+w_k,\end{array}\right.---\left(1\right)$

其中x_k为系统状态变量；u_k为控制输入；y_k为测量输出；f_k为故障信号；w_k为测量噪声，v_k为系统过程噪声，w_k和v_k为零均值且互相独立；g(·)表示系统函数，h(·)表示测量函数，g(·)和h(·)为μ+1阶连续可微。

对于任意向量a，将其m阶Kronecker积记作a^[m]。

系统初始状态x₀和噪声彼此独立，对于i=1,2,…,2μ均有，数学期望且范数有界：

$\left|\right|x_0^{\left[i\right]}\left|\right|\≤s_{0,i},\left|\right|v_k^{\left[i\right]}\left|\right|\≤s_{v,k,i},\left|\right|w_k^{\left[i\right]}\left|\right|\≤s_{w,k,i}.$

对于和利用Taylor展开式和中值定理将公式（1）转化为多项式逼近模型：

$\left\{\begin{array}{c}{x}_{k+1}^{\left[m\right]}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{G}_{m,i}({\tilde{x}}_k,{\mathrm{\μ}}_k,v_k,f_k){(x_k-{\tilde{x}}_k)}^{\left[i\right]}+G_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,u_k,v_k,f_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}\\ y_k^{\left[m\right]}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{H}_{m,i}({\tilde{x}}_k,w_k){(x_k-{\tilde{x}}_k)}^{\left[i\right]}+H_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,w_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}\end{array}\right.---\left(2\right)$

可以看出，公式（2）表示的多项式逼近模型由μ阶多项式和μ+1阶多项式余项构成。其中，为第k+1时刻系统状态变量m阶kronecker幂，为第k时刻测量输出变量的m阶kronecker幂，而θ∈[0,1]，x_k为系统状态变量，为系统状态变量x_k的估计值。公式（2）中，

$G_{m,i}({\tilde{x}}_k,u_k,v_k,f_k)=\frac{1}{i!}({\▿}_x^{\left[i\right]}\⊗{(g+v+f)}^{\left[m\right]}){|}_{(x={\tilde{x}}_k,u=u_k,f=f_k)},$

$H_{m,i}({\tilde{x}}_k,w_k)=\frac{1}{i!}({\▿}_x^{\left[i\right]}\⊗{(h+w)}^{\left[m\right]}){|}_{(x={\tilde{x}}_k,w=w_k)},$

运算符定义为

${\▿}_x^{\left[0\right]}\⊗x=x,$

${\▿}_x^{[i+1]}\⊗x={\▿}_x\⊗{\▿}_x^{\left[i\right]}\⊗x,i>0,$

其中 ${\▿}_x=[\∂/{\∂x}_1,...,\∂/2x_n].$

在无故障情况下，将所述多项式逼近模型中的μ+1阶多项式余项表示为系数为比例矩阵和不确定矩阵的积的低阶多项式，所述不确定矩阵的范数有界；

具体的，在离散非线性系统无故障的情况下，若系统连续，可得：

$G_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,u_k,v_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}=L_{g,m,k}{\mathrm{\Δ}}_{g,m,k}{(x_k-{\tilde{x}}_k)}^{\left[\mathrm{\μ}\right]}$

以及

$H_{m,\mathrm{\μ}+1}({{\tilde{x}}_{\mathrm{\θ}}}_k,u_k,v_k){(x_k-{\tilde{x}}_k)}^{[\mathrm{\μ}+1]}=L_{h,m,k}{\mathrm{\Δ}}_{h,m,k}{(x_k-{\tilde{x}}_k)}^{\left[\mathrm{\μ}\right]}$

其中L_g,m,k和L_h,m,k为比例矩阵，△_g,m,k和△_h,m,k为不确定矩阵，且有‖△_g,m,k‖≤1和‖△_h,m,k‖≤1。

基于所述比例矩阵L_g,m,k和L_h,m,k和不确定矩阵△_g,m,k和△_h,m,k获得增广状态演化方程，包括以下两个子步骤：

首先，在有故障的情况下，依据所述比例矩阵和不确定矩阵针对所述非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别构建系统演化方程；

具体的，可将公式(2)转化为系统演化方程，用下式表示：

其中

$A_{m,j,k}=\underset{i=j}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(({\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]})\⊗{\mathrm{\ζ}}_{v,k,m-p})M_i^j(I_{n^j}\⊗{(-{\tilde{x}}_k)}^{[i-j]}),$

$C_{m,j,k}=\underset{i=j}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(({\▿}_x^{\left[i\right]}\⊗h^{\left[p\right]})\⊗{\mathrm{\ζ}}_{w,k,m-p})M_i^j(I_{n^j}\⊗{(-{\tilde{x}}_k)}^{[i-j]}),$

$v_{m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(\left(({\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]}){(x_k-{\tilde{x}}_k)}^{\left[i\right]}\right)\⊗(v_k^{[m-p]}-{\mathrm{\ζ}}_{v,k,m-p})),$

$w_{m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}{M}_m^p(\left(({\▿}_x^{\left[i\right]}\⊗p^{\left[p\right]}){(x_k-{\tilde{x}}_k)}^{\left[i\right]}\right)\⊗(w_k^{[m-p]}-{\mathrm{\ζ}}_{w,k,m-p})),$

其次，定义由非线性系统的系统状态变量和测量输出变量的kronecker幂构成增广状态：

$x_k^{\left(\mathrm{\μ}\right)}=\left[\begin{array}{c}1\\ x_k\\ x_k^{\left[2\right]}\\ .\\ .\\ .\\ x_k^{\left[\mathrm{\μ}\right]}\end{array}\right],y_k^{\left(\mathrm{\μ}\right)}=\left[\begin{array}{c}1\\ y_k\\ y_k^{\left[2\right]}\\ .\\ .\\ .\\ y_k^{\left[\mathrm{\μ}\right]}\end{array}\right],$

根据公式（3），则可得增广状态演化方程，用下式表示：

其中，为增广状态下第k+1时刻的系统状态变量，为增广状态下第k时刻的测量输出变量，为增广状态下的过程噪声，为增广状态下的测量噪声，

$v_k^{\left(\mathrm{\μ}\right)}={(v_{0,k}^T,v_{1,k}^T,...,v_{\mathrm{\μ},k}^T)}^T,w_k^{\left(\mathrm{\μ}\right)}={(w_{0,k}^T,w_{1,k}^T,...,w_{\mathrm{\μ},k}^T)}^T,$ Υ_k＝[Υ_0，k，Υ_1，k，...，Υ_μ，k]，

L_g,k=diag{L_g,0,k,L_g,1,k,…,L_g,μ,k},L_h,k=diag{L_h,0,k,L_h,1,k,…,L_h,μ,k},

△_g,k=diag{△_g,0,k,△_g,1,k,…,△_g,μ,k},△_h,k=diag{△_h,0,k,△_h,1,k,…,△_h,μ,k}

；

综上，在步骤S101中，对公式（1）定义的非线性系统进行阶次为μ的多项式逼近。考虑系统中的故障为加性故障，当非线性环节g(·)和h(·)μ+1阶可微时，可以利用Taylor展开式以及中值定理在无故障情况下对非线性系统进行多项式逼近。对于其中的高阶多项式余项，可以表示为系数为范数有界不确定矩阵的低阶多项式。

步骤S102，滤波器设计步骤，针对所述增广状态演化方程（4）构建滤波器函数：

$\left\{\begin{array}{c}{\tilde{x}}_{k+1}^{\left(\mathrm{\μ}\right)}=A_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)}+K_k{r}_k\\ =A_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)}+K_k(y_k^{\left(\mathrm{\μ}\right)}-C_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)})\end{array}\right.---\left(5\right)$

其中为是对的估计，且有而K_k为滤波增益系数，在下述步骤中确定K_k的取值，为第k时刻的残差信号，为增广状态下第k时刻的系统状态变量，为增广状态下第k+1时刻的系统状态变量的估计值，

由（5）式看出，滤波器函数根据增广状态下第k时刻的系统状态变量残差信号，以及滤波增益系数K_k获得第k+1时刻的系统状态变量

以下步骤确定所述滤波增益系数以使得增广状态估计误差的方差的上界在第k+1时刻最小；

记增广状态估计误差和增广状态估计误差的方差

设γ_1k,γ_2k,γ_3k,ε均为正实数，若初始值为的如下递推矩阵方程

有正定解，且满足如下条件

其中，

$Y_k=(1+\mathrm{\ϵ})C_k{({\stackrel{\‾}{P}}_k^{-1}-{\mathrm{\γ}}_{1,k}{Y}_k^T{Y}_k)}^{-1}{C}_k^T+2({\mathrm{\γ}}_{1,k}^{-1}(1+\mathrm{\ϵ})+{\mathrm{\γ}}_{2,k}^{-1}(1+{\mathrm{\ϵ}}^{-1}))L_{h,k}{L}_{h,k}^T+{\stackrel{\‾}{\mathrm{\ψ}}}_k^w,$

$Z_k=(1+\mathrm{\ϵ})C_k{({\stackrel{\‾}{P}}_k^{-1}-{\mathrm{\γ}}_{1,k}{Y}_k^T{Y}_k)}^{-1}{A}_k^T,$

$\begin{array}{c}{\stackrel{\‾}{\mathrm{\ψ}}}_{a,b,k}^v=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{j=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{a}{\mathrm{\Σ}}}\underset{q=0}{\overset{b}{\mathrm{\Σ}}}\frac{1}{i!j!}{M}_a^p(\left(({\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]}){\stackrel{\‾}{\mathrm{\Ω}}}_{i,j,k}{({\▿}_x^{\left[j\right]}\⊗g^{\left[q\right]})}^T\right)\⊗({\mathrm{st}}_{n^{a-p},n^{b-q}}^{-1}({\mathrm{\ζ}}_{v,k,a-p+b-q})\\ -{\mathrm{\ζ}}_{v,k,a-p}{\mathrm{\ζ}}_{v,k,b-q}^T\left)\right){\left(M_b^q\right)}^T,\end{array}$

$\begin{array}{c}{\stackrel{\‾}{\mathrm{\Ψ}}}_{a,b,k}^w=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{j=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{a}{\mathrm{\Σ}}}\underset{q=0}{\overset{b}{\mathrm{\Σ}}}\frac{1}{i!j!}{M}_a^p(\left(({\▿}_x^{\left[i\right]}\⊗h^{\left[p\right]}){\stackrel{\‾}{\mathrm{\Ω}}}_{i,j,k}{({\▿}_x^{\left[j\right]}\⊗h^{\left[q\right]})}^T\right)\⊗({\mathrm{st}}_{n^{a-p},n^{b-q}}^{-1}({\mathrm{\ζ}}_{w,k,a-p+b-q})\\ -{\mathrm{\ζ}}_{w,k,a-p}{\mathrm{\ζ}}_{w,k,b-q}^T\left)\right){\left(M_b^q\right)}^T,\end{array}$

${\stackrel{\‾}{\mathrm{\Ω}}}_{i,j,k}=\underset{p=0}{\overset{i}{\mathrm{\Σ}}}\underset{q=0}{\overset{j}{\mathrm{\Σ}}}{M}_i^p\left((I_{n^p}\⊗{(-{\tilde{x}}_k)}^{[i-p]}){\stackrel{\‾}{X}}_{p,q,k}(I_{n^q}\⊗{(-{\tilde{x}}_k)}^{[j-q]})\right){\left(M_j^q\right)}^T,$

${\stackrel{\‾}{X}}_{k+1}=A_k{({\stackrel{\‾}{X}}_k^{-1}-{\mathrm{\γ}}_{3,k}{Y}_k^T{Y}_k)}^{-1}{A}_k^T+{\mathrm{\γ}}_{3,k}^{-1}{L}_{g,k}{L}_{g,k}^T+{\stackrel{\‾}{\mathrm{\Ψ}}}_k^v,$

其中，为矩阵因数，可利用Francesco等人1996年“离散时间非高斯线性系统多项式滤波”（Polynomial filtering for linear discrete time non-Gaussian systems）一文中的定理得到，是维数为nj的单位矩阵，

且有：

则由下式设计的滤波增益系数：

$K_k=Z_k^T{Y}_k^{-1}---\left(10\right)$

可保证对于任意时刻k，均有即增广状态估计误差的方差存在上界且公式(10)中的滤波增益系数K_k使得第k+1时刻的上界最小，下文给出证明过程；

至此为止，通过设计滤波增益系数K_k使得在每一时刻增广状态估计误差的方差的上界最小，实现了最小二乘意义下的滤波器设计；

综上，在步骤S102中，对于公式（1）表示的非线性系统转化成的含有不确定矩阵的多项式系统，利用最小二乘准则设计滤波器，用于估计系统状态并检测系统故障。通过计算得到增广状态估计误差方差的某一个上界，并通过滤波器增益设计，使得该上界在每一时刻达到最小。

步骤S103，故障检测步骤，根据所述增广状态演化方程（4）确定第k时刻的检测阈值：

其中，

${\mathrm{\δ}}_{v,m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}\left|\right|M_m^p\left|\right|\left|\right|{\▿}_x^{\left[i\right]}\⊗g^{\left[p\right]}\left|\right|{\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}(s_{v,k,m-p}+|\left|{\mathrm{\ζ}}_{v,k,m-p}\right),$

${\mathrm{\δ}}_{w,m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}\underset{p=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{i!}\left|\right|M_m^p\left|\right|\left|\right|{\▿}_x^{\left[i\right]}\⊗h^{\left[p\right]}\left|\right|{\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}(s_{w,k,m-p}+|\left|{\mathrm{\ζ}}_{w,k,m-p}\right|\left|\right),$

${\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}=\underset{q=0}{\overset{i}{\mathrm{\Σ}}}\left|\right|M_i^q(I_{n^q}\⊗{(-{\tilde{x}}_k)}^{[i-q]})\left|\right|\left(\right|\left|{\tilde{x}}_{q,k}^{\left(\mathrm{\μ}\right)}\right||+{\mathrm{\δ}}_{e,k}),$

${\mathrm{\δ}}_{v,k}=\sqrt{\underset{m=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{\left({\mathrm{\δ}}_{v,m,k}\right)}^2},$

${\mathrm{\δ}}_{w,k}=\sqrt{\underset{m=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{\left({\mathrm{\δ}}_{w,m,k}\right)}^2},$

δ_e,k的初始值： ${\mathrm{\δ}}_{e,i,0}=s_{0,i}+\left|\right|{\mathrm{\ζ}}_{0,i}\left|\right|,{\mathrm{\δ}}_{e,0}=\sqrt{\underset{i=0}{\overset{\mathrm{\μ}}{\mathrm{\Σ}}}{\mathrm{\δ}}_{e,i,0}^2}.$

根据残差信号定义如下的故障检测策略：

对于第k时刻，当所述残差信号的范数小于所述检测阈值时，即||r_k||≤δ_r,k时，系统正常；当所述残差信号的范数大于所述检测阈值时，即||r_k||>δ_r,k时，系统故障。

本实施例确定的检测阈值δ_r,k为在线更新的时变阈值，与定常阈值相比，时变检测阈值具有更好的跟踪和在线调整能力。

证明过程一

以下说明比例矩阵L_g，m，k和L_h，m，k的推导过程：

L_g,m,k可利用如下方式确定，记 ${\▿}_x^{\left[\mathrm{\μ}\right]}\⊗{(g+v)}^{\left[m\right]}={[{\mathrm{\ϵ}}_{1,\mathrm{\μ},m}(x,u,v),...,{\mathrm{\ϵ}}_{n^m,\mathrm{\μ},m}(x,u,v)]}^T,$ 其中对所有i=1,…,n^m，均有ε_i,μ,m(x,u,v):则可得：

其中

考虑如下情形：系统真实状态x_k与估计状态之间存在如下关系：其中E_k为分布矩阵，z_k为范数小于1的不确定矢量。由于函数对连续，因此一定存在极值使得 $\left|\right|{\▿}_x^T\⊗{\mathrm{\ϵ}}_{i,\mathrm{\μ},m}^T({\tilde{x}}_{{\mathrm{\θ}}_k},u_k,v_k)\left|\right|\≤{\stackrel{\‾}{m}}_{i,k},$ 且考虑到 $x_k={\tilde{x}}_k+E_k{z}_k,$ 则有 $\left|\right|{(x_k-{\tilde{x}}_k)}^T({\▿}_x^T\⊗{\mathrm{\ϵ}}_{i,\mathrm{\μ},m}^T({\tilde{x}}_{{\mathrm{\θ}}_k},u_k,v_k))\left|\right|\≤{\stackrel{\‾}{m}}_{i,k}\left|\right|E_k\left|\right|.$ 因此存在向量使得‖ζ_i,k‖≤1，并且

${(x_k-{\tilde{x}}_k)}^T({\▿}_x^T\⊗{\mathrm{\ϵ}}_{i,\mathrm{\μ},m}^T({\tilde{x}}_{{\mathrm{\θ}}_k},u_k,v_k))={\stackrel{\‾}{m}}_{i,k}\left|\right|E_k\left|\right|{\mathrm{\ζ}}_{i,k}^T.;$

记 ${\mathrm{\Ξ}}_k={[{\mathrm{\ζ}}_{1,k},...,{\mathrm{\ζ}}_{n^m,k}]}^T,$ 则有

由Ξ_k定义可知，因此可得

$L_{g,m,k}=\frac{\sqrt{n^m}}{(\mathrm{\μ}+1)!}\left|\right|E_k\left|\right|\mathrm{diag}\{{\stackrel{\‾}{m}}_{1,k},...,{\stackrel{\‾}{m}}_{n^m,k}\}.$

L_h,m,k可利用类似方法获得。

证明过程二

以下证明当滤波增益系数时，增广状态估计误差方差上界最小。

首先可证明 和

由矩阵不等式相关知识可得：

$\begin{array}{c}{P}_{k+1}\≤(1+\mathrm{\ϵ})((A_k-K_k{C}_k){({\stackrel{\‾}{P}}_k^{-1}-{\mathrm{\γ}}_{1,k}{Y}_k^T{Y}_k)}^{-1}{(A_k-K_k{C}_k)}^T+{2\mathrm{\γ}}_{1,k}^{-1}(L_{g,k}{L}_{g,k}^T+K_k{L}_{h,k}{L}_{h,k}^T{K}_k^T))\\ +2{\mathrm{\γ}}_{2,k}^{-1}(1+{\mathrm{\ϵ}}^{-1})(L_{g,k}{L}_{g,k}^T+K_k{L}_{h,k}{L}_{h,k}^T{K}_k^T)+{\stackrel{\‾}{\mathrm{\Ψ}}}_k^v+K_k{\stackrel{\‾}{\mathrm{\Ψ}}}_k^w{K}_k^T\end{array}$

考虑到Y_k和Z_k的定义，可得

$\begin{array}{c}{P}_{k+1}\≤K_k{Y}_k{K}_k^T-K_k{Z}_k-Z_k^T{K}_k^T+(1+\mathrm{\ϵ})A_k{({\stackrel{\‾}{P}}_k^{-1}-{\mathrm{\γ}}_{1,k}{Y}_k^T{Y}_k)}^{-1}{A}_k^T\\ +2({\mathrm{\γ}}_{1,k}^{-1}(1+\mathrm{\ϵ})+{\mathrm{\γ}}_{2,k}^{-1}\left({1+\mathrm{\ϵ}}^{-1}\right))L_{g,k}{L}_{g,k}^T+{\stackrel{\‾}{\mathrm{\Ψ}}}_k^v\end{array}$

注意到对上式的K_k进行配方，可得

$\begin{array}{c}{P}_{k+1}\≤(K_k-Z_k^T{Y}_k^{-1})Y_k{(K_k-Z_k^T{Y}_k^{-1})}^T-Z_k^T{Y}_k^{-1}{Z}_k+(1+\mathrm{\ϵ})A_k{({\stackrel{\‾}{P}}_k^{-1}-{\mathrm{\γ}}_{1,k}{Y}_k^T{Y}_k)}^{-1}{A}_k^T\\ +2({\mathrm{\γ}}_{1,k}^{-1}(1+\mathrm{\ϵ})+{\mathrm{\γ}}_{2,k}^{-1}(1+{\mathrm{\ϵ}}^{-1}))L_{g,k}{L}_{g,k}^T+{\stackrel{\‾}{\mathrm{\Ψ}}}_k^v.\end{array}$

因此，增广状态估计误差的方差存在上界，显然当时，增广状态估计误差方差上界最小。

实施例二

本实施例提供一种离散非线性系统：

$\left\{\begin{array}{c}{x}_{k+1}=0.1x_k\sin x_k-0.8x_k+u_k+v_k+f_k,\\ y_k=0.1x_k+w_k,\\ u_k=y_k.\end{array}\right.$

其中过程噪声满足：

$\left\{\begin{array}{c}\mathrm{Pr}(v_k=-1\×{10}^{-3})=0.6,\\ \mathrm{Pr}(v_k=0)=0.2,\\ \mathrm{Pr}(v_k=3\×{10}^{-3})=0.2,\end{array}\right.$

测量噪声满足：

$\left\{\begin{array}{c}\mathrm{Pr}(w_k=-7\×{10}^{-3})=0.3,\\ \mathrm{Pr}(w_k=3\×{10}^{-3})=0.7.\end{array}\right.$

其中Pr(·)表示某个随机事件发生的概率。

从第k=20时刻开始加入幅度为0.4的故障，如下形式

$f_k=\left\{\begin{array}{cc}0.4,& k>20,\\ 0,& k\≤20.\end{array}\right.$

在本实施例的离散非线性系统中，使用实施例一中提供的故障检测方法，得到图2和图3所示的结果。

图2所示为真实状态与估计状态的曲线图，当在第k=20时刻之前时，估计状态和真实状态匹配，不会产生故障信号；当第k=20时刻出现故障时，估计状态与真实状态之间出现较大差距，该差距可显示在图3所示的残差信号的范数中。图3所示为残差信号范数与故障检测阈值的曲线图。由图3可以看出，当第k=20时刻出现故障时，残差信号的范数大于故障检测阈值，可以及时检测到故障信号的发生。同时，时变检测阈值的设置可以准确判断出故障信号的发生。

综上所述，本发明针对一类非线性系统，给出了一种新颖的状态估计与故障诊断的多项式方法。对非线性系统进行阶次为μ的多项式逼近，其中高阶多项式余项则利用系数为范数有界不确定矩阵的低阶多项式来刻画。通过状态增广，得到由原系统状态Kronecker幂组成的增广状态演化规律，针对增广之后的状态设计估计及故障检测器。滤波器参数设计问题转化为最小二乘滤波问题，通过求解一系列线性矩阵递推方程，可得到滤波器的参数设计方案，使得在无故障情况下，增广状态估计误差方差上界最小。通过分析估计状态和多项式系数，计算故障检测阈值，并设计对应的检测策略。

对于基于估计的故障检测方法，残差评价在故障检测中起着非常重要的作用，包含阈值确定和检测策略两部分。与定常阈值相比，本发明中利用测量信号在线更新的时变阈值具有更好的跟踪和在线调整能力。

在多项式方法中，高阶余项表达式形势极为复杂，本发明在考虑高阶多项式余项的情况下计算了估计误差方差的上界，从而设计出更加合理的滤波器，结合时变检测阈值精确检测故障。

虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (3)

1.一种非线性系统的多项式滤波故障检测方法，其特征在于：包括：

多项式逼近步骤，将非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别用多项式表示，获得由μ阶多项式和μ+1阶多项式余项构成的多项式逼近模型；在无故障情况下，将所述多项式逼近模型中的μ+1阶多项式余项表示为系数为比例矩阵和不确定矩阵的积的低阶多项式，所述不确定矩阵的范数有界；基于所述比例矩阵和不确定矩阵获得增广状态演化方程；其中，1≤μ,1≤m≤μ，μ和m均为大于等于1的正整数，μ是人为设定的逼近阶次；

滤波器设计步骤，针对所述增广状态演化方程构建滤波器函数，所述滤波器函数根据增广状态下第k时刻的系统状态变量，残差信号，以及滤波增益系数获得第k+1时刻的系统状态变量的估计值；确定所述滤波增益系数以使得增广状态估计误差的方差的上界在第k+1时刻最小；

故障检测步骤，根据所述增广状态演化方程确定第k时刻的检测阈值，根据所述检测阈值确定故障检测策略，

所述多项式逼近步骤中，所述多项式逼近模型为：

$\left\{\begin{array}{c}{x}_{k+1}^{\[m\]}=\underset{i=0}{\overset{\mathrm{\μ}}{\Σ}}{G}_{m,i}({\tilde{x}}_k,u_k,v_k,f_k){(x_k-{\tilde{x}}_k)}^{\[i\]}+G_{m,\mathrm{\μ}+1}({\tilde{x}}_{{\mathrm{\θ}}_k},u_k,v_k,f_k){(x_k-{\tilde{x}}_k)}^{\[\mathrm{\μ}+1\]}\\ y_k^{\[m\]}=\underset{i=0}{\overset{\mathrm{\μ}}{\Σ}}{H}_{m,i}({\tilde{x}}_k,w_k){(x_k-{\tilde{x}}_k)}^{\[i\]}+H_{m,\mathrm{\μ}+1}({\tilde{x}}_{{\mathrm{\θ}}_k},w_k){(x_k-{\tilde{x}}_k)}^{\[\mathrm{\μ}+1\]}\end{array}\right.$

其中，为第k+1时刻系统状态变量m阶kronecker幂，为第k时刻测量输出变量的m阶kronecker幂，θ∈[0,1]，x_k为系统状态变量，为系统状态变量x_k的估计值，u_k为控制输入，y_k为测量输出变量，f_k为故障信号，w_k为测量噪声，v_k为系统过程噪声，w_k和v_k为零均值且互相独立；

$G_{m,i}({\tilde{x}}_k,u_k,v_k,f_k)=\frac{1}{i!}({\▿}_x^{\[i\]}\⊗{\left(g+v+f\right)}^{\[m\]}){|}_{(x={\tilde{x}}_k,u=u_k,v=v_k,f=f_k)},$

$H_{m,i}({\tilde{x}}_k,w_k)=\frac{1}{i!}({\▿}_x^{\[i\]}\⊗{(h+w)}^{\[m\]}){|}_{(x={\tilde{x}}_k,w=w_k)},$

运算符定义为：

${\▿}_x^{\[0\]}\⊗\mathrm{\χ}=\mathrm{\χ},$

${\▿}_x^{\[i+1\]}\⊗\mathrm{\χ}={\▿}_x\⊗{\▿}_x^{\[i\]}\⊗\mathrm{\χ},i>0,$

其中，g表示系统函数，h表示测量函数，

所述多项式逼近步骤中，在无故障情况下，将所述多项式逼近模型中的μ+1阶多项式余项表示为：

$G_{m,\mathrm{\μ}+1}({\tilde{x}}_{{\mathrm{\θ}}_k},u_k,v_k){(x_k-{\tilde{x}}_k)}^{\[\mathrm{\μ}+1\]}=L_{g,m,k}{\mathrm{\Δ}}_{g,m,k}{(x_k-{\tilde{x}}_k)}^{\[\mathrm{\μ}\]}$

$H_{m,\mathrm{\μ}+1}({\tilde{x}}_{{\mathrm{\θ}}_k},u_k,v_k){(x_k-{\tilde{x}}_k)}^{\[\mathrm{\μ}+1\]}=L_{h,m,k}{\mathrm{\Δ}}_{h,m,k}{(x_k-{\tilde{x}}_k)}^{\[\mathrm{\μ}\]}$

其中，L_g,m,k和L_h,m,k为比例矩阵，Δ_g,m,k和Δ_h,m,k为不确定矩阵，且有||Δ_g,m,k||≤1和||Δ_h,m,k||≤1，

所述多项式逼近步骤中，基于所述比例矩阵和不确定矩阵获得增广状态演化方程包括：

在有故障的情况下，依据所述比例矩阵和不确定矩阵针对所述非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别构建系统演化方程；

定义由非线性系统的系统状态变量和测量输出变量的kronecker幂构成的增广状态，基于所述系统状态转移方程获得增广状态演化方程，

所述在有故障的情况下，依据所述比例矩阵和不确定矩阵针对所述非线性系统的系统状态变量和测量输出变量的m阶kronecker幂分别构建故障状态时的系统演化方程为：

其中，

$A_{m,j,k}=\underset{i=j}{\overset{\mathrm{\μ}}{\Σ}}\underset{p=0}{\overset{m}{\Σ}}\frac{1}{i!}{M}_m^p\left(\left({\▿}_x^{\[i\]}\⊗g^{\[p\]}\right)\⊗{\mathrm{\ξ}}_{v,k,m-p}\right)M_i^j\left(I_{n^j}\⊗{(-{\tilde{x}}_k)}^{\[i-j\]}\right),$

$C_{m,j,k}=\underset{i=j}{\overset{\mathrm{\μ}}{\Σ}}\underset{p=0}{\overset{m}{\Σ}}\frac{1}{i!}{M}_m^p\left(\left({\▿}_x^{\[i\]}\⊗h^{\[p\]}\right)\⊗{\mathrm{\ξ}}_{w,k,m-p}\right)M_i^j\left(I_{n^j}\⊗{(-{\tilde{x}}_k)}^{\[i-j\]}\right),$

$v_{m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\Σ}}\underset{p=0}{\overset{m}{\Σ}}\frac{1}{i!}{M}_m^p\left(\right(\left({\▿}_x^{\[i\]}\⊗g^{\[p\]}\right){\left(x_k-{\tilde{x}}_k\right)}^{\[i\]})\⊗(v_k^{\[m-p\]}-{\mathrm{\ξ}}_{v,k,m-p}\left)\right),$

$w_{m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\Σ}}\underset{p=0}{\overset{m}{\Σ}}\frac{1}{i!}{M}_m^p\left(\right(\left({\▿}_x^{\[i\]}\⊗h^{\[p\]}\right){\left(x_k-{\tilde{x}}_k\right)}^{\[i\]})\⊗(w_k^{\[m-p\]}-{\mathrm{\ξ}}_{w,k,m-p}\left)\right),$

其中，为矩阵因数，是维数为n^j的单位矩阵，

所述增广状态定义为：

$x_k^{\left(\mathrm{\μ}\right)}=\left[\begin{array}{c}1\\ x_k\\ x_k^{\[2\]}\\ .\\ .\\ .\\ x_k^{\[\mathrm{\μ}\]}\end{array}\right],y_k^{\left(\mathrm{\μ}\right)}=\left[\begin{array}{c}1\\ y_k\\ y_k^{\[2\]}\\ .\\ .\\ .\\ y_k^{\[\mathrm{\μ}\]}\end{array}\right],$

所述增广状态演化方程由下式表示：

其中，为增广状态下第k+1时刻的系统状态变量，为增广状态下第k时刻的测量输出变量，为增广状态下的过程噪声，为增广状态下的测量噪声，

L_g,k＝diag{L_g,0,k,L_g,1,k,...,L_g,μ,k},L_h,k＝diag{L_h,0,k,L_h,1,k,...,L_h,μ,k},

Δ_g,k＝diag{Δ_g,0,k,Δ_g,1,k,...,Δ_g,μ,k},Δ_h,k＝diag{Δ_h,0,k,Δ_h,1,k,...,Δ_h,μ,k}

所述滤波器设计步骤中，所述滤波器函数表示为：

$\begin{array}{c}{\tilde{x}}_{k+1}^{\left(\mathrm{\μ}\right)}=A_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)}+K_k{r}_k\\ =A_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)}+K_k(y_k^{\left(\mathrm{\μ}\right)}-C_k{\tilde{x}}_k^{\left(\mathrm{\μ}\right)})\end{array}$

其中，为第k时刻的残差信号，为增广状态下第k时刻的系统状态变量，为第k+1时刻的系统状态变量的估计值，K_k为滤波增益系数；

所述增广状态估计误差为

2.如权利要求1所述方法，其特征在于，所述故障检测步骤中，所述检测阈值为：

其中，

${\mathrm{\δ}}_{v,m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\Σ}}\underset{p=0}{\overset{m}{\Σ}}\frac{1}{i!}\left|\right|M_m^p\left|\right|\left|\right|{\▿}_x^{\[i\]}\⊗g^{\[p\]}\left|\right|{\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}(s_{v,k,m-p}+|\left|{\mathrm{\ξ}}_{v,k,m-p}\right|\left|\right),$

${\mathrm{\δ}}_{w,m,k}=\underset{i=0}{\overset{\mathrm{\μ}}{\Σ}}\underset{p=0}{\overset{m}{\Σ}}\frac{1}{i!}\left|\right|M_m^p\left|\right|\left|\right|{\▿}_x^{\[i\]}\⊗h^{\[p\]}\left|\right|{\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}(s_{w,k,m-p}+|\left|{\mathrm{\ξ}}_{w,k,m-p}\right|\left|\right),$

${\stackrel{\‾}{\mathrm{\ρ}}}_{i,k}=\underset{q=0}{\overset{i}{\Σ}}\left|\right|M_i^q(I_{n^q}\⊗{\left(-{\tilde{x}}_k\right)}^{\[i-q\]})\left|\right|\left(\right|\left|{\tilde{x}}_{q,k}^{\left(\mathrm{\μ}\right)}\right||+{\mathrm{\δ}}_{e,k}),$

${\mathrm{\δ}}_{v,k}=\sqrt{\underset{m=0}{\overset{\mathrm{\μ}}{\Σ}}{\left({\mathrm{\δ}}_{v,m,k}\right)}^2},$

${\mathrm{\δ}}_{w,k}=\sqrt{\underset{m=0}{\overset{\mathrm{\μ}}{\Σ}}{\left({\mathrm{\δ}}_{w,m,k}\right)}^2},$

δ_e,k的初始值：

3.如权利要求2所述的方法，其特征在于，所述故障检测步骤中的故障检测策略为：对于第k时刻，当所述残差信号的范数小于所述检测阈值时，系统正常；当所述残差信号的范数大于所述检测阈值时，系统故障。