基于貝葉斯公式的地下水污染源及含水層參數同步反演

張雙圣,劉漢湖,強 靜,劉喜坤,朱雪強

基于貝葉斯公式的地下水污染源及含水層參數同步反演

張雙圣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 quanti?cation and Bayesian inversion of a regional groundwater ?ow 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