張文鴿,侯勝玲,殷會娟
(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ū)動因素。
河套灌區(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ū)所處位置
本文基于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.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.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)部門的重視,盡早采取行動,以防地下水超采。
驅(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ū)動因子與地下水埋深變化情況圖
本文以內(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)境的目的。