CN101158873A - 一种非线性过程故障诊断方法 - Google Patents

Abstract

一种非线性过程故障诊断方法，包括采集数据、进行相似性分析、利用核主元分析对数据进行白化处理，求解白化后的观测变量z、利用修正独立元分析ICA提取独立元、利用T²和SPE统计量和LS－SVM进行故障检测与诊断步骤。本发明提出了非线性动态过程故障诊断技术，结合Kernel、ICA和LS－SVM三者的优点，即，发挥Kernel对非线性的表达能力，同时发挥ICA对动态特性的把握能力和LS－SVM的分类能力。

Description

一种非线性过程故障诊断方法

技术领域

本发明属于故障检测与诊断技术领域，特别涉及一种基于改进的核独立元分析和支持向量机相结合的非线性过程诊断方法。

背景技术

随着计算机和电子测量技术的飞速发展，现代工业过程大都具有完备的传感测量装置，可以在线获得大量的过程数据。对这些数据的统计分析可以帮助操作人员及时发现过程故障，避免重大事故的发生，这就促进了人们研究基于数据分析的过程监测方法。目前该类方法中使用最多的是主元分析(PCA)、偏最小二乘(PLS)及独立元分析(ICA)等。然而这些方法都是基于线性变换的统计方法，即假定过程变量之间满足线性相关关系。但是对于复杂的工业过程，像化工业的田纳西-伊斯曼过程、废水处理过程以及矿业中的连续采煤机采煤过程等等，变量之间往往呈现出非线性的相关关系特性，利用传统的统计方法将给过程监测带来较大的误差，增大过程故障误报、漏报的概率。

针对这种情况，Kramer在《AIChE Journal》期刊中提出了一种基于自联想的五层(输入，映射，瓶颈，解映射，及输出层)神经网络的非线性PCA。Dong和McAvoy在《Computers and ChemicalEngineering》期刊中提出了一种基于主曲线和神经网络的非线性PCA，并把它应用到非线性过程监测。基于遗传编程(genetic programming)以及输入训练神经网络的非线性PCA方法也已经分别在《Computers and Chemical Engineering》和《AIChE Journal》期刊中被提出。但是大部分现有的非线性PCA方法都是基于神经网络，必须离线训练，要计算主元的话也必须解决非线性最优化问题，并且在训练神经网络之前必须事先明确主元的个数。而且这些方法使用的非线性变换函数一般难以获得，神经网络的训练也较为困难。目前，核理论在非线性工业过程中的应用逐渐增多。对于非线性过程监测和故障诊断，引入了核主元分析(KPCA)和核独立元分析(KICA)。KICA首次是由Lee和Qin等人提出的。它是一种KPCA和ICA结合的新的非线性扩展。KPCA在一个与输入空间非线性相关的高维特征空间中计算主元。ICA提取高阶统计量以及分解观测数据。因此，独立元(ICs)比主元(PCs)从观测数据中展露了更多动态信息。同其它的非线性方法相比，KICA结合了KPCA和ICA的优点，成为一种在线监测故障的非线性动态方法。目前，把KICA方法用于工业过程存在的问题是：形成的核矩阵的维数是采样次数的平方，这样对于过程非常不利。ICA方法是用于动态过程中提取独立元，但是由于输入空间和特征空间不能像PCA那样自如地转换，使得故障诊断变得困难。

发明内容

针对现有技术存在的问题，本发明提供一种非线性过程故障诊断方法。对核独立元分析(KICA)方法进行了改进，并同最小二乘支持向量机(LS-SVM)结合，提出了一种基于改进的KICA和LS-SVM的非线性过程故障诊断方法。

针对故障检测与诊断，本发明在KICA的基础上进行了如下改进：

1)能够表达所有采样数据信息的维数适当的新的核矩阵的形成技术，提出了一个组群中数据的相似性概念。经过相似性分析之后，提取了特征空间中非线性信息来表达所有观测数据，这样当训练样本个数很大时，就降低了KPCA的计算量。

2)基于LS-SVM提出了一种故障分类策略。

核理论已经越来越多的应用于非线性工业过程，表明了核函数在工业过程对观测数据非线性化的优越性。建模之前必须进行数据相似性分析，原因如下：

1)映射数据进入特征空间变得线性冗余；

2)假如数据是线性的，那么当使用核技巧非线性化时将会产生更大的误差；

3)在KPCA训练过程中，核矩阵的大小是样本个数的平方。当样本数变大，特征值和特征向量的计算将会十分耗费时间。

本发明通过分别探索输入空间和特征空间样本的线性相关性来分析观测数据。并引入了相似性概念来定义新数据和节点集之间的相似性因子。

改进的KICA和LS-SVM均利用mercer核，两者很容易进行结合。本发明用LS-SVM进行故障诊断，一旦检测到故障，来自改进的KICA的核变换得分将直接被引入，作为LS-SVM的输入对故障进行分类。

本发明的故障诊断方法包括如下步骤：

步骤一、采集数据

采集过程中相关变量的数据，对于每个故障，产生两组数据，即训练数据和实时工况数据。训练数据用于建立模型，实时工况数据用于在线监测。并用均值和标准偏差规范化采集的数据。

步骤二、进行相似性分析

相似性分析包括输入空间的相似性分析和特征空间的相似性分析。本发明根据相似性因子对数据进行分析，去除掉相似性较强的数据。既解决了特征空间中数据的线性冗余，又降低了计算负载，还减小了使用核技巧带来的误差。

1、输入空间的相似性分析

在输入空间中，假如三个数据点x，y，z满足

$\cos \mathrm{\θ}=\frac{\mathrm{cov}(x-y,x-z)}{\mathrm{std}(x-y)\mathrm{std}(x-z)}=1---\left(1\right)$

那么x，y，z之间是线性关系。

学习的开始，只有两个数据点，即 $N_1=\{{\tilde{x}}_1,{\tilde{x}}_2\}=\{x_1,x_2\}.$ 第i个节点集被表示为 $N_i=\{{\tilde{x}}_1,\·\·\·,{\tilde{x}}_{n_i}\},$ i＝1，...，t₁，且满足n_i＜t₁，t₁为训练数据的个数。根据方程(1)定义相似性因子为

$S_1={\left(\frac{\mathrm{cov}(x_{i+1}-{\tilde{x}}_{k_1},x_{i+1}-{\tilde{x}}_{l_1})}{\mathrm{std}(x_{i+1}-{\tilde{x}}_{k_1})\mathrm{std}(x_{i+1}-{\tilde{x}}_{l_1})}\right)}^2---\left(2\right)$

式中  x_i+1-训练数据中的第i+1个的数据点；

-节点集N_i中任意两个数据点；

S₁满足0≤S₁≤1。当S₁＝1时相似性最强。因此，将x_i+1同节点集N_i中所有数据都两两进行相似性计算，假如所有 $S_1<\sqrt{{\mathrm{\&gamma;}}_0},$ 其中γ₀是满足 $0\&le;\sqrt{{\mathrm{\&gamma;}}_0}\&le;1$ 的预定义小变量，新数据x_i+1将被引入建立新扩展节点集N_i+1＝{N_i，x_i+1}；否则，这个数据点将被拒绝，即N_i+1＝N_i。在输入空间相似性分析后，获得第一个数据子集Ω₁。

2、特征空间中的相似性分析

得到子集Ω₁后，将Ω₁中的数据通过非线性函数Ф映射到特征空间，进行特征空间的相似性分析。令x，y及z是Ω₁中输入空间的三个数据点。映射三个数据点到特征空间，他们变为Ф(x)，Ф(y)及Ф(z)。

假如

$\cos \mathrm{\&theta;}=\frac{(\mathrm{\&Phi;}\left(z\right)-\mathrm{\&Phi;}\left(x\right),\mathrm{\&Phi;}\left(z\right)-\mathrm{\&Phi;}\left(y\right))}{\left|\right|\mathrm{\&Phi;}\left(z\right)-\left(x\right)\left|\right|\left|\right|\mathrm{\&Phi;}\left(z\right)-\mathrm{\&Phi;}\left(y\right)\left|\right|}=1---\left(3\right)$

那么

θ＝0

也就是，Ф(x)，Ф(y)及Ф(z)在特征空间都是线性的。

学习的开始，只有两个数据，即 $L_1=\{\mathrm{\&Phi;}\left({\tilde{x}}_1\right),\mathrm{\&Phi;}\left({\tilde{x}}_2\right)\}=\{\mathrm{\&Phi;}\left(x_1\right),\mathrm{\&Phi;}\left(x_2\right)\}.$ 特征空间第j个节点集被表示为 $L_j=\{\mathrm{\&Phi;}\left({\tilde{x}}_1\right),\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Phi;}\left({\tilde{x}}_{n_j}\right)\},$ j＝1，...，t₂，且满足n_j＜t₂，t₂为子集Ω₁中数据的个数。根据方程(3)定义特征空间的相似性因子为：

$S_2={\left(\frac{\mathrm{cov}(\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{k_2}\right),\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{l_2}\right))}{\mathrm{std}(\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{k_2}\right))\mathrm{std}(\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{l_2}\right))}\right)}^2$

式中Ф(x_j+1)-数据子集Ω₁中第j+1个数据点x_j+1映射后的数据；

-节点集L_j中的任意两个数据；

因而，S₂满足0≤S₂≤1。当S₂＝1时相似性是最强的。

S₂的可选择相似性因子能按下面过程来决定，方程(3)等价于方程(4)：

$\frac{f(x,y,z)}{g(x,y,z)}=1---\left(4\right)$

其中

f(x，y，z)＝(Ф(x)，Ф(y))-(Ф(y)，Ф(z))-(Ф(x)，Ф(z))+(Ф(z)，Ф(z))

${g(x,y,x)=((\mathrm{\&Phi;}\left(x\right),\mathrm{\&Phi;}\left(x\right))-2(\mathrm{\&Phi;}\left(x\right),\mathrm{\&Phi;}\left(z\right))+(\mathrm{\&Phi;}\left(z\right),\mathrm{\&Phi;}\left(z\right)))}^{\frac{1}{2}}((\mathrm{\&Phi;}\left(y\right),\mathrm{\&Phi;}\left(y\right))-$

$2(\mathrm{\&Phi;}\left(y\right),\mathrm{\&Phi;}\left(z\right))+{(\mathrm{\&Phi;}\left(z\right),\mathrm{\&Phi;}\left(z\right)))}^{\frac{1}{2}}---\left(5\right)$

方程(5)利用核技巧

k(a，b)＝(Ф(a)，Ф(b))    (6)

其中a，b为输入空间的两个数据点，将方程(4)变为

$\frac{k(x,y)-k(y,z)-k(x,z)+k(z,z)}{{(k(x,x)-2k(x,z)+k(z,z))}^{\frac{1}{2}}{(k(y,y)-2k(y,z)+k(z,z))}^{\frac{1}{2}}}=1---\left(7\right)$

简化方程(6)，本文采用径向基核函数 $k(x,y)=\exp \left(\frac{{\left|\right|a-b\left|\right|}^2}{{\mathrm{\&sigma;}}^2}\right),$ 有

$\frac{k(x,y)-k(y,z)-k(x,z)+1}{{(1-k(x,z))}^{\frac{1}{2}}{(1-k(y,z))}^{\frac{1}{2}}}=2---\left(8\right)$

方程(7)等价于方程(8)

$\frac{k(x,y)-k(y,z)-k(x,z)+1}{{(1-k(x,z)-k(y,z)+k(x,z)k(y,z))}^{\frac{1}{2}}}=2---\left(9\right)$

令

h＝1-k(y，z)-k(x，z)

$\frac{h+k(x,y)}{{(h+k(x,z)k(y,z))}^{\frac{1}{2}}}=2---\left(10\right)$

因此假如数据点x，y，z满足

$k(x,y)=2{(h+k(x,z)k(y,z))}^{\frac{1}{2}}-h---\left(11\right)$

那么θ＝0，也就是，数据点x，y，z在特征空间是线性的。

且假如

k(x，z)＝k(x，y)    (12)

k(y，z)＝4k(x，y)-3

那么满足方程(10)。

当核函数选取径向基核函数，那么当数据点x，y，z满足方程(12)时x，y，z在特征空间是线性的。因此，特征空间中相似性因子可以转化为在输入空间进行计算，定义如下：

$S_{21}=1-\frac{|k(x_{j+1},x_{k_2})-k(x_{k_1},x_{l_2})|}{k(x_{j+1},x_{k_2})}---\left(13\right)$

且

$S_{22}=1-\frac{|k(x_{j+1},x_{l_2})-4k(x_{k_2},x_{l_2})+3|}{k(x_{j+1},x_{l_2})}---\left(14\right)$

其中S₂₁，S₂₂满足0≤S₂₁≤ 1，0≤S₂₂ ≤ 1，当S₂₁＝1且S₂₂＝1时相似性是最强的。因此不需要对子集Ω₁中的数据进行非线性映射就可在特征空间对数据进行相似性分析，重新定义第j个节点集 $L_j=\{{\tilde{x}}_1,\&CenterDot;\&CenterDot;\&CenterDot;,{\tilde{x}}_{n_j}\}.$ 假如 $S_{21}<\sqrt{{\mathrm{\&gamma;}}_1}$ 且 $S_{22}<\sqrt{{\mathrm{\&gamma;}}_1},$ 其中γ₁是一个满足 $0\&le;\sqrt{{\mathrm{\&gamma;}}_1}\&le;1$ 的预定义小变量，x_j+1将被引入建立新扩展节点集L_j+1＝{L_j，x_j+1}。否则这个数据点将被拒绝，即L_j+1＝L_j。在特征空间对Ω₁中的数据进行相似性分析后，保留的数据点个数表示为N；且获得了第二个数据集Ω₂。

步骤三、利用核主元分析(KPCA)对数据做白化处理，求解白化后的观测变量z

通过非线性映射将输入空间映射到一个特征空间接着在此特征空间对观测数据进行白化处理，得到白化后的观测变量z。通过白化处理后，减少了ICA需要估计的参数，从而简化计算过程。

子集Ω₂中包含m个变量的观测数据x_k∈R^m，k＝1，...，N，其中N是观测的个数。利用非线性映射Ф：R^m→F，原始空间中的观测数据就扩展到高维特征空间F，Ф(x_k)∈F。此特征空间的协方差结构是一个单位矩阵。因而，特征空间中协方差矩阵将为

$C^F=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(x_j\right)\mathrm{\&Phi;}{\left(x_j\right)}^T---\left(15\right)$

其中Ф(x_j)，j＝1，...，N，假定为零均值和单位方差。令Θ＝[Ф(x₁)，...，Ф(x_N)]，因而C^F能表示为 $C^F=\frac{1}{N}\mathrm{\&Theta;}{\mathrm{\&Theta;}}^T.$ 定义一个N×N维的Gram核矩阵K

[K]_ij＝K_ij＝<Ф(x_i)，Ф(x_j)>＝k(x_i，x_j)    (16)

有K＝Θ^TΘ。核函数k(x_i，x_j)的应用可以在不执行非线性映射的情况下在F中计算内积。即可以避免执行非线性映射，并且在特征空间通过引入一个核函数k(x，y)＝<Ф(x)，Ф(y)>计算内积。常用的核函数有径向基核函数： $k(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{{\mathrm{\&sigma;}}^2}),$ 多项式核函数：k(x，y)＝<x，y>^r，Sigmoid(S形)核函数：k(x，y)＝tanh(β₀<x，y>+β₁)。核函数的选择潜在的决定了映射Ф以及特征空间F。从核矩阵K可知，高维空间中Ф(x_k)的中心化能按以下方法执行，即将Ф(x_k)的中心化转化为对K的中心化处理。中心化核矩阵

可由下式获得

$\tilde{K}=K-1_{N}K-K{1}_N+1_{N}K{1}_N---\left(17\right)$

其中

特征值分解

$\mathrm{\&lambda;\&alpha;}=\tilde{K}\mathrm{\&alpha;}---\left(18\right)$

式中λ为

的特征根，α为λ对应的特征向量。由式(18)能够获得同

的d个最大正特征值λ₁≥λ₂≥...≥λ_d相应的标准正交特征向量α₁，α₂，...，α_d。因而C^F的d个最大正特征值是 $\frac{{\mathrm{\&lambda;}}_1}{N},\frac{{\mathrm{\&lambda;}}_2}{N},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{{\mathrm{\&lambda;}}_d}{N},$ 相应的特征向量v₁，v₂，...，v_d能够表示为

$V_j=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_j}}\mathrm{\&Phi;}{\mathrm{\&alpha;}}_j$ j＝1，...，d    (19)

特征向量矩阵V＝[v₁，v₂，...，v_d]可以简单的由下式表达

V＝ΘHΛ^-1/2    (20)

其中Λ＝diag(λ₁，λ₂，...，λ_d)和H＝[α₁，α₂，...，α_d]分别为

的d个最大特征值的对角阵及对应的

特征向量矩阵。特征空间中的白化矩阵P和映射数据的白化变换如下

z＝P^TФ(x)    (21)

具体为，

$z=P^T\mathrm{\&Phi;}\left(x\right)=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T{\mathrm{\&Theta;}}^T\mathrm{\&Phi;}\left(x\right)=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T{[\mathrm{\&Phi;}\left(x_1\right),...,\mathrm{\&Phi;}\left(x_N\right)]}^T\mathrm{\&Phi;}\left(x\right)$

$=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T{[\tilde{k}(x_1,x),...,\tilde{k}(x_N,x)]}^T$

$=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T\tilde{k}---\left(22\right)$

式中

$\tilde{k}=k-1_{t}K-K{1}_N+1_t+1_{t}K{1}_N,$ k＝[k(x₁，x)，...，k(x_N，x)]^T，1_t＝(1/I)[1，...，1]∈R^1×N；

x₁，...，x_N—Ω₂中的数据；

x-需要白化处理的数据；

步骤四、利用修正独立元分析(ICA)提取独立元

在KPCA变换空间提取独立元s。修正ICA进一步从白化后的观测变量提取一组独立元。并使独立元各变量之间相互统计独立。找出p(≤d)个独立元，即s＝{s₁，s₂，...，s_p}，满足E(ss^T)＝D＝diag{λ₁，...，λ_p}以使s的元素变得尽可能彼此独立，利用

s＝C^Tz    (23)

其中C∈R^d×p为得分转换矩阵和C^TC＝D。

定义归一化的独立元s_n为

s_n＝D^-1/2s＝D^-1/2C^Tz＝C_n ^Tz    (24)

C_n为标准得分转换矩阵，很明显D^-1/2C^T＝C_n ^T，C_n ^TC_n＝I，以及E(s_ns_n ^T)＝I。因此作为正交的结果，把原始ICA寻找一个任意分离矩阵的问题简化为寻找更少参数估计的矩阵C_n的简单问题。因而，所提出算法的目标改为从z∈R^d寻找s_n∈R^p和C_n，使s_n的元素变得尽可能彼此独立且满足E(s_ns_n ^T)＝I。设定z的第一个数据中p个元为s_n的初始元。令C_n ^T＝[I_p0]，其中I_p是p维的单位矩阵且0是p×(d-p)的零矩阵。根据修正ICA算法计算标准得分转换矩阵C_n：

(1)选择p，估计独立元的个数，设定计数器i←1；

(2)获取初始向量c_i；

(3)令c_i←E{zg(c_i ^Tz)}-E{g′(c_i ^Tz)}c_i，其中g′是g的一阶导数，Hyvrinen给出了三种g的表达函数为g₁(u)＝tanh(a₁u)，g₂(u)＝u exp(-a₂u²/2)，g₃(u)＝u³；

(4)执行正交化： $c_i\&LeftArrow;c_i-\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}\left({c_i}^T{c}_i\right)c_i.$ 这种正交化排除了现有方案所包含的信息；

(5)归一化 $c_i\&LeftArrow;\frac{c_i}{\left|\right|c_i\left|\right|};$

(6)假如c_i还没有收敛，返回步骤(3)；

(7)假如c_i收敛，输出向量c_i。若i≤p，则设定i←i+1且返回步骤(2)；一旦找到C_n，从下式可得到独立元：

s＝D^1/2C_n ^Tz    (25)

步骤五、利用T²和SPE统计量以及LS-SVM进行故障检测与诊断

采用T²和SPE统计量进行在线故障检测，当观测数据的统计量没有超出统计量规定的控制限，则属于正常数据，反之属于异常数据，表明出现故障。一旦故障发生，将实时工况数据的独立元作为LS-SVM的输入进行故障诊断，根据训练数据建立的决策函数，即可判断故障的类型。

从Lee和Qin(2006)那，知道T²统计量和SPE统计量定义如下：

T²＝s^TD^-1s    (26)

$\mathrm{SPE}=e^{T}e={(z-\hat{z})}^T(z-\hat{z})=z^T(I-C_n{C_n}^T)z,$ z＝P^TФ(x)    (27)

其中 $e=z-\hat{z},$

能够由下式得来：

$\hat{z}{C}_n{D}^{-1/2}s=C_n{C}_n^{T}z---\left(28\right)$

因为s不服从高斯分布，T²的控制限通过F分布来决定。

SPE的控制限由下面权重χ²分布来计算(Qin，2003)

SPE～μχ_h ²，μ＝b/2a，h＝2a²/b    (29)

其中a和b分别是标准操作数中SPE的估计均值和方差。

要对故障进行分类，首先要建立各个故障的决策函数，具体如下：

LS-SVM分类器在F中构造一个最大边界(margin)超平面。在输入空间，这与个非线性决策边界相对应。

$\left\{\begin{array}{cc}\&lang;w,\mathrm{\&Phi;}\left(s_i\right)\&rang;+b\&GreaterEqual;+1-{\mathrm{\&xi;}}_i,& \mathrm{if}{y}_i=+1\\ \&lang;w,\mathrm{\&Phi;}\left(s_i\right)\&rang;+b\&le;-1+{\mathrm{\&xi;}}_i,& \mathrm{if}{y}_i=-1\end{array}\right.---\left(30\right)$

其中w表示权重向量且b为偏置(bias)。两种情况都能合成为单一框架如下：

y_i[<w，Ф(s_i>+b]≥+1-ξ_i

其中ξ_i≥0，i＝1，2，...，N是正松弛变量，用于允许不等式集的错误分类。因为两个超平面之间的边界宽是2/‖w‖，目标函数被定义为‖w‖²和惩罚项的最小值：

$\begin{array}{cc}\min & \frac{1}{2}{w}^{T}w+C\underset{i=i}{\overset{n}{\mathrm{\&Sigma;}}}\end{array},{\mathrm{\&xi;}}_i^2---\left(31\right)$

$s.t.\left\{\begin{array}{c}{y}_i[w^T\mathrm{\&Phi;}\left(s_i\right)+b]\&GreaterEqual;+1-{\mathrm{\&xi;}}_i,i=\mathrm{1,2},...,N\\ {\mathrm{\&xi;}}_i>0,i=\mathrm{1,2},...,N\\ C>0\end{array}\right.$

其中C∈R⁺是一整定超参数，它决定边界最大化和训练误差最小化之间的权衡；松弛变量ξ_i≥0，i＝1，2，...，N是用于放宽硬边界LS-SVM分类器的约束以便允许存在一些被错误分类的数据；在LV-SVM中，当引入核技巧后，上式(31)经过一定的运算后，上式最优化问题能重新描述为：

$\left[\begin{array}{cc}0& y^T\\ y& \mathrm{\&Omega;}+V_f\end{array}\right]\left[\begin{array}{c}b\\ a\end{array}\right]=\left[\begin{array}{c}0\\ I_f\end{array}\right]---\left(32\right)$

式中I_f-I_f ^T＝[1，...，1]∈R^1×N；

Ω-Ω_ij＝y_iy_jФ(s_i)^TФ(s_j)＝y_iy_jk(s_i，s_j)；

α-拉格朗日因子；

根据方程(32)可以求出偏置b和拉格朗日因子α，根据这两个参数可以确定相应的决策函数即分类器如下：

$f\left(s\right)=\mathrm{sgn}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_{i}k(s,s_i)+b)---\left(33\right)$

式中s-需要判断故障类型数据的独立元；

s_i-训练数据的独立元，i＝1，...，N；

决策函数建立后，当检测到故障，输入实时工况数据的独立元，根据对应的f(s)的值判断故障的类型。选取其中一个故障的决策函数，当f(s)＝+1时，发生本决策函数对应的故障，反之f(s)＝-1时，没有发生本决策函数对应的故障，继续根据其他决策函数进行故障诊断。

本发明的具体操作过程如图11所示的流程图，其中定义训练数据的子集Ω₂中的观测变量为x_k∈R^m，白化后的观测变量为z，独立元为s。实时工况数据为x_t∈R^m，白化后的观测变量为z_t，独立元为s_t。上面步骤中的计算公式(22)，(25)，(26)以及(33)等属于通用公式，均可用于训练数据和实时工况数据的计算。

本发明——基于改进的核独立元分析和支持向量机相结合的非线性过程故障诊断方法的优越性表现为：

提出了非线性动态过程故障诊断技术，结合Kernel、ICA和LS-SVM三者的优点，即，发挥Kernel对非线性的表达能力，同时发挥ICA对动态特性的把握能力和LS-SVM的分类能力。

附图说明

图1三个数据点在特征空间的角度；

图2田纳西-伊斯曼过程示意图；

图中：1，2，3，4，5，6，7，8，9，10，11，12，13分别表示为：流1，流2，流3，流4，流5，流6，流7，流8，流9，流10，流11，流12，流13；FC：流量控制；FI：流量指示器；PI：压力指示器；PHL：压力控制；SC：同步回旋加速器；TC：温度控制；TI：温度指示器；LI：液位指示器；LC：液位控制；XC：成分控制；XA，XB，XD，XE分别是成分A分析，成分B分析，成分D分析，成分E分析；

图3 a)，b)分别是田纳西-伊斯曼过程中故障8基于KPCA的监测结果和两个主元(nPC)；图4特征空间中的成分分析；

图5 a)，b)分别是田纳西-伊斯曼过程中故障8基于KICA的监测结果和两个独立元(IC)；

图6 a)，b)分别是田纳西-伊斯曼过程中故障8基于改进的KICA的监测结果和两个IC；

图7废水处理过程WWTP示意图；

图中：Q_in：流入量；Q_a：内部循环流量；Q_e：流出量；Q_f：进料层流量；Q_r：返回污泥的流量；Q_u：沉淀池底部流流量；Q_w：处理后的水的输入流量；

图8标准操作条件数据的变量模式；

图9 a)，b)分别是WWTP在两场暴雨情况下KICA的监测结果和KICA的两个IC；

图10 a)，b)分别是WWTP在两场暴雨情况下基于γ₀＝0.05的改进的KICA的监测结果和两个IC；

图11本发明方法实现的流程图。

以美国ROCKWELL公司的可编程控制器(PLC)实现基础控制，监控程序用RSView32提供的VBA应用软件编制。监控软件在单独的计算机上运行，该计算机上装有RSLinx通讯软件，负责与PLC和上位机进行数据通讯，RSLinx与监控程序之间通过DDE方式进行双向通讯。把监控结果输出到计算机的系统管理画面，同时把监控结果保存到实时数据库中，为操作者或相关技术工人进行监督操作提供参考指导作用。

具体实施方式

实例1、田纳西-伊斯曼(Tennessee Eastman)过程

本发明提出的方法应用到了田纳西-伊斯曼过程仿真数据中，并同KPCA及原KICA的监测结果进行对比。田纳西-伊斯曼过程是一个复杂非线性过程，它是由Eastman化学品公司所创建，其目的是为评价过程控制和监测方法提供一个真实的工业过程。控制结构如图2所示。过程包括五个主要单元：反应器，冷凝器，压缩机，分离器，汽提塔；而且，它包含了八种成分：A，B，C，D，E，F，G和H。四种反应物A，C，D以及E和惰性B一起被加进反应器里，形成产品G和H，还有副产品F。田纳西-伊斯曼过程包括21个预设定的故障，如表1所示。包括22个连续过程测量，12个控制变量，和19个成分测量。如表2所示。除了反应器的搅拌器的搅动速度(因为没有对它进行控制)，共52个观测变量用于本研究的监测。

表1.田纳西-伊斯曼过程的过程故障描述

No.	描述	类型	No.	描述	类型
						1234567891011	A/C进料比率，B成分不变(流4)B成分，A/C进料比率不变(流4)D的进料温度(流2)反应器冷却水的入口温度冷凝器冷却水的入口温度A进料损失(流1)C存在压力损失-可用性降低(流4)A，B，C进料成分(流4)D的进料温度(流2)C的进料温度(流4)反应器冷却水的入口温度	阶跃阶跃阶跃阶跃阶跃阶跃阶跃随机变量随机变量随机变量随机变量	12131415161718192021	冷凝器冷却水的入口温度反映动态反应器冷却水阀门冷凝器冷却水阀门未知未知未知未知未知流4的阀门固定在稳态位置	随机变量慢偏移粘住粘住恒定位置

表2.田纳西-伊斯曼过程中的监测变量

No.	过程测量	No.	成分测量	No.	控制变量
						12345678910111213141516171819202122	A进料(流1)D进料(流2)E进料(流3)总进料(流4)再循环流量(流8)反应器进料速度(流6)反应器压力反应器等级反应器温度排放速度(流9)产品分离器温度产品分离器液位产品分离器压力产品分离器塔底低流量(流10)汽提器等级汽提器压力汽提器塔底流量(流11)汽提器温度汽提器流量压缩机功率反应器冷却水出口温度分离器冷却水出口温度	23242526272829303132333435363738394041	成分A(流6)成分B(流6)成分C(流6)成分D(流6)成分E(流6)成分F(流6)成分A(流9)成分B(流9)成分C(流9)成分D(流9)成分E(流9)成分F(流9)成分G(流9)成分H(流9)成分D(流11)成分E(流11)成分F(流11)成分G(流11)成分H(流11)	424344454647484950515253	D进料量(流2)E进料量(流3)A进料量(流1)总进料量(流4)压缩机再循环阀排放阀(流9)分离器灌液流量(流10)汽提器液体产品流量(流11)汽提器水流阀反应器冷却水流量冷凝器冷却水流量搅拌速度

步骤一、采集数据

对于训练数据和实时工况数据采用了三分钟的采样间隔采集数据。每个故障的训练数据由480个观测数构成，实时工况数据由960个观测数构成。开始时都是没有故障的，训练数据中，故障都是在第20次采样时引入，实时工况数据中，故障都是在第160次采样引入的。训练数据和实时工况数据中的数据都包含了52个观测变量。本实例主要针对故障8进行建模分析，随机选取了故障8的训练数据和实时工况数据中控制变量的5组数据分别如表3和表4所示：

表3.故障8训练数据中控制变量的五组数据

表4.故障8实时工况数据中控制变量的五组数据

步骤二、进行相似性分析

对训练数据的数据进行输入空间的相似性分析，取γ₀＝0.95。学习的开始只有两个数据点，即节点集N₁＝{x₁，x₂}。当引入数据点x₃时，根据公式(1)进行计算，当 $S_1<\sqrt{0.95}$ 时，

x₃加入节点集N₂中，得到N₂＝{x₁，x₂，x₃}，否则N₂＝N₁。当引入数据点x₄时，也根据公式(1)进行计算，当 $S_1<\sqrt{0.95}$ 时，x₄加入节点集N₃中，得到N₃＝{N₂，x₄}，否则N₃＝N₂。

依次类推，每引入一个新的数据点，都根据公式(1)进行计算，直到引入x₄₈₀为止，得到第一个数据子集Ω₁。

再对数据子集Ω₁中的数据进行特征空间的相似性分析，取γ₁＝0.95，核函数选径向基核函数。由于Ω₁中数据的编号不再连续所以要对Ω₁中的数据重新编号。学习的开始也只有两个数据点，即节点集L₁＝{X₁，x₂}。当引入数据点x₃时，根据公式(13)、(14)进行计算，当 $S_{21}<\sqrt{0.95}$ 且 $S_{22}<\sqrt{0.95}$ 时，x₃加入节点集L₂，得到L₂＝{x₁，x₂，x₃}，否则L₂＝L₁。当引入数据点x₄时，根据公式(13)、(14)进行计算，当 $S_{21}<\sqrt{0.95}$ 且 $S_{22}<\sqrt{0.95}$ 时，x₄加入节点集L₃，得到L₃＝{L₂，x₄)，否则L₃＝L₂。依次类推，直到Ω₁中最后一个数据点进行相似性分析，得到第二个数据子集Ω₂。

步骤三、利用KPCA对数据做白化处理

选取径向基核函数作为核函数，其中σ＝0.5。对第二个数据子集Ω₂中的数据根据公式(16)进行核计算，得到Gram核矩阵K。接着根据公式(17)对Gram核矩阵K进行中心化处理，得到中心化的核矩阵

。通过计算出了它的特征值，并选取了12个最大的正特征值为λ₁≥λ₂≥...≥λ₁₂。根据公式(18)得到12个特征值对应的特征向量矩阵H。最后在特征空间进行白化，根据公式(22)依次对Ω₂中的数据进行白化处理得到白化后的观测变量z。

根据公式(22)对实时工况数据进行白化处理得到白化后的观测变量z，。

步骤四、利用修正ICA提取独立元

取p＝8，根据ICA算法计算训练数据的标准得分转换矩阵C_n，其中g(u)＝tanh(1.5u)。接着利用上一步得到白化后的观测变量z，根据公式s＝D^1/2C_n ^Tz计算得到训练数据的独立元s。再利用上一步得到的白化后的观测变量z_t，根据公式s_t＝D^1/2C_n ^Tzt计算得到实时工况数据的独立元s_t。

步骤五、故障检测与诊断

将训练数据的独立元

作为LS-SVM的输入，确定相应的决策函数。建立训练样本S：{s_i，y}，i＝1，2，…，n，其中y_i∈{+1，-1)，当s_i为故障8的数据时，y_i＝+1，否则y_i＝-1。核函数仍然选取径向基核函数。由此就建立了故障8相应的决策函数为 $f\left(s_t\right)=\mathrm{sgn}(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_{i}k(s_t,s_i)+0.1)$ 。决策函数建立后，根据公式(26)、(27)计算实时工况数据的T²和SPE统计量。从统计量图中看出检测到故障，此时将实时工况数据的独立元s_t代入决策函数进行分类，根据f(s_t)的值判断实时工况数据中发生故障的故障类型。具体监测结果见对比实例。

对比实例

对所有的数据都在KPCA、原始KICA和改进的KICA应用之前进行规范化。对于统计量，每幅图中的点线是控制限，控制限上方的统计量表示检测到一个故障(统计量用实线表示)。对于故障8，基于KPCA的监测结果如图3所示。非线性主元分析如图4所示。基于原始KICA的监测结果如图5所示。基于改进的KICA的监测结果如图6所示。在故障8的情况下，流4的A，B，C原料成分是随机改变的。在训练的过程中LS-SVM需要正常训练样本和故障训练样本。对于本次故障，KPCA能够在大约第185个样本检测到故障(图3)且延时时间是25s。然而，尽管发生故障，但是有些样本(样本210-220，370-390，630-650，750-760，790-820以及895-905)仍低于95％的控制限，这就给过程操作员带来一幅错误的过程状态图。原始KICA监测图显示统计量成功地从第180个样本检测到故障(图5)延时时间是20s。且比起低于95％的控制限的KPCA需要更少的样本。从改进的KICA监测图中可以看出，它和原始的KICA监测结果基本相同，但通过相似性分析，改进的KICA降低了计算负载。由表5可知，改进的KICA的整个监测率明显高于KPCA的。其中故障3，9和15的检测率比所有的故障检测率都要低很多。对于故障3，9，15，任何一个统计量都将导致低的检测率，换句话说，故障3，9以及15从实时工况数据中是不可观测的。

表5.田纳西-伊斯曼过程中每种方法的故障诊断率

故障	KPCA		改进的KICA
					T2	SPE	T2	SPE
	12345678910	1009849259910095443	1009881002710010097451	100989813010010097581					1009831002010010098378

11121314151617181920

249794795307490341

81989510075295904955

85999510078096918065

77999510056595857550

分类器的LS-SVM输入变量的选择对它的分类性能有着很重要的影响。输入变量是根据过程知识和操作特性进行选定的。(Chiang，2001)尽管考虑了有着不同类型和位置的三种故障，这三种故障是：故障4(反应器冷水入口温度)，故障6(A进给丢失(流1))和故障11(反应器冷水入口温度)。然而，只有变量1对故障6是重要的，只有变量9和变量51对故障4和故障11是重要的。其余的49个变量就分类而言不提供重要信息，也就是，对于故障4，6和11的训练数据和实时工况数据观测不到这49个变量有什么不同。故障4的数据集同故障11重叠。从表6和7可以看出，通过使用主元nPCs和独立元ICs，从原始信息中提取出了有用的特征。且使用ICs的分类率高于nPCs的分类率。原因是改进的KICA中的负熵能比KPCA更好地解释高阶原始输入信息。从表8中，发现改进的KICA和LS-SVM结合的训练时间比原始KICA与LS-SVM结合的更短。

表6.田纳西-伊斯曼过程中LS-SVM利用nPCs后的故障非类率

故障	训练(％)LS-SVM中的核参数	测试(％)LS-SVM中的核参数	SVs的个数
				σ＝0.5	σ＝0.5
	4611	53233				494543	333

表7.田纳西-伊斯曼过程中LS-SVM利用ICs后的故障非类率

故障

训练(％)LS-SVM中的核参数

测试(％)LS-SVM中的核参数

SVs的个数

	σ＝0.5	σ＝0.5
				4611	63537	534846	333

表8.田纳西-伊斯曼过程的训练时间

LS-SVM中的核参数	原始KICA加LS-SVM(s)	改进的KPCA加LS-SVM(s)
			σ＝0.5σ＝0.6σ＝0.7	238256289	157178195

实例2、废水处理过程(WWTP)

把基于改进的KICA的监测方法用于废水处理过程。它包括活性污泥模型No.1(ASMl)和一个十层沉淀池模型(settler model)，活性污泥模型No.1和一个十层沉淀池模型分别用于仿真生物反应和沉淀过程。WWTP系统的工艺布局如图7所示。生物反应器的第一二个分室不充气，其它的充气。所有分室都被认为是理想混合的，而二级沉淀池用一系列一维十层进行建模。对于过程监测，选定了8个变量如表9所列，因为他们都很重要且在实际WWTP系统监测中比较典型。这种过程的变量往往一个周期有很大的波动，他们的均值和方差不能够一直保持为常量。根据这种情况，常用的多变量统计过程监测(MSPM)，隐含地假定了一个稳定的潜在过程，会导致大量的错误警报及故障丢失。所以提出了本发明方法对此过程进行监测。

表9.WWTP基准的监测变量

No.	变量	描述	稳态值
				12345678	SsXbhXsXiSnhSndXndQin	易生物降解的底物活性异养生物的数量慢生物降解的底物颗粒惰性有机物NH4 ++NH3氮可溶性生物降解的有机氮颗粒可生物降解的有机氮输入流量	69.9g/m328.17g/m3202.32g/m351.2g/m331.56g/m36.95g/m310.59g/m318446m3/day

步骤一、数据采样

采集如表9所示变量的数据。包括一周的标准数据(训练数据)和两周的实时工况数据，其中采样周期为15分钟，由此获得了672组标准数据和1344组实时工况数据。在实时工况数据中，此过程经历了长期干燥的天气后突然下了两场暴雨，第一场暴雨：样本850-865，第二场暴雨：样本1050-1110。因为WWTP受到每日流量大的波动及进给流的成分的影响，数据有一定的周期特性如图8所示。分别随机选取了标准数据和实时工况数据中的十组数据如表10和表11所示：

表10.WWTP标准数据中的十组数据

表11.WWTP实时工况数据中的十组数据

步骤二、进行相似性分析

对标准数据进行输入空间的相似性分析，取γ₀＝0.95 。学习的开始只有两个数据点，即节点集N₁＝{x₁，x₂}。当引入数据点x₃时，根据公式(1)进行计算，当 $S_1<\sqrt{0.95}$ 时，x₃加入节点集N₂中，得到N₂＝{x₁，x₂，x₃}，否则N₂＝N₁。依次类推，每引入一个新的数据点，都根据公式(1)进行计算，直到引入x₆₇₂为止，得到第一个数据子集Ω₁。

再对数据子集Ω₁中的数据进行特征空间的相似性分析，取γ₁＝0.95，核函数选径向基核函数。学习的开始也只有两个数据点，即节点集L₁＝{x₁，x₂}。当引入数据点x₃时，根据公式(13)、(14)进行计算，当 $S_{21}<\sqrt{0.95}$ 且 $S_{22}<\sqrt{0.95}$ 时，x₃加入节点集L₂，得到L₂＝{x₁，x₂，x₃}，否则L₂＝L₁。依次类推，直到Ω₁中最后一个数据点进行相似性分析，得到第二个数据子集Ω₂。

步骤三、利用KPCA对数据做白化处理

选取径向基核函数作为核函数，其中σ＝0.5。对第二个数据子集Ω₂中的数据根据公式(16)进行核计算，得到Gram核矩阵K。接着根据公式(17)对Gram核矩阵K进行中心化处理，得到中心化的核矩阵通过计算出了它的特征值，并选取了4个最大的正特征值为λ₁≥λ₂≥...≥λ₄。根据公式(18)得到4个特征值对应的特征向量矩阵H。通过标准数据得到公式(22)所需参数后，就可以在特征空间利用公式(22)对实时工况数据进行白化处理，得到白化后的观测变量z_t。

步骤四、利用修正ICA提取独立元

取p＝3，根据ICA算法计算训练数据的标准得分转换矩阵C_n，其中g(u)＝tanh(1.5u)。接着利用上一步得到的白化后的观测变量z_t，根据公式s_t＝D^1/2C_n ^Tz_t计算得到实时工况数据的独立元s_t。

步骤五、故障检测

利用上一步得到的实时工况数据的独立元s_t根据公式(26)和(27)计算实时工况数据的T²统计量和SPE统计量。得到改进的KICA监测结果如图10所示，其中点线为它们的控制限，实线为统计量的分布曲线。通过实时工况数据的统计量分布，可以检测到故障发生。对于原始的KICA，得到如图9所示的监测结果。从图9和图10可以看出，对于这两场暴雨，原始的KICA能够成功的对其进行监测，改进的KICA也成功地对两场暴雨进行了监测。但是改进的KICA降低了计算负载，而且达到的监测结果与原始的KICA基本一致。由此表明改进的KICA比原始的KICA更有优势。此实例由于故障就是两场暴风雨的引入，所以不再对故障的种类进行诊断。

本发明——基于改进的核独立元分析和支持向量机相结合的非线性过程故障诊断方法还能对其他的非线性过程像连续采煤机采煤过程、轧钢过程等进行故障诊断。通过对他们的监测，及时检测出故障，避免了系统崩溃及由此造成的物质损失和人员伤亡，提高了工业过程的生产效益。

Claims (5)

1.一种非线性过程故障诊断方法，其特征在于该方法包括以下步骤：

步骤一、采集数据

采集过程中相关变量的数据，对于每个故障，产生两组数据，即训练数据和实时工况数据；训练数据用于建立模型，实时工况数据用于在线监测；并用均值和标准偏差规范化采集的数据；

步骤二、进行相似性分析

相似性分析包括输入空间的相似性分析和特征空间的相似性分析，根据相似性因子对数据进行分析，去除掉相似性较强的数据；

步骤三、利用核主元分析对数据进行白化处理，求解白化后的观测变量z

通过非线性映射将输入空间映射到一个特征空间接着在此特征空间对观测数据进行白化处理，得到白化后的观测变量z；

步骤四、利用修正独立元分析ICA提取独立元

在核主元分析KPCA变换空间提取独立元s；利用修正ICA方法从白化后的观测变量z中提取一组独立元；并使独立元各变量之间相互统计独立；

步骤五、利用T²和SPE统计量和LS-SVM进行故障检测与诊断

采用T²和SPE统计量进行在线故障检测，当观测数据的统计量没有超出统计量规定的控制限，则属于正常数据，反之属于异常数据，表明出现故障；一旦故障发生，将实时工况数据的独立元作为LS-SVM的输入进行故障诊断，根据训练数据建立的决策函数，即可判断故障的类型。

2.根据权利要求1所述的一种非线性过程故障诊断方法，其特征在于所述的步骤二中的相似性分析的具体过程如下：

1)输入空间的相似性分析

对采集的数据进行相似性分析；

学习的开始，只有两个数据点，即 $N_1=\{{\tilde{x}}_1,{\tilde{x}}_2\}=\{x_1,x_2\};$ 假设第i个节点集被表示为 $N_i=\{{\tilde{x}}_1,\&CenterDot;\&CenterDot;\&CenterDot;,{\tilde{x}}_{n_i}\},$ i＝1，...，t₁，且满足n_i＜t₁，t₁为训练数据的个数；定义相似性因子为

$S_1={\left(\frac{\mathrm{cov}(x_{i+1}-{\tilde{x}}_{k_1},x_{i+1}-{\tilde{x}}_{l_1})}{\mathrm{std}(x_{i+1}-{\tilde{x}}_{k_1})\mathrm{std}(x_{i+1}-{\tilde{x}}_{l_1})}\right)}^2---\left(2\right)$

式中x_i+1-训练数据中的第i+1个数据点；

-节点集N_i中任意两个数据点；

S₁满足0≤S₁≤1；当S₁＝1时相似性最强；将x_i+1同节点集N_i中所有数据都两两进行相似性计算，假如所有 $S_1<\sqrt{{\mathrm{\&gamma;}}_0},$ 其中γ₀是满足 $0\&le;\sqrt{{\mathrm{\&gamma;}}_0}\&le;1$ 的预定义小变量，新数据x_i+1将被引入建立新扩展节点集N_i+1＝{N_i，x_i+1}；否则，这个数据点将被拒绝，即N_i+1＝N_i；在输入空间对所有训练数据进行相似性分析后，获得第一个数据子集Ω₁；

2)特征空间的相似性分析

得到子集Ω₁后，将Ω₁中的数据通过非线性函数Ф映射到特征空间，进行特征空间的相似性分析；学习的开始，只有两个数据，即 $L_1=\{\mathrm{\&Phi;}\left({\tilde{x}}_1\right),\mathrm{\&Phi;}\left({\tilde{x}}_2\right)\}=\{\mathrm{\&Phi;}\left(x_1\right),\mathrm{\&Phi;}\left(x_2\right)\};$ 假设特征空间第j个节点集被表示为 $L_j=\{\mathrm{\&Phi;}\left({\tilde{x}}_1\right),...,\mathrm{\&Phi;}\left({\tilde{x}}_{n_j}\right)\},$ j＝1，...，t₂，且满足n_j＜t₂，t₂为子

集Ω₁中数据的个数；定义特征空间的相似性因子为：

$S_2={\left(\frac{\mathrm{cov}(\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{k_2}\right),\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{l_2}\right))}{\mathrm{std}(\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{k_2}\right))\mathrm{std}(\mathrm{\&Phi;}\left(x_{j+1}\right)-\widetilde{\mathrm{\&Phi;}}\left(x_{l_2}\right))}\right)}^2$

式中Ф(x_j+1)-数据子集Ω₁中第j+1个数据点x_j+1映射后的数据；

-节点集L_j中的任意两个数据；

由于非线性函数Ф难以直接确定，利用核技巧k(a，b)＝<Ф(a)，Ф(b)>，其中a，b为输入空间的数据，对此问题进行了解决，核函数采用径向基核函数 $k(x,y)=\exp (-\frac{{\left|\right|a-b\left|\right|}^2}{{\mathrm{\&sigma;}}^2});$

经过一定转换，特征空间中相似性因子转化为在输入空间进行计算，定义如下：

$S_{21}=1-\frac{|k(x_{j+1},x_{k_2})-k(x_{k_1},x_{l_2})|}{k(x_{j+1},x_{k_2})}---\left(13\right)$

且 $S_{22}=1-\frac{|k(x_{j+1},x_{l_2})-4k(x_{k_2},x_{l_2})+3|}{k(x_{j+1},x_{l_2})}---\left(14\right)$

其中S₂₁，S₂₂满足0≤S₂₁≤1，0≤S₂₂≤1，当S₂₁＝1且S₂₂＝1时相似性最强；因此不需要对子集Ω₁中的数据进行非线性映射就可在特征空间对数据进行相似性分析，重新定义第j个节点集 $L_j=\{{\tilde{x}}_1,\&CenterDot;\&CenterDot;\&CenterDot;,{\tilde{x}}_{n_j}\};$ 假如 $S_{21}<\sqrt{{\mathrm{\&gamma;}}_1}$ 且 $S_{22}<\sqrt{{\mathrm{\&gamma;}}_1},$ 其中γ₁是一个满足 $0\&le;\sqrt{{\mathrm{\&gamma;}}_1}\&le;1$ 的预定义小变量，x_j+1将被引入建立新扩展节点集L_j+1＝{L_j，x_j+1}；否则这个数据点将被拒绝，

即L_j+1＝L_j；在特征空间对Ω₁中的数据进行相似性分析后，保留的数据点个数表示为N；且获得了第二个数据集Ω₂。

3.根据权利要求1所述的一种非线性过程故障诊断方法，其特征在于所述的步骤三中求解白化后的观测变量z的具体过程如下：

子集Ω₂中包含m个变量的观测数据x_k∈R^m，k＝1，...，N，其中N是观测的个数；利用非线性映射Ф：R^m→F，原始空间中的观测数据就扩展到高维特征空间F，Ф(x_k)∈F；

此特征空间的协方差结构是一个单位矩阵；因而，特征空间中协方差矩阵将为

$C^F=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(x_j\right)\mathrm{\&Phi;}{\left(x_j\right)}^T---\left(15\right)$

其中Ф(x_j)，j＝1，...，N，假定为零均值和单位方差；令Θ＝[Ф(x₁)，...，Ф(x_N)]，因而C^F能表示为 $C^F=\frac{1}{N}\mathrm{\&Theta;}{\mathrm{\&Theta;}}^T;$ 定义一个N×N维的Gram核矩阵K

[K]_ij＝K_ij＝<Ф(x_i)，Ф(x_j)>＝k(x_i，x_j)    (16)

有K＝Θ^TΘ；核函数k(x_i，x_j)的应用可以在不执行非线性映射的情况下在F中计算内积；即可以避免执行非线性映射，并且在特征空间通过引入一个核函数k(x，y)＝<Ф(x)，Ф(y)>计算内积；从核矩阵K可知，高维空间中Ф(x_k)的规范化能按以下方法执行，即将Ф(x_k)的中心化转化为对K的中心化处理；中心化核矩阵

可有下式获得

$\tilde{K}=K-1_{N}K-K{1}_N+1_{N}K{1}_N---\left(17\right)$

其中

把

特征值分解

$\mathrm{\&lambda;\&alpha;}=\tilde{K}\mathrm{\&alpha;}---\left(18\right)$

式中λ为

的特征根，α为λ对应的特征向量；由式(18)能够获得同

的d个最大正特征值λ₁≥λ₂≥...≥λ_d相应的标准正交特征向量α₁，α₂，...，α_d；因而C^F的d个最大正特征值是 $\frac{{\mathrm{\&lambda;}}_1}{N},\frac{{\mathrm{\&lambda;}}_2}{N},...,\frac{{\mathrm{\&lambda;}}_d}{N},$ 相应的特征向量v₁，v₂，...，v_d能够表示为

$v_j=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_j}}\mathrm{\&Theta;}{\mathrm{\&alpha;}}_j$ j＝1，...，d    (19)

特征向量矩阵V＝[v₁，v₂，...，v_d]可以简单的由下式表达

V＝ΘHΛ^-1/2    (20)

其中Λ＝diag(λ₁，λ₂，...，λ_d)和H＝[α₁，α₂，...，α_d]分别为

的d个最大特征值的对角阵及对应的特征向量；特征空间中的白化矩阵P和映射数据的白化变换如下

z＝P^TФ(x)    (21)

具体为，

$z=P^T\mathrm{\&Phi;}\left(x\right)=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T{\mathrm{\&Theta;}}^T\mathrm{\&Phi;}\left(x\right)=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T{[\mathrm{\&Phi;}\left(x_1\right),\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Phi;}\left(x_N\right)]}^T\mathrm{\&Phi;}\left(x\right)$

$=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T{[\tilde{k}(x_1,x),\&CenterDot;\&CenterDot;\&CenterDot;,\tilde{k}(x_N,x)]}^T---\left(22\right)$

$=\sqrt{N}{\mathrm{\&Lambda;}}^{-1}{H}^T\tilde{k}$

式中

$\tilde{k}=k-1_{t}K-K{1}_N+1_t+1_{t}K{1}_N,$ k＝[k(x₁，x)，...，k(x_N，x)]^T，1_t＝(1/I)[1，...，1]R^1×N；

x₁，...，x_N，—Ω₂中的数据；

x-需要白化处理的数据。

4.根据权利要求1所述的一种非线性过程故障诊断方法，其特征在于所述的步骤四求解独立元s的具体过程如下：

找出p(≤d)个独立元，即s＝{s₁，s₂，...，s_p}，满足E(ss^T)＝D＝diag{λ₁，...，λ_p}以使s的元素变得尽可能彼此独立，利用

s＝C^Tz    (23)

其中C∈R^d×p为得分转换矩阵和C^TC＝D；定义归一化的独立元为

s_n＝D^-1/2s＝D^-1/2C^Tz＝C_n ^Tz    (24)

C_n为标准得分转换矩阵，很明显D^-1/2C^T＝C_n ^T，C_n ^TC_n＝I，以及E(s_ns_n ^T)＝I；因此算法的目标改为从z∈R^d寻找s_n∈R^p和C_n，使s_n的元素变得尽可能彼此独立且满足E(s_ns_n ^T)＝I；设定z的第一个数据中p个元为s_n的初始元；令C_n ^T＝[I_p0]，其中I_p是p维的单位矩阵且0是p×(d-p)的零矩阵；根据修正ICA算法计算标准得分转换矩阵C_n：

(1)选择p，估计独立元的个数，设定计数器i←1；

(2)获取初始向量c_i；

(3)令c_i←E{zg(c_i ^Tz)}-E{g′(c_i ^Tz)}c_i，其中g′是g的一阶导数，Hyvrinen给出了三种g的表达函数为g₁(u)＝tanh(a₁u)，g₂(u)＝uexp(-a₂u²/2)，g₃(u)＝u³；

(4)执行正交化： $c_i\&LeftArrow;c_i-\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}\left({c_i}^T{c}_i\right)c_i;$ 这种正交化排除了现有方案所包含的信息；

(5)归一化 $c_i\&LeftArrow;\frac{c_i}{\left|\right|c_i\left|\right|};$

(6)假如c_i还没有收敛，返回步骤(3)；

(7)假如c_i收敛，输出向量c_i，若i≤p，则设定i←i+1且返回步骤(2)；

一旦找到C_n，根据下式得到独立元

s＝D^1/2C_n ^Tz    (25)。

5.根据权利要求1所述的一种非线性过程故障诊断方法，其特征在于所述的步骤五中检测统计量的计算以及决策函数的建立过程如下：

T²统计量和SPE统计量定义如下：

T²＝s^TD^-1s    (26)

$\mathrm{SPE}=e^{T}e={(z-\hat{z})}^T(z-\hat{z})=z^T(I-C_n{C_n}^T)z,$ z＝P^TФ(x)    (27)

其中 $e=z-\hat{z},$

能够由下式得来：

$\hat{z}{C}_n{D}^{-1/2}s=C_n{C}_n^{T}z---\left(28\right)$

因为s不服从高斯分布，T²的控制限通过F分布来决定；SPE的控制限由下面权重χ²分布来计算(Qin，2003)

SPE～μχ_h ²，μ＝b/2a，h＝2a²/b    (29)

其中a和b分别是标准操作数中SPE的估计均值和方差；

要对故障进行分类，首先要建立各个故障的决策函数，具体如下：

LS-SVM分类器在F中构造一个最大边界(margin)超平面；在输入空间，这与一个非线性决策边界相对应；

$\left\{\begin{array}{cc}\&lang;w,\mathrm{\&Phi;}\left(s_i\right)\&rang;+b\&GreaterEqual;+1-{\mathrm{\&xi;}}_i,& \mathrm{if}{y}_i=+1\\ \&lang;w,\mathrm{\&Phi;}\left(s_i\right)\&rang;+b\&le;-1+{\mathrm{\&xi;}}_i,& \mathrm{if}{y}_i=-1\end{array}\right.---\left(30\right)$

其中w表示权重向量且b为偏置(bias)；两种情况都能合成为单一框架如下：

y_i[<w，Ф(s_i)>+b]≥+1-ξ_i

其中ξ_i≥0，i＝1，2，...，N是正松弛变量，用于允许不等式集的错误分类；因为两个超平面之间的边界宽是2/‖w‖，目标函数被定义为‖w‖²和惩罚项的最小值：

$\begin{array}{cc}\min & \frac{1}{2}{w}^{T}w+c\underset{j=i}{\overset{n}{\mathrm{\&Sigma;}}}\end{array},{\mathrm{\&xi;}}_i^2---\left(31\right)$

$s.t.\left\{\begin{array}{c}{y}_i[w^T\mathrm{\&Phi;}\left(s_i\right)+b]\&GreaterEqual;+1-{\mathrm{\&xi;}}_i,i=\mathrm{1,2},...,N\\ {\mathrm{\&xi;}}_i>0,i=\mathrm{1,2},...,N\\ C>0\end{array}\right.$

其中C∈R⁺是一整定超参数，它决定边界最大化和训练误差最小化之间的权衡；松弛变量ξ_i≥0，i＝1，2，...，N是用于放宽硬边界LS-SVM分类器的约束以便允许存在一些被错误分类的数据；在LV-SVM中，当引入核技巧后，方程(31)经过一定的运算后，能重新描述为：

$\left[\begin{array}{cc}0& y^T\\ y& \mathrm{\&Omega;}+V_f\end{array}\right]\left[\begin{array}{c}b\\ a\end{array}\right]=\left[\begin{array}{c}0\\ I_f\end{array}\right]---\left(32\right)$

式中I_f-I_f ^T＝[1，...，1]∈R^1×N；

Ω-Ω_ij＝y_iy_jФ(s_i)^TФ(s_j)＝y_iy_jk(s_i，s_j)；

α-拉格朗日因子；

根据方程(32)求出偏置b和拉格朗日因子α，根据这两个参数确定相应的决策函数即分类器如下：

$f\left(s\right)=\mathrm{sgn}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_{i}k(s,s_i)+b)---\left(33\right)$

式中s-需要判断故障类型数据的独立元；

s_i-训练数据的独立元，i＝1，...，N；

决策函数建立后，当检测到故障，输入实时工况数据的独立元，根据f(s)的值判断故障的类型，选取其中一个故障的决策函数，当f(s)＝+1时，发生本决策函数对应的故障，反之f(s)＝-1时，没有发生本决策函数对应的故障，继续根据其他决策函数进行故障诊断。

Images