基于贝叶斯公式的地下水污染源及含水层参数同步反演

张双圣,刘汉湖,强 静,刘喜坤,朱雪强

基于贝叶斯公式的地下水污染源及含水层参数同步反演

张双圣1,3,刘汉湖1,强 静2*,刘喜坤3,朱雪强1

(1.中国矿业大学环境与测绘学院,江苏 徐州 221116；2.中国矿业大学数学学院,江苏 徐州 221116；3.徐州市城区水资源管理处,江苏 徐州 221018)

监测方案优化；贝叶斯公式；信息熵；Kriging替代模型；差分进化自适应Metropolis算法；参数后验均值偏离率

地下水污染源及含水层参数反演是指通过建立地下水水流及溶质运移模型,运用已有的污染物浓度监测数据反求污染源位置、污染源释放强度、释放时间等污染源信息以及含水层参数等.其本质是运用监测数据对地下水水流及溶质运移模型的源汇项及控制方程参数进行反演.地下水污染源及含水层参数的反演识别方法主要包括解析法、模拟优化方法和不确定性分析方法等.解析法形式简单,而且计算效率较高,但该方法主要适用于水文地质条件比较单一,而且必须满足一定的假定条件的地下水系统.由于反演问题均可转化为优化问题,因此模拟优化方法成为地下水污染源识别及含水层参数反演的常用方法.求解模拟优化问题的算法主要有模拟退火算法[1]、遗传算法[2]、禁忌搜索算法[3]等,这些算法均属于启发式智能搜索优化算法,具有全局搜索性能,保证全局最优,有效克服传统优化算法不收敛或局部最优的缺点.但是模拟优化算法未考虑参数取值、测量误差等各种不确定性因素,得到参数的一个确定的估计值,容易丢掉真值.

近年来,不确定性分析方法在研究参数反演问题中应用较多,主要包括广义似然不确定性估计方法[4]、地质统计学方法[5]、卡尔曼滤波方法[6]以及贝叶斯统计方法[7]等.贝叶斯统计方法应用最为广泛,常借助独立抽样的蒙特卡罗方法(MC)[8]进行近似求解,其中马尔科夫链蒙特卡罗方法(MCMC)[9-17]作为一种经典抽样方法得到广泛应用.单链MCMC算法[9-14]容易存在反演不收敛,或者反演结果局部最优的问题[17].而多链MCMC算法[15-16]适用于参数维度高,有多个局部最优值点,搜索量大的参数空间,能够更好地解决马尔科夫链局部收敛的问题[15-16].差分进化Metropolis算法(DE- MC)[15]和差分进化自适应Metropolis算法(DREAM)[16]是常用多链MCMC算法,其中 DREAM算法是DE-MC算法的改进版本,采用自适应随机子空间采样技术及能够自适应调整的交叉概率,并且运用IQR统计方法[16]去除无用链,可显著提高搜索效率和求解结果的精度.

在污染源及含水层参数反演过程中,往往出现监测数据与反演参数关联性较弱的问题,因此对监测方案(包括监测井的位置、数量及监测频率等)进行优化设计成为参数反演的关键.通常做法是定义某个目标函数,对监测方案的信息含量进行量化.常用的目标函数有信噪比(SNR)[18],以及基于贝叶斯公式的相对熵[19-20].但信噪比仅考虑监测误差对监测数据的干扰影响,而忽略了监测数据对参数后验分布的影响;相对熵未考虑参数先验分布对后验分布的影响.为此将信息熵概念[21-22]引入到地下水监测方案优化设计中.

另外在污染源及含水层参数反演过程中均需反复调用地下水水流及溶质运移数值模拟模型,导致巨大的计算负荷.通常采用构建数值模拟模型的替代模型的方法减少计算量.常用构建替代模型方法有多项式回归法[23],径向基函数神经网络法[24], Kriging法[25-26]等.其中Kriging替代模型具有兼顾局部和全局的统计特性,适用于解决非线性程度较高的问题.

本文充分考虑模拟模型参数不确定性、替代模型误差不确定性以及测量误差的不确定性,运用地下水模拟软件GMS(Groundwater Modeling System)建立二维非均质各向同性含水层水流及溶质运移数值模型.在确定初始监测时刻、监测间隔时间及监测次数条件下,利用最优拉丁超立方抽样方法及Kriging法建立数值模拟模型的替代模型,采用基于信息熵最小累进加井的方式制定监测方案,并根据此优化监测方案监测数据,基于贝叶斯统计方法运用DREAM算法进行污染源及含水层参数的同步反演.

1 研究方法

1.1 贝叶斯统计方法

贝叶斯公式如下:

本研究采用MCMC (Markov Chain Monte Carlo)算法对(5)式的后验分布进行求解.

1.2 基于贝叶斯公式与信息熵的监测方案优化设计

式(11)的求解算法较为复杂,难以得出显示表达式,本文运用文献[20,22]中蒙特卡罗方法近似求解.

1.3 DREAM算法及收敛性判定

DREAM算法具体步骤可查阅文献[16].采用Gelman-Rubin收敛诊断方法[27]对DREAM算法后50%采样过程的收敛性进行判断,判断指标为:

如果马尔科夫链满足Gelman-Rubin收敛准则,则计算终止,否则继续进化平行序列.

1.4 Kriging替代模型和抽样方法

为了减少监测方案优化设计及参数反演过程中多次调用地下水溶质运移数值模拟模型产生的计算量,运用Kriging法[25-26]构造数值模拟模型的替代模型.

替代模型的建立与DREAM算法初始样本的选取均需采取一定的抽样方法抽取样本.在替代模型构建过程中,为保证替代模型能够捕捉到研究对象函数的变化趋势,保证抽取的样本均匀分布在先验分布的全空间,以最优拉丁超立方抽样方法[28]进行样本抽样.而DREAM算法是多链MCMC算法,可采用拉丁超立方抽样方法[25]在参数先验分布范围内相对均匀的抽取随机样本作为初始样本,降低了随机选取样本对反演结果产生的影响.

2 算例应用

2.1 模型建立及问题概述

研究区内共有10眼监测井,地下水水质背景值较好,含水层污染物的初始浓度为零,某日研究区下游发现污染物X,初步断定污染源S在上游的某一区域范围内,且在某时间段内以注水井的形式(200m3/d)持续恒定地向含水层中排放污染物.研究区平面示意图如图1.

图1 算例模型示意

表1 研究区域已知水文地质参数

以西南角为坐标原点建立坐标系,根据研究区水文地质条件,建立地下水水流数值模型:

建立的地下水水流及溶质运移模型运用GMS软件中的MODFLOW和MT3DMS模块进行计算求解.为保证每个网格中心对应一个潜在污染源(污染源样本)的位置,将研究区域剖分为100行250列的正方形有限差分网格,基本单元格边长为4m.

表2 10眼监测井的坐标位置

表3 待求参数的先验分布范围

2.2 Kriging替代模型的建立

表4 从待求参数的先验分布范围中得到400组训练输入样本

(2)分别对10眼监测井建立Kriging替代模型将表4中400组参数分别代入到GMS软件中,得到10眼监测井在[800d,1027d]内每隔30d的污染物浓度计算值,作为Kriging替代模型输出值.将400组输入值与输出值作为训练样本代入MATLAB软件中,运用DACE工具箱建立各备选监测井的Kriging替代模型.

表5 从待求参数的先验分布范围中得到30组检验输入样本

图2 Kriging替代模型与数值模拟模型输出结果对比

表6 10眼备选监测井Kriging替代模型的R2,RMSE及WRE

由表6数据可知, 10眼备选监测井Kriging替代模型的平均决定系数2=0.9874,表明Kriging替代模型模拟精度较高.

2.3 监测方案优化设计

首先优化单井监测方案,即求(11)式的最小值,问题可概化如下:

表7 从参数的先验分布中得到20组“参数真值”

常用的参数反演效果评价指标有相对均方根误差和相对误差均值[30],指标具体公式如下:

②参数后验均值估计MP与“参数真值”之间相对误差均值

其中:表示第组“参数真值”,表示参数的第个分量.

表8 单井监测方案,

由表8可知,6#备选监测井的信息熵最小,故把其作为多井监测方案选定的第1眼监测井.后续选取多井监测方案中的第2眼监测井,步骤为:将选定6#井和其余9眼备选监测井分别进行组合,得到9种组合形式,并计算这9种组合监测方案的信息熵.此问题可概化如下:

由表9可知,信息熵随着监测井数量的增加而减小,其中4眼井组合监测方案信息熵、5眼井组合监测方案信息熵以及10眼井组合监测方案信息熵之间差距很小,故6眼井及以上组合监测方案信息熵不再赘述.兼顾反演精度及监测成本,并满足每个参数分区至少有1眼监测井的要求,本算例最优监测方案为5眼井组合监测方案(6,5,1,2,8)(定义为方案1).为了验证所选监测方案参数反演效果,与10眼井组合监测方案(1,2,3,4,5,6,7,8,9,10)(定义为方案2)的参数反演结果进行对比分析.

表9 多井监测方案优化设计中各监测方案的E(MP)

2.4 污染源及含水层参数同步反演

以表7中第1组参数真值为例,采用2.3节选定的2种监测方案,运用DREAM算法(初始平行链数为11,反演稳定后平行链数记为)进行参数反演,其中每条马尔科夫链长度是50000.

由表10可知,方案1中11个参数后验均值偏离率的平均值为10.3%,方案2中11个参数后验均值偏离率的平均值为9.1%,两种监测方案下参数后验均值偏离率结果相近,即2种方案的反演效果类似.但是方案2用了10眼监测井,方案1仅用了5眼监测井,方案2所需监测成本是方案1所需监测成本的2倍.对于本文案例,基于监测成本及反演精度双重考虑,认为方案1为最佳的监测方案.

表10 2种监测方案下模型参数后验统计结果及收敛性判断指标

2.5 灵敏度分析

由表10可知,污染源位置参数S和S的后验分布范围较先验分布范围显著减小(其中S的后验分布范围仅为先验分布范围的25%,S的后验分布范围为先验分布范围的10%),而其他参数的后验分布范围未见明显减小,这主要是由于各参数对污染物浓度监测值的灵敏度不同造成的.

为避免局部灵敏度分析方法没有考虑参数之间相互作用对输出结果影响的缺陷,运用全局灵敏度分析方法(Sobol’法)[31]及2.2节Kriging替代模型,得到参数对10眼备选监测井10个监测数据的1阶灵敏度系数,见表11.参照参数灵敏度分级[31](表12),参数S和S对于各监测井监测数据来讲属于中等灵敏参数或灵敏参数,而其余参数对于各监测井监测数据来讲多为不灵敏参数.

本文算例采用理想地下水水流及溶质运移模型进行了监测方案优化设计及参数反演.但这些方法应用到实际案例中时,首先需要保证建立的水流模型符合实际问题,需对水流模型各水文地质参数进行反演识别,然后再进行监测方案优化设计及溶质运移模型参数反演.

表11 参数a的1阶灵敏度系数

表12 参数灵敏度分级

3 结论

3.1 运用最优拉丁超立方抽样方法和Kriging法建立数值模拟模型的替代模型,其精度较高,能以较小的计算量得到和地下水数值模拟模型相近的输入输出关系,显著降低监测方案优化设计与污染源及含水层参数识别过程中反复调用地下水溶质运移数值模拟模型产生的计算负荷.

3.2 参数后验均值偏离率有效避免了参数真值取值对反演结果评价的影响,适用性更广.参数反演结果的后验均值偏离率与监测方案信息熵呈现较好的正相关关系.参数后验均值偏离率越小,说明参数反演效果越好.信息熵是参数后验分布不确定性的有效量度,信息熵越小,参数后验范围越小,基于贝叶斯公式与信息熵的监测方案优化设计方法是确定地下水污染监测方案的有效方法.

3.3 多井监测方案优化设计过程中,不同数量监测井组合产生的备选监测方案数量较大,为寻求最优监测方案,需求解所有备选监测方案的信息熵,导致计算量较大.而基于信息熵最小的累进加井的多井监测方案优化设计方法,可有效降低优化设计过程中求解信息熵的计算量, 监测方案优化效率得到显著提高.

3.4 DREAM算法采用多条马尔科夫链并行计算策略及自适应随机子空间采样技术,并且能够自适应调整交叉概率,可大大提高后验概率分布的搜索效率,能够有效解决局部收敛问题,而且能够实现多参数同步反演,计算效率显著提高.

[1] Dougherty D E, Marryott R A. Optimal groundwater management: 1. Simulated annealing [J]. Water Resources Research, 1991,27(10): 2493-2508.

[2] Mahinthakumar G, Sayeed M. Hybrid genetic algorithm-local search methods for solving groundwater source identification inverse problems [J]. Journal of Water Resources Planning and Management, 2005,131(1):45-57.

[3] Glover F, Laguna M. Tabu Search [M]. Boston, MA: Kluwer Academic Publisher, 1997:125-151.

[4] Hassan A E, Bekhit H M, Chapman J B. Uncertainty assessment of a stochastic groundwater flow model using GLUE analysis [J]. Journal of Hydrology, 2008,362(1/2):89-109.

[5] Snodgrass M F, Kitanidis P K. A geostatistical approach to contaminant source identification [J]. Water Resources Research, 1997, 33(4):537-546.

[6] Rasmussen J, Madsen H, Jensen K H, et al. Data assimilation in integrated hydrological modeling using ensemble Kalman filtering: evaluating the effect of ensemble size and localization on filter performance [J]. Hydrology and Earth System Sciences, 2015,12(2): 2267-2304.

[7] Chen M, Izady A, Abdalla O A, et al. A surrogate-based sensitivity quantification and Bayesian inversion of a regional groundwater flow model [J]. Journal of Hydrology, 2018,557:826-837.

[8] Roberts C P, Casella G. Monte Carlo statistical methods (Second edition) [M]. Berlin: Springer, 2004:79-122.

[9] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of state calculations by fast computing machines [J]. The Journal of Chemical Physics, 1953,21(6):1087-1092.

[10] Hastings W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika, 1970,57(1):97-109.

[11] Haario H, Saksman E, Tamminen J. An adaptive Metropolis algorithm [J]. Bernoulli, 2001,7(2):223-242.

[12] Tierney L, Mira A. Some adaptive Monte Carlo methods for bayesian inference [J]. Statistics in Medicine, 1999,18:2507-2515.

[13] Mira A. Ordering and improving the performance of Monte CarloMarkov Chains [J]. Statistical Science, 2002,16:340-350.

[14] Haario H, Laine M, Mira A. DRAM: Efficient adaptive MCMC [J]. Statistics and Computing, 2006,16(4):339-354.

[15] Ter Braak C J F. A Markov chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces [J]. Statistics and Computing, 2006,16(3):239-249.

[16] Vrugt J A, Ter Braak C J F, Diks C G H, et al. Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling [J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2009,10(3):273- 290.

[17] Vrugt J A. Markov chain Monte Carlo simulation using the DREAM software package: Theory, consepts, and Matlab implementation [J]. Enviromental Modeling & Software, 2016,75:273-316.

[18] Czanner G, Sarma S V, Eden U T, et al. A Signal-to-Noise Ratio Estimator for Generalized Linear Model Systems [J]. Lecture Notes in Engineering & Computer Science, 2008,2171:1063-1069.

[19] Lindley D V. On a measure of the information provided by an experiment [J]. The Annals of Mathematical Statistics, 1956,27(4): 986-1005.

[20] Huan X, Marzouk, Y M. Simulation-based optimal Bayesian experimental design for nonlinear systems [J]. Journal of Computational Physics, 2013,232(1):288-317.

[21] Shannon C E. A mathematical theory of communication [J]. The Bell System Technical Journal, 1948,27(3):379-423.

[22] 张双圣,强 静,刘汉湖,等.基于贝叶斯公式的地下水污染源识别 [J]. 中国环境科学, 2019,39(4):1568-1578. Zhang S, Qiang J, Liu H, et al. Identification of groundwater pollution sources based on Bayes’ theorem [J]. China Environmental Science, 2019,39(3):1568-1578.

[23] Knill D L, Giunta A A, Baker C A, et al. Response surface models combining linear and Euler aerodynamics for supersonic transport design [J]. Journal of Aircraft, 1999,36(1):75-86.

[24] Li J, Chen Y, Pepper, D. Radial basis function method for 1-D and 2-D groundwater contaminant transport modeling [J]. Computational Mechanics, 2003,32(1):10-15.

[25] 安永凯,卢文喜,董海彪,等.基于克里格法的地下水流数值模拟模型的替代模型研究 [J]. 中国环境科学, 2014,34(4):1073-1079. An Y, Lu W, Dong H, et al. Surrogate model of numerical simulation model of groundwater based on Kriging [J]. China Environmental Science, 2014,34(4):1073-1079.

[26] Lophaven S N, Nielsen H B, Sondergaard J. Dace: A MATLAB Kriging toolbox [R]. Kongens Lyngby: Technical University of Denmark, Technical Report No. IMM-TR-2002-12.

[27] Grlman A G, Rubin D B. Inference from iterative simulation using multiplesequences [J]. Statistical Science, 1992,7:457-472.

[28] Hickernell, F.A. A generalized discrepancy and quadrature error bound [J]. Mathematics of Computation of the American Mathematical Society, 1998,67(221):299-322.

[29] 初海波.基于自适应替代模型的DNAPLs污染地下水修复优化设计及其不确定性分析 [D]. 长春:吉林大学, 2015.Chu H. The optimization design and uncertainty analysis of DNAPLs- contaminated groundwater remediation based on the adaptive surrogate model [D]. Changchun: Jilin University, 2015.

[30] 张江江.地下水污染源解析的贝叶斯监测设计与参数反演方法 [D]. 杭州:浙江大学, 2017. Zhang J. Bayesian monitoring design and parameter inversion for groundwater contaminant source identification [D]. Hangzhou: Zhejiang University, 2017.

[31] Lenhart L, Eckhardt K, Fohrer N, et al. Comparison of two different approaches of sensitivity analysis [J]. Physics and Chemistry of the Earth, 2002,27:645-654.

Synchronous inversion of groundwater pollution source and aquifer parameters based on Bayesian formula.

ZHANG Shuang-sheng1,3, LIU Han-hu1, QIANG Jing2*, LIU Xi-kun3, ZHU Xue-qiang1

(1.School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China；2.School of Mathematics, China University of Mining and Technology, Xuzhou 221116, China；3.Xuzhou City Water Resource Administrative Office, Xuzhou 221018, China)., 2019,39(7)：2902~2912

Aiming at the optimization of monitoring schemes in the process of the identification of pollution source and the inversion of aquifer parameters in the heterogeneous underground aquifer, this paper proposes an optimization method for the multi-well monitoring schemes based on Bayesian formula and progressive addition of wells with minimum information entropy. Firstly, the two-dimensional heterogeneous isotropic subsurface groundwater flow and solute transport models under hypothetical case were constructed, and the numerical simulation models were solved by GMS software. The surrogate model of the numerical simulation model was established by the optimal Latin hypercube sampling method and Kriging method. Then Taking the minimum information entropy of the parameter posterior distribution as the objective function, the optimization design of multi-well monitoring schemes was carried out by means of progressive addition of wells. Finally, the differential evolution adaptive Metropolis algorithm was used to inverse the pollution source and aquifer parameters synchronously according to the optimized monitoring scheme. The case study results showed that: The 5combination monitoring scheme (6, 5, 1, 2, 8) under the condition of taking into account the inversion accuracy and monitoring cost and ensuring that there was at least one monitoring well in each parameter section was the optimal monitoring scheme. Compared with the 10combined monitoring scheme (1, 2, 3, 4, 5, 6, 7, 8, 9, 10) with the smallest information entropy, the 11parameters posterior mean deviation rate increased by 1.2%, but the monitoring cost was reduced by 50%.

monitoring well optimization；bayesian formula；information entropy；kriging surrogate model；differential evolution adaptive metropolis algorithm；parameter posterior mean deviation rate

X523

A

1000-6923(2019)07-2902-11

张双圣(1983-),男,山东省昌邑市人,工程师,博士,主要从事地下水污染控制研究.发表论文20余篇.

2018-11-27

国家水体污染控制与治理科技重大专项基金资助项目(2015ZX07406005)

* 责任作者, 讲师, 57591340@qq.com