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

田鼎文

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

?

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

田鼎文

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

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

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

干旱區(qū)水資源短缺，地下水資源常常成為干旱區(qū)綠洲的重要供水水源[1]。地下水不恰當(dāng)開采引起的地下水位下降和地下水水質(zhì)惡化是干旱區(qū)綠洲水資源利用中的突出問題，會對綠洲農(nóng)業(yè)生產(chǎn)及生態(tài)保護產(chǎn)生不利影響。塔里木河是中國第一大內(nèi)陸河，地下水資源是當(dāng)?shù)鼐G洲生態(tài)健康的重要保障，了解流域內(nèi)地下水位的空間分布規(guī)律及其動態(tài)變化，對于地下水資源的可持續(xù)開發(fā)利用具有重要意義。

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

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

1計算模型選用

1.1模型簡介

回歸分析法是研究兩個或多個隨機變量間關(guān)聯(lián)性的方法，利用回歸分析法導(dǎo)出的因變量與自變量之間的數(shù)學(xué)關(guān)系式稱為回歸方程或回歸模型。塔里木河下游地下水位線性回歸模型如下：

(1)

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

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

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

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

1.2變異函數(shù)

變異函數(shù)γ(h)是區(qū)域化變量空間結(jié)構(gòu)的一種形式化表達，是地統(tǒng)計學(xué)理論的基礎(chǔ)，其數(shù)學(xué)表示為兩個隨機變量Z(x)和Z(x+h)之間增量的方差的一半。在滿足二階平穩(wěn)性假設(shè)下(期望值不變，協(xié)方差或變異函數(shù)僅與空間間隔距離h相關(guān))，原始變異函數(shù)γx(h)為單純的空間間隔距離h的函數(shù)γ(h)，而與空間位置無關(guān)，即

(2)

常用的變異函數(shù)模型有指數(shù)模型、球形模型、高斯模型等。變異函數(shù)有三個重要參數(shù)：變程、塊金值、基臺值。變異函數(shù)會隨間隔距離h的增大而增大，最后達到一個穩(wěn)定的常數(shù)，該穩(wěn)定的常數(shù)稱為基臺值；此時的距離h為變程；當(dāng)距離為0時，變異函數(shù)值為塊金值?；仔?yīng)為塊金值與基臺值的比值，表示空間變異程度。

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

1.3克里金法

1.3.1普通克里金法

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

(3)

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

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

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

n——用于估計過程的測量值的個數(shù)。

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

(4)

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

(5)

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

(6)

1.3.2普通克里金插值平滑效應(yīng)的Yamamoto修正法

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

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

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

(7)

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

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

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

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

(8)

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

1.3.3泛克里金法

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

(9)

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

(10)

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

(11)

式中ak——未知系數(shù)；

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

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

(12)

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

(13)

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

泛克里金最終估計式為：

(14)

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

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

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

3樣本變異函數(shù)隨時間變化的特點

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

4地下水位基本統(tǒng)計量

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

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

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

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

6結(jié)語

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

參考文獻

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

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

[3]李金榮,楊振放,李云峰.兩種方法在地下水位估值中的應(yīng)用[J].水文地質(zhì)工程地質(zhì),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