滑坡灾害具有隐蔽性和突发性等特点, 不仅威胁人类生命财产安全, 还会破坏区域生态系统的结构、功能、安全和健康, 从而影响社会经济、生态文明以及人类发展的自然基础^[1]。生态风险评估是一种估计危险环境下可能发生的不利生态影响的方法^[2—3]。量化生态风险为政府制定风险防范、生态系统管理和可持续发展政策提供科学依据。在区域生态风险评估中, 生态扰动的原因通常包括污染物排放^[4]、土地利用变化^[5—6]等人类活动, 但也可能包括自然灾害如地质灾害^[1]、干旱^[7]和洪水^[8]等。其中, 地质灾害是生态风险的重要风险源之一, 已成为国内外研究热点^[9—10]。

对滑坡灾害的生态风险评估, 尽管也会考虑灾害对人类社会的影响, 但通常生态系统是风险的主要承灾体^[11]。目前, 相关的生态风险评估模型主要包括相对风险模型^[12]、概率损失模型^[13—14]和生态风险景观指数^[15]等。其中, 概率损失模型因其有效性和效率而被广泛应用^[16]。该方法通过考虑地质灾害发生的概率、对地质灾害的响应以及承灾体的潜在损失来构建二维或三维模型^{[10—11, 17]}。“危险性-易损性”二维模型重点关注滑坡发生的概率及其对生态系统可能造成的扰动, 缺乏对生态系统的暴露响应表征^[11]。“危险性-脆弱性-潜在损失”三维模型可体现信息间的交互状态, 增强风险评价时效性^[18]。许嘉慧等^[19]同样基于该三维评价框架, 进行格网、行政、子流域多尺度下的滑坡灾害生态风险评价。基于三维评价框架, 可为生态风险评价与管理提供更为系统全面的参考。然而, 在潜在损失方面, 学者通常关注生态系统服务的损失, 这些损失可以基于多种模型或指数来量化, 但对直接暴露于灾害环境的人群的关注仍然不足^[20]。

长江上游地区的地质构造活动强烈、地形地貌复杂、气候条件多变, 是滑坡等地质灾害的高风险区域。尽管此前多项研究评估了长江上游地区地质灾害的危险性^[21—22], 但这些研究主要集中在灾害孕育环境的危险性, 鲜有研究关注滑坡等地质灾害生态系统风险的综合评估。本文结合遥感和GIS技术, 建立了基于滑坡灾害、生态脆弱性和潜在损失的评估技术框架, 对长江上游滑坡灾害造成的生态风险进行评估, 为政府部门提供防灾减灾科学依据。

1 数据与方法 1.1 研究区概况

长江上游地区包括我国西南地区(90°32′—111°27′E, 24°27′—35°45′N之间), 具体包括四川省、重庆市、云南省、贵州省、湖北省、甘肃省等部分地区, 面积约100万km²(图 1)。该地区横跨我国一、二级地貌阶梯, 地形复杂, 地势西高东低, 西部以草原为主(平均海拔超过3000 m), 东部则以四川盆地平原和丘陵地带为主, 分布灌溉农田、灌木、阔叶林^[23]。该地区年平均降水量约1000 mm, 地区差异较大, 降水量自东南向西北递减。多年平均气温介于16—18 ℃之间, 在空间上也呈现出东南高、西北低的特征^[24]。主要为高原山地地貌, 地形破碎、山高坡陡, 新构造活动强烈和地震频发, 加之降雨量充沛, 导致其成为我国滑坡等地质灾害最为严重的地区^[25]。

1.2 数据源

本研究使用的数据包括数字高程模型(Digital Elevation Model, DEM)、土地利用、滑坡清单、归一化植被指数(Normalized Difference Vegetation Index, NDVI)、降水、河流、岩性、断层、人口、净初级生产力(Net Primary Production, NPP)、植被蒸散发(Evapotranspiration, ET)和土壤数据(表 1)。所有空间数据均采用三次卷积法重采样至1 km分辨率, 投影统一为正轴等积割圆锥投影。

1.3 研究方法 1.3.1 研究框架

滑坡灾害生态风险可被定义为灾害发生的概率与其发生后对生态系统和人类社会造成的损失的乘积^[30]。本文基于概率损失框架提出了一个集成滑坡危险性、生态脆弱性和潜在损失的评估框架(图 2), 生态风险根据滑坡危险性、生态脆弱性和潜在损失进行量化, 评价结果按下式计算：“生态风险=危险性×脆弱性×潜在损失”, 三个指标的具体释义见表 2。考虑到三个指标无统一量纲, 需要先对它们进行归一化处理, 然后再相乘。

1.3.2 基于随机森林的滑坡危险性评价

滑坡灾害危险性评价总体方法流程图如图 2所示。首先制作研究区完整的滑坡清单及风险评价指标体系, 并建立相应的空间数据库。其次, 利用滑坡的训练集提取每个条件因子的数值作为随机森林模型的输入特征。第三步, 利用经过训练的模型生成滑坡危险性地图。此外, 通过受试者工作特征(Receiver Operating Characteristic, ROC)曲线来评估模型的准确性。

随机森林(Random Forest, RF)是一种利用多个决策树进行分类或回归的机器学习方法, 能够提供贝叶斯分类器的近似结果^[30]。与其他模型不同的是, 随机森林能够提供多种变量重要性的度量, 其中最可靠的度量是基于树中的变量值被随机置换时分类精度的降低^[31]。在解决分类问题时, 随机森林的预测被认为是未加权的多数类别投票^[32]。bagging技术用于从变量的随机样本中选择训练数据集以进行模型校准。对于每个变量, 如果该变量的值在袋外观测值中进行排序, 则该函数将确定模型预测误差^[33]。

在滑坡危险性评价中, 随机森林模型利用树之间的差异性, 让每棵树对类别进行投票, 并根据多数票确定最终的类别。这种集成方法在复杂的数据集上表现出稳健和准确的性能, 几乎不需要调整参数, 而且可以容忍许多噪声变量^[34]。随机森林模型采用基于R的数据分析工具包“randomForest”库, 通过网格搜索进行参数优化, 并使用五折交叉验证来评估最终精度。在本研究中, 树的数量设置为100, 特征的数量设置为4。

1.3.3 基于景观格局指数的生态脆弱性评估

生态系统脆弱性反映了生态系统因外部干扰、结构、组成和其他生态特征而发生变化的可能性。景观格局反映了人类对自然生态系统的影响方式和程度, 这些生态影响具有区域性和累积性, 能够通过生态系统的结构和组成来反映。景观格局指数(Landscape Pattern Index, LPI)量化了景观异质性与生态过程之间的相互作用, 适用于更大空间尺度的脆弱性测绘^[35]。

本研究参考了景观指数选择的基本原则和技术方法, 选取了斑块密度(Patch Density, PD)、景观分离指数(DIVISION)和景观干扰指数(Landscape disturbance index, LDI)^[36], 以构建综合指数衡量区域生态系统在地质灾害干扰下的脆弱性。对三个景观指数进行计算和标准化, 得到其平均值, 具体公式如下所示：

(1)

其中, 综合脆弱性指数(Vulnerability Index, VI)表示区域生态系统的脆弱程度, 其数值越高, 说明生态系统的稳定性和对外界扰动的抵抗能力越差。PD代表景观生态系统的破碎程度, 破碎程度越高, 生态系统对干扰的恢复能力就越差。DIVISION是揭示景观类型中斑块个体分离程度的指标。PD和DIVISION两个指标可以通过专业景观指数软件Fragstats 4.3进行计算。LDI能够反映不同景观生态系统的干扰程度, 可以根据景观破碎度指数(C_i)、S_i和景观优势度指数(D_i)计算得出。

(2)

(3)

(4)

(5)

式中, a、b和c的权重分别被指定为0.5、0.3和0.2；N_i是第i个景观的斑块数量；A_i是第i个景观的面积(km²)；A为所有景观的总面积。

在计算景观生态风险时, 景观分析单元的大小应为平均斑块大小的2—5倍^[37]。因此, 将研究区划分为1 km × 1 km网格, 共991163个网格单元, 作为脆弱性评价单元。

1.3.4 潜在损失计算

生态风险的最终承受者不仅包括人类, 也涉及生态系统的组成结构。本研究将生态系统和暴露人群视为生态风险的受体, 重点关注生态系统服务, 如水土保持、水源涵养和生物多样性保护等。特定的生态系统格局维持着这些生态服务, 一旦被破坏, 服务功能将降低, 增加生态风险。

考虑到研究区的生态环境状况, 选取水源涵养、水土保持和植被NPP三种生态系统服务来代表潜在的生态损失。本研究对这三种生态服务系统进行标准化处理, 然后在空间上叠加, 作为生态系统受到地质灾害胁迫后的潜在损失。此外, 利用人口分布数据来表征地质灾害发生后的人口暴露程度。通过对潜在生态系统损失和人口暴露进行标准化处理, 借助ArcGIS软件中的空间叠加工具得到最终的潜在损失。

(1) 水源涵养

水源涵养的计算基于水平衡法^[38]：

(6)

式中, Q_wr代表节水量(m³/A)；P_i代表降水量(mm/a)；R_i代表地表径流量(mm/a)；ET_i代表蒸散量(mm/a)；i代表生态系统服务类型。

(2) 水土保持

水土保持的计算基于修正通用土壤流失方程(Revised Universal Soil Loss Equation, RUSLE), 通过计算潜在水土流失与实际水土流失之间的差值来确定水土保持量^[39]。

(7)

(8)

(9)

式中, SOR为土壤保持量(t hm^-2 a^-1)；SE_p代表潜在的侵蚀；SE_r代表实际侵蚀量；R代表降雨侵蚀力因子(MJ mm hm^-2 h^-1 a^-1)；K为土壤侵蚀因子(t h MJ^-1 mm^-1), LS为地形因子, 可以通过DEM计算^[40]；C是植被因子；P是管理因子, 可参考以前的研究^[41]对不同的土地利用进行赋值。

2 结果与分析 2.1 滑坡危险性评估

绘制完整的滑坡清单是进行滑坡危险性评估的基础。通过收集和整理各省地质环境监测总站和中国科学院资源环境科学与数据中心的滑坡数据, 本文获取了长江上游地区包括59212处滑坡的完整滑坡清单, 并使用核密度分析工具对其进行可视化(图 3)。根据密度分析结果, 研究区滑坡主要分布在龙门山断裂带(四川省广元市-成都市西北部)、大渡河峡谷一带和长江三峡库区一带。这些地区具有复杂的地质结构、活跃的地壳运动、频繁的地震活动、严重的河流侵蚀以及集中且高降雨量等特点, 为大量滑坡的形成和诱发提供了环境。

研究选取了59212个滑坡地质灾害作为正样本, 并在滑坡点500 m缓冲区外随机选择了等量的点作为负样本。在所有样本中, 随机选取70%的点用于模型训练, 其余30%的点用于验证。

滑坡的发生是各种影响因素综合作用的结果^[42]。选择合适的影响因素是建立滑坡评价模型的关键, 但目前还没有统一的标准和方法^[43]。根据地质环境特征、前人文献的分析以及数据可获取性, 确定了9个滑坡影响因子：坡度、坡向、高程、距河流距离、距断层距离、NDVI、工程地质岩组和降水。利用ArcGIS软件, 基于DEM数据中生成了坡度、坡向和剖面曲率三个评价因子, 并制作了河流和断层的缓冲分析图(图 4)。

在滑坡危险性评价中, 因子的选取是关键步骤, 既要考虑各因子对滑坡的影响程度, 又要避免各因子之间的强相关性。为评估因子之间的相关性, 进行了皮尔逊相关系数分析。图 5显示了各因子间的相关系数, 除降水和高程外, 所有因子间的相关系数都小于0.8, 大部分因子间的相关系数小于0.3。这表明选取的9个因子之间没有较强的相关性, 适合用于滑坡危险性评价。

在滑坡危险性评估中, 采用自然断点法将滑坡敏感性指数划分为五个等级：极低危险性、低危险性、中危险性、高危险性、极高危险性(图 6)。其中, 高和极高危险区占总面积的25%, 主要集中在四川省雅安市、四川省巴中市至达州市一带以及三峡库区附近。中危险区占17%, 主要分布在贵州省北部, 低危险区占14%, 极低危险区占44%。模型的分类性能通过ROC曲线评估^[44], 曲线下面积(Area Under the Curve, AUC)为0.86, 远高于0.5, 说明模型具有很高的准确性。

2.2 生态脆弱性和潜在损失空间分布特征 2.2.1 生态脆弱性空间分布特征

根据1.3.3提出的生态脆弱性评估方法, 对各单元进行了脆弱性值的计算。然后, 基于自然断点法将脆弱性值划分为五个等级：极低、低、中、高、极高。生态脆弱性的空间分布如图 7所示, 极高和高等级脆弱区占总面积的40%, 中等级脆弱区占总面积的28%, 低和极低等级脆弱区占总面积的32%。高生态脆弱地区主要集中在四川盆地的平原和丘陵地区、云贵高原、三峡库区以及黄土高原, 这些地区土地利用强度高、景观破碎化、生态系统稳定性受到威胁, 导致生态系统高度脆弱。

2.2.2 潜在损失空间分布特征

本研究通过将生态损失层和人口暴露层进行叠加来计算潜在损失, 所有层均采用自然断点法划分为五个级别(极低、低、中、高、极高), 其空间分布特征如图 8所示。研究区的潜在生态损失主要为中、低度。生态损失值最高的区域集中在东部的三峡库区和西南部的横断山区, 这些地区具备良好的水热条件, 植被茂盛, 净初级生产力强劲, 生态系统服务功能较强。然而, 西北地区海拔较高, 气候寒冷, 植被覆盖度低, 生态系统服务能力弱, 因而生态损失价值较低(图 8)。在人口暴露方面, 空间分布差异较大：高暴露和极高暴露区域集中在成渝地区双城经济圈, 该区域人口密集, 经济发达, 因此在滑坡灾害发生时面临较高的人口暴露风险；相比之下, 其他地区人口较少, 暴露风险相对较低(图 8)。

从潜在损失总额来看, 等级较高的地区主要分布在三峡库区和横断山区, 这些地区拥有较高的植被覆盖度, 生态系统服务功能强大(图 8)。然而, 这些地区在面临滑坡灾害时潜在损失巨大, 丰富的生态系统服务和高植被覆盖度使它们更容易受到生态破坏和损失, 因此需要高度关注和防范。

2.3 滑坡灾害生态风险评估 2.3.1 生态风险空间格局

根据1.3.1的计算方法, 得到了在网格尺度下的长江上游地区滑坡灾害生态风险。采用自然断点法将滑坡灾害生态风险划分为五个等级：极高风险、高风险、中风险、低风险和极低风险(图 9)。

长江上游地区生态风险等级高和极高的区域比例较低, 分别占总面积的2.83%和0.38%(表 3), 总体上呈现局部聚集、全局分散的分布格局。在具体分布上, 研究区的中西部和东部均属于高风险级别的区域, 形成了较为突出的分布特点。聚集性较高的高和极高风险区主要分布在岷江、大渡河、长江流域地区, 这些区域受地质环境和河流侵蚀影响较大。中风险区域主要位于高风险区域的外围以及丘陵地带, 低风险区主要为研究区的西部地区。总体而言, 生态风险的空间分布呈现出明显的地域性特点, 反映出地形地貌、生态环境和人类活动等多方面的综合作用。

2.3.2 生态风险影响机制分析

滑坡灾害生态风险成因分析对于风险防控至关重要。本文以两个生态风险较高的地区为典型案例, 结合生态风险评估结果和自然条件, 分析其对生态风险的影响机制。高风险区域(图 9)包括：(1)四川省雅安市-四川省乐山市一带；(2)重庆市万州市-湖北省宜昌市一带。

(1) 四川省雅安市-四川省乐山市一带：该区域位于四川盆地与横断山区的交界地带, 地形复杂多样, 山地和丘陵占据主导。区域内山脉纵横, 新构造活动强烈, 复杂的地质背景和脆弱的地质环境使滑坡灾害频繁发生^[45]。该区域的植被覆盖率较高, 这有助于维持生态系统的稳定性。然而, 由于土地利用复杂, 涉及多种类型的土地利用形式, 如农业、林业和城市建设等, 导致区域内的生态系统面临较高的脆弱性。根据评估结果, 该区域的滑坡危险性、生态脆弱性和潜在生态损失都比较高, 这三个因素共同导致了较高的生态风险。

(2) 重庆市万州市-湖北省宜昌市一带：该区域位于三峡库区, 地形多为山地和丘陵, 地质条件复杂。季节性降雨和急剧变化的水库蓄水对土壤和岩层的稳定性产生显著影响, 导致滑坡等地质灾害频繁发生^[46]。此外, 万州市和宜昌市经济较为发达, 人口密集, 人类活动强度高, 对生态环境造成明显扰动, 从而加剧了该区域的生态风险。

总体而言, 研究区高生态风险的成因和机制具有显著的空间异质性。除了孕育灾害的客观环境外, 包括农业和工程在内的人类活动也在高生态风险中发挥了重要作用。

3 讨论 3.1 滑坡危险性评价因子选取

滑坡危险性评价因子选取至关重要, 在图 5中的评价因子之间的皮尔逊相关系数分析中, 降水和高程的相关系数超过了0.8, 这表明这两个因子之间存在较强的相关性, 应将降水评价因子剔除。然而, 本文所使用的降水数据集是基于薄板样条插值法对全球气象站点进行插值生成的, 并使用海拔作为协方差。因此, 在本研究中, 高程和降水因子显示出强烈的相关性, 这主要是因为两者在数据生成过程中存在一定的交互作用和依赖关系。

为进一步检验评价因子的协方差, 使用容差和方差膨胀因子(Variance inflation factor, VIF)进行了分析。当方VIF大于或等于5或容差小于或等于0.2时, 该指标被视为存在协方差。如表 4所示, 高程和降水的VIF值分别为3.8和3.6, 均未达到5, 表明这两个因子之间的协方差尚未达到显著水平, 因此两者都可以参与危险性评价。

3.2 风险控制与展望

不同高风险地区影响灾害发生的因素有所不同, 因此提出的风险控制措施应更加针对性和因地制宜。未来的政策制定应同时兼顾经济发展和生态保护。对于滑坡危险性、生态脆弱性、潜在损失较大的地区, 应提高地质灾害监测频次, 并在部分地区设立禁止开发区。对于生态脆弱性较高的地区, 应制定相应的保护措施, 如对不稳定边坡进行加固、植被恢复等。此外, 对于潜在生态损失较大的地区, 应以保护为主, 尽量避免过度的土地开发和利用。

目前的研究存在一定的局限性, 在未来的研究工作中应认识到这些局限性并寻求解决之道。本文提出的“危险性-脆弱性-潜在损失”滑坡灾害生态风险评价框架在计算滑坡造成的潜在损失时引入了人口暴露这一概念, 充分考虑了滑坡灾害对人口的威胁。但是, 在生态系统方面, 潜在损失不仅仅包含水土保持、水源涵养和生物多样性保护这三类生态系统服务。由于详细空间数据的缺乏, 本研究的潜在损失计算仅限于这三类生态系统服务。此外, 本研究中采用的风险评估方法可以通过纳入不确定性和敏感性分析来改进, 以更好地理解结果的空间异质性。在未来的研究中, 应考虑更广泛的生态系统服务, 以改进潜在损失的计算, 进行更全面的不确定性和敏感性分析, 以更好地了解评估结果的有效性。

4 结论

(1) 基于9个滑坡条件因素, 使用随机森林模型对长江上游地区进行了滑坡危险性评价, 区域内滑坡危险性等级主要为极低和低, 面积占比分别为44%和14%。较高和极高危险性区域分别占总面积的17%和8%, 集中在四川省雅安市、四川省巴中市至达州市一带以及三峡库区。

(2) 极高和高脆弱区占总面积的40%, 中脆弱区占总面积的28%, 低和极低脆弱区占总面积的32%, 四川盆地、云贵高原、三峡库区和黄土高原的平原和丘陵地区存在较高的生态脆弱性, 生态系统稳定性受到威胁。潜在损失较高区域主要为三峡库区和横断山区, 这些地区具有强大的生态系统服务功能, 但在发生地质灾害时, 潜在损失较高。

(3) 研究区域内的生态风险总体较低, 较高的生态风险区域呈现局部聚集, 占总面积的3.21%, 主要集中在岷江、大渡河、长江流域地区。高生态风险的成因和机制存在显著的空间异质性, 人类活动显著影响生态风险。