CN106599271A - 一种燃煤机组排放监测时序数据异常值检测方法 - Google Patents

Abstract

本发明涉及一种燃煤机组排放监测时序数据异常值检测方法，具体是运用k‑means聚类算法，根据机组运行状态和环保设施运行方式对测点全工况历史样本数据集进行工况划分，然后针对每个工况历史样本集，基于统计假设检验方法，建立异常值检测模型，根据工况样本集的数据分布特征选择不同的假设检验方法对测点时序数据进行在线辨识，及时定位与历史分布和当前变化趋势不一致的异常点。本发明使用假设检验方法对现场监测数据进行异常值在线检测，动态辨识与历史分布或当前变化趋势不一致的异常点，为更为复杂的故障监测提供早期预警，提高了报警工作的精细化程度，提升在线监测管理水平和快速响应能力。

Description

一种燃煤机组排放监测时序数据异常值检测方法

技术领域

本发明属于能源领域，具体涉及一种燃煤机组排放监测时序数据异常值检测方法。

背景技术

随着我国节能减排工作的不断深入开展，燃煤机组污染物排放标准要求越来越严苛。2014年，国家发改委、环保部等部门联合制定了《煤电节能减排升级与改造行动计划(2014-2020年)》，提出了燃煤机组达到燃气轮机组污染物排放标准的要求，即所谓的超低排放，要求在基准氧含量6％条件下，烟尘、二氧化硫、氮氧化物排放浓度分别不高于10、35、50mg/Nm³。

为了响应国家政策的号召，江苏省燃煤机组于2014年开始大力开展超低排放改造，江苏方天电力技术有限公司在政府部门和有关政策的支持下，致力于开展超低排放系统数据在线监控和超低排放电价考核管理工作，将燃煤机组超低排放系统相关数据实施联网集成，其中包括负荷工况参数和烟囱排口SO2、NOx、烟尘浓度、氧量、温度、湿度、压力以及相关环保系统过程辅助参数。

随着超低排放改造工作的不断开展，燃煤机组污染物排放浓度数据的正常与否成为发电企业关注的焦点问题之一。燃煤机组超低排放监测数据随时间变化，是典型的时间序列数据。因此，超低排放时序数据中的异常值有2种类型：1、远离正常值分布范围的孤立点，即与历史分布不一致的异常点；2、由于时间序列数据变化的惯性，异常值会表现为不同于之前的变化趋势，即与当前趋势不一致的异常点。此外，超低排放监测数据以快至10秒、慢至1分钟的时间间隔进行采集，因此一个测点一年的数据记录条数约在52.6～315万之间。因此，超低排放时序数据异常值检测，就是要在这上百万条记录中发现数据正常分布的特征，建立高效的在线检测方法及时发现时序数据变化过程中的异常点，给出报警。

传统的监测数据异常值检测应用最多的是3西格玛法，它是根据测点在正常状态下的历史数据，分别计算出算术平均值μ和标准差σ，然后将正常值的范围设定在区间[μ-3σ,μ+3σ]内，超出这个范围的数据就是异常值。但是在燃煤机组超低排放异常值在线检测场景下，这种固定报警阈值的方法存在以下缺陷：

1)机组和超低排放设施的运行工况和运行方式是不断变化的，在不同的状态下测点会表现出不同的变化特征，进而导致不同状态下的测点数据的均值和方法表现出差异，使得基于全工况样本集得到的固定阈值不能很好地适应持续变化的运行工况；

2)由于固定阈值是在全工况样本集上得到的，为了适应所有工况其取值范围必然相对宽松，这就使得落在阈值范围内、但是其变化趋势异常的数据点无法被有效识别，导出出现漏检情况；

3)3西格玛法是基于数据呈正态分布这个假设的，而燃煤机组超低排放测点的监测数据即使在稳态工况下也不一定服从正态分布，因此3西格玛法不完全适用于燃煤机组超低排放应用。

发明内容

本发明的针对现有技术中的不足，提供一种燃煤机组排放监测时序数据异常值检测方法。

为实现上述目的，本发明采用以下技术方案：

一种燃煤机组排放监测时序数据异常值检测方法，其特征在于，包括以下步骤：

采用k-means聚类算法，对测点全工况历史样本数据集S进行工况划分，将相同工况的数据点归为一类；

针对每个工况样本集S_C，基于统计假设检验方法，建立异常值检测模型；

根据工况样本集S_C的数据特征选择不同的假设检验方法对测点时序数据进行辨识，定位与历史分布和当前变化趋势不一致的异常点。

为优化上述技术方案，采取的具体措施还包括：

所述工况划分包括：

随机选取全工况历史样本数据集S中的k个数据点作为初始聚类中心点；

对于全工况历史样本数据集S中每一个数据点，计算该数据点与各个初始聚类中心点之间的欧式距离，得到距离最近的初始聚类中心点，将这个初始聚类中心点的类作为该数据点的类；

对于每一个聚类，计算该聚类中所有数据点的均值，并将均值作为新聚类中心点；

计算每一个聚类中心点的变化，当最大变化值不超过预定的阈值，终止算法，得到k个聚类中心，表示k个不同的工况划分。

当工况样本集S_C的样本数量>200时，使用KDE核密度估计方法，计算样本集S_C数据分布的概率密度函数，得到95％置信区间，并根据新样本x是否落入95％置信区间来判断是否为异常点。

所述概率密度函数为其中f为概率密度，x₁,x₂...,x_N为独立同分布F的N个样本点，K为满足对称性及∫K(x)dx＝1的核函数，h(h＞0)为平滑参数，σ是样本标准差。

当工况样本集S_C的样本数量≤200时，通过W检验判定工况样本集S_C是否满足正态分布。

将N个独立观察值按从小到大的次序排列，记为x₍₁₎,x₍₂₎...,x_(N)；

计算W检验统计量其中为样本均值，a_i(W)的值根据样本容量N和样本排序下标i从W检验的系数a_i(W)表中得到，

根据显著水平α和样本容量N，基于W检验统计量的p分位数Z_p表得到统计量W的α分位数W_α；

通过比较W和W_α来判断总体是否服从正态分布。

当工况样本集S_C满足正态分布时，根据分布的平均值μ和标准差σ，计算出95％置信区间[μ-1.96σ,μ+1.96σ]，并根据新样本x是否落入95％置信区间来判断是否为异常点。

将容量分别为N₁、N₂的源分布样本和检测分布样本中的所有数据混合，按照大小顺序排列；

计算源分布样本和检测分布样本的秩和，分别记为W₁、W₂；

计算源分布样本和检测分布样本的U检验统计量U₁、U₂：

将U₁、U₂中的最小值作为检验的统计量U＝min{U₁,U₂}；

根据源分布样本和检测分布样本的容量N₁、N₂和置信水平α，基于U检验的临界值表得到统计量U的临界值U_α；

通过比较U和U_α来判断源分布样本和检测分布样本是否来自于同一分布。

本发明的有益效果是：基于统计假设检验的异常值检测方法，对海量运行数据进行工况划分，不同的工况使用不同的样本集，提高了报警工作的精细化程度；这种方法没有使用基于正态分布假设的固定报警阈值，使用了无需先验知识的假设检验方法，提高了方法的通用性；同时，这种方法不仅可以对已有数据集进行异常值检测，也可以用于对现场数据进行异常值在线检测，前者可以用于数据发掘之前的数据与处理过程，后者可以用于工业控制过程，有着广阔的应用前景。

附图说明

图1是本发明的方法流程图。

图2是测试机组全年烟尘排放浓度监测数据。

图3是测试样本集的趋势图。

图4是异常值检测图。

具体实施方式

现在结合附图对本发明作进一步详细的说明。

本发明是采取以下技术方案进行实现的。

一、工况划分

采用k-means聚类算法对测点全工况历史样本数据集S进行工况划分，将相同工况的数据点归为一类。k-means算法的处理过程如下：

1、随机选取数据集S中的k个数据点作为初始聚类中心点；

2、while循环：

1)对于数据集S中每一个数据点，计算其与各个聚类中心点之间的欧式距离，找到距离最近的聚类中心点，将这个中心点的类作为该数据点的类；

2)对于每一个聚类，计算该聚类中所有数据点的均值，并将其作为该聚类新的中心点；

3)计算每一个聚类中心点的变化，如果最大变化值不超过预定的阈值，说明聚类趋于稳定，则终止算法。

在对数据点进行聚类的过程中，一旦聚类数量k确定，k-means算法将在若干次迭代之后得到k个聚类中心，表示k个不同的工况划分。

二、异常值检测方法

1、异常值检测总体思路

基于假设检验的时序数据异常值检测方法处理流程如附图1所示，假设检验过程如下：

1)当工况样本集S_C的样本数量>200时，使用KDE核密度估计方法，计算该样本集数据分布的概率密度函数，进而得到95％的置信区间。如果新样本x落入这个区间内，则判断该点为正常，否则为异常；

2)当样本数量≤200时，由于KDE核密度估计方法要求样本数量要足够多，因此首先通过Shapiro-Wilk检验(W检验)判定工况样本集S_C是否满足正态分布：如果满足正态分布，则利用正态分布的性质，根据分布的平均值μ和标准差σ，计算出95％置信区间[μ-1.96σ,μ+1.96σ]，根据新样本x是否落入该区间来判断其是否为异常点；如果不满足正态分布，则转入下一步；

3)当样本数量≤200，且工况样本集S_C不满足正态分布时，使用Mann-Whitney检验(U检验)来判断新到来的样本x是否来自于工况样本集S_C的数据分布。

2、基于KDE的异常值检测方法

基于KDE的异常值检测方法，首先是使用核密度估计方法，根据已知的历史样本数据求得未知总体的概率密度函数；然后，利用概率密度函数计算95％置信区间，当数值落在置信区间内则认为其正常，否则判断其为异常点。

1)概率密度函数估计

核密度估计(Kernel Density Estimation，KDE)是一类基于概率密度函数的非参数估计法，它从数据样本本身出发研究数据分布的特征，其基本思想是：基于已知的历史数据，对未知总体的概率密度函数进行估计，使估计到的密度函数与真正的密度函数间的均方积分误差最小。

核密度估计方法，就是借助一个移动的单元格(核函数)放在每一个数据点的位置上；然后，将核函数的作用效果叠加起来，获得一条光滑的曲线。而核函数的选择条件为单个峰值下的函数面积为1。

假设x₁,x₂...,x_N为独立同分布F的N个样本点，设其概率密度为f，则其核密度函数估计为：

公式(1)中的K被称为核函数，通常满足对称性及∫K(x)dx＝1，核函数是一种加权函数，利用数据点x_i到x的距离(x-x_i)来决定x_i在估计点x的密度时所起的作用，一般选择标准正态函数作为核函数，如公式(2)所示，则离x点越近的样本点其加权也越大。

公式(1)中的h(h＞0)是一个平滑参数，称为带宽(band width)，带宽值的选择对估计f(x)的影响很大。当h很小时，只有特别接近x的点才起较大作用；随着h的增大，则距离x远一些点的作用也随之增加。标准正态核函数的带宽h可由Silverman拇指法则得到：

式中σ是样本标准差。

2)95％置信区间估计

对于一个给定的数据分布(概率密度函数已知)，95％置信度区间是指服从该分布的任意随机变量将有95％的可能性落在该区间内。对于具有对称性的概率密度函数，则是概率为2.5％至97.5％之间的样本数值。

特别地，对于服从N(μ,σ²)正态分布的样本，其95％置信区间[μ-1.96σ,μ+1.96σ]，其中μ、σ分别为样本均值和标准差。

3、正态性检验-Shapiro-Wilk检验(W检验)

为了处理小样本情况，本方法采用Shapiro-Wilk检验，即W检验，来确定工况样本是否服从正态分布。在满足正态性(通过W检验)的前提下，直接通过正态分布的95％置信区间估计方法得到正常数据的取值范围，不在这个取值范围内的数值就被认为是异常点。

检验步骤如下：

1)将N个独立观察值按从小到大的次序排列，记为：x₍₁₎,x₍₂₎...,x_(N)；

2)依次将x₍₁₎,x₍₂₎...,x_(N)的值带入公式(4)计算W检验统计量：

式中为样本均值；a_i(W)的值根据样本容量(等同于独立观察值)N和样本排序下标i从W检验的系数a_i(W)表中得到；[N/2]根据公式(5)计算：

3)根据预先设定的显著水平α和样本容量N，查W检验统计量的p分位数Z_p表得到统计量W的α分位数W_α；如果W＜W_α，则拒绝原假设H₀，认为总体不服从正态分布；否则，接受假设H₀，认为总体服从正态分布。

4、同分布检验-Mann-Whitney检验(U检验)

对于不服从正态分布的样本，可以采用Mann-Whitney检验(U检验)来判断两组样本数据是否来自于同一分布。在异常值在线检测应用场景中，将工况历史样本作为源分布，将在工况历史样本集中加入待检测的样本点后形成的样本作为检测分布，则如果以不低于95％的概率来源于同一分布的话，就认为待检测点为正常点，将其历史样本集后并没有改变源数据分布。

其检验步骤如下：

1)将容量分别为N₁、N₂的两组样本中所有数据混合，然后按照大小顺序排列；最小值的秩为1，第二小的值秩为2，以此类推。如果有数值相等的情形，则取这几个数值排序下标的平均值作为其秩，如：样本(3,5,5,9)中，数值5的秩为(2+3)/2＝2.5，所有样本的秩为(1,2.5,2.5,4)；

2)根据上一步中得到的秩，分别求出两个样本的秩和，即样本集中所有数值的秩之和，记为W₁、W₂；

3)按照公式(6)分别计算两个样本的U检验统计量U₁、U₂：

4)将U₁、U₂中的最小值作为本次检验的统计量U＝min{U₁,U₂}；

5)根据两个样本的容量N₁、N₂和置信水平α，查U检验的临界值表得到统计量的U的临界值U_α；

6)比较U和U_α，如果U＜U_α，则拒绝原假设H₀，认为两组样本来源于不同分布；否则，接受原假设H₀，认为两组样本分布没有差异。

以下结合算法的实施过程对本发明作具体的介绍：

以一台装机容量为600MW的燃煤机组(以下称为测试机组)，2015年7月至2016年6月1年的除尘超低排放设施运行数据为样本建立模型，以2016年7月1日数据为测试样本，以烟尘排放浓度为目标测点，检测其时序数据中的异常值。

测试机组的除尘超低排放设施采用脱硫前干式除尘+脱硫后湿式除尘工艺，其中干式除尘器有4个电室4个电场，湿式除尘器有8个电场。采集测点包括：机组负荷、所有电场的二次电流和二次电压、排口的烟尘排放浓度，其中机组负荷代表机组运行工况，决定了除尘设施入口烟气参数；所有电场的运行状态(运行、停运)的不同组合代表了除尘设施的运行方式。

按照本方法的处理流程，首先要使用聚类方法对测试机组烟尘排放浓度全工况历史样本集进行工况划分。由于烟尘排放浓度受机组负荷和除尘设施运行方式的影响，因此选择这些测点作为工况划分参数，对历史样本及进行分组，将相似工况下的排放浓度数据放在一起形成工况样本集。为了便于展示，表1中仅列出了以机组负荷为工况参数进行聚类得到的工况统计数据。

表1测试机组历史样本集工况划分结果

由于进行了超低排放改造，烟尘排放浓度降至烟尘超低排放标准10mg/m³以下。因此，为了不影响样本均值和标准差的计算结果，将烟尘排放浓度数值>10的样本剔除。测试机组全年烟尘排放浓度监测数据如图2所示。

以烟尘排放浓度全年监测数据为历史样本集，其均值μ＝4.08、标准差σ＝0.41，则根据3西格玛法其95％置信区间为[2.85,5.31]。图3是测试样本集的趋势图，图中用横线标注出了95％置信区间的上下限，从图中可以看出所有数据都落在了置信区间内，本方法没有发现任何异常值。

根据工况划分的结果，通过计算每个测试样本到每个工况聚类中心(用机组负荷平均值、所有电场运行状态组成的向量表示)的距离，然后将距离最近的聚类作为该数据点所属工况，并以该工况的样本作为源样本，通过假设检验逻辑来判断数据点是否异常。

使用本方法在测试样本集上进行测试得到的异常值检测，结果如图4所示，异常值分别出现在样本序号50和256处，图中用圆圈标出。这两处异常都是趋势图中出现突变的地方，使得之前的变化趋势被破坏，但是其取值又没有超出3西格玛法的95％置信区间[2.85,5.31]。由此可见，本方法对与当前趋势不一致的异常点比较敏感，可以高效地从时间序列数据中发现异常值。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (8)

1.一种燃煤机组排放监测时序数据异常值检测方法，其特征在于，包括以下步骤：

采用k-means聚类算法，对测点全工况历史样本数据集S进行工况划分，将相同工况的数据点归为一类；

针对每个工况样本集S_C，基于统计假设检验方法，建立异常值检测模型；

根据工况样本集S_C的数据特征选择不同的假设检验方法对测点时序数据进行辨识，定位与历史分布和当前变化趋势不一致的异常点。

2.如权利要求1所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：所述工况划分包括：

随机选取全工况历史样本数据集S中的k个数据点作为初始聚类中心点；

对于全工况历史样本数据集S中每一个数据点，计算该数据点与各个初始聚类中心点之间的欧式距离，得到距离最近的初始聚类中心点，将这个初始聚类中心点的类作为该数据点的类；

对于每一个聚类，计算该聚类中所有数据点的均值，并将均值作为新聚类中心点；计算每一个聚类中心点的变化，当最大变化值不超过预定的阈值，终止算法，得到k个聚类中心，表示k个不同的工况划分。

3.如权利要求1所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：当工况样本集S_C的样本数量>200时，使用KDE核密度估计方法，计算样本集S_C数据分布的概率密度函数，得到95％置信区间，并根据新样本x是否落入95％置信区间来判断是否为异常点。

4.如权利要求3所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：所述概率密度函数为其中f为概率密度，x₁,x₂...,x_N为独立同分布F的N个样本点，K为满足对称性及∫K(x)dx＝1的核函数，h(h＞0)为平滑参数，σ是样本标准差。

5.如权利要求1所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：当工况样本集S_C的样本数量≤200时，通过W检验判定工况样本集S_C是否满足正态分布。

6.如权利要求5所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：

将N个独立观察值按从小到大的次序排列，记为x₍₁₎,x₍₂₎…,x_(N)；

计算W检验统计量其中为样本均值，a_i(W)的值根据样本容量N和样本排序下标i从W检验的系数a_i(W)表中得到，

根据显著水平α和样本容量N，基于W检验统计量的p分位数Z_p表得到统计量W的α分位数W_α；

通过比较W和W_α来判断总体是否服从正态分布。

7.如权利要求5所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：当工况样本集S_C满足正态分布时，根据分布的平均值μ和标准差σ，计算出95％置信区间[μ-1.96σ,μ+1.96σ]，并根据新样本x是否落入95％置信区间来判断是否为异常点。

8.如权利要求5所述的一种燃煤机组排放监测时序数据异常值检测方法，其特征在于：

将容量分别为N₁、N₂的源分布样本和检测分布样本中的所有数据混合，按照大小顺序排列；

计算源分布样本和检测分布样本的秩和，分别记为W₁、W₂；

计算源分布样本和检测分布样本的U检验统计量U₁、U₂：

$U_1=N_1{N}_2+\frac{N_1(N_1+1)}{2}-W_1$

$U_2=N_1{N}_2+\frac{N_2(N_2+1)}{2}-W_2;$

将U₁、U₂中的最小值作为检验的统计量U＝min{U₁,U₂}；

根据源分布样本和检测分布样本的容量N₁、N₂和置信水平α，基于U检验的临界值表得到统计量U的临界值U_α；

通过比较U和U_α来判断源分布样本和检测分布样本是否来自于同一分布。