基于卡尔曼滤波和TITAN路径的雷电预警

胡 鹏,伍瑞林,伍光胜,张志坚

(1.广州市突发事件预警信息发布中心，广州 511430; 2.广州市荔湾区气象局，广州 510140)

0 引言

雷电是一种常见的天气现象，常伴随着强对流等天气过程出现，容易造成破坏性极大的雷灾损失。广东省地势北高南低，珠三角三面环山，南朝南海，呈一个半封闭低谷，山地的向阳坡和迎风坡有利于海洋水汽抬升，冷空气来临时容易出现锋面雷暴天气，尤其广州市地处珠江三角洲，城市范围大，人口分布稠密，是广东省雷击风险高值区之一[1]。经统计，2019年广州市共发生云对地闪电约23万次，平均地闪回击密度值为30.64次/(平方公里·年)，平均地闪强度为10.44 kA，发生雷灾事故37起，占全省雷灾总数的14.5%，较2018年增加48%，因雷电灾害造成的经济损失215.41万元，其中电力核电行业是雷灾的重灾区，雷灾起数最多，造成的经济损失最大。因此，加强对雷电灾害产生的分析和研究，尤其是对雷暴运动发展过程的研究，从而实现高精度的雷电预警，对进一步完善城市雷电灾害防御等方面有着重要意义。

目前，国内外均主要利用闪电定位资料、大气电场资料和雷达回波资料开展雷电预警相关技术研究。如美国、加拿大、巴西、法国、英国、日本等均建立了闪电监测定位网，其中美国于20世纪七八十年代就已开始建设雷电预警系统[2]。在我国，中国气象科学研究院近几年也开发出了综合利用雷达、卫星、闪电定位、大气电场及探空等多种资料的雷电临近预警系统[3]，并投入业务运行；梁巧倩、林良勋根据雷电发生的天气学分型和完全预报方法研究出了一种广州地区雷电短期预报方案，对雷电潜势预报具有较高参考价值[4]；王凯等通过利用黄山周边地区雷暴过程大气电场仪和闪电定位仪数据，选取预报因子建立了雷电预警的预报方程，获得较好的应用效果[5]；路郁等对闪电定位数据进行聚类分析，利用“膨胀-腐蚀”算法还原雷暴云，并外推未来雷电落区[6]。以上雷电预警预报研究都取得了一定进展，但雷电发生时具有一定的随机性、瞬间性、地域相关性等特征，增加了其规律研究的难度。近年，随着计算机技术和模式识别技术的快速发展，又出现了利用闪电定位和雷达回波数据识别、追踪雷暴云，从而采用外推法进行雷暴运动趋势临近预报的方法，并取得了一定成果[7]。如R.E.Rinchert提出利用模式识别方法识别和匹配相邻强回波单体；J.T.Johnson提出改进的SCIT算法[8]，使雷暴跟踪的准确率提高到90%；文献[9]提出了利用闪电定位数据基于DBScan聚类算法进行雷暴云识别从而进行其运动趋势预测的方法，可实现预测未来15 min内雷暴位置和覆盖区域预测。

本文提出了一种新型的基于clique聚类识别[10]和卡尔曼滤波[11-12]进行雷电识别及外推的方法，并在雷电预警监测系统中分别利用本方法和传统的TITAN风暴模型[13]路径预测方法实现了未来1小时逐6分钟的雷暴路径外推，基于2020年广州闪电定位数据对两种方法进行了有效性检验评估。

1 雷电识别外推算法设计

1.1 基于卡尔曼滤波的追踪外推算法

卡尔曼滤波是一种最优化自回归数据处理算法，是一种针对线性系统的高效率状态滤波技术，在对时间序列分析、轨迹最优化等方面都有较好的效果[14-15]。它是一种递归的滤波过程，即只要获知上一时刻状态的估计值以及当前状态的观测值就可计算出当前状态的估计值，因此不需要记录观测或者估计的历史信息(如图1)。

图1 卡尔曼滤波过程

在雷暴预测应用中，可以首先利用Clique空间网格聚类算法对雷暴单体进行识别，然后利用卡尔曼滤波算法，对识别出的雷暴单体运动进行线性外推。其中每6分钟设为一个时间段，在时间段结束时，采集融合该6分钟内的闪电数据，按图2算法流程执行。

图2 卡尔曼滤波路径外推流程图

1.1.1 基于Clique聚类算法的雷暴识别

由于闪电基本上由雷暴云团生成，在空间上具有连续性，雷暴云团的移动意味着闪电位置的移动，因此在进行雷电外推时，可通过聚类算法，将在空间上具有一定连续性的离散闪电聚合成独立的雷电单体，从而简化雷电外推过程。此处选择了运算速度较快，抗干扰性较好，可进行任意形状聚合的Clique网格空间聚类算法[16-18]。其基本步骤如下。

步骤1：闪电坐标网格化，将6分钟内的所有闪电根据经纬度坐标放入网格化的平面直角坐标系；

步骤2：根据设定好的阈值参数，扫描符合闪电频数阈值和雷暴面积阈值的网格并合并；

步骤3：记录符合阈值要求的网格为闪电簇区域；

步骤4：合并相邻的闪电簇区域并编号；

步骤5：扫描所有闪电簇区域，利用椭圆法包络闪电簇，返回识别的雷暴区域。

1.1.2 卡尔曼预测

把雷暴的移动过程看作是系统的变化，应用卡尔曼滤波中的时间更新方程对前一时间段的雷暴进行预测，从而得到现时间段的预测值。通过前一时间的系统状态Xn-1和协方差矩阵Pn-1预测最新时间的系统状态和协方差，并把预测结果作为预测值保存下来。该过程使用时间更新方程：

(1)

(2)

其中:F是状态转移矩阵，代表系统变化的方向和幅度，Q是状态转移噪声矩阵，代表该过程可能产生的误差变量。

1.1.3 雷暴追踪

步骤1：获取在前一时间段的N1个雷暴，在雷暴识别步骤中获取现时间段识别的N2个雷暴，把两个相邻时间段的雷暴作为二分图中的两个集合A和B；

步骤2：计算A中每个雷暴i对应B中每个雷暴j的代价函数Cij，其中0

Cij=dp+ds

(3)

其中:dp是雷暴中心点的位置差；ds是雷暴面积的平方根的差。将Cij作为雷暴二分图之间边上的权值；

步骤3：针对带权值的雷暴二分图，利用二分图最小权值匹配寻找到能使代价函数的和∑Cij的匹配方案；

步骤4：完成匹配后，把雷暴分为四个集合，分别为前一时间段的已匹配和未匹配雷暴集合Mn-1和UMm-1，以及现时间段的已匹配和未匹配雷暴集合Mn和UMn；

步骤5：结合卡尔曼预测所得的预测值，再分别在重合的基础上匹配Mn-1和UMn、UMm-1和Mn，以此判断是否产生了雷暴分裂和合并，并对分裂和合并的雷暴进行方向和速度上的校准；

步骤6：分析处理剩余雷暴，添加新生雷暴，舍弃追踪失败的雷暴，并把雷暴追踪的结果作为观察值。

1.1.4 卡尔曼更新

应用卡尔曼滤波中的状态更新方程对雷暴状态进行更新，从而得到现时间段的雷暴状态值。首先利用预测过程得到的协方差预测值结合观察矩阵计算卡尔曼增益Kn，再利用系统的预测值Xn-1和最新时间的观察值Yn更新系统的最新状态。该过程使用状态更新方程：

(4)

(5)

(6)

其中:Xn是更新后的后验状态矩阵，Pn是更新后的后验协方差矩阵，H为观察模型，R是观察噪声矩阵。

1.1.5 外推预测路径

根据更新后得到的雷暴状态，通过线性外推得到雷暴预测路径，以此外推未来1小时逐6分钟(即6、12、……、60 min)的雷暴状态，得到一条雷暴外推路径；在每次时间段结束都可以进行迭代外推，从而不断更新且精确化外推路径。

1.2 基于TITAN风暴路径的追踪外推算法

基于TITAN风暴路径进行雷电追踪和外推则是目前应用较多的另一种雷电临近预报方法。其主要过程也包括：闪电融合、聚类识别、追踪外推[19-20]。

1.2.1 闪电融合

由于TITAN风暴移动路径时间间隔和雷达回波探测间隔均为6 min，首先需要将每次采集到的零散闪电定位数据融合成逐6 min过程，与雷达的探测间隔相对应，融合过程如图3所示。

图3 闪电融合过程图

当前时间采集到闪电，计算出起始和截止时间，将其对应06/12/18/24/30/36/42/48/54/00分钟时段内，按照时段的起止时间，从数据表中检索出已存在的闪电进行融合，从而形成一个逐6 min的闪电集合，作为后续进行聚类识别和外推的雷电数据。

1.2.2 聚类识别

此处采用DBSCAN (density-based spatial clustering of applications with noise)密度聚类算法，用距离作为闪电密度的聚合因素，将多个离散的闪电聚合成若干个雷电单体，离散的闪电聚合成雷电单体后，利用雷电单体内的闪电位置，加权计算出雷电单体的质心位置，再以质心位置与风暴移动路径进行空间关联和外推。

1.2.3 外推预测

利用成熟的TITAN雷暴识别、追踪、分析和临近预报算法，从前后两个时次的雷达回波图中识别出雷暴单体，利用两个单体之间的位置偏移，计算出雷暴的移动方向和移动速度，从而得到风暴的移动路径，将其作为雷电的移动路径。由于TITAN数据接口提供的是未来1小时逐10 min的风暴移动路径，还需要进行差值运算生成未来逐6 min的风暴移动路径以匹配雷电预报路径。

通常雷电单体会随着风暴移动，但雷电单体和风暴位置往往并不重叠，同时由于闪电定位仪和雷达探测设备之间数据处理存在时间差，因此需利用时间和空间关联方法，根据聚类出来的雷电单体的时间和质心位置，匹配到最近时间和空间距离最近的风暴移动路径，将风暴预报路径平移后，作为雷电单体的预报路径，对雷电单体进行逐6 min的位置外推(图4中虚线部分)，生成未来1小时逐6 min的雷电预报产品。

图4 移动路径外推图

2 雷电预警模型和平台设计

基于上述两种雷电识别外推方法和广州范围内闪电定位仪、大气电场仪、雷达等雷电实况监测数据，广州市气象局开发了雷电监测预警系统，可在GIS平台上有效提供当前雷电实况数据和未来1小时雷电预报产品，并对两种识别外推方法进行有效性检验评估，可广泛应用于防雷减灾、雷灾调查、预警预报服务等领域。其整体方案如图5所示。

图5 系统结构图

1)雷电实况监测资料：包括EN全闪、粤港澳闪电、风暴路径和雷达回波及广州新建大气电场仪等。

2)雷电预警模型：包括采集、融合、聚类分析、外推、预警和网格化等过程。①采集：从广东省气象探测数据中心的通用气象数据访问接口采集雷达回波、EN全闪、粤港澳闪电、风暴路径等雷电实况资料，入库存储；②融合：将逐分钟的闪电定位融合成逐6 min的数据，与雷达观测时间进行匹配；③聚类分析：将离散的闪电定位聚合成雷电单体；④外推：利用风暴移动路径，对雷电单体进行未来1小时的外推；⑤预警和网格化：将地图划分成1 * 1 km分辨率的网格，设置每个格点的雷电红色和黄色预警标准，利用未来1小时逐6 min的雷电单体位置，按照距离权重法，插值成未来1小时逐6分钟的雷电精细化格点预警产品。

3)雷电数据库：采用关系型数据库和NetCDF文件库相结合的方式来存储雷电实况监测和预警数据。其中关系型数据库用于保存闪电和风暴路径等结构化雷电数据；NetCDF文件库用于存储雷达回波和雷电精细化格点预警产品等网格数据。

4)雷电预警模型检验系统：在WebGIS地理信息系统上，实现雷达回波、闪电、风暴移动路径和雷电精细化格点预警产品的综合显示，通过对雷电数据的可视化分析，检验雷电预警模型的输出结果。

系统界面如图6，其中黑色路径为基于TITAN风暴路径生成的未来1小时雷电路径预报，灰色路径为基于卡尔曼滤波算法生成的未来1小时雷电路径预报。

图6 系统界面

3 实验结果与分析

3.1 数据来源

为了对两种外推算法结果的有效性进行检验评估，本文使用了2020年5月1日至2020年10月27日广东省气象局数据接口所提供的广东EN全闪观测网的全省闪电定位数据作为检验结果的数据来源。

3.2 检验方法

定义命中率和提前量为两个检验指标。

其中命中率是针对每个预报路径点，读取对应时刻发生的闪电，如闪电定位数据显示半径5 km范围内有地闪，则认为命中，否则认为空报，即：

(7)

其中:N命中为命中次数，N空报为空报次数。

提前量则是针对每个预报路径点，其对应时刻的未来1小时内发生闪电，则认为命中，计算第一次闪电和当前预报之间的时间差，作为预报的提前量。

统计时，以间隔6分钟为一个时间片，从起始时间至结束时间整个时间段内逐时间片计算每个预测路径的命中率和提前量。

3.3 结果分析

统计后的检验结果如图7可见，两种外推算法的命中率均随预测时间的增加而衰减，越是临近当前时刻，预测命中率越高，未来6分钟预测命中率最高可达到87.5%，未来60分钟预测命中率最高27.2%。且由于2020年8月份对基于卡尔曼滤波的外推算法进行了参数优化，将聚类识别过程中的网格密度增大后，其预测命中率在2020年9月至10月的检验中已明显优于基于TITAN路径的外推算法。而两种外推算法在各时次的预测提前量上差别不大，平均相差1分钟左右，整体表现比较接近。但在命中次数上，基于TITAN路径的外推算法要明显高于基于卡尔曼滤波的外推算法，其主要原因还是因为后者所采用的Clique聚类算法对分散的较小型雷暴单体识别率低于前者采用的DBSCAN密度聚类算法，需进一步进行优化改进。

图7 检验结果对比

4 结束语

本文提出了一种基于clique聚类识别和卡尔曼滤波算法进行雷电识别及外推的方法，并利用该方法和传统的基于TITAN风暴路径外推的两种算法，研究开发了广州雷电监测预警系统，该系统接入实时雷电气象数据，可生成未来1小时逐6分钟1*1 km分辨率的雷电精细化格点预警产品，并对预报产品进行检验，可为雷电预警业务化提供有效平台。实际数据检验表明，两种算法均能有效识别追踪雷暴并进行未来1小时路径预测，经过参数优化后，基于卡尔曼滤波的算法在预测命中率上已优于传统的基于TITAN风暴路径的外推算法，下一步还需要在聚类识别算法上进行研究优化，进一步提高雷电预警的准确性、及时性。