地下水污染監(jiān)測井優(yōu)化設(shè)計(jì)及污染源識別

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

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

關(guān)鍵詞：監(jiān)測井優(yōu)化;污染源識別;貝葉斯方法;信息熵;最優(yōu)拉丁超立方抽樣;差分進(jìn)化自適應(yīng)Metropolis算法;克里金

中圖分類號：X523? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻(xiàn)標(biāo)志碼：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

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

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

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

算法的改進(jìn)版本，相比于DE-MC算法，DREAM算法采用自適應(yīng)隨機(jī)子空間采樣技術(shù)及能夠自適應(yīng)調(diào)整的交叉概率，并且運(yùn)用IQR統(tǒng)計(jì)方法去除無用鏈，這幾個(gè)方面提高了DREAM算法的搜索效率和解的精度[21].

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

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

1? ?研究方法

1.1? ?貝葉斯公式

貝葉斯公式如下：

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

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

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

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

式（10）左端項(xiàng)含有監(jiān)測值d，但在進(jìn)行監(jiān)測井方案設(shè)計(jì)時(shí)，并未真正得到d，可認(rèn)為監(jiān)測值d是隨機(jī)變量，其概率密度函數(shù)為p（d|MP）. 為了得到一個(gè)僅包含監(jiān)測方案MP的函數(shù)，對式（10）兩側(cè)乘以

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

望為：

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

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

1.3? ?DREAM算法

1.3.1? ?DREAM算法具體步驟

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

1）在模型未知參數(shù)α先驗(yàn)范圍內(nèi)隨機(jī)產(chǎn)生NP

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

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

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

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

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

對？坌i∈{1，2，…，NP}，當(dāng)j = 1時(shí)，令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是區(qū)間[0，1]上的隨機(jī)數(shù).

5）根據(jù)IQR統(tǒng)計(jì)方法[21]去除無用鏈.

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

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

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

DREAM算法后50%采樣過程的收斂性進(jìn)行判斷，判斷指標(biāo)為：

1.4? ?最優(yōu)拉丁超立方抽樣方法和Kriging替代模型

1.4.1? ?最優(yōu)拉丁超立方抽樣方法

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

1.4.2? ?Kriging替代模型

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

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

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

2? ?算例應(yīng)用

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

2.1.1? ?模型建立

假定研究區(qū)域?yàn)榫匦螀^(qū)域，長1 000 m，寬600 m，含水層為厚度35 m的砂質(zhì)潛水含水層（水文地質(zhì)參數(shù)見表1），西部邊界Γ1與東部邊界Γ3為給定水頭邊界，其中東部水頭為25 m，西部水頭為30 m，北部邊界Γ2與南部邊界Γ4為隔水邊界.天然狀態(tài)下地下水流為自西向東流動(dòng)的二維均質(zhì)各向同性的非承壓穩(wěn)定流.

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

以西南角為坐標(biāo)原點(diǎn)建立坐標(biāo)系，根據(jù)研究區(qū)水文地質(zhì)條件，建立地下水水流數(shù)值模型：

式中：Dx、Dy分別為水動(dòng)力彌散系數(shù)在x、y方向上的分量，m2/d;vx、vy分別為x、y方向上地下水的滲流速度，m/d;n為含水層介質(zhì)的孔隙度，無量綱;c為污染物質(zhì)量濃度，mg/L;Qinj為向含水層的注入液體量，m3/d;Cinj為隨液體進(jìn)入含水層的污染物質(zhì)量濃度，mg/L;f（x，y，t）表示僅在水動(dòng)力彌散作用下，單位時(shí)間通過實(shí)際過水?dāng)嗝鎲挝幻娣e的溶質(zhì)質(zhì)量.

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

2.1.2? ?問題概述

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

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

如下：

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）運(yùn)用最優(yōu)拉丁超立方抽樣方法從污染源未

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

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

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

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

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

為進(jìn)一步檢驗(yàn)替代模型的精度，分別運(yùn)用決定系數(shù)，平均絕對誤差及均方根誤差3個(gè)指標(biāo)對7眼監(jiān)測井的替代模型進(jìn)行檢驗(yàn)評價(jià)，結(jié)果如表5所示.

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

2.3? ?監(jiān)測方案優(yōu)化設(shè)計(jì)

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

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

式中：MP表示監(jiān)測方案.

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

式中： j為第j組真實(shí)參數(shù)α;k表示參數(shù)α的第k個(gè)分量. 由式（20）得出Ci7種監(jiān)測方案的反演結(jié)果MRE（MP），如表7所示.

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

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

檢驗(yàn)假設(shè)H0 ： E（MP） 和RT（MP）之間沒有真正的線性關(guān)系，H1 ： E（MP）和RT（MP）之間有線性關(guān)系;顯著性水平α = 0.05.

經(jīng)計(jì)算可得F檢驗(yàn)統(tǒng)計(jì)量F = 40.739，而Fα（1，35-2） = Fα（1，33） = 4.139，因此F > Fα（1，33），拒絕H0 . F檢驗(yàn)表明E（MP） 和RT（MP）之間存在線性關(guān)系. 另外圖3（c）中實(shí)線表示擬合直線，兩側(cè)虛線標(biāo)出其95%置信區(qū)間，只有2個(gè)點(diǎn)落在置信區(qū)間外面，也充分說明了E（MP） 和RT（MP）之間具有正線性關(guān)系.

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

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

E（MP）隨監(jiān)測井?dāng)?shù)量的變化曲線，如圖4所示.

由圖4可知，在各監(jiān)測類型的最優(yōu)方案條件下，信息熵隨著監(jiān)測井?dāng)?shù)量的增加而減小，其中兩眼井的信息熵顯著小于一眼井的信息熵，但是與其他數(shù)量的監(jiān)測井的信息熵相差不大.

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

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

2.4? ?基于優(yōu)化監(jiān)測方案的污染源識別

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

由表9和圖5可知，方案1到方案3參數(shù)后驗(yàn)分布范圍逐漸減小，與相應(yīng)監(jiān)測方案信息熵的變化趨勢一致（圖4），且參數(shù)真值均在參數(shù)后驗(yàn)分布范圍內(nèi). 進(jìn)一步表明監(jiān)測方案信息熵越小，參數(shù)后驗(yàn)分布不確定性越小.

為進(jìn)一步說明3種監(jiān)測方案的參數(shù)反演效果，對馬爾科夫鏈穩(wěn)定后的剩余2 000個(gè)樣本點(diǎn)進(jìn)行后驗(yàn)統(tǒng)計(jì)分析，結(jié)果見表10.

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

監(jiān)測井?dāng)?shù)越多所需費(fèi)用越多，單井監(jiān)測費(fèi)用最少. 方案2與方案3信息熵差距很小，5個(gè)參數(shù)后驗(yàn)樣本均值相對誤差的平均值也近似相等，卻增加了1眼井進(jìn)行監(jiān)測，費(fèi)用明顯增加. 對于本文案例，如果需要綜合考慮監(jiān)測成本及參數(shù)后驗(yàn)分布范圍大小，認(rèn)為方案2為最佳的監(jiān)測井方案.

3? ?結(jié)? ?論

1）與拉丁超立方抽樣方法相比，最優(yōu)拉丁超立方抽樣可有效提高參數(shù)先驗(yàn)分布抽樣的均勻性，避免抽樣結(jié)果的隨機(jī)性及對整個(gè)抽樣空間的覆蓋填充程度的差異性.

2）Kriging替代模型屬于黑箱模型，能以較小的計(jì)算量得到和地下水?dāng)?shù)值模型相近的輸入輸出關(guān)系，在保證模擬精度的條件下，顯著降低計(jì)算負(fù)荷.

3）參數(shù)反演結(jié)果的相對均方根誤差與監(jiān)測方案信息熵呈現(xiàn)較好的正相關(guān)關(guān)系. 信息熵是參數(shù)后驗(yàn)分布不確定性的有效量度，信息熵越小，參數(shù)后驗(yàn)范圍越小，基于貝葉斯公式與信息熵的監(jiān)測方案優(yōu)化設(shè)計(jì)方法是確定地下水污染監(jiān)測方案的有效方法.

4）并非監(jiān)測井?dāng)?shù)量越多越有利于污染源的反

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

參考文獻(xiàn)

[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]? 張江江. 地下水污染源解析的貝葉斯監(jiān)測設(shè)計(jì)與參數(shù)反演方法 [D]. 杭州：浙江大學(xué)環(huán)境與資源學(xué)院，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]? 肖傳寧，盧文喜，趙瑩，等.基于徑向基函數(shù)模型的優(yōu)化方法在地下水污染源識別中的應(yīng)用 [J]. 中國環(huán)境科學(xué)，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 代理模型的優(yōu)化設(shè)計(jì)方法及其在注塑成型中的作用 [D]. 大連：大連理工大學(xué)運(yùn)載工程與力學(xué)學(xué)部，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]? 周圣武，李金玉，周長新. 概率論與數(shù)理統(tǒng)計(jì)[M].2版. 北京：煤炭工業(yè)出版社，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）