杨 阳,窦艳星,王云强,安韶山

1 中国科学院地球环境研究所 黄土与第四纪地质国家重点实验室,西安 710061 2 中国科学院第四纪科学与全球变化卓越创新中心,西安 710061 3 陕西黄土高原地球关键带国家野外科学观测研究站,西安 710061 4 西北农林科技大学林学院, 杨凌 712100 5 西北农林科技大学黄土高原土壤侵蚀与旱地农业国家重点实验室, 杨凌 712100

生态系统服务是人类从生态系统中直接/间接地获得的收益,它侧重于人类从大自然中所获取的各项服务,影响着生态环境质量和生态系统的可持续性发展；开展生态系统服务的研究能够帮助我们有效促进生态环境质量的改善,对于资源短缺等问题提供很好的解决办法[1—3]。近百年来,随着人类对大自然开发利用程度的加重,加上自然环境的恶化,全球范围内的生态系统结构和功能发生了巨大的变化,同时伴随着空气质量下降、水污染等一系列的环境问题,各项生态系统服务出现了急剧降低,这给人类的生产生活带来了极大的困扰,已成为全球性问题[4—5]。21世纪初,在联合国的组织下,“生态系统服务与人类福祉”计划(MA)得以实施,该计划表明,近60%以上的生态系统服务出现了不同程度的下降,并且这种下降趋势很可能会持续,这严重阻碍了人类的可持续发展[6—8]。据相关统计,生态系统功能能够为人类创造高达15万亿英镑的价值[9]；然而在人类不合理的开发利用下,这种价值损失了近2/3,诸如水土流失、生物多样性大幅下降等问题,在导致生态功能下降的同时,也破坏了生态环境的可持续发展[7—8],如何评估和改善生态系统服务显得刻不容缓[10—11]。

对于生态系统服务的评估方法,出现了一些以复杂理论和研究成果为基础的评估模型,如InVEST(Integrated Valuation of Ecosystem Services and Tradeoffs)、ARIES(Artificial Intelligence for Ecosystem Services)和SolVES(Social Values for Ecosystem Services)等模型[12]。对比而言,InVEST模型发展更为成熟,且处于不断更新和完善阶段,具有“空间化、精细化、定量化”的特征,因此得到广泛应用,并在全球范围内均取得较好的反响[13—16]。然而,由于生态系统服务种类多样、空间分布不均衡、以及人类的选择性使用,使得各种服务之间呈现出相互作用、联系及交织的动态关系,具体表现为协同和权衡关系。当前,已有较多的研究探讨了多种尺度(如国家、区域、流域等尺度)上生态系统服务之间的关系,这些生态系统服务的协同和权衡关系不尽一致[13—15]。例如,陈登帅等[17]采用生产可能性边界法定量化了渭河流域生态系统服务之间的关系,发现生物多样性与产水量呈权衡关系,而固碳与生物多样性呈协同关系；对于延安市生态系统服务,孙艺杰等[18]的结果表明食物供给与土壤保持表现为权衡关系,净初级生产力(Net primary productivity,NPP)与土壤保持、水源涵养表现为协同关系[19]。张瑜[20]等利用InVEST模型评估了黄土高原生态系统服务价值的动态变化特征。结合InVEST模型与情景分析方法,Ouyang等[21]对多种生态系统服务在不同土地利用方式下的相互作用关系,及流域土地利用变化趋势进行了定量研究,结果表明,科学地使用土地能够有效权衡各生态系统服务直接的关系,促进区域的可持续发展。傅伯杰等[22]总结生态系统服务的权衡很大部分是由于人类的干扰和社会的发展,拉大了不同生态系统服务此消彼长的趋势,因此需要从时空尺度合理评估生态系统服务。尽管众多的学者研究了不同尺度下生态系统服务之间的关系,但这些研究仅仅静态的分析了生态系统服务之间的关系；还缺乏对生态系统服务整体的差异以及长时间序列的动态观测研究。

我国是世界上生态工程实施最广泛的国家,贡献了全球植被变绿面积的25%[23]。黄土高原是我国植被变绿增幅最大的区域之一,也是退耕还林还草、治沟造地、固沟保塬、淤地坝建设等重大生态工程实施的典型区域；在这些生态工程作用下,黄土高原植被覆盖度明显增加,其带来的各项生态系统服务(储水、固碳、生物多样性等生态系统服务)也显著提升[24—27]。有研究显示,我国各类生态系统在2000—2010年期间,水源涵养、土壤保持与固碳、防风固沙、洪水调蓄和食物生产功能得以改善和提升,而生物多样性功能表现为下降趋势,其中退耕还林还草工程和自然区保护政策对生态系统服务起到了关键的作用[28]。随着植被的快速恢复,人类对生态系统作用力的不断增强,生态系统服务过度开发,使人类福扯受到威胁,大部分生态系统和生态系统服务正在或者趋于退化,生态与环境问题大量涌现,由此引发了一系列的水土流失、沙漠化严重、灌丛入侵、水资源平衡破坏、多样性降低、保护区生境质量下降等,因此,生态系统服务权衡和供给也成为研究的前沿与热点领域[29—32]。当前在“一带一路”、“退耕还林/草工程”、“美丽中国建设”、“黄河高质量发展”、“双碳战略”等一系列重大战略的推动下,黄土高原作为我国天然的西部生态屏障,该区域的生态环境建设与植被恢复显得极为重要。因此,合理评估黄土高原植被恢复条件下生态系统服务,是推动我国西部生态环境建设的重大国家需求和科学保障。

以黄土高原植被恢复为研究背景,在黄土丘陵沟壑区选取两个不同治理模式的典型小流域(纸坊沟流域,人工治理模式；坊塌流域,工程治理模式)；其中纸坊沟和坊塌流域都是延河流域的支流,两个流域处于不同治理阶段和模式,地理位置基本相毗邻,隶属于中国科学院水土保持研究所野外站台,2个小流域是黄土丘陵区退耕还林/草工程以来植被恢复和治理模式的典范代表。借助CASA(Carnegie-Ames-Stanford approach)模型,综合1998、2008和2018年的遥感影像图、土壤数据、气象数据、归一化差值植被指数(NDVI)等多源数据,利用实测数据、遥感解译数据、退耕还林数据、降雨数据、社会经济等数据,采用InVEST模型对生态服务(碳储量、土壤保持量、生境质量和产水量)进行评估,探究植被恢复过程中各生态系统服务的协同效应和权衡效应,这对于科学评价人类活动对黄土高原地区生态环境和可持续发展具有重要的现实和理论意义。

1 材料与方法

1.1 研究区概况

纸坊沟是黄土高原延河流域的支流之一,主要流经地形复杂的黄土丘陵沟壑区,该流域分布在纸坊沟村等3个自然村,位于延安市安塞县境内；坊塌流域分布于安塞县境内的沿河湾镇坊塌自然村,与纸坊沟流域一样同属于延河流域的重要支流之一。对于这2个小流域而言,由于环境治理时间及模式等存在差异,其生态系统服务具有一定的差异性,二者均位于安塞县东南区域,属于中科院水保所的重点流域研究观测对象,2个小流域是黄土丘陵区退耕还林/草工程以来,黄土丘陵沟壑区植被恢复和治理模式的典范代表。

纸坊沟流域呈现出明显的南北狭长分布特点,最低海拔为1010 m,最高海拔为1432 m,面积达到8.2 km2；坊塌流域面积略大于纸坊流域,达到8.6 km2,海拔与纸坊流域较为接近,最低海拔1023 m,最高海拔1381 m；受地理分布因素的制约,二者呈现出典型的暖温带半干旱季风气候,通过该地区20年气象资料分析,其年均温达到8.8℃,年日照达2400 h[33—34]。

1.2 植被调查和样品采集

2016年7月—10月进行在纸坊沟和坊塌流域植被调查、土壤样品采集等工作。在2个小流域选择6种不同土地利用方式,其中包括耕地、退耕草地、人工草地、人工灌丛、人工林地、自然灌丛；如表1所示,结合不同的土地类型采取平行样地的设置方法,在各样地内布设3个样方；植被调查参照方精云等人[35]的植被清查统计方法,林地、灌丛样方分别设置为10 m和3 m的正方形,而对于草本及裸地设置1 m×1 m的样方；记录植物群落结构特点。利用五点混合法收集土壤样品,在记录海拔及坡度等信息后,将土壤样品带回实验室测定理化性质。其他基础数据的搜集来源于文献和站点监测数据。

坡向以顺时针方向,每45°作为一个等级(范围在0—8之间)。数字1、2、3分别代表上坡位、中坡位和下坡位

1.3 遥感数据来源

(1)遥感影像获取

纸坊沟和坊塌流域2018年遥感影像资料的获得,通过黄土高原科学数据中心提取,对于1998年和2008年的遥感影像采取购买方式；空间分辨率为20 m,解译精度大于84%。

(2)DEM数据

纸坊沟和坊塌流域内DEM数据的获取采取矢量化及等高线转换的方式,利用ArcGIS (Arc Hydro模块)填挖处理DEM,使DEM表面平滑,坡度数据借助于空间分析模块计算,通过流域边界的提取获得高程数据。

(3)土地利用数据

结合纸坊沟和坊塌流域的遥感影像资料及Quick Bird数据,借助于墨卡托投影系统,在结合实地调查统计的基础上,进行对比分析。

(4)气象数据

通过黄土高原安塞国家气象站获取相应的气象资料,主要是气温、降雨量及太阳辐射等方面,此外,还结合统计年鉴及相关文献资料。

1.4 InVEST模型评估

采用InVEST模型和ArcGIS 10.2软件计算典型小流域的生境质量、碳储量、产水量、土壤保持服务,以及空间分布特征。

(1)NPP(Net Primary Production)

NPP指的是植被的净初级生产力,代表碳储存服务。采用CASA模型对NPP进行计算[36]：

NPP (x,t)=APAR (x,t) ×ε(x,t)

式中：x表示栅格中的具体像元；t表示月份；NPP(x,t)、APAR(x,t)、ε(x,t)分别表示t月份时,x像元上的植被净初级生产力、植被吸收的光合有效辐射、植被的实际光能利用率。

(2)碳储量

碳储量的估算借助于InVEST模型,碳储量主要包括大气碳储量、地上碳储量、土壤碳储量、地下碳储量、死亡有机碳。由于大气碳难以进行有效衡量,因此,本研究过程中对此不予估算。碳库模型的计算建立在不同土地利用类型的基础上,首先对碳密度进行测算,之后再结合研究区域面积,进而获取总的碳储量。具体的计算公式如下[37]：

C总碳=C地上碳+C土壤碳+C地下碳+C死亡碳

(3)产水模块

利用水循环作用机理测算产水量,具体的参数包括降雨量、径流、蒸发蒸腾、植被用水及土壤深度等影响因素,其中的Water Yield模块运行模式采取的是栅格单元,对于每个栅格单元而言,首先测定降雨量与蒸发蒸腾效应之差,而流域尺度采取模块加权的方式获取相应的产水量。

Yxj=(1-AETxj/Px)×Px

式中,Yxj、AETxj分别表示第j类土地利用类型在x单元格的年均产水量、年实际蒸散量；Px是x单元格的年降雨量。

(4)土壤保持模块

在InVEST模型中,采用土壤流失方程(USLE)计算土壤保持量,USLE最早由美国农业部提出,1980年引入中国,我国大量学者对此进行了一系列的改进和验算；主要原理是潜在土壤侵蚀与实际土壤侵蚀之差。计算公式如下[38]：

ΔUSLE=RKLS-USLE=R×(1-C-Ps)

式中,ΔUSLE是年平均土壤保持量(t hm-2a-1)；USLE是植被覆盖下的土壤流失量；RKLS是裸地土壤流失量；LS是坡长与坡度因子,根据DEM模型获得；C代表植被与管理因子；Ps代表土壤保持措施因子。

R是降水径流侵蚀力因子(MJ mm hm-2h-1a-1)),计算模型表达如下[38]：

式中,Pi表示月平均降雨量；P表示年平均降雨量；单位均为mm；

R的单位是100 ac h a-1,转换成国际单位MJ mm hm-2h-1a-1需乘以17.02的转换系数。

K指的是土壤可侵性因子(t h MJ-1mm-1),与土壤质地、降雨径流、有机碳、土壤渗透等有关,通过EPIC(Erosion/Productivity Impact Calculator)模型计算[38]：

{1.0-0.7×(1-SAN/100)/[1-SAN/100+exp[22.9×(1-SAN/100)-5.51]]}

式中,SIL、SAN、CLA分别指的是粉粒、沙粒和黏粒含量(%)；SOC指的是有机碳含量(%)。

1.5 生态系统服务权衡值的量化

逐像元的提取4种生态系统服务后,计算两两之间的相关性系数,当系数为负值时,认为两者属于权衡关系,当为正值时,则为协同关系。在空间上,通过MATLAB软件开展逐像元的长时间序列的相关性分析。计算公式如下[39]：

式中：R是两种服务之间的相关系数；i为像元的行号；j为像元的列号；n为栅格数据的时间序列；ES1、ES2为生态系统服务。

1.6 生态系统服务的生产可能性边界

生产可能性边界又可以称为效率曲线(Efficiency frontiers)或者帕累托曲线(Pareto frontiers)。在坐标轴上绘制出的图形称为生产可能性曲线(Production Possibility Frontier, PPF),或转换线[40—42]。PPF是对经济社会所能达到的最大产量组合的描述,适合用于筛选各种生产组合。PPF曲线内的点表示仍有资源未得到完全利用,存在资源闲置,说明生产还有潜力；而曲线之外的点,则表示现有的技术和资源条件是达不到的；只有曲线上的点,才是资源配置效率最高的点。本研究采用PPF来分析黄土高原不同生态系统服务之间的关系。先对不同生态系统服务两两图层进行比值,按照比值大小将每个单元格从小到大排列,然后依次对单元格对应的地理位置累计求和,绘制出相应的曲线,即为PPF曲线。

2 结果分析

2.1 NPP变化趋势

图1显示了纸坊沟流域1998—2018年平均值空间分布特征。1998—2018年平均植被NPP自中部向北部和南部递增,北部和南部植被NPP变化幅度较大,中部处于轻微变化区域；1998年纸坊沟流域植被NPP分布不均匀,此时NPP大多集中在15—40 g/m2之间,2008年NPP较均匀,此时NPP大多在20—80 g/m2之间,2018年NPP分布较为均匀,此时NPP大多在50—100 g/m2之间。

图1 纸坊沟流域1998—2018年NPP分布图Fig.1 NPP from 1998 to 2018 in Zhifanggou

图2显示了坊塌流域1998—2018年平均值空间分布特征。1998—2018年平均植被NPP自中部向四周递增,东部和南部植被NPP变化幅度较大,中部处于轻微变化区域；1998年NPP分布较为均匀,此时NPP大多集中在20—40 g/m2之间,2008年NPP大多在20—60 g/m2之间,2018年NPP大多在50—110 g/m2之间。

2.2 土壤保持功能评估

通过图3可知,纸坊沟流域1998—2018年单位面积土壤保持量呈上升趋势,不同区域存在较为明显的差异,其中变化较大的是南部区域,而中部变化较小；1998年单位面积土壤保持量并没有较为明显的变化,其低值为45 t/hm2,高值为232 t/hm2；2008年呈现出中部高的特点,南北区域则变化小,最低值为53 t/hm2,最高值为299 t/hm2；2018年则保持较为均匀的变化趋势,最低值为78 t/hm2,最高值为320 t/hm2。

通过图4得知,坊塌流域单位面积土壤保持量存在一定的时空差异,整体呈上升趋势,且变化幅度较小,呈较为均匀的分布特点；1998年单位面积土壤保持量并没有较为明显的变化,最低值为22 t/hm2,最高值为190 t/hm2；2008年呈现出中部高的特点,南北区域则变化小,最低值为35 t/hm2,最高值为221 t/hm2；2018年则保持较为均匀的变化趋势,最低值为46 t/hm2,最高值为267 t/hm2。

图2 坊塌流域1998—2018年NPP分布图Fig.2 NPP from 1998 to 2018 in Fangta

图3 纸坊沟流域1998—2018年单位面积土壤保持量Fig.3 Soil conservation per unit area from 1998 to 2018 in Zhifanggou

图4 坊塌流域1998—2018年单位面积土壤保持量Fig.4 Soil conservation per unit area from 1998 to 2018 in Fangta

2.3 产水量功能评估

通过图5得知,纸坊沟流域1998—2018年单位面积产水量呈逐步增加的趋势,南部区域波动不大,而北部地区变化较大；1998年单位面积产水量南北差异较大,南部区域产水量较高(43 mm/hm2),北部产水量较低(317 mm/hm2)；2008年单位面积产水量南北差异较大,最低值为33 mm/hm2,最高值为326 mm/hm2；2018年单位面积产水量变化趋势呈现出典型的南北小、中部大特点,最低值为26 mm/hm2,最高值为356 mm/hm2。

通过图6得知,坊塌流域1998年单位面积产水量空间分布差异较大,南部、北部区域较高,而中部较低；2008年呈现出北低南高的特点,最低值为17 mm/hm2,最高值为298 mm/hm2；2018年由北部向南部区域的产水量呈现上升趋势,最低值为23 mm/hm2,最高值为324 mm/hm2。

图5 纸坊沟流域1998—2018年单位面积土壤产水量Fig.5 Soil water yield per unit area from 1998 to 2018 in Zhifanggou

图6 坊塌流域1998—2018年单位面积土壤产水量Fig.6 Soil water yield per unit area from 1998 to 2018 in Fangta

2.4 碳储量功能评估

纸坊沟流域单位面积碳储量在不同年份存在较大的差异(图7)。1998年单位面积碳储量较高的是北部区域,而南部较低,最低值为1.59 t/hm2,最高值为13.26 t/hm2；2008年单位面积碳储量呈现出典型的南北高、中部低的变化特点,最低值为3.02 t/hm2,最高值为16.99 t/hm2；2018年碳储量也呈现出典型的南北高、中部低的变化特点,最低值为2.47 t/hm2,最高值为27.26 t/hm2。

从图8可知,坊塌流域1998年单位面积碳储量分布具有明显的地域差异,南部及北部区域含量较高,而中部低,最低值为1.02 t/hm2,最高值为1.53 t/hm2；对于2008年,由北向南递减的变化趋势尤为突出,最低值为1.23 t/hm2,最高值为14.63 t/hm2,呈现较大的差异；而对于2018年,由北向南递减,最低值为2.01 t/hm2,最高值为19.52 t/hm2。

图7 纸坊沟流域1998—2018年单位面积碳储量Fig.7 Soil carbon storage per unit area from 1998 to 2018 in Zhifanggou

图8 坊塌流域1998—2018年碳储量Fig.8 Soil carbon storage per unit area from 1998 to 2018 in Fangta

2.5 生境质量评估

结合生境质量的分布特点,纸坊沟流域生境质量呈逐渐上升趋势(图9)；1998年生境质量指数较高的是中部区域,最低值为0.06,最高值为0.81；2008年,由北向南递减趋势,最低值为0.05,最高值为0.92；2018年呈两端高、中部较低的特点。坊塌流域1998年生境指数变化并不大(图10),最低值为0.09,最高值为0.77；对于2008和2018年,生境指数最低值分别为0.16、0.15,最高值分别为0.79、0.87,呈现出较为均匀的分布特点。

图9 纸坊沟流域1998—2018年生境质量Fig.9 Habitat quality from 1998 to 2018 in Zhifanggou

2.6 不同土地利用方式间权衡与协同特征

由图11可知,纸坊沟流域林地、草地、梯田和灌丛,土壤保持量和碳储量以协同关系为主,协同关系比例超过了50%；道路、建筑用地、耕地和裸地,空间方面的权衡关系比例超过了50%。对于裸地和道路,产水量与土壤保持呈权衡关系,在空间方面的权衡关系比例超过了50%。

由图12可知,坊塌流域林地、草地、梯田和灌丛,土壤保持和碳储量表现为协同关系；所有土地利用类型下产水量和土壤保持表现为权衡关系；对于建筑用地和道路,土壤保持和生境质量以权衡关系为主,其他土地利用类型更多的呈现协同关系。对于建筑用地和道路,碳储量、产水量呈现权衡关系。

2.7 生态系统服务优化

使用Python语言绘制产水量、碳储量、生境质量、土壤保持之间的生产可能性边界(PPF)曲线,然后从定量角度研究生态系统服务之间的关系。由图13可知,纸坊沟流域土壤保持量与碳储量之间协同关系明显,从点a到b和c到d,土壤保持增加量分别为43.49 t/hm2和37.53 t/hm2,碳储量分别增加了2.97 t/hm2和3.16 t/hm2。产水量与碳储量、土壤保持量的PPF曲线均是“凹”趋势(协同关系)。也即：随着产水量的增加,碳储量、土壤保持量增加的更明显。而生境质量与产水量、土壤保持量、碳储量的PPF曲线均是“凸”趋势(权衡关系),也即：随着生境质量逐渐提升,产水量、土壤保持量、碳储量逐渐减小。

由图14可知,坊塌流域土壤保持量和碳储量协同关系表现明显,从点a到b和c到d,土壤保持增加量分别为44.53 t/hm2和60.59 t/hm2,碳储量分别增加了2.72 t/hm2和5.40 t/hm2,这说明碳储量累计越多,对土壤保持越有利。同样地,产水量与碳储量、生境质量与土壤保持量、生境质量与产水量之间的PPF曲线表现为向内“凹”的趋势,也即随着碳储量的增加,产水量和土壤保持量呈同增同减的协同关系。然而,产水量与土壤保持量、生境质量与碳储量之间的PPF曲线均是“凸”趋势(权衡关系),也即随着碳储量累积的增加,生境质量逐渐减小；随着土壤保持量逐渐累积,产水量也逐渐减小。

图11 纸坊沟流域不同土地利用权衡协同统计Fig.11 Tradeoff and synergy of different land use types in Zhifanggou

图12 坊塌流域不同土地利用权衡协同统计Fig.12 Tradeoff and synergy of different land use types in Fangta

图13 纸坊沟流域生态系统服务帕累托效率曲线Fig.13 Pareto efficiency curve of ecosystem services in Zhifanggou图中从点a到点b的变化量表示生态系统服务的机会成本,从点c到点d的变化量表示生态系统服务所需要付出的计划成本

图14 坊塌流域生态系统服务帕累托效率曲线Fig.14 Pareto efficiency curve of ecosystem services in Fangta

3 讨论

近些年,随着“两山”、“双碳”、“美丽中国”和“黄河高质量发展”等重大战略的提出,我国对黄土高原生态环境的重视程度越来越高,与之相关的生态系统服务受到了更多的关注[43—46]。本文通过InVEST模型,对黄土丘陵沟壑区1998—2018年纸坊沟和坊塌流域的产水、固碳、水土保持和生境质量4种服务进行分析,探究植被恢复过程中各生态系统服务的协同效应和权衡关系。研究发现,纸坊沟与坊塌流域NPP在1998—2018年间逐渐增加,其中前十年增加微弱,后十年增加急剧,这一结果与前人对黄土高原NPP时间变化的结果一致[47—48],表明退耕还林后,土地利用结构发生了改变,植被覆盖度明显提高,NPP也随之增加；其次,植被生长具有一定的周期性,耕地面积的减小并不能瞬时带来NPP的增加,因此NPP的增加与植被恢复过程也具有一定的滞后性,其中前十年为植被缓慢恢复期,后十年为植被显著恢复期。

降水与蒸发是决定产水量的两个主要因子,气象数据表明,自1998年以后黄土高原地区气温波动上升使得蒸发量增强,且1998—2018年降水量增加显著,两者的增加导致了产水量的增加[47],与本研究的结果相一致。此外,碳储量、土壤保持量、生境质量1998—2018年均逐渐增加,增加幅度有所不同,主要是由于随着退耕还林还草工程的实施,土壤养分和质量得到了明显的提升；对于1998—2008年,作为退耕的第一阶段,大量耕地转变为林地和草地,植被得到逐步恢复,而后十年是植被的快速恢复期[49],林地和草地面积明显增加,生态系统服务得到了大幅提升,尤其是在碳储量和土壤保持功能方面。此外,产水量和生境质量的提高促进了多项生态系统服务的提升[49]。对于2008—2018年,作为退耕的后十年,耕地面积的下降幅度有所减缓,相应地,林地和草地的增加速度减慢,生态系统服务的提升速度减缓,尤其是碳储量与土壤保持,这与前十年相比出现较大差异；通过生境质量分析得知,随着退耕还林工程的不断实施,生境质量得以提升。从空间上看,纸坊沟流域产水、碳储量、土壤保持、生境质量呈现出南北高、中部低,而坊塌流域各生态系统服务由北向南递减,主要是由于纸坊沟流域南北部和坊塌流域的南部森林、灌木分布广泛,植被覆盖度高,大量的枯枝落叶对水分起到了拦截左右,这样可以促进水源涵养的提高；与此同时,庞大的植物根系不断传输水分给土壤,使土壤结构更加稳固[50]；本研究还发现：森林、灌丛和草地碳储量较大,主要是由于植被的光合作用固碳、以及根系固碳使得生态系统碳汇增加,伴随着土壤碳储量的增加,土壤肥和水分均得以提升,进而减少了土壤侵蚀和增加了土壤保持量[51]。与此同时,这些区域大部分属于保护区,受到人类的管控,生态系统类型多样化,生境质量和多样性都维持在较高水平。

生态系统服务的权衡和协同关系可知,对于纸坊沟流域,不同年份的生态系统服务存在较大差异,其中2008和2018年呈现出明显的权衡关系,而在1998年表现出协同关系。对于坊塌流域,1998和2018年呈现出权衡关系,而2008年表现为协同关系。低协同是一种比较稳定的状态,对于每一种服务而言,独立或者综合作用结果都较差,1998年纸坊沟流域基本表现出低的协同关系,主要是由于其脆弱的生态环境与人类的活动对环境的破坏作用导致。高协同关系一般各项生态系统服务较好,各项生态系统服务综合表现为互利协同模式[50—51],比如2018年各项生态系统服务以高协同为主。相对比而言,生态系统服务权衡关系较弱的区域,主要是由于部分生态系统服务较高导致其他生态系统服务较低,生态系统服务波动较大,彼此之间存在激烈竞争的状态[52—53]。无论是1998年的坊塌流域,还是2008年纸坊沟流域,高权衡关系的区域分布较为广泛,高权衡关系主要是因为地形起伏不定,林分质量较差以及人类活动的干扰,使得各生态系统服务呈明显的权衡关系[38—39]。针对这些区域,需加强对林分的管理以提升其固碳量、土壤保持量和生境质量,以缓解较强的权衡关系；此外,还可以结合当地经济发展,适当的改变粗放的农业管理模式,用以提升森林资源的生态系统服务。

纸坊沟、坊塌流域生态系统服务之间的关系存在明显的空间差异性。协同区域面积比例虽大,但其中又混合着权衡区域,权衡和协同区域呈絮状交错分布[54]。胡秀芳等[55]对黑河流域的研究表明：生态系统服务之间的关系南北区域性差异很大,但仍然有30%的区域没有权衡和协同关系。张静静等[39]研究了伏牛山森林生态系统服务权衡与协同关系,结果发现：高权衡关系的森林面积约占61%,高协同关系的森林面积约占29%。在荒漠草原,碳储量与土壤保持呈协同关系,而碳储量与水源涵养的协同关系均不明显[52]。在青藏高原,气候调节与碳储量表现为协同关系,而产水量与碳储量表现为权衡关系[56]。这些研究与本研究不相一致,可能不同地区的地形类型、气候类型和人类活动等影响不尽一致[57-58]。此外,本研究处于黄土丘陵沟壑区,受复杂的植被、地形地貌及气候条件的影响,纸坊沟和坊塌流域生态系统服务之间的协同和权衡关系呈明显的空间差异分布。由此,在退耕还林还草工程实施过程中,要充分结合不同区域的生态系统服务开展针对性的保护,在协同或权衡效应的基础上,提出针对性保护措施,并维持生态系统的平衡和稳定性发展。

本项研究从静态的时间点上探讨了1998年、2008年和2018年黄土丘陵沟壑区典型流域生态系统服务之间的权衡与协同关系,具有明显的区域分异和时空分异特征,然而,全面理解生态系统服务权衡和协同区域和时空尺度仍是一大难题[59]。生态系统服务之间的关系是随时间动态变化的,未来更应加强权衡与协同的连续性和动态转变特征研究。另外,生态系统服务权衡和协同关系中蕴含了许多生物过程、物理机理以及化学机制,如果缺乏对生物、物理、化学过程的深入剖析,这将会影响对生态系统服务的驱动机制的认知。因此,在研究生态系统服务之间的关系时,应更多的涉及生物、物理、化学过程的解释,以提升生态系统的平衡和稳定,促进可持续发展。

4 结论

(1)随着植被的恢复,坊塌和纸坊沟NPP和各项生态服务均得到了大幅的提升,其中1998—2008年增幅较大(生态系统服务显著提高),而2008—2018年增幅较小(生态系统服务提高)。

(2)在植被恢复20年(1998—2018),坊塌和纸坊沟流域受自然因素及人为因素的制约,各生态系统服务之间主要表现为协同关系,局部区域表现出权衡关系；这主要是由于植被恢复造成了各生态系统服务的增加,然而在时间和空间上增幅不同,从而表现出以协同关系为主的分布模式；因此,在退耕还林的生态环境保护过程中,需要基于各生态系统服务之间的关系,因地制宜的划分和制定保护策略。