基于ConvLSTM的雷达回波外推预报

高云，郭艳萍，周建慧，张叶娥，杨泽民

山西大同大学计算机与网络工程学院，山西大同，037009

0 引言

随着全球气候变暖，工业污染加重，全球生态遭到严重破坏，引发的气候异常问题也越来越严重，中国的自然灾害是世界上最严重的国家之一。在各种重大自然灾害中，冰雹、暴雨、雷暴等突发性强的气象灾害，发生频率最高；又因其程度强而持续时间短而难以提前预知和防范，导致危害程度较大。如2005年07月，印度西部连日特大暴雨引发多起洪水和泥石流等灾害，造成至少900人死亡。2005年06月，黑龙江省宁安市洪灾死亡人数达106人，给人民的生命财产带来很大危害。

强对流天气是具有重大破坏性的灾害性天气之一。目前，通常使用雷达回波外推的方法对强对流天气进行监测和预警。雷达回波外推是指对回波运行的速度、方向、强度和形态变化等特征进行跟踪和预测[1-2]。传统雷达回波外推方法存在预测回波变化较快的天气过程时准确率大幅下降，预测时效较短，随着预测时长的增加准确性快速下降等不足。

早期的外推方法最常使用TREC技术，即交叉相关法。Hilst等先后在文献[3-6]中应用TREC技术对整体移动的风暴簇进行了跟踪，结果表明该技术能够对风暴簇的内部运动进行较好的反演。后经过多位学者的改进，徐亚钦等在文献[7]中提出了多重动态区域TREC算法，对TREC方法进行改进。单体质心法是较常用的一种回波外推法。乔春贵等在文献[8]中对整幅雷达回波图通过单体质心法进行线性外推，得出对于稳定的层状云该方法虽然有较好的外推能力，但对其他类型的云则外推结果较差。1950年Gibson在文献[9]中首先提出了光流法，张蕾等在文献[10]中对传统光流法进行了改进。陈家慧等在文献[11]中使用BP神经网络进行混合型降水回波的外推，表现出了较好的外推能力。

雷达回波外推方法近年来虽然取得了较大的进步，外推预报的准确性逐渐提高，但目前的雷达回波外推算法，对于1小时以上的外推预报，准确性下降很快。雷达回波外推方法还需要进一步地探索和改进。

1 ConvLSTM模型

ConvLSTM是LSTM结构的变体，通过将W的权值计算变成了卷积运算来提取出图像的特征。图1所示即为其单元结构图。

ConvLSTM的思路就是使用卷积来代替矩阵乘法。它在输入到状态以及状态到状态转换中都采用了卷积结构，通过堆叠多个ConvLSTM层并形成编码预测结构来建立更一般的时空序列预测模型。

该设计的一个显著特点是所有输入X1,...,Xt，细胞输出C1,...,Ct，隐藏状态H1,...,Ht，和ConvLSTM的3个门it，ft，ot都是3维张量，2维网格向3维张量的转换如图2所示。

ConvLSTM内部结构如图3所示。

其对应公式如下：

其中星号*表示卷积，小圆圈○表示哈达玛乘积。可以看出，这里在输入门、遗忘门和输出门这3个门的输入上，都考虑了细胞状态Ct-1。

为了确保细胞状态Ct-1具有与输入相同的行数和列数，在应用卷积运算之前需要先进行padding。将边界点上隐藏状态的填充视为使用外部世界的状态进行计算。通常，在第一个输入到来之前，将LSTM的所有状态初始化为零，这对应于对于未来的“完全无知”。

ConvLSTM模型用于预测网络时其编解码结构如图4所示。

ConvLSTM也可以作为更复杂结构的构建块。对于时空序列预测问题，使用图 4所示的结构。它包括编码网络和预测网络2个网络，将编码网络的最后状态复制到预测网络的初始状态和单元输出，2个网络都是通过堆叠多个ConvLSTM层形成的。预测目标与输入具有相同的维度，将预测网络中的所有状态连接起来并发送到1×1卷积层以生成最终预测。

2 实验过程

2.1 数据准备

本算法进行分析的原始数据序列为山西省大同市2019-6-1～2019-6-6共1191幅雷达站基本反射率图。雷达站名：大同；雷达扫描间隔时间为5′47"；数据范围200km；显示仰角：0.5/1.5/2.4。以原始基数据的命名方式进行命名，方便后期进行按顺序索引。其中2019-6-5 16:14:55BJT基本反射率图如图5所示。

（1）切割

从图 5可以看出，雷达站大同扫描范围为200km，除覆盖大同地区以外，还覆盖了山西省其他地区以及内蒙古自治区和河北省等部分地区。本实验只对大同地区的强对流气候进行预测，排除其他区域的干扰，先设定图片输出分辨率和初始化最终数组，处理的第一步是对原始雷达回波图像进行切割，只切割原始图像（200,150, 310, 350）的区域作为本次实验的数据。切割后的图像如图 6所示，只保留了大同市周边地区。同时减少比色卡、坐标、图片标题等干扰项，增加模型训练的精度。

（2）灰度转换

第二步是将三通道的彩色图片转换为单通道灰度图片。由于计算机硬件限制，原始绘制的300dpi高分辨率图像不适合输入到神经网络中直接进行训练，首先要将裁切后的图像降低分辨率（压缩）至100*100。原始雷达回波图像为RGB三通道图像，为减少训练时间可将原始图片转换为灰度图像，每个像素用8个bit表示，RGB转换为灰度图像的公式为：

将图像抗锯齿，取消字体平滑，抗混叠。每个雷达站每个时刻的体扫数据均可得到一个100*100分辨率大小的灰度图像，转换为灰度的图像如图7所示。

（3）得到训练集

增加一个维度备用，每幅图像载入内存后即为一个100*100*1大小的三维数组。处理完单个体扫数据之后再依次遍历所有图片，将每个图片得到的数组连接成一个四维的数组序列。为了排除杂点的干扰，再次将数组转化为2维（图像数量，100*100），依次遍历数组中每个点的值，如果该点的值＜50，则置为0。得到所需的图像数组序列。最后将所得数组保存为.npz格式的numpy数组方便后续取用，这就是最终可以输入到网络中进行训练的样本集。

2.2 ConvLSTM模型建立与分析

使用ConvLSTM模型进行雷达回波图像外推步骤如下：

（1）训练集和测试集划分

ConvLSTM模型样本数据的输入尺寸为如图8所示的 5D 张量（samples，time, rows，cols，channels），要提前将训练集和测试集reshape成如上形式的tensor张量，即样本总数、各样本帧数、图像宽度、图像高度和颜色通道数。输出尺寸为，如果返回全部序列，则返回5D张量，即（samples, timesteps, output_row, output_col, filters）；否则，返回4D张量，即（samples,output_row, output_col, filters）。o_row和o_col取决于filter和padding的尺寸。

假如上一层是ConvLSTM2D layer，那么其输出为以上形式的4D张量或5D张量，当后面再接另外一个layer时，就要考虑该layer是否能接受4D张量或5D张量，即要考虑ConvLSTM2D的输出能否作为该layer的输入。

数据准备生成的.npz格式的numpy数组是2维的，首先将其读入并reshape为4维数据，即（NUMBER, WIDTH, HEIGHT, 1），样本集总数为NUMBER，数据集为灰度图像，将颜色通道数设为1。将数据集设置为2个5D张量的数组，BASIC_SEQUENCE和NEXT_SEQUENCE，分别存放当前时刻图像和下一时刻图像作为训练集和结果集，数据为5维(NUMBER - FRAMES,FRAMES, WIDTH, HEIGHT, 1)，样本总量=原样本总数-各样本帧数。

训练时将整个样本集切割成长度统一的样本，16（FRAMES）帧图像为一个单元。共分为（NUMBER –FRAMES）个单元，即可直接训练。

（2）ConvLSTM模型建立

卷积长短时记忆神经网络自提出以来逐步完善已经成为了较为成熟的图像序列预测模型，目前Python中也有不少对其模型的封装，方便研究人员直接构建网络对自己的数据进行训练和预测。

本文所建模型为如图 9所示的Keras中的序贯模型，即逐层嵌套依次连接，数据流由输入端输入，逐层流动，反向传播时更新参数，逐步降低损失函数。该网络由四层ConvLSTM2D网络层堆叠而成，最后为一个三维卷积层用以格式化输出数组以便求取损失函数或获得预测结果。

本文所建模型由4层ConvLSTM2D网络层+1层ConvLSTM3D网络层堆叠而成，每层卷积核数目不同。4层ConvLSTM2D网络层卷积核数目分别为40、40、60、40；卷积核大小均为3*3；补0策略为′same′，即输出shape与输入shape相同；使用线性激活函数；返回输出序列中的全部序列。ConvLSTM3D网络层卷积核数目分别为1；卷积核大小为3*3*3；补0策略为′same′，即输出shape与输入shape相同；使用′sigmoid′激活函数；输出中维度的顺序为通道在后面。

为了保证在线性向非线性转变时，权重的尺度不变，所以使用BatchNormalization在激活函数前对输入进行了批标准化操作，优化神经网络，将分散的数据统一。

由于sigmod函数具有梯度饱和现象的缺点，所以本实验使用的是学习率为0.01的随机梯度下降优化器optimizers.SGD，将所有参数梯度裁剪到数值范围为[-0.5,0.5]的区间内。

2.3 ConvLSTM模型训练与预测

（1） 模型训练

在训练模型之前，需要通过compile来对学习过程进行配置。compile接收三个参数，分别为优化器、损失函数和指标列表。优化器optimizer为一个已预定义的优化器名，或者一个类型为Optimizer的对象；损失函数loss为目标函数，目的模型优化，可以使用预定义的损失函数名或一个自定义的损失函数；指标列表metrics的设置是针对分类问题，一般设metrics值为′accuracy′。指标可以是一个预定义指标的名字，也可以是一个用户定制的函数。指标函数应该返回单个张量，或一个完成metric_name - ＞ metric_value映射的字典。

本模型的损失函数为对数损失函数′binary_crossentropy′，以adadelta为优化方法进行编译。为了降低训练时间，在Keras中创建一个多GPU模型，设定4GPU并行处理。输入的数据即为2.1中处理的BASIC_SEQUENCE和NEXT_SEQUENCE作为训练集和结果集序列；batch_size值为10，即计算梯度所需的样本数量为10；epochs值为10，即代表样本集内所有的数据经过了10次训练；validation_split=0.05，即样本集中5%的数据为验证集。训练样本集，取其中一回合的部分训练结果如图10所示。

（2）预测

模型训练完成后存入“nice_model.h5”文件中，取样本集中第600组图片中的12帧进行预测，预测16次，训练结果连接形成预测集。

（3）可视化结果

将预测结果与训练所用样本集及其对应的结果集对比形成可视化结果。结果分为16组图片，前8组图片显示初始轨迹与标记数据的对比结果，如图11所示为第4幅图片。

后8组图片显示预测结果与标记数据对应结果的对比结果，如图 12所示为第12幅图片。

2.4 算法评价

从图10所示的模型训练结果来看，损失函数趋于减小的趋势，最终稳定在一个较小的值，说明模型的拟合程度较好。从图11的初始轨迹对比来看，轨迹运动与标记数据吻合度很好，从图12的预测结果来看，预测图像与标记数据的结果集吻合度基本一致，说明预测结果正确性好。

值得注意的是，模型图像的处理对后期结果的影响很大，由于不同时段雷达回波反射率因子分布区间不同，在进行灰度转换时，若使用原始数据绘制图像可能使比色卡（color map）的范围不同，这将导致不同图像中的相同强度的反射率因子值对应的RGB色彩不同，若不进行处理将极大影响训练效果。

另外，为了得到更为准确的外推结果，可在模型训练完成后，进行3次以上的外推，可以提高预测精度。

3 结语

本算法使用大同市历史气象雷达数据，重点介绍了ConvLSTM模型在雷达回波外推短期降水预测的应用方法，并详细介绍了ConvLSTM模型在雷达回波外推即时空序列分析建模的整个过程，同时对相应的算法和整个建模流程给出了实现及结果说明。本算法可以从以下方面进行拓展思考：

（1）本实验在建模时使用的损失函数是对数损失函数，损失函数有很多种，实际应用中均方误差函数也是常用的一种损失函数。后续实验中应考虑如何优化损失函数，选择使模型拟合性更好的损失函数。

（2）ConvLSTM是Xingjian Shi等提出的，后面他针对雷达图的天气预测又提出了TrajGRU，基于运行轨迹对图像做更精准的捕捉。本实验如何使用TrajGRU模型进行更好的预测，这些问题都需要进一步地进行探讨。