雷达回波反射率垂直剖面图的冰雹识别方法

路志英，朱俊秀，田 硕，贾惠珍

(1. 天津大学电气与自动化工程学院，天津 300072；2. 天津市气象局，天津 300074)

雷达回波反射率垂直剖面图的冰雹识别方法

路志英1，朱俊秀1，田 硕1，贾惠珍2

(1. 天津大学电气与自动化工程学院，天津 300072；2. 天津市气象局，天津 300074)

为了更加准确地识别冰雹天气，对某地区2005—2010年间24个冰雹过程和19个暴雨过程的825个雷达样本数据进行了分析处理，建立了基于雷达回波反射率垂直剖面图的冰雹自动识别的客观模型．在分析冰雹云体形成机理和结构的基础上确定最佳剖线，并用插值法生成雷达回波反射率垂直剖面图，通过图像处理方法，提取特征数据(强回波(45,dBZ以上)与0,℃和－20,℃温度层的高度差、弱回波区和有界弱回波区的宽度和高度)．然后采用粗未选集理论数据挖掘方法对相关特征数据进行处理，建立了自动识别冰雹天气的客观模型．测试结果表明：该识别模型的判别规则对28个冰雹天气过程的383个样本的正确识别命中率是82.77%，可有效地识别和预报冰雹，有助于减轻冰雹灾害天气造成的损失．

图像处理；弱回波区；有界弱回波区；粗糙集；数据挖掘

冰雹作为一种强对流天气会对人们的生命财产造成严重的威胁．因此，准确地识别和预报冰雹天气成为气象研究人员的重点目标之一．随着雷达探测技术的发展，现已能够获得大量的实时数据，通过对这些数据的分析，能够对雷达回波反射率图和径向速度图进行相关的特征提取，如强度特征、形态特征及梯度特征的提取[1-2]．另外，Montanya等[3]和李磊等[4]还研究了闪电活动对冰雹云的影响．但事实表明仅仅通过这些特征来进行预报并不能达到十分准确的效果．因此一些气象人员通过对典型的冰雹天气的雷达回波剖面图进行分析，总结出了雷达反射率垂直剖面图(RCS)产品反映冰雹天气过程的主要特征，如弱回波区和有界弱回波区(穹窿)等特征[5-6]．赵俊荣等[7]分析了冰雹过程的中气旋、液态水含量和有界弱回波区的特征，并介绍了这些特征与冰雹落点的关系，为人工消雹提供依据．上述研究表明综合运用反射率图和垂直剖面上的一些特征来进行雹云单体的识别更为准确有效．

本文在Visual Studio 2010开发环境下，利用MFC类库实现了在生成雷达回波反射率剖面图基础上提取单体强回波(45,dBZ以上)高度、弱回波区和有界弱回波区以及探空得到的0,℃和－20,℃温度层高度等相关特征数据，并对特征数据进行了基于粗集理论的数据挖掘[8-11]，获得了区分强天气过程冰雹和暴雨的客观模型——判别规则，然后利用这些规则对天气过程进行自动识别，结果表明本文的方法能有效地识别冰雹和暴雨过程．

1 基于剖面图的特征提取

1.1 剖面图生成算法

在了解冰雹的形成机理和雹云结构特征的基础上，通过对雷达回波反射率图中回波单体的形态特征分析，确定剖线的位置，生成剖面图．本文综合3种不同的方法确定最佳剖线生成剖面图．

1.1.1 剖线的确定方法

方法1 低仰角下梯度最强处与高仰角强回波中心连线.

以0.5°仰角下的反射率图为研究基础，如图1所示为一冰雹回波，ABCD为单体所在区域的正外接矩形．对其进行扫描，获得冰雹单体最强回波的区域位置(矩形A1,B1,C1,D1)，并将其几何中心作为该冰雹单体的强回波中心(图中点O)．求出点O到矩形ABCD各边的距离，将其中最大距离(图中OE的长度)的2倍为边长，以点O为几何中心，做正方形A2,B2,C2,D2．然后以OB2作为扫描的起始位置，每10°为单位对正方形A2,B2,C2,D2进行扫描．统计每一个区域中反射率因子不小于30,dBZ的像素点数，数目最小的区域即为反射率变化梯度最大区域，分别将该区域的每一点与包含30,dBZ以上回波的最高仰角的强回波中心的水平投影位置的连线做剖线，对其剖面进行分析，特征最明显的剖面所对应的剖线即为此方法的最终剖线．

图1 剖线算法1示意Fig.1 Schematic diagram for algorithm 1

方法2 各个仰角下的强回波中心点的线性拟合.

雷达VCP21体扫示意(如图2所示)的基数据可以生成9个仰角的雷达回波反射率．

图2 雷达体扫示意Fig.2 Schematic diagram for radar scanning

检测同一个单体在不同仰角下的最强回波中心(各个仰角下单体的最大回波值可能不一样)．然后将所获得的点进行直线拟合(一般情况下只能获得前6个仰角的强回波中心，高仰角的一般无强回波)．

拟合过程：设直线方程为ykxb=+，根据直线拟合的性质，可利用公式分别求得斜率及截距，从而获得最佳剖线，其计算式为

方法3 最高仰角和最低仰角的强回波中心的连线.

找到包含30,dBZ以上强度回波反射率图中的最高仰角和最低仰角的强回波中心，连接这两点并延伸到整个单体就得到剖线．

1.1.2 剖面的生成

本文中的反射率图都包含15种颜色，分别描述了雷达回波的15个强度区间(对应关系见图3).

由于雷达基数据只包含9个仰角对应的锥面上的回波强度，且每个仰角上只存储整数方位角和整数径向距离交叉点上的回波强度信息，数据点是不连续的，要生成连续的剖面图，必须采用插值算法[12]．

雷达回波反射率的示意如图4所示，中心点O为雷达位置，线段AB为剖线，点C为剖线上任意点，图中s表示点到雷达的径向距离，θ表示方位角．

图3 反射率强度与颜色的对应关系Fig.3 Corresponding relation between colors and reflectivity intensity

图4 雷达回波反射率示意Fig.4 Schematic of radar reflectivity

首先，从剖线的左侧端点A开始，依次计算剖线上每个点的径向距离和方位角．然后，设点M(如图5所示)是剖面中与剖线上任意点所对应的空间竖直线上的点，它在雷达空间坐标系中的坐标为（r,α,h），其中r为斜距，α为方位角，h为点M在剖面图中的高度．点M位于其上下相邻的仰角为e1和e2的扫描锥面之间．M1和M2分别是经过点M的垂线与2个扫描锥面的交点，它们在雷达数据库中的坐标分别是(r1,α,e1)和(r2,α,e2)，其高度计算式为

图5 线性内插示意Fig.5 Schematic of liner interpolation

式中：0h为雷达天线的高度；ir为第i个仰角上的径向距离为等效地球半径．

点M的反射率强度的计算式为

根据上述插值算法获得的剖面如图6所示．

图6 图3中剖线对应的剖面(RCS)Fig.6 RCS for cutting line in Fig.3

1.2 强回波(45,dBZ以上)高度

雹块的增长发展过程是相当复杂的，形成冰雹的雨滴半径为4～5,mm，甚至更大．云体必须垂直发展到－20,℃温度层以上才能形成大冰雹．另一个影响降落到地面的冰雹尺寸的因素是雹块从冻结层(环境湿球温度0,℃层高度)落向地面过程中的融化效应．因此，冰雹识别首先要分析雹云单体中强回波(45,dBz以上)与相关温度层之间的关系．

在0.5°仰角图上过最强回波区域的几何中心点每隔1°做一条剖线生成剖面．检测每个剖面中45,dBZ以上强回波区域的高度，最终从180个高度数据中取得最大值作为当前时刻雹云单体强回波高度，并将其与0,℃和－20,℃温度层高度(探空测量数据)写入MySQL数据库．经过对训练样本的统计得出45,dBZ高度、0,℃和－20,℃温度层高度之间的关系，如表1所示.

表1 冰雹暴雨天气中强回波高度与0,℃和－20,℃高度之间的关系Tab.1 Relationship between strong echo height，0,℃ and －20,℃ heights in hail and heavy rain process

1.3 弱回波区和有界弱回波区

冰雹增长为大雹块的一个必要条件是有能支撑雹块的强上升气流，如果低层上升气流较强，使该处形成的降水质点被携带上升，加上风暴顶的辐散和环境风的影响，形成了低层无回波或回波强度很弱的弱回波区(weak echo recognition，WER)，当其被上方的中高层悬垂回波包围时，就形成了有界弱回波区(bounded weak echo recognition，BWER)[13]，如图7所示．

图7 WER和BWER示意Fig.7 Sketch for WER and BWER

在满足0,℃层到地面的距离比较适宜的情况下，如果回波形态再呈现出弱回波和回波悬垂特征，则产生大冰雹的可能性会明显增加，若出现有界弱回波区，则出现大冰雹的概率几乎是100%[14]．

因此本文通过图像处理[15]方法自动提取单体的WER和BWER特征，写入数据库，作为自动识别冰雹的重要特征．WER和BWER特征提取的具体操作步骤如下.

步骤1 划分区域，改变颜色．30,dBZ以下为背景色(黑色)，30～40,dBZ为黄色，40,dBZ以上为红色．

步骤2 去除干扰区域．在改变颜色的基础上，利用八邻域连通法标识各个连通域，统计每个连通域的大小，然后根据预先设定的阈值，滤除小于这一阈值的区域．去除干扰区域后的效果如图8所示．

图8 去除干扰区域后的效果Fig.8 Profile after remove interference area

步骤3 确定弱回波区或穹窿位置.

(1) 从回波墙开始，分别向两边检测30,dBZ以上区域和40,dBZ以上区域的下边界，设为RWER1(右侧30,dBZ以上区域的下边界)、RWER2 (右侧40,dBZ以上区域的下边界)、LWER1(左侧30,dBZ以上区域的下边界)、LWER2(左侧40,dBZ以上区域的下边界)，将两边边界进行比较，对较长一侧的2个边界进行特征分析．图8(a)中弱回波区在回波墙右侧，RWER1、RWER2如图所示，LWER1、LWER2为零，故取RWER1、RWER2进行分析；图8(b)中弱回波区在回波墙左侧，LWER1、LWER2如图所示，RWER1、RWER2为零，故取LWER1、LWER2进行分析．

(2) 检测右侧2个层次区域(30,dBZ以上区域和40,dBZ以上区域)下边界的宽度和高度，如果宽度大于等于4,km(经过多次调整观察确定)，则认为该边界下方区域是弱回波区．

(3) 在存在弱回波区的基础上，对边界做进一步分析：如果弱回波区上方边界在上升过程中存在下降的趋势(即存在被强回波包围的弱回波区)，并且下降部分的高度大于0.5,km、宽度大于2,km，则认为存在有界弱回波区．

(4) 如果右侧的2个层次区域都存在弱回波区或者都存在有界弱回波区，则标注30,dBZ区域所显示出的位置为弱回波区或者有界弱回波区，并计算其宽度和高度；如果只有1个区域存在有界弱回波区或者弱回波区，则按照该区域边界计算弱回波区和有界弱回波区的高度和宽度数据．

本文从3种剖线所生成的剖面中选择有界弱回波区或者弱回波区形态最明显的剖面图在屏幕上显示，并标出边界位置(如图9所示)．图9中弱回波区和有界弱回波区的高度和宽度数据计算结果如表2所示．

表2 图9中WER和BWER的特征数据Tab.2 Feature data for WER and BWER of Fig.9 km

图9 弱回波区和有界弱回波区Fig.9 WER and BWER

2 特征数据库生成

本文从某地区雷达站2005—2010年的历史资料中的24个冰雹天气过程和19个暴雨天气过程共825个样本数据进行特征提取，并将特征数据和天气类型写入数据库，得到825条数据记录，其中10条数据如表3所示．

表3 部分特征数据Tab.3 Part of the feature data

3 冰雹识别客观模型的建立

由于同一个特征的数据在不同强天气过程之间并没有严格的阈值界限，在部分特征相同的条件下可能会出现不同的天气过程或者多种天气过程同时发生，因此天气过程的识别具有不确定性[16]．本文采用可以处理不确定性数据系统的粗糙集理论对特征数据库中的特征(WER宽度、WER高度、BWER宽度、BWER高度、0,℃层高度、−20,℃层高度、强回波(45,dBZ以上)高度与0,℃以及−20,℃层高度之间的高度差和天气类型)进行分析．

粗糙集理论不需要任何附加条件就可以直接对由数据构成的决策表进行推理，它作为一种处理信息不确定、不精确、不完善系统的新的数学工具，是目前使用较多的一种归纳学习方法．

粗糙集理论认为知识是一种根据特征属性对现实中抽象的对象进行分类的能力．为了处理智能数据，粗糙集将感兴趣的研究领域对象的集合表示为知识表达系统．如果知识表达系统的属性集合分为条件属性和决策属性，那么这样的知识表达系统就称为决策表，本文所讨论的问题(见表3)即可以用决策表来表达，其中第2～9列为条件属性，第10列的天气类型为决策属性．

3.1 数据离散化

粗糙集理论只能处理离散数据，因此，在数据准备阶段首先应该完成数据的离散化，就是对连续属性的取值范围进行划分，把一个大区间划分为若干个小区间．本文采用等频离散化方法进行数据离散，该方法是一种简单的无监督单变量离散化方法．它是根据设定的频数K利用每个特征属性将测试样本分为K个子区间，每个子区间包含相同的样本数．本文将每个特征属性值划分为4个区间．

3.2 属性约简

为了获取识别模型(即分类规则)，需要将离散化的决策表进行约简．所谓约简就是化简决策表中的条件属性和属性值，使得决策表在保持原有决策能力的同时，具有较少的条件属性和属性值．这里所说的决策能力实际上是指分类能力，即依据条件属性值去判别对象类别的能力．由于决策表可看作是分类知识，因此决策表约简就是知识约简，即对知识的过滤、压缩和提炼．本文采用Johnson算法进行决策表的约简．

3.3 客观模型建立

本文根据决策表的约简结果，最终生成分类规则. 通过对43个天气过程(冰雹24个，暴雨19个)的825个(冰雹462个，暴雨363个)数据样本进行处理，得到825条特征数据，针对这些数据进行基于粗糙集理论的数据挖掘和规则提取，最终得到89条分类规则，其中8条规则如表4所示．

表4 规则表Tab.4 Rules

规则表中第1条规则表示：如果某个样本的特征同时满足WER宽度小于2,km，−20,℃温度层高度在［7.5，8.0)km区间内，并且强回波高度与−20,℃温度层高度之间的高度差小于−0.6,km的条件，则该关注对象属于冰雹的支持票数增加13票，属于暴雨的支持票数增加27票，最终按照各自支持票数与各自总样本数的比值(支持率)的大小对样本进行分类．

4 冰雹的自动识别

利用本文中获得的分类规则建立冰雹识别的客观模型，对新的风暴回波进行类型识别时，首先计算出该实况所有关注对象的特征数据，然后采用投票决策法对样本进行识别，具体步骤如下．

步骤1 目标回波区域的选择．采用连通域标记法搜索反射率强度大于等于40,dBZ的目标区域．

步骤 2 特征提取.

(1) 分别按3种剖线确定方法确定剖线位置，生成反射率垂直剖面图．

(2) 分别提取(1)中3个剖面图的弱回波区和有界弱回波区特征，将有弱回波区或者有界弱回波区特征最明显的剖面图作为最佳剖面，并将其特征数据写入特征数组feature[]．

步骤 3 规则匹配.

(1) 变量初始化．bao表示冰雹支持度变量，初始值为0；yu表示暴雨支持度变量，初始值为0；RuleNum表示分类规则总数；循环变量i＝1．

(2) 将特征数组中的特征数据依次与第i条规则匹配，如果匹配成功，则将该条规则中的冰雹支持度和暴雨支持度分别加到变量bao和yu中，并令i＝i＋1．若匹配不成功，直接令i＝i＋1．

(3) 判断i是否满足大于RuleNum，如果满足，执行步骤4，否则执行步骤3中的(2)．

步骤 4 类型识别.

将变量bao和yu分别与冰雹和暴雨的总支持票数作比值，比值大的类型就是当前回波所属的类型．

5 识别结果及分析

本文将冰雹样本作为正例，暴雨样本作为反例，经过经验丰富的气象专家采用本文冰雹识别方法对未参加训练的28个冰雹过程的383条历史样本和17个暴雨过程的262个历史样本的测试，所得识别结果如表5和表6所示．通过对表5识别结果的分析、计算，本文冰雹识别方法的命中率POD为82.77%、误警率FAR为3.06%、成功指数CSI为80.66%，其中：

表5 样本识别结果Tab.5 Recognition results of sample

表6 单个过程中的正确识别率统计Tab.6 Correct recognition rate for single monomer

24个冰雹天气过程中空报为暴雨的几乎都是在冰雹单体初步形成或者消亡阶段间断的时刻．2009年6月16日发生在某地区的冰雹(最大直径9,mm)天气过程，0,℃和−20,℃高度分别为4.2,km和7.6,km，强回波高度基本上维持在−20,℃高度层以下(如图10所示)，这些条件使得暴雨的支持度很大，导致误判全部误识别为暴雨．该过程持续时间为30,min，降雹时间只有3,min．

另外，本文方法对冰雹的短时预报的平均时效达21.6,min，具体情况如表7所示．

图10 漏报的冰雹实况Fig.10 Characteristics of hail omission

表7 预报时间提前量Tab.7 Ahead time of the forecast

6 结 论

(1) 用本文提出的确定剖线的算法能使剖面准确有效地确定出反映单体的弱回波区和有界弱回波区结构特征的剖面图．

(2) 本文提取的WER和BWER以及强波高度等特征数据的方法对于冰雹自动识别是合理有效的．

(3) 通过基于粗糙集理论的数据挖掘方法建立的客观模型和投票决策法的综合应用，能有效地将冰雹和暴雨天气区分开，实现了冰雹和暴雨的识别和短时临近预报．

(4) 本文实现了WER和BWER的自动检测，为更准确自动识别冰雹创造了条件，为获取冰雹消雹点以减轻冰雹造成的灾害奠定了基础．

［1］ 王 萍，董晓凯，贾惠珍. 雹云雷达回波反射率图像梯度及形态特征的提取[J]. 计算机应用，2007，27(增1)：372-377.

Wang Ping，Dong Xiaokai，Jia Huizhen. Hail cloud radar echo reflectance image gradient and morphological characteristics of extraction[J]. Computer Application， 2007，27(Suppl 1)：372-377(in Chinese).

［2］ Lu Zhiying，Wang Lei，Ma Hongmin，et al. Hailstone detection based on image mining[C] //Fifth International Conference on Fuzzy Systems and Knowledge Discovery. Jinan，China，2008，10(4)：39-43.

［3］ Montanya J，Soula S，Pineda N，et al. Study of the total lightning activity in a hailstorm[J]. Atmospheric Research，2009，91(24)：430-437.

［4］ 李 磊，朱 莉. “20110629”云南特大区域性冰雹天气过程分析[J]. 云南大学学报，2013，35(2)：199-207.

Li Lei，Zhu Li. Yunnan large regional hail weather proc-ess analysis for 20110629[J]. Journal of Yunnan University，2013，35(2)：199-207(in Chinese).

［5］ 江玉华，刘 娟，何 跃，等. 一次伴随冰雹的超级单体风暴特征[J]. 气象科技，2011，39(2)：172-181.

Jiang Yuhua，Liu Juan，He Yue，et al. Characteristics of supercell storm accompanied by hail[J]. Meteorological Science and Technology，2011，39(2)：172-181(in Chinese).

［6］ 覃 靖，潘 海. 具有双有界弱回波区的超级单体雹暴雷达回波分析[J]. 安徽农业科学，2009，37(5)： 2088-2090.

Qin Jing，Pan Hai. Analysis of supercell hailstorms radar echo with double bounded weak echo[J]. Journal of Anhui Agricultural Sciences，2009，37(5)：2088-2090(in Chinese).

［7］ 赵俊荣，郭金强，杨景辉，等. 一次致灾冰雹的超级单体风暴雷达回波特征分析[J]. 高原气象，2011， 30(6)：1681-1689.

Zhao Junrong，Guo Jinqiang，Yang Jinghui，et al. A super monomer of disastrous hail storm radar echo characteristics analysis[J]. Plateau Meteorology，2011， 30(6)：1681-1689(in Chinese).

［8］ Lu Zhiying，Yang Le，Li Yanying，et al. Sand-dust storm ensemble forecast model based on rough set[J].Transactions of Tianjin University，2007，13(6)：441-446.

［9］ 张 静. 基于粗糙集理论的数据挖掘算法研究[D]. 西安：西北工业大学计算机学院，2006

Zhang Jing. Theory Research of Data Mining Algorithm Based on Rough Set[D]. Xi’an：Computer College， Northwestern Polytechnical University，2006(in Chinese).

［10］ 王光宏. 数据挖掘综述[J]. 同济大学学报：自然科学版，2004，32(2)：246-252.

Wang Guanghong. Summarize for data mining[J]. Journal of Tongji University：Natural Science Edition， 2004，32(2)：246-252(in Chinese).

［11］ 张学民，安利平. 动态数据库规则获取和更新的粗糙集方法[J]. 天津大学学报，2001，34(4)：529-532.

Zhang Xuemin，An Liping. Obtaining and updating rules based on rough set theory in dynamic databases[J]. Journal of Tianjin University，2001，34(4)：529-532(in Chinese).

［12］ 潘留杰，朱伟军，周毓荃，等. 天气雷达资料的三维格点化及拼图初探[J]. 安徽农业科学，2010， 38(5)：2517-2519，2570.

Pan Liujie，Zhu Weijun，Zhou Yuquan，et al. Three dimensional attunements and puzzles for weather radar data[J]. Journal of Anhui Agricultural Sciences， 2010，38(5)：2517-2519，2570(in Chinese).

［13］ 章国材. 强对流天气分析与预报[M]. 北京：气象出版社，2011.

Zhang Guocai. Analysis and Forecasting for Strong Convective Weather[M]. Beijing：China Meteorological Press，2011(in Chinese).

［14］ 俞小鼎，姚秀萍，熊廷南，等. 多普勒天气雷达原理与业务应用[M]. 北京：气象出版社，2006.

Yu Xiaoding，Yao Xiuping，Xiong Tingnan，et al. Doppler Weather Radar Principle and Business Applications[M]. Beijing：China Meteorological Press， 2006(in Chinese).

［15］ 申俊琦，冯胜强，朱莉娜，等. 基于数学形态学的焊缝图像边缘提取[J]. 天津大学学报，2010，43(4)： 373-377.

Shen Junqi，Feng Shengqiang，Zhu Lina，et al. Edge detection of seam image based on mathematical morphology[J]. Journal of Tianjin University，2010， 43(4)：373-377(in Chinese).

［16］ 穆 穆，陈博宇. 气象预报的方法与不确定性[J]. 气象，2011，37(1)：1-13.

Mu Mu，Chen Boyu. Weather forecast method and its uncertainty[J]. Meteorological Monthly，2011，37(1)： 1-13(in Chinese).

(责任编辑：孙立华)

Hail Recognition Using Radar Echo Reflectivity Cross Section

Lu Zhiying1，Zhu Junxiu1，Tian Shuo1，Jia Huizhen2
(1. School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China；2. Tianjin Meteorological Bureau，Tianjin 300074，China)

To recognize the hail weather more accurately，852 radar base data were analyzed. These data were collected from 24 hail processes and 19 heavy rain processes in a certain area between 2005 and 2010. Automatic identification model for hail was built by using radar reflectivity cross section. The optimal cutting lines were determined based on the analysis on the forming mechanism and structure of hail cloud. Radar reflectivity cross section was obtained by using interpolation algorithm and was analyzed with image processing methods. Features data such as the height differences between strong echo(over 45,dBZ)and 0,℃ as well as －20,℃ isothermal layer，and the height and width of weak echo recognition and bounded weak echo recognition were obtained. Rules for hail and heavy rain recognition were acquired by using data mining based on rough set，and automatic identification model was built. 385 hail samples from 28 hail processes were distinguished with this model. The test results show that the accuracy of recognition rate is 82.77%. This provides an effective method for identification and short-time forecast for hail，and it is helpful for reducing the loss caused by hail.

image processing；weak echo recognition；bounded weak echo recognition；rough set；data mining

P412.25

A

0493-2137(2015)08-0742-08

10.11784/tdxbz201401022

2014-01-08；

2014-03-16.

天津市自然科学基金资助项目(14JCYBJC21800)．

路志英（1964— ），女，博士，副教授.

路志英，luzy@tju.edu.cn.

时间：2014-03-14.

http://www.cnki.net/kcms/doi/10.11784/tdxbz201401022.html.