楊瀅嘉, 張丹蓉, 莊會波, 劉 江, 管儀慶, 邵廣文
(1.河海大學 水文水資源學院,江蘇 南京 210098;2.山東省水文局,山東 濟南250000)
水資源是社會經濟發展和人民生產生活的主要影響因素,地下水作為水資源的重要組分,具有水質良好、水量穩定、調節力強、便于開采等優點,已成為城鎮發展建設、人民飲水安全、生態環境改善的重要保障[1]。針對目前很多地區出現地下水開采利用過度的現象,研究地下水變化規律和分布特征,對了解區域地下水量、制定水資源管理方針、合理開采利用地下水至關重要[2-3]。
地下水觀測數據是進行地下水運動研究和模擬的基礎,主要是對地下水監測井進行水位、水質狀況的觀測。地下水監測井是一群空間離散分布點,因此需要在離散型數據的基礎上采用空間插值法來探討地下水埋深時空變化特征,其中地統計學中的Kriging插值法是普遍使用的一種內插法。Machiwal等[4]將地統計學與GIS相結合,利用50眼監測井數據,研究了印度Ahar河流域36個月間地下水位的時空變化;Ma等[5]采用普通Kriging法模擬了Kansas中南部地下水位變化;Sun等[6]采用簡單Kriging法,分析了民勤綠洲48眼監測井過去22年的地下水埋深狀況,揭示了其變化規律及影響因素。地下水開采量過大的地區易形成漏斗區,引發地面沉降、咸水入侵、水質變差等一系列環境惡化的問題,明確地下水開采現狀及超采區劃分非常重要,李文體等[7]基于地下水動態資料,對河北平原區的地下水開采區進行不同開采程度的劃分,分析了超采情況;黃曉燕等[8]運用水位動態法等多種方法,開展江蘇省超采區評價并進行對比分析;陳國浩[9]運用水位動態法,分析了濟寧市236眼監測井水位動態變化速率,對淺層孔隙水進行超采區劃分,提出了相應整治對策。目前對地下水埋深時空變化特征的研究較多,而在此基礎上科學合理地評價地下水開采適宜度在很多區域尚未實現。
山東省濰坊市人均水資源量僅為全國人均的15.9%,是我國極度資源型缺水城市之一[10]。濰坊中部平原高產糧食作物,多數傳統農業灌溉嚴重消耗水資源,隨著城市化建設、居民生活、服務業需水量不斷增加,過度開采利用地下水導致流域水循環的失衡,危及流域的生態環境健康,目前已有一些學者對濰坊市超采現狀和降落漏斗成因進行了分析[11-13]。本文根據濰坊市中部彌河流域山前平原區1985-2015年地下水觀測井五日觀測資料,分別對地下水埋深年際年內變化趨勢、埋深變幅速率、不同埋深段面積變化特征進行分析,淺析其影響因素,并劃分流域超采區,評價地下水開采適宜性。研究成果有助于了解彌河流域平原區地下水埋深變化情況,對有效控制局部降落漏斗及合理開采利用地下水資源有重要的參考價值。
彌河是山東省15條骨干河道之一,發源于山東中部魯山和沂山之間,自南向北流經臨朐、青州、壽光后入渤海,彌河全長206 km,流域面積3 847.5 km2,如圖1所示。
研究區處于魯中以北的沿海平原區,坐標位置為北緯118°22′~119°10′,東經36°27′~37°19′,總面積為3 387.2 km2,屬于暖溫帶季風區大陸性氣候,多年平均降雨量為561.5 mm,年均蒸發量為1 003.6 mm,年均氣溫11.7~14.3℃,年日照時數在2 450 h左右,全年多盛行南偏東風。地勢自南向北減小(圖1(a)),地下水含水巖組主要是松散巖類第四系孔隙含水巖組,巖性以中粗砂、礫石為主,區域總體富水性較強(圖1(b)),主要補給來源有降水補給、回灌補給、側向徑流補給等。
收集整理所需的基本資料,本文所采用的數據信息如表1所示。

圖1 彌河流域地理位置、高程分布及水文地質圖

表1 研究所需基本資料數據信息及來源
Mann-Kendall是進行時間序列趨勢分析的一種常用的非參數檢驗法,以使用約束少、人為干擾小、適用范圍寬等優點被廣泛使用,其計算結果可信度高[14]。在趨勢檢驗中,統計變量Z>0,表示序列呈上升趨勢;Z<0,表示呈下降趨勢。其中|Z|>1.96,說明上升或下降趨勢顯著。
相關性分析法是對兩個或多個存在關聯的變量因子進行分析計算,通過相關系數來衡量變量因子之間的密切度[15-16]。
本文使用該方法計算地下水位數據及其影響因子之間的相關系數,分析影響地下水動態變化的因素,相關系數對應相關程度如表2。
Kriging法是ArcGIS地統計學中的一種對空間未知樣點的區域化變量的取值進行線性無偏最優內插的方法,該方法最大程度地利用了已知樣點包括的所有信息[17-18]。其基本思路是首先通過分析區域化變量在空間位置上的分布,選擇最佳理論變異函數模型,其次確定對未知樣點進行插值的影響范圍,在范圍內根據已知樣點與未知樣點的空間位置和相關性的不同,賦予每個已知樣點相應的權重,進行加權平均后估計出待測樣點屬性值。本文對數據的空間趨勢走向和殘差項綜合分析后,采用Geostatistical Analyst工具中的普通Kriging插值法對研究區56眼監測井地下水位值進行空間插值,轉為90 m分辨率的地下水位柵格圖,將90 m分辨率的DEM與同分辨率的水位柵格圖相減,得到1985、1995、2002、2009、2012、2015年地下水埋深分布圖。
水位動態法是指以地下水水位年均下降速率和現狀年的地下水埋深作為判斷指標進行超采區劃分的方法[19]。柵格計算是對多個柵格數據進行空間函數分析運算;重分類是對原有柵格數據按照新的分類標準重整輸出一組新值[20]。本文利用ArcGIS中的重分類和柵格計算兩個功能實現對研究區地下水是否超采的劃分。
首先利用ArcGIS中的重分類將研究區2015年地下水埋深空間分布圖分為3類,即埋深<6 m,賦值1;6 m≤埋深≤25 m,賦值2;埋深>25 m,賦值3;
其次利用柵格計算器計算區域研究期內地下水埋深年均變化速率,重分為3類,即小于0.2 m/a的區域,賦值4(非超采區);0.2~1 m/a的區域,賦值5(一般超采區);大于1 m/a的區域,賦值6(嚴重超采區);
最后對前兩個重分類得到的新柵格數據進行計算(2015年地下水埋深重分類數據×地下水開采區重分類數據),對新的結果再次重分類為4類,得到研究區地下水可開采情況建議圖,即值為4、5的可大量開采區,值為6、8的可適宜開采區,值為10、12的可限量開采區,值為15、18的禁止開采區。
對1985-2015年彌河流域平原區的年降水量和蒸發量進行趨勢分析和M-K檢驗法分析,結果如圖2所示。
由圖2(a)可以看出,1985-2015年彌河流域多年平均降水量為561.5 mm,年最大降水量為904.3 mm(1990年),最小降水量為367.5 mm(2002年);由圖2(b)可知,年蒸發量多年平均值為1 003.6 mm,在898.3(1990年)~1 173.6 mm(1986年)間波動。兩者M-K檢驗統計量|Z|<1.96,均未通過0.05置信水平檢驗,變化不明顯,年降水量呈極緩慢上升趨勢,年蒸發量呈緩慢下降趨勢。
4.2.1 年際變化 計算彌河流域平原區56個監測井1985-2015年地下水埋深平均值,繪制年際變化曲線,結果如圖3。
由圖3可以看出,地下水埋深總體呈逐漸增大的趨勢,結合M-K趨勢檢驗,|Z|=7.04,通過了0.05置信水平檢驗,埋深從1985年的12.38 m增加到2015年的22.32 m,增速為0.33 m/a。1985-2000年,彌河流域平原區地下水埋深在12~21 m之間,即地下水埋深從12.38 m增加到20.15 m,平均速率為0.52 m/a。其中1990年平原區蒸發量較小,地下水埋深較上年增加速率明顯緩和,由于1990年降雨量較大且1991年流域蒸發量偏小,水位出現回升。2000-2015年,地下水埋深在19~23 m之間,埋深增加速率逐漸變緩,2000年開始實際開采量逐漸減小,如2003、2007年,降雨量偏大且蒸發量較小,地下水位開始有所上升,特別在遇豐水年,如2004年,地下水位明顯有大幅度的上升,此期間地下水位有升有降逐漸趨于平緩,2010年后地下水埋深維持在22 m左右。

表2 地下水位數據其及影響因素之間的相關系數值與關聯程度對應表

圖2 1985-2015年彌河流域平原區年降水量及蒸發量的動態變化

圖3 1985-2015年彌河流域平原區平均地下水埋深年際變化
4.2.2 年內變化 繪制年徑流頻率曲線,取25%、50%、75%頻率對應的年份為豐水年(1995年)、平水年(1999年)、枯水年(2006年),分析豐、平、枯3年的地下水埋深年內變化情況結果如圖4。
由圖4可知,彌河流域平原區各觀測井地下水位在豐、平、枯水年內的變化趨勢基本一致,2-3月達到最高水位,然后水位逐漸開始下降,在5-6月左右降至年內最低水位,受6-8月較多的降水補給,地下水位逐漸上升。年降水量大小不同導致年內水位回升幅度不同,年降水量越小,地下水位的回升幅度越緩。年降水量越小,年內水位受年內蒸發量、地下水開采等其他因素影響的變化越明顯,在枯水年,如果降雨偏小,地下水位會持續下降。
4.3.1 地下水埋深縱向變化 利用Kriging插值法繪制出1985、1995、2002、2009、2012、2015年6個年份地下水埋深圖,結果如圖5所示。
從圖5的地下水埋深空間分布可看出,1985年埋深自西南山前傾斜區向東北近海區逐漸變淺,自1985年起,埋深總體逐年增大,東部近海區和中部平原區埋深增加速率更快,其中中部平原區逐漸形成一個明顯的漏斗區。
1985年研究區地下水埋深較淺,平均埋深12.38 m,最小埋深0.27 m,最大埋深36.34 m,東北近海區埋深普遍低于10 m,青州市與壽光市交界上小范圍出現10~20 m的埋深區;1995年研究區中部埋深增加明顯,流域平均埋深17.56 m,青州市與壽光市交界處出現降落漏斗跡象,中心地下水埋深為25.43 m;2002年研究區平均埋深為20.74 m,最大埋深43.34 m,中部漏斗區已較明顯,中心埋深達39.43 m;2009年東北近海區埋深基本都增加至15~20 m,中部降落漏斗迅速擴張,漏斗中心埋深為46.08 m;2012年埋深下降速率趨于緩慢,平均埋深21.80 m,最大埋深51.83 m,降落漏斗范圍基本穩定不變,中心埋深增加至51.23 m;2015年平均埋深22.32 m,最小埋深0.63 m,最大埋深53.87 m,漏斗中心埋深51.8 m。
4.3.2 不同埋深段面積變化 通過ArcGIS中的計算幾何對各年不同埋深段所占的面積進行統計,得到6個年份不同地下水埋深的面積變化情況,如表3所示。

圖4 研究區地下水埋深豐、平、枯年內變化

圖5 研究區不同年份的地下水埋深空間分布

表3 不同年份不同地下水埋深段面積統計表 km2
由表3可知,彌河流域平原區在1985-2015年間,地下水埋深不同分段所占面積在1985-2009年間變化較劇烈,自2010年開始變化趨勢基本平穩。小于5 m的埋深段所占面積呈減小趨勢,減少幅度為24.38 km2/a;埋深大于5 m小于15 m范圍所占面積年均下降幅度先減小后增大,其中1985-1995年間年均面積減少幅度較大,約50.4 km2/a,2009年前后面積減少幅度較小,約18.86 km2/a,2012-2015年間年均面積減少幅度最大,約53.2 km2/a;埋深大于15 m小于25 m范圍所占面積呈先減小后增大的趨勢增加,其中1995年前年均面積增幅較大,約61.07 km2/a,2012-2015年間面積增加65.07 km2/a,增幅最大,2009-2012年間增加0.58 km2/a,增幅最小;埋深大于25 m范圍所占面積增幅較小,2009年以前逐漸增大,2009年以后基本保持穩定不變,面積增幅最大為23.23 km2/a,最小為8.73 km2/a。
地下水動態變化受氣象因素和人為活動共同影響。伴隨國民經濟快速發展和水資源供需不平衡加劇,開始開發利用地下水,由《山東省水資源公報》可知,濰坊市開采量自20世紀70年代后不斷增加,從1985年10.79×108m3增至2000年15.73×108m3,地下水位不斷下降;社會經濟發展也間接地影響地下水位的變化,《山東省水利統計年報》 中表明,研究區近15年來,GDP由218×108元增至1 538.9×108元,增長了7倍。
選擇壽光136#、184#、121#和青州115#、100#共5個典型觀測井(見圖1),利用SPSS軟件,對地下水位數據與研究區的降水量、蒸發量、開采量、GDP進行雙變量相關性分析,結果見表4。
由表4可知,觀測井136#、184#地下水位變化與降水量相關程度很弱,其余3個均為強相關,說明該研究區地下水動態變化明顯受到降水量的影響。5個觀測井地下水位變化與蒸發量和開采量的皮爾遜相關系數0.4<|γ|<0.8,均達到了中等相關或強相關,且為負相關,說明研究區蒸發量和人工開采量對地下水動態變化影響很大,即隨著蒸發量和開采量的增大,地下水位越低。研究區地下水位與GDP存在負相關,隨著社會經濟的飛速發展,水資源供不應求,過度地開采使用地下水,導致地下水位持續下降,出現局部地下水漏斗區。
利用ArcGIS,以1985-2015年為研究期,2015年為現狀水平年,通過水位動態法分析研究期內地下水埋深的分布及變化速率情況,繪制出地下水埋深、研究區內地下水埋深變化及地下水可開采情況建議圖,具體分析結果見圖6~8,研究區不同分區的面積見表5。

表4 雙變量相關性分析結果
注:**為通過1%的顯著性水平檢驗;* 為通過5%的顯著性水平檢驗。

圖6彌河流域平原區現狀水平年(2015年)地下水埋深分布圖 圖7彌河流域平原區研究期內地下水埋深變化圖 圖8彌河流域平原區地下水可開采情況建議圖

表5 研究區不同分區的面積
從圖6可明顯的看出,彌河流域平原區地下水埋深幾乎均超過6 m,淺埋區僅占區域面積的1.8%,埋深大于25 m的區域主要集中在壽光與青州交界處以及山前傾斜區,由圖7和表5可知,彌河流域平原區的超采區面積是非超采區面積的3倍,分別為2549.7和829.5 km2,一般超采區范圍較大,主要分布于研究區西北至東南走向。由圖8和表5可知,可適宜開采區主要位于研究區北部和南部,其中壽光市、廣饒縣、青州市、昌樂縣、臨朐縣均有分布,總面積為710.2 km2,占研究區面積21.0%,可限量開采區幾乎遍布研究區,面積高達2 058.3 km2,而在青州市與壽光市交界線處的降落漏斗處應禁止開采地下水。
本文以彌河流域平原區56眼監測井五日地下水數據(1985-2015年)為基礎,結合當地氣象、地下水開采、社會經濟發展情況,利用ArcGIS和SPSS統計方法,對研究區地下水埋深時空變化特征和可開采情況進行分析,得出以下結論:
(1)彌河流域平原區1985-2015年的地下水埋深逐漸增加,2000年后增幅平緩,年內地下水位一般在2-3月出現最高值,在5-6月降至最低值;
(2)彌河流域平原區1985、1995、2002、2009、2012、2015年地下水埋深空間分布相似,從西南到東北,地下水埋深隨地勢大體降低,也存在一定空間變異性,埋深小于15 m的范圍逐漸縮小,埋深大于15 m的區域面積不斷增加,其中1995年左右青州市與壽光市交界處開始出現小范圍降落漏斗,漏斗區不斷擴大,至2012年漏斗范圍基本不再增大,但漏斗中心埋深仍不斷增大;
(3)彌河流域平原區可適宜開采區僅占總面積的21.0%,主要位于壽光、青州、臨朐等市縣,研究區大多為可限量開采區(面積高達2 058.3 km2),山前一帶傾斜區和降落漏斗區屬于禁止開采區,應當采取措施合理保護利用地下水;
(4)5個典型監測井的水位數據大部分與各影響因子皮爾遜相關系數|γ|>0.4,通過1%顯著性水平檢驗,相關性高,表明氣候變化和人類活動對地下水的動態變化均有影響。
參考文獻::
[1] 王 瑗,盛連喜,李 科, 等.中國水資源現狀分析與可持續發展對策研究[J].水資源與水工程學報,2008,19(3):10-14.
[2] 胡汝驥. 中國天山自然地理[M]. 北京:中國環境科學出版社,2004.
[3] ALLEY W M,TAYLOR C J . The value of long-term ground water level monitoring[J]. Ground Water,2010,39(6):801-801.
[4] MACHIWAL D, MISHRA A, JHA M K, et al. Modeling short-term spatial and temporal variability of groundwater level using geostatistics and GIS[J]. Natural Resources Research, 2012, 21(1): 117-136.
[5] MA T S, SOPHOCLEOUS M, YU Y S. Geostatistical applications in ground-water modeling in south-central Kansas[J]. Journal of Hydrologic Engineering, 1999, 4(1): 57-64.
[6] SUN Yue,KANG Shaozhong,LI Fusheng et al. Comparison of interpolation methods for depth to groundwater and its temporal and spatial variations in the Minqin oasis of northwest China[J]. Environmental Modelling and Software,2009,24(10):1163-1170.
[7] 李文體,劉向華,馮謙誠.河北省地下水超采區劃分及超采狀況分析研究[J]. 水資源保護,2001,17(4):26-30+72.
[8] 黃曉燕,馮志祥,李 朗,等.江蘇省地下水超采區劃分方法對比研究[J].水文地質工程地質,2014,41(6):26-31.
[9] 陳國浩.濟寧市地下水超采區劃分研究[D]. 南京:河海大學,2007.
[10] 高樹東. 濰坊市地下水資源評價[D]. 南京:河海大學,2005.
[11] 曹 濱. 濰坊北部地下水漏斗成因機理與演變規律研究[D].濟南:山東農業大學,2017.
[12] 徐月琴,孫景林,王沁宇.濰坊市濰北區地下水降落漏斗成因分析[J].山東水利,2018(6):15-16.
[13] 隋 偉,苗乃華,陳吉賢,等.濰坊市地下水超采現狀及對策[J].地下水,2010,32(3):51-52.
[14] 曹潔萍,遲道才,武立強,等. Mann-Kendall檢驗方法在降水趨勢分析中的應用研究[J].農業科技與裝備,2008(5):35-37+40.
[15] 汪冬華. 多元統計分析與SPSS應用[M]. 上海:華東理工大學出版社,2010.
[16] 孫逸敏. 利用SPSS軟件分析變量間的相關性[J]. 新疆教育學院學報,2007,23(2):120-123.
[17] 黃耀裔. 基于普通Kriging的地下水空間插值研究——以福建晉江市為例[J]. 廊坊師范學院學報(自然科學版), 2013,13(6):12-17.
[18] 曾懷恩,黃聲享. 基于Kriging方法的空間數據插值研究[J]. 測繪工程, 2007,16(5):5-8+13.
[19] 郭秀紅,趙 輝.淺談地下水超采區劃分[J].中國水利,2015(1):41-43.
[20] 吳秀芹. ArcGIS 9 地理信息系統應用與實踐[M].北京:清華大學出版社,2007.