CN103824137B - 一种复杂机械设备多工况故障预测方法 - Google Patents

Abstract

本发明涉及一种复杂机械设备多工况故障预测方法，其步骤：1）针对多工况过程建立多PCA模型，对每个PCA模型计算出相应的检测指标T²统计量和SPE；2）对多PCA模型中的各个T²统计量和SPE两检测指标进行优化，并对机械设备进行故障检测，检测得到过渡过程机械设备的故障数据；3）对优化后的两检测指标T²统计量和SPE检测到的过渡过程机械设备故障数据进行故障重构，得到使重构后SPE最小化的幅值估计值；4）对过渡过程中不同工况下的同一故障重构后的幅值估计值进行一致性幅值估计；5）根据进行一致性幅值估计后的幅值估计值

Description

一种复杂机械设备多工况故障预测方法

技术领域

本发明涉及一种机械设备故障预测方法，特别是关于一种复杂机械设备多工况故障预测方法。

背景技术

随着科学技术和工业发展，机械设备向着大型化、高速化、复杂化发展。因此，现在企业生产对设备及系统的可靠性、连续性、经济性等要求日益提高，在以往对设备及系统故障进行有效诊断和提出解决方案的基础上，进一步要求，在故障只出现微小异常征兆时即可实现对故障进行预报并提出相应紧急处理措施。故障预测的方法有多种多样，其中统计过程监控技术已经发展20余年，并且广泛应用于工业过程的故障检测、诊断和估计等。在最近的研究中，基于主元分析（Principal Component Analysis，PCA）的故障估计方法被成功用于故障预测中，但对于多工况过程的数据，基于主元分析的故障估计方法却不能够非常准确地进行故障预测，当复杂系统运行在多个工况时，变量之间的关系会根据系统当前的工作模式发生相应的变化，如果按照单工况的数据模型描述系统，会造成大量的误报和漏报。因此需要针对多工况过程的数据提出新的故障预测方案。

发明内容

针对上述问题，本发明的目的是提供一种复杂机械设备多工况故障预测方法，该方法的预测精度较高，能有效避免故障误报和漏报情况，避免了不同工况下同一种故障的估计误差。

为实现上述目的，本发明采取以下技术方案：一种复杂机械设备多工况故障预测方法，其包括以下步骤：1）针对多工况过程建立多PCA模型，并针对每个PCA模型计算出其相应的检测指标Hotelling’s T²统计量和SPE，其中Hotelling’s T²统计量简称T²统计量；（1）假设x∈R^m表示具有m个测量变量的样本向量，正常运行时的样本有n个；数据矩阵X∈R^n×m由n个样本组成，其中每一行代表了一个样本，每一列代表了一个测量变量共有n次采样；将数据矩阵X各列经过标准化处理成零均值和单位方差的变量，得到进行标准化后的样本x的协方差矩阵S，并对该协方差矩阵特征值分解并按大小降序排列；协方差矩阵S为：

$S=\mathrm{cov}\left(x\right)\≈\frac{1}{n-1}{X}^{T}X;$

（2）PCA模型将样本矩阵X∈R^n×m分解成建模部分和残差部分E两个部分：，其中，表示被建模部分；E表示残差部分；P∈R^m×A为负载矩阵，是由S的前A个特征向量组成的，A表示主元的个数；T∈R^n×A为得分矩阵，T=XP；（3）根据已有的测量变量在每个稳态工况下的历史数据建立相对应的主元模型，组成包含所有工况的多主元模型：，其中，q表示稳定工况个数；（4）在多PCA模型中，针对每个PCA模型需要计算出其相应的检测指标T²统计量和SPE，用SPE指标衡量样本向量在残差空间投影的变化，用T²统计量衡量测量变量在主元空间中的变化；2）采用随时间加权算法对多PCA模型中的各个T²统计量和SPE两检测指标进行优化，并根据优化后的检测指标T²统计量和SPE对机械设备进行故障检测，检测得到过渡过程机械设备的故障数据：（1）将多工况过渡过程中的上一工况的权重w₁分和下一工况的权重w₂分别取为：

$w_1=\frac{t-t_1}{\mathrm{\Δt}},$

$w_2=\frac{t_2-t}{\mathrm{\Δt}},$

式中，t₁表示过渡过程开始的时刻，t₂表示过渡过程结束的时刻，Δt=t₂-t₁；w₁+w₂=1；（2）假设t=k时刻，分别表示上一工况和下一工况的测量变量均值，σ₁和σ₂分别表示上一工况和下一工况的测量变量标准差，则此时过渡过程的均值和标准差σ_t＝k分别为：

${\stackrel{\‾}{x}}_{t=k}=w_1\×{\stackrel{\‾}{x}}_1+w_2\×{\stackrel{\‾}{x}}_2,$

σ_t＝k＝w₁×σ₁+w₂×σ₂；

（3）样本向量通过加权后的均值和标准差σ_t＝k进行预处理，求出样本向量的协方差矩阵S_t＝k和过渡过程的负载矩阵P_t＝k，进而得到优化后的SPE指标和T²统计量：

${\mathrm{SPE}}_{t=k}={\left|\right|(I-P_{t=k}{P_{t=k}}^T)x_{t=k}\left|\right|}^2\≤{\mathrm{\δ}}_{\mathrm{\αt}=k}^2,$

$T_{t=k}^2=x_{t=k}^T{P}_{t=k}{\mathrm{\Λ}}_{t=k}^{-1}{P}_{t=k}^T{x}_{t=k}\≤T_{\mathrm{\α}=k}^2;$

3）对步骤2）中优化后的两检测指标T²统计量和SPE检测到的过渡过程机械设备故障数据进行故障重构，得到使重构后SPE最小化的幅值估计值；4）对过渡过程中不同工况下的同一故障重构后的幅值估计值进行一致性幅值估计：（1）假设，式中，表示同一故障在上一工况和下一工况下的幅值估计值，Ξ_{i_1}、Ξ_{i_2}表示各自故障在残差空间的投影矩阵，假设维数为a；（2）根据步骤（1）得到同一故障在不同工况下的幅值估计值为：求解使最小化的d和e，得到同一故障在不同工况下的一致估计5）根据进行一致性幅值估计后的幅值估计值，利用支持向量机预测模型对故障幅值进行趋势预测。

所述步骤3）中，所述故障重构的过程如下：（1）假设检测得到的故障数据用x表示，正常数据用x^*表示，x^*，x∈[x₁,x₂,...,x_m]，m为每一个测量样本的维数，故障幅值的大小用f表示，则：x=x^*+Ξf，式中，Ξ表示故障子空间；（2）重构正常数据以消除故障的影响，得到正常数据的估计值为：x_x ^*=x-Ξ_if_i，式中，f_i表示对故障幅值大小f的估计；（3）根据步骤（2）重构后的正常数据的估计值得到使重构后SPE最小化的幅值估计值，式中，表示故障子空间Ξ到残差子空间的投影，C表示残差子空间。

所述步骤5）中，所述支持向量机预测方法如下：（1）对输入的测量变量采用降噪平滑处理进行数据的预处理，根据FPE准则计算嵌入维数，然后分别形成训练和测试样本；（2）对支持向量机模型进行初始化，将拉格朗日乘子α_i和以及阈值b赋以随机初始值；（3）利用训练样本建立目标函数，得到拉格朗日乘子α_i和以及阈值b，最后选取最优的惩罚因子C^*和核函数参数g；（4）将步骤（3）得到的参数值代入预测模型，用测试样本计算在未来某时刻的预测值；（5）计算误差函数，当误差的绝对值小于预先设定值时，结束学习过程，否则返回步骤（1）继续学习。

所述步骤5）中，所述支持向量机预测模型对故障幅值进行趋势预测精度采用均方误差和平均相对预测误差来评价。

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于采用基于多PCA故障重构的方法对多工况过程的故障进行预测，可以很好地解决过程数据的多工况问题，从隐含故障的数据中挖掘出故障方向并估计出故障幅值，也考虑到了故障的多维特性，可以获得更加精确的故障预测结果。2、由于复杂机械设备运行在多个工况时，测量变量之间的关系会根据设备当前的工作模式发生相应的变化。如果按照单工况的数据模型描述设备状态，会造成大量的误报和漏报。因此，本发明从过程监控的角度采用多主元分析模型，利用过渡过程的加权检测方法，提出了基于多主元分析模型的重构、故障估计和故障一致性估计，避免了不同工况下同一种故障的估计误差。3、本发明采用基于重构的多工况故障预测方法，当故障幅值被估计出来以后，采用支持向量机对故障幅值进行趋势预测。由于故障重构的目的是沿故障的方向用故障数据构造正常数据，消除故障对正常数据的影响，因此故障重构后进行预测的结果更为准确。4、本发明采用支持向量机对故障幅值进行直接预测，解决了传统预测方法无法利用多维有效数据的问题，在长期预测中，支持向量机模型具有很好的预测精度，具有很好的泛化性能。本发明可以广泛在机电设备故障预测中应用。

附图说明

图1是本发明在多PCA模型下的SPE指标检测仿真曲线示意图；

图2是本发明在多PCA模型下的T²统计量指标检测仿真曲线示意图；

图3是本发明的两个工况下的分别故障估计曲线示意图；

图4是本发明的两个工况下的分别故障一致估计曲线示意图；

图5是本发明基于重构后的故障预测曲线示意图。

具体实施方式

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

如图1所示，本发明提供一种机械设备多工况故障预测方法，其包括以下步骤：

1）针对多工况过程建立多PCA模型（主元分析模型），并针对每个PCA模型计算出其相应的检测指标Hotelling’s T²统计量（以下简称T2统计量）和SPE（平方预测误差，也称Q统计量）；

（1）假设x∈R^m表示具有m个测量变量的样本向量（即m为样本x的维数），正常运行时的样本有n个。数据矩阵X∈R^n×m由n个样本组成，其中每一行代表了一个样本，每一列代表了一个测量变量共有n次采样。将数据矩阵X各列经过标准化处理成零均值和单位方差的变量，可以得到进行标准化后的样本x的协方差矩阵S，并对该协方差矩阵特征值分解并按大小降序排列。协方差矩阵S为：

$S=\mathrm{cov}\left(x\right)\≈\frac{1}{n-1}{X}^{T}X,---\left(1\right)$

其中，将数据矩阵X各列经过标准化处理成零均值和单位方差的变量的方法是将数据矩阵X的每一列减去相应的变量均值并且除以相应的变量标准差。

（2）根据PCA模型将测量变量空间分成主元子空间和残差子空间两个正交且互补的子空间，任意一个样本向量均可分解成为在主元子空间和残差子空间上的投影，即PCA模型将样本矩阵X∈R^n×m分解成建模部分和残差部分E两个部分：

$X=\hat{X}+E{=\mathrm{TP}}^T+E,---\left(2\right)$

其中，表示被建模部分；E表示残差部分；P∈R^m×A为负载矩阵，是由S的前A个特征向量组成的，A表示主元的个数；T∈R^n×A为得分矩阵，T=XP。

（3）根据已有的测量变量在每个稳态工况下的历史数据建立相对应的主元模型，组成包含所有工况的多主元模型：

$X^{\left(q\right)}={\hat{X}}^{\left(q\right)}+E^{\left(q\right)}---\left(3\right)$

其中，q表示稳定工况个数。

（4）在多PCA模型中，针对每个PCA模型需要计算出其相应的检测指标T²统计量和SPE，即T^2(q)和SPE^(q)，用SPE指标衡量样本向量在残差空间投影的变化，用T²统计量衡量测量变量在主元空间中的变化：

其中，SPE指标表达式为：

$\mathrm{SPE}={\left|\right|(I-{\mathrm{PP}}^T)x\left|\right|}^2\≤{\mathrm{\δ}}_{\mathrm{\α}}^2---\left(4\right)$

式中，I为单位矩阵；表示置信水平为α时SPE的控制限。当SPE在控制限内时，认为当前运转过程处于正常状态。当SPE值超出了统计控制限时，代表当前运转过程发生了故障，SPE值的变化代表着数据间相关性的变化。该控制限的计算公式为：

${\mathrm{\δ}}_{\mathrm{\α}}^2={\mathrm{\θ}}_1{(\frac{c_{\mathrm{\α}}\sqrt{{2\mathrm{\θ}}_2{h}_0^2}}{{\mathrm{\θ}}_1}+1+\frac{{\mathrm{\θ}}_2{h}_0(h_0-1)}{{\mathrm{\θ}}_1^2})}^{1/h_0}---\left(5\right)$

式中，，λ_j为样本矩阵X的协方差矩阵∑的特征值，c_α为标准正态分布在置信水平α下的阈值，m是样本x的维数。

T²统计量表达式为：

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

其中，Λ＝diag{λ₁,...,λ_A}，表示置信度为α的T²统计限。当T²位于控制限内时，认为当前运转过程处于正常工作状态。

由此可知，主元分析是根据数据变化的方差大小来确定变化方向的主次地位，按主次顺序得到各个主元，这些主元彼此是相互独立的。本发明采用SPE和T²统计量两种指标来检测机械设备运转过程是否发生异常，T²统计量指标衡量的是样本距离主元子空间原点的距离，而SPE指标主要体现了正常过程变量之间的相关性被改变的程度，显示了异常的过程状态。主元子空间主要反映正常数据变化的测试，残差子空间主要反映噪声情况，因此T²统计量定义的控制限比SPE定义的要大得多，较小的故障可能不会超过T²统计量控制限但却可以超过SPE控制限。假设在正常的状态下，各个过程变量之间的关系是不变的，系统的动态变化主要表现在主元子空间中。当有故障发生时，残差子空间将会发生变化。因此采用SPE指标来检测残差空间中是否发生异常，T²统计量衡量测量变量在主元空间中的变化。

2）采用随时间加权算法对多PCA模型中的各个T²统计量和SPE两检测指标进行优化，并根据优化后的检测指标T²统计量和SPE对机械设备进行故障检测，检测得到过渡过程机械设备的故障数据，通过优化后的检测指标进行故障检测可以有效避免在工况过渡过程中出现故障误报；

其中，对T²统计量和SPE两检测指标进行优化的方法包括以下步骤：

（1）由于在多工况的过渡过程中每一个样本向量点的权重值都是不同的，都会随时间逐渐增加，上一工况的权重值从1逐渐到0，下一工况的权重值从0逐渐到1。因此，将多工况过渡过程中的上一工况的权重w₁分和下一工况的权重w₂分别取为：

$w_1=\frac{t-t_1}{\mathrm{\Δt}}---\left(7\right)$

$w_2=\frac{t_2-t}{\mathrm{\Δt}}---\left(8\right)$

式中，t₁表示过渡过程开始的时刻，t₂表示过渡过程结束的时刻，Δt=t₂-t₁；w₁+w₂=1。

当一个稳定运行的工况刚刚向下一个工况进行过渡时，过渡过程的特性受该工况的影响较大，随着时间的推移，越靠近新的工况的时间点上，过渡过程的特性越贴近下一个工况，直到进入下一工况，因此过渡过程中的每一个点的权重值都是不同的。

（2）假设t=k时刻，分别表示上一工况和下一工况的测量变量均值，σ₁和σ₂分别表示上一工况和下一工况的测量变量标准差，则此时过渡过程的均值和标准差σ_t＝k分别为：

${\stackrel{\‾}{x}}_{t=k}=w_1\×{\stackrel{\‾}{x}}_1+w_2\×{\stackrel{\‾}{x}}_2---\left(9\right)$

σ_t=k＝w₁×σ₁+w₂×σ₂。 (10)

（3）样本向量通过步骤（2）中求出的加权后的均值和标准差σ_t=k进行预处理，求出样本向量的协方差矩阵S_t=k和过渡过程的负载矩阵P_t=k，进而得到优化后的SPE指标和T²统计量：

${\mathrm{SPE}}_{t=k}={\left|\right|(I-P_{t=k}{P_{t=k}}^T)x_{t=k}\left|\right|}^2\≤{\mathrm{\δ}}_{\mathrm{\αt}=k}^2,---\left(11\right)$

$T_{t=k}^2=x_{t=k}^T{P}_{t=k}{\mathrm{\Λ}}_{t=k}^{-1}{P}_{t=k}^T{x}_{t=k}\≤T_{\mathrm{\α}=k}^2.---\left(12\right)$

如图1所示，SPE指标的检测仿真曲线，横坐标表示样本数，纵坐标表示SPE值。虚线代表各个统计量99%的控制限，若超过控制限则表示系统出现故障。由图可知，基于多PCA方法的SPE在第166个样本点超过控制限。

如图2所示，T²统计量指标的检测仿真曲线，横坐标表示样本数，纵坐标表示T²值。虚线代表各个统计量99%的控制限，若超过控制限则表示系统出现故障。由图可知，基于多PCA方法的T²在第164个样本点超过控制限。

3）对步骤2）中优化后的两检测指标T²统计量和SPE检测到的过渡过程机械设备故障数据进行故障重构，得到使重构后SPE最小化的幅值估计值，以使故障因素对故障数据内正常部分的影响降至最小，进而构造出正常状态下的过程数据。

故障重构的具体过程如下：

（1）假设检测得到的故障数据用x表示（即x为样本），正常数据用x^*表示，x^*，x∈[x₁,x₂,...,x_m]，m代表传感的个数，即每一个测量样本的维数，故障幅值的大小用f表示，则：

x=x^*+Ξf， (13)

式中，Ξ表示故障子空间（也称为故障方向矩阵）。

（2）重构正常数据以消除故障的影响，得到正常数据的估计值为：

x_i ^*=x-Ξ_if_i (14)

式中，f_i表示对故障幅值大小f的估计，从几何意义上讲，是将样本x沿着故障所在的子空间拉回到主元子空间。

（3）根据步骤（2）重构后的正常数据的估计值，得到使重构后SPE最小化的幅值估计值

${\hat{f}}_i={\widetilde{\mathrm{\Ξ}}}^{+}{x}_i,---\left(15\right)$

式中，,表示故障子空间Ξ到残差子空间的投影，C表示残差子空间。

其中，重构后的SPE指标为：

$\mathrm{SPE}\left(x^{*}\right)={\left|\right|{\tilde{x}}^{*}\left|\right|}^2={\left|\right|\tilde{x}-{\widetilde{\mathrm{\Ξf}}}_i\left|\right|}^2---\left(16\right)$

重构就是根据公式（16）搜索，使得：

${\hat{f}}_i=\arg \min \left|\right|\tilde{x}-{\widetilde{\mathrm{\Ξf}}}_i{\left|\right|}^2={\left({\widetilde{\mathrm{\Ξ}}}_i^T{\widetilde{\mathrm{\Ξ}}}_i\right)}^{-1}{\widetilde{\mathrm{\Ξ}}}_i^T\tilde{x}---\left(17\right)$

式（17）的最优解为：

${\hat{f}}_i={\widetilde{\mathrm{\Ξ}}}^{+}{\tilde{x}}_i={\widetilde{\mathrm{\Ξ}}}^{+}{x}_i.$

其中，表示的Moore-Penrose伪逆，表示使重构后SPE最小化的故障幅值的估计。

4）对过渡过程中不同工况下的同一故障重构后的幅值估计值进行一致性幅值估计：当故障出现在工况一的过程中，而且延续到了工况二的过程中，这时的故障方向矩阵Ξ发生了变化，但是故障的幅值没有发生变化。由于是同一个故障，即使在不同的工况下被重构估计，理论上也应该会得到一致估计，可是在实际工业过程中，数据的噪声和不同工况下机器的干扰等因素都会使结果产生偏差，因此需要对不同工况下的同一故障进行一致性幅值估计。

（1）假设

$A={\hat{f}}_1{\mathrm{\Ξ}}_{i\_1}^{\′}---\left(18\right)$

$B={\hat{f}}_2{\mathrm{\Ξ}}_{i\_2}^{\′}---\left(19\right)$

式中，表示同一故障在上一工况和下一工况下的幅值估计值（即不同工况下的幅值估计值），Ξ_{i_1}、Ξ_{i_2}表示各自故障在残差空间的投影矩阵，假设维数为a。

（2）根据公式（18）、（19）可以得到同一故障在不同工况下的幅值估计值为：

$f_1^{*}=\mathrm{Ad}---\left(20\right)$

$f_2^{*}=\mathrm{Be}---\left(21\right)$

求解使最小化的d和e，其求解方法如下：

令矩阵G=[A,B]，G∈R^n×2k对G^T进行奇异值分解（SVD）：

U^HG^TV=Σ(22)

其中，U和V都是酉矩阵维数分别为2k和n，Σ是对角阵，对角线上元素为奇异值。按照由大到小对奇异值进行排序，取左奇异向量U的后a列特征向量，此时阶数是2k×k，前k行即为d的估计，后k行即为e的估计，得到d和e后代入公式（20）和公式（21）中，得到的即为同一故障在不同工况下的一致估计。

如图3所示，假设工况一先发生故障，然后转换工况，故障进入工况二。如图所示，样本点从0到900为工况一，故障是在样本点160处发生，样本点900到1100是过渡过程，从样本点1100开始进入工况二。故障估计的曲线理论上应该是在工况一和工况二时误差较小，但如图3所示，工况一的故障估计和工况二的故障估计误差较大。

如图4所示，假设工况一先发生故障，然后转换工况，故障进入工况二。如图所示，样本点从0到900为工况一，故障是在样本点160处发生，样本点900到1100是过渡过程，从样本点1100开始进入工况二。经过一致幅值估计的算法后，得到的故障估计如图4所示，这时的故障估计达到一致，经过过渡过程后曲线更平滑。

5）根据进行一致性幅值估计后的幅值估计值，利用支持向量机(SVM)预测模型对故障幅值进行趋势预测。其预测方法如下：

（1）对输入的测量变量采用降噪平滑处理进行数据的预处理，根据FPE准则计算嵌入维数，然后分别形成训练和测试样本。

（2）对支持向量机模型进行初始化，将拉格朗日乘子α_i和以及阈值b赋以随机初始值。

（3）利用训练样本建立目标函数，得到拉格朗日乘子α_i和以及阈值b，最后选取最优的惩罚因子C^*和核函数参数g。

（4）将步骤（3）得到的参数值代入预测模型，用测试样本计算在未来某时刻的预测值。

（5）计算误差函数，当误差的绝对值小于预先设定值时，结束学习过程，否则返回步骤（1）继续学习。

例如，对于给定的时间序列{x₁,x₂,...,x_n},i＝1,2,...,n，x_n是预测的目标值，将先前的{x_n-1,x_n-2,...,x_n-m}作为输入值，建立自相关输入与输出y_n＝{x_n}之间的映射关系f：，m为嵌入维数，经过变换之后，得到用于支持向量机的学习样本：

$X=\left[\begin{array}{cccc}{x}_1& x_2& ...& x_m\\ x_2& x_3& ...& x_{m+1}\\ ...& ...& ...& ...\\ x_{n-m}& x_{n-m+1}& ...& x_{n-1}\end{array}\right],Y=\left[\begin{array}{c}{x}_{m+1}\\ x_{m+2}\\ ...\\ x_n\end{array}\right]---\left(23\right)$

对向量机进行训练的回归函数为：

$y_t={\mathrm{\Σ}}_{i=1}^{n-m}({\mathrm{\α}}_i-{\mathrm{\α}}_i^{*})K(x_i,x_t)+b,t=m+1,...,n---\left(24\right)$

其中，α_i和是拉格朗日乘子，b是阈值，K(x_i,x_t)是核函数，本发明采用径向基核函数（Radial Basis Function）：

$K(x,x\′)=\exp (-\frac{{\left|\right|x-x\′\left|\right|}^2}{{2\mathrm{\σ}}^2})---\left(25\right)$

则一步预测模型为：

$y_{n+1}={\mathrm{\Σ}}_{i=1}^{n-m}(\mathrm{\α}-{\mathrm{\α}}_i^{*})K(x_i^{*},x_{n-m+1}^{*})+b---\left(26\right)$

由此得到一个样本，表示第n+1个数据的预测值。

进而有

$y_{n+p}={\mathrm{\Σ}}_{i=1}^{n-m}({\mathrm{\α}}_i-{\mathrm{\α}}_i^{*})K(x_i^{*},x_{n-m+p}^{*})+b---\left(27\right)$

式（27）表示第p步的预测模型，得x_n表示第n个数据的实际值，表示第n个数据的预测值。

上述步骤5）中，SVM预测模型对故障幅值进行趋势预测精度可以采用均方误差(MSE)和平均相对预测误差(ARE)来评价，MSE和ARE越小预测精度越精确，本发明优选将MSE和ARE控制在5%以内。其中，均方误差和平均相对预测误差分别为：

均方误差为：

$\mathrm{MSE}=\frac{1}{n}\×\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-{\hat{y}}_i)}^2}---\left(28\right)$

平均相对预测误差为：

$\mathrm{ARE}=\frac{1}{n}\×\mathrm{\Σ}\left|\frac{y_i-{\hat{y}}_i}{y_i}\right|.---\left(29\right)$

例如，如图5所示，选取2000个样本数据用于故障检测，故障在第166个样本点被SPE指标检测出来。当故障被检测出来以后，利用故障重构对故障幅值进行估计，得到相应的故障估计。采用重构后得到的故障估计进行故障预测，利用支持向量机预测模型对未来一步和十步的故障估计进行预测。图5中实线表示实际曲线，点划线表示一步预测曲线，星型线表示十步预测曲线，则一步预测和十步预测的相对预测误差如表1所示。

表1不同预测步数下的相对预测误差

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

Claims (4)

1.一种复杂机械设备多工况故障预测方法，其包括以下步骤：

1)针对多工况过程建立多PCA模型，并针对每个PCA模型计算出其相应的检测指标Hotelling’s T²统计量和SPE，其中Hotelling’s T²统计量简称T²统计量；

(1)假设x∈R^m表示具有m个测量变量的样本向量，正常运行时的样本有n个；数据矩阵X∈R^n×m由n个样本组成，其中每一行代表了一个样本，每一列代表了一个测量变量共有n次采样；将数据矩阵X各列经过标准化处理成零均值和单位方差的变量，得到进行标准化后的样本x的协方差矩阵S，并对该协方差矩阵特征值分解并按大小降序排列；协方差矩阵S为：

$S=\mathrm{cov}\left(x\right)\≈\frac{1}{n-1}{X}^{T}X;$

(2)PCA模型将样本矩阵X∈R^n×m分解成建模部分和残差部分E两个部分：

$X=\hat{X}+E={\mathrm{TP}}^T+E,$

其中，表示被建模部分；E表示残差部分；P∈R^m×A为负载矩阵，是由S的前A个特征向量组成的，A表示主元的个数；T∈R^n×A为得分矩阵，T＝XP；

(3)根据已有的测量变量在每个稳态工况下的历史数据建立相对应的主元模型，组成包含所有工况的多主元模型：

$X^{\left(q\right)}={\hat{X}}^{\left(q\right)}+E^{\left(q\right)},$

其中，q表示稳定工况个数；

(4)在多PCA模型中，针对每个PCA模型需要计算出其相应的检测指标T²统计量和SPE，用SPE指标衡量样本向量在残差空间投影的变化，用T²统计量衡量测量变量在主元空间中的变化；

2)采用随时间加权算法对多PCA模型中的各个T²统计量和SPE两检测指标进行优化，并根据优化后的检测指标T²统计量和SPE对机械设备进行故障检测，检测得到过渡过程机械设备的故障数据：

(1)将多工况过渡过程中的上一工况的权重w₁分和下一工况的权重w₂分别取为：

$w_1=\frac{t-t_1}{\mathrm{\Δ}t},$

$w_2=\frac{t_2-t}{\mathrm{\Δ}t},$

式中，t₁表示过渡过程开始的时刻，t₂表示过渡过程结束的时刻，Δt＝t₂-t₁；w₁+w₂＝1；

(2)假设t＝k时刻，分别表示上一工况和下一工况的测量变量均值，σ₁和σ₂分别表示上一工况和下一工况的测量变量标准差，则此时过渡过程的均值和标准差σ_t＝k分别为：

${\stackrel{\‾}{x}}_{t=k}=w_1\×{\stackrel{\‾}{x}}_1+w_2\×{\stackrel{\‾}{x}}_2,$

σ_t＝k＝w₁×σ₁+w₂×σ₂；

(3)样本向量通过加权后的均值和标准差σ_t＝k进行预处理，求出样本向量的协方差矩阵S_t＝k和过渡过程的负载矩阵P_t＝k，进而得到优化后的SPE指标和T²统计量：

${\mathrm{SPE}}_{t=k}=\left|\right|(I-P_{t=k}{P_{t=k}}^T)x_{t=k}|{|}^2\≤{\mathrm{\δ}}_{\mathrm{\α}t=k}^2,$

$T_{t=k}^2=x_{t=k}^T{P}_{t=k}{\mathrm{\Λ}}_{t=k}^{-1}{P}_{t=k}^T{x}_{t=k}\≤T_{\mathrm{\α}=k}^2;$

3)对步骤2)中优化后的两检测指标T²统计量和SPE检测到的过渡过程机械设备故障数据进行故障重构，得到使重构后SPE最小化的幅值估计值

4)对过渡过程中不同工况下的同一故障重构后的幅值估计值进行一致性幅值估计：

(1)假设

$A={\hat{f}}_1{\mathrm{\Ξ}}_{i\_1}^{\′},$

$B={\hat{f}}_2{\mathrm{\Ξ}}_{i\_2}^{\′},$

式中，表示同一故障在上一工况和下一工况下的幅值估计值，Ξ_{i_1}、Ξ_{i_2}表示各自故障在残差空间的投影矩阵，假设维数为a；

(2)根据步骤(1)得到同一故障在不同工况下的幅值估计值为：

$f_1^{*}=Ad,$

$f_2^{*}=Be,$

求解使最小化的d和e，得到同一故障在不同工况下的一致估计f₁ ^*、

5)根据进行一致性幅值估计后的幅值估计值利用支持向量机预测模型对故障幅值进行趋势预测。

2.如权利要求1所述的一种复杂机械设备多工况故障预测方法，其特征在于：所述步骤3)中，所述故障重构的过程如下：

(1)假设检测得到的故障数据用x表示，正常数据用x^*表示，x^*，x∈[x₁,x₂,...,x_m]，m为每一个测量样本的维数，故障幅值的大小用f表示，则：

x＝x^*+Ξf，

式中，Ξ表示故障子空间；

(2)重构正常数据以消除故障的影响，得到正常数据的估计值为：

x_i ^*＝x-Ξ_if_i，

式中，f_i表示对故障幅值大小f的估计；

(3)根据步骤(2)重构后的正常数据的估计值得到使重构后SPE最小化的幅值估计值

${\hat{f}}_i={\widetilde{\mathrm{\Ξ}}}^{+}{x}_i,$

式中，表示故障子空间Ξ到残差子空间的投影，C表示残差子空间。

3.如权利要求1或2所述的一种复杂机械设备多工况故障预测方法，其特征在于：所述步骤5)中，所述支持向量机预测方法如下：

(1)对输入的测量变量采用降噪平滑处理进行数据的预处理，根据FPE准则计算嵌入维数，然后分别形成训练和测试样本；

(2)对支持向量机模型进行初始化，将拉格朗日乘子α_i和以及阈值b赋以随机初始值；

(3)利用训练样本建立目标函数，得到拉格朗日乘子α_i和以及阈值b，最后选取最优的惩罚因子C^*和核函数参数g；

(4)将步骤5)中的步骤(3)得到的参数值代入预测模型，用测试样本计算在未来某时刻的预测值；

(5)计算误差函数，当误差的绝对值小于预先设定值时，结束学习过程，否则返回步骤(1)继续学习。

4.如权利要求1或2所述的一种复杂机械设备多工况故障预测方法，其特征在于：所述步骤5)中，所述支持向量机预测模型对故障幅值进行趋势预测精度采用均方误差和平均相对预测误差来评价。