基于替代模型的地下水溶質運移不確定性分析

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

?

基于替代模型的地下水溶質運移不確定性分析

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

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

關鍵詞：地下水溶質運移；蒙特卡羅；不確定性分析；替代模型；風險評估；克里格

* 責任作者, 教授, Luwx999@163.com

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

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

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

1　研究方法

1.1 模型建立

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

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

表1　參數取值范圍Table 1　Ranges of model parameters

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

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

1.2 Sobol’全局靈敏度分析

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

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

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

圖3　各參數總靈敏度計算成果Fig.3　Total sensitivity of each parameter

1.3 模擬模型的Kriging替代模型

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

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

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

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

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

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

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

2　結果與討論

2.1 統計分析

首先,畫出3口觀測井污染物濃度的頻數直方圖,如圖5所示,以直觀地分析輸出結果的不確定性.然后分別對各個觀測井的輸出值的統計指標進行計算(表2).

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

利用SPSS軟件的K-S檢驗分別對3口井的輸出值進行分布檢驗(包括指數分布,均勻分布以及正態分布).檢驗結果顯示,井3的輸出近似服從均值為76.24,標準差為5.19的正態分布;1,2號井的輸出均不服從以上分布,因此暫且認為其分布未知.

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

表2　各觀測井污染物輸出值統計指標(mg/L)Table 2　Output statistic indexes of observation wells (mg/L)

2.2 區間估計

區間估計在概率學中指在一定置信度下對未知參數真值存在范圍的估計.本文將其引申,把各個觀測井的輸出值當做是要估計的參數,計算其在不同置信水平下的區間范圍.

切比雪夫不等式[18]是俄國的數學家Chebyshev在研究隨機變量的統計規律時,利用標準差構建的一個不等式.式(7)是它的一種形式.

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

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

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

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

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

3　結論

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

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

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

參考文獻：

[1] 申 升.基于貝葉斯理論的地下水溶質運移的不確定性研究[D]. 長沙:湖南大學, 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] 束龍倉,朱元生.地下水資源評價中的不確定性因素分析 [J].水文地質工程地質, 2000,(6):6-8.

[4] 陸 樂,吳吉春.地下水數值模擬不確定性的貝葉斯分析 [J].水利學報, 2010,41(3):264-271.

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

[6] 吳吉春,陸 樂.地下水模擬不確定性分析 [J]. 南京大學學報(自然科學), 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] 張質明,王曉燕,李明濤.基于全局敏感性分析方法的WASP模型不確定性分析 [J]. 中國環境科學, 2014,34(5):1336-1346.

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

[10] 束龍倉,朱元生,孫慶義,等.地下水允許開采量確定的風險分析[J]. 水利學報, 2000,(3):77-81.

[11] 李如忠,汪家權,錢家忠.地下水允許開采量的未確知風險分析[J]. 水利學報, 2004,(4):54-60.

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

[13] 周 振,喬衛敏,蔣路漫,等.基于數學模型的AAO系統階躍響應特性分析 [J]. 中國環境科學, 2015,35(2):442-447.

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

[15] 張 巍,鄭一,王學軍.水環境非點源污染的不確定性及分析方法 [J]. 農業環境科學學報, 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] 安永凱,盧文喜,董海彪,等.基于克里格法的地下水流數值模擬模型的替代模型研究 [J]. 中國環境科學, 2014,34(4):1073-1079.

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

主要從事地下水數值模擬與優化管理方面的研究.發表論文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-),男,湖南洞口縣人,吉林大學博士研究生,

基金項目：國家自然科學基金(41072171)

收稿日期：2015-09-21

中圖分類號：X523

文獻標識碼：A

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