CN105182955B - 一种多变量工业过程故障识别方法 - Google Patents

Abstract

本发明涉及一种多变量工业过程故障识别方法，含有以下步骤：(一)收集历史数据库的正常操作工况数据集X和K类已知的故障模式数据集，计算正常操作工况数据集的均值mean(X)和标准差std(X)，对已知的故障模式数据集进行标准化处理获得新故障模式数据集。(二)在各个故障模式数据集下构造数据窗，计算六种统计量变量。(三)检测过程故障，收集实时故障数据S，进行标准化处理。(四)在步骤(三)的基础上执行统计量主元相异度分析，计算待识别故障数据集和已知故障模式数据集之间的故障识别指数FRI。(五)对故障识别指数FRI进行排序，获得故障识别结果。本发明基于统计量主元相异度分析，在相异度分析中，提取主元信息，摒弃次要数据信息，抑制噪声的影响，能够充分挖掘数据高阶统计信息。

Description

一种多变量工业过程故障识别方法

技术领域

本发明属于工业过程故障识别技术领域，具体地说，涉及一种基于统计量主元相异度分析的多变量工业过程故障识别方法。

背景技术

随着高集成、大规模化的现代工业系统日益形成，故障诊断技术成为确保现代工业系统安全稳定运行的一项关键技术，并且由于工业过程数据库中存储了丰富的运行数据，基于数据驱动的故障诊断技术已经引起了工业过程工程师和学术研究人员的广泛关注。研究人员提出了一系列基于数据驱动的故障诊断方法，包括主元分析(PCA)、独立元分析(ICA)、偏最小二乘(PLS)。目前的故障诊断方法研究多数集中在故障检测问题(即如何快速有效的发现过程故障)，然而故障源诊断问题(即检测到故障后识别故障的类型和原因)的研究相对较少，后者是一种更具有挑战性的研究问题。

当工业过程数据库中存在一些已知的故障模式数据时，故障模式识别方法就成为故障源诊断的一种有效手段。对于复杂的故障模式识别问题，可以利用的方法有费舍尔判别子空间方法、支持向量机方法进行故障辨识。Kano等人首先提出相异度分析(参考文献：ManabuKano,ShinjiHasebe,IoriHashimoto,HiromuOhno,Statisticalprocessmonitoringbasedondissimilarityofprocessdata,AIChEJournal,48(2002)1231-1240)。近几年来，相异度分析(Dissimilarityanalysis)作为一种故障诊断和识别的有效技术，引起了广泛的研究关注。该方法通过比较未知故障数据和已知故障模式数据之间的相异度指数来判断故障模式。虽然该方法取得初步的成功应用，但是其缺陷在于：(1)在相异度指数计算过程中没有考虑噪声抑制问题，但是噪声信息会恶化识别效果；(2)相异度分析直接针对原始测量变量，无法充分利用数据的高阶统计信息。

发明内容

本发明针对现有技术存在的上述不足，提供一种简单、精确有效的多变量工业过程故障识别方法，该方法基于统计量主元相异度分析，在相异度分析中，提取主元信息，摒弃次要数据信息，抑制噪声的影响，能够充分挖掘数据高阶统计信息。

本发明的技术方案是：一种多变量工业过程故障识别方法，含有以下步骤：

(一)收集历史数据库的正常操作工况数据集X和K类已知的故障模式数据集{H_o1,H_o2,...,H_oK}，计算正常操作工况数据集的均值mean(X)和标准差std(X)，对已知的故障模式数据集{H_o1,H_o2,...,H_oK}进行标准化处理获得新故障模式数据集{H₁,H₂,...,H_K}。

(二)在各个故障模式数据集下构造数据窗，计算低阶统计量的变量均值ε_i(t)、低阶统计量的方差v_i(t)、高阶统计量的3阶中心距高阶统计量的4阶中心距高阶统计量的偏斜度γ_i(t)和高阶统计量的峰度κ_i(t)六种统计量变量。

(三)检测过程故障，收集实时故障数据S，进行标准化处理。

(四)在步骤(三)的基础上执行统计量主元相异度分析，计算待识别故障数据集和已知故障模式数据集之间的故障识别指数FRI。

(五)对故障识别指数FRI进行排序，获得故障识别结果。

进一步的，所述步骤(一)中，通过公式(1)对故障模式数据进行标准化处理，公式(1)的表达式如下：

H_i＝(H_oi-mean(X))/std(X),(i＝1,2,...,K)(1)

经上述公式(1)标准化处理后即可获得新故障模式数据集{H₁,H₂,...,H_K}。

进一步的，所述步骤(二)中，构造数据窗的步骤为：对任意t时刻的数据窗记为：

式中，w为数据窗宽度。

进一步的，所述六种统计变量分别由计算公式(3)、(4)、(5)、(6)、(7)、(8)计算获得，公式(3)、(4)、(5)、(6)、(7)、(8)的表达式如下所示：

${\mathrm{\ϵ}}_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{x}_i(t-h)---\left(3\right)$

$v_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2---\left(4\right)$

$c_i^{\left(3\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3---\left(5\right)$

$c_i^{\left(4\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4---\left(6\right)$

${\mathrm{\γ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i\left(t-h\right)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i\left(t-h\right)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^{3/2}}---\left(7\right)$

${\mathrm{\κ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i\left(k-h\right)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i\left(k-h\right)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^2}-3---\left(8\right)$

根据上面所述的六个公式(3)至(8)获得数据窗的6类变量，将其堆积为一个行向量[ε(t)v(t)c⁽³⁾(t)c⁽⁴⁾(t)γ(t)κ(t)]，其中，ε(t)＝[ε₁(t)…ε_m(t)]，v(t)＝[v₁(t)…v_m(t)]， $c^{\left(3\right)}\left(t\right)=\[c_1^{\left(3\right)}\left(t\right)...c_m^{\left(3\right)}\left(t\right)\],c^{\left(4\right)}\left(t\right)=\[c_1^{\left(4\right)}\left(t\right)...c_m^{\left(4\right)}\left(t\right)\],$ γ(t)＝[γ₁(t)…γ_m(t)]，κ(t)＝[κ₁(t)…κ_m(t)]；

将各个故障模式数据集分割为一系列如公式(9)所示的数据窗H_it(w≤t≤n)，公式(9)表达式如下：

$H_t=\left[\begin{array}{c}{H}_{iw}\\ .\\ .\\ .\\ H_{it}\\ .\\ .\\ .\\ H_{in}\end{array}\right]---\left(9\right)$

针对各故障模式数据集的每一个数据窗H_it(w≤t≤n)计算其统计量，进而获得各个故障模式数据集的统计量矩阵，记为：

$X_i=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right],(i=1,2,\mathrm{...},K)---\left(10\right)$

进一步的，所述步骤(四)中，执行统计量主元相异度分析的步骤为：首先构造实时数据S的窗宽为w的数据窗，根据公式(3)至(8)计算待识别故障数据集的统计量，记为：

$X_S=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right]---\left(11\right)$

对于已知故障模式数据集的统计量矩阵(i＝1,2,...,K)和待识别故障数据统计量分别包含m个变量的n_i(i＝1,...,K)和n_s个样本，其协方差阵计算公式为：

$R_i=\frac{1}{n_i-1}{X}_i^T{X}_i,R_S=\frac{1}{n_s-1}{X}_S^T{X}_S---\left(12\right)$

对第i个故障模式，构造总体统计量数据集的协方差矩阵，记为：

$R=\frac{1}{n_i+n_s-1}{X}^{T}X---\left(13\right)$

对其开展特征值分解，获得：

RP＝PΛ(14)

式中，P是特征向量矩阵，Λ(λ₁≥λ₂≥...≥λ_m)是特征值矩阵；

保留L个特征值和对应的特征向量p₁,p₂,...,p_L，对矩阵X_i和X_S进行主元线性变换，得到Y_i和Y_s，记为：

$Y_i=\sqrt{\frac{n_i-1}{n_i+n_s-1}}{X}_i{P}_L{\mathrm{\Λ}}_L^{-1/2}---\left(15\right)$

$Y_S=\sqrt{\frac{n_s-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(16\right)$

式中，P_k＝[p₁,p₂,...,p_k]，Λ_k为{λ₁,λ₂,...,λ_k}构成的对角矩阵；

Y_i和Y_s的协方差矩阵记为：

$C_i=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_i^T{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(17\right)$

$C_S=\frac{1}{n_1+n_s-1}{A_L}^{-1/2}{P}_L^T{X}_S^T{X}_S{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(18\right)$

两个数据集的协方差矩阵满足下式：

$C_i+C_S={{\mathrm{\Λ}}_k}^{-1/2}{P}_k^T{\mathrm{RP}}_k{{\mathrm{\Λ}}_k}^{-1/2}={{\mathrm{\Λ}}_k}^{-1/2}{P}_k^T{P}_k{\mathrm{\Λ}}_k{{\mathrm{\Λ}}_k}^{-1/2}=I---\left(19\right)$

对协方差矩阵应用特征值分解，得到

$C_i{w}_j^{\left(i\right)}={\mathrm{\λ}}_j^{\left(i\right)}{w}_j^{\left(i\right)},j=1,2,\mathrm{...},k---\left(20\right)$

$C_S{w}_j^{\left(s\right)}={\mathrm{\λ}}_j^{\left(s\right)}{w}_j^{\left(s\right)}---\left(21\right)$

由于和C_S的特征向量与C_i的特征向量相同，仅仅是次序相反，对相同大小的数据集X_i和X_S，如果具有相似的数据分布和关联关系，则特征值λ_j和1-λ_j均接近0.5；相反的，如果数据集描述了不同的数据分布和数据关系，则特征值分别趋近1和0；

计算各故障模式的统计量变量与待识别故障数据统计量的故障识别指数FRI的步骤为：相异度指数D可以用来评价数据集之间的相似程度，定义为：

$D=\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(22\right)$

式中，λ_j是C_i或者C_s的特征值，如果两个数据集相似程度高，则特征值接近0.5，D接近0；

为了方便故障识别，定义故障识别指数FRI为S_DA，其计算公式为：

$S_{DA}=1-D=1-\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(23\right)$

当两个数据集相似程度高时，S_DA接近1。

进一步的，所述步骤(三)中，过程检测到故障后，收集实时故障数据集S_o，使用公式S＝(S_o-mean(X))/std(X)进行标准化，上式中，X为正常操作工况数据集。

进一步的，所述步骤(五)中，对故障识别指数进行排序，最大的故障识别指数为故障识别结果X_imax，其中，故障识别结果的序号为：

$i_{max}=\arg \underset{i}{max}\{S_{SPDA}\left(i\right),i=1,2,\mathrm{...},K\}---\left(24\right)$

真实故障模式的序号为i_real，如果i_real＝i_max，则故障识别结果是正确的。

进一步的，所述步骤(五)中，获得故障识别结果后，为了评价不同方法的故障识别结果，通过识别误差指数REI和识别对比度RCD两个故障识别性能指标进行不同方法的故障识别效果对比。

进一步的，所述识别误差指数REI的计算公式为：

$REI=\frac{1}{K}\underset{i=1}{\overset{K}{\Σ}}\left|{\mathrm{re}}_i\right|+\mathrm{\ρ}---\left(25\right)$

式中，re_i是识别误差，其计算公式为：

${\mathrm{re}}_i=\left\{\begin{array}{ccc}{S}_{SPDA}\left(i\right)-1,& for& i=i_{real}\\ S_{SPDA}\left(i\right)-0,& for& i\≠i_{real}\end{array}\right.---\left(26\right)$

式中，ρ是惩罚因子，错误的识别结果为1，正确的识别结果为0。

进一步的，所述识别对比度的计算公式为:

$RCD=\mathrm{\γ}\frac{S_{SPDA}\left(i_{max}\right)-S_{SPDA}\left(i_{submax}\right)}{S_{SPDA}\left(i_{max}\right)}---\left(27\right)$

式中，S_SPDA(i_max)是最大的FRI，S_SPDA(i_submax)是次大的FRI，γ是故障识别错误参数，对于错误的识别结果，γ是-1，否则是1。

本发明的有益效果是：(1)本发明基于统计量主元相异度分析，在传统相异度分析方法的基础上引入噪声抑制机制，在相异度分析中，提取主元信息，摒弃次要数据信息，抑制工业噪声对故障辨识结果的影响，防止噪声信息恶化识别结果。(2)本发明引入统计模式分析技术，首次计算原始测量变量的统计量变量，故障识别时构建原始测量变量的六个统计量变量，对统计量数据集进行主元相异度分析，能够充分挖掘数据高阶统计信息，改善故障识别结果。

附图说明

图1为本发明统计量主元相异度分析的故障识别流程图。

图2为本发明SPDA方框图。

图3为本发明故障D03T下的规范化FRI指数。

图4为本发明对CSTR系统各故障REI性能指标对比。

图5为本发明对CSTR系统各故障RCD性能指标对比。

图6为本发明TE化工过程的流程图。

具体实施方式

下面结合一个仿真实例并结合附图对本发明作出进一步说明。

实施例一：以CSTR系统为例，CSTR系统为一个化学反应器，物料A的溶液进入反应器，发生一阶不可逆化学方应，生成物料B，反应为放热反应，因此需要通过外面的夹套冷却剂取走反应热量，为保证过程正常运行，一般带有串级控制系统。

根据过程机理，建立CSTR系统的动态机理模型如下：

$\frac{{\mathrm{dc}}_A}{dt}=-k_0{e}^{-E/RT}+\frac{Q_F{c}_{AF}-Q_F{c}_A}{Ah}---\left(28\right)$

$\frac{dT}{dt}=\frac{{k_0}^{-E/RT}{c}_A(-\mathrm{\Δ}H)}{{\mathrm{\ρC}}_p}+\frac{Q_F{T}_F-Q_{F}T}{Ah}+\frac{{\mathrm{UA}}_C(T_C-T)}{{\mathrm{\ρC}}_{p}Ah}---\left(29\right)$

$\frac{{\mathrm{dT}}_C}{dt}=\frac{Q_C(T_{CF}-T_C)}{V_C}+\frac{{\mathrm{UA}}_C(T-T_C)}{{\mathrm{\ρ}}_C{C}_{pC}{V}_C}---\left(30\right)$

$\frac{dh}{dt}=\frac{Q_F-Q}{A}---\left(31\right)$

式中，A是反应器截面积，c_A是反应器内物料A的浓度，c_AF是物料A在进料中的浓度，C_p是反应物比热，C_pC是冷却剂比热，E是活化能，h是反应器液位，k₀是反应因子，Q_F进料流量，Q_C是冷却剂流量，R是气体常数，T是反应器内温度，T_C是冷却剂出口温度，T_CF是冷却剂入口温度，T_F是反应器进料温度，U是换热系数，A_C是总的换热面积，ΔH是反应热，ρ是反应物密度，是冷却剂密度。

根据动态机理模型，对CSTR系统进行仿真。在仿真过程中，采集反应器进料流量、反应器进料温度、进料中物料A浓度、反应器内温度、反应器液位、反应器出料流量、反应器出料中物料A浓度、冷却剂入口温度、冷却剂出口温度和冷却剂流量10个测量变量。

CSTR仿真过程包含了1个正常操作模态和8个故障操作模态，各包含500个数据样本。故障模态列表见表1。表1中包含两组故障参数，参数1用于产生历史故障模式数据D01H～D08H.参数2用于产生测试数据D01T～D08T.

表1

故障	描述	参数1	参数2
				D01	冷却剂温度缓慢改变	-0.1K/min	+0.1K/min
D02	进料浓度缓慢改变	+6.0×10-4(molL)min	+8.0×10-4(molL)min
				D03	进料温度缓慢改变	+0.05K/min	+0.07K/min
D04	冷却剂阀门粘滞	29.5％	31.5％
				D05	催化剂失活	+2K/min	+4K/min
D06	换热器结垢	-120(J/(min.K))/min	-100(J/(min.K))/min

D07	冷却剂出口温度测量偏差	+3K	-2K
				D08	反应器出口流量测量偏差	+3L/min	+5L/min

上述CSTR系统的多变量工业过程故障识别方法，含有以下步骤：

(一)收集CSTR系统历史数据库的正常操作工况数据集X和K类已知的故障模式数据集{H_o1,H_o2,...,H_oK}，计算正常操作工况数据集的均值mean(X)和标准差std(X)，对已知的故障模式数据集{H_o1,H_o2,...,H_oK}进行标准化处理获得新故障模式数据集{H₁,H₂,...,H_K}。

通过公式(1)对故障模式数据进行标准化处理，公式(1)的表达式如下：

H_i＝(H_oi-mean(X))/std(X),(i＝1,2,...,K)(1)

经上述公式(1)标准化处理后即可获得新故障模式数据集{H₁,H₂,...,H_K}。

(二)在各个故障模式数据集下构造数据窗，计算低阶统计量的变量均值ε_i(t)、低阶统计量的方差v_i(t)、高阶统计量的3阶中心距高阶统计量的4阶中心距高阶统计量的偏斜度γ_i(t)和高阶统计量的峰度κ_i(t)六种统计量变量。

构造数据窗的步骤为：对任意t时刻的数据窗记为：

式中，w为数据窗宽度。

所述六种统计变量分别由计算公式(3)、(4)、(5)、(6)、(7)、(8)计算获得，公式(3)、(4)、(5)、(6)、(7)、(8)的表达式如下所示：

${\mathrm{\ϵ}}_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{x}_i(t-h)---\left(3\right)$

$v_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2---\left(4\right)$

$c_i^{\left(3\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3---\left(5\right)$

$c_i^{\left(4\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4---\left(6\right)$

${\mathrm{\γ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i\left(t-h\right)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i\left(t-h\right)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^{3/2}}---\left(7\right)$

${\mathrm{\κ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(k-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(k-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^2}-3---\left(8\right)$

根据上面所述的六个公式(3)至(8)获得数据窗的6类变量，将其堆积为一个行向量[ε(t)v(t)c⁽³⁾(t)c⁽⁴⁾(t)γ(t)κ(t)]，其中，ε(t)＝[ε₁(t)…ε_m(t)]，v(t)＝[v₁(t)…v_m(t)]， $c^{\left(3\right)}\left(t\right)=\[c_1^{\left(3\right)}\left(t\right)...c_m^{\left(3\right)}\left(t\right)\],c^{\left(4\right)}\left(t\right)=\[c_1^{\left(4\right)}\left(t\right)...c_m^{\left(4\right)}\left(t\right)\],$ γ(t)＝[γ₁(t)…γ_m(t)]，κ(t)＝[κ₁(t)…κ_m(t)]；

将各个故障模式数据集分割为一系列如公式(9)所示的数据窗H_it(w≤t≤n)，公式(9)表达式如下：

$H_t=\left[\begin{array}{c}{H}_{iw}\\ .\\ .\\ .\\ H_{it}\\ .\\ .\\ .\\ H_{in}\end{array}\right]---\left(9\right)$

针对各故障模式数据集的每一个数据窗H_it(w≤t≤n)计算其统计量，进而获得各个故障模式数据集的统计量矩阵，记为：

$X_i=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right],(i=1,2,\mathrm{...},K)---\left(10\right)$

(三)检测过程故障，收集实时故障数据S，进行标准化处理。

过程检测到故障后，收集实时故障数据集S_o，使用公式S＝(S_o-mean(X))/std(X)进行标准化，上式中，X为正常操作工况数据集。

(四)在步骤(三)的基础上执行统计量主元相异度分析，计算待识别故障数据集和已知故障模式数据集之间的故障识别指数FRI。

执行统计量主元相异度分析的步骤为：首先构造实时数据S的窗宽为w的数据窗，根据公式(3)至(8)计算待识别故障数据集的统计量，记为：

$X_S=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right]---\left(11\right)$

对于已知故障模式数据集的统计量矩阵(i＝1,2,...,K)和待识别故障数据统计量分别包含m个变量的n_i(i＝1,...,K)和n_s个样本，其协方差阵计算公式为：

$R_i=\frac{1}{n_i-1}{X}_i^T{X}_i,R_S=\frac{1}{n_s-1}{X}_S^T{X}_S---\left(12\right)$

对第i个故障模式，构造总体统计量数据集的协方差矩阵，记为：

$R=\frac{1}{n_i+n_s-1}{X}^{T}X---\left(13\right)$

对其开展特征值分解，获得：

RP＝PΛ(14)

式中，P是特征向量矩阵，Λ(λ₁≥λ₂≥...≥λ_m)是特征值矩阵；

大的特征值意味着重要的数据变化，而小的特征值意味着微弱的数据波动。前者是故障识别的主要因素，而后者往往是数据噪声导致的偶然变化。因此，保留L个特征值和对应的特征向量p₁,p₂,...,p_L，这前L个特征向量作为主元载荷向量，尽可能保留主要的数据变化，压制数据噪声。

对矩阵X_i和X_S进行主元线性变换，得到Y_i和Y_s，记为：

$Y_i=\sqrt{\frac{n_i-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(15\right)$

$Y_S=\sqrt{\frac{n_s-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(16\right)$

式中，P_k＝[p₁,p₂,...,p_k]，Λ_k为{λ₁,λ₂,...,λ_k}构成的对角矩阵；

Y_i和Y_s的协方差矩阵记为：

$C_i=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_i^T{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(17\right)$

$C_S=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_S^T{X}_S{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(18\right)$

两个数据集的协方差矩阵满足下式：

$C_i+C_S={{\mathrm{\Λ}}_k}^{-1/2}{P}_k^T{\mathrm{RP}}_k{{\mathrm{\Λ}}_k}^{-1/2}={\mathrm{\Λ}}_k^{-1/2}{P}_k^T{P}_k{\mathrm{\Λ}}_k{{\mathrm{\Λ}}_k}^{-1/2}=I---\left(19\right)$

对协方差矩阵应用特征值分解，得到

$C_i{w}_j^{\left(i\right)}={\mathrm{\λ}}_j^{\left(i\right)}{w}_j^{\left(i\right)},j=1,2,...,k---\left(20\right)$

$C_S{w}_j^{\left(s\right)}={\mathrm{\λ}}_j^{\left(s\right)}{w}_j^{\left(s\right)}---\left(21\right)$

由于和C_S的特征向量与C_i的特征向量相同，仅仅是次序相反，对相同大小的数据集X_i和X_S，如果具有相似的数据分布和关联关系，则特征值λ_j和1-λ_j均接近0.5；相反的，如果数据集描述了不同的数据分布和数据关系，则特征值分别趋近1和0。

计算各故障模式的统计量变量与待识别故障数据统计量的故障识别指数FRI的步骤为：相异度指数D可以用来评价数据集之间的相似程度，定义为：

$D=\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(22\right)$

式中，λ_j是C_i或者C_s的特征值，如果两个数据集相似程度高，则特征值接近0.5，D接近0；

为了方便故障识别，定义故障识别指数FRI为S_DA，其计算公式为：

$S_{DA}=1-D=1-\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(23\right)$

当两个数据集相似程度高时，S_DA接近1。

(五)对故障识别指数FRI进行排序，获得故障识别结果。

对故障识别指数进行排序，最大的故障识别指数为故障识别结果X_imax，其中，故障识别结果的序号为：

$i_{max}=\arg \underset{i}{max}\{S_{SPDA}\left(i\right),i=1,2,\mathrm{...},K\}---\left(24\right)$

真实故障模式的序号为i_real，如果i_real＝i_max，则故障识别结果是正确。

最理想的识别结果是：第i_max个FRI接近1，其他FRI均接近0。实际识别过程中，有可能会出现多个FRI都具有比较大的数值，这样就增加故障识别的难度，因为，即使i_real＝i_max，多个较大的FRI也会降低了故障识别结果的可信度。

获得故障识别结果后，为了评价不同方法的故障识别结果，通过识别误差指数REI和识别对比度RCD两个故障识别性能指标进行不同方法的故障识别效果对比。

所述识别误差指数REI的计算公式为：

$REI=\frac{1}{K}\underset{i=1}{\overset{K}{\Σ}}\left|{\mathrm{re}}_i\right|+\mathrm{\ρ}---\left(25\right)$

式中，re_i是识别误差，其计算公式为：

${\mathrm{re}}_i=\left\{\begin{array}{ccc}{S}_{SPDA}\left(i\right)-1,& for& i=i_{real}\\ S_{SPDA}\left(i\right)-0,& for& i\≠i_{real}\end{array}\right.---\left(26\right)$

式中，ρ是惩罚因子，错误的识别结果为1，正确的识别结果为0。

所述识别对比度的计算公式为:

$RCD=\mathrm{\γ}\frac{S_{SPDA}\left(i_{max}\right)-S_{SPDA}\left(i_{submax}\right)}{S_{SPDA}\left(i_{max}\right)}---\left(27\right)$

式中，S_SPDA(i_max)是最大的FRI，S_SPDA(i_submax)是次大的FRI，γ是故障识别错误参数，对于错误的识别结果，γ是-1，否则是1。

很显然，大的RCD数值意味着清晰正确的识别结果，可信度高；小的RCD数值意味着识别结果不正确，可信度低。

本实施例中，采用相异度分析DA法和主元相异度分析PDA法两种方法与本发明基于统计量主元相异度分析SPDA的故障识别方法对CSTR系统进行故障模式识别。

统计量数据窗宽度选为30，综合对比上述三种方法，以故障D03为例说明故障识别效果。

故障测试数据集为D30T，为了比较不同方法，分别计算上述三种故障识别指数S_DA,S_PDA,S_SPDA，结果如表2所示。

表2

模式号	D01H	D02H	D03H	D04H	D05H	D06H	D07H	D08H
									SDA	0.8000	0.8005	0.9790	0.8607	0.8164	0.8020	0.8077	0.8184
SPDA	0.0471	0.0447	0.9747	0.7770	0.6447	0.0459	0.5400	0.7052
									SSPDA	0.0023	0.0015	0.9000	0.0433	0.0054	0.0023	0.0058	0.0050

为了便于比较，对数据进行规范化处理，将表2中的每行数据均除以其最大值，获得规范化故障识别指数S_DA,S_PDA,S_SPDA，结果如表3所示，同时绘制于图3中。

表3

模式号	D01H	D02H	D03H	D04H	D05H	D06H	D07H	D08H
									SDA	0.8172	0.8177	1.0000	0.8792	0.8339	0.8192	0.8251	0.8360
SPDA	0.0483	0.0459	1.0000	0.7972	0.6614	0.0471	0.5540	0.7235
									SSPDA	0.0025	0.0016	1.0000	0.0481	0.0060	0.0025	0.0065	0.0055

由图3可以看出，三种方法都可以正确识别故障数据集D03T。然而，DA方法结果中的第1、2、4、5、6、7、8故障识别指数都非常大，均在0.8以上，对比度不清晰。在应用PDA方法时，第4、5、7、8故障识别指数仍然比较大。应用SPDA方法时取得最好的故障识别效果，不但第3个故障识别指数最大，而且其他的识别结果都非常小，均小于0.1。从该故障案例的分析中可以看出，PDA方法可以改进DA方法，而SPDA方法能够进一步改进PDA方法。

为了有一个整体的衡量，计算测试数据集D03T识别过程中的性能指标REI和RCD，列于表4中。

表4

Methods	DA	PDA	SPDA
				REI	0.7285	0.3597	0.0091
RCD	0.1208	0.2028	0.9519

由表4可以看出，SPDA方法具有最小的REI结果，而且具有最大的RCD结果，这说明SPDA方法的故障识别误差最小，故障识别对比度最高，故障识别效果最佳。

计算全部8个测试故障数据集的REI和RCD指标，列于表5和表6中，其中，表5为对CSTR系统各故障REI性能指标的对比，表6为对CSTR系统各故障RCD性能指标的对比，为了有一个直观的表现可见图4和图5。

表5

表6

从上述表5、表6以及图4、图5中可以看出，PDA能够有效改善REI指标，SPDA则进一步减小REI指标；同时PDA的RCD指标大于DA方法，SPDA的RCD指标大于PDA方法。因此，PDA和SPDA具有更好的故障识别性能，以本发明SPDA的故障识别效果最佳。

实施例二：以一个著名的TE化工过程为例，TE化工过程已经广泛应用于故障诊断和过程控制方法研究中。TE化工过程由五个单元组成，共有52个过程变量，其包含11个操作变量，19个成分变量和22个连续测量变量，TE化工过程的流程图如图6所示。关于TE化工过程的其它细节详见参考文献(J.J.Downs,E.F.Vogel,Aplant-wideindustrialprocesscontrolproblem,Computers&ChemicalEngineering,17(1993)245-255)。

TE化工过程的仿真数据包可在网站http://web.mit.edu/braatzgroup/ index.html上下载。数据包中含有1个正常工况和21个故障工况数据IDV(1)～IDV(21)。正常工况数据有500个采样点52个变量。故障数据集分为两个类别，一类是故障模式数据集F01H～F21H，另一类是故障测试数据集F01T～F021T，各有480个采样点和52个变量值。

上述TE化工过程的多变量工业过程故障识别方法，含有以下步骤：

(一)收集TE化工过程历史数据库的正常操作工况数据集X和K类已知的故障模式数据集{H_o1,H_o2,...,H_oK}，计算正常操作工况数据集的均值mean(X)和标准差std(X)，对已知的故障模式数据集{H_o1,H_o2,...,H_oK}进行标准化处理获得新故障模式数据集{H₁,H₂,...,H_K}。

通过公式(1)对故障模式数据进行标准化处理，公式(1)的表达式如下：

H_i＝(H_oi-mean(X))/std(X),(i＝1,2,...,K)(1)

经上述公式(1)标准化处理后即可获得新故障模式数据集{H₁,H₂,...,H_K}。

(二)在各个故障模式数据集下构造数据窗，计算低阶统计量的变量均值ε_i(t)、低阶统计量的方差v_i(t)、高阶统计量的3阶中心距高阶统计量的4阶中心距高阶统计量的偏斜度γ_i(t)和高阶统计量的峰度κ_i(t)六种统计量变量。

构造数据窗的步骤为：对任意t时刻的数据窗记为：

式中，w为数据窗宽度。

所述六种统计变量分别由计算公式(3)、(4)、(5)、(6)、(7)、(8)计算获得，公式(3)、(4)、(5)、(6)、(7)、(8)的表达式如下所示：

${\mathrm{\ϵ}}_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{x}_i(t-h)---\left(3\right)$

$v_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2---\left(4\right)$

$c_i^{\left(3\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3---\left(5\right)$

$c_i^{\left(4\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4---\left(6\right)$

${\mathrm{\γ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^{3/2}}---\left(7\right)$

${\mathrm{\κ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(k-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(k-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^2}-3---\left(8\right)$

根据上面所述的六个公式(3)至(8)获得数据窗的6类变量，将其堆积为一个行向量[ε(t)v(t)c⁽³⁾(t)c⁽⁴⁾(t)γ(t)κ(t)]，其中，ε(t)＝[ε₁(t)…ε_m(t)]，v(t)＝[v₁(t)…v_m(t)]， $c^{\left(3\right)}\left(t\right)=\[c_1^{\left(3\right)}\left(t\right)...c_m^{\left(3\right)}\left(t\right)\],c^{\left(4\right)}\left(t\right)=\[c_1^{\left(4\right)}\left(t\right)...c_m^{\left(4\right)}\left(t\right)\],$ γ(t)＝[γ₁(t)…γ_m(t)]，

κ(t)＝[κ₁(t)…κ_m(t)]；

将各个故障模式数据集分割为一系列如公式(9)所示的数据窗H_it(w≤t≤n)，公式(9)表达式如下：

$H_t=\left[\begin{array}{c}{H}_{iw}\\ .\\ .\\ .\\ H_{it}\\ .\\ .\\ .\\ H_{in}\end{array}\right]---\left(9\right)$

针对各故障模式数据集的每一个数据窗H_it(w≤t≤n)计算其统计量，进而获得各个故障模式数据集的统计量矩阵，记为：

$X_i=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right],(i=1,2,\mathrm{...},K)---\left(10\right)$

(三)检测过程故障，收集实时故障数据S，进行标准化处理。

过程检测到故障后，收集实时故障数据集S_o，使用公式S＝(S_o-mean(X))/std(X)进行标准化，上式中，X为正常操作工况数据集。

(四)在步骤(三)的基础上执行统计量主元相异度分析，计算待识别故障数据集和已知故障模式数据集之间的故障识别指数FRI。

执行统计量主元相异度分析的步骤为：首先构造实时数据S的窗宽为w的数据窗，根据公式(3)至(8)计算待识别故障数据集的统计量，记为：

$X_S=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right]---\left(11\right)$

对于已知故障模式数据集的统计量矩阵(i＝1,2,...,K)和待识别故障数据统计量分别包含m个变量的n_i(i＝1,...,K)和n_s个样本，其协方差阵计算公式为：

$R_i=\frac{1}{n_i-1}{X}_i^T{X}_i,R_S=\frac{1}{n_s-1}{X}_S^T{X}_S---\left(12\right)$

对第i个故障模式，构造总体统计量数据集的协方差矩阵，记为：

$R=\frac{1}{n_i+n_s-1}{X}^{T}X---\left(13\right)$

对其开展特征值分解，获得：

RP＝PΛ(14)

式中，P是特征向量矩阵，Λ(λ₁≥λ₂≥...≥λ_m)是特征值矩阵。

大的特征值意味着重要的数据变化，而小的特征值意味着微弱的数据波动。前者是故障识别的主要因素，而后者往往是数据噪声导致的偶然变化。因此，保留L个特征值和对应的特征向量p₁,p₂,...,p_L，这前L个特征向量作为主元载荷向量，尽可能保留主要的数据变化，压制数据噪声。

对矩阵X_i和X_S进行主元线性变换，得到Y_i和Y_s，记为：

$Y_i=\sqrt{\frac{n_i-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(15\right)$

$Y_S=\sqrt{\frac{n_s-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(16\right)$

式中，P_k＝[p₁,p₂,...,p_k]，Λ_k为{λ₁,λ₂,...,λ_k}构成的对角矩阵；

Y_i和Y_s的协方差矩阵记为：

$C_i=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_i^T{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(17\right)$

$C_S=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_S^T{X}_S{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(18\right)$

两个数据集的协方差矩阵满足下式：

$C_i+C_S={{\mathrm{\Λ}}_k}^{-1/2}{P}_k^T{\mathrm{RP}}_k{{\mathrm{\Λ}}_k}^{-1/2}={{\mathrm{\Λ}}_k}^{-1/2}{P}_k^T{P}_k{\mathrm{\Λ}}_k{{\mathrm{\Λ}}_k}^{-1/2}=I---\left(19\right)$

对协方差矩阵应用特征值分解，得到

$C_i{w}_j^{\left(i\right)}={\mathrm{\λ}}_j^{\left(i\right)}{w}_j^{\left(i\right)},j=1,2,...,k---\left(20\right)$

$C_S{w}_j^{\left(s\right)}={\mathrm{\λ}}_j^{\left(s\right)}{w}_j^{\left(s\right)}---\left(21\right)$

由于和C_S的特征向量与C_i的特征向量相同，仅仅是次序相反，对相同大小的数据集X_i和X_S，如果具有相似的数据分布和关联关系，则特征值λ_j和1-λ_j均接近0.5；相反的，如果数据集描述了不同的数据分布和数据关系，则特征值分别趋近1和0。

计算各故障模式的统计量变量与待识别故障数据统计量的故障识别指数FRI的步骤为：相异度指数D可以用来评价数据集之间的相似程度，定义为：

$D=\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(22\right)$

式中，λ_j是C_i或者C_s的特征值，如果两个数据集相似程度高，则特征值接近0.5，D接近0；

为了方便故障识别，定义故障识别指数FRI为S_DA，其计算公式为：

$S_{DA}=1-D=1-\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(23\right)$

当两个数据集相似程度高时，S_DA接近1。

(五)对故障识别指数FRI进行排序，获得故障识别结果。

对故障识别指数进行排序，最大的故障识别指数为故障识别结果X_imax，其中，故障识别结果的序号为：

$i_{max}=\arg \underset{i}{max}\{S_{SPDA}\left(i\right),i=1,2,...,K\}---\left(24\right)$

真实故障模式的序号为i_real，如果i_real＝i_max，则故障识别结果是正确。

最理想的识别结果是：第i_max个FRI接近1，其他FRI均接近0。实际识别过程中，有可能会出现多个FRI都具有比较大的数值，这样就增加故障识别的难度，因为，即使i_real＝i_max，多个较大的FRI也会降低了故障识别结果的可信度。

获得故障识别结果后，为了评价不同方法的故障识别结果，通过识别误差指数REI和识别对比度RCD两个故障识别性能指标进行不同方法的故障识别效果对比。

所述识别误差指数REI的计算公式为：

$REI=\frac{1}{K}\underset{i=1}{\overset{K}{\Σ}}\left|{\mathrm{re}}_i\right|+\mathrm{\ρ}---\left(25\right)$

式中，re_i是识别误差，其计算公式为：

${\mathrm{re}}_i=\left\{\begin{array}{ccc}{S}_{SPDA}\left(i\right)-1,& for& i=i_{real}\\ S_{SPDA}\left(i\right)-0,& for& i\≠i_{real}\end{array}\right.---\left(26\right)$

式中，ρ是惩罚因子，错误的识别结果为1，正确的识别结果为0。

所述识别对比度的计算公式为:

$RCD=\mathrm{\γ}\frac{S_{SPDA}\left(i_{max}\right)-S_{SPDA}\left(i_{submax}\right)}{S_{SPDA}\left(i_{max}\right)}---\left(27\right)$

式中，S_SPDA(i_max)是最大的FRI，S_SPDA(i_submax)是次大的FRI，γ是故障识别错误参数，对于错误的识别结果，γ是-1，否则是1。

很显然，大的RCD数值意味着清晰正确的识别结果，可信度高；小的RCD数值意味着识别结果不正确，可信度低。

本实施例中，采用相异度分析DA法和主元相异度分析PDA法两种方法与本发明基于统计量主元相异度分析SPDA的故障识别方法分别对TE化工过程进行故障识别。

统计量数据窗宽度选为200，综合上述三种方法对全部21个故障测试数据集说明故障模式识别的效果。

上述三种方法的两个性能指标REI和RCD如表7和表8所示。

表7

表8

表7中列出的为REI指标的结果，表8中列出的为RCD指标的结果，在表中最好的识别结果用粗体字标出。从表7中可以看出，除了故障模式F15T和F16T，SPDA有最小的REI指标，而从表8中可以得出，除了故障模式F09T、F15T和F16T，SPDA有最大的RCD指标。通过全局分析可以看出，REI和RCD的均值进行对比得到PDA和SPDA的故障识别效果优于传统的DA方法。因此，得到结论是在TE化工过程中PDA和SPDA方法具有更加显著的故障识别效果，且SPDA方法的故障识别效果最佳。

因为识别一个故障得到一个非常小的RCD指标说服力不强，因此一个小RCD指标的故障识别可能被归为一个不清晰的故障识别结果。这里对一个0.1的阈值，也就是说，只有当RCD大于0.1，故障识别结果被认为是清晰的。否则，可以认为故障识别结果没有说服力。因此，定义显著故障识别率(CFRR)为RCD大于0.1的故障数据集数目占总的故障数据集数目的比例。分析表8中的RCD指标，对显著故障识别率(CFRR)进行对比如表9所示。

表9

	DA	PDA	SPDA
				CFRR	4.76％	47.62％	85.71％

由表9可以看出，SPDA有85.71％的CFRR指标，PDA有47.62％的CFRR指标，而DA的CFRR指标仅有4.76％。由此可见，采用本发明基于统计量主元相异度分析SPDA的故障识别方法进行故障识别的效果最好。

以上所举实施例仅用为方便举例说明本发明，并非对本发明保护范围的限制，在本发明所述技术方案范畴，所属技术领域的技术人员所作各种简单变形与修饰，均应包含在以上申请专利范围中。

Claims (8)

1.一种多变量工业过程故障识别方法，其特征在于：含有以下步骤：

(一)收集历史数据库的正常操作工况数据集X和K类已知的故障模式数据集{H_o1,H_o2,...,H_oK}，计算正常操作工况数据集的均值mean(X)和标准差std(X)，对已知的故障模式数据集{H_o1,H_o2,...,H_oK}进行标准化处理获得新故障模式数据集{H₁,H₂,...,H_K}；

(二)在各个故障模式数据集下构造数据窗，计算低阶统计量的变量均值ε_i(t)、低阶统计量的方差v_i(t)、高阶统计量的3阶中心距高阶统计量的4阶中心距高阶统计量的偏斜度γ_i(t)和高阶统计量的峰度k_i(t)六种统计量变量；所述六种统计变量分别由计算公式(3)、(4)、(5)、(6)、(7)、(8)计算获得，公式(3)、(4)、(5)、(6)、(7)、(8)的表达式如下所示：

${\mathrm{\ϵ}}_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{x}_i(t-h)---\left(3\right)$

$v_i\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2---\left(4\right)$

$c_i^{\left(3\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3---\left(5\right)$

$c_i^{\left(4\right)}\left(t\right)=\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4---\left(6\right)$

${\mathrm{\γ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^3}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(t-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^{3/2}}---\left(7\right)$

${\mathrm{\κ}}_i\left(t\right)=\frac{\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(k-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^4}{{\left(\frac{1}{w}\underset{h=0}{\overset{w-1}{\Σ}}{\[x_i(k-h)-{\mathrm{\ϵ}}_i\left(t\right)\]}^2\right)}^2}-3---\left(8\right)$

根据上面所述的六个公式(3)至(8)获得数据窗的6类变量，将其堆积为一个行向量

[ε(t)v(t)c⁽³⁾(t)c⁽⁴⁾(t)γ(t)κ(t)]，其中，ε(t)＝[ε₁(t)…ε_m(t)]，

v(t)＝[v₁(t)…v_m(t)]，

γ(t)＝[γ₁(t)…γ_m(t)]，κ(t)＝[κ₁(t)…κ_m(t)]；

将各个故障模式数据集分割为一系列如公式(9)所示的数据窗H_it(w≤t≤n)，公式(9)表达式如下：

$H_t=\left[\begin{array}{c}{H}_{iw}\\ .\\ .\\ .\\ H_{it}\\ .\\ .\\ .\\ H_{in}\end{array}\right]---\left(9\right)$

针对各故障模式数据集的每一个数据窗H_it(w≤t≤n)计算其统计量，进而获得各个故障模式数据集的统计量矩阵，记为：

$X_i=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right],(i=1,2,...,K)---\left(10\right);$

(三)检测过程故障，收集实时故障数据S，进行标准化处理；

(四)在步骤(三)的基础上执行统计量主元相异度分析，计算待识别故障数据集和已知故障模式数据集之间的故障识别指数FRI；执行统计量主元相异度分析的步骤为：首先构造实时数据S的窗宽为w的数据窗，根据公式(3)至(8)计算待识别故障数据集的统计量，记为：

$X_S=\left[\begin{array}{cccccc}\mathrm{\ϵ}\left(w\right)& v\left(w\right)& c^{\left(3\right)}\left(w\right)& c^{\left(4\right)}\left(w\right)& \mathrm{\γ}\left(w\right)& \mathrm{\κ}\left(w\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(t\right)& v\left(t\right)& c^{\left(3\right)}\left(t\right)& c^{\left(4\right)}\left(t\right)& \mathrm{\γ}\left(t\right)& \mathrm{\κ}\left(t\right)\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \mathrm{\ϵ}\left(n\right)& v\left(n\right)& c^{\left(3\right)}\left(n\right)& c^{\left(4\right)}\left(n\right)& \mathrm{\γ}\left(n\right)& \mathrm{\κ}\left(n\right)\end{array}\right]---\left(11\right)$

对于已知故障模式数据集的统计量矩阵和待识别故障数据统计量分别包含m个变量的n_i(i＝1,...,K)和n_s个样本，其协方差阵计算公式为：

$R_i=\frac{1}{n_i-1}{X}_i^T{X}_i,R_S=\frac{1}{n_s-1}{X}_S^T{X}_S---\left(12\right)$

对第i个故障模式，构造总体统计量数据集的协方差矩阵，记为：

$R=\frac{1}{n_i+n_s-1}{X}^{T}X---\left(13\right)$

对其开展特征值分解，获得：

RP＝PΛ(14)

式中，P是特征向量矩阵，Λ(λ₁≥λ₂≥...≥λ_m)是特征值矩阵；

保留L个特征值和对应的特征向量p₁,p₂,...,p_L，对矩阵X_i和X_S进行主元线性变换，得到Y_i和Y_s，记为：

$Y_i=\sqrt{\frac{n_i-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(15\right)$

$Y_S=\sqrt{\frac{n_s-1}{n_i+n_s-1}}{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(16\right)$

式中，P_k＝[p₁,p₂,...,p_k]，Λ_k为{λ₁,λ₂,...,λ_k}构成的对角矩阵；

Y_i和Y_s的协方差矩阵记为：

$C_i=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_i^T{X}_i{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(17\right)$

$C_S=\frac{1}{n_1+n_s-1}{{\mathrm{\Λ}}_L}^{-1/2}{P}_L^T{X}_S^T{X}_S{P}_L{{\mathrm{\Λ}}_L}^{-1/2}---\left(18\right)$

两个数据集的协方差矩阵满足下式：

$C_i+C_S={{\mathrm{\Λ}}_L}^{-1/2}{P}_k^T{\mathrm{RP}}_k{{\mathrm{\Λ}}_L}^{-1/2}={{\mathrm{\Λ}}_L}^{-1/2}{P}_k^T{P}_k{\mathrm{\Λ}}_k{{\mathrm{\Λ}}_L}^{-1/2}=I---\left(19\right)$

对协方差矩阵应用特征值分解，得到

$C_i{w}_j^{\left(i\right)}={\mathrm{\λ}}_j^{\left(i\right)}{w}_j^{\left(i\right)},j=1,2,...,k---\left(20\right)$

$C_S{w}_j^{\left(s\right)}={\mathrm{\λ}}_j^{\left(s\right)}{w}_j^{\left(s\right)}---\left(21\right)$

由于和C_S的特征向量与C_i的特征向量相同，仅仅是次序相反，对相同大小的数据集X_i和X_S，如果具有相似的数据分布和关联关系，则特征值λ_j和1-λ_j均接近0.5；相反的，如果数据集描述了不同的数据分布和数据关系，则特征值分别趋近1和0；

计算各故障模式的统计量变量与待识别故障数据统计量的故障识别指数FRI的步骤为：相异度指数D可以用来评价数据集之间的相似程度，定义为：

$D=\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(22\right)$

式中，λ_j是C_i或者C_s的特征值，如果两个数据集相似程度高，则特征值接近0.5，D接近0；

为了方便故障识别，定义故障识别指数FRI为S_DA，其计算公式为：

$S_{DA}=1-D=1-\frac{4}{L}\underset{j=1}{\overset{L}{\Σ}}{({\mathrm{\λ}}_j-0.5)}^2---\left(23\right)$

当两个数据集相似程度高时，S_DA接近1；

(五)对故障识别指数FRI进行排序，获得故障识别结果。

2.根据权利要求1所述的多变量工业过程故障识别方法，其特征在于：所述步骤(一)中，通过公式(1)对故障模式数据进行标准化处理，公式(1)的表达式如下：

H_i＝(H_oi-mean(X))/std(X),(i＝1,2,...,K)(1)

经上述公式(1)标准化处理后即可获得新故障模式数据集{H₁,H₂,...,H_K}。

3.根据权利要求1所述的多变量工业过程故障识别方法，其特征在于：所述步骤(二)中，构造数据窗的步骤为：对任意t时刻的数据窗记为：

式中，w为数据窗宽度。

4.根据权利要求1所述的多变量工业过程故障识别方法，其特征在于：所述步骤(三)中，过程检测到故障后，收集实时故障数据集S_o，使用公式S＝(S_o-mean(X))/std(X)进行标准化，上式中，X为正常操作工况数据集。

5.根据权利要求1所述的多变量工业过程故障识别方法，其特征在于：所述步骤(五)中，对故障识别指数进行排序，最大的故障识别指数为故障识别结果X_imax，其中，故障识别结果的序号为：

$i_{max}=\arg \underset{i}{max}\{S_{SPDA}\left(i\right),i=1,2,...,K\}---\left(24\right)$

真实故障模式的序号为i_real，如果i_real＝i_max，则故障识别结果是正确的。

6.根据权利要求1所述的多变量工业过程故障识别方法，其特征在于：所述步骤(五)中，获得故障识别结果后，为了评价不同方法的故障识别结果，通过识别误差指数REI和识别对比度RCD两个故障识别性能指标进行不同方法的故障识别效果对比。

7.根据权利要求6所述的多变量工业过程故障识别方法，其特征在于：所述识别误差指数REI的计算公式为：

$REI=\frac{1}{K}\underset{i=1}{\overset{K}{\Σ}}\left|{\mathrm{re}}_i\right|+\mathrm{\ρ}---\left(25\right)$

式中，re_i是识别误差，其计算公式为：

${\mathrm{re}}_i=\left\{\begin{array}{cc}{S}_{SPDA}\left(i\right)-1,& fori=i_{real}\\ S_{SPDA}\left(i\right)-0,& fori\≠i_{real}\end{array}\right.---\left(26\right)$

式中，ρ是惩罚因子，错误的识别结果为1，正确的识别结果为0。

8.根据权利要求6所述的多变量工业过程故障识别方法，其特征在于：所述识别对比度RCD的计算公式为:

$RCD=\mathrm{\γ}\frac{S_{SPDA}\left(i_{max}\right)-S_{SPDA}\left(i_{submax}\right)}{S_{SPDA}\left(i_{max}\right)}---\left(27\right)$

式中，S_SPDA(i_max)是最大的FRI，S_SPDA(i_submax)是次大的FRI，γ是故障识别错误参数，对于错误的识别结果，γ是-1，否则是1。