工業園區污染廢氣的動態關聯分析方法

朱 濤， 王曉凱， 衛曉旭， 凌德森

(山西大學 物理電子工程學院， 山西 太原 030006)

0 引 言

隨著改革開放推動國民經濟的發展， 工業園區逐漸建立起來， 它推動著國家工業和制造業的發展， 國家GDP的60%都是由各級各類園區占據[1]. 盡管工業園區給國家的經濟帶來大幅度的推動， 但是在環境污染方面卻出現了很多問題. 2019 年， 全國337個地級及以上城市中， 仍有180個城市環境空氣質量超標[2]， 全國環境污染問題仍然嚴重. 如何快速地監測和追蹤廢氣污染源頭是目前環境治理的重中之重.

工業園區內企業眾多， 地形復雜， 再加上天氣、 風向、 季節、 壓強等眾多影響廢氣濃度分布的因素， 尋找污染因素的工作不好開展. 目前， 國內外對于環境污染治理采用的代表性方法是氣體擴散模型和關聯分析方法. 氣體擴散模型是一種以概率模型為基礎， 模擬氣體擴散的模型， 以高斯煙羽模型[3]和高斯煙團模型[4]為主要代表. 張成才等[5]利用高斯擴散模型建立廢氣污染擴散預測系統， 模擬工業廢氣污染對周邊空氣質量的影響. 李萬莉[6]通過建立改進的高斯煙羽模型， 在天然氣泄漏時預測天然氣的濃度分布有很好的效果. 擴散模型主要通過污染源信息得到污染物的擴散和分布， 然而在前向分析情況下， 由于假設的參數過量， 高斯模型很難反演計算， 當設置的參數不準確時， 反向模型會輸出不同的推理結果. 所以， 在變量影響的推斷和源頭追蹤方面， 擴散模型存在很大的不足.

現代工業園區廢氣擴散過程中會涉及到氣體之間的相互作用以及氣體與外界之間的相互作用， 這導致工業園區周邊環境的每一種物質會因為一定的相互作用而聯系在一起[7 ]， 所以廢氣污染物關聯分析方法是確定污染因素、 推理污染源頭的有效方法. 目前， 關聯分析屢見不鮮， 主成分分析法、 灰度關聯分析法、 斯皮爾曼相關系數(Spearman)、 最大信息系數(Maximal Information Coefficient, MIC)等在各行各業都有廣泛應用. Qiao Z等[8]利用偏相關和層次聚類分析方法， 首次在全國范圍內探討了空氣污染指數(API)與多個氣象參數之間的關系, 提出了空氣質量對氣象條件敏感性的時空變化見解. Spearman相關系數和最大信息系數[9]都是度量兩個變量之間的統計學關聯系數， 均可用于線性和非線性數據； 范秋香[10]采用Spearman相關系數， 分析了臨沂市經濟技術開發區4年內污染排放量的濃度分布變化， 得出O3濃度呈現顯著上升趨勢， 污染顆粒物PM10， PM2. 5和大氣主要污染物SO2的年均值呈顯著的下降趨勢； 朱青等[11]通過最大信息系數等分析2014年鄱陽湖生態環境質量空間分布特征， 得出影響原始遙感生態指數(RSEI0)的關鍵影響因素， 為湖區環境污染治理提供了關鍵信息.

由于影響空氣污染的不確定因素很多， 需要進行全面的關聯分析， 相關分析方法層出不窮， 選用合適的相關分析方法作為本文的關聯方法非常重要. 傳統的關聯分析只能看出廢氣污染物與影響因素之間的長期關聯度， 無法從時間維度上分析廢氣污染物與其影響因素的時序特征， 并且效率差， 分析效果不佳. 關聯分析不僅要做到分析準確， 還要做到高效率. 為了解決此問題， 本文重點對動態分析進行研究， 為了做到動態分析， 本文提出含自適應滑窗的動態分析方法， 配合Spearman相關系數和最大信息系數加權來分析廢氣污染物與各個影響因素之間的動態關聯度.

1 動態關聯分析總體框架

目前， 在國家環境治理的大力監督管理下， 各城市已經形成了網格化的監測系統， 像工業園區這種工業廢氣主要排放區內通常安裝有不同位置的工業廢氣濃度傳感器， 組成了工業污染監測站點， 環保部門通過各個站點采集數據， 監測著工業園區的環境質量. 本文通過分析某站點的數據， 采用動態關聯分析方法來分析與異常污染物關聯度高的其他污染物. 圖 1 為動態關聯分析方法總框架， 主要分為3個步驟： 數據獲取、 動態關聯分析和篩選異常指標和時刻.

圖 1 動態關聯分析方法總體框架Fig.1 The general framework of dynamic correlationanalysis method

當某監測站某污染物濃度異常(工業主要廢氣污染物濃度數據超標或者即將超標)， 獲取該污染物歷史濃度數據， 一般獲取一個月內的歷史數據， 并獲取該站點其他污染物歷史濃度數據； 利用含滑窗的動態關聯分析方法計算目標污染物與該站點其他污染物之間的關聯度， 該關聯度為動態指標， 具有時序性； 通過對關聯度的分析， 獲取與該站點異常污染物具有高關聯性的一種或幾種廢氣污染物及其時刻； 根據分析結果對工業園區進行針對性的監督和管理.

2 動態關聯分析方法

2.1 含變異系數的動態滑窗設計

變異系數是衡量數據離散大小的統計量， 它的計算方法是標準差和平均值的比值， 它可以消除量綱帶來的影響[12]. 以變異系數為基礎， 可以計算一段時序數據的變異系數變化率， 從而反應該段數據的波動情況， 通過設置合理的閾值來判斷該段數據的波動范圍， 并以反饋的形式來調整需要關聯分析的數據長度.

在分析每段時間序列的關聯性時， 為了加快算法效率但不丟失關鍵特征， 需要考慮數據的波動情況， 采用自適應滑窗的方法來處理數據， 當數據波動過于平穩時， 自適應滑窗將加長截取窗口的長度， 擴大關聯分析的范圍； 當數據波動過于劇烈時， 自適應滑窗將縮短截取窗口的長度， 減小關聯分析的范圍， 提高瞬時特性辨識度. 在此思路下， 提出基于變異系數的自適應滑窗確定方法， 具體包括以下步驟：

Step 1： 根據設置的初始窗長L0截取數據段， 并求取截取的數據段的變異系數C

(1)

式中：σ0是該段數據的標準差；m0是該段數據的平均值.

Step 2： 將該數據段再分成s段， 每段長度為L0/s；

Step 3： 求取每小段的變異系數Ci

(2)

式中：σi是第i小段數據的標準差；mi是第i小段數據的平均值.

Step 4： 求取該數據段的變異系數變化率

(3)

Step5： 滑窗長度調整比例

(4)

且

(5)

式中：kmax，kmin為穩定性判別閾值；L為更新后的窗口長度.

2.2 Spearman相關系數和最大信息系數

Spearman相關系數和最大信息系數都是度量兩個變量之間的統計學關聯系數， 均可用于線性和非線性數據. 而且兩者計算復雜度低， 魯棒性高， 對于樣本的數量沒有過多的要求， 均可以處理小樣本數據， 更適用于污染廢氣動態關聯分析方法.

2.2.1 Spearman相關系數

Spearman相關系數可以用來量度兩個波形的相關程度， 取值在-1和1之間. 若兩個波形在一定范圍內呈正相關， 存在高度的相似性， 則取值為正值， 絕對值越高， 關聯度越高； 若兩個波形在一定范圍內呈負相關， 波形存在相反的趨勢， 則取值為負值， 絕對值越高， 負相反的程度越高. Spearman相關系數具體思想為：

將時序數據x和y分別按照升序或者降序來排列數據， 將數據x和y內的每個元素在排列中的位置記作該元素的秩次， 從而得到數據x和y的秩次數列r和e， 將數列r和數列e內每個元素對應相減得到秩次差數列d={d1，d2，…，dn}， 再將其帶入斯皮爾曼相關系數公式[13]

(6)

式中：n為數據樣本量；ρ為斯皮爾曼相關系數；i為第i個樣本.

2.2.2 最大信息系數

最大信息系數也可以用來反應兩個波形之間的關聯程度， 以互信息為基礎， 采用網格劃分的方法， 相較于互信息而言有更高的準確度， 具有普適性、 公平性和對稱性. 其計算過程如下：

對于給定的數據集D={(Xi,Yi),i= 1,2,…,n}， 在直角坐標系中， 如果將X軸劃分為f個格子，Y軸劃分為g個格子， 可以得到一個f×g的網絡劃分G， 將落入G的點的數量占數據集D數量的比值看作是其概率密度D|G， 而根據不同的網格劃分情況得到的概率分布D|G.變量X和Y的最大互信息為[14]

(7)

式中：D是給定的數據集；f,g是對這個數據集的某種劃分數量；p(X,Y)是聯合概率密度；p(X) 和p(Y)是邊緣概率密度.

相同f×g的網絡劃分有很多種劃分方式， 會得到不同大小的互信息值， 記錄不同劃分情況下最大的互信息值為MI(D,X,Y).之后再進行歸一化， 使其取值在[0,1]之間

(8)

取不同網格劃分下的最大的歸一化后的互信息值作為最大信息系數的值. 假設樣本數量為n， 最大信息系數的定義為

(9)

式中：fg

2.2.3 Spearman相關系數和MIC加權

Spearman相關系數和最大信息系數可以衡量數據變量之間的線性關系和非線性關系， 但是最大信息系數只能表現出兩變量的某種函數關聯度， 無法表征出負相關的程度， 而Spearman相關系數具備這一點， 但是Spearman相關系數的魯棒性沒有最大信息系數的高. 所以鑒于兩種系數的優點， 采用Spearman相關系數和最大信息系數加權的方式作為相關方法： 當Spearman相關系數和最大信息系數都為正值時， 取兩者的平均值； 當Spearman相關系數為負值時， 取兩者的絕對值的平均值， 然后均值再加負號.

2.3 動態關聯分析算法

動態關聯分析方法是一種能對時序數據在時間維度上進行動態關聯分析的算法， 依靠含變異系數的動態滑窗設計， 根據數據的波動情況， 實時調整窗口的長度， 配合Spearman相關系數和最大信息系數的加權相關法， 將研究對象的長期關聯度切割成若干個瞬時關聯度， 反應出數據的動態特性， 工業園區污染廢氣動態關聯分析算法流程如圖 2 所示. 大致流程如下：

Step 1： 獲取異常污染物濃度數據和其他污染物濃度數據.

圖 2 工業園區污染廢氣動態關聯分析算法流程Fig.2 Dynamic correlation analysis algorithm flow ofindustrial park pollution exhaust gas

Step 2： 截取T時刻及以前一定長度L0異常污染物歷史數據，L0為初期調研實驗后得出的最佳值.

Step3： 通過變異系數及其變化率判斷數據是否穩定， 具體判斷規則如2.1節所示，kmax，kmin為初期調研實驗后得出的最佳值.

Step 4： 如果數據穩定， 使用Spearman相關系數和MIC加權相關方法進行關聯分析； 若數據不穩定， 則通過調整規則調整滑窗長度， 然后再進行關聯分析.

Step 5： 判斷所有數據是否全部完成動態關聯分析， 沒有的話， 截取T+1時刻重復step 2到step 4的步驟.

Step 6： 如果數據已全部完成動態關聯分析， 輸出動態關聯結果.

3 實驗及結果

3.1 實驗數據及實驗設置

以某工業園區CO濃度數據作為本文的主要研究對象， 以2017年4月到6月的數據(時均值)作為實驗數據， 其中包括SO2, NO2, CO, O3、 溫度、 濕度等指標(因為示例數據是非化工園區采集的數據， 所以污染物主要以碳硫化物、 氮氧化物為主). 本文共設置3組實驗， 實驗1： 取2017年5月和6月的CO數據進行含變異系數的自適應滑窗實驗， 數據長度為660 h， 觀察數據波動情況與窗長的變化. 實驗2： 取2017年6月的CO數據分別與NO2和SO2數據進行含滑窗的動態關聯實驗， 數據長度為 660 h， 用不同的相關方法(皮爾遜相關系數(Pearson)、 Spearman相關系數、 最大信息系數、 肯德爾系數(Kendall)、 本文中的加權相關方法)進行比較. 實驗3： 以2017年4月的CO數據為研究對象， NO2, SO2為關聯對象進行本文中的動態關聯分析實驗， 數據長度為660 h(時間尺度約為1個月)， 初始窗長L0均設置為48 h，kmax設置為0.45，kmin設置為0.2，s設置為 8 h.

3.2 實驗1

取2017年5月和6月的CO數據進行含變異系數的動態滑窗實驗， 觀察數據波動情況與窗長的變化. 2017年5月的CO數據定義為數據1， 2017年6月的CO數據定義為數據2.

圖 3 為數據1數據波動情況和窗口長度變化情況， 圖 4 為數據2數據波動情況和窗口長度變化情況.

圖 3 窗口變化情況(數據1)

圖 4 窗口變化情況(數據2)Fig.4 Window changes (data 2)

圖 3(a)和圖4(a)為CO數據原圖， 圖3(b)和圖4(b)為2.1節介紹的滑動窗口調整比例k值的變化， 圖3(c)和圖4(c)為更新后窗口長度的變化. 由圖可知在數據波動較大的區域， 窗口可以及時地變小， 在數據相對平穩的區域， 窗口可以適當的變大. 數據1相較于數據2比較平穩， 數據2中CO濃度最高峰值達到3 mg/m3， 相應的數據1的窗口長度加長的時刻比較多， 數據2的窗口長度縮短的時刻比較多.

3.3 實驗2

取2017年6月的CO濃度數據分別與NO2和SO2濃度數據進行含滑窗的動態關聯實驗， 為了方便觀察只截取前100 h， 用不同的相關方法(Pearson系數、 Spearman相關系數、 MIC、 Kendall系數、 本文中的加權相關方法(Spearman-MIC))進行比較. 圖 5 為CO濃度與NO2濃度的動態關聯分析， 圖 6 為CO濃度與SO2濃度的動態關聯分析， 從圖中可以看出Pearson系數忽高忽低， Kendall系數普遍偏低， 本文中的Spearman相關系數和MIC加權的相關方法始終處于中間的位置， 更能準確地表征兩變量的關聯特征.

圖 5 CO-NO2關聯度Fig.5 CO-NO2 correlation degree

圖 6 CO-SO2關聯度Fig.6 CO-SO2 correlation degree

3.4 實驗3

以2017年4月的CO數據為研究對象， NO2, SO2為關聯對象進行本文中的動態關聯分析實驗. 圖 7、 圖 9 分別展示了CO與NO2, SO2的動態關聯分析結果. 并且展示了CO與這兩種不同污染物的長期靜態關聯度. 兩個不同的指標都有較為突出的高關聯時間段， 說明此時此刻該指標對于CO濃度有著較大的影響. 圖7中大約有7個以上峰值超過0.8， 觀測點550往后有著較高的關聯度， 說明這些時刻CO和NO2關系密切， 在圖9中大部分的區域都在長期關聯度以上， 說明CO濃度和SO2有著很大的關系. 圖 8、 圖10分別展示了高關聯度時間段內， CO與NO2, SO2的波形變化圖(因不同污染物數據量綱不同， 為了方便觀察故作歸一化處理). 通過波形圖可以看出， 當NO2濃度和SO2濃度升高后， CO濃度也隨之升高. 根據這些信息進行污染的專項防治， 為CO濃度異常提供了關鍵信息.

圖 7 CO與NO2動態關聯度Fig.7 Dynamic correlation degree between CO and NO2

圖 8 CO與NO2高關聯度區間波形變化Fig.8 Variation of high correlation interval waveformbetween CO and NO2

圖 9 CO與SO2動態關聯度Fig.9 Dynamic correlation degree between CO and SO2

圖 10 CO與SO2高關聯度區間波形變化Fig.10 Variation of high correlation interval waveformbetween CO and SO

4 結束語

針對工業園區廢氣污染治理問題， 本文基于變異系數的動態滑窗提出了一種動態關聯分析方法. 該方法使用變異系數及其變化率判斷數據段波動性自適應地調整窗口大小， 借助Spearman相關系數和最大信息系數加權來反映關聯特征， 最后使用某工業園區的實際監測數據中的CO濃度與其他污染氣體濃度進行動態關聯分析實驗. 結果表明， 該方法可以有效提取數據的動態關聯特性， 能夠很好地反應出CO濃度異常的其他影響因子和時間信息， 為工業園區廢氣污染溯源提供了重要的參考信息. 文中選用的CO氣體類型是示例， 也可以用于其他氣體以及其他影響因素的計算. 未來可以分析多站點之間動態關聯度， 在空間維度上將多個站點聯系起來， 提取數據的時空特征, 推理污染路線.