基于MCMC和ES-MDA方法的地下水?dāng)?shù)值模型非均質(zhì)參數(shù)場及開采量的反演研究

劉墉達，陳 喜，高 滿，孟詳博，劉維翰，黃日超

(天津大學(xué) 地球系統(tǒng)科學(xué)學(xué)院 表層地球系統(tǒng)科學(xué)研究院，天津 300072)

1 研究背景

對于復(fù)雜的多層含水層，根據(jù)有限的地下水水位等觀測數(shù)據(jù)，反演水文地質(zhì)參數(shù)、開采量等通常存在不唯一性、不確定性問題，且在調(diào)用地下水?dāng)?shù)值模型進行參數(shù)反演時，隨著調(diào)用次數(shù)和參數(shù)維度增加，反演計算成本變高。解決這些問題對提高地下水?dāng)?shù)值模擬精度、增加地下水資源開發(fā)利用評估的可靠性具有重要意義。

模型參數(shù)調(diào)試和反演方法早期較多采用試錯法。隨著計算機發(fā)展，基于最優(yōu)化理論的確定性和隨機優(yōu)化方法得到快速發(fā)展。相對而言，梯度下降法等確定性尋優(yōu)方法對參數(shù)初值要求較高，優(yōu)化結(jié)果容易陷入局部最優(yōu)解。隨機優(yōu)化方法推求參數(shù)的概率分布，特別適用于解決多參數(shù)反演問題[1-2]。如采用蒙特卡羅方法在參數(shù)變化范圍隨機抽樣，構(gòu)建馬爾科夫鏈生成參數(shù)樣本集合，即馬爾科夫鏈蒙特卡羅方法(MCMC)，可以處理高維參數(shù)空間和模擬非線性問題，因此被廣泛應(yīng)用于地下水?dāng)?shù)值模型參數(shù)反演。Yan等[3]基于MCMC算法對理想均質(zhì)含水層非穩(wěn)定流下的污染源位置和污染物泄露時間進行了反演；李雪利等[4]以二維非均質(zhì)含水層區(qū)域為研究算例，基于MCMC方法反演了地下水污染源，且對比了MCMC不同建議分布對參數(shù)反演結(jié)果的影響。MCMC方法需要多次調(diào)用地下水?dāng)?shù)值模型，隨著待反演參數(shù)增多，地下水?dāng)?shù)值模型調(diào)用次數(shù)隨之增加，產(chǎn)生較高的計算代價。

近年來發(fā)展的替代模型(Surrogate model)，通過構(gòu)建一個擬合函數(shù)(如克里金法[3，5](Kriging)、徑向基函數(shù)法[6-7]、支持向量機[8-10]、多項式混沌展開[11]以及貝葉斯網(wǎng)絡(luò)[12]等)，模仿復(fù)雜的物理模型輸入-輸出過程，可大幅減少計算成本。其中，Kriging法采用回歸多項式的確定性分析和不確定隨機分析，描述全局和局部的統(tǒng)計特性，建模方便、精度較高。Yan等[3]針對理想均質(zhì)含水層，采用Kriging法建立了非穩(wěn)定流數(shù)值模型的替代模型，驗證結(jié)果表明建立的替代模型具有較高精度。替代模型與MCMC結(jié)合可提高參數(shù)優(yōu)選的計算效率，如張江江[13]采用自適應(yīng)稀疏格子插值法構(gòu)造了替代模型，并與MCMC結(jié)合優(yōu)選監(jiān)測井分布、反演模型參數(shù)。

另一種參數(shù)反演途徑是數(shù)據(jù)同化，通過融合觀測數(shù)據(jù)與模型預(yù)測值信息，自動調(diào)整模型參數(shù)、狀態(tài)和運行軌跡，實現(xiàn)對模型參數(shù)的反演估計[14]。常用的數(shù)據(jù)同化方法包括集合卡爾曼濾波[15](EnKF)、集合平滑器[13](ES)和多重數(shù)據(jù)同化集合平滑器[16](ES-MDA)等。EnKF使用當(dāng)前時刻觀測值對上一時刻的參數(shù)和狀態(tài)進行更新的實時算法，需要多次調(diào)用地下水模型，計算過程較為繁瑣，且EnKF對模型參數(shù)和狀態(tài)同時更新，可能導(dǎo)致更新后的參數(shù)和狀態(tài)的不一致性。ES作為一種批處理算法，一次性使用所有時刻的觀測值更新模型參數(shù)和狀態(tài)，且只更新參數(shù)以避免參數(shù)和狀態(tài)的不一致性。ES-MDA通過多次迭代進行數(shù)據(jù)同化，可有效地避免了觀測數(shù)據(jù)過擬合，具有調(diào)用模型次數(shù)少、計算成本低的優(yōu)點。陳沖等[16]基于MODFLOW構(gòu)建的黑河流域中游地下水模型，使用ES-MDA方法反演了模型中13個參數(shù)最優(yōu)取值，對比了ESM-DA不同超參數(shù)對反演結(jié)果的影響。周念清等[17]利用ES-MDA方法融合高密度電阻率法采集的觀測數(shù)據(jù)，實現(xiàn)了對污染源參數(shù)和滲透系數(shù)場的聯(lián)合反演。然而ES-MDA也存在不足，如敏感度低的模型參數(shù)反演精確度低。

目前常用的地下水?dāng)?shù)值模擬軟件一般包括GMS、Visual MODFLOW和FEFLOW等。這些軟件具有較為實用的交互界面，可以直接在圖形界面構(gòu)建模型，由于它們無法與優(yōu)化算法直接耦合，在過去的研究中更多傾向于先優(yōu)化后模擬，即先選出待優(yōu)化參數(shù)的可能取值，再將參數(shù)輸入地下水模型進行驗證，還有研究采用松散耦合進行參數(shù)反演[18]。這些方法無法利用計算機內(nèi)存在優(yōu)化算法和模型之間快速共享數(shù)據(jù)，極大地降低了計算速度和求解效率，F(xiàn)loPy[19]是由Mark Bakker等基于Python編寫的第三方庫，可以和其他優(yōu)化算法緊密耦合反演地下水?dāng)?shù)值模型參數(shù)。

針對地下水?dāng)?shù)值模型和參數(shù)反演各方法的優(yōu)缺點，本文以復(fù)雜三維含水層(潛水、承壓水)為案例，假設(shè)某一含水層滲透系數(shù)為異質(zhì)性空間變化場，在各含水層開采量未知情形下，建立地下水?dāng)?shù)值模型和基于Kriging方法的替代模型，模擬含水層分層地下水水頭變化，建立基于替代模型的MCMC、替代模型和數(shù)值模型相結(jié)合的兩階段MCMC以及ES-MDA反演方法，對比各方法對參數(shù)和隨機場反演精度，為復(fù)雜含水層參數(shù)反演方法的選擇提供借鑒。

2 研究方法

2.1 地下水?dāng)?shù)值模擬及參數(shù)反演問題對于潛水、承壓為水多層含水層，地下水流運動的一般表達式為：

(1)

根據(jù)監(jiān)測井地下水水位的觀測值，記為d=(d1，d2，…，dn)，以模型模擬值F(θ)與觀測值d之間誤差(ε=d-F(θ))為目標(biāo)函數(shù)，推求ε最小時的模型參數(shù)θ(如K等)，即為參數(shù)反演問題。

2.2 隨機參數(shù)場的Karhunen-Loève展開對于某一非均質(zhì)孔隙介質(zhì)含水層，如滲透系數(shù)K值空間變化，形成參數(shù)反演的高維問題。相比于裂隙和巖溶介質(zhì)顯著的非均質(zhì)性，孔隙介質(zhì)在一定尺度上可以認為相對連續(xù)并具有自相關(guān)性。對此，可利用Karhunen-Loève(K-L)展開[20]，將K值的空間變化表述為確定性函數(shù)和隨機場的線性疊加，以達到降低維度，即：

(2)

在地下水?dāng)?shù)值模擬中，K的自然對數(shù)(lnK)的協(xié)方差函數(shù)C(S1，S2)服從指數(shù)分布，即任意兩點(S1和S2)滿足：

(3)

式中：η1和η2分別為x和y方向上的自相關(guān)長度；σ2為方差。

(4)

式中：G(θ)為θ的已知全局函數(shù)；Z(θ)為高斯隨機函數(shù)，為系統(tǒng)局部偏差，其滿足E(Z)=0，Var(Z)=σ2，且任一兩點Z(θ)的相關(guān)函數(shù)可表示為：

cov[Z(θ(i))，Z(θ(j))]=σ2R(θ(i)，θ(j))

(5)

式中：R(θ(i)，θ(j))為任意兩個樣本之間的協(xié)方差函數(shù)(常用有線性函數(shù)、高斯函數(shù)、指數(shù)函數(shù)等)。

Kriging替代模型實現(xiàn)步驟如下：

(1)建立模型輸入(參數(shù)θ)-輸出(水頭y)集合。在各參數(shù)先驗分布上抽取N組樣本作為輸入集合，調(diào)用地下水?dāng)?shù)值模型，模擬每個監(jiān)測井不同時段水頭值y作為樣本輸出集合。

(2)Kriging替代模型訓(xùn)練。將步驟(1)樣本的輸入和輸出集合代入替代模型進行訓(xùn)練，建立每個監(jiān)測井的水頭替代模型

2.4 參數(shù)反演方法

2.4.1 貝葉斯方法 從觀測數(shù)據(jù)中推斷出未知參數(shù)的可能取值，即后驗概率分布p(θ∣d)表述為：

(6)

式中：θ為待反演的參數(shù)；d為地下水水頭觀測值；p(θ)為θ的先驗概率密度函數(shù)；p(d∣θ)為似然函數(shù)，表示模擬值與觀測值d的接近程度；p(d)為歸一化積分常數(shù)。

通常根據(jù)野外調(diào)查、試驗等信息得出參數(shù)某一取值范圍，并假設(shè)其服從某一分布，即參數(shù)的先驗分布p(θ)。假設(shè)模型中未知參數(shù)共有n個，即θ=(θ1，θ2，…，θn)，服從均勻分布，則參數(shù)θi的p(θ)為：

(7)

式中ai、bi分別為參數(shù)θi先驗分布的下限和上限。

如果各參數(shù)獨立，總先驗分布P(θ)可定義為：

(8)

根據(jù)水頭觀測值d與模擬值F(θ)之間的總誤差(ε=d-F(θ))，假設(shè)ε=(ε1，ε2，…，εn)服從n維高斯分布，得出ε的條件概率密度函數(shù)(似然函數(shù))p(d∣θ)為：

(9)

由于p(d∣θ)計算中需要調(diào)用正演模型F(θ)，當(dāng)參數(shù)維度較高時，一般無法得到p(θ∣d)的解析解。可采用基于MCMC的概率抽樣方法和ES-MDA數(shù)據(jù)同化方法，推求p(θ∣d)的統(tǒng)計量(如均值、標(biāo)準(zhǔn)差等)，當(dāng)作參數(shù)反演問題的解。

2.4.2 基于MCMC的參數(shù)后驗分布求解 MCMC通過隨機采樣構(gòu)建一個馬爾科夫鏈，當(dāng)馬爾科夫鏈?zhǔn)諗亢蟮牟蓸泳拖喈?dāng)于在后驗分布p(θ∣d)上抽取樣本，利用樣本的統(tǒng)計量就可以得到參數(shù)均值和標(biāo)準(zhǔn)差。

常用實現(xiàn)MCMC方式是Metropolis-Hastings(M-H)算法[21]，其基本思想是通過構(gòu)造一個接受-拒絕準(zhǔn)則，以保證采到的樣本服從參數(shù)后驗分布p(θ∣d)，具體步驟如下：

(1)初始化參數(shù)。采用蒙特卡羅方法，隨機生成一組參數(shù)初值θi，i=0；

(2)生成參數(shù)的候選值。從建議分布(如正態(tài)分布)q(θ)中隨機生成一個候選值θ*。

(3)計算接受概率：

(10)

式中：log(p(d∣θ))為對數(shù)似然函數(shù)；在建議分布為正態(tài)分布時，q(θi，θ*)=q(θ*，θi)。

(4)接受或拒絕候選狀態(tài)。從均勻分布U(0，1)中抽取隨機數(shù)u，若u≤α(θi，θ*)，則令θi+1=θ*，否則令θi+1=θi。

(5)令i=i+1，返回步驟(2)，重復(fù)多次，直到采樣得到足夠的參數(shù)樣本。

2.4.3 基于ES-MDA的參數(shù)后驗分布求解 ES-MDA在每次迭代中對觀測誤差引進一個膨脹系數(shù)αi，以較低的計算代價解決高維非線性的參數(shù)反演問題，便于實現(xiàn)多次數(shù)據(jù)同化。具體步驟如下：

(1)給定數(shù)據(jù)同化次數(shù)Na，集合樣本N，第i次數(shù)據(jù)同化時的膨脹系數(shù)αi，其中αi滿足：

(11)

(2)根據(jù)各參數(shù)的先驗分布p(θ)抽取N個樣本，組成初始參數(shù)集合θ1：

θ1=[θ1，1，θ1，2，…，θ1，N]

(12)

(3)初始i=1，對集合中每個樣本進行正演模型計算，得到模擬值集合：

di，j=F(θi，j)，j=1，2，…，N

(13)

(4)對樣本參數(shù)集合進行更新：

(14)

θi+1，j=θi，j+CMD(CDD+R)-1(duc，j-di，j)

(15)

式中：dobs為觀測值；duc為添加αi擾動后的觀測值；z為修正引入服從高斯分布的隨機誤差；CMD為模型參數(shù)θi與模擬值di之間的協(xié)方差矩陣；CDD與R分別為模擬值和觀測誤差的自協(xié)方差矩陣。

(5)i=i+1，重復(fù)執(zhí)行步驟(3)和(4)直到完成Na次數(shù)據(jù)同化，得到最終樣本集合θNa=[θNa，1，θNa，2，…，θNa，N]，以及參數(shù)的后驗估計。

2.4.4 反演參數(shù)的精度評估 根據(jù)各參數(shù)的后驗分布，統(tǒng)計各參數(shù)的均值及標(biāo)準(zhǔn)差。對比各參數(shù)均值與真值，兩者越接近說明反演的參數(shù)誤差越??；同時對比不同反演方法各參數(shù)的標(biāo)準(zhǔn)差，標(biāo)準(zhǔn)差越小，說明參數(shù)反演結(jié)果的不確定性越小。

對于反演的參數(shù)空間分布場，采用均方根誤差(RMSE)和結(jié)構(gòu)相似性指標(biāo)(SSIM)作為評估依據(jù)，其定義如下：

(16)

(17)

式中：u和v分別為參數(shù)標(biāo)準(zhǔn)場和反演參數(shù)場；Np為網(wǎng)格個數(shù)；μ和σ分別為參數(shù)場均值和標(biāo)準(zhǔn)差；C1=0.01和C2=0.03是穩(wěn)定計算的常數(shù)。SSIM的閾值范圍為[-1，1]，SSIM值越接近1，反演參數(shù)場與參數(shù)標(biāo)準(zhǔn)場越趨近一致。

3 研究算例

3.1 研究區(qū)及模型構(gòu)建

3.1.1 地下水?dāng)?shù)值模型及參數(shù)設(shè)置 假定研究區(qū)為邊長1000 m、厚度50 m的5層立方體，其中第1層為潛水含水層，第3、5層為承壓含水層，第2、4層為隔水層。西部邊界與東部邊界為定水頭邊界，水頭分別為10 m和5 m；北部邊界與南部邊界為隔水邊界(圖1)。

圖1 研究區(qū)平面和剖面示意圖

將研究區(qū)剖分為50行、50列、5層有限差分網(wǎng)格，單元邊長為20 m、高為10 m，模擬期10 d，計算時間步長1 d，所構(gòu)建的模型具有一定的代表性。在本案例中模型的不確定性來自第一層承壓含水層非均質(zhì)滲透系數(shù)場和5個開采井的開采量，其他參數(shù)給定，潛水含水層水平滲透系數(shù)(5 m/d)、第二層承壓含水層水平滲透系數(shù)(4 m/d)、隔水層水平滲透系數(shù)(1×10-5m/d)、給水度(0.2)、入滲補給強度(1×10-6m/d)和承壓含水層儲水系數(shù)(1×10-4)，假設(shè)垂向與水平向滲透系數(shù)比值為0.1，初始水頭分布為開采量是0時的穩(wěn)定流水頭。

研究區(qū)內(nèi)共有5個開采井(W1—W5)和10個監(jiān)測井(O1—O10)，其中W1—W3開采第一承壓含水層，W4開采潛水含水層，W5開采第二承壓含水層；O1監(jiān)測潛水含水層水位，O2—O9監(jiān)測第一承壓含水層水位，O10監(jiān)測第二承壓含水層水位。本次研究開采井分散布設(shè)，以減少各開采井的相互影響。監(jiān)測井布設(shè)在開采井附近，旨在捕捉開采井附近水頭變化情況。

使用Flopy調(diào)用地下水?dāng)?shù)值模型MODFLOW，計算逐網(wǎng)格單元水頭，選擇10個監(jiān)測井所在位置10 d逐日水頭，作為參數(shù)反演的觀測值。后續(xù)將10個監(jiān)測井的模擬值與觀測值的誤差最小作為目標(biāo)函數(shù)，使用基于MCMC和ES-MDA方法反演模型參數(shù)。

表1 參數(shù)的先驗分布區(qū)間和真值

表2 開采井開采量 單位：m3/d

3.2 基于MCMC方法的參數(shù)反演為了權(quán)衡模型計算成本與反演精度，本研究對比基于替代模型的MCMC參數(shù)反演與兩階段的MCMC參數(shù)反演結(jié)果。前者在參數(shù)整個變化范圍內(nèi)抽取大數(shù)據(jù)樣本(鏈條長度為50 000，且建議分布方差較大)，基于替代模型的MCMC確定參數(shù)的后驗分布；后者在前者基礎(chǔ)上根據(jù)替代模型參數(shù)后驗分布縮小參數(shù)取值范圍(鏈條長度為10 000，且建議分布的方差較小)，使用基于地下水?dāng)?shù)值模型的MCMC反演參數(shù)后驗分布。假設(shè)誤差ε服從正態(tài)分布ε～N(0，0.1)，使用馬爾科夫鏈接受采樣中最后的1 000組樣本，反演的各參數(shù)后驗分布。

基于替代模型MCMC和兩階段的MCMC參數(shù)反演結(jié)果統(tǒng)計值見表3，以q1、q2、q3為例，反演的參數(shù)后驗分布如圖3。可以看出，兩種反演方法得出的各參數(shù)均值與表1中參數(shù)真值十分接近。相對而言，兩階段的MCMC各參數(shù)均值與真值更為接近，且標(biāo)準(zhǔn)差更小，反演的后驗分布峰值也更集中于真值。

圖3 基于替代模型的MCMC和兩階段的MCMC反演參數(shù)的后驗分布

圖4分別為(a)參數(shù)K3標(biāo)準(zhǔn)場、(b)反演前隨機生成的初始參數(shù)場、(c)基于替代模型和(d)兩階段的MCMC反演K3參數(shù)場?；谔娲Ｐ秃蛢呻A段的MCMC反演參數(shù)場與參數(shù)標(biāo)準(zhǔn)場之間的RMSE分別為8.3×10-2和2.6×10-2；SSIM分別為0.88和0.98，兩種方法反演的參數(shù)場都較為準(zhǔn)確，但兩階段的MCMC反演的參數(shù)場精度更高。說明僅使用替代模型的MCMC參數(shù)反演，雖然計算費時少，但替代模型的誤差影響反演參數(shù)值的可靠性。而兩階段的MCMC通過第一階段縮小參數(shù)范圍，第二階段再調(diào)用數(shù)值模型，可在實現(xiàn)減小計算成本的同時提高反演參數(shù)的精度。

3.3 基于ES-MDA方法的參數(shù)反演在使用ES-MDA方法反演地下水模型參數(shù)過程中，存在某些超參數(shù)的選擇(如同化次數(shù)Na、集合樣本N、膨脹系數(shù)αi)，直接影響同化算法對模型參數(shù)的更新效果。本文借鑒前人經(jīng)驗[16]，取Na=2，N=100，α1=α2=2，采用模型參數(shù)的先驗分布(表1)，使用ES-MDA方法反演模型參數(shù)的后驗分布。

ES-MDA兩次同化反演參數(shù)結(jié)果的對比見表4，以q1、q2、q3為例，反演參數(shù)后驗分布見圖5。第一、二次同化反演結(jié)果各參數(shù)均值與表1中參數(shù)真值十分接近。但第二次同化各參數(shù)均值與真值更為接近，且標(biāo)準(zhǔn)差均小于第一次同化標(biāo)準(zhǔn)差，反演的后驗分布峰值位于真值附近(圖5)。

表4 ES-MDA兩次同化模型參數(shù)統(tǒng)計值

圖5 基于ES-MDA方法兩次同化反演參數(shù)的后驗分布

圖6分別為(a)參數(shù)K3標(biāo)準(zhǔn)場、(b)反演前隨機生成的初始參數(shù)場、(c)基于ES-MDA第一次、(d)第二次同化反演的參數(shù)場。兩次同化反演參數(shù)場與參數(shù)標(biāo)準(zhǔn)場之間的RMSE分別為1.8×10-2和9.1×10-3，SSIM分別為0.93和0.96。兩者都較為準(zhǔn)確，且第二次同化反演參數(shù)場精度更高。說明增加數(shù)據(jù)同化次數(shù)可以增加ES-MDA方法的參數(shù)反演精度，第二次同化反演的后驗分布更為精確。

圖6 滲透系數(shù)場K3

3.4 正演模型驗證在非均質(zhì)參數(shù)場和開采量的校驗后，還要進行正演模擬，用來驗證預(yù)報水頭場的能力。通過將基于兩階段的MCMC和ES-MDA所反演出的非均質(zhì)參數(shù)場和開采量輸入到地下水?dāng)?shù)值模型中，與監(jiān)測井真實水頭值和真實水頭場進行對比。如圖7所示，其中紅色實線為真實水頭，綠色虛線和藍色點劃線分別為基于兩階段的MCMC和基于ES-MDA參數(shù)反演后正演模擬得到的監(jiān)測井水頭。

圖7 基于兩階段的MCMC和基于ES-MDA反演后正演模擬得到的監(jiān)測井水頭與真實值的水頭過程線

根據(jù)圖7可以看出，基于兩種方法反演的非均質(zhì)參數(shù)場和開采量進行正演模擬的監(jiān)測井水頭過程線與真實水頭過程線十分接近，兩種方法都充分捕捉到了監(jiān)測井位置的水頭變化情況，說明監(jiān)測井水頭模擬效果較好。

圖8—12展示了基于兩種方法反演的非均質(zhì)參數(shù)場和開采量在進行正演模擬不同時間節(jié)點的第一承壓含水層水頭場以及與真實水頭場之間的誤差。其中圖8表示在同一個應(yīng)力周期3個不同時間節(jié)點(3 d，6 d，9 d)的真實水頭場。在相同時間節(jié)點，圖9表示基于兩階段的MCMC反演后正演模擬得到的水頭場，圖10表示基于兩階段的MCMC反演后正演模擬得到的水頭場與真實水頭場的之間誤差場。圖11表示基于ES-MDA反演后正演模擬得到的水頭場，圖12表示基于ES-MDA反演后正演模擬得到的水頭場與真實水頭場的之間誤差場。

圖8 第一承壓含水層真實水頭場

圖9 基于兩階段的MCMC反演后正演模擬得到的水頭場

圖10 兩階段的MCMC正演場與真實場之間的誤差場

圖11 基于ES-MDA反演后正演模擬得到的水頭場

圖12 ES-MDA正演場與真實場之間的誤差場

由基于兩階段的MCMC和ES-MDA正演水頭場與真實水頭場的誤差場可看出，基于兩階段的MCMC反演后在正演模擬水頭場誤差較小，誤差值基本控制在0.1 m以內(nèi)?；贓S-MDA反演后在正演模擬水頭場效果一般，在開采井附近誤差較大，誤差值可達到0.4 m，局部水頭場模擬與真實水頭場相差較大，這與ES-MDA對開采量反演精度較低有關(guān)。兩階段的MCMC反演準(zhǔn)確度更高，而且更為均勻，在不同時期對水頭場的模擬都相較于ES-MDA更為準(zhǔn)確。

3.5 不同反演方法結(jié)果對比分析對于開采量(q1—q5)反演結(jié)果，不同方法反演的參數(shù)后驗分布，均值及峰值均接近真值；從非均質(zhì)產(chǎn)生的高維滲透系數(shù)K3場反演精度來看，與K3實際分布相比，不同方法都能較好地反演出參數(shù)隨機場的空間分布，RMSE都較?。灰許SMI為評價指標(biāo)時，基于替代模型的MCMC方法、ES-MDA方法反演的參數(shù)場精度低于兩階段的MCMC反演參數(shù)場精度。說明先使用替代模型的MCMC反演縮小參數(shù)場分布范圍，再調(diào)用地下水?dāng)?shù)值模型，可調(diào)和復(fù)雜模型參數(shù)反演中降低計算成本與提高反演精度之間的矛盾。

對比計算成本，基于替代模型的MCMC方法需要先調(diào)用1000次數(shù)值模型用于訓(xùn)練Kriging替代模型，再調(diào)用50 000次替代模型進行后驗分布采樣；兩階段的MCMC方法還需再調(diào)用10 000數(shù)值模型進行MCMC精確采樣；而ES-MDA方法僅設(shè)定2次同化次數(shù)，集合樣本為100，總共調(diào)用數(shù)值模型次數(shù)為300(1個初始集合加2個同化集合)。因此，ES-MDA方法計算效率更高，而兩階段的MCMC方法計算效率較低。由此可見，在進行地下水?dāng)?shù)值模型參數(shù)反演時，可以根據(jù)參數(shù)反演精度和計算效率要求，選擇合適的反演方法。

4 結(jié)論

基于設(shè)定的多層含水層水文地質(zhì)條件，建立了地下水?dāng)?shù)值模型和Kriging替代模型?；谪惾~斯理論，使用隨機抽樣的MCMC算法和數(shù)據(jù)同化的ES-MDA算法，反演模型非均質(zhì)滲透系數(shù)場及開采量。對比替代模型-MCMC方法、替代模型-地下水?dāng)?shù)值模型的兩階段MCMC以及ES-MDA同化算法反演的模型參數(shù)精度及計算效率。得出如下主要結(jié)論：(1)Kriging替代模型可以有效地模仿地下水?dāng)?shù)值模型模擬的水頭，由此可減小模型參數(shù)反演中計算量。但對于高維滲透系數(shù)場，如僅根據(jù)建立的有限觀測值(如本文中10個監(jiān)測井地下水水頭)，基于替代模型MCMC反演的滲透系數(shù)場誤差相對較大。(2)在基于替代模型MCMC反演的參數(shù)后驗分布基礎(chǔ)上縮小參數(shù)變化范圍，再采用MCMC降低地下水?dāng)?shù)值模型調(diào)用次數(shù)，這種兩階段MCMC反演的參數(shù)精度更高，不確定性小。(3)ES-MDA算法反演的參數(shù)后驗分布計算效率最高，但對高維參數(shù)場和開采量的反演效果不及兩階段的MCMC。

本文僅基于設(shè)定的含水層水文地質(zhì)條件，并采用地下水?dāng)?shù)值模型模擬的各含水層中地下水水頭作為真值，實際水文地質(zhì)條件更為復(fù)雜，且觀測的地下水水頭存在誤差；因此，未來可利用實際區(qū)域觀測資料，權(quán)衡不同反演方法的精度及計算效率來選擇合適的方法，本研究為大范圍地下水?dāng)?shù)值模擬和地下水資源評價等提供適宜的方法。