1 引言

飑线是天气雷达进行探测时，在雷达回波图像上呈现出的明显的带状或类弧状图像，该形状区域的回波强度明显高于该区域以外的周边区域。通常飑线是由多个雷暴单体或超级单体近似一条直线或弧线排列组成的云雨带，其宽度约为几十km，长度约为几十至数百km，其初步形成至基本消散的持续时间可达数小时至十几小时^[1-2]。

以往的研究表明，飑线的形成与锋面系统、低层辐合线、出流边界以及地形或热力抬升等作用息息相关，同时，大气内部的扰动和不稳定也是诱因之一^[3]。当飑线出现时，往往会出现气压陡升、温度骤降、风向突变和风力加强等强对流性天气现象，同时还可能伴有雷暴、闪电、冰雹或龙卷等灾害性天气。虽然飑线不是造成灾害性天气的直接原因，但飑线的出现往往预示着破坏性大风、龙卷、冰雹和短时强降水的发生，且一旦发生，其破坏力是非常巨大的。例如，2009年6月3日20时30分左右(北京时间，下同)，河南商丘地区突遭强对流天气袭击，风力最大达11级，并伴有暴雨、冰雹，整个过程持续约50分钟。据灾后统计，全市172个乡镇中有137个不同程度受灾，倒塌房屋4 800多间，损毁树木769.3万株，倒伏小麦309.84万亩，20人失去生命，直接经济损失高达15亿元^[4]。灾害性天气每年都在发生，且随着人们生产活动的不断扩展，因其造成的影响也越来越严重^[5-7]。

对于有经验的气象专业人员而言，通过对雷达回波图像等气象资料的主观分析，能够准确地判读出飑线的形成和发展过程。而对于计算机程序而言，如何准确地识别出雷达回波图像中的飑线特征，尽可能避免漏判和误判，却是一项复杂的算法问题。Kelly等^[8]在其专利中提出利用二维傅立叶变换、基于反射率值的边缘检测和形态分析，对雷达回波图像进行处理，实现对灾害性天气的识别和图像展示。程凌舟等^[9]利用小波变换模极大边缘检测方法和Hu矩理论方法提取飑线中心强度较强且连接成线这一回波特征，实现飑线雷达回波的自动识别，从个例检验的效果来看该方法具有良好的效果。李哲等^[10]利用计算机分析反射率因子、长轴长度、系统面积和线性程度，并择优选取线性程度高且包含较强回波的飑线作为飑线最优识别结果，结果表明，其方法能够同时识别强、弱飑线，识别命中率比单一参数法分别提高了26.3%、10%、40%和89%。

由于各类高危高影响天气在雷达等探测资料中所呈现出的复杂性和多样性，使得对其识别追踪和预报的自动化变得极为困难。事实上，近5~ 10年来气象业务人员能够实时获得的气象观测和探测数据爆炸式增长，这些实时推送、更新的海量资料，其信息量已超出人们能够及时进行主观分析、判读的能力。因此，亟需借助计算机对包括飑线在内的各类短时强对流性天气进行自动化、智能化的识别，并能够满足识别的高准确率要求。

2 飑线在雷达上的回波特征分析

经过数十年的发展，国内外很多专家学者通过雷达探测资料、数值模拟和地面实况等资料，分析总结出一些有关飑线生成发展过程中的重要特征。

根据飑线的形态特征，美国国家强风暴实验室(NSSL)将飑线系统划分为断线型、后部建立型、碎块型和嵌入型^[11]。Rotunno等^[12]最早提出较系统、完备的飑线发展维持机理，并以水平涡度的方式刻画了飑线维持和传播的物理机制，其成果对后续研究奠定了重要基础。田荟君等^[13]对发生在合肥附近的一次暖区飑线过程进行了分析，雷达分析显示，飑线呈弓形回波，具有明显的对流区、层云区和过渡带，线尾涡旋位于弓形回波北端。万夫敬等^[14]利用雷达和地面观测资料对发生在山东东部地区的一次强雷暴大风天气过程进行了个例分析，提出此次飑线的形成与中尺度边界辐合线持续触发新的对流单体密切相关，同时，近地层的径向辐散和中层辐合对于大风的识别诊断也有重要意义。侯淑梅等^[15]利用自动站和雷达资料，对一次山东境内飑线发展过程进行了深入分析，并总结了一系列有关飑线发展过程中对流单体之间和弓状回波与对流单体之间的发展变化特征。张萍等^[16]分析了一次发生在江西境内的强飑线天气过程，提出弓状回波带由若干个对流回波核组成，回波长度超过100 km。吴瑞姣等^[17]通过对以往35条飑线的分析，发现江淮地区的飑线多为向东南方向移动，移速集中在8~16 m/s，飑线长度为200~250 km，回波中心强度为60~65 dBZ，持续时间为3~4 h。

通过天气雷达能够清楚地探测到飑线形成、发展的过程以及飑线大小和空间分布情况。在其形成阶段，雷达回波图像上一般表现为若干个相对独立的对流单体，各个单体的轮廓清晰，单体中心区域的回波强度超过40 dBZ，单体的水平尺度为15~20 km，各个单体空间分布尚较分散^[18]。这些对流单体在不断发展的过程中，体积逐步增加，单体顶高逐步抬升到10 km以上，并逐渐聚拢、合并形成带状的空间分布。在其旺盛阶段，雷达回波图像上出现明显的长条状或长弧状高值回波区域，回波强度达到50 dBZ甚至60 dBZ以上，回波顶高上升到15 km以上。高值回波区域的外沿回波强度梯度较大，表现为强超级单体或强对流单体群。此时，在雷达回波图像上即可表现出典型的飑线特征^[19-20]。

3 雷达回波飑线特征识别方法 3.1 拟解决的技术问题

随着对飑线形成机理和物理过程研究的深入，人们对飑线所表征的天气现象做出主观分析的经验也越来越丰富。然而，对于计算机自动识别飑线而言，却需要解决以下几点技术问题。

(1) 飑线所对应的高值回波区域在形态上往往是非连续的，仅凭选定某一阀值，并过滤掉低于该阀值的区域，得到的往往是若干个独立的高值小区域，计算机难以直接通过这些小区域判定飑线存在与否。

(2) 飑线所对应的回波强度难以给出客观定量的指标，不同的天气个例其飑线区域的回波强度和空间分布形态差异很大，难以采用传统图像模式识别技术做出判定。

(3) 由于飑线与地面自动观测站等气象设备记录到的天气现象(如大风、雷暴、短时强降水)之间缺少清晰的数学关系，所以，若直接将雷达基数据或回波图像与气象观测设备记录到的天气现象作为神经网络的输入，通过大规模历史个例数据集训练来寻求模型的最优参数集，其训练结果不论是采用传统神经网络模型或SVM，还是采用深度神经网络模型，往往数据集越大，模型的超参数集越难收敛，识别或预测的结果往往偏向于训练数据集中的最大概率的天气现象。这种情况下，即便检验时的统计误差很小，但空报率或漏报率至少有一项会非常高(往往超过90%)。

3.2 总体算法流程

本方法的主要目标是通过计算机自动识别出对流发展到旺盛阶段时，在雷达回波图像上显现出的长条状或长弧状高值回波区域。为了解决上述技术问题，本方法通过数值预处理、高通滤波、二值化降噪、图像特征提取、目标物的中轴线提取，以及飑线形态分析等一系列步骤，构建对雷达飑线特征的自动识别，计算流程如图 1所示。

步骤1：首先，定义变量α∈[0，max_α]，α表示雷达体扫的仰角编号，max_α表示该雷达体扫的最大仰角编号，max_α由雷达的硬件型号和体扫参数确定。α初始值为0。再定义变量Filt_A，用于记录高通滤波的阈值。Filt_A的初始值通常取[60,100]之间的正整数。

步骤2：读取待分析的雷达基数据文件，从该文件中提取α仰角面上的基本反射率数据，记为R(α)。为了便于计算，将R(α)由极坐标形式转换为平面直角坐标形式的数据，记为RP(α，x，y)。

步骤3：将RP(α，x，y)中的基本反射率数值进行高通滤波，滤波阈值为步骤1中定义的变量Filt_A。经高通滤波后的基本反射率记为RPF(α，x，y)，其中，x和y分别表示雷达探测点的横坐标和纵坐标，每个雷达探测点由(x，y)确定，以下简称探测点或像素。

步骤4：对RPF(α，x，y)进行二值化处理，将步骤4中所有被过滤的像素值设为0，所有未被过滤的像素值设为1，由此得到二值化后的RPF(α，x，y)，记为RPFB(α，x，y)。

步骤5：借鉴高斯滤波^[21]的基本思想，对RPFB(α，x，y)进行“降噪”处理，同时剔除孤立杂波点和被高通滤波过滤掉的孤立点，计算结果记为RPFBD(α，x，y)。具体计算方法详见3.3节。

步骤6：从RPFBD(α，x，y)中识别目标物的中轴线。为了更清楚地表述，可将二维数组形式的RPFBD(α，x，y)视为一张二值化图像，RPFBD(α，x，y)中像素值为1的点对应该图像中的雷达基本反射率高值区，即待识别的目标物，像素值为0的点视作背景图像。目标物在图像上表现为若干个大小不均的色块，该步骤的目的是从这些色块中提取出色块的中轴线，记为RPFBD_kn(α，x，y)。具体计算过程详见3.4节。

步骤7：由于飑线在雷达回波图像上呈现为带状或类弧状图形，因此，可用一条直线来近似带状或类弧状图形的中轴线。换言之，如果存在飑线，就应当能从RPFBD_kn(α，x，y)中找到一条或多条直线，使得RPFBD_kn(α，x，y)中大量凌乱无章、长短不一的线段总体上围绕这些直线分布。该步骤是为了判别RPFBD_kn(α，x，y)中是否存在带状或类弧状图形的中轴线。具体计算过程详见3.5节。

步骤8：上述步骤如果识别出飑线，则将一组数据[α，Filt_A，FN]记录到数据集Rlt中，进入步骤9。其中，FN是一个计数器，其初值为0。如果没有识别出飑线，则先令FN = FN + 1，再令Filt_A = Filt_A - FN×Filt_B，其中，Filt_B是一个定值参数，一般取1~10的正整数。此时，如果Filt_A ≥ Filt_N，则将α重置为0，返回步骤2；反之，则进入步骤9。其中，Filt_N是一个定值参数，0＜Filt_N＜Filt_A，Filt_N的物理意义是，小于Filt_N的基本反射率往往不可能发生飑线。

步骤9：令α=α+1，Filt_A重置为初始值，返回步骤2，继续进行步骤2~步骤8的计算。直到α > max_α，进入步骤10。

步骤10：分析数据集Rlt，该数据集有可能为空，也可能由1组或多组[α，Filt_A，FN]记录构成。Rlt的实践意义在于，当其为空，表示所分析的雷达数据没有识别出飑线；当其不为空，表示识别出了飑线，且[α，Filt_A，FN]的组数越多，表示在多个雷达探测仰角面均识别出了飑线，飑线特征很强烈；同时，FN的值从0开始，数值越小，表明该仰角面上飑线特征越明显，反之，飑线特征相对较弱。

3.3 二值化回波图像降噪算法

由于雷达回波图像经过高通滤波、二值化处理等预处理计算后，图像中产生大量孤立杂波点和被高通滤波过滤掉的孤立点，这些点对于后续的强回波色块特征提取造成不可忽略的干扰。为了改善这一问题，需要对二值化后的回波图像进行降噪处理。算法的核心思想与高斯滤波^[22-23]相近，但方法和效果上又有显著不同。

传统高斯滤波过程中常用的高斯函数如下式所示：

$ G\left( {x,y} \right) = \frac{1}{{2\pi {\sigma ^2}}}{{\mathop{\rm e}\nolimits} ^{ - \frac{{{x^2} + {y^2}}}{{2{\sigma ^2}}}}} $

(1)

式中，σ表示高斯核，其大小决定了高斯函数的宽度。通过高斯核卷积或进行傅立叶变换，能够有效降低孤立点噪声，而对于被高通滤波过滤掉的略小于Filt_A的孤立点，传统高斯滤波却无能为力。本算法的改进和计算过程如下。

首先，从RPFB(α，0，0)开始，划定一个边长为SL的正方形，分别统计该正方形区域内像素值为0和为1的像素数量，分别记为Stat_0和Stat_1。其中，SL为正整数且SL∈[3,50]。如果Stat_0 > k× Stat_1，则将该正方形区域内所有像素的值设为0；反之，则将该正方形区域内所有像素的值设为1，其中，k为经验阈值，k∈[0.1，10]。再以SL为步长，向右或向下滑动上述正方形区域，重复上述计算过程。直到所有像素完成一次遍历计算。

如果RPFB(α，x，y)中x或y的最大值不是SL的整数倍，会存在这样一种特殊情况，当边长为SL的正方形以SL为滑动步长，滑动到RPFB(α，x，y)最右侧或最下端，则会发生边长“溢出”，此时将正方形“溢出”的像素默认取值为0，再按照上述方法进行处理。

3.4 目标物中轴线的提取算法

为了从上述步骤6所述的若干个大小不均的色块中提取出色块的中轴线，本算法采用图像外周层层剥离的思想，具体计算过程如下。

(1) 从一个固定方向，依次分析RPFBD(α，x，y)中每个像素点的值。如果值为0，则跳过，继续分析下一个像素；如果值为1，进入步骤(2)。为了便于表述，不妨定义当前像素值为1的点的x = i，y = j。

(2) 以坐标(i，j)为中心，统计与之紧相邻的8个像素中，像素值为1的像素数量C1，如果C1=3，即RPFBD(α，i，j)的周边存在3个雷达基本反射率为高值的像素，则将当前位置的值置为0，即 RPFBD(α，i，j)= 0。根据步骤(1)、(2)，完成所有像素点的遍历，输出结果记为RPFBD_k1(α，x，y)。

(3) 依次分析RPFBD_k1(α，x，y)中每个像素点的值，如果值为0，则跳过；如果值为1，则统计与之紧相邻的8个像素中，像素值为1的像素数量C2，如果C2∈[3,4]，即当前像素的周边存在3个或4个雷达基本反射率为高值的像素，则将当前位置的值置为0。本步完成所有像素的遍历，输出结果记为RPFBD_k2(α，x，y)。

(4) 与步骤(3)类似，分析RPFBD_k2(α，x，y)中每个值为1像素，统计与之紧相邻的8个像素中，像素值为1的像素数量C3，如果C3∈[3,4,5]，则将当前位置的值置为0。完成所有像素的遍历，输出结果记为RPFBD_k3(α，x，y)。

(5) 与步骤(3)类似，分析RPFBD_k3(α，x，y)中每个值为1像素，统计与之紧相邻的8个像素中，像素值为1的像素数量C4，如果C4∈[3,4,5,6]，则将当前位置的值置为0。完成所有像素的遍历，输出结果记为RPFBD_k4(α，x，y)。

(6) 与步骤(3)类似，分析RPFBD_k4(α，x，y)中每个值为1像素，统计与之紧相邻的8个像素中，像素值为1的像素数量C5，如果C5∈[3,4,5,6,7]，则将当前位置的值置为0。完成所有像素遍历，输出结果记为RPFBD_k5(α，x，y)。

(7) 经过上述计算，可将雷达基本反射率高值区所表征的各个色块向其中心位置收缩，从图像上直观地看，各个色块的位置没有变化，但面积在缩小。

(8) 为了得到各个色块的中轴线，需要进行多轮步骤(1)~(7)的迭代计算，并将上一轮的RPFBD_k5(α，x，y)作为下一轮中步骤(1)的RPFBD(α，x，y)，直到不再有新的像素值被置为0，停止迭代计算。输出结果记为RPFBD_kn(α，x，y)，图像呈现为大量凌乱无章、长短不一的线段。

3.5 基于霍夫变换的特征提取算法

霍夫变换(Hough Transform)是Paul Hough于1959年公开的一种图像处理方法，主要用于识别照片中的复杂图案，并于1962年获得美国专利。该专利对直线采用斜截距的参数化，即$y = ax + b $的表达式，由于斜率可能变成无穷大，进而导致无限变换空间，为了解决这一问题，Richard Duda等^[24]于1972年提出了广义霍夫变换，该思想一直沿用至今。

霍夫变换能够用于检测目标图像中的直线、线段和圆弧等图形。为了判别上述步骤6中的RPFBD_kn(α，x，y)中是否存在带状或类弧状图形的中轴线，这里采用广义霍夫变换的基本思想，具体计算过程如下。

首先，采用法线式定义一条直线：

$ r\left( {x', y', w} \right) = x' \cdot \cos w + y' \cdot \sin w $

(2)

其中，r表示从坐标原点到直线的最短距离，理论上r ∈ [-∞, ∞]，实际计算时r的最大值可取为RPFBD_kn(α，x，y)所对应的图像的对角线长度；ω表示直线与X轴正方向的夹角，ω ∈ [0 °, 180 °)，( x', y')表示平面直角坐标系中的点，对于平面直角坐标系中的任一条直线而言，其(r, ω)是固定不变的，不同的(r, ω)构成不同的直线。

然后，以ω_i为横向的单位步长，并以r_i为纵向的单位步长，构建一个二维矩阵MAT[p，q]，其中，p表示行，$p \in \left[ { - n \cdot r\_i, - \left( {n - 1} \right) \cdot r\_i, \cdots \cdots , - 2 \cdot r\_i, 0, r\_i, 2r\_i, \cdots \cdots , \left( {n - 1} \right) \cdot r\_i, n \cdot r\_i} \right] $，q表示列，$q \in \left[ {0, w\_i, 2w\_i, 3w\_i, \cdots \cdots , 180 - w\_i} \right] $。MAT[p，q]中所有单元的初值均为0。其中，ω_i用于定义直线与X轴正方向的夹角划分的粒度，ω_i值越小，横向(列)上的刻度数量越多，夹角被划分的粒度越细，ω_i通常取一个可整除180的正整数，如1、2、6、20或45等。同理，r_i用于定义直线到原点之间最短距离被划分的粒度，r_i值越小，纵向(行)上的刻度数量越多，最短距离被划分的粒度越细，n∙r_i和-n∙r_i分别表示MAT纵向(行)上的最大刻度和最小刻度。

接着，依次遍历RPFBD_kn(α，x，y)中各像素，如果值为0，则跳过不处理，接着遍历下一个像素，直到所有像素全部遍历一遍；如果该像素的值为1，则需进一步分析：假设当前像素的坐标为(i，j)，即x =i，y = j，按式(2)，将i代入 x'，j代入y'，再将MAT横向的每一个q，$q \in \left[ {0, w\_i, 2w\_i, 3w\_i, \cdots \cdots , 180 - w\_i} \right] $，逐个代入式(2)中的ω，由此得到一系列r(i，j，q)。对于每一个r(i，j，q)，从p中找到与之绝对值差值最小的任一项，记为p′，再将MAT[p′，q]单元的数值自增1。

最后，分析MAT[p，q]中各单元的数值大小，按式(3)筛选出的所有单元，这些单元的数量即为本方法识别出的飑线的中轴线的数量。

$ \frac{{{\rm{MAT}}\left[ {p,q} \right]}}{{{\rm{Max}}\left( {{\rm{MAT}}} \right)}} > Trd $

(3)

其中，Max(MAT)表示MAT[p，q]中所有单元数值的最大值，Trd∈(0，1)，其取值与雷达基本反射率图像中像素的数量、易发生强对流天气的季节、雷达探测的地理空间等气象学因子息息相关。Trd取值越小，识别飑线的漏报率越低，但误报率会增加；反之会增加漏报率，并降低误报率。可看出MAT[p，q]中各个单元所对应的(p，q)分别表征了一条直线，RPFBD_kn(α，x，y)中大量凌乱无章、长短不一的线段刚好围绕这些直线分布。

若式(3)计算得到的中轴线数量为0，则表明当前分析的雷达数据中不存在飑线；反之，则表明识别到飑线的存在。

4 个例检验与分析 4.1 数据说明

为了检测上述方法的有效性，这里首先选用2019年4月9日11:35:05(北京时，下同)南京地区的多普勒天气雷达基数据按上述计算流程进行分析。当日08—14时，安徽江南地区及江苏西部出现强对流性天气，自西向东相继出现中到大雨，局地出现暴雨并伴有雷雨大风的天气，最大小时雨强达到40~60 mm。11: 35: 05的雷达回波图像(0.5 °仰角)如图 2所示。

4.2 实施过程

按照3.2节所述计算流程，对上述雷达数据实施计算，各项参数在本实施例中取值如表 1所示。

首先，从最低一个仰角开始对雷达基数据进行预处理，经坐标转换得到的256阶灰度图像如图 3a所示。

使用Filt_A=60 dBZ进行高通滤波，得到RPF (0，x，y)，所绘制的图像如图 3b所示。经二值化处理得到RPFB(α，x，y)，所绘制的图像如图 3c所示。再进行降噪处理，同时剔除孤立杂波点和被高通滤波过滤掉的孤立点，得到RPFBD(0，x，y)，所绘制的图像如图 3d所示，可看出图 3d显现为若干个大小不均的白色块。接着，从这些白色块中提取出色块的中轴线，得到RPFBD_kn(0，x，y)，所绘制的图像如图 3e所示。为了清楚地显示图像细节，图 3f是对图 3e主体区域的放大显示，可看出图 3d所示的色块经上述步骤6的迭代计算，得到一些看似杂乱无章长短不一的线段，这些线段正是各个色块的中轴线。

接着，按照步骤7，定义一个法线式直线表达式，$r\left( {x', y', w} \right) = x' \cdot \cos w + y' \cdot \sin w $。本实施例中，上述RPFBD_kn (0，x，y)所对应的图像的边长为810，r_i=20，ω_i=5，$n \cdot {r_i} = \sqrt {{{810}^2} + {{810}^2}} \approx 1146 $，这里取r_i的整数倍，令n∙r_i =1 140，则n=57。依次遍历RPFBD_kn(α，x，y)中每一个像素，如果该像素的值为1，则将当前像素坐标(i，j)代入(x′，y′)，再将所构建的矩阵MAT横向的每一个q逐个代入法线表达式中的ω，由此得到一系列r(i，j，q)。举例来说，假设p ∈ [-10，-8，-6，……，6，8，10]，q∈ [ 0,45,90,135 ] ，如果r(i，j，0)= 7.5，则p′ =8；再如r(i，j，45)= -9，则p'取-10或-8均可，实际应用中，应当统一“就高”取-8，或统一“就低”取-10。然后，分别将MAT[8,0]和MAT[-8, 45](或MAT[-10, 45])的单元值加1。

进一步分析计算$\frac{{{\rm{MAT}}\left[ {p,q} \right]}}{{{\rm{Max}}\left( {{\rm{MAT}}} \right)}} > Trd$，此处Trd取值为0.75。当前，由RPFBD_kn(0，x，y)计算得到有且仅有MAT[70,240]单元的数值满足上述不等式，将ω=70，r=240代入$r\left( {x', y', w} \right) = x' \cdot \cos w + y' \cdot \sin w $，在图 3f的基础上绘制出的直线如图 3g所示。可看出该直线所表征的方向和位置与图 2中飑线的方向和位置是一致的。将由该步骤计算得到的一组数据[α，Filt_A，FN]记录到Alt中，即Alt = {[0,60,0]}。

接着，进入步骤9，再按照上述流程计算下一个仰角，直到α > max_α。

最后，分析数据集Rlt，由于本实施例中Rlt包含2组Alt，且FN = 0，因此，可判定当前雷达资料中存在飑线特征。进一步通过Rlt和MAT计算出的法线式直线参数，可对雷达回波飑线区域做出标识，本实施例将上述直线贯穿区域紧相邻20 km范围内基本反射率值大于55 dBZ的回波做出高亮标识，结果如图 3h所示。数据集Rlt所记录的信息，既能够直接向气象专业人员提供客观、量化的飑线评价指标，也可集成到相关业务系统中，以更丰富的展现形式，实现飑线识别的功能目标。

采用上述相同方法及参数值，对2019年4月9日10时28分—11时23分共12个雷达基数据进行计算分析，各数据0.5 °仰角基本反射率图像及识别出的飑线的参数化直线如图 4所示。

从图 4中可看出，各条识别出的飑线的参数化直线与实际飑线在位置上存在1~10像素的偏差，这是因为识别过程中，参数r_i取值为20，即MAT矩阵的步长为20个像素。如果适当减小r_i，则能够降低这种位置偏差，但运算量及运算时长将有所增加。其中，图 4a还出现了一条夹角约45 °的斜线，实际回波图像中在两条直线相交处偏西南方向的确存在超过60 dBZ的回波，呈细条状，但空间尺度较小，而非气象意义上的飑线，这里属于误判。图 4c中，实际飑线的形态和空间尺度与图 4b和图 4d相近，这里却发生漏判，可能是由于图 4c提取出的中轴线映射到MAT矩阵的规律性较差，使用$\frac{{{\rm{MAT}}\left[ {p,q} \right]}}{{{\rm{Max}}\left( {{\rm{MAT}}} \right)}} > Trd $条件不满足，从而无法识别出来。后经测试，将参数Trd调整为0.72，或将Filt_A降低到55 dBZ，又能识别出与图 4b和图 4d相近的一组[p，q]参数。图 4各幅图像仅显示了Filt_A = 60 dBZ，FN =0时所识别出的[p，q]参数，再由此参数绘制出的直线。按照上述计算流程，当FN增加，Filt_A随之减小时，被高通滤波后的色块面积有所增加，可识别的区域及中轴线数量也将相应增加，这可在一定程度上降低单个雷达数据，其飑线特征被漏判的机率，如图 4c出现的这种情况。同时，当雷达体扫的仰角面增加，最终形成的Alt数据集中[α，Filt_A，FN]记录组数也可能会增加，进而能够起到降低漏判率，提高识别成功率的目的。

4.3 检验结果分析

前文以2019年4月9日发生在江苏地区的一次强对流性天气过程为例，对本文所述雷达回波飑线识别的方法进行实例检验，对于单仰角面上形态特征显著的飑线，本方法能够较准确地做出识别，但也存在个别漏判和误判。通过组合不同仰角面上的回波，以及降低高通滤波的阈值等手段，可提高识别的成功率。

前文图 3g和图 4中的直线是根据上述计算步骤中二维矩阵MAT计算出的法线式直线参数绘制而成。由于判定飑线存在与否的关键性指标是经式(3)筛选出的MAT单元的数量，也就是本方法识别出的飑线的中轴线的数量，当该数量大于0，则表明识别到有飑线存在，当该数量为0，则表明未识别到飑线。因此，在实际应用中，如图 3g和图 4所示的直线并不会也不需要被绘制出来。这里绘制的目的是为了以图形化的方式来检验本方法实施过程中参数[p，q]的计算结果。可看出图 3g和图 4中绘制的大多数直线的方向和位置与实际飑线的方向和位置基本一致。

4.4 进一步的检验与分析

为进一步检验上述方法的普适性，我们收集了2017年3—6月南京和合肥两部雷达的全部基数据用以检验分析。检验借鉴了钟敏等^[25]提出的“雷达识别跟踪产品的评估分析”方法，采用击中率(POD)、虚警率(FAR)和临界成功指数(CSI)对识别的效果进行评估。

由于雷达设备存在停机维护等特殊情况，所收集到的基数据时间上有少量缺失，文件数分别为21 049和29 938。为了检验本方法的效果，需要先从这些文件中筛选出的确存在飑线特征的基数据，并以此作为检验基准。考虑到文件数量较大，我们对这些数据先进行了简单筛选，筛选的指标为：(1)基本反射率最大值超过55 dBZ；(2)大于50 dBZ的探测点数量大于50个；(3)雷达探测时间的前后各1 h内，该省内地面站观测风速有超过10 m/s的记录且同一时刻这样的站点数量超过3个。同时满足这3项指标的基数据文件数量分别减少到5 642和4 206。再通过人工判读这些基数据所绘制的最低2个仰角的回波图像，统计得到回波图像中存在显著飑线天气特征(一条连续的长条状或长弧状高值回波区)的文件数分别为191和87。将这些文件作为“存在飑线”的判定基准，并以此评估识别的POD指标。

实际上，对于飑线特征的筛选存在主观性，为了更客观反映识别的效果，我们又以地面站的观测记录来评估FAR指标。具体过程是，如果由本方法识别出存在飑线，则再分析相应地理空间范围内的地面站是否记录到大风或短时强降水等天气现象，如果有，则评价为识别成功，反之，则判为虚警。此外，考虑到地面站观测到的天气现象往往滞后于雷达回波形态所表征的现象，这里以每次雷达探测时间为基准，后续2 h内的观测记录作为判定参考。

仍采用上述方法及表 1所示的参数值，但其中的参数k和Trd 增加为多组，其中，k分别取值为[0.10，0.15，0.20，0.25，0.30，0.35，0.40]，Trd分别取值为[0.60，0.65，0.70，0.75，0.80，0.85，0.90]，对上述全部基数据进行识别。得到POD指标如图 5所示。

图中横坐标为不同取值的参数Trd，纵坐标为不同取值的参数k，图中不同颜色所表示的是POD指标。可看出总体上Trd取值与POD呈现负相关，k取值对POD的影响相较Trd要小，这些参数组合中，(k，Trd)=(0.25，0.60)表现出相对较优的性能，POD指标最高达到95%。

表 2列举了几组(k，Trd)不同取值组合时三项评价指标的得分情况。

与图 5的结论相同，表 2中Trd取值越小，POD得分相对越高，但CSI在Trd=0.75时表现相对最好。可看出FAR和CSI得分普遍较差，经分析，其中一个重要原因是检验数据中存在一些如图 6所示的形态“异常”的回波图像，本方法未能将其与飑线加以准确区分，产生了误判。另一方面，由于近地面层结稳定度状况的差异，即使出现飑线形态的雷达回波，地面也不一定会观测到大风，采用大风作为记录飑线可能漏掉许多线状对流，这也是导致FAR居高不下的原因之一。

总体而言，本方法对于雷达回波图像中飑线形态的识别，具体较高的识别率，POD达到95%，但由于算法对于异常回波图像分析的不完备，以及检验方法的缺陷，当前FAR指标过高，进而也影响了CSI的得分。

5 小结

本文提出一种雷达回波图像中飑线特征自动识别的方法，以多普勒天气雷达探测资料为主要数据源，对雷达探测到的基本反射率的空间分布和强度进行分析，利用数值预处理、高通滤波、二值化降噪等技术手段解决高值回波区域不连续的问题，采用图像特征提取、目标物的中轴线提取等技术手段解决飑线所对应的回波强度难以客观定量判定的问题。

该方法的应用能够将以往需要由气象专业人员主观分析、判读雷达回波图像的工作自动化、客观化，提高了飑线识别、强对流天气预警相关业务的准确性和时效性。此外，经本方法识别出飑线后，将判定结果作为神经网络的输入，将有助于提高基于机器学习技术的短临天气预报相关业务模型训练的可靠性。

在实践业务中，气象专业人员仅通过研读某一个仰角面的雷达回波图像，往往无法准确地判定是否存在飑线。本方法能够将不同仰角面的基本反射率数据加以综合分析，进一步提高了自动识别飑线特征的准确性和可靠性。后续研究计划将该方法作为前置基础，运用深度神经网络，将雷达回波图像和本方法的结果分别作为深度网络的输入，并以自动站等地面观测记录为天气现象发生的依据，结合光流场和外推算法，实现雷雨、大风、冰雹等灾害性天气的智能识别和预警。