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

1(中國(guó)科學(xué)院大學(xué),北京 100049)

2(中國(guó)科學(xué)院 沈陽(yáng)計(jì)算技術(shù)研究所,沈陽(yáng) 110168)

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

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

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

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

1 水動(dòng)力反演

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

1.1 重金屬污染物的時(shí)空變化

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

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

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

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

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

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

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

1.2 水動(dòng)力情景模擬

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

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

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

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

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

2 溯源模型建立

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

2.1 時(shí)空溯源模型

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

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

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

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

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

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

2.2 污染物排放量模型

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

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

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

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

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

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

3 溯源模型求解

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

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

3.1 改進(jìn)步長(zhǎng)和視野

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

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

3.2 改進(jìn)覓食行為

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

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

3.3 改進(jìn)的AFSA算法流程

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

圖1 改進(jìn)的AFSA算法流程圖

4 污染源的排查

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

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

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

5 算法驗(yàn)證

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

5.1 背景

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

圖2 河流地圖

表1 自動(dòng)監(jiān)測(cè)站b的水文信息

5.2 溯源結(jié)果

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

圖3 目標(biāo)函數(shù)1

5.3 溯源企業(yè)分析

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

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

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

圖4 目標(biāo)函數(shù)2

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

圖5 溯源后的河流地圖

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

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

6 結(jié)論

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