☰
机器学习干旱预测实战:从SPI计算到XGBoost建模评估
2026/9/25 5:24:43 网站建设 项目流程

简介:面向气候科学的机器学习管道,源自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(window=timescale, min_periods=timescale).sum() def fit_gamma_params(precip_rolled: pd.Series): """ 对滚动累计降水做 Gamma 拟合。 只取非零值拟合,零值单独计算经验频率。 """ vals = precip_rolled.dropna() positive = vals[vals > 0] # floc=0 固定位置参数为 0,Gamma 分布只允许大于 0 的输入 shape, loc, scale = gamma.fit(positive, floc=0) 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, loc=loc, scale=scale) # 避免 log(0) 和正态逆变换的极端值 cdf = np.clip(cdf, 1e-6, 1 - 1e-6) return pd.Series(norm.ppf(cdf), index=precip_rolled.index) # 使用示例 precip = pd.Series(...) # 你的月降水序列,index 为 DatetimeIndex rolled_3 = aggregate_precip(precip, timescale=3) params_3 = fit_gamma_params(rolled_3) spi3 = calc_spi(rolled_3, params_3)

这段代码的关键在于fit_gamma_params里的参数处理。gamma.fit(positive, floc=0)固定了 Gamma 分布的位置参数,如果放开floc让它自由拟合,分布会被零值拉偏,SPI 结果在干旱端会失真。另一个注意点是rolling(min_periods=timescale),如果序列中间有缺失月份,不要自动填充,宁可让这一段不参与拟合,否则滚动聚合会把缺失值“吞掉”变成错误的累计值。

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), timestamps

lead_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_splits=5, gap=6) # gap=6 表示训练/验证之间空出 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_estimators=2000, learning_rate=0.02, max_depth=4, subsample=0.8, colsample_bytree=0.8, scale_pos_weight=pos_weight, eval_metric="auc", early_stopping_rounds=50, random_state=42 ) model.fit( X_train, y_train, eval_set=[(X_valid, y_valid)], verbose=False ) 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])

这段代码有两个参数值得单独说明。gap=6是给训练集和验证集之间留出 6 个月的缓冲。为什么需要这个缓冲?因为 SPI 的自相关长度会跨好几个月,训练集最后一年的状态会通过滞后特征影响到验证集第一年的标签,gap 能切断这条泄漏路径。scale_pos_weight的作用是把干旱样本的损失权重调高,因为干旱事件本来就比较少,如果不平衡设置,模型会倾向把所有样本都预测为“不干旱”,整体准确率很高,但对预警业务毫无用处。

4.2 评估指标:CSI、POD、FAR 才是业务关心的数字

分类准确率在干旱预警里没有意义。假设干旱事件只有 15% 的频率,模型什么都不做全预测“不干旱”,准确率也是 85%,但决策者拿到这个模型等于没有模型。气象预警领域通行的指标是 CSI(Critical Success Index)、POD(Probability of Detection)和 FAR(False 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, threshold=0.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, label="Observed SPI-3", color="black") plt.plot(timestamps_valid, pred_prob, label="Predicted drought prob", color="red") plt.axhline(-1.0, linestyle="--", color="gray") plt.legend() plt.show()

曲线贴合不代表好,如果预测概率曲线看起来像 SPI 曲线的左右平移,说明模型只是在延迟复述观测。好的预测是概率峰值出现在 SPI 谷值之前,这才是“预测”而不是“跟唱”。

5. 四个必踩的坑:数据泄漏、样本不平衡与气候外推

5.1 坑一:随机切分导致“用未来预测过去”,验证集虚高

现象:用默认train_test_split训练模型,验证集 AUC 高达 0.85,CSI 也很漂亮。换成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 年用最新数据重新拟合一次模型,给旧样本设衰减权重;三是监控特征分布漂移,比如在周度、月度的模型评估中加上 PSI(Population 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_names=feature_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 时间序列和特征分布图,而不是急着跑模型——数据要是脏的,模型再花哨也是白搭。希望这些从数据到评估再到解释的流程和踩坑记录能帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询