简介这是一份以机器学习应对干旱预测的开源研究型代码包面向气候科学、遥感或数据科学领域有一定Python基础的研究者与学生。项目将不同机器学习方法组织成可统一训练与测试的端到端管道从原始数据格式整理、特征与样本构建到模型训练、相互比较与误差评估均有模块化实现代码按功能拆分为src下的不同类方便替换或新增算法复现实验。包体为zip压缩格式大小约49.31MB文件总数显示为0平台未提供文件类型明细结合项目描述和JupyterNotebook标签核心内容预计以Python脚本、Notebook演示及conda环境配套YAML为主。资源还提供三个程序入口并附带设计思路和配置说明可通过environment.yml快速创建名为esowc-drought的独立环境。已有142人学习查看对希望把机器学习引入气候预测的入门到进阶学习者而言是一份可直接运行、便于二次开发的完整范例。1. 机器学习预测干旱ml_drought 项目想解决的核心问题ml_drought 是 GitHub 上 ml-clim 组织下的一个项目标题直译就是“用机器学习更好地预测和了解干旱”。第一次刷到这个项目时我没急着看模型架构而是先翻了它的数据目录和特征定义——因为做过几个干旱预测项目之后我越来越确信这个方向真正的门槛不在模型深度而在于怎么把降水、温度、土壤湿度这些时序观测变成一份不会泄漏未来信息的训练样本。机器学习能补足传统统计和物理模型的两个短板一是从高维气象变量里自动找非线性组合二是给出因子贡献度让预测不止于概率数字。适合读这篇文章的人气象、农业、水文领域想引入机器学习的工程师做灾害预警的学生以及想评估 GitHub 上这类开源项目值不值得深挖的开发者。2. 先解决“用什么数据”构造干旱预测数据集的四个关键选择2.1 干旱指标选型SPI、SPEI 到底用哪个做干旱预测第一件事是定义“旱”。气象观测里最常用的标准化降水指数SPI只吃降水按每个地区的历史分布做标准化所以对不同气候区天然可比。SPI-3 代表近 3 个月降水偏少程度适合反映季节尺度的气象干旱SPI-6 或 SPI-12 则偏向水资源和长期干旱。如果关心农业或生态就需要考虑温度高的地方蒸散发也大这时候用标准化降水蒸散发指数SPEI更合理。SPEI 在 SPI 的算法基础上加入了气温驱动的潜在蒸散发项能捕捉升温对干旱的放大效应。在 ml_drought 项目里如果目标是“了解干旱”我建议同时算两个指标一个做预测目标一个做特征或对照。这样你还能看出模型到底在靠哪套指标学习。选择上有一个简单的经验只有降水资料用 SPI有日最高/最低温度和降水优先 SPEI业务系统面向农业保险和作物生长季直接做 SPEI-3 或 SPEI-6。注意无论选哪个窗口长度对标签的影响远大于模型选择。2.2 特征目录降水、温度、植被和遥相关因子类别常用特征时间尺度说明气象降水距平、月均温、最高/最低温度、相对湿度、风速月距平需要减去当地历史均值水文土壤湿度、蒸散发、径流月可用再分析产品或陆面模式输出植被NDVI、EVI、VCI旬/月干旱发生后期植被会响应遥相关Nino3.4、SOI、IOD、AO月大尺度气候模式对区域干旱有前兆性静态经度、纬度、海拔、气候区静态帮助模型区分区域差异这些特征不是全都要收。刚开始做先用降水距平、温度距平、土壤湿度和两三个遥相关指数即可。特征太多会让模型在样本量少时学到站点特有规律迁到新地区就翻车。NDVI 这种响应指标虽然对“了解干旱”有用但在做未来预测时要格外小心如果目标是提前 3 个月预警而植被往往滞后一个月才出现枯黄那它在预测任务里就不是前兆信号。我的习惯是先建立一套最小特征集跑出基线再用 SHAP 决定加什么。2.3 自己写一个 SPI-3 计算函数别用黑匣子import numpy as np import pandas as pd from scipy.stats import gamma, norm def calc_spi(precip, scale3, min_nonzero30): 用滑动累积降水计算 SPI。 precip 必须按月排序单位 mm缺失值先填充。 scale: 累积月数3 表示 SPI-3。 min_nonzero: 非零样本最小数量。 # 1. 滑动累积形成 3 个月累计降水序列 cum pd.Series(precip).rolling(scale, min_periods1).sum().dropna() # 2. 用非零累积值拟合 Gamma 分布 nonzero cum[cum 0] if len(nonzero) min_nonzero: raise ValueError(f非零样本 {len(nonzero)} 太少无法稳定拟合) shape, loc, scale_param gamma.fit(nonzero.values, floc0) # 3. 计算每个值的累积概率 cdf gamma.cdf(cum.values, shape, locloc, scalescale_param) # 4. 零值样本单独处理赋予小于所有正样本的概率 zero_mask cum.values 0 if zero_mask.any(): cdf[zero_mask] np.linspace(0, 0.001, zero_mask.sum()) # 5. 转成标准正态分位数即 SPI cdf np.clip(cdf, 1e-6, 1 - 1e-6) return pd.Series(norm.ppf(cdf), indexcum.index)第一步rolling(scale)生成 3 个月累积降水SPI 的本质是对不同尺度累积降水做概率标准化。第二步用 Gamma 拟合是因为月降水累积量通常服从偏态分布而 Gamma 族能描述这种偏态。第三步得到每个累计值的百分位。floc0固定位置参数为 0让分布从 0 开始否则拟合可能产生负值支撑不符合降水特征。第四步处理零降水如果某地常年有零降水月Gamma 拟合对零值不友好这里给它们一个接近 0 的累积概率避免算出负无穷 SPI。最后用norm.ppf把累积概率映射到标准正态SPI 的负数代表比常年偏干正数代表偏湿。提示SPI 公式看起来简单但窗口、概率分布、零值处理都会影响结果。你在 GitHub 项目评估时如果看到代码里没有对零值做处理就要怀疑标签质量。2.4 时间对齐与滞后特征避免“用未来预测过去”特征矩阵里每一行对应位置目标月份。假设我们在做 T 月预测 T3 月是否干旱那么 X 只能使用截面时间 T 及之前的信息。真实项目中有人会把 T 到 T3 的降水距平都放进特征结果 AUC 漂亮得不可思议。另外如果标签是 SPI-3而特征里也放 T 月的 3 个月累积降水那这个特征和标签几乎同源模型不再预测只是在复述当前状态。def build_samples(features, target, lead3, feature_lag1): features: 按 (时间, 地点) 展开的特征表含 year, month target: T 月的干旱标签序列 lead: 提前几个月预测 feature_lag: 特征截止到目标月之前的第几个月 X, y [], [] for t in range(max(lead, feature_lag) 1, len(target)): # 特征截止于 t-feature_lag 月预测 tlead 月的干旱状态 X.append(features.iloc[t - feature_lag].values) y.append(target[t lead]) return np.array(X), np.array(y)这段代码是逻辑示例实际还会带上地区和季节列。核心思路是特征时间窗和目标时间窗要有明确间隔。lead越大任务越难但业务价值越高因为预警提前量更大。3. 模型选型与训练从随机森林到 LSTM 的三步走3.1 先跑随机森林把基线钉死数据到位后不要急着上深度学习。干旱事件样本量通常小一个站点十年数据才 120 个月正样本可能只有几十条。随机森林对缺失值不敏感能捕捉非线性交互训练速度也快非常适合做基线。如果随机森林做不出合理效果问题八成在特征或标签不在模型。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split, cross_val_score X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, shuffleFalse, stratifyNone ) model RandomForestClassifier( n_estimators500, max_depth6, min_samples_leaf10, max_featuressqrt, class_weightbalanced, random_state42, n_jobs-1 ) scores cross_val_score(model, X_train, y_train, cv5, scoringroc_auc) print(f交叉验证 AUC: {scores.mean():.3f} ± {scores.std():.3f})n_estimators500是树的数量再高收益递减max_depth6限制单棵树深度防止过拟合min_samples_leaf10让叶节点至少 10 个样本不让模型记住个别站点max_featuressqrt随机抽取特征子集增加树之间多样性class_weightbalanced自动按类别比例加权缓解干旱样本少的问题。这里shuffleFalse是因为时序数据必须保持顺序。随机森林默认的 KFold 会打乱时间顺序所以这个得分偏乐观后面第 6 章会专门讲时序交叉验证。3.2 XGBoost 调参类别不平衡和缺失值处理随机森林跑通之后XGBoost 是更好的第二选择。它对特征单调性约束更灵活自带缺失值学习路径而且能通过scale_pos_weight直接控制不平衡比率。干旱预测里模型往往要面对“几十年才有几次严重干旱”的数据这个参数比准确率重要得多。import xgboost as xgb neg sum(y_train 0) pos sum(y_train 1) model_xgb xgb.XGBClassifier( n_estimators300, max_depth4, learning_rate0.05, subsample0.8, colsample_bytree0.8, scale_pos_weightneg / pos, reg_lambda1.0, eval_metricauc, early_stopping_rounds20, use_label_encoderFalse ) model_xgb.fit( X_train, y_train, eval_set[(X_test, y_test)], verboseFalse )max_depth4是因为 XGBoost 深度太大容易让叶子权重过拟合干旱特征交互通常不深。subsample0.8和colsample_bytree0.8分别控制行采样和列采样增加鲁棒性。scale_pos_weightneg / pos把少数类别权重放大到多数类别与少数类别的比值实际效果通常比class_weight更直接。reg_lambda1.0是 L2 正则防止高维特征下权重爆炸。early_stopping_rounds20在验证集上连续 20 轮不涨就停。如果你的 XGBoost 版本低于 1.6use_label_encoderFalse可以删掉。3.3 LSTM 做序列预测滑窗、归一化和早停树模型看不到“连续几个月降水偏低”这样的动态模式LSTM 可以这也是标题里“机器学习预测干旱”最常被想到的路线。但 LSTM 的坑也很明显样本量少、归一化不当、滑窗随机抽样都会导致翻车。我一般只在树模型已经稳定之后才试 LSTM。import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout from tensorflow.keras.callbacks import EarlyStopping # 假设已经按时间顺序归一化过特征 def make_sequences(X, y, timesteps12): Xs, ys [], [] for i in range(len(X) - timesteps): Xs.append(X[i:i timesteps]) ys.append(y[i timesteps]) return np.array(Xs), np.array(ys) X_seq, y_seq make_sequences(X_scaled, y, timesteps12) model Sequential([ LSTM(32, return_sequencesTrue, input_shape(X_seq.shape[1], X_seq.shape[2])), Dropout(0.2), LSTM(16), Dense(1, activationsigmoid) ]) model.compile(optimizeradam, lossbinary_crossentropy, metrics[AUC]) early EarlyStopping(monitorval_loss, patience5, restore_best_weightsTrue) model.fit(X_seq, y_seq, epochs50, batch_size64, validation_split0.2, callbacks[early], verbose0)timesteps12意味着用过去 12 个月的特征预测未来包含完整年周期。第一层LSTM(32, return_sequencesTrue)输出序列给第二层LSTM(16)继续提取时序特征。Dropout(0.2)防过拟合。validation_split0.2直接切最后 20%但这里有个隐患如果数据不是按时间顺序排列验证集会穿越时间。实际执行时建议先把数据排序再手动切分不要依赖默认validation_split。EarlyStopping的patience5对干旱任务常用 5~10太小会欠拟合。如果你的样本量少于几百LSTM 很可能打不过 XGBoost这不是能力问题是数据量问题。3.4 评估指标选 CSI 和 HSS而不是准确率指标公式关注点CSIhits / (hits misses false_alarms)干旱事件捕获质量忽略正确无旱样本HSS(hits - expected_hits) / (N - expected_hits)相对随机分类的提升对不平衡敏感AUC排序能力概率排序但不反映概率校准Accuracy(hits correct_no_rain) / N不平衡时无意义干旱预测里“无旱”是绝大多数样本全预测无旱准确率可能超过 95%但业务完全不可用。CSI 只看该报警的样本有多少被报出来HSS 则衡量模型是否比瞎猜强。比如验证集上模型漏掉一半干旱事件AUC 可能还有 0.82但 CSI 只有 0.35。所以评估一个 GitHub 上的干旱预测项目先看验证集报的是 CSI/HSS 还是准确率就能判断是否经得起业务考验。4. 让模型“说出”干旱的原因SHAP 归因的实操4.1 特征重要性为什么不够用随机森林和 XGBoost 都可以打印feature_importances_但那只是一个总分某个特征在所有树分裂中带来的纯度提升和频率。它看不到特征与干旱概率的方向性也不回答“什么条件下这个特征才起作用”。比如“3 个月降水距平”很重要模型预测干旱概率升高但可能只是在某些站点、某些季节里才有效。SHAP 给每个样本的每个特征一个贡献值能逐条解释模型为什么判断某个月份会旱这正是标题里“了解干旱”的含义。4.2 用 SHAP 在随机森林上做全局解释import shap explainer shap.TreeExplainer(model) # model 为已训练的 XGBoost 或 RandomForest shap_values explainer.shap_values(X_test) shap.summary_plot(shap_values, X_test, feature_namesfeature_names, max_display15)TreeExplainer对树模型有精确算法比KernelExplainer快得多。summary_plot把样本按 SHAP 值排序横轴表示该特征对预测的贡献方向颜色代表原始特征值高低。如果特征多设置max_display15只看前 15 个。随机森林的shap_values在部分版本中可能是列表形式报错时查看 shap 版本 API通常建议用shap.TreeExplainer(model).shap_values(X)。4.3 分季节和分区域看关键驱动因子干旱驱动机制有明显的季节和区域差异。把测试集按月份拆成 DJF、MAM、JJA、SON 四组分别画 SHAP 图经常能看到相反的符号。例如季风区夏季干旱与“前 3 个月降水距平”负相关明显冬季则与“温度距平”正相关高温加速积雪融化。如果不分组这些细节会被平均掉。import matplotlib.pyplot as plt for season in [DJF, MAM, JJA, SON]: mask (test_month // 3) % 4 {DJF: 0, MAM: 1, JJA: 2, SON: 3}[season] sv explainer.shap_values(X_test[mask]) plt.figure() shap.summary_plot(sv, X_test[mask], feature_namesfeature_names, showFalse) plt.title(f{season} SHAP) plt.show()test_month // 3 % 4只是示例映射你需要按实际月份列调整。分季节后样本量减少解释结果可能波动大要看多个年份是否一致不要只凭一季判断。如果某个季节只有十几次干旱样本SHAP 图画出来的排序很可能包含噪声。4.4 把归因结果写进干旱成因分析报告SHAP 不只给模型找补它可以产出业务报告。比如从 SHAP 依赖图上看到“当 3 个月降水距平低于 -30mm 且温度距平高于 1.5℃ 时干旱概率超过 0.7”你可以在预警产品里给出类似“高风险降水严重偏少叠加高温”的成因标签而不只输出概率。这是机器学习从黑匣子走向决策辅助的关键一步。具体做法是基于 SHAP 值把每个预测样本归因到一个主要驱动因子生成“时间、地点、概率、主导因子”的表格交给业务侧去核对。5. 干旱预测项目常见的 5 个坑现象、原因与解决5.1 标签泄漏训练时 AUC 高得反常现象训练集 AUC 0.999测试集 AUC 0.95上线后预报成功率却低得离谱。原因最常见的泄漏有两个一是把目标月份之后的降水、温度数据当作特征二是在特征计算时用了全样本的均值与方差做标准化。比如StandardScaler.fit在整份数据上调用而不是只 fit 训练集。解决先按时间顺序划分训练、验证、测试再做特征工程和归一化。检查那些“预测超前”的气候因子比如 ENSO 指数的发布时间是否真的早于预测时刻。如果你发现模型输入里有时刻 T 之后的数据整个工作流的结论都要推倒重来。5.2 时间错位用同期降水预测同期干旱现象模型觉得自己找到了规律当月降水距平一低当月干旱标签就为 1。原因SPI-3 标签本身是过去 3 个月降水累计的结果如果用当月降水距平当特征实际上是拿结果解释结果。解决要么把特征滞后至少一个完整窗口比如预测 SPI-3 时特征只到 T-3 月要么改用目标月前的气候态距平比如历史同期均值距平并删除所有和标签窗口重叠的累积变量。这个坑在时间序列机器学习项目里非常普遍GitHub 仓库的 README 一般不会写得靠代码审查发现。5.3 空间自相关随机切分把验证集变成“同地域复制”现象用 KFold 随机切分时 CSI 有 0.6按站点分组切分后掉到 0.35。原因相邻站点的气象数据高度相关随机切分会把同一片气候区的站点同时塞进训练和验证模型相当于见过“邻居答案”。解决用GroupKFold把 group 设为站点 ID 或网格 ID如果担心时间因素用GroupShuffleSplit按站点分组。做迁移评估时还可以按气候区划分训练和验证这样更接近真实部署。5.4 类别不平衡模型永远预报“无旱”现象看分类报告无旱的精确率和召回率都很高干旱的召回率接近 0但 AUC 还挺高。原因AUC 只看正样本排在负样本前面的概率即使所有样本输出概率都在 0.05 以下只要正样本稍微高于负样本AUC 也能到 0.8。业务需要的是概率可比较的预警AUC 不能反映这一点。解决训练时用scale_pos_weight或class_weight预测后做阈值校准评估用 CSI 和 HSS。另外严重干旱样本极少可以考虑把标签放宽为“轻度干旱及以上”增加正样本数量。但这样会把目标从灾害预警变成状态监测业务上要讲清楚。5.5 迁移失效换个气候区就崩现象在甲地区训练好的模型直接拿到乙地区预测CSI 从 0.6 降到 0.2。原因模型学到了地区特有的气候条件。湿润区干旱主要由连续无雨日数决定半干旱区则更依赖降水距平百分比如果特征里没有气候区或标准化方式不一致模型无法外推。解决把气候区作为特征或者对每个气候区做单独的标准化。更稳健的做法是加入静态地理特征纬度、海拔等让模型按区域自适应。如果有多个气候区的数据用区域分组交叉验证评估迁移能力。很多开放数据集上效果好的模型最后都是在这一步被刷下来的。6. 让模型真正可用时序交叉验证、决策阈值校准与重训练节奏6.1 用 TimeSeriesSplit 替代 KFoldfrom sklearn.model_selection import TimeSeriesSplit from sklearn.metrics import roc_auc_score tscv TimeSeriesSplit(n_splits5, gap3) fold_scores [] for train_index, test_index in tscv.split(X): X_train_fold, X_test_fold X[train_index], X[test_index] y_train_fold, y_test_fold y[train_index], y[test_index] model_fold RandomForestClassifier( n_estimators300, max_depth6, min_samples_leaf10, class_weightbalanced, random_state42, n_jobs-1 ) model_fold.fit(X_train_fold, y_train_fold) fold_scores.append(roc_auc_score(y_test_fold, model_fold.predict_proba(X_test_fold)[:, 1])) print(f时序交叉验证 AUC: {np.mean(fold_scores):.3f} ± {np.std(fold_scores):.3f})gap3让训练集最后 3 个月不参与训练防止自相关导致信息穿越。这套打分比普通 KFold 更接近真实上线场景。做干旱预测如果时序交叉验证和随机切分的结果相差很大先别急着调模型回去检查特征是否泄漏。6.2 根据代价矩阵校准预警阈值干旱预警中漏报代价远高于误报所以不能把 0.5 作为决策边界。用验证集概率画 Precision-Recall 曲线挑选 F1 最高或业务代价最小的阈值。from sklearn.metrics import precision_recall_curve prec, rec, thr precision_recall_curve(y_val, proba_val) f1_scores 2 * prec * rec / (prec rec 1e-6) best_thr thr[np.argmax(f1_scores)]如果业务里漏报代价是误报的 3 倍可以按代价定义cost misses * cost_miss false_alarms * cost_fa再找到最低 cost 的阈值。这个阈值要每个月根据最新验证集重新校准。6.3 重训练与数据漂移监控气候非平稳模型会随时间漂移。我一般每季度或半年重训练一次具体频率看业务更新周期。同时监控特征分布变化用 PSI 这类指标。当某个特征分布漂移超过阈值比如降水距平均值变化 0.5 个标准差需要把新数据并入训练集并重训。这个环节在开源 GitHub 项目里通常没有但实际落地必须自己补。我现在做干旱预测项目习惯先把数据构造代码和验证方式当作第一优先级模型清单排第二。很多开源项目模型很漂亮但数据泄漏和随机切分的问题一眼就能看到。把这个思路固定下来后我自己的预测系统上线翻车的次数少了很多。希望帮到你。本文还有配套的精品资源点击获取
机器学习预测干旱全指南：从特征工程到模型落地
简介这是一份以机器学习应对干旱预测的开源研究型代码包面向气候科学、遥感或数据科学领域有一定Python基础的研究者与学生。项目将不同机器学习方法组织成可统一训练与测试的端到端管道从原始数据格式整理、特征与样本构建到模型训练、相互比较与误差评估均有模块化实现代码按功能拆分为src下的不同类方便替换或新增算法复现实验。包体为zip压缩格式大小约49.31MB文件总数显示为0平台未提供文件类型明细结合项目描述和JupyterNotebook标签核心内容预计以Python脚本、Notebook演示及conda环境配套YAML为主。资源还提供三个程序入口并附带设计思路和配置说明可通过environment.yml快速创建名为esowc-drought的独立环境。已有142人学习查看对希望把机器学习引入气候预测的入门到进阶学习者而言是一份可直接运行、便于二次开发的完整范例。1. 机器学习预测干旱ml_drought 项目想解决的核心问题ml_drought 是 GitHub 上 ml-clim 组织下的一个项目标题直译就是“用机器学习更好地预测和了解干旱”。第一次刷到这个项目时我没急着看模型架构而是先翻了它的数据目录和特征定义——因为做过几个干旱预测项目之后我越来越确信这个方向真正的门槛不在模型深度而在于怎么把降水、温度、土壤湿度这些时序观测变成一份不会泄漏未来信息的训练样本。机器学习能补足传统统计和物理模型的两个短板一是从高维气象变量里自动找非线性组合二是给出因子贡献度让预测不止于概率数字。适合读这篇文章的人气象、农业、水文领域想引入机器学习的工程师做灾害预警的学生以及想评估 GitHub 上这类开源项目值不值得深挖的开发者。2. 先解决“用什么数据”构造干旱预测数据集的四个关键选择2.1 干旱指标选型SPI、SPEI 到底用哪个做干旱预测第一件事是定义“旱”。气象观测里最常用的标准化降水指数SPI只吃降水按每个地区的历史分布做标准化所以对不同气候区天然可比。SPI-3 代表近 3 个月降水偏少程度适合反映季节尺度的气象干旱SPI-6 或 SPI-12 则偏向水资源和长期干旱。如果关心农业或生态就需要考虑温度高的地方蒸散发也大这时候用标准化降水蒸散发指数SPEI更合理。SPEI 在 SPI 的算法基础上加入了气温驱动的潜在蒸散发项能捕捉升温对干旱的放大效应。在 ml_drought 项目里如果目标是“了解干旱”我建议同时算两个指标一个做预测目标一个做特征或对照。这样你还能看出模型到底在靠哪套指标学习。选择上有一个简单的经验只有降水资料用 SPI有日最高/最低温度和降水优先 SPEI业务系统面向农业保险和作物生长季直接做 SPEI-3 或 SPEI-6。注意无论选哪个窗口长度对标签的影响远大于模型选择。2.2 特征目录降水、温度、植被和遥相关因子类别常用特征时间尺度说明气象降水距平、月均温、最高/最低温度、相对湿度、风速月距平需要减去当地历史均值水文土壤湿度、蒸散发、径流月可用再分析产品或陆面模式输出植被NDVI、EVI、VCI旬/月干旱发生后期植被会响应遥相关Nino3.4、SOI、IOD、AO月大尺度气候模式对区域干旱有前兆性静态经度、纬度、海拔、气候区静态帮助模型区分区域差异这些特征不是全都要收。刚开始做先用降水距平、温度距平、土壤湿度和两三个遥相关指数即可。特征太多会让模型在样本量少时学到站点特有规律迁到新地区就翻车。NDVI 这种响应指标虽然对“了解干旱”有用但在做未来预测时要格外小心如果目标是提前 3 个月预警而植被往往滞后一个月才出现枯黄那它在预测任务里就不是前兆信号。我的习惯是先建立一套最小特征集跑出基线再用 SHAP 决定加什么。2.3 自己写一个 SPI-3 计算函数别用黑匣子import numpy as np import pandas as pd from scipy.stats import gamma, norm def calc_spi(precip, scale3, min_nonzero30): 用滑动累积降水计算 SPI。 precip 必须按月排序单位 mm缺失值先填充。 scale: 累积月数3 表示 SPI-3。 min_nonzero: 非零样本最小数量。 # 1. 滑动累积形成 3 个月累计降水序列 cum pd.Series(precip).rolling(scale, min_periods1).sum().dropna() # 2. 用非零累积值拟合 Gamma 分布 nonzero cum[cum 0] if len(nonzero) min_nonzero: raise ValueError(f非零样本 {len(nonzero)} 太少无法稳定拟合) shape, loc, scale_param gamma.fit(nonzero.values, floc0) # 3. 计算每个值的累积概率 cdf gamma.cdf(cum.values, shape, locloc, scalescale_param) # 4. 零值样本单独处理赋予小于所有正样本的概率 zero_mask cum.values 0 if zero_mask.any(): cdf[zero_mask] np.linspace(0, 0.001, zero_mask.sum()) # 5. 转成标准正态分位数即 SPI cdf np.clip(cdf, 1e-6, 1 - 1e-6) return pd.Series(norm.ppf(cdf), indexcum.index)第一步rolling(scale)生成 3 个月累积降水SPI 的本质是对不同尺度累积降水做概率标准化。第二步用 Gamma 拟合是因为月降水累积量通常服从偏态分布而 Gamma 族能描述这种偏态。第三步得到每个累计值的百分位。floc0固定位置参数为 0让分布从 0 开始否则拟合可能产生负值支撑不符合降水特征。第四步处理零降水如果某地常年有零降水月Gamma 拟合对零值不友好这里给它们一个接近 0 的累积概率避免算出负无穷 SPI。最后用norm.ppf把累积概率映射到标准正态SPI 的负数代表比常年偏干正数代表偏湿。提示SPI 公式看起来简单但窗口、概率分布、零值处理都会影响结果。你在 GitHub 项目评估时如果看到代码里没有对零值做处理就要怀疑标签质量。2.4 时间对齐与滞后特征避免“用未来预测过去”特征矩阵里每一行对应位置目标月份。假设我们在做 T 月预测 T3 月是否干旱那么 X 只能使用截面时间 T 及之前的信息。真实项目中有人会把 T 到 T3 的降水距平都放进特征结果 AUC 漂亮得不可思议。另外如果标签是 SPI-3而特征里也放 T 月的 3 个月累积降水那这个特征和标签几乎同源模型不再预测只是在复述当前状态。def build_samples(features, target, lead3, feature_lag1): features: 按 (时间, 地点) 展开的特征表含 year, month target: T 月的干旱标签序列 lead: 提前几个月预测 feature_lag: 特征截止到目标月之前的第几个月 X, y [], [] for t in range(max(lead, feature_lag) 1, len(target)): # 特征截止于 t-feature_lag 月预测 tlead 月的干旱状态 X.append(features.iloc[t - feature_lag].values) y.append(target[t lead]) return np.array(X), np.array(y)这段代码是逻辑示例实际还会带上地区和季节列。核心思路是特征时间窗和目标时间窗要有明确间隔。lead越大任务越难但业务价值越高因为预警提前量更大。3. 模型选型与训练从随机森林到 LSTM 的三步走3.1 先跑随机森林把基线钉死数据到位后不要急着上深度学习。干旱事件样本量通常小一个站点十年数据才 120 个月正样本可能只有几十条。随机森林对缺失值不敏感能捕捉非线性交互训练速度也快非常适合做基线。如果随机森林做不出合理效果问题八成在特征或标签不在模型。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split, cross_val_score X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, shuffleFalse, stratifyNone ) model RandomForestClassifier( n_estimators500, max_depth6, min_samples_leaf10, max_featuressqrt, class_weightbalanced, random_state42, n_jobs-1 ) scores cross_val_score(model, X_train, y_train, cv5, scoringroc_auc) print(f交叉验证 AUC: {scores.mean():.3f} ± {scores.std():.3f})n_estimators500是树的数量再高收益递减max_depth6限制单棵树深度防止过拟合min_samples_leaf10让叶节点至少 10 个样本不让模型记住个别站点max_featuressqrt随机抽取特征子集增加树之间多样性class_weightbalanced自动按类别比例加权缓解干旱样本少的问题。这里shuffleFalse是因为时序数据必须保持顺序。随机森林默认的 KFold 会打乱时间顺序所以这个得分偏乐观后面第 6 章会专门讲时序交叉验证。3.2 XGBoost 调参类别不平衡和缺失值处理随机森林跑通之后XGBoost 是更好的第二选择。它对特征单调性约束更灵活自带缺失值学习路径而且能通过scale_pos_weight直接控制不平衡比率。干旱预测里模型往往要面对“几十年才有几次严重干旱”的数据这个参数比准确率重要得多。import xgboost as xgb neg sum(y_train 0) pos sum(y_train 1) model_xgb xgb.XGBClassifier( n_estimators300, max_depth4, learning_rate0.05, subsample0.8, colsample_bytree0.8, scale_pos_weightneg / pos, reg_lambda1.0, eval_metricauc, early_stopping_rounds20, use_label_encoderFalse ) model_xgb.fit( X_train, y_train, eval_set[(X_test, y_test)], verboseFalse )max_depth4是因为 XGBoost 深度太大容易让叶子权重过拟合干旱特征交互通常不深。subsample0.8和colsample_bytree0.8分别控制行采样和列采样增加鲁棒性。scale_pos_weightneg / pos把少数类别权重放大到多数类别与少数类别的比值实际效果通常比class_weight更直接。reg_lambda1.0是 L2 正则防止高维特征下权重爆炸。early_stopping_rounds20在验证集上连续 20 轮不涨就停。如果你的 XGBoost 版本低于 1.6use_label_encoderFalse可以删掉。3.3 LSTM 做序列预测滑窗、归一化和早停树模型看不到“连续几个月降水偏低”这样的动态模式LSTM 可以这也是标题里“机器学习预测干旱”最常被想到的路线。但 LSTM 的坑也很明显样本量少、归一化不当、滑窗随机抽样都会导致翻车。我一般只在树模型已经稳定之后才试 LSTM。import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout from tensorflow.keras.callbacks import EarlyStopping # 假设已经按时间顺序归一化过特征 def make_sequences(X, y, timesteps12): Xs, ys [], [] for i in range(len(X) - timesteps): Xs.append(X[i:i timesteps]) ys.append(y[i timesteps]) return np.array(Xs), np.array(ys) X_seq, y_seq make_sequences(X_scaled, y, timesteps12) model Sequential([ LSTM(32, return_sequencesTrue, input_shape(X_seq.shape[1], X_seq.shape[2])), Dropout(0.2), LSTM(16), Dense(1, activationsigmoid) ]) model.compile(optimizeradam, lossbinary_crossentropy, metrics[AUC]) early EarlyStopping(monitorval_loss, patience5, restore_best_weightsTrue) model.fit(X_seq, y_seq, epochs50, batch_size64, validation_split0.2, callbacks[early], verbose0)timesteps12意味着用过去 12 个月的特征预测未来包含完整年周期。第一层LSTM(32, return_sequencesTrue)输出序列给第二层LSTM(16)继续提取时序特征。Dropout(0.2)防过拟合。validation_split0.2直接切最后 20%但这里有个隐患如果数据不是按时间顺序排列验证集会穿越时间。实际执行时建议先把数据排序再手动切分不要依赖默认validation_split。EarlyStopping的patience5对干旱任务常用 5~10太小会欠拟合。如果你的样本量少于几百LSTM 很可能打不过 XGBoost这不是能力问题是数据量问题。3.4 评估指标选 CSI 和 HSS而不是准确率指标公式关注点CSIhits / (hits misses false_alarms)干旱事件捕获质量忽略正确无旱样本HSS(hits - expected_hits) / (N - expected_hits)相对随机分类的提升对不平衡敏感AUC排序能力概率排序但不反映概率校准Accuracy(hits correct_no_rain) / N不平衡时无意义干旱预测里“无旱”是绝大多数样本全预测无旱准确率可能超过 95%但业务完全不可用。CSI 只看该报警的样本有多少被报出来HSS 则衡量模型是否比瞎猜强。比如验证集上模型漏掉一半干旱事件AUC 可能还有 0.82但 CSI 只有 0.35。所以评估一个 GitHub 上的干旱预测项目先看验证集报的是 CSI/HSS 还是准确率就能判断是否经得起业务考验。4. 让模型“说出”干旱的原因SHAP 归因的实操4.1 特征重要性为什么不够用随机森林和 XGBoost 都可以打印feature_importances_但那只是一个总分某个特征在所有树分裂中带来的纯度提升和频率。它看不到特征与干旱概率的方向性也不回答“什么条件下这个特征才起作用”。比如“3 个月降水距平”很重要模型预测干旱概率升高但可能只是在某些站点、某些季节里才有效。SHAP 给每个样本的每个特征一个贡献值能逐条解释模型为什么判断某个月份会旱这正是标题里“了解干旱”的含义。4.2 用 SHAP 在随机森林上做全局解释import shap explainer shap.TreeExplainer(model) # model 为已训练的 XGBoost 或 RandomForest shap_values explainer.shap_values(X_test) shap.summary_plot(shap_values, X_test, feature_namesfeature_names, max_display15)TreeExplainer对树模型有精确算法比KernelExplainer快得多。summary_plot把样本按 SHAP 值排序横轴表示该特征对预测的贡献方向颜色代表原始特征值高低。如果特征多设置max_display15只看前 15 个。随机森林的shap_values在部分版本中可能是列表形式报错时查看 shap 版本 API通常建议用shap.TreeExplainer(model).shap_values(X)。4.3 分季节和分区域看关键驱动因子干旱驱动机制有明显的季节和区域差异。把测试集按月份拆成 DJF、MAM、JJA、SON 四组分别画 SHAP 图经常能看到相反的符号。例如季风区夏季干旱与“前 3 个月降水距平”负相关明显冬季则与“温度距平”正相关高温加速积雪融化。如果不分组这些细节会被平均掉。import matplotlib.pyplot as plt for season in [DJF, MAM, JJA, SON]: mask (test_month // 3) % 4 {DJF: 0, MAM: 1, JJA: 2, SON: 3}[season] sv explainer.shap_values(X_test[mask]) plt.figure() shap.summary_plot(sv, X_test[mask], feature_namesfeature_names, showFalse) plt.title(f{season} SHAP) plt.show()test_month // 3 % 4只是示例映射你需要按实际月份列调整。分季节后样本量减少解释结果可能波动大要看多个年份是否一致不要只凭一季判断。如果某个季节只有十几次干旱样本SHAP 图画出来的排序很可能包含噪声。4.4 把归因结果写进干旱成因分析报告SHAP 不只给模型找补它可以产出业务报告。比如从 SHAP 依赖图上看到“当 3 个月降水距平低于 -30mm 且温度距平高于 1.5℃ 时干旱概率超过 0.7”你可以在预警产品里给出类似“高风险降水严重偏少叠加高温”的成因标签而不只输出概率。这是机器学习从黑匣子走向决策辅助的关键一步。具体做法是基于 SHAP 值把每个预测样本归因到一个主要驱动因子生成“时间、地点、概率、主导因子”的表格交给业务侧去核对。5. 干旱预测项目常见的 5 个坑现象、原因与解决5.1 标签泄漏训练时 AUC 高得反常现象训练集 AUC 0.999测试集 AUC 0.95上线后预报成功率却低得离谱。原因最常见的泄漏有两个一是把目标月份之后的降水、温度数据当作特征二是在特征计算时用了全样本的均值与方差做标准化。比如StandardScaler.fit在整份数据上调用而不是只 fit 训练集。解决先按时间顺序划分训练、验证、测试再做特征工程和归一化。检查那些“预测超前”的气候因子比如 ENSO 指数的发布时间是否真的早于预测时刻。如果你发现模型输入里有时刻 T 之后的数据整个工作流的结论都要推倒重来。5.2 时间错位用同期降水预测同期干旱现象模型觉得自己找到了规律当月降水距平一低当月干旱标签就为 1。原因SPI-3 标签本身是过去 3 个月降水累计的结果如果用当月降水距平当特征实际上是拿结果解释结果。解决要么把特征滞后至少一个完整窗口比如预测 SPI-3 时特征只到 T-3 月要么改用目标月前的气候态距平比如历史同期均值距平并删除所有和标签窗口重叠的累积变量。这个坑在时间序列机器学习项目里非常普遍GitHub 仓库的 README 一般不会写得靠代码审查发现。5.3 空间自相关随机切分把验证集变成“同地域复制”现象用 KFold 随机切分时 CSI 有 0.6按站点分组切分后掉到 0.35。原因相邻站点的气象数据高度相关随机切分会把同一片气候区的站点同时塞进训练和验证模型相当于见过“邻居答案”。解决用GroupKFold把 group 设为站点 ID 或网格 ID如果担心时间因素用GroupShuffleSplit按站点分组。做迁移评估时还可以按气候区划分训练和验证这样更接近真实部署。5.4 类别不平衡模型永远预报“无旱”现象看分类报告无旱的精确率和召回率都很高干旱的召回率接近 0但 AUC 还挺高。原因AUC 只看正样本排在负样本前面的概率即使所有样本输出概率都在 0.05 以下只要正样本稍微高于负样本AUC 也能到 0.8。业务需要的是概率可比较的预警AUC 不能反映这一点。解决训练时用scale_pos_weight或class_weight预测后做阈值校准评估用 CSI 和 HSS。另外严重干旱样本极少可以考虑把标签放宽为“轻度干旱及以上”增加正样本数量。但这样会把目标从灾害预警变成状态监测业务上要讲清楚。5.5 迁移失效换个气候区就崩现象在甲地区训练好的模型直接拿到乙地区预测CSI 从 0.6 降到 0.2。原因模型学到了地区特有的气候条件。湿润区干旱主要由连续无雨日数决定半干旱区则更依赖降水距平百分比如果特征里没有气候区或标准化方式不一致模型无法外推。解决把气候区作为特征或者对每个气候区做单独的标准化。更稳健的做法是加入静态地理特征纬度、海拔等让模型按区域自适应。如果有多个气候区的数据用区域分组交叉验证评估迁移能力。很多开放数据集上效果好的模型最后都是在这一步被刷下来的。6. 让模型真正可用时序交叉验证、决策阈值校准与重训练节奏6.1 用 TimeSeriesSplit 替代 KFoldfrom sklearn.model_selection import TimeSeriesSplit from sklearn.metrics import roc_auc_score tscv TimeSeriesSplit(n_splits5, gap3) fold_scores [] for train_index, test_index in tscv.split(X): X_train_fold, X_test_fold X[train_index], X[test_index] y_train_fold, y_test_fold y[train_index], y[test_index] model_fold RandomForestClassifier( n_estimators300, max_depth6, min_samples_leaf10, class_weightbalanced, random_state42, n_jobs-1 ) model_fold.fit(X_train_fold, y_train_fold) fold_scores.append(roc_auc_score(y_test_fold, model_fold.predict_proba(X_test_fold)[:, 1])) print(f时序交叉验证 AUC: {np.mean(fold_scores):.3f} ± {np.std(fold_scores):.3f})gap3让训练集最后 3 个月不参与训练防止自相关导致信息穿越。这套打分比普通 KFold 更接近真实上线场景。做干旱预测如果时序交叉验证和随机切分的结果相差很大先别急着调模型回去检查特征是否泄漏。6.2 根据代价矩阵校准预警阈值干旱预警中漏报代价远高于误报所以不能把 0.5 作为决策边界。用验证集概率画 Precision-Recall 曲线挑选 F1 最高或业务代价最小的阈值。from sklearn.metrics import precision_recall_curve prec, rec, thr precision_recall_curve(y_val, proba_val) f1_scores 2 * prec * rec / (prec rec 1e-6) best_thr thr[np.argmax(f1_scores)]如果业务里漏报代价是误报的 3 倍可以按代价定义cost misses * cost_miss false_alarms * cost_fa再找到最低 cost 的阈值。这个阈值要每个月根据最新验证集重新校准。6.3 重训练与数据漂移监控气候非平稳模型会随时间漂移。我一般每季度或半年重训练一次具体频率看业务更新周期。同时监控特征分布变化用 PSI 这类指标。当某个特征分布漂移超过阈值比如降水距平均值变化 0.5 个标准差需要把新数据并入训练集并重训。这个环节在开源 GitHub 项目里通常没有但实际落地必须自己补。我现在做干旱预测项目习惯先把数据构造代码和验证方式当作第一优先级模型清单排第二。很多开源项目模型很漂亮但数据泄漏和随机切分的问题一眼就能看到。把这个思路固定下来后我自己的预测系统上线翻车的次数少了很多。希望帮到你。本文还有配套的精品资源点击获取
