基于地統(tǒng)計(jì)學(xué)的伊犁-鞏乃斯河谷滲透系數(shù)空間變異性研究

宋 浩 ， 張 琛， 劉 明, 曾 磊， 李 瑛， 馬稚桐

(1. 長安大學(xué) 旱區(qū)地下水與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710061； 2.長安大學(xué) 環(huán)境科學(xué)與工程學(xué)院， 陜西 西安 710061； 3.中國地質(zhì)調(diào)查局 西安地質(zhì)調(diào)查中心， 陜西 西安 710054)

1 研究背景

水文地質(zhì)參數(shù)的空間變異性是影響飽和及非飽和帶地下水水流和溶質(zhì)運(yùn)移不確定性的主要因素[1]，是開展區(qū)域數(shù)值模擬、地下水資源形成演化以及水資源評(píng)價(jià)的基礎(chǔ)，而這種空間變異性大多都具有很強(qiáng)的結(jié)構(gòu)性[2]。20世紀(jì)70年代以來，自DAGAN率先提出了隨機(jī)理論研究方法，使得空間變異性得以在數(shù)學(xué)模型的基礎(chǔ)上進(jìn)行解析，隨后，空間變異性便成了研究領(lǐng)域的熱點(diǎn)問題。國內(nèi)外眾多學(xué)者[3-16]注重于研究河床滲透系數(shù)、表土飽和滲透系數(shù)以及農(nóng)業(yè)小流域的滲透系數(shù)，如Fu Tonggang等[7]和Botros等[8]運(yùn)用實(shí)驗(yàn)變差函數(shù)擬合了喀斯特地貌下土壤飽和滲透系數(shù)的最優(yōu)模型來表征該地區(qū)的空間變異性，并發(fā)現(xiàn)最優(yōu)模型均是高斯模型。黃冠華等[11]通過測量試驗(yàn)地土壤相關(guān)的水文地質(zhì)參數(shù)，發(fā)現(xiàn)淺層土壤的滲透系數(shù)大致呈對(duì)數(shù)正態(tài)分布。而施小清等[12]和商瑋麟等[13]通過研究發(fā)現(xiàn)，承壓含水層的滲透系數(shù)樣本數(shù)據(jù)大多服從正態(tài)分布或者對(duì)數(shù)正態(tài)分布。

已有的研究中鮮有使用正態(tài)性更好的Box-Cox變換，承壓含水層的空間變異性研究也比較少，并且傳統(tǒng)的水文地質(zhì)參數(shù)分區(qū)大多是以地形地貌作為參考依據(jù)。本文旨在對(duì)伊犁-鞏乃斯河谷基于鉆孔的滲透系數(shù)進(jìn)行相關(guān)分析，探討伊犁河谷潛水含水層滲透系數(shù)的空間分布規(guī)律，尋求一種利用有限水文地質(zhì)資料獲取區(qū)域滲透系數(shù)的方法，并為野外水文地質(zhì)工作提供一個(gè)可供參考的滲透系數(shù)分區(qū)。

2 研究區(qū)概況

伊犁河谷區(qū)隸屬于伊犁哈薩克自治州。地理坐標(biāo)范圍為80°9′42″～91°01′45″E，40°14′16″～49°10′45″N，總面積為19 900 km2。本文主要以伊犁-鞏乃斯河谷作為研究區(qū)，研究區(qū)總面積約為10 000 km2，見圖1。

研究區(qū)西端總寬度約為90 km，向東延伸，寬度逐漸減小。地層結(jié)構(gòu)呈現(xiàn)單層、多層交替變化，巖層顆粒也粗細(xì)交替，并在最東端形成斷陷谷地。在這些地層中，賦存有大量的淺層潛水以及深層的承壓水，通過大范圍的水文地質(zhì)調(diào)查，山前沖洪積扇區(qū)地下水位大于100m，而河谷區(qū)兩側(cè)地下水位只有不足15 m。伊犁河谷區(qū)的地下水補(bǔ)給來源主要為北部以及南部山區(qū)的冰雪消融，整個(gè)谷地地下水總流向與伊犁河一致，為自東向西徑流。

研究區(qū)的第四層地質(zhì)大體可以描述為下更新統(tǒng)西域礫巖主要分布于南部庫克蘇河源頭，中更新統(tǒng)冰水沉積主要分布在鞏乃斯谷地兩側(cè)，上更新統(tǒng)的風(fēng)積層主要分布于霍城以北，沖洪積層主要分布于烏孫山北麓以及特克斯以東，全新統(tǒng)風(fēng)積層分布于伊寧市以西的伊犁河兩側(cè)。

圖1 研究區(qū)位置與抽水試驗(yàn)點(diǎn)分布圖

3 數(shù)據(jù)來源及方法

通過收集整理伊犁河谷地區(qū)的抽水試驗(yàn)資料，共獲取52個(gè)潛水含水層和31個(gè)承壓含水層的滲透系數(shù)數(shù)據(jù)。抽水試驗(yàn)點(diǎn)分布圖見圖1。本文用到的主要理論及方法包括Box-Cox變化、變差函數(shù)理論等。

3.1 Box-Cox變換

本文中由于潛水含水層滲透系數(shù)既不滿足正態(tài)分布，也不滿足對(duì)數(shù)正態(tài)分布，所以特引入Box-Cox變換從而提高了樣本的正態(tài)性。 變換[17]就是對(duì)反應(yīng)變量y進(jìn)行變換,變換的公式為：

(1)

式中：λ為待定參數(shù),λ>0;y為反應(yīng)變量。

由公式(1)可以看出對(duì)于不同的λ會(huì)產(chǎn)生不一樣的變換，很明顯，λ=0就是對(duì)數(shù)變換，可見Box-Cox變換就是一種比對(duì)數(shù)變換更高級(jí)的變換。λ一般采用最大似然估計(jì)求得。

3.2 變差函數(shù)理論

變差函數(shù)是一個(gè)分析區(qū)域變量隨機(jī)性與結(jié)構(gòu)性特別有效的地質(zhì)統(tǒng)計(jì)學(xué)方法。該方法假設(shè)區(qū)域變化量Z(x)是二階平穩(wěn)[18]的，則其變差函數(shù)定義為：

(2)

假設(shè)Z(xi)是平穩(wěn)區(qū)域化變量Z(x)在n個(gè)點(diǎn)上的觀測值，數(shù)據(jù)對(duì){Z(xi),Z(xi+h)}為在某一方向上相隔|h|的點(diǎn)對(duì)(xi,xi+h)上的觀測值，則定義Z(x)的實(shí)驗(yàn)變差函數(shù)為：

(3)

式中:γ*(h)為h點(diǎn)的實(shí)驗(yàn)變異函數(shù);N(h)為數(shù)據(jù)對(duì){Z(xi),Z(xi+h}的對(duì)數(shù)。

3.3 變差函數(shù)的理論模型

變差函數(shù)的理論模型[19]主要分為有基臺(tái)值和無基臺(tái)值模型兩大類，本文主要研究的是有基臺(tái)值的情況，所以就只羅列有基臺(tái)值的三大模型，見表1。

表1中，各公式中a，C，C0分別為各自變差函數(shù)模型中的數(shù)值。

4 結(jié)果分析與討論

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

對(duì)獲取的樣本數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，結(jié)果分別見表2以及圖2、3。

從表2可以看出，該地區(qū)的潛水含水層滲透系數(shù)分布極度不均，滲透系數(shù)值分布在0.46～40.68 m/d之間，進(jìn)一步表明研究區(qū)的潛水含水介質(zhì)變化較大，既有類似砂礫等粗顆粒介質(zhì)，也有類似于黏土等截水能力較強(qiáng)的細(xì)顆粒介質(zhì)。承壓含水層的滲透系數(shù)統(tǒng)計(jì)分析結(jié)果表明，滲透系數(shù)分布在0.54～20.40 m/d之間，相比于潛水含水層，這個(gè)值較穩(wěn)定，主要原因是承壓含水介質(zhì)較潛水含水層來說，在整個(gè)地區(qū)的更替變化較小，地層結(jié)構(gòu)更加穩(wěn)定。從計(jì)算結(jié)果來看，潛水和承壓水的原始K值以及變換后的K值都在0.1～1.0范圍，但是很明顯的區(qū)別是進(jìn)行Box-Cox變換后，變異系數(shù)達(dá)到了最小，可見Box-Cox變換可以減小樣本的離散度。

由圖2、3可以很直接地看出來，首先是不論潛水還是承壓水，原始K值直方圖擬合的鐘型曲線效果都不好，正態(tài)性比較差，主要還是由于K的樣本離散度太大。此外，由Q-Q圖(描述正態(tài)分布的散點(diǎn)圖)也可以很明顯的看出來，K值的正態(tài)概率分布穩(wěn)定性很差，體現(xiàn)在樣本點(diǎn)在斜率為標(biāo)準(zhǔn)差、截距為均值的直線上比較分散。而對(duì)K值進(jìn)行對(duì)數(shù)轉(zhuǎn)換和Box-Cox轉(zhuǎn)換后，很明顯地可以看出來鐘型曲線擬合效果有所提升，尤其是Box-Cox變換，達(dá)到最優(yōu)效果，這正體現(xiàn)出了Box-Cox變換對(duì)于樣本正態(tài)性提升的優(yōu)越性。由公式(1)來看，對(duì)數(shù)轉(zhuǎn)換是Box-Cox變換的一種特殊形式，對(duì)于離散度較高的樣本數(shù)據(jù)，Box-Cox變換則能顯示出它的優(yōu)勢，這點(diǎn)從圖2可以直接觀察出來。但是在類似于承壓含水層滲透系數(shù)穩(wěn)定性較好情況下，對(duì)數(shù)變換能達(dá)到最優(yōu)轉(zhuǎn)換的效果。從表2可以看出，對(duì)數(shù)變換和Box-Cox變換的結(jié)果基本上是一致的，由圖3也可以看出，兩者的正態(tài)性相似。

運(yùn)用K-S檢驗(yàn)和S-W檢驗(yàn)對(duì)K、lnK、以及進(jìn)行Box-Cox變換后的K值進(jìn)行正態(tài)性檢驗(yàn)，結(jié)果見表3。

由表3可知，潛水含水層滲透系數(shù)只有進(jìn)行Box-Cox變換后，樣本顯著性Sig值才大于0.05，說明潛水含水層滲透系數(shù)只有進(jìn)行Box-Cox變換后，才會(huì)在95%的置信水平下服從于正態(tài)分布；如前面所言，由于承壓含水層滲透系數(shù)離散度低，對(duì)數(shù)變換就可以服從正態(tài)分布。因此，對(duì)于潛水含水層滲透系數(shù)，分析過程中需要對(duì)原有數(shù)據(jù)進(jìn)行Box-Cox變換，而承壓含水層滲透系數(shù)，分析過程中只需要對(duì)其做對(duì)數(shù)變換即可。

4.2 實(shí)驗(yàn)變差函數(shù)的計(jì)算與分析

前面的研究發(fā)現(xiàn)伊犁-鞏乃斯河谷的潛水含水層滲透系數(shù)服從于Box-Cox變換正態(tài)分布，承壓含水層滲透系數(shù)服從對(duì)數(shù)正態(tài)分布。運(yùn)用公式(3)計(jì)算實(shí)驗(yàn)變差函數(shù)[20]，以Box-Cox變換后的潛水含水層K值和對(duì)數(shù)變換后的承壓含水層K值作為區(qū)域變化量，在進(jìn)行基本的統(tǒng)計(jì)之后，發(fā)現(xiàn)樣本值較為平均，符合平穩(wěn)假設(shè)的條件。運(yùn)用GS+軟件進(jìn)行擬合，潛水含水層樣本數(shù)據(jù)選用滯后距為143 412 m，間距為13 200 m，承壓含水層樣本數(shù)據(jù)選用滯后距為123 987 m，間距為10 300 m，方向容限都為22.5°的情況下，計(jì)算該地區(qū)的潛水和承壓含水層樣本值的實(shí)驗(yàn)變差函數(shù)，得到變差函數(shù)特征參數(shù)的計(jì)算結(jié)果與理論模型擬合結(jié)果，見表4。

表2 滲透系數(shù)K統(tǒng)計(jì)分析結(jié)果

圖2 潛水含水層滲透系數(shù)直方圖與Q-Q圖

圖3 承壓含水層滲透系數(shù)直方圖與Q-Q圖

從表4可以看出來，潛水含水層和承壓含水層滲透系數(shù)的擬合結(jié)果均表明球狀、指數(shù)和高斯模型的擬合優(yōu)度相差不大，但是三者的變程相差比較大，很明顯，高斯模型的變程a都是最小的，并且C0/(C0+C)都是最高的，所以選擇高斯模型作為擬合結(jié)果的最優(yōu)模型。依據(jù)Iqbal等[21]所闡述的觀點(diǎn)，如果C0/(C0+C)<25%，說明系統(tǒng)的空間相關(guān)性很強(qiáng)烈。對(duì)比計(jì)算結(jié)果，可以看出，潛水含水層和承壓含水層滲透系數(shù)的空間相關(guān)性都很強(qiáng)烈，依據(jù)相關(guān)研究[22]可知，這些變異主要由結(jié)構(gòu)性因素所引起，也就是受地層結(jié)構(gòu)的影響較大。圖4給出了實(shí)驗(yàn)變差函數(shù)及理論模型的擬合圖。

表3 滲透系數(shù)正態(tài)性檢驗(yàn)結(jié)果

表4 潛水與承壓水K值變差函數(shù)特征參數(shù)的計(jì)算與擬合

圖4 實(shí)驗(yàn)變差函數(shù)及理論模型擬合圖

結(jié)合表4與圖4來看，潛水含水層的高斯模型擬合優(yōu)度達(dá)到了0.378，C0/(C0+C)達(dá)到了18.85%，說明潛水含水層K值的Box-Cox變換空間相關(guān)性較高；變程a=10 046 m，所以當(dāng)兩個(gè)樣本點(diǎn)的距離小于10 046 m時(shí)是空間相關(guān)的。承壓含水層的高斯模型擬合優(yōu)度達(dá)0.441，C0/(C0+C)達(dá)到了9.20%，說明承壓含水層K值的對(duì)數(shù)變換空間相關(guān)性較好；變程a=9 180 m，所以說當(dāng)兩個(gè)樣本點(diǎn)的距離小于9 180 m時(shí)是空間相關(guān)的。

4.3 滲透系數(shù)參數(shù)分區(qū)

4.3.1 潛水含水層滲透系數(shù)參數(shù)分區(qū) 利用克里金插值法可以得到該地區(qū)潛水含水層K值的等值線圖，插值結(jié)果見圖5。

圖5中由于潛水含水層只有Box-Cox值服從正態(tài)分布，所以潛水含水層插值數(shù)據(jù)采用的是滲透系數(shù)的Box-Cox變換值，相應(yīng)地，承壓含水層采用的是K值的對(duì)數(shù)變換值。

由圖5的插值結(jié)果可以看出來，研究區(qū)潛水含水層的滲透系數(shù)呈現(xiàn)出霍城縣以西的滲透系數(shù)較高，而在伊寧縣附近滲透系數(shù)突然變得很小，產(chǎn)生了整個(gè)研究區(qū)的最小的滲透系數(shù)區(qū)。察縣和喀什河與伊犁河交界點(diǎn)偏南部的滲透系數(shù)又很高，而在整個(gè)鞏乃斯河附近的滲透系數(shù)變化較穩(wěn)定，滲透系數(shù)呈現(xiàn)出西部偏大而東部偏小的特點(diǎn)。結(jié)合研究區(qū)的地貌類型見圖6。

由地貌類型圖6可以看出，該地區(qū)的地層巖性從西向東主要是呈現(xiàn)出顆粒越來越細(xì)的特點(diǎn)。有研究表明，大尺度的滲透系數(shù)的空間變異性都是以地貌類型條件作為控制性因素，正是由于地貌類型的千差萬別，導(dǎo)致河流的水動(dòng)力條件發(fā)生了改變，從而使得沉積物巖性在空間上分布不均，直接導(dǎo)致滲透系數(shù)的大小也分布不均。研究區(qū)的礫質(zhì)細(xì)土平原區(qū)占地面積較大，也就是說，較粗顆粒的地貌類型占研究區(qū)比例較大，正因?yàn)槿绱?，研究區(qū)總體應(yīng)該呈現(xiàn)出西部地區(qū)的滲透系數(shù)大于東部地區(qū)的滲透系數(shù)的特點(diǎn)，這與克里金插值結(jié)果一致。圖5中滲透系數(shù)最小的區(qū)域?yàn)閳D6中的山前黃土丘陵區(qū)，圖5中的3個(gè)滲透系數(shù)最大值，全都在礫質(zhì)平原區(qū)。東部地區(qū)由于主要地貌類型為山前黃土丘陵以及含礫細(xì)土平原區(qū)，所以滲透系數(shù)相比于西部地區(qū)普遍較小。結(jié)合圖5、6，對(duì)研究區(qū)的潛水含水層的滲透系數(shù)進(jìn)行參數(shù)的分區(qū)，分區(qū)結(jié)果見圖7。

圖7的分區(qū)主要是依據(jù)圖5所示的滲透系數(shù)的插值結(jié)果以及研究區(qū)的地形地貌和水文地質(zhì)條件，R1～R8表示潛水含水層滲透系數(shù)各個(gè)不同的分區(qū)，各分區(qū)的滲透系數(shù)值是按照插值結(jié)果的平均值進(jìn)行確定的，結(jié)果見表5，由于研究區(qū)的樣本數(shù)據(jù)比較少，只能分出一個(gè)大概的框架。如果要進(jìn)行比較細(xì)致的水文地質(zhì)工作，需要收集更詳細(xì)的資料來對(duì)研究區(qū)進(jìn)行更加細(xì)致地劃分。

表5 潛水含水層滲透系數(shù)K分區(qū)結(jié)果 m/d

4.3.2 承壓含水層滲透系數(shù)參數(shù)分區(qū) 利用克里金插值法可以得到該地區(qū)承壓含水層K值的等值線圖，見圖8。從圖8可以看出，整個(gè)區(qū)域滲透系數(shù)最大值點(diǎn)分布在察縣東部以及新源縣北部偏西處，最小值區(qū)域分布在霍城縣西部。承壓含水層的滲透系數(shù)主要受地區(qū)深部地層結(jié)構(gòu)以及沉積物類型的影響較大。經(jīng)過調(diào)查得知，研究區(qū)內(nèi)上游主要為卵礫石層，山前洪積扇中部以下和各大河流中下游逐漸變?yōu)樯暗[石、黏土等互層結(jié)構(gòu)。由此可見，該地區(qū)的深層沉積物主要表現(xiàn)為從東至西顆粒逐漸變細(xì)，且沉積物厚度越來越大，從而導(dǎo)致了該地區(qū)的含水介質(zhì)滲透系數(shù)呈現(xiàn)出東部大于西部的特點(diǎn)。對(duì)研究區(qū)的承壓含水層滲透系數(shù)進(jìn)行參數(shù)分區(qū)，分區(qū)結(jié)果見圖9。

圖5 潛水含水層滲透系數(shù)插值結(jié)果圖 圖6 研究區(qū)主要地貌類型圖

圖7 潛水含水層滲透系數(shù)分區(qū)圖 圖8 承壓水含水層滲透系數(shù)插值結(jié)果圖

由于承壓含水層的沉積物類型沒有可供參考的資料，所以承壓含水層的分區(qū)主要是依據(jù)插值結(jié)果圖以及滲透系數(shù)的大小進(jìn)行分區(qū)，R1-R6表示承壓含水層滲透系數(shù)不同的分區(qū)，分區(qū)結(jié)果見表6。

圖9 承壓含水層滲透系數(shù)分區(qū)圖

表6 承壓含水層滲透系數(shù)K分區(qū)結(jié)果m/d

5 結(jié) 論

采用地質(zhì)統(tǒng)計(jì)學(xué)對(duì)伊犁河谷的潛水和承壓水含水層滲透系數(shù)進(jìn)行了正態(tài)性檢驗(yàn)，用實(shí)驗(yàn)變差函數(shù)理論進(jìn)行了模型擬合，結(jié)合區(qū)域地貌地形以及克里金插值結(jié)果對(duì)該地區(qū)的滲透系數(shù)進(jìn)行了分區(qū)，取得了如下認(rèn)識(shí)：

(1)伊犁-鞏乃斯河谷地區(qū)滲透系數(shù)樣本數(shù)據(jù)經(jīng)Box-Cox變換后符合正態(tài)性分布；而承壓含水層滲透系數(shù)樣本經(jīng)對(duì)數(shù)變換即可滿足正態(tài)性分布。

(2)伊犁-鞏乃斯河谷地區(qū)的滲透系數(shù)顯示出了強(qiáng)烈的空間相關(guān)性，并且該地區(qū)滲透系數(shù)的空間變異性主要是由巖層的結(jié)構(gòu)本身引起的。

(3)基于克里金插值方法，結(jié)合研究區(qū)地形地貌特征，開展了研究區(qū)滲透系數(shù)分區(qū)工作，分區(qū)結(jié)果表明：研究區(qū)總體應(yīng)該呈現(xiàn)出西部地區(qū)的滲透系數(shù)大于東部地區(qū)的滲透系數(shù)的特點(diǎn)，與克里金插值結(jié)果一致。分區(qū)結(jié)果顯示該地區(qū)總體滲透系數(shù)較大，最大值在R3區(qū)，達(dá)到了31 m/d；最小值在R5區(qū)，僅為0.7 m/d。

(4)根據(jù)調(diào)查結(jié)果，承壓含水層的地層巖性基本表現(xiàn)為從東至西逐漸變細(xì)的特點(diǎn)，相應(yīng)地其滲透系數(shù)也在逐漸變小，這與插值結(jié)果也比較符合。分區(qū)結(jié)果顯示承壓含水層滲透系數(shù)總體偏小，最大值在R6區(qū)，可以達(dá)到13.5 m/d；最小值在R1區(qū)，僅為1.1 m/d。