基于时序栅格的暴雨事件提取与追踪方法

杨光辉， 薛存金， 刘敬一， 邬群勇， 伍程斌

1.福州大学卫星空间信息技术国家地方联合工程中心，福州350116

2.中国科学院数字地球重点实验室，北京100094

3.中国科学院空天信息研究院，北京100094

综合对地观测和众源技术为获取长时序、高时空分辨率的降水时序栅格数据集提供了技术支撑[1].高时空分辨率的暴雨事件集可以为获取暴雨影响区域的时空分布、挖掘暴雨事件的时空模式、开展面向个体的精准暴雨灾害预警预报服务提供重要数据基础[2].因此，基于时序栅格数据集并结合暴雨事件的时空认知方法提出暴雨事件提取方法，进而构建高时空分辨率的暴雨事件集，这对暴雨灾害及城市内涝的监测和预测具有重要意义[3].

目前，暴雨事件的提取方法主要分为3 类：1）基于站点日降水数据，即根据气象上暴雨定义直接提取阈值，如文献[4]以站点日降水数据提取暴雨天数，分析了中国近40年暴雨发生频率的年代际时空变化特征；文献[5]根据站点日降水数据分析了暴雨过程和暴雨日的年、季、月气候分布及变化特征；文献[6]基于日降水数据统计年代际暴雨雨量、雨日、雨强来分析中国1951—2010年的年代际暴雨时空变化格局.2）对栅格数据进行日累加处理，设定特定阈值进行暴雨事件提取，如文献[7]对全球降水测量(global precitation measurement,GPM)与热带降雨测量任务(tropical rainfall measuring mission, TRMM)降雨数据进行日累加后提取暴雨栅格；文献[8]以日累加方式处理GPM 降雨数据后进行适用性评价.3）基于雷达反射率阈值开展时刻状态的暴雨事件提取，并利用暴雨的移动特征实现暴雨事件的追踪，如文献[9]提出用暴雨识别、跟踪、分析和预报(thumderstorm identification, tracking,analysis, and nowcasting, TITAN)方法将空间上反射率超过一定阈值的连续区域作为暴雨区域进行识别和匹配追踪；文献[10]基于双阈值迟滞法识别暴雨，避免了只有几个像素的暴雨区域；文献[11]基于TITAN 在识别暴雨过程中使用多阈值分割技术（hierarchical threshold segmentation, HTS），减少了暴雨过程识别中错误的暴雨对象合并情况.

上述方法在暴雨提取分析方面取得了一系列成果.基于站点的降水数据在局部尺度上（站点位置）具有较高的观测精度，但在偏远地区或地形复杂地区，因为观测站点分布稀疏，所以很难提供精度较高的降水空间分布特征[12].将时间分辨率高的栅格降水数据先进行日累加再进行提取分析则会损失原有数据精度，以致无法提取整个暴雨事件的发展细节.基于TITAN 的暴雨事件识别与追踪算法利用了暴雨事件的移动特征，为暴雨事件的提取与重构提供了思路，但在追踪过程中只考虑了暴雨路径起止位置的分裂合并，因而无法识别暴雨路径中间分裂合并行为特征.针对以上问题，本文设计了一种基于时序栅格的暴雨事件提取与追踪方法，实现了对暴雨事件的完整提取.

1 暴雨事件提取与追踪方法

暴雨事件是指降水量超过一定阈值且在时空中连续变化的降雨事件.暴雨事件提取与追踪方法（rainstorm event extracting and tracking method, REETM）基于长时间序列栅格数据集，依据暴雨在时空上连续渐变的过程特性[13-15]开展暴雨事件提取与追踪.REETM 主要包括以下3 个步骤：

步骤1在时间维度上提取暴雨栅格，依据降水数据划分暴雨等级；

步骤2在空间维度上提取暴雨对象，连接空间上孤立的暴雨栅格并统计暴雨信息；

步骤3暴雨事件时空追踪，在时空维度上追踪暴雨并重构暴雨事件.

1.1 时间维度暴雨提取

基于长时间序列降水量栅格数据集进行时间维度的暴雨提取处理，提取出暴雨级别信息，获取栅格暴雨分级数据集，如图1 所示.

图1 时间维度上的暴雨提取Figure 1 Rainstorm extraction in a temporal domain

具体步骤如下：

步骤1取出每个栅格像元所在位置的所有图像栅格值，形成一个按照时间顺序排列的一维数组.

步骤2遍历这个一维数组，寻找极大值点，如果数组任意两个不为0 的位置时间差值小于阈值t，则认为降水连续.以极大值点为中心，向前、向后搜寻降水连续的序列.

步骤3统计降水连续的序列范围内的累计降水量，如果序列时间长度超过24 h，则以极大值点前12 h 这一时刻为统计起始点，以极大值点后12 h 这一时刻为统计结束点.依据中国气象局制定的《降水量等级》[16]对该暴雨序列进行分级，分级依据如表1 所示.

表1 暴雨等级划分Table 1 Rainstorm classification

1.2 空间维度暴雨提取

在时间维度暴雨提取结果的基础上进行空间维度暴雨提取.在时间维度上的暴雨提取保证了暴雨事件在时间维度的连续性；在空间维度上的暴雨提取不但保证了暴雨事件在空间维度的相邻性，而且能提取出暴雨对象及其信息，如图2 所示.

图2 空间维度上的暴雨提取Figure 2 Rainstorm connection in a spatial domain

具体步骤如下：

步骤1采用结点搜索法[17]连接孤立的暴雨栅格，将栅格暴雨分级数据集转换为矢量暴雨数据集.

步骤2遍历矢量暴雨数据集中每一个内部连通矢量多边形即暴雨对象，从原始栅格数据集和栅格暴雨分级数据集中寻找出该暴雨对象覆盖范围内的所有栅格像元.

步骤3统计暴雨对象覆盖范围内所有栅格像元信息以获取暴雨对象信息，包括时间信息、空间信息、属性信息.时间信息即该暴雨对象的发生时刻，空间信息包括暴雨对象空间范围、面积和周长，属性信息包括暴雨强度、降水体积、平均降水量、最大降水量、最小降水量.

步骤4保存所有暴雨对象，得到完整的矢量暴雨数据集.

1.3 暴雨事件时空追踪

暴雨事件时空追踪的核心思想如下：一个暴雨事件在时空中是连续变化的，它由多个时刻状态组成，且每个时刻状态包含一个或多个暴雨对象，于是可以通过追踪相邻时刻空间范围拓扑相交的暴雨对象来实现暴雨事件的追踪.

暴雨事件时空追踪的流程与步骤如下：同一暴雨事件在短时间内移动距离较短，空间范围存在重叠.因此，对于相邻时刻暴雨对象，若其空间范围拓扑相交，则认为暴雨对象属于同一个暴雨事件，如图3 所示.

图3 暴雨事件时空追踪Figure 3 Rainstorm tracking in spatiotemporal domain

在图3 中，绿色多边形表示暴雨对象，T1时刻暴雨对象1 与T2时刻暴雨对象3 在空间上拓扑相交，则认为暴雨对象1 与暴雨对象3 存在关系1→3，即暴雨对象1 发展为暴雨对象3，且暴雨对象1 与暴雨对象3 属于同一个暴雨事件；同理找出所有时刻属于同一个暴雨事件的暴雨对象，并统计暴雨事件信息.具体算法步骤如下：

步骤1读取矢量暴雨数据集，取出每一个暴雨对象并将其存储到暴雨对象集合P 中.

步骤2从第1 个时刻开始追踪暴雨对象，令t=1.

步骤3遍历暴雨对象集合P 中t 时刻每一个暴雨对象pi，如果集合P 中t+1 时刻存在暴雨对象pj，使得pi的空间范围与pj的空间范围拓扑相交，则记录暴雨对象关系pi→pj，并将该暴雨对象关系存储到暴雨对象关系集合R 中；若pi与t+1 时刻多个暴雨对象空间范围拓扑相交，则记录pi与每个相交暴雨对象的关系并存储到集合R 中.

步骤4令t 自增1，如果t 为最后时刻，则执行步骤5；否则执行步骤3.

步骤5若P 中不存在暴雨对象，则执行步骤8；否则创建暴雨对象集合S 以存储属于同一个暴雨事件的所有暴雨对象，创建暴雨对象关系集合R'以存储该暴雨事件中所有暴雨对象关系.创建暴雨对象集合P'和P''，从暴雨对象集合P 中取出一个暴雨对象存储到P'中.

步骤6遍历P'中每一个暴雨对象p，若R 中存在关系p →P'，则从R 中取出该暴雨关系储存到暴雨对象关系集合R'中，并从P 中取出暴雨对象P'存储到P''中.

步骤7将P'中所有暴雨对象取出后存储到S 中.若P''中存在暴雨对象，则将P''中所有暴雨对象取出储存到P'中，执行步骤6；否则保存暴雨对象S 中的所有暴雨对象并统计该暴雨事件信息（暴雨级别、平均降水量、降水体积等）后予以保存.保存属于该暴雨事件的对象关系集合R'中的所有暴雨对象关系，执行步骤5.

步骤8所有暴雨事件提取完毕，则算法结束.

在以上算法中需要特别指出的是：暴雨事件是一个在时空上变化的实体，其平均降水量尚未有定义，这里暂且采用的计算公式为

式中，Ravg为暴雨事件平均降水量，n 为暴雨事件中的暴雨对象数目，Ai为第i 个暴雨对象的面积，Ri为第i 个暴雨对象的平均降水量.

2 实例验证

本文以GPM IMERG Final Products（简称为GPM 产品数据集）作为实例数据开展暴雨事件提取工作.GPM 产品数据集来源于NASA（disc.gsfc.nasa.gov），其空间分辨率为0.1◦，时间分辨率为0.5 h；验证数据集包括暴雨预警信息（www.weather.com.cn）、49 个气象台站的日降水数据集和2016-07-01 的天气雷达数据集（data.cma.cn）；研究区域为70◦E∼140◦E,0◦∼60◦N，时间范围为2016-06-01∼2016-09-30.

将时间维度暴雨提取过程中最大降水连续间隔时间t 设置为3 h，提取出49 003 个暴雨事件，其中暴雨级别事件有48 937 个，大暴雨级别事件有66 个，特大暴雨级别事件有0 个.

2.1 与日累加法的对比分析

采用探测率（probability of detection, POD）nPOD统计指标评价卫星降水数据对暴雨事件的探测能力，该指标表示探测正确的数目占实际数目的比例，其值越高表示暴雨事件提取方法的漏报率越低.表2 对比了日累加法与本文提出的REETM 方法的暴雨探测率，其中D0为暴雨真实天数，D1为采用日累加方法探测得到的暴雨天数，D2为采用本文方法探测得到的暴雨天数.通过计算得出本文方法的探测率平均比日累加法高12.97%，原因是日累加法只统计每日固定起止时间累计降水量来识别暴雨，而本文方法统计累计降水量的起止时间是动态变化的，可以探测出固定时间段内探测不到的暴雨，所以探测率较高.

表2 日累加法与REETM 的暴雨探测率对比Table 2 Comparison of rainstorm POD between daily accumulation method and REETM

2.2 与暴雨预警信息的对比分析

选取关键时刻展示一个暴雨事件的演变过程如图4 所示，紫色多边形为本方法提取出的暴雨事件覆盖区域.该暴雨事件的起止时间是2016-07-18∼2016-07-21，起源于湖南和陕西，经过发展合并分裂等过程影响全国多个地区，由东北地区离开中国国境并逐渐消亡.2016-07-18 8:00 陕西出现暴雨，由西向东移动并逐渐扩大，途径陕西北部、山西大部分区域、河北大部分区域与河南北部.同时，湖南区域暴雨，由西南向东北移动并逐渐扩大，途径湖南北部、湖北中部和河南南部，于2016-07-20 0:00 在河南中部与起源于陕西的暴雨合并，然后继续向东北移动，于2016-07-20 2:00 分裂为两部分，一部分逐渐缩小，于2016-07-20 22:00 在河南境内消亡；另一部分继续向东北移动，途径河南、河北、山东、北京、天津和辽宁等地，于2016-07-21 13:00 逐渐离开中国国境并消亡.将本暴雨事件与多地气象台预警信息进行对比，可以发现本文提出的REETM 方法所提取的暴雨事件在发展过程中与气象台发布暴雨预警在时间和空间上基本吻合.

2.3 与雷达数据集的对比分析

图5 将本文方法提取的暴雨事件图像与雷达反射率图像进行对比，其中上侧4 幅图为本方法提取出的暴雨事件，下侧4 幅图为对应时刻雷达反射率图像.可以看出本文方法提取出的暴雨事件与雷达图像中基本反射率较高的区域基本吻合，故本文方法提取的暴雨事件位置与实际暴雨发生位置基本吻合.

3 结 论

针对长时间序列、高时空分辨率的降水量栅格数据集，本文从时间维度、空间维度和时空追踪这3 个层次设计了一种暴雨事件提取与追踪方法.时间维度的暴雨提取保证了暴雨在时间上的连续性，空间维度的暴雨提取保证了暴雨在空间维度的相邻性，暴雨时空追踪保证了暴雨在时空上的连续性和一致性.通过实验例证与传统提取方法的对比分析，验证了本文提取方法REETM 的正确性和可行性，并得出了以下3 个主要结论：

图4 2016-07-18∼2016-07-21 某暴雨事件Figure 4 A rainstorm event from July 18 to 21, 2016

图5 暴雨事件图像与雷达反射率图像Figure 5 Rainstorm events images and radar reflectivity images

1）本文方法从时间、空间、时空3 个维度上进行设计，实现了暴雨从产生发展到消亡事件的提取与追踪；

2）相较于日累加暴雨提取方法，本文方法的探测率平均提高12.97%；

3）依据暴雨预警信息和雷达图像验证了本文方法提取暴雨事件的有效性和可行性.

本文提取的暴雨事件清楚地体现了暴雨在整个发展过程中的空间分布以及动态变化过程，对暴雨的分析及预测有着重要的意义.比如，暴雨在某个时刻发生了分裂，可能由大气气流突然变化造成.然而，对于暴雨对象运动过快以及前后时刻空间覆盖范围没有重叠的情况，采用本文方法提取的准确度会受到影响.因此，今后的研究重点在于结合暴雨运动路径对暴雨位置进行预测，以完善优化算法.