CN106156434A - 基于局部时滞重构的滑动窗时间差‑高斯过程回归建模方法 - Google Patents

Abstract

本发明涉及一种基于局部时滞重构的滑动窗时间差‑高斯过程回归建模方法，适合应用于具有时滞、非线性、时变特性的化工过程。本发明所述方法能够通过滑动窗口策略逐步跟踪最新的过程时变动态，同时，在滑动窗口中采用模糊曲线分析方法对过程阶段性时滞特征进行参数提取，并用于局部模型训练样本和测试样本的时滞重构；然后采用时间差高斯过程回归(TDGPR)模型来描述局部重构滑动窗口上的变量漂移特征。本发明为工业过程提供了一种有效的实时预测和控制的技术支持手段，有利于提高产品质量，控制生产成本，规避安全隐患。

Description

基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法

技术领域

本发明涉及基于局部时滞重构的滑动窗口时间差高斯过程回归(LTR-MWTDGPR)的软测量建模方法，属于复杂工业过程建模和软测量领域。

背景技术

随着现代工业过程对于产品质量控制和优化要求的不断加强，对过程变量的在线测量技术也相应地提出了更高的要求。在许多的实际应用中，关乎质量的变量(如气体流量浓度值，产品中某种化学成分含量等)往往很难在线进行测量，但亟待及时地进行获取。在这一背景下，软测量技术应运而生，它通过构建辅助变量集与主导变量的数学模型来推断难测主导变量的实时值，被广泛地运用于工业过程中。

在某些情况下，即使现场有安装好的主导变量在线分析仪表，但是由其测量分析得到的主导变量序列与采样得到的辅助变量时间序列值并不是一致的，呈现出显著的时间滞后特征，造成这一特征的原因有系统的容积时延，信号和物料的传输时延以及测量仪表的分析周期带来的时延等。如果在建模时忽略这样的时滞特性，那么模型精度和系统的控制质量就会大大降低。时延的增加不仅带来控制性能的下滑，而且建立的模型不能够解释实时过程状态，甚至带来灾难性的后果。因此，为了对化工生产过程进行优化控制，对辅助变量集与主导变量之间的时滞参数进行可靠的估计是非常必要的。

目前，相比于机理建模方法，基于数据驱动理论的方法由于不需要深入的先验过程知识，同时具有更好的泛化性能而广泛应用于软测量建模领域。数据驱动模型，顾名思义，基于过程大批量的历史采样数据而建立，随着DCS系统的快速发展，过程能够大批量地进行同步数据采集，由于过程仪表时空分布上的差异，过程变量与质量变量间的时滞依然存在。因此，实时采集的数据集中包含着十分有用的时滞信息，为建立时滞软测量模型提供了可能。

为了提取过程时滞信息，国内外文献中多采用过程硬件仪表的设计参数(如反应器容积，管道长度等)估计装置的大致时延范围，或者利用输入输出变量间的相关系数，将滞后的样本信息引入软测量模型，抑或是利用互信息指标选择合适的输入变量以提高软测量模型的可靠性。然而，在上述方法中，模型结构中的滞后样本个数选择依赖于反复试验方法，易出现不稳定的模型性能。基于相关系数分析的时延估计方法只对线性系统描述较好，基于互信息的估计算法往往具有很高的计算复杂度。尽管时延参数有时候能够通过先验知识或对过程机理的深入预分析而事先确定，但是对于较为复杂的系统，这样的方法往往伴随着很大的随机性和不确定性。因此，非常希望找到一种既能够有效追踪过程非线，又具有较低的计算复杂度的时滞估计方法。

当软测量模型建立并投入使用，很难避免模型品质下降的问题。为了降低性能下降程度。近年来自适应学习机制被广泛提出，使模型长期地维护其预测精度。最常用的方法有迭代法(RM)，滑动窗口法(MW)，实时学习方法(JITL)，以及时间差方法(TD)。在这些方法中，基于MW和RM的模型能够同时处理过程变量和质量变量的缓慢漂移，基于JITL策略的模型适用于过程变量的漂移情况。然而，在许多应用中，单独地使用以上每一种策略都需要反复地重建模型，而TD模型结构不仅可以同时对输入输出漂移进行建模具有高稳定性，同时不存在反复模型更新的问题。如同所有的离线模型一样，全局的TD模型也会随着时间而老化。

综上所述，在实际工业过程建模时，时变特性、过程非线性以及时滞特性都是需要着重处理的问题，任何一种特性都不能被忽视，为了对过程进行实时控制和优化，建立考虑以上三种特征的高效、高精度软测量模型是意义深远且至关重要的。

有鉴于上述的缺陷，本设计人积极加以研究创新，以期创设一种基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法，使其更具有产业上的利用价值。

发明内容

为解决上述技术问题，本发明的目的是提供一种不需要等待新输入数据的到来，仅需要利用相应的时滞重构样本作为模型输入进行TDGPR建模，从而提前获得主导变量的实时预测值，一定程度上解决了主导变量获得具有时间滞后的问题的基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法。

本发明基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法，包括：

S1估计过程主导变量y与辅助变量集X＝[x₁,x₂,...,x_m]间存在的最大时延参数T_max；

S2获取过程连续均匀采样输入输出变量集，建立包含L组连续时间样本的初始滑动窗口W_ini,W_ini＝[X(t),y(t)]_t＝1,...,L；

其中X(t)＝[x₁(t),x₂(t),...,x_m(t)]_t＝1,...,L,m为辅助变量的个数，L为正整数；

S3对初始滑动窗口内m个辅助变量进行扩展，得到m(T_max+1)维的时滞扩展变量集，扩展方式为：

采用模糊曲线分析方法对时滞扩展变量集的各扩展变量进行分析，得到模糊覆盖范围最大的m个时滞变量，将各自的最优时滞参数记为d₁,d₂,...,d_m；

S4利用最优时滞参数对初始滑动窗口内L组实时对应的辅助变量进行时滞校正，重构后的窗口数据个数变为L-T_max，重构窗口记为W_rec，其中

$W_{rec}={\[x_1(t-d_1),x_2(t-d_2),...,x_m(t-d_m),y\left(t\right)\]}_{t=T_{max}+1,...,L};$

S5当t_new采样时刻到来，无需采样辅助变量x(t_new)的情况，采用滑动窗口内的估计好的时滞参数对将到来的实时输入变量重构到d_i时刻之前作为查询样本进行建模，对x(t_new)进行时滞重构后得到的查询样本记为x_d(t_new)，这里的x_d(t_new)＝[x₁(t_new-d₁),x₂(t_new-d₂),...,x_m(t_new-d_m)]；

S6对窗口中重构的训练样本及重构查询数据求取j次时间差分，然后对时间差分输入输出训练样本建立GPR模型，即建立了TDGPR模型结构，通过该模型结构得到主导变量的动态漂移预测值Δy_j,pred(t_new),在历史主导变量实时测量值基础上求出x_d(t_new)对应的实时输出y_j,pred(t_new)，预测完成后，丢弃当前时滞重构的MWTDGPR模型，参见下式：

$\left\{\begin{array}{c}\mathrm{\Δ}{x}_{d,j}\left(t_{new}\right)=x_d\left(t_{new}\right)-x_d(t_{new}-j)\\ \mathrm{\Δ}{y}_{j,pred}\left(t_{new}\right)=f_{GPR}\left(\mathrm{\Δ}{x}_{d,j}\right(t_{new}\left)\right)\\ y_{j,pred}\left(t_{new}\right)=y_j(t_{new}-j)+\mathrm{\Δ}{y}_{j,pred}\left(t_{new}\right)\end{array}\right.$

S7将最新的样本更新至初始滑动窗口，剔除初始滑动窗口中最旧的样本，即窗口向前滑动一步，然后，重复上述步骤S1至S6。

进一步地，采用模糊曲线分析方法得到各自的最优时滞参数的具体方法为：

对于扩展的时滞窗口输入变量{x_i(t-λ),λ＝0,1,…,T_max}，假设初始滑动窗口内输入输出对应关系为x_i(t-λ)→y(t)，t＝T_max+1,...,L。那么，可将x_i(t-λ)的模糊隶属度函数定义为：

${\mathrm{\Φ}}_{it}\left(x_i\right(t-\mathrm{\λ}\left)\right)=\exp \[-{\left(\frac{I_n{x}_{it,\mathrm{\λ}}-x_i(t-\mathrm{\λ})}{b}\right)}^2\]$

Φ_it(x_i(t-λ))为变量x_i(t-λ)关于第t个数据点的输入变量模糊隶属度函数，b取变量x_i(t-λ)值域范围的20％，I_n＝[1,1,...,1]^T,n为序列x_i(t-λ)的长度，n＝L-T_max。x_it,λ表示x_i(t-λ)序列在t时刻的采样值，x_i(t-λ)为时延变量序列；

通过下式对扩展后的(T_max+1)个窗口新变量质心去模糊化，得第i个变量时延值为λ条件下的模糊曲线C_i,λ

$C_{i,\mathrm{\λ}}\left(\mathrm{\λ}\right)=\frac{{\mathrm{\Σ}}_{t=T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]\·y\left(t\right)}{{\mathrm{\Σ}}_{T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]}$

其中，C_i,λ(λ)_max和C_i,λ(λ)_min为模糊曲线上点值域的最大值和最小值，覆盖范围即为模糊曲线最大值与最小值之差，d_i为使模糊曲线C_i,λ覆盖范围最大的λ，也即

借由上述方案，本发明至少具有以下优点：

本发明搜集足够多的均匀采样数据，组成历史数据库,根据过程先验知识或者预分析方法确定过程参数，如各变量最大时延参数，滑动窗口长度等。利用滑动窗口策略来增加新样本丢弃旧样本，达到实时跟踪过程动态的目的，组建初始的离线滑动窗口，并对该初始滑动窗口利用模糊曲线分析方法进行时滞分析，提取窗口各变量的最优时滞参数进行时序重匹配。然后，采用时间差高斯过程回归模型结构来拟合重匹配后滑动数据库上的漂移特性，当新的采样时刻到来前，可利用窗口经过时滞重组的输入作为预测模型输入，来获得新采样时刻的主导变量预测值，这对于高精度的产品质量控制和实时过程监控都是有益的。

本发明采用模糊曲线分析(Fuzzy Curve Analysis,FCA)的方法将变量时滞信息引入到软测量模型中，这种方法的特点是计算复杂度较低同时易于理解，能够直观有效地确定输入变量的重要性程度。同时，考虑到过程控制中操作条件的时变特性，过程采集数据也呈现了显著的阶段特性。因此，在考虑时滞参数估计的同时，不同操作条件下的时滞和漂移特征都需要进行对症处理，以更好地进行局部时序匹配，同时也进一步提高软测量模型的可靠性。

为了确保全局TD模型能够不断适应过程缓变和突变，本发明采用滑动窗口、时间差分方法两种不同的自适应机制结合使用以提高模型对于时变非线性动态的可靠性。另外，局部模型的选择对于局部非线性特征的描述程度而言亦具有很大意义。一些常用的模型如偏最小二乘(PLS)、主成分分析(PCA)等能够很好地处理输入变量和输出变量之间的线性关系，人工神经网络(ANN)、支持向量机(SVM)、最小二乘支持向量机(LS-SVM)能够有效地处理过程的非线性关系。近年来，高斯过程回归(GPR)作为一种非参数概率模型，不仅可以给出预测值，还可以得到预测值的不确定程度。故本发明选择GPR模型描述局部动态特征，结合TD思想来有效地处理过程输入输出的漂移。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，并可依照说明书的内容予以实施，以下以本发明的较佳实施例并配合附图详细说明如后。

附图说明

图1为本发明基于局部时滞重构的滑动窗-时间差高斯过程回归(LTR-MWTDGPR)建模方法建模示意图；

图2为硫回收装置过程示意图；

图3为硫回收装置过程TDGPR,MWTDGPR及本发明方法在不同L情况下j与RMSE指标的关系；

图4为硫回收装置过程本发明在不同滑动窗口长度下RMSE指标变化趋势；

图5为脱丁烷塔过程示意图；

图6为脱丁烷塔过程TDGPR,MWTDGPR及本发明在不同L情况下j与RMSE指标的关系；

图7为脱丁烷塔过程本发明在不同滑动窗口长度下RMSE指标变化趋势。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

结合图1所示的建模流程示意图和2个实际工业过程案例，对本发明作进一步详述：

工业案例一：硫回收装置过程

硫回收装置过程(Sulfur recovery unit，SRU)是精炼厂处理系统的重要环节，负责对含硫气体(如H₂S和SO₂)的处理，以避免其对环境造成巨大的危害，具体的反应过程流程简要示意图如图2所示，过程有5个辅助变量，2个主导变量，过程数据来源中，7个变量均可由在线传感器或者在线分析仪表获得，本发明中选择H₂S浓度作为过程需要进行实时估计的主导变量进行建模。

步骤1：采集过程的前1000组连续时间样本，采样时间间隔为1min，由于对该过程输入输出间的时滞情况缺乏充分的先验认知，假设在均匀采样得到的数据充分多的情况下认为过程时滞特性越稳定，对该数据进行时滞预分析，设置最大时延参数分别为5～15间的整数，对800组数据通过模糊曲线分析方法(见步骤3)来得到相应的变量时滞参数估计结果，发现当T_max设置为15时，无论怎样地改变训练样本数量，时滞参数不再局部收敛，其估计结果最大值也不再超过15，因此将T_max设置为15；

步骤2：将滑动窗口的长度L分别设置为30，50，70，90，110。以滑动窗口长度为30为例，建立初始离线滑动窗口，记为W_ini＝[X(t),y(t)]_t＝1,...,30,其中输入矩阵记为X(t)＝{x₁(t),x₂(t),...,x₅(t)}_t＝1,...,30；

步骤3：在长度为30的滑动窗口上对原始5维变量进行扩展，得到80维的输入变量集，扩展方式为：

对于扩展的时滞窗口输入变量{x_i(t-λ),λ＝0,1,…,T_max}，进行模糊曲线分析，假设窗口内输入输出对应关系为x_i(t-λ)→y(t),x_i(t-λ)的模糊隶属度函数定义为：

${\mathrm{\Φ}}_{it}\left(x_i\right(t-\mathrm{\λ}\left)\right)=\exp \[-{\left(\frac{I_n{x}_{it,\mathrm{\λ}}-x_i(t-\mathrm{\λ})}{b}\right)}^2\]---\left(2\right)$

Φ_it(x_i(t-λ))为变量x_i(t-λ)关于第t个数据点的输入变量模糊隶属度函数，b取变量x_i(t-λ)值域范围的20％，I_n＝[1,1,...,1]^T,n为序列x_i(t-λ)的长度，n＝L-T_max＝15。x_it,λ表示x_i(t-λ)序列在t时刻的采样值，x_i(t-λ)为时延变量序列。通过式(3)对扩展后的16个窗口新变量质心去模糊化，可得第i个变量时延值为λ条件下的模糊曲线C_i,λ，d_i为使模糊曲线C_i,λ覆盖范围最大的λ，C_i,λ(λ)_max和C_i,λ(λ)_min为模糊曲线上点值域的最大值和最小值，覆盖范围即为模糊曲线最大值与最小值之差；

$C_{i,\mathrm{\λ}}\left(\mathrm{\λ}\right)=\frac{{\mathrm{\Σ}}_{t=T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]\·y\left(t\right)}{{\mathrm{\Σ}}_{T_{max}+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]}---\left(3\right)$

$d_i=\underset{\mathrm{\λ}}{\mathrm{argmax}}\[C_{i,\mathrm{\λ}}{\left(\mathrm{\λ}\right)}_{max}-C_{i,\mathrm{\λ}}{\left(\mathrm{\λ}\right)}_{\min }\]---\left(4\right)$

步骤4：利用上一步分析得到的时滞参数对窗口内30组实时对应的输入输出样本进行时滞校正，重构后的窗口数据个数变为L-T_max，即15组样本，重构窗口可记为W_rec：

W_rec＝[x₁(t-d₁),x₂(t-d₂),...,x₅(t-d₅),y(t)]_{t＝16,...,30}；

步骤5：当t_new采样时刻到来，无需采样辅助变量x(t_new)的情况，可采用滑动窗口内的估计好的时滞参数对将到来的实时输入变量重构到d_i时刻之前作为查询样本进行建模，对x(t_new)进行时滞重构后得到的查询样本可记为x_d(t_new)，这里的x_d(t_new)＝[x₁(t_new-d₁),x₂(t_new-d₂),...,x₅(t_new-d₅)]；

步骤6：对窗口中重构的训练样本及重构查询数据求取j次时间差分，然后对时间差分输入输出训练样本建立GPR模型，即建立了TDGPR模型结构，通过该模型结构可得到主导变量的动态漂移预测值Δy_j,pred(t_new)，最后在历史主导变量实时测量值基础上可以求出x_d(t_new)对应的实时输出y_j,pred(t_new)由(5)给出，预测完成后，丢弃当前时滞重构的MWTDGPR模型；

$\left\{\begin{array}{c}{\mathrm{\Δx}}_{d,j}\left(t_{new}\right)=x_d\left(t_{new}\right)-x_d(t_{new}-j)\\ {\mathrm{\Δy}}_{j,pred}\left(t_{new}\right)=f_{GPR}\left({\mathrm{\Δx}}_{d,j}\right(t_{new}))\\ y_{j,pred}\left(t_{new}\right)=y_j(t_{new}-j)+{\mathrm{\Δy}}_{j,pred}\left(t_{new}\right)\end{array}\right.---\left(5\right)$

步骤7：将最新的样本更新至原滑动窗口，剔除窗口中最旧的样本，进行下一次时滞估计，同时，等待新采样时刻的查询样本到来；不断重复上述步骤，可以提前实现H₂S浓度的实时预测。

当L分别设置为30,50,70,90,110大小时，将本发明方法分别与不考虑局部时滞的MWTDGPR方法，传统全局的TDGPR方法进行了不同时间差大小的对比，j＝1～9，如图3所示，可以发现，当j由1到9变化时，传统单一的TDGPR模型随着时间差增加性能恶化显著，而加入滑动窗口更新自适应机制的TDGPR模型稳定性和预测性能较TDGPR模型显著增强；

同时由图3可知，当L＝30时，本发明方法仅在j＝1和j＝2时有着最佳的性能，随着时间差增加，重构时滞样本对于局部动态的描述能力下降。这是因为此时窗口长度较短，在时间差阶次增大时，可供后续TDGPR模型训练的样本逐渐减少，而硫回收过程又伴随较明显的过程突变，导致建模数据信息量匮乏，因此此时TD模型精度下降。此外，当滑动窗口窗口长度逐渐增大时，可以注意到本发明方法的精度在时间差阶次增大时有了逐步提高；在L＝70和90时，由于既考虑了局部时滞，同时又能够给TDGPR模型带来较多有用信息的数据量，因此相对于传统的MWTDGPR模型可靠性和精度都有了明显的改善。最后，当窗口长度增加到110时，由于窗口长度较长，导致时滞估计的局部描述能力减弱，尽管TDGPR模型提供了较多的时序匹配样本，也不能显著地优于传统MWTDGPR方法。

5种窗口大小情况下的硫回收过程软测量模型RMSE指标趋势图，如图4所示。对于硫回收装置过程，采样间隔为1min，发现当j＝1～5时，取L＝30这样小的窗口大小预测精度最佳，这时既可以确保时滞估计涵盖局部特征，又能够追踪到变量动态漂移；若实际情况需要根据6～8分钟之前历史分析值来预测当前时刻H₂S浓度情况时，则需要选择如50～70左右的窗口大小，这样可以更好地在局部时滞和非线性动态中取得平衡，同时也证明了本发明方法对于实时预测H₂S浓度的有效性和精度。

工业案例二：脱丁烷塔过程

脱丁烷塔过程是石油炼制生产过程中脱硫和石脑油分离装置的重要组成部分，过程示意图如图5所示。发明所用的验证数据集为标准数据集，共包含样本2394组，过程有7个辅助变量，分别描述为：x₁塔顶温度；x₂塔顶压力；x₃塔顶回流量；x₄塔顶产品流出量；x₅第6层塔板温度；x₆塔底温度1；x₇塔底温度2，1个主导变量为塔底丁烷浓度，其值不能直接检测，需要通过在线仪表分析获得，在线仪表的测量周期和安装位置使得各辅助变量与主导变量间存在大约45～90min的滞后。为了方便起见，发明中将2个塔底温度变量取平均后作为1个辅助变量。

步骤1：以数据采样6min为间隔，取过程前800组数据用于软测量建模，由于各辅助变量与主导变量间存在大约45～90min的滞后，基于以上过程先验知识，设置T_max参数为19；

步骤2：将滑动窗口的长度L分别设置为30，50，70，90，110。以滑动窗口长度为30为例，建立初始离线滑动窗口，记为W_ini＝[X(t),y(t)]_t＝1,...,30,其中输入矩阵记为X(t)＝{x₁(t),x₂(t),...,x₆(t)}_t＝1,...,30；

步骤3：在长度为30的滑动窗口上对原始6维变量进行扩展，得到120维的输入变量集，扩展方式为：

对于扩展的时滞窗口输入变量{x_i(t-λ),λ＝0,1,…,T_max}，进行模糊曲线分析，假设窗口内输入输出对应关系为x_i(t-λ)→y(t),x_i(t-λ)的模糊隶属度函数定义为：

${\mathrm{\Φ}}_{it}\left(x_i\right(t-\mathrm{\λ}\left)\right)=\exp \[-{\left(\frac{I_n{x}_{it,\mathrm{\λ}}-x_i(t-\mathrm{\λ})}{b}\right)}^2\]---\left(7\right)$

Φ_it(x_i(t-λ))为变量x_i(t-λ)关于第t个数据点的输入变量模糊隶属度函数，b取变量x_i(t-λ)值域范围的20％，I_n＝[1,1,...,1]^T,n为序列x_i(t-λ)的长度,n＝11。x_it,λ表示x_i(t-λ)序列在t时刻的采样值，x_i(t-λ)为时延变量序列。通过式(8)对扩展后的20个窗口新变量质心去模糊化，可得第i个变量时延值为λ条件下的模糊曲线C_i,λ，d_i为式(9)使模糊曲线C_i,λ覆盖范围最大的λ，C_i,λ(λ)_max和C_i,λ(λ)_min为模糊曲线上点值域的最大值和最小值，覆盖范围即为模糊曲线最大值与最小值之差；

$C_{i,\mathrm{\λ}}\left(\mathrm{\λ}\right)=\frac{{\mathrm{\Σ}}_{t=T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]\·y\left(t\right)}{{\mathrm{\Σ}}_{T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]}---\left(8\right)$

$d_i=\underset{\mathrm{\λ}}{\mathrm{argmax}}\[C_{i,\mathrm{\λ}}{\left(\mathrm{\λ}\right)}_{max}-C_{i,\mathrm{\λ}}{\left(\mathrm{\λ}\right)}_{\min }\]---\left(9\right)$

步骤4：利用上一步分析得到的时滞参数对窗口内30组实时对应的输入输出样本进行时滞校正，重构后的窗口数据个数变为L-T_max，即11组样本，重构窗口可记为W_rec：

W_rec＝[x₁(t-d₁),x₂(t-d₂),...,x₆(t-d₆),y(t)]_{t＝20,...,30}；

步骤5：当t_new采样时刻到来，无需采样辅助变量x(t_new)的情况，可采用滑动窗口内的估计好的时滞参数对将到来的实时输入变量重构到d_i时刻之前作为查询样本进行建模，对x(t_new)进行时滞重构后得到的查询样本可记为x_d(t_new)，这里的x_d(t_new)＝[x₁(t_new-d₁),x₂(t_new-d₂),...,x₆(t_new-d₆)]；

步骤6：对窗口中重构的训练样本及重构查询数据求取j次时间差分，然后对时间差分输入输出训练样本建立GPR模型，即建立了TDGPR模型结构，通过该模型结构可得到主导变量的动态漂移预测值Δy_j,pred(t_new),最后在历史主导变量实时测量值基础上可以求出x_d(t_new)对应的实时输出y_j,pred(t_new)由(10)给出,预测完成后,丢弃当前时滞重构的MWTDGPR模型；

$\left\{\begin{array}{c}{\mathrm{\Δx}}_{d,j}\left(t_{new}\right)=x_d\left(t_{new}\right)-x_d(t_{new}-j)\\ {\mathrm{\Δy}}_{j,pred}\left(t_{new}\right)=f_{GPR}\left({\mathrm{\Δx}}_{d,j}\right(t_{new}))\\ y_{j,pred}\left(t_{new}\right)=y_j(t_{new}-j)+{\mathrm{\Δy}}_{j,pred}\left(t_{new}\right)\end{array}\right.---\left(10\right)$

步骤7：将最新的样本更新至原滑动窗口，剔除窗口中最旧的样本，进行下一次时滞估计，同时，等待新采样时刻的查询样本到来；不断重复上述步骤，可以提前实现塔底丁烷浓度的实时预测。

图6为基于不同TD模型进行丁烷浓度预测的可靠性对比结果，从图6可以看出，引入局部时滞的LTR-MWTDGPR模型相比于传统的MWTDGPR模型，以及传统的全局TDGPR模型在时间差间隔增加时，预测RMSE指标下降程度十分显著，意味着预测精度提高，对于变量漂移的处理更加有效。对于脱丁烷塔过程而言，该过程没有太大的突然波动，过程辅助变量和主导变量都有着很明显的变量漂移情况，而且过程伴随着较长时间的时滞，在不同的L取值情况下，考虑局部时滞和局部变量漂移的本发明方法都具有最低的预测RMSE，显著地提高了传统基于滑动窗口策略的TDGPR方法的可靠性。

接着，在L＝30,50,70,90,110时绘出了脱丁烷塔过程各方法的预测RMSE指标走势曲线，如图7所示。由图7可知，窗口大小为30时，预测性能最佳，可靠性最强，这是因为过程局部时滞充分被引入，校正的时序训练样本使得过程动态漂移的学习效果增强。在不同的窗口大小下，当j逐渐增大的时候，三种方法的RMSE指标变化呈现出不同的趋势，TDGPR方法建立的软测量模型性能下降急剧，而MWTDGPR方法在考虑时变特征和局部漂移后RMSE增加趋势放缓，而本发明方法随着j增加性能变差程度不明显，整体上较其它两种方法取得了令人满意的预测精度改善。窗口大小的选择对于脱丁烷塔过程软测量建模有着较大的影响，在图7里的5种窗口大小RMSE对比结果中，大致呈现了窗口大小越大，模型性能改善程度越低的状况。这充分地说明了脱丁烷塔过程，在窗口长度较大的时候，过程时滞变化越难以进行局部跟踪，而窗口较小时，能够充分考虑局部时滞特性，所得的MWTDGPR模型也更加精确，基于历史测量值的精度显著提升。经过一系列的仿真对比，进一步论证了考虑局部时滞的LTR-MWTDGPR方法对于塔底丁烷浓度实时预测的精度和可靠性，该发明方法即使基于很长时间前的主导变量值也相较于传统方法有着更佳的动态捕捉能力和预测精度。

以上所述仅是本发明的优选实施方式，并不用于限制本发明，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (2)

1.一种基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法，其特征在于，包括：

S1估计过程主导变量y与辅助变量集X＝[x₁,x₂,...,x_m]之间存在的最大时延参数T_max；

S2获取过程连续均匀采样输入输出样本集，建立包含L组连续时间样本的初始滑动窗口W_ini,W_ini＝[X(t),y(t)]_t＝1,...,L；

其中，X(t)＝[x₁(t),x₂(t),...,x_m(t)]_t＝1,...,L,m为辅助变量的个数，L为初始滑动窗口长度，为正整数；

S3对初始滑动窗口内m个辅助变量分别进行扩展，得到m(T_max+1)维的时滞扩展变量集，扩展方式为：

采用模糊曲线分析方法对时滞扩展变量集的各扩展变量进行分析，得到模糊覆盖范围最大的m个时滞变量，将各自的最优时滞参数记为d₁,d₂,...,d_m；

S4利用最优时滞参数对初始滑动窗口内L组实时对应的输入输出样本进行时滞校正，重构后的窗口数据个数变为L-T_max，重构窗口记为W_rec，其中

$W_{\mathrm{rec}}=[x_1(t-d_1),x_2(t-d_2),...,x_m(t-d_m),y\left(t\right){]}_{t=T_{\max }+1,...,L};$

S5当t_new采样时刻到来，无需采样辅助变量x(t_new)的情况，采用滑动窗口内的估计好的时滞参数对将到来的实时输入变量重构到d_i时刻之前作为查询样本进行建模，对x(t_new)进行时滞重构后得到的查询样本记为x_d(t_new)，这里的x_d(t_new)＝[x₁(t_new-d₁),x₂(t_new-d₂),...,x_m(t_new-d_m)]；

S6对重构窗口中重构的训练输入输出样本集及重构查询数据求取j次时间差分，然后对时间差分输入输出训练样本建立GPR模型，即建立了TDGPR模型结构，通过该模型结构得到主导变量的动态漂移预测值Δy_j,pred(t_new),在历史主导变量实时测量值基础上求出x_d(t_new)对应的实时输出y_j,pred(t_new)预测完成后，丢弃当前时滞重构的MWTDGPR模型，参见下式：

$\left\{\begin{array}{c}{\mathrm{\Δx}}_{d,j}\left(t_{new}\right)=x_d\left(t_{new}\right)-x_d(t_{new}-j)\\ {\mathrm{\Δy}}_{j,pred}\left(t_{new}\right)=f_{GPR}\left({\mathrm{\Δx}}_{d,j}\right(t_{new}\left)\right)\\ y_{j,pred}\left(t_{new}\right)=y_j(t_{new}-j)+{\mathrm{\Δy}}_{j,pred}\left(t_{new}\right)\end{array}\right.;$

S7将最新的样本更新至初始滑动窗口，剔除初始滑动窗口中最旧的样本，即窗口向前滑动一步，然后，重复上述步骤S1至S6。

2.根据权利要求1所述的基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法，其特征在于，采用模糊曲线分析方法得到滑动窗口内各辅助变量的最优时滞参数的具体方法为：

对于每个辅助变量扩展得到的窗口时滞输入变量集，记为{x_i(t-λ),λ＝0,1,…,T_max}，λ为引入的时延参数变量，假设初始滑动窗口内输入输出对应关系为x_i(t-λ)→y(t)，t＝T_max+1,...,L，那么，可将x_i(t-λ)的模糊隶属度函数定义为：

${\mathrm{\Φ}}_{it}\left(x_i\right(t-\mathrm{\λ}\left)\right)=\exp \[-{\left(\frac{I_n{x}_{it,\mathrm{\λ}}-x_i(t-\mathrm{\λ})}{b}\right)}^2\]$

Φ_it(x_i(t-λ))为变量x_i(t-λ)关于第t个数据点的输入变量模糊隶属度函数，b取变量x_i(t-λ)值域范围的20％，I_n＝[1,1,...,1]^T,n为x_i(t-λ)变量采样序列的长度，n＝L-T_max，x_it,λ表示x_i(t-λ)序列在t时刻的采样值，x_i(t-λ)为时延变量序列；

通过下式对扩展后的(T_max+1)个窗口新变量质心去模糊化，得第i个变量时延值为λ条件下的模糊曲线C_i,λ

$C_{i,\mathrm{\λ}}\left(\mathrm{\λ}\right)=\frac{{\Σ}_{t=T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]\·y\left(t\right)}{{\Σ}_{T_{\max }+1}^L{\mathrm{\Φ}}_{it}\[x_i(t-\mathrm{\λ})\]}$

其中，C_i,λ(λ)_max和C_i,λ(λ)_min为模糊曲线上点值域的最大值和最小值，覆盖范围即为模糊曲线最大值与最小值之差，d_i为使模糊曲线C_i,λ覆盖范围最大的λ，也即