受日益加剧的气候变化(如气候变暖和极端气候事件)和人类活动(如城市化和耕地扩张)的影响, 生态系统退化已成为全球范围内最严重的生态环境问题之一^[1—2]。生态脆弱性作为表征生态系统质量的一个关键指标, 能综合反映一个生态系统面对外界干扰时的承受程度, 以及适应新环境的能力^[3]。定量评估生态脆弱性及其驱动机制, 可加深对全球变化背景下生态系统退化风险的认识, 进而指导区域生态系统的合理开发与保护, 因此成为近年来全球变化和可持续发展研究的热点^[4—8]。

生态脆弱性评估主要包括综合指标法和单一指标法两种评价方法。综合指标法主要依据“暴露-敏感性-适应力”^[9—10]、“压力-状态-响应”^[11]或“压力-敏感性-恢复力”^[12]等评估框架, 从资源、环境、社会和经济等方面选取多个影响因子构建评价指标体系, 进而采用层次分析(Analytic Hierarchy Processes, AHP)、主成分分析(Principal Component Analysis, PCA)或熵权法等加权计算得到生态脆弱性^{[11, 13—14]}。该方法实际是通过各因子某一年份综合状态来反映当年的生态脆弱性状况。由于指标选取具有明确的脆弱性内涵^[15—16], 综合指标法被大多数研究所青睐。但该方法存在指标数据获取困难、指标间相关性高和指标权重分配主观性强等问题^[17—20]。对比之下, 单一指标法以能够综合反映研究区生态脆弱性的单一因子为评价指标^{[3, 21]}, 通过指标的多年波动状况^[22]表征生态脆弱性。该方法具有数据依赖性较低和计算相对简单的优势^{[17, 19, 22]}, 但可能存在指标代表性不足的问题。由于两种方法的指标体系和计算原理完全不同且各具优势, 厘清二者评估结果的一致性是深入理解生态脆弱性格局及其不确定性的关键。然而, 现有研究大多选择一种方法对某个区域的生态脆弱性进行评估^{[3, 9—14, 17—19]}, 尚缺乏针对两种方法评估结果一致性的定量分析。

黄土高原是世界上面积最大的黄土堆积区, 也是我国典型生态环境脆弱区和水土保持的重点防治区^[23—25]。在气候变化(如干旱)和长期不合理人类活动(如过度开垦和放牧)的影响下, 水土流失和土地资源退化严重, 已严重威胁区域生态环境和社会经济可持续发展^[24—27]。本研究以黄土高原地区为例, 基于依据“暴露-敏感性-适应力”模型框架的综合评价法和以生态系统净初级生产力(Net Primary Production, NPP)为指标的单一评价法, 对比分析了该区生态脆弱性空间分布格局的差异, 并结合地理探测器探讨了两种方法所估算的脆弱性与驱动因子关系的一致性。本研究旨在探寻不同评估方法下的生态脆弱性分布格局及其与驱动因子关系的一致性, 研究结果不仅可以加深对脆弱性评估结果不确定性的认识, 还可以为后续研究方法的合理选择提供参考。

1 研究区概况

黄土高原(33°43′—41°16′N, 100°54′—114°33′E)位于黄河流域中部, 包括太行山以西、日月山以东、秦岭以北、鄂尔多斯高原以南的广大地区(图 1), 是我国乃至世界上水土流失最严重、生态环境最脆弱的地区。行政区域跨青海省、甘肃省、宁夏回族自治区、内蒙古自治区、陕西省、山西省及河南省七省区, 东西绵延约1000km, 南北地跨750km左右, 总面积约64万km²。黄土高原地区地势西北高、东南低, 区域海拔为84—5206m。该区属暖温带大陆性季风气候, 年均气温3.6—14.3℃, 年均降水量300—800mm。降水量在空间分布上自西北向东南逐渐增多, 在季节动态上主要集中在夏季, 7—9月降水量占年降水量的60%—70%。植被由东南向西北可划分为森林带、森林-草原带、草原带和草原-荒漠带。该区基本上覆盖了100—300m厚的细腻黄土微粒, 土质疏松, 土壤侵蚀强烈, 水土流失面积高达45万km²。

2 研究方法 2.1 数据来源与预处理

本研究主要关注不同方法估算生态脆弱性分布格局差异, 涉及的数据类型多, 因此本研究主要采用已公开发布的或经模型运算得到的网格化数据产品, 未采用数量有限的野外实测数据。具体数据主要包括遥感、气象、社会经济和其他四类。遥感数据包括2000—2018年净初级生产力(Net Primary Production, NPP)和归一化植被指数(Normalized Difference Vegetation Index, NDVI), 分别由美国地质调查局网站(https://lpdaac.usgs.govg11)提供的MOD17A3H产品(年尺度, 500m)和MOD13A3产品(月尺度, 1km)获取。气象数据主要包括2000—2018年月均气温、月均相对湿度和年降水量, 基于研究区及其周围148个气象站点的观测数据, 利用ANUSPLIN插值软件, 结合数字高程模型(数据来源于美国奋进号航天飞机的雷达地形测绘SRTM), 采用局部薄盘光滑样条函数插值获取得到。湿润度指数采用伊万诺夫湿润度计算公式, 由年降水量与年蒸发力的比值计算得到^[28]。社会经济数据包括2000—2018年GDP、人口密度和路网密度。其中GDP和人口密度(空间分辨率为1km)由中国科学院资源环境科学数据中心(http://www.resdc.cn/)下载获取；路网密度基于对中国交通地图册矢量化获取得到的道路空间分布数据, 由县域内道路长度与县域面积的比值计算得到。其他数据包括坡度、起伏度、水蚀模数、风蚀模数、水源涵养量和香农多样性指数。坡度和起伏度由数字高程模型计算得到。2000—2018年水蚀模数、风蚀模数和水源涵养量分别采用修正的通用土壤流失方程(RUSLE)^[29]、修正风蚀方程(RWEQ)^[30]和降水贮存量法^[31]估算。香农多样性指数(Shannon′s diversity index, SHDI)基于研究区2000—2018年土地利用数据(由中国科学院资源环境科学数据中心下载获取, 空间分辨率为1km)中的土地利用类型信息, 利用Fragstats 4.2软件计算得到, 算式如下：

(1)

式中, P_i为第i种土地利用类型所占的面积比例；n为土地利用类型数量。

本研究将所有指标数据的空间分辨率统一重采样为1km, 再进行后续分析。

2.2 综合指标评价法和单一指标评价法

综合指标评价法包括指标体系构建、指标权重确定和生态脆弱性指数计算三个步骤。首先, 基于“暴露-敏感性-适应力”的生态脆弱性概念与分析框架, 结合黄土高原地区气候干旱和土壤侵蚀严重等主要生态问题, 同时考虑数据的时效性和可获取性, 选取了13个指标(表 1)。在综合评价指标体系中, 暴露反映了生态系统所处环境中外界干扰的程度, 包括年降水量、湿润度指数、大于10℃积温、水蚀模数、风蚀模数和人口密度6个指标；敏感性反映了生态系统容易受到外界干扰胁迫的程度, 包括年最大NDVI、坡度、起伏度、香农多样性指数和水源涵养量5个指标；适应能力反映了生态系统面对外界干扰时的内在调节能力或恢复潜力, 包括人均GDP和公路密度2个指标。

其次, 采用AHP计算各评价指标的权重。首先参考相关研究^[32—33]建立生态脆弱性评价准则层的判断矩阵, 进而计算得到暴露、敏感性和适应力的权重分别为0.43、0.43和0.14。然后将表征暴露、敏感性和适应力的所有指标项标准化后乘以该准则层的权重。其中, 正相关指标和负相关指标采取不同的标准化处理方法, 以保证处理之后的指标值越大反映该指标对生态脆弱性影响越显著。

最后, 基于PCA计算得到的累计贡献率大于85%的前4个主成分结果, 估算生态脆弱性指数(Ecological Vulnerability Index, EVI)：

(2)

式中, Y和r分别表示各主成分及对应的贡献率, EVI的值越高, 生态系统越脆弱。

单一指标法通过计算研究时段内NPP的年际波动以及年际变异的变化趋势来估算生态脆弱性^[34]：

(3)

(4)

(5)

式中, S代表敏感性, 反映了生态系统受外界扰动的影响, 由NPP的年际波动情况计算得到；A代表适应性, 表征了生态系统在面对外界干扰时的内在调节和恢复能力, 由NPP年际变率线性拟合趋势线的斜率计算得到；F_i表示黄土高原地区第i年的NPP值；F表示该地区研究时段内NPP的平均值；x为自然数1、2、3…, 表示研究时段内的年份时间序列；y为NPP年际变率的绝对变化量, 即|F_i-F|。由公式(4)和(5)计算得到的S和A, 经极差标准化后(统一量纲), 进而输入公式(3)计算V。

由于仅2000年、2005年、2010年、2015年和2018年具备综合指标评价法所需的全部指标因子数据, 为了确保两种方法评估结果的可比性, 单一指标法研究时段选为2000—2018年, 而综合指标法所估算的脆弱性分布格局利用五个年份脆弱性指数的平均值表征。

2.3 生态脆弱性分级

本研究采用自然断点分级法对生态脆弱性结果进行分级。自然断点分级法主要依据数值的统计分布规律, 结合数据直方图和标准差进行数据分级^[35]。相比其他分级方法, 该方法能使各级内部方差之和最小且各级之间差异最大化, 较客观地划分生态脆弱性等级。本研究将生态脆弱性划分为微度脆弱、轻度脆弱、中度脆弱、重度脆弱和极度脆弱5级。其中, 综合指标法2000年、2005年、2010年、2015年和2018年结果采用相同的生态脆弱性指数分级阈值, 以保证五期结果之间的可比性。

2.4 驱动因子分析

本文采用地理探测器^[36]定量分析了基于两种方法所估算的生态脆弱性空间格局与各影响因子关系的一致性。各因子影响力的高低, 通过地理探测器计算得到的q统计值来度量。q值越大说明驱动因子对于生态脆弱性空间分布的解释力越高, 反之越低。首先, 选取EVI作为因变量, 选取综合指标评价体系中的13个评价指标作为自变量因子, 在15km×15km格网尺度, 统计EVI和各自变量因子平均值。其次, 利用自然断点分级法将各自变量因子分为5级, 由数值量转变为类型量。最后, 利用地理探测器计算q统计值。

3 结果 3.1 基于单一与综合指标法的黄土高原地区生态脆弱性空间格局差异

整体上, 单一指标法与综合指标法得到的生态脆弱性空间格局一致(图 2)。两种方法均表明黄土高原地区生态脆弱性处于较高水平, 中度及以上脆弱等级面积占总面积的58%—63%, 整体呈现出从东南到西北递增的趋势。其中, 东南部地区, 如渭河和黄河流域构成的河谷平原区以及太行山区主要为微度或轻度脆弱分布区；西北部地区, 如库布齐沙漠和毛乌素沙地地区主要为重度或极度脆弱区。

从像元尺度来看, 尽管单一指标法与综合指标法得到的生态脆弱性指数相关性较高(r=0.80, P < 0.001), 但具体数值和等级划分差异较大(图 3)。单一指标法估算的生态脆弱性指数数值范围为-0.8—1.2, 综合指标法结果的数值范围为0.4—1.0。两种方法生态脆弱性等级结果不一致区域占全区总面积的45%, 主要表现为相邻等级间差异(占全区面积的39.61%), 且单一指标法得到的生态脆弱性等级比综合指标法结果高一个等级的区域占比相对较大(占全区总面积的24.39%)(图 3)。在空间分布上, 两种方法等级划分结果不一致区域主要分布在植被覆盖较好区域, 包括东部边缘的太行山区、东部的吕梁山区、东南部的渭河和汾河河谷平原区、西北边缘的黄河农灌区以及西部的洮河流域(图 2)。

3.2 基于单一与综合指标法的黄土高原地区生态脆弱性控制因子差异

地理探测器结果显示, 单一和综合指标法获得的生态脆弱性分布与驱动因子的关系整体一致(图 4)。就单因子而言, 解释力较高的均为年降水量、年最大NDVI和湿润度指数, 单一指标法相应的q值为0.59—0.71, 综合指标法相应的q值为0.62—0.84；其次为风蚀模数、水源涵养和坡度, 两种方法的q均值在0.29—0.44之间；相对而言, 水蚀模数、大于10℃积温、起伏度和香农多样性指数的解释力较弱, 两种方法的q均值都小于0.09。同时, 不同控制因子之间存在明显的双因子增强交互作用和非线性增强交互作用。其中, 年降水量和年最大NDVI的双因子增强交互作用具有最高的解释力, 其对单一指标法和综合指标法结果的解释力分别为78%和92%, 平均解释力高达85%。但各单因子对单一和综合指标法评估结果的解释力数值及排序存在一定差异。例如, 对单一指标法结果解释力最高的单因子为年降水量(71%), 其次为年最大NDVI(69%)；对综合指标法结果解释力最高的单因子为年最大NDVI(84%), 其次为年降水量(82%)。

4 讨论

从生态脆弱性程度和空间格局来看, 单一指标法和综合指标法均表明黄土高原地区整体处于较高生态脆弱性水平, 且呈现出从东南到西北递增的空间格局, 这与过去关于该地区的脆弱性评估结果一致^{[27, 32, 37—38]}。这主要是因为研究区大部分地区属于半干旱气候且为黄土堆积区, 植被条件较差, 水土流失严重^{[23, 39]}, 对气候变化(如干旱)和不合理人类活动(如过度开垦和放牧)等外界干扰较为敏感^[40—42], 整体上生态系统比较脆弱。相对而言, 在降水量较高、植被覆盖度较好且水土流失较弱的东南部区域^[43], 生态脆弱性级别较低。然而, 从像元尺度来看, 虽然两种方法估算的生态脆弱性指数相关性较高, 但在植被覆盖较好区域的生态脆弱性等级存在差异。这可能是由于两种方法的原理不同所致。例如, 综合指标法依据所选指标的现状来估算生态脆弱性(表 1), 植被覆盖度较好的区域生态脆弱性级别较低；而单一指标法通过计算植被生产力的年际波动状况来估算生态脆弱性, 波动剧烈的区域生态脆弱性级别较高。研究时段内, 受生态工程和气候变化的影响, 部分植被条件较好区域的植被生产力亦表现出较剧烈的年际波动态势^{[42, 44]}, 可能是导致单一指标法估算生态脆弱性等级偏高的重要原因。

从控制因子来看, 两种方法估算的生态脆弱性均与气候湿润程度(年降水量和湿润度指数)和植被状况(年最大NDVI)的相关性最强, 与其他学者的研究结果相一致^{[27, 32]}。这主要是因为研究区大部分区域属于干旱半干旱区, 加之土壤侵蚀严重, 土壤含水量低, 蓄水工程少, 植被生长对降水比较敏感^{[42, 44]}；同时, 研究区具有连续的第四纪黄土堆积, 容易发生水土流失, 生态系统服务功能受植被条件影响较大^{[39, 45—47]}。在研究区西北部的库布齐沙漠和毛乌素沙地地区, 气候长期持续干旱少雨, 植被稀疏, 地表可蚀性颗粒充足, 且大风频繁发生, 土壤风蚀严重^{[41, 43—44]}, 生态脆弱性程度较高, 属于重度或极度脆弱区。而在研究区东南部的河谷平原和丘陵区, 年降水量相对较高, 为森林和森林-草原分布区, 植被覆盖度较高, 生态系统水土保持功能较好^{[43, 46]}, 为微度或轻度脆弱区。

研究结果表明单一指标法所估算的生态脆弱性程度、空间格局及其与驱动因子的关系与综合指标法均具有较高的相似性。这说明, 单一指标法虽然计算方法简单, 但可以快速有效地实现区域生态脆弱性评估。相对而言, 尽管以往研究大多采用综合指标法, 但该方法需要从资源、环境、社会和经济等多个方面搜集数据, 数据的时效性、可获取性和精度难以得到保证。研究选取的13个指标中, 部分指标(如水蚀模数、风蚀模数和水源涵养量)可获取的数据仅到2018年, 补充计算则需要大量实测数据和遥感数据的支持。综合指标法基于某个特定年份的气象、植被、土壤和社会经济状况进行生态脆弱性评估, 其评估结果可能受所选年份的影响较大。如本研究综合指标法所估算的不同年份生态脆弱性指数数值整体呈先增后降再增再降的波动变化态势(图 5), 五个年份脆弱性等级完全一致的面积仅占全区的5.89%。本研究仅对比了基于“暴露-敏感性-适应力”框架综合指标法与单一指标法的一致性, 基于不同评估框架(如“压力-状态-响应”和“压力-敏感性-恢复力”)和不同指标因子所估算的生态脆弱性格局差异，有待深入研究。

5 结论

采用基于“暴露-敏感性-适应力”模型框架的综合指标法和基于生态系统净初级生产力的单一指标法, 定量分析了黄土高原地区生态脆弱性空间变异特征及其驱动因子的一致性。整体上, 单一和综合指标法评估结果具有较高的一致性, 其结果均表明该地区生态脆弱性程度较高, 且呈现出东南低而西北高的空间格局；但像元尺度上, 生态脆弱性等级划分在近一半区域不一致, 主要表现为单一指标法结果高一个等级, 尤其在东南部、中部及西部植被覆盖较好的区域。从生态脆弱性空间格局的控制因子来看, 两种方法均表明主控因子为年最大NDVI、年降水量和湿润度指数, 其中年降水量和年最大NDVI的双因子增强交互作用解释力最强。综合考虑两种方法的数据需求、计算原理与评估结果的一致性, 单一指标法因其数据依赖性低、计算简单及与综合指标法结果相似性高等优势, 是生态脆弱性空间格局评价的更有效方法。