內(nèi)蒙古河套灌區(qū)地下水埋深時空變化及其驅(qū)動因素

張文鴿，侯勝玲，殷會娟

(1. 黃河水利委員會黃河水利科學(xué)研究院，鄭州 450003；2. 鄭州大學(xué)水利科學(xué)與工程學(xué)院，鄭州 450001)

地下水作為一種穩(wěn)定的水源且具有地表水所不具備的優(yōu)越性[1]，如分布廣、水質(zhì)好、可持續(xù)開發(fā)利用等，成為干旱地區(qū)最佳的供水選擇[2]，重視地下水合理開發(fā)與污染防治已經(jīng)成為世界各國提高社會與經(jīng)濟(jì)效益的一項重要戰(zhàn)略任務(wù)[3]。因此，地下水埋深時空變化規(guī)律的分析對區(qū)域水資源可持續(xù)利用和地下水開發(fā)有重要的指導(dǎo)作用。

近幾年來，國內(nèi)學(xué)者對我國地下水時空變異規(guī)律開展了諸多研究，其中，地統(tǒng)計學(xué)法是較為常用的方法[4]。張喜風(fēng)等[5]利用遙感和地統(tǒng)計學(xué)方法，分析敦煌綠洲土地利用變化對地下水位時空變異的影響，得出人為因素對地下水位空間異質(zhì)性變化有較大影響并且地下水資源的連通性和脆弱性也在增加。鄧寶山等[6]運用經(jīng)典統(tǒng)計學(xué)方法，探討克里雅綠洲地下水埋深與土壤鹽分的時空變化。喬雨[7]以水文地質(zhì)學(xué)、地下水動力學(xué)、地質(zhì)統(tǒng)計學(xué)等理論為指導(dǎo)，總結(jié)了吉林省中部高平原區(qū)地下水資源時空演化特征。也有一些學(xué)者提供了新思路，王大康[8]利用土地利用轉(zhuǎn)移矩陣及退縮度和擴(kuò)張度模型，從各地類和整體的角度分析綠洲空間演變規(guī)律。賴喬楓等[9]采用ArcGIS空間分析及主成分分析法對大安市地下水動態(tài)變化規(guī)律進(jìn)行分析，運用BP神經(jīng)網(wǎng)絡(luò)對地下水位進(jìn)行預(yù)測。但目前研究主要集中在單井單地區(qū)的多年觀測資料分析，觀測序列也大多為2015年前的數(shù)據(jù)，對多口井和多灌域的研究較少，數(shù)據(jù)實時更新較慢，難以建立更具有時效性和區(qū)域連續(xù)性的時空分布模型。

河套灌區(qū)作為黃河中游的大型灌區(qū)，同時也是中國設(shè)計灌溉面積最大的灌區(qū)，但灌區(qū)自產(chǎn)地表水資源很少，主要是依靠引黃河水[10]，地下水資源是主要供水資源[11]。本文根據(jù)河套灌區(qū)不同灌域的地下水埋深，采用地統(tǒng)計學(xué)、趨勢分析和灰色關(guān)聯(lián)度等方法，研究河套灌區(qū)地下水埋深時空變化及其驅(qū)動因素。

1 資料來源與研究方法

1.1 研究區(qū)概況

河套灌區(qū)位于內(nèi)蒙古自治區(qū)西部巴彥淖爾市境內(nèi)，處于東經(jīng) 106°20′～109°19′、 北緯 40°19′～41°18′之間，研究區(qū)所處位置如圖1所示。河套灌區(qū)南北寬50 km，東西長250 km，總土地面積約112 萬hm2，設(shè)計灌溉面積為73 萬hm2，2017年實際灌溉面積為60.7 萬hm2，其中井渠雙灌面積為7.5 萬hm2。灌區(qū)分為5個灌域：烏蘭布和灌域、義長灌域、永濟(jì)灌域、解放閘灌域和烏拉特灌域。

圖1 研究區(qū)所處位置

1.2 研究方法

本文基于2008-2018年河套灌區(qū)224眼地下水觀測井逐月資料，地下水觀測井位置如圖2所示。運用五點三次平滑法分析灌區(qū)時間序列特征，并利用ArcGIS10.2 軟件，對灌區(qū)地下水位高程及其相應(yīng)的地下水位埋深進(jìn)行克里格( Kriging) 插值并繪制灌區(qū)地下水位等值線圖，從而分析河套灌區(qū)地下水位時空分布特征。

圖2 灌區(qū)地下水監(jiān)測井位置

1.2.1 五點三次平滑

五點三次平滑法有降低通濾波器的作用，以便展現(xiàn)出變量的變化趨勢，與滑動平均相比五點三次平滑的優(yōu)點在于不會削弱過多的波幅，它可以很好地反映序列變化的實際趨勢，特別適合于作相對短時期變化的趨勢的分析。本文利用五點三次平滑分析法對灌區(qū)2008-2018年地下水埋深進(jìn)行時間變化的趨勢分析。

對時間序列x，在其每個數(shù)據(jù)點前后各取兩相鄰數(shù)據(jù)，用三次多項式擬合：

(1)

式中：a、b、c、d可以根據(jù)最小二乘法原理確定，以此確定五點三次平滑計算公式。

五點三次平滑公式如下：

(2)

(3)

(4)

(5)

(6)

對序列的開始兩點用公式(2)和(3)平滑，最后兩點用公式(5)和(6)進(jìn)行平滑，其余各點均按式(4)進(jìn)行平滑。

1.2.2 灰色關(guān)聯(lián)度

灰色關(guān)聯(lián)度分析，是20世紀(jì)80年代鄧聚龍教授創(chuàng)立的灰色系統(tǒng)理論中的重要內(nèi)容之一[12]。通過關(guān)聯(lián)度表征地下水系統(tǒng)各因素的密切相關(guān)程度，從而決定影響系統(tǒng)的主要因素，為進(jìn)一步進(jìn)行系統(tǒng)分析提供指導(dǎo)[13]。

本文通過灰色關(guān)聯(lián)度法計算降水量、蒸發(fā)量和引黃水量三種影響因素與地下水位埋深的關(guān)聯(lián)度，分析地下水位埋深主要驅(qū)動因素[14]。關(guān)聯(lián)度分析的計算和步驟如下：

設(shè)k(k=1,2,…,n)為時間編號序列；Xi(k)為驅(qū)動因子Xi關(guān)于第k年的觀測數(shù)據(jù)，則得到驅(qū)動因子Xi的橫向序列為Xi(k)|k=1,2,…,n，(i=1,2,…,m)，m為驅(qū)動因子的個數(shù)。

(1)原始數(shù)據(jù)變換。消除原始數(shù)據(jù)的量綱，轉(zhuǎn)換為可比較的數(shù)據(jù)序列，可采用均值化變換、初值化變換和標(biāo)準(zhǔn)化變換，本文采用均值化變換。

(2)求差序列。

Δi(k)=|L′(k)-X′(k)|

(7)

Δi=Δi(k)|k=1,2,…,n(i=1,2,…,m)

(8)

式中：X′(k)為Xi(k)的均值；L′(k)為時間k的折減系數(shù)均值。

(3)求關(guān)聯(lián)系數(shù)。

(9)

式中：η為分辨系數(shù)，其作用在于提高關(guān)聯(lián)系數(shù)間的差異顯著性，η∈(0,1)，一般取η=0.5。

(4)計算關(guān)聯(lián)度。

(10)

1.2.3 克里格( Kriging) 插值法

克里格(Kriging)插值又稱空間局部插值法，在有限區(qū)域內(nèi)對區(qū)域化變量進(jìn)行無偏最優(yōu)預(yù)估的一種方法，是地統(tǒng)計學(xué)的主要內(nèi)容之一[15]。法國著名統(tǒng)計學(xué)家G Matheron將該方法理論系統(tǒng)化，并命名為 Kriging，即克里格方法[16]。利用此方法對灌區(qū)地下水位埋深進(jìn)行空間插值并分析灌區(qū)地下水埋深空間分布特點。

克里格插值法的公式如下：

(11)

式中：Z(x0)為未知樣點的值；Z(xi)為未知樣點周圍的已知樣本點的值；ωi為第i個已知樣本點對未知樣點的權(quán)重；n為已知樣本點的個數(shù)。

地質(zhì)統(tǒng)計學(xué)方法的理論基礎(chǔ)是半方差函數(shù)，它是描述隨機(jī)變量空間變異結(jié)構(gòu)的一個函數(shù)。半方差函數(shù)的計算式為：

(12)

式中：r(h)為變異函數(shù)；h為分離間距；Nh是間距為h的樣本“對”數(shù)，它的下標(biāo)h表示Nh是分離距離的函數(shù)。

樣本變異函數(shù)僅僅是一個數(shù)據(jù)的概括技術(shù)，還需借助理論模型來擬合變異函數(shù)曲線。常用的理論模型為有基臺值模型，其中包括球狀模型、高斯模型和指數(shù)模型等，其一般公式見表1。

2 結(jié)果與分析

2.1 灌區(qū)地下水埋深時間變化

2.1.1 年際變化

河套灌區(qū)現(xiàn)有地下水位觀測孔224個，觀測系列較長，選取灌區(qū)2008-2018年烏蘭布和灌域、義長灌域、解放閘灌域、永濟(jì)灌域和烏拉特灌域5個灌域地下水埋深數(shù)據(jù)，利用五點三次平滑法對灌區(qū)和各灌域地下水動態(tài)規(guī)律進(jìn)行分析，并繪制各灌域和灌區(qū)多年地下水位變化過程見圖3和圖4。

表1 理論模型的一般公式

圖3 河套灌區(qū)各灌域2008-2018年地下水埋深變化趨勢

由圖3可見，各灌域的地下水位變幅有逐年下降的趨向，因為灌溉面積、人口規(guī)模、取水條件和經(jīng)濟(jì)發(fā)展水平的差異，各灌域地下水埋深差別較大。2008-2010年之間，各灌域的年變幅在1.92 m左右，到2011年各灌域地下水位大幅度降落。其中，烏蘭布和灌域地下水埋深最淺，多年平均值為1.82 m；義長灌域地下水埋深最大，多年平均值為2.08 m，兩者相差0.27 m。永濟(jì)、解放閘和烏蘭布和灌域地下水位年變化幅度較小，義長灌域與烏拉特的地下水位變化幅度最大。

圖4 灌區(qū) 2008-2018年地下水埋深變化與平滑趨勢

地下水動態(tài)受到降水、蒸發(fā)和融凍等氣象因素的影響，表現(xiàn)為季節(jié)性變化，受引黃灌溉水影響表現(xiàn)為周期變換[17]。從內(nèi)蒙古河套灌區(qū)總體來看，地下水埋深的年際動態(tài)為灌溉周期，基本處于一個波動狀態(tài)，如圖4所示，河套灌區(qū)2008-2018年的年平均地下水埋深歷年變化幅度較大，年均地下水埋深在1.85～2.25 m之間變動。2011-2013年，灌區(qū)地下水埋深有小幅度的上升趨勢，但總體來看，河套灌區(qū)地下水埋深年際變化具有人類活動因素所造成的開采型下降趨勢[11]。

2.1.2 年內(nèi)變化

分析灌區(qū)2008-2018年全灌區(qū)的月平均地下水埋深數(shù)據(jù)，灌區(qū)月均地下水位變化趨勢見圖5。

圖5 灌區(qū) 2008-2018年地下水埋深月均變化

由于河套灌區(qū)地下水埋深的年內(nèi)變化受凍融、降雨和蒸發(fā)等氣象條件的影響，年內(nèi)動態(tài)變化表現(xiàn)為季節(jié)性變化。每年4月下旬開始春灌，灌后水位大幅回升，升幅在0.5 m以上；5月上旬至8月下旬的作物生長需水期，潛水平均埋深在1.8 m上下，蒸發(fā)排泄水量較多；9月上旬夏灌結(jié)束，地下水位逐漸下降，直至9月中旬冬灌前，埋深下降到2.5 m左右；9月下旬冬灌開始，水位再次回升至接近地表，見圖5。

2.2 灌區(qū)地下水埋深空間分布

2.2.1 地下水埋深的插值分析

為了準(zhǔn)確直觀地對河套灌區(qū)地下水埋深的時空分布進(jìn)行分析，選取2008、2010、2012、2014、2016和2018年地下水?dāng)?shù)據(jù)，運用克里格插值繪制出河套灌區(qū)各年地下水埋深的空間分布圖(見圖6)，直觀反映了研究區(qū)地下水埋深空間分布及其變異特征、變異程度。

圖6 灌區(qū)各年地下水埋深空間分布圖

從圖6能夠看出，整體上，河套灌區(qū)的地下水位埋深都較淺，主要是因為關(guān)于地下水開采量較小[10]。地下水埋深>10 m的區(qū)域主要分布在義長灌域東北部以及城鎮(zhèn)密集的區(qū)域；埋深4～10 m的區(qū)域主要分布在解放閘灌域北部和南部以及永濟(jì)灌域北部，區(qū)域面積變化較為顯著；埋深<4 m的區(qū)域面積較大，并且是連貫的分布，大致分布在灌區(qū)中部。2008-2012年，地下水埋深為2.05～4.10 m的區(qū)域面積開始減少，大部分轉(zhuǎn)變?yōu)?.03～5.24 m。隨著地下水的不斷開發(fā)使用，灌區(qū)地下水埋深逐步加深，連貫區(qū)域被逐漸增大的深埋深區(qū)域割裂為不連貫的小塊。從各灌域來看，烏蘭布和灌域大部分區(qū)域為沙區(qū)，地下水開采較為困難，烏蘭布和各年地下水埋深下降趨勢較不明顯；解放閘灌域地下水的主要補(bǔ)給是灌溉水和降水入滲[18]，潛水蒸發(fā)是地下水的主要排泄出路，在降雨與蒸發(fā)的影響下，解放閘地下水開采量增加，地下水埋深一直下降；永濟(jì)灌域位于河套灌區(qū)中間，無越流補(bǔ)給發(fā)生，地下水補(bǔ)給受地表灌溉和降雨等要素影響，該灌域地下水埋深變化較大；義長灌域處于河套灌區(qū)中部位置，是灌區(qū)生物活動頻繁之處，且受東北低西南高的地勢影響，導(dǎo)致南部地下水埋深年際變化相對平緩。在人為活動與地勢的影響下，在空間上該灌域多年地下水埋深呈現(xiàn)連續(xù)性分布；烏拉特灌域距三盛公樞紐較遠(yuǎn)，引水條件較差，地下水用量逐年增加，導(dǎo)致地下水埋深下降。

綜上所述，灌區(qū)地下水埋深較淺但逐年下降，考慮到河套灌區(qū)的輸配水效率、灌溉節(jié)水能力和水資源量等條件，地下水位可能面臨超采風(fēng)險。在灌區(qū)各年地下水埋深分布圖的基礎(chǔ)上，利用ArcGIS軟件中的重分類將地下水埋深面積分為<2 m、2～5 m和>5 m三類，計算不同地下水埋深所占面積，得到灌區(qū)不同地下水埋深所占面積統(tǒng)計圖(圖7)。

圖7 河套灌區(qū)不同地下水埋深所占面積統(tǒng)計圖

2.2.2 地下水埋深的空間變異分析

由圖7可知，地下水埋深<2 m的區(qū)域在2008年、2010年以及2012年占比較大，分別為64.8%、59.1%和62.3%，然而從2014年開始減少至53.0%。隨著埋深<2 m的區(qū)域面積減少，埋深在2～5 m的區(qū)域面積開始不斷增加，從2008年埋深2～5 m區(qū)域面積占比為32.4%上升到2018年的54.8%。埋深范圍>5 m的區(qū)域面積變化較小，2008-2018年面積占比圍繞3.0%波動。對比2008-2018年灌區(qū)地下水埋深數(shù)據(jù)，埋深<2 m的區(qū)域不斷向埋深2～5 m區(qū)域轉(zhuǎn)化，但>5 m地下水埋深無顯著變化，說明灌區(qū)淺層地下水被不斷開采而較深地下水還未受到影響。當(dāng)灌區(qū)超過一半面積為2 m以上的埋深時，灌區(qū)種植的作物和自然植被的生長就會受到嚴(yán)重影響[3]，2016年與2018年灌區(qū)地下水埋深面積已有超過一半為2 m以上，這一現(xiàn)象應(yīng)當(dāng)引起相關(guān)部門的重視，盡早采取行動，以防地下水超采。

2.3 地下水埋深的影響因素

驅(qū)動河套灌區(qū)地下水位變化的因素主要有降雨量、蒸發(fā)量和引黃水量。地下水位埋深及主要影響因素觀測資料見表2，采用灰色關(guān)聯(lián)度對灌區(qū)2008-2018年均地下水位變化影響因素進(jìn)行分析，結(jié)果見表3。并繪制2008-2018年蒸發(fā)量、降雨量、引黃水量和地下水埋深變化情況見圖8。

表2 地下水位埋深及主要影響因素觀測資料

由表3可以得到：2008-2018年均地下水埋深變化驅(qū)動因素的灰色關(guān)聯(lián)度分別為：降水量0.67，蒸發(fā)量0.78，引黃水量0.75，因此，3個主要因素與地下水位關(guān)聯(lián)度由大到小排列為：蒸發(fā)量(0.78)>引黃水量(0.75)>降水量(0.67)?？梢钥闯?，蒸發(fā)量和引黃水量對灌區(qū)年均地下水位埋深影響最大。主要因為河套灌區(qū)位于我國干旱區(qū)，蒸發(fā)量大，降水量少，灌區(qū)采取引黃灌溉，地下水減少受灌排和蒸發(fā)的影響，主要靠灌溉入滲和降水補(bǔ)給。

表3 各相關(guān)影響因素關(guān)聯(lián)度結(jié)果

由圖8可以看出灌區(qū)蒸發(fā)量對地下水位埋深影響較為顯著，在2009-2011與2015-2017之間，灌區(qū)蒸發(fā)量處于較高水平，相應(yīng)地下水位埋深也呈現(xiàn)出下降趨勢。2012年6月26日開始，河套灌區(qū)連續(xù)降雨3天，據(jù)氣象、水文部門測定為年一遇特大暴雨，降雨幾乎覆蓋了整個灌區(qū)[19]，灌區(qū)降水量達(dá)到287.60 mm，從而使得當(dāng)年灌區(qū)蒸發(fā)量和引黃水量下降，同時地下水埋深也上升至1.87 m。引黃水量作為灌區(qū)重要的水資源來源，對地下水進(jìn)行補(bǔ)給，影響地下水位埋深，這是因為河套灌區(qū)具有自產(chǎn)水少，過境水多的特點，境內(nèi)唯一過境水資源為黃河。

圖8 各驅(qū)動因子與地下水埋深變化情況圖

3 結(jié) 語

本文以內(nèi)蒙古河套灌區(qū)2008-2018年224 眼地下水觀測井逐月資料為基礎(chǔ)，運用五點三次平滑法分析灌區(qū)地下水年際年內(nèi)趨勢變化，利用Arcgis軟件對灌區(qū)地下水埋深做插值分析并繪制地下水埋深圖和不同埋深所占面積統(tǒng)計圖，并結(jié)合灰色關(guān)聯(lián)度探究地下水埋深的主要驅(qū)動因素。根據(jù)趨勢圖與數(shù)據(jù)計算結(jié)果，灌區(qū)整體地下水埋深較淺，但在人為活動與自然條件影響下，灌區(qū)地下水開采量逐年增加，造成地下水埋深逐年下降，其中義長灌域位于灌區(qū)中部，由于植被分布密集、人類活動頻繁，地下水埋深變化最為明顯。從關(guān)聯(lián)度的角度來看，蒸發(fā)量和引黃水量對灌區(qū)地下水埋深變化影響較大。河套灌區(qū)空間變異大，政府部門需制定相關(guān)政策達(dá)到合理開發(fā)地下水資源，有效提高灌區(qū)水資源的利用率，提高各灌域地下水開采治理力度，保護(hù)灌區(qū)水生態(tài)環(huán)境的目的。