張雙圣,強(qiáng) 靜,劉漢湖,劉喜坤,朱雪強(qiáng)

?

基于貝葉斯公式的地下水污染源識(shí)別

張雙圣1,3,強(qiáng) 靜2*,劉漢湖1,劉喜坤3,朱雪強(qiáng)1

(1.中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院,江蘇 徐州 221116；2.中國(guó)礦業(yè)大學(xué)數(shù)學(xué)學(xué)院,江蘇 徐州 221116；3.徐州市城區(qū)水資源管理處,江蘇 徐州 221018)

監(jiān)測(cè)井優(yōu)化；污染源識(shí)別；貝葉斯公式；信息熵；延遲拒絕自適應(yīng)Metropolis算法；拉丁超立方抽樣；多目標(biāo)優(yōu)化模型

由于地下水污染的發(fā)現(xiàn)具有滯后性,事件發(fā)生的時(shí)間、地點(diǎn)和污染源的強(qiáng)度等污染源信息往往是未知的,因此能否準(zhǔn)確地得到污染源的相關(guān)信息,對(duì)于地下水污染的治理具有重要的現(xiàn)實(shí)意義.

建立地下水水質(zhì)運(yùn)移模型,利用監(jiān)測(cè)井處的污染物濃度監(jiān)測(cè)值反求污染源位置、污染源釋放強(qiáng)度及釋放時(shí)間等信息,即為地下水污染源識(shí)別反問(wèn)題,其本質(zhì)是運(yùn)用監(jiān)測(cè)值對(duì)溶質(zhì)運(yùn)移模型參數(shù)進(jìn)行反演與識(shí)別.目前,求解此反問(wèn)題的方法主要包括貝葉斯統(tǒng)計(jì)方法[1]、地質(zhì)統(tǒng)計(jì)學(xué)方法[2]微分進(jìn)化算法[3]、遺傳算法[4]和模擬退火算法[5]等.這些算法可分為確定性方法和不確定性方法,其中確定性方法未考慮誤差因素,得到參數(shù)的一個(gè)確定的估計(jì)值,但容易丟掉“真值”;不確定性方法具有較強(qiáng)的隨機(jī)性,得到的是參數(shù)估計(jì)值的一個(gè)范圍,但保證了反演結(jié)果的可靠性.其中,貝葉斯統(tǒng)計(jì)方法以從監(jiān)測(cè)值中獲取參數(shù)信息為目標(biāo),將參數(shù)先驗(yàn)概率密度函數(shù)與樣本似然函數(shù)結(jié)合在一起,形成一套非常靈活、直觀、易于理解的統(tǒng)計(jì)方法,應(yīng)用越來(lái)越廣泛.

在對(duì)模型參數(shù)進(jìn)行反演識(shí)別時(shí),貝葉斯統(tǒng)計(jì)方法中經(jīng)常需要求解參數(shù)的后驗(yàn)估計(jì)值或者后驗(yàn)分布,在參數(shù)維數(shù)不是特別高時(shí)可以采用數(shù)值積分或者正態(tài)近似的方法求解[6].但是,隨著參數(shù)維數(shù)的增加,數(shù)值積分算法的計(jì)算量將呈指數(shù)增長(zhǎng),求解過(guò)程復(fù)雜而且難度較大,往往需要借助獨(dú)立抽樣的蒙特卡羅方法(MC方法)[7]進(jìn)行近似求解,其中馬爾科夫鏈蒙特卡羅方法(MCMC方法)[8-13]作為一種經(jīng)典抽樣方法得到廣泛應(yīng)用.近些年,比較流行的MCMC算法主要包括經(jīng)典Metropolis算法[8-9]、延遲拒絕算法(DR)[10],自適應(yīng)Metropolis算法(AM)[11-12],延遲拒絕自適應(yīng)Metropolis算法(DRAM)[13]等,但是這些算法均存在反演結(jié)果不收斂,或者局部最優(yōu)的問(wèn)題.

另一方面,受限于監(jiān)測(cè)經(jīng)費(fèi)及客觀條件,監(jiān)測(cè)方案[14-15](包括監(jiān)測(cè)井的位置、數(shù)量及監(jiān)測(cè)頻率)往往難以達(dá)到理想要求,造成監(jiān)測(cè)數(shù)據(jù)包含的信息量不充分或者監(jiān)測(cè)值和未知參數(shù)的關(guān)聯(lián)性較弱,且含有誤差,因此在對(duì)模型參數(shù)進(jìn)行反演識(shí)別時(shí),常常導(dǎo)致反問(wèn)題不適定性[16].為達(dá)到理想的模型參數(shù)反演效果,需對(duì)監(jiān)測(cè)井方案進(jìn)行優(yōu)化設(shè)計(jì).首先需要定義目標(biāo)函數(shù),對(duì)監(jiān)測(cè)井方案的信息量進(jìn)行量化[15].監(jiān)測(cè)方案優(yōu)化中常用的目標(biāo)函數(shù)有基于信息矩陣的A-最優(yōu),D-最優(yōu)和E-最優(yōu)等準(zhǔn)則[15],信噪比(SNR)[17],基于貝葉斯公式的相對(duì)熵[14,15,18]等.其中基于信息矩陣的A-最優(yōu)、D-最優(yōu)和E-最優(yōu)等準(zhǔn)則適用于線性模型及參數(shù)分布假設(shè)為正態(tài)分布情形;信噪比、相對(duì)熵這2個(gè)指標(biāo)適用于非線性模型及參數(shù)分布假設(shè)為非高斯情形[14,18].在現(xiàn)實(shí)案例中,地下水模型多為非線性模型,但是信噪比僅考慮監(jiān)測(cè)誤差對(duì)監(jiān)測(cè)數(shù)據(jù)的干擾影響,相對(duì)熵未考慮參數(shù)先驗(yàn)分布對(duì)后驗(yàn)分布的影響. 為此引入信息熵概念[19],信息熵是信息不確定性的度量,不確定性越大,信息熵越大.

1 研究方法

1.1 貝葉斯公式

貝葉斯公式[21]如下:

在實(shí)測(cè)數(shù)據(jù)固定的條件下,(6)式是關(guān)于參數(shù)的函數(shù).通過(guò)積分求解參數(shù)的后驗(yàn)分布很難得出顯式表達(dá)式,本研究采用MCMC方法對(duì)(6)式進(jìn)行求解.

1.2 基于貝葉斯公式與信息熵的監(jiān)測(cè)井優(yōu)化設(shè)計(jì)

(12)

式(12) 的求解算法較為復(fù)雜,難以得出顯示表達(dá)式,本文運(yùn)用蒙特卡羅方法[14]近似求解.

首先利用式(8)對(duì)式(12)改寫(xiě)如下:

運(yùn)用蒙特卡羅方法求解(14)式:

1.3 基于拉丁超立方抽樣的改進(jìn)型多鏈延遲拒絕自適應(yīng)Metropolis算法

Metropolis算法是一種經(jīng)典的MCMC 方法,應(yīng)用較為廣泛,但在在運(yùn)用貝葉斯公式對(duì)模型參數(shù)進(jìn)行反演識(shí)別時(shí),Metropolis算法容易出現(xiàn)局部最優(yōu),或者難以收斂,而且抽樣效率低下的問(wèn)題,本研究提出一種基于拉丁超立方抽樣的改進(jìn)型多鏈延遲拒絕自適應(yīng)Metropolis算法.

1.3.1 DRAM算法 延遲拒絕自適應(yīng)Metropolis算法(DRAM)由Haario等[13]在2006年首次提出,是一種高效自適應(yīng)MCMC方法,其本質(zhì)是把DR算法和AM算法2者組合起來(lái),既保證了馬爾科夫鏈的局部自適應(yīng),又保證了鏈條的全局自適應(yīng)調(diào)整策略,因此DRAM算法的抽樣效率優(yōu)于DR算法和AM算法各自單獨(dú)使用的情況[13].DRAM算法具體步驟如下:

(1)非自適應(yīng)階段

③重復(fù)過(guò)程②,直至達(dá)到迭代次數(shù)0.

(2)自適應(yīng)階段

③重復(fù)過(guò)程②,直至達(dá)到迭代次數(shù).

1.3.2 基于拉丁超立方抽樣的改進(jìn)型多鏈DRAM算法 為防止反演結(jié)果局部最優(yōu),或者產(chǎn)生難以收斂的問(wèn)題,且保證樣本初始點(diǎn)的隨機(jī)性和均勻性,本研究采用基于拉丁超立方抽樣[22]方法優(yōu)化抽樣過(guò)程.拉丁超立方抽樣是一種多維分層隨機(jī)抽樣方法,具有良好的散布均勻性和代表性.假設(shè)要在維模型參數(shù)先驗(yàn)范圍[A,B](=1,2,…,)內(nèi)抽取組樣本,具體步驟如下:

(1)把維模型參數(shù)的個(gè)先驗(yàn)范圍[A,B] (=1,2,…,)都均分成個(gè)小區(qū)間,小區(qū)間可記為[A,B](=1,2,…,,=1,2,…,),總共產(chǎn)生′個(gè)小區(qū)間;

(4) 矩陣中每個(gè)列向量是1組樣本,共抽取得到組樣本.

基于拉丁超立方抽樣的改進(jìn)型多鏈DRAM算法具體步驟:

第1步:運(yùn)用拉丁超立方抽樣方法,在模型參數(shù)先驗(yàn)范圍內(nèi)隨機(jī)抽取組初始樣本.

第2步:分別以組初始樣本作為初始點(diǎn)采用DRAM算法進(jìn)行迭代運(yùn)算,得到條Markov Chains.

第3步:對(duì)上面條Markov Chains求平均值作為最終的結(jié)果.

此方法改進(jìn)之處在于盡量削弱樣本初始點(diǎn)對(duì)算法結(jié)果的影響,符合Monte Carlo模擬方法的思想.

1.3.3改進(jìn)型多鏈DRAM算法的收斂性判斷 本文采用Gelman-Rubin收斂診斷方法[23]對(duì)多鏈DRAM算法后50%采樣過(guò)程的收斂性進(jìn)行判斷.收斂性判斷指標(biāo)為:

2 算例應(yīng)用

2.1 算例概述

假設(shè)研究區(qū)域內(nèi)含水層為二維均質(zhì)各向同性含水層,建立坐標(biāo)系,在S范圍內(nèi)向無(wú)限平面的均勻穩(wěn)定流中注入某污染物,瞬時(shí)注入的質(zhì)量為,為水流方向(圖1),則此時(shí)的對(duì)流-彌散方程[24]為:

式中:C為監(jiān)測(cè)點(diǎn)在t時(shí)刻的污染物濃度, ;t為污染物排放后經(jīng)過(guò)的時(shí)間,d;Dx、Dy分別為縱向、橫向的彌散系數(shù),;u為含水層水流平均流速,.

方程的解析解可表示為:

式中:為污染物的排放強(qiáng)度,g;為含水層的厚度,m;、分別為距污染源的縱向、橫向距離,m;為含水層的有效孔隙率.

監(jiān)測(cè)井位置為.區(qū)域水文地質(zhì)參數(shù)如表1所示.

表1 研究區(qū)域水文地質(zhì)參數(shù)

2.2 監(jiān)測(cè)方案優(yōu)化

表2 研究區(qū)域待求參數(shù)取值范圍

由1.2節(jié)可知,監(jiān)測(cè)井位置優(yōu)化問(wèn)題可概化如下:

圖2 遺傳算法中每一代目標(biāo)函數(shù)的最小值和平均值

式中:表示第組真實(shí)參數(shù),表示參數(shù)的第個(gè)分量.

表3 11種監(jiān)測(cè)方案,和

表4 從參數(shù)的先驗(yàn)分布中得到20組真實(shí)參數(shù)

2.2.2 監(jiān)測(cè)方案多目標(biāo)優(yōu)化模型 在現(xiàn)實(shí)生活中,污染源反演問(wèn)題常常不僅需要優(yōu)化監(jiān)測(cè)方案以使反演誤差最小(即信息熵最小),同時(shí)為了盡快找到污染源還要求監(jiān)測(cè)方案耗時(shí)最少,需要在信息熵和監(jiān)測(cè)方案耗時(shí)之間進(jìn)行權(quán)衡.設(shè)定監(jiān)測(cè)次數(shù)仍為5次,以信息熵最小和監(jiān)測(cè)耗時(shí)最短建立多目標(biāo)優(yōu)化模型,數(shù)學(xué)公式如下:

信息熵最小化目標(biāo):

監(jiān)測(cè)耗時(shí)最短目標(biāo):

一般來(lái)說(shuō),多目標(biāo)優(yōu)化問(wèn)題的最優(yōu)解并不唯一,它的最優(yōu)解集由那些任一個(gè)目標(biāo)函數(shù)值的減少都必須以增加其他目標(biāo)函數(shù)值為代價(jià)的解組成的集合,稱之為Pareto域.本文采用帶精英策略的非支配排序遺傳算法(NSGA-Ⅱ)解此多目標(biāo)優(yōu)化問(wèn)題.NSGA-Ⅱ算法是由Deb等[20]在NSGA的基礎(chǔ)上提出的.

圖4 目標(biāo)函數(shù)的Pareto前沿

NSGA-Ⅱ算法中參數(shù)如下設(shè)定:種群大小為100,最大進(jìn)化代數(shù)為50,交叉系數(shù)為0.8,變異系數(shù)為0.2.運(yùn)用Matlab中g(shù)amultiobj函數(shù)求解此多目標(biāo)優(yōu)化問(wèn)題,得到目標(biāo)函數(shù)的Pareto前沿,見(jiàn)圖4.

2.3 基于優(yōu)化監(jiān)測(cè)方案的污染源識(shí)別

由表5可知,2種監(jiān)測(cè)方案下,4個(gè)參數(shù)反演中0,02個(gè)參數(shù)反演效果均較好,均值相對(duì)誤差都在4.0%內(nèi);而參數(shù),0反演均值相對(duì)誤差在10%~ 20%,主要原因在于和0先驗(yàn)分布范圍較0,0大,導(dǎo)致和0的后驗(yàn)分布信息熵較大,其不確定性相應(yīng)增大,進(jìn)而反演誤差增大,因此為降低參數(shù)的反演誤差,需根據(jù)調(diào)查盡量縮小參數(shù)的先驗(yàn)分布范圍.

另一方面,與基于單目標(biāo)優(yōu)化的監(jiān)測(cè)方案反演結(jié)果相比,基于多目標(biāo)優(yōu)化的監(jiān)測(cè)方案條件下,污染源參數(shù)的反演均值相對(duì)誤差雖分別增加了0.4%、0.2%、0.3%、2.9%,但監(jiān)測(cè)時(shí)間卻顯著縮短了55.6%,因此基于多目標(biāo)優(yōu)化的監(jiān)測(cè)方案對(duì)于污染源的快速識(shí)別更具有現(xiàn)實(shí)意義.

3 結(jié)論

3.1 信息熵是參數(shù)反演結(jié)果精度的有效量度,信息熵越小,參數(shù)反演結(jié)果精度越高.參數(shù)反演結(jié)果的誤差與其先驗(yàn)分布范圍的選取有直接關(guān)系,參數(shù)先驗(yàn)分布范圍越大,其后驗(yàn)分布的信息熵較大,導(dǎo)致反演誤差增大.為降低參數(shù)的反演誤差,需根據(jù)調(diào)查盡量縮小參數(shù)的先驗(yàn)分布范圍.

3.2 與單目標(biāo)優(yōu)化監(jiān)測(cè)方案相比,多目標(biāo)優(yōu)化監(jiān)測(cè)方案雖然增加了反演結(jié)果的誤差,但是卻能顯著縮短監(jiān)測(cè)時(shí)間,對(duì)于污染源的快速識(shí)別更具有現(xiàn)實(shí)意義.

3.3 基于拉丁超立方抽樣的DRAM算法在保證反演結(jié)果準(zhǔn)確性的前提下,可顯著提高反演效率,其可靠性和穩(wěn)定性均優(yōu)于經(jīng)典Metropolis算法.

3.4 本文設(shè)計(jì)的案例為污染物在含水層中遷移的理想模型,具有解析解,在進(jìn)行監(jiān)測(cè)井的優(yōu)化設(shè)計(jì)及污染源識(shí)別時(shí),可將解析解公式通過(guò)matlab編程的方式嵌入到相關(guān)算法中,實(shí)現(xiàn)解析公式的反復(fù)調(diào)用.但在現(xiàn)實(shí)生活中,地下含水層的結(jié)構(gòu)及排污過(guò)程較為復(fù)雜,地下水在含水介質(zhì)中的遷移運(yùn)動(dòng)是一個(gè)十分復(fù)雜的系統(tǒng),一般通過(guò)建立數(shù)學(xué)模型(沒(méi)有解析解),運(yùn)用商業(yè)軟件(如GMS、Visual MODFLOW等)進(jìn)行數(shù)值模擬求解.在運(yùn)用貝葉斯公式和信息熵進(jìn)行監(jiān)測(cè)井優(yōu)化設(shè)計(jì)及污染源識(shí)別時(shí),需多次調(diào)用數(shù)值模擬模型,計(jì)算量巨大.本文采用地下水污染理想模型,既驗(yàn)證了基于貝葉斯公式與信息熵的監(jiān)測(cè)井優(yōu)化設(shè)計(jì)方法的可行性,又強(qiáng)有力的減少了計(jì)算工作量.

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

[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] 李久輝,盧文喜,常振波,等.基于不確定性分析的地下水污染超標(biāo)風(fēng)險(xiǎn)預(yù)警 [J]. 中國(guó)環(huán)境科學(xué), 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] 張江江.地下水污染源解析的貝葉斯監(jiān)測(cè)設(shè)計(jì)與參數(shù)反演方法 [D]. 杭州:浙江大學(xué), 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] 戴英彪.基于拉丁超立方試驗(yàn)設(shè)計(jì)的事故再現(xiàn)結(jié)果不確定性分析[D]. 長(zhǎng)沙:長(zhǎng)沙理工大學(xué), 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] 陳崇希,李國(guó)敏.地下水溶質(zhì)運(yùn)移理論及模型 [M]. 武漢:中國(guó)地質(zhì)大學(xué)出版社, 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] 周圣武,李金玉,周長(zhǎng)新.概率論與數(shù)理統(tǒng)計(jì)(第二版) [M]. 北京:煤炭工業(yè)出版社, 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

國(guó)家自然科學(xué)基金資助項(xiàng)目(51774270);國(guó)家水體污染控制與治理科技重大專項(xiàng)基金資助項(xiàng)目(2015ZX07406005)

*責(zé)任作者, 講師, 57591340@qq.com

張雙圣(1983-),男,山東省昌邑市人,中國(guó)礦業(yè)大學(xué)博士研究生,主要從事地下水污染控制研究.發(fā)表論文20余篇