基于靈敏度分析和替代模型的地下水污染風險評價方法

常振波,盧文喜*,辛 欣,顧文龍,崔尚進(1.吉林大學,地下水與資源環境教育部重點實驗室,吉林長春 130012；2.吉林大學環境與資源學院,吉林 長春 130012)

基于靈敏度分析和替代模型的地下水污染風險評價方法

常振波1,2,盧文喜1,2*,辛 欣1,2,顧文龍1,2,崔尚進1,2(1.吉林大學,地下水與資源環境教育部重點實驗室,吉林長春 130012；2.吉林大學環境與資源學院,吉林 長春 130012)

采用蒙特卡洛方法,借助隨機模型進行地下水污染風險評價,模型中隨機變量利用靈敏度分析的方法確定,使地下水風險評價結果更為可靠,并借助一個假想例子來說明評價過程.結果表明,模擬輸出數據符合正態分布規律,對正態分布概率密度函數積分可以得到污染風險,井1、井2和井3的污染風險分別為0%、78.52%和100%;根據整個模擬區的污染風險分布圖可以劃分出具有不同污染風險程度的子區域,藉此能夠定量評價模擬區不同子區域的污染風險程度.

地下水污染風險評價；靈敏度分析；蒙特卡洛；替代模型；概率統計

地下水人為污染事件時有發生[1],相比地表水污染,地下水污染修復所需費用更高.地下水污染具有發生的隱蔽性和發現的滯后性[2-3],污染一經發現難以修復.為了緩解地下水質惡化的現狀,可靠的地下水污染風險評價工作是十分必要的.決策者可以根據風險評價結果,通過人為規劃和調控降低污染發生的可能性[4-5].

風險評價實質就是不確定性分析,沒有不確定性,就沒有風險[6].因此,對地下水溶質運移中不確定性因素進行定量分析,并在結果中體現風險程度,可以使地下水污染風險評價結果更加科學.對此國內外展開了大量的研究,Goodrich等[7]應用蒙特卡洛模擬,研究由于參數不確定性導致地下水流和溶質運移過程的不確定性;Bennett等[8]將蒙特卡洛方法得到的一系列污染物遷移結果用于暴露評價;梁婕[9]利用貝葉斯方法推斷出滲透系數隨機場的相關參數的后驗分布,基于參數的后驗分布研究地下水溶質運移規律;申升[10]進一步考慮了縱向彌散度對地下水溶質運移的影響.

前人研究中隨機變量多為事先確定,然而對于不同的模型參數對結果的影響是不同的,引入靈敏度分析方法,篩選對結果影響較大的參數作為隨機變量,利用蒙特卡洛方法對地下水溶質運移進行不確定性分析,并以濃度值超過某一限定值的概率作為污染風險,對整個模擬區的污染風險進行定量評價.為了減少大量重復調用模擬模型而產生的計算負荷,本文借助替代模型完成蒙特卡洛過程.

1 研究方法

1.1 蒙特卡洛法

蒙特卡洛法作為最常用的不確定性分析方法之一被廣泛地應用于地下水數值模擬中[11-13].蒙特卡洛方法就是通過設置一系列平行的確定性實驗來模擬隨機問題.在設置平行實驗時選擇待求因素作為隨機變量,研究隨機變量變化對結果的影響.蒙特卡洛法將參數的不確定性轉化為溶質運移的不確定性,能夠很好地表征污染風險[14-15].本文以蒙特卡洛方法作為框架,具體研究思路如圖1所示.

1.2 靈敏度分析

靈敏度分析可以定量地評價參數不確定性對結果造成的影響,本文利用靈敏度分析的方法分析參數對輸出濃度的影響,挑選兩個靈敏度系數較大的參數作為隨機變量,進行地下水污染風險評價,這樣會使得風險評價更為可靠.分析方法可分為全局靈敏度分析和局部靈敏度分析,其中局部靈敏度分析能得到單一參數變化對結果的影響.靈敏度利用結果對參數求偏導數,反映參數變化對結果的影響程度[16-17],即:

式中:Sk表示當參數αk變化時對輸出結果y的影響程度,也就是靈敏度系數.

對于某一特定參數靈敏度系數的求解,可通過以下方法近似,使該參數的值由αk變化為αk+Δαk,結果則由 yi( αk)變化為 yi( αk+Δαk),可用以下公式求得近似值,即:

1.3 替代模型

少量的數據無法支撐地下水污染風險評價,而大量數據的獲取需要大量的時間,其中時間主要花費在模擬模型的求解上.為了減少模擬模型的使用次數,引入替代模型.替代模型能夠在某些行為上逼近模擬模型,且調用簡便,運行時間短.蒙特卡羅模擬與替代模型結合能有效減小計算負荷[18-19].

克里格方法利用協方差的變化來表達空間的變化,現被延伸為一種替代模型的建立方法,廣泛地應用于多個行業,是一種黑箱模型[20].

克里格模型的方程可以表示成以下形式:

上述對 y( x)的估計值y?( x)可以分為兩部分,其中前一部分為線性回歸部分,后一部分為隨機部分.其中 f( x) =[f1( x) , f2(x) ,… , fk(x )]為已知回歸模型的基函數;待定參數 β = [β1, β2,… ,βk],可以通過訓練數據求得;z( x)為隨機部分,其方差為σ2,均值為0,協方差:

R( xi, xj)為采樣點xi和點xj的關聯函數,有多種類型可供選擇,本文采用較為常用的高斯模型:

式中:θk為待定參數,通過求解優化模型求得;xki為第i個樣本的k維坐標.

根據克里格模型,在預測點 x處的響應值y( x)的預測估計值為:

式中:f為基函數,為方便起見,本文選擇常數型,f為一常數列向量;r( x)為點x與n個訓練樣本采樣點 (x1, x2,… ,xn)之間的相關向量,;y為與n個采樣點對應的響應值,為 n×1的向量;β為線性回歸部分的待定參數,可以通過最優線性無偏估計求得:

R為n個采樣點的相關矩陣:

方差σ2的估計值為:

待定參數 θk可以通過一個無約束優化問題求得:

通過MATLAB軟件實現上述替代模型的建立過程.

2 假想算例

2.1 概況

圖1 技術路線流程Fig.1 The flow chart of technical route

模擬區內存在二維均質各向同性潛水含水層,含水層信息如圖1所示,含水層厚為12m,模擬區初始水位為8m,水力梯度約為0.007,年平均降水量730mm,降雨入滲系數為0.2,選擇滲透系數、給水度、縱向彌散度和孔隙度4個參數進行靈敏度分析,根據經驗給出其取值范圍和概率分布,如表1所示,其中因橫向彌散度與縱向彌散度的比值為 0.1,故只考慮縱向彌散度.區內共設三口抽水井,定流量抽水,各井抽水量均為300m3/d,持續時間 1a,同時將抽水井作為水質觀測井.現假設在模擬區內建造一座工廠,工廠每天向含水層排放400m3的污水,污染物濃度500mg/L,污染質為保守物質,遷移過程中不發生物理化學反應.含水層污染物的初始濃度為零,定水頭邊界可以視為零濃度邊界,隔水邊界可以視為零通量邊界.

表1 參數的概率分布及取值情況Table 1 The probability distribution and values information of parameters

2.2 隨機模型建立

2.2.1 數學模型 地下水溶質運移模型是在地下水水流模型的基礎上建立的,研究區地下水為二維均質各向同性潛水非穩定流,可以建立數學模型:

地下水溶質運移數學模型:

地下水水流模型和溶質運移模型通過達西定律聯系起來:

式中:K為含水層的滲透系數,m/d;H為潛水水位, m;B為潛水含水層底板高程,m;P 為人工開采強度,m/d;R為降水入滲補給量,m/d;μ為含水層給水度,無量綱;G為模擬區范圍;S1、 S2為含水層邊界,在水流模型中為定水頭邊界,在溶質運移模型中為定濃度邊界;S3、 S4為含水層邊界,在水流模型中為零流量邊界,在溶質運移模型中為零彌散通量邊界;Kn為邊界法向量上的滲透系數,m;n為含水層介質的孔隙度,無量綱;c為污染質濃度,mg/L;Dx、Dy為水動力彌散系數在x、y方向的分量,m2/d;vx、 vy為滲透流速v在x、y方向上的分量,m/d.

圖2 研究算例平面示意Fig.2 Plan view of the aquifer configuration

數學模型是借助GMS軟件的MODFLOW模塊和MT3DMS模塊進行求解的.

2.2.2 隨機變量確定 本文利用靈敏度分析的方法確定隨機變量,參與分析的參數有滲透系數K、給水度μ、縱向彌散度 αL和孔隙度n.由圖3可知,4個參數的靈敏度系數由大到小依次為SK、 Sα、Sn、 Sμ,滲透系數和縱向彌散度的靈敏度系數較大,選擇這兩個參數作為隨機變量,按照其經驗概率分布進行取值,剩下的兩個參數作為確定值輸入模型.

2.3 替代模型的建立

利用拉丁超立方抽樣的方法對隨機變量滲透系數和縱向彌散度按照其概率分布進行抽樣,得到50組輸入數據集,將50組參數組合分別輸入模擬模型,利用GMS求解,得到50組輸出結果.將參數組合和輸出結果中三口觀測井的末時刻濃度值分別作為輸入輸出,利用40組輸入和輸出建立克里格替代模型,并利用 10組數據(每一組數據均包括3口觀測井的濃度值)驗證替代模型的可靠性.圖4為模擬模型輸出結果和替代模型輸出結果對比圖,可以看出二者的相對誤差在2%以下,可認為替代模型滿足誤差要求,并能夠替代模擬模型進行調用.

圖3 參數靈敏度曲線Fig.3 Sensitivity curves of parameters

2.4 蒙特卡羅模擬

利用蒙特卡羅方法模擬500次試驗,對結果進行統計.即首先根據滲透系數和縱向彌散度的概率分布,按照拉丁超立方抽樣的方法,得到500組輸入數據集;然后將500組參數組合輸入替代模型,得到 500組輸出數據;對輸出數據進行分析.

3 結果分析

3.1 單井污染風險分析

利用SPSS軟件中的單樣本K-S檢驗分析500組輸出,假設三口井的濃度值為正態分布,可以得到結果如表 2所示,其中當 P值大于 0.05時認為假設成立,得到三口井的分布均為正態分布.由此畫出三口井濃度值的累積頻率直方圖和其對應的正態分布概率密度函數,如圖5(a)～(c)所示.

圖4 模型預測結果對比Fig.4 Comparison of different models’ result

表2 輸出數據單樣本K-S分析結果Table 2 The results of K-S single sample analysis on output data

圖5 觀測井污染物濃度值直方圖及正態分布概率密度函數Fig.5 Observation wells pollutant concentration histogram and normal probability density function

分析三口供水井的污染風險,以污染物濃度值超過100mg/L的地下水為是被污染的,將超過100mg/L概率作為污染風險大小,利用如下公式可計算三口井污染風險:

式中:P為單井污染風險大小;()f x為單井濃度值的概率密度函數;σ為正態分布的標準差;μ為概率分布的均值.井1、井2和井3的污染風險分別為0%、78.52%和100%.由此可以得到,井2和井3水質均有較大的污染風險,為此可能要對工廠排放污水的濃度做出一定的限制或者尋找新的供水井,下面通過對區域的污染風險進行評價來為尋找新的供水井提出一些建議.

3.2 區域污染風險分析

圖6 模擬區污染風險分布Fig.6 the pollution risk distribution map of the whole simulation area

為更加全面地分析整個研究區的污染風險,計算出每個離散網格所對應位置的污染風險,并繪制等值線圖見圖 6.可以看出污染風險在不同空間位置上具有差異性,研究區上方區域污染風險較小且范圍較小,下方區域污染風險較大且范圍較大.井1位于白色區域即污染低風險區域,在此區域的地下水受到污染的風險較小,低于20%;井2位于淺色區域即污染中等風險區域,此區域地下水受到污染的風險在20%到80%;井3位于深色區域即污染高風險區域,此區域地下水受到污染的風險在80%以上.井2和井3均有較高的污染風險,為了使地下水水質能夠維持使用標準,需要尋找新的供水井.圖中白色區域風險小,且靠近定水頭邊界水量充足,可以在此區域內選擇新供水井.

4 結論

4.1 在文中假設情境下靈敏度分析方法篩選出滲透系數和縱向彌散度作為隨機參數,利用這兩個參數作為隨機變量得到的地下水污染風險評價結果更加科學和可靠.

4.2 利用本文的方法進行地下水風險評價,能夠得到整個模擬區污染風險分布圖,根據此圖可以確定污染風險小的新供水井位置.

4.3 本文利用蒙特卡羅的方法與克里格替代模型,大大減少了計算負荷,而且克里格替代模型預測結果與模擬模型模擬結果相比,相對誤差均小于 2%,說明替代模型可以較好地模仿模擬模型的行為.

[1] 董文福,傅德黔.近年來我國環境污染事故綜述 [J]. 環境科學與技術, 2009,32(7):75-77.

[2] 王 昊.我國地下水污染成因及防治措施研究 [J]. 資源節約與保護, 2013(5):88-97.

[3] 王焰新.地下水污染與防治 [M]. 北京:高等教育出版社, 2007.

[4] 張麗君.地下水脆弱性和風險性評價研究進展綜述 [J]. 水文地質工程地質, 2006,33(6):113-119.

[5] 劉增超,董 軍,何連生,等.基于過程模擬的地下水污染風險評價方法研究 [J]. 中國環境科學, 2013,33(6):1120-1126.

[6] 馬祿義,許學工,徐麗芬.中國綜合生態風險評價的不確定性分析 [J]. 北京大學學報:自然科學版, 2011(5):893-900.

[7] Goodrich M T, Mccord J T. Quantification of uncertainty in exposure assessments at hazardous waste sites [J]. Ground Water, 1995,33(5):727-732.

[8] Bennett D H, James A L, Mckone T E, et al. On uncertainty in remediation analysis: variance propagation from subsurfacetransport to exposure modeling [J]. Reliability Engineering and System Safety, 1998,62(1/2):117-129.

[9] 梁 婕.基于不確定理論的地下水溶質運移及污染風險研究[D]. 長沙:湖南大學, 2009.

[10] 申 升.基于貝葉斯理論的地下水溶質運移的不確定性研究[D]. 長沙:湖南大學, 2012.

[11] 王 宇,盧文喜,卞建民,等.基于小波神經網絡的地下水流數值模擬模型的替代模型研究 [J]. 中國環境科學, 2015,35(1):139-146.

[12] 歐陽琦,盧文喜,侯澤宇,等.基于替代模型的地下水溶質運移不確定性分析 [J]. 中國環境科學, 2016,36(4):1119-1124.

[13] 曾獻奎,王 棟,吳吉春.地下水流概念模型的不確定性分析 [J].南京大學學報(自然科學), 2012,48(6):746-752.

[14] 姚磊華.白楊水源地潛水水質的 Monte—Carlo隨機模擬 [J].長春科技大學學報, 1998(3):296-302.

[15] Hamed M M, Conte J P, Bedient P B. Probabilistic screening tool for ground-water contamination assessment [J]. Journal of Environmental Engineering, 1995,121(11):767-775.

[16] 束龍倉,劉佩貴,劉 波,等.傍河水源地數學模型的參數靈敏度分析—以遼寧省北票市某傍河水源地為例 [J]. 工程勘察, 2006,(8):29-31.

[17] 束龍倉,王茂枚,劉瑞國,等.地下水數值模擬中的參數靈敏度分析 [J]. 河海大學學報(自然科學版), 2007,35(5):491-495.

[18] 安永凱,盧文喜,董海彪,等.基于克里格法的地下水流數值模擬模型的替代模型研究 [J]. 中國環境科學, 2014,34(4):1073-1079.

[19] 顧文龍,盧文喜,馬洪云,等.地下水數值模擬分析中降水入滲補給強度及滲透系數不確定性評價 [J]. 水電能源科學, 2015(11): 45-48.

[20] 吳義忠.多領域物理系統的仿真優化方法 [M]. 北京:科學出版社, 2011.

Groundwater contamination risk assessment method based on sensitivity analysis and surrogate model.

CHANG Zhen-bo1,2, LU Wen-xi1,2*, XIN Xin1,2, GU Wen-long1,2, CUI Shang-jin1,2
(1.Key Laboratory of Groundwater Resources and Environment Ministry of Education, Jilin University, Changchun 130012, China；2.College of Environment and Resources, Jilin University, Changchun 130012, China). China Environmental Science, 2017,37(1)：167～173

Stochastic model with Monte Carlo method was used to assess the risk of groundwater pollution. The random variables in the model were determined by the sensitivity analysis, which made the result of the groundwater risk assessment more reliable, and the evaluation process was illustrated by a hypothetical example. The results showed that the simulated output data was in accord with the normal distribution law, and normal probability density function of the integrator was the risk of contamination. The pollution risks of well 1, 2 and 3 were 0%, 78.52% and 100%, respectively. Sub regions with different risk of pollution were divided according to the pollution risk distribution map of the whole simulation area, so as to quantitatively evaluate the pollution risk of different sub area in the simulation area.

groundwater contamination risk assessment；sensitivity analysis；Monte Carlo；surrogate model；probability statistics

X703

A

1000-6923(2017)01-0167-07

常振波(1993-),男,河南開封人,吉林大學碩士研究生,主要從事地下水溶質運移的不確定性研究.

2016-05-15

國家自然科學基金項目(41372237);吉林大學研究生創新基金資助項目(2016207,2016100)

* 責任作者, 教授, luwenxi@jlu.edu.cn