李久輝,盧文喜*,常振波,李明彧,苗添升,趙 瑩,張將偉(1.吉林大學(xué)地下水與資源環(huán)境教育部重點實驗室,吉林 長春 130012；2.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130012)

基于不確定性分析的地下水污染超標(biāo)風(fēng)險預(yù)警

李久輝1,2,盧文喜1,2*,常振波1,2,李明彧1,2,苗添升1,2,趙 瑩1,2,張將偉1,2(1.吉林大學(xué)地下水與資源環(huán)境教育部重點實驗室,吉林 長春 130012；2.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130012)

地下水污染預(yù)測可以通過地下水污染質(zhì)運移數(shù)值模擬模型予以實現(xiàn),為分析模型參數(shù)取值不確定性對模型輸出結(jié)果的影響,本文運用蒙特卡洛方法對模型輸出結(jié)果進行不確定性分析.為降低數(shù)值模擬模型復(fù)雜程度,運用靈敏度分析方法篩選對模型影響較大參數(shù)作為模型中隨機變量;為減少重復(fù)調(diào)用數(shù)值模擬模型產(chǎn)生的計算負(fù)荷,在保證一定精度前提下,運用克里格方法建立模擬模型的替代模型完成模擬過程.結(jié)果表明：應(yīng)用概率密度函數(shù)積分可以估計地下水遭受污染風(fēng)險與不同置信程度下污染物濃度區(qū)間.污染羽分布圖與分級污染超標(biāo)風(fēng)險預(yù)警圖可以分別對研究區(qū)不同等級污染覆蓋面積和研究區(qū)不同污染風(fēng)險對應(yīng)污染羽分布進行估計.基于污染質(zhì)運移數(shù)值模擬不確定性分析的地下水污染超標(biāo)風(fēng)險預(yù)警可以更加客觀地對地下水污染問題進行預(yù)測.

地下水污染；數(shù)值模擬模型；蒙特卡洛方法；克里格替代模型；不確定性分析

地下水污染預(yù)測可以運用地下水污染質(zhì)運移數(shù)值模擬模型予以實現(xiàn).過去地下水中溶質(zhì)運移模擬多采用確定性模型求解,沒有將水文地質(zhì)參數(shù)獲取的不確定性考慮在內(nèi),其結(jié)果參考價值不高,基于其結(jié)果的地下水污染預(yù)測與評價亦不夠可靠,嚴(yán)重影響污染預(yù)防與污染治理相關(guān)措施實施,使污染得不到及時處理與改善,導(dǎo)致地下水環(huán)境及生態(tài)環(huán)境污染進一步惡化.故研究參數(shù)的不確定性對模型輸出結(jié)果的影響對于地下水污染預(yù)測與風(fēng)險預(yù)警是十分必要的.

近年來不確定性分析研究在國內(nèi)外發(fā)展迅速.水文地質(zhì)參數(shù)不確定性與地下水研究結(jié)合層出不窮,但多用于水資源量評價,Kuczera等[1]將降雨量模型不確定性分析結(jié)果與降雨數(shù)據(jù)評估相結(jié)合,具有實際意義;曾獻奎等[2]將地下水流概念模型不確定性分析應(yīng)用于模型預(yù)測精度研究,使與地下水相關(guān)研究的不確定性分析進一步發(fā)展;束龍倉等[3]將水文地質(zhì)參數(shù)不確定性分析與地下水補給量研究相結(jié)合,來評價地下水補給量離散程度,不確性分析多用于水量評價,在水質(zhì)問題研究領(lǐng)域卻很少見,直到歐陽琦等[4], 苗添升等[5]將不確定性分析結(jié)果與地下水遭受污染風(fēng)險相結(jié)合,水質(zhì)預(yù)測與不確定性分析結(jié)合才有所進展.

根據(jù)前人經(jīng)驗,將不確定性分析與地下水污染預(yù)測相結(jié)合,運用靈敏度分析方法篩選出對數(shù)值模擬模型輸出結(jié)果影響較大的參數(shù),作為模型中隨機變量以降低模型復(fù)雜程度;為減少重復(fù)調(diào)用模擬模型產(chǎn)生的計算負(fù)荷,在保證一定精度前提下,運用克里格方法建立數(shù)值模擬模型的替代模型完成蒙特卡洛模擬,最后對蒙特卡洛結(jié)果進行統(tǒng)計分析.利用概率密度函數(shù)估計不同置信程度的污染物濃度范圍;運用數(shù)理統(tǒng)計方法分析研究區(qū)各個有限差分網(wǎng)格(運用 GMS軟件剖分離散得到的差分網(wǎng)格)污染物濃度值服從分布規(guī)律,實現(xiàn)對研究區(qū)地下水不同程度污染的污染羽分布預(yù)測與分級污染超標(biāo)風(fēng)險預(yù)警,為地下水污染預(yù)測提供更多選擇依據(jù).

1 研究方法

1.1 蒙特卡洛方法

本文以蒙特卡洛方法作為不確定性分析的的總體研究思路.蒙特卡洛方法的主要思想是當(dāng)待求解的問題是某種隨機事件出現(xiàn)的概率,或者是某個隨機變量的期望值時,通過“實驗”統(tǒng)計的方法,用這個隨機事件在實驗中出現(xiàn)的頻率對隨機事件發(fā)生的概率進行估計,或者利用實驗中隨機事件的數(shù)字特征,作為研究問題的解.在研究過程中無論隨機事件結(jié)果滿足何種分布,在模擬次數(shù)足夠多的情況下,都可以得到一個比較精確的概率分布[6].文中將污染質(zhì)運移數(shù)值模型視為“實驗”發(fā)生器,將模型輸出的污染物濃度結(jié)果視為研究區(qū)污染問題這一隨機事件出現(xiàn)的結(jié)果,最后利用輸出的污染物濃度結(jié)果估計地下水污染事件總體.

1.2 模型建立

建立一個東西長 1000m,南北寬 550m,含水層厚度為20m的研究區(qū),含水層為均質(zhì)各向同性潛水非穩(wěn)定流,初始水位為16m,水流方向由西至東;區(qū)域內(nèi)地表平坦,南側(cè) Г2及北側(cè) Г4均為隔水邊界;西側(cè) Г1有一河流,切穿潛水含水層且與其具有較好的水力聯(lián)系,河流水位穩(wěn)定為18m;東側(cè)Г3有一湖泊,水位穩(wěn)定為16m.研究區(qū)在垂向上均勻的接受大氣降水補給,年平均降水量為580mm;在忽略蒸發(fā)騰散消耗量前提下,地下水排泄方式主要是人工開采,井1、2,3抽水量均為500m3/d.區(qū)內(nèi)一化工廠污水未經(jīng)處理泄漏到地下含水層中,污水污染物濃度為 400mg/L,泄漏量為600m3/d.將污染物視為不會發(fā)生化學(xué)變化及生物轉(zhuǎn)化與遷移的穩(wěn)定污染物.潛水含水層中污染物的初始濃度為零,定水頭邊界Г1,Г3與隔水邊界Г2,Г4分別視為零濃度邊界和零通量邊界.根據(jù)研究區(qū)概況,為更好地對研究區(qū)污染質(zhì)運移情況進行研究,利用GMS將研究區(qū)剖分離散為220個有限差分網(wǎng)格,區(qū)內(nèi)剖分情況及觀測井與污染源位置詳見圖1：

圖1 地下水溶質(zhì)運移模型剖分及污染源,觀測井位置Fig.1 Model mesh of the groundwater solute transport and the location of the pollution source and observation wells

根據(jù)研究區(qū)水文地質(zhì)條件建立研究區(qū)的地下水水流數(shù)值模型：

式中：K為含水層滲透系數(shù),m/d;H為潛水水位,m; B為潛水含水層底板高程,m;Q為人工開采強度,m/d;R為降水入滲補給量,m/d;μ為含水層給水度,無量綱;S為模擬區(qū)范圍;Г1,Г3為定水頭邊界;Г2,Г4為零流量邊界;Kn為邊界法向量上的滲透系數(shù),m.

污染質(zhì)隨水流運移,在建好的地下水水流數(shù)值模型基礎(chǔ)之上,建立地下水污染質(zhì)運移數(shù)值模型：

式中：S為模擬區(qū)范圍;Г1,Г3為定濃度邊界;Г2,Г4為零彌散通量邊界;n為含水層介質(zhì)的孔隙度,無量綱;c為污染質(zhì)濃度,mg/L;Dx,Dy為水動力彌散系數(shù)在x、y方向的分量,m2/d;xyv v、 為滲透流速v在x、y方向上的分量, m/d; I為單位時間單位液相體積內(nèi)溶質(zhì)質(zhì)量的增減量,mg/(d?m3).

在建立好研究區(qū)的污染質(zhì)運移數(shù)值模型基礎(chǔ)之上,利用GMS軟件將研究區(qū)剖分為220個有限差分網(wǎng)格,然后利用其MODFLOW工具箱和MT3DMS工具箱模擬15年后研究區(qū)污染質(zhì)運移情況.

為了簡化模型計算,可以將對模型影響較小的隨機參數(shù)作為定值輸入模型,影響較大的參數(shù)作為隨機變量輸入模型,為實現(xiàn)這一研究過程,應(yīng)該運用有效方法對模型中多個隨機參數(shù)進行逐一篩選.

1.3 靈敏度分析方法篩選隨機參數(shù)

運用靈敏度分析方法篩選出對污染質(zhì)運移數(shù)值模型輸出結(jié)果影響較大的參數(shù),可以在提高研究精度同時減少計算負(fù)荷.靈敏度分析常用的兩種分析方法是全局靈敏度分析方法和局部靈敏度分析方法[7].靈敏度分析應(yīng)用廣泛,束龍倉,王茂枚等[8]對靈敏度分析方法進行了闡述,并將局部靈敏度分析和全局靈敏度分析引入地下水?dāng)?shù)值模型,來模擬塔里木河下游地下水水位變化情況.文中采用局部靈敏度分析方法篩選隨機參數(shù).局部靈敏度分析即利用參數(shù)求偏導(dǎo)數(shù)反映參數(shù)變化對結(jié)果影響程度的大小,來選擇影響較大參數(shù),公式如下：

式中：Sk是靈敏度系數(shù);αk是輸入模型參數(shù)值;△αk是輸入模型參數(shù)變化量;yi(αk+△αk)是參數(shù)變化△αk時,模型輸出結(jié)果;yi(αk)是參數(shù)為 αk時模型輸出結(jié)果.靈敏度系數(shù)越大,說明此參數(shù)對模型輸出結(jié)果影響越大.時,先將模型中參數(shù)取波動均值見表 1,輸出污染物濃度,然后分別將模型中某個隨機參數(shù)加減20%,10%,其他參數(shù)不變,利用式(3)篩選出2個靈敏度系數(shù)較大的參數(shù)作為隨機變量.

圖2 參數(shù)靈敏度分析圖Fig.2 Sensitivity analysis of parameters

根據(jù)圖 2可知,對模型輸出結(jié)果影響較大的2個隨機參數(shù)分別是滲透系數(shù)和縱向彌散系數(shù).將滲透系數(shù)和縱向彌散系數(shù)作為隨機變量,其它參數(shù)取波動均值(定值)輸入模型.

為實現(xiàn)滲透系數(shù)和縱向彌散系數(shù)兩個隨機變量在參數(shù)變化范圍內(nèi)多組不同取值,應(yīng)該選擇適當(dāng)?shù)姆椒▽蓚€隨機參數(shù)分別抽樣取值.

1.4 拉丁超立方方法

運用拉丁超立法抽樣方法對篩選出的兩個隨機參數(shù)分別抽樣取值,可以有效的降低參數(shù)的抽樣誤差,使參數(shù)取值具有代表性.

拉丁超立法抽樣方法是將服從某個分布的參數(shù)在其波動范圍內(nèi)進行等概率分層,并在每一層中抽取相同數(shù)量的參數(shù)組成參數(shù)樣本的一種抽樣方法.根據(jù)前人經(jīng)驗總結(jié),滲透系數(shù),縱向彌散系數(shù)服從分布如表1所示[9].運用MATLAB利用拉丁超立方方法對滲透系數(shù),縱向彌散系數(shù)分別抽樣取值.

表1 參數(shù)概率分布及取值情況Table 1 Values and probability distribution of model parameters

對滲透系數(shù),縱向彌散系數(shù)分別抽樣 50組,運用污染質(zhì)運移數(shù)值模擬模型求解,輸出井1、2、3的50組污染質(zhì)濃度結(jié)果.

因為少量的數(shù)據(jù)沒有辦法支撐地下水污染預(yù)測研究,多次重復(fù)調(diào)用模擬模型又會產(chǎn)生大量的計算負(fù)荷,且耗時較長,故利用模擬模型的 50組隨機參數(shù)與輸出結(jié)果數(shù)據(jù)集建立克里格替代模型代替模擬模型求解.

1.5 建立替代模型

在保證一定精度前提下,為減少運用污染質(zhì)運移數(shù)值模擬模型的計算負(fù)荷建立數(shù)值模擬模型的替代模型完成蒙特卡洛過程.此次研究選擇克里格內(nèi)插法建立替代模型.克里格方法由南非地質(zhì)學(xué)家 Krige[10]提出的一種插值方法,它從變量的相關(guān)性和變異性出發(fā),利用方差的變化表達空間的變化,在引入似然函數(shù)估計誤差最小的情況下,通過已知觀察點估計未知觀察點.目前克里格方法被延伸為一種建立替代模型的方法,替代模型做為一種黑箱模型[11],應(yīng)用于多個工程領(lǐng)域.替代模型建立可以按以下步驟實現(xiàn)：

式中：R(xi,xj)為任意兩采樣點 xi,xj之間的空間相關(guān)關(guān)系方程,也稱關(guān)聯(lián)函數(shù),文中采用 EXP型的關(guān)聯(lián)函數(shù)：

式中：θk為待定參數(shù);xi為第i個樣本的k維坐標(biāo).k

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

式中：r( x)為點x與n個采樣點 (x1,x2,…, xn)的相關(guān)向r( x) =[R( x, x1),R( x, x2),… ,R( x, xn)];y為n個采樣點對應(yīng)的響應(yīng)值,為 n×1階向量;β為線性回歸部分的待定參數(shù),通過最優(yōu)線性無偏估計可以求得：

式中：R為n個采樣點相關(guān)系數(shù)組成的n×n階相關(guān)矩陣：

方差2σ 估計值用(10)式確定：

替代模型的建立可以通過求解上面的非線性無約束優(yōu)化問題[12]來實現(xiàn).待定參數(shù)θk求出后,通過建立的替代模型可獲得待求響應(yīng)值. θk可以通過無約束優(yōu)化(11)式求得[13-14]：

圖3 替代模型擬合相對誤差直方圖Fig.3 Histogram of relative error of the surrogate model

由圖 3可知替代模型輸出的污染物濃度與檢驗樣本的相對誤差均小于 1%,誤差較小,故認(rèn)為替代模型精度較高,可以代替模擬模型求解.

再次利用拉丁超立方方法對滲透系數(shù)和縱向彌散系數(shù)分別抽樣1000組,輸入克里格替代模型代替模擬模型求解,得到井1、2、3的各1000組污染物濃度值.對井1、2、3的各1000組污染物濃度值進行統(tǒng)計與分析.

2 結(jié)果與討論

2.1 觀測井污染預(yù)測與分析

2.1.1 污染物濃度統(tǒng)計分析 對 3口觀測井的各1000組污染物濃度值各項統(tǒng)計指標(biāo)進行統(tǒng)計,繪制3口觀測井污染濃度直方圖.結(jié)合圖4和表2可知,井2污染物波動范圍較大,分布較井1,3分散,可能原因是井2位于地下水流向下方,與污染源距離較遠(yuǎn)且相對來講靠近定水頭邊界,這些因素的綜合影響導(dǎo)致井 2中污染物濃度波動較大.且濃度與井1,3相比偏高,原因可能井2位于地下水流向一側(cè),污染質(zhì)在水流攜帶作用下運移加快,導(dǎo)致其濃度偏高.井 1,3與隔水邊界距離相等,與污染源距離有微小差距.故2口井污染濃度相差較小,波動區(qū)間也比較相近.

2.1.2 污染超標(biāo)風(fēng)險預(yù)警 為了對觀測井中污染物濃度進行更具體的統(tǒng)計分析,應(yīng)用偏度·峰度檢驗法[15]對替代模型輸出濃度值服從分布進行檢驗.偏度·峰度檢驗法是檢驗數(shù)據(jù)是否服從正態(tài)分布的一種數(shù)學(xué)方法,方法原理是利用數(shù)據(jù)樣本計算出樣本峰度,偏度與正態(tài)分布總體峰度,偏度的近似程度來判斷數(shù)據(jù)是否來自正態(tài)分布總體.

圖4 各觀測井污染物濃度頻數(shù)直方圖Fig.4 Frequency histogram of pollutant concentration of observation wells

表2 各觀測井污染物輸出值統(tǒng)計指標(biāo)(mg/L)Table 2 Statistic indexes of pollutant concentration of observation wells (mg/L)

圖5 污染物濃度概率密度擬合曲線圖Fig.5 The probability density fitting curve of pollutant concentration

利用峰度·偏度檢驗法檢驗出3口觀測井中污染物濃度值在滿足顯著水平為α=1%、置信程度為 1-α=99%時來自正態(tài)分布總體.井 1、2、3污染物濃度值的擬合情況見圖5：

在評價井中地下水遭受污染風(fēng)險時可以對1000組污染物濃度值進行正態(tài)分布擬合,求出擬合的正態(tài)概率密度函數(shù),對正態(tài)概率密度函數(shù)進行積分求出地下水受污染風(fēng)險,公式見(12)：

式中：P是觀測井污染風(fēng)險大小;f(x)是污染物濃度值服從的概率密度函數(shù);σ為正態(tài)分布的標(biāo)準(zhǔn)差;μ是正態(tài)分布的均值;m是污染物超標(biāo)標(biāo)準(zhǔn).將污染物濃度超過 120mg/L視為輕度污染,超過 200mg/L 視為較嚴(yán)重污染,超過280mg/L視為重度污染,利用(12)式進行積分求出井1、2、3為輕度污染風(fēng)險分別為99.94%、98.47%、74.50%;為較嚴(yán)重和嚴(yán)重污染風(fēng)險均為0%,根據(jù)風(fēng)險評價結(jié)果可知3口井相比井3受污染風(fēng)險較小,相比于井1和井2,井3更適合做水源.

對正態(tài)概率密度函數(shù)積分可以得到井1、2、3中污染物濃度在置信程度分別為80%、60%、40%, 20%時地下水污染物濃度波動區(qū)間,見表3.

表3 不同置信程度污染物濃度分布區(qū)間(mg/L)Table 3 The interval estimation of pollutant concentration of observation wells under different confidence (mg/L)

2.2 研究區(qū)污染預(yù)測與分析

2.2.1 污染羽范圍分布預(yù)測與分析 為預(yù)測研究區(qū)污染羽分布范圍,利用替代模型輸出研究區(qū)220個剖分網(wǎng)格的各1000組污染物濃度的均值繪制不同程度污染的污染羽分布圖.污染物濃度超過 120mg/L的區(qū)域為輕度污染區(qū)域,超過200mg/L區(qū)域為較嚴(yán)重污染區(qū)域,超過 280mg/L為嚴(yán)重污染區(qū)域,利用均值繪制圖6,根據(jù)圖6可知 15年后研究區(qū)受污染總面積約為 175000m2占研究區(qū)總面積的 32%,其中嚴(yán)重污染覆蓋面積最約為90000m2占研究區(qū)總面積16%.

圖6 污染羽分布Fig.6 The distribution of contaminant plume

根據(jù)研究區(qū)被污染面積大小,決策者可以提前準(zhǔn)備污染治理工作,根據(jù)研究區(qū)發(fā)展情況規(guī)劃需要治理的污染面積.

2.2.2 研究區(qū)污染分級風(fēng)險預(yù)警 為對研究區(qū)遭受污染風(fēng)險進行預(yù)警,統(tǒng)計替代模型輸出的研究區(qū)220個剖分網(wǎng)格的各1000組污染物濃度輸出值服從的分布規(guī)律,根據(jù)統(tǒng)計結(jié)果,分析各剖分網(wǎng)格受不同程度污染風(fēng)險,對研究區(qū)不同程度污染的不同污染超標(biāo)風(fēng)險對應(yīng)污染羽分布進行預(yù)測.由污染超標(biāo)風(fēng)險預(yù)警圖(圖 7～圖 9)可知在污染超標(biāo)風(fēng)險相同情況下,污染羽覆蓋范圍隨污染級別增大反而有減小態(tài)勢.對于相同程度污染,隨污染風(fēng)險置信程度增高,污染羽分布有減小態(tài)勢,若選擇污染風(fēng)險為 99%時的區(qū)域作為重點污染治理區(qū),需要治理面積可由污染超標(biāo)風(fēng)險預(yù)警圖計算得到.根據(jù)資金,設(shè)備,人員情況,也可以選擇相信其它污染風(fēng)險對應(yīng)的污染面積進行治理.

除此之外,可以根據(jù)研究區(qū)能夠承受的污染風(fēng)險上限對污染問題進行提前防治工作.也可以將污染超標(biāo)風(fēng)險預(yù)警圖與研究區(qū)內(nèi)人口密度及生活情況相結(jié)合,確定更經(jīng)濟更清潔的水源位置.結(jié)合污染風(fēng)險預(yù)警圖可以看出相同污染程度下,井3受污染風(fēng)險最小,與井1、2比較,井3更適合做水源.

圖7 輕度污染風(fēng)險預(yù)警Fig.7 The early warning of Light pollution risk

圖8 較嚴(yán)重污染超標(biāo)風(fēng)險預(yù)警Fig.8 The early warning of moderate pollution risk

圖9 嚴(yán)重污染超標(biāo)風(fēng)險預(yù)警Fig.9 The early warning of Serious pollution risk

3 結(jié)論

3.1 運用靈敏度分析方法篩選出滲透系數(shù)與縱向彌散系數(shù)作為數(shù)值模擬模型中隨機變量,在保證一定精度前提下建立模擬模型的替代模型完成蒙特卡洛模擬過程,可大幅度減小計算負(fù)荷,對于工程領(lǐng)域的研究工作具有實踐意義.

3.2 利用概率密度函數(shù)可以對地下水遭受污染風(fēng)險與不同置信程度的污染物濃度范圍進行估計;污染羽分布圖可以對研究區(qū)不同等級污染的污染羽覆蓋面積進行估計,分級污染風(fēng)險預(yù)警圖可以對研究區(qū)不同污染等級的污染超標(biāo)風(fēng)險進行預(yù)警.

3.3 基于污染質(zhì)運移數(shù)值模擬不確定性分析的地下水污染超標(biāo)風(fēng)險研究與確定性的地下水污染問題研究相比,可以對地下水污染進行更加客觀的預(yù)測,為決策者提供更多選擇依據(jù),對于地下水污染研究有重大意義.

[1]Kuczera G, Kavetski D, Franks S, et al. Towards a Bayesian total error analysis of conceptual rainfall-runoff models：Characterising model error using storm-dependent parameters [J]. Journal of Hydrology, 2006,331(1)：161-177.

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

[3]束龍倉,陶玉飛,劉佩貴.考慮水文地質(zhì)參數(shù)不確定性的地下水補給量可靠度計算 [J]. 水利學(xué)報, 2008,(3)：346-350.

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

[5]苗添升,盧文喜,歐陽琦,等.地下水?dāng)?shù)值模擬的不確定性分析在水質(zhì)預(yù)測中的應(yīng)用 [J]. 水電能源科學(xué), 2016,(8)：20-23+44.

[6]吳吉春,陸 樂.地下水模擬不確定性分析 [J]. 南京大學(xué)學(xué)報(自然科學(xué)版), 2011,(3)：227-234.

[7]韓林山,李向陽,嚴(yán)大考.淺析靈敏度分析的幾種數(shù)學(xué)方法 [J].中國水運(下半月), 2008,(4)：177-178.

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

[9]陳 彥,吳吉春.含水層滲透系數(shù)空間變異性對地下水?dāng)?shù)值模擬的影響 [J]. 水科學(xué)進展, 2005,16(4)：482-487.

[10]Krige D G. A Statistical Approach to Some Mine Valuations and Allied Problems at the Wit—watersrandr [D]. Wilwatersrand：University of Wit—watersrand, 1951.

[11]吳義忠,陳立平.多領(lǐng)域物理系統(tǒng)的仿真優(yōu)化方法 [M]. 科學(xué)出版社, 2011.

[12]徐亞棟,錢林方,陳龍淼.基于Kriging方法的復(fù)合材料身管結(jié)構(gòu)近似分析 [J]. 中國機械工程, 2007,18(8)：988-990.

[13]王曉鋒,席 光,王尚錦.Kriging與響應(yīng)面方法在氣動優(yōu)化設(shè)計中的應(yīng)用 [J]. 工程熱物理學(xué)報, 2005,26(3)：423-425.

[14]張 崎.基于 Kriging方法的結(jié)構(gòu)可靠性分析及優(yōu)化設(shè)計 [D].大連理工大學(xué), 2005.

[15]田 禹.基于偏度和峰度的正態(tài)性檢驗 [D]. 上海交通大學(xué), 2012.

Risk prediction of groundwater pollution based on uncertainty analysis

LI Jiu-hui1,2, LU Wen-xi1,2*, CHANGZhen-bo1,2, LI Ming-yu1,2, MIAO Tian-sheng1,2, ZHAO Ying1,2, ZHANG Jiang-wei1,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(6)：2270～2277

Groundwater contamination prediction could be achieved through the numerical simulation model of groundwater contaminant transport.In order to analyze the influence of the uncertainty of parameters in the model, Monte Carlo method was used to analyze the uncertainty of the model output in this paper,In order to reduce the complexity of the numerical simulation model, the sensitivity analysis method was used to select the random variables in the model.To reduce the calculated load generated by the numerical simulation model called repeated, set up Krig surrogate model instead of the numerical simulation model to complete the simulation process. The results showed that the probability density function integral could be used to estimate the groundwater contamination risk and the contamination concentration range under different confidence level. The distribution of contaminant plume could estimate different level pollution cover acreage of the study area and the classification contamination risk early warning could estimate the contamination plume distribution of the study area in different contamination risk; Groundwater contamination risk warning, which was based on the uncertainty analysis of groundwater contaminant transport numerical simulation could forecast the groundwater problem more objectivily.

groundwater contamination；numerical simulation model；Monte Carlo；Krig surrogate models；uncertainty analysis

X523

A

1000-6923(2017)06-2270-08

李久輝(1992-),女,內(nèi)蒙古赤峰人,碩士研究生,主要從事污染質(zhì)運移數(shù)值模擬的不確定性研究.

2016-11-16

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

* 責(zé)任作者, 教授, luwenxi@jlu.edu.cn