气候与土地利用变化对中国黄河流域水生态系统服务时空演变的相互作用 刘焕彩、 霍星宇、 李曼、 马慧强

《Land》：Interactive Effects of Climate and Land-Use Changes on the Spatiotemporal Evolution of Water Ecosystem Services in the Yellow River Basin, China Huancai Liu, Xingyu Huo, Man Li and Huiqiang Ma

　　摘要：黄河流域（YRB）正面临着来自气候变率和土地利用强度变化的日益增加的压力；然而，这些因素如何共同塑造水生生态系统服务（WES）尚未得到充分的量化。本研究使用InVEST模型评估了2000至2020年间黄河流域的水量产出、土壤保护和水净化能力。通过情景模拟和地理检测方法，研究了包括气候和土地利用变化在内的各种驱动因素对WES的影响。在2000年至2020年间，水量产出和土壤保护显示出上升趋势，分别增加了26.24%和30.19%，而氮和磷的输出量略有下降，分别为4.82%和3.08%，表明水净化能力有一定程度的改善。所有指标在东南亚地区值较高，在西北地区值较低，三个地形区域之间存在明显差异。气候变化对水量产出、土壤保护、氮输出和磷输出的贡献率分别为97.4–99.3%、94.5–98.3%、87.2–96.0%和85.7–95.2%，表明气候因素在流域尺度上的水相关生态过程中起着核心调节作用。然而，随着时间的推移，土地利用因素的贡献率增加，并对水净化能力产生了更大的影响。WES在空间上表现出显著的异质性。气候、地形及其他因素之间的相互作用显著增强了水量产出和土壤保护服务的空间变异性，而土地利用与其他驱动因素的相互作用对水净化能力的空间变异性影响更大。因此，需要制定区域差异化的政策策略。我们的研究结果为黄河流域及世界其他类似地区在水资源管理中区分气候干预和土地利用干预提供了定量依据。

1. 引言
水生生态系统服务（WESs）是指人类通过水生态系统的结构、功能和过程获得的直接和间接利益，是水生态学领域中的一种特定形式的生态系统服务[1]。作为水生态系统的核心功能[2]，水量产出、土壤保护和水净化服务对于维护区域水资源安全、增强生态屏障和确保环境稳定性具有关键作用。目前，极端气候事件的频率上升，加上中国快速的城市化进程，加剧了环境挑战，包括水资源短缺和水生态系统退化，给其功能表现带来了巨大压力。因此，迫切需要研究WES变化的影响因素和驱动机制。
WES的评估已从单一指标估计发展到多指标综合评估。早期研究主要关注量化个别服务，如水量产出和土壤侵蚀。随着生态系统服务概念框架的完善，水量产出、土壤保护和水净化被确定为水生态系统的三个核心功能。相应地，综合评估方法也从基于统计的经验模型发展到基于过程的物理模型[2]。最近，由于计算技术的广泛应用，如InVEST和SWAT等模型因能够空间表征生态系统服务的形成过程而被广泛采用[3]。例如，有研究基于InVEST模型分析了东非高地[4]、斯洛伐克的尼特拉地区[5]和印度的通加巴德拉流域[6]的WES功能特征。中国学者主要关注流域和区域尺度，评估了鄱阳湖[7]、洞庭湖[8]、太湖[9]、卢安河[10]、吴江河[11]、北盘河[12]和重庆龙溪河[13]等流域的水量产出、土壤保护和水净化服务，以及长江经济带[14]、中国西南部[15]和京津冀地区[16]等区域的水资源状况。然而，模型参数的本地化以及多项服务之间权衡和协同作用的量化是当前评估方法的主要挑战。因此，对于生态脆弱的黄河流域（YRB）来说，缺乏一种系统的方法来准确描述核心WES的时空异质性。
气候变化和土地利用变化是WES变化的核心驱动因素[17]。气候变化通过影响植被分布、水循环和营养循环来影响水生态系统[18]。相比之下，土地利用变化通过改变水和营养物质的迁移过程、土壤结构特征以及非点源污染的形成过程来影响生态系统服务功能[19]。值得注意的是，土地利用变化是短期内影响生态系统服务供给的最动态因素。不同类型的土地利用在提供生态系统服务的能力上存在根本差异，次级土地利用类别进一步放大了这些差异。同时，土地利用的集约化产生了一种“隐藏”效应，即同一土地利用类型的不同管理实践可能导致服务输出的显著变化[20]。土地利用的空间分布模式通过控制生态系统内的物质和能量流动来影响区域生态过程。先前的研究表明，不同地区的主导驱动因素存在显著差异：在干旱和半干旱地区，气候变化对水量产出服务的贡献超过60%[21,22,23]；而在城市化快速发展的地区，土地利用变化是水质退化的首要原因[21,24]。值得注意的是，气候变化和土地利用变化并非独立作用，而是通过植被-土壤-水文反馈过程形成复杂的非线性相互作用[17]。然而，目前的研究主要基于统计相关性分析来确定驱动因素，尚不足以定量解析气候和土地利用之间的相互作用以及识别阈值效应。
黄河流域（YRB）具有相对脆弱的自然生态基础，对干扰的抵抗力低且敏感性高。目前，在气候变迁和人为干扰的共同作用下，水系统中出现了生态和环境问题，如径流急剧减少、土壤侵蚀严重和水质恶化。因此，迫切需要科学评估YRB内WES的时间和空间变化及其背后的驱动力量。目前该流域的研究主要集中在生态系统服务评估（涵盖多种服务功能[25,26,27]、生态安全格局构建[28]、权衡和协同作用[29]等方面。例如，Fan等人[30]系统分析了河南省YRB的三种一般服务（碳储存、水量产出和土壤保护），得出这三种服务具有类似的空间分布模式。这些研究大大加深了我们对YRB特定区域生态系统服务功能的理解。然而，仍存在一些不足之处。从研究主题来看，大多数研究集中在一般服务或静态评估上，很少有系统性的动态评估，针对流域尺度上的三种核心水生态系统服务（水量产出、土壤保护和水净化）。从机制分析来看，大多数研究仅定性地描述了驱动因素，缺乏对气候和土地利用变化各自贡献及其交互效应的定量区分。这些不足极大地限制了我们对YRB内WES时空演变模式和因果机制的系统性理解。
基于这些认识，我们将解决两个关键科学问题：(1) 2000至2020年间YRB核心WES的时空演变特征是什么？(2) 气候变化和土地利用变化对这些变化的贡献及其交互机制是什么？为此，本研究使用InVEST模型全面评估了2000至2020年间YRB核心WES指标（水量产出服务、土壤保护功能和水净化能力）的时间和空间演变模式。此外，还利用情景模拟和地理检测方法定量分析了气候变化和土地利用变化等驱动因素对WES变化的影响程度和相互作用。本研究旨在为YRB的生态保护和高质量发展提供科学支持。

2. 研究区域和研究方法
2.1. 研究区域
黄河发源于青藏地区的巴颜哈拉山脉。以内蒙古托克托县的河口镇和河南省郑州市的桃花峪为分水岭，它从西向东流经青海、四川、甘肃、宁夏、内蒙古、山西、陕西、河南和山东九个省份，最终通过山东省垦利县注入渤海，全长5464公里，是世界上沉积物浓度最高的河流[31]（图1）。图1. 研究区域示意图。YRB（北纬32°–42°，东经96°–119°）覆盖面积795,000平方公里。由于水生态功能主要受地形因素（海拔和地形起伏）的限制，本研究采用了一种不同于传统的上游、中游和下游划分的区域划分方法。根据地形特征和县级行政边界，该流域被划分为三个区域。兰州上游的源头地区（区域I）以海拔超过4000米、年降水量超过400毫米的众多山脉以及集中在安颜哈拉山脉和巴颜哈拉山脉的广泛冰川和永久积雪为特征，这些是黄河的重要水源。主要的地表覆盖类型为草原、林地和未利用土地（在中国土地分类系统中主要是高山沙漠和裸露岩石）。兰州下游和中游地区（区域II）以黄土高原为主，海拔在1000至2000米之间，年降水量在200至600毫米之间。耕地、草原和林地是主要的土地利用类型。下游地区（区域III）海拔普遍低于500米，属于温带季风气候，年降水量超过600毫米，土地利用以耕地、建设用地（包括城市用地、农村居民点和基础设施）和水体为主。

2.2. 数据来源
本研究主要使用了五个时期的气象和土地利用数据集（2000年、2005年、2010年、2015年和2020年），以及土壤、地形、植被和社会经济数据集。这些数据集的具体信息见表1。需要注意的是，由于气象数据的可靠性较低，高海拔地区的研究结果可能存在不确定性。为了实现数据标准化并便于综合分析，所有栅格数据被重新采样为1公里×1公里的网格，并将投影坐标系统标准化为Krasovsky_1940_Albers。

2.3. 研究方法
2.3.1. 评估水生态系统服务的方法
综合生态系统服务和权衡评估（InVEST）模型（版本3.14.1 [37]）是由斯坦福大学（美国加利福尼亚州斯坦福）、大自然保护协会（TNC）（美国弗吉尼亚州阿灵顿）和世界自然基金会（WWF）（美国华盛顿特区）共同开发的开源评估系统。该模型的核心功能是模拟生态系统服务的动态演变。该模型具有明显的优势，包括驱动数据易于获取、定量评估精度高，以及评估过程和结果的清晰空间表示[37]。近年来，InVEST模型在生态系统服务研究领域的学术认可度和应用普及度方面处于领先地位[38]。因此，本研究主要基于InVEST模型来计算和分析YRB的水生态系统核心功能。该模型分别于2000年、2005年、2010年、2015年和2020年以5年为间隔运行。
(1) 水量产出
水量产出服务代表生态系统的供水能力，定义为特定区域内单位面积单位时间的水资源供应量。在InVEST模型中，年水量产出评估模块基于Budyko水-热耦合平衡原理使用栅格网格定量估算供水能力。计算公式如下 [14]：
（1）
（2）
（3）
（4）
其中，Yx 代表网格单元 x 的年平均产水量（毫米）；AETx 对应于像素 x 处的年实际蒸散量（毫米）；Px 表示像素 x 接收的年平均降水量（毫米）；ωx 表示修正后的植物可用水量与年平均降水量的比值；Rx 表示 Budyko 干燥系数；Z 是 Zhang 系数；PAWCx 表示像素 x 的植物可用水储存量（毫米）；Kx 代表植被蒸散系数；ETox 是年平均参考蒸散量（毫米）。InVEST 模型中的年水资源评估模块需要输入每种土地利用类型的根深度和 Kx 值（表 2），这些数据来源于先前的研究 [39] 以及 InVEST 模型手册中推荐的参数 [37]。此外，Zhang 系数是模型中的一个重要输入参数，根据之前的研究 [39]，该系数被设定为 3.6。表 2. 产水量计算参数表。

**土壤保护功能**
土壤保护功能指的是生态系统调节和防止土壤因风蚀、径流和其他侵蚀过程而损失的能力。根据通用土壤流失方程（USLE），InVEST 模型中的沉积物输送模块通过测量潜在土壤侵蚀与实际土壤侵蚀之间的差异来量化土壤保护效果。计算公式如下 [40]：
（5）
（6）
（7）
其中，SD 代表土壤保护能力（吨·公顷^-2）；RKLS 表示潜在土壤侵蚀量（吨·公顷^-2）；USLE 表示实际土壤侵蚀量（吨·公顷^-2）；R 表示降雨侵蚀因子；K 表示土壤可侵蚀性因子；L 和 S 分别表示地形因子（坡长因子和坡度因子）；C 表示植被覆盖因子；P 表示土壤保护措施因子。在模型中，C 和 P 是两个固定参数，用于调整土壤保护估计的实际偏差，它们的具体值基于相关文献确定 [41]（表 3）。表 3. 土壤保护计算参数表。

**水质净化能力**
水质净化能力是指通过植被、微生物和水生生物拦截和吸附径流中的营养物质（如氮和磷）来实现的水系统的自然自我净化能力。在 InVEST 模型中，营养物质输送模块计算营养物质向水体的输出，从而反映流域水质。较低的氮和磷输出表明水质净化能力更强。营养物质通过地表径流和地下水流动传输，但由于缺乏特定路径的数据，本研究主要关注地表径流。计算公式如下 [42]：
（8）
（9）
其中，Xexp,i 表示流域内单元 i 的营养物质排放量（氮和磷）（千克·公顷^-2）。loadsurf,i 表示像素 i 处的地表营养物质负荷（氮和磷）（千克·公顷^-2），通过将不同土地覆盖类型的平均营养物质负荷与径流潜力指数（由降水量表示）进行校准得出。NDRsurf,i 是网格 i 的地表营养物质输送率，根据最大营养物质滞留效率（NDR0,i）、地形指数（ICi, IC0）和校准参数 k 确定，这些参数决定了像素 i 与河流之间的关系。在模型的输入参数中，氮和磷营养物质负荷系数（load_n 和 load_p）用于计算地表营养物质负荷 loadsurf,i。此外，植被拦截效率（eff_n 和 eff_p）对应于最大营养物质滞留效率（NDR0,i）；所有上述参数均基于 InVEST 模型手册 [37] 确定（表 4）。表 4. 水质净化能力计算参数表。

**2.3.2. 场景模拟方法**
使用 5 年为时间间隔的场景模拟 [39,43]，评估了气候和土地利用变化对 YRB 中水资源服务（WES）变化的贡献和差异。为了避免单一固定基线的累积偏差，我们采用了滚动基线方法，即每个 5 年间隔的起始年的气候和土地利用类型作为该间隔的参考。每个间隔的场景设置如下表 5 所示。以 2000 年至 2005 年为例，场景 1 假设土地利用类型不变，仅发生气候变化；而场景 2 假设气候不变，仅发生土地利用变化。随后使用 InVEST 模型基于控制变量概念对每个场景设置的结果进行了比较分析。因此，场景模拟涵盖了从 2000 年到 2020 年的整个研究时期，通过一系列 5 年的滚动比较来完成。贡献率计算如下 [23]：
（10）
（11）
（12）
（13）
其中，ES1 表示参考场景下的 WES，EScli 和 ESlu 分别代表场景 1 和 2 的模拟值。ΔEScli 和 ΔESlu 分别是由气候变化和土地利用变化引起的变化。Rcli 和 Rlu 表示气候和土地利用因素对这些变化的相对贡献。表 5. 土地和气候变化场景设置。

**2.3.3. Geodetector 方法**
上述场景分析主要用于量化气候与土地利用对 WES 随时间变化的贡献，而 Geodetector（2024 年 10 月 30 日从 http://www.geodetector.cn 访问）用于识别各个因素及其交互作用对 WES 模式的空间解释能力。这两种方法解决了不同的科学问题，因此是互补的。为了识别 YRB 中 WES 空间异质性的关键驱动因素，本研究采用了“基本因素 + 新增因素”相结合的方法进行因素选择。基本因素是服务评估过程中的关键输入变量，主要是降水量和土地利用类型。在地理检测中选择因素并非“越多越好”。引入与因变量空间异质性相关性较弱的过多因素会引入“噪声”，干扰关键驱动因素的识别。本研究旨在揭示关键决定因素，而不是不加区分地增加变量数量。因此，基于自然因素和人为活动因素的组合，我们选择了具有代表性的影响因素：降水量和年平均温度作为气候因素；数字高程模型（DEM）和地形起伏作为地形因素；归一化植被指数（NDVI）作为植被因素；以及土地利用类型、地区 GDP 和人口密度 [23,40,44] 作为人为活动因素。Jenks 自然断裂法用于对上述连续因素变量进行分类。显著性测试表明，大多数因素通过了 0.05 的显著性水平，而 GDP 和人口密度在某些年份未能通过测试。受行政边界的限制，GDP 和人口密度表现出阶梯式城乡分布，使得自然断裂法难以识别有效的空间分层。这导致层内方差过大且 q 统计量不显著。通过因素检测和交互作用分析，我们研究了这些因素如何解释 WES 的空间变异性及其交互作用的强度。因素检测的表达如下 [45]：
（14）
其中，q 的范围从 0 到 1，较高值表示该因素的解释能力更强。h 表示因素的层次等级，L 表示分层总数；Nh 和 N 分别表示层次 h 内和整个地区的单元计数； 和 分别表示层次 h 内和整个地区的 WES 方差。在此框架基础上，评估了交互作用，以确定任何配对因素的组合效应是否增强了或减弱了对 WES 空间异质性的解释能力。交互作用被分为五个不同的类别（表 6）。表 6. 因素交互作用分类。

**3. 结果**
**3.1. 水生生态系统服务的时空模式**
**3.1.1. 产水量**
从 2000 年到 2020 年，YRB 的总产水量呈波动上升趋势，从 2000 年的 301.56 毫米增加到 2020 年的 380.70 毫米，增加了 26.24%（图 2a）。这表明该流域的水资源服务能力显著提高。从空间上看，产水量服务功能的排序如下：第三区 > 第一区 > 第二区，多年平均值分别为 458.87 毫米、434.50 毫米和 244.42 毫米。图 2. 2000 年至 2020 年 YRB 中水生生态系统服务的时间序列。注：(a) 产水量；(b) 土壤保护；(c) 氮输出；(d) 磷输出。产水量反映了某个地区的水资源供应，其生成能力与降水量和蒸散量等因素密切相关。如图 3a 所示，位于高海拔地区的第一区降雨量丰富，蒸散量低，总体产水量能力强，范围在 393.45 至 484.25 毫米之间。这一能力的数值在南部地区较高，在北部地区较低，南部的高海拔地区的产水量超过 800 毫米。第二区位于半湿润温带季风气候区和半干旱温带大陆气候区的过渡地带，由于降雨不足和蒸散量较高，加上黄土高原严重的土壤侵蚀，导致总体产水量在 204.36 至 323.20 毫米之间，从东南向西北逐渐减少。第三区位于平原地区，由于人口密度和工业化加速，自然植被被农田和城市发展所取代，这在一定程度上抑制了产水量。然而，由于该地区位于中国东部，降水量相对较高，其总体产水量仍然较强，空间变化较小，范围在 395.97 至 511.54 毫米之间。值得注意的是，这一结果部分受到 InVEST 产水量模块的影响，该模块通过从降水量中减去实际蒸散量来计算产水量。在具有高不透水表面的地区（例如第三区的建筑用地），尽管水文连通性和基流受损，但蒸发量减少导致产水量增加。因此，模拟出的高产水量并不一定表示水资源丰富。图 3. 2000 年至 2020 年 YRB 中水生生态系统服务的空间分布。注：(a) 产水量；(b) 土壤保护；(c) 氮输出；(d) 磷输出。图 4a 显示了产水量的时间变化特征。从 2000 年到 2005 年，第一区的产水量在大部分地区增加，尤其是在中部（主要是青海省中部），增加了 100–200 毫米。相比之下，第二区和第三区的产水量总体呈下降趋势，减少幅度在 ?100 至 0 毫米之间。从 2005 年到 2010 年，第一区和第三区的产水量普遍下降，特别是在第一区的西南部（主要是青海省西南部），减少了超过 100 毫米。与此同时，第二区的产水量呈上升趋势，特别是在东北部（主要是内蒙古自治区中部）增加了 100–200 毫米。从 2010 年到 2015 年，整个 YRB 的产水量下降，特别是在第二区的南部（主要是陕西省），减少了超过 100 毫米。从 2015 年到 2020 年，整个流域的产水量增加，特别是在第一区的南部（主要是四川省）和第二区的中部（主要是山西省），增加了超过 200 毫米。图 4. 2000 年至 2020 年 YRB 中水生生态系统服务的变化。
**3.1.2. 土壤保护**
从 2000 年到 2020 年，YRB 的土壤保护能力从 2000 年的 88.18 吨·公顷^-2 增加到 2020 年的 114.80 吨·公顷^-2，增加了 30.19%（图 2b）。这表明该流域的总体土壤保护能力显著提高。鉴于植被生长与降水模式密切相关，并直接影响土壤稳定性，第一区的波动特征在生态政策（如三江源生态保护和发展计划及甘肃省水土保持条例）和近二十年上游降水量变化的共同作用下变得尤为明显。在空间分布上，土壤保持服务功能的排名如下：第一区 > 第三区 > 第二区，多年平均值分别为171.90吨/公顷、118.91吨/公顷和61.64吨/公顷。土壤保持能力反映了一个地区抵抗侵蚀的能力，这与地形坡度、土壤结构、水流速度、植被覆盖度和人类活动等因素密切相关。如图3b所示，第一区广阔的森林和草地面积显著减少了降水的地表径流，加上有限的人类活动，使得该地区的土壤保持能力很强，范围在146.20至198.14吨/公顷之间。相比之下，第二区的土壤保持水平普遍较低，并且存在显著的空间差异。例如，第二区西北部土壤疏松、植被稀少且风速较大，由于气候条件、土壤和植被覆盖的共同影响，该地区已成为黄河流域（YRB）中“低土壤保持区”，其土壤保持能力通常低于25吨/公顷。第二区东南部的黄土高原地区受到“退耕还林”和“陕西省水土保持条例”等法规的影响，实施了大量的土壤保持项目，包括梯田和淤泥坝建设，显著提高了该地区的土壤保持能力。值得注意的是，汾河和渭河沿岸的土壤保持能力达到了150吨/公顷。第三区的人口相对密集，尽管长期的农业开发和建设用地的扩张削弱了土壤的抗侵蚀能力，但低平的地形导致了大量的沉积物积累，从而保持了107.87至124.10吨/公顷的土壤保持量。

图4b显示了土壤保持的时间变化特征。从2000年到2005年，第一区的土壤保持能力在大范围内增加，尤其是在中部地区（主要是青海省中部和甘肃省西南部），增加了超过100吨/公顷。相比之下，第二区和第三区的土壤保持能力增减幅度不大，变化范围在-50至50吨/公顷之间。从2005年到2010年，第一区和第三区的土壤保持能力下降，降幅大多在-50至0吨/公顷之间。与此同时，第二区的土壤保持能力呈现上升趋势，增幅范围在0至50吨/公顷之间。从2010年到2015年，整个黄河流域的土壤保持能力下降，特别是在第一区中部（主要是甘肃省），减少了超过50吨/公顷。从2015年到2020年，整个流域的土壤保持能力增加，尤其是在第一区中南部（主要是甘肃省）和第二区东部（主要是山西省），增幅超过了100吨/公顷。

**3.1.3 水质净化**
从2000年到2020年，黄河流域的氮和磷输出量略有下降，分别减少了4.82%和3.08%，表明在水污染控制方面取得了一定进展（图2c,d）。在空间分布上，氮和磷的输出量排名分别为第三区 > 第二区 > 第一区，多年平均氮输出量分别为6.31千克/公顷、2.72千克/公顷和1.98千克/公顷，多年平均磷输出量分别为0.65千克/公顷、0.27千克/公顷和0.24千克/公顷。氮和磷的输出量是水质净化能力的基本指标，受到自然背景条件和人类活动的影响。在空间上，氮和磷的输出量在东南部较高，在西北部较低（图3c,d）。第三区作为一个传统的农业生产区和快速社会经济发展的地区，由于过度使用化肥、土地扩张和不可渗透表面的增加，氮和磷的输出量明显高于第一区和第二区，分别为5.612千克/公顷和7.138千克/公顷以及0.591千克/公顷和0.729千克/公顷。类似地，兰州市（第一区）东北部和第二区东南部这些人口密集的城市化地区，由于快速的工业化和城市化进程导致了大量的土地扩张，成为整个流域氮和磷的高排放区。低排放区主要位于第二区西北部（降雨量少且土壤贫瘠）和第一区西南部（有广阔的森林和草地）。这些地区的农业基础较弱，社会经济发展水平较低，导致化肥施用量减少，城市化程度适中，地表硬化程度较低，这些因素有助于维持生态系统的污染物滞留能力，有效减少了氮和磷的输出量。如图4c,d所示，研究期间氮和磷输出量的变化范围分别在-5至5千克/公顷和-0.1至0.1千克/公顷之间。从2000年到2005年，第一区的氮和磷输出量增加，而第二区和第三区则减少。从2005年到2010年，第一区和第三区的氮和磷输出量减少，但第二区有所增加。从2010年到2015年，黄河流域中部的氮和磷输出量下降，尤其是第一区中部（主要是甘肃省）。从2015年到2020年，整个流域的氮和磷输出量增加，特别是在第一区中南部（主要是甘肃省）和第二区东部（主要是山西省），增幅超过了100吨/公顷。

**3.2 驱动因素对水生态系统服务的影响**
**3.2.1 气候的影响**
2000年至2020年间，气候对黄河流域不同水生态系统服务（WES）的影响显著不同（图5）。模拟结果显示，气候因素占这一时期水量变化的97.4%-99.3%。值得注意的是，情景1显示从2015年到2020年，水量增加了38.2%，表明气候变化显著增强了水生产能力。气候对土壤保持的贡献率下降了3.26%。这种下降发生在2005-2010年和2010-2015年，因为这些时期的气候特征削弱了土壤稳定能力，而2000-2005年和2015-2020年的气候波动显著改善了土壤保持服务。气候对氮和磷输出量的贡献率继续下降，分别减少了9.17%和9.98%。比较不同时期情景1和参考情景下的氮和磷输出量模拟结果显示，气候变化在过去二十年中一直提升了水质净化能力。总之，气候变化主导了所有WES的变化。虽然气候因素对水量和土壤保持的影响大于对氮和磷输出的影响，但其对这些服务的影响随时间逐渐减弱。

**3.2.2 土地利用的影响**
随着中国对生态保护的重视，如将农田转化为森林和将牧场转化为草地等举措，2000年至2020年间黄河流域内的耕作土地和未利用土地面积减少，而森林和草地覆盖面积扩大。同时，快速的 Urbanization（城市化）导致建设用地大幅扩展。这些变化加剧了土地利用变化对流域内WES的影响，观察到了明显的空间差异（图6）。就水量变化而言，土地利用的贡献率从东向西逐渐减少。值得注意的是，下游地区（如河南省和山东省）由于建设用地的扩张，土地利用的贡献率高于中上游地区。就土壤保持服务而言，宁夏和内蒙古扩大了的森林和草地面积显著增强了当地的抗侵蚀能力，使得土地利用对土壤保持功能的贡献远远超过其他省份。对于水质净化，土地利用因素显著影响了各省氮和磷输出量的变化。相比之下，其对水量和土壤保持的影响相对较小。

**3.2.3 驱动因素对水生态系统服务的影响**
**3.2.3.1 气候的影响**
2000年至2020年间，气候对黄河流域不同水生态系统服务的影响差异显著（图5）。模拟结果显示，气候因素解释了这一时期水量变化的97.4%-99.3%。值得注意的是，情景1记录从2015年到2020年水量增加了38.2%，表明气候变化显著增强了水生产能力。气候对土壤保持的贡献率下降了3.26%。这种下降是因为2005-2010年和2010-2015年的气候特征削弱了土壤稳定能力，而2000-2005年和2015-2020年的气候波动显著改善了土壤保持服务。气候对氮和磷输出量的贡献率继续下降，分别减少了9.17%和9.98%。比较不同时期情景1和参考情景下的氮和磷输出量模拟结果发现，气候变化在过去二十年里一致提升了水质净化能力。总之，气候变化主导了所有WES的变化。虽然气候因素对水量和土壤保持的影响大于对氮和磷输出的影响，但其对这些服务的影响随时间逐渐减弱。

**3.2.2 土地利用的影响**
随着中国对生态保护的重视，从2000年到2020年间，黄河流域内将农田转化为森林和将牧场转化为草地的举措减少了耕作土地和未利用土地的面积，同时森林和草地覆盖面积扩大。然而，快速的 Urbanization（城市化）也导致建设用地大幅扩展。这些变化加剧了土地利用变化对流域内WES的影响，观察到了明显的空间差异（图6）。就水量变化而言，土地利用的贡献率总体上从东向西减少。值得注意的是，下游地区（如河南省和山东省）由于建设用地的扩张，土地利用的贡献率高于中上游地区。就土壤保持服务而言，宁夏和内蒙古扩大的森林和草地面积显著增强了当地抵抗土壤侵蚀的能力，使得土地利用对土壤保持功能变化的贡献远超过其他省份。对于水质净化，土地利用因素显著影响了各省氮和磷输出量的变化。相比之下，其对水量和土壤保持的影响相对于化肥施用和建设用地扩张的影响较小。

**3.3 主要驱动因素的识别**
通过因子分析研究了不同驱动因素对黄河流域水量、土壤保持和水质净化服务空间差异的解释力度（图7）。研究发现，气候条件、地形特征和植被因素是水量服务空间变化的主要驱动因素。相比之下，人为影响的作用较小。降水分别发挥了80.5%和77.1%的最高解释力度，其次是温度。降水和地形对土壤保持服务的解释力度也相对较高。这表明气候和地形因素是水量和土壤保持服务空间异质性的关键决定因素。相比之下，人类活动对水质净化服务的影响更大。通过实施将农田转化为森林、采用精准施肥技术和改进污水处理设施等措施，土地利用因素分别解释了氮和磷输出量空间异质性的55.8%和61.3%。相比之下，污水处理设施分别解释了45.4%和52.2%。这些发现表明，在农业和城市非点源污染的共同作用下，土地利用仍然是驱动黄河流域（YRB）养分（氮和磷）输出空间异质性的主要因素。图7显示了水生态系统服务变化可归因于YRB中的关键驱动因素。注意：x1：年平均降水量；x2：年平均温度；x3：数字高程模型（DEM）；x4：地形起伏；x5：归一化差异植被指数（NDVI）；x6：土地利用类型；x7：国内生产总值（GDP）；x8：人口密度。通过交互检测方法确定了这些因素对WES空间异质性的影响（图8）。所有交互作用都显示出增强特性（二元或非线性），证实了综合因素对空间模式的控制作用比任何单一因素都要强。例如，2000年，年平均温度和DEM之间的交互作用对水资源服务空间异质性的解释能力达到了54.7%，超过了任一单独因素的最大解释能力（52.9%）。地形起伏和NDVI之间的交互作用将土壤保持服务的解释能力提高了31.7%（2020年），这比单独因素解释能力之和（14.5% + 13.6% = 28.1%）高出3.6%。图8展示了影响YRB水生态系统服务空间模式的因素交互作用。值得注意的是，降水量与其他因素之间的交互作用对水资源服务空间异质性的解释能力增强最为显著。此外，降水量与地形起伏的交互作用分别增强了土壤保持服务的空间异质性解释能力，而土地利用与其他因素的交互作用则增强了氮和磷输出空间异质性的解释能力。尽管GDP和人口密度对每种服务的单独解释能力较弱，但它们与其他因素的交互作用增强了解释能力，其q值通常超过10%。这表明这些因素通过与自然因素的协同效应影响WES的空间异质性，而不是独立作用。因此，气候和地形因素是水资源服务和土壤保持服务空间异质性的“增强器”。

4. 讨论
4.1. 空间模式和驱动机制
总体而言，本研究中确定的水资源产量、土壤保持以及氮和磷的输出显示出明显的空间模式：东南部值较高，西北部值较低（图3），这与以往的研究结果一致。例如，Fang等人[46]表明2015年黄河流域的水资源产量和净化能力从东南向西北递减。Fang等人[41]还发现，2000年至2016年间，黄河流域的高土壤保持值主要分布在西南部山区，最高值出现在陡峭的山地地区。Geng等人[47]发现，2000年至2018年间，流域内的主要土壤保持区域集中在流域中部和上游的高原地区。尽管本研究中获得的绝对数值与前述文献中的数值存在一定差异，但它们仍处于同一数量级范围内。更重要的是，包括“东南部高，西北部低”的空间梯度以及三个地形区之间的相对排名在内的主要时空模式与以往的研究结果高度一致，从而支持了本研究采用的方法的可靠性。绝对数值的差异应归因于数据来源、处理方法（如重采样分辨率和参数定位）以及模型参数设置的差异，而不是相互矛盾的结果。此外，本研究的地理检测分析与Zhu等人[40]的结论一致，表明气候和地理因素是影响水资源服务和土壤保持服务的重要因素，而土地利用类型是影响净化能力的最重要驱动因素。本研究证实气候变化主导了黄河流域水资源服务的变化，这与多尺度研究的结果高度一致。例如，Yang等人[43]将1995年至2018年间流域水资源变化量的96.68%-100.48%归因于气候变化，Wu等人[48]发现1985年至2019年期间降水对水资源的影响远大于土地利用。这些研究共同支持了气候因素在流域尺度上对水生态过程起着核心调节作用的科学观点。本研究将时间序列延长至2020年，并同时涵盖了三项核心服务，从而增强了结论的普遍性和时效性。关于土地利用的影响，研究显示土地利用对各种WES的贡献率持续上升，其对净化能力的影响随着时间的推移显著增强，这与Hou等人[21]在山西黄河流域和Wang等人[24]在吴柳河流域的研究结果相符。三项研究均指出，农业非点源污染（化肥施用）和建设用地扩张导致的土地利用变化是氮和磷输出变化的关键驱动因素。这一共同结论强调了土地利用变化对不同流域水环境的一致干扰，突显了“通过优化土地利用结构控制非点源污染”成为流域生态管理中的核心挑战。

从驱动机制的角度来看，气候变化在黄河流域WES中的主导作用源于其全面而长期的影响。对于水资源服务而言，降水是区域水资源的主要补给来源，其“东南多，西北少”的空间分布直接决定了流域内的水资源产量梯度。对于土壤保持而言，气候变化通过气候、植被和土壤之间的三维相互作用影响服务能力。极端降水事件的频率和强度增加增强了地表径流对土壤的侵蚀作用，提高了侵蚀风险。然而，气候变暖加上黄土土壤的松散质地导致土壤开裂，加速了有机物的分解，从而降低了土壤的抗侵蚀能力。干旱压力导致植被覆盖减少，根系-土壤稳定性减弱。在净化过程中，土地利用改变了氮和磷的来源（例如化肥施用），而气候则从根本上决定了氮和磷迁移的规模以及通过水文过程的净化效率。从上述分析中我们可以看出，尽管气候变化在黄河流域WES中起主导作用，但土地利用的影响也不容忽视。土地利用不仅是WES变化的直接驱动因素，其影响日益增强，而且还是气候变化影响的调节因素（例如通过植被-土壤-水文反馈过程）；因此，它在气候变化的主导作用中起着重要作用。这一点从其与其他因素的交互作用中的非线性增强作用可以得到证明。例如，2000年降水量与土地利用类型之间的交互作用解释了85.1%的水资源服务空间异质性，超过了它们单独解释能力之和（降水80.5%，土地利用3.8%，合计84.3%），表明了非线性的增强效应。

4.2. 职献与局限性
在本研究中，我们克服了以往研究的局限性，以往研究主要关注特定区域的生态系统服务（如水资源产量、土壤保持、栖息地质量、碳储存、粮食生产等）[30]。通过关注WESs，我们系统地评估了黄河流域的三大关键服务：水资源产量、土壤保持和水净化。我们建立了一个“水量-水质-水生态”的综合评估框架，提供了对流域内水生态功能进展的更全面理解。这是首次量化气候和土地利用变化如何影响整个流域各种WES的时间动态，揭示了一个新的模式：土地利用贡献率持续上升，对水质的影响日益显著。与气候变化相比，土地利用变化更受人类生产和生活方式选择以及土地利用政策的驱动，因此更具可控性。因此，考虑到对该流域水生态系统管理的研究，应进一步开展关于土地利用的研究，例如探讨同一土地利用类型内不同管理强度的影响，包括控制措施（如化肥施用量）。此外，使用地理检测工具进行的交互检测揭示了“气候-地形-土地利用”多因素协同作用机制，放大了服务的空间差异，弥补了以往单一因素分析的不足。

然而，本研究也存在一些需要未来进一步探讨的局限性。首先，尽管InVEST模型在生态系统服务研究中被广泛使用，但在描述复杂地形方面存在局限性。为了满足模型的运行要求，我们将所有数据集重采样至1公里分辨率。虽然这种方法适用于流域尺度分析，但不可避免地引入了聚合效应。例如，在山区和地形破碎严重的地区，一些重要的地形特征会丢失，这会影响与坡度相关的计算以及土壤保持估算。同样的问题也适用于水文路径和养分传输，粗分辨率可能会扭曲水流路径。此外，模型参数（如植被蒸散系数）依赖于流域范围内的平均值，忽略了可能影响模拟精度的次区域差异。模型还将水资源产量与气候变量直接关联起来，并明确指出它不区分地表流、地下流和基流，而是假设所有来自像素的水资源产量都通过这些路径之一到达目标点。这意味着本研究可能低估了土地利用对水资源产量的贡献，因为模型没有明确考虑通过特定路径（如地下水补给或渗透）的水文效应。然而，这一局限性超出了水资源产量模型的分析范围。因此，“气候起主导作用”的结论存在一定的不确定性。未来的研究必须集成土壤与水分评估工具（SWAT）模型（专门用于模拟详细水文过程）进行多模型耦合和验证，并通过敏感性分析优化参数空间区分设置。其次，本研究的情景模拟采用了既定方法来区分气候和土地利用的影响，重点关注两者的相对贡献率。因此，计算贡献率的方法忽略了变化的方向。未来的研究应进一步探讨变化方向对贡献率的影响。第三，目前的水净化评估仅量化了通过地表径流的非点源氮和磷输出。尽管这种方法能够充分捕捉磷的传输，但由于地下流动性，可能会低估氮的负荷。此外，本研究没有考虑点源污染（工业和生活污水）的影响，也没有考虑水生生物多样性对净化能力的影响。后续研究应加入地下水模块以解决氮的传输途径问题，整合点源因素，建立双源评估框架，并将这些与水质模型（如溶解氧）相结合，以扩展评估维度。虽然现有研究已经确定了气候和土地利用的主导作用，但它们通常没有量化不同土地利用转换（如将农田转换为森林或草地转换为建设用地）对服务的边际影响，也没有明确气候因素（如降水量和温度）与土地利用因素之间的协同作用阈值。未来的研究可以采用双尺度细化方法（从次流域到县级别），利用强度分析模型量化土地利用转换的服务效应，并通过阈值检测确定因素协同作用的临界阈值。第四，为了进行定量验证，我们将模拟输出与关键水文站的观测数据进行了比较（图S1）。模拟的年水资源产量与梨津站（黄河河口）的实测径流呈现出强相关性（R = 0.85，p < 0.1）。模拟的实际土壤侵蚀（USLE）与同一站点的实测沉积物负荷高度相关（R = 0.94，p < 0.05）。对于水净化，模拟的总氮输出显示出与2005-2015年卢口站测量的氨氮浓度一致的下降趋势。这些对比结果证实，尽管存在上述提到的固有不确定性，模型模拟仍然能够捕捉到黄河流域（YRB）水资源产量、土壤保护和水净化的主要时间动态变化。未来的研究应使用更高频率的数据进行更详细的逐站点校准。值得注意的是，在本研究中，我们使用了五个离散年份（2000年、2005年、2010年、2015年和2020年）来代表二十年的变化情况。虽然这5年的时间间隔能够捕捉到十年间的趋势，但它可能忽略了年际变化和极端事件（如干旱或洪水），这些事件会对水资源产量和土壤保护产生不成比例的影响。此外，InVEST模型的养分输送模块采用年度时间步长，无法捕捉到化肥施用和径流的季节性动态变化。未来的研究应考虑更高的时间分辨率（例如月度或每日），并纳入气候极端情况，以更好地捕捉水资源生态系统的全部变异范围。

4.3 政策启示
黄河流域三个区域的水资源生态系统的显著空间异质性要求采取针对不同区域的管理策略。第一区域（源头地区）具有较高的水资源产量和良好的土壤保护能力，但养分输出极低。该地区具有生态脆弱性和较高的气候敏感性（例如冰川退缩、局部城市污染等）。因此，第一区域的调整策略应优先考虑适应气候的生态保护措施，包括严格限制高山草甸和湿地的放牧活动及基础设施扩张，建立长期的冰川-永久冻土监测网络，并通过改进污水处理和限制不透水表面的增长来控制城市非点源污染。第二区域（黄土高原）的水资源产量和土壤保护能力最低，养分输出量适中。该地区的主要问题是严重的土壤侵蚀和明显的农业非点源污染。因此，第二区域的调整策略应主要侧重于持续加强土壤和水资源保护项目（如梯田建设、挡土坝修建、植树造林等），在易侵蚀地区实施精准施肥和使用控释肥料，通过将陡坡农田转换为森林或草地来优化土地利用结构，并在支流沿线建设河岸缓冲带以拦截农业污染物。第三区域（下游平原）的水资源产量和养分输出量最高。该地区的特点是城市化进程加快、不透水表面面积扩大，以及水净化压力不断增加。因此，第三区域的推荐策略应包括严格控制建设用地的扩张，推动海绵城市的开发，提升污水处理厂的排放标准，推广精准施肥和使用控释肥料，以及在黄河及其支流沿线恢复或建设湿地公园以提高自然净化能力。在流域层面，为应对气候变化导致的极端降雨和干旱风险，需要完善整个流域的防洪和抗旱工程体系。此外，还需要建立跨区域协调机制。上游地区（第一和第二区域）通过实施土壤保护和污染控制措施，可以通过水质交易或生态补偿机制获得下游第三区域的补偿。

5. 结论
基于InVEST模型、情景模拟和地理检测方法，本研究分析了2000年至2020年间黄河流域水资源产量、土壤保护和水净化能力的时空演变，并量化了气候和土地利用变化对这些水资源生态系统的单独及交互影响。我们的研究结果为区分水资源管理中的气候因素和土地利用因素提供了定量依据，有助于为黄河流域及全球类似地区制定具有针对性的生态恢复政策。这些地区对气候变化敏感，跨越多个气候带，并且人类活动频繁。主要结论如下：
(1) 2000年至2020年间，黄河流域的水资源产量和土壤保护能力总体呈上升趋势，年际波动分别为26.24%和30.19%；相比之下，氮和磷的输出量分别下降了4.82%和3.08%，表明水净化能力有所改善。从空间上看，这三个区域的服务指标在东南部较高，在西北部较低，三个地形区域之间存在明显差异。第一、第二和第三区域的多年平均水资源产量分别为434.50毫米、244.42毫米和458.87毫米；土壤保护量分别为171.90吨/公顷、61.64吨/公顷和118.91吨/公顷；氮输出量分别为1.98千克/公顷、2.72千克/公顷和6.31千克/公顷；磷输出量分别为0.24千克/公顷、0.27千克/公顷和0.65千克/公顷。
(2) 气候是黄河流域水资源生态系统变化的重要驱动因素。气候变化对水资源产量、土壤保护、氮输出和磷输出变化的贡献率分别为97.4–99.3%、94.5–98.3%、87.2–96.0%和85.7–95.2%。进一步分析表明，尽管气候变化对流域服务的影响正在减弱，但随着时间的推移，土地利用的影响却在增加，尤其是在氮和磷的输出方面表现尤为明显。这突显了未来研究中需要更加关注土地利用问题。
(3) 不同影响因素之间的相互作用增强了单个因素在解释水资源生态系统空间异质性方面的解释力。具体来说，气候和地形与其他因素结合时，会显著加剧水资源产量和土壤保护服务的空间差异。土地利用与其他因素的相互作用进一步加剧了流域内水净化能力的空间异质性。
(4) 需要制定差异化的区域政策策略。在第一区域（源头地区），应优先考虑适应气候的生态保护措施，包括限制放牧、监测冰川-永久冻土以及控制城市非点源污染。在第二区域（黄土高原），应重点推进土壤和水资源保护项目、精准施肥和土地利用结构优化以及河岸缓冲带建设。在第三区域（下游平原），管理策略应包括控制建设用地扩张、推动海绵城市的建设、提升污水处理水平、实施精准施肥和使用控释肥料，以及沿黄河及其支流恢复湿地公园。在整个流域范围内，需要建立跨区域协调机制（如水质交易或生态补偿）和完善的防洪抗旱工程体系以应对气候变化引起的极端事件。

补充材料
相关支持信息可下载自：https://www.mdpi.com/article/10.3390/land15050791/s1
图S1：使用实测数据验证水资源生态系统服务：径流、沉积物负荷和氨氮的比较。
注释：
(a) 水资源产量；
(b) 实际土壤侵蚀量；
(c) 氮输出量。
表S1：2000年至2020年不同情景下的生态系统服务情况。