引言

雷暴活动中始终伴随着强烈的放电现象，闪电活动能很好地反映雷暴活动的强弱变化与移动趋势^[1-3]。近年来，随着国家雷电监测网络的建立与完善，使全国范围的云-地闪电 (简称闪电) 监测成为可能。相对于雷达与卫星数据，闪电观测数据具有更高的观测实时性与更低的传输时延，对于实时监测快速生消的中小尺度对流系统具有非常重要的意义^[4-5]。目前，国内外基于雷达的对流识别、追踪与外推的算法研究较多 (如TITAN，SCIT等^[6-9])，产品也得到较为广泛的应用，对强对流天气的临近预报有较好的指示性作用。

基于闪电数据的雷暴识别、追踪与外推算法也有一些尝试。侯荣涛等^[10]构建了江苏地区基于DBSCAN (Density-Based Spatial Clustering of Applications with Noise) 聚类算法的闪电簇识别、追踪与临近预报模型；Tuomi等^[11]开发了利用时空阈值区分闪电簇进而识别并追踪雷暴的算法，分析了芬兰地区的雷暴活动特征。Betz等^[12]利用基于闪电密度判别的雷暴追踪算法分析了LINET (VLF/LF lightning detection network) 闪电定位网络探测的闪电数据。Kohn等^[13]将闪电定位数据引入预警决策支持系统 (WDSS-Ⅱ)，利用闪电密度阈值完成雷暴识别与追踪，取得较好的短临预报效果。Bonelli等^[14]利用雷达与闪电定位数据，通过设置时空阈值实现雷暴路径的识别、追踪与外推，并能够实现30 min或更长时间的外推。姚叶青等^[15]通过叠加分析雷电集中区与风暴，发现雷电集中区一般对应有风暴，且雷电集中区中心在风暴中心附近，二者移动方向基本一致，移速接近；通过匹配雷电集中区与风暴，采用雷电集中区将随其所匹配的风暴一起移动的思路，利用风暴追踪技术，对雷电集中区进行临近外推，从而实现雷电的临近落区预报。吕伟涛等^[16]利用雷达、卫星、闪电等数据，开发了区域识别、追踪和外推算法AITEA，可对已发生闪电的区域 (如回波强度超过某个阈值或云顶亮温低于某个阈值的区域等) 进行识别, 利用一段时间的监测资料便可进行跟踪，采用Holt双参数线性指数平滑方法，对区域的中心位置坐标进行预测。

本文介绍一种基于全国地闪数据的雷暴识别、追踪与外推方法。方法中使用全新、高效的闪电聚类方法实现雷暴识别，利用Kalman滤波技术实现追踪并完成对雷暴路径0~60 min外推，可实现对全国范围雷暴活动的监测与临近预报。

1 雷暴识别、追踪与外推 1.1 雷暴的识别

由于雷暴过程中闪电往往集中发生于雷暴云中对流最旺盛的区域，因此，雷暴在闪电数据上的表现形式即为一簇密集的闪电。利用聚类算法可实现雷暴识别，本文中聚类算法采用密度极大值快速搜索算法 (clustering by fast search and find of density, CFSFD)^[17]。相对于传统K-means与DBSCAN等经典聚类算法而言，CFSFD不需要预指定聚类中心，同时对非球面形状的雷暴能有更好地识别效果。该算法只考虑点与点之间的距离，因此，不需要将点映射到一个向量空间中，算法复杂度较传统算法有较大改进。

闪电密度计算。

(1)

式 (1) 中，；n为聚类的闪电数量；d_i，j为两个闪电i，j之间的距离；d_c为间隔阈值，根据文献[17]，d_c的取值应使d_i，j < d_c的闪电数量约占总数的2%。

利用式 (1) 计算每个闪电的密度，即计算每个闪电方圆距离d_c内的闪电数量。闪电密度值是雷暴中心识别的重要依据，闪电密度越大，其闪电越密集，表征越强烈的放电过程。

(2)

实际聚类过程中，也可使用高斯核 (Gaussian Kernel) 函数计算闪电密度 (式 (2))。高斯核函数从中心到外围根据距离指数衰减，更易确定唯一的雷暴中心点。

闪电距离计算。

(3)

利用式 (3)，对每个闪电，计算所有其他闪电中密度比其大的闪电距离该闪电的最小距离 (对于闪电密度最大的闪电，其δ=max(d_ij))。对于δ越大的闪电，其周围散乱点越少，某一区域上簇状独立性越高。

雷暴聚类中心的确认。假设图 1中有两个雷暴，图 1中所有闪电的密度值按照由高到低排列，“1”表示密度最高的点，“2”次之，以此类推。图 1b给出了每个闪电的密度与归一化距离分布，纵坐标为相对距离比值。雷暴中心闪电的确认可以通过给定的δ_min和ρ_min筛选出同时满足ρ > ρ_min和δ>δ_min条件的点作为距离中心点。图 1中闪电1与闪电10可以作为雷暴闪电中心。闪电2~8虽然密度大于闪电10，但其δ < δ_min，因此，不能作为雷暴中心闪电。闪电26~28虽有较大的δ值，但其ρ < ρ_min，同样不能作为雷暴聚类中心。

雷暴属于中小尺度天气系统，综合考虑雷暴尺度与实际聚类效果，本文取ρ_min=1.5，δ_min=20 km，对各类雷暴单体闪电簇有较好的识别效果。

其余闪电的分配。当雷暴中心闪电确定后，剩余闪电的类别标签按照以下原则指定：闪电的类别标签与高于该闪电密度的最近闪电的类别一致。雷暴边界的确定。定义雷暴边界区：某一雷暴的雷暴边界区，由该雷暴中与其他雷暴任意闪电距离小于d_c的闪电构成。寻找雷暴边界区中密度最大的闪电，将其密度记为ρ_max。将该雷暴中ρ < ρ_max的闪电作为噪声去除，从而确定雷暴边界。

1.2 雷暴的追踪与外推

Kalman滤波器能在线性无偏最小方差估计准则下对路径做出最优估计^[18]。在实现雷暴识别的前提下，假设雷暴线性移动，可利用Kalman滤波算法对雷暴展开线性外推。

对雷暴的移动与发展过程作如下假设^[6]：雷暴移动路径为直线, 雷暴强度呈线性增强或减弱。

1.2.1 Kalman追踪与外推模型

关于Kalman滤波器原理可参见文献[18]。根据雷暴移动特征，建立Kalman状态滤波模型。

系统状态方程和观测方程后利用卡尔曼滤波方程式通过递推方法，不断预测目标在下一时刻的位置，并实现0~60 min的路径外推。

1.2.2 雷暴的追踪

利用Kalman滤波器输出的下一时刻雷暴位置，同时通过设置相应时间阈值ΔT与空间阈值ΔD，实现对雷暴的追踪。具体追踪过程如下：

① 对于t₀时刻已有雷暴，Kalman滤波器输出其下一时刻的预测位置 (t₁，P₁)。②搜索新识别的雷暴中对 (t₁，P₁) 满足时空阈值 (ΔT，ΔD) 的雷暴，如有多个雷暴满足条件，选择距离P₁最近的雷暴，认为其为雷暴路径的延续。③若新识别的雷暴未能匹配到路径，则认为是新生的雷暴。④雷暴的合并表示几个雷暴合成一个雷暴的过程 (图 2a)。对于t₁时刻的单体1与单体2，如果t₂时刻只监测到1个单体3，若单体1与单体2在t₂时刻的Kalman外推落在单体3的单体范围之内，且单体3、单体1与单体2均满足时空阈值限定，则认为单体3为单体1与单体2的合并。⑤雷暴的分裂表示1个大的雷暴分裂成几个小的雷暴的过程 (图 2b)。类似于雷暴合并的处理，若单体1下一时刻分裂为单体2与单体3，且其中心位置落在单体1的Kalman外推区域之内，则认为单体2与单体3为单体1的分裂。

1.2.3 雷暴的外推

在雷暴历史路径的基础上，根据1.2.2节所述的建模过程，利用Kalman滤波器不难得到其0~60 min的外推矢量。将雷暴面积作为雷暴发展与消亡的特征参数，可线性外推雷暴的强度变化。综合以上雷暴识别、追踪与外推的步骤，总流程如图 3所示。

2 数据来源

闪电定位数据来自国家雷电监测网提供的2013年地闪观测数据。国家雷电监测网2013年投入业务考核的站点达到347站，较好地覆盖了除青藏高原、内蒙古中西部地区外的全国大部分地区。闪电定位信息包括地闪回击二维空间信息、发生时间、闪电强度、闪电陡度、定位方式、定位误差等丰富信息。整套系统网内定位精度优于300 m，探测效率不低于80%^[19]。经过长期的业务化运行，整套系统较为稳定，能够不间断地对覆盖范围进行地闪监测。

3 数据处理结果

利用上述雷暴识别、追踪方法，处理2013年全国的雷电监测数据，数据处理时间间隔为10 min (即将每10 min累积的闪电进行识别、追踪与外推)，最终得到雷暴70840个。

3.1 雷暴的识别结果

图 4为2013年3月19日21:20—21:30 (北京时，下同) 位于江西东北部、浙江西部的闪电聚类 (即雷暴识别) 效果。图 4a显示对于不同形状分布、不同闪电数量的闪电簇，CFSFD聚类算法均能按照其空间位置分布实现较好地聚类。将识别结果与雷达回波 (图 4b) 进行对比发现，由闪电聚类、识别的雷暴活动区与雷达回波图中的对流单体位置具有较好的一致性。

由此可见，本文的聚类算法能够对闪电进行聚类，进而有效识别雷暴活动范围。

3.2 雷暴追踪结果

图 5给出了2013年3月20日华南地区6个雷暴的路径分布。雷暴A起始于20：40，向东偏南方向发展，持续至22:50，雷暴生命史为2 h。雷暴B起始于12:50，开始向北发展，而后向东移动。雷暴在13:20分裂，主单体继续向东移动，共持续4 h，移动距离超过350 km；分裂出的小单体，向东北方向移动，持续至13:50。雷暴C起始于12:40，雷暴生成后向东偏北方向发展，13:50雷暴出现分裂，分裂子雷暴向东北方向移动，生命史较短，迅速消散；雷暴主体继续向东移动，维持至15:50。雷暴D具备局地性较强、迅速生成并消散的特征，从生成至消散，整个过程约30 min，雷暴移动距离不超过15 km。利用闪电识别的雷暴E, F持续时间短，仅维持20 min，为较弱的雷暴过程。

总体而言，闪电分布呈较明显的线状，雷暴路径与闪电线状分布能有较好地重合，表明雷暴路径的追踪具有较好的效果。

为了验证雷暴的识别与追踪效果，图 6利用雷达回波对其进行验证分析，由图 6可清晰地观察到对流单体的发生、发展特征。

12:50雷达回波显示，雷暴B, C对应的对流单体已生成，此时对流的最大反射率因子与面积均较小，明显处于新生阶段。此后雷暴不断发展，雷暴面积增加，中心最大反射率因子增大。13:30雷达回波显示雷暴B对应的单体出现分裂，子单体位置与B1，B2位置一致，B2继续东移，B1向东北方向移动；与此同时，单体C出现分裂迹象。雷暴E对应到明显的对流单体。14:00雷暴B1只持续至13:50，此时雷达回波显示B1对应的对流单体回波强度减小，结构松散，呈消散状态，表明二者基本一致；雷暴追踪结果显示雷暴C于13:50分裂，与此相对应，雷达回波显示14:00雷暴C对应的单体分裂已完成，子单体位置与C1，C2位置一致。雷暴F对应的对流单体出现。14:40—17:50雷暴B，C分裂后的主体继续向东移动，雷暴与对流单体，二者移动轨迹基本重叠。在雷暴追踪结束之后，雷达回波显示对流单体仍然存在，但很快消失。

值得注意的是，12:50—13:10雷达反射率因子拼图显示较多的对流单体。相比之下，利用闪电仅识别了4个雷暴 (雷暴B，C，E，F)。究其原因，其余单体最大反射率因子较小，对流层顶高度较低 (图略)，单体的持续时间短、移动距离小，表明较弱的对流活动不利于闪电活动的发生^[20-21]。

雷暴路径与雷达反射率因子拼图中的最大反射率因子中心能较好地重合，表明定位到的闪电多发生于雷暴中对流最旺盛的阶段；根据闪电数据识别并追踪得到的两个雷暴的移动路径与雷达回波显示的雷暴路径有较好的一致性，表明方法有效。

雷暴B，C的分裂的追踪结果与雷达回波显示的单体分裂具有很好的一致性，合并的原理与分裂相似，表明方法对雷暴单体分裂与合并的有效性。雷暴提前于对流单体消失，时间提前量约为10~60 min。这说明闪电对对流单体的消散具有一定指示意义。

较弱的对流单体中，闪电活动较弱或无闪电活动。雷暴E，F实际上是同一对流系统的不同发展阶段，然而由于该对流单体中闪电活动较弱，导致被识别为两个独立的雷暴。由此可见，雷暴与对流单体仍存在一定差异。

3.3 雷暴外推结果评估

为了评估外推算法的性能，利用类似Dixon等^[6]的外推评估方法，即将路径与外推结果按5 km×5 km的网格化处理，如果有闪电出现在某网格内，则认为该网格为雷暴活动区。定义以下术语^[6]：命中指外推路径与实际路径均为雷暴活动区；漏报指实际路径为雷暴活动区，外推路径为非雷暴活动区；空报指外推路径为雷暴活动区，实际路径为非雷暴活动区。

命中率

(4)

漏报率

(5)

临界成功指数

(6)

其中，N为样本量。

为了验证雷暴外推的效果，本文选取了雷暴活动活跃的2013年3月19—20日作为检验时段，该时段内共识别、追踪和外推雷暴623个，外推检验结果如表 1所示。

检验结果显示：类似于TITAN等外推算法，基于Kalman的雷暴外推可靠性随时间迅速衰减。这可能是因为风暴的生命史太短, 在预报期间内生长和消亡的变化太快导致的^[22]。表 1显示，10 min的外推结果其外推命中率约为0.7，30 min则降至0.6，60 min继续降低；与此同时，虚报率随着时间增加从0.5以下增加至0.7；临界成功指数从0.4降低至0.1~0.2。

总体而言，对比TITAN等算法，本文开发的外推算法在各项检验参数上整体表现接近，某些方面有更好的表现。

4 结论与讨论

本文介绍了一种基于全国地闪数据的雷暴识别、追踪与外推的新方法。该方法利用密度极大值快速搜索聚类算法实现雷暴识别，利用Kalman滤波算法实现雷暴追踪与外推。数据处理结果表明：

1) 基于闪电数据，本文提出的方法能有效识别雷暴，并对雷暴进行实时追踪。该方法能有效识别并追踪各类雷暴，同时也能有效识别雷暴的分裂与合并的情况。

2) Kalman滤波算法具有较好的0~60 min的临近外推预报能力，各项性能指标整体接近TITAN算法，对临近预报具备一定的参考价值。

3) 闪电往往出现在较强的对流系统中，对对流系统的发展、消散具有一定指示性意义。

闪电活动对于各类强对流天气具有不同活动特征，研究上述变化特征，进而实现利用闪电变化特征实现强对流天气的临近预报是目前的重要研究方向。未来将以本文工作为基础，以单个雷暴单体为载体，进一步结合闪电变化特征，在雷暴外推的基础上，实现分类强对流天气的预警预报，进一步挖掘闪电数据的使用价值。