CN108958226A - 基于生存信息势—主成分分析算法的te过程故障检测方法 - Google Patents

Abstract

本发明涉及非高斯随机系统故障诊断领域，具体为基于生存信息势—主成分分析算法的TE过程故障检测方法，解决TE过程不服从高斯分布且变量个数较多，传统的方法结果不精确的问题，步骤：一、采集正常数据归一化处理；二、计算平均重构误差；三、构建SIP‑PCA算法；四、估计SPE及T²统计量的控制限；五、采集故障数据与控制限对比，检测是否发生故障。优点：计算复杂度低，迭代不涉及指数运算，收敛速度快；无分布假设，更具一般性；控制限的计算过程中采用核密度估计法，充分考虑非高斯噪声的影响，提高了故障检测的精度；使用SIP使变量的随机性变小，在非高斯噪声干扰下的估计结果具有更好的平滑性。

Description

基于生存信息势—主成分分析算法的TE过程故障检测方法

技术领域

本发明涉及非高斯随机系统故障诊断领域，具体为基于生存信息势(SurvivalInformation Potential，SIP)—主成分分析(Principal Component Analysis,PCA)算法的田纳西—伊斯曼(Tennessee Eastman，TE)过程故障检测方法。

背景技术

化工过程具有技术密集、工艺复杂、生产过程联系紧密等特点。化工过程如果轻视安全的重要性，就可能造成重大事故，这样不仅会影响生产而且会造成巨大的经济损失甚至是人员的伤亡，如果大量的有毒有害的物质发生泄漏将还会对环境造成重大的污染。为了确保生产过程运行正常，降低设备的维修率，保证产品的生产质量，需要对生产过程进行全程监测，所以故障诊断成为有实际意义的研究课题。

随着故障检测和诊断在内的多变量统计方法得到了迅速发展。许多于从海量数据集中提取过程信息并解释这些信息的统计技术已在各个领域得到发展。这种方法不需要准确的经验知识也不需要复杂数学模型的建立，而且釆用的数据是工业过程中的第一手资料，因此更接近真实情况。附加成本低，假定前提少，且极易维护，因此适用于复杂的工业过程。已证明其中特别有效的一种方法是使用PCA结合T²图表和SPE图表。

但是传统的基于PCA的监测都假设数据服从高斯分布，实际上，并非所有的过程变量都是高斯分布的，其中一些可能遵循不同类型的非高斯分布，例如TE过程，该过程是依据实际化工联合反应过程开发的模型仿真平台，这个过程模型包括42个测量变量和12个控制变量，其产生的数据具体时变、强耦合、非线性特征，且数据分布不服从高斯分布。在这种情况下，传统的基于PCA的监控方法可能无法很好地运行，可能会导致过程监控的不准确。因此，为了解决以上这些问题，本发明提出一种基于生存信息势—主成分分析算法的TE过程故障检测方法。

发明内容

本发明旨在解决TE过程不服从高斯分布且变量个数较多，采取传统的主成分分析方法会出现监测结果不精确的问题，提出了一种基于生存信息势—主成分分析算法的TE过程故障检测方法。

本发明是通过以下操作步骤实现的：基于生存信息势—主成分分析算法的TE过程故障检测方法，包括以下操作步骤：

一、采集正常工况下的数据X∈R^N×m(x₁,x₂,...,x_m)作为正常数据，并且对其进行归一化处理，其中N表示样本数，m表示变量个数；

二、从重构误差的角度考虑，PCA就是寻找一个投影矩阵，使得平均重构误差最小，平均重构误差的计算公式如下：

即

其中x_i为任一样本点，为样本的投影点；t_i＝P^Tx_i是m×1维的得分向量，即为主成分向量；P为加载矩阵，即为主成分的投影方向；

三、基于SIP—PCA的算法：

1.构建SIP—PCA算法：SIP的本质是在传统的信息势表达式中，用生存函数S(x)＝P{X＞x}＝1-F(x)取代概率密度函数；由此，对于一个向量其α-阶(α>0)生存信息势的定义为：

其中，是随机向量|X|的多元生存函数；|X|代表随机向量X的绝对值，是一个包含元素|X₁|,|X₂|,K,|X_m|的m维随机向量；

为了计算简便，引入了经验SIP，假设有N个样本X(1)，X(2),...,X(N)，在每个样本点的经验生存函数为：

因此，经验SIP可以表示为：

当样本数据为标量数据时，即m＝1，对于N个样本，0≤X(1)≤X(2)≤…≤X(N)，可以得到：

假设X(0)＝0，可以得到一种更一般的形式：

其中：

由于式7在x(i)＝0处不平滑，故使用有序样本的平方0≤X²(1)≤X²(2)≤…≤X²(N)来计算经验SIP，即：

然后，我们将式2带入式9，得到SIP-PCA的性能指标：

其中，P要求为正交矩阵，最终得出：

其中，t_j＝P^Tx_j，P＝[p₁,p₂,...,p_m]；

2.采用梯度下降法对性能指标进行优化，得到所有方向的方向矩阵P：

PCA是一种统计投影技术，其本质是坐标变换，由上述建立SIP—PCA算法，当重构误差最小时说明已经在保证数据不失真的情况下，将原始数据投影到了新的坐标系中；

SIP是衡量估计结果随机性的度量，即当SIP最小时，原始数据X与重构数据之间的误差最小；

因此，通过优化性能指标即可得到最优的投影矩阵P；由于性能指标式11是带约束的优化问题，本发明使用拉格朗日因子法，将上述带约束问题转换成无约束优化问题，如式12所示：

然后，采用梯度下降法对该性能指标进行优化，主要分为以下两步：

b1、对加载矩阵P和拉格朗日因子求偏导：

b2、根据递推公式求得未知参数P和

其中μ和η为学习率，且μ∈[0,1]，η∈[0,1]；其存在是为了均衡偏导部分在幅值上对所求最小值的影响；

3.使用交叉验证法确定PCA的主元个数：

使用PCA构造主元模型时，必须确定主元个数，而主元个数的确定需要考虑两个方面的因素：原始数据维数的降低和原始测量数据信息的丢失。主元个数的选择会直接影响故障诊断的结果，可见，主元个数的选择十分重要。常用的主元个数的选择方法有累计方差贡献率法、未重构方差法、交叉检验法等。因为前两种方法都需要人为设置一个阈值来确定要选取得主元个数，所以本发明采用交叉验证方法求取主元，具体步骤如下：

c1、将步骤一采集的原始数据X∈R^N×m分为两部分，去掉第i行得到X^(-i)∈R^(N-1)×m；以及只有第i行的x⁽ⁱ⁾∈R^1×m；

c2、使用数据X^(-i)∈R^(N-1)×m执行SIP-PCA算法，得到加载矩阵P^(-i)；

c3、将数据x⁽ⁱ⁾∈R^1×m投影到加载矩阵P^(-i)上求得其得分：

t⁽ⁱ⁾＝x⁽ⁱ⁾P^(-i) 式15

c4、确定x⁽ⁱ⁾∈R^1×m的残差：

重复执行c1-c4，直至得到所有的残差；

c5、计算预测残差平方和：

PRESS最小值所对应的j即为所求的最佳主元个数；

四、将置信度设置为α，使用核密度估计法来估计SPE统计量以及T²统计量的概率密度函数，通过对两个统计量的概率密度函数积分得到故障检测控制限，包括如下步骤：

1.T²统计量表示得分向量的平方和，其具体计算为：

T²＝x^TPP^Tx 式18

2.SPE统计量表示每个采样数据与统计模型之间的误差，具体的计算为：

SPE＝||(I-PP^T)x||² 式19

3.由于数据是非高斯的，仍使用传统基于高斯假设的控制限的求法，可能会导致结果不准确，因此使用核密度估计法估计SPE和T²的概率密度函数：

4.根据置信度α，对T²和SPE统计量的概率密度函数分别求积分，得到对应的控制限和CL_SPE：

五、采集故障工况下的数据作为故障数据，对其进行归一化处理，随后将故障数据投影到加载矩阵上，计算得到故障数据模型的SPE统计量以及T²统计量，通过与式22和式23计算出的控制限作对比，检测系统是否发生了故障。

与现有的技术相比本发明具有以下几个优势：实际化工过程中，并非所有的过程变量都是高斯分布的，其中一些可能遵循不同类型的非高斯分布。本方法使用了一种新的衡量估计结果随机性的性能指标—SIP，其本质是在传统的信息势表达式中，用生存函数取代概率密度函数。SIP不涉及核计算，使得其计算复杂度低，每次迭代不涉及指数运算，收敛速度快。SIP没有任何的分布假设，更具一般性；与传统的求取控制限的方法不同，本发明充分考虑到非高斯噪声的影响，在计算控制限的过程中采用核密度估计法，有效解决了传统PCA算法基于高斯分布的假设，提高了故障检测的精度。此外，使用SIP不仅可以使变量的随机性变小，还能使误差接近于零，且在非高斯噪声干扰下的估计结果具有相对较好的平滑性。

附图说明

图1为本发明SIP—PCA算法的故障诊断结构框图；

图2a为其中一个变量的正常数据分布图

图2b为图2a变量加入故障后的数据分布图；由图2b与图2a对比可以看出，未出现故障时，数据的分布范围为[-4,4]，当给定一个阶跃信号作为故障数据时，数据的分布范围变成了[0,8]，偏离了正常范围；

图3a为T²统计量的概率密度函数图；

图3b为SPE统计量的概率密度函数图；已知曲线与X轴所围成的面积等于1，这里通过设置统计阈置信度α来具体确定置信限：例如，此时置信度设为0.99，那么对应的控制限就是一个T²与SPE的值，在该值右边曲线与X轴围城的面积占总面积的0.01，即占总面积的α；此时，便求出了SIP-PCA的控制限与CL_SPE；

图4a为传统PCA算法的主元分析统计量变化图(故障检测结果图)；

图4b为SIP-PCA算法的主元分析统计量变化图(故障检测结果图)；由图4a与图4b对比可知，在250秒时人为加入故障，传统的PCA算法不能及时检测出故障，而本申请提出的SIP-PCA算法则在检测故障的快速性和准确性方面均显著优于传统PCA算法。

具体实施方式

以下结合具体实施例对本发明作进一步说明：在某化工生产厂的整个TE过程主要有五个操作单元组成：反应器、产品冷凝器、气液分离器、循环压缩机和汽提塔。气态的反应物进入到反应器中，生成液态产品气相的反应是在一种不挥发的气相催化剂的作用下进行的。反应器内置有冷凝包用来移除反应产生的热量。产品以气态的形式出来，并夹杂有一些未反应物。催化剂仍然滞留在反应器中。

该TE生产过程主要由四种气态物料参与反应，分别为A、C、D和E，生产出两种产品G、H，并伴有一种副产品F，此外在产品的进料中含有少量的惰性气体B。整个过程主要由四种反应组成，反应方程式如下：

此TE过程主要有12个控制变量和41个测量变量，在确定过程监测变量时，通常选取常用的与过程运行状况紧密相关的16个过程变量，如表1所示；作为正常数据对系统进行故障诊断，数据集分为两部分：训练数据和测试数据，其中采样时间设置为1s，训练样本个数为500，测试数据在250s引入故障。

表1 TE过程检测变量列表

表2所示为TE过程的故障列表，其列举了多种类型的故障，可以体现实际应用中可能会出现的各种故障。

表2 TE过程故障列表

该实施例采用基于生存信息势—主成分分析算法的TE过程故障检测方法，包括以下操作步骤：

一、采集正常工况下的数据X∈R^N×m(x₁,x₂,...,x_m)作为正常数据，并且对其进行归一化处理，其中N表示的样本数在本实施例中为16，m表示的变量个数为500；

二、计算平均重构误差：

即

其中x_i为任一样本点，为样本的投影点；t_i＝P^Tx_i是m×1维的得分向量，即为主成分向量；P为加载矩阵，即为主成分的投影方向；

三、基于SIP—PCA的算法：

1.构建SIP—PCA算法：

2.采用梯度下降法对性能指标进行优化，得到所有方向的方向矩阵P：

b1、对加载矩阵P和拉格朗日因子求偏导：

b2、根据递推公式求得未知参数P和

其中μ和η为学习率，且μ∈[0,1]，η∈[0,1]；

3.使用交叉验证法确定PCA的主元个数，num(num<16)，则原始数据的维度降低X∈R⁵⁰⁰×¹⁶→X∈R⁵⁰⁰×^num，具体步骤如下：

c1、将步骤一采集的原始数据X∈R^N×m分为两部分，去掉第i行得到X^(-i)∈R^(N-1)×m；以及只有第i行的x⁽ⁱ⁾∈R^1×m；

c2、使用数据X^(-i)∈R^(N-1)×m执行SIP-PCA算法，得到加载矩阵P^(-i)；

c3、将数据x⁽ⁱ⁾∈R^1×m投影到加载矩阵P^(-i)上求得其得分：t⁽ⁱ⁾＝x⁽ⁱ⁾P^(-i) 式15

c4、确定x⁽ⁱ⁾∈R^1×m的残差：

重复执行c1-c4，直至得到所有的残差；

c5、计算预测残差平方和：

PRESS最小值所对应的j即为所求的最佳主元个数；至此，完成了SIP-PCA算法的构建；

四、将置信度设置为α，使用核密度估计法来估计SPE统计量以及T²统计量的概率密度函数，通过对两个统计量的概率密度函数积分得到故障检测控制限，包括如下步骤：

1.T²统计量表示得分向量的平方和，其具体计算为：

T²＝x^TPP^Tx 式18

2.SPE统计量表示每个采样数据与统计模型之间的误差，具体的计算为：

SPE＝||(I-PP^T)x||² 式19

3.使用核密度估计法估计SPE和T²的概率密度函数：

4.根据置信度α，对T²和SPE统计量的概率密度函数分别求积分，得到对应的控制限和CL_SPE：

五、采集故障工况下的数据作为故障数据，对其进行归一化处理，随后将故障数据投影到加载矩阵上，计算得到故障数据模型的SPE统计量以及T²统计量，通过与式22和式23计算出的控制限作对比，检测系统是否发生了故障。

例如，给定一个数据表达式如下：

X＝UY+2E 式25

其中X∈R^500×7；Y是利用高斯混合模型得到的非高斯分布的数据；E是服从Beta分布的噪声，在250s时刻，给数据X的第二个变量一个阶跃信号作为故障数据，其故障诊断结果如图4b所示，从结果可以看到，在250s时，T²和SPE统计量的值均超过了控制限，说明系统出现了故障。

Claims (1)

1.一种基于生存信息势-主成分分析算法的TE过程故障检测方法，其特征在于：包括以下操作步骤：

一、采集正常工况下的数据X∈R^N×m(x₁，x₂，...，x_m)作为正常数据，并且对其进行归一化处理，其中N表示样本数，m表示变量个数；

二、计算平均重构误差：

即

其中x_i为任一样本点，为样本的投影点；t_i＝P^Tx_i是m×1维的分向量，即为主成分向量；P为加载矩阵，即为主成分的投影方向；

三、基于SIP-PCA的算法：

1.构建SIP-PCA算法：

2.采用梯度下降法对性能指标进行优化，得到所有方向的方向矩阵P：

b1、对加载矩阵P和拉格朗日因子求偏导：

b2、根据递推公式求得未知参数P和

其中μ和η为学习率，且μ∈[0，1]，η∈[0，1]；

3.使用交叉验证法确定PCA的主元个数，具体步骤如下：

c1、将步骤一采集的原始数据X∈R^N×m分为两部分，去掉第i行得到X^(-i)∈R^(N-1)×m；以及只有第i行的x⁽ⁱ⁾∈R^1×m：

c2、使用数据X^(-i)∈R^(N-1)×m执行SIP-PCA算法，得到加载矩阵P^(-i)；

c3、将数据x⁽ⁱ⁾∈R^1×m投影到加载矩阵P^(-i)上求得其得分：

t⁽ⁱ⁾＝x⁽ⁱ⁾P^(-i)

c4、确定x⁽ⁱ⁾∈R^1×m的残差：

重复执行c1-c4，直至得到所有的残差；

c5、计算预测残差平方和：

PRESS最小值所对应的j即为所求的最佳主元个数；至此，完成了SIP-PCA算法的构建；

四、将置信度设置为α，使用核密度估计法来估计SPE统计量以及T²统计量的概率密度函数，通过对两个统计量的概率密度函数积分得到故障检测控制限，包括如下步骤：

1.T²统计量表示得分向量的平方和，其具体计算为：

T²＝x^TPP^Tx

2.SPE统计量表示每个采样数据与统计模型之间的误差，具体的计算为：

SPE＝||(I-PP^T)x||²

3.使用核密度估计法估计SPE和T²的概率密度函数：

4.根据置信度α，对T²和SPE统计量的概率密度函数分别求积分，得到对应的控制限和CL_SPE：

五、采集故障工况下的数据作为故障数据，对其进行归一化处理，随后将故障数据投影到加载矩阵上，计算得到故障数据模型的SPE统计量以及T²统计量，通过与步骤四、4计算出的控制限作对比，检测系统是否发生了故障。