机器学习干旱预测实战：从SPI计算到XGBoost建模评估

简介面向气候科学的机器学习管道源自ml_drought项目为研究干旱预测及相关气候问题提供端到端解决方案。资源面向具备一定Python基础、希望将机器学习引入气候数据分析的研究者与开发者包含创建、相互比较和评估机器学习方法的完整流程内置多个任务类统一数据格式并支持自定义模块。压缩包约49.31MB内容以Python源码、Jupyter Notebook示例、环境配置文件如environment.yml为主按src目录划分功能类并提供三个入口点便于快速上手。管道设计清晰模块划分合理便于二次开发。通过配置环境文件可一键创建conda环境免去繁琐依赖安装。已有141人学习适合作为入门气候机器学习的实战参考。通过该资源可掌握从数据预处理、特征工程到模型评估的管道搭建思路并可基于示例快速构建自己的干旱预测研究工作。1. ml_drought 是什么把干旱预测变成监督学习问题而不是继续堆物理模型做干旱预测的人过去一提到改进第一反应是换更复杂的陆面过程模式、加更多的同化数据。但 ml_drought 这个项目走的是完全相反的路子把干旱预测当成一个监督学习问题用历史气象观测和再分析数据训练机器学习模型让模型自己从降水、温度、蒸散、土壤湿度的时空序列里找出干旱发生的规律。它的核心贡献不只是“用 ML 替代传统模型”而是把“预测”和“理解”两件事放在一起做。这个方案解决的是实际业务里最头疼的两件事一是短期干旱预警的准确率上不去传统统计模型和物理模式在 1 到 3 个月的尺度上经常翻车二是模型即使预报对了你也说不清是哪个因子在起主导作用。ml_drought 的做法是用可解释性工具把黑匣子打开告诉你这次干旱预测是海温异常驱动的还是本地土壤缺水驱动的。适合正在做干旱监测预警、农业保险定价、水资源调度决策的从业者也适合想上手“气象 机器学习”这个方向的开发者。仓库在 GitHub 的 ml-clim 组织下代码和数据管线是完整开源的照着跑一遍就能在自己的数据集上复现。2. 从气象数据到干旱标签SPI 计算与数据集构建2.1 为什么先算 SPI干旱的定义本身就是标签工程做监督学习第一步是定标签。干旱和降雨预测不同它不是“明天下不下雨”这种直接观测得到的事件而是一个基于统计定义的缓慢过程。最常见的做法是用标准化降水指数 SPI把某段时间的累计降水量放到历史分布里做标准化负值代表偏旱-1 以下算中度干旱-1.5 以下算重度干旱。ml_drought 选择 SPI 而不是 PDSI 或土壤湿度百分位有几个实际考虑。SPI 只需要降水数据不依赖复杂的陆面参数它的时间尺度可调SPI-3 反映季度尺度水分亏缺SPI-12 反映长期干旱而且 SPI 在气象学、水文学界被广泛使用业务系统愿意认。把预测目标定为“未来某个月的 SPI-3 是否低于 -1”这个问题就从回归变成了二分类和下游预警流程衔接更顺。SPI 的计算本身不是随便调个库就完事的。零降水月份在干旱区非常常见而 Gamma 分布不允许输入为 0。很多新手在这里翻车直接对含零的序列做拟合得到的分布参数是错的。我一般把序列拆成两部分处理零降水月份用经验频率单独计算正降水月份用 Gamma 拟合 CDF然后按概率加权合并最后做正态逆变换。2.2 用降水与温度数据计算 SPI-3 / SPI-6核心脚本这里给一份可以直接跑的 SPI 计算脚本输入是某站点或某格点过去 30 年以上的月降水序列。数据源用 ERA5 或者 CPC 全球降水产品都行先把 NetCDF 转成 pandas DataFrame一列是时间一列是降水。import numpy as np import pandas as pd from scipy.stats import gamma, norm def aggregate_precip(monthly_precip: pd.Series, timescale: int) - pd.Series: 把月降水滚动聚合到目标时间尺度例如 SPI-3 就是近 3 个月累计 return monthly_precip.rolling(windowtimescale, min_periodstimescale).sum() def fit_gamma_params(precip_rolled: pd.Series): 对滚动累计降水做 Gamma 拟合。 只取非零值拟合零值单独计算经验频率。 vals precip_rolled.dropna() positive vals[vals 0] # floc0 固定位置参数为 0Gamma 分布只允许大于 0 的输入 shape, loc, scale gamma.fit(positive, floc0) p_zero (vals 0).mean() return shape, loc, scale, p_zero def calc_spi(precip_rolled: pd.Series, params) - pd.Series: 用拟合好的 Gamma 分布把累计降水映射为标准正态值。 零降水的 CDF 直接取 p_zero非零值用 Gamma CDF 换算。 shape, loc, scale, p_zero params cdf_zero p_zero # 非零部分Gamma CDF 乘以 (1 - p_zero)加上零值概率 cdf cdf_zero (1 - cdf_zero) * gamma.cdf(precip_rolled, shape, locloc, scalescale) # 避免 log(0) 和正态逆变换的极端值 cdf np.clip(cdf, 1e-6, 1 - 1e-6) return pd.Series(norm.ppf(cdf), indexprecip_rolled.index) # 使用示例 precip pd.Series(...) # 你的月降水序列index 为 DatetimeIndex rolled_3 aggregate_precip(precip, timescale3) params_3 fit_gamma_params(rolled_3) spi3 calc_spi(rolled_3, params_3)这段代码的关键在于fit_gamma_params里的参数处理。gamma.fit(positive, floc0)固定了 Gamma 分布的位置参数如果放开floc让它自由拟合分布会被零值拉偏SPI 结果在干旱端会失真。另一个注意点是rolling(min_periodstimescale)如果序列中间有缺失月份不要自动填充宁可让这一段不参与拟合否则滚动聚合会把缺失值“吞掉”变成错误的累计值。2.3 组建时序样本集滑动窗口、滞后特征与训练/验证/测试切分SPI 算出来后下一步是把数据组织成监督学习格式每个样本是“过去 N 个月的观测特征 未来 M 个月的目标标签”。这里有两个关键约束。第一个约束是特征和标签必须严格错开时间。举个例子如果用 t 年 1 月到 3 月的降水特征去预测 t 年 1 月到 3 月的 SPI-3那就是在泄漏答案——SPI-3 本身就是这三个月累计降水算出来的。ml_drought 这类项目里最常见的错误就是把同期变量当特征用。正确的是用 t-1、t-2、t-3 个月的降水、温度、SPI 等做特征预测 t 月及以后的 SPI。第二个约束是切分方式必须按时间顺序不能随机 shuffle。具体做法我会在第 4 章展开这里先给出样本构建的思路def build_samples(spi_series: pd.Series, precip: pd.Series, temp: pd.Series, feature_window: int 12, lead_time: int 1): 用过去 feature_window 个月的观测构造特征 预测 lead_time 个月后的 SPI-3 是否为负干旱事件。 features, labels, timestamps [], [], [] for t in range(feature_window, len(spi_series) - lead_time): # 特征窗口t-feature_window 到 t-1 f np.concatenate([ precip.values[t-feature_window:t], temp.values[t-feature_window:t], spi_series.values[t-feature_window:t] ]) # 标签未来第 lead_time 个月是否发生干旱 y 1 if spi_series.values[t lead_time] -1.0 else 0 features.append(f) labels.append(y) timestamps.append(spi_series.index[t]) return np.array(features), np.array(labels), timestampslead_time是核心业务参数它决定了你是做“下个月是否干旱”还是“三个月后是否干旱”。通常 lead_time 越大预测难度越高因为大气过程的可预报性在 2 周以上就快速衰减。实际项目里训练集覆盖 30 到 50 年就够用了过长的历史序列反而可能引入气候背景不一致的样本比如 1960 年的温度平均值和 2010 年完全不同模型会去拟合这种不想要的趋势。3. 特征工程与模型选型为什么树模型在干旱预测里常比深度学习更稳3.1 特征怎么造滞后降水、温度异常与大尺度气候指数标签定了决定模型上限的是特征。ml_drought 的思路和传统统计预报很接近只是不再手工指定“哪个月降水最重要”而是把一组有物理依据的特征丢给模型自己去学。我通常把特征分成三层。第一层是本地变量包括过去 12 个月的逐月降水、平均气温、SPI-3 和 SPI-12 历史值。SPI 的滞后值尤其重要因为干旱有持续性上个月的 SPI-3 对下个月有很强的指示意义。第二层是季节性信息把月份编码成周期变量让模型知道“这是在季风季节还是干季”。第三层是大尺度气候指数比如 Nino 3.4 海温异常、印度洋偶极子 IOD、北大西洋涛动 NAO这些指数是跨区域干旱的重要驱动因子尤其是 ENSO 对东南亚、澳大利亚、非洲南部干旱的影响非常显著。特征工程的粒度上也值得推敲。站点级预测用逐点观测没问题但如果是区域网格预测我一般会对周边格点做空间聚合比如取目标点周围 2 度范围内的平均降水而不是只取单格点。这样做有两个好处一是平滑掉单格点的噪声二是让模型学到空间梯度信息。3.2 XGBoost 与 LSTM 的选型边界不是越深越好很多人一看到“机器学习预测干旱”就默认要上 LSTM、Transformer。我做过对比实验结论是在月尺度的干旱预测任务上经过调参的 XGBoost 通常能打赢结构复杂的深度学习模型而且训练成本低一个数量级。原因不复杂——月尺度气象数据的样本量本身就少一个站点 50 年逐月数据只有 600 个样本切掉验证集和测试集后训练样本可能不到 400 条这种规模撑不起大模型的容量。LSTM 有它的价值场景当输入是日尺度或小时尺度的连续观测序列长度在数百步以上且你需要捕获降水的日内演变过程时序列模型有优势。但 ml_drought 这类以 SPI 为目标的任务输入已经被滚动聚合成了月尺度特征时间依赖信息其实已经被 SPI 的滞后项吸收得差不多了再用序列模型属于重复建模。我一般会跑一个轻量对比来确定选型同样的特征XGBoost、随机森林、逻辑回归、LSTM 各出一版结果比较验证集上的 CSI 和 AUC。如果 LSTM 没有明显超过 XGBoost 5% 以上生产环境就选 XGBoost理由是可解释性好、调参快、不容易在样本少的时候过拟合。3.3 基准必须做持续预测和气候态预测是两条及格线在评估任何 ML 模型之前先跑两个基准。第一个是持续预测persistence直接用上个月的 SPI-3 作为下个月的预测值。因为干旱有惯性这个朴素基准的成绩往往不低如果模型打不过它说明学到的东西还不如一个“干旱会延续”的假设。第二个是气候态预测climatology用历史同期的 SPI 平均值作为预测这个基准检验的是季节性规律的作用。这两个基准不是摆设它们决定了你的 ML 模型有没有实际增量。我在一个区域水旱风险项目里就碰到过这种情况XGBoost 的 AUC 做到 0.78看着不错但 persistence 基准是 0.74模型的真实增益只有 4 个百分点。这个结果不能说没有价值但它提醒你收敛预期也提醒你把精力花在找真正有效的新特征上而不是继续调参。4. 训练与评估最小闭环TimeSeriesSplit、早停与 CSI 指标4.1 训练代码用 TimeSeriesSplit 避免“用未来预测过去”时序数据切分有个容易忽略的细节sklearn的train_test_split默认随机切分这对时间序列是致命的。随机切分会让训练集里出现测试集时间附近的样本由于气象序列有强自相关模型相当于提前看过了“答案附近的天气状态”验证集的分数会虚高 0.1 到 0.2 个 AUC。正确的做法是用TimeSeriesSplit它保证训练集永远在验证集之前而且可以通过gap参数在训练集和验证集之间留出空隙避免训练集末尾的样本与验证集开头样本因为时间相邻而产生信息流向。import xgboost as xgb import numpy as np from sklearn.model_selection import TimeSeriesSplit X, y ..., ... # 来自 build_samples 的输出 tscv TimeSeriesSplit(n_splits5, gap6) # gap6 表示训练/验证之间空出 6 个月 cv_scores [] for fold, (train_idx, valid_idx) in enumerate(tscv.split(X)): X_train, X_valid X[train_idx], X[valid_idx] y_train, y_valid y[train_idx], y[valid_idx] # 处理干旱样本不均衡正样本权重与负样本比例对齐 pos_weight np.sum(y_train 0) / np.sum(y_train 1) model xgb.XGBClassifier( n_estimators2000, learning_rate0.02, max_depth4, subsample0.8, colsample_bytree0.8, scale_pos_weightpos_weight, eval_metricauc, early_stopping_rounds50, random_state42 ) model.fit( X_train, y_train, eval_set[(X_valid, y_valid)], verboseFalse ) best_iter model.best_iteration pred_prob model.predict_proba(X_valid)[:, 1] # 这里先存概率CSI 等指标统一在评估阶段计算 cv_scores.append((y_valid, pred_prob, best_iter)) print(各折最佳迭代轮数:, [s[2] for s in cv_scores])这段代码有两个参数值得单独说明。gap6是给训练集和验证集之间留出 6 个月的缓冲。为什么需要这个缓冲因为 SPI 的自相关长度会跨好几个月训练集最后一年的状态会通过滞后特征影响到验证集第一年的标签gap 能切断这条泄漏路径。scale_pos_weight的作用是把干旱样本的损失权重调高因为干旱事件本来就比较少如果不平衡设置模型会倾向把所有样本都预测为“不干旱”整体准确率很高但对预警业务毫无用处。4.2 评估指标CSI、POD、FAR 才是业务关心的数字分类准确率在干旱预警里没有意义。假设干旱事件只有 15% 的频率模型什么都不做全预测“不干旱”准确率也是 85%但决策者拿到这个模型等于没有模型。气象预警领域通行的指标是 CSICritical Success Index、PODProbability of Detection和 FARFalse Alarm Ratio。POD 是“干旱事件被正确预报出来的比例”漏报率过高会让用户不信任模型。FAR 是“预报了干旱但实际没发生”的比例空报率过高会让决策者麻木狼来了喊多了就失效。CSI 把漏报和空报都纳入考虑分数越高越好。下面这段代码用验证集概率和一组阈值计算这三个指标。def calc_skill_scores(y_true: np.ndarray, y_prob: np.ndarray, threshold: float 0.5): pred (y_prob threshold).astype(int) hit np.sum((pred 1) (y_true 1)) miss np.sum((pred 0) (y_true 1)) false_alarm np.sum((pred 1) (y_true 0)) correct_negative np.sum((pred 0) (y_true 0)) pod hit / (hit miss) if (hit miss) 0 else 0.0 far false_alarm / (hit false_alarm) if (hit false_alarm) 0 else 0.0 csi hit / (hit miss false_alarm) if (hit miss false_alarm) 0 else 0.0 return {POD: pod, FAR: far, CSI: csi} # 对第 4.1 节的每个 fold 计算指标 for y_valid, pred_prob, _ in cv_scores: scores calc_skill_scores(y_valid, pred_prob, threshold0.3) print(scores)注意这里我把阈值设成了 0.3 而不是默认的 0.5。原因是干旱样本比例低模型输出的概率普遍偏低卡 0.5 会让漏报率高得无法接受。阈值怎么选要看业务权衡防汛抗旱部门宁可信其有不可信其无阈值就调低接受更高的空报率换取更低的漏报率农业保险定价相反需要压低 FAR。这个阈值不应该拍脑袋而是遍历 0.1 到 0.6画出 POD 和 FAR 的权衡曲线再定。4.3 时间序列回检模型有没有滞后一眼就能看出来指标算完还不能直接交付必须把预测结果拉回时间轴上看一眼。做法很简单把验证集的预测概率按时间顺序画成曲线和真实的 SPI 序列叠在一起。这个可视化的价值在于它能暴露两类指标看不出的问题。第一类问题是系统性滞后。如果模型给出的干旱概率峰值总是比真实 SPI 谷值晚一两个月说明模型过度依赖 SPI 滞后项基本就是把 persistence 抄了一遍。这时候需要检查特征里是否加入了当期变量或者 lead_time 太短导致模型“来不及反应”。第二类问题是季节错位。如果模型在某个固定季节频繁报高概率比如每年 9 月都报干旱但实际只有一半年份干旱说明模型把季节性当成了干旱信号需要检查特征里月份编码的影响是否过大。顺手贴一个检查滞后的代码思路import matplotlib.pyplot as plt # pred_prob 已经按时间排好序 plt.figure(figsize(12, 4)) plt.plot(timestamps_valid, spi_valid, labelObserved SPI-3, colorblack) plt.plot(timestamps_valid, pred_prob, labelPredicted drought prob, colorred) plt.axhline(-1.0, linestyle--, colorgray) plt.legend() plt.show()曲线贴合不代表好如果预测概率曲线看起来像 SPI 曲线的左右平移说明模型只是在延迟复述观测。好的预测是概率峰值出现在 SPI 谷值之前这才是“预测”而不是“跟唱”。5. 四个必踩的坑数据泄漏、样本不平衡与气候外推5.1 坑一随机切分导致“用未来预测过去”验证集虚高现象用默认train_test_split训练模型验证集 AUC 高达 0.85CSI 也很漂亮。换成TimeSeriesSplit重跑AUC 掉到 0.74直接缩水 0.1 以上。原因SPI 序列本身有很强的自相关。随机切分时验证集中的样本和训练集中的某些样本在时间上只差一两个月它们的气象状态高度相似。模型相当于在考试时看到了“同类型题目的标准答案”。这不是模型变强了是信息从训练集流向了测试集。解决所有涉及时间的过程都必须用TimeSeriesSplit并且在折与折之间设置gap。gap 的大小参考 SPI 的自相关长度月尺度数据我一般至少空 3 到 6 个月。如果数据是日尺度gap 需要对应扩大到 90 天以上。5.2 坑二把同期观测数据放进特征模型在“读答案”现象模型对训练集和验证集的拟合都好得不可思议AUC 接近 0.95但换到新数据上一塌糊涂。原因这是最隐蔽的泄漏。假设你要预测 2024 年 3 月的 SPI-3特征里却包含了 2024 年 3 月的降水量。SPI-3 本来就是 1 月到 3 月降水的函数模型根本不需要学任何气象规律直接根据 3 月降水反推 SPI 即可。我自己早期做类似项目时也在这上面吃过亏当时的“高精度”模型其实就是个查表器。解决严格检查每个特征的“信息截止时间”。构造特征时只能使用t - lead_time之前的数据。一个简单粗暴的验证方法是把目标月的降水特征从特征列表里删掉重新跑一遍训练如果成绩大幅下降别急着高兴想想删掉的到底是什么。5.3 坑三干旱样本太少模型变成“永远不报干旱”的复读机现象验证集准确率 85%但 POD 接近 0模型一个干旱事件都没报出来。原因SPI 在 -1 以下的样本比例理论上只有 15.9%重度干旱更是稀罕事。类别不平衡时模型学到的最优策略是全部预测为负类因为这样总损失最小。解决先算正负样本比用scale_pos_weight或者class_weight调整损失权重。注意不要用 SMOTE 这类过采样方法因为气象时间序列的合成样本很可能违背天气演变的物理过程造出“1 月高温 2 月极端降水”这种现实中不会出现的组合。更稳妥的做法是调整决策阈值把阈值从 0.5 降到 0.3 甚至 0.2用更高的 FAR 换回 POD。5.4 坑四气候非平稳性让训练集和测试集分布不一致模型外推失效现象模型在验证集上成绩不错但部署后第一年效果尚可第二年开始频繁空报。原因气候变化让训练期的气候态和现在的气候态已经不同。比如模型在 1990 到 2015 年的数据上学到的“温度异常与干旱的关系”在 2020 年后全球升温背景下被打破了。模型的输入分布发生漂移输出自然不可靠。解决没有完美的答案但有缓解手段。一是把大尺度气候指数和长期趋势项作为特征给模型提供气候背景信息二是训练时不要用全部历史数据而是做滚动训练每 1 到 2 年用最新数据重新拟合一次模型给旧样本设衰减权重三是监控特征分布漂移比如在周度、月度的模型评估中加上 PSIPopulation Stability Index如果 PSI 超过 0.2就要安排重训。6. 用 SHAP 打开黑匣子把预测结果接进预警流程指标告诉你模型“准不准”SHAP 告诉你模型“看到了什么”。ml_drought 项目把可解释性和预测并列这个设计不是锦上添花——业务方真正要的是一个能辅助决策的预报依据不是一个只吐概率的黑盒子。用 SHAP 的 TreeExplainer 可以直接对 XGBoost 输出归因值import shap explainer shap.TreeExplainer(model) shap_values explainer.shap_values(X_valid) # 全局特征重要性按平均绝对值排序 shap.summary_plot(shap_values, X_valid, feature_namesfeature_names)我常用的两个图是 summary plot 和 dependence plot。前者能看到哪些特征对干旱预测影响最大通常 SPI-3 的滞后值和 Nino 3.4 海温异常会排在最前面这符合物理预期。后者看单个特征的边际效应比如 Nino 3.4 指数异常升高时模型输出的干旱概率如何变化能直接验证模型是否学到了“厄尔尼诺导致东南亚干旱”这类已知机制。把解释性接进业务流程的做法是每次预报输出时附带 top-3 贡献特征。比如“2024 年 3 月干旱概率 68%主要由 SPI-3 持续偏低、Nino 3.4 异常、过去 3 个月降水低于历史第 20 百分位驱动”。这种格式的预报反馈比单一概率值更容易让跨部门的人接受也方便事后复盘模型到底哪里看错了。我现在做这类项目拿到新数据集的第一件事一定是先画 SPI 时间序列和特征分布图而不是急着跑模型——数据要是脏的模型再花哨也是白搭。希望这些从数据到评估再到解释的流程和踩坑记录能帮到你。本文还有配套的精品资源点击获取