基于地統計學的伊犁-鞏乃斯河谷滲透系數空間變異性研究

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

(1. 長安大學 旱區地下水與生態效應教育部重點實驗室， 陜西 西安 710061； 2.長安大學 環境科學與工程學院， 陜西 西安 710061； 3.中國地質調查局 西安地質調查中心， 陜西 西安 710054)

1 研究背景

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

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

2 研究區概況

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

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

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

圖1 研究區位置與抽水試驗點分布圖

3 數據來源及方法

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

3.1 Box-Cox變換

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

(1)

式中：λ為待定參數,λ>0;y為反應變量。

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

3.2 變差函數理論

變差函數是一個分析區域變量隨機性與結構性特別有效的地質統計學方法。該方法假設區域變化量Z(x)是二階平穩[18]的，則其變差函數定義為：

(2)

假設Z(xi)是平穩區域化變量Z(x)在n個點上的觀測值，數據對{Z(xi),Z(xi+h)}為在某一方向上相隔|h|的點對(xi,xi+h)上的觀測值，則定義Z(x)的實驗變差函數為：

(3)

式中:γ*(h)為h點的實驗變異函數;N(h)為數據對{Z(xi),Z(xi+h}的對數。

3.3 變差函數的理論模型

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

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

4 結果分析與討論

4.1 滲透系數特征統計分析

對獲取的樣本數據進行統計分析，結果分別見表2以及圖2、3。

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

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

運用K-S檢驗和S-W檢驗對K、lnK、以及進行Box-Cox變換后的K值進行正態性檢驗，結果見表3。

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

4.2 實驗變差函數的計算與分析

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

表2 滲透系數K統計分析結果

圖2 潛水含水層滲透系數直方圖與Q-Q圖

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

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

表3 滲透系數正態性檢驗結果

表4 潛水與承壓水K值變差函數特征參數的計算與擬合

圖4 實驗變差函數及理論模型擬合圖

結合表4與圖4來看，潛水含水層的高斯模型擬合優度達到了0.378，C0/(C0+C)達到了18.85%，說明潛水含水層K值的Box-Cox變換空間相關性較高；變程a=10 046 m，所以當兩個樣本點的距離小于10 046 m時是空間相關的。承壓含水層的高斯模型擬合優度達0.441，C0/(C0+C)達到了9.20%，說明承壓含水層K值的對數變換空間相關性較好；變程a=9 180 m，所以說當兩個樣本點的距離小于9 180 m時是空間相關的。

4.3 滲透系數參數分區

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

圖5中由于潛水含水層只有Box-Cox值服從正態分布，所以潛水含水層插值數據采用的是滲透系數的Box-Cox變換值，相應地，承壓含水層采用的是K值的對數變換值。

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

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

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

表5 潛水含水層滲透系數K分區結果 m/d

4.3.2 承壓含水層滲透系數參數分區 利用克里金插值法可以得到該地區承壓含水層K值的等值線圖，見圖8。從圖8可以看出，整個區域滲透系數最大值點分布在察縣東部以及新源縣北部偏西處，最小值區域分布在霍城縣西部。承壓含水層的滲透系數主要受地區深部地層結構以及沉積物類型的影響較大。經過調查得知，研究區內上游主要為卵礫石層，山前洪積扇中部以下和各大河流中下游逐漸變為砂礫石、黏土等互層結構。由此可見，該地區的深層沉積物主要表現為從東至西顆粒逐漸變細，且沉積物厚度越來越大，從而導致了該地區的含水介質滲透系數呈現出東部大于西部的特點。對研究區的承壓含水層滲透系數進行參數分區，分區結果見圖9。

圖5 潛水含水層滲透系數插值結果圖 圖6 研究區主要地貌類型圖

圖7 潛水含水層滲透系數分區圖 圖8 承壓水含水層滲透系數插值結果圖

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

圖9 承壓含水層滲透系數分區圖

表6 承壓含水層滲透系數K分區結果m/d

5 結 論

采用地質統計學對伊犁河谷的潛水和承壓水含水層滲透系數進行了正態性檢驗，用實驗變差函數理論進行了模型擬合，結合區域地貌地形以及克里金插值結果對該地區的滲透系數進行了分區，取得了如下認識：

(1)伊犁-鞏乃斯河谷地區滲透系數樣本數據經Box-Cox變換后符合正態性分布；而承壓含水層滲透系數樣本經對數變換即可滿足正態性分布。

(2)伊犁-鞏乃斯河谷地區的滲透系數顯示出了強烈的空間相關性，并且該地區滲透系數的空間變異性主要是由巖層的結構本身引起的。

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

(4)根據調查結果，承壓含水層的地層巖性基本表現為從東至西逐漸變細的特點，相應地其滲透系數也在逐漸變小，這與插值結果也比較符合。分區結果顯示承壓含水層滲透系數總體偏小，最大值在R6區，可以達到13.5 m/d；最小值在R1區，僅為1.1 m/d。