塔里木河下游地下水位空間變異性研究

田鼎文

(巴州新宇水利水電工程咨詢有限責任公司， 新疆 庫爾勒　841000)

?

塔里木河下游地下水位空間變異性研究

田鼎文

(巴州新宇水利水電工程咨詢有限責任公司， 新疆 庫爾勒841000)

【摘要】塔里木河下游干旱缺水，地下水資源是生態健康的基本保障。本文采用線性回歸法對塔里木河下游地下水位進行回歸分析，并采用地統計學方法對當地地下水位的空間變異規律進行分析。結果表明：塔里木河下游各月地下水位與其測點空間坐標呈線性關系，具有很強的空間相關性； 對于局部變化劇烈的地下水位，采用克里金平滑效應修正法可明顯減小克里金插值所帶來的平滑效應，但整體插值精度略低于普通克里金與泛克里金法。

【關鍵詞】塔里木河下游； 地下水位； 空間變異； 研究

干旱區水資源短缺，地下水資源常常成為干旱區綠洲的重要供水水源[1]。地下水不恰當開采引起的地下水位下降和地下水水質惡化是干旱區綠洲水資源利用中的突出問題，會對綠洲農業生產及生態保護產生不利影響。塔里木河是中國第一大內陸河，地下水資源是當地綠洲生態健康的重要保障，了解流域內地下水位的空間分布規律及其動態變化，對于地下水資源的可持續開發利用具有重要意義。

地下水位是具有隨機性與空間相關性雙重屬性的隨機化區域變量，預測地下水位空間分布的統計學方法有回歸分析法、反距離加權插值法、多項式插值法等。陳志宏[2]運用多元線性回歸分析法對地下水位的變化規律及其與相關因素的關系進行研究，發現多元線性回歸分析法可以充分利用觀測資料，反映地下水位變化規律。在其他眾多研究地下水位空間插值的方法中，大都認為克里金插值法預測地下水位空間分布效果最優。例如：李金榮等[3]在地下水位插值計算方法的比較中，認為克里金法優于趨勢面法。Yao等[4]以西北干旱區石羊河流域為例，對比八種空間插值方法，發現普通克里金法優于其他方法，能更真實反映地下水位變化趨勢。Varouchakis[5]在預測希臘某流域地下水位時，也發現克里金法優于反距離加權法和最小曲率法。

預測地下水位空間分布的統計學方法有回歸分析法、反距離加權插值法、多項式插值法等。本文采用線性回歸分析法，對塔里木河下游2000—2013年間地下水位的空間變化趨勢進行分析，并采用地統計學中幾種克里金插值法對地下水位進行插值分析，對比分析普通克里金法、克里金平滑效應修正法以及泛克里金法在該地區地下水位空間估值中的適用性，探討該地區地下水位的空間變異規律及其動態變化。

1計算模型選用

1.1模型簡介

回歸分析法是研究兩個或多個隨機變量間關聯性的方法，利用回歸分析法導出的因變量與自變量之間的數學關系式稱為回歸方程或回歸模型。塔里木河下游地下水位線性回歸模型如下：

(1)

式中Hi——第i個月份的地下水位，m；

X，Y——測點的投影直角坐標，m(以下稱為空間坐標)；

ai，bi，ci——回歸系數，回歸系數a或b大于0，說明地下水位與空間坐標X或Y呈正相關，反之，則呈負相關。

回歸方程的擬合優度檢驗可以采用判定系數R2統計量，判定系數R2越接近于1，說明回歸方程對地下水位樣本數據點的擬合優度越高，該月份的地下水位與空間坐標的線性回歸關系較強；反之，越接近于0，說明回歸方程對地下水位的擬合優度越低，該月份的地下水位與空間坐標的線性關系較弱。

1.2變異函數

變異函數γ(h)是區域化變量空間結構的一種形式化表達，是地統計學理論的基礎，其數學表示為兩個隨機變量Z(x)和Z(x+h)之間增量的方差的一半。在滿足二階平穩性假設下(期望值不變，協方差或變異函數僅與空間間隔距離h相關)，原始變異函數γx(h)為單純的空間間隔距離h的函數γ(h)，而與空間位置無關，即

(2)

常用的變異函數模型有指數模型、球形模型、高斯模型等。變異函數有三個重要參數：變程、塊金值、基臺值。變異函數會隨間隔距離h的增大而增大，最后達到一個穩定的常數，該穩定的常數稱為基臺值；此時的距離h為變程；當距離為0時，變異函數值為塊金值。基底效應為塊金值與基臺值的比值，表示空間變異程度。

本文的區域化變量為地下水位，當觀測點之間的距離大于變程時，表示研究區地下水位的空間相關性消失。

1.3克里金法

1.3.1普通克里金法

在滿足二階平穩假設時，地下水位的普通克里金(Kriging)法估計公式為:[6]

(3)

式中Z*(x0)——x0位置的地下水位估計值；

Z(xi)——xi位置的地下水測量值；

λi——分配給Z(xi)的權重；

n——用于估計過程的測量值的個數。

普通克里金法估計值誤差的方差為

(4)

在求估計方差的極小值時須引入拉格朗日乘數μ，使估計方差最小，得到了計算權重λi的線性方程組：

(5)

求解式(5)，得到拉格朗日乘數μ和權重λi，將權重λi代入式(3)，即可求得地下水位的普通克里金估計值。此時，估計值誤差的方差計算式為

(6)

1.3.2普通克里金插值平滑效應的Yamamoto修正法

為克服普通克里金插值時產生的平滑效應，Yamamoto提出了一套對普通克里金估計值進行后處理的方法，可以較好地修正平滑效應的影響。普通克里金估計值在經過修正后，還必須檢查修正后的估計值是否落在給定范圍內實測點的數值范圍(Zmin，Zmax)內。如果超出這個范圍，還需要對修正后的估計值做二次修正。

若Z**(x0)>Zmax或Z**(x0)

Z**(x0)=Z**(x0)+δ×f

(7)

式中Z**(x0)——修正后的普通克里金估計值；

Ns0——插值標準差數，為負的普通克里金估計誤差與插值標準差(插值方差的平方根)的比值；當Ns0<0時，δ=-[Z*(x0)-Zmin]；當Ns0>0時，δ=Zmax-Z*(x0)；

f——優化因子，以保證修正后值的方差盡可能接近樣本的方差。

將經過上述二次修正的估計值與所有的實測樣本點進行對比，檢查該值是否落在全部樣本點數值范圍(Smin，Smax)內。若超出該范圍，將再次進行修正，即

(8)

上式保證了修正后普通克里金估計值的平均值與樣本的平均值盡量保持一致。

1.3.3泛克里金法

當區域變量不滿足二階平穩假設時，一般情況下空間變量的均值不再是一個常數，而是一個空間變量。假定在位置x上的區域變量z(x)可用確定性漂移m(x)和殘差部分r(x)來代表，即：

(9)

通過漂移的定義，z(x)在x的期望為m(x)，即:

(10)

假定漂移m(x)可以用下式表示：

(11)

式中ak——未知系數；

pk(x)——已知的x的函數，可以用xk-1的形式來表示。

假定殘差部分r(x)滿足二階平衡條件，用σ(h)來表示它的協方差函數，如果xx′是Ω中的兩點x和x′間的距離，r(x)和r(x′)的協方差只是距離xx′的函數：

(12)

假定σ(h)是已知的，泛克里金法的主要問題是確定漂移部分的最佳系數ak。漂移的估計值是樣本值z(xi)的線性組合：

(13)

式中，權重bi是xi的函數，根據普遍性、最佳性選擇權重，從而得到漂移與殘差的估計值。

泛克里金最終估計式為：

(14)

2下游地下水位線性回歸分析

利用測點的空間坐標對塔里木河下游2010—2013年各月地下水位進行二元線性回歸分析。結果表明：2010—2013年，線性回歸方程的判定系數均在0.80左右，在此期間，地下水位與空間坐標有良好的線性關系。根據已有研究，塔里木河下游的農作物覆蓋度自2000年之后相對穩定，表明塔里木河下游地下水位只有在植被覆蓋穩定時期與空間坐標有較好的線性回歸關系。

在2010—2013年，每年4月或5月以及10月或11月的判定系數小于0.50。說明在此時間段，每年這兩個月份的地下水位與其空間坐標無良好的線性關系。4月、5月為農田作物的播種期，10月、11月為農田作物的收獲期，這期間地下水開采頻繁且沒有規律，導致局部地下水位嚴重下降，間接擾亂了地下水位與空間坐標的線性關系。然而，影響地下水位變化的因素除空間坐標外，還包括降雨量、地下水開采量等很多復雜因素，在研究地下水位變化時還應考慮其他因素的影響。

3樣本變異函數隨時間變化的特點

根據地統計學理論和方法，利用塔里木河下游2010—2013年各月地下水位計算它們的樣本方差函數值，做出實測樣本變異函數圖，對實測樣本變異函數進行擬合，并采用交叉驗證計算各擬合模型的方差與相關系數，計算結果顯示高斯模型擬合效果最好。采用高斯模型對各月地下水位的變異函數進行擬合，發現其時間變程均在1400km以下，變異函數變程可以反映區域化變量的影響范圍，可知本研究區地下水位最大相關距離為1400km。塊金值與基臺值比值可以反映空間變異程度，該比值高，說明由隨機部分引起的空間變異性程度較大，相反則由結構性因素引起的空間變異程度較大。由區域化變量空間相關程度的分級標準(小于25%為強烈空間相關)可知：本研究區在2010—2013年間地下水位的空間隨機性很小，具有很強的空間相關性。

4地下水位基本統計量

由于2010年地下水位觀測數據較全面，因此本文對2010年各月地下水位進行基本統計分析(見下表)。由表中可以看出：在不同月份，地下水位的最小值年內波動較小。樣本(地下水位)各月變異系數均在7%以下，按照變異系數Cv的等級劃分(弱變異性，Cv<10%；中等變異性，10%100%)可知本研究區內地下水位樣本各月的變異性較弱。

5克里金法在塔里木河下游地下水位空間插值的應用比較

利用塔里木河下游2010年觀測井的年均地下水位觀測數據，分別采用普通克里金插值法、普通克里金Yamamoto修正法以及泛克里金插值法，對研究區地下水位進行空間插值預測。發現三種方法的所有預測結果與實測值對比均方根誤差(RMSE)相差不大，其中普通克里金方法最小，Yamamoto修正法最大；三種方法的相關系數(R2)均接近于1，表明用三種方法預測的地下水位值與實測水位之間具有良好的一致性和相關性。

為進一步驗證Yamamoto修正法對于平滑效應的修正程度，選取實測地下水位極端值與其預測值進行對比，三種方法的均方根誤差(RMSE)以Yamamoto修正法最小，三種插值方法得到的研究區地下水位空間分布趨勢大致相同，但在一些局部區域，使用Yamamoto修正法修正后的插值結果比修正前的值小，而在另一些區域修正后的值會變大，這在一定程度上消除了克里金法的平滑效應。總體分析，三種方法均可應用于塔里木河下游地下水位的空間預測，普通克里金法可以保證整體的插值精度；而在側重于局部較高或較低地下水位研究時，Yamamoto修正法可以消除平滑效應，更接近于極端水位的真實值。

6結語

本文采用多元線性回歸分析法對塔里木河下游地下水位與其空間坐標進行回歸分析，并根據地統計學中的變異函數與克里金插值法對塔里木河下游地下水位分布空間變異性進行了分析。結果表明：在植被覆蓋穩定時期各月地下水位與其測點空間坐標呈顯著的線性相關，而地下水無規律的頻繁開采可能擾亂此線性關系；地下水位分布的變異函數模型選用高斯模型較為合適，各個月份地下水位均具有很強的空間相關性。對塔里木河下游進行地下水位空間分布插值預測時，Yamamoto修正法可以消除克里金插值所帶來的平滑效應，在地下水位極高值與極低值點的預測上，可以根據實測點數據較真實地重現地下水位的空間分布，但整體插值精度要略小于普通克里金法與泛克里金法。

參考文獻

[1]張鴻義.中國干旱區地下水資源及開發潛力分析[J].干旱區研究,2009,26(3):149-161.

[2]陳志宏.多元線性回歸方法在地下水位預測中的應用[J].北京地質,1999(3):20-26.

[3]李金榮,楊振放,李云峰.兩種方法在地下水位估值中的應用[J].水文地質工程地質,2003(3):42- 46.

[4]YAO L Q, HUO Z L, FENG S Y, et al. Evaluation of spatial interpolation methods for groundwater level in an arid inland oasis, northwest China[J]. Environmental earth sciences，2014，71(4): 1911-1924.

[5]VAROUCHAKIS E A, Hristopulos D T.Comparison of stochastic and deterministic methods for mapping groundwater level spatial variability in sparsely monitored basins[J]. Environmental Monitoring and Assessment, 2013,185(1):1-19.

[6]JANG C S, CHEN S K, LIN C C. Using multiple-variable indicator kriging to assess groundwater quality for irrigation in the aquifers of the Choushui River alluvial fan[J]. Hydrological Processes, 2008, 22(22): 4477- 4489.

Research on underground water level spatial variability on the lower reaches of the Tarim River

TIAN Dingwen

(BazhouXinyuWaterConservancyandHydropowerEngineeringConsultingCo.,Ltd.,Korla841000,China)

Abstract：The lower reaches of the Tarim River is dry and lack of water. Groundwater resources are the basic guarantee of ecological health. In the paper, linear regression method is adopted for regression analysis on underground water level on lower reaches of the Tarim River. Statistical method is adopted to analyze spatial variation law of local underground water level. The results show that underground water level on the lower reaches of the Tarim River has linear relation with spatial coordinates at the measuring points each month, which has strong spatial correlation. Kriging smoothing effect correction method is adopted for prominently reducing the smoothing effect brought by Kriging interpolation aiming at underground water level with sharp partial change. However, overall interpolation precision is slightly lower than ordinary Kriging and universal Kriging method.

Key words：the lower reaches of the Tarim River; underground water level; spatial variation; research

DOI:10.16616/j.cnki.10-1326/TV.2016.02.020

中圖分類號：TV211

文獻標志碼：A

文章編號：1005- 4774(2016)02- 0070- 05