基于植被生产力的黄土高原地区生态脆弱性及其控制因子分析

2022-12-09杨艳周德成宫兆宁刘子源张良侠

生态环境学报 2022年10期

杨艳，周德成，宫兆宁，刘子源，张良侠*

1. 南京信息工程大学应用气象学院，江苏南京 210044；2. 首都师范大学资源环境与旅游学院，北京 100048；3. 资源环境与地理信息系统北京市重点实验室，北京 100048

生态环境是人类赖以生存的基本条件。近年来，随着全球气候变化和人类活动加剧，生态系统的自我调节能力快速下降，呈现出脆弱化的趋势，严重威胁着区域生态安全和经济社会的可持续发展（Hu et al.，2021；徐广才等，2009；徐兴良等，2022）。依据IPCC第五次评估报告，生态脆弱性被定义为生态系统容易受到不利影响的倾向或者习性，它包括生态系统对不利影响的敏感性以及对新的外界条件的适应性两个方面（IPCC，2014；苏胜亮等，2022）。科学合理的对区域生态脆弱性进行量化评估，对理解生态系统的功能和过程以及因地制宜的进行生态系统修复和保护具有重要意义。

生态脆弱性评价的关键在于构建科学的评价指标体系。当前相关研究采用的评价指标体系大致可以分为综合指标体系和单一指标体系两种（杨飞等，2019）。其中，综合指标体系大多基于不同的评估模型框架，从自然风险、人为干扰和社会经济等方面选取多个指标构建而成（张学玲等，2018）。同时，在进行脆弱性评估过程中，往往需要依据不同的方法（如主成分分析、层次分析和专家打分法）对选取的多个指标进行权重分配（鄢继尧等，2020）。例如，张学渊等（2021）基于“压力-状态-响应”评估模型框架，从气候、土壤、植被、地形地貌和社会经济等方面选取了 13个指标，利用空间主成分分析法对 2000—2018年西北干旱区生态脆弱性进行了评估；朱琪等（2021）基于“成因及结果表现”评估模型框架，从地形、地表、气象和社会经济等方面选取 13个指标，结合空间主成分分析法评估了 2005—2015年东北森林带生态脆弱性状况；金丽娟等（2022）基于“敏感-恢复-压力”评估模型框架，选取了 17个涵盖地形、地表、气象、结构、功能和压力等多方面的评价因子，并基于层次分析法对 2005—2018年四川省进行生态脆弱性评价。然而，由于考虑范围广和指标数量多，综合指标体系存在指标因子权重分配具有主观性以及因子之间相关性高等问题。相较而言，采用能够综合反映区域生态脆弱性的主要因子建立单一指标评价体系具有较强的可靠性和可操作性（Zhang et al.，2017；肖桐等，2010）。

生态系统生产力是生态系统功能的重要指标和生态系统承载力的物质基础，能够有效反映出生态系统的健康状态，对外界环境变化具有高度敏感性，是衡量生态脆弱性的重要依据（于贵瑞等，2020）。国内外学者基于生态系统净初级生产力（net primary productivity，NPP）和总初级生产力（gross primary productivity，GPP）已成功地对不同地区的生态脆弱性进行了评估。如以NPP作为评价指标，分别对1989—2007年的三江源地区、1961—2000年长江中下游地区、2000—2010年岷江上游地区进行了生态脆弱性评估（Zhang et al.，2017；肖桐等，2010；於琍等，2012）；以 GPP为评价指标，分别对2000—2015年西南地区以及2000—2014年天山—塔里木河绿洲地区生态脆弱性进行了评估（何敏等，2019；王鹤松等，2021）。这些研究工作有助于提高对生态脆弱性问题的理解，但在探究生态脆弱性控制因子时仍存在一定局限性。如通过对各控制因子（如温度、降水、高程、坡度）作分区统计，来分析其对生态脆弱性的影响（肖桐等，2010；王鹤松等，2021）；分别基于沿各控制因子的变异梯度随机选取的40个采样点数据和基于整个研究区的所有数据利用线性回归分析方法，分析每个控制因子与生态脆弱性的相关关系（Zhang et al.，2017；何敏等，2019）。受研究方法的限制，这些研究仅初步分析了区域尺度上各控制因子对生态脆弱性的单一影响，难以剔除各控制因子间的相互作用，更不能定量表达各控制因子对生态脆弱性的贡献率。地理探测器是探测空间分异性并揭示其背后驱动力的一组统计学方法，该方法不需要对影响因子预设限定，且无线性假设，计算得出的q值具有明确的物理含义，能够客观反映出自变量对因变量的解释力，因此能够较全面的探究生态脆弱性时空分布格局的控制因子（王劲峰等，2017）。

黄土高原处于干旱、半干旱和半湿润区的过渡带，植被退化和水土流失问题严重，相对落后的经济、不合理的土地利用以及煤炭资源的开发更加剧了该地区生态环境的退化，属于中国典型的生态脆弱区（贺鹏等，2022；张家政等，2022）。本研究以黄土高原为研究区，选择NPP作为评价指标，依据IPCC中有关生态脆弱性的定义首先计算敏感性和适应性，进而对 2001—2020年黄土高原生态脆弱性的空间分布格局进行了评估，并结合地理探测器定量分析了生态脆弱性的控制因子及其交互作用。研究结果有助于深化对该地区生态系统对气候变化和人类活动响应的认识，并为当地采取针对性措施以减少外部胁迫对生态系统产生的不良影响、促进可持续发展和规划提供科学参考和理论支持。

1 材料和方法

1.1 研究区概况

黄土高原（ 33°43′— 41°16′N ， 100°54′—114°33′E）地处中国中部偏北，位于黄河中上游和海河上游地区，总面积约6.40×105km2（刘静等，2020），覆盖了青海、甘肃、宁夏、内蒙古、陕西、山西和河南等7个省区（图1）。该区具有冬季寒冷干燥、夏季湿润温暖的暖温带大陆性季风气候特征，年均温度4—12 ℃；降水分布不均匀，自西北向东南逐渐增多，年均降水量400—600 mm（陈剑南等，2022）。区域内以黄土覆盖为主，土体疏松，抗侵蚀能力较弱，从而形成了沟壑纵横支离破碎的典型黄土地貌。土地利用类型以草地、耕地和林地为主，面积占比分别为41%、32%和15%。植被由东南向西北可划分为森林带、森林草原带、典型草原带、荒漠草原和草原化荒漠带（修丽娜等，2019）。能矿资源丰富，是全国能源安全保障基地和生态安全屏障的重要组成部分（王佳宁等，2020）。

图1 2010年黄土高原地区土地利用情况Figure 1 Land use in the Loess Plateau in 2010

1.2 数据来源与预处理

本文使用的数据主要包括遥感、气象、高程、土地利用/覆被、人均GDP、人口密度和坡度数据。遥感数据包括NPP和归一化植被指数（normalized difference vegetation index，NDVI）两种，时间跨度为 2001—2020年，均从美国地质调查局网站（https://lpdaac.usgs.gov）下载获取。其中，NPP来自MODIS数据集中的MOD17A3H产品，时间分辨率为年，空间分辨率为500 m；NDVI来自MODIS数据集中的MOD13A2产品，时间分辨率为16 d，空间分辨率为 1 km。本研究利用最大值合成法将NDVI数据由16 d尺度转换为年尺度。气象数据为覆盖研究区及其周围148个气象站点的2001—2020年的年均气温和年均降水量，来源于中国气象科学数据共享网（http://data.cma.cn）。本研究基于气象观测数据，利用AUSPLIN插值软件结合数字高程模型（digital elevation model，DEM）进行空间插值，生成1 km×1 km空间分辨率的气象要素栅格图像。DEM、2010年土地利用情况、人均GDP和人口密度数据均来自于中国科学院资源环境科学数据中心(http://www.resdc.cn)。其中，DEM数据空间分辨率为90 m，土地利用、人均GDP和人口密度数据空间分辨率为1 km。坡度数据由DEM计算得到。为了便于分析不同指标因子与生态脆弱性之间的关系，本研究将所有数据统一重采样为Albers投影下1 km×1 km的栅格数据。

1.3 脆弱性计算方法

选取NPP作为生态脆弱性评价指标，根据IPCC第五次评估报告中有关脆弱性的概念进行生态脆弱性的计算（IPCC，2014）。生态脆弱性被定义为生态系统容易受到外界破坏或伤害的程度，可由敏感性和适应性的差值计算得到：

式中：

V——生态脆弱性；

S——敏感性，表征了生态系统受气候变化或其他外界扰动的影响；

A——适应性，表征了生态系统受自身调节和恢复能力的制约。

S由NPP的年际波动情况来表示，计算公式如下：

式中：

Fi——该地区第i年的NPP值（n=20）；

——黄土高原NPP的平均值。

A由NPP年际变率线性拟合趋势线的斜率来表示，计算公式如下：

式中：

x——自然数 1、2、3…，对应于 2001—2020年的时间序列；

y——NPP年际变率的绝对变化量。

值得注意的是，由公式（2）和（3）计算得到的S和A数值不在同一量纲，难以直接用于数据分析和生态脆弱性计算。因此，本研究基于极差标准化方法（Hu et al.，2021）对S和A数据进行了标准化处理。

为了进一步揭示黄土高原地区生态脆弱性空间变异特征，需要对生态脆弱性指数进行分级。常用的数据分级方法包括自定义间隔、等距、标准差和自然断点等，其中自然断点分级法是一种根据数值统计分布规律分级和分类的统计方法，可使得各级内部方差之和最小且差异最大化（Liu et al.，2017）。本研究采用自然断点分级法确定生态脆弱性指数分级阈值，将其划分为微度脆弱、轻度脆弱、中度脆弱、重度脆弱和极度脆弱5级。

1.4 控制因子分析

为进一步探究黄土高原生态脆弱性的控制因子，选取年均温、年均降水量、海拔、坡度、人均GDP、人口密度和NDVI 7个指标（表1），利用地理探测器模型定量分析了探测器包括因子探测、交互作用探测、风险区探测和生态各因子对生态脆弱性空间分布的影响。利用地理探测器模型定量分析了各因子对生态脆弱性空间分布的影响。地理探测器包括因子探测、交互作用探测、风险区探测和生态探测四个部分（王劲峰等，2017）。其中因子探测利用因子解释力衡量自变量因子对因变量变化产生的贡献力大小，可以反映出黄土高原生态脆弱性的主要控制因子；交互作用探测通过比较两个指标的单因子解释力之和及其交互作用的解释力，可用于探测自变量因子两两之间的交互后对生态脆弱性的影响程度。地理探测器模型主要用q值度量解释力结果，其值越大说明自变量对于因变量空间分异的解释力越强，反之则越弱。首先，将生态脆弱性指数V作为因变量Y，将选取的7个控制因子指标作为自变量X，借助ArcGIS 10.7软件在15 km×15 km格网尺度统计脆弱性和各自变量因子平均值；其次，利用自然断点法将各自变量因子重分类为5级，由数值量转变为类型量，并获取脆弱性数据与因子分类数据的空间关联属性表；最后将属性表导入地理探测器运行得到q值，对影响黄土高原生态脆弱性空间分布的控制因子进行探测分析。

表1 影响因子Table 1 Impact Factor

2 结果

2.1 黄土高原地区生态脆弱性空间分布格局

2001—2020年黄土高原生态脆弱性整体较高，中度及以上脆弱等级地区面积占比为61%（表2），且呈现出西北高而东南低的空间分布格局（图2c）。对比而言，适应性空间分布与脆弱性刚好相反，呈现西北低而东南高的格局，而敏感性呈现中部高而东部及南部低的格局（图2a、b）。具体而言，微度和轻度脆弱区所占面积分别为 1.05×105km2和1.46×105km2，主要分布在黄土高原东南部，如南部的六盘山和子午岭地区，以及东部的吕梁太行山区，具有敏感性低且适应性极高的特点；重度和极度脆弱区所占面积则分别为 1.49×105km2和 6.80×104km2，主要分布在黄土高原中部的丘陵沟壑区及西南部的高塬沟壑区，敏感程度高且适应性较低。此外，不同土地利用类型间的脆弱性差异也较大（图3）。林地整体脆弱性较低，以微度、轻度和中度脆弱为主，面积占比超过 85%；其次为建设用地、耕地和草地，中度及以下脆弱区面积占比分别为76%、66%和 59%；未利用地的脆弱性程度最高，以中度、重度和极度脆弱为主，面积占比超过86%。

图2 2001—2020年黄土高原地区敏感性（a）、适应性（b）和生态脆弱性（c）空间分布特征Figure 2 Spatial distributions of sensitivity, adaptability and ecological vulnerability in the Loess Plateau in 2001-2020

图3 2001—2020年黄土高原不同土地利用类型下各生态脆弱性等级面积占比Figure 3 Area proportion of each ecological vulnerability under different land use types in the Loess Plateau in 2001-2020

表2 黄土高原不同脆弱性等级所占面积及占比Table 2 Area and proportion of different vulnerability levels in the Loess Plateau

2.2 黄土高原地区生态脆弱性空间格局的控制因子分析

7个影响因子的q值大小排序为NDVI＞降水＞坡度＞温度＞海拔＞人均GDP＞人口密度（表3）。其中NDVI解释力最大，且远高于其他指标，q值为0.59；其次为降水，q值为0.48；坡度次之，q值为0.18；海拔、人均GDP和人口密度的影响最小，q值均为0.09。这表明植被覆盖度和年降水量变化对2001—2020年黄土高原生态脆弱性变化的解释程度相对较大，其他因子对生态脆弱性变化的解释程度相对较低。

表3 黄土高原地区各影响因子对生态脆弱性的解释力Table 3 Explanatory power of different drivers on ecological vulnerability in the Loess Plateau in 2001-2020

生态脆弱性因子交互作用值均大于单个因子的最大值（表 4），这说明各指标对于黄土高原生态脆弱性的影响并不是相互独立，而是相互作用的结果。海拔与降水、海拔与人均GDP、海拔与人口密度、坡度与降水和坡度与温度之间表现为非线性增强，其余都表现为双因子增强。虽然单因子探测结果显示海拔对黄土高原地区生态脆弱性影响较小，但是当它与NDVI交互作用时，对生态脆弱性具有最强的解释力（q=0.66）。综上，NDVI和降水分别与其他因子交互作用后产生的解释力都远高于其他因子之间的相互作用（q＞0.51）。

表4 黄土高原地区生态脆弱性各影响因子间的交互作用Table 4 Interactions among drivers of ecological vulnerability in the Loess Plateau

3 讨论

本研究表明，2001—2020年黄土高原地区生态脆弱性水平整体较高，且呈现出明显的西北高而东南低的分布格局，该结果与张良侠等（2022）基于“暴露-敏感性-适应力”评估模型框架，结合层次分析和空间主成分分析方法得出的黄土高原生态脆弱性评估结果相似，且与刘会军等（2015）在全国尺度上界定的生态脆弱区相一致。极度脆弱区主要集中在黄土高原西南部的高塬沟壑区，该地区气候干旱且植被覆盖稀疏，水土流失严重，生态脆弱性程度最高，与陈枫等（2018）和苏胜亮等（2022）的研究结果相一致。重度和部分极度脆弱区主要分布在陕北、晋北和鄂尔多斯东部地区，其中包含一些重要的煤炭基地。该区位于干旱和半干旱地带，煤炭能源的开采不仅破坏了地表原本稀疏的植被且进一步加剧了水资源的短缺，从而导致生态环境高度脆弱，这与王佳宁等（2020）的研究结果相一致。中度脆弱区分布在黄土高原中北部的毛乌素沙地和库布齐沙漠地区，该地区沙化草原面积占比大，生态系统对气候变化敏感且抵御外界干扰能力差（Sun et al.，2021）。然而，本研究结果与部分前人研究存在一些不同，例如杨雯娜等（2021）发现黄河流域中的极度脆弱区集中在宁夏北部以及杨晴青等（2019）发现陕西佳县的生态脆弱性为中等脆弱，产生差异的原因可能是受研究尺度、指标选取和评价方法等因素的影响，生态脆弱性评价结果是相对性的，区域之间不具有可比性，这也说明了在区域尺度上建立系统客观的评价体系的重要性。

本研究发现黄土高原地区生态脆弱性主要的两个控制因子为植被覆盖度和降水，这与张良侠等（2022）和王丽霞等（2021）在黄土高原区的研究结果具有一致性。降水是影响干旱和半干旱区植被生长发育的重要因素，植被覆盖度高意味着其生态功能较强，二者的增加客观上都有利于生态脆弱性的降低（Li et al.，2021；陈玉兰等，2022）。例如，在黄土高原东南部，较高的年降水量适宜植被生长，生态脆弱性等级整体较低。而在中部半干旱区，植被覆盖度较低，生态系统较脆弱。研究时段内黄土高原地区年降水量整体呈增加趋势（刘洋洋等，2019），同时大规模实施的生态恢复工程显著提高了植被覆盖度（Xiao，2014；杨丹等，2022），这两者均会降低黄土高原地区的生态脆弱性。然而，何敏等（2019）发现西南地区生态脆弱性主导因子为温度和海拔，李路等（2021）发现气温、地形和植被覆盖度是新疆喀什地区生态脆弱性的主要控制因子，张学渊等（2021）发现西北干旱区生态脆弱性受土壤和地形影响最大，产生分歧的原因主要在于所选研究区不同，自然条件存在差异。相较于西南地区或者新疆喀什地区，黄土高原所在区域整体热量条件的变化以及生态系统受地形的制约程度相对较低，因而温度和海拔对于黄土高原生态脆弱性的影响较小。本研究表明，海拔与植被覆盖度相互作用对黄土高原生态脆弱性具有最大的解释力。

相较于以往的研究（张学渊等，2021；朱琪等，2021；金丽娟等，2022），本研究采用的将NPP作为单一指标的生态脆弱性评价方法，回避了综合指标体系中的指标选取具有重叠性和权重分配主观性高等问题，能够有效评估黄土高原地区生态脆弱性的空间分布格局。基于该方法得到的结果与前人研究结果相一致，证明本研究方法是科学可行的。地理探测器的引入也为定量识别生态脆弱性的控制因子提供了评估方法。然而，本研究仍存在一些不足。首先，本研究基于2001—2020年NPP的动态变化特征来计算生态脆弱性，结果只能表征研究区生态脆弱性的空间变异特征，无法反映其时间变异特征。其次，受数据限制，本研究时段相对较短，且所选生态脆弱性驱动因子没有考虑其他人为活动（如生态工程、耕作活动和牲畜数量），退耕还林和退牧还草等生态政策对黄土高原脆弱性的影响还有待更深入的研究。未来还需进一步完善生态脆弱性的计算方法，利用长时序高空间分辨率数据，并结合研究区的社会经济发展状况，选取系统综合的指标进行生态脆弱性的主控因子分析，以期为黄土高原生态保护和建设提供具体的科学参考和理论支持。