基于M-ACO算法的空氣污染源反演①

祁柏林, 崔英杰, 王 帥, 武 暕

1(中國科學院 沈陽計算技術研究所, 沈陽 110168)

2(中國科學院大學, 北京 100049)

3(遼寧省沈陽生態環境監測中心, 沈陽 110167)

4(遼寧省生態環境監測中心, 沈陽 110165)

隨著工業的發展, 空氣污染已經給人民群眾的生活質量和我國經濟的快速發展造成了很大的影響. 對于人類的影響包括造成了一系列疾病, 甚至縮短了人類的預期壽命[1,2]. 對于城市經濟發展的影響包括削弱了城市的集聚效應與規模, 減緩了城市的生產活動, 進而抑制了城市化進程[3]. 在空氣污染超標事故發生時首要解決的就是如何尋找污染源. 傳統的對于空氣污染源反演的方法主要是通過對企業的源排放清單進行排查. 但在源清單的獲取過程中, 由于監測區域空氣污染物的開放性, 使得排放清單的獲取難度過大, 數據不夠完全. 即便是能夠獲取到完備的數據, 整個污染源清單的數據量也會過于龐大, 從而導致單純通過人力排查十分困難且低效. 為了解決源清單排查法具有的低效和盲目性的問題, 國內外學者紛紛開始尋找更優質的方法來應用到污染源的反演工作中.

根據近些年來國內外學者對于污染源反演工作的相關研究, 可以把污染源反演的方法劃分為兩類: 基于數理概率統計理論的反演方法和基于全局優化搜索的方法. (1)基于數理統計理論的方法是利用概率論的推理來進行污染源問題的描述, 根據推理結果的分布情況來反映污染源的概率統計[4–6]. 鄭情文[5]利用伴隨概率定位法的計算策略, 以2 km×2 km的理想開放空間為研究區域, 理想化風向和風速, 利用風洞模型實現了污染源的追蹤. 這種方法對于高維問題無法快速降維導致求解時間過長. (2)基于全局優化搜索的方法是根據氣體擴散的特征、氣象等條件, 利用氣體擴散模型,建立目標函數, 采取啟發式搜索算法對目標函數進行求解, 從而得到污染源的相關參數[7–11]. 沈澤亞等[7]通過對遺傳算法、粒子群算法和單純形算法的兩兩耦合,分別對污染源進行了反演, 最終得到遺傳算法和單純形算法耦合時, 污染源反演的效率最高, 但準確性無法保證, 粒子群與單純形算法耦合時, 污染源反演的準確性最高, 但效率無法保證; Thomson等[8]采用模擬退火算法, 結合隨機搜索算法對污染源進行了反演, 其反演效率較低; Ma等[9]提出了一種將粒子群算法同吉洪諾夫規則耦合的污染源反演方法, 但是目前無法在實際中投入使用.

相比于其他群智能算法, 蟻群算法采用的是正反饋并行機制, 具有魯棒性強, 對初值依賴性不大等優點,在最短路徑規劃問題上得到了很廣泛的應用. 本文創新性地將蟻群算法應用到空氣污染源反演的領域中來.

雖然蟻群算法相比于其他群智能算法存在著一些優勢, 但仍然存在著一些缺陷, 包括初期收斂速度較慢以及容易陷入局部最優值. 本文將針對這兩個缺陷進行改進, 從而對空氣污染源反演的準確性和高效性進行保證.

1 反演情景

對于污染源的參數反演需要確定污染源和監測站之間的污染物擴散關系, 這就需要一種合適的氣體擴散模型來將污染物的擴散規律進行量化, 進而建立目標函數, 將空氣污染源的反演問題轉化為全局最優化問題.

1.1 氣體擴散模型

為了和先前學者的工作保持一定的延續性, 在污染源的下風向處, 可以利用點源連續擴散的高斯煙羽模型模擬出空氣在三維空間中的擴散濃度變化情況[12].其表現形式如下:

其中,C(x,y,z)表示下風向處任意一點的氣體模擬濃度( g/m3);Q0表示污染源的強度 (g /s );u表示平均風速(m /s) ; σy和 σz分別表示平面和豎直垂直于下風向的擴散系數;x0和y0分 別表示污染源的平面坐標;x和y分別表示監測站的平面坐標;z表示污染源下風處任意一點的垂直坐標;H0表示污染源的高度. 在實際的空氣污染監測中, 我們一般最關注的焦點是污染物在近地面的濃度分布狀況, 因此本文中忽略所有的監測站的高度,也就是令z=0, 可以得到:

式(2)即為我們建立污染源反演模型的基礎, σy和σz分別可由Pasquill-Gifford模型擴散系數方程確定.其一般表現形式為:

其中, γ1, λ1, γ2, λ2均由不同環境以及大氣穩定度等級來確定[13].

1.2 反演情景模擬

在式(2)的基礎上, 假定和分別為第個監測站的實際監測污染物濃度和通過式(2)計算模擬到的第i個監測站所處位置的污染物濃度, 那么以此可以建立污染源參數反演的目標函數, 以平方損失函數為基礎, 對每個監測站的監測濃度進行加和, 建立目標函數即反演模型[7,9,10]表現形式如下:

根據式(4)可知, 源參數反演的問題可以轉化為對目標函數Ltar進行最優化的問題. 目標函數Ltar越小, 表明通過反演得到的污染源的相關參數越接近真實污染源參數, 即反演結果越接近真實污染源, 反演的效果越好. 上述污染源的相關參數包括污染源的橫坐標x0、縱坐標y0、 源強Q0和 污染源的高度H0.

以沈陽市渾南區部署的若干空氣質量監測站的數據為本文的數據基礎. 在t時刻, 有若干個監測站監測到某污染物濃度超標, 記錄所有監測到污染物濃度超標的監測站的污染物濃度值(i=1,2,···,N). 與此同時,包括風速、風向等常規氣象監測信息也是已知的, 從而可以確定大氣穩定度等級, 然后參照文獻[13]可以確定水平以及縱向的擴散系數σy和 σz.

為了簡化擴散模擬場景, 將整個監測的區域平面抽象成一個矩形平面, 并建立二維平面坐標系, 對每個監測站進行坐標轉化, 每個監測站坐標可以以平面坐標(xi,yi)(i=1,2,···,N)來表示. 將監測到污染物濃度超標的監測站的坐標、監測的污染物濃度和當前的平均風速代入反演模型中, 即可進行污染源的反演, 最終得到污染源的4個參數結果.

2 反演模型求解

對于上述最優化問題, 本文采用改進后的ACO算法進行迭代求解, 即選取式(4)作為目標函數, 通過迭代求其最小值, 最終得到污染源的參數, 包括源強Q0,源坐標x0和y0以 及源釋放高度H0.

2.1 傳統蟻群算法

傳統的ACO算法是是由意大利學者Dorigo等[14]提出的一種人工仿生智能優化算法, 螞蟻在覓食過程中, 會在經過的路徑留下信息素, 信息素濃度高的路徑會更大概率被其他螞蟻選擇, 整個種群會根據信息素的交流和傳遞來進行不斷調整下一步的移動方向, 進而找到最佳的覓食路徑.

然而傳統ACO算法存在著容易陷入局部極值和收斂速度較慢的缺點. 本文設計一種機制用來對基本蟻群算法的正反饋機制進行制衡, 使其能夠跳出局部極值. 同時對信息素更新機制進行調整, 使其在迭代初期能夠進行更快速地收斂, 最終將其歸納為M-ACO算法.

2.2 引入選擇交叉思想避免陷入局部極值

傳統的蟻群算法容易陷入局部極值是因為其正反饋機制會導致如果當前解不是最優解時, 依舊會吸引其他螞蟻更大概率選擇當前解, 從而逐漸喪失了種群的多樣性, 導致算法陷入局部極值. 因此, 本文類比遺傳算法中的選擇淘汰和交叉操作, 在每次所有螞蟻當前這代結果求解完成后, 對整個蟻群按照目標函數值從小到大進行排序, 將值最大的1/4個體舍去, 從剩下的3/4個體中隨機選擇值處于中間位置的1/3個體分為兩個部分, 分別從兩個部分中隨機選擇兩個個體P1和P2進行交叉操作, 參照遺傳算法中的部分映射交叉形式將P1和P2的一部分進行交換, 從而形成兩個新的個體, 反復操作直到新個體的數量達到原種群數量的1/4, 然后填充到之前舍去的1/4個體處, 達到維持種群多樣性的目的, 具體交叉方式可如下表示:

其中,P1、P2表示所選螞蟻個體, 包含反演源的解組成的矩陣,表示交叉后形成新的螞蟻個體,A和B表示交叉點, 其取值范圍為(0, 1)且A＜B.

2.3 提出獎懲因子機制加快收斂速度

傳統蟻群算法中, 螞蟻的信息素越大, 表明該螞蟻代表的可行解越優質, 在后續的迭代過程中能夠吸引劣質螞蟻向其靠近. 在每次迭代結束后, 每只螞蟻都會根據解的優劣程度來更新自己的信息素, 信息素更新機制可如下表示:

其中,n為當前迭代的次數; ρ表示信息素揮發因子;τi(n) 為第i只螞蟻在n次迭代開始時包含的信息素,τi(n+1) 為第i只螞蟻在n次迭代后更新完畢的信息素,同時也作為其在第n+1次迭代開始時的信息素;Δτi(t)為第i 只螞蟻在第n 次迭代過程中獲得的信息素增量, 其表現形式如下:

由式(6)可知, 在傳統蟻群算法中, 每只螞蟻的信息素更新雖然會根據其解的優劣程度來區分更新, 但是在解的優劣程度相差不明顯的時候, 較優解被發現的概率下降, 劣質螞蟻無法向優質螞蟻進行靠攏, 使得收斂速度降低. 為了加快收斂速度, 本文對信息素更新機制進行改進, 在式(7)的基礎上, 提出獎懲因子機制m, 對信息素增量的計算進行修改, 可表示為:

其中, n 為當前迭代的次數, Δ τi(t)為 第i 只 螞蟻在第n 次迭代過程中獲得的信息素增量; μ為獎懲系數, 其值由式(9)決定, 對于第i 只螞蟻, 如果其n -1代的目標函數值n-1) 大于其第n 代的目標函數值n), 那么m則會大于1; 反之, m 則會小于1.

式(8)和式(9)共同決定了第i只螞蟻在第n 次迭代過程中能獲得的信息素增量Δ τi(n). 式(8)表示了目標函數值n)越小的螞蟻可以獲得越多的信息素增量,式(9)通過第i只螞蟻相鄰兩代之間的目標函數值的比值, 即獎懲系數μ 來動態地放大或者縮小每只螞蟻各自對應的信息素增量. 這樣會使得螞蟻在根據式(6)更新各自的信息素后, 形成兩極分化的狀態, 優質螞蟻信息素更多而劣質螞蟻信息素更少. 優質螞蟻在后續的迭代過程中能以更大的概率被發現, 吸引劣質螞蟻向其靠攏, 從而加快了算法的收斂速度.

2.4 M-ACO算法流程

M-ACO算法流程如圖1所示, 其與傳統蟻群算法流程中最大的區別在于: (1)每次迭代完畢后都需要進行選擇淘汰交叉來維持種群多樣性; (2)每次更新信息素時需要利用獎懲因子對每只螞蟻的信息素增量進行放縮.整個算法過程中需要初始化的算法參數包括: 蟻群規模 m , 迭代次數T , 信息素揮發因子ρ, 信息素和啟發函數的重要程度因子α 和β.

圖1 M-ACO算法流程圖

3 實驗驗證

3.1 實例背景

以沈陽市渾南區某區域內在2019年某日因某工廠違規排放引起的顆粒物濃度超標事件為實驗實例.該區域內存在小范圍部署的11個空氣監測站, 其實際部署情況如圖2所示.

圖2 監測站點位實際部署情況

此次污染物濃度超標事件一共使得6個監測站監測到顆粒物濃度超標, 根據當前其他常規氣象參數確定大氣穩定度等級為E. 那么根據文獻[13]可以得知該時刻的σy和 σz的取值如式(8)所示:

已知風向為東偏北風, 為了方便確定污染源范圍,另給出上風向處一個未監測到顆粒物濃度超標的監測站作為參考監測站.

為了方便實驗, 根據上述11個常規監測站以及另給出的1個參考監測站的實際部署情況, 以地圖的某個點為原點, 按照監測站之間的實際距離等比例進行描點, 將其抽象成為二維平面圖, 并且用圓圈框住上述監測到顆粒物濃度超標的6個監測站, 則當前整個區域內所有監測站的狀態如圖3所示.

圖3 顆粒物濃度超標時各監測站情況

將圖3圓圈內的6個監測站的坐標以及此時的顆粒物監測濃度記錄下來, 并將其代入式(4)中, 用來進行污染源的反演.

3.2 反演結果

通過上述實驗過程, 利用Matlab軟件編寫ACO算法和M-ACO算法分別進行求解. 經過反復調試, 結合程序的整體運行效率和運算結果的準確性進行綜合考慮, 最終決定設置迭代次數為1 000, 蟻群數量為200, 信息素揮發因子為0.8, 信息素重要程度因子為1,啟發函數重要程度因子為5. 最終可以得到兩種算法運算的目標函數值對比變化情況如圖4所示.

由圖4可以看出, 目標函數在ACO算法中在615代后達到收斂狀態, 并且此收斂狀態并不是全局最優狀態; 而在M-ACO算法中, 142代時, 目標函數的收斂速度急劇下降, 目標函數即將陷入局部極值, 此時通過不斷地交叉操作擴大搜索范圍, 在243代后進一步收斂, 在336代后基本達到全局最優.

圖4 目標函數值對比變化情況

從收斂速率上來看, 目標函數在M-ACO算法的迭代過程中, 收斂效率比ACO算法要高許多, 這是由于利用獎懲因子機制更新信息素的策略, 使得優質螞蟻個體能夠積累更多的信息素, 從而放大了其被選擇的概率; 而在ACO算法中, 由于各個螞蟻之間的信息素濃度差異不夠明顯, 導致優質螞蟻無法被選擇到, 因此收斂速度較慢.

3.3 結果評估

通過對該次顆粒物濃度超標事件的歷史記錄和相關資料的進一步調查了解, 最終獲得了真實污染源的4項相關參數數據, 如表1所示.

表1 真實污染源相關參數信息

此時, 根據兩種算法分別進行求解得到的反演源的4項參數值的變化情況, 同表1中真實污染源進行對比, 可以繪制出反演源中4項參數值的變化趨勢, 進而直觀地對比兩種算法尋優結果的準確性和效率.

由圖5可知, 在使用ACO算法迭代過程中, 目標函數收斂速度慢, 并且最終陷入了局部極值, 最終得到的反演源的4項參數結果與真實源4項參數結果偏差較大, 無法保證空氣污染源反演的準確性和效率; 而在M-ACO算法迭代過程中, 最終得到的反演源4項參數結果更準確. 更接近真實污染源的4項參數值.

圖5 反演源各項參數對比變化情況

同時圖5也驗證了使用M-ACO算法在本次顆粒物濃度超標事件中, 可以快速準確地求出污染源的位置, 源強和排放高度等信息, 達到了空氣污染源反演的目的, 證明了該算法在之后的污染源反演工作中可以進行實際應用.

由上述實驗求得的每個參數的最終結果如表2所示.

表2 兩種算法反演結果對比

根據表2中分別利用兩種算法得到的反演源的4項參數值, 與真實源的4項參數值做對比, 計算出相對誤差率, 如表3所示.

從表3可以看出, 使用M-ACO算法得到的反演結果的誤差率明顯低于使用ACO算法的誤差率, 污染源反演的準確性得到了保證.

表3 反演結果相對誤差率對比

4 結論

通過上述實例, 可以發現M-ACO算法在污染源反演的準確性和效率上都比傳統的ACO算法要優秀,能夠把單個空氣污染源反演結果的相對誤差率控制到2%以內; 同時解決了傳統ACO算法因容易陷入局部極值而無法保證結果準確性的問題, 為實際的污染源定位排查工作提供有力的數據支持, 為空氣污染的治理起到了關鍵的引導性作用, 促進了空氣污染治理由憑感覺、憑經驗管理轉向精準化、靶向化管理, 保護了我們的空氣環境.