基于规范变量残差的化工过程微小故障检测与诊断

商亮亮 陆智林 文传博 邱爱兵

(1.南通大学电气工程学院,江苏南通 226019;2.上海电机学院电气学院,上海 201306)

1 引言

高度相关的子系统与控制单元导致工业过程变得越来越复杂[1].大多数工业过程对可靠性和安全性有更加迫切的需求[2-3].随着数据驱动技术的迅速发展,多变量统计过程监测(multivariate statistic process monitoring,MSPM)方法在工业和学术界引起了极大的关注[4-6].3种最常用的数据驱动监测方法是主成分分析(principal component analysis,PCA)[7-8],偏最小二乘(partial least squares,PLS)[9-10]和规范变量分析(canonical variate analysis,CVA)[11-13].其中,基于CVA的监测方法在动态过程监测中显示了它的优势.然而,大多数已发表的文献集中于检测幅值相对较高且变化较大的突变故障[13].微小故障的主要特征是其变化发展缓慢[14],且用常规故障检测方法不容易检测到.微小故障的诊断方法始于贝尔实验室提出的故障树分析法[15].李娟等[16]进行了微小故障诊断方法综述,将微小故障诊断分为定性诊断、定量诊断和半定性半定量诊断的方法.定性诊断方法主要由图论方法和专家系统,故障树分析法便是图论方法的一种;定量诊断方法分为基于模型的方法和数据驱动的方法.

许多学者对数据驱动检测方法进行了改进.一类基于残差指标增强检测灵敏度的方法:Pilario和Cao等[13]提出了规范变量差异分析(canonical variate dissimilarity analysis,CVDA),并用于动态过程中的早期故障检测,其中证明了基于相异度的统计量优于传统的统计量.此外,针对非线性动态过程微小故障,Pilario和Cao等[1]还提出了一种混合核CVDA的监测方法,增强了检测的灵敏度,并在(continuous stirred tank reactor,CSTR)仿真中验证.为了检测和诊断缓慢变化的故障,Li等[17]提出了结合规范变量残差和T2,Q的组合监测统计量及贡献图,该方法的有效性在工业旋转机械的应用中得到了验证.Chen和Peng等[18]提出了基于累积规范变量相关分析的传感器精度下降检测方法,且在文章中提出了一种阈值学习方法来计算控制限,并在实际的(continuous stirred tank heater,CSTH)过程中通过3种方法比较验证了该方法的可行性和有效性.针对高速自动机多行程运动形态的特点,王宝祥等[19]提出一种分时段基于规范变量残差的动态特性监测方法,在12.7 mm高速自动机上验证了算法的有效性.邓佳伟等[20]针对传统统计局部核主元分析在构造改进残差时未考虑样本的差异性,进而导致故障样本信息易被掩盖的问题提出了基于加权统计局部核主元分析的非线性化工过程微小故障诊断方法.另一类是借助散度对微小变化的灵敏性:Zhang和Diallo等[3]通过测量概率分布之间的微小变化,提出了基于Jensen-Shannon散度的微小故障检测和诊断方案.尹晨[21]针对微小故障容易被系统噪声掩盖的问题,提出了基于K-L(Kullback-Leibler,K-L)散度的微小故障诊断的方法.此外,Ji和Zhou等[22]基于MSPM方法提出了微小故障检测的策略和指标,介绍了两种有代表性的平滑技术,并应用于CSTR过程和一个实际的多单元电动空气制动过程.Shang和Zhou等[23]提出了递归变换分量统计分析用于监测微小故障检测,在田纳西-伊斯曼(Tennessee-Eastman,TE)过程中的仿真应用说明了其检测性能较好.针对闭环系统中的传感器故障,陶松兵等[24]提出了结合K-L距离与快速移动窗口主成分分析的微小故障在线检测与估计模型.Chen等[25]提出了基于K-L散度和PCA的微小故障检测方法,并应用于含未知噪声的电驱动力系统.

通过上述文献分析,与基于规范变量差异分析中统计量D不同.本文不仅考虑了数据的时间序列特性,还利用了残差空间中数据投影的差异性,提出了两个基于规范变量残差的权重平均统计量及其二维贡献图,来实现工业过程中的微小故障检测与诊断.本文的主要贡献如下:首先,提出了两种基于规范变量残差的加权平均统计量;其次,采用CSTR过程中的反应器传感器故障和催化剂失活两个微小故障来验证所提出的检测与诊断方法的可行性和有效性;然后,计算出两个故障中统计量的归一化贡献并绘制二维贡献图以识别故障变量.最后,绘制了两种故障的故障检测率和诊断率对比图.

本文的其余部分概述如下:第2部分回顾了基于CVA的故障检测与诊断方法;提出了两种基于规范变量残差的加权平均统计量和绘制二维贡献图的方法,并给出了离线训练和在线检测与诊断步骤.第3部分将该方法应用在CSTR过程中,画出了各个变量的二维贡献图,并与传统统计量T2,Q以及CVDA方法中D统计量的故障检测率和诊断率进行了比较.第4部分给出了结论.

2 基于规范变量残差的故障监测与诊断

2.1 规范变量分析

无输入的多元状态空间模型可以表述为

其中:x(k)∈Rn和y(k)∈Rm分别是状态和观测向量;k表示采样时刻;A ∈Rn×n,C ∈Rmy×m分别是系统矩阵和输出矩阵;w(k)∈Rn为随机扰动,v(k)∈Rm为测量噪声,假设零均值高斯白噪声.

过去和未来的观测向量yp(r),yf(r)可以通过使用过去和未来的测量数据组成

其中:r表示一般索引,p表示过去观测向量的窗口长度.

根据文献[26]中的方程,过去和未来Hankel矩阵Yp和Yf定义为

给定一个具有m个测量变量的M个观测值的数据集,公式(3)中yp(p+1)的最后一个元素是y1,公式(4)中yf(p+N)的最后一个元素为yl.这些Hankel矩阵中列的维数是N=M −2p+1.随后,协方差矩阵Σpp,Σff和互协方差矩阵Σfp就可以通过下列公式计算出:

CVA试图取得yp(r)和yf(r)之间的最优线性组合.这可以转化成如下一个约束优化的问题:

优化问题可以通过对Hankel矩阵进行奇异值分解来解决

规范变量sr可以通过对公式(12)进行奇异值分解得到

其中J=是转换矩阵,它将过去测量数据转换到规范变量空间.规范变量空间由主元空间和残差空间组成,主要的奇异值可以用来确定主元空间的维度.前n个最大的奇异值可以确定主元空间x(r)∈Rn的维数,剩余的奇异值可以确定残差空间d(r)∈Rmp−n的维数,其中n ＜mp.n的数值应该是已知的.因此,规范变量可以写成如下形式:

状态向量不仅是估计的规范变量的一部分,也是基于过去观测向量的线性组合.因此,状态向量可以被定义为一个基于过去观测向量的线性组合

其中Vn由矩阵V中的前n列构成.预测的误差在过去观测空间中表示如下:

为了检测过程异常,在常规CVA检测方法中考虑了来自状态空间和噪声空间的T2,Q两个检测统计量.T2统计量代表了模型内部的变化,而Q统计量代表了模型外部的变化.

而基于状态空间的T2检测统计量是由文献[11]中提出的

置信度α的相应控制限(UCL)如下推导:

3 基于规范变量残差的统计量和控制限

与仅使用过去矢量数据的传统统计量T2和Q不同,基于规范变量残差的统计量同时利用未来矢量数据和过去矢量数据.在获得残差信号后,基于残差统计量的故障检测包括两个主要步骤:第1步是建立监测统计模型,第2步是根据数据分布确定合适的控制限.两个残差信号可以由下面公式产生:

其中矩阵Σr1,Σr2分别是残差d1,d2的协方差矩阵.受参考文献[18]的启发,提出了加权平均残差统计量.因考虑到连续样本的序列相关性,所以采用加权平均残差统计量来提高故障检测率和诊断率.这两个加权平均残差统计量由下列公式给出:

实际上,由于周围的噪声和有限的采样,过程变量并不严格遵循高斯分布[18].当数据不遵循正态分布时,通常采用概率密度函数估计(probability density function,PDF)的方法.因为核密度估计(kernel density estimation,KDE)是一种估计PDF的非参数方法[11],所以KDE也被考虑在内.假设x1是随机变量,p(x1)表示密度函数

因此,如果知道p(x1),可以通过用WD1和WD2统计量的概率密度函数估计特定置信度α的控制限.两个统计量的控制限可由以下方法估算

如果出现WD1＞的情况,则检测到异常状况或故障.

在检测出微小故障后,需要对产生的故障进行诊断,进而确定工业过程中是哪一个或几个变量发生故障导致的微小故障.依次计算k时刻每个变量对于故障的贡献,贡献最大的变量则会被认为是最有可能的故障原因.本文不仅采用了T2,Q统计量的贡献,也采用了加权平均统计量WD的贡献.第i个变量在k个采样时刻的统计量J,由表示,其中J ∈{T2,Q,WD1,WD2},各统计量的贡献由下列公式进行计算

在计算贡献时,首先利用公式(29)-(31)获取来自正常训练数据的贡献.随后,对于在线操作,相对的贡献可以通过将减去对应的正常训练数据贡献的均值,并且通过正常训练数据贡献的标准差进行缩放计算

将公式(29)-(31)的数据代入公式(32)计算,便可以得到每个变量缩放后的贡献值,更利于绘制贡献图进行比较.

3.1 离线训练和在线检测与诊断程序

该方法包括两个阶段:1)离线训练监测模型;2)在线监测与诊断模型.利用来自正常运行条件的数据进行CVA模型的离线训练.采用包含微小故障的测试数据集来验证所提出的检测与诊断方法.所提出的故障检测与诊断方法的具体步骤总结如下:

离线训练

1) 获取正常运行条件下的训练数据矩阵X;

2) 使用公式(5)-(6)组装矩阵Yp和Yf;

3) 用公式(7)-(9)来计算互协方差矩阵Σpp,Σff和协方差矩阵Σfp;

4) 对公式(12)中的Hankel矩阵H进行奇异值分解;

5) 用公式(15)-(16)计算状态向量和残差向量;

6) 先用公式(21)-(22)计算出两个残差信号后,继续用公式(23)-(24)计算残差统计量,最后利用(25)-(26)求出权重平均残差统计量WD1和WD2;

7) 使用公式(28)确定统计量WD1和WD2的相应控制限.如果其中一个统计量超出了它们相应的控制限,则表明检测到故障发生;

8) 将正常运行的训练数据代入公式(29)-(31),计算出各统计量的贡献,用于在线检测与诊断.

在线检测与诊断

1) 在时间点2p采集一个观测值,并使用训练数据集Xn的均值和方差对样本进行标准化;

2) 使用公式(15)-(16)将其投影到状态向量和残差向量空间;

3) 使用公式(25)-(26)计算新观测值的统计量WD1和WD2.如果一个或两个统计数据超过了KDE确定的相应控制限,则检测到故障;

4) 然后进行故障诊断,将含有故障的测试矩阵代入公式(29)-(31)计算出结果,再代入公式(32)进行缩放计算,最后进行归一化后绘制出贡献图,通过贡献图对故障进行诊断,最终确定故障变量.

4 在连续搅拌反应器系统中的应用

人们可以利用连续搅拌釜反应器(CSTR)系统[12,26]来验证提出的化工过程监测方法是否有效.CSTR模拟系统的原理图如图1所示.

溶剂和反应物a的入口流在一级不可逆反应的条件下产生单一组分B,即出口流.夹套的冷却流用来进行传递反应产生的热量.反应器的液位和温度可以通过串级控制策略调节.根据能量、质量和组分平衡的基本原理,CSTR模拟系统的动力学模型可以表述为

其中:A是反应器的横截面积,Ca表示反应器中物质a的浓度,Cp是内容物的热容量,Caf是进料流中物质a的浓度,Cpc是冷却剂的热容量,h是反应器的液位,E是活化能,k0是预膨胀因子,Qc是冷却剂流速,Qf是进入反应器的进料流速,R是气体常数,Qo是反应器的出口流速,T表示反应器的温度,Tcf是冷却剂进料的温度,Tc是冷却夹套中冷却剂的温度,Tf是进料流的温度,U是传热系数,ΔH是反应热,ρ是内含物密度,Ac是总传热面积,ρc是冷却剂的密度.

CSTR仿真系统在正常操作下以10 s为采样时间产生训练和测试数据集.表1列出了10个测量变量和模拟过程中使用的相应高斯噪声标准偏差值.

表1 测量变量和噪声标准变量差值Table 1 The measured variables and standard deviation of the noise

通过使用训练数据集的300个样本来建立CVA模型.初始故障由两种类型的漂移变化来描述,从k=301到k=700[26-27]:(1)反应器温度传感器T上升的故障T=0.05 K·min−1;(3)催化剂失活的初始故障E/R为3 K·min−1.

4.1 案例研究1

在本案例研究中,考虑了反应器温度传感器T的微小故障,前300 个采样点是正常数据,在k=300到k=700 时引入了反应器温度上升的故障T=0.05 K·min−1.图2显示了常规的检测统计量T2,Q的监测效果,如图所示其故障检测效果较差,故障检测率分别为81.59%和39.39%.图3中所示为基于残差的统计量WD1和WD2的检测效果,分别给出了相对较高的98.98%和72.02%的故障检测率,图3中基于统计量D的故障检测时间较长.通过图2-3中对比最先检测到故障的顺序,表明图3中WD1统计量的故障检测效果最好.

图2 案例1中T2,Q统计量的监测图表Fig.2 The monitoring chart of statistics T2,Q in case 1

图3 案例1中WD1,WD2,D统计量的监测图表Fig.3 The monitoring chart of statistics WD1, WD2, D in case 1

图4中,给出了T2,Q,WD1,WD2,D5种统计量的贡献图与FIR图.如图4中贡献图所示,案例1的故障与变量3,5,10有关,其根本原因是因为反应器温度T的变化.而反应器温度T会影响冷却套中冷却剂温度Tc,进而影响到冷却剂流量Qc,这样最终导致3个变量的贡献较大.

图4 案例1中T2,Q,WD1,WD2,D 5种统计量的贡献图与FIR图Fig.4 Contribution plot of the five statistics T2, Q, WD1,WD2,D and FIR in case 1

为了准确比较统计量的检测和诊断性能,采用了故障检测率和故障诊断率.故障诊断率比较如图9所示,基于WD1的贡献图故障诊断率最高为91.71%,而CVDA中基于D统计量的故障诊断率为91%,稍低于WD1.案例中的故障检测率(FDR)与故障识别率(FIR)的具体公式如下:

其中:Fv是故障样本的数量,Fn是故障样本的总数量,Nv是具有最大贡献的故障变量个数,Nn是从故障开始时刻的样本总数.

4.2 案例研究2

在本案例研究中,从采样时刻k=301到k=700引入催化剂失活E/R为3 K·min−1的缓慢变化微小故障.图5显示了传统CVA方法的检测统计量T2,Q的检测结果,其故障检测率分别为70.33%,58.57%.图6中所示基于残差的统计量WD1和WD2的故障检测率分别为91.82%和60.62%要远高于前者,且WD1和WD2在采样时刻k=299,k=381检测到故障比T2,Q,D检测到故障的时间更早,因此故障检测性能更好.

图5 案例2中T2,Q统计量的监测图表Fig.5 The monitoring chart of statistics T2,Q in case 2

图6 案例2中T2,Q,WD1,WD2,D统计量的监测图表Fig.6 The monitoring chart of statistics WD1, WD2, D in case 2

故障的贡献图表明此次故障与变量5有着最大的关系,由于故障的根本原因是催化剂失活影响反应器温度T变化,进而导致冷却夹套中冷却剂的温度Tc的变化.图7是T2,Q,WD1,WD2,D这5种统计量的贡献图与FIR图,WD1和WD2的故障诊断率分别为77.20%,75.39%.经过图7中故障诊断率的对比,表明所提出的方法要优于传统的T2,Q,D统计量.

图7 案例2中T2,Q,WD1,WD2,D五种统计量的贡献图与FIR图Fig.7 Contribution plot of the five statistics T2, Q, WD1,WD2,D and FIR in case 2

上述所研究的两种故障的检测率和诊断率比较如图8-9所示.在这两种情况下,统计量WD1的所有FDR都远高于统计量T2.统计量WD2的FDR高于统计量Q的FDR.在故障案例1中统计量WD1的FIR与统计量D的FIR基本一致,但在故障案例2中高于统计量D的FIR.仿真结果表明了基于残差的统计量的故障检测与诊断的有效性.

图8 两个案例中统计量T2,Q,WD1,WD2,D的FDR比较Fig.8 FDR of the statistics T2,Q,WD1,WD2,D in the two case studies

图9 两个案例中统计量T2,Q,WD1,WD2,D的FIR比较Fig.9 FIR of the statistics T2, Q, WD1, WD2, D in the two case studies

综上所述,基于规范变量残差的两个加权平均统计量能够有效提高微小故障的故障检测率与诊断率.

5 结论

本文将提出的基于规范变量残差的化工过程微小故障检测和诊断方法应用于CSTR过程进行仿真研究.该方法提出了两种基于规范变量残差的加权平均统计量,并用核密度估计来确定其控制限,最后画出二维贡献图进行故障诊断.在CSTR过程仿真应用中,通过故障检测率和故障诊断率的比较,验证了该方法较好的检测与诊断性能.仿真结果表明,与传统的统计量T2,Q以及CVDA中的D统计量相比,这两种基于规范变量残差的加权平均统计量不仅能够及时检测到微小故障,而且不同程度地提高了故障检测率与诊断率,是一种有效的微小故障诊断方法.下一步工作将研究基于规范变量的微小故障预测.