CN103714255A - 一种基于非线性故障重构的故障预测方法 - Google Patents

Abstract

本发明涉及一种基于非线性故障重构的故障预测方法，其包括如下步骤：1）利用KPCA模型对旋转机械系统运行时监测的过程数据进行离线非线性建模，并进行异常检测，提取故障信息，其中KPCA模型为核主元分析模型；2）通过步骤1）中的KPCA模型和异常检测实现对故障程度进行定量描述，并采用最优化方法求解故障重构下的故障估计；3）用多层递阶的方法对估计出的故障幅值f的发展趋势完成预测。本发明能有效提高故障预报效率，减少误报和漏报率，为进一步研究复杂机械系统的预测维护技术奠定基础。本发明可以广泛应用在石化、治金、煤炭等诸多企业的大型机电设备的在线监控系统中。

Description

一种基于非线性故障重构的故障预测方法

技术领域

本发明涉及一种故障预测方法，特别是关于一种机械故障检测技术领域中的基于非线性故障重构的故障预测方法。

背景技术

随着科学技术和现代工业的飞速发展，国民经济的机械、能源、石化、运载和国防等行业的装备日趋大型化、高速化、集成化和自动化，这类复杂系统一旦发生事故，便会造成巨大的财产损失和人员伤亡。因此，防范因部件性能退化导致的安全问题已成为复杂系统日常运行所要面临的重要难题。故障预测作为预测维护的重要环节，现已得到国内外专家的广泛关注。

在复杂工业过程环境下，传统的基于过程机理模型的故障诊断技术难以发挥作用，而数据驱动的过程监控技术在近年来得到了较好的发展。在数据驱动的方法中，研究内容和应用案例数量最多的是基于统计的方法。多元统计过程监控技术在诸多不同的工业过程控制领域得到了成功的应用，包括化工过程，高分子聚合物，微电子制造等等。多元统计过程监控主要由多元投影方法和多元控制图构成。其中多元投影方法的基本思想是将被监测的样本向量从高维的变量空间投影到由较少潜变量所张成的潜空间中去。多元控制图则是基于这种空间分解的结构，分别构造出反映所在空间变化的统计量，用以监测过程的异常。但是传统的统计模型是线性的，而在实际过程中，数据之间的关系一般是非线性的。如果这种非线性特性较强时，用线性模型进行过程监控容易导致大量误报和漏报。

对于一般的线性故障过程，可以采用自回归建模或是指数平滑模型，因为它们易于计算和使用。但前提都是假定故障过程可以被直接观测，但是实际上这一点是很难满足的。现有技术中采用将故障强度建模成一个离散的随机系统，虽然不能被测量，却能影响被测变量,即隐马尔科夫模型(HMM)；还有采用半隐马尔科夫模型（HSMM），将状态停留时间也考虑进去。尽管基于HMM和HSMM的方法能够处理隐含故障过程，但是它只能考虑单个部件的剩余有效寿命预测。对于一个大型的工业过程，对整个过程或设备的健康情况进行预测是更为重要的问题。由此可知，现有复杂系统的隐含性能退化故障预测问题的主要难点有以下两个方面：（1）多变量预测问题。目前人们主要利用单变量信号进行故障诊断，采用单一传感器所提供的信息进行研究有很大的局限性。采用多变量的故障预测方法才能进行整个系统级别的故障预测。（2）非线性问题。传统的将非线性因素忽略掉的结果往往并不理想，随着对故障检测和预测要求日益提高，原可忽略的非线性问题越来越突出。

发明内容

针对上述问题，本发明的目的是提供一种基于非线性故障重构的故障预测方法，该故障预测方法能提高故障预报效率，减少误报和漏报率，为进一步研究复杂机械系统的预测维护技术奠定基础。

为实现上述目的，本发明采取以下技术方案：一种基于非线性故障重构的故障预测方法，其包括如下步骤：1）利用KPCA模型对旋转机械系统运行时监测的过程数据进行离线非线性建模，并进行异常检测，提取故障信息，其中KPCA模型为核主元分析模型；2）通过步骤1）中的KPCA模型和异常检测实现对故障程度进行定量描述，并采用最优化方法求解故障重构下的故障估计；所述最优化方法求解故障重构下的故障估计包括以下步骤：（1）确定故障方向的提取：若故障方向未知，需要从已有的历史故障数据中提取故障方向，令X_f=[x₁,x₂,...,x_nf]^T代表故障方向Ξ下的故障数据，Ξ∈R^m×d；并采用滑动平均滤波法消除非线性故障数据中的正常部分

对非线性故障数据x_i的影响，滤波后故障数据矩阵

为：

${\underset{\‾}{X}}_f^T=\mathrm{\Ξ}[{\underset{\‾}{f}}_1,{\underset{\‾}{f}}_2,...,{\underset{\‾}{f}}_{\mathrm{nf}}],$

式中，Ξ为故障方向，f _nf为滤波后的故障幅值；并对故障数据矩阵进行奇异值分解，得到：

式中，D为对角矩阵，U为左奇异矩阵，V为右奇异矩阵，将对角阵D的非零奇异值按照降序排列，选故障方向Ξ＝U(d)，d为故障维数；（2）根据步骤（1）中的故障方向，进行非线性故障重构及估计：①将SPE指标作为重构检测指标，根据SPE指标进行故障幅值的大小的确定，则第i个样本的故障幅值f_i为：

f_i＝argminSPE(φ(x_i-Ξf_i))，i=1，2，…，n_f，

式中，argmin表示对目标函数最小化所对应的参数，φ表示样本在特征空间的映射向量；②取故障幅值f最小值，故障重构就是沿故障方向减去纯故障因素的部分，构造出数据的正常因素部分，将重构后的正常数据假设为z，则

z=x-Ξf；   （1）

将SPE指标对f^T求偏导，在0点处获得极值，得到故障幅值f最小值，根据步骤①得到SPE指标为：

$\begin{array}{c}{\mathrm{SPE}}_i=\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)-\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T{P}_{\stackrel{\‾}{\mathrm{\φ}}}{P_{\stackrel{\‾}{\mathrm{\φ}}}}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)\\ =\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)-\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T{\stackrel{\‾}{\mathrm{\Φ}}}^T{P}_{\stackrel{\‾}{k}}{M}^{-1}{P_{\stackrel{\‾}{k}}}^T\stackrel{\‾}{\mathrm{\Φ}}\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right),\\ =\stackrel{\‾}{k}(z_i,z_i)-\stackrel{\‾}{k}{\left(z_i\right)}^T{M}_s\stackrel{\‾}{k}\left(z_i\right)\end{array}---\left(2\right)$

式中，i=1，2，…，n_f， $M^{-1}=\mathrm{diag}(\frac{1}{{\mathrm{\μ}}_1},\frac{1}{{\mathrm{\μ}}_2},...,\frac{1}{{\mathrm{\μ}}_p}),P_{\stackrel{\‾}{k}}=[{\mathrm{\β}}_1,{\mathrm{\β}}_2,...,{\mathrm{\β}}_p],M_s=P_{\stackrel{\‾}{k}}{M}^{-1}{P_{\stackrel{\‾}{k}}}^T,$ Μ_s是对称矩阵，μ_p为

的特征值，β_p为与特征值μ_p对应的经过施密特标准正交化后的特征向量，

为均值中心化处理后的核矩阵；由式（1）得到SPE指标对f^T求偏导，进而得到故障大小f^T为：

$f^T=\frac{[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})]B\left(z\right)\mathrm{\Ξ}}{[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})]k\left(z\right)},$

式中，k(z)为核函数， $B\left(z\right)=\left[\begin{array}{c}k(x_1,z){(x-x_1)}^T\\ k(x_2,z){(x-x_2)}^T\\ \·\\ \·\\ \·\\ k(x_n,z){(x-x_n)}^T\end{array}\right];$

将f^T转化为f的形式： $f=\frac{{\mathrm{\Ξ}}^T{B\left(z\right)}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}\left(z\right)]}{k{\left(z\right)}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}\left(z\right)]},---\left(3\right)$

联立式（1）和式（3），得到递推公式如下，

$f^{(j+1)}=\frac{{\mathrm{\Ξ}}^{T}B{(x-{\mathrm{\Ξf}}^{\left(j\right)})}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}(x-{\mathrm{\Ξf}}^{\left(j\right)})]}{k{(x-{\mathrm{\Ξf}}^{\left(j\right)})}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}(x-{\mathrm{\Ξf}}^{\left(j\right)})]};$

3）用多层递阶的方法对估计出的故障幅值f的发展趋势完成预测。

所述步骤1）中，异常检测故障提取过程中，采用衡量过程数据主要变化的Hotelling T²统计量和衡量过程残差部分的SPE指标两类检测指标进行过程检测，其中SPE指标为平方预测误差指标，Hotelling T²统计量、SPE指标分别计算如下：

$T^2=k{\left(x\right)}^T{\mathrm{P\Λ}}^{-2}{P}^{T}k\left(x\right)\≤T_{\mathrm{\α}}^2,$

$\mathrm{SPE}=k(x,x)-k{\left(x\right)}^T{\mathrm{P\Λ}}^{-1}{P}^{T}k\left(x\right)\≤{\mathrm{\δ}}_{\mathrm{\α}}^2,$

式中，k(x)^T为核矩阵的转置矩阵；P为负载矩阵；Λ为协方差矩阵；

为Hotelling T²统计量的控制限值；k(x,x)为样本x在高维特征空间的内积范数；

为SPE指标的控制限值；当两种检测指标分别超过了各自的控制限值时，则表明过程数据中出现了异常；若在控制限值内，则表明过程数据正常。

本发明由于采取以上技术方案，其具有以下优点：本发明由于采用KPCA模型对旋转机械系统运行时监测的过程数据进行离线非线性建模，在非线性数据模型中，定义故障幅值并估计，进而通过故障幅值进行故障预测，因此能有效提高故障预报效率，减少误报和漏报率，为进一步研究复杂机械系统的预测维护技术奠定基础。本发明可以广泛应用在石化、治金、煤炭等诸多企业的大型机电设备的在线监控系统中。

附图说明

图1是基于KPCA的故障检测结果示意图；

图2是基于PCA模型的故障检测结果示意图；

图3是KPCA模型下，故障重构前后的SPE示意图；

图4（a）是KPCA模型下，故障重构估计出来的故障幅值f示意图；

图4（b）是KPCA模型下，故障重构估计出来的故障幅值绝对值|f|示意图；

图5是基于多层递阶方法的KPCA故障幅值的趋势预测。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提供一种基于非线性故障重构的故障预测方法，其包括以下步骤：

1）利用核主元分析（KPCA）模型对旋转机械系统运行时监测的过程数据进行离线非线性建模，并进行异常检测，提取故障信息；

假设x₁，x₂，…，x_n∈R^m是供核主元分析学习的n个m维的列向量样本，非线性映射为φ，原始空间R通过φ映射到高维特征空间F^high中去。原始数据x_i在高维特征空间F^high中的像为φ(x_i)，记Φ＝[φ(x₁),φ(x₂),...,φ(x_n)]^T。KPCA建模过程如下：

C_φv＝λv/n，   （1）

其中，C_φ＝Φ^TΦ/n为高维特征空间样本数据φ(x_i)的协方差矩阵，λ为特征值，v是与特征值λ对应的特征向量。

定义核矩阵K＝ΦΦ^T，特征向量α＝Φ^Tv，则有：

Kα＝λα，   （2）

其中，核矩阵K=（k_ij）中的元素按如下定义，

k_ij＝k(x_i,x_j)＝φ(x_i)^Tφ(x_j)，   （3）

根据公式（3），对高维特征空间样本数据φ(x_i)进行均值中心化处理后得到：

$\left\{\begin{array}{c}\stackrel{\‾}{\mathrm{\φ}}\left(x_i\right)=\mathrm{\φ}\left(x_i\right)-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\φ}\left(x_i\right)\\ \stackrel{\‾}{k}\left(x_{i,}{x}_j\right)=\stackrel{\‾}{\mathrm{\φ}}{\left(x_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(x_j\right)\end{array}\right.,$

函数k(x_i,x_j)即为核函数，函数

为均值中心化处理后的核函数，本申请采用的核函数为径向基函数，则有：

$k(x_i,x_j)=e^{-\frac{{\left|\right|x_i-x_j\left|\right|}^2}{c}},---\left(4\right)$

其中，c是径向基函数的带宽。

对式（2）求解得到特征值和特征向量（λ_i,α_i），根据式（4）可得：

v_i=λ^-1Φ^Tα_i   （5）

为了保证v_i是单位向量，需要对特征向量α_i进行归一化处理：

${\mathrm{\α}}_i\sqrt{{\mathrm{\λ}}_i}{\mathrm{\α}}_i^o,---\left(6\right)$

其中

表示式（2）中得到的单位长度的特征向量。令P_f=[v₁,...,v_l]为KPCA的负载矩阵，其中l为主元个数，由累计方差贡献率确定。对于一个新样本x，可以计算其非线性主元t如下：

$t=P_f^T\mathrm{\φ}\left(x\right)={\mathrm{\Λ}}^{-1}{P}^{T}k\left(x\right)---\left(7\right)$

其中，核矩阵k(x)＝[k(x₁,x),...,k(x_m,x)]^T，负载矩阵P=[α₁,...,α_l]，协方差矩阵Λ＝diag{λ₁,...,λ_l}。由此可知，式（7）中的变量都是可以计算的量，因此可以用于监测。

在基于核主元分析模型的多变量统计过程监控当中，采用衡量过程数据主要变化的Hotelling T²统计量和衡量过程残差部分的平方预测误差（SPE，SquaredPrediction Error）指标两类检测指标进行过程检测，实现异常检测故障提取。在KPCA模型下，两种检测指标都可以把故障点成功检测出来，相对于线性PCA而言，KPCA更适合非线性相关的过程监控（如图1、图2所示）。在KPCA模型中，Hotelling T²统计量、SPE指标可以分别计算如下：

$T^2=k{\left(x\right)}^T{\mathrm{P\Λ}}^{-1}{P}^{T}k\left(x\right)\≤T_{\mathrm{\α}}^2,---\left(8\right)$

$\mathrm{SPE}=k(x,x)-k{\left(x\right)}^T{\mathrm{P\Λ}}^{-1}{P}^{T}k\left(x\right)\≤{\mathrm{\δ}}_{\mathrm{\α}}^2,---\left(9\right)$

式中，k(x)^T为核矩阵的转置矩阵；P为负载矩阵；Λ为协方差矩阵；

为HotellingT²统计量的控制限值；k(x,x)为样本x在高维特征空间的内积范数；

为SPE指标的控制限值。

当两种检测指标分别超过了各自的控制限值时，则表明过程数据中出现了异常；若在控制限值内，则表明过程数据正常。

2）通过步骤1）中的KPCA模型和异常检测实现对故障程度进行定量描述，并采用最优化方法求解故障重构下的故障估计；其中，故障重构就是重新构造去除了故障信息的过程数据，使故障因素对该过程数据内正常部分的影响减至最小，从而构造出正常状态下的数据；故障估计就是故障重构成功后，估计出故障的大小（如图3、图4所示）。

最优化方法求解故障重构下的故障估计包括以下步骤：

（1）确定故障方向的提取。若故障方向未知，需要从已有的历史故障数据中提取故障方向。令X_f=[x₁,x₂,...,x_nf]^T代表故障方向Ξ下的故障数据，Ξ∈R^m×d。由于正常数据

是零均值的，因此可以采用滑动平均滤波法消除其对非线性故障数据x_i的影响，滤波后故障数据矩阵

为：

${\underset{\‾}{X}}_f^T=\mathrm{\Ξ}[{\underset{\‾}{f}}_1,{\underset{\‾}{f}}_2,...,{\underset{\‾}{f}}_{\mathrm{nf}}],---\left(10\right)$

式中，Ξ为故障方向，f _nf为滤波后的故障幅值。

对故障数据矩阵进行奇异值分解，得到：

$X_f^T={\mathrm{UDV}}^T,---\left(11\right)$

式中，D为对角矩阵，U为左奇异矩阵，V为右奇异矩阵。将对角阵D的非零奇异值按照降序排列，选故障方向Ξ＝U(d)，d为故障维数。实际应用中，故障维数d的取值是使得故障样本的重构检测指标处于正常范围的最小值。

（2）根据步骤（1）中的故障方向，进行非线性故障重构及估计：

①将SPE指标作为重构检测指标，根据SPE指标进行故障幅值的大小的确定，则第i个样本的故障幅值f_i为：

f_i＝argminSPE(φ(x_i-Ξf_i))，（i=1，2，…，n_f），   （12）

式中，argmin表示对目标函数最小化所对应的参数，φ表示样本在特征空间的映射向量。

尽管上式中的重构检测指标也可以用Hotelling T²来替代，但是由于某些最优算法（如牛顿最优算法）不能保证Hotelling T²的收敛性，所以故障诊断时一般不采用此检测指标。

②取故障幅值f最小值，故障重构就是沿故障方向减去纯故障因素的部分，构造出数据的正常因素部分。将重构后的正常数据假设为z，则：

z=x-Ξf。   （13）

将SPE指标对f^T求偏导，在0点处获得极值，得到故障幅值f最小值。根据式（9）对SPE指标进一步化简为如下：

$\begin{array}{c}{\mathrm{SPE}}_i=\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)-\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T{P}_{\stackrel{\‾}{\mathrm{\φ}}}{P_{\stackrel{\‾}{\mathrm{\φ}}}}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)\\ =\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)-\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T{\stackrel{\‾}{\mathrm{\Φ}}}^T{P}_{\stackrel{\‾}{k}}{M}^{-1}{P_{\stackrel{\‾}{k}}}^T\stackrel{\‾}{\mathrm{\Φ}}\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right),\\ =\stackrel{\‾}{k}(z_i,z_i)-\stackrel{\‾}{k}{\left(z_i\right)}^T{M}_s\stackrel{\‾}{k}\left(z_i\right)\end{array}---\left(14\right)$

在式（14）中，i=1，2，…，n_f， $M^{-1}=\mathrm{diag}(\frac{1}{{\mathrm{\μ}}_1},\frac{1}{{\mathrm{\μ}}_2},...,\frac{1}{{\mathrm{\μ}}_p}),$ 负载矩阵 $P_{\stackrel{\‾}{k}}=[{\mathrm{\β}}_1,{\mathrm{\β}}_2,...,{\mathrm{\β}}_p]$ 对称矩阵

μ_p为

的特征值，β_p为与特征值μ_p对应的经过施密特标准正交化后的特征向量，

为均值中心化处理后的核矩阵。

则由式（14）可得SPE指标对f^T求偏导为：

$\frac{\∂\mathrm{SPE}}{{\∂f}^T}=\frac{\∂\stackrel{\‾}{k}(z,z)}{{\∂f}^T}-\frac{\∂\left[\stackrel{\‾}{k}{\left(z\right)}^T{M}_s\stackrel{\‾}{k}\left(z\right)\right]}{{\∂f}^T}.---\left(15\right)$

根据公式（3）得到 $\stackrel{\‾}{k}(x_i,x_j)=k(x_i,x_j)-k{\left(x_i\right)}^T{1}_n-k{\left(x_j\right)}^T{1}_n+{1_n}^T{K1}_n,$ 则式（15）等号右边的第一部分可写为：

$\begin{array}{c}\frac{\∂\stackrel{\‾}{k}(z,z)}{{\∂f}^T}=\frac{\∂k(z,z)}{{\∂f}^T}-{{21}_n}^T\frac{\∂k\left(z\right)}{{\∂f}^T}+\frac{\∂\left[{1_n}^T{K1}_n\right]}{{\∂f}^T}\\ =-{{21}_n}^T\frac{\∂k\left(z\right)}{{\∂f}^T}\end{array},---\left(16\right)$

由于本申请选择的核函数是径向基函数

则式（16）第一个等号右边第一项实际上为0。

式（16）中向量对故障幅值求偏导，可以展成如下形式

$\frac{\∂k\left(z\right)}{{\∂f}^T}=\frac{\∂k\left(z\right)}{{\∂z}^T}\frac{\∂z}{{\∂f}^T}---\left(17\right)$

由核矩阵k(z)的定义知，k(z)＝[k(x₁,z),k(x₂,z),...,k(x_n,z)]^T。由于向量核矩阵k(z)中的任意一个元素可以写成

因此其对求偏导：

$\begin{array}{c}\frac{\∂k(x_i,z)}{{\∂z}^T}=e^{-\frac{{\left|\right|x_i-z\left|\right|}^2}{c}}\×[-\frac{1}{c}\frac{\∂({x_i}^T{x}_i-{{2x}_i}^{T}z+z^{T}z)}{{\∂z}^T}\\ =-\frac{1}{c}k(x_i,z)[-{{2x}_i}^T\frac{\∂z}{{\∂z}^T}+\frac{\∂\left(z^{T}z\right)}{{\∂z}^T},(i=\mathrm{1,2},...,n)\\ =\frac{2}{c}k(x_i,z)({{-x}_i}^T+z^T)\end{array}---\left(18\right)$

而 $\frac{\∂z}{{\∂f}^T}=[\frac{\∂z}{{\∂f}_1},\frac{\∂z}{{\∂f}_2},...,\frac{\∂z}{{\∂f}_d}]=-[{\mathrm{\ξ}}_1,{\mathrm{\ξ}}_2,...,{\mathrm{\ξ}}_d]=-\mathrm{\Ξ},$ 因此式（17）可以写成：

$\begin{array}{c}\frac{\∂k\left(z\right)}{{\∂f}^T}=[\frac{\∂k(x_1,z)}{{\∂z}^T},\frac{\∂k(x_2,z)}{{\∂z}^T},...,\frac{\∂k(x_n,z)}{{\∂f}^T}{]}^T\frac{\∂z}{{\∂f}^T}\\ =-\left[\begin{array}{c}\frac{2}{c}k(x_1,z)(-{x_1}^T+z^T)\\ \frac{2}{c}k(x_2,z)(-{x_2}^T+z^T)\\ \·\\ \·\\ \·\\ \frac{2}{c}k(x_n,z)(-{x_n}^T+z^T)\end{array}\right]\mathrm{\Ξ}\\ =-\frac{2}{c}\left[\begin{array}{c}k(x_1,z){(x-x_1)}^T+k(x_1,z)f^T{\mathrm{\Ξ}}^T\\ k(x_2,z){(x-x_2)}^T+k(x_2,z)f^T{\mathrm{\Ξ}}^T\\ \·\\ \·\\ \·\\ k(x_n,z){(x-x_n)}^T+k(x_n,z)f^T{\mathrm{\Ξ}}^T\end{array}\right]\mathrm{\Ξ}\end{array},---\left(19\right)$

若假设

$B\left(z\right)=\left[\begin{array}{c}k(x_1,z){(x-x_1)}^T\\ k(x_2,z){(x-x_2)}^T\\ \·\\ \·\\ \·\\ k(x_n,z){(x-x_n)}^T\end{array}\right],$

则式（19）可以写为：

$\frac{\∂k\left(z\right)}{{\∂f}^T}=-\frac{2}{c}[B\left(z\right)\mathrm{\Ξ}-k\left(z\right)f^T],---\left(20\right)$

则式（15）等号右边第二项可化为

$\begin{array}{c}\frac{\∂\left[\stackrel{\‾}{k}{\left(z\right)}^T{M}_s\stackrel{\‾}{k}\left(z\right)\right]}{{\∂f}^T}=\frac{\∂\stackrel{\‾}{k}{\left(z\right)}^T}{{\∂f}^T}{M}_s\stackrel{\‾}{k}\left(z\right)+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s\frac{\∂\stackrel{\‾}{k}\left(z\right)}{{\∂f}^T}\\ =2\stackrel{\‾}{k}{\left(z\right)}^T{M}_s\frac{\∂\stackrel{\‾}{k}\left(z\right)}{{\∂f}^T}\end{array},---\left(21\right)$

由式（21）可得，

$\frac{\stackrel{\‾}{k}\left(z\right)}{{\∂f}^T}=(I-1_{n\×n})\frac{\∂k\left(z\right)}{{\∂f}^T},$ 其中，I为单位矩阵。   （22）

把（16）、（20）、（21）、（22）代入（15），可得，

$\begin{array}{c}\frac{\∂\mathrm{SPE}}{{\∂f}^T}-{{21}_n}^T\frac{\∂k\left(z\right)}{{\∂f}^T}-2\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})\frac{\∂k\left(z\right)}{{\∂f}^T}\\ =\frac{4}{c}[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})][B\left(z\right)\mathrm{\Ξ}-k\left(z\right)f^T]\end{array},---\left(23\right)$

令（23）等于零向量，则可求得故障大小f^T为：

$f^T=\frac{[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})]B\left(z\right)\mathrm{\Ξ}}{[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})]k\left(z\right)},---\left(24\right)$

由此可知，式（24）中右边的分母其实是一个标量，而非矩阵。当此标量不为0时，则能达到可重构的条件；若为0，则故障幅值不可被重构。

为求取方便，将f^T转化为f的形式：

$f=\frac{{\mathrm{\Ξ}}^{T}B{\left(z\right)}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}\left(z\right)]}{k{\left(z\right)}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}\left(z\right)]}.---\left(25\right)$

由式（25）可知，故障幅值f无法直接求取，需要通过递推的方式确定。联立式（13）和式（25），可得递推公式如下，

$f^{(j+1)}=\frac{{\mathrm{\Ξ}}^{T}B{(x-{\mathrm{\Ξf}}^{\left(j\right)})}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}(x-{\mathrm{\Ξf}}^{\left(j\right)})]}{k{(x-{\mathrm{\Ξf}}^{\left(j\right)})}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}(x-{\mathrm{\Ξf}}^{\left(j\right)})]},---\left(26\right)$

其中上标括号内的j代表递推的次数。

3）最后用多层递阶的方法对估计出的故障幅值f的发展趋势完成预测，由于多层递阶方法把系统的时变特性考虑了进去，因此预报精度有所提升（如图5所示）。

上述各实施例仅用于说明本发明，其中各部件的结构和连接方式等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (2)

1.一种基于非线性故障重构的故障预测方法，其包括如下步骤：

1）利用KPCA模型对旋转机械系统运行时监测的过程数据进行离线非线性建模，并进行异常检测，提取故障信息，其中KPCA模型为核主元分析模型；

2）通过步骤1）中的KPCA模型和异常检测实现对故障程度进行定量描述，并采用最优化方法求解故障重构下的故障估计；所述最优化方法求解故障重构下的故障估计包括以下步骤：

（1）确定故障方向的提取：若故障方向未知，需要从已有的历史故障数据中提取故障方向，令X_f=[x₁,x₂,...,x_nf]^T代表故障方向Ξ下的故障数据，Ξ∈R^m×d；并采用滑动平均滤波法消除非线性故障数据中的正常部分

对非线性故障数据x_i的影响，滤波后故障数据矩阵

为：

${\underset{\‾}{X}}_f^T=\mathrm{\Ξ}[{\underset{\‾}{f}}_1,{\underset{\‾}{f}}_2,...,{\underset{\‾}{f}}_{\mathrm{nf}}],$

式中，Ξ为故障方向，f _nf为滤波后的故障幅值；并对故障数据矩阵

进行奇异值分解，得到：式中，D为对角矩阵，U为左奇异矩阵，V为右奇异矩阵，将对角阵D的非零奇异值按照降序排列，选故障方向Ξ＝U(d)，d为故障维数；

（2）根据步骤（1）中的故障方向，进行非线性故障重构及估计：

①将SPE指标作为重构检测指标，根据SPE指标进行故障幅值的大小的确定，则第i个样本的故障幅值f_i为：

f_i＝argminSPE(φ(x_i-Ξf_i))，i=1，2，…，n_f，

式中，argmin表示对目标函数最小化所对应的参数，φ表示样本在特征空间的映射向量；

②取故障幅值f最小值，故障重构就是沿故障方向减去纯故障因素的部分，构造出数据的正常因素部分，将重构后的正常数据假设为z，则

z=x-Ξf；   （1）

将SPE指标对f^T求偏导，在0点处获得极值，得到故障幅值f最小值，根据步骤①得到SPE指标为：

$\begin{array}{c}{\mathrm{SPE}}_i=\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)-\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T{P}_{\stackrel{\‾}{\mathrm{\φ}}}{P_{\stackrel{\‾}{\mathrm{\φ}}}}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)\\ =\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right)-\stackrel{\‾}{\mathrm{\φ}}{\left(z_i\right)}^T{\stackrel{\‾}{\mathrm{\Φ}}}^T{P}_{\stackrel{\‾}{k}}{M}^{-1}{P_{\stackrel{\‾}{k}}}^T\stackrel{\‾}{\mathrm{\Φ}}\stackrel{\‾}{\mathrm{\φ}}\left(z_i\right),\\ =\stackrel{\‾}{k}(z_i,z_i)-\stackrel{\‾}{k}{\left(z_i\right)}^T{M}_s\stackrel{\‾}{k}\left(z_i\right)\end{array}---\left(2\right)$

式中，i=1，2，…，n_f， $M^{-1}=\mathrm{diag}(\frac{1}{{\mathrm{\μ}}_1},\frac{1}{{\mathrm{\μ}}_2},...,\frac{1}{{\mathrm{\μ}}_p}),P_{\stackrel{\‾}{k}}=[{\mathrm{\β}}_1,{\mathrm{\β}}_2,...,{\mathrm{\β}}_p],M_s=P_{\stackrel{\‾}{k}}{M}^{-1}{P_{\stackrel{\‾}{k}}}^T,$ Μ_s是对称矩阵，μ_p为的特征值，β_p为与特征值μ_p对应的经过施密特标准正交化后的特征向量，

为均值中心化处理后的核矩阵；

由式（1）得到SPE指标对f^T求偏导，进而得到故障大小f^T为：

$f^T=\frac{[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})]B\left(z\right)\mathrm{\Ξ}}{[{1_n}^T+\stackrel{\‾}{k}{\left(z\right)}^T{M}_s(I-1_{n\×n})]k\left(z\right)},$

式中，k(z)为核函数， $B\left(z\right)=\left[\begin{array}{c}k(x_1,z){(x-x_1)}^T\\ k(x_2,z){(x-x_2)}^T\\ \·\\ \·\\ \·\\ k(x_n,z){(x-x_n)}^T\end{array}\right];$

将f^T转化为f的形式： $f=\frac{{\mathrm{\Ξ}}^T{B\left(z\right)}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}\left(z\right)]}{k{\left(z\right)}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}\left(z\right)]},---\left(3\right)$

联立式（1）和式（3），得到递推公式如下，

$f^{(j+1)}=\frac{{\mathrm{\Ξ}}^{T}B{(x-{\mathrm{\Ξf}}^{\left(j\right)})}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}(x-{\mathrm{\Ξf}}^{\left(j\right)})]}{k{(x-{\mathrm{\Ξf}}^{\left(j\right)})}^T[1_n+(I-1_{n\×n})M_s\stackrel{\‾}{k}(x-{\mathrm{\Ξf}}^{\left(j\right)})]};$

3）用多层递阶的方法对估计出的故障幅值f的发展趋势完成预测。

2.如权利要求1所述的一种基于非线性故障重构的故障预测方法，其特征在于：所述步骤1）中，异常检测故障提取过程中，采用衡量过程数据主要变化的HotellingT²统计量和衡量过程残差部分的SPE指标两类检测指标进行过程检测，其中SPE指标为平方预测误差指标，Hotelling T²统计量、SPE指标分别计算如下：

$T^2=k{\left(x\right)}^T{\mathrm{P\Λ}}^{-2}{P}^{T}k\left(x\right)\≤T_{\mathrm{\α}}^2,$

$\mathrm{SPE}=k(x,x)-k{\left(x\right)}^T{\mathrm{P\Λ}}^{-1}{P}^{T}k\left(x\right)\≤{\mathrm{\δ}}_{\mathrm{\α}}^2,$

式中，k(x)^T为核矩阵的转置矩阵；P为负载矩阵；Λ为协方差矩阵；

为HotellingT²统计量的控制限值；k(x,x)为样本x在高维特征空间的内积范数；

为SPE指标的控制限值；当两种检测指标分别超过了各自的控制限值时，则表明过程数据中出现了异常；若在控制限值内，则表明过程数据正常。

Images