三江平原典型区地下水流场变化及主要影响因素分析

束龙仓，徐丽丽，袁亚杰，吕 岩，鲁程鹏，刘 波

（1.河海大学 水文水资源学院，江苏 南京 210098；2.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；3.黑龙江省水文水资源中心，黑龙江 哈尔滨 150001；4.黑龙江省佳木斯水文水资源中心，黑龙江 佳木斯 154002）

1 研究背景

三江平原土壤肥沃，是我国九大商品粮生产基地之一。其所属的三江盆地属于一个大型地下水汇水盆地，地下水在重力势能驱动下，由山前地带高势能区向黑龙江河谷的低势能区运动，进而形成三江平原区域地下水流动系统。1980年代前，三江平原地下水开发利用程度较低，地下水流场处于天然状态，地下水位的年际变化与降水及地表水的年际变化相似［1］。但随着社会经济的快速发展及水田范围的不断扩大，地下水开采尤其是农田灌溉用水大幅增加，人工开采的影响日趋明显［2］。

三江平原地下水开发利用引发的问题一直备受关注。郭龙珠［3］应用偏最小二乘法研究了三江平原宝清县地下水位动态变化规律，并建立地下水位动态模型分析了地下水位变化的转折点；Yang等［4］采用自适应神经模糊推导系统对建三江地区胜利农场的水位进行了分析和预测；刘东等［5-6］采用分形理论与小波理论有机结合方法、近似熵方法对建三江分局地下水埋深序列复杂性进行测度分析；危润初等［7］采用双向回归突变分析方法研究了黑龙江建三江地区（1992—2011）地下水埋深序列的趋势突变点；刘伟朋等［8］分析了三江平原地下水流场时空演化特征，并提到不同因素对地下水位会产生的影响；束龙仓等［9］对三江平原典型区地下水与河流代表性河段河水之间水量交换的时空变化规律进行了深入分析。虽然目前已有一些学者对三江平原地下水位、流场变化问题进行了相关研究，但对该区地下水位下降幅度及是否处于超采状态结论不一［10］。一些学者研究显示三江平原由于灌溉用水增多，部分区域地下水位大幅下降，出现地下水超采问题［8］，甚至有研究表明有些区域出现严重的地下水降落漏斗，对农业生产和居民生活用水构成了威胁［11-12］。另一些研究则表明整个三江平原的潜水区可判断为无超采，还有一定的开采潜力［13-14］，甚至发现2005—2010年间三江平原地下水位没有明显的下降，仅是随着降雨的丰平枯年份动态变化［15］。不同学者对于引起研究区地下水位变化的主要原因也没有一致结论，王韶华等［16］通过分析得出1997—2002年三江平原地区地下水位下降是由大规模种植水稻、大量开采地下水所引起；赵孟芹等［17］对一些代表性区域地下水位动态进行了分析，认为气候因素对地下水位变化起到了决定性的作用；刘东等［5-6］则研究认为地下水开采是影响地下水位动态变化的关键因素。

综上，受时间、资料等限制，目前对三江平原地下水流场变化情况和主要影响因素还未能得出统一定论。此外，已有对三江平原地下水位变化原因的研究仅仅是停留在定性分析层面，不同影响因素对地下水位的影响程度有待进一步研究揭示。交叉小波变换能够有效诊断不同信号间的相关性、时延性和相位结构，适宜于分析地下水位和各影响因素间的时延相关特征和时频相位关系。在交叉小波分析得出时频相位关系的基础上，辅以皮尔逊相关分析，能够较好地研究各影响因素对地下水位的影响。灰色关联度分析可以用来分析两个因素之间的关联程度，可以定量比较各影响因素对地下水位影响程度的相对大小。为此，本文选取三江平原典型区域作为研究区，进一步探明2001—2019年地下水流场变化情况，并分析地下水位变化的空间差异。基于地下水位变化空间差异选取代表性站点，采用交叉小波分析、皮尔逊相关分析及灰色关联度分析等方法，研究了降水及地下水开采这两大主要影响因素对不同区域地下水位变化的影响程度，对于深入探究区域地下水位变化和指导该区地下水资源合理开发利用具有重要意义。

2 研究区概况及数据来源

2.1 研究区概况三江平原由黑龙江、乌苏里江及松花江冲、洪积作用所形成，面积约10.89万km2。研究区位于三江平原腹地，是三江平原地下水资源开发利用问题最为突出的地区。同时兼顾以河流划分的自然边界和行政区边界，确定了包括富锦市、同江市、抚远市、饶河县部分区域作为本文的研究区，北部、东部、南部、西北部分别以黑龙江、乌苏里江、挠力河、松花江为界，总面积约为2.21万km2（图1）。该区域为温带湿润、半湿润大陆季风气候，多年平均气温为2.8℃，多年平均降水量为383.5～886.1 mm，其中70%的降水多发生在6—9月，多年平均蒸发量约为800 mm。研究区地势西南高于东北，除少数山丘外，绝大部分是平原沼泽地带，地面坡度为1/10 000～1/5000。

图1 研究区位置

研究区含水层结构单一，区内堆积有厚达200余米的第四系砂及砂砾等松散沉积物，形成巨厚的孔隙含水岩组，孔隙水在上覆连续较厚黏土层分布的地区为微承压水及层间水，其他地区为潜水。松散岩类孔隙水具有厚而稳定的含水层，地下水蕴藏丰富，是研究区地下水开采利用的主要对象。根据单井涌水量Q，区内松散岩类孔隙水分布可划分为水量极丰富（Q＞5000 m3/d），水量丰富（1000 m3/d＜Q＜5000 m3/d）、水量中等（100 m3/d＜Q＜1000 m3/d）、水量贫乏（Q＜100 m3/d）4类区域。此外，研究区还存在基岩裂隙水，其富水性较差，在区内分布范围较小。

大气降水入渗补给是研究区松散岩类孔隙水最为主要的补给方式，灌溉入渗量和侧向径流补给也是地下水补给量的重要组成部分。地下水排泄途径以人工开采为主，同时还存在侧向径流排泄和少量蒸发排泄。此外，研究区内的沿河地区在汛期可接受来自河流的地表水补给，而枯水期则是地下水向河流排泄。

2.2 数据来源本文气象、水文、地下水位动态观测资料和水资源开发利用资料收集于中国气象科学数据共享服务网、黑龙江流域水文年鉴、黑龙江省水文局和项目监测、统测及实验数据。其中地下水位监测井共计100眼，所采用的地下水监测井均为潜水或弱承压水监测井，主要为2000年7月—2019年3月序列，部分站点水位资料缺失不连续；降水资料主要采用研究区10个雨量站2000年—2018年降水量数据。地下水开采数据为各级水资源分区结合行政区的地下水开采量，为2010—2018年序列，对此根据灌区分布等实际情况进行整理得到研究区地下水开采量及各灌区地下水开采量等数据。地下水位监测站和雨量站位置分布见图1。

3 研究方法

本文采用空间插值及栅格代数方法研究三江平原典型区2001—2019年地下水流场变化情况，采用交叉小波分析、皮尔逊相关性分析方法，分析降水及地下水开采变化对不同区域地下水位变化程度的影响，并采用灰色关联度分析方法，定量比较降水及地下水开采对不同区域地下水位影响程度的相对大小。

3.1 交叉小波分析交叉小波变换（CrossWavelet Transform，XWT）是一种信号分析技术，可以对两个时间序列在不同时频域中的相互关系进行分析研究［18］。本研究中采用交叉小波变换分析地下水位变化与降水及地下水开采之间的关系，两个时间序列xn和yn的交叉小波变换可以定义为WX、WY。其交叉小波谱为WXY＝WXWY*，其中WY*代表WY的复共轭。它们的交叉小波功率谱密度被定义为，其值越大，二者相关程度越高。复数辐角arg（WXY）可以看作为时频空间中xn和yn之间的局部相对相位。对连续交叉小波功率谱的检验方法为：假定时间序列xn和yn的期望谱为红色噪音谱pXk和pYk，两个时间序列的交叉小波功率谱理论分布可参考文献［19］，其表示如下：

式中：σX、σY分别为时间序列xn、yn的标准差；v为自由度；Zv（p）为置信度。

3.2 皮尔逊相关分析皮尔逊相关系数（Pearson Correlation Coefficient），可以反映两个随机变量之间的线性相关程度。给定两个随机变量X、Y，其皮尔逊相关系数r的计算公式如下：

式中：r为皮尔逊相关系数，无量纲；n为样本数量；Xi、Yi为变量X、Y对应的i点观测值；分别为X、Y样本平均数。

r的取值在-1与1之间，其绝对值越大表示相关程度越高。r值大于0表示正相关，取值小于0表示负相关，取值等于0表示线性不相关。

3.3 灰色关联度分析灰色理论系统是指部分信息已知，部分信息未知的系统［20］。灰色关联度分析（Grey Relation Analysis，GRA），是一种以灰色系统为基础，根据因素之间发展趋势的相似或相异程度来衡量因素间关联程度的一种方法。

本研究以研究区地下水水位变化为母序列，记为x0（t），降水序列x1（t）和地下水开采量序列x2（t）为子数序列，则在t＝k时刻地下水水位变化与其相关影响因素之间的关联系数γ（x0（k），xi（k））为：

式中：γ（x0（k），xi（k））为关联系数；δ为分辨系数，研究中取0.5；i为子序列下标；k为某时刻；分别为两级最小差和两级最大差。

母序列地下水水位变化与子序列降水、地下水开采量等的关联度ri的计算公式为：

式中各参数意义同式（3）。

4 结果与分析

4.1 2001—2019年三江平原典型区地下水流场变化分析根据相关研究［21］，研究区潜水完整的水文周期划分自当年3、4月至翌年3、4月较合适，且农田灌溉用水主要集中在4—8月份，故取每年3月地下水位作为每年地下水位变化依据，以进行地下水流场年际变化分析。由于地下水位监测站分布较为均匀，利用较为常用的反距离权重法对2001年3月和2019年3月地下水位数据进行插值分析，分别绘制流场分布如图2所示。

图2 研究区2001年3月和2019年3月地下水流场对比

由图可见2001—2019年间地下水整体流向未发生变化，为自西南流向东北，但区域内地下水位明显下降。2001年3月研究区地下水位平均值为47.51 m，水位最大值为60.30 m，最小值为37.56 m；而到2019年3月地下水位平均值为46.16 m，水位最大值为59.30 m，最小值为36.22 m，均有所减小。研究期间地下水流场形态发生了较为明显的变化，2001年3月西南部富锦市内水位较高，地下水水力梯度比较大，地下水位整体自西南向东北方向逐渐减小，最终在东北方向排泄入河流或在河流下方以地下水侧向排泄的方式流出区外。而到2019年3月富锦市附近区域水力梯度变小，地下水径流变缓，中东部地区地下水位整体下降较为明显，在中部甚至形成了一定范围的降落漏斗，改变了流网的形态，说明研究区中部地下水开采较多，引起地下水流动路径发生变化。从水田面积变化情况来看（图3），1980—2018年水田面积呈上升趋势，1980—2000年间呈小幅度增长，而2000年以后呈现大幅度上升趋势，由2000年的742.56 km2上升为2018年的10 507.87 km2，增长了13.15倍，成为研究区内的第一大用地类型，占总土地面积的47.68%［22］。而水田的发展依赖于灌溉，区内灌溉引用部分地表水，但整体以抽取地下水灌溉为主，地下水灌溉用水量占农业灌溉总用水量的70%以上（图4），农业地下水灌溉开采量的增加导致地下水位的整体下降。因此，研究期内地下水位整体下降与水田的迅速扩张有关。

图3 研究区1980—2018年水田面积变化

图4 研究区2010—2018年农业灌溉用水情况

采用栅格代数方法，将2001年3月和2019年3月地下水位进行代数运算，得到地下水位变幅如图5所示。根据栅格统计，研究区2001—2019年期间，地下水位下降区域面积为1.86万km2，占整个研究区的83.97%，即大部分地区地下水位呈现下降趋势，仅富锦市附近等小部分范围内有水位抬升现象。地下水位变幅平均值为-1.37 m，区内平均地下水位年均降幅为0.07 m，中部和东南部地区水位下降较为严重，水位降幅最大可达8.29 m。

研究区地下水位下降区域主要分布在中东部。从水田和灌区分布来看（图6），研究区水田主要分布在中东部区域，灌区依随水田也多分布在中东部地区，区内分布有青龙山灌区、八五九灌区、乌苏镇灌区、勤得利灌区、三村灌区以及临江灌区。从农场分布来看，研究区农场同样多集中分布在中东部区域（图5），并且农场集中分布范围与地下水位降幅1 m以上区域大体吻合。由此认为研究区中东部地区地下水位的明显下降，主要是由于该区域是水田、农场集中分布区，因农业灌溉需要开采大量地下水，导致该区域地下水位出现大幅下降。研究区西部存在地下水位小幅下降区域，区内同样有零星水田分布（图6），水田发展造成该区域地下水开采量增加，由于该区域地势较高，地下水沿东北方向排泄，地下水经开采后补给恢复较慢，故地下水位出现轻微下降现象。对于富锦市附近的地下水位上升区域，认为是由气象因素造成，2000—2018年间研究区降水量增加了49.60～259.19 mm，且在富锦市与同江市附近降水量增幅最大［22］。

图5 2001年3月—2019年3月地下水位变幅图

图6 研究区土地利用类型及灌区分布

根据2001—2019年地下水位变化情况，在地下水位变化程度不同区域选取4个代表性地下水位监测站，各站点所在位置如图5所示。

根据各站点2001—2019年地下水位过程线（图7），可见4个代表站点仅菜咀子站地下水位有小幅上升，上升2.74 m，其他3个站点水位有不同程度下降，前锋农场站下降6.06 m，别拉洪站下降1.06 m，临江站下降0.39 m。4个地下水位监测站能够展现研究区内2001—2019年不同水位变化情况，具有代表性，以下研究将基于此4个代表站进行。

图7 代表站2001—2019年地下水位变化过程线

4.2 降水变化对区域地下水位影响分析

4.2.1 交叉小波分析 为了探究降水变化对不同区域地下水位变化的影响程度，采用4个代表性地下水位监测站2000年7月—2018年12月逐月地下水位和附近雨量站所测得的逐月降水量进行交叉小波分析（图8）。图中，实线内区域通过95%置信度的红噪声标准谱检验，细弧线区域表示小波影响椎内的有效谱值。

图8 代表站降水和地下水位交叉小波功率谱

由菜咀子站降水与地下水位变化的交叉小波谱图8（a）可知，在11～13月周期内，菜咀子站降水和地下水位序列在整个2000—2018年期间表现出通过95%置信度检验的显著共振关系。箭头基本指向右下方，表示地下水位与降水同相位，二者为正相关关系，且地下水位滞后于降水。这些特征表明，在整个研究期间，降水对年尺度上的地下水位变化有很强的影响。根据平均相位角，地下水位相对于降水变化的滞后时间约为11.52 d。此结果说明菜咀子站所在的地下水位上升区域地下水位受气候因素影响显著。

图8（b）和图8（c）分别为临江站和别拉洪站降水与地下水位变化的交叉小波谱，由图可见箭头基本指向右下方，说明临江站和别拉洪站降水和地下水位同相位。二者在整体上同样存在11～13月尺度上的显著共振关系，但相关性特征存在不稳定现象，如临江站2005年1月—2008年7月、别拉洪站2000年7月—2006年5月在11～13月尺度上共振关系不显著，分析认为可能是人为开采等因素在这些时间段内对地下水位的影响较大，从而使地下水位受降水影响程度减小。

图8（d）为前锋农场站降水与地下水位变化的交叉小波谱，可见此站在11～12月周期尺度上虽然通过了红噪声检验，但降水与地下水位交叉小波功率谱能量较低，未表现出较强的关联性，并且箭头指向左上方，降水与地下水位变化呈反相位，二者为负相关关系。该结果说明前锋农场所在的农场中部地下水位大幅下降区域，研究期间地下水位变化受降水影响程度非常微小。

4.2.2 皮尔逊相关分析 对4个代表性地下水位监测站地下水位变幅与其附近雨量站降水数据分别进行皮尔逊相关性分析显示（表1），菜咀子站地下水位与降水存在相关性，其相关系数为0.485，通过0.01级别的双尾检验；临江站和别拉洪站地下水位变幅与降水呈现不同程度的微弱相关性；前锋农场站地下水位变幅与降水相关性较差。该结果进一步说明：地下水位上升区域地下水位受降水影响较大，在地下水位降幅较小区域，地下水位在一定程度上受降水变化影响，而地下水位下降较严重区域降水对地下水位的影响程度较小。

表1 代表站水位变化与降水相关性分析结果

4.3 地下水开采对区域地下水位影响分析

4.3.1 交叉小波分析 为了探究地下水开采对不同区域地下水位变化的影响程度，采用4个代表性地下水位监测站2010—2018年地下水位月降幅和其所在区域地下水月开采量进行交叉小波分析（图9），其中临江站、别拉洪站以及前锋农场站位于灌区内，地下水开采量数据采用其所在灌区的地下水开采量，菜咀子站位于灌区外，站点及附近地区无地下水灌溉开采，其地下水开采量序列采用其所在行政区生活用水等其他地下水开采量。

由菜咀子站地下水开采与地下水位变化之间的交叉小波谱图9（a）可知，二者交叉小波功率谱能量较低，未表现出较强的关联性，也不存在稳定的周期性相关。该结果说明菜咀子站所在的地下水位上升区域地下水受开采影响程度较小，原因是由于菜咀子站位于灌区外，该区域地下水开采主要为生活用水开采，无灌溉开采，地下水开采量较少。

图9（b）和图9（c）分别为临江站和别拉洪站地下水开采与地下水位变化的交叉小波谱，由图可见两站地下水开采与地下水位变化整体上主要存在11～13月尺度的显著共振关系，但是该尺度上的相关性特征不稳定。临江站地下水开采与地下水位变化在2010—2013年间无明显共振关系，2014年之后两者共振关系显著，该结果说明在2010—2013年间，该区域地下水位受地下水开采影响较小，2014年后影响逐渐增大。此变化与其所在区域水田发展有关，由于该区域水田范围扩大，地下水开采量渐增，地下水开采对地下水位的影响逐渐增大，2014年之后开始对地下水位有显著影响。

由前锋农场站地下水开采与地下水位变化的交叉小波功率谱图9（d）可知，在11～13月周期内，地下水开采和地下水位变化在整个2010—2018年期间表现出通过95%置信度检验的显著共振关系，箭头基本指向右下方，表示地下水开采与地下水降幅同相位，二者为正相关关系，且地下水位变化滞后于地下水开采。根据平均相位角计算得，地下水位相对于地下水开采变化的滞后时间约为10.8 d。这些特征表明，前锋农场站所在的中部地下水位大幅下降区域，在整个研究期间，地下水开采对年尺度上的地下水位变化有很强的影响，地下水位受地下水开采影响显著。

图9 代表站地下水开采和地下水位变化交叉小波功率谱

4.3.2 皮尔逊相关分析 对4个代表性地下水位监测站地下水位变幅与其所在地区的地下水开采量数据分别进行皮尔逊相关性分析显示（表2），仅前锋农场站地下水位变化与地下水开采间存在一定相关性，其相关系数为-0.409，通过0.01级别的双尾检验，其他站点地下水位变化与地下水开采之间的相关性较差。此结果进一步说明在地下水位大幅下降区域，地下水位受地下水开采的影响较大，而在地下水位降幅较小及水位上升区域，地下水位受地下水开采影响相对较小。

表2 代表站水位变化与地下水开采相关性分析结果

4.4 降水和地下水开采对区域地下水位相对影响程度分析地下水位变化受到多种因素影响，而降水变化和地下水开采是影响研究区地下水位变化的两个主要因素。菜咀子站由于处于灌区外，受地下水开采影响较小，上文此站降水与地下水位的交叉小波结果显示此区域地下水位变化受降水影响较大，故该区域降水对地下水位的影响程度明显大于地下水开采对地下水位的影响程度。为了进一步分析灌区内地下水开采与降水变化对于不同区域地下水位变化影响程度的相对大小，对前锋农场站、临江站、别拉洪站3个代表性站点分别进行灰色关联度分析，得到各站点地下水位与降水及地下水开采量之间的关联度。

由表3可见，位于研究区地下水位下降区域的3个代表站地下水开采与地下水位变化间的关联度均大于降水与地下水位变化间的关联度，表明其所在区域地下水开采对地下水位的影响程度均大于降水对地下水位变化的影响程度。前锋农场站位于灌区中部地下水位大幅下降区域，受地下水开采影响较大，其地下水开采与地下水位变化间关联度明显大于降水与地下水位变化间的关联度。而别拉洪站和临江站地下水位和降水及地下水开采之间的关联度差别不大，表明此两站所在的地下水位小幅下降区域虽然受地下水开采的影响程度更大，但同时也很受降水情况的影响。

表3 灰色关联度分析结果

5 结论

本文选取三江平原典型区域作为研究区，采用空间插值方法分析了区内2001—2019年地下水流场变化情况。继而利用栅格代数方法计算了区内地下水位变化的空间差异，基于地下水位变化空间差异选取代表性站点，采用交叉小波变换分析了地下水位和降水及地下水开采间的时延相关特征和时频相位关系，并在小波分析得出时频相位关系的基础上，辅以皮尔逊相关分析，使结果更具可靠性。最后，采用灰色关联度分析方法，比较了降水和地下水开采这两大主要影响因素对不同区域地下水位影响程度的相对大小。相比于已有研究，本研究采用了更充足的研究数据和更系统的研究方法，定量地阐明了降水和地下水开采对于三江平原地下水位及其流场的影响。通过研究得出以下具体结论：

（1）研究区2001—2019年地下水大体流向未发生明显变化，但地下水位呈现整体下降趋势，下降区域面积为1.86万km2，占整个研究区的83.97%，地下水位变幅平均值为-1.37 m。中部和东南部地区水位下降较为严重，在中部甚至形成了一定范围的降落漏斗，水位降幅最大可达8.29 m；富锦市附近小部分范围内有水位抬升现象。

（2）研究区地下水位的整体下降与研究期间水田面积迅速扩张有关。研究区2000—2018年间水田面积增长了13.15倍。中部地下水位大幅下降区域主要为水田、农场集中分布区，水位降幅1 m以上区域与农场集中分布范围大体吻合；西部地下水位小幅下降区域内同样有零星的水田分布。

（3）研究区地下水位上升区域内，地下水位受降水因素影响显著，受地下水开采变化影响较小。水位上升区地下水位与降水在年周期上存在显著共振关系，二者皮尔逊相关系数为0.485，而地下水位与地下水开采之间不存在稳定的周期性相关，相关性较差。

（4）研究区中东部地下水位大幅下降区域内，地下水位受人类开采的影响较大，受降水影响较小。中东部水位大幅下降区地下水位与降水间不存在稳定的周期性相关，相关性较差，而地下水位与地下水开采之间在年周期上存在显著共振关系，二者皮尔逊相关系数为-0.409。

（5）研究区地下水位下降区域中，中部水位大幅下降区域地下水位受开采影响明显大于受降水影响，灌区内水位小幅下降区域地下水位与二者之间关联度相差不大，受地下水开采影响程度略大于受降水影响。