基于最小响应单元的地下水埋深时空分异特征研究

牛 腾 岳德鹏 于佳鑫 于 强 王朋冲 胡雅慧

(1.北京林业大学精准林业北京市重点实验室， 北京 100083； 2.北京林业大学林学院， 北京 100083)

0 引言

地下水埋深状态是水文地质要素的综合反映，它不仅反映了地下水补给、径流、排泄的随机变化过程，同时也是判断地下水是否出现环境问题及其严重程度的重要指标。地下水的变化会引起整个生态环境系统天然平衡状态的失衡和破坏。过量开采地下水将造成地下水埋深下降。浅层地下水埋深大幅度下降，原有的绿洲会变成沙漠，原有的景观会遭到破坏。地下水埋深是地下水资源量最直接的表现形式，是地下水管理最重要的控制性指标[1]。

翁牛特旗是内蒙古自治区严重缺水地区之一，缺水面积占全旗的30%以上。由于特定的自然地理和水文地质条件，地下水一直是该地区的重要水源，为地区经济的发展发挥了重要的作用[2-3]。但是，由于严重缺水，加之长期超量开采，地下水资源匮乏已成为制约该地区经济发展的生态“瓶颈”。改进水文响应单元是根据流域内土壤、降水、植被等因素划分的具有相同水文特性的最小水文单元，基于最小水文响应单元的研究方法可在水资源的合理调配、环境风险评估等领域发挥重要作用。将最小水文响应单元理念运用于地下水埋深的模拟，能将栅格单元转化为矢量单元，将局部的数据模拟至全局，并能较为有效地反映实际情况[4-7]。空间自相关法以空间关联测度为核心，研究地下水埋深的空间分布特征，发现其空间关联模式，能较好地从区域角度分析土壤中地下水埋深的空间自相关特征。

采用协同克里金插值方法，使插值结果不仅基于地下水埋深数据的空间地理特性，也对影响地下水埋深的多种因素进行拟合，并对结果进行综合模拟，但协同克里金插值结果具有一定的范围局限性[8]。利用水文响应单元理论将空间内的栅格单元转化为矢量单元，使之更加贴近实际、符合地下水埋深的扩散和汇流状态。本文采用协同克里金插值方法，以2005—2017年每隔4 a、共4 a地下水埋深数据为主变量，对应年份的NDVI(归一化植被指数)、降水量和河网密度为协变量，将地下水埋深在研究区空间范围内呈现。利用空间自相关的方法将赋予地下水埋深插值结果的最小响应单元进行空间分析[9-10]，研究翁牛特旗时间、空间双尺度上的分异规律及特性，以期对地下水资源的合理开发利用、当地土地规划利用和农业发展规划提供参考。

1 材料与方法

1.1 研究区概况

翁牛特旗地处内蒙古自治区赤峰市中部，南与敖汉旗和赤峰市松山区相连，西接克什克腾旗，北与林西县、巴林右旗、阿鲁科尔沁旗为邻，东与通辽市奈曼旗接壤。位于东经117°49′～120°43′、北纬42°26′～43°25′之间，西部是高寒台地，中部是起伏平缓、山川交错的丘陵区；东部为科尔沁沙地。全旗既有沙海和草甸相间的低洼沼泽，也有老哈河与西拉沐沦河交汇地带流水冲击形成的平原，构成全旗“五沙四山一分田”的复杂地形。全旗总面积11 889 km2。研究区属于温带半干旱大陆性季风气候。翁牛特旗历年平均降水量在300～400 mm之间。西部山区降水量为400 mm，中部丘陵区降水量360 mm，东部荒漠区降水量300 mm。区内地下水时空分布极不均匀，主要通过大气降水、农田灌溉、干渠入渗和水库渗漏等进行补给[11]。

1.2 数据来源与处理

地下水埋深样点共40个，分布于翁牛特旗监测井，如图1所示。监测井均从2005年开始观测，观测至2017年12月，有较为完整的资料。本文选取2005、2009、2013、2017年夏季地下水埋深数据作为空间插值的原始数据。选取翁牛特旗夏季且少云的2005年7月和2009年7月Landsat TM 影像和2013年7月和2017年7月的Landsat OLI影像，利用ENVI 5.1软件对影像进行波段合成、图像增强和几何校正处理，通过建立解译标志与实地调查得到翁牛特旗2005、2009、2013、2017年土地利用现状图以及NDVI分布图。从空间地理数据云下载翁牛特旗DEM，用于水文分析中水文响应单元的划分。区域降水量根据巴林左旗、林西、翁牛特旗、宝国图和赤峰5个降水测站点的数据插值获取。根据全国分级水网数据，通过ArcGIS软件密度分析工具，构建翁牛特旗水网密度。

图1 研究区区划及地下水测站分布图Fig.1 Division of study area and distribution map of groundwater stations

1.3 协同克里金插值模型

如图1所示，地下水埋深是以点数据的形式分布于研究区内，为研究整个空间范围内地下水埋深的时空分异规律，并对其进行定性定量分析，首先需对地下水埋深数据在空间上进行呈现。因此，选取协同克里金插值，相对于其他插值方式，协同克里金插值(co-Kriging)在区域范围内进行插值的主变量与多个属性变量之间存在很强的空间相关性，协同克里金插值是将区域空间估计值的参考变量在一个主变量的基础上增加了多个区域化协同变量的方法，通过模型的拟合，分析主变量与协同变量之间的相关性[12-14]，有利于更好地估算整个区域范围内的因变量分布情况。

本文结合影响地下水埋深的静态特征(如植被覆盖程度、降水量和水网分布等因子)作为协变量，模拟地下水埋深在研究区范围内的分布情况。地表植被覆盖情况和地形的分布特点对于地下水埋深的扩散范围、扩散方式与扩散程度有重要的影响，植被覆盖情况能反映土壤的固水能力，植被覆盖率高的地方普遍地下水埋深较为稳定，多年间变化不大，影响地下水埋深的时间变化特征，且辐射能力较弱，地形则影响地下水的汇流过程，是地下水扩散的主要途径，降水是地下水径流的主要来源，水网是地下水扩散的重要载体，二者共同决定地下水埋深的空间变化特征。因此，本文选取NDVI(归一化植被指数)、区域降水量和水网密度作为协变量，对近10年间地下水埋深进行协同克里金插值。

协同区域化的协变量K所对应的影响因子数据{Zk(x)，k=1，2，…，x}，k0为主变量地下水埋深变化速率，在各地下水测站点上埋深变化速率Zk0(xi)(i=1,2,…,n)的数学期望存在并为一个定值mk0。

ZVk0对中心点x0在待估测区域V(x0)内区域化变量地下水埋深变化速率k0估测的平均值为

(1)

承载Vαk上的平均值为

(2)

(3)

式中μαk——协同克里金插值各因素权重系数

Zαk——协同化区域变量

协同化区域变量在地下水测站点{x1，x2，…,xn}的范围内符合二阶平稳假设，二阶平稳假设是指区域化变量Z(x)的任一n维分布函数不因空间点x发生位移而改变，具有数学期望存在且平稳，方差和协方差存在且平稳的性质。在本文中协同化区域变量为降水量、河网密度和NDVI，分别用ui、vj、ps表示，即

(4)

式中ai、bj、cs——降水量、河网密度和NDVI的权重系数

1.4 最小水文响应单元

依据协同克里金插值对地下水埋深的处理结果在空间上呈现为栅格数据，但在实际情况下，因多重因素空间差异和土地利用类型的不规则化，使得地下水埋深的划分单元不是类似于栅格的规则形状。因此，利用划分水文响应单元的方法，构建最小水文响应单元[15-17]，使栅格数据转换为矢量数据，使后面的空间分析更贴合于研究区的实际情况。

最小水文响应单元基于水文响应单元(Hydrological response units，HRU)，原定义为SWAT模型中模拟的基本单元，SWAT模型中的基本单元HRU通过阈值划分由土地利用、土壤类型和坡度共同定义。结合适应指数法，构建指数阈值体系，能在DEM的基础上通过水文分析更好地划分研究区的子流域。在生成子流域的基础上，依据土地利用数据作为划定标准，确定水文响应单元，并以此尺度开展相关研究[18-21]。地下水的来源主要有降水入渗、灌溉水入渗、地表水入渗补给、越流补给和人工补给。综合分析，主要由降水、地表径流和所处环境决定。地表径流入渗是地下水补给的重要途径，但本文不仅依据地表径流和地形划分的水文响应单元来作为分析单元，而是在此基础上依据研究区的地类信息，考虑不同土壤特性，不同土地利用情况的因素，并以此构建最小响应单元。因地下水的响应单元划分无法完全实地获取，相对于栅格单元和单纯地表汇流划分的水文响应单元，最小响应单元更具备科学性。

流域计算根据ArcGIS软件中的水文分析功能，利用D8单流向算法原理对DEM每个栅格单元的周围8个方向进行定义，并假定雨水降落在地形中某一个格子上，该格子的水流将会流向周围8个格子中地形最低的格子。如果多个像元格子的最大下降方向都相同，则会扩大相邻像元范围，直到找到最陡下降方向为止。由于DEM数据并不均匀分布，存在因周边像元过高产生不合理的凹陷点而带来基础数据的误差，故将不合理的凹陷点填充至合理值[22-24]。

因此，用水文响应单元来表达地下水埋深的变化规律具备可行性，将协同克里金插值结果在ArcGIS软件中利用水文分析工具对水文响应单元进行划定，水文响应单元依据研究区的地形进行划分，对基于DEM划分的流域进行二次划分，根据DEM划分的流域在空间内仅能代表地形一种特点，在同一流域内也存在多种地类，存蓄的地下水含量也各不相同。因此，利用土地覆盖数据对流域进行二次划分，并对边界的小斑块进行合并，得到最终的最小水文响应单元，每个最小水文响应单元的土地覆盖与流域划分各不相同，具备一定的独特性，以此作为地下水埋深评价的最小水文响应单元能够最大限度地分析地下水埋深的时空分异规律。

对划分好的全新最小水文响应单元进行赋值，因协同克里金插值结果是在整个研究区的栅格数据，而水文响应单元是空间数据，因此，通过分析工具中的zonal模块对插值结果统计空间范围内各水文响应单元的均值，并通过属性连接在空间上进行呈现和分析。

1.5 空间自相关

定量研究地下水的时空分异特征，需要将划分好的水文响应单元和地下水埋深协克里金插值结果耦合分析，将地下水埋深数值通过水文响应单元进行呈现，对其结果进行空间自相关分析。空间自相关是指一些变量在同一个分布区内的观测数据之间潜在的相互依赖性。空间自相关就是研究空间中，某空间单元与其周围单元间，就某种特征值，通过统计方法，进行空间自相关性程度的计算，以分析这些空间单元在空间上分布现象的特性[25-28]。描述空间自相关的参数包括全局自相关和局部自相关，本文主要通过局部自相关来表达地下水的扩散及汇流过程。

局部空间自相关是用来度量局部空间单元对于整体研究范围空间自相关的影响程度的相关性指数，代表特定水文响应单元与邻近水文响应单元的地下水埋深的相关程度，计算公式为

(5)

式中m——与水文响应单元i相邻接的水文响应单元个数

Wij——n阶矩阵的元素，代表第i个单元格与第j个单元格的空间拓扑关系

xi——第i个单元格地下水埋深

xj——第j个单元格地下水埋深

结合多年地下水埋深数据，根据最小水文响应单元划分的区域进行空间自相关分析。可得到近20年间地下水埋深的变化趋势，并以空间地块的方式进行呈现，每个最小水文响应单元分别具备不同的景观类型，分析实际情况下研究区内地下水埋深分布的聚类程度和相关情况。具备空间高聚集特性的区域，地下水汇流过程较为通畅，由结果反推过程，根据明显的自相关特性分析得出在范围内地下水流通的过程不具备较大的阻碍。

2 结果与分析

2.1 地下水埋深协同克里金插值

基于ArcGIS软件中的地统计向导工具，以2005、2009、2013、2017年地下水埋深为主变量，以对应年份夏季的NDVI、降水量和研究区内水网密度为协变量，进行协同克里金插值分析。在协同克里金插值中选择不同插值方式和参数应用，确定最好的插值效果图像。对比普通克里金插值、简单克里金插值、泛克里金插值和析取克里金插值4种插值方法，普通克里金插值多用于单个变量的无偏估计，地下水埋深结果具有较强的聚类特点，所以普通克里金插值不适用于地下水埋深研究；泛克里金插值需要知道整个插值空间的整体变化趋势，在本研究中，插值前并没有地下水变化趋势的整体估计；析取克里金插值是非线性的插值模型，在协同克里金插值结合协变量的基础上非线性插值增加了插值结果的不确定性，故此，普通克里金插值、泛克里金插值和析取克里金插值并不适用于本文的研究，研究地表形变及其影响因子的耦合关系分析，确定简单克里金插值效果最好。选择不同协方差函数进行函数拟合精度分析，环形模型、球面模型、三球模型、K-贝塞尔模型、J-贝塞尔模型和稳定模型等11个函数模型的对比如表1所示，在对协方差模型的精度验证中，将标准平均值作用选择最优模型的标准，在11个模型中稳定模型的效果最好。

表1 不同半方差函数模型插值精度评价Tab.1 Evaluation of interpolation accuracy of different semivariance function models m

图2 普通克里金插值结果Fig.2 Ordinary Kriging interpolation results

图3 协同克里金插值结果Fig.3 Co-Kriging interpolation results

对比用单一的地下水埋深作为插值因素的克里金插值(图2)，协同克里金插值结果(图3)没有突兀的点状聚集特性，插值结果在空间上更平滑，能更好体现地下水在空间上的分布规律。在整个研究区范围内，插值结果分布空间分异较为明显，其中，东西两侧的地下水埋深差距明显，且在所有年份中地下水埋深受土地利用类型影响很大。在4 a时间里，东侧地下水埋深较大，主要分布在西拉木伦河下游以胜利村为中心的农田和杜勒本吉附近的沙丘荒地；翁牛特旗中部及东部沙丘主要是平原沙丘，主要为沙漠，零星有荒草等植被覆盖，连年的干旱缺水环境和无植被根系固水的情况下，在这个范围地下水埋深较低；而大片的农田地带虽然紧邻西拉木伦河，且存在一定的农作物覆盖，但在农作物密集种植的情况下，其需水量巨大，当地居民普遍在农田周围附近打井就近灌溉，每年巨大的需水量导致在大范围农田分布区地下水埋深也普遍较深，老哈河下游农田集中区也与其类似，河流下游的农业开发活动导致地下水埋深情况不佳。地下水埋深较好的区域主要分布在西部山区，综合4 a插值结果来看，西部山区植被覆盖程度较高，高大乔灌木较多，也更为密集，植物的多样化也有利于土壤环境的维持与改善，发达的根系和良好的土壤环境有利于土壤水分的维持，每次降雨有大量的水分被植物固定在土壤中，因此，研究区东部普遍地下水埋深较浅；此外，河流附近的地下水埋深与周围相比也呈现出较好的特点，其中尤以研究区北部西拉木伦河和南部羊肠子河、红山水库最为明显，地下水埋深都是通过就近汇入河流进行扩散，在一些大型河流的附近由于小河流的汇入和水分的扩散，致使地下水埋深较高。

对比2005、2009、2013、2017年的地下水埋深插值结果，以4 a为界限，分析地下水埋深在时间上的变化规律。2005年地下水埋深在2.731 67～4.003 64 m之间，整体地下水埋深分布较为均匀，地下水埋深较大值主要分布在研究区的东北部沙丘地带，且研究区内大部分区域地下水埋深在3 m以上，说明当年地下水扩散和汇流较为顺畅，没有太多阻碍因素，地下水埋深低值沿河流分布特征明显。2009年地下水埋深在2.603 1～4.107 38 m之间，2009年地下水空间分异不明显，所有地下水埋深高值区域集中分布在东北部和南部农田范围内，且地下水埋深达到3 m以上，说明该年地下水埋深受环境影响较大，且降水较为稀少，河流区与农田区地下水埋深较为接近，农田需要过多抽取地下水进行农业灌溉。2013年地下水埋深在3.086 97～4.643 16 m之间，地下水埋深差异性很小，地下水在土壤中扩散与汇流过程通畅，不受太多阻碍，但整体地下水埋深较低，整体分布趋势与2009年类似，降水较为平均但不够充沛，地面径流起到主要调节作用。2017年地下水埋深在2.324 57～4.520 02 m之间，地下水埋深在研究区内分布较为均匀，在西拉木伦河、老哈河和羊肠子河附近地下水埋深较低，河流对于地下水的汇集有突出作用。综合以上4 a结果可得出，短期内地下水埋深的平均值不会有太大波动，但局部地区，如：东北部农田、中东部沙丘等，变化明显，且地下水埋深受当年降水等因素的影响较大，在时间序列上，地下水埋深受河流影响越来越大，地下水在土壤中的汇流过程起到了越来越重要的作用。因此，对翁牛特旗地下水资源的开发、利用与保护应充分利用这一点，以改善水资源现状，维持水环境的平衡。

2.2 最小响应单元构建

依据ArcGIS软件中的水文分析，根据DEM划分得到水文响应单元6 614个，如图4所示，但传统的水文响应单元在边界处存在大量的细碎小斑块，通过消除工具对小于75 hm2的细碎小斑块进行消除，得到水文响应单元83个。但此水文响应单元仅仅是根据研究区内地形划分的水文响应单元，地下水埋深在空间上的分区划分还受土地利用类型的影响，例如耕地、林地和建筑用地对于地下水扩散和汇流的影响各不相同，不同的土地类型具备不同的土壤类型、植被覆盖程度，土壤含水率和用水程度等多种特征的差别。因此，将依据地形划分的水文响应单元结合土地利用情况进行二次划分，划分后进行小斑消除，最终得到最小响应单元551个。最小响应单元代表不同的土地利用类型和地下水汇流方向，以此为载体在空间上呈现地下水埋深的扩散和汇流情况具有重要意义。

图4 水文响应单元分布图Fig.4 Distribution map of hydrological response units

研究区内最小响应单元的分布如图5所示，不同类型的最小响应单元具备一定的分布规律。其中林地的最小响应单元50个，均分布在翁牛特旗的西部和西南部，并存在明显的以国道为界限，国道以西以林地为主，且高程普遍在1 000 m以上，林地最小响应单元分布较为集中，且面积较大，在林地范围内，最小响应单元的聚集特点突出，范围内的地下水埋深扩散与汇流特性较为一致。研究区内水体分布也较为集中，最小响应单元主要分布在北部西拉木伦河流域以及老哈河下游的红山水库，共26个，面积均较小。耕地最小响应单元158个，面积普遍不大且分布较为零散，具体特点不显著，分布在翁牛特旗东北角，西拉木伦河下游平原区、西部林地中河流两侧、中部林草-沙漠交界带以及老哈河下游冲击平原。建筑用地最小响应单元1个，位于翁牛特旗市中心。草地最小响应单元230个，为6种地类中最多的，也是分布最为零散的一类最小响应单元，其特点是环绕沙漠地区分布。裸地最小响应单元86个，主要分布在研究区的中部和东部，其中有5块面积较大，较为集中，地下水在这类最小响应单元中扩散的规律最为明显。

图5 最小响应单元分布图Fig.5 Distribution map of minimum response unit

2.3 空间自相关分析

协同克里金插值的结果在空间上呈现为栅格数据，完成最小响应单元构建后，将协同克里金插值结果赋给最小响应单元，再利用空间统计工具中的空间自相关下的局部自相关分析4 a地下水埋深的分布及演变规律。空间自相关结果如图6所示。

图6 地下水埋深空间自相关特征Fig.6 Spatial autocorrelation characteristics of groundwater depth

由图6可知，研究区东北部的135个最小响应单元的地下水埋深结果在2005—2017年均呈现出高聚集特点，其大部分为耕地，说明在长时间序列上，耕地的地下水埋深均处在一个水位较低的区间内，且聚集特点明显，说明在大片农田范围内，地下水的扩散机制较为通畅，农作物的需水量使得这部分连片区域呈现为高值聚集区；研究区中东部的沙丘平原区也呈现出大范围的高值聚集区，在个别年份呈现出非较大值聚集特点，说明平原沙丘区的地下水埋深空间分布很不规律，不够稳定，空间内地下水埋深较低，但是存在个别斑块地下水埋深较高，地下水交互存在障碍，聚集特征不如耕地明显；研究区西部有114个最小响应单元在4 a中均为低值聚集区，西部依据林地和地形划分的最小响应单元偏大，同样面积内，低值聚集区的最小响应单元个数会少很多，由低值聚集区的分布可明显发现其主要分布于林地的范围内，而不是在水体等地下水埋深较高的最小响应单元，这说明西部山区范围内地下水埋深呈现大范围的聚集，森林中树木发达的根系能够促进地下水的产流与汇流。对比4 a地下水埋深的空间分异特性，研究区西部山区林地最小响应单元内的低值聚集区变化差异不大，但在微弱的年际变化间，低值聚集区呈现明显的聚集特点，部分低值次聚集区和低值分散区逐渐减少甚至消失，致使研究区西部的聚集特点逐年增加；在研究区东部，2005—2017年整体分布和变化规律与西部类似，但除2013年外东部的变化更加微弱，耕地的聚集特性较为稳定。这说明近年来整体地下水埋深呈现聚集特点，地下水的汇流过程呈现一种逐渐有序的状态。

3 结论

(1)以2005—2017年4期地下水埋深为主变量，以NDVI、降水量和河网密度等在空间上影响地下水埋深的数据作为协变量，采用协同克里金插值方法将地下水埋深数据在空间呈现。在4 a插值结果中，地下水埋深平均变化不大，地下水埋深整体在2～5 m之间浮动，但在中东部沙丘区变化较为明显，整体呈现出向河流方向汇集的趋势。

(2)运用水文响应单元模式，在空间上将协同克里金插值后的地下水埋深栅格数据转换为矢量数据。通过水文分析划定水文响应单元，结合研究区土地利用类型，对其重新划分为551个最小水文响应单元。在对最小水文响应单元赋值的基础上，进行空间自相关分析，4 a数据中，持续高值聚集区135个，低值聚集区114个，整体地下水埋深呈现聚集特点，地下水的汇流过程呈现一种逐渐有序的状态。

(3)对研究区内地下水埋深时空分异特性分析表明，在空间上，地下水埋深东西分异规律明显，这主要源于西部山区森林带、东部农田和沙丘带的土地利用以及地形地貌的差异，地下水位呈现西高、东低的态势，且地下水埋深情况较好的区域逐渐沿翁牛特旗的主要河流分布，受降水和河流影响逐渐变大；在时间序列上，地下水埋深平均变化不大，但逐渐呈现聚集的趋势，这说明阻碍地下水扩散的障碍正在逐渐消除。因此，建议对于翁牛特旗地下水资源的开发与利用应结合以上特点，合理规划、重点开发、有效利用。