复现Nature Geoscience的VPD+SM AutoML+XAI归因解耦

干旱对植被的胁迫，其实有两条路：一条在脚下——土壤水分（SM）不够，根系吸不上水；另一条在头顶——大气饱和水汽压差（VPD）太高，空气"抽水"太狠，叶片只能关气孔。两者在干旱时往往同时出现，高度共线，传统相关分析很难把它们分开。

最近 Nature Geoscience 的一篇文章给出了一个很有冲击力的结论：

🌈
Li, L., Fang, Y., Mao, J., Xu, S., Zheng, Z., Li, W., Shi, M., Koven, C. D., Chambers, J., Leung, L. R. & McDowell, N. Vegetation increasingly controlled by atmospheric dryness rather than soil moisture.Nature Geoscience (2026). https://doi.org/10.1038/s41561-026-02107-w^[1]

多数地区植被仍由 SM 主导，但过去四十年 VPD 的影响在增强、SM 在减弱，而且"SM 主导→VPD 控制"的转移明显多于反向。

这篇文章的方法核心是"可解释人工智能（XAI）+ 自动机器学习（AutoML）"。我们从原始数据开始，把整套流程独立重跑了一遍，代码和全部图件已开源：

GitHub：https://github.com/GISWLH/vpd-sm-automl-xai

一、整体技术路线

一句话概括：每个 0.5° 格点单独训练一个模型，让模型告诉我们 NDVI 的波动有多少来自 SM、多少来自 VPD。

气象驱动：ERA5 逐日 0.25° 数据，聚合到逐月 0.5°，得到 VPD、SM、气温（Tas）、短波辐射（SWR），时间范围 1982–2022（01_prep_era5.py）。
植被信号：PKU GIMMS-MODIS 拼接 NDVI，半月合成→逐月 0.5°，按 QC 剔除云和雪（03_prep_anomalies.py）。
植被类型：MCD12C1 2001 年 IGBP 分类，0.5° 取众数并合并为原文的 9 类，比如北方灌丛、热带常绿阔叶林等（02_prep_landcover.py）。
去季节、去趋势：对所有序列做 STL 分解，取残差（距平），并且只保留生长季月份。这一步很关键：我们要解释的是"年际波动"，不是季节循环。

二、第一级归因：AutoML + TreeSHAP 拆解"谁在驱动"

每个格点上：

样本：中心格点加 3×3 邻域的生长季月份距平，自变量是 SM、VPD、Tas、SWR，因变量是 NDVI 距平；
AutoML：用 FLAML 在 Random Forest 和 XGBoost 之间自动选模型、调超参，每个格点 8 秒预算，80/20 划分，用 R² 评估；
TreeSHAP：对中心格点逐月计算 SHAP 值 φ，取 |φ| 的均值归一化成相对重要性 RI(%)；
时段对比：分 1982–2001 和 2002–2022 两段，分别算 RI，得到重要性的相对变化 δRI；
主导指数：D = RI_VPD − RI_SM，D<0 为 SM 主导；ΔD>0 表示向 VPD 控制转移。

工程上，3 个计算节点各开 48 进程（共 144 个 rank）并行（04_worker.py、run_node.sh），每个 rank 支持断点续算。最后筛掉 R² < 0.2 的格点，剩下 33,988 个有效格点。

三、复现结果①：SM 和 VPD 的相对重要性

面积加权平均 RI_SM = 26.5%，RI_VPD = 27.2%，两者旗鼓相当；
图 e 中蓝色为 SM 主导，占 42.3% 面积，主要在半干旱草地、稀树草原和澳大利亚、非洲南部；
图 g 按植被类型看：非北方灌丛（SHU_NB）、草地（GRA）偏 SM 主导，常绿阔叶林（EBF）、农田（CRO）偏 VPD；
图 h–l 是第二级 SHAP：背景 SM 均值解释了 D 空间格局的 43.5%，越湿的地方越偏 VPD 控制。

论文中原图：

四、复现结果②：重要性怎么变

两个时段对比：71.8% 的面积 VPD 重要性上升，45.3% 的面积 SM 重要性下降。VPD 在"变强"这个结论，在复现中非常稳健。

论文原图：

五、复现结果③：主导性转移（ΔD）

把"原来谁主导"和"往哪边变"交叉成四象限（图 a、b、c）：

象限	面积占比
SM 主导 → 更偏 VPD	26.9%
SM 主导 → 更偏 SM	15.4%
VPD 主导 → 更偏 VPD	36.0%
VPD 主导 → 更偏 SM	21.7%

总体 62.9% 的面积 ΔD>0。SM 主导区向 VPD 转移的面积（26.9%）明显大于留在 SM 一侧的（15.4%），这和原文说的"不对称转移"一致。

论文原图：

六、第二级归因：为什么会转向？

"解耦"的第二步，是把第一级得到的 D、δRI_SM、δRI_VPD、ΔD 当作因变量，用各格点的背景气候特征（VPD/SM 的均值和变异系数 CV，以及它们在两个时段间的变化 Δ）去解释它们的空间差异。为了降低空间自相关和噪声，我们先把 0.5° 面积加权聚合成 2.5°（1,761 个块，05b_coarsen.py），再跑一次 AutoML + TreeSHAP（06_attribution.py）。四个目标最终都选中了 XGBoost：

目标	test R²	前三位驱动
D	0.23	SM 均值 43.5%、VPD CV 27.3%、VPD 均值 18.5%
δRI_SM	0.43	ΔCV_SM 29.0%、ΔCV_VPD 19.1%、SM 均值 11.8%
δRI_VPD	0.55	ΔVPD 均值 28.1%、ΔSM 均值 15.0%、ΔCV_VPD 15.0%
ΔD	0.59	ΔVPD 均值 21.3%、ΔCV_SM 17.7%、ΔSM 均值 17.7%

结论很直观：**VPD 均值涨得越多、土壤变得越干的地方，VPD 的控制权就增强得越多。**而且起作用的不只是均值，变率（CV）的变化同样排在前列，对应原文"均值和变率共同改变植物水分胁迫响应"的说法。

七、SHAP 可视化

SHAP 图是这类文章的门面。同一份 δRI_VPD 归因结果，我们画了三种风格：

条形重要性 + 蜂群 + LOWESS 依赖曲线

07_shap_style.py：左边是 mean|SHAP| 条形图，右边是按特征值着色的散点和 LOWESS 趋势。这种排版信息密度高，适合在正文里做主图。

SHAP可视化风格

09_shap_final.py 融合了前两版，按 IPCC AR6 综合报告的版式来画：

标题直接写结论（VPD importance grows most where VPD rises most and soils dry most）；
a 栏：重要性条形；b 栏：格点级蜂群，用箱体标出 25–75% 区间，并加了注释箭头；c 栏：用分位数做横轴的响应曲线，带 25–75% 带；
底部的分组色带显示 VPD 组和 SM 组的总贡献（60.3% vs 39.7%）；
配色经过色盲友好（CVD）校验。

另外三个目标的最终版：

八、几点复现心得

传感器切换要不要担心？ 1982–2022 横跨 AVHRR 和 MODIS 两代传感器。我们检查了 NDVI 残差标准差：MODIS 时代约为 AVHRR 时代的 1.76 倍，但这个比值与 ΔD 的 Spearman 相关只有 −0.004，可以排除 ΔD 格局是传感器切换造成的。
R² 门槛：逐格点模型的平均 test R² 约 0.39，用 R² ≥ 0.2 筛选后再解释，否则 SHAP 解释的就是噪声。
第二级 D 模型 R² 只有 0.23，解读"背景态"驱动要谨慎；变化量模型（R² 0.43–0.59）更可靠。
AutoML 的价值：3.4 万个格点不可能手动调参，FLAML 在 8 秒预算内自动在 RF 和 XGBoost 之间做选择，加上 TreeSHAP 的精确解，让"每个格点一个可解释模型"变得可行。

九、代码与数据

全部代码、结果图（PNG/PDF/SVG）、SHAP 绘图数据、Excalidraw 手绘流程图源文件都在：

https://github.com/GISWLH/vpd-sm-automl-xai^[3]

原始 ERA5、NDVI 和 MODIS 数据体量大，没有放进仓库，脚本里写明了各自的来源和处理方式。欢迎 Star、Fork，也欢迎拿去换成自己的研究区和变量。

参考文献：Li L. et al. Vegetation increasingly controlled by atmospheric dryness rather than soil moisture. Nature Geoscience, 2026. doi:10.1038/s41561-026-02107-w

END

声明： 欢迎转载、转发。气象学家公众号转载信息旨在传播交流，其内容由作者负责，不代表本号观点。文中部分图片来源于网络，如涉及内容、版权和其他问题，请联系小编(微信：qxxjgzh)处理。