一种基于简略核主元分析的精馏塔异常状态识别方法

本发明涉及一种化工过程运行状态监测方法，特别涉及一种基于简略核主元分析的精馏塔异常状态识别方法。

由于计算机技术与先进测量仪表技术广泛应用于化工生产，化工生产过程中的温度，压力，流量等数据信息可以被实时测量并存储起来，这些海量样本数据为当前化工智能制造提供了坚实的数据基础。近十几年来，利用采样数据来检测化工过程运行中出现的异常状态在安全化工生产领域受到了越来越多的重视。经过十几年的发展，化工过程异常状态监测已经建立起了一套以主成分分析(Principal Component Analysis，缩写：PCA)与偏最小二乘(Partial Least Square，缩写：PLS)等多变量分析算法为基础的异常状态识别方法。这些数据驱动的方法技术的核心本质在于对采样数据进行潜在特征的挖掘。换句话说，所建立的数据驱动模型都是旨在提取采样数据数据中潜藏的特征。

精馏塔设备是石油化工生产中应用极为广泛的一种传质传热装置，其主要作用通过利用混合物中各组分具有不同的挥发度，实现对物质的分离目的。因此，实时检测精馏塔设备是否出现异常运行状态对于整个化工生产具有重要的意义。精馏塔中各层塔板的温度，回流流量等数据信息都可以反应出精馏塔设备是否工作在期望工况上。考虑到精馏塔设备采样数据的非线性变化特征，可使用基于核PCA(Kernel PCA，缩写：KPCA)模型的异常状态识别方法。KPCA方法通过一个隐性的非线性映射函数将原始数据投影至高维空间后，使得非线性特征在高维空间可直接按照PCA方法的思路进行特征分析与提取。虽然KPCA方法能很好的处理精馏塔采样数据的非线性特征，并通过监测非线性特征的变化来实现异常状态识别目的，但是KPCA方法技术本身的缺陷同样是非常明显的。

数据驱动的方法技术有一个共同特点，就是训练数据越多越好。换句话说，KPCA方法在训练模型时理应使用尽可能多的样本数据。然而，KPCA的离线建模阶段完成后，对新采样数据实施特征变换时，其在线的计算量是和离线建模阶段所用数据量直接成正比例关系的。此外，在线为各个样本数据计算核向量时，都是涉及到了非线性的指数运算。若是采样间隔短，利用KPCA实施精馏塔异常状态的在线识别就会面临着严峻的实时性效率问题。换句话说，KPCA方法正在利用当前采样时刻进行运算与异常状态检测时，下一个新采样得到的样本数据已经出现了。因此，如何提升KPCA方法的在线计算效率是一个很值得关注的问题，尤其是在采样间隔很短的精馏塔设备异常状态识别上。

本发明所要解决的主要技术问题是：如何提升利用KPCA方法实时检测精馏塔设备异常状态的效率。本发明方法是通过降低离线建模时的数据量来实现对在线识别效率的提升。简单的来讲，与其使用尽可能多的样本数据训练KPCA模型，不如使用典型具有代表性的样本数据来训练KPCA模型。

本发明方法解决上述问题所采用的技术方案为：一种基于简略核主元分析的精馏塔异常状态识别方法，包括以下步骤：

步骤(1)：利用精馏塔设备所安装的测量仪表，在精馏塔正常运行状态时采集N个样本数据x1，x2，…，xN，其中第i个采样时刻的样本数据xt∈Rm×1由m个采样数据组成，具体包括：塔釜液位、塔釜压力，塔釜底部产品流量，进料流量，进料温度，顶部回流流量，冷凝器液位，和各层塔板的温度，i∈{1，2，…，N}。

步骤(2)：对N个样本数据x1，x2，…，xN实施标准化处理，得到N个m×1维的数据向量

步骤(3)：对N个数据向量进行边缘点的筛选，从而保留n个数据向量z1，z2，…，zn建立核主成分分析模型，其中n小于N，具体的实施过程如下所示。

步骤(3.1)：初始化i＝1。

步骤(3.2)：按照欧式距离的大小，从N个数据向量中为搜寻出k个欧式距离最小的数据向量，并记录成其中i∈{1，2，…，n}。

步骤(3.3)：根据如下所示公式计算的法向量fi。

其中c＝1，2，…，k，表示计算与之间的欧式距离。

步骤(3.4)：根据如下所示公式计算边缘点指标ζi：

上式中，θci为二进制数，其取值规律如下所示：

其中，上标号T表示矩阵或向量的转置。

步骤(3.5)：判断是否满足条件i＜N；若是，则设置i＝i+1后，返回步骤(3.2)；若否，则按照数值大小对边缘点指标ζ1，ζ2，…，ζN进行降序排列，并将最大的n个边缘点指标所对应的数据向量依次记录为z1，z2，…，zn。

步骤(4)：利用z1，z2，…，zn建立核主元分析模型，具体的实施过程如下所示。

步骤(4.1)：根据如下所示公式计算核矩阵K∈Rn×n中第a行第b列的元素K(a，b)：

其中，δ为核参数，通常可取δ＝100m，a∈{1，2，…，n}，b∈{1，2，…，n}，Rn×n表示n×n维的实数矩阵。

步骤(4.2)：根据如下所示公式对核矩阵K进行中心化处理得到矩阵

其中，矩阵IIn∈Rn×n中全部元素都是1。

步骤(4.3)：求解矩阵所有特征值所对应的特征向量，并按照大小对特征值进行降序排列得到λ1≥λ2≥…≥λn，与特征值λ1，λ2，…，λn分别对应的特征向量为p1，p2，…，pn。

步骤(4.4)：计算特征值的均值并确定出特征值λ1，λ2，…，λn中大于均值的个数，记作参数d。

步骤(4.5)：根据如下所示公式计算核矩阵J∈RN×n中第i行第b列元素J(i，b)：

其中，i∈{1，2，…，n}，b∈{1，2，…，n}。

步骤(4.6)：根据如下所示公式对J实施中心化处理得到

上式中，矩阵IIN∈RN×n中所有元素都是1，RN×n表示N×n维的实数矩阵。

步骤(4.7)：根据公式与分别计算主得分矩阵S1与副得分矩阵S2，再根据公式A＝(S1TS1)-1计算矩阵A，其中，主变换矩阵P1＝[p1，p2，…，pd]，副变换矩阵P2＝[pd+1，pd+2，…，pn]。

步骤(4.8)：根据公式D＝diag{S1AS1T}与Q＝diag{S2S2T}分别计算指标向量D与Q，其中，diag{}表示将大括号内的矩阵对角线元素转变成向量的操作运算。

步骤(4.9)：利用核密度估计法估计出指标向量D与Q的控制上限Dlim与Qlim。

从步骤(4.5)至步骤(4.8)中可以看出，虽然使用了边缘点指标筛选出了部分数据向量用于求解KPCA的特征值向量，但是计算指标向量时，用的却是全部的采样数据；一方面，在线计算核向量的样本数据大量减少，这也是为何本发明方法叫做简略核主元分析；另一方面，所有采集的样本数据在确定控制上限时都派上用场了。

本发明方法的离线建模过程至此已全部结束，接下来就是利用在线采样数据不间断地实施对多精馏塔设备运行状态的实时监测，从而及时识别出异常状态。

步骤(5)：在最新采样时刻t，利用精馏塔设备所安装的测量仪表测量得到由m个采样数据组成的数据向量xt∈Rm×1，并对其实施与步骤(2)中相同的标准化处理，得到新数据向量

步骤(6)：根据如下所示公式计算核向量kt∈R1×n中的第b个元素kt(b)：

上式中，b∈{1，2，…，n}，R1×n表示1×n维的实数向量。

从上式⑧中可以看出，本发明方法针对在线测量到的每个数据，只需计算n次即可得到核向量；相比之下，传统的KPCA方法需要计算N次才能得到核向量；因此本发明方法的在线计算量得到了明显的降低。

步骤(7)：根据如下所示公式对核向量kt实施中心化处理得到

上式中，向量IIt∈R1×n中所有元素都为1，矩阵IIn∈Rn×n中全部元素都是1。

步骤(8)：根据公式与分别计算主得分向量s1与副得分向量s2后，再分别根据公式Dt＝s1As1T与Qt＝s2s2T计算识别指标Dt与Qt。

步骤(9)：判断是否满足条件：Dt≤Dlim且Qt≤Qlim；若是，则当前采样时刻精馏塔运行状态正常，返回步骤(5)继续实施对下一最新采样时刻的异常状态识别；若否，则执行步骤(10)从而决策是否识别出异常状态。

步骤(10)：返回步骤(5)继续实施对下一最新采样时刻样本数据的异常状态识别，若连续3个采样时刻的检测指标都不满足步骤(9)中的判断条件，则精馏塔进入异常工作状态；否则，返回步骤(5)继续实施对下一最新采样时刻的异常状态识别。

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

首先，本发明方法旨在提高KPCA用于精馏塔在线异常状态识别时的计算效率，通过边缘点指标筛选出部分具有代表性的数据训练KPCA模型，从而极大的降低了核向量的在线计算量；其次，本发明方法虽有筛选部分数据，但是并未全盘否定未选中的数据，而是通过构造识别指标的方式将其考虑进在线异常状态识别中；最后，通过具体的实施案例，验证本发明方法在计算效率上的提升，与此同时对于异常状态的识别能力却未曾降低。

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

图2为精馏塔设备的实景图和组成结构示意图。

图3为本发明方法在识别多种精馏塔异常状态的能力对比。

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

如图1所示，本发明公开了一种基于简略核主元分析的精馏塔异常状态识别方法，下面结合一个具体应用实例来说明本发明方法的具体实施方式。

从图2中的精馏塔实景图可以看出，精馏塔设备不单只是单独的一个精馏设备，在底部还有配套的再沸器，顶部配套有冷凝器。从图2中的结构示意图中可以看出，该精馏塔设备测量仪表包括：流量仪表，温度仪表，液位仪表三类，对应测量的变量有17个，具体包括：塔釜液位、塔釜压力，塔釜底部产品流量，进料流量，进料温度，顶部回流流量，冷凝器液位，以及每层塔板的温度(总计有10层)。

步骤(1)：利用精馏塔设备所安装的测量仪表，在精馏塔正常运行状态时采集N＝1000个样本数据x1，x2，…，x1000。

步骤(2)：对1000个样本数据x1，x2，…，x1000实施标准化处理，对应得到1000个17×1维的数据向量

步骤(3)：按照上述步骤(3.1)至步骤(3.5)对这1000个数据向量进行边缘点的筛选，从而保留n＝400个数据向量z1，z2，…，z400建立核主成分分析模型，样本个数得到了大幅的降低。

由于样本个数的大幅度降低，在线计算核向量时(公式⑧)所需的计算时间成比例减小。因为原来的KPCA需要针对每个在线测量样本调用公式⑧1000次，而本发明方法只需调用400次即可。很显然，本发明方法的计算效率高。

步骤(4)：利用z1，z2，…，z400建立核主元分析模型，保留主副变换矩阵P1与P2，矩阵A，以及控制上限Dlim与Qlim。

离线建模阶段完成后，即可按照如下所示步骤不间断的对精馏塔实施在线状态监测，需要利用每个新采样时刻的样本数据。

步骤(5)：在最新采样时刻t，利用精馏塔设备所安装的测量仪表测量得到由m个采样数据组成的数据向量xt∈R17×1，并对其实施与步骤(2)中相同的标准化处理，得到新数据向量

步骤(6)：根据上述公式⑧计算核向量kt∈R1×400中的第i个元素kt(i)。

步骤(7)：根据上述公式⑨对核向量kt实施中心化处理得到

步骤(8)：根据公式与分别计算主得分向量s1与副得分向量s2后，再分别根据公式Dt＝s1As1T与Qt＝s2s2T计算识别指标Dt与Qt。

步骤(9)：判断是否满足条件：Dt≤Dlim且Qt≤Qlim；若是，则当前采样时刻精馏塔运行状态正常，返回步骤(5)继续实施对下一最新采样时刻的异常状态识别；若否，则执行步骤(10)从而决策是否识别出异常状态。

步骤(10)：返回步骤(5)继续实施对下一最新采样时刻样本数据的异常状态识别，若连续3个采样时刻的检测指标都不满足步骤(9)中的判断条件，则精馏塔进入异常工作状态；否则，返回步骤(5)继续实施对下一最新采样时刻的异常状态识别。

最后，为了验证本发明方法在提升计算效率的同时，不会对状态监测结果造成负面影响，针对精馏塔在冷凝水温度异常、再沸器温度异常、和回流阀门粘滞的异常状态识别问题，对比分析了本发明方法与传统KPCA的识别几率。从图3中可以看出，本发明方法与传统KPCA方法几乎差异甚微。但是，前面已经分析过本发明方法的计算效率由原来的需要计算1000次降低到现在计算400次，计算耗时能够很明显的降低。