基于贝叶斯公式的地下水污染源识别

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



基于贝叶斯公式的地下水污染源识别

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

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

监测井优化；污染源识别；贝叶斯公式；信息熵；延迟拒绝自适应Metropolis算法；拉丁超立方抽样；多目标优化模型

由于地下水污染的发现具有滞后性,事件发生的时间、地点和污染源的强度等污染源信息往往是未知的,因此能否准确地得到污染源的相关信息,对于地下水污染的治理具有重要的现实意义.

建立地下水水质运移模型,利用监测井处的污染物浓度监测值反求污染源位置、污染源释放强度及释放时间等信息,即为地下水污染源识别反问题,其本质是运用监测值对溶质运移模型参数进行反演与识别.目前,求解此反问题的方法主要包括贝叶斯统计方法[1]、地质统计学方法[2]微分进化算法[3]、遗传算法[4]和模拟退火算法[5]等.这些算法可分为确定性方法和不确定性方法,其中确定性方法未考虑误差因素,得到参数的一个确定的估计值,但容易丢掉“真值”;不确定性方法具有较强的随机性,得到的是参数估计值的一个范围,但保证了反演结果的可靠性.其中,贝叶斯统计方法以从监测值中获取参数信息为目标,将参数先验概率密度函数与样本似然函数结合在一起,形成一套非常灵活、直观、易于理解的统计方法,应用越来越广泛.

在对模型参数进行反演识别时,贝叶斯统计方法中经常需要求解参数的后验估计值或者后验分布,在参数维数不是特别高时可以采用数值积分或者正态近似的方法求解[6].但是,随着参数维数的增加,数值积分算法的计算量将呈指数增长,求解过程复杂而且难度较大,往往需要借助独立抽样的蒙特卡罗方法(MC方法)[7]进行近似求解,其中马尔科夫链蒙特卡罗方法(MCMC方法)[8-13]作为一种经典抽样方法得到广泛应用.近些年,比较流行的MCMC算法主要包括经典Metropolis算法[8-9]、延迟拒绝算法(DR)[10],自适应Metropolis算法(AM)[11-12],延迟拒绝自适应Metropolis算法(DRAM)[13]等,但是这些算法均存在反演结果不收敛,或者局部最优的问题.

另一方面,受限于监测经费及客观条件,监测方案[14-15](包括监测井的位置、数量及监测频率)往往难以达到理想要求,造成监测数据包含的信息量不充分或者监测值和未知参数的关联性较弱,且含有误差,因此在对模型参数进行反演识别时,常常导致反问题不适定性[16].为达到理想的模型参数反演效果,需对监测井方案进行优化设计.首先需要定义目标函数,对监测井方案的信息量进行量化[15].监测方案优化中常用的目标函数有基于信息矩阵的A-最优,D-最优和E-最优等准则[15],信噪比(SNR)[17],基于贝叶斯公式的相对熵[14,15,18]等.其中基于信息矩阵的A-最优、D-最优和E-最优等准则适用于线性模型及参数分布假设为正态分布情形;信噪比、相对熵这2个指标适用于非线性模型及参数分布假设为非高斯情形[14,18].在现实案例中,地下水模型多为非线性模型,但是信噪比仅考虑监测误差对监测数据的干扰影响,相对熵未考虑参数先验分布对后验分布的影响. 为此引入信息熵概念[19],信息熵是信息不确定性的度量,不确定性越大,信息熵越大.

1 研究方法

1.1 贝叶斯公式

贝叶斯公式[21]如下:

在实测数据固定的条件下,(6)式是关于参数的函数.通过积分求解参数的后验分布很难得出显式表达式,本研究采用MCMC方法对(6)式进行求解.

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

(12)

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

首先利用式(8)对式(12)改写如下:

运用蒙特卡罗方法求解(14)式:

1.3 基于拉丁超立方抽样的改进型多链延迟拒绝自适应Metropolis算法

Metropolis算法是一种经典的MCMC 方法,应用较为广泛,但在在运用贝叶斯公式对模型参数进行反演识别时,Metropolis算法容易出现局部最优,或者难以收敛,而且抽样效率低下的问题,本研究提出一种基于拉丁超立方抽样的改进型多链延迟拒绝自适应Metropolis算法.

1.3.1 DRAM算法 延迟拒绝自适应Metropolis算法(DRAM)由Haario等[13]在2006年首次提出,是一种高效自适应MCMC方法,其本质是把DR算法和AM算法2者组合起来,既保证了马尔科夫链的局部自适应,又保证了链条的全局自适应调整策略,因此DRAM算法的抽样效率优于DR算法和AM算法各自单独使用的情况[13].DRAM算法具体步骤如下:

(1)非自适应阶段

③重复过程②,直至达到迭代次数0.

(2)自适应阶段

③重复过程②,直至达到迭代次数.

1.3.2 基于拉丁超立方抽样的改进型多链DRAM算法 为防止反演结果局部最优,或者产生难以收敛的问题,且保证样本初始点的随机性和均匀性,本研究采用基于拉丁超立方抽样[22]方法优化抽样过程.拉丁超立方抽样是一种多维分层随机抽样方法,具有良好的散布均匀性和代表性.假设要在维模型参数先验范围[A,B](=1,2,…,)内抽取组样本,具体步骤如下:

(1)把维模型参数的个先验范围[A,B] (=1,2,…,)都均分成个小区间,小区间可记为[A,B](=1,2,…,,=1,2,…,),总共产生´个小区间;

(4) 矩阵中每个列向量是1组样本,共抽取得到组样本.

基于拉丁超立方抽样的改进型多链DRAM算法具体步骤:

第1步:运用拉丁超立方抽样方法,在模型参数先验范围内随机抽取组初始样本.

第2步:分别以组初始样本作为初始点采用DRAM算法进行迭代运算,得到条Markov Chains.

第3步:对上面条Markov Chains求平均值作为最终的结果.

此方法改进之处在于尽量削弱样本初始点对算法结果的影响,符合Monte Carlo模拟方法的思想.

1.3.3改进型多链DRAM算法的收敛性判断 本文采用Gelman-Rubin收敛诊断方法[23]对多链DRAM算法后50%采样过程的收敛性进行判断.收敛性判断指标为:

2 算例应用

2.1 算例概述

假设研究区域内含水层为二维均质各向同性含水层,建立坐标系,在S范围内向无限平面的均匀稳定流中注入某污染物,瞬时注入的质量为,为水流方向(图1),则此时的对流-弥散方程[24]为:

式中:C为监测点在t时刻的污染物浓度, ;t为污染物排放后经过的时间,d;Dx、Dy分别为纵向、横向的弥散系数,;u为含水层水流平均流速,.

方程的解析解可表示为:

式中:为污染物的排放强度,g;为含水层的厚度,m;、分别为距污染源的纵向、横向距离,m;为含水层的有效孔隙率.

监测井位置为.区域水文地质参数如表1所示.

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

2.2 监测方案优化

表2 研究区域待求参数取值范围

由1.2节可知,监测井位置优化问题可概化如下:

图2 遗传算法中每一代目标函数的最小值和平均值

式中:表示第组真实参数,表示参数的第个分量.

表3 11种监测方案,和

表4 从参数的先验分布中得到20组真实参数

2.2.2 监测方案多目标优化模型 在现实生活中,污染源反演问题常常不仅需要优化监测方案以使反演误差最小(即信息熵最小),同时为了尽快找到污染源还要求监测方案耗时最少,需要在信息熵和监测方案耗时之间进行权衡.设定监测次数仍为5次,以信息熵最小和监测耗时最短建立多目标优化模型,数学公式如下:

信息熵最小化目标:

监测耗时最短目标:

一般来说,多目标优化问题的最优解并不唯一,它的最优解集由那些任一个目标函数值的减少都必须以增加其他目标函数值为代价的解组成的集合,称之为Pareto域.本文采用带精英策略的非支配排序遗传算法(NSGA-Ⅱ)解此多目标优化问题.NSGA-Ⅱ算法是由Deb等[20]在NSGA的基础上提出的.

图4 目标函数的Pareto前沿

NSGA-Ⅱ算法中参数如下设定:种群大小为100,最大进化代数为50,交叉系数为0.8,变异系数为0.2.运用Matlab中gamultiobj函数求解此多目标优化问题,得到目标函数的Pareto前沿,见图4.

2.3 基于优化监测方案的污染源识别

由表5可知,2种监测方案下,4个参数反演中0,02个参数反演效果均较好,均值相对误差都在4.0%内;而参数,0反演均值相对误差在10%~ 20%,主要原因在于和0先验分布范围较0,0大,导致和0的后验分布信息熵较大,其不确定性相应增大,进而反演误差增大,因此为降低参数的反演误差,需根据调查尽量缩小参数的先验分布范围.

另一方面,与基于单目标优化的监测方案反演结果相比,基于多目标优化的监测方案条件下,污染源参数的反演均值相对误差虽分别增加了0.4%、0.2%、0.3%、2.9%,但监测时间却显著缩短了55.6%,因此基于多目标优化的监测方案对于污染源的快速识别更具有现实意义.

3 结论

3.1 信息熵是参数反演结果精度的有效量度,信息熵越小,参数反演结果精度越高.参数反演结果的误差与其先验分布范围的选取有直接关系,参数先验分布范围越大,其后验分布的信息熵较大,导致反演误差增大.为降低参数的反演误差,需根据调查尽量缩小参数的先验分布范围.

3.2 与单目标优化监测方案相比,多目标优化监测方案虽然增加了反演结果的误差,但是却能显著缩短监测时间,对于污染源的快速识别更具有现实意义.

3.3 基于拉丁超立方抽样的DRAM算法在保证反演结果准确性的前提下,可显著提高反演效率,其可靠性和稳定性均优于经典Metropolis算法.

3.4 本文设计的案例为污染物在含水层中迁移的理想模型,具有解析解,在进行监测井的优化设计及污染源识别时,可将解析解公式通过matlab编程的方式嵌入到相关算法中,实现解析公式的反复调用.但在现实生活中,地下含水层的结构及排污过程较为复杂,地下水在含水介质中的迁移运动是一个十分复杂的系统,一般通过建立数学模型(没有解析解),运用商业软件(如GMS、Visual MODFLOW等)进行数值模拟求解.在运用贝叶斯公式和信息熵进行监测井优化设计及污染源识别时,需多次调用数值模拟模型,计算量巨大.本文采用地下水污染理想模型,既验证了基于贝叶斯公式与信息熵的监测井优化设计方法的可行性,又强有力的减少了计算工作量.

[1] 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.

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

[3] Ruzek B, Kvasnicka M. Differential evolution algorithm in the earthquake hypocenter location [J]. Pure and Applied Geophysics, 2001,158:667-693.

[4] 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.

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

[6] Tanner M A. Tools for statistical inference: methods for the expectation of posterior distribution and likelihood functions [M]. New York: Springer-Verlag, 1996:1-52.

[7] 李久辉,卢文喜,常振波,等.基于不确定性分析的地下水污染超标风险预警 [J]. 中国环境科学, 2017,37(6):2270-2277. Li J, Lu W, Chang Z, et al. Risk prediction of groundwater pollution based on uncertainty analysis [J]. China Environmental Science, 2017, 37(6):2270-2277.

[8] 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.

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

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

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

[12] Xu W, Jiang C, Yan L, et al. An adaptive metropolis-hastings optimization algorithm of Bayesian estimation in non-stationary flood frequency analysis [J]. Water Resources Management, 2018,32(4): 1343-1366.

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

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

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

[16] Carrera J, Neuman S P. Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, stability, and solution algorithms [J]. Water Resources Research, 1986,22(2):211- 227.

[17] 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(1).

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

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

[20] Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-Ⅱ[J]. IEEE Transactions on Evolutionary Computation, 2002,6(2):182-197.

[21] Berger J O. Recent development and applications of Bayesian analysis [D]. Proceedings 50th ISI, Book I. 1995.

[22] 戴英彪.基于拉丁超立方试验设计的事故再现结果不确定性分析[D]. 长沙:长沙理工大学, 2011. Dai Y. Uncertainty analysis of vehicle accident reconstruction results based on Latin hypercube sampling [D]. Changsha: Changsha University of Science and Technology, 2011.

[23] Gelman A G, Rubin D B. Inference from iterative simulation using multiple sequences [J]. Statistical Science, 1992,7:457-472.

[24] 陈崇希,李国敏.地下水溶质运移理论及模型 [M]. 武汉:中国地质大学出版社, 1996:46-49. Chen C, Li G. Theory and model of solute transport in groundwater [M]. Wuhan: China University of Geosciences Press, 1996:46-49.

[25] 周圣武,李金玉,周长新.概率论与数理统计(第二版) [M]. 北京:煤炭工业出版社, 2007:208-215. Zhou S, Li J, Zhou C. Probability theory and mathematical statistics (Second edition) [M]. Beijin: China Coal Industry Publishing House, 2007:208-215.

Identification of groundwater pollution sources based on Bayes’ theorem.

ZHANG Shuang-sheng1,3, QIANG Jing2*, LIU Han-hu1, 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(4)：1568~1578

monitoring well optimization；contamination source identification；Bayes’ Theorem; information entropy；delayed rejection adaptive Metropolis algorithm；Latin hypercube sampling；multi-objective optimization model

X523

A

1000-6923(2019)04-1568-11

2018-08-20

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

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

张双圣(1983-),男,山东省昌邑市人,中国矿业大学博士研究生,主要从事地下水污染控制研究.发表论文20余篇