简介这份PDF文献面向气象预报、雷达外推与短临预报方向的研究人员和业务人员聚焦强对流天气临近预报中光流法及其与机器学习结合的应用。资源为单文件PDF压缩包约1.51MB内容完整便于直接阅读与引用。文献系统梳理了光流法追踪雷达回波运动矢量场的原理说明其相较传统交叉相关法在强对流追踪中的优势并引入半拉格朗日外推方案保持回波旋转效应提升0至2小时外推精度。在此基础上作者采用ConvLSTM空间深度学习模型结合Z-R关系开展1小时降水预测结果表明其空漏报率低于SWAN业务系统QPF结果。文中还结合2017年福建强对流个例进行对比分析展示了从方法原理到实验验证的完整思路。目前已有161人学习适合希望了解光流法、机器学习与雷达回波外推结合路径的读者参考。1. 光流法做临近预报从两张雷达图到未来两小时降水临近预报要回答的问题很具体未来 0 到 2 小时某块区域会不会下雨、雨多大、往哪走。这个时间尺度太短数值天气预报模式跑一次都不够所以业内长期靠外推法撑场面。外推法的核心假设是降水系统在短时间内形态变化不大主要做平移运动于是问题就变成了——怎么从前后两帧雷达回波图里把运动场估出来。光流法就是干这个的经典工具它逐像素计算位移矢量得到一张稠密的运动场再把这张运动场外推到未来时刻就得到了预报场。这套方案适合谁做短临预报业务的气象工程师、做雷达外推算法研究的同学、以及想把机器学习引入气象外推但不知道从哪下手的开发者。纯光流法能跑出一个基线但它有天花板回波的生消、合并、分裂它处理不了因为光流只建模平移。于是近几年的主流思路变成——用光流法提供运动场先验用 CNN、ConvLSTM 这类模型去学形态演变两者结合。这篇就按这个脉络从光流法的最小实现讲到机器学习怎么接进去把参数、坑和验证方法都摊开。2. 光流法估运动场Horn-Schunck 与 Lucas-Kanade 怎么选2.1 光流的两条约束到底约束了什么光流法的出发点是亮度恒常假设同一个降水粒子在相邻两帧之间移动它的回波强度不变。写成数学式就是 I(x, y, t) I(xu, yv, t1)对时间做一阶泰勒展开得到 I_x·u I_y·v I_t 0。这一个方程有两个未知数 u 和 v欠定所以必须再加约束才能解。Horn-Schunck 加的是全局平滑约束假设整个场的运动是平滑的最小化 (I_x·u I_y·v I_t)² α²(|∇u|² |∇v|²)。它给出稠密光流每个像素都有矢量适合雷达回波这种连续场。Lucas-Kanade 加的是局部窗口内运动一致的假设在一个小窗口里用最小二乘解给出稀疏但更稳的光流适合特征点跟踪。对临近预报来说你要的是一张覆盖整个雷达图的运动场所以 Horn-Schunck 或者它的变体是更自然的选择。但 Horn-Schunck 对噪声敏感雷达图里的地物杂波、超折射回波会污染运动场估计这是第一个要处理的坑。2.2 用 OpenCV 跑通最小可用的光流外推先不急着上机器学习把纯光流基线跑出来后面才有对比的锚点。下面这段代码读两帧雷达回波图算稠密光流然后做半拉格朗日外推。import cv2 import numpy as np # 读两帧雷达回波灰度图0-255 对应反射率因子 frame1 cv2.imread(radar_t0.png, cv2.IMREAD_GRAYSCALE) frame2 cv2.imread(radar_t1.png, cv2.IMREAD_GRAYSCALE) # Farneback 稠密光流参数含义见下方说明 flow cv2.calcOpticalFlowFarneback( frame1, frame2, None, pyr_scale0.5, # 金字塔缩放0.5 表示每层缩小一半 levels3, # 金字塔层数3 层能抓大位移 winsize15, # 平均窗口越大越平滑但越糊 iterations3, # 每层迭代次数 poly_n5, # 像素邻域多项式展开大小 poly_sigma1.2, # 高斯标准差 flags0 ) # 半拉格朗日外推对每个输出像素沿光流反向找源像素 def advect(frame, flow, steps1): h, w frame.shape yy, xx np.mgrid[0:h, 0:w].astype(np.float32) # 反向追踪所以减号 src_x xx - flow[..., 0] * steps src_y yy - flow[..., 1] * steps # 双线性插值采样 out cv2.remap(frame, src_x, src_y, cv2.INTER_LINEAR, borderModecv2.BORDER_REPLICATE) return out # 外推未来 6 帧每帧间隔 10 分钟 pred frame2.copy() for i in range(1, 7): pred advect(pred, flow, steps1) cv2.imwrite(fpred_t{i}.png, pred)这段代码的逻辑是Farneback 算法用多项式展开近似每个像素邻域的亮度变化从而估计稠密光流advect 函数做半拉格朗日外推对每个输出像素沿光流反向追踪到源位置用双线性插值取值。参数里最需要调的是 winsize 和 levels。winsize 太小运动场噪声大回波边缘会出现碎裂winsize 太大运动场过度平滑强对流单体的旋转结构会被抹掉。levels 决定能捕捉多大位移雷达图相邻帧间隔 6 到 10 分钟强回波移动可能超过 20 像素levels 至少设 3。提示Farneback 的输入最好是归一化后的反射率因子不要直接喂原始 dBZ 值。原始值动态范围大弱回波和强回波的梯度差异悬殊光流会被强回波主导弱降水区的运动场估计不准。2.3 光流场的后处理三个必做的滤波原始光流场不能直接用至少要做三步后处理。第一步是杂波掩膜把地物杂波和超折射区域抠掉这些区域的光流是伪运动。第二步是中值滤波去掉孤立的异常矢量用 cv2.medianBlur 对 u、v 分量分别做 5×5 中值。第三步是散度约束真实降水场的运动场散度不应该太大如果某个区域散度超过阈值说明光流估计失败可以用邻域均值替换。这三步做完运动场的视觉质量会明显提升外推的回波不会出现大面积的撕裂和漂移。但即便如此纯光流外推的预报时效通常也就 30 到 60 分钟再往后回波强度衰减和形态变化的误差就盖过平移误差了。这就是为什么要把机器学习接进来。3. 机器学习接进光流CNN 和 ConvLSTM 各管什么3.1 为什么不让 CNN 直接从两帧预测未来帧一个自然的想法是既然 CNN 这么强为什么不直接把前后两帧堆成通道让 CNN 回归未来帧这样做在实验里能跑出结果但泛化性差。原因是 CNN 的卷积核是空间不变的它学到的是局部纹理到局部纹理的映射而降水系统的演变是高度依赖运动方向的。同一个回波形态往东移动和往东北移动未来演变完全不同但 CNN 在卷积层面区分不了这个差异除非你把运动信息显式喂给它。所以更稳的架构是双分支一个分支用光流法或可变形卷积估计运动场另一个分支用 CNN 提取回波形态特征然后在特征层面做运动补偿。ConvLSTM 则适合处理时序它把 LSTM 的门控机制换成卷积操作能在保留空间结构的同时建模时间依赖。对临近预报ConvLSTM 的输入通常是最近 5 到 10 帧输出未来 6 到 12 帧。3.2 用 PyTorch 搭一个光流引导的 ConvLSTM 外推网络下面是一个最小可跑的架构光流分支用预计算的 Farneback 光流作为输入形态分支用 CNN两者融合后送进 ConvLSTM 做时序外推。import torch import torch.nn as nn class ConvLSTMCell(nn.Module): def __init__(self, in_ch, hid_ch, kernel_size3): super().__init__() padding kernel_size // 2 self.conv nn.Conv2d(in_ch hid_ch, 4 * hid_ch, kernel_size, paddingpadding) self.hid_ch hid_ch def forward(self, x, h, c): combined torch.cat([x, h], dim1) gates self.conv(combined) i, f, o, g torch.split(gates, self.hid_ch, dim1) i, f, o, g torch.sigmoid(i), torch.sigmoid(f), \ torch.sigmoid(o), torch.tanh(g) c_next f * c i * g h_next o * torch.tanh(c_next) return h_next, c_next class FlowGuidedNowcaster(nn.Module): def __init__(self, in_ch2, hid_ch64): super().__init__() # 形态分支提取回波纹理 self.form_cnn nn.Sequential( nn.Conv2d(in_ch, 32, 3, padding1), nn.ReLU(), nn.Conv2d(32, 32, 3, padding1), nn.ReLU() ) # 光流分支2 通道 u,v self.flow_cnn nn.Sequential( nn.Conv2d(2, 16, 3, padding1), nn.ReLU(), nn.Conv2d(16, 16, 3, padding1), nn.ReLU() ) self.fuse nn.Conv2d(48, hid_ch, 1) self.convlstm ConvLSTMCell(hid_ch, hid_ch) self.head nn.Conv2d(hid_ch, 1, 1) def forward(self, frames, flows): # frames: (B, T, C, H, W), flows: (B, T, 2, H, W) B, T, _, H, W frames.shape h torch.zeros(B, 64, H, W, deviceframes.device) c torch.zeros(B, 64, H, W, deviceframes.device) outputs [] for t in range(T): f_form self.form_cnn(frames[:, t]) f_flow self.flow_cnn(flows[:, t]) x self.fuse(torch.cat([f_form, f_flow], dim1)) h, c self.convlstm(x, h, c) outputs.append(self.head(h)) return torch.stack(outputs, dim1)这个架构的关键设计点光流分支和形态分支在通道维拼接后用 1×1 卷积融合1×1 卷积的作用是学习两个分支的加权组合而不是简单相加。ConvLSTM 的隐藏状态在时间步之间传递所以它能记住前几帧的运动趋势。输出头用 1×1 卷积把隐藏状态映射回单通道回波强度。参数上hid_ch 设 64 是精度和显存的折中如果你做 512×512 的雷达图batch size 8显存大概占 6 到 8 GB。如果显存不够把 hid_ch 降到 32但预报的细节会损失。输入帧数 T 建议 5 到 10太少学不到运动趋势太多训练慢且容易过拟合。3.3 训练时的损失函数MSE 不够要加梯度损失临近预报的评估不只看像素误差还看回波结构。只用 MSE 训练模型倾向于输出模糊的平均场强回波中心会被抹平。常见做法是 MSE 加梯度损失梯度损失用 Sobel 算子算预测场和真值场的梯度差逼模型保留边缘。def gradient_loss(pred, target): # Sobel 核 kx torch.tensor([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtypetorch.float32).view(1, 1, 3, 3) ky kx.transpose(2, 3) kx kx.to(pred.device) ky ky.to(pred.device) pred_gx torch.abs(torch.nn.functional.conv2d(pred, kx, padding1)) pred_gy torch.abs(torch.nn.functional.conv2d(pred, ky, padding1)) tgt_gx torch.abs(torch.nn.functional.conv2d(target, kx, padding1)) tgt_gy torch.abs(torch.nn.functional.conv2d(target, ky, padding1)) return torch.mean(torch.abs(pred_gx - tgt_gx)) \ torch.mean(torch.abs(pred_gy - tgt_gy)) # 总损失 loss mse_loss(pred, target) 0.5 * gradient_loss(pred, target)梯度损失的权重 0.5 是经验值调大到 1.0 会让边缘更锐利但可能引入噪声调小到 0.1 则退化成接近纯 MSE。训练时先用 MSE 预热几个 epoch再引入梯度损失收敛更稳。4. 数据、评估与避坑临近预报落地时最容易翻车的地方4.1 雷达数据预处理的四个关键步骤原始雷达基数据不能直接喂模型标准流程是解码、坐标变换、质量控制、归一化。解码用 pyart 或 wradlib 读 CINRAD 或 NEXRAD 格式。坐标变换把极坐标的方位角-距离网格插值到笛卡尔网格常用最近邻或双线性。质量控制包括地物杂波抑制和超折射剔除pyart 的 dealias 和 filter 能处理一部分。归一化把反射率因子从 -10 到 70 dBZ 线性映射到 0 到 1。注意坐标变换的插值方法会影响光流估计。最近邻插值会产生块状伪影光流在这些边界上会算出虚假的大位移。建议用双线性插值虽然会轻微平滑回波但对光流更友好。4.2 评估指标CSI、FSS 和光流场的端点误差临近预报的评估不能只看 RMSE。业务上常用 CSI 和 FSS。CSI 针对某个阈值算命中率公式是 hits / (hits misses false_alarms)。FSS 是邻域空间检验对高分辨率预报更公平因为它允许小位移误差。光流场本身的评估用端点误差就是估计的 (u, v) 和真实 (u, v) 的欧氏距离但真实光流很难获取通常用合成数据或者人工标注的少数样本做验证。指标适用场景阈值建议注意点CSI单阈值命中评估20/30/40 dBZ对位移误差敏感FSS高分辨率邻域评估20/30/40 dBZ邻域窗口选 5×5 或 9×9RMSE整体强度误差全量程会被弱回波主导端点误差光流场质量像素/帧需要真值难获取4.3 避坑清单五条血泪经验现象光流场在回波边缘出现大量异常矢量。原因回波边缘梯度大亮度恒常假设不成立Farneback 的多项式展开在边缘处失效。解决对光流场做中值滤波或者用加权中值滤波权重用回波强度弱回波区的光流矢量降权。现象模型在训练集上 CSI 很高测试集上掉 20 个点。原因雷达数据有强烈的季节和地域差异训练集如果只覆盖夏季对流测试集遇到层状云降水就翻车。解决训练集要覆盖不同季节、不同降水类型至少包含对流性、层状云和混合性三类。现象ConvLSTM 预测到第 6 帧以后回波强度整体偏弱。原因MSE 损失对弱回波的惩罚小模型倾向于输出保守的平均值。解决加梯度损失或者对强回波区域加权权重可以设为回波强度的平方。现象光流外推的回波整体偏移方向对但移动速度偏慢。原因Farneback 的 winsize 太大运动场被过度平滑位移矢量被低估。解决减小 winsize 到 9 或 11同时增加 levels 到 4让金字塔捕捉更大位移。现象训练 loss 震荡不收敛。原因雷达数据的动态范围大未归一化或者归一化不当。解决确认输入在 0 到 1 之间用 BatchNorm 或者 GroupNorm 稳定训练学习率从 1e-4 开始用 cosine 退火。5. 把预报时效从 60 分钟推到 120 分钟一个可验证的改进技巧纯光流外推的时效瓶颈在 60 分钟左右原因是它假设运动场不变。但真实降水系统的运动场是随时间演变的强对流单体会发展、合并、分裂。要把时效推到 120 分钟核心思路是让运动场也变成可预测的。我一般会用一个轻量的运动场预测头输入最近几帧的光流场输出未来几帧的光流场然后用预测的光流场做外推而不是用最后一帧的光流场反复外推。具体做法在 FlowGuidedNowcaster 的基础上加一个辅助任务让 ConvLSTM 的隐藏状态同时预测下一帧的光流。损失函数里加一项光流预测的 MSE权重设 0.3。这样模型在学回波演变的同时也被迫学运动场的演变。训练完后推理时用预测的光流场做半拉格朗日外推而不是复用第一帧的光流。验证方法上我习惯做分时效评估分别算 0-30 分钟、30-60 分钟、60-90 分钟、90-120 分钟的 CSI。如果 60 分钟以后的 CSI 比纯光流基线高 5 个点以上说明运动场预测头起了作用。如果没提升检查光流预测的端点误差可能是辅助任务的权重太小或者 ConvLSTM 的隐藏状态容量不够。还有一个技巧是多尺度融合。雷达回波里既有大尺度的层状云也有小尺度的对流单体。单一尺度的光流场要么抓大尺度丢小尺度要么反过来。我一般会做两个尺度的光流估计一个用大 winsize 抓大尺度运动一个小 winsize 抓小尺度然后在预报时按回波强度加权融合。强回波区用小尺度光流弱回波区用大尺度光流。这个技巧在春季混合性降水里提升明显CSI 能涨 3 到 5 个点。最后说一个我踩过的坑不要用未来帧的信息去算光流。听起来是废话但在做数据管道时很容易犯。比如你用 t 和 t1 帧算光流然后外推 t2这没问题。但如果你用 t1 和 t2 算光流去外推 t3在实时业务里 t2 还没到这就是数据泄漏。训练时如果管道写错验证指标会虚高上线就翻车。我的习惯是在数据加载器里显式检查时间戳确保输入帧的时间严格早于预测帧。这套方案值不值得做如果你只需要 30 分钟内的外推纯光流加后处理就够了成本低、可解释。如果你要做 60 到 120 分钟的预报并且有历史雷达数据积累光流加机器学习的混合方案是目前性价比最高的路径。算力上一张 12 GB 显存的卡就能训数据上至少需要 1 到 2 年的雷达体扫数据。先从单站点跑通再考虑多站点融合。希望帮到你。本文还有配套的精品资源点击获取
