CN101308385B - 基于二维动态核主元分析的非线性过程故障检测方法 - Google Patents

Abstract

一种基于二维动态核主元分析的非线性过程故障检测方法，属于故障检测技术领域，包括以下步骤，步骤一：确定采样参数，判断执行过程。确定采样参数，选择影响故障的数据参数，判断是进行训练还是检测；步骤二：训练，即采集正常工作的数据，用二维动态核主元分析方法提取训练数据的非线性主元，计算训练数据的平方预测误差，并确定控制限；步骤三：检测，即采集在线观测数据，利用二维动态核主元分析方法提取在线观测数据的非线性主元，并计算实时在线观测数据的平方预测误差，比较实时在线观测数据的平方预测误差与训练数据平方预测误差的控制限，如果超出控制限，则显示并报警。本发明能够及时检测到生产过程中故障，减少工业生产过程中的损失。

Description

基于二维动态核主元分析的非线性过程故障检测方法

技术领域

本发明属于故障检测技术领域，特别涉及一种基于二维动态核主元分析(2D-DKPCA，Two Dimensional Dynamic Kernel Principal Components Analysis)的非线性过程故障检测方法。

背景技术

随着经济的不断发展，工业企业不断增多，生产部门之间的竞争也日益激烈。要提高企业的效益就要提高产品质量，减少不合格产品的数量。除了对产品生产工艺的要求越来越高，对生产过程的故障检测方法能力的要求也越来越高。这是因为在复杂非线性生产过程中，一些故障是难以避免的，比如外界干扰导致生产条件偏离初始设定的最优条件。所以及时地检测到生产过程中发生的故障，尽早实施补救措施是十分重要的。

许多非线性过程数据具有动态模式，且数据之间具有相关性，所以静态的或线性的检测方法就不能精确检测到故障。

对于静态非线性过程的故障检测方法也有相关的研究。核主元分析(KPCA，KernelPrincipal Components Analysis)方法也用来提取非线性关系。KPCA可以在高维特征空间中通过积分算子和非线性核函数有效的计算主元(PCs，Principal Components)。以上方法都适用于静态非线性过程的故障检测中。

所以，为了解决动态非线性过程的故障检测问题，应当发展一个基于核主元分析法的故障检测方法时应该考虑过程的动态性。

发明内容

本发明提出了一种基于二维动态核主元分析的非线性过程故障检测方法。用二维动态核主元分析方法从正常运行的过程数据中提取一些捕捉非线性和动态性的主元，并通过平方预测误差(SPE，Squared Prediction Error)统计变量的分布检测是否有故障发生。

使用二维动态核主元分析方法进行在线监测批处理过程的方法的详细步骤如下：

步骤一：确定采样参数

确定采样参数，选择影响故障的数据作为采样参数，采样参数的选择依据每个具体生产过程可能导致故障的数据；然后人工选择，判断是进行训练还是检测，训练则转入步骤二，检测则转入步骤三。

训练数据是指正常工作过程的历史数据，检测数据是指实际工作过程中的在线观测数据。

步骤二：训练，即采集正常工作的数据，用二维动态核主元分析方法提取训练数据的非线性主元，计算训练数据的平方预测误差，并确定其控制限。

1.用变量的均值和标准偏差来规范化正常情况下的训练数据，经过处理后，数据值都在0～1之间。

2.构造二维增广矩阵，规范化以后的训练数据与滞后区域的数据构成二维增广矩阵。

使用KPCA方法通过对角化协方差矩阵来对一组给定的非线性相关的数据集进行退耦，应该在线性特征空间中表示协方差矩阵而不是在非线性输入空间中表示，所以首先把输入向量从输入空间投影到线性特征空间中。

设Φ(·)是把输入向量从输入空间投影到特征空间F的非线性映射函数。观测变量数据由I个连续的批，J个变量和K次采样产生。设x_j(i，k)为第i个批中在第k次采样时变量j的测量值，i＝1…I；j＝1…J；k＝1…K。在二维动态过程的批处理中，投影到特征空间后变量的当前值不仅仅依赖于时间方向的过去值x_j(i，k-1)，x_j(i，k-2)，...，x_j(i，k-n+1)；也取决于批次方向的过去值x_j(i-1，k)，x_j(i-2，k)，...，x_j(i-m+1，k)，甚至交叉方向的过去值x_j(i-1，k-1)，...，x_j(i-m+1，k-n+1)，其中m和n是两个方向的自回归阶数。观测变量投影到特征空间后，由其当前值和时间、批次、交叉方向的滞后值构成的二维增广矩阵定义如下：

$\mathrm{\Φ}\left(X\right)=\left[\begin{array}{c}\mathrm{\Φ}\left(x_{m,n}\right)\\ \·\\ \·\\ \·\\ \mathrm{\Φ}\left(x_{i,k}\right)\\ \·\\ \·\\ \·\\ \mathrm{\Φ}\left(x_{l,k}\right)\end{array}\right]=\left[\begin{array}{c}\mathrm{\Φ}\left(y_1\right)\\ \·\\ \·\\ \·\\ \mathrm{\Φ}\left(y_l\right)\\ \·\\ \·\\ \·\\ \mathrm{\Φ}\left(y_L\right)\end{array}\right]---\left(1\right)$

其中，滞后区域的数据指在本次采样之前的数据，比如时间方向的自回归系数是2，那么就往前看两步，为了便于之后的运算，用一个脚标的x_l替换两个脚标的x_m，n，它们之间的换算关系式如下

l＝(i-m)(K-n+1)+k-n+1，l＝1…L    (2)

L＝(I-m+1)(K-n+1)                 (3)

x_l＝y_l＝[x(i，k)，x(i，k-1)，…，

x(i，k-n+1)，x(i-1，k)，x(i-1，k-1)，…，

x(i-1，k-n+1)，…，x(i-m+1，k)， (4)

x(i-m+1，k-1)，…，x(i-m+1，k-n+1)]

x(i，k)＝[x₁(i，k)，x₂(i，k)…，x_J(i，k)](5)

其中，x_j(i，k)为第i个批中在第k次采样时变量j的测量值，x(i，k)是由第i个批第k次采样时得到的J个观测变量构成的向量。x_i，k是由当前的x(i，k)和x(i，k)分别在时间方向上、批次方向上和交叉方向上的滞后值构成的向量，为了便于计算，用y_l替换x_i，k，即y_l＝x_i，k。Φ(·)是把输入向量从输入空间投影到特征空间F的非线性映射函数。Φ(y_l)表示向量y_l在特征空间上的投影。L为常数，是增广矩阵Φ(X)的行数。

3.计算中心化核矩阵，求取协方差矩阵的特征值和特征向量，确定保留主元的个数。

在线性特征空间中，映射后的观测变量的协方差矩阵定义如式(6)所示

$C^F=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\Σ}}}\mathrm{\Φ}\left(y_d\right)\mathrm{\Φ}{\left(y_d\right)}^T---\left(6\right)$

其中，Φ(y_d)代表由映射到特征空间后的观测变量增广矩阵的第d个行向量，d＝1，…，L。

为了对角化协方差矩阵，下面求解协方差矩阵的特征值。

λv＝C^Fv    (7)

其中，特征值λ≥0，且v∈F\{0}是λ对应的非零特征向量。C^Fv表示如下：

$C^{F}v=\left(\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\Σ}}}\mathrm{\Φ}\left(y_d\right)\mathrm{\Φ}{\left(y_d\right)}^T\right)v$

(8)

$=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}<\mathrm{\&Phi;}\left(y_d\right),v>\mathrm{\&Phi;}\left(y_d\right)$

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

λ<Φ(y_d)，v>＝<Φ(y_d)，C^Fv>       (9)

在KPCA中，协方差矩阵的特征向量是Φ(y_d)的线性组合，则一定存在系数 ${\mathrm{\&alpha;}}_d=({\mathrm{\&alpha;}}_d^1,{\mathrm{\&alpha;}}_d^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&alpha;}}_d^L),$ d＝1，…，L，使得

$v=\left[\begin{array}{c}{v}^1\\ v^2\\ .\\ .\\ .\\ v^L\end{array}\right]=\left[\begin{array}{c}\underset{d=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^1\mathrm{\&Phi;}\left(y_d\right)\\ \underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^2\mathrm{\&Phi;}\left(y_d\right)\\ .\\ .\\ .\\ \underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^L\mathrm{\&Phi;}\left(y_d\right)\end{array}\right]---\left(10\right)$

结合方程(9)和(10)，得到

$\mathrm{\&lambda;}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d<\mathrm{\&Phi;}\left(y_g\right),\mathrm{\&Phi;}\left(y_d\right)>$

(11)

$=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d<\mathrm{\&Phi;}\left(y_g\right),\underset{h=1}{\overset{L}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(y_h\right)><\mathrm{\&Phi;}\left(y_h\right),\mathrm{\&Phi;}\left(y_d\right)>$

其中，g＝1，…，L。

尽管映射向量Φ(·)存在，但通常计算上并不容易处理。然而，通常不需要明确计算Φ(·)，仅仅需要知道特征空间中两个向量的点积。为了避免执行非线性映射Φ(·)及在特征空间中的点积运算，下面通过引入核函数k(x，y)＝<Φ(x)，Φ(y)>来简化计算。

定义一个L×L维的矩阵K，对K执行高维空间中的中心化。很容易得到中心化核矩阵

$\tilde{K}=K-1_{L}K-K{1}_L+1_{L}K{1}_L---\left(12\right)$

其中， $1_L=\frac{1}{L}\left[\begin{array}{ccccc}1& .& .& .& 1\\ .& .& & & .\\ .& & .& & .\\ .& & & .& .\\ 1& .& .& .& 1\end{array}\right]\&Element;R^{\mathrm{LL}}.$

则方程(11)可以转换为：

$\mathrm{\&lambda;}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d{\tilde{k}}_{\mathrm{gd}}=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d\underset{h=1}{\overset{L}{\mathrm{\&Sigma;}}}{\tilde{k}}_{\mathrm{gh}}{\tilde{k}}_{\mathrm{hd}}---\left(13\right)$

其中， ${\tilde{k}}_{\mathrm{gd}}=<\mathrm{\&Phi;}\left(y_g\right),\mathrm{\&Phi;}\left(y_d\right)>,$ ${\tilde{k}}_{\mathrm{gh}}=<\mathrm{\&Phi;}\left(y_g\right),\mathrm{\&Phi;}\left(y_h\right)>,$ ${\tilde{k}}_{\mathrm{hd}}=<\mathrm{\&Phi;}\left(y_h\right),\mathrm{\&Phi;}\left(y_d\right)>.$

方程左边为

$\mathrm{\&lambda;}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d{\tilde{k}}_{\mathrm{gd}}=\mathrm{\&lambda;}\tilde{K}\mathrm{\&alpha;}---\left(14\right)$

方程右边为

$\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d\underset{h=1}{\overset{L}{\mathrm{\&Sigma;}}}{\tilde{k}}_{\mathrm{gh}}{\tilde{k}}_{\mathrm{hd}}=(1/L){\tilde{K}}^2\mathrm{\&alpha;}---\left(15\right)$

左边等于右边，得到

$\mathrm{\&lambda;L}\tilde{K}\mathrm{\&alpha;}={\tilde{K}}^2\mathrm{\&alpha;}---\left(16\right)$

为了解方程(16)，我们求解下式的特征值的问题

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

现在，在特征空间执行KPCA等价于解决方程(17)的特征值问题。通过仅保留前A个特征向量来降低问题的维数，A≤L。并且按特征空间中相应向量归一化的要求对前A个特征向量进行归一化。

4.计算训练数据的非线性动态主元

设x为规范化后的变量，

x＝[x(i，k)，x(i，k-1)，…，

x(i，k-n+1)，x(i-1，k)，x(i-1，k-1)，…，(18)

x(i-1，k-n+1)，…，x(i-m+1，k)，

x(i-m+1，k-1)，…，x(i-m+1，k-n+1)]

x(i，k)＝[x₁(i，k)，x₂(i，k)…，x_J(i，k)](19)

t_k表示x的第k个非线性主元，利用下式计算t_k，

$t_k=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_k}}<v^k,\mathrm{\&Phi;}\left(x\right)>$

$=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_k}}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^k<\mathrm{\&Phi;}\left(y_d\right),\mathrm{\&Phi;}\left(x\right)>---\left(20\right)$

$=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_k}}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^k\tilde{k}(y_d,x)$

其中，k＝1，…，A。Φ(x)表示映射到特征空间后的观测变量。v^k表示协方差矩阵对应的第k个特征向量，α_d ^k表示第k个特征向量的第d个变量对应的系数。

5.计算训练数据的平方预测误差统计量，并确定其控制限

通过平方预测误差SPE的分布可以检测到是否有故障发生。当观测数据的统计量超出统计量规定的上限时属于异常数据，表明发生了故障。

特征空间中的平方预测误差SPE定义为：

$\mathrm{SPE}={\left|\right|e\left|\right|}^2=\tilde{k}(x,x)-t^{T}t---\left(21\right)$

其中，t＝(t₁，t₂，…，t_A)

当 $\mathrm{SPE}\&le;{\mathrm{\&delta;}}_{\mathrm{\&alpha;}}^2$ 时，表示没有故障发生。其中，δ_α ²代表置信度为α的SPE的控制上线，由权重χ²分布来计算。

这里采用Nomikos和MacGregor给出的δ_α ²的计算式：

${\mathrm{\&delta;}}_{\mathrm{\&alpha;}}^2=g{\mathrm{\&chi;}}_{h,\mathrm{\&alpha;}}^2---\left(22\right)$

其中，g＝v_s/2m_s，h＝2(m_s)²/v_s，m_s是K次采样的SPE值的平均值，v_s是相应的方差。

根据正常情况建模具体算法如下：

(1)获得观测数据子集；

(2)用每个变量的均值和标准偏差来规范化数据子集；

(3)根据规范化的数据计算核矩阵，并由式(12)计算中心化核矩阵

(4)根据式(17)计算协方差矩阵的特征向量，并选取要保留的特征向量的个数。

(5)根据式(20)提取动态非线性主元；

(6)由式(21)计算正常运行的监测统计量(SPE)；

(7)根据等式(22)确定SPE图表的控制界限。

步骤三：检测，即采集在线观测数据，利用二维动态核主元分析方法提取训练数据的非线性主元，计算在线观测数据的SPE，与训练数据SPE的控制限相比较，如果超出控制限，则显示并报警。

在线监测的具体算法如下：

(1)在线采集观测变量数据；

(2)用每个变量的均值和标准偏差来规范化数据子集；

(3)计算规范化后的在线观测数据的核向量k_new∈R^1×L

[k_new]_g＝[k(x_new，y_g)]         (23)

其中，x_new∈R^mnJ是规范化后的在线观测数据向量，y_g∈R^mnJ是正常情况的建模数据向量，[k_new]_g表示核向量的第g个元素，g＝1，…，L；

(4)根据式(24)计算中心化核向量

${\tilde{k}}_{\mathrm{new}}=k_{\mathrm{new}}-1_{\mathrm{new}}K-k_{\mathrm{new}}{1}_L+1_{\mathrm{new}}K{1}_L---\left(24\right)$

其中， $1_{\mathrm{new}}=\frac{1}{L}[1,\&CenterDot;\&CenterDot;\&CenterDot;,1]\&Element;R^{1\&times;L},$ K是由建模数据计算得到的核矩阵， $1_L=\frac{1}{L}\left[\begin{array}{ccccc}1& .& .& .& 1\\ .& .& & & .\\ .& & .& & .\\ .& & & .& .\\ 1& .& .& .& 1\end{array}\right]\&Element;R^{\mathrm{LL}}.$

(5)根据式(20)提取在线观测数据的非线性主元。

(6)根据式(21)计算在线观测数据的SPE监测统计变量。

(7)将在线观测数据的SPE与训练数据平方预测误差的SPE相比较，如果超出控制限，则显示并报警。

本发明具有以下优势：

1.本发明首次提出了一种基于2D-DKPCA的非线性过程故障检测方法，实现对复杂过程的监测。

2.本发明能够全面的捕捉批处理过程中的非线性特征和动态性特征，从而更有效的监测生产过程。

3.本发明能够及时的检测到生产过程中故障，从而减少工业生产过程中的损失。

本发明的研究成果还可以应用到其他生产过程中。

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

附图说明

图1.本发明的一个实施例从第5批开始具有过程漂移的2D-DPCA方法SPE检测图；

图2.本发明的一个实施例从第5批开始具有过程漂移的2D-DKPCA方法SPE检测图；

图3.本发明的一个实施例从第5批开始相关结构发生变化的2D-DPCA方法SPE检测图；

图4.本发明的一个实施例从第5批开始相关结构发生变化的2D-DKPCA方法SPE检测图；

图5.青霉素的发酵过程示意图；其中，FC：流量控制；T：温度指示器；pH：pH值指示器；

图6.本发明的另外一个实施例从第2批的第100次采样到400次采样之间，酶的供给率减少为正常的80％的二维动态核主元分析(2D-DKPCA)法的故障检测图；

图7.本发明的另外一个实施例从第三批的100小时开始，酶的供给率以-0.005的斜率线性递减的2D-DKPCA方法SPE检测图；

图8.本发明的另外一个实施例在第五批的100小时处，酶的供给率减少为正常的80％的2D-DKPCA方法SPE检测图；

图9.本发明的2D-DKPCA故障检测方法流程图。

具体实施方式

下面参照附图，说明本发明在实例中的应用。

实例1.

为了模拟具有二维动态非线性特征的批处理过程，本例创建了一个非线性系统来检验2D-DKPCA方法的性能。系统模型如下：

x₁(i，k)＝0.5sin(x₁(i，k-1)+0.8x₁(i-1，k)-0.33x₁(i-1，k-1)+w₁

$x_2(i,k)=0.67\sin \left(x_2(i,k-1)\right)+0.44x_2(i-1,k)-0.11\cos \left(x_2(i-1,k-1)\right)+\sin \left(\frac{\mathrm{k\&pi;}}{0.3K}\right)+w_2$

x₃(i，k)＝0.65x₁(i，k)+0.35x₂(i，k)+w₃

x₄(i，k)＝-1.26x₁(i，k)+0.33x₂(i，k)+w₄

其中，x₁和x₂是两个独立的变量，x₃和x₄是x₁与x₂的线性组合；w_j(j＝1，2，3，4)是方差为0.02的高斯型随机变量。总共有10个批次，每个批次有50次采样，即I＝10，J＝4，K＝50。对于每批的第一次采样，设变量x₁和x₂分别为x₁(i，1)＝x₁(i-1，1)+w₁和x₂(i，1)＝x₂(i-1，1)+w₂。变量的2-D自回归阶数都是(1，1)，二维增广数据向量由映射后的x(i，k)，x(i，k-1)，x(i-1，k)，x(i-1，k-1)构成，其中x(i，k)＝[x₁(i，k)，x₂(i，k)，x₃(i，k)，x₄(i，k)]。在这里，我们讨论两种故障，即过程的漂移和相关结构的变化。

(1)小的过程漂移

为了产生一些带有故障的批次，从第5批开始给变量x₂加上一个随着时间和批次缓慢增加的信号。我们保留4个主元。在这个例子中，设定参数c＝6000。图1和2分别为使用2D-DPCA方法和2D-DKPCA方法的监测结果。在图2中，SPE值随着时间和批次的增加而增加，这表明有小的过程漂移。

(2)结构变化

另外一个故障是通过改变相关结构产生的。变量x₂从第5批开始变化，即

$x_2(i,k)=0.8\cos \left(x_2(i,k-1)\right)+0.67x_2(i-1,k)-0.47\cos \left(x_2(i-1,k-1)\right)+\sin \left(\frac{\mathrm{k\&pi;}}{0.35K}\right)$

应用2D-DPCA方法进行故障监测，保留4个主元，监测结果如图3所示。再使用2D-DKPCA方法，设定参数c＝500，我们得到了明显反应了相关结构变化的SPE图表，如图4所示。

通过这两种故障的对比，可以很明显的看出2D-DKPCA方法在非线性动态过程故障检测上的优势。

实例2.

在这个例子中，2D-DKPCA方法应用于一个著名的标准过程的监测——青霉素发酵过程中。图5给出了青霉素的发酵过程流程图。青霉素的产生过程具有非线性特征。在青霉素发酵的初始阶段，在发酵罐中放入必要的细胞物质。当原有的酶被微生物消耗掉的时候，开始添加酶。由于代谢阻遏物影响，发酵罐中低浓度的培养基是高发酵率的重要保证。在发酵开始连续的供给葡萄糖。引入小的变动，监测青霉素发酵过程在故障条件下的运行情况。监测的变量如表1所示。每个批的持续时间是400小时，包含约45小时的开始阶段和约355小时的反馈阶段。在第2批引入故障，即从100小时到400小时，酶的供给率减少为正常的80％。

表1.青霉素发酵的监测变量

No.

变量

1	通风率(L/h)
		2	搅拌功率(W)
3	底物流加速率(L/h)
		4	底物流温度(k)
5	溶解氧的浓度(mmol/L)
		6	培养容积(L)
7	CO2浓度(mmol/L)
		8	pH值
9	产生的热量(kcal)

步骤一：数据采集

批2的观测变量数据集由400组数据构成，包含了9个监测变量。本实例主要针对霉的供给速率故障进行实例分析，对故障批2随机选取了10组观测变量数据如表2所示。

表2.观测变量数据

通风率	搅拌功率	底物流加速率	底物流温度	溶解氧的浓度	培养容积	CO2浓度	pH值	产生的热量
									8.602	29.79	0.040613	295.83	95.6857	101.4995	2.0156	5.0349	67.1035
8.6035	29.78	0.040613	295.835	96.3238	101.5152	2.2144	5.0304	67.162
									8.606	29.775	0.040565	295.84	95.8797	101.5309	2.0999	5.0261	67.2191
8.608	29.775	0.040518	295.845	96.3414	101.5457	2.1732	5.0219	67.2746
									8.61	29.765	0.040494	295.8475	95.9682	101.5606	2.248	5.0179	67.329
8.6115	29.76	0.040518	295.85	96.0833	101.5761	2.0803	5.014	67.3821
									8.6135	29.76	0.040541	295.85	95.8888	101.5915	2.1131	5.0102	67.4339

8.6155	29.755	0.040589	295.8525	96.2489	101.6066	2.0992	5.0065	67.4853
									8.617	29.75	0.040613	295.8525	95.9859	101.6223	2.1171	5.0028	67.5369
8.6195	29.755	0.040636	295.8525	95.961	101.6381	2.1951	5.2386	67.5897

步骤二：数据处理

根据计算观测变量均值和标准偏差，并用变量的均值和标准偏差来规范化观测变量数据；

步骤三：提取动态非线性主元素

首先把输入向量从输入空间投影到线性特征空间中。设Φ(·)是把输入向量从输入空间投影到特征空间的非线性映射函数。在二维动态过程的批处理中，投影到特征空间后观测变量的当前值依赖于时间方向的过去值、批次方向的过去值和交叉方向的过去值。所以在线观测变量投影到特征空间后，就构成了形如式(18)所示的向量。

在离线部分计算协方差矩阵C^F时选取径向基核函数 $k(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{c})$ 作为核函数，其中c＝5000，并保留前5个特征向量。对于在线观测数据，使用2D-DKPCA方法对第2批运行情况进行在线监控。根据公式(23)计算核向量k_new。根据公式(24)计算中心化核向量

根据式(20)提取在线观测数据的非线性主元。

步骤四：利用统计变量SPE进行故障检测

通过统计变量SPE的分布可以检测到是否有故障发生。当观测数据的统计量超出统计量规定的上限时属于异常数据，表明发生了故障。根据式(21)计算在线观测数据的SPE监测统计变量。

2D-DKPCA方法的检测结果如图6所示。由图可见，从发酵开始到45小时，由于每个批的启动运行条件有偏差，所以SPE上下波动。对应的SPE在第100小时处明显超出控制线，表明有故障发生，说明此方法能够检测到生产过程中的故障。

再次测试2D-DKPCA方法，在第3和5批中分别引入故障。在第三批中，酶的供给率从100小时开始，以-0.005的斜率线性递减，持续到结束；在第五批中，酶的供给率从100小时开始，减少为正常的80％，持续到结束。核函数选用的都是径向基核函数 $k(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{c}),$ 其中c＝500。

检测结果分别如图7和如图8所示。如图7所示，2D-DKPCA方法在180小时检测到故障，比故障产生时间延迟了80小时。这是由于酶的供给率以-0.005的斜率线性递减，由于这种变化太微小了，所以在很长时间内酶的供给率仍在正常范围内，没有形成故障。如图8所示，2D-DKPCA方法在100小时检测到故障，即在故障发生时立即检测到故障。

本发明提出了一个基于二维动态核主元分析(2D-DKPCA)的故障检测方法。这个方法能够全面的检测系统的故障。本方法适用于非线性动态生产过程，比如发酵过程，制药过程等。

Claims (2)

1.一种基于二维动态核主元分析的非线性过程故障检测方法，该方法包括以下步骤：步骤一：确定采样参数，判断执行过程

选择影响故障的数据作为采样参数，判断是进行训练还是检测，训练转入步骤二，检测转入步骤三；

步骤二：训练，即采集正常工作的数据，用二维动态核主元分析方法提取训练数据的非线性主元，计算训练数据的平方预测误差，并确定其控制限；

步骤三：检测，即采集在线观测数据，利用二维动态核主元分析方法提取在线观测数据的非线性主元，并计算实时在线观测数据的平方预测误差，比较实时在线观测数据的平方预测误差与训练数据平方预测误差的控制限，如果超出控制限，则显示并报警；其特征在于步骤二中所述的训练过程如下：

1)用变量的均值和标准偏差来规范化采集数据，经过处理后，数据值都在0～1之间；

2)构造二维增广矩阵，规范化以后的数据与滞后区域的数据构成二维增广矩阵，二维增广矩阵定义如下：

$\mathrm{\&Phi;}\left(X\right)=\left[\begin{array}{c}\mathrm{\&Phi;}\left(x_{m,n}\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&Phi;}\left(x_{i,k}\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&Phi;}\left(x_{I,K}\right)\end{array}\right]=\left[\begin{array}{c}\mathrm{\&Phi;}\left(y_1\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&Phi;}\left(y_l\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&Phi;}\left(y_L\right)\end{array}\right]---\left(1\right)$

l＝(i-m)(K-n+1)+k-n+1，l＝1…L            (2)

L＝(I-m+1)(K-n+1)                         (3)

x_l＝y_l＝[x(i，k)，x(i，k-1)，…，

x(i，k-n+1)，x(i-1，k)，x(i-1，k-1)，…，

                                          (4)

x(i-1，k-n+1)，…，x(i-m+1，k)，

x(i-m+1，k-1)，…，x(i-m+1，k-n+1)]

x(i，k)＝[x₁(i，k)，x₂(i，k)…，x_J(i，k)] (5)

设Φ(·)是把输入向量从输入空间投影到特征空间F的非线性映射函数，观测变量数据由I个连续的批，J个变量和K次采样产生，设x_j(i，k)为第i个批中在第k次采样时变量j的测量值，i＝1…I；j＝1…J；k＝1…K，其中m和n是两个方向的自回归阶数，

x_i，k是由当前的x(i，k)和x(i，k)分别在时间方向上、批次方向上和交叉方向上的滞后值构成的向量，y_l＝x_i，k，Φ(y_l)表示向量y_l在特征空间F上的投影，L为常数；

3)计算中心化核矩阵，求取协方差矩阵的特征值和特征向量，确定保留主元的个数在特征空间F中，映射后的观测变量的协方差矩阵定义如式(6)所示

$C^F=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(y_d\right)\mathrm{\&Phi;}{\left(y_d\right)}^T---\left(6\right)$

其中，Φ(y_d)代表由映射到特征空间F后的观测变量增广矩阵的第d个行向量，d＝1，…，L；下面求解协方差矩阵的特征值，

λv＝C^Fv            (7)

其中，特征值λ≥0，且v∈F\{0}是λ对应的非零特征向量，C^Fv表示如下：

$C^{F}v=\left(\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(y_d\right)\mathrm{\&Phi;}{\left(y_d\right)}^T\right)v$

                   (8)

$=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}<\mathrm{\&Phi;}\left(y_d\right),v>\mathrm{\&Phi;}\left(y_d\right)$

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

λ<Φ(y_d)，v>＝<Φ(y_d)，C^Fv>        (9)

协方差矩阵的特征向量是Φ(y_d)线性组合，存在系数

d＝1，…，L，使得

$v=\left[\begin{array}{c}{v}^1\\ v^2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ v^L\end{array}\right]=\left[\begin{array}{c}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^1\mathrm{\&Phi;}\left(y_d\right)\\ \underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^2\mathrm{\&Phi;}\left(y_d\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^L\mathrm{\&Phi;}\left(y_d\right)\end{array}\right]---\left(10\right)$

结合方程(9)和(10)，得到

$\mathrm{\&lambda;}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d<\mathrm{\&Phi;}\left(y_g\right),\mathrm{\&Phi;}\left(y_d\right)>$

                 (11)

$=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d<\mathrm{\&Phi;}\left(y_g\right),\underset{h=1}{\overset{L}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(y_h\right)><\mathrm{\&Phi;}\left(y_h\right),\mathrm{\&Phi;}\left(y_d\right)>$

其中，g＝1，…，L；引入核函数k(x，y)＝<Φ(x)，Φ(y)>，

定义一个L×L维的矩阵K，对K执行高维空间中的中心化，得到中心化核矩阵

$\tilde{K}=K-1_{L}K-K{1}_L+1_L{K1}_L---\left(12\right)$

其中，

则方程(11)转换为：

$\mathrm{\&lambda;}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d{\tilde{k}}_{\mathrm{gd}}=\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d\underset{h=1}{\overset{L}{\mathrm{\&Sigma;}}}{\tilde{k}}_{\mathrm{gh}}{\tilde{k}}_{\mathrm{hd}}---\left(13\right)$

其中， ${\tilde{k}}_{\mathrm{gd}}=<\mathrm{\&Phi;}\left(y_g\right),\mathrm{\&Phi;}\left(y_d\right)>,$ ${\tilde{k}}_{\mathrm{gh}}=<\mathrm{\&Phi;}\left(y_g\right),\mathrm{\&Phi;}\left(y_h\right)>,$ ${\tilde{k}}_{\mathrm{hd}}=<\mathrm{\&Phi;}\left(y_h\right),\mathrm{\&Phi;}\left(y_d\right)>;$

$\mathrm{\&lambda;}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d{\tilde{k}}_{\mathrm{gd}}=\mathrm{\&lambda;}\tilde{K}\mathrm{\&alpha;}---\left(14\right)$

$\frac{1}{L}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d\underset{h=1}{\overset{L}{\mathrm{\&Sigma;}}}{\tilde{k}}_{\mathrm{gh}}{\tilde{k}}_{\mathrm{hd}}=(1/L){\tilde{K}}^2\mathrm{\&alpha;}---\left(15\right)$

$\mathrm{\&lambda;L}\tilde{K}\mathrm{\&alpha;}={\tilde{K}}^2\mathrm{\&alpha;}---\left(16\right)$

为了解方程(16)，求解下式的特征值

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

保留前A个特征向量来降低问题的维数，A≤L，对前A个特征向量进行归一化；

4)计算训练数据的非线性动态主元

设x为规范化后的变量，

x＝[x(i，k)，x(i，k-1)，…，

x(i，k-n+1)，x(i-1，k)，x(i-1，k-1)，…，

                                                        (18)

x(i-1，k-n+1)，…，x(i-m+1，k)，

x(i-m+1，k-1)，…，x(i-m+1，k-n+1)]

x(i，k)＝[x₁(i，k)，x₂(i，k)…，x_J(i，k)]               (19)

t_k表示x的第k个非线性主元，利用下式计算t_k，

$t_k=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_k}}<v^k,\mathrm{\&Phi;}\left(x\right)>$

$=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_k}}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^k<\mathrm{\&Phi;}\left(y_d\right),\mathrm{\&Phi;}\left(x\right)>---\left(20\right)$

$=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_k}}\underset{d=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d^k\tilde{k}(y_d,x)$

其中，k＝1，…，A，v^k表示协方差矩阵对应的第k个特征向量，

表示第k个特征向量的第d个变量对应的系数；

5)计算训练数据的平方预测误差统计量，并确定其控制限

通过平方预测误差SPE的分布可以检测到是否有故障发生，当观测数据的统计量超出统计量规定的上限时属于异常数据，表明发生了故障；

特征空间F中的平方预测误差SPE定义为：

$\mathrm{SPE}={\left|\right|e\left|\right|}^2=\tilde{k}(x,x)-t^{T}t---\left(21\right)$

其中，t＝(t₁，t₂，…，t_A)

当

时，表示没有故障发生，其中，代表置信度为α的SPE的控制限；

的计算式如下：

${\mathrm{\&delta;}}_{\mathrm{\&alpha;}}^2=g{\mathrm{\&chi;}}_{h,\mathrm{\&alpha;}}^2---\left(22\right)$

其中，g＝v_s/2m_s，h＝2(m_s)²/vs，m_s是K次采样的SPE值的平均值，v_s是相应的方差。

2.按照权利要求1所述的一种基于二维动态核主元分析的非线性过程故障检测方法，其特征在于步骤三中所述的检测过程，

(1)在线采集观测变量数据；

(2)用每个变量的均值和标准偏差来规范化数据子集；

(3)计算规范化后的在线观测数据的核向量k_new∈R^1×L

[k_new]_g＝[k(x_new，y_g)]        (23)

其中，x_new∈R^mnJ是规范化后的在线观测数据向量，y_g∈R^mnJ是正常情况的建模数据向量，[k_new]_g表示核向量的第g个元素，g＝1，…，L；

(4)根据式(24)计算中心化核向量

${\tilde{k}}_{\mathrm{new}}=k_{\mathrm{new}}-1_{\mathrm{new}}K-k_{\mathrm{new}}{1}_L+1_{\mathrm{new}}K{1}_L---\left(24\right)$

其中，

K是由建模数据计算得到的核矩阵，

(5)提取在线观测数据的非线性主元；

(6)计算在线观测数据的平方预测误差SPE监测统计变量；平方预测误差SPE定义为：

$\mathrm{SPE}={\left|\right|e\left|\right|}^2=\tilde{k}(x,x)-t^{T}t$

其中，t＝(t₁，t₂，…，t_A)

(7)将在线观测数据的平方预测误差SPE与训练数据平方预测误差的控制限相比较，如果超出控制限，则显示并报警。

Images