基于改進AFSA算法的河流突發水污染溯源①

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

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

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

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

4(阜新市生態環境保護服務中心,阜新 123000)

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

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

1 水動力反演

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

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

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

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

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

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

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

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

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

1.2 水動力情景模擬

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

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

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

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

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

2 溯源模型建立

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

2.1 時空溯源模型

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

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

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

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

函數的約束條件為:

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

2.2 污染物排放量模型

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

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

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

函數的約束條件為:

為了方便求解,同時空溯源模型一樣將目標函數變為求最大值,目標函數的公式為:

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

3 溯源模型求解

在上述過程中,已經把突發性水污染溯源的兩個模型轉化成了關于兩個求最值的優化模型,本文選用改進后的AFSA算法進行模型求解.

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

3.1 改進步長和視野

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

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

3.2 改進覓食行為

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

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

3.3 改進的AFSA算法流程

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

圖1 改進的AFSA算法流程圖

4 污染源的排查

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

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

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

5 算法驗證

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

5.1 背景

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

圖2 河流地圖

表1 自動監測站b的水文信息

5.2 溯源結果

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

圖3 目標函數1

5.3 溯源企業分析

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

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

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

圖4 目標函數2

表2 溯源結果分析

圖5 溯源后的河流地圖

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

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

6 結論

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