基于卷积光流生成网络的雷达回波图像预测算法研究

张毓杉，周 晓

(武汉理工大学 机电工程学院，湖北 武汉 430070)

天气预报分为长期预报，中期预报，短期预报，短时预报和临近预报5种[1]。本文研究的临近预报是指的2 h以内的预报。临近预报主要是通过雷达回波图像的外推技术来实现的，对于雷达回波历史图像，经过算法分析处理后，预测生成未来时刻的雷达回波图像。雷达回波图像是灰度图像，其灰度值大小反映了当前的降雨情况[2]，图像上灰度值越大的区域，表示此处观测到的回波值越大，也说明了此处降雨程度越严重。1978年Rinehart[3]提出了使用交叉相关法来跟踪雷达回波。1993年，Dixon[4]提出了基于雷达回波数据对云团进行识别、追踪和预报的TITAN(thunderstorm identification tracking analysis and nowcasting)算法。后来Johnson[5]提出了SCIT(storm cell identification and tracking)算法，对TITAN算法进行了改进。因为光流法的兴起，Bowler[6-7]等提出利用光流算法来外推雷达回波图像来实现降雨预测的目的。

笔者以雷达回波降雨预测为研究目标，结合传统的光流法[8]和卷积神经网络提出了基于卷积神经网络的光流生成网络(convolutional neural networks-flow, CNN-Flow)生成雷达回波图像间的光流场。最后再利用光流场和半拉格朗日算法预测出未来时刻的图像序列。

1 卷积光流生成网络

卷积神经网络(convolutional neural networks, CNN)拥有各种各样的结构，相对于全连接神经网络，由于其权重共享的特性，大大减少了训练时的参数数量，因此选择CNN作为基础结构来生成光流场。假设雷达回波图像的宽与长分别为w和h，且由于雷达回波图像是灰度图像，其通道数为1，则两个输入图像的尺寸均为h×w×1。而光流场中每个点代表了图像中每个像素点的横纵移动速度分量，因此网络输出的光流场尺寸为h×w×2，网络的总体工作流程如图1所示。

图1 CNN-Flow网络工作流程

从图1可知，输入网络的两帧图像首先进行通道融合，在通道维度上进行叠加，再输入网络，经过编码解码器的上下采样U型结构后，最后生成对应的光流场。

1.1 网络结构

该网络由编码器和解码器组成。编码器的作用不仅是为了减轻训练压力，也压缩了图像的尺寸，使得卷积层可以获得更大的感受野，提取更加深层的特征。解码器的作用是恢复图像的尺寸，因为雷达回波各个部分移动的速度可能不一样，有的区域位移过大，有的区域位移过小，为了能够同时处理不同移动幅度的回波，网络在解码过程的同时在不同尺寸的图像上利用步长为1的卷积层提取不同尺度的光流场，然后将不同尺度的光流场同时与经过缩放的真实光流场进行损失值的计算，最后网络输出与原图像尺寸一样大小的光流场，CNN-Flow的具体结构如图2所示。

图2 CNN-Flow网络结构

从图2可知，由于网络层数过少会使得提取的特征不足，层数过多对训练结果也不会有太大的提升，反而提高了训练难度。因此该网络一共设有6层编码层和6层解码层，其中Conv代表卷积层，Deconv代表反卷积层，Concat代表连接层(用于将多个特征图按通道维度进行融合)，flow1-6代表不同尺寸大小的光流场，这些光流场均是在解码过程中利用卷积层对不同分辨率的特征图进行提取而生成的，由于光流场的通道数是2(代表两个速度分量)，因此提取光流场的卷积层的卷积核数量均为2。

CNN-Flow的编码器由10个卷积层叠加而成，每一层的作用是向下一层或者同层解码层传递特征，解码器由6个反卷积层及6个连接层组成，解码器的每一层输入的特征都由上一层解码层输出的特征以及上一层解码层提取的光流场和同一层的编码层输出的特征通过跳层连接融合而成的。网络训练过程中将所有提取的不同尺寸的光流场与相对应尺寸的真实光流场进行损失函数计算，然后将所有损失值进行加权和。由于LeakyReLU[9]激活函数解决了神经元死亡的问题，因此网络选取了LeakyReLU作为激活函数。

1.2 多权重组合损失函数

由于光流场的每一个格点均由两个速度分量u和v组成，为了使得训练时输出的光流场接近真实的光流场，使用如式(1)所示的指标EPE(endpoint error,端点误差)作为损失函数的基本计算公式。

(1)

网络中一共输出了6个不同尺度的光流场作为输出，最后需要计算6个不同尺度的EPE，然后获得EPE加权损失函数。尺寸越大的光流场所计算出的EPE的权重越大，加权损失函数计算公式如式(2)所示。

Loss=a1×EPE1+a2×EPE2+a3

×EPE3+a4×EPE4+a5*EPE5+

a6×EPE6

(2)

式中：EPE1为尺度最小的光流场；EPE6为尺度最大的光流场。经过多次尝试比较，最后权重设置如下：a1=0.004，a2=0.01，a3=0.02，a4=0.06，a5=0.1，a6=0.31。

2 雷达回波图像预测

通过CNN光流网络获取光流场之后，还需要利用光流场以及计算出光流场的原始两帧图像进行图像外推计算。此处外推计算使用的是半拉格朗日外推法，一维的半拉格朗日算法为：现假设某一粒子的F值由其坐标x和时间t所决定，即F(x,t)，则半拉格朗日的基本方程如式(3)和式(4)所示：

(3)

(4)

式(3)和式(4)表明F在沿着粒子运动轨迹时是一个恒定不变的值。U(x,t)代表着当前点的速度。

图3为半拉格朗日外推方法的大致模型。图3表示了粒子的移动情况，AC这条粗线表示了粒子从t-Δt时刻到t+Δt时刻的真实移动路径；aC这条细线表示了粒子从t-Δt时刻到t+Δt时刻的近似移动路径。A和a作为粒子的起始点；B是粒子的中间点；C是粒子的终点。一般终点在图像格点上，而起点和中间点都不会刚好在网格点上。Δt是粒子一次移动的时间间隔；2α是粒子从t-Δt时刻到t+Δt时刻的移动距离。

图3 半拉格朗日外推方法的大致模型

由图3可以得到式(5)和式(6)：

(5)

α=ΔtU(x-α,t)

(6)

式(6)中α值可由式(7)进行迭代求出。

(7)

这样就能在t-Δt时刻找到t+Δt时刻对应的F值。

对于雷达回波的外推问题，因为是图像外推，所以需要将一维的半拉格朗日扩展成二维，现假设粒子(图像中的像素点)的F值(即灰度值)由其坐标x和y以及时间t决定，二维半拉格朗日的方程变为如下形式：

(8)

(9)

(10)

式(8)～式(10)中U为像素点x方向的速度；V为像素点y方向上的速度。于是像素点的运动轨迹可由式(11)～式(13)决定：

(11)

α=ΔtU(x-α,y-β,t)

(12)

β=ΔtV(x-α,y-β,t)

(13)

式(12)和式(13)依然可用迭代方式求出。由于起始点和中间点不可能会刚好落在图像整数坐标上，因此F(x-2α,y-2β,t-Δt)，U(x-α,y-β,t)，V(x-α,y-β,t)需要通过图像插值的方式求出。假设(i,j)，(i,j+1)，(i+1,j)，(i+1,j+1)4个点分别为邻近的4个格点，而起始点和中间点会落在(i+u,j+v)处，u和v均为浮点数，可采用双线性插值法，即通过式(14)求出。

f(i+u,j+v)=(1-u)×(1-v)×f(i,j)+

u×(1-v)×f(i+1,j)+(1-u)×v×

f(i,j+1)+u×v×f(i+1,j+1)

(14)

通过上述方法就能利用光流场以及原始的两帧图像推出第三帧的图像。由于CNN-Flow结合半拉格朗日外推算法后只能通过输入两帧图像预测出第三帧图像，笔者采用递归调用的方式生成10帧预测图像，通过输入最新生成的图像以及上一帧图像，重新计算光流场并推算下一帧图像。

3.1 数据集

采用的数据集来自于深圳市气象局提供的粤港澳大湾区的雷达回波图像数据集(SRAD)，SRAD共包含20 000个样本，其中每个雷达数据个案样本覆盖时长为4小时，样本中每帧图像采集时间间隔6分钟，图像长和宽均为256(即大约覆盖255 km×255 km的区域)。

设计的CNN-Flow需要两帧连续的雷达回波图像作为输入，而输出是一帧两通道的光流场。网络要求的数据集的每一个样本应该包含两张雷达回波图像以及一个光流场文件，但是由于目前还没有能够直接检测到光流的传感器，无法建立标准的雷达回波图像的光流数据集，因此笔者首先将SRAD进行滑动窗口切分，保持每个样本中包含有两帧连续的雷达回波图像，再对其打乱并随机筛选切分成了8 000个训练样本及2 000个测试样本，最后用Farneback算法处理每一个样本的两帧图像，生成对应的光流场文件加入对应的样本中，然后将新的数据集称为FSRAD(flow SRAD)。

3.2 评价指标

实验采取的是降水临近预报评估方法，其中包含有CSI(critical success index,临界成功指数)，FAR(false alarm rate,误报率)，POD(probebility of detection,检测概率)[10]3个指标来评判预报的准确性。首先通过设定一个降雨阈值，将预测图像与真实图像中在阈值内的格点设置为下雨点(即标记为1)，阈值外的格点设置为无雨点(即标记为0)，然后按像素遍历预测图像及真实图像，当预测值为1，真实值为1时记为一个hits；当预测值为1，真实值为0时记为一个alarms；当预测值为0，真实值为1时记为一个misses。最后统计hits，alarms以及misses的个数。式(15)为3个指标的计算公式。

(15)

式中：nhits为hits点的个数；nalarms为alarms点的个数；nmisses为misses点的个数。其中CSI与POD的值越高代表预测越准确，FAR的值越低代表预测误报率越低。

3.3 模型训练参数设置

网络训练超参数设置如表1所示。

表1 网络的超参数

3.4 实验结果分析

由于EPE指标可用于验证生成的光流场及真实光流场之间的差距，且EPE值越小代表生成的光流场越接近真实值，因此该网络在训练过程中记录了迭代过程中在FSRAD训练集及测试集上的平均EPE数值变化，结果如图4所示。

图4 EPE变化曲线

从图4可知，迭代了160次之后，训练集和测试集的平均EPE已经开始收敛，并且训练集的EPE收敛于0.9201，测试集的EPE收敛于0.724 4。表2为CNN-Flow和传统光流法在FSRAD测试集上生成一帧光流场所需的平均时间对比。

表2 生成光流场的平均时间消耗

从表2可知，CNN-Flow速度明显快于传统的光流法，因此基于深度学习的光流法在运行效率上有了较大的提升。为了验证CNN-Flow在预测图像上的降雨预测精度，将生成的光流场结合半拉格朗日外推法，通过两帧雷达回波图像预测未来10帧的图像，实验中将回波阈值设为30 dB，以此阈值在CNN-Flow和传统光流法预测的10帧图像以及实际的10帧图像上用3个降雨指标进行对比分析。表3为两种算法在测试集上所计算的3种指标的均值，可以看出CNN-Flow在CSI指标和FAR指标上均略差于传统光流法，在POD指标上强于传统光流法，结合表2和表3分析可得，CNN-Flow网络在总体精度上和传统光流法差不多的情况下，运行速度远远超过了传统光流法。

表3 算法在测试集上的平均指标

为了测试算法在个例样本上的预测效果，选取了某个时刻的个例样本绘制了两种算法的CSI，FAR，POD指标随着预测时间增加的变化曲线，如图5所示，其中横坐标代表输入前两帧图像，预测3至12帧图像。从图5可知，CNN-Flow和传统光流法在预测时间较短时，精度都较高，但是当预测时间较长时，精度都严重下降，这是由于CNN-Flow和传统光流法都是仅仅分析了两帧图像而预测后续图像，随着时间增加，不确定性会逐渐增加，例如回波的突然消失和诞生，这种情况导致了精度的严重下降。CNN-Flow在CSI指标上和传统光流法在预测前期几乎相等，随着时间增加，逐渐略强于传统光流法；在FAR指标上，略差于传统光流法；在POD指标上始终略强于传统光流法。

图5 评价指标曲线图

为了在视觉角度分析本文算法在预测图像上的效果，图6是CNN-Flow和传统光流法在不同样本上预测的10帧图像中第1、4、7、10帧图像的对比。从图6可知，两个算法均能成功预测出回波图像的正确移动趋势，并且在预测时间较短时，预测的回波和真实的回波几乎一样，但是当预测时间增加时，两种算法对于图像细节的预测能力均有减弱。

图6 算法预测图像对比

4 结论

笔者在卷积神经网络的基础上结合传统的光流法提出了使用CNN-Flow网络及半拉格朗日算法进行雷达回波图像预测，该算法采用了SRAD数据集作为原数据进行实验分析，并在时间、降雨指标以及图像上对算法进行了效果的验证，在实验对比中发现该算法能够较准确地预测雷达回波图像，并在计算速度上远远超过了传统光流法。