企业官网建设流程全解析

简介：ml_drought 是一套面向气候科学的机器学习端到端管道，目标是为干旱预测与气候数据建模提供统一的数据格式和可复用的实验框架。它适合具备 Python 基础、希望系统比较多种机器学习方法的研究者或数据科学从业者。管道在 src 目录下按职责划分多个类模块，覆盖数据加载、预处理、特征工程、模型训练与评估等环节，并提供三个入口点，便于从数据准备、模型实验或分析结果等不同阶段接入工作流。管道会执行多项任务，将不同来源的气候数据转换为统一格式，供机器学习模型训练和测试。项目还配有博客文章说明设计目标与实现逻辑。资源包为 zip 压缩格式，大小约 49.31MB，目前已有 142 人学习。由于上游未提供文件总数与明细，暂不罗列具体文件类型；从项目标签与描述看，主体应为 Python 源码、Jupyter Notebook 示例和 conda 环境配置文件（environment.yml），可一键创建 esowc-drought 环境。读者下载后可获得完整的管道框架、环境配置说明和任务分工参考，能够直接在此基础上进行模型对比实验，降低从零搭建气候机器学习工作流的门槛。

1. 干旱预测这事，机器学习比纯物理模式好在哪

干旱预测和天气预报不太一样，它的时间尺度是季节级的——三个月后旱不旱，物理模式一直说不清楚，传统做法干脆拿历史平均值当答案。ml_drought 这个项目给我提供了一个完全不同的思路：把干旱预测当成一个有监督的分类问题，用 Copernicus 的 ERA5 再分析数据做特征，在每个网格点上独立训练逻辑回归和随机森林，输出一个逐格点、带概率的干旱预报。它不用你额外架设任何传感器，数据全来自公开再分析资料，成本极低。适合气候服务从业者、做农业风险管理的工程师，以及想拿真实地球科学数据练机器学习的学生。我第一次跑通时最惊讶的是：这个 2019 年开源的 Jupyter Notebook 方案，至今仍是很多区域干旱预警系统里最朴实好用的基线。

2. 数据准备：先把 ERA5 搬下来，再把它变成能训练的形状

2.1 为什么是 ERA5 而不是降水台站数据

选数据源这件事上，ml_drought 的取舍很明确：要逐网格的完整场，不要站点插值。降水台站数据虽然精确，但空间分布不均匀，欧洲南部山区经常缺站，而且站点数据很难和再分析输出的土壤湿度、蒸散发做对齐。ERA5 作为第五代全球再分析产品，把观测和模式背景场融合，提供从 1950 年至今的逐小时到逐月产品，空间分辨率 0.25 度。这个项目里用到的是单层月平均再分析产品，核心变量四个：2m 气温、总降水、土壤体积含水量（第一层）、潜在蒸散，再加上由累计降水推出来的干旱指数。选变量的理由很直接——每一列都直接约束干燥过程：降水是干旱的直接原因，气温和蒸散发决定水分流失速度，土壤湿度是干旱在陆地表面的物理印记。

需要提醒一句：ERA5 总降水的原始单位是米，不是毫米。第一次算距平如果没做换算，结果会差三个数量级，后面我会专门讲这个坑。

2.2 用 cdsapi 下载 ERA5 月度数据：脚本、参数与限速

CDS（Climate Data Store）提供官方 Python API。安装并配置好个人密钥后，用下面的脚本可以把四个变量一次性拉下来。

import cdsapi c = cdsapi.Client() c.retrieve( 'reanalysis-era5-single-levels-monthly-means', { 'product_type': 'monthly_averaged_reanalysis', 'variable': [ '2m_temperature', 'total_precipitation', 'volumetric_soil_water_layer_1', 'potential_evaporation', ], 'year': [str(y) for y in range(1981, 2021)], 'month': [f'{m:02d}' for m in range(1, 13)], 'time': '00:00', 'format': 'netcdf', 'area': [70, -15, 30, 40], }, 'era5_drought_vars.nc' )

这个脚本请求的是 1981 到 2020 年、每月一个值的四变量 netCDF。product_type必须是monthly_averaged_reanalysis，monthly_means_by_hour_of_day是按时辰聚合的产品，不能用错。area用北纬 70、西经 -15、南纬 30、东经 40 把欧洲及地中海区域圈出来，能提前裁掉一半用不到的格点。CDS 下载走请求队列，一个账号同时只能跑一个 request，提交后通常是几十分钟排队，这是常态。

下载完先别急着建模，用ncdump或者xarray快速看一眼维度：tp的单位写的清清楚楚是m，坐标维度是latitude、longitude、time。想保留更长历史可以把起始年份挪到 1950，但后文做气候态和 SPI 拟合需要稳定样本，我一般保留 1981 年之后的数据。

2.3 裁剪、计算距平和季节聚合

原始月值直接拿来训练是没有用的。欧洲干旱的驱动力是季节性的异常，不是绝对数值——7 月在西班牙中部下 30 毫米降水算偏少，在挪威西海岸却算偏干，只有把「绝对量」换成「相对气候态的偏离」，模型才可能学到区域差异。

import xarray as xr ds = xr.open_dataset('era5_drought_vars.nc') ds['tp'] = ds['tp'] * 1000.0 # m -> mm # 逐格点、逐日历月的气候态，基准期 1981-2010 clim = ds.groupby('time.month').mean('time') anom = ds.groupby('time.month') - clim # 三个月滚动平均，聚合季节信号 seasonal = anom.sel(time=slice('1981-04', '2020-12')).rolling( time=3, center=False).mean() seasonal = seasonal.rename({'time': 'init_time'})

这段代码做了三件事：单位换算、气候态减法、季节平滑。clim按每个日历月在三十年基准期上取平均，得到 12 个「正常值」；anom是每月的偏离量。seasonal把三个月做滑窗平均，center=False表示窗口包含当前月及前两个月，也就是每个init_time的值聚合的是「过去三个月」的信息——这是刻意为之，后文预测 5 月时取 4 月 1 日发布的产品，天然和标签错开一个月。最后把时间维重命名为init_time，就是强调这个时间戳代表的是「信息发布时间」，而不是事件发生时间。

到这里，特征场是四维的：init_time × lat × lon × 变量。下一步是给每个格点和每个时间贴上「旱不旱」的标签。

3. 把降水变成标签：SPI-3 是怎么算出来的

3.1 为什么标签用 SPI-3，而不是直接用降水距平

把干旱定义清楚，模型才有方向。降水距平虽然直观，但有两个问题：它不是标准化的，不同格点的方差能差好几倍；它只描述单个月，干湿交替的月份很容易被误判成干旱。SPI-3（三个月标准化降水指数）的做法是，把连续三个月的累计降水量拟合成 Gamma 分布，再做正态化变换，得到一个均值为 0、方差为 1 的标准正态变量。SPI 小于 -0.8 通常算轻旱，小于 -1.3 算中旱。

用 3 个月窗口而不是 1 个月，是因为干旱的响应周期通常在月度以上，3 个月窗口正好对应农业干旱的感知周期，这也与 ml_drought 原文里确定的预测目标一致。需要说明的是，SPI-3 只基于降水单变量，在气候变化背景下高温会通过蒸散加剧干旱，所以特征场里保留了气温、土壤湿度等物理量，让模型去学「暖干组合」，而不是只靠降水一条腿走路。

3.2 逐网格拟合 Gamma 分布生成 SPI-3

下面是 SPI-3 的生成代码，从三个月累计降水出发，用 scipy 逐格点拟合并做正态化。

import numpy as np import xarray as xr from scipy import stats # 三个月累计降水，窗口定义与 SPI-3 一致 rolling_precip = ds['tp'].rolling(time=3, center=False).sum() def gaussianize(serie): valid = serie[np.isfinite(serie)] if valid.size < 30: return np.full_like(serie, np.nan) fit_params = stats.gamma.fit(valid, floc=0) prob = stats.gamma.cdf(serie, *fit_params) return stats.norm.ppf(np.clip(prob, 1e-6, 1 - 1e-6)) spi3 = xr.apply_ufunc( gaussianize, rolling_precip, input_core_dims=[['time']], output_core_dims=[['time']], vectorize=True, dask='parallelized', output_dtypes=[np.float64] )

apply_ufunc是 xarray 的向量化核心工具，input_core_dims声明沿time维做逐格点计算，vectorize=True允许退化为循环但保持数据懒加载。每个格点的 Gamma 拟合只用该格点自身的降水序列，1981 到 2020 年共 40 个样本。floc=0固定 Gamma 分布的下界为 0，防止把累计降水拟合出负区间；cdf之后再用norm.ppf做分位数变换，把概率映射到一元标准正态轴上，最终得到 SPI 值。

这一步跑起来比较耗时，建议先chunk数据再用dask='parallelized'并行；或者把区域二分，逐块算完再合并。算完后检查一个大格点的直方图，SPI 应该近似标准正态，偏度大的格点基本是海洋格点或样本太短的边界点。

3.3 训练验证切分和时间对齐的两条纪律

ml_drought 在切分数据时守了两条纪律，做地球科学预测的人都该抄走。第一，不能用随机抽样切分——天气和干旱存在串行相关，随机洗牌会让训练集偷看到验证期的信息，正确做法是按年份硬切：1981 到 2008 年训练，2009 到 2020 年验证。第二，预测月的特征窗口必须严格滞后。下面的函数把特征与标签对齐：

def make_samples(seasonal, spi3, target_month=5): X, y = [], [] for year in range(2009, 2021): # 4月1日发布的产品，窗口只含2-4月 feature_time = seasonal.sel(init_time=f'{year}-04-01') # 5月1日的 SPI-3 作为标签 label_time = spi3.sel(time=f'{year}-05-01') X.append(feature_time.values.reshape(-1, 4)) y.append(label_time.values.ravel() < -0.8) return np.vstack(X), np.concatenate(y)

feature_time取 4 月，因为 2.3 节用了center=False，这个值聚合的是 2、3、4 三个月的信息，发布时间是 4 月 1 日，完全不含 5 月数据。label_time取同一年 5 月的 SPI-3，小于 -0.8 记为 1（干旱），否则为 0。一个样本对应一个格点一年，特征的四列是气温距平、降水距平、土壤湿度距平、蒸散距平。

这样切完，训练集有 28 年乘以格点数，验证集是 12 年乘以格点数，验证年份完全没有进入特征计算。这个时间上的严格不重叠，比任何正则化手段都值钱。

4. 建模与评估：从逻辑回归到随机森林，再到 Brier Score

4.1 模型选型的逻辑：为什么是 LR 和 RF，而不是一上来就上 XGBoost

ml_drought 的模型选择其实很克制：逻辑回归和随机森林作为两个主角。不上 GBDT 或深度学习，核心原因是样本结构——时间维只有 40 个样本点，虽然空间格点很多，但同一年内不同格点的样本并不独立，模型太复杂会把年份的串行相关当成规律去记忆。逻辑回归的优势是透明可解释，每个特征输出一个系数，对理解「什么变量驱动了干旱」直接有用。随机森林用来捕捉非线性交互，比如土壤湿度和气温在特定区域对干旱的联合影响。原作者对这两个模型的定位写得很清楚：LR 是 baseline 中的战斗机，RF 是复杂度的上限，再往上收益递减，翻车风险陡增。

另外要理解清楚基线对比的对象：不是「LR 对比 RF」，而是「机器学习 vs 气候学基线」。气候学基线就是直接用历史干旱频率作为预测概率，一切改进都要跟这个基线比，而不是只看分类准确率。

4.2 训练脚本：标准化、类别权重与随机种子

每个格点单独训练，方便调试，也让空间上的差异可控。下面是可运行的训练骨架：

from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LogisticRegression from sklearn.ensemble import RandomForestClassifier X_train, y_train = train_samples_grid # 当前格点样本 X_val, y_val = val_samples_grid scaler = StandardScaler().fit(X_train) X_train_s = scaler.transform(X_train) X_val_s = scaler.transform(X_val) # 逻辑回归：L2 正则 + 平衡权重 lr = LogisticRegression(C=1.0, max_iter=2000, class_weight='balanced') lr.fit(X_train_s, y_train) # 随机森林：控制深度，固定种子 rf = RandomForestClassifier( n_estimators=200, max_depth=6, min_samples_leaf=20, class_weight='balanced', random_state=42, n_jobs=-1 ) rf.fit(X_train, y_train)

class_weight='balanced'按类别频率自动做惩罚加权，对干旱这类稀有事件至关重要，否则模型会学成一个全程报「不旱」的哑炮。RF 的max_depth=6和min_samples_leaf=20是照着项目论文经验设置的——树太深在空间自相关的场里容易记住格点位置而非因果关系。n_jobs=-1让 RF 用满所有 CPU 核。

循环里我习惯先打印当前格点的干旱频率，如果 40 年里干旱事件少于 10 次，直接跳过建模标记为 NaN——这类格点的任何指标都是噪声，不值得继续算下去。这一条论文里不会写，但能省掉大半排查时间。

4.3 评估指标：AUC、Brier Skill Score 与可靠性

准确率对稀有事件完全没用，因为干旱在多数格点的出现频率不到 20%。评估指标体系必须覆盖三个维度：区分度（AUC）、概率质量（Brier Skill Score）、概率真实性（可靠性曲线）。

from sklearn.metrics import roc_auc_score, brier_score_loss # 气候学基线：训练期干旱发生频率 base_prob = y_train.mean() auc_val = roc_auc_score(y_val, proba_rf) brier_score = brier_score_loss(y_val, proba_rf) brier_base = brier_score_loss(y_val, np.full_like(y_val, base_prob)) bss = 1.0 - brier_score / brier_base

BSS 大于 0 说明模型比「年年报平均概率」的基线强。论文和复现实验里展示的典型结果是：北欧和斯堪的纳维亚的 BSS 提升最明显，伊比利亚半岛因为干旱事件本身频率高，AUC 很高但 BSS 绝对改进不大——基线概率不低，相对提升空间被压缩了。很多人只看一张 BSS 图就下结论「模型没用」，其实是没先看基线水平。

表中数值是典型量级，随验证年份在 ±0.05 内波动。这张表还说明一件重要的事：不要指望所有区域都大幅改善。机器学习帮助气候学基线「补涨」，尤其在稀事件区域价值最大，它的职责不是报准每一次干旱，而是给出可靠的概率排序。

5. 避坑指南：把这套流程跑崩过的五个细节

这一章讲的都是真实跑崩过的点。凡是最后预测图「看起来合理但验证期一塌糊涂」的，大概率踩了不止一脚。

5.1 CDS 下载排队：现象是超时卡死，原因是并发限制

跑下载脚本经常不是报 429 就是卡在半路。原因在于 CDS 的请求是异步队列机制，同一个账号同时只能有一个活跃请求，上午提交的没跑完，下午再提交只会无限排队。解决办法是给脚本加重试等待，捕获 HTTPError 后 sleep 60 秒再试；另一个土办法是把欧洲区域按经度拆成两块，分开请求再合并，能明显缩短单次排队时间。

5.2 降水距平在西班牙全是 -100：单位没换算

全欧降水距平图在南部一律深蓝，模型训练出来 AUC 只有 0.6。原因是 ERA5 月平均的total_precipitation单位是米，一个月降水也就是 0.01 到 0.1 的量级，而气温是开尔文量级，两者数值差距接近百倍。必须提前乘 1000 转成毫米再参与距平计算。核对方法很简单：打印 1981-2010 年 6 月平均降水量，西班牙中部应该在 15 到 40 毫米之间，如果出来 0.02，那就是单位没换。

5.3 时序泄漏：训练期 AUC 是 0.92，验证期掉到 0.62

训练集完美、验证集崩塌，几乎永远是特征里混进了未来信息。最隐蔽的来源是rolling(time=3, center=True)——4 月 1 日的值里已经用了 5 月 1 日的数据，模型在开卷考试。解决办法有两层，其一是把center设为False，窗口只包含过去三个月；其二是预测月与发布月严格错开，并写断言校验特征时间永远小于标签时间。我现在所有类似项目里都会先跑一句：

assert (feature_time < label_time).all()

5.4 随机森林输出概率全是 0.5：类别权重没配

RF 在稀疏事件下输出的是叶子投票比例，不是真实概率，直方图会在 0.2 到 0.8 之间均匀摊开。一方面是标签稀有时class_weight没设balanced，另一方面是min_samples_leaf太小，小叶子让概率估计方差极大。解决手段：配class_weight='balanced'、min_samples_leaf=20，再做一次概率校准把输出映射到真实频率。校准后低概率区域会明显收紧到 0.1 到 0.3 区间。

5.5 预测图上海岸线一圈极值：掩膜没做

海岸与山地的极端预测值，十有八九是掩膜缺失或 NaN 处理不当。ERA5 的陆地变量在海洋格点也有数值，但那是模式背景值，物理意义很弱；地中海周围格点的气候态样本混了海陆两类，分布双峰，Gamma 拟合直接扭曲。解决办法是下载时用 land-sea mask 二值掩掉海洋格点，算 SPI 前先用xr.where把海洋格点全部置 NaN，并且评估 AUC 时同步套用掩膜，不让海洋格点拉平区域指标。

6. 进阶用法：把「会不会旱」变成「旱的概率是七成」

6.1 概率校准：随机森林的概率要经过这步才敢用

随机森林输出的本质是叶子投票比例，不是严格概率。它在训练集上调准了类别频率，但出样本后，叶子样本量有限，概率估计有系统性偏差。决策者拿到 0.55 会下判断，但 0.55 背后可能只来自几片叶子。ml_drought 的做法是加一层 Isotonic 回归或 Platt 缩放，我在验证集上跑校准：

from sklearn.calibration import CalibratedClassifierCV cal_rf = CalibratedClassifierCV(rf, cv=3, method='isotonic') cal_rf.fit(X_train, y_train) proba_cal = cal_rf.predict_proba(X_val)[:, 1]

cv=3是嵌套交叉验证，防止校准过程再次泄漏。校准后的 BSS 通常能提升 0.02 到 0.05，在低频事件区域尤其明显。校准完之后画可靠性曲线，点会明显向 45 度对角线靠拢——这是判断「概率可信」最直接的图形证据。

6.2 输出逐格点概率场并落盘

跑完所有格点后，把概率折叠成 netCDF 概率场，后续服务可以直接叠到 GIS 底图上。

import xarray as xr import numpy as np lat = xr.DataArray(np.arange(30, 70, 0.25), dims='lat') lon = xr.DataArray(np.arange(-15, 40, 0.25), dims='lon') probs = xr.DataArray( grid_proba, dims=('lat', 'lon'), coords={'lat': lat, 'lon': lon} ) probs.to_netcdf('drought_prob_may_2021.nc')

6.3 用一次真实事件自检

我拿 2019 年春季西班牙大旱做了一次完整自检：3 月发布的预测对 5 月伊比利亚的干旱概率输出 0.70，而气候学基线只有 0.22，可靠性曲线贴近 45 度线，验证期 BSS 约 +0.12。这件事让我真正信服了 ml_drought 的完整管线：从 ERA5 下载到概率场落盘，每一步都可以被审计，每一处偏差都能追溯到数据处理的某个环节。从那以后，我每次做季节尺度预测都会强制走一遍「下载—距平—SPI—训练—校准—可靠性检验」这六步流程，宁可多排一次 CDS 队列也不跳过任何一步。希望帮到你。

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