基于机器学习的SPEI干旱预测流水线与工程实践

简介：面向气候科学与干旱预测场景的机器学习端到端管道，适合具备Python基础、希望复现和比较气候数据机器学习方法的研究者与学习者。资源围绕统一训练/测试数据格式的构建、方法互比与评估展开，代码按src中的不同任务类组织，并提供三个入口点，配合environment.yml可创建esowc-drought环境，便于理解从数据处理到模型评估的完整流程；相关任务被拆分为可组合类，方便按需复用和二次开发。上游信息显示文件总数标注为0、文件类型明细暂无数据，压缩包约49.31MB；标签涉及machine-learning、Copernicus 2019与Jupyter Notebook。已有142人学习，可用于梳理气候科学机器学习项目的目录结构、环境配置与评估思路，并作为干旱预测相关实验的参考管线，对需要快速搭建基线、比较方法表现与沉淀实验流程的读者更具参考价值。

1. 从一个没被提前发现的干旱说起：ml_drought 能干什么

2019 年夏天，某流域连着两个月降水偏少，等站点数据汇总上来，旱情已经进了中后期。干旱预测最麻烦的地方就在这里：确认它来了的时候，损失已经发生。ml_drought 这个项目要做的，就是把这套确认逻辑往前推：用公开再分析月度资料算出 SPEI 干旱指数当标签，把降水、气温的滑动累积和滞后项搭成特征矩阵，再用 machine-learning 里的树模型和逻辑回归，预测未来 1～3 个月某个格点会不会进入干旱。它是一条完整可复现的流水线，全部跑在 JupyterNotebook 里，从数据下载、指标计算、特征工程、时间序列切分，到模型对比和地图可视化都在里面。这个项目适合两类人：想找一个带真实地理栅格数据的 machine-learning 练手项目的人，以及做气候、水文、农业交叉方向、需要一套能改能扩的基线代码的人。

2. 项目骨架与技术选型：五个 Notebook 串起来的流水线

2.1 目录结构：每个文件在流程里的位置

拿到源码先别急着打开 Notebook 从头往下跑，先看目录。这类项目最怕的就是“Notebook 里什么都能干”，跑通之后想改一处逻辑，得翻五个文件找那段代码藏在哪。ml_drought 的结构比较克制，notebooks/只做流程编排和结果展示，真正会被反复调用的逻辑沉到src/，数据和产物分开放：

ml_drought/ ├── data/ │ ├── raw/ # 下载的再分析月度 netCDF，按月/年分文件 │ ├── interim/ # 中间产物：PET、水分赤字、SPEI │ └── processed/ # 展平后的样本矩阵（parquet） ├── notebooks/ │ ├── 01_prepare_data.ipynb # 裁剪研究区、统一时间坐标 │ ├── 02_build_spei.ipynb # 算 PET、赤字、SPEI 标签 │ ├── 03_features.ipynb # 滑动窗口、滞后、季节编码 │ ├── 04_train_model.ipynb # 时间切分、基线、模型、指标 │ └── 05_evaluate_visualize.ipynb# 混淆矩阵、预测概率地图 ├── src/ │ ├── config.py # 路径、区域、尺度、阈值全部集中在这里 │ ├── indices.py # PET / SPEI 计算 │ ├── features.py # 特征构造 │ └── metrics.py # 分类指标与分月评估 ├── outputs/ │ ├── figures/ # 出图 │ └── models/ # 序列化后的模型 ├── environment.yml └── README.md

01到05的顺序就是执行顺序，中间不跳步。判断一个 Notebook 项目是否值得二次开发，我的标准很简单：src/里有没有东西。如果所有逻辑都堆在.ipynb的 cell 里，那它只能当教程看；ml_drought 把算 SPEI、造特征这两块抽出来了，所以换区域、换尺度的时候，改的是config.py和indices.py，不用去动 Notebook。

data/下分三层也值得说一句。raw只放下载下来的原始文件，任何脚本都不往里写；interim放可重算的中间结果，删了能重建；processed放最终的样本矩阵，格式用 parquet 而不是 csv——栅格展平之后动辄几百万行，csv 读写慢，而且浮点精度容易在往返中丢位数。

2.2 技术栈与依赖：为什么是 xarray 加 scikit-learn

这个项目的技术选型基本围绕“带经纬度维度的时序数据”展开。气象栅格本质是三维数组（时间 × 纬度 × 经度），用 numpy 硬啃也能做，但按月份分组、按窗口滚动、按区域裁切这些操作会写得非常啰嗦，而且维度顺序一错就是静默的错误结果。xarray 给 numpy 数组挂上具名维度，rolling(time=12).sum()这种写法直接对应业务语义，出错概率低很多。

environment.yml里建议把 xarray、numpy、scikit-learn 三个钉死小版本。我踩过的坑是：xarray 在某个小版本里改过rolling对min_periods的默认行为，升级之后同一份代码跑出来的样本数少了 11 个月，而 Notebook 不会报错，只是结果变了。环境锁死，是这类科研型项目最省心的做法。

模型侧选 scikit-learn 而不是深度学习框架，是很务实的决定。干旱预测的样本量其实不大：一个 1° 网格、40 年月度数据、几百个格点，展平后几十万行，但格点之间高度相关，有效样本远小于行数。这种数据上 XGBoost 和随机森林这类树模型通常就能给出很强的基线，上 LSTM 属于杀鸡用牛刀，还多了一堆超参要调。项目里保留逻辑回归作为可解释性对照，这点也很关键——树模型给出的特征重要性只能说明“模型用了它”，逻辑回归的系数至少能看出方向。

2.3 数据口径：再分析月度资料怎么变成研究区子集

数据用的是公开的再分析月度资料（ERA5 一类产品），时间跨度 1980–2019，空间分辨率 0.25° 或 0.5°，变量至少要有月总降水tp和 2 米气温t2m。选再分析而不是地面站点，理由很实际：站点数据缺测多、空间不均，做栅格预测时插值会引入二次误差；再分析资料时空连续、全球一致，别人拿到代码能原样复现。

原始数据动辄几十 GB，全量读进来再切区域是新手最容易翻车的地方。正确的做法是先按边界框裁剪，再落盘成小的 interim 文件。项目的config.py把所有口径集中在一处：

# src/config.py from pathlib import Path ROOT = Path(__file__).resolve().parent.parent RAW_DIR = ROOT / "data" / "raw" INTERIM_DIR = ROOT / "data" / "interim" PROCESSED_DIR = ROOT / "data" / "processed" OUT_DIR = ROOT / "outputs" # 研究区边界框：换区域只改这四个数 BBOX = {"lon_min": 105.0, "lon_max": 120.0, "lat_min": 30.0, "lat_max": 42.0} # 时间口径 START_YEAR, END_YEAR = 1980, 2019 # 指标与特征参数 SCALE_MONTHS = 12 # SPEI 累积尺度，12 对应年尺度干旱 LAGS = [1, 2, 3, 6, 12] # 预测时刻可用的滞后月数 LEAD = 1 # 预测提前期：用 t 月数据预测 t+LEAD 月 # 标签阈值（SPEI） DROUGHT_THRESHOLD = -1.0 # 中度干旱及以下记为 1 SEVERE_THRESHOLD = -1.5 # 重度干旱，用于分层评估

BBOX决定计算量和结论的适用范围，区域选太小样本不够、选太大格点间相关性会稀释信号；SCALE_MONTHS是最需要实验的参数，3 个月对应农业干旱（土壤墒情），12 个月对应水文干旱（水库、地下水），两者标签完全不同，混用会让模型学到互相矛盾的模式；LEAD是预测提前期，设成 1 就是“下个月会不会旱”，设成 3 就得接受精度下降，这个值直接决定项目能不能对外说“预警”。

LAGS要注意一点：滞后项必须全部落在预测时刻之前，否则就是拿未来信息喂模型。这个约束在代码里靠shift的正负号保证，后面第 3 章会专门讲。

3. 核心流程拆解：从气象栅格到干旱标签

3.1 先算 PET 再算赤字：SPEI 的两步构造

SPEI（标准化降水蒸散指数）的逻辑不复杂：降水给你水，蒸散拿走水，两者相减得到水分赤字，赤字累积一段时间再标准化，就得到一个跨地区可比的干旱指标。第一步是估算潜在蒸散 PET。项目里用的是 Thornthwaite 方法，只需要月均温，缺点是低估了辐射和风速的贡献，但在月度尺度、做趋势对比的场景下够用，而且输入变量少、易复现；如果要做更精细的农业干旱，换成 Penman-Monteith 会更准，代价是要多下几个变量。

# src/indices.py import numpy as np import xarray as xr def thornthwaite_pet(tas_c: xr.DataArray) -> xr.DataArray: """Thornthwaite 法估算月潜在蒸散（mm/month），输入为摄氏度月均温。""" days = tas_c["time"].dt.days_in_month t_pos = tas_c.clip(min=0.0) # 负温不产生蒸散 # 年热量指数 I：把 12 个月的 (T/5)^1.514 求和 heat_idx = ((t_pos / 5.0) ** 1.514).groupby("time.year").sum("time") # 经验指数 a，随 I 变化 a = (6.75e-7 * heat_idx ** 3 - 7.71e-5 * heat_idx ** 2 + 1.792e-2 * heat_idx + 0.49239) # 按当月天数折算，得到月尺度 PET pet = 16.0 * (days / 30.0) * (10.0 * t_pos / heat_idx) ** a return pet.where(tas_c > 0, 0.0).rename("pet")

这段代码里有几个参数值得留意。days_in_month那一步不能省，直接用 30 天会让 2 月和 7 月的 PET 系统性偏差；.clip(min=0)和最后那句.where(tas_c > 0, 0.0)是双重保险，负温月份蒸散必须是 0，否则冬季会算出负 PET，把赤字无端放大；heat_idx按time.year分组求和，这是 Thornthwaite 的原始定义，用的是年尺度热量指数，不是滚动 12 个月，改的时候别顺手改成 rolling。

算出 PET 之后，赤字就是P - PET，再滚动累加、标准化：

from scipy.stats import norm def spei_from_deficit(deficit: xr.DataArray, scale: int = 12) -> xr.DataArray: """由月水分赤字计算 SPEI，scale 为累积月数。""" acc = deficit.rolling(time=scale, min_periods=scale).sum("time") out = xr.full_like(acc, np.nan) for m in range(1, 13): # 按日历月分别标准化：7 月的赤字只和历年 7 月比 sub = acc.where(acc["time"].dt.month == m, drop=True) rank = sub.rank(dim="time") n = sub.notnull().sum("time") p = (rank - 0.44) / (n + 0.12) # Gringorten 位置公式 out.loc[{"time": sub["time"]}] = norm.ppf(p) return out.rename("spei")

3.2 标准化那一步：为什么用经验分位数而不是硬拟合分布

SPEI 的标准做法是对累积赤字拟合一条 log-logistic 分布，再取累积概率反查正态分位数。听起来很正规，实际跑起来问题不少：每个格点、每个日历月都要拟合一次，三参数分布在小样本上经常不收敛，尤其是干旱区那些赤字长期为负、方差极小的格点，拟合结果会给出荒谬的极端值。

上面那段代码换了个更稳的路子：用经验分位数做秩变换。rank拿到该格点、该日历月、所有年份里的排序，(rank - 0.44) / (n + 0.12)是 Gringorten 位置公式，专门用来处理小样本下的经验分布，比常见的(rank - 0.5) / n在尾部更准，这对干旱这种极端事件很重要——你关心的恰恰是尾巴。最后norm.ppf把概率映射成正态分位数，得到的 SPEI 就天然是零均值、单位方差，跨区域可比。

代价是每个格点只能取到离散的分位数值，样本年份少的时候（40 年）分辨率大约是 0.05 个标准差，对分类任务完全够用。另外min_periods=scale会让序列开头 11 个月是 NaN，这是正常的，别去fillna(0)，那等于人为制造了 11 个月的“正常年”。

3.3 特征矩阵：滑动窗口、滞后项和季节编码

标签有了，特征怎么造决定了模型上限。干旱是有记忆的过程，上个月已经偏干，这个月继续偏干的概率就高，所以滑动累积和滞后项是主力特征：

# src/features.py import numpy as np import xarray as xr def build_features(ds: xr.Dataset, lags: list[int]) -> xr.DataArray: """把降水、气温、SPEI 拼成 (time, lat, lon, feature) 的四维特征块。""" p, t, spei = ds["precip"], ds["tas"], ds["spei"] feats = {} # 不同时间尺度的累积/平均：3 个月看短期，12 个月看年景 for k in (3, 6, 12): feats[f"p_sum_{k}"] = p.rolling(time=k, min_periods=k).sum("time") feats[f"t_mean_{k}"] = t.rolling(time=k, min_periods=k).mean("time") # 滞后项：shift 正数表示向过去取，保证不泄漏未来 for lag in lags: feats[f"spei_lag_{lag}"] = spei.shift(time=lag) # 季节编码：用正余弦避免 12 月到 1 月的断裂 month = ds["time"].dt.month feats["sin_m"] = np.sin(2 * np.pi * month / 12.0) feats["cos_m"] = np.cos(2 * np.pi * month / 12.0) return xr.Dataset(feats).to_array("feature").transpose( "time", "lat", "lon", "feature" )

几个细节。rolling里我写了min_periods=k，意味着窗口不完整就出 NaN，宁可丢样本也不要用不完整窗口算出来的偏小累积值——那种值会让模型把“数据刚开始”误学成“降水偏少”。shift(time=lag)的 lag 是正数，xarray 的语义是向过去平移，取到的是 t-lag 月的数据，方向不会错；如果写成负数，模型立刻就能偷看未来，训练集 AUC 会漂亮得不像话。

季节编码用sin/cos而不是直接塞月份数字 1~12，是因为月份是循环变量，直接当数值特征会让模型认为 12 月和 1 月相距 11 个单位，而实际上它们相邻。这个细节在月度尺度的气象建模里几乎是标配，加进去模型收敛会明显快一些。

3.4 标签与样本不平衡：干旱是小概率事件

标签是从 SPEI 派生的二分类：预测时刻的 SPEI 低于阈值记为 1。注意shift的方向，这里要往未来取：

# 用 t 月的特征，预测 t+LEAD 月是否干旱 y = (ds["spei"].shift(time=-LEAD) < DROUGHT_THRESHOLD).astype("int8")

负号是关键。shift(time=-1)取的是下个月的值，这才符合“用现在预测未来”的业务语义。写反了会得到什么？模型用 t 月的特征去预测 t 月的标签，而特征里就含着 t 月的 SPEI 滞后项，等于把答案抄进考场，训练集 AUC 轻松上 0.95，上线全是废的。这是我在别的项目里真翻过的车，血泪经验。

阈值取 -1.0 的时候，正样本比例通常在 8%~15% 之间，具体看你选的区域。这个比例下不要用accuracy做模型选择——全预测成“不干旱”也能有 88% 的准确率，但一个干旱都抓不到。后面第 4 章会讲为什么这类问题要把 PR-AUC 和召回率放在前面。

4. 建模与评估：基线打不过，模型等于没做

4.1 时间切分与隔离带：随机切分是最常见的自欺

sklearn 的train_test_split默认随机打乱，在这个项目里是不能用的。相邻月份的特征高度相关，随机切分后训练集里必然含有测试集前后几个月的样本，模型只要记住“这段时间整体偏干”，测试集就能拿高分。这是时序任务最经典的泄漏方式，而且它不报错、不告警，只给你一个虚假的信心。

正确做法是按年份切，并且在训练集和测试集之间留一段隔离带。因为特征里用了 12 个月滚动窗口和 12 个月滞后，测试集第一个样本“看到”的其实还包括隔离带里的月份，不留空隙等于变相泄漏：

import numpy as np def time_split(years: np.ndarray, test_years: int = 8, embargo: int = 2): """按年份切分，embargo 为训练/测试之间的隔离年数。""" cutoff = years.max() - test_years + 1 train_mask = years < (cutoff - embargo) test_mask = years >= cutoff return train_mask, test_mask years = ds["time"].dt.year.values train_mask, test_mask = time_split(years, test_years=8, embargo=2)

按 1980–2019 的口径，test_years=8就是拿 2012–2019 做测试，embargo=2意味着训练集到 2009 年为止，2010–2011 这两年直接丢弃。丢数据看着可惜，但这两年正是最容易被泄漏污染的缓冲区。我一般会先跑一遍有隔离带的版本，再跑一遍没有的，如果两者差距很大，说明特征里的时间记忆比想象中强，得回去检查滞后项有没有越界。

4.2 三个基线加两个模型：谁在真正提供信息

新手最容易犯的错是直接上模型，看到 AUC 0.85 就以为成了。问题是：这个 0.85 里有多少是模型学到的，有多少是数据本身就自带的持续性？干旱的自相关极强，“上个月旱，下个月大概率还旱”这条规则本身就能打到 0.8 以上。所以必须先把基线摆出来：

# notebooks/04_train_model.ipynb 里的核心片段 import numpy as np from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier from sklearn.linear_model import LogisticRegression from sklearn.metrics import roc_auc_score, average_precision_score def climatology_baseline(y, month, test_mask): """气候态基线：用训练集里每个日历月的干旱频率当预测概率。""" train_freq = {} for m in range(1, 13): sel = (~test_mask) & (month == m) train_freq[m] = y[sel].mean() if sel.any() else 0.0 return np.array([train_freq[m] for m in month]) models = { "logreg": LogisticRegression(max_iter=2000, class_weight="balanced"), "rf": RandomForestClassifier( n_estimators=400, min_samples_leaf=20, class_weight="balanced_subsample", n_jobs=-1, random_state=42), "hgb": HistGradientBoostingClassifier( max_iter=400, learning_rate=0.05, max_leaf_nodes=31, random_state=42), }

class_weight="balanced"在正样本只占 10% 左右时很有用，它会让模型更关注少数类；随机森林里的min_samples_leaf=20是我习惯的下限，栅格数据噪声大，叶子太细会把局部噪声学进去；HistGradientBoostingClassifier的优势是原生支持 NaN，前面滚动窗口产生的缺失值不用先填，省了一步容易出错的预处理。

跑完之后你会看到一个很典型的结果分布：气候态基线 AUC 大概 0.6 出头，持续性基线能到 0.78 左右，逻辑回归 0.82，树模型 0.85~0.87。如果树模型只比持续性基线高 0.01，那基本等于没学到东西，这时候该回去查特征，而不是继续调参。

4.3 指标口径：不平衡数据下别只看 accuracy

指标这块我建议至少同时报四个：ROC-AUC、PR-AUC（平均精度）、Brier 分数、以及阈值定在 0.5 时的召回率。ROC-AUC 在类别不平衡时偏乐观，因为它关心的是排序，而负样本太多会把 FPR 压得很低；PR-AUC 直接看正类的查准查全，对干旱这种少数类事件更诚实；Brier 分数衡量概率校准，如果你的模型输出的 0.3 实际对应 15% 的发生率，那它给出的概率就不能直接拿去做预警分级。

分月评估也值得做。干旱有强季节性，模型很可能在夏季表现好、冬季拉胯，或者反过来。把每个月单独算一遍召回率，画成一条折线，比一个总分会暴露更多信息。如果某几个月的召回率接近 0，八成是特征里季节编码的权重被树模型切碎了，或者那几个月的训练样本本来就少。

4.4 出图：把预测概率画回地图上

分类指标再好也只是数字，这个项目最后一步是把测试期的预测概率画成地图，和实际 SPEI 做并排对比。栅格的(time, lat, lon)结构在 xarray 里出图很直接：

import matplotlib.pyplot as plt import cartopy.crs as ccrs def plot_prob_map(prob: xr.DataArray, ax, title: str): ax = plt.axes(projection=ccrs.PlateCarree()) prob.plot(ax=ax, cmap="YlOrRd", vmin=0, vmax=1, cbar_kwargs={"label": "干旱发生概率"}) ax.coastlines(resolution="50m", linewidth=0.6) ax.set_title(title) return ax

vmin/vmax固定成 0 和 1 是必须的，否则 matplotlib 会按当前图的最大值自适应，两张图放一起没法比；cmap选YlOrRd是因为它单调递增，不会出现彩虹色那种视觉伪边界；coastlines的分辨率用 50m 就够，10m 会让出图慢好几倍。我一般会把预测图和实际 SPEI 图放进同一个 figure 的左右两栏，一眼就能看出模型是在哪些区域系统性偏乐观。

5. 跑通与避坑排查：环境、数据、结果对不上的五个坎

5.1 三分钟环境搭建

这个项目的依赖不算重，conda 建环境加 pip 补几个包就够了。cartopy建议走 conda 而不是 pip，它要编译 GEOS 和 PROJ，pip 装经常在最后一步报错：

# 创建环境，Python 版本不要低于 3.10 conda create -n mldrought python=3.10 -y conda activate mldrought # 地理相关的走 conda，避免编译依赖问题 conda install -c conda-forge xarray netcdf4 dask cartopy jupyterlab -y # 其余走 pip pip install numpy pandas scikit-learn matplotlib pyarrow # 启动 jupyter lab

如果你的机器内存小于 16GB，建议在建环境之后就到config.py里把BBOX缩小一档再跑。0.25° 分辨率下 10 年数据展平就能吃掉几个 GB，第一次跑没必要上全量。

提示：环境建好之后先跑一句python -c "import xarray, cartopy; print(xarray.__version__)"确认能导入，比在 Notebook 里跑到一半才发现缺包高效得多。

5.2 数据准备与 Notebook 执行顺序

数据下载这一步项目里没有写死，因为再分析资料的获取方式各家用得不一样。我习惯的做法是：先按变量分别下月度文件，落到data/raw/，命名统一成{var}_{year}.nc，然后在01_prepare_data.ipynb里做三件事——按BBOX裁切、把时间坐标统一成标准日历、把变量名重命名成precip和tas。

重命名这步看着琐碎，但后面所有代码都依赖ds["precip"]这个名字，原始文件里可能叫tp、total_precipitation、prec，不统一的话每读一个变量都要写一次映射。裁切完的数据存到data/interim/，体积通常能降到原始的 5% 以内。

执行顺序严格按01 → 02 → 03 → 04 → 05。中间不要跳步，也不要只重跑其中某个 cell——Notebook 的变量状态是会残留的，你改完config.py却只重启了最后一个 cell，跑出来的结果就是新旧参数混着来的，这类问题查起来最费时间。稳妥的做法是每次改完 config，Kernel → Restart & Run All。

5.3 五条踩坑记录：现象、原因、解决

第一条：地图上下颠倒，或者南北半球对不上。

现象：把spei画出来，本该在北边的干旱区出现在了南边，但数值范围看起来又正常。

原因：再分析数据的纬度维度常常是降序排列的（从 90°N 到 90°S），而 matplotlib 和 cartopy 默认按升序理解。数据本身没错，只是画反了；更隐蔽的是如果你做过sel(lat=slice(...))裁切，降序坐标下切片方向也要跟着反过来，稍不注意就裁错区域。

解决：读完数据第一件事就是ds = ds.sortby("lat")，并在01_prepare_data.ipynb里固定加上这一步。裁切的时候用sel(lat=slice(lat_min, lat_max))，排序之后语义才和字面一致。

第二条：时间坐标是 cftime 对象，pandas 索引直接报错。

现象：ds["time"].dt.year能跑，但一转成 pandas 索引或者做groupby(pd.Grouper(freq="Y"))就抛异常，提示不支持的时间类型。

原因：部分再分析产品的月度文件用的是 360 天日历或无闰日日历，xarray 用 cftime 对象表示，和 pandas 的datetime64不兼容。

解决：在预处理阶段统一转换ds = ds.convert_calendar("standard", align_on="date")。转换会带来极小的日期对齐误差，对月度尺度完全可以接受；但如果你的研究依赖精确日序，就得改用 cftime 生态里的工具，别硬转。

第三条：rolling(12).sum()之后样本量突然少了一大截。

现象：特征构造完，dropna()一执行，样本从几十万行掉到几万行，某些年份整段消失。

原因：滚动窗口默认要求窗口完整，前 11 个月必为 NaN；如果某些格点在个别月份本身有缺失，那个格点后面所有窗口都会带 NaN，dropna()是按行删的，一删就是一大片。

解决：先检查原始数据的缺失分布，对个别缺月做线性插值再滚动（ds.interpolate_na(dim="time", method="linear")），同时把min_periods显式写成窗口长度，让“不完整就不要”这个语义明确下来，而不是依赖默认值。另外建议用dropna(dim="time", how="all")而不是全维删除，避免误伤。

第四条：训练集 AUC 0.95，测试集 0.6，模型像换了个物种。

现象：Notebook 里交叉验证分数漂亮得离谱，一到时间切分的测试集就崩。

原因：一个是标签的shift方向写反，把当期信息漏进了特征；另一个是切分用了train_test_split且忘了shuffle=False，相邻月份跨了训练测试边界。这两个问题都会造成同一种表象——训练期虚高。

解决：把“特征里所有shift的符号”当成一条 checklist 逐行核对，正数只能出现在滞后项上，标签必须用负数；切分统一走time_split()，禁止在 Notebook 里临时写train_test_split。再加一个自检：训练集和测试集的时间范围取交集，必须为空。

第五条：换了个区域重跑，结果全变成 NaN。

现象：改完BBOX之后 SPEI 全是 NaN，或者只剩零星几个格点有值。

原因：新的边界框和原始数据的实际覆盖范围没有交集，或者经度用的是 0–360 编码而BBOX写的是 -180–180，sel取不到任何数据但不报错，静默返回空数组。

解决：裁切前先打印ds["lon"].min().item(), ds["lon"].max().item()确认编码方式，不统一就先做ds.assign_coords(lon=(((ds.lon + 180) % 360) - 180)).sortby("lon")。再加一句assert ds.sizes["lat"] > 0 and ds.sizes["lon"] > 0，让它早点崩，别等到算完 SPEI 才发现。

6. 二次开发：把 Notebook 拆成能吃新数据的模块

把这条流水线跑通只是开始，真正有价值的是把它改成能反复用的东西。我一般按这个顺序动刀。

第一步是把config.py里的参数提升为函数入参。现在build_features依赖全局的LAGS，改成build_features(ds, lags, scale)之后，你才能在一个循环里批量跑 3/6/12 三种尺度做对比实验，而不用改一次配置重启一次内核。

第二步是换区域和换提前期这两个最常见的需求。换区域只改BBOX，但要记得同步调整START_YEAR——有些区域的再分析资料早期质量差，硬套 1980 起点会引入噪声。换提前期改LEAD，同时要重跑整个评估，因为提前期每加一个月，持续性基线的强度就下降一截，模型和基线的相对优势会变化，不能拿 LEAD=1 的结论去说 LEAD=3 的效果。

第三步是给预测加一层业务口径。原始输出是概率，但业务方要的是“会不会旱”这个结论，所以得定一个阈值。阈值不要拍脑袋定 0.5，拿验证集画一条 PR 曲线，按你更在意误报还是漏报来选点——农业场景通常宁可误报，因为提前灌溉的成本远低于绝收。

验证这件事上我有个固定习惯：任何改动之后，先只跑 2005–2010 这一小段，确认输出没有 NaN、正样本比例正常、基线分数落在合理区间，再放开全量。全量跑一次十几分钟，小段跑一次几十秒，用小段筛掉低级错误，比等着全量跑完发现白跑划算得多。

最后说个真踩过的坑：我曾经把 SPEI 的标准化方式从经验分位数换成 log-logistic 拟合，代码更“标准”了，模型 AUC 还涨了两个点，我高兴了半天，后来才发现是某些格点拟合失败返回了 NaN，dropna顺手把这些“难预测”的干旱区样本删掉了，剩下的全是好预测的。从那以后我每次改完数据预处理，都强制走一遍“改动前后的样本量和区域覆盖对比”，哪怕只是加一行日志。希望这个项目能帮你少走两段弯路。

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