CN104714537A - 一种基于联合相对变化分析和自回归模型的故障预测方法 - Google Patents

Abstract

本发明公开了一种基于联合相对分析和自回归建模的故障预测方法，首先在主成分分析的基础上，提出了一种基于联合相对分析方法分解故障影响，确定故障方向。随后，根据确定的故障方向，基于重构技术估计故障幅度，进行正常部分的数据恢复。定义新监测统计量D²涵盖正常的数据波动，进而建立事件报警限Elt_D。最后，基于新监测统计量D²建立自回归模型用于在线监测统计量的预测，实现故障情况的提前报警。本发明简便有效，不依赖先验过程知识和假设。故障预测结果对后续的故障诊断和修复具有重要意义，有助于过程工程师对过程运行状态做出及时判断，从而保证工业生产的安全可靠运行和产品的高质量追求。

Description

一种基于联合相对变化分析和自回归模型的故障预测方法

技术领域

本发明属于复杂工业过程故障建模与诊断领域，特别是涉及一种在线的过程故障预测方法。

背景技术

现代工业过程生产设备众多、工艺原理日趋复杂，实时故障检测及诊断技术对于保证操作安全和提升质量具有重要作用。随着市场竞争的加剧和对生产过程安全可靠运行的强烈需求，过程在线监测、故障建模和诊断技术受到了越来越多的关注。在过去的几十年中，诸如主成分分析、偏最小二乘为代表的多变量统计分析技术由于具有不依赖过程知识、操作简便等优良特性，已在大量工业实际中得到了应用，如钢铁冶金、精细化工、聚合物反应等等。这些方法从测量数据中获取潜在的过程特性，并定义特定的监测模型和统计量来判断当前过程是否处于“统计控制状态”来指示过程是否正常进行。“统计控制状态”即涵盖所有正常波动的操作区域，当过程数据不在此区域时，认为过程正在发生异常变化。

通常，故障检测及诊断包括异常过程行为的检测和故障原因的辨识。然而，监测统计量只能指示是否发生过程异常，而对于是否会有严重的故障发生或是故障将在什么时候发生不能做到有效的预测。对于一些故障，一旦监测统计量超出其置信限会导致严重的问题，甚至会带来灾难性的后果。因此，期望在监测统计量报警后依靠工程师来进行过程调整和补救将为时已晚。对于一些故障呈缓慢变化的复杂工业过程，如果能提前预测故障影响、准确获悉事故发生时间以便采取补救措施具有重大意义。

一般来说，故障预测有两个重要问题需要解决：1)如何量化故障影响；2)故障影响将于何时超越置信限。目前，部分方法假设过程故障可以直接观测得到，与大部分实际情况不符；一些方法针对此问题，利用统计方法提取故障的特征子空间或方向，但是未能很好的评价过程影响。可以说，目前已有的方法还没有形成对故障过程的有效评测和预测技术，从而直接或间接的影响后续的故障诊断的性能和采取补救措施的可能性。

发明内容

本发明的目的在于针对现有故障建模及预测方法的不足，提供一种基于联合相对变化分析和自回归模型的故障预测方法。该方法能及时捕捉故障数据下的过程变化，有效确定故障变化方向。基于故障方向，构建预测模型进行监测统计量的预测，以便进行及时可靠的在线过程诊断，并最终应用于实际工业生产现场，确保复杂工业过程的安全可靠运行以及产品的高质量追求。

本发明的目的是通过以下技术方案来实现的：一种基于联合相对变化分析和自回归模型的故障预测方法，该方法包括以下步骤：

步骤1：过程分析与数据准备

设过程有J个测量变量和K个采样点，则可以形成一个K×J矩阵；用X_r(N_r×J)和X_a(N_a×J)分别描述正常数据集和故障数据集；N_r和N_a分别代表正常数据集和故障数据集的样本个数；对数据集X_r和X_a按式(1)进行中心化和归一化处理，使其均值为0，标准差为1；

$x_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_j}{{\mathrm{\δ}}_j}---\left(23\right)$

其中x_k,j是数据集中第k行第j列的数据，k∈[1,K]是采样时刻，j∈[1,J]是变量索引；是矩阵第j个变量的均值，δ_j是第j个变量的标准差，满足以下关系：

${\stackrel{\‾}{x}}_j=\frac{1}{K}\underset{k}{\overset{K}{\mathrm{\Σ}}}{x}_{k,j}$

${\mathrm{\δ}}_j=\sqrt{\frac{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(x_{k,j}-{\stackrel{\‾}{x}}_j)}^2}{K-1}}---\left(24\right)$

步骤2：联合相对变化分析

(2.1)对正常数据集X_r进行主成分分析(PCA)获取监测模型：

$\left\{\begin{array}{c}{T}_r=X_r{P}_r\\ E_r=X_r{P}_r^e{P}_r^{\mathrm{eT}}\end{array}\right.---\left(25\right)$

其中T_r(N_r×R_r)和P_r(J×R_r)分别是主元成分和相应的主元负载；E_r(N_r×J)和是PCA分解得到的残差矩阵和相应的负载；R_r是通过累积解释方差保留的主元空间的主元个数；是残差空间中保留的主元方向个数，有

根据公式(3)计算获取的主元得分T_r和残差矩阵E_r，计算两个监测统计指标：T²统计量和SPE统计量，计算公式如下：

${T_r}^2={(t_r-\stackrel{\‾}{t_r})}^T{{\mathrm{\Σ}}_r}^{-1}(t_r-\stackrel{\‾}{t_r})---\left(26\right)$

SPE_r＝e_r ^Te_r

其中，t_r(R×1)是来自T_r(N_r×R_r)的主元向量；是T_r(N_r×R_r)的均值向量，由于归一化的原因其值为0；Σ_r是T_r(N_r×R_r)的协方差矩阵；e_r(J×1)是E_r(N_r×J)的残差向量；

T²统计限和SPE统计限分别可以通过因子为α的F分布和χ²分布获取：

$T^2~\frac{R(N^2-1)}{N(N-R)}{F}_{R,N-R,\mathrm{\α}}---\left(27\right)$

SPE～gχ_h,α ²   (28)

其中g＝ν/2m、h＝2(m)²/ν；m是公式(4)中计算得到的SPE平均值，v是SPE的方差；

(2.2)利用正常监测模型检测故障数据

将故障数据X_a向P_r投影，获取故障数据的主元T_a；将X_a向投影，获得故障数据的残差E_a，如公式(7)所示：

$\left\{\begin{array}{c}{T}_a=X_a{P}_r\\ E_a=X_a{P}_r^e{P}_r^{\mathrm{eT}}=\underset{i=1}{\overset{R_e}{\mathrm{\Σ}}}{p}_{r,i}^e{p_{r,i}^e}^T\end{array}\right.---\left(29\right)$

沿J个不同的监测方向，分别对比故障数据集和正常数据集在主元空间和残差空间的波动变化：

$\begin{array}{c}{\mathrm{Ratio}}_i=\frac{\mathrm{var}\left(T_{a,i}\right)}{\mathrm{var}\left(T_{r,i}\right)}(i=\mathrm{1,2},...,R_r)\\ {\mathrm{\Δ}}_i={\left|\right|X_a{p}_{r,i}^e{p_{r,t}^e}^T\left|\right|}^2-{\left|\right|X_r{p}_{r,i}^e{p_{r,i}^e}^T\left|\right|}^2(i=\mathrm{1,2},...,R_r^e)\end{array}---\left(30\right)$

其中，var(*)表示距离中心点(从正常数据集中获取，此处为0)的主元方差；下标i代表矩阵的第i列；Ratio是由Ratio_i组成的R_r维向量；||||是欧氏距离，Δ是由Δ_i构成的维向量；

(2.3)故障相关方向排序

对Ratio和Δ进行排序，找出所有Ratio_i>1和Δ_i>0的方向，这些方向在主元空间和残差空间中具有更大的波动变化，是可能的故障方向；

(2.4)更新监测统计量

移除步骤2.3找到的故障方向，并重新计算正常数据集和故障数据集新的T²和SPE监测统计量；定义正常数据集新的监测统计限，并与故障数据的监测统计量进行对比；如果故障数据的监测统计量在正常数据的监测统计限内，认为关键故障方向已经得到移除；否则表明故障数据中仍存在故障方向；

(2.5)迭代更新

重复步骤(2.2)-(2.4)直至所有故障方向都被移除；

将提取的故障方向构成数据矩阵和作为两个不同监测子空间的故障方向，保留的主元方向个数分别是和

步骤3：故障方向集成

沿和X_a的波动可以建模为：

$\begin{array}{c}{X}_{a,f}=X_a{P}_r^{*}{P}_r^{*T}\\ X_{a,f}^e=X_a{P}_r^{e*}{P}_r^{e*T}\end{array}---\left(31\right)$

对X_a,f和进行PCA分解，获得系统主要变化P_f(J×R_f)和其中，R_f和分别为主元个数，且有 $R_f=\mathrm{rank}\left(X_a{P}_r^{*}\right),P_f^{*}=\mathrm{rnak}\left(X_a{P}_r^{e*}\right);$ 集成P_f和为统一的联合模型 $P_{\mathrm{fc}}(J\×R_{\mathrm{fc}})=[P_f,P_f^{*}](R_{\mathrm{fc}}=R_f^{*}+R_f);{}^{}$

步骤4：基于重构的故障幅度估计

(4.1)正常数据的重构

由步骤3分解得到故障方向P_fc，不受故障影响的数据能够恢复为：

x^*＝x_f-P_fcf   (32)

其中，x_f是故障样本，x^*是从x_f中恢复得到的正常数据，f代表系统子空间的故障幅度；

(4.2)故障幅度的计算

x^*的最优估计通过最小化x^*到监测子空间(主元空间或是残差空间)的距离来得到；故障幅度f计算如下：

$\begin{array}{c}f=\arg \min {\left|\right|\mathrm{\Ω}{x}^{*}\left|\right|}^2=\arg \min {\left|\right|\mathrm{\Ω}(x-P_{\mathrm{fc}}f)\left|\right|}^2=\arg \min {\left|\right|\tilde{x}-{\tilde{P}}_{\mathrm{fc}}f\left|\right|}^2\\ ={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}\end{array}---\left(33\right)$

其中，Ω表示监测空间的投影算子，在主元空间Ω＝P_r(P_r ^TP_r)^-1P_r ^T，在残差空间Ω＝I-P_r(P_r ^TP_r)^-1P_r ^T；I是J×J的单位矩阵；

当Ω＝P_r(P_r ^TP_r)^-1P_r时：

$\begin{array}{c}{\tilde{P}}_{\mathrm{fc}}=\mathrm{\Ω}{P}_{\mathrm{fc}}=P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r{P}_{\mathrm{fc}}=P_r{\left({P_r}^T{P}_r\right)}^{-1}{P_r}^T[P_f,P_f^{*}]=P_f\\ f={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({P_f}^T{P}_f\right)}^{-1}{P_f}^T{P}_r{\left({P_r}^T{P}_r\right)}^{-1}{P_r}^T={P_f}^{T}x\end{array}---\left(34\right)$

当Ω＝I-P_r(P_r ^TP_r)^-1P_r时

$\begin{array}{c}{\tilde{P}}_{\mathrm{fc}}=\mathrm{\Ω}{P}_{\mathrm{fc}}=(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)P_{\mathrm{fc}}=(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)[P_f,P_f^{*}]=P_f^{*}\\ f^{*}={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left(P_f^{*T}{P}_f^{*}\right)}^{-1}(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{p}_r)x=P_f^{*T}x\end{array}---\left(35\right)$

主元空间和残差空间的综合故障幅度f_c计算为：

$f_c=\left[\begin{array}{c}f\\ f^{*}\end{array}\right]=\left[\begin{array}{c}{P_f}^T\\ P_f^{*T}\end{array}\right]x={P_{\mathrm{fc}}}^{T}x---\left(36\right)$

步骤5：联合故障监测指标的确定

(5.1)联合故障监测统计指标D²实际值的计算

对于多变量系统故障，通过典型的马氏距离统计量计算得到故障影响：

$D^2={(f_c-\stackrel{\‾}{f_c})}^T{\mathrm{\Λ}}^{-1}(f_c-\stackrel{\‾}{f_c})---\left(37\right)$

其中，表示由正常数据计算得到的故障幅度平均值；Λ为对角矩阵，其对角元素是由正常数据计算得到的每个方向上的主元方差，对于一些非常接近于0的方差值赋予1；

(5.2)联合故障监测统计限Ctr_D的建立

在正常操作条件下，假设过程数据服从多变量正态分布，当样本数N足够大时，可以利用χ²分布计算D²的置信限：

$D^2\≤{\mathrm{Ctr}}_D\≡g{\mathrm{\χ}}_{h,a}^2---\left(38\right)$

其中，g＝ν/2m、h＝2(m)²/ν,m是在正常训练数据下根据公式(15)计算得到的D²平均值，ν是相应的方差；

步骤6：故障预测

(6.1)预测模型输入输出数据整理

记L代表输入矩阵的时间长度，在每一个采样时刻，具有一个R_fc维的故障幅度向量，输入矩阵X(N×R_fcL)表示为：

$X(N\×R_{\mathrm{fc}}L)=\left[\begin{array}{c}{f_{c,1}}^T,{f_{c,2}}^T,...,{f_{c,L}}^T\\ {f_{c,2}}^T,{f_{c,3}}^T,...,{f_{c,L+1}}^T\\ .\\ .\\ .\\ {f_{c,N}}^T,{f_{c,N+1}}^T,...,{f_{c,L+N-1}}^T\end{array}\right]---\left(39\right)$

其中，N为输入矩阵的样本数，每一个x_i＝[f_c,i ^T,f_c,i+1 ^T,...,f_c,i+L-1 ^T]是一个R_fcL维的向量，包含了从时刻i到i+L-1的故障幅度信息；

预测时域PH表示从当前时刻对PH步后的故障幅度进行预测；输出矩阵Y(N×R_fc)的每一行对应于x_i ^T(1×R_fcL)PH步后的故障幅度；

(6.2)预测模型的建立

基于步骤6.1中的输入矩阵X和输出矩阵Y利用多变量最小二乘方法估计模型系数和建立经验预测模型：

Θ＝(X^TX)^-1X^TY   (40)

$\hat{Y}=\mathrm{X\Θ}$

其中Θ(R_fcL×R_fc)是回归系数矩阵，是得到的预测值；

(6.3)故障影响估计

由预测值与实际值Y得到的预测误差E(N×R_fc)为：

$E=Y-\hat{Y}---\left(41\right)$

在线应用时，通过结合当前数据和之前L步的历史数据，每一时刻的预测向量x_new ^T(1×R_fcL)可以表述为[f_c,new-L+1 ^T,f_c,new-L ^T,...,f_c,new ^T]；通过模型(18)，得到提前PH步的故障幅度预测值当实际值y_new可用时，预测误差e_new计算如下：

$\begin{array}{c}{\hat{y}}_{\mathrm{new}}={x_{\mathrm{new}}}^T\mathrm{\Θ}\\ e_{\mathrm{new}}=y_{\mathrm{new}}-{\hat{y}}_{\mathrm{new}}\end{array}---\left(42\right)$

基于预测的故障幅度，故障影响通过式(21)计算得到：

${{\hat{D}}_{\mathrm{new}}}^2={({\hat{y}}_{\mathrm{new}}-\stackrel{\‾}{f_c})}^T{\mathrm{\Λ}}^{-1}({\hat{y}}_{\mathrm{new}}-\stackrel{\‾}{f_c})---\left(43\right)$

定义事件报警限Elt_D，Elt_D略大于联合故障监测统计限Ctr_D，对于某一特定的PH，比较预测得到的与定义的事件报警限Elt_D，如果表明故障影响已经非常严重，如果表明故障在容许范围内；

(6.4)计算剩余时间RT

剩余时间RT是采取有效故障修复措施的时间，定义k^*为超出事件报警限Elt_D的时间；k^#为D²实际值超出事件报警限Elt_D的时间，得到剩余时间RT：

RT＝PH-(k^*-k^#)   (44)

k*-k^#表示监测统计量预测值相对于实际值的时延，用于评价故障预测结果的敏感性，由于提前PH步预测，如果PH大于时延，RT将大于0，有助于故障修复。

与现有技术相比，本发明的有益效果是：本发明所提出的方法可以及时预测故障演化趋势，在造成严重后果前及时报警，有效增加操作工程师的反应时间，为故障诊断和消除提供了可能。所提出的方法适用于一类具有慢时变动态演化特征的故障过程，且无需任何前提假设和过程先验知识。

附图说明

图1是不同预测步长下的过程故障监测统计量预测结果，(a)PH＝5时的训练数据(b)PH＝5时的测试数据(c)PH＝10时的训练数据(d)PH＝10时的测试数据；

图2是不同预测时域下PH的RT结果，(a)训练数据(b)测试数据；

图3是不同预测长度L下的预测性能(PH＝6)，(a)训练数据(b)测试数据；

图4是L＝15,PH＝6下的在线故障预测结果，(a)训练数据(b)测试数据。

具体实施方式

下面结合附图及具体实例，对本发明作进一步详细说明。

田纳西-伊斯曼化学工业过程是一个复杂的非线性过程，由伊斯曼化学品公司创建，其目的是为评价过程控制和监测方法提供一个真实的工业过程。这个过程共包含4个反应，有4种气体进料A、C、D、E，生成两种产品G、H，此外在进料中含有少量的惰性组分B和副产品F。

反应方程式如下：

A+C+D＝G

A+C+E＝H     (45)

A+E＝F

3D＝2F

过程包含41个测量变量及12个控制变量，变量如表1所示。

表1田纳西-伊斯曼过程测量变量表

序号	变量名称	序号	变量名称	序号	变量名称
						1	A组分进料流量	2	D组分进料流量	3	E组分进料流量
4	A和C组分进料流量	5	循环流量	6	反应器进料流量
						7	反应器压力	8	反应器液位	9	反应器温度
10	放空流量	11	产品分离温度	12	产品分离器液位
						13	产品分离器压力	14	产品分离器下部出料	15	汽提塔液位
16	汽提塔压力	17	汽提塔下部出料	18	汽提塔温度
						19	反应器冷却水出口温度	20	压缩机功率	21	汽提塔蒸汽流量
22	汽提塔冷却水出口温度	23	成分A(反应器进料)	24	成分B(反应器进料)
						25	成分C(反应器进料)	26	成分D(反应器进料)	27	成分E(反应器进料)
28	成分F(反应器进料)	29	成分A(放空气体分析)	30	成分B(放空气体分析)
						31	成分C(放空气体分析)	32	成分D(放空气体分析)	33	成分E(放空气体分析)
34	成分F(放空气体分析)	35	成分G(放空气体分析)	36	成分H(放空气体分析)
						37	成分D(产品分析)	38	成分E(产品分析)	39	成分F(产品分析)
40	成分G(产品分析)	41	成分H(产品分析)

本发明一种基于联合相对分析和自回归建模的故障建模及预测方法，包括以下步骤：

步骤1：过程分析与数据准备

设过程有J个测量变量和K个采样点，则可以形成一个K×J矩阵。其中测量变量为温度、速度压力、位移等过程运行的状态参数。X_r(N_r×J)和X_a(N_a×J)分别表示正常数据集和故障数据集，N_r和N_a分别代表正常数据集和故障数据集的样本个数。本实例的田纳西-伊斯曼过程有测量变量为41个，正常数据样本150个，故正常数据集为X_r(150×41)。此外，为验证所提方法的在线故障预测性能，考虑了该工业过程存在的一种缓慢偏移的动态特性故障，建模样本为150个。因此，故障数据集为X_a(150×41)。对数据集X_r和X_a按照式(2)进行中心化和归一化处理，使其均值为0标准差为1。

$x_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_j}{{\mathrm{\δ}}_j}---\left(46\right)$

其中x_k,j是数据集中第k行第j列的数据，k∈[1,K]是采样时刻，j∈[1,J]是变量索引，是矩阵第j个变量的均值，δ_j是第j个变量的标准差，满足以下关系：

${\stackrel{\‾}{x}}_j=\frac{1}{K}\underset{k}{\overset{K}{\mathrm{\Σ}}}{x}_{k,j}$

${\mathrm{\δ}}_j=\sqrt{\frac{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(x_{k,j}-{\stackrel{\‾}{x}}_j)}^2}{K-1}}---\left(47\right)$

步骤2：联合相对变化分析

(2.1)对正常数据集X_r进行主成分分析(PCA)获取监测模型：

$\left\{\begin{array}{c}{T}_r=X_r{P}_r\\ E_r=X_r{P}_r^e{P}_r^{\mathrm{eT}}\end{array}\right.---\left(48\right)$

其中T_r(N_r×R_r)和P_r(J×R_r)分别是主元成分和相应的主元负载；E_r(N_r×J)和是PCA分解得到的残差矩阵和相应的负载；R_r是通过累积解释方差(在本发明中选取为90％)保留的主元空间的主元个数。主元是正常数据集中最大的几个波动方向同时也是监测方向。是残差空间中保留的主元方向个数，有

根据公式(4)计算获取的主元得分T_r和残差矩阵E_r，计算两个监测统计指标：T²统计量和SPE统计量，计算公式如下：

${T_r}^2={(t_r-\stackrel{\‾}{t_r})}^T{{\mathrm{\Σ}}_r}^{-1}(t_r-\stackrel{\‾}{t_r})---\left(49\right)$

SPE_r＝e_r ^Te_r

其中，t_r(R×1)是来自T_r(N_r×R_r)的主元向量；是T_r(N_r×R_r)的均值向量，由于归一化的原因其值为0；Σ_r是T_r(N_r×R_r)的协方差矩阵；e_r(J×1)是E_r(N_r×J)的残差向量。

T²统计限和SPE统计限分别可以通过因子为α(α＝0.01)的F分布和χ²分布获取：

$T^2~\frac{R(N^2-1)}{N(N-R)}{F}_{R,N-R,\mathrm{\α}}---\left(27\right)$

SPE～gχ_h,α ²   (51)

其中g＝ν/2m、h＝2(m)²/ν；m是公式(5)中计算得到的SPE平均值，v是SPE的方差。

(2.2)利用正常监测模型检测故障数据

将故障数据X_a向P_r投影，获取故障数据的主元T_a；将X_a向投影，获得故障数据的残差E_a，如公式(8)所示：

$\left\{\begin{array}{c}{T}_a=X_a{P}_r\\ E_a=X_a{P}_r^e{P}_r^{\mathrm{eT}}=\underset{i=1}{\overset{R_e}{\mathrm{\Σ}}}{p}_{r,i}^e{p_{r,i}^e}^T\end{array}\right.E---\left(52\right)$

沿J个不同的监测方向，分别对比故障数据集和正常数据集的在主元空间和残差空间的波动变化：

$\begin{array}{c}{\mathrm{Ratio}}_i=\frac{\mathrm{var}\left(T_{a,i}\right)}{\mathrm{var}\left(T_{r,i}\right)}(i=\mathrm{1,2},...,R_r)\\ {\mathrm{\Δ}}_i={\left|\right|X_a{p}_{r,i}^e{p_{r,t}^e}^T\left|\right|}^2-{\left|\right|X_r{p}_{r,i}^e{p_{r,i}^e}^T\left|\right|}^2(i=\mathrm{1,2},...,R_r^e)\end{array}---\left(53\right)$

其中，var(*)表示距离中心点(从正常数据集中获取，此处为0)的主元方差；下标i代表矩阵的第i列。Ratio是由Ratio_i组成的R_r维向量。||||是欧氏距离，Δ是由Δ_i构成的维向量。

(2.3)故障相关方向排序

对Ratio和Δ进行排序，如果Ratio_i>1或是Δ_i>0，意味着在第i个方向上，故障数据比正常数据的波动大。因此，该方向对于监测统计量的超限起到了重要作用。找出所有Ratio_i>1和Δ_i>0的方向，这些方向在主元空间和残差空间中具有更大的波动变化，是可能的故障方向。

(2.4)更新监测统计量

移除步骤2.3中的故障方向，并重新计算正常数据和故障数据集新的T²和SPE监测统计限。定义正常数据新的监测统计限，并与故障数据的监测统计量进行对比。如果故障数据的监测统计量在正常数据的监测统计限内，认为关键故障方向已经得到移除；否则表明故障数据中仍存在一定的故障方向。

(2.5)迭代更新

重复步骤(2.2)-(2.4)直至所有故障方向都被移除。

将提取的故障方向构成数据矩阵和作为两个不同监测子空间的故障方向，保留的主元方向个数分别是和

步骤3：故障方向集成

沿和X_a的波动建模为：

$\begin{array}{c}{X}_{a,f}=X_a{P}_r^{*}{P}_r^{*T}\\ X_{a,f}^e=X_a{P}_r^{e*}{P}_r^{e*T}\end{array}---\left(54\right)$

可知，它们表示从正常数据到故障数据的增量变化，也就是受故障影响的信息。

对X_a,f和进行PCA分解，获得系统主要变化P_f(J×R_f)和其中，R_f和分别为主元个数，且有基于此，上述两个模型涵盖了受故障影响的方向，揭示了故障情况下的相对变化。这两个模型本质分别上是原PCA监测空间(主元空间P_r和残差空间)的一部分。集成到P_f和为统一的联合模型 $P_{\mathrm{fc}}(J\×R_{\mathrm{fc}})=[P_f,P_f^{*}](R_{\mathrm{fc}}=R_f^{*}+R_f),$ 易知联合模型的所有方向都是正交的。

步骤4：基于重构的故障幅度估计

(4.1)正常数据的重构

对于故障样本，监测统计量可能会超出之前设定的置信区域，产生报警信号。故障重构的目标是在消除故障影响后估计数据，也就是说移除T²和SPE这两个监测统计量的报警信号。由步骤3分解得到故障方向P_fc，不受故障影响的数据恢复为：

x^*＝x_f-P_fcf   (55)

其中，x_f是故障样本，x^*是从x_f中恢复得到的正常数据，f代表系统子空间的故障幅度；

(4.2)故障幅度的计算：

x^*的最优估计通过最小化x^*到监测子空间(主元空间或是残差空间)的距离来得到；故障幅度f计算如下：

$\begin{array}{c}f=\arg \min {\left|\right|\mathrm{\Ω}{x}^{*}\left|\right|}^2=\arg \min {\left|\right|\mathrm{\Ω}(x-P_{\mathrm{fc}}f)\left|\right|}^2=\arg \min {\left|\right|\tilde{x}-{\tilde{P}}_{\mathrm{fc}}f\left|\right|}^2\\ ={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}\end{array}---\left(56\right)$

其中，Ω表示向监测空间的投影算子，在主元空间Ω＝P_r(P_r ^TP_r)-1P_r ^T，在残差空间Ω＝I-P_r(P_r ^TP_r)^-1P_r ^T；I是J×J的单位矩阵；

当Ω＝P_r(P_r ^TP_r)^-1P_r时：

$\begin{array}{c}{\tilde{P}}_{\mathrm{fc}}=\mathrm{\Ω}{P}_{\mathrm{fc}}=P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r{P}_{\mathrm{fc}}=P_r{\left({P_r}^T{P}_r\right)}^{-1}{P_r}^T[P_f,P_f^{*}]=P_f\\ f={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({P_f}^T{P}_f\right)}^{-1}{P_f}^T{P}_r{\left({P_r}^T{P}_r\right)}^{-1}{P_r}^T={P_f}^{T}x\end{array}---\left(57\right)$

当Ω＝I-P_r(P_r ^TP_r)^-1P_r时

$\begin{array}{c}{\tilde{P}}_{\mathrm{fc}}=\mathrm{\Ω}{P}_{\mathrm{fc}}=(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)P_{\mathrm{fc}}=(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)[P_f,P_f^{*}]=P_f^{*}\\ f^{*}={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left(P_f^{*T}{P}_f^{*}\right)}^{-1}(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{p}_r)x=P_f^{*T}x\end{array}---\left(58\right)$

主元空间和残差空间的综合故障幅度f_c计算为：

$f_c=\left[\begin{array}{c}f\\ f^{*}\end{array}\right]=\left[\begin{array}{c}{P_f}^T\\ P_f^{*T}\end{array}\right]x={P_{\mathrm{fc}}}^{T}x---\left(59\right)$

步骤5：联合故障监测指标的确定

(5.1)联合故障监测统计指标D²实际值的计算

对于多变量系统故障，通过典型的马氏距离统计量计算得到故障影响：

$D^2={(f_c-\stackrel{\‾}{f_c})}^T{\mathrm{\Λ}}^{-1}(f_c-\stackrel{\‾}{f_c})---\left(60\right)$

其中，表示由正常数据计算得到的故障幅度平均值；Λ为对角矩阵，其对角元素是由正常数据计算得到的每个方向上的主元方差。对于一些非常接近于0的方差值赋予1以避免病态问题的出现。建立新的监测统计量D²沿不同方向和权重综合了这些波动。

(5.2)联合故障监测统计Ctr_D的建立

在正常操作条件下，假设过程数据服从多变量正态分布，当样本数N足够大时，可以利用χ²分布计算D²的置信限：

$D^2\≤{\mathrm{Ctr}}_D\≡g{\mathrm{\χ}}_{h,a}^2---\left(61\right)$

其中，g＝ν/2m、h＝2(m)²/ν,m是在正常训练数据下根据公式(16)计算得到的D²平均值，ν是相应的方差。

步骤6：故障预测

(6.1)预测模型输入输出数据整理

联合故障幅度构成的向量随着过程的进行不断变化，并由于向量间的相关性形成多变量时间序列。这些向量包含丰富的故障幅度变动信息和数据相关信息。通过与当前故障幅度的关系，潜在的故障演化特性可以由和过去、将来故障幅度关系来反应。因此，故障动态特性(多变量时间序列相关关系)建模过程可以看成关于故障幅度预测值(输出)和历史故障测量值(输入)的过程。记L代表输入矩阵的时间长度，在每一个采样时刻，具有一个R_fc维的故障幅度向量。基于以上分析，输入矩阵X(N×R_fcL)可以表示为：

$X(N\×R_{\mathrm{fc}}L)=\left[\begin{array}{c}{f_{c,1}}^T,{f_{c,2}}^T,...,{f_{c,L}}^T\\ {f_{c,2}}^T,{f_{c,3}}^T,...,{f_{c,L+1}}^T\\ .\\ .\\ .\\ {f_{c,N}}^T,{f_{c,N+1}}^T,...,{f_{c,L+N-1}}^T\end{array}\right]---\left(62\right)$

其中，N为输入矩阵的样本数，每一个x_i＝[f_c,i ^T,f_c,i+1 ^T,...,f_c,i+L-1 ^T]是一个R_fcL维的向量，包含了从时刻i到i+L-1的故障幅度信息；

预测时域PH表示从当前时刻对PH步后的故障幅度进行预测，如PH＝6表示对6个时刻后的故障值进行预测，输出矩阵Y(N×R_fc)的每一行由对应于x_i ^T(1×R_fcL)PH步后的故障幅度构成；

(6.2)预测模型的建立

利用多变量最小二乘方法估计模型系数和建立经验预测模型：

Θ＝(X^TX)^-1X^TY   (63)

$\hat{Y}=\mathrm{X\Θ}$

其中Θ(R_fcL×R_fc)是回归系数矩阵，是得到的预测值。

(6.3)故障影响估计

由预测值与实际值Y的得到的预测误差E(N×R_fc)为：

$E=Y-\hat{Y}---\left(64\right)$

在线应用时，通过结合当前数据和之前L步的历史数据，每一时刻的预测向量_xnew ^T(1×R_fcL)可以表述为[f_c,new-L+1 ^T,f_c,new-L ^T,...,f_c,new ^T]。通过模型(19)，得到提前PH步的故障幅度预测值当实际值y_new可用时，预测误差e_new计算如下：

$\begin{array}{c}{\hat{y}}_{\mathrm{new}}={x_{\mathrm{new}}}^T\mathrm{\Θ}\\ e_{\mathrm{new}}=y_{\mathrm{new}}-{\hat{y}}_{\mathrm{new}}\end{array}---\left(65\right)$

基于预测的故障幅度，故障影响通过式(22)计算得到：

${{\hat{D}}_{\mathrm{new}}}^2={({\hat{y}}_{\mathrm{new}}-\stackrel{\‾}{f_c})}^T{\mathrm{\Λ}}^{-1}({\hat{y}}_{\mathrm{new}}-\stackrel{\‾}{f_c})---\left(66\right)$

定义事件报警限Elt_D，Elt_D略大于联合故障监测统计限Ctr_D。随着故障过程进行，故障幅度持续增大。对于某一特定的PH，比较预测得到的与定义的事件监测限Elt_D。如果时，表明故障影响已经非常严重；如果表明故障在容许范围内。

图1的(a)-(d)分别展示了PH＝5和PH＝10时的故障预测结果，图中同时画出了原始的D²监测统计量及其预测值。从图中可以看出，D²预测值和其实际值之间具有很好的一致性，表明预测模型具有较高的精度，可以有效指示过程故障的发生。在D²预测的波峰处，会有较大的预测误差。同时，可以看出随着PH的增大，事件报警也会相应的增大。

(6.4)计算剩余时间RT

在故障影响变得严重之前，获知故障何时到达某一程度十分重要。剩余时间RT是采取有效故障修复措施的时间。定义k^*为超出事件报警限Elt_D的时间；k^#为D²实际值超出事件报警限Elt_D的时间。得到剩余时间RT：

RT＝PH-(k^*-k^#)   (67)

k*-k^#表示监测统计量预测值相对于实际值的时延，用于评价故障预测结果的敏感性。由于提前PH步预测，如果PH大于时延，RT将大于0，有助于故障修复。综上，RT越大，报警信息越可以提早预测，工程师会有更多的故障处理时间。当然，PH的选择会对RT产生一定的影响。一般来说，PH偏大时，预测精度将降低，时延也将变大；相反，PH偏小时，预测精度将升高，时延减小。

图2展示了预测步长PH对剩余时间RT的影响关系。对于反应数据拟合能力的训练数据，随着PH的增加，RT也将增加，表明PH的增幅快于时延。对于反应数据泛化能力的测试数据，随着PH的增加，除了少部分PH值外，RT基本保持为3不变。这可能是时延具有和PH相同的增幅。

图3给出了剩余寿命RT和不同预测长度L之间的关系。对于训练数据，较大的L具有较小的拟合误差；但是对于测试数据，较大的L的拟合误差也较大。同时，训练数据的RT受L变化影响较小，而测试数据受影响较大。当L处于10到13时，RT的值小于0，表明时延比PH值大。L＝15时，RT为45，比PH大很多。从图4(b)也可以看出，测试数据的D²预测值在第113个采样点处超限，而其实际值直到第152个采样点处才报警。

在采取有效的校正措施之前，首先需要确定预测结果是否可靠。假设式(20)得到的数据遵从多变量正态分布，计算统计量SPE基于之前的预测结果：

SPE＝e^Te(68)

其中，e(R_fc×1)是从式(20)计算得到的残差矩阵E的行向量。

SPE的100(1-α)％控制限服从权重卡方分布，定义如下：

SPE～gχ_h,α ²   (69)

其中，g＝ν/2m、h＝2(m)²/ν；m是式(22)中所有SPE的平均值，ν是SPE的方差。

SPE统计量提供了一种有效的预测结果评估手段，如果SPE统计量超出式(25)的置信限，意味着当前预测得到的故障幅度是不可信的，对应计算出的D²预测值应该舍弃。

本发明一种基于联合相对分析和自回归建模的故障预测方法，阐述了故障预测中的两个重要问题：1)如何定量地衡量这些故障过程变化对监测性能的影响；2)预测故障影响将于何时超越置信限。首先，在主成分分析的基础上，提出了一种基于联合相对分析算法分解故障影响，确定故障方向。随后，根据确定的故障方向，基于重构技术估计故障幅度，进行正常部分数据的恢复。定义新监测统计量D²涵盖正常的数据波动，建立置信限。最后，建立了基于新监测指标的自回归模型用于在线监测统计量D²的提前预测，实现故障情况的提前报警。本发明简便有效，不依赖先验过程知识和假设，在典型的复杂工业过程田纳西-伊斯曼化工过程中的实验效果良好。故障预测结果对后续的故障诊断和修复具有重要意义，有助于过程工程师对过程运行状态做出及时判断，从而保证工业生产的安全可靠运行和产品的高质量追求。

应该理解，本发明并不局限于上述具体实施例的田纳西-伊斯曼过程，凡是熟悉本领域的技术人员在不违背本发明精神的前提下还可做出等同变形或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。

Claims (1)

1.一种基于联合相对变化分析和自回归模型的故障预测方法，其特征在于，该方法包括以下步骤：

步骤1：过程分析与数据准备

设过程有J个测量变量和K个采样点，则可以形成一个K×J矩阵；用X_r(N_r×J)和X_a(N_a×J)分别描述正常数据集和故障数据集；N_r和N_a分别代表正常数据集和故障数据集的样本个数；对数据集X_r和X_a按式(1)进行中心化和归一化处理，使其均值为0，标准差为1；

$x_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_j}{{\mathrm{\δ}}_j}---\left(1\right)$

其中x_k,j是数据集中第k行第j列的数据，k∈[1,K]是采样时刻，j∈[1,J]是变量索引；是矩阵第j个变量的均值，δ_j是第j个变量的标准差，满足以下关系：

${\stackrel{\‾}{x}}_j=\frac{1}{K}\underset{k}{\overset{K}{\mathrm{\Σ}}}{x}_{k,j}$

${\mathrm{\δ}}_j=\sqrt{\frac{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(x_{k,j}-{\stackrel{\‾}{x}}_j)}^2}{K-1}}---\left(2\right)$

步骤2：联合相对变化分析

(2.1)对正常数据集X_r进行主成分分析(PCA)获取监测模型：

$\left\{\begin{array}{c}{T}_r=X_r{P}_r\\ E_r=X_r{P}_r^e{P}_r^{e^T}\end{array}\right.---\left(3\right)$

其中T_r(N_r×R_r)和P_r(J×R_r)分别是主元成分和相应的主元负载；E_r(N_r×J)和是PCA分解得到的残差矩阵和相应的负载；R_r是通过累积解释方差保留的主元空间的主元个数；是残差空间中保留的主元方向个数，有

根据公式(3)计算获取的主元得分T_r和残差矩阵E_r，计算两个监测统计指标：T²统计量和SPE统计量，计算公式如下：

${T_r}^2={(t_r-{\stackrel{\‾}{t}}_r)}^T{{\mathrm{\Σ}}_r}^{-1}(t_r-{\stackrel{\‾}{t}}_r)---\left(4\right)$

SPE_r＝e_r ^Te_r

其中，t_r(R×1)是来自T_r(N_r×R_r)的主元向量；是T_r(N_r×R_r)的均值向量，由于归一化的原因其值为0；Σ_r是T_r(N_r×R_r)的协方差矩阵；e_r(J×1)是E_r(N_r×J)的残差向量；

T²统计限和SPE统计限分别可以通过因子为α的F分布和χ²分布获取：

$T^2~\frac{R(N^2-1)}{N(N-R)}{F}_{R,N-R,\mathrm{\α}}---\left(5\right)$

SPE～gχ_h,α ²                         (6)

其中g＝ν/2m、h＝2(m)²/ν；m是公式(4)中计算得到的SPE平均值，v是SPE的方差；

(2.2)利用正常监测模型检测故障数据

将故障数据X_a向P_r投影，获取故障数据的主元T_a；将X_a向投影，获得故障数据的残差E_a，如公式(7)所示：

$\left\{\begin{array}{c}{T}_a=X_a{P}_r\\ E_a=X_a{P}_r^e{P_r^e}^T=\underset{i=1}{\overset{R_e}{\mathrm{\Σ}}}{P}_{r,i}^e{p_{r,i}^e}^T\end{array}\right.---\left(7\right)$

沿J个不同的监测方向，分别对比故障数据集和正常数据集在主元空间和残差空间的波动变化：

${\mathrm{Ratio}}_i=\frac{\mathrm{var}\left(T_{a,i}\right)}{\mathrm{var}\left(T_{r,i}\right)}(i=\mathrm{1,2},...,R_r)---\left(8\right)$

${\mathrm{\Δ}}_i={\left|\right|X_a{p}_{r,i}^e{p_{r,i}^e}^T\left|\right|}^2-{\left|\right|X_r{p}_{r,i}^e{p_{r,i}^e}^T\left|\right|}^2(i=\mathrm{1,2},...,R_r^e)$

其中，var(^*)表示距离中心点(从正常数据集中获取，此处为0)的主元方差；下标i代表矩阵的第i列；Ratio是由Ratio_i组成的R_r维向量；‖‖是欧氏距离，Δ是由Δ_i构成的维向量；

(2.3)故障相关方向排序

对Ratio和Δ进行排序，找出所有Ratio_i>1和Δ_i>0的方向，这些方向在主元空间和残差空间中具有更大的波动变化，是可能的故障方向；

(2.4)更新监测统计量

移除步骤2.3找到的故障方向，并重新计算正常数据集和故障数据集新的T²和SPE监测统计量；定义正常数据集新的监测统计限，并与故障数据的监测统计量进行对比；如果故障数据的监测统计量在正常数据的监测统计限内，认为关键故障方向已经得到移除；否则表明故障数据中仍存在故障方向；

(2.5)迭代更新

重复步骤(2.2)-(2.4)直至所有故障方向都被移除；

将提取的故障方向构成数据矩阵和作为两个不同监测子空间的故障方向，保留的主元方向个数分别是和

步骤3：故障方向集成

沿和X_a的波动可以建模为：

$X_{a,f}=X_a{P}_r^{*}{P}_r^{*T}$

                                               (9)

$X_{a,f}^e=X_a{P}_r^{e*}{P}_r^{e*T}$

对X_a,f和进行PCA分解，获得系统主要变化P_f(J×R_f)和其中，R_f和分别为主元个数，且有集成P_f和为统一的联合模型 $P_{\mathrm{fc}}(J\×R_{\mathrm{fc}})=[P_f,P_f^{*}](R_{\mathrm{fc}}=P_f^{*}+R_f);$

步骤4：基于重构的故障幅度估计

(4.1)正常数据的重构

由步骤3分解得到故障方向P_fc，不受故障影响的数据能够恢复为：

x^*＝x_f-P_fcf                      (10)

其中，x_f是故障样本，x^*是从x_f中恢复得到的正常数据，f代表系统子空间的故障幅度；

(4.2)故障幅度的计算

x^*的最优估计通过最小化x^*到监测子空间(主元空间或是残差空间)的距离来得到；故障幅度f计算如下：

$\begin{array}{c}f=\arg \min {\left|\right|\mathrm{\Ω}{x}^{*}\left|\right|}^2=\arg \min {\left|\right|\mathrm{\Ω}(x-P_{\mathrm{fc}}f)\left|\right|}^2=\arg \min {\left|\right|}^2=\arg \min {\left|\right|\tilde{x}-{\tilde{P}}_{\mathrm{fc}}f\left|\right|}^2\\ ={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}\end{array}---\left(11\right)$

其中，Ω表示监测空间的投影算子，在主元空间Ω＝P_r(P_r ^TP_r)^-1P_r ^T，在残差空间Ω＝I-P_r(P_r ^TP_r)^-1P_r ^T；I是J×J的单位矩阵；

当Ω＝P_r(P_r ^TP_r)^-1P_r时：

${\tilde{P}}_{\mathrm{fc}}=\mathrm{\Ω}{P}_{\mathrm{fc}}=P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r{P}_{\mathrm{fc}}=P_r{\left({P_r}^T{P}_r\right)}^{-1}{P_r}^T[P_f,P_f^{*}]=P_f$

                                              (12)

$f={\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({P_f}^T{P}_f\right)}^{-1}{P_f}^T{P}_r{\left({P_r}^T{P}_r\right)}^{-1}{P_r}^{T}x={P_f}^{T}x$

当Ω＝I-P_r(P_r ^TP_r)^-1P_r时

${\tilde{P}}_{\mathrm{fc}}=\mathrm{\Ω}{P}_{\mathrm{fc}}=(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)P_{\mathrm{fc}}=(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)[P_f,P_f^{*}]=P_f^{*}$

                                          (13)

$f^{*}{\left({{\tilde{P}}_{\mathrm{fc}}}^T{\tilde{P}}_{\mathrm{fc}}\right)}^{-1}{{\tilde{P}}_{\mathrm{fc}}}^T\tilde{x}={\left({P_f^{*}}^T{P}_f^{*}\right)}^{-1}{P_f^{*}}^T(I-P_r{\left({P_r}^T{P}_r\right)}^{-1}{P}_r)x={P_f^{*}}^{T}x$

主元空间和残差空间的综合故障幅度f_c计算为：

$f_c=\left[\begin{array}{c}f\\ f^{*}\end{array}\right]=\left[\begin{array}{c}{P_f}^T\\ {P_f^{*}}^T\end{array}\right]x={P_{\mathrm{fc}}}^{T}x---\left(14\right)$

步骤5：联合故障监测指标的确定

(5.1)联合故障监测统计指标D²实际值的计算

对于多变量系统故障，通过典型的马氏距离统计量计算得到故障影响：

$D^2={(f_c-{\stackrel{\‾}{f}}_c)}^T{\mathrm{\Λ}}^{-1}(f_c-{\stackrel{\‾}{f}}_c)---\left(15\right)$

其中，表示由正常数据计算得到的故障幅度平均值；Λ为对角矩阵，其对角元素是由正常数据计算得到的每个方向上的主元方差，对于一些非常接近于0的方差值赋予1；

(5.2)联合故障监测统计限Ctr_D的建立

在正常操作条件下，假设过程数据服从多变量正态分布，当样本数N足够大时，可以利用χ²分布计算D²的置信限：

$D^2\≤{\mathrm{Ctr}}_D\≡g{\mathrm{\χ}}_{h,a}^2---\left(16\right)$

其中，g＝ν/2m、h＝2(m)²/ν,m是在正常训练数据下根据公式(15)计算得到的D²平均值，ν是相应的方差；

步骤6：故障预测

(6.1)预测模型输入输出数据整理

记L代表输入矩阵的时间长度，在每一个采样时刻，具有一个R_fc维的故障幅度向量，输入矩阵X(N×R_fcL)表示为：

$X(N\×R_{\mathrm{fc}}L)=\left[\begin{array}{c}{f_{c,1}}^T,{f_{c,2}}^T,...,{f_{c,L}}^T\\ {f_{c,2}}^T,{f_{c,3}}^T,...,{f_{c,L+1}}^T\\ .\\ .\\ .\\ {f_{c,N}}^T,{f_{c,N+1}}^T,...,{f_{c,L+N-1}}^T\end{array}\right]---\left(17\right)$

其中，N为输入矩阵的样本数，每一个x_i＝[f_c，i ^T，f_c，i+1 ^T，…，f_c，i+L-1 ^T]是一个R_fcL维的向量，包含了从时刻i到i+L-1的故障幅度信息；

预测时域PH表示从当前时刻对PH步后的故障幅度进行预测；输出矩阵Y(N×R_fc)的每一行对应于x_i ^T(1×R_fcL)PH步后的故障幅度；

(6.2)预测模型的建立

基于步骤6.1中的输入矩阵X和输出矩阵Y利用多变量最小二乘方法估计模型系数和建立经验预测模型：

Θ＝(X^TX)^-1X^TY                 (18)

$\hat{Y}=\mathrm{X\Θ}$

其中Θ(R_fcL×R_fc)是回归系数矩阵，是得到的预测值；

(6.3)故障影响估计

由预测值与实际值Y得到的预测误差E(N×R_fc)为：

$E=Y-\hat{Y}$

在线应用时，通过结合当前数据和之前L步的历史数据，每一时刻的预测向量x_new ^T(1×R_fcL)可以表述为[f_c,new-L+1 ^T,f_c,new-L ^T,...,f_c,new ^T]；通过模型(18)，得到提前PH步的故障幅度预测值当实际值y_new可用时，预测误差e_new计算如下：

${\hat{y}}_{\mathrm{new}}={x_{\mathrm{new}}}^T\mathrm{\Θ}$

                                    (20)

$e_{\mathrm{new}}=y_{\mathrm{new}}-{\hat{y}}_{\mathrm{new}}$

基于预测的故障幅度，故障影响通过式(21)计算得到：

${{\hat{D}}_{\mathrm{new}}}^2={({\hat{y}}_{\mathrm{new}}-{\stackrel{\‾}{f}}_c)}^T{\mathrm{\Λ}}^{-1}({\hat{y}}_{\mathrm{new}}-{\stackrel{\‾}{f}}_c)---\left(21\right)$

定义事件报警限Elt_D，Elt_D略大于联合故障监测统计限Ctr_D，对于某一特定的PH，比较预测得到的与定义的事件报警限Elt_D，如果表明故障影响已经非常严重，如果表明故障在容许范围内；

(6.4)计算剩余时间RT

剩余时间RT是采取有效故障修复措施的时间，定义k^*为超出事件报警限Elt_D的时间；k^#为D²实际值超出事件报警限Elt_D的时间，得到剩余时间RT：

RT＝PH-(k^*-k^#)                   (22)

k^*-k^#表示监测统计量预测值相对于实际值的时延，用于评价故障预测结果的敏感性，由于提前PH步预测，如果PH大于时延，RT将大于0，有助于故障修复。