朱 濤， 王曉凱， 衛(wèi)曉旭， 凌德森

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

0 引 言

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

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

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

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

1 動態(tài)關(guān)聯(lián)分析總體框架

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

圖 1 動態(tài)關(guān)聯(lián)分析方法總體框架Fig.1 The general framework of dynamic correlationanalysis method

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

2 動態(tài)關(guān)聯(lián)分析方法

2.1 含變異系數(shù)的動態(tài)滑窗設計

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

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

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

(1)

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

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

Step 3： 求取每小段的變異系數(shù)Ci

(2)

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

Step 4： 求取該數(shù)據(jù)段的變異系數(shù)變化率

(3)

Step5： 滑窗長度調(diào)整比例

(4)

且

(5)

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

2.2 Spearman相關(guān)系數(shù)和最大信息系數(shù)

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

2.2.1 Spearman相關(guān)系數(shù)

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

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

(6)

式中：n為數(shù)據(jù)樣本量；ρ為斯皮爾曼相關(guān)系數(shù)；i為第i個樣本.

2.2.2 最大信息系數(shù)

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

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

(7)

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

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

(8)

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

(9)

式中：fg

2.2.3 Spearman相關(guān)系數(shù)和MIC加權(quán)

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

2.3 動態(tài)關(guān)聯(lián)分析算法

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

Step 1： 獲取異常污染物濃度數(shù)據(jù)和其他污染物濃度數(shù)據(jù).

圖 2 工業(yè)園區(qū)污染廢氣動態(tài)關(guān)聯(lián)分析算法流程Fig.2 Dynamic correlation analysis algorithm flow ofindustrial park pollution exhaust gas

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

Step3： 通過變異系數(shù)及其變化率判斷數(shù)據(jù)是否穩(wěn)定， 具體判斷規(guī)則如2.1節(jié)所示，kmax，kmin為初期調(diào)研實驗后得出的最佳值.

Step 4： 如果數(shù)據(jù)穩(wěn)定， 使用Spearman相關(guān)系數(shù)和MIC加權(quán)相關(guān)方法進行關(guān)聯(lián)分析； 若數(shù)據(jù)不穩(wěn)定， 則通過調(diào)整規(guī)則調(diào)整滑窗長度， 然后再進行關(guān)聯(lián)分析.

Step 5： 判斷所有數(shù)據(jù)是否全部完成動態(tài)關(guān)聯(lián)分析， 沒有的話， 截取T+1時刻重復step 2到step 4的步驟.

Step 6： 如果數(shù)據(jù)已全部完成動態(tài)關(guān)聯(lián)分析， 輸出動態(tài)關(guān)聯(lián)結(jié)果.

3 實驗及結(jié)果

3.1 實驗數(shù)據(jù)及實驗設置

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

3.2 實驗1

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

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

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

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

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

3.3 實驗2

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

圖 5 CO-NO2關(guān)聯(lián)度Fig.5 CO-NO2 correlation degree

圖 6 CO-SO2關(guān)聯(lián)度Fig.6 CO-SO2 correlation degree

3.4 實驗3

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

圖 7 CO與NO2動態(tài)關(guān)聯(lián)度Fig.7 Dynamic correlation degree between CO and NO2

圖 8 CO與NO2高關(guān)聯(lián)度區(qū)間波形變化Fig.8 Variation of high correlation interval waveformbetween CO and NO2

圖 9 CO與SO2動態(tài)關(guān)聯(lián)度Fig.9 Dynamic correlation degree between CO and SO2

圖 10 CO與SO2高關(guān)聯(lián)度區(qū)間波形變化Fig.10 Variation of high correlation interval waveformbetween CO and SO

4 結(jié)束語

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