一种适用于多采样率化工过程的故障检测方法

本发明涉及一种化工过程运行状态监测方法，特别涉及一种适用于多采样率化工过程的故障检测方法。

由于计算机技术与先进测量仪表技术在现代化工工业过程中的广泛应用，化工生产过程中的温度，压力，流量等数据信息可以被实时测量并存储起来，这些海量样本数据为当前化工“大数据”建设提供了坚实的数据基础。近十几年来，利用采样数据来检测化工过程运行中出现的故障在安全化工生产领域受到了越来越多的重视，其思想在于通过采样数据来实时检测化工过程运行是否出现故障，从而保证安全生产与维持产品质量的稳定性。经过十几年的发展历程，适用于化工过程对象的数据驱动故障检测方法已经从早期的单变量思路转变成以主成分分析(Principal Component Analysis，缩写：PCA)与偏最小二乘(Partial Least Square，缩写：PLS)等多变量分析算法为基础的多变量故障检测方法。这些主流的故障检测方法实施的核心主要关注于数据潜在特征的挖掘。换句话说，所建立的数据驱动模型都是旨在提取采样数据中潜藏的特征。

一个典型的化工过程对象的核心生产单元都包含有：反应设备、分离设备、冷却设备、和回流装置。这些设备之间相互链接，运行机理与控制系统回路较为复杂。此外，由于测量仪器仪表购置批次不同或对测量物理量的实时性要求不同，导致不同测量变量会存在不同的采样率。举个简单的例子：化工过程中的温度、压力等物理量可在较短的时间内得到测量，其测量周期可按照分种或秒来测量；然而，产品的浓度等检验指标，因测量仪表或人员的限制，其测量周期会相对较长，通常是按照小时来录入。如此一来，在相同采样时间段内，化工过程各个测量仪器采集到的样本数据的个数是不一致的。

一方面，从检测故障发生与否的角度出发，任何测量变量的异常变化都能够反映出化工过程运行中出现了故障。因此，为了保证故障检测的全面性，所有的测量数据都应该得到有效的应用。另一方面，若是按照不同采样率将测量变量分成多个变量块或者集合，然后对每种采样率数据分开进行建模，会忽略变量块或集合之间的相关性。这主要是因为化工过程的运行是一个整体，从运行机理上讲各个测量变量间是存在耦合关系的，反应在数据上就是各测量变量间的相关性了。若是变量的相关性发生了变化，化工过程同样是进入了故障运行状态的。因此，研究发明出一种适用于多采样率化工过程的故障检测方法是一个还未得到有效解决的问题。

本发明所要解决的主要技术问题是：发明出一种适用于多采样率化工过程的故障检测方法。由于需要考虑多采样率导致测量样本数据不一致问题，本发明方法解决的技术难点则在于：如何在充分考虑测量变量间相关性的基础上实施多采样率的样本数据的特征提取。

本发明方法解决上述问题所采用的技术方案为：一种适用于多采样率化工过程的故障检测方法，包括以下步骤：

步骤(1)：确定化工过程对象中各测量仪器仪表的采样频率，并将采样频率相同的测量仪器仪表所对应的测量变量归属为一个变量集合，得到B个变量集合φ1，φ2，…，φB，并将各变量集合所对应的采样频率分别记录为f1，f2，…，fB。

步骤(2)：利用化工过程对象中的测量仪器仪表采集该化工过程正常运行状态下的样本数据，并根据步骤(1)中的B个变量集合φ1，φ2，…，φB对应地将采样数据分别组成B个数据矩阵X1，X2，…，XB，其中，由第b个变量集合对应的Nb个样本数据组成、Mb为第b个变量集合中的测量变量个数、表示Mb×Nb维的实数矩阵、下标号b＝1，2，…，B。

步骤(3)：分别对B个数据矩阵X1，X2，…，XB进行标准化处理后，对应得到新矩阵再根据公式Tb＝1/fb分别计算出各变量集合所对应的采样周期T1，T2，…，TB，并出T1，T2，…，TB的最小公倍数T0后，再根据公式nb＝T0/Tb计算出参数n1，n2，…，nB。

由于新矩阵之间列向量的个数不一致，需要进行适当的变换处理。若将新矩阵中Nb个列向量表示成后，可按照如下所示公式①分别构造延时矩阵Z1，Z2，…，ZB：

其中，表示新矩阵中的第i列向量，表示Mb×1维的列向量，下标号i＝1，2，…，Nb。

上述步骤(1)至步骤(3)涉及根据不同采样率区分出的多个变量集合，这些变量集合之间是存在相关性的。因此，为了充分考虑不同采样率下变量块间的相关性，在实施特征提取时，本发明方法涉及的多采样率建模算法需满足如下所示目标函数：

上式②中，为投影向量，为系数向量，表示kronecker乘积，下标号c＝1，2，…，B。若设则的运算结果如下所示：

其中，为系数向量βb中的各个元素。

上式②的求解需要用到经典的拉格朗日乘子法，即通过引入拉格朗日乘子λb与γb后构造如下所示函数L：

函数L分别对wb与βb求解偏导数，并将偏导数置零后，可得到如下所示关系等式：

从上式⑤与⑥中可以发现，这是个特征值与特征向量的求解问题。然而，由于wb的求解需要已知βb，而βb的求解需要已知wb，因此上式的特征值问题求解相互耦合。为此，本发明方法发明出如下所示的循环迭代步骤。

步骤(4)：利用多采样率建模算法求解出投影矩阵W1，W2，…，WB、载荷矩阵P1，P2，…，PB、以及得分矩阵S1，S2，…，SB，具体的实施过程如下所示。

步骤(4.1)：设置d＝1后，初始化投影矩阵W1，W2，…，WB、载荷矩阵P1，P2，…，PB、以及得分矩阵S1，S2，…，SB为空矩阵。

步骤(4.2)：将新矩阵中Nb个列向量表示成后，按照如下所示公式⑦分别构造延时矩阵Z1，Z2，…，ZB：

其中，表示新矩阵中的第i列向量，表示Mb×1维的列向量，下标号i＝1，2，…，Nb。

步骤(4.3)：初始化w1，w2，…，wB为任意Mb×1维的实数列向量，并初始化β1，β2，…，βB为任意nb×1维的列向量。

步骤(4.4)：根据公式wb＝wb/||wb||与公式βb＝βb/||βb||分别对wb与βb实施单位化处理，其中||wb||表示计算wb的长度，||βb||表示计算βb的长度。

步骤(4.5)：根据公式如下所示公式⑧计算矩阵Θb：

上式⑧中，表示Mb×Mb维的单位矩阵；当b＝c时，Hbc＝0；当b≠c时，Hbc＝1；上标号T表示矩阵或向量的转置符号，下标号c＝1，2，…，B。

步骤(4.6)：判断是否满足条件：d＞1；若否，则求解特征值问题中最大特征值λb所对应的特征向量若是，则求解特征值问题中最大特征值λb所对应的特征向量

步骤(4.7)：根据公式对特征向量实施单位化处理从而得到w1，w2，…，wB，再根据如下所示公式⑨计算矩阵Φb

上式⑨中，表示nb×nb维的单位矩阵。

步骤(4.8)：分别计算矩阵Φ1，Φ2，…，ΦB的最大特征值所对应的特征向量ξ1，ξ2，…，ξB，再分别根据公式βb＝ξb/||ξb||对特征向量实施单位化处理得到β1，β2，…，βB。

步骤(4.9)：判断w1，w2，…，wB是否收敛，收敛的标准为元素不再发生变化；若否，则返回步骤(4.5)；若是，则执行下一步骤(4.10)。

步骤(4.10)：根据公式分别计算出得分向量s1，s2，…，sB后，再分别根据公式计算出载荷向量p1，p2，…，pB。

步骤(4.11)：分别将投影矩阵W1，W2，…，WB中的第d列的列向量设置成w1，w2，…，wB，分别将载荷矩阵P1，P2，…，PB中的第d列的列向量设置成p1，p2，…，pB，再分别将得分矩阵S1，S2，…，SB中的第d行的行向量设置成s1，s2，…，sB。

步骤(4.12)：根据公式分别更新后，判断是否满足条件：d≥min{M1，M2，…，MB}或min{λ1，λ2，…，λB}≤ζ；若否，则设置d＝d+1后返回步骤(4.2)；若是，则得到最终的B个投影矩阵W1，W2，…，WB、B个载荷矩阵P1，P2，…，PB、以及B个得分矩阵S1，S2，…，SB后，再根据公式θb＝W(PTW)-1计算特征提取矩阵θ1，θ2，…，θB；其中，min{M1，M2，…，MB}表示M1，M2，…，MB中的最小值，min{λ1，λ2，…，λB}表示λ1，λ2，…，λB中的最小值，阀值ζ的取值范围为：B/15≤ζ≤B/5。

步骤(5)：判断在b分别等于1，2，…，B时，是否满足条件：Mb＝d；若否，则对步骤(4.11)中更新的矩阵实施奇异值分解：并设置分解矩阵其中对角矩阵Ab对角线上的元素为的非零奇异值，Ub与Vb是奇异值分解得到的两个酉矩阵；若是，则不对矩阵实施任何操作。

步骤(6)：判断在b分别等于1，2，…，B时，是否满足条件：Mb＝d；若否，则根据公式与分别计算监测指标向量Db与Qb后，利用核密度估计法确定出控制上限Db，lim与Qb，lim；若是，则只需根据公式计算监测指标向量Db，并利用核密度估计法确定出控制上限Db，lim；其中diag{}表示将大括号内矩阵对角线元素转变成列向量的操作，协方差矩阵

离线建模过程至此已全部结束，接下来就是利用在线采样数据不间断地实施对多采样率化工过程运行状态的实时监测。

步骤(7)：根据B个不同的采样频率f1，f2，…，fB在线采集多采样率化工过程对象的样本数据x1，x2，…，xB，其中，为第b个变量集合φb所对应的样本数据。

由于存在B个不同的采样频率f1，f2，…，fB，样本数据x1，x2，…，xB被采集到的频率也不相同，只要x1，x2，…，xB中任何一个样本数据被采集到，即刻开始实施步骤(8)。

步骤(8)：对xb实施与步骤(3)中相同的标准化处理得到向量后，根据公式tb＝θbTxb与分别计算特征向量tb与残差向量eb。

步骤(9)：判断是否满足条件：Mb＝d；若否，则根据公式与计算第b个变量集合的监测指标与若是，则只需根据公式计算第b个变量集合的监测指标

步骤(10)：判断是否满足条件：且若是，则当前采样时刻该化工该过程运行状态正常，返回步骤(7)继续实施对下一采样时刻的运行状态监测；若否，则执行步骤(11)从而决策是否触发故障警报。

步骤(11)：返回步骤(7)继续实施对下一采样时刻样本数据的运行状态监测，若连续3个采样时刻的监测指标都不满足步骤(10)中的判断条件，则触发故障警报；否则，不触发故障警报，并返回步骤(7)继续实施对下一采样时刻的运行状态监测。

与本发明方法的优势与特点如下所示。

首先，本发明方法不仅解决了多采样率化工过程的数据建模问题，而且多采样率建模算法充分地考虑了不同采样频率测量变量之间的相关性。其次，本发明方法将不同采样频率测量变量间的相关性特征提取出来后，单独进行监测从而可以起到对故障的发生进行预测监测的效果的目的；最后，通过具体的实施案例，对比验证了本发明方法在监测多采样率化工过程运行状态上的可行性。

图1为本发明方法的实施流程示意图。

图2为TE化工过程的组成结构示意图。

图3为本发明方法在监测多采样率的TE化工过程时的在线监测曲线图。

下面结合附图和具体实施方式对本发明进行详细说明。

如图1所示，本发明公开了一种适用于多采样率化工过程的故障检测方法，下面结合一个具体应用实例来说明本发明方法的具体实施方式。

如图2所示，应用对象是来自于美国田纳西-伊斯曼(TE)化工过程对象，是该化工厂化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于化工过程状态监测研究。整个TE过程包括22个测量变量、12个操作变量、和19个产品成分测量变量。该TE化工过程对象涉及三种不同的采样频率：3分钟每次、9分钟每次、和15分钟每次，是一个非常典型的多采样率化工过程。

步骤(1)：确定TE化工过程对象中各测量仪器仪表的采样频率，并将采样频率相同的测量仪器仪表所对应的测量变量归属为一个变量集合，得到如表1所示的B＝3个变量集合φ1，φ2，φ3，并将各变量集合所对应的采样频率分别记录为f1＝1/3，f2＝1/9，f3＝1/15。

表1：不同采样频率所对应的变量详情。

步骤(2)：利用化工过程对象中的测量仪器仪表采集该化工过程正常运行状态下的样本数据，并根据步骤(1)中的B＝3个变量集合φ1，φ2，φ3对应地将采样数据分别组成3个数据矩阵X1，X2，X3。

步骤(3)：分别对B＝3个数据矩阵X1，X2，X3进行标准化处理后，对应得到新矩阵再根据公式Tb＝1/fb分别计算出各变量集合所对应的采样周期T1＝3，T2＝9，T3＝15，并出T1，T2，T3的最小公倍数T0＝45后，再根据公式nb＝T0/Tb计算出参数n1＝15，n2＝5，nB＝2。

步骤(4)：利用多采样率建模算法求解出投影矩阵W1，W2，W3、载荷矩阵P1，P2，P3、得分矩阵S1，S2，S3、以及特征提取矩阵θ1，θ2，θ3，具体的实施过程详见前述步骤(4.1)至步骤(4.12)。

步骤(5)：判断在b分别等于1，2，3时，是否满足条件：Mb＝d；若否，则对步骤(4.11)中更新的矩阵实施奇异值分解：并设置分解矩阵若是，则不对矩阵实施任何操作。

步骤(6)：判断在b分别等于1，2，…，B时，是否满足条件：Mb＝d；若否，则根据公式与分别计算监测指标向量Db与Qb后，利用核密度估计法确定出控制上限Db，lim与Qb，lim；若是，则只需根据公式计算监测指标向量Db，并利用核密度估计法确定出控制上限Db，lim。

离线建模阶段至此完成，接下来进入在线异常状态检测阶段，需要实时利用TE化工过程对象的在线采样数据。

刚开始实施在线状态监测时，TE化工过程对象运行在正常状态下，一段时间后，冷凝器冷却水的温度发生了阶跃跳变的故障。

步骤(7)：根据3个不同的采样频率，在线采集多采样率TE化工过程对象的样本数据x1，x2，x3。

由于存在3个不同的采样频率，样本数据x1，x2，x3被采集到的频率也不相同，只要x1，x2，x3中任何一个样本数据被采集到，即刻开始实施步骤(8)。

步骤(8)：对xb实施与步骤(3)中相同的标准化处理得到向量后，根据公式tb＝θbTxb与分别计算特征向量tb与残差向量eb。

步骤(9)：判断是否满足条件：Mb＝d；若否，则根据公式与计算第b个变量集合的监测指标与若是，则只需根据公式计算第b个变量集合的监测指标

步骤(10)：判断是否满足条件：且若是，则当前采样时刻该化工该过程运行状态正常，返回步骤(7)继续实施对下一采样时刻的运行状态监测；若否，则执行步骤(11)从而决策是否触发故障警报。

步骤(11)：返回步骤(7)继续实施对下一采样时刻样本数据的运行状态监测，若连续3个采样时刻的监测指标都不满足步骤(10)中的判断条件，则触发故障警报；否则，不触发故障警报，并返回步骤(7)继续实施对下一采样时刻的运行状态监测。

如图3所示，本发明方法在监测该多采样率化工过程运行状态的实施监测曲线图。此外，利用传统PCA方法对TE过程实施监测，对于不同采样频率采取复制数据的手段从而保证PCA方法的可适用性。从图3的对比中可以发现，本发明方法能持续不断的出发警报，效果优越于传统PCA方法。

上述实施案例只用来解释说明本发明的具体实施，而不是对本发明进行限制。在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改，都落入本发明的保护范围。