一种基于嵌套迭代费舍尔判别分析的故障变量隔离方法

摘要

本发明公开了一种基于嵌套迭代费舍尔判别分析的故障变量隔离方法，该方法充分发掘了过程数据所包含的潜在信息，能有效区分对故障起不同程度作用的变量，从而根据不同类别分别建立判别模型，改善了模型精度，大大提高了在线故障诊断的性能，增强了实际在线故障诊断的可靠性和可信度，有助于工业工程师对故障进行准确的修复，从而保证实际生产的安全可靠运行和产品的高质量追求。

说明书

技术领域

本发明属于化工过程统计监测领域，特别是涉及一种基于嵌套迭代费舍尔判 别分析的故障变量隔离方法。

背景技术

作为工业生产中一种重要的生产方式，化工过程与人们的生活息息相关，已 被广泛应用于冶金、炼油、造纸、制革等领域。如何保证化工过程生产安全，提 高产品质量和经济效益是人们关注的焦点。随着化工过程日益复杂，在线故障检 测及诊断也越来越重要。故障诊断是指在检测到故障发生后，进一步判断发生了 哪种故障。而通过故障变量隔离则可以很好的识别发生的故障类别，保障生产的 安全可靠运行以及产品的高质量，从而可以避免重大安全事故，减少人员伤亡以 及提高经济效益。随着技术的发展，工业现场可以获得越来越多的数据，基于数 据的在线故障诊断策略越来越受到研究人员和现场工程师的青睐。化工过程的统 计建模、在线监测、故障诊断及质量预测已成为广泛的研究课题。

前人对此已经作了相应的研究与探讨，基于不同的角度提出了相应的在线故 障诊断办法。归纳起来有如下几种：基于重构指标的方法、基于相似度计算的模 式匹配方法、基于模型的改进隔离方法以及降维的方法。其中，降维的方法如主 元分析、费舍尔判别分析能有效处理高维度、高相关性的数据，它们通过构造潜 变量将高维度原始测量数据投影到低维度的监测空间，从而提高了故障诊断精 度，被广泛应用于工业过程在线故障诊断。总体来说，上述几种基于降维的故障 诊断方法各有各的适用场合与优缺点。相对而言，费舍尔判别分析方法侧重于区 分具有不同特性的数据，在故障诊断方面更有优势。但是，基于传统费舍尔判别 分析的在线故障诊断方法应用于实际化工过程时存在三方面的问题：首先，化工 过程数据往往是高度耦合的，这可能导致类内散布矩阵是奇异的，从而无法进行 奇异值分解提取过程数据潜在信息。其次，由于类间散布矩阵奇异，有可能导致 判别成分的个数小于类别的个数，从而使得散布矩阵无法提供充足的过程信息。 最后，在每个类中，所提取的判别成分是线性相关的，这就导致所提取的过程信 息冗余。针对传统费舍尔判别分析方法的不足，研究人员提出了一系列的改进方 法。总的来说，这些方法均采用两步法解决高耦合数据带来的散步矩阵奇异性问 题，其关键在于如何在进行费舍尔判别分析之前进行数据降维。然而，前人的方 法在解决奇异性问题时均存在一定程度的问题，如数据压缩不当而导致无法提取 出过程数据的关键潜在信息，或者过程重要信息缺失等等，从而导致故障诊断精 度欠缺。为了克服上述问题，研究人员提出一种迭代嵌套的费舍尔判别分析方法， 提高了故障诊断性能。然而，该方法没有考虑到不同的变量对故障的影响程度不 同，有些变量对故障影响很大而另一些变量与故障不甚相关，使用所有变量数据 建立诊断模型，降低了模型精度，影响了在线故障诊断的性能。

本发明的内容深入考虑了化工过程的复杂性、数据的高维度高耦合性以及不 同变量对故障作用的程度不同，提出了一种基于嵌套迭代费舍尔判别分析的故障 变量隔离方法。该方法充分发掘了故障相对于正常数据的变化情况，区分了对故 障有重要影响的故障变量和没有影响的普通变量，并分别建立了故障诊断模型， 大大提高了在线故障诊断性能。到目前为止，尚未见到与本发明相关的研究报道。

发明内容

本发明的目的在于针对现有针对化工生产过程的在线故障诊断技术的不足， 提供一种基于嵌套迭代费舍尔判别分析的故障变量隔离方法。

本发明的目的是通过以下技术方案来实现的：一种基于嵌套迭代费舍尔判别 分析的故障变量隔离方法，该方法包括以下步骤：

(1)获取过程分析数据：设一个化工生产过程具有J个测量变量和操作变量， 则每一次采样可得到一个1×J的向量，采样K次后得到的数据表述为一个二维矩 阵X(K×J)，所述测量变量为运行过程中可被测量的状态参数，包括流量、温度、 速率；所述操作变量包括给料量、阀门开度；分别获取正常数据二维矩阵X_n(K×J) 和故障数据二维矩阵X_f,m(K×J)，其中，下标n表示正常数据，下标f表示故障数 据，m表示故障的类别；将正常数据和故障数据统一标示为X_i(K×J)，其中下标 i表示数据的类别；

(2)选取正常数据样本和一类故障数据样本作为总样本其中， 由X_i(i＝1,2)从上到下排列组成；

(3)数据准备：分别计算总样本均值向量每类样本均值向量总类内 散布矩阵S_w和类间的散布矩阵S_b，计算公式如下：

$\left(\begin{array}{c}{S}_i=\underset{x_i\in X_i}{\Sigma }(x_i-{\bar{x}}_i){(x_i-{\bar{x}}_i)}^T\\ S_w=\underset{i=1}{\overset{2}{\Sigma }}{S}_i\\ S_b=\underset{x_i\in X_i}{\Sigma }{K}_i({\bar{x}}_i-\bar{x}){({\bar{x}}_i-\bar{x})}^T\end{array}\right)---\left(1\right)$

其中，S_i是每个类的散布矩阵；

(4)提取初始判别成分，该步骤由以下子步骤来实现：

(4.1)最大化类间离散度：求取使类间离散度最大的权重向量w，即相当于 求取类间散布矩阵S_b的最大特征值所对应的特征向量w，所述类间离散度为 w^TS_bw，获取w后，按公式(2)求取相应的总样本初始判别成分t；

$t=\bar{X}w---\left(2\right)$

其中，是减均值中心化后的总样本，那么对于每一类样本，其所对应的类 判别成分为可知，t由t_i从上到下依次排列构成；

(4.2)数据压缩：对减均值中心化后的总样本根据下式进行数据压缩：

$\left(\begin{array}{c}{p}^T={\left(t^{T}t\right)}^{-1}{t}^T\bar{X}\\ \bar{E}=\bar{X}-t{p}^T\end{array}\right)---\left(3\right)$

其中：p表示总样本的负载向量，表示总样本中与t无关的残差；

同理，对于每类样本都可以通过公式(4)得到与t_i无关的残差且由从上到下排列组成：

${\bar{E}}_i={\bar{X}}_i-t_i{p}^T={\bar{X}}_i-{\bar{X}}_{i}w{p}^T---\left(4\right)$

最后，用上述数据压缩关系wp^T更新每一个类的信息，以保证判别成分的正 交性：

E_i＝X_i-X_iwp^T                 (5)

(4.3)迭代更新过程数据

(4.3.1)用步骤(4.2)中获得的E_i代替步骤(3)中的X_i，按步骤(3)重 新计算总样本均值每类样本均值向量总类内散布矩阵S_w和类间的散布矩 阵S_b，按步骤(4.1)、(4.2)再次提取初始判别成分；

(4.3.2)重复步骤(4.3.1)直到所提取的初始判别成分的个数等于S_w的阶数 N；那么，同时可以得到由权重向量w组成的权重矩阵W(J×N)和相应的负载向 量p组成的负载矩阵P(J×N)、总样本初始判别成分t组成的总样本的初始判别成 分矩阵其中，T由T_i按从上至下排列构成，T_i是每个类的判别成分 矩阵；最后，求取初始判别成分的系数矩阵R＝W(P^TW)^-1，且T和T_i可直接由系 数矩阵根据公式(6)求出：

$\left(\begin{array}{c}T=\bar{X}R=\bar{X}W{\left(P^{T}W\right)}^{-1}\\ T_i={\bar{X}}_{i}R={\bar{X}}_{i}W{\left(P^{T}W\right)}^{-1}\end{array}\right)---\left(6\right)$

(5)提取最终判别成分，该步骤通过以下子步骤来实现：

(5.1)过程数据预处理：使用X_iR代替每类初始数据集合X_i，按步骤(3) 重新计算每类样本均值总样本均值总类内散布矩阵S_w^*以及类间散布矩阵 S_b^*；

(5.2)确定最终判别成分：最终判别成分通过以下步骤来确定：

(5.2.1)求取最优判别成分方向向量w^*，使得类间散布矩阵与类内散布矩阵 的比值J(θ)最大；其中，w^*则可通过公式(7)求 取矩阵S_w^*-1S_b^*最大特征值所对应的特征向量得到：

S_w^*-1S_b^*w^*＝λw^*                 (7)

(5.2.2)求取每类的最终判别成分向量t_i^*：

$\left(\begin{array}{c}{t_i}^{*}=X_{i}R{w}^{*}=X_i\theta \\ \theta =R{w}^{*}\end{array}\right)---\left(8\right)$

(5.2.3)将t_i^*从上至下依次排列构成总样本的最终判别成分向量

(5.3)压缩过程数据：为了保证每类样本的判别成分之间是正交的，进行如 下处理：

$\left(\begin{array}{c}{p_i}^{*T}={\left({t_i}^{*T}{t_i}^{*}\right)}^{-1}{t_i}^{*T}{X}_i\\ {E_i}^{*}=X_i-{t_i}^{*}{p_i}^{*T}\end{array}\right)---\left(9\right)$

其中，是每类的负载向量，E_i^*是与t_i^*无关的残差；

(6)迭代更新过程数据，该步骤包括以下子步骤：

(6.1)用步骤(5.3)中E_i^*代替步骤(3)中的X_i，按步骤(3)重新计算每 类样本均值、总样本均值，总类内散布矩阵S_w^*以及类间散布矩阵S_b^*，按步骤(4) 和步骤(5)再次提取最终判别成分向量t_i^*；

(6.2)重复步骤(6.1)直至获得足够的最终判别成分t_i^*并构成最终的判别 成分矩阵T_i^*，T_i^*所保留的最终判别成分个数为R，所述R通过交叉检验的方法确 定；相应的，同时可以获得权重矩阵Θ(J×R)和负载矩阵P_i^*(J×R)；其中，Θ(J×R) 和P_i^*(J×R)分别由θ(J×1)和p_i^*(J×1)构成；

(6.3)求取最终系数矩阵R_i^*(J×R)：

R_i^*＝Θ(P_i^*TΘ)^-1                (10)

至此，步骤(2)中选取的该类故障最终系数矩阵被求取出来；

(7)选取正常数据和另一类故障数据作为总样本，重复步骤(4)-(6)， 获得该类故障样本的最终系数矩阵；

(8)重复步骤(7)直到M类故障的最终系数矩阵都被求 取出来；

(9)故障数据衡量；该步骤由以下子步骤完成：

(9.1)求取每类故障X_f,m及相应正常数据X_n的类内成分矩阵T_f,m^*和T_n,m^*：

$\left(\begin{array}{c}{T_{f,m}}^{*}=X_{f,m}{R_{f,m}}^{*}\\ {T_{n,m}}^{*}=X_n{R_{f,m}}^{*}\end{array}\right)---\left(11\right)$

(9.2)求取每类故障X_f,m及相应正常数据X_n的D²指标：

$\left(\begin{array}{c}{t}_{n,m,k}={x_{n,k}}^T{R_{f,m}}^{*}\\ {D_{n,m,k}}^2={(t_{n,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}(t_{n,m,k}-{\bar{t}}_{n,m})\\ t_{f,m,k}={x_{f,k}}^T{R_{f,m}}^{*}\\ {D_{f,m,k}}^2={(t_{f,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}(t_{f,m,k}-{\bar{t}}_{n,m})\end{array}\right)---\left(12\right)$

其中，x_n,k是X_n中一个样本，t_n,m,k是x_n,k的判别成分，是公式(11)中T_n,m^*的均值向量，Σ_n则表示T_n,m^*的协方差矩阵，如果方差过小而趋近于0，那么它们 将被置为1；同理，x_f,m是X_f,m中一个样本，t_n,m,k是x_f,m的判别成分；

(9.2)建立正常数据基于D²指标的控制限：由于过程数据服从多变量正态 分布，那么可知D²统计量服从带权重的χ²分布，则可依据带权重的χ²分布建立 D²统计量的控制限

(10)比较每类故障数据的D²指标与相应的正常数据控制限若没有 超限报警则说明该类故障与正常数据相似，无需进行以下步骤；若有超限报警， 则执行以下步骤辨识和分离故障变量；

(11)变量贡献度衡量，该步骤包括以下子步骤；

(11.1)按公式(13)计算每类故障及相应正常数据的变量贡献度指标

$\left(\begin{array}{c}{t}_{·,m,k}={x_{·,k}}^T{R_{f,m}}^{*}\\ C_{D^2,·,k,j}={(t_{·,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}{r}_{m,j}(x_{·,k,j}-{\bar{x}}_{n,j})\end{array}\right)---\left(13\right)$

其中，下标·表示正常数据或者故障数据，x_·,k表示相应正常或者故障数据中 的一个样本，t_·,m,k表示相应样本x_·,k的判别成分，x_·,k,j是相应数据中的第k个样本 的第j个变量，是正常数据中第j个变量的均值，r_m,j是第j个变量的权重系数；

(11.2)因为正常数据服从一定的分布规律，根据分布规律确定正常数据贡 献度指标的控制限；

(11.3)按公式(14)计算故障样本相对于正常样本的变量贡献度的比值

${\mathrm{RC}}_{D^2,f,k,j}=\frac{C_{D^2,f,k,j}}{\mathrm{ctr}\left(C_{D^2,n,j}\right)}=\frac{{(t_{·,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}{r}_{m,j}(x_{·,k,j}-{\bar{x}}_{n,j})}{\mathrm{ctr}\left(C_{D^2,n,j}\right)}---\left(14\right)$

其中，表示正常数据的控制限；

(11.4)对于每个变量，将通过公式(14)所计算得到的M类故障对于正常 样本的变量贡献度的比值构成向量并计算的均值

(12)故障变量选择；将步骤(11.4)所得到的排序，选择最 大所对应的变量作为故障变量，将该变量从故障数据和正常数据 中移除，得到新的故障数据和正常数据

(13)模型更新；对更新后的故障数据和正常数据执行步骤(3)～ (8)，求取新的最终系数矩阵，即代表故障方向；

(14)基于更新后的最终系数矩阵，重复执行步骤(9)～(12)直到所有故 障变量都被选取出来；对于每类故障数据，故障变量构成矩阵其 中，J_f,m表示故障变量的个数，N_f表示故障样本的个数；则每类故障被分为两部 分数据：对故障有重大影响的以及对故障没有影响的其中，J_n,m＝J-J_f,m；同样，对于正常数据也被分为两部分和 其中，N_n表示正常数据的个数，J_n,m表示普通变量的个数；那么， 对于每一对正常数据X_n和故障数据X_f,m，可以得到由故障变量构成的矩阵 和以及由剩余的普通变量构成的矩阵和

(15)建立故障诊断模型，该步骤包括以下子步骤：

(15.1)用和分别代替步骤(2)中获取的正常样 本数据和故障数据，重复步骤(3)至(8)提取出新的最终系数矩阵即表 示有重要区分度的故障方向；

(15.2)计算每类故障子数据的判别成分和故障衡量统计指标

$\left(\begin{array}{c}{\tilde{T}}_{f,m}={\tilde{X}}_{f,m}{\tilde{R}}_{f,m}\\ {{\tilde{D}}_{f,m}}^2={({\tilde{t}}_{f,m}-{\overline{\stackrel{~}{t}}}_{f,m})}^T{{\tilde{\Sigma }}_{f,m}}^{-1}({\tilde{t}}_{f,m}-{\overline{\stackrel{~}{t}}}_{f,m})\end{array}\right)---\left(15\right)$

其中，表示判别成分矩阵，是的行向量的转置，是的均值 向量，是的协方差矩阵；

(15.3)建立每类故障变量衡量指标的控制限；由于过程数据服从多维高斯 分布，且样本数量足够大，故变量衡量指标服从加权χ²分布，则可依据带权重的 χ²分布建立D²统计量的控制限

(16)将步骤(14)中普通变量构成的矩阵和分 别进行减均值除标准差的预处理后，构成大样本矩阵对进行PCA 分解建模，得到负载矩阵表示所监测的波动方向；

(17)计算普通变量矩阵的故障衡量指标

其中，其中，表示判别成分矩阵，是的行向量的转置，是的均 值向量，是的协方差矩阵；

(18)建立每类普通变量衡量指标的控制限；由于过程数据服从多维高斯分 布，且样本数量足够大，故变量衡量指标服从加权χ²分布，则可依据带权重的χ² 分布建立D²统计量的控制限

(19)在线故障诊断；当发生故障后，需要检测发生故障的变量以及故障的 类型，该步骤包括以下子步骤：

(19.1)按照步骤(1)获取新数据x_new(J×1)，x_new(J×1)可以被分为两部分， 故障变量部分和普通变量部分即

(19.2)调用故障诊断模型衡量故障变量：

$\left(\begin{array}{c}{\tilde{t}}_{\mathrm{new},m}={{\tilde{x}}_{\mathrm{new},m}}^T{\tilde{R}}_{f,m}\\ {{\tilde{D}}_{\mathrm{new},m}}^2={({\tilde{t}}_{\mathrm{new},m}-{\overline{\stackrel{~}{t}}}_{f,m})}^T{{\tilde{\Sigma }}_{f,m}}^{-1}({\tilde{t}}_{\mathrm{new},m}-{\overline{\stackrel{~}{t}}}_{f,m})\end{array}\right)---\left(17\right)$

其中，是新数据中故障变量部分，是的判别成分，和均 由步骤(15.2)求得，为故障变量衡量指标；

(19.3)调用故障诊断模型衡量普通变量：

其中，是新数据中普通变量部分，是的判别成分，和均由 步骤(17)求得，是普通变量衡量指标；

(19.4)依次比较和其控制限和其控制限若均没 有超限报警，则表示所调用的故障诊断模型能很好地诊断出故障，则判断该故障 属于调用的故障诊断模型所表示的类别；若有一个以上指标超限，则继续调用该 故障诊断模型对接下来的新故障数据进行诊断，按步骤(19.1)～(19.3)计算和若有连续新故障数据的和有一个以上指标超限，则表示故 障不属于该故障诊断模型所表示的故障类别；考虑到误报率和漏报率，需要对连 续多个新故障数据进行统计分析来判断故障的类别，记录新故障数据在调用每类 故障诊断模型时两个D²指标的未超限个数，则新故障数据属于D²指标的未超限 个数最大的故障类别。

本发明的有益效果是：该方法克服了传统方法存在的类内散布矩阵奇异性问 题、判别成分的个数局限问题、判别成分线性相关问题，使得方法可以应用于高 维度高耦合数据的复杂化工过程。该方法将传统的判别分析方法与变量选择方法 有效结合，区分了故障变量和普通变量，建立了更准确的诊断模型，能有效地区 分生产过程中故障类别，提高了在线故障诊断的性能，有助于工程师准确修复故 障，确保生产的安全可靠运行以及产品的高质量追求。

附图说明

图1是本发明基于嵌套迭代费舍尔判别分析的故障变量隔离方法的流程图；

图2是本发明方法调用故障1的诊断模型针对故障2的在线故障诊断结果图， (a)为普通变量，(b)为故障变量；

图3是本发明方法调用故障2的诊断模型针对故障2的在线故障诊断结果图， (a)为普通变量，(b)为故障变量；

图4是本发明方法调用故障9的诊断模型针对故障2的在线故障诊断结果图 (a)为普通变量，(b)为故障变量；

图5是本发明方法调用故障15的诊断模型针对故障2的在线故障诊断结果 图，(a)为普通变量，(b)为故障变量。

具体实施方式

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

以田纳西-伊斯曼过程为例，田纳西-伊斯曼过程是一个典型的复杂化工生产 过程，该过程由五个主要运行机构构成，分别是反应器、产品冷凝器、蒸汽液体 分离器、循环压缩机以及产品剥离器。在整个过程中可以采集到两部分过程变量： 41个测量变量以及11个操作变量。变量见表1和表2。

表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(产品分析)    

表2田纳西-伊斯曼过程操作变量表

序号 变量名称 序号 变量名称 序号 变量名称 1 A组分进料 2 D组分进料 3 E组分进料 4 A和C组分进料 5 循环阀门 6 放空阀门 7 分离器液相流量 8 汽提塔液相流量 9 汽提塔蒸汽阀门 10 反应器冷却水流量 11 冷凝器冷却水流量    

如图1所示，本发明一种基于嵌套迭代费舍尔判别分析的故障变量隔离方法， 包括以下步骤：

(1)获取过程分析数据：设一个化工生产过程具有J个测量变量和操作变量， 则每一次采样可得到一个1×J的向量，采样K次后得到的数据表述为一个二维矩 阵X(K×J)，所述测量变量为运行过程中可被测量的状态参数，包括流量、温度、 速率；所述操作变量包括给料量、阀门开度；分别获取正常数据二维矩阵X_n(K×J) 和故障数据二维矩阵X_f,m(K×J)，其中，下标n表示正常数据，下标f表示故障数 据，m表示故障的类别；将正常数据和故障数据统一标示为X_i(K×J)，其中下标 i表示数据的类别；

表3田纳西-伊斯曼过程故障表

序号 故障变量 发生类型 序号 变量名称 发生类型 1 APC进料变化 跃变 9 进料2温度变化 随机 2 成分B变化 跃变 10 进料4温度变化 随机 3 进料2温度变化 跃变 11 反应器冷却水温度变化 随机 4 反应器冷却水速度变化 跃变 12 冷凝器冷却水温度变化 随机 5 冷凝器冷却水速度变化 跃变 13 反应动力学特性变化 缓慢漂移 6 进料1损失 跃变 14 反应器冷却水阀门 粘滞 7 成分C进料压力下降 跃变 15 冷凝器冷却水阀门 粘滞 8 进料4中ABC组分变化 随机      

(2)选取正常数据样本和一类故障数据样本作为总样本其中， 由X_i(i＝1,2)从上到下排列组成；

(3)数据准备：分别计算总样本均值向量每类样本均值向量总类内 散布矩阵S_w和类间的散布矩阵S_b，计算公式如下：

$\left(\begin{array}{c}{S}_i=\underset{x_i\in X_i}{\Sigma }(x_i-{\bar{x}}_i){(x_i-{\bar{x}}_i)}^T\\ S_w=\underset{i=1}{\overset{2}{\Sigma }}{S}_i\\ S_b=\underset{x_i\in X_i}{\Sigma }{K}_i({\bar{x}}_i-\bar{x}){({\bar{x}}_i-\bar{x})}^T\end{array}\right)---\left(1\right)$

其中，S_i是每个类的散布矩阵；

(4)提取初始判别成分，该步骤由以下子步骤来实现：

(4.1)最大化类间离散度：求取使类间离散度最大的权重向量w，即相当于 求取类间散布矩阵S_b的最大特征值所对应的特征向量w，所述类间离散度为 w^TS_bw，获取w后，按公式(2)求取相应的总样本初始判别成分t；

$t=\bar{X}w---\left(2\right)$

其中，是减均值中心化后的总样本，那么对于每一类样本，其所对应的类 判别成分为可知，t由t_i从上到下依次排列构成；

(4.2)数据压缩：对减均值中心化后的总样本根据下式进行数据压缩：

$\left(\begin{array}{c}{p}^T={\left(t^{T}t\right)}^{-1}{t}^T\bar{X}\\ \bar{E}=\bar{X}-t{p}^T\end{array}\right)---\left(3\right)$

其中：p表示总样本的负载向量，表示总样本中与t无关的残差；

同理，对于每类样本都可以通过公式(4)得到与t_i无关的残差且由从上到下排列组成：

${\bar{E}}_i={\bar{X}}_i-t_i{p}^T={\bar{X}}_i-{\bar{X}}_{i}w{p}^T---\left(4\right)$

最后，用上述数据压缩关系wp^T更新每一个类的信息，以保证判别成分的正 交性：

E_i＝X_i-X_iwp^T                 (5)

(4.3)迭代更新过程数据

(4.3.1)用步骤(4.2)中获得的E_i代替步骤(3)中的X_i，按步骤(3)重 新计算总样本均值每类样本均值向量总类内散布矩阵S_w和类间的散布矩 阵S_b，按步骤(4.1)、(4.2)再次提取初始判别成分；

(4.3.2)重复步骤(4.3.1)直到所提取的初始判别成分的个数等于S_w的阶数 N；那么，同时可以得到由权重向量w组成的权重矩阵W(J×N)和相应的负载向 量p组成的负载矩阵P(J×N)、总样本初始判别成分t组成的总样本的初始判别成 分矩阵其中，T由T_i按从上至下排列构成，T_i是每个类的判别成分 矩阵；最后，求取初始判别成分的系数矩阵R＝W(P^TW)^-1，且T和T_i可直接由系 数矩阵根据公式(6)求出：

$\left(\begin{array}{c}T=\bar{X}R=\bar{X}W{\left(P^{T}W\right)}^{-1}\\ T_i={\bar{X}}_{i}R={\bar{X}}_{i}W{\left(P^{T}W\right)}^{-1}\end{array}\right)---\left(6\right)$

(5)提取最终判别成分，该步骤通过以下子步骤来实现：

(5.1)过程数据预处理：使用X_iR代替每类初始数据集合X_i，按步骤(3) 重新计算每类样本均值总样本均值总类内散布矩阵S_w^*以及类间散布矩阵 S_b^*；

(5.2)确定最终判别成分：最终判别成分通过以下步骤来确定：

(5.2.1)求取最优判别成分方向向量w^*，使得类间散布矩阵与类内散布矩阵 的比值J(θ)最大；其中，w^*则可通过公式(7)求 取矩阵S_w^*-1S_b^*最大特征值所对应的特征向量得到：

S_w^*-1S_b^*w^*＝λw^*                   (7)

(5.2.2)求取每类的最终判别成分向量t_i^*：

$\left(\begin{array}{c}{t_i}^{*}=X_{i}R{w}^{*}=X_i\theta \\ \theta =R{w}^{*}\end{array}\right)---\left(8\right)$

(5.2.3)将t_i^*从上至下依次排列构成总样本的最终判别成分向量

(5.3)压缩过程数据：为了保证每类样本的判别成分之间是正交的，进行如 下处理：

$\left(\begin{array}{c}{p_i}^{*T}={\left({t_i}^{*T}{t_i}^{*}\right)}^{-1}{t_i}^{*T}{X}_i\\ {E_i}^{*}=X_i-{t_i}^{*}{p_i}^{*T}\end{array}\right)---\left(9\right)$

其中，是每类的负载向量，E_i^*是与t_i^*无关的残差；

(6)迭代更新过程数据，该步骤包括以下子步骤：

(6.1)用步骤(5.3)中E_i^*代替步骤(3)中的X_i，按步骤(3)重新计算每 类样本均值、总样本均值，总类内散布矩阵S_w^*以及类间散布矩阵S_b^*，按步骤(4) 和步骤(5)再次提取最终判别成分向量t_i^*；

(6.2)重复步骤(6.1)直至获得足够的最终判别成分t_i^*并构成最终的判别 成分矩阵T_i^*，T_i^*所保留的最终判别成分个数为R，所述R通过交叉检验的方法确 定；相应的，同时可以获得权重矩阵Θ(J×R)和负载矩阵P_i^*(J×R)；其中，Θ(J×R) 和P_i^*(J×R)分别由θ(J×1)和p_i^*(J×1)构成；

(6.3)求取最终系数矩阵R_i^*(J×R)：

R_i^*＝Θ(P_i^*TΘ)^-1                (10)

至此，步骤(2)中选取的该类故障最终系数矩阵被求取出来；

(7)选取正常数据和另一类故障数据作为总样本，重复步骤(4)-(6)， 获得该类故障样本的最终系数矩阵；

(8)重复步骤(7)直到M类故障的最终系数矩阵都被求 取出来；

(9)故障数据衡量；该步骤由以下子步骤完成：

(9.1)求取每类故障X_f,m及相应正常数据X_n的类内成分矩阵T_f,m^*和T_n,m^*：

$\left(\begin{array}{c}{T_{f,m}}^{*}=X_{f,m}{R_{f,m}}^{*}\\ {T_{n,m}}^{*}=X_n{R_{f,m}}^{*}\end{array}\right)---\left(11\right)$

(9.2)求取每类故障X_f,m及相应正常数据X_n的D²指标：

$\left(\begin{array}{c}{t}_{n,m,k}={x_{n,k}}^T{R_{f,m}}^{*}\\ {D_{n,m,k}}^2={(t_{n,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}(t_{n,m,k}-{\bar{t}}_{n,m})\\ t_{f,m,k}={x_{f,k}}^T{R_{f,m}}^{*}\\ {D_{f,m,k}}^2={(t_{f,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}(t_{f,m,k}-{\bar{t}}_{n,m})\end{array}\right)---\left(12\right)$

其中，x_n,k是X_n中一个样本，t_n,m,k是x_n,k的判别成分，是公式(11)中T_n,m^*的均值向量，Σ_n则表示T_n,m^*的协方差矩阵，如果方差过小而趋近于0，那么它们 将被置为1；同理，x_f,m是X_f,m中一个样本，t_n,m,k是x_f,m的判别成分；

(9.2)建立正常数据基于D²指标的控制限：由于过程数据服从多变量正态 分布，那么可知D²统计量服从带权重的χ²分布，则可依据带权重的χ²分布建立 D²统计量的控制限

(10)比较每类故障数据的D²指标与相应的正常数据控制限若没有 超限报警则说明该类故障与正常数据相似，无需进行以下步骤；若有超限报警， 则执行以下步骤辨识和分离故障变量；

(11)变量贡献度衡量，该步骤包括以下子步骤；

(11.1)按公式(13)计算每类故障及相应正常数据的变量贡献度指标

$\left(\begin{array}{c}{t}_{·,m,k}={x_{·,k}}^T{R_{f,m}}^{*}\\ C_{D^2,·,k,j}={(t_{·,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}{r}_{m,j}(x_{·,k,j}-{\bar{x}}_{n,j})\end{array}\right)---\left(13\right)$

其中，下标·表示正常数据或者故障数据，x_·,k表示相应正常或者故障数据中 的一个样本，t_·,m,k表示相应样本x_·,k的判别成分，x_·,k,j是相应数据中的第k个样本 的第j个变量，是正常数据中第j个变量的均值，r_m,j是第j个变量的权重系数；

(11.2)因为正常数据服从一定的分布规律，根据分布规律确定正常数据贡 献度指标的控制限；

(11.3)按公式(14)计算故障样本相对于正常样本的变量贡献度的比值

${\mathrm{RC}}_{D^2,f,k,j}=\frac{C_{D^2,f,k,j}}{\mathrm{ctr}\left(C_{D^2,n,j}\right)}=\frac{{(t_{·,m,k}-{\bar{t}}_{n,m})}^T{{\Sigma }_n}^{-1}{r}_{m,j}(x_{·,k,j}-{\bar{x}}_{n,j})}{\mathrm{ctr}\left(C_{D^2,n,j}\right)}---\left(14\right)$

其中，表示正常数据的控制限；

(11.4)对于每个变量，将通过公式(14)所计算得到的M类故障对于正常 样本的变量贡献度的比值构成向量并计算的均值

(12)故障变量选择；将步骤(11.4)所得到的排序，选择最 大所对应的变量作为故障变量，将该变量从故障数据和正常数据 中移除，得到新的故障数据和正常数据

(13)模型更新；对更新后的故障数据和正常数据执行步骤(3)～ (8)，求取新的最终系数矩阵，即代表故障方向；

(14)基于更新后的最终系数矩阵，重复执行步骤(9)～(12)直到所有故 障变量都被选取出来；对于每类故障数据，故障变量构成矩阵其 中，J_f,m表示故障变量的个数，N_f表示故障样本的个数；则每类故障被分为两部 分数据：对故障有重大影响的以及对故障没有影响的其中，J_n,m＝J-J_f,m；同样，对于正常数据也被分为两部分和 其中，N_n表示正常数据的个数，J_n,m表示普通变量的个数；那么， 对于每一对正常数据X_n和故障数据X_f,m，可以得到由故障变量构成的矩阵 和以及由剩余的普通变量构成的矩阵和

(15)建立故障诊断模型，该步骤包括以下子步骤：

(15.1)用和分别代替步骤(2)中获取的正常样 本数据和故障数据，重复步骤(3)至(8)提取出新的最终系数矩阵即表 示有重要区分度的故障方向；

(15.2)计算每类故障子数据的判别成分和故障衡量统计指标

$\left(\begin{array}{c}{\tilde{T}}_{f,m}={\tilde{X}}_{f,m}{\tilde{R}}_{f,m}\\ {{\tilde{D}}_{f,m}}^2={({\tilde{t}}_{f,m}-{\overline{\stackrel{~}{t}}}_{f,m})}^T{{\tilde{\Sigma }}_{f,m}}^{-1}({\tilde{t}}_{f,m}-{\overline{\stackrel{~}{t}}}_{f,m})\end{array}\right)---\left(15\right)$

其中，表示判别成分矩阵，是的行向量的转置，是的均值 向量，是的协方差矩阵；

(15.3)建立每类故障变量衡量指标的控制限；由于过程数据服从多维高斯 分布，且样本数量足够大，故变量衡量指标服从加权χ²分布，则可依据带权重的 χ²分布建立D²统计量的控制限

(16)将步骤(14)中普通变量构成的矩阵和分 别进行减均值除标准差的预处理后，构成大样本矩阵对进行PCA 分解建模，得到负载矩阵表示所监测的波动方向；

(17)计算普通变量矩阵的故障衡量指标

其中，其中，表示判别成分矩阵，是的行向量的转置，是的均 值向量，是的协方差矩阵；

(18)建立每类普通变量衡量指标的控制限；由于过程数据服从多维高斯分 布，且样本数量足够大，故变量衡量指标服从加权χ²分布，则可依据带权重的χ² 分布建立D²统计量的控制限

(19)在线故障诊断；当发生故障后，需要检测发生故障的变量以及故障的 类型，该步骤包括以下子步骤：

(19.1)按照步骤(1)获取新数据x_new(J×1)，x_new(J×1)可以被分为两部分， 故障变量部分和普通变量部分即

(19.2)调用故障诊断模型衡量故障变量：

$\left(\begin{array}{c}{\tilde{t}}_{\mathrm{new},m}={{\tilde{x}}_{\mathrm{new},m}}^T{\tilde{R}}_{f,m}\\ {{\tilde{D}}_{\mathrm{new},m}}^2={({\tilde{t}}_{\mathrm{new},m}-{\overline{\stackrel{~}{t}}}_{f,m})}^T{{\tilde{\Sigma }}_{f,m}}^{-1}({\tilde{t}}_{\mathrm{new},m}-{\overline{\stackrel{~}{t}}}_{f,m})\end{array}\right)---\left(17\right)$

其中，是新数据中故障变量部分，是的判别成分，和均 由步骤(15.2)求得，为故障变量衡量指标；

(19.3)调用故障诊断模型衡量普通变量：

其中，是新数据中普通变量部分，是的判别成分，和均由 步骤(17)求得，是普通变量衡量指标；

(19.4)依次比较和其控制限和其控制限若均没 有超限报警，则表示所调用的故障诊断模型能很好地诊断出故障，则判断该故障 属于调用的故障诊断模型所表示的类别；若有一个以上指标超限，则继续调用该 故障诊断模型对接下来的新故障数据进行诊断，按步骤(19.1)～(19.3)计算和若有连续新故障数据的和有一个以上指标超限，则表示故 障不属于该故障诊断模型所表示的故障类别；考虑到误报率和漏报率，需要对连 续多个新故障数据进行统计分析来判断故障的类别，记录新故障数据在调用每类 故障诊断模型时两个D²指标的未超限个数，则新故障数据属于D²指标的未超限 个数最大的故障类别。

本发明通过基于迭代嵌套费舍尔判别分析的故障变量隔离方法根据历史故 障数据建立故障诊断模型，工程师可以实时获得新过程采样数据的在线故障诊断 结果，判断故障的类别从而进行相应的修复措施。当基于历史故障数据建立的诊 断模型能很好容纳新故障样本所求得的D²指标，即调用该类故障的诊断模型所 得的新样本D²指标产生的报警信号显著少于调用其他类故障诊断模型所产生的 报警信号，则说明新样本属于该类故障，那么当前故障被正确诊断出来。考虑到 误报率和漏报率，统计了新样本D²在每类故障诊断模型下的未超限信号的个数， 新样本则属于未超限个数最多的那类故障。为了描述故障诊断的性能，我们定义 了一个指标—正确诊断指标，若故障被正确诊断，则该指标为1；否则，该指标 值为0。表4总结了本发明所提出的方法与传统方法针对15种故障类型的正确诊 断指标的比较结果，并且计算均值(Mean)以及标准差(STD)用于综合评估故 障诊断的性能。

表4本方法与传统方法故障诊断性能对比(衡量指标：正确诊断指标)

从表中可以看出，基于本发明方法的故障诊断方法提高了实际在线故障诊断 的可靠性和可信度，对于15种故障均能正确诊断出来，其性能优于传统方法。 此外，图2至图5展示了本发明所提出方法针对于故障2的在线故障诊断结果。 是实线表示每类故障基于D²指标的控制限，点线表示D²指标。可以看出，而本 发明则能正确诊断出故障2，说明了本发明方法的有效性。总体来说，本发明的 方法有优越的在线故障诊断性能，有助于工业工程师对故障进行准确修复，保证 实际生产过程的安全可靠运行。

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