CN112560924B

CN112560924B - 基于动态内部慢特征分析的丙烯精馏塔状态监控方法

Info

Publication number: CN112560924B
Application number: CN202011437212.8A
Authority: CN
Inventors: 杜文莉; 王振雷; 周钊; 田洲; 叶贞成; 钱锋
Original assignee: East China University of Science and Technology
Current assignee: East China University of Science and Technology
Priority date: 2020-12-07
Filing date: 2020-12-07
Publication date: 2024-01-19
Anticipated expiration: 2040-12-07
Also published as: CN112560924A

Abstract

本发明公开了一种基于动态内部慢特征分析算法(dynamic inner slow feature analysis，DiSFA)的丙烯精馏塔运行状态监控方法，所述方法包括：对丙烯精馏塔正常运行状态下的生产过程数据进行标准化处理；利用标准化后的数据拟合DiSFA模型，确定模型参数；利用PCA和ICA方法对DiSFA方法提取到的特征进行进一步分解；构造用于监控的统计量，确定对应的控制限；在线实时采集样本，利用得到的模型参数，计算对应该样本的统计量，并与相应的控制限进行比较，如统计量超出控制限则说明工作点发生了变化或离线建立的动态内模型不再适用。

Description

基于动态内部慢特征分析的丙烯精馏塔状态监控方法

技术领域

本发明属于化工过程运行状态监控技术领域，具体涉及一种丙烯精馏塔运行状态监控方法。

背景技术

精馏是一种广泛应用于化工、石油、轻工等生产过程中的传质工艺。它利用混合材料中各组分挥发性的差异，即相同温度下各组分的饱和蒸汽压差异，来分离其组分，高挥发性组分称为轻组分，低挥发性组分称为重组分。

精馏塔是精馏过程中的关键设备，它将原料分成不同的馏分，每个馏分的产品可以作为后续馏分的原料出售，也可以直接作为混合后的成品出售。生产过程如下：原料经加热炉加热，加热至所需温度后，原料板进入精馏塔，蒸馏后分离成不同的组分。在一段时间内，各组分的含量受多种因素的影响，如进料成分、温度、流量、各塔板的温度和压力、各侧线的产量、回流流量等。

丙烯精馏系统由2个大小相差很大的精馏塔串联组成。来自甲烷汽提塔塔釜的碳三馏份，进入丙烯精馏系统主塔(精馏段)的第76块塔板。经主塔高回流精馏后，主塔塔顶产出聚合级高纯度丙烯产品，同时将绝大多数丙烯产品回流到主塔内，进行精馏；主塔塔釜流出液经泵加压后进入副塔(提馏段)的塔顶进行再分离，同时副塔塔顶蒸汽直接进入主塔塔釜，这样两塔处于串联操作方式。

为了保证系统的长期自动化安全运行，实时监控系统的运行状态至关重要。这有利于及时有效发现故障，排除故障，避免造成恶劣的社会影响和重大经济损失。因此，精馏塔状态监控作为精馏研究的一个热点领域，具有重要的理论和实践意义。

精馏过程机理复杂，不适合基于机理模型的方法进行过程监控。另一方面，受控制系统的作用和设备的惯性和滞后特征，以及内外环境随机干扰因素的影响，系统中采集数据包含丰富的自相关和互相关关系，使得数据驱动的精馏塔状态监控研究具有一定挑战性。

发明内容

为了能够更有效的进行精馏过程监测，本发明考虑反馈控制作用对过程中采集数据的影响，通过同时对采集到数据间的动态关系和自相关关系进行建模，对丙烯精馏塔运行状态进行监控，这对于丙烯精馏塔长期自动化安全运行具有极其重要的实用价值。本发明解决了传统的多变量统计过程监控(MSPM)方法不能很好地拟合采集样本间动态特性的问题。本发明结合慢性先验，同时关注数据间的互相关性和自相关性，提取的特征更有利于过程故障检测。

本发明一方面提供一种基于动态内部慢特征分析算法的丙烯精馏塔状态监控模型的建模方法，所述方法包括：

(1)将丙烯精馏塔正常运行状态下的生产过程数据作为训练集；

(2)对步骤(1)得到的训练集数据进行标准化处理；

(3)利用步骤(2)标准化之后的数据拟合动态内部慢特征分析模型，确定模型参数，得到预提取数据特征，所述动态内部慢特征分析模型的目标函数如式(1)所示：

式中，X＝[x₁,x₂,...,x_N+s]^T表示标准化后的数据矩阵， X_i＝[x_i,x_i+1,...,x_N+i-1]^T(i＝1,2,...,s+1)，Z_s＝[X₁,X₂,...,X_s]，w为投影向量，β＝[β₁ β₂ … β_s]^T为不同自回归模型的拟合系数，为Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵；

(4)利用步骤(3)得到的特征进行进一步分解，并构建监控统计量，确定对应的监控限。

在一个或多个实施方案中，步骤(1)中，为保证采集到的样本数据能够反应生产过程的全部状态信息，样本数量不宜太少，采集的丙烯精馏塔生产过程数据的样本数一般不少于1000个。

在一个或多个实施方案中，步骤(2)中，数据标准化的目的是为了避免因原始数据中不同维度间数据量纲不同，导致对后续模型拟合产生的不利影响。标准化后的数据各维度数据均值为0，方差为1。

在一个或多个实施方案中，步骤(2)中，数据标准化为以训练集数据的均值μ和方差σ作为模型参数，并对训练集进行归一化处理。

在一个或多个实施方案中，步骤(2)中，所述特征包括动态潜变量和残差。

在一个或多个实施方案中，步骤(2)中，所述特征包括动态潜变量矩阵T、残差矩阵E、投影矩阵W、重构矩阵P和系数矩阵Θ^*。

在一个或多个实施方案中，步骤(3)中，所述基于慢特征分析算法拟合数据具体步骤如下：

初始化模型超参数，所述模型超参数包括动态阶次s和动态潜变量个数l，使用动态内部慢特征分析算法对步骤(2)标准化后的数据进行特征提取；

令X＝[x₁,x₂,...,x_N+s]^T表示标准化后的数据矩阵，定义 X_i＝[x_i,x_i+1,...,x_N+i-1]^T(i＝1,2,...,s+1)，Z_s＝[X₁,X₂,...,X_s]，则动态内部慢特征分析算法的目标函数定义如下：

式中，w为投影向量，β＝[β₁ β₂ … β_s]^T为不同自回归模型(AR模型)的拟合系数，为Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵。求解该目标函数，确定最优解w^*和β^*，即可确定一维动态潜特征t＝x^Tw^*；

为了继续提取下一维动态潜特征，需对数据进行如下处理：

X＝X-tp^T (2)

其中，

将根据公式(2)得到的X代入公式(1)中，迭代提取动态潜变量，直到提取潜变量的个数达到设定的超参数动态潜变量个数l；

设最终提取到的动态潜变量矩阵为T，T由l个动态潜变量组成，可以证明T中不同维的潜变量间相互正交，但它们在不同时滞间可能仍可能具有相关性；

参照根据X矩阵定义X_i矩阵的方式，根据T矩阵定义矩阵T_i(i＝1,2,...,s+1)，并定义建立如下目标函数，建模不同时滞间的相关性：

这是一个最小二乘问题，解析解为

计算残差其中P为l个p组成的矩阵。

步骤(3)中，所述慢特征分析算法的超参数需在拟合数据前确定，可依据经验进行确定，也可按照如下方法确定：

(a)初始化动态阶次s为1；

(b)设原始数据的维度为m，初始化为选择动态潜变量个数l为m-1；

(c)按照动态内部慢特征分析算法确定模型参数，所述模型参数包括残差和动态潜变量；

(d)计算残差矩阵E中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_ei(k)，k为相位差，对于所有维度的autocurr_ei(k)，若存在 autocurr_ei(k)>0.5，则l＝l+1；计算潜变量T中各维度的自相关系数，其中第i 维数据的自相关系数为autocurr_Ti(k)，k为相位差，对于所有维度的autocurr_Ti(k)，若存在autocurr_Ti(k)<0.2，则l＝l-1；

(e)重复步骤(c)-(d)，直到确定合适的参数l；

(f)在验证集上检验残差矩阵E是否包含明显的自相关关系，若是，修改动态阶次s＝s+1，所述验证集为丙烯精馏塔正常运行状态下的生产过程数据；

(g)重复步骤(c)-(f)，直到确定合适的参数s。

在一个或多个实施方案中，将丙烯精馏塔正常运行状态下的生产过程数据的一部分作为训练集，另一部分作为验证集。

在一个或多个实施方案中，步骤(4)中，进一步特征分解和统计量构造包含如下两个步骤。

第一步：进一步分解经动态内部慢特征分析处理后得到的特征；为避免T 中数据包含自相关性对分解过程的影响，与根据X矩阵定义X_i矩阵类似，根据T矩阵定义矩阵T_i(i＝1,2,...,s+1)，并定义则动态潜变量的一步超前预测值为/>可得预测残差为V＝T-T；为了充分挖掘其中蕴含的信息，按照式(4)-(7)，利用主成分分析(PCA)和独立成分分析(ICA) 方法分别从不同的角度进行分解V和E，分离各特征维度间的相关性：

E＝W_EICAE_{non_gaussICA}+E_noiseICA (6)

V＝W_VICAV_{non_gaussICA}+V_noiseICA (7)

其中，T_EPCA和E_noisePCA分别为由主成分分析算法得到的残差矩阵E对应的负载矩阵、得分矩阵和未建模残差矩阵；W_EICA、E_{non_gaussICA}和E_noiseICA分别为由独立成分分析算法得到的残差矩阵E对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；T_VPCA和V_noisePCA分别为由主成分分析算法得到的预测残差矩阵V对应的负载矩阵、得分矩阵和未建模残差矩阵，W_VICA、 V_{non_gaussICA}和V_noiseICA分别为由独立成分分析算法得到的预测残差矩阵V对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；

第二步：监控统计量的设计；

对于正常工况下的生产过程数据经归一化后的数据x_k，令R＝W(P^TW)^-1，其中，W为l个动态潜变量所对应的l个投影向量组成的矩阵，P为l个p组成的矩阵，按如下方式进行分解x_k：

为了充分利用得到的各部分特征，设计如下统计量：

ISPE_e＝(e_k-W_EICAe_{knon_gauss})(e_k-W_EICAe_{knon_gauss})^T (14)

ISPE_v＝(v_k-W_VICAv_{knon_gauss})(v_k-W_VICAv_{knon_gauss})^T (16)

其中，Λ_e为由残差矩阵E进行协方差分解得到的特征值对角阵，即n为E中样本个数；Λ_v为由残差矩阵V进行协方差分解得到的特征值对角阵，即n为V中样本个数；

按照核密度估计方法确定(9)-(17)中各个统计量的监控限，并标记所对应的监控限分别为SPE_{e_Lim}、/>SPE_{v_Lim}、/>ISPE_{e_Lim}、/>ISPE_{v_Lim}和/>

在一个或多个实施方案中，为了计算每个统计量的监控阈值，需要通过已知样本点估计对应统计量的概率分布密度，优选设置置信限99％处对应的值为监控阈值。

在一个或多个实施方案中，所述核密度估计方法利用公式进行数据拟合，其中z_j为第j个样本点，n位样本点的个数， h为窗口宽度，K(·)为核函数，优选为高斯密度函数。

在一个或多个实施方案中，由于各个统计量包含的物理意义不尽相同，将 (9)-(17)中的统计量分组融合，最终监控用统计量为：

对于融合后的统计量SIS和DIS，由其构造方式可确定两者的控制限均为 1。

本发明的另一个方面提供一种基于动态内部慢特征分析算法的丙烯精馏塔状态监控方法，所述方法包括离线建模和在线检测两个阶段，所述离线建模阶段包括本文任一实施方案所述的基于动态内部慢特征分析算法的丙烯精馏塔状态监控模型的建模方法的步骤(1)-(4)，所述在线检测阶段包括以下步骤：

(5)实时采集丙烯精馏塔的生产过程数据，记为x_new，利用训练集的均值μ和方差σ对x_new进行归一化处理；

(6)利用步骤(3)和(4)确定的模型参数，按式(8)对归一化后的x_new进行分解；

(7)利用步骤(6)的分解结果按式(9)-(17)计算基于实时数据的相应的统计量，并利用步骤(4)确定的各个统计量的监控限按式(18)、(19) 计算基于实时数据的融合后的统计量SIS_new和DIS_new：

对于融合后的统计量SIS和DIS，由其构造方式可确定两者的控制限均为1；

(8)判断SIS_new＞1、DIS_new＞1是否成立，若SIS_new＞1成立，则说明工作点发生了变化；若DIS_new＞1成立，则说明离线建立的动态内模型不再适用；若SIS_new＞1、 DIS_new＞1均不成立，则返回步骤(5)。

在一个或多个实施方案中，所述离线建模阶段包括：

(1’)采集正常工作状态下的数据作为历史数据集，并划分为训练集和验证集；

(2’)计算步骤(1’)训练集数据的均值μ和方差σ作为模型参数，并对训练集进行归一化处理；

(3’)初始化动态阶次s为1；

(4’)设原始数据的维度为m，初始化为选择动态潜变量个数l为m-1；

(5’)利用步骤(2’)归一化处理后的数据，按照动态内部慢特征分析算法确定模型参数；

(6’)计算残差矩阵E中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_ei(k)，k为相位差，对于所有维度的autocurr_ei(k)，若存在 autocurr_ei(k)>0.5，则l＝l+1；计算潜变量T中各维度的自相关系数，其中第i 维数据的自相关系数为autocurr_Ti(k)，k为相位差，对于所有维度的autocurr_Ti(k)，若存在autocurr_Ti(k)<0.2，则l＝l-1；

(7’)重复步骤(5’)-(6’)，直到确定合适的参数l；

(8’)在验证集上检验残差矩阵E是否包含明显的自相关关系，若是，修改动态阶次s＝s+1；

(9’)重复步骤(5’)-(8’)，直到确定合适的参数s；

(10’)利用步骤(5’)得到的模型处理步骤(2’)产生的归一化数据，得到残差矩阵E和预测残差矩阵V，并依PCA和ICA方法进行分解，确定参数W_EICA、/>W_VICA、Λ_e和Λ_v；

(11’)按照公式(9)-(17)计算相应的统计量，并基于核密度估计方法确定对应的监控限；

在一个或多个实施方案中，所述在线监测阶段包括：

(1”)实时采集到的丙烯精馏塔数据，设为_xnew，利用离线阶段步骤(2’) 确定的参数μ和σ对_xnew进行归一化处理；

(2”)利用离线阶段步骤(5’)确定的模型参数，计算对应的量 t_new,e_new,v_new；

(3”)利用离线阶段步骤(11’)确定的各统计量控制限，按公式(9)-(17) 计算相应的统计量，并按公式(18)、(19)计算融合后的统计量；

(4”)判断SIS_new＞1，DIS_new＞1是否成立，若SIS_new＞1成立，则说明工作点发生了变化；若DIS_new＞1成立，则说明离线建立的动态内模型不再适用，可能是因为发生了故障，导致系统控制质量下降，也可能是因为系统转移到了新的工作点，但该工作点与训练用历史数据所在工作点存在不一样的动态特性，需进一步判断处理，否则，返回步骤(1”)。

附图说明

图1为本发明的基于动态内部慢特征分析算法的丙烯精馏塔状态监控方法的处理流程示意图。

图2为实施例1利用本发明的基于动态内部慢特征分析算法的丙烯精馏塔状态监控方法得到的丙烯精馏塔装置状态监控结果，其中，上图展示了各采样点的SIS统计量，下图展示了各采样点的DIS统计量。

图3为对比例1采用经典PCA方法得到的丙烯精馏塔装置状态监控结果，其中，上图展示了各采样点的T²统计量，下图展示了各采样点的SPE统计量。

具体实施方式

为使本领域技术人员可了解本发明的特点及效果，以下谨就说明书及权利要求书中提及的术语及用语进行一般性的说明及定义。除非另有指明，否则文中使用的所有技术及科学上的字词，均为本领域技术人员对于本发明所了解的通常意义，当有冲突情形时，应以本说明书的定义为准。

本文中，为使描述简洁，未对各个实施方案或实施例中的各个技术特征的所有可能的组合都进行描述。因此，只要这些技术特征的组合不存在矛盾，各个实施方案或实施例中的各个技术特征可以进行任意的组合，所有可能的组合都应当认为是本说明书记载的范围。

本发明提出一种基于动态内部慢特征分析算法(dynamic inner slow featureanalysis，DiSFA)的故障检测方法。本发明提出的DiSFA算法通过内部自回归(auto-regression,简称AR)模型和外部协方差最大化模型进行建模，以保证可以同时提取数据间的自相关和互相关关系，将动态特性和工作点信息分开监控，有利于提供当前工况更准确的信息，将实际故障与操作条件变化区分开来。本发明的基于DiSFA的故障检测方法可应用到丙烯精馏实际过程中，实现对动态过程故障的有效检测。本发明的基于DiSFA的丙烯精馏塔状态监控方法的处理流程示意图如图1所示。

本发明中，整个方法实施流程包括离线建模和在线监测两个阶段。离线建模部分利用正常工况下从丙烯精馏系统中采集的数据作为训练数据，进行标准化处理，在此基础上拟合DiSFA模型，提取能够表征过程工作状态的特征，并利用主成分分析(PCA)和独立成分分析(ICA)方法从不同的角度对特征进行再分解，构建相应的统计量，利用核密度估计方法确定控制限，最后按统计量的物理意义分组融合，用于进行故障检测。在线检测阶段，实时从丙烯精馏塔过程中采集数据，利用离线阶段获得的正常工况下各维度数据的均值和方差对实时数据进行标准化处理，而后利用离线确定的DiSFA、PCA和ICA模型参数，对特征进行分解，最后根据分解得到的特征计算SIS和DIS统计量，若超出控制限，则可以认为当前监控数据反应精馏塔过程内可能发生了故障。

本发明所述的DiSFA算法是指一种利用标准化后数据进行特征提取的算法，具体是指公式(1)所述的目标函数及其求解过程：

式中，X＝[x₁,x₂,...,x_N+s]^T表示标准化后的数据矩阵， X_i＝[x_i,x_i+1,...,x_N+i-1]^T(i＝1,2,...,s+1)，Z_s＝[X₁,X₂,...,X_s]，w为投影向量，β＝[β₁ β₂ … β_s]^T为不同自回归模型的拟合系数，为Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵。

本发明所述的DiSFA模型是指利用DiSFA算法建立的模型，包括式(1) 所示的目标函数、以及该目标函数求解过程中涉及到的公式和函数。

本发明中，使用DiSFA算法建立特征预提取模型，基于PCA、ICA方法进行特征再分解，以及相关统计量的计算统称为基于DiSFA的过程监控方法。本发明中，提取的特征包括动态潜变量矩阵T、残差矩阵E、投影矩阵W、重构矩阵P和系数矩阵Θ^*中的一个或多个或优选全部。

本发明将现场采集到的丙烯精馏塔在正常工况下运行时产生的生产过程数据作为训练数据。本文中，正常数据具有本领域常规的含义，是指系统运行平稳，无异常现象发生的工况。本文中，所使用的丙烯精馏塔生产过程数据包含多个维度，同一时刻采集得到的数据组成一个向量，称为一个样本。数据集以矩阵形式呈现，矩阵中每行数据表示一个样本，每列数据(列向量)对应同一个物理量在不同时刻的值。矩阵中的行，即每个样本的采集时刻，应按照时间顺序依次排列。本发明中，为充分描述正常工况包含的特征，用于训练模型的数据不宜太少，一般不少于1000个。若样本的采样间隔为1min，则训练集样本数优选为1000-2000组。

本发明利用实时采集的丙烯精馏塔过程数据(简称在线数据、实时数据) 作为测试数据验证模型的有效性。作为测试数据的在线数据可能包含异常工况下的数据。在线数据的采集方式应保持与建立模型时所用离线数据的采集方式保持相同，包括数据采样周期，各个维度数据的数据源。在线进行过程监测时，每次针对单个数据进行计算，判断其对应的监控统计量是否超出了监控范围，若超出监控范围，则模型指示当前数据点受到了非正常因素的影响，即丙烯精馏塔过程可能发生了故障。通过比较本发明所给方法的判断结果与实际运行状态是否一致，可以确认是否出现检测错误，进而计算检测率和误报率。

本发明中，用作训练数据、测试数据的生产过程数据的测量变量优选包括进料量、进料组分1、进料组分2、塔压差、中沸器流量、中沸器液位、裂解气入口温度、塔顶甲烷浓度、采出量、不凝气量、冷凝器1液位、冷凝器2液位、丙烯冷剂流量、丙烯入口温度、凝液罐液位和塔釜采出量。

本发明将离线建模阶段计算的训练数据集的均值μ和方差σ作为参数对在线数据进行标准化处理。本发明中，基于平均值和方差对数据进行标准化处理的方法是本领域常规的。训练数据集在标准化后各个维度的数据均满足均值为0、方差为1的关系。数据标准化的目的是为了避免因原始数据中不同维度间数据量纲不同，导致对后续模型拟合产生的不利影响。对训练数据集进行标准化处理可以是，设原始数据矩阵中，第i维的数据表示为则对于第i维数据的标准化公式为：/>其中，/> 为标准化后的第i维数据。

本发明利用标准化后的数据拟合DiSFA模型，确定模型参数(超参数)，进行特征提取。DiSFA模型中存在两个重要的超参数：模型动态阶次s和提取潜变量(又称动态潜特征、动态潜变量)的个数l，在s和l给定的情况下，对标准化后的数据矩阵X＝[x₁,x₂,...,x_N+s]^T，定义 X_i＝[x_i,x_i+1,...,x_N+i-1]^T(i＝1,2,...,s+1)，Z_s＝[X₁,X₂,...,X_s]，则DiSFA方法的目标函数是：

式中，w为投影向量，β＝[β₁ β₂ … β_s]^T为不同自回归(autoregressive，简称 AR)模型的拟合系数，为Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵。求解该目标函数，即确定最优解w^*和β^*。由w^*确定一维动态潜特征t＝x^Tw^*。为了继续提取下一维动态潜特征，需对数据进行如下处理：

X＝X-tp^T (2)

其中，迭代提取动态潜变量，直到提取潜变量的个数达到设定超参数动态潜变量个数l。设最终提取到的动态潜变量矩阵为T，可以证明T中不同维的潜变量间相互正交，但它们在不同时滞间可能仍可能具有相关性。采用根据X矩阵定义X_i矩阵类似的方法，根据T矩阵定义矩阵T_i(i＝1,2,...,s+1)，并定义/>建立如下目标函数，建模不同时滞间的相关性：

这是一个最小二乘问题，解析解为

计算残差本文中，P是指l个前述p组成的矩阵。

为描述方便，本发明关于求解DiSFA模型参数的计算过程伪代码给出如下：

输入：数据矩阵X，提取潜变量个数l，动态阶数s。

输出：动态潜变量矩阵T，残差矩阵E，投影矩阵W，重构矩阵P，系数矩阵Θ^*。

归一化原始数据矩阵X，使处理后数据矩阵各变量满足零均值和单位方差。

for i＝1:l

随机初始化一个单位向量w

do{

β_i＝G^-1([t₁,t₂,...,t_s]^Tt_s+1)

while(w与β均收敛)

t＝Xw

X＝X-tp^T

W＝[W,w_i]

End

构造矩阵T，T_i(i＝1,2,...,s+1)，计算/>

计算残差：

本发明中，超参数s和l可由操作人员依据经验进行确定，也可按照如下方法确定：

(a)初始化动态阶次s为1；

(c)按照DiSFA算法确定模型的动态潜变量和残差；

(d)计算训练集中得到的残差矩阵E的自相关系数，若存在自相关系数大于预设阈值，则l＝l+1；计算训练集中得到的动态潜变量T的自相关系数，若存在噪声序列，则l＝l-1；

(e)重复步骤(c)-(d)，直到确定合适的参数l；

(f)在验证集上检验残差矩阵E是否包含明显的自相关关系，若是，修改动态阶次s＝s+1；

(g)重复步骤(c)-(f)，直到确定合适的参数s。

本文中，自相关可以理解为同一个信号在不同时间点的互相关，可以通过计算时间序列与其自身在不同时间点的互相关系数进行评估，一般若计算得到的相关系数均处于[0,0.1]范围内，可以认为不存在相关性；处于[0.1,0.3]，可以认为存在弱相关关系；处于[0.3,0.5]，可以认为存在中等相关关系；处于[0.5,1]，可以认为存在强相关关系。对于时间序列来讲，噪声的特征在于其没有明显的自相关性，同样可以通过计算时间序列与其自身在不同时间点的互相关系数进行评估。本发明中，前述自相关系数的预设阈值可以是0.5。本发明中，当自相关系数小于0.2时可认为存在噪声序列。

在一些实施方案中，步骤(d)包括：计算残差矩阵E中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_ei(k)，k为相位差，对于所有维度的autocurr_ei(k)，若存在autocurr_ei(k)>0.5，则l＝l+1；计算潜变量T中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_Ti(k)，k为相位差，对于所有维度的autocurr_Ti(k)，若存在autocurr_Ti(k)<0.2，则l＝l-1。本发明中，相位差k一般选为1。

残差矩阵E中第i维数据的自相关系数autocurr_ei(k)的计算公式为：

其中，k为相位差，一般选为1；为残差矩阵E中第i维数据中的第t个数据；/>为残差矩阵E中第i维数据中的第t-k个数据；/>为残差矩阵E中第i维数据的平均值。

潜变量T中第i维数据的自相关系数autocurr_Ti(k)的计算公式为：

其中，k为相位差，一般选为1；为潜变量T中第i维数据中的第t个数据；/>为潜变量T中第i维数据中的第t-k个数据；/>为潜变量T中第i维数据的平均值。

本发明中，验证集(验证数据)为现场采集到的丙烯精馏塔在正常工况下运行时产生的生产过程数据，且与训练集不重复。在一些实施方案中，采集丙烯精馏塔正常运行状态下的生产过程数据，并划分为训练集和验证集。

本发明中，经DiSFA方法初步提取特征后，需进一步进行特征分解和统计量构造两个步骤。

第一步：进一步分解经动态内部慢特征分析处理后得到的特征。

为避免T中数据包含自相关性对分解过程的影响，采用根据X矩阵定义X_i矩阵类似的方法，根据T矩阵定义矩阵T_i(i＝1,2,...,s+1)，并定义则动态潜变量的一步超前预测值为/>可得预测残差为V＝T-T。

本发明中，进一步分解经动态内部慢特征分析处理后得到的特征，是指进一步分解残差矩阵E和预测残差矩阵V。

为了充分挖掘其中蕴含的信息，利用PCA和ICA方法按照公式(4)-(7) 分别从不同的角度进行分解V和E，分离各特征维度间的相关性：

E＝W_EICAE_{non_gaussICA}+E_noiseICA (6)

V＝W_VICAV_{non_gaussICA}+V_noiseICA (7)

其中，T_EPCA和E_noisePCA分别为由主成分分析算法得到的残差矩阵E对应的负载矩阵、得分矩阵和未建模残差矩阵；W_EICA、E_{non_gaussICA}和E_noiseICA分别为由独立成分分析算法得到的残差矩阵E对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；T_VPCA和V_noisePCA分别为由主成分分析算法得到的预测残差矩阵V对应的负载矩阵、得分矩阵和未建模残差矩阵，W_VICA、 V_{non_gaussICA}和V_noiseICA分别为由独立成分分析算法得到的预测残差矩阵V对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵。利用PCA算法计算负载矩阵、得分矩阵和未建模残差矩阵以及利用ICA算法计算独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵的方法是已知的。

进一步分解经动态内部慢特征分析处理后得到的特征的目的在于求取公式(4)-(7)中的W_EICA、/>和W_VICA。

第二步：监控统计量的设计。

对于训练数据集经归一化处理后形成的数据x_k，令R＝W(P^TW)^-1，其中， W为l个动态潜变量所对应的l个投影向量组成的矩阵，按如下方式进行分解 x_k：

为了充分利用得到的各部分特征，设计如下统计量：

ISPE_e＝(e_k-W_EICAe_{knon_gauss})(e_k-W_EICAe_{knon_gauss})^T (14)

ISPE_v＝(v_k-W_VICAv_{knon_gauss})(v_k-W_VICAv_{knon_gauss})^T (16)

本发明中，式(9)-(17)中各个统计量的控制限是基于训练数据集利用核密度估计方法求得的，并标记所对应的监控限(监控阈值)分别为SPE_{e_Lim}、/>SPE_{v_Lim}、ISPE_{e_Lim}、/>ISPE_{v_Lim}和/>核密度估计是求取监控限的成熟方法，该方法利用公式/>进行数据拟合，其中z_j为第j个样本点，n位样本点的个数，h为窗口宽度，K(·)为核函数，通常为高斯密度函数。为了计算每个统计量的监控阈值，需要通过已知样本点估计对应统计量的概率分布密度。本发明优选设置置信限为99％处对应的值为监控阈值。

由于各个统计量包含的物理意义不尽相同，将式(9)-(17)中的统计量分组融合，最终本发明的监控用统计量为：

对于融合后的统计量SIS和DIS，由其构造方式可确定两者的控制限均为1。

本发明中，利用训练数据得到模型的所有参数后，可利用得到的这些模型参数(例如包括动态潜变量矩阵T、残差矩阵E、投影矩阵W、重构矩阵P、系数矩阵Θ^*、W_EICA、和W_VICA)进行在线过程监测。在线过程监测时，实时采集丙烯精馏塔过程数据，利用训练集的均值μ和方差σ进行标准化处理，利用基于训练数据得到的模型参数根据式(8)对标准化处理后的在线数据进行分解，按照公式(9)-(17)计算基于在线数据的各个统计量，再利用基于训练数据得到的各个统计量的监控限按照公式(18)-(19)计算基于在线数据的融合后的SIS和DIS统计量，并判断SIS和DIS是否超出控制限1，如果超出则记为异常状态，由此实现对丙烯精馏塔过程状态的监控。

本发明中，异常状态可能是系统发生故障(未正常运转)，也可能是系统工作点发生转移(变化)。具体而言，若_SISnew＞1成立，则说明工作点发生了变化；若_DISnew＞1成立，则说明离线建立的动态内模型不再适用。本文中，“动态内模型不适用”是指根据训练数据得到的“动态内模型”无法描述当前时刻的数据变化趋势，因此导致DIS统计量发生显著变化，并非指DiSFA方法对于监控当前丙烯精馏塔过程状态不适用。

本文中，SIS和DIS统计量的值所反应的丙烯精馏塔过程状态如表1所示。

表1：SIS和DIS统计量意义

	DIS_new＜1	DIS_new＞1
			SIS_new＜1	正常工况	工作点未变，但系统运行不平稳
SIS_new＞1	工作点变化，但系统仍正常运行	工作点变化，且系统运行不平稳

本文中，静态特性是指系统正常运行时过程变量间的互相关关系，关注的是同一时刻不同变量间的互相关性。工作点发生变化是指在线数据所在的工作点与训练用历史数据所在工作点存在不一样的静态特性。动态特性指的是包括控制器在内的整个闭环系统内，采集到的相关物理量中表现出的变化趋势，动态特性的主要关注点在于时间序列的自相关特性。SIS统计量是由残差矩阵E 通过分解后计算得到的，E表示的潜变量表示原始数据的能力，如果E发生了显著变化，说明得到的潜变量不能作为采集到的数据的低维表示。因此，SIS 统计量描述建模阶段使用数据的工作点信息，因操作变量变化或发生故障都会使得工作点变化。DIS统计量是由矩阵V通过分解后计算得到的，V表示的是由原始数据直接投影得到的潜变量和由过去时刻的原始数据投影得到的潜变量预测的当前时刻潜变量的残差，即表示的是由过去预测当前的能力，V发生了显著变化则说明变化趋势发生了变化，不能由离线建模时得到的自回归模型描述，即动态特性发生了变化。DIS统计量监控的是序列间的自相关关系，DIS 统计量未超出控制限说明各潜变量间的自相关关系与建模阶段相似，系统仍在正常运转，若DIS统计量超出控制限，可能是故障导致，也可能是因为转移到的新工作点与旧工作点之间的动态性差别较大的原因导致。

本发明的主要优点在于：

(1)本发明能够更好的处理具有动态特征的过程数据；

(2)本发明能够同时建模数据间的自相关和互相关关系；

(3)本发明提出的监控统计量融合了从不同角度得到特征的表达能力，有利于提高故障检测率。

下面通过实施例对本发明进行具体描述。有必要在此指出的是，以下实施例只用于对本发明作进一步说明，不能理解为对本发明保护范围的限制，本领域技术人员根据本发明的内容作出的一些非本质的改进和调整，仍属于本发明的保护范围。

实施例1

下面以本发明在某一石化企业具体实施的情况，给出采用本发明的基于动态内部慢特征分析(DiSFA)方法的丙烯精馏塔状态监控方法对丙烯精馏塔进行状态监控的详细操作流程。

1、将现场采集到的丙烯精馏塔正常工况下1000组样本数据作为训练数据，对数据集中的每维变量进行标准化处理，设原始数据矩阵中，第i维的数据表示为则对于第i维数据的标准化公式为：/>其中，/> 为标准化后的第 i维数据。本实施例中选取的测量变量包括进料量、进料组分1、进料组分2、塔压差、中沸器流量、中沸器液位、裂解气入口温度、塔顶甲烷浓度、采出量、不凝气量、冷凝器1液位、冷凝器2液位、丙烯冷剂流量、丙烯入口温度、凝液罐液位、塔釜采出量共计16个变量；

2、利用标准化后的数据拟合DiSFA模型，确定模型参数，DiSFA模型存在两个超参数，即提取潜变量个数l和动态阶数s，根据文本所述的确定超参数s和l的步骤(a)到(g)，通过选择多组不同的参数进行拟合，对比相应的DiSFA模型提取到的动态潜变量的预测残差和剩余残差部分的自相关性，从而确定合适的超参数，本实施例中，最终选择参数潜变量个数l＝5，动态阶数 s＝3。

DiSFA方法的目标函数是：

式中，w为投影向量，β＝[β₁ β₂ … β_s]^T为不同AR模型的拟合系数，为 Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵。求解该目标函数，确定最优解w^*和β^*，进而确定一维动态潜特征t＝x^Tw^*。

为了继续提取下一维动态潜特征，需对数据进行如下处理：

X＝X-tp^T (2)

其中，

迭代提取动态潜变量，直到提取潜变量的个数达到设定超参数动态潜变量个数l。设最终提取到的动态潜变量矩阵为T，可以证明T中不同维的潜变量间相互正交，但它们在不同时滞间可能仍具有相关性。与根据X矩阵定义X_i矩阵类似，根据T矩阵定义矩阵T_i(i＝1,2,...,s+1)，并定义建立如下目标函数，建模不同时滞间的相关性：

这是一个最小二乘问题，解析解为

残差为

求解DiSFA模型参数的计算过程详细如下：

输入：数据矩阵X，提取潜变量个数l，动态阶数s。

for i＝1:l

随机初始化一个单位向量w

do{

β_i＝G^-1([t₁,t₂,...,t_s]^Tt_s+1)

while(w与β均收敛)

t＝Xw

X＝X-tp^T

W＝[W,w_i]

End

构造矩阵T，T_i(i＝1,2,...,s+1)，计算/>

计算残差：

3、经DiSFA方法初步提取特征后，进一步进行特征分解和统计量构造两个步骤。

为避免T中数据包含自相关性对分解过程的影响，与根据X矩阵定义X_i矩阵类似，根据T矩阵定义矩阵T_i(i＝1,2,...,s+1)，并定义则动态潜变量的一步超前预测值为/>可得预测残差为V＝T-T。

为了充分挖掘其中蕴含的信息，利用PCA和ICA方法分别从不同的角度进行分解V和E，分离各特征维度间的相关性，分解方式如式(4)-(7)所示：

E＝W_EICAE_{non_gaussICA}+E_noiseICA (6)

V＝W_VICAV_{non_gaussICA}+V_noiseICA (7)

其中，T_EPCA和E_noisePCA分别为由主成分分析算法得到的残差矩阵E对应的负载矩阵、得分矩阵和未建模残差矩阵；W_EICA、E_{non_gaussICA}和E_noiseICA分别为由独立成分分析算法得到的残差矩阵E对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；T_VPCA和V_noisePCA分别为由主成分分析算法得到的预测残差矩阵V对应的负载矩阵、得分矩阵和未建模残差矩阵，W_VICA、 V_{non_gaussICA}和V_noiseICA分别为由独立成分分析算法得到的预测残差矩阵V对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵。

第二步：监控统计量的设计。

对于归一化后的数据x_k，令R＝W(P^TW)^-1，其中，W为l个动态潜变量所对应的l个投影向量组成的矩阵，按如下方式进行分解x_k：

为了充分利用得到的各部分特征，设计如下统计量：

ISPE_e＝(e_k-W_EICAe_{knon_gauss})(e_k-W_EICAe_{knon_gauss})^T (14)

ISPE_v＝(v_k-W_VICAv_{knon_gauss})(v_k-W_VICAv_{knon_gauss})^T (16)

其中，Λ_e为由残差矩阵E进行协方差分解得到的特征值对角阵，即n为E中样本个数；Λ_v为由残差矩阵V进行协方差分解得到的特征值对角阵，即/>n为V中样本个数。

本实施例中，式(9)-(17)中各个统计量的控制限是基于训练数据集利用核密度估计方法求得的，核密度估计方法利用公式进行数据拟合，其中z_j为第j个样本点，n为样本点的个数，h为窗口宽度，K(·) 为高斯密度函数，将置信限为99％处对应的值为监控限。将各个统计量所对应的监控限分别记为/>SPE_{e_Lim}、/>SPE_{v_Lim}、ISPE_{e_Lim}、/>ISPE_{v_Lim}和/>由于各个统计量包含的物理意义不尽相同，将式(9)-(17) 中的统计量分组融合，最终监控用统计量为：

本发明中，得到模型的所有参数后，可利用得到的这些模型参数进行在线过程监测。在线过程监测时，实时采集丙烯精馏塔过程数据，并进行标准化处理，先后经DiSFA、PCA和ICA算法进行特征提取，计算基于在线数据的SIS 和DIS统计量，并判断是否超出控制限1，如果SIS或DIS统计量超出1则记为异常状态，由此实现对丙烯精馏塔过程状态的监控。

5、丙烯精馏塔生产过程状态在线监控

训练完成后，利用在线采集的丙烯精馏塔的生产过程数据作为测试数据进行验证，即参照上述步骤1对数据进行标准化处理，参照上述步骤2确定DiSFA 模型参数，参照上述步骤3进一步分解经DiSFA方法提取到的特征，并计算相应的SIS和DIS统计量，判断是否超出控制限1，如果超出则记为异常状态。本实施例共验证了1000组在线数据，监控结果如图2所示，图中与x轴平行的虚线表示统计量的控制限，与y轴平行的虚线表示故障数据与正常数据的分割线，虚线右边即从第300个数据开始为进料量发生变化，导致工作点变化，但系统仍在运转。

本发明的监控方法提供的SIS统计量反应同一时刻采集到变量间的互相关关系，可以反映出当前工作点是否改变，而DIS统计量反应根据历史潜变量预测当前潜变量的能力，由于用于预测的自回归模型的模型参数是根据离线数据训练得到的，如果DIS统计量发生显著变化，则说明此自回归模型不能正确描述当前数据的变化趋势(即文中所述动态特性)，也即系统的平稳程度与正常工况不同，反应系统内发生了某些故障，导致平稳性发生变化。从图2中可以看出，SIS统计量能够及时检测到这一工作点变化的工况，但DIS统计量只有一段时间超出控制限，这是因为实施例1中引入的变化是由于进料量改变，并非真正意义的故障，因此，DIS统计量并未发生显著变化，系统工作正常。相比于传统方法，本方法将工作点信息和动态特性分开监控，在实际生产过程中，有利于提供当前工况更准确的信息，帮助从业者根据监测结果将实际故障与操作条件变化这种正常工况区分开来。

对比例1

利用与实施例1相同的训练数据采用经典的主成分分析(PCA)方法建立监控模型，并使用与实施例1相同的测试数据进行验证。

图3为经典PCA方法得到的监控图，图中与x轴平行的虚线表示统计量的控制限，与y轴平行的虚线表示故障数据与正常数据的分割线，从图中可以看出，两个统计量在进料量发生变化后，均超出控制限，即PCA方法只能监控静态工作点信息，而不能提供动态信息的表征。

Claims

1.一种基于动态内部慢特征分析算法的丙烯精馏塔状态监控模型的建模方法，其特征在于，所述方法包括：

(2)对步骤(1)得到的训练集数据进行标准化处理；

(3)利用步骤(2)标准化之后的数据拟合动态内部慢特征分析模型，确定模型参数，提取数据特征，所述特征包括动态潜变量和残差，所述动态内部慢特征分析模型的目标函数如式(1)所示：

式中，X＝[x₁，x₂，...，x_N+s]^T表示标准化后的数据矩阵，X_i＝[x_i，x_i+1，...，x_N+i-1]^T(i＝1，2，...：s+1)，Z_s＝[X₁，X₂，...，X_s]，w为投影向量，β＝[β₁ β_２ … β_s]^T为不同自回归模型的拟合系数，为Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵；

(4)采用主成分分析法和独立成分分析法对步骤(3)得到的特征进行进一步分解，并构建监控统计量，确定对应的监控限；

步骤(3)中，基于慢特征分析算法拟合数据的步骤如下：

初始化模型超参数，所述模型超参数包括动态阶次s和动态潜变量个数l，使用动态内部慢特征分析算法对步骤(2)标准化后的数据进行特征提取；令X＝[x₁，x₂，...，x_N+s]^T表示标准化后的数据矩阵，定义X_i＝[x_i，x_i+1，...，x_N+i-1]^T(i＝1，2，...，s+1)，Z_s＝[X₁，X₂，...，X_s]，则动态内部慢特征分析算法的目标函数定义如下：

式中，w为投影向量，β＝[β₁ β₂ … β_s]^T为不同自回归模型的拟合系数，为Kronecker积符号，ΔZ_s为Z_s的一阶差分矩阵，求解该目标函数，确定最优解w^*和β^*，即可确定一维动态潜特征t＝x^Tw^*；

为了继续提取下一维动态潜特征，对数据进行如下处理：

X＝X-tp^T (2)

其中，

迭代提取动态潜变量，直到提取潜变量的个数达到设定的超参数动态潜变量个数l；

设最终提取到的动态潜变量矩阵为T，与根据X矩阵定义X_i矩阵类似，根据T矩阵定义矩阵T_i(i＝1，2，...，s+1)，并定义建立如下目标函数，建模不同时滞间的相关性：

这是一个最小二乘问题，解析解为

计算残差其中P为l个p组成的矩阵；

步骤(3)中，所述慢特征分析算法的超参数s和l在拟合数据前依据经验进行确定，或按照如下方法确定：

(a)初始化动态阶次s为1；

(c)按照动态内部慢特征分析算法确定动态潜变量和残差；

(d)计算残差矩阵E中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_ei(k)，k为相位差，对于所有维度的autocurr_ei(k)，若存在autocurr_ei(k)＞0.5，则l＝l+1；计算潜变量T中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_Ti(k)，k为相位差，对于所有维度的autocurr_Ti(k)，若存在autocutr_Ti(k)＜0.2，则l＝l-1；

(e)重复步骤(c)-(d)，直到确定合适的参数l；

(f)以丙烯精馏塔正常运行状态下的生产过程数据作为验证集，在验证集上检验残差矩阵E是否包含明显的自相关关系，若是，修改动态阶次s＝s+1；

(g)重复步骤(c)-(f)，直到确定合适的参数s；

步骤(4)中，进一步特征分解和统计量构造包含如下两个步骤：

第一步：进一步分解经动态内部慢特征分析处理后得到的特征；

动态潜变量的一步超前预测值为可得预测残差为/>利用主成分分析和独立成分分析方法分解预测残差V和残差E，分离各特征维度间的相关性，分解结果如式(4)-(7)所示：

E＝W_EICAE_{non_gaussIC}A+E_noiseICA (6)

V＝W_VICAV_{nan_gaussICA}+V_noiseICA (7)

其中，T_EPCA和E_noisePCA分别为由主成分分析算法得到的残差矩阵E对应的负载矩阵、得分矩阵和未建模残差矩阵；W_EICA、E_{man_gaussICA}和E_noiseICA分别为由独立成分分析算法得到的残差矩阵E对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；/>T_VPCA和V_noisePCA分别为由主成分分析算法得到的预测残差矩阵V对应的负载矩阵、得分矩阵和未建模残差矩阵，W_VIC，A、V_{nan_gaussICA}和V_noiseICA分别为由独立成分分析算法得到的预测残差矩阵V对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；

第二步：监控统计量的设计；

设计如下统计量：

ISPE_e＝(e_k-W_EICAe_{knon_gauss})(e_k-W_EICAe_{non_gauss}) (14)

ISPE_v＝(v_k-W_VICAv_{know_gauss})(v_k-W_VICAv_{know_gauss})^T (16)

采用核密度估计方法确定式(9)-(17)中各个统计量的监控限，并将各个统计量所对应的监控限分别记为SPE_{e_Lim}、/>SPE_{v_Lim}、/>ISPE_{e_Lim}、/>ISPE_{v_Lim}和

2.如权利要求1所述的基于动态内部慢特征分析算法的丙烯精馏塔状态监控模型的建模方法，其特征在于，步骤(1)中，采集的丙烯精馏塔的生产过程数据的样本数不少于1000组。

3.如权利要求1所述的基于动态内部慢特征分析算法的丙烯精馏塔状态监控模型的建模方法，其特征在于，步骤(2)中，标准化后的数据各维度数据均值为0，方差为1。

4.一种基于动态内部慢特征分析算法的丙烯精馏塔状态监控方法，其特征在于，所述方法包括：

(2)对步骤(1)得到的训练集数据进行标准化处理；

(3)利用步骤(2)标准化之后的数据拟合动态内部慢特征分析模型，确定模型参数，提取数据特征，所述特征包括动态潜变量和残差，步骤如下：

令X＝[x₁，x₂，...，x_N+s]^T表示标准化后的数据矩阵，定义X_i＝[x_i，x_i+1...，x_N+i-1]^T(i＝1，２，...，s+1)，Z_s＝[X₁，X₂，...，X_s]，则动态内部慢特征分析算法的目标函数定义如下：

为了继续提取下一维动态潜特征，对数据进行如下处理：

X＝X-tp^T (2)

其中，其中P为l个p组成的矩阵；

这是一个最小二乘问题，解析解为

计算残差

(4)采用主成分分析法和独立成分分析法对步骤(3)得到的特征进行进一步分解，并构建监控统计量，确定对应的监控限，包含如下两个步骤：

E＝W_EICAE_{non_gaussICA}+E_noiseICA (6)

V＝W_VICAV_{non_gaussICA}+V_noiseICA (7)

其中，T_EPCA和E_noisePCA分别为由主成分分析算法得到的残差矩阵E对应的负载矩阵、得分矩阵和未建模残差矩阵；W_EICA、E_{non_gaussICA}和E_noiseICA分别为由独立成分分析算法得到的残差矩阵E对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；/>T_VPCA和V_noisePCA分别为由主成分分析算法得到的预测残差矩阵V对应的负载矩阵、得分矩阵和未建模残差矩阵，W_VICA、V_{non_gaussICA}和V_noiseICA分别为由独立成分分析算法得到的预测残差矩阵V对应的独立因子混合矩阵、非高斯成分矩阵和剩余高斯成分矩阵；

第二步：监控统计量的设计；

设计如下统计量：

ISPE_e＝(e_k-W_EICAe_{know_gauss})(e_k-W_EICAe_{know_gauss})^T (14)

ISPE_v＝(v_k-W_VICAv_{knon_gauss})(v_k-W_VICAv_{know_gauss}) (16)

其中，Λ_e为由残差矩阵E进行协方差分解得到的特征值对角阵，即n为E中样本个数；Λ_v为由残差矩阵V进行协方差分解得到的特征值对角阵，即/>n为V中样本个数；

采用核密度估计方法确定式(9)-(17)中各个统计量的监控限，并将各个统计量所对应的监控限分别记为SPE_eLim、/>SPE_{v_Lim}、/>ISPE_{e_Lim}、/>ISPE_{v_Lim}和

(7)利用步骤(6)的分解结果按式(9)-(17)计算基于实时数据的相应的统计量，并利用步骤(4)确定的各个统计量的监控限按式(18)、(19)计算基于实时数据的融合后的统计量SIS_new和DIS_new：

(8)判断SIS_new＞1、DIS_new＞1是否成立，若SIS_new＞1成立，则说明工作点发生了变化；若DIS_new＞1成立，则说明离线建立的动态内模型不再适用；若SIS_new＞1、DIS_new＞1均不成立，则返回步骤(5)；

(a)初始化动态阶次s为1；

(c)按照步骤(3)中的动态内部慢特征分析算法确定动态潜变量和残差；

(d)计算残差矩阵E中各维度的自相关系数，其中第i维数据的自相关系数为autocurr_ei(k)，k为相位差，对于所有维度的autocurr_ei(k)，若存在autocurr_ei(k)＞0.5，则l＝l+1；计算潜变量T中各维度的自相关系数，其中第i维数据的自相关系数为autocrr_Ti(k)，k为相位差，对于所有维度的autocurr_Ti(k)，若存在autocurr_Ti(k＜0.2，则l＝l-1；

(e)重复步骤(c)-(d)，直到确定合适的参数l；

(f)以丙烯精馏塔正常工作状态下的生产过程数据作为验证集，在验证集上检验残差矩阵E是否包含明显的自相关关系，若是，修改动态阶次s＝s+1；

(g)重复步骤(c)-(f)，直到确定合适的参数s。

5.如权利要求4所述的基于动态内部慢特征分析算法的丙烯精馏塔状态监控方法，其特征在于，步骤(1)中，采集的丙烯精馏塔的生产过程数据的样本数不少于1000组。

6.如权利要求4所述的基于动态内部慢特征分析算法的丙烯精馏塔状态监控方法，其特征在于，步骤(2)中，标准化后的数据各维度数据均值为0，方差为1。