CN112100130B - 一种基于数据立方体模型的海量遥感变量多维聚合信息的计算方法 - Google Patents

Abstract

本发明公开了一种基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，属于数据处理领域。本发明的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，采用预先聚合遥感变量信息，将大部分聚合计算转换成为已有聚合信息的查询过程，通过预先计算的多个小粒度聚合信息自动组合成任意查询条件对应的聚合结果，从而避免了传统方法所面临的计算性能问题。因此，本发明能够在人机交互场景中以近实时的响应速度返回遥感变量聚合信息，同时节省大量计算资源。本发明支持时间聚合和空间聚合，也同时支持时空独立聚合与时空耦合聚合等复杂形式。

Description

一种基于数据立方体模型的海量遥感变量多维聚合信息的计 算方法

技术领域

本发明属于数据处理领域，尤其是一种基于数据立方体模型的海量遥感变量多维聚合信息的计算方法。

背景技术

当前，高时/空/谱卫星遥感技术与地理模拟系统的迅速发展，产生了海量高维的对地观测数据及其再分析数据集，并被广泛应用于气候灾害、环境生态等各领域。在全球气候变化与人地关系日益耦合的背景下，自然观测数据与人文统计数据的联合分析已成为多学科交叉研究的重要方向。时空聚合(Spatiotemporal Aggregation)是实现自然-人文多元地理数据集成分析的重要手段，它通过对数据集进行信息抽取与汇总，有效去除冗余信息，进而将多源异构数据纳入到统一的时空分析尺度空间中。根据聚合维度的组合方式，时空聚合分为时间聚合、空间聚合、时空独立聚合、时空耦合聚合，其分析范围往往涉及复杂的空间几何(如土地利用边界)、不均等时间区域(如冷暖季时长不同)、相互耦合的时空边界(如不同区域作物类型具有不同物候)以及多重属性值域条件(如“植被指数>0and土壤类型＝[1,2]”)。

当前大部分GIS软件或遥感图像软件不具备时空一体化的地图代数运算能力，仅能解决时间聚合与空间聚合问题。对于时空独立聚合、时空耦合聚合以及时空多维聚合，现有研究一般采用多维数组/多维地图代数运算实现，但是需先将复杂时空区域转换为多维矩阵，存在着高内存消耗与高计算延迟问题。在面对海量、高维的遥感信息时，时空聚合的查询性能低下，尤其是在交互式环境下的实时响应就更加困难。

发明内容

本发明的目的在于解决时空独立聚合、时空耦合聚合以及时空多维聚合，需先将复杂时空区域转换为多维矩阵，存在着高内存消耗与高计算延迟的问题，提供一种基于数据立方体模型的海量遥感变量多维聚合信息的计算方法。

为达到上述目的，本发明采用以下技术方案予以实现：

一种基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，包括以下步骤：

1)采用MultiDim E-R方法设计时间维度层级、空间维度层级以及遥感变量维度层级，以FreeMan-Diaconis频率直方图作为数据立方体模型的度量，建立面向遥感变量聚合信息的数据立方体模型，并生成数据立方体模型的配置文件；

以分块压缩的NetCDF文件格式存储不同粒度的数据立方体模型子集和原始遥感变量；

所述配置文件包括数据立方体模型描述、数据立方体模型与原始数据之间的映射关系；

2)根据输入的查询时空区域建立数据立方体模型维度层级拓扑包含的最小深度聚合树，所述最小深度聚合树中的每个叶子节点为完全处于查询时空区域内的自适应时空单元；

根据所述自适应时空单元对应的编码，获取自适应时空单元中的直方图信息，利用直方图信息推算出时空单元内遥感变量的聚合信息，之后利用对应的分布聚合公式计算得到查询时空区域内部的遥感变量的近似聚合结果；

对于与查询时空区域边界相交的时空单元，利用多维地图代数进行数据单元内的时空聚合计算，得到时空聚合结果；

3)将时空聚合结果与所述近似聚合查询结果进行合并，得到精确的聚合结果。

进一步的，步骤1)中，NetCDF格式文件存储在本地服务器，以单机形式来支持访问；

或者以对象存储形式部署于分布式云环境中。

进一步的，步骤1)中的所述配置文件记录有维度配置信息、度量配置信息、不同维度粒度组合对应的NetCDF文件存储路径、原始遥感变量的NetCDF文件存储路径；

所述维度配置信息包括维度类型、维度范围和层级结构；

所述度量配置信息为直方图统计区间的粒度。

进一步的，步骤1)中的所述配置文件以XML格式部署于本地服务器或分布式云环境中。

进一步的，步骤2)中输入查询时空区域的方式为，在构建交互式地图界面上：

通过绘制面状图形或上传shapefile/GeoJSON格式多边形文件指定空间区域；

通过时间滑动轴指定时间区域；

通过上传带有起止时间字段的shapefile/GeoJSON格式多边形文件来指定与空间耦合的时间区域；

通过多选框指定遥感变量，通过对应的值域滑动轴指定不同遥感变量的值域。

进一步的，步骤2)建立数据立方体模型维度层级拓扑包含的最小深度聚合树，具体操作为：

根据输入的查询时空区域，从数据立方体模型时空维度的最高层级逐层进行时空拓扑计算，获取与查询时空相交的时空单元，对于与区域边界相交的时空单元向下层细化，直至时空单元完全处于查询时空区域内，得到最小深度聚合树，并获取查询范围内自适应时空单元集合。

进一步的，步骤2)中所述的聚合信息，具体为求和、计数、均值、众数、中位数和方差。

进一步的，其特征在于，步骤2)中利用直方图信息推算出聚合信息的具体操作为：

利用查询条件中遥感变量的阈值，对各时空单元的直方图信息进行查询，过滤掉阈值范围外的时空单元，并在最小深度聚合树中更新满足查询条件的时空单元直方图信息。

与现有技术相比，本发明具有以下有益效果：

本发明的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，采用预先聚合遥感变量信息，将大部分聚合计算转换成为已有聚合信息的查询过程，通过预先计算的多个小粒度聚合信息自动组合成任意查询条件对应的聚合结果，从而避免了传统方法所面临的计算性能问题。因此，本发明能够在人机交互场景中以近实时的响应速度返回遥感变量聚合信息，同时节省大量计算资源。本发明支持时间聚合和空间聚合，也同时支持时空独立聚合与时空耦合聚合等复杂形式；以交互式在线地图界面为用户提供遥感变量的各种聚合功能，屏蔽了算法处理的专业性，增强了时空聚合功能的实用性，提升了用户使用体验。

附图说明

图1为基于数据立方体模型的遥感变量聚合信息多尺度组织方法的流程图；

图2为多尺度数据空间时空划分示意图；

图3为面向遥感变量聚合信息的数据立方体模型示意图；

图4为海量遥感变量时空聚合交互式查询方法流程图；

图5为遥感变量精确聚合查询流程图；

图6为第一个测试项目的结果图。

具体实施方式

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

下面结合附图对本发明做进一步详细描述：

本发明的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，可分解成基于数据立方体模型的遥感变量聚合信息多尺度组织方法和海量遥感变量时空聚合交互式查询方法进行。

参见图1，图1为基于数据立方体模型的遥感变量聚合信息多尺度组织方法的流程图，包括以下实施步骤：

S1.利用时间维度层级、空间格网层级和遥感变量维度对数据空间进行多尺度划分，具体地：

本发明实施例中时间维度层级定义为“日<月<年<ENSO周期”的树形结构(采用10位字符编码，如A120000918代表ENSO暖季2000年9月18日)；空间格网层级定义为“0.2°<0.4°<0.8°<1.6°”多级全球四叉树剖分格网(采用Z曲线四叉树QuadTile编码，如04accd代表0.2°格网级别与编码为accd的空间单元)；遥感变量维度为气温、降雨、湿度、风速等天气变量类型(采用英文名称编码，如temperature代表温度)。

根据时空维度层级上时间粒度和空间粒度的排列组合，将数据空间划分为不同尺度上的时空单元，如图2所示，时空单元按照粒度粗化方向存在着<0.2°,day>，<0.2°,month>，<0.2°,day>，<0.4°,day>，<0.4°,month>，<0.4°,year>，<0.8°,day>，<0.8°,month>，<0.8°,year>等多种粒度组合，上层粒度可由下层粒度按照固定的维度层级方向聚合而得到。

遥感变量维度按照变量类型默认保持一级分类，也可以对每个变量类型扩展二级或多级分类，如基于风速值定义风速的语义分级区域“[无风,微风,中风,大风]”。

S2.以多尺度数据空间划分单元为数据组织单位，以单元内遥感变量直方图为度量类型建立数据立方体模型。

数据立方体模型通过在数据库或内存中构建多维交互模型和多粒度聚合机制，为海量多维信息的洞察分析提供了强有力的数据引擎。本发明提出建立面向遥感变量聚合信息的数据立方体模型。

采用MultiDim E-R方法设计时间维度层级、空间维度层级以及遥感变量维度，同时以FreeMan-Diaconis频率直方图作为数据立方体模型的度量，如图3所示，空间维层级(Spatial Grid)之间是空间包含拓扑关系，Geocode记录了空间格网编码；时间维层级(Observed Time)之间是时间包含拓扑关系，ID记录时间编码；遥感变量维度(Multivariate)仅有单层级，Name记录变量类型编码。事实表(Cube Fact)中存储了各维度组合对应数据单元所统计的FreeMan-Diaconis频率直方图。

根据维度粒度组合的不同，对于每种数据空间划分，以时间编码T、空间编码S、变量编码V为三个维度建立支持标签的多维数组[T][S][V]，以存储数据立方体模型的一个子集。其中，每个数组单元内存储的值为H，其中，H为数据空间单元内变量值的FreeMan-Diaconis频率直方图，采用一维列表[<x_i,f_i>]结构组织，f_i为遥感变量在直方图区间x_i上的频数。细粒度划分单元内的直方图可以通过基于求和计算(SUM)的聚合过程直接得到粗粒度单元内的直方图。

S3.以压缩后的多维数组数据格式存储不同时空尺度的数据立方体模型，以数据立方体模型配置文件建立数据立方体模型描述及其与原始数据之间的关系映射。

以分块压缩的NetCDF文件格式存储不同粒度的数据立方体模型子集。其中，NetCDF格式文件可存储在数据聚合引擎所在的本地服务器，以单机形式来支持访问，也可以对象存储形式部署于分布式云环境中。原始遥感变量同样以分块压缩的NetCDF文件格式进行存储，但是一个NetCDF文件只存储单一变量。

数据立方体模型配置文件以XML格式部署于本地服务器或分布式云环境中。该配置文件记录了维度配置信息(维度类型、维度范围、层级结构)、度量配置信息(直方图统计区间粒度)。记录了不同维度粒度组合对应的NetCDF文件存储路径。记录了原始遥感变量的NetCDF文件存储路径。

数据聚合引擎在内存中加载数据立方体模型配置文件，建立数据立方体模型虚拟对象，然后以Lazy load方式读取对应的NetCDF文件分块部分，生成多维数组，从而执行后续的聚合计算。

参见图4，图4为海量遥感变量时空聚合交互式查询方法的流程图，包括以下实施步骤：

S1.计算覆盖用户通过交互式地图界面指定时空区域范围内数据立方体模型维度层级所对应的最小深度聚合树，从而获取多粒度时空划分单元集合，以自适应描述用户指定的任意时空区域。

构建交互式地图界面，用户通过手动绘制面状区域或上传shapefile/GeoJSON格式多边形文件指定空间区域，通过时间滑动轴选取统一时间区域，或通过上传带有起止时间字段的shapefile/GeoJSON格式多边形文件来指定与空间耦合的时间区域。通过多选框指定遥感变量，通过对应的值域滑动轴指定不同遥感变量的值域。

对于用户自定义的时空范围，查询算法将自动从数据立方体模型时空维度的最高层级逐层进行时空拓扑计算，获取与指定区域相交的时空单元，对于与区域边界相交的单元向下层细化，直至单元处于完全被包含关系，最终获取指定区域覆盖的自适应时空单元集合{<T,S>}，以近似描述其时空边界范围，作为时空条件记录在查询请求JSON中，其中，T为时间单元编码，S为空间单元编码。对于自定义的多个时空聚合区域，将以<ZS,ZT,V>编码形式记录统计区域，其中ZS为空间聚合区域，ZT为时空聚合区域，V为变量类型。

对于用户自定义的遥感变量阈值范围{<V,[p₁,p₂]>}，p₁和p₂是起止值，作为变量条件也记录在查询请求JSON中。

S2.根据查询区域的多维单元编码集合获取各单元存储的直方图信息，利用直方图信息推算单元内遥感变量的求和、计数、均值等各类聚合信息，然后通过对应的分布聚合公式，实时计算得到查询时空区域内部的遥感变量近似聚合结果。

对于上述查询请求，数据聚合引擎将首先获取数据立方体模型对应编码位置的时空单元，生成多维数组。基于单元内直方图信息，利用遥感变量阈值范围对各单元进行过滤，同时对过滤后的直方图信息进行更新，随后执行下述统计。

假定某数据单元内的统计直方图为[<x_i,f_i>]，由频数分布直方图信息可以推导求解该数据样本对应的各类聚合算子值，如表1所示。其中众数、中位数和方差是由直方图分组数据计算出来的近似值，且其精度随着组距变小而提高。

表1时空单元直方图的各类聚合信息计算方法

根据S1结果，查询区域覆盖多个时空单元，根据每个单元内的直方图分布聚合信息统计整体区域的聚合值。上述聚合算子根据可分布聚合性判定可分为三类：①可分布聚合(Distributive)，可通过子集分布式函数直接计算全集聚合值，例如Sum和Count；②可代数聚合(Algebraic)，即通过子集代数函数组合计算全集聚合值，具有间接可分布聚合性，例如Avg＝Sum/Count；③整体聚合(Holistic)，即只能基于全集计算聚合值，例如中位数Med。因此，整体区域的总数、总和与均值可以直接由分布聚合公式统计而来。整体区域的众数、中位数和方差需要基于多个时空单元内直方图的最终聚合结果而计算得到。

S3.对于与时空区域边界相交的数据单元，利用多维地图代数进行单元内的时空聚合计算，并将计算结果与上述S2中的近似聚合查询结果进行合并，最终得到精确的聚合结果。

由于时空查询区域覆盖了区域内部的时空单元，同时也包含与区域边界相交的时空单元。对于区域内部的时空单元，通过上述S2步骤可统计出覆盖时空查询区域的近似聚合信息<Q₁,Q₂,…,Q_n>，其中Q_n为第n个时空单元的聚合值，如图5所示。

对于区域边界上的时空单元，需要通过数据立方体模型配置文件获取对应时空位置的原始数据。以该时空单元内的原始数据为统计目标，以区域边界为统计区域，基于多维地图代数Zonal函数计算该单元内的分布直方图信息和聚合值<C₁,C₂,…,C_k>，其中C_k为第k个时空单元的涉及区域聚合值。若统计区域具有m个部分，则该时空单元分解为m个子单元，对应统计得到m个分布直方图。

最后，根据统计区域编码，对上述区域内部时空单元统计结果，和区域边界时空单元统计结果进行分组汇总。汇总聚合方法同样参考B2中的分布聚合公式，从而得到用户查询区域内遥感变量的精确聚合值。

实施例

为了验证本发明的数据组织与聚合方法对于遥感变量聚合查询应用是否满足交互式环境下近实时响应的需求，采用光纤局域网环境、服务主机(配置：CPU i7-7700 3.6G，内存32G，存储1TB SSD)和若干客户端进行测试。测试原始数据包含3个变量(最大气温，最小气温，降雨量)，空间范围为全球区域，空间分辨率为0.04°，时间跨度为58年，时间分辨率为月，总数据量在350GB左右。根据本发明数据模型建立数据立方体，空间层级为”0.7°<1.4°<2.8°<5.6<11.2°”，时间层级为”月<年<ENSO周期”。测试方法采用3种，本发明方法简称为Cube，两个传统方法分别为：①在同一测试机上编写ArcGIS python脚本依次执行空间聚合、时间聚合与变量阈值过滤，以模拟时空聚合(简称ArcGIS)。②基于Python Numpy包和Geopandas包实现多维地图代数Zonal函数，执行单次时空聚合计算(简称MMA)。

利用以上三种方法分别执行以下测试项目3次：①通过在GIS交互界面依次绘制多个任意多边形(面积约为10，20，40,80，160，300万平方公里)，指定10年的时间区域，查询每个时空区域内3种遥感变量的年均值与中位数值；②一次性查询全球范围所有国家30年ESNO平均降水量。

测试过程中，记录并按均值统计服务器内存变化量与查询延迟。如图6所示，图6为第一个测试项目的结果，数据立方体响应延迟、内存利用率以及CPU利用率最小，且随查询面积增加而增长缓慢；ArcGIS响应延迟最长，内存利用率呈线性增长，CPU利用率增幅不大；多维地图代数在初始状态下响应延迟较小，但随着查询面积增大，响应延迟、内存利用率和CPU利用率急剧增加。如表2所示，表2为第二个测试项目的数据，第二个测试项目结果显示在大量数据聚合情况下，多维地图代数出现严重性能问题，导致内存溢出，程序无法顺利执行；ArcGIS采取多批次执行模式，因而能完成计算任务，但是耗时较长；本发明提出的方法响应延迟为45s，能够较好完成查询任务，具有明显优势。

表2第二个测试项目的测试结果

测试方法	响应延迟(s)	内存利用率％	CPU利用率％
				ArcGIS	9188	45.7	30.5
多维地图代数	/	/	/
				数据立方体(本发明)	45	38.6	35.2

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (8)

1.一种基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，包括以下步骤：

1)采用MultiDim E-R方法设计时间维度层级、空间维度层级以及遥感变量维度层级，以FreeMan-Diaconis频率直方图作为数据立方体模型的度量，建立面向遥感变量聚合信息的数据立方体模型，并生成数据立方体模型的配置文件；

以分块压缩的NetCDF文件格式存储不同粒度的数据立方体模型子集和原始遥感变量；

所述配置文件包括数据立方体模型描述、数据立方体模型与原始数据之间的映射关系；

2)根据输入的查询时空区域建立数据立方体模型维度层级拓扑包含的最小深度聚合树，所述最小深度聚合树中的每个叶子节点为完全处于查询时空区域内的自适应时空单元；

根据所述自适应时空单元对应的编码，获取自适应时空单元中的直方图信息，利用直方图信息推算出时空单元内遥感变量的聚合信息，之后利用对应的分布聚合公式计算得到查询时空区域内部的遥感变量的近似聚合结果；

对于与查询时空区域边界相交的时空单元，利用多维地图代数进行数据单元内的时空聚合计算，得到时空聚合结果；

所述最小深度聚合树由以下方式得到：对于用户自定义的时空范围，从数据立方体模型时空维度的最高层级逐层进行时空拓扑计算，获取与指定区域相交的时空单元，对于与区域边界相交的单元向下层细化，直至单元处于完全被包含关系，最终获取指定区域覆盖的自适应时空单元集合{<T,S>}，作为时空条件记录在查询请求JSON中，其中，T为时间单元编码，S为空间单元编码；对于自定义的多个时空聚合区域，将以<ZS,ZT,V>编码形式记录统计区域，其中ZS为空间聚合区域，ZT为时空聚合区域，V为变量类型；

3)将时空聚合结果与所述近似聚合查询结果进行合并，得到精确的聚合结果。

2.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤1)中，NetCDF格式文件存储在本地服务器，以单机形式来支持访问；

或者以对象存储形式部署于分布式云环境中。

3.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤1)中的所述配置文件记录有维度配置信息、度量配置信息、不同维度粒度组合对应的NetCDF文件存储路径、原始遥感变量的NetCDF文件存储路径；

所述维度配置信息包括维度类型、维度范围和层级结构；

所述度量配置信息为直方图统计区间的粒度。

4.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤1)中的所述配置文件以XML格式部署于本地服务器或分布式云环境中。

5.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤2)中输入查询时空区域的方式为，在构建交互式地图界面上：

通过绘制面状图形或上传shapefile/GeoJSON格式多边形文件指定空间区域；

通过时间滑动轴指定时间区域；

通过上传带有起止时间字段的shapefile/GeoJSON格式多边形文件来指定与空间耦合的时间区域；

通过多选框指定遥感变量，通过对应的值域滑动轴指定不同遥感变量的值域。

6.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤2)建立数据立方体模型维度层级拓扑包含的最小深度聚合树，具体操作为：

根据输入的查询时空区域，从数据立方体模型时空维度的最高层级逐层进行时空拓扑计算，获取与查询时空相交的时空单元，对于与区域边界相交的时空单元向下层细化，直至时空单元完全处于查询时空区域内，得到最小深度聚合树，并获取查询范围内自适应时空单元集合。

7.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤2)中所述的聚合信息，具体为求和、计数、均值、众数、中位数和方差。

8.根据权利要求1所述的基于数据立方体模型的海量遥感变量多维聚合信息的计算方法，其特征在于，步骤2)中利用直方图信息推算出聚合信息的具体操作为：

利用查询条件中遥感变量的阈值，对各时空单元的直方图信息进行查询，过滤掉阈值范围外的时空单元，并在最小深度聚合树中更新满足查询条件的时空单元直方图信息。

Images