生态系统服务功能关乎人类福祉, 是社会和自然可持续发展的基础, 受到众多管理者和研究学者的关注。水源涵养是一项重要的生态系统调节服务, 指生态系统对降水的拦截和维持, 能够减少降水峰值流量、调节地表径流, 同时增加可用水资源, 是水资源得以持续的基础和保障^[1—2]。鉴于此, 近年来水源涵养功能受到了学者们相当大的关注和重视^[2—3]。开展水源涵养功能定量评估和变化的驱动力分析, 有助于政府部门作出科学决策, 对于保障区域生态安全、实现水资源和社会经济的可持续发展意义重大^[4]。

对现有水源涵养相关研究进行回顾, 发现在区域尺度上, 大多数研究还停留在定量评估、空间格局关键因素识别、时空演化特征^[5—7], 以及对单一驱动因素的响应^[8—9]。而且多采用相关性分析、主成分或聚类分析等传统方法探究水源涵养变化驱动力^[10—13], 或通过情景假设模拟水源涵养量对气候和土地利用变化的响应^[14—16], 缺少从地理分异的角度对水源涵养变化的定量归因和驱动因素的空间差异性进行研究。气候和人类活动等因素对水源涵养变化的影响常常是相互作用和非线性的, 量化它们的贡献具有挑战性^{[4, 17]}。地理探测器是一种强大的、能直接量化驱动因素及其交互作用的影响的方法, 它不必受限地遵循传统统计方法的假设, 且不涉及复杂的参数设置过程^[18—19]。目前该方法已成功用于量化驱动因素对水土保持、碳储量和生境质量等生态系统服务变化的影响^{[6, 20—23]}。而地理加权回归(Geographically weighted regression, GWR)作为一种局部回归分析方法, 能更好地揭示驱动因素的区域差异, 已被应用于探讨生态系统服务的影响因素^{[20, 24]}。地理探测器和GWR在分析驱动因素方面优势互补, 但目前鲜有将两者结合应用以深入分析区域水源涵养功能变化的驱动因素的研究, 缺少从地理空间的角度探讨水源涵养功能对多个驱动因素的综合响应及其空间差异。

黑龙江省涉及大小兴安岭等重点国有林区和三江平原、松嫩平原国家重要的粮食生产基地, 是我国重要的水源涵养区和东北“大粮仓”, 在保障国家生态安全、国土安全、粮食安全以及促进经济绿色增长中具有重要的战略地位^[25]。InVEST模型产水模块基于水量平衡原理, 通过输入降水、蒸散、土壤深度、土地利用等参数可进行大时空尺度建模, 对数据和水文专业知识的要求相对较低, 可实现空间数据量化和可视化, 能在ArcGIS下进行制图和空间分析, 因此适合评估气候和土地利用变化下区域产水量和水源涵养服务功能变化^{[5, 7, 12, 16, 26]}。本研究基于InVEST模型产水模块, 在区域尺度上开展黑龙江省水源涵养功能的定量评估和时空动态分析, 通过地理探测器进行水源涵养功能变化的定量归因, 同时结合地理加权回归模型分析水源涵养变化与驱动因素之间的地理空间关系, 科学地分析水源涵养变化驱动因素的空间分异性, 为管理部门制定有针对性的生态环境保护决策和优化生态系统管理提供科学指导。

1 研究区概况

黑龙江省位于中国东北部, 121°11′—135°05′E, 43°26′—53°33′N。地势大致呈西北、北部和东南部高, 东北、西南部低, 西北部为东北—西南走向的大兴安岭山地, 北部为西北—东南走向的小兴安岭山地, 东南部是由长白山脉的张广才岭、老爷岭、太平岭和完达山组成的东北—西南走向的山地, 东部为三江平原(包括兴凯湖平原), 西部是松嫩平原^[27—29]。属寒温带与温带大陆性季风气候，春季低温干旱, 夏季温热多雨, 秋季易涝早霜, 冬季寒冷漫长, 干燥少雨。平均年降水量为250—700 mm。1月份平均气温-32—-17℃, 7月份平均气温16—23℃。地带性植被为寒温带针叶林和温带针阔叶混交林, 主要树种有兴安落叶松、樟子松、红皮云杉、红松、白桦、山杨、蒙古栎和水曲柳等^[30], 主要分布在大小兴安岭、张广才岭、老爷岭和完达山等山地^[29]。全省植被在植物区系上分别属长白植物区系、大兴安岭植物区系和蒙古植物区系。

2 方法 2.1 数据来源

本研究收集土地利用、气候、归一化植被指数(NDVI)、土壤、数字高程模型(DEM)和流域边界等基础数据进行水源涵养功能评估和变化驱动分析, 数据说明和来源详见表 1。

将以上数据转换成统一的投影坐标系, 空间分辨率设置为250m。需说明的是, 本研究中黑龙江省大兴安岭地区的边界包括了具有行政管辖权的松岭区和加格达奇。

2.2 产水量模拟

本研究采用InVEST模型产水模块评估黑龙江省水源涵养功能空间分布及其变化。该模块是一种基于水量平衡的估算方法, 某栅格单元的降水量减去实际蒸散发后的水量即为产水量^[31], 包括地表产流、土壤含水量、枯落物持水量和冠层截留量。模块根据Zhang等^[32]基于Budyko水热耦合平衡假设提出的算法计算实际蒸散。

模型主要计算公式如下：

(1)

(2)

(3)

(4)

(5)

(6)

式中, Y_xj为第j种土地利用类型栅格单元x的产水量; AET_xj为第j种土地利用类型栅格单元x的每年实际蒸散量; P_x为栅格单元x的年降水量; R_xj为无量纲的Budyko干燥指数; ω_x为自然气候—土壤性质的非物理参数; AWC_x是栅格单元x的植物可利用水量; Z为Zhang系数, 是季节常数; K_j为第j种土地利用类型的植被蒸散系数; ET₀是潜在蒸散; soil_depth是土壤深度; root_depth是根系深度; PAWC是植物可利用水量; sand%、silt%、clay%和OM%分别为土壤质地的砂粒含量、粉砂粒含量、粘粒含量和土壤有机质含量。

模型输入变量包括：降水量(栅格数据)、潜在蒸散(栅格数据)、土壤深度(栅格数据)、植物可利用水量(栅格数据)、土地利用(栅格数据)、流域(矢量数据)、流速系数、植物蒸散系数和Z参数, 所有输入数据要具有统一投影坐标系, 且栅格数据空间分辨率一致。植物可利用水量根据土壤质地的砂粒含量、粉砂粒含量、粘粒含量和土壤有机质含量计算得出。流速系数和植物蒸散系数结合相关文献和InVEST说明手册获得^[33—34], Z参数根据黑龙江省产水量情况并结合水资源公报中地表水资源量和相关文献确定为9.433。

2.3 水源涵养量估算

用InVEST模型计算产水量之后, 考虑不同土地利用类型地表径流流速系数的影响, 结合地形指数和土壤饱和导水率再计算水源涵养量^[7], 计算公式如下：

(7)

(8)

式中, WR为单位面积水源涵养量(mm); V为流速系数; TI为地形指数, 无量纲; K_s为土壤饱和导水率(cm/d); Y为产水量(mm)。D_area为集水区栅格数量, soil_depth为土壤深度(mm), P_slope为百分比坡度。

本研究在子流域层面分析水源涵养量变化对驱动因素的综合响应。以研究期间子流域水源涵养量变化值作为因变量, 例如2020年水源涵养量减去2000年水源涵养量为2000—2020年的因变量。使用ArcGIS 10.3的“空间分析(Spatial Analyst Tools)”-“区域分析(Zonal)”-“以表格显示分区统计(Zonal Statistics as Table)”对以上评估得到的全省水源涵养量栅格数据进一步处理, 获取各子流域水源涵养量变化量的平均值。

2.4 驱动因素选择

本研究从气候因素、植被因素和人为因素三个方面选择与水源涵养功能变化相关的具有代表性且易量化、易获取的驱动因素(表 2)。降水和蒸散直接影响产水量, 气温通过影响蒸散来间接影响产水量及水源涵养功能; 植被能够减缓地表径流流速, NDVI、草地覆盖率和林地覆盖率用于反映植被状况; 耕地覆盖率和城乡工矿用地覆盖率可直接反映人类活动的影响。将2000—2020年间各个驱动因素的变化值作为自变量来进行水源涵养量变化的驱动分析, 使用ArcGIS 10.3区域分析(Zonal)工具(同2.3节)计算子流域各自变量的平均值。由于地理探测器模型只能处理离散变量^[18], 需要对8个连续型自变量重分类, 以转换为离散型自变量。可利用ArcGIS中内置的自然间断法(Natural Break)进行重分类, 该方法根据数据的固有属性进行聚类, 可以减少组内方差, 增加组间方差, 是一种常用的连续型自变量离散化处理方法^{[19, 21]}。

2.5 地理探测器

本研究使用地理探测器(Geographical detector)并以子流域作为空间分异单元来分析水源涵养量变化的主要驱动因素。主要使用地理探测器模型中的因子探测器和交互探测器。因子探测器用于分析自变量对因变量的解释程度; 交互作用探测器是探索两个自变量的联合效应是否会增加、减少对因变量的解释力或对解释力没有影响^{[18, 35]}。解释力的强弱通过比较q值的大小反映^{[18, 23]}。将因变量和离散化后的自变量读入到用Excel编制的地理探测器软件(http://www.geodetector.cn/), 然后运行。q统计量的表达式如下：

式中, q表示离散化后的自变量(表 2)对因变量(水源涵养量变化量)的解释能力, 取值在0—1之间; h=1, 2, 3, …, L是自变量的分类或分层; N_h和σ_h²分别为h层的样本数和方差; N和σ²分别是总样本量(子流域的总数目)和方差。

2.6 地理加权回归

本研究利用GWR模型来识别不同驱动因素对水源涵养量变化影响的空间分布。GWR在局部尺度上建立自变量和因变量之间的回归关系, 有效避免了变量空间差异导致的误差^[20]。某个驱动因素的回归系数反映了水源涵养量变化对该驱动因素的空间响应。较高的调整R²值表示较高的解释力和模型适合度。具体以子流域为单元, 基于ArcGIS 10.3的“空间统计(Spatial Statistics Tools)”-“空间关系建模工具集(Modeling Spatial Relationships)”, 采用“地理加权回归(Geographically Weighted Regression)”工具进行分析。回归模型的表达式如下：

(10)

式中, y_i是第i个子流域的因变量(水源涵养量的变化量); (u_i, v_i)为第i个子流域的空间位置; n是子流域的数量; β₀(u_i, v_i)是截距; p是驱动因素的数量; x_ij是第i个子流域的第j个自变量(表 2); β_j(u_i, v_i)表示第i个子流域第j个驱动因子的回归系数; ε_i是误差项。

模型参数的设置：因变量字段为水源涵养量的变化量, 自变量字段为表 2中8个驱动因子的变化量(连续型)。核的类型选择“固定距离法”, 核带宽选择“AIC(最小信息准则)”。

3 结果 3.1 土地利用变化

黑龙江省土地利用以林地、耕地、草地为主, 其次是沼泽、水域和城乡工矿。近二十年来, 黑龙江省耕地面积增加19447.25 km², 林地面积增加4780.20 km², 草地面积减少5971.14 km², 水域面积减少539.77 km², 城乡工矿面积增加了4604.25 km², 沼泽和其他用地面积分别减少19210.12 km²和3110.67 km²(表 3)。

3.2 水源涵养功能时空变化

黑龙江省水源涵养量高的地方主要集中在西北部和北部的大小兴安岭山地、南部的张广才岭、老爷岭和太平岭以及东部的完达山脉, 其次是东部的三江平原, 西部松嫩平原(半干旱区)最低(图 1)。全省单位面积水源涵养量呈现先降后升的变化趋势, 2000年为13.81 mm/m², 2010年为13.03 mm/m², 2020年达到20.29 mm/m²。全省水源涵养总量2000年为6.44×10⁹ m³; 2010年为6.07×10⁹ m³, 相比2000年减少了5.68%;2020年达到9.46×10⁹ m³, 相比2010年提高了55.76%(表 4)。

整体上来看, 2000—2010年间黑龙江省大部分区域水源涵养空间变化不明显, 只有小兴安岭南部包括伊春市乌翠区、友好区、伊美区、铁力市和绥化市绥棱县平均水源涵养量明显减少(变化量≤-5 mm/m², 下同)。2010—2020年间大部分地方的水源涵养功能呈增加趋势, 特别是小兴安岭南部、三江平原西部和南部以及张广才岭、老爷岭等区域水源涵养量明显增加(变化量≥5 mm/m², 下同), 而三江平原东北部的抚远市和大兴安岭北部的漠河市、塔河县等局部区域水源涵养量明显减少。近二十年间的水源涵养变化趋势和近十年的基本一致。

从单位面积水源涵养量来看, 各土地利用类型均呈现先下降再明显上升的变化趋势, 2020年单位面积水源涵养量比2000年提升0.57 mm/m²(沼泽)— 8.51 mm/m²(林地)。2000年以来, 草地和沼泽的水源涵养总量不断减少, 而耕地、林地、城乡工矿用地和其他用地的水源涵养总量呈先减少后增加的变化趋势。根据水源涵养总量, 林地水源涵养能力最高, 其次是耕地和草地, 然后依次是城乡工矿、沼泽和其他用地, 水域最低(表 4)。

3.3 水源涵养变化的驱动分析

地理探测器结果表明, 8个驱动因子对黑龙江省水源涵养功能的变化具有显著性影响(P值＜0.05), 不同驱动因子对水源涵养变化的解释能力存在差异性。降水和蒸散的q值最大(表 5), 表明它们是决定水源涵养量变化的空间分异特征的主导因素。草地覆盖率对水源涵养量变化的影响大于林地覆盖率。耕地覆盖率对水源涵养量变化的影响大于城乡工矿用地覆盖率。

驱动因素之间的交互作用对水源涵养量变化的影响超过任何单一因素的影响, 且主要呈非线性增强(表 5)。近二十年来, 降水和蒸散之间的耦合对水源涵养量变化的影响最为显著, 其次是降水与草地覆盖率、降水与林地覆盖率的交互作用。蒸散与其他驱动因素之间的交互作用也显著影响水源涵养量的变化, 尤其是蒸散与气温、与草地、林地和耕地覆盖率的交互作用(表 5)。

3.4 水源涵养量对驱动因素响应的空间差异性

地理加权回归模型分析表明, 本研究选择的8个驱动因素可以解释水源涵养量变化的86.13%(调整R² =0.8613)。水源涵养量和驱动因素之间的地理空间关系存在空间异质性。近二十年来, 黑龙江省水源涵养量变化与降水正相关, 与实际蒸散发负相关, 并且在大兴安岭、小兴安岭北部和三江平原东北部地区相关系数较高(图 2)。水源涵养量与草地覆盖率、林地覆盖率正相关, 并且在小兴安岭、三江平原和黑龙江省南部山地相关系数较高。耕地和城乡工矿用地覆盖率与水源涵养量变化之间的关系在不同区域表现为不同的正负相关关系, 其中在大兴安岭、小兴安岭北部和三江平原东部北部地区负相关性明显(图 2)。

4 讨论

InVEST产水量模型对降水参数非常敏感, 降水变化与水源涵养量变化趋势往往直接相关^[36—39]。本研究中黑龙江省水源涵养总量在2000—2010年间减少, 而在2010—2020年间明显增加, 与降水量波动趋势一致。地理探测器结果也显示降水因子对水源涵养量变化的解释力最高, 表明水源涵养功能变化的主要驱动力是降水。国内学者在黄河流域^{[33, 37]}、三江源^[13]、青海湖流域^[40]、青藏高原^[41]和疏勒河流域^[42]等地区也得出了类似的结论。根据水量平衡原理, 降水与实际蒸散之间的平衡直接影响区域产水量和水源涵养量, 因此实际蒸散也与水源涵养功能直接相关^[12—13]。Nahib等对西塔鲁姆河流域^[17]和王辉源等对秦岭^[16]的研究指出, 蒸散是除了降水以外对水源涵养量影响最大的因子。本研究也发现蒸散对水源涵养量变化的解释力仅次于降水。大小兴安岭地区的水源涵养量变化与实际蒸散明显负相关, 特别是漠河市、塔河县、新林区、呼中区和爱辉区实际蒸散明显增加, 导致水源涵养量明显减少。因此, 近二十年来黑龙江省水源涵养功能变化的主要影响因素是气候变化。

土地利用对水源涵养量的影响小于气候因子, 影响过程更为复杂。土地利用变化会改变土壤条件、土壤侵蚀和下垫面等情况, 影响入渗、蒸散、土壤保水和地表径流等水文过程, 进而影响水源涵养^{[37, 40, 43]}。林地对降水的涵养能力较强, 由于植被结构和根系深度不同, 草地和耕地对降水的调节作用弱于林地。而林地的蒸散系数＞耕地＞草地, 导致林地产水量低于草地和耕地, 但耕地的流速系数＞林地＞草地, 最终草地的单位面积水源涵养量＞林地＞耕地。建设用地通常覆盖有沥青、水泥和混凝土, 形成不透水层, 减少了水的渗透和聚集时间, 导致涵养水源的能力低下^{[40, 44]}。近二十年黑龙江省草地面积变化强度比林地大, 耕地面积变化强度比城乡工矿用地大, 结合地理探测器q值的大小, 可知草地覆盖率的变化对水源涵养量的影响大于林地, 而耕地覆盖率的变化对水源涵养量的影响大于城乡工矿用地。土地利用变化对区域水源涵养功能的影响取决于转换面积的大小和土地覆被的类型。局部地区土地利用变化对水源涵养功能的影响可能会大于气候因子, 例如三江平原东北部的抚远市水源涵养量明显下降, 与林地转化为耕地和城乡工矿用地密切有关。

区域水源涵养功能的变化是一个多因素耦合驱动的生态过程^{[6, 45]}。近二十年, 国家在大小兴安岭重点国有林区实施了天然林资源保护一期和二期工程, 在老爷岭和张广才岭山地丘陵区进行退耕还林^[46], 工程区森林覆盖率明显提高, 黑龙江省林地面积增加, 林地水源涵养功能提升, 表明生态工程的实施取得了显著的生态效益。虽然黑龙江省生态环境总体而言向好发展, 但有研究发现由于气候变暖变干、人口增加和牲畜增长等因素, 大兴安岭地区草地退化^[47—49]。而该地区年降水量不高, 近二十年林地面积增加导致蒸散明显提高, 局部水源涵养量明显下降。因此在大兴安岭地区增加草地面积可能比增加林地面积更有利于提升水源涵养功能。此外, 大兴安岭地区水源涵养量与草地覆盖率的相关系数明显大于与林地的相关系数, 也说明了这一点。本研究发现大小兴安岭和三江平原地区水源涵养变化与耕地和城乡工矿用地覆盖率高度负相关, 而黑龙江省南部山地和三江平原地区水源涵养变化与林草覆盖率高度正相关。因此改善土地利用结构和合理恢复植被有利于提升水源涵养功能^{[2, 20, 50]}。在大小兴安岭、三江平原东北部以及老爷岭和张广才岭等水源涵养功能重要区域, 管理部门需要密切关注气候变化和人类活动, 警惕可能导致水源涵养功能降低的不利因素, 定期评估服务功能变化情况; 在老爷岭、张广才岭和三江平原东北部地区持续推进生态恢复工程, 增加自然植被, 提高林草覆盖率, 在大小兴安岭地区注重优化人类活动, 合理控制耕地和城乡工矿等建设活动强度, 以保障黑龙江省水源涵养服务的稳定性和可持续性。

5 结论

探究水源涵养功能变化的驱动因素及其在地理空间中的关系是保障水资源可持续发展的基础和前提。本研究基于InVEST产水量模块评估了黑龙江省2000—2020年水源涵养量变化, 从地理空间的角度分析了气候和土地利用变化等因素对水源涵养变化的影响。结果表明：大小兴安岭、三江平原东北部以及老爷岭和张广才岭等南部山地是黑龙江省水源涵养功能的重要区域, 近二十年来全省水源涵养量先减少后增加。气候变化是水源涵养功能变化的主要驱动因素, 只有三江平原东北部局部区域土地利用变化对水源涵养量的影响大于气候变化。草地覆盖率的变化对水源涵养量的影响大于林地, 而耕地覆盖率的变化对水源涵养量的影响大于城乡工矿用地。大兴安岭局部地区水源涵养量降低, 与草地面积减少, 林地、耕地和城乡工矿用地面积增加密切相关。在大小兴安岭地区合理控制耕地和城乡工矿等人类建设活动强度, 在老爷岭、张广才岭和三江平原东北部地区持续推进生态恢复工程, 增加林草植被, 有助于提升黑龙江省水源涵养功能、维护区域生态安全。