基于温度层强回波区域面积的天气雷达雷暴自动识别方法

解妍琼 张云 杨波 王佳 王亚 张鹏

(1 国防科技大学 气象海洋学院，南京 211101；2 中国气象局 云雾物理环境重点实验室；北京 100081；3 陆军工程大学 教研保障中心，南京 210007)

引 言

天气雷达雷暴自动识别与跟踪是自动天气分析和雷电临近预报中的一个基本问题。雷暴一经识别，其空间位置、移动路径等发展特征可进一步确定。利用天气雷达回波数据进行雷暴(或风暴)识别、跟踪和临近预报的研究已经进行了几十年[1-5]。天气雷达反射率因子是雷暴识别中最常用的参数，识别风暴的最简单方法是设定一个合理的雷达反射率因子对图像进行阈值处理，由于该方法简单实用，在许多风暴识别方法中得到了应用。两种著名的方法是TITAN(雷暴识别、跟踪、分析和临近预报)[2]和SCIT(风暴单体识别和跟踪)[3]。还有很多关于雷暴雷达回波特征的研究[6-10]，主要以雷达的反射率因子为研究对象。如Trapp， et al[4]研究了1998—2000年发生的3 828次龙卷的雷达反射率因子图像，以确定准线性对流系统引发的龙卷风的百分比。

在对天气雷达回波雷暴特性的研究中，许多研究给出了与不同温度层和反射率因子相关的统计结果[11-16]。Martinez， et al[17]发现，当强度超过40 dBZ的回波对应的高度大于7 km时才会发生闪电。LIU， et al[18]利用TRMM(Tropical Rainfall Mea￣suring Mission)1998—2010年共13 a的观测资料研究了雷达反射率因子垂直剖面与闪电频度(Lightning Flash Rates, LFR)之间的关系，检验了20、30和40 dBZ回波顶温度与闪电活动的相关性。结果表明，回波顶温度与闪电活动的相关性并不高，而混合相态的强回波区大小与闪电发生率及闪电发生区域有更高的相关性。文献[19—21]都将-10 ℃高度层对应的反射率因子值达到40 dBZ作为地闪活动的预报因子。

YANG， et al[22]提出了一种基于双偏振雷达资料识别雷暴的算法，利用双偏振雷达资料进行粒子类型识别，然后综合识别出的霰粒子面积与雷达反射率因子进行雷暴识别。结果表明，在CAPPI 0 ℃高度层强回波区域反射率因子识别阈值为30 dBZ，霰粒子区筛选阈值为2 km2时，得到了最优的雷暴识别结果，识别准确度达到91%以上，识别虚警率(False Alarm Rate, FAR)仅为6.9%，临界成功指数为85.3%，这一方法能够达到很好的雷暴识别效果，但必须基于双偏振雷达。虽然我国业务雷达正在进行双偏振的升级改造，但相当长一段时间内，单偏振雷达仍是业务主体。如何使单偏振天气雷达具备相应的雷暴自动识别能力一直是气象业务的热点和难点问题，本文在前期研究基础上建立一种基于温度层强回波区域面积的天气雷达雷暴识别方法，解决单偏振天气雷达雷暴自动识别问题，并与文献[22]提出的算法进行比较，分析本文算法的识别效果，为提高单偏振天气雷达雷暴自动分析能力提供技术手段。

1 数据来源

本章利用S波段双偏振多普勒天气雷达2014年和2015年雷暴探测实验观测资料结合江苏省闪电定位资料、地面大气电场资料和国家气候基准站南京站点的探空数据研究雷暴的雷达回波特征。观测使用的S波段双偏振多普勒天气雷达，性能指标可参考文献[22-23]，不再赘述。江苏省闪电定位系统为ADTD型雷电定位系统，该系统综合使用磁定向法以及多站时差定位法探测地闪，定位精度约为500 m，探测效率约为90%～95%。

表1 样本资料简要说明Table 1 The brief description of sample data

为研究雷暴过程的双偏振雷达回波特征，选取2014—2015年共17个雷暴天气过程进行分析，具体样本日期见表1。为保证数据的有效性，选取大气电场仪或闪电定位系统监测区域内的雷达回波进行分析，采用了地面大气电场资料、江苏省闪电定位资料对对流单体进行筛选，区分雷暴单体与非雷暴单体。

2 雷暴识别算法

由于雷暴是一种中小尺度强对流天气系统，在水平范围内具有较强的雷达回波反射率因子，对流高度能够达到几公里甚至十几公里，且对流越旺盛越有利于雷暴起电，因此本文雷暴识别算法的基本思想是利用雷达回波体扫数据首先在水平范围内搜索强回波区域；再根据某一温度层上，达到设定阈值的雷达回波反射率因子积分面积，判断水平范围强回波区域是否符合雷暴特征；当强回波区域符合雷暴特征时，判定该强回波区域为雷暴区域。根据雷暴的雷达回波特征，采用不同温度层的雷达回波反射率因子积分面积实现雷暴识别。

2.1 算法总体流程

识别算法总体流程如图1所示，主要包括：

(1) 将雷达体扫数据作为输入量，完成数据解析；

(2) 利用体扫回波数据各层的反射率因子生成组合反射率数据；

(3) 通过设定的反射率因子阈值Zth1，利用组合反射率数据进行雷达回波强中心搜索；

(4) 通过设定的区域尺度阈值A1，对组合反射率数据雷达回波强中心进行区域尺度筛选，筛选出区域尺度大于区域尺度阈值A1的回波区域，筛选结果构成待定雷暴区域集合C1；

(5) 利用步骤(1)解析得到的体扫数据，生成各温度层CAPPI反射率因子数据；

(6) 利用某温度层CAPPI反射率因子阈值Zth2进行雷达回波强中心搜索；

(7) 通过设定的区域尺度阈值A2，对CAPPI反射率因子雷达回波强中心进行区域尺度筛选，筛选出区域尺度大于区域尺度阈值A2的回波区域，构成待定雷暴区域集合C2；

(8) 将第4步得到的待定雷暴区域集合C1与第7步得到的待定雷暴区域集合C2进行匹配，当集合C2中某区域中心的水平投影落在集合C1中某区域的范围内时，判断为匹配成功，将C1对应的待定雷暴区域确定为雷暴单体。

其中，步骤(5)中的CAPPI的高度根据探空数据的温度确定，本文主要采用了0 ℃、-10 ℃和-15 ℃对应的高度，由于探空数据一般情况不会出现0 ℃、-10 ℃和-15 ℃这几个高度，因此，温度层对应的高度在实际处理过程中采用插值方法得到。

图1 雷暴识别算法流程Fig.1 The overall process of thunderstorm identification algorithm

2.2 雷暴强中心提取

(1)强回波点提取

图2 强回波点提取效果：(a)原始雷达回波；(b)强回波点提取结果Fig.2 Effect of intense echo point extraction: (a) the original radar echo image;(b) the intense echo point extraction results

设组合反射率雷达回波图和CAPPI反射率因子回波图为H×V像素的图像，则图像共包含H×V个像素点，由ZH(H，V)表示各像素点的雷达回波反射率因子值，当ZH(i，j)>ZTh时，ZH(i，j) 即为强回波点，其中，i为回波图横坐标，j为回波图纵坐标，0

(2)强回波段合成

在横向或纵向上，逐点对上一步骤得到的强回波点集合ZHC进行相邻点合并，构成横向或纵向强回波段集合，表示为ZHS(N)，N为强回波段集合大小。以纵向强回波段为例，强回波段ZHS(x)是由多个横坐标相同，在纵向坐标位置连续的点构成的一维回波数据点集合。当某强回波点ZHC(i，j)为孤立点时，则判定该回波点为无效点，直接删除。

(3)强回波区域合成

对强回波段合成后得到的强回波段集合ZHS(N)进行相邻强回波段合成，即可构成强回波区域集合ZHA(M)，M为强回波区域集合大小。当某强回波段ZHS(x)为孤立强回波段时，则判定该回波段为无效回波段，直接删除。

强回波区域合成后，可利用边界点识别算法提取各强回波区域的边界点，以便于CAPPI识别出的强回波区域和组合反射率识别出的强回波区域进行水平投影位置匹配。本文采用遍历各强回波点的算法提取边界点，如图3所示。红色区域表示强回波区域，每一个格点代表一个回波点，遍历各格点数据，当某数据点的相临点个数小于4时，该数据点为边界点，图中A、B两点的相邻点的个数为4，判定为非边界点，除这两点外的各数据点的相邻点个数都小于4，判定为边界点。采用该算法即可提取出强回波区域的边界。

图3 强回波区域边界点提取算法示意Fig.3 Schematic diagram of boundary point extractionalgorithm in intense echo region

图4 强回波区域边界提取效果Fig.4 Effect of intense echo region boundary extraction

根据上述算法，采用合适的参数即可提取出强回波区域，提取效果如图4所示。强回波区域包含强回波段，强回波段包含强回波点，因此强回波区域包含区域内各强回波点的坐标信息和反射率因子信息。各强回波区域的边界可以从回波反射率因子图中提取出来。实际处理算法中可对提取出来的每个强回波区域进行编号，再通过遍历的方式对这些强回波区域进行区域尺度筛选和位置匹配。

2.3 区域尺度筛选

由于雷暴单体的尺度一般在几公里到上百公里，对应尺度过小的区域，可通过区域尺度筛选将其判断为非雷暴区域。雷暴识别算法总体流程中，组合反射率回波图和CAPPI回波图中都采用了区域尺度筛选，对尺度过小的区域进行筛选。在组合反射率回波图中以区域尺度阈值A1，CAPPI回波图中以区域尺度阈值A2进行筛选，当强回波区域的面积小于对应的区域尺度阈值时，强回波区域集合中删除该区域。将组合反射率筛选得到的待定雷暴区域集合C1，CAPPI在某高度层上筛选得到的待定雷暴区域集合C2。

2.4 筛选结果匹配

最后将待定雷暴区域集合C1和C2进行水平投影位置匹配，若匹配成功，即判定对应的组合反射率区域为雷暴区域。水平投影匹配基本思路是首先利用2.2节中提取的各强回波区域边界，求取CAPPI各强回波区域的几何中心点Pt，当Pt的水平投影落在组合反射率的某个待定雷暴区域时，则判定匹配成功。

3 结果检验与分析

为检验上述雷暴识别算法的有效性，利用表1所示天气过程的雷达回波数据和对应的闪电定位数据进行验证，为便于雷暴识别准确率统计，将闪电定位数据通过经纬度映射到雷达回波图像中。

根据上述雷暴识别算法流程，需要在组合反射率和CAPPI中分布搜索反射率因子强回波中心。由于本算法雷暴识别的关键在于CAPPI不同高度层上反射率因子阈值参数选择，为减小统计工作量，将组合反射率强回波中心搜索阈值设为固定值40 dBZ，区域尺度阈值A1设为固定值1 km2。在此基础上，识别过程中还需设置CAPPI图像的强回波中心搜索反射率因子阈值Z2，温度层HT以及反射率因子面积阈值A2，构成识别算法参数组合(Z2，HT，A2)。

在温度层HT的处理上本文算法做了近似，由于雷暴是强对流系统，其内部存在强烈的上升气流和下沉气流，云内某一温度层不可能在固定的高度，考虑到环境温度是雷暴起电的重要影响因素，在大量相关的研究中都采用探空数据给出的温度层(固定高度层)作为宏观特征量对雷暴进行研究，因此本文算法也将温度层近似认为是一个固定的高度，利用最近的探空站点最接近的时间点探空数据判定温度层的高度，然后在CAPPI中得到不同高度的反射率因子图像。

为检验各种参数的识别效果，采用探测概率 (Probability of Detection，POD)、漏警率 (Miss Rate，MR)、虚警率 (False Alarm Rate，FAR)和临界成功指数 (Critical Success Index，CSI)指数来衡量检测结果，分别表示为：

(1)

(2)

(3)

(4)

其中：R1，R2和R3分别为识别准确数、漏警数和虚警数。

检验过程中，利用上述识别雷暴算法结合识别参数组合，生成识别结果回波图像，再将雷达回波数据对应时段的闪电定位结果叠加到识别结果回波图像上，如图5所示。其中椭圆内的强回波区为通过上述识别算法自动计算得到的雷暴区域；当闪电定位结果落在雷暴识别的区域时，判定雷暴识别结果准确，相应R1增加1；当某区域有闪电发生，但未识别为雷暴区域时，判定为漏警，相应R2增加1；当某区域被识别为雷暴区域，但该区域没有闪电发生时，判定为虚警，相应R3增加1。

现有的研究表明，将不同的温度层的反射率因子用于雷暴预测时，会得到不同的预测效果，且多项研究表明在-10 ℃高度层回波强度达到40 dBZ以上，能够较好地区分雷暴和非雷暴单体[17-19]。文献[21]研究了0、-10和-15 ℃高度层霰粒子分布对雷暴识别结果的影响，为对比有无双偏振信息对结果的影响，本文对应地也采用这三个高度层进行研究。

图5 雷达回波雷暴识别结果检验(“-”代表负闪电，“+”代表正闪电)：(a) 系统性雷暴识别结果；(b)孤立单体雷暴识别结果；(c)飑线雷暴识别结果；(d)混合性雷暴识别结果Fig.5 Radar echo thunderstorm identification result test chart( “-” represents negative lightning and “+” represents positive lightning): (a) identification results of systematic thunderstorm; (b) identification results of isolated thunderstorm cell;(c) identification results of squall line thunderstorm;(d) identification results of mixed thunderstorm

检验过程中，首先考虑Z2为30 dBZ时的参数组合，研究HT为0、-10和-15 ℃层高度上区域尺度阈值A2分别是1、2和3 km2时的情况，此时，共有9种参数组合，利用这9种参数组合对表1所示天气过程的雷达回波数据进行雷暴识别，识别统计结果如表2所示，统计的有闪电数据对应的雷暴单体共有312个。其中，识别结果较好的参数组合是HT为0 ℃高度的3种参数组合，第1种参数组合(HT为0 ℃高度，A2=1 km2)能够得到最高的探测概率，POD达到89.7%，但虚警率略高，FAR为41.1%，临界成功指数为55.2%；从CSI看，第5种参数组合(HT为-10 ℃，A2=2 km2)能够得到最佳的识别效果，POD虽然低于第一种参数组合，仅达到了71.8%，但虚警率相对较低，FAR为20.6%，CSI为60.5%。此外，HT为-10 ℃时的三种组合CSI都优于0 ℃和-15 ℃的几种组合。

再研究Z2为35 dBZ时的参数组合，与Z2为30 dBZ时的分析方法相同，HT设为0、-10和-15 ℃，A2设为1、2和3 km2，共9种参数组合，识别统计结果如表3所示，可发现在0 ℃高度层上的识别效果最好，(35 dBZ， 0 ℃， 2 km2)参数组合能够获得最好的识别效果，POD达到87.5%，虚警率为32.9%，临界成功指数为61.2%。0 ℃高度层上的3种参数组合POD和CSI相对于其余两种高度层的参数组合，有较大的优势，但FAR则相对较高。

最后考虑ZTh为40 dBZ时的参数组合，采用相同的(HT，A2)参数组合进行检验，识别统计结果如表4所示，可以发现对应的9种参数组合都仅能得到很低的POD，最高仅达到了43.2%，CSI最高仅为42.2%。

表2 Z2=30 dBZ时的雷暴识别结果Table 2 Thunderstorm identification results (Z2=30 dBZ)

表3 Z2=35 dBZ时的雷暴识别结果Table 3 Thunderstorm identification results (Z2=35 dBZ)

由以上统计结果，在某温度层上，对反射率因子大于阈值的区域进行面积求和，将总面积作为雷暴识别的基本指标合理且能够得到较好的识别结果。又由于不同雷暴的对流发展不同，导致不同雷暴在不同的温度层的超过反射率因子阈值的面积不同，对流越旺盛，在较高高度上的反射率因子越大，越符合雷暴的特征，雷暴的闪电活动也越剧烈。在将-10 ℃或-15 ℃高度层的反射率因子阈值作为识别参数时，尽管识别POD相对较低，但识别到的雷暴区域是所有雷暴区域中闪电活动相对较为剧烈的区域，因此，以强雷暴筛选为目的时，可以选择相应的参数组合。从CSI结果来看，27种参数组合中，(35 dBZ， 0 ℃， 2 km2)参数组合能够达到最优的识别效果。

4 结论

本文提出了基于温度层强回波区域面积的天气雷达雷暴自动识别方法，该方法将温度层、反射率因子阈值和强回波区域面积相结合，包括组合反射率强反射率因子区域筛选、尺度筛选，某温度层上CAPPI强反射率因子区域筛选，组合反射率识别区域和CAPPI识别区域匹配几个步骤，该方法中的关键在于选择恰当的CAPPI筛选温度层、强回波区反射率因子阈值和强回波区域面积尺度阈值。研究中，采用了17个雷暴天气过程的312个雷暴区域检验识别算法，结果表明：当 CAPPI数据处理中采用(35 dBZ， 0 ℃， 2 km2)参数组合能够达到最优的识别效果，识别概率可达87.5%，虚警率为32.9%，临界成功指数CSI为61.2%。该方法可以用于天气雷达的雷暴计算机自动识别业务。此外，虽然对比文献[22]将双偏振参数用于雷暴识别，本文的雷暴识别能力要弱一些，但本方法实现了只有单偏振雷达时的雷暴自动识别，为目前气象业务提供重要指导。