宗成元，康學遠，施小清，吳吉春

(南京大學地球科學與工程學院/表生地球化學教育部重點實驗室，江蘇 南京 210023)

地下水模型有助于地下水資源的開發(fā)利用和污染修復，在地下水資源管理以及地下水環(huán)境風險性評估、預測等方面具有重要意義。模型模擬、預測效果的好壞取決于模型參數(shù)(如滲透系數(shù)、儲水系數(shù)等)精度的高低。然而，復雜的水文地質(zhì)環(huán)境與觀測資料的有限性造成無法預先精確刻畫模型參數(shù)[1]。

地下水模型包括地下水流模型、溶質(zhì)運移模型等，目前模型參數(shù)估計主要使用非線性回歸法和貝葉斯方法，例如，Bachmaf等[2]利用非線性參數(shù)估計軟件PEST[3]推估水文地球化學反應(yīng)的相關(guān)參數(shù)；Carniato等[4]利用貝葉斯方法對反應(yīng)運移模型進行參數(shù)估計，并對殘差相關(guān)性及模型誤差進行討論。模型參數(shù)個數(shù)多、觀測數(shù)據(jù)類型多樣的特點，給傳統(tǒng)估計方法帶來了困難與挑戰(zhàn)。鑒于傳統(tǒng)方法計算效率低和估計效果差，地下水模型參數(shù)估計中引入了數(shù)據(jù)同化方法。數(shù)據(jù)同化方法廣泛用于水文地質(zhì)、石油工程等領(lǐng)域，通過融合多源觀測數(shù)據(jù)，模擬值與觀測值之間的差盡可能小，進而估計模型參數(shù)值。例如，通過同化地下水水位和示蹤劑濃度觀測數(shù)據(jù)推估滲透系數(shù)場[5-6]。參數(shù)推估方法很多，其中包括Kalman filter方法[7](KF，Kalman)，ES方法[8](Ensemble Smoother，VanLeeuwen and Evensen)，集合卡爾曼濾波法[9](Ensemble Kalman Filter，EnKF，Evensen)，iES方法[10](iterative Ensemble Smoother，Chen)。相比于其他方法，Emerick和Reynolds在2013年提出來的Ensemble Smoother with Multiple Data Assimilation(ESMDA)方法可以獲得更好的數(shù)據(jù)擬合、參數(shù)估計效果，且計算成本更低，但其僅可處理服從高斯分布的滲透系數(shù)場[11-14]。

在沖積含水層中，由于巖相的非均質(zhì)分布，滲透系數(shù)一般呈現(xiàn)出明顯的非高斯特性，非高斯特性給參數(shù)推估帶來了困難與挑戰(zhàn)[15]。例如，當滲透系數(shù)的分布受到具有高滲優(yōu)勢流動管道的古河道的影響；或含有連續(xù)的高滲透帶的沖積含水層的影響[16]。多點地質(zhì)統(tǒng)計方法(Multiple-point geo-statistics，MPS)廣泛用于模擬非高斯場，保持參數(shù)場的非高斯特性，通過訓練圖像描述實際含水層的結(jié)構(gòu)。MPS中的Direct Sampling(DS)方法結(jié)合逐點模擬和模式模擬方法的優(yōu)點，利用順序模擬的策略，不需要預掃描訓練圖像和存儲地質(zhì)模式，計算效率較高[17]，但動態(tài)觀測數(shù)據(jù)(如水頭和濃度等)無法融入DS方法的模擬過程。耦合數(shù)據(jù)同化方法與DS方法既可保持參數(shù)場的非高斯特性，又通過融合多源動態(tài)觀測數(shù)據(jù)提高非高斯場的推估精度。Cao等[18]構(gòu)建了iES-DS數(shù)據(jù)同化框架，融合單一類型觀測數(shù)據(jù)(地下水位)，較好地刻畫了非高斯參數(shù)場。

本文基于DS與ESMDA，構(gòu)建一種新的數(shù)據(jù)同化框架(ESMDA-DS)。ESMDA-DS方法結(jié)合了ESMDA與DS的各自優(yōu)勢，通過融合多源動態(tài)觀測數(shù)據(jù)推估非高斯參數(shù)場。構(gòu)建3個理想算例，包括僅融合水位，同時融合水位與濃度，同時融合水位、濃度以及對數(shù)滲透系數(shù)，分析融合不同類型觀測數(shù)據(jù)對推估精度的影響，并對比基于不同類型觀測數(shù)據(jù)推估場的水位及濃度預測效果。

1 研究方法

1.1 水流和溶質(zhì)運移模型

多孔介質(zhì)飽和承壓含水層中地下水非穩(wěn)定流運動的控制方程為[19]：

(1)

式中：K——滲透系數(shù)；

h——水位；

Ss——貯水率；

W*——源匯項；

t——時間。

本文采用MODFLOW[20]求解上述水流問題。

飽和承壓含水層中溶質(zhì)運移的對流—彌散方程為[19]：

(2)

式中：c——溶液中某種組分的濃度；

t——時間；

D——水動力彌散系數(shù)；

u——實際平均流速矢量。

本文基于MT3DMS[21]求解上述溶質(zhì)運移問題。

1.2 ESMDA-DS方法

ESMDA[12]方法是Emerick和Reynolds在2013年提出的一種參數(shù)推估方法；DS方法作為多點地質(zhì)統(tǒng)計方法的一種，具有重構(gòu)諸如非高斯參數(shù)場等復雜地質(zhì)模式方面的優(yōu)勢。本文提出一種新的參數(shù)推估方法，通過構(gòu)建ESMDA與DS的數(shù)據(jù)同化框架，融合多源動態(tài)觀測數(shù)據(jù)解決非高斯場的參數(shù)推估問題。在推估過程中，隨機產(chǎn)生一定數(shù)量的向?qū)c，ESMDA方法基于水位、濃度更新向?qū)c處的參數(shù)值，DS方法基于向?qū)c處更新的參數(shù)值與對數(shù)滲透系數(shù)觀測值更新全場參數(shù)值，保持對數(shù)滲透系數(shù)場的非高斯特性。

ESMDA-DS方法更新參數(shù)場的主要步驟見圖1。

步驟1：基于DS方法生成參照場與樣本場。

步驟2：確定觀測井與向?qū)c數(shù)量及位置。

步驟3：構(gòu)建狀態(tài)向量X：

(3)

式中：p——向?qū)c與觀測點處對數(shù)滲透系數(shù)的集合；

h——運行MODFLOW后的模擬水位或濃度。

步驟4：確定樣本數(shù)量Ne，迭代次數(shù)Na，系數(shù)αi，每次數(shù)據(jù)同化滿足以下限制條件：

(4)

ESMDA中每個數(shù)據(jù)同化的系數(shù)取α1=9.333，α2=7.0，α3=4.0，α4=2.0[12]。

步驟5：每次迭代時，利用更新的參數(shù)從零時刻開始運行耦合模型：

(5)

式中：i——樣本編號；

f——預測；

a——分析。

通過式(5)，根據(jù)模型參數(shù)獲得不同時刻的狀態(tài)變量，同時該變量可在步驟5中用于更新。

步驟6：更新樣本：

(6)

式中：CYD——狀態(tài)向量與模擬值之間的互協(xié)方差矩陣；

CDD——模擬值的自協(xié)方差矩陣；

l——ESMDA的迭代次數(shù)，l=1,2，…，Na；

CD——觀測誤差的協(xié)方差矩陣；

dobs——具有協(xié)方差αlCD噪聲的擾動觀測值；

d——預測數(shù)據(jù)。

步驟7：DS方法基于向?qū)c處的參數(shù)更新值與對數(shù)滲透系數(shù)觀測值更新全場對數(shù)滲透系數(shù)場分布。

ESMDA-DS中，向?qū)c發(fā)揮重要作用，通過向?qū)c概念耦合ESMDA與DS。若不設(shè)置向?qū)c，ESMDA-DS方法中僅利用DS方法，此時僅能保持參數(shù)場的非高斯特性，動態(tài)水位或濃度數(shù)據(jù)無法融合到多點地質(zhì)統(tǒng)計學模擬中。如果模型所有網(wǎng)格均為向?qū)c，ESMDA-DS方法中僅利用ESMDA方法，僅能融合多源觀測數(shù)據(jù)，不會保持對數(shù)滲透系數(shù)場的非高斯特性。通過耦合ESMDA方法和DS方法，既可保持參數(shù)場的非高斯特性，又通過融合多源觀測數(shù)據(jù)更好地推估非高斯含水層，刻畫模型中高滲優(yōu)勢通道。

圖1 ESMDA-DS方法的主要步驟流程圖Fig.1 Flow chart of the proposed ESMDA-DS method

2 算例

構(gòu)建理想算例驗證ESMDA-DS方法推估非高斯參數(shù)場的效果。理想算例為二維平面矩形承壓含水層非穩(wěn)定流模型。模型面積為0.04 km2，x與y方向長度均為200 m，z方向的含水層厚度為5 m，含水層底板高程為0 m。網(wǎng)格均勻剖分，邊長為5 m×5 m，網(wǎng)格數(shù)為40×40。模型初始水位為30 m，南北邊界為隔水邊界，東側(cè)邊界是水位為30 m的給定水頭邊界。為提供同化價值高的觀測數(shù)據(jù)，流場水位波動較大，將西側(cè)邊界設(shè)定為分段賦值的給定流量邊界。模型模擬時長為40 d，前30 d用于數(shù)據(jù)擬合，為進一步評價參數(shù)估計效果，后10 d用于預測，模型預測精度越高，說明模型參數(shù)估計效果越好。模型具體參數(shù)設(shè)置，見圖 2與表 1。如圖 2所示，在研究區(qū)域內(nèi)均勻布設(shè)9口觀測井，提供水位、濃度觀測數(shù)據(jù)，其中1～6號井還提供滲透系數(shù)觀測值。

圖2 模型設(shè)置Fig.2 Model settings

本文共設(shè)置3個理想算例(表2)，Case 1僅同化水位數(shù)據(jù)，Case 2同時同化水位與濃度數(shù)據(jù)，Case 3同時同化水位、濃度數(shù)據(jù)以及滲透系數(shù)觀測值(lnK)。在模型東側(cè)位置設(shè)置線狀持續(xù)泄露的污染源泄露點，設(shè)定單個網(wǎng)格點的源強濃度為1 mg/L?；趨⒄請鲞\行模型，計算得到水位與濃度模擬值，分別加上噪聲作為水位觀測值(h)與濃度觀測值(c)。假設(shè)水位、濃度、對數(shù)滲透系數(shù)的觀測誤差均服從高斯分布，均值均為0，標準差分別為模型最大降深的5%、模擬值的10%、0.1%[22-23]。在研究區(qū)內(nèi)隨機布設(shè)160個向?qū)c，DS方法將基于向?qū)c處參數(shù)值更新全場參數(shù)分布。

表2 算例情景設(shè)置Table 2 Case settings

DS方法是一種基于模式識別的方法，直接根據(jù)訓練圖像獲取其結(jié)構(gòu)特征。與傳統(tǒng)的兩點地質(zhì)統(tǒng)計學方法相比，多點法能夠更好地描述含水層的非高斯結(jié)構(gòu)。本文利用DS方法生成含水層中砂和黏土兩相分布情況，基于表 3中的隨機函數(shù)產(chǎn)生一系列的連續(xù)性變量。

圖3(a)為描述含水層特征的訓練圖像，圖3(b)所示的參照場在訓練圖像的基礎(chǔ)上利用DS方法產(chǎn)生，每一相都以連續(xù)性變量填充。

在推估參數(shù)場過程中，首先，基于DS方法隨機生成200個對數(shù)滲透系數(shù)場作為初始樣本。與產(chǎn)生參照場的方法類似，利用DS方法和圖 3(a)所示的訓練圖像，基于與產(chǎn)生參照場相同的隨機函數(shù)參數(shù)隨機生成200個參數(shù)場集合，從而構(gòu)成lnK場的初始參數(shù)集合。其次，ESMDA融合動態(tài)觀測數(shù)據(jù)更新向?qū)c處的參數(shù)值。最后，DS根據(jù)向?qū)c與滲透系數(shù)觀測點處的參數(shù)值對參數(shù)場進行插值，從而計算對數(shù)滲透系數(shù)的全場分布。

表3 生成砂、黏土兩相的參數(shù)設(shè)置Table 3 Parameter settings for sand and clay

圖3 理想算例的訓練圖像與參照場Fig.3 Training image and reference field of the ideal case

3 結(jié)果與討論

3.1 非高斯對數(shù)滲透系數(shù)場的推估

圖4 lnK場的集合均值及集合方差分布Fig.4 Mean and variance distribution of lnK field注：Case 1融合h，Case 2融合h與c，Case 3融合h、c與lnK。

未融合lnK數(shù)據(jù)的初始場、融合lnK數(shù)據(jù)的初始場以及融合不同類型觀測數(shù)據(jù)lnK場推估結(jié)果的集合均值及集合方差分布見圖4。未融入lnK時，缺乏滲透系數(shù)觀測值，隨機產(chǎn)生的初始樣本場集合均值分布圖為類均質(zhì)特征，基于此初始場構(gòu)建Case 1與Case 2。Case 1僅推估出上側(cè)的一條高滲優(yōu)勢通道，無法準確推估出兩條高滲優(yōu)勢通道?？赡苁怯捎谒徽`差較大，觀測數(shù)據(jù)量不足，同時存在異參同效現(xiàn)象導致推估精度不高。Cao等[18]基于耦合的iES-DS方法融合誤差較小的水位數(shù)據(jù)較好地刻畫了非高斯參數(shù)場。觀測誤差大小代表觀測水平的高低，其直接影響觀測數(shù)據(jù)的可用性，隨著觀測誤差增大，推估精度越低。當觀測數(shù)據(jù)量較少時，同樣無法得到準確的推估結(jié)果，導致推估場與真實場相差較大[24]。Case 2可大致推估出兩條高滲優(yōu)勢通道，表明同時融合多源數(shù)據(jù)可提高參數(shù)場的推估精度[25-26]。融入lnK時，在推估過程中觀測點處對數(shù)滲透系數(shù)值保持不變，初始樣本場集合均值分布圖可刻畫出上側(cè)的高滲優(yōu)勢通道。與未融入lnK相比，其更接近參照場。同時集合方差分布圖下側(cè)顯示較大的不確定性，表明融合對數(shù)滲透系數(shù)觀測值可生成較為精確的初始場?；诖顺跏紙鰳?gòu)建Case 3，與Case 1、Case 2推估結(jié)果相比，Case 3的推估效果更好，表明融合對數(shù)滲透系數(shù)觀測值較好的約束非高斯場的推估精度。

為定量分析滲透系數(shù)場的推估精度，可通過計算均方根誤差(Root Mean Square error，RMSE)和集散離合度(Ensemble Spread，ES)評價數(shù)據(jù)同化效果，RMSE越小說明參數(shù)估計結(jié)果與參照場越接近，ES值下降越快說明同化過程收斂越快：

(7)

(8)

式中：N——模型中網(wǎng)格數(shù)量；

Ne——迭代次數(shù)；

Yi——每個格點的參照值；

var(xi)——每個格點處不同樣本的方差。

圖5為不同算例lnK場推估結(jié)果直方圖，圖6為不同算例lnK場的RMSE與ES基于同化步數(shù)的變化曲線?；趌nK場推估結(jié)果直方圖，驗證了基于ESMDA-DS方法可以保持參數(shù)場的非高斯特性。此外，對比融合單一類型數(shù)據(jù)與多源觀測數(shù)據(jù)的lnK場推估結(jié)果直方圖，RMSE與ES變化曲線圖，表明同時同化水位、濃度及對數(shù)滲透系數(shù)時推估精度更高，同時收斂速度更快。

圖5 lnK場直方圖Fig.5 Histogram of the lnK field注：Case 1融合h，Case 2融合h與c，Case 3融合h、c與lnK。

圖6 lnK場的RMSE與ES隨同化步數(shù)變化曲線Fig.6 RMSE and ES of the lnK field注：Case 1融合h，Case 2融合h與c，Case 3融合h、c與lnK。

3.2 觀測數(shù)據(jù)的擬合與模型預測

利用水位與濃度的擬合、預測精度進一步評價推估效果。

圖7為1、2號井的水位與濃度擬合、預測圖。前300個時間步長(即前30天)為擬合期，后100個時間步長(即后10天)為預測期。在水位預測結(jié)果中，Case 1中1、2號井的水位預測值與真實值相差較大，預測效果較差。Case 2與Case 3相比，兩口井的水位預測精度大致相當，但同時融合3種觀測數(shù)據(jù)可降低推估結(jié)果的不確定性。在濃度預測結(jié)果中，與Case 1、Case 2相比，Case 3聯(lián)合多源數(shù)據(jù)可更準確地刻畫非高斯參數(shù)場，獲得更高的預測精度[27-28]。

Case 2中，融合水位、濃度數(shù)據(jù)后參數(shù)場推估效果與預測精度不理想的原因是：(1)溶質(zhì)運移范圍與時間相關(guān)，在模擬期內(nèi)，高滲優(yōu)勢通道外觀測井內(nèi)溶質(zhì)濃度變化很小或者無變化，數(shù)據(jù)價值不高，反而由于引入額外誤差使同化效果變差。9口井中，1、2、6、7號井位于高滲優(yōu)勢通道內(nèi)，其僅在模擬后期提供濃度信息，有價值的數(shù)據(jù)量較少。(2)在模擬期內(nèi)水位變化幅度較濃度變化幅度大，與濃度數(shù)據(jù)價值相比，水位變化對數(shù)據(jù)同化效果影響更大。由于水位存在較大的觀測誤差，因此即使融合了濃度數(shù)據(jù)，參數(shù)推估效果與預測精度仍未得到較大改善。

4 結(jié)論與展望

本文構(gòu)建了ESMDA-DS耦合框架，基于理想算例，驗證了其對非高斯參數(shù)場的推估效果，進一步探討了不同類型觀測數(shù)據(jù)對推估效果、水位與濃度預測精度的影響。研究表明：

(1)ESMDA-DS方法吸收了DS方法保持對數(shù)滲透系數(shù)場非高斯特性以及ESMDA方法在同化多種數(shù)據(jù)方面的優(yōu)勢，既保留了參數(shù)推估場的非高斯特性，又準確推估出高滲優(yōu)勢通道。

(2)單獨融合水位數(shù)據(jù)時，由于水位誤差較大、數(shù)據(jù)量較少，導致僅能推估出參照場的局部特征，水位、濃度預測精度較低。同時融合水位與濃度數(shù)據(jù)可改善參數(shù)場的推估效果、提高水位與濃度的預測精度，但存在較大水位誤差導致濃度數(shù)據(jù)對參數(shù)場的推估效果提高有限。而同時融合水位、濃度與滲透系數(shù)時可有效提高參數(shù)場的推估精度、水位與濃度的預測精度，同時降低推估結(jié)果的不確定性。

水文地質(zhì)模型包含許多參數(shù)，如滲透系數(shù)、孔隙度、彌散度等，本文為驗證ESMDA-DS方法的有效性，算例僅推估滲透系數(shù)。后續(xù)可進一步融合多源觀測數(shù)據(jù)(如溫度[29])，聯(lián)合同時推估多個水文地質(zhì)參數(shù)，從而有效提高模型反演精度。