基于SHAP方法识别气象驱动因子对水稻单产的非线性阈值及交互主导性——以宁波地区单季稻为例

《Atmosphere》：Identifying Nonlinear Thresholds and Interaction Dominance of Meteorological Drivers on Rice Yield: A SHAP-Based Approach

　　

量化作物系统对气象驱动因子的非线性响应是农业气象学的核心难题。尽管可解释人工智能（Explainable Artificial Intelligence, XAI）提供了新途径，但现有基于SHAP（Shapley Additive exPlanations）的阈值识别方法多局限于效应方向反转点的判定，且缺乏因子间交互效应的统一量化分级标准。为探究影响水稻单产的气象因子阈值及交互作用强度，研究人员利用宁波市9个区县1995—2024年水稻单位面积产量与气象数据，基于LASSO及6种机器学习算法构建水稻产量预测模型，采用基于SHAP算法的递归特征消除（Recursive Feature Elimination, RFE）筛选出11个核心气象因子。在此基础上提出两项创新指标：一是引入导数极值阈值（Derivative Extrema Threshold, DET）作为零穿越阈值（Zero-Crossing Threshold, ZCT）的补充，通过对平滑后SHAP依赖图曲线一阶导数的极值点定位，识别效应强度发生质变但未发生方向逆转的临界位置；二是提出交互主导比（Interaction Dominance Ratio, IDR），将交互变异度归一化至总效应框架内，建立强、中、弱三级交互分级标准。研究发现特征优化后的LightGBM模型预测效果最优（R2= 0.833）。ZCT识别出置信区间极窄的方向逆转点，如8月累计降水量210.6 mm、6月平均气温24.5 ℃；DET进一步揭示了效应强度突变特征，如5月累计降水量64.9 mm处"减产效应减弱"。IDR分析准确刻画了由8—9月平均气温、6月最低气温及8月累计降水量构成的"交互主导三角网络"。本研究构建的"ZCT?DET?IDR"框架推动了从方向判别、强度预警到多因子协同分析的范式进阶，为区域农业气象灾害精细化预警提供了可量化的新视角。

论文解读：基于SHAP方法识别气象驱动因子对水稻单产的非线性阈值及交互主导性

研究背景与意义

水稻是我国最重要的口粮作物之一，其稳产直接关系到国家粮食安全。准确预测水稻产量并阐明气象因子对产量的影响机制，对指导农业生产及制定防灾减灾措施具有重要意义。作物生长发育及产量形成本质上是非线性自然过程，传统统计回归或相关分析模型多局限于特定时空条件，难以处理复杂的非线性问题；而机器学习与深度学习等非线性预测模型虽显著提升了精度，却因"黑箱"特性限制了对决策过程的理解，阻碍其在产量预测中的深入应用。可解释人工智能（Explainable Artificial Intelligence, XAI），特别是SHAP（Shapley Additive exPlanations）方法，为打开模型黑箱提供了理论依据。现有基于SHAP依赖图（Dependence Plot）的阈值识别多采用零穿越阈值（Zero-Crossing Threshold, ZCT），仅能捕捉效应方向由正转负或由负转正的关键点，无法识别效应方向不变但强度发生质变（如增产效应由快速上升转为饱和、减产效应加速恶化）的位置，且缺乏因子间交互效应的归一化定量分级标准。为此，研究人员以宁波市1995—2024年单季稻为研究对象，构建"ZCT?DET?IDR"耦合框架，系统阐明气象因子对水稻产量的非线性影响及协同机制，该研究成果发表于《Atmosphere》。

主要关键技术方法

研究人员选用宁波市9个区县1995—2024年水稻单位面积产量数据（来源于《宁波统计年鉴》）及同期对应站点逐日气象观测资料（平均气温、最高气温、最低气温、相对湿度、降水量、最大风速），将实际产量分解为趋势产量与气象产量（采用Hodrick–Prescott滤波器提取），对6类基础气象变量进行时间窗口扩展生成216个时序特征，经LASSO初步筛选保留15个特征。分别以多元线性回归（Multiple Linear Regression, MLR）及支持向量回归（Support Vector Regression, SVR）、Bagged Trees、随机森林（Random Forest, RF）、反向传播神经网络（Back Propagation Neural Network, BPNN）、LightGBM构建气象产量预测模型，采用基于时间序列分割的5折交叉验证网格搜索/Bayesian优化调参，按时间顺序划分训练集（1995—2018）与测试集（2019—2024）。以基于SHAP值均值绝对大小的递归特征消除（Recursive Feature Elimination, RFE）确定最终11个核心气象因子输入LightGBM最优模型。在此基础上：(1) 对SHAP依赖图做移动平均平滑，通过相邻点SHAP值异号的线性插值计算零穿越阈值（ZCT），并以Bootstrap（n=500）估计95%置信区间；(2) 对平滑曲线三次样条拟合求一阶导数，取其一阶导数列的局部极值点为导数极值阈值（Derivative Extrema Threshold, DET），并结合箱线图四分位距（InterQuartile Range, IQR）过滤分布尾部伪阈值；(3) 基于SHAP交互值计算交互主导比（Interaction Dominance Ratio, IDR）＝(Φ_ij的P₉₀–P₁₀) / 特征对总效应的标准差，按IDR＜0.30（弱交互）、0.30≤IDR＜0.50（中等交互）、IDR≥0.50（强交互）三级标准分级，并构建"绝对交互幅度–IDR"二维分类空间。

研究结果

3.1 初始数据分析

研究人员经HP滤波分解得到气象产量变幅为?2.12～+1.49 t·hm^?2，极端减产区与增产区不对称，表明气象灾害对产量的负面冲击大于有利天气的正向增益。11个气象因子与气象产量的Pearson相关系数均低于|r|＜0.08且不显著（p＞0.05），证实单因子与产量间呈弱线性关系，需采用非线性机器学习建模。温度因子变异系数（Coefficient of Variation, CV）低且分布紧凑，降水因子CV高且右偏，说明温度是生长条件基线，降水是年产波动主因，亦解释温度衍生阈值置信区间较降水更窄。

3.2 不同产量模型预测性能比较

六类模型在测试集上的评估显示LightGBM表现最优（R²=0.809，RMSE=0.352 t·hm^?2，MAE=0.258 t·hm^?2），优于SVR（R²=0.774）、Bagged Trees和RF，BPNN（R²=0.709）与MLR（R²=0.641）精度最低。梯度提升树框架在处理气象因子—水稻产量非线性关系上具明显优势，选定LightGBM开展后续分析。

3.3 模型优化与气象因子贡献

基于SHAP重要性值的RFE过程表明保留11个特征时5折时间序列交叉验证R²达峰值，剔除4个冗余特征后优化LightGBM在独立测试集上R²提升至0.833（RMSE=0.330 t·hm^?2，MAE=0.252 t·hm^?2）。SHAP蜂群图（Beeswarm Plot）与饼图显示核心驱动因子依次为：8月累计降水量R2020_8（相对贡献17.4％）、3月平均气温TAVG_3（13.1％）、6月平均气温TAVG_6（—降至第5位）、8—9月平均气温TAVG_8-9（11.0％）、6月最低气温TMIN_6（升至第3位，10.9％），前四位累计贡献约51.5％，降水和春夏气温共同主导宁波水稻产量预测。

3.4 水稻产量气象阈值识别

3.4.1 ZCT识别结果

11个因子SHAP依赖图均呈典型单零点穿越行为。ZCT分别为：TAVG_3为11.6 ℃，TAVG_6为24.5 ℃，TAVG_8-9为26.2 ℃，TMIN_6为21.6 ℃，R2020_8为210.6 mm，R2020_5为113.1 mm等。置信区间宽度显示温度因子属窄过渡带（如TMIN_6仅0.05 ℃），降水因子属宽过渡带（如R2020_4-7达26.08 mm），反映土壤—作物系统的水分调节缓冲能力。窄过渡带因子可设明确单值预警阈值，宽过渡带宜设过渡预警区间。

3.4.2 DET识别结果

经IQR过滤后各因子DET点显示：TAVG_8-9的DET位于ZCT右侧（26.86～26.88 ℃），SHAP均值为显著负值，代表乳熟期高温超限后负效应加速恶化区；R2020_5的DET位于ZCT左侧（64.9～79.5 mm），表征移栽返青期干旱缓解加速区（减产效应减弱）；TAVG_3具双侧DET分布。长时段累计降水变量DET落于数据稀疏极端尾部者视为样条外推不确定标记予以说明，不赋予农学意义。DET弥补了ZCT无法捕获非方向逆转性质变点的不足。

3.5 因子交互效应分析

SHAP交互矩阵分解显示11个因子纯主效占比76.99％～98.06％，但交互效应占比存在差异，TMIN_6与R2020_8存在较明显效应抵消（分别约9.16％和5.58％），说明传统均值绝对SHAP会低估此类因子真实潜力。IDR热图表明55个因子对中仅3对（5.5％）达强交互主导（IDR≥0.50）：TAVG_8-9×TMIN_6（IDR=0.622）、TAVG_8-9×R2020_8（IDR=0.549）、TMIN_6×R2020_8（IDR=0.525），三者恰好构成以TAVG_8-9、TMIN_6、R2020_8为顶点的交互主导三角网络（Interaction Dominance Triangular Network）；6对（10.9％）为中等交互，46对（83.6％）为弱交互。二维分类空间中TAVG_8-9×TMIN_6落入高IDR–高绝对幅度象限，为最具农学意义的协同组合；TAVG_6-7×R2020_8虽绝对交互强度排第三但IDR=0.390（中等），因两因子分属不同生育阶段缺乏条件依赖性，以主效驱动为主。三维交互曲面图直观印证上述差异。

讨论与结论总结

讨论部分指出：ZCT与本地农气实践吻合——3月均温11.6 ℃接近宁波早稻安全播种温度指标（10～12 ℃），低于此芽腐苗死率升高；8—9月均温26.2 ℃与8月降水210.6 mm构成交互主导三角节点，对应宁波单季稻抽穗—灌浆期高温（日均温≥27 ℃降低结实率与千粒重）叠合台风致持续淹水（月降水＞200 mm引发田间渍害）的复合胁迫；10月降水83.3 mm与晚稻灌浆期秋雨致根系缺氧相符。DET揭示ZCT之外的效应强度演变特征（如高温超阈值后无缓冲加速恶化、5月适中降水使减产效应衰减）。IDR量化交互主导程度，发现6月低温背景调制8—9月高温敏感性（物候窗口偏移）及高温—过湿并发放大损失，提示传统单因子阈值预警存在系统性低估风险，应构建含前期温度异常修正项的多因子复合气象指数。方法学上ZCT?DET?IDR框架较传统偏依赖图（Partial Dependence Plots, PDP）/累积局部效应（Accumulated Local Effects, ALE）及普通SHAP依赖图更能定位强度质变点与归一化交互强度，但需注意DET对样条平滑参数敏感、IDR的分级阈值借用了Cohen's效应量大中小基准属类比推广、SHAP归因为条件相关性而非因果机制、结果需在受控试验中外推验证及跨区域校准。研究局限含时间窗口扩展致多重共线性仅部分消除、DET稀疏区伪阈值说明、模型R²提升未做Bootstrap显著性检验、未做分区交互异质性分析及未做跨作物跨气候区外部验证。

结论：经LASSO与SHAP?RFE双重筛选，含11个核心气象因子的LightGBM模型对宁波地区水稻产量预测效果最优（R²=0.833）。提出的导数极值阈值（DET）弥补了传统零穿越阈值（ZCT）仅能识别方向逆转点的局限，成功提取效应强度突变特征；交互主导比（IDR）实现交互效应强度的归一化横向分级，准确刻画了由8—9月平均气温、6月最低气温与8月累计降水量组成的交互主导三角网络。"ZCT?DET?IDR"框架识别的核心阈值（3月均温11.6 ℃、8—9月均温26.2 ℃、8月降水210.6 mm等）与宁波本地农气试验及灾情记录高度一致，证明所识别阈值对应具明确物理与农学意义的灾害发生边界。该框架为区域农业气象灾害精细化预警提供了兼具数学严谨性与实用性的决策基础，且具有跨XAI解释场景的可移植潜力，未来需结合因果推断与多源数据集完善为标准化的"阈值?交互"双指标XAI工具箱。