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

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

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

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

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

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

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

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

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

1 研究方法

1.1 蒙特卡洛方法

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

1.2 模型建立

建立一個(gè)東西長 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ū)在垂向上均勻的接受大氣降水補(bǔ)給,年平均降水量為580mm;在忽略蒸發(fā)騰散消耗量前提下,地下水排泄方式主要是人工開采,井1、2,3抽水量均為500m3/d.區(qū)內(nèi)一化工廠污水未經(jīng)處理泄漏到地下含水層中,污水污染物濃度為 400mg/L,泄漏量為600m3/d.將污染物視為不會(huì)發(fā)生化學(xué)變化及生物轉(zhuǎn)化與遷移的穩(wěn)定污染物.潛水含水層中污染物的初始濃度為零,定水頭邊界Г1,Г3與隔水邊界Г2,Г4分別視為零濃度邊界和零通量邊界.根據(jù)研究區(qū)概況,為更好地對(duì)研究區(qū)污染質(zhì)運(yùn)移情況進(jìn)行研究,利用GMS將研究區(qū)剖分離散為220個(gè)有限差分網(wǎng)格,區(qū)內(nèi)剖分情況及觀測(cè)井與污染源位置詳見圖1：

圖1 地下水溶質(zhì)運(yùn)移模型剖分及污染源,觀測(cè)井位置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為人工開采強(qiáng)度,m/d;R為降水入滲補(bǔ)給量,m/d;μ為含水層給水度,無量綱;S為模擬區(qū)范圍;Г1,Г3為定水頭邊界;Г2,Г4為零流量邊界;Kn為邊界法向量上的滲透系數(shù),m.

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

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

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

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

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

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

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

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

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

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

1.4 拉丁超立方方法

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

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

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

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

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

1.5 建立替代模型

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

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

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

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

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

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

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

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

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

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

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

2 結(jié)果與討論

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

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

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

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

表2 各觀測(cè)井污染物輸出值統(tǒng)計(jì)指標(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

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

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

式中：P是觀測(cè)井污染風(fēng)險(xiǎn)大小;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)式進(jìn)行積分求出井1、2、3為輕度污染風(fēng)險(xiǎn)分別為99.94%、98.47%、74.50%;為較嚴(yán)重和嚴(yán)重污染風(fēng)險(xiǎn)均為0%,根據(jù)風(fēng)險(xiǎn)評(píng)價(jià)結(jié)果可知3口井相比井3受污染風(fēng)險(xiǎn)較小,相比于井1和井2,井3更適合做水源.

對(duì)正態(tài)概率密度函數(shù)積分可以得到井1、2、3中污染物濃度在置信程度分別為80%、60%、40%, 20%時(shí)地下水污染物濃度波動(dòng)區(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ù)測(cè)與分析

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

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

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

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

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

3 結(jié)論

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

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

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

[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]曾獻(xiàn)奎,王 棟,吳吉春.地下水流概念模型的不確定性分析 [J].南京大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012,(6)：746-752.

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

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

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

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

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

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

[9]陳 彥,吳吉春.含水層滲透系數(shù)空間變異性對(duì)地下水?dāng)?shù)值模擬的影響 [J]. 水科學(xué)進(jìn)展, 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]. 中國機(jī)械工程, 2007,18(8)：988-990.

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

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

[15]田 禹.基于偏度和峰度的正態(tài)性檢驗(yàn) [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ì)運(yùn)移數(shù)值模擬的不確定性研究.

2016-11-16

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

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