CN103412557A - 一种适于非线性过程在线监控的工业故障检测与诊断方法 - Google Patents

Abstract

本发明公开一种适于非线性过程在线监控的工业故障检测与诊断方法，其涉及模式识别领域，本质是对工业故障来源的分类识别，实现步骤如下：a)样本数据采集；b)信号预处理；c)数据特征提取；d)模式分类与故障诊断试验；e)结果输出；最后对临界故障状态进行报警及应急处理。本发明提供的工业故障检测与诊断方法在现有的技术基础上有效地提高了运行效率。其适用于非线性过程在线监控，能对故障源做出准确的检测与诊断。

Description

一种适于非线性过程在线监控的工业故障检测与诊断方法

技术领域

本发明是一种工业故障检测与诊断方法，具体地说，是涉及一种适于非线性过程在线监控的基于QR分解的核判别分析法(IKDR/QR)的工业故障检测与诊断方法。

背景技术

随着近年来计算机系统、数据库系统的普及应用，工厂拥有了相当丰富的生产数据资源。在工业生产过程中，不可避免的会发生设备损耗和故障的情况，不但影响产品的质量，有时还会发生重大生产事故，完全依靠人力的传统监测方法己不足以解决复杂的质量控制问题。

故障是指使系统至少一个特性或参数出现较大的偏差，并超出了可接受的范围。此时系统的性能明显低于正常水平，难以完成其预期的功能。故障对整个系统的影响，小到如某个变量值不能正常显示，大到导致整个系统的崩溃，甚至产生人身事故。故障类型一般分为三类：(1)被控过程(对象)的故障；(2)仪表器件的故障；(3)软件故障。

实际生产过程中的生产数据有线性的、非线性的、静态的、动态的等，对于不同的生产过程所具有的特征，应选用不同的故障检测方法，这样才能有效地检测到故障。

经专利查询统计，国内已有一些故障检测与诊断方面的专利：例如，基于非线性频谱分析的故障检测与诊断方法（99115737.0）、故障诊断服务器、车辆故障检测与诊断方法、装置及系统（200910090461.1）、基于无人飞行器的接触网故障检测与诊断方法（201110037895.2）等。上述发明专利很少处理非线性数据且不适于在线学习应用。本发明针对这两方面问题，提出一种适于非线性过程在线监控的工业故障检测与诊断方法。该方法针对于工业系统中的非线性生产数据。近年来，对于普遍存在的非线性问题，核方法的出现大大推动了非线性过程研究的发展。主要的核方法有核主成分分析(kernel principal component analysis，KPCA)；核Fisher判别分析(kernel fisher discriminator analysis，KFDA)；核独立成分分析(kernel independent component analysis，KICA)等。本发明本质上也是一种核判别分析方法。

发明内容

本发明要克服现有技术无法处理非线性数据且不适于在线学习应用的缺点，提供一种高效的故障检测与诊断方法，充分利用基于QR分解的核判别分析法的先降维后特征提取思想，将核空间映射到低维空间进行计算，减少了构造核矩阵的计算量，降低了核矩阵的存储空间，同时引入了增量计算思想。其具有非线性处理能力且能够快速有效地处理高维数据，有效地解决了在线学习中冗余计算的问题。

本发明解决其技术问题所采用的技术方案是：

一种适于非线性过程在线监控的工业故障检测与诊断方法，包括以下步骤：

a)样本数据采集：利用TE过程的仿真数据，仿真数据包括15组不同类别的故障数据，每组故障数据包含301个样本和51种故障测量值，构成4515×51的数据集；

b)信号预处理：首先利用小波变换软阈值法对样本数据进行小波去噪，其中小波变换软阈值法采用三层小波分解，并将小波重构后的信号作为原始样本输入；

c)数据特征提取：运用基于QR分解的增量式核判别分析法计算最佳投影矩阵，根据最近投影矩阵对样本对象进行特征降维，后续根据新的训练样本数据及时更新最佳投影矩阵；

d)故障源分类判别：采用欧式近邻法进行判别分类；

e)结果输出：输出分类信息以确定故障源状态；

最后对故障源状态进行报警及应急处理。

具体来讲，在本发明的工业故障检测与诊断方法中，所述数据特征提取过程具体包括下述步骤：

a1确定前一次基于QR分解的增量式核判别分析法计算中类间散度矩阵B^Φ中的矩阵部分，确定前一次基于QR分解的增量式核判别分析法计算中总体散度矩阵T^Φ中的

矩阵部分，确定每类样本数量n_i(i=1,2,…,c)，确定新插入的样本x；

a2根据步骤a1中的

矩阵部分、每类样本数量n_i以及新插入的样本x更新类间散度矩阵B^Φ；

a3根据步骤a1中的

矩阵部分、每类样本数量n_i以及新插入的样本x更新总体散度矩阵T^Φ；

a4利用步骤a2、步骤a3得到的类间散度矩阵B^Φ和总体散度矩阵T^Φ，对(T^Φ+μI)^-1B^Φ进行特征分解得到c个特征向量

令低维最佳投影矩阵 $V^{\mathrm{\Φ}}=[{\mathrm{\φ}}_1^{\mathrm{\Φ}},{\mathrm{\φ}}_2^{\mathrm{\Φ}},\·\·\·,{\mathrm{\φ}}_c^{\mathrm{\Φ}}];$

a5通过G^Φ=C^Φ(R^Φ)^-1V^Φ计算更新后的最佳投影矩阵G，其中C^Φ为核空间H中的质心矩阵，R^Φ为对C^Φ进行QR分解的R矩阵。

本发明的技术构思：复杂的非线性过程面临在线学习问题，在实际应用中，采集一组完整且具有代表性的数据是昂贵、费时的。在很多应用的初始阶段只有少量可供使用的样本数据，而随着时间的推移，新增样本会陆续呈现，需根据新的训练数据及时更新最佳投影矩阵。如果采用传统批量式的非线性特征提取方法，每次新样本到来都需要重新构造核矩阵K，并且需要将新样本与原有样本一起重新进行一次完整的投影矩阵计算，这样严重影响了特征提取的效率。因此本文提出了一种专门针对非线性系统在线学习效率问题的增量式非线性特征提取方法。这种方法有效地提高了更新投影矩阵的效率，同时减少核矩阵K的存储空间，使其更适用于非线性过程在线监控。基于QR分解的核判别分析(kerneldiscriminant analysis via QR decomposition，KDR/QR)是由Xiong等人提出的。它是一种利用核思想对非线性数据进行特征降维的方法。该方法的主要思想是，首先根据核方法将样本原始空间中的类间散度矩阵S_b和总体散度矩阵S_t映射到高维核空间中，然后在高维核空间中以类间散度矩阵

最大为目标再将和

映射到一个低维空间中，并在低维空间中应用Fisher判别分析求取低维空间的最佳投影矩阵，再利用QR分解求得原始空间的最佳投影矩阵。KDR/QR的优势是将核空间H映射到一个低维空间L后再进行判别分析，减小了特征分解的规模以提升特征提取效率。本文提出的基于QR分解的增量式核判别分析法(incrementalkernel discriminant analysis via QR decomposition，IKDR/QR)利用了这种先降维后特征提取的思想，在增量过程中充分利用低维空间L降低构造核矩阵K的规模，减少核矩阵K的存储空间。

本发明的优点是：具有非线性处理能力且能够快速有效地处理高维数据，有效地解决了在线学习中冗余计算的问题。

附图说明

图1是本发明的三层小波分解信号图；

图2是本发明的小波高频信号软阈值处理图；

图3是本发明的小波去噪后的采样单位振动信号图；

图4是本发明的高斯核分类性能实验结果图；

图5是本发明的多项式核分类性能实验结果图；

图6是本发明的特征提取效果比较图；

图7是本发明的训练时间比较图；

图8是本发明的流程图。

具体实施方式

下面对本发明作进一步描述。参照附图：

一种适于非线性过程在线监控的工业故障检测与诊断方法，包括以下步骤：

a)样本数据采集：Eastman化学公司开发的Tennessee Esatman过程模拟器提供了一个实际工业过程的仿真平台，近年来被广泛用于验证各种故障检测与诊断方法，成为一种国际上通用的标准仿真模型。利用TE过程的仿真数据进行仿真实验，仿真数据包括15组不同类别的故障数据，每组故障数据包含301个样本和51种故障测量值，构成4515×51的数据集。TE过程包括5个主要操作单元、4种气体进料、2个气液放热反应生成的2种主产品、2个衍生放热反应生成的2种副产品等，过程机理复杂，变量较多，包括阶跃、随机变化、慢漂移、粘滞和恒定位置等多种数据故障类型。其中TE过程的仿真数据可以从http://depts.washington.edu/control/下载得到；

b)信号预处理：对传感器采集的时域信号进行小波去噪，时频域转换以及提取适当频率上的幅值构成输入样本对象：

1)本发明采用三层小波分解，小波基选用db3。小波分解后得到一组低频信号成分和三组高频信号成分。如附图1所示，其中第一行是低频信号成分，其余三行是高频信号成分；

2)对三层高频成分运用软阈值方法进行处理，选取每层小波幅值的均值作为每层阈值δ，软阈值处理后的效果如附图2；

3)附图3为小波去噪后的采样单位振动信号图。

c)数据特征提取：运用基于QR分解的增量式核判别分析法计算最佳投影矩阵，根据最近投影矩阵对样本对象进行特征降维，后续根据新的训练样本数据及时更新最佳投影矩阵：

1)确定前一次基于QR分解的增量式核判别分析法计算中类间散度矩阵B^Φ中的矩阵部分，确定前一次基于QR分解的增量式核判别分析法计算中总体散度矩阵T^Φ中的

矩阵部分，确定每类样本数量n_i(i=1,2,…,c)，确定新插入的样本x；

2)根据步骤1)中的

矩阵部分、每类样本数量n_i以及新插入的样本x更新类间散度矩阵B^Φ，更新类间散度矩阵B^Φ时直接更新其中

矩阵部分，包括两个步骤：

a1更新M^TKM中非t类样本与t类样本的核关系：

${\tilde{A}}_{\mathrm{pt}}=\frac{1}{n_p{\tilde{n}}_t}\underset{a_i\∈c_p,a_j\∈{\tilde{c}}_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ij}}\right)$

$=\frac{1}{n_p(n_t+1)}(\underset{a_i\∈c_p,a_j\∈c_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ij}}\right)+\underset{a_i\∈c_p}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ix}}\right))$

$=\frac{n_t}{(n_t+1)}{A}_{\mathrm{pt}}+\frac{1}{n_p(n_t+1)}\underset{a_i\∈c_p}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ix}}\right)$

其中Φ(a)为映射Φ将数据a从原始空间映射到更高维的核空间H中，A_pt为前一次增量计算的结果，p≤c，p≠t。

a2更新M^TKM中t类样本的类内关系：

${\tilde{A}}_{\mathrm{tt}}=\frac{1}{{\tilde{n}}_t^2}\underset{a_{i,j}\∈{\tilde{c}}_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ij}}\right)$

$=\frac{1}{{(n_t+1)}^2}(\underset{a_{i,j}\∈c_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ij}}\right)+2\underset{a_i\∈c_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ix}}\right)+\mathrm{\Φ}\left(a_{\mathrm{xx}}\right))$

$=\frac{n_t^2}{{(n_t+1)}^2}{A}_{\mathrm{tt}}+\frac{2}{{(n_t+1)}^2}\underset{a_i\∈c_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ix}}\right)+\frac{1}{{(n_t+1)}^2}\mathrm{\Φ}\left(a_{\mathrm{xx}}\right)$

其中A_tt为前一次增量计算的结果。

3)根据步骤1)中的

矩阵部分、每类样本数量n_i以及新插入的样本x更新总体散度矩阵T^Φ，更新总体散度矩阵T^Φ时随着新的样本x的插入，矩阵KM产生两种改变：

a1更新每个样本与t类样本的核关系：

${\tilde{A}}_{\mathrm{pt}}=\frac{1}{{\tilde{n}}_t}\underset{a_j\∈{\tilde{c}}_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{pj}}\right)$

$=\frac{1}{n_t+1}(\underset{a_j\∈c_t}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{pj}}\right)+\mathrm{\Φ}\left(a_{\mathrm{px}}\right))$

$=\frac{n_t}{n_t+1}{A}_{\mathrm{pt}}+\frac{1}{n_t+1}\mathrm{\Φ}\left(a_{\mathrm{px}}\right)$

其中p≤c，p≠t。

a2插入代表新样本x核空间关系向量：

$\tilde{K}\tilde{M}=(\·\·\·;\frac{1}{{\tilde{n}}_1}\underset{a_i\∈c_1}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ix}}\right),\·\·\·,\frac{1}{{\tilde{n}}_c}\underset{a_i\∈c_c}{\mathrm{\Σ}}\mathrm{\Φ}\left(a_{\mathrm{ix}}\right);\·\·\·)$

其中向量

即为x样本与每类样本的核空间关系。

4)利用步骤2)、步骤3)得到的类间散度矩阵B^Φ和总体散度矩阵T^Φ，对(T^Φ+μI)^-1B^Φ进行特征分解得到c个特征向量

令低维最佳投影矩阵 $V^{\mathrm{\Φ}}=[{\mathrm{\φ}}_1^{\mathrm{\Φ}},{\mathrm{\φ}}_2^{\mathrm{\Φ}},\·\·\·,{\mathrm{\φ}}_c^{\mathrm{\Φ}}];$

5)通过G^Φ=C^Φ(R^Φ)^-1V^Φ计算更新后的最佳投影矩阵G，其中C^Φ为核空间H中的质心矩阵，R^Φ为对C^Φ进行QR分解的R矩阵。

d)故障源分类判别：采用最近邻分类器对故障源部位类型进行判别分类；

e)结果输出：输出分类信息以确定故障源状态。

最后对故障源状态进行报警及应急处理。

分类性能实验选择KFDA算法以及KDR/QR算法在不同数量的训练数据下比较降维效果以反映在线学习问题中训练样本数量的变化对算法性能的影响。从每类样本中选择20、50、100、200个样本，作为四组实验的训练样本，并选择每类50，共750个样本作为测试数据。核函数分别选用高斯核函数k(x,y)=exp(-||x-y||²/σ)，核参数σ取400，以及多项式核函数k(x,y)=(<x,y>+1)^d，d取2。通过5-fold交叉验证法选择适合的降维维度为30。选用欧式近邻法作为分类算法，每组实验各进行10次，取平均值作为结果。实验结果如图4，图5所示。从图4，图5可知，对于两种不同的核函数，KDR/QR相比KFDA有着更好的降维效果。随着训练样本的增多，KFDA的降维效果慢慢接近KDR/QR，但在训练样本较少时KDR/QR的优势明显。这一优势特别适用于在线学习中初始训练样本较少的情况，由于本发明提出的IKDR/QR每一次更新投影矩阵的计算都建立在前一次计算上，初始的训练效果会影响到后续的增量结果，所以KDR/QR适合增量化。

以下分析比较IKDR/QR与KDR/QR在线学习问题上的算法性能，包括特征提取效果和运算速度两方面。在线学习问题的特点是初始训练数据较少，随着训练的进行会有新的样本到来，因此我们从每类样本中选取20个，共300个作为初始训练样本，同时选取每类80个样本，共1200个作为待插入样本。核函数选择高斯核函数，将1200个样本依次插入训练样本，实验结果如图7，图8所示。从图7可知，增量式算法IKDR/QR与批量式算法KDR/QR的识别正确率相近，这证明了IKDR/QR算法在非线性系统在线学习过程中对样本空间特性的有效保持，其在逐次更新的过程中相比进行完整降维计算的批量式方法没有损失较多的信息。另一方面，从图8可知随着新插入的样本数量增多，IKDR/QR算法的时间开销一直维持在一个较低水平，而批量式方法每次计算投影矩阵的耗时快速递增，并随着训练样本的增加，增量式方法在速度上的优势越明显。为了更好反映出IKDR/QR在计算效率上的优越性，表1记录了插入新的训练样本后训练时间开销的变化。由表可知，当插入150个样本时，两者训练时间的差距只有0.5秒，而当插入1200个样本时，两者训练时间的差距达到了16.6秒。

表1训练时间开销（秒）

Claims (4)

1.一种适于非线性过程在线监控的工业故障检测与诊断方法，其特征在于，包括以下步骤：

a)样本数据采集：利用TE过程的仿真数据，仿真数据包括15组不同类别的故障数据，每组故障数据包含301个样本和51种故障测量值，构成4515×51的数据集；

b)信号预处理：利用小波变换软阈值法对样本数据进行小波去噪，其中小波变换软阈值法采用三层小波分解，并将小波重构后的信号作为原始样本输入；

c)数据特征提取：运用基于QR分解的增量式核判别分析法计算最佳投影矩阵，根据最近投影矩阵对样本对象进行特征降维，后续根据新的训练样本数据及时更新最佳投影矩阵；

d)故障源分类判别：采用欧式近邻法进行判别分类；

e)结果输出：输出分类信息以确定故障源状态；

最后对故障源状态进行报警及应急处理。

2.根据权利要求1所述的工业故障检测与诊断方法，其特征在于，在所述步骤c)中，运用基于QR分解的增量式核判别分析法计算最佳投影矩阵包括下述五个步骤：

a1确定前一次基于QR分解的增量式核判别分析法计算中类间散度矩阵B^Φ中的

矩阵部分，确定前一次基于QR分解的增量式核判别分析法计算中总体散度矩阵T^Φ中的

矩阵部分，确定每类样本数量n_i(i=1,2,…,c)，确定新插入的样本x；

a2根据步骤a1中的

矩阵部分、每类样本数量n_i以及新插入的样本x更新类间散度矩阵B^Φ；

a3根据步骤a1中的矩阵部分、每类样本数量n_i以及新插入的样本x更新总体散度矩阵T^Φ；

a4利用步骤a2、步骤a3得到的类间散度矩阵B^Φ和总体散度矩阵T^Φ，对(T^Φ+μI)^-1B^Φ进行特征分解得到c个特征向量

令低维最佳投影矩阵 $V^{\mathrm{\&Phi;}}=[{\mathrm{\&phi;}}_1^{\mathrm{\&Phi;}},{\mathrm{\&phi;}}_2^{\mathrm{\&Phi;}},\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&phi;}}_c^{\mathrm{\&Phi;}}];$

a5通过G^Φ=C^Φ(R^Φ)^-1V^Φ计算更新后的最佳投影矩阵G，其中C^Φ为核空间H中的质心矩阵，R^Φ为对C^Φ进行QR分解的R矩阵。

3.根据权利要求2所述的工业故障检测与诊断方法，其特征在于，在所述步骤a2中更新类间散度矩阵B^Φ时直接更新其中

矩阵部分，包括两个步骤：

a1更新M^TKM中非t类样本与t类样本的核关系：

${\tilde{A}}_{\mathrm{pt}}=\frac{1}{n_p{\tilde{n}}_t}\underset{a_i\&Element;c_p,a_j\&Element;{\tilde{c}}_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ij}}\right)$

$=\frac{1}{n_p(n_t+1)}(\underset{a_i\&Element;c_p,a_j\&Element;c_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ij}}\right)+\underset{a_i\&Element;c_p}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ix}}\right))$

$=\frac{n_t}{(n_t+1)}{A}_{\mathrm{pt}}+\frac{1}{n_p(n_t+1)}\underset{a_i\&Element;c_p}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ix}}\right)$

其中Φ(a)为映射Φ将数据a从原始空间映射到更高维的核空间H中，A_pt为前一次增量计算的结果，p≤c，p≠t。

a2更新M^TKM中t类样本的类内关系：

${\tilde{A}}_{\mathrm{tt}}=\frac{1}{{\tilde{n}}_t^2}\underset{a_{i,j}\&Element;{\tilde{c}}_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ij}}\right)$

$=\frac{1}{{(n_t+1)}^2}(\underset{a_{i,j}\&Element;c_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ij}}\right)+2\underset{a_i\&Element;c_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ix}}\right)+\mathrm{\&Phi;}\left(a_{\mathrm{xx}}\right))$

$=\frac{n_t^2}{{(n_t+1)}^2}{A}_{\mathrm{tt}}+\frac{2}{{(n_t+1)}^2}\underset{a_i\&Element;c_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ix}}\right)+\frac{1}{{(n_t+1)}^2}\mathrm{\&Phi;}\left(a_{\mathrm{xx}}\right)$

其中A_tt为前一次增量计算的结果。

4.根据权利要求2所述的工业故障检测与诊断方法，其特征在于，在所述步骤a3中更新总体散度矩阵T^Φ时随着新的样本x的插入，矩阵KM产生两种改变：

a31更新每个样本与t类样本的核关系：

${\tilde{A}}_{\mathrm{pt}}=\frac{1}{{\tilde{n}}_t}\underset{a_j\&Element;{\tilde{c}}_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{pj}}\right)$

$=\frac{1}{n_t+1}(\underset{a_j\&Element;c_t}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{pj}}\right)+\mathrm{\&Phi;}\left(a_{\mathrm{px}}\right))$

$=\frac{n_t}{n_t+1}{A}_{\mathrm{pt}}+\frac{1}{n_t+1}\mathrm{\&Phi;}\left(a_{\mathrm{px}}\right)$

其中p≤c，p≠t。

a32插入代表新样本x核空间关系向量：

$\tilde{K}\tilde{M}=(\&CenterDot;\&CenterDot;\&CenterDot;;\frac{1}{{\tilde{n}}_1}\underset{a_i\&Element;c_1}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ix}}\right),\&CenterDot;\&CenterDot;\&CenterDot;,\frac{1}{{\tilde{n}}_c}\underset{a_i\&Element;c_c}{\mathrm{\&Sigma;}}\mathrm{\&Phi;}\left(a_{\mathrm{ix}}\right);\&CenterDot;\&CenterDot;\&CenterDot;)$

其中向量

即为x样本与每类样本的核空间关系。

Images