基于高斯過程回歸的地下水模型結(jié)構(gòu)不確定性分析與控制

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

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

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

模型參數(shù)不確定性已經(jīng)得到水文及地下水研究者的廣泛認(rèn)同，如束龍倉等[6]考慮了水文地質(zhì)參數(shù)不確定性對地下水補(bǔ)給量計算的影響。目前模型參數(shù)不確定性研究方法主要包括矩方法和蒙特卡洛(MC)方法[7-9]；矩方法通過求解地下水隨機(jī)偏微分方程獲得變量(如水頭)的均值、方差等統(tǒng)計量，該方法計算量小，但一般難以處理復(fù)雜條件的地下水模型(如不規(guī)則的邊界條件、非高斯分布的模型參數(shù))，且不能獲得變量(如滲透系數(shù)、水頭)的完整概率分布信息。MC方法通過在模型參數(shù)的概率分布空間內(nèi)隨機(jī)抽樣控制模擬不確定性，主要研究方法包括通用似然不確定性分析(GLUE)[10]和馬爾科夫鏈蒙特卡洛模擬(MCMC)[11]。Hassan等[12]利用GLUE方法對地下水流模型參數(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)不確定性表示概念模型對真實對象的刻畫誤差，模型校準(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ū)動的結(jié)構(gòu)誤差統(tǒng)計學(xué)習(xí)方法。BMA通過建立多個地下水模型來描述地下水系統(tǒng)，對這些模型的模擬結(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在實際應(yīng)用過程中存在一些問題，如所建立的模型具有主觀性、難以考慮模型間的相關(guān)性、模型先驗權(quán)重的賦值具有主觀性等，從而影響了BMA進(jìn)行不確定性分析的效果。

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

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

1 研究方法

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

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

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ù)θ的后驗分布：

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

(2)

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

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

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

(3)

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

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

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

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

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

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

(4)

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

ε——測量誤差；

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

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

(5)

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

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

Φ——超參數(shù)；

λ——特征長度的超參；

(6)

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

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

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

(7)

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

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

(8)

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

(9)

(10)

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

(11)

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

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

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

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

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

(12)

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

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

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

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

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

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

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

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

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

2.1 海水入侵模型概況

本文建立一個理想條件下三維巖溶裂隙介質(zhì)海水入侵模型，模型裂隙結(jié)構(gòu)見圖1(a)。研究區(qū)在平面上為一個矩形，長2 000 m，寬1 600 m，垂向上厚度為40 m。裂隙網(wǎng)絡(luò)包括4條主裂隙和若干次級裂隙。采用等效多孔介質(zhì)方法對該區(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，兩個觀測井。研究區(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ù)。

考慮實際場地條件下的巖溶裂隙分布復(fù)雜，難以對其空間分布進(jìn)行準(zhǔn)確刻畫，通常對裂隙網(wǎng)絡(luò)進(jìn)行簡化處理，從而導(dǎo)致模型結(jié)構(gòu)誤差。如圖1(b)所示，假設(shè)野外勘察工作僅識別出4條主裂隙，而忽略了次級裂隙分布。基于該簡化的模型結(jié)構(gòu)來描述海水入侵過程時，模擬結(jié)果將會受到模型結(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ù)值模擬不確定性分析

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

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

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

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

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

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

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

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

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

如圖4所示，分別為在不考慮和考慮結(jié)構(gòu)誤差條件下的O1觀測井濃度預(yù)測曲線，可以發(fā)現(xiàn)考慮結(jié)構(gòu)誤差之后，預(yù)測區(qū)間明顯變寬。根據(jù)對模型預(yù)測性能的定量評價指標(biāo)(表2)，可以發(fā)現(xiàn)，在模型識別期，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測ME絕對值和RMSE值均顯著降低，預(yù)測區(qū)間覆蓋率有所提升(由42.9%提升至100%)。在模型驗證期，考慮模型結(jié)構(gòu)誤差之后，模型預(yù)測ME絕對值和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ù)值模擬算例分析

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

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

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

圖5 實驗裝置圖和觀測曲線Fig.5 Experimental appurratus and observation curve

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

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

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

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

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

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

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

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

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

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

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

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

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

圖8 不考慮(a)和考慮(b)結(jié)構(gòu)誤差時上端流出溶液相對濃度預(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對模型結(jié)構(gòu)誤差進(jìn)行統(tǒng)計學(xué)習(xí)，并將DREAMzs算法與GPR算法相結(jié)合，有效提高了地下水?dāng)?shù)值模擬及預(yù)測結(jié)果的可靠性。

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

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