李欣欣,王 寧,姜秋俚,劉 樞,魏建勛,張 楠

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

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

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

4(阜新市生態(tài)環(huán)境保護服務中心,阜新 123000)

近年來,國內(nèi)外突發(fā)性重金屬水污染事件屢屢發(fā)生,由于其發(fā)生的突然性和危害的嚴重性會威脅到廣大人民的人身和財產(chǎn)安全,給社會經(jīng)濟的發(fā)展帶來不可估量的損失,因此國內(nèi)外對此事件的關(guān)注極高,并不斷地尋求解決辦法[1].在我國水資源最大的特點是地域分布不均,而我國又是人口大國,缺水狀況也是我國一直在解決的問題.我國有跨國界河流40多條,涉及19個主要國家[2],一旦發(fā)生嚴重的事件,重金屬污染物不只對水環(huán)境和人體的危害極大,甚至會影響社會治安和國際友好關(guān)系.經(jīng)調(diào)查分析,突發(fā)性重金屬水污染的發(fā)生主要是由企業(yè)違規(guī)排放廢水和工廠事故泄漏造成的.因此在自動監(jiān)測站監(jiān)測到河水中重金屬污染物超標時,能否快速精確的定位污染源企業(yè)或工廠,對于后續(xù)采取有效措施具有現(xiàn)實意義.本文研究的是由河流附近企業(yè)違規(guī)排放廢水導致的突發(fā)性重金屬水污染事故的追蹤溯源.

水動力學反演法,根據(jù)河流的水文特征選取合適的水動力學方程,再根據(jù)真實的污染物濃度監(jiān)測數(shù)據(jù)、污染源的分布數(shù)據(jù)等,通過構(gòu)造相關(guān)模型并利用相關(guān)方法反演求解得到污染源的排放位置、排放時間和排放質(zhì)量[2].目前,主要應用了遺傳算法、粒子群算法、模擬退火算法等[3,4],但存在容易得到局部最優(yōu)和搜索時間較長的問題.基于此類問題,本文選用改進后的AFSA算法,算法具有更強的跳出局部極值的能力和更快的搜索能力.

1 水動力反演

重金屬污染源的追蹤溯源需要水動力學反演,因此需要選擇合適的水動力方程進行模擬,從而得到水動力模擬庫,建立相關(guān)模型,快速準確的定位重金屬污染源.

1.1 重金屬污染物的時空變化

在水動力學中,模擬污染物在一維河道中的遷移變化規(guī)律可以用以下公式表示:

式(1)中,C表示河流在x斷面在t時刻污染物的濃度(單位是m g/L),Dx表示河流縱向的彌散系數(shù)(單位是m2/min),ux表示河流縱向平均流速(單位m2/min),K表示污染物的衰減系數(shù)(單位是s?1)[5].

本文根據(jù)實際情況研究的是瞬時點源排放重金屬污染物,故假設在河流的一個斷面x處投入質(zhì)量為M的重金屬污染物,此時污染物的衰減系數(shù)K為0,河流的縱向平均流速為ux.此時污染物與斷面的河水混合,那么方程的初始條件和邊界條件是:

式(2)中,污染物的初始濃度為C0=M/Q,M表示瞬時投放到河流中的污染物的質(zhì)量(單位是kg),Q表示單位時間內(nèi)河水的流量(單位是m3/s)[5],A表示河流的斷面面積(單位是m2).

通過函數(shù) δ的特性和L aplace的變換得到式(1)在式(2)條件下的解析解,即一維河流中重金屬污染物的時空變化方程:

式(3)中,C(x,t)表示污染物在下游的污染物濃度分布(單位是mg/L),x0表示污染源的位置(單位是m),t0表示污染源的排放時刻(單位是min)[5],其余同上.

通過式(3)可知,對于突發(fā)性重金屬水污染的追蹤溯源問題,需要確定污染源的排放位置、排放時間和排放量這3個因素.

1.2 水動力情景模擬

如果重金屬水污染事故突然發(fā)生,自動監(jiān)測站報警,此時已知自動監(jiān)測站處的水文信息,包括監(jiān)測站處河流的斷面面積A、河流的縱向彌散系數(shù)Dx、河流的縱向平均流速ux以及監(jiān)測站監(jiān)測到的某段離散時間內(nèi)污染物濃度變化數(shù)據(jù)等.記自動監(jiān)測站處得到的隨時間變化的污染物濃度數(shù)據(jù)為Ci(i=1,2,···,n).

根據(jù)自動監(jiān)測站處的重金屬污染物的污染范圍信息得重金屬污染源的位置范圍,即根據(jù)情景,假設干流自動監(jiān)測站a(xa)從未監(jiān)測到污染物,而在a下游的干流自動監(jiān)測站b(xb)在時刻T第一次報警,監(jiān)測到了污染物,故得到污染源的排放位置范圍是(xa,xb),而污染源的排放時間與自動監(jiān)測站b(xb)監(jiān)測到污染物的時間間隔范圍是(0,tmax),其中tmax可以由以下公式表示:

通過上述過程已知污染源的位置范圍和排放時間范圍,由式(3)進行水動力情景模擬.為了便于情景模擬,可以將式(3)進行變形得到以下公式:

式(5)中,C(x,t)表示污染物在自動監(jiān)測站a(xa)下游的污染物濃度分布,單位是g /L;x表示污染源距離自動監(jiān)測站b(xb)的距離長度,單位是m;t表示污染源在自動監(jiān)測站b(xb)監(jiān) 測到污染物時刻前t分鐘排放的污染物,單位是min,其余同上.

在位置和時間的范圍內(nèi)進行突發(fā)性重金屬污染源的情景模擬.假設在重金屬污染源的參數(shù)是[x′,t′,M′],代入式(5)得到與自動監(jiān)測站b(xb)相同的一段離散時間內(nèi)一系列的污染物的濃度變化數(shù)據(jù),記為Ci′(i=1,2,···,n).

2 溯源模型建立

重金屬污染源的追蹤溯源需要位置、時間和質(zhì)量3個參數(shù),若將重金屬污染源追蹤溯源問題的3個參數(shù)作為一個整體未知參量進行求解,需要構(gòu)造復雜的模型,模型的復雜會出現(xiàn)溯源時間過長,效率比較低等問題,故不能快速準確地找到污染源.因而本文構(gòu)建兩個模型進行重金屬污染物的溯源,將位置和時間這兩個參數(shù)作為一個整體未知參量建立時空溯源模型,將質(zhì)量作為一個未知參量建立污染物排放量模型.

2.1 時空溯源模型

假定重金屬污染物的排放質(zhì)量是M′,將上述的自動監(jiān)測站的離散時間的重金屬污染物濃度變化序列值Ci(i=1,2,···,n)與水文情景模擬得到位置和時間范圍內(nèi)一系列的在相同離散時間內(nèi)的重金屬污染物濃度變化序列值Ci′(i=1,2,···,n)進行擬合度分析,利用相關(guān)系數(shù)R來判斷真實值與模擬值的相關(guān)性.

在時空溯源模型中,未知參數(shù)量是位置和時間,質(zhì)量設為M′.故若情景模擬的重金屬污染源的位置和時間極為接近于真實重金屬污染源的位置和排放時間那么相關(guān)系數(shù)R的值也極為接近1[6].根據(jù)相關(guān)系數(shù)R的值從大到小排列,R越大越大代表是重金屬污染源的可能性越大,排序越靠前.相關(guān)系數(shù)的公式如下:

式中,Ci表示離散時間內(nèi)時刻t的污染物濃度,C表示離散時間內(nèi)污染物濃度的平均值,Ci′表示情景模擬的離散時間內(nèi)時刻t的污染物濃度,表示情景模擬的離散時間內(nèi)污染物濃度的平均值.

根據(jù)相關(guān)系數(shù)R存在最大值的性質(zhì),為了從一系列的[x′,t′,Ci′] (i=1,2,···,n)中篩選出準確的結(jié)果(x′,t′),故構(gòu)建以下目標函數(shù):

函數(shù)的約束條件為:

由式(4)～式(10)構(gòu)成了時空溯源模型,通過時空溯源模型得到重金屬污染源的位置和排放時間,再通過構(gòu)建污染物排放量模型.

2.2 污染物排放量模型

通過時空溯源模型,重金屬污染源的排放位置和排放時間已經(jīng)確定,即可在此模型中作為已知條件,此時重金屬污染物溯源的未知參數(shù)只有污染源排放污染物的質(zhì)量.故由式(11)可以大致確定源強度范圍.

式中,M表示所求重金屬污染物的質(zhì)量(單位是g),表示上述過程中假設的重金屬污染物的質(zhì)量(單位是g)[5],其余同上.

為了確定污染源排放該重金屬的總量,通過構(gòu)建以下函數(shù)求最小值:

函數(shù)的約束條件為:

為了方便求解,同時空溯源模型一樣將目標函數(shù)變?yōu)榍笞畲笾?目標函數(shù)的公式為:

由式(11)～式(14)構(gòu)成了污染物排放量模型,通過模型得到重金屬污染源的排放量.

3 溯源模型求解

在上述過程中,已經(jīng)把突發(fā)性水污染溯源的兩個模型轉(zhuǎn)化成了關(guān)于兩個求最值的優(yōu)化模型,本文選用改進后的AFSA算法進行模型求解.

改進的AFSA 求解時空溯源模型時,選取F=F1作為適應度目標函數(shù),人工魚個體狀態(tài)是X=(x,t),其當前所在位置的食物濃度函數(shù)用Y=F(X)表示[7],求最大值.人工魚群算法求解污染物排放量模型時,選取F=F2作為適應度目標函數(shù),人工魚個體狀態(tài)是X=M,其當前所在位置的食物濃度函數(shù)用Y=F(X)表示[7],求最大值.

3.1 改進步長和視野

步長和視野范圍是算法的關(guān)鍵參數(shù),步長和視野范圍是算法的關(guān)鍵參數(shù),若希望算法的全局搜索能力更強并且收斂迅速,則需要設置較大的視野;若希望算法的局部搜索能力較強,則需要設置較小的視野.步長大則收斂快,但是步長過大會發(fā)生震蕩現(xiàn)象,步長小則收斂慢但是求解精度高[8].故本文參考文獻[9]通過式(15)動態(tài)調(diào)整視野Visual和步長Step增強算法的搜索能力和精確度[9]:

式中,Visual=xmax/4,S tep=visual/8,Visualmin=0.001、S tepmin=0.0002分別為人工魚的視野、步長的最小值,gen表示當前的迭代次數(shù),maxGen表示最大的迭代次數(shù),xmax表示搜索范圍的最大值[9].

3.2 改進覓食行為

在滿足前進條件的前提下,只有當隨機選擇的狀態(tài)比當前的狀態(tài)好,原覓食行為才會選擇向該方向移動一步,這樣存在搜索速度慢的問題.故為了加快搜索速度,當前人工魚i的狀態(tài)為Xi,通過式(16)隨機得到狀態(tài)Xj,如果滿足前進條件Yj>Yi,則將當前人工魚i直接移動到Xj狀態(tài),如果不滿足前進條件則重新按照式(16)隨機得到狀態(tài)Xj.反之,當達到try_number次后仍然不滿足前進條件,則通過式(17)隨機生成狀態(tài)Xj,將當前人工魚i直接移動到Xj狀態(tài).

式中,i,j=1,2,···,fishnum,rand()為[0,1]之間的隨機數(shù).

3.3 改進的AFSA算法流程

圖1所示的算法流程中初始化設置算法參數(shù)包括:確定人工魚群規(guī)模fishnum、迭代次數(shù)gen、最大迭代次數(shù)maxGen、擁擠度因子δ、嘗試次數(shù)try_number、個體間的距離di,j[8].

圖1 改進的AFSA算法流程圖

4 污染源的排查

通過調(diào)查分析,根據(jù)實際情況,通過上述模型得到污染源在一維河道的位置,排放時間以及排放量后篩選得到的涉銅企業(yè)名單,當企業(yè)距離河流越近并且企業(yè)的信譽越低,則是污染源企業(yè)的可能性越高.故用L表示企業(yè)到納污河流的距離,用E表示企業(yè)的信譽度,其等級分為優(yōu)/良(用0表示),差(用1表示).通加入權(quán)值后得到公式:

式中,w1=0.9975、w2=0.0025,Li表示企業(yè)i到河流的距離,Ei表示企業(yè)i的信譽度,Fi表示加入權(quán)值后企業(yè)i為污染源的函數(shù)值.

當上述企業(yè)名單有k個企業(yè)時,通過概率計算,將企業(yè)名單根據(jù)概率從大到小排列,從而得到相關(guān)工作人員的排查名單.概率計算公式為:

5 算法驗證

本文研究的是瞬時點源突然超標排放重金屬銅導致的河流重金屬水污染,并利用已有的GIS進行結(jié)果展示[10].

5.1 背景

以自動監(jiān)測站b處銅突發(fā)污染超標報警為例,2018年的12月16日上午11點,自動監(jiān)測站b第一次銅超標報警,此時監(jiān)測值為7.5 mg/L.而在同一時刻圖2(黑點代表企業(yè),三角形代表涉銅企業(yè))標注的正常運行的自動監(jiān)測站a、自動監(jiān)測站c以及其他在圖2中未標注的正常運行的自動監(jiān)測站并沒有發(fā)生銅突發(fā)污染超標報警,故污染源存在于自動監(jiān)測站a和自動監(jiān)測站b之間.將該時間作為初始監(jiān)測時間,然后每間隔1小時監(jiān)測一次,監(jiān)測到下午14點結(jié)束.自動監(jiān)測站b處的水文信息如表1所示.

圖2 河流地圖

表1 自動監(jiān)測站b的水文信息

5.2 溯源結(jié)果

通過上述兩個模型得到運行結(jié)果圖3和圖4.圖3中的尋優(yōu)過程圖中圓的圓心坐標數(shù)據(jù)為表2中溯源結(jié)果的污染源位置和污染源排放時間,圖4中的尋優(yōu)過程圖中圓的圓心坐標數(shù)據(jù)為表2中溯源結(jié)果的污染源排放量.表2中,污染源位置即河流中污染源的位置距離自動監(jiān)測站b的距離,污染源排放時間即污染源排放時刻與自動監(jiān)測站b第一次報警的時刻的時間差,污染源排放量即污染源的排放質(zhì)量.

圖3 目標函數(shù)1

5.3 溯源企業(yè)分析

通過上述結(jié)果可以看出上述兩個模型能夠快速準確地得到污染源的3個參數(shù),之后根據(jù)這3個參數(shù)并利用GIS 地圖得到圖5.其中標注圓的半徑通過式(20)確定.

其中,Li代表第i個涉銅企業(yè)到其納污河流的距離(i=1,2,···,n).

圖5中紅色圓內(nèi)涉銅污染源企業(yè)有6個,表明這6個企業(yè)可能為污染源企業(yè),再由式(18)和式(19)計算這6個企業(yè)為污染源的概率并按照概率從小到大排列得到概率排查清單,即表3.

圖4 目標函數(shù)2

表2 溯源結(jié)果分析

圖5 溯源后的河流地圖

表3 概率排查清單(單位:%)

當時發(fā)生重金屬銅污染后,相關(guān)工作人員通過走訪排查自動監(jiān)測站b的上游相關(guān)涉銅企業(yè)進行污染源企業(yè)的追蹤溯源,最終確定為企業(yè)6導致此次河流重金屬銅污染事故的發(fā)生,與表3中概率最大的企業(yè)一致.

6 結(jié)論

通過實驗發(fā)現(xiàn)概率最高的企業(yè)6就是人工排查得到的污染源企業(yè).目前已有的研究工作主要是求得污染源的排放位置,排放時間以及污染物排放量這3個參數(shù)后再通過人工排查確定污染源企業(yè),沒有明確的人工排查順序,故排查過程耗時耗力,不能及時的找到污染源企業(yè)從而采取相應的措施處理問題.實際上這3個參數(shù)值總因為各種因素導致存在相應的誤差,本文提出的方法能在水文信息較少的情況下能更加快速準確地求得污染源的3個參數(shù),并且考慮到由于各種因素導致這3個參數(shù)的誤差無法完全排除,構(gòu)造相關(guān)公式并考慮到企業(yè)信譽度因素,最終為相關(guān)工作人員提供排查的概率清單,進而快速準確地找到污染源企業(yè),及時采取措施解決突發(fā)性重金屬銅水污染問題,保護好我們的水環(huán)境.