顾及计算强度的地质矢量大数据计算域 均衡分解方法研究

朱月琴，郭明强，黄 颖，吴 亮,5，谢 忠,5

(1.自然资源部地质信息技术重点实验室，北京 100037； 2.中国地质调查局发展研究中心，北京 100037； 3.中国地质大学(武汉)信息工程学院，湖北 武汉 430074； 4.武汉中地数码科技有限公司，湖北 武汉 430074； 5.国家地理信息系统工程技术研究中心，湖北 武汉 430074)

0 引 言

随着科学研究方式第四范式的诞生，即数据密集型的知识发现，地质学领域具有多源、异构等复合型数据已经被列入大数据的范畴，被称之为“地质大数据”[1-2]。大数据量地理信息数据在线服务应用提供的矢量服务，其响应速度和数据量、服务器性能，网络环境相关性大，目前矢量大数据服务应用存在的瓶颈有以下四点：数据发布准备工作费时费力[3-4]、时效性难以保证[5-6]、无法在线编辑处理、难以支持大数据量应用[7-9]，迫切需要研究适应矢量大数据实时渲染的并行处理方法。

1 地质矢量大数据实时渲染系统总体框架

面向地质大数据实时渲染的计算域均衡分解系统层次架构的设计包括：数据准备层、前端应用层、业务逻辑层、GIS服务层。数据准备层为独立的桌面软件，与后面层级没有直接的关系，但后面层级在运行之前必须首先通过数据准备层生成必要的准备数据。前端应用层通过ajax向业务逻辑层发送http请求，业务逻辑层通过多线程并行的向多个MapGIS IGServer服务节点发送REST服务请求获得矢量地图图像，如图1所示。

图1 系统架构图

数据准备层：该层为一个独立的桌面软件，即桌面自动采集工具，主要的功能是对矢量大数据进行样本自动采集，并通过对采集的样本数据进行统计分析，得到矢量大数据渲染时间的计算公式。同时实现了对矢量大数据自动网格化，为后面矢量大数据优化做数据准备工作。

前端应用层：该层为面向地质大数据实时渲染的计算域均衡分解系统的前端部分，主要实现了地图显示、地图基本操作，包括放大，缩小，复位等功能，以及动态的监视并显示每一个服务器节点渲染矢量数据所耗费的时间。

业务逻辑层：该层为本系统的核心部分，主要实现了最佳网格级别查找、可视域范围内的网格单元检索、可视域范围内矢量要素渲染时间统计、可视域范围的均衡划分、通过多线程对矢量数据并行渲染、无缝拼接渲染结果等一系列功能。

GIS服务层：即MapGIS IGServer，该部分主要负责地质矢量大数据服务的管理与服务发布，提供了丰富的GIS功能服务以及完善的二次开发接口。

2 地质矢量大数据实时渲染关键技术

2.1 矢量地图服务时间样本自动采集

矢量地图服务时间样本的自动采集主要实现了在矢量数据的空间范围内随机生成大小、位置各异的矩形范围，利用每次生成的矩形范围对矢量数据进行空间查询，记录当前矩形范围内要素的个数。并且利用以当前随机矩形范围作为矢量地图的渲染范围，从服务器获取图像，同时记录获取图像所耗费的时间。将每一个随机矩形范围对应的要素总数及矢量地图的渲染时间记录到样本文件并保存至指定的文件夹下。实现流程图如2所示。

图2 矢量地图服务时间样本自动采集实现流程图

2.2 矢量地图服务时间公式自动统计

矢量地图服务时间公式自动统计的功能是在矢量地图服务时间样本自动采集完成的基础上进行的，其目的是对采集的样本数据通过一元线性回归的方法，计算出在矢量要素渲染的过程中要素数量与渲染时间之间的线性关系CI=a×count+b中的系数a与截距b，以及R2。其主要的实现思路为：首先从采集的样本数据中获取多组观察值(count1，CI1)，(count2，CI2)，…，然后调用已封装好的最小二乘法计算接口，求出系数a、截距b以及R2。

实现效果图如3所示。

图3 矢量地图服务时间公式自动统计实现效果图

2.3 矢量图层时间网格自动生成

矢量图层时间网格自动生成主要实现了将矢量图层使用四叉树的方式进行网格划分，划分的实现思路为：网格组织方式按照四叉树(4n)组织，第0级一个网格，第1级4个网格，依次类推。算法：首先获取图层的外包矩形范围[Xmin，Ymin，Xmax，Ymax]，然后将其扩展为一个正方形范围[XRmin，YRmin，XRmax，YRmax]，扩展分为两种情况：

第一种，外包矩形宽度(Xmax-Xmin)大于高度(Ymax-Ymin)时，扩展计算公式见式(1)～(4)。

XRmin=Xmin

(1)

(2)

XRmax=Xmax

(3)

(4)

第二种情况，外包矩形宽度(Xmax-Xmin)小于高度(Ymax-Ymin)时，扩展计算公式见式(5)～(8)。

(5)

YRmin=Ymin

(6)

(7)

YRmax=Ymax

(8)

计算第0～n级网格的分辨率，计算公式见式(9)～(11)。

(9)

(10)

…

(11)

计算指定级别l下第r行第c列网格单元的矩形范围[XCmin，YCmin，XCmax，YCmax],规定左下角(XRmin，YRmin)为计算原点，具体计算公式见式(12)～(15)。

XCmin=XRmin+c.Rl

(12)

YCmin=YRmin+r.Rl

(13)

XCmax=XCmin+Rl

(14)

YCmax=YCmin+Rl

(15)

在MapGISLocal数据源指定数据库中创建区图层，将计算的每个网格单元的矩形范围作为几何参数，该网格单元的要素个数、渲染时间、行号、列号等作为属性参数添加到上述新建的区图层中。

2.4 矢量地图服务均匀划分

矢量地图服务均匀划分首先需要在均衡分解模块的全局配置文件中配置结点数量c和一个均衡参数p。如现有部署方式为：1+2模式，1台机器为示例网站部署机器，2台MapGIS IGServer服务器，则c=2。可视域范围采用网格划分方法，则可视域到底在哪个分辨率的网格下划分最均衡至关重要。为计算最佳网格分辨率，设置了均衡参数p。设现有可视域范围为[XVmin，YVmin，XVmax，YVmax],最佳网格分辨率Rp计算公式见式(16)。

(16)

计算得到最佳分辨率Rp之后，通过遍历所有级别的网格的分辨率，得到分辨率最接近Rp的网格的级别数lp。然后利用当前可视域范围对级别数为lp的网格区图层进行空间查询，查询该范围内的所有计算时间网格单元矩形要素。若可视域范围不能完全覆盖正方形范围，则边界求占比，修改可视域范围覆盖不完整的网格单元的渲染时间并且修改网格单元的矩形范围。经过研究分析，将视域范围不能完全覆盖的网格单元按照其所在不同的位置分别计算其边占比，参与边占比计算的网格单元的分布情况如图4所示，具体的算法如下所述。

图4 参与边占比计算的网格单元分布图

为了便于阅读下列公式，此处给出公式相关参数含义：网格单元矩形范围为[XCmin，YCmin，XCmax，YCmax],可视域矩形范围为[XVmin，YVmin，XVmax，YVmax],修改后单元格范围为[XNmin，YNmin，XNmax，YNmax],边占比为Ratio，修改后网格单元的渲染时间为CInew。

左上角单元格计算公式为式(17)。

(17)

修改XCmin为XVmin，YCmax为YVmax。

左边单元格计算公式式(18)。

(18)

修改XCmin为XVmin。

左下角单元格计算公式为式(19)。

(19)

修改XCmin为XVmin，YCmin为YVmin。

下边单元格计算公式为式(20)。

(20)

修改YCmin为YVmin。

右下角单元格计算公式为式(21)。

(21)

修改XCmax为XVmax，YCmin为YVmin。

右边单元格计算公式为式(22)。

(22)

修改XCmax为XVmax。

右上角单元格计算公式为式(23)。

(23)

修改XCmax为XVmax，YCmax为YVmax。

上边单元格计算公式为式(24)。

(24)

修改YCmax为YVmax。

修改渲染时间的公式为式(25)。

CInew=CI×Ratio

(25)

网格数据检索完成之后，遍历所有的网格单元，计算当前可视域范围内网格单元内矢量要素的渲染时间总和。根据服务器节点数量求出平均时间。利用可视域范围内的网格单元进行可视域划分。划分思路可简述为：将可视域范围内的所有的单元格按照其行号进行分行，遍历单元格行，累加单元格渲染时间同时累加每一个单元格的空间范围，直至渲染时间总和大于等于平均渲染时间为止，此处分为两种情况：第一种，渲染时间总和等于平均渲染时间，则直接将空间范围的累加结果分配给第一个服务器；第二种，渲染时间总和大于平均渲染时间，则计算超出平均时间部分与平均时间的占比，根据该占比调整空间范围，然后再将该范围分配给第一个服务器节点，多出的空间范围则累计至下一个服务器节点的空间范围，依此类推。实现流程图如图5所示。

3 结 语

本文在地质大数据爆发式增长的背景下，结合目前在线大数据地图服务在实际应用中的瓶颈。面向地质大数据实时渲染的计算域均衡分解系统采用了一系列的优化方案，包括矢量地图服务时间样本自动采集、矢量地图服务时间工时自动统计、矢量地图网格化、基于网格的可视域划分、多线程并行渲染，渲染结果合并等，在矢量大数据在线服务的响应速率提升方面获得了比较明显的效果。其中基于网格的可视域划分是该系统的核心部分，结合矢量大数据的特征，采用网格单元将矢量要素的渲染时间拆分，再根据服务器节点的数量按照平均渲染时间进行网格单元合并的划分方案，保证了每个服务器节点的渲染时间基本相等，负载均衡效率提升明显。

图5 矢量地图服务均匀划分实现流程图