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

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

1(中國科學(xué)院 沈陽計(jì)算技術(shù)研究所, 沈陽 110168)

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

3(遼寧省沈陽生態(tài)環(huán)境監(jiān)測中心, 沈陽 110167)

4(遼寧省生態(tài)環(huán)境監(jiān)測中心, 沈陽 110165)

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

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

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

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

1 反演情景

對于污染源的參數(shù)反演需要確定污染源和監(jiān)測站之間的污染物擴(kuò)散關(guān)系, 這就需要一種合適的氣體擴(kuò)散模型來將污染物的擴(kuò)散規(guī)律進(jìn)行量化, 進(jìn)而建立目標(biāo)函數(shù), 將空氣污染源的反演問題轉(zhuǎn)化為全局最優(yōu)化問題.

1.1 氣體擴(kuò)散模型

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

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

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

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

1.2 反演情景模擬

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

根據(jù)式(4)可知, 源參數(shù)反演的問題可以轉(zhuǎn)化為對目標(biāo)函數(shù)Ltar進(jìn)行最優(yōu)化的問題. 目標(biāo)函數(shù)Ltar越小, 表明通過反演得到的污染源的相關(guān)參數(shù)越接近真實(shí)污染源參數(shù), 即反演結(jié)果越接近真實(shí)污染源, 反演的效果越好. 上述污染源的相關(guān)參數(shù)包括污染源的橫坐標(biāo)x0、縱坐標(biāo)y0、 源強(qiáng)Q0和 污染源的高度H0.

以沈陽市渾南區(qū)部署的若干空氣質(zhì)量監(jiān)測站的數(shù)據(jù)為本文的數(shù)據(jù)基礎(chǔ). 在t時(shí)刻, 有若干個(gè)監(jiān)測站監(jiān)測到某污染物濃度超標(biāo), 記錄所有監(jiān)測到污染物濃度超標(biāo)的監(jiān)測站的污染物濃度值(i=1,2,···,N). 與此同時(shí),包括風(fēng)速、風(fēng)向等常規(guī)氣象監(jiān)測信息也是已知的, 從而可以確定大氣穩(wěn)定度等級, 然后參照文獻(xiàn)[13]可以確定水平以及縱向的擴(kuò)散系數(shù)σy和 σz.

為了簡化擴(kuò)散模擬場景, 將整個(gè)監(jiān)測的區(qū)域平面抽象成一個(gè)矩形平面, 并建立二維平面坐標(biāo)系, 對每個(gè)監(jiān)測站進(jìn)行坐標(biāo)轉(zhuǎn)化, 每個(gè)監(jiān)測站坐標(biāo)可以以平面坐標(biāo)(xi,yi)(i=1,2,···,N)來表示. 將監(jiān)測到污染物濃度超標(biāo)的監(jiān)測站的坐標(biāo)、監(jiān)測的污染物濃度和當(dāng)前的平均風(fēng)速代入反演模型中, 即可進(jìn)行污染源的反演, 最終得到污染源的4個(gè)參數(shù)結(jié)果.

2 反演模型求解

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

2.1 傳統(tǒng)蟻群算法

傳統(tǒng)的ACO算法是是由意大利學(xué)者Dorigo等[14]提出的一種人工仿生智能優(yōu)化算法, 螞蟻在覓食過程中, 會(huì)在經(jīng)過的路徑留下信息素, 信息素濃度高的路徑會(huì)更大概率被其他螞蟻選擇, 整個(gè)種群會(huì)根據(jù)信息素的交流和傳遞來進(jìn)行不斷調(diào)整下一步的移動(dòng)方向, 進(jìn)而找到最佳的覓食路徑.

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

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

傳統(tǒng)的蟻群算法容易陷入局部極值是因?yàn)槠湔答仚C(jī)制會(huì)導(dǎo)致如果當(dāng)前解不是最優(yōu)解時(shí), 依舊會(huì)吸引其他螞蟻更大概率選擇當(dāng)前解, 從而逐漸喪失了種群的多樣性, 導(dǎo)致算法陷入局部極值. 因此, 本文類比遺傳算法中的選擇淘汰和交叉操作, 在每次所有螞蟻當(dāng)前這代結(jié)果求解完成后, 對整個(gè)蟻群按照目標(biāo)函數(shù)值從小到大進(jìn)行排序, 將值最大的1/4個(gè)體舍去, 從剩下的3/4個(gè)體中隨機(jī)選擇值處于中間位置的1/3個(gè)體分為兩個(gè)部分, 分別從兩個(gè)部分中隨機(jī)選擇兩個(gè)個(gè)體P1和P2進(jìn)行交叉操作, 參照遺傳算法中的部分映射交叉形式將P1和P2的一部分進(jìn)行交換, 從而形成兩個(gè)新的個(gè)體, 反復(fù)操作直到新個(gè)體的數(shù)量達(dá)到原種群數(shù)量的1/4, 然后填充到之前舍去的1/4個(gè)體處, 達(dá)到維持種群多樣性的目的, 具體交叉方式可如下表示:

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

2.3 提出獎(jiǎng)懲因子機(jī)制加快收斂速度

傳統(tǒng)蟻群算法中, 螞蟻的信息素越大, 表明該螞蟻代表的可行解越優(yōu)質(zhì), 在后續(xù)的迭代過程中能夠吸引劣質(zhì)螞蟻向其靠近. 在每次迭代結(jié)束后, 每只螞蟻都會(huì)根據(jù)解的優(yōu)劣程度來更新自己的信息素, 信息素更新機(jī)制可如下表示:

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

由式(6)可知, 在傳統(tǒng)蟻群算法中, 每只螞蟻的信息素更新雖然會(huì)根據(jù)其解的優(yōu)劣程度來區(qū)分更新, 但是在解的優(yōu)劣程度相差不明顯的時(shí)候, 較優(yōu)解被發(fā)現(xiàn)的概率下降, 劣質(zhì)螞蟻無法向優(yōu)質(zhì)螞蟻進(jìn)行靠攏, 使得收斂速度降低. 為了加快收斂速度, 本文對信息素更新機(jī)制進(jìn)行改進(jìn), 在式(7)的基礎(chǔ)上, 提出獎(jiǎng)懲因子機(jī)制m, 對信息素增量的計(jì)算進(jìn)行修改, 可表示為:

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

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

2.4 M-ACO算法流程

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

圖1 M-ACO算法流程圖

3 實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)例背景

以沈陽市渾南區(qū)某區(qū)域內(nèi)在2019年某日因某工廠違規(guī)排放引起的顆粒物濃度超標(biāo)事件為實(shí)驗(yàn)實(shí)例.該區(qū)域內(nèi)存在小范圍部署的11個(gè)空氣監(jiān)測站, 其實(shí)際部署情況如圖2所示.

圖2 監(jiān)測站點(diǎn)位實(shí)際部署情況

此次污染物濃度超標(biāo)事件一共使得6個(gè)監(jiān)測站監(jiān)測到顆粒物濃度超標(biāo), 根據(jù)當(dāng)前其他常規(guī)氣象參數(shù)確定大氣穩(wěn)定度等級為E. 那么根據(jù)文獻(xiàn)[13]可以得知該時(shí)刻的σy和 σz的取值如式(8)所示:

已知風(fēng)向?yàn)闁|偏北風(fēng), 為了方便確定污染源范圍,另給出上風(fēng)向處一個(gè)未監(jiān)測到顆粒物濃度超標(biāo)的監(jiān)測站作為參考監(jiān)測站.

為了方便實(shí)驗(yàn), 根據(jù)上述11個(gè)常規(guī)監(jiān)測站以及另給出的1個(gè)參考監(jiān)測站的實(shí)際部署情況, 以地圖的某個(gè)點(diǎn)為原點(diǎn), 按照監(jiān)測站之間的實(shí)際距離等比例進(jìn)行描點(diǎn), 將其抽象成為二維平面圖, 并且用圓圈框住上述監(jiān)測到顆粒物濃度超標(biāo)的6個(gè)監(jiān)測站, 則當(dāng)前整個(gè)區(qū)域內(nèi)所有監(jiān)測站的狀態(tài)如圖3所示.

圖3 顆粒物濃度超標(biāo)時(shí)各監(jiān)測站情況

將圖3圓圈內(nèi)的6個(gè)監(jiān)測站的坐標(biāo)以及此時(shí)的顆粒物監(jiān)測濃度記錄下來, 并將其代入式(4)中, 用來進(jìn)行污染源的反演.

3.2 反演結(jié)果

通過上述實(shí)驗(yàn)過程, 利用Matlab軟件編寫ACO算法和M-ACO算法分別進(jìn)行求解. 經(jīng)過反復(fù)調(diào)試, 結(jié)合程序的整體運(yùn)行效率和運(yùn)算結(jié)果的準(zhǔn)確性進(jìn)行綜合考慮, 最終決定設(shè)置迭代次數(shù)為1 000, 蟻群數(shù)量為200, 信息素?fù)]發(fā)因子為0.8, 信息素重要程度因子為1,啟發(fā)函數(shù)重要程度因子為5. 最終可以得到兩種算法運(yùn)算的目標(biāo)函數(shù)值對比變化情況如圖4所示.

由圖4可以看出, 目標(biāo)函數(shù)在ACO算法中在615代后達(dá)到收斂狀態(tài), 并且此收斂狀態(tài)并不是全局最優(yōu)狀態(tài); 而在M-ACO算法中, 142代時(shí), 目標(biāo)函數(shù)的收斂速度急劇下降, 目標(biāo)函數(shù)即將陷入局部極值, 此時(shí)通過不斷地交叉操作擴(kuò)大搜索范圍, 在243代后進(jìn)一步收斂, 在336代后基本達(dá)到全局最優(yōu).

圖4 目標(biāo)函數(shù)值對比變化情況

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

3.3 結(jié)果評估

通過對該次顆粒物濃度超標(biāo)事件的歷史記錄和相關(guān)資料的進(jìn)一步調(diào)查了解, 最終獲得了真實(shí)污染源的4項(xiàng)相關(guān)參數(shù)數(shù)據(jù), 如表1所示.

表1 真實(shí)污染源相關(guān)參數(shù)信息

此時(shí), 根據(jù)兩種算法分別進(jìn)行求解得到的反演源的4項(xiàng)參數(shù)值的變化情況, 同表1中真實(shí)污染源進(jìn)行對比, 可以繪制出反演源中4項(xiàng)參數(shù)值的變化趨勢, 進(jìn)而直觀地對比兩種算法尋優(yōu)結(jié)果的準(zhǔn)確性和效率.

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

圖5 反演源各項(xiàng)參數(shù)對比變化情況

同時(shí)圖5也驗(yàn)證了使用M-ACO算法在本次顆粒物濃度超標(biāo)事件中, 可以快速準(zhǔn)確地求出污染源的位置, 源強(qiáng)和排放高度等信息, 達(dá)到了空氣污染源反演的目的, 證明了該算法在之后的污染源反演工作中可以進(jìn)行實(shí)際應(yīng)用.

由上述實(shí)驗(yàn)求得的每個(gè)參數(shù)的最終結(jié)果如表2所示.

表2 兩種算法反演結(jié)果對比

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

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

表3 反演結(jié)果相對誤差率對比

4 結(jié)論

通過上述實(shí)例, 可以發(fā)現(xiàn)M-ACO算法在污染源反演的準(zhǔn)確性和效率上都比傳統(tǒng)的ACO算法要優(yōu)秀,能夠把單個(gè)空氣污染源反演結(jié)果的相對誤差率控制到2%以內(nèi); 同時(shí)解決了傳統(tǒng)ACO算法因容易陷入局部極值而無法保證結(jié)果準(zhǔn)確性的問題, 為實(shí)際的污染源定位排查工作提供有力的數(shù)據(jù)支持, 為空氣污染的治理起到了關(guān)鍵的引導(dǎo)性作用, 促進(jìn)了空氣污染治理由憑感覺、憑經(jīng)驗(yàn)管理轉(zhuǎn)向精準(zhǔn)化、靶向化管理, 保護(hù)了我們的空氣環(huán)境.