一种基于地面实况的降雹风暴体客观标识方法

刘伯骏 张亚萍* 黎中菊 韩 潇 芦 华 张 勇

1)(重庆市气象台，重庆 401147) 2)(重庆市气象科学研究所，重庆 401147)

引 言

近年随着计算机科学的发展和大数据处理能力的提升，应用深度学习等人工智能技术进行强对流风暴识别和预报，成为越来越多研究者的选择[1-3]。深度学习一般分为监督学习和无监督学习，当带标识的样本较少时，也会采用半监督学习。无论监督学习还是半监督学习，标识工作是建立深度学习数据集的关键基础，这在观测数据匮乏的灾害性天气智能预报中尤为重要。标识正确、质量可靠的数据集，对深度学习算法准确率有决定性影响。由于强对流事件发生概率低、过程时间短、空间尺度小，导致观测难度大，观测数据往往不完备，制约了深度学习在相关领域的应用。其中，冰雹的地面实况观测最少，观测时间不明也很常见。将冰雹实况信息与天气雷达监测和跟踪的风暴体进行匹配，可以充分利用地面实况，提高冰雹观测数据完备性。

本文将无地面实况情况下，根据天气雷达观测到的风暴体特征进行降雹风暴体(下文简称雹暴)判识称为雹暴识别，而将标明导致地面出现冰雹实况的风暴体称为雹暴标识，即将天气雷达观测到的强对流风暴体与地面实况相匹配，算法核心在于根据地面冰雹事件发生的时空信息，找出导致该事件的风暴体(雹暴)。为弥补地面观测不足，Gange等[2]将雹暴识别结果(WSR-88D中冰雹探测算法[4])用于建立冰雹实况数据集，并指出气候统计证实冰雹识别结果与冰雹地面报告的空间分布基本一致[5]。张秉祥等[6]根据地面站的冰雹观测，利用其相对准确的时间，搜索冰雹发生时刻及前两个体扫的风暴体，通过综合指数标识雹暴。王瑾等[7]在雹暴标识上，采用以下规则：①风暴单体识别和跟踪(the storm cell identification and tracking，简称SCIT[8-9])算法识别的风暴体要在人工影响天气炮点5 km范围内；②风暴体在连续两个体扫时间内达到一定强度；③对于有准确时间记录的冰雹事件，识别出风暴体的时间必须在地面记录时刻前15 min到记录时刻后5 min范围内；④对于无准确时间的降雹记录，考虑距离炮点最近的风暴体为雹暴；⑤风暴体必须在雷达静锥区外且在观测范围内。若满足上述条件的风暴体不唯一，以距离为判断依据。然而炮点观测有限，对于如何利用更普遍的、时空信息更模糊的地面灾情报告进行雹暴标识，仍没有完整的解决方案。

在以上雹暴识别和标识方法的基础上，本文利用重庆地区灾情信息中有准确观测时间的冰雹过程实况数据，选取空间上有利于描述地面观测与风暴体匹配程度的判别因子，采用模糊逻辑算法，提出一种融合地面冰雹实况信息和雷达观测的客观标识方法，以解决时间信息模糊时的雹暴标识问题。

1 数 据

1.1 雷达数据及产品

利用中国气象局灾害性天气短时临近预报预警业务系统SWAN(Severe Weather Automatic Nowcast System)进行雷达回波三维拼图及SCIT产品计算。由于重庆及其周边雷达站陆续建成，不同年份拼图所使用的原始雷达数量不同：2008—2009年仅使用重庆雷达数据，2010—2011年使用重庆、万州2部雷达数据拼图，2012年使用重庆、万州、黔江3部雷达拼图，2013—2014年使用重庆、万州、黔江、永川4部雷达拼图，2015—2019年使用重庆、万州、黔江、永川、南充、宜宾、成都、达州、遵义、恩施、宜昌、怀化共12部雷达拼图。SCIT算法使用雷达拼图数据的水平分辨率为1 km×1 km，时间分辨率为6 min。

1.2 灾情报告

灾情报告是一种描述气象灾害详细情况的信息。灾情信息通常由气象信息员观测收集，或通过民政部、应急管理部等相关单位共享得到，经各基层气象部门汇总后形成灾情报告，上报到中国气象局气象灾害管理系统，供相关用户下载使用。气象信息员行政村(社区)覆盖率一直稳定在较高水平，这一数据被认为有较好一致性和较高空间密度[10]。但应注意到灾情报告通常为事后收集，对冰雹等灾害性天气的时空描述往往比较模糊，仅有少量过程能提供相对精确的降雹时间点。

本文将冰雹事件分为两类：有相对准确时间信息的过程(简称时间确定过程)和无相对准确时间信息的过程(简称时间不确定过程)。重庆地区灾情报告显示：2008—2019年共发生冰雹事件64次，27个区县受冰雹灾害影响。其中时间确定过程19次，时间不确定过程45次，但受到观测数据限制(如雷达数据缺失或雹暴处于雷达观测盲区)，可用于标识算法研究的时间确定过程和时间不确定过程分别仅有13次和22次(时间确定过程见表1，其中时间为北京时，下同)。为此，基于时间确定过程中2008—2018年的8次过程(简称参照集)构建标识算法，将2019年5次过程(简称验证集)用以检验。

表1 2008—2019年重庆地区有相对准确时间和雷达观测数据的冰雹过程Table 1 Hail cases with exactly time and radar observation in Chongqing during 2008-2019

2 雹暴客观标识方法

本文基于2008—2019年重庆地区冰雹灾情报告，确定冰雹发生的大致时空信息。根据对应时段的SCIT产品确定备选风暴体，并基于三维拼图数据计算风暴体雷达观测特性信息。利用模糊逻辑算法，综合时空信息和风暴体雷达观测特性进行雹暴标识(图1)。

图1 雹暴客观标识方法流程Fig.1 Diagram of the object-labeling algorithm for hailstorm

2.1 指标与隶属函数

由于强对流天气的发生发展需要满足多个物理条件，且不同类型强对流天气的物理量统计结果表明，很难找到一个完全明确的、单一物理量阈值表征某类天气发生发展的物理条件[11-12]，因此需要能够综合应用代表不同物理条件的多个物理量的方法，模糊逻辑方法正是其中之一[13]。1965年美国数学家Zadeh[14]首先提出模糊集合的概念，标志着模糊数学的诞生。模糊逻辑善于表达界线不清晰的知识与经验，借助隶属函数的概念，处理模糊关系，在气象领域应用广泛[6，15-21]。本文根据灾情报告给出的时空信息，将距离也引入隶属函数，结合以往研究雹暴识别的常用指标，采用模糊逻辑算法，形成雹暴客观标识方法。

结合前人研究成果[22-24]与重庆本地经验总结[25]，选取风暴体最大反射率因子(CRM，单位：dBZ)、风暴体最大回波顶高(ETM，单位：km)、最大垂直积分液态水含量(VILM，单位:kg·m-2)及45 dBZ所在最大高度(H，单位：km)作为雷达观测指标，选用风暴体质心与冰雹灾害初始猜测位置(简称初猜位置，精确到秒)间的距离(DIST，单位：km)作为地面观测指标。

2.1.1 风暴体质心与初猜位置间的距离(DIST)

初猜位置由灾情报告中冰雹灾害的空间位置信息确定，进而根据SCIT算法识别的风暴质心位置计算DIST。Witt等[4]认为，当最近的风暴体与地面灾情报告位置相距50 km或更远时不合理，因此DIST的上限被设置为50 km；作为冰雹母体的雷暴是γ中尺度系统，根据Orlanski[26]尺度划分理论，其水平特征空间尺度为20 km。若假设风暴体近似为圆形，则其半径为10 km。因此，当DIST不大于10 km时，可认为在距离上已达最优情况(图2)。

2.1.2 风暴最大反射率因子(CRM)

多普勒天气雷达是监测对流系统的基本手段，强对流天气系统常伴有强烈的上升运动和显著的水汽相态变化[27]，CRM较大是表征强对流的重要指标。本文设置CRM上限时参考张秉祥等[6]在利用模糊逻辑进行雹暴判识时的阈值60 dBZ；考虑到重庆地形复杂、遮挡严重，将CRM隶属函数下限设置为30 dBZ(图2)。

2.1.3 风暴最大回波顶高(ETM)

ETM指风暴体投影范围内，18 dBZ所能达到的最大高度。相比于雷暴和短时大风等强对流天气，冰雹对应的风暴体回波顶更高[28]。张秉祥等[6]将回波顶上下限分别设置为8 km和6 km。本研究中，为了将更多风暴体纳入备选范围，将ETM下限降至5 km。结合重庆本地强对流统计结果，强对流天气的ETM主要在10 km以上，以此确定ETM上限为10 km(图2)。

2.1.4 风暴45 dBZ所在最大高度(H)

H指风暴体中45 dBZ回波所到达的最大高度。Smith等[29]指出，云内最初冰雹凝结体直径中位数约为0.4～0.5 cm。根据Marshall-Palmer雨滴谱指数分布关系式推导出，冰雹凝结体初期等效雷达反射率因子约为45 dBZ。判断冰雹简单有效的方法是比较强回波高度相对于0℃和-20℃等温线高度的位置[6]，当45 dBZ超过0℃高度时，可能出现冰雹；当45 dBZ回波所在高度超过-20℃层时，出现冰雹概率极大。根据2014—2018年重庆沙坪坝站探空数据统计，3—9月重庆0℃层和-20℃层平均高度如表2所示。据此，将H的下限和上限分别设置为3 km和8 km(图2)。

表2 2014—2018年3—9月重庆沙坪坝站平均0℃和-20℃层高度Table 2 Average heights of 0℃ and -20℃ at Shapingba sounding station from Mar to Sep during 2014-2018

图2 降雹风暴体客观标识方法中各指标的隶属函数Fig.2 Membership functions of each index of object-labeling algorithm for hailstorm

2.1.5 风暴最大垂直积分液态水含量(VILM)

VILM定义为风暴体投影范围内最大垂直积分液态水含量，值越大，对流性越强。该指标对于冰雹过程有较好指示意义。张秉祥等[6]将VILM的下限和上限分别设置为35 kg·m-2和55 kg·m-2。考虑重庆地区雷达波束有时受地形遮挡影响，本文对上述指标略作下调，分别为25 kg·m-2和50 kg·m-2(图2)。

2.2 权重设置

本文对DIST，CRM，ETM，H，VILM的权重分别设置为0.4，0.2，0.2，0.1和0.1。其中，对CRM，ETM，H，VILM，本文参考已有研究成果设置权重。而DIST作为标识算法中独有的指标，将测试不同距离权重对标识算法结果的影响。为此，设定以下敏感性试验：当距离权重取值自0.1逐步递增至0.9时，其余物理量按比例分配剩余权重，保持各权重之和为1。将以上描述可表示为

wDIST=i/10，i=1，2，3，…，9；

(1)

(2)

(3)

式(1)～式(3)中，wDIST，wCRM，wETM，wH，wVILM分别表示DIST，CRM，ETM，H，VILM的权重。

当距离选取不同权重时，算法在参照集中标识结果变化如图3所示。图中横坐标为距离权重，纵坐标为不同过程中，算法标识时间与地面灾情报告中记录的冰雹发生时间之差，其正(负)值代表算法标识时间较晚(早)。当wDIST较小，算法挑选发展更旺盛、距离稍远的对流风暴体；反之，当wDIST过大，则标识更多参考距离的影响。由图3可见，对于部分过程(如2008年4月10日、2011年7月23日、2014年4月17日、2015年7月25日的冰雹过程)，客观标识算法结果与距离权重无关，这是因为雹暴最强时刻与地面观测到的降雹时刻相对应。但对于其余过程，降雹时刻与雹暴最强时刻不同，此时引入距离信息可以改善标识效果(如2008年6月5日、2010年5月6日、2018年3月17日的冰雹过程)。总体上，随着DIST权重增加至0.4，时间差的方差(图略)收敛至最小，平均值接近于0。因此，后续计算中，取DIST权重为0.4。

图3 距离权重选取对标识结果的影响Fig.3 The influence of distance weight on labeling results

3 标识算法结果

3.1 总体效果

由于雷达拼图时间以6 min为间隔，标识结果很难与地面灾情报告时间一一对应。若标识结果与观测时间之差处于[-18 min,12 min]的时间窗内，则认为标识正确。总体而言，客观算法可以准确给出冰雹发生时间及对应的风暴体，针对8次过程，有7次能够正确标识，其中有5次标识时间与灾情报告描述的时间相差在1个体扫(6 min)内(表3)。对于2008年6月5日丰都小冰雹过程，因SWAN中的SCIT算法未识别出风暴体，所以造成客观算法未能标识成功。表3同时给出客观算法对不同过程所标识出的风暴体的匹配程度，即模糊逻辑中的综合隶属度。匹配程度由各判别因子乘各自对应权重并求和得到。表3表明：参照集中风暴体匹配程度普遍较大(大于0.78)，平均匹配程度为0.924，标准差为0.096。匹配程度与冰雹尺寸的相关性较差，如对大冰雹而言，既有匹配程度达到1.000的2014年沙坪坝冰雹，也有匹配程度为0.805 的2018年彭水冰雹。分析发现，绝大部分过程中导致匹配程度下降的原因是垂直积分液态水含量较小，而垂直积分液态水含量主要与冰雹发生地点与雷达站的距离有关。

表3 雹暴标识客观算法在参照集的标识效果Table 3 Results of objective labeling algorithm for hailstorm based on reference set

图4展示冰雹标识算法给出的标识时刻垂直积分液态水含量、组合反射率因子以及SCIT算法给出的风暴路径信息。为了更好地测试算法特性，从以下几个方面进行试验：①针对上述过程，输入数据时间长度随机分布，使每个过程计算时间长度为4.5～8 h。结果显示，输入数据时长对标识结果没有显著影响。②上述过程共有大冰雹过程3次，小冰雹过程3次，冰雹大小不明过程2次。若SWAN中SCIT算法可以识别出风暴体，则冰雹尺寸对标识结果无显著影响。③冰雹灾害所发生的月份不同、时间不同(午后或夜间)，未对算法的准确性产生影响。④算法对风暴体生命史长度不敏感。无论对2010年5月6日的长生命史风暴体，还是2015年7月25日丰都太平坝乡的局地生消雹暴，标识算法均能给出较好标识结果。

图4 各冰雹过程中算法标识时刻的垂直积分液态水含量(填色)、组合反射率因子(实线,单位：dBZ)及SCIT识别的风暴体路径(三角，红色为标识时刻)Fig.4 Vertical integrated liquid water(the shaded),composite reflectivity(the line，unit:dBZ) and storm path from SCIT at labeling time for each hail case(triangles，the red one means labeling time)

续图4

3.2 标识失败过程分析

虽然客观标识算法整体标识效果较好，但在2008年6月5日南川冰雹过程中标识效果并不理想：地面灾情报告中提及的时间为15:20，而标识结果为15:06，时间偏差为14 min；对于丰都冰雹，灾情报告中记录时间约为15:30，而标识算法未能给出结果。

图5给出此次过程降雹临近时刻SWAN拼图组合反射率因子及SCIT识别风暴结果，且标明地面观测指出的降雹位置和时间，15:24时刻数据缺失。由图5可见，对于丰都冰雹，SCIT算法在临近时刻均未识别出风暴体，算法标识失败。对于南川冰雹，SCIT算法在15:18，15:30也未能识别出风暴体，因此标识算法仅能从15:06，15:12，15:42中挑选。风暴强度15:42最强，15:06次之，15:12最弱。结合距离因素，客观算法最终将冰雹发生时间定为15:06。以上分析表明：客观标识算法对风暴体识别结果较为敏感。

图5 2008年6月5日南川、丰都冰雹过程组合反射率因子演变及风暴路径信息Fig.5 Composite reflectivity evolution with storm path information in the hail case at Nanchuan and Fengdu on 5 Jun 2008

3.3 初猜位置对标识结果的影响

由标识方法构建过程及权重分配可知，在标识过程中，初猜位置对标识结果影响较大。以2008年4月10日南川冰雹为例分析初猜位置偏差的影响，并利用有SCIT识别结果的参照集过程，分析初猜扰动对标识结果的影响。

对2008年4月10日南川冰雹过程，客观算法标识时刻为21:42，而灾情报告指出南川德隆镇21：30 开始降雹，时间偏差为12 min。图6为21:30—21:42组合反射率因子演变及SCIT识别单体情况。由于临近降雹时刻，仅有1个单体在德隆镇上空，因此可确定该单体造成降雹。追溯历史位置，该风暴体以约40 km·h-1的速度向东移动，并于21:30—21:36分进入德隆镇边界范围，与灾情报告时间基本对应。标识偏差相对较大与初猜位置定位有关。初猜位置由灾情报告中乡镇名称定位而得，因此定位位置常为乡镇政府所在地或人口密集区。对于面积较大(或行政边界中，任一方向距离较远)的乡镇，这种定位方式会导致误差，对标识结果产生不利影响。由此可以推断，若降雹风暴体尺度较大或呈带状，可能导致风暴识别算法对风暴定位与降雹位置产生一定偏移，影响标识结果。

图6 2008年4月10日南川冰雹组合反射率因子演变及风暴路径信息Fig.6 Composite reflectivity evolution with storm path information in the hail case at Nanchuan on 10 Apr 2008

为了定量探讨初猜位置对标识结果影响，通过改变参照集(除去无SCIT产品的2008年6月5日丰都冰雹)的初猜位置进行扰动试验：以原始初猜位置为圆心，10 km为半径，每隔30°方位构建新的初猜位置，并基于新位置进行冰雹标识。标识结果平均时间偏差和均方根误差(图7)显示，初猜位置的10 km偏差会对标识结果产生平均值为[-8 min,12 min]的影响，这在多数情况下不影响标识准确率。绝大部分过程中，时间偏差的均方根误差小于14 min，与平均偏差量级相当。但对于部分过程(如2008年6月5日南川冰雹)，初猜位置的偏移会引起较大标识结果离散度，这可能与当天回波较为分散、备选风暴体较多有关。

图7 初猜位置扰动对标识结果的影响Fig.7 Influence of initial guess position disturbance on labeling results

4 标识算法结果检验

为了检验客观标识方法的适用性，本文分别对验证集和时间不确定过程进行标识。对于验证集，在标识过程中去除准确时间信息，使算法在标识过程中将时间搜索范围扩大到冰雹发生前后的几个小时，从而模拟灾情报告没有准确时间信息描述的情形。如果算法标识时间处于地面观测时段或其邻近时间窗范围内，则认为算法标识正确。对于时间不确定过程，文中将标识算法给出的时间与经验丰富的短临预报员给出的标识结果进行比较，以判断客观方法是否与预报员标识结果接近。

4.1 验证集检验结果

对于已建立的客观标识方法，采用验证集进行测试，结果显示：对于验证集5个过程，客观算法标识正确率为100%，匹配程度分布在0.887至1.000之间(表4)，这与标识算法基于参照集的分布基本相当。对各过程中匹配程度丢分项分析表明：除2019年4月24日巫溪冰雹过程因风暴体质心距离初猜位置较远导致匹配程度低之外，其余过程均因风暴体垂直积分液态水含量较低而导致匹配程度降低。以上结果初步证实客观标识算法具有较好的泛化能力，可以在时间信息模糊的情况下进行正确标识，为冰雹数据集的扩展增加可能性。

表4 雹暴标识客观算法在验证集的标识效果Table 4 Performance of objective labeling algorithm for hailstorm based on validation set

4.2 时间不确定过程标识结果检验

时间不确定过程中，灾情报告无法提供精确的时空信息，但可以将标识算法得到的时间与经验丰富的短临预报员给出的标识结果进行比较，判断二者是否一致。针对22次时间不确定过程，预报员给出的标识结果与客观算法基本一致：18个过程标识正确；2个过程回波偏弱，因预报员未给出标识结果无法比较；2个过程时间相差较远，分别为30 min和36 min，经查验，由于灾情报告给出的空间范围较大，主客观标识过程中对初猜位置的认定有出入，导致标识时间相差较远。此外，比较主客观标识结果显示：预报员更倾向于将标识时间提前1～2个体扫，标识出的风暴体较灾情报告描述地点稍远、强度稍弱，在此情况下，主客观往往标识的是同一个风暴体，因此对数据集构建影响不大。

以2018年4月30日江津中山镇冰雹过程为例，考察标识算法在时间不确定过程中的标识效果。针对此次过程，灾情报告给出如下描述：“在低层低涡和冷空气的共同作用下，4月30日20:00—5月1日14:00我区各地出现明显雷阵雨天气，部分地区伴有雷电、短时强降水、阵性大风、冰雹等强对流天气。具体雨量……受局部大风和冰雹影响，我区中山镇遭受风雹灾害。”据此，测算初猜位置为28.79°N,106.28°E，并确定搜索时间段为4月30日20:00—5月1日14:00。将以上信息输入到标识算法，得到标识结果：4月30日20:12中山镇西部风暴体导致降雹。图8a给出4月30日20:12 SWAN组合反射率因子、初猜位置信息。风暴体最强回波超过65 dBZ，且距初猜位置较近。对风暴体沿图8a中AB两点进行剖面(图8b)，可见最强回波中心约处于5 km高度，超过当日0℃层所在高度(30日20:00邻近探空站测得0℃层高度为4.38 km)，且45 dBZ延伸至12 km左右(-20℃层高度为7.63 km)。风暴体东北侧低层有弱回波区，对应低层较强的上升气流。以上分析可知，识别算法给出降雹时间的雷达回波符合传统降雹风暴体特征，该风暴体降雹可能性很大，且风暴体距离灾情报告中致灾位置较近，标识结果合理。此次过程，人工标识降雹时间为20:06，比客观方法提前6 min。此时刻风暴体已从西面接近中山镇，但距离稍远，强度基本相当。

图8 2018年4月30日20：12雷达组合反射率因子(a)及其沿图8a中AB所作剖面(b)Fig.8 Composite reflectivity at 2012 BT 30 Apr 2018(a) and profile along AB in Fig.8a(b)

5 结论与讨论

本文选取2008—2019年重庆地区冰雹天气过程，并将灾情报告与雷达观测的强对流风暴体进行合理匹配，发展基于地面实况观测的雹暴标识算法。对于参照集，8次冰雹过程，算法正确标识7次，其中有5次标识时间与灾情报告描述的时间相差1个体扫。通过验证集和时间不确定过程的检验，证实标识方法可在时间信息模糊的情况下进行标识。具体分析表明：

1) 冰雹大小对标识结果没有显著影响。冰雹灾害发生的月份和时间(午后或夜间)，未对标识算法的准确性产生影响。冰雹生命史长短对标识结果无显著影响，无论2010年5月6日的长生命史雹暴，还是2015年7月25日丰都太平坝乡的局地生消雹暴，算法均能给出较好标识效果。

2) 标识结果对SCIT风暴体识别算法较为敏感，对风暴识别算法能否识别出风暴体对标识结果有直接影响，其给出的风暴体质心位置(或风暴体大小等特性)对标识结果有间接影响。另外，初猜位置也会影响标识结果。

3) 验证集检验和时间不确定过程主客观标识结果表明，标识算法可在无准确时间信息情况下对雹暴进行合理标识。

以上研究所选用降雹过程有限，且均为重庆地区冰雹过程，后续研究中将丰富降雹过程，进一步验证标识算法的泛化能力。此外，在冰雹数据集中负样本(即该风暴体确定没有降雹)的标识方法有待研究。