許增光, 王 雪, 柴軍瑞, 覃 源
(西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710048)
目前,潔凈水資源短缺已成為一個(gè)世界性問(wèn)題,如何準(zhǔn)確掌握地下水的運(yùn)動(dòng)規(guī)律,利用地下水的天然可再生性,并通過(guò)人工調(diào)控方法對(duì)地下水資源進(jìn)行有效管理和科學(xué)利用已成為當(dāng)前水文地質(zhì)學(xué)研究的前沿及今后的發(fā)展方向[1-3]。因滲流產(chǎn)生的滲透變形、滑坡、潰壩以及地下水開(kāi)采引起地面下沉、庫(kù)區(qū)滲漏,庫(kù)岸滑坡等情形,一旦發(fā)生,往往后果嚴(yán)重,因此掌握滲流場(chǎng)分布規(guī)律尤為重要。滲流分析通常都是基于傳統(tǒng)的確定性模型,但實(shí)際滲流場(chǎng)常常受到諸多隨機(jī)性因素的影響,具有很強(qiáng)的不確定性。造成不確定性的因素主要包括:滲流參數(shù)、初始條件、邊界條件、源匯項(xiàng)等[4-5]。而滲透系數(shù)是影響滲流場(chǎng)的重要因素,因此獲取較為精確的滲透系數(shù)的分布規(guī)律至關(guān)重要。
滲透系數(shù)是表征含水層特性的一個(gè)關(guān)鍵要素。通過(guò)壓水試驗(yàn)計(jì)算滲透系數(shù)是實(shí)際工程的一貫做法,可以獲取直接的滲透系數(shù),精度高,但同時(shí)也存在成本高及操作繁瑣的問(wèn)題。加入間接測(cè)量數(shù)據(jù)電阻率來(lái)獲取滲透系數(shù),與壓水試驗(yàn)得到直接的滲透系數(shù)結(jié)合,可以在一定程度上減少打井鉆孔的數(shù)量,在操作和成本上都是有益的。特別是在大范圍區(qū)域情況下,效果會(huì)更加明顯。電阻率與滲透系數(shù)之間的關(guān)系一直在被研究,有很多經(jīng)驗(yàn)、半經(jīng)驗(yàn)公式。電阻率可由電測(cè)深、電剖面、高密度電阻率儀測(cè)量,操作簡(jiǎn)單且數(shù)據(jù)量大。文獻(xiàn)[6]利用ERT數(shù)據(jù)推求出了非均質(zhì)多孔介質(zhì)滲透系數(shù)的分布。諸多國(guó)外學(xué)者對(duì)電阻率與滲透系數(shù)之間的關(guān)系進(jìn)行了深入的研究[7-18],最后皆得到了二者之間的聯(lián)系。
本文采用兩種方法,利用電阻率測(cè)量數(shù)據(jù)求滲透系數(shù)。一種是經(jīng)驗(yàn)公式法,Archie公式利用電阻率計(jì)算孔隙率,Kozeny-Carman公式利用孔隙率計(jì)算滲透系數(shù);另一種是人工神經(jīng)網(wǎng)絡(luò)法,利用MATLAB訓(xùn)練電阻率與滲透系數(shù)的關(guān)系。然后借助Surfer軟件畫(huà)出滲透系數(shù)分布的等值線圖。最后對(duì)這兩種方法進(jìn)行對(duì)比分析。
經(jīng)驗(yàn)公式法是將Archie公式與Kozeny-Carman公式結(jié)合,充分利用了電阻率、孔隙率和滲透系數(shù)之間的關(guān)系。人工神經(jīng)網(wǎng)絡(luò)法是利用數(shù)學(xué)方法直接建立電阻率與滲透系數(shù)的數(shù)學(xué)關(guān)系式。
Archie公式:
ρ0=aρwφ-m
(1)
式中,ρ0和ρw分別是地層電阻率和地層流體電阻率(Ω·m),φ是介質(zhì)的孔隙率,a和m是與介質(zhì)有關(guān)的常數(shù)值。Archie公式只有在純砂巖地層或者泥層含量小于5%的地層中適用[19]。
Kozeny-Carman公式:
(2)
式中,k是滲透系數(shù)(m/s),δw是水的密度(1 000 kg/m3),g是重力加速度(9.81 m/s2),μ是運(yùn)動(dòng)粘滯系數(shù)(kg/m·s),d是地層巖性顆粒粒徑值(m)。
通過(guò)公式(1)求出孔隙率,代入公式(2)就可得到滲透系數(shù)。
傳輸函數(shù):
a=f(wp+b)
(3)
隱含層節(jié)點(diǎn)的輸出:
(4)
輸出節(jié)點(diǎn)的計(jì)算輸出:
(5)
誤差函數(shù):
(6)
在上述公式中,p是輸入數(shù)據(jù),w和v是權(quán)重,b是閾值,t是目標(biāo)數(shù)據(jù)。誤差函數(shù)分別對(duì)隱含層節(jié)點(diǎn)、輸出層節(jié)點(diǎn)和閾值求導(dǎo)。求導(dǎo)的目的是修正變量Δvlj和Δwji。根據(jù)方程(4),沿著負(fù)梯度方向修正w和b。訓(xùn)練直到滿足精度要求為止。
本文案例取自于文獻(xiàn)[8],研究區(qū)域位于德國(guó)波鴻魯爾大學(xué)南部的魯爾河谷,在凱姆納德湖西南部河岸,覆蓋面積大約占199 016 m2。這個(gè)區(qū)域有6口已經(jīng)存在的鉆井,6個(gè)滲透系數(shù)值作為實(shí)測(cè)值。用電測(cè)深法測(cè)量20個(gè)點(diǎn)的電阻率。已知的地質(zhì)和水文地質(zhì)資料詳見(jiàn)文獻(xiàn)[8],測(cè)量資料見(jiàn)表1[8]。
表1 測(cè)量資料Tab.1 Survey data
根據(jù)已有6口井的滲透系數(shù)反推常數(shù)a和m,最終求得a=0.8,m=1.4。結(jié)合式(1)(2)求出20個(gè)滲透系數(shù)。計(jì)算的滲透系數(shù)與電阻率的關(guān)系見(jiàn)圖1。
圖1 經(jīng)驗(yàn)公式法求解的滲透系數(shù)與電阻率的關(guān)系Fig.1 Relationship between hydraulic conductivity and resistivity calculated by empirical formula method
地層因子F=ρ/ρw。圖2是計(jì)算滲透系數(shù)與地層因子的關(guān)系。
圖2 經(jīng)驗(yàn)公式法計(jì)算的滲透系數(shù)與地層因子的關(guān)系Fig.2 Relationship between hydraulic conductivity and formation factor calculated by empirical formula method
由圖1和圖2可以看出用地層因子F表達(dá)滲透系數(shù)更為合適。經(jīng)驗(yàn)公式法滲透系數(shù)的計(jì)算值與滲透系數(shù)的觀測(cè)值之間的關(guān)系見(jiàn)圖3。因?yàn)閣6與ves1相差很大,可能是由于測(cè)量上的誤差或者是這一點(diǎn)上變異大非均質(zhì)性強(qiáng),所以在圖中未畫(huà)出點(diǎn)。
圖3 經(jīng)驗(yàn)公式法滲透系數(shù)的計(jì)算值與滲透系數(shù)的觀測(cè)值比較Fig.3 Comparison between the calculated values and the observed values of hydraulic conductivity from the empirical formula
最后滲透系數(shù)與地層因子的關(guān)系式為:
k=129.44F-2.445+1.8522
(7)
平均誤差百分比為:
(8)
式中,N是滲透系數(shù)觀測(cè)數(shù)據(jù)的個(gè)數(shù),obsi是滲透系數(shù)觀測(cè)值;compi是滲透系數(shù)計(jì)算值。
最后用6個(gè)點(diǎn)的測(cè)量數(shù)據(jù)和計(jì)算數(shù)據(jù)計(jì)算經(jīng)驗(yàn)公式法的平均誤差百分比。PE=10.82%。除去ves1的值,計(jì)算的PE=5.9%。誤差較小,證明預(yù)測(cè)的效果是較為理想的。
利用BP神經(jīng)網(wǎng)絡(luò)對(duì)電阻率進(jìn)行訓(xùn)練。公式(3)是訓(xùn)練的基礎(chǔ)公式,變量由輸入數(shù)據(jù)、權(quán)重和閾值組成。輸入變量為ρ和ρw,輸出變量是6口井觀測(cè)的滲透系數(shù)。訓(xùn)練出合適的函數(shù),然后再求20個(gè)電阻率對(duì)應(yīng)的滲透系數(shù)。
人工神經(jīng)網(wǎng)絡(luò)法滲透系數(shù)的計(jì)算值與滲透系數(shù)的觀測(cè)值之間的關(guān)系見(jiàn)圖4。兩種方法得到的20個(gè)滲透系數(shù)數(shù)據(jù)見(jiàn)表2。
圖4 人工神經(jīng)網(wǎng)絡(luò)法滲透系數(shù)的計(jì)算值與滲透系數(shù)的觀測(cè)值比較Fig.4 Comparison between the calculated values and the observed values of hydraulic conductivity from artificial neural network
點(diǎn)數(shù)經(jīng)驗(yàn)公式法k/(10-2m/s)人工神經(jīng)網(wǎng)絡(luò)法k/(10-2m/s)1292283220522932352544334268536525762052247196215825727194553551024326911199220122262301331824214259238152132441621018317226276182341841920322020438242
按公式(8)計(jì)算人工神經(jīng)網(wǎng)絡(luò)法的平均誤差百分比,PE=9.23%。除去ves1的值,計(jì)算的PE=3.2%。這個(gè)誤差百分比也比較小,預(yù)測(cè)效果比較理想。
利用地質(zhì)統(tǒng)計(jì)法中的克里格插值分別對(duì)這兩種方法計(jì)算出的滲透系數(shù)分布。Surfer軟件有很強(qiáng)大的插值功能,并且可以繪制等值線圖。圖5和圖6分別是經(jīng)驗(yàn)公式法和人工神經(jīng)網(wǎng)絡(luò)得到的滲透系數(shù)分布圖。
圖5 經(jīng)驗(yàn)公式法得到的滲透系數(shù)分布Fig.5 Hydraulic conductivity distribution obtained by empirical formula method
圖6 神經(jīng)網(wǎng)絡(luò)法得到的滲透系數(shù)分布Fig.6 Hydraulic conductivity distribution obtained by artificial neural network method
圖5與圖6分布圖,證明兩種預(yù)測(cè)方法是可行的,基于實(shí)測(cè)的電阻率數(shù)據(jù)可以得到滲透系數(shù)的分布。不僅操作簡(jiǎn)單還可以獲取大量的數(shù)據(jù),效果很理想。從中也可以看出經(jīng)驗(yàn)公式法得到的滲透系數(shù)值比人工神經(jīng)網(wǎng)絡(luò)法大,而人工神經(jīng)網(wǎng)絡(luò)法分布較為均勻。在本工程案例上,人工神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的效果最佳。計(jì)算值與觀測(cè)值相差較小,其平均誤差小于經(jīng)驗(yàn)公式法的平均誤差。
通過(guò)對(duì)兩種預(yù)測(cè)方法的對(duì)比分析,經(jīng)驗(yàn)公式法是依據(jù)電阻率、孔隙率和滲透系數(shù)之間的關(guān)系而定,預(yù)測(cè)值比較大,這與a和m的確定和地質(zhì)條件都有關(guān)系。神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)是依賴目標(biāo)數(shù)據(jù),由于觀測(cè)數(shù)據(jù)分布均勻(1.83~2.8 cm/s)而且數(shù)量較少(6個(gè)),這導(dǎo)致神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的整體效果也比較均勻。但是總體的趨勢(shì)是一致的。從平均誤差來(lái)看,本案例預(yù)測(cè)效果最佳的是人工神經(jīng)網(wǎng)絡(luò)法。
本文主要是采取兩種方法利用電阻率數(shù)據(jù)得到滲透系數(shù)數(shù)據(jù)。從而減少壓水試驗(yàn)帶來(lái)的高成本和復(fù)雜操作問(wèn)題??梢越?jīng)濟(jì)高效地獲取滲透系數(shù)的分布。
經(jīng)驗(yàn)公式法是結(jié)合Archie公式和Kozeny-Carman公式,結(jié)合前人的經(jīng)驗(yàn),理論方法簡(jiǎn)單。但是在a和m的確定上存在著誤差,a和m的值影響著最后的計(jì)算結(jié)果。Archie公式不適合在黏土中使用,所以經(jīng)驗(yàn)公式法是不適合黏土地層地質(zhì)情況。同時(shí)也不適合太復(fù)雜地層,影響滲透系數(shù)的因素不僅僅是孔隙率,影響的因素越多,預(yù)測(cè)精度就會(huì)越低。
人工神經(jīng)網(wǎng)絡(luò)法是把工程問(wèn)題數(shù)學(xué)化,采用人工神經(jīng)網(wǎng)絡(luò),利用計(jì)算機(jī)龐大的計(jì)算能力訓(xùn)練出電阻率與滲透系數(shù)之間的關(guān)系,但人工神經(jīng)網(wǎng)絡(luò)法不適合在沒(méi)有已知的滲透系數(shù)情況下使用。
兩種方法得到的ves1的滲透系數(shù)與測(cè)量值相差都很大,可能這一點(diǎn)上滲透系數(shù)的觀測(cè)值存在著誤差或者是此處的非均質(zhì)變異性較其他地方強(qiáng)。根據(jù)平均誤差百分比的計(jì)算,兩種方法的平均誤差都較小。根據(jù)最后的插值等值線圖,預(yù)測(cè)的結(jié)果很理想。
在這個(gè)工程案例下兩種預(yù)測(cè)方法都是可行的,而且預(yù)測(cè)的效果都是很理想的。雖然兩者之間有一些偏差,可能是由于研究區(qū)域的面積大,已知數(shù)據(jù)太少造成的,再加上測(cè)量出現(xiàn)一些誤差等情況都會(huì)導(dǎo)致結(jié)果的不準(zhǔn)確。
參考文獻(xiàn):
[1] JORDI B A. Importance of stream infiltration data for modelling surface water-groundwater interactions [J]. Journal of Hydrology, 2015, 528: 683-693.
[2] 張博, 李國(guó)秀, 程品, 等. 基于隨機(jī)理論的地下水環(huán)境風(fēng)險(xiǎn)評(píng)價(jià) [J]. 水科學(xué)進(jìn)展, 2016, 27(1): 100-106.
ZHANG Bo, LI Guoxiu, CHENG Pin, et al. Groundwater environment risk assessment based on stochastic theory [J]. Advances in Water Science, 2016, 27(1): 100-106.
[3] 王林, 徐青. 基于蒙特卡羅隨機(jī)有限元法的三維隨機(jī)滲流場(chǎng)研究 [J]. 巖土力學(xué), 2014, 35(1): 287-292.
WANG Lin, XU Qing. Analysis of three dimensional random seepage field based on Monte Carlo stochastic finite element method [J]. Rock and Soil Mechanics, 2014, 35(1): 287-292.
[4] 覃榮高, 曹廣祝, 仵彥卿. 非均質(zhì)含水層中滲流與溶質(zhì)運(yùn)移研究進(jìn)展 [J]. 地球科學(xué)進(jìn)展, 2014, 29(1): 30-41.
QIN Ronggao, CAO Guangzhu, WU Yanqing. Review of the study of groundwater flow and solute transport in heterogeneous aquifer [J]. Advances in Earth Science, 2014, 29(1): 30-41.
[5] 梁婕, 曾光明, 郭生練, 等. 滲透系數(shù)的非均質(zhì)性對(duì)地下水溶質(zhì)運(yùn)移的影響 [J]. 水利學(xué)報(bào), 2008, 39(8): 900-906.
LIANG Jie, ZENG Guangming, GUO Shenglian, et al. Effect of hydraulic conductivity heterogeneity on solute transport in groundwater [J]. Journal of Hydraulic Engineering, 2008, 39(8): 900-906.
[6] 喻永祥, 吳吉春. 利用ERT數(shù)據(jù)推求非均質(zhì)多孔介質(zhì)滲透系數(shù)初探 [J]. 水文地質(zhì)工程地質(zhì), 2006, 33(2): 41-44.
YU Yongxiang, WU Jichun. Application of electrical resistivity tomography data to estimate hydraulic conductivity in porous medium [J]. Hydrogeology and Engineering Geology, 2006, 33(2): 41-44.
[7] SUN K, GOLTZ M N. Direct estimation of hydraulic parameters relating to steady state groundwater flow [J]. Environmental Modelling &Software, 2016, 86: 50-55.
[8] NIWAS S, CELIK M. Equation estimation of porosity and hydraulic conductivity of Ruhrtal aquifer in Germany using near surface geophysics [J]. Journal of Applied Geophysics, 2012, 84(9): 77-85.
[9] OLATUNJI S, MUSA A. Estimation of aquifer hydraulic characteristics from surface geoelectrical methods: case study of the Rima Basin, north western Nigeria [J]. Arabian Journal for Science and Engineering, 2014, 39(7): 5475-5487.
[10] OKIONGBO K, MEBINE P. Estimation of aquifer hydraulic parameters from geoelectrical method-a case study of Yenagoa and environs, Southern Nigeria [J]. Arabian Journal of Geosciences, 2015, 8(8): 6085-6093.
[11] MONDAL N C, DEVI A B, RAJ P A, et al. Estimation of aquifer parameters from surfacial resistivity measurement in granitic area in Tamil Nadu [J]. Current Science, 2016, 111(3): 524-534.
[12] SINGH S, SINGH V S. Estimation of hydraulic characteristics from electrical resistivity data in coastal aquifers of southern India [J]. Journal Geological Society of India, 2016, 88(1): 77-86.
[13] TIZRO A T, VOUDOURIS K, BASAMI Y. Estimation of porosity and specific yield by application of geoelectrical method-A case study in western Iran [J]. Journal of Hydrology, 2012, 454-455(4): 160 -172.
[14] KAZAKIS N, VARGEMEZIS G, VOUDOURIS K. Estimation of hydraulic parameters in a complex porous aquifer system using geoelectrical methods [J]. Science of the Total Environment,2016,550:742 -750.
[15] PERDOMO S, AINCHIL J, KRUS E. Hydraulic parameters estimation from well logging resistivity and geoelectrical measurements [J]. Journal of Applied Geophysics, 2014, 105(2): 50-58.
[16] BOUDERBALA A, REMINI B, HAMOUDI S. Geoelectrical investigation of saline water intrusion into freshwater aquifers: A case study of Nador coastal aquifer, Tipaza, Algeria [J]. Geofísica Internacional, 2016, 55(4): 239-253.
[17] TIZRO A T, VOUDOURIS K, SALEHZADE M, et al. Hydrogeological framework and estimation of aquifer hydraulic parameters using geoelectrical data: a case study from West Iran [J]. Hydrogeology Journal, 2010, 18(4): 917-929.
[18] ASFAHANI J. Hydraulic parameters estimation by using an approach based on vertical electrical soundings (VES) in the semi-arid Khanasser valley region, Syria [J]. Journal of African Earth Sciences, 2016, 117: 196-206.
[19] 楊飛. 地球物理資料綜合解釋 [M]. 北京: 石油工業(yè)出版社, 1981: 158-159.