CN109034252B - 空气质量站点监测数据异常的自动化识别方法 - Google Patents

Abstract

一种空气质量站点监测数据异常的自动化识别方法，包括：步骤A：接收站点监测数据；步骤B：对所接收的站点监测数据进行初级检查，识别显著异常的观测数据；步骤C：对完成初级检查后的数据进行时空一致性检查，识别时空不一致异常；步骤D：步骤C后，使用四项补充检查识别未被初级检查和时空一致性检查识别的异常；步骤E：输出经过质控后的观测数据、时空一致性估计值、及各项检查中的概率值；所述方法用以缓解现有技术自动化异常识别方法中难以识别特有的周期性异常、延滞异常、以及对于正定的，更接近于对数正态的空气质量监测数据(PM_2.5，PM₁₀，SO₂，NO₂，CO和O₃)的异常识别效果较差，难以识别数值较低或观测误差小于观测标准差的异常数据等技术问题。

Description

空气质量站点监测数据异常的自动化识别方法

技术领域

本公开涉及大气污染领域，尤其涉及一种空气质量站点监测数据异常的自动化识别方法。

背景技术

大范围、准确的常规污染物观测数据是衡量空气质量的重要依据和相关研究的基础。然而，由于仪器故障、恶劣环境、以及监测方法的局限异常观测数据出现不可避免。在实际应用过程中，通常需要人工对监测数据进行审核和质控，以剔除异常的观测数据。这种方法通常能有效剔除大气污染监测中异常数据。其主要缺点是非常繁杂，需要耗费大量人力和时间，很难快速获得大量的质控数据，制约了数据的快速应用。此外，不同人的质控标准具有一定主观性，难以完全一致，从而可能给质控数据集引入一定的偏差。因此，有必要建立一种客观的、具有统一标准的质控技术方法。

气象观测数据的自动化质控较为成熟，已成为各气象数据集归档时不可或缺的一部分，在质控中，会根据观测变量特有的异常特征针对性地设计算法识别风向、风速、降雨、降雪等变量的异常观测数据；同时，会依据观测数据在时间和空间上的一致性判断数据合理性；在海洋观测(温度，深度，盐度)和土壤观测(温度，湿度)中，质控也有一些较为成熟的研究，但在大气污染领域，美国国家环境保护局、欧洲环境署、英国环境部都制定了观测质控的规范手册，但其核心是观测操作规范和仪器维护，大范围的大气污染监测数据自动化质控研究仍然非常少，在实际应用过程中，通常以人工审核的方式，或借用其他领域的通用自动化质控方法识别异常监测数据。

现有的主流空气质量监测数据的异常识别方法有两种，一种是以人工审核的方式识别异常数据，另一种借用其他领域的通用自动化质控方法识别异常监测数据。人工审核的方法依赖质控员的经验，通过人眼从监测的时间序列或空间分布等其他图件或表格中，挑出可以的观测数据。对可疑的监测数据可以通过组织相关人员进行站点周边的实地考核，进一步确认该监测数据是否异常。通用的质控方法可参考气象数据异常识别中常用的z-score方法，该方法分三步。第一步计算监测值平均，第二步计算监测值的标准差，第三部将偏离均值几倍标准差的数据标记为异常，现有的两种异常识别方法都存在各自的缺陷，其中人工审核效率低下，难以适用于实时或大规模的监测数据应用，例如在线监测数据发布、将在线监测数据同化进空气质量预报系统以改进预报效果、以及构建全国多年的空气质量再分析场，于此同时，人工审核的方法缺少原理支撑，其结果易受质控员的主观经验的影响；另一种异常识别方法，即通用的自动化异常识别方法未针对中国环境空气质量监测网设计，难以识别其特有的周期性异常、延滞异常(异常定义在本公开技术方案中介绍)，并且通用方法常隐含监测数据正态分布的假设，对于正定的，更接近于对数正态的空气质量监测数据(PM_2.5，PM₁₀，SO₂，NO₂，CO和O₃)的异常识别效果较差，难以识别数值较低的异常数据。因此，迫切需要发展一种针对大气环境监测网络的常规大气污染物监测数据的典型异常类型的自动化质控新方法。

公开内容

(一)要解决的技术问题

本公开提供了一种空气质量站点监测数据异常的自动化识别方法，以缓解现有技术自动化异常识别方法中难以识别其特有的周期性异常、延滞异常、以及对于正定的，更接近于对数正态的空气质量监测数据(PM_2.5，PM₁₀，SO₂，NO₂，CO和O₃)的异常识别效果较差，难以识别数值较低或观测误差小于观测标准差的异常数据等技术问题。

(二)技术方案

本公开提供一种空气质量站点监测数据异常的自动化识别方法，包括：步骤A：接收站点监测数据；步骤B：对步骤A所接收的站点监测数据进行初级检查，识别显著异常的观测数据；步骤C：对步骤B完成初级检查后的数据进行时空一致性检查，识别时空不一致异常；步骤D：经步骤C后，使用四项补充检查识别未被初级检查和时空一致性检查识别的异常；以及步骤E：输出经过质控后的观测数据、时空一致性估计值、及各项检查中的概率值。

在本公开实施例中，所述步骤B中的初级检查，包括：步骤B₁：完整性检查；步骤B₂：超量程检查，对监测数据进行上下限检查，将超出仪器量程的错误记录剔除；以及步骤B₃：大观测误差检查，剔除超出合理值很多的观测，以减弱其对时空连续性检查性能的影响。

在本公开实施例中，步骤C中所述时空一致性检查，包括：时间一致性回归；以及空间一致性回归。

在本公开实施例中，所述时间一致性回归，利用检验点邻近时刻的观测数据，计算检验点的时间回归值，回归方法采用低通滤波，即：

其中F_t为滤波估计值，i为检验点的时次，k代表滤波时间窗口从检验点往前和往后的时间长度，f为原始观测，h为滤波系数。

在本公开实施例中，所述空间一致性回归，是结合邻近空间范围内的观测值计算得到检验点的估计值，具体计算公式如下：

其中F_s(i)为目标站点在i时刻的空间一致性估计值。f_r为第r个参考站点的观测值。a_r为检验站点与参考站点间的一致性指标，采用以下方法进行计算：

其中f_r(i+k)为参考站点在i+k时刻的观测值，

为滑动窗口内的观测平均值。

在本公开实施例中，所述空间一致性回归，权重c_r采用Gaspari-Cohn(高斯-康恩)方案计算：

其中d为目标站点与参考站点之间的距离，d_c为截止距离。

在本公开实施例中，根据所述时间和空间的一致性估计值F_t和F_s，计算检验点的归一化估计残差Z_t和Z_s，再计算残差相关系数：

进而计算残差概率：

其中i为目标时次，ρ为时空残差的相关系数，Z_t，Z_s分别为归一后的时间和空间的回归残差，

分别为滑动窗口内时间和空间的归一残差平均，i-n和i+n分别为滑动窗口起始和结束时间。

在本公开实施例中，步骤D所述的四项补充检查，包括：

小变化异常检查，观测值呈现出长时间常值或过于缓慢的异常时段，与实际大气污染变化特征不吻合，所述异常时段数据剔除；

周期性异常检查，识别周期出现的异常并进行剔除；

PM₁₀＜PM_2.5异常检查，当PM_2.5与PM₁₀浓度出现倒挂时，将PM₁₀观测数据剔除；以及

有效数据量检查，即对每个观测数据，统计其前后12小时内的有效数据，若有效数据少于5个则对其进行剔除。

在本公开实施例中，所述小变化异常检查中异常时段的残差概率P_a计算如下式：

其中R_a、S_a、Z_a为延滞时段的回归残差、回归残差标准差和归一化的回归残差，R_s、S_s为上面内容中计算的空间回归残差及其标准差，b和e分别为延滞时段的开始和结束时次，计算得到小变化异常残差概率P_a，残差概率P_a小于设定阈值时的观测值识别为异常并剔除，所述阈值为10^-3～10^-9。

在本公开实施例中，所述周期性异常检查，首先，以24小时为间隔，对原始观测f进行滑动平均计算，如下：

其中i为待检查时次，然后，对f_c进行中值滤波，如下：

F_c(i)＝M(f_c(i+k)，k∈[-1，1])

其中M为集合的中位数，计算得到回归值F_c，滑动窗口长度取3个时次，接着，通过f_c和F_c计算得到回归残差R_c，并以94百分位的回归残差作为回归残差的标准差，使得得到的标准差σ_c大于一天中第二大的回归残差，公式如下：

S_c(i)＝g(R_c(i+k)，k∈[-72，72])

其中g为集合的94百分位，最后，将R_c和S_c一起代入下列公式：

计算得到周期性异常残差概率P_c，残差概率P_c小于阈值的观测识别为异常并进行剔除，所述阈值为10^-2～10^-4。

(三)有益效果

从上述技术方案可以看出，本公开空气质量站点监测数据异常的自动化识别方法至少具有以下有益效果其中之一或其中一部分：

(1)可以在一分钟内识别十万个观测数据中的异常数据。

(2)结果具有确定性，不受质控员主观因素的影响。

(3)可以应用于实时或大规模的监测数据应用。

(4)能进行在线监测数据发布、将在线监测数据同化进空气质量预报系统以改进预报效果、以及构建全国多年的空气质量再分析场。

(5)通过基于时空连续性等污染监测数据特征设计正常观测数据的回归算法，可以对每个观测数据计算其概率，从而定量地评估每个观测数据的合理性，效果完全取决于算法中的参数设置和观测数据，其结果具有确定性，不受质控员主观因素的影响。

(6)针对中国环境空气质量监测网的异常数据特征设计，可以识别其特有的监测数据异常。

(7)设计正常观测的回归算法，使异常识别方法更适用于正定的，更接近于对数正态分布的空气质量监测数据，其置信区间更小，能更好地识别浓度较低的异常数据。

(8)基于概率理论，可将多项检查有机结合。

附图说明

图1为本公开实施例空气质量站点监测数据的异常情况分类示意图。

图2为本公开实施例空气质量站点监测数异常的自动化识别方法的流程框架示意图。

图3为本公开实施例空气质量站点监测数据异常的自动化识别方法的具体流程示意图。

图4为本公开实施例超量程异常数据示意图。

图5为本公开实施例中国环境空气质量国控站所用仪器的量程检测范围示意图。

图6为本公开实施例滤波系数的计算结果示意图。

图7为本公开实施例时间一致性检查和时空一致性检查的质控效果对比示意图。

图8为本公开实施例周期性异常检查的质控效果示意图。

图9为本公开实施例z-score方法与基于回归残差概率方法的监测数据异常识别效果对比示意图。

具体实施方式

本公开提供了一种空气质量站点监测数据异常的自动化识别方法，通过将大气污染监测数据异常分类后，利用计算机的快速运算性能可以实时或大规模的监测数据应用，实现监测数据异常的自动化识别，以缓解现有技术中即自动化异常识别方法中难以识别其特有的周期性异常、延滞异常、以及对于正定的，更接近于对数正态的空气质量监测数据(PM_2.5，PM₁₀，SO₂，NO₂，CO和O₃)的异常识别效果较差，难以识别数值较低或观测误差小于观测标准差的异常数据等技术问题。

为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。

图1为空气质量监测数据的异常情况分类示意图，如图1所示，所述空气质量监测数据异常情况分类包括：

(1)数据在时空上的一致性异常(ST异常)，大气污染的平流、扩散等过程使其浓度在时间和空间上呈现出连续性，当某一观测点数据出现足够大的误差时，结合图1(a)和图1(b)所示，其观测值会与邻近时次和周边站点出现较大差异，表现出离群特征，这种异常我们分别称之为数据的时间一致性异常和空间一致性异常。

(2)数据的长时间小变化异常(LV异常)，结合图1(c)和图1(d)所示，大气污染观测数据另外一种异常使观测值呈现出长时间常值或过于缓慢的变化，与实际大气污染变化特征不吻合，导致这种异常的原因可能是观测仪器的采样泵卡死，纸带耗尽等，此外，CO监测仪器依据两个吸收室的压力差计算浓度，当两个吸收室红外辐射源的衰老不一致时，会使观测数据产生持续增强的漂移。

(3)周期性异常(P异常)，如图1(e)所示，大气污染监测数据中还存在周期性出现的异常数据，这类异常数据通常连续地在每天的固定时刻出现，大气污染物监测中，由于仪器发光元件老化、环境变化等影响，仪器测量值会产生漂移，需要对仪器定期进行校准。

(4)PM_2.5和PM₁₀观测值的“倒挂”异常(LP异常)，如图1(f)所示，实际观测的PM_2.5测量值可能大于PM₁₀的测量，这与理论不相符，导致这种异常的主要原因是PM₁₀和PM_2.5监测方法的差异，中国很多城市在2012年以前只开展PM₁₀业务化监测，而没有PM_2.5监测，当时PM₁₀监测主要是采用恒温加热的β射线法和震荡天平法进行测量，其中恒温加热的β射线法的不足之处在于：环境温度过低时可能将半挥发性有机物过度蒸发；震荡天平法的不足之处在于：加热采样气体进行除湿时，可能将半挥发性有机物过度蒸发，2012年以后，中国开始对PM_2.5进行全面监测，对方法进行了改进，采用动态加热的β射线法或联用膜动态补偿的震荡天平法，可以防止挥发性有机物过度挥发，或对其进行补偿，因此，当同一站点对PM₁₀和PM_2.5监测仪器的方法原理不同时，可能会出现PM_2.5和PM₁₀观测值的“倒挂”异常。

在本公开实施例中，提供一种空气质量站点监测数据异常的自动化识别方法，图2为所述方法的流程框架示意图，图3为所述方法的具体流程示意图。时空一致性异常是空气质量监测异常数据的主要特征，本方法以时空一致性检查为核心，为保障时空一致性检查的效果，在其之前设计初级检查以剔除显著异常的干扰数据；在其之后，针对尚未被识别的异常数据，根据其特征，设计四项其他检查作为补充。结合图2和图3所示，所述空气质量站点监测数据异常的自动化识别方法，包括：

步骤A：接收站点监测数据；

步骤B：对步骤A所接收的站点监测数据进行初级检查，识别显著异常的观测数据；

步骤C：对步骤B完成初级检查后的数据进行时空一致性检查，识别时空不一致异常；

步骤D：经步骤C后，使用四项补充检查识别未被初级检查和时空一致性检查识别的异常；以及

步骤E：输出经过质控后的观测数据、时空一致性估计值、及各项检查中的概率值。

在本公开实施例中，所述步骤E中时空一致性估计值，用于缺测值的补充；所述各项检查中的概率值，为监测数据可靠性的量化指标，便于与人工检查相结合，步骤E后，即完成了空气质量站点监测数据异常的自动化识别。

在本公开实施例中，所述步骤B中初级检查，包括：

步骤B₁：完整性检查，以保障获得的观测记录在传输存储过程中保持完整可靠，完整性检查最先进行，在数据读取中就加以考虑；

步骤B₂：超量程检查，完成完整性检查后，对监测数据进行上下限检查，将超出仪器量程的错误记录剔除；以及

步骤B₃：大观测误差检查，剔除超出合理值很多的观测，以减弱其对时空连续性检查性能的影响。

中国环境空气质量国控站对全国各地的空气质量进行联网实时监控，各地检测站的观测数据直接传送到中国环境监测总站，对于这类观测与数据收集处理不在一起的远程观测，因数据传输引起的错误需不可避免(Gandin，1988)，在传输、存储和读取中，记录中可能出现缺失、重复、乱码、断行等问题，针对这一问题，质控系统在进行其他检查前，首先对观测数据记录的完整性进行检查，这项检查中，只有当某条记录拥有确定的字符长度，在特定位置满足特定字符类型，并且通过冗余检查的数据才在第一轮检测中确认其有效，第二轮检测中对已经识别的、可修复的异常记录进行修复。这类异常有多余空白、记录断行、空白代替数字，第三轮检测中对记录的唯一性进行检查，确保每个时次每个站点只有一条完整记录。

图4为本公开实施例超量程异常数据示意图，图5为本公开实施例中国环境空气质量国控站所用仪器的量程检测范围示意图。如图4所示，监测数据中有时会出现超出仪器量程的不合理数值，可以完全确定地判断其为异常数据，这类异常可能来源于仪器参数的不合理设置，也可能是恶劣天气下仪器来不及达到平衡导致比如红外光源不稳定会导致CO的观测产生异常高值和负值，暴雨后采样膜的水汽蒸发量可能高于颗粒物累积量使观测浓度出现负值。这类异常可能单独出现，也可能连续出现，在连续出现时若只看起相对大小可能难以判断其为异常，如图4所示，其中大部分时次的观测大于10000的仪器观测上限，故而有必要对其进行单独检查，质控系统在完成数据完整性检查后对这部分异常数据进行识别，仪器量程参考中国环境空气质量国控站所用仪器的检测范围(图5所示)，将超出各污染物仪器量程的数据剔除，注：量程范围来自国标(HJ653-2013，HJ654-2013)，颗粒物监测有两套量程对应两种监测方法。

Z-score方法是识别气象观测数据异常最常用的方法，其计算方法如公式(1)所示，它利用观测数据的均值和标准差，将原始观测归一化得到Z，最后将Z超过设定阈值的观测识别为异常数据。

式中f为原始观测，

为观测平均，σ为观测标准差，i为检验点时次。

Z-score方法能识别显著偏离大部分观测的异常数据，但它不适用于显著偏离正态分布的观测，大气污染物非负的浓度数据更接近对数正态分布，直接使用Z-score方法可能将部分高浓度的正常观测数据识别为异常数据。针对这个问题，同时为了增强异常识别能力，我们对Z-score方法做如下修改。

使用更符合正态分布的回归残差代替原始观测公式，将回归残差的Z值作为评判依据以检验观测的合理性。式中F为回归值，R为回归残差。每项异常检查中，根据检查的内容设计相应的回归方法计算得到F。回归方法设计时，竭力使异常数据的回归残差最大程度地大于正常数据的回归残差，并且使正常数据的回归残差接近均值为0的正态分布。

R(i)＝f(i)-F(i) (2)

使用滑动窗口代替整个研究时段计算单个时刻的标准差，如公式(3)，i-n和i+n分别为滑动窗口起始和结束时间，S为回归残差的标准差，i+k为滑动窗口计算时的采样点，通过使用滑动窗口，可以使回归残差的标准差S随观测浓度的变化规律自适应地改变。

基于上述两点修改，用回归残差R及其标准差S代替原始观测f和标准差σ，故将公式1改为为公式(4)。

由于理论上回归残差的均值

为0，且后续的概率计算不受正负的影响，故而式中分子可由

简化为R(i)。

进一步计算Z的概率。使用公式(5)得到的Z服从均值为0，标准差为1的正态分布，其概率P的计算公式列于式公式(6)，使用概率的优点是可以依据概率理论更好地将多项检查结合起来。

修改后的回归残差概率方法(PRR，Probability of Regression Residual)将应用于识别监测数据的异常情况。

在本公开实施例中，PRR方法使用观测点(检验点)邻近时空的观测值，计算得到观测点的估计值，然后根据观测值和估计值的差异大小对异常数据进行判断和识别，这个过程中，估计值的精度与邻近时空观测值的精度直接相关，如果其中包含很大观测误差的异常数据，可能使估计值产生偏差，进而影响数据质控的效果，因此，为了保障数据质控的效果，本公开其他使用PRR方法的检查之前，对显著异常的观测数据进行了预剔除，其方法如下：

使用中值滤波对作为PRR中的回归方法，对测量值进行估计，中值滤波的公式如下：

F_m(i)＝M(f(i+k)，k∈[-n，n]) (7)

其中，F_m(i)为i时刻的滤波估计值，M为中值函数，f为原始观测，i-n和i+n分别为滑动窗口起始和结束时间。

S_m(i)＝1.4826M(|R_m(i+k)|，k∈[-n，n]) (8)

其中S_m是使用估计残差绝对值的中位数间接估计，R_m为将F_m带入公式(2)所得的残差。

与直接使用原始数据相比，这种方法估计的标准差更为稳健，不容易受大误差数据的影响。最后利用公式(5)和(6)，计算得到中值滤波估计值的残差概率P_m。通过敏感性测试，将10^-15作为概率限制，将P_m小于10^-15的检验点作为显著异常数据进行剔除。

观测数据在时空上的离群异常是观测中最常见的异常，其呈现出来的特征也最复杂，需要结合一定时空范围内的观测信息才能有效识别。针对这一数据异常问题，本公开中介绍的基于概率判断的异常识别方法，结合检验点邻近时空范围内的观测值，计算二元正态分布假设下时空回归估计值的残差概率，再根据概率对异常数据进行判断和识别。

在本公开实施例中，步骤C中所述时空一致性检查，包括：

时间一致性回归；以及

空间一致性回归。

第一步利用检验点邻近时刻的观测数据，计算检验点的时间回归值。本回归方法采用低通滤波：

其中F_t为滤波估计值，i为检验点的时次，k代表滤波时间窗口从检验点往前和往后的时间长度，f为原始观测，h为滤波系数，滤波系数由通行频率，截止频率，和计算方法决定。本公开的通行和截止频率分别为1/8和1/24小时，使用equiripple FIR filter(等波纹有限脉冲响应滤波器)作为滤波系数的计算方法，滤波系数的计算结果如图6所示，(注：滤波系数关于时移对称分布，表中只列出非负部分)，使用低通滤波可以有效的抑制原始数据的瞬时变化，并保留日以上周期变化，相比滑动平均，它赋予临近时次观测更大的权重，对正常观测数据有更小的回归残差，可以更好地区分正常和异常观测，相比中值滤波，其估计的残差更接近正态分布。

本公开实施例质控方法不仅利用检验点邻近时间窗口的观测信息来判断，同时还利用检验点邻近空间范围内的观测信息来判断。因此，第二步是结合邻近空间范围内的观测值计算得到检验点的估计值，具体计算公式如下：

其中F_s(i)为目标站点在i时刻的空间一致性估计值。f_r为第r个参考站点的观测值。a_r为检验站点与参考站点间的一致性指标，采用以下方法进行计算：

其中f_r(i+k)为参考站点在i+k时刻的观测值，

为滑动窗口内的观测平均值。一致性指标常用于评估模拟效果，在质控中也有应用。相比于相关系数，它受奇异值的影响更小，能更好地评估站点间观测的一致性。其不足之处在于，对于两组完全不相关的序列，计算得到的一致性系数也不为0，这使得估计结果会受不相关站点观测值的干扰。针对这一问题，我们借鉴同化中局地化的思路，减小距检验站点较远的参考站点权重，权重c_r采用Gaspari-Cohn(高斯-康恩)方案计算：

其中d为目标站点与参考站点之间的距离，d_c为截止距离。

通过前述公式(9)和(10)分别获得时间和空间的一致性估计值F_t和F_s后，估计残差符合二元正态分布。首先利用公式(2)-(5)计算检验点的归一化估计残差Z_t和Z_s，再计算残差相关系数：

进而计算残差概率：

其中i为目标时次，ρ为时空残差的相关系数，Z_t，Z_s分别为归一化后的时间和空间的回归残差。

分别为滑动窗口内时间和空间的归一残差平均。i-n和i+n分别为滑动窗口起始和结束时间。概率过低的检验点观测值与临近时空范围观测值有显著差异，具备“数据在时空上的离群异常”的特征，识别为异常观测数据。这一方法的优点在于通过二元正态分布，可以将一定时间和空间内的观测信息协同评估。相比于只考虑单个站点在一定时间范围内的观测信息或只考虑某一个时次周边站点的观测信息，该方法更易识别在时空上离群的异常观测数据。

图7是时间一致性检查和时空一致性检查的质控效果对比。其中浅色阴影区域为时间一致性置信区间，区间外的观测为时间不一致性异常，深色阴影表示时空一致性置信区间。图中观测1、2、3处与周边时次和周边站点有明显差异，判断其为时空不一致异常。观测4、5处虽然与周边时次差异较大，但与周边站点的变化保持一致，判断其为正常观测。时间一致性检查虽然能有效识别1、2，但在未考虑空间一致性的情况下，异常3没有被有效剔除，而正常观测4、5被错误地识别为异常。而在时空一致性检查中，通过综合考虑前后时次和周边站点，置信区间在大部分区域有很好的压缩，能更好的识别跟周围站点和时次变化不一致的异常观测。

在本公开实施例中，步骤D中所述的四项补充检查，包括：小变化异常检查。

与数据在时空上的“离群”特征不同，这类异常观测值呈现出长时间常值或过于缓慢的变化，单个观测值往往在正常测量值范围内，与周边时次的观测值也较为接近。采用常规的时空一致性检查方法往往难以直接识别和剔除。针对这类异常，本项检查先将小变化的时段识别出来，再对其合理性进行判断：

首先，通过观测值随时间变化的一阶、二阶导数，识别出连续相同或变化很小的观测时段。在一些比较清洁和稳定的时段，由于观测仪器的分辨率有限，识别出的部分小变化时段的测量值可能为正常观测。故而本项检查结合空间连续性，进一步判断小变化时段的合理性。

其次，将识别出的时段作为整体，计算延滞时段的归一化残差，(公式15)，时段的回归残差为空间回归残差的平均，时段的标准差也依据正态分布中单样本与均值间标准差的关系相应改变。将归一化残差Z_a带入公式(6)得到延滞时段的残差概率Pa，残差概率小于10^-6的时段识别为异常并进行剔除。

其中R_a、S_a、Z_a为延滞时段的回归残差、回归残差标准差和归一化的回归残差。R_s、S_s为上面内容中计算的空间回归残差及其标准差，b和e分别为延滞时段的开始和结束时次。

本检查能有效识别出观测值变化小且空间一致性差的异常观测。相比于剔除连续相同值的方法，它能识别出不严格相同的缓慢变化时段，同时也能保留变化虽小但空间一致性高的合理观测。

在本公开实施例中，步骤D中所述的四项补充检查还包括：周期性异常检查。

大气污染观测数据存在另外一种异常，这类异常通常间隔24小时周期性出现，针对这一特征，我们根据观测数据的日变化规律对这类异常进行识别。

首先，以24小时为间隔，对原始观测f进行滑动平均计算，由公式(16)，得到新序列f_c作为本项检查的“观测值”。新的观测序列为临近五天的日变化平均，可以增强观测的日变化，使得周期出现的异常更为突出。

其中i为待检查时次。

然后，对f_c进行中值滤波，由公式(17)，得到回归值F_c。因为自校准异常单独出现，故而滑动窗口长度取3个时次。

F_c(i)＝M(f_c(i+k)，k∈[-1，1]) (7)

其中M表示集合的中位数。

接着，由f_c和F_c计算得到回归残差R_c，并以94百分位的回归残差作为回归残差的标准差，公式(18)，使得得到的标准差σ_c大于一天中第二大的回归残差(第二大残差处于93.75百分位)，从而确保在本项检查中一天最多剔除一个时次。

S_c(i)＝g(R_c(i+k)，k∈[-72，72]) (8)

其中g为集合的94百分位。

最后，将R_c和S_c一起代入公式(5)和(6)，得到概率P_c，概率小于10^-4的观测识别为异常并进行剔除。

图8给出第三类异常即周期性异常检查的质控效果示意图。以武汉汉口花桥站的数据监测数据为例，如图所示，在每天凌晨4时臭氧观测值突然跃升，与其日变化规律不符。部分第三类异常的时空一致性较差，可以在时空一致性检查中被识别出来，另一部分异常的波动幅度小于日变化波动幅度，需要通过本项检查加以识别。

在本公开实施例中，步骤C中所述的四项补充检查还包括：PM₁₀＜PM_2.5异常检查。

同一站点同一时次观测的PM_2.5浓度值大于PM₁₀浓度值。因为中国绝大部分PM_2.5监测都是2012年以后开始，采用了原理更加先进的监测仪器，因此，当观测PM_2.5与PM₁₀浓度出现倒挂时，我们认为PM_2.5测量结果更加可信，而将PM₁₀观测数据剔除。由于各类PM_2.5浓度异常也可能导致“倒挂”问题的出现，所以本项检查放在其它检查之后，以减少错误剔除。

在本公开实施例中，步骤D所述的四项补充检查，还包括有效数据量检查，经上述系列方法进行质控后，能剔除大部分监测数据异常，但并不能保证将所有异常进行剔除，对一段时间内连续异常数据会由个别异常没有识别出来，对于时间上孤立的观测，尽管其可能与真值差异较小，但在周边没有数据进行佐证情况下，可信度较低。故而本公开中对每个数据，统计其前后12小时内的有效数据，若有效数据少于5个则对其进行剔除。

至此，已经结合附图对本公开实施例进行了详细描述。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：

(1)在观测数据严格遵守数据格式时，或直接调用数据库时可省略初级检查中的完整性检查；

(2)小变化异常检查和周期性异常检查的顺序可以互换；

(3)围绕时空一致性检查的不足设计了完整性检查、超量程检查、大观测误差检查、小变化检查、周期性异常检查、PM₁₀＜PM_2.5异常检查、有效数据量检查，若省略其中的一项或多项，异常识别效果将降低，但仍能识别大部分异常数据；

(4)在时间一致性回归中，使用了低通滤波，该滤波系数可以根据特定需要进行微调。

依据以上描述，本领域技术人员应当对本公开空气质量站点监测数据异常的自动化识别方法有了清楚的认识。

综上所述，本公开提供了一种空气质量站点监测数据异常的自动化识别方法，通过将大气污染监测数据异常分类后，利用计算机的快速运算性能可以实时或大规模的监测数据应用，实现监测数据异常的自动化识别，以缓解现有技术中即自动化异常识别方法中难以识别其特有的周期性异常、延滞异常、以及对于正定的，更接近于对数正态的空气质量监测数据(PM_2.5，PM₁₀，SO₂，NO₂，CO和O₃)的异常识别效果较差，难以识别数值较低或观测误差小于观测标准差的异常数据等技术问题。

还需要说明的是，实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本公开的保护范围。贯穿附图，相同的元素由相同或相近的附图标记来表示。在可能导致对本公开的理解造成混淆时，将省略常规结构或构造。

并且图中各部件的形状和尺寸不反映真实大小和比例，而仅示意本公开实施例的内容。另外，在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的限制。

除非有所知名为相反之意，本说明书及所附权利要求中的数值参数是近似值，能够根据通过本公开的内容所得的所需特性改变。具体而言，所有使用于说明书及权利要求中表示组成的含量、反应条件等等的数字，应理解为在所有情况中是受到「约」的用语所修饰。一般情况下，其表达的含义是指包含由特定数量在一些实施例中±10％的变化、在一些实施例中±5％的变化、在一些实施例中±1％的变化、在一些实施例中±0.5％的变化。

再者，单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。

说明书与权利要求中所使用的序数例如“第一”、“第二”、“第三”等的用词，以修饰相应的元件，其本身并不意味着该元件有任何的序数，也不代表某一元件与另一元件的顺序、或是制造方法上的顺序，该些序数的使用仅用来使具有某命名的一元件得以和另一具有相同命名的元件能做出清楚区分。

此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。

类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (9)

1.一种空气质量站点监测数据异常的自动化识别方法，包括：

步骤A：接收站点监测数据；

步骤B：对步骤A所接收的站点监测数据进行初级检查，识别显著异常的观测数据；

步骤C：对步骤B完成初级检查后的数据进行时空一致性检查，识别时空不一致异常；

步骤D：经步骤C后，使用四项补充检查识别未被初级检查和时空一致性检查识别的异常；以及

步骤E：输出经过质控后的观测数据、时空一致性估计值、及各项检查中的概率值；

步骤D所述的四项补充检查，包括：

小变化异常检查，观测值呈现出长时间常值或过于缓慢的异常时段，与实际大气污染变化特征不吻合，所述异常时段数据剔除；

周期性异常检查，识别周期出现的异常并进行剔除；

PM₁₀＜PM_2.5异常检查，当PM_2.5与PM₁₀浓度出现倒挂时，将PM₁₀观测数据剔除；以及

有效数据量检查，即对每个观测数据，统计其前后12小时内的有效数据，若有效数据少于5个则对其进行剔除。

2.根据权利要求1所述的自动化识别方法，其中，所述步骤B中的初级检查，包括：

步骤B₁：完整性检查；

步骤B₂：超量程检查，对监测数据进行上下限检查，将超出仪器量程的错误记录剔除；以及

步骤B₃：大观测误差检查，剔除超出合理值很多的观测，以减弱其对时空连续性检查性能的影响。

3.根据权利要求1所述的自动化识别方法，其中，步骤C中所述时空一致性检查，包括：

时间一致性回归；以及

空间一致性回归。

4.根据权利要求3所述的自动化识别方法，其中，所述时间一致性回归，利用检验点邻近时刻的观测数据，计算检验点的时间回归值，回归方法采用低通滤波，即：

其中F_t为滤波估计值，i为检验点的时次，k代表滤波时间窗口从检验点往前和往后的时间长度，f为原始观测，h为滤波系数。

5.根据权利要求3所述的自动化识别方法，其中，所述空间一致性回归，是结合邻近空间范围内的观测值计算得到检验点的估计值，具体计算公式如下：

其中F_s(i)为目标站点在检验点的时次i的空间一致性估计值，f_r为第r个参考站点的观测值，a_r为检验站点与参考站点间的一致性指标，采用以下方法进行计算：

其中f_r(i+k)为参考站点在i+k时刻的观测值，

为滑动窗口内的观测平均值。

6.根据权利要求5所述的自动化识别方法，其中，所述空间一致性回归，权重c_r采用Gaspari-Cohn(高斯-康恩)方案计算：

其中d为目标站点与参考站点之间的距离，d_c为截止距离。

7.根据权利要求4或5所述的自动化识别方法，其中，根据所述时间和空间的一致性估计值F_t和F_s，计算检验点的归一化估计残差Z_t和Z_s，再计算残差相关系数：

进而计算残差概率：

其中i为检验点的时次，ρ为时空残差的相关系数，Z_t，Z_s分别为归一后的时间和空间的回归残差，

分别为滑动窗口内时间和空间的归一残差平均，i-n和i+n分别为滑动窗口起始和结束时间。

8.根据权利要求1所述的自动化识别方法，其中，所述小变化异常检查中异常时段的残差概率P_a计算如下式：

其中R_a、S_a、Z_a为延滞时段的回归残差、回归残差标准差和归一化的回归残差，R_s、S_s为上面内容中计算的空间回归残差及其标准差，b和e分别为延滞时段的开始和结束时次，i为检验点的时次，计算得到小变化异常残差概率P_a，残差概率P_a小于设定阈值时的观测值识别为异常并剔除，所述阈值为10^-3～10^-9。

9.根据权利要求1所述的自动化识别方法，其中，所述周期性异常检查，首先，以24小时为间隔，对原始观测f进行滑动平均计算，如下：

其中i为待检查时次，然后，对f_c进行中值滤波，如下：

F_c(i)＝M(f_c(i+k)，k∈[-1，1])

其中M为集合的中位数，k代表滤波时间窗口从检验点往前和往后的时间长度，计算得到回归值F_c，滑动窗口长度取3个时次，接着，通过f_c和F_c计算得到回归残差R_c，并以94百分位的回归残差作为回归残差的标准差，使得得到的标准差S_c大于一天中第二大的回归残差，公式如下：

S_c(i)＝g(R_c(i+k)，k∈[-72，72])

其中g为集合的94百分位，最后，将R_c和S_c一起代入下列公式：

计算得到周期性异常残差概率P_c，残差概率P_c小于阈值的观测识别为异常并进行剔除，所述阈值为10^-2～10^-4。

Images