深圳市气象局临近预报技术进展

陈元昭 兰红平 刘琨

（深圳市气象局，深圳 518030）

0 引言

基于雷达数据的雷暴识别追踪和外推预报方法是最早出现的临近预报技术[1]。20世纪50年代，Ligda[2]首先开始了雷达回波外推进行临近预报研究，此后临近预报技术得到不断发展，主要体现在一些预报方法的提出。单体质心法的核心是首先将整个雷暴单体视为一个三维的单体，然后对单体进行识别、分析后，再进行追踪，采用拟合外推法来做雷暴的临近预报[3-5]。交叉相关法是通过计算整块回波不同区域的连续时次的相关系数，采用最优相关法，来获得回波的运动矢量场，然后通过这些运动矢量场来预报回波未来的位置和形状。交叉相关法的优势在于它的计算方法相对简单，不但可以跟踪对流降水，还可以对层状云降水进行追踪[6-8]，在气象业务部门一度得到广泛应用[9-10]。但交叉相关法对强度和形状随时间变化很快、局地生成的回波，常会出现跟踪失败的情况[11-12]。近年来，基于雷达的临近预报方法发展迅速，中国气象局在完成灾害性天气短时临近预报系统（SWAN）[13]的基础上进一步完善了强对流天气的监测技术，研发了基于闪电和卫星产品的临近预报产品[14]。光流法是计算机视觉领域中的重要方法[15]，近年来该方法在天气临近预报领域得到广泛应用[16]，但对于西风带系统，尤其是飑线的预报，光流法存在一定的局限性[17]；而未来几年以雷达外推技术为主的临近预报将继续在业务中应用[18-19]。近年来，随着计算能力、信息技术和智能算法技术的不断突破，人工智能技术备受社会各界的高度重视，发展迅猛[20-21]，基于深度学习的智能天气临近预报方法也受到气象领域的高度关注。

本文主要从技术角度对深圳市气象局的临近预报技术的发展历程进行回顾总结。第一部分主要介绍交叉相关法、光流法、粒子滤波临近预报方法以及雷达基数据质量控制方法，第二部分主要介绍基于生成对抗网络（generative adversarial networks，GAN）[22]的人工智能临近预报方法，最后对不同方法的优缺点进行总结和讨论。

1 基于运动矢量场的临近预报方法研究进展

深圳市气象局临近预报方法的研究始于2006年，在引进和吸收了国内外对流风暴追踪技术算法的基础上，2008年研发了交叉相关临近预报方法，并应用到2011年深圳世界大学生运动会气象服务中。随后逐渐研发了局部约束光流法、粒子滤波融合法等临近预报技术。以上述方法为核心搭建了临近预报决策支持平台（PONDS）系统，并在业务中实时运行，对广东地区的强对流天气过程进行实时监测和临近预报，成为深圳市气象局三大核心业务系统之一。

1.1 雷达资料处理

深圳临近预报采用广东省12部S波段天气雷达（分别位于广州、深圳、韶关、珠海、清远、阳江、河源、汕尾、汕头、梅州、湛江、肇庆）2.5 km高度的反射率因子雷达拼图资料。雷达的虚假回波经常出现在低层，选取2.5 km高度的回波拼图作为预报的初始场。这一层的雷达回波既能代表华南地区对流的水平分布特征，同时在很大程度上避免了虚假回波和地物杂波的影响。为计算方便，雷达资料从极坐标格式利用最近邻居法和垂直方向的线性内插法相结合的插值法插值到三维直角坐标系中。拼图的时间间隔为6 min。拼图前对雷达基数据资料进行质量控制，剔除了超折射、地物等的影响。

1.2 交叉相关临近预报方法

交叉相关法的核心就是通过计算连续时次雷达回波不同区域的最优空间相关系数，来确定回波在过去的移动矢量特征。交叉相关法把整个数据区域划分成若干小区域，然后在相邻时刻雷达回波图像的小区域之间计算相关系数，通过最大相关系数确定相邻时刻图像中的区域对应关系，进而确定回波区域的平均运动。一个简单的示意图如图1所示。将相邻时刻（t1，t2）雷达反射率因子回波图像（图1），分别划分成若干大小相同的子区域。对于t1时刻的一个子区域A，对应t2时刻回波图像的候选区域（大小由平均风速决定），将所有可能的子区域分别与A做相关计算，相关系数R表示为：

图1 交叉相关法示意图 Fig. 1 Cross-correlation approach

式中，1和2分别是1和2时刻的反射率因子，是一个子区域内象素点的数量。t2时刻相关系数最大的子区域B即为与A匹配的子区域。将A和B的中心位置连接起来，即为回波运动矢量。

交叉相关法仅利用最近两个时次的雷达回波获得的运动矢量场是不稳定的，用这种风场外推的预报结果也具有较大误差。雷达标定、地物及异常传播的电磁波受到阻碍物遮挡、衰减、雷达波束的不完全充塞等随机因素是出现这种误差的主要原因。在进行运动矢量场计算时，需要抑制风场中的这些噪声干扰，以提高回波预报稳定性和准确性。在交叉相关法中利用卡尔曼滤波法改进雷达回波交叉相关算法[23]，具体算法见[24-25]，计算中考虑了过去0.5～1h的回波趋势。从效果上看，提升了回波预报的稳定性[23]。

在变化较为平缓的层状云降水系统中，从雷达回波图像上观察到的雷暴的外形和移动速度变化都不是很剧烈，此时交叉相关法的效果较好[11]。但对强对流天气过程，尤其是在雷暴的外形和移动速度随时间变化很快的情况下，通过简单的计算相关系数的方法难以保证追踪的准确性，交叉相关法的预报效果就会明显降低，跟踪失败的情况显著增加，并最终影响预报的结果[11]。所以，对变化较为剧烈的强对流降水系统，交叉相关法有先天的缺陷，需要引入新的方法。

1.3 局部约束光流法临近预报方法

1.3.1 Lucas-Kanade 法计算光流场

光流法的基本原理是以图像亮度变化来识别运动目标和观测器之间的相对运动产生的瞬时位移。图像中所有像素点的瞬时位移就构成了图像的光流场。而光流法的核心就是通过识别出的连续图像系列计算光流场。简单的来说，光流场就是类刚体物体的速度矢量场。光流方程如下：

计算光流场的约束方法通常有Lucas-Kanade局部约束[26]和Horn-Schunk全局约束[27]。局部约束方法是对某个点周围给定小区域的光流给定限定条件，而全局约束方法是在整个图像区域范围内满足一定的约束条件。临近预报关注回波的局地变化，光流法中采用Lucas-Kanade局部约束法作为计算光流的约束条件。用获得的光流场进行临近预报。

与传统的交叉相关法相比，光流法的优势在于立足于变化，而不是选定不变特征再跟踪不变特征移动的方式。

1.3.2 质量控制——中值滤波

在光流法中，对雷达数据的质量控制采用中值滤波法。中值滤波是一种非线性平滑技术，它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值[28]。

2013年4月30日，受锋面低槽影响，一条飑线自西北向东南橫扫广东全省。本次过程回波滤波前后对比可以看出（图2），滤波前后飑线的长度、形状、强回波中心等特征大体一致。但可以看出，滤波前（图2a），飑线后部南端的层状云回波比较零碎，杂乱。滤波后（图2b）去噪效果好，图像清晰，较好地保留了回波的边缘特征，没有出现明显失真。

图2 2013年4月30日14时（北京时，下同）广东雷达拼图光流法追踪得到的回波：（a）中值滤波前；（b）中值滤波后 Fig. 2 Optical flow method based on the Guangdong multi-radar mosaic at 14:00 BT 30 April 2013: (a) before median filtering; (b) after median filtering

1.3.3 预报效果个例对比

2014年5月20日中午，一条近东西向的弓形回波从珠江口西侧向东移动，追上原位于其东侧的正在减弱的回波。这两条回波合并后，迅速增加，在珠江口附近发展成为倒“Y”字形的强回波。给珠江口一带带来了局地强降水和短时雷雨大风。

光流法对这次强降水过程30 min外推预报效果较好。从图3可知，光流法预报出了回波合并后强度加强，预报的30 min后回波的位置与实况基本相符，也基本预报出了回波的倒“Y”字形结构，只是西南部的回波预报比实况偏弱。

对华南地区回波预报评估结果表明[16]，光流法整体预报效果优于交叉相关法，尤其对移动型局地生成的回波及强度和形状随时间变化很快的雷暴。

光流法对西风带系统，尤其是飑线的预报仍存在一定的局限性[19]，需要研制新方法来弥补雷达临近预报技术在0～1 h的预报能力方面的缺陷。

1.4 粒子滤波临近预报方法

1.4.1 粒子滤波法基本原理

粒子滤波（Particle filter）由Carpenter等[29]首次提出，主要是通过非参数化的蒙特卡洛（Monte Carlo）模拟来实现递推贝叶斯（Bayes）滤波估计。

在进行雷达外推预报前需获得雷达回波的运动矢量场。和获得唯一运动矢量场不同，粒子滤波算法是利用同一时次的雷达回波采用两种不同的估测方法得到的相同时次的两组回波运动矢量场，然后采用粒子滤波算法对两组运动矢量场进行融合，再利用得到的运动矢量场来进行回波的外推预报。在粒子滤波算法中，回波矢量场可表示为：

图3 光流法对2014年5月20日珠江口附近迅速加强回波预报：（a）14时实况；（b）14时的30 min外推预报；（c）14：30雷达回波实况 Fig. 3 Optical flow method based on the Guangdong multi-radar mosaic on 20 May 2014:(a) real radar echo at 14:00 BT; (b) 30 min extrapolation forecast at 14:00 BT; (c) real radar echo at 14:30 BT

式中，f(x)、h(x)分别为状态方程和观测方程。x(t)、y(t)、m(t)、n(t)分别表示系统的状态、观测值、过程噪声和观测噪声。其中f(x)为采用Lucas-Kanade约束的光流法得到的运动矢量场，而h(x)采用基于Harris角点算法[30]追踪回波得到运动矢量场。

图4 粒子滤波法流程图 Fig. 4 Flow chat of particle filtering fusion algorithm

图4为粒子滤波法流程图。粒子滤波计算方法的核心是用一组加权随机样本（粒子）来近似表示后验概率密度函数。递推贝叶斯滤波理论给出了该问题的严格求解。在进行求解时引入了优化算法对粒子滤波参数进行寻优，使算法对整个参数空间进行高效搜索以获得最优解。相比其他方法，粒子滤波可以得到更优化的回波运动矢量场。粒子滤波利用一系列带权值的空间随机采样，来逼近后验概率密度函数。

1.4.2 质量控制—双边滤波

1998年，Overton等[31]最先提出了双边滤波（bilateral filtering）的概念，它是在基于空间分布的高斯滤波函数基础上提出的。传统高斯滤波仅仅关注像素的空间距离而忽略了像素值的变化，因此会在滤除噪声的同时造成边缘的模糊。双边滤波在传统高斯滤波器的基础上添加了考虑像素值变化程度的权重。因此，在滤除图像噪声的同时考虑到灰度相似度的信息，使得权重系数随着图像灰度的变化而改变，有效地保持了图像边缘结构，最终能够进行自适应的滤波。

研究选取2016年4月13日飑线过程的回波进行滤波试验。从滤波前的回波（图5a）可以看出，强回波的边界、强回波区的对流回波及部分层状云回波比较零碎、杂乱，部分回波呈锯齿状。从滤波后回波对比可以看出，中值滤波和双边滤波后（图5b和5c）回波的形状、强回波中心以及强回波后的层状云回波和实况大致一样，强回波的边缘、强回波中心，后部的层状云回波变得更平滑、有序。对比滤波后飑线前沿回波，中值滤波后飑线前沿的回波被削弱，而双边滤波较好保留了前沿回波。中值滤波可以改善回波的质量，但是会削弱回波的边沿，尤其是飑线的边沿，在进行回波运动矢量场计算时，边沿的风矢量场只能由估计分析得来，从而影响运动估计的质量。而双边滤波在提高回波质量的同时，也很好地保留了回波边沿，在进行运动矢量计算时可以直接获得该处的风矢量场，质量相对来说更可靠。

图5 2016年4月13日雷达回波滤波前后对比：（a）滤波前；（b）中值滤波后；（c）双边滤波后 Fig. 5 Comparison between the Guangdong multi-radar mosaic on 13 April 2016: (a) before filtering; (b) after median filtering; (c) after bilateral filtering

1.4.3 运动矢量场对比

获得真实、平滑的运动矢量场是进行外推预报的关键。分别用交叉相关法、光流法、粒子滤波法追踪2016年4月13日06：30飑线过程回波的运动矢量场，其中交叉相关法和光流法均采用中值滤波法对雷达基数据进行质量控制，粒子滤波法采用双边滤波进行质量控制。

该飑线在06：30—07：30的1 h内大致呈西北—东南方向移动，粒子滤波融合法预报的移动路径与实况大致相同。在实况中飑线的不同位置移动方向有所差别，飑线南段向偏东南方向移动，北段向偏东移方向移动，东南段介于两者之间。粒子滤波融合算法追踪的飑线运动矢量风场（图6a）中，飑线南段大致为西北气流控制，北段偏西气流控制为主，东南段为西北气流控制，与飑线相应位置的移动方向基本一致，对矢量场的刻画也更精细。交叉相关法追踪的矢量场大致为偏西风（图6b），预报了飑线在1 h内朝偏东方向移动。而光流法追踪的飑线南段为偏西风，北段为西南风（图6c），预报飑线在1 h内朝偏东北方向移动。交叉相关法和光流法在进行运动矢量追踪时，由于质量控制的原因，强回波中心或者旋转处运动矢量会被周边的运动矢量所平滑，导致预报的移动路径出现偏差。由此可知，粒子滤波融合法获得了更接近实况、更真实的、更精细的回波运动矢量。

图6 2016年4月13日06：30运动矢量场对比：（a）粒子滤波融合法；（b）交叉相关法；（c）光流法 Fig. 6 Comparison between motion vector field at 06:30 BT 13 April 2016: (a) particle filtering fusion algorithm;(b) cross-correlation approach; (c) optical flow method

检验结果表明[19]，粒子滤波法能对各中天气类型的回波进行较好的预报，尤其是飑线预报效果好于光流法和交叉相关法。

2 人工智能临近预报方法研究

2.1 人工智能临近预报方法研究进展

20世纪80年代开始，人工智能技术已经在气象领域得到探索性的应用，从气候资料的处理分析到预报产品的制作都有涉足，其中主要有[32]基于人工神经网络、基于遗传计算及其融合算法、基于支持向量机方法、基于贝叶斯方法等气象预报方法。由于计算能力、数据存储分析能力、以及人工智能算法等基础条件的限制，尽管当时的一些成果已显示出较好的应用前景，人工智能技术在气象领域的进展仍受到了制约[20]。

面对信息技术和智能预报技术的快速发展，国内气象界也在积极探索将智能预报技术与气象业务科技的发展相结合，并开展了积极的尝试。如中央气象台、福建省气象台等气象部门通过合作开展基于卷积神经网络、决策树模型等算法的智能临近预报方法研究。香港天文台从2014年开始和香港科技大学合作开发基于卷积神经网络预测模型的深度学习降水临近预报方法[33]。通过机器深度学习，模拟雷达回波未来2 h的移动路径，较传统基于光流法预测雷达回波移动具有了新的优势。

2.2 基于生成对抗网络GAN 的智能预报技术

基于生成对抗网络GAN的智能预报技术是从一系列雷达观测资料中，运用图像识别等方法提取知识和信息建立预测模型，并通过损失函数训练模型，输出产品，进行预报的人工智能临近预报新技术。该方法可快速和智能化地识别回波，从而预报未来一段时间内对流天气发生和影响区域（图7）。

图7 基于GAN的人工智能临近预报方法示意图 Fig. 7 Artificial intelligent nowcasting approach based on GAN

GAN由谷歌大脑团队的Goodfellow等[22]提出，近几年在学术界受到广泛的关注，但目前研究主要集中在传统的视觉领域，即利用GAN模型生成传统视觉领域图片，如鸟、猫、狗等。深圳市气象局首次使用GAN方法来做回波临近预报。

2.2.1建 模模型时的，建需立对及原训始练雷达图片进行5次卷积[34]。其中第一到第四次卷积分别对输入图片采用卷积层和纠正线性函数ReLU[35]进行卷积，分别获取图片特征，第五次卷积采用卷积层和双曲正切函数[36]进行卷积，获得第五特征图，再对5次获取的卷积进行堆叠完成模型。建模应用了2015—2017年共3年的雷达数据。

通过堆叠5个卷积层从原始图像中提取的特征，建立预报模型，输出预报目标，通过计算预报目标与真实目标图像之间的损失函数，再通过损失函数对模型进行训练，促使模型向正确的方向优化。

严格来说，一个GAN框架，最少拥有两个组成部分，一个是生成器G（Generator），一个是判别器D（Discriminator）。生成器G基于模型生成预报回波图像序列。判别器则需要对生成器生成的图片进行真假判别。在训练过程中，会把生成器生成的样本和真实样本信息传递给判别器D。判别器D的目标是尽可能正确地识别出真实样本，和尽可能正确地揪出生成的样本，也就是假样本。而生成器的目标则和判别器相反，就是尽可能最小化判别器揪出它的概率。

生成器G 和判别器D 组成一个最小最大游戏（min-max game），在训练过程中双方都不断优化自己，直到达到平衡——双方都无法变得更好，也就是假样本与真样本完全不可区分，从而完成模型的优化训练。

GAN方法本身没有明确的“预报”过程。在模型的训练过程中，使用历史序列雷达图片和1 h内预测目标的数据对作为训练样本。这样模型在训练中学到的输出就是未来1 h的预报。对抗学习模型见图8。

图8 基于GAN的智能临近预报对抗学习外推模型 Fig. 8 Adversarial learning model based on GAN

2.2.2 模型的效果

受南海西北部热带低压和西南气流的共同影响，2018年6月5日广东大部分地区出现了暴雨降水。降水主要出现在5日上午后。低压外围的对流云团不断涌现，给广东带来持续的间歇性降水。

选取当日12时的回波进行60 min的外推预报试验。对比这次降水过程雷达反射率因子拼图13时的实况（图9c）和智能预报方法在12时起报的13时外推预报（图9b），智能预报方法预报出了主要的四块回波：粤西一带的东北—西南向的条状回波、珠江口西侧的偏西北—东南向的条状回波、珠江口东侧的回波及粤东沿海的条状回波。预报的回波位置与实况大致相当，但强度偏弱，强回波的范围偏大。

图9 2018年6月5日反射率因子拼图实况和预报：（a）16：12实况；（b）16：12的60 min预报；（c）17：12实况 Fig. 9 Comparison between the Guangdong multi-radar mosaic on June 5, 2018: (a) real radar echo at 16:12 BT;(b) 60 min forecast at 16:12 BT; (c) real radar echo at 17:12 BT

相对于天气过程的非线性、多因子作用的特点，只用雷达数据做预报显然是不够的。在进行模型训练时，除了雷达数据，还应该包含大气动力信息、热力信息，甚至地形信息等，以提高模型的适应能力，提高智能预报水平。对算法的改进，可以考虑在模型中加入Dense Net结构、Batch Normalization批归一化层等，以提高模型的性能。

3 预报效果评估

多年临近预报方法的研究，极大提升了深圳市气象局临近预报预警水平。对流天气预警时间提前量由30 min左右提高到2018年的56 min。为了对临近预报效果进行定量检验，分别计算粒子滤波法、光流法和GAN智能临近预报方法的击中率（POD）、空报率（FAR）、临界成功指数（CSI），具体算法见文献[16]。选取2018年的18个降水过程，包括西风带、西南季风、台风降水等，将回波强度分为25～80，35～80，45～80 dBz三个等级。结果如图10。

图10 粒子滤波法、光流法及GAN方法1 h预报试验结果击中率、空报率和临界成功指数对比 Fig. 10 Comparation of hit rate, false forecasting rate and critical success index (CSI) in 1 hour later forecasting experiment via particle filtering method, optical flow and GAN

由图10可知，三种方法对25 dBz以上中等强度回波的POD都在55%以上，其中光流法最高，GAN方法次之，略低于光流法。而成功临近指数CSI中，光流法和GAN方法大致相当。由此可知，总体优势上，GAN和传统的光流法大致相当。对光流法，三种强度回波的预报效果均好于GAN方法和粒子滤波法，目前业务应用以光流法为主。但对不同的过程，三种方法预报效果不尽相同，光流法对局地生成及强度随时间快速变化的回波预报效果较好，粒子滤波法对飑线的预报效果较好，而GAN方法对飑线和台风的预报也具有较大的参考价值。

4 结论和讨论

近些年来，深圳市气象局逐渐开发了基于运动矢量场的交叉相关法、光流法、粒子滤波融合法临近预报技术，并与哈尔滨工业大学深圳研究生院合作开展了基于生成对抗网络（GAN）的人工智能临近预报方法研究。

基于运动矢量场的临近预报方法中，三种预报方法对于不同天气过程各有优势，业务中很难有一种方法能很好地做出所有天气过程的临近预报，对于不同的天气过程，只能通过预报员的判断来选择采用哪种预报方法。基于生成对抗网络（GAN）的人工智能临近预报方法大致可以预报1 h的回波，预报结果有一定的参考价值。

目前，基于雷达回波的外推临近预报技术是基于雷达反射率因子单个因子，并没有考虑影响风暴生消的动力学、热力学等因素，几乎没有能力预报回波强度的变化，对于生命史小于预报时效的雷暴，预报意义不大。而对于风暴生命史较长的过程，比如超级单体风暴、持续几个小时的强飑线、由锋面造成的强降水雨带以及台风等大尺度天气过程，外推临近预报具有一定的现实意义。

对于人工智能临近预报方法，灾害性天气事件如冰雹、龙卷、台风等的气象观测数据目前尚未达到深度学习需要的样本数量，所以基于深度学习的智能临近预报方法的研究是一个持续的过程，需要持续的投入研究。

对于临近预报未来发展主要依赖以下几方面：一是数值模式的发展，将先进的外推预报方法同快速更新循环的高时空分辨率数值模式预报以及二者的融合，是未来强对流天气短时临近预报的重要发展方向[13]；二是大数据的挖掘，发展基于人工智能的临近预报方法，希望可以寻找到一些规则并应用到业务服务中；三是风暴追踪算法的不断改进；四是预报员能力的逐渐提高。