退耕还林还草背景下陕北地区生态系统服务动态及热点识别

王芸， 赵鹏祥， 李治国， 夏冰， 高健洲， 马鹏宇

(1.河北开放大学 教务处，河北 石家庄 050080； 2.国家开放大学 农林医药教学部，北京 100039；3.西北农林科技大学 林学院，陕西 杨凌 712100； 4.河北博嘉农业有限公司，河北 石家庄 050000；5.北京林业大学 园林学院，北京 100083； 6.河南职业技术学院 环境艺术工程学院，河南 郑州 450046；7.北京铁汉生态环境科学研究院有限公司，北京 100083)

生态系统服务是指生态系统和生态过程所提供的维持人类生命的自然条件和效用[1]。千年生态系统评估结果表明，大约60%的生态系统服务出现了恶化[2]。全球多个国家已经实施了许多生态恢复措施以改善区域生态系统服务。我国实施的退耕还林还草工程是发展中国家中最大的生态修复计划，也是世界上最大的生态建设项目[3-4]。了解退耕还林还草背景下生态系统服务的时空变化特征及相互作用，确定多个生态系统服务的热点区域，对于制定新一轮的退耕还林还草规划和土地资源合理配置有重要的参考意义。

自千年生态系统评估之后，生态系统服务研究受到了地理学、生态学等相关学科领域专家的广泛关注，研究成果丰硕，研究内容主要涉及生态系统服务的概念和内涵[1-2,5-6]、生物物理量化和价值化[7-8]、时空变化和驱动机制[9-11]、权衡与协同作用[10,12-13]、供需平衡与空间流动[14-15]。如：李燕和周亮广[16]从服务强度、区位指数视角对安徽省市域范围的生态系统服务价值时空变化进行了分析；ZHANG Z等[17]以福建省为例，研究了多种生态系统服务之间的相互作用，并结合地理加权回归模型确定了城市化、气候和地形等对权衡与协同作用的空间影响力；YANG S等[18]通过构建不同的土地利用变化情景，对延河流域的生态系统服务进行量化，揭示了土地利用变化与生态系统服务之间的相互关系和影响。上述研究成果为区域生态环境建设和生态管理决策的制定提供了重要依据。然而，鉴于生态系统的复杂性，对典型地区生态系统服务的解析在生态系统管理领域中仍需进一步加强。

陕北地区作为黄土高原乃至中国退耕还林还草工程实施的重要区域之一，拥有脆弱的生态环境和高强度的人类活动，是开展生态系统服务相关研究的理想区域和典型区域。以陕北地区为例，分析退耕还林还草工程实施以来土地利用的变化特征，结合多种模型量化陕北地区4种关键生态系统服务(产水量、粮食产量、植被净初级生产力(net primary productivity，NPP)和土壤保持)的时空变化和相互作用，识别生态系统服务的热点区域，以期为该地区进一步的生态恢复治理提供理论依据。

1 研究数据与方法

1.1 研究区概况

陕北地区，位于35°21′～39°34′N、107°28′～111°15′E，辖榆林市和延安市的所有区域，面积为78 819.69 km2，海拔为417～1 907 m，如图1所示。陕北地区位于黄土高原中部，以典型黄土丘陵和沟壑为主要地形特征，气候以半干旱大陆性季风气候为主，年平均气温为8～12 °C，年平均降水量为250～600 mm，降水量年内分布不均，大约有60%～70%的降水集中在7—9月份[19]。土壤以黄绵土分布范围最广，其次为草原风沙土，土壤质地疏松，易受侵蚀。在退耕还林还草工程启动之前，集中降水和不适当的人为土地利用活动是导致该地区生态环境脆弱，特别是水土流失严重的主要原因[20]。经过大规模的退耕还林还草，该地区的生态环境得到了明显改善，生态系统服务发生了巨大变化。

图1 研究区位置

1.2 数据基础

以退耕还林还草工程实施的关键时间节点(2000年和2015年)之间的时间段为研究期，进行相关研究数据的收集与处理。涉及的数据主要有：

1)土地利用数据，来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)，空间分辨率为30 m。在统计土地利用变化时借助ArcGIS 10.6重分类工具将其划分为耕地、林地、草地、水域、建设用地和未利用地6类，结果如图2所示。

图2 2000年和2015年陕北地区土地利用类型空间分布示意图

2)数字高程模型(Digital Elevation Model，DEM)数据，来源于地理国情监测云平台 (http://www.dsac.cn),空间分辨率为30 m。

3)月尺度的归一化植被指数(Normalized Difference Vegetation Index，NDVI)数据，来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)，空间分辨率为1 000 m。

4)土壤数据，来源于世界土壤数据库(Harmonized World Soil Database，HWSD)的中国土壤数据集(v1.1)(http://westdc.westgis.ac.cn),空间分辨率为1 000 m。

5)月尺度的气象数据，来源于国家气象科学数据中心(http://data.cma.cn/)，通过克里金插值法生成1 000 m分辨率的栅格数据。

6)粮食产量数据，来自于相关县级行政单元的统计年鉴(https://data.cnki.net/)。

7)人口和GDP数据，来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)，空间分辨率为1 000 m。

1.3 生态系统服务量化方法

1.3.1 产水量

采用生态系统服务和权衡的综合评估模型(Integrated Valuation of Ecosystem Services and Trade-offs，InVEST)模型中的产水量模块对陕北地区2000年和2015年的产水量进行量化。产水量模块是基于水量平衡原理，以像元内的降水量与实际蒸散量之差作为产水量的计算结果[9]，其基本原理如下：

WYx=(1-AETx/Px)Px，

(1)

(2)

(3)

PETx=KxETx。

(4)

式中：WYx为像元x的年产水量；AETx为像元x的年实际蒸散量；Px为像元x的年均降水量；PETx为像元x的年潜在蒸散量；wx为自然气候-土壤性质的非物理参数；Z为季节参数，其数值与降水次数呈正相关；AWCx为像元x的植被可利用含水量；Kx为像元x的植被蒸散系数；ETx为像元x的参考作物蒸散量。

1.3.2 粮食产量

由于粮食产量与NDVI之间存在着显著的线性关系[21]，采用全部为正值的NDVI对粮食产量进行空间化处理，得到粮食产量栅格数据(Gx)。具体的计算公式为：

(5)

式中：Gsum,i为县级行政单元i的粮食产量；NDVIx为像元x的NDVI；NDVIsum,i为县级行政单元i耕地的NDVI总和。

1.3.3 植被净初级生产力(NPP)

采用基于过程的卡内基—埃姆斯坦福方法(Carnegie-Ames-Stanford Approach，CASA)模型计算陕北地区2000年和2015年的NPP，CASA模型是一种陆地生态系统评估模型，在NPP的评估研究中应用广泛[10,12]，模型如下：

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

(6)

APAR(x,t)=0.5SOL(x,t)FPAR(x,t)，

(7)

ε(x,t)=Tε1Tε2Wεε*。

(8)

式中：APAR(x,t)为像元x在t月吸收的光合有效辐射；ε(x,t)为像元x在t月份的光能转化率；SOL(x,t)为像元x在t月份的太阳总辐射；FPAR(x,t)为像元x在t月份植被对光合有效辐射的吸收比例；ε*为植被最大光能利用率；Tε为对ε*产生影响的温度胁迫因子；Wε为水分胁迫因子。

1.3.4 土壤保持量

采用修正的通用土壤流失方程(Revised Universal Soil Loss Equation，RUSLE)模型以潜在土壤侵蚀量与实际土壤侵蚀量之差来表示生态系统中的土壤保持量[22]。其基本原理为:

As=Ap-Ac=RKLS(1-CP)。

(9)

式中：As为土壤保持量；Ap为潜在土壤保持量；Ac为实际土壤侵蚀量；R为降雨侵蚀力因子；K为土壤可侵蚀性因子；L为坡长因子；S为坡度因子；C为植被覆盖和管理因子；P为水土保持措施因子。

1.4 协同作用与权衡判定依据

生态系统服务之间的协同指的是生态系统服务同时增加或减少的情况，而权衡是指一种生态系统服务的增加导致另一种生态系统服务减少的情况[23]。本研究在SPSS 24.0软件平台上结合Spearman相关系数法分析了陕北地区生态系统服务之间的协同作用与权衡。当Spearman相关系数为正值且通过显著性检验时，说明两种生态系统服务之间存在着协同作用；当Spearman相关系数为负值且通过显著性检验时，说明两种生态系统服务之间存在着权衡作用。

1.5 服务热点识别方法

将产水量、粮食产量、NPP或土壤保持量超过其平均值的区域定义为该生态系统服务的热点区域[24-26]，将这4种生态系统服务的热点区域在ArcGIS 10.6内进行空间叠加，得到多个生态系统服务的热点区域。如果某个像元中的这4种生态系统服务的量化值均不超过对应生态系统服务在整个研究区的平均值，则该像元为无热点区域。同样，如果某个像元中，只有1种生态系统服务的量化值超过了该生态系统服务在整个研究区的平均值，则将该像元划分为一类服务热点(1类热点)；如果某个像元中，只有2种生态系统服务的量化值超过了对应生态系统服务在整个研究区的平均值，则将该像元划分为二类服务热点(2类热点)，以此类推，划分为三类服务热点(3类热点)和四类服务热点(4类热点)。

2 结果与分析

2.1 退耕还林还草后土地利用类型的变化

退耕还林还草工程的实施对陕北地区的土地利用产生了重大影响， 2000年和2015年陕北地区各种土地利用类型的面积占比如图3所示。

图3 2000年和2015年陕北地区各种土地利用类型的面积占比

图2和图3表明：陕北地区的草地分布最广，在2000年和2015年分别占总面积的43.53%和45.61%；其次为耕地和林地的，耕地和林地在2000年的面积占比分别为35.82%和13.86%，在2015年的面积占比分别为31.90%和15.38%；最后依次为未利用地、建设用地和水域，他们的面积占比非常小，2000年和2015年两个年份的未利用地面积占比均不足5.8%，建设用地和水域的面积占比均不足1.0%。自退耕还林还草工程实施以来，草地、林地、耕地和建设用地面积变化很大，存在着草地、林地、建设用地面积增加和耕地面积减少的变化趋势。

陕北地区实施退耕还林还草以来土地利用类型相互转化的情况见表1。

表1 陕北地区各种土地利用类型面积转移矩阵 km2

由表1可知：研究期内，各种土地利用类型均存在向其他土地利用类型转化的变化特征。具体来说，耕地主要存在着向草地、林地和建设用地转化的特征，转化面积分别为3 699.65、1 028.70、176.43 km2；林地向其他土地利用类型的转化面积非常小，仅有不足400.00 km2的林地面积转化为其他土地利用类型；草地存在着向耕地、林地和建设用地转化的轻微趋势，转化面积对于草地总面积来说占比很小；2015年，建设用地面积增加一倍多，主要是耕地和草地向建设用地的转化，转化面积分别为176.43、196.97 km2；未利用地和水域的面积变化较小，对于其他土地利用类型面积变化的贡献非常微小。

2.2 生态系统服务的时空变化

利用ArcGIS 10.6绘出4种生态系统服务的空间分布情况，如图4所示。

图4 2000年和2015年陕北地区生态系统服务的空间分布

由图4可以看出，陕北地区的4种生态系统服务在研究期内各自的空间分布存在着明显的分异性。尽管部分生态系统服务(产水量、粮食产量和NPP)在2000年和2015年的空间分布存在着较小的差异，但4种生态系统服务在2个年份的总体分布特征是相似的。产水量在2000年和2015年均存在着自南向北递减的空间分布特征，其中在榆林西北部地区的产水量相对较高，特别在2015年，该地区的产水量达到最大值；粮食产量的空间分布比较分散，这与粮食产量的栅格化过程有关，即粮食产量的空间分布与耕地的空间分布一致，粮食产量的分布面积随耕地面积的减少而减少；NPP呈现出由南向北递减的空间格局，而且随时间的推移，陕北大部分地区的NPP得到了改善；土壤保持量在2000年和2015年这2个研究年份的空间分布比较相似，即自南向北土壤保持量逐渐减少，低值区位于陕北地区的西北部，呈条带状分布。

利用ArcGIS 10.6比较了陕北地区4种生态系统服务区域在2000年和2015年的整体水平，发现产水量、粮食产量和NPP呈现出增加趋势，分别由2000年的6.75 mm、17.37 t/(km2·年)、322.21 g·C/(m2·年)增加至2015年的21.35 mm、27.02 t/(km2·年)、439.49 g·C/(m2·年)，增幅分别为216.30%、55.56%和36.40%；而土壤保持量呈现出减少趋势，由2000年的36 637.67 t/(km2·年)减少至2015年的24 319.74 t/(km2·年)，减少了33.62%。

2.3 权衡与协同作用

采用Spearman相关系数描述产水量、粮食产量、NPP和土壤保持量两两之间的相互作用，结果见表2。

表2 2000年和2015年陕北地区4种生态系统服务之间的Spearman相关系数

由表2可以看出:2000年，4种生态系统服务两两之间的关系均呈现为显著的协同作用，其中以产水量与NPP之间的协同作用(0.576)最为明显，其次为NPP与土壤保持量(0.470)之间的、NPP与粮食产量(0.460)之间的，再次为产水量与粮食产量(0.420)之间的、产水量与土壤保持量(0.390)之间的，粮食产量与土壤保持量之间的协同作用最弱(0.103)。以上括号中的数值为Spearman相关系数值。随着退耕还林还草工程的实施，陕北地区4种生态系统服务两两之间的关系2015年的与2000年的相比发生了明显的变化，主要表现为同一生态系统服务关系对协同作用的减弱或逆转为权衡作用。具体来说，产水量、粮食产量与NPP这3种生态系统服务两两之间的协同作用出现衰减；土壤保持量与产水量、粮食产量之间的协同作用逆转为2015年的权衡作用；此外，土壤保持量与NPP之间的协同作用在2015年得到了加强，Spearman相关系数值由2000年的0.470增大至2015年的0.515。

2.4 多个生态系统服务热点区域

结合1.5节的热点区域识别方法，在ArcGIS 10.6中将研究区划分为无热点、1类热点、2类热点、3类热点和4类热点区域，各热点区域的分布和面积占比分别如图5和图6所示。由图5和图6可知：

1)从空间分布上来看，2000年和2015年各类型生态系统服务热点区域的空间分布情况整体上比较相似，4类热点区域和3类热点区域主要分布在研究区的南部地区，呈现出比较明显的空间集聚特征，2类热点区域主要分布在研究区的中部和中南部地区，1类热点区域和无热点区域主要分布在研究区的北部和西北部地区，呈相对分散的空间分布格局。

2)通过对比2000年和2015年各类型生态系统服务热点区域的空间分布可以发现，4类热点、3类热点和无热点区域的分布范围缩小，2类热点和1类热点区域的空间分布范围增大。

3)具体来说(如图6所示)，4类热点区域的面积由2000年的占研究区总面积的2.17%减少至2015年的1.35%，3类热点区域的面积占比由2000年的13.17%减少至2015年的10.64%，无热点区域的面积占比由2000年的27.80%减少至2015年的21.72%，2类热点区域和1类热点区域的面积占比分别由2000年的22.67%、34.19%增加至2015年的29.21%、37.08%。

图5 2000年和2015年各类型生态系统服务热点区域的分布情况

图6 2000年和2015年各类型生态系统服务热点区域的面积占比

3 结语

文中以退耕还林还草工程实施的重点区域——陕北地区为研究区，分析了退耕还林还草工程实施以来土地利用类型的变化特征，量化了该背景下产水量、粮食产量、NPP和土壤保持量的时空变化特征和这4种生态系统服务之间的权衡与协同作用，确定了多个生态系统服务的热点区域。研究发现：

1)2015年，陕北地区的土地利用类型主要以草地、耕地和林地为主，占研究区总面积的比例分别为45.61%、31.90%和15.38%。2000—2015年，研究区呈现出耕地面积减少(-4 976.03 km2)且主要转化为草地(+3 699.65 km2)和林地(+1 028.70 km2)的变化特征，这些变化趋势与该地区土地利用的其他研究结果相一致[27]，这也符合退耕还林还草工程实施的初衷。

2)陕北地区4种生态系统服务在研究期内的空间分布存在着明显的分布差异性。产水量、NPP和土壤保持量总体呈现出自南向北递减的空间分布特征；粮食产量的空间分布由于耕地分布分散的原因呈现出比较分散的空间分布格局。研究期内产水量、粮食产量和NPP呈现出增加趋势，而土壤保持量呈现出减少趋势，这些结果也与退耕还林还草背景下LIU Y等[28]、SU C等[29]的研究结果一致。退耕还林还草的实施，在改善地表植被覆盖情况[30]的同时也增强了土壤的保持能力，这不仅改善了NPP的供给能力，还抑制了潜在土壤侵蚀量的减少和实际土壤保持量的下降[31]。此外，在此背景下，耕地面积的减少是显而易见的，但由于21世纪以来国家对农业的政策支持和粮食单产能力的提高[4]，粮食产量呈现出了不减反增的趋势[4,21]。已有研究发现，退耕还林还草工程的实施减少了产水量[12,19]，这与本研究的结果相反，产水量与不同研究期内的降水情况有关，降水量的明显增加导致研究区的产水量增加[32]。

3)2000—2015年，陕北地区4种生态系统服务两两之间的关系发生了显著变化，表现为2000年所有生态系统服务关系对之间的协同作用在2015年出现衰减甚至逆转为权衡作用。退耕还林还草工程的实施，增强了土壤保持量与NPP之间的协同作用，同时也导致了土壤保持量与产水量、粮食产量之间的权衡，这与PAN J等[13]、潘竟虎等[25]的研究结果一致。需要注意的是，本研究证实了生态系统服务之间的相互作用会随时间的变化而变化，特别是在人类活动密集的地区[33-34]，单一时间节点上的生态系统服务权衡与协同作用的静态研究可能不足以为人类活动日益密集背景下的生态管理决策提供充足依据。

4)确定多个生态系统服务热点区域是确定生态保护区域优先级和优化土地利用结构的重要途径，这对政策制定者平衡资源分配和领土空间规划具有重要意义[25-26]。本研究基于生态系统的多服务功能评估，识别了生态系统服务热点区域，发现研究区自南向北呈现出4类热点区向无热点区逐渐过渡的空间分布格局，这与LIU L等[34]的研究结果类似。然而随着时间的推移，4类热点区域和3类热点区域的面积在缩小，这表明了生态系统服务之间的权衡在增加或者协同在减弱，难以同时获得多种生态系统服务的大量供给，这与文献[24]和[26]的研究结果一致，也与本研究得出的生态系统服务权衡与协同作用的时间变化特征相符合。

然而，文中尚存在一些不足。首先，各类数据的分辨率不一致，这会导致结果存在误差。例如，粮食产量的估算是基于县域粮食产量进行的，若能获得镇域甚至村域粮食产量数据，结果可能会更准确。其次，本研究仅在栅格尺度上考虑了生态系统服务之间的相互作用，已有研究发现，生态系统服务权衡与协同作用会随空间尺度的变化而发生变化[35]，多尺度下的生态系统服务相互作用的比较分析需要进一步加强。最后，热点区域识别方法是以平均值为分界线的，即高于生态系统服务平均值的区域被定义为热点区域，尽管很多研究也是采用这个平均值作为分界线的[24-26,34]，未来应结合区域发展的重点、因地制宜地考虑设置热点识别的阈值。