观潮 · 模式预报与气象 AI 观察
15321处暴雨滑坡样本喂进XGBoost后 SHAP把主控因子指认为高程坡度NDVI AUC 0.828
> 2023 年 7 月 29—31 日台风「杜苏芮」外围水汽叠加副热带高压，北京西部山区 83 小时累计降水约 400–700 mm，门头沟全域、房山北部、延庆南部发生大量浅层降雨型滑坡。这篇把灾后实测的 15321 处滑坡栅格作为正样本、随机抽 15186 个非滑坡栅格做负样本，在约 2800 km²、30 m 分辨率、共约 311 万个栅格上用 XGBoost 做易发性建模，测试集 AUC 0.828；再用 SHAP 从全局贡献、非线性阈值、双因子交互、单样本误判溯源四个层面拆解模型。结论是大尺度地形主控——高程（平均绝对 SHAP 0.61）、坡度（0.43）、NDVI（0.39）排前三，极高与高易发区集中在门头沟西部、房山北部海拔 500–800 m、坡度 20°–35° 的破碎山地；作者同时明确交代了空间自相关泄漏、单次暴雨外推、SHAP 不能当因果这三条边界。
缺的是哪一块
滑坡易发性评价的老路是频率比、支持向量机、随机森林这类模型输出一张概率图，停在「变量重要性排序」——但降雨型浅层滑坡里坡度、高程、NDVI、降雨从来不是各自单独起作用，而是通过地形—水文—覆被耦合共同控制坡体稳定。只排重要性问题没解决：主控因子有没有稳定的非线性阈值？交互效应怎么刻画？模型解释在空间自相关和变量共线性条件下还能信多少？
传统统计模型在线性假设、变量独立性、因子交互刻画上都吃力；CNN 之类的深度学习模型能提空间特征，可解释性不足又卡住了它在灾害机理分析里的直接应用。这篇的路线是梯度提升加可解释机器学习：XGBoost 负责判别，SHAP 负责把每个预测拆成基准值加各因子贡献。
原文图 1：研究区概况。左上为北京市行政区划（研究区为门头沟全域、房山北部、延庆南部），左下为土壤类型分布，右为研究区高程与滑坡点叠加——红色高值集中在西北部中低山区。
模型怎么搭的
研究区 115°25′–116°15′E、39°30′–40°10′N，约 2800 km²，暖温带半湿润季风气候，年均降水 550–650 mm、七八月占七成左右，海拔 100–2303 m。九类影响因子统一采样到 30 m 栅格：
因子 · 类型 · 说明
高程、坡度、坡向 · 地形 · 由 DEM 提取，坡向 0°–360°
地形位指数 T · 地形复合 · T=LN[E/(Ē)·S/(S̄)]，由高程与坡度组合
NDVI · 植被 · 近红外与红光波段运算，取值 −1 到 1
土地利用类型 · 覆被 · 分类，one-hot 编码
暴雨期累计降水量 · 气象 · 2023-07-29 20 时至 07-31 13 时共 83 h 站点累计雨量插值到 30 m
土壤类型、土壤侵蚀强度 · 下垫面 · 土壤侵蚀按 RUSLE 分级，SE_11=微度、SE_12=轻度
样本与建模：与滑坡起始点重合的 15321 个栅格为正样本（来自 RLBJ v1.0 野外实测滑坡数据库），非滑坡栅格随机下采样 15186 个，合计 30507 个建模样本、正负近似均衡，按 8∶2 随机分层划分训练与测试。XGBoost 网格搜索定参 n_estimators=400、max_depth=4，靠叶节点数量惩罚、L2 正则化和采样策略压过拟合。
共线性先做体检：坡向、NDVI、暴雨雨量、土壤类型、土壤侵蚀、土地利用的 VIF 都低于 5；但高程、坡度、地形位指数的 VIF 分别为 15.500、14.290、54.705，远超阈值 10——原因很直接，地形位指数就是用高程和坡度算出来的，三者结构性相关。作者的处理是让 VIF 结果用于诊断结构性相关、并限定后续 SHAP 解释的适用范围，而不是把这三个因子踢掉。
原文图 2：九类影响因子空间分布。上排为 DEM、坡度、地形位指数，中排为坡向、NDVI、土地利用，下排为土壤类型、土壤侵蚀强度、暴雨期累计降水量——降雨高值区与西部中高山区高度重叠。
结果数字
判别性能与因子排序
测试集 AUC 0.828（随机分层划分、当前参数设置下）。平均绝对 SHAP 值排序：
排名 · 因子 · 平均绝对 SHAP 值
1 · 高程 · 0.61
2 · 坡度 · 0.43
3 · NDVI · 0.39
4 · 地形位指数 · 0.29
5 · 暴雨期累计降水量 · 0.25
6 · 土壤侵蚀强度（微度） · 0.21
7 · 坡向 · 0.20
8–9 · 土地利用类型、土壤类型 · 更低
原文图 3：平均绝对 SHAP 值条形图。Altitude（高程）+0.61 居首，Slope +0.43、NDVI +0.39、地形位指数 T +0.29、暴雨累计雨量 Rain_cum +0.25 依次排后；土地利用与土壤类型的各子类贡献都在 0.13 以下。
原文图 4：SHAP 蜂群图。横轴为对模型输出的影响方向，颜色为因子取值高低——高程、坡度上低值样本（蓝）大量落在负贡献区，高值样本（红）偏向正贡献；NDVI 的高值样本在正负两侧都有分布，形状明显比地形因子散。
非线性阈值：模型响应区间和实际样本集中区间不是一回事
高程的 SHAP 响应分段明显：0–300 m 以负贡献为主，超过 300 m 转为正向，500–1200 m 是明显的正向贡献区间。但实际滑坡样本和高易发单元集中在 500–800 m——作者专门区分了「模型正向贡献区间」与「实际样本集中区间」这两个概念，二者不是同一回事。坡度的优势区间是 20°–35°：更缓的坡下滑分力小，更陡的坡土层薄、基岩出露，浅层土质滑坡反而少。坡向 0°–100° 整体正贡献、150°–350° 多为负贡献，且响应随高程变化。
原文图 5：海拔的 SHAP 依赖图、按 NDVI 着色。300 m 以下大量深蓝样本压在负贡献区（最低到 −4），300–1200 m 抬升到 0 附近及以上，1300 m 以上重新回落；高 NDVI（红）样本集中在中高海拔段——这正是后面要说的共现问题。
原文图 6：坡向的 SHAP 依赖图、按高程着色。0°–100°（北偏东向坡）整体正贡献，150°–350° 逐步转负、300° 附近最负；红（高海拔）与蓝（低海拔）样本的响应走势不同，提示坡向与高程存在交互，作者同时提醒坡向还与日照、蒸发、植被和土体水分共变。
一个假阴性样本的完整溯源
作者从测试集里挑了一处「实际发生滑坡、但模型判为低概率」的样本做瀑布图拆解：基准值 E[f(x)]=0.004 是对数几率空间的平均输出，对应约 0.501 的概率——不是 0.004 的发生概率。该样本各因子贡献从大到小：
因子（取值） · SHAP 贡献 · 方向
坡度（6.227°） · −1.35 · 最强负贡献
地形位指数（0.746） · −0.59 · 负
耕地利用 · −0.48 · 负
高程（159 m） · −0.37 · 负
NDVI（0.566） · +0.30 · 正
坡向（86.037°） · +0.16 · 正
微度水力侵蚀 · −0.13 · 弱负
林地 / 轻度侵蚀 · −0.09 / −0.05 · 弱负
各贡献累加 f(x)=−2.594，对应预测概率约 0.070。一眼能看出问题在哪：这处滑坡发生在 6.2° 的缓坡、159 m 的低海拔——全都落在模型学到的「安全区间」里，是个典型的缓坡低海拔例外样本。
原文图 7：测试集中一处实际滑坡但被判为低概率样本的 SHAP 瀑布图。从基准 E[f(X)]=0.004 出发，坡度 −1.35 是最大下拉项，逐项累加到 f(x)=−2.594（概率约 0.070）。
原文图 8：滑坡与非滑坡样本对比。(a) 海拔箱线图，滑坡样本显著更高（P<0.001）；(b) 坡度箱线图，同样 P<0.001；(c) NDVI 核密度，滑坡样本偏 0.6–0.8；(d) 土地利用类型样本量柱状图；(e) 因子相关性热力图；(f) 样本空间分布，红点（滑坡）集中在西北部山地。
空间分布
对全部约 311 万个有效栅格逐一计算易发性、自然断点法分五级：极低、低、中、高、极高。极高与高易发区集中在门头沟西部、房山北部海拔 500–800 m、坡度 20°–35°、地形位指数较高的破碎山地带；低海拔平缓河谷和山前台地多为极低或低易发区。整体格局是「大尺度地形主控、局部多因子耦合调控」——中易发单元沿山间村镇和乡村道路呈线状或零散分布，反映土地利用扰动与局地微地貌的耦合。
作者承认的限制
- 空间自相关可能抬高测试精度。模型评价以随机分层划分为主，相邻栅格可能同时进训练集和测试集；XGBoost 的正则化并不能消除这种训练—测试泄漏，作者明确说需要补空间分块交叉验证。
- SHAP 不是因果。SHAP 反映的是模型关联结构，不能区分相关关系与因果关系；相关特征较强时贡献还会在特征间重新分配。典型例子是 NDVI：高 NDVI 的正向贡献主要反映它与中高海拔、残坡积物条件的空间共现，即空间混杂效应，不能据此认定植被覆盖直接促进滑坡。坡向同理，它与日照、蒸发、土体水分共变。
- 单次暴雨。降雨因子只对应 2023 年 7 月 29—31 日这一场暴雨（83 h 累计雨量），结论主要适用于该事件，向其他暴雨过程的外推未检验。
- 因子体系不完整。未包含断层距离、岩性分组、地下水和边坡工程等精细变量。
- 样本有漏测。滑坡样本依托灾后野外实测点位构建，小型隐蔽滑坡受勘察条件限制存在少量漏测，正负样本配比待无人机遥感普查补充后优化。
一句话判断
可解释机器学习在暴雨滑坡上的价值不在那 0.828 的 AUC，而在它能指出「高程 300 m 是贡献由负转正的门槛、20°–35° 是优势坡度带、那个被判错的样本错在缓坡低海拔」——但作者把边界画得很清楚：随机划分的精度有空间自相关水分，SHAP 排序是关联不是因果，单场暴雨的阈值换一场雨未必成立。