一种基于VLF/LF闪电数据的多算法融合雷暴系统识别技术

2024-01-12伍华丽植耀玲卢炳夫曾鹏

热带气象学报 2023年5期

伍华丽，植耀玲，卢炳夫，曾鹏

(广西壮族自治区气象灾害防御技术中心，广西南宁 530022)

1 引言

雷电作为一种能引起严重灾害的自然现象，能造成人员伤亡，损毁电气及电子通信设备、建筑物，诱发火灾和爆炸事故，威胁人们的生命和财产安全，是联合国国际减灾委员会公布的最严重的十种自然灾害之一，是国家防灾减灾工作重点关注的领域。因此，开展雷电临近预警预报关键技术研究工作对保护人民的生命财产安全有着非常重要的意义。雷电临近预警预报常通过雷暴识别、追踪及外推来实现，雷暴系统识别的准确与否直接影响后续追踪及趋势外推的效果，是做好雷电临近预警预报的关键。

学者们常利用雷达数据开展灾害性天气系统的识别技术研究，吕伟涛等[1]利用雷达资料对可能发生或已经发生闪电的区域进行识别，并基于识别结果进行追踪和外推，实现雷电临近预警。王萍等[2]采用分类决策树法，基于多雷达参数构建50 km以内雷暴系统分类识别模型，较好地实现了小尺度雷暴系统中短时降水和冰雹之间的分类识别。Kelly 等[3]利用傅立叶变换算法、基于雷达反射率图像的边缘检测以及形态分析，识别灾害性天气。

相对于雷达数据，闪电定位数据具有更高的监测频率和更低时延，更有利于识别快速生消的中小尺度雷暴系统，也有一些学者基于二维(ADTD)闪电定位系统资料利用聚类算法开展雷暴系统的识别研究，侯荣涛等[4]基于改进DBSCAN聚类算法(IDBSCAN)对江苏地区一次雷暴过程中的闪电定位数据进行聚类分析，识别地闪聚类簇及其中心并基于聚类中心的变化规律实现了对闪电移动趋势的预报。周康辉等[5]利用密度极大值快速搜索聚类法，基于闪电定位数据实现了对雷暴的识别，有效处理雷暴分裂与合并的情况。侯荣涛等[6]、路郁[7]运用改进的OPTICS 聚类算法对闪电定位资料进行分时段聚类分析并构建了雷电临近预警模型。上述研究充分印证了基于闪电定位数据开展雷暴系统识别的可行性，为基于闪电定位数据的雷暴系统识别技术提供了理论支撑，同时也指出了目前闪电定位数据在雷暴系统识别应用方面存在的一些问题，如：传统聚类法(如DBSCAN、K-Means 等)对参数敏感，无法有效处理密度差异较大的雷电数据集，不适用于多尺度的雷暴系统；闪电定位数据的数据量及密集程度影响着聚类效果的好坏，仅针对强雷暴及强对流天气的情形识别效果较佳；当被探测到的闪电定位数据稀少时，无法通过聚类来识别雷暴系统。

近年来随着VLF/LF 三维闪电监测网的高速建设[8]，源自该监测网的闪电定位数据在北京[9]、湖北[10]、四川[11]、江苏[12]、浙江[13]、福建[14]、山西[15]、贵州[16-17]、云南[18]、海南[19]、广西[20]等多地气象部门得到业务化应用。针对VLF/LF 三维闪电监测网闪电定位数据(下文简称VLF/LF闪电定位数据)的研究，均集中在闪电活动特征分析、与二维(ADTD)闪电定位系统资料对比分析等方面，研究结果表明VLF/LF 三维闪电定位系统对地闪的探测效率及定位精度均高于ADTD 系统，探测的地闪数是ADTD 系统的2～3 倍[10,13-16]；两套系统监测到的雷电活动时空分布趋势比较一致。目前，针对VLF/LF 闪电定位数据的雷暴系统识别及预警应用研究开展较少，鉴于数据量级和密度分布上的差异性，现有的闪电定位数据的聚类识别技术研究成果无法适应于VLF/LF闪电定位数据，迫切需要开展基于VLF/LF 闪电定位数据的雷暴灾害性天气识别研究，深入挖掘VLF/LF三维闪电监测网在雷电临近预警中的应用效益。

针对上述需求，本文提出一种多算法融合雷暴系统识别技术，该技术首创一种基于密度特性自适应的改进的OPTICS 聚类法，改善传统OPTICS 聚类法无法处理影响雷暴分布识别的稀疏点、无法识别区分密度差异不大的相邻雷暴体的缺陷；同时，与K-Means 聚类法融合应用，解决OPTICS 聚类法[6,21]等无法准确定位雷暴体的中心点的难题，实现任意尺度形状雷暴体的识别及其几何中心点的精准计算，可为后续雷暴追踪及趋势外推工作提供良好的雷暴系统识别结果，也能为雷暴灾害天气的预报及雷电预警模型的训练提供有效的知识库支撑。

2 资料及来源

2.1 VLF/LF闪电定位数据

本文采用的闪电定位资料来源于广西气象部门VLF/LF 三维闪电监测定位系统监测的闪电定位数据[22]。该系统由中国科学院电工研究所制造[8]，包含20个三维闪电探测站(ADTD-2C 型闪电定位仪)和1 个数据处理中心，基线距离大于70 km，子站密度高，有效探测范围覆盖广西全域。系统采用三维时差定位法(3D-TOA)，能对云地闪回击、云闪闪击进行定位，时间精度优于10-4s，闪电强度相对误差小于10%(在10～25 kA 范围内)，10 kA 以上闪电探测效率大于90%，网内平均水平定位误差小于300 m，垂直定位误差小于500 m[8,23-24]，可探测闪电位置、高度、强度、时间、类型以及极性等重要参数[8,25]。

2.2 雷达数据

本文所用的雷达数据，是中国气象局研制的SWAN 系统(Severe Weather Automatic Nowcast System)[26]提供的雷达组合反射率拼图数据。组合反射率因子表示的是垂直空气柱内的最大反射率因子在笛卡尔坐标水平面上的投影，是各层较强回波的迭加反映，能较明显地反映对流天气的强度及发生位置。SWAN系统对广西境内10部雷达(南宁、柳州、桂林、百色、河池、崇左、梧州、防城港、北海、玉林)基数据进行实时采集解析，经过数据质量控制、坐标系转换(极坐标系与笛卡尔坐标系插值转换)、重叠区域融合拼接等流程后，生成网格化的雷达拼图数据集，包括组合反射率拼图、回波顶高拼图、VIL 拼图等。其中，组合反射率拼图数据以dBZ 为单位存放，水平分辨率为0.01 °×0.01 °(约1 km×1 km)，时间分辨率为6 min，空间范围是102～113 °E，19～28 °N。

3 多算法融合雷暴系统识别技术

3.1 主要技术路线

多算法融合的雷暴系统识别技术的主要技术路线如图1 所示。首先对传统的OPTICS 法进行参数自适应改进，即基于VLF/LF 闪电定位数据的空间分布特征实时计算初始输入参数Min_ρ，使其自适应于待处理的闪电数据，有效解决固定的初始输入参数因缺乏普适性而影响聚类效果的问题；而后与K-Means 聚类法融合应用，计算雷暴系统的中心点，利用K-Means 法通过迭代计算聚类中心的优势，改进仅凭OPTICS法无法直接准确定位雷暴体的中心点(聚类中心)的缺陷。

图1 多算法融合雷暴系统识别技术算法流程图

3.2 OPTICS聚类法简介

OPTICS(Ordering Points to identify the clustering structure)算法是一种基于密度的空间聚类算法。该算法将具有足够密度的区域划分为簇，并在含有噪声的空间数据库中发现任意形状的簇。

3.2.1 算法涉及的基本概念

ε邻域：给定对象半径为ε内的区域称为该对象的ε邻域。

核心点：如果给定对象ε邻域内的样本点数大于等于Min_ρ，则称该对象为核心点。

直接密度可达：对于样本集合D，如果样本点q在p的ε邻域内，并且p为核心点，那么对象q从对象p直接密度可达。

密度可达：对于样本集合D，给定一串样本点p1，p2，……，pn，p=p1，q=pn，假如对象pi从pi-1直接密度可达，那么对象q从对象p密度可达。

密度相连：存在样本集合D 中的一点o，如果对象o到对象p和对象q都是密度可达的，那么p和q密度相联。

核心距离：是使一个样本点p0 成为核心对象的最小半径，在给定邻域半径ε和Min_ρ参数的前提下，核心距离cd(p0) ≤邻域半径ε。

可达距离：核心点p0 与其邻域内的点p1 的距离d(p0，p1)与核心距离cd(p0)中较大的称为p1 到p0 的可达距离rd(p1，p0)。易知，核心距离内邻域点到核心点的可达距离为固定值cd(p0)，核心距离外邻域点到核心点可达距离为真实距离d(p0，p1)。

3.3 OPTICS 聚类法改进方案

由3.2 节中的OPTICS 聚类法原理可知，OPTICS 聚类法初始参数的设置主要用于筛选参与聚类的点(核心点)，剔除离散点。当筛选出参与聚类的点后，依靠Min_ρ来确定扩张半径，若该点处于密集区，则扩张半径就很小，若处于稀疏区，扩张半径就大，基于Min_ρ确定的扩张半径不是一个特定的值，而是一个范围，邻域ε的大小决定了这个范围值的上限，则聚类结果对于邻域ε不再敏感，邻域ε只要在合理的范围内设定足够大，就能取得不错的聚类效果。

由于OPTICS聚类法对于初始输入参数邻域ε不再敏感，则聚类效果的好坏主要取决于Min_ρ的选取。改进的OPTICS聚类法与传统的OPTICS聚类法相比，区别在于初始输入参数Min_ρ不会取固定值，而是根据闪电空间分布特性实时计算获取，具有自适应的特性，随实时闪电空间分布特性的变化而变化，自适应于待处理的闪电数据。

改进算法的流程如图3所示。

图3 改进的OPTICS算法流程图

密度特征参量ρ0根据闪电密度空间分布特性实时提取，计算步骤如下。

(1) 对VLF/LF 闪电定位数据进行网格化处理，即将广西区域划分为10 km×10 km 的网格，计算每个网格内VLF/LF闪电定位数据的频数。

(2) 提取闪电频数大于0 的网格，将网格按闪电频数大小降序排列，得到序列Kn。其中n为降序排列后给每个网格标注的序号(n=1，2，3，……，N，N为所含闪电频数大于0 的网格的总数)，Kn为序号n对应的网格内VLF/LF 闪电定位数据的总频数。

(3) 对排序前n的网格数量进行累加得到X(n)，即：

(4) 将排序前的n个网格内所含的闪电数量进行累加得到Y(n)，即：

其中，m1为序号1的网格内VLF/LF闪电定位数据的总数，以此类推，mn为序号n的网格内VLF/LF闪电定位数据的总数。

(5)以X(n)为横坐标，Y(n)为纵坐标绘制曲线，采用双总体t检验法，计算该曲线的曲率突变点(X(m)，Y(m))，进而得到对应的Km，并通过公式ρ0=πε2×Km/(10×10)计算密度特征参量ρ0（图4）。

图4 双总体t检验法计算流程图

其中，选用的t检验统计公式为：

3.4 基于K-Means法的雷暴中心点提取

K-Means 法是一种基于划分的聚类算法，以距离作为数据对象间相似性度量的标准(即数据对象间的距离越小，则它们的相似性越高)，将相似的对象归于一个簇中，簇中心通过簇中所有点的均值来计算，通过迭代计算寻找最佳的分类结果(簇集合)及其聚类中心。

本文利用K-Means 法通过迭代计算聚类中心的优势来计算OPTICS法得到的簇类的几何中心，算法的流程描述：(1) 根据OPTICS 聚类法的输出结果初始化簇数量值k及其对应的初始聚类中心(x1，y1)，(x2，y2)，……，(xk，yk)；(2) 针对VLF/LF 闪电定位数据集中每个闪电计算得到k个聚类中心的距离并将其分到距离最小的聚类中心所对应的簇类中；(3) 针对每个簇，重新计算它的聚类中心；(4) 重复(2)和(3)，直到聚类中心不再变化，这时就产生了最后的分类结果；(5) 输出聚类簇集合以及相应的聚类中心。

4 结果分析

本文提出的多算法融合雷暴系统识别技术，旨在解决两个主要问题，一是探索适用于VLF/LF闪电定位数据的聚类算法；二是基于聚类结果研制雷暴系统识别产品，为雷暴系统的追踪外推提供产品支撑，深入挖掘VLF/LF三维闪电监测网在雷电临近预警中的应用效益。综上，本小节分别从闪电数据的聚类效果及雷暴识别结果的优劣两方面进行检验分析。

4.1 实验数据

为了检验本文算法对不同雷暴天气过程的适应性，选取不同类型雷暴天气过程的VLF/LF闪电定位数据作为检验数据。从聚类效果、单个时段的识别结果、连续时段雷暴系统运行轨迹识别结果三个方面进行检验。

根据雷暴过程雷达回波发展形态，参考相关研究[27-29]，可将广西常见的雷暴天气过程分为三类：(1) 区域性雷暴天气过程，即为与大尺度对流系统相关的雷暴天气过程(如飑线、超级单体等系统)；(2) 局地性雷暴天气过程，即为局地突发的小尺度雷暴天气系统；(3) 分散性雷暴天气过程，即为分散性多单体的中小尺度雷暴天气过程。选取上述三类雷暴天气过程的个例进行检验实验，详情如表1所示(表中时间为北京时间，下同)。

表1 雷暴天气过程个例属性表

4.2 聚类效果对比

分别利用传统OPTICS 聚类法及本文提出的多算法融合雷暴系统识别技术，对个例1～3 的检验数据进行处理，结果如图5～图10 所示，其中图5a～图10a 为OPTICS 聚类法输出的距离可达图，图中的每一个凹陷部分代表一个聚类簇，横坐标代表经过聚类后被识别为簇的散点个数，纵坐标为散点之间的距离；图5b～图10b 为基于可达图进一步提取的雷暴聚类结果的GIS 展示图。特别说明的是，对于常规OPTICS 聚类法，本文参考了相关文献[4, 6, 21]的阈值确定方法，经多次实验后，选取初始输入参数Min_ρ=15、ε=10。

图5 区域性雷暴天气过程(2020年7月10日19:20—30)VLF/LF闪电数据基于本文算法的处理结果a的横坐标为排序好的样本，下同。

由图5a 和图6a、图7a 和图8a、图9a 和图10a对比可知，对于不同类型雷暴灾害性天气，利用传统OPTICS聚类法得到的可达图形态不规则，不利于制定簇群提取规则，进而影响簇群提取的效果，得到不佳的聚类结果，例如图6b所示的结果，将本该属于一个雷暴体的闪电数据识别得四分五裂，输出多个雷暴体簇，与实际严重不符；例如图8b所示的结果，不仅将闪电数据过度过滤，还将本应属于两个雷暴体的数据识别为一个雷暴体；例如图10b所示的结果，将本应属于一个雷暴体的数据识别为两个雷暴体。而本文提出的多算法融合雷暴系统识别技术的输出结果如图5、图7 以及图9 所示，初始输入参数Min_ρ及ε由实时VLF/LF 闪电定位数据估算动态获取，得到的可达图形态规则，凹陷和凸起的位置分明(图5a、图7a、图9a)，最后的簇群提取效果良好(图5b、图7b、9b)，有效滤除影响雷暴体识别结果的离散闪电数据，很好地识别区分密度差异不大的相邻雷暴体。

图6 区域性雷暴天气过程(2020年7月10日19:20—30)VLF/LF闪电数据基于OPTICS法的处理结果

图7 分散性雷暴天气过程(2020年7月2日16:40—50)VLF/LF闪电数据基于本文算法的处理结果

图8 分散性雷暴天气过程(2020年7月2日16:40—50)VLF/LF闪电数据基于OPTICS法的处理结果

图9 局地性雷暴天气过程(2020年7月9日17:40—50)VLF/LF闪电数据基于本文算法的处理结果

图10 局地性雷暴天气过程(2020年7月9日17:40—50)VLF/LF闪电数据基于OPTICS法的处理结果

4.3 雷暴识别结果分析

利用雷达组合反射率拼图数据，分别从某个时刻的识别结果以及一段时间内识别得到的雷暴运行轨迹两个角度，对本文算法得到的雷暴识别结果进行检验。

4.3.1 单个时刻的识别结果分析

将本文算法的识别结果与同时刻的雷达组合反射率拼图数据叠加展示（图11、图12 及图13），针对不同类型的雷暴天气过程，本文提出的算法能实现任意尺度形状雷暴体的识别，识别出的雷暴体均与雷达回波图中对流单体位置具有较好的一致性。

图11 区域性雷暴天气过程(2020年7月10日19:20—30)本文算法输出的雷暴聚类识别结果与雷达数据的叠加图

图12 分散性雷暴天气过程(2020年7月2日16:40—50)本文算法输出的雷暴聚类识别结果与雷达数据的叠加图

图13 局地性雷暴天气过程(2020年7月9日17:40—50)本文算法输出的雷暴聚类识别结果与雷达数据的叠加图

4.3.2 雷暴运行轨迹识别结果分析

图14 给出了2020 年7 月10 日18:10—21:10广西中南部的一次雷暴天气过程中，由本文算法识别出的雷暴体移动路径分布。雷暴A、B均起始于18:10，均向东北方向移动，在19:10 的时候，雷暴A分裂成单体A1和A2。单体A1继续向东北方向移动。单体A2 向南移动，并于20:20 与一直持续向东北方向移动的雷暴B 合并成雷暴体C。此后直至21:10 雷暴体A1 与C 均持续向东北方向移动。

图14 雷暴运行轨迹识别结果与2020年7月10日18:10—21:10闪电数据的对比

总体而言，整个雷暴天气过程中监测到的闪电散点呈较明显的折线状分布，与雷暴路径能较好地重合，表明本文算法能识别得到的雷暴运行轨迹，并且能有效识别雷暴的分裂与合并的情况，能对雷暴未来移动路径的追踪外推提供较好的识别产品支撑。

将分裂前(18:20)、分裂后(19:20)以及合并后(20:20)三个时刻的雷达组合反射率拼图与雷暴运行轨迹识别结果进行叠加对比，检验雷暴体运行轨迹识别效果，结果如图15～17 所示。本文算法识别出的雷暴运行轨迹与雷达强回波中心能较好重合；雷暴运行轨迹和雷达回波均显示了雷暴系统整体持续向东北方向运动的趋势；本文算法识别出的雷暴A 的分裂以及雷暴A2 与B 的合并情况，与雷达回波显示的单体先分裂后合并在时间和空间上均具有很好的一致性，表明本文算法能有效识别雷暴的分裂及合并。

图15 2020年7月10日18:20雷达组合反射率拼图与雷暴运行轨迹识别结果的叠加对比

图16 2020年7月10日19:20雷达组合反射率拼图与雷暴运行轨迹识别结果的叠加对比

图17 2020年7月10日20:20雷达组合反射率拼图与雷暴运行轨迹识别结果的叠加对比

5 结语

本文为深入挖掘VLF/LF 闪电数据在雷暴灾害性天气识别预警方面的天然优势，针对目前传统密度聚类法(如OPTICS 法等)不适用于VLF/LF闪电数据聚类处理、无法识别区分密度差异不大的相邻雷暴体、无法准确定位雷暴体的中心点等缺陷，提出一种多算法融合雷暴系统识别技术，并利用不同尺度的雷暴天气系统过程数据对其进行检验，结论如下。

(1) 本文算法利用双总体t检验法基于实时闪电定位数据密度分布特征对OPTICS 法进行初始输入参数自适应改进，较常规OPTICS 法，能输出形态规则、无毛刺感，凹陷和凸起的位置分明的可达图，有利于雷暴簇的自动提取，滤除影响雷暴体识别结果的离散闪电数据，识别区分密度差异不大的相邻雷暴簇。

(2) 本文提出的多算法融合雷暴系统识别技术能较好地适应不同类型雷暴系统，实现任意尺度形状雷暴体的识别及几何中心的较精准定位，可为雷暴系统未来移动路径的追踪外推提供较好的识别产品支撑，一方面有助于填补由于山脉和其他遮挡物造成的雷达回波衰减区的监测空白，另一方面可作为常规风暴识别产品的补充，与基于雷达的风暴识别产品联合应用，有利于实现针对雷暴系统的风暴分类识别。

(3) 本文旨在探索适用于VLF/LF 闪电数据的聚类方法，仅针对聚类算法的创新性和可行性进行讨论，后续将考虑对该算法的应用效益进行深入挖掘，充分利用VLF/LF 闪电数据的三维属性，按高度进行分层聚类进而反演雷暴系统三维结构，研制具备三维属性的雷暴系统识别产品。