基于正交局部慢性特征的故障检测方法
张展博1, 王振雷1, 王昕2    
1. 华东理工大学 化工过程先进控制和优化技术教育部重点实验室, 上海 200237;
2. 上海交通大学 电工与电子技术中心, 上海 200240
摘要:为提高化工行业中数据驱动故障检测的效果,该文针对实际工业系统中闭环控制导致的过程动态特性和数据流形中蕴含的局部信息,提出了一种基于局部时空正则慢特征提取(local time-space regularized slow feature extraction,LTSS)的方法进行故障检测。首先,构造基于局部时空正则的目标函数得到投影矩阵,进而得到预提取特征S,则S张成的空间中包含了静态信息,而S的一阶差分张成的空间中包含了动态信息。其次,基于独立成分分析(independent components analysis,ICA)方法,分别为2个空间构建对应的S2和SPE统计量进行监控,用于实时故障检测。在TE(Tennessee Eastman)过程上的案例研究可以证明所提方法的有效性。
关键词故障检测    动态特性    局部相似性    慢性先验    独立成分分析(ICA)    TE过程    
Fault detection based on orthogonal local slow features
ZHANG Zhanbo1, WANG Zhenlei1, WANG Xin2    
1. Key Laboratory of Advanced Control and Optimization for Chemical Processes, East China University ofScience and Technology, Shanghai 200237, China;
2. Center of Electrical & Electronic Technology, Shanghai Jiao Tong University, Shanghai 200240, China
Abstract: A local time-space regularized slow feature extraction method was developed to improve data-driven fault detection in the chemical industry based on the process dynamics of closed-loop control systems and the local information contained in the data manifold. An objective function was defined based on the local time-space term to obtain a projection matrix and the pre-extraction feature, S. The span of S contains the static information, while the first derivative of the span of S contains the dynamic information. An independent component analysis was used to obtain statistics for S2 and SPE for both spaces for real-time fault detection. A case study on the Tennessee Eastman process shows the validity of this method.
Key words: fault detection    process dynamics    local similarity    slow feature    independent component analysis (ICA)    Tennessee Eastman (TE) process    

为了提升工业产品的质量和生产效率,现代工业生产过程日益复杂,依靠操作人员监测整个过程的运行状态变得越来越困难; 另一方面,现代工业过程系统结构也日趋复杂,很难建立精确的机理模型,因此基于机理模型的故障检测方法受到很大限制。随着计算机技术在工业过程中的广泛应用,大量的过程数据被采集储存,基于数据驱动的多元统计方法被广泛地应用在工业过程中故障检测和故障诊断任务中,特别是基于主元分析(principal components analysis, PCA)、偏最小二乘(partial least squares, PLS)和独立成分分析(independent component analysis, ICA)的技术及其改进方法[1-5]得到了大量的研究。

PCA方法是一种提取过程数据信息的有效方法,但该方法存在数据满足高斯分布的假设条件,而这在实际工业数据中往往不成立,由此引入了ICA方法,考虑非高斯特征。但PCA和ICA本质上是一个静态方法,忽略了过程数据中的时间维度,这会在特征提取过程中造成动态信息的损失。为解决该问题,Qin等[6-7]提出动态PCA(dynamic principal components analysis, DPCA)方法,将过程数据矩阵进行增广,但这种方法大大增加了数据向量的维度,容易受到“维度诅咒”效应的影响。为此,Zhao[8]在动态图嵌入(dynamic graph embedding, DGE)方法中,使用了二维矩阵分解策略,其结果表明,基于矩阵的特征提取算法比基于矢量的特征提取算法计算效率更高,但这些方法都并未显式地对时间信息进行建模。2012年,慢特征分析(slow feature analysis, SFA)方法被引入故障检测领域[9]。文[10]提出概率慢特征分析(probabilistic slow feature analysis, PSFA)方法,Guo等[11]在故障检测任务中引入这种方法,将PSFA表达为状态空间模型,阐明其实质是一个隐变量符合一阶Markov过程的动态模型,显式对时间信息进行建模,可以用来提取历史数据中包含的动态信息。这非常适合于数据依时间顺序产生的故障检测任务。为了充分利用提取到的动态信息,Gu等[12]将PCA方法与SFA方法融合,提出SFPCA (slow feature principal component analysis)方法,并得到一定效果。

SFA方法是一个全局算法,没有考虑数据的局部信息。实际上,数据的局部特征在故障检测任务中至关重要,LLE[13](locally linear embedding)、LPP[14](locality preserving projections)、NPE[15](neighborhood preserving embedding)等是经典的能够捕捉到数据中局部信息的方法,Miao等[16]考虑时间邻域信息,扩展了NPE方法,提出时间邻域保持嵌入算法,类似地,Tong等[17]利用更多的局部信息,提出NoMNPE (nonlocal and multiple neighborhoods preserving embedding)方法,都取得了较好的效果,但这些方法中并未对历史数据集中的动态信息显式建模,这在时序数据处理任务中是不合理的。

为了得到更佳的故障检测效果,本文同时对数据集中的局部信息和动态信息进行建模,提出一种基于LTSS的故障检测方法,该方法首先基于局部时空正则构造优化目标,从历史数据中提取慢特征S,所提取到的特征中包含历史数据中的局部信息和动态信息,可以认为慢特征S张成的空间中主要包含数据集中的工作点信息,而慢特征S的一阶差分$\dot S$张成的空间中主要包含数据集中的动态信息,最后分别为2个空间建立S2、SPE及Sd2、SPEd统计量,监视2个空间中的数据分布,检测当前工况是否发生故障。

1 LTSS-ICA故障检测方法

Sprekeler[18]阐明了SFA与Laplace图嵌入(Laplacian eigenmaps, LEMs)方法的联系,而LEMs方法是基于流形假设下的一种重要降维方法,且已经在化工行业数据的分析中得到了大量研究,取得了不错的实验结果。本文的目标在于构建慢特征的同时,引入局部时空正则,以提取到原始数据空间中的更多信息。

1.1 基于局部时空正则的慢特征提取 1.1.1 数据的局部几何特征

化工行业的数据大都通过布置在现场的传感器获得,依时间顺序变化,描述所监测量的变化趋势,因此时间上相邻的数据蕴含系统动态信息,另一方面,良好的控制系统中在某一小段时间内采集到的数据在相应的维度上应具有紧凑的聚集特性,因此,本文考虑数据在空间和时间上的局部相似性。存在多种不同的方法可以构建数据内的局部相似性,本文采用文[8]使用的度量方法:

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{ij}}( {\rm{Dist}}{ _{ij}};{\beta _1},{\beta _2}) = \\ \left\{ \begin{array}{l} {\rm{exp}}\left( {\frac{{ - {\rm{Dist}}{ _{ij}}}}{{{\beta _1}}}} \right){\rm{exp}}\left( {\frac{{ - {{(i - j)}^2}}}{{{\beta _2}}}} \right),{\rm{ 如果 }}{\mathit{\boldsymbol{x}}_i},{\mathit{\boldsymbol{x}}_j}{\rm{互为}}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k{\kern 1pt} {\kern 1pt} {\rm{近邻;}}\\ {\rm{0,其他情况}}{\rm{.}} \end{array} \right. \end{array} $ (1)

其中,Distij=‖xixj22xi表示第i个数据,i, j分别为数据序号,表示当前数据在数据序列中的位置信息; β1β2为超参数,用于调节距离邻域结构和时间邻域结构。

1.1.2 基于局部时空正则的慢特征提取

慢性先验假设采集到的数据由隐含的变化较慢的因素驱动产生,设时刻t采集到的m维数据表示为x(t)=[x1(t), x2(t), …, xm(t)]T, t∈[t0, t1],LTSS方法的目标是找到映射函数g(x)=[g1(x), g2(x), …, gJ(x)]T,使得输入数据经过该J维向量函数变换后得到的J维特征为s(t)=[s1(t), s2(t), …, sJ(t)]T,其中si(t)=gi(x(t))。期望得到的J维特征s变化尽量缓慢,该目标可表示为

$ {J_1} = \sum\limits_t {\left\| {\mathit{\boldsymbol{\dot s}}(t)} \right\|_2^2} . $ (2)

其中: ‖·‖2表示向量的二范数,$\mathit{\boldsymbol{\dot s}}$(t)(表示t时刻数据对应的慢特征)为s(t)的一阶差分,用来衡量其变化趋势,$\mathit{\boldsymbol{\dot s}}$(t)=s(t)-s(t-1)。

另一个合理的期望是,由原始数据中较为相似数据得到的慢特征s的变化速率仍较为相近,该期望强调了对原始数据中动态信息的描述,可表示为如下的形式:

$ {J_2} = \sum\limits_{{t_i}} {\sum\limits_{{t_j}} {{C_{ij}}} } \left\| {\mathit{\boldsymbol{\dot s}}({t_i}) - \mathit{\boldsymbol{\dot s}}({t_j})} \right\|_2^2. $ (3)

其中,Cij由式(1)得到。

为了能够提取到原始数据中蕴含的尽可能多的信息,希望得到的慢特征s也能够体现出原始数据中蕴含的局部信息,目标是由原始数据中较为相似数据得到的慢特征s仍较为相似,可表示为如下形式:

$ {J_3} = \sum\limits_{{t_i}} {\sum\limits_{{t_j}} {{C_{ij}}} } {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\| {\mathit{\boldsymbol{s}}({t_i}) - \mathit{\boldsymbol{s}}({t_j})} \right\|_2^2. $ (4)

整合以上3个目标,得到最终的目标函数:

$ J = {\rm{arg}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \mathop {{\rm{min}}}\limits_g ({J_1} + \lambda {J_2} + \gamma {J_3}). $ (5)

目标函数式(5)是一个泛函问题,难以求解,因此这里考虑函数g(x)为线性函数的简单情况。假设由原始数据组成的矩阵为X,则在线性情况下,数据矩阵X对应的慢特征矩阵S可通过S=UTX得到,可以重写目标函数为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} J = {\rm{arg}}{\kern 1pt} {\kern 1pt} \mathop {{\rm{min}}}\limits_U (\sum\limits_t {\left\| {{\mathit{\boldsymbol{U}}^{\rm{T}}}{\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{\dot x}}(t)} \right\|_2^2 + } \\ \lambda \sum\limits_{{t_i}} {\sum\limits_{{t_j}} {{C_{ij}}} } \left\| {{\mathit{\boldsymbol{U}}^{\rm{T}}}{\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{\dot x}}({t_i}) - {\mathit{\boldsymbol{U}}^{\rm{T}}}{\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{\dot x}}({t_j})} \right\|_2^2 + \\ \left. {\gamma \sum\limits_{{t_i}} {\sum\limits_{{t_j}} {{C_{ij}}} } \left\| {{\mathit{\boldsymbol{U}}^{\rm{T}}}{\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{x}}({t_i}) - {\mathit{\boldsymbol{U}}^{\rm{T}}}{\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{x}}({t_j})} \right\|_2^2} \right) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{arg}}{\kern 1pt} {\kern 1pt} \mathop {{\rm{min}}}\limits_U {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{tr}} \left( {{\mathit{\boldsymbol{U}}^{\rm{T}}}\left[ {\mathit{\boldsymbol{\dot X}}(\mathit{\boldsymbol{I}} + \lambda (\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{C}})){{\mathit{\boldsymbol{\dot X}}}^{\rm{T}}} + } \right.} \right.\\ \left. {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \gamma \mathit{\boldsymbol{X}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{C}}){\mathit{\boldsymbol{X}}^{\rm{T}}}} \right]\mathit{\boldsymbol{U}}} \right). \end{array} $ (6)

设目标函数式(6)的最优解为U*,则只需为U*乘上小于1的标量,就能使目标函数值更小,因此需引入约束条件来约束U,使得存在最优解U*,考虑到矩阵D度量了每个点与邻域内其他点的相似程度,若一个点与越多的点相似,则该点会起到越重要的作用,因此引入约束条件如下:

$ \mathit{\boldsymbol{u}}_i^{\rm{T}}(\mathit{\boldsymbol{\dot XD\dot X}}){\mathit{\boldsymbol{u}}_i} = 1. $ (7)

其中,ui为矩阵U的第i列。

另外,文[19]表明保证矩阵U中不同列间的正交性,更有利于提取原始数据中的信息。因此,引入另一个约束条件:

$ \mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{u}}_1} = 0, \cdots ,\mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{u}}_{i - 1}} = 0. $ (8)

为了方便,约定以下表达式:

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{Q}} = \mathit{\boldsymbol{\dot XD}}{{\mathit{\boldsymbol{\dot X}}}^{\rm{T}}},{\mathit{\boldsymbol{U}}^{i - 1}} = [{\mathit{\boldsymbol{u}}_1},{\mathit{\boldsymbol{u}}_2}, \cdots ,{\mathit{\boldsymbol{u}}_{i - 1}}],}\\ {{\mathit{\boldsymbol{B}}^{i - 1}} = {{[{\mathit{\boldsymbol{U}}^{i - 1}}]}^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{ - 1}}{\mathit{\boldsymbol{U}}^{i - 1}},}\\ {\mathit{\boldsymbol{H}} = \mathit{\boldsymbol{\dot X}}(\mathit{\boldsymbol{I}} + \lambda (\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{C}})){{\mathit{\boldsymbol{\dot X}}}^{\rm{T}}} + \gamma \mathit{\boldsymbol{X}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{C}}){\mathit{\boldsymbol{X}}^{\rm{T}}},}\\ {{\mathit{\boldsymbol{M}}^i} = \{ \mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{Q}}^{ - 1}}{\mathit{\boldsymbol{U}}^{i - 1}}{{({\mathit{\boldsymbol{B}}^{i - 1}})}^{ - 1}}{{({\mathit{\boldsymbol{U}}^{i - 1}})}^{\rm{T}}}\} {\mathit{\boldsymbol{Q}}^{ - 1}}\mathit{\boldsymbol{H}}.} \end{array}} \right. $ (9)

则该优化问题的最优解为:u1Hu=λQu的最小特征值对应的特征向量,ui(i>1)为Mi的最小特征值对应的特征向量,详细推导过程参见文[19],这里不再给出。根据S=UTX,即可得到预提取特征矩阵。

1.2 基于预提取慢特征的故障检测方法 1.2.1 基于ICA的分解

n维连续型随机变量M=[m1, m2, …, mn]T,且M的概率密度分布函数为fM(m),存在n阶方阵A可逆,使得Y=AM,则Y的概率密度函数与M的概率密度函数满足以下关系:

$ {f_Y}(y) = \frac{1}{{| {\rm{det}} (\mathit{\boldsymbol{A}})|}}{f_X}({\mathit{\boldsymbol{A}}^{ - 1}}y). $ (10)

即,可以通过分析随机变量Y的概率分布间接得到随机变量M的概率分布特征,因此引入ICA方法,对特征SS张成的空间进行再分解,构造统计量进行监测,达到近似监测原数据空间静态信息和动态信息分布的目的。

由于不同维度的特征所包含信息的重要程度有所不同,变化较慢的特征一般被认为包含更本质的信息,而变化相对较快的特征中可能混有采集数据中的噪声,特征的慢度可近似由对应的广义特征值所衡量,为此,本文对不同慢度的特征进行加权,加权策略如下:

$ {w_i} = \frac{{{{\rm{e}}^{ - h{\lambda _i}}}}}{{\sum\limits_j {{{\rm{e}}^{ - h{\lambda _j}}}} }}. $ (11)

即,加权后的特征为s′=sw。其中: wi为对第i个慢特征的加权系数; h为调节参数,调节权重的分散程度; λi为第i个慢特征对应的特征值; w=[w1, w2, …, wn]T,⊙表示Hadamard积。

最终对s′进行分解得到

$ {\mathit{\boldsymbol{s}}^\prime } = {\mathit{\boldsymbol{W}}_1}{\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{s}}_{{\rm{ non\_gauss }}}} + \mathit{\boldsymbol{e}}. $ (12)

其中,snon_gauss为非高斯独立成分,e为残差,W1为混合矩阵,由ICA算法确定,对表征动态特性的$\mathit{\boldsymbol{\dot s}}$张成的空间进行类似分解,得到

$ {\mathit{\boldsymbol{\dot s}}^\prime } = {\mathit{\boldsymbol{W}}_2}{\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{s}}_{{\rm{ non\_gauss\_d }}}} + {\mathit{\boldsymbol{e}}_{\rm{d}}}. $ (13)

其中,$\mathit{\boldsymbol{\dot s}}$′为s′的一阶差分,特别的对于t时刻的特征$\mathit{\boldsymbol{\dot s}}$′(t)可计算为:$\mathit{\boldsymbol{\dot s}}$′(t)=s′(t)-s′(t-1)。同上所述,W2为混合矩阵,由ICA算法确定。

1.2.2 故障检测策略

为监视当前工作状态,定义S2和SPE统计量:

$ {{S^2} = {{\left\| {{\mathit{\boldsymbol{s}}_{{\rm{ non\_gauss }}}}} \right\|}_2},} $ (14)
$ {{\rm{SPE}} = {{\left\| \mathit{\boldsymbol{e}} \right\|}_2}.} $ (15)

当上述2个统计量发出故障警告时,表示当前系统内监控变量间的互相关关系可能发生了变化,而提取到的动态信息可以给出当前工况的进一步信息,以帮助现场操作人员快速定位到故障源,定义如下监测当前系统动态信息的统计量:

$ {S_{\rm{d}}^2 = {{\left\| {{\mathit{\boldsymbol{s}}_{{\rm{non\_gaus}}{{\rm{s}}_{\rm{ - }}}{\rm{d}}}}} \right\|}_2}.} $ (16)
$ {{\rm{SP}}{{\rm{E}}_{\rm{d}}} = {{\left\| {{\mathit{\boldsymbol{e}}_{\rm{d}}}} \right\|}_2}.} $ (17)

为了能够在线实时监测过程工作点是否发生变化,需要给出S2和SPE统计量的监控限,另外,为了监测系统动态信息的变化,Sd2和SPEd统计量也需要给出相应监控限。然而,由于工业过程中控制系统内物理量间的关系错综复杂,很难保证测得的过程数据符合某一特定分布,因此这里对原始数据的分布不做任何假定,最终也无法确定这些统计量所满足的分布。核密度估计(kernel density estimation, KDE)方法是估计某一随机变量的概率密度分布函数的强有力工具,这里使用KDE方法确定上述统计量的监控限[20]

该方法的执行如大多数故障检测方法一样分为离线建模阶段和在线检测阶段,具体工作流程如图 1所示。在离线建模阶段使用正常工况下的数据拟合LTSS-ICA方法中的参数,并基于KDE方法确定相应控制限(本文置信区间设为95%); 在线监测阶段,根据式(14)—(17)计算实时采集到的样本的相关统计量,并与设定好的控制限进行对比,此时可分为2种情况考虑:

图 1 LTSS-ICA流程图

1) 若表征系统工作状态的任一统计量连续超出控制限,且表征系统动态信息的任一统计量也连续超出控制限,则认为当前系统中发生了难以控制的故障。

2) 若表征系统工作状态的任一统计量连续超出控制限,但表征系统动态信息的统计量在一段连续时间内超出控制限,而后又恢复到控制限以下,则认为当前系统中可能发生了输入小幅阶跃变化等不影响系统正常工作的情况。

2 案例研究

TE过程实验平台最早由Downs和Vogel开发,模拟Eastman公司的实际化工过程(如图 2所示)。该仿真过程包含21类预定义故障类型,涉及阶跃、随机、缓慢漂移、阀门黏滞等故障类型,还有一系列未知类型故障,详细故障描述可参考文[21]。本文使用数据可从http://web.mit.edu/braatzgroup/links.html下载,其中包含22种工况,即正常工况和预定义的21种故障工况,每种工况下包含一个训练集和一个测试集,正常工况下的训练集包含500个样本,用于训练模型,而在故障工况下,训练集包含480个样本,测试集包含960个样本,用于测试得到的模型,其中故障工况测试集在第160个样本后引入故障。

图 2 TE流程示意图

为证明所提方法的有效性,本文与Zhao所提的DGE[8]方法进行对比。与本文所提方法相比,DGE方法也同时考虑了局部信息和动态信息,具有一定可比性。

故障5为阶跃型故障,故障描述为冷凝器冷却水的入口温度发生变化,本文所提方法在特征空间和特征一阶差分空间建立的控制图分别如图 3所示,从图中可以看出,S2和SPE统计量持续超出控制限,而S2和SPE统计量在发生故障时迅速给出警报信息,但一段时间后又恢复到控制限以下,这说明该阶跃发生后,系统经过闭环控制作用最终得以稳定,而DGE方法不能给出这样的信息。

图 3 故障5下监测结果

故障19为未知类型故障,从图 4可以看出,在引入故障后,LTSS-ICA方法中S2统计量能够很好地检测到系统内静态信息的变化,而Sd2和SPEd统计量在发生故障后也表现出与故障前完全不一样的特征,这说明系统内可能发生了闭环难以控制的故障。而DGE方法提出的2种统计量均没有实现良好的故障检测。

图 4 故障19下监测结果

在TE过程所有预定义故障工况下,使用ICA、DGE和LTSS-ICA方法,得到的相应统计量的故障检测率(fault detection ratio,FDR)如表 1所示,其值越大表现越好(LTSS-ICA方法中,监测动态信息的统计量并未列出),从表中可以看出,故障3、9、15在3种方法中都没有被很好地检测到,这是因为这3种故障发生后,测得的数据发生的变化幅度过于微小,淹没在测量噪声之中,实际上现已提出的方法中,除了专门针对微小故障实施类似累积和控制图策略的方法[22]能够检测到,其他方法大都无法捕获到这3种情况下的故障信息。对于其他故障情况,LTSS-ICA方法总体上都有同样优秀或更好的表现,且LTSS-ICA方法在检测出静态信息发生变化后,可以根据Sd2和SPEd统计量获得关于系统动态特性变化的信息,有益于操作人员更快地洞悉引起故障警报的原因。

表 1 相关统计量在TE数据集上的FDR
故障编号 ICA DGE LTSS-ICA
S2 SPE S2 SPE S2 SPE
1 0.9975 0.9900 0.995 0.9937 0.9975 0.9988
2 0.9888 0.9862 0.9889 0.9837 0.9912 0.9725
3 0.0488 0.0088 0.0625 0.0463 0.0362 0.0213
4 0.9975 0.9225 0.3283 0.9962 1 0.7500
5 1 0.9988 0.9937 0.5963 1 0.9988
6 1 1 1 0.9987 1 1
7 1 0.3887 0.6613 1 1 1
8 0.9825 0.7300 1 0.9825 0.9838 0.8287
9 0.0100 0.0200 0 0.0225 0.0175 0.0625
10 0.9087 0.5825 0.9262 0.7300 0.9287 0.7138
11 0.7588 03725 0.3125 0.9525 0.8137 0.5188
12 0.9988 0.9575 0.9550 0 0.9988 0.9788
13 0.9563 0.9250 1 1 0.9625 09425
14 1 0.8712 0.0725 1 1 1
15 0.1975 0.0413 0.3975 0.5338 0.2387 0.1025
16 0.9100 0.6412 0.1800 0.2600 0.9237 0.7113
17 0.9587 0.8838 0.6338 0.9612 0.9625 0.8925
18 0.9100 0.9050 0.8887 0.9262 0.9213 0.9063
19 0.9063 0.6038 0.0250 0.0925 0.9387 0.6425
20 0.9125 0.6075 0.8613 0.8475 0.9150 0.8063
21 0.5787 0.3525 0.6038 0.4550 0.6262 0.2225

3 结论

本文提出了局部时空正则的慢特征提取-独立成分分析(LTSS-ICA)方法用于故障检测,并在TE过程数据集下进行了验证。首先基于慢性先验项和时空邻域正则项构建了一个无监督的优化目标,用于初步特征提取,在提取到的特征SS的一阶差分$\dot S$张成的空间中使用ICA方法,建立相应的S2和SPE统计量并根据KDE方法确定每一统计量的控制限,用于近似监控原始数据中工作点信息和动态信息的分布,超出控制限时及时给出警报。实验结果表明,所提出的LTSS-ICA方法具有很好的表现,能够提升故障检测的性能,且能提供其他2种方法不能提供的对动态信息的监测,相比DGE方法有更好的适用性。但在具有高度非线性的过程中,LTSS-ICA方法可能会表现出不太满意的效果,且本文尚未给出故障诊断的方法,这正是下一步研究中所要实验和探索的方向。

参考文献
[1]
DENG X G, TIAN X M, CHEN S, et al. Nonlinear process fault diagnosis based on serial principal component analysis[J]. IEEE Transactions on Neural Networks and Learning Systems, 2018, 29(3): 560-572.
[2]
YIN S, XIE X C, SUN W. A nonlinear process monitoring approach with locally weighted learning of available data[J]. IEEE Transactions on Industrial Electronics, 2017, 64(2): 1507-1516. DOI:10.1109/TIE.2016.2612161
[3]
YIN S, ZHU X P, KAYNAK O. Improved PLS focused on key-performance-indicator-related fault diagnosis[J]. IEEE Transactions on Industrial Electronics, 2015, 62(3): 1651-1658. DOI:10.1109/TIE.2014.2345331
[4]
CAI L F, TIAN X M, CHEN S. Monitoring nonlinear and non-Gaussian processes using Gaussian mixture model-based weighted kernel independent component analysis[J]. IEEE Transactions on Neural Networks and Learning Systems, 2017, 28(1): 122-135.
[5]
CHEN H T, JIANG B, LU N Y, et al. Deep PCA based real-time incipient fault detection and diagnosis methodology for electrical drive in high-speed trains[J]. IEEE Transactions on Vehicular Technology, 2018, 67(6): 4819-4830. DOI:10.1109/TVT.2018.2818538
[6]
LI W H, QIN S J. Consistent dynamic PCA based on errors-in-variables subspace identification[J]. Journal of Process Control, 2001, 11(6): 661-678. DOI:10.1016/S0959-1524(00)00041-X
[7]
LI G, LIU B S, QIN S J, et al. Dynamic latent variable modeling for statistical process monitoring[J]. IFAC Proceedings Volumes, 2011, 44(11): 12886-12891.
[8]
ZHAO H T. Dynamic graph embedding for fault detection[J]. Computers & Chemical Engineering, 2018, 117(2): 359-371.
[9]
DENG X G, TIAN X M, HU X Y. Nonlinear process fault diagnosis based on slow feature analysis[C]//Proceedings of the 10th World Congress on Intelligent Control and Automation. Beijing, China: IEEE, 2012: 3152-3156.
[10]
TURNER R, SAHANI M. A maximum-likelihood interpretation for slow feature analysis[J]. Neural Computation, 2007, 19(4): 1022-1038. DOI:10.1162/neco.2007.19.4.1022
[11]
GUO F H, SHANG C, HUANG B, et al. Monitoring of operating point and process dynamics via probabilistic slow feature analysis[J]. Chemometrics and Intelligent Laboratory Systems, 2016, 151(15): 115-125.
[12]
GU S M, LIU Y L, ZHANG N. Fault diagnosis in Tennessee Eastman process using slow feature principal component analysis[J]. Recent Innovations in Chemical Engineering, 2016, 9(1): 49-61.
[13]
ROWEIS S T, SAUL L K. Nonlinear dimensionality reduction by locally linear embedding[J]. Science, 2000, 290(5500): 2323-2326. DOI:10.1126/science.290.5500.2323
[14]
HE X F, NIYOGI P. Locality preserving projections[J]. Advances in Neural Information Processing Systems, 2003, 16(1): 186-197.
[15]
HE X F, CAI D, YAN S C, et al. Neighborhood preserving embedding[C]//Tenth IEEE International Conference on Computer Vision. Beijing, China: IEEE, 2005: 1208-1213.
[16]
MIAO A M, GE Z Q, SONG Z H, et al. Time neighborhood preserving embedding model and its application for fault detection[J]. Industrial & Engineering Chemistry Research, 2013, 52(38): 13717-13729.
[17]
TONG C D, LAN T, SHI X H, et al. Statistical process monitoring based on nonlocal and multiple neighborhoods preserving embedding model[J]. Journal of Process Control, 2018, 65: 34-40. DOI:10.1016/j.jprocont.2017.10.009
[18]
SPREKELER H. On the relation of slow feature analysis and Laplacian eigenmaps[J]. Neural Computation, 2011, 23(12): 3287-3302. DOI:10.1162/NECO_a_00214
[19]
CAI D, HE X F, HAN J W, et al. Orthogonal laplacianfaces for face recognition[J]. IEEE Transactions on Image Processing, 2006, 15(11): 3608-3614. DOI:10.1109/TIP.2006.881945
[20]
BOTEV Z I, GROTOWSKI J F, KROESE D P. Kernel density estimation via diffusion[J]. The Annals of Statistics, 2010, 38(5): 2916-2957. DOI:10.1214/10-AOS799
[21]
LYMAN P R, GEORGAKIS C. Plant-wide control of the Tennessee Eastman problem[J]. Computers & Chemical Engineering, 1995, 19(3): 321-331.
[22]
DU Y C, DU D P. Fault detection and diagnosis using empirical mode decomposition based principal component analysis[J]. Computers & Chemical Engineering, 2018, 115(12): 1-21.