2. 江西省气象台 江西省气象局天气预报开放实验室, 南昌 330096
2. Jiangxi Meteorological Observatory, Weather Forecast open laboratory of Jiangxi Meteorological Bureau, Nanchang 330096
强对流天气是导致气象灾害的重要天气类型,与中小尺度天气系统关系密切(俞小鼎和郑永光,2020),天气雷达则是探测中小尺度天气系统的有力工具(罗义等,2021),在雷暴大风(朱平和俞小鼎,2019;黄小彦,2020;杨晓亮等,2020)、冰雹(何炳文等,2020;杨吉等,2020;曾勇等,2020)、暴雨(何钰等,2021;苏爱芳等,2021;张晓辉等,2022)、龙卷(周后福等,2014;黄先香等,2018;邱阳阳等,2021)等强对流天气的预报预警中发挥了巨大作用。近年来,随着人工智能技术的迅猛发展,采用机器学习、深度学习等对雷达回波进行智能识别、短临外推等方法越来越多。张文海和李磊(2019)通过贝叶斯分类法对正、负样本数据集进行机器学习,检验结果表明人工智能算法比传统概念模型法命中率高9%。顾建峰等(2020)采用TrajGRU技术(Zhang et al.,2021),初步建立了三维雷达回波智能临近预报模型,结合地面观测资料和U-NET技术(Ronneberger et al.,2015),对雷暴大风和冰雹进行智能识别。这些方法需依赖大量历史雷达图像数据集进行学习建模,才有可能使人工智能学习到历史数据多种特征,最终建立泛化性能较好的分类、回归模型。但如果历史雷达图像中存在较多的错误数据干扰,甚至一些仅有1个或数个像素点大小的杂波,都可能导致特征学习出现异常(Elsayed et al.,2018)。因此雷达数据集的质量至关重要。
雷达在实际业务观测中,经常会遇到一些非气象回波的干扰,如地物回波、超折射回波、同波长干扰回波、飞机船只等回波、海浪回波、由天线辐散特性造成的虚假回波等(张培昌等,2001)。这些虚假回波会严重影响历史观测数据集质量,因此有必要在使用之前对雷达回波进行质量控制,尽量剔除虚假回波的干扰。马中元等(2010)采用非雷达资料(云图、自记雨量和闪电数据)排除地物杂波,采用统计“突变率”、两层仰角叠加分析和改进型中值滤波算法去除数据杂波。仰美霖等(2018)等将5×5区域内的中心点周围有效值少于某一阈值时,将中心点赋为无效值;利用有效回波的连续性,记录总库数的最大、最小值,判断是否为条幅状干扰回波等。魏鸣等(2019)采用支持向量机进行地物杂波的识别,并与人工神经网络识别结果进行对比,发现支持向量机的效果更好。文浩等(2020)采用模糊逻辑方法对径向分布的条幅状干扰回波进行了识别,通过建立隶属函数和阈值,对径向干扰回波和降水回波的识别和消除效果较好。
以上雷达回波质量控制方法虽然可以有效消除多种非气象回波的干扰,但也存在较多的问题。首先,需要使用多种资料来协助剔除虚假回波,这需要所有观测数据全部到达后才能进行,会增加雷达回波图像到达用户手中的滞后时间;当其它资料有缺失时,质控效果也会有所降低。其次,采用统计或者机器学习方法需要提供已经区分好的历史数据进行学习和统计,而这些历史数据也需要事先进行人工整理。在人工智能数据集制作过程中,面对大量雷达反射率图像时(周康辉等,2021;程文聪等,2020),如何快速提取所需雷达回波并进行质量控制,也是一个问题。第三,之前的雷达反射率因子质量控制方法大多从基数据入手,采用时间、空间特征分析方法对雷达反射率因子进行质量判断,计算较为繁琐,耗时也长,考虑的因素也多。但实际工作中发现,预报员从单张雷达反射率因子回波图像中,用较短时间就能确定非气象回波、弱回波、降水回波等多种特征回波,表明单张雷达回波图像上确实包含可以区分这些回波的特征信息。而深度学习核心技术的卷积技术正是寻找图像特征最有利的工具之一。本文利用不同形状的卷积核,通过改进后的卷积技术,实现对任意雷达反射率因子图像的质量控制,能够根据需求快速剔除任意雷达反射率因子图像中多种形状、大小的非气象回波和弱回波。该方法经过在江西省气象台长期运行,可以满足大部分业务和科研需求。
1 资料来源与研究方法 1.1 资料来源本文所用资料为南昌(115.91°E,28.67°N)多普勒雷达基本反射率因子图像,江西省及邻省雷达组合反射率因子拼图。其中基本反射率因子以1.5°仰角为例挑选具有多种回波图像进行质量控制说明。
1.2 研究方法本文采用自定义卷积核和改进卷积方法来对雷达反射率拼图进行质量控制。常规卷积核在图像处理上有较广泛的用途,Prewitt算子(Seung-Ju et al.,2009)、Sobel算子(Qu et al.,2005)等一系列滤波算子,很早就被用于边缘检测等工作中,也被称为Filter,其主要计算公式如下
$ (f * g)(u, v)=\sum\limits_i \sum\limits_j f(i, j) g(u-i, v-j)=\sum\limits_i \sum\limits_j a_{i, j} b_{u-i, v-j} $ | (1) |
其中,f(i,j)为卷积核,i、j分别为卷积核x和y坐标方向上的格点数量;g(u,v)为输入的2D图像,u、v分别为图像x和y坐标方向上的格点数量,*指卷积核在2D图像上的卷积操作。其中u>i, v>j, a是卷积核f格点上的值,b是图像g格点上的值。
有时为了卷积后的图像与之前的图像尺寸大小保持一致,需要将输入的2D图像用0值对图像进行填充(Padding)操作(图 1a)。因此公式(1)更改为
$ (f * g)(U, V)=\sum\limits_i \sum\limits_j f(i, j) g(U-i, V-j)=\sum\limits_i \sum\limits_j a_{i, j} b_{U-i, V-j} $ | (2) |
其中,U和V分别是进行填充操作以后的图像长度和高度,分别为U=u+i-1,V=v+j-1,带入公式(2)可得
$ (f * g)(U, V)=\sum\limits_i \sum\limits_j f(i, j) g(u-1, v-1)=\sum\limits_i \sum\limits_j a_{i, j} b_{u-1, v-1} $ | (3) |
从公式(3)可见,原始图像经过0值填充后(图 1a),再利用卷积核(图 1b)进行卷积,可以得到一个新的图像,且图像大小与原输入图像保持一致(图 1c)。该新图像每个格点的值为原始图像格点与卷积核卷积后的格点值之和。因此,a×b后的值再进行相加,可以有效放大不同特征之间的区别,有利于卷积操作中寻找不同特征,这也是卷积操作在深度学习中的核心地位原因之一。本文仅需要剔除异常格点值,保留正常值,因此不需要执行a×b后再相加这一步骤。设原图像某格点值为b1,则该格点卷积后的值p为
$ p= \begin{cases}b_1, & z>t \\ 0, & z \leq t\end{cases} $ | (4) |
其中,t为根据不同卷积核设定的阈值,在本文中是雷达反射率因子值,根据不同卷积核的作用分别给定相应的阈值t。计算步骤为:首先根据公式(3),从图像左上角开始与卷积核进行
在雷达反射率因子回波上,非气象类回波经常呈现孤立点状和射线状(张培昌等,2001)。在2016年4月15日01时26分(世界时,下同)南昌雷达1.5°仰角基本反射率图(图 2)中,A1、A2、A3均为孤立点状回波,其中A1由大量孤立点状回波构成的孤立点回波群,包含一个或多个像素点组成的离散孤立点回波;A2是较为典型的孤立点状回波,虽然由数个像素组成,但并不是真正的天气回波;A3面积较大,严格来说已经不再是孤立点状回波,本文称为准孤立点状回波。射线状回波B也是雷达反射率因子回波中经常见到的非气象类回波,在雷达信号处理器异常,信号受到外部电磁干扰时,雷达会观测到沿径向分布的条幅状干扰回波(文浩等,2020)。地物回波也是较为常见的非气象类回波,且这种回波经常覆盖范围较大,无论是在平常预报业务当中,还是在制作雷达回波数据集时,都容易造成误判,是必须要剔除的杂波,详见2.3节。此外,如果研究对象仅关注较强回波,那弱回波C也需要剔除。
对于1个像素点大小的孤立回波,可以采用图 3a的卷积核(以下简称孤立卷积核Ⅰ);对于9个像素点以内(长、宽均不超过3个像素)的孤立回波,可以采用图 3b的卷积核(以下简称孤立卷积核Ⅱ)进行卷积。经过测试,孤立卷积核Ⅰ的阈值取10 dBz,孤立卷积核Ⅱ的阈值取30 dBz,消除孤立点状回波效果最好。从图 3b中也可知,如果孤立回波的长或宽超过了3个像素,采用孤立卷积核Ⅱ是无法消除的,需要经过多次卷积操作,才可能消除,甚至无法完全消除。长、宽均超过3个像素的像素群则无法被消除。长、宽均超过5个像素的像素群已经不再是孤立点状回波,也不能被这2个卷积核给消除。利用孤立卷积核Ⅰ对图 4a进行卷积后可以看到,A1孤立回波群中1个像素大小的回波全部被剔除(图 4b),A2、A3及左下方的高回波区均得以保留。利用孤立卷积核Ⅱ进行卷积后发现,A1孤立回波群中更多的多像素孤立回波被剔除;A2孤立点状回波也被剔除了4个像素点,仅剩3个像素点(图 4c)。对图 4c用孤立卷积核Ⅱ再进行一次卷积,A1的零散孤立点状回波进一步被剔除,A2则完全被剔除(图 4d)。几次操作中可以看到,虽然A3的面积比A2略大一些,但没有受到卷积核卷积的影响,仍保持相同大小,证明这两种卷积核造成的图像信息丢失程度较低。
上述两种卷积核能够较好地解决孤立点回波的剔除,但对其他较大面积的线状回波则无法胜任,如果需要剔除的虚假回波形状较为复杂,例如射线状虚假回波,基本上无法剔除。为此,需要针对线状回波等线性非气象类回波设计新的卷积核方案。线状非气象类回波在雷达反射率图上经常表现为一条等宽(图 5a—d)或射线状的线状回波(图 5e、f),其中等宽线状回波的长或宽不会有明显变化,射线状回波越远离图像中央,线宽越大。图 5a、b分别是水平和垂直等宽线状回波,图 5c、d分别是不同方向的倾斜等宽线状回波,图 5e、f分别是不同方向的射线状回波。
针对图 5a的水平线状回波,可以采用上下为1,中间为0的卷积核来进行卷积操作(图 6a);针对图 5b的垂直线状回波,可以采用左右为1,中间为0的卷积核来进行卷积操作(图 6b)。图 5c、d都是倾斜线状回波,仅方向不一样,需要分别制作一个卷积核(图 6c、d),来剔除不同方向的等宽线状回波。以上4种卷积核可以有效剔除横截面积不超过3个像素点的线性回波。由于回波图中经常同时存在水平、垂直线状回波,还有两个方向的倾斜线状回波,因此经常两两混合使用。
在实际回波图中应用可以看到,在利用水平、垂直卷积核对图 7a进行卷积后(图 7b),A1孤立回波群的非气象回波大量减少,仅剩余4个较大孤立回波团,A2回波面积大幅度减少,表明水平、垂直卷积核较好的对这种孤立小回波进行了质控。利用倾斜卷积核卷积后的回波图与图 7b类似(图 7c),在A1孤立回波群中剩余了4个较大的回波团,A2回波也完全消失。两种卷积方案结果相似的原因在于,倾斜线状回波也可以分解为水平和垂直方向的水平线状回波。但在回波复杂时,用倾斜卷积核效果比水平、垂直卷积核效果要好些。将水平、垂直卷积核和2种倾斜卷积核先后卷积后可以发现(图 7d),A1孤立回波群中最左侧的较大回波团基本被剔除,剩余三个回波团的面积也略有减少。这表明经过多次卷积后,部分非气象类杂波可以逐渐被剔除,同时也可能会增加图像信息的丢失。
利用图 6a—d的卷积核对射线状回波进行卷积后发现,如果射线横截面没超过3个像素,则射线会被剔除。但随着射线远离图像中心,宽度会逐渐增加,一旦横截面超过3个像素,则不能再被剔除(图 8b)。如果将图 6a—d的卷积核从5×5扩大为7×7,0值区扩大为5×5,虽然可以剔除更多的射线状回波,但也可能会剔除更多的正确回波。为此,本文采用图 6e、f的射线卷积核来对射线状回波进行剔除。虽然卷积核大小仍然是5×5,但倾斜的0值区域按照3-4-5-4-3交错排列,使得该卷积核能够对横截面为5像素的射线状回波进行剔除。虽然射线状卷积核与前4种线状卷积核形状相似,但这种卷积核的阈值要设置的较低才有效果。如本文种射线卷积核的阈值仅为3 dBz,最大不超过5 dBz。只有这样才不会出现将正常成片回波错误剔除。如果射线横截面宽度达到了7像素大小,可以采用图 6g—h的射线状卷积核来处理。极端条件下,可以采用图 6i的射线状卷积核进行卷积,该卷积核能处理横截宽度达到11个像素的射线,但图像信息丢失程度可能也会相对较高。
图 8a中的射线状回波在经过2个不同方向的倾斜等宽卷积核先后卷积后,在图像边缘仍然存在较长的残存回波(图 8b);同样,经过水平、垂直和2个倾斜等宽卷积核先后卷积,也无法彻底解决射线状回波,在图像边缘残存了一段较短的回波(图 8c)。但利用两个射线卷积核卷积进行卷积后,射线状回波已经消失,A1孤立回波群、A2孤立点状回波也进行了剔除(图 8d)。因此,射线状回波可以同时完成射线状回波和孤立点回波的剔除,且效果也较好,消除弱回波区域比孤立点卷积核消除的面积稍大,图像信息丢失程度中等。
经过实验发现,图 3和图 6的卷积核均无法剔除地物回波及其他一些弱回波,图 8b—d也表明,图 2的弱回波C和孤立回波A并没有完全消除。实际上,这些卷积核也不能剔除地物回波,这是因为这些弱回波区的回波都是成片连续的,容易超过阈值导致无法剔除,必须再设计新的卷积核来处理。为此,本文设计了地物回波和弱回波卷积核(图 9),能有效剔除地物回波及其他弱回波。
当采用图 9的弱回波卷积核时,阈值设置为40 dBz,即斜线方向上3个格点的dBz值之和不超过40 dBz,就可以将中间格点值赋为0 dBz,否则仍然取原值。如果设置为30 dBz,弱回波保留的会略多一些;如果设置为50 dBz,弱回波保留的会略少一些,可以根据需要自行调整。利用弱回波卷积核对图 10a进行卷积,可以看到绝大部分弱回波区域、孤立回波区域和射线状回波全部被剔除,图像基本上只保留了20 dBz以上的较强回波区域(图 10b);在A1和A3区域还残留着很少的几个孤立像素,再用孤立点卷积核或线状卷积核卷积一遍即可全部去除。图 10c是2016年4月15日01时26分江西及周边省份多部雷达组合反射率因子回波拼图图像,可以看到在图像北部有大量晴空湍流回波,其面积比南部的较强回波面积还要大。这些晴空湍流回波无法用孤立点回波和线状回波进行消除,但采用弱回波卷积核进行卷积后可以看到,几乎所有的晴空湍流回波全部被剔除,仅保留了南部较强回波部分,质量控制效果明显(图 10d)。与前几个卷积核相比,该卷积核的卷积计算丢失的信息相对是最多的。
表 1是本文所有卷积核的名称、作用和建议阈值,建议阈值是利用江西8部雷达5 a反射率图像进行数据质量控制后经过统计、均衡卷积效果后得到(图 6i卷积核出现次数极少,故没有统计在内)。可以利用多种卷积核进行多次卷积,以消除局地孤立回波(建议先进行水平、垂直、倾斜射线、弱回波等卷积操作,再利用孤立卷积核消除残余的孤立像素点)。在实际操作中,可以根据需要,构造其他形态的卷积核核阈值,以满足不同要求。
本文基于卷积计算原理,对常规卷积方法进行了改进,构造孤立点状卷积核、线状卷积核和弱回波卷积核,可以实现对雷达单站回波、组合反射率拼图中的孤立点回波、线状回波、地物回波等非气象类回波和弱回波区域进行剔除。主要结论如下:
(1) 定义改进卷积方法,即对图像上的某一个格点进行卷积操作
(2) 构造孤立卷积核Ⅰ和孤立卷积核Ⅱ,可以有效剔除1个像素和9个像素(长、宽均不超过3个像素)内的孤立回波;构造水平、垂直卷积核、2个方向的倾斜卷积核,可以有效剔除孤立点回波、水平方向、垂直方向和倾斜方向的线状回波(横截面宽度不超过3个像素);构造2个方向的射线卷积核Ⅰ和射线卷积核Ⅱ,可以分别有效剔除横截面宽度不超过5和7个像素的射线状回波;射线卷积核Ⅲ则可以剔除横截面宽度不超过11像素的射线状回波;构造弱回波卷积核,可以有效剔除地物回波和弱回波。
(3) 不同卷积核的回波剔除能力不一样,图像信息丢失程度也有较大区别,可以根据需要使用一种或多种卷积核相互配合进行卷积以剔除不需要的回波。同时卷积核构造简单,使用方便,可以非常容易自定义新的卷积核,以便对更复杂的雷达反射率因子回波进行质量控制操作。
本文使用的改进卷积方法仍然有一定的不足,该方法只能剔除多余的错误回波,对缺失(例如遮挡区域)的回波不能补齐,也不能判断弱回波是不是真实回波,例如利用弱回波卷积核进行卷积时可能会错误地剔除阵风锋等一些有意义的弱回波区,因此采用这种方法要充分考虑图像信息丢失的接受程度。若对弱回波有需求,则尽量采用非弱回波卷积核进行质量控制;若只关注图像中的强回波形态变化,则可以采用弱回波卷积核进行质量控制。
通过对江西8部雷达5 a反射率图像的计算和统计,基本确立了在预报业务和科研工作中的质控方案,即在日常预报业务中,为了减少图像信息的丢失程度,主要采取多个非弱回波卷积核组合卷积的方式对雷达回波进行质量控制,这种应用场景下要求剔除非气象类杂波的时候尽量减少对正常回波的影响。例如最常用的卷积组合为:“孤立卷积核Ⅱ”+ 2个“倾斜卷积核”+“垂直卷积核”+“水平卷积核”+“孤立卷积核Ⅰ”+“孤立卷积核Ⅱ”,并至少执行2遍卷积操作,即可在最大限度保留图像信息的同时,消除部分非气象杂波,效果与图 8c类似。在人工智能深度学习等科研工作中,科研人员通过剔除弱反射率回波,让深度学习模型更关注短时强降水等强对流回波的演变特征,从而提高短时强降水等强对流天气的预报准确率。通过同一个基于卷积长短期记忆神经网络方法进行学习预测,经过图像质控的数据集性能比未经图像质控的应用效果更好。
程文聪, 史小康, 张文军, 等. 2020. 基于深度学习的数值模式降水产品降尺度方法[J]. 热带气象学报, 36(3): 307-316. |
顾建峰, 周国兵, 刘伯骏, 等. 2020. 人工智能技术在重庆临近预报业务中的初步研究与应用[J]. 气象, 46(10): 1286-1296. |
何炳文, 胡振菊, 高伟, 等. 2020. 湘西北地区强冰雹的多普勒天气雷达旁瓣回波统计分析[J]. 暴雨灾害, 39(3): 269-275. |
何钰, 陈小华, 李耀孙, 等. 2021. 云南省副热带高压外围类短时强降水的雷达回波特征[J]. 气象, 47(4): 450-462. |
黄先香, 俞小鼎, 炎利军, 等. 2018. 广东两次台风龙卷的环境背景和雷达回波对比[J]. 应用气象学报, 29(1): 70-83. |
黄小彦, 孙继松, 刘文婷. 2020. 地形作用下低空急流的演变与强降水对流风暴系统的相互作用[J]. 气象学报, 78(4): 551-567. |
罗义, 梁旭东, 王刚, 等. 2021. 雷达反演的多尺度风场在临近预报中的应用研究[J]. 暴雨灾害, 40(4): 401-409. |
马中元, 朱春巧, 刘熙明, 等. 2010. CINRAD雷达数据质量控制方法初探[J]. 气象, 36(8): 134-141. |
邱阳阳, 程向阳, 杨祖祥, 等. 2021. 登陆台风"温比亚"引发的龙卷过程分析[J]. 暴雨灾害, 40(5): 531-540. |
苏爱芳, 吕晓娜, 崔丽曼, 等. 2021. 郑州"7.20"极端暴雨天气的基本观测分析[J]. 暴雨灾害, 40(5): 445-454. |
魏鸣, 管理, 梁学伟, 等. 2019. 基于支持向量机的雷达地物回波识别研究[J]. 大气科学学报, 42(4): 631-640. |
文浩, 张乐坚, 梁海河, 等. 2020. 基于模糊逻辑的新一代天气雷达径向干扰回波识别算法[J]. 气象学报, 78(1): 116-127. |
杨吉, 郑媛媛, 徐芬. 2020. 江淮地区一次冰雹过程的双线偏振雷达观测分析[J]. 气象学报, 78(4): 568-579. |
仰美霖, 江源, 刘黎平, 等. 2018. 北京SA雷达电磁干扰回波特征及质控算法初探[J]. 干旱气象, 36(5): 802-812. |
杨晓亮, 杨敏, 隆璘雪, 等. 2020. 冷涡背景下河北雷暴大风环境条件与对流风暴演变个例分析[J]. 暴雨灾害, 39(1): 52-62. |
俞小鼎, 郑永光. 2020. 中国当代强对流天气研究与业务进展[J]. 气象学报, 78(3): 391-418. |
曾勇, 万雪丽, 李丽丽, 等. 2020. 一次多单体冰雹天气过程的雷达回波与闪电特征分析[J]. 暴雨灾害, 39(3): 250-258. |
张培昌, 杜秉玉, 戴铁丕, 等. 2001. 雷达气象学(第二版)[M]. 北京: 气象出版社.
|
张文海, 李磊. 2019. 人工智能在冰雹识别及临近预报中的初步应用[J]. 气象学报, 77(2): 82-291. |
张晓辉, 高艳春, 易永力, 等. 2022. 承德市冰雹时空分布特征及雷达预警指标分析[J]. 沙漠与绿洲气象, 16(2): 16-22. |
周后福, 施丹平, 刁秀广, 等. 2014. 2013年7月7日苏皖龙卷环境场与雷达特征分析[J]. 干旱气象, 9(3): 415-423. |
周康辉, 郑永光, 王婷波. 2021. 利用深度学习融合NWP和多源观测数据的闪电落区短时预报方法[J]. 气象学报, 79(1): 1-14. |
朱平, 俞小鼎. 2019. 青藏高原东北部一次罕见强对流天气的中小尺度系统特征分析[J]. 高原气象, 38(1): 1-13. |
Elsayed G F, Shankar S, Cheung B. 2018. et al. Adversarial Examples that Fool both Human and Computer Vision [J]. DOI: 10.48550/arXiv.1802.08195
|
Qu Y D, Cui C S, Chen S B, et al. 2005. A fast subpixel edge detection method using Sobel-Zernike moments operator. Image Vis Comput[J]. Image & Vision Computing, 23(1): 11-17. DOI:10.1016/j.imavis.2004.07.003 |
Ronneberger O, Fischer P, Brox T. 2015. U-net: Convolutional networks for biomedical imagesegmentation [C]//Proc. 18th Int. Conf. on Medical Image Computing and Computer Assisted Intervention, Munich, Germany, Springer. DOI: 10.1007/978-3-662-54345-0_3
|
Seung J J, Hyung B K, Sang B Y. 2009. Implementation of image filter for computer graphic software[J]. International journal of computer science and network security, 9(3): 79-86. |
Zhang L, Huang Z, Liu W, et al. 2021. Weather radar echo prediction method based on convolution neural network and long short-term memory networks for sustainable e-Agriculture[J]. Journal of Cleaner Production, 298: 126776. DOI:10.1016/j.jclepro.2021.126776 |