CN110009020B - 一种基于多核主元分析模型的非线性过程监测方法 - Google Patents

Abstract

本发明公开一种基于多核主元分析模型的非线性过程监测方法，旨在建立并融合多个核函数所对应的KPCA模型，从而避免核函数的选择问题，并在此基础上实施有效地非线性过程监测。首先，本发明方法将常用的核函数类型全部考虑进来，避免了核函数的选择问题。因此本发明方法的通用性较强。其次，本发明方法因使用多个核函数分别建立多个不同的非线性过程监测模型，充分发挥了多模型建模的优势。可以说，本发明方法的故障检测效果不会弱于任何一个使用单个核函数的过程监测模型。综合这两点优势，本发明方法克服了传统基于KPCA的过程监测方法的不足，是一种更为优选的非线性过程监测方法。

Description

一种基于多核主元分析模型的非线性过程监测方法

技术领域

本发明涉及一种工业过程监测方法，尤其是涉及一种基于多核主元分析模型的非线性过程监测方法。

背景技术

现代工业过程对象规模的大型化与生产的高效化对实时监测过程运行状态提出了越来越高的要求，及时地检测出过程对象运行过程中出现的故障工况是保证产品质量的根本途径。可以说，针对以故障检测为核心任务的过程监测技术的研究一直伴随着工业发展的历程。而今，数据驱动的过程监测方法是最为主流的实施技术手段，这主要得益于测量仪表与计算机技术的广泛应用，可测量与可存储的工业过程采样数据为数据驱动的过程监测方法奠定了坚实的数据基础。考虑到现代工业过程对象的非线性特性，使用非线性的数据建模与监测方法通过能保证其有效性。

作为最为广泛使用的非线性过程建模与监测方法，核主元分析(KernelPrincipal Component Analysis，KPCA)算法通过利用核学习技巧实现了对采样数据非线性特征的描述。核学习技巧的巧妙之处在于通过一个未知的非线性函数将原采样数据映射至高维空间，映射后的数据在高维空间中只存在线性关系特征，但是KPCA却不需要知道非线性函数的具体形式。正是得益于这种优势，在现有科研文献与专利技术材料中，基于KPCA算法或者类似核学习思想的建模与监测方法层出不穷。然而，KPCA算法虽然不需要知道非线性函数的具体形式，但是需要人为预先给定一个核函数，而常用核函数的选择有多种可选形式。一般来讲，基于KPCA算法的过程监测模型都是指定其中一种通用性较好的核函数(即高斯核函数)，鲜有考虑使用其他两种核函数的研究成果。可是，核函数直接决定着非线性映射函数的映射结果，对基于KPCA的非线性过程监测性能有着直接影响。

此外，另一个值得指出的问题是，KPCA算法如何选择合适的核函数至今仍未出现一个被广泛接受的标准。这主要是因为KPCA算法所建立的模型是一种单分类模型，建模过程只需要利用正常工况下的采样数据即可。在没有历史故障数据的前提下，无法从提升故障检测效果的角度出发选择合适于被监测对象的核函数。即使某些过程对象历史数据库中可以提供些许故障工况下的采样数据，但是故障工况的类型不可计数，利用这些有限的故障数据来帮助选择核函数所得到的KPCA模型有可能只针对这些故障类型，而对提升未知类型故障的检测可能是没有任何作用的。因此，基于KPCA的非线性过程监测方法中核函数的选择是个值得关注的问题，但是至今却未得到深入研究。

近十几年来，多模型的建模思路已被较多领域验证为一种能有效提升模型精度的策略，多模型的泛华能力通常也优越于单个模型。多模型建模实施的关键在于如何对多个子模型进行有效融合，这种多模型的思路或许能为KPCA算法中核函数的选择问题提供可借鉴的方案。也就是说，放弃选择单个的核函数，转而建立多个核函数的非线性模型，这样可以避免核函数的选择问题，同时又能发挥出多模型的优势。

发明内容

本发明所要解决的主要技术问题是：如何建立并融合多个核函数所对应的KPCA模型，从而避免核函数的选择问题，并在此基础上实施有效地非线性过程监测。本发明解决上述技术问题所采用的技术方案为：一种基于多核主元分析模型的非线性过程监测方法，包括以下步骤：

步骤(1)：在生产过程正常运行状态下，采集n个样本数据组成训练数据矩阵X∈R^{n ×m}，并对X中各列实施标准化处理，得到标准化后的矩阵

其中，

表示标准化后的第i个样本数据，下标号i＝1，2，…，n，m为测量变量个数，上标号T表示矩阵或向量的转置，R为实数集，R^n×m表示n×m维的实数矩阵。

步骤(2)：设置高斯核函数的参数g后，按照下式①计算高斯核矩阵K_G∈R^n×n：

上式中，K_G(i，j)表示高斯核矩阵K_G中的第i行、第j列元素，下标号i＝1，2，…，n与j＝1，2，…，n，exp表示以自然常数e为底的指数函数，符号|| ||表示计算向量的长度，参数g的建议值为g＝10m。

步骤(3)：设置多项式核函数的参数p后，根据下式②计算多项式核矩阵K_P∈R^n×n：

上式中，K_P(i，j)表示多项式核矩阵K_P中的第i行、第j列元素，参数p的建议值为p＝3。

步骤(4)：根据如下所示公式③计算Sigmoid核矩阵K_S∈R^n×n：

上式中，K_S(i，j)表示Sigmoid核矩阵K_S中的第i行、第j列元素，tanh表示双曲正切函数。

步骤(5)：根据如下所示公式④计算线性核矩阵K_L∈R^n×n：

上式中，K_L(i，j)表示线性核矩阵K_L中的第i行、第j列元素。

步骤(6)：依据如下公式分别对核矩阵K_G、K_P、K_S、和K_L实施中心化处理，对应得到中心化后的核矩阵

和

上式中，下标号c∈{G，P，S，L}分别表示高斯核、多项式核、sigmoid核、以及线性核，方阵Θ∈R^n×n中各元素都等于1。

步骤(7)：计算高斯核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_G为非零特征值的个数。

步骤(8)：计算多项式核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_P为非零特征值的个数。

步骤(9)：计算Sigmoid核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_S为非零特征值的个数。

步骤(10)：计算线性核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_L为非零特征值的个数。

步骤(11)：保留高斯核主元变换矩阵

多项式核主元变换矩阵

Sigmoid核主元变换矩阵

以及线性核主元变换矩阵

步骤(12)：设置核主元个数为d后，先根据公式

计算核主元得分矩阵Φ_G、Φ_P、Φ_S、和Φ_L，其中A_c(1：d)表示矩阵A_c中的第1列至第d列向量组成的矩阵，再根据公式Λ_c＝Φ_c ^TΦ_c/(n-1)计算各核主元得分矩阵的协方差矩阵Λ_G、Λ_P、Λ_S、和Λ_L。

步骤(13)：根据如下所示公式分别计算核主元监测统计量的控制上限D_lim和残差监测统计量的控制上限Q_G，lim、Q_P，lim、Q_S，lim、以及Q_L，lim：

上式中，g_c与h_c分别为卡方分布的权重系数与自由度，F_d，n-d，σ表示自由度为d与n-d的F分布在置信度σ＝99％下的取值，

表示自由度为h_c的卡方分布在置信度σ＝99％下的取值。

以上所述步骤为本发明方法的离线建模阶段，其中步骤(3)至步骤(6)中分别使用了四种典型的核函数计算和矩阵。本发明方法的在线过程监测阶段包括如下所示步骤(14)至步骤(18)。

步骤(14)：在线采集最新采样时刻的样本数据x_new∈R^1×m，并对其进行标准化处理得到

步骤(15)：根据如下所示公式分别计算高斯核向量z_G∈R^1×n、多项式核向量z_P∈R^{1 ×n}、Sigmoid核向量z_S∈R^1×n、以及线性核向量z_L∈R^1×n：

上式中，z_G(i)、z_P(i)、z_S(i)、和z_L(i)分别为核向量z_G、z_P、z_S、和z_L中的i个元素。

步骤(16)：依据下式计算对各核向量实施中心化处理得到

和

其中，行向量θ＝[1，1，…，1]∈R^1×n中各元素都等于1。

步骤(17)：根据公式

计算核主元得分向量t_G、t_P、t_S、和t_L后，根据如下所示公式分别计算监测统计量D_c与Q_c：

上式中，t_c(1：d)表示向量t_c中第1列至第d列元素组成的向量，t_c(d+1：N_c)表示向量t_c中第d+1列至第N_c列元素组成的向量，下标号c∈{G，P，S，L}。

步骤(18)：对于下标号c的所有四个取值，判断是否满足条件：D_c≤D_lim且Q_c≤Q_c，lim？若否，则当前采样时刻已进入故障工况；若是，则过程对象处于正常运行状态，返回步骤(14)实施对下一采样时刻的过程监测。

与现有方法相比，本发明方法的优势在于：

首先，本发明方法将常用的核函数类型全部考虑进来，避免了核函数的选择问题。因此本发明方法的通用性较强。其次，本发明方法因使用多个核函数分别建立多个不同的非线性过程监测模型，充分发挥了多模型建模的优势。可以说，本发明方法的故障检测效果不会弱于任何一个使用单个核函数的过程监测模型。综合这两点优势，本发明方法克服了传统基于KPCA的过程监测方法的不足，是一种更为优选的非线性过程监测方法。

附图说明

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

图2本发明方法与传统KPCA方法在TE过程21类故障上的平均故障检测率对比。

具体实施方式

下面结合附图对本发明方法进行详细的说明。

如图1所示，本发明公开一种基于多核主元分析模型的非线性过程监测方法。现结合一个具体的实施案例来陈述本发明方法的具体实施方式。

所测试的过程对象为TE过程，该过程原型是伊斯曼化工生产车间的一个实际工艺流程。目前，TE过程因其流程的复杂性，已作为一个标准实验平台被广泛用于故障检测研究。整个TE过程包括22个测量变量、12个操作变量、和19个成分测量变量。所采集的数据分为22组，其中包括1组正常工况下的数据集与21组故障数据。而在这些故障数据中，有16个是已知故障类型，如冷却水入口温度或进料成分的变化、阀门粘滞、反应动力学漂移等，还有5个故障类型是未知的。为了对该过程进行监测，选取如表1所示的33个过程变量，接下来结合该TE过程对本发明具体实施步骤进行详细的阐述。

表1：TE过程监测变量。

序号	变量描述	序号	变量描述	序号	变量描述
						1	物料A流量	12	分离器液位	23	D进料阀门位置
2	物料D流量	13	分离器压力	24	E进料阀门位置
						3	物料E流量	14	分离器塔底流量	25	A进料阀门位置
4	总进料流量	15	汽提塔等级	26	A和C进料阀门位置
						5	循环流量	16	汽提塔压力	27	压缩机循环阀门位置
6	反应器进料	17	汽提塔底部流量	28	排空阀门位置
						7	反应器压力	18	汽提塔温度	29	分离器液相阀门位置
8	反应器等级	19	汽提塔上部蒸汽	30	汽提塔液相阀门位置
						9	反应器温度	20	压缩机功率	31	汽提塔蒸汽阀门位置
10	排空速率	21	反应器冷却水出口温度	32	反应器冷凝水流量
						11	分离器温度	22	分离器冷却水出口温度	33	冷凝器冷却水流量

步骤(1)：在生产过程正常运行状态下，采集n＝500个样本数据组成训练数据矩阵X∈R^500×33，并对X中各列实施标准化处理，得到标准化后的矩阵

步骤(2)：设置高斯核函数的参数g＝10×33＝330后，根据上述公式①计算高斯核矩阵K_G∈R^500×500：

步骤(3)：设置多项式核函数的参数p＝3后，根据上述公式②计算多项式核矩阵K_P∈R^500×500：

步骤(4)：根据上述公式③计算Sigmoid核矩阵K_S∈R^500×500：

步骤(5)：根据上述公式④计算线性核矩阵K_L∈R^500×500：

步骤(6)：依据上述公式⑤对核矩阵K_G、K_P、K_S、和K_L实施中心化处理，对应得到中心化后的核矩阵

和

步骤(7)：计算高斯核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_G为非零特征值的个数。

步骤(8)：计算多项式核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_P为非零特征值的个数。

步骤(9)：计算Sigmoid核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_S为非零特征值的个数。

步骤(10)：计算线性核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_L为非零特征值的个数。

步骤(11)：设置核主元个数为d＝36，从而得到高斯核主元变换矩阵

多项式核主元变换矩阵

Sigmoid核主元变换矩阵

以及线性核主元变换矩阵

步骤(12)：先根据公式

计算核主元得分矩阵Φ_G、Φ_P、Φ_S、和Φ_L，再根据公式Λ_c＝Φ_c ^TΦ_c/(n-1)计算各核主元得分矩阵的协方差矩阵Λ_G、Λ_P、Λ_S、和Λ_L。

步骤(13)：根据上述公式⑥与公式⑦计算核主元监测统计量的控制上限D_lim和残差监测统计量的控制上限Q_G，lim、Q_P，lim、Q_S，lim、以及Q_L，lim，F分布与卡方分布的具体数值可通过查概率表获得。

采集TE过程21种不同故障工况下的采样数据，每组采样数据包含960个样本，其中前160个样本为正常工况下采样，后800个样本才是故障工况下的采样。然后，执行如下所示步骤(14)至步骤(18)实施在线故障检测。

步骤(14)：对各个样本数据x_new∈R^1×m进行标准化处理得到

步骤(15)：根据前述公式⑧至公式

分别计算高斯核向量z_G∈R^1×500、多项式核向量z_P∈R^1×500、Sigmoid核向量z_S∈R^1×500、以及线性核向量z_L∈R^1×500：

上式中，z_G，i、z_P，i、z_S，i、和z_L，i分别为核向量z_G、z_P、z_S、和z_L中的i个元素。

步骤(16)：依据前述公式

计算对各核向量实施中心化处理得到

和

步骤(17)：根据公式

计算核主元得分向量t_G、t_P、t_S、和t_L后，再根据前述公式

分别计算监测统计量D_c与Q_c。

步骤(18)：对于下标号c的所有四个取值，判断是否满足条件：D_c≤D_lim且Q_c≤Q_c，lim？若否，则当前采样时刻已进入故障工况；若是，则过程对象处于正常运行状态。

计算21组不同类型故障工况下测试数据的故障检测成功率，然后求取均值后将本发明方法的故障检测成功率与传统使用单个核函数(高斯核函数)的KPCA方法的故障检测成功率对比于图2中。从图2可以发现，本发明方法的故障检测成功率明显高于传统KPCA方法。

上述实施例仅是对本发明的优选实施方式，在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改和改变，不应排除在本发明的保护范围之外。

Claims (1)

1.一种基于多核主元分析模型的非线性过程监测方法，其特征在于，包括以下步骤：

首先，离线建模阶段包括如下所示步骤(1)至步骤(13)；

步骤(1)：在生产过程正常运行状态下，采集n个样本数据组成训练数据矩阵X∈R^n×m，并对X中各列实施标准化处理，得到标准化后的矩阵

其中，

表示标准化后的第i个样本数据，下标号i＝1，2，…，n，m为测量变量个数，上标号T表示矩阵或向量的转置，R为实数集，R^n×m表示n×m维的实数矩阵；

步骤(2)：设置高斯核函数的参数g后，按照下式①计算高斯核矩阵K_G∈R^n×n：

上式中，K_G(i，j)表示高斯核矩阵K_G中的第i行、第j列元素，下标号i＝1，2，…，n与j＝1，2，…，n，exp表示以自然常数e为底的指数函数，符号|| ||表示计算向量的长度；

步骤(3)：设置多项式核函数的参数p后，根据下式②计算多项式核矩阵K_P∈R^n×n：

上式中，K_P(i，j)表示多项式核矩阵K_P中的第i行、第j列元素；

步骤(4)：根据如下所示公式③计算Sigmoid核矩阵K_S∈R^n×n：

上式中，K_S(i，j)表示Sigmoid核矩阵K_S中的第i行、第j列元素，tanh表示双曲正切函数；

步骤(5)：根据如下所示公式④计算线性核矩阵K_L∈R^n×n：

上式中，K_L(i，j)表示线性核矩阵K_L中的第i行、第j列元素；

步骤(6)：依据如下公式分别对核矩阵K_G、K_P、K_S、和K_L实施中心化处理，对应得到中心化后的核矩阵

和

上式中，下标号c∈{G，P，S，L}分别表示高斯核、多项式核、Sigmoid核、以及线性核，方阵Θ∈R^n×n中各元素都等于1；

步骤(7)：计算高斯核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_G为非零特征值的个数；

步骤(8)：计算多项式核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_P为非零特征值的个数；

步骤(9)：计算Sigmoid核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_S为非零特征值的个数；

步骤(10)：计算线性核矩阵

所有非零特征值

所对应的特征向量

此步骤要求各个特征向量皆为单位长度，N_L为非零特征值的个数；

步骤(11)：保留高斯核主元变换矩阵

多项式核主元变换矩阵

Sigmoid核主元变换矩阵

以及线性核主元变换矩阵

步骤(12)：设置核主元个数为d后，根据公式

计算核主元得分矩阵Φ_G、Φ_P、Φ_S、和Φ_L，其中A_c(1：d)表示矩阵A_c中的第1列至第d列向量组成的矩阵，再根据公式Λ_c＝Φ_c ^TΦ_c/(n-1)计算各核主元得分矩阵的协方差矩阵Λ_G、Λ_P、Λ_S、和Λ_L；

步骤(13)：根据如下所示公式分别计算核主元监测统计量的控制上限D_lim和残差监测统计量的控制上限Q_G，lim、Q_P，lim、Q_S，lim、以及Q_L，lim：

上式中，g_c与h_c分别为卡方分布的权重系数与自由度，F_d，n-d，σ表示自由度为d与n-d的F分布在置信度σ＝99％下的取值，

表示自由度为h_c的卡方分布在置信度σ＝99％下的取值；

其次，完成上述离线建模阶段后，即可实施在线故障监测，具体包括如下所示步骤(14)至步骤(18)；

步骤(14)：在线采集最新采样时刻的样本数据x_new∈R^1×m，并对其进行标准化处理得到

步骤(15)：根据如下所示公式分别计算高斯核向量z_G∈R^1×n、多项式核向量z_P∈R^1×n、Sigmoid核向量z_S∈R^1×n、以及线性核向量z_L∈R^1×n：

上式中，z_G(i)、z_P(i)、z_S(i)、和z_L(i)分别为核向量z_G、z_P、z_S、和z_L中的i个元素；

步骤(16)：依据下式计算对各核向量实施中心化处理得到

和

其中，行向量θ＝[1，1，…，1]∈R^1×n中各元素都等于1；

步骤(17)：根据公式

计算核主元得分向量t_G、t_P、t_S、和t_L后，根据如下所示公式分别计算监测统计量D_c与Q_c：

上式中，t_c(1：d)表示向量t_c中第1列至第d列元素组成的向量，t_c(d+1：N_c)表示向量t_c中第d+1列至第N_c列元素组成的向量，下标号c∈{G，P，S，L}；

步骤(18)：对于下标号c的所有四个取值，判断是否满足条件：D_c≤D_lim且Q_c≤Q_c，lim；若否，则当前采样时刻已进入故障工况；若是，则过程对象处于正常运行状态，返回步骤(14)实施对下一采样时刻的过程监测。

Images