基于替代模型的地下水溶质运移不确定性分析

欧阳琦,卢文喜*,侯泽宇,顾文龙,辛 欣(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