面向大规模地理矢量线数据的多层级实时可视化技术*

刘泽邦,陈 荦,马梦宇,杨岸然,钟志农,景 宁

(国防科技大学 电子科学学院, 湖南 长沙 410073)

地理矢量数据模型通常以点、线、面等离散对象来描述地理表面空间要素的形状和位置,以及刻画地理空间现象的范围与分布,是全空间地理信息系统描述时空实体的核心数据模型[1]。地理矢量线数据作为地理矢量数据的重要组成部分,通过一系列顺序排列的点来实现对现实世界的表达,被广泛用于道路、行政区划、河流水系、建筑物轮廓等多种线状地物的描述和表达。随着地理空间数据采集技术的快速发展,地理矢量线数据的规模呈现爆炸性增长态势。根据2021年交通运输行业发展统计公报数据显示,截至2021年底全国公路总里程已达5.280 7×107km,涉及的道路矢量数量规模达到了千万级[2];截至2022年,OSM(OpenStreetMap)中的地理矢量线要素的规模已达到亿级[3]。这些体量庞大且不断增长的地理矢量线数据不仅数据量巨大,而且绝大部分的数据仅仅以表格或文本形式呈现,用户需耗费大量的时间和精力进行数据的处理和分析。

空间数据可视化作为空间分析的重要手段,是探索挖掘空间数据的基础[4-5],有助于用户直观地对空间数据进行总结、分析和推理[6-7]。随着网络信息技术的快速发展,基于瓦片金字塔技术[8]的网络瓦片地图服务被提出并用于在互联网浏览地理矢量数据,大规模地理矢量数据的高效可视化成为当前地理信息科学领域的研究热点[9-13]。

在可视过程中,需要遵循经典的可视分析探索流程[14]:既能以鹰眼视角浏览数据的总体分布,又可进行缩放浏览数据的个体细节。为提升可视化效率,单机上的可视化方法通过数据采样、聚合、简化等一系列减少数据规模的方法来实现快速可视化响应[15-20]。一些研究利用GPU的图形处理优势来加速地理矢量数据可视化处理过程[21-23]。矢量瓦片技术作为近期流行的可视化技术[24-25],可支持在客户端直接对数据进行样式渲染。文献[26]实现了大规模卫星点数据的多分辨率的热区分布图展示。Eldawy等[27]实现了散点图、路网图、热力图等多种可视化效果。YU等[28-29]实现了亿级规模空间数据的地图可视化。文献[30-31]采用预先生成部分可视化结果的思想实现了高分辨率的地图实时可视化。然而,上述可视化方法采用数据导向的计算方式,以单个矢量要素作为计算单元进行可视化计算。随着数据规模的急剧增长,会产生地图瓦片数据集生成耗时长,规模大,样式难以交互更改等问题,即使采用高性能并行计算技术也难以满足无缝数据浏览、高效实时渲染等可视化需求[32]。

1 可视化模型描述

当前主流可视化方法以数据为计算导向,如图1所示,其核心思想是以单个矢量要素作为计算单元,将其转换为便于屏幕显示的栅格对象,最后合并得到可视化结果。

图1 传统数据导向的可视化模型Fig.1 Traditional data-driven visualization model

显示导向的可视化模型[33-36]如图2所示。该模型以屏幕像元为计算单元,将依赖矢量要素规模的计算复杂度转换为依赖像元规模的计算复杂度,大幅提升了计算效率。

图2 显示导向的可视化模型Fig.2 Display-driven visualization model

在可视浏览时,可视瓦片的属性不尽相同,如显示层级或高或低,要素数量或多或少,导致绘制效率差异较大。为实现各显示层级上的实时可视化,提出面向多层级瓦片绘制的自适应可视化模型。如图3所示,模型对瓦片绘制任务进行自适应判定,选择最优方式绘制可视结果。当采用显示导向可视化算法,进行像元生成判定;当采用数据导向可视化算法,检索瓦片内矢量要素进行栅格化。在模型中,地理矢量数据需转换成特定的空间索引用以支撑多层级瓦片绘制,该索引一方面能为瓦片绘制算法的自适应选取提供支撑,另一方面能兼顾显示导向与数据导向两种可视化算法的计算需求。最后,在自适应判定器中设定判定因子,根据瓦片绘制任务的判定结果来选取可视化算法。判定因子包括:瓦片显示层级,即显示层级越低,则线要素越多,可视化耗时长;瓦片内矢量线要素数量,即线要素越多,可视化耗时越长;瓦片内矢量线要素总长度,即线要素总长度越长,可视化耗时越长。

图3 面向多层级瓦片绘制的自适应可视化模型Fig.3 Adaptive visualization model for multilevel tile drawing

2 像元四叉R树索引

2.1 索引结构设计

在面向多层级瓦片绘制的自适应可视化模型中,为绘制任意显示层级上的一张瓦片图,数据导向的可视化方法需快速检索与该瓦片对应空间范围拓扑相交的矢量线要素;显示导向的可视化方法需要判断像元对应空间范围是否与任一矢量线要素拓扑相交。

为实现快速可视化所需的拓扑相交和数据检索的计算需求,设计了一种PQR树空间索引。如图4所示,PQR树是由多个索引结点连接的类四叉树结构,有别于常规四叉树以数据集最小外包框作为根结点索引范围,PQR树将全球空间范围作为根结点索引范围,再不断向下进行递归四分得到索引结点,索引结点共五种类型:根结点、左上结点、右上结点、左下结点、右下结点。在索引结点中,记录结点的空间范围、地理编码、结点属性与结点指针等信息。当达到预先设定的最高层级后则不再进行划分,在PQR树的最高层级上,每个索引结点连接一个R树,对与索引结点索引范围拓扑相交的矢量线要素进行高效组织。

图4 PQR树索引结构设计Fig.4 PQR-tree index layout

PQR树索引结点说明如下:①基于Geohash的编码方法来获取索引结点编码。如图4所示,从根结点向下编码,每向下递归划分一层:左上、右上、左下、右下结点的编码分别为“00”,“10”,“01”和“11”,不同层级结点的编码长度不一致,最终得到二进制形式的编码作为索引结点编码。②对于矢量线数据,选取要素数量count、显示层级level、要素长度length三个判定因子作为结点属性。③结点指针包括四个子结点指针LU-ptr、RU-ptr、LB-ptr、RB-ptr和一个R树指针R-ptr。非叶子结点的四个子结点指针分别指向该结点索引范围四分后的四个子区域所对应的子索引结点,R树指针为空;叶子结点的四个子结点指针均为空,R树指针指向该索引结点所连接的R树。

PQR树的索引结点的索引范围与瓦片/像素对应空间范围是唯一映射的,这一特性使其可同时满足自适应判定、数据导向与显示导向可视化方法的计算需求:①只需快速找到瓦片对应的索引结点,并对其中记录的判定因子进行比较,即可完成自适应判定;②对于瓦片对应的索引结点,对该结点连接的R树执行空间检索获取瓦片对应空间范围内的矢量要素,即可完成数据导向可视化方法的计算需求;③对于瓦片对应的索引结点,只需从该结点向下判断像素对应的索引结点是否存在,即可完成显示导向可视化方法的计算需求。

2.2 索引构建

地理矢量线数据的PQR树索引构建算法的核心思想是从根结点“自顶而下”地将矢量线要素与全球地理空间范围递归四分的结果进行拓扑相交判断(“位于”和“叠置”)来生成索引结点,在达到最高层级后,将矢量线要素组织生成R树结构。由于像元对应空间范围是由瓦片对应空间范围递归划分八次得到,为满足显示导向可视化方法的计算需求,PQR树的层级不低于8。

在创建索引结点时,最常采用近似表达的方式将矢量线要素的最小外包框(minimum boundary box, MBB)与空间范围框进行判定。考虑当在低显示层级时,矢量线要素完全位于空间范围框内,而在高显示层级时,矢量要素只与部分空间范围框相交。为提高计算效率,将线要素表示为双层外包框结构,如图5所示,其包含了线要素的MBB和线要素中每条线段的MBB。在相交判断时,依次采用线要素的MBB和每条子段的MBB用于计算,两层过滤后再判断每条子段与空间范围框是否空间拓扑相交。

图5 地理矢量线要素的双层级最小外包框结构Fig.5 Double MBB structure of geographic vector linestring features

PQR树索引构建流程包括:①创建根结点,对于任一地理矢量数据集,根结点的空间范围均为全球地理空间范围。②逐一遍历数据集中地理矢量线要素,将其表示为双层级外包框结构。③将根结点的空间范围划分为四个子空间,分别与线要素的MBB进行空间关系判定,若线要素位于某个子空间内,则创建相应的索引结点,并从该索引结点继续向下进行递归判定,直至生成包含要素MBB的最小索引结点。④依次选取要素的每条子段的MBB,将其分别与最小索引结点的空间范围划分后所得的四个子空间进行空间关系判定。⑤若子段的MBB位于某个子空间内,则创建相应的索引结点;若子段的MBB与某些子空间相叠置,则再精细判定子段与这些子空间是否相叠置来创建相应的索引结点;继续从新创建的索引结点向下递归判定,直至设定的最高层级。⑥分别将要素子段插入与其相交的索引结点所连接的R树中。⑦当所有的矢量线要素判定完毕后,构建得到最终的PQR树索引。

3 基于像元四叉R树的自适应可视化

在面向多层级瓦片绘制的自适应可视化模型中,为快速绘制各个显示层级上的瓦片图,设计基于PQR树的自适应可视化(PQR-tree-based adaptive visualization, PAV)算法。PAV对判定因子的判定流程进行了设定:当采用显示导向可视化方法,核心是判定像元对应的索引结点在PQR树中是否存在;当采用数据导向可视化方法,核心是从PQR树中快速检索瓦片对应空间范围内的矢量线要素。

根据2.2节,对线数据集构建PQR树后,与线要素存在空间拓扑相交的空间范围框会生成对应的索引结点。所以对于待绘制的瓦片任务,需要判断瓦片对应的索引结点在PQR树中是否存在来决定瓦片是否需要被绘制。同时,在显示导向的可视化绘制时,也需要判断像元对应的索引结点在PQR树中是否存在来决定像元值是否需要被生成。基于此,设计了索引结点查找函数SEARCHNODE,如算法1所示,“自顶而下”地从设定索引结点向下查找对应的索引结点。函数输入包括开始查询的索引结点sNode和瓦片/像元对应空间范围框的地理编码rCode,当查询到对应的索引结点时进行输出,反之输出为空。算法主要步骤为:①逐个遍历初始结点的子结点指针,判断子结点cNode是否存在;②当子结点存在时,获取cNode的地理编码cCode;③判断cCode与rCode长度是否相等;若两者不等,取两者的前缀对应位数进行异或操作,称为前缀匹配;④若操作结果为0,将cNode与rCode作为函数输入,递归执行索引结点生成函数;⑤直至当rCode与cCode完全相等,输出cNode为与空间范围框相对应的索引结点,其余情况下均输出为空。

算法1 索引结点查找函数Alg.1 Searching function of index node

当瓦片对应的索引结点存在时,则说明瓦片需要被绘制。此时,PAV获取到该瓦片对应的索引结点中的属性信息(显示层级level、线要素数量count、线要素总长度length),并与三个判定因子的阈值进行比较,具体的判定流程如下:①判断层级信息level与δlevel的大小,当level大于δlevel时,则选择数据导向的可视化方法;反之继续判断。②判断count与δcount的大小,当count大于δcount时,选择显示导向的可视化方法;反之继续判断。③判断length与δlength的大小,当length大于δlength时,选择显示导向的像素生成策略;反之选择数据导向的可视化方法。

PAV的具体实现主要分为五步:①调用算法1的索引结点查找函数,判断待绘制瓦片对应的索引结点在索引中是否存在,若不存在则直接输出空白的瓦片图;反之继续执行。②获取索引结点的结点属性信息(level,count,length),与设定的判定因子阈值进行判定,自适应地选取可视化方法。③若选取显示导向的可视化方法,则遍历瓦片的每个像元,继续调用算法1的索引结点查找函数,判断像元对应的索引结点是否存在,若存在则直接在瓦片图上生成像元值。④若选取数据导向的可视化方法,则对该索引结点连接的R树执行空间范围检索操作,得到地理矢量线要素,并逐一进行栅格化处理生成像元值。⑤最后输出得到瓦片图。

PAV的优势在于:采用显示导向的可视化方法时,查找像素对应的索引结点是否存在只有唯一的判断路径,采用简单的二进制运算即可快速生成瓦片图;采用数据导向的可视化方法时,在查找瓦片对应的索引结点是否存在同样只有唯一的判断路径,进而快速检索到矢量要素并进行栅格化生成瓦片图。最终可实现各个显示层级上瓦片图的快速绘制。

4 实验与分析

实验数据集如表1所示,来自开源志愿者网络地图服务平台OSM。从L1至L7,数据规模逐渐增大。L7包含OSM全球所有线要素数据,其中线要素的规模达到了10亿级别。实验环境:内存为256 GB,核数为32的2.6 GHz的Inter处理器,且处理系统为Ubuntu 20.04。

表1 实验数据集Tab.1 Datasets used in the experiments

为确立PAV中判定因子的阈值,实验一选取L1～L6作为测试数据集,对比显示导向的可视化方法(display-driven visualization, DiDV)和数据导向的可视化方法(data-driven visualization, DaDV)在不同的判定因子取值下的瓦片绘制性能,为PAV中的自适应判定提供依据。实验二在所有线数据集上分析PAV算法的可视化性能,其中L7作为测试集用以验证算法的泛化能力,最终验证PAV算法可以支撑大规模地理矢量数据集的多层级实时交互式可视化。实验三将PAV与当前主流可视化方法进行比较,验证PAV算法的优势。在实验中所有数据集均被预先转换成PQR索引以支撑PAV算法。同时,由于每个瓦片的绘制任务间相互独立,故采用多进程消息传递框架(message passing interface, MPI)与OpenMp混合并行的方式来提升瓦片绘制的性能,在所有实验中均启动了8个MPI进程(每个进程中包含4个OpenMp线程)用于绘制瓦片。

4.1 实验一:判定因子阈值设定

本实验用于获取判定因子的阈值。采用DiDV与DaDV绘制瓦片图,统计每张瓦片判定因子的值,评估瓦片的显示层级、线要素数量、线要素总长度对两种方法性能的影响。

为分析显示层级的影响,在各显示层级随机选取2 000张非空瓦片,计算两种方法在各显示层级上的绘制速率,实验结果如图6所示,当显示层级较低时,DiDV的可视化性能明显优于DaDV,并且当数据规模较大时的性能优势更明显,在L4数据集上当显示层级为6时,DiDV的绘制速率为DaDV的35倍;随着显示层级升高,DaDV的效率显著提升并开始超越DiDV,在L6数据集上的最大绘制速率达到563.8张/s。综上可知,当显示层级较低时应选择DiDV,反之选择DaDV。在设定阈值时,将完全包含数据集的最小瓦片层级作为其空间尺度值,例如全球空间的数据集的空间尺度则为0。当数据集的空间尺度不同时,显示层级阈值δlevel也不同,对于L1～L6,δlevel

(a) L1

(b) L2

(d) L4

(e) L5

(f) L6图6 瓦片绘制速率随显示层级变化的对比Fig.6 Comparison of tile drawing speed at different zoom levels

的值比空间尺度的值大8至11,为满足所有数据集在各个显示层级上的快速可视化,以数据集的空间尺度值+11后的值设为PAV算法中的δlevel。

为分析瓦片空间范围内线要素数量的影响,设定六个不同的范围区间:0～400、400～800、800～1 200、1 200～1 600、1 600～2 000、>2 000,对于每个范围区间随机选取2 000张非空瓦片,计算DiDV和DaDV的瓦片绘制速率,实验结果如图7所示。DaDV在线要素数量较小时表现出优于DiDV的可视化性能,然而其性能随线要素数量的增长而显著下降,而DiDV的瓦片绘制速率维持着相对稳定的水平。对于数据集L1～L6,实验结果显示在线要素数为400～800范围内DiDV的瓦片绘制速率开始超越DaDV,由此以800作为PAV算法中的线要素数量阈值δcount。

(a) L1

(b) L2

(c) L3

(d) L4

(e) L5

(f) L6图7 瓦片绘制速率随线要素数量变化的对比Fig.7 Comparison of tile drawing speed varies with the number of linestring features

为分析瓦片空间范围内线要素的影响,设定六个不同的范围区间:0～200 km、200～400 km、400～600 km、600～800 km、800～1 000 km、>1 000 km。对于每个范围区间随机选取2 000张非空瓦片,计算DiDV和DaDV的瓦片绘制速率,实验结果如图8所示,DaDV在线要素总长度较小时表现出优于DiDV的可视化性能,然而其性能随线要素总长度的增长而显著下降,而DiDV的瓦片绘制速率维持着相对稳定的水平。对于数据集L1～L6,实验结果显示在线要素总长度为200～400 km的范围区间上DiDV的瓦片绘制速率开始超越DaDV,由此以400 km作为PAV算法中的线要素总长度阈值δlength。

(a) L1

(b) L2

(c) L3

(d) L4

(e) L5

(f) L6图8 瓦片绘制速率随线要素总长度变化的对比Fig.8 Comparison of tile drawing speed varies with the total length of linestring features

4.2 实验二:支撑多层级实时可视化性能

为验证PAV算法具备支撑大规模地理矢量线数据的多层级实时可视化的能力,将实验一确定的阈值(δlevel=11、δcount=800、δlength=400 km)代入PAV算法中,对于所有数据集L1～L7,在0～20显示层级上随机选择5 000张非空瓦片,分别统计生成5 000张瓦片的总耗时和生成每张可视化瓦片的绘制时间。实验结果见图9,对于所有数据集,PAV算法均维持着相对稳定高效的可视化性能,瓦片绘制速率最低也能达到366.5张/s,远高于实时浏览数据所需的瓦片请求速率(通常不超过100张/s)。在10亿规模测试数据集L7上,瓦片生成总耗时仅为13.6 s,瓦片绘制速率达到366.5张/s。同时,从L1～L7数据集规模显著增长,而算法的瓦片生成总耗时并无明显的增长,且瓦片绘制速率并无显著的下降趋势,说明其具有对数据集规模不敏感的特性,非常适用于对大规模地理矢量线数据的多层级实时可视化。图10为生成5 000张瓦片的时间分布箱型图,其中“”表示瓦片的平均生成时间,从结果可知,对于所有数据集,绝大部分的瓦片在30 ms内可完成绘制,在10亿规模测试数据集L7上,瓦片的最长生成时间仅为44 ms。假设当前浏览器发送了100个瓦片请求,所有任务通过13轮(100÷8)绘制,最差情况也能在0.57 s((13×44) ms)内完成任务,完全可以满足实时可视化的需求。综上,结合图9和图10结果分析可知,PAV算法可以高效地支撑10亿级规模的地理矢量线数据的多层级实时可视化。

图9 生成5 000张瓦片的总耗时和绘制速率Fig.9 Total time and tile drawing speed of drawing 5 000 tiles

图10 生成5 000张瓦片的时间分布箱型图Fig.10 Box diagram of time distribution of drawing 5 000 tiles

4.3 实验三:与主流可视化方法的性能对比

为验证PAV算法的优势,将其与当前两个主流的空间数据可视化方法GeoSparkViz和HadoopViz进行比较,对于所有数据集,分别采用3种方法生成0～8显示层级上的所有可视化瓦片,得到三种方法的总耗时对比,如图11所示,对于所有数据集,GeoSparkViz生成瓦片的耗时少于HadoopViz,GeoSparkViz的瓦片生成性能相较HadoopViz更优,然而随着数据规模增长,这两种方法的瓦片生成性能急剧下降。而PAV生成可视化瓦片的耗时均远少于GeoSparkViz与HadoopViz,表现出优异的可视化性能,对于体量超过亿级的数据集L4～L7,PAV的瓦片生成耗时分别仅为GeoSparkViz的13.2%、11.6%、8.1%、7.5%,且数据规模越大,PAV相较GeoSparkViz的瓦片生成性能优势越明显。从L1～L7,随着数据规模的增长,PAV的耗时并未出现大幅增长,验证了PAV具有对数据规模不敏感的特性,能够支撑大规模地理矢量线数据的多层级实时可视化。

图11 生成可视化瓦片的总耗时Fig.11 Total time to render all tiles in zoom levels 0～8

为进一步对比三种方法在不同显示层级上的可视化性能,图12对比了三种方法在亿级规模数据集L4～L7上第1、3、5、7、9显示层级上的瓦片绘制速率,从实验结果可知:GeoSparkViz和HadoopViz在较高的显示层级上具有较高的瓦片绘制速率,其中GeoSparkViz的瓦片绘制速率优于HadoopViz,GeoSparkViz在L4数据集的第9显示层级上的瓦片绘制速率能达到269.77张/s。然而,随着显示层级的降低,GeoSparkViz与HadoopViz的瓦片绘制速率大幅下降,按照可视化浏览从整体到精细的原则,GeoSparkViz与HadoopViz无法实现大规模地理矢量线数据的多层级实时可视化。而PAV在所有数据集的各个显示层级上均具有较高的瓦片绘制速率,特别当显示层级较低时,其瓦片绘制速率最差在L7数据集上也能达到115.29张/s,大幅优于GeoSparkViz(0.01张/s)与HadoopViz(0.001张/s)的可视化性能。

(a) L4

(b) L5

(c) L6

(d) L7图12 各显示层级上的瓦片生成速率对比Fig.12 Comparison of tile drawing speed on each zoom level

从上述结果可知:对于反映可视化效率的两个主要指标瓦片生成总耗时和各显示层级的瓦片绘制速率,PAV均大幅优于主流可视化方法GeoSparkViz和HadoopViz。原因在于:GeoSparkViz和HadoopViz采用数据导向的计算方式,从L1至L7,随着数据规模的增长,这两种方法的瓦片生成总时间必然会大幅增长;并且对于同一数据集,当在较高的显示层级上,瓦片内的矢量线要素的规模较小,所以GeoSparkViz和HadoopViz具备相对较高的瓦片绘制速率,随着显示层级的降低,数据规模急剧增长,某张瓦片内的矢量线要素可能达到千万甚至上亿量级,从而使得GeoSparkViz和HadoopViz在低显示层级上的瓦片绘制速率显著下降。PAV采用显示导向与数据导向相结合的自适应可视化方法,当数据规模较小或显示层级较高时使用数据导向的可视化方法能快速绘制瓦片,反之采用显示导向的可视化方法直接计算像元点的方式能大幅提升可视化效率,由此在较低显示层级上仍具有较高的瓦片绘制速率,所以PAV在所有数据集的各个显示层级上均具有较高的瓦片绘制速率,并且总的瓦片生成总耗时大幅降低。综上所述,PAV相较现有主流的数据导向的可视化方法GeoSparkViz与HadoopViz具有明显优势,能够支撑大规模地理矢量线数据的多层级实时可视化。

5 结论

针对大规模地理矢量数据可视化这一当前GIS领域研究的热点问题,以大规模地理矢量线数据为研究对象提出高效的多层级可视化技术。该技术采用瓦片金字塔模型来展示地理矢量线数据,以实现高效的瓦片绘制性能为计算目标,提出面向多层级瓦片绘制的自适应可视化模型,在模型中通过自适应判定来选择最优的显示导向或数据导向的可视化方法,同时设计空间索引结构和可视化算法来支撑模型中涉及的数据组织和可视化流程。实验证明该技术在各个显示层级上均具有高效的瓦片绘制性能,可在单机条件下支撑亿级规模地理矢量线数据的多层级实时可视化,在空间数据的可视化探索上具有较好的应用前景。