CN103838217B - 一种基于mica-ocsvm的发酵过程故障监测方法 - Google Patents

Abstract

本发明公开了一种实现对青霉素发酵过程进行实时故障监测的新方法。为了保证青霉素发酵过程的安全性和稳定性，建立一种有效的过程监测方案来及时地检测异常现象是十分有必要的。本发明包括“离线建模”和“在线监测”两个阶段。“离线建模”包括：首先对发酵过程的三维数据进行处理；之后采用独立成分分析（ICA）提取数据的独立成分信息；最后利用单类支持向量机（OCSVM）进行建模构造监控统计量，并利用核密度估计方法确定控制限。“在线监测”包括：对新采集的数据按照模型进行处理，计算其统计量并与控制限进行比较判断发酵过程是否运行正常。本发明无需假设发酵过程变量服从高斯或非高斯等具体的分布，故障监测的准确率较高。

Description

一种基于MICA-OCSVM的发酵过程故障监测方法

技术领域

本发明涉及基于数据驱动的故障诊断技术领域，特别是涉及一种针对间歇过程的故障诊断技术。本发明的基于数据驱动的方法即是在典型间歇过程——青霉素发酵过程故障监测方面的具体应用。

背景技术

近几十年来，间歇过程因其可以满足生产高附加值产品的需求而受到广泛关注。但其机理复杂、操作复杂度高、产品质量易受不确定性因素的影响。青霉素发酵过程作为典型的间歇过程，具有较强的非线性、动态性、混合高斯分布等特点，为了保证发酵过程操作系统的安全性和稳定性，建立一种有效的过程监控方案来及时地检测异常现象是十分有必要的。

目前，多元统计技术已经广泛应用到间歇过程（包括发酵过程）的过程监控。其中应用较多的方法有多向主成分分析(Multi-way Principal Component Analysis，MPCA)和多向偏最小二乘分析(Multi-way Partial Least Squares，MPLS)，这些方法需要假设过程变量服从或近似服从多元高斯分布，应用在实际间歇过程监控存在着不足。MICA可以利用高阶统计信息提取过程数据中的非高斯信号，显现出比MPCA更优的监控效果。然而，实际的工业过程多为高斯和 非高斯的混合分布，因此以上方法因需假设过程变量服从具体的分布而受到应用限制。一些学者采用PCA和ICA分步建模来解决这一问题，提出了组合方法MICA-PCA对非高斯信息和高斯信息分别监控，但该方法不能明确选择主要独立成分的个数，无法做到有效区分非高斯信息和高斯信息。

发明内容

为了克服上述现有技术的不足，本发明提供了一种基于MICA-OCSVM的发酵过程故障监测方法。将MICA提取的独立成分用于OCSVM建模并构造非线性的监控统计量。OCSVM是基于统计学习理论提出的，因此其无需假设过程变量服从具体的分布。而且OCSVM可以确定正常工况下潜隐变量的非线性边界，有效减少过程监测中误报、漏报的发生，提高故障监测的准确性。

本发明采用了如下的技术方案及实现步骤：

A.离线建模阶段：

1）采集发酵过程正常工况下的历史数据，所述的历史数据X由离线测试得到的同一发酵过程相同工艺下的I批次数据构成，X=(X₁,X₂,...,X_I)^T，其中X_i表示第i批次数据，每个批次包含K个采样时刻，每个采样时刻采集J个过程变量，即X_i=(X_i,1,X_i,2,...,X_i,K)，其中X_i,k表示第i批次第k采样时刻采集的数据，X_i,k=(x_i,k,1,x_i,k,2,...,x_i,k,J)，其中x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值；

2）对历史数据进行标准化处理，处理方式如下：

首先计算历史数据X的所有时刻上所有过程变量的均值和标准 方差，其中第k采样时刻的第j个过程变量的均值x_k,j的计算公式为，x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值，k=1,...,K，j=1,...,J；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为， $s_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{i,k.j}-{\stackrel{\‾}{x}}_{k,j})}^2},k=1,...,K,j=1,...,J.$

然后对历史数据X进行标准化，其中第i批次中第k采样时刻的第j个过程变量的标准化计算公式如下：

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

其中，i=1,...,I，j=1,...,J，k=1,...,K；

3）将步骤2）标准化后的数据重新构造成二维矩阵X'，该矩阵共有J个列向量，即X'=(X₁',X'₂,...,X'_J)，其中第j个列向量X'_j=(X'_j,1,...,X'_j,K)^T，X'_j,k=(X'_j,k,1,...,X'_j,k,I)^T，其中X'_j,k,i表示经过步骤2）标准化处理后的第j个过程变量第k个采样时刻在第i个批次中对应的值，其中i=1,...,I，j=1,...,J，k=1,...,K；

4）利用独立成分分析（ICA）方法提取独立成分：

4.1）对二维矩阵X'利用主成分分析（PCA）方法得到白化矩阵Q，Q=Λ^-1/2U^T，其中U和Λ分别为X'的协方差矩阵对应的特征向量矩阵和特征值矩阵，之后将X'进行白化，白化公式为：Z=QX'；

4.2）利用快速ICA算法（Fast ICA）从Z中估计出新的混合矩阵B和分离矩阵W，再根据S=B^TZ得到独立成分S，S=(S₁',...,S'_K)，S'_k=(S'_k,1,...,S'_k,I)，S'_k,i=(S′_1,k,i,...,S'_J,k,i)^T，其中i=1,...,I，k=1,...,K；

5）利用单类支持向量机（OCSVM）对所有S'_k分别进行训练， k=1,...,K，依次得到S在每个采样时刻的决策超平面函数，其中F_k(·)表示利用OCSVM对S'_k进行训练得到的第k采样时刻的决策超平面函数，然后计算建模数据相应的D统计量，D=(D₁',...,D'_K)，D'_k=(d'_k,1,...,d'_k,I)，其中，d'_k,i=-F_k(S'_k,i)；

6）利用核密度估计方法估计上述求得的D统计量在预设置信限时的估计值，并将其作为模型的控制限；

B.在线监测阶段：

7）采集当前发酵过程第k采样时刻的J个过程变量的数据x_k，并根据步骤2）中得到的k时刻的均值和标准方差对其进行标准化得到其中第k采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

其中，x_k,j为当前第k采样时刻所采集发酵数据中的第j个过程变量，为第k采样时刻的第j个过程变量的平均值，s_k,j为第k采样时刻的第j个过程变量的标准方差，j=1,...,J，k=1,...,K；

8）提取标准化后的k时刻采集数据的独立成分s_k，计算公式如下：

$s_k=W{\tilde{x}}_k$

其中，W为离线建模阶段步骤4.2）中所确定的分离矩阵；

9）计算当前发酵过程k时刻采集数据的监控统计量计算公式如下：

${\stackrel{\‾}{D}}_k=-F_k\left(s_k\right)$

其中，F_k为离线建模阶段步骤5）所确定的第k时刻的决策超平面函数，s_k为标准化后的k时刻采集数据的独立成分；

10）将上述计算得到的监控统计量与建模阶段的步骤6）确定的控制限进行比较，如果超限则认为发生故障，进行报警；否则即为正常。

11）若发酵过程完毕，则终止监测；否则采集下一时刻的数据，返回步骤7），继续进行过程监测。

有益效果

与现有技术相比，本发明将MICA提取的独立成分直接用于OCSVM建模并构造非线性的监控统计量。OCSVM是基于统计学习理论提出的，因此其无需假设过程变量服从具体的分布。而且OCSVM可以确定正常工况下潜隐变量的非线性边界，有效利用独立成分的结构信息。本发明方法可以减少过程监测中误报、漏报的发生，提高故障监测的准确性。

附图说明

图1为本发明的方法流程图；

图2为正常数据X_{30 × 10 × 400}组成形式示意图；

图3为标准化后重新构造的二维矩阵X'组成形式示意图；

图4为提取的独立成分S的组成形式示意图；

图5（a）为现有MICA方法对正常批次数据的独立成分平方和I²监测图；

图5（b）为现有MICA方法对正常批次数据的预测误差平方和SPE监测图；

图6为本发明方法对正常批次数据的故障监测效果图；

图7（a）为现有MICA方法对阶跃故障批次数据的独立成分平方和I²监测效果图；

图7（b）为现有MICA方法对阶跃故障批次数据的预测误差平方和SPE监测图；

图8为本发明方法对阶跃故障批次数据的故障监测效果图；

图9（a）为现有MICA方法对斜坡故障批次数据的独立成分平方和I²监测效果图；

图9（b）为现有MICA方法对斜坡故障批次数据的预测误差平方和SPE监测图；

图10为本发明方法对斜坡故障批次数据的故障监测效果图。

其中，图1中的虚线表示“在线监测”与“离线建模”各步骤之间存在的联系。

具体实施方式

青霉素是一种高效、低毒、临床应用广泛的重要抗生素，其生产过程是一个典型的动态、非线性、多阶段的间歇生产过程。美国Illinois州立理工学院过程监控与技术小组开发的青霉素仿真平台PenSim2.0，为青霉素间歇生产过程的监测、故障诊断以及控制提供了一个标准平台。在此平台上可以实现青霉素发酵过程的一系列仿真，相关研究已表明该仿真平台的实用性与有效性，其已经成为国际 上较有影响的青霉素仿真平台。

本实验即以PenSim2.0为仿真研究对象，设置采样时间间隔为1h，选取10个过程变量监视过程运行状况，如表1所示。仿真了31批正常数据，2批故障数据。其中选择30批正常数据X_{30 × 10 × 400}用于建立模型，另外1批正常数据和2批故障数据作为测试数据，用于验证方法的有效性。设置的2种故障类型、幅度、起止时间，见表2。

表1建立模型所用变量

表2故障设置情况

将本发明方法应用到上述发酵过程仿真对象包括离线建模和在线监测两个大步骤，具体陈述如下：

A.离线建模阶段：

步骤1：将上述30批正常数据X_{30 × 10 × 400}沿批次方向展开，具体形 式见图2。可以看到一共有400个矩形框（即400个时刻），每个矩形框都由30行10列构成（即30个批次，10个变量）；

步骤2：对X_{30 × 10 × 400}进行标准化处理。首先按公式计算第k采样时刻的第j个过程变量在所有批次上的均值，其中x_i,k,j为X_{30 × 10 × 400}第i批次中第k采样时刻的第j个过程变量的测量值，k=1,...,400，j=1,...,10；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为， $s_{k,j}=\sqrt{\frac{1}{30-1}\underset{i=1}{\overset{30}{\mathrm{\Σ}}}{(x_{i,k.j}-{\stackrel{\‾}{x}}_{k,j})}^2},k=1,...,400,j=1,...,10.$

然后对X_{30 × 10 × 400}进行标准化，其中第i批次中第k采样时刻的第j个过程变量的标准化计算公式如下：

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

其中，i=1,...,30，j=1,...,10，k=1,...,400；

步骤3：将步骤2）标准化后的数据重新构造成二维矩阵X'，见图3，可以看到一共有400个矩形框（即400个时刻），每个矩形框都由30行10列构成（即30个批次，10个变量）；

步骤4：利用独立成分分析（ICA）方法提取独立成分：

4.1：利用Matlab中的主成分分析函数（princomp）对二维矩阵X'进行分解，得到X'的协方差矩阵对应的特征向量矩阵U和特征值矩阵Λ：[U,Tr,Λ]=princomp(X')，其中的得分矩阵Tr没有使用；之后构造白化矩阵Q，Q=Λ^-1/2U^T；最后将X'进行白化，白化公式为：Z=QX'；

4.2：利用Matlab中的Fast ICA工具箱，从Z中估计出新的 混合矩阵B和分离矩阵W，再根据S=B^TZ得到独立成分S，S的具体形式见图4，可以看到一共有400个矩形框（即400个时刻），每个矩形框都由10行30列构成（即10个独立成分，30个批次）；

步骤5：利用Matlab中libsvm工具箱的“svmtrain”函数，对每个时刻的独立成分分别进行训练，依次得到S在第k采样时刻的决策超平面函数F_k(·)，k=1,...,400，然后使用“svmpredict”函数计算建模数据相应的D统计量，D=(D₁',...,D'_K)，D'_k=(d'_k,1,...,d'_k,I)，其中，d'_k,i=-F_k(S'_k,i)，i=1,...,30，S'_k,i为S在第k采样时刻的第i列；

步骤6：利用Matlab中的核密度估计函数“ksdensity”估计上述求得的D统计量在置信限为0.99时的值，并将其作为模型的控制限；

B.在线监测阶段：

步骤7：采集当前发酵过程第k采样时刻的10个过程变量的数据x_k，并根据步骤2中得到的k时刻的均值和标准方差对其进行标准化得到其中第k采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

其中，x_k,j为当前第k采样时刻所采集发酵数据中的第j个过程变量，为第k采样时刻的第j个过程变量的平均值，s_k,j为第k采样时刻的第j个过程变量的标准方差，j=1,...,10；

步骤8：提取标准化后的k时刻采集数据的独立成分s_k，计算公式如下：

$s_k=W{\tilde{x}}_k$

其中，W为离线建模阶段步骤4.2中所确定的分离矩阵；

步骤9：利用“svmpredict”函数按照步骤5中“svmtrain”训练得到的第k时刻的模型，计算当前发酵过程s_k的监控统计量D_k，即：

${\stackrel{\‾}{D}}_k=-F_k\left(s_k\right)$

其中，F_k为离线建模阶段步骤5所确定的第k时刻的决策超平面函数，s_k为标准化后的k时刻采集数据的独立成分；

步骤10：将上述计算得到的监控统计量与建模阶段的步骤6确定的控制限进行比较，如果超限则认为发生故障，进行报警；否则即为正常。

步骤11：若发酵过程完毕，则终止监测；否则采集下一时刻的数据，返回步骤7，继续进行过程监测。

上述步骤即为本发明方法在青霉素发酵仿真平台故障监测领域的具体应用。为了验证本方法的有效性，对3批测试数据分别进行了在线监测阶段的实验。得到的实验结果见图5至图10，每幅图分别包括与横坐标平行的线和曲线，其中与横坐标平行的线为通过核密度估计方法确定的控制限，曲线为实时的监测值。如果曲线的值大于控制限的值，说明在此时刻发酵过程发生了故障；否则说明发酵过程运行正常。

图5和图6分别为现有MICA方法和本发明方法对正常批次数据的监测效果图。图5（a）中与横坐标平行的线为控制限，曲线为实时的I²监测值；图5（b）中与横坐标平行的线为控制限，曲线为实 时的SPE监测值；图6中与横坐标平行的线为控制限，曲线为实时的D监测值。可以发现，图5（a）中的I²监测图在第3、30、39时刻发生了3次超限，为误报警；图5（b）中的SPE监测图在第108、139、167、315、355、371时刻发生了6次误报警。而本发明方法的D监测图没有发生任何误报警，监测效果较好。图7和图8分别为现有MICA方法和本发明方法对阶跃故障批次数据的监测效果图。图7（a）中的横线为控制限，曲线为实时的I²监测值；图7（b）中的横为控制限，曲线为实时的SPE监测值；图8中的横线为控制限，曲线为实时的D监测值。由于故障1为较大的阶跃型故障，两种监测方法均可有效的检测出故障的发生。但是现有方法依然存在较多的误报警，比如图7（a）中的I²统计量在第1、4、7、12、18、19、28、30、32、34、37、39时刻均发生了误报警，而本发明方法没有任何误报警，显示出更好的监测效果。图9和图10分别为现有MICA方法和本发明方法对斜坡故障批次数据的监测效果图。图9（a）中的横线为控制限，曲线为实时的I²监测值；图9（b）中的横线为控制限，曲线为实时的SPE监测值；图10中的横线为控制限，曲线为实时的D监测值。由于此故障为缓慢变化的斜坡故障，两种监测方法均不能及时检测出在第200时刻发生的故障，有一定的延时性。其中，图9（a）中的I²监测图在第215时刻开始检测到故障，延时性较长；图9（b）中的SPE监测图和图10中的D监测图均在第211时刻开始检测到故障。相比之下，本发明方法在故障检测的及时性方面略优于现有MICA方法。

为了更形象地对比现有方法与本发明方法应用于发酵过程故障 监测的有效性，对上述三批测试数据的监测效果列表对比如下：

表3故障检测情况

注：1)未发生故障时不产生警报的样本数目加上故障发生后报警的样本数目总和与总的样本数目比值；2)未发生故障时产生报警的样本与总的样本数目比值;3)发生故障后未报警的样本与总的样本数目的比值 。

从上表3不难发现，无论是在误报率、漏报率方面还是在准确率方面，本发明方法的统计量均比现有方法的统计量有所提升，提高了青霉素发酵过程故障监测效果。

Claims (2)

1.一种基于MICA-OCSVM的发酵过程故障监测方法，其特征包括“离线建模”和“在线监测”两个阶段，具体步骤如下：

A.离线建模阶段：

1)采集发酵过程正常工况下的历史数据，所述的历史数据X由离线测试得到的同一发酵过程相同工艺下的I批次数据构成，X＝(X₁,X₂,...,X_I)^T，其中X_i表示第i批次数据，每个批次包含K个采样时刻，每个采样时刻采集J个过程变量，即X_i＝(X_i,1,X_i,2,...,X_i,K)，其中X_i,k表示第i批次第k采样时刻采集的数据，X_i,k＝(x_i,k,1,x_i,k,2,...,x_i,k,J)，其中x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值；

2)对历史数据X进行标准化处理，处理方式如下：

首先计算历史数据X的所有时刻上所有过程变量的均值和标准方差，其中第k采样时刻的第j个过程变量的均值的计算公式为，x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为，

然后对历史数据X进行标准化，其中第i批次中第k采样时刻的第j个过程变量的标准化计算公式如下：

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

3)将步骤2)标准化后的数据重新构造成二维矩阵X'，该矩阵共有J个列向量，即X'＝(X′₁,X'₂,...,X'_J)，其中第j个列向量X'_j＝(X'_j,1,...,X'_j,K)^T，X'_j,k＝(X′_j,k,1,...,X'_j,k,I)^T，其中X'_j,k,i表示经过步骤2)标准化处理后的第j个过程变量第k个采样时刻在第i个批次中对应的值；

4)利用独立成分分析ICA方法提取独立成分：

4.1)对二维矩阵X'利用主成分分析PCA方法得到白化矩阵Q，Q＝Λ^-1/2U^T，其中U和Λ分别为X'的协方差矩阵对应的特征向量矩阵和特征值矩阵，之后将X'进行白化，白化公式为：Z＝QX'；

4.2)利用快速ICA算法Fast ICA从Z中估计出新的混合矩阵B和分离矩阵W，再根据S＝B^TZ得到独立成分S，S＝(S₁',...,S'_K)，S'_k＝(S'_k,1,...,S'_k,I)，S'_k,i＝(S₁'_,k,i,...,S'_J,k,i)^T，其中i＝1,...,I，I表示批次数据；

5)利用单类支持向量机OCSVM对所有S'_k分别进行训练，k＝1,...,K，依次得到S在每个采样时刻的决策超平面函数，其中F_k(·)表示利用OCSVM对S'_k进行训练得到的第k采样时刻的决策超平面函数，然后计算建模数据相应的D统计量，D＝(D₁',...,D'_K)，D'_k＝(d'_k,1,...,d'_k,I)，其中，d'_k,i＝-F_k(S'_k,i)；

6)利用核密度估计方法估计上述求得的D统计量在预设置信限时的估计值，并将其作为模型的控制限；

B.在线监测阶段：

7)采集当前发酵过程第k采样时刻的J个过程变量的数据x_k，并根据步骤2)中得到的k时刻的均值和标准方差对其进行标准化得到其中第k采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

其中，x_k,j为当前第k采样时刻所采集发酵数据中的第j个过程变量，为第k采样时刻的第j个过程变量的平均值，s_k,j为第k采样时刻的第j个过程变量的标准方差，j＝1,...,J，k＝1,...,K；

8)提取标准化后的k时刻采集数据的独立成分s_k，计算公式如下：

$s_k=W{\tilde{x}}_k$

其中，W为离线建模阶段步骤4.2)中所确定的分离矩阵；

9)计算当前发酵过程k时刻采集数据的监控统计量计算公式如下：

${\stackrel{\‾}{D}}_k=-F_k\left(s_k\right)$

其中，F_k为离线建模阶段步骤5)所确定的第k时刻的决策超平面函数，s_k为标准化后的k时刻采集数据的独立成分；

10)将上述计算得到的监控统计量与建模阶段的步骤6)确定的控制限进行比较，如果超限则认为发生故障，进行报警；否则即为正常；

11)若发酵过程完毕，则终止监测；否则采集下一时刻的数据，返回步骤7)，继续进行过程监测。

2.根据权利要求1所述的一种基于MICA-OCSVM的发酵过程故障监测方法，其特征在于：步骤5)所述的决策超平面函数F_k(·)计算公式如下：

F_k(·)＝w_kφ_k(·)+b_k

其中w_k、φ_k(·)、b_k为使用OCSVM对S'_k进行训练时得到的参数。