2000—2020年祁连山生态系统服务时空分异研究

赵文鹏， 吕荣芳， 庞吉丽， 张建明， 王乃昂

（兰州大学 资源环境学院，甘肃 兰州 730000）

0 引言

生态系统服务（Ecosystem Services， ES）是人类为了满足自身需求而直接或间接地从生态系统中获得的各项产品和服务［1］。它是连接自然生态过程和人类活动的纽带，其对自然资源的合理利用和配置，是实现生态环境保护和区域可持续发展的前提和基础［2-3］。生态系统服务的形成依赖于一定的时间和空间尺度上的生态系统结构和生态过程，在不同的时空尺度上其表现出特定的主导作用和效果，从而密切影响着人类从生态系统中所获利益大小［4-6］。不同区域之间的地形地貌特征［7-8］、气候变化［9］、人类活动［10-11］以及土地覆盖变化［12］等因素会影响生态系统结构及其空间特征，进而导致了生态系统服务中的时空分布与差异［5］。定量评估与全面分析特定时空内的各项生态系统服务的时空分异，能够促进自然资源的保护管理和针对性的生态恢复，对区域生态系统的科学管理与可持续发展具有重要意义［13-15］。

祁连山以其独特地理环境成为全球典型生态环境脆弱区，是我国西北干旱区、半干旱区重要的生态屏障，被誉为河西走廊的“母亲山”［16］，对于祁连山地区的生态系统服务研究也一直是学术界关注的重点内容。赵亚茹等［17］对石羊河上游产水量进行定量评估与驱动因素识别；刘洋等［18］探究疏勒河流域碳储量时空变化与土地覆盖之间的关系；Wei等［19］评估了气候与土地覆盖变化对疏勒河流域产水量的影响；卿苗等［20］等利用InVEST 模型和FLUS模型研究了石羊河流域在过去40年间和未来不同情景下土地覆盖变化对碳储量的影响；李芳等［21］基于InVEST 模型和时空地理加权回归方法探索了黑河流域上游产水服务时空演变特征；Xu等［22］基于InVEST 模型定量计算了2000—2019年祁连山自然保护区的水源涵养量并分析了其变化趋势的原因。上述研究多在祁连山地区内的局部尺度上对单项生态系统服务进行多年的定量计算，并探索其时空变化。陈童尧等［23］对2015 年祁连山自然保护区土壤保持服务在不同海拔、坡度和土地覆盖等进行定量评估；杜可心等［24］基于定量计算的2018年的多项生态系统服务，并结合海拔、坡度、地形起伏度和地形位指数因子分析其梯度效应；吕荣芳等［2］基于2019 年的六项生态系统服务，从社会-生态环境角度分析祁连山地区各项生态系统服务间权衡与协同关系，这些研究多基于一年计算的生态系统服务在局部尺度上对其梯度效应以及权衡协同关系等展开分析。总体来看，祁连山地区生态系统服务评估与分析的研究多基于流域、自然保护区等局部尺度，在祁连山地区全域进行多年尺度多项服务的时空分异性探索有所不足，不同的时空尺度在一定程度上会影响其生态系统管理策略的制定［6］，该区域复杂的生态环境也亟须从土地覆盖以及地形梯度等方面对其生态系统服务变化进行全面的分析［25］。

因此，本文基于祁连山2000—2020年5期的30 m分辨率的土地覆盖、气象等数据，利用InVEST 模型和CASA 模型评估祁连山产水量、碳储存、土壤保持、水质净化、碳固定和生境质量六项生态系统服务的时空分异特征，利用分区统计得到六项生态系统服务的土地覆盖以及梯度分布，识别其热点区域，以期能够深入了解祁连山地区六项生态系统服务分布的空间特征以及时间演变方向，探索不同土地覆盖与海拔梯度之间生态系统供给能力的差异性，从而为地理环境复杂多样的祁连山地区的生态管理与可持续发展提供科学指导和决策依据。

1 研究区概况

祁连山地区（93°56′～103°90′ E， 35°84′～39°97′ N）地处青藏、蒙新、黄土三大高原的交会处，内有黑河、疏勒河、石羊河三大内陆河流域，是中国西北干旱半干旱地区最重要的水源地之一，也是中国西北地区重要的生态安全屏障（图1）。祁连山地区总面积约为19.47×104km2，主要由一系列东南—西北走向的高山与谷底组成，存在着多条现代冰川，海拔在1 965～5 792 m 之间，其中约有80%的区域海拔高于3 000 m。该区域气候条件差异较大，属于典型的山地半干旱气候，气候特征以寒、旱为主，且自东向西逐渐变冷变干。年均降水量约为300～700 mm，且降水总量的80%集中于夏秋两季［19］。年均气温较低，约为-2.2～2.3 ℃，气温年较差、日较差较大，降水和气温随海拔高度变化呈现出明显的垂直梯度。此外，祁连山地区太阳辐射强度较大，大部分地区一年日照时数大于2 800 h。这种复杂的自然条件也造就了该区域不同地区显著的水热条件差异，在垂直与水平方向上均形成了多样的植被与土壤类型，土地覆盖类型复杂多样，主要为草原、草甸、裸地和森林［26］。

图1 祁连山地区的位置和海拔Fig. 1 The location of Qilian Mountains in China and its elevation

2 数据来源与研究方法

2.1 数据来源

本文所涉及的高程海拔数据、土壤类型均来自于地理空间数据云（http：//www.gscloud.cn/）；人口数据来源于中国科学院资源环境科学数据中心（https：//www.resdc.cn/）；2000年、2005年、2010年、2015 年、2020 年5 期土地覆盖数据来源于国家青藏高原数据中心（http：//data. tpdc. ac. cn/zh-hans/）与中国科学院资源环境科学与数据中心（http：//www.resdc.cn）［27-29］；气象站点数据下载于中国气象科学数据共享服务网的《中国地面气候资料日值数据集》（http：//data. cma. cn/），共45 个站点，并利用ANUSPLIN 模型进行空间插值［30］；水文站点河流年径流量数据来自于《甘肃省水资源公报》和《青海省水资源公报》，共9个站点。上述数据分辨率均统一为30 m。

2.2 生态系统服务计算

本文基于InVEST 模型与CASA 模型等生态系统服务定量评估模型，估算祁连山地区产水量、土壤保持、水质净化、碳固定、生境质量与碳储存服务时空分异特征。产水量服务基于InVEST 模型的Water Yield 模块进行计算，该模块基于水量平衡原理并利用Budyko曲线根据降水、植物蒸腾和地表蒸发等计算所得，相关参数通过研究区文献进行修正［19，31-32］，基于水文站点实测数据进一步对模型结果进行分流域校准，总体精度达到94%以上（表1）；土壤保持服务基于InVEST 模型的SDR 模块进行计算，模型中所需降雨侵蚀因子（R）运用Wischmeier基于月尺度的经验模型计算所得［33］，植被覆盖因子（C）参考蔡崇法进行计算［34］，其余参数则参考研究区相关文献［23］获得；水质净化服务基于InVEST模型的NDR 模块进行评估，模型仅考虑面源污染，根据相关研究所获取的污染物拦截系数与过滤量以及相关参数［35-37］，对植被与土壤对径流中的N、P营养盐保持量进行估算；碳固定服务基于InVEST模型的Carbon 模块进行评估，碳库数据参考研究区相关文献［2，18，38-42］所得；生境质量服务基于In-VEST 模型的Habitat Quality 模块，由相关文献［43-44］以及结合研究区实际情况确定不同土地覆盖的威胁因子敏感度、威胁因子的空间权重与影响距离与外界威胁强度；碳储存则利用CASA 模型基于生态系统碳循环过程估算植被净初级生产力，具体模型原理及参数确定详见文献［45］。以上生态系统服务计算中相关生物物理参数详见表2与表3。

表1 2000—2020年祁连山地区子流域水文站点年径流量与模拟年产水量对比Table 1 Comparison of annual runoff and simulated annual water production at hydrological stations in the Qilian Mountains sub-basin from 2000 to 2020

表2 产水量、碳固定、土壤保持与生境敏感性因子生物物理参数表Table 2 Table of biophysical parameters of water yield， carbon sequestration， soil retention and habitat sensitivity factors

表3 威胁因子参数表Table 3 Threat factor parameter

2.3 地形数据处理

（1） 地形起伏度

地形起伏度是基于ArcGIS Pro 的“邻域分析工具”进行计算，其计算公式为：

式中：D为地形起伏度；Cmax为窗口内的海拔最高值；Cmin为窗口内的海拔最低值。

（2） 地形位指数

地形位指数是综合描述空间内任一点海拔和坡度属性信息的指标［46］，其计算公式：

式中：T为地形位指数；E与E0分别代表研究内任意海拔与平均海拔；S与S0分别代表研究区内任一坡度与平均坡度。坡度大、海拔高的地区地形位指数越大，反之越小。

（3） 地形数据分级标准

为分析不同地形梯度中生态系统服务分布特点，采用分位数法［47］将坡度、海拔、地形起伏度与地形位指数分为6级（表4）。

表4 海拔、坡度、地形起伏度和地形位指数分级标准Table 4 Classification of elevation， slope， landform relief gradients and terrain niche gradients

2.4 热点分析方法

生态系统服务热点区往往是具有高服务多样性、高供给能力和高生物物理等优势的区域，能够进一步体现生态系统服务的空间供给能力强弱［4］。本研究将各生态系统服务大于其各自平均值的栅格视为该服务的热点区域，并利用ArcGIS Pro“叠加分析工具”，对祁连山地区多重生态系统服务热点区进行识别。根据栅格单元中热点区的频次，将各项服务供给值均未超过各自平均值的栅格定义为非热点区，若仅有一种服务供给值超过其平均值定义为Ⅰ类热点区，以此类推分别定义为Ⅱ类、Ⅲ类、Ⅳ类和Ⅴ类热点区。

3 结果与分析

3.1 2000—2020年祁连山地区土地覆盖变化

2000—2020 年祁连山地区各土地覆盖类型变化总体较为稳定，但不同土地覆盖之间的转入转出较为复杂（图2）。草地与裸地分布广泛，约占总面积的86%，各土地覆盖之间的流动主要集中在冰川积雪、裸地与草地之间。其中，在2015 年之前冰川积雪转入裸地与草地一直占据主导地位；2015 年至2020 年之间，裸地大量转入冰川积雪，导致冰川积雪用地面积百分比增加。祁连山地区的水热条件变化，也在一定程度上影响了土地覆盖的转入与转出。在这21年间，祁连山地区的气温和降水均呈上升趋势，平均气温与平均降水量分别为-2.1 ℃与406.6 mm。气温倾向率为0.18 ℃·（5a）-1，而降水倾向率为32 mm·（5a）-1，气候逐渐向暖湿转型［48］，这也导致冰川积雪面积的减少与草地、裸地面积的增加，总体呈现植被面积占比的波动增加趋势。

图2 2000—2020年祁连山土地覆盖变化流动桑基图Fig.2 Sankey diagram of land use changes in Qilian Mountains from 2000 to 2020

3.2 祁连山地区生态系统服务时空分异特征

在时间上，2000—2020 年祁连山地区的产水量、碳固定、土壤保持、水质净化、碳储存和生境质量六项生态系统服务供给量均呈上升趋势，增量较高的区域主要分布在中西部的山区（图3）。其中，产水量在研究时段内平均增长速率为0.98 mm·hm-2·a-1，西部山区为该服务增长高值区；碳储存的单位面积供给量由2000 年的118.2 gC·m-2增加到2020 年的140.5 gC·m-2；土壤保持单位面积供给量增加显著，增幅达111%；水质净化服务则整体变化不大，平均增幅为3%左右；碳固定在5 个时期的单位面积供给量分别为144.7、146.6、150.1、150.4、158.5 t·hm-2，呈波动增长趋势；生境质量呈现稳步增长趋势，平均值增幅为1%，且西部增长幅度要高于东部。

图3 2000—2020年祁连山生态系统服务空间演变（从左至右分别为2000年、2020年以及2000—2020年均值）Fig. 3 Spatial evolution of ESs in Qilian Mountains， 2000—2020 （from left to right：the average values in 2000， 2020 and 2000—2020， respectively）

在空间上，各项服务供给量年均值存在明显的空间差异化规律（图3），供给量高值区主要集中在气候更为暖湿的东部地区，总体呈现“东高西低”的分布格局。产水量的空间分布海拔差异明显，供给量高值区集中在研究区内山区，总体上自东部至西部产水量逐渐降低，在研究区东南部高值面积出现了小幅退缩。土壤保持供给量高值区分布在北部、东部以及南部边缘的山区，总体呈现东北向西南减弱趋势。碳储存、碳固定、水质净化以及生境质量空间分布格局较为相似，其供给量高值区集中在祁连山东部地势平坦地带，低值区则多分布在西部裸地、山地冰川与柴达木盆地西北边缘地区。

3.3 基于土地覆盖的祁连山地区生态系统服务及变化

如图4 所示，不同土地覆盖类型上的生态系统服务供给量具有明显差异，随着时间增加，各土地覆盖的生态系统服务单位面积供给能力总体上呈增加趋势。其中，林地、草地和灌丛既是研究区主要分布土地覆盖，也是生态系统服务供给量和单位面积生态系统服务供给量最高的土地覆盖类型，其单位面积供给量约占所有土地覆盖单位面积总供给量的60%以上。其中，耕地与建设用地的单位面积的水质净化供给量最大；林地、灌丛与草地的单位面积碳固定供给量最大，三种土地覆盖的单位面积碳固定供给量约占全部土地覆盖的86%。单位面积碳储存与生境质量供给量年度差异不明显，其高值区均主要集中在林地、草地与湿地上。各土地覆盖的单位面积产水量、碳储存以及土壤保持供给量总体呈波动增长趋势，林地、草地的单位面积供给量较高；由于冰川积雪和裸地的蒸散发能力较弱，其单位面积产水量供给量较高，二者约占各土地覆盖单位面积总供给量的50%。

图4 2000—2020年基于土地覆盖的祁连山地区生态系统服务均值及变化Fig. 4 Mean and trend of ESs in Qilian Mountains based on land class from 2000 to 2020： water yield （a）， carbon sequestration （b）， soil retention （c）， nutrient retention （d）， carbon storage （e） and habitat quality （f）

3.4 基于地形梯度的祁连山地区生态系统服务分布

对2000—2020 年各生态系统服务的年均值进行梯度分析之前，为了消除不同服务均值之间的量级影响，对各服务均值进行z-score 标准化处理（图5）。除海拔外，其余各项地形梯度上各生态系统服务均值的分布基本相似。其中，产水量在不同的地形梯度上都呈增加趋势，而碳固定、碳储存与生境质量在各个梯度上均呈现着近似的分布特点，除在海拔上具有先降低后增加再降低的起伏波动特点外，均在4级梯度处达到最高值。土壤保持与产水量在坡度、地面起伏度以及地形位指数上呈相似的增加趋势，而在海拔中却与碳固定等服务相似，呈现先减后增再减的趋势，在1 级、3 级和4 级梯度处达到高值。水质净化服务总体呈现随梯度等级增加而服务降低的趋势，2 级梯度处是其在海拔、坡度以及地面起伏度上的转折点，其高值多分布在海拔较低且坡度较小的地区，也是人类活动较为频繁的区域。

图5 祁连山地区生态系统服务在各地形梯度上的分布特点Fig. 5 Distributions of ESs in Qilian Mountains in topographic gradient zones：elevation （a）， slope （b）， landform relief （c） and terrain niche （d）

3.5 热点区分析

如图6 所示，2000—2020 年祁连山生态系统服务热点区时空变化整体较为稳定。在时间变化上，各热点区平均每年变化面积约占研究区面积的0.5%，Ⅲ类和Ⅳ类热点区面积呈增加趋势，而其余热点区面积呈减少趋势，其中Ⅲ类热点区增加最为显著，平均每年增加2 184 km2，非热点区减少最显著，平均每年减少1 709.5 km2，Ⅰ类、Ⅱ类和Ⅴ类热点区则平均每年减少523.2 km2。在空间分布上，研究区的热点区呈现自东向西逐渐由强变弱的特征趋势，非热点区和Ⅲ类热点区分布范围最广，平均约占研究区总面积的22%，Ⅴ类热点区分布范围最小，多集中在东部地区的林地、灌丛和草地中。研究区中西部热点区增加明显，以Ⅰ类和Ⅲ类热点区增加最为明显。总体来看，研究区整体环境趋于改善，但在局部区域有退化趋势，因此需要加强Ⅰ类、Ⅱ类和Ⅴ类热点区的管理，预防其生态环境退化，同时推动Ⅲ类热点区向Ⅳ类热点区转化，实现区域生态系统服务供给的协同最优。

4 讨论

生态系统服务的时空分异特征往往是自然和人类活动因素共同驱动的结果［2］，总体来看，2000—2020 年祁连山地区六项目标生态系统服务的供给能力呈明显增加趋势，各项服务之间的协同程度有所提升，也在一定程度上说明该区域的生态环境在逐渐改善和提高，但在部分地区仍旧存在着脆弱甚至倒退的生态环境状况。研究区的中部和西南部目标生态系统服务供给量提升明显，而东南部和西北部等人类活动相对频繁的河谷地带、耕作区的生态脆弱性特征明显，部分服务供给量甚至出现了下降。祁连山地区生态系统服务的时空分异性与气候、人类活动、地形、土地覆盖以及生态措施等密切相关［49］。

祁连山地区的暖湿化倾向与人类活动强度的增加是造成目标服务发生时空变化的最主要原因［50-51］。一方面区域内气温与降水呈现暖湿变化趋势，另一方面退耕还林、天然林保护与多项水利补偿工程也在不断实施与干预地区自然环境，从而共同导致祁连山地区植被覆盖度整体增加［52-53］，林地、草地以及灌丛等土地覆盖类型占比提升，进一步促使目标服务供给量在2000—2020 年间整体上都处于上升趋势。西部山区由于地形对降水的影响以及人为干扰较少，供给增值更为明显。产水量、碳储存和土壤保持服务对自然环境的变化较为敏感，其单位面积供给均值与各土地覆盖单位面积供给量增加显著，但产水量在西北部疏勒河流域的山谷中却出现了减少趋势，这是由于河流与湖泊水量的增加以及人类在昌马灌区的农业灌溉日益频繁，导致该区的土壤含水量增高、实际蒸散发量增加而产水量下降［54］，所以在该区应推行滴灌、井灌井排等节水农业，提升水源涵养能力。水质净化、碳固定和生境质量服务受到祁连山地区东部的冰川积雪退缩较快的影响［55］，冰川退缩形成的裸地、草地成为了服务增长的高值区。然而由于模型基于土地覆盖的计算方式，各土地覆盖的单位面积供给能力基本无变化［56］。生境质量仅在中东部河流沿岸受人类活动影响呈减弱趋势，整体生境质量水平在植被覆盖度增加的影响下，呈现提升趋势。总体来看，祁连山地区西部在暖湿气候条件影响下非热点区面积显著减少，Ⅲ类和Ⅳ类热点区在2000年前后提出的退耕还草、围栏放牧等生态措施影响下，其面积不断增加，但在东部地区，主要分布着林地的V类热点区受到人类活动影响而呈减少趋势,在生态管理中要加强该区域林地资源的保护，协调人与自然的发展。

祁连山地区六项生态系统供给量在空间分布格局上均呈现东部高于西部的特征。这是由于东部地区受到东亚季风带来丰沛水汽的影响，水热条件优于主要受西风带影响的西部地区［57］，从而导致植被覆盖度东部大于西部。同时东部地区是林地、灌丛等生态系统供给能力较强的土地覆盖的主要分布区，影响了祁连山地区生态系统服务在空间上的整体分布。其中，产水量的空间分布与降水量有着较好的耦合关系，冰川积雪区由于其实际蒸散发量小，单位面积产水供给能力较强；土壤保持年均值的高值区则多分布在研究区东部降水较多且坡度平缓的河流谷地，该区域适宜林草生长，利于当地的土壤保持［23］；碳储存、碳固定和生境质量的林地、草地单位面积供给量较大；而水质净化由于耕地与建设用地的氮、磷污染物来源量大，导致其单位面积供给能力最强［58］。在海拔梯度上，产水量与土壤保持受山地降水的影响，随着海拔的升高而增加，而碳储存、碳固定与生境质量则呈单峰分布，在海拔3 650～3 950 m 处供给量最大，这一海拔区间也与祁连山地区中西部的最大降水带高度相当［59］；随着海拔梯度级数增加，水热条件变差，植被盖度减少，服务的供给能力降低，这与Ma 等［26］的研究结论一致。人类活动集中在较低的海拔、坡度和地面起伏度处，水质净化服务在这些区域形成高值。因此在生态管理中，对于研究区海拔低于3 350 m 以下的区域应以协调人地关系为重点，注重污染物的处理；对高海拔地区划定保护红线，防止潜在土壤侵蚀风险进一步提高；而在坡度大于15°的地区，应坚决落实退耕还林还草政策，坡度大于25°的地区只能用作林业用地，以防范土壤侵蚀，提升土壤保持能力［47］。山区相对于其周边区域，由于其林地、草地的分布热点区较强，在研究区西部分布的Ⅰ类、Ⅱ类热点区在Ⅲ类热点区的增加趋势下，呈减少变化，未来需要预防该区域的生态环境退化，促使其转为多项生态系统服务协同变化。

本文也存在以下不足：一是本文所引用的土地覆盖数据是生态系统服务计算的基础，而其分类精度在部分年份仅为80%左右，可能会带来一些系统误差。随着高分辨率土地覆盖数据集的推出，这类误差在今后的研究中有望得到改善。二是由于研究区实测的水文、气象数据有限，InVEST 模型校准中，部分数据结果缺乏校准，影响到模型计算结果的精度，未来应通过实地调查、采样等准确校准模型参数［2，6］。另外，需要延长研究年份，定量评估生态系统服务的连续变化及其与土地覆盖之间的联系，更为准确地反映出生态工程建设所带来的影响。

5 结论

本文基于InVEST 模型、CASA 模型对2000—2020 年祁连山地区产水量、碳储存、土壤保持、水质净化、碳固定和生境质量六项主要生态系统服务进行定量评估，并结合自然和人类活动影响分析其时空演变特征，提出了生态管理的建议，对研究区生态恢复和保护有一定的借鉴作用。研究结论如下：

在时间上，祁连山地区土地覆盖变化面积较小，但变化类型复杂，六项生态系统服务整体空间格局较稳定，各项服务供给量总体呈上升趋势，各土地覆盖供给能力也有所提升，山区是服务供给量增加的高值区，但在局部地区也存在受人类活动影响而出现生态环境退化现象，西部的整体生态环境趋向好转，Ⅲ类和Ⅳ类热点区面积呈增加趋势，其余热点区则相反。在进行生态管理时须因地制宜，制定针对性政策。

在空间上，祁连山地区六项生态系统服务均呈“东高西低”的分布格局，林地与草地单位面积供给能力较高，耕地与建设用地的单位面积供给能力对水质净化影响最大，3、4 级梯度处的碳储存、碳固定和生境质量均值较高，产水量和土壤保持高值区5级梯度处，而水质净化则在1级梯度处均值最大，Ⅴ类热点区多分布在东部山区林地与林草过渡带，要注意协调该区域的人地关系，注重林地资源保护。

致谢：本研究工作得到国家重点研发计划项目（2019YFC0507402）资助和兰州大学超算中心支持。