一种基于阵风锋时空变化特征的大风识别算法研究

吴奕霄，王珂清,程 婷，项 瑛，周学东，徐 芬

（1.江苏省气候中心，江苏南京 210041；2.江苏省人工影响天气中心，江苏南京 210041；3.江苏省气象信息中心，江苏南京 210041；4.南京气象科技创新研究院，江苏南京 210041）

0 引言

阵风锋过境伴有气压升高、风速突增和温度陡降等现象，阵风锋及其母雷暴经常会产生地面大风，引发人员伤亡或者经济损失，因此阵风锋研究对防灾、减灾有重要意义［1］。提前识别和预警阵风锋有利于采取有效的措施避免危险事故发生，进而降低和减少经济损失和人员伤亡。雷达目标检测能够判断回波信号中目标存在与否，并提取目标位置信息［2-3］。因此，利用雷达进行易引发大风的阵风锋自动检测技术研究引起了广泛关注。模板匹配技术是在设定一个目标模板后，与具有相似目标的图像进行比对，以获取目标模板在图像中的位置，以此达到图像识别的目的，其优点是匹配的算法更加简单，在图像变化不大的情况下，识别率理想［4-5］。

在天气雷达反射率因子图像上阵风锋表现为一条高度较低（低仰角层观测）、形态呈弧线结构和具有移动特点这3 个特征的窄带回波。在雷达径向速度数据场上呈现速度辐合特征。Uyeda等［6］首先提出阵风锋自动探测算法（Automatic Detection of Gust Fronts,ADGF），该算法只使用沿雷达径向的速度辐合，将径向切变超过预设阈值的段组合在一起，用曲线拟合出每个段上的最大切变值，这条曲线就代表探测到的阵风锋。后来提出了改进的阵风锋探测算法（Gust Front Detection Algorithm,GFDA）［7-9］，利用两个低仰角垂直关联减少了虚假预警，该算法被用于美国TDWR 雷达业务系统。美国强风暴实验室和麻省理工学院林肯实验室联合开发了改进的阵风锋自动探测算法［10-11］，包含了阵风锋雷达回波细线探测，通过寻找局地最大反射率因子与阵风锋反射率因子和宽度一致性探测阵风锋回波细线。Delanoy等［12］开发了机器智能识别阵风锋算法（Machine Intelligent Gust Front Algorithm,MIGFA），该算法首次使用模板函数计算雷达回波图和速度标准偏差得分，并设计了相应的质量控制算法抑制虚假信号，最后在组合得分图像上提取阵风锋。随着需求的提升和技术的发展，一些新的或改进的阵风锋识别技术也得到全面发展，徐芬等［13］通过统计江苏3种窄带回波分布特征和径向波形特征，在此基础上设计了基于回波平坦性计算和阵风锋径向波形寻找算法的阵风锋自动识别算法。冷亮等［14］基于阵风锋的3 个特征，提出了基于数字图像处理技术中的数学形态学算法的阵风锋自动识别算法，通过建立雷达观测范围内的二维特征量场，对其处理后得到一条反映阵风锋骨架特征的曲线。徐芬等［15］改进和调整了MIGFA 算法，设计了较高仰角（1.5°/2.4°）反射率的阵风锋细线函数模板，替代原空间差分反射率函数模板，新的算法对于CINRAD/SA 雷达探测到的江苏沿江地区阵风锋回波的识别更有针对性。HWANG 等［16］使用模糊逻辑理论开发了新的神经模糊阵风锋检测算法（Neuro-Fuzzy Gust Front Detection Algorithm,NFGDA），该算法通过使用神经网络的训练进行优化，对于阵风锋的识别率高达92%。Yuan等［17］基于统计特征分析和图像处理技术，提出一种基于本地二进制双模板（Local Binary with Dual-Template,LBDT）算法改进的雷达数据自动检测阵风锋的方法。该算法对反射率值和阵风锋宽度的变化不敏感，在自动识别阵风锋过程中检测概率高、误报率低。徐月飞等［18］基于Faster RCNN 算法和InceptionV2 模型，提出了改进的深度卷积网络阵风锋识别算法（d-CNNGFDA），利用窄带弱回波的绕雷达中心旋转不变性的特性，只对垂直走向的窄带弱回波进行标记，扩充了数据量用以提升复杂回波中阵风锋引起的窄带弱回波的自动识别能力。

上述这些方法均需选取阵风锋样本符合各自算法设定条件上才可以较好的识别，各有使用局限性。本文试图融合易引发大风的阵风锋通用特征:雷达回波图像所在高度较低、图像特征呈弧线结构和具有移动这3 个特征，设计能表征这3个特征的函数模板，再采用风暴回波特征抑制技术和其他弧线回波特征抑制技术以降低检测虚警率，最终实现更具普适性的一套阵风锋检测技术。

1 识别方法研究

1.1 算法流程概述

根据观测事实以及学者们总结［12］的易引发大风的阵风锋雷达反射率因子图像具有以下特征:1）阵风锋雷达反射率因子图像呈细长线状，其数值一般小于25 dBz；2）阵风锋图像高度较低；3）呈现移动的特点，如图1方框内弧线状图像所示。根据这3个时空变化特征，采用相应预处理方法提高雷达反射率因子图像质量后，借鉴模板匹配方法［12］得到3种不同的得分值结果。

图1 阵风锋天气雷达反射率因子图像（方框内弧线状图像）

融合3种时空变化特征匹配结果后，采用区域增长法提取阵风锋弧线图像结构，达到自动检测阵风锋的目的。基于阵风锋时空变化特征的大风识别算法流程如图2所示，由该图可见通过对包含有不同特征的雷达回波图像进行相应的模板匹配后，融合多组得分值，经过两次有效提取，最终实现对易引发大风的阵风锋雷达反射率因子图像的提取。

在利用函数模板提取阵风锋图像特征前，为提高图像提取成功率，对数据质量进行如下预处理:将0.5°和1.5°仰角的雷达反射率因子图像采用9 点平滑滤除孤立点、奇异点和回波图像中的空点。将这两个仰角的数据从极坐标下转换到直角坐标下，空间分辨率为1 km×1 km。后续基于函数模板和得分函数的计算均基于直角坐标系下的雷达回波图像。

1.2 函数模板设计

根据阵风锋图像在不同仰角、仰角间图像特征差异、同一仰角前后时刻图像差异等设计相应函数模板提取能够反映阵风锋的图像特征。根据阵风锋在空间变化上的特征设计了阵风锋线状函数模板；根据阵风锋在空间上线状差异特征分别设计了风暴边缘函数模板和放射状回波图像函数模板；根据阵风锋在时间变化上的特征设计了阵风锋差分反射率因子图像函数模板。以上几种函数模板和得分函数分别设计如下:由于阵风锋天气雷达反射率因子图像特征呈细长弧线结构，根据这一特征设计函数模板和得分函数如图2所示，利用相应的得分函数可得到高的阵风锋图像得分值，抑制其他图像得分值。该设计思路是依据尽管阵风锋图像长度和弯曲弧度不等，但将其分为若干段，每一段均近似为直线，因此图3（a）中的模板被设计为直线段主要基于假设阵风锋图像的一个段为近似直线。黑色和灰色格点代表参与计算的点，灰色代表模板中阵风锋回波图像上的点，对应图3（b）中的得分函数1。因为阵风锋回波一般在25 dBz 以下，所以给25 dBz 以下回波一个较高的得分值，特别是9 dBz 附近［10］，是阵风锋最常出现的反射率因子值。给大于35 dBz 的回波一个缓慢下降的得分值，因为阵风锋在某些情况下会与风暴边缘结合在一起，而阵风锋一般不会出现负值，所以给负的反射率因子值一个负得分值。黑色的点代表阵风锋回波两侧的点，对应图3（b）中的得分函数0，给得分函数0 低反射率因子一个高的分值，而给高反射率因子一个低得分值。需要说明的是，函数模板和得分函数的长宽是可调的，因为各地发生的阵风锋长度、宽度和反射率因子值可能不尽相同，笔者结合江苏地区阵风锋反射率因子值和相关文献确定的得分函数值仅供参考。

因为阵风锋回波图像出现的高度较低，一般为3 km 以下，可利用高低仰角层差分反射率因子来识别阵风锋，这样做的好处是排除较多其他回波干扰，特别是贯穿了高低层的线状对流风暴图像。将上下两个仰角的反射率因子图像经上述数据预处理方法质控后，上下仰角同一格点数据相减得到差分反射率因子图像（图4（a））。如图4（a）所示，图中大部分回波已被消除，阵风锋回波图像基本被保留。根据阵风锋差分反射率因子图像特征，设计图5中的函数模板和得分函数对差分反射率因子图像进行特征提取。与图3 阵风锋线状函数模板和得分函数对比可见，采用的函数模板不变，得分函数根据阵风锋差分反射率图像特征作了一定调整。同理，因为易引发大风的阵风锋具有移动的特征，所以可以利用前后时刻差分反射率因子图像特征来识别阵风锋。将前后连续两个时刻的反射率因子图像经上述数据预处理方法质控后，前后时刻同一格点数据相减得到阵风锋差分反射率因子图像（图4（b）），仍旧采用图5的函数模板和得分函数对前后时刻阵风锋差分反射率图像特征进行提取。

图4 阵风锋天气雷达反射率因子图像（方框内弧线状图像）

图5 阵风锋差分反射率因子图像函数模板

图4（b）中的前后时刻差分反射率因子排除了较多的回波，却保留了风暴边缘图像，该风暴边缘呈弧线状，很容易被误识别为阵风锋图像。本文设计风暴函数模板和相应得分函数用以削弱线状风暴边缘对识别结果的影响。如图6（a）所示，函数模板被设计为一个近似的圆形，图中的每个格点代表1 km 距离，该函数能够抑制线状回波图像得分值，突出呈片状的风暴和层云图像得分值。图6(b)为风暴边缘得分函数图。

图6 风暴边缘函数模板

另外，在图1 和图4 的图像中都受到一种放射状回波的干扰。如图4（a）第一象限中能看到一个沿雷达扫描方向的线状图像（图中箭头所示），这是由于放射状回波干扰所致。为了排除放射状回波干扰，设计如图7所示的函数模板和得分函数，图7（a）中每个格点代表1 km，得分函数根据放射状回波图像特征设计。函数1 代表在放射状回波图像上，函数0代表在放射状回波图像两侧。因为放射状回波干扰图像一般只会出现在雷达扫描方向，所以该模板仅在极坐标反射率图像上沿雷达扫描方向进行计算。再分别利用图3、图5 函数模板和得分函数处理每幅雷达反射率图像得到的得分值结果中减去使用该模板计算得到的结果，即可排除放射状回波干扰。

图7 放射状回波图像函数模板

为了让得分值结果便于提取，本文还设计了如图8所示的函数模板和得分函数对融合的得分值结果进行平滑处理。如图8（a）所示，以黑色圆点为中心点，灰色方格代表阵风锋得分值上的点，黑色方格代表阵风锋得分值两侧的点，模板中心点两侧设计了较多的灰色格点，目的是让阵风锋得分值更加连续。

图8 平滑图像函数模板

由于阵风锋图像的朝向是不确定的，在某个格点上，上述图3、图5、图8 中的函数模板以黑色小圆点为中心，让其以10°的间隔旋转180°，得到不同角度下模板匹配的得分值，最后在这个格点上只保留最高得分值。然后以1 个格点的步长移过所有格点，并在所有格点上重复以上步骤。旋转公式如下:

式中，θ代表旋转角度，x'和y'代表旋转后新坐标，x和y代表原坐标。

经过上述多个函数模板和得分函数处理后，对多组得分值进行融合处理，并对得分值结果作二值化处理，得分值大于0.5 赋值为1，小于0.5 则赋值为0。为提取属于阵风锋的图像，本文借鉴区域增长法［11］对二值化图像进行处理，得到的候选阵风锋图像如果长度小于6 km，则排除这个候选阵风锋。然后再采用一个有条件的区域增长法对初次提取的阵风锋进行处理:如果一个格点满足与原候选阵风锋包含的点相邻，得分值不低于0.2且该点得分值角度与原候选阵风锋平均角度相差不超过30°，则认为该点是候选阵风锋中的一个点，并入候选阵风锋中。如果候选阵风锋长度超过20 km，则认为是阵风锋，否则排除。

2 识别方法检测

以2012年5月16日盐城新一代多普勒天气雷达探测到的阵风锋图像为例（图1方框内）。如图1（a）中0.5°仰角雷达站东南方向的阵风锋图像被风暴边缘掩盖，图1（b）中1.5°仰角雷达站东南方向的阵风锋图像较明显，基本未受影响，但是一般阵风锋图像出现高度较低，抬高仰角至1.5°时，随着离开雷达站距离增加，探测高度也在增加，在雷达站东偏北约25°方向第二个距离圈处已探测不到阵风锋图像。可见0.5°雷达资料探测高度较低，在近距离处阵风锋图像易被地物回波所掩盖；1.5°仰角探测时，受地物回波影响较小，却由于远距离处探测高度较高而不能探测到阵风锋。由此可以看出，单仰角的雷达资料并不一定能探测到完整的阵风锋图像，这也是文中利用了两个仰角的图像进行阵风锋自动识别的原因。

采用本文图3 设计的函数模板和得分函数分别对0.5°仰角和1.5°仰角雷达图像进行处理后，用灰度值表达计算得到的得分值如图9（a）和图9（b）所示，图像中越高的灰度值代表得分值越高。将图9 与图1 实况图对比可以看出，阵风锋回波处明显得到了较高的得分值，说明模板函数匹配技术能够在阵风锋识别中发挥重要作用。图9（a）和图9（b）中部分非阵风锋图像也得到了较高的得分值，这也是阵风锋自动识别中最困难的问题所在——相似回波较多且难以区分，所以文中基于多种阵风锋特征来对其进行自动识别。图9（c）是由图9（a）和图9（b）综合后的组合得分值。综合两个仰角的探测优势，将得分值图像组合起来，图9（c）中得到较完整、连续的阵风锋得分值，但同时也组合了两个得分值图像中的干扰回波图像得分值。

图9 阵风锋线状函数模板处理的0.5°和1.5°仰角雷达图像

再利用图5的函数模板和得分函数分别对图4（a）和图4（b）进行得分计算，得到图10两个得分值图像，从图10 中可以看出阵风锋图像处得分值明显较高，且具有较好的连续性，与图4 中阵风锋图像所在位置和形状有较高的一致性。

图10 阵风锋差分反射率因子图像函数模板计算的两个得分值图像

根据图2 的算法流程，利用图6 函数模板和得分函数处理前后时刻差分反射率图像（图4（b））得到风暴边缘图像得分值（图略），用图10（b）得分值结果减去风暴边缘图像得分值，降低风暴边缘图像对识别结果的影响。再分别对图9（c）、图10（a）和图10（b）得分值减去放射状回波图像得分值结果后，采用加权平均的方式融合3种得分值得到综合得分值结果，如图11（a）所示。由该图可见，雷达站南偏东约20°方向的放射状回波图像已得到较好的抑制。利用图8 平滑图像函数模板和得分函数处理图11（a），得到图11（b）最终的得分值结果。由图11（b）可见，本文设计的检测技术成功突出了阵风锋得分值，而抑制其他得分值。在实际中可根据算法在本地运行情况调整权重，进一步提高识别成功率。最后，利用区域增长法处理图11（b），得到结果如图11（c），对比实况图1 可以看出，算法能够较好地识别出阵风锋，且识别结果与原阵风锋位置和形状比较一致。

图11 多次处理后得分值及自动识别结果

3 结束语

本文根据我国易引发大风的阵风锋雷达回波图像高度较低、呈弧线结构和具有移动特点的时空变化特征，在3 种不同反射率因子图像上，采用模板匹配技术分别对其进行处理，得到3种不同得分值结果。对3 种不同得分值结果进行加权平均后，采用能够使得分值平滑细化的函数模板技术进行处理，最后采用两次不同的区域增长法提取阵风锋，并排除虚假预警。在得到3种得分值结果的过程中，还采用了风暴边缘图像抑制模板技术和放射状回波图像模板抑制技术。结果表明，算法能够有效识别阵风锋，识别结果与原阵风锋位置和形状比较一致。

作为会对生活和生产带来严重影响的天气因素，阵风锋常年带来损失不容忽视。本文所设计的阵风锋模板匹配技术为能准确地监测和预警阵风锋导致的灾害性天气提供多一种可能。目前本文中涉及的算法函数模板以及得分函数中的阈值集均来自江淮地区多年阵风锋图像特征统计结果，适用于该地区的阵风锋特征识别，若要在其他地区进行推广，还需要结合其他地区阵风锋图像特征进行相应阈值的调试；另外本文设计的算法仅利用了雷达反射率因子图像，阵风锋在速度图像上还具有辐合特征，但是雷达仅能探测到径向速度，且算法对径向速度图像质量有较高要求，导致在目前阵风锋图像特征识别中雷达径向速度图像数据难以利用。后期作者将尝试结合天气雷达径向速度图像特征以及地面自动站气象要素变化特征，以进一步提升基于阵风锋的大风识别算法。