生态系统服务是人类从生态系统中所获得的各种惠益^[1], 它是通过生态系统的功能直接或间接得到的产品和服务^[2], 是连接自然环境与人类福祉的桥梁^[3]。不同的生态系统服务之间存在密切而复杂的关系, 主要表现为此消彼长的权衡关系和相互增益的协同关系^[4-5], 加强对生态系统服务之间的权衡与协同关系的理解, 有利于提高生态系统总体效益、保证生态系统服务的可持续供给^[6]。

生态系统服务的时空变化以及不同生态系统服务之间的关系是生态系统服务研究的重要内容, 国内外学者已对此做了很多研究。Yang等^[7]利用玫瑰图以及生产可能性边界(PPF)的方法, 对关中-天水地区的碳固存、产水和土壤保持之间的权衡与协同效应进行了评价, 结果表明碳固存与土壤保持之间是显著权衡关系, 碳固存与产水量以及产水量与土壤保持之间是协同关系。Lautenbach等^[8]比较分析了3种(地图对比法、情景分析法、基于优化景观权衡分析法)生态系统服务权衡的方法, 结果指出地图对比法和情景分析法是目前研究生态系统服务权衡关系的有效方法。Howe等^[9]通过分析给人类带来福祉的全球生态系统服务的案例发现权衡记录几乎是协同增效的3倍, 权衡比协同更能带来理想结果。傅伯杰等^[10]以黄土高原地区为例探索了生态系统服务权衡及区域集成方法, 研究结果表明, 土地利用变化与土壤保持、碳固定具有正效应, 与产水量之间存在负效应。陈登帅等^[11]利用CA-markov模型预测了渭河流域2050年在最优生态系统服务下的土地利用格局。王鹏涛等^[12]基于逐像元偏相关的时空统计制图方法对汉江流域的土壤保持服务、产水服务、NPP进行空间制图, 并分析了三者之间的权衡与协同关系。但是目前对生态系统服务权衡与协同关系的研究多集中在局地尺度, 对宏观尺度的研究较少^[3]。

生态恢复是通过人为干预、启动或促进退化生态系统恢复进程的活动, 被认为是应对生态退化, 改善生态系统服务的有效手段^[13]。生态恢复能够从物种组成、生态系统功能、景观环境和生态系统稳定性4个方面对生态系统产生作用, 改变生态系统的结构、功能及自我维持能力从而对生态系统服务产生重要影响。近年来, 对于生态恢复的研究也逐渐聚焦到生态系统服务上来^[14]。已有研究对全球范围内89项生态恢复案例进行meta分析, 结果显示生态恢复与生态系统服务呈正相关关系, 能够促使生态系统服务恢复25%, 但尚未达到退化前水平。大量研究表明，造林对生态系统服务有很大影响, 它可以保持土壤且增加碳汇, 但是会使产水量降低^[15-29]。Jia等^[17]研究表明在退耕还林工程的影响下供给服务和调节服务之间表现出权衡关系, 调节服务之间表现出协同关系。Koch等^[18]分析证明了西澳洲实施的红柳桉树林恢复项目使得该区域的碳固存和水源涵养得到明显改善。Marton等^[19]证实了美国湿地保护项目和保育休耕项目恢复了湿地和滨岸缓冲带由农田开发而丧失的水质净化服务。孙艺杰等^[20]探究了退耕还林政策和地形坡度对生态系统服务权衡与协同关系动态变化的影响。虽然目前已有很多研究来证明生态恢复对生态系统服务的改善作用, 但是目前还没有系统地评价恢复行动在提高生态系统服务方面的效力^[21]。过去50年来, 我国在国家和地方层面上开展了大量的生态恢复技术研发工作以及大量的生态工程建设, 这些工程有的已经显现出一定的效益, 有的还处在摸索和调整阶段。由于生态恢复是一个长期的动态过程, 对其进行科学、系统的评价能够为生态恢复提供反馈信息, 进而为改善生态系统提高生态系统服务功能提供有益的信息参考。因此探究生态恢复对生态系统服务权衡与协同关系的影响对进一步调整生态恢复工程具有重要意义。

我国是世界上生态环境脆弱区分布面积最大、脆弱生态类型最多、生态脆弱性表现最明显的国家之一。据统计, 中度以上生态脆弱区面积约占陆地总面积的55%^[22], 并且大多位于生态过渡区和植被交错区, 是典型的农牧、林牧、农林等复合交错带, 也是我国目前生态问题突出、经济相对落后和人民生活贫困区, 同时也是我国环境监管的薄弱地区。目前对我国重点脆弱生态区的生态脆弱性已有学者做过不同程度的研究。郭兵等^[23-24]分别对青藏高原区和西北干旱荒漠区的生态系统脆弱性做了研究, 结果表明青藏高原区的生态系统脆弱性整体处于中度脆弱状态, 西北干旱区则处于中-重度脆弱状态, 并且其生态系统脆弱性时空变化分异格局与地形、气候(气温、降水等)、人口密度存在显著相关性。彭建等^[25]以云南省永胜县为例用生态系统服务功能衡量了滇西北生态脆弱区土地利用类型的相对生态价值。高旺盛等^[26]通过分析黄土高原生态系统服务价值发现生态环境脆弱区的生态系统服务价值依然存在, 但远不足以维持系统的自我循环, 才会导致生态环境恶化趋势日趋严重。目前, 多数研究主要集中在生态脆弱区的脆弱性实证评价研究^[27], 对整个重点脆弱区的生态系统服务权衡与协同关系研究较少。生态脆弱区生态恢复是保证我国及生态脆弱区经济可持续发展的必然选择。因此, 探析我国重点脆弱生态区的生态系统服务的时空变化格局以及不同生态系统服务之间的权衡与协同关系, 有利于科学制定与实施我国重点脆弱生态区的生态系统管理政策^[21], 对与实现生态系统服务效益最大化^[28]和保障国家生态安全的重大战略需求^[29]具有重要的意义。

基于此, 本文采用RUSLE(The Revised Universal Soil Loss Equation)模型、CASA(Carnegie Ames-Stanford Approach)模型和InVEST产水模型分别评估我国重点脆弱生态区1990—2015年土壤保持、植被净初级生产力(Net primary production, NPP)和产水服务的时空分布格局特征, 采用皮尔逊积矩(Pearson)相关系数法量化不同服务空间上的权衡与协同关系, 并基于植被分布场数据(Vegetation continuous fields, VCF)分析了生态恢复对生态系统服务权衡与协同关系动态变化的影响, 以期为大尺度的国家重点脆弱生态区生态系统服务优化提供参考, 为实现人类福祉和人类社会与生态系统“双赢”的目标而努力^[6]。

1 研究区概况与研究方法 1.1 研究区概况

我国重点脆弱生态区包括北方风沙区、西北干旱区、黄土高原区、青藏高原区和西南喀斯特区(图 1)。其中北方风沙区位于我国北部的农牧交错带上, 总面积约68万km², 以半干旱的草原生态系统为主。西北干旱区分布于贺兰山以西, 祁连山、阿尔金山、昆仑山以北和新疆全境的广大区域, 面积约176万km², 以山地和荒漠生态系统为主。黄土高原区位于我国中部偏北部, 面积约64万km², 高原内大部分为厚层黄土覆盖, 经流水长期强烈侵蚀, 逐渐形成千沟万壑、地形支离破碎的特殊自然景观, 严重的水土流失是该区域主要生态环境问题。青藏高原区位于我国地形的第一阶梯, 平均海拔4000 m以上, 面积约258万km², 是我国最大的生态脆弱区；该区生态系统结构简单、抗干扰能力弱并且易受全球环境变化的影响。西南喀斯特区是我国喀斯特地貌的主要分布区, 面积约193万km², 该区生态环境脆弱易受人类活动影响。

1.2 数据来源

气象数据包括降水、气温、日照时数等来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn)；遥感影像数据来源于地理空间数据云Landsat TM、MODIS数据；DEM数据亦来源于地理空间数据云ASTER GDEM V2版数据, 分辨率为30 m；土壤数据采用基于世界土壤数据库(HWSD)的中国土壤数据集(v1.1)(http://westdc.westgis.ac.cn)；植被分布场数据来自Song等人的研究成果^[30](https://doi.org/10.1038/s41586-018-0411-9)。该数据是将地表表示为植被功能类型的分类组合, 所采用的数据集主要来源于高级超高分辨率辐射计(AVHRR), 该数据集组合了来自多个卫星传感器的光学观测, 如中分辨率成像光谱仪、Landsat增强型专题成像仪等。

1.3 研究方法 1.3.1 土壤保持服务

本文的土壤保持量是以TM、MODIS等遥感数据为基础, 结合气象、DEM、植被和土壤类型数据, 生成模型中涉及的各项侵蚀因子、坡长坡度因子和人为管理因子, 基于RUSLE模型计算单个像元内的土壤保持量^[31]。该模型适用范围广, 可操作性强。

式中, SC为土壤侵蚀指数(t hm^-2 a^-1), R为降雨侵蚀力因子, K为土壤可蚀性因子, LS为坡长坡度因子, C为植被覆盖因子, P为水土保持措施因子。

1.3.2 植被净初级生产力

本文主要基于Potter等建立的Carnegie Ames-Stanford Approach(CASA)模型, 来计算重点脆弱生态区的植被净初级生产力。其估算方法如下^[32]：

式中, NPP为植被净初级生产力(gC m^-2 a^-1), APAR (x, t)为像元x在t月时, 植被吸收的光合有效福射(MJ/m²), ε(x, t)表示像元x在t月时, 植被实际的光能利用率(gC/MJ)。

1.3.3 产水服务

本文生态系统产水服务, 采用区域水量平衡法进行评价, 利用区域水分的输入量降水与输出量蒸散发的差值, 得到区域生态系统的产水量^[33]。该方法在全球不同气候类型的流域内进行了校正^[34], 所得结果较为科学, 计算公式如下：

式中, Q为产水量(mm), P为年降水量(mm), ET为蒸散量(mm), PET为潜在蒸散量(mm), 可采用Hamon潜在蒸散法求得^[35]。ω为植被耗水系数, 根据前人已有校正结果来设置^[36-38]。

1.3.4 生态恢复区的识别确定

植被分布场数据相比较于基于分类系统的土地利用数据能够更加具体量化地表的变化, 它表示每个陆地像素被TC(tree canopy), SV(short vegetation)和BG(bare ground)覆盖的百分比, 以此代表当地高峰生长季节时的植被组成。其中, TC指的是树冠垂直投影覆盖的地面比例。SV表征除树木以外的植被覆盖的地面比例, 包括灌木, 草本植被和苔藓, 而BG表示未被植被覆盖的地表比例。树木被定义为高度超过5米的所有植被^[30]。本文基于Matlab R2016a平台, 采用Theil-Sen斜率(TS_slope)估计^[39-40]和Mann-Kennal趋势检验^[41]的方法获取重点脆弱生态区植被分布场的趋势变化情况, 如果M-K检验无统计学意义(P>=0.05), 净变化为0, 如果趋势检验显著(P < 0.05), 则用趋势分析估计每个像元的净变化。在人工种植植被改善生态环境的背景下将TC、SV、BG的变化组合进行分类(表 1), 将显著变化的像元分类(P < 0.05), 提取出表 1中前4种类型为植被连续变化显著的区域, 即生态恢复区。其他变化不明显、没有变化以及植被覆盖度减少的区域为非生态恢复区。其中生态恢复区的面积约为214万km², 非生态恢复区的面积约536万km²(图 2)。

1.3.5 生态系统服务时空变化及相关性量化方法

本文在分析1990—2015年重点脆弱生态区3种服务变化趋势的基础上, 采用基于Matlab R2016a平台的皮尔逊积矩(Pearson)相关系数法对生态系统服务权衡与协同关系进行量化。并利用可视化制图来量化其空间差异。利用t检验的方法, 将生态系统服务之间的权衡与协同关系分为7类, 即显著协同(r>0, P≤0.05)；较显著协同(r>0, 0.05 < P≤0.1)；协同(r>0, P>0.1)；显著权衡(r < 0, P≤0.05)；较显著权衡(r < 0, 0.05 < P≤0.1)；权衡(r < 0, P>0.1)；其余的为无相关。

2 结果与分析 2.1 生态系统服务时空变化分析

1990—2015年重点脆弱生态区的生态系统服务整体上呈上升趋势, 其中NPP上升趋势最为明显, 产水的年际变化波动较大(图 3：(a)、(b)、(c))。具体而言, 1990—2015年重点脆弱生态区土壤保持服务的年平均值为1426.24 t km^-2 a^-1, 最高值在1998年(1932.60 t km^-2 a^-1), 最低值在1992年(1157.60 t km^-2 a^-1)；NPP的年均值为376.43 gC m^-2 a^-1, 最高值在2013年(419.37 gC m^-2 a^-1), 最低值在2000年(352.71 gC m^-2 a^-1), 多年来平均以2.04 gC m^-2 a^-1的速率增长。产水服务的年均值为225.97 mm a^-1, 年际变化波动较大, 最高值在2015年(301.58 mm a^-1), 最低值在1992年(172.19 mm a^-1)。从生态系统服务在时间上的变化关系可以看出, 产水量与土壤保持在1990—2015年呈现出比较明显的协同变化规律。NPP与其他两种服务在1990—2000年间均呈现出权衡变化规律。

从图 3((d)、(e)、(f)))可以看出, 土壤保持量在各个子区域中的年际变化较缓和, 西南喀斯特区的土壤保持量明显高于其他区域, 也是年际变化量最大的区域, 其他从高到低依次为青藏高原区、黄土高原区、北方风沙区和西北干旱区。NPP在各个子区域中的年际变化都不显著, 西南喀斯特的NPP服务量明显高于其他区域, 其次是黄土高原区, 最低的是西北干旱区, 北方风沙区和青藏高原区不相上下, 但北方风沙区的年际变化大于青藏高原区, 这主要是因为北方风沙区受人类活动的影响较大。对于产水服务而言, 各区域年际波动变化较其他两种服务大, 其中变化较大的是西南喀斯特区和黄土高原区, 西北干旱区的均为负值且变化不大, 北方风沙区的值在-100~100之间波动且小于0的年份要多于大于0的年份。

在空间特征上(图 4), 3种服务都呈现出南多北少的特征, 变化趋势在不同区域表现不同, 图中的变化量均通过显著性检验(P < 0.05)。土壤保持服务各个子区差异较小, 相对高值区位于青藏高原区的东南部和西南喀斯特区；土壤保持量明显增加的地区几乎没有, 减少较明显的地区除了青藏高原西北部连片分布以外, 其他分布较散。NPP的值明显呈现出南多北少分布特征, 高值区和低值区集中连片的分布；黄土高原区和喜马拉雅山脉东南部的小范围地区是变化量的高值区, 大部分地区基本不变。产水服务的空间差异相对较小, 西北干旱区、和北方风沙区是相对低值区, 高值区分布在青藏高原区南部和西南喀斯特区的东部；整体上以基本不变为主, 主要变化的区域在青藏高原区和西南喀斯特区。这种分布特征与我国水热分布和地理环境的差异有很大关系。

2.2 生态系统服务权衡与协同时空动态变化

本文对整个重点脆弱生态区1990—2015年两两服务间的相关系数进行了动态趋势分析, 图 5的相关系数均经过显著性检验(P < 0.05)。从整个重点脆弱生态区的年平均水平来看生态系统服务之间的相关系数均大于0, 说明整个重点脆弱区生态系统服务以协同为主, 并且协同关系均呈减弱的趋势。其中协同性降低速度最快的是产水与NPP(平均下降速率为5.7%)。变化较慢的是产水与土壤保持, 并且产水与NPP整体上的协同性较其他两者较高, 平均相关系数为0.32。总之, 26年来重点脆弱生态区的生态系统服务整体表现为协同减弱的趋势。

在空间上(图 6), 就土壤保持与产水而言, 整体上以显著协同为主, 西北干旱区部分地区出现权衡关系；对于土壤保持与NPP, 西北干旱区的大部分和青藏高原的北部以及北方风沙区的大部分地区为无相关关系, 在青藏高原南部和西南喀斯特区较显著权衡和权衡关系交叉分布, 黄土高原区以协同关系为主；NPP与产水在西北干旱区的大部分和青藏高原的北部以及北方风沙区的大部分地区为无相关关系, 在黄土高原地区以协同关系为主, 而西南喀斯特区以显著权衡以及权衡关系为主。

就各个子区域的具体差异而言, 在土壤保持与产水方面, 各子区域的显著协同关系占该区总面积的比例都在65%以上(图 7), 其中黄土高原最高(86.04%), 三种等级的权衡关系在各区的所占比例都非常小(1%左右)。在土壤保持与NPP方面, 各子区中无相关关系所占比例较大, 最高的是西北干旱区(83.01%), 在黄土高原区域内协同关系的分布面积占50.54%。这与黄土高原极易被侵蚀的土壤类型有关, 该区的土壤保持和NPP表现出比其他区域较明显的相互增益的协同关系。此外, 西南喀斯特区的权衡关系也较明显, 占该区的47.01%。在NPP与产水方面, 除了无相关关系分布较广以外, 黄土高原的协同关系较突出(占38.48%), 西南喀斯特的显著权衡关系和权衡关系较明显, 分别占该区的30.84%和46.60%。总之, 5个子区整体的生态系统服务权衡与协同关系与整个重点脆弱生态区的一致, 都是以协同关系为主。

2.3 生态恢复对生态系统服务权衡与协同的影响

植被分布场数据的显著变化归因于直接的人类活动和自然驱动因素, Song等人研究结果表明中国北方区域VCF的变化中受人类活动直接影响的超过了75%^[30]。为了维护国家生态安全, 我国采取不同的生态工程来改善一些脆弱的生态环境, 比如退耕还林还草工程、三北防护林工程和天然林保护工程等等, 这些一系列工程均影响我国VCF的变化, 从而对区域生态系统服务功能产生一定的影响。因此, 本文在大规模人工种植改变植被分布的背景下来提取生态恢复区和非生态恢复区各类生态系统服务, 分析其权衡与协同关系的动态变化(图 8), 探讨生态恢复对生态系统服务权衡与协同关系变化的影响。

在生态恢复过程中, 生态系统服务之间的协同是指不同生态系统服务的变化趋势相同, 共同增长或减弱；权衡是指不同服务之间的变化趋势相异^[13]。2000年前, 恢复区的产水与NPP为协同关系(相关系数为0.21), 而非恢复区为权衡关系(相关系数为-0.03), 2000年之后恢复区和非恢复区的权衡与协同关系差异减小。土壤保持与NPP的权衡与协同关系在恢复区和非恢复区的差异最大, 恢复区的土壤保持与NPP的协同程度(相关系数为0.25)要远大于非恢复区(相关系数为-0.004), 这表明生态恢复工程增强了土壤保持与NPP之间的协同关系, 伴随着NPP的增加土壤保持也不断增加, 彼此之间相互促进。土壤保持与产水在生态恢复区与非生态恢复区均为协同关系, 但两者的差异较小(相关系数均约为0.2)。整体上恢复区的生态系统服务的协同程度要大于非生态恢复区, 且协同关系总体呈减弱的趋势, 但因为非恢复区的生态系统比较稳定, 而恢复区的生态系统处于恢复阶段, 自我调节能力和稳定性较差, 难以在短时间内形成稳定且完整的生态系统^[42]。所以在相对较短的时间范围内恢复区的生态系统服务的协同程度与非恢复区的并没有太大的差异。

就各个子区而言(图 9), 每个子区的土壤保持与NPP无论在恢复区还是非恢复区都是协同关系, 除了西北干旱区内恢复区的协同程度低于非恢复区外, 其他子区内恢复区的土壤保持与NPP的协同程度均高于非恢复区。这说明了生态恢复工程增强了重点脆弱生态区的土壤保持和NPP的协同关系。对于产水与NPP而言, 除了北方风沙区基本为无相关以外, 其他子区恢复区内均为权衡关系, 并且权衡关系均成减弱趋势, 说明生态恢复减弱了各子区产水与NPP的权衡关系, 北方风沙区、青藏高原区、西北干旱区的土壤保持与产水在恢复区的协同程度均高于非恢复区；黄土高原区在恢复区和非恢复区的土壤保持与产水的相关系数均在±0.01范围内, 但恢复区在2000年以后有向协同转变的趋势；西南喀斯特区的土壤保持与产水表现为权衡关系, 从动态变化来看, 恢复区权衡程度减弱的趋势明显高于非恢复区, 这表明生态恢复减弱了该区土壤保持与产水的冲突关系。总之, 从整体上看, 生态恢复区生态系统服务的协同程度大于非生态恢复区, 权衡关系减弱的趋势高于非生态恢复区, 即生态恢复对于生态系统加强协同减弱权衡有一定的促进作用, 这符合生态系统服务加强协同、减弱权衡的生态管理目标^[43]。但有些区域仍需进一步调整生态工程加强对生态系统的管理。

3 结论与讨论 3.1 讨论

生态系统具有尺度效应^[44], 不同等级时空尺度上的生态系统服务的效益差异较大^[45-46]。目前对于生态系统服务权衡与协同关系研究多集中在局地尺度上的时空动态分析及空间表达上, 对大尺度上生态系统服务及其权衡与协同关系在长时间序列上的时空变化的分析较少^[3]。本文定量测算了1990—2015年我国重点脆弱生态区的土壤保持、产水、NPP 3种服务, 分析了这3种生态系统服务间的权衡与协同关系及其生态恢复工程对其造成的影响。结果表明我国的生态恢复工程对于改善生态环境起到了一定的作用, 但个别区域效果依然不理想, 例如西南喀斯特区, 生态恢复工程还需调整, 以此来保障生态恢复的可持续发展。尽管生态恢复工程通过改变生态系统格局对生态系统服务有一定的影响, 但生态系统服务还受其他多种复杂因素的影响, 例如全球变化^[47]人类社会相应政策^[48]、人口、城市化和工业化^[49]等因素的影响。MA千年生态系统评估将生态系统服务划分为四大类, 包括支持服务, 供给服务, 调节服务以及文化服务^[1]。本文仅涉及到供给、调节和支持服务, 并不能完全代表重点脆弱生态区的全部生态系统服务, 比如北方风沙区以防风固沙服务为主, 可以减缓风沙对京津冀地区造成的危害。而且从本文的分析结果可以看出各子区域上的生态系统服务权衡与协同关系与整个重点脆弱生态区上的并不一致, 例如整个重点脆弱生态区的NPP和产水整体上表现为协同关系, 但在西南喀斯特区却是权衡关系为主。这主要是因为各个子区的自然环境条件差异较大, 所以在宏观研究的基础上还需要探讨不同子区各种生态系统服务变化的深层次原因。因此在多尺度上进行多种类型生态系统服务权衡与协同关系研究并深入探究影响生态系统服务权衡与协同关系的其他因素将成为未来研究的重要方面。

3.2 结论

(1) 1990—2015年间, 我国重点脆弱生态区的土壤保持、NPP及产水呈现出南多北少的格局, 并整体上有增加的趋势。3种生态系统服务之间总体上呈协同关系, 且协同关系呈减弱的趋势, 在各子区域上, 土壤保持、NPP及产水三者之间以协同关系为主。但不同区域之间存在着空间异质性。我国重点脆弱生态区的生态系统服务在进行生态系统管理时, 需要同时考虑到整体与局部空间的不同, 以制定更为合理的生态保护政策。

(2) 恢复区和非恢复区的生态系统服务之间总体上表现出协同关系, 并且恢复区的相关系数总体上高于非恢复区, 说明生态恢复工程对生态系统服务有一定的影响。不过生态恢复工程对生态系统服务的影响是一个长期的过程, 在制定生态恢复工程时, 应考虑长远的可持续发展前景。

(3) 各子区域的统计结果表明, 生态恢复区生态系统服务的协同程度大于非生态恢复区, 动态变化上权衡关系减弱的趋势高于非生态恢复区。但各个区域内又存在差异, 所以在制订或完善生态恢复工程时应在把握整体特征的基础上结合各子区的特征, 让生态恢复工程更好的推动生态环境的可持续发展。