基于ConvLSTM的改进雷达回波外推方法研究

赵玉娟，李宗飞，陈凯华，王 彦，朱男男，李祥海，姜罕盛

(1.天津市气象信息中心，天津 300074；2.天津市人工影响天气办公室，天津 300074；3.天津海洋中心气象台，天津 300074)

0 引言

临近预报通常指未来0～2 h的天气预报。强对流天气包括雷雨大风、冰雹等天气，因其历时短、破坏性强、演变规律复杂，是临近预报的重点和难点。海洋强对流天气临近预报在海洋气象灾害防御中具有重要地位。海上大风及其引发的次生灾害是导致海洋气象灾害的主要因素。为了提升海上大风等海上强对流天气预报能力，气象工作者持续在开展大风天气特征及预报方法等方面的探索。多普勒天气雷达是气象部门用于大气监测的重要设备，能够提供高时空分辨率的精细监测产品，在灾害性、突发性天气监测预警中是极为重要的参考指标，因此，许多学者开展了分析雷达回波在雷暴大风天气条件下特征的研究，为相关领域天气预报能力提升提供支撑。王彦等[1]利用天津地区46次雷暴大风过程统计分析了雷达回波在雷暴大风天气过程下的特征，得出了影响渤海西部的雷暴大风在雷达回波形态方面有弓状、带状、阵风锋等四种类型，弓状回波对应的雷暴大风天气最强烈等结论；郭庆利等[2]通过对烟台北部沿海5年的雷雨大风天气个例分析，得到了渤海海峡雷雨大风天气下雷达反射率因子的回波强度多在45 dbz以上，形状大致分带状、弓形等5类的结论。郭鸿鸣等[3]综合运用6部天气雷达拼图、WRF模式物理场等资料分析了强对流系统入海前后的时空变化规律。王亚南等[4]利用海岛、平台、浮标等站点加密观测资料，分析了渤海西部雷雨大风的统计特征。陈明轩等[5]、俞小鼎等[6]、程丛兰[7]等则论述了强对流天气临近预报方法。

雷达回波外推结果是临近预报主要参考数据之一。如何快速、准确地生成雷达回波预测数据是近年来气象领域研究热点之一。传统雷达外推方法包括交叉相关法(COTREC)[8]、光流法等[9-10]。交叉相关法和光流法均是假定雷达反射率因子的运动符合拉格朗日守恒规律，在稳定性降水预报中效果较好[11]，但对于局地强对流天气，因雷达回波演变快，不满足守恒条件，预报效果则会随着时间快速下降[12-13]。

探索雷达回波外推新方法，高效、准确地生成雷达回波预测数据，对提升海上强对流天气临近预报和服务能力具有重要意义。深度学习算法可自动学习海量数据中蕴含规律，且无需较多先验知识，因此在气象领域应用日益广泛，诸多学者开展了相关方法在短临天气预报方面研究。郭尚瓒等[14]开展了多层感知器在短时降雨预测方面的探索。Shi等[15]利用带卷积的长短时记忆单元构建RNN，将其在雷达外推预报方面应用取得了较好效果。施恩[16]提出了基于输入的动态卷积神经网络模型，卷积核含有当前输入的特征，网络模型测试期间还可基于输入图像变化，输入、输出图像的强相关性得到保证，利用南京、杭州、厦门三地雷达CAPPI图像数据作为样本试验的结果表明，所提方法较传统雷达回波外推方法，预测图像准确率和外推时效均有所提高。郭瀚阳[17]等借助基于自编码的卷积GRU网络，利用雷达拼图数据训练得到了可利用历史0.5 h数据预测未来1 h回波的雷达回波外推模型。试验结果证明其所提方法在预测精度上明显优于传统方法。黄兴友[18]等采用Causal-LSTM单元构建神经网络模型实现雷达回波外推，并使用带权重的损失函数进行模型训练，测试集及个例检验表明其构建的模型在强回波预报方面优于光流法。

目前基于深度学习的雷达回波外推研究多面向降水预报，针对海上大风天气的研究较少。本研究面向海上大风临近预报需求选取雷达数据，并从输入数据格式、损失函数两方面进行改进，构建了基于自编码的ConvLSTM网络，利用沧州雷达站4年的历史观测数据对其进行训练，得到了雷达回波外推模型，可利用近1 h雷达反射率因子数据预测未来1 h雷达反射率因子数据。测试集及典型天气个例对模型预测效果进行检验的结果表明，改进后模型较传统ConvLSTM模型在强回波预测方面效果更好。

1 方法

1.1 网络模型

本研究拟通过学习历史雷达数据时空变化规律，预测未来雷达回波序列，属于时间序列预测问题，循环神经网络在相关领域应用较多。不同于一般神经网络不同层神经元节点互相独立，循环神经网络各隐藏层节点不仅依赖当前输入，还依赖前一时刻中间状态，处理新数据时，也可记忆历史计算结果。卷积长短时记忆单元ConvLSTM是应用较广的循环神经网络模型之一，通过输入门、遗忘门及输出门实现信息流动控制，可防止有价值信息因为预测序列长度的增大而被丢，还能选择性地实现“更新”和“遗忘”。ConvLSTM工作原理如式(1)～(5)所示，it表示输入门，ft表示遗忘门，ot表示输出门，算子(o)表示矩阵对应元素相乘，“*”表示卷积操作，“σ”表示Sigmoid函数，Wx-，Wh-是二维卷积核，输入X1，…，Xt和单元状态C1，…，Ct，隐藏状态H1，…，Ht及it、ft、ot均为3维张量。

it=σ(Wxi*Xt+Whi*Ht-1+WciοCt-1+bi)

(1)

ft=σ(Wxf*Xt+Whf*Ht-1+WcfοCt-1+bf)

(2)

Ot=σ(Wxo*Xt+Who*Ht-1+WcoοCt+bo)

(3)

Ct=ftοCt-1+itοtanh(Wxc*Xt+Whc*Ht-1+bc)

(4)

Ht=Otοtanh(Ct)

(5)

循环神经网络还能作为基本单元构建更复杂网络。本研究借鉴了郭瀚阳等[17]研究思路，亦采用自编码模型构建网络，不同之处在于本文采用基于自编码的ConvLSTM网络进行雷达回波序列预测，通过多层堆叠ConvLSTM，增强模型学习能力。自编码模型包含编码、解码两阶段。编码阶段先用最后输入的隐藏层状态代表所有输入序列信息，再将编码最后一步隐藏层状态用于初始化解码阶段隐藏层的状态。解码阶段先用编码输入的最后一帧作为第一个输入得到第一个预测输出，再用第一个预测输出作为输入得到第二个预测输出，持续迭代此过程得到所有预测输出，该方法具备可产生变长预测序列的优势。本研究的编码和解码阶段均采用三层堆叠ConvLSTM来学习数据特征。

1.2 损失函数改进

不同天气状况下，雷达反射率因子强度具有明显差异。天气雷达实时观测到的回波强度是判断强对流天气的重要参考数据。预报经验和已有对海上强对流大风天气的雷达回波特征分析研究表明，雷达反射率因子包含较多强回波情况下，发生海上强对流大风天气的概率相对更高[19-22]，王福侠等[21]利用天气雷达和自动站资料研究统计2006～2008年河北省中南部地区28次雷暴大风过程中出现地面大风的262个观测站上空的雷达回波特征发现，一般雷暴大风的反射率因子都在50 dbz以上，但干对流雷暴大风的反射率因子一般只有40 dbz左右。强回波预测准确与否对预报效果影响更大，提高强回波预测准确率是提升预报效果的关键。

ConvLSTM模型常规采用的是均方差损失函数MSE，均方差损失函数先计算真实数据与预测数据所有对应点误差平方的总和，再计算其均值，对于不同大小的真实值权重一样，对于雷达回波数据则表现为对不同强度回波的权重一样。对于雷达回波预测应用而言，希望较大的回波值能够具有更高的预测准确率，为提升强回波预测能力，本文改进了损失函数，将损失函数构造为原始MSE损失函数和加权重的MSE函数两部分的组合，两部分的系数各为0.5，对于加权重的MSE损失函数，在利用均方差损失函数计算真实值与预测值差值基础上增加回波真实值作为权重，强回波的值较大权重会更高，弱回波值较小则权重会更低，从而实现提升强回波预测准确率的目标。改进后损失函数如式(6)所示，其中pred代表预测值集合，label代表实际观测值集合，N代表系数，改进后损失函数在label值较大即观测值为强回波时将产生更大影响，lossfunction为MSE函数，计算公式如式(7)所示，其中yi代表实际观测值集合label中的第i个真实观测值，yip代表第i个观测值在预测值集合pred中对应的预测值。基于公式(6)、(7)可推导得出公式(8)，由公式(8)可见，改进后的损失函数增加了N2yi2作为真实值与预测值差值的权重系数，当预测值与真实值存在差异且真实观测值较大时，改进后损失函数计算得出的数值较改进前更大，因此能够实现增大强回波权重，提升强回波预测效果的目的。

loss=0.5*lossfunction(pred，label)+

0.5*lossfunction(N*label*label，N*pred*label)

(6)

(7)

(8)

1.3 输入数据改进

气象部门长期存储的主要是二进制格式的雷达基数据，该类数据无法直接用于深度学习训练，已有研究多采用基数据生成图像文件方式作为输入数据，而单个图像文件仅能保存1个仰角的观测数据。天气雷达是对一定空间范围内降水回波的观测，一个体扫包含多个仰角观测数据。大气流体运动变化过程是在三维空间中进行，每个高度层均与其附近高度层有较强相关性[23]。在每个时刻输入相邻多个仰角观测数据有助于解决单层雷达回波图像外推局限性。数组是在程序设计中，为了处理方便，把具有相同数据类型的若干元素按有序的形式组织起来的一种数据组织方式，是用于储存多个相同类型数据的集合。每个观测时次，天气雷达单个仰角对一定空间范围的观测数据可抽象为二维数组，多个仰角的雷达观测数据可抽象为三维数组，npy文件是Python语言针对多维数组(Ndarray)的科学计算库NumPy专用的二进制文件格式，能够保存任意维度的NumPy数组，可满足存储多仰角雷达观测数据的需求，而且NumPy库提供了save、load函数为便捷地将数组数据存储到npy文件和从npy文件读取数组内容提供了有力支撑，因此，本文对深度学习模型的雷达输入数据存储格式进行改进，将同一个观测时次多个仰角的雷达反射率观测数据定义为(n，x，y)形式的三维数组进行存储，n代表观测仰角个数，x代表经度方向的观测数据点个数，y代表纬度方向的观测数据点个数，并利用numpy库的save函数将其存储在同一个npy文件中，N个观测时次的观测数据转化为N个npy文件存储，为后续实现多仰角输入数据训练打下基础，实际训练时可按需灵活提取npy文件中的1个或多仰角的雷达观测数据用于训练。

2 数据集构建

深度学习雷达数据集构建主要包括雷达数据筛选、数据预处理、深度学习样本组构造三个步骤，流程如图1所示。

图1 雷达数据深度学习数据集构建流程

2.1 雷达数据筛选

雷达数据时段选取强对流大风发生较多的5～9月。首先使用渤海西部有代表性的地面自动气象站、海上平台站逐小时观测资料，根据王亚男等[4]提出的指标选取强对流天气导致的雷暴大风过程，然后依据大风过程日期挑选雷达数据。强对流天气导致雷暴大风过程选取标准如下：

代表站瞬时风速(或最大风速)≥17 m/s，相应海域出现雷电天气(为消除系统性大风过程，对于瞬时风速(或最大风速)≥17 m/s且持续时间大于3小时的过程予以去除)。当多个观测站在不超过12小时范围内先后监测到雷雨大风并受同一天气系统影响时，记为一次过程。

2.2 数据预处理

雷达基数据预处理包括数据解码、坐标转换等步骤。CINRAD-SA型多普勒天气雷达基数据文件为二进制格式，存储的数据采用极坐标系，存储内容包括反射率、速度、谱宽等。数据解码步骤负责按照雷达基数据文件格式，完成雷达观测描述信息和反射率观测数据的提取并按不同信息的数据类型完成格式转换，坐标转换步骤首先建立极坐标和经纬度坐标的对应关系，然后将以极坐标形式存储的反射率等数据，投影到经纬度坐标系，得到反射率数据的格点矩阵。坐标转换步骤中建立极坐标和经纬度坐标对应关系的过程如下：雷达观测径向剖面图如图2所示，∠φ代表雷达观测仰角，d代表雷达站点与观测位置的径向距离。建立以地球中心为原点的三维直角坐标系，雷达站点r的坐标为(x1，y1，z1)，雷达观测任意点位置为p，其在地球表面投影为p’(x2，y2，z2)，Δlon，Δlat为投影点与雷达站点的经纬度差，观测点p在地球表面投影p’坐标利用式(9)计算得到。通过式(10)，计算得到任意p’位置(lon，lat)的极坐标值(∠A，∠φ，d)，其中∠A，d通过计算得到，∠φ为已知值，∠A代表雷达观测方位角。

图2 雷达观测径向剖面图

(9)

(10)

坐标转换后单个仰角的观测数据为400*400的格点矩阵，直接用于训练计算开销较大，为提高处理效率，本研究对数据进行了抽样处理，提取以雷达站点为中心±2度经纬度范围的格点数据，抽样处理后单个仰角观测数据的格点分辨率为100*100，数据量降低为抽样前的1/16。

在特定环境下，雷达观测受地物杂波、晴空回波等诸多因素影响，将产生非降水回波干扰。地物杂波已在雷达RPG中进行了处理，为进一步去除无效回波干扰，质量控制方面，本研究主要通过去除10 dbz以下晴空回波以保障模型训练效果。

2.3 深度学习样本组构造

雷达数据因观测设备故障、维护等原因存在缺测，导致数据存在时间间隔不等的问题，为去除时间不连续数据，按以下步骤构造观测数据文件时间间隔相同的深度学习数据集：

1)针对预处理得到的文件生成其文件名、观测时间索引信息。

2)利用索引信息，按照时间连续原则，筛选构造深度学习样本组。对符合时间连续性检查的数据文件使用长度为20的滑动窗口以步长为1进行滑动采样，得到所有样本组。每组样本包含20个观测时间连续的雷达数据文件{x1，x2，x3，…，x10，y1，y2，y3，…，y10}，x1到x10为输入数据，对应近一小时10个观测时间的雷达回波，y1到y10为输出数据，代表下一小时10个观测时间的雷达回波，每个文件时间间隔为6 min。

3 实验设计及结果分析

为验证本文方法有效性，利用Pytorch实现了深度学习模型，按照本研究方法，基于沧州雷达站2016～2019年5～9月观测数据，利用Python语言开发软件构建了雷达数据深度学习数据集，构建的训练集包含10 640个样本，测试集包含2 000个样本。激活函数是深度学习神经网络模型的重要组成部分，旨在使神经网络模型能更好地拟合数据分布，输出更准确的结果，其选择对神经网络性能、模型收敛速度有很大影响，本研究激活函数使用LeakyReLU.LeakyReLU是神经网络常用激活函数ReLU的变体。对于ReLU函数，当输入x>0时，输出为x，当输入x≤0时，输出始终为0，导致神经元不能更新。为了解决ReLU函数这一问题，在ReLU函数的负半区引入一个非常小的常数leak，即当x≤0时，输出为leak*x，使得输入信息小于0时，信息没有被完全丢掉，仍有很小的梯度。由于导数总是不为零，能减少静默神经元的出现，允许基于梯度的学习，解决了ReLU函数进入负区间后，导致神经元不学习的问题。训练过程采用反向传播算法计算误差，网络参数利用Adam算法更新，Adam算法本质上是带有动量项的RMSprop，它利用梯度的一阶矩估计和二阶矩估计动态调整每个参数的学习率，该算法的优点主要在于经过偏置校正后，每一次迭代学习率都有确定的范围，参数比较平稳，对内存需求少，适用于大数据集和高维空间。学习率参数的设置很重要，学习率取值过小会使网络收敛缓慢，过大则又会使训练易陷入局部最优解。本研究通过引入动量项加快网络收敛速度，减少震荡，使网络训练不易陷入局部最优解。损失函数采用原始MSE函数和带权重的MSE函数组合的方式，详见1.2节，实验过程对比了N设置为2、3、4情况下和输入数据为1个仰角、2个仰角时的模型预测效果，对比结果表明输入为1个仰角、N=2时预测效果最优，最终选定的损失函数参见式(11)。其他训练参数如表1所示，初始学习率采用0.000 1，模型的训练采用批训练方式，batch_size设置为20。训练采用早停策略(Early Stop)方式停止以避免过拟合。最大迭代次数设置为400次，当10次循环迭代后训练集精度提高不超过0.000 1，则训练结束。训练完成后，模型能够根据雷达近1小时10个观测时次的回波数据作为输入，预测出未来1小时10个时次的雷达回波数据。

表1 雷达数据深度学习训练参数

loss=0.5*lossfunction(pred，label)+

0.5*lossfunction(2*label*label，2*pred*label)

(11)

3.1 模型评价指标

采用分阈值评估方式对测试集进行检验，阈值分别选取15、20、30、40 dbz，预测时长包括0.5 h和1 h，预报时间间隔为6 min，评价指标采用气象领域常用指标临界成功指数(CSI)、命中率(POD)和虚警率(FAR)。通过逐点对比预测值与观测值，得出各预测点所属分类(分类标准参见表2)，进而计算得出评估指标。a、b、c分别代表预测数据命中数、空报数和漏报数，评估指标公式如式(12)～(14)：

表2 雷达回波像素点分类标准

(12)

(13)

(14)

3.2 测试集检验

表3对比了原模型和本文改进模型对于测试集在15、20、30、40 dBz共4个反射率阈值0.5 h和1 h的预测结果。由表3可知，本文改进模型在4个阈值的CSI、POD指标均明显优于改进前，FAR指标在阈值较小时略有增长，但在阈值较高时优于改进前。

表3 模型改进前后在测试集检验指标

POD指标方面，阈值为30 dbz及以下时，0.5 h、1 h的预测较改进前分别提高了9%～17%和13%～17%；阈值为40 dbz时，0.5 h、1 h的预测分别提高了8%和2%。CSI指标方面，阈值为30 dbz及以下时，0.5 h、1 h的预测较改进前分别提高了4%～12%和7%～10%；阈值为40 dbz时，0.5 h、1 h的预测较改进前分别提高了6%和1%。FAR指标方面，0.5 h的预测较改进前略有增长，1 h的预测在30、40 dbz较改进前降低。0.5 h的预测在阈值为30 dbz及以下时，较改进前增幅不超过6%，阈值为40 dbz的预测和改进前相当；1 h的预测在阈值为20 dbz及以下时，较改进前增幅不超过4%，在阈值为30、40 dbz时的预测较改进前分别降低1%和13%，优于改进前。

3.3 个例分析

文献[21]研究表明，产生雷暴大风的回波主要为弓形回波、带状回波和块状回波，其中带状回波是产生雷暴大风的主要回波。为更直观分析预测效果，从测试集和非测试集共选取2019年7月29日、2020年05月21日、2022年06月12日3组发生强对流大风的个例日期进行雷达回波外推分析，通过回波图像和预测指标两方面对比原始ConvLSTM模型和本文改进模型未来1 h雷达回波预测效果，对比图形包括四行，每行图像之间的时间间隔为6分钟，第一行为输入的近1小时雷达观测真实回波图像，第二行为未来1小时雷达实际观测的真实回波图像，第三、四行分别为原始ConvLSTM模型和改进模型基于第一行真实雷达观测预测的未来1小时雷达回波图像。

图3为测试集2019年07月29日09点06分到10点00分雷达回波个例预测效果对比，该个例是一次强带状回波引发的海上大风过程。带状回波通常是由多个对流回波单体相连排列成带状的回波，回波长度远大于回波宽度，有强回波时，传播方向与回波带垂直[22]。由图3可见，未来1 h的真实回波图像中在右下方一直存在一条形态较为明显、强度较高的带状回波，原始模型的预测，尤其在35 min之后的预测，对于带状回波的覆盖区域、强度方面均与实况有较大差异，而改进后模型预测的回波图像则整体上与观测实况更为相似，而且对带状强回波的预测效果有更为明显的提升，能较准确地预测出右下角带状回波形态、强度和移动位置，尤其是后30分钟的预测效果相比原始模型改进更为明显，较改进前更完整、清晰地预测了强回波形态、强度。该个例在阈值为30 dbz时，0.5 h、1 h的预测命中率分别为73%和62%，较改进前分别提升29%、32%；阈值为40 dbz时，0.5 h、1 h的预测命中率分别为37%和17%，较改进前分别提升36%、17%，虚警率皆低于改进前。不过也应注意到，改进后模型虽然在回波细节预测上优于改进前，但与观测实况仍有一定差距，而且预测细节与实际图像的差异随着预测时间的变长逐渐加大，外推时间越长，细节丢失也越来越多。

图4展示了非测试集个例2020年05月21日07点00分到07点54分雷达回波预测效果对比情况，该个例是一次小弓形回波引发的海上大风过程。弓形回波是指快速移动的向着运动方向凸起的，形如弓的强对流回波。弓形回波的空间尺度大小不一，小的弓形回波长度仅几十公里，有的可达上百公里.弓形回波是由后侧强烈的下沉气流造成的。显著弓形回波在低层反射率因子图上除了形如弓形外，弓形回波前沿存在高的反射率因子梯度，在较强回波带后侧有弱回波通道或者后侧入流缺口[21-22]。由图4中未来1 h的真实回波可见，图像上方持续存在一小弓形回波，并逐渐自左向右移动，改进后模型对回波整体形态预测更为准确，对于图像上部的回波强度、形态和移动趋势预测效果更优，能够更为明显地呈现出小弓形回波的形态特征和强度。该个例在阈值为30 dbz时，0.5 h、1 h的预测命中率分别为41%和29%，较改进前分别提升13%、11%；阈值为40 dbz时，0.5 h、1 h的预测命中率分别为27%和18%，较改进前分别提升15%、13%，虚警率皆低于改进前。

图4 模型改进前后2020年05月21日07点00分至07点54分预测效果对比

图5为非测试集个例2022年06月12日15点00分至15点54分雷达回波预测效果对比，该个例在未来1 h的真实雷达回波中存在两个小带状回波。由图5中可见，未来1 h的真实雷达回波中左侧带状回波形态较为稳定，且维持了较高的强度，右侧的带状回波则逐渐变弱、消散，原始模型预测左侧带状回波逐渐消散、右侧回波维持，而改进后模型则较准确地预测了未来1 h左下角回波一直保持较高强度、右下角回波逐渐减弱的变化趋势，与实况更为一致。该个例在阈值为30 dbz时，0.5 h、1 h的预测命中率分别为76%和73%，较改进前分别提升18%、36%；阈值为40 dbz时，0.5 h、1 h的预测命中率分别为30%和14%，较改进前皆提升6%，0.5 h预测的虚警率与改进前相当，1 h预测的虚警率较改进前分别降低18%、19%。

图5 模型改进前后2022年06月12日15点00分至15点54分预测效果对比

4 结束语

本文基于深度学习模型ConvLSTM，提出了从输入数据格式和损失函数两方面改进的雷达回波外推模型构建思路，通过沧州雷达观测站4年的历史数据对模型进行训练得到了可利用近1 h雷达观测数据预测未来1 h雷达回波的模型。测试集和非测试集典型个例检验结果表明，改进模型对于强回波的预测能力明显提升，预测的雷达强回波形态、强度、位置相比改进前与实况具有更高的相似性，能够为临近预报提供一定的参考。但改进模型在回波强度最大值、形态细节等方面的预测结果距离雷达真实观测仍有一定差距，同时，随着预测时间的增长，预测值与真实观测值的差异也逐渐加大，而且在外推预报的雷达数据产品方面，本研究仅探索了对反射率因子数据的预测，对于雷达观测数据中对大风预报也有较好指示意义的径向速度数据的预测尚未涉及。训练数据的质量和用于预测的模型均会对预测效果的提升产生影响，更丰富的预报产品也将为预报能力提升提供更好支撑，因此，下一步将继续开展基于更丰富和更高质量训练数据、增加注意力机制的深度学习模型以及关于雷达径向速度观测数据产品预测等方面的研究，进一步优化预测模型设计，丰富预测产品种类，提升预测效果，以便为大风临近预报提供更好支撑。