1. 引言
龙卷是地球大气中一种具有局地性、小尺度、突发性特征的强对流天气,是一种最为剧烈且最具有破坏性的天气现象,世界范围内每年发生龙卷约为2000个,美国超过1200多个 [1] ,中国的龙卷发生频率约为美国十分之一 [2] ,由于龙卷小尺度的对流系统瞬变信息复杂,预报预警难度大,往往造成极大的人员伤亡和经济财产损失 [3] ,因此,进行龙卷的有效识别是一个亟待解决的问题。
20世纪70年代,美国强风暴实验室第一次观测到龙卷的涡旋特征(TVS),并进一步发现相距900 m的两个距离库之间存在较大切变以及对应着区域内的最大正负速度值 [4] 。之后通过进一步实验及模拟验证了TVS验证龙卷发生的可行性 [5] 。20世纪80年代,使用MDA算法大大提升了涡旋识别 [6] ,针对TVS阈值识别低探测率的问题,20世纪90年代,美国开发出一种新的龙卷识别方法(tornado detection algorithm, TDA),由此,龙卷识别范围由超级单体龙卷扩大至非超级单体龙卷 [7] 。进入21世纪,随着单偏振雷达的升级,双偏振雷达识别龙卷也有了长足的进展。Ryzhkov等发现在钩状回波区域,龙卷碎片特征对应的相关系数CC (低于0.8)及差分反射率ZDR (小于0.5 dB)都很低 [8] ,Umeyama等根据龙卷中碎片和水凝物的速度的差异性,依靠双极化谱密度,识别出碎片和水凝物,最后将碎片产生的速度偏差滤除 [9] 。Suzuki [10] 利用KSR雷达观测的数据,发现超级单体龙卷拥有ZDR弧,ZDR柱以及KDP柱的特征。在国内,戴建华等发现超级单体回波的有界弱回波区、回波悬垂、回波墙以及中气旋等典型特征 [11] 。郑永光等(2020)分析了2019年7月3日辽宁开原EF4级强龙卷的形成与消亡机制 [12] 。管理等(2022)基于经典概念模型和机器学习方法,开发了KDP足和ZDR弧的自动识别算法 [13] 。
上述工作表明,针对龙卷的识别,识别出中气旋,钩状回波、龙卷碎片、ZDR弧、ZDR柱、KDP柱等典型特征非常重要。本文提出在中气旋识别基础上,利用反射率因子Z和谱宽W及速度V之间的相关性以及在体扫的低仰角区域依托具有垂直相关性的ZDR弧、KDP足的特征,进而识别出龙卷的算法。针对2016年6月23日阜宁以及2016年5月9日广州三水地区的灾害性天气过程,进行龙卷识别结果的验证。
2. 资料选取
本文选取2016年6月23日阜宁地区以及2016年5月9日广州三水地区的灾害性天气过程,雷达资料分别是盐城S波段单偏振雷达以及广州S波段双偏振雷达,分辨率为0.25 km,0.25 km,时间分辨率皆为6 min。
3. 识别方法
针对单偏振雷达龙卷识别,Basedon Burgessand Magsig’s指出,龙卷发生区域内,反射率因子Z和谱宽W及速度V之间存在负相关性 [14] ,但在实际中,单偏振雷达准确识别龙卷发生区域很难,故扩大识别范围,在龙卷发生区域周围,往往反射率因子Z和谱宽W及速度V之间存在正相关性,利用这一特性,叠加中气旋识别算法,识别出龙卷。针对双偏振雷达,识别龙卷分为三步。第一步进行中气旋的识别。第二步进行强ZH单体,强KDP单体,强ZDR单体的识别,第三步进行质心计算,确认龙卷发生区域。
3.1. 指标相关性(ZH, V, W)
龙卷发生区域附近,具有较大的反射率因子,速度以及谱宽。且反射率因子和谱宽W以及速度V之间存在正相关性。具体方法如下:
选取5*5的模板,在该模板范围内,速度V平均值需满足−5 m/s至5 m/s之间,速度V绝对值均值需大于15 m/s,反射率因子Z均值需大于30 dBZ,谱宽W均值需大于3 m/s。ZW,ZV的相关性公式如下:
(1)
(2)
选取反射率因子,速度,速度绝对值满足条件的以及相关性PZW,PZV均大于0的区域,作为龙卷待选择地,之后进行中气旋的识别。
3.2. 中气旋识别
龙卷发生时,龙卷形成的涡旋接近几公里甚至十几公里,回波图上会有明显正负速度对,即TVS (龙卷涡旋特征)。具体算法如下:
首先进行一维旋转风矢量的提取,对低层仰角同一距离圈内的速度V进行筛选,保留速度连续增长的部分,储存元素分别有初始方位角Φb,结束时的方位角Φe,起始径向速度Vb,结束径向速度Ve,以及距雷达的距离r。之后进行三次测试,前两次须满足的公式如下:
(3)
(4)
第一次测试筛选出与切变有关的区域,所计算的角动量必须大于最低多普勒角动量阈值(50 ms−1∙km),第二次测试筛选出大于背景切变值2 × 10−3 s−1的旋转风矢量,第三次测试分为两步,首先将第一次测试计算得到的旋转风矢量的角动量与高角动量阈值Hm (缺省值150 ms−1∙km)比较,若大于则保留;若小于则将其与Hs (缺省值4 ms−1∙km−1)比较,若大于即可保存该旋转风矢量,进行下一步处理。
对保留下来的一维旋转风矢量,进行二维特征的合并,计算每个旋转风矢量的质心,判断相邻风矢量的质心距离Ir (缺省值1 km)是否在1 km之内,同时判断相邻方位角之差是否在Iaz (缺省值2.2˚)之内,若两者都满足,则合并成一个二维特征,并不断重复该过程,得到多个风矢量组成的二维特征。对每组二维特征进行筛选,若二维特征存在大于等于4个的一维风矢量,则将其保留。将保留下来的二维特征所对应的反射率因子进行反射率因子阈值筛选,剔除小于等于15 dBZ的数据,减少噪声及旁瓣回波的干扰。
最后计算二维特征的切变,计算公式如下:
(5)
若切变大于15 ms−1∙km−1,则判定为潜在TVS,再做垂直相关性分析,若有垂直相关特征,则将其判定为龙卷风。若切变阈值不能满足,则进行下一步单体联合识别。
3.3. 单体联合识别
单体联合识别涉及雷达反射率因子ZH,差传播相位常数KDP和差分反射率ZDR。首先对三种数据进行线性插值处理。最低层仰角的三种数据全部插值到300 × 300 (km)规则网格上,网格精度为0.25 × 0.25 (km),之后使用DBSCAN聚类算法进行聚类。
基于DBSCAN聚类算法的单体识别
在灾害性天气环境下,ZH,KDP等强单体数量未知,目标单体个数难以确定,相较于需要确定个数的K-Means聚类算法,DBSCAN算法是一种基于密度的聚类算法,无需确定个数,只需要输入搜索半径eps以及最少可达目标数minPoints。其次,算法对于离群点的识别也能减小杂波的影响。
DBSCAN算法首先在数据集中任选一个初始点划分到集群1,然后统计该点搜索半径eps范围内点的个数,若大于最少可达目标数minPoints,就保留,并将所有邻域内的点分配到集群1。然后,针对该集群内的其他点,重复搜寻该点搜索半径eps内点的操作并统计个数,若个数大于最少可达目标数minPoints,则将这些点拓展到集群1中,以此类推,直到集群1无法扩充为止。之后在集群1外围选取任意一点放入集群2,重复前述操作,直到数据集中的所有点都被分配到某个集群或者是杂波中。
针对反射率因子ZH,以平均强度超过45 dBz和单体面积超过20 km2两个约束条件进行单体的识别,即DBSCAN聚类算法的参数eps(半径)设置为0.5 km,minPoints (最少个数)设置为320。针对差传播相位常数KDP,在对应ZH大于35 dBz的部分,以大于1.5˚/km、单体面积大于2 km2作为约束条件识别强单体,即DBSCAN聚类算法的参数eps (半径)设置为0.5 km,minPoints (最少个数)设置为32,针对ZDR,使用DBSCAN聚类算法进行聚类,针对大于3 dB的部分,以面积大于1km2作为约束条件识别强单体,即DBSCAN聚类算法的参数eps (半径)设置为0.5 km,minPoints (最少个数)设置为16。为避免近地面层地物回波的干扰,还需引入ZDR单体平均CC大于0.88的补充约束条件。若强ZH单体与强ZDR单体的质心距离小于15 km,则识别强单体为ZDR弧。
对于处理完毕的KDP足和ZDR弧,对其相邻仰角数据进行同样方法处理,确保KDP足以及ZDR弧存在垂直结构,之后计算识别出的中气旋,KDP足和ZDR弧的质心距离,若质心间距小于5公里,则将质心连线的中心处识别为龙卷发生地。
4. 案例分析
4.1. 阜宁龙卷–盐城S波段单偏振雷达
2016年6月23日,江苏省盐城市阜宁县出现了一次强对流天气,在涡旋系统的影响下,江苏省阜宁市出现了雷暴,冰雹,大风等灾害性天气。本文采用盐城SA波段单偏振雷达观测数据,选取时间范围在2016年6月23日5时12分至2016年6月23日7时0分内的雷达资料,该部分雷达数据距离分辨率为0.25 km,时间分辨率为6 min。图1中,图1(a)是0.5˚仰角上6:26:10的反射率因子图,强雷暴中心位于雷达中心西北部55 km处,有明显的钩状回波,钩状回波凹形中心的最大反射率因子超过60 dBZ,沿着钩状边缘线凹处有明显的弱回波区域。该处对应的同时段的速度如图1(b)所示,该处有明显的涡旋结构,从速度色值上看出,正负速度对的风速相差巨大,判断该处有剧烈的大风灾害。
Figure 1. (a) Radar reflectivity factor (ZH) plot of SA-band single-polarization radar in Yancheng on June 23, 2016 at 06:26; (b) Velocity (V) plot of SA-band single-polarization radar in Yancheng on June 23, 2016 at 06:26
图1. (a) 盐城SA波段单偏振雷达2016年6月23日6时26分的雷达反射率因子(ZH)图;(b) 盐城SA波段单偏振雷达2016年6月23日6时26分的速度(V)图
对此次数据进行中气旋的识别,图2是识别出的中气旋结果图:该处的风切变高达21.36 ms−1∙km−1,大于15 ms−1∙km−1,识别为TVS,并对其进行垂直相关性分析,1.45˚仰角上识别出的TVS质心位置与仰角0.5˚上识别出的TVS质心位置距离相差3.57 km,具有垂直特征结构,该处识别为龙卷发生地。
Figure 2. Mid-cyclone identification map of SA-band single-polarization radar in Yancheng at 06:26 on June 23, 2016
图2. 盐城SA波段单偏振雷达2016年6月23日6时26分的中气旋识别图
对连续时间段的数据进行龙卷的中气旋识别,最终识别出在6:26分和6:48分出现具有垂直结构的中气旋,结合ZV,ZW相关性识别出来的区域,识别结果一致,将其判断为龙卷,图3为6:26分用指标相关性识别出来的区域,表1显示了两种方法在不同时间段的识别情况。
Table 1. Identification judgments of different time periods using indicator correlation and mesocyclone identification
表1. 利用指标相关性和中气旋识别对不同时间段的识别判断
Figure 3. Identified regions using indicator correlation
图3. 利用指标相关性识别出的区域
4.2. 三水龙卷–广东佛山S波段双偏振雷达
2016年5月9日,广州三水范围内出现了一次强对流天气,本文采用广州S波段双偏振雷达观测数据,选取时间范围在2016.5.9 7:30~2019.5.9 10:48内的雷达资料,该部分雷达数据距离分辨率为0.25 km,时间分辨率为6 min。图4中,图4(a)是0.51˚仰角上9:00的反射率因子图,强雷暴中心位于雷达站西北部75 km处,形状为钩状回波,中心区域的反射率因子大于60 dBZ,周围强回波向东北延伸,呈带状分布,图4(b)是0.51˚仰角上9:00的速度图,径向速度反映出偏西风气流,西北部75 km处以及北部100 km处具有一些蜗旋特征,系统有发展潜力。图4(c)是0.51˚仰角上9:00的差分反射率ZDR图,图4(d)是0.51˚仰角上9:00的差传播相位常数KDP图:
Figure 4. S-band dual polarization radar reflectivity factor (a) plot, velocity (b) plot, differential reflectivity (c) plot, differential propagation phase constant (d) plot in Guangzhou on May 9, 2016
图4. 2016年5月9日广州S波段双偏振雷达反射率因子(a)图,速度(b)图,差分反射率(c)图,差传播相位常数(d)图
在图4速度(b)图中明显看到在方位角340˚方向,100 km处以及西北部75 km处有正负速度对,且核心区域对应的反射率因子超过60 dBZ,对速度数据进行中气旋识别并展示TVS (龙卷涡旋特征):
Figure 5. Results of medium cyclone identification algorithm
图5. 中气旋识别算法结果
根据中气旋识别算法,图5识别出在雷达正北方向100 km处两侧,方位角315˚方向,75 km处及方位角290˚,100 km处存在中气旋,由于中气旋识别算法算出的切速度小于15 ms−1∙km−1,需要进行下一步单体联合识别:
Figure 6. Single-unit joint identification map (yellow markers are tornado areas)
图6. 单体联合识别图(黄色标记为龙卷区域)
根据进一步的单体联合识别显示,如图6在黄色标记(−45.27, 58.08) km处,即方位角319˚方向,78 km处有龙卷发生,该处位于强Z,强KDP,强ZDR交界处,接下来对过该处的径向及相邻径向进行剖面分析。
图7(a)中反映出在78 km处,有单体生成,40 dBZ回波顶高高达8 km,6 km高度处出现大于70 dBZ的强反射率。图7(b)中78 km处有高达四公里的强速度垂直区域,与反射率因子60 dBZ区域相对应。在该径向两侧,即图7(c)和图7(d)中,同78 km处,出现强ZDR柱及KDP柱。该特征是由于龙卷内部的粒子分选机制造成的,通过四处径向的分析,判断方位角319˚方向,78 km应为龙卷的发生地。
Figure 7. (a) Radial profile of the reflectivity factor at azimuth 319.8˚, (b) radial profile of the velocity at azimuth 318.6˚, (c) radial profile of the differential reflectivity (ZDR) at azimuth 317.6˚, and (d) radial profile of the differential propagation phase constant KDP
图7. (a) 方位角319.8˚的反射率因子径向剖面图,(b) 方位角318.6˚速度径向剖面图,(c) 方位角317.6˚的差分反射率(ZDR)径向剖面图,(d) 差传播相位常数KDP的径向剖面图
4.3. 识别效果评估
针对灾害的识别效果评估,给出四个指标,分别是临界成功指数ICS,命中率RH,漏报率RM和虚警率RFA,其对应的计算方式如(6)式所示:
,
,
,
(6)
上式中,x为灾害案例中识别成功的样本个数,y为实际存在气象灾害但未识别出来的样本个数,z为识别出灾害但实际并未发生的样本个数。
在本次龙卷识别实验中,三水龙卷案例的临界成功指数和命中率为0.67,漏报率为0.33,阜宁龙卷案例的临界成功指数为0.36,漏报率为0.64,见表2,三水龙卷案例的龙卷识别效果明显优于阜宁龙卷案例,证明在识别出中气旋的基础上,引入双偏振参量进行龙卷的识别能够较大提高龙卷的识别预警率。总的来说,本次龙卷识别实验的临界成功指数和命中率为0.5,漏报率为0.5,未来算法仍需改进。
Table 2. Tornado identification effect evaluation table
表2. 龙卷识别效果评估表
5. 结论与讨论
本文研究了分别利用单偏振以及双偏振雷达数据基于ZV、ZW的相关性,中气旋识别算法叠加DBSCAN聚类算法的龙卷识别方法,并通过江苏阜宁的一次强对流天气过程以及广东三水地区的一次龙卷实例验证了该识别方法的有效性。
1) 目前,对于龙卷的识别比较艰难,很难找到准确的识别方法,因此将传统的中气旋识别算法叠加双偏振参量算法能够提高识别准确率。
2) 对于龙卷灾害性天气而言,ZDR弧具有警示作用,该特征在龙卷发生前10~15分钟发生,持续至龙卷结束,因此,对ZDR弧的识别能够提前预警龙卷。同时,KDP足和ZDR弧的识别能够较好表征低层入流强度和切变的情况。
3) 对于DBSCAN聚类算法,该算法能够很好的识别单体并且不用提前设置单体数量,同时,任意形态的单体都能识别出来,这对于ZDR弧,带状ZH单体具有很好的识别作用。
但本文仍有一些不足,1) 本文使用个例较少,该方法是否具有普适应性还有待考证。2) 垂直相关的方法还可以进一步改进。3) 本文使用的ZDR数据并未进行预处理,未来可根据实际情况探究是否需要对该数据进行预处理。