引言

冰雹灾害作为我国主要自然灾害之一，对农业生产和人民生命财产安全有着严重的影响，人工防雹是目前防灾减灾工作的一项重要措施。随着我国新一代天气雷达组网建设，其在冰雹等强对流天气事件中的监测和预警作用越来越重要(俞小鼎等，2005；李德俊等，2011；王秀玲等，2012；周泓等，2014)，雷达探测的高时空分辨率特性，使我们能获得更多有效的天气信息，提高了冰雹识别和预警的准确性，并逐渐成为我国人工防雹的主要参考依据(王瑾等，2011；张秉祥等，2014)。如何充分利用雷达探测到的对流性单体的结构特征，了解对流单体中的云水分布，对人工影响天气有着重要的影响。从多普勒天气雷达组网以来，人们通过日常工作的积累和研究，得出了一些冰雹云识别的重要指标。Winston等(1986)研究表明，垂直积分液态水含量对冰雹的存在有较好的指示作用；李红斌等(2010a)利用各类雹云在初期和发展阶段的回波强度、30 dBz强回波中心高度和顶高、VIL等雷达主要参数特征及跃增变化，由初始回波和雷达参数的跃增变化总结了识别雹云的技术方法; 付双喜等(2004)利用雷达产品VIL进行了冰雹云识别，发现VIL值越大预示着降水强度越大，当VIL值>7.0 kg·m^-2时，将预示有强对流天气特别是冰雹天气的出现。以上研究表明，利用雷达探测资料，分析对流性单体的结构特征及其演变，有利于进一步提高人工防雹作业和效果评估的科学性和准确性。

近年来，一些省、市经过多年的努力，也逐步建立了适合当地条件的人工防雹作业指挥系统(刘春文等，2009；李红斌等，2010b；邹书平，2011)。但从目前的实际应用来看，人工防雹作业主要还是利用雷达回波的特征参数，结合经验给出识别指标，这些方法简单实用，但存在一定的局限性，即经验性太强，而且主要基于雷达的二次产品数据，受某一因子的影响比重较大，自动化程度低，大多靠预报员主观性识别得到，工作量大，对防雹作业的敏感区认识不够，造成作业效果差强人意。如何将雷达探测到的回波信息与冰雹云的结构特征有机结合起来，提高冰雹云的识别能力，提供相对准确的作业部位，为本文研究的重点。冰雹天气不同于大风、强降水等强对流天气，在雷达回波上往往表现为反射率因子较强，强回波中心的相对高度较高(Witt et al, 1984; 胡胜等，2015)。因此，本文以风暴识别为核心，利用风暴识别算法，提取风暴各类结构特征参数，并以其特征参数建立决策树模型，对识别出来的风暴按强度进行分类，降低部分发展相对不旺盛的对流单体造成冰雹云识别的空报率，提高作业选择的可决策性，最后集成人工影响天气最新技术研究成果，实现了具有高度自动化、科学决策的人工防雹作业指挥系统。通过对2006—2014年重庆、辽宁大连和河南三门峡的风暴个例识别效果检验，证明该方法可以提高强风暴的命中率，降低其空报率，能较有效提高人工防雹作业指挥的效率，并在大连、新疆、黑龙江和河南等地的人工防雹作业指挥系统中得到了实际应用。

1 资料与研究方法

本文选取辽宁大连、河南三门峡和重庆地区2006—2014年的多次灾害性天气过程，收集整理得到31次较为典型的冰雹天气过程(同时具有雷达观测和冰雹观测记录)，累计182站次冰雹样本。雷达资料采用辽宁大连、河南三门峡、重庆主城和万州4部多普勒天气雷达基数据，其中辽宁大连为C波段雷达，其余3部雷达为S波段。由于探空资料在时间和空间上分辨率都较粗，所以分析时自动选取距离雷达最近的测站，探空资料的时间为距离雷达探测时间最近的一次探空。冰雹资料采用的是辽宁市、河南省和重庆市人工影响天气办公室提供的炮点观测记录和有冰雹观测的测站观测记录，包括冰雹出现的时间、地点和最大冰雹直径等数据。

目前风暴跟踪识别方法主要分为两类，一种是利用单阈值的TITAN(Thunderstorm Identification, Tracking, Analysis, and Nowcasting)(Dixon et al, 1993)算法，一种是利用多阈值的SCIT算法(Johnson et al, 1998)。TITAN侧重于风暴的形态识别，能处理风暴的分裂和合并，而SCIT侧重于风暴单体质心位置和风暴体结构的识别。考虑到目前SCIT算法为我国CINRAD-PUP的主要算法，已在业务中大量使用，所以本文首先选择SCIT算法，对每个风暴单体的特征参数进行识别和跟踪，确定每个风暴单体的生命周期、移动路径、风暴顶高、风暴体积、强中心强度和高度、基于风暴单体质心的垂直累计液态含水量(VIL)等风暴结构主要特征参数；然后采用决策树分类模型方法(Kang et al, 1997)，根据不同强度风暴单体的概念模型(俞小鼎等，2006)及当地日常业务总结的风暴识别指标，选择风暴结构特征参数作为决策因子，建立决策树模型，采用回归分析方法将风暴单体分为雷雨云、单体风暴、多单体风暴和强风暴(俞小鼎等，2006)，风暴分类的定义原则上依赖于风暴的强度和形态等物理特征，例如风暴体积、质量等。雷雨云一般强度较弱、质量较小，不考虑冰雹存在，而单体风暴、多单体风暴和强风暴均有可能是普通风暴或强风暴，只是风暴的体积、质量和强度、垂直切变、强中心位置的高度等有一定的差别，其出现冰雹的概率相对单体风暴，强风暴概率更大；最后利用GIS技术，结合风暴移动路径及其预报，对下游地区作业点根据其距离风暴体的远近和风暴体的位置、强度和作业敏感区的高度，自动输出预警或作业参数，指挥作业站点人员进行防雹作业。

2 风暴自动识别技术

正确识别冰雹云是进行人工防雹作业的前提。如何将冰雹云和雷雨云区分开来是解决选择人工增雨或防雹作业的关键问题，从现有的研究来看，雷达识别冰雹云的方法大致包括：通过回波形态主观识别冰雹云；通过回波参数客观识别冰雹云；通过回波形态和参数相结合的主客观方法识别冰雹云。

考虑到冰雹天气总是伴随着较强的风暴单体出现，本研究以风暴单体识别为基础，改进传统回波形态和特征参数识别的方法，把重心转移到可能出现冰雹天气的风暴识别上来。目前我国业务雷达软件使用的风暴识别算法采用的是SCIT算法，该算法由四个子功能组成：风暴单体段、风暴单体质心、风暴单体跟踪和风暴位置预报。SCIT算法输出的单体属性在一定程度上能反映风暴的发生、发展过程，但对于冰雹云的识别，尤其是人工防雹作业指挥来说还不够，吴林林(2006)详细分析了美国NEXTRAD的冰雹探测算法，重新对算法进行了实现。本文在此算法的基础上，结合人工防雹作业敏感参数，增加了0℃层高度和-20℃层高度对应的风暴中心强度、风暴反射率因子垂直递减率、回波强度30和45 dBz所在高度、中高层闭合强中心所在高度等客观参数。

2.1 0℃层高度和-20℃层高度对应的风暴中心强度

在探空分析中，0℃层和-20℃层是冰雹分析中最引人关注的两个特征高度层，以往简单有效地判断有无大冰雹的方法是根据强回波区相对于0℃层和-20℃层等温线高度的位置。强回波区必须扩展到0℃等温线以上才能对强降雹的潜势有所贡献(俞小鼎，2014；濮文耀等，2015)。当强回波区扩展到-20℃等温线高度之上时，对强降雹的潜势贡献最大(Witt et al, 1998), 为表征这一特征，在此定义Z_0℃和Z_-20℃分别表示0℃层高度和-20℃层高度对应的回波强度，由于雷达探测的限制，随着距离增大，水平分辨率和垂直分辨率降低，在距离雷达100~200 km的区域需要对数据进行插值，本文采用双线性插值的方法进行插值。

2.2 风暴反射率因子垂直递减率

为了表征风暴的垂直结构，特定义风暴反射率因子垂直递减率 ${\nabla _Z}Z$ ：

${\nabla _Z}Z = \left( {{Z_i} - {Z_{i - 1}}} \right)/\left( {{H_i} - {H_{i - 1}}} \right)$

(1)

式中，Z_i为风暴第i层分量中心回波强度，单位为：dBz；H_i为第i层分量中心回波强度对应的高度，单位为：km。该参数反映了风暴的垂直结构特征，可以很好地表征风暴的悬垂结构。回波强度垂直廓线的位置分别为风暴二维分量中回波强中心所在的位置，当该位置处强度垂直递减率出现突变，而且变化较大时，往往是风暴悬垂结构所在的位置。

2.3 回波强度30和45 dBz所在高度

从统计数据来看，30和45 dBz是风暴分别达到雷雨云和冰雹云的临界指标，Smith等(1975)提出，云内最初冰雹直径约为0.4~0.5 cm，利用雷达气象方程关系式推到得出，冰雹云初期等效雷达反射率因子约为45 dBz。为提取这一特征参数，在风暴识别算法中分别定义H₃₀和H₄₅来表征回波强度30和45 dBz所在高度，通过风暴各层水平分量的强度进行线性插值得到，该指标作为风暴分类的一个特征样本。

2.4 中高层闭合强中心

回波顶高和强回波中心是人工防雹作业的重要参考特征参数，对于单体回波来说，强回波中心往往只有一个，但多个单体组成的雹云回波，其垂直高度上往往存在多个强回波中心。而回波顶高是雷达探测到的以最高仰角为基础的回波顶高度，并不一定能反映冰雹云的中高层特征。由于防雹作业关心的多数是冰雹生长过程中，中高层的强回波中心(许焕斌等，2008)，为科学定位作业的方位和仰角，提高作业效率，在风暴识别的过程中，特定义Z_Hm和H_m分别表示中高层闭合强中心的强度和高度，以区分低层强度较强的回波。其中中高层闭合强中心的定义为：从风暴的顶部到底部对风暴各高度层的强中心进行搜索，当风暴各层强中心强度出现先增大，然后逐渐减小的过程，则认为该转折点为中高层闭合强中心位置，即当风暴体垂直高度上各层强回波中心强度随高度呈双峰或多峰分布时，则高度最高的强中心位置为中高层闭合强中心。

3 基于决策树分类模型的风暴分类方法

数据分类是数据挖掘中的一种分析方法，通过学习训练集构造一个分类函数或分类模型，该函数或模型能够把数据记录映射到给定类别中的某一个，从而可以应用于数据预测。数据集由一组数据库记录构成，记录形式可以表示为(V₁, V₂, …, V_n, C)，其中V_i表示样本的属性值，C表示类别。分类模型的构造方法包括统计方法(如贝叶斯方法)、机器学习方法(如决策树方法)、神经网络、粗集算法和遗传算法等。本文采用基于二分法的CE4.5决策树分类算法(Ruggieri, 2002)进行风暴的自动分类。算法中利用风暴识别输出的各类风暴特征参数作为样本属性值和风暴类型作为类别，构建决策树模型(如图 1)，最后根据决策树模型按风暴属性对识别出的风暴单体进行自动分类。如下所示：

$样本:X = \{ {x_1},{x_2}, \ldots ,{x_n}\} $

(2)

$结果:Y = \{ {y_1},{y_2}, \ldots ,{y_m}\} $

(3)

式中，x_n为风暴形状、回波强度、动力学特征和热力学特征等18个风暴特征参数，并根据不同风暴的概念模型及重庆2009—2014年雷达观测统计得到决策参考值，如表 1，y_m为雷雨云、单体风暴、多单体风暴、强风暴等4类风暴类型。

分类过程按照图 1，采用CE4.5算法进行逐级判别，直到所有风暴单体得到确定的分类结果。

CE4.5算法的主要原理：

(1) 初始化决策树T为只含一个树根(X, Q)，其中X是全体样本集，Q为全体属性集。

(2) if(T中所有叶节点(A, X’)都满足X属于同一类或Q’为空)then算法停止；

(3) else{任取一个不具有(2) 中所述状态的叶节点(X’，Q’)}；

(4) for each Q’中的属性A do计算信息增益gain(A，X’)；

(5) 选择具有最高信息增益的属性B作为节点(X’，Q’)的测试属性；

(6) for each B的取值b_i，do{从该节点(X’，Q’)伸出分支，代表测试输出B=b；求得X中B值等于b_i的子集x_i，并生成相应的叶节点(X_i’，Q’-{B})}；

(7) 转(2)。

4 防雹预警和作业输出

当风暴识别成功后，根据风暴移动路径和预报结果，考虑风暴预报路径存在一定的偏差，在实际作业决策指挥过程中，暂定以风暴移动方向为中心轴，左右60°夹角所包含的扇形区域内所有防雹作业站点作为预警和作业相关联的站点(如图 2)，其中距离风暴强中心距离10 km的范围为需要立即作业的站点所属区域，即作业区，距离风暴强中心10~30 km的范围作为准备作业的站点所属区域，即预警区，利用GIS区域识别技术，即根据炮点的经纬度位置，判断炮点是否落在作业区或预警区内，并对落在该区域内的站点通过短信平台输出预警信息或作业参数。

4.1 风暴跟踪和预警

本文主要利用雷达对风暴进行实时跟踪，并对风暴的移动路径做出预报，从而对有冰雹作业条件的风暴、最佳作业站点及作业时机等作出实时预警分析和判断。为了能较准确地得到风暴移向、移速等信息，

需要对在两个相继体扫中识别的风暴单体进行时间相关处理以确定每一个被识别的风暴单体的路径。

(1) 检查体扫之间的时间间隔。如果该时间间隔大于阈值TIME(缺省值=20 min)，则将不进行体扫间的时间相关。

(2) 利用前一个体扫单体质心的位置及移动方向、移动速度确定目前体扫中该单体的初猜位置。如果该单体在前一个体扫中首次被识别，则利用前一个体扫所有风暴单体运动向量的平均值作为该单体的移动方向和速度。

$\begin{array}{l} X = {X_0} + V \times {\rm{cos}}\alpha \\ Y = {Y_0} + V \times {\rm{sin}}\alpha \end{array}$

(2)

式中，X、Y为单体的初猜位置，X₀、Y₀为单体的历史位置，V为单体的移动速度，α为单体的移动方向。

(3) 计算在当前体扫中识别的每个单体与前一体扫识别出的单体的相关性。首先，计算当前体扫中识别的每个单体与前一体扫的所有单体初猜位置的距离，如果距离小于一个阈值，则与初猜位置相对应的旧单体被储存下来作为一个与当前单体的可能匹配。然后对当前单体与其所有可能匹配的旧单体，按照强度分布进行相关分析，相关性最大的被认为是同一单体。如果是同一单体，就被赋予与前次体积扫相同的ID(单体识别号)。如果不是同一单体，就被赋予新的ID。

$r = \frac{{\sum ({Z_{i,j}} - \bar Z)(Z_{i,j}^\prime - {{\bar Z}^\prime })}}{{\sqrt {\sum {{({Z_{i,j}} - \bar Z)}^2}\sum {{(Z_{i,j}^\prime - {{\bar Z}^\prime })}^2}} }}$

(3)

式中，Z为当前单体各点的强度，Z′为可能匹配的旧单体各点强度。

(4) 计算所有至少已经在相继两个体扫中被识别的目前单体的新的运动向量。单体运动向量的计算主要根据该单体历史体扫位置，使用线性最小二乘拟合法确定单体预报时效内每一时刻对应的预测位置。

4.2 作业判别与指挥

以提示预警的风暴为中心，寻找周围的作业站点，根据风暴的类型(单体风暴、多单体风暴、强风暴)和作业工具的类型(火箭、高炮)，设计了六套防雹作业技术方案(李红斌等，2011)。系统自动依据风暴单体所在位置选择周边合适的火箭或高炮进行作业，并自动输出作业的方位角、仰角和用弹量等作业参数，同时根据风暴移动的方向和速度，对未来可能作业的作业区域进行提前预警，输出预警参数。其中图 3a为雷达实时监测结果，包括风暴识别结果，作业点位置和作业工具；图 3b为实时输出的作业参数，包括防雹预警信息、防雹作业信息和增雨作业信息。

5 算法检验 5.1 检验标准

从重庆、辽宁大连和河南三门峡等地的冰雹发生地来看，多数冰雹天气发生在山区，有冰雹观测的测站比较少，其实况记录主要来源于国家级测站的报告和人工防雹作业点上报的作业记录，因此本文做如下规定：(1) 实况检验范围：以有冰雹观测项目的站点和有作业记录的人工影响天气作业站点为参与检验的站点。(2) 识别单体检验范围：以参与检验的实况观测站点为中心，在其30 km×30 km的范围内选择识别出来的单体风暴、多单体风暴和强风暴作为检验对象。(3) 时间匹配窗口：由于冰雹发生往往滞后于雷达观测，将实况与识别结果检验定义在冰雹观测时间前的1~2个体扫，即6~12 min作为时间匹配的窗口。

通过对比被识别的风暴单体和冰雹实况出现位置和时间，分析冰雹识别检验效果。根据灾害性天气评分原则(俞小鼎等，2006)，采用命中率(POD)、虚警率(FAR)和临界成功指数(CSI)来计算识别评分。如果在时间匹配窗口范围内，两者对应，则认为识别正确；如果识别到的风暴单体没有对应冰雹出现，则认为虚警；如果冰雹出现的站点附近没有识别到风暴单体，则认为漏报。其中虚警率计算以作业参数输出的站点作为样本计数，漏报以实况有冰雹记录的站点作为样本计数。

5.2 风暴分类识别评分

统计2006—2014年重庆、辽宁大连和河南三门峡冰雹天气过程，共搜寻到有雷达观测和冰雹报告的典型冰雹过程31次(重庆18次，主要发生在东北和东南部地区；大连6次；三门峡7次)，累计出现冰雹182站次。表 2给出了累计识别出的单体风暴、多单体风暴和强风暴的次数，实况观测到的冰雹站数，及冰雹识别评分信息。

考虑到雷达连续观测，一个风暴体会多次被识别到，在计算评分时，在同一检验区域只计一次，所以实际参与评分的风暴体小于累计识别到的风暴体次数。从评分效果来看，单体风暴的命中率最低、虚警率最高，而强风暴的命中率最高、虚警率最低。分析其原因，单体风暴虚警率高的原因主要是单体风暴强度不强，出现冰雹的概率较小，而命中率较低的主要原因，一是当冰雹出现区域距离雷达位置较远时，雷达能探测到的层次较少，识别出来的风暴厚度等参数远小于实际风暴体的厚度；二是可能与0℃和-20℃层高度有关。强风暴命中率较高和虚警率较低的原因主要是强风暴的判别参数较为严格，当风暴被分类为强风暴时，对应出现冰雹的概率较高。这样的分类结果，为实际人工防雹作业时也提供了重要的参考依据，当出现多单体风暴及强风暴时应考虑人工防雹作业，而出现单体风暴时，可以根据实际情况选择性作业。

5.3 主要风暴参数的演变分析

为进一步了解风暴分类过程中风暴特征参数的演变情况，检验冰雹作业预警和输出的时机选择，本文选择2010年5月6日发生在重庆垫江—梁平的强风暴天气过程和2014年7月14日发生在河南三门峡市卢氏县的多单体风暴天气过程，从风暴顶高、风暴底高、风暴最大反射率因子、风暴累积液态水含量和风暴质量等主要风暴参数分析风暴的演变过程，结果表明风暴各项指标的演变均与风暴的发生、发展和消亡过程相对应，其预警或作业信息的输出时间均在风暴发展的中前期，通过作业能有效抑制冰雹云中冰晶的凝结，发挥防雹作业的作用。

5.3.1 “5.6”大风冰雹天气过程

2010年5月6日00：30，风暴分类算法开始对重庆境内的强风暴进行预警。从7号强风暴单体来看(表 3)，00:30，识别出来的7号风暴为两个，分别位于垫江县的高峰镇和沈家镇，分类结果为雷雨云；00:36，7号风暴发展并合并，位置移动到垫江县黄沙乡，强度发展为强风暴；00:42，7号风暴又开始分裂，一个位于垫江县曹回乡，维持原有强度，为强风暴，一个位于垫江县永平乡为雷雨云；01:00，曹回乡的风暴体开始向梁平县的周嘉镇移动，风暴强度仍为强风暴(图 4)；01:24，风暴移动到梁平县的回龙镇，强度有所减弱，变为雷雨云；01:30，风暴在原地加强为强风暴；01:36，风暴向东北方向移动，此后逐渐减弱。

从风暴的主要结构特征参数时间演变来看, 风暴最大反射率因子在01:06左右达到最大，然后逐渐递减，与风暴识别的结果基本一致(图 5a)；风暴垂直累积液态水含量在01:15左右达到最大，然后呈波动状态逐渐递减(图 5b)；风暴顶高在00:54左右到达最大，随后呈阶梯状递减，而风暴底高呈线性递增，风暴厚度逐渐减小(图 5c)；风暴质量在01:12左右到达最大，然后迅速降低，风暴呈减弱状态(图 5d)。

综上分析，风暴从00:30分开始处于发展阶段，01：00左右风暴处于成熟阶段，然后风暴开始逐渐衰减，从风暴分类结果来看，风暴基本上也是从雷雨云逐渐发展为强风暴，然后再变化为单体风暴，最后为雷雨云。这一过程与实况(在垫江沙坪、梁平回龙分别于6日01：12和01：23出现了瞬时最大风速达11级的大风，并伴有冰雹出现)基本一致，说明风暴在00:48—01：30处于最强阶段，而大风冰雹正好也发生于这个时段。从作业时机的选择来看，00：36开始预警信息和作业参数信息均有输出，比实际降雹时间提前了近40 min。

5.3.2 “7.14”冰雹天气过程

“7.14”冰雹天气过程为2014年7月14日河南省三门峡市雷达观测到的一次多单体风暴产生的冰雹天气过程。这次过程时间较短，从风暴发展到消亡，约30 min，回波底较高，在0.5°仰角的雷达PPI上基本无明显的回波出现。11：25—11：49风暴分类算法均在卢氏县境内识别到多单体风暴，系统向东南偏东方向移动(表 4、表 5)。11：25，留书岭和耿家岭炮点提示作业并输出作业参数；11：31，留书岭炮点提示作业，机场炮点出现预警；11：37留书岭和耿家岭炮点提示作业，涧北和机场两个炮站出现预警；11：49机场炮点提示作业。其中11：37为风暴发展最旺盛的时候(图 6)，出现两站作业和两站预警，而11：40左右在风暴单体偏南区域出现冰雹。从实况和作业输出情况来看，此次作业的位置和输出时间与实况较吻合，能有效对冰雹作业进行决策指挥。

从此次风暴的主要结构特征参数时间演变来看, 风暴最大反射率因子在11:31左右达到最大，然后逐渐递减，与风暴识别的结果基本一致(图 7a)；风暴垂直累积液态水含量在11:31左右达到最大，然后呈波动状态逐渐递减(图 7b)；11:25—11：37风暴顶高较高，风暴低高较低，风暴较为深厚，随后风暴减弱(图 7c)；风暴质量在11:31左右到达最大，然后迅速降低，风暴呈减弱状态(图 7d)。

6 结论与讨论

利用雷达、常规探空和冰雹灾害报告资料，基于风暴识别算法，对风暴结构特征参数建立决策树模型进行分类，并对不同强度风暴体按防雹作业敏感区域进行自动预警或作业输出。通过对2006—2014年期间三地发生的31次冰雹天气过程识别检验，得出了以下结论：

(1) 通过改进风暴识别算法，增加0℃层高度和-20℃层高度对应的风暴中心强度、风暴反射率因子垂直递减率、回波强度30和45 dBz所在高度、中高层闭合强中心所在高度等结构特征指数，能有效地提高风暴的识别准确度。

(2) 建立基于决策树模型的分类识别算法，能有效地将雷达识别出的风暴体按强度进行分类，并按不同强度的风暴设计相应的作业方案，有效地输出防雹作业敏感区所在的位置和用弹量等作业参数，提高人工防雹决策的自动化程度和效率。

(3) 在人工防雹作业中，根据风暴位置及移动方向，将下游30 km范围划定为作业区和预警区，有效提高了作业预警时间，作业时机选择也基本上在风暴发展的中前期，能提前通过作业抑制冰雹云的生长。

综上所述，基于决策树模型的风暴分类识别算法自动化程度高，能提高人工防雹作业的效率。但该方法仅应用了单多普勒雷达资料，受限于雷达探测的局限性，无法完全识别远距离的风暴特征，下一步将利用雷达拼图资料予以完善。