黄河横跨我国三大地理阶梯, 构成了中国重要的生态屏障。随着人类活动干扰不断增强, 加之生态本底脆弱^[1], 水土流失^[2]、水资源短缺^[3]、土壤侵蚀严重^[4]等生态问题日趋受到学界的关注, 推动黄河流域的高质量发展已成为维护国家生态安全的重要研究课题之一。2021年10月中共中央、国务院印发《黄河流域生态保护和高质量发展规划纲要》标志着黄河流域高质量发展上升为国家战略, 纲要指出要加强对黄河流域的生态环境保护与建设, 加强黄河流域水源涵养生态服务功能的保护力度。《中华人民共和国黄河保护法》2023年4月1日正式实施也开启了黄河依法保护的新篇章, 明确要求对九曲黄河实现水资源、水环境、水生态三位一体协同治理, 推动流域高质量发展, 实现人与自然和谐共生。鉴于此, 开展黄河流域水源涵养服务功能定量评估及驱动因素探究对实现黄河流域水资源、水生态的可持续发展具有重要参考价值。

水源涵养作为生态系统水量调节服务功能的核心^[5], 具有调节地表径流、影响养分循环、增加可用水资源等多项功能^[6], 具有极其重要的研究价值。水源涵养服务功能评价是实现黄河流域水资源有效管理的重要需求。由于水源涵养具有复杂性与动态性^[5], 因此准确评估黄河流域长时间序列时空尺度的水源涵养服务功能及识别影响其变化的关键因素仍是当下研究的难点与重点。目前对水源涵养服务功能的评估基本通过InVEST模型、SWAT模型、元胞自动机模型等计算得出的水源涵养量来量化实现^[7], 其中InVEST模型凭借可视化生态系统服务功能、数据量少、输入简单、评估结果精确等优点, 成为目前评估水源涵养服务功能应用最广泛的模型之一。国内外学者已将其用于评估森林^[8]、流域^[9—11]、丘陵^[12]、河谷^[13]、生态保护区^[14—15]的水源涵养服务功能。目前, 对水源涵养变化驱动机制研究大多采用相关性分析^[16]、聚类分析^[17]等传统方法, 从地理分异角度并考虑空间异质性与多因子交互作用的水源涵养变化定量归因研究较少。对黄河流域水源涵养相关研究回溯发现, 虽有部分学者对黄河全流域及部分子流域产水服务与水源涵养服务进行研究^[18—21], 但总体研究时间尺度相对较短, 长时间序列细致分析黄河全流域水源涵养服务功能演变规律相对较少。

鉴于此, 本文以黄河流域为研究区, 以1980—2020年气象、土壤、土地利用等数据为数据源, 运用InVEST模型量化近40年来黄河流域水源涵养服务功能, 使用空间自相关分析黄河流域水源涵养服务功能空间分布模式, 采用地理探测器方法对黄河流域水源涵养服务功能变化进行归因分析, 以期揭示黄河流域水源涵养功能的时空演变特征, 为推进黄河流域生态保护与高质量发展提供重要参考。

1 研究区概况

黄河流域地处32°—42°N, 96°—119°E, 全长5464 km, 发源于巴颜喀拉山脉, 途径九省汇入渤海^[22](图 1)。黄河上、中、下游分别以内蒙古自治区托克托县河口镇、河南省郑州市桃花峪为界。黄河流域地势悬殊、西高东低, 跨越三大阶梯、四个地貌单元、三个温度带, 气候类型自西向东包含高原山地气候、温带大陆性气候、温带季风气候。上游地区以兰州为分界点, 其上主要为亚寒带, 气候高原寒冷, 冰川地貌发达, 以下主要包含中温带, 平均海拔3000 m以上, 最高峰超4000 m, 河道曲折, 上游流域多分布湖泊、沼泽、草地, 水流稳定河水清, 土地覆被以草地为主, 主要地貌包含祁连山、宁夏平原、河套平原, 河湖湿地众多, 是重要的水源补给地; 中游地区主要受中温带与暖温带的控制, 由土壤侵蚀严重的黄土地貌组成, 河流多经高山峡谷, 水流湍急, 坡度大, 流经黄土高原、汾渭平原, 携带大量泥沙, 水土流失严重, 为根治水害的关键河段; 下游地区以暖温带为主体^[23], 主要地貌为滨海平原区, 海拔较低, 泥沙长期淤积形成“地上悬河”。黄河流域主要土地利用类型为草地^[18]。黄河流域年均降水量约480 mm, 总体呈现为由东南向西北逐渐减少^{[22, 24]}。

2 数据与方法 2.1 数据来源

本文流域数据、土地利用数据、高程数据均来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/), 于ArcGIS掩膜提取获得。气象数据集来源于中国气象数据网(http://www.nmic.cn/)中1980—2020年研究区域各站点月值数据, 经ArcGIS软件中的克里金插值得气温与降水数据, 以5年为时间节点, 通过彭曼公式计算潜在蒸散发数据。土壤数据于寒区旱区科学数据中心的1 ∶ 100万土壤数据库获取, 在ArcGIS中获取相应土壤属性栅格数据。数据说明和来源详见表 1。

模型其余变量包括植物蒸散系数(Kc)、最大根系深度(Maximum root depth)、流速系数(V)、经参考他人研究成果获得^[25—26]; 季节常数(Z参数)经参考他人研究并根据实际情况调整确定为3.6^[18], 模型所需植被可利用含水率(PAWC)数据计算参考相关文献^[26]。

将以上数据转换为统一的投影坐标系, 空间分辨率为1 km。

2.2 研究方法 2.2.1 InVEST模型

InVEST模型(Integrated Valuation of Ecosystem Services and Tradeoffs)是对生态系统服务提供量化和价值评估功能的平台^[27]。本文通过产水模块估算黄河流域产水量并与地形、土壤物理性质和流速结合进行校正得研究区栅格尺度水源涵养量, 方程式如下：

(1)

(2)

(3)

式中, WC、Y、TI、D_area、soil_dep、P_slope、K_s、V、S、C分别表示单位面积水源涵养量(mm)、产水量(mm)、地形指数、集水区栅格数、土层深度(mm)、百分比坡度(mm)、土壤饱和导水率(cm/d)、流速系数、砂粒、粘粒含量(%)。产水量详细算法参考InVEST模型官方用户手册。

2.2.2 空间自相关

空间自相关是通过计算某一位置方差与邻近位置方差的关系来判断其间是否存在相互依赖性的地统计学方法^[28]。本文用Moran′s I来反映黄河流域水源涵养量变化的空间关联与差异, 公式如下：

(4)

式中, I为全局莫兰指数, 取值为[-1, 1], 正值表示空间聚集, 负值表示空间分散, 0表示空间上随机分布; n为观测总数; W_ij为空间权重, 若斑块i与j相邻, 则W_ij为1, 否则为0。

2.2.3 地理探测器

地理探测器是一种探测空间分异性及定量揭示其背后驱动因素的统计学方法^[29]。本文主要采用地理探测器模型中的因子探测器与交互探测器来分析黄河流域水源涵养量变化的主要驱动因子。因子探测器探测自变量对因变量的解释程度^[29], 而交互作用探测器用于分析两自变量间的交互作用, 即两因子共同作用的联合效应是否会增加或减少对因变量的解释力, 或对因变量的影响是相互独立的^[30]。解释程度的强弱通过q值度量^[31], 表达式为：

(5)

式中, q的取值范围为[0, 1], 值越大, 表明该因子对因变量的解释力越强; h=1, 2, …, L为变量的分层; N_h、N、σ²、σ_h²分别为层h的单元数、总单元数、总方差、层h的方差。

本研究以黄河流域1980—2020年水源涵养量为因变量, 选取降水、蒸散发、高程(DEM)、坡度、土地利用5个驱动力因子运用地理探测器进行影响因子探究。鉴于地理探测器模型对输入变量的离散要求^[29], 将降水、蒸散发、高程(DEM)、坡度按ArcGIS内置自然间断法^[32](Natural Break)重分类为6类, 该方法以数据固有属性进行聚类, 增加组间方差, 减少组内方差, 为一种连续型自变量常用的离散化处理方法^[33]。

3 结果与分析 3.1 黄河流域水源涵养量时空分布特征

本文参考已有研究^[18]及模拟计算, 最终将季节常数Z值确定为3.6, 得到多年黄河流域平均产水体积为626.08×10⁸ m³, 与《黄河水资源公报》公布的1956—2000年平均水资源总量638.37×10⁸ m³相比误差最小, 表明模型运算结果较为可靠。

1980—2020年黄河流域平均水源涵养能力为33.3393 mm, 水源涵养量为265.0475×10⁸ m³, 水源涵养量变化趋势与产水量变化一致(图 2), 年际涵养总量波动较大, 波动幅度分别为57.68%、-35.80%、-23.66%、75.13%、-17.95%、-27.22% 和47.81%, 整体呈现波动上升趋势。2020年水源涵养能力较1980年上升了19.1457 mm, 水源涵养量增加152.2086亿m³。2020年水源涵养量为历年峰值, 相较2015年增加了47.81%;2000年水源涵养量较1995年减小35.80%, 为历年最低值(图 3)。

1980—2020年黄河流域各年水源涵养量在空间分布上具有一致性, 在空间尺度上表现为西南、东南高西北低的格局(图 4), 这一格局与杨洁等人^[18]研究黄河流域产水量时空分布格局一致, 表明其空间分布的准确性。上游地区包括黄河源头区到青藏高原峡谷区与草甸区以及洛河峡谷为高水源涵养分布区, 低值区域主要集中于“几”字湾西北部, 最低值位于鄂尔多斯高原西北部, 地处库布齐沙漠, 降水偏低, 植被稀少。1980—2020年黄河流域上游与中下游地区对整体水源涵养贡献度变化幅度不大, 流域东部贡献幅度整体表现出先升后降再升的态势。2000年后, 渭河与汾河下游涵养值有所下降与该区城镇化水平提高, 植被减少有关^[34]。

3.2 不同土地利用类型的水源涵养量

土地利用分布情况会影响区域内气候、降水、土壤蓄水能力进而影响黄河流域不同地类水源涵养情况。不同土地利用类型单位面积水源涵养量可反映不同地类的水源涵养能力, 单位面积水源涵养量越大表示水源涵养能力越强。植被丰富的地区水源涵养能力相对较高, 1980—2020年, 黄河流域灌木林水源涵养能力最强, 为71.6661 mm, 其次为有林地, 其年均单位面积水源涵养量为65.6988 mm, 尽管灌木林与有林地蒸散作用较强, 但林冠层截流降水保护了土壤结构, 枯枝落叶层会加强对水分的吸收, 故水源涵养能力较强; 建设用地与未利用地集中的地区水源涵养能力相对较弱, 分别为16.3685 mm与12.1308 mm, 主要原因为建设用地无植被截留降水, 降水不易吸收, 故水源涵养能力较弱。黄河流域1980—2020年不同地类平均水源涵养能力排序：灌木林>有林地>高覆盖草地>其他林地>中覆盖草地>低覆盖草地>旱地>建设用地>未利用地>水田>水域(图 5)。

不同土地利用类型的水源涵养总量不仅与其水源涵养能力有关, 也与其土地利用类型的面积有很大关系。1980—2020年黄河流域广泛分布着各种覆被类型的草地, 总占地面积382904 km², 水源涵养能力为122.653 mm, 水源涵养量为145.3425亿m³, 贡献度为52.94%;林地占总面积的38.34%, 水源涵养能力为179.9541mm, 水源涵养量为66.63亿m³, 贡献度为24.27%(表 2)。草地与林地为黄河流域水源涵养量的主要贡献地类。

3.3 黄河流域水源涵养量空间相关性

黄河流域1980—2020年水源涵养量的全局Moran′s I值分别为0.841、0.813、0.849、0.814、0.856、0.826、0.814、0.828, P值均小于0.001, 且均通过了置信水平95%检验, 表明多年来黄河流域水源涵养服务在空间分布上存在一定的空间依赖性及较强的空间正相关性, 以高-高聚集和低-低聚集类型为主, 即水源涵养量高值区与较高区靠拢, 而水源涵养量低值区与较低区相邻(图 6)。

纵观40年来“高-高”聚集型主要分布在黄河流域东南部与西南部, 与水源涵养量高值分布区一致(图 4)。“低-低”聚集型主要分布在“几”字湾北部, 地处库布齐沙漠, 水源涵养能力弱, 叠加降水量低、蒸散量高的缘故, 故而呈“低-低”聚集分布。“低-高”与“高-低”聚集类型所占比例均较小。

3.4 黄河流域水源涵养量影响因子的地理探测

本文选取降水、蒸散发、高程、坡度、土地利用5个驱动因素, 对不同年份黄河流域水源涵养量进行了地理探测。基于因子探测的结果表明：1980—2020年, 5个驱动因子对黄河流域水源涵养变化的解释力均通过显著性检验(P＜0.05)且对水源涵养量的解释力存在差异性, 从大到小依次为降水、蒸散发、高程、土地利用类型、坡度。降水q值最高、蒸散发q值次之(表 3), 二者是决定水源涵养空间变化分异的主导因素。土地利用类型也会在一定程度上影响流域内各地的水源涵养能力, 且其因子解释力随着年份变化波动增强, 但影响力不及降水, 坡度对水源涵养的影响力最弱。

交互作用探测结果表明：驱动因子间的交互作用对水源涵养量的解释力均强于任何单一因素的影响(图 7)。近四十年来(以2020年为例, 其余年份类似), 降水与土地利用类型间的耦合对水源涵养的解释力最强, 其对水源涵养空间分布的解释力为土地利用类型单因子的3.93倍; 其次为降水与高程间的交互作用, 表明在同一高程下, 由于降水量不同导致水源涵养量表现差异, 主要表现为黄河流域中部地区南北水源涵养空间分布差异巨大。鉴于降水对水源涵养能力影响最为显著, 导致降水与其他因子的交互作用解释力均大于其余因子间的交互作用。降水与土地利用类型间的交互作用为高程与土地利用类型间交互作用的1.49倍, 即在同一土地利用类型下, 降水对水源涵养量的影响要显著强于海拔。

4 讨论

本研究显示黄河流域水源涵养服务功能时空分异特征明显。从时间尺度上看, 1980—2020年黄河流域水源涵养量呈现明显波动增加态势, 呈现2个波谷3个波峰, 谷值出现在2000年与2015年, 峰值分别出现在1990年、2005年与2020年, 这与以往黄河流域的相关研究结果基本一致。方露露等^[35]研究2000—2016年间黄河流域产水服务表明产水量最大值出现在2003年, 而本文峰值出现在2005年, 些许差别可能由于本文时间选取以5年为间隔导致。从空间尺度来看, 黄河流域水源涵养服务空间分布极不均衡并呈现明显空间异质性, 且不同时期水源涵养服务空间分异格局变化较小, 整体表现具有一定规律性, 水源涵养服务较高区域主要集中在黄河流域东南部与西南部, 较低区域主要集中表现在流域西北部沙漠。此空间分布格局与杨洁等^[18]关于1995—2015年间黄河流域产水量的研究结果一致。故而在制定黄河流域生态保护政策时, 一方面应注意强化水源涵养高值区, 另一方面可针对低值区重点修复, 以期实现黄河流域水源涵养能力整体提升。

降水和蒸散发是黄河流域水源涵养量变化的重要指标。降水量变化对水源涵养量变化趋势具有明显的直接影响^[36—39], 本文地理探测器分析结果也显示降水量对水源涵养量的解释力最高。根据水量平衡原理, 水源涵养量直接受制于降水与蒸散量间的平衡^[15]。地理探测器分析结果也显示蒸散量对水源涵养量的解释能力仅次于降水。黄河流域源头区位于青藏高原东部, 地属高原山地气候, 流经湟水谷底, 又因地势作用, 多降雨, 故水源涵养量高; 流域中上游地处西北干旱区, 属温带大陆性气候, 如“几”字湾西北部鄂尔多斯高原的库布齐沙漠, 蒸散量远远大于降水量, 降水稀少, 蒸散量高, 故流域西北部水源涵养量值最低; 而流域东部, 属温带季风性湿润气候, 降雨充沛, 故而同为水源涵养量高值区。再次表明水源涵养能力很大程度上取决于气候变化。这一点与江西东江源^[14]、三江源^[15]、黑龙江省^[32]、海南岛^[40]等地区的研究结论一致。

除了气候变化外, 土地利用类型对水源涵养能力的影响过程更为复杂^[32], 下垫面严重变化会导致有效拦截降水、调节径流、净化水质、调节生态等功能降低, 从而导致水源涵养能力下降^[7]。林地与草地在黄河流域水源涵养服务功能中扮演着重要角色。因其覆盖面积大、林冠层截流降水保护土壤结构^[41], 枯枝落叶层加强了对水分的吸收而水源涵养能力强。建设用地因其一般覆盖混凝土、水泥与沥青等不透水层, 导致水分渗透减少而水源涵养能力低下^[8]。2020年黄河下游三角洲地区水源涵养量较高, 部分原因归功于近几年实施的湿地生态恢复工程^[42]。本文地理探测器分析结果也指出降水与土地利用类型的交互作用对水源涵养的解释力最强, q值达到了0.7598。因此, 当地政府应将水源涵养服务功能的保护与建设生态林地紧密结合, 同时严格贯彻落实退耕还林还草政策以维持水源涵养服务功能的稳定, 构建九曲黄河健康水生态。

InVEST模型产水模块虽为量化水源涵养提供了良好的动态评估结果, 但该模块未将地形因素考虑在内, 未能贴切反映复杂下垫面环境下的水量平衡过程, 故而存在某些不确定性^[15]。此外, 模型中输入的最大根系深度、植被蒸散系数等均参考他人研究成果获得, 可能会对模型精度存在一定影响。未来应加强对参数本地化的研究, 同时对影响因子的选取层面应对人口、GDP等社会经济因子加以考虑。

5 结论

本文基于InVEST模型产水量模块及相关系数校正模拟量化1980—2020年黄河流域水源涵养能力, 采用空间自相关分析流域水源涵养服务功能空间分布模式, 并运用地理探测器分析降水、蒸散发、高程、坡度、土地利用类型5个驱动力因子对研究区水源涵养空间分布影响。研究结论如下：

(1) 1980—2020年黄河流域水源涵养量为174.8639亿m³—378.4538亿m³, 多年平均水源涵养量为265.0475亿m³, 呈现上下游地区较高, 中游地区较低分布格局。

(2) 1980—2020年黄河流域不同地类平均水源涵养能力排序：灌木林>有林地>高覆盖草地>其他林地>中覆盖草地>低覆盖草地>旱地>建设用地>未利用地>水田>水域。

(3) 黄河流域多年平均水源涵养量在空间分布上具有空间依赖性及较强的空间正相关性, 呈现聚集分布且以低-低聚集与高-高聚集为主。

(4) 降水量是影响黄河流域水源涵养量变化的主要驱动因子; 降水量与土地利用类型的交互作用最为显著。