0 引言

随着探测手段的不断丰富和加工处理能力的不断增强, 气象数据产生的频次和时增量已分别达到了秒级和TB级, 数据总量也呈指数级增长。传统气象数据主要是以站点或者格点为表现形式的结构化数据, 且大多采用Oracle、MySQL等关系型数据库, 以及常用于科学计算领域的NetCDF或HDF等文件数据形式存储, 并根据时间形成序列化数据进行存储和管理。

面对海量增长的历史气候资料, 无论是传统的关系型数据库还是文件系统, 数据的检索和定位性能均较差, 多维数据查询和分析的效率均较低, 目前已难于满足大规模、多样化的气象数据分析挖掘的性能需求。近年, NASA (National Aeronautics and Space Administration) 气候模拟中心利用MapReduce分布式计算框架开展了地球气候模拟数据的高性能分析查询研究^[1-2], 对比传统的NetCDF数据集查询有明显的计算性能提升。广东省气象局李永生等^[3]利用HBase实现了格点类数值天气预报数据的分布式存储和管理, 其访问效率比传统方法也有明显提升。然而, 在气象业务中将分布式存储技术与分布式计算有效结合的相关研究还未见报道。

本文利用近年来积累的海量多普勒天气雷达数据, 展开基于分布式计算框架的风暴三维追踪的方法研究, 进而研究气象业务背景下的分布式存储与计算技术相结合的方法, 以帮助业务技术人员能够更高效地实现局部地区强对流天气的历史气候特征分析。最后，以珠三角超大城市群在不同季节对于风暴的生消起到的促进与抑制作用为例, 验证了采用本文方法获得的分析结果的有效性。

1 三维风暴追踪方法 1.1 TITAN算法

TITAN (Thunderstorm Identification, Tracking, Analysis and Nowcasting) 算法是目前临近天气预报的主要方法之一, 主要用于对强对流天气的识别、跟踪和预报, 它使用组合最优算法对风暴体进行追踪, 使用几何算法来实现对风暴的合并、分裂的识别^[4], 本文以基于TITAN算法思想实现三维风暴的追踪, 并将其与分布式存储管理技术相结合实现高效计算。

1.2 三维风暴的识别

三维风暴的识别与追踪其本质可视为在立体网格数据中目标像素点的近似邻近像素点的搜索问题。其基本步骤包括以下几步:

1) 从给定的立体网格数据的最下层开始进行行扫描, 当遇到第一个回波强度大于预先设定的阈值像素点 (x, y, 0) 时, 则以该回波为起点, 搜索与该像素点相邻的其他相连像素, 并检查是否有其他相邻新像素点的回波强度达到阈值以上, 如果有则进行标记, 一方面实现覆盖率的提高, 另一方面可避免像素点的重复, 进而提高搜索效率。

2) 进行垂直搜索, 作适当的“膨胀”处理。当回波强度从阈值以上下降到阈值以下, 像素点仍然被标记为阈值以上, 并进行下一步的扩展搜索, 搜索过程中, 如未发现“连通”的其他网格点, 则停止继续膨胀。

3) 对所有搜索到的新像素点, 重复1) 中操作, 直到搜索不到任何新的达到阈值的像素点为止, 所得所有像素点的集合, 便是一个强度在阈值范围内的强回波区, 形成风暴单体。

4) 从像素 (x+1, y, 0) 开始, 重复过程1)~3), 直到立体数据扫描完毕完成扫描, 并对所有风暴单体进行标记。

经过以上过程, 即可实现立体网格数据中的风暴单体的确定, 进而可以基于确定的各风暴单体中包含的像素点数据计算获得风暴中心强度、风暴伸展高度、风暴单体厚度等其他风暴体相关的一系列属性信息。在该过程中, 可利用MapReduce方法将各风暴单体的搜索过程进行任务切分, 实现搜索过程的分布式, 以及属性信息计算的分布式。

1.3 三维风暴的追踪算法

雷暴、冰雹和强降水等恶劣天气在雷达强度图像上表现为强回波区, 即由彼此连通且回波强度大于某一阈值的像素点的集合。因此, 识别和追踪强烈的恶劣天气, 就是追踪大气三维立体网格数据中的强回波区域。本文将雷达回波反射率的阈值设定为35 dBZ, 三个垂直正交平面18连通, 跟踪搜索风暴体位置和边界 (见图 1)。

基于大气三维立体网格数据实现风暴的追踪算法如下所示。

输入:三维立体网格数据。

输出:风暴中心强度、风暴伸展高度、风暴单体体积、风暴单体厚度。

1) 雷达回波反射率的阈值=35 dBZ;

2) If (回波强度>阈值)

3) Then起点=该回波强度的像素 (x, y, 0);

4) 搜索该起点邻近相连的12个像素;

5) If (有新像素点的回波强度>阈值)

6) Then打上记号;

7) 垂直搜索时, 作适当的“膨胀”处理;

8) If (回波强度从阈值以上下降到阈值以下)

9) Then该点仍然标记为阈值以上;

10) 进行下一步的扩展搜索;

11) If (仍然没有发现“连通”的网格点)

12) Then放弃搜索;

13) 对所有搜索到的新像素点;

14) If (回波强度像素点 < 阈值)

15) Then停止搜索;

16) Else GOTO 2);

17) For (起点=像素 (x+1, y, 0), 起点in立体数据, 起点的下一个数据)

18) 对风暴单体进行标记;

2 基于大规模多普勒天气雷达数据的分布式实现 2.1 数据的分布式存储管理 2.1.1 数据存储管理模型设计

将三维多普勒天气雷达数据解析后, 利用HBase实现描述雷达数据属性的元数据和雷达回波数据实体的存储, 并以带时间戳的Row Key进行标记。考虑到三维风暴追踪算法和未来三维可视化的需求, 没有对数据垂直分层存储, 而将三维数据整体存储入HBase的列存储中, 通过带时间的Row Key管理长时间序列的雷达历史资料集。由于大多数时间大气是处在没有降水的状态下, 因此, 三维数据中存在数据稀疏问题。因此在存入HBase的列存储之前, 对数据将作GZIP压缩, 一方面可以有效地解决数据稀疏问题；另一方面可实现数据存储空间的降低, 进而间接提升网络数据访问的延迟。在服务接口设计中, 提供分层读取数据的接口, 满足二维平面数据访问需求。具体存储方式如图 2所示。

2.1.2 数据服务接口开发

针对雷达数据量多、访问频次高的问题, 本文设计了针对网络应用开发的基于文本流的Restful接口 (实现方式如图 3所示), 该接口用于对外提供数据服务, 主要包括数据查询、时间序列数据接入与分析处理、三维雷达反射率场数据的接入与分析处理、二维平面和三维切片数据的接入与分析处理。此外, 接口的数据标准遵循JSON格式, 并可以支持可视化软件Surfer格式文件下载的接口服务, 进而实现满足主流数据格式的具体应用服务的接入。

2.2 算法的分布式实现方法

通过实施算法的并行化可以显著提高计算效率, 可以从算法本身的并行化处理以及基于数据分割的分布式计算两个方面进行考虑。由于三维风暴追踪算法原理复杂度不高, 所需数据是时间 (Row Key, RADAR:YYYYMMDDHHmm) 相关的, 每次风暴单体搜索过程是独立的, 主要的计算量源自长时间序列海量数据带来的密集型任务。因此, 设计主要采取基于数据分割的分布式算法实现, 这也是MapReduce编程模型的典型应用。如图 4所示, 启动MapReduce计算过程时^[5], HBase资源管理按照Row Key向Map任务提供数据片段, Map通过识别Row key所对应的雷达元数据和回波反射率实体数据进行分布式计算。计算结果通过Reduce任务汇总, 可以回写HBase, 也可直接连接可视化组件或者进入统计分析流程。

3 实验与分析 3.1 实验环境

采用HDFS技术, 建立适合气象业务实际的大容量数据分布式存储系统, 基于MapReduce技术探讨建立分布式计算模型, 选取典型类型的气象资料在云计算环境下设计开发大容量数据业务应用程序。目前以数值预报产品和三维雷达拼图进行实验。基于MapReduce分布式计算框架如图 5所示。

本文的实验平台在广东省气象局数值预报产品分布式管理服务平台基础上进行拓展^[6]。平台的基础架构由分布式存储设施层、基于分布式架构的数据环境与计算平台层和数据应用服务的接口层组成。

基于Hadoop技术架构, 构建了分布式存储设施层^[7], 为各类气象数据提供可伸缩、高容错的弹性存储资源。服务器配置为:CPU是Intel Xeon CPU E5-2603@ 1.80 GHz×4, 内存3 658 MB, 硬盘2.5 TB, 共7台实验机器, 网络连接1 Gb/s。

分布式数据环境与计算平台层由数据处理模块、分布式数据存储和管理模块以及分布式计算引擎三部分组成, 其中数据处理模块对自定义格式的雷达基数据进行解码和三维格式转换, 并开发实现了将数据实时接入平台的功能; 分布式数据存储和管理模块对HDFS和HBase进行集成^[8], 实现了高时空分辨率气象雷达数据的高效管理和读写操作; 在分布式计算引擎模块中MapReduce通过把对数据集的大规模操作分发给网络上的每个节点进行处理^[9], 每个节点会周期性地把完成的工作和状态同步更新到管理节点上。

服务接口层包括分布式数据服务接口和分布式算法服务接口, 为了更好地与其他业务系统实现对接, 对外提供了基于Restful接口的标准化数据查询和访问接口^[10]。

3.2 实验数据

多普勒天气雷达通过发射电磁波, 对三维大气在不同仰角上作360°的立体扫描 (Volume Scan), 得到球坐标系的三维基数据 (Base Data), 其中的雷达回波反射率表征大气中降水的强度。设定不同的高度, 应用标准大气下的雷达测高公式, 获取该高度上的数据。如果该高度上没有数据, 则在最邻近的上下2个仰角之间利用双线性插值获得数据。通过设定高度值就可以获得等高平面上的雷达回波反射率值 (Constant Altitude Planar Position Indicator, CAPPI)。循环进行该算法, 就可以获得从近地层到对流层上部三维立体的雷达回波反射率。本文使用广州及周边地区多普勒天气雷达基数据进行三维立体CAPPI的网格化处理, 获得500×400×40层的雷达回波反射率产品, 约21 MB, 三维网格雷达回波反射率数据描述如表 1所示。多普勒天气雷达每6 min完成一次体扫描, 以2007~2009年3年周期作为实验分析, 原始数据文件262800个, 共约5.52 TB数据。

实验测试包括雷达数据的访问效率测试以及执行雷达三维风暴追踪算法的性能测试, 并分别与NetCDF数据接口和单机读取文件序列方式进行性能对比。

3.3 雷达数据访问效率

雷达数据的访问效率测试主要考察在并发条件下, 分别读取三维雷达资料场数据、等高的二维平面、任意方向的三维切片的性能效率, 并和传统的利用科学计算通用格式NetCDF^[11]的访问进行对比, 如图 6所示。利用通过分布式存储管理系统访问三维雷达数据, 一般控制在100 ms以内, 同时随着并发数的增长能够保持高效率。由于HBase采用了数据压缩, 因此在三维数据和数据切片读取效率上, 较NetCDF有明显优势。NetCDF原本并不支持并发访问, 随着用户数的增加, 出现访问锁死现象, 所以没有给出多用户访问数据。

3.4 算法性能测试

采用分布式MapReduce/Hadoop架构, 针对2007~2009年的雷达数据, 开展三维风暴追踪算法的3年回算, 分析珠江三角洲 (广州为中心方圆200 km为界) 发生的强对流天气 (雷达回波反射率>35) 的活动情况, 共费时193 min, 获取517时次风暴活动。此前, 采用单机读取26万个文件, 逐一计算的方式, 扣除期间中断时次, 根据程序日志, 共耗时219 h。

3.5 案例分析

利用风暴活动位置数据, 绘制不同月份的热点图, 如图 6。前汛期 (4月为代表), 广州及周边地区, 强对流风暴单体多发区在城市范围, 而周边则是风暴活动的低值区, 特别是广州市区一带, 5 km之内的网格一个月平均有30~45个风暴单体的活动。由于风暴的移动特性, 这种分布也往往表示着广州市城市群在前汛期 (湿季) 对于风暴的发生发展有促进作用。

后汛期 (7月为代表), 风暴活动较前汛期更加频繁, 但活动的主力区域是珠江三角洲的沿海地区, 这也反映了后汛期风暴单体主要源于主要是副高南侧边缘的对流活动以及东风波动系统西传所致。从图 7所示的广州市区的风暴活动分布图上看, 广州的西边 (芳村一带) 是风暴活动的活跃区, 但进入市区则风暴活动明显减少, 而东边移出广州市区, 风暴又重新活跃。这也反映了在后汛期 (相对的汛期干季), 大城市效应对于风暴的活动是有抑制和减弱作用的。这也和国外的一些研究以及开展的外场实验结果存在相似之处。

4 结语

本文基于分布式计算架构搭建了海量气象数据存储管理和算法应用平台, 实现了TB级别的雷达三维格点数据的分布式存储管理, 具有较好的扩展性; 采用MapReduce计算框架设计和实现了风暴追踪算法的分布式处理, 通过Restful API接口可以实现对数据的高效访问和算法的便捷调用。通过业务应用测试和典型场景的应用, 可以看到基于分布式架构的大数据分析应用模型, 对于未来高分卫星的精细遥感资料的实时和长时间序列的应用挖掘带来助力, 为业务人员更好分析局地天气气候事件和环境生态演变提供更强的计算能力和效率工具。