CN107423503A - 基于核慢特征分析和时滞估计的gpr建模 - Google Patents

Abstract

本发明公开了一种基于核慢特征分析和时滞估计的GPR建模方法。用于具有时延和非线性的化工过程。该方法首先通过模糊曲线分析充分挖掘工业数据中的时延信息，求出数据中的最优时滞，并进行建模数据的重构；进一步通过核慢特征分析方法对重构数据进行非线性的特征提取；最后，基于提取后的特征建立高斯过程回归模型，实现对关键变量的精确预测，从而提高产品质量，降低生产成本。

Description

基于核慢特征分析和时滞估计的GPR建模

技术领域

本发明涉及核慢特征分析和时滞估计的GPR建模，属于复杂工业过程建模和软测量领域。

背景技术

在实际工业过程中，一些关键变量的测量对于生产高质量的产品至关重要。但是在现有技术条 件和经济代价等问题的约束下，使得直接获取关键变量变得十分困难。

基于这样的背景下，软测量技术应运而生。它通过构建过程易测变量和关键变量的函数关系， 来推断和估计难测的关键变量，因此得到了广泛应用。常见的软测量建模方法如偏最小二乘、神经网络、 最小二乘支持向量机等都可以得到良好的预测效果。高斯过程回归(Gaussian process regression，GPR)由于 其具有预测精度高，超参数容易确定，输出具有概率意义等优点，近几年被广泛应用于软测量建模。因此 选取高斯过程回归作为建模方法。

实际工业过程建模中，一般都会对过程数据预先进行特征提取。现有的线性特征提取方法，如 主成分分析和独立主元分析，以及非线性特征提取方法核主元分析，都是从变量的角度构造新的成分，并 未分析样本之间的变化情况。慢特征分析(slow featureanalysis，SFA)作为一种新的特征提取方法，能够从 样本的角度出发，提取出一部分最缓慢的特征代表过程的本质特征，基于本质特征建立精度更高模型。近 年来，SFA得到了各个研究领域的关注，并且取得了广泛的应用，例如盲信号分离、软测量建模、行为识 别、故障诊断以及遥感。然而现有的慢特征分析方法往往是基于二阶多项式扩展，其处理非线性的能力较 弱。当数据的维度较大时，二阶多项式扩展会引起维数灾难，也有可能会造成数据的过拟合。针对工业过 程中的强非线性和高维特征，通过引入核函数扩展的核慢特征分析(kernel slow feature analysis，KSFA)代替 二阶多项式扩展的慢特征分析，加强了对非线性数据的特征提取能力。

实际的生产过程中输入输出数据会存在显著的时间滞后，如物料流动引起的传输时延、信号传 输带来的纯时延等。若对这些时滞置之不理，会导致输入和输出变量之间的因果关系发生变化，从而导致 建模精度的下降。因此，如何准确估计输入和输出之间的时延变得尤为重要。因此，如何准确估计输入和 输出之间的时延变得尤为重要。由于实际工业过程工艺复杂，通过过程工艺估计时延变得不太可能。Zhang 和Komulaine通过构建输入变量和输出变量的相关系数来估计时延，但是由于该方法只考虑了变量之间的 线性关系，对于非线性情况可能会得到错误的结论。阮宏镁等通过差分估计(differenceestimation，DE)算 法优化过程变量间的联合互信息来确定时延，但是智能优化算法容易陷入局部最优，且联合互信息分析相 关性时计算复杂度较高，运行时间较长。模糊曲线分析(fuzzy curve analysis，FCA)方法最早用于确定对输 出变量重要的输入变量，这种方法的优点在于计算复杂度较低，能够直观有效的确定出输入变量对输出变 量的重要性程度。因此本文将采用该方法确定最优时滞。通过引入不同的时滞信息扩展输入变量，构造输入变量和输出变量之间的模糊逻辑，从而确定对输出变量最重要的输入变量所对应的时滞值，此时的时滞 值则为该输入变量的最优时滞。

一种基于核慢特征分析和时滞估计的高斯过程回归建模方法(简称FCA-KSFA-GPR)能够很好的 处理时延和非线性，提高软测量模型的性能。通过FCA挖掘出各辅助变量的最优时滞，能够确定输入变 量对输出变量的重要性程度，因此使输入变量和输出变量的因果性更强，能够更好的解决时延问题。相对 与传统的特征提取方法，SFA从样本的角度出发，更能挖掘出本质特征，建立精度更高模型。通过KSFA 对重构数据进行特征提取，更能加强对非线性数据的特征提取能力。

发明内容

针对工业过程中存在的时延和非线性特性，提出了一种考虑过程时滞信息，并基于核慢特征分 析的高斯过程回归建模方法。

通过模糊曲线分析充分挖掘工业数据中的时延特征，求出数据中的最优时滞，并进行建模数据 的重构；进一步通过核慢特征分析方法对重构数据进行非线性的特征提取；最后，基于提取后的特征建立 高斯过程回归模型，实现对化工过程产品质量的在线估计，从而提高产量，降低生产成本。

本发明的目的是通过以下技术方案实现的：

基于核慢特征分析和时滞估计的GPR建模方法，所述方法包括以下过程：针对于时延和非线性 的化工过程，用模糊曲线分析对过程数据进行分析，确定出各输入变量的最优时滞，并进行数据的重构。

通过核慢特征分析对重构数据进行非线性特征提取，通过交叉验证，确定出最优的核慢特征特 征向量阵和对应的核慢特征，并通过核慢特征建立高斯过程回归模型。

对于新来的样本，通过训练样本得到的核慢特征特征向量阵求取新来样本的核慢特征，代入已 建立好的高斯过程回归模型中，即可实现对关键变量的精确预测，从而提高产品质量，降低生产成本。

附图说明

图1慢特征分析流程图；

图2重构后的数据；

图3基于核慢特征分析和时滞估计的GPR建模流程图；

图4前6个输入变量不同时延所对应的RC值；

图5训练样本和测试样本提取的四个变化缓慢的核慢特征；

图6是四种建模方法的预测结果图

图7四种建模方法的散点图

具体实施方式

下面结合图3所示，对本发明做进一步详述：

以常见的化工过程——脱丁烷塔过程为例。实验数据来自于脱丁烷塔过程，对脱丁烷塔底部的 丁烷浓度进行预测。

步骤1：收集输入输出数据组成历史训练数据库。

步骤2：标准化训练样本数据，并将每个输入变量进行T_max+1维扩展。其中，T_max为最大时延。

步骤3：通过模糊曲线分析(fuzzy curve analysis，FCA)确定各输入变量的最优时滞，定义为 d₁,d₂,...,d_m。其中m为样本的维度。所述的模糊曲线分析算法为：

FCA用于选择重要的输入变量。通过构建输入和输出变量之间的模糊逻辑，从而确定对输出变 量重要的输入变量。假设有m个输入变量{x_i,i＝1,2,...,m}和输出变量y，采集到的样本个数为n，令x_ik为 第k个样本的第i个变量，y_k为对应的输出变量。对于每个输入变量x_i,i＝1,2,...,m在第k个样本处的模 糊隶属度函数φ_ik(x_i)定义为：

式中，φ_ik为变量x_i所对应的第k个样本点的输入变量模糊隶属函数。φ_ik这里选取为高斯函数，b选取为 输入变量x_i范围的20％。因此对于n个训练样本对应的每个输入变量都有n条模糊规则，其中第k个规则 描述为：R^k:if x_i isφ_ik(x_i),then y is y_k。

通过式(2)对输入变量x_i进行重心去模糊化，计算得到其对应的模糊曲线为c_i。进一步通过式 (3)，可以计算出模糊曲线c_i的变化范围Rc_i。若得到的Rc_i越大，表明此时输入变量对于输入变量越重要。

Rc_i＝max(c_i)-min(c_i) (3)

由于共有m个输入变量，故能得到m个Rc_i，将Rc_i大小按降序排列，可以得到各输入变量的重要性。

通过引入时滞信息，原有的输入变量x_i转变为T_max+1维： {x_i(t),x_i(t-1),…,x_i(t-d),…,x_i(t-T_max)}，其中t为输出变量的采样时刻，d为待确定的时滞变量，T_max为最大的变量时延值。根据模糊曲线法原理，将{x_i(t),x_i(t-1),…,x_i(t-d),…,x_i(t-T_max)}转化为T_max+1个 待选输入变量，从中选取一个最重要的变量x_i(t-d)，则d即为该输入变量所对应的最优时滞。求取出各 输入变量的最优时滞分别为d₁,d₂,...,d_m后，则进行样本数据的重构。

对重构后的训练样本进行标准化处理，并通过核慢特征分析(kernel slowfeature analysis，KSFA) 确定其核慢特征和相应特征向量矩阵。慢特征分析算法为：

给定一个p维的时序输入信号X(t)＝{x₁(t),x₂(t),...,x_p(t)}，SFA算法旨在找到一个q维的变换 函数G(x)＝{g₁(X),g₂(X),...g_q(X)}，使得q维的输出信号S(t)＝{S₁(t),S₂(t),...,S_q(t)}变化尽可能缓慢。 其中S_j(t):＝g_j(X(t))，j∈[1,2,...,q]。SFA算法最终转化为一个优化问题：

且满足约束条件：

<S_j>_t＝0 (5)

<(S_j)²>_t＝1， (6)

其中为S_j的一阶导数，<·>_t为时间平均。

式(4)为优化目标函数，使得输出信号变化最为缓慢；式(5)表示的零均值约束，仅仅是为了简化 问题的求解；单位方差约束(6)不但对输出信号的尺度进行了归一化，还避免了常值解的出现，使输出信号 必须包含一定的信息；式(7)保证了输出信号的各个分量之间是不相关的，避免了冗余信号的出现。除此之 外，慢特征根据其变化快慢进行排序，若按升序排序，即S₁(t)是最缓慢的特征，S₂(t)是次缓慢的特征， 以此类推。

当变换函数g_j(·)为线性时，此时每个慢特征S_j都是输入变量的线性组合，即：

S_j＝g_j(X)＝X(t)w_j (8)

由于X(t)已经预先进行了均值化处理，因此满足式(5)的约束条件。式(4)优化目标和约束条件(6)可以 重写为：

式中，B＝<X(t)^TX(t)>_t，进而公式(9)的优化求解问题等价于如下的广义特征值分解 问题：

AW＝BWΛ (10)

其中，Λ＝(λ₁,...λ_q)为广义特征值构成的对角矩阵且λ₁＜λ₂＜…＜λ_q，W＝(w₁,...w_p)为对应的广义特征 向量矩阵。式(4)的优化目标正是对角矩阵Λ中的主对角元素。

当变换函数g_j(·)为非线性时，此时的变换函数可以等效为线性变换函数在非线性空间的扩展， 多项式扩展为常用的扩展方式，以二阶多项式扩展为例，定义如下：

经过二阶扩展后，的维度由p维变成了L维，其中L＝p+p(p+1)/2。将进行均值化处理 后得V(x)，将V(x)替换式(8)中的X(t)，就将非线性情况转化成了线性情况。

核慢特征分析就是将核函数扩展方法代替多项式扩展从而加强对非线性数据的处理，从而实现 非线性特征提取。核函数的选取有多种选择，这里选择泛化和平滑估计能力强的高斯核函数。

将经过核函数扩展后的数据K代替二阶多项式扩展的V(t)，优化目标就转换为：

其中，式(10)的优化目标求解问题等价于如下的广义特征值分解问题：

A_kW＝B_kWΛ (15)

其中，表示的协方差阵，B_k＝<K^TK>_t表示K的协方差阵。

通过交叉验证确定核慢特征个数，并将此时的核慢特征向量矩阵和对应的核慢特征分别记为 W，S_tr。然后基于(S_tr,Y)建立高斯过程回归模型，其中Y表示重组训练样本中的输出变量。建立的GPR 模型为：

给定训练样本集X∈R^D×N和y∈R^N，其中X＝{x_i∈R^D}_i＝1…N，y＝{y_i∈R}_i＝1…N分别代表D 维的输入和输出数据，输入和输出之间的关系由公式(16)产生：

y＝f(x)+ε (16)

其中f是未知的函数形式，ε是均值为0，方差为的高斯噪声。对于一个新的输入x^*，相应 的概率预测输出y^*也满足高斯分布，其均值和方差如式(11)和(12)所示：

y^*(x^*)＝c^T(x^*)C^-1y (17)

式中c(x^*)＝[c(x^*,x₁),…,c(x^*,x_n)]^T是训练数据和测试数据之间的协方差矩阵。是训练数据之间 的协方差矩阵，I是N×N维的单位矩阵。c(x^*,x^*)是测试数据的自协方差。

GPR可以选择不同的协方差函数c(x_i,x_j)产生协方差矩阵Σ，只要选择的协方差函数能保证产 生的协方差矩阵满足非负正定的关系。这里选择高斯协方差函数：

式中，v控制协方差的量度，ω_d代表每个成分x^d的相对重要性。

对式(4)中的未知参数v,ω₁,…,ω_D和高斯噪声方差的估计，一般最简单的方法就是通过极大 似然估计得到参数

为了求得参数θ的值，首先将参数θ设置为一个合理范围内的随机值，然后用共轭梯度法得到优 化的参数。获得最优参数θ后，对于测试样本x^*，可以用式(17)和(18)来估计GPR模型的输出值。

对计算所得的所有相似系数s_i进行降序排列，建立局部模型时，只选择前L个相似系数较大的 数据。为了选择合适比例的建模数据，针对于TE化工过程，数据比例选择从10％逐渐增大到100％，最 后得到最佳的数据比例为70％。JITL进行数据选择时，不同的比例下CPU耗时和预测的精度如图2和图 3所示。

根据JITL选择的数据建立的局部GPR模型为：

给定训练样本集X∈R^D×N和y∈R^N，其中X＝{x_i∈R^D}_i＝1…N，y＝{y_i∈R}_i＝1…N分别代表D维的 输入和输出数据。输入和输出之间的关系由公式(10)产生：

y＝f(x)+ε (10)

其中f是未知的函数形式，ε是均值为0，方差为的高斯噪声。对于一个新的输入x^*，相应 的概率预测输出y^*也满足高斯分布，其均值和方差如式(11)和(12)所示：

y^*(x^*)＝c^T(x^*)C^-1y (11)

式中c(x^*)＝[c(x^*,x₁),…,c(x^*,x_n)]^T是训练数据和测试数据之间的协方差矩阵。是训练数据之间 的协方差矩阵，I是N×N维的单位矩阵。c(x^*,x^*)是测试数据的自协方差。

GPR可以选择不同的协方差函数c(x_i,x_j)产生协方差矩阵Σ，只要选择的协方差函数能保证产 生的协方差矩阵满足非负正定的关系。这里选择高斯协方差函数：

式中v控制协方差的量度，ω_d代表每个成分x^d的相对重要性。

对式(4)中的未知参数v,ω₁,…,ω_D和高斯噪声方差的估计，一般最简单的方法就是通过极大 似然估计得到参数

为了求得参数θ的值，首先将参数θ设置为一个合理范围内的随机值，然后用共轭梯度法得到优 化的参数。获得最优参数θ后，对于测试样本x^*，可以用式(11)和(12)来估计GPR模型的输出值。

对于新来的样本，通过训练样本得到的核慢特征特征向量阵求取新来样本的核慢特征，代入已 建立好的高斯过程回归模型中，即可实现对关键变量的精确预测，从而提高产品质量，降低生产成本。

图6是脱丁烷塔底部的丁烷浓度预测值和实际值拟合曲线，并且与GPR、FCA-GPR、FCA-SFA-GPR所建立的软测量模型进行了比较。由图可知，基于核慢特征分析和时滞估计的GPR建模能 够有效地预测脱丁烷塔底部的丁烷浓度。

Claims (2)

1.基于核慢特征分析和时滞估计的GPR建模方法，其特征在于，该方法步骤为：

步骤1：收集输入输出数据组成历史训练数据库

步骤2：标准化训练样本数据，并将每个输入变量进行T_max+1维扩展。其中，T_max为最大时延

步骤3：通过模糊曲线分析(fuzzy curve analysis，FCA)确定各输入变量的最优时滞，定义为d₁,d₂,...,d_m。其中m为样本的维度，所述的模糊曲线分析算法为：

通过构建输入和输出变量之间的模糊逻辑，从而确定对输出变量重要的输入变量；假设有m个输入变量{x_i,i＝1,2,...,m}和输出变量y，采集到的样本个数为n，令x_ik为第k个样本的第i个变量，y_k为对应的输出变量；对于每个输入变量x_i,i＝1,2,...,m在第k个样本处的模糊隶属度函数φ_ik(x_i)定义为：

式中，φ_ik为变量x_i所对应的第k个样本点的输入变量模糊隶属函数；φ_ik这里选取为高斯函数，b选取为输入变量x_i范围的20％；因此对于n个训练样本对应的每个输入变量都有n条模糊规则，其中第k个规则描述为：R^k:if x_i isφ_ik(x_i),then y is y_k；

通过式(2)对输入变量x_i进行重心去模糊化，计算得到其对应的模糊曲线为c_i；进一步通过式(3)，可以计算出模糊曲线c_i的变化范围Rc_i；若得到的Rc_i越大，表明此时输入变量对于输入变量越重要；

Rc_i＝max(c_i)-min(c_i) (3)

由于共有m个输入变量，故能得到m个Rc_i，将Rc_i大小按降序排列，可以得到各输入变量的重要性；

通过引入时滞信息，原有的输入变量x_i转变为T_max+1维：{x_i(t),x_i(t-1),…,x_i(t-d),…,x_i(t-T_max)}，其中t为输出变量的采样时刻，d为待确定的时滞变量，T_max为最大的变量时延值；根据模糊曲线法原理，将{x_i(t),x_i(t-1),…,x_i(t-d),…,x_i(t-T_max)}转化为T_max+1个待选输入变量，从中选取一个最重要的变量x_i(t-d)，则d即为该输入变量所对应的最优时滞；求取出各输入变量的最优时滞分别为d₁,d₂,...,d_m后，则进行样本数据的重构

步骤4：对重构后的训练样本进行标准化处理，并通过核慢特征分析(kernel slowfeature analysis，KSFA)确定其核慢特征和相应特征向量矩阵；慢特征分析算法为：

给定一个p维的时序输入信号X(t)＝{x₁(t),x₂(t),...,x_p(t)}，SFA算法旨在找到一个q维的变换函数G(x)＝{g₁(X),g₂(X),...g_q(X)}，使得q维的输出信号S(t)＝{S₁(t),S₂(t),...,S_q(t)}变化尽可能缓慢；其中S_j(t):＝g_j(X(t))，j∈[1,2,...,q]；SFA算法最终转化为一个优化问题：

且满足约束条件：

<S_j>_t＝0 (5)

<(S_j)²>_t＝1， (6)

其中为S_j的一阶导数，<·>_t为时间平均

式(4)为优化目标函数，使得输出信号变化最为缓慢；式(5)表示的零均值约束，仅仅是为了简化问题的求解；单位方差约束(6)不但对输出信号的尺度进行了归一化，还避免了常值解的出现，使输出信号必须包含一定的信息；式(7)保证了输出信号的各个分量之间是不相关的，避免了冗余信号的出现；除此之外，慢特征根据其变化快慢进行排序，若按升序排序，即S₁(t)是最缓慢的特征，S₂(t)是次缓慢的特征，以此类推；

当变换函数g_j(·)为线性时，此时每个慢特征S_j都是输入变量的线性组合，即：

S_j＝g_j(X)＝X(t)w_j (8)

由于X(t)已经预先进行了均值化处理，因此满足式(5)的约束条件，式(4)优化目标和约束条件(6)可以重写为：

式中，B＝<X(t)^TX(t)>_t，进而公式(9)的优化求解问题等价于如下的广义特征值分解问题：

AW＝BWΛ (10)

其中，Λ＝(λ₁,...λ_q)为广义特征值构成的对角矩阵且λ₁＜λ₂＜…＜λ_q，W＝(w₁,…w_p)为对应的广义特征向量矩阵，式(4)的优化目标正是对角矩阵Λ中的主对角元素。

当变换函数g_j(·)为非线性时，此时的变换函数可以等效为线性变换函数在非线性空间的扩展，多项式扩展为常用的扩展方式，以二阶多项式扩展为例，定义如下：

经过二阶扩展后，的维度由p维变成了L维，其中L＝p+p(p+1)/2，将进行均值化处理后得V(x)，将V(x)替换式(8)中的X(t)，就将非线性情况转化成了线性情况；核慢特征分析就是将核函数扩展方法代替多项式扩展从而加强对非线性数据的处理，从而实现非线性特征提取；核函数的选取有多种选择，这里选择泛化和平滑估计能力强的高斯核函数：

将经过核函数扩展后的数据K代替二阶多项式扩展的V(t)，优化目标就转换为：

其中，式(10)的优化目标求解问题等价于如下的广义特征值分解问题：

A_kW＝B_kWΛ (15)

其中，表示的协方差阵，B_k＝<K^TK>_t表示K的协方差阵

步骤5：通过交叉验证确定核慢特征个数，并将此时的核慢特征向量矩阵和对应的核慢特征分别记为W，S_tr；然后基于(S_tr,Y)建立高斯过程回归模型，其中Y表示重组训练样本中的输出变量，建立的GPR模型为：

给定训练样本集X∈R^D×N和y∈R^N，其中X＝{x_i∈R^D}_i＝1...N，y＝{y_i∈R}_i＝1...N分别代表D维的输入和输出数据，输入和输出之间的关系由公式(16)产生：

y＝f(x)+ε (16)

其中f是未知的函数形式，ε是均值为0，方差为的高斯噪声，对于一个新的输入x^*，相应的概率预测输出y^*也满足高斯分布，其均值和方差如式(17)和(18)所示：

y^*(x^*)＝c^T(x^*)C^-1y (17)

式中c(x^*)＝[c(x^*,x₁),…,c(x^*,x_n)]^T是训练数据和测试数据之间的协方差矩阵，是训练数据之间的协方差矩阵，I是N×N维的单位矩阵，c(x^*,x^*)是测试数据的自协方差。

GPR可以选择不同的协方差函数c(x_i,x_j)产生协方差矩阵Σ，只要选择的协方差函数能保证产生的协方差矩阵满足非负正定的关系，这里选择高斯协方差函数：

式中v控制协方差的量度，ω_d代表每个成分x^d的相对重要性；

对式(19)中的未知参数v,ω₁,…,ω_D和高斯噪声方差的估计,采用极大似然估计得到参数

为了求得参数θ的值，首先将参数θ设置为一个合理范围内的随机值，然后用共轭梯度法得到优化的参数，获得最优参数θ后，对于测试样本x^*，可以用式(17)和(18)来估计GPR模型的输出值。

2.根据权利要求1所述的基于核慢特征分析和时滞估计的GPR建模方法，其特征在于，先通过模糊曲线分析，求出数据中的最优时滞；进一步通过核慢特征分析方法对重构数据进行非线性的特征提取；最后，基于提取后的特征建立高斯过程回归模型，实现对关键变量的精确预测，从而提高产品质量，降低生产成本。