鐘樂樂，曾獻(xiàn)奎，吳吉春

(南京大學(xué)地球科學(xué)與工程學(xué)院/表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023)

近年來地下水?dāng)?shù)值模擬已在科學(xué)研究與生產(chǎn)實(shí)踐中得到廣泛應(yīng)用，如地下水資源評(píng)價(jià)、海水入侵治理、地下水污染修復(fù)等方面[1]。實(shí)際地下水系統(tǒng)十分復(fù)雜，限于認(rèn)知程度一般對(duì)其進(jìn)行簡化[2-3]，由此導(dǎo)致了模擬及預(yù)測結(jié)果的不確定性。根據(jù)地下水?dāng)?shù)值模擬過程，其不確定性的來源可以分為模型參數(shù)，模型結(jié)構(gòu)及觀測數(shù)據(jù)的不確定性[4-5]。由于模擬結(jié)果的不確定性直接影響了相關(guān)決策的科學(xué)性，因此有必要對(duì)地下水模型進(jìn)行不確定性的定量分析與控制，從而提高模擬結(jié)果的可靠性。

模型參數(shù)不確定性已經(jīng)得到水文及地下水研究者的廣泛認(rèn)同，如束龍倉等[6]考慮了水文地質(zhì)參數(shù)不確定性對(duì)地下水補(bǔ)給量計(jì)算的影響。目前模型參數(shù)不確定性研究方法主要包括矩方法和蒙特卡洛(MC)方法[7-9]；矩方法通過求解地下水隨機(jī)偏微分方程獲得變量(如水頭)的均值、方差等統(tǒng)計(jì)量，該方法計(jì)算量小，但一般難以處理復(fù)雜條件的地下水模型(如不規(guī)則的邊界條件、非高斯分布的模型參數(shù))，且不能獲得變量(如滲透系數(shù)、水頭)的完整概率分布信息。MC方法通過在模型參數(shù)的概率分布空間內(nèi)隨機(jī)抽樣控制模擬不確定性，主要研究方法包括通用似然不確定性分析(GLUE)[10]和馬爾科夫鏈蒙特卡洛模擬(MCMC)[11]。Hassan等[12]利用GLUE方法對(duì)地下水流模型參數(shù)及模型預(yù)測不確定性進(jìn)行分析。陸樂等[13]利用單分量自適應(yīng)Metropolis(SCAM)采樣算法，在貝葉斯框架下進(jìn)行地下水流模型的參數(shù)不確定分析。Vrugt等[14-16]在原有MCMC算法的基礎(chǔ)上提出DREAM、DREAMzs、DREAM(ABC)、MT-DREAM等算法，提高了參數(shù)不確定性分析的效率，并在水環(huán)境模擬不確定性分析中得到廣泛應(yīng)用。

模型結(jié)構(gòu)不確定性表示概念模型對(duì)真實(shí)對(duì)象的刻畫誤差，模型校準(zhǔn)過程中若忽略模型結(jié)構(gòu)不確定性，將導(dǎo)致參數(shù)補(bǔ)償(參數(shù)不確定性來補(bǔ)償結(jié)構(gòu)誤差)，從而出現(xiàn)參數(shù)過矯正問題，使得模型預(yù)測能力較差[17-18]。當(dāng)前的處理方法主要有貝葉斯模型平均(BMA)[19]和基于數(shù)據(jù)驅(qū)動(dòng)的結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)方法。BMA通過建立多個(gè)地下水模型來描述地下水系統(tǒng)，對(duì)這些模型的模擬結(jié)果進(jìn)行加權(quán)平均處理模型結(jié)構(gòu)的不確定性。目前基于BMA的模型結(jié)構(gòu)不確定性研究較多，Rojas等[20]利用BMA定量分析了地下水?dāng)?shù)值模擬不確定性的組成，Parrish等[21]提出了將粒子濾波算法(PF)與BMA結(jié)合來減小模型結(jié)構(gòu)誤差。然而，BMA在實(shí)際應(yīng)用過程中存在一些問題，如所建立的模型具有主觀性、難以考慮模型間的相關(guān)性、模型先驗(yàn)權(quán)重的賦值具有主觀性等，從而影響了BMA進(jìn)行不確定性分析的效果。

數(shù)據(jù)驅(qū)動(dòng)法(DDM)主要利用機(jī)器學(xué)習(xí)方法(如神經(jīng)網(wǎng)絡(luò)，支持向量基回歸，隨機(jī)森林，高斯過程回歸)，對(duì)模型結(jié)構(gòu)誤差進(jìn)行統(tǒng)計(jì)模擬，即使用某種統(tǒng)計(jì)模型來擬合模型結(jié)構(gòu)誤差。Demissie[22]等通過一個(gè)地下水流數(shù)值模擬案例，分別將人工神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)、決策樹和基于實(shí)例權(quán)重法用于模型結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)，并對(duì)比了這四種方法進(jìn)行模型結(jié)構(gòu)不確定性分析的表現(xiàn)，結(jié)果表明通過數(shù)據(jù)驅(qū)動(dòng)方法進(jìn)行結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)能夠提高模型預(yù)測能力，且人工神經(jīng)網(wǎng)絡(luò)具有較優(yōu)的模型結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)表現(xiàn)。然而，Demissie[22]等將參數(shù)識(shí)別與結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)分開獨(dú)立進(jìn)行，參數(shù)識(shí)別階段未考慮結(jié)構(gòu)誤差會(huì)導(dǎo)致參數(shù)識(shí)別誤差(過矯正)，同時(shí)影響結(jié)構(gòu)誤差的統(tǒng)計(jì)學(xué)習(xí)。因此，需要將模型結(jié)構(gòu)與模型參數(shù)不確定性分析聯(lián)合進(jìn)行。Xu等[23-24]將高斯過程回歸(GPR)分別用于理想與實(shí)際的地下水流模型的結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)，通過DREAM算法同時(shí)識(shí)別地下水流模型和統(tǒng)計(jì)模型參數(shù)。結(jié)果表明，GPR能夠考慮結(jié)構(gòu)誤差時(shí)空相關(guān)性，避免參數(shù)補(bǔ)償，提高模型預(yù)測能力。

本次研究擬使用高斯過程回歸方法識(shí)別模型結(jié)構(gòu)誤差，基于DREAMzs算法的MCMC模擬同時(shí)識(shí)別地下水模型和統(tǒng)計(jì)模型參數(shù)，實(shí)現(xiàn)地下水?dāng)?shù)值模擬不確定性的控制及其定量分析。此外，分別通過理想條件下基巖裂隙區(qū)海水入侵過程與室內(nèi)柱體實(shí)驗(yàn)溶質(zhì)運(yùn)移過程的數(shù)值模擬模型，驗(yàn)證不確定性分析方法的效率與可靠性。

1 研究方法

1.1 貝葉斯參數(shù)不確定性分析

對(duì)于未考慮結(jié)構(gòu)誤差的模型，通常將結(jié)構(gòu)誤差與測量誤差作為一個(gè)整體：

y=M(x,θ)+δ

(1)

式中：y(y1,…,yn)——狀態(tài)向量(如水頭、溫度、濃度等)；

M(M1,…,Mn)——模型函數(shù)；

x(x1,…,xm)——模型輸入向量(m為輸入向量維數(shù))；

θ(θ1,…,θd)——參數(shù)向量；

d——參數(shù)維數(shù)；

δ——誤差向量包括觀測誤差、結(jié)構(gòu)誤差及參數(shù)誤差。

根據(jù)貝葉斯原理，通過觀測數(shù)據(jù)y，推導(dǎo)參數(shù)θ的后驗(yàn)分布：

p(θ|y)∝p(y|θ)p(θ)

(2)

式中：p(θ|y)——參數(shù)θ后驗(yàn)分布；

p(θ)——參數(shù)θ先驗(yàn)分布。

p(y|θ)為參數(shù)θ的似然函數(shù)[25]：

(3)

式中：∑——誤差向量δ的協(xié)方差矩陣。

由于p(θ|y)很難直接求得，一般采用MCMC模擬，抽樣得到模型參數(shù)θ和模型輸出的后驗(yàn)分布，MCMC模擬基本步驟如下：

(1)確定待識(shí)別的模型參數(shù)θ及其先驗(yàn)邊緣分布；

1.2 基于高斯回歸的模型結(jié)構(gòu)不確定性分析

對(duì)于考慮結(jié)構(gòu)誤差的模型，在式(1)上引入了結(jié)構(gòu)誤差項(xiàng)，將結(jié)構(gòu)誤差與測量誤差分開處理[26]：

y=M(x,θ)+b(x,Φ)+ε

(4)

式中：b——結(jié)構(gòu)誤差；

ε——測量誤差；

Φ——高斯回歸超參數(shù)。

高斯過程回歸(GPR)[27]是一種非參數(shù)貝葉斯核函數(shù)回歸方法，該方法能夠?qū)φ`差進(jìn)行統(tǒng)計(jì)學(xué)習(xí)。在高斯過程回歸中，隨機(jī)變量b(x,Φ)服從多元高斯分布即b(x,Φ)～N (μ,C)，μ為均值函數(shù)，C為協(xié)方差函數(shù)。本文均值函數(shù)μ=0，協(xié)方差函數(shù)采用平方指數(shù)型：

(5)

式中：n——觀測數(shù)據(jù)的個(gè)數(shù)；

Ι——指示函數(shù)，若i=j，則Ι=1，否則Ι=0；

Φ——超參數(shù)；

λ——特征長度的超參；

(6)

式(6)中，右側(cè)第一項(xiàng)代表模型輸出對(duì)觀測數(shù)據(jù)的擬合效果，第二項(xiàng)代表對(duì)模型M的復(fù)雜度懲罰，第三項(xiàng)為標(biāo)準(zhǔn)化常數(shù)。

利用GPR進(jìn)行模型結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)時(shí)，需要同時(shí)識(shí)別模型參數(shù)θ與超參數(shù)Φ。假設(shè)θ與Φ相對(duì)獨(dú)立，根據(jù)貝葉斯原理，通過觀測數(shù)據(jù)y，推導(dǎo)出參數(shù)θ與Φ的后驗(yàn)分布，即

p(θ,Φ|y)∝p(y|θ,Φ)p(θ)p(Φ)

(7)

其中，p(y|θ,Φ)為似然值，可由式(6)求得；p(θ,Φ|y)為后驗(yàn)分布；p(θ)為θ先驗(yàn)分布；p(Φ)為Φ先驗(yàn)分布。由于參數(shù)后驗(yàn)p(θ,Φ|y)難以直接求得，本文使用DREAMzs算法結(jié)合GPR識(shí)別θ，Φ的后驗(yàn)分布。

識(shí)別出θ，Φ后驗(yàn)分布之后，可以對(duì)結(jié)構(gòu)誤差進(jìn)行模擬預(yù)測。由于y-M與b*先驗(yàn)聯(lián)合分布服從多元正態(tài)分布，即：

(8)

在已有訓(xùn)練數(shù)據(jù)的條件下，得到b*的后驗(yàn)分布：

(9)

(10)

Cbb(b*)=C**-C*TC-1C*

(11)

μ*——預(yù)測點(diǎn)結(jié)構(gòu)誤差的先驗(yàn)均值；

C*、C**——先驗(yàn)協(xié)方差矩陣；

Cbb(b*)——后驗(yàn)協(xié)方差矩陣。

進(jìn)行多元正態(tài)隨機(jī)抽樣得到b*、ε*，可得模型模擬預(yù)測值：

y*=M*(x*)+b*+ε*

(12)

式中：x*——模型預(yù)測輸入向量(即預(yù)測點(diǎn))；

b*——預(yù)測點(diǎn)結(jié)構(gòu)誤差；

ε*——預(yù)測點(diǎn)觀測誤差；

M*——預(yù)測點(diǎn)模型輸出值；

y*——預(yù)測值。

綜上，結(jié)合GPR進(jìn)行模型結(jié)構(gòu)不確定性分析的基本步驟如下：

(1)確定待識(shí)別的模型參數(shù)θ與超參數(shù)Φ及其先驗(yàn)邊緣分布；

2 海水入侵?jǐn)?shù)值模擬算例分析

針對(duì)理想條件下的巖溶裂隙介質(zhì)海水入侵問題，利用簡化模型結(jié)構(gòu)的地下水模型描述海水入侵過程時(shí)，分別進(jìn)行考慮和不考慮模型結(jié)構(gòu)誤差的海水入侵預(yù)測不確定性分析，驗(yàn)證基于高斯回歸方法模型結(jié)構(gòu)不確定性分析的可靠性。

2.1 海水入侵模型概況

本文建立一個(gè)理想條件下三維巖溶裂隙介質(zhì)海水入侵模型，模型裂隙結(jié)構(gòu)見圖1(a)。研究區(qū)在平面上為一個(gè)矩形，長2 000 m，寬1 600 m，垂向上厚度為40 m。裂隙網(wǎng)絡(luò)包括4條主裂隙和若干次級(jí)裂隙。采用等效多孔介質(zhì)方法對(duì)該區(qū)域進(jìn)行概化，裂隙區(qū)水平滲透系數(shù)為50 m/d，縱向彌散度為30 m，孔隙度為0.2，非裂隙區(qū)水平滲透系數(shù)為5 m/d，縱向彌散度為10 m，孔隙度為0.1，垂直水平滲透系數(shù)之比為0.1，橫縱彌散度比為0.1。模型東側(cè)為海洋邊界即定水頭H=0 m與定濃度邊界C=20 g/L，南、北、西邊界為隔水邊界。研究區(qū)內(nèi)共有4口抽水井，流量均為350 m3/d，兩個(gè)觀測井。研究區(qū)表面均勻接受降水補(bǔ)給，降水入滲系數(shù)為0.001 m/d。初始流場為抽水前穩(wěn)定狀態(tài)下的地下水流場。利用三維、變密度海水入侵?jǐn)?shù)值模擬程序—SEAWAT4[28]建立數(shù)值模型，模擬期設(shè)為1 600 d。將觀測井處的模擬濃度加上高斯白噪聲觀測誤差(均值為0，方差為1)作為本次算例分析的觀測數(shù)據(jù)。

考慮實(shí)際場地條件下的巖溶裂隙分布復(fù)雜，難以對(duì)其空間分布進(jìn)行準(zhǔn)確刻畫，通常對(duì)裂隙網(wǎng)絡(luò)進(jìn)行簡化處理，從而導(dǎo)致模型結(jié)構(gòu)誤差。如圖1(b)所示，假設(shè)野外勘察工作僅識(shí)別出4條主裂隙，而忽略了次級(jí)裂隙分布?；谠摵喕哪Ｐ徒Y(jié)構(gòu)來描述海水入侵過程時(shí)，模擬結(jié)果將會(huì)受到模型結(jié)構(gòu)和模型參數(shù)不確定性的影響。

圖1 理想條件(a)和簡化的(b)研究區(qū)裂隙網(wǎng)絡(luò)Fig.1 Virtual real fracture network and simplified fracture network in the study area

2.2 海水入侵?jǐn)?shù)值模擬不確定性分析

針對(duì)簡化模型結(jié)構(gòu)的地下水模型，分別在考慮和不考慮模型結(jié)構(gòu)誤差的條件下，進(jìn)行海水入侵?jǐn)?shù)值模擬的不確定性分析。假設(shè)地下水模型中的非裂隙區(qū)滲透系數(shù)K1、彌散度α1、裂隙區(qū)滲透系數(shù)K2、彌散度α2為未知參數(shù)，而模型邊界條件、初始條件、模型其它參數(shù)和源匯項(xiàng)與真實(shí)模型相同。模型識(shí)別期為0～1 300 d，模型驗(yàn)證期為1 301～1 600 d。4個(gè)未知模型參數(shù)(K1、K2、α1、α2)和3個(gè)高斯回歸超參數(shù)(σs、λ、σε)的先驗(yàn)分布，見表1。其中，對(duì)于特征長度λ，基于Brynjarsdóttir與O’Hagan的經(jīng)驗(yàn)[29]，其先驗(yàn)分布為Gamma分布，且特征長度λ過小將會(huì)降低模型預(yù)測能力，λ過大則難以學(xué)習(xí)到誤差的統(tǒng)計(jì)特征，因此需要對(duì)λ進(jìn)行限制[27]。此外，本次研究采用DREAMzs算法識(shí)別未知的模型參數(shù)(包括統(tǒng)計(jì)模型參數(shù))時(shí)，設(shè)置4條平行的馬爾科夫鏈。

表1 模型參數(shù)的邊緣先驗(yàn)分布與地下水模型參數(shù)真實(shí)值

不考慮結(jié)構(gòu)誤差時(shí)，MCMC單鏈運(yùn)行510次達(dá)到收斂，取收斂后的8 000次樣本統(tǒng)計(jì)其參數(shù)邊緣后驗(yàn)分布，見圖2(a)。

圖2 不考慮和考慮模型結(jié)構(gòu)誤差識(shí)別得到的各模型參數(shù)邊緣后驗(yàn)分布Fig.2 Identified marginal posterior distributions of groundwater model parameters without and with considering the model structural error

考慮結(jié)構(gòu)誤差時(shí)，MCMC單鏈運(yùn)行1 220次達(dá)到收斂，取收斂后的8 000次樣本統(tǒng)計(jì)其參數(shù)邊緣后驗(yàn)分布，如圖2(b)、圖3所示。

圖3 識(shí)別得到的各高斯回歸超參數(shù)邊緣后驗(yàn)分布Fig.3 Identified marginal posterior distributions of the hyper parameters with Gaussian progress regression

由圖2可以看出，對(duì)于參數(shù)K1，考慮結(jié)構(gòu)誤差時(shí)得到的K1邊緣后驗(yàn)分布在真實(shí)值處具有較大的預(yù)測概率，且考慮結(jié)構(gòu)誤差時(shí)K1真值更加靠近后驗(yàn)分布峰值；對(duì)于參數(shù)K2，不考慮結(jié)構(gòu)誤差時(shí)邊緣后驗(yàn)分布不能包含K2真值，而考慮結(jié)構(gòu)誤差后K2真值落在邊緣后驗(yàn)分布內(nèi)；對(duì)于參數(shù)α1，不考慮結(jié)構(gòu)誤差時(shí)邊緣后驗(yàn)分布同樣不能包含α1的真值，但考慮結(jié)構(gòu)誤差時(shí)α1的真值落在邊緣后驗(yàn)分布內(nèi)，且較靠近峰值；對(duì)于參數(shù)α2，考慮與不考慮結(jié)構(gòu)誤差時(shí)，后驗(yàn)分布均包含α2真值，但考慮結(jié)構(gòu)誤差時(shí)在真實(shí)值處具有更大的預(yù)測概率，且峰值更加靠近真實(shí)值。因此，基于高斯過程結(jié)構(gòu)誤差驅(qū)動(dòng)模型顯著降低了參數(shù)K1、K2、α1、α2補(bǔ)償模型結(jié)構(gòu)誤差的程度。

基于MCMC得到的參數(shù)后驗(yàn)分布樣本，獲得相對(duì)應(yīng)的模型模擬結(jié)果，通過頻率統(tǒng)計(jì)可得到模擬值(或預(yù)測值)的置信區(qū)間與平均值等，如將95%水平下的置信區(qū)間(2.5%～97.5%)作為本次研究的預(yù)測區(qū)間。使用預(yù)測區(qū)間覆蓋率(即觀測數(shù)據(jù)位于預(yù)測區(qū)間的比例)，ME值(即觀測值與預(yù)測值之差的均值)和RMSE值(即觀測值與預(yù)測值之差的均方根)來表征模型的預(yù)測性能。預(yù)測區(qū)間包含率越高，ME值的絕對(duì)值越小，RMSE值越小，預(yù)測值越靠近觀測值，表示預(yù)測性能越好。在考慮和不考慮結(jié)構(gòu)誤差情況下的各指標(biāo)的統(tǒng)計(jì)分析見表2。

表2 模型預(yù)測性能指標(biāo)統(tǒng)計(jì)

如圖4所示，分別為在不考慮和考慮結(jié)構(gòu)誤差條件下的O1觀測井濃度預(yù)測曲線，可以發(fā)現(xiàn)考慮結(jié)構(gòu)誤差之后，預(yù)測區(qū)間明顯變寬。根據(jù)對(duì)模型預(yù)測性能的定量評(píng)價(jià)指標(biāo)(表2)，可以發(fā)現(xiàn)，在模型識(shí)別期，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測ME絕對(duì)值和RMSE值均顯著降低，預(yù)測區(qū)間覆蓋率有所提升(由42.9%提升至100%)。在模型驗(yàn)證期，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測ME絕對(duì)值和RMSE值均有所降低，預(yù)測區(qū)間覆蓋率顯著提升(由0提升至100%)。因此，通過考慮模型結(jié)構(gòu)誤差，能夠提高模型預(yù)測能力。

圖4 不考慮(a)和考慮(b)結(jié)構(gòu)誤差觀測井O1預(yù)測濃度曲線Fig.4 Concentration prediction curve of the observation well (O1) without and with considering the model structural error

3 溶質(zhì)運(yùn)移柱體實(shí)驗(yàn)數(shù)值模擬算例分析

本文以LV等[30]一維溶質(zhì)運(yùn)移柱體實(shí)驗(yàn)為例，建立溶質(zhì)運(yùn)移數(shù)值模型，分別在考慮結(jié)構(gòu)誤差與不考慮結(jié)構(gòu)誤差的條件下，對(duì)水文地質(zhì)參數(shù)進(jìn)行反演并進(jìn)行不確定性分析。

3.1 溶質(zhì)運(yùn)移實(shí)驗(yàn)概況

圖5(a)為該實(shí)驗(yàn)裝置圖，裝置為圓形砂柱，直徑2.6 cm，砂柱高20.0 cm，砂柱由兩部分組成：RegionⅠ直徑為920 μm的石英砂均勻填充，RegionⅡ由直徑為550 μm的石英砂均勻填充。實(shí)驗(yàn)開始前75.4 min，以1 mL/min的流速從砂柱下端注入濃度為C0的KNO3溶液，之后以1 mL/min的流速通入純凈水。從注入KNO3時(shí)開始計(jì)時(shí)，每隔4 min檢測上端流出溶液濃度C，并得到相對(duì)濃度C/C0。在圖5(a)裝置條件下，共得到71個(gè)時(shí)刻(即第1、4、…、281 min時(shí)刻)的KNO3的相對(duì)濃度，以此作為觀測值，得到相對(duì)濃度觀測曲線(圖5b)。

圖5 實(shí)驗(yàn)裝置圖和觀測曲線Fig.5 Experimental appurratus and observation curve

3.2 溶質(zhì)運(yùn)移數(shù)值模擬不確定性分析

對(duì)上述KNO3運(yùn)移實(shí)驗(yàn)過程進(jìn)行概化，利用ModFlow2005[31]和MT3DMS[32]分別建立地下水流模型和溶質(zhì)運(yùn)移模型。模型在平面上剖分為24×24個(gè)格柵，模型垂向上剖分為40層。假設(shè)地下水模型中RegionⅠ滲透系數(shù)K1、縱向彌散度α1和RegionⅡ滲透系數(shù)K2、縱向彌散度為α2為未知參數(shù)。模型圓柱的內(nèi)壁設(shè)為隔水邊界，上下底部均設(shè)為定流量邊界。此外，整個(gè)實(shí)驗(yàn)柱體孔隙度設(shè)為0.33、橫向彌散度與縱向彌散度之比取0.1。

考慮到實(shí)際溶質(zhì)運(yùn)移過程的復(fù)雜性，難以對(duì)含水層結(jié)構(gòu)進(jìn)行精細(xì)刻畫，未考慮吸附化學(xué)反應(yīng)等過程及常規(guī)對(duì)流彌散方程的缺陷[33]，從而產(chǎn)生結(jié)構(gòu)誤差。利用概化的模型模擬溶質(zhì)過程時(shí)，模擬結(jié)果將會(huì)受到模型結(jié)構(gòu)和模型參數(shù)不確定性的影響。

分別在考慮和不考慮模型結(jié)構(gòu)誤差的條件下，進(jìn)行柱體實(shí)驗(yàn)溶質(zhì)運(yùn)移數(shù)值模擬的不確定性分析。模型識(shí)別期為0～117 min，模型驗(yàn)證期為118～201 min。表3所示為4個(gè)未知模型參數(shù)(K1、K2、α1、α2)和3個(gè)高斯回歸超參數(shù)(σs、λ、σε)的先驗(yàn)分布。本次研究采用DREAMzs算法識(shí)別未知的模型參數(shù)(包括統(tǒng)計(jì)模型參數(shù))，設(shè)置4條平行的馬爾科夫鏈。

不考慮結(jié)構(gòu)誤差時(shí)，MCMC單鏈運(yùn)行260次達(dá)到收斂，取收斂后的8 000次樣本統(tǒng)計(jì)其參數(shù)邊緣后驗(yàn)分布，如圖6(a)所示。

考慮結(jié)構(gòu)誤差時(shí)，MCMC單鏈運(yùn)行990次達(dá)到收斂，取收斂后的8 000次樣本統(tǒng)計(jì)其參數(shù)邊緣后驗(yàn)分布(圖6b、圖7)。

表3 模型參數(shù)的邊緣先驗(yàn)分布

對(duì)于真實(shí)的地下水模型，模型參數(shù)的真值未知，因此無法驗(yàn)證參數(shù)后驗(yàn)分布的正確性。對(duì)比圖6(a)、(b)，可以看出兩種方法得到的參數(shù)邊緣后驗(yàn)分布有一定的重合，但是對(duì)于參數(shù)K1與K2，考慮結(jié)構(gòu)誤差時(shí)參數(shù)邊緣后驗(yàn)分布的范圍及峰態(tài)與不考慮結(jié)構(gòu)誤差相比存在較大的差異。同時(shí)對(duì)于參數(shù)α1與α2，兩種方法得到的參數(shù)后驗(yàn)分布的峰態(tài)也存在差異。這表明，是否考慮結(jié)構(gòu)誤差，會(huì)導(dǎo)致模型參數(shù)識(shí)別的差異。

基于MCMC得到的參數(shù)后驗(yàn)分布樣本，計(jì)算相對(duì)應(yīng)的模型模擬結(jié)果，通過頻率統(tǒng)計(jì)可得的模擬值的95%置信區(qū)間(即預(yù)測區(qū)間)與平均值(即預(yù)測值)。對(duì)預(yù)測區(qū)間覆蓋率，ME值和RMSE值進(jìn)行統(tǒng)計(jì)分析，見表4。

表4 模型預(yù)測性能指標(biāo)統(tǒng)計(jì)

圖6 不考慮(a)和考慮(b)模型結(jié)構(gòu)誤差識(shí)別得到的各模型參數(shù)邊緣后驗(yàn)分布Fig.6 Identified marginal posterior distributions of groundwater model parameters without and with considering the model structural error

圖7 識(shí)別得到的各高斯回歸超參數(shù)邊緣后驗(yàn)分布Fig.7 Identified marginal posterior distributions of the hyper parameters with Gaussian progress regression

在不考慮和考慮結(jié)構(gòu)誤差條件下實(shí)驗(yàn)柱體端口處平均濃度的穿透曲線，如圖8所示。可以看出，考慮結(jié)構(gòu)誤差條件下，預(yù)測區(qū)間覆蓋區(qū)域更大。結(jié)合模型預(yù)測性能的定量評(píng)價(jià)指標(biāo)(表4)，可以發(fā)現(xiàn)，在模型識(shí)別期，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測ME絕對(duì)值和RMSE值均明顯降低，預(yù)測區(qū)間覆蓋率顯著提升(由33.3%提升至100%)。在模型驗(yàn)證期，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測ME絕對(duì)值和RMSE值均有所降低，預(yù)測區(qū)間覆蓋率也有所提升(由42.9%提升至100%)。因此，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測值更加接近觀測值，模型具有更好的預(yù)測性能。

圖8 不考慮(a)和考慮(b)結(jié)構(gòu)誤差時(shí)上端流出溶液相對(duì)濃度預(yù)測曲線Fig.8 Relative concentration prediction curve of the solution flowing from the upper without and with considering the model structural error

4 結(jié)論及展望

地下水模型結(jié)構(gòu)受到地層巖性、地質(zhì)構(gòu)造、化學(xué)反應(yīng)過程等多種不確定因素的影響。但囿于人類認(rèn)知與勘測資料，這些不確定性因素?zé)o法避免，因此存在模型系統(tǒng)偏差。本次研究提出一種基于高斯過程回歸(GPR)的模型結(jié)構(gòu)不確定性分析方法。該方法利用GPR對(duì)模型結(jié)構(gòu)誤差進(jìn)行統(tǒng)計(jì)學(xué)習(xí)，并將DREAMzs算法與GPR算法相結(jié)合，有效提高了地下水?dāng)?shù)值模擬及預(yù)測結(jié)果的可靠性。

基于兩個(gè)案例分析—理想巖溶裂隙區(qū)海水入侵模型與室內(nèi)溶質(zhì)運(yùn)移柱體實(shí)驗(yàn)，通過對(duì)比在考慮與不考慮結(jié)構(gòu)誤差的條件下，獲得的參數(shù)邊緣后驗(yàn)分布及預(yù)測分布區(qū)間。結(jié)果發(fā)現(xiàn)：考慮模型結(jié)構(gòu)誤差之后，模型參數(shù)邊緣后驗(yàn)分布更加靠近真實(shí)值，在參數(shù)識(shí)別過程減輕了參數(shù)補(bǔ)償?shù)挠绊?，同時(shí)模型預(yù)測區(qū)間對(duì)觀測數(shù)據(jù)的包含率顯著變大，預(yù)測結(jié)果的平均誤差(ME)絕對(duì)值與均方根誤差(RMSE)值也明顯降低，模型具有更好的預(yù)測性能。因此，通過數(shù)據(jù)驅(qū)動(dòng)模型考慮結(jié)構(gòu)誤差，可以減輕參數(shù)識(shí)別過程中的參數(shù)補(bǔ)償，提高模型的模擬及預(yù)測性能。

在觀測數(shù)據(jù)有限的條件下，考慮結(jié)構(gòu)誤差的數(shù)據(jù)驅(qū)動(dòng)模型方法會(huì)導(dǎo)致預(yù)測區(qū)間寬度一定程度的增加，通過融合其他類型的觀測數(shù)據(jù)，可進(jìn)一步降低模型預(yù)測不確定性。