基于可持续发展的太行山区生态系统服务垂直分类管理*

高会，刘金铜，朱建佳，李晓荣

①中国科学院遗传与发育生物学研究所 农业资源研究中心/中国科学院农业水资源重点实验室/河北省节水农业重点实验室，石家庄 050022；②中国科学院大学，北京 100049；③河北科技师范学院 园艺科技学院，河北 秦皇岛 066004

中国是一个多山的国家，山地面积约占陆地总面积的70 %， 超过40 %的人口生活在山区[1-2]。山地为人们提供着多种直接或者间接的生态系统服务[3]，这些生态系统服务来源于人类和自然相互作用形成的山地生态系统，而山地生态系统对于快速的全球变化在生态和经济上都表现出敏感性[4]。在全球气候变化、经济一体化的条件下，山区的发展面临着众多挑战。其中，山区可持续发展成为山区研究的重要关注点[5-7]，优化调控管理是实现山区可持续发展的根本措施[8-9]。

不同空间尺度的生态系统服务的评估与空间分布已有不少研究[10-11]，Hayha等[12]评估了山区森林生态系统服务，并给出了其空间分布。但是，在山区生态系统服务研究中，人们忽略了山地垂直地带性[13]，只对其空间分布进行了一些描述[14]，并没有阐明垂直梯度上生态系统服务的特征。因此，在山区生态系统服务评估的基础上，明确其垂直梯度生态系统服务功能格局是非常有必要的。

垂直地带性分异是山地特有的一种现象，形成地貌 、植被和土壤等的垂直带以及自然景观垂直带[13]。目前垂直地带性分异的研究主要有垂直带谱和垂直带界线研究[15-16]。以前研究中的垂直带谱中垂直带和垂直带界线，通常以各类植被和土壤为主要标志，并结合水热条件和地貌特点来划分，不能全面反映山区生态系统的垂直分异特征。

太行山区东西横跨4个经度，南北纵跨6个纬度，是黄土高原和华北平原的天然分界线。黄土高原的部分河流经太行山流入华北平原，因此太行山区具有重要的地理意义和生态意义。太行山区由于海拔高度并不足够高，在垂直梯度上还没有表现出明显的垂直带谱分布。但是，根据对大量调研考察资料以及相关资料的整合归纳，发现太行山区在垂直梯度上表现出较明显的生态特征分异规律。研究太行山区生态系统分异特征，确定其垂直分异的关键线带，分析垂直梯度上生态系统服务分布特征，可以为太行山区生态系统服务梯度分类优化管理、实现区域可持续发展提供理论依据。

借鉴垂直带谱中的研究理论，本文以太行山区水、土、气、生的典型生态要素指标的高程变化为依据，确定了太行山区生态系统垂直分异的关键线带，进而将太行山区划分为亚高山区、中山区和低山丘陵区3类生态系统空间格局。根据生态系统服务梯度变化的差异性，界定太行山3类生态系统的生态系统服务垂直梯度上管理的可持续发展最优模式为：亚高山区生态系统服务以调节和支持功能为主；中山区生态系统服务同时兼顾调节、支持与生产功能；低山丘陵区生态系统服务以生产功能为主，兼顾其他功能，但要注重各功能的协调。

1 研究区概况

太行山区是华北平原的天然屏障和重要的水源地，位于34°36′N～40°47′N，110°42′E～116°35′E，包括北京、河南、山西和河北的101个县，面积约为1 400万hm2(1 hm2=10 000 m2)。山势东坡陡、西坡缓，北高南低，最高海拔约3 000 m，最低海拔约为10 m，平均海拔为860 m(图1)。太行山区属于暖温带半湿润大陆性季风气候，多年年平均降雨量为570 mm，年均温为10℃。太行山区主要土地利用类型为耕地、林地、草地、建筑用地，占总面积的90 %以上。

图1 太行山区地理位置图

2 数据处理

太行山区生态系统垂直分异关键线带的确定选择土地利用、植被NPP、产水量和气温作为太行山区生态系统要素土、生、水和气的表征性指标。不同的土壤会有不同的土地利用类型[17]，不同的土地利用类型会产生不同的生态过程，进而反应不同的生态系统结构[18]。生态系统要素包括生态系统中的一切生物。而植被是其中最重要的一种生物。植被NPP表征生态系统的质量状况和生产能力[19]，是评价生态系统结构和功能的重要指标[10-11]。山区素有水塔之称[20]，太行山区为山区和京津冀地区提供淡水。生态系统服务功能能够反应生态系统结构[21]。温度变化影响生态系统的结构、功能和生产力[22]，不同的气温条件下具有不同的生态系统结构。

综上所述，土地利用、植被NPP、产水量和气温作为太行山区生态系统要素土、生、水和气的表征性指标，具有其科学性和合理性。

2.1 土地利用

太行山区土地利用数据是利用30 m分辨率的TM影像(2010年)，通过目视判读而得，经过实地考察证明其分类精度在80 %以上。土地资源分类系统依据国家二调分类标准而定，包括6个一级类型(耕地、林地、草地、水域、城乡居民点和工矿用地、未利用土地)，25个二级类型(水田、旱作；有林地、灌木林地、疏林地，其他林地；高覆盖草地、中覆盖草地、低覆盖草地；河渠、湖泊、水库坑塘、冰川积雪、海涂、滩地；城镇用地、农村居民点、工交建设用地；沙地、戈壁、盐碱地、沼泽地、裸土地、裸岩、其他)。在此基础上对土地利用类型进行了重新分类，分为7个类型(农田、森林、草地、建设用地、湿地、荒漠、水域)。

2.2 植被NPP

植被NPP数据空间分辨率为1 km(https://ipdaac.usgs.gov)。该产品基于MODIS/TERRA卫星遥感数据，通过BIOME-BGC模型计算得到植被NPP数据，其空间分辨率为1 km。用一元线性回归分析法，基于像元分析太行山区2000—2014年植被NPP变化，线性回归的斜率k即为NPP变化速率，公式为：

式中：k为NPP变化速率；N=15为年数；NPPi为第i年的NPP值，t为年份。基于ArcGIS的栅格计算功能，可得到植被NPP在2000—2014年间的变化趋势图。当k＞0，说明植被NPP呈增加趋势，反之为减少趋势。

2.3 产水量

基于inVEST模型中产水量子模块[23-24]原理，本研究中的产水量为降雨量减去实际蒸散发量。降雨数据是由太行山区102个多年平均降水气象站数据为基础进行插值，插值方法采用(普通)克里金算法，生成栅格分辨率1 km。太行山区实际蒸散发是以MOD16蒸散发为数据源，空间分辨率为1 km(https://ntsg.umt.edu)，该产品基于卫星遥感数据，通过Penman-Monteith公式[25]计算得到实际蒸散发数据。

2.4 气温

太行山区年均气温数据由太行山区102个多年平均气温气象站数据为基础进行插值，插值方法采用(普通)克里金算法，生成栅格分辨率1 km。

2.5 高程

高程数据来源于ASTER GDEM(advanced spaceborne theemal emission and reflection radiometer global digital elevation model)，即先进星载热发射和反射辐射仪全球数字高程模型(http://gdem.ersdac.jspacesystems.or.jp)，其空间分辨率为30 m(图1)。

3 研究方法

3.1 关键线带确定

利用ArcGIS软件将土地利用、植被NPP、产水量和气温4个指标数据在垂直梯度上以100 m为一个统计单元分别进行统计。由于高海拔地区的面积比例较少，其中的土地利用和植被NPP以及产水量指标将海拔大于2 500 m的区域作为一个统计单元。

在确定各个生态系统要素表征指标垂直分异的关键线的基础上，采用加权的方法综合各个关键线，确定太行山区生态系统垂直分异的关键线。由于生态系统要素表征指标对于生态系统的区划都具有重要的意义[26-27]，定义每个指标的权重值为1/3(由于年均气温指标没有明显的垂直分异特征，不能确定其垂直分异的关键线)。

3.2 生态系统服务评估

计算生态系统服务的公式如下：

其中，AK是某种生态系统类型面积，VCK是某种生态系统类型单位面积生态系统服务。VCK是一个当量的价值量乘以某种当量，本研究一个当量的价值为3 406.5元。

3.3 当量表修正

根据谢高地等[11]单位面积生态系统服务当量表，对太行山区的生态系统服务进行计算。由于其当量表是针对全国的一个平均值，为了更加准确地计算研究区域的生态系统服务，本文利用NPP(MODIS)对其进行修正。其修正系数公式如下：CCi=NPPi/NPPeve，其中，CCi是当量表的修正系数，NPPi是研究区域某种土地利用类型的净初级生产力的平均值，NPPeve是全国相对应的土地利用类型的净初级生产力的平均值。

3.4 生态系统服务功能格局指数

太行山区各个分区的面积有较大的差异，因此生态系统服务的总价值不能说明太行山区各分区生态系统服务功能格局。另外，在谢高地等[11]的当量表中，单位面积的生态系统服务绝对价值存在调节服务＞支持服务＞供给服务＞文化服务的趋势，因此，各个分区的生态系统服务的主体功能都是调节服务，这与太行山区实地情况不符。为了确定各分区生态系统服务功能格局，这里引入了功能格局指数，该指数以太行山区每种单位面积生态系统价值为基准，与各分区进行比较，确定其分区的主体服务功能，公式为：FPD=ESVij/ESVej，其中i=1、2、3，分别代表低山丘陵区、中山区和亚高山区，j=(1、2、3、4)，分别代表供给服务、调节服务、支持服务和文化服务。FPD为功能格局指数，ESVij为太行山区i分区j服务单位面积价值，ESVej为太行山区j服务单位面积价值，服务功能指数越大说明该项服务功能在分区中越重要。

4 结果与分析

4.1 太行山区生态系统垂直分异特征

(1)土地利用的垂直分异特征

以海拔高度和土地利用中农田、森林、草地、建设用地面积百分比作为聚类指标，通过聚类分析将所有指标分为3类，高程＜600 m，高程600～1 600 m和高程＞1600 m，即600 m、1 600 m为太行山区土地利用垂直分异的突变线。太行山区各类分区土地利用类型主要以农田、森林、草地、建设用地为主，3个分区随着海拔的升高，农田、建设用地面积占所在分区总面积的比例逐渐减少，森林、草地的面积比例逐渐增加。

(2)植被NPP垂直分异特征

由图2可知，太行山区垂直分布上NPP年际间变化速率有两个明显的峰值。利用回归分析将NPP年际间变化速率分段进行拟合，二次多项式拟合较好(R2＞0.8)，拟合曲线分别在海拔约600 m、1 300 m处出现峰值，表明在两个海拔处NPP年际间变化速率最大。NPP表征生态系统的质量状况和生产能力[19]，是评价生态系统结构和功能的重要指标[10-11]。600 m、1 300 m处NPP年际间变化速率大说明该海拔处生态系统结构和功能出现分异的特征。

(3)产水量垂直分异特征

由图3可知太行山区垂直分布上产水量有两个峰值。利用回归分析将产水量分段进行拟合，二次多项式拟合较好(R2≈0.8)，拟合曲线分别在海拔约500 m、1 400 m处出现峰值。太行山区每公顷产水量约500 m3，产水量随着海拔的升高呈现减少的趋势，且呈现波动性的变化。太行山区垂直分布上产水量分别在海拔约500 m、1 400 m处出现峰值，说明对于提供淡水这项服务功能，生态系统结构在这两处达到最佳状态，其余的海拔高度处生态系统结构处于可优化状态。因此500 m、1 400 m为产水量垂直分异的关键线。

图2 太行山区植被NPP年际间变化速率

图3 太行山区产水量

(4)年均气温

太行山区多年年均气温为8.99 ℃。随着海拔的升高年均气温稳定降低(图4)。年均气温在垂直梯度上没有明显的分异特征。但是，气温是生态系统不可缺少的要素，在太行山区生态系统要素垂直分异的关键线/带的确定中不可或缺。

4.2 太行山区生态系统垂直分异关键线带确定

每个生态系统要素表征指标都有其垂直分异的关键线，且关键线有一定的差异。利用加权方法综合各个要素指标的关键线，得到太行山区生态系统垂直分异的两个关键线，即567 m、1 433 m。由于多数指标都是以100 m为统计单元的，所以其关键线也应该以100 m为单元，即太行山区生态系统垂直分异的两个关键带为(567±50) m和(1433±50) m，即500～600 m和1 400～1 500 m。两个关键带的面积约占太行山区总面积的9.07 %，其中，500～600 m关键带的面积比例为4.85 %，主要位于太行山区东部，跨越太行山区南北；1 400～1 500 m关键带的面积比例为4.22 %，零散分布于太行山区中西部 (图5)。

图4 太行山区年均气温

图5 太行山区生态系统垂直分异的关键带空间分布

4.3 太行山区生态系统服务垂直分异特征

4.3.1 太行山区垂直梯度生态系统服务评估

太行山区各个分区，即低山丘陵区、中山区、亚高山区，1公顷单位面积服务总价值分别为19 788元、21 706元、82 502元，4种生态系统服务单位面积服务价值也是亚高山区＞中山区＞低山丘陵区(图6)，这是因为在单位面积服务价值高的主要土地利用类型(除了水域、湿地)，如林地和草地，随着海拔的升高其面积比例增加。由于在价值当量表中所有的土地利用类型的4种生态系统服务当量都是调节服务＞支持服务＞供给服务＞文化服务，因此，在3个分区单位面积的4种生态系统服务具有相同的趋势，即调节服务＞支持服务＞供给服务＞文化服务。

4.3.2 太行山区垂直梯度生态系统服务功能格局

表1列出了太行山各个分区的服务功能指数。在低山丘陵区，供给服务的服务功能指数最高，为主要服务功能，约是其他服务功能指数的2倍。在中山区，4种服务功能指数都接近于0.5，说明4种服务功能基本是均衡的。在高山区，除了供给服务功能格局指数小于2，其余3种服务功能指数都大于2，说明除了供给功能，其余3种服务功能都是高山区的主体功能。

综上所述，太行山区垂直梯度生态系统服务功能格局为：随着海拔的升高，以供给服务功能为主转变为4种生态系统服务功能均衡，再由4种生态系统服务功能均衡转变为调节、支持和文化服务功能为主。

图6 太行山区各分区单位面积生态系统服务

表1 服务功能指数

4.4 基于可持续发展的太行山区生态系统分类管理模式

太行山低山丘陵区，居住人口约占太行山区总人口的66 %，耕地资源在该区域较多。山区耕地所带来的负生态系统服务也是显而易见的，如耕地降低土壤的有机质含量[28]、造成水污染、减少生物多样性[29]、消耗地下水资源[30]。太行山低山丘陵区存在人多地少的矛盾，为了实现生存的需要，必须实现山坡地高产出，同时还要保护环境，防止耕地对环境造成的负服务。农林复合系统是在综合考虑社会、经济和生态因素的前提下，将乔木和灌木等有机结合到农牧生态系统中，为社会提供供给功能的同时保障自然资源的可持续生产[31]。太行山区低山丘陵区应该发展农林复合生态系统，找到农林最佳组合模式，即保证供给服务的同时，又确保生态系统可持续发展。

太行山中山区，农村土地承包经营权流转逐渐成为农村经济发展中的热点问题。在太行山中山区经过土地流转发展了一批因地制宜的开发模式，如井陉县的万亩连翘基地、平山县葫芦峪的现代农业产业园、武安市农林牧畜相结合的模式。在开发利用初期，土壤发生扰动，土壤微生物群落发生改变[32]，土地利用变化会改变区域的水文过程[33]，造林会降低土壤的总磷含量[34]，伐林会降低土壤有机质的含量[35]；植被覆盖率较低，水土流失较严重[36]。但随着开发利用模式形成，土壤固定，植被覆盖率增加，在初期出现的负生态系统服务消失，并且在不同程度上增加了供给服务，促进经济发展的同时保护了环境。

太行山亚高山区，土地利用类型主要以森林和草地为主，森林和草地的主要功能为调节气候、水文，保持土壤、维持养分循环[37-39]。为保护这些森林和草地免受人类活动破坏，高山区建立了多个自然保护区，同时适当地发展旅游业，增加了文化服务。

5 结论

山区可持续发展是山区研究的重要关注点，生态系统服务优化调控管理是实现山区可持续发展的根本措施。太行山区因海拔高度较低，虽然在垂直梯度上没有表现出明显的垂直带谱特征，但仍呈现较明显的生态特征垂直分异规律。

选取太行山区产水量、土地利用、温度和植被NPP作为生态系统主要要素的表征指标，利用ArcGIS软件以100 m海拔梯度为一个单位分别研究其垂直梯度变化特征，并综合4个指标确定了太行山区生态系统500～600 m和1 400～1 500 m两个垂直分异关键线带，这两个线带将太行山区划分为＞1 500 m的亚高山区、600～1 400 m的中山区和＜500 m的低山丘陵区3类生态系统空间格局。

在确定了“二带三区”的基础上，分析得到3个生态类型区单位面积生态系统服务由大到小的梯度变化依次为亚高山区＞中山区＞低山丘陵区；分析3个生态类型区生态系统服务功能指数后表明：低山丘陵区供给服务最高；中山区供给、调节、支持和文化4种服务基本相同；亚高山区除供给服务较低外其余3种服务均较高。由此说明：随着海拔升高，太行山区垂直梯度生态系统服务由以供给服务为主转变为4种服务均衡，再转变为以非供给服务为主。

太行山区3个生态类型区生态系统服务梯度调控管理的可持续发展模式和目标为：亚高山区以非供给服务为主，调控功能目标为生态系统服务外溢区；中山区以4种功能均衡发展，调控功能目标为生态系统服务外溢廓道和供给服务人工调控提升区；低山丘陵区以供给服务为主导，兼顾协调其他功能，调控功能目标为生态系统服务过程控制与调整、水土耦合与农林结构优化配置区。

(2018年1月9日收稿)

[1] 邓伟, 唐伟. 试论中国山区城镇化方向选择及对策[J]. 山地学报,2013, 31(2): 168-173.

[2] 邓伟, 熊永兰, 赵纪东, 等. 国际山地研究计划的启示[J]. 山地学报,2013, 31(3): 377-384.

[3] KHAN S M, PAGE S E, AHMAD H, et al. Sustainable utilization and conservation of plant biodiversity in montane ecosystems: the western Himalayas as a case study [J]. Annals of Botany, 2013, 112(3): 479-501.

[4] MUNTEANU C, KUEMMERLE T, BOLTIZIAR M, et al. Forest and agricultural land change in the Carpathian region—a meta-analysis of long-term patterns and drivers of change [J]. Land Use Policy, 2014,38: 685-697.

[5] BORSDORF A，BRAUN V．The European and Global dimension of mountain research [EB/OL]. 2009 [2018-01-09]. http://rga.revues.org/index630.

[6] BRUNNER S H, GRÊT-REGAMEY A. Policy strategies to foster the resilience of mountain social-ecological systems under uncertain global change [J]. Environmental Science & Policy, 2016, 66: 129-139.

[7] IMPERIALE A J, VANCLAY F. Using social impact assessment to strengthen community resilience in sustainable rural development in mountain areas [J]. Mountain Research and Development, 2016, 36(4):431-442.

[8] 赵东娟, 齐伟, 赵胜亭, 等. 基于GIS的山区县域土地利用格局优化研究[J]. 农业工程学报, 2008, 24(2): 101-106.

[9] REGER B, GÖTTLEIN A, KATZENSTEINER K, et al. Assessing the sensitivity of mountain forests to site degradation in the northern limestone alps, Europe [J]. Mountain Research and Development,2015, 35(2): 139-151.

[10] COSTANZA R, D’ARGE R, DE GROOT R S, et al. The value of the world’s ecosystem services and natural capital [J]. Nature, 1997, 387:253-260.

[11] 谢高地, 张彩霞, 张雷明, 等. 基于单位面积价值当量因子的生态系统服务价值化方法改进[J]. 自然资源学报, 2015, 30(8): 1243-1254.

[12] HÄYHÄ T, FRANZESE P P, PALETTO A, et al. Assessing, valuing,and mapping ecosystem services in Alpine forests [J]. Ecosystem Services, 2015, 14: 12-23.

[13] 伍光和. 自然地理学[M]. 北京: 高等教育出版社, 2000.

[14] 李国伟, 赵伟, 魏亚伟, 等. 天然林资源保护工程对长白山林区森林生态系统服务功能的影响[J]. 生态学报. 2015, 35(4): 984-992.

[15] 侯学煜. 论中国各植被区的山地垂直植被带谱的特征, 中国植物学会,中国植物学会三十周年年会论文摘要汇编[C]. 北京：[出版者不详], 1963: 254-258.

[16] 姚永慧, 张百平, 韩芳, 等. 横断山区垂直带谱的分布模式与坡向效应[J].山地学报, 2010, 28(1): 11-20.

[17] 傅伯杰, 张立伟. 土地利用变化与生态系统服务：概念、方法与进展[J]. 地理科学进展, 2014, 33(4): 441-446.

[18] FU B J, CHEN L D, MA K M, et al. The relationships between land use and soil conditions in the hilly area of the Loess Plateau in northern Shaanxi, China [J]. Catena, 2000, 39(1): 69-78.

[19] WANG J, MENG J J, CAI Y L. Assessing vegetation dynamics impacted by climate change in the southwestern karst region of China with AVHRR NDVI and AVHRR NPP time-series [J]. Environmental Geology, 2008, 54: 1185-1195.

[20] CALDWELL P V, MINIAT C F, ELLIOTT K J, et al. Declining water yield from forested mountain watersheds in response to climate change and forest mesophication [J]. Global Change Biology, 2016, 22(9):2997-3012.

[21] CHUAI X W, HUANG X J, WU C Y, et al. Land use and ecosystems services value changes and ecological land management in coastal Jiangsu, China [J]. Habitat International, 2016, 57: 164-174.

[22] 初小静, 韩广轩. 气温和降雨量对中国湿地生态系统CO2交换的影响[J]. 应用生态学 报, 2015, 26 (10 ): 2978 -2990.

[23] REDHEAD J W, STRATFORD C, SHARPS K, et al. Empirical validation of the InVEST water yield ecosystem service model at a national scale [J]. Science of the Total Environment, 2016, 569: 1418-1426.

[24] SHARPS K, MASANTE D, THOMAS A, et al. Comparing strengths and weaknesses of three ecosystem services modelling tools in a diverse UK river catchment [J]. Science of the Total Environment,2017, 584/585: 118-130.

[25] MONTEITH J L. Evaporation and environment [J]. Symposia of the Society for Experimental Biology, 1965, 19: 205-234.

[26] 郑度, 葛全胜, 张雪芹, 等. 中国区划工作的回顾与展望[J]. 地理研究, 2005, 24(3): 330-344.

[27] 傅伯杰, 刘世梁, 马克明. 生态系统综合评价的内容与方法[J]. 生态学报, 2001, 21(11): 1886-1892

[28] GASPAR Q L, FERRER L G, LZQUIERDO N A. Lateral and depth patterns of soil organic carbon fractions in a mountain Mediterranean agrosystem [J]. Journal of Agricultural Science, 2016, 154(2): 287-304.

[29] SWINTON S M, LUPI F, ROBERTSON G P, et al. Ecosystem services and agriculture: cultivating agricultural ecosystems for diverse benefits[J]. Ecological Economics, 2007, 64(2): 245-252.

[30] FOSTER S S D, PERRY C J. Improving groundwater resource accounting in irrigated areas: a prerequisite for promoting sustainable use [J]. Hydrogeology Journal, 2010, 18: 291-294.

[31] LAUDGREN. The defination of agroforestry [J]. Agroforestry System.1982, 1(1): 4-7.

[32] FRENCH K E, TKACZ A, TURNBULL L A. Conversion of grassland to arable decreases microbial diversity and alters community composition [J]. Applied Soil Ecology, 2017, 110: 43-52.

[33] WOLDESENBET T A, ELAGIB N A, RIBBE L, et al. Hydrological responses to land use/cover changes in the source region of the Upper Blue Nile Basin, Ethiopia [J]. Science of the Total Environment, 2017,575: 724-741.

[34] DENG Q, MCMAHON D E, XIANG Y, et al. A global meta-analysis of soil phosphorus dynamics after afforestation [J]. New Phytologist,2017, 213(1): 181-192.

[35] VILLARINO S H, STUDDERT G A, BALDASSINI P, et al.Deforestation impacts on soil organic carbon stocks in the Semiarid Chaco Region, Argentina [J]. Science of the Total Environment, 2017,575: 1056-1065.

[36] PARDINI G, GISPERT M, EMRAN M, et al. Rainfall/runoff/erosion relationships and soil properties survey in abandoned shallow soils of NE Spain [J]. Journal of Soils and Sediments, 2017, 17(2): 499-514.

[37] BUMA B, WESSMAN C A. Forest resilience, climate change, and opportunities for adaptation: A specific case of a general problem [J].Forest Ecology and Management, 2013, 306: 216-225.

[38] GARCIA-NIETO A P, GARCIA-LLORENTE M, INIESTAARANDIA I, et al. Mapping forest ecosystem services: From providing units to beneficiaries [J]. Ecosystem Services, 2013, 4: 126-138.

[39] WINOWIECKI L, VAGEN T G, MASSAWE B, et al. Landscapescale variability of soil health indicators: effects of cultivation on soil organic carbon in the Usambara Mountains of Tanzania [J]. Nutrient Cycling in Agroecosystems, 2016, 105(3): 263-274.