1 前言

间歇过程作为一种重要的工业生产方式，具有操作性强和生产效率高等优点，在生物制药、食品饮料和精细化工等小批量、高附加值产品的生产制备中得到广泛应用^[1]。强非线性、多模态特性是间歇过程的固有特征，是实际工业过程和在线监控中面临的问题，而寻求一种有效解决途径对降低经济损失和避免人员伤亡具有重要意义^[2]。保障生产过程的安全稳定运行和可靠的生产质量一直是过程监控的重要任务，对质量相关的故障检测和诊断方法提出了更高的要求^[3]。

传统多元统计过程监控方法只能解决线性、单模态的生产过程，间歇过程由于初始条件和不同产品质量的要求不同，批次之间不可避免地出现多模态，为间歇过程多模态问题提出了挑战^[4]。Deng等^[5]通过提出局部近邻相似度分析的方法解决了连续过程多模态特性，但并没有考虑间歇过程；Zhao^[6]提出有限批次的阶段分析建模方法实现了多模态间歇过程在线监控，但没有考虑过程变量与质量变量的关系。偏最小二乘(partial least squares, PLS)是典型的质量变量预测方法，通过容易得到的过程变量来监测难以得到的质量变量，更有助于对产品质量指标的波动进行实时、在线监控和预测^[7]。针对间歇过程监测与质量预测，Yu^[8]提出了一种自适应核偏最小二乘方法；Luo等^[9]基于多元线性回归提出了多线性PLS和滑动窗技术相结合的方法；Wang等^[10]通过一种特征向量的选取方法简化了双核多向偏最小二乘，实现了更高的预测精度。这些方法虽然已成功应用于间歇过程质量预测，但PLS构建的模型主元数多，不利于对高维数据过程监控，并且主元空间依然含有与质量预测无关的质量变量正交成分。Qin等^[11]为了对过程变量和质量变量更加充分地检测与预测，提出了并发潜结构投影(concurrent projection to latent structures，CPLS)模型，但没有考虑实际工业过程的强非线性和间歇过程多模态特性。

为了解决非线性问题，Peng等^[12]提出了全影核偏最小二乘(total kernel partial least squares, TKPLS)方法，通过引入核函数实现线性可分，提高了与质量相关故障的监控效果。近年来，由于高斯混合模型(Gaussian mixture model, GMM)可以更好地处理复杂过程数据而被用于估计复杂数据分布，并成功应用于非线性、多模态特性的生产过程中^[13-15]。本文针对间歇过程的非线性、多模态特性，以及质量相关问题，提出了多向高斯混合模型-并发核熵潜结构投影(multi-way Gaussian mixture model-concurrent kernel entropy projection to latent structures, MGMM-CKEPLS)算法。将非线性采样数据投影到高维核特征空间，实现变量线性化；通过计算Renyi熵贡献的大小对数据降维处理，弥补了PLS选取主元时的不足；根据GMM获得最优高斯分布和统计参数，对多模态间歇过程实现模态聚类；分别对不同模态建立CPLS模型，确定各模态监控统计量并集成统一监控量，实现非线性、多模态过程检测及质量预测。

2 并发潜结构投影(CPLS)

PLS模型通过最大化过程变量和质量变量的协方差进行质量预测，对过程变量方差并没有降序排列提取主元，不利于过程监控，并且单独对得分空间进行监控，忽略了不能被过程变量预测的质量变量。Qin等^[11]提出的CPLS将过程空间和质量空间分成了五个子空间，对每个子空间分别进行监控，实现了对过程空间和质量空间的完整监控。

PLS在低维空间对过程变量和质量变量分解如下：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{X}}{\bf{ = }}\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{P}}^{\rm T}}{\bf{ + }}\mathit{\boldsymbol{E}}\\ \mathit{\boldsymbol{Y}}{\bf{ = }}\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{Q}}^{\mathop{\rm T}\nolimits} }{\bf{ + }}\mathit{\boldsymbol{F}} \end{array} \right. $

(1)

对可以被预测的质量变量$\mathit{\boldsymbol{\hat Y}}{\bf{ = }}\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{Q}}^{\mathop{\rm T}\nolimits} }$分解奇异值：

$ \mathit{\boldsymbol{\hat Y}}{\bf{ = }}{\mathit{\boldsymbol{U}}_{\rm{c}}}{\mathit{\boldsymbol{D}}_{\rm{c}}}\mathit{\boldsymbol{V}}_{\rm{c}}^{\mathop{\rm T}\nolimits} \equiv {\mathit{\boldsymbol{U}}_{\rm{c}}}\mathit{\boldsymbol{Q}}_{\rm{c}}^{\mathop{\rm T}\nolimits} $

(2)

其中，${\mathit{\boldsymbol{Q}}_{\rm{c}}}{\bf{ = }}{\mathit{\boldsymbol{V}}_{\rm{c}}}{\mathit{\boldsymbol{D}}_{\rm{c}}}$由降序排列的非零奇异值和对应的右奇异向量得到，由于$\mathit{\boldsymbol{V}}$为正交阵，则${\mathit{\boldsymbol{U}}_{\rm{c}}}$可由式(3)计算：

$ {\mathit{\boldsymbol{U}}_{\rm{c}}} = \mathit{\boldsymbol{XR}}{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{V}}_{\rm{c}}}\mathit{\boldsymbol{D}}_{\rm{c}}^{ - 1} \equiv \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{R}}_{\rm{c}}} $

(3)

其中，${\mathit{\boldsymbol{R}}_{\rm{c}}} = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{V}}_{\rm{c}}}\mathit{\boldsymbol{D}}_{\rm{c}}^{ - 1}$，T、Q、R由PLS模型得到。

不可被预测的质量变量：

$ {{\mathit{\boldsymbol{\tilde Y}}}_{\rm{c}}}{\bf{ = }}{\mathit{\boldsymbol{T}}_{\rm{y}}}\mathit{\boldsymbol{P}}_{\rm{y}}^{\rm{T}}{\bf{ + }}\mathit{\boldsymbol{\tilde Y}} $

(4)

与质量变量无关的过程变量：

$ {{\mathit{\boldsymbol{\tilde X}}}_{\rm{c}}}{\bf{ = }}{\mathit{\boldsymbol{T}}_{\rm{x}}}\mathit{\boldsymbol{P}}_{\rm{x}}^{\rm{T}}{\bf{ + }}\mathit{\boldsymbol{\tilde X}} $

(5)

因此，构建CPLS模型：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{X}}{\bf{ = }}{\mathit{\boldsymbol{U}}_{\rm{c}}}\mathit{\boldsymbol{R}}_{\rm{c}}^{ - 1}{\bf{ + }}{\mathit{\boldsymbol{T}}_{\rm{x}}}\mathit{\boldsymbol{P}}_{\rm{x}}^{\rm{T}}{\bf{ + }}\mathit{\boldsymbol{\tilde X}}\\ \mathit{\boldsymbol{Y}}{\bf{ = }}{\mathit{\boldsymbol{U}}_{\rm{c}}}\mathit{\boldsymbol{Q}}_{\rm{c}}^{\rm{T}}{\bf{ + }}{\mathit{\boldsymbol{T}}_{\rm{y}}}\mathit{\boldsymbol{P}}_{\rm{y}}^{\rm{T}}{\bf{ + }}\mathit{\boldsymbol{\tilde Y}} \end{array} \right. $

(6)

CPLS将过程变量和质量变量同时投影到五个子空间：联合输入输出的协方差子空间、输出主子空间、输出残差子空间、输入主子空间和输入残差子空间。CPLS只是从全局监控的角度对过程变量和质量变量进行了进一步分解，计算量还是比较复杂，并且只是一种线性处理方法，对于非线性和多模态的间歇过程不能直接应用。

3 多向高斯混合模型-并发核熵潜结构投影(MGMM-CKEPLS) 3.1 间歇过程数据预处理

间歇过程是一种多批次生产过程，其实时监测得到的样本集为三维矩阵：$\mathit{\boldsymbol{X}}\left({I \times J \times K} \right)$，其中，I为批次个数，J为变量个数，K为采样时间。在实际的生产过程中，由于每个批次的采样时间不同等原因而不可能做到完全重复同批次的生产，使得样本数据的批次长度不一致，出现不等长数据问题。本文采用“最短长度法”，以最短的一个批次数据作为标准而将剩余批次数据进行硬性截取，使得所有批次数据都具有相等长度。对三维数据进行建模，需要把已处理成相等长度的数据展开成与连续过程相类似的二维数据。由于批次展开可以降低数据非线性，变量展开可以消除数据预估问题，因此本文结合批次-变量展开方法进行展开，其展开方式如图 1所示。这种展开方式可以消除在线监控时的数据预估问题，提高了在线应用时的监控精度^[16]。

3.2 并发核熵潜结构投影(CKEPLS)

Jenssen^[17]于2010年提出了核熵成分分析(kernel entropy component analysis, KECA)算法用于数据特征的提取和分析，具有主元个数少、运算量低和非线性处理能力强等优点。本文结合KECA对数据处理的优越性和CPLS对质量预测的优点，提出了并发核熵潜结构投影(concurrent kernel entropy projection to latent structures，CKEPLS)算法，减少标准PLS建模时所需的主元数，降低计算量，有效处理非线性过程数据。

假设经过n次采样的m个过程变量和p个质量变量分别构成过程数据$\mathit{\boldsymbol{X}} = {\left[ {{x_1}, \cdots, {x_n}} \right]^{\rm{{\rm T}}}} \in {R^{n \times m}}$和质量数据$\mathit{\boldsymbol{Y}} = {\left[ {{\mathit{y}_1}, \cdots, {\mathit{y}_\mathit{n}}} \right]^{\rm{{\rm T}}}} \in {\mathit{R}^{\mathit{n} \times \mathit{p}}}$。若n维样本X的概率密度函数是p(x)，则其Renyi熵为：

$ H\left( p \right) = - {\rm{lg}}\int {{p^2}\left( x \right){\rm{d}}x} $

(7)

令$V\left(p \right) = \int {{p^2}\left(x \right){\rm{d}}x} $，引入Parzen窗方法对$V\left(p \right)$进行概率密度估计，对核矩阵 K进行特征值分解：

$ \mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \mathit{\boldsymbol{ED}}{\mathit{\boldsymbol{E}}^{\rm{T}}} $

(8)

其中，$\mathit{\boldsymbol{E}} = \left({{\mathit{\boldsymbol{e}}_1}, \cdots, {\mathit{\boldsymbol{e}}_\mathit{n}}} \right)$为特征向量矩阵，$\mathit{\boldsymbol{D}} = {\rm{diag}}\left({{\lambda _1}, \cdots, {\lambda _\mathit{n}}} \right)$为特征值矩阵，通过计算可得：

$ \hat V\left( p \right) = \frac{1}{{{n^2}}}{\left( {\sqrt {{\lambda _i}} \mathit{\boldsymbol{e}}_i^{\rm{{\rm T}}}{\bf{1}}} \right)^2} $

(9)

其中1是每个元素均为1的n维列向量，${{\lambda _i}}$、${\mathit{\boldsymbol{e}}_\mathit{i}}$为核矩阵进行特征值分解得到的特征值和对应的特征向量。根据Renyi熵贡献大小，选取前k个较大的特征值$\mathit{\boldsymbol{D}} = {\rm{diag}}\left({{\lambda _1}, \cdots, {\lambda _\mathit{k}}} \right)$和特征向量$\mathit{\boldsymbol{E}} = \left({{\mathit{\boldsymbol{e}}_1}, \cdots, {\mathit{\boldsymbol{e}}_\mathit{k}}} \right)$，由$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \mathit{\boldsymbol{D}}_k^{\frac{1}{2}}\mathit{\boldsymbol{E}}_k^{\rm{{\rm T}}}$计算核矩阵$\mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{{\rm T}}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$，即：

$ \mathit{\boldsymbol{K}}\left( {{x_i},{x_j}} \right) = \left\langle {\phi \left( {{x_i}} \right),\phi \left( {{x_j}} \right)} \right\rangle $

(10)

对数据进行中心化处理：

$ {\mathit{\boldsymbol{K}}_{\text{new}}} = \left( {\mathit{\boldsymbol{I}} - \frac{1}{n}\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{E}}^{\rm{{\rm T}}}}} \right)\mathit{\boldsymbol{K}}\left( {\mathit{\boldsymbol{I}} - \frac{1}{n}\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{E}}^{\rm{{\rm T}}}}} \right) $

(11)

根据CPLS，对$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$、Y进行如下分解：

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {\mathit{\boldsymbol{U}}_{\rm{c}}}\mathit{\boldsymbol{R}}_{\rm{c}}^{ - 1} + {\mathit{\boldsymbol{T}}_{\rm{x}}}\mathit{\boldsymbol{P}}_x^{\rm{{\rm T}}} + \mathit{\boldsymbol{ \boldsymbol{\tilde \varPhi} }} $

(12)

$ \mathit{\boldsymbol{Y}} = {\mathit{\boldsymbol{U}}_{\rm{c}}}\mathit{\boldsymbol{Q}}_{\rm{c}}^{\rm{{\rm T}}} + {\mathit{\boldsymbol{T}}_{\rm{y}}}\mathit{\boldsymbol{P}}_{\rm{y}}^{\rm{{\rm T}}} + \mathit{\boldsymbol{\tilde Y}} $

(13)

式中，${\mathit{\boldsymbol{U}}_{\rm{c}}}$为$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$中与质量预测${\mathit{\boldsymbol{\hat Y}}}$相关的协方差部分，${\mathit{\boldsymbol{T}}_{\rm{x}}}$为$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$中与质量预测${\mathit{\boldsymbol{\hat Y}}}$无关的方差部分，${\mathit{\boldsymbol{T}}_{\rm{y}}}$为${\mathit{\boldsymbol{Y}}}$中不能被预测的方差部分，${\mathit{\boldsymbol{P}}_{\rm{x}}}$、${\mathit{\boldsymbol{U}}_{\rm{c}}}$分别为$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}$和${\mathit{\boldsymbol{Y}}}$的负载矩阵，${\mathit{\boldsymbol{R}}_{\rm{c}}}$为满足${\mathit{\boldsymbol{T}}_{\rm{x}}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{R}}_{\rm{c}}}$的权重矩阵。

3.3 多向高斯混合模型-并发核熵潜结构投影(MGMM-CKEPLS)

将三维间歇过程数据$\mathit{\boldsymbol{X}}\left({I \times J \times K} \right)$按批次-变量方向展开得到二维数据$\mathit{\boldsymbol{X}}\left({IK \times J} \right)$，其高斯混合模型^[18]可表示为：

$ p\left( \mathit{\boldsymbol{X}} \right) = \sum\limits_{k = 1}^C {{\omega _k}g\left( {\mathit{\boldsymbol{X}}\left| {{\mu _k},{\Sigma _k}} \right.} \right)} $

(14)

其中C、${{\omega _k}}$、${{\mu _k}}$和${{\Sigma _k}}$分别为高斯成分数、第 k个高斯成分的权重、均值和方差，第 k个多变量高斯密度函数$g\left({\mathit{\boldsymbol{X}}\left| {{\mu _k}, {\Sigma _k}} \right.} \right)$为：

$ g\left( {\mathit{\boldsymbol{X}}\left| {{\mu _k},{\Sigma _k}} \right.} \right) = \frac{1}{{{{\left( {2\pi } \right)}^{J/2}}{{\left| {{\Sigma _k}} \right|}^{1/2}}}}\exp \left[ { - \frac{1}{2}{{\left( {x - {\mu _k}} \right)}^{\rm{{\rm T}}}}\Sigma _k^{ - 1}\left( {x - {\mu _k}} \right)} \right] $

(15)

通过K-means得到初始值${C_0}$、${\omega _0}$、${\mu _0}$和${\Sigma _0}$，由EM迭代过程逐步去掉权重为零的高斯成分，自动确定高斯成分参数${\theta ^{\left({s + 1} \right)}}\left({{C^{\left({s + 1} \right)}}, {\omega ^{\left({s + 1} \right)}}, {\mu ^{\left({s + 1} \right)}}, {\Sigma ^{\left({s + 1} \right)}}} \right)$。

Step 1  (E-Step):

$ {p^{\left( s \right)}}\left( {{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right) = \frac{{\omega _k^{\left( s \right)}g\left( {\mathit{\boldsymbol{X}}\left| {\mu _k^{\left( s \right)},\Sigma _k^{\left( s \right)}} \right.} \right)}}{{\sum\limits_{k = 1}^C {\omega _k^{\left( s \right)}g\left( {\mathit{\boldsymbol{X}}\left| {\mu _k^{\left( s \right)},\Sigma _k^{\left( s \right)}} \right.} \right)} }} $

(16)

式中${p^{\left(s \right)}}\left({{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right)$表示在第s步迭代中，二维样本X的第k个高斯成分的后验概率。

Step 2  (M-Step):

高斯混合模型的迭代公式：

$ \omega _k^{\left( {s + 1} \right)} = \frac{{\sum\limits_{i = 1}^{IK} {{p^{\left( s \right)}}\left( {{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right)} }}{{IK}} $

(17)

$ \mu _k^{\left( {s + 1} \right)} = \frac{{\sum\limits_{i = 1}^{IK} {{p^{\left( s \right)}}\left( {{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right)\mathit{\boldsymbol{X}}} }}{{\sum\limits_{i = 1}^{IK} {{p^{\left( s \right)}}\left( {{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right)} }} $

(18)

$ \Sigma _k^{\left( {s + 1} \right)} = \frac{{\sum\limits_{i = 1}^{IK} {{p^{\left( s \right)}}\left( {{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right)\left( {\mathit{\boldsymbol{X}} - \mu _k^{\left( {s + 1} \right)}} \right){{\left( {\mathit{\boldsymbol{X}} - \mu _k^{\left( {s + 1} \right)}} \right)}^{\rm{{\rm T}}}}} }}{{\sum\limits_{i = 1}^{IK} {{p^{\left( s \right)}}\left( {{C_k}\left| \mathit{\boldsymbol{X}} \right.} \right)} }} $

(19)

其中$\omega _k^{\left({s + 1} \right)}$、$\mu _k^{\left({s + 1} \right)}$和$\Sigma _k^{\left({s + 1} \right)}$分别表示第s+1步迭代中第k个高斯成分的权重、均值和方差。

假设通过GMM聚类后，${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k} = \left[ {{x_1}, \cdots, {x_m}} \right] \in {R^{{{\rm{c}}_k} \times {\rm{m}}}}$和${\mathit{\boldsymbol{Y}}_k} = \left[ {{y_1}, \cdots, {y_p}} \right] \in {R^{{{\rm{c}}_k} \times {\rm{p}}}}$分别为其中一模态的核空间样本数据和质量数据，${{{\rm{c}}_k}}$为第 k个模态的样本数。对${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_k}$和${\mathit{\boldsymbol{Y}}_k}$运用CKEPLS算法，得到${\mathit{\boldsymbol{U}}_{{{\rm{c}}_k}}}$、${\mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}}$和${\mathit{\boldsymbol{T}}_{{{\rm{y}}_k}}}$，然后在每个子空间计算监控统计量。

质量相关的监控指标：

$ \mathit{\boldsymbol{T}}_{{{\rm{c}}_k}}^2 = \left( {n - 1} \right)\mathit{\boldsymbol{U}}_{{{\rm{c}}_k}}^{\rm T}{\mathit{\boldsymbol{U}}_{{{\rm{c}}_k}}} $

(20)

$ {\mathit{\boldsymbol{Q}}_{{{\rm{x}}_k}}} = \mathit{\boldsymbol{ \boldsymbol{\tilde \varPhi} }}_k^{\rm{{\rm T}}}{{\mathit{\boldsymbol{ \boldsymbol{\tilde \varPhi} }}}_k} - 2\mathit{\boldsymbol{ \boldsymbol{\tilde \varPhi} }}_k^{\rm{{\rm T}}} \times {\mathit{\boldsymbol{P}}_{{{\rm{x}}_k}}}{\mathit{\boldsymbol{T}}_{{x_k}}} + \mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}^{\rm{{\rm T}}}\mathit{\boldsymbol{P}}_{{{\rm{x}}_k}}^{\rm{{\rm T}}}{\mathit{\boldsymbol{P}}_{{{\rm{x}}_k}}}{\mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}} $

(21)

质量无关的监控指标：

$ \mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}^2 = \mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}^{\rm{{\rm T}}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{{\rm{x}}_k}}^{ - 1}{\mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}} $

(22)

其中，${\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{{\rm{x}}_k}}} = \frac{1}{{{{\left({n - 1} \right)}^2}}}\mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}^{\rm{{\rm T}}}{\mathit{\boldsymbol{T}}_{{{\rm{x}}_k}}}$。

不能被预测的质量变量监控指标：

$ \mathit{\boldsymbol{T}}_{{{\rm{y}}_k}}^2 = \mathit{\boldsymbol{T}}_{{{\rm{y}}_k}}^{\rm{{\rm T}}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{{\rm{y}}_k}}^{ - 1}{\mathit{\boldsymbol{T}}_{{{\rm{y}}_k}}} $

(23)

$ {\mathit{\boldsymbol{Q}}_{{{\rm{y}}_k}}} = {\left( {{\mathit{\boldsymbol{Y}}_k} - {\mathit{\boldsymbol{Q}}_{{{\rm{c}}_k}}}{\mathit{\boldsymbol{U}}_{{{\rm{c}}_k}}}} \right)^{\rm{{\rm T}}}}\left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{P}}_{{{\rm{y}}_k}}}\mathit{\boldsymbol{P}}_{{{\rm{y}}_k}}^{\rm{{\rm T}}}} \right) \times \left( {{\mathit{\boldsymbol{Y}}_k} - {\mathit{\boldsymbol{Q}}_{{{\rm{c}}_k}}}{\mathit{\boldsymbol{U}}_{{{\rm{c}}_k}}}} \right) $

(24)

其中，${\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{{{\rm{y}}_k}}} = \frac{1}{{{{\left({n - 1} \right)}^2}}}\mathit{\boldsymbol{T}}_{{{\rm{y}}_k}}^{\rm{{\rm T}}}{\mathit{\boldsymbol{T}}_{{{\rm{y}}_k}}}$。

最后，根据GMM得到的各模态权重，对所有模态的监控指标加权求和集成统一的监控量，其监控统计量和控制限如表 1所示，其中，$\chi _{{l_{\rm{c}}}, \alpha }^2$表示自由度为${l_{\rm{c}}}$、置信度为$\alpha $的${\chi ^2}$分布，${g_{\rm{x}}}\chi _{{{\rm{h}}_{\rm{x}}}, \alpha }^2$表示尺度因子为${g_{\rm{x}}}$、自由度为${{{\rm{h}}_{\rm{x}}}}$、置信度为$\alpha $的${\chi ^2}$分布的临界值。由于本文算法考虑了间歇过程模态之间的差异性，更接近于实际工业生产过程，提高了检测精度。

4 基于MGMM-CKEPLS的监控策略 4.1 离线建模

(1) 在无故障工况下获取若干批次过程数据$\mathit{\boldsymbol{X}}\left({I \times J \times K} \right)$，采用3.1节的方法进行数据预处理，最终得到已标准化处理的二维数据$\mathit{\boldsymbol{X}}\left({IK \times J} \right)$；

(2) 根据给定核函数计算核矩阵并通过式(11)中心化；

(3) 对核矩阵进行特征值分解，通过式(9)得到每个特征值所对应的Renyi熵，按照Renyi熵的大小，选取前k个较大的特征值和特征向量，得到数据Φ；

(4) 通过GMM对多模态过程进行聚类，得到 K 个模态的权重、均值和方差；

(5) 根据式(12)、(13)对每一个模态进行分解，并根据式(20)~(24)计算监控统计量，由GMM得到的权值对所有模态的控制限进行加权求和集成统一的监控量。

4.2 在线监测

(1) 将测试数据按3.1节的方法进行数据预处理；

(2) 计算核矩阵、Renyi熵、得到数据Φ；

(3) 根据GMM进行多模态聚类，根据式(12)、(13)对每一个模态进行分解，并根据式(20)~(24)计算监控统计量；

(4) 如果监控统计量没有超过控制限则无故障发生，重复步骤(1)~(3)，或者进行质量预测，否则，产生故障进行报警或生产过程结束。

4.3 质量预测

通过式(25)进行质量预测：

$ \mathit{\boldsymbol{\hat Y}} = {\mathit{\boldsymbol{K}}_{{\rm{new}}}}\mathit{\boldsymbol{U}}{\left( {{\mathit{\boldsymbol{T}}^{\rm{{\rm T}}}}{\mathit{\boldsymbol{K}}_{{\rm{new}}}}\mathit{\boldsymbol{U}}} \right)^{ - 1}}{\mathit{\boldsymbol{T}}^{\rm{{\rm T}}}}\mathit{\boldsymbol{Y}} $

(25)

其中，${\mathit{\boldsymbol{K}}_{{\rm{new}}}}$由式(11)得到。

5 案例研究及分析 5.1 故障监测

青霉素生产制备过程是典型的间歇过程，具有非线性、多模态等特性。本文基于美国Illinois州立理工学院开发的Pensim2.0标准青霉素仿真平台^[19]产生出间歇过程数据，对青霉素发酵过程设置不同初始条件和故障类型进行仿真实验。设定400 h为每批次的反应时间，0.5 h为采样时间，在正常范围设置不同初始条件，不引入故障的情况下共产生35个批次正常工况数据，从产生的18个变量数据中选择其中10个过程变量(通风速率、搅拌功率、补料温度、溶解氧浓度、反应器体积、排气二氧化碳浓度、pH值、发酵罐温度、产生热、冷水流加速率)和2个质量变量(菌体浓度、产物浓度)作为监控变量，构成三维数据X(35×10×800)。

Pensim2.0提供了三种故障类型(通风速率、搅拌功率、底物流加速率)。本文引入通风速率和搅拌功率作为故障，详细的故障参数设置如表 2所示。仿真采用高斯核函数，核参数选取32，多模态聚类结果如图 2所示，得到3个模态批次。分别用MKPLS和本文提出的MGMM-CKEPLS算法进行监控，本文给出故障1 (在时间200~300 h，即采样点400~600，加入10%的阶跃信号作为故障信号)的监控图，如图 3~8所示，其余监控结果见表 3。

图 3、图 4为基于传统方法的MKPLS监控图。从图 3、图 4可以看出，MKPLS在间歇过程监控中有较大的误报、漏报产生。虽然其 Q统计量的漏报率为1%，但误报率达到25%，其T²统计量的误报和漏报率分别为8.67%和15%，明显高于本文算法的监控统计量。由表 3对不同故障的监测结果可知，针对间歇过程的多模态和非线性，MGMM-CKEPLS优于MKPLS方法。

图 5~8为所提算法的监控图。图 5为质量相关变量的监控图，图 6为质量不相关的过程变量监控图，从图 5、图 6可以看出，针对质量相关的变量，发生故障后会向正常工况靠近，而质量不相关的过程变量可以全部检测出故障，误报率低；图 7为过程变量残差的监控图，误报、漏报率分别为0.83%和0.5%，表明本文提出的算法对过程空间变量进行分解的有效性；图 8为不能被预测的质量变量的监控图，明显看出发生故障后又很快回到正常工况。基于以上分析，针对间歇过程的非线性和多模态问题，本文提出的算法能快速高效地检测出故障，误报、漏报率低，提高了检测精度。

5.2 质量预测

本文选取青霉素发酵过程产生的两个质量变量(菌体浓度、产物浓度)进行产品质量预测，与传统方法MKPLS进行比较，质量预测结果如图 9~12所示，根据式(26)所示的均方误差(RMSE)计算预测精度，结果如表 4所示。

$ {\rm{RMSE}} = \sqrt {\frac{1}{\mathit{K}}\sum\limits_{i = 1}^\mathit{K} {{{\left( {{{\mathit{\tilde y}}_\mathit{i}} - {\mathit{y}_\mathit{i}}} \right)}^2}} } $

(26)

其中，K为采样时间，${{{\mathit{\tilde y}}_\mathit{i}}}$和${{\mathit{y}_\mathit{i}}}$分别为第 i个时刻的预测值和实际值。

由式(26)可知，RMSE越大，预测精度越低。图 9~10分别为MKPLS对菌体浓度和产物浓度的质量预测，可以明显看出预测值跟实际值有较大偏差，其RMSE分别为0.4836和0.3677，均方误差较大，预测精度低；图 11~12为本文算法的质量预测，其RMSE分别为0.0634和0.0298，均方误差小，预测精度高。因此，本文算法相比MKPLS算法，提高了产品质量预测精度。

6 结论

强非线性、多模态是间歇过程的固有特性。本文提出的基于MGMM-CKEPLS过程监测及质量预测方法，通过核熵投影解决了非线性数据的同时降低了主元数；利用GMM对复杂数据处理的优点对多模态过程进行聚类，突出了各模态的差异性；然后在不同模态构建CPLS模型进行故障监测和质量预测。通过仿真实验表明本文方法比MKPLS有更好的监控效果，降低了误报率、漏报率，提高了故障监控的快速性和精确性，对产品质量有更高的预测精度。