基于社區網絡的大氣污染源定位算法

宿兵暢 張亞娟 張素琪

1(河北工業大學人工智能與數據科學學院 天津 300401) 2(河北省大數據計算重點實驗室 天津 300401) 3(天津商業大學信息工程學院 天津 300134)

0 引 言

隨著工業的發展和城鎮化速度的不斷加快，我國大氣污染問題由間歇性突發問題逐漸演變成常態性多發問題[1]，只有從復雜的污染數據中分析出污染源，才能從根本上解決大氣污染問題，因此污染源定位對于治理大氣污染具有重要意義。

傳統的污染源定位算法是基于氣體濃度衰減模型的定位方法。這類方法以氣體擴散理論為基礎，通過氣體擴散過程中污染物濃度和擴散距離的關系，計算出污染源的位置[2]。2006年匡興紅等[3]基于氣體污染源濃度衰減模型，利用最小二乘法和極大似然預估算法，對污染源定位進行了研究。2009年朱勇等[4]提出一種基于無線傳感器網絡的加權組合三邊算法來對氣體源進行定位。2013年楊一帆等[5]通過構建突發型大氣污染源位置識別的空間反演算法對污染源位置進行了識別。2014年吳正佳等[6]假設各監測點為PM2.5擴散源，求解每個監測點擴散到其他站點的PM2.5濃度，最后根據各監測點位的PM2.5濃度和，定量分析確定污染源位置。上述基于氣體濃度衰減模型的污染源定位算法是在不考慮風向、風速、污染物釋放速率、大氣壓強、溫度和地表狀況等環境因素影響的理想化條件下進行的，然而環境因素對氣體濃度衰減模型的結果會產生很大影響。

為了解決氣體濃度衰減模型在污染源定位中易受環境因素影響的問題，近年來，部分學者將社區網絡分析方法(Social Network Analysis，SNA)應用于研究大氣污染問題。2015年薛安等[7]以161座城市作為節點，以城市間PM2.5質量濃度的相關性與距離的比值作為邊的權重，構建加權網絡，并采用Girvan Newman算法對網絡進行劃分，得到了不同季節中國的污染區域分布情況。2017年孫亞男等[8]通過分析長三角地區的空氣質量指數AQI，證明了大氣污染城市聯動網絡的存在，并從時空的角度揭示了55個城市在大氣污染聯動機制中表現出的內在科學規律。2018年馬博宇等[9]通過構建污染網絡，并根據節點重要性綜合評價結果從整個大網絡中提取出聯系最為緊密、相關性最強的一個霧霾污染網絡。2019年肖琴等[10]通過分析由66座城市組成的污染網絡的社區結構、模體等性質得出了網絡中的主要污染區域。上述利用社區網絡的方法對大氣污染的研究僅停留在發現污染區域的層面，對具體污染源的定位研究還相對較少。

針對上述問題，本文提出了基于社區網絡分析的大氣污染源定位算法。算法首先采用Granger因果檢驗方法[11]分析各監測點處AQI的時間序列，得出任意兩個監測點的AQI之間的影響關系，并且以這些影響關系作為兩個監測點之間的關聯關系，監測點作為節點，構建由監測點組成的污染網絡[12]；然后，通過改進的標簽傳播算法(ILPA)將污染網絡劃分成多個污染區域；最后，采用社區網絡分析方法對劃分后的每個污染子網絡進行中心性分析，根據分析結果預測污染源的位置。本文以273個監測點的AQI作為樣本數據，構建了由監測點組成的污染網絡，并通過對網絡的聚類劃分，以及對各子網絡中節點的度中心度、出度和入度的計算與分析，將污染源的定位范圍由污染區域進一步縮小到了具體的監測點附近，提高了定位的精確性，實驗表明本文方法可準確定位污染源位置。

本文算法流程如圖1所示。

圖1 算法流程

1 模型建立

1.1 關聯關系矩陣的構建

從大氣污染的擴散性角度分析，一些先被污染的監測點處的污染物濃度變化很可能是引起其他監測點處污染物濃度變化的原因，即一個監測點處空氣質量的變化可能會影響其他監測點處空氣質量的變化。因此，分析AQI時間序列之間因果關系，可以得出哪些監測點污染的變化引起了另外一些監測點污染的變化。本文通過Granger因果關系檢驗方法來計算任意兩個監測點處AQI變化的因果關系。

Granger因果關系檢驗是用于判斷兩個時間序列之間是否存在因果關系的一種方法。在時間序列的情形下，兩個變量a1和a2之間的Granger因果關系定義如下：

定義若在包含了變量a1和a2的過去信息的條件下，對a2的預測效果要優于只單獨由a2的過去信息對a2進行的預測效果，即變量a1有助于解釋a2的將來的變化，則認為變量a1是a2的Granger原因[13]。

進行Granger因果關系檢驗的一個前提條件是：參與檢驗的時間序列是平穩序列，對非平穩時間序列需要進行差分處理直到時間序列平穩，否則可能會出現虛假回歸問題，因此在進行格蘭杰因果關系檢驗之前首先應對各時間序列的平穩性進行檢驗。

平穩性指的是監測點v處的AQI的時間序列{AQIv,t}的特性不隨時間t的變化而變化，即均值、方差和協方差都是與時間t無關的常數。常用ADF檢驗來分別對各指標序列的平穩性進行檢驗。具體的Granger因果檢驗步驟如下：

構建兩個監測點v1和v2處AQI的關聯關系方程，表示為：

(1)

式中：AQIv1,t和AQIv2,t分別表示監測點v1和監測點v2在t時刻的具有平穩性的AQI的數據值；αj、βj、γj(j=1，2)是待測參數；{εj,t}(j=1，2)為殘差項；m、n、p、q為滯后階數。

根據參數γ1,i和參數γ2,i是否全為0，檢驗結果存在以下四種可能的情況：

(1) 若γ1,i至少有一項不為0，γ2,i全部為0，則表示：{AQIv2,t}是{AQIv1,t}的Granger原因，且{AQIv1,t}不是{AQIv2,t}的Granger原因，兩監測點之間的污染影響情況關聯關系表現為{AQIv2,t}→{AQIv1,t}，即v2監測點處的大氣污染對v1監測點處的大氣污染有單向影響。

(2) 若γ1,i全部為0，γ2,i至少有一項不為0，則表示：{AQIv1,t}是{AQIv2,t}的Granger原因，且{AQIv2,t}不是{AQIv1,t}的Granger原因，兩監測點之間的污染影響情況關聯關系表現為{AQIv1,t}→{AQIv2,t}，即v1監測點處的大氣污染對v2監測點處的大氣污染有單向影響。

(3) 若γ1,i和γ2,i全部滿足至少有一項不為0，則表示：{AQIv1,t}是{AQIv2,t}的Granger原因，且{AQIv2,t}是{AQIv1,t}的Granger原因，兩監測點之間的污染影響情況關聯關系表現為{AQIv1,t}→{AQIv2,t}且{AQIv2,t}→{AQIv2,t}，即v1監測點處的大氣污染和v2監測點處的大氣污染存在有雙向影響。

(4) 若γ1,i和γ2,i全部為0，則表示：{AQIv1,t}不是{AQIv2,t}的Granger原因，且{AQIv2,t}也不是{AQIv1,t}的Granger原因，即v1監測點處的大氣污染和v2監測點處的大氣污染不存在互相影響。

根據Granger檢驗結果，如果v1監測點處的大氣污染對v2監測點處的大氣污染之間的關聯關系表現為{AQIv1,t}→{AQIv2,t}，則在監測點間的污染影響情況鄰接關系矩陣中Rv1,v2的值為1，否則在監測點間的污染影響情況的鄰接關系矩陣中Rv1,v2的值為0。因此可以構建如式(2)所示的監測點之間的污染影響情況的關聯關系矩陣R。

(2)

式中：Ri,j表示第vi個監測點和第vj個監測點之間的關聯關系值，Ri,j等于1或0。規定當i=j時，Ri,j的值為0。

1.2 污染網絡的構建與劃分

如圖2所示，以監測點作為網絡的節點，以監測點間的污染影響情況的鄰接關系矩陣R作為網絡中節點間的關聯關系連邊，構建監測點間的污染網絡結構圖G。

圖2 污染網絡構建過程

一個污染源產生污染的影響范圍是有限的，在這個有限的范圍內，各個監測點的污染變化情況互相影響較大，而不同的污染范圍之間的污染情況的相互影響相對較小。在污染網絡中表現為：在一個污染區域的內部，監測點之間的具有直接的、較強的、相對緊密的關聯關系，而不同的污染區域之間的關聯關系相對較弱。根據這一性質，可利用標簽傳播的思想將污染網絡劃分，如圖3所示。

圖3 污染網絡的劃分圖

在實際情況中，污染網絡結構圖中所涵蓋的污染區域存在兩類情況，如圖4所示，其中，單一污染區域情況指所構建的污染網絡中僅存在一個污染區域；多污染區域情況指所構建的污染網絡中存在多個污染區域。

(a) 單一污染區域 (b) 多污染區域圖4 兩類污染區域圖

傳統的標簽傳播方法(LPA)僅從網絡的結構特征來決定標簽的傳播[14]，在標簽傳播的過程中，節點根據式(3)選擇其鄰居節點中出現次數最多的標簽更新自身的標簽。

(3)

式中：Ni(l)表示節點i的標簽為l的鄰居節點的集合；Ui(l)表示節點i的鄰居節點中標簽l出現的次數；li表示節點i的標簽，此處Ri,j取1。

然而LPA在真實的污染網絡中并不完全適用。如圖5所示，若網絡中有標簽分別為A、B、C的三個節點1、2、3，其中節點2和節點1在現實的情況里距離較近，節點2與節點3的距離較遠。由于距離較遠，節點3處的污染對節點2處的污染的影響可以忽略不計，因此節點2和節點1應該在同一個污染區域，節點3單獨處于一個污染區域。然而，傳統的LPA并沒有考慮實際情況中的距離因素對于污染擴散的影響。因此，更新節點2的標簽時，可選擇標簽A或標簽C中的任意一個，若選擇將節點2的標簽更新為C時，最終將會造成與實際情況不符的污染網絡的劃分，即將節點2和節點3劃分到一個污染區域，從而影響對于真實的污染網絡的劃分的準確度。

圖5 LPA標簽傳播過程示意圖

因此，本文引入污染影響因子I(i,j)來優化標簽的更新過程。針對污染網絡的劃分，本文基于傳統的LPA算法，提出了一種基于污染影響因子的改進ILPA算法。污染影響因子計算式為：

(4)

式中：(xi，yi)和(xj，yj)分別表示節點vi和vj的緯度坐標及經度坐標。

基于污染影響因子I(i,j)優化后的更新節點標簽的計算式為：

(5)

ILPA算法具體步驟如下：

Step1為每一個節點初始化一個唯一的標簽。

Step2隨機對網絡中的節點進行排序。

Step3根據Step 2中的節點排列順序更新節點，根據式(5)在待更新節點的所有鄰居節點中確定更新后的標簽li。當li不唯一時，隨機選擇一個標簽作為更新后的標簽。

Step4當網絡中的所有節點的標簽穩定不變時，算法結束，否則返回Step 3繼續。

1.3 污染子網絡中污染源的判斷

在一個社交網絡中，描述節點特征的屬性一般有出度、入度等。節點的出度表示的是從該節點出發的指向其他節點的有向線段的條數，它在污染監測網絡中描述了該監測點處污染的變化情況影響的監測點數量，若其出度越大，表示該監測點處的污染情況影響的監測點數量越多。節點的入度表示的是從其他節點出發指向該節點的有向線段的條數，它在污染監測網絡中，描述了對監測點處污染產生影響的監測點數量，若其入度越大，表示對該監測點污染產生影響的其他監測點數量越多。此外，在社交網絡分析中常常采用度中心度這一屬性來描述整體網絡中一個節點與網絡中其他節點的關聯程度。其中，度中心度表示的是與該節點直接相連的節點的個數[15]。在污染網絡圖中，度中心度越大，表示該監測點處污染情況與網絡中其他監測點處污染情況的關聯程度越大。

根據網絡中節點的特征屬性，在一個污染區域內污染源的判斷標準如下：若一個監測點的度中心度較大，且該站點的出度大于入度，則該監測點在污染網絡處于中心位置，與其他監測點的關聯程度最大，且該監測點處的污染情況對其他站點處的污染情況的影響程度大于其他站點處的污染情況對該站點處污染情況的影響，那么該監測點附近存在污染源。如圖6所示，節點6的度中心度相對較大，且出度大于入度，因此在節點6的附近存在污染源。

圖6 污染源分析示意圖

2 實 驗

實驗所采用的原始數據為2018年3月273個監測點的污染物數據，該污染數據主要包括站點名稱、經度、緯度、時間和AQI。

2.1 數據預處理

對各監測點每小時監測到的AQI的數據進行規范化處理，得到各監測點的AQI的時間序列：

(name,x,y,AQIt1,AQIt2,…,AQItn，…，AQItN)

其中，name是監測點的名稱；x和y分別是監測點的緯度坐標值和經度坐標值；AQItn是該監測點在tn時刻監測到的AQI的值。

2.2 活躍監測點的篩選

為了降低實驗的復雜度，提高污染源定位的效率，需要將一些明顯不會對其他站點處的污染產生影響的監測點排除，因此需要對規范化后的大氣污染物數據的進行分類處理。

根據下式分別計算n個監測點在各時刻所監測到的AQI的均值AQItE：

(7)

將AQItE與相應時刻第k個監測點監測到的空氣質量指數AQIk,t比較，根據結果將監測點劃分為不活躍監測點和活躍監測點，劃分標準如表1所示。

表1 不活躍與活躍監測點

如圖7所示，在每一個時刻，由于不活躍監測點的污染水平都低于所有監測點的污染的平均水平，所以這些監測點處的污染對其他監測點處的污染影響能力較小，而活躍監測點的污染水平高于所有監測點的污染平均水平，這些監測點對其他監測點處的污染影響能力較大，活躍監測點中存在大氣污染源的可能性較大，因此利用篩選出的活躍監測點的大氣污染物信息的AQI的時間序列數據進行的實驗。

圖7 不同類型監測點AQI對比

實驗在273個監測點中篩選出210個活躍監測點及63個不活躍監測點。

2.3 實驗結果分析

本文對篩選出的210個活躍監測點的AQI時間序列兩兩之間進行Granger因果關系檢驗，其中對式(1)中的回歸項系數的顯著性檢驗的原假設如下：

原假設1：γ1,i全為0。

原假設2：γ2,i全為0。

對于上述假設的檢驗結果，根據顯著性水平α來確定關聯關系。對于原假設1，若式(1)中第一個方程回歸項系數的顯著性檢驗結果P值小于α，則表示在α的顯著性水平下拒絕了原假設1，否則接受原假設。對于原假設2，若式(1)中第二個方程回歸項系數的顯著性檢驗結果P值小于α，則表示在α的顯著性水平下拒絕了原假設2，否則接受原假設。若在α的顯著性水平下拒絕原假設，則關聯關系為1，否則為0。

因此，根據α=5%的顯著性水平確定參與檢驗的兩個監測點之間關聯關系，部分節點的關聯關系矩陣，如表2所示。

表2 關聯關系矩陣表

通過社會網絡分析軟件Ucinet對關聯關系矩陣進行了可視化，得到污染網絡圖，圖8為上述7個節點組成的污染網絡圖，圖9為全部節點組成的污染網絡與局部放大圖。圖中的節點表示監測點，數字表示監測點的編號，有向線段表示監測點間污染的影響關系。

圖8 部分節點污染網絡圖

圖9 污染網絡與局部放大圖

以圖8中節點為例，用ILPA算法對其進行社區劃分，節點的緯度和經度信息如表3所示，節點的更新順序為節點31、節點37、節點67、節點122、節點121、節點187、節點28。

表3 節點經、緯度信息表

續表3

根據式(5)確定各節點更新后的新標簽的計算過程如表4所示。

表4 標簽更新計算過程表

標簽傳播過程如表5所示，根據結果最終的社區劃分結果為{31，37，67}和{28，121，122，187}。

表5 部分節點標簽傳播過程表

根據ILPA對整體污染網絡的劃分結果，污染網絡中含有4個內部聯系緊密污染區域，即整個污染網絡被劃分成4個污染子網絡，此時模塊度的值最大為0.647。圖10至圖13分別表示污染子網絡1至4的分布圖。

圖10 污染子網絡1的分布圖

圖11 污染子網絡2的分布圖

圖12 污染子網絡3的分布圖

圖13 污染子網絡4的分布圖

結合污染監測網絡中各站點的經度和緯度坐標值，對劃分的結果進行分析可知，4個污染子網絡的主要特征呈現區域性分布，統計結果如表6所示。

表6 污染子網絡特征情況表

根據污染源判斷標準，分別對四個子網絡計算其節點的度中心度、出度和入度。表7反映了污染子網絡1中度中心度排名前10位的節點特征情況。根據污染源判斷標準可知，監測點“(28)東海特鋼西南”的度中心度較大且出度明顯大于入度，說明在子網絡1中監測點“(28)東海特鋼西南”處于污染的中心位置，與其他監測點的交互影響最為直接，且該站點處的污染情況對其他站點處的污染情況影響最大，則在監測點“ (28)東海特鋼西南”附近存在大氣污染源。

表7 子網絡1中節點特征情況表

表8反映了污染子網絡2中度中心度排名前10位的節點特征情況。根據污染源判斷標準可知，監測點 “(138)三興化工西北”“(83)冀東水泥三友2號”、“(137)三興化工東南”的度中心度較大且出度明顯大于入度，說明在子網絡2中，上述3個監測點處于污染的中心位置，與其他監測點的交互影響最為直接，且這3站點處的污染情況對其他站點處的污染情況影響最大，則在以上3個監測點附近存在大氣污染源。

表8 子網絡2中節點特征情況表

表9反映了污染子網絡3中度中心度排名前10位的節點特征情況。根據污染源判斷標準可知，監測點“(146)十二中正南”“(61)國豐鋼鐵東南”“(72)環城北路與銀河路交口東南”的度中心度較大且出度明顯大于入度，說明在子網絡3中，上述3個監測點處于污染的中心位置，與其他監測點的交互影響最為直接，且這3站點處的污染情況對其他站點處的污染情況影響最大，則在以上3個監測點附近存在大氣污染源。

表9 子網絡3中節點特征情況表

表10反映了污染子網絡4中度中心度排名前10位的節點特征情況。根據污染源判斷標準可知，監測點“(42)豐潤污水處理廠東南”“(9)曹雪芹西道與G102國道交口”“(127)仁潤印刷東南”“(208)中車西北”的度中心度較大且出度明顯大于入度，說明在子網絡4中，上述4個監測點處于污染的中心位置，與其他監測點的交互影響最為直接，且這4站點處的污染情況對其他站點處的污染情況影響最大，則在以上4個監測點附近存在大氣污染源。

表10 子網絡4中節點特征情況表

為了驗證實驗結果的準確性，查證上述11個監測點處的PM2.5的均值及平均空氣質量等級情況，如表11所示。

表11 11個監測點的PM2.5情況表

結合表11可知，實驗結果較為準確，編號為9、28、42、61、72、83、127、137、138、146和208的11處監測點處的PM2.5值普遍較高，且空氣質量等級多為重度污染，因此這些監測點附近存在污染源，在污染治理的過程中，應該著重在這些監測點附近排查污染原因，并采取相關措施進行治理，從而達到從根本上解決大氣污染問題的目的。

3 結 語

本文針對多發的大氣污染環境問題，提出一種基于社區網絡的大氣污染源定位算法，該算法利用Granger因果檢驗法挖掘各監測點之間污染變化的相互影響，然后將影響關系轉化為監測點網絡結構圖，利用社區分析的方法在構建的監測點網絡結構圖中進行污染區域的劃分和污染源的定位。實驗表明，本文提出的大氣污染源定位算法，不僅能夠在復雜的環境中找到相應的污染區域，而且能夠將污染源的定位范圍進一步精確到具體的監測點附近，從而進行高效、準確的大氣污染源定位。