1 前言

在过程工业中，如化工、冶金、生物、制药、污水处理等领域，很多关键的质量参数无法实时在线测量，严重制约了过程的监测、控制及优化水平。软测量技术为此类难测参数的在线估计提供了有效途径^[1-2]。随着工业生产过程的现代化，大量隐含有关生产状况的操作数据得以保存下来。因此，基于数据驱动的软测量建模方法日益受到青睐^[3-6]。此类方法的典型代表有：偏最小二乘^[7]、人工神经网络^[8]、支持向量回归^[9]、高斯过程回归^[10]等。

近年来，集成学习作为一种有效的机器学习与数据挖掘工具，在软测量领域得到了广泛应用^[11-15]。一般而言，集成学习软测量建模包含2个关键步骤：基模型的生成和融合。大量研究表明，集成学习成败的关键在于能否建立一组准确性高、多样性强的基模型^[16-17]。其中，基模型的多样性生成机制尤为重要。然而，常规集成学习软测量建模方法主要通过样本扰动方式生成基模型，如聚类^[18]、移动窗口^[19]、自助采样^[20]、序列采样^[21]等策略，忽略了输入变量选择的多样性。实际上，输入变量也是影响软测量模型性能的关键因素^[22-25]。然而，在实际应用中，难以获得单一的最优变量组合，而获得一系列次优的输入变量组合则相对容易得多。因此，本文重点探索基于输入特征扰动的集成学习软测量建模方法。

此外，构建高性能集成学习软测量模型的另一个关键步骤在于基模型的融合。根据ZHOU等^[26]的研究，集成学习中会出现“Many Could Be Better Than All”的现象，即在基模型生成之后，通过移除一些基模型来获得较小的集成，不仅有助于减小模型的存储开销和计算负担，而且有可能提升模型的泛化性能。该操作称为集成修剪，是构建高性能选择性集成学习模型的关键所在。究其根本，集成修剪的目的在于确保基模型更好地满足准确性和多样性的要求，本质上这是一个优化问题。但是，现有的基于优化的修剪方法仅考虑单一的优化目标，无法有效确保准确性和多样性目标的平衡。

针对上述问题，提出一种基于进化多目标优化的选择性集成高斯过程回归(evolutionary multi-objective optimization based selective ensemble Gaussian process regression，EMO-SEGPR)软测量建模方法。该方法首先引入输入特征扰动的多样性产生机制，结合bootstrapping随机重采样和偏互信息(partial mutual information，PMI)相关分析构建一组多样性的输入特征子集，并据此建立一组多样性的高斯过程回归(Gaussian process regression，GPR)基模型。然后，基于进化多目标优化(evolutionary multi-objective optimization，EMO)算法实现集成修剪，从而获得一组规模较小而且同时满足多样性和准确性要求的GPR基模型。最后，采用Stacking集成策略实现多样性基模型的融合。在应用案例部分，通过对青霉素发酵过程和Tennessee Eastman (TE)化工过程的建模研究，验证了EMO-SEGPR软测量方法的有效性和优越性。

2 算法介绍 2.1 高斯过程回归

高斯过程(Gaussian process, GP)是任意有限个随机变量均具有联合高斯分布的集合。对于输入输出样本集$\boldsymbol{D}=\{\boldsymbol{X}, \boldsymbol{y}\}=\left\{\boldsymbol{x}_{i}, y_{i}\right\}_{i=1}^{n}$，考虑如下回归模型

$y = f\left( \boldsymbol{x} \right) + \varepsilon , \varepsilon \sim N\left( {0, \sigma _\text{n}^2} \right)$

(1)

其中，x为输入向量，n为输入样本点的总数，f(·)为未知函数，y是受到噪声干扰的观测值，$\varepsilon$是均值为0、方差为$\sigma _\text{n}^2$的高斯噪声。在函数f空间上，一个高斯过程完全由均值函数m(x)和协方差函数C(x, x')确定，即

$m\left( \boldsymbol{x} \right) = \text{E}\left[ {f\left( \boldsymbol{x} \right)} \right]$

(2)

$C\left( {\boldsymbol{x}, \boldsymbol{x}'} \right) = \text{E}\left[ {\left( {f\left( \boldsymbol{x} \right) - m\left( \boldsymbol{x} \right)} \right)\left( {f\left( {\boldsymbol{x}'} \right) - m\left( {\boldsymbol{x}'} \right)} \right)} \right]$

(3)

如果数据已经进行了适当的预处理，可以假设训练样本集产生于一个零均值高斯过程，即输出观测值服从

$\boldsymbol{y} \sim N\left( {0, {\boldsymbol{C}}} \right)$

(4)

其中C为一个n×n对称正定的协方差矩阵，$\boldsymbol{x}'$为$\boldsymbol{x}$的转置计算。

对于一个新的测试输入${\boldsymbol{x}_ * }$，训练输出y和测试输出${y_ * }$的联合先验分布为

$\left[ {\begin{array}{*{20}{c}} {\boldsymbol{y}} \\ {{y_ * }} \end{array}} \right] = N\left( {0, \left[ {\begin{array}{*{20}{c}} {\boldsymbol{C}}&{{\boldsymbol{k}_ * }} \\ {\boldsymbol{k}_ * ^\text{T}}&{C\left( {{\boldsymbol{x}_ * }, {\boldsymbol{x}_ * }} \right)} \end{array}} \right]} \right)$

(5)

其中，${\boldsymbol{k}_ * } = {\left[ {C\left( {{{\boldsymbol{x}}_ * }, {{\boldsymbol{x}}_1}} \right), \cdots , C\left( {{{\boldsymbol{x}}_ * }, {{\boldsymbol{x}}_n}} \right)} \right]^\text{T}}$为测试点${{\boldsymbol{x}}_ * }$与训练集输入之间的n×1阶协方差矩阵；$C\left( {{{\boldsymbol{x}}_ * }, {{\boldsymbol{x}}_ * }} \right)$为测试点${{\boldsymbol{x}}_ * }$自身的协方差。给定了协方差函数，测试输出${y_ * }$的后验分布为

${y_ * }\left| {{\boldsymbol{X}}, {\boldsymbol{y}}, {{\boldsymbol{x}}_ * }} \right. \sim N\left( {{{\hat y}_ * }, \sigma _ * ^2} \right)$

(6)

这里，${\hat y_ * }$和$\sigma _ * ^2$按下式计算

${\hat y_ * } = \text{E}\left( {{y_ * }} \right) = {\boldsymbol{k}}_ * ^\text{T}{{\boldsymbol{C}}^{ - 1}}{\boldsymbol{y}}$

(7)

$\sigma _ * ^2 = {\rm{Var}}\left( {{y_ * }} \right) = C\left( {{{\boldsymbol{x}}_ * }, {{\boldsymbol{x}}_1}} \right) - {\boldsymbol{k}}_ * ^\text{T}{{\boldsymbol{C}}^{ - 1}}{{\boldsymbol{k}}_ * }$

(8)

其中，$E\left( \cdot \right)$和${\rm{Var}}\left( \cdot \right)$分别表示均值和方差算子。${\hat y_ * }$和$\sigma _ * ^2$分别表示GPR模型测试样本点${{\boldsymbol{x}}_ * }$的预测值和预测方差。

GPR模型的建立首先要确定协方差函数，所选的协方差函数要能在任何数据集上产生非负定协方差矩阵，同时希望相邻的输入产生相邻的输出。因此，选择带噪声项的Mateŕn协方差

$C\left( {{{\boldsymbol{x}}_i}, {{\boldsymbol{x}}_j}} \right) = \sigma _\text{f}^2\left( {1 + \frac{{\sqrt 3 \parallel {{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}\parallel }}{l}} \right) = \exp \left( { - \frac{{\sqrt 3 \parallel {{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}\parallel }}{l}} \right) + \sigma _n^2{\delta _{ij}}$

(9)

其中，$\sigma _\text{f}^2$为输出尺度，l为输入尺度，$\sigma _n^2$为噪声方差；i = j时，${\delta _{ij}}{\rm{ = }}\sigma _\text{f}^2$，否则$\delta$_ij = 0。超参数集合$\mathit{\Theta} = \left\{ {\sigma _\text{f}^2, 1, \sigma _n^2} \right\}$通过极大似然法求得。训练样本上的对数似然函数为

$\log p\left( {{\boldsymbol{y}}\left| {\boldsymbol{X}} \right.} \right) = - \frac{1}{2}{{\boldsymbol{y}}^\text{T}}{{\boldsymbol{C}}^{ - 1}}{\boldsymbol{y}} - \frac{1}{2}\log \left| {\boldsymbol{C}} \right| - \frac{n}{2}\log \left( {2\pi } \right)$

(10)

对超参数求偏导：

$\frac{{\partial \left[ {\log p\left( {{\boldsymbol{y}}\left| {\boldsymbol{X}} \right.} \right)} \right]}}{{\partial \theta }} = - \frac{1}{2}{\rm{tr}}\left( {{{\boldsymbol{C}}^{ - 1}}\frac{{\partial {\boldsymbol{C}}}}{{\partial \theta }}} \right) + \frac{1}{2}{{\boldsymbol{y}}^\text{T}}{{\boldsymbol{C}}^{ - 1}}\frac{{\partial \boldsymbol{C}}}{{\partial \theta }}{{\boldsymbol{C}}^{ - 1}}{\boldsymbol{y}}$

(11)

其中，$\frac{{\partial {\boldsymbol{C}}}}{{\partial \theta }}$可以从协方差函数求得，${\rm{tr}}\left( \cdot \right)$表示计算矩阵的迹。

GPR建模技术由于存在非线性处理能力强、能提供预测不确定信息、参数设定简单等独特优势，在本文中将用于构建集成学习基模型。

2.2 偏互信息

互信息(mutual information，MI)是信息论里的一种非参数、非线性相关性评价指标。MI能评估输入变量与输出变量之间的相关程度，但是忽略了输入变量之间的独立性。而PMI^[27]不仅可以评价候选变量与输出变量之间的相关性，同时还考虑了输入变量之间的冗余度问题。PMI的计算式如下：

${\rm{PMI}}\left( {X, Y\left| Z \right.} \right) = \iint {{f_{X', Y'}}\left( {{{x'}_i}, {{y'}_i}} \right)}\ln \left[ {\frac{{{f_{X', Y'}}\left( {{{x'}_i}, {{y'}_i}} \right)}}{{{f_{X'}}\left( {{{x'}_i}} \right){f_{Y'}}\left( {{{y'}_i}} \right)}}} \right]\text{d}x'\text{d}y'$

(12)

其中，

$x' = x - \text{E}\left[ {x\left| Z \right.} \right], y' = y - \text{E}\left[ {y\left| Z \right.} \right]$

(13)

式中，x'，y'表示样本数据集的冗余成分，${f_{X'}}\left( {X'} \right)$和${f_{Y'}}\left( {Y'} \right)$是边缘概率密度，f_X'Y'(x', y')是联合概率密度。给定已选输入变量集Z，PMI满足：(ⅰ)对称性，即${\rm{PMI}}\left( {X, Y\left| Z \right.} \right) = {\rm{PMI}}\left( {Y, X\left| Z \right.} \right)$；(ⅱ)非负性，${\rm{PMI}}\left( {X, Y\left| Z \right.} \right) \geqslant 0$，在Z条件下X与Y相互独立时等号成立。$Z = \mathit{\Phi} $时，${\rm{PMI}}\left( {X, Y\left| Z \right.} \right) = {\rm{MI}}\left( {X, Y} \right)$。

基于离散样本的PMI值可估计如下：

${\rm{PMI}}\left( {X, Y\left| Z \right.} \right) = \frac{1}{n}\sum\nolimits_{i = 1}^n {\ln \left[ {\frac{{{f_{X', Y'}}\left( {{{x'}_i}, {{y'}_i}} \right)}}{{{f_{X'}}\left( {{{x'}_i}} \right){f_{Y'}}\left( {{{y'}_i}} \right)}}} \right]} $

(14)

其中，${x'_i}, {y'_i}$表示数据集中第$i$个样本的剩余信息。

3 EMO-SEGPR软测量建模 3.1 多样性GPR基模型构建

在集成学习中，增强基模型的多样性能有效提升集成模型的泛化性能。本文通过挖掘输入特征扰动机制，从而构建一组多样性的基模型。常见的特征扰动方法如随机子空间法，更适用于特征较多的样本数据，若特征较少或冗余特征很少，反而不利于模型性能的提升^[28]。而且，输入特征的随机选择将难以保证基模型的预测性能。一种更为有效的特征选择方法是，采用某种相关性准则评估待选变量与输出变量之间的相关性。但是，其选择结果易受建模样本选择的影响。因此，提出一种样本采样与相关分析相结合的多样性子空间构建方法。其基本思想是：首先通过bootstrapping重采样获得一组多样性的建模样本子集$\left\{ {{{\boldsymbol{X}}_m}, {{\boldsymbol{y}}_m}} \right\}_{m = 1}^M$，然后对每个子集进行PMI相关分析，从而获得一组多样性的输入子空间$\left\{ {{{\boldsymbol{S}}_m}} \right\}_{m = 1}^M$。

为了选择与输出变量相关性较强的变量，同时解决输入变量之间的冗余问题，选用PMI相关性准则。首先，采用k最近邻(k-nearest neighbor，kNN)法^[27]对PMI值进行估计，其中使用k折交叉验证和置换检验方法^[29]确定最优最近邻数目。然后，确定PMI阈值(决定何时停止输入变量选择)，本文使用一种统计置信限来判断输入与输出变量是否相关。对于一个候选变量，重复p次bootstrapping算法产生多个随机化变量，然后从中确定第α百分位，以此作为变量相关性的重要性判别阈值。若原始变量的PMI值大于阈值，则认为该变量与输出变量相关。在本研究中，p = 100，α = 95%。

对于一个bootstrapping重采样子集，PMI变量选择步骤如下：

(1) 初始化输入变量集$S = \left\{ {{X_i}} \right\}_{i = 1}^d$，已选变量集$S' = \mathit{\Phi }$；Y为输出变量；

(2) 确定kNN估计器的最佳最近邻数目，并计算每个输入变量与输出变量的${\rm{PMI}}\left( {{X_i}, Y\left| {S'} \right.} \right)$；

(3) 选出(2)中PMI得分最高的变量$X'$，如果该PMI值高于第95百分位随机化PMI阈值，则将其加入$S'$，即$S \leftarrow S - \left\{ {{X_i}} \right\}, S' \leftarrow S' + \left\{ {{X_i}} \right\}$，否则终止变量选择；

(4) 重复(2)和(3)直至选出所有重要变量。

通过上述方法获得一组多样性的输入子空间，由对应的变量索引提取原始训练样本，从而构建一组多样性的GPR基模型$\left\{ {{\rm{GPR}}} \right\}_{m = 1}^M$。GPR模型是一种概率性的非参数模型，其优势在于模型输出具有概率意义，可以同时提供预测均值和预测方差。其中，预测方差可以用来评估预测值的可信度。${\left\{ {{\rm{GPR}}} \right\}_m}$预测模型可以描述如下：

$\left\{ \begin{array}{l} {{\hat y}_{m, {\rm{new}}}} = \text{E}\left( {{y_{m, {\rm{new}}}}} \right) = {\boldsymbol{k}}_{m, {\rm{new}}}^\text{T}{\boldsymbol{C}}_m^{ - 1}{\boldsymbol{y}_m} \\ \sigma _{m, {\rm{new}}}^2 = {\rm{Var}}\left( {{y_{m, {\rm{new}}}}} \right) = C\left( {{\boldsymbol{x}_{m, {\rm{new}}}}, {\boldsymbol{x}_{m, {\rm{new}}}}} \right) - {\boldsymbol{k}}_{m, {\rm{new}}}^\text{T}{\boldsymbol{C}}_m^{ - 1}{{\boldsymbol{k}}_{m, {\rm{new}}}} \end{array} \right.$

(15)

式中，${{\boldsymbol{k}}_{m, {\rm{new}}}} = {\left[ {C\left( {{\boldsymbol{x}_{{\rm{new}}}}, {\boldsymbol{x}_{m, {\rm{1}}}}} \right), \cdots , C\left( {{\boldsymbol{x}_{{\rm{new}}}}, {\boldsymbol{x}_{m, n}}} \right)} \right]^\text{T}}$，${\hat y_{m, {\rm{new}}}}$和$\sigma _{m, {\rm{new}}}^2$分别为基模型${\left\{ {{\rm{GPR}}} \right\}_m}$的预测均值和方差。

3.2 基于进化多目标优化的集成修剪

在GPR基模型的构建过程中，难免存在部分GPR基模型相关性强、预测性能不佳的现象。若将所有基模型用于集成建模，势必造成集成模型复杂度的提升，甚至恶化预测性能。解决这一问题的有效方法是移除部分基模型，即集成修剪^[30]。

对初始GPR基模型执行集成修剪的核心问题在于，如何挑选一组较小规模的GPR基模型，而且尽可能保证基模型的预测性和多样性。然而，这两个指标往往存在一定程度的冲突，是典型的多目标优化问题。因此，本文将基模型的集成修剪问题转化为如下多目标优化问题：

$\max \left[ {{f_1}, {f_2}} \right]$

(16)

f₁，f₂为目标函数，分别衡量基模型的预测精度和多样性。

为求解上述多目标优化问题，首先需要确定决策向量、约束条件、目标函数的具体形式。假设生成了M个GPR基模型，对每一个基模型进行二进制编码，其中每一位编码表示是否选中相应的GPR基模型，1表示选中该模型，0表示未选中，然后将这一组二进制变量作为决策向量。为了有效控制选择到的模型个数，将模型选择数目M_select作为约束条件。

对于基模型的预测精度指标，先估计单个模型的预测精度，再将所有基模型的估计值平均作为整体预测精度。为此，将基模型均方根误差(root-mean-square error, RMSE)作为预测精度指标，RMSE越小，则基模型预测精度越高。给定验证样本${{\boldsymbol{D}}_{{\rm{val}}}} = \left\{ {{{\boldsymbol{X}}_{{\rm{val}}}}, {{\boldsymbol{y}}_{{\rm{val}}}}} \right\}$，精度目标函数定义为：

${\rm{RMS}}{{\rm{E}}_{{\rm{avg, val}}}} = \frac{{\sum\nolimits_{i = 1}^M {{\rm{RMSE}}_{{\rm{val}}}^i} }}{M}$

(17)

其中，${\rm{RMSE}}_{{\rm{val}}}^i$为基于验证样本的单个${\rm{GP}}{{\rm{R}}_m}$模型的预测均方根误差。

对于基模型的多样性指标，目前仍没有统一的定义。本文提出使用相关系数作为多样性衡量指标，其基本思想是：两个基模型之间的差异越大，预测误差系列之间的相关系数指标就越小。为衡量集成多样性，先估计任意两个基模型之间的差异性，即相关性系数：

$r\left( {{\boldsymbol{e}_i}, {\boldsymbol{e}_j}} \right) = \frac{{{\rm{Cov}}\left( {{\boldsymbol{e}_i}, {\boldsymbol{e}_j}} \right)}}{{\sqrt {{\rm{Var}}\left( {{\boldsymbol{e}_i}} \right){\rm{Var}}\left( {{\boldsymbol{e}_j}} \right)} }}$

(18)

其中，${\boldsymbol{e}_i}, {\boldsymbol{e}_j}$分别表示任意两个GPR基模型的预测误差，${\rm{Cov}}\left( \cdot \right)$用于计算任意两个误差之间的协方差，${\rm{Var}}\left( \cdot \right)$表示方差算子。最后将(18)式得出的所有相关系数值平均作为集成多样性指标：

${r_{{\rm{avg, val}}}} = \frac{{\sum\nolimits_{j = i + 1}^{{M_\text{select}}} {\sum\nolimits_{i = 1}^{{M_\text{select}} - 1} {r\left( {{\boldsymbol{e}_i}, {\boldsymbol{e}_j}} \right)} } }}{{{{\left( {M_\text{select}^2 - {M_\text{select}}} \right)} \mathord{\left/ {\vphantom {{\left( {M_\text{select}^2 - {M_\text{select}}} \right)} 2}} \right. } 2}}}$

(19)

综上分析可知，最大化基模型的预测精度和多样性相当于最小化式(17)和式(19)两个衡量指标。因此，式(16)中的最大化多目标优化问题转化为如下最小化优化问题：

$\min \left[ {{\rm{RMS}}{{\rm{E}}_{{\rm{avg, val}}}}, {r_{{\rm{avg, val}}}}} \right]$

(20)

接下来对上述多目标优化问题进行求解。传统最优化方法对目标函数、约束函数要求较高，而且最优性达到的条件太苛刻。相比而言，进化优化通过模拟生物进化过程与机制求解，其优势在于通用性较强、对优化问题没有过多的特殊要求。因此，进化多目标优化算法在实际应用中日益受到青睐^[31-33]。本文选用的是带精英策略的非支配排序的遗传算法(non-dominated sorting genetic algorithmⅡ，NSGA-Ⅱ)，由2002年DEB等^[34]对其算法NSGA的改进，是目前为止最优秀的进化多目标优化算法之一。

采用NSGA-Ⅱ优化实现集成修剪的主要步骤如下：

(1) 对基模型的选择问题进二进制编码，随机初始化每个个体的染色体，生成大小为P的初始种群${P_t} = 0$。

(2) 对种群P_t中的每一个个体进行评估。从染色体上解码决策向量，确定被选中的GPR基模型，然后基于验证样本集评估每个被选中基模型的预测性能，并计算两个优化目标。

(3) 重复以下步骤，直到满足终止条件。

3a) 对种群进行非支配排序，并计算所有非支配解的拥挤距离。

3b) 利用二进制锦标赛选择、交叉和变异，生成一个与P_t相同大小的子代Q_t。

3c) 重复步骤(2)对Q_t进行评估。

3d) 结合${R_t} = {P_t} \cup {Q_t}$，保留精英父代，根据非支配排序方法进行排序。

3e) 从R_t中获取P个解来创建新的一代P_t+1，增加进化代数t + 1 = t。

(4) 使用非支配排序找出在上一代合并后的个体中的最优解。

通过设置合适的种群数和迭代数进行优化，得到Pareto最优解集，其中任意一个Pareto解对应一个基模型选择的二进制变量组合。然后对二进制染色体串解码，由此选出参与集成的GPR基模型，从而实现基模型的集成修剪。

3.3 多样性GPR基模型集成

通过3.2节的集成修剪操作，最终选出M_select个GPR基模型参与集成。给定测试样本X_new，相应可获得M_select个预测值。最后融合局部预测值得到最终估计值。要达到满意的预测效果，需选择合理的集成策略。由于基模型的预测性能往往存在差异，因此加权集成策略是一种合理选择。当训练数据很多时，一种强大的集成策略是学习法，即通过另一个学习器进行结合。Stacking^[35]是学习法的典型代表，它先从初级数据集训练初级学习器，然后生成一个新数据集用于训练次级学习器。在新的训练数据中，初级学习器的预测输出被用作次级学习器的输入。

本文在Stacking框架下采用偏最小二乘法(partial least squares，PLS)对GPR基模型进行集成，最佳主成分个数由交叉验证确定。假设$\left\{ {{{\hat y}_m}} \right\}_{m = 1}^{{M_\text{select}}}$、$\left\{ {\sigma _m^2} \right\}_{m = 1}^{{M_\text{select}}}$为基模型的预测输出和方差，则基于PLS的Stacking集成输出可以表示为：

$\hat y{\rm{ = }}{w_0}{\rm{ + }}{w_1}{\hat y_1} + {w_2}{\hat y_2} + \cdots + {w_{{M_\text{select}}}}{\hat y_{{M_\text{select}}}}$

(21)

根据不确定度合成原理，集成预测方差${\sigma ^2}$可计算为^[36]：

${\sigma ^2}{\rm{ = }}{\sum\nolimits_{i = 1}^{{M_\text{select}}} {\left( {\frac{{\partial y}}{{\partial {y_i}}}} \right)} ^2}\sigma _i^2 + 2\sum\nolimits_{1 \leqslant i \leqslant j}^{{M_\text{select}}} {\frac{{\partial y}}{{\partial {y_i}}}\frac{{\partial y}}{{\partial {y_j}}}} {\rho _{ij}}{\sigma _i}{\sigma _j}$

(22)

其中，${\sigma _i}$和${\sigma _j}$为任意两个基模型的预测输出不确定度；${\rho _{ij}}$为不确定度变量${\sigma _i}$和${\sigma _j}$的之间的相关系数。

根据式(21)可知，${w_i} = {{\partial y} \mathord{\left/ {\vphantom {{\partial y} {\partial {y_i}}}} \right. } {\partial {y_i}}}$，因此式(22)所示的集成预测方差可重写为：

${\sigma ^2} = \sum\nolimits_{i = 1}^{{M_\text{select}}} {w_i^2} \sigma _i^2 + 2\sum\nolimits_{1 \leqslant i \leqslant j}^{{M_\text{select}}} {{w_i}{w_j}} {\rho _{ij}}{\sigma _i}{\sigma _j}$

(23)

由式(23)可以看出，集成预测方差取决于基模型自身的预测方差水平和预测方差之间的相关程度。如果基模型完全独立，则${\rho _{ij}} = 0$，集成预测方差退化为${\sigma ^2} = \sum\nolimits_{i = 1}^{{M_\text{select}}} {w_i^2} \sigma _i^2$。但在实际应用中，基模型之间难免存在某种程度的相关性，即${\rho _{ij}} \ne 0$。由于真实的${\rho _{ij}}$值无法获得，本文按以下方法对其进行估计。给定验证样本集${{\boldsymbol{D}}_{{\rm{val}}}} = \left\{ {{{\boldsymbol{X}}_{{\rm{val}}}}, {{\boldsymbol{y}}_{{\rm{val}}}}} \right\}$，可获得任意两个GPR基模型在${{\boldsymbol{D}}_{{\rm{val}}}}$上的预测均值向量和方差向量$\left\{ {{{{\boldsymbol{\hat y}}}_{i, {\rm{val}}}}, {{\boldsymbol{\sigma }}_{i, {\rm{val}}}};{{{\boldsymbol{\hat y}}}_{j, {\rm{val}}}}, {{\boldsymbol{\sigma }}_{j, {\rm{val}}}}} \right\}$。此时，${\rho _{ij}}$可估计为：

${\rho _{ij}} = \frac{{{\rm{Cov}}\left( {{{\boldsymbol{\sigma }}_{i, {\rm{val}}}}, {{\boldsymbol{\sigma }}_{j, {\rm{val}}}}} \right)}}{{\sqrt {{\rm{Var}}\left( {{{\boldsymbol{\sigma }}_{i, {\rm{val}}}}} \right){\rm{Var}}\left( {{{\boldsymbol{\sigma }}_{j, {\rm{val}}}}} \right)} }}$

(24)

除了“学习法”，GPR基模型还可以通过其他方式进行融合。例如，一种应用较为广泛的方法是通过GPR的预测方差信息实现融合^[37-38]。

3.4 实施原理

EMO-SEGPR软测量方法的基本原理如图 1所示。具体实施步骤如下：

(1) 采集输入输出样本，将样本分为训练集、验证集和测试集；

(2) 通过bootstrapping随机重采样，生成多个训练样本子集；

(3) 采用PMI准则对每个训练样本子集进行相关分析，从而构建多样性输入子空间，然后提取对应的原始训练样本，建立GPR基模型；

(4) 通过进化多目标优化算法对(3)得到的GPR基模型实施集成修剪得到Pareto最优解集；

(5) 选取一个Pareto最优解进行解码，确定最终选中的GPR基模型，并建立PLS集成模型；

(6) 对新的测试样本，首先给出GPR基模型的预测输出和预测方差，然后采用PLS stacking集成融合，最终得到测试样本的集成预测输出和预测方差。

4 应用案例研究

为了验证EMO-SEGPR软测量建模方法的有效性，本文以青霉素发酵过程和TE化工过程作为应用案例。实验中比较了如下方法：(1) GPR：使用单一的输入变量集建立全局模型；(2) EGPR：仅使用多样性的输入子空间进行集成建模；(3) EMO-SEGPR：本文提出的方法，使用多样性的输入子空间进行集成建模，并考虑集成修剪。此外，EGPR、EMO-SEGPR方法对比了3种集成策略：简单平均法(simple averaging rule，SA)、预测方差集成(prediction variance based ensemble, PVE)^[20]、PLS stacking集成。

模型预测性能评价采用均方根误差RMSE和决定系数${{\rm{R}}^2}$：

${\rm{RMSE}} = \sqrt {\frac{1}{{{N_{{\rm{test}}}}}}\sum\nolimits_{i = 1}^{{N_{{\rm{test}}}}} {{{\left( {{{\hat y}_i} - {y_i}} \right)}^2}} } $

(25)

${{\rm{R}}^2} = 1 - \frac{{\sum\nolimits_{i = 1}^{{N_{{\rm{test}}}}} {{{\left( {{{\hat y}_i} - {y_i}} \right)}^2}} }}{{\sum\nolimits_{i = 1}^{{N_{{\rm{test}}}}} {{{\left( {{{\hat y}_i} - {{\bar y}_i}} \right)}^2}} }}$

(26)

其中，${N_{{\rm{test}}}}$为测试样本的数目，${\hat y_i}$, ${y_i}$和${\bar y_i}$分别为估计值、实际值和实际输出的均值。

4.1 青霉素发酵过程

青霉素发酵过程被广泛用于间歇过程的建模、监控和控制的研究。图 2给出了青霉素发酵过程的工艺流程图。本文以青霉素浓度作为主导变量，以表 1所列过程变量作为辅助变量建立软测量模型。建模数据来源于PenSim 2.0^[39]平台，发酵周期为400 h，采样时间间隔为2 h，共收集20批青霉素发酵过程数据。将发酵数据分为训练集(8批，50%)、验证集(4批，25%)、测试集(4批，25%)。其中，训练集用于模型训练，验证集用于多目标优化适应度函数的计算，测试集用于软测量模型的性能评估。

图 3显示了EMO-SEGPR建模中NSGA-Ⅱ优化产生的Pareto前沿图。其中，进化的种群数、迭代数、选择的基模型数分别取500、100和20。实验中，仅选择其中一个Pareto解进行解码，并用于EMO-SEGPR建模。本案例研究中，初始构建了66个GPR基模型，经过NSGA-Ⅱ优化最终选择其中20个基模型参与集成建模。图 4显示了基模型的二进制选择结果。

表 2比较了不同软测量建模方法在青霉素发酵过程中的预测误差。可以看出，GPR的预测精度远低于EGPR和EMO-SEGPR方法，这是因为GPR并未考虑输入变量选择的多样性。相比而言，EGPR和EMO-SEGPR通过结合输入特征扰动和集成学习，获得了显著的性能提升。此外，实验比较了3种集成策略，可以看出，EMO-SEGPR方法有效剔除了冗余子模型，从而大幅度降低了集成建模的复杂度。而且，EMO-SEGPR在3种集成策略下均表现出最佳预测性能，其中PLS stacking集成效果最佳。因此，本文提出的EMO-SEGPR方法表现出了较好的预测性能。此外，图 5显示了基于EMO-SEGPR (PLS stacking)方法的青霉素浓度预测趋势曲线。由图可见，EMO-SEGPR模型预测值与实际值高度吻合，由此进一步说明EMO-SEGPR软测量方法的良好性能。

4.2 TE化工过程

TE化工过程^[40]是基于实际工业过程的模拟仿真，其工艺流程如图 6所示。TE生产过程主要有A、C、D和E四种气态物料参与反应，生产出两种产品G、H，并伴有一种副产品F。TE过程主要包括41个测量变量和12个操纵变量。在本实验中，选用如表 3中所示的22个连续测量变量和12个操纵变量作为原始输入，Stream 9中的E成分浓度作为软测量模型的输出。输入和输出数据在5个不同操作模式下以12 min的采样间隔采集获得，并将其分为训练集(50%)、验证集(25%)、测试集(25%)。

图 7给出了由NSGA-Ⅱ优化的Pareto前沿。其中，进化的种群数、迭代数、选择的GPR基模型数分别取300、100和30。对TE过程，通过输入特征扰动构建了85个多样性GPR基模型，然后进行优化选出30个最佳基模型，集成修剪结果如图 8所示。

表 4给出了不同软测量建模方法的预测结果。同时，图 9给出了基于EMO-SEGPR (PLS stacking)的E成分预测趋势曲线。由表 4可知，基于输入特征扰动的集成学习软测量模型EGPR和EMO-SEGPR明显优于常规的全局模型GPR。相比于EGPR，集成修剪策略的引入使得EMO-SEGPR实现了模型复杂度的大幅下降，同时维持了较高的预测性能。对比不同的集成策略可看出，基于PLS stacking的集成效果显著优于其他两种集成方法。总体上，所提出的EMO-SEGPR软测量建模方法表现出较高的预测性能。从图 9中也可以看出，基于EMO-SEGPR方法的E成分预测输出与实际输出高度吻合，说明该方法可以获得满意的预测性能。

5 结论

本文提出了一种基于进化多目标优化的选择性集成学习软测量建模方法EMO-SEGPR。该方法通过bootstrapping随机重采样和PMI准则产生多样性子空间，然后采用GPR技术构建局部模型，有效增强了基模型的多样性。随后将多样性基模型的集成修剪转化为多目标优化问题，并通过EMO算法实现基模型的精简，最后将得到的基模型进行融合。通过青霉素发酵过程和TE化工过程验证了该方法的有效性。EMO-SEGPR方法充分考虑了输入变量选择的多样性，能有效保证基模型的多样性和预测精度。而且，集成修剪的引入有效克服了传统集成学习融合所有局部模型的缺陷，既显著降低了集成建模的复杂度，又提升了集成预测性能。此外，相比于传统的集成学习方法，EMO-SEGPR方法不仅能有效提高模型的估计精度，而且具有很强的通用性，可以根据不同的需求选择相应的集成策略、EMO技术、基学习器等。