利用卫星嵌入技术和多维环境因素对怒江中下游地区的滑坡进行测绘与易发性评估（2017–2025年）

《Remote Sensing》：Landslide Mapping and Susceptibility Assessment in the Middle and Lower Reaches of the Nujiang River (2017–2025) Using Satellite Embedding and Multidimensional Environmental Factors

　　**亮点**
- 主要发现是什么？
本研究开发了一个基于卫星嵌入的框架，用于2017年至2025年间怒江中下游地区的年度滑坡清单制图。
滑坡表现出强烈的空间聚集性、显著的年际变化以及持续存在的热点区域，揭示了滑坡活动的明显时空模式。

**主要发现的意义是什么？**
该工作流程提高了在复杂山区进行年度滑坡识别和时空特征分析的效率。
这些发现为区域滑坡监测、热点区域追踪、危险区划和基于风险的管理提供了科学支持。

**摘要**
滑坡制图和易发性评估对于危险识别、基础设施保护和风险管理至关重要。怒江中下游地区地形起伏较大，地貌变化迅速，景观条件脆弱，这增加了滑坡发生的风险并阻碍了及时检测。为了改进滑坡活动的时空特征分析，我们开发了一个多源地球观测框架，用于年度滑坡制图和易发性评估。
首先，生成了年际嵌入变化强度图，以辅助对滑坡相关地表扰动的视觉解释；其次，通过实地验证和视觉解释收集了年度滑坡和非滑坡样本；第三，利用Google Earth Engine上的随机森林生成了2017-2025年的年度10米滑坡地图；最后，将24个多维环境因素纳入滑坡易发性建模中。
滑坡主要集中在怒江沿岸及相邻的高起伏峡谷坡面上，具有显著的年际变化，但热点区域相对稳定。SHAP分析进一步确定BSI_mean是最重要的预测因子，其平均绝对SHAP值为0.116，其次是NDVI_mean和地形相关变量，表明裸露地表、植被状况和地形切割与滑坡发生密切相关。
本研究提供了年度滑坡清单和易发性信息，有助于灾害缓解和基础设施规划。

**1. 引言**
滑坡是指在重力作用下岩石、土壤或松散碎屑向下滑动的现象，对人类生命、关键基础设施和山地生态系统构成严重威胁[1,2,3,4]。2004年至2016年间，全球共记录了4862起非地震性滑坡事件，导致约55,997人死亡[5]。在全球变化背景下，由于气候变化和人为干扰的增加，滑坡灾害可能会加剧，从而提高易发地区的滑坡频率[3,6,7]。因此，及时准确地识别现有滑坡和滑坡易发区域对于山区灾害评估、风险缓解和灾害预防至关重要。
遥感已成为滑坡调查的主要工具，能够在偏远和难以到达的山区进行重复的大规模观测。光学图像有助于识别地表沉积物、植被扰动和光谱异常，但在持续云层覆盖、地形阴影和明显的季节性地表变化条件下效果较差[4,8]。合成孔径雷达（SAR）能够在全天候、昼夜条件下进行地球观测，在光学数据难以获取的地区特别有价值[4]。然而，在山区，其有效使用受到阴影、几何畸变和密集植被影响下的可解释性降低的制约[9,10]。星载干涉合成孔径雷达（InSAR）能够直接测量地面变形，对于检测缓慢移动或逐渐失稳的斜坡非常有效[11,12,13,14]。但在时间相关性严重、植被覆盖密集或斜坡失稳速度过快而无法有效捕捉的情况下，其能力会减弱[10,15,16]。这些限制表明，没有单一传感器能够同时满足及时检测、大面积覆盖和连续时间监测的需求。在这种情况下，卫星嵌入数据集通过将多源地球观测数据整合到一个紧凑的特征空间中，提供了比单一传感器观测更有效的复杂地表扰动特征描述方法[17]。对于多云的山区，年度嵌入表示可以减少对无云光学观测的依赖，并在频繁云层覆盖的情况下提供更稳定的地表条件描述。然而，由于这些数据集只是最近才与大规模地理空间基础模型一起出现，其在滑坡检测中的应用仍有限。

**2. 材料与方法**
**2.1. 研究区域**
研究区域如图1所示。云南省位于中国西南部（图1a），研究区域位于该省的西部（图1b），包括十个县级行政单位：贡山德荣和怒族自治县、福贡县、泸水市、云龙县、龙阳区、石甸县、龙陵县、芒市、永德县和镇康县，总面积约为32,800平方公里。怒江自北向南流经该区域。该地区海拔超过3000米，形成了深切的高山峡谷地貌，坡度陡峭，山谷深邃（图1c）。这些地形条件加上频繁的地质灾害和脆弱的生态环境，使得该地区极易发生滑坡。由于自然和社会经济限制，怒江干流尚未开发大规模水电站，水电开发受到限制[29,30]。此外，该地区还受到广泛的人为干扰，包括小型水电站、道路网络和沿河谷及交通走廊的局部工程活动，这些活动可能改变坡面几何形状并降低局部坡面稳定性。在这种条件下，及时和全面的滑坡调查对于灾害缓解、社区韧性和未来大规模基础设施的可持续规划至关重要。尽管具有实际意义，但该地区的系统滑坡调查仍较为有限。

**2.2. 数据与预处理**
本研究使用的数据集如表1所示。Google卫星嵌入数据主要用于年度滑坡识别，而Shuttle Radar Topography Mission（SRTM）数字高程模型（DEM）、MERIT Hydro、WorldCover 2021数据、气候数据和人为指标则用于时空分析和滑坡易发性建模。基于这些数据集，推导出24个多维环境因素，并将其分为四类：地形-地质因素、气候-水文因素、土地覆盖因素和人类活动因素。此外，通过实地调查和HR光学卫星图像的视觉解释，建立了2022-2025年的滑坡样本，用于模型训练和准确性评估。

**2.2.1. 卫星嵌入数据集**
本研究使用了GEE中的年度Google卫星嵌入数据集。该数据集由大规模Google AlphaEarth基金会模型生成，提供10米空间分辨率的全球地理空间嵌入。每个像素由一个64维向量表示，总结了来自多个地球观测数据的年度地表条件和时间动态[17]。由于嵌入是基于多年观测的年度表示，因此比传统的基于图像的解释方法更少依赖于无云光学场景，特别适用于受季风影响的山区。与传统光谱带相比，这些学习到的嵌入特征能更有效地捕捉复杂地表变化，因此非常适合分类和变化检测任务[32]。使用2017-2025年的年度嵌入图像生成了年际地表变化强度层，作为视觉解释、样本选择和滑坡制图的辅助信息。

**2.2.2. 地形-地质因素**
地形和地质因素用于表征与滑坡发生相关的地形形态、地表粗糙度和岩石学背景条件。地形变量包括海拔（图1c）、坡度、朝向、总曲率、地形起伏和地形粗糙度（以海拔的局部标准差表示）。所有这些变量均来自30米SRTM DEM，提供1弧秒空间分辨率的近全球高程数据[33]。
还包括了一个重新分类的地质图（GeologyReclass），用于表示研究区域的岩石学特征。具体来说，将原始1:50,000比例尺的区域地质图扫描后数字化，并根据系列、系统和岩石学特征重新分类为12种地质类型（图1d）。为确保后续建模的空间一致性，所有地形和地质因素均重新采样为30米的空间分辨率。

**2.2.3. 气候-水文因素**
气候和水文因素用于表征与滑坡发生相关的外部触发条件和水文环境。气候变量包括年平均降水量（MeanAnnualPrecip，图1e）、平均雨季降水量（MeanWetSeasonPrecip）、年平均最大日降水量（MeanAnnualMaxDaily）和极端降水事件的平均频率（MeanExtremeFreq）。本研究将雨季定义为5月至10月，极端降水事件定义为降水量超过50毫米的天气。这些变量来自2000-2025年的CHIRPS日降水量数据集，用于表示研究区域的长期降水条件、季节性降雨集中程度、降雨强度和极端降雨活动[34,35]。
水文变量包括上游汇水面积（FlowAcc_UPA）、上游坡度相关的流量积累（FlowAcc_UPG）、距最近排水口的距离（HAND）[36]、河流能量指数（SPI）[37]、地形湿度指数（TWI）[38]、河流密度（RiverDensity）和距河流的距离（DistRiver）。这些变量主要来自MERIT Hydro数据集，用于描述径流积累、局部湿度条件和排水影响以及潜在的河流侵蚀效应[39]。其中，UPA和UPG表征地表径流的潜在集中程度；HAND反映了位置相对于最近排水通道的相对垂直位置；SPI表示集中水流的侵蚀能力；TWI表示水分积累和地表湿润的趋势。河流密度和距河流的距离进一步反映了排水网络发展和河流接近度对坡面不稳定性的影响。

**2.2.4. 土地覆盖因素**
土地覆盖因素用于表征地表覆盖状况、植被状况和地表暴露程度。本研究考虑的变量包括土地覆盖类型（WorldCover）、平均归一化植被指数（NDVI_mean）[40]、NDVI振幅（NDVI_amplitude）和平均裸露地表指数（BSI_mean）[41]。WorldCover来自ESA WorldCover 2021产品，用于表示研究区域的主要土地覆盖类型的空间分布[42]。NDVI_mean、NDVI_amplitude和BSI_mean来自2019-2025年的Sentinel-2地表反射率数据。

**2.2.5. 人类活动因素**
人类活动因素用于表征人为干扰对滑坡发生的影响。本研究考虑了两个变量：距离道路的距离（DistRoad）和距离居民点的距离（DistSettlement）。DistRoad数据来源于OpenStreetMap的道路数据，而DistSettlement数据来源于ESA WorldCover 2021产品的建成区类别。这两个变量都被转换为距离栅格，以表示每个像素与主要人类活动特征的接近程度。2.2.6. 野外调查和高分辨率（HR）影像样本和验证数据来源于野外调查和高分辨率光学卫星影像。野外调查从2022年持续到2025年，主要在研究区域的国家和省级道路上进行，特别关注了交通走廊附近的斜坡、河谷、居民点以及工程扰动区域。野外调查点记录在滑坡体的中心附近，而不是滑坡边缘，坐标精确到小数点后五位，大约相当于1米的精度。2022年、2023年、2024年和2025年分别调查了259个、240个、256个和300个野外观测点，其中在野外验证了45个、55个、62个和70个完整的滑坡。为了补充野外数据，从中国资源卫星数据与应用中心获取了2022年至2025年间的无云高分辨率卫星影像，覆盖了整个研究区域。这些影像包括高分一号/二号/六号（GF-1/2/6）和资源一号/三号（ZY-1/3）的数据，并被锐化到2米的空间分辨率。每年使用1月至6月期间获取的无云高分辨率影像，而野外调查则在7月和8月进行。当同一区域有多个高分辨率影像时，优先选择云量较低、地形可见度较高且拍摄日期接近野外调查的影像。野外调查数据和高分辨率影像是样本解释和准确性评估的主要参考。构建滑坡和非滑坡样本的详细程序在第2.3.2节中描述。本研究使用的高分辨率卫星影像的详细信息见表2。表2. 本研究使用的高分辨率卫星影像。2.3. 方法论本研究的方法论框架如图2所示。该框架围绕四个功能组件设计：嵌入变化分析用于解释指导、年度特定样本构建、年度滑坡制图、易发性评估和时空分析。首先，使用卫星嵌入数据集生成年际表面变化强度层，为后续的手动解释提供辅助空间指导。其次，通过整合嵌入变化强度图、高分辨率光学影像和野外调查数据，构建了2022年至2025年的年度滑坡和非滑坡样本集。然后，使用2022年至2024年收集的样本训练随机森林（RF）模型，并将其应用于生成2017年至2025年的年度滑坡图。最后，利用多年滑坡图来表征滑坡表面的空间分布，并将24个多维环境变量纳入滑坡易发性建模中，以生成易发性地图。图2. 本研究的技术流程图。2.3.1. 嵌入变化分析用于视觉解释指导在卫星嵌入特征空间中测量年际表面变化强度，为视觉解释提供辅助指导。滑坡通常会导致表面形态、物质暴露和植被状况的突然和局部扰动。因此，受滑坡活动影响的像素预计在嵌入特征上表现出比稳定区域更大的年际偏差[17,43,44]。据此，基于嵌入的变化分析用于突出显示地表扰动的异常区域，从而在复杂的山区地形中缩小手动解释的搜索范围。如图3a所示，2023年和2024年稳定区域的样本在UMAP投影的嵌入空间中大部分重叠[45]，表明它们的嵌入特征在相邻年份间相对稳定。相比之下，受滑坡影响区域的样本在2023年的事件前状态和2024年的事件后状态之间显示出明显的位移，表明滑坡发生会导致嵌入特征的显著变化。这种对比表明年度嵌入特征对滑坡相关的地表扰动敏感。图3. 2023年和2024年稳定区域与受滑坡影响区域之间的嵌入空间特征比较。(a) 2023年和2024年稳定区域的年度卫星嵌入特征以及滑坡发生前后的受滑坡影响区域的UMAP可视化；(b) 2023年和2024年嵌入向量之间的余弦距离的箱线图。在此基础上，对2017年至2025年相邻年份的年度嵌入产品进行了逐像素比较。具体来说，使用64维嵌入向量计算的余弦相似度比较了2017–2018年、2018–2019年……和2024–2025年的年度嵌入向量对。采用这种度量是因为嵌入被设计为一个综合特征表示，所有维度共同编码年度表面特征。对于给定的像素i，让 和 表示两个相邻年份的嵌入向量。余弦相似度按照公式(1)计算：(1) 相应的嵌入变化强度定义为ΔE：(2) ΔE的统计显著性在图3b中展示。与稳定区域相比，受滑坡影响区域在2023年和2024年之间的余弦距离显著增大，证实了通过余弦距离测量的年度嵌入变化可以有效捕捉滑坡相关的地表扰动。对所有相邻年份对重复这一过程后，得到了一系列年度嵌入变化图，这些图随后被用作识别可能受滑坡活动影响区域的辅助证据。连续的ΔE图没有通过阈值处理转换为二进制滑坡掩膜，而是作为辅助嵌入变化层来指导高分辨率影像的视觉检查。高ΔE值也可能反映非滑坡地表扰动；因此，滑坡归属依赖于高分辨率影像和野外观察的地形证据，而不仅仅是嵌入变化。相反，低ΔE值并不自动表示非滑坡状况，当高分辨率影像或野外记录中存在可识别的地形证据时，即使是微小的滑坡也会被解释出来。然而，那些在ΔE图和高分辨率影像中都缺乏可检测表面表现的低对比度滑坡可能在当前工作流程中仍然无法识别。2.3.2. 从卫星影像和野外调查中构建训练样本训练样本是根据野外调查和手动解释的高分辨率光学影像构建的，因为现有的滑坡清单数据集并不完全适合10米空间分辨率的年度滑坡制图。在某些情况下，现有的清单包含位置不准确的信息，没有及时更新，可能包括已经稳定或重新植被化的古老或早期滑坡[26,27]。因此，使用这些清单进行监督训练可能会引入大量的标签噪声并导致显著的制图错误。为了解决这些限制，专门使用高分辨率影像和野外数据构建了2022年至2025年的滑坡样本数据集（图4a）。在样本解释过程中，连续的嵌入变化强度层与高分辨率影像和野外调查记录一起使用，以指导视觉检查，而不是作为基于阈值的筛选掩膜。最终样本标签主要根据高分辨率影像和野外证据确定，ΔE值低的区域并没有被自动排除。滑坡样本主要是从显示可识别滑坡特征的位置的高分辨率影像和野外证据中手动识别的，包括弧形或不规则的疤痕边界、裸露的物质、下坡移动痕迹以及堆积或沉积模式。非滑坡样本是从周围地形单元中选择的，这些单元没有滑坡失败的视觉证据。与滑坡无关的表面变化，如道路建设、农业活动和裸露的非滑坡表面，根据地形、坡度位置和野外证据被排除。非滑坡样本是使用相应年份的高分辨率影像分别解释的，并且不会在年份之间重复使用作为固定标签。具有模糊地形证据的区域，包括疑似旧沉积物、重新植被化的疤痕、道路切割区和建筑工地，被排除在非滑坡样本集之外。最终样本集旨在代表当代滑坡状况，而不是来自过时或混合时间序列清单的信息。图4. 手动解释样本集的分布。(a) 高分辨率影像；(b) 基于黄色六边形网格的视觉解释；(c) 样本集的空间分布。采用基于六边形网格的采样策略来减少空间自相关并提高样本集的空间代表性[46,47]。整个研究区域被划分为10平方公里的六边形单元，在每个单元内进行随机采样。对于非滑坡类别，每个单元随机选择五个样本点；而对于滑坡类别，根据滑坡规模选择两个到五个样本点（图4b）。这种策略允许样本覆盖更广泛的地形、土地覆盖和扰动条件，同时减少在有限数量的滑坡热点内的过度聚集。总共获得了3368个滑坡样本和65,325个非滑坡样本，如图4c和表3所示。表3. 本研究收集的样本。2.3.3. 使用GEE和RF算法绘制滑坡地图年度滑坡制图是在GEE中实现的，GEE提供了一个集成环境，用于多源数据访问和基于云的计算，从而能够在区域尺度上进行多年分类[48,49]。在本研究中，GEE用于处理候选滑坡区域、管理解释的样本数据，并在统一的工作流程中训练RF分类器。绘制的类别代表每个年度影像记录中可通过遥感识别的滑坡表面，包括新暴露的疤痕、活动表面以及早期或重新激活的滑坡的仍然可见的表面。年度地图描绘了可检测的滑坡表面，并不构成严格的年度新事件清单。此外，绘制的类别没有进一步细分为特定类型的滑坡，因为类型特定的分类超出了这项区域年度制图研究的范围。2022年至2024年的年度样本集用于训练，而2025年的样本集用于时间独立但空间不重叠的验证。将训练好的分类器应用于2017年至2021年涉及在同一研究区域和年度嵌入特征空间内的时间转移。这种转移得到了可检测滑坡表面的共同影像特征的支持，包括暴露的物质、植被扰动和疤痕形态模式。在分类工作流程中，每个目标年的年度卫星嵌入影像作为预测数据集，解释的样本用于提取相应的嵌入特征以进行RF训练。RF分类器配置了200棵树，以平衡多年GEE制图中的分类稳定性和计算效率。树的数量敏感性测试表明，F1分数和OOB分数在大约170棵树后趋于稳定，支持使用200棵树作为稳定设置。其他参数保持默认值，以保持年份间模型结构的一致性。然后将训练好的分类器应用于年度嵌入影像，将研究区域内的像素分类为滑坡或非滑坡。分类器输出是一个二进制地图，其中滑坡像素被赋予值1，非滑坡像素被赋予值0。2.3.4. 使用多维因素进行滑坡易发性建模滑坡易发性建模使用了多年滑坡制图结果和第2.2节中描述的24个多维环境因素。为了构建易发性训练数据集，使用2017年至2025年年度滑坡地图中的滑坡斑块作为正样本的来源，而滑坡斑块外的区域作为负样本的来源。为了减少噪声，在提取样本之前移除了小于1000平方米的碎片化映射斑块。然后从保留的滑坡斑块的中心部分选择正样本，而不是从斑块边界选择，因为斑块边界处更可能出现混合像素和分类不确定性。总共选择了5000个滑坡样本和10,000个非滑坡样本，并将数据集随机分为70%的训练子集和30%的验证子集。采用RF分类器进行易发性建模，其超参数通过网格搜索结合五折交叉验证进行优化[24,50]。根据优化结果，RF模型配置了以下参数：树的数量=400，最大深度=20，最小样本分割=2，每个叶子的最小样本=1。检查了24个候选因素之间的相关性，以评估输入变量之间的潜在冗余。如图5a所示，几组变量表现出相对较强的成对相关性，尤其是在水文和降水相关变量之间，这表明直接包含所有因素可能会导致冗余。因此，应用了递归特征消除与交叉验证（RFECV）[51]方法，并使用得到的最优特征子集进行最终的易发性建模。根据RFECV的结果（图5b），移除了三个变量——FlowAcc_UPA、FlowAcc_UPG和SPI，而其余因素被保留用于模型构建。图5显示了用于滑坡易发性建模的24个环境因素的相关性结构和基于递归特征消除与交叉验证（RFECV）的特征选择结果。(a) 24个因素的相关性热图；(b) RFECV曲线显示了随着保留特征数量的变化模型性能的变化，其中最优特征子集对应于最高的曲线下面积（AUC）。红色虚线表示最优特征数量。

2.3.5. 性能评估指标
年度滑坡制图和滑坡易发性建模的性能使用基于混淆矩阵的指标进行评估。具体来说，使用生产者准确率（PA）、用户准确率（UA）、总体准确率（OA）、卡帕系数（kappa coefficient）和平均交并比（mIoU）来评估年度滑坡地图，而接收者操作特征（ROC）曲线和ROC曲线下面积（AUC）用于评估易发性模型的区分能力。公式如下：
(3) (4) (5) (6) (7) (8) (9) (10)
其中TP、TN、FP和FN分别表示真正例、真负例、假正例和假负例的数量，N是样本总数。

3. 结果
3.1. 准确性评估
3.1.1. 滑坡制图的准确性
首先使用2025年的验证样本集评估了基于嵌入的滑坡制图的准确性，该样本集在时间上是独立的，但在空间上没有重叠，如表4所示。共有16,300个非滑坡样本和831个滑坡样本被正确分类，其中5个为假正例，92个为假负例。滑坡类别的PA、UA、OA、F1分数和kappa系数分别为90.03%、99.40%、99.44%、94.49%和0.9419，表明基于点的分类性能很强。尽管样本不平衡，但高的滑坡PA、UA和F1分数表明在当前的验证设置下，基于嵌入的分类器并没有明显偏向于多数非滑坡类别。

表4. 2025年滑坡制图结果的混淆矩阵和准确性指标。
为了进一步检验卫星嵌入数据集的有效性，我们还使用原始的10个多光谱波段训练并评估了一个基于Sentinel-2的RF模型。相同的2025年独立现场样本集被用于验证。如表5所示，基于Sentinel-2的结果获得了96.49%的OA和75.60%的UA。然而，其滑坡PA和F1分数仅为51.03%和60.93%，表明有大量的滑坡样本被遗漏。

表5. 2025年基于Sentinel-2的滑坡制图结果的混淆矩阵和准确性指标。
尽管如此，基于点的验证主要反映了样本级别的分类准确性，并不能完全描述边界划分或斑块级别的空间一致性。为了解决这一限制，使用2017-2025年的手动解释的真实图像块进行了代表性的对象级别比较。每年选择了20个代表性块，以覆盖不同的县、地形、滑坡大小和图像条件。这些块不仅来自视觉上明显的滑坡；每个块都包括手动解释的滑坡表面和周围的非滑坡区域，以评估遗漏和误报错误。对于2022-2025年，参考解释得到了高分辨率（HR）图像和现场证据的支持。对于2017-2021年，由于缺乏HR图像，参考块主要基于Sentinel-2图像和地貌证据进行解释。因此，2017-2021年的基于Sentinel-2的指标并不完全独立，而是作为比较参考值报告，而不是严格的独立准确性估计。

表6显示，基于嵌入的结果在F1分数、kappa系数和mIoU方面始终优于基于Sentinel-2的结果。平均F1分数从Sentinel-2的50.44%增加到基于嵌入方法的73.60%，而平均mIoU从63.65%增加到77.53%。这些结果表明，卫星嵌入特征比原始的Sentinel-2多光谱波段提供了更有效的滑坡相关表面扰动的表示，特别是在对象级别的滑坡划分方面。2021年基于嵌入的结果的相对较低的UA和F1分数表明该年的误报错误增加，这可能与参考斑块的碎片化、植被恢复以及与其他暴露或受扰动表面的混淆有关。2021年的异常值也突显了年际间对象级别制图不确定性的变化。尽管2017-2021年的Sentinel-2指标的独立性有限，但这种比较提供了有用的证据，表明卫星嵌入特征改善了遥感可识别滑坡表面的划分，特别是在2022-2025年这些年份得到了HR图像的支持。这种基于块的评估改进了斑块级别空间一致性的评估，但仍是一种代表性的对象级别比较，而不是完全随机的验证设计。

3.1.2. 滑坡易发性模型的性能
图6显示了滑坡易发性模型测试集验证的ROC曲线和精确度-召回率曲线。共有1500个滑坡样本和3000个非滑坡样本用于评估模型的性能。模型获得了0.9697的AUC，表明在基于清单的易发性数据集中，滑坡样本和非滑坡样本之间的区分度很高。从精确度-召回率曲线得出的平均精确度为0.9426，表明模型在类别不平衡的情况下保持了高精确度。此外，随着召回率的增加，模型保持了高精确度，表明在识别滑坡样本方面表现稳健。这些指标反映了基于映射的清单易发性数据集中的区分能力，相关不确定性在4.4节中进行了讨论。

图7显示了研究区域内10米分辨率下映射的滑坡表面的空间分布。映射的滑坡表面显示出明显的空间聚集，而不是在整个流域内的均匀分布。根据经度和纬度分布曲线，滑坡区域主要集中在98.5°至99.0°E和25.6°至26.4°N之间。这些高密度区域主要集中在怒江大峡谷内，特别是在贡山、福贡和泸水。这些空间统计数据来自映射的滑坡表面，因此作为描述性基于映射的估计报告，而不是经过误差调整的面积估计。

图8显示了怒江中下游地区滑坡易感性的空间分布和县级组成。滑坡易感性使用固定的等间隔RF易感性分数阈值分为五个相对等级：非常低（0-0.2）、低（0.2-0.4）、中（0.4-0.6）、高（0.6-0.8）和非常高（0.8-1.0）。这些等级代表相对的易感性水平，而不是精确的绝对概率。采用等间隔方案是为了提供透明且可复制的分类规则，并避免由自然断裂或分位数分类产生的分布依赖性阈值。在区域尺度上，非常低易感性等级占据了研究区域的最大比例（87.8%），而低、中、高和非常高等级分别占5.0%、2.5%、2.2%和2.5%。这种分布表明RF易感性分数强烈偏向低值，研究区域中的大多数像素具有相对较低的模型易感性。广泛的非常低易感性区域主要位于人口稀少的山区。相比之下，虽然高和非常高的易感性等级所占面积相对较小，但它们主要分布在怒江的主干和支流沿线。这些河谷走廊与居民点、道路和其他人类活动更为紧密相关，意味着人口和基础设施面临更大的滑坡风险。

图9显示了研究区域内不同距离带与道路和居民点之间滑坡易感性的比例组成。在0-500米的范围内，高和非常高易感性等级的比例相对较低；在500-1000米范围内显著增加；然后随着距离的增加逐渐减少。这种模式表明滑坡易感性与人活动密切相关，同时也表明最高易感性并不直接出现在道路和居民点附近。相反，最易感的区域集中在中等距离带，而靠近道路和居民点的斜坡可能通过工程措施得到了更好的保护或稳定。

4. 讨论
4.1. 滑坡发生的主要控制因素
确定滑坡发生的主要控制因素在易感性分析中至关重要，因为它有助于确定模型的预测行为是否与潜在的环境过程一致。大量证据表明，滑坡发生受到地形、水文、地质、土地覆盖、气候和人为干扰的影响[52,53,54,55]。严格解释这些因素对于将统计或机器学习模型与基于过程的理解联系起来非常重要[23,25,56]。最近的研究应用了可解释性方法，如Shapley Additive Explanations（SHAP），来识别与滑坡发生相关的主要因素，因为这些方法可以量化预测变量的相对重要性及其对模型输出的影响方向[57,58]。

SHAP分析表明，研究区域的滑坡发生主要与表面条件和地形相关变量有关。如图10所示，BSI_mean和NDVI_mean具有最大的平均绝对SHAP值，其次是RoughnessStddev、Aspect和TerrainRelief，而RiverDensity和TotalCurvature等变量的贡献相对较小。较高的BSI_mean值增加了模型输出，而较高的NDVI_mean值减少了模型输出，表明与裸露表面暴露和植被损失有很强的关联。然而，BSI_mean和NDVI_mean可能同时反映了滑坡前的表面条件和滑坡后的情况，包括暴露的物质、植被扰动和随后的恢复；因此，它们的SHAP贡献表明了预测关联，而不仅仅是纯粹的因果控制。这些关联代表了综合的滑坡类别，而不是针对个别失效机制的特定类型控制。地形粗糙度、起伏度、坡度和靠近道路和河流的距离也增加了易感性，表明了地形切割和人为干扰的影响。坡向显示出相对较高的贡献度，且滑坡像素主要集中在朝南至东南方向的斜坡上，这可能反映了与坡向相关的土壤湿度、植被状况、风化强度和局部山谷形态的差异。相比之下，与降水相关的变量贡献较弱，这可能是因为它们是通过长期统计指标来表示的，而不是事件规模的触发降雨。CHIRPS的相对较低分辨率可能会进一步平滑深切割峡谷地形中的局部地形降雨信号，从而减弱了与降雨相关变量的明显贡献。图10. 基于SHAP的解释，显示了2017-2025年研究区域内滑坡发生的主要控制因素。(左图)：根据平均绝对SHAP值对前15个预测因子进行排名。(右图)：相同预测因子的SHAP蜜蜂群图，显示了每个因素对滑坡易发性的影响方向和程度。点的颜色代表特征值，蓝色表示低值，红色表示高值。SHAP值表明一个因素是使预测向更高的滑坡易发性驱动还是向更低的滑坡易发性驱动，灰色的垂直线表示SHAP值为零。本研究的结果与以往的文献大体一致，但也揭示了当前滑坡清单的一个独特特征。裸露表面指标的正面贡献和与植被相关的变量的负面影响与先前的研究结果一致[16,55]。同样，坡度、地形起伏以及与道路和河流的接近程度的重要性也与早期研究结果一致，这些因素表明陡峭的地形、道路切割和河流侵蚀是滑坡发生的常见控制因素[59,60]。然而，与许多传统的研究不同，在这些研究中降雨、岩性或坡度是主要预测因子，而本研究中BSI_mean和NDVI_mean的排名最高。这种差异可能是由于训练数据集的时空特异性，该数据集来源于2017-2025年的10米分辨率滑坡测绘。因此，清单主要由近期发生的和地貌活跃的滑坡组成，这提高了表面暴露和植被损失指标的诊断敏感性。

4.2 滑坡的空间分布特征
滑坡的发生表现出明显的空间选择性，而不是随机分布。先前的研究表明，其空间分布通常受多种环境因素的控制，包括地形、水文、植被、岩性和人为干扰[10,22,52]。因此，分析这些因素对滑坡的梯度响应对于识别滑坡更可能发生的环境条件非常重要。为了研究滑坡的空间模式，对SHAP分析确定的15个最具影响力的因素进行了梯度分析，如图11所示。滑坡主要集中在这几个变量的相对狭窄范围内。大多数滑坡像素发生在海拔1000-3000米、坡度15-45°以及地形起伏200-600米的区域。就坡向而言，朝南至东南方向的斜坡占所有滑坡像素的约60%。此外，滑坡像素在河流（0-800米）和道路（0-1000米）附近高度集中。大约86%的滑坡像素发生在草地和树木覆盖类别中，平均NDVI值为0.2-0.6，NDVI振幅值为0.2-0.4。滑坡像素还主要分布在年平均降水量为800-1600毫米且河流密度值在0.04到0.12之间的区域。地质和曲率类别也显示出不均匀的分布，表明滑坡与研究区域内的特定岩性和地貌环境有偏好关系。图11. SHAP确定的15个最重要预测因子对应的滑坡像素的梯度分布。直方图显示了每个因素不同区间或类别内的滑坡像素比例。在“坡向”图中，标签如N、NE和E表示坡度方向；在“WorldCover”图中，BSV表示裸露/稀疏植被，PW表示永久性水体，ML表示苔藓/地衣；在“地质类别”图中，缩写表示冲积-洪积松散沉积物（ADLS）、坡积-残余松散沉积物（CRLS）、软质碎屑岩（WCR）、中等硬度碎屑岩（MHCR）、碳酸盐岩（CR）、低至中等级别变质软岩（LMGMSR）、火山岩（VR）、中等硬度变质岩（MHMR）和侵入岩（IR）。这些结果表明，研究区域的滑坡主要发生在具有中等至陡峭坡度和中等强烈地表切割的中等海拔斜坡上。这种空间模式与河流侵蚀、坡脚侵蚀、道路切割和排水相关干扰的综合影响一致，这些都是山区斜坡稳定性的已知控制因素[61,62]。此外，不同地质类别之间的对比表明，在陡峭且高度切割的地形条件下，岩性背景对滑坡发生有重要控制作用。

4.3 绘制滑坡影响斜坡单元的时间动态
在斜坡单元尺度上的分析为复杂山区环境中绘制滑坡影响斜坡单元的持续性和时空变异性提供了重要见解。与基于像素的方法相比，基于斜坡单元的分析更符合地貌单元，因此提供了更具物理意义的滑坡过程的空间组织和时间演变的表示[63,64]。因此，应用r.slopeunits工具来分析怒江中下游大于1000平方米的滑坡的时空动态[65]。如图12所示，在研究区域内共划分了103,736个斜坡单元。该过程将像素级别的滑坡测绘结果汇总到对象级别，从而便于评估斜坡稳定性。滑坡影响斜坡单元的年度计数和比例直接从测绘结果计算得出，因此代表了基于测绘的描述性指标。图12. 滑坡影响斜坡单元的空间分布。左侧面板显示了2025年研究区域内的斜坡单元分布。右侧三个子面板（顶部、中部和底部）显示了斜坡单元边界与坡向图、高分辨率影像和滑坡测绘结果的叠加。右侧子面板中的彩色线条表示斜坡单元边界，底部子面板中的红色斑块表示绘制的滑坡表面。图13显示了2017-2025年间10个县滑坡影响斜坡单元的年际变化。从绝对数量来看，贡山、泸水、云龙和富贡县的滑坡影响斜坡单元数量始终最多，而芒市、镇康和永德县的滑坡影响斜坡单元数量相对较低。几个县在年份间显示出明显的滑坡峰值，表明时间变化非单调。泸水和富贡县在大多数年份保持了最高的滑坡比例，贡山也持续显示出较高的比例。相比之下，龙陵、石甸、芒市、镇康和永德县的滑坡比例在整个研究期间一直较低。云龙县显示出明显的年际变化，2017年和2019年的比例最高，随后几年显著下降。图13. 2017-2025年10个县滑坡影响斜坡单元的县级时间变化。条形图显示了10个县每年滑坡影响斜坡单元的数量，带圆圈标记的黑色线条显示了每个县滑坡影响斜坡单元与总斜坡单元的年度比例。在斜坡单元尺度上，绘制的滑坡既表现出空间持续性也表现出年际变化。在高活动性的县，易发生滑坡的区域的空间分布随时间相对稳定，表明这些区域主要由持久的地貌和环境因素控制。同时，绘制的滑坡影响斜坡单元数量的明显年际变化表明，尽管存在空间热点，但可检测的滑坡扰动仍然具有动态性。总体而言，绘制的滑坡影响斜坡单元主要集中在少数高起伏的峡谷县，尽管这些县的基于测绘的年度比例仍表现出明显的年际变化。

4.4 局限性和未来改进
尽管滑坡测绘和易发性评估结果的准确性很高，但仍存在一些局限性。年度样本主要基于高分辨率（HR）光学影像和现场证据进行解释。低对比度的滑坡，如移动缓慢、部分植被覆盖、重新激活或形态微妙的滑坡，在ΔE层和高分辨率影像中缺乏可识别的表面表现时可能会被遗漏。当前的易发性模型本质上是静态的。SRTM衍生的地形变量和WorldCover 2021土地覆盖层提供了地形和土地覆盖条件的静态表示，因此无法完全捕捉2017-2025年期间的局部地形变化、土地覆盖变化或植被恢复情况。这可能会在易发性建模中引入时间不一致性和不确定性。此外，研究区域的滑坡发生还受到动态和时变控制因素的影响，如降雨异常、短期地表扰动和渐进的斜坡变形，这些因素无法通过多年平均环境变量充分捕捉。不确定性也可能从滑坡测绘传播到易发性分区。移除小于1000平方米的斑块和基于中心的正样本提取减少了测绘噪声的影响，但并未完全消除与测绘相关的不确定性。由于这种不确定性量化需要基于概率的验证和误差调整的面积估计，因此没有为年度测绘面积、斜坡单元计数和热点持续性估计正式的置信区间。尽管InSAR提供了关于滑坡前兆和缓慢移动斜坡变形的宝贵信息，但由于本研究的主要目标是年度滑坡相关表面扰动的清单测绘，因此没有将其整合进来。在怒江中下游，陡峭的地形、密集的植被、时间去相关性和快速的斜坡变形给大规模InSAR应用带来了额外的挑战。未来的工作将受益于整合高分辨率光学影像、InSAR、现场调查、更新的DEM和年度土地覆盖信息，以构建更完整的滑坡清单，并更好地描述动态地形和表面变化。还需要一个动态的滑坡易发性框架来描述长期环境背景和短期触发条件[66,67]。此外，引入结合不确定性意识和可解释性的建模策略可以提高预测的可靠性并加强基于过程的解释[68,69]。

5. 结论
(1) 本研究开发了一种卫星嵌入辅助的工作流程，用于怒江中下游的年度滑坡测绘和易发性评估。通过整合卫星嵌入数据、高分辨率影像、现场证据、地理信息工程（GEE）和射频（RF）建模，生成了2017-2025年可遥感识别的滑坡表面的年度10米地图和基于清单的易发性地图。结果表明，卫星嵌入技术在复杂山区地形中的区域滑坡测绘具有潜力。
(2) 绘制的滑坡表面沿怒江走廊和相邻的高起伏峡谷斜坡显示出明显的空间聚集，特别是在贡山、富贡、泸水和云龙县。在斜坡单元尺度上，绘制的滑坡影响单元表现出明显的年际变化，而主要热点区域总体上保持稳定。
(3) 易发性分析确定BSI_mean、NDVI_mean、RoughnessStddev、Aspect和TerrainRelief是综合绘制滑坡类别的主要预测因子。这些结果表明，裸露表面暴露、植被状况、地形粗糙度、坡度方向和地形起伏与滑坡发生密切相关。然而，由于BSI_mean和NDVI_mean可能部分反映了滑坡后的光谱响应，它们的SHAP贡献表明了预测关联而非纯粹的因果控制。
(4) 结果是在基于遥感的年度滑坡测绘范围内得出的。年度地图描述了可检测的滑坡表面，而不是严格的新事件清单，而易发性地图代表了基于清单的背景易发性结果。因此，应用这些结果时需要考虑与滑坡后光谱响应、静态环境预测因子和基于清单的标签相关的不确定性。进一步整合变形监测、更新的地形和土地覆盖数据以及不确定性意识建模将改善动态滑坡监测和易发性评估。