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

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

(1.吉林大學(xué) 地下水與資源環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,吉林 長春 130012；2.吉林大學(xué) 環(huán)境與資源學(xué)院,吉林 長春 130012)

0 引 言

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

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

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

1 研究方法

1.1 模型建立

建立一個(gè)三維的非穩(wěn)定潛水含水層的假想研究區(qū)，研究區(qū)兩側(cè)高，中間低，河流自南側(cè)流入，北側(cè)流出。

1.1.1 概念模型

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

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

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

表1 研究區(qū)地下水水文地質(zhì)參數(shù)值Tab.1 Groundwater Hydrogeology parameter initial values of study area

(3)源匯項(xiàng)。研究區(qū)接受大氣降水補(bǔ)給，研究區(qū)的主要排泄方式為蒸發(fā)排泄和人工開采?，F(xiàn)在研究區(qū)內(nèi)設(shè)置3口抽水井，2個(gè)地表水地表水觀測(cè)點(diǎn)和21口地下水觀測(cè)井，對(duì)地表徑流量與地下水水位進(jìn)行數(shù)據(jù)的采集，各井的分布情況如圖1所示。模擬時(shí)長為3個(gè)月，為降水較為充沛的5，6，7月，本次研究對(duì)模型運(yùn)行3個(gè)月后的結(jié)果進(jìn)行分析計(jì)算。

表2 研究區(qū)地表水力參數(shù)值Tab.2 Initial values of surface hydraulic parameter

1.1.2 數(shù)學(xué)模型

采用三維Richard方程和擴(kuò)散波近似的二維平均深度圣維南(Saint Venant)方程描述地下水水流和地表水水流運(yùn)動(dòng)，采用雙重節(jié)點(diǎn)法對(duì)地表水與地下水模型進(jìn)行耦合，通過達(dá)西流關(guān)系描述地表水系統(tǒng)與地下水系統(tǒng)的交換量，得到如下微分方程，定解條件和耦合方程：

圖2 地表水與地下水的參數(shù)分區(qū)圖Fig.2 Distribution map of hydrogeology parameter of groundwater and surface water

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

(1)

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

1.2 Kriging方法

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

Kriging模型的數(shù)學(xué)表達(dá)式是：

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

(2)

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

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

(3)

式中：R(xi,xj)為空間采樣點(diǎn)xi,xj的相關(guān)函數(shù)，可以采用以下函數(shù)形式進(jìn)行表示：

(4)

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

(5)

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

1.3 Monte Carlo模擬

Monte Carlo方法又稱為統(tǒng)計(jì)模擬方法，是一種常用的數(shù)值模擬的不確定性分析方法，在解決非線性的隨機(jī)問題的過程中具有一定的優(yōu)勢(shì)[7]。

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

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

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

2 結(jié)果與討論

2.1 Kriging替代模型

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

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

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

2.2 不確定性分析

2.2.1 統(tǒng)計(jì)分析

使用SPSS軟件的K-S檢驗(yàn)對(duì)各地表水觀測(cè)點(diǎn)和地下水觀測(cè)井的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，得到地表徑流量與地下水位的頻率分布直方圖，結(jié)果顯示，2號(hào)地表水觀測(cè)點(diǎn)的徑流量的分布較集中，服從均值為0.243，標(biāo)準(zhǔn)差為0.026的正態(tài)分布，徑流量波動(dòng)范圍是0.235～0.246，地下水觀測(cè)井的水位分布不服從正態(tài)分布，認(rèn)為地下水水位的分布未知。

表3 不同觀測(cè)井的統(tǒng)計(jì)結(jié)果Tab.3 Statistics of different observation wells

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

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

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

2.2.2 區(qū)間估計(jì)

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

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

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

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

2.2.3 風(fēng)險(xiǎn)評(píng)估

對(duì)地下水的水位進(jìn)行統(tǒng)計(jì)分析，并據(jù)此繪制累積概率分布圖如圖7所示。

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

地下水觀測(cè)井沿著河流兩側(cè)自上而下分布，當(dāng)1號(hào)和2號(hào)地下水觀測(cè)井同時(shí)滿足地下水水位分別低于62 m和37 m的條件時(shí)，該研究區(qū)的淺層地下水處于疏干狀態(tài)，繼續(xù)抽水會(huì)造成研究區(qū)水生態(tài)的惡化。

根據(jù)各地下水觀測(cè)井水位的累積概率分布圖可以看出，1號(hào)地下水觀測(cè)井水位低于62 m的風(fēng)險(xiǎn)是0.6，2號(hào)地下水觀測(cè)井水位低于37 m的風(fēng)險(xiǎn)為0.2，沿河分布的兩口觀測(cè)井的風(fēng)險(xiǎn)值明顯下降，地下水受到地表水的補(bǔ)給，水位上升，地表水于地下水之間存在水力聯(lián)系。綜合考慮研究區(qū)的用水需求，認(rèn)為地下水開采造成研究區(qū)水環(huán)境失衡的風(fēng)險(xiǎn)為0.2。

2.3 運(yùn)算負(fù)荷分析

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

3 結(jié) 論

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

(2)結(jié)合累積概率密度曲線對(duì)較高置信度下地下水位區(qū)間進(jìn)行風(fēng)險(xiǎn)評(píng)估，發(fā)現(xiàn)隨著地表徑流和降水入滲的補(bǔ)給，河口地表水水位抬升，對(duì)地下水進(jìn)行補(bǔ)給，地表水與地下水之間存在水力聯(lián)系和相互轉(zhuǎn)換作用。

(3)在模擬模型替代模型的基礎(chǔ)上，對(duì)地表水地下水耦合模擬模型的源匯項(xiàng)進(jìn)行不確定性分析，可以定性地對(duì)模型輸出結(jié)果的可靠性進(jìn)行分析，為決策者提供參考，對(duì)于區(qū)域水資源開發(fā)利用提供科學(xué)依據(jù)。