胡玉珊，丁 濤，張舒羽，顏丹佳

(1. 浙江同濟科技職業(yè)學(xué)院，浙江 杭州 311231；2. 中國計量大學(xué) 質(zhì)量與安全工程學(xué)院，浙江 杭州 310018)

面對頻繁發(fā)生的水環(huán)境污染事件，如何快速和精準地進行污染物溯源，盡快獲知引發(fā)事故的污染源信息(排放位置、排放時間和排放量)，對實現(xiàn)水污染事故的靶向治理和制定有效的應(yīng)急處置方案具有重要意義。近幾十年來，在突發(fā)流域性水污染事故時，水域的水質(zhì)得不到保障[1]，國內(nèi)外已有很多學(xué)者在研究未知水體污染事件的溯源及污染源識別方面取得了一定的成果。這些研究大部分圍繞確定性方法和不確定性方法的思路對河流水環(huán)境污染事件進行溯源研究與討論[2-3]。水污染事故污染源識別的確定性理論方法包括直接求解法和啟發(fā)式優(yōu)化法。直接求解法是指將參數(shù)作為因變量，利用空間域和時間域上的導(dǎo)數(shù)信息，應(yīng)用基于水質(zhì)-水動力模型的反向計算法進行污染源信息求解的一種方法[4-10]。

水污染事故污染源識別的不確定性理論方法常采用隨機方法。隨機方法用概率分布函數(shù)來描述客觀事故的隨機性，包括統(tǒng)計歸納法[11]、貝葉斯推理法[12]以及最小相對熵(MRE)等。統(tǒng)計歸納法的優(yōu)點在于能基于大量數(shù)據(jù)作不確定性分析，如Huang等[13]在收集錢塘江大量監(jiān)測點數(shù)據(jù)的基礎(chǔ)上，通過統(tǒng)計歸納分析得到錢塘江污染物來源的重點區(qū)域。最小相對熵(MRE)的優(yōu)點在于對待求問題進行不確定性分析，即它能基于待求問題的先驗分布獲取待求問題的二次估計，如Woodbury等[14]率先將MRE應(yīng)用于追溯地下水污染源信息以及重構(gòu)大氣污染歷史。許多學(xué)者在基于傳統(tǒng)貝葉斯推理的基礎(chǔ)上加入MCMC抽樣算法，通過以樣本頻率代替概率的方法得到污染源信息，通常稱這種方法為Bayesian inference-MCMC算法。如Jiang等[15]通過直接MC抽樣法，在反演污染源信息的基礎(chǔ)上和溯源過程中公認的不確定性分析框架下，采用貝葉斯推理方法對擴散系數(shù)、流速等不確定性參數(shù)進行了敏感度分析及反演。Wei等[16]運用Bayesian-MCMC抽樣算法，在污染源溯源過程中考慮算法中敏感參數(shù)的不確定性，消除了算法中敏感參數(shù)所帶來的誤差，得到了更加準確的污染源溯源結(jié)果。

從上述文獻可知，目前水污染事故溯源研究中無論是確定性方法還是不確定性方法，其研究對象主要針對突發(fā)水污染事件的瞬時污染源，鮮見對持續(xù)排放的連續(xù)性污染源溯源方法的研究工作。而在現(xiàn)實環(huán)境中，水環(huán)境污染事件可分為兩種，一種是瞬時污染事故，即在極短時間內(nèi)通過大口徑排放管道把大量污水排放到自然水體中，如運輸化學(xué)品車輛傾覆或企業(yè)污水未經(jīng)處理短時間大流量地直接排入自然水體；另一種是持續(xù)排放一段時間的污染事故，如化工企業(yè)在夜間的持續(xù)偷排或管道破裂引起的一段時間小流量污水泄漏事件等。因此，筆者針對連續(xù)排放的水污染事故溯源問題，基于二維對流擴散積分模型，研究提出了水污染溯源方法的理論公式，彌補和擴展了現(xiàn)有水污染事故溯源的方法體系。為使該理論公式能夠應(yīng)用于水污染溯源實踐，同時研究構(gòu)建了兩點和四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)的布設(shè)方式，通過構(gòu)建的應(yīng)急監(jiān)測網(wǎng)絡(luò)來監(jiān)測和計算目標污染物濃度在空間域和時間域上的差分信息，得到優(yōu)化布設(shè)方案，從而解決溯源方法向應(yīng)用轉(zhuǎn)化的關(guān)鍵問題。

1 二維對流擴散連續(xù)源積分模型

連續(xù)排放污染源溯源算法基于二維對流擴散連續(xù)源積分模型。在原點(0,0)處從t=0開始以速率為Q連續(xù)排放污染物，濃度為C0，因為是連續(xù)排放，可以視為一系列瞬時排放的疊加。

設(shè)在原點(0,0)處于t=τ時刻，經(jīng)過dτ，注入的污染物排放量為dM=C0Qdτ，在(x,y)處，t時刻的污染物濃度為

(1)

式(1)中：C(x,y,t)為河流點的污染物濃度,mg/L；C0為連續(xù)源排放濃度,mg/L；Q為連續(xù)源排放速率,L/s；Dx，Dy分別為x，y方向湍流擴散系數(shù),m2/s；t為時間,s；x為與排放污染源的水平距離,m；y為與排放污染源的垂直距離,m；u，v分別為水流在x，y方向的速度分量,m/s。

2 二維連續(xù)排放污染源溯源算法設(shè)計

式(1)通常用于二維對流擴散連續(xù)源信息已知的情況下，分析污染源排入河流后導(dǎo)致的污染物濃度變化情況。而在溯源反問題中，污染源信息(排放位置、排放時刻、排放速率)都是未知的。河流某固定監(jiān)測點處的污染物濃度C，x和y方向上的彌散系數(shù)Dx和Dy，x方向上的流速u，y方向上的流速v均可由河流實際監(jiān)測情況獲得。通過這些可知量，利用式(1)的逆運算，可求得源強信息。因此，本文基于上面的思路，對二維對流擴散連續(xù)源積分模型式(1)進行一系列的逆向推導(dǎo)，來反向計算連續(xù)污染源的源強信息，那么式(1)可寫為

(2)

(3)

其中，C′是關(guān)于x,y,t的函數(shù)關(guān)系式，對C′(x,y,t)的x,y,t分別求導(dǎo)，得到

(4)

(5)

(6)

為了方便計算，將式(4)～式(6) 中的已知量設(shè)為

(7)

(8)

(9)

將式(7)～式(9)代入式(4)～式(6)中，求得

x=(u-a)×(t-τ)

(10)

y=-b(t-τ)

(11)

(12)

3 兩點和四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)方式

3.1 兩點應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)

圖1為兩點應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)示意圖，圖中以污染源排放點為原點(0，0)，河流流向平行方向為x軸，河流流向垂直方向為y軸。在x方向上距離固定監(jiān)測點P1的Δx處布置一個應(yīng)急監(jiān)測點P2，在y軸方向上距離固定監(jiān)測點P1的Δy處布置一個應(yīng)急監(jiān)測點P3。

圖1 兩點應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)示意

此時

(13)

式(13)中：CP1(T)為P1在T時刻的監(jiān)測濃度；CP1(T+τ)為P1在T+τ時刻的監(jiān)測濃度；CP2(T)為P2在T時刻的監(jiān)測濃度；CP2(T+τ)為P2在T+τ時刻的監(jiān)測濃度。

(14)

式(14)中：CP3(T)為P3在T時刻的監(jiān)測濃度；CP3(T+τ)為P3在T+τ時刻的監(jiān)測濃度。

(15)

式(15)中：CP1(T+Δt)為P1在T+Δt時刻的監(jiān)測濃度；CP1(T+τ+Δt)為P1在T+τ+Δt時刻的監(jiān)測濃度。

連續(xù)排放污染源溯源算法工作步驟如下：

1)假設(shè)(0,0)處為污染源排放位置，河流下游固定監(jiān)測點(X,Y)處監(jiān)測到污染物濃度超標。

2)及時在x方向上與固定監(jiān)測點P1相距Δx處布置應(yīng)急監(jiān)測點P2，在y方向上與固定監(jiān)測點P1相距Δy處布置應(yīng)急點P3。

3)在T時刻，P1，P2，P3同時測量河流水質(zhì)，得到這三個點的污染物濃度，分別為CP1(X,Y,T),CP2(X+Δx,Y,T)，CP3(X,Y+Δy,T)。

4)在T+τ時刻，P1，P2，P3同時測量河流水質(zhì)，得到這三個點的污染物濃度，分別為CP1(X,Y,T+τ)、CP2(X+Δx,Y,T+τ)、CP3(X,Y+Δy,T+τ)。

5)在T+Δt時刻，P1再次監(jiān)測河流污染物濃度，為CP1(X,Y,T+Δt)。

6)在T+Δt+τ時刻，P1再次測量河流污染物濃度，為CP1(X,Y,T+Δt+τ)。

7)將上述步驟中測得的污染物濃度以及Δx，Δy，Δt，τ代入式(4)～式(15)，解得連續(xù)源信息的計算值x，y，t。將計算出來的x，y，t代入式(1)，可解得污染源排放速率C0Q。

3.2 四點應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)

在x和y方向上設(shè)置對稱應(yīng)急監(jiān)測點，即四點應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)方法，示意圖如圖2所示。

圖2 四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)示意

此時

(16)

(17)

(18)

式(18)中:CP2(T),CP4(T)分別為同一時間x方向上與P1相距Δx處對稱的兩個應(yīng)急監(jiān)測點在T時刻的污染物監(jiān)測值;CP3(T)，CP5(T)分別為同一時間y方向上與P1相距Δy處對稱的兩個應(yīng)急監(jiān)測點在T時刻的污染物監(jiān)測值。CP1(T)和CP1(T+Δt)分別為P1在T時刻以及T+Δt時刻的污染物監(jiān)測值。

溯源工作步驟如下：

1)假設(shè)(0,0)處發(fā)生突發(fā)水污染事故，河流下游固定監(jiān)測點(X,Y)處監(jiān)測到污染物濃度超標。

2)及時在x方向上與P1相距±Δx處布設(shè)應(yīng)急監(jiān)測點P2,P4，在y方向上與P1相距±Δy處布設(shè)應(yīng)急監(jiān)測點P3,P5。

3) 在T時刻，P1,P2,P3,P4,P5同時測量河流污染物濃度，得到它們在T時刻的污染物監(jiān)測濃度值，分別為CP1(X,Y,T),CP2(X+Δx,Y,T),CP4(X-Δx,Y,T),CP3(X,Y+Δy,T),CP5(X,Y-Δy,T)。

4)在T+Δt時刻，P1再次監(jiān)測污染物濃度，為CP1(X,Y,T+Δt)。

5)將CP1(X,Y,T),CP2(X+Δx,Y,T),CP4(X-Δx,Y,T),CP3(X,Y+Δy,T),CP5(X,Y-Δy,T)，CP1(X,Y,T+Δt)和Δx、Δy、Δt代入式(16)、式(18)以及式(7)、式(12)，獲得污染源信息(X,Y,T)的計算值(x,y,t)，將計算出來x，y，t的值代到式(1)中，可解得污染源排放量。

4 研究實踐與分析

筆者提出的河流水污染事故快速溯源方法是基于恒定流下污染物二維對流擴散方程的解析解，是在水文理想條件下的純理論概化推導(dǎo)，目前研究階段僅適用于人工輸水河渠背景下的污染物溯源。筆者取溫州地區(qū)瑞安市環(huán)城河馬嶼鎮(zhèn)段作為研究對象，主要原因為：1) 瑞安市環(huán)城河是人工輸水河渠，且河流灣道少，相對較直，符合研究需要的基本河道特性。2) 基于GIS的浙江省工業(yè)企業(yè)空間分布特征與水污染關(guān)系研究[17]，馬嶼鎮(zhèn)工業(yè)發(fā)達，河流水污染事故發(fā)生概率較高，具有比較高的研究價值。

以瑞安市環(huán)城河馬嶼鎮(zhèn)段河道為研究對象，討論了兩點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布局參數(shù)Δx、Δy和T對溯源算法精度的影響，以及對比分析兩點和四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)形式。假設(shè)污染源以C0Q=0.01 kg/s的速率排入河渠的上游，參見圖1，P1在x和y方向上與污染源位置相距X=5 400 m，Y=0 m，擴散系數(shù)Dx=50 m2/s、Dy=0.1 m2/s；河流流速u=0.5 m/s，v=0 m/s。P1在監(jiān)測到污染物濃度超標時報警，系統(tǒng)需要通過應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)和污染源溯源方法，進行溯源任務(wù)，得到連續(xù)源的排放速率Q，排放位置X和Y以及排放時間T的計算值。在溯源算法中，水深h，彌散系數(shù)Dx和Dy，流速u和v均為恒定值，而應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)法中的參數(shù)T，Δx，Δy，Δt的值都是可變因子。

4.1 溯源計算與分析

從數(shù)學(xué)角度看，當(dāng)Δx→0、Δy→0、Δt→0時，溯源結(jié)果最為準確。但由于現(xiàn)實條件限制，Δx,Δy,Δt不可能無限趨向于0。因此，設(shè)置了Δx,Δy的不同因子水平值。在不同因子水平下，通過比較溯源計算值和設(shè)計值之間的誤差，得到優(yōu)化的應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)參數(shù)。

4.1.1Δx對溯源方法精度的影響分析

Δx取6個水平值，分別為1，5，10，30，100,300 m。其余參數(shù)取值如下：Dx=50 m2/s，Dy=0.1 m2/s，u=0.5 m/s，v=0 m/s，T=3.0 h，Δy=1 m，Δt=30 s，τ=30 s。通過溯源算法得到污染源信息的計算值，并計算溯源誤差，結(jié)果如圖3所示。

由圖3可以看出，當(dāng)Δx設(shè)置為1～300 m時，X的誤差分別為3.25%，4.25%，5.29%，9.25%，20.89%,42.13%;T的誤差分別為0.08%，1.16%，2.27%，6.51%，18.99%,41.74%;C0Q的誤差分別為2.90%，4.00%，5.13%，9.41%，21.88%,44.14%。分析上述數(shù)據(jù)可得，Δx越小溯源誤差越小。Δx在30 m以內(nèi)的溯源誤差小于10%，在可接受范圍內(nèi)。而隨著Δx的增大，誤差越來越大，會嚴重影響溯源精度。

圖3 應(yīng)急監(jiān)測網(wǎng)絡(luò)參數(shù)Δx對溯源算法精度的影響

4.1.2Δy對溯源算法精度的影響分析

Δy取5個水平值，分別為1，5，10，30,100 m。其余參數(shù)取值如下：Dx=50 m2/s，Dy=0.1 m2/s，u=0.5 m/s，v=0 m/s，T=3.0 h，Δx=1 m，Δt=30 s，τ=30 s.通過溯源算法得到污染源信息的計算值，并計算溯源誤差，對溯源算法精度的影響計算結(jié)果如圖4所示。

圖4 應(yīng)急監(jiān)測網(wǎng)絡(luò)參數(shù)Δy對溯源算法精度的影響

由圖4可以看出，當(dāng)Δy設(shè)置為1～100 m時，X的誤差分別為3.25%，3.37%，3.75%，6.97%,10.91%;T的誤差分別為0.08%，0.21%，0.60%，3.92%，7.97%;C0Q的誤差分別為2.90%，2.91%，2.91%，3.02%,3.31%。分析上述數(shù)據(jù)可得，Δy越小溯源誤差越小。當(dāng)Δy為100 m時，其X值的相對誤差為10.91%，在可接受范圍內(nèi)。因此，對于應(yīng)急監(jiān)測網(wǎng)絡(luò)中的Δy可以不做太高要求。

4.1.3 T對溯源算法精度的影響

T分別取1.0 h(代表污染團遠未到達P1點)、2.0 h(污染團前端到達P1點)、3.0 h(代表污染團中心到達P1點)、4.0 h(代表污染團尾翼到達P1點)、6.0 h(代表污染團遠離P1點)。其余參數(shù)取值如下：Dx=50 m2/s，Dy=0.1 m2/s，u=0.5 m/s，v=0 m/s，Δx=1 m，Δy=1 m，Δt=30 s,τ=30 s。通過溯源算法得到連續(xù)源信息的計算值，并計算溯源誤差，結(jié)果如圖5所示。

圖5 T對溯源算法精度的影響

由圖5看出，在T=1.0 h時，污染還未到達監(jiān)測點，此時溯源誤差非常大;在T=2.0 h時，監(jiān)測點受污染的影響較小，此時溯源誤差仍比較大;在T=3.0 h時，溯源誤差是最小的，此時污染團中心接近固定監(jiān)測點，對監(jiān)測點的影響最大，溯源結(jié)果最為可靠。根據(jù)二維對流擴散連續(xù)方程，污染物濃度在后期將保持不變。因此在T=4.0 h之后，應(yīng)急監(jiān)測網(wǎng)絡(luò)的濃度差將縮小，導(dǎo)致溯源誤差增加。

4.2 應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)形式對算法精度的影響分析

以應(yīng)急監(jiān)測網(wǎng)絡(luò)參數(shù)Δx為例，分析四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布置法對溯源誤差的影響，并將結(jié)果與二點應(yīng)急監(jiān)測網(wǎng)絡(luò)布置法進行對比。Δx取值分別為1，5，10，30，100，300 m。其余參數(shù)取值如下：Dx=50 m2/s，Dy=0.1 m2/s，u=0.5 m/s，v=0 m/s，T=3.0 h，Δy=1 m，Δt=30 s，τ=30 s。通過溯源算法得到污染源信息的計算值，并計算溯源誤差，結(jié)果如圖6所示。

圖6 應(yīng)急監(jiān)測網(wǎng)絡(luò)布置方式對溯源精度的影響

由圖6可以看出，四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)法相比二點法，溯源誤差都要小，說明設(shè)置對稱應(yīng)急監(jiān)測點可以降低溯源誤差。且在同一T下，溯源誤差隨著Δx的增大變化不大，說明設(shè)置對稱應(yīng)急監(jiān)測點可以提高溯源誤差的穩(wěn)定性。因此，在條件允許的情況下，采用四點對稱應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)法明顯優(yōu)于兩點應(yīng)急監(jiān)測網(wǎng)絡(luò)布設(shè)，從而降低溯源誤差。

筆者提出的河流水污染事故快速溯源方法是在恒定流污染物二維擴散方程解析解(式(1)的基礎(chǔ)上逆向推導(dǎo)得到，是水文理想條件下的求解方法。因此，目前溯源方法的研究成果僅適用于滿足恒定流條件的人工順直輸水渠道，不適于水深、流量和流速等水力要素隨時間和地點變化的自然河道。為促進該方法由理論分析向應(yīng)用轉(zhuǎn)化，應(yīng)進一步開展如下研究工作：應(yīng)結(jié)合室內(nèi)瑞安市環(huán)城河馬嶼鎮(zhèn)段水槽模型進一步開展算法的驗證性工作；加強污染物測量誤差對溯源精度的影響分析。目前方程中僅考慮了污染物隨水流的遷移和擴散作用，未考慮污染物的沉降、再懸浮、吸附和降解等作用，應(yīng)加強該部分的研究工作。同時，將研究成果應(yīng)用到瑞安市環(huán)城河馬嶼鎮(zhèn)段的實際監(jiān)測中，并不斷矯正參數(shù)，達到可用可塑的目的。

5 結(jié) 論

1)提出了河流二維連續(xù)排放污染源的溯源算法，實現(xiàn)了對河流二維對流擴散方程的反向求解，解決了污染源位置、釋放強度和釋放時間等源強信息計算的理論公式問題。筆者提出的兩點對稱和四點對稱的河流應(yīng)急監(jiān)測網(wǎng)絡(luò)的布設(shè)方法，解決了目標污染物濃度在空間域和時間域的差分計算問題，從而為河流二維連續(xù)污染源的溯源算法的應(yīng)用奠定了基礎(chǔ)。

2)應(yīng)急監(jiān)測網(wǎng)絡(luò)影響因子Δx，Δy，T等對溯源誤差的影響趨勢相同。相對而言，Δx和T對溯源精度的影響更大，Δy對溯源精度的影響較小。固定監(jiān)測點與污染團中心的相對位置對溯源精度影響很大。在污染團中心接近或到達常設(shè)監(jiān)測點時，運用該溯源方法的相對誤差最小。應(yīng)急監(jiān)測網(wǎng)絡(luò)中應(yīng)急監(jiān)測點的位置離常設(shè)監(jiān)測點越近，相對誤差越小。四點對稱應(yīng)急網(wǎng)絡(luò)布置法與二點應(yīng)急網(wǎng)絡(luò)布置法相比，相對誤差更小，且當(dāng)布置條件發(fā)生改變時，四點對稱應(yīng)急網(wǎng)絡(luò)布置法的溯源誤差穩(wěn)定性更高。

3)通過將瑞安市環(huán)城河馬嶼鎮(zhèn)段作為實際研究對象，設(shè)計的污染場景以及不同應(yīng)急監(jiān)測網(wǎng)絡(luò)布局形式下污染源設(shè)計值與計算值的相對誤差計算，驗證了溯源算法的可行性。該溯源方法可提高對水污染源事件的應(yīng)急響應(yīng)能力，有效并快速地為污染處理提供科學(xué)的理論依據(jù)。

4)在實踐應(yīng)用中，應(yīng)開展算法的驗證性工作，分析水環(huán)境條件、模型參數(shù)和測量誤差等要素對溯源的影響，進一步完善和論證溯源方法的可行性和可靠性。