基于替代模型的地表水地下水耦合模型不確定性分析

韓 玉，盧文喜，張將偉，安永凱

(1.吉林大學 地下水與資源環境教育部重點實驗室,吉林 長春 130012；2.吉林大學 環境與資源學院,吉林 長春 130012)

0 引 言

地表水和地下水是水資源的重要組成部分，兩者之間相互作用，具有密切的水力聯系。在對區域水資源的評價管理過程中，兩者之間的聯系常常被簡化，為解決這一問題，本文建立了地表水地下水耦合模擬模型對地表水和地下水的水流運動及兩者之間的水力聯系進行描述。與此同時，在對源匯項的概化過程中，由于主觀與客觀因素的影響，使得人工開采的不確定性較大，對研究區地表水地下水的耦合模型輸出結果會產生一定的影響。為評估模擬模型的輸出結果的可靠性，考慮地下水開采量的不確定性及可能造成的風險對于區域水資源的開發利用具有重要的意義。

近年來，模型的不確定性分析已經得到了人們的廣泛關注。1968年，Yen等人對水文模型的不確定性進行分類[1]，為不確定性的研究提供理論依據；Lopes等人將降雨的空間差異性與流域徑流量的變化情況進行結合，研究降雨空間分布的不確定性[2]，具有實際意義；2012年，Harp和Vesselinov采用水文地質接受概率的方法對水文地質統計模型的不確定性進行綜合分析，為不確定性分析提供一個新的方向[3]；苗添升在2016年結合實際例子對研究區地下水水流模擬模型的不確定性進行分析[4]。但不確定性分析大多應用于地下水水流的數值模擬中[5]，對于地表水地下水的耦合模擬方面的應用較少。在耦合模擬模型的不確定分析過程中，模擬模型的調動次數過多，計算時間較長。為解決這一問題，應用Kriging模型建立地表水地下水耦合模擬模型的替代模型，解決了不確定分析過程中計算負荷較大的問題。

在建立地表水地下水耦合模擬數學模型的基礎上，對地下水開采量進行不確定性分析。采用拉丁超立方抽樣抽取均勻分布的地下水抽水量數據集，將其帶入已建立的地下水地表水耦合模擬模型中，運行HydroGeoSphere軟件得到相應的地表水徑流量和地下水水位，通過輸入輸出數據集對Kriging模型進行學習訓練，從而得到地表水地下水耦合模擬模型的替代模型。按照相同的方式抽取樣本容量足夠大的地下水井開采量數據，借助替代模型得到相應的地表徑流量和地下水水位輸出數據，借助切比雪夫不等式，計算不同置信度下地表徑流量與地下水水位的置信區間。最后通過繪制累積概率分布圖，結合地下水生態水位閾值對研究區的地下水水位進行風險評估。

1 研究方法

1.1 模型建立

建立一個三維的非穩定潛水含水層的假想研究區，研究區兩側高，中間低，河流自南側流入，北側流出。

1.1.1 概念模型

(1)邊界概化。含水層均質各向同性，北側Γ1與下游湖泊相連接，水位穩定為28 m，南側Γ3與上游河流連同，地下水水位保持在68 m左右，東側Γ4和西側Γ2為隔水邊界；對于地表水水流，上游河口S1為定流量邊界，下游河口S2為臨界水深邊界條件。

圖1 研究區邊界示意圖Fig.1 The schematic diagram of study area boundary

(2)參數分區及賦值。對研究區的地表水水流和地下水水流進行參數分區，分區情況如圖2所示，具體賦值情況見表1。

表1 研究區地下水水文地質參數值Tab.1 Groundwater Hydrogeology parameter initial values of study area

(3)源匯項。研究區接受大氣降水補給，研究區的主要排泄方式為蒸發排泄和人工開采?，F在研究區內設置3口抽水井，2個地表水地表水觀測點和21口地下水觀測井，對地表徑流量與地下水水位進行數據的采集，各井的分布情況如圖1所示。模擬時長為3個月，為降水較為充沛的5，6，7月，本次研究對模型運行3個月后的結果進行分析計算。

表2 研究區地表水力參數值Tab.2 Initial values of surface hydraulic parameter

1.1.2 數學模型

采用三維Richard方程和擴散波近似的二維平均深度圣維南(Saint Venant)方程描述地下水水流和地表水水流運動，采用雙重節點法對地表水與地下水模型進行耦合，通過達西流關系描述地表水系統與地下水系統的交換量，得到如下微分方程，定解條件和耦合方程：

圖2 地表水與地下水的參數分區圖Fig.2 Distribution map of hydrogeology parameter of groundwater and surface water

圖3 研究區抽水井與觀測井分布圖Fig.3 Distribution of pumping and observation wells in the study area

(1)

式中：K為滲透系數，[LT-1]；kr為相對滲透率，用以表征包氣帶對地表水地下水交換的影響；H為地下水的水頭，[L]；h1為已知地下水水位，[L]；h為地表水水頭，[L]；qgs為地表水與地下水的交換量，[T-1]；QG為地下水的源匯項，[L3L-3T-1]；Sw為飽和度；SS為儲水系數，[L-1]；θs為孔隙度；Kox,Koy為x,y方向上的地表傳導系數，[LT-1]；ho為地表水面高程，[L]；Q為地表水的源匯項，[LT-1]；φ為地表水等效孔隙度；Kn為邊界法方向的地表傳導系數，[LT-1]；qo1為第二類邊界條件，[LT-1]；qo2為臨界深度邊界條件，[LT-1]；h為邊界水頭，[L]；S1,S2為地表水邊界；Γ1,Γ2,Γ3,Γ4為地下水邊界；lexch為地表水地下水耦合長度，[L]。

1.2 Kriging方法

替代模型也稱為“元模型”，能夠以簡單的程序算法替代復雜繁瑣的系統運行過程，從而大幅度的減少實際系統的運算負荷。Kriging方法是一種通過已知點的相關信息來預測未知點的插值方法[6]，能夠較好的反映模擬模型的輸入輸出關系。Kriging方法的輸入為抽水井的開采量，輸出則為地表觀測點的徑流量與地下觀測井的水位。

Kriging模型的數學表達式是：

y(x)=βfT(xi)+ε(xi)

(2)

式中：f(xi)=[f1(x),…,fk(x)]T為已知的回歸函數的基函數，β=(β1,…,βk)T為相應的未知參數，ε(xi)為高斯隨機過程，其均值為0，方差為σ2，協方差為：

Cov[ε(xi),ε(xj)]=σ2R(xi,xj)

(3)

式中：R(xi,xj)為空間采樣點xi,xj的相關函數，可以采用以下函數形式進行表示：

(4)

由已知n個樣本點[X]=({X1}T,…, {Xn}T)T的響應[y]=(y1,…,yn)T來預測任一點{x*}的響應：

(5)

式中：rT(x)=[R(x*,x1),…,R(x*,xn)]，[R]為矩陣元素為R(xi,xj)的(n×n)矩陣。

1.3 Monte Carlo模擬

Monte Carlo方法又稱為統計模擬方法，是一種常用的數值模擬的不確定性分析方法，在解決非線性的隨機問題的過程中具有一定的優勢[7]。

(1)采用拉丁超立方抽樣方法對地下水開采量進行隨機抽樣，將樣本數據作為輸入數據集；

(2)在地表水地下水耦合模擬模型的基礎上，通過HydroGeoSphere軟件進行并行同步求解，得到輸出量(如地表水徑流量和地下水水位)的數據集；

(3)對輸出結果進行統計分析，得到方差，均值，極大值，極小值等估計量，并用統計結果繪制圖表，根據概率分布情況進行不確定性分析。

2 結果與討論

2.1 Kriging替代模型

采用拉丁超立方抽樣方法得到50組均勻分布的開采量樣本數據，將其作為輸入值，運行HydroGeoSphere軟件，得到相對應的地表水徑流量和地下水位降深作為輸出值，選取前40組輸入-輸出樣本作為訓練樣本，得到耦合模擬的替代模型；將后10組的輸入-輸出樣本對替代模型進行檢驗，結果見圖4。

圖4 替代模型檢驗樣本誤差Fig.4 The relative error of the test sample for the surrogate model

從圖4可以看出，1號地表水觀測點和1號地下水觀測井的相對誤差波動幅度較小，2號地下水觀測井的相對誤差波動幅度較大，認為在河流入口處地下水開采對于地表徑流的影響較小，出口處的地下水降深受地下開采的影響較為明顯。但總體變化較小，一般小于0.1%。說明替代模型的結果能夠很好地逼近地表水地下水耦合模擬模型的運行結果，且精度較高。

2.2 不確定性分析

2.2.1 統計分析

使用SPSS軟件的K-S檢驗對各地表水觀測點和地下水觀測井的數據進行統計分析，得到地表徑流量與地下水位的頻率分布直方圖，結果顯示，2號地表水觀測點的徑流量的分布較集中，服從均值為0.243，標準差為0.026的正態分布，徑流量波動范圍是0.235～0.246，地下水觀測井的水位分布不服從正態分布，認為地下水水位的分布未知。

表3 不同觀測井的統計結果Tab.3 Statistics of different observation wells

圖5 地表水觀測點徑流量的頻率分布直方圖Fig.5 Frequency distribution histogram of the runoff at the observation point

圖6 地下水觀測井水位降深的頻率分布直方圖Fig.6 Frequency distribution histogram of water level deepens in groundwater observation wells

由表3可以看出，地表觀測點徑流量的均值偏向極大值，地表水水位降深的均值偏向于極小值，說明在相同的水利條件下，地表徑流量的變化情況比地下水位降深的變化情況更為明顯；結合圖5和圖6可以看出，地下水的方差較大，數據離散化較明顯，故地下水的不確定性較大。

2.2.2 區間估計

置信區間或稱置信間距，是指在某一置信度時，總體參數所在的區域距離或區域長度[8]。置信度又稱顯著性水平，是指估計總體參數落在某一區間時，可能發生錯誤的概率。由于地表水徑流的分布已知，可以采用概率密度函數對其進行進一步的分析計算，而地下水的分布情況未知，故采用切比雪夫不等式對其分布情況進行確認和計算。

切比雪夫不等式是概率不等式的一種，利用隨機變量的數學期待和方差計算不同置信度下的區間范圍，適用于概率分布未知但方差已知的隨機變量的區間估計。

表4 不同置信水平下地表徑流量和地下水位的區間估計Tab.4 Interval estimation of surface runoff and groundwater depth at different confidence levels

在估計總體參數時，一般都會給出一個較高的置信度。本文通過替代模型獲得樣本容量為1 000的地表水徑流量和地下水水位數據，通過表4可以看出置信度越高，置信區間就越大，即估計的參數的相對精度就會越低。反之，置信度越低，則精度相對就會越高，并且區間集中于均值附近。

2.2.3 風險評估

對地下水的水位進行統計分析，并據此繪制累積概率分布圖如圖7所示。

圖7 各地下水觀測井水位累積概率分布圖Fig.7 Probability distribution of water level in groundwater observation wells

地下水觀測井沿著河流兩側自上而下分布，當1號和2號地下水觀測井同時滿足地下水水位分別低于62 m和37 m的條件時，該研究區的淺層地下水處于疏干狀態，繼續抽水會造成研究區水生態的惡化。

根據各地下水觀測井水位的累積概率分布圖可以看出，1號地下水觀測井水位低于62 m的風險是0.6，2號地下水觀測井水位低于37 m的風險為0.2，沿河分布的兩口觀測井的風險值明顯下降，地下水受到地表水的補給，水位上升，地表水于地下水之間存在水力聯系。綜合考慮研究區的用水需求，認為地下水開采造成研究區水環境失衡的風險為0.2。

2.3 運算負荷分析

研究區地表水地下水聯合模擬模型在CPU為Intel Core i5-4430，內存為4GB，操作系統為Windows 10 的計算機上，運行一次需要的時間大約為12 600 s(約3.5 h)。假設后續對地表水地下水耦合模擬模型進行不確定性分析時，需要反復調用該模型1 000次，則運行消耗的總時長約146 d。如果通過所建立的替代模型來進行計算，僅需要在前期調用50次耦合模擬模型來建立模擬模型的替代模型，耗時約7天，然后利用替代模型來進行不確定分析，耗時為5 s，總耗時為11 d。結果表明，替代模型能夠大幅度減小計算負荷和計算時間。

3 結 論

(1)采用Kriging方法建立模擬模型的替代模型，運行結果能夠充分逼近模擬模型的輸入輸出關系，能夠大幅度的減少不確定分析過程中的計算負荷。

(2)結合累積概率密度曲線對較高置信度下地下水位區間進行風險評估，發現隨著地表徑流和降水入滲的補給，河口地表水水位抬升，對地下水進行補給，地表水與地下水之間存在水力聯系和相互轉換作用。

(3)在模擬模型替代模型的基礎上，對地表水地下水耦合模擬模型的源匯項進行不確定性分析，可以定性地對模型輸出結果的可靠性進行分析，為決策者提供參考，對于區域水資源開發利用提供科學依據。