歐陽琦,盧文喜*,侯澤宇,顧文龍,辛 欣(1.吉林大學(xué)地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,吉林長春 130021；2.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130021)

?

基于替代模型的地下水溶質(zhì)運(yùn)移不確定性分析

歐陽琦1,2,盧文喜1,2*,侯澤宇1,2,顧文龍1,2,辛 欣1,2(1.吉林大學(xué)地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,吉林長春 130021；2.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130021)

摘要：為分析參數(shù)的不確定性對地下水溶質(zhì)運(yùn)移數(shù)值模型的影響,采用蒙特卡羅(Monte Carlo)模擬對一算例進(jìn)行分析,并從風(fēng)險(xiǎn)評估的角度對不確定性分析的結(jié)果進(jìn)行了闡釋.為減小計(jì)算負(fù)荷,利用Sobol’法對模型參數(shù)進(jìn)行了靈敏度分析,篩選出較為敏感的參數(shù)作為隨機(jī)變量,建立了模擬模型的克里格(Kriging)替代模型,進(jìn)而實(shí)現(xiàn)Monte Carlo模擬.結(jié)果表明:置信度為80%時(shí),井1,2,3濃度值的置信區(qū)間分別為23.46～42.06,47.99～66.73,69.54～82.94mg/L;結(jié)合風(fēng)險(xiǎn)評估,計(jì)算出地下水受污染的風(fēng)險(xiǎn)為0.54,可為地下水污染物防控與修復(fù)提供科學(xué)依據(jù).

關(guān)鍵詞：地下水溶質(zhì)運(yùn)移；蒙特卡羅；不確定性分析；替代模型；風(fēng)險(xiǎn)評估；克里格

* 責(zé)任作者, 教授, Luwx999@163.com

地下水溶質(zhì)運(yùn)移數(shù)值模型可以模擬污染物在地下水中的遷移轉(zhuǎn)化規(guī)律,預(yù)測污染物的分布情況,是地下水環(huán)境管理者常用的科技工具.由于地下水系統(tǒng)本身的復(fù)雜性,模擬時(shí)對模擬系統(tǒng)的簡化及其他一些無法預(yù)測的因素[1-2],在利用地下水溶質(zhì)運(yùn)移數(shù)值模型時(shí)存在著諸多的不確定性因素.束龍倉等[3]根據(jù)不確定性因素產(chǎn)生的原因,將其劃分為客觀不確定性因素和主觀不確定性因素兩大類.由于確定性模型只能得到唯一的結(jié)果[4],單純依靠確定性模型進(jìn)行地下水污染物運(yùn)移模擬,其結(jié)果的可靠性有待驗(yàn)證.不確定性分析對于提高模型的可靠性具有十分重要的意義,還可以滿足管理者對于模型極端情況預(yù)測的需求[5].

Monte Carlo(蒙特卡羅)模擬又稱統(tǒng)計(jì)試驗(yàn)法,是國內(nèi)外常見的分析地下水?dāng)?shù)值模擬不確定性的方法[4,6-7].它適用性廣,方法簡單,能將模型參數(shù)的不確定性轉(zhuǎn)化為模擬結(jié)果的不確定性[4],適用于討論多變量的模型中不確定性的問題[8].在進(jìn)行Monte Carlo模擬時(shí)需要進(jìn)行大量的統(tǒng)計(jì)試驗(yàn)(通常是成百上千次),若直接調(diào)用模擬模型,計(jì)算成本比較大.替代模型是模擬模型輸入輸出響應(yīng)關(guān)系的代替,它能夠以較小的計(jì)算量得到和模擬模型相近的輸入輸出關(guān)系[9].因此基于替代模型的Monte Carlo模擬能夠大幅減小計(jì)算負(fù)荷,現(xiàn)階段文獻(xiàn)還未見相關(guān)應(yīng)用.

吳吉春等[6]認(rèn)為將地下水模擬不確定性分析結(jié)果與風(fēng)險(xiǎn)評估相結(jié)合,具有重大的實(shí)踐意義.目前將不確定性分析結(jié)果與風(fēng)險(xiǎn)評估相結(jié)合的做法多用于地下水資源評價(jià)方面,如束龍倉等[10],李如忠等[11]曾將不確定分析與風(fēng)險(xiǎn)評估相結(jié)合,用于地下水允許開采量的確定,而在地下水溶質(zhì)運(yùn)移模擬方面并不多見.本次研究借助替代模型進(jìn)行Monte Carlo模擬,又在地下水溶質(zhì)運(yùn)移模擬方面將不確定性分析的結(jié)果與風(fēng)險(xiǎn)評估相結(jié)合.首先根據(jù)一假想算例建立了地下水溶質(zhì)運(yùn)移數(shù)值模擬模型,然后利用Sobol’全局靈敏度分析方法對模型參數(shù)進(jìn)行靈敏度分析,篩選出較為敏感的參數(shù)作為隨機(jī)變量建立模擬模型的Kriging替代模型,進(jìn)而借助替代模型實(shí)現(xiàn)Monte Carlo模擬,最后將不確定性分析的結(jié)果與風(fēng)險(xiǎn)評估相結(jié)合.

1　研究方法

1.1 模型建立

建立了一個(gè)二維潛水含水層的溶質(zhì)運(yùn)移模型,用以模擬假想的實(shí)際含水層.采用確定性地下水流數(shù)值模型GMS求解該地下水溶質(zhì)運(yùn)移問題.假定該模型中含水層水平等厚,平面上為正方形,邊長為1000m,厚度為30m,網(wǎng)格剖分如圖1.溶質(zhì)運(yùn)移模型是在水流模型的基礎(chǔ)上建立的,首先定義水流模型,含水層西側(cè)和東側(cè)設(shè)定為定水頭邊界,水頭值分別為27m和24m,北部和南部為隔水邊界.對于溶質(zhì)運(yùn)移模型,西側(cè)為已知濃度邊界(濃度為0mg/L),其余邊界為零溶質(zhì)通量邊界.含水層概化為均質(zhì)各向同性,通常情況下,模型參數(shù)的先驗(yàn)分布難以確定,通常采用均勻分布來描述參數(shù)的先驗(yàn)分布[12],本文亦采用均勻分布來描述溶質(zhì)運(yùn)移模型參數(shù)的先驗(yàn)分布(參數(shù)取值情況見表1).

圖1　地下水溶質(zhì)運(yùn)移模型剖分和注水井,觀測井位置Fig.1　Model mesh of the groundwater solute transport and the location of the injection and observation wells

表1　參數(shù)取值范圍Table 1　Ranges of model parameters

現(xiàn)在模擬區(qū)內(nèi)加入1口注水井以及3口抽水井(同時(shí)也作為觀測井),注水井流量為100m3/d,抽水井1,2,3流量分別為100,100, 50m3/d,且注水井以200mg/L的濃度向含水層持續(xù)排放污染物,各井位置見圖1.含水層在垂向上接受大氣降水補(bǔ)給,區(qū)域年降水量設(shè)為450mm(降水量年內(nèi)均勻分布).含水層主要排泄方式為人工抽水,蒸發(fā)排泄可忽略.模擬時(shí)間為12年,模型運(yùn)行4,8,12年后污染物運(yùn)移情況見圖2(本研究以模型運(yùn)行12年后的結(jié)果進(jìn)行計(jì)算).

圖2　模型污染物運(yùn)移情況Fig.2　Transport of pollutant in model (a):4年,(b):8年,(c):12年

1.2 Sobol’全局靈敏度分析

靈敏度分析是一種研究系統(tǒng)各種輸入變化對輸出的影響程度的方法[13-14].它雖然可以單獨(dú)作為不確定性分析的一種方法,但是它對系統(tǒng)的輸出缺少統(tǒng)計(jì)學(xué)上的分析,因此本次研究只將其作為不確定性分析中的一個(gè)步驟.相對于局部靈敏度分析法,全局靈敏度分析法可以分析多個(gè)參數(shù)的變化對模型運(yùn)行結(jié)果總的影響,并能分析每一個(gè)參數(shù)及其參數(shù)之間相互作用對模型結(jié)果的影響[15],得到的結(jié)果更加接近于實(shí)際情況.Sobol’方法是目前最常用的基于方差的全局敏感性分析方法.其原理可參考文獻(xiàn)[8].

根據(jù)Sobol’法的原理在MATLAB上編寫程序,計(jì)算了5個(gè)參數(shù)總靈敏度,各個(gè)參數(shù)的總靈敏度包括了多個(gè)參數(shù)的協(xié)同作用.各變量的總敏感度結(jié)果如圖3所示.

敏感度分析結(jié)果表明:滲透系數(shù),降水入滲系數(shù)以及給水度對模型輸出的影響較大,孔隙度以及縱向彌散度影響較小.因此,可將滲透系數(shù),降水入滲系數(shù)以及給水度作為隨機(jī)變量,著重考察這3個(gè)參數(shù)的不確定性對模型輸出的影響,而將孔隙度,縱向彌散度取其平均值作為確定值.這樣既降低了替代模型的輸入維度,也減小了不確定性分析的工作負(fù)荷.

圖3　各參數(shù)總靈敏度計(jì)算成果Fig.3　Total sensitivity of each parameter

1.3 模擬模型的Kriging替代模型

替代模型是模擬模型輸入輸出響應(yīng)關(guān)系的代替,它能夠以較小的計(jì)算量得到和模擬模型相近的輸入輸出關(guān)系.本研究中的輸入即是通過靈敏度分析篩選出來的較為敏感的參數(shù)的取值(即滲透系數(shù),降水入滲系數(shù)以及給水度),輸出則是模型在不同參數(shù)取值情況下運(yùn)行12年后,各觀測井污染物的濃度值.

Luo等[16]、安永凱等[17]認(rèn)為,Kriging(克里格)法是建立地下水?dāng)?shù)值模擬模型的替代模型的有效方法.因此本次研究采用該方法建立替代模型.

1.3.1 替代模型的建立 首先利用拉丁超立方抽樣方法(LHS)在隨機(jī)變量的取值范圍內(nèi)(表1)進(jìn)行抽樣,共抽取70組參數(shù)作為輸入值,然后將其帶入到模擬模型中獲得各個(gè)樣本所對應(yīng)的輸出值,輸入與輸出構(gòu)成一組樣本.選取前60組樣本作為訓(xùn)練樣本,利用Kriging法建立替代模型, 后10組樣本作為檢驗(yàn)樣本對所建立的替代模型的精度進(jìn)行檢驗(yàn).替代模型檢驗(yàn)情況如圖4所示.

從圖4來看,各個(gè)井的大部分相對誤差都非常小(小于0.005),只有極個(gè)別誤差稍大,最大的也不超過0.02.說明替代模型的精度非常高,能充分逼近模擬模型的輸入輸出關(guān)系,其結(jié)果是可信的.

1.3.2 Monte Carlo模擬 Monte Carlo模擬假定隨機(jī)變量的概率分布函數(shù)已知, 通過隨機(jī)抽樣方法得到多組輸入變量,每組變量相當(dāng)于一次統(tǒng)計(jì)試驗(yàn),然后把隨機(jī)變量帶入模型中從而得到大量的輸出,通過對輸出(如地下水模擬中的水頭或溶質(zhì)濃度)的統(tǒng)計(jì)可以得到均值,方差等統(tǒng)計(jì)估計(jì)量,并擬合其輸出結(jié)果的概率分布情況.

圖4　替代模型檢驗(yàn)樣本相對誤差箱線圖Fig.4　Box plot of relative error of the surrogate model sample

本次研究利用拉丁超立方抽樣法(LHS)抽取1000組參數(shù)樣本,將其帶入建立好的替代模型,得到1000組響應(yīng)輸出(即污染物濃度值),然后對這1000組響應(yīng)輸出進(jìn)行不確定性分析.

2　結(jié)果與討論

2.1 統(tǒng)計(jì)分析

首先,畫出3口觀測井污染物濃度的頻數(shù)直方圖,如圖5所示,以直觀地分析輸出結(jié)果的不確定性.然后分別對各個(gè)觀測井的輸出值的統(tǒng)計(jì)指標(biāo)進(jìn)行計(jì)算(表2).

由圖4以及表2可以看出,3口井的輸出結(jié)果較為分散,從1號井到3號井,極差和標(biāo)準(zhǔn)差有所增大,說明輸出的結(jié)果越來越分散,可能與觀測井距污染源的位置越來越遠(yuǎn)有關(guān),離污染源越近輸出的不確定性越小,分布越集中.可能的原因是,距離污染源越近溶質(zhì)運(yùn)移的距離越短,受含水層各參數(shù)的影響就越小.

利用SPSS軟件的K-S檢驗(yàn)分別對3口井的輸出值進(jìn)行分布檢驗(yàn)(包括指數(shù)分布,均勻分布以及正態(tài)分布).檢驗(yàn)結(jié)果顯示,井3的輸出近似服從均值為76.24,標(biāo)準(zhǔn)差為5.19的正態(tài)分布;1,2號井的輸出均不服從以上分布,因此暫且認(rèn)為其分布未知.

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

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

2.2 區(qū)間估計(jì)

區(qū)間估計(jì)在概率學(xué)中指在一定置信度下對未知參數(shù)真值存在范圍的估計(jì).本文將其引申,把各個(gè)觀測井的輸出值當(dāng)做是要估計(jì)的參數(shù),計(jì)算其在不同置信水平下的區(qū)間范圍.

切比雪夫不等式[18]是俄國的數(shù)學(xué)家Chebyshev在研究隨機(jī)變量的統(tǒng)計(jì)規(guī)律時(shí),利用標(biāo)準(zhǔn)差構(gòu)建的一個(gè)不等式.式(7)是它的一種形式.

利用它可以計(jì)算出不同置信水平下的區(qū)間范圍.它對于分布未知,而已知平均值和方差的情況比較適用,但是其估計(jì)結(jié)果比較粗糙.對于已知分布的隨機(jī)變量,可以通過其概率密度函數(shù)來進(jìn)行確切的計(jì)算.

圖6　各井污染風(fēng)險(xiǎn)評估Fig.6　Map of contamination risk assessment of each observation well圖中曲線上的某點(diǎn)(x,y)表示濃度大于x的風(fēng)險(xiǎn)(即概率)是y

表3　不同置信水平下各觀測井濃度值的區(qū)間估計(jì)(mg/L)Table 3　Concentration interval estimation of each observation well under different confidence (mg/L)

因此,對1,2號井輸出結(jié)果采用切比雪夫不等式進(jìn)行區(qū)間估計(jì),對3號井輸出結(jié)果利用其概率密度函數(shù)進(jìn)行區(qū)間估計(jì).由表3可以看出,置信水平越高,區(qū)間范圍越大;置信水平越低,區(qū)間范圍越小,也越集中在均值附近.

2.3 風(fēng)險(xiǎn)評估假設(shè)當(dāng)同時(shí)滿足井1,井2,井3的濃度分別大于30,60,70mg/L時(shí),可認(rèn)為研究區(qū)地下水已受到污染.根據(jù)做出各個(gè)觀測井輸出值的累積分布函數(shù)畫出風(fēng)險(xiǎn)評估圖(圖6).從圖6中可以看出,井1濃度大于30mg/L的風(fēng)險(xiǎn)是0.68,井2濃度大于55mg/L的風(fēng)險(xiǎn)是0.66,井3濃度大于75mg/L的風(fēng)險(xiǎn)是0.54.因此,該研究區(qū)地下水受到污染的風(fēng)險(xiǎn)是0.54,風(fēng)險(xiǎn)較大.

3　結(jié)論

3.1 靈敏度分析結(jié)果顯示:滲透系數(shù),降水入滲系數(shù)以及給水度對模型輸出的影響較大,孔隙度以及縱向彌散度影響較小.通過靈敏度分析可以篩選出對模型輸出影響較大的參數(shù),能減少所考察的變量的個(gè)數(shù),減小了不確定性分析的工作負(fù)荷,且降低了替代模型的輸入維數(shù).

3.2 利用Kriging法建立的替代模型精度較高,能充分逼近模擬模型的輸入輸出關(guān)系.因此基于替代模型的Monte Carlo模擬是可行的.

3.3 由于不確定性是客觀存在的,相對于確定的輸出,一定置信水平下的取值區(qū)間往往更有指導(dǎo)意義,切比雪夫不等式可以很好地應(yīng)用于觀測井濃度的區(qū)間估計(jì).

參考文獻(xiàn)：

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

[2] Feyen L, Caers L. Quantifying geological uncertainty for flow and transport modeling in multi-modal heterogeneous formations [J]. Advances in Water Resources, 2006,29:912-929.

[3] 束龍倉,朱元生.地下水資源評價(jià)中的不確定性因素分析 [J].水文地質(zhì)工程地質(zhì), 2000,(6):6-8.

[4] 陸 樂,吳吉春.地下水?dāng)?shù)值模擬不確定性的貝葉斯分析 [J].水利學(xué)報(bào), 2010,41(3):264-271.

[5] 劉悅憶,趙建世,黃躍飛,等.基于蒙特卡洛模擬的水質(zhì)概率預(yù)報(bào)模型 [J]. 水利學(xué)報(bào), 2015,46(1):51-57.

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

[7] Wu C M, Yeh T C J, Zhu J F, et al. Traditional analysis of comparing apples to oranges [J]. Water Resources Research, 2005,41(9)W09402.

[8] 張質(zhì)明,王曉燕,李明濤.基于全局敏感性分析方法的WASP模型不確定性分析 [J]. 中國環(huán)境科學(xué), 2014,34(5):1336-1346.

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

[10] 束龍倉,朱元生,孫慶義,等.地下水允許開采量確定的風(fēng)險(xiǎn)分析[J]. 水利學(xué)報(bào), 2000,(3):77-81.

[11] 李如忠,汪家權(quán),錢家忠.地下水允許開采量的未確知風(fēng)險(xiǎn)分析[J]. 水利學(xué)報(bào), 2004,(4):54-60.

[12] 林 青,徐紹輝.基于GLUE方法的飽和多孔介質(zhì)中溶質(zhì)運(yùn)移模型參數(shù)不確定性分析 [J]. 水利學(xué)報(bào), 2012,43(9):1017-1024.

[13] 周 振,喬衛(wèi)敏,蔣路漫,等.基于數(shù)學(xué)模型的AAO系統(tǒng)階躍響應(yīng)特性分析 [J]. 中國環(huán)境科學(xué), 2015,35(2):442-447.

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

[15] 張 巍,鄭一,王學(xué)軍.水環(huán)境非點(diǎn)源污染的不確定性及分析方法 [J]. 農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào), 2008,27(4):1290-1296.

[16] Luo J, Lu W. Comparison of surrogate models with different methods in groundwater remediation process [J]. J. Earth Syst. Sci., 2014,123(7):1579–1589.

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

[18] 盛 驟,謝式千,潘承毅.概率論與數(shù)理統(tǒng)計(jì)(第四版) [M]. 北京:高等教育出版社, 2011.

主要從事地下水?dāng)?shù)值模擬與優(yōu)化管理方面的研究.發(fā)表論文2篇.

Uncertainty analysis of groundwater solute transport based on surrogate model.

OUYANG Qi1,2, LU Wen-xi1,2*, HOU Ze-yu1,2, GU Wen-long1,2, XIN Xin1,2(1.Key Laboratory of Groundwater Resources and Environment, Ministry of Education, Jilin University, Changchun 130021, China；2.College of Environment and Resources, Jilin University, Changchun 130021, China). China Environmental Science, 2016,36(4)：1119～1124

Abstract：In order to analyze the influence that parameter uncertainty has on groundwater solute transport numerical simulation, this study adopted Monte Carlo simulation to conduct an uncertainty analysis on an example and illustrated its results from the perspective of risk assessment. For the purpose of reducing computation load, Sobol’ method was used to analyze the sensitivity of model parameters, which helped to select the more sensitive parameters as random variables and to construct the Kriging surrogate model of the simulation model. This surrogate model provided further help in achieving the Monte Carlo simulation. Results showed: when the confidence was 80%, the confidence intervals of well 1, 2, 3 were 23.46～42.06,47.99～66.73,69.54～82.94mg/L, respectively. Combined with risk assessment, the risk of groundwater contamination was calculated as 0.54, which could serve as a scientific guide for the prevention and control of groundwater pollution.

Key words：groundwater solute transport；Monte Carlo；uncertainty analysis；surrogate model；risk assessment；Kriging

作者簡介：歐陽琦(1991-),男,湖南洞口縣人,吉林大學(xué)博士研究生,

基金項(xiàng)目：國家自然科學(xué)基金(41072171)

收稿日期：2015-09-21

中圖分類號：X523

文獻(xiàn)標(biāo)識碼：A

文章編號：1000-6923(2016)04-1119-06