GeoCube：面向大规模分析的多源对地观测时空立方体

高凡，乐鹏，2，姜良存，曹志鹏，梁哲恒，上官博屹，胡磊，赵帅锋

1.武汉大学遥感信息工程学院,武汉 430079;

2.湖北省空间信息智能处理工程技术研究中心,武汉 430079;

3.广东南方数码科技股份有限公司,广州 510665

1 引言

随着对地观测体系的建立，遥感大数据不断累积（张立福等，2017；廖小罕，2021；Wagemann等，2018；Baumann 等，2018）。据统计目前全球在轨卫星达172 颗，包括已持续监测地球近50年的Landsat 卫星系列，欧洲航天局哥白尼项目下的哨兵卫星系列以及中国自主研发的高分卫星系列等（CEOS，2020a）。按照预估，仅哨兵系列的数据体量在2022年就将达到30 PB （Wagemann 等，2018）。然而，受限于对地观测领域传统数据处理模式，这些遥感大数据的价值难以得到充分发挥。其中涉及到若干关键技术问题包括：（1）海量对地观测数据的管理，包括数据监管、发现以及访问；（2）海量对地观测数据的处理，包括在线处理以及高性能计算；（3）海量对地观测数据的应用分析，包括传统分析以及基于深度学习的智能分析（李德仁等，2014；Sudmanns等，2020；Yue等，2016）。

20世纪90年代到21世纪初期间，对地观测领域发生了一次从面向本地分析到面向产品定制的服务模式变革。面向本地分析的模式中，用户需从数据门户下载数据，利用本地资源进行处理和分析，该模式受限于网络传输带宽、本地存储空间和计算能力，且不同门户的数据组织模式与数据存储格式多样，不利于数据发现和互操作（Wulder和Coops，2014；Kopp等，2019）。随着网络技术的发展，衍生了面向产品定制的服务模式，该模式基于网络空间服务标准，提供统一的数据发现、访问和处理接口，将数据处理迁移到远端服务器上执行，解决了互操作以及本地资源受限的问题（Yue 等，2015）。其中地理信息标准联盟OGC（Open Geospatial Consortium）提供了一系列领域认可的服务标准，比如网络覆盖服务、网络处理服务、网络地图服务等。虽然面向产品定制的服务模式通过计算迁移缓解了本地计算的压力，但该模式仍未解决海量数据的管理和分析问题。近年来，随着网络基础设施的快速发展，对地观测领域正经历着二次模式变革，从面向产品定制转变为面向大规模分析的服务模式。数据立方基础设施作为新一代模式中的核心技术，打破了传统模式技术屏障，其基于高性能存储/计算设施，采用统一的时空基准，将对地观测数据组织为分析就绪型数据ARD（Analysis Ready Data），向用户提供海量对地观测数据发现、访问和大规模分析的功能。其中，分析就绪型数据可定义为按照最低要求进行预处理并组织为可立即用于分析的卫星数据（Lewis 等，2018）。目前，对地观测数据立方体定义局限于遥感影像数据，本文在此基础上对立方体中的数据源进一步扩充，给出多源对地观测时空立方体定义：时间序列上的多维对地观测产品，其在空间上对齐并组织为适用于大规模计算的分析就绪数据，数据源可包括遥感影像、矢量数据、轨迹数据以及点云数据等。

随着数据立方体模型在对地观测领域掀起的热潮，国际上陆续开展了大量相关研究，比如地球服务数据立方体（EarthServer Data Cube）、开源数据立方体ODC（Open Data Cube）、地球系统数据立方体ESDC（Earth System Data Cube）（Baumann等，2018；CEOS，2020b；Mahecha 等，2020）。表1分别从数据源、组织模式以及计算策略等方面对比了地球服务数据立方体、开源数据立方体、地球系统数据立方体以及本文所提出的面向大规模分析的多源对地观测时空立方体GeoCube。EarthServer数据立方体数据源主要为海洋对地观测数据和气候数据，支持洲际尺度的数据管理与分析。EarthServer 数据立方体底层基于Rasdaman 阵列数据库系统（Baumann 等，1999），采用OGC 地理标记语言覆盖标准统一管理元数据。用户可通过OGC 网络覆盖服务和网络覆盖处理服务访问和处理立方体中数据，服务请求会被解析为Rasdaman 查询语言并执行。EarthServer 支持丰富的语言接口包括C++、Java、R以及JavaScript。

表1 不同立方体基础设施对比Table 1 Comparison of different cube infrastructures

开源数据立方体是由卫星对地观测委员会CEOS（Committee on Earth Observation Satellites）推动的项目，目标是为了降低数据立方构建过程中的技术屏障，提供一种构建数据立方体的开源解决方案（CEOS，2020b）。开源数据立方面向卫星遥感数据，其通过对遥感影像进行重采样、重投影等处理，以瓦片的组织模式将数据纳入到统一的时空基准下进行管理。开源数据立方采用5个简单的数据表并结合JSONB 格式来管理瓦片元数据，采用NetCDF 格式存储瓦片数据。开源数据立方同样提供了OGC 网络覆盖服务和网络处理服务等接口，并且支持用户在Jupyter Notebook 中自定义开发。

地球系统数据立方体是一个面向多源栅格数据的基础设施，该立方体将数据源作为一个维度，将多源栅格数据组织在一个统一分辨率和坐标系统的立方体中，有助于多源栅格数据的联合分析。该立方体主要采用Julia 语言进行开发，同时也支持Python和R语言。数据存储方面采用Zarr数据格式在云上存储（https://zarr.readthedocs.io［2020-12-15］），通过文件目录的形式组织。针对每一个立方体实例有一个配置文件，描述了立方体的时空分辨率等信息。用户可通过文件传输FTP（File Transfer Protocol）的方式浏览目录并访问数据，也可通过OGC 网络覆盖服务来获取数据，或通过JupyterNotebook访问并分析数据。

本文提出了一个面向多源数据管理和大规模分析的立方体基础设施GeoCube，相比于已有的立方体方法：（1）GeoCube 支持统一时空基准下的多源数据管理，包括栅格数据和矢量数据。基于商业智能领域事实星座模型管理元数据，可灵活扩展支持更多数据源；（2）GeoCube 结合瓦片组织模式和云计算优势，支持多源数据的长时序大规模分析；（3）GeoCube 引入人工智能领域技术，优化负载均衡性能，提高并行计算效率。

2 方法和设计

本文提出了一个新的地理空间数据立方体。首先，该立方体兼容多源地理空间数据包括栅格和矢量。针对数据特点，为立方体设计了4个维度包括空间、时间、产品以及波段。立方体中数据以空间对齐的方式组织，形成时间序列产品，因而该立方体又称为多源时空立方体。其次，立方体采用云计算技术来支持瓦片的高效处理，提出了一套面向立方体数据的分布式内存对象。目前支持矢量和栅格数据，可扩展支持轨迹等其他数据类型。最后，本文采用人工智能机器学习从数据科学的角度捕捉瓦片的计算特征，优化大规模并行计算过程的负载均衡性能，提高硬件的利用率。

2.1 立方体维度设计

2.1.1 立方体多维组织模型

对地观测数据通常由多维地理空间数据组成，维度包括时空维度以及其他主题维度。本文针对栅格和矢量数据特点，为立方体设计了4个维度如下（图1）：

图1 立方体模型Fig.1 The cube model

（1）产品维：数据的产品信息，包含产品名称、产品类型、所属传感器、层级信息以及参考坐标系。其中产品名称描述了某一类产品，比如Landsat 8 分析就绪数据和OSM 水体数据。产品类型用于判断产品属于栅格还是矢量数据类型。产品维和波段维通过一个产品&波段附加表关联，可获取该产品包含的波段信息。矢量产品没有波段信息，因此该附加表中只有影像产品。

（2）空间维：每个网格的空间信息，包含网格编码、行列号、层级、所属城市和省份。其中网格编码采用Z空间填充曲线以提高查询速度；层级属性存储了全球网格基准信息，支持用于分析的瓦片分辨率信息存储，比如1°×1°瓦片包含6250×6250个像素，同时支持用于可视化的金字塔瓦片的不同层级分辨率信息存储；城市和省份属性可满足用户在空间维度上的语义查询。

（3）时间维：瓦片的时间维度信息，包含成像时间，处理时间和年、月、日等信息。其中年、月、日的层级设计可以方便用户在时间维度上做聚集操作。

（4）波段维：记录不同传感器的所有波段信息，以及高级别产品的自定义波段信息，比如NDWI 产品的波段可定义为NDWI 波段。并且针对SAR 影像提供了波段和极化方式的组合。矢量数据没有波段信息，因此不和波段维关联。

（5）栅格瓦片事实：由度量信息和4个维度成员包括产品、空间、时间、波段组成，这里的度量信息包括瓦片ID 和栅格瓦片计算强度。瓦片ID作为瓦片的数据地址，用于访问栅格瓦片数据表中存储的瓦片数据。瓦片计算强度用于优化并行计算过程中负载均衡性能，细节将在2.3 节中介绍。

（6）矢量瓦片事实：和栅格瓦片事实共享3个维度成员包括产品、空间、时间，度量信息包括矢量瓦片计算强度以及每个矢量瓦片包含的矢量数据ID 集合。其中矢量数据ID 集合即瓦片内矢量数据的地址，指向矢量数据表中的数据。矢量瓦片计算强度用途将在2.3节中介绍。

立方体模型设计基于商业智能领域事实星座模型，该模型支持多个事实表共享不同的维度表，契合多源数据立方体模型设计理念（Chaudhuri 和Dayal，1997）。如图1所示，由于矢量数据没有波段维度，因此矢量和栅格数据共享产品、时间、空间3个维度，栅格数据独享波段维度。

2.1.2 立方体数据访问

基于立方体模型，用户可以在立方体上进行联机分析处理，比如切块和切片。这两种操作也是对地观测领域最常用的两种查询操作，其中切块操作是沿着不同的维度选择多个值得到子立方的过程，切片操作则是在维度上只选择一个值得到降维立方的过程（Baumann，2017；Voidrot 和Percivall，2020）。为方便用户访问数据，本文分别针对栅格瓦片和矢量数据设计了查询接口getRasterTiles（）和getFeatures（），接口以QueryParams类对象作为参数，该类通过多维语义查询参数初始化。图2展示了3个例子。

图2 立方体切块和切片操作Fig.2 Cube trim and slice operation

（1）切块：查询时间范围在2016年7月21日至2018年11月13日，空间范围在32.0°N—35.0°N，113.0°E—115.0°E 中的数据，返回落在该时空范围内的子数据立方；

（2）沿着一个维度切片：查询波段为蓝波段的数据，返回不包括波段维的降维数据立方；

（3）沿着两个维度切片：查询波段为蓝波段且时间为2017年6月11日的数据，返回不包括波段维和时间维的降维数据立方。

2.2 基于云计算的立方体数据处理

地理空间数据的海量、多源等特点，为立方体基础设施下的数据处理带来了挑战。事实上，本文立方体模型中的瓦片组织模式为大规模并行计算提供了便利，通过结合云计算技术，可以实现快速高效的数据处理。在立方体模型中，数据都沿着维度进行了切片，即立方体每一个单元都唯一对应着各个维度的一个成员。在空间维度上，栅格数据都物理上被切为瓦片，矢量数据采用逻辑切片的方式保证其拓扑完整性，即不真正对矢量切片，而是根据矢量数据与网格的空间关系对矢量数据进行组织，形成矢量逻辑瓦片。在时间维度上，每个瓦片对应着一个时刻。在产品维度上，每个瓦片属于某一种产品类别。在波段维度上，每个瓦片采用单波段的组织方式，即对应着某一个波段。

云计算作为兴起的大数据处理技术，已被广泛应用于各领域。其中Apache Spark作为云计算中的代表技术之一，是一种专为大规模数据处理而设计的高效计算引擎，采用弹性分布式数据集RDD（Resilient Distributed Dataset）作为数据容器来支持分布式计算（Zaharia 等，2010）。本文对Spark RDD 进行了扩展，提出了一套分布式立方体对象。一方面可以继承Spark RDD 分布式计算能力，另一方面可以兼容立方体多源数据类型。

如图3为设计的分布式立方体对象。位于最顶层的分布式内存立方体对象为CubeRDD［T］，其继承于原生Spark Abstract RDD，CubeRDD［T］是一个Key-Value 类型的RDD，其中Key 类型固定为SpaceTimeProductKey 类，该类包含地理空间数据最基础的3个维度信息即时间、空间和产品，其中时间信息采用Unix 时间戳表示，空间信息由Z 空间填充曲线编码生成，产品则直接采用字符串表示。比如用户查询Landsat 8 在武汉市2016年5月1日至6月1日的数据，则查询结果中的每一个SpaceTimeProductKey 由产品信息Landsat 8，时间信息为转换后的Unix 时间戳范围1462032000 至1464710400 中的某个值，空间信息为武汉市范围对应的Z 空间填充曲线编码中的某个值组成。CubeRDD［T］中的Value 类型为泛型T，满足不同来源数据的集成需求，可接受比如栅格瓦片对象RasterTile、矢量对象集合FeatureCollection 等。CubeRDD［T］还提供了其他属性包括网格空间布局、每个维度的范围以及参考坐标系。

图3 分布式内存立方体对象设计Fig.3 The design of distributed cube object

本文基于CubeRDD［T］，分别为栅格和矢量设计了RasterRDD 和FeatureRDD 分布式内存对象。RasterRDD 以CubeRDD［RasterTile］为构造参数，并可在时—空—产品键的基础上添加其他维度键比如波段信息，形成时间—空间—产品—波段键SpaceTimeProductBandKey。 RasterRDD 的值为RasterTile 类对象，RasterTile 类为ArrayTile 类的父类，ArrayTile 类存储了瓦片的像素数组，针对不同数据类型的像素数组，分别有ByteArrayTile、IntArrayTile、FloatArrayTile、DoubleArrayTile 等继承ArrayTile。针对矢量数据设计了分布式内存对象FeatureRDD，该类以CubeRDD［FeatureCollection］为构造参数，采用时—空—产品信息作为键，值为FeatureCollection 类对象，该类为Feature 类对象集合，每个Feature 类包含了几何信息Geometry 以及属性信息Attribute，FeatureRDD 可进一步转换为只包含几何信息的SpatialRDD。由于Geometry 可转换为点、线、面几何对象，因此SpatialRDD 又可转换为PointRDD、PolyLineRDD 和PolygonRDD。设计的分布式内存立方体对象具有良好的可扩展性，除了目前已支持的矢栅数据类型，可进一步扩展支持其他数据类型比如轨迹数据。

CubeRDD 的设计模式便于用户沿着不同维度进行分布式分析，比如可对RasterRDD 沿着时间维度即TimeKey做聚合，从而实现分布式地时间序列分析，同样可沿着波段维度即BandKey做聚合，从而实现分布式地波段融合分析。针对FeatureRDD，可沿着产品维度即ProductKey做聚合，从而可实现两个产品或图层的空间几何分析。另一方面，矢栅数据处于统一的时空基准中，因此可对RasterRDD和FeatureRDD 采用Spark 连接操作Join 基于时空键合并，使得具有相同时空键的栅格瓦片和矢量瓦片落在同一个时空立方体单元中，实现多源数据的联合分析。

2.3 人工智能优化的立方体计算

立方体模型采用全球空间网格基准对数据进行切片，以瓦片的形式来组织数据，便于实现分布式并行计算，充分调度起底层硬件资源。然而，地理空间数据尤其是矢量数据在空间上存在分布不规律、异质性突出等问题，易造成并行过程中负载失衡，难以充分利用硬件资源（Armstrong 和Densham，1992）。针对此问题，目前研究主要从两个角度去解决，包括域划分和动态调度。其中，动态调度适用于存在充足并行粒度的计算场景，这里的并行粒度即立方体中的瓦片。因此，动态调度方案可用于大尺度分析计算，比如以大量瓦片为输入的全球尺度分析，然而针对中尺度或小尺度分析比如洲际尺度、国家尺度，甚至省级或市级尺度，动态调度方案则难以实现较好的负载均衡性能（Bansal 等，2006；Ren 等，2017）。另一种策略，域划分方案实际属于调度前的一个阶段，常和静态调度策略结合，该方案通过评估并行粒度的计算强度，将这些并行粒度均匀地组合并分发到各个计算单元以实现负载均衡。域划分方案的关键问题是并行任务的计算强度评估，如评估不准确同样易造成负载失衡问题（Wang 和Armstrong，2009）。

针对域划分方案中的计算强度评估问题，作者已有相关研究发表于（Yue 等，2020）。研究从GIS 计算人工智能优化的角度出发，通过数据驱动的机器学习方法来刻画数据/算法特征，实现对地理空间域的合理表征，完成自动化预测并行地理计算任务计算强度。研究基于两个典型的地理空间分析案例，验证了提出方法的可行性与高效性，并提供了一个将机器学习方法用于地理空间分析任务计算强度评估的参考框架。具体来讲，首先需要根据经验知识，结合数据与算法特征，总结可能影响计算强度的候选特征。比如在矢量空间相交中，多边形的数量、多边形顶点数量、多边形分布方差等可作为候选特征；然后，框架通过结合多种机器学习特征选择算法和回归模型训练算法，根据精度评价指标从候选特征中选择出最优特征子集，同时得到最优预测模型；最后，在并行计算中便可通过模型预测计算强度，进行负载分配。需要注意的是，该方法中计算强度预测过程被嵌入到了并行计算过程中，存在一定的额外时间代价。

本文将该方法集成到立方体模型中，并通过提前预测划分粒度的计算强度，存储到数据库中，可消除原方法中的额外时间代价。方法集成过程如图4所示，首先，和原提出方法相同，通过结合机器学习特征选择技术以及训练算法，生成计算强度预测模型；然后，查询数据库中存储的瓦片数据，读取预先训练的模型，计算瓦片特征并预测计算强度，存储于图1中栅格或矢量瓦片事实表中（见2.1 节）。本文采用JSON 格式存储瓦片计算强度信息，并针对单图层分析和双图层分析分别设计了存储方法。其中单图层分析即对单个特征数据集进行的单层分析比如缓冲区分析，双图层分析则是对两个特征数据集进行的双层分析比如相交分析。单图层分析中，每个瓦片在不同的分析算法中具有不同的计算强度，因此分别存储不同分析算法对应的计算强度即可；双图层分析中，瓦片计算强度不仅依赖于分析算法，而且依赖于和其进行分析的另一个瓦片，因此需要根据分析算法和瓦片来区分计算强度，如图4所示。最后，智能并行计算阶段，计算强度连同瓦片数据一起获取，通过自定义分区器，即根据分区数目，将瓦片按照计算强度值均匀地分配给各个分区。分配过程中首先对所有瓦片的计算强度求和，然后按照分区数量得到分区平均计算强度，最后将瓦片依次放入到分区，直到该分区的累计计算强度大于平均计算强度。

图4 人工智能优化的立方体计算模式Fig.4 Artificial-intelligence optimized cube computing

3 实现与结果分析

3.1 GeoCube框架

本文提出了一种面向大规模分析的多源数据对地观测时空立方体基础设施GeoCube，为对地观测领域中的多源数据管理、多维语义查询、大规模分析等问题提供了一种解决方案。如图5所示为GeoCube 整体框架，分为基础设施层、数据接入层、数据存储层、数据视图层、计算引擎层和访问接口层。

图5 GeoCube框架Fig.5 The framework of GeoCube

（1）基础设施层：该层为底层硬件基础，GeoCube建立在一个由高性能计算集群和大规模存储阵列组成的高性能基础设施之上。

（2）数据接入层：GeoCube目前支持Landsat 8、高分一号、哨兵一号SAR、OSM 和专业领域矢量数据的接入，其中影像数据预处理为分析就绪型数据级别。

（3）数据存储层：GeoCube 采用统一的时空基准，基于全球尺度定义的空间网格对矢栅两种数据进行切分。根据自定义的元数据规范，基于多维数据组织模型，分别采用关系数据库和分布式数据库存储维度属性数据与原数据。

（4）数据视图层：GeoCube 支持从空间维度、时间维度、产品维度、波段维度向用户提供多维语义查询接口。

（5）计算引擎层：该层结合分布式计算和人工智能领域技术，优化大规模计算和异质性较强的并行计算，提供地图代数运算、时间序列分析、空间几何分析和联合分析处理等功能。

（6）访问接口层：GeoCube 通过提供Jupyter Notebook在线脚本控制台、标准网络空间服务以及Web平台的方式，方便用户访问GeoCube API。

3.2 数据与环境

GeoCube中的遥感影像主要覆盖2013年—2019年中国部分区域，包括Landsat 8 OLI 30 m 分辨率数据、高分一号PMS 2/8 m 分辨率和WFV 16 m 分辨率数据，以及哨兵1 号SAR 数据。矢量数据主要来自于OSM提供的数据。

硬件环境方面，GeoCube 建立在3 台物理服务器上，包括2 台Inspure NF5280M4 服务器，配有Intel（R）Xeon（R）CPU E5-2690 v4@2.60 GHz，Centos 7.3.1611操作系统，28个CPU 核以及128 GB内存，另外一台为Inspure NF5270M3服务器，配有Intel（R）Xeon（R）CPU E5-2692 v2@2.20 GHz，Ubuntu 14.04.1 操作系统，24 个CPU 核以及94 GB内存。3 台服务器和一个PB 级存储阵列连接，该存储阵列由84块16 TB的存储磁盘组成。

软件环境方面，GeoCube 采用PostgreSQL 9.3.24 存储维度信息以及事实信息，Hadoop 2.7.4和HBase 1.4.3 存储瓦片数据，采用Spark 2.4.3 和Geotrellis 3.0.0 实现分布式计算，基于Springboot 2.3.1 将GeoCube 数据和算法发布为OGC 覆盖服务和处理服务，采用Almond 实现Scala 语言Jupyter Notebook，采用Cesium1.71实现数据可视化。

3.3 效率评估

3.3.1 查询效率

GeoCube 支持对数据进行多维度语义查询，如图6所示，展示了对数据的查询界面与可视化效果。其中图6（a）所示为GeoCube 目前已导入的产品信息，包括各类栅格产品与矢量产品；图6（b）所示为对某一景产品叠加到Cesium 三维地球上的可视化效果；图6（c）为在Jupyter Notebook 中基于GeoCube内核对瓦片数据进行多维语义查询，返回结果可以是一个瓦片数组也可以是一个CubeRDD［T］。

图6 GeoCube中数据查询与可视化Fig.6 Data retrieval and visualization in GeoCube

在大规模分析场景中，用户往往需要使用大量的瓦片数据，GeoCube采用了混合并行的方式来提高数据访问速度。如图7所示为数据查询流程，包含元数据查询以及瓦片数据访问两个阶段，其中元数据查询阶段即获取瓦片地址以及属性信息，瓦片数据访问阶段即将HBase中的瓦片数据反序列化为瓦片内存对象。在元数据查询阶段，首先根据用户提供的维度语义参数，访问各个维度表得到维度键，然后到事实表中查询瓦片地址即瓦片ID，得到瓦片ID 集合。在瓦片数据访问阶段，瓦片ID 集合被分发到各个节点，各节点根据分配的瓦片ID 分布式地反序列化HBase 数据表中的瓦片。在该阶段，瓦片分辨率过高会导致单张瓦片反序列化时间过长，因此GeoCube在每个任务中进一步采用多线程技术，每个线程负责反序列化瓦片的部分数据。

图7 GeoCube瓦片查询流程Fig.7 Data access flow in GeoCube

图8 对比了地理立方GeoCube 与开源数据立方OpenDataCube 的查询效率，统计了将数据从磁盘读取到内存中的时间。OpenDataCube 在数据量较小时查询效率略优于GeoCube，但随着数据量增大，GeoCube 在查询效率上和OpenDataCube 基本持平。造成该现象原因为GeoCube采用的分布式存储模式，相比于OpenDataCube 采用的数据集中存储在文件系统中的方式，GeoCube在分布式查询过程中存在更多的网络通信，比如某一并行任务所查询的数据不在当前节点，则需要从其他节点获取数据。分布式存储模式虽然给查询效率带来了一定影响，但却可以大大提升计算效率，因为数据已经提前分布在了分布式节点中，减少了节点间的数据通信，从而整体上具有优势，计算效率评估见3.3.2节。

图8 地理立方和开源数据立方查询效率对比Fig.8 Comparison of data query performance between GeoCube and OpenDataCube

3.3.2 处理效率

GeoCube的一个重要目标就是向用户提供快速的大规模分析。本文以水体观测频率分析为例，验证GeoCube的高效计算能力。水体观测频率分析属于一种时间序列分析，它通过计算每个像素在时间序列上被检测为水体的频率，来反映不同区域在一段时间内水体覆盖的整体情况（Mueller等，2016）。如图9所示为该分析的并行执行流程，由于待处理的瓦片数据量可能较大导致无法一次性读入到内存，本文采用了按空间进行分批的方式，解决内存不足问题的同时保证每个空间网格上具有查询时间范围内全时间序列数据。每个批次的瓦片数据已分布在各个节点上，因此可直接基于空间关系聚集时间序列瓦片数据，分布式地进行水体观测频率分析。但由于进行了空间分批，每个批次内的瓦片数量有限，会导致无法充分调度起硬件资源。针对该问题，GeoCube在分布式处理过程中，在每个分布式任务里结合了多线程技术，每个线程处理时间序列上部分空间范围内瓦片，从而进一步提高处理速度。

图9 GeoCube水体观测频率分析流程Fig.9 Water observation frequency analysis in GeoCube

图10 对比了地理立方GeoCube 与开源数据立方OpenDataCube 在水体观测频率分析上的性能。OpenDataCube 采用和GeoCube 相同的硬件环境，通过Dask 技术实现高性能并行处理（DASK，2020）。如图10所示，处理数据量从GB 到TB 级别，GeoCube全程处于优势，并且随着数据量的增多优势逐渐增大。GeoCube具有处理优势的原因一方面得益于上文所提的混合并行的处理模式；另一方面是GeoCube采用的分布式存储模式，相比于OpenDataCube 将数据集中存储在文件系统的方式，GeoCube在并行处理前，数据已分布在各节点，导致处理过程具有较小的通信代价。图11 展示了湖北省部分区域在2013年—2018年的水体观测频率分析结果，深色部分表示该区域大部分时间都表现为水体，而浅色部分表示该区域常年干旱。

图10 地理立方和开源数据立方处理效率对比Fig.10 Comparison of processing performance between GeoCube and OpenDataCube

图11 2013年—2018年湖北省部分区域的水体观测频率Fig.11 Water observation frequency result in some regions of Hubei province from 2013 to 2018

3.3.3 负载均衡性能

相较于栅格数据，矢量数据具有空间分布不均、异质性等特点，极易造成矢量处理算法在并行化过程出现负载失衡问题。本节将以矢量数据空间相交算法为例，给出2.3 节基于人工智能优化的方法在GeoCube 中负载均衡问题上的性能表现。

根据立方体存储模型，立方体中矢量数据已基于全球统一空间网格进行了划分，可直接用于分布式并行空间相交，这也是目前SpatialSpark 等主流框架采用的方法（You 等，2015），即空间网格划分后直接进行并行计算。但这种方法忽略了矢量切片的计算强度差异，完全依赖于Spark 自身调度机制来实现负载均衡，当切片计算强度相差较大或切片数量不够充足时，容易造成处理单元负载失衡。本节将对比这种未经优化的方法和经过2.3 节基于人工智能优化后的方法在不同数据集上的负载均衡性能。

实验数据来自于OSM 提供的2020年中国和印度区域矢量数据，均已按照空间网格切分并导入GeoCube。实验分为3组测试，第1组测试中两个待相交图层分别包含169万和49万个多边形，第2组测试中两图层分别包含250 万和82 万个多边形，第3组测试中两图层分别包含300万和87万个多边形。如图12所示，两种方法在数据查询阶段的耗时基本相同，然而在数据处理阶段，随着数据量的增大，基于人工智能优化后的方法相较于未经优化的方法优势逐渐增大。图13 给出了两种方法在3 组测试中的负载均衡表现，实验统计了8 个执行器上的执行时间，可见基于人工智能优化后的方法负载均衡性能明显优于未经优化的方法。

图12 未经优化和人工智能优化方法计算效率对比Fig.12 Comparison of computing performance between the AI-optimized approach and the traditional approach

图13 未经优化和人工智能优化方法负载均衡对比Fig.13 Comparison of load balance performance between the AI-optimized approach and the traditional approach

3.4 应用案例

3.4.1 长时序光谱指数产品

光谱指数计算作为遥感领域广泛应用的分析方法之一，通过对两个相关波段进行操作来监测陆表地物覆盖情况。常用的光谱指数包括归一化差异水体指数NDWI（Normalized Difference Water Index）（Gao，1996；McFeeters，1996）、改进归一化差异水体指数MNDWI （Modified Normalized Difference Water Index）（Xu，2006）、归一化差异植被指数NDVI（Normalized Difference Vegetation Index）（Crippen，1990）、增强植被指数EVI（Enhanced Vegetation Index）（Huete等，2002）、土壤调整植被指数SAVI（Soil Adjusted Vegetation Index）（Huete，1988）等。

NDWI是目前应用最广泛的一种水体指数，包括两种计算方式，分别为通过组合近红外和绿波段消除土壤和地表植被要素背景来提取水体（McFeeters，1996），以及通过使用近红外和短波红外来提取植物液态水（Gao，1996）。其中用于水体提取的NDWI 在建筑物背景较多时效果较差，会扩大实际水体面积。MNDWI 在此基础上做了改进，可区分阴影和水体，并且可提取水体悬浮沉积物分布等细微特征（Xu，2006）。NDVI 是植被监测中最常用的光谱指数，它可以部分消除太阳高度角、传感器观测角和大气等影响，但对土壤背景变化敏感。目前对植被指数的改进方法较多，比如顾虑到土壤亮度影响的SAVI，针对冠层结构变化设计的EVI等。

GeoCube 目前已实现了NDWI，MNDWI，NDVI以及SAVI，用户可基于GeoCube JupyterNotebook自定义实现更多光谱指数。如图14 为使用GeoCube 生成的湖北省某区域2015年—2018年的NDWI部分专题产品。

图14 2015年—2018年湖北省某区域NDWI部分专题产品Fig.14 NDWI product in some regions of Hubei province from 2015 to 2018

3.4.2 多源数据联合分析

本节将以2016年中国南部某省洪涝灾害为案例，验证GeoCube 的多源数据联合分析能力。该案例背景为2016年台风电母登录中国南部，造成大规模降雨，导致诸多城市出现洪涝灾害，约50 万人受到此次灾害影响，大量的交通设施遭到破坏。

本节以该案例为基础，通过矢栅联合分析来提取受灾基础设施。数据包括灾前灾后高分一号数据以及部分交通设施矢量数据，比如桥梁、道路等。受灾要素提取流程为：基于灾前灾后的高分遥感影像，生成水体提取产品；对灾前灾后水体提取产品做水体变化检测，得到洪涝受灾区域；将洪涝受灾区域产品与矢量数据做叠置分析，提取受灾要素。如图15所示，红色区域为受到洪涝淹没的区域，蓝色箭头为该地区受到洪涝影响的桥梁。图16 给出了在Jupyter Notebook 中该案例的分析流程，涉及GeoCube 内核导入，多维数据查询，受灾要素提取以及专题产品结果展示等流程。

图15 中国南部某地区受灾桥梁提取结果Fig.15 The affected bridge extraction in a region of Southern China

图16 Jupyter Notebook中案例实现流程Fig.16 The case implemented in Jupyter Notebook

4 结论

在对地观测领域，数据立方基础设施正作为一项突破性的技术推动着新一代的模式变革。本文自主研发了一套时空立方基础设施GeoCube，提供了针对海量多源数据存储、访问、分析和高性能计算的方案，其涉及一种面向多源数据的多维度时空组织策略，一套分布式内存立方体对象，一种基于人工智能优化的并行计算方法。本文设计了多组实验用以验证GeoCube 的可行性与高效性，得出以下结论：

（1）基于事实星座模型的多源数据多维度时空组织策略，在功能上一方面支持多维度语义查询，包括产品、时间、空间、波段维度；另一方面支持矢量栅格等多源数据管理；

（2）基于分布式计算引擎的分布式内存立方体对象设计，在功能上具有良好的可扩展性，不仅面向遥感影像数据，同时兼容矢量数据映射；在性能上通过结合更细粒度的线程级并行，可实现高效的处理分析；

（3）基于人工智能优化的并行计算方法，在功能上面向异质性较强的算法/数据，比如矢量空间相交；在性能上通过预先预测的计算强度信息，可实现较好的负载均衡。

目前，GeoCube 存储的数据量和管理的数据类型仍然有限，未来将接入更多的数据，并增加支持数据类型，比如轨迹数据。另外，数据立方间的互操作问题是目前对地观测领域面临的主要挑战之一，未来将结合GeoCube 和其他框架的立方体，探索互操作问题的解决方案（Nativi 等，2017；Voidrot和Percivall，2020）。