CN112634113B - 一种基于动态滑窗的污染废气关联分析方法 - Google Patents

Abstract

本发明公开了一种基于动态滑窗的污染废气关联分析方法，属于数值分析与数据挖掘技术领域。本发明首先使用变异系数及其变化率来定义滑窗窗口长度的变化，其次将斯皮尔曼相关系数Spearman和最大信息系数MIC加权来反映关联特征，最后将两种方法结合起来通过滑动窗口自适应的提取数据，辅助斯皮尔曼相关系数和最大信息系数加权相关法来分析数据间的动态关联度。本发明具有动态特性，可将两个研究对象在时间维度上，将恒定关联度切割成若干瞬时关联性。本发明中的Spearman相关系数和MIC加权相关方法较其他相关方法准确度高，更能准确的表征两变量的关联特征，能较好的分析出与异常污染物有关的污染因素和时间信息，为大气污染溯源提供了重要的参考信息。

Description

一种基于动态滑窗的污染废气关联分析方法

技术领域

本发明属于数值分析与数据挖掘技术领域，具体涉及一种基于动态滑窗的污染废气关联分析方法。

背景技术

随着改革开放推动国民经济的发展，工业园区逐渐建立起来，它是众多工业企业的聚集地，维持和推进着各个城市的工业发展，推动着国家工业和制造业的发展，是国家产能的发动机。截至2018年，我国共有国家级经济技术开发区219家，国家级高新技术开发区168家，省级及其他开发区2000多家，各级各类园区的经济总量占全国GDP的60％以上。尽管工业园区给国家的经济带来大幅度的推动，但是在环境污染方面却出现了很多问题。由于工业园区里百分之八十都是工业工厂，工业废气的排放量也随着工业园区的兴盛逐渐增加，大部分工业园区虽然建立在离市区很远的地方，但是工业废气在季节和风向的影响下依然影响着城市上空的空气质量。2019年全国337个地级及以上城市中，180个城市环境空气质量超标，占53.4％，106个城市优良天数比例在50％～80％之间，337个城市累计发生严重污染452天，全国环境污染问题仍然严重。而环境空气治理的源头除了汽车废气的排放外，其次就是工业园区的废气排放，据统计，工业园区二氧化硫排放占全国二氧化硫排放总量的90％，氮氧化物占70％，烟尘占85％。如何快速的监测和追踪废气污染源头是目前环境治理的重中之重，废气污染物关联分析方法是确定污染因素，推理污染源头的有效方法。

发明内容

针对上述问题本发明提供了一种基于动态滑窗的污染废气关联分析方法。

为了达到上述目的，本发明采用了下列技术方案：

一种基于动态滑窗的污染废气关联分析方法，包括以下步骤：

步骤1，含变异系数的动态滑窗；

步骤2，Spearman相关系数和MIC的加权；

步骤3，动态关联分析算法。

进一步，所述步骤1含变异系数的动态滑窗还包括以下步骤：

步骤S1.1：根据设置的初始窗长L₀截取数据段，并求取截取的数据段的变异系数C：

其中σ₀是该段数据的标准差，m₀是该段数据的平均值；

步骤S1.2：将该数据段再分成s段，每段长度为L₀/s；

步骤S1.3：求取每小段的变异系数C_i：

其中σ_i是第i小段数据的标准差，m_i是第i小段数据的平均值。

步骤S1.4：求取该数据段的变异系数变化率R：

步骤S1.5：滑窗长度调整比例k为：

且：

其中k_max，k_min为稳定性判别阈值，L为更新后的窗口长度。

在分析每段时间序列的关联性时，为了加快算法效率但不丢失关键特征，需要考虑数据的波动情况，采用滑窗的方法来处理数据，当数据波动过于平稳时，动态滑窗将加长截取窗口的长度，扩大关联分析的范围；当数据波动过于剧烈时，动态滑窗将缩短截取窗口的长度，减小关联分析的范围，提高瞬时特性辨识度。

进一步，所述步骤2Spearman相关系数和MIC的加权还包括以下步骤：

步骤S2.1，首先把时序数据x按升序或者降序排列得到排序数列a＝{a₁，a₂，…，a_n}，将数据x内每个元素在数列a中的位置记为r_i，称其为该元素的秩次，从而得到数据x对应的秩次数列r；

将另一时序数据y按升序或者降序排列得到排序数列b＝{b₁，b₂，…，b_n}，继而可以得到时序数据y对应的秩次数列e；

将数列r和数列e内每个元素对应相减得到秩次差数列d＝{d₁，d₂，…，d_n}，再将其带入斯皮尔曼相关系数公式：

式中：n为数据样本量，ρ为斯皮尔曼相关系数，i为第i个样本；

Spearman相关系数用来量度两个波形相关程度的一种方法。若两个波形在一定范围内呈正相关，存在高度的相似性，则取值为正值，绝对值越高，关联度越高；若两个波形在一定范围内呈负相关，波形存在相反的趋势，则取值为负值，绝对值越高，相反的程度越高。Spearman相关系数对波形的分布不作要求，适用范围广，不限样本大小。

步骤S2.2，对于给定的数据集D＝{(X_i,Y_i),i＝1,2,…,n}，如果将X轴划分为f个格子，Y轴划分为g个格子，那么我们可以得到一个f×g的网络划分G，将落入G的点的数量占数据集D数量的比例看作是其概率密度D|_G，而根据不同的网格划分情况得到的概率分布D|_G，变量X和Y的最大互信息I为：

其中，D是给定的数据集；f,g是对这个数据集的某种划分数量；p(X,Y)是联合概率密度，p(X)和p(Y)是边缘概率密度，

相同f×g的网络划分有很多种情况，会得到不同的互信息值，记录其中最大的互信息值为MI(D,X,Y)，之后在进行归一化，使其取值在[0,1]之间：

取不同网格划分下的最大的归一化后的互信息值作为最大信息系数的值；假设样本数量为n，最大信息系数MIC(D)的定义为：

其中fg<B(n)为网格划分的上界，一般情况下B(n)＝n^0.6效果最佳。

最大信息系数经常用来衡量两个基因之间的关联程度，以互信息为基础，采用网格划分的方法，相较于互信息而言有更高的准确度，具有普适性、公平性和对称性。

斯皮尔曼相关系数和最大信息系数可以衡量数据变量之间的线性关系和非线性关系，但是最大信息系数只能表现出两变量的某种函数关联度，无法表征出负相关的程度，而斯皮尔曼具备这一点，但是斯皮尔曼相关系数的鲁棒性没有最大信息系数的高。所以鉴于两种系数的优点，采用斯皮尔曼相关系数和最大信息系数加权的方式作为相关方法：当斯皮尔曼系数和最大信息系数都为正值时，取两者的平均值；当斯皮尔曼系数为负值时，取两者的绝对值的平均值，然后均值再加负号。

进一步，所述步骤3动态关联分析算法还包括以下步骤：

步骤S3.1，获取异常污染物浓度数据和其他污染物浓度数据；

步骤S3.2，截取T时刻及T时刻前一定长度L₀异常污染物历史数据，L₀为初期调研实验后得出的最佳值。

步骤S3.3，通过步骤S1的变异系数和变异系数变化率判断数据是否稳定，k_max，k_min为初期调研实验后得出的最佳值。

步骤S3.4，当数据稳定时，再进行Spearman相关系数和最大信息系数加权相关方法进行关联分析；若数据不稳定，则通过调整规则调整滑窗长度，然后再进行关联分析。

步骤S3.5，判断所有数据是否全部完成动态关联分析，没有全部完成，则截取T+1时刻重复步骤S3.2到步骤S3.4。

步骤S3.6，当数据已全部完成动态关联分析，输出动态关联结果。

所述步骤3动态关联分析算法，将上述两个步骤结合，是一种能对时序数据在时间维度上进行动态关联分析的算法，依靠含变异系数的动态滑窗设计，根据数据的波动情况，实时调整窗口的长度，配合Spearman相关系数和最大信息系数的加权相关法，将研究对象的长期关联度切割成若干个瞬时关联度，反应出数据的动态特性。

与现有技术相比本发明具有以下优点：较于长期关联方法，本发明具有动态特性，可以将两个研究对象在时间维度上，将恒定关联度切割成若干瞬时关联性。本发明中的Spearman相关系数和MIC加权相关方法较其他相关方法准确度高，更能准确的表征两变量的关联特征，能较好的分析出与异常污染物有关的污染因素和时间信息，为大气污染溯源提供了重要的参考信息。

附图说明

图1大气污染物动态关联分析算法流程；

图2数据波动情况和窗口变化情况；

图3为CO与NO₂关联度不同相关方法比较；

图4为CO与NO₂动态关联度；

图5为CO与SO₂动态关联度。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本发明实施例给出一种基于动态滑窗的污染废气关联分析方法，如附图1所示，其为本发明实施例中提出的大气污染物动态关联分析算法流程，该方法具体包括：

以某工业园区CO浓度数据作为本文的主要研究对象，以2017年4月到6月的数据(时均值)作为实验数据，其中包括SO₂、NO₂、CO、O₃、温度、湿度等指标。

取2017年6月的CO浓度数据进行含变异系数的动态滑窗实验。分析该段时间序列的关联性时，为了加快算法效率但不丢失关键特征，需要考虑数据的波动情况，采用滑窗的方法来处理数据，当数据波动过于平稳时，动态滑窗将加长截取窗口的长度，扩大关联分析的范围；当数据波动过于剧烈时，动态滑窗将缩短截取窗口的长度，减小关联分析的范围，提高瞬时特性辨识度。在此思路下，用该段数据进行动态滑窗确定方法，具体包括以下步骤：

步骤S1.1：根据设置的初始窗长48截取数据段，并求取截取的数据段的变异系数C：

其中σ₀是该段数据的标准差，m₀是该段数据的平均值。

步骤S1.2：将该数据段再分成8段，每段长度为6；

步骤S1.3：求取每小段的变异系数C_i：

其中σ_i是第i小段数据的标准差，m_i是第i小段数据的平均值。

步骤S1.4：求取该数据段的变异系数变化率R：

步骤S1.5：滑窗长度调整比例k为：

且：

其中k_max＝0.45，k_min＝0.2为稳定性判别阈值。

附图2为该数据波动情况和窗口长度变化情况。在附图2中，图a为CO数据原图，图b为滑动窗口调整比例k值的变化，图c为窗口长度的变化。由图可知在数据波动较大的区域，窗口可以及时的变小，在数据相对平稳的区域，窗口可以适当的变大。

用该CO浓度数据与NO₂浓度进行含滑窗的动态关联分析实验，步骤S2.1首先把该CO浓度数据x按升序或者降序得到排序数列a＝{a₁，a₂，…，a_n}，将数据x内每个元素在数列a中的位置记为r_i，称其为该元素的秩次，从而得到数据x对应的秩次数列r。将NO₂浓度数据y按同样方式排列得到排序数列b＝{b₁，b₂，…，b_n}，继而可以得到时序数据y对应的秩次数列e。将数列r和数列e内每个元素对应相减得到秩次差数列d＝{d₁，d₂，…，d_n}，再将其带入斯皮尔曼相关系数公式：

式中：n为数据样本量，ρ为斯皮尔曼相关系数，i为第i个样本。

步骤S2.2，对于上述给定的数据集D＝{(X_i,Y_i),i＝1,2,…,n}，X_i为上述CO浓度，Y_i为上述NO₂浓度。如果将X轴划分为f个格子，Y轴划分为g个格子，那么我们可以得到一个f×g的网络划分G，将落入G的点的数量占数据集D数量的比例看作是其概率密度D|_G，而根据不同的网格划分情况得到的概率分布D|_G。变量X和Y的最大互信息I为：

其中，D是给定的数据集；f,g是对这个数据集的某种划分数量；p(X,Y)是联合概率密度，p(X)和p(Y)是边缘概率密度。

相同f×g的网络划分有很多种情况，会得到不同的互信息值，记录其中最大的互信息值为MI(D,X,Y)。之后在进行归一化，使其取值在[0,1]之间：

取不同网格划分下的最大的归一化后的互信息值作为最大信息系数的值。假设样本数量为n，最大信息系数MIC(D)的定义为：

其中fg<B(n)为网格划分的上界，一般情况下B(n)＝n^0.6效果最佳。

将上述取得的Spearman相关系数ρ和最大信息系数MIC(D)进行加权，加权规则为：当斯皮尔曼系数和最大信息系数都为正值时，取两者的平均值；当斯皮尔曼系数为负值时，取两者的绝对值的平均值，然后平均值再加负号。用不同的相关方法(Pearson系数、斯皮尔曼相关系数、最大信息系数、Kendall系数、本文中的加权相关方法)进行比较。附图3为CO浓度与SO₂浓度的不同相关方法动态关联分析比较，从图中可以看出Pearson系数普遍偏高，Kendall系数普遍偏低，本文的Spearman相关系数和MIC加权的相关方法始终处于中间的位置，能更准确的表征两变量的关联特征。

以CO为研究对象，NO₂、SO₂浓度数据为关联对象完整的进行本文中的动态关联分析实验。附图4、5为CO浓度和NO₂、SO₂浓度的动态关联度，并且展示了CO与这两种不同指标的长期静态关联度。图4中大约有7个以上峰值超过0.8，观测点550往后有着较高的关联度，说明这些时刻CO和NO₂关系密切，在图5中大部分的区域都在长期关联度以上，说明CO浓度和SO₂有着很大的关系。根据这些信息进行污染的专项防治，为CO浓度异常提供了关键信息。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (1)

1.一种基于动态滑窗的污染废气关联分析方法，其特征在于：包括以下步骤：

步骤1，含变异系数的动态滑窗；

步骤2，Spearman相关系数和MIC的加权；

步骤3，动态关联分析算法；所述步骤1含变异系数的动态滑窗还包括以下步骤：

步骤S1.1：根据设置的初始窗长L₀截取数据段，并求取截取的数据段的变异系数C：

其中σ₀是该段数据的标准差，m₀是该段数据的平均值；

步骤S1.2：将该数据段再分成s段，每段长度为L₀/s；

步骤S1.3：求取每小段的变异系数C_i：

其中σ_i是第i小段数据的标准差，m_i是第i小段数据的平均值；

步骤S1.4：求取该数据段的变异系数变化率R：

步骤S1.5：滑窗长度调整比例k为：

且：

其中k_max，k_min为稳定性判别阈值，L为更新后的窗口长度；所述步骤2Spearman相关系数和MIC的加权还包括以下步骤：

步骤S2.1，首先把时序数据x按升序或者降序排列得到排序数列a＝{a₁，a₂，…，a_n}，将数据x内每个元素在数列a中的位置记为r_i，称其为该元素的秩次，从而得到数据x对应的秩次数列r；

将另一时序数据y按升序或者降序排列得到排序数列b＝{b₁，b₂，…，b_n}，继而可以得到时序数据y对应的秩次数列e；

将数列r和数列e内每个元素对应相减得到秩次差数列d＝{d₁，d₂，…，d_n}，再将其带入斯皮尔曼相关系数公式：

式中：n为数据样本量，ρ为斯皮尔曼相关系数，i为第i个样本；

步骤S2.2，对于给定的数据集D＝{(X_i,Y_i),i＝1,2,…,n}，如果将X轴划分为f个格子，Y轴划分为g个格子，那么我们可以得到一个f×g的网络划分G，将落入G的点的数量占数据集D数量的比例看作是其概率密度D|_G，而根据不同的网格划分情况得到的概率分布D|_G，变量X和Y的最大互信息I为：

其中，D是给定的数据集；f,g是对这个数据集的某种划分数量；p(X,Y)是联合概率密度，p(X)和p(Y)是边缘概率密度，

相同f×g的网络划分有很多种情况，会得到不同的互信息值，记录其中最大的互信息值为MI(D,X,Y)，之后在进行归一化，使其取值在[0,1]之间：

取不同网格划分下的最大的归一化后的互信息值作为最大信息系数的值；假设样本数量为n，最大信息系数MIC(D)的定义为：

其中fg<B(n)为网格划分的上界，B(n)＝n^0.6效果最佳；所述步骤3动态关联分析算法还包括以下步骤：

步骤S3.1，获取异常污染物浓度数据和其他污染物浓度数据；

步骤S3.2，截取T时刻及T时刻前一定长度L₀异常污染物历史数据，L₀为初期调研实验后得出的最佳值；

步骤S3.3，通过步骤S1的变异系数和变异系数变化率判断数据是否稳定，k_max，k_min为初期调研实验后得出的最佳值；

步骤S3.4，当数据稳定时，再进行Spearman相关系数和最大信息系数加权相关方法进行关联分析；若数据不稳定，则通过调整规则调整滑窗长度，然后再进行关联分析；

步骤S3.5，判断所有数据是否全部完成动态关联分析，没有全部完成，则截取T+1时刻重复步骤S3.2到步骤S3.4；

步骤S3.6，当数据已全部完成动态关联分析，输出动态关联结果。