白驹，乔国华，高宇，陈广雨

1. 国家基础地理信息中心，北京 100830；

2. 泰安市不动产登记和交易中心，泰安 271001；

3. 自然资源部信息中心，北京 100830；

4. 黑龙江铁汇测绘工程有限公司，哈尔滨 150000

1 引 言

生态系统服务及驱动机制研究一直是地理学、生态学和可持续科学领域研究前沿（傅伯杰等，2020）。运用多源自然地理时空信息开展长序列生态系统服务时空动态监测及驱动机制研究，可为确定优化多目标生物多样性保护、服务组合交付和关键利益相关者偏好的管理实践提供支撑，也是时空信息赋能我国当前高质量发展的主要形式之一（单卫东等，2023；陈军等，2023）。由于生态系统内部组成、结构迥异，在复杂气候变化、人类活动和其他因素影响下，各种生态系统服务呈现时空动态的异质性，且存在驱动因素的不确定性，尤其是在生态脆弱的高寒地带（Ma 等，2021）。

受高寒、干旱的极端自然环境影响，高寒地带草地生态系统内部结构、组成对气候变化和人类活动的响应尤为敏感，并显著反馈于生态系统服务功能（姜翠红等，2016）。气候变化通过温度、降水等气候因子直接影响生态系统服务，通过影响土地利用/覆盖的强度间接改变生态系统服务，并且气候变化、植被动态等引发的多类型生态系统服务供给量的变化，都可能导致生态系统服务驱动机制的复杂变化（张明珠等，2021）。部分研究选取不同地域开展了相关分析。潘梅等（2020）对2000~2015年京津冀地区生态系统服务研究发现，受气候暖湿和生态修复工程双重因素影响，该区域水源涵养量、防风固沙量均呈不显著上升趋势。云南省生态系统产水量在1992~2019 年呈先增大后减小的波动变化趋势，但是各年度的产水量空间分布格局趋于一致，且证实降水量和实际蒸散量是导致产水服务空间分异的主要驱动因子（黄欣等，2022）。此外，有研究还选取特定的生态系统展开了生态系统服务时空动态监测及驱动机制分析。景海超等（2022）对青藏高原地区那曲市生态系统服务的时空动态及自然驱动因素研究发现，产水量与土壤保持服务呈现下降趋势，生境质量趋于稳定，而固碳服务增长5.6%。降水量、归一化植被指数（normalized difference vegetation index，NDVI）是决定产水量和固碳空间异质性的关键因素，土壤保持主要受坡度因子影响。相比土地利用和人类活动直接干扰强烈地区，高寒地带草地生态系统受到自然因素影响更为显著（姜翠红等，2016）。但生态过程控制下不同生态系统服务对自然因子的响应机制尚不明晰，极有必要进一步揭示高寒地区长期生态系统服务时空动态及自然因子驱动机制。

黄河流域青海段是黄河的发源地，也是重要的水源涵养区与补给区，具有很强的水源涵养、水土保持及碳固存的功能，是维系黄河中下游地区生态安全的天然屏障。近年来，虽然国家在该区域实施了区划轮牧、草原修复等生态保护措施，草原退化和河水流量剧减的趋势得到有效遏制。但是，由于气候变化和过度人类活动的制约，该区域生态环境脆弱的局面尚未从根本上扭转，生态系统服务时空动态仍不明确，自然因素对生态系统服务的驱动机制仍不明晰。因此，本文以黄河流域青海段为研究区域，选取高寒草地生态系统产水量、土壤保持、防风固沙和固碳四项典型生态系统服务功能作为研究对象，分析长时间序列下高寒生态系统产水量、土壤保持、防风固沙和固碳四项典型生态系统服务的时空动态，揭示生态系统服务时空动态的自然驱动机制，以为制定黄河流域上游高寒地区生态系统服务管理对策和可持续发展提供科学支撑。

2 研究数据与方法

2.1 研究区概况

黄河流域青海段位于青藏高原东南部，包括黄河河源区和上游区域，约占黄河流域总面积的23.8%，有重要水源涵养、土壤保持、防风固沙、固碳和生物多样性保护功能，是中国乃至亚洲重要生态安全屏障区。研究区地势西高东低，平均海拔3854 m，属高原大陆性气候，常年平均气温为2.8 ℃，多年平均降水量441.0 mm，呈东南向西北递减趋势。区域太阳辐射较强，年辐射量140~177 kJ/cm2，日照时长2400~2800 h。根据全国水资源及黄河流域水资源分区可知，该区域划分为河源区、玛曲段、龙羊峡段、兰州段和湟水流域（图1）。河源区为黄河发源地区，海拔4200~4800 m，属高原湖泊沼泽地貌，地势较平缓；玛曲段为重要水源涵养区，海拔3500~4000 m，地形起伏大，多高山深谷；龙羊峡段多高山峡谷地貌，间有谷地及高山草地。兰州段峡谷盆地相间，河道高程降至2670 m 以下，分布耕地、林地。东北部为支流湟水流域；子流域北部多高山河沼泽，以畜牧业为主；南部谷地、黄土丘陵地貌，以种植农业为主。

图1 研究区地表覆盖类型分布Fig.1 Surface cover type distribution in the study area

2.2 数据源

本研究中生态系统服务核算和自然驱动因子以气象、地形、土壤等数据集为基础。

（1）气象数据。地表温度、参考蒸散发数据分别采用美国国家航空航天局（National Aeronautics and Space Administration，NASA）发布的MOD11A2、MOD16A2 8 d 合成产品，分辨率1 km；降水数据为NASA 发布雷达与辐射计组合的网格降水量估计数据，包括逐月的热带降水测量计划（tropical rainfall measurement mission，TRMM）和全球降水观测（global precipitation measurement，GPM）两代全球卫星降水产品，分辨率为0.25º×0.25º；风速数据采用欧洲中期天气预报中心（European Centre for Medium-range Weather Forecasts，ECWMF）提供的哥白尼气候变化服务局（Copernicus Climate Change Service，C3S）的能源数据集中10 m 高度处二维水平空气速度，分辨率为0.25°×0.25°，单位为m/s。

（2）地形数据包括海拔、坡度。海拔采用先进星载热发射和反射辐射仪全球数字高程模型（advanced spaceborne thermal emission and reflection radiometer global digital elevation model，ASTER GDEM）数据，分辨率为30 m，数据来源于地理空间数据云（https://www.gscloud.cn/）；坡度数据基于海拔数据，通过ArcGIS 的Surface 模块计算获取。

（3）土壤数据来源于国家冰川冻土沙漠科学数据中心（National Cryosphere Desert Data Center，NCDC）1∶100 万土壤数据库，含土壤机械组成、有机碳含量、碳酸钙含量、土壤最大根系深度、植被可利用水分信息等。土壤湿度数据以增强植被指数（enhanced vegetation index，EVI）和地表温度数据为基础，选用温度植被干旱指数（temperature vegetation dryness index，TVDI）数据模型生产。土地利用（覆被）采用2000 年、2010 年、2020 年GlobeLand30 数据，源于全国地理信息资源目录服务系统（https://www.webmap.cn/），分辨率为30 m×30 m。EVI 数据采用NASA 发布的MOD13A2 16 d合成产品，分辨率为1000 m。植被覆盖度可采用EVI，通过像元二分法计算获得（Jiapaer 等，2011）。净初级生产力（net primary productivity，NPP）数据采用NASA 提供的MOD17A3HGF 年数据，分辨率为500 m。

2.3 研究方法

2.3.1 生态系统服务核算模型

本研究针对黄河流域上游高寒草地生态系统的典型生态系统服务功能，分别选取生态系统服务核算模型，测算产水量、土壤保持、防风固沙和固碳四项生态系统服务物理量。

产水量主要基于水量平衡原理，结合区域气候因素、地形及土壤特征等因素，以降水作为系统的输入，蒸散量作为系统的输出，两者差值即为区域产水量（Sharp 等，2016）：

式中，Y(x)为生态系统类型栅格x的年产水量，mm；AET(x)为栅格x的年实际蒸散量，mm；P(x)为栅格x上的年降水量，mm；PET(x)为年内潜在蒸散量，mm；ω为表征自然气候–土壤性质的非物理经验参数值。

土壤保持是生态系统通过其结构与过程减少由于降水所导致的土壤侵蚀的作用，是生态系统重要调节服务之一。生态系统土壤保持服务通过土壤保持量来衡量，为无植被覆盖情况下的潜在土壤侵蚀模数与现实土壤侵蚀模数的差值，土壤侵蚀模数采用修正通用水土流失方程（revised universal soil loss equation，RUSLE）计算得到（Renard 等，1997）。

防风固沙主要针对生态系统防治风蚀的服务功能进行测算。选用生态系统中气候因素、植被覆盖状况、土壤质地等因素综合评价生态系统减少风蚀量、防风固沙的能力，为裸土条件下土壤风蚀和地表覆盖植被条件下的土壤风蚀量的差值，通过修正风蚀模型（revised wind erosion equation，RWEQ）计算风蚀区土壤风蚀模数（Gregory 等，2004）：

式中，G为防风固沙量，t/ha；SLQ为潜在风蚀量，t/ha；SL为实际风蚀量，t/ha；Qmax_Q为潜在转运量，t/ha；Qmax为实际转运量，t/ha；SQ为潜在关键地块长度，m；z为下风向最大风蚀出现距离，m，本次计算选取50 m；WF 为气候因子，t/km，表征了在考虑降水、温度、日照及雪盖等因素的条件下风力对土壤颗粒的搬运能力；EF 为土壤可蚀性因子，为一定土壤理化条件下土壤受风蚀影响大小；SCF 为土壤结皮因子，为一定土壤理化条件下土壤结皮抵抗风蚀能力的大小；K′为地表粗糙度因子，由地形所引起的地表粗糙程度对风蚀影响的反映；C为植被覆盖因子，植被条件对防止风蚀发生程度的贡献。

植被通过光合作用吸收大气中的二氧化碳是陆地生态系统固碳的主要方式之一。NPP 作为地表碳循环的重要组成部分，不仅直接反映了植被群落在自然环境条件下的生产能力，而且是判定生态系统碳源汇和调节生态过程的主要因子，在全球变化及碳平衡中扮演着重要的作用。CASA （Carnegie-Ames-Stanford approach）模型可以通过植被吸收的光合有效辐射（absorbed photosynthetic active radiation，APAR）和光能利用率(ε)两个变量估算植被NPP。由CASA 模型可知，植被每生产1 kg 有机物，能固定1.62 kg CO2。因此，植被固碳量：

式中，GCS为植被年固碳量，t/ha；RC为二氧化碳中碳的含量，为27.27%。

2.3.2 冗余分析法

冗余分析（redundancy analysis，RDA）是响应变量矩阵与解释变量矩阵之间多元多重线性回归的拟合值矩阵的主成分分析（principal component analysis，PCA）方法，也是多响应变量回归分析的拓展，广泛应用于多变量响应数据建模（董旭辉等，2007）。Canoco5.0 中RDA 分析可以显示每个环境因子解释率和显著性检验的结果（图2）。响应变量与解释变量箭头之间夹角的余弦值可以判断它们之间的相关性。Rv1 和Rv2 是响应变量，Ev1 、Ev2 和Ev3 为解释变量。∠a接近90°，即接近正交，表明Rv1 和Ev1 之间的相关性很小。∠b小于90°，表明Rv1 和Ev2 之间存在正相关；角度越小，则正相关性越大。∠c大于90°，表明Rv1 和Ev3之间存在负相关；角度越大，则负相关性越大。本文中采用随机样点1000 个作为解释变量。

图2 冗余分析二维排序Fig. 2 Double sequence diagram of RDA

3 结果与分析

3.1 生态系统服务空间格局特征

研究期内黄河流域青海段产水量、土壤保持、防风固沙和固碳四种生态系统服务物理量均值分别为340.3 mm、1328.2 t/ha、23.73 t/ha、1025.0 kg/ha（表1）。由图3 看出，产水量整体呈东南向西北递减趋势，高值区集中分布在玛曲段、龙羊峡段南部和湟水流域北部，年产水量深度超过500 mm；低值区则分布在河源区、兰州段和龙羊峡段北部，大部分地区年产水量深度低于200 mm。研究区内76.5%地表土壤保持值低于平均值，高值区集中分布在龙羊峡段、兰州段东北部和湟水流域下游地区，其他水资源分区有零星分布。防风固沙空间分布整体较为均匀，88.01%区域低于均值。高值核心区位于研究区东北部边界区域，并呈放射状向外衰减。固碳格局整体呈东部向西部递减趋势，局部区域呈现空间异质性。

表1 2000 年、2010 年、2020 年黄河流域青海段典型生态系统服务物理量均值统计Tab.1 Average value of typical ecosystem services in Qinghai section of Yellow River basin in 2000, 2010, and 2020

图3 2000~2020 年黄河流域青海段生态系统服务物理量均值空间格局Fig.3 Spatial pattern of typical ecosystem services in Qinghai section of Yellow River basin from 2000 to 2020

3.2 生态系统服务时空动态特征

从时空动态来看，流域尺度下产水量、土壤保持、防风固沙的年物理量均值都呈现增加趋势（表1）。流域内2010 年产水量比2000 年整体提升 8.9%，2020 年产水量比 2010 年整体提升44.6%，提升明显区域集中在南部、东部区域（图4）。2010 年土壤保持物理量均值比2000 年增加12.8%，占总面积的72.7%，增加明显区域集中分布在龙羊峡段东北部和河源区西部，而湟水流域下游地区土壤保持略有降低。2020 年土壤保持物理量均值比2010 年增加8.1%，增加面积占流域总面积34.8%，主要分布在龙羊峡段东北部。2010年防风固沙物理量均值比2000 年增加10.7%，增加区域占总面积54.5%；2020 年防风固沙物理量比2010 年增加10.1%，增加区域占总面积34.3%。两期防风固沙物理量增加高值区都集中在湟水流域东部边界区域。流域内2010 年固碳物理量均值呈现比2000 年升高29.2%，增加明显区域集中在龙羊峡段、兰州段和湟水流域下游区域；2020 年固碳物理量均值比2010 年呈现略微降低，比例为1.7%，降低区广泛分布在玛曲段。

图4 黄河流域青海段生态系统服务时空动态格局Fig.4 Spatial-temporal dynamic of typical ecosystem services in Qinghai section of Yellow River basin

3.3 生态系统服务与自然因子冗余分析

根据RDA 分析结果，2000 年、2010 年、2020年地表温度、风速、海拔、植被覆盖状况、土壤湿度、坡度和降水量七个自然因子解释变量对四项典型生态系统服务物理量时空动态贡献率分别为47.0%、38.1%、41.2%，存在年际波动。2000年、2010 年植被覆盖状况、海拔是两个贡献率最大的解释变量，总贡献率分别为29.3%、26.5%。土壤湿度、降水的解释变量仅为15.4%、10.3%，相对植被状况和海拔较低。2020 年海拔和地表温度是影响生态系统服务最关键两个自然因子，解释变量的贡献率合计为29.9%，植被状况和降水的贡献率仅为10.6%。风速、坡度对生态系统服务的贡献率最小。

各自然因子对不同类型生态系统服务物理量时空动态贡献率不同，且对单一生态系统服务物理量的贡献率也存在年际差异（图5）。降水是研究期内决定产水量的首要正向因子，即降水量的多少直接决定该区域的产水量。土壤湿度与产水量表现出正相关性，地表温度与产水量负相关，即流域内地表温度越高的地区产水量越低。坡度对产水量的影响呈现降低趋势。2010 年、2020 年，海拔由负向指标转变成为除降水外控制产水量的关键正向指标。在对土壤保持时空动态解释中，地表温度是唯一贯穿整个研究期的正向驱动因子。海拔是决定土壤保持物理量的主要逆向驱动因子，其次为植被状况、土壤湿度。风速、降水、坡度对土壤保持服务的贡献具有不确定性。土壤湿度、植被状况和坡度是控制防风固沙物理量的正向驱动因子，温度是关键的逆向驱动因子。2000 年风速对防风固沙的影响并不明显，而2010 年、2020 年风速与防风固沙服务呈现明显的负相关性。2000 年、2010 年，植被状况、海拔、降水、坡度、土壤湿度是控制生态系统固碳服务的正向驱动因子，温度与固碳服务物理量呈现负相关性；而2020 年温度与固碳服务物理量呈现正相关，降水、海拔和土壤湿度成为逆向驱动因子。

图5 生态系统服务与自然因子冗余分析二维排序Fig.5 Double sequence diagram of RDA between ecosystem services and natural factors

4 讨 论

4.1 自然因子对生态系统服务时空动态的影响

研究期内该区域产水量年均值为278.6~438.9 mm，整体呈现东南向西北递减趋势，年际变化呈增强态势。根据水量平衡原理，降水量、实际蒸散量以及二者之间的平衡是直接控制区域产水量能力、格局的关键因素。气候变化趋势下，近20 年以来黄河流域青海段降水量整体呈现增加趋势（姬兴杰等，2022）。实际蒸散量表现出显著降低是黄河流域青海段产水量增加的原因。

研究区土壤保持物理量呈现明显年际增加趋势，土壤保持服务强度高值区在中北部和东部。根据修正通用水土流失方程，除受降水等气候因子影响外，地形地貌、植被覆盖情况都是影响土壤保持物理量的自然因素。龙羊峡段东北部地区土壤质地以土质疏松沙地为主，增加了潜在土壤侵蚀量，利于植被发挥更明显的土壤保持服务功能。兰州段东北部和湟水下游多峡谷、盆地交错分布，河道高程落差较大，复杂地貌条件是该区域潜在土壤侵蚀量高于其他区域主要诱因，但区域多耕地和林地，有较高植被覆盖度和较好土壤侵蚀防治措施，可以减少地表径流，降低实际土壤侵蚀量（陈朝良等，2021）。在土壤侵蚀防治措施增强的情况下，降水对土壤保持物理量年际变化影响显著。

研究区防风固沙服务物理量呈现明显年际增加趋势，高值区位于东北部边界区域。气候因素方面风速对固沙功能有根本影响，决定潜在风蚀量，而不同土地覆盖类型通过防风功能发挥固定土壤作用，决定实际风蚀量大小，防风功能上呈现森林>农田>草地>未利用地（Li 和Xu，2019）。该区域主要草地覆盖，少数耕地、森林，土地覆盖类型的差异导致防风固沙服务能力差异。东北部边界区域为祁连山南麓，森林和山体分布发挥较强的防风固沙功能。

研究区固碳能力呈现年际变动趋势，龙羊峡段、兰州段和湟水流域下游区域是固碳高值区。近20 年来，青藏高原地区暖湿气候变化趋势利于植被生长茂盛，NPP 呈现增长趋势，直接反映固碳的年际动态（夏冰等，2023）。不同植被类型之间的NPP存在不同程度的差异，决定了该区域森林>草地>耕地>未利用地的固碳能力强弱，因此，植被类型也是决定固碳物理量空间格局的重要因素（Wei 等，2022）。流域内东部多森林、西北多高寒草原分布决定了固碳物理量呈现整体呈现东部向西北递减趋势。

4.2 自然因素对生态系统服务时空动态驱动的复杂性

生态系统服务时空动态是一个由自然条件、人类活动等多种因素直接或间接共同驱动的复杂过程。本研究选取的七个自然因子对2000 年、2010年、2020 年生态系统服务物理量时空动态解释变量仅为38.1%~47.0%，说明可能存在其他驱动因素控制着黄河流域青海段的生态系统服务物理量时空动态。

本研究发现植被覆盖状况在对四项生态系统服务时空动态解释中较其他自然因子具有较高贡献率。植被通过蒸腾作用、调节小气候等改变生态系统产水量，通过对降水的截留和再分配作用改良土壤结构的方式调节土壤保持服务（赵苗苗等，2017；李鑫等，2018）。植物可以利用根茎叶削弱地表风力、固定土壤减缓风蚀作用，通过光合作用固定碳元素。水热条件、土壤特性是控制陆地生态系统植被类型和生长状况的核心因素。海拔是影响水热条件的主要因素，通过对水热条件的再分配影响土壤特性。同样地，坡度成为各类生态系统服务物理量时空格局的重要解释变量。

各自然因子对不同类型生态系统服务物理量时空动态贡献率不同，这基本符合自然因子对生态系统服务物理量时空动态的贡献与各类生态系统服务物理量核算中涉及的生态过程具有一致性的规律（苏常红和王亚璐，2018）。产水量是区域水循环过程中综合考虑降水量与实际蒸散量收入支出平衡的产物。降水是自然界中土壤水分最直接、主要的补给方式，其必然也与产水量呈现正相关性。在排除水分可获得性限制前提下，温度升高导致实际蒸散量增加，为间接控制产水量的负向因素。土壤保持物理量与降水之间的关系并不是简单的线性关系，而与降水量是否处于植被生长的水分胁迫阶段密切相关（Zhao 等，2022）。在水分充足的前提下，温度利于植被生长、恢复，表现出对土壤保持时空动态具有显著积极效应。地形通过改变气候、植被因子引起土壤保持物理量的年际变化。

在黄河流域青海段，除温度外，风速和降水未对防风固沙表现出稳定的贡献。潜在风蚀量决定防风固沙服务物理量的潜力，其产生前提是风速大于临界起沙风速。黄河流域青海段不属于大风天气核心区，模型采用的平均风速可能将大于临界的瞬间起沙风速（5 m/s）过滤掉，结果将弱化风速对黄河流域青海段防风固沙服务的贡献（王式功等，2003）。黄河流域青海段降水与大风天气的时间不重合，削弱了降水在生态系统防风固沙功能中的效用。土壤湿度、海拔、坡度和植被覆盖情况对防风固沙时空动态起着重要积极作用。首先，土壤湿度增加土壤颗粒之间的黏聚力以抵抗风蚀，间接改善地表植被覆盖情况以减小实际风蚀量；其次，海拔和坡度差异引起的多风大风天气和土壤颗粒的迁移利于潜在风蚀，且多样地貌类型提供的土壤和气候条件影响植被的生长、分布和形态，影响防风固沙能力；最后，植被通过直接覆盖地表使其免受吹蚀、分散风动能减少气流与土壤颗粒之间动能传递、直接阻挡土壤颗粒运动的方式发挥防风固沙功能（牛丽楠等，2022）。这与本研究EVI 与防风固沙呈现显著正相关结果一致。

本研究的固碳核算方法仅考虑了植被通过光合作用自然的碳封存过程，因此植被特征变化是决定固碳时空动态的直接因素，EVI 与固碳时空动态呈现正相关性。此外，气候、地形和土壤都会通过改变植被特征或光合作用直接或间接影响固碳时空动态（林枫等，2021）。在近20 年来青藏高原暖湿气候变化趋势下，降水、海拔、土壤湿度、温度对固碳作用亦存在变动性。根据2000 年、2010 年和2020 年三期研究结果，大部分自然因子对生态系统服务时空动态的驱动机制比较明晰，但产水量和坡度、土壤保持和风、固碳和风等并未显示出显著的相关性或它们的相关性随着时间的推移并不稳定，这表明气候因素对生态系统服务的变化是波动的、复杂的，需要更长时间序列研究揭示真实趋势。

5 结 论

黄河流域上游地区高寒草地生态系统是我国重要的生态屏障地区。本研究通过对该区域水源涵养、土壤保持、防风固沙和固碳生态系统服务功能的时空动态监测，探明了四项典型生态系统服务功能的时空动态特征，并分析了其自然驱动机制。

黄河流域青海段产水量、土壤保持、防风固沙、固碳四项种生态系统服务物理量均值分别为340.3 mm、1328.2 t/ha、23.73 t/ha、1025.0 kg/ha，并呈现四种不同的地理空间格局。从时空动态来看，2000 年、2010 年、2020 年流域的产水量、土壤保持、防风固沙物理量呈现增加趋势。虽然2020年固碳物理量比2010 年略微降低，但相对2000 年整体仍然呈现增加态势。2000 年、2010 年、2020年，降水、地表温度、风速、土壤湿度、植被覆盖状况、海拔、坡度七个自然因子解释变量对四项典型生态系统服务物理量时空动态贡献率分别为47.0%、38.1%、41.2%，存在年际波动。植被覆盖状况、海拔是2000 年、2010 年两个贡献率最大的解释变量，而海拔、地表温度是2020 年最关键两个自然因子，风速和坡度对四项生态系统服务时空动态影响较弱。降水是决定产水量的主要正向因子，地表温度为负向因子。在对土壤保持时空动态解释中，地表温度是唯一贯穿整个研究期的正向驱动因子。海拔是决定土壤保持物理量的主要逆向驱动因子，风速、降水、坡度对土壤保持服务的贡献具有不确定性。土壤湿度、植被覆盖状况和坡度是控制防风固沙物理量的正向驱动因子，温度是关键的逆向驱动因子。温度是控制固碳的重要因子，但正反向贡献的变化表明温度与固碳服务之间存在非线性关系。