弥河流域平原区地下水埋深时空变化特征及可开采性分析

杨滢嘉， 张丹蓉， 庄会波， 刘 江， 管仪庆， 邵广文

(1.河海大学 水文水资源学院，江苏 南京 210098；2.山东省水文局，山东 济南250000)

1 研究背景

水资源是社会经济发展和人民生产生活的主要影响因素，地下水作为水资源的重要组分，具有水质良好、水量稳定、调节力强、便于开采等优点，已成为城镇发展建设、人民饮水安全、生态环境改善的重要保障[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年地下水观测井五日观测资料，分别对地下水埋深年际年内变化趋势、埋深变幅速率、不同埋深段面积变化特征进行分析，浅析其影响因素，并划分流域超采区，评价地下水开采适宜性。研究成果有助于了解弥河流域平原区地下水埋深变化情况，对有效控制局部降落漏斗及合理开采利用地下水资源有重要的参考价值。

2 研究区概况

弥河是山东省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))，主要补给来源有降水补给、回灌补给、侧向径流补给等。

3 数据来源和研究方法

3.1 数据来源

收集整理所需的基本资料，本文所采用的数据信息如表1所示。

图1 弥河流域地理位置、高程分布及水文地质图

表1 研究所需基本资料数据信息及来源

3.2 M-K检验法

Mann-Kendall是进行时间序列趋势分析的一种常用的非参数检验法，以使用约束少、人为干扰小、适用范围宽等优点被广泛使用，其计算结果可信度高[14]。在趋势检验中，统计变量Z>0，表示序列呈上升趋势；Z<0，表示呈下降趋势。其中|Z|>1.96，说明上升或下降趋势显著。

3.3 相关性分析法

相关性分析法是对两个或多个存在关联的变量因子进行分析计算，通过相关系数来衡量变量因子之间的密切度[15-16]。

本文使用该方法计算地下水位数据及其影响因子之间的相关系数，分析影响地下水动态变化的因素，相关系数对应相关程度如表2。

3.4 Kriging插值法

Kriging法是ArcGIS地统计学中的一种对空间未知样点的区域化变量的取值进行线性无偏最优内插的方法，该方法最大程度地利用了已知样点包括的所有信息[17-18]。其基本思路是首先通过分析区域化变量在空间位置上的分布，选择最佳理论变异函数模型，其次确定对未知样点进行插值的影响范围，在范围内根据已知样点与未知样点的空间位置和相关性的不同，赋予每个已知样点相应的权重，进行加权平均后估计出待测样点属性值。本文对数据的空间趋势走向和残差项综合分析后，采用Geostatistical Analyst工具中的普通Kriging插值法对研究区56眼监测井地下水位值进行空间插值，转为90 m分辨率的地下水位栅格图，将90 m分辨率的DEM与同分辨率的水位栅格图相减，得到1985、1995、2002、2009、2012、2015年地下水埋深分布图。

3.5 超采区划分方法

水位动态法是指以地下水水位年均下降速率和现状年的地下水埋深作为判断指标进行超采区划分的方法[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的禁止开采区。

4 结果分析与讨论

4.1 气候因子变化

对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 地下水埋深动态变化特征

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 地下水埋深时空变化特征分析

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。

4.4 地下水动态变化特征成因分析

地下水动态变化受气象因素和人为活动共同影响。伴随国民经济快速发展和水资源供需不平衡加剧，开始开发利用地下水，由《山东省水资源公报》可知，潍坊市开采量自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存在负相关，随着社会经济的飞速发展，水资源供不应求，过度地开采使用地下水，导致地下水位持续下降，出现局部地下水漏斗区。

4.5 地下水可开采性探讨

利用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，而在青州市与寿光市交界线处的降落漏斗处应禁止开采地下水。

5 结 论

本文以弥河流域平原区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.