基于地统计分析的河套灌区地下水埋深与矿化度时空变异规律研究

2020-09-05孙贯芳杨金忠

灌溉排水学报 2020年8期

姚玲，杨洋，孙贯芳，朱焱，杨金忠

（武汉大学水资源与水电工程科学国家重点实验室，武汉 430072）

0 引言

灌区土壤水分和盐分是制约作物生长与灌区长久发展的2 个重要因素，土壤盐分不仅受到灌溉制度、种植结构、降雨蒸发等因素的影响，还与地下水的埋深以及地下水矿化度密切相关[1]。河套灌区作为一首制引黄灌区，随着近年来引黄灌溉水量的减少以及灌区节水改造配套工程的实施，灌区原有的水文循环规律随之改变，农田水盐的动态平衡也发生变化[2-4]。地下水是灌区水盐大环境的重要组成成分，地下水在灌溉水和降雨的入渗补给下，其埋深和矿化度均会发生变化，地下水流动过程中与含水层介质会发生一系列的地球化学反应，影响着其化学组成和变化[5]，而土壤盐渍化与地下水文密切相关，因此，地下水埋深与地下水矿化度研究是灌区水盐大环境研究中不可或缺的一部分，其对灌区土壤盐渍化的改善，水资源的合理开发利用和管理及其有关的生态环境保护和建设都具有一定的现实意义[6-7]。

目前，针对河套灌区的地下水研究多以短时间或田间尺度地下水变化特征为研究重点[8]，而利用大空间尺度、长时间序列的地下水埋深和矿化度数据分析其时空变化规律以及相互关系的研究较少。因此，本文借助前人的地统计学方法，以土壤盐渍化情况严重的河套灌区为研究区，从长时间序列分析灌区地下水埋深和矿化度的时空变化特征，为灌区地下水埋深和地下水矿化度之间的关系提供理论依据。

1 材料与方法

1.1 研究区概况

河套灌区是我国三大特大型灌区之一，位于内蒙古自治区的西部，地理坐标为东经106°20′—109°19′，北纬40°15′—41°18′，东西长250 km，南北宽50 km，总土地面积1.94×109m2。灌区自西向东由乌兰布和灌域、解放闸灌域、永济灌域、义长灌域、乌拉特灌域组成[9]。灌区具体地理位置如图1 所示。灌区地处干旱、半干旱、半荒漠草原地带，年平均降水量为130～215 mm，自西向东递增。各地区蒸发能力极强，蒸发量与降雨量之比可达11～16 倍。河套灌区含水层水平方向上的分布规律为：由于坳陷深度自东向西、由南向北递增，含水层厚度沿此方向增厚。地下水以潜水为主，潜水含水层由细沙和中细沙组成。少量半承压水存在于地势低洼，黏质土覆盖层厚度较大的地区。灌区潜水主要来源于田间灌溉入渗及渠道渗漏，其次是山洪水和降水。

图1 河套灌区地理位置 Fig.1 Location map of Hetao Irrigation District

1.2 数据采集

河套灌区5 个灌域共有215 个地下水埋深观测点，取样时间为1998—2017 年，地下水埋深数据每5 d采集1 次，分别是在每个月的1、6、11、16、21、26日，1 年共采集72 批水样。灌区5 个灌域共有117个地下水矿化度观测点，取样时间为1998—2017 年，每年采集7 批次水样，分别是在1、3、5、7、9、10、11 月。

1.3 研究方法

本文在研究地下水埋深，矿化度的空间变异性时主要采用了地统计学的方法，半变异函数是地统计学的基本工具，是一个关于数据点的半变异值与数据点间距离的函数，公式为：

式中：h 为滞后距，或称步长、距离段；E 表示数学期望，Z(x)为在位置x 处的变量值，Z(x,h)为在与位置x 偏离h 处的变量值。随着距离段的变化，可计算出一系列的变异函数值。以h 为横坐标，V(h)为纵坐标作图，便得到了变异函数图。从计算公式可见，变异函数实际上是一个协方差函数，是同一个变量在一定相隔距离上差值平方的期望值，差值越小，说明在此距离段上该变量值的相关性越好；反之亦然。

式（1）是较为严格的数字定义，同时适用于空间上连续分布的变量。但在实际工作中，采样点常常是离散的，为此将式（1）改写为；

式中:N(h)为距离等于h 的点对数;Z(xi)为点xi处变量的实测值；Z(xi+h)为与点xi偏离h 处变量的实测值。式（2）被称为实验变异函数。根据实验变差函数图的特征，需要将实验变差函数拟合成相应的理论变差函数模型，本次研究多为指数函数和球面函数模型：

指数函数模型：

球面函数模型：

式中：C0为块金值；[C0+C]为基台值；a 为变程。块金值C0表示随机变异的大小，基台值[C0+C]是半变异函数达到的极限值。块金值与基台值的比值[C0/C0+C]称为块金系数，反映了随机变异占总变异的大小，一般认为块金系数<0.25 时，空间相关性强；块金系数在0.25～0.75 时，空间相关性中等；块金系数>0.75 时，空间相关性弱。变程a 表示具有相似性质斑块的空间连续性范围，变程以内，空间变量具有空间自相关性或空间依赖性，变程以外则不存在空间依赖性[10-11]。

模型拟合的精度采用决定系数R2和残差平方和RSS 来判断，残差平方和RSS 越小、决定系数R2越大，说明理论模型对实验半变异函数的拟合效果越好[12]。

2 结果与分析

2.1 地下水埋深和矿化度地统计特征

地下水埋深和矿化度按照每5 年（即1998—2002年、2003—2007 年、2008—2012 年、2013—2017 年）分4 个阶段进行研究。计算其埋深和矿化度的统计特征值如表1 所示，变异系数小于0.1 为弱变异性，0.1～1为中等变异性，大于1 为强变异性[13]。在4 个研究阶段，地下水埋深的平均值分别为1.99、2.23、2.27、2.47 m，变异系数分别为62.31%、73.54%、83.62%、89.88%，处在中等变异程度，数据的平均值和变异系数逐渐增大，说明灌区地下水埋深及其空间变异性逐渐增大。地下水矿化度的平均值分别为4 114.65、4 279.31、4 100.30、3 814.23 mg/L，变异系数分别为98.78%、134.47%、125.81%、122.40%，呈强变异性，数据的平均值和变异系数均先增大后减小，说明灌区地下水矿化度及其空间变异性呈现先增大后减小的趋势，地下水矿化度存在较多异常值，异常高值可能是由于生物富集或局部的土壤有机质，或者是局部施肥、污染等人工因素造成。

地统计学中半变异函数的成立必须建立在本征假设或二阶平稳假设基础上，这就要求样本点的观测值必须符合正态分布或近似正态分布[14]。呈现正态分布的数据其偏度系数和峰度系数均应接近于0，由表1 可知样本数据并不符合正态分布，因为原始观测数据中的特异值相对较多，利用ArcGIS 地统计模块中的探索数据功能，结合域处理法及数据直方图对特异值进行筛选，域处理法是将数据的平均值m加减3 倍的标准差σ 作为界限，将大于或小于[m±3σ]的数据视为特异值[15]。

特异值的存在影响了半变异函数拟合的精度[16]，对于分析地下水埋深及矿化度的空间相关性有一定干扰，故在保持数据真实性的情况下，剔除了少数特异值。对剔除特异值之后的数据进行对数转换，转换之后的数据，偏度接近于0，峰度接近于0，说明转换后的数据基本接近于正态分布，可以用于半变异函数的拟合。

表1 地下水埋深和矿化度参数特征统计表 Table 1 Statistics eigenvalue of groundwater depth and salinity parameters

利用GS+软件，采用全灌区去除特异值后的198个经过对数转换后的地下水埋深数据和108 个地下水矿化度数据，选取不同的半变异函数模型对数据进行拟合，各阶段埋深及矿化度数据的最优拟合模型及相关参数汇总如表2 所示，半变异函数模型如图2、图3 所示。

表2 地下水埋深和矿化度最优拟合模型及参数表 Table 2 The optimal fitting model and parameters of groundwater depth and salinity

4 个研究阶段埋深的最优拟合模型均为指数模型，块金系数均小于0.25，说明灌区地下水埋深的空间相关性较强，影响地下水埋深变化主要为灌区的环境因素。从1998—2017 年变程逐渐增大，说明灌区地下水埋深的空间自相关性距离变长、空间连续性变强。有3 个研究阶段地下水矿化度的最优拟合模型均为指数函数，仅2008—2012 年为球面函数，由块金系数可以看出灌区1998—2007 年地下水矿化度空间相关性较强，而2008—2017 年地下水矿化度在空间上呈中等相关性。2008—2017 年块金系数的增大说明灌区地下水矿化度的空间结构性变差，受人为因素影响造成的随机变异性增强。变程逐年增大，说明灌区地下水矿化度空间自相关性距离逐渐增大，空间连续性逐渐增强，地下水矿化度在空间上分布的均匀性增强。在4 个研究阶段，地下水矿化度的块金系数皆大于地下水埋深的块金系数，由于地下水埋深受到地形、地貌、气候等大尺度因子影响，因而其随机变异较小，空间相关性较为强烈，而地下水矿化度受到土壤性质、降水、灌溉、施肥等小尺度因子影响，随机变异性较强，空间相关性没有地下水埋深那么强烈。

2003—2017 年各地下水埋深观测点的空间相关距离在15 000～117 400 m 之间，大于实际采样距离的均值13 947 m，各地下水矿化度观测点的空间相关距离在21 600～45 300 m 之间，大于实际采样距离的均值18 063 m，故各取样点均在空间变异范围内，具有代表性。1998—2002 年的地下水埋深的变程为10 200 m，地下水矿化度的变程为17 100 m，小于各自的采样距离均值，说明试验采样间距较大，为了达到反映埋深和矿化度空间结构的目的，需要增加采样数，减小采样距离[17]。模型拟合残差极小，且决定系数在0.556～0.934 之间，均为显著水平，表明所选模型可以较准确地描述其空间结构特征。

图2 地下水埋深数据半变异函数模型 Fig.2 Semi-variogram model of groundwater depth data

图3 地下水矿化度数据半变异函数模型 Fig.3 Semi-variogram model of groundwater salinity data

2.2 地下水埋深和矿化度空间分布

运用ArcGIS 软件地统计模块的普通克里金插值法对灌区4 个研究阶段的地下水埋深和矿化度进行插值如图4、图5 所示。统计不同阶段不同埋深区间、不同矿化度区间的面积所占比例如图6 所示。灌区地下水埋深在4 个研究阶段逐渐增大，埋深在1.5～2 m的低值区域面积从76.5%降低到33.6%，埋深>3 m的高值区域面积所占比例从0.6%升高到8.2%。地下水矿化度处于2 000～4 000 mg/L 中等水平的区域面积较大，1998—2017 年灌区地下水矿化度较低值区域（<2 000 mg/L）的面积逐渐扩大，所占比例从9.4%上升到22.74%后又降低到18.53%，矿化度高值区域（>5 000 mg/L）面积从2003 年以后逐年减少，所占比例从18.58%降低至13.14%。

在空间分布上，灌区西南部沿黄河附近埋深相对较浅，基本在2 m 以下，而西北部和东北部沿狼山山前相对较深，部分区域埋深可达10 m 以上，空间分布表现出明显的条带状和局部斑块状的格局。结合图1 中灌区机电井的分布，可以看出机电井分布区域与地下水埋深高值区域的分布相似，机电井多为解决由于2003 年干旱导致农业灌溉水不足问题，作为农田灌溉的取水井，随地下水不断被开采利用，在义长灌域西南部和西北部以及乌拉特灌域三湖河地区地下水埋深在2003—2017 年呈现逐渐增大趋势。解放闸灌域和永济灌域中北部、义长灌域中部、乌拉特灌域西部无机电井覆盖，因此这些区域地下水埋深无明显增大趋势，依旧保持在<2.5 m 的较浅水平。

地下水矿化度空间方向上的分布规律受地貌、古地理和构造条件的控制，由于河套灌区整体地势西南高（高程为1 042 m），西北、东南相对较低（高程为1 034、1018 m），灌区地下水水平迁移主要是由西南向西北和东南，地下水埋深逐渐变浅，蒸发浓缩作用加剧，地下水矿化度升高，因此矿化度较高的区域分布在灌区西北部和东南部地区，西南及中部局部地区地下水矿化度较低。矿化度的空间分布还受到地下水径流的强烈影响，在灌区西南部，地下水受地形影响，水力梯度较大，水循环交替作用强烈，因此矿化度较小。灌区东部，地形变缓，地下水径流相对缓慢，地下水离子组分沿水流路径逐渐富集，矿化度增大。第一含水层组下段在乌拉特灌域的西部一带底板隆起，岩层厚度较小，阻碍了地下水的流动，导致该地地下水埋深较浅，在常年土壤蒸发强烈条件下，土壤盐碱化严重，地下水矿化度较高。

图4 地下水埋深空间分布图 Fig.4 Spatial distribution of groundwater depth

图5 地下水矿化度空间分布图 Fig.5 Spatial distribution of groundwater salinity

图6 不同阶段不同地下水埋深和矿化度区间面积占比图 Fig.6 Area proportion of different groundwater depth and salinity interval in different stages

2.3 地下水埋深和矿化度的关系

地下水埋深影响着地下水总量和土壤盐分的变化，而地下水矿化度也受地下水总量和土壤盐分的影响，因此地下水矿化度的变化与地下水埋深的变化之间必然存在一定的联系。经过比较图4 和图5 地下水埋深和矿化度的空间分布图可以明显看出：地下水矿化度较高的区域地下水埋深相应的较小（如图4、图5 中的西北部及东南部），地下水矿化度相对较低的区域地下水埋深相应较大（如图4、图5 中西南部及东北部分区域），尤其体现在灌区南部永济灌域左下方的区域，矿化度低值与埋深高值呈现鲜明的对比。

利用ArcGIS 中的栅格分割功能，将典型年2002年（降雨量和引黄灌溉水量均处于多年平均水平）的地下水埋深和矿化度空间插值图像分割成48 个区域，每个区域埋深和矿化度的平均值见图7，在绝大多数区域（编号1-8，11-16，20-24，30-37，40-48）地下水埋深和矿化度呈现一高一低的分布。地下水中的离子组分主要来源于土壤中被灌溉和降雨所淋洗的盐分，还有部分来源于灌溉水中的盐分以及地下水中岩石矿物的风化溶解。灌溉水量以及降雨量的大小很大程度上决定了地下水中离子量的多少也即矿化度的大小，同时，灌溉水量及降雨量的多少也决定了其对地下水的补给程度，影响着地下水位的变化。因此灌溉量较大的地区通过灌溉淋洗到地下水的离子量大，矿化度相应较大，而由于灌溉水的入渗补给，地下水位升高，埋深相应地较小；河套灌区目前的井灌区主要分布在狼山山前和三河湖地区（狼山山前区域编号为22、31、32；三河湖区域编号为19）由于抽取地下水灌溉要求地下水矿化度相对较小，且地下水的利用会导致这些地区埋深的增大，因此这些地区地下水埋深大，矿化度小。

图7 不同分割区域地下水埋深和矿化度均值图 Fig.7 Average groundwater depth and salinity in different segmentation areas

2.4 降雨和灌溉对地下水埋深和矿化度的影响

河套灌区1998—2017 年的年降雨量及灌溉水量如图8 所示。选取降雨量较小的2005 年作为枯水年（年降雨量为77 mm），降雨量较大的2012 年为丰水年（年降雨量为285 mm），通过对比枯水年与丰水年地下水矿化度与埋深的分布情况，分析降雨对地下水埋深和矿化度的影响。2005 年和2012 年地下水埋深和矿化度的空间分布如图9 所示。不同地下水埋深区间面积所占比例如表3 所示。不同地下水矿化度区间面积所占比例如表4 所示。结合图表可知，相比2005年，2012 年地下水埋深高值区域面积有所减小，低值及中值区域面积显著增大；地下水矿化度高值区域面积明显缩小，低值区域面积明显增大。埋深<1.6 m的低值区域面积所占比例从2005 年的2.6%增加到2012 年的11.6%，埋深处于1.6～2.0 m 的较低值区间面积所占比例从2005 年的57.41%增加到2012 年的62.73%，埋深处于2.0～2.5 m 的较高值区间的面积所占比例从2005 年的34.37%减少到2012 年的19.6%。矿化度<2 000 mg/L 的低值区间面积所占比例从2005年的12.51%增加到2012 年的31.82%，矿化度>5 000 mg/L 的高值区间面积所占比例从2005 年的20.2%减少到2012 年的11.35%。

丰水年大量的降雨对灌区整体地下水的补给作用，导致了地下水位的上升，地下水得到淡化，尽管降雨过程对土壤盐分也有一定的淋洗作用，但降雨所覆盖面积为整个灌区的耕作及未耕作区域，因此对于地下水的淡化作用强于依靠灌溉水对土壤盐分淋洗而造成的地下水矿化度上升作用。而枯水年降雨量稀少，相应地引黄灌溉水增加，在引黄灌溉水对耕作土壤的淋洗过程中，土壤中的盐分随之下渗到地下水中，使地下水矿化度增大。

图8 1998—2017 年河套灌区年降雨量及灌溉水量 Fig.8 Annual rainfall and irrigation amount in Hetao Irrigated District from1998—2017

表3 枯水年与丰水年不同地下水埋深区间面积占比 Table 3 Area proportion of different groundwater depth interval in dry year and wet year %

表4 枯水年与丰水年不同地下水矿化度区间面积占比 Table 4 Area proportion of different groundwater interval salinity in dry year and wet year %

图9 枯水年与丰水年地下水埋深和矿化度空间分布图 Fig.9 Spatial distribution of groundwater depth and salinity in dry year and wet year

灌区典型观测点地下水埋深和矿化度的年内变化趋势如图10 所示。选取平水年2002 年（年降雨量为162 mm）不同月份（3、7、11 月）分析灌溉对地下水埋深和矿化度的影响，不同月份地下水埋深和地下水矿化度空间分布如图11 所示。统计灌区不同地下水埋深区间面积占比如表5 所示，不同地下水矿化度区间面积占比如表6 所示。结合图表可知，地下水埋深在年内的变化规律十分明显，埋深低值区域（<1.4 m）面积从3—7 月增长了13.73%，从7—11月增长了39.81%，埋深较高值区域（2.4～3 m）面积则从3 月的36.97%减小至11 月的3.08%。矿化度低值区域（<2 000 mg/L）的面积从3 月的23.02%减小到11 月的21.25%，矿化度处于2 000～3 000 mg/L的区域面积从3 月的26.32%增加到11 月的29.38%，矿化度较高值区域（3000～5 000 mg/L）和高值区域（>5 000 mg/L）的面积从3—7 月再到12 月仅有微小波动。

图10 典型观测点地下水埋深和矿化度年内变化 Fig.10 Annual variation of groundwater depth and salinity at typical observation sites

图11 典型年不同灌溉期地下水埋深和矿化度空间分布图 Fig.11 Spatial distribution of groundwater depth and salinity in different irrigation periods in typical year

表5 典型年不同灌溉期不同地下水埋深区间面积占比 Table 5 Area proportion of different groundwaterdepth interval in different irrigation periods in typical year %

表6 典型年不同灌溉期不同地下水矿化度区间面积占比 Table 6 Area proportion of different groundwatersalinity interval in different irrigation periods in typical year %

3 月为土壤封冻期末，地下水埋深处于全年最小时段，4 月底5 月初春灌开始，灌区引入大量的黄河水灌溉，这些引黄水通过渠系和田间入渗补给地下水，作物生育期内水分蒸腾作用强烈，灌区处于灌溉-补给-蒸发的循环过程中，盐分通过蒸腾作用聚集在土壤与浅层地下水中，因此7 月地下水埋深减小，地下水矿化度增大；9—10 月灌区进行大面积秋浇，秋浇的形式主要是在已收割作物的土地上进行大水漫灌，起到冲盐保墒的作用，没有作物吸收水分，因此在短时间内地下水得到大量补给，地下水埋深显著减小，土壤中的盐分被淋洗到地下水中，但由于灌溉水对地下水整体的补给淡化作用，地下水矿化度相较于7 月并未呈明显增大。

3 讨论

河套灌区地下水埋深逐年下降，为保证作物的正常生长以及灌区生态环境的可持续发展，地下水埋深的下降应控制在一定范围内；河套灌区地下水矿化度先增大后减小，地下水在2010—2017 年有轻微的淡化趋势。地下水埋深在年内的变化规律十分明显，而地下水矿化度年内变化规律不明显，且特异值较多，因地下水矿化度受地形、地貌、包气带岩性、地下水的径流和气候以及地下水埋深多重条件的控制，虽然蒸发作用及灌溉水的溶滤作用对地下水矿化度的影响为主导作用，但其变化受外界局部环境干扰较大，外界污染物对地下水质也会产生直接影响，因此地下水矿化度的年际变化与年内变化趋势甚微。

降雨和灌溉水入渗均会补充地下水量，因此会导致地下水埋深的减小，然而对于矿化度的影响比对地下水埋深的影响更为复杂，一方面，降雨和灌溉入渗的水量会对地下水起到淡化的作用，从而使得地下水矿化度减小；另一方面，在降雨和灌溉水入渗的过程中，会将土壤中的盐分淋洗到浅层地下水中，从而导致地下水矿化度增大，因此降雨和入渗究竟会造成地下水矿化度的增大还是减小还需要综合考虑降雨量、灌溉入渗水量的大小以及对土壤盐分淋洗量的大小。

地下水矿化度较高的区域地下水埋深相对较小，这与前人的研究结果一致[18]，利用地下水埋深与矿化度之间的这种响应关系，可以根据地下水埋深的空间变化情况来预估地下水矿化度的空间变化，从而合理调控灌区地下水位，控制地下水矿化度和土壤盐渍化程度。