简介:ml_drought是一套面向气候科学的机器学习端到端管道,聚焦干旱预测与模型对比研究,适合气候科研人员、环境数据分析者及有一定Python基础的开发者。管道通过src目录下的多个任务类,把数据格式统一、特征构建、模型训练与评估等环节标准化,并设置三个独立入口点,用户可按需调用或替换模块。压缩包大小约49.31MB,核心内容以Jupyter Notebook和Python脚本为主,随附environment.yml环境配置,可用Anaconda一键创建运行环境。已有142人学习下载,可据此快速上手:既能复现实验流程,也能将模块改造后接入自己的干旱数据集,其清晰的模块边界也适合作为科研项目的起点,或用于教学演示与算法基准测试。整体设计兼顾可读性与扩展性,是一份实用且便于二次开发的气候机器学习工具包。
1. 机器学习预测干旱:从数据到决策的最后一公里
ml_drought 是 GitHub 上 ml-clim 组织下的一个开源项目,核心就一句话:用机器学习把“干旱预测”这件事从统计描述推进到空间网格上的定量估计。很多做气象和农业的朋友第一反应是,干旱不是有 SPI、SPEI 这些指数了吗?但那些指数算的是“过去一段时间的干湿状况”,属于事后描述。真正难的是提前几个月知道哪个区域会偏旱、偏旱到什么程度。ml_drought 要解决的就是这个提前量问题,它把干旱预测包装成一个有监督的时空回归任务,输入历史气候场,输出未来干旱指数。这个项目适合三类人:搞气候服务的技术人员、做农业保险和粮食安全评估的算法工程师,以及想找一个非计算机视觉场景练手时空建模的研究生。先说清楚,机器学习不会直接告诉你“哪里有干旱”,它预测的是一个可验证的目标量,这个目标量怎么定义、怎么准备数据,才是整个项目真正花时间的地方。
2. 干旱预测为什么需要机器学习:从指数到时空模型的跨越
2.1 传统干旱指数的局限:SPI 是记账,不是天气预报
做干旱这件事,第一步通常是计算干旱指数。最常见的标准化降水指数(SPI)和标准化降水蒸散发指数(SPEI),本质都是把降水或降水减蒸散序列做概率分布标准化。计算流程很简单:对某个站点过去几十年的逐月降水序列,分别对 1 个月、3 个月、6 个月窗口做滑动累加,再对累计值做 Gamma 分布拟合,最后映射到标准正态分布。下面这段代码就是用 scipy 算 SPI 的常规做法。
import numpy as np import pandas as pd from scipy import stats def compute_spi(precip: pd.Series, timescale: int = 3): # 滑动累加,得到过去 timescale 个月的累计降水 roll = precip.rolling(window=timescale, min_periods=timescale).sum() # 剔除 NaN,用最大似然估计拟合 Gamma 分布参数 valid = roll.dropna() fit_alpha, fit_loc, fit_beta = stats.gamma.fit(valid) # 对每个累计值求累计概率,再转成标准正态分位数 cdf = stats.gamma.cdf(roll, fit_alpha, loc=fit_loc, scale=fit_beta) spi = stats.norm.ppf(cdf.clip(lower=1e-6, upper=1 - 1e-6)) return spi这段代码的逻辑不复杂,但问题也在这。SPI 用的是“已经发生的降水”,滚动窗口结尾落在当前时刻,所以它回答的是“过去三个月够不够干”。你想要“未来三个月会不会干旱”,用 SPI 直接外推只能靠自回归,或者把气候模式的降水预报喂进去再算一个预报版 SPI。这就是 ml_drought 这类项目的切入点:与其手工把数值模式输出转成指数,不如让模型直接学习“过去气候场 → 未来指数”的映射。这里的要点是,ml_drought 的标签仍然可以用 SPI 或 SPEI 构建,但输入特征不再局限于降水,而是温度、气压、土壤湿度、海温等多个变量组成的网格场。
从模型角度看,气候网格数据是典型的二维空间结构,加上时间维度就是三维张量。GraphCast、PanguWeather 这类大模型用的也是类似思路,只是尺度更大。ml_drought 更聚焦干旱场景,输入输出都是区域网格,分辨率通常在 0.25 到 1 度之间。很多人听到卷积神经网络就以为那是计算机视觉专用,其实 CNN 处理网格气候场非常自然,卷积核就是在捕捉空间上的局地相关性。机器学习与计算机视觉的边界在这里很有意思:同样的 U-Net 结构,CV 里用来分割道路,干旱模型里用来把多个气候变量编码成对未来干旱指数的预测场,问题定义比网络结构本身更关键。
2.2 把干旱预测变成有监督回归问题:输入、输出与损失
ml_drought 的建模套路,我按自己做类似项目的经验拆解一下。输入是一个历史时间窗口,比如过去 12 个月的月均变量场,每个变量是一张 H×W 的图,多变量叠加后变成 C×H×W,时间维可以并到通道维,也可以用 ConvLSTM 保留时间结构。输出是未来某个时间尺度的 SPEI 或 SPI 场,比如未来 3 个月平均的 SPEI-3。这样就不需要做站点级别的预测再插值,而是直接学习网格到网格的映射。
# 数据张量形状设计,以 TensorFlow 为例 # inputs: [batch, time, channels, height, width] # 12个月,3个变量(降水、温度、蒸散发),0.5度网格 # inputs.shape = (32, 12, 3, 128, 256) # targets: [batch, height, width],未来3个月SPEI-3 # targets.shape = (32, 128, 256)这个设计里有三个坑,直接影响模型能不能训练起来。第一个是时间窗和预测窗不能重叠。如果你用第 t 个月往前 12 个月的资料预测第 t+1 到 t+3 个月的指数,那么样本划分时要避免把标签所在时段也放进输入特征,否则数据泄漏会让验证集指标一片祥和。第二个是目标变量的空间平滑性。SPEI 在网格上本身有很强空间自相关,模型很容易学到“抄邻居”的捷径,也就是输出的场和输入最后几帧几乎一样,这在气象上叫持续性预测。所以损失函数不能只看整体 RMSE,还要看异常区域的命中率。第三个是多变量尺度差异。降水可能是 0 到 200 毫米,气温是 -30 到 40 摄氏度,不先做标准化,模型训练会不稳定。Ridge 回归里那个惩罚项就相当于隐式标准化,但在深度学习里最好显式处理。
标准的回归损失可以直接用 MSE,但对干旱这种极端事件建模,我更推荐在损失里加一个权重项,让模型更关注尾部样本。比如对每个像素,如果真实 SPEI 小于 -1(达到轻旱),就把该像素的 MSE 损失乘以 2.0。实现很简单,就是在 NumPy 和 TensorFlow 里做个 mask 相乘。
import tensorflow as tf def weighted_mse(y_true, y_pred, drought_threshold=-1.0, weight=2.0): squared_error = tf.square(y_true - y_pred) # 对干旱像素提高损失权重,让模型不把精力全花在“正常”区域 drought_mask = tf.cast(y_true < drought_threshold, dtype=tf.float32) return tf.reduce_mean(squared_error * (1.0 + weight * drought_mask))这个加权损失要配合训练日志里的分桶指标看,否则容易变成模型把所有像素都预测成一个接近干旱的值,因为那样也能降低加权损失。后面在评估章节再展开讲。总之,ml_drought 的建模逻辑本身并不神秘,它就是给物理气候问题套上标准的机器学习壳子,但壳子里面装的数据组织和目标定义,才是决定效果上限的地方。
3. 跑通 ml_drought 的最小环境与数据准备
3.1 环境依赖:先解决“能跑”和“能复现”
从 GitHub 上把仓库克隆下来只是第一步。ml_drought 这类研究型项目对 Python 环境比较敏感,最常见的翻车就是本地 Python 版本和项目依赖冲突。我一般会先用 conda 建一个独立环境,再按仓库里的 requirements 或 environment.yml 安装依赖。注意,不能直接用 base 环境去试,因为 xarray、dask、tensorflow 这几个库的版本矩阵很容易打架。
# 创建独立环境,Python 3.9 是目前气候机器学习项目兼容性最好的版本 conda create -n ml_drought python=3.9 -y conda activate ml_drought # 安装基础地理和科学计算库 pip install xarray netcdf4 dask scipy pandas matplotlib # 安装深度学习框架,CPU 机器先装上能跑通调试 pip install tensorflow-cpu这里特别说一下为什么要 xarray。气候数据几乎都是 NetCDF 格式,包含经纬度坐标和时间坐标,普通 NumPy 数组处理起来非常容易搞错维度顺序。xarray 的 DataArray 把坐标和变量绑在一起,晚点做按区域切片、按时间聚合、与观测对齐都省心。如果你跑的是 GPU 版本的 TensorFlow,避免和 CPU 版混装,卸载干净再装。一个实用经验:先装 xarray 和 dask,再装 tensorflow,最后装项目里的其他依赖,这样 pip 不会因为顺序问题把某个包降级。
依赖装完之后,最好用 conda list 把关键包版本记录到一个 requirements-lock.txt 里。气候数据处理的库更新很快,半年后你再跑同一个项目,可能因为 NumPy 版本变化导致随机种子对不上,结果无法复现。这个锁文件就是你的后悔药。
3.2 数据源与预处理:ERA5 和 CMIP6 怎么接
ml_drought 需要两个层面的数据:训练模型用的历史再分析资料,以及验证预测用的气候模式输出。历史资料最常用的是欧洲中心的 ERA5,空间分辨率 0.25 度,时间上从 1940 年到现在,提供降水、温度、蒸散发、土壤湿度等变量。对于中国区域的干旱研究,有人用中国区域高分辨率数据集,但刚开始跑通流程时用 ERA5 就够了,因为下载接口成熟、格式统一。
import xarray as xr # 读取单个 NetCDF 文件,预览变量和坐标 ds = xr.open_dataset("era5_monthly_2010_2020.nc") print(ds.variables) # 常见变量名:tp(总降水)、t2m(2米气温)、pev(蒸散发) # 需要检查经纬度范围、时间频率,以及是否有缺失值 # 选择中国区域,并降采样到 0.5 度 ds_region = ds.sel(lat=slice(55, 15), lon=slice(70, 140)) ds_down = ds_region.resample(time="1M").mean()注意,ERA5 的降水变量是累积量,单位是米,逐月的值要累加起来。如果直接当瞬时量用,模型输入的数值范围会非常乱。气温和蒸散发则是月平均。把这些变量统一成月尺度后,再按多年气候态做标准化。常见做法是计算每个网格点每个月份的长期均值和标准差,然后用 (x - mean) / std 得到标准化异常。这里千万不能把全年数据混在一起算一个均值和标准差,否则季节信号会把异常信号淹没。
def compute_anomaly(ds, var_name): clim = ds[var_name].groupby("time.month").mean(dim="time") std = ds[var_name].groupby("time.month").std(dim="time") # 对每个月份减去对应的气候态均值 anomaly = (ds[var_name].groupby("time.month") - clim) # 再除以标准差,得到标准化异常 anomaly_std = anomaly / std return anomaly_std这段代码里有个隐藏边界:用全时段计算的气候态来标准化,在训练验证时会轻微泄漏未来信息。严格做法是用训练期数据计算气候态,再应用到验证期。但实际操作中,如果时间序列足够长(超过 30 年),这种泄漏影响很小。不过你要在论文里讲清楚的。
关于 CMIP6 气候模式数据,我的建议是先不用。ml_drought 如果要评估未来气候情景下的干旱风险,才需要 CMIP6 输出。日常训练用 ERA5 就够,验证模型泛化性时,再引入一套独立的再分析产品(比如日本 JRA-55)做外部验证。如果一开始就把 CMIP6 不同模式的输出混进训练集,模式之间的系统偏差会干扰模型学习真正的干旱机制,这是很多项目踩坑的地方。
3.3 用配置驱动实验:YAML 里的关键字段解读
研究型项目最容易出现“改一个参数改一次代码”的烂摊子。ml_drought 这种项目从设计上必须跑大量消融实验,所以很自然会用 YAML 配置来管理参数。我自己看一个深度学习项目,先不看模型代码,直接看配置文件,能看出作者对数据切分、损失函数、训练策略的思考。下面是一个典型的干旱预测项目配置片段。
data: era5_path: "./data/era5_monthly.nc" target_var: "spei3" input_vars: ["tp", "t2m", "pev"] start_year: 1980 valid_year: 2010 test_year: 2015 seq_len: 12 lead_time: 3 grid_size: [128, 256] normalize: "monthly_clim" model: arch: "unet" input_channels: 36 # 12个月 * 3个变量 hidden_channels: [64, 128, 256] dropout: 0.1 train: batch_size: 16 epochs: 50 learning_rate: 0.001 loss: "weighted_mse" drought_weight: 2.0 optimizer: "adam"这些字段看似简单,但每一个都对应一个决策边界。例如lead_time: 3是说预测未来 3 个月,如果改成 6,就要考虑输入窗口和标签窗口的距离拉长后样本数量减少的问题。input_channels: 36是 seq_len 和 input_vars 的乘积,在你设计模型时先不用手动计算,可以从配置读取。
import yaml with open("config/experiment.yaml", "r") as f: config = yaml.safe_load(f) seq_len = config["data"]["seq_len"] input_vars = config["data"]["input_vars"] model_input_channels = seq_len * len(input_vars) print("模型输入通道数 =", model_input_channels)这样做的好处是,实验记录里只要保存一份 YAML,就能完整还原这次跑的数据范围、模型结构、训练策略。我见过太多人用“exp1.py”这种文件来管理实验,最后自己都分不清哪个文件对应哪次结果。配置文件就是实验的身份证,这一点从第一天跑模型就要养成习惯。
4. 训练与评估:损失函数、评价指标与可视化
4.1 模型选型:U-Net、ConvLSTM 还是别的新花样
ml_drought 这类网格到网格的预测任务,最常见的主干是 U-Net 和 ConvLSTM,两者解决的问题略有不同。U-Net 通过编码器下采样捕获大尺度气候模式,再用解码器上采样还原到原始分辨率,非常适合像素级稠密预测。ConvLSTM 则把循环结构嵌在卷积里,更适合保留时间依赖。但实际跑下来,如果你只有 12 个月的输入序列,把时间维直接堆成通道,用 U-Net 也能达到相当好的效果,而且训练稳定得多。LSTM 类模型对序列长度和步长的敏感度高,调参成本也高。
import tensorflow as tf def conv_block(x, filters, skip=True): x = tf.keras.layers.Conv2D(filters, 3, padding="same", activation="relu")(x) x = tf.keras.layers.BatchNormalization()(x) x = tf.keras.layers.Conv2D(filters, 3, padding="same", activation="relu")(x) return x def build_unet(input_shape): inputs = tf.keras.Input(shape=input_shape) # 编码器部分 c1 = conv_block(inputs, 64) p1 = tf.keras.layers.MaxPooling2D(pool_size=2)(c1) c2 = conv_block(p1, 128) p2 = tf.keras.layers.MaxPooling2D(pool_size=2)(c2) c3 = conv_block(p2, 256) # 瓶颈层 bottleneck = conv_block(c3, 512) # 解码器部分 u1 = tf.keras.layers.UpSampling2D(size=2)(bottleneck) u1 = tf.keras.layers.Concatenate()([u1, c2]) c4 = conv_block(u1, 256) u2 = tf.keras.layers.UpSampling2D(size=2)(c4) u2 = tf.keras.layers.Concatenate()([u2, c1]) c5 = conv_block(u2, 128) # 输出层,回归到单通道 outputs = tf.keras.layers.Conv2D(1, 1, activation="linear")(c5) model = tf.keras.Model(inputs, outputs) return model这个 U-Net 是精简版,已经能跑通基本流程。但要说明几个设计选择和参数。首先是 BatchNormalization,在气候网格回归任务里,BN 能加速收敛,但对批量大小比较敏感。如果 batch size 只有 4 或 8,BN 的统计量会抖动,建议改用 LayerNormalization 或者干脆去掉。其次是上采样方式,这里用 UpSampling2D 加 Concatenate,最近几年流行的转置卷积和深度可分离上采样各有优劣。我的经验是,气候场本身比较平滑,上采样造成的高频伪影不多,简单 UpSampling 就够。最后是输出层用 linear,因为预测的 SPEI 是连续值,不是分类,千万不要在输出层加 sigmoid 或 softmax。
4.2 训练参数与调参优先级
训练一个干旱预测模型,不同参数的敏感度差异很大。优先级最高的是学习率。气候数据标准化后数值范围比较规整,初始学习率 0.001 通常能正常下降,但如果损失不降或者震荡,优先调低到 0.0005。第二个要调的是 batch size,它直接影响 GPU 显存和 BN 层表现。第三个才是网络宽度和深度。很多人一上来就把 U-Net 的 hidden_channels 加倍,结果训练时间翻倍,效果提升不显著。真正值得花时间的是数据侧:目标变量计算方式、预测提前期、区域划分。
# Keras 训练参数示意 model.compile( optimizer=tf.keras.optimizers.Adam(learning_rate=0.001), loss=weighted_mse, metrics=[mse, drought_hit_rate], ) callbacks = [ tf.keras.callbacks.ReduceLROnPlateau(patience=5, factor=0.5), tf.keras.callbacks.EarlyStopping(patience=10, restore_best_weights=True), ] history = model.fit( train_dataset, epochs=config["train"]["epochs"], validation_data=valid_dataset, callbacks=callbacks, )训练日志里除了 loss 和 mse,我建议加一个反映干旱异常区域命中率的自定义指标。干旱预测最怕的是模型把所有像素都预测成 0(正常),MSE 还能看过去,但业务上完全没用。比如定义 hit_rate:当真实 SPEI 小于 -0.5 时,预测也小于 -0.5 的像素占比。这样训练过程中就能直观看到模型是在学大尺度信号还是在小区域精细捕捉。
4.3 评估指标:不看总分,看分桶和极端事件
对干旱预测模型做评估,单独看 RMSE 或 MAE 不够。SPEI 在空间上大部分区域处于正常范围(-0.5 到 0.5),模型即使一直输出接近 0 的平均场,RMSE 也不会特别差。真正把一个预测系统拉出差距的,是对干旱区域的空间位置和强度峰值的把握。所以要按干旱强度分桶评估。常见的做法是把真实值和预测值都分成三档:湿润(SPEI > 0.5)、正常(-0.5 到 0.5)、干旱(SPEI < -0.5),然后看混淆矩阵。
import numpy as np from sklearn.metrics import confusion_matrix def evaluate_by_category(y_true_flat, y_pred_flat, thresholds=[-0.5, 0.5]): # 分桶:0=干旱,1=正常,2=湿润 y_true_cat = np.digitize(y_true_flat, thresholds, right=True) y_pred_cat = np.digitize(y_pred_flat, thresholds, right=True) cm = confusion_matrix(y_true_cat, y_pred_cat, labels=[0, 1, 2]) # 返回干旱类别命中率,即真实干旱中有多少被预测为干旱 drought_hit = cm[0, 0] / (cm[0, :].sum() + 1e-6) return cm, drought_hit从这张分桶表能看出很多问题。比如模型如果经常把干旱预测成正常,那就是漏报,对农业决策危害很大。如果经常把正常预测成干旱,那是虚报,会浪费救济资源。漏报和虚报在社会影响上不对等,所以评估时不要只看一个总分,最好分别计算干旱类别的召回率和命中率(精确率)。再深入一点,可以做空间尺度分析。干旱通常覆盖一定面积,单点命中率难以反映空间形态。把模型输出做空间平滑(比如高斯滤波)再与真实事件比较,更接近业务上的区域预警需求。很多项目在论文里声称精度高,但落到实际预警,紧急程度和空间范围才是决策者关心的。
5. 干旱预测避坑:数据泄漏、天气噪声与时空尺度
5.1 数据泄漏:滑动窗口切分惹的祸
现象:训练集损失和验证集损失都很好,但一到测试集或者换一个年份就崩溃。原因很可能是样本构造时滑动窗口的输入和目标有重叠。比如输入是 1 到 12 月,目标是 12 月当月的 SPEI,那模型只需要“抄”最后一个月就能拿到一个不错的分数。更隐蔽的是标准化泄漏:用全时间段计算的均值标准差去标准化输入,验证集的分布信息提前进入了训练。解决方法是构建样本时强制输入窗口结束月份早于目标窗口起始月份至少一个月,并在时间维度上按年份切块而不是随机采样。
# 错误示例:输入窗口结束与目标开始重叠 # inputs: 2000-01 到 2000-12,targets: 2000-12,模型直接可抄 # 正确:lead_time 设为 1 # inputs: 2000-01 到 2000-12,targets: 2001-01 到 2001-03在实践中,我最早做这个项目时就是没注意重叠,验证集 ACC 到了 0.9 以上,后来把时间轴画出来才发现目标月和输入最后一个月重合了。从那以后我养成了一个习惯:模型训练好之后,随机挑几个样本,把输入的最后几个月和输出目标画在同一张图上,肉眼检查有没有信息泄露的迹象。
5.2 天气尺度噪声:气候模型预测干旱是“信号弱、背景吵”
现象:模型预测的空间分布很破碎,像椒盐噪声一样,干旱中心不明显,甚至出现相邻网格预测值跳变。原因是输入的月平均降水本身就包含很多天气尺度随机性,而干旱是低频累积过程,模型在逐像素回归时容易拟合高频噪声。解决方法是显式约束输出空间平滑性。加一个空间平滑损失项,或者对目标变量做空间平滑后再训练。比如把 SPEI 场用高斯核做平滑,让模型学到平滑的目标分布。更工程化的方法是使用一个损失函数同时包含 MSE 和空间梯度惩罚。
def smoothness_loss(y_true, y_pred): # 计算预测场在水平和垂直方向的差分,惩罚剧烈跳变 dy = tf.abs(y_pred[:, 1:, :] - y_pred[:, :-1, :]) dx = tf.abs(y_pred[:, :, 1:] - y_pred[:, :, :-1]) return tf.reduce_mean(dy) + tf.reduce_mean(dx)这个平滑项在训练时可以给一个较小的权重,比如 0.1,否则模型会过度平滑,把干旱区域的峰值磨平。调这个参数的过程有些玄学,我的做法是先固定平滑权重,跑一轮看预测场的空间自相关,如果旱区中心不明显就调低。
5.3 样本不平衡:干旱事件少,模型很容易“躺平”
现象:模型预测结果几乎全为正常状态,干旱区域完全没有。原因和药物副作用数据不平衡一样,正常年份占绝大多数,模型只要保守预测,整体损失就能很低。解决这个问题有三个层面。一是数据层,把样本按目标年是否发生大型干旱进行加权采样,让干旱年份每个 epoch 出现的频率更高。二是损失函数层,用前面提到的加权 MSE 或自定义干旱权重。三是评价层,训练时用按干旱区域加权的指标作为早停依据,而不是原始 MSE。
# 加权采样示意:根据标签场是否大面积干旱决定是否重复采样 weights = [] for i in range(len(dataset)): label = dataset[i]["target"] drought_area = np.mean(label < -0.5) # 干旱面积越大,权重越高 weights.append(1.0 + 5.0 * max(drought_area - 0.1, 0.0))注意不要太激进,把权重拉到很高可能导致模型过拟合到几个严重干旱年份,失去泛化能力。一般干旱面积超过 30% 的样本权重到 5 倍就够。
5.4 区域和季节偏差:模型可能只学会了某个季节
现象:总体指标还可以,但逐季节分析发现夏季干旱预测能力显著弱于冬季。原因在于气候模式对不同季节可预报性的物理基础不同。热带以外地区冬季的大尺度环流更稳定,模型的输入特征信息量大,所以效果好。解决方法是把评估结果按月份或季节拆分,并在训练时考虑季节编码。常见做法是给输入张量增加一个“月份正弦/余弦编码”通道,让模型知道当前预测的起点是哪个月份,而不是让它自己从数据里猜。
# 给输入拼接月份编码 import numpy as np month = 5 # 假设输入最后一个月是5月 sin_enc = np.sin(2 * np.pi * month / 12) cos_enc = np.cos(2 * np.pi * month / 12) # 将两个标量编码成与空间网格同形状的通道,拼到输入最后这个做法能有效防止模型把不同季节的预测混为一谈。特别是做提前 3 个月的预测时,1 月出发的预测和 7 月出发的预测,环境背景完全不同。
6. 进阶:把预测结果变成可解释的区域干旱预警信号
6.1 归因图与类激活图:让模型告诉你它看了哪块云
模型训练好只是第一步,要用它做决策,必须知道预测的干旱中心是由哪个变量、哪个区域驱动的。对 U-Net 这类卷积模型,最容易上手的是梯度加权类激活图(Grad-CAM)。它用输出通道对最后一层特征图的梯度作为权重,加权求和得到空间热力图,能直观看到模型重点关注的输入区域。
def grad_cam(model, sample_input, layer_name="conv2d_7"): grad_model = tf.keras.models.Model( inputs=model.input, outputs=[model.get_layer(layer_name).output, model.output] ) with tf.GradientTape() as tape: feature_map, prediction = grad_model(sample_input) # 取预测场中干旱区域平均强度作为归因目标 drought_score = tf.reduce_mean(prediction[0, ...]) grads = tape.gradient(drought_score, feature_map) # 对梯度的空间维求均值,得到每个通道的重要性 channel_weights = tf.reduce_mean(grads, axis=(1, 2)) cam = tf.reduce_sum(tf.multiply(channel_weights, feature_map[0]), axis=2) cam = tf.maximum(cam, 0) # 只保留正向贡献 cam = (cam - tf.reduce_min(cam)) / (tf.reduce_max(cam) - tf.reduce_min(cam) + 1e-6) return cam.numpy()归因图不要只看单个样本,要按事件合成。比如把过去 20 年里面最大的 5 个干旱事件样本分别做 CAM,再平均,得到的平均归因图更有代表性。如果平均图中干旱中心附近总是对应着当时的暖异常区域,那说明模型学到的是温度偏高加剧蒸散发的机制,这在物理上是合理的,可信度就高。
6.2 从预测输出到干旱等级:业务侧的分类决策
模型输出的 SPEI 是连续值,但实际预警需要的是等级。常见的四分法是:轻旱(SPEI 在 -1 到 -0.5)、中旱(-1.5 到 -1)、重旱(-2 到 -1.5)、特旱(小于 -2)。从连续输出映射到等级时,不要直接用固定阈值,因为模型的预测有不确定性。更稳健的办法是同时预测均值和方差,或者用多个训练种子(比如 5 个不同随机种子)得到集成结果,然后看集成成员落在某个等级的比例。
# 多成员集成预测,输出每个格点落入干旱等级的概率 def ensamble_to_alert(ensemble_preds, thresholds=[-2, -1.5, -1, -0.5]): # ensemble_preds: [5, H, W] 多个模型对同一时次的预测 probs = [] for t in thresholds: prob = np.mean(ensemble_preds < t, axis=0) probs.append(prob) # 综合多个等级概率,输出最高等级 alert = np.argmax(np.stack(probs), axis=0) return alert这个集成方法虽然简单,但对业务报告非常实用。干旱预警最怕的是“确定性的单点输出”给决策者带来虚假信心。用集成概率表达,每一个格点都能给出“未来三个月进入中旱的可能性 80%”之类的话。我现在每跑一个实验,都会额外保存每个测试样本每个格点的 ensemble 输出,这样一旦上级或客户问“你这个置信度怎么来的”,能立刻给出一张空间概率图,而不是只甩一个数值。
做机器学习预测干旱,有时候难的不是模型有多先进,而是能不能把模型输出转化为经得起追问的决策依据。我自己的习惯是训练后第一步不是看 loss,而是随机抽三个月的时间切片,把输入、预测、真实值和归因图画在一张图上,先肉眼看有没有明显不合理的错位,再调参数。数据泄漏、季节偏差、空间破碎这些问题,单靠指标很难发现,但图像一眼就能看出问题。希望这些经验能帮你少走几步弯路。
本文还有配套的精品资源,点击获取