基于指示Kriging法的地下水埋深与土地利用变化关系研究

2021-08-23许越越

人民珠江 2021年8期

许越越，苏涛，雷波，王蕾

(1.安徽理工大学空间信息与测绘工程学院，安徽淮南 232001；2.中国水利水电科学研究院，北京 100000；3.中国水利水电科学研究院发展战略与政策研究室，北京 100000)

作为区域生态环境的主要影响因子之一，地下水与地表植被演变之间有着复杂的关系,灌域地下水的动态变化对土地利用/覆被变化 (Land use/land cover change,LUCC)具有强烈的影响及其响应，其变化往往会引起区域整个生态环境系统的改变和转化，合理的地下水埋深的确定在一定程度上影响着地下水、土壤以及植物之间的动态平衡。

近些年,一些专家学者开始从水文地质或景观生态学等学科方向入手,对地下水在环境中的生态调控作用或对植物变化的影响进行过研究[1]。俞斌传等[2]对抚州市临川区的2005—2015年土地利用与景观格局变化进行了相关的分析。叶红梅[3]以疏勒河流域中下游区域为研究区,分析该区土地覆盖动态变化及其与地下水的相关性,发现天然植被与地下水埋深有较大的相关性。很多研究是根据行政区划来分析的，李晓岚等[4]提出了一种跨行政区化的土地利用动态度空间融合方法。目前还有很多研究以一些代表性地区为主，比如有关绿洲、草原、灌区、长江流域等地区。张喜风等[5]通过克里格插值得到敦煌绿洲1987年和2008年这两年地下水位的时空变异性特征，并结合这两年的土地利用数据，叠加分析了地下水埋深与土地利用变化之间的关系，得到了绿洲地下水下降的主要原因。

近30年来，由于一些自然因素以及人为社会经济因素的影响，河套灌区的土地利用类型已经产生了转变，对河套灌区的土地资源进行动态监测，是解决土地资源配置问题的主要方法。地下水水位的变化在维持河套灌溉区的生态环境中起着重要作用，随着中国大型灌区节水改造工程的实施，节水灌溉等措施可能会导致区域灌溉制度和土地利用类型发生改变。地下水埋深过浅就会容易发生土壤盐碱化，地下水埋深过深就会容易导致耕地荒漠化。将河套灌区土地利用与地下水方面结合，并分析地下水埋深与土地利用变化的关系，对灌区的土地利用的管理及可持续发展有着重要的意义。因此对河套灌区土地利用变化进行动态监测分析，研究节水改造后地下水的变化和土地利用之间的关系，对河套灌区的土地利用的管理及可持续发展有着重要的意义。

1 材料与方法

1.1 研究区概况

解放闸灌域位于内蒙古河套灌区上游,东南紧邻黄河，西北地处乌拉山脚，灌域属中温带高原、大陆性气候特征。灌域地处北纬40°34′～41°14′，东经106°43′～107°27′，南北长约87 km,东西宽约81 km，海拔高度为1 030～1 046 m，西南高，东北低，总体上呈三角形状。灌域总控制面积21.56万hm2，其中灌溉面积14.21万hm2，约占河套灌区总农田灌溉面积的27%，非灌溉面积7.35万hm2。灌域多年平均地下水埋深1.72 m,年平均降水138.2 mm,年平均蒸发量为2 096.4 mm,年平均风速为2～3 m/s。灌域内设有地下水日常观测井57眼,其中有27眼同步观测地下水水质。灌域有干渠3条，分干渠16条，干沟3条，分干沟12条，灌域内各级渠道和排水沟构成了灌排配套的渠系供排水网络，根据当地种植结构，灌区年灌水7次，分别是夏灌(4—6月)3次、秋灌(7—9月)3次、秋浇(10—11月)1次，11月到次年3月为冻融期。该区风大雨少,气候干燥,蒸发量大，日照时间长，地下水埋深较浅，是典型的土壤盐渍化的灌域。

1.2 数据资料及处理

本文所使用的影像数据是Landsat5卫星影像，获取时间分别为1989、2000、2012年这3个时期的7、8、9月份，该时期研究区的植被特征明显。根据解放闸灌域的自然、社会经济特征，并结合灌域的植被和土地利用特点，对三期遥感影像进行图像处理，将研究区域的土地利用类型分为耕地、林草地、建筑用地、盐荒地、水域和其他未利用地这6大类，然后再结合各种相关专题信息及实际调查资料，对分类结果进行验证和修改，最终形成三期土地利用/土地覆被分类结果图。

地下水数据来自整个河套灌区1989—2012年各年观测井的数据，研究区所需地下水数据均为孔隙潜水类型，其中包括观测井的位置数据，观测井的地下水埋深，观测井的地下水埋深在0.22～0.40 m,通过对这些数据的整理得到观测井的年平均地下水埋深数据。为了客观并真实地反映灌区地下水位的时空变异性，依据该灌域内的地下水监测资料，选择解放闸灌域及灌域周边的53眼长期观测井1989、2000、2012年的地下水埋深数据，输入 ArcGIS中形成灌域地下水观测井分布(图1)。

图1 解放闸灌域地下水观测井位置分布

1.3 分析方法

本文以河套灌区解放闸灌域为研究区域，根据研究区1989—2012年地下水埋深观测资料，利用地统计学中的半变异函数模型、指示Kriging插值方法分析1989年、2000年以及2012年3、7、11月3个不同特征季节的地下水埋深时空分布规律，并结合遥感解译生成的三期土地利用/覆被类型分布图，叠加分析得到各土地利用/覆被类型的转变趋势。

1.3.1地统计学方法

地统计学方法能够较好地描述环境的空间变异规律，揭示自然现象的空间异质性和空间格局，在生态学领域应用十分广泛，其中半变异函数模型分析与克里格插值法是最常用的地统计分析方法[6]。指示Kriging法(Indicator Kriging)是Journel提出的一种非参数估计方法，是一种对区域不确定性估计的合理性成为处理有偏数据的有力工具[7]。

本文采用SPSS 25软件进行地下水位埋深数据的描述性分析，利用K-S检验对地下水位埋深进行正态分布检验；地下水位埋深的空间分析采用地统计学方法，首先利用GS+ 9.0软件确定不同阈值条件下地下水位埋深的变异函数模型进行半方差函数分析，然后调用ArcGIS10.2中地统计分析模块建立球状模型，插值得到地下水埋深分布图，进而分析灌域地下水埋深的空间分布特征及变异性。

其中半变异函数 (Semivairogram)是一个关于数据点的变异性与数据点间距离的函数[8]。当定量描述区域的变异特征时，需要建立变异函数的理论模型，再按照实验值进行理论模型曲线的确定和最优曲线的选择。球体模型(Spherical Model)是最常用的拟合模型之一，本文正是采用球体模型对该灌域的地下水埋深进行建模，模型的公式如下：

(1)

式中R——半变异函数；C0——块金值；C——偏基台值；C0+C——基台值；a——变程；r——步长。

在地统计学中通常采用基台值、块金值和变程来描述空间异质性程度[9]。块金值是由实验误差和小于取样尺度上的随机因素共同引起的差异,较大的块金方差表明较小尺度上存在的某种因素不容忽视；基台值为区域化变量最大变异,其值越大表明变量的空间变异程度越大；变程反映变量自相关范围的大小[10]。块金值/基台值称为基底效应，用来表示空间变异的程度，该比值越高，说明由随机部分引起的空间变异程度越大；反之，则由结构性因素引起的空间变异性程度越大。当基底效应小于25% 时，变量具有很强的空间相关性；位于25% ～75%之间，变量具有较强的空间相关性；超过75% 时,空间相关性则较弱。

1.3.2土地利用/覆被时空变异分析

首先利用1989、2000和2012年这三期的土地利用/覆被图初步对比分析解放闸灌域20多年间土地利用的变化，然后运用交叉分析方法建立逐年间土地利用类型的转移矩阵，从而定量分析各土地利用类型间相互转换的方向和强度。单一土地利用类型的动态度可用来反映该土地利用类型随时间的动态变化，可用来描述土地利用变化的区域差异分析和预测未来土地利用趋势[11]。本文利用土地动态度对各土地利用类型在时间上的变化进行定量化的描述，见式(2)：

Si=[(Lb-La)/La]/T×100%

(2)

式中Si——研究时段内某一土地利用类型的动态度；La、Lb——研究初期和末期某一种土地利用类型的面积；T——研究时段长度，a。

2 结果与分析

2.1 地下水埋深随时间的变化特征

2.1.1年际变化规律

灌域内1989—2012年3、7和11月这3个典型时间段地下水埋深变化见图2。

a)年均及线性

b)3、7、11月

从图2a可以看出，1989—2012年解放闸灌域地下水埋深总体上变化较为平缓，呈波动增大趋势。其中1989—2000年均地下水埋深呈先减小后增大的趋势，地下水埋深在1.42～1.78 m变化，2000—2012年地下水埋深在1.66～2.17 m变化；从图2b中可以看出，1989—2004年灌域内3、7、11月的地下水埋深在1.39～2.01 m变化，整体上变化一致，均呈增长趋势；2004—2012年该研究区3月(非灌溉期)地下水埋深在2.02～2.75 m、7月(秋灌期)地下水埋深在1.64～2.14 m、11月(秋浇期)地下水埋深在1.16～2.01 m，处于秋灌期和秋浇期的地下水埋深较非灌溉期要浅。

2.1.2年内变化规律

利用解放闸灌域1989—2012年逐月的地下水埋深资料，统计分析其年内变化特征，见图3。

图3 解放闸灌域地下水埋深的年内变化

从图3可以看出，1989—2012年该研究区地下水埋深整体上变化趋势一致，均随着灌溉时段呈现周期性涨跌变化规律。其中1989、2000年地下水埋深变化幅度较为缓慢，2012年地下水埋深年内变化呈较为明显的双峰型曲线，地下水埋深最小值出现在5、6月，最大埋深值出现于8、9月。这是因为灌域内夏灌期的气温相对较低，蒸发量小，渠系间的引水量已经满足作物的需水量，地下水开采量小，潜水位持续上升，使地下水埋深较小；秋灌期由于气温较高，蒸发作用强烈，作物耗水量大，地下水大量开采，使水位持续下降，故8、9月埋深最大。秋浇期由于气温相对较低，蒸发量小，作物的耗水量较小，地下水埋深呈减小趋势。

2.2 地下水埋深的空间变异性分析

地下水埋深的原始样本数为53，利用SPSS 25.0进行Log变换，剔除偏离正态分布插值点后，统计中值、均值、偏度、峰值及 K-S检验值进行评价，并采用交叉验证法进行最优模型及参数选择(表1)。由表 1 可知，该模型拟合的地下水埋深精度基本满足要求，经过对数变换后，地下水位埋深具有近似正态分布特征，中值与均值近似相等，偏度和峰值均接近于0，K-S检验值均大于0.2，较好，验证精度高，可进一步进行空间插值。

表1 地下水埋深正态分布评价指标

利用 ArcGIS10.2的地统计模块计算半变异函数模型参数，并采用平均误差(M)、均方根误差(S)、平均标准误差(δ)、标准化平均误差(Δ)以及标准化均方根误差(Sδ)对模型进行检验，最终得到最佳半变异拟合模型及其相关参数(表2)。

表2 地下水埋深半变异函数模型相关参数

由表2可知:总体上，1989年地下水埋深的块金值较于2000年和2012年大，而基台值有明显的增大，说明地下水在随机尺度上产生了一定的变化，同时也在结构性尺度上发生了改变。1989年3个月份的块金值和基台值较2000、2012年变化相对较大，总体上季节性因素影响不显著。1989、2000以及2012年3个月份地下水位埋深的块金值较小，基台值相对来说较大，这说明解放闸灌域在这3年间地下水埋深总的空间异质性程度较高，而随机性因素引起的空间异质性较小，结构因素引起的空间变异程度逐渐增大。同时，地下水埋深总体上的基底效应呈现逐年下降的趋势，按照区域化变量空间相关性程度的分级标准，均小于25%，这3年地下水埋深均呈现出较强的空间相关性，但2012年比1989、2000年空间相关程度减弱，空间异质性在逐渐增强；1989—2012年变程呈现着先减小后增大的变化趋势，说明地下水埋深在中尺度上存在着较强的变异。综上可见，近20年来解放闸灌域地下水埋深的变化在随机和结构性尺度上均发生着变化，空间异质性在逐渐增强，其变异更多是由于空间自相关引起，主要是与该研究区的地形、地貌、气候、土壤类型等自然因素有关，城建、农垦、种植制度、灌排强度等人为活动等随机因素影响较小。

经多次拟合结果确定解放闸灌域内地下水埋深变异函数值的拟合模型为球形模型，由表3可知，由该模型拟合的地下水埋深精度基本满足要求，1989、2000年以及2012年3、7、11月的地下水埋深拟合效果较好，平均误差和标准化平均误差均接近于0，均方根误差和平均标准误差近似相等，标准化平方根误差接近于1，说明该半变异函数模型能较好地反映地下水埋深的空间特征，验证精度较高。

表3 半变异函数模型交叉验证结果

2.3 地下水埋深空间分布特征

利用球形模型变异函数对解放闸灌域内的地下水埋深数据进行指示克里格插值，进而分析地下水埋深在空间上的分布特征，放闸灌域在1989、2000、2012年3、7、11月的地下水埋深在空间上的分布见图4。

从图4可以看出整个灌域20多年的地下水埋深处于动态变化过程中，7月(秋灌期)以及11月(秋浇期)地下水埋深比3月(非灌溉期)要深，灌溉期水位有所下降，非灌溉期水位逐渐上升，地下水埋深随灌溉呈周期性的涨落变化；灌域内地下水埋深总体上有增大趋势，不同地区地下水埋深变幅较大，部分地区地下水埋深增幅相较其他地区偏大，除受节水措施的影响外，也与少量的井灌区域分布有关；同时该研究区不同时期地下水埋深从西南向东北逐渐变深，具有浅埋区面积增大、深埋区面积减小的特征。1989年不同时期地下水埋深在空间上的变幅不大，但从局部上看，3—11月份东北部地下水埋深呈现逐渐减小的趋势，这说明季节性的气候变化存在着一定的影响，不同灌溉期排水强度也会有所影响，中东部的地下水埋深相对于西南部较大；2000年地下水埋深总体上较1989年相比明显增大，这与2000年后各干渠月均引水量和排水方式的改变有关，2000年后，灌域内地下水开采量有所缩减，块金值变小，说明此时由地下水补给、地势地貌等结构性因素引起的空间变异程度变大。一方面受节水灌溉的影响，2000—2012年灌域内部分地下水被消耗掉，地下水对蒸散发的贡献率在一定程度呈上升趋势，其中西南部地下水埋深增大趋势尤为明显。此外，大量的渠系水运移过程中渗漏补给地下水，蒸散发消耗地下水时使得这部分水分被再次消耗利用。灌域内蒸散发量高，为满足耗水需求，部分地下水被消耗，地下水位有所下降，这是导致灌域地下水埋深持续增大的主要原因。

a)1989年3月

b)1989年7月

c)1989年11月

e)2000年7月

f)2000年11月

g)2012年3月

h)2012年7月

i)2012年11月

2.4 土地利用/覆被变化特征

1989—2012年解放闸灌域各土地利用类型总体上变化比较明显，见图5。其中1989年林草地占比最大，其次是水域和其他未利用地，建筑用地和耕地占比较小，盐荒地占比最少，土地盐渍化问题不显著；2000年林草地和未利用地有所减少，水域和耕地面积逐渐扩大，建筑用地面积也有所减少，同时盐荒地占比也呈现增大趋势，土地逐渐盐渍化；2012年林草地和水域面积均有所减小，未利用地占比增大，盐荒地面积呈增大趋势，盐渍化问题加剧。

a)1989年

b)2000年

c)2012年

2.4.1随时间变化特征分析

解放闸灌域在20多年间各土地利用类型的数量变化较为显著，6类土地利用/覆被的面积均发生了不同程度的变化(图6)。其中，耕地、建筑用地、盐荒地和其他未利用地4种地类的面积在增加，林草地的面积在逐年减少，水域的面积也在大幅度减少。根据单一土地利用类型的动态度计算出解放闸灌域6种土地利用类型的年变化率(表4)。从表4可以看出，解放闸灌域土地利用变化速度比较快，其中以耕地和盐荒地的变化速度最快，1989—2012年耕地的面积增加了270.99 km2，动态度达到6.76%，这说明该灌域是以耕地为主要用地的社会经济结构，盐荒地的面积增加了99.24 km2，动态度为9.47%，土地盐渍化问题显著；其次是林草地和水域，1989—2012年林草地总体面积减少了449.67 km2，水域减少了204.18 km2；建筑用地和其他未利用地的年变化率较小，其中建筑用地的面积增加了26.53 km2，动态度只有0.67%，说明在20多年间,该地区的城市化进程相对较慢,城镇扩张不明显，其他未利用地面积增加了257.1 km2，动态度达到2.75%。

图6 解放闸灌域3个时期土地利用/覆盖类型面积

表4 解放闸灌域土地利用/覆被随时间动态变化

2.4.2随空间变化特征分析

通过对1989、2000、2012年3期土地利用/覆被图像进行叠加分析,求得解放闸灌域内不同土地利用类型面积变化的转移矩阵(表5、6)。

表5 解放闸灌域1989—2000年土地利用/覆被转移矩阵单位:km2

表6 解放闸灌域2000—2012年土地利用/覆被转移矩阵单位:km2

通过表5、6可以看出，灌域内各土地利用类型在20多年间各土地利用类型相互之间均存在着不同程度的转化。1989—2000年该灌域内的耕地主要转化为林草地和未利用地，其中转化为林草地的面积最多，转出大于转入，退耕还林效果显著；耕地与未利用地之间也存在相互转换，且转出均大于转入,表明灌域内土地荒漠化现象日益加剧。其他未利用地主要和林草地、耕地进行转化，其中耕地面积的增加主要由未利用地转变而来,该灌域农业发展较为快速；未利用地向林草地转化的面积占比较大，转出大于转入，植被呈现着良好的生长状态与覆盖度。水域和建筑用地之间的相互转化较为紧密,水域向建筑用地的转化趋势较强，建筑用地呈扩大趋势，城镇化发展迅速。盐荒地与其他土地利用类型之间的相互转化不明显，但土地盐渍化现象仍然存在。2000—2012年该灌域的耕地主要转化为其他未利用地，且转出大于转入，耕地有所退化，转化的面积较大，土地荒漠化问题严重；林草地和水域的面积存在大幅度减小趋势，主要转化为未利用地；建筑用地面积有所增加，城镇发展较为稳定，部分转化为盐荒地，同时盐荒地又在向其他未利用地转化，说明该灌域的土地盐渍化问题较为突出。其他未利用地主要转化为耕地，且转出大于转入，说明该灌域是以农业为主的社会经济结构，耕地是主要的土地利用类型。

2.5 地下水与土地利用/覆被变化的时空变异分析

利用指示克里金法对解放闸灌域地下水埋深进行插值，得到1989、2000以及2012年的地下水埋深空间分布，见图7；解放闸灌域的地下水埋深20多年间地下水埋深的空间变化见图8。

a)1989年

b)2000年

c)2012年

a)1989—2000年

b)2000—2012年

由图8可以看出：1989—2012年地下水埋深时空变异性发生了显著的变化，地下水埋深变化最大的区域由0.34～0.81 m的东北部向0.67～1.89 m的东南部转移，地下水埋深变化最小的区域由-0.74～-0.29 m的南部向-1.53～-0.59 m的西南部和北部转移。解放闸灌域是干旱半干旱地区，地下水埋深较浅，此时地下水埋深越大，耕地退化的现象会有所减少，灌域内土地利用类型就会形成一个相对良性的转化，极大程度上有利于灌域的生态可持续发展。将1989、2000、2012年解放闸灌域的土地利用/覆被(图5)与对应插值的地下水埋深分布(图7)叠加分析可得：地下水埋深在1.82～1.99 m变化，有利于该研究区内耕地、林草地的稳定生长；从局部区域的生态环境变化来看，当地下水埋深大于2 m，各地区均存在土地退化的不良转化。解放闸灌域东北部地区的土地利用/覆被类型变化主要是其他未利用地和林草地有所减少，由于城镇化发展较为快速稳定，人口聚集，耕地与建筑用地均有所增加，水域和盐荒地变化不明显；西南部地区是典型的农林地种植区域，由于不同季节的定期灌溉，地下水埋深变化较快，加上人类社会活动的影响，土地利用类型相应发生了变化，主要土地利用类型为水域、林草地和耕地。在该阙值范围内，该灌域内的土地利用类型和地下水埋深在空间上的变化相似，并且地下水位下降趋势越大的地区，植被自然变化程度越高，耕地退化越严重，生态系统不稳定。