基于Monte Carlo方法的地表水地下水耦合模擬模型不確定分析

張將偉，盧文喜，曲延光 ，安永凱

（1.吉林大學(xué) 地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，吉林 長(zhǎng)春 130021；2.吉林大學(xué) 環(huán)境與資源學(xué)院，吉林 長(zhǎng)春 130021；3.吉林省水文水資源局，吉林 長(zhǎng)春 130021）

1 研究背景

地表水與地下水之間的相互轉(zhuǎn)換是水循環(huán)過程中的重要特征［1］，運(yùn)用數(shù)學(xué)模型對(duì)地表水地下水進(jìn)行耦合模擬是研究地表水地下水相互作用的常用手段。由于模型輸入、模型結(jié)構(gòu)以及觀測(cè)數(shù)據(jù)等多個(gè)方面存在不確定性，耦合模型的計(jì)算結(jié)果與真實(shí)值之間往往存在一定偏差，這是耦合模型中不確定性的實(shí)際體現(xiàn)［2-5］。單純運(yùn)用確定性模型研究地表水與地下水的運(yùn)動(dòng)與轉(zhuǎn)化規(guī)律，只能得到唯一結(jié)果，難以確定其可靠程度。對(duì)地表水地下水耦合模擬模型進(jìn)行不確定性分析，可以分析模型的可靠程度，滿足管理者對(duì)于模型極端情況預(yù)測(cè)的需求［6］，對(duì)地表水地下水聯(lián)合調(diào)度有重要意義。

模型參數(shù)的不確定性是模型不確定性研究的重要內(nèi)容之一［7］。目前，國(guó)內(nèi)外學(xué)者在數(shù)學(xué)模擬模型的參數(shù)不確定性分析方面開展了大量工作。在水文模型方面，Beven等人提出水文模型的異參同效性，并發(fā)展了GLUE方法（generalized likelihood uncertainty estimation），對(duì)水文模型的參數(shù)不確定性進(jìn)行估計(jì)［8-10］。在地下水模型方面，按照求解原理可以將不確定性分析方法分為：蒙特卡羅法、矩方程法和貝葉斯法三大類以及其它方法（如條件模擬、敏感度分析、一次二階矩法等）［11］。Hassan等［12］應(yīng)用馬爾科夫鏈-蒙特卡洛方法（Markov Chain Monte Carlo，MCMC）對(duì)Alaska，Amchitka Island的Milrow試驗(yàn)場(chǎng)的地下水模型參數(shù)進(jìn)行了不確定性評(píng)價(jià)，證明了MCMC方法能夠有效的探索參數(shù)分布空間的高概率密度區(qū)域，并將其反饋到參數(shù)后驗(yàn)概率分布的構(gòu)造過程；姚磊華等［13］在矩方程法基礎(chǔ)上，將Taylor展開、攝動(dòng)技術(shù)、待定系數(shù)法和有限元法相結(jié)合，提出了待定系數(shù)攝動(dòng)隨機(jī)有限元法，推導(dǎo)出地下水水頭均值和方差的表達(dá)式，并選取二維承壓地下水水流進(jìn)行隨機(jī)模擬；梁婕等［14］應(yīng)用貝葉斯推斷方法和隨機(jī)模擬技術(shù)，定量研究了滲透系數(shù)的非均質(zhì)性對(duì)地下水溶質(zhì)運(yùn)移的影響，并進(jìn)行了二維理想算例的分析計(jì)算。但以往的研究只局限于地表水或地下水，針對(duì)流域尺度的地表水地下水耦合模型進(jìn)行不確定分析還是有待解決的問題。

本文以石頭口門水庫(kù)上游飲馬河匯水流域?yàn)槔\(yùn)用Monte Carlo方法對(duì)地表水地下水耦合模擬模型進(jìn)行不確定性分析。首先，根據(jù)研究區(qū)地表水地下水的實(shí)際條件建立地表水地下水耦合模擬模型，用以描述地表水地下水各自的運(yùn)動(dòng)特征及二者之間的水力聯(lián)系和相互轉(zhuǎn)換作用。運(yùn)用HGS軟件對(duì)耦合模型進(jìn)行求解，并根據(jù)實(shí)測(cè)數(shù)據(jù)對(duì)耦合模型進(jìn)行校正和檢驗(yàn)。應(yīng)用局部靈敏度分析方法對(duì)地表水地下水耦合模擬模型參數(shù)進(jìn)行靈敏度分析，篩選出對(duì)模型輸出結(jié)果影響顯著的參數(shù)。將靈敏度較高的參數(shù)作為隨機(jī)變量，進(jìn)行隨機(jī)試驗(yàn)，對(duì)地表水地下水耦合模擬模型進(jìn)行不確定性分析，考察參數(shù)的不確定性對(duì)模型輸出結(jié)果的影響。為減小不確定性分析過程中的計(jì)算負(fù)荷，應(yīng)用Kriging方法，建立了地表水地下水耦合模擬模型的替代模型。最后，根據(jù)不確定性分析結(jié)果，對(duì)研究區(qū)地表水和地下水生態(tài)惡化風(fēng)險(xiǎn)進(jìn)行風(fēng)險(xiǎn)評(píng)估。

2 研究區(qū)概況

研究區(qū)位于吉林省中部長(zhǎng)春市西側(cè)，為石頭口門上游飲馬河匯水流域，總面積約為1197 km2，如圖1所示。研究區(qū)地勢(shì)總體上呈南高北低，全區(qū)的平均海拔為418 m。石頭口門上游飲馬河匯水流域?qū)儆跍貛Т箨懶约撅L(fēng)性氣候類型，四季氣候變化明顯，多年平均降水量642.9 mm，多年平均蒸發(fā)量1339 mm。研究區(qū)內(nèi)主要河流為飲馬河，由南至北流動(dòng)，最終匯入石頭口門水庫(kù)。該河段河道較彎曲，河槽寬淺，兩岸為平緩的低丘和臺(tái)地，中下游多為洼地和沼澤，河岸兩側(cè)多為砂土和砂壤土，河床內(nèi)沉積物為砂或者細(xì)砂，多年平均徑流量約為3 m3/s。

研究區(qū)北部為山前平原區(qū)，地層巖性以第四系沖積洪積層粗砂和砂礫石為主。南部和西北部為低山丘陵區(qū)，地層巖性主要為二疊系拉溪組的淺變質(zhì)巖和侏羅系上統(tǒng)沙河子組和營(yíng)城組碎屑巖。中部為波狀臺(tái)地與丘陵?duì)钆_(tái)地區(qū)，地層巖性主要為為下更新統(tǒng)白土山組砂礫石、中更新統(tǒng)黃土狀亞黏土。區(qū)內(nèi)人工開采地下水主要用于生活用水和農(nóng)田灌溉，開采量約為1.3億m3/a。

圖1 研究區(qū)概況圖

3 研究方法

3.1 HydrogeosphereHydroGeoSphere（HGS）是一款由Aquanty公司開發(fā)的地表水地下水耦合模擬軟件。它從水文循環(huán)物理過程出發(fā)，充分考慮地表水地下水之間水力聯(lián)系和相互轉(zhuǎn)化過程，是一款能夠在流域尺度下進(jìn)行地表水地下水耦合模擬的軟件。

根據(jù)研究區(qū)實(shí)際情況對(duì)研究區(qū)地表水系統(tǒng)和地下水系統(tǒng)進(jìn)行概化。將研究區(qū)地表水系統(tǒng)的上游入流斷面概化為已知流量邊界，研究區(qū)下游為石頭口門水庫(kù)，故將下游斷面概化為臨界深度邊界，將匯水流域分水嶺概化為零深坡度邊界。區(qū)內(nèi)源項(xiàng)主要包括大氣降水，匯項(xiàng)主要包括蒸發(fā)蒸騰。不同的土地利用類型導(dǎo)致不同的下墊面條件，因此，根據(jù)不同土地利用類型將研究區(qū)分為耕地、林地、草地、灘地、居民點(diǎn)和城鎮(zhèn)六個(gè)類別，賦予不同的地表水水力參數(shù)，如圖2所示。

圖2 研究區(qū)地表水邊界與土地利用類型

圖3 研究區(qū)地下水邊界與參數(shù)分區(qū)

研究區(qū)地下水系統(tǒng)頂部邊界為地表面，與農(nóng)業(yè)灌溉回歸水和城鎮(zhèn)、農(nóng)村生活用水回滲等存在水力聯(lián)系（大氣降水補(bǔ)給僅在地表水模型中給予考慮，在地下水模型中不再重復(fù)計(jì)算。）。研究區(qū)底部為透水性較差的淤泥質(zhì)亞黏土，作為隔水邊界。研究區(qū)側(cè)向邊界如圖3所示。研究區(qū)北部為山前平原區(qū)，地層巖性主要粗砂和砂礫石，大部分區(qū)域介質(zhì)透水性較好，富水性較強(qiáng)；南部低山丘陵區(qū)，地層巖性主要為碎屑巖、變質(zhì)巖和花崗巖，發(fā)育在其中的網(wǎng)狀裂隙部分被填充，連通性差，地下水富水性較弱；中部為波狀臺(tái)地與丘陵?duì)钆_(tái)地區(qū)，地層巖性主要為砂礫石和黃土狀亞黏土，地下水富水性較弱。

根據(jù)概化結(jié)果建立數(shù)學(xué)模擬模型，見式（1），使用二維Saint Venant方程來描述地表水運(yùn)動(dòng)，使用三維Richard方程描述地下水水流運(yùn)動(dòng)，使用達(dá)西定律描述地表水與地下水交換過程，并在達(dá)西定律的基礎(chǔ)上引入相對(duì)滲透率，反映包氣帶對(duì)地表水地下水水量交換的影響。

式中：do為地表徑流深度，m；qo為地表水流速，m/s；qgs為地下水和地表水之間單位體積流量，s-1；ho為地表水水位，m；Q為地表水的源匯項(xiàng)，為單位面積體積流量，m/s；φo為等效地表孔隙度；h（x，y，0）為（x，y）點(diǎn)初始水位，m；Γ1為上游段面，Γ2為分水嶺，Γ3為下游斷面；H為地下水水頭，m；W為地下水的源匯項(xiàng)，為單位體積流量，s-1；φ為孔隙度，無量綱；Sw為飽和度，無量綱；Ss為儲(chǔ)水系數(shù)，m-1；qG為邊界上的水流通量，m/s；H0為地下水初始水頭，m；S1已知流量邊界，S2為隔水邊界；Γgs為地表水與地下水之間交換通量，s-1；Kzz為含水介質(zhì)滲透系數(shù)，m/s；krw為相對(duì)滲透率，無量綱；lexch為地表水與地下水之間的耦合長(zhǎng)度，m。

通過Van Genuchten關(guān)系來描述相對(duì)滲透率（krw）與飽和度（Sw）和壓力水頭（ψ）之間的關(guān)系，其關(guān)系式如下：

式中：lp為孔隙連接指數(shù)；γ為經(jīng)驗(yàn)參數(shù)；Se為有效含水量；Swr為殘余含水率；α、β為VG方程參數(shù)。

在平面上，采用三角網(wǎng)格剖分法對(duì)模型進(jìn)行剖分，共得到節(jié)點(diǎn)5539個(gè)，三角形單元10448個(gè)。在垂向上，根據(jù)地層巖性，將含水層剖分為10層，弱透水層剖分為1層。同時(shí)，為研究地表水與地下水的相互作用，在地下水地表水交界面附近進(jìn)行加密剖分，模型空間離散圖如圖4所示。

圖4 模型空間離散示意圖

基于研究區(qū)已有的水文及水文地質(zhì)資料，選取2016年7月25日作為初始時(shí)刻，使用觀測(cè)資料插值獲得研究區(qū)地下水初始水位和地表水初始水深，非飽和帶水分分布可由式（2）、式（3）計(jì)算得出。將2016年7月25日—8月25日作為模型的校正時(shí)段，校正期內(nèi)逐月觀測(cè)數(shù)據(jù)作為校正依據(jù)，采用試估校正法來校正模型各參數(shù)和源匯項(xiàng)。在此基礎(chǔ)上將2016年8月26日—9月26日作為模型的檢驗(yàn)時(shí)段，檢驗(yàn)期內(nèi)逐月觀測(cè)數(shù)據(jù)作為檢驗(yàn)依據(jù)。

3.2 靈敏度分析靈敏度分析可以定量地評(píng)價(jià)輸入?yún)?shù)變化對(duì)輸出結(jié)果的影響，通過確定模型中靈敏度較高參數(shù)可以減少M(fèi)onte Carlo模擬所考察的參數(shù)個(gè)數(shù)，從而簡(jiǎn)單、高效地對(duì)模型進(jìn)行不確定性分析。按分析方法分類，靈敏度可以分為局部靈敏度分析和全局靈敏度分析［15-17］。由于全局靈敏度分析時(shí)需要多次調(diào)用模擬模型，耗費(fèi)時(shí)間較長(zhǎng)，尤其是本次研究過程中所建立的地表水地下水耦合模擬模型采用牛頓迭代方法求解，在保證研究結(jié)果精度的前提下，運(yùn)行時(shí)間較長(zhǎng)。因此本文采用局部靈敏度分析方法對(duì)模擬模型參數(shù)進(jìn)行靈敏度分析。

從數(shù)學(xué)意義上來講，參數(shù)靈敏度可理解為變量x的微小變化所引起的函數(shù)F（x）的變化程度。局部靈敏度計(jì)算公式如式（4）所示：

式中：S為參數(shù)靈敏度；xi代表不同的模型參數(shù)；Δxi代表對(duì)參數(shù)xi的微小擾動(dòng)。

3.3 Monte Carlo方法Monte Carlo方法是一種隨機(jī)分析方法，它是一種通過隨機(jī)實(shí)驗(yàn)來求解隨機(jī)問題的技術(shù)，即采用統(tǒng)計(jì)抽樣理論近似求解數(shù)學(xué)問題或物理問題的一種方法，是進(jìn)行不確定性分析的一種有效方法［18-21］。

蒙特卡羅模擬方法假定隨機(jī)變量的概率分布函數(shù)已知，通過隨機(jī)抽樣方法獲得多組輸入變量，每組變量都相當(dāng)于一次統(tǒng)計(jì)試驗(yàn)，然后將隨機(jī)變量帶入模型中從而得到大量的輸出，根據(jù)對(duì)輸出結(jié)果的統(tǒng)計(jì)可以得到均值、方差等統(tǒng)計(jì)估計(jì)量，并擬合其輸出結(jié)果的概率分布情況。

它的主要思路如下：（1）在各隨機(jī)變量的可行域內(nèi)根據(jù)其服從的概率分布進(jìn)行隨機(jī)抽樣。在抽樣過程中，可以采用拉丁超立方抽樣法，使得樣本的覆蓋率更高，減少抽樣次數(shù)。（2）將隨機(jī)變量的抽樣結(jié)果進(jìn)行隨機(jī)組合后，通過運(yùn)轉(zhuǎn)模型獲得相應(yīng)的輸出結(jié)果。（3）統(tǒng)計(jì)分析模型輸出結(jié)果，根據(jù)統(tǒng)計(jì)結(jié)果的均值、方差等估計(jì)量，擬合輸出結(jié)果的概率分布情況，定量描述模型的不確定性。

為了減小Monte Carlo方法應(yīng)用過程中反復(fù)調(diào)用耦合模擬模型所導(dǎo)致的大量計(jì)算負(fù)荷，本文應(yīng)用Kriging方法建立耦合模擬模型的替代模型。

3.4 Kriging方法Kriging方法是建立在變異函數(shù)理論分析基礎(chǔ)上，對(duì)有限區(qū)域內(nèi)的區(qū)域化變量取值進(jìn)行無偏最優(yōu)估計(jì)的一種空間局部?jī)?nèi)插法。近年來，被廣泛應(yīng)用于建立替代模型［22-25］。Kriging模型表達(dá)式如下：

式中：f（x）=［f1（x），f2（x），…，fk（x）］為已知回歸模型的基函數(shù)；β=［β1，β2，…，βk］為待定參數(shù)；z（x）為隨機(jī)部分，均值為0，方差為σ2，協(xié)方差為：

式中R（xi，xj）為點(diǎn)和點(diǎn)的關(guān)聯(lián)系數(shù)，可用高斯模型進(jìn)行計(jì)算：

式中：n為設(shè)計(jì)變量個(gè)數(shù)；θk和pk為待定系數(shù)。則有已知在點(diǎn)x處的y（x）估計(jì)值為：

式中：rT（x）為點(diǎn)x與m個(gè)采樣點(diǎn)之間的相關(guān)向量；y為m×1向量，包含m個(gè)采樣點(diǎn)對(duì)應(yīng)的響應(yīng)值；β可由最優(yōu)線性無偏估計(jì)得到：

式中R為n個(gè)采樣點(diǎn)相關(guān)系數(shù)組成的相關(guān)矩陣：

方差σ2的估計(jì)值可用σ2=（y-fβ）TR-1（y-fβ），參數(shù)θk則根據(jù)極大似然估計(jì)給出［26］。

4 結(jié)果與討論

4.1 模型的校正與檢驗(yàn)結(jié)果地表水地下水耦合模擬模型校正和檢驗(yàn)?zāi)┛蹋饔^測(cè)井地下水水頭的實(shí)測(cè)值與計(jì)算值如圖5所示，各測(cè)站地表徑流量的實(shí)測(cè)值與計(jì)算值如表1所示。

圖5 地下水觀測(cè)井實(shí)測(cè)水位與計(jì)算水位

表1 地表水測(cè)站實(shí)測(cè)流量與計(jì)算流量

由圖5可知，校正和檢驗(yàn)階段地表水地下水耦合模擬模型中地下水水位的計(jì)算值、實(shí)測(cè)值與y=x間的擬合度分別為0.9997和0.9993。

使用相對(duì)誤差和NSE（納什效率系數(shù)）誤差評(píng)價(jià)地表水流量的校正檢驗(yàn)的精度。NSE誤差計(jì)算方法如下：

式中：T為地表水測(cè)站數(shù)量；Qo為地表水流量的實(shí)測(cè)值；Qm為計(jì)算值。

由表1可知，各測(cè)站地表徑流量的計(jì)算值與實(shí)測(cè)值的相對(duì)誤差均小于20%，且NSE誤差均大于0。

綜上所述，所建立的地表水地下水耦合模擬模型能夠反映研究區(qū)的地下水水位變化和地表水流量變化的實(shí)際情況。研究區(qū)水文地質(zhì)參數(shù)和地表水力參數(shù)的取值如表2、表3所示。

4.2 靈敏度分析結(jié)果根據(jù)研究區(qū)水文地質(zhì)參數(shù)分區(qū)，選取井5、井13、井14（圖1）作為不同水文地質(zhì)區(qū)域的典型研究對(duì)象，地表水測(cè)站選取小金屯測(cè)站為研究對(duì)象。

對(duì)地表水地下水耦合模擬模型中的滲透系數(shù)、孔隙度、彈性釋水系數(shù)、曼寧系數(shù)洼地儲(chǔ)存量和耦合長(zhǎng)度進(jìn)行局部靈敏度分析，分析結(jié)果如圖6所示。由圖6可知對(duì)耦合模擬模型中地下水水位影響較大的2個(gè)隨機(jī)參數(shù)分別是滲透系數(shù)和孔隙度，對(duì)耦合模擬模型中地表徑流量影響較大的2個(gè)隨機(jī)參數(shù)分別是滲透系數(shù)和曼寧粗糙系數(shù)，滲透系數(shù)的變化不僅對(duì)會(huì)耦合模型中地下水水位影響較大，也對(duì)地表徑流量影響較大。

表2 研究區(qū)含水層水文地質(zhì)參數(shù)校正結(jié)果

表3 研究區(qū)地表水水力參數(shù)校正結(jié)果

圖6 典型觀測(cè)井與測(cè)站靈敏度分析結(jié)果

故將滲透系數(shù)、孔隙度和曼寧粗糙系數(shù)作為隨機(jī)變量，其他參數(shù)取校正值作為確定值輸入模型。4.3 Kriging模型將預(yù)報(bào)期設(shè)為2016年10月31日，利用地表水地下水耦合模擬模型進(jìn)行預(yù)報(bào)。根據(jù)靈敏度分析結(jié)果，將曼寧粗糙系數(shù)、孔隙度和滲透系數(shù)作為隨機(jī)變量，其取值范圍為校正值±20%，其他參數(shù)及源匯項(xiàng)設(shè)定為校正值。根據(jù)前人經(jīng)驗(yàn)，曼寧粗糙系數(shù)、孔隙度和滲透系數(shù)服從的參數(shù)分布如表4所示。

運(yùn)用拉丁超立方方法抽樣100組隨機(jī)變量，作為100組預(yù)報(bào)方案。統(tǒng)計(jì)100組預(yù)報(bào)方案下地表水地下水耦合模擬模型中井5、井13、井14與小金屯測(cè)站的預(yù)報(bào)結(jié)果，以此為基礎(chǔ)建立Kriging替代模型代替耦合模擬模型。

運(yùn)用MATLAB軟件建立替代模型。將耦合模擬模型的100組預(yù)報(bào)結(jié)果中前90組作為替代模型的訓(xùn)練樣本來訓(xùn)練替代模型。后10組作為替代模型的檢驗(yàn)樣本，來檢驗(yàn)替代模型的擬合精度，檢驗(yàn)結(jié)果見圖7。

由圖7可知，替代模型的地下水水頭輸出結(jié)果與耦合模擬模型地下水水頭的變化趨勢(shì)基本一致，且相對(duì)誤差均小于0.2%；替代模型的地表徑流量輸出結(jié)果與耦合模擬模型地表徑流量的變化趨勢(shì)基本一致，相對(duì)誤差小于4%。擬合相對(duì)誤差較小，說明替代模型精度較高，可以代替耦合模擬模型求解。

4.4 不確定性分析再次利用拉丁超立方方法抽取1000組隨機(jī)變量樣本，作為輸入數(shù)據(jù)集。將輸入數(shù)據(jù)集輸入Kriging替代模型，求解得到1000組井5、井13、井14和小金屯測(cè)站的輸出值。

對(duì)所獲取的1000組輸出值進(jìn)行統(tǒng)計(jì)分析，統(tǒng)計(jì)結(jié)果如表5所示。由表5可知井14、井5、井13標(biāo)準(zhǔn)差依次變大，說明輸出結(jié)果越來越分散，不確定性越來越大。這可能與含水介質(zhì)不同有關(guān)，井5、井13、井14分別為碎屑巖-花崗巖含水層、黃土狀亞黏土含水層和第四系沖積洪積物含水層地區(qū)典型觀測(cè)井，故認(rèn)為耦合模擬模型中亞黏土含水層地區(qū)輸出結(jié)果的不確定性大于碎屑巖-花崗巖含水層和第四系沖積洪積物含水層。

表4 參數(shù)的概率分布與取值

圖7 替代模型精度檢驗(yàn)

假設(shè)井5、井13、井14和小金屯測(cè)站的輸出數(shù)值均符合正態(tài)分布，利用SPSS軟件中的S-K檢驗(yàn)分別對(duì)各組輸出數(shù)值進(jìn)行假設(shè)檢驗(yàn)。根據(jù)前人經(jīng)驗(yàn)認(rèn)為當(dāng)雙側(cè)漸進(jìn)顯著性大于等于0.05時(shí)假設(shè)成立，故小金屯測(cè)站輸出數(shù)值符合均值為11.21、標(biāo)準(zhǔn)差為0.68的正態(tài)分布，井5輸出數(shù)值符合均值為230.77、標(biāo)準(zhǔn)差為0.09的正態(tài)分布。井13、井14輸出數(shù)值分布未知。

表5 輸出數(shù)據(jù)單樣本K-S分析結(jié)果

4.5 區(qū)間估計(jì)區(qū)間估計(jì)指在概率學(xué)中一定置信度下對(duì)未知參數(shù)真值存在范圍的估計(jì)。本文將井5、井13、井14和小金屯測(cè)站的輸出值作為未知參數(shù)，計(jì)算其在不同置信水平下的區(qū)間范圍。

切比雪夫不等式適用于分布未知，而已知平均值和方差的情況，故對(duì)井13、井14運(yùn)用切比雪夫不等式進(jìn)行區(qū)間估計(jì)。井5、小金屯測(cè)站的輸出數(shù)值均符合正態(tài)分布，運(yùn)用概率密度函數(shù)計(jì)算。井5、井13、井14和小金屯測(cè)站區(qū)間估計(jì)結(jié)果如表6所示，由表可以看出，置信水平越高，區(qū)間范圍越大；置信水平越低，區(qū)間范圍越小，也越集中在均值附近。

表6 不同置信水平下典型觀測(cè)井與測(cè)站的區(qū)間估計(jì)

4.6 生態(tài)風(fēng)險(xiǎn)評(píng)估生態(tài)風(fēng)險(xiǎn)評(píng)估是描述和評(píng)估人為活動(dòng)、環(huán)境污染和自然災(zāi)害對(duì)生態(tài)系統(tǒng)結(jié)構(gòu)與功能產(chǎn)生負(fù)面作用的可能性及其大小的過程［27］。通常使用指數(shù)法，模擬法和統(tǒng)計(jì)法進(jìn)行生態(tài)風(fēng)險(xiǎn)評(píng)估。本文采用模擬法，根據(jù)統(tǒng)計(jì)不確定分析過程中獲取的1000組輸出結(jié)果和地表水流量與地下水水位的生態(tài)閾值，對(duì)地表水流量與地下水水位變化超出生態(tài)閾限的狀況進(jìn)行了風(fēng)險(xiǎn)評(píng)估。

根據(jù)經(jīng)驗(yàn)，假設(shè)預(yù)報(bào)期末刻，井5、井13、井14水位同時(shí)低于230.65、221.80、184.8 m時(shí)存在地下水生態(tài)環(huán)境惡化風(fēng)險(xiǎn)；小金屯水文站流量小于10.5 m3/s會(huì)引起地表水水生態(tài)環(huán)境惡化。根據(jù)各個(gè)觀測(cè)井、水文測(cè)站輸出值累積分布函數(shù)畫出風(fēng)險(xiǎn)評(píng)估圖（圖8）可以看出，在當(dāng)前開采狀態(tài)下，井5、井13、井14水位低于230.65、221.80、184.8 m分別為11%、6%和12%，故存在地下水生態(tài)環(huán)境惡化風(fēng)險(xiǎn)概率為6%。同理，在當(dāng)前開采狀態(tài)下，地表水生態(tài)環(huán)境惡化風(fēng)險(xiǎn)為15%。

圖8 各觀測(cè)井和測(cè)站風(fēng)險(xiǎn)評(píng)估圖

4.7 計(jì)算負(fù)荷分析在CPU為Intel Xeon E5 2.4G Hz，內(nèi)存為64 GB，操作系統(tǒng)為Windows 7的計(jì)算機(jī)上，運(yùn)行一次石頭口門水庫(kù)上游飲馬河匯水流域地表水地下水耦合模擬模型需要的時(shí)間約為加粗4.7 s（5.5 h），而運(yùn)行一次替代模型需要的時(shí)間約為10 s。說明建立替代模型可以有效減小Monte Carlo過程中反復(fù)調(diào)用耦合模擬模型所產(chǎn)生的計(jì)算負(fù)荷。

5 結(jié)論

本文以石頭口門水庫(kù)上游飲馬河匯水流域?yàn)檠芯繀^(qū)，建立研究區(qū)地表水地下水水流耦合數(shù)學(xué)模擬模型。通過局部靈敏度分析法獲取了3個(gè)敏感性較大參數(shù)作為隨機(jī)變量。運(yùn)用Monte Carlo方法對(duì)耦合模擬模型進(jìn)行不確定性分析。

靈敏度分析結(jié)果表明，滲透系數(shù)，孔隙度和曼寧粗糙系數(shù)對(duì)地表水地下水耦合模擬模型的輸出結(jié)果影響較大。其中，滲透系數(shù)的變化不僅會(huì)對(duì)耦合模型中地下水輸出結(jié)果產(chǎn)生影響，也會(huì)通過改變地表水地下水之間的交換量，對(duì)地表水的輸出結(jié)果產(chǎn)生影響。

利用Kriging方法建立的替代模型能夠充分逼近地表水地下水耦合模擬模型的輸入響應(yīng)關(guān)系。在不確定性分析過程中運(yùn)用替代模型能夠大幅度反復(fù)調(diào)用耦合模擬模型所產(chǎn)生的計(jì)算負(fù)荷。

風(fēng)險(xiǎn)評(píng)估結(jié)果表明，在當(dāng)前開采狀態(tài)下，研究區(qū)地下水生態(tài)環(huán)境惡化風(fēng)險(xiǎn)概率為6%，地表水生態(tài)環(huán)境惡化風(fēng)險(xiǎn)為15%。

運(yùn)用區(qū)間估計(jì)方法可以計(jì)算不同置信水平下地表水流量和地下水水位的取值區(qū)間，相較于確定性模型所得到的唯一結(jié)果更有實(shí)際指導(dǎo)意義。