CN101158693B - 基于多核独立元分析的批量生产过程故障检测方法 - Google Patents

Abstract

一种基于多核独立元分析的批量生产过程故障检测方法，包括采集数据、数据处理、利用核主元分析对数据进行白化、利用修正ICA提取独立元及利用T²和SPE统计量进行故障检测步骤。本发明首次提出一种基于MKICA的批量生产过程监测方法，实现对复杂过程进行监测。能较早的检测到故障，从而减少工业生产过程中的损失。

Description

基于多核独立元分析的批量生产过程故障检测方法

技术领域

本发明属于故障诊断技术领域，特别涉及一种基于多核独立元分析MKICA(MultiwayKernel Independent Component Analysis)的批量生产过程监测方法。

背景技术

随着工业规模的日益扩大，生产产品数量的增加，越来越多的工业领域使用批量生产方式生产。比如啤酒行业、味精行业和制药行业等许多领域。尤其是在复杂非线性工业生产过程中，批量生产方式应用的更为普遍。由于非线性生产过程中决定产品质量的因素较多，对条件要求比较苛刻，所以条件的变化大都会降低产品的质量，比如发酵过程。发酵过程是一个非常复杂的非线性过程，有关过程操作条件优化方面的报道不多。对于一个发酵装置、菌体及其相应的发酵培养基均确定的发酵过程而言，影响发酵过程的条件是pH、温度及溶解氧等要素。对于像生产木糖醇这类依靠氧化还原酶实现生物转化的过程来说，pH和溶解氧的因素尤为重要。

在许多采用批量方式生产的化学制品、医药制品、生物制品和其它产品制造行业中，生产过程中关键时期的微小变动都会降低产品质量。因此应用有效的故障检测方法能在生产过程中较早的检测到故障，可以减少不合格批次的数量，明显提高产品质量。批量生产过程的故障主要发生在进料泵和反应装置等部位。

目前，基于多元统计分析的几项技术被广泛应用于批量生产的故障检测。通过扩展多元统计过程控制方法，MacGregor在《Technometrics》等期刊中、Gallagher在《Comput.Chem.Eng.》中分别提出多主元分析方法。在《IEEE Transactions on semiconductor manufacturing》和《J.of Process Control》中，Qin提出了基于递归最小二乘方法的多偏最小二乘法，这种方法可用于监测过程数据和产品质量。在《Chem.Intell.Lab.Sys.》中，Nomikos和MacGregor提出了对批量生产过程的丢失数据的估计方法。这些方法在实际中很起作用。然而，如果应用不需要补充丢失数据的多偏最小二乘监测方法会更好。在《AIChE J.》中，Dong和McAvoy使用基于主曲线和神经网络相结合的非线性主元分析法来监测批量生产。非线性主元分析法是基于神经网络的。由于神经网络需要离线和在线学习近似的动态变化，所以是很慢的，不适用于实际的过程监测中。在《Chem.Intell.Lab.Sys.》中，R

nnar等提出一个有适应性的批量生产过程监测方法。这种方法使用等级主元分析方法，克服了在主元分析方法中对丢失数据的需要。最近，核主元分析方法也用来提取非线性关系。已经证实这种方法能通过分析历史数据很有效的检测故障，同时在线监测新批次的效果也很好。但是以上提出的方法应用在非高斯性的工业过程，如发酵过程中表现不很理想。

发明内容

本发明首次提出基于MKICA的批量生产过程的故障检测方法。针对批量生产过程所具有的非线性特征，提出一种基于多核独立元分析方法的故障检测方法，在标准操作条件模式下建模，从在线监测得到的实时工况数据中提取独立元，利用统计变量T²和SPE进行故障检测。此方法与基于多独立元分析(Multiway Independent Component Analysis，MICA)的方法比较，证实了基于多核独立元分析方法的过程监测方法能较早的检测到故障，从而提高产品的质量和生产率。

使用多核独立元分析方法进行批量生产过程的故障检测方法的详细步骤如下：

步骤一：数据采集

采集观测变量数据。对于每个故障，采集两组数据，即训练数据和实时工况数据。其中，训练数据用于建立模型，实时工况数据用于故障检测。例如在发酵过程中，进料泵的故障导致耗氧率的改变，需要采集温度、压力和PH值等相关变量进行故障检测。

步骤二：数据处理

在数据采集过程中，会丢失一些数据。对于丢失数据的补充有以下几个方法：

①用零补充丢失的数据；

②用当前值补充丢失的数据；

③用平均值补充丢失的数据；

④用主元分析方法补充丢失的数据；

⑤用变量之间冗余的关系补充丢失的数据。

依据批量生产的特征选择最合适的方法，方法②和③应用的比较普遍。本发明采用方法②。

批量生产过程数据是关于批次、变量和时间的数据。首先把数据放入三维矩阵X(I×J×K)中，其中I是批次，J是变量的数目，K是每个批采样的次数，然后将矩阵X(I×J×K)变换成二维矩阵X(I×JK)，如图1所示。例如，一批JK维的标准操作数据，x_k∈R^JK，k＝1，...，I，x_k为观测数据，R^JK为JK维空间。用均值和标准偏差规范化X(I×JK)。经过转化后，运用MKICA方法就等效于在大的二维矩阵X上运用普通的核独立元分析(Kernel Independent Component Analysis，KICA)方法。

步骤三：用核PCA对数据进行白化处理

白化的作用是减少独立元分析(ICA)方法需要估计的参数，从而简化计算过程。这一步的目的就是通过非线性映射把观测数据从输入空间映射到特征空间中，然后在特征空间进行白化处理。

首先，进行非线性映射。设观测数据x_k∈R^JK，通过非线性映射Ф：R^JK→F，F为特征空间，把原始空间的观测数据映射到高维特征空间，Ф(x_k)∈F。特征空间中协方差矩阵为

$C^F=\frac{1}{I}\underset{k=1}{\overset{I}{\mathrm{\Σ}}}\mathrm{\Φ}\left(x_k\right)\mathrm{\Φ}{\left(x_k\right)}^T$

其中Ф(x_k)，k＝1，...，I，为映射到高维特征空间的数据。假定Ф(x_k)为零均值和单位方差。令Θ＝[Φ(x₁)，...，Ф(x_I)]，因而协方差矩阵C^F能表示为 $C^F=\frac{1}{I}\mathrm{\&Theta;}{\mathrm{\&Theta;}}^T.$ 为了避免在特征空间中执行非线性映射和内积计算，定义了一个I×I维的Gram核矩阵K：[K]_ij＝K_ij＝<Ф(x_i)，Ф(x_j)>＝k(x_i，x_j)    (1)

有K＝Θ^TΘ，k(x_i，x_j)为核函数，x_i和x_j为观测数据，1≤i，j≤I。应用核函数k(x_i，x_j)，可以不通过非线性映射就能在F中计算内积，这样就避免了在特征空间中执行非线性映射和内积计算。常用的核函数是径向基核函数 $k(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{c}),$ 多项式核函数k(x，y)＝<x，y>^r和S形核函数k(x，y)＝tanh(β₀<x，y>+β₁)。核函数的选择决定了映射Ф和特征空间F。由核矩阵K可知，高维空间中Ф(x_k)的中心化能按以下方法执行，即将Ф(x_k)的中心化转化为对K的中心化处理。从下式可以得到中心化核矩阵

$\tilde{k}=K-1_{I}K-K{1}_I+1_{I}K{1}_I---\left(2\right)$

其中，

对进行特征值分解：

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

可以得到

的d个最大正特征值λ₁≥λ₂≥...≥λ_d和相应的标准正交的特征向量α₁，α₂，...，α_d。在本发明中，按经验选取满足 $\frac{{\mathrm{\&lambda;}}_i}{\mathrm{sum}\left({\mathrm{\&lambda;}}_i\right)}>0.001$ 的特征值。那么，协方差矩阵C^F的d个最大正特征值是 $\frac{{\mathrm{\&lambda;}}_1}{I},\frac{{\mathrm{\&lambda;}}_2}{I},...,\frac{{\mathrm{\&lambda;}}_d}{I},$ 相应的标准正交特征向量v₁，v₂，...，v_d可以表示为

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

C^F的特征矩阵V＝[v₁，v₂，...，v_d]可以简单的表示如下：

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

其中Λ＝diag(λ₁，λ₂，...，λ_d)是由d个特征值构成的对角矩阵，H＝[α₁，α₂，...，α_d]是对应的特征向量构成的矩阵。用C^F的特征矩阵V把协方差矩阵C^F化为对角矩阵：

$C^F=\mathrm{Vdiag}(\frac{{\mathrm{\&lambda;}}_1}{I},\frac{{\mathrm{\&lambda;}}_2}{I},...,\frac{{\mathrm{\&lambda;}}_d}{I})V^T=\frac{1}{I}{\mathrm{V\&Lambda;V}}^T---\left(6\right)$

令

$P=V{\left(\frac{1}{I}\mathrm{\&Lambda;}\right)}^{-1/2}=\sqrt{I}\mathrm{\&Theta;H}{\mathrm{\&Lambda;}}^{-1}---\left(7\right)$

那么，

P^TC^FP＝E    (8)

E为单位矩阵。这样，得到了白化矩阵P。

映射到特征空间中的数据Ф(x)经过如下转换可以被白化：

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

z为白化后的观测变量。具体的，

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

$=\sqrt{I}{\mathrm{\&Lambda;}}^{-1}{H}^T{[\tilde{k}(x_1,x),...,\tilde{k}(x_I,x)]}^T---\left(10\right)$

其中，x是需要白化的数据。

步骤四：使用修正ICA提取独立元

从KPCA变换空间中提取独立元。应用Lee和Qin.在  AIChE J.》中提出的修正ICA方法从观测数据中提取独立元。这种方法的优点是不但可以为过程监测提取较少的主要独立元，还可以根据独立元的方差进行排列。

修正ICA方法可以找到p个独立元，y＝{y₁，...，y_p}。为使元素彼此之间能尽可能的独立，y需满足E{yy^T}＝D＝diag{λ₁，...，λ_p}，λ₁，...，λ_p是p个独立元所对应的特征值。设

y＝C^Tz    (11)

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

定义归一化的独立元y_n为：

y_n＝D^-1/2y＝D^-1/2C^Tz＝C_n ^Tz    (12)

C_n为标准得分转换矩阵。显然，D^-1/2C^T＝C_n ^T，C_n ^TC_n＝E，E{y_ny_n ^T}＝E。这样一旦找到C_n，由式(11)和(12)就可以得到独立元y。所以下面就来求解C_n。由于z是观测数据经过白化后得到的，具有不相关性，所以可以选择z的第一个数据中p个独立元作为y_n的初始值。令C_n ^T的初始矩阵设置为：

C_n ^T＝[I_p

0]        (13)

其中I_p是p维单位矩阵，0是p×(d-p)维零矩阵。

用修正ICA算法计算标准得分转换矩阵C_n。C_n的详细的算法如下：

1)选择要估计的独立元的个数p。设置计数器i←1；

2)取初始向量c_n，i；

3)最大化近似负熵：使c_n，i←E{zg(c_n，i ^Tz)}-E{g′(c_n，i ^Tz)}c_n，i，其中g是G的一阶导数，g′是G的二阶导数。Hyv

rinen给出了三种g的表达函数：g₁(u)＝tanh(a₁u)，g₂(u)＝u exp(-a₂u²/2)，g₃(u)＝u³；

4)为了排除已包含的信息做正交化： $c_{n,i}\&LeftArrow;c_{n,i}-\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}\left({c_{n,i}}^T{c}_{n,j}\right)c_{n,j};$

5)归一化 $c_{n,i}\&LeftArrow;\frac{c_{n,i}}{\left|\right|c_{n,i}\left|\right|};$

6)如果c_n，i没有收敛，返回第3)步；

7)如果c_n，i收敛了，输出向量c_n，i。如果i≤p，则i←i+1，返回2)。这样就找到C_n。

因此独立元可由下式获得：

y＝D^1/2C_n ^Tz＝D^1/2C_n ^TP^TФ(x)    (14)

得到的目标函数y_n，i，  y_n，i＝(c_n，i)^Tz，i＝1，...，p是统计独立的，具有最大化的非高斯性。

步骤五：利用统计变量T²和SPE进行故障检测。

通过统计变量T²和SPE的分布可以检测到是否有故障发生。当观测数据的统计量超出统计量规定的控制限时属于异常数据，表明发生了故障。

T²统计定义如下：

T²＝y^TD^-1y， $y=D^{1/2}{C}_n^T{P}^T\mathrm{\&Phi;}\left(x\right)---\left(15\right)$

SPE统计定义如下：

$\mathrm{SPE}=e^{T}e={(z-\hat{z})}^T(z-\hat{z})=z^T(E-C_n{C_n}^T)z,$ z＝PTФ(x)    (16)

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

可由下式得到：

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

因为y不服从高斯分布，所以T²的控制上界可由F分布确定。

如果大部分非高斯性包含在提取的独立元中，剩余的子空间将包含大部分的高斯噪声，这个噪声可以认为是正态分布的。假设预测误差是正态分布的，SPE的控制界限可由下式加权x²分布计算：

SPE～μχ_h ²

μ＝b/2a，h＝2a²/b    (18)

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

具体算法如下：

本发明方法的实现过程包括标准操作条件模式的确定和在线监测两部分。

1.离线分析(标准操作条件模式的确定)

(1)在正常的生产中得到观测变量；

(2)用当前值补充丢失的数据，把观测变量矩阵从X(I×J×K)变化到X(I×JK)；

(3)用均值和标准偏差对X(I×JK)进行规范化；

(4)用等式(1)计算核矩阵K∈R^I×I。用等式(2)计算中心化核矩阵

用式(3)对进行特征值分解，得到与最大正特征值λ₁≥λ₂≥...≥λ_d对应的的正交特征向量(α₁，α₂，...，α_d)；

(5)根据公式(7)计算白化矩阵P；

(6)根据修正ICA中给出的算法计算C_n；

2.在线监测

(1)用当前值补充丢失的数据。把数据X_t(K×J)，伸展为X_t ^T(1×JK)；(2)由[k_t]_k＝[k_t(x_t，x_k)]计算核向量k_t∈R^1×I。其中，x_k是被规范化的训练数据，x_k∈R^JK，k＝1，2，...，I；

(3)中心化k_t如下：

${\tilde{k}}_t=k-1_{t}K-K{1}_I+1_{t}K{1}_I---\left(17\right)$

其中，K是由标准操作条件模式的第(4)步确定的。1_t＝(1/I)[1，...，1]∈R^1×I，

(4)对于实时工况数据x_t，通过下式计算白化后的观测变量z_t

$z_t=P^T\mathrm{\&Phi;}\left(x_t\right)=\sqrt{I}{\mathrm{\&Delta;}}^{-1}{H}^T[\tilde{k}(x_1,x_t),...,\tilde{k}(x_1,x_t){]}^T---\left(20\right)$

(5)根据标准操作条件模式计算的C_n，通过下式计算独立元y_t

y_t＝D^1/2C_n ^Tz_t    (21)

其中，D＝diag{λ₁，...，λ_p}。

(6)计算实时工况数据的监测统计量(T²和SPE)；

本发明方法具有以下优势：

1、本发明首次提出基于MKICA的批量生产过程的故障检测方法。

2、本发明的控制方法能较早的检测到故障，从而减少工业生产过程中的损失。

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

附图说明

图1三维数据转换为二维数据；

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

图3.从10小时引进故障的青霉素监测过程a)MICA监测结果；b)MKICA监测结果；

图4.从1小时引进故障的诺西肽监测过程a)MICA方法的监测结果，b)MKICA方法的监测结果；

图5.从45小时引进故障的诺西肽监测过程a)MICA方法的监测结果，b)MKICA方法的监测结果；

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

具体实施方式

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

实例1.

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

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

No.

变量

1	通风率(l/h)
		2	煽动功率(W)
3	培养基供给率
		4	培养基温度
5	溶解氧的浓度
		6	培养容积(l)
7	CO2浓度
		8	pH值
9	产生的热量(kcal)

步骤一：数据采集

批15的训练数据由400组数据构成，  实时工况数据由800组数据构成。训练数据和实时工况数据都包含了9个监测变量。本实例主要针对霉的供给速率故障进行实例分析，对故障批15随机选取了10组训练数据和10组实时工况数据，所监测的变量数据分别如表2和表3所示。

表2训练数据

通风率(l/h)	煽动功率(W)	培养基供给率	培养基温度	溶解氧的浓度
					8.6015	3.04E+01	0.042325	295.9975	96.3164
8.6025	3.04E+01	0.0423	295.9975	96.0497
					8.6045	3.04E+01	0.0423	295.995	96.2852
8.6065	3.04E+01	0.042275	295.9925	96.3152
					8.609	30.36	0.042225	295.9875	96.3498
8.612	3.03E+01	0.042225	295.98	96.4273
					8.6145	3.03E+01	0.042275	295.97	96.0804
8.6165	3.03E+01	0.04235	295.965	96.0652
					8.6175	3.03E+01	0.0424	295.965	96.6134

8.6175

3.03E+01

0.04245

295.965

96.7248

培养容积(l)	CO2浓度	pH值	产生的热量(kcal)
				99.808	1.8438	5.1092	52.7822
99.8255	1.7173	5.0712	53.0277
				99.8433	1.8334	5.0366	53.2699

99.861	1.8646	5.018	53.5109
				99.8781	1.7888	5.0003	53.7506
99.8954	1.8471	5.2167	53.9915
				99.9131	1.7185	5.1698	54.2273
99.9309	1.8684	5.1279	54.4581
				99.9483	1.7137	5.0898	54.687
99.9661	1.7321	5.0551	54.9139

表3实时工况数据

通风率(l/h)	煽动功率(W)	培养基供给率	培养基温度	溶解氧的浓度
					8.6015	30.365	0.042325	295.9975	96.3164
8.6025	30.365	0.0423	295.9975	96.0497
					8.6045	30.37	0.0423	295.995	96.2852
8.6065	30.37	0.042275	295.9925	96.3152
					8.609	30.36	0.042225	295.9875	96.3498
8.612	30.345	0.040114	295.98	96.4273
					8.6145	30.325	0.040161	295.97	96.1555
8.6165	30.3	0.040233	295.965	96.1635

8.6175	30.28	0.04028	295.965	96.719
					8.6175	30.25	0.040328	295.965	96.8329

培养容积(1)	CO2浓度	pH值	产生的热量(kcal)
				99.808	1.8438	5.1092	52.7822
99.8255	1.7173	5.0712	53.0277
				99.8433	1.8334	5.0366	53.2699
99.861	1.8646	5.018	53.5109
				99.8781	1.7888	5.0003	53.7506
99.8954	1.8471	5.2167	53.9915
				99.911	1.7181	5.1714	54.2117
99.9267	1.8673	5.1325	54.4109
				99.942	1.7117	5.0973	54.6038
99.9577	1.7293	5.065	54.7938

步骤二：数据处理

把训练数据放入三维矩阵X(I×J×K)中，其中I是批次，J是变量的数目，K是每个批采样的次数。变换X(I×J×K)到X(I×JK)，如图1所示。用均值和标准偏差规范化X(I×JK)，得到训练数据x_k∈R^JK，k＝1，...，I；把实时工况数据X_t(k×J)，伸展为X_t ^T(1×Jk)，得到实时工况数据x_t∈R^JK，k＝1，2，...，I。

步骤三：用KPCA对数据进行白化处理

(1)利用训练数据计算白化矩阵。选取径向基核函数作为核函数，把数据x_k∈R^JK代入公式(1)计算核函数，得到核矩阵K。接着根据公式(2)对Gram核矩阵K进行中心化处理，得到中心化核矩阵根据公式(3)计算出了

的特征值，根据经验选取了5个最大的正特征值为λ₁≥λ₂≥...≥λ₅。从而得到协方差矩阵C^F的5个最大正特征值是 $\frac{{\mathrm{\&lambda;}}_1}{I},\frac{{\mathrm{\&lambda;}}_2}{I},...,\frac{{\mathrm{\&lambda;}}_5}{I},$ 并且根据公式(4)得到相关的标准正交特征向量v₁，v₂，...，v₅。  由公式(5)-(8)得到白化矩阵P。

(2)对实时工况数据进行白化处理。把数据x_t∈R^JK，代入式[k_t]_k＝[k_t(x_t，x_k)]，k＝1，2，...，I中计算核向量k_t∈R^1×I。根据式(19)和由训练数据得到的核矩阵K中心化核向量kt。中心化后的核向量为由训练数据得到的白化矩阵P，根据公式(20)白化映射到特征空间中的监测变量数据x_tk，白化后的观测变量为z_t。

步骤四：使用修正ICA进一步处理z

本步骤利用修正ICA方法找到需要的p个独立元。首先确定要估计的独立元的个数p，取p＝4。根据修正ICA算法计算标准操作条件模式C_n，其中g(u)＝tanh(u)。然后利用上一步得到的白化后的观测变量z_t，根据公式(21)计算得到实时工况数据的独立元y_t。

步骤五：利用统计变量T²和SPE进行故障检测

统计变量T²和SPE可以监测到是否有故障发生。根据公式(15)、(16)计算统计量T2和SPE的控制限。为了检测本方法的效果，将本方法与MICA方法进行对比，其中MICA方法选取的独立元个数为3个。给第一个批加一个故障，就是由进料泵导致的故障使耗氧率从3.5到0.5非线性的减少直到运行结束。故障的起始时间是10小时。监测结果如下：如图3.a)所示，MICA方法的T²图没有检测到故障，MICA方法的SPE图在52小时检测到故障；如图3.b)所示MKICA方法的T²图没有检测到故障，MKICA方法的SPE图在45小时检测到故障，比MICA方法早了7小时。对比可见，MKICA方法能较早的检测到故障，因此是比MICA好的方法

实例2.

本发明提出的方法应用到了诺西肽过程的数据中。诺西肽是一种由链霉菌产生的双环缩氨酸抗生素。它的分子式是C₅₁H₄₃N₁₃Q₁₂S₆。由于诺西肽能促进牲畜成长并且在牲畜体内没有残余，所以它主要用作饲料添加剂。一些植物能产生诺西肽，可以作为添加剂。诺西肽的生产过程是一个需要96个小时的批量需氧发酵过程。在发酵物到达成熟阶段之前，培养菌株的过程是在种子发酵罐中进行的。在到达成熟阶段之后，将成熟发酵物转到末级发酵罐中。为了最优化诺西肽的合成，这些发酵罐中的发酵物是在正常条件下采用批量模式发酵的。在诺西肽的发酵过程中，在线的测量变量包括物理和化学变量，如表4所示。

表4.诺西肽发酵过程变量

No.	变量	单位	描述	最小值	最大值
						1	T	℃	温度	5	50
2	Do	％set.	溶解的氧气量	0	100
						3	P	Pa	压力	0	35
4	O2	％	废气中O2的浓度	12	22
						5	CO2	％	废气中CO2的浓度	0	7
6	pH	pH	pH值	3	10
						7	OUR	mol	O2吸收率	0	5

8

CER

mol

CO2产生速率

0

5

在研究中，8个监测变量中有噪声。为了测试运行过程中是否有故障，对30个标准操作条件模式下的批执行离线分析。这个数据集分别用MICA和MKICA方法分析。一般的，由于KICA选择的独立元是从高维特征空间中提取的，所以KICA选择的独立元数量多于ICA的。通过交叉确认分别为MICA和MKICA模型选择3个和4个独立元。选择径向基核函数 $k(x,y)=\exp (-\frac{{\left|\right|x-y\left|\right|}^2}{c})$ 作为核函数，其中c＝0.5。

接下来用MICA和MKICA方法进行在线监测。基于离线分析，对30个批的实时工况过程应用MICA和KMICA方法进行在线监测。用99％的置信限度分别检验MICA和KMICA方法故障检测的效果。为了检测MICA和MKICA的性能，给第一个批加一个故障，就是由进料泵导致的故障使耗氧率从4.5到0.5非线性的减少直到运行结束。故障的起始时间是1小时。监测结果如图4.a)和4.b)所示。在MICA的情况下，如图4.a)，T²的图没有检测到一个明显的偏差而SPE图能从77小时检测到。与MICA方法相比，MKICA方法的在线监测图较早的检测到了故障。如图4.b)所示，T²图在60小时检测到了故障，而SPE图在30小时就检测到故障。总之，MKICA方法的检测时间比MICA方法早了47小时。第二个批在45小时耗氧率非线性的减少，一直减少到发酵过程结束。对于这个故障批，如图5.a)所示，MICA方法在60小时检测到故障。然而，与MICA方法对比，如图5.b)所示，可见MKICA方法的T²图在55小时显示了一个明显的偏差，表明检测到故障。而MKICA方法的SPE图在60小时检测到故障，比故障发生晚了10小时。这个延迟是耗氧率的影响通过相关变数慢慢传播作用的结果。

这种基于MKICA的批量生产过程故障检测方法可以有效的捕捉变量间的非线性关系，它在过程监测中的应用表现出比MICA更好的性能。这种新方法能较早的检测到故障，从而减少损失，提高产品的质量。本发明的研究成果还可以应用其他生产过程中，例如：啤酒，味精，制药等生产过程中。

Claims (1)

1.一种基于多核独立元分析的批量生产过程故障检测方法，其特征在于将该方法应用于青霉素发酵过程，包括以下步骤：

步骤一、采集数据

采集过程中相关变量的数据，对于每个故障，采集两组数据，即标准操作模式的训练数据和在线监测的实时工况数据；训练数据由400组数据构成，实时工况数据由800组数据构成；训练数据和实时工况数据都包含了9个监测变量，即：通风率、煽动功率、培养基供给率、培养基温度、溶解氧的浓度、培养容积、CO₂浓度、PH值和产生的热量；其中，训练数据用于建立模型，实时工况数据用于故障检测；

步骤二、数据处理

用当前值补充丢失的数据，用均值和标准偏差规范化采集的观测数据，批量生产过程数据是关于批次、变量和时间的数据，首先把数据放入三维矩阵X(I×J×K)中，其中I是批次，J是变量的数目，K是每个批采样的次数，然后将矩阵X(I×J×K)变换成二维矩阵X(I×JK)；

步骤三、利用核主元分析对数据进行白化处理

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

具体过程如下：

首先进行非线性映射，设观测数据x_k∈R^JK，k＝1，...，I，k是观察数据的个数，I是批次，J是变量的数目，K是每个批采样的次数，R^JK为JK维空间，通过非线性映射Φ：R^JK→F，F为特征空间，把原始空间的观测数据映射到高维特征空间中，Φ(x_k)∈F，特征空间中协方差矩阵为

$C^F=\frac{1}{I}\underset{k=1}{\overset{I}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}\left(x_k\right)\mathrm{\&Phi;}{\left(x_k\right)}^T$

其中Φ(x_k)，k＝1，...，I，为映射到高维特征空间的数据，假定Φ(x_k)为零均值和单位方差，令Θ＝[Φ(x₁)，...，Φ(x_I)]，因而协方差矩阵C^F能表示为定义一个I×IGram维的核矩阵K：

[K]_ij＝K_ij＝〈Φ(x_i)，Φ(x_j)〉＝k(x_i，x_j)(1)

有K＝Θ^TΘ，k(x_i，x_j)为核函数，x_i和x_j为观测数据，1≤i，j≤I，应用核函数k(x_i，x_j)，可以不通过非线性映射就能在F中计算内积，这样就避免了在特征空间中执行非线性映射和内积计算，常用的核函数是径向基核函数

多项式核函数k(x，y)＝〈x，y〉^r和S形核函数k(x，y)＝tanh(β₀〈x，y〉+β₁)，核函数的选择决定了映射Φ和特征空间F；由核矩阵K可知，高维空间中Φ(x_k)能按以下方法进行中心化，即将Φ(x_k)的中心化转化为对K的中心化处理，从下式可以得到中心化核矩阵

$\tilde{K}=K-1_{I}K-{K1}_I+1_I{K1}_I---\left(2\right)$

其中，

对

进行特征值分解：

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

可以得到

的d个最大正特征值λ₁≥λ₂≥...≥λ_d和相应的标准正交的特征向量α₁，α₂，...，α_d，按经验选取满足

的特征值，i＝1，...，d，则协方差矩阵C^F的d个最大正特征值是

相应的标准正交特征向量v₁，v₂，...，v_d可以表示为

$v_j=\frac{1}{\sqrt{{\mathrm{\&lambda;}}_j}}\mathrm{\&Theta;}{\mathrm{\&alpha;}}_j,j=1,...,d---\left(4\right)$

C^F的特征矩阵V＝[v₁，v₂，...，v_d]可以简单的表示如下：

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

其中Λ＝diag(λ₁，λ₂，...，λ_d)是由d个特征值构成的对角矩阵，H＝[α₁，α₂，...，α_d]是对应的特征向量构成的矩阵，用协方差矩阵C^F的特征矩阵V把C^F化为对角矩阵：

$C^F=\mathrm{Vdiag}(\frac{{\mathrm{\&lambda;}}_1}{I},\frac{{\mathrm{\&lambda;}}_2}{I},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{{\mathrm{\&lambda;}}_d}{I})V^T=\frac{1}{I}\mathrm{V\&Lambda;}{V}^T---\left(6\right)$

令

$P=V{\left(\frac{1}{I}\mathrm{\&Lambda;}\right)}^{-1/2}=\sqrt{I}\mathrm{\&Theta;H}{\mathrm{\&Lambda;}}^{-1}---\left(7\right)$

那么，

P^TC^FP＝E    (8)

P为白化矩阵，E为单位矩阵，这样就得到了白化矩阵P；

所述的求解白化后的观测变量z的具体过程如下：

映射到特征空间中的数据Φ(x)经过如下转换可以被白化：

z＝P^TΦ(x)(9)

z为白化后的观测变量，

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

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

其中，x是需要白化的数据，

j＝1，2，...，I，为中心化核函数；

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

利用修正ICA算法把白化后的观测变量z转换成独立元，并使独立元各变量之间尽可能相互统计独立；

求解独立元y的具体过程如下：

修正ICA方法可以找到p个独立元，y＝{y₁，...，y_p}，为使元素彼此之间能尽可能的独立，y需满足E{yy^T}＝D＝diag{λ₁，...，λ_p}，λ₁，...，λ_p是p个独立元所对应的特征值，利用

y＝C^Tz    (11)

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

定义归一化的独立元y_n为：

y_n＝D^-1/2y＝D^-1/2C^Tz＝C_n ^Tz    (12)

C_n为标准得分转换矩阵，显然，D^-1/2C^T＝C_n ^T，C_n ^TC_n＝E，E{y_ny_n ^T}＝E，由于z是观测数据经过白化后得到的，具有不相关性，所以可以选择z的第一个数据中p个独立元作为y_n的初始值，因此可以把C_n ^T的初始矩阵设置为：

其中I_p是p维单位矩阵，0是p ×(d-p)维零矩阵；

利用修正ICA算法计算标准得分转换矩阵C_n，由式(11)和(12)就可以得到独立元y，独立元可由下式获得：

y＝D^1/2C_n ^Tz＝D^1/2C_n ^TP^TΦ(x)(14)；

步骤五、利用T²和SPE统计量进行故障检测

采用T²和SPE统计量进行在线故障检测，当观测数据的统计量没有超出统计量规定的控制限时，则属于正常数据，反之属于异常数据，表明出现故障；

批量生产过程统计变量的T²和SPE的计算：

T²统计定义如下：

T²＝y^TD^-1y， $y=D^{1/2}{C}_n^T{P}^T\mathrm{\&Phi;}\left(x\right)---\left(15\right)$

SPE统计定义如下：

$\mathrm{SPE}=e^{T}e={(z-\hat{z})}^T(z-\hat{z})=z^T(E-C_n{C_n}^T)z,z=P^T\mathrm{\&Phi;}\left(x\right)---\left(16\right)$

其中，

可由下式得到：

$\hat{z}=C_n{D}^{-1/2}y={C_{n}C}_n^{T}z---\left(17\right)$

因为y不服从高斯分布，所以T²的控制上界可由F分布确定；

如果大部分非高斯性包含在提取的独立元中，剩余的子空间将包含大部分的高斯噪声，这个噪声可以认为是正态分布的，假设预测误差是正态分布的，SPE的控制界限可由下式加权χ²分布计算：

$\mathrm{SPE}~{\mathrm{\&mu;\&chi;}}_h^2$

μ＝b/2a，h＝2a²/b(18)

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

Images