简介：这份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(x+u, y+v, t+1)，对时间做一阶泰勒展开，得到 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_scale=0.5, # 金字塔缩放，0.5 表示每层缩小一半 levels=3, # 金字塔层数，3 层能抓大位移 winsize=15, # 平均窗口，越大越平滑但越糊 iterations=3, # 每层迭代次数 poly_n=5, # 像素邻域多项式展开大小 poly_sigma=1.2, # 高斯标准差 flags=0 ) # 半拉格朗日外推：对每个输出像素，沿光流反向找源像素 def advect(frame, flow, steps=1): 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, borderMode=cv2.BORDER_REPLICATE) return out # 外推未来 6 帧，每帧间隔 10 分钟 pred = frame2.copy() for i in range(1, 7): pred = advect(pred, flow, steps=1) cv2.imwrite(f'pred_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_size=3): super().__init__() padding = kernel_size // 2 self.conv = nn.Conv2d(in_ch + hid_ch, 4 * hid_ch, kernel_size, padding=padding) self.hid_ch = hid_ch def forward(self, x, h, c): combined = torch.cat([x, h], dim=1) gates = self.conv(combined) i, f, o, g = torch.split(gates, self.hid_ch, dim=1) 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_ch=2, hid_ch=64): super().__init__() # 形态分支：提取回波纹理 self.form_cnn = nn.Sequential( nn.Conv2d(in_ch, 32, 3, padding=1), nn.ReLU(), nn.Conv2d(32, 32, 3, padding=1), nn.ReLU() ) # 光流分支：2 通道 u,v self.flow_cnn = nn.Sequential( nn.Conv2d(2, 16, 3, padding=1), nn.ReLU(), nn.Conv2d(16, 16, 3, padding=1), 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, device=frames.device) c = torch.zeros(B, 64, H, W, device=frames.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], dim=1)) h, c = self.convlstm(x, h, c) outputs.append(self.head(h)) return torch.stack(outputs, dim=1)这个架构的关键设计点：光流分支和形态分支在通道维拼接后用 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]], dtype=torch.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, padding=1)) pred_gy = torch.abs(torch.nn.functional.conv2d(pred, ky, padding=1)) tgt_gx = torch.abs(torch.nn.functional.conv2d(target, kx, padding=1)) tgt_gy = torch.abs(torch.nn.functional.conv2d(target, ky, padding=1)) 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×9 |
| RMSE | 整体强度误差 | 全量程 | 会被弱回波主导 |
| 端点误差 | 光流场质量 | 像素/帧 | 需要真值，难获取 |
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 和 t+1 帧算光流，然后外推 t+2，这没问题。但如果你用 t+1 和 t+2 算光流去外推 t+3，在实时业务里 t+2 还没到，这就是数据泄漏。训练时如果管道写错，验证指标会虚高，上线就翻车。我的习惯是在数据加载器里显式检查时间戳，确保输入帧的时间严格早于预测帧。
这套方案值不值得做？如果你只需要 30 分钟内的外推，纯光流加后处理就够了，成本低、可解释。如果你要做 60 到 120 分钟的预报，并且有历史雷达数据积累，光流加机器学习的混合方案是目前性价比最高的路径。算力上，一张 12 GB 显存的卡就能训，数据上，至少需要 1 到 2 年的雷达体扫数据。先从单站点跑通，再考虑多站点融合。希望帮到你。
本文还有配套的精品资源，点击获取
