考虑参数不确定性的地下水污染源识别

摘要：为分析参数不确定性对地下水污染源识别的影响,本文通过模拟-优化方法、灵敏度分析方法、蒙特卡罗方法和克里格方法的综合运用,建立了描述渗透系数与污染物质释放强度之间关系的推算模型,进行了考虑参数不确定性的地下水污染源识别研究.研究结果表明,推算模型具有较高的精度,确定性系数和平均相对误差分别为0.9895和4.51%;运用推算模型推算了8000组渗透系数影响下的污染源识别结果,节省了约99%的计算负荷和时间;对8000组污染源识别结果进行了定量的统计与分析,得到了概率密度最大的污染源识别结果和置信水平分别为80%、60%、40%和20%对应的污染源识别结果置信区间.本研究改善了应用模拟-优化方法进行地下水污染源识别时,难以考虑参数不确定性的缺点,可以为决策者提供更多的参考依据.

关键词：地下水污染源；模拟-优化；不确定性；替代模型；推算模型

地下水污染具有存在的隐蔽性和发现的滞后性特点,致使人们对于地下水污染源的特征缺乏了解和掌握.这给地下水污染修复方案的合理设计、污染责任认定和污染风险评估都带来了很大的困难[1-3].因此,关于地下水污染源识别的研究就显得格外重要.

地下水污染源识别兴起于20世纪80年代,发展到今天应用于地下水污染源识别的方法包括直接方法、概率和地统计模拟方法、模拟-优化方法和地球物理探测法等[4].其中,模拟-优化方法,近些年来被广泛应用于地下水污染源识别[5-7].江思珉等[8]运用单纯形模拟退火混合算法求解优化模型,识别地下水污染源强度.随后,又将模拟-优化方法与卡尔曼滤波方法结合,进行基于污染羽形态对比的地下水污染源识别研究[9].肖传宁等[10]应用基于径向基函数替代模型的模拟-优化方法,识别地下水污染源的污染物质泄漏量.侯泽宇等[11]应用基于核极限学习机替代模型的模拟-优化方法,对地下水重非水相流体(DNAPLs)污染源及含水层参数的进行同步识别.

尽管应用模拟-优化方法进行地下水污染源识别,取得了丰硕的研究成果.但是,应用模拟-优化方法进行地下水污染源识别研究,只会得到唯一的污染源识别结果(某一组参数取值影响下的污染源特征)[12-13].因此,基于模拟-优化方法,考虑参数不确定性的污染源识别研究,实施起来尤为困难.而参数的不确定性客观存在[14-15],会影响污染源的识别结果.

本研究提出将模拟-优化方法、灵敏度分析方法、蒙特卡罗方法和克里格方法结合,进行考虑参数不确定性的地下水污染源识别研究.首先根据研究区的具体条件,建立地下水污染物质运移数值模拟模型.运用灵敏度分析方法筛选出对模拟模型输出结果影响最大的参数.为减少调用模拟模型耗费的计算时间和负荷,基于克里格方法建立了模拟模型的替代模型.然后,确定决策变量、目标函数和约束条件,建立识别污染源的优化模型.将替代模型做为等式约束连接到优化模型中.对筛选出的参数抽样90组,把所有参数都依次做为约束条件赋值到优化模型中并求解优化模型,得到所有参数影响下的污染源识别结果.最后,基于多组参数及其影响下的污染源特征,利用克里格方法建立描述参数与污染源特征之间关系的推算模型,运用推算模型计算成千上万组参数影响下的污染源识别结果(蒙特卡罗模拟),并对识别结果进行定量的统计与分析.

1 研究方法

1.1 局部灵敏度分析方法

局部灵敏度分析方法,是用来筛选模型参数对模型输出结果影响大小的方法.它的原理是利用模型输出结果对输入模型的参数求偏导数,利用偏导数大小,来判断输入模型的参数对模型输出结果影响程度的大小(式1),为了方便计算,式1可以转换为式2:

式中:Sk是灵敏度系数;xk是输入模型的第k个参数;xk是输入模型的第k个参数的变化量;是参数变化时,模型输出结果;是参数为xk时模型输出结果;n是区域观测井的数量.灵敏度系数越大,说明该参数对模型输出结果影响越大.式2中符号单位根据具体分析参数而定.

1.2 克里格方法

克里格方法是地统计学的一种插值方法.近些年来,克里格方法被延伸为一种建立替代模型的方法,被应用于多个工程领域[16-17].克里格方法的原理如下:

式中:qk为待定参数;和分别是第i和第j个样本的k维取值.

b基函数的待定参数,可以通过最优线性无偏估计可以求得:

1.3 蒙特卡罗方法

蒙特卡罗方法,又称随机抽样或统计试验方法.蒙特卡罗方法的基本思想是通过“试验”的方法,统计某个随机事件发生的频率,或是某个随机变量的一些数字特征,以这种事件出现的频率估计这一随机事件的概率,或者将某些数字特征,作为随机事件的解.

蒙特卡罗方法的实现一般包括以下几步:(1)基于需要解决的问题,构造易于实现的概率统计模型或者随机过程.(2)确定输入模型的随机变量和随机变量抽样方法.(3)基于概率统计模型或随机过程,进行多次模拟试验,对试验的结果进行统计与分析,将某一结果的发生频率或者某些数字特征(包括均值和方差、标准差等)作为问题的解.蒙特卡罗方法的基本原理,详见尹增谦等[18]和朱辉等[19]文章.

2 案例研究

2.1 研究区概况

图1 研究区域概况

Fig.1 Overview of the study area

S1和S2分别代表第1和第2个污染源;O1和O7分别代表第1到第7口观测井

本文借鉴文献[4,20]的案例进行研究.研究区是具有5个参数分区,边界不规则的二维非均质各项同性承压含水层,地下水流为非稳定流,水流方向由AB边界指向CD边界(图1);AB和CD边界是已知水头边界,水头分别为100 和80m;AC和BD边界为隔水边界;研究区在垂直方向上接受均匀的水量补给,补给率为0.0000864m/d.含水层的水文地质参数见表1.地下水中污染物质的初始浓度为100mg/L.污染物质迁移的总模拟时间为10a,共有20个模拟期(每6个月为1个模拟期,每月30d).假设污染物质是不经过生物转化或化学变化的保守污染物.污染源的位置已知,在第1个模拟期初至第4个模拟期末,污染物被持续释放到地下水中,在后来的模拟期不再释放污染物质(表2).区域有7口观测井.含水层中5个参数区的分布,观测井和污染源的位置见图1.

表1 含水层参数

Table 1 Aquifer parameters

表2 污染物质在各释放时段的释放强度

Table 2 Contaminant release intensity in each release periods

基于以上的研究案例,本研究的技术路线见图2.

2.2 数值模拟模型

根据研究区的具体条件,建立水文地质概念模型.在概念模型的基础上,建立水流和溶质运移数值模拟模型.描述二维承压含水层系统中非稳态流的地下水水流的控制偏微分方程如下:

式中:K是渗透系数,m/d;H是水头,m;W是水流模拟模型的源汇项,m/d;ms是贮水率或释水率,m-1;t是时间,d;x代表横向方向,y代表纵向方向(长度单位为m).

图2 技术路线图

Fig.2 Technology Roadmap

图3 水流空间分布

Fig.3 Spatial distribution map of water flow

描述溶质运移的控制偏微分方程如下:

(12)

式中:q是孔隙度,无量纲;C是污染物质浓度,mg/L;ux是横向水流实际平均流速(纵向水流实际平均流速与横向计算方法相同),m/d;D是弥散度,m2/d;R是溶质运移模拟模型的源汇项,mg/d.

图4 污染物质空间分布

Fig.4 Spatial distribution map of contaminant

达西定律可用于计算式12中的ui,如式13所示:

(13)

建立水流和溶质运移模拟模型后,使用GMS软件中MODFLOW和MT3DMS工具箱来模拟水流和污染物质运移的过程.