天气雷达反射率因子图像中三体散射自动识别

王 萍，杜雪峰，徐考基

(天津大学电气与自动化工程学院，天津 300072)

天气雷达反射率因子图像中三体散射自动识别

王 萍，杜雪峰，徐考基

(天津大学电气与自动化工程学院，天津 300072)

三体散射的准确识别可以有效降低强冰雹的误警率、提高预报强冰雹的时效性．在详细分析三体散射产生机理的基础上，将三体散射按照其图像特点定义成 2种类型．设计标准模板，在将搜索区域进行标准化变换的配合下成功解决了楔状的、短的三体散射的自动识别；用链码技术将走向一致的区域边界的比例定量地描述出来，以此解决长直的三体散射的自动识别．实验表明，该方法对短三体散射的区域断裂和长三体散射的远端粘连具有不敏感性，对三体散射的总击中率达到94.8%，误警率不足1%．

三体散射；楔状识别；链码特征提取；冰雹预警

天气雷达是对强对流天气进行监测和预警的主要工具之一．这种雷达按照预定的几个仰角分别依次旋转 360°向空中发射脉冲形式的电磁波，当电磁波遇到雪花、雨滴或冰雹等降水物质时，大部分能量会继续前进，而有一小部分能量被降水物质向四面八方散射，其中，向后散射的能量被雷达接收[1]从而探测到这些降水物质的存在．为方便观测，雷达系统提供一种反射率图像以直观展现探测到的所有数据信息．1987年，Zrnic[2]在雷达反射率图像中首次发现了“三体散射(three-body scatter spike，TBSS)”现象，即一种与强冰雹回波相关的、出现在强回波外侧、沿径向向外延伸的、弱的、窄带状的回波，并称其为“三体散射特征”，后来，Wilson等[3]称其为“火焰回波(flare rcho)”、Lemon[4]称其为“三体散射长钉”．根据 Zrnic[2]的理论，雷达发出的电磁波遇到强冰雹类降水物质时，向四周散射的能量会很强，于是使向下散射到地面的部分能量有能力返回到原降水物质再回射到雷达，从而使雷达误认为在冰雹的径向外侧存在着弱的降水物，因此在强冰雹回波的径向外延成虚假的、弱的回波．

国内外学者对三体散射与冰雹的关系做了深入的研究，廖玉芳等[5]对伴随产生TBSS现象的中国11个强对流过程中的 23个强风暴进行了统计，每一个强风暴都产生了直径超过 2,cm的强冰雹，其中有 4个强风暴产生了直径 10,cm 以上的冰雹．因此得出“如果 S波段新一代天气雷达产生三体散射现象，则降强冰雹的可能性几乎是 100%，同时大多伴随地面灾害性大风”的重要结论．据湖南省北部的统计，在2002—2005年间的19个雹暴中，有16次出现了TBSS，占84%，且 16次所降冰雹的最大直径均超过2,cm[5]．更重要的，TBSS的持续时间多为 30～60,min，其中的大部分时段位于地面降雹之前，因而使TBSS预警强冰雹具有时效性．另外，Hubbert等[6]计算了三体散射区域的偏差反射率，进而又通过对地面的反射建立了统计和解析模型，他们还获得了冰雹单体的反射系数并对冰雹单体进行估计. Smallcomb[7]对三体散射速度图进行了分析，说明了三体散射区域常常伴有较大的速度值，并且其并不是真正速度，有时会导致中气旋算法的误警，影响预报质量．Zrnic等[8]还证明三体散射与冰雹尺度有关．如此可见，及时扑捉到TBSS有利于提高强冰雹预警的水平和质量．

我国新一代天气雷达的冰雹指数产品采用的是由美国开发的强冰雹预警算法，该算法主要借助扩展到0,℃层以上、达到或超过45,dBz的强反射率因子，累计计算一种“冰雹动能通量”[9]，再形成基于发生于美国的若干强冰雹实况的冰雹指数[10]．据统计，该冰雹指数在对强冰雹平均命中70%的同时，平均误警率高达 49%[5]．究其原因，是因为这个冰雹指数考虑的因素并不主要为强冰雹所独有，许多强降水回波也会因跃于 0°层甚至-20°高度以上而获得很高的冰雹指数．

TBSS应该是强冰雹雷达回波所特有的现象．如何通过计算机自动捕捉三体散射现象在国内外都未找到相关报道，笔者从TBSS的形成机理和回波特点出发，利用图像识别和计算机视觉技术构建并实现专用的TBSS识别算法，以代替人工凭借经验的低效率查找，为提高具有TBSS现象的强冰雹的识别准确性和预警时效性提供重要的客观依据．

1 三体散射现象及其图像特点

根据 Zrnic[2]对三体散射的定义，三体散射的形成经过了3次反射，设雷达发射的电磁波在离开雷达径向距离r后遇到强冰雹A，如图1所示，此后，电磁波的大部分能量继续向前传输，剩余部分向周围散射，其中，散射波的一小部分返回雷达，以示探测到目标 A，另有一部分散射向地面，再形成返回波到达冰雹A，然后有部分能量经过A的反射又回到雷达．容易看出，当冰雹 A的强度、位置一定时，源于冰雹 A的散射波到达地面的路径越长，回到冰雹 A的波的能量就越弱，设只有落到图1中以BC或BD为半径的圆形区域内的返回波剩余能量才能够支撑其行进至雷达并不被解读为噪声，这时，雷达会将这部分途径地面的返回波信号误读为在沿电磁波发出方向、远离冰雹 A的一定距离范围内存在着弱的降水物，而该距离理论上大于等于 A到地面距离 h，即有

图1 TBSS形成原理示意Fig.1 Forming principle of TBSS

在雷达反射率因子图像上，TBSS的特征是明显的和有规律的，具体情况如下所述．

1.1 产生TBSS的风暴具有高反射率强度的核区

Lemon[4]指出要产生能够分辨的 TBSS，风暴核反射率因子必须大于 63,dBz，文献[5]通过对我国境内499例TBSS的统计，发现产生TBSS的风暴核反射率因子大于等于60,dBz的比例高达99.2%．

1.2 TBSS出现的位置具有规律性

根据图1所示TBSS的形成机理，TBSS一定出现在雷达与风暴核区连线的延长线上，另据我国夏季多为西风或偏西风，产生TBSS的风暴大多位于探测雷达的西侧或偏西侧[3]．例如，图2为2张实测的显现TBSS的反射率图像，其中，分别在图2(a)第2象限(西北)和图 2(b)第 3象限(西南)径向射线上、高反射率核区外侧的椭圆标记处出现了TBSS．

1.3 TBSS具有特定的形态

TBSS具有直-带状(见图 2(a))或楔-带状(见图2(b))的形态，其长度与风暴核区的反射率强度、面积以及地面条件有关，据文献[5]统计，TBSS的长度通常小于15,km，超过3,km的比例高达91.6%．

图2 TBSS图像特点Fig.2 Image feature of TBSS

1.4 TBSS具有弱的反射率强度

由于雷达探测波能量历经雷达-冰雹物-地面、再到冰雹物-雷达的较长路途的消耗，其反射率强度一般很低，据资料统计，TBSS的反射率强度不超过25,dBz，更多的低于20,dBz．

简言之，TBSS是一块位于强反射率风暴核区径向外侧的、直-带状或楔-带状的弱回波区域．

1.5 TBSS在径向速度图的特性

Lemon[4]提出沿着三体散射回波伸展方向的径向速度是关于冰雹单体 A径向和垂直运动的函数，这个函数接近 sinw uθ+ ，w为单体的垂直运动速度，u为散射位置的径向速度分量，但是由于地面上物体的运动，三体散射的径向速度图会受到污染．冰雹的形成在零度层以上有大量的温度很低的水汽时，强大的上升气流促进了大冰雹的形成．

多数情况 w是负数，只有当有强大上升气流时其才为正数，这样，根据速度图上三体散射的信息，可以看出此时冰雹是上升还是下沉，文献[11]对这 2种情况进行了详细的讨论，并且说明了三体散射可能产生速度模糊．一般情况下，三体散射长钉上都会形成比较大的负速度，图3(b)为2012年4月16日广东的基数据中三体散射的径向速度图，可以看到部分核区及部分三体散射区域相对周围有明显的负速度．

图3 TBSS的雷达反射率、径向速度及其剖面Fig.3 Reflectivity，radial velocity and their sections of TBSS

1.6 TBSS在径向速度图和反射率图上剖面图特征

图 3(c)为图 3(a)反射率图的剖面图，可以看到三体散射在4,km高度以上朝单体外侧斜向上伸出，并且随着与核中心距离的增加，其反射率值逐渐降低.

图 3(d)为图 3(b)径向速度图的剖面图，可以发现图中三体散射区域速度跳变相对较大，且有较大的负速度点(浅绿色)，三体散射也同样在 4,km以上伸出来．

2 TBSS识别算法

算法设计采用“在可能存在TBSS的区域内寻找满足TBSS图像特点的弱回波”的指导思想．

2.1 关注区域的确定

设天气雷达采用九仰角的降水探测模式，鉴于6°及其以下探测仰角的反射率图上均有发现 TBSS的报道，以下方法适合于 6°及其以下探测仰角的反射率图．

分别在这些仰角(α)的反射率图上，找到达到或超过60,dBz的区域，计算该区域中心的径向距离r和距地面的高度h，用该区域中心径向外侧h之后的2段圆弧和向径向线2侧各偏离10°的2条射线圈定一个搜索 TBSS的扇形区域，如图 4所示，其中，参数 h为风暴核区中心高度(nn× )、扇区厚度Δr依据统计值设定．

图4 可能存在TBSS的范围Fig.4 Area that TBSS probably exist

2.2 定义“楔-带状”和“直-带状”

TBSS大体形如带状．当其长度介于 3～10,km之间时，其根部相对较粗、头部相对较尖(见图 2(b))的特点突出出来，本文将其定义为“楔-带状”；而长度超过10,km的TBSS，其“直”且“等宽”的部分占主导(见图 2(a))，本文称其为“直-带状”．据文献[3]统计，前者形状的 TBSS占 54.5%，后者占37.1%，其余 8.4%的长度小于 3,km，这时，图像上TBSS仅为 2个像素左右，“带”状特点已经消失，且易混淆于具有起伏特点的回波区域边界之中，因此，本文从长度超过3,km以上TBSS的特点出发进行算法设计．

2.3 “楔-带状”TBSS识别算法

匹配法是检测具有特定形状的对象的常用方法，如模板匹配法[12]、基于特征点的频率匹配法[13]等．其中，模板匹配通过计算相似度函数搜索图像中对模板的最佳匹配点，鉴于气象雷达的大数据量特点且快速运行的要求，不宜直接搬用已有模板或方法．为此，本文提出一种基于搜索区域标准化处理的楔状模板匹配算法用于快速找到匹配点，判断是否存在“楔-带状”的TBSS．

2.3.1 构建“楔-带状”标准模板

根据TBSS的长度范围，设定正方形模板(nn× )的尺度为

图5 标准模板Fig.5 Standard temp template

2.3.2 标准化搜索区域

为适应图5所示的标准模板，首先将扇形的搜索区域二值化，即对区域内在0～25,dBz范围内取值的像素点记为 1，其余的记为-1，形成区域 S的函数fS(i, j)；再对 fS(i, j)进行旋转变换得到区域S′，描述函数为 fS′(i , j )，旋转角度等于区域 S的主轴到正西向射线的夹角，如图 6(a)所示；最后按照图 6(b)方法将搜索区域S′标准化．

图6 搜索区域旋转变换及标准化搜索区域Fig.6 Target area transform to the standard searching

2.3.3 模板匹配

令模板 (,)g i j在区域BS′上按照式(2)进行遍历搜索运算，得到判定值 (,)R x y，若存在满足式(3)所述条件的判定值，则判定该区域出现“楔-带状”的TBSS．

由于可能存在 TBSS的区域被二值化和模板的三值化布局使得该模板遇到“楔-带状”的 TBSS时，对象部分和背景部分被同时反向加权，可获得搜索效率的提升，且 (,)R x y越大，越接近理想的“楔-带状”，其中， (,)R x y的最大值为29，阈值M在理论分析和大量样本统计的基础上确定．

2.4 “直-带状”TBSS识别算法

当 TBSS较长时，由于其形状趋于平直，使其区域边界的多数元素呈现出与 TBSS方向相关的主流走向，因此，可以借助对边界线的描述和对边界线各元素“同方向数”的统计，实现对“直-带状”TBSS的识别．

2.4.1 提取区域边界

首先对关注区域进行二值化和补空洞处理再用Canny算子提取出关注区域(0～25,dBz)联通域的边缘，这些边缘将组成区域的边界L．

2.4.2 描述区域边界

取图7所示的8方向基元组描述边界线段，即顺次扫描搜索区域BS′得到其边界线的一个点，令该点为起点进行边界线的追踪，同时，以级联的方式将相关基元不断追加到L中，最终形成该区域边界线的8方向链码．

图7 8方向链码基元Fig.7 Chain code of eight directions

2.4.3 计算边界链码的方向数

设边界线的链码总长度为 NL，链码中，设方向 i和j上的基元数分别为num(i)和num(j)，其中，方向i和j满足关系式

则边界线的i-j方向数ijN 定义为

式(5)中，方向i的选择与搜索区域的主轴方向β(正北为0°)有关，即

当ijN 大于阈值(0.4)时表示 i-j方向的边界元素占主导，这时判定该区域存在“直-带状”的TBSS．

2.5 算法流程

本文根据 TBSS的形状特点提出 2种识别算法．考虑到“楔-带状”更具普遍性，采用首先使用“楔-带状”算法、再补充使用“直-带状”算法的 2级识别策略，如图 8所示．如此，对于根部具有“楔”状特点的“直-带状”TBSS，也只需经一级识别即可．

图8 TBSS识别算法流程Fig.8 Flow chart of identification algorithm of TBSS

3 测试与分析

按照本文提出的方案和算法，对广东、江西、天津共 144个含有三体散射的实例进行测试．在 144个实例中，涉及 0.5°～6.0° 6种仰角的反射率图像，径向范围从25,km到200,km，包含“直-带状”TBSS共27个，占18.75%，“楔-带状”TBSS共117个，占81.25%．又随机选择了含有强核区但未出现TBSS的反射率图像 150张，对本文算法展开反例样本测试．总测试结果如表1所示，根据表1可得

识别TBSS的CSI评分：

图9～图12给出了测试示例．

图9中的3组图片是“楔-带状”算法的识别结果，小窗口展示的是寻找到目标区域并判别为三体散射的情况．图 9(a)和(b)属于楔状 TBSS，前者较长、后者较短，均被正确识别出来；而图 9(c)属于“直-带状”的 TBSS，但从小窗口展示的关注区域看，其根部也呈现出楔状，于是使模板匹配时获得高值，因此对此类样本无需再启用“直-带状”算法，即可完成TBSS的识别．

表1 算法测试统计结果Tab.1 Statistics results of the test

图9 “楔-带状”TBSS识别算法的识别结果示例Fig.9 Identification results of the cuneiform-banding TBSS

基于专用标准模板的匹配法除了能应对较短TBSS、根部呈楔状的“直-带状”TBSS以外，对断裂噪声也具有一定的抵御能力，例如，图10中的TBSS在尖端附近出现了断裂，但其“楔-带状”未遭破坏，算法依然能将其识别出来．

图10 有断裂的“楔-带状”TBSS的识别结果Fig.10 Identification results of the breaking cuneiformbanding TBSS

图 11反射率图中出现的 TBSS显然是“直-带状”的，由于其过长，并未完全被包含到所提取的关注区域中，其根据式(6)求出借助图11给出的区域边界链码不难算出：于是判断存在 TBSS，识别结果用方形框标注于图 11链码下面的图中．

图11 “直-带状”TBSS识别算法的识别结果示例Fig.11 Identification results of straight-banding TBSS

在9张漏识的TBSS中，基本涉及3种情况．其一，也是最主要的情况(约占6成)，就是TBSS较短且过粗，致使模板匹配时无法获得高值，如图 12(a)所示，或是 TBSS短到 3,km甚至 3,km以下，如图12(b)所示；其二，就是TBSS与远端降水回波发生混合粘连，导致识别失败，如图 12(c)、图 12(d)所示；第3种情况仅出现了一例，就是TBSS不规范，既不像“楔-带状”、也不是“直-带状”，如图 12(e)所示．但是在以上各例中仰角变化时，能够找到可被算法识别的三体散射，产生预警，所以基于多仰角的算法探测有益于弥补上述单仰角对三体散射的漏识问题，被本文算法误警为 TBSS的样本如图 12(f)所示，这是使用“楔-带状”算法的误警结果．同时还可以通过其雷达反射率剖面图以及径向速度图来辅助判断是否为三体散射现象，见图 13，发现其形状向左伸出且又向上生长，与图 3(c)相差较大，伸出的低反射率回波高度位于 3,km以下，不满足三体散射雷达反射率剖面图特性，并且其速度图中，误警区域速度不存在较大变化或绝对值很大的负向速度，可见其为误识，所以通过对反射率剖面图及速度图的分析，可以协助进一步确认识别三体散射的正确性．

图12 漏识和误警Fig.12 Failure identification and false alarm

图13 误警产生时三体散射的反射率剖面和径向速度Fig.13 Reflectivity section and radial velocity of the false alarm TBSS

4 结 语

TBSS是当天气雷达探测到即将发生的强冰雹可能出现在反射率图像中的特殊现象，其可能出现的条件是明确的，可能出现的部位是确定的，但其形状、长短、强弱又是不严格一致的．本文就其自动识别问题做了以下工作．

将 TBSS归纳成两类形状，即“楔-带状”和“直-带状”，根据各自的特点设计不同的识别算法．但在使用时，并不首先做形状识别，而是采用一种“主-协”识别策略，即“楔-带状”算法做“主”识别，这时，凡根部到中部呈现“楔状”的TBSS，即便因 TBSS很长，其根部的楔状特征已不突出的情况，也可将其识别出来，用作协助识别的“直-带状”算法旨在“楔-带状”算法无功而返的时候投入工作，将那些细长的TBSS识别出来．

为识别“楔-带状”TBSS，设计出专用的、反映这类 TBSS特点的、用于匹配搜索的标准模板，配合搜索区域的标准化处理，能够 100%地将 3,km以上、未与其他回波发生粘连的“楔-带状”TBSS识别出来，同时能够在一定程度上抵御断裂、长度变化等干扰因素的影响．

当TBSS很长时，TBSS区域边界的多数成分将与TBSS的取向一致，本文借助边界链码将这一特征定量地计算出来，由于该方法侧重的是取向一致的比例，因此在识别“直-带状”类TBSS的时候，能够抵御远端处粘连上其他回波的干扰．

本文的方法计算速度快，时间复杂度为 ()O n，其中，模板和链码的算法均为线性复杂度，且问题规模不大，即使在计算多个仰角的情况下，对一个时刻的处理也在毫秒级．这种综合考虑各种情况的方法比单一使用模板判别或链码判别有更高的识别率和鲁棒性，且依据三体散射出现位置对其进行旋转，简单有效地解决了识别形状的旋转问题．

为说明本文方法的有效性，展开了一定规模的正反例样本的算法测试．测试结果表明，本文方法能够将不与其他回波粘连的“较短”TBSS及直的“较长”TBSS识别出来，满足算法条件的个别降水回波形状有可能被识别为三体散射，这时通过对雷达反射率和速度剖面图的分析可以进一步判断其真伪．

人工识别三体散射时，因为易受到预报员经验水平、工作状态及各种人为因素干扰，而造成识别不准确，工作效率低；自动识别的技术具有准确率高、工作效率高且节省人力等诸多优势．

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

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

［2］ Zrnic D S. Three-body scattering produces precipitation signature of special diagnostic value[J]. Radio Science，1987，22(1)：76-86.

［3］ Wilson J B，Reum D. The flare echo：Reflectivity and velocity signature[J]. Journal of Atmospheric and Oceanic Technology，1988，5：197-205.

［4］ Lemon L R. The radar “three-body scatter spike”：An operational large-hail signature[J]. Weather and Forcast ing，1988，13：327-340.

［5］ 廖玉芳，俞小鼎，吴林林，等. 强雹暴雷达三体散射统计与个例分析[J]. 高原气象，2007，26(4)：812-820.

Liao Yufang，Yu Xiaoding，Wu Linlin，et al. Statistic and case studies on radar three body scattering of severe hailstorm [J]. Plateau Meteorology，2007，26(4)：812-820(in Chinese).

［6］ Hubbert J C，Bringi V N. The effects of three-body scattering on differential reflectivity signatures[J]. Journal of Atmospheric and Oceanic Technology，2000，17(1)：51-61.

［7］ Smallcomb C. Hail spike impacts on Doppler radial velocity data during several recent lower Ohio Valley convective events[C]//23rd Conf on Severe Local Storms. St Louis，MO，USA，2006：1-22.

［8］ Zrnic D S，Zhang G. Three-body scattering and hail size[J]. Journal of Applied Meteorology and Climatology，2010，49(4)：687-700.

［9］ Waldvogel A，Schmid W，Federer B. The kinetic energy of hailfalls(part I)：Hailstone spectra[J]. Journal of Applied Meteorology，1978，17：515-520.

［10］ Donavon R A，Jungbluth K A. Evaluation of a technique for radar identification of large hail across the upper[J]. Weather and Forecasting，2007，22(10)：244-254.

［11］ Stan-Sion A. Three body scatter signature producing uncertainties on Doppler radial velocity structures [C]// 33rd Conference on Radar Meteorology. Australia，2007：5. 3.

［12］ 余立功，王 强，陈 纯. 多尺度模板匹配算法[J].工程图学学报，2005(3)：80-83.

Yu Ligong，Wang Qiang，Chen Chun. On the multilevel template match algorithm[J]. Journal of Engineering Graphics，2005(3)：80-83(in Chinese).

［13］ 刘 睿，姚素英，史再峰. 基于特征点的图像配准子系统的研究与实现[J]. 电路与系统学报，2012，17(3)：10-14.

Liu Rui，Yao Suying，Shi Zaifeng. Image matching circuit subsystem based on the interest points[J]. Journal of Circuits and Systems，2012，17(3)：10-14(in Chinese).

(责任编辑：孙立华)

Automatic Identification of Three-Body Scatter Spike in the Doppler Radar Reflectivity Image

Wang Ping，Du Xuefeng，Xu Kaoji
(School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China)

The accurate identification of three-body scatter spike(TBSS) can reduce the false alarm rate of strong hail and improve the timeliness of strong hail forecast. This paper provides a detailed analysis of the formation of the threebody scatter spike and then defines it into two classes. Designing the standard template，normatively transforming the target area solve the problem of short，cuneiform-shaped TBSS identification；and then the chain code is used to calculate the main direction of the boundary for identifying the long and straight TBSS. Test shows that the method is not good at the identification of the short TBSS with broken region or the long TBSS whose end is connected with others. Identification rate of TBSS is 94.8% and false alarm rate is less than 1%.

three-body scatter spike(TBSS)；cuneiform-shaped identification；feature extraction from the chain code；hail warning

TP391

A

0493-2137(2014)08-0711-08

10.11784/tdxbz201305037

2013-05-17；

2013-06-19.

中国气象局气象探测中心2012年新一代天气雷达建设业务软件系统开发及应用(ZQC-H12085).

王 萍（1955— ），女，教授，博士生导师，wangps@tju.edu.cn.

杜雪峰，dxf_friendship@sina.com.

时间：2013-09-11.

http：//www.cnki.net/kcms/detail/12.1127.N.20130911.0923.001.html.