基于卷积神经网络的雷达回波外推方法

施 恩，李 骞，顾大权，赵章明

(国防科学技术大学 气象海洋学院，南京 211101)

0 引言

临近预报主要指0～3h的高时空分辨率的天气预报，主要预报对象包括强降水、大风、冰雹等灾害性天气[1]。目前，临近预报的主要手段是基于天气雷达资料的雷达回波外推技术，即根据当前时刻雷达观测结果，推测雷达回波未来的位置和强度，以实现对强对流系统的跟踪预报[2-3]。目前常用的雷达回波外推方法是质心跟踪法和交叉相关法(Tracking Radar Echoes by Correlation, TREC)[4]。质心跟踪法通过预测雷达回波的质心来推测下一时刻的回波位置，该方法对采集到的信息的利用更加充分，但它依赖阈值来识别风暴单体，仅适用于对风暴的追踪，难以应用于预测大范围降水的回波变化[5]; TREC通过对雷达回波划分区域，并求得各个区域当前时刻与前一时刻的相关系数，相关系数最大的区域即是回波移动矢量的终点[6]。研究人员在TREC算法的基础上，进一步发展了COTREC(Continuity Of TREC vectors)[7]和DITREC(Difference Image-based TREC)[8]等方法，但是此类方法仅根据几个时刻的回波特征推测下一时刻的回波分布，数据利用率较低，并且假设回波是线性演变的，而实际情况下回波的变化更为复杂。

针对上述传统雷达回波外推方法存在的问题，本文在传统卷积神经网络(Convolutional Neural Network, CNN)的基础上，提出了一种基于输入的动态卷积神经网络结构(Dynamic Convolutional Neural Networks based on Input，DCNN-I)，该网络从历史的时序雷达回波图像集中学习回波变化规律，从而实现对雷达回波图像的预测与外推。DCNN-I通过改进传统CNN的结构和网络训练过程中的参数，使网络中直接影响图像变换的卷积核在网络测试阶段仍然能够随着输入的不同而变化，更加适用于雷达回波外推这一类输入图像与输出图像之间存在较强相关性的问题，能够取得较好的雷达回波图像预测结果。

本文首先介绍DCNN-I的结构和特点，进而介绍利用该网络进行雷达回波外推的全过程。将大量有序的雷达回波强度等高平面显示 (Constant Altitude Plan Position Indicator, CAPPI) 图像集构造成训练样本集和测试样本集，利用训练样本集训练网络，使网络收敛，并用测试样本集测试网络，得到雷达回波外推图像。在对比实验中采用降水预报业务中常用的COTREC算法和DITREC算法作为对比方法，并从预测图像的准确率和外推时效两方面进行对比。

1 基于输入的动态卷积神经网络

卷积神经网络是深度学习的重要分支，由于CNN具有特殊的网络结构，在图像处理领域具有深远的研究前景，近年来得到高度重视并引起广泛研究[9-12]。CNN是一种具有多层结构的网络，其低层由卷积层和下采样层交替组成，是网络提取图像特征的重要环节，高层为分类器，一般为全连接层[12]。

考虑到雷达回波外推问题的输入与输出之间存在较强相关性，本文在传统卷积神经网络基础上提出一种基于输入的动态卷积神经网络结构(DCNN-I)，增加了动态子网络(Dynamic Sub-Network, DSN)和概率预测层(Probability Prediction Layer, PPL)，PPL中的卷积核由DSN计算得出，训练阶段结束之后卷积核可根据输入图像序列的不同而动态变化。DCNN-I的整体结构如图1所示。

图1 DCNN-I的结构 Fig. 1 Architecture of DCNN-I

其中DSN与传统的CNN具有相同的结构，用于计算概率向量，PPL中将输入图像序列中的最后一幅雷达回波图像与卷积核进行两次卷积操作得到最终的外推图像。DSN中包含DCNN-I的全部权值参数，其输出作为PPL的卷积核。

在传统的卷积神经网络中，卷积核在网络训练阶段通过反向传播算法不断得到更新，而在测试阶段保持不变。与传统的CNN不同，DCNN-I中PPL的卷积核是DSN对输入图像序列进行处理的结果，因此在网络测试阶段仍会随着输入样本的不同而变化，使网络具有动态特性，对于雷达回波外推这一类输出图像与输入图像强相关的问题具有更好的适用性。

1.1 动态子网络DSN

在DCNN-I中，输入的雷达回波图像首先输入动态子网络DSN，经过卷积、采样以及激活函数处理后，得到一维列向量VPV和一维行向量HPV。DSN可以视为一个独立的卷积神经网络，其低层包含5个卷积层和4个下采样层，高层为1个全连接层，DSN的结构如图2所示。

DSN中包含C1、C2、C3、C4和C5共5个卷积层，图2中给出了各个卷积层包含的卷积核个数(在数值上与输出特征图数量相同)以及输出特征图的分辨率，其中C5层的卷积核大小为7×7，其余各个卷积层的卷积核大小均为9×9。在卷积层中，输入特征图与卷积核相卷积，卷积结果加上一个偏置项后作为激活函数的输入，经过激活函数处理后得到该层的输出特征图。DSN中的卷积层的计算公式为：

(1)

图2 DSN的结构 Fig. 2 Architecture of DSN

DSN中包含S1、S2、S3和S4共4个下采样层，图2中给出了各个下采样层输出特征图的数量和分辨率，各层中采样核的大小为2×2，采样步长为2。在下采样层中，对输入特征图降分辨率，在提取图像特征的同时降低计算复杂度，采样时以步长为2进行连续采样，采样过后输出特征图的分辨率降为原输入特征图的1/4。DSN中的下采样层的计算公式为:

(2)

在DSN的全连接层F1中，首先将C5层的32个分辨率为16×16的输出特征图按列顺序展开，得到大小为521×1的列向量。分别计算该向量与垂直参数矩阵WV和水平参数矩阵WH的外积，计算结果分别加上一个偏置向量，再经过Softmax函数处理之后，得到一维列向量VPV和一维行向量HPV。DSN中的F1层的计算公式为：

(3)

其中：aC5表示C5层输出的大小为512×1的特征向量，WV和WH均为41×512的参数矩阵，BV和BH均为41×1的偏置向量，计算VPV时需要将计算结果转置，最终得到大小为1×41的VPV和41×1的HPV。由于Softmax函数的特性，得到的一维向量VPV和HPV中所有元素均为正值且元素之和为1，因此可以将VPV视为垂直概率向量，将HPV视为水平概率向量。

1.2 概率预测层PPL

PPL包含两个输入，输入图像序列中的最后一幅图像作为PPL的输入特征图，DSN的输出概率向量VPV和HPV作为PPL的卷积核，在PPL中，将对其输入特征图进行两次卷积操作，得到外推图像，PPL的结构如图3所示。

PPL包含DC1和DC2两层结构，垂直概率向量VPV为DC1层的卷积核，水平概率向量HPV为DC2层的卷积核。在PPL中进行两次卷积操作，在DC1层中，将输入图像序列中的最后一幅雷达回波图像与VPV相卷积，得到分辨率为240×280的输出特征图；在DC2层中，将DC1层的输出特征图与HPV相卷积，得到分辨率为240×240的预测图像。由于概率向量的特殊性(元素均为正，元素和为1)，可以将DC1层和DC2层中的卷积操作视为垂直方向的预测和水平方向的预测。

图3 PPL的结构 Fig. 3 Architecture of PPL

2 数据集构造

实验中采用南京、杭州、厦门三市的2016年7至9月份CINRAD-SA型多普勒天气雷达资料训练DCNN-I，数据源中大约包含60 000个样本，可满足神经网络训练需求。由于雷达系统中为用户提供可视化界面的PUP软件没有批量生成雷达回波图像的功能，因此本文需要研究雷达回波强度CAPPI图像生成方法，并构造数据集。

CINRAD-SA型多普勒天气雷达基数据为三维极坐标下的散乱数据，需要依次经过坐标转换、数据插值、水平采样和绘制才能够得到一幅回波强度CAPPI图像， 图4为多普勒天气雷达基数据预处理过程。

图4 多普勒天气雷达基数据预处理过程 Fig. 4 Preprocessing process of Doppler weather radar data

首先根据文献[13]中提供的方法，将三维极坐标下的数据转换到三维笛卡尔直角坐标系中，采用反距离加权法[14]进行数据插值，得到三维笛卡尔直角坐标系下的规整网格数据。然后对数据进行水平采样，提取某一高度下的二维平面数据，将数据映射至0～255，便能够得到回波强度CAPPI灰度图像，每幅图像的分辨率为2 000×2 000。为了降低CNN的训练复杂度，不对图像进行颜色映射。

雷达基数据经过数据预处理之后得到回波强度CAPPI灰度图像，此时图像的边缘部分包含大量的空白信息，需要进行裁剪和压缩处理，得到分辨率为280×280的灰度图。然后再根据图像之间的时间顺序进行采样，构造网络训练和测试需要的样本，得到的每一组样本包含5幅图像{x1,x2,x3,x4,y}，每幅图像之间的时间间隔恒定(6 min)。其中{x1,x2,x3,x4}是分辨率为280×280的输入样本图像序列，{y}是再次经过裁剪、分辨率为240×240的对照标签(网络的期望输出)，此时对第5幅图像进行裁剪的目的是保证对照标签与DCNN-I的输出图像的分辨率相一致。最终得到包含48 000幅样本图像的训练样本集和包含4 800幅样本图像的测试样本集，图5为数据集中的一组样本图像。

图5 数据集中的一组样本图像 Fig. 5 A group of sample images in the dataset

3 实验结果与分析

为了验证本文方法的有效性，在2.40 GHz CPU、内存大小为4 GB的PC上实现了算法。其中卷积核和参数矩阵(VPV和HPV)的初始化依据Xavier初始化方法[15]，将偏置初始化为0向量。训练过程采用反向传播算法计算误差，并采取带动量的梯度下降法更新网络参数，学习率为0.000 1，动量系数为0.5。网络训练阶段采取批训练的方式，每次输入10组样本，网络迭代次数为40。

对比实验中，对2016年7月7日南京地区、2016年7月16日杭州地区以及2016年8月11日厦门地区三个大范围降水过程进行回波强度的外推，并分别从回波图像和预测指标两方面分析不同雷达回波外推方法的优劣。

3.1 回波图像对比分析

对比实验中分别对3个降水过程进行雷达回波强度外推，结果如图6所示，图中的第1行为输入的雷达回波图像，第2行为实际观测到的雷达回波图像，第3行和第4行分别为采用DCNN-I和COTREC的外推结果，每一行中相邻的图像之间的时间间隔为6 min。

由图6可知，与COTREC算法相比，基于DCNN-I的方法外推图像与观测到的实况图像相似程度更高；尤其是在图像边缘部分，采用COTREC算法很难准确计算出观测区域边缘部分的矢量场，而DCNN-I可以在训练过程中学习观测区域边缘部分的回波变化特征，从而给出较为可靠的预测。此外，图6中采用DCNN-I方法得到的外推图像都较为模糊，这是因为雷达回波外推本身存在的不确定性，其输入和输出并不是简单的对应关系，DCNN-I在训练过程中学习了各种天气过程回波变化规律，网络在预测时几乎不可能对整个观测区域内雷达回波作出精准的预测，导致外推的回波图像较为模糊。

3.2 预测指标对比分析

在3.1节中，本文从图像层面分析外推的雷达回波图与观测到的雷达回波实况图的相似性，对比分析了外推结果。除了根据图像的相似性来判断外推准确性以外，还可以通过预测指标来评价外推结果。本文中采用的预测指标包括: 临界成功指数(Critical Success Index，CSI)、误报率(False Alarm Rate，FAR)以及探测概率(Probability Of Detection，POD)[16]。

图6 三个降水过程的预测结果 Fig. 6 Prediction results of the three example

首先设定一个回波强度阈值(对于回波图像，即代表某一灰度阈值)，若雷达回波图像中某个像素的灰度值超过该阈值，则判定该像素对应区域的预报值是活跃的，反之，则判定该像素的预报值为不活跃。因此对于单个像素的预测分为成功(预报值与实际值均活跃，标记为S)、漏报(预报值不活跃，实际值活跃，标记为M)和空报(预报值活跃，实际值不活跃，标记为F)。本文以nS、nM和nF分别表示预测图像中成功、漏报和空报的格点数，则定义CSI、FAR和POD的计算公式如下：

表1 不同雷达回波外推方法的结果对比Tab. 1 Comparison of mean and variance for different method over 15 prediction steps

(4)

CSI、FAR和POD能够较好地评估预测的雷达回波图像的预报效果，但是仍然不能直观地反映预测的雷达回波的准确率，因此本文进一步采用均方误差(Mean Square Error, MSE)来衡量雷达回波外推的准确率，其计算形式为：

(5)

进行计算之前，需要先根据Z-R关系[17]，将图像中的像素每个像素的值换算成每个格点的降水强度。上式中，N为观测区域内Ω的格点总数，实际上就是指图像包含的像素点的总数，F为实际观测的雷达回波图像，F′为预测的雷达回波图像，t0为预报时刻，τ为预测的时效，F′(t0+τ,x)表示在t0+τ时刻，预测的雷达回波在x格点处的降水强度。

实验中，采用COTREC算法和DITREC算法作为对比方法，此外还将输入图像中的最后一幅( Last Image)作为对比的依据。对2016年8月11日厦门地区大范围降水过程进行连续15次外推预测，并分别计算4种方法预测结果的MSE、CSI、FAR和POD预测指标，对比结果如图7所示。

实验中回波强度阈值设置为10 dBZ，对应图像的灰度值为25。由图7可知，随着预测时间的增长，实况回波图像与Last Image之间的差别越来越大，COTREC算法与DITREC算法的预测结果较为接近，基于DCNN-I的雷达回波外推方法预测结果最好，尤其是在准确率(以MSE表示)和误报率方面，采用DCNN-I优化效果明显。该降水过程的15次外推过程中，各项预测指标的平均值和方差如表1所示。

从表1可得，采用DCNN-I的雷达回波外推方法能够得到更好的预测结果，预测结果也较为稳定，预测准确率约有10%的提升。

本文通过去相关时间L来描述不同雷达回波外推方法的最大预测时效[18]。去相关时间L与预测的雷达回波图像与观测到的实况回波图像之间的相关系数c的值有关，相关系数c的计算公式为：

(6)

其中，G与G′分别表示实况雷达回波图像和预测的雷达回波图像中某个像素点的灰度值。相关系数的取值为c∈[0,1]，c的值越大，表明两幅图像越相似，雷达回波外推的效果越好。由于相关系数c符合指数规律，因此将相关系数下降到1/e时对应的时间定义为去相关时间L。当相关系数c低于1/e时，可以判定预测的雷达回波不再具有有效性，因此去相关时间L能够作为衡量雷达回波外推时效的依据。图8为分别采取COTREC、DITREC和DCNN-I方法对三个降水过程进行外推时，相关系数关于预测时间的变化。

图7 不同方法外推结果的预测指标对比 Fig. 7 Comparison of four precipitation nowcasting metrics on different method

图8 三次降水过程中相关系数的变化 Fig. 8 Correlation coefficient changes over forecast time during three precipitation processes

由图8可知：在三次降水过程中，相关系数随着预测时间呈指数递减；与COTREC算法和DITREC算法相比，采用DCNN-I的外推方法对应的相关系数随时间下降更为缓慢。从去相关时间L来看，采用传统方法进行雷达回波外推得到的L值约为2.4～3 h，采用DCNN-I的雷达外推方法得到的L则大于4 h，证明基于DCNN-I的雷达回波外推方法能够有效延长外推时效。

4 结语

卷积神经网络在二维图像的识别和预测上具有深远的应用前景，目前卷积神经网络在天气预报领域的应用较少，本文尝试利用卷积神经网络进行雷达回波外推，在传统卷积神经网络的基础之上，通过增加DSN和PPL，提出了基于输入的动态卷积神经网络DCNN-I，该网络中PPL的卷积核由DSN计算得出，卷积核中包含当前输入图像的特征，在网络测试阶段仍然能够根据输入图像的不同而变化，从而保证网络输入图像和输出图像之间的强相关性。与传统的雷达回波外推方法相比，基于DCNN-I的雷达回波外推方法通过网络训练从大量的雷达数据中学习雷达回波变化特征，有效提高了数据的利用率，对于“陌生”的雷达回波具有较强的适应性。实验中利用南京、杭州、厦门三个地区的雷达回波图训练DCNN-I，并将外推结果与降水预报业务中常用的COTREC算法和DITREC算法相比，结果表明，基于DCNN-I的雷达回波外推方法能够提高外推回波图像的准确率，并有效提高外推时效。

在下一步的研究工作中将进一步优化DCNN-I，降低预测图像的 “模糊” 程度，并将其应用于其他的图像识别与预测领域。

参考文献(References)

[1] DOVIAK R J, ZRNIC D S. Doppler Radar & Weather Observations [M]. Waltham, MA: Academic Press, 2014: 1-9.

[2] 张沛源, 杨洪平, 胡绍萍. 新一代天气雷达在临近预报和灾害性天气警报中的应用[J].气象,2008,34(1):3-11.(ZHANG P Y, YANG H P, HU S P. Applications of new generation weather radar to nowcasting and warning of severe weather [J]. Meteorological Monthly, 2008, 34(1): 3-11.)

[3] OTSUKA S, TUERHONG G, KIKUCHI R, et al. Precipitation nowcasting with three-dimensional space-time extrapolation of dense and frequent phased-array weather radar observations [J]. Weather and Forecasting, 2016, 31(1): 329-340.

[4] RINEHART R E, GARVEY E T. Three-dimensional storm motion detection by conventional weather radar [J]. Nature, 1978, 273(5660): 287-289.

[5] 李英俊,韩雷.基于三维雷达图像数据的风暴体追踪算法研究[J].计算机应用,2008,28(4):1078-1080.(LI Y J, HAN L. Storm tracking algorithm developmentbased on the three-dimensional radar image data [J]. Journal of Computer Applications, 2008, 28(4): 1078-1080.)

[6] LIANG Q Q, FENG Y R, DENG W J, et al. A composite approach of radar echo extrapolation based on TREC vectors in combination with model-predicted winds [J]. Advances in Atmospheric Sciences, 2010, 27(5): 1119-1130.

[7] FLETCHER T D, ANDRIEU H, HAMEL P. Understanding, management and modelling of urban hydrology and its consequences for receiving waters: a state of the art [J]. Advances in Water Resources, 2013, 51(1): 261-279.

[8] 张亚萍,程明虎,夏文梅,等.天气雷达回波运动场估测及在降水临近预报中的应用[J].气象学报,2006,64(5):631-646.(ZHANG Y P, CHENG M H, XIA W M, et al. Estimation of weather radar echo motion field and its application to precipitation nowcasting [J]. Acta Meteor Sinica, 2006, 64(5): 631-646.)

[9] KRIZHEVSKY A, SUTSKEVER I, HINTON G E. Imagenet classification with deep convolutional neural networks [C]// Proceedings of the 25th International Conference on Neural Information Processing Systems. [S.l.]: Curran Associates Inc., 2012,1: 1097-1105.

[10] 常亮,邓小明,周明全,等.图像理解中的卷积神经网络[J].自动化学报,2016,42(9):1300-1312.(CHANG L, DENG X M, ZHOU M Q, et al. Convolutional neural networks in image understanding [J]. Acta Automatic Sinica, 2016, 42(9): 1300-1312.)

[11] 李倩玉,蒋建国,齐美彬.基于改进深层网络的人脸识别算法[J].电子学报,2017,45(3):619-625.(LI Q Y, JIANG J G, QI M B. Face recognition algorithm based on improved deep networks [J]. Acta Electronica Sinica, 2017, 45(3): 619-625.)

[12] LeCUN Y, HUANG F J. Large-scale learning with SVM and convolutional for generic object categorization [C]// CVPR ’06: Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Washington, DC: IEEE Computer Society, 2006, 1: 284-291.

[13] 张志强,刘黎平,王红艳.三维可视化技术在雷达三维组网产品显示中的运用[J].气象科技,2010,38(5):605-608.(ZHANG Z Q, LIU L P, WANG H Y. Application of 3D visualization technology to display of Doppler radar networking products [J]. Metrological Science and Technology, 2010, 38(5): 605-608.)

[14] LU G Y, WONG D W. An adaptive inverse-distance weighting spatial interpolation technique [J]. Computers & Geosciences, 2008, 34(9): 1044-1055.

[15] GLOROT X, BENGIO Y. Understanding the difficulty of training deep feedforward neural networks [EB/OL]. [2017- 03- 01]. http://www.weblio.jp/redirect?url=http%3A%2F%2Fjmlr.org%2Fproceedings%2Fpapers%2Fv9%2Fglorot10a%2Fglorot10a.pdf&etd=3e03a9174d723be0.

[16] 王改利, 赵翠光, 刘黎平,等. 雷达回波外推预报的误差分析[J]. 高原气象, 2013, 32(3):874-883. (WANG G L, ZHAO C G, LIU L P, et al. Error analysis of radar echo extrapolation [J]. Plateau Meteorology, 2013, 32(3):874-883.)

[17] UIJLENHOET R, POMEROY J H. Raindrop size distribution and radar reflectivity-rain rate relationships for radar hydrology [J]. Hydrology & Earth System Sciences & Discussions, 2001, 5(4): 3012-3018.

[18] SHI X, CHEN Z, WANG H, et al. Convolutional LSTM network: a machine learning approach for precipitation nowcasting [C]// Proceedings of the 28th International Conference on Neural Information Processing Systems. Cambridge, MA: MIT Press, 2015,1: 802-810.

This work is partially supported by the National Natural Science Foundation of China (41305138, 61473310).

SHIEn, born in 1992, M. S. candidate. His research interests include neural network, pattern recognition.

LIQian, born in 1980, Ph. D., associate professor. His research interests include pattern recognition.

GUDaquan, born in 1959, professor. His research interests include intelligent information system, neural network.

ZHAOZhangming, born in 1993, M. S. candidate. His research interests include pattern recognition, distributed network.