一种确定复杂工业过程中关键变量因果关系的方法

摘要

本发明公开了一种确定复杂工业过程中关键变量因果关系的方法，通过对工业过程中待确定因果关系的关键变量的历史数据采用伪近邻思想计算每个关键变量的最佳时序嵌入维度；针对两两关键变量，假定因果关系，以假定因变量的最佳时序嵌入维度为标准，构造两个关键变量的时序重构流形，利用收敛交叉映射算法计算两者间的收敛交叉映射能力；基于蒙特卡洛模拟确定CCM的能力判定阈值，以此判定关键变量间假定因果关系的正确性，从而构建工业过程中关键变量的初步因果关系网络；利用时滞检测方法修正初步因果关系网络，得到最终关键变量因果关系网络。本发明充分利用生产离线数据，对生产过程无干扰作用，提高了安全性和经济效益。

说明书

技术领域

本发明涉及工业过程异常工况诊断与溯因分析中的关键变量因果关系分析等技术领域，具体涉一种确定复杂工业过程中关键变量因果关系的方法。

背景技术

在石油炼制、钢铁冶炼等大型复杂工业系统，流程中运行变量较多，各变量之间耦合严重，集成度高。例如，加氢裂化过程是石油炼制流程中一个子过程，是氢气在较高压力和温度下经催化剂作用使重质油发生加氢、裂化和异构化反应转化为轻质油的加工过程，主要包括加氢精制反应、加氢裂化反应、高低压分离以及分馏系统这四个重要部分。加氢裂化的工业过程作为这类典型的庞大而复杂的流程工业，具有流程长、变量多、工况多以及非线性的特点。当流程中一个关键变量出现异常时，可能会引起相关变量出现异常，从而导致整个流程出现波动，然而操作人员根据经验得到的结果往往无法及时并准确地定位到异常发生的原因，使工业过程长时间处于非优状态，造成资源能源的浪费。

在异常工况诊断研究中，有关学者提出了基于定性模型、基于定量模型及基于数据驱动的方法。其中，基于定量模型和定性模型是从对流程物理特性的基本理解发展而来，依靠模型的先验知识，但它只适用于机理明确的工业过程。目前，在流程工业中存有海量的历史数据，因而基于数据驱动的诊断方法对比其他两类方法具有明显的优势，而通过这类方法分析异常传播路径、定位异常原因的核心是寻找变量之间的因果关系，但是由于工业过程中变量具有非线性、耦合严重等特性，变量之间即使存在相关性也不能保证变量之间存在因果关系，给流程中变量间的因果关系识别带来了挑战。

针对上述挑战，目前已有不少学者从不同角度出发提出了一些解决方法。格兰杰提出的Granger因果关系检验主要适合于变量之间存在线性相关性，同时还要求测试的两个变量相互独立，但是该方法对于耦合动态系统不适用，虽然一些学者改进并提出了非线性Granger因果关系分析方法，但是该方法在不可分离非线性动态系统中测试失败。互相关函数(Cross Correlation Function)通过计算变量之间最大CCF值得到时间延迟来确定信号的因果关系，但是这个方法适用于线性系统，对于非线性系统计算结果不能保证准确性。Schreiber提出的传递熵(Transfer Entropy)算法来判断因果关系虽然对非线性和线性系统都能适用，但是不仅对参数改变敏感，而且计算量大、实现相对困难，同时结果计算的准确性与概率密度函数的估计密切相关，因此要求数据具有时间稳定性。此外，还有学者提出通过输入激励信号来监测模型的变化，但是这类方法对于实际工业流程来说，不仅难以设计激励信号，而且加入的激励信号本身会影响工业流程的正常运行。

针对上述方法的局限性以及复杂工业过程的实际特点，亟需提出一种用于确定复杂工业过程中关键变量间因果关系的方法。考虑Sugihara等提出的收敛交叉映射(CCM)方法不仅适用于非线性、中强耦合的复杂系统，而且计算方便参数选择量少的特点，本专利采用收敛交叉映射(CCM)的思想确定关键变量间的因果关系。然而由于CCM初始应用于生态系统，并没有对时序嵌入维度的选择进行讨论，而且对于具有高维空间的复杂工业过程，依据人工经验选择嵌入维度和CCM能力判定阈值无法确保结果的准确性，并且很难实现自动监测。因此本专利通过改进CCM方法确定复杂工业过程中关键变量间因果关系的方法，最终形成关键变量的因果关系网。

发明内容

本发明的目的在于提供一种计算方便，参数选择量少，充分考虑了关键变量之间的耦合问题，能解决具有中强耦合的非线性复杂动态系统中关键变量间的因果关系识别，为工业过程异常工况的溯因分析的方法。

为了实现上述发明目的，本发明提供技术方案如下：一种确定复杂工业过程中关键变量因果关系的方法，所述方法包含以下步骤：

步骤一：收集工业过程关键变量的生产历史数据，对样本进行数据预处理：

选取待确定因果关系的n个关键变量，收集这些关键变量一段连续的时间序列的生产历史数据组成一个变量维度为n，样本长度为L的数据集作为样本；然后利用小波变换对所述的样本进行去噪处理；所述的n为大于1的自然数；

步骤二：确定重构流形的最佳时序嵌入维度E：

对步骤一预处理后的样本中每个关键变量，采用伪近方法计算其最佳时序嵌入维度E；即对所述的样本中的关键变量的第i个重构流形，设嵌入维度E＝k的最近邻点的距离为D_i(k)，嵌入维度为E＝k+1的最近邻点的距离为D_i(k+1)，计算两者的几何平均，通过对数形式比较两者的值，使其结果更具有鲁棒性；则有

其中，L表示时间序列的样本长度，τ为选取的时间间隔；当流形中的最近邻点的距离随着嵌入维度的增加基本不变时，即D_i(k)无限接近于D_i(k+1)，有G(k)的值无限接近于0，称为饱和状态。当G(k)的值基本不变时，选取最先进入饱和状态时刻对应的嵌入维度E为重构流形的最佳时序嵌入维度；

步骤三：计算两两关键变量间的收敛交叉映射能力：

对于两两关键变量，假定因果关系，并以假定因变量的最佳时序嵌入维度E为标准，构造两关键变量的重构流形，然后采用k近邻算法计算得到最近重构流形坐标，利用CCM算法计算两者间的收敛交叉映射能力；

步骤四：基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络：

采用蒙特卡洛抽样，选取多对不相关的随机时间序列，分别计算它们在不同样本长度L下收敛交叉映射能力的均值和标准差；然后，曲线拟合得到与样本长度L相关的均值函数和标准差函数σ_r＝1.17L^-0.46；接着，根据确定的样本长度L在拟合函数中找到对应的均值与标准差；以此按照3σ准则得到收敛交叉映射能力的判定阈值r₀＝-3.8L^-1.08+3.51L^-0.46；若步骤三中计算的某两关键变量收敛交叉映射能力的收敛值超出该判定阈值，则认为假定的关键变量因果关系成立，遍历所有关键变量，判定所有关键变量两两间的因果关系，构建初步因果关系网络。

进一步地，还包括步骤五，采用时滞检测方法修正初步因果关系网络，得到最终的工业过程中关键变量因果关系网络：

首先，从步骤四中得到的初步因果关系网络中挑选出双向因果关系的关键变量：根据工业过程的工艺大致确定关键变量之间可能存在的时间延迟范围，设为Δt₁，Δt₂，其中Δt₁≤0，Δt₂≥0，假设时间序列X_i和X_j，利用CCM算法通过X_i(t)来反映X_j(t-Δt)，其中Δt∈[Δt₁，Δt₂]，计算不同时滞下的收敛交叉映射能力；比较不同时滞下的收敛交叉映射能力，得到的最大收敛交叉映射能力为最优结果，此最优结果对应的时滞认为是在流程中两关键变量之间的时间延迟；如果最优收敛交叉映射能力对应的时滞为正，则表示时间序列X_i的未来值更能准确的反映X_j的过去值，因此存在从X_j到X_i方向上的因果关系，反之，时滞为负则表示不存在；

然后，在消除了假象双向因果关系的网络中挑选出存在相同原因且两者之间存在因果关系的关键变量，如果变量X对变量Z的影响先于X对Y的影响，此时可理解为Z的过去值更能准确的预测Y的未来值，此时判断最优收敛交叉映射能力对应的时滞为负；从步骤四中得到的初步因果关系网络中经过以上分析后，对步骤四中得到的初步因果关系网络进行修正，得到最终的工业过程中关键变量因果关系网络。

进一步地，步骤三中所述的“利用CCM算法计算关键变量两两间的收敛交叉映射能力”的过程为：

对于两两含有时序的关键变量X_i和X_j，假定X_j到X_i存在因果关系；以假定因变量X_j的最佳时序嵌入维度为标准，构造变量X_i和X_j的重构流形；对于长度为L的样本时间序列X_i和X_j，其中i，j＝1，2，3，......，n且i≠j，则时间序列在t时刻的重构流形坐标分别为：

M_i(t)＝[x_i(t)，x_i(t-τ)，......，x_i(t-(E-1)τ)]

M_j(t)＝[x_j(t)，x_j(t-τ)，......，x_j(t-(E-1)τ)]

由上式可得时间序列X_i和X_j的重构流形分别为M_i＝{M_i(t)}，M_j＝{M_j(t)}；从流形M_i中找到距离M_i(t)最近的E+1个点M_i，t＝{M_i(t₁)，M_i(t₂)，M_i(t₃)，......，M_i(t_E+1)}，为了降低计算量，提高计算速度，采用k近邻算法来寻找这E+1个点；根据时间序列X_i和X_j存在从X_j到X_i的因果关系的假设，通过M_i，t来估计时间序列X_j的值，得到估计的时间序列其具体计算公式如下：

D＝∑d_k>

计算估计时间序列与原时间序列X_j的相关系数得到收敛交叉映射能力r_i→j如下：

其中，表示时间序列X_j的均值，表示估计时间序列的均值。根据CCM算法原理可得，如果存在从X_j到X_i的因果关系，则随着样本时间序列长度L的增加，相关系数r_i→j将收敛于一个大于0的值，此时，称该相关系数为收敛交叉映射能力；如果收敛交叉映射能力收敛于大于0的值，则认为两个关键变量之间存在因果关系，否则不存在。

进一步地，步骤一中所述的“一段连续的时间序列”为相对稳定的、无人工调整的一段时间。

进一步地，所述的L为大于或等于200的自然数。

进一步地，所述的工业过程为石油炼制的加氢裂化过程，所述的关键变量包括但不限于出口温度、入口温度、一床层顶温、二床层顶温、三床层顶温、冷氢注入一二床层量、冷氢注入二三床层量。

进一步地，所述的工业过程为烟草制丝的烘丝过程过程，所述的关键变量包括但不限于入口水分、热风温度、热风风量、一区筒壁压力、二区筒壁压力、出口水分这六个关键变量。

与现有技术相比，本发明的有益效果在于：一是通过改进CCM的方法，能解决非线性中强耦合的复杂动态系统中关键变量间的因果关系识别问题，为分析复杂工业过程中异常工况的溯因分析提供重要条件。二是采用伪近邻的思想确定CCM算法中的重要调优参数(时序嵌入维度)，解决了由于人工选择工业过程中多个关键变量的调优参数而导致算法准确性降低的问题。三是基于蒙特卡洛模拟方法确定CCM的判定阈值，减少因果分析中的人为判断。四是通过时滞检测消除假象双因果关系以及部分同源因果造成的伪因果关系，降低了工业过程关键变量因果关系网络的干扰项。依据本专利构建的工业过程中关键变量因果关系网络，对比现有方法，适用于复杂的工业实际过程，计算更加简便，充分利用工业过程的离线数据，且对工业过程无干扰作用，提高了安全性和经济效益，对回溯异常传播路径、定位异常原因效果良好。

附图说明

图1为本发明的一个实施例的算法流程图。

图2为本发明的一个实施例加氢裂化流程精制反应装置图。

图3为本发明的一个实施例嵌入维度曲线图。

图4为本发明的一个实施例交叉映射能力曲线。

图5为本发明的一个实施例初步因果关系网络图。

图6为本发明的一个实施例最终因果关系网络图。

图7为本发明的另一个实施例的嵌入维度曲线图。

图8为本发明的另一个实施例的交叉映射能力曲线.。

图9为本发明的另一个实施例的初步因果关系网络。

图10为本发明的另一个实施例的最终因果关系网络。

具体实施方式

为了更好地说明本发明的研发目的、技术方案及优势，以及充分公开本发明技术方案，下面对本发明的技术方案做进一步说明。

本发明提供的技术适用于石油炼制、烟草制丝、钢铁冶炼等复杂工业过程中关键变量间因果关系的确定。

一种确定复杂工业过程中关键变量因果关系的方法，包含以下步骤：

步骤一：收集工业过程关键变量的生产历史数据，对样本进行数据预处理：

选取待确定因果关系的n(n大于等于1)个关键变量，收集这些关键变量一段连续的时间序列的生产历史数据组成一个变量维度为n，样本长度为L的数据集作为样本；然后利用小波变换对所述的样本进行去噪处理。

为了使选取的样本更具稳定性，生产历史数选取相对稳定的、无人工调整的一段时间产生的数据。

步骤二：确定重构流形的最佳时序嵌入维度E：

对步骤一预处理后的样本中每个关键变量，采用伪近方法计算其最佳时序嵌入维度E；即对所述的样本中的关键变量的第i个重构流形，设嵌入维度E＝k的最近邻点的距离为D_i(k)，嵌入维度为E＝k+1的最近邻点的距离为D_i(k+1)，计算两者的几何平均，通过对数形式比较两者的值，使其结果更具有鲁棒性。则有

其中，L表示时间序列的样本长度，τ为选取的时间间隔；当流形中的最近邻点的距离随着嵌入维度的增加基本不变时，即D_i(k)无限接近于D_i(k+1)，有G(k)的值无限接近于0，称为饱和状态。当G(k)的值基本不变时，选取最先进入饱和状态时刻对应的嵌入维度E为重构流形的最佳时序嵌入维度。

步骤三：计算两两关键变量间的收敛交叉映射能力：

对于两两关键变量，假定因果关系，并以假定因变量的最佳时序嵌入维度E为标准，构造两变量的重构流形，然后采用k近邻算法计算得到最近重构流形坐标，利用CCM算法计算两者间的收敛交叉映射能力；

对于两两含有时序的关键变量X_i和X_j，假定X_j到X_i存在因果关系；以假定因变量X_j的最佳时序嵌入维度为标准，构造关键变量X_i和X_j的重构流形；对于长度为L的样本时间序列X_i和X_j，其中i，j＝1，2，3，......，n且i≠j，则时间序列在t时刻的重构流形坐标分别为：

M_i(t)＝[x_i(t)，x_i(t-τ)，......，x_i(t-(E-1)τ)]

M_j(t)＝[x_j(t)，x_j(t-τ)，......，x_j(t-(E-1)τ)]

由上式可得时间序列X_i和X_j的重构流形分别为M_i＝{M_i(t)}，M_j＝{M_i(t)}；从流形M_i中找到距离M_i(t)最近的E+1个点M_i，t＝{M_i(t₁)，M_i(t₂)，M_i(t₃)，......，M_i(t_E+1)}，为了降低计算量，提高计算速度，采用k近邻算法来寻找这E+1个点；根据时间序列X_i和X_j存在从X_j到X_i的因果关系的假设，通过M_i，t来估计时间序列X_j的值，得到估计的时间序列其具体计算公式如下：

D＝∑d_k>

计算估计时间序列与原时间序列X_j的相关系数得到收敛交叉映射能力r_i→j如下：

其中，表示时间序列X_j的均值，表示估计时间序列的均值。根据CCM算法原理可得，如果存在从X_j到X_i的因果关系，则随着样本时间序列长度L的增加，相关系数r_i→j将收敛于一个大于0的值，此时，称该相关系数为收敛交叉映射能力；如果收敛交叉映射能力收敛于大于0的值，则认为两个关键变量之间存在因果关系，否则不存在。

步骤四：基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络：

采用蒙特卡洛抽样，选取多对不相关的随机时间序列，分别计算它们在不同样本长度L下收敛交叉映射能力的均值和标准差；然后，曲线拟合得到与样本长度L相关的均值函数和标准差函数σ_r＝1.17L^-0.46；接着，根据确定的样本长度L在拟合函数中找到对应的均值与标准差；以此按照3σ准则得到收敛交叉映射能力的判定阈值r₀＝-3.8L^-1.08+3.51L^-0.46；若步骤三中计算的某两关键变量收敛交叉映射能力的收敛值超出该判定阈值，则认为假定的关键变量因果关系成立，遍历所有关键变量，判定所有关键变量两两间的因果关系，构建初步因果关系网络。

为了消除假象双因果关系以及部分同源因果造成的伪因果关系，降低了工业过程关键变量因果关系网络的干扰项，对前述技术方案的一种改进为：还包括：

步骤五，采用时滞检测方法修正初步因果关系网络，得到最终的工业过程中关键变量因果关系网络：

首先，从步骤四中得到的初步因果关系网络中挑选出双向因果关系的关键变量：根据工业过程的工艺大致确定关键变量之间可能存在的时间延迟范围，设为Δt₁，Δt₂，其中Δt₁≤0，Δt₂≥0，假设时间序列X_i和X_j，利用CCM算法通过X_i(t)来反映X_j(t-Δt)，其中Δt∈[Δt₁，Δt₂]，计算不同时滞下的收敛交叉映射能力；比较不同时滞下的收敛交叉映射能力，得到的最大收敛交叉映射能力为最优结果，此最优结果对应的时滞认为是在流程中两关键变量之间的时间延迟；如果最优收敛交叉映射能力对应的时滞为正，则表示时间序列X_i的未来值更能准确的反映X_j的过去值，因此存在从X_j到X_i方向上的因果关系，反之，时滞为负则表示不存在；

然后，在消除了假象双向因果关系的网络中挑选出存在相同原因且两者之间存在因果关系的关键变量，如果变量X对变量Z的影响先于X对Y的影响，此时可理解为Z的过去值更能准确的预测Y的未来值，此时判断最优收敛交叉映射能力对应的时滞为负；从步骤四中得到的初步因果关系网络中经过以上分析后，对步骤四中得到的初步因果关系网络进行修正，得到最终的工业过程中关键变量因果关系网络。

具体实施例1：

下面将结合我国某炼油厂加氢裂化过程的实际运行数据举例说明，针对本发明中提出的基于改进收敛交叉映射算法确定加氢裂化过程中关键变量因果关系，根据图1中流程图执行，并做出详细说明。应该强调的是，下述说明仅仅是示例性的，且下述中实例只是加氢裂化流程中一部分关键变量，而不是为了限制本发明的范围及其应用。所提供方法包括如下步骤：

步骤一，收集加氢裂化流程关键变量的生产历史数据，进行样本数据预处理。即，对于待确定因果关系的加氢裂化流程加氢精制反应部分7个关键变量，其中加氢精制反应装置简化示意图如图2所示，从历史数据库中选取同一工况下非人工干预的连续时间序列来组成数据集长度为1000的样本，并利用小波变换对样本进行去噪处理。表1为待测的加氢裂化流程中部分关键变量的描述。

表1关键变量描述

步骤二，对步骤一得到的预处理后的每个关键变量，采用伪近邻思想得到重构流形的最佳时序嵌入维度E。由于提取加氢裂化流程关键变量的时间序列本身是离散的点，采样间隔为每分钟，因此选择时间间隔为τ＝1。计算得到各关键变量的G(k)图形如图3所示。根据G(k)曲线，当其变化速率小于0.02时，认为G(k)的值趋于饱和得到最佳嵌入维度。表2为通过伪近邻计算方法得到各关键变量的重构流形的最佳嵌入维度。

表2关键变量最佳嵌入维度

步骤三，计算两两关键变量间的收敛交叉映射能力。对于两两关键变量，假定因果关系，并以假定因变量的最佳嵌入维度E为标准，构造两变量的重构流形，得到重构流形，采用k近邻算法得到最近重构流形坐标。利用CCM算法计算关键变量两两间的收敛交叉映射能力，在不同样本长度下计算收敛交叉映射能力得到的收敛曲线，然后对比曲线是否收敛得到收敛结果，图4给出两两关键变量间收敛的CCM结果曲线。需注意的是，由于计算的两关键变量之间可能存在嵌入维度不一致的情况，此时选择原因变量的最佳时序嵌入维度来重构流形。表3为关键变量两两之间计算得到的收敛交叉映射能力，其中行代表原因变量，列代表结果变量，不同关键变量之间的“/”表示计算的CCM值不收敛。

表3收敛交叉能力

步骤四，基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络。采用假设检验的方式，使用多对不相关的随机时间序列，分别计算在不同样本长度下收敛交叉映射能力的均值和标准差，曲线拟合得到与样本长度相关的均值和标准差函数，并根据3σ准则得到收敛交叉映射能力的阈值r₀＝-3.8L^-1.08+3.51L^-0.46。在实例中，样本长度L＝1000，因而得到收敛交叉映射能力的阈值为r₀＝0.1441，判断表3中的数据，得到加氢裂化流程关键变量的初步因果关系网络如图5所示。

步骤五，采用时滞检测修正初步因果关系网络。在加氢裂化流程中，前后变量之间存在时间延迟，如果上游变量的当前值改变，经过一段时间后才能影响到下游变量的值。因此可通过时滞检测区分出由于强单向因果关系引起的同步现象以及部分同源因果引起的伪因果关系。根据加氢裂化流程工艺大致确定两关键变量之间可能存在的时间延迟范围，从初步因果关系网络中可得到存在双向因果关系的关键变量，经过时滞检测消除假象双向因果关系后得到进一步的因果关系网络，再从中挑选出存在相同原因且两者之间存在因果关系的关键变量，经过时滞检测后消去变量之间的伪因果关系。关键变量之间的时滞正负结果如表4所示。依据时滞检测结果修正初步因果关系网络后得到最终的因果关系网络如图6所示。

表4关键变量时滞检测结果

经过验证，得到的最终因果关系网络与加氢裂化流程工艺了解的经验一致，符合实际情况。该结果与人工通过流程经验知识判别因果关系相比，不仅具有快速性而且降低了操作人员的工作量。本专利完成的加氢裂化流程因果关系网络，对比现有方法，计算更加简便，充分利用加氢裂化流程的离线数据，且对加氢裂化流程无干扰作用，提高了安全性和经济效益，进而为回溯异常传播路径、定位异常原因打下良好基础。

具体实施例2：

再以烟草制丝生产为例，烟草制丝生产是典型的间歇生产过程，烘丝是最关键的一道工序，它的目标是把烟叶等原材料加工成符合卷烟质量标准以及下一阶段卷制生产要求的叶丝。该生产过程较为复杂，且叶丝受到水分、温度等方面因素的影响，因此，以烘丝过程为例，分析该过程中关键变量之间的因果关系，主要步骤如下：

步骤一，选取烘丝过程中入口水分、热风风量、一区筒壁压力、二区筒壁压力、出口水分这五个关键变量，样本的长度为300；

步骤二，确定重构流形的最佳时序嵌入维度E，各关键变量的G(k)图形如图7所示，因而得到各关键变量的最佳时序嵌入维度分别为7、6、6、6、5；

步骤三，计算得到关键变量两两间的收敛交叉映射能力，判断是否收敛，得到收敛的曲线图如图8所示，表5为关键变量两两之间计算得到的收敛交叉映射能力，其中行代表原因变量，列代表结果变量，不同关键变量之间的“/”表示计算的CCM值不收敛。

表5收敛交叉能力

步骤四，判定阈值判定得到初步因果关系网络图如图9所示；

步骤五，进行时滞检测结果，得到最终因果关系网络如图10所示，结果符合实际情况。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的保护范围内。