基于高斯過程回歸的地下水模型結構不確定性分析與控制

鐘樂樂，曾獻奎，吳吉春

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

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

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

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

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

本次研究擬使用高斯過程回歸方法識別模型結構誤差，基于DREAMzs算法的MCMC模擬同時識別地下水模型和統計模型參數，實現地下水數值模擬不確定性的控制及其定量分析。此外，分別通過理想條件下基巖裂隙區海水入侵過程與室內柱體實驗溶質運移過程的數值模擬模型，驗證不確定性分析方法的效率與可靠性。

1 研究方法

1.1 貝葉斯參數不確定性分析

對于未考慮結構誤差的模型，通常將結構誤差與測量誤差作為一個整體：

y=M(x,θ)+δ

(1)

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

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

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

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

d——參數維數；

δ——誤差向量包括觀測誤差、結構誤差及參數誤差。

根據貝葉斯原理，通過觀測數據y，推導參數θ的后驗分布：

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

(2)

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

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

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

(3)

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

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

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

1.2 基于高斯回歸的模型結構不確定性分析

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

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

(4)

式中：b——結構誤差；

ε——測量誤差；

Φ——高斯回歸超參數。

高斯過程回歸(GPR)[27]是一種非參數貝葉斯核函數回歸方法，該方法能夠對誤差進行統計學習。在高斯過程回歸中，隨機變量b(x,Φ)服從多元高斯分布即b(x,Φ)～N (μ,C)，μ為均值函數，C為協方差函數。本文均值函數μ=0，協方差函數采用平方指數型：

(5)

式中：n——觀測數據的個數；

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

Φ——超參數；

λ——特征長度的超參；

(6)

式(6)中，右側第一項代表模型輸出對觀測數據的擬合效果，第二項代表對模型M的復雜度懲罰，第三項為標準化常數。

利用GPR進行模型結構誤差統計學習時，需要同時識別模型參數θ與超參數Φ。假設θ與Φ相對獨立，根據貝葉斯原理，通過觀測數據y，推導出參數θ與Φ的后驗分布，即

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

(7)

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

識別出θ，Φ后驗分布之后，可以對結構誤差進行模擬預測。由于y-M與b*先驗聯合分布服從多元正態分布，即：

(8)

在已有訓練數據的條件下，得到b*的后驗分布：

(9)

(10)

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

(11)

μ*——預測點結構誤差的先驗均值；

C*、C**——先驗協方差矩陣；

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

進行多元正態隨機抽樣得到b*、ε*，可得模型模擬預測值：

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

(12)

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

b*——預測點結構誤差；

ε*——預測點觀測誤差；

M*——預測點模型輸出值；

y*——預測值。

綜上，結合GPR進行模型結構不確定性分析的基本步驟如下：

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

2 海水入侵數值模擬算例分析

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

2.1 海水入侵模型概況

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

考慮實際場地條件下的巖溶裂隙分布復雜，難以對其空間分布進行準確刻畫，通常對裂隙網絡進行簡化處理，從而導致模型結構誤差。如圖1(b)所示，假設野外勘察工作僅識別出4條主裂隙，而忽略了次級裂隙分布。基于該簡化的模型結構來描述海水入侵過程時，模擬結果將會受到模型結構和模型參數不確定性的影響。

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

2.2 海水入侵數值模擬不確定性分析

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

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

不考慮結構誤差時，MCMC單鏈運行510次達到收斂，取收斂后的8 000次樣本統計其參數邊緣后驗分布，見圖2(a)。

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

考慮結構誤差時，MCMC單鏈運行1 220次達到收斂，取收斂后的8 000次樣本統計其參數邊緣后驗分布，如圖2(b)、圖3所示。

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

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

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

表2 模型預測性能指標統計

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

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

3 溶質運移柱體實驗數值模擬算例分析

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

3.1 溶質運移實驗概況

圖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 溶質運移數值模擬不確定性分析

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

考慮到實際溶質運移過程的復雜性，難以對含水層結構進行精細刻畫，未考慮吸附化學反應等過程及常規對流彌散方程的缺陷[33]，從而產生結構誤差。利用概化的模型模擬溶質過程時，模擬結果將會受到模型結構和模型參數不確定性的影響。

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

不考慮結構誤差時，MCMC單鏈運行260次達到收斂，取收斂后的8 000次樣本統計其參數邊緣后驗分布，如圖6(a)所示。

考慮結構誤差時，MCMC單鏈運行990次達到收斂，取收斂后的8 000次樣本統計其參數邊緣后驗分布(圖6b、圖7)。

表3 模型參數的邊緣先驗分布

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

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

表4 模型預測性能指標統計

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

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

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

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

4 結論及展望

地下水模型結構受到地層巖性、地質構造、化學反應過程等多種不確定因素的影響。但囿于人類認知與勘測資料，這些不確定性因素無法避免，因此存在模型系統偏差。本次研究提出一種基于高斯過程回歸(GPR)的模型結構不確定性分析方法。該方法利用GPR對模型結構誤差進行統計學習，并將DREAMzs算法與GPR算法相結合，有效提高了地下水數值模擬及預測結果的可靠性。

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

在觀測數據有限的條件下，考慮結構誤差的數據驅動模型方法會導致預測區間寬度一定程度的增加，通過融合其他類型的觀測數據，可進一步降低模型預測不確定性。