夏普利值揭示了地形因素对裸露基岩与砾石分异过程的控制作用

《Geoderma》：Shapley values reveal geomorphic controls on exposed bedrock-gravel differentiation

准确区分裸露基岩（EB）与砾石表面对高原地貌区土壤资源评估、侵蚀过程解析及生态保护至关重要。本研究提出融合地形数据、Sentinel-2光谱特征及可解释机器学习（RFE-RF）的分类框架，在西藏雅鲁藏布江流域4120km2区域验证。通过7,798个样本构建训练集，采用递归特征消除优化31个预测变量，随机森林模型在验证集和测试集分别达到94.28%和95.64%总体精度，显著优于传统单源模型（OA=85.78-88.91%）和TNLRS全国土地调查图（OA=88%）。Shapley值分析揭示坡高（SlpHeight＞25m）、地形位置指数（TPI）极性、短波红外波段（B11/B12）及黏土指数（ClaI）是核心驱动因子，形成"地形主导、光谱辅助"的解译机制，为高海拔地区精细地表分类提供可移植范式。

　　在高海拔地区，精确区分裸露岩面（EB）和砾石地表对于评估土壤资源、理解侵蚀过程对土壤形成的影响以及指导生态保护策略至关重要。传统的手动视觉判读方法不仅费时费力，而且成本较高，通常将两者归为一类“混合岩面-砾石地表”，这导致对土壤资源的误判。因此，建立一种新的分类方法，能够准确区分EB和砾石地表，成为该领域研究的重要课题。本文提出了一种融合地形特征、遥感光谱指数和可解释机器学习（ML）的方法，用于对中国西藏南部高山地区（平均海拔超过4,500米）的EB和砾石地表进行分类。该方法在高分辨率Google Earth Pro影像基础上生成了7,798个样本，并通过结合Sentinel-2光谱波段、土壤和植被相关指数以及数字高程模型（DEM）导出的地形变量，构建了一个递归特征消除-随机森林（RFE–RF）模型，取得了95.64%的整体准确率，显著优于传统方法（整体准确率为88%）。独立的实地验证进一步确认了预测结果的可靠性。通过Shapley分析，我们发现坡高和地形位置指数（TPI）是EB和砾石区分的主要驱动因素，反映了山脊地区的侵蚀过程和山谷中的沉积作用。短波红外波段（B11, B12）及其衍生指数（如黏土指数、地质指数）也进一步提高了区分效果。最终生成的分类图与高分影像和人工判读结果高度一致，表明该方法在高精度地表分类方面具有广泛的应用前景。

### 1. 引言

在高山生态系统中，EB和砾石地表的空间差异对理解地表过程具有重要意义。它直接影响土壤侵蚀动态、地质灾害风险以及生态系统恢复潜力。传统遥感（RS）分类方法主要依赖光谱特征或单一地形参数，这使得EB和砾石的区分变得困难，因为它们在光谱特性上存在相似性，例如Sentinel-2红边波段中的重叠反射。现有的土地资源调查，包括中国的第三次全国土地资源调查（TNLRS），将两者归为“EB和砾石地表”这一类别，这限制了对其生态功能和灾害易感性的准确评估。在这些调查中，“EB和砾石地表”被定义为地表土壤砾石覆盖超过70%的区域。然而，地表下仍可能存在未被利用的土壤。因此，改进能够区分EB和砾石的高精度空间数据集对于推动土壤特性研究和风化层厚度分析至关重要。尽管TNLRS地图提供了对EB和砾石地表空间分布的概述，但其在高山地区的可靠性尚未得到系统验证。然而，该地图仍是开发改进的EB和砾石制图方法的重要资源。

### 2. 材料与方法

#### 2.1. 研究区域

研究区域位于中国西藏自治区林芝市朗县，地理坐标为东经92°43′–93°60′，北纬28°67′–29°49′，地处喜马拉雅山脉北麓和雅鲁藏布江中下游。研究区域总面积约为4,120平方公里，主要土壤类型包括淡色黑钙土、冰冻淋溶土、普通淋溶土和钙质黑钙土。这些土壤主要由花岗岩、片麻岩、混合岩和泥质岩形成，质地较轻，易受水蚀影响。研究区域的地形特征包括高原山地和V形山谷，北侧和南侧较高，中部较低。景观主要由开阔的山谷、坡地和山地组成，海拔范围从2,918米到6,110米，平均约为4,500米。5,000米以上的区域广泛分布着裸露岩面、砾石和冰川。坡度范围从0°到85°，平均约为28°。该地区属于青藏高原的温暖、半湿润气候区，年均温为12.20°C，年均降水约为400毫米。

#### 2.2. 训练、验证和测试数据集生成及实地调查

为了评估传统TNLRS地图的准确性，我们从Google Earth Pro高分辨率影像中生成了EB、砾石和OLT样本。样本生成结合了视觉判读和分层空间采样，遵循标准化方法：（1）标准化分类定义：根据土地利用现状分类（GB/T 21010–2017）和岩石-地貌规格，定义了三种土地覆盖类别：EB（无植被覆盖的连续岩面）、砾石（地表主要由砾石覆盖，≥70%的面积）和OLT（植被、水体和人工地表）。（2）空间分层策略：根据海拔梯度将研究区域划分为五个子区域，并在每个子区域中进行平衡的随机采样，以确保空间覆盖和地貌多样性的均匀性。（3）质量控制机制：所有样本均与2米分辨率的高分影像进行了比对，最终获得了7,798个验证点：2,733个EB、2,499个砾石和2,566个OLT。

在2025年5月，我们进行了独立的实地调查，以验证EB和砾石分类的准确性。在代表性地点，使用GPS定位和摄影记录了38个EB和33个砾石样本。这些样本为评估EB和砾石分类图的准确性提供了独立的验证数据集。

#### 2.3. 环境协变量准备

##### 2.3.1. 地形导出变量

使用12.5米的Advanced Land Observation Satellite（ALOS）DEM，并通过立方样条插值将其重采样到10米分辨率，以匹配Sentinel-2数据集。通过Automated Geoscientific Analyses（SAGA）GIS v9.7.2，我们导出了22个地形预测变量。尽管尺度效应会影响地形导出变量，但所有协变量均使用SAGA的默认算法计算，确保了一致的10米分辨率作为基本映射单元。每个导出变量的有效空间尺度根据其计算原理有所不同（见表1）。局部尺度导出变量（如坡度、坡向、TPI）在原生10米分辨率下计算。多尺度指数（如多分辨率谷底平坦度指数和山脊平坦度指数）则通过内部平滑程序进行调整，以适应地貌特征。基于流的指标（如地形湿润指数、流功率指数）则通过水文路由确定其有效尺度。表1列出了本研究中使用的所有地形变量。

##### 2.3.2. Sentinel-2影像

我们从Google Earth Engine获取了Sentinel-2地表反射率数据（COPERNICUS/S2_SR_HARMONIZED），这是经过大气校正的Level-2A数据集，无需额外预处理。我们共融合了285个云覆盖小于10%的场景，这些场景的时间跨度从2018年1月至2024年12月。保留了10个高分辨率波段（B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12），并统一重采样到10米分辨率。从这些波段中计算了21个与土壤和植被反射率相关的光谱指数（见表3）。

#### 2.4. 方法

本研究分为四个主要步骤：（i）从Google Earth Pro中生成EB、砾石和OLT样本，以评估传统TNLRS地图的准确性；（ii）收集和预处理地形和Sentinel-2数据；（iii）实施四种特征工程策略，使用RF分类器生成高分辨率的EB和砾石地图——MODEL A：仅地形变量；MODEL B：仅Sentinel-2波段和导出指数；MODEL C：所有变量；MODEL D：通过RFE–RF优化后的所有变量；（iv）评估模型的准确性，并应用Shapley值以识别EB和砾石分类结果的主要驱动因素。

##### 2.4.1. RFE

RFE是一种特征选择算法，通过排名预测变量并逐步剔除最不重要的变量，直到找到最优子集（Guyon et al., 2002）。本研究采用RFE–RF方法，具体步骤如下：（1）初始化候选特征集（V1, V2, …, Vm）；（2）在当前特征子集上训练RF分类模型，并使用10折交叉验证（CV）重复10次，以Kappa系数和准确率作为评估指标；（3）根据重要性排名剔除最不重要的特征；（4）重复步骤（2）–（3），每次迭代剔除一个特征，直到只剩下一个特征。在每次迭代中，记录CV的Kappa系数和准确率，并选择最大化这两个指标的子集作为最优。如图4所示，具有约10个特征的模型表现出稳定的性能（Kappa > 0.88）。然而，具有最优31个特征子集的模型在CV中表现出最佳性能（Kappa = 0.92；准确率 = 94.30%）。因此，选择31个特征子集用于MODEL D，以最大化对EB和砾石的分类准确率。

##### 2.4.2. RF分类器

RF分类器被用于区分EB和砾石，因为它能够处理高维的RS和地形数据，对过拟合具有较强的鲁棒性，并且在大规模空间预测中具有计算效率（Breiman, 2001）。为了优化模型性能，我们系统调整了三个关键参数：树的数量（ntree）、每个分裂使用的变量数（mtry）和最小终端节点数（nodesize）。通过观察Out-of-Bag误差的收敛趋势，确定了最优的ntree值（0–1000）。通过10折CV评估了mtry（1–10）对分类准确率的影响，并通过网格搜索评估了nodesize（1–20）对模型复杂度的影响。参数选择结果见图S1。最终的参数设置为：MODEL A：ntree = 194，mtry = 8，nodesize = 3；MODEL B：ntree = 665，mtry = 6，nodesize = 5；MODEL C：ntree = 183，mtry = 3，nodesize = 2；MODEL D：ntree = 514，mtry = 5，nodesize = 2。

基于这些参数，我们使用RF模型进行了空间外推。分类概率用于量化分类的不确定性（Malone et al., 2017）。所有实验均在R 4.3.2中进行，模型训练使用了randomForest包（Liaw and Wiener, 2002）。CV通过caret包（Kuhn, 2008）结合空间分层采样，以确保参数优化结果的地貌代表性。

##### 2.4.3. 模型评估

对于每个模型，我们将7,798个点的完整数据集随机分为训练、验证和测试子集，比例为7:2:1。验证和测试集主要用于评估模型性能和生成分类图的准确性。关键评估指标包括整体准确率（OA）、Kappa系数、精确率、召回率和F1得分。这些指标的计算公式如下（Taghizadeh-Mehrjardi et al., 2020）：

- OA = (TP + TN) / (TP + TN + FP + FN)
- Kappa = 1 - (1 - OA) / (1 - Pe)
- 精确率 = TP / (TP + FP)
- 召回率 = TP / (TP + FN)
- F1得分 = (精确率 × 召回率) / (精确率 + 召回率)

其中，N表示样本数量；TP：真阳性；FP：假阳性；TN：真阴性；FN：假阴性。较高的OA、Kappa、精确率、召回率和F1得分表示更好的模型性能。Kappa系数的解释如下：0–0.20（差），0.21–0.40（一般），0.41–0.60（中等），0.61–0.80（显著），0.81–1.00（近乎完美）（Landis and Koch, 1977）。这些指标均来自混淆矩阵。

此外，模型性能还通过接收者操作特征（ROC）曲线和曲线下面积（AUC）进行评估。ROC曲线表示在不同分类阈值下真阳性率和假阳性率的关系。AUC量化模型性能：0.90–1.00表示优秀预测能力；0.80–0.90表示非常好；0.70–0.80表示好；0.60–0.70表示一般；0.50–0.60表示差（Kantardzic, 2011）。

最后，使用实地调查的点对最优的EB和砾石分类图进行了独立验证。

##### 2.4.4. 模型可解释性

为了识别影响EB和砾石空间分布的主要环境驱动因素，本研究采用基于Shapley值的局部可解释框架。Shapley值源于合作博弈理论（Kuhn and Tucker, 1953），用于量化每个预测变量对模型输出的边际贡献，从而为原本不透明的模型提供可解释性（Wadoux and Molnar, 2022）。Shapley值的计算公式如下（Lundberg and Lee, 2017）：

- ?_i(f) = ∑_{S?N/i} [S! × (N - S - 1)! / N!] × [f(S ∪ i) - f(S)]

其中，N是所有特征的完整集合；S是任何不包含特征i的子集；f(S)表示使用特征S进行模型预测；f(S ∪ i)表示特征i被包含时的模型预测；|S|表示集合S的大小；|N|是特征的总数。

Shapley值大于0表示对预测的正向贡献；小于0表示负向影响。通过计算所有样本的Shapley值的绝对值平均，得到了平均绝对Shapley（MAS）值。较高的MAS值表示特征的整体重要性（Wadoux et al., 2023）。全局特征重要性通过Shapley总结图可视化，而Shapley依赖图则揭示了个体特征如何影响模型输出，并帮助识别EB和砾石分类之间的环境阈值。此外，通过成对的Shapley相互作用分析，我们研究了关键变量（如坡高和TPI）对分类结果的综合影响。

所有基于Shapley的分析均在Python 3.12中使用shap模块进行。总结图、依赖图和相互作用值分析用于量化多源预测变量的全局和局部效应，从而为EB和砾石分布的环境控制提供机制性见解。

### 3. 结果

#### 3.1. 模型在验证和测试集上的性能

混淆矩阵提供了模型性能的详细见解（见图5）。在验证和测试数据集中，所有四个模型的分类结果都集中在对角线上，表明大多数样本被正确分类。MODEL C和MODEL D在对角线上的值更高，反映了它们在分类准确性方面优于MODEL A和MODEL B。此外，对于EB和砾石类别，MODEL C和MODEL D实现了更平衡的预测，减少了误分类。

#### 3.2. EB和砾石的空间分布模式

使用MODEL D生成的10米分辨率EB和砾石分布图（见图7）与传统TNLRS地图一致，但同时也识别出未被TNLRS地图分类的EB和砾石区域，这些区域已被高分辨率影像验证。在图7b中，可以观察到MODEL D更准确地代表了EB和砾石的空间特征，特别是在高山和陡坡区域，其输出与高分辨率卫星观测结果一致。EB主要分布在高海拔的山脊和陡坡区域，而砾石则集中在相对较低的地形区域。这种分布模式与高山环境中EB和砾石的自然分布相符，为模型的可靠性提供了支持。此外，预测地图显示出高可靠性和低不确定性，如图7e–f所示。

独立的实地验证进一步确认了这些预测的可靠性：38个EB样本中有37个（97.37%）和33个砾石样本中有32个（96.97%）被正确分类。

#### 3.3. Shapley值揭示EB和砾石差异的驱动因素

Shapley值分析提供了一个定量且可解释的框架，用于解释优化后的RF模型（MODEL D）所捕捉的关系。该分析不仅对变量的重要性进行了排序，还揭示了预测变量对EB和砾石分类概率的方向和大小影响。此外，它还识别了关键的环境阈值和地形与光谱驱动因素之间的协同作用。

Shapley值的分布（见图8、9、10、11、12）显示，地形变量如坡高和TPI在分类中的贡献显著。特别是，在研究中发现的25米坡高阈值（见图11、12）反映了侵蚀与沉积之间的动态平衡。当坡高超过25米时，陡坡侵蚀占主导地位，显著增加了EB出现的概率（正Shapley值）。这种模式对应于高能侵蚀条件下岩石风化阻力的降低和坡地不稳定性的增加。相反，当坡高低于25米时，沉积物运输和堆积成为主导因素，增加了砾石分类的概率（正Shapley值）。

同样，TPI的正负划分（见图11、12、13）展示了地形位置如何调节物质运输。负TPI区域（如山谷底部或凹陷地带）由于沉积物堆积而促进了砾石的沉积（正Shapley值），而正TPI区域（如山脊和高地形区域）则在持续侵蚀下有利于EB的形成。

光谱响应在短波红外波段（B11 = 1610 nm；B12 = 2190 nm）（见图9、10、11）突显了分类的矿物学基础。B11波段的高反射率与砾石正相关，这是由于次生黏土矿物在沉积区的吸湿性。黏土矿物的羟基吸收特征在B11/B12波段（黏土指数）产生了独特的光谱信号，而岩层由于其硅酸盐晶格的稳定性，表现出较低的B12反射率。相比之下，红边波段（B5–B7）的光谱模糊性（见图14）源于风化EB表面的氧化铁涂层，这些涂层的反射率模式与砾石沉积物相似。

Shapley分析还揭示了环境变量的多尺度协同效应（见图12）。在宏观尺度上，海拔和坡度定义了地貌框架：EB主要分布在高海拔和陡坡区域，而砾石则集中在较低的、坡度较缓的区域。在微观尺度上，导出参数如坡高和TPI捕捉了更精细的空间对比。例如，EB在5,000米以上的区域受到低温的限制，而砾石的形成在5,000米左右达到峰值。此外，TPI < 0且坡高 < 25米的区域砾石形成更为频繁。这种嵌套的“多尺度协同”有效地描述了地貌表面的异质性。

总体而言，这些结果表明Shapley值能够独特地量化EB–砾石差异的阈值、驱动因素和相互作用。特别是，25米坡高阈值区分了侵蚀和沉积主导的区域，短波红外波段（B11/B12）及其导出指数（如黏土指数、地质指数）解决了光谱模糊性，而TPI与坡高、海拔与坡度之间的协同效应则进一步明确了这些因素对分类结果的影响。尽管Shapley值已广泛应用，但其在多源RS和DEM数据与ML框架结合中的应用为高海拔景观的地表过程控制提供了新的、空间明确的理解。这种方法推动了可解释ML在地貌研究中的应用，产生了传统统计或实地方法难以获得的机制性见解。

### 4. 讨论

#### 4.1. 不同预测模型的性能比较

本研究开发了四种预测模型（MODEL A–D），以探讨多源数据融合和特征优化在区分EB和砾石地表中的机制。模型性能的显著差异不仅反映了输入数据集的物理特性差异，还反映了研究区域的空间异质性以及ML算法捕捉特征相互作用的能力。通过RFE–RF特征选择，MODEL D取得了最高的准确率（OA = 95.64%），超过了单源模型（OA = 85.78–88.91%）和传统TNLRS地图（OA = 88%）。独立的实地验证进一步确认了EB和砾石分类图的可靠性。

#### 4.2. Shapley值揭示的模型可视化和地理驱动机制

Shapley值的分布（见图8、9、10、11、12）揭示了地形变量如坡高和TPI在分类中的显著空间差异，反映了关键的地貌响应模式。特别是在本研究中发现的25米坡高阈值（见图11、12），代表了侵蚀和沉积主导的区域之间的动态平衡。当坡高超过25米时，陡坡侵蚀占主导地位，显著增加了EB出现的概率（正Shapley值）。这种模式对应于高能侵蚀条件下岩石风化阻力的降低和坡地不稳定性的增加。相反，当坡高低于25米时，沉积物运输和堆积成为主导因素，增加了砾石分类的概率（正Shapley值）。

同样，TPI的正负划分（见图11、12、13）展示了地形位置如何调节物质运输。负TPI区域（如山谷底部或凹陷地带）由于沉积物堆积而促进了砾石的沉积（正Shapley值），而正TPI区域（如山脊和高地形区域）则在持续侵蚀下有利于EB的形成。

光谱响应在短波红外波段（B11 = 1610 nm；B12 = 2190 nm）（见图9、10、11）突显了分类的矿物学基础。B11波段的高反射率与砾石正相关，这是由于次生黏土矿物在沉积区的吸湿性。黏土矿物的羟基吸收特征在B11/B12波段（黏土指数）产生了独特的光谱信号，而岩层由于其硅酸盐晶格的稳定性，表现出较低的B12反射率。相比之下，红边波段（B5–B7）的光谱模糊性（见图14）源于风化EB表面的氧化铁涂层，这些涂层的反射率模式与砾石沉积物相似。

Shapley分析还揭示了环境变量的多尺度协同效应（见图12）。在宏观尺度上，海拔和坡度定义了地貌框架：EB主要分布在高海拔和陡坡区域，而砾石则集中在较低的、坡度较缓的区域。在微观尺度上，导出参数如坡高和TPI捕捉了更精细的空间对比。例如，EB在5,000米以上的区域受到低温的限制，而砾石的形成在5,000米左右达到峰值。此外，TPI < 0且坡高 < 25米的区域砾石形成更为频繁。这种嵌套的“多尺度协同”有效地描述了地貌表面的异质性。

总体而言，这些结果表明Shapley值能够独特地量化EB–砾石差异的阈值、驱动因素和相互作用。特别是，25米坡高阈值区分了侵蚀和沉积主导的区域，短波红外波段（B11/B12）及其导出指数（如黏土指数、地质指数）解决了光谱模糊性，而TPI与坡高、海拔与坡度之间的协同效应则进一步明确了这些因素对分类结果的影响。尽管Shapley值已被广泛使用，但其在多源RS和DEM数据与ML框架结合中的应用为高海拔景观的地表过程控制提供了新的、空间明确的理解。这种方法推动了可解释ML在地貌研究中的应用，产生了传统统计或实地方法难以获得的机制性见解。

#### 4.3. 在其他山地地区的适用性与局限性

很少有研究将高分辨率地形、Sentinel-2影像和可解释ML结合，用于同时分类EB和砾石地表。大多数先前的研究仅关注特定环境背景下的EB，并报告了中等准确率（OA = 79–88%）。相比之下，我们的框架在青藏高原（TP）应用，该地区以稀疏的高山植被和强烈的地形对比为特征，通过融合光谱指数（如NDVI、SAVI、红边波段）和地形导出指标（如TPI、坡高、TRI），取得了95.64%的整体准确率。Sentinel-2的短波红外波段（B11/B12）在区分硅酸盐岩层和富含黏土的砾石方面表现出色，而TP的低植被覆盖率和活跃的侵蚀和沉积过程进一步增强了区分效果。

然而，在植被茂密或经常有云的山地地区，植被指数（如NDVI、SAVI）可能主导预测变量的重要性，掩盖了由SWIR波段衍生的矿物学信号。同样，在以碳酸盐为主的或高度风化的热带地形中，岩层和砾石之间的光谱差异可能减少，从而降低分类准确率。尽管如此，该框架的物理驱动因素具有广泛的应用前景：植被指数能够捕捉冠层覆盖；SWIR波段反映了不同岩石类型的矿物学变化；地形指标则代表了控制岩层暴露和沉积物积累的基本地貌过程。适应性的RFE–RF算法能够自动识别新环境下的最相关预测变量，而Shapley分析则提供了对特征贡献的透明、定量解释。尽管本研究仅在TP进行了验证，但物理可解释的预测变量集和数据驱动的设计表明该框架在多样化的山地环境中具有较强的可转移性。未来应在欧洲阿尔卑斯山、北美洛基山脉和南美安第斯山脉等地区进行有针对性的验证活动，以确定其性能边界并优化模型配置，以适应不同的高山生态系统。

为了克服光学RS在植被覆盖或云层多的区域的光谱限制，未来的研究应结合合成孔径雷达（SAR）数据（如Sentinel-1），这可以提供全天候成像和部分冠层穿透能力。SAR后向散射对地表粗糙度和土壤湿度敏感，使其能够通过物理结构（如密集岩层与松散砾石）区分EB和砾石（Bousbih et al., 2019）。极化信息（VV/VH）可以进一步表征微地形和湿度梯度（Nativel et al., 2022；Ullmann and Stauch, 2020）。多时相的Sentinel-1和Sentinel-2数据的协同分析也可以捕捉季节性风化和沉积动态，从而提高时间和空间分辨率（Vuolo et al., 2018）。

本研究旨在展示一个可解释、高精度的分类框架，而不是基准算法。RF模型在RS和地球科学领域已被广泛验证，具有较强的鲁棒性和计算效率，能够实现超过95%的高准确率。然而，我们并未测试更先进的算法，如XGBoost、LightGBM或深度学习模型，这些模型在某些条件下可能优于RF。未来的工作应探索系统比较和混合架构，如将RF与三维卷积神经网络或Transformer结合，以捕捉复杂的空间、光谱和地形相互作用（Liu et al., 2020）。

尽管RF模型表现出色（OA > 95%），但仍需承认一些局限性。首先，训练样本来源于Google Earth Pro影像的视觉判读，即使经过质量控制，仍可能引入主观偏差。其次，数据集是随机划分的，未考虑空间自相关性。第三，Sentinel-2和DEM数据的10米分辨率无法完全表示细尺度异质性；单个像素可能包含EB和砾石，导致混合像素效应。第四，地形导出变量（如TPI、坡度）是使用固定DEM分辨率和默认算法计算的，未优化邻域大小（如3×3与15×15像素），这可能无法捕捉最相关的地貌背景。最后，SHAP衍生的阈值需要与实地观测和地貌测量进行验证。未来的研究应通过以下方式解决这些问题：（1）应用空间CV方法（如k最近邻距离匹配，KNNDM）以严格测试泛化能力；（2）纳入更高分辨率或多源数据集以减少混合像素效应；（3）进行多尺度敏感性分析以确定地形协变量的最佳邻域大小；（4）扩大实地调查以收集更具代表性的样本，并评估SHAP衍生的阈值是否对应实际的地貌过程。这些改进将增强所提出框架的稳健性、严谨性和在多样化山地环境中的可转移性。

### 5. 结论

本研究开发了一个整合地形导出变量、Sentinel-2影像和可解释机器学习（ML）的框架，用于区分高山地区的EB和砾石地表。该研究的主要贡献如下：

1. **提高分类精度**：通过融合多源数据（地形和光谱指数）并使用RFE–RF进行特征优化，显著提升了预测性能。MODEL D取得了最高的准确率（OA = 95.64%），超过了单源模型（OA = 85.78–88.91%）和传统TNLRS地图（OA = 88%）。独立的实地验证进一步确认了EB和砾石分类图的可靠性。

2. **机制性见解**：可解释的Shapley分析定量识别了EB–砾石差异的主要驱动因素。地形阈值（如坡高 > 25米、TPI极性）和光谱特征（SWIR波段B11/B12、黏土指数、地质指数）成为关键预测变量。Shapley依赖图揭示了关键的环境阈值，如坡高 > 25米的EB预测和TPI < 0的砾石预测，以及预测变量之间的协同效应（如TPI与坡高、海拔与坡度之间的协同效应）。这些发现为高山地区地表差异的机制性、可量化的理解提供了支持。

3. **实际应用价值**：该框架生成了更精细的空间数据集，有助于地质灾害评估、生态系统恢复规划和土壤侵蚀建模。通过解决传统地图（如TNLRS）将EB和砾石归为一类的局限性，它提供了更准确的地表特征。由于该方法依赖于全球可访问的预测变量，包括Sentinel-2影像和地形导出指标，它在多种山地地区具有应用潜力，前提是进行本地校准和验证。