一种基于马尔可夫随机场模型的工业场景工况聚类方法、系统、设备及存储介质

技术领域

本发明涉及一种工业场景下工况的构建方法、系统、设备及存储介质，属于工业数据分析技术领域。

背景技术

在工业大数据分析领域，实际工业信号往往来自许多生产系统传感器，且这些数据多随时间变化。工况定义为将该高维时序信号划分的若干可能的状态。如热电厂的运行状态大致可分为启动段、运行段、停止段、异常段等。工况挖掘和识别通过机器学习算法挖掘高维数据中隐藏的工况信息，这对实际生产有一定的指导意义，同时也是许多后续算法(预测、分类)重要的前处理手段。

现有的工况聚类算法主要有基于距离度量、基于数据概率分布层次聚类等方法。基于距离度量的聚类算法的主要问题在于，无论是采用何种距离度量以及聚类算法，都需要对原始信号进行切割，且要求切割后的子序列长度一致。过长的子序列可能包含不止一个的隐藏工况，而过短的子序列则可能包含较少的数据特征，都对聚类算法的结果有不利影响，且人为对序列进行分割也不利于保留原时间序列的时间一致性。基于数据分布的层次聚类方法则直观地根据数据在某些维度上的分布特征确定分割点进行聚类，优势在于算法直观，可解释性强，不足之处在于用于层次聚类的特征往往是分布十分集中的类别特征，不平稳或波动性大的特征则不能用于该类方法，使用较为受限。

因此亟需一种工业场景下的工况聚类方法，以挖掘工业场景下的潜在工况。

发明内容

本发明为了高精度、高效率的实现工业场景下的潜在工况挖掘，在下文中给出了关于本发明的简要概述，以便提供关于本发明的某些方面的基本理解。应当理解，这个概述并不是关于本发明的穷举性概述。它并不是意图确定本发明的关键或重要部分，也不是意图限定本发明的范围。

本发明提供了一种基于马尔可夫随机场模型的工业场景工况聚类方法，包括以下步骤：

步骤S100、从发电厂众多特征中初步筛选出能够影响发电厂运行工况变化的特征，并对这些能够影响发电厂运行工况变化的特征进行初步处理；

步骤S200、将步骤S100初步处理后的影响发电厂运行工况变化特征进行降噪处理；

步骤S300、将步骤S200中进行降噪处理后的影响发电厂运行工况变化特征数据进行降维处理；

步骤S400、将降维处理后的数据进行聚类，以得到运行工况。

进一步地：步骤S100中，在燃煤锅炉发电厂已有所有测点数据中，采集以5秒为时间周期，总时间长度为4个月的数据中，选择影响锅炉工况变化的特征。

进一步地：步骤S200中，对影响发电厂运行工况变化特征进行降噪处理方法是：首先将异常数据进行修正，再对高频噪声进行抑制，具体方式为：

异常数据的处理方法：

自定义时间窗口值w，从X

噪声的处理方法是：

自定义时间窗口值w，从X

通过高斯滤波、小波降噪算法或数字低通滤波器将输入数据的高频噪声进行抑制，通过小波软阈值降噪方法将输入数据的数据异常进行修正。

进一步地：步骤S300中，对影响发电厂运行工况变化特征数据进行降维处理具体方法如下：

步骤S310，计算任意两个影响发电厂运行工况变化特征数据的相关系数，当相关系数大于0.95时，剔除这两个特征中方差较小的特征，即剔除每个磨煤机的电机转速反馈信号、风粉一次风速、风粉混合物温度、混合风量特征，燃烧室上下端部的二次风层风门开度特征，此时数据特征降记为X∈R

步骤S320，将X∈R

进一步地：步骤S400中，将步骤S300中影响发电厂运行工况变化特征数据进行降维处理后的数据X∈R

步骤S410、将输入数据X∈R

步骤S420、根据聚类结果P及参数Θ，得到能够表征聚类结果的值，当表征值的变化小于一定阈值(根据表征聚类结果的具体要求所设定阈值)，或迭代次数达到最大迭代次数，则停止迭代，此时的聚类结果为最终输出结果。

进一步地：还包括步骤S430、根据整体聚类优化目标对现有聚类及其参数进行迭代更新，参数更新分为两步，首先在Θ不变情况下更新P，在Θ不变情况下，改写优化目标，采用动态规划算法求解优化目标，对P进行更新；其次在P不变情况下更新聚类Θ,在P不变情况下，改写优化目标，采用交替方向乘子算法求解优化目标，对Θ、U、Z进行更新。

进一步地：所述步骤S410中，进行数据X∈R

步骤S411、输入数据形式是维度为n，长度为T的多元时间序列，可表示为x

根据实际需要自定义聚类数量K、子序列分割的窗长w、最大迭代次数max_iter、正则项系数β、π；

步骤S412、将x

X

步骤S413、根据自定义聚类数量K，进行随机初始聚类，其结果表示为：

步骤S414、对每个聚类定义其高斯逆协方差矩阵Θ

其中A

步骤S415、对每个Θ

因此，参数初始化模块将输入数据映射为子序列数组X

进一步地：在步骤S420中，高斯逆协方差矩阵Θ优化问题遵循以下3个条件：

(1)假定X

(2)增加惩罚项，鼓励连续的时间点分类给相同的聚类；

(3)避免过拟合，希望Θ尽可能稀疏；

基于以上三点，可以得出问题的优化目标：

其中ll(X

ll(X

因此聚类结果表征值表达式为：

进一步地：步骤S430中根据整体聚类优化目标对现有聚类及其参数进行迭代更新，作用是根据优化目标：

对P＝{P

更新分为两步，第一步：在Θ不变情况下更新P＝{P

更新第一步：

在Θ不变情况下，优化目标为：

解该优化问题可以等价于：最大化对数似然和聚类在时间上尽可能连续之间的权衡。时间连续性的重要性由参数β控制。

该优化问题可以采用动态规划算法，动态规划算法有向图如下图所示；

在动态规划算法视角下，该最优化问题等价于：求解有向图中从1到T时刻权值最短的路径，图中结点(i，j)的值为当X

因此，通过以上方法，可在Θ不变情况下，采用动态规划算法求解优化问题，对P＝{P

更新第二步：

在P不变情况下，极大似然项可以改写成：

由于Θ

求解该优化问题，每次都需要至少重复K次，总体算法复杂程度较高。为降低计算复杂度，采用具有并行化计算特性的交替方向乘子算法(ADMM)以解决该优化问题。

采用ADMM算法，引入一致性变量Z以及对偶变量U，则优化问题可改写为：

其对应的增广拉格朗日函数为：

其中ρ＞0，是ADMM算法的惩罚系数。

ADMM算法执行以下三个步骤重复执行直至收敛，停止条件设定为Θ和Z在两次迭代结束时的结果残差小于阈值：

(a)

(b)

(c)U

上述迭代步骤中，步骤(a)和(b)属于Θ和Z优化函数，解优化函数能够得到新的Θ和Z。根据新的Θ和Z，可以通过步骤(c)得到新的U。Θ和Z的更新方法如下所述。

Θ更新：Θ参数的更新与Z和U变量无关，可以写作下式：

该优化问题有一个已知的解析解：

其中QDQ

Z更新：Z参数更新可以写作下式：

定义变量

令

可以得到以下Z的封闭解：

因此，通过更新第一步和更新第二步，能够得到新的P＝{P

一种基于马尔可夫随机场模型的工业场景工况聚类系统，包括初步筛选模块、噪声抑制模块、特征降维模块和工况聚类模块；

所述初步筛选模块用于从发电厂众多特征中初步筛选出能够影响发电厂运行工况变化的特征，并对这些特征进行初步处理；

所述噪声抑制模块用于对电厂运行数据中存在的异常数据进行筛选，并对噪声进行修正；

所述特征降维模块用于对进行降噪处理后的影响发电厂运行工况变化特征数据进行降维处理；

所述工况聚类模块将降维处理后的数据进行聚类，以得到运行工况。

进一步地：所述工况聚类模块包括参数初始化子模块、聚类赋值子模块以及参数更新子模块；

所述参数初始化子模块用于将输入数据进行初始化，并自定义聚类参数，形成特定形式的数据，以用于后续聚类运算；

所述聚类赋值子模块是根据聚类结果参数，得到能够表征聚类结果的值；

所述参数更新子模块根据整体聚类优化目标对现有聚类及其参数进行迭代更新。

一种基于马尔可夫随机场模型的工业场景工况聚类设备，包括存储器和处理器，存储器存储有计算机程序，所述的处理器执行所述计算机程序时实现前述公开的一种基于马尔可夫随机场模型的工业场景工况聚类方法的步骤。

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现前述公开的一种基于马尔可夫随机场模型的工业场景工况聚类方法。

本发明具有以下有益效果：

1、本发明能够平稳的实现工业生产机组负荷稳态工况的聚类，通过对工业数据的筛选、抑制、降维和聚类方式对机组海量数据的稳态工况的提取，为后续的机组性能指标计算和设备运行状态诊断分析提供了有效的，有价值的数据支撑；

2、本发明进一步提高了处理和分析机组海量数据的效率，实现了人工处理海量机组数据时所无法达到的高效性和经济性。

附图说明

图1是一种基于马尔可夫随机场模型的工业场景工况聚类系统图；

图2是一种基于马尔可夫随机场模型的工业场景工况聚类方法流程图；

图3是动态规划算法有向图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图中示出的具体实施例来描述本发明。但是应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

具体实施方式一：

本实施方式的目的在于提供一种基于马尔可夫随机场模型的工业场景工况聚类方法与系统。通过在发电厂各个位置布置传感器从而获得的发电厂实时运行数据，经过人工筛选尽可能选出能够影响锅炉工况变化的特征数据，经过数据清洗、噪声抑制、特征降维等方法对数据进行处理，之后对数据进行聚类，每个类别视为不同工况。本实施方式主要分为初步筛选模块；噪声抑制模块；特征降维模块；工况聚类模块；

初步筛选模块

本模块作用是从发电厂众多特征中初步筛选出能够影响发电厂运行工况变化的特征，并对这些特征进行初步处理。

在燃煤锅炉发电厂已有所有测点采集的一定的时间周期(例如5秒)、总时间长度为2个月的数据中，选择能够影响锅炉工况变化的特征。特征包括总给煤量、总给水量、总给风量、锅炉负荷、磨煤机数据、燃烧室数据。磨煤机数据包括每台磨煤机的电机转速反馈信号、给煤量反馈信号、风粉一次风速、风粉混合物温度、混合风量，燃烧室数据包括sofa风层、上下端部一次风层、二次风层的四个角的风门开度，共106个特征。

对筛选出来的特征数据进行初步处理，包括根据总给水量和总给煤量的比值得到水煤比，根据总给风量和总给煤量的比值得到风煤比，将燃烧室每个风层的四角风门开度数据取平均值得到每个风层的平均风门开度。之后剔除总给水量、总给风量、每个风层的四角风门开度。保存为初步筛选数据。初步筛选数据维度为51个。

保存初步筛选数据，作为本模块的输出，并作为噪声抑制模块的输入。

噪声抑制模块

本模块作用是将输入数据中存在的噪声进行降噪。实际电厂运行数据中存在各类噪声，噪声可分为由于工业传感器自身波动、由于工业传感器偶然失真产生的数据异常。

本模块根据标准差筛选出异常数据，并用异常数据周围数据均值对异常数据进行替换；采用高斯滤波方法对噪声进行修正。

经过噪声抑制的数据作为本模块的输出，并作为特征降维模块的输入。

特征降维模块

本模块作用是将输入数据进行降维处理，以降低计算难度，节省计算时间。

本模块通过相关系数法、方差选择法将呈现正线性相关的冗余特征进行剔除，可剔除每个磨煤机的电机转速反馈信号、风粉一次风速、风粉混合物温度、混合风量特征，还可以剔除燃烧室上下端部的二次风层风门开度特征，将维度将至19个；再通过主成分分析(PCA)、奇异值分解(SVD)算法对特征进行降维，将维度将至8个。经PCA和SVD处理的数据的特征不再具有确切的物理含义。

经过特征降维的数据作为本模块的输出，并作为工况聚类算法的输入。

工况聚类模块

本模块作用是将输入数据进行聚类。根据聚类结果得到每个类别的时间戳，即可将锅炉运行数据分为不同工况。

本模块分为参数初始化子模块、聚类赋值子模块以及参数更新子模块。

参数初始化子模块作用是将输入数据进行初始化，并自定义聚类参数，形成特定形式的数据，以用于后续聚类运算。根据自定义聚类参数，生成聚类结果P；对聚类P定义了具有托普利兹矩阵性质的高斯逆协方差矩阵Θ；根据Θ随机初始化托普利兹矩阵U、Z，作为后续计算的参数。

聚类赋值子模块作用是根据聚类结果P及参数Θ，得到能够表征聚类结果的值。当表征值的变化小于一定阈值，或迭代次数达到最大迭代次数，则停止迭代，此时的聚类结果为最终输出结果。

参数更新子模块作用是根据整体聚类优化目标对现有聚类及其参数进行迭代更新。参数更新分为两步，首先在Θ不变情况下更新P，其次在P不变情况下更新聚类Θ。在Θ不变情况下，改写优化目标，采用动态规划算法求解优化目标，对P进行更新。在P不变情况下，改写优化目标，采用交替方向乘子算法(ADMM)求解优化目标，对Θ、U、Z进行更新。每一次更新聚类结果以及参数，即为一次迭代。更新的聚类结果以及参数传入聚类赋值子模块。

具体实施方式二：

本实施方式是将通过在发电厂各个位置布置传感器从而获得的发电厂实时运行数据，经过数据清洗、特征降维等方法将数据进行基本处理，之后输入工况聚类模块，将锅炉数据进行聚类，每个类别视为不同工况。其流程如图1和图2所示，其中：

一、初步筛选模块

初步筛选模块作用是从发电厂众多特征中初步筛选出能够影响发电厂运行工况变化的特征，并对这些特征进行初步处理。

在燃煤锅炉发电厂已有所有测点采集的一定的时间周期(例如5秒)、总时间长度为4个月的数据中，选择能够影响锅炉工况变化的特征。特征包括总给煤量、总给水量、总给风量、锅炉负荷、磨煤机数据、燃烧室数据。磨煤机数据包括每台磨煤机的电机转速反馈信号、给煤量反馈信号、风粉一次风速、风粉混合物温度、混合风量，燃烧室数据包括sofa风层、上下端部一次风层、二次风层的四个角的风门开度，共106个特征。

对筛选出来的特征数据进行初步处理，包括根据总给水量和总给煤量的比值得到水煤比，根据总给风量和总给煤量的比值得到风煤比，将燃烧室每个风层的四角风门开度数据取平均值得到每个风层的平均风门开度。之后剔除总给水量、总给风量、每个风层的四角风门开度。保存为初步筛选数据。初步筛选数据形式为X＝{X

初步筛选数据作为本模块的输出，并作为噪声抑制模块的输入。

二、噪声抑制模块

噪声抑制模块作用是将输入数据中存在的噪声进行降噪。实际电厂运行数据中存在各类噪声，噪声可分为由于工业传感器自身波动、由于工业传感器偶然失真产生的数据异常。因此需要对输入数据进行降噪处理。

输入数据有M＝51维特征，每维特征是长度为T的时间序列，任一维特征数据记为X

异常数据的处理方法是：

自定义时间窗口值w(可选w＝100)，从X

噪声的处理方法是：

自定义时间窗口值w(可选w＝100)，从X

通过以上方法，首先将异常数据进行修正，再对高频噪声进行抑制。经过噪声抑制的数据有M维特征，每维特征是长度为依然记为T。处理后的数据进行保存，作为本模块的输出，并作为特征降维模块的输入。

三、特征降维模块

特征降维模块作用是将输入数据进行降维处理，以降低计算难度，节省计算时间。

降维具体方法如下：

步骤一，计算任意两个特征数据的相关系数。当相关系数大于0.95时，剔除这两个特征中方差较小的特征。通过此步骤，可剔除每个磨煤机的电机转速反馈信号、风粉一次风速、风粉混合物温度、混合风量特征，还可以剔除燃烧室上下端部的二次风层风门开度特征。此时数据特征降为k＝19，长度为T，记为X∈R

步骤二，将X∈R

经过特征降维的数据维度为n＝8，长度为T。降维后的数据为本模块的输出，并作为工况聚类算法的输入。

四、工况聚类模块

工况聚类模块作用是将输入数据进行聚类。根据聚类结果得到每个类别的时间戳，即可将锅炉运行数据分为不同工况。

本模块的输入数据维度为n＝8，长度为T，记为X∈R

工况聚类模块可分为参数初始化子模块、聚类赋值子模块以及参数更新子模块。

参数初始化子模块：

参数初始化子模块输入数据为特征降维模块输出数据，作用是将输入数据进行初始化，形成特定形式的数据用于后续计算。

输入数据形式是维度为n，长度为T的多元时间序列，可表示为x

x

其中x

步骤一，需根据实际需要自定义聚类数量K、子序列分割的窗长w、最大迭代次数max_iter、正则项系数β、ρ。

步骤二，将x

X

步骤三，根据自定义聚类数量K，进行随机初始聚类，其结果表示为

步骤四，对每个聚类定义其高斯逆协方差矩阵Θ

其中A

步骤五，对每个Θ

因此，参数初始化模块将输入数据映射为子序列数组X

聚类赋值子模块：

聚类赋值子模块作用是，根据聚类结果P以及参数Θ，得到能够表征聚类结果的值。当表征值的变化小于一定阈值，或迭代次数达到最大迭代次数max_iter，则停止迭代，此时的聚类结果P＝{P

Θ优化问题遵循以下3个条件：1，假定X

其中ll(X

ll(X

因此聚类结果表征值表达式为：

参数更新子模块：

参数更新子模块作用是根据优化目标：

对P＝{P

更新分为两步：1在Θ不变情况下更新P＝{P

更新步骤1：

在Θ不变情况下，优化目标为：

解该优化问题可以等价于：最大化对数似然和聚类在时间上尽可能连续之间的权衡。时间连续性的重要性由参数β控制。

该优化问题可以采用动态规划算法，动态规划算法有向图如图3所示。

在动态规划算法视角下，该最优化问题等价于：求解有向图中从1到T时刻权值最短的路径，图中结点(i，j)的值为当X

因此，通过以上方法，可在Θ不变情况下，采用动态规划算法求解优化问题，对P＝{P

更新步骤2：

在P不变情况下，极大似然项可以改写成：

由于Θ

求解该优化问题，每次都需要至少重复K次，总体算法复杂程度较高。为降低计算复杂度，采用具有并行化计算特性的交替方向乘子算法(ADMM)以解决该优化问题。

采用ADMM算法，引入一致性变量Z以及对偶变量U，则优化问题可改写为：

其对应的增广拉格朗日函数为：

其中ρ＞0，是ADMM算法的惩罚系数。

ADMM算法执行以下三个步骤重复执行直至收敛，停止条件设定为Θ和Z在两次迭代结束时的结果残差小于阈值：

(a)

(b)

(c)U

上述迭代步骤中，步骤(a)和(b)属于Θ和Z优化函数，解优化函数能够得到新的Θ和Z。根据新的Θ和Z，可以通过步骤(c)得到新的U。Θ和Z的更新方法如下所述。

Θ更新：Θ参数的更新与Z和U变量无关，可以写作下式：

该优化问题有一个已知的解析解：

其中QDQ

其中

Z更新：Z参数更新可以写作下式：

定义变量

令

可以得到以下Z的封闭解：

因此，通过更新步骤1和更新步骤2，能够得到新的P＝{P

具体实施方式三：

本领域内的技术人员通过上述实施例提及的系统及方法，本实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式，模块之间也可根据计算机逻辑结构进行重新组织。而且，本实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

根据本实施例的方法、设备(系统)、和计算机程序产品的流程图或方框图来描述的。应理解可由计算机程序指令实现流程图或方框图中的每一流程或方框、以及流程图或方框图中的流程或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程或方框图一个方框或多个方框中指定的功能的步骤。

在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。

存储器可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。存储器是计算机可读介质的示例。

计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。

本实施方式只是对本专利的示例性说明，并不限定它的保护范围，本领域技术人员还可以对其局部进行改变，只要没有超出本专利的精神实质，都在本专利的保护范围内。