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

摘要：为分析参数不确定性对地下水污染源识别的影响,本文通过模拟-优化方法、灵敏度分析方法、蒙特卡罗方法和克里格方法的综合运用,建立了描述渗透系数与污染物质释放强度之间关系的推算模型,进行了考虑参数不确定性的地下水污染源识别研究.研究结果表明,推算模型具有较高的精度,确定性系数和平均相对误差分别为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工具箱来模拟水流和污染物质运移的过程.

与实际问题不同,假想例子没有实际监测数据.所以将真实的污染源特征,输入污染物质运移模拟模型,正向运行模拟模型得到观测井的污染物质浓度,作为反向识别污染源特征时的实测数据.水流和污染物质在第1800和3600d的空间分布情况分别见图3和图4.观测井的污染物质浓度见图5.

图5 观测井的污染物质浓度实测数据

Fig.5 Measured data of contaminant concentration in observation wells

2.3 灵敏度分析

参数不确定性会影响污染物质运移模拟模型的输出结果.考虑模型所有参数对输出结果的影响,又会导致模型过于复杂.因此,应用局部灵敏度分析方法[21],筛选出对模拟模型输出结果影响最大的参数,作为输入模拟模型的随机变量,其它参数作为模拟模型的确定性变量.

筛选随机变量时,首先,将输入模拟模型的参数取均值(表3)输入模型,运行模型,得到观测井的污染物质浓度.然后,将输入模型的某个参数加减10%,20%,其它参数保持不变.再次运行模型,得到某一参数变化后,对应的观测井污染物质浓度.利用式(2)计算出各参数的灵敏度系数,选择灵敏度系数最大的参数,作为输入模型的随机变量.灵敏度分析结果见图6.

表3 参数服从的概率分布及取值范围

Table 3 Probability distribution and value range of parameters

由图6可以看出,灵敏度系数最大的参数是渗透系数.因此,将渗透系数作为输入模型的随机变量,其它参数作为确定性变量,取定值输入模型.根据前人经验[24],渗透系数服从的概率分布见表3.利用拉丁超立方抽样方法[22-23]对渗透系数在给定范围内抽样600组,每组参数样本具有5个值,分别对应研究区的5个参数分区.同时,利用该方法对污染物质释放强度抽样600组,每组样本具8个值,分别是2个污染源在4个释放时段的污染物质释放强度.使渗透系数和释放强度随机组合,得到600组由渗透系数和释放强度组成的输入样本.

图6 参数灵敏度分析

Fig.6 Parameter sensitivity analysis graph

应用模拟-优化方法识别地下水污染源时,污染物质运移模拟模型会作为等式约束连接到优化模型中,求解优化模型时,成百上千次的调用模拟模型,会产生大量的计算负荷,浪费大量的计算时间.建立模拟模型的替代模型,能有效解决这一问题.

2.4 建立替代模型

应用克里格方法建立污染物质运移模拟模型的替代模型.参考克里格方法的原理,利用MATLAB编写克里格替代模型的训练程序.然后,建立污染物质运移模拟模型的替代模型.建立替代模型的步骤如下:

1.将600组输入样本依次输入污染物质运移模拟模型,运行模拟模型,可以得到与600组输入样本一一对应的600组观测井的污染物质浓度.

2.利用540组输入样本作为替代模型的输入,污染物质的浓度作为替代模型的输出,训练替代模型.

3.将余下的60组输入样本,依次输入替代模型,得到替代模型输出的60组观测井的污染物质浓度.利用60组替代模型的输出和模拟模型的输出,检验替代模型的精度.

应用确定性系数和平均相对误差检验替代模型的精度,它们的计算方法见式14和15:

式中:yi是模拟模型输出的第i个时段末的污染物质浓度,mg/L;是替代模型输出的第i个时段末的污染物质浓度,mg/L;是模拟模型输出的各时段末的污染物质浓度均值,mg/L;n是观测样本的总数量.

2.5 建立优化模型

建立替代模型后,需要建立识别地下水污染源特征的非线性优化模型.非线性优化模型包括目标函数、决策变量和约束条件3部分.将观测井的污染物质实测浓度与模拟计算浓度拟合(替代模型的输出浓度)最小平方误差和,作为目标函数;将污染物质的释放强度,作为优化模型的决策变量;建立识别地下水污染源特征的非线性优化模型.优化模型的约束条件包括污染源释放强度不等式约束、5个分区渗透系数和地下水溶质运移规律等式约束;优化模型表示如下:

式中:Ki是某个参数分区的渗透系数,m/d;qm是污染物质释放强度,mg/d;是在观测井处污染物质的模拟计算浓度,mg/L;是在观测井处的污染物质实测浓度,mg/L.是替代模型.

利用遗传算法[25]求解某渗透系数组约束下的优化模型,只能得到与该组渗透系数影响下的污染源识别结果.而研究渗透系数不确定性对污染源识别结果的影响,需要统计多组渗透系数样本影响下的污染源识别结果.因此,应用蒙特卡罗方法进行考虑渗透系数不确定性的污染源识别.

2.6 建立推算模型

应用蒙特卡罗的思想,将优化模型作为“试验”发生器,将某渗透系数样本影响下的污染源识别结果,作为某次随机“试验”的结果.通过对成千上万组随机“试验”结果的统计与分析,考虑渗透系数不确定性对污染源识别结果的影响.

建立优化模型后,再次应用拉丁超立方方法,对渗透系数抽样90组.将第1组渗透系数,赋值到优化模型中,然后求解优化模型(整个求解过程中,渗透系数一直保持不变),得到与第1组渗透系数影响下的污染物质释放强度(包括S1的4个污染物质释放强度值和S2的4个污染物质释放强度值).对第2组渗透系数做与第1组渗透系数相同处理,直到第90组渗透系数.这样就依次求解了90个优化模型,获得了90组渗透系数影响下的污染物质释放强度.

但是,仅仅考虑90组渗透系数对污染源识别结果的影响,很难具有代表性.而求解成千上万个优化模型,又会产生海量的计算负荷,花费大量的计算时间.因此提出了应用克里格方法,建立能够描述渗透系数与污染物质释放强度之间关系的推算模型,来推算不同渗透系数取值影响下的污染物质释放强度.将70组渗透系数和与之影响下的污染物质释放强度,分别作为推算模型的输入与输出,训练推算模型.利用余下的20组渗透系数影响下的污染物质释放强度作为检验数据(每组包括S1的4个污染物质释放强度值和S2的4个污染物质释放强度值),检验推算模型的精度.推算模型与替代模型的精度评价指标相同.

3 结果与讨论

3.1 替代模型的精度评价

表4 替代模型的R2和MRE

Table 4R2and MRE of surrogate model

图7 7口观测井的污染物质浓度拟合

Fig.7 Fitting graph of contaminant concentration of seven observation wells

(a)~(g)分别表示第1到第7口观测井

对基于克里格方法替代模型的精度进行了评价.通常情况下,如果替代模型的R2高于0.9,MRE在10%以内,则认为替代模型精度满足研究需要.结合表4和图7,图8可以看出,替代模型与模拟模型输出浓度的平均相对误差均小于2%,确定性系数都大于0.99,接近1,替代模型与模拟模型的输出结果拟合程度非常好,替代模型的精度很高,可以代替模拟模型连接到优化模型中.

由图8可知,利用7口观测井浓度数据进行替代模型的精度检验时,出现了个别的异常值,但是异常值的相对误差均未超过4%,在10%以内,可见即使对于异常值,替代模型对模拟模型的拟合程度仍然很高.出现异常值的可能原因是,建立替代模型的训练数据,在个别的输入样本取值及取值领域内覆盖的不够全面,导致替代模型在该输入样本取值处泛化性能欠佳.因此,在获取训练数据时,尽量去覆盖输入样本的取值范围是很有必要的.

图8 7口观测井的污染物质浓度相对误差箱线图

Fig.8 The relative error box plot of the contaminant concentration of seven observation wells

3.2 识别结果的不确定性分析

图9 推算模型拟合曲线

Fig.9 Fitting graph of the inference model

应用克里格方法,建立了描述渗透系数与污染物质释放强度之间关系的推算模型.由图9可知, 推算模型的确定性系数达到了0.9895,与1非常接近;平均相对误差为4.51%,在10%以内;推算模型的精度很高,可以用来推算任意渗透系数影响下的污染物质释放强度.

考虑渗透系数不确定性的地下水污染源识别,需要统计成百上千组渗透系数影响下的污染源识别结果.具体做法是:再次对渗透系数抽样8000组,将其全部输入推算模型,可以得到各组渗透系数影响下的污染物质释放强度.对8000组污染物质释放强度进行统计与分析.分析因为渗透系数不确定性,导致的污染源识别结果的不确定性.

由图10可以看出,受渗透系数不确定性影响,污染源的识别结果也具有不确定性.应用切比雪夫不等式[22]估计不同置信水平对应的污染源识别结果置信区间.由表5可知,不同置信程度对应的置信区间,均包含2个污染源的实际释放强度.随着置信程度增大,置信区间变大,反之则反.决策者可以根据不同的管理目标,选择相信不同置信程度对应的置信区间,做为污染源特征的识别结果.也可以选择概率密度最大的释放强度作为污染源识别结果(识别结果见表6).

图10 不确定性分析

Fig.10 Uncertainty analysis diagram

(a)~(d)是S1各个释放时段的污染物质释放强度;(e)~(h)是S2各个释放时段的污染物质释放强度

表5 不同置信水平对应的污染源识别结果置信区间

Table 5 Confidence interval of contamination source identification results corresponding to different confidence levels

由表6可知,识别得到的S1和S2的释放强度与真实的污染源特征虽然有一定的误差,但是整体上对污染源的真实特征逼近程度较好.其中,S1第4时段的释放强度,与其真实释放强度差距最大,导致这样结果的原因可能有以下2方面:(1)替代模型和推算模型均具有误差,会在一定程度上影响识别结果的准确性.(2)抽样的过程具有随机性,抽样的样本不一定会完美覆盖到真实的含水层渗透系数.含水层渗透系数偏离真值,会导致污染源的识别结果偏离污染源特征的真实结果.

求解一个优化模型,进行1000次迭代计算,需要的时间约为0.45h,求解8000个优化模型需要3600h.而建立推算模型需要约40.5h,应用推算模型推算8000组渗透系数影响下的污染源特征,仅需要不到1s(瞬间输出).应用推算模型可以节省约99%的计算时间和计算负荷.

表6 概率密度最大的污染源识别结果

Table 6 Identification results of contamination sources with the highest probability density

4 结论

4.1 基于克里格方法建立的替代模型具有较高的精度,确定性系数高于0.99,平均相对误差小于1.2%.将克里格替代模型做为等式约束条件连接到优化模型中,供求解优化模型调用,在保证一定精度的前提下,可以减少大量的负荷和时间.

4.2 应用克里格方法,建立了描述渗透系数与污染物质释放强度之间关系的推算模型.建立的推算模型具有较高的精度,确定性系数和平均相对误差分别为0.9895和4.51%.运用克里格推算模型计算了8000渗透系数影响下的污染源特征,节省了约99%的计算负荷和时间.

4.3 对8000组污染源特征进行定量的统计与分析,统计得到了置信水平为80%、60%、40%和20%对应的污染源识别结果置信区间;分析出S1在4个释放时段概率密度最大的释放强度分别是74.22×105、53.67×105、40.35×105和16.36 ×105mg/d,S2在4个释放时段概率密度最大的释放强度分别是55.08× 105、50.50×105、39.41×105和21.01×105mg/d,将概率密度最大的释放强度做为污染源识别结果,对污染源的真实特征有较好的逼近程度.

4.4 将模拟-优化方法、灵敏度分析方法、蒙特卡罗方法和克里格方法结合,进行考虑渗透系数不确定性的污染源识别研究,更加符合实际情况,而且能够得到更加丰富的污染源识别结果,为决策者提供更多的参考依据.