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

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

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

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

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

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

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

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

1 研究方法

1.1 蒙特卡洛法

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

1.2 靈敏度分析

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

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

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

1.3 替代模型

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

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

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

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

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

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

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

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

R為n個采樣點的相關(guān)矩陣:

方差σ2的估計值為:

待定參數(shù) θk可以通過一個無約束優(yōu)化問題求得:

通過MATLAB軟件實現(xiàn)上述替代模型的建立過程.

2 假想算例

2.1 概況

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

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

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

2.2 隨機模型建立

2.2.1 數(shù)學模型 地下水溶質(zhì)運移模型是在地下水水流模型的基礎(chǔ)上建立的,研究區(qū)地下水為二維均質(zhì)各向同性潛水非穩(wěn)定流,可以建立數(shù)學模型:

地下水溶質(zhì)運移數(shù)學模型:

地下水水流模型和溶質(zhì)運移模型通過達西定律聯(lián)系起來:

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

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

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

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

2.3 替代模型的建立

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

圖3 參數(shù)靈敏度曲線Fig.3 Sensitivity curves of parameters

2.4 蒙特卡羅模擬

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

3 結(jié)果分析

3.1 單井污染風險分析

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

圖4 模型預(yù)測結(jié)果對比Fig.4 Comparison of different models’ result

表2 輸出數(shù)據(jù)單樣本K-S分析結(jié)果Table 2 The results of K-S single sample analysis on output data

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

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

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

3.2 區(qū)域污染風險分析

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

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

4 結(jié)論

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

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

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

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

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

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

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

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

[6] 馬祿義,許學工,徐麗芬.中國綜合生態(tài)風險評價的不確定性分析 [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] 梁 婕.基于不確定理論的地下水溶質(zhì)運移及污染風險研究[D]. 長沙:湖南大學, 2009.

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

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

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

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

[14] 姚磊華.白楊水源地潛水水質(zhì)的 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] 束龍倉,劉佩貴,劉 波,等.傍河水源地數(shù)學模型的參數(shù)靈敏度分析—以遼寧省北票市某傍河水源地為例 [J]. 工程勘察, 2006,(8):29-31.

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

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

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

[20] 吳義忠.多領(lǐng)域物理系統(tǒng)的仿真優(yōu)化方法 [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-),男,河南開封人,吉林大學碩士研究生,主要從事地下水溶質(zhì)運移的不確定性研究.

2016-05-15

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

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