地下水污染監測井優化設計及污染源識別

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

摘? ?要：在地下水污染源識別過程中，針對監測井監測值信息量不充分或監測值與模型參數關聯性較弱的問題，提出一種基于貝葉斯公式與信息熵的監測井優化方法. 構建二維地下水溶質運移模型，并運用GMS軟件進行數值求解. 為減少監測井優化設計及污染源識別過程中反復調用數值模型的計算負荷，采用克里金法建立數值模型的替代模型. 以信息熵作為優化指標，篩選出不同監測類型的最優監測方案，并以監測成本和反演精度為參考因素，選定相應監測方案，最后運用差分進化自適應Metropolis算法進行污染源識別. 算例研究表明： 7口監測井的克里金替代模型的決定系數均大于0.98，可較好地替代原數值模型. 基于監測成本最小的方案1（3號單井），其信息熵為12.772;兼顧監測成本和反演精度的方案2（井（2，3）組合），其信息熵為9.723;基于反演精度較高的方案3（3井（2，3，5）組合），其信息熵為9.377.方案1到方案3參數后驗分布范圍及標準差均逐漸減小，驗證了信息熵是參數后驗分布不確定性的有效量度.

關鍵詞：監測井優化;污染源識別;貝葉斯方法;信息熵;最優拉丁超立方抽樣;差分進化自適應Metropolis算法;克里金

中圖分類號：X523? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻標志碼：A

Optimization Design of Groundwater Pollution Monitoring

Wells and Identification of Contamination Source

ZHANG Shuangsheng1，3，LIU Hanhu1，QIANG Jing2，LIU Xikun3，ZHU Xueqiang1

（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）

Abstract： In the process of identifying groundwater pollution sources，a monitoring well optimization method based on Bayesian formula and information entropy is proposed for the problem that the monitoring value of monitoring wells is insufficient or the correlation between monitoring values and model parameters is weak. The two-dimensional groundwater contaminant transport model was numerically solved by GMS software. To reduce the computational load of the numerical model repeatedly in the optimization design of the monitoring wells and the identification process of the pollution source， the Kriging method was used to establish the surrogate model of the numerical model. As an optimization index， the optimal monitoring schemes of different monitoring types were selected， and the monitoring cost and inversion accuracy were taken as reference factors for the corresponding monitoring schemes. Then， the differential evolution adaptive Metropolis algorithm was used to identify the pollution source. The case study results show that： The determination coefficient of the Kriging surrogate models of the 7 monitoring wells was greater than 0.98， which indicated that the Kriging surrogate models can well replace the original numerical model. The scheme 1（single well No. 3） based on the lowest monitoring cost has an information entropy of 12.772;The scheme 2 （the combination of well No.2 and No.3） taking the monitoring cost and inversion accuracy into account has an information entropy of 9.723;The scheme 3（the combination of well No.2，3 and 5） with higher inversion precision has an information entropy of 9.377. Both the posterior distribution ranges and the standard deviation of model parameters from scheme 1 to scheme 3 were gradually reduced， which verifies that the information entropy is an effective measure of the uncertainty of the posterior distribution of the parameters.

Key words： monitoring well optimization;contamination source identification;Bayesian approach;information entropy;optimal Latin hypercube sampling;differential evolution adaptive Metropolis algorithm;Kriging

地下水污染具有隱蔽性、發現滯后性及修復難度大、費用高的特點，能否準確地得到污染源的相關信息，對于地下水污染的治理具有重要的現實意義.地下水污染源識別反問題是指通過建立地下水溶質運移模型，利用監測井處的污染物濃度監測值反求污染源位置、污染源釋放強度及釋放時間等信息. 由于建設監測井成本高昂，且存在監測數據包含的信息量不充分或者監測值和未知參數的關聯性較弱的問題，需要對現有監測方案進行優化設計.通常以信噪比（Signal-noise Ratio，SNR）[1]，基于貝葉斯公式的相對熵[2-3]作為監測井方案信息量的量度指標.信噪比僅考慮監測誤差對監測數據的干擾影響，相對熵未考慮參數先驗分布對后驗分布的影響.為此引入信息熵概念[4]，信息熵是信息不確定性的度量，不確定性越大，信息熵越大.

地下水污染源識別的求解方法主要包括貝葉斯統計方法[5-6]、地質統計學方法[7]，微分進化算法[8]、遺傳算法[9]和模擬退火算法[10]等. 其中，貝葉斯統計方法應用較為廣泛，在運用該方法對模型參數進行反演識別時，經常需要求解參數的后驗估計值或者后驗分布，在參數維數不是特別高時可以采用數值積分或者正態近似的方法求解[11]. 但是，隨著參數維數的增加，數值積分算法的計算量將呈指數增長，

求解過程復雜而且難度較大，往往需要借助獨立抽樣的蒙特卡羅方法（Monte Carlo 方法，簡稱MC方法）[12]進行近似求解，其中馬爾科夫鏈蒙特卡羅方法（Markov chain Monte Carlo方法，簡稱MCMC方法）[13-18] 作為一種經典抽樣方法得到廣泛應用. 近些年，比較流行MCMC算法主要包括經典Metropolis算法[13-14]、延遲拒絕算法（delayed rejection ，DR）[15-16]，自適應Metropolis算法（AM）[17]，延遲拒絕自適應Metropolis算法（DRAM）[18]等. 但是這些算法均是單鏈MCMC算法，容易出現反演結果不收斂，或者局部最優的問題.多鏈MCMC算法適用于參數維度高，有多個局部最優值點，搜索量大的參數空間，能夠更好地解決Markov chain局部收斂的問題[19]. 常用多鏈MCMC算法有DE-MC算法（Differential Evolution Markov Chain）[20]，DREAM算法（Differential Evolution Adaptive Metropolis）[21]等. DREAM算法是DE-MC

算法的改進版本，相比于DE-MC算法，DREAM算法采用自適應隨機子空間采樣技術及能夠自適應調整的交叉概率，并且運用IQR統計方法去除無用鏈，這幾個方面提高了DREAM算法的搜索效率和解的精度[21].

此外在監測井優化設計及污染源識別過程中需要多次調用地下水溶質運移數值模型，計算代價非常高，而替代模型的應用能夠有效地減少計算量.常用構造替代模型方法有多項式回歸法[22]、徑向基函數法[23-24]、人工神經網絡法[25]、克里金法（Kriging）[26-28]等，其中Kriging法作為多項式回歸分析的一種改進方法，包含多項式和隨機過程兩部分，同時具有局部和全局的統計特性，是一種監督式學習算法. 而且在MATLAB軟件中可用專門的DACE工具箱[29]建立Kriging替代模型，方便實用，因此得到廣泛應用.

綜合上述研究進展，本文建立二維地下水溶質運移數值模型，在確定初始監測時刻，監測間隔時間及監測次數條件下利用最優拉丁超立方抽樣方法及Kriging法建立數值模型的替代模型，以信息熵為評價指標分別計算不同監測類型下的最優監測井方案，并篩選出兼顧監測成本與監測精度的監測方案，然后以篩選出的監測方案運用DREAM算法進行污染源識別，為地下水污染源識別及監測方案優化研究提供借鑒.

1? ?研究方法

1.1? ?貝葉斯公式

貝葉斯公式如下：

在實測數據d固定的條件下，式（5）是關于參數α的函數.通過積分求解參數α的后驗分布很難得出顯式表達式，本研究采用MCMC （Markov Chain Monte Carlo）算法[21]對式（5）的后驗分布進行求解.

1.2? ?基于貝葉斯公式與信息熵的監測井優化設計

監測方案Monitoring Proposal （MP）的優化設計主要包括對監測井的數量、位置以及監測頻率的優化.假設初始監測時間為t1（固定值），由監測方案MP所得監測值仍記為d，此時貝葉斯公式可以改寫成：

運用污染物濃度監測值d反演未知參數α，α后驗概率密度函數p（α|d，MP），α的后驗分布的信息熵類似定義如下：

式（10）左端項含有監測值d，但在進行監測井方案設計時，并未真正得到d，可認為監測值d是隨機變量，其概率密度函數為p（d|MP）. 為了得到一個僅包含監測方案MP的函數，對式（10）兩側乘以

p（d|MP），再對d積分，求出信息熵H（MP，d）的期

望為：

式（11）中E（H（MP，d）只受到監測方案MP的影響，可簡記為E（MP）. 通過求解E（MP）的最小值，就可以獲得最優監測井設計方案MP*.根據信息熵概念，利用由MP*獲得的污染物濃度監測值d*反演未知參數α，此時α后驗分布的信息熵最小，表示 α的不確定性也最小，反演效果最優.

因此，只要監測設計方案MP固定，通過式（12）（14）（15）就可以得到此種監測方案的信息熵近似值.

1.3? ?DREAM算法

1.3.1? ?DREAM算法具體步驟

DREAM算法[21]是在DE-MC算法[20]的基礎上提出的. DREAM算法步驟如下：

1）在模型未知參數α先驗范圍內隨機產生NP

個初始樣本Xi（t）=（Xi，1（t），Xi，2（t），…，Xi，m（t））T（i=1，2，…，NP，t=0），作為NP條馬爾科夫鏈的起始點.

2）針對第i（i = 1，2，…，NP）條馬爾科夫鏈，運用DE方法產生參數的變異樣本Zi（t）.

式中： I表示m階單位矩陣;e表示m階方陣，其對角線上元素服從均勻分布U（-b，b），b為自定義的極小值;δ表示用于產生候選樣本的平行鏈數的二分之一，r1（j），r2（k）∈{1，2，…，NP}為隨機選取的平行鏈編號，且當j，k = 1，2，…，δ 時，r1（j）≠r2（k）≠i;

γ（δ，d′）為比例因子;ε 表示m維向量，其m個元素均服從正態分布N（0，b*），b*為自定義的極小值.

3）引入交叉概率Cr，交叉混合Xi（t）和Zi（t）得到候選樣本Vi如下：

對？坌i∈{1，2，…，NP}，當j = 1時，令d′ = d;對j = 1，2，…，m，有

Xi，j（t） = Xi，j（t）? ? if? U≤1 - Cr，d′ = d′ - 1Zi，j（t）? ? 其他 （17）

其中0≤Cr≤1，U是區間[0，1]上的隨機數.

5）根據IQR統計方法[21]去除無用鏈.

6）判斷收斂性.如果馬爾科夫鏈滿足Gelman-Rubin收斂準則[30]，若滿足條件則計算終止，否則繼續進化平行序列.

1.3.2? ?DREAM算法的收斂性判斷

本文采用Gelman-Rubin收斂診斷方法[30]對

DREAM算法后50%采樣過程的收斂性進行判斷，判斷指標為：

1.4? ?最優拉丁超立方抽樣方法和Kriging替代模型

1.4.1? ?最優拉丁超立方抽樣方法

拉丁超立方抽樣是一種多維分層隨機抽樣方法，但抽樣結果有很大地隨機性，每一次抽樣得到的結果差別很大.假設要在m維參數α取值范圍[0，1]m內運用拉丁超立方抽樣方法抽取q組樣本，q組樣本數據用矩陣Φqm表示;符合要求的矩陣Φqm共有（q！）m種，總共可以得到（q！）m種拉丁超立方抽樣方案，（q！）m種抽樣方案Φqm組成的集合記為？撰.這些抽樣方案對整個抽樣空間的覆蓋填充程度有很大地差異[28]. 本文運用中心化L2偏差[31]作為優化指標尋求覆蓋填充程度最優的抽樣方案，中心化L2偏差表達式如下：

1.4.2? ?Kriging替代模型

Kriging法[26]是由Matheron等人發明的一種優化插值方法，其具體工作原理可查閱文獻[27]，基本步驟如下：

Kriging替代模型響應值y與自變量k之間的關系可用下式表示：y（k） = f T B + Z，式中k∈IRmk （mk表示自變量k維數），f（k）=（f1（k），f2（k），…，fp（k））T，其中fi（·）（i = 1，2，…，p）表示事先確定的多項式函數，常常是0階，1階或者2階多項式;β = （β1，β2，…，βp）T為待定參數;Z為誤差隨機變量，E（Z） = 0，Var（E） = σ2.

2002年，Lophaven等人[29]基于Matlab平臺創建了DACE工具箱，用于生成Kriging替代模型，本文中Kriging替代模型的建立就是基于此工具箱.

2? ?算例應用

2.1? ?模型建立及問題概述

2.1.1? ?模型建立

假定研究區域為矩形區域，長1 000 m，寬600 m，含水層為厚度35 m的砂質潛水含水層（水文地質參數見表1），西部邊界Γ1與東部邊界Γ3為給定水頭邊界，其中東部水頭為25 m，西部水頭為30 m，北部邊界Γ2與南部邊界Γ4為隔水邊界.天然狀態下地下水流為自西向東流動的二維均質各向同性的非承壓穩定流.

研究區內共有7眼監測井，地下水水質背景值較好，含水層污染物的初始濃度為零，某日研究區下游發現污染物X，初步斷定污染源S在上游的某一區域范圍（先驗范圍）內，且在某時間段內以注水井的形式（200 m3/d）持續恒定地向含水層中排放污染物. 研究區平面示意圖如圖1所示.

以西南角為坐標原點建立坐標系，根據研究區水文地質條件，建立地下水水流數值模型：

式中：Dx、Dy分別為水動力彌散系數在x、y方向上的分量，m2/d;vx、vy分別為x、y方向上地下水的滲流速度，m/d;n為含水層介質的孔隙度，無量綱;c為污染物質量濃度，mg/L;Qinj為向含水層的注入液體量，m3/d;Cinj為隨液體進入含水層的污染物質量濃度，mg/L;f（x，y，t）表示僅在水動力彌散作用下，單位時間通過實際過水斷面單位面積的溶質質量.

建立的地下水水流及溶質運移模型運用GMS軟件中的MODFLOW和MT3DMS模塊進行計算求解. 為保證每個網格中心對應一個潛在污染源（污染源樣本）的位置，將研究區域剖分為150行200列的正方形有限差分網格，基本單元格邊長為4 m.

2.1.2? ?問題概述

針對潛在的污染源范圍，要求在現有7眼監測井的基礎上，制定優化的監測井方案，以此進行污染源的識別（包括污染源的位置，污染物投放和結束時間，以及污染物的投放質量濃度），即求解污染源未知參數α = （Xs，Ys，T1，T2，Cs），其中（Xs，Ys）為污染源位置，m;T1、T2分別為污染源開始排放和結束排放的時刻，d;Cs是注入污染物的質量濃度，mg/L.表2為7眼監測井的坐標位置.

以未發生污染某確定時刻作為初始時間，此時t=0.要求從第750 d至第900 d內完成污染源識別任務（期間每10 d進行一次監測，每眼井共監測16次），假設上述5個參數的先驗分布為均勻分布α = （Xs，Ys，T1，T2，Cs）5個待求參數的先驗范圍分別

如下：

60 m≤Xs≤220 m，240 m≤Ys≤400 m，8 d≤T1≤12 d，18 d≤T2≤22 d，2 800 mg/L≤Cs≤3 300 mg/L.

2.2? ?Kriging替代模型的建立

1）運用最優拉丁超立方抽樣方法從污染源未

知參數α的先驗范圍內抽得均勻分布的100組參數作為Kriging替代模型訓練輸入值（此時CL2（Φ100，5）=0.002 462），如表3所示.

2）分別對7眼監測井建立Kriging替代模型. 將表3中100組參數分別代入GMS軟件中，得到7眼監測井不同監測時間點的污染物濃度值，作為

Kriging替代模型輸出值.將100組輸入值與輸出值作為訓練樣本代入MATLAB軟件，利用DACE工具箱對各監測井的Kriging替代模型進行訓練.

在未知參數的先驗范圍內運用拉丁超立方抽樣方法重新得到10組參數作為檢驗樣本的輸入值（表4），再次代入GMS軟件中，得到7眼監測井不

以7號監測井為例，以檢驗樣本的輸出值作為橫坐標，替代模型的輸出值作為縱坐標，繪制替代模型輸出值與檢驗樣本輸出值的對比圖（圖2）. 由圖2可知，數據散點均集中分布于y = x直線上，表明替代模型可較好地對數值模型進行替代.

為進一步檢驗替代模型的精度，分別運用決定系數，平均絕對誤差及均方根誤差3個指標對7眼監測井的替代模型進行檢驗評價，結果如表5所示.

由表5數據可知，Kriging替代模型預測精度較高，且7眼備選監測井平均絕對誤差的平均值為0.14，平均決定系數R2 = 0.996 8. 在參數反演及監測井優化問題中，由于運用Kriging替代模型，整體誤差包括替代模型誤差及測量誤差. 假設第i（i=1，2，…，7）眼備選監測井的Kriging替代模型誤差εi′滿足正態分布N（0，（σi′）2），且均值E（εi′）=0，均方差σi′ =RMSEi;假設測量誤差ε″滿足正態分布N（0，（σ″）2），且均值E（ε″） = 0，均方差σ″ = 0.01. 由于替代模型誤差εi′及測量誤差ε″相互獨立，故第i（i=1，2，…，7）口監測井整體誤差εi = εi′ + ε″滿足正態分布N（0，（σi′）2+（σ″）2），據此可以確定1.1節中式（4）和式（5）中協方差矩陣C（ε）.

2.3? ?監測方案優化設計

7眼井可得到7組監測數據，任取i（i=1，2，…，7）組監測數據作為參數反演的監測數據d，共有Ci7（i=1，2，…，7）種組合形式，每種組合形式代表一種監測方案，從而監測方案共有Ci7（i=1，2，…，7）種. 由于i=1，2，…，7，按照選取的監測井數量進行劃分，監測類型共有7類. 根據1.2節可知，每類監測方案優化問題其實就是在Ci7（i=1，2，…，7）種監測方案里面選取信息熵最小的監測方案，信息熵最小的監測方案可視為最優監測方案. 問題可概化如下：

min E（MP） = -ln p（α） + min U（MP）? （19）

式中：MP表示監測方案.

根據式（12）（14）（15）分別求得不同監測方案的信息熵E（MP）.為了驗證基于貝葉斯公式與信息熵的監測井優化設計效果，以信息熵E（MP）與反演結果相對誤差均值MRE（MP）為指標對Ci7種監測方案進行評價.從參數α的先驗分布里運用拉丁超立方抽樣方法隨機并且均勻地得到20組參數作為真實參數（表6），對應于Ci7種監測方案，20組真實參數通過Kriging替代模型產生了20×Ci7組濃度監測值.利用產生的監測值，運用DREAM算法（初始平行鏈數為10）反演參數α，其中每條馬爾科夫鏈長度為12 000，在馬爾科夫鏈長度10 000時所有參數的收斂性判斷指標■<1.2.為了保證精度，只將馬爾科夫鏈趨于穩定后的最后2 000組樣本進行后驗統計，得出參數后驗均值估計MMP . 并將表6中的真實參數α一并代入MRE（MP）表達式，如下：

式中： j為第j組真實參數α;k表示參數α的第k個分量. 由式（20）得出Ci7種監測方案的反演結果MRE（MP），如表7所示.

分別將表7中1眼井監測方案，2眼井組合監測方案以及3眼井組合監測方案下的MRE（MP）與 E（MP）進行線性擬合，結果見圖3和表8. 由于4眼井組合監測方案，5眼井組合監測方案以及6眼井組合監測方案下E（MP）數值接近，同時MRE（MP）數值也十分接近，并且考慮到計算誤差的影響，不再對這3種類型監測方案下MRE（MP） 與E（MP）進行線性擬合.

由于圖3（c）中決定系數R2=0.552比較小，運用 F檢驗對圖3（c）中線性方程的顯著性進行檢驗[32].

檢驗假設H0 ： E（MP） 和RT（MP）之間沒有真正的線性關系，H1 ： E（MP）和RT（MP）之間有線性關系;顯著性水平α = 0.05.

經計算可得F檢驗統計量F = 40.739，而Fα（1，35-2） = Fα（1，33） = 4.139，因此F > Fα（1，33），拒絕H0 . F檢驗表明E（MP） 和RT（MP）之間存在線性關系. 另外圖3（c）中實線表示擬合直線，兩側虛線標出其95%置信區間，只有2個點落在置信區間外面，也充分說明了E（MP） 和RT（MP）之間具有正線性關系.

綜上，由圖3和表8可知，MRE（MP）與E（MP）呈較好的正相關關系，說明E（MP）是參數反演結果精度的有效量度，E（MP）越小，參數反演結果精度越高. 但表7中E（MP）取得最小值的監測方案的MRE（MP）未必是最小值（如1眼監測井情況下，選取3監測方案時E（MP）為12.772，MRE（MP）為0.069;選取2監測方案時E（MP）為13.004，MRE（MP）為0.063），主要原因在于選取20組參數真值（表6）時，雖然運用拉丁超立方抽樣方法使20組參數真值盡量均勻分布在參數先驗范圍內，但由于數據較少，無法在真正意義上均勻分布在先驗范圍內，而信息熵的求解運用蒙特卡洛方法（MC方法）在參數先驗范圍內運用拉丁超立方抽樣方法抽取樣本20 000次，二者相比，因此認為以E（MP）最小值作為選取最優監測方案的指標更加可信，而不能以MRE（MP）最小值作為選取最優監測方案的指標.

另由表7可知，并不是任何條件下監測井數量越多，信息熵越小，如兩眼監測井情況下，（2，3）組合方式下的E（MP）為9.723，3眼監測井情況下，（1，5，6）組合方式下的E（MP）為11.666. 表7顯示每種監測類型下均存在信息熵最小的最優監測井方案.因此從表7中各監測類型中篩選出信息熵最小的監測方案作為相應監測井組合方案中的最優方案，并繪制不同監測類型中最優方案的信息熵

E（MP）隨監測井數量的變化曲線，如圖4所示.

由圖4可知，在各監測類型的最優方案條件下，信息熵隨著監測井數量的增加而減小，其中兩眼井的信息熵顯著小于一眼井的信息熵，但是與其他數量的監測井的信息熵相差不大.

通常情況下，對于利用監測井進行污染源識別，既要考慮反演精度，又要考慮監測成本. 因此本算例選取3種監測方案進行污染源識別. 方案1：監

測成本最小的1眼監測井方案（3號單井）. 方案2：兼顧反演精度與監測成本的監測井方案：由于2眼井的信息熵與3～7眼井的信息熵相差不大，但隨著監測井數量的增加，監測費用顯著增加，因此該方案選2眼井監測方案（監測井（2，3）組合）. 方案3：基于較高反演精度的監測井方案，由于3～7眼井的信息熵相差不大，認為反演精度相近，因此該方案選3眼井監測方案（監測井（2，3，5）組合）.

2.4? ?基于優化監測方案的污染源識別

以表6中第1組參數真值α=（Xs，Ys，T1，T2，Cs）=（81.4，266.3，11.2，19.9，3 011.6）為例，分別利用2.3節選取的3種監測方案，運用DREAM算法（初始平行鏈數為10，反演穩定后平行鏈數記為n）進行污染源參數反演，其中每條馬爾科夫鏈長度是12 000.

由表9和圖5可知，方案1到方案3參數后驗分布范圍逐漸減小，與相應監測方案信息熵的變化趨勢一致（圖4），且參數真值均在參數后驗分布范圍內. 進一步表明監測方案信息熵越小，參數后驗分布不確定性越小.

為進一步說明3種監測方案的參數反演效果，對馬爾科夫鏈穩定后的剩余2 000個樣本點進行后驗統計分析，結果見表10.

由表10可知，從方案1到方案3，5個參數的標準差均逐漸減小;方案2與方案3下5個參數后驗均值相對誤差數值接近，與方案1相比，Xs，Ys以及T1的后驗均值相對誤差均大幅減少，但是T2、Cs 2個參數的后驗均值相對誤差卻有所增大.結合圖3、圖4及表10分析，進一步說明參數后驗分布的信息熵與5個參數后驗均值相對誤差的平均值成正比，但不能保證參數后驗分布的信息熵與每個參數后驗均值相對誤差成正比，即運用參數后驗分布的信息熵作為監測井優化設計的指標，在參數整體反演精度最高的情況下，不能保證每個參數的反演效果均達到最佳.

監測井數越多所需費用越多，單井監測費用最少. 方案2與方案3信息熵差距很小，5個參數后驗樣本均值相對誤差的平均值也近似相等，卻增加了1眼井進行監測，費用明顯增加. 對于本文案例，如果需要綜合考慮監測成本及參數后驗分布范圍大小，認為方案2為最佳的監測井方案.

3? ?結? ?論

1）與拉丁超立方抽樣方法相比，最優拉丁超立方抽樣可有效提高參數先驗分布抽樣的均勻性，避免抽樣結果的隨機性及對整個抽樣空間的覆蓋填充程度的差異性.

2）Kriging替代模型屬于黑箱模型，能以較小的計算量得到和地下水數值模型相近的輸入輸出關系，在保證模擬精度的條件下，顯著降低計算負荷.

3）參數反演結果的相對均方根誤差與監測方案信息熵呈現較好的正相關關系. 信息熵是參數后驗分布不確定性的有效量度，信息熵越小，參數后驗范圍越小，基于貝葉斯公式與信息熵的監測方案優化設計方法是確定地下水污染監測方案的有效方法.

4）并非監測井數量越多越有利于污染源的反

演識別，必須以信息熵為篩選指標制定各監測類型下的最優監測井組合方案，而且可以以監測成本、監測效率、反演精度等為限制條件進行具體工況條件下最優監測方案的選擇.

參考文獻

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

[2]? ? HUAN X，MARZOUK Y M. Simulation-based optimal Bayesian experimental design for nonlinear systems [J]. Journal of Computational Physics，2013，232（1）：288—317.

[3]? ? LINDLEY D V. On a measure of the information provided by an experiment [J]. The Annals of Mathematical Statistics，1956，27（4）：986—1005.

[4]? ? SHANNON C E. A mathematical theory of communication [J]. The Bell System Technical Journal，1948，27（3）：379—423.

[5]? ? SOHN M D，SMALL M J，PANTAZIDOU M. Reducing uncertainty in site characterization using bayes monte carlo methods [J]. Journal of Environmental Engineering-asce，2000，126（10）：893—902.

[6]? ? 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.

[7]? ? SNODGRASS M F，KITANIDIS P K. A geostatistical approach to contaminant source identification [J]. Water Resources Research，1997，33（4）：537—546.

[8]? ? RUZEK B，KVASNICKA M. Differential evolution algorithm in the earthquake hypocenter location [J]. Pure and Applied Geophysics，2001，158：667—693.

[9]? GIACOBBO F，MARSEGUERRA M，ZIO E. Solving the inverse problem of parameter estimation by genetic algorithms：the case of a groundwater contaminant transport model [J]. Annals of Nuclear Energy，2002，29（8）：967—981.

[10] DOUGHERTY D E，MARRYOTT R A. Optimal groundwater management：simulated annealing [J]. Water Resources Research，1991，27（10）：2493—2508.

[11]? TANNER M A. Tools for statistical inference：methods for the expectation of posterior distribution and likelihood functions [M]. Berlin：Springer，1996：14—39.

[12] ROBERTS C P，CASELLA G. Monte carlo statistical methods [M] .2nd ed. Berlin：Springer，2004：79—122.

[13]? 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.

[14] HASTINGS W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika，1970，57（1）：97—109.

[15]? TIERNEY L，MIRA A. Some adaptive Monte Carlo methods for bayesian inference [J]. Statistics in Medicine，1999，18：2507—2515.

[16]? MIRA A. Ordering and improving the performance of Monte Carlo Markov Chains [J]. Statistical Science，2002，16：340—350.

[17]? HAARIO H，SAKSMAN E，TAMMINEN J. An adaptive metropolis algorithm [J]. Bernoulli，2001，7（2）：223—242.

[18]? HAARIO H，LAINE M，MIRA A. DRAM：efficient adaptive MCMC [J]. Statistics and Computing，2006，16（4）：339—354.

[19]? 張江江. 地下水污染源解析的貝葉斯監測設計與參數反演方法 [D]. 杭州：浙江大學環境與資源學院，2017.

ZHANG J J. Bayesian monitoring design and parameter inversion for groundwater contaminant source identification [D]. Hangzhou：College of Environmental and Resource Sciences，Zhejiang University，2017. （In Chinese）

[20]? 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.

[21]? 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.

[22]? 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.

[23]? 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.

[24]? 肖傳寧，盧文喜，趙瑩，等.基于徑向基函數模型的優化方法在地下水污染源識別中的應用 [J]. 中國環境科學，2016，36（7）：2067—2072.

XIAO C N，LU W X，ZHAO Y，et al. Optimization method of identification of groundwater pollution sources based on radial basis function model [J]. China Environmental Science，2016，36（7）：2067—2072. （In Chinese）

[25] CHRISTIE M，DEMYANOV V，ERBAS D. Uncertainty quantification for porous media flows [J]. Journal of Computational Physics，2006，217（1）：143—158.

[26] MATHERTON G. Principles of geostatistics [J]. Economic Geology，1963，58（8）：1246—1266.

[27]? SACKS J，WELCH W J，MITCHELL T J，et al. Design and analysis of computer experiments [J]. Statistical Science，1989，4（4）：409—435.

[28]? 高月華.基于Kriging 代理模型的優化設計方法及其在注塑成型中的作用 [D]. 大連：大連理工大學運載工程與力學學部，2008.

GAO Y H. Optimization methods based on Kriging surrogate model and their application in injection molding [D]. Dalian：Faculty of Vehicle Engineering and Mechanics，Dalian University of Technology，2008. （In Chinese）

[29]? 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.

[30] GELMAN A G，RUBIN D B. Inference from iterative simulation using multiple sequences [J]. Statistical Science，1992，7：457—472.

[31? HICKERNELL F J. A generalized discrepancy and quadrature error bound [J]. Mathematics of Computation of the American Mathematical Society，1998，67（221）：299—322.

[32]? 周圣武，李金玉，周長新. 概率論與數理統計[M].2版. 北京：煤炭工業出版社，2007：208—215.

ZHOU S W，LI J Y，ZHOU C X. Probability theory and mathematical statistics [M].2nd ed. Beijin：China Coal Industry Publishing House，2007：208—215.（In Chinese）