張 飛,溫忠輝,戴鳳君,李　偉

(1.河海大學水文水資源學院,江蘇 南京　210098;2.江蘇水文水資源勘測局鎮(zhèn)江分局,江蘇 鎮(zhèn)江　212003;3. 江蘇省地質調查研究院,江蘇 南京　210018)

基于滲透系數(shù)序貫高斯模擬的水庫滲漏量不確定性分析

張 飛1,溫忠輝1,戴鳳君2,李偉3

(1.河海大學水文水資源學院,江蘇 南京210098;2.江蘇水文水資源勘測局鎮(zhèn)江分局,江蘇 鎮(zhèn)江212003;3. 江蘇省地質調查研究院,江蘇 南京210018)

摘要：水庫庫底地層滲透系數(shù)是影響水庫滲漏量的主要因素之一,研究滲透系數(shù)的空間變異特征可為水庫滲漏量計算結果的可靠性分析提供科學依據(jù)。以天津市北塘水庫庫底地層滲透系數(shù)為研究對象,先利用克里金插值對該參數(shù)系列進行了插補,再對插補后的滲透系數(shù)系列進行了500次序貫高斯模擬,最后分析了其空間變異性對庫區(qū)滲漏量計算值不確定性的影響。結果表明:庫底地層垂向滲透系數(shù)空間結構可用球狀模型描述,在變程1 210 m范圍內具有高度空間相關性。以2005年7月16日滲漏量計算值總體為例,最小值與最大值分別為1 389.49 m3/d、1 897.30 m3/d,且落在區(qū)間1 600～1 850 m3/d的頻率高達77.6%。

關鍵詞：滲透系數(shù);條件模擬;空間變異性;水庫滲漏量;不確定性

地表水和地下水的相互作用是自然界中普遍存在的一種現(xiàn)象。正確認識地表水與地下水之間的相互作用,尤其是水量和水質的時空變化規(guī)律研究及交換水量的準確計算,對水資源的綜合管理及合理開發(fā)利用具有重要的科學意義[1]。進入21世紀,國內外許多組織和機構都將兩者的相互作用作為水文學及水文地質學方面研究的熱點問題[2]。

水文地質參數(shù)的空間變異性是影響地表水與地下水相互作用強度的主要因素之一。其中,滲透系數(shù)在野外實際條件下的空間變異性較其他水文地質參數(shù)要大得多,最大值與最小值相差可達13個數(shù)量級[3]。這就決定了滲透系數(shù)必須作為隨機變量處理。同時滲透系數(shù)又具有一定程度上的結構性,比如,含水介質是在特定地質條件下沉積形成的,具有一定成層性[4]。因此,地表水與地下水之間交換量計算的可靠度在很大程度上依賴于參數(shù)隨機場的模擬是否反映出實際地區(qū)參數(shù)的空間變異性。

胡伏生等[5]利用條件模擬對內蒙古霍林河市主要含水層滲透系數(shù)的分布特征進行了分析;陳彥等[4]采用Monte Carlo法來研究含水介質滲透系數(shù)空間變異性對地下水數(shù)值模擬的影響;彭偉等[6]進行了滲透系數(shù)的序貫高斯模擬對污染物運移的不確定性研究;束龍倉等[7]在對北塘水庫庫底地層滲透系數(shù)進行Monte Carlo模擬的基礎上,進行了地表水與地下水交換量的不確定性研究。綜合前人研究的成果,與Monte Carlo法相比,條件模擬既考慮了參數(shù)的隨機性,又考慮了參數(shù)的結構性,更能客觀地反映參數(shù)實際變化特征[6]。筆者在文獻[7]研究的基礎上,運用地質統(tǒng)計學的原理和方法,分析參數(shù)的隨機性和空間相關性,進一步探討參數(shù)條件模擬對交換量不確定性的影響。

1序貫高斯模擬法

對于某一特定的含水層而言,其空間上每一點的滲透系數(shù)都是確定的。但由于滲透系數(shù)難于準確獲取,并且獲取所有點的滲透系數(shù)也不太現(xiàn)實,因此如何合理利用有限個參數(shù)樣本來推斷含水層參數(shù)在空間上的分布特征成為難點之一。條件模擬是地質統(tǒng)計學一個重要的隨機模擬手段,它是以區(qū)域化變量為基礎、借助變異函數(shù)來研究既具有隨機性又具有結構性的自然現(xiàn)象的手段。20世紀90年代我國在水土資源領域開始引入條件模擬進行應用研究,并受到廣泛關注[8]。在條件模擬的多種方法中,序貫高斯模擬法由于靈活、簡單、高效,成為條件模擬最常用的方法之一。

顧名思義,序貫高斯模擬就是將順序模擬和高斯模擬相結合的一種隨機模擬方法。序貫高斯模擬的適用前提是模擬變量必須符合正態(tài)分布,若不符合,則須先將其進行正態(tài)變換。該方法前人已作了大量的敘述,限于篇幅,這里不再贅述。序貫高斯模擬方法原理及算法的實現(xiàn)可參見文獻[6]。

2地表水與地下水交換量計算方法

天津市北塘水庫庫區(qū)垂向滲漏過程中,滲漏水量受下伏相對隔水層的阻擋,庫水先垂向滲漏補給庫底含水層,之后在側向水力梯度作用下以側向滲流為主。因此,在滲漏量計算過程中,可以將上述過程概化為側向滲流。

庫區(qū)滲漏量按照式(1)～(3)進行計算:

(1)

(2)

(3)

式中:Q為地下水與地表水交換量;Qi為各個分區(qū)的交換量;n為分區(qū)個數(shù);Ki為各分區(qū)水量交換層滲透系數(shù);Ai為各分區(qū)的面積;Ji為地下水與地表水之間的水力梯度,無量綱;Li為各分區(qū)的形心到對應地下水水位觀測井的距離;ΔHi為庫水位與地下水的水位差。

利用序貫高斯模擬法生成滲透系數(shù)隨機場,滲流方向的確定仍沿用文獻[7]提出的D8算法。主要計算步驟如下:

a. 將庫底地層在平面上剖分成多個大小相同的正方形網(wǎng)格,逐一計算每個網(wǎng)格中心點到周圍所有觀測井的水力梯度,取其最大值作為該點滲流方向,得到任一網(wǎng)格的Ji;

b. 將序貫高斯模擬法生成的參數(shù)模擬值分別賦值到相應的網(wǎng)格上,得到任一網(wǎng)格所代表的Ki;

c. 根據(jù)式(2)計算得到任一網(wǎng)格Ai上地表水與地下水交換量Qi;

d. 將所有網(wǎng)格各自的Qi相加即得水庫滲漏總量Q;

e. 重復步驟a～c,得到水庫滲漏量Q的多次模擬值。

3實例研究

3.1研究區(qū)概況

北塘水庫位于天津市塘沽區(qū)北塘鎮(zhèn)西北約2 km,是一座已建成的中型平原水庫。地理坐標為117°38′E～117°41′E,39°06′N～39°08′N。2003年12月北塘水庫增容擴建后,總庫容達到3 977萬m3,設計水位7.00 m,蓄水面積7.31 km2,圍堤總長約10.40 km。庫底地層巖性以黏土、粉質黏土為主,夾粉砂、粉土透鏡體,整體上含水不豐富,埋藏較淺,地下徑流微弱。庫區(qū)地下水為孔隙潛水和微承壓水,分別分布于上部陸相層和第一海相層上部的粉砂、粉土層中。庫區(qū)周邊地下水埋深0.55～2.56 m,平均1.42 m,北側和東北側地下水水位較高。北塘水庫地理位置、平面形狀及庫底地層巖芯采樣點位置見圖1。

圖1　北塘水庫地理位置、平面形狀及庫底地層巖芯采樣點位置

庫區(qū)共布置了巖芯采樣點14個,在對土樣進行室內顆粒分析試驗后,選取ALYAMANI&SEN公式最終計算得到這14個測點的滲透系數(shù)。

本文選取2005年7月16日,8月31日,9月15日地表水與地下水水位動態(tài)觀測資料進行分析計算,地表水、地下水水位動態(tài)資料見表1。

表1　地表水、地下水水位動態(tài)觀測資料　m

3.2序貫高斯模擬

滲透系數(shù)半方差函數(shù)分析及序貫高斯模擬均用GS+9.0軟件進行計算和模擬,研究區(qū)參數(shù)模擬的空間分布圖用ArcMap10.2軟件繪制,參數(shù)統(tǒng)計和P-P圖在SPSS12.0軟件中進行。

3.2.1滲透系數(shù)統(tǒng)計特征分析

為了對研究區(qū)滲透系數(shù)分布有個總體的認識,對這14個測點樣本進行了統(tǒng)計分析及非參數(shù)檢驗[9],結果表明:樣本分布類型為對數(shù)正態(tài),均值為4.56 cm/d,標準差為8.67,最大值為24.71 cm/d,最小值為8.35×10-5cm/d,變異系數(shù)為190%，屬高強度變異[10]。插值精度在很大程度上受到數(shù)據(jù)變異性的影響,在相同的插值精度條件下,數(shù)據(jù)變異性越大所需的樣本數(shù)越多。若僅對這14個測點樣本進行地質統(tǒng)計分析,則資料不夠充分,并且?guī)彀短幉荚O點數(shù)較少。為了解決此問題,可以利用這14個樣本進行克里金插值,得到如圖2所示的庫區(qū)庫底地層滲透系數(shù)等值線,然后在該圖中較均勻地取點,并結合庫底地層地質勘探資料,最終獲得20個點來插補資料系列[11]。根據(jù)國外相關研究成果[12],雖然該方法獲取的數(shù)據(jù)有一定近似性,但對數(shù)據(jù)的半定量及定性分析不會產(chǎn)生明顯影響。

圖2　庫區(qū)庫底地層滲透系數(shù)等值線(單位：cm/d)

對插補后的34個參數(shù)樣本分別繪制正態(tài)(圖3(a))和對數(shù)正態(tài)(圖3(b))P-P概率圖,相對于正態(tài)分布而言,該滲透系數(shù)系列更接近對數(shù)正態(tài)分布。將該系列取自然對數(shù)后,再對該參數(shù)系列進行K-S檢驗[13]。檢驗結果表明,在0.05顯著性水平下,lnK系列的P>0.05,接受正態(tài)分布的原假設,故lnK系列服從正態(tài)分布,滿足了序貫高斯模擬需以數(shù)據(jù)場呈正態(tài)分布的前提[14]。

圖3　K系列的正態(tài)及對數(shù)正態(tài)P-P概率圖

3.2.2半方差函數(shù)的結構分析

半方差函數(shù)是研究含水層滲透系數(shù)空間變異的一種有效方法。利用GS+9.0軟件可以很方便地獲得庫底地層滲透系數(shù)的最佳理論半方差函數(shù)(圖4)。

圖4　含水層滲透系數(shù)半方差函數(shù)

圖4中半方差函數(shù)為球狀模型,C0=2.12,C0+C=39.52,a=1 210 m,R2=0.879,RSS=29.53。其中:C0表示塊金效應;C0+C表示基臺值;C0/(C0+C)表示參數(shù)空間變異性程度,研究區(qū)滲透系數(shù)C0/(C0+C)值為0.05%,屬高度空間相關性[10]。a表示參數(shù)變程,即在該值范圍內,觀測值之間存在一定的相關性,觀測值之間的距離超過該值時則相互獨立。研究區(qū)參數(shù)變程為1 210 m,大于采樣距離。綜合以上分析,研究區(qū)滲透系數(shù)在變程1 210 m范圍內具有高度的空間相關性。

3.2.3滲透系數(shù)條件模擬

利用圖4半方差函數(shù)理論模型對庫底地層滲透系數(shù)進行了500次序貫高斯模擬。為避免模擬過程中模擬值出現(xiàn)負值[15],在GS+9.0軟件中將模擬值的上限和下限分別設為參數(shù)實測值的最大值和最小值。設定剖分網(wǎng)格大小為20 m×20 m,即每個網(wǎng)格面積為400 m2。將模擬的參數(shù)隨機場導入ArcMap中轉換成點位圖,剔除研究區(qū)外的數(shù)據(jù)。最終得到插值數(shù)據(jù)18 048個,用于滲漏量的計算。圖5列出模擬中隨機抽取的兩次條件模擬的參數(shù)空間分布圖。

圖5　滲透系數(shù)條件模擬實現(xiàn)

由圖5可知,滲透系數(shù)具有明顯的空間變異性,滲透系數(shù)較大的區(qū)域位于庫區(qū)中北部,也是水庫滲漏風險最大的區(qū)域,其他區(qū)域滲透系數(shù)較小。這是由于庫區(qū)中北部第四系全新統(tǒng)上部陸相層在水庫增容擴建期間破壞較為嚴重,下伏的粉砂和粉土層出露庫底,而其他區(qū)域地層多為黏土類沉積物,透水性差[16]。

3.3滲流方向的確定

依據(jù)文獻[7]提出的D8算法,分別計算每個網(wǎng)格中心點到周邊5口觀測井的水力梯度,取其最大值作為滲流方向,最終得到每個網(wǎng)格的Ji。圖6為不同時刻滲流方向示意圖。

圖6　不同時刻滲流方向示意圖

3.4計算滲漏量

在給每一網(wǎng)格Ji、Ki賦值的基礎上,根據(jù)式(1)～(3)分別計算3個時刻的滲漏量,這就得到每個時刻一次滲透量模擬值。按照一共模擬500次滲透系數(shù)隨機場計算,可得到每個時刻樣本容量為500的滲漏量總體。滲漏量計算結果見表2，其頻率分布見圖7。

表2　水庫滲漏量計算結果　m3/d

綜合表2及圖7分析,基于滲透系數(shù)序貫高斯模擬的滲漏量計算結果,最大值與最小值相差約1.4倍,頻率分布呈偏右型。以7月16日計算的滲漏量總體為例,最小值與最大值分別為1 389.49 m3/d、1 897.30 m3/d,且落在區(qū)間1 600～1 850 m3/d的頻率較高,其頻率高達77.6%。

圖7　2005年滲漏量頻率分布

4結論

a. 庫區(qū)滲透系數(shù)較大的區(qū)域位于中北部,也是水庫滲漏風險最大的區(qū)域,其他區(qū)域滲透系數(shù)較小。與Monte Carlo法相比,基于序貫高斯模擬的滲透系數(shù)隨機場既考慮了參數(shù)空間隨機性,又考慮了參數(shù)的結構性,其模擬結果更能反映出實際庫底地層地質條件。

b. 影響滲漏量的主要因素是滲透系數(shù)的空間變異性,同一時刻基于不同滲透系數(shù)隨機計算出的滲漏量相差約1.4倍。此外,不同時刻空間滲流方向的變化也會對滲漏量的計算產(chǎn)生影響。

c. 針對其他的水文地質參數(shù),如孔隙度、給水度、導水系數(shù)等,在對其進行正態(tài)檢驗或正態(tài)變換后,均可以通過多次序貫高斯模擬,產(chǎn)生大量不同插值結果的參數(shù)隨機場,進而進行不確定性分析。若參數(shù)樣本偏少,利用克里金插值來插補資料不失為一種有效的方法。

參考文獻：

[1] 胡立堂,王忠靜,趙建世,等.地表水和地下水相互作用及集成模型研究[J].水利學報,2007,38(1):54-59.(HU Litang,WANG Zhongjing,ZHAO Jianshi,et al.Advances in the interactions and integrated model between surface water and groundwater [J].Journal of Hydraulic Engineering,2007,38(1):54-59.(in Chinese))

[2] 胡俊鋒,王金生,滕彥國.地下水與河水相互作用的研究進展[J].水文地質工程地質,2004 (1):108-113.(HU Junfeng,WANG Jinsheng,TENG Yanguo.Study progress of interaction between stream and groundwater [J].Hydrogeology and Engineering Geology,2004 (1):108-113.(in Chinese))

[3] FREEZE R A,CHERRY J A.Groundwater [M].Upper Saddle River,NJ:Prentice-Hall,1979.

[4] 陳彥,吳吉春.含水層滲透系數(shù)空間變異性對地下水數(shù)值模擬的影響[J].水科學進展,2005,16(4):482-487.(CHEN Yan,WU Jichun.Effect of the spatial variability of hydraulic conductivity in aquifer on the numerical simulation of groundwater [J].Advances in Water Science,2005,16(4):482-487.(in Chinese))

[5] 胡伏生,陳連竹,萬力,等.含水層參數(shù)不確定性與地下水開采量可靠性評價[J].南京大學學報(自然科學版),1999,35(6):675-682.(HU Fusheng,CHEN Lianzhu,WAN Li,et al.Uncertainty of aquifer property and assessment of groundwater yield reliability [J].Journal of Nanjing University (Natural Sciences),1999,35(6):675-682.(in Chinese))

[6] 彭偉,吳劍鋒,閻婷婷.滲透系數(shù)的條件模擬對污染物運移的不確定性分析[J].南京大學學報(自然科學版),2008,44(3):280-288.(PENG Wei,WU Jianfeng,YAN Tingting.Uncertainty analysis of the contaminant transport fate using conditional simulation of hydraulic conductivity [J].Journal of Nanjing University (Natural Sciences),2008,44(3):280-288.(in Chinese))

[7] 束龍倉,魯程鵬,李偉.考慮參數(shù)不確定性的地表水與地下水交換量的計算方法[J].水文地質工程地質,2008 (5):68-71.(SHU longcang,LU Chengpeng,LI Wei.Calculation method of the exchange volume between surface water and groundwater based on uncertainty of parameters [J].Hydrogeology and Engineering Geology,2008 (5):68-71.(in Chinese))

[8] 陳亞新,徐英,史海濱,等.條件模擬原理和技術進展及在水資源系統(tǒng)中的應用[J].水科學進展,2002,13(2):246-252.(CHEN Yaxin,XU Ying,SHI Haibin,et al.Theory,development and application of the conditional simulation in water resources system [J].Advances in Water Science,2002,13(2):246-252.(in Chinese))

[9] 陳元芳.統(tǒng)計試驗方法及應用[M].哈爾濱:黑龍江人民出版社,2000:25-65.

[10] 楊奇勇,楊勁松,姚榮江,等.耕地土壤有機質空間變異性的隨機模擬[J].農(nóng)業(yè)工程學報,2010,26(12):324-329.(YANG Qiyong,YANG Jinsong,YAO Rongjiang,et al.Stochastic simulation of cultivated soil organic matter spatial variability [J].Transactions of the Chinese Society of Agricultural Engineering,2010,26(12):324-329.(in Chinese))

[11] 郭建清,母敏霞,鄭麗萍,等.淺層含水層水文地質參數(shù)的統(tǒng)計分析與空間相關性[J].勘察科學技術,2002(4):9-12.(GUO Jianqing,MU Minxia,ZHENG Liping,et al.Analysis on statistic distribution and spatial correlativity of hydrogeological parameters of shallow aquifer [J].Site Investigation Science and Technology,2002(4):9-12.(in Chinese))

[12] ROBERT J H,PETER K K.An application of geostatistical approach to the inverse problem in two-dimensional groundwater modeling [J].Water Resource Research,1984,20(7):1003-1020.

[13] 束龍倉,李偉.北塘水庫庫底地層滲透系數(shù)的隨機特性分析[J].吉林大學學報(地球科學版),2007,37(2):216-220.(SHU Longcang,LI Wei.Analysis on the stochastic characteristic of the hydrauic conductivity of bottom stratum,Beitang Reservoir,Tianjin City [J].Journal of Jilin University(Earth Science Edition),2007,37(2):216-220.(in Chinese))

[14] 李保國,胡克林,陳德立,等.農(nóng)田土壤表層飽和導水率的條件模擬[J].水利學報,2002(2):36-46.(LI Baoguo,HU Kelin,CHEN Deli,et al.Conditional simulation of soil surface saturated hydraulic conductivity at field scale [J].Journal of Hydraulic Engineering,2002(2):36-46.(in Chinese))

[15] 張澤譜,王學軍.土壤微量元素含量空間分布的條件模擬[J].土壤學報,1998,35(3):424-429.(ZHANG Zepu,WANG Xuejun.Conditional simulation for spatial distribution of trace elements in soils [J].Acta Pedologica Sinica,1998,35(3):424-429.(in Chinese))

[16] 李偉.北塘水庫對地下水的滲漏補給特性研究[D].南京:河海大學,2007.

Uncertainty analysis of reservoir leakage based on sequential Gaussian simulation of hydraulic conductivity

ZHANG Fei1,WEN Zhonghui1,DAI Fengjun2,LI Wei3

(1.CollegeofHydrologyandWaterResources,HohaiUniversity,Nanjing210098,China;2.ZhenjiangBranchofJiangsuProvinceHydrologyandWaterResourcesInvestigationBureau,Zhenjiang212003,China;3.GeologicalSurveyofJiangsuProvince,Nanjing210018,China)

Abstract：Hydraulic conductivity of the reservoir bottom stratum is one of the main factors affecting reservoir leakage. Studying the spatial variability of hydraulic conductivity can provide a scientific basis for the reliability analysis of the calculation results of reservoir leakage. Using the bottom stratum of the Beitang Reservoir as a case study, the parameter series were extended through Kriging interpolation. Then, the interpolated hydraulic conductivity was estimated by sequential Gaussian simulation which was carried out 500 times. Finally, the effect of the spatial variability of hydraulic conductivity on the uncertainty of calculated leakage was analyzed. The results show that the spatial structure of vertical hydraulic conductivity can be described by a spherical model, and it has a high spatial correlation within 1 210 m. As for the leakage calculated on July 16, 2005, the smallest value was 1 389.49 m3/d, the largest value was 1 897.30 m3/d, and the frequency of the calculated value falling into the range of 1 600 to 1 850 m3/d reached 77.6%.

Key words：hydraulic conductivity; conditional simulation; spatial variability; reservoir leakage; uncertainty

DOI：10.3880/j.issn.1004-6933.2016.03.013

基金項目:國家自然科學基金(41172203,41201029,41301017);中國水利水電科學研究院流域水循環(huán)模擬與調控國家重點實驗室開放研究基金(IWHR-SKL-201502)

作者簡介:張飛(1991—),男,碩士研究生,研究方向為地下水數(shù)值模擬及開發(fā)利用。E-mail:jeff91214@sina.com 通信作者:溫忠輝,副教授。E-mail:wenzh2812@sina.com.cn

中圖分類號：P641.2

文獻標志碼：A

文章編號：1004-6933(2016)03-0069-05

(收稿日期：2015-06-02編輯：徐娟)