间歇过程生产方式已被广泛应用于精细化工、生物制药、食品加工等诸多领域[1-2]．由于涉及非常复杂的物理化学反应，并受到外部不确定因素的极大干扰，间歇过程生产的最终产品质量难以保持稳定，生产现场须要及时调整工艺参数．然而，最终产品质量只能在生产结束后进行离线测量，质量反馈严重滞后，且需要较高的人力和时间成本，设计准确性高、解释性强的软测量方法已成为必然选择．软测量是以最终产品质量作为输出，并选取与其相关且易测的变量作为输入，构造某种数学关系来实现对最终产品质量的预测和估计[3]．现有研究工作大致分为第一性原理方法和数据驱动方法．前者通过物理、化学实验或仿真，并借助统计学模型来推理实际工业过程的因果关系．例如，文献[4]针对注塑成型这一典型的间歇过程，通过分析聚合物熔体特性来在线监测产品重量的变化．显然，只有深刻理解其中的物理化学规律，该类方法才能很好地发挥作用．随着数据采集和机器学习等技术的广泛应用，数据驱动方法已成为软测量研究的主流[5]，其重点在于挖掘数据中的关联关系，而不试图反映实际工业过程的因果关系．例如，偏最小二乘(PLS)、主成分分析(PCA)、支持向量回归(SVR)等浅层学习模型，以及多层感知器(MLP)、长短期记忆网络(LSTM)、卷积神经网络(CNN)等先进的深度学习模型[6]．文献[7]指出数据和特征决定了机器学习的上限，而模型和算法只是逼近这个上限．基于机器学习的软测量建模，其本质上是对数据中模式的识别，只有数据中包含最终产品质量的因果信息，算法的准确性和解释性才能得以保证．因此，必须进行充分的数据分析尤其是因果分析，以获得关于工业过程本身的感性见解，进而选择最终产品质量的因果特征作为输入来支撑软测量建模．否则，盲目改进机器学习算法，带来的是复杂的模型结构和大量的超参数，这与奥卡姆剃刀原理和工业领域的可靠性需求相悖．鉴于此，本研究提出基于非线性因果效应的间歇过程软测量方法，对间歇过程变量与最终产品质量之间的因果效应进行量化，并结合间歇过程的多时段特点提出堆叠集成AdaBoost算法，实现间歇过程最终产品质量的软测量．1 间歇过程数据驱动的建模策略间歇过程具有多时段特性，例如注塑成型过程，由合模、注射、保压、熔胶、冷却、开模等多个时段组成，过程变量在同一时段的数据分布是相似的，在不同时段则具有不同的变量轨迹、运行模式及相关性特征[8]．考虑采样频率的不同，间歇过程历史数据集中包含两种类型的原始数据，如图1所示．a．时段无关变量Vu：包含L个变量(l =1，2，…，L)和I个批次运行(i=1，2，…，I)的二维数组．这类变量在批次方向上是非线性动态变化的，但在同一批次的时间方向上是不变的，如某一关键时刻的温度和压力，某一关键时段的持续时间，设备状态的设定参数等．b．时段相关变量Vr，k：包含J个变量( j=1，2，…，J)、K个时间间隔(k=1，2，…，K)和I个批次运行(i=1，2，…，I)的三维数组．这类变量如温度、压力、流量等随生产过程在时间和批次方向上都是非线性动态变化的．10.13245/j.hust.220605.F001图1间歇过程历史数据集最终产品质量QFP取决于时间方向上的整体运行水平，是每个时间点对最终产品质量影响的累积，可以表示为时段无关变量以及时段相关变量在每个时间点产生影响的加权叠加，QFP=ω0Q0(Vu)+∑k=1KωkQk(Vr,k)，(1)式中：QFP为某一批次的最终产品质量；Q0，Qk分别为Vu和Vr对最终产品质量的局部影响；ω0，ωk分别为局部影响的权值．本研究提出的间歇过程数据驱动的建模策略，首先在每个时段提取时段相关变量的均值、方差、极差、中位数、偏度和峰度等时序特征，如表1所示．10.13245/j.hust.220605.T001表1时序特征表示时序特征计算公式均值∑k=1Nxj(k)/N方差∑k=1N(xj(k)-x¯j)2/(N-1极差max(xj)-min(xj)中位数median(xj)偏度∑k=1N(xj(k)-x¯j)3/N∑k=1N(xj(k)-x¯j)2/(N-1)1.5峰度∑k=1N(xj(k)-x¯j)4/N∑k=1N(xj(k)-x¯j)2/(N-1)2然后，从时段无关变量和每个时段的时序特征中选择最终产品质量的因果特征集，分别记为Su，Sr，p(p=1，2，…，P)，其中P为时段．如图2所示，为表述方便，将时段无关变量产生的影响记为时段0．由此，式(1)可简化为QFP=ω0Q0(Su)+∑p=1PωpQp(Sr,p ．(2)10.13245/j.hust.220605.F002图2间歇过程抽象描述示意图最后，通过堆叠集成AdaBoost算法从数据中学习得到每个时段的局部影响Q，并获得局部影响的权值ωk，代入式(2)中实现最终产品质量的预测．2 非线性因果效应与特征选择特征选择是从给定的候选特征集中选择一个特征子集作为机器学习算法的输入，以帮助理解数据，实现特征降维，降低算法的复杂度[9]．总的来说，主要涉及以下两个关键环节．a. 子集搜索策略．给定一个包含M个输入特征的集合F = {X1，X2，…，XM}，则有2M个候选输入特征子集S ⊇ F，其目标是从F中选择一个最优特征子集S[10]．先用空集初始化S，在每次迭代中，根据子集评价标准从F中选择一个特征加入S，直至达到停止阈值，这种策略称为前向搜索，即S=S⋃{EC(S⋃Xi,Y)TS,Xi∈F\S}，(3)式中：EC(·)为子集评价标准；TS为停止阈值；Xi (i = 1，2，…，M)为候选输入特征；Y为输出特征．也可以先初始化S = F，每次迭代从S中删除一个特征，直至达到停止阈值，这种策略称为后向搜索，即S=S\{EC(S\Xi,Y)TS,Xi∈S}．(4)显然，上述策略是贪婪的，若非穷举搜索则只能实现局部最优．此外，难以确定具有良好解释性和理论基础的最佳评价标准和停止阈值．b. 子集评价标准．在每次迭代中，用来判断是否保留或去除候选特征，如发散程度、关联程度及机器学习性能等．将子集搜索策略与子集评价标准相结合，即可得到不同的特征选择方法，大致分为过滤式、包裹式和嵌入式[11-12]三类．过滤式特征选择采用方差、皮尔逊相关系数(PCC)或最大信息系数(MIC)作为子集评价标准．方差可以用来衡量特征本身的发散程度，但不考虑输入输出之间的关联关系；PCC可以度量输入输出之间的线性关联，而MIC可以检测两个特征之间的非线性关联，但其需要大量的样本[13]．这类方法不依赖于机器学习算法，本身具有较强的解释性．包裹式特征选择采用某些机器学习算法本身所具有的参数作为子集评价标准，如回归系数和决策树的信息增益等可以衡量特征的权重或重要程度．这类方法依赖于机器学习的训练过程，计算成本高，且其本质上是基于关联关系．嵌入式特征选择以机器学习性能的最大化作为子集评价标准，结合智能优化算法可以自动选择最优特征组合[14- 15]，其计算成本高，在工业应用中容易产生过拟合现象．本研究受后非线性因果模型(PNM)[16-17]这一因果推理算法的启发，结合信息论定义了变量间的非线性因果效应，进而实现全自动的因果特征选择．2.1　后非线性因果模型如图3所示，若X与Y具有因果效应，且X是原因变量，Y是结果变量，则PNM引入的基本假设为Y=f2(f1(X)+ε)，(5)式中：f1为原因变量X的非线性作用；f2为观测时的非线性变形(可逆)；ε为噪声或干扰，且ε ⊥ X．10.13245/j.hust.220605.F003图3后非线性因果模型示意图若给定一组原因变量{X1，X2，…，Xk}和结果变量Y，式(5)可推广为Y=f2(f1(X1,X2,…,Xk)+εk)．(6)为了检测另一个变量Xk+1是否为Y的原因变量，将Xk+1加入到式(6)，得到Y=f2(f1(X1,X2,…,Xk,Xk+1)+εk+1)．(7)若Xk+1的加入使噪声项减弱，则表示其包含Y的因果信息，Xk+1对Y的因果效应可以定义为[18]EXk+1→Y=12logσ2(εk)σ2(εk+1)，(8)式中σ2(·)为方差．该定义利用回归残差的方差定义了因果效应，缺点在于须要建立两次非线性回归模型，时间复杂度较高，且假设后非线性因果模型的数据生成机制比较强．2.2　基于信息论的非线性因果效应度量在信息论中，一个离散随机变量X所包含的不确定性可以用香农熵来度量，即H(X)=-∑xP(x)logP(x)，(9)式中：H(·)为香农熵；P(x)为概率质量函数．两个离散随机变量X和Y所包含的总的不确定性可以用联合熵来度量，即H(X,Y)=-∑x,yP(x,y)logP(x,y)．(10)若X中的信息可以使Y中的不确定性降低，采用条件熵来度量给定X后，则Y中的剩余不确定性为H(Y|X)=H(X,Y)-H(X)．(11)类似地，给定一组原因变量{X1，X2，…，Xk}，结果变量Y中的剩余不确定性可以表示为      H(Y|X1,X2, …,Xk)=H(X1,X2, …,Xk,Y)-  H(X1,X2,…,Xk)．(12)进一步，将Xk+1加入到式(12)，得到H(Y|X1,X2, …,Xk,Xk+1)=                  -H(X1,X2, …,Xk,Xk+1)．(13)若Xk+1的加入使得Y中的剩余不确定性进一步减小，则Xk+1对Y具有因果效应，可得到式(8)的改进形式EXk+1→Y=H(Y|X0,X1,…,Xk)- H(Y|X0,X1, …,Xk,Xk+1)．(14)非线性因果效应度量方法其思路来源于后非线性因果模型，推导过程也可以用图4所示的非线性因果效应示意图表示．10.13245/j.hust.220605.F004图4非线性因果效应示意图2.3　基于非线性因果效应的特征选择当进行特征选择时，本研究采用前向搜索策略，以非线性因果效应作为子集评价标准，提出基于非线性因果效应的特征选择方法，形式化表示为S=S⋃{EXi→Y≠0,Xi∈F\S}. (15)与现有的特征选择方法不同，本文方法只须要按照一定的顺序遍历所有候选特征Xi，不须要设置停止阈值，自动地选择因果效应不为零的输入特征组合，算法1主要流程伪代码如下．算法1 基于非线性因果效应的特征选择方法输入 包含M个输入特征的集合F = {X1，X2，…，XM}，输出特征Y．输出 因果特征集S初始化S = ∅；特征离散化处理for i = 1 to M do计算Xi对Y的因果效应EXi→Yif EXi→Y≠0，then S=S⋃{Xi}end for返回S．3 堆叠集成AdaBoost软测量决策树回归是一类常见的机器学习方法，主要指分类回归树(CART)算法，可用于解决分类或回归问题．在此以决策树回归作为基本的学习器，设计堆叠集成AdaBoost软测量算法．3.1　决策树回归模型给定训练数据集D = {(x1，y1)，(x2，y2)，…，(xN，yN)}，当应用CART算法解决回归问题时，基于二分递归切分思想，采用平方误差最小化准则选出最优切分变量j及切分点s，即求解minj,sminc1∑xi∈R1(j,s)(yi-c1)2+minc2∑xi∈R2(j,s)(yi-c2)2. (16)然后通过该变量及切分点将输入空间划分为两个区域R1和R2，即由该节点生成两个子节点，分别含有N1和N2个样本，R1(j,s)={x|x(j)≤s};R2(j,s)={x|x(j)s}. (17)并确定这两个区域上的最优输出值ĉ1和ĉ2，ĉ1=1N1∑xi∈R1(j,s)yi; 　ĉ2=1N2∑xi∈R2(j,s)yi. (18)依次递归该过程直至满足结束条件，最终将输入空间划分为W个区域R1，R2，…，RW，生成决策树f(x)=∑w=1WĉwI(x∈Rw), (19)式中I(·)为指示函数，若x∈Rw，则I=1，否则I=0．生成回归树后，从底部开始由下往上依次进行剪枝，直至根节点．对于每一种剪枝情况都会生成一颗子树，由此形成一个子树序列f1(x)，f2(x)，…，fn(x)；然后在独立的验证数据集上使用交叉验证的方法，比较各子树对于验证集的平方误差，选出最优决策树fα(x)．3.2　堆叠集成AdaBoost算法在批次方向上，由于多种工况或产品规格的存在，批产过程通常具有多模态特性，使用单一的机器学习模型很难提高预测性能．AdaBoost集成策略的核心思想在于使每一个单一机器学习模型在训练过程中依次更加关注之前的错误样本，以保证模型之间的差异性，使得最终预测结果更加准确、可靠，算法主要流程如算法2所示．算法 2 堆叠集成AdaBoost软测量算法输入 含有N个样本的训练集D = {(x1，y1)，(x2，y2)，…，(xN，yN)}, 基本的学习器(本研究采用决策树回归模型)，迭代次数T．输出 软测量模型f (x)初始化：W1 = {w1(1)，w1(2)，…，w1(N)}，w1(i) = 1/N，i = 1，2，…，N．for t = 1 to T do用分布Wt = {wt(1)，wt(2)，…，wt(N)}从D 中取一个样本集Rt；用Rt训练基本的学习器，并构建回归模型ft(x)；计算每个训练样本的损失函数Lt(i)；计算平均损失L¯t；设置αt，将Wt更新为Wt+1={wt+1(1)，wt+1(2)，…，wt+1(N)}．end for返回 f (x)．给定D={(x1，y1)，(x2，y2)，…，(xN，yN)}为训练集，Wt={wt(1)，wt(2)，…，wt(N)} (t=1，2，…，N)表示第t次迭代时D上的权值分布．在以后的迭代中，Wt将通过增加性能较差样本的权重、降低性能较好样本的权重进行更新．其中性能可采用平均损失函数来评估，L¯t=∑i=1NLt(i)wt(i) , (20)式中Lt为范围为0到1的损失函数．文献[19]给出了3个候选Lt，本研究采用指数形式，即Lt(i)=1-exp{lt(i)/max(lt(i))}(i=1,2,…,N), (21)式中lt(i) = |ft(xi)-yi|为每个训练样本的损失．权值更新程序如下，wt+1(i)=wt(i)αt1-Lt(i)/Zt, (22)式中：αt=L¯t/(1-L¯t)为权值更新参数；Zt为使得Wt+1成为一个概率分布的规范化因子．最终，AdaBoost回归的结果表示为f(x)=-∑t=1Tft(x)logαt . (23)4 实验验证与结果分析为了验证上述方法在批产过程软测量中的有效性，在此以注塑成型过程为案例进行对比实验．4.1　数据集描述与实验设置注塑成型过程是利用注塑机将塑料原料在高温下熔化，通过高速、高压将塑料熔体注射到模具中，在恒定压力下经历复杂的物理化学变化，最终得到塑料制品．通过该过程的反复操作，可以生产出大量相同的产品．在此过程中，最终产品质量的测量存在较大的延迟，严重影响了确保质量稳定的及时决策，亟需准确可靠的软测量技术[20]．本研究利用第四届工业大数据创新竞赛中富士康集团提供的1.66×104个生产批次样本(数据下载网址：https://www.industrial-bigdata.com)，包括86个时段无关变量和22个时段相关变量，以产品尺寸作为最终产品质量进行实验验证，并以其中的9 960个样本作为训练集，其余样本作为测试集．为了验证提出的软测量方法的有效性，在分时段因果特征选择的基础上，以PLS，SVR和决策树作为基准方法进行对比验证．在上述模型的超参数选择过程中，本研究仅在默认参数的基础上进行微调：PLS的成分数量设为20；SVR均采用默认参数；决策树的最大深度设为10，最小样本分割设为5；AdaBoost含有20颗决策树，学习率设为1.3．针对本研究所采用的数据集和所提方法，上述模型参数对软测量性能影响较小，读者可自行验证．采用均方根误差(ERMS)和决定系数R2作为性能评价指标，ERMS=∑i=1NT(yi-ŷi)2/NT　 ; (24)R2=1-∑i=1NT(yi-ŷi)2/∑i=1NT(yi-y¯)2, (25)式中：NT为测试集中的样本数；yi为第i个样本的真实值；ŷi为软测量模型的估计值；y¯为所有估计值的平均值．本研究的所有代码在Python3.7中运行，计算机配置为Intel(R)Core(TM)i7-8700CPU@3.20 GHz 32.00 GiB RAM．4.2　对比结果与讨论在注塑成型过程的11个关键生产时段(合模、注射、保压、熔胶、后松退、冷却、开模、顶进、顶退、中子进和中子退)，分别提取22个时段相关变量的6个时序特征(见表1)，即每个生产时段衍生出132个候选输入特征，进而分时段计算候选输入特征对产品尺寸的因果效应并进行特征选择．在此基础上，采用不同时段的输入特征集训练得到11个软测量模型，再加上基于时段无关变量训练得到的软测量模型，采用线性回归确定这12个软测量模型各自的权重，最后基于式(2)得到注塑成型产品尺寸的最终估计结果．不同软测量建模方法的ERMS和R2对比结果见表2，可以看出：本文方法(AdaBoost)所得到的ERMS最低，R2最大，获得最佳性能；PLS的所得到的ERMS最大，R2最小，即性能最差．综合图5和图6的对比结果也可以看出，SVR、决策树和AdaBoost这三者之间的性能差异较小，这是因为基于非线性因果效应进行特征选择，在输入特征相同的情况下，不同机器学习模型的性能差异较小．10.13245/j.hust.220605.F005图5不同软测量方法的估计结果对比10.13245/j.hust.220605.T002表2不同软测量方法的ERMS和R2对比方法ERMS/mmR2/%PLS0.02576.7SVR0.02381.2决策树0.02085.5AdaBoost0.02085.610.13245/j.hust.220605.F006图6不同软测量方法估计结果的散点图及估计误差的概率密度曲线对比5 结语a. 在后线性因果模型的启发下，基于信息论定义了非线性因果效应度量方式，进而提出了基于非线性因果效应的特征选择方法．本方法能实现间歇过程因果特征的自动选择，无须设置停止阈值，避免信息损失；在实际软测量应用中，该方法具有良好的准确性，能够降低软测量结果的不确定性．b. 考虑多时段特性建立了间歇过程的最终产品质量模型，基于AdaBoost集成决策树方法实现多时段模型集成的质量在线预测．对比实验表明：本文方法与决策树方法的ERMS和R2最佳，PLS和SVR效果略差．本文方法模型结构简单、超参数较少，具有广阔的应用前景．