探索快速城市化地区生态系统服务相互作用的时空异质性：权衡/协同变化及其驱动机制

《Ecological Indicators》：Exploring spatio-temporal heterogeneity of ecosystem service interactions in rapidly urbanizing areas: Trade-offs/synergistic changes and their driving mechanisms

黄河下游区生态系统服务时空演变与驱动机制研究 摘要：基于InVEST和RUSLE模型，量化了1990-2020年黄河下游区水 yield（WY）、碳 storage（CS）、土壤 conservation（SC）和食物 production（FP）的时空特征，结合Spearman与地理加权回归（GWR）揭示权衡与协同（TOS）关系，并运用随机森林（RF）与偏最小二乘结构方程模型（PLS-SEM）解析驱动机制。结果表明：WY空间异质性显著，与CS呈强权衡（2010年峰值-0.42），而SC与FP保持相对稳定；TOS存在显著时空异质性和尺度效应，Spearman与GWR结果存在不一致；驱动因素中LUCC、Pre、DEM和PET对WY影响最显著，LUCC主导CS，DEM主导SC，NDVI和POP显著影响FP，且驱动因素作用方向与强度随时间动态变化。研究为生态管理提供了科学依据。

　　随着经济的快速发展和城市化进程的加快，黄河流域下游地区（LYRR）的自然生态系统及其提供的生态系统服务（ESs）经历了不可逆的破坏，加剧了生态保护与经济社会发展的矛盾。本研究利用1990年至2020年的多源空间数据，结合InVEST和RUSLE模型，量化了水产出（WY）、碳储存（CS）、土壤保持（SC）和粮食生产（FP）。通过Spearman相关性和地理加权回归（GWR）分析了这些服务之间的权衡与协同关系，同时应用随机森林（RF）和偏最小二乘结构方程模型（PLS-SEM）识别了驱动因素及其作用路径。研究结果揭示了WY空间格局的显著变化，而其他三种ESs则相对稳定。ESs之间的相互作用表现出显著的时空异质性和尺度效应，导致Spearman与GWR之间出现差异。其中，WY与CS之间的最强权衡关系在2010年达到-0.42***。驱动机制表明，土地利用/覆盖变化（LUCC）、年降水量（Pre）、地形高程（DEM）和潜在蒸散发（PET）主导了WY的变化；LUCC主要驱动了CS的变化；DEM对SC有显著影响；而LUCC、归一化植被指数（NDVI）和人口（POP）则对FP有重要影响。30年间，驱动因素对ESs的影响方向和强度存在显著差异。例如，在1990年，Pre（0.734***)对WY产生了最强的正向影响，而LUCC（-0.934***）对CS则有最显著的负向影响。然而，它们通过中介路径产生的间接效应仍较弱。这些发现为快速城市化地区生态管理与可持续发展提供了科学依据。

### 1. 引言

生态系统服务（ESs）是指生态系统及其生态过程为人类生存提供的条件和利益，这些服务涵盖人类从生态系统中直接或间接获得的所有优势（Gray et al., 2024）。这些服务可以分为四大类：供给服务、调节服务、支持服务和文化服务（Reid et al., 2005）。例如，水产出（WY）、碳储存（CS）、土壤保持（SC）和粮食生产（FP）等，都与人类的生存和生活质量密切相关（Carpenter et al., 2009）。ESs研究的目的在于理解这些服务对人类福祉的贡献（Chen et al., 2023）。在整个形成和发展过程中，ESs受到自然因素、社会经济发展、人类干预等多方面的影响。不同的ESs之间表现出不同的互动关系，主要体现为权衡关系，即在某些情况下存在冲突，而在其他情况下则可能实现协同（Wu and Li, 2019）。近年来，人类对自然生态系统的过度开发加剧了生态服务功能的退化。这种趋势在快速城市化和建设用地扩张的背景下变得尤为明显（Li et al., 2025；Wei et al., 2024）。因此，在快速城市化的背景下，有必要开展相关研究，以深入理解ESs的变化及其背后的驱动因素。

### 2. 材料与方法

#### 2.1. 研究区域

研究区域覆盖了中国河南省和山东省的16个地级市，地理坐标为34°11′至38°16′N，112°33′至119°18′E，总面积约为1.22×10^5 km2（图1）。下游黄河流域主要由深且肥沃的冲积平原构成，是重要的粮食生产基地和国家生态安全的重要组成部分。该地区具有温和的大陆性季风气候，年均降雨量在473至1209 mm之间。地形以平原为主，其次是山地和丘陵，海拔范围从-23米到1680米。LYRR独特的地貌特征导致历史上经常发生河岸决口和洪水，同时带来了土地退化、土壤盐碱化、水系统污染和排水不畅等问题。此外，人类的不合理开发也加剧了该地区生态系统的脆弱性和不稳定性。研究区域面临诸如人地矛盾突出、建设用地扩张、耕地侵占、水资源短缺和土地需求多样化等挑战，这对生态安全构成了严重威胁。因此，在LYRR开展ESs及其影响机制的研究，对于实现区域生态环境保护和高质量发展具有重要意义。

#### 2.2. 数据来源

本研究的数据主要包括遥感数据、土壤和地形数据、气象数据以及社会经济数据。1990年至2020年的土地利用/覆盖变化（LUCC）数据来源于中国科学院资源与环境科学数据中心，根据中国科学院LUCC的一级分类系统，将其重新划分为六类：耕地、林地、草地、建设用地、水域和未利用地。为确保跨时期动态分析的可靠性和科学性，不同时间段的数据在分类语义和空间表达上进行了严格标准化处理。粮食产量数据来源于河南省和山东省的省级统计年鉴，涵盖LYRR内16个地级市的总粮食产量，通过耕地网格的归一化植被指数（NDVI）评估粮食供应能力的空间分布。数据资源及其应用详见表1。表中所有地理信息空间数据均使用ArcGIS 10.8进行预处理，包括根据LYRR边界进行掩膜提取以统一空间范围，对1 km×1 km网格进行重采样以统一地理空间分辨率，并转换为Krasovsky_1940_Albers投影坐标系以统一空间参考。这些步骤确保了多时期数据在空间参考和像素对齐方面的完全一致性，为长期序列分析提供了可靠的基础。

#### 2.3. 方法

本研究提出了一种结合Spearman和GWR方法与机器学习和路径建模算法的研究框架，以探索快速城市化区域中的权衡与协同关系（TOS）及其驱动机制。该研究主要分为四个部分（图2）：第一部分涉及数据收集和预处理，具体细节详见第2.2节“数据来源”；第二部分使用InVEST和RUSLE模型对WY、CS、SC和FP进行定量评估，并分析其空间分布模式和演变趋势；第三部分结合Spearman和GWR揭示了TOS的时空特征及其演变特性；第四部分则通过集成RF和PLS-SEM模型深入探讨了ESs的驱动机制。

#### 2.3.1. 生态系统服务量化

本研究选择了四个ESs进行评估：WY、CS、SC和FP。选择这四个ESs主要基于以下三个考虑：首先，这些服务类型涵盖了供给、调节和支持三大类生态系统服务，能够较为全面地表征区域生态系统的多功能性；其次，这些指标对土地利用变化和人类干扰敏感，特别是在快速城市化地区；最后，每种服务在LYRR中均具有重要的生态功能和战略意义，直接关系到区域的可持续发展目标，如水资源安全、粮食安全、碳汇增强和生态屏障稳定性。因此，对这些服务的全面评估能够系统揭示城市化过程中生态系统结构和功能变化的机制。

WY的量化基于InVEST模型，其计算方法如下：

$$
WY(x) = (1 - \frac{ET(x)}{PRE(x)}) \times PRE(x)
$$

其中，WY(x)代表每年的水产出（mm），ET(x)表示每年的实际蒸散发（mm），PRE(x)表示年降水量（mm）。PET0(x)表示潜在蒸散发（mm），Cet(x)代表植被蒸散发系数，Rct(x)表示参考植被蒸散发。

PAWC（植物可用水含量）的计算基于HWSD数据和Zhou et al. (2005)的方法，具体公式如下：

$$
PAWC = 54.509 - 0.132 \times Sand - 0.003 \times Sand^2 - 0.055 \times Silt - 0.06 \times Silt^2 - 0.738 \times Clay + 0.007 \times Clay^2 - 2.688 \times OM + 0.501 \times OM^2
$$

其中，PAWC代表土壤中的植物可用水含量；Sand、Silt和Clay分别表示土壤中砂、粉砂和黏土的比例（%）；OM表示土壤有机质含量；OC表示有机碳含量，OM = 1.724 × OC。E(x)是经验参数，Z表示张系数，即所谓的“季节因子”。该参数作为调整因子，通过经验数据得出，用于考虑季节变化对WY的影响。根据Yang et al. (2020)在YRB上进行的WY多期运行结果，结合实际情况，Z值被确定为3.6。

CS的量化基于InVEST模型，该模型结合了LUCC和四个基本碳库作为输入数据。计算公式如下：

$$
C_{total} = C_{above} + C_{below} + C_{soil} + C_{dead}
$$

其中，C_total代表总CS（t/ha），C_above和C_below分别代表地上和地下生物碳库，C_soil表示土壤碳库，C_dead代表死有机碳库。碳密度值参考了Yang et al. (2021a)的研究，结合了YRB中的降水量和年平均温度，对不同土地利用类型进行了修订。

SC的估算基于RUSLE模型（Benavidez et al., 2018）。其公式如下：

$$
A = R \times K \times L \times S \times C \times P
$$

其中，A表示年均土壤侵蚀量；R表示降雨侵蚀因子；K表示土壤可蚀性因子；L和S分别代表地形因子（坡长和坡度）；C表示植被覆盖和管理因子；P表示土壤-水保护措施因子。由于潜在土壤侵蚀A_p不考虑植被和水保措施，其P和C分别取值为1。

FP的计算基于NDVI和LUCC数据，通过NDVI对耕地栅格进行空间分配。其计算公式如下：

$$
G_i = G_{sum} \times \frac{NDVI_i}{NDVI_{sum}}
$$

其中，G_i表示第i个栅格单元的作物产量，G_sum表示LYRR内的总FP，NDVI_i表示第i个栅格单元的NDVI值，NDVI_sum表示LYRR内的总NDVI值。

#### 2.3.2. 坡度趋势分析

研究使用MATLAB对1990年至2020年的WY、CS、SC和FP进行了趋势分析。计算了每个像素30年间的变动坡度。结合显著性水平（p值），将变动趋势客观地划分为五个类别：基本稳定、显著增长、轻微增长、显著下降和轻微下降（详见表2）。表中展示了坡度趋势分类的标准。

#### 2.3.3. 权衡与协同分析

本研究基于领域内常用的统计标准和相关研究中的知识，通过相关系数（r）来判断关系的强度。具体标准如下：|r| ≥ 0.6表示强相关，0.3 ≤ |r| < 0.6表示中等相关，0 < |r| < 0.3表示弱相关。系数的正负区分了影响的方向：负值表示权衡，正值表示协同。这些标准在本研究中统一应用于解释Spearman和GWR模型的结果，确保TOS分类的一致性。

Spearman相关性的计算公式如下：

$$
ρ = 1 - \frac{6 \sum d_i^2}{n(n^2 - 1)}
$$

其中，ρ表示相关系数，其值在[-1, 1]范围内，d_i^2表示第i对数据的排序差值的平方，n为样本量。正的ρ表示ESs之间的协同关系，负的ρ表示权衡关系。如果ρ不显著或为0，则表示两个ESs之间无相关性。

GWR模型的计算公式如下：

$$
Y_i = β_0(u_i, v_i) + \sum_{j=1}^p β_j(u_i, v_i) X_{ij} + ε_i
$$

其中，Y_i是第i个位置的因变量值，X_ij（j = 0, 1, …, p）是第i个位置的自变量值，p是自变量的数量，(u_i, v_i)表示第i个位置的坐标，β_0(u_i, v_i)是截距项，β_j(u_i, v_i)（j = 0, 1, …, p）是回归系数，ε_i表示随机误差。

#### 2.3.4. 驱动机制分析

RF和PLS-SEM在方法论上实现了互补整合：RF提供了客观的变量重要性排序和筛选标准，而PLS-SEM揭示了驱动因素对ESs及其权衡/协同关系的直接和间接路径、影响方向和作用程度。

### 3. 结果分析

#### 3.1. 生态系统服务的时空变化

研究结果显示，LYRR中ESs的时空变化存在显著差异。1990年，高WY区域主要分布在东部，低WY区域集中在西南部，整体呈现从西南向东部逐渐增加的趋势（图4a）。与1990年相比，2000年的WY空间分布特征表现出相反的趋势：西南部WY较高，而北部和东部WY较低，整体呈现从西南向东北逐渐减少的趋势（图4b）。2000年，水产出服务的空间分布出现了更明显的转变，这可能与降水作为直接因素的影响有关。2010年，低值区域主要分布在西南部和水域，而高值区域则没有明显的聚集，表现出不均匀分布，整体上东部的WY高于西部（图4c）。2020年，LYRR的总体WY呈现出“东高西低”的空间格局，高值区域主要分布在东部和东南部，而低值区域的分布范围较广，主要集中在北部、西部和西南部（图4d）。从坡度趋势（图4e）来看，1990年至2020年间，LYRR中部的WY出现了显著下降，而东部和南部的小面积区域则表现出轻微增长。内湖以及北部和东南部的水域WY基本保持稳定，而大部分区域则呈现出轻微下降的趋势。

#### 3.2. 生态系统服务之间的相互作用

从1990年至2020年，ESs之间的TOS表现出显著的时空异质性，主要表现为弱权衡和弱协同（图9）。总体来看，除了CS-FP以外，其他ESs之间均以弱协同为主。WY-CS的权衡关系在2010年达到最强（-0.42***），而CS-FP的权衡关系则在2020年转变为协同关系（-0.21），与Liu et al. (2024a)和Yang et al. (2021b)的研究结果一致。然而，Spearman与GWR分析之间出现了显著差异，揭示了明显的尺度效应，这与Zhang et al. (2022)的研究结论不同。这种差异主要源于局部空间异质性：Spearman作为全局统计方法，难以捕捉空间非平稳性，而GWR则通过局部加权回归有效地揭示了TOS关系的空间异质性。结合这两种方法，能够提供全局和局部尺度的互补分析视角。空间相关性通常揭示了调节和维持服务之间的协同效应，而供给服务往往与其它服务之间存在权衡（Yang et al., 2021b；Zhang et al., 2022）。研究结果表明，除了CS-FP以外，其他ESs之间的相互作用主要表现为弱权衡。这种现象与供给服务与其他服务之间的权衡模式相吻合。然而，CS-FP之间的弱协同关系和CS-SC之间的弱权衡关系，与调节和支持服务之间通常存在的协同效应模式相矛盾。这种差异可能源于LYRR独特的生态系统结构、功能和复杂的动态过程，同时也反映了全局统计与局部建模在识别空间异质性方面的固有差异。

#### 3.3. 驱动机制分析

RF分析结果显示，1990年至2020年间，各驱动因素对不同ESs的影响程度存在显著差异。例如，LUCC、Pre、DEM和PET是WY的主要驱动因素；LUCC和SLP是CS的主要驱动因素；DEM是SC的主要驱动因素；而LUCC、NDVI、PET和POP则显著影响FP。此外，不同驱动因素对ESs的影响方向和强度在不同年份表现出显著变化。例如，在1990年，Pre（0.734***)对WY产生了最强的正向影响，而LUCC（-0.934***）对CS的影响则最显著。然而，它们通过中介路径产生的间接效应仍较弱。这种变化趋势表明，驱动因素对ESs的影响不仅取决于其直接作用，还受到区域生态结构和功能变化的影响。

### 4. 讨论

#### 4.1. 生态系统服务的时空动态及其驱动因素

研究发现，从1990年至2020年，LYRR中WY的空间异质性最为显著，整体呈现下降趋势。Pre和LUCC是影响WY的两个关键驱动因素。以往研究表明，该地区降水在过去70年中持续减少（Wang et al., 2022），这表明气候因素为WY下降提供了重要的自然背景。尽管城市化增加了不透水地表，从而在一定程度上提高了城市区域的径流能力（Li et al., 2020），但随着城乡用水需求的持续增长以及农业灌溉用水的大幅增加，区域总水资源消耗量稳步上升。同时，植被恢复政策如“退耕还林”项目（Yang et al., 2021b）显著提高了生态用地的植被覆盖率，增强了冠层截留和植被蒸腾，从而减少了地表径流的生成。总体而言，城市化带来的径流增加远小于生态恢复和气候条件对径流减少的影响，导致整个区域的WY能力下降。

在1990年至2020年间，LYRR中CS和SC的空间格局表现出相似性。高CS区域主要分布在东部和西南部的山地和丘陵地区，这与Yang et al. (2021b)的结论一致。驱动机制分析表明，LUCC和DEM分别是CS和SC的主要驱动因素，部分与Zhang et al. (2021)的结论相吻合。CS的增加主要源于森林和草地的广泛分布。高植被覆盖率通过光合作用促进碳汇，并通过植物残体的长期存储增强土壤中的有机碳储量（Niu et al., 2022）。SC能力主要受DEM控制，植被冠层截留降水和根系土壤稳定化有效缓解了侵蚀，从而增强了SC功能。而低SC区域则主要集中在地形平坦的平原地区，这些区域由于植被覆盖率低、土壤有机碳储存能力弱以及人类活动带来的高碳排放强度，导致SC能力较低。

FP的空间分布模式在1990年至2020年间保持稳定。高FP区域主要集中在平原地区的耕地，而低FP区域则主要分布在水域、建设用地以及东部和西南部的高地形山地地区。这种分布模式表明FP高度依赖于高质量的耕地资源，并受到城市化和地形条件的限制。这种空间差异主要由几个关键因素驱动：LUCC、NDVI和POP。具体而言，耕地面积的变化是影响粮食供应稳定性的关键因素。NDVI作为植被生长的指标，能够有效反映作物生长和粮食产量的潜在变化。人口增长直接增加了粮食需求，高人口密度进一步提升了整体粮食需求。因此，平原地区受益于集中的耕地和有利的生产条件，粮食供应能力较强。相比之下，水域、建设用地和山地地区则因自然条件和土地利用模式受到限制，粮食生产能力普遍较低。

#### 4.2. 权衡与协同关系的尺度效应

本研究评估的TOS关系在一定程度上与相关研究（Yang et al., 2021b）存在相似性和差异性。总体相关性显示，WY-SC、WY-FP和CS-SC之间存在协同关系，而WY-CS和SC-FP之间则表现出权衡关系。CS-FP关系从权衡转变为协同，这与Liu et al. (2024a)和Yang et al. (2021b)的研究结果一致。然而，Spearman与GWR分析之间存在显著差异，揭示了明显的尺度效应，这与Zhang et al. (2022)的研究结论不同。这种差异主要源于局部空间异质性：Spearman作为全局统计方法，难以捕捉空间非平稳性，而GWR则通过局部加权回归有效地揭示了TOS关系的空间异质性。结合这两种方法，能够提供全局和局部尺度的互补分析视角。空间相关性通常揭示了调节和维持服务之间的协同效应，而供给服务往往与其它服务之间存在权衡（Yang et al., 2021b；Zhang et al., 2022）。研究结果表明，除了CS-FP以外，其他ESs之间的相互作用主要表现为弱权衡。这种现象与供给服务与其他服务之间的权衡模式相吻合。然而，CS-FP之间的弱协同关系和CS-SC之间的弱权衡关系，与调节和支持服务之间通常存在的协同效应模式相矛盾。这种差异可能源于LYRR独特的生态系统结构、功能和复杂的动态过程，同时也反映了全局统计与局部建模在识别空间异质性方面的固有差异。

#### 4.3. 局限性与未来展望

尽管本研究在量化评估ESs方面力求严谨，但在数据收集和参数确定过程中仍不可避免地存在误差。例如，在使用InVEST模型估算CS时，碳库的碳密度值参考了Yang et al. (2021a)的研究，考虑了不同土地利用类型下的降水量和年平均温度。同样，RUSLE模型中使用的土壤侵蚀防治因子P是基于相关研究的综合评估（Wang et al., 2024b；Xiao et al., 2020；Zhang et al., 2024）。此外，RUSLE模型在建设用地的应用中存在固有的不确定性。本研究中使用的C因子主要基于植被覆盖率，这种方法在自然和农业用地中被证明是有效的。然而，建设用地的侵蚀主要发生在建设阶段，此时地表受到严重扰动，土壤松散度高。建设完成后，硬化地表几乎完全阻止侵蚀。由于土地利用数据无法区分建设阶段和稳定阶段，本研究中使用的年度平均C因子值可能高估了稳定建设用地的侵蚀量，因此可能对估算结果产生一定影响。其次，受限于时空数据的分辨率，本研究仅依赖1990年至2020年的多时期数据，这在时间跨度上显得不足。未来的研究应结合更高分辨率的遥感产品与长期地面观测数据，以分析更长时间跨度和更精细空间与时间尺度的趋势。

最后，当前FP的计算仅考虑了耕地中的粮食产量。未来的研究可以纳入更多粮食来源，如林地、草地和水域中的水果、蔬菜、肉类、蛋类、奶制品和水产品。此外，在使用PLS-SEM时，研究选择了WY-CS、CS-FP和SC-FP作为三个观测变量，以代表整体TOS的潜在变量，这可能削弱了驱动因素通过ESs中介路径产生的间接效应。本研究主要依赖单因素分析；未来的研究可以进一步探讨多因素的相互作用，以增强对ESs动态的理解。同时，随着气候变化和人类活动的持续演变，未来的研究可以将情景预测与ESs变化相结合，以评估在不同条件下的适应性和可持续性。

#### 4.4. 政策建议

基于对ESs时空动态、TOS和驱动机制的全面分析，本研究提出了以下具体的政策建议：

（1）基于TOS的时空格局实施生态红线控制。在权衡区域（如WY-CS），应采取水-碳协同的植被恢复策略，与双碳目标相协调。在低WY区域，应避免高耗水的造林项目，优先选择耗水低、固碳能力强的本地物种。在高CS区域，应避免盲目造林，以减少水压，确保生态恢复与水源保护之间的平衡。在协同区域（如CS-FP），应优先纳入生态红线，推广生态农业，采用混合种植模式，以提升粮食生产与碳汇的协同效应。

（2）优化以LUCC调控为中心的土地利用分区。首先，严格保护东部泰山和西南部太行山麓的山地生态源区（高CS和SC区域），限制开发和耕作，保持其水源保护、土壤保持和碳汇的核心功能。其次，优化中部平原的主产区（高CS区域），通过推广保护性耕作和生态沟渠，减少非点源污染和土壤侵蚀。同时，建立农田林网，增强农业生态系统中的碳汇，实现“藏粮于地”和“藏碳于田”的协同效应。此外，应通过设立城市扩展边界，严格保护高产农田和城市周边的关键生态节点，对城市扩张实施严格限制。

（3）建立基于空间显式评估的动态监测与预警机制。构建一个整合多源数据的区域ESs监测与预警系统，定期评估服务之间的平衡。根据监测结果，动态调整生态修复与保护项目，实施“一区一策”的精准管理，提高治理能力以适应变化。

（4）基于定量ESs评估，完善生态补偿政策。在关键生态功能区，应建立与ESs数量和价值变化相匹配的精准补偿标准，鼓励地方政府和农民参与生态修复与保护项目，以实现互利共赢。研究发现，LYRR中CS-FP的关系从权衡转变为协同，因此建议推进跨区域和多功能的生态补偿政策，以维持这种协同发展的趋势。