基于核主元分析贡献图的非线性过程故障辨识方法

摘要

一种基于核主元分析贡献图的非线性故障辨识方法，包括采集数据、提取特征空间的基、提取非线性主元、故障检测和故障辨识五个步骤，本发明提出了一种新的特征空间基的提取方法，通过基的提取来减少了数据映射到特征空间后的线性冗余，以及在训练样本数量很大的时候减少KPCA的计算量。本发明方法采用贡献图进行故障辨识，对过程的每个变量计算其贡献图和控制限，通过贡献图和控制限的关系判断出各变量失控对故障的发生所承担的责任，此辨识方法克服了输入空间和特征空间不能自由转换所带来的故障辨识的困难。

说明书

技术领域

本发明属于故障检测与诊断技术领域，提出了一种基于核主元分析贡献图的非线性过程 故障辨识方法。

背景技术

近十年里，多变量统计方法在包括故障检测和故障诊断的过程监控中已经得到了迅速的 发展。如主元分析(PCA)，偏最小二乘(PLS)，独立元分析(ICA)和支持向量基(SVM) 等已经在工业过程中得到了广泛的应用。但是这些方法主要解决线性变化问题，基于这些方 法的过程监控一般假定过程为线性的，而实际的工业过程都存在不同程度的非线性变化。对 于非线性变化较强的过程，应用这些方法进行在线监控将产生较高的误报率。为了解决观测 数据的非线性问题，人们就提出了一系列的非线性方法。如一种基于主曲线和神经网络的非 线性PCA、一种基于自联想的五层神经网络的非线性PCA以及基于输入训练神经网络的非线 性PCA方法等。但是大部分现有的非线性PCA方法都是基于神经网络，必须离线训练，要 计算主元的话也必须解决非线性最优化问题。

最近，由Schlkopf等人首次在《Neural Computation》提出一种新的PCA的非线性扩展， 即核主元分析(KPCA)。核主元分析是首次通过非线性函数把输入空间映射成特征空间，接 着在特征空间中提取主元。它通过引入一个核函数避免执行非线性映射和在特征空间中计算 内积。同其他非线性方法相比，它主要的优点是不涉及非线性最优化问题，本质上仅需要线 性代数，使用起来同PCA一样简单。但是在KPCA的训练过程中核矩阵的大小是样本数量的 平方，当样本数变大，特征值和特征向量的计算就变的很费时间。而且对于现有的KPCA，故 障变量的分离和故障的识别仍然是个问题。在PCA方法中，故障变量和监控量之间存在一定 的线性关系，可以比较容易的计算变量贡献用于绘制贡献图。如Qin在《JOURNAL OF CHEMOMETRICS》中所提出的在输入空间基于重构的贡献图方法。而对于KPCA，由于输 入空间到特征空间的非线性映射函数未知，因此类似PCA中的贡献图无法直接用于KPCA的 故障变量辨识。

发明内容

针对上述现有技术存在的不足，本发明提供了一种基于改进的核主元分析的非线性过程 故障辨识方法。它在原始KPCA的基础上进行了以下两方面的改进：

1)提出了特征空间基的概念，通过提取基来减少了数据映射到特征空间后的线性冗余， 以及在训练样本数量很大的时候减少KPCA的计算量；

2)提出了基于重构的故障辨识策略，首次通过基于重构的贡献图来辨识过程中各变量失 控对故障产生所承担的责任。

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

步骤一、采集数据

采集过程相关变量的数据，并对其进行中心化和归一化，即规范化处理。

步骤二、提取特征空间的基

根据特征因子依次判断规范化后的数据与基中数据的线性关系，并去除线性关系较强的 数据，其余数据加入基中，最后得到基Ω。为了获得一个能够代表特征空间所有信息的特征 空间的基，本发明采用基于几何学思想的方法。假定${\Phi }_B=(\Phi \left({\tilde{x}}_1\right),\Phi \left({\tilde{x}}_2\right),···,\Phi \left({\tilde{x}}_{B_S}\right))$是一个特 征空间的基，其中Φ(·)表示输入空间映射到特征空间的非线性映射函数。特征空间的任意数 据点能够表示如下：

$\hat{\Phi }\left(x_i\right)={\Phi }_B{\beta }_i---\left(1\right)$

其中${\beta }_i={({\beta }_i^1,{\beta }_i^2,···,{\beta }_i^{B_S})}^T$是系数向量。本发明通过判断数据点Φ(x_i)与差值的大小来寻 找一个特征空间的基，差值表示为

${\eta }_i={\left|\right|\Phi \left(x_i\right)-\hat{\Phi }\left(x_i\right)\left|\right|}^2---\left(2\right)$

当差值η_i等于0时，则$\Phi \left(x_i\right)=\hat{\Phi }\left(x_i\right)={\Phi }_B·{\beta }_i,$表明数据点Φ(x_i)可由基Φ_B线性表示。

通过对方程(2)的最小化求解得到

学习的开始，基中只有一个数据点，即${\Phi }_1=\left\{\Phi \left({\tilde{x}}_1\right)\right\}=\left\{\Phi \left(x_1\right)\right\}.$第i个特征空间的基定义 为${\Phi }_i=\{\Phi \left({\tilde{x}}_1\right),···,\Phi \left({\tilde{x}}_{n_i}\right)\}=\{\Phi \left(x_1\right),···,\Phi \left(x_{n_i}\right)\},$i＝1，2，…，t₁，且满足n_i＜t₁，t₁为采集数据的个 数。根据方程(3)得到

${\beta }_{i+1}={\left({\Phi }_i^T{\Phi }_i\right)}^{-1}{\Phi }_i^T\Phi \left(x_{i+1}\right)---\left(4\right)$

根据方程(1)和(2)定义特征因子

$S_{i+1}=\frac{{\left|\right|\Phi \left(x_{i+1}\right)-\hat{\Phi }\left(x_{i+1}\right)\left|\right|}^2}{{\left|\right|\Phi \left(x_{i+1}\right)\left|\right|}^2}---\left(5\right)$

其中$\hat{\Phi }\left(x_{i+1}\right)={\Phi }_i{\beta }_{i+1},$当S_i+1＝0时，Φ(x_i+1)可由基Φ_i中数据点线性表示。

本发明引入核技巧k(x，y)＝(Φ(x)，Φ(y))，从而特征空间的特征提取转化为在输入空间进 行计算，方程(4)和方程(5)分别转化为

$S_{i+1}=\frac{k(x_{i+1},x_{i+1})-\left(\begin{array}{ccc}k(x_{i+1},x_1)& ···& k(x_{i+1},x_{n_i})\end{array}\right){\beta }_{i+1}}{k(x_{i+1},x_{i+1})}---\left(7\right)$

因此，通过方程(6)和(7)，即可在输入空间确定最终的基。学习的开始基中只有一个数 据即${\Phi }_1=\left\{{\tilde{x}}_1\right\}=\left\{x_1\right\}.$当$S_{i+1}>\sqrt{{\gamma }_0},$其中γ₀是满足$0\le \sqrt{{\gamma }_0}\le 1$且预先确定的较小常数，新数 据x_i+1将被引入建立新基，得到新基Φ_i+1＝{Φ_i，x_i+1}；否则，这个数据点将被拒绝，没有基扩 展，即Φ_i+1＝Φ_i。按此方法依次对采集的数据进行特征提取，去除与基中数据线性相关性较 强的数据，最后得到基Ω，它保留了N个数据。

步骤三、提取非线性主元

KPCA把标准PCA扩展到了非线性数据分布。KPCA的基本思想是通过非线性映射把输 入空间的观测数据x映射到特征空间F，接着在特征空间F执行线性PCA。在执行PCA之前， 将Ω中的观测数据x_k∈R^m，k＝1，2，…，N，映射到一个高维特征空间F

x_k→Φ(x_k)                       (8)

假定${\Sigma }_{j=1}^N\Phi \left(x_j\right)=0,$则映射后数据的协方差矩阵可表示为

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

对协方差进行特征分解

λv＝C^Fv                        (10)

其中特征值λ≥0，v∈F≠{0}是λ对应的特征向量。这里，C^Fv能够表示如下：

$C^{F}v=\left(\frac{1}{N}\underset{j=1}{\overset{N}{\Sigma }}\Phi \left(x_j\right){\Phi \left(x_j\right)}^T\right)v=\frac{1}{N}\underset{j=1}{\overset{N}{\Sigma }}⟨\Phi \left(x_j\right),v⟩\Phi \left(x_j\right)---\left(11\right)$

其中<x，y>表示x和y之间的内积。因此λv＝C^Fv等价于

λ<Φ(x_k)，v>＝<Φ(x_k)，C^Fv>，k＝1，2，…，N    (12)

且存在系数α_i(i＝1，2，…，N)使得

$v=\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i\Phi \left(x_i\right)---\left(13\right)$

结合方程(12)和(13)，得到

$\lambda \underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i⟨\Phi \left(x_k\right),\Phi \left(x_i\right)⟩=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i⟨\Phi \left(x_k\right),\underset{j=1}{\overset{N}{\Sigma }}\Phi \left(x_i\right)⟩⟨\Phi \left(x_j\right),\Phi \left(x_i\right)⟩---\left(14\right)$

其中所有k＝1，2，…，N。通常，尽管映射函数Φ(·)存在，但计算上并不容易处理。然而，它不 需要明确被计算，仅仅需要知道特征空间中两个向量的内积。通过引入核函数 k(x，y)＝<Φ(x)，Φ(y)>避免了执行非线性映射以及在特征空间计算两者的内积。

现在，定义一个N×N维的核矩阵K

[K]_ij＝K_ij＝(Φ(x_i)，Φ(x_j))＝k(x_i，x_j)  i＝1，2，…，N，j＝1，2，…，N    (15) 那么方程(14)的左边能够表示为

$\lambda \underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i⟨\Phi \left(x_k\right),\Phi \left(x_i\right)⟩=\lambda \underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i{K}_{\mathrm{ki}}---\left(16\right)$

由于k＝1，2，…，N，方程(16)转化为λKα。方程(14)的右边能够表示为

$\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i⟨\Phi \left(x_k\right),\underset{j=1}{\overset{K}{\Sigma }}\Phi \left(x_i\right)⟩⟨\Phi \left(x_j\right),\Phi \left(x_i\right)⟩=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i\underset{j=1}{\overset{N}{\Sigma }}{K}_{\mathrm{kj}}{K}_{\mathrm{ji}}---\left(17\right)$

由于k＝1，2，…，N，方程(17)转化为(1/N)K²α。结合方程(16)和(17)，得到

λNKα＝K²α                                (18)

其中系数α＝[α₁，…，α_N]^T。方程(18)的求解转化为下式特征向量的求解

Nλα＝Kα                                  (19)

在应用KPCA之前，应先执行高维空间中的中心化。核矩阵K中心化如下

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

其中

因此方程(19)转化为

$\mathrm{N\lambda \alpha }=\tilde{K}\alpha ---\left(21\right)$

现在，在特征空间执行PCA等价于求解方程(21)。求解后得到特征值λ₁≥λ₂≥…≥λ_N对应 的特征向量α₁，α₂，…α_N。通过仅保留前p个特征向量来降低问题的维数。按特征空间F中相 应向量归一化的要求对α₁，α₂，…α_p进行归一化，即，

<v_k，v_k>＝1，k＝1，2，…，p              (22)

利用$v_k=\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i^k\Phi \left(x_i\right),$方程(22)改为

$1=⟨\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i^k\Phi \left(x_i\right),\underset{j=1}{\overset{N}{\Sigma }}{\alpha }_j^k\Phi \left(x_j\right)⟩$

$=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\alpha }_i^k{\alpha }_j^k⟨\Phi \left(x_i\right),\Phi \left(x_j\right)⟩$

$=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\alpha }_i^k{\alpha }_j^k{K}_{\mathrm{ij}}=⟨{\alpha }_k,{\mathrm{K\alpha }}_k⟩={\lambda }_k⟨{\alpha }_k,{\alpha }_k⟩---\left(23\right)$

式中α_i^k，α_j^k——分别表示第k个特征向量α_k的第i行和第j行的值。 通过把Φ(x)投影到特征空间F中的特征向量v_k上，其中k＝1，2，…，p，任意数据x的非线性 主元t计算如下

$t_k=⟨v_k,\tilde{\Phi }\left(x\right)⟩=\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i^k⟨\tilde{\Phi }\left(x_i\right),\tilde{\Phi }\left(x\right)⟩=\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i^k\tilde{k}(x_i,x)---\left(24\right)$

式中k-k＝1，2，…，p；

x-要进行非线性主元提取的数据；

x_i-训练数据中的第i个数据，i＝1，2，…，N；

是中心化核向量的第i个值，

其中$\tilde{k}=k-1_{t}K-{k1}_N+1_t{K1}_N,$k＝[k(x，x₁)，…，k(x，x_N)]。

步骤四、故障检测

利用Hotelling的T²统计和平方预测误差(SPE)统计进行故障检测，通过控制图即可判 断出是否发生故障，当T²统计和SPE统计超出各自的控制限时，认为有故障发生，反之，整 个过程正常。

Hotelling的T²统计量和平方预测误差(SPE)统计量的计算如下

T²＝t^TΛ^-1t＝[t₁，t₂，…，t_p]Λ^-1[t₁，t₂…，t_p]^T    (25)

其中t_k(k＝1，2，…，p)由方程(24)决定，Λ^-1是训练数据前p个最大正特征值的对角阵的逆。

T²的100α％控制限是利用F分布获得的

$T_{\lim }^2=\frac{p(n-1)}{n-p}F(p,n-1,\alpha )---\left(26\right)$

其中F(p，n-1，α)是自由度为p和n-p的F分布的上限100α％的临界点。

$\mathrm{SPE}=\tilde{K}(x,x)-t^{T}t---\left(27\right)$

SPE的100(1-α)％控制限是

${{\mathrm{SPE}}_{\lim }=\theta }_1{(\frac{c_{\alpha }\sqrt{{2\theta }_2{h}_0^2}}{{\theta }_1}+1+\frac{{\theta }_2{h}_0(h_0-1)}{{\theta }_1^2})}^{1/h_0}---\left(28\right)$

其中

${\theta }_i=\underset{j=i+1}{\overset{m}{\Sigma }}{\lambda }_j^i,i=\mathrm{1,2,3}---\left(29\right)$

$h_0=1-\frac{{2\theta }_1{\theta }_3}{{3\theta }_2^2}---\left(30\right)$

且c_α是上限100(1-α)％相应的正常偏差。m表示特征空间有效的维数。

步骤五、故障辨识

本发明在特征空间提出了一种KPCA技术的故障识别方法，即贡献图方法。一旦检测到 故障发生，采用贡献图进行故障辨识，对过程的每个变量计算其贡献图和控制限，通过贡献 图和控制限的关系即可判断出各变量失控对故障的发生所承担的责任。

按照每个已知故障的方向重构故障。首次提出基于重构的非线性贡献图。当Ξ_i被假定为 输入空间的故障方向，Φ(Ξ_i)是特征空间的故障方向，他们并不一定是真实的故障，可能是 也可能不是，因而，按照Φ(Ξ_i)方向定义重构的测量向量如下

$z_i=\bar{\Phi }\left(x\right)-\Phi \left({\Xi }_i\right)f_i---\left(31\right)$

其中均值和修正系数f_i为

$\bar{\Phi }\left(x\right)=\Phi \left(x\right)+\underset{i=1}{\overset{N}{\Sigma }}\Phi \left(x_i\right)---\left(32\right)$

$f_i=\arg \min \mathrm{SPE}\left(z_i\right)=\arg \min {\left|\right|\tilde{\Phi }\left(x\right)-\tilde{\Phi }\left({\Xi }_i\right)f_i\left|\right|}^2---\left(33\right)$

其中

$\tilde{\Phi }\left(x\right)=(I-{\mathrm{vv}}^T)\bar{\Phi }\left(x\right)---\left(34\right)$

$\tilde{\Phi }\left({\Xi }_i\right)=(I-{\mathrm{vv}}^T)\Phi \left({\Xi }_i\right)---\left(35\right)$

式中$v=\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right)$为训练数据协方差矩阵的特征向量，α_j为训练数据的中心化核矩阵的特 征向量。

因此

$f_i={\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left(x\right)---\left(36\right)$

且重构的残差向量变为

${\tilde{z}}_i=[I-\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)]\tilde{\Phi }\left(x\right)---\left(37\right)$

故障重构的残差向量根据故障自由残差测量可以表示为

${\tilde{z}}_i=[I-\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)]{\tilde{\Phi }}^{*}\left(x\right)+$

$[I-\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)]\tilde{\Phi }\left({\Xi }_j\right)f_i---\left(38\right)$

其中是假定的实际故障，即$\tilde{\Phi }\left({\Xi }_i\right)=\tilde{\Phi }\left({\Xi }_j\right),$那么

${\tilde{z}}_i=[I-\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)]{\tilde{\Phi }}^{*}\left(x\right)---\left(39\right)$

和

$\left|\right|{\tilde{z}}_i\left|\right|=\left|\right|[I-\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)]{\tilde{\Phi }}^{*}\left(x\right)\left|\right|\le \left|\right|{\tilde{\Phi }}^{*}\left(x\right)\left|\right|---\left(40\right)$

因此SPE指标就被引入了控制限δ²，得到

$\mathrm{SPE}\left|\right|z_i\left|\right|={\left|\right|{\tilde{z}}_i\left|\right|}^2\le {\left|\right|{\tilde{\Phi }}^{*}\left(x\right)\left|\right|}^2=\mathrm{SPE}\left({\Phi }^{*}\left(x\right)\right)\le {\delta }^2---\left(41\right)$

从方程(31)中，残余子空间的大量修正是因而，通过SPE的重构定义贡献图为

$\mathrm{SPE}\_t={\left|\right|\tilde{\phi }\left({\Xi }_i\right)f_i\left|\right|}^2=f_i^2{\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)---\left(42\right)$

方程(42)取代了方程(36)，得到

$\mathrm{SPE}\_t={\left[{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left(x\right)\right]}^2{\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)$

令$\tilde{C}=I-{\mathrm{vv}}^T$因为${\tilde{C}}^2=\tilde{C},$${\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left(x\right)={\tilde{\Phi }}^T\left({\Xi }_i\right)\Phi \left(x\right).$因此，

$\mathrm{SPE}\_t={\Phi }^T\left(x\right)\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\Phi \left(x\right)---\left(43\right)$

其中

${\Phi }^T\left(x\right)\tilde{\Phi }\left({\Xi }_i\right)={\Phi }^T\left(x\right)(I-{\mathrm{vv}}^T)\Phi \left({\Xi }_i\right)={\Phi }^T\left(x\right)(I-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right){\left(\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right)\right)}^T)\Phi \left({\Xi }_i\right)$

$={\Phi }^T\left(x\right)(I-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right){\left(\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right)\right)}^T)\Phi \left({\Xi }_i\right)$

$=k(x,{\Xi }_i)-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_{j}k(x,x_j)\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i)---\left(44\right)$

${\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)={\left((I-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right){\left(\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_l\Phi \left(x_l\right)\right)}^T)\Phi \left({\Xi }_i\right)\right)}^T(I-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right){\left(\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_l\Phi \left(x_l\right)\right)}^T)\Phi \left({\Xi }_i\right)$

$=({\tilde{\Phi }}^T\left({\Xi }_i\right)-\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i){\left(\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right)\right)}^T)(\tilde{\Phi }\left({\Xi }_i\right)-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right)\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i))$

$=1-2\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i)\left(\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_{j}k(x_j,{\Xi }_i)\right)+$

$\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i)\underset{j=1}{\overset{P}{\Sigma }}\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_j{\alpha }_{l}k(x_l,x_j)\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i)---\left(45\right)$

${\tilde{\Phi }}^T\left({\Xi }_i\right)\Phi \left(x\right)={\left((I-{\mathrm{vv}}^T)\Phi \left({\Xi }_i\right)\right)}^T\Phi \left(x\right)$

$={(\Phi \left({\Xi }_i\right)-\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_j\Phi \left(x_j\right){\left(\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_l\Phi \left(x_l\right)\right)}^T\Phi \left({\Xi }_i\right))}^T\Phi \left(x\right)$

$=k(x,{\Xi }_i)-\underset{l=1}{\overset{P}{\Sigma }}{\alpha }_{l}k(x_l,{\Xi }_i)\underset{j=1}{\overset{P}{\Sigma }}{\alpha }_{j}k(x_j,x)---\left(46\right)$

考虑到自由故障测量变量服从正态分布且基于重构的分布在二次型中定义，则基于重构 的SPE贡献图的控制限定义如下

${\delta }_i^2\equiv g_i^{\mathrm{SPE}}{\chi }_{(h_i^{\mathrm{SPE}},\alpha )}^2---\left(47\right)$

其中

$g_i^{\mathrm{SPE}}=\frac{\mathrm{tr}{\left\{S\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\right\}}^2}{\mathrm{tr}\left\{S\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\right\}}---\left(48\right)$

$h_i^{\mathrm{SPE}}=\frac{\mathrm{tr}{\left\{S\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\right\}}^2}{\mathrm{tr}{\left\{S\tilde{\Phi }\left({\Xi }_i\right){\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}{\tilde{\Phi }}^T\left({\Xi }_i\right)\right\}}^2}---\left(49\right)$

式中S是x的协方差矩阵。

利用矩阵轨迹特性，可以把方程(48)和(49)重新排列为

$g_i^{\mathrm{SPE}}=\frac{\mathrm{tr}{\left\{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)\right\}}^2}{\mathrm{tr}\left\{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)\right\}}---\left(50\right)$

$h_i^{\mathrm{SPE}}=\frac{{\left[\mathrm{tr}\right\{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)\left\}\right]}^2}{\mathrm{tr}{\left\{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)\right\}}^2}---\left(51\right)$

方程(50)和(51)括号里的项都是标量且能够根据迹的运算求出。由以上得出

$g_i^{\mathrm{SPE}}=\frac{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-2}{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)}^2}{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-1}\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)}=\frac{{\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)}{{\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)}---\left(52\right)$

$h_i^{\mathrm{SPE}}=\frac{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-2}{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)}^2}{{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)\right)}^{-2}{\left({\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)\right)}^2}=1---\left(53\right)$

因此，δ_i²表示为：

${\delta }_i^2=\frac{{\tilde{\Phi }}^T\left({\Xi }_i\right)S\tilde{\Phi }\left({\Xi }_i\right)}{{\tilde{\Phi }}^T\left({\Xi }_i\right)\tilde{\Phi }\left({\Xi }_i\right)}{\chi }_{(1,\alpha )}^2---\left(54\right)$

通过各变量的贡献图SPE_t和对应的控制限δ_i²，就可判断出各变量对故障所承担的责任。当 $\mathrm{SPE}\_T\le {\delta }_i^2$时，则第i个变量对故障的发生要承担一定的责任。根据差值的大小判断责任的大 小。

具体操作步骤如下：

一、离线训练

(1)采集数据，并进行规范化处理；

(2)对规范化后的数据进行特征提取，得到基Ω；

(3)根据方程(15)对基Ω中的数据计算核矩阵K；

(4)根据方程(20)在特征空间执行中心化，得到中心化核矩阵

(5)根据方程(21)求解特征向量，并根据方程(23)对前p个最大正特征值所对应的特征 向量进行归一化；

二、在线监测

(6)获得新数据，并对它进行规范化得到x_t∈R^m，由[k_t]_i＝[k_t(x_t，x_i)]计算核向量k_t∈R^1×N，

其中x_i是被规范化的标准操作数据。x_i∈R^m，i＝1，2，…，N；

(7)计算中心化核向量如下：

${\tilde{k}}_t=k_t-1_{t}K-k_t{1}_N+1_t{K1}_N---\left(55\right)$

其中K由离线建模的第(3)步确定，1_t＝(1/I)[1，…，1]∈R^1×N；

(8)计算非线性主元

$t_{k,\mathrm{new}}=(v_k,\tilde{\Phi }\left(x_t\right))=\underset{i=1}{\overset{N}{\Sigma }}{\alpha }_i^k{\tilde{k}}_t(x_t,x_i)---\left(56\right)$

(9)计算新数据的检测统计量

$T^2=t_{\mathrm{new}}^T{\Lambda }^{-1}{t}_{\mathrm{new}}=[t_{1,\mathrm{new}},t_{2,\mathrm{new}},···,t_{p,\mathrm{new}}]{\Lambda }^{-1}{[t_{1,\mathrm{new}},t_{2,\mathrm{new}},···,t_{p,\mathrm{new}}]}^T---\left(57\right)$

$\mathrm{SPE}=\tilde{K}(x_{\mathrm{new}},x_{\mathrm{new}})-t_{\mathrm{new}}^T{t}_{\mathrm{new}}---\left(58\right)$

(10)确定统计量的控制限；

(11)计算SPE贡献图；

(12)确定SPE贡献图的控制限。

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

本发明——基于改进的核主元分析的非线性过程故障辩识方法的优越性表现为：

1)提出一种新的特征空间基的提取方法；

2)基于重构的故障辩识方法克服了输入空间和特征空间不能自由转换所带来的故障辩识的 困难。

附图说明

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

图中：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分析。

图2田纳西-伊斯曼过程中故障10基于原始KPCA的监测结果；

图3田纳西-伊斯曼过程中故障10基于改进的KPCA的监测结果；

图4田纳西-伊斯曼过程中故障10基于735次采样的改进的KPCA的贡献图；

图5田纳西-伊斯曼过程中故障10基于785次采样的改进的KPCA的贡献图；

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

具体实例

田纳西-伊斯曼过程

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

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

    No.    过程测量     No.     成分测量     No.   控制变量     1     2     3     4     5     6     7     8     9     10     11     12    A进料(流1)    D进料(流2)    E进料(流3)    总进料(流4)    再循环流量(流8)    反应器进料速度(流6)    反应器压力    反应器等级    反应器温度    排放速度(流9)    产品分离器温度    产品分离器液位     23     24     25     26     27     28     29     30     31     32     33     34     成分A(流6)     成分B(流6)     成分C(流6)     成分D(流6)     成分E(流6)     成分F(流6)     成分A(流9)     成分B(流9)     成分C(流9)     成分D(流9)     成分E(流9)     成分F(流9)     42     43     44     45     46     47     48     49     50     51     52     53   D进料量(流2)   E进料量(流3)   A进料量(流1)   总进料量(流4)   压缩机再循环阀   排放阀(流9)   分离器灌液流量(流10)   汽提器液体产品流量(流11)   汽提器水流阀   反应器冷却水流量   冷凝器冷却水流量   搅拌速度

 13  14  15  16  17  18  19  20  21  22   产品分离器压力   产品分离器塔底低流量(流10)   汽提器等级   汽提器压力   汽提器塔底流量(流11)   汽提器温度   汽提器流量   压缩机功率   反应器冷却水出口温度   分离器冷却水出口温度  35  36  37  38  39  40  41 成分G(流9) 成分H(流9) 成分D(流11) 成分E(流11) 成分F(流11) 成分G(流11) 成分H(流11)

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

    No.     描述   类型  No. 描述   类型     1     2     3     4     5     6     7     8     9     10     11     A/C进料比率，B成分不变(流4)     B成分，A/C进料比率不变(流4)     D的进料温度(流2)     反应器冷却水的入口温度     冷凝器冷却水的入口温度     A进料损失(流1)     C存在压力损失-可用性降低(流4)     A，B，C进料成分(流4)     D的进料温度(流2)     C的进料温度(流4)     反应器冷却水的入口温度   阶跃   阶跃   阶跃   阶跃   阶跃   阶跃   阶跃   随机变量   随机变量   随机变量   随机变量  12  13  14  15  16  17  18  19  20  21 冷凝器冷却水的入口温度 反映动态 反应器冷却水阀门 冷凝器冷却水阀门 未知 未知 未知 未知 未知 流4的阀门固定在稳态位置   随机变量   慢偏移   粘住   粘住             恒定位置

步骤一、采集数据

对于训练数据和实时工况数据采用了三分钟的采样间隔采集数据。故障10的训练数据由 480个观测数构成，实时工况数据由960个观测数构成。开始时都是没有故障的，训练数据 中，故障都是在第20次采样时引入，实时工况数据中，故障都是在第160次采样引入的。训 练数据和实时工况数据中的数据都包含了52个观测变量。随机选取了故障10的训练数据和 实时工况数据中控制变量的5组数据分别如表3和表4所示。采集数据后再对其进行规范化 处理。

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

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

步骤二、提取特征空间的基

对规范化后的训练数据进行特征空间基的提取，本例核函数采用3次多项式核函数 k(x，y)＝(x·y+1)³，γ₀＝0.05。学习的开始基中只有一个数据点Φ₁＝{x₁}，引入数据x₂，根据 方程(6)和(7)计算S₂，当$S_2>\sqrt{0.05}$时，数据x₂被加进基中，得到新基Φ₂＝{x₁，x₂}，否 则，数据x₂被拒绝，Φ₂＝Φ₁。接着引入x₃，根据方程(6)和(7)计算S₃，当$S_3>\sqrt{0.05}$时， 数据x₃被加入基中，得到新基Φ₃＝{Φ₂，x₃}，否则，数据x₃被拒绝，Φ₃＝Φ₂。按此方法依次 计算，直到计算完480组数据。最终得到基Ω。

步骤三、提取非线性主元

选取径向基核函数$k(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{5^2}),$对基Ω中的数据根据方程(15)进行核计算， 得到Gram核矩阵K。接着根据公式(20)对Gram核矩阵K进行中心化处理，得到中心化 的核矩阵通过计算出了它的特征值，选取了最大的12个特征值以达到降维的目的。 根据方程(21)计算出了选取的12个特征值所对应的特征向量，并根据方程(23)对特征向 量进行归一化得到归一化后的特征向量α_k，k＝1，2，…，12。

参数都确定之后，就可以对规范化后的新数据进行非线性主元提取，根据公式 [k_t]_i＝[k_t(x_t，x_i]计算新数据的核向量，其中x_t是规范化后的新数据，x_i是训练数据， i＝1，2，…，316。然后根据方程(55)对其进行了中心化处理，得到中心化核向量最后根 据方程(56)计算得到新数据的非线性主元。

步骤四、故障检测

首先通过SPE统计量对田纳西-伊斯曼过程进行检测，根据方程(58)计算SPE统计， 接着根据方程(28)确定它的控制限，得到改进的KPCA对田纳西-伊斯曼的监测结果如图3 所示，从图中看出SPE统计量在大约170s检测到故障。本发明还利用原始KPCA对故障进 行检测，原始KPCA也是在大约170s检测到故障，延时时间是10s。如图2所示。由此可以 得出，改进的KPCA能达到原始KPCA同样的效果，但改进的KPCA检测结果基于了更少的 样本，减少了计算量。

步骤五、故障辨识

由步骤四的结果得知检测到故障发生，因此采用SPE贡献图确定哪些变量引起了故障的 发生。针对每个变量分别计算它们的贡献图和控制限，本发明随机选取了第735次采样的数 据对各变量进行计算，首先，对第一个变量A进料(流1)进行计算，看看它对故障所承担的责 任，假定故障方向Ξ₁＝(1，0，0，…，0)^T∈R⁵²，根据方程(43)和(54)分别计算得到A进料(流 1)的SPE_t和控制限δ₁²；接着对第二个变量D进料(流2)进行计算，看看它对故障所承担的责 任，假定故障方向Ξ₂＝(0，1，0，…，0)^T∈R⁵²，根据方程(43)和(54)分别计算得到D进料(流 2)的SPE_t和控制限δ₂²；按此方法依次对田纳西-伊斯曼的52个变量进行计算，得到贡献图 如图4所示。当各变量的SPE_t＞δ²时，那么此变量要对故障的发生承担一定的责任，根据 SPE_t-δ²的大小判断责任的轻重。从图中可以看出第1-6，8-13，15-17，19，21，22，24，28，30， 32，46-48，52个变量失控。接着又随机选取了第785次采样的数据对个变量进行计算，用上面 同样的方法，得到结果如图5所示。从图中知道第1-17，19，21，22，24-26，28，30-33，38，43， 45-48，52个变量失控。因此，通过此方法就可以很明显的从图中看出各变量对故障的发生所 要承担的责任，就不用盲目的对整个过程的所有变量进行检查。

通过上面的实例，表明了本发明——基于改进的核主元分析的非线性过程故障辨识方法的 有效性以及同原始KPCA相比的优越性。它不仅实现了故障的检测，还辨识出各变量失控对 故障所承担的责任，减少了设备维修的盲目性。