2. 武汉大学 测绘遥感信息工程国家重点实验室，湖北 武汉 430079
2. State Key Laboratory for Information Engineering in Surveying, Mapping and Remote Sensing; Wuhan University, Wuhan 430079, China
1 引 言
近年来，基于矢量图形信息的三维可视化和基于实景影像信息的三维可视化技术应用越来越广泛，以LiDAR技术为代表的空间数据获取手段正支撑其快速发展[1]。地面激光雷达已经成为在近地空间获取精细三维图形信息的主要手段之一。地面激光雷达获取的数据是按阵列方式存放的、带反射强度的三维点坐标。在此定义按阵列存放的三维点坐标序列为点图像。经过配准的多视(站)点图像称为点图像模型。点图像模型具有数据量大(海量性)、按扫描线排列(栅格性)、数据未经任何处理(原始性)等特点。有学者利用这些特性进行了很有意义的工作，采用球面有向搜索和球面LOP局部优化算法对单站点图像进行不规则三角网模型的构建，提高了构网时间效率，避免了由于数据噪声引起的构网错误和空洞现象[2]；采用一种基于距离依赖的球面重采样方法进行点云数据的精简工作，有效地降低了配准操作的时间成本[3]。点图像模型能够从整体上表达地物的真实形态，是必须永久保存的原始数据。但其数据量很大，且分布不规则，单站数据就有几百万至几千万个三维点。一般地面激光雷达站数较多，如故宫太和殿上下里外扫描达上百站。如何对这种海量数据进行管理是当前面临的最大问题。现有的激光雷达数据处理软件一般将其转换为非结构化随机分布的散乱点云，对海量点云数据进行分块处理，只选取有用的数据，丢弃其他数据，失去了数据的原始性和栅格性。特别是对于故宫建筑物这样的古建筑，其具有对象部件多样化、数据查询操作要求效率高等特点，空间数据的组织、存储与管理具有较大难度。
利用数据库来存储海量空间和属性数据是发展方向，空间数据和属性数据可采用对象关系数据库技术来存储[4, 5]；同时，有针对性地设计简单高效的空间索引也是关键问题之一。纯粹二维索引效率很低，因为三维空间分布的点云在不同平面单元里差异很大，传统三维索引一般又没有顾及不同的对象。目前适合点云的三维空间索引方法主要有格网、四叉树、八叉树、k-d树[5, 6, 7]、R树及其变种[8, 9, 10]和一些很具代表性的混合索引[11, 12]等。这些索引中很多是点对象和空间实体的最小外包矩形建立的，但其索引粒度基本上都索引到了单个点对象，不适用于数据量特别大的点云三维空间索引。现有的空间索引往往支持内存或外存的细节层次技术，在生成空间索引和细节层次数据时通常是复制点云数据备份，采用“以空间换时间”的方法来选取所需层次的数据进行处理与可视化。目前还没有一种空间索引方法明显优于其他方法[5, 7]，往往需要根据具体的应用来设计合适的空间索引。
为了便于整体点云模型的场景表达，提高海量点云数据处理和可视化的效率，保证原始数据信息不被破坏和丢失，在充分利用点图像模型高信息量的基础上，针对上述问题提供一套新的海量点图像数据库管理算法，为点云数据的后处理 工作奠定基础。对于单站点图像模型提出一种二维与三维集成的空间索引进行数据组织，对于多站点图像数据采用该集成空间索引与三维R树空间索引的混合索引来组织。将空间索引与点云数据分别组织，只生成一份点云数据，直接利用空间索引树结构生成细节层次模型，通过建立一种可视化调度策略针对不同视距采用不同的细节层次可视化。
2 海量点图像数据管理算法原理
针对地面激光雷达扫描获取的单站和多站点图像数据，本文采用的算法主要包含以下几个方面的内容：二维与三维集成空间索引的设计、混合空间索引的构建与数据存储和数据细节层次调度，算法原理流程如图1所示;利用点云数据的栅格性在单站数据读取的同时采用自上而下再自下而上的方式建立二维四叉树和三维最小外包盒树空间索引结构(quad-MBB，简称QMBB树)，针对多站点图像数据则采用3D-R树空间索引管理各单站数据的根节点MBB，将海量原始数据分别组织成序列化空间索引数据和经压缩后的点云数据，并将组织好的两类数据和其他属性数据存储到大型商用数据库中；基于自适应细节层次技术和数据调度策略实现单站或多站点图像数据的数据库调度与可视化。
|
| 图 1 算法原理总流程图 Fig. 1 The total flow of the algorithm principle |
利用原始扫描点云的栅格性，按点云的行数和列数先以自上而下的方式构建基于二维四叉树空间索引，每次按照实际扫描的横向步进角度和纵向步进角度分别进行二分，对点图像数据进行二维分割，每个分割都形成一个二维最小外包矩形(minimal bounding rectangle,MBR)，如图2左列所示，表示基于激光雷达二维深度图的四叉树分割示意图。根据每个MBR对应的点云在三维空间中计算每个MBR对应的最小外包盒(minimal bounding box,MBB),以二维四叉树的形式构建3D-MBB树，如图2右列所示。将二维四叉树与3D-MBB树集成起来构成2D-3D空间索引QMBB树，它兼具二维和三维特征，是一种集成空间索引。
|
| 图 2 QMBB树空间索引的数据分割 Fig. 2 Data segmentation using QMBB-tree spatial index |
QMBB树空间索引的索引粒度为一个2D-MBR和对应的3D-MBB，实际是索引一个点集。普通四叉树节点只含二维信息，而QMBB树则同时具备二维属性和三维属性，可进行二维和三维空间查询，分割后的点图像数据块没有重叠点，而且QMBB树非叶子节点内的采样点均来自原始数据，没有经过插值运算。QMBB树三维空间索引的节点数据结构如图3所示。它由指向子节点的标识children和指向父节点的标识pParent、三维坐标偏移CoordOffset和三维点的ID标识偏移IDOffset、三维点数据ID标识IdList、节点的最小外包矩形体MBBox、节点的行列位置最小外包矩形MBR、节点的名称编码strCode等组成。节点的名称编码采用线性四叉树编码方法，每个节点名称编码的字符串长度代表它在QMBB空间索引树中的深度。
|
| 图 3 QMBB树空间索引节点数据结构 Fig. 3 Node data structure of Q-MBB tree spatial index |
QMBB树空间索引的构建流程为，在每个点坐标读取的过程中依次定义点的ID值并根据ID值将其插入QMBB树索引的叶子节点中，如果某一个节点内存储的点数超过阈值就将该节点进行分裂，直至所有的空间点均插入QMBB树。QMBB树构建完毕后计算每个叶子节点的MBB，利用均匀采样的方法以自下而上的方式填充空间索引树的非叶子节点数据，直到根节点，该QMBB树的非叶子节点也只存储坐标数据的索引ID。其构建方法如图4所示。首先获取单视点图像的行号列号，根据行列号530×789从上至下构建QMBB树，接下来将每个三维点数据按顺序插入对应叶子节点，如红虚线的点插入“第三层次叶子节点1”，蓝虚线的点插入“第三层次叶子节点2”，当每个有效三维点均插入QMBB树叶子节点后，再从叶子节点自下而上依次重采样填充各层次非叶子节点的点云数据，如图绿实线箭头所示。这样QMBB树在数据读取的过程中一次性快速实时生成，点图像在逻辑上被分割成多分辨率细节层次数据块。
|
| 图 4 单站点图像QMBB树构建及数据插入方法示意图 Fig. 4 QMBB-tree index building of single-station grid point-cloud |
QMBB树空间索引结构建立后，多分辨率点云数据的层次通过QMBB树的深度来表达，QMBB树的每一层代表不同分辨率的点云数据，分辨率的大小取决于QMBB树空间索引构建时点云的均匀采样率。QMBB树空间索引结构中只存储三维点的ID号，作为真实点坐标及灰度纹理信息的指针，而且点坐标数据只存储单个拷贝，检索时通过ID号换算为偏移指针查询数据，从而避免多层次数据的存储冗余。利用3D-R树空间索引管理每站点图像数据，3D-R树的每个叶子节点均存储QMBB树的根节点MBB，同时3D-R树的每个节点也存储各单站点图像的均匀采样数据，采用二级空间混合索引的方式来管理所有点图像模型数据。以深度优先遍历的方式将空间索引序列化并入库，将细节层次点云数据和一些必要的属性数据，如扫描时间、扫描地点、扫描人员、扫描仪型号等分别存储到已设计好的数据库表中，按空间对象、点图像数据对象和存储的数据块表3级模式来存储海量点图像数据，再次检索序列化索引数据后可快速恢复该点云的空间索引结构。受文献[13]的双字节存储法启发，对于三维浮点坐标值采用空间位置量化的方法进行压缩，并将点的ID、灰度和纹理等属性信息分别进行压缩，压缩率可达37.5%，便于传输和现有微机显存一次性绘制较大规模的数据。
3 点图像数据可视化近几年国际国内研究较多的是利用点模型[14, 15, 16]来进行大规模点云的可视化工作。它能更好地表达形状复杂且不规则的物体，能够兼顾绘制速度与绘制质量[17, 18, 19]，但其绘制的点数最多也只有千万级。而且在点模型的数据预处理阶段往往需要很大的时间和空间消耗来计算模型表面的几何属性，如求解主曲率、主方向，计算法矢和法锥面等，而原始的扫描数据往往是没有这些几何信息的，需要进行求解与计算，且花费时间大都比较长，而且超过一定容量的数据不能被预处理。如QSplat软件[14]的.qs格式转换工具Laser Splat就对数据量大小有限制，又如XGRT软件[20]在构建基于外存的数据文件时生成的二进制文件约为原始数据文件的3倍，这些都降低了海量点云后处理工作的效率。
针对可视化时海量点图像数据的可视化问题，本算法采用的方案是首先对单站或多站点云的MBB进行视锥裁剪，然后采用基于视点远近的细节层次(LOD)技术来选择需要显示的指定级别数据，最后将选出的MBB进行遮挡剔除，确定实际需要显示到屏幕上的点云数据。在内存内设置检索查找表，判断哪些数据块不用再到数据库内调度，只检索不在内存中的坐标数据。通过选出的MBB获取点ID，再根据ID值在数据库中直接检索对应的真实点图像模型数据参与可视化，提高了时间和空间利用率，可视化数据调度选择流程如图5所示。
|
| 图 5 栅格点云模型可视化 Fig. 5 Visualization of grid point-cloud model |
利用3D-R树和QMBB树空间索引结构中节点的MBB与三维场景视锥体进行视域裁剪，判断MBB与视锥体是否相交，排除不在视锥体内的节点。进行单站或多站之间的细节层次数据的选择，在图6所示的自适应细节层次选取方案示意图中，对于每个存在视锥体内的MBB的8个顶点，逐个计算屏幕视点到各个顶点的距离，记录距离的最大值MaxDist和最小值MinDist以及对应的2个顶点，将排除最大最小两个顶点后的6个顶点分别投影到视锥体的近平面上，计算投影的六边形(q1′q6′q3′q7′q4′q5′)在近平面上的投影面积Area,定义当前每个点在屏幕上显示的大小cursize，将Area的值乘以cursize和定权因子a即为实际投影面积ScreenArea。其中定权因子a的意义为点集中所有点的平均显示大小的权值，由于点的固定大小在显示的时候会有部分重叠或完全叠置的现象，用a来调节显示的平均大小，提高LOD级别的判断准确性。最终采用ScreenArea与MinDist的比值定义比例因子RateofArea2Dist来保守估算细节层次级别数据。利用不同细节层次下RateofArea2Dist取值范围的不同来自适应确定需要显示到计算机屏幕上的数据，比例因子越小，则对应较为粗略的层次数据；比例因子越大，则对应较为精细的层次数据。本文实现的细节层次模型的自适应性不仅表现在依据视点远近关系判断LOD层次上，而且表现在依据单个MBB的视距分别判断其各自的LOD层次，对不同的MBB树节点自动地判定其在当前三维场景中的LOD层次上，也即不仅是判断整个点图像数据的LOD层次，而且是针对每个点图像分割数据块分别判定其LOD层次，不同的MBB树节点内的点图像分块数据具有不同的LOD层次。
|
| 图 6 单站时细节层次选取 Fig. 6 single-station LOD selection |
层次数据选定后，根据整个场景的视点和MBB间的分布情况判断多站间哪些MBB内的数据应该保留下来参与绘制，统计需要绘制的点的总数，如果超过设定阈值，则将各个数据块的LOD级别降为较粗略的一层数据显示，直到点数满足阈值要求为止，保证内存中需要绘制点的数量和内存消耗量在合理的水平，最终达到有较高的显示帧率和较好的用户体验。记录保留下来参与绘制的MBB节点存储的三维坐标点的ID值，利用ID值从数据库内查询点云数据，最后将一定个数的点数据置于GPU中进行绘制。
本文采用配置为Intel(R) Core(TM)2 CPU6600 2.4GHz,3.00GB内存、GeForce8500GT显卡的微机，以Oracle11g作为数据库管理平台，利用数据库系统设计工具PowerDesigner设计数据库对象模型、数据库表，使用PL/SQL语言进行数据库建库和管理，利用Microsoft C#/C++.NET 2005实现了具有快速浏览和检索查询功能的交互式可视化管理原型系统。
4 试验与分析图7是点图像(单视)的LOD显示效果，扫描后得到的数据文件大小为18.52 M，含有彩色有效点个数为358 062，由于三维激光扫描仪的获取深度图像的特性，文件内还含有以某种格式记录的空间无效点，但其不参与可视化。每次屏幕内显示的点数约为150 000左右，在表2所示配置的计算机上采用OpenGL三维图形库实现了本文的算法，平均显示帧率约为120。图7中按箭头所示方向点云的LOD级别逐渐减低，表示数据细节层次分割粒度的MBB线框随之变大，但如图所见每一帧上点云可视化效果并没有减低。
|
| 图 7 单视点图像LOD显示效果 Fig. 7 LOD display of grid point-cloud |
图8是点图像数模型的LOD显示效果，数据文件大小为269.14M，含有12块栅格数据，彩色有效点总数为5 347 765，每次屏幕内显示的点数约为300 000左右，在同样配置的计算机上平均显示帧率约为60左右。该图则是展示点云的LOD级别逐渐增大的过程，红色的线框表示内部分块节点数据的MBB，在不同的分辨率下查询出的MBB的大小不同，分辨率越高，MBB线框越密集，LOD级别越高，数据显示越精细。需要指出的是，本文采用的帧率计算方法为每次强制执行显卡缓冲区内数据，而不用等到显卡缓冲区满再执行指令，采用此方法计算的帧率不受计算机屏幕刷新率的限制，更能反映真实的显示速度。
|
| 图 8 点图像模型数据LOD显示效果 Fig. 8 LOD display of grid point-cloud model |
表1显示利用本文算法调度和绘制海量点图像模型的效率。根据LOD选择后计算获取的点集的数量的大小，不同总点数绘制的帧率也有所不同，在保证基本绘制质量的前提下，经多次试验统计，其中绘制速度达到或超过10 000 000点/s(实际绘制点数×实时帧率)，比文献[19]所用的方法快2倍以上，较文献[17]报告的绘制速度快6倍以上。LOD数据检索总时间包括多站视锥裁剪、LOD选择等点的ID选定过程耗用时间之和。
本文方法能同时调度多站点图像数据进行绘制，如果某一层次的多站数据绘制时点的总数超过了某一阈值，则重新从3D-R树和QMBB树混合空间索引中检索上一层次的数据显示，确保每次屏幕显示的点的总数不超过规定的阈值，基本不影响绘制质量的前提下保持较高的帧率。如图9所示，左图为故宫太和门整体点图像模型，交换格式文件总大小为1.54 G，它由241 596 602个有效点组成，实际绘制的点数为10 711 195个；右图为中和殿整体点图像模型，交换格式文件总大小为4.74 G，它由2 038 308 028个有效点组成，实际绘制的点数为3 940 690个，绘制帧率平均保持在4 fps左右,能够基本满足流畅显示的要求。
| 交换格式文件
总大小/M |
数据块数 | 点的行列数 | 有效点总数 | 实际绘制点的个数 | 显示帧率
/fps |
LOD数据检索
总时间/s |
| 18.52 | 1 | 530×789 | 358 062 | 40 000~200 000 | 110~458 | 0.043 18 |
| 269.14 | 12 | 121×95
2963×607 … |
5 347 765 | 200 000~350 000 | 30~80 | 0.787 31 |
| 421.60 | 7 | 1682×865
568×2760 … |
7 774 564 | 550 000~1 340 000 | 14~50 | 0.748 63 |
|
| 图 9 多站点图像数据整体调度与可视化 Fig. 9 The whole visualization of multi-station grid point-cloud models |
本文的所提出的算法可以在海量点图像数据读取的同时快速建立QMBB树集成空间索引，利用3D-R树和QMBB树多级混合索引方便地存取自组织的空间索引数据和点云数据，避免了外存的不必要LOD数据备份与采样点冗余，节约外存和内存的空间，数据结构简单有效，可视化操作时查询效率较高，同时利用大型商用数据库来存储数据保证了数据的安全性与并发性，具有较高的实用价值。
下一步研究的主要方向有：构建基于网络的应用平台，实现实时远程数据处理应用；在绘制质量与绘制效率上考虑新型的GPU可编程绘制手段和自适应点模型绘制模式，提高可视化质量与效率；利用本文的空间数据模型和数据库存储的数据，开发数据处理与可视化的算法与应用系统，为大型LiDAR数据处理与管理平台提供技术支撑。
| [1] | LI Deren. 3D Visualization of Geospatial Information: Graphics Based or Imagery Based[J]. Acta Geodaetica et Cartographica Sinic. 2010, 39(2): 111-114. (李德仁. 论地球空间信息的三维可视化:基于图形还是基于影像[J]. 测绘学报, 2010, 39(2): 111-114.) |
| [2] | ZHANG Fan, HUANG Xianfeng, LI Deren. Spherical Projection Based Triangulation for One Station Terrestrial Laser Scanning Point Cloud[J]. Acta Geodaetica et Cartographica Sinic. 2009, 38(1): 48-54. (张帆, 黄先锋, 李德仁. 基于球面投影的单站地面激光扫描点云构网方法[J]. 测绘学报, 2009, 38(1): 48-54.) |
| [3] | MANDOW A, MARTíNEZ J L, REINA A, et al. Fast Range-independent Spherical Subsampling of 3D Laser Scanner Points and Data Reduction Performance Evaluation for Scene Registration[J]. Pattern Recognition Letters. 2010, 31(11): 1239-1250. |
| [4] | ZHU Qing. A Survey of Three Dimensional GIS Technologies [J]. Geomatics World, 2004, 2(3): 8-12. (朱庆. 三维地理信息系统技术综述[J]. 理信息世界, 2004, 2(3): 8-12.) |
| [5] | ZHENG Kun, ZHU Liangfeng, WU Xincai, et al. Study on Spatial Indexing Techniques for 3D GIS[J]. Geography and Geo-Information Science, 2006, 22(4): 35-39. (郑坤, 朱良峰, 吴信才, 等. 三维GIS空间索引技术研究[J]. 地理与地理信息科学, 2006, 22(4): 35-39.) |
| [6] | LI Qingquan, YANG Bisheng, SHI Wenzhong, et al. The Theory and Methods of Spatial Data Acquisition, Modeling, Visualization[M]. Wuhan: Wuhan University Press, 2003. (李清泉,杨必胜,史文中,等. 三维空间数据的实时获取、建模与可视化[M]. 武汉:武汉大学出版社, 2003.) |
| [7] | SHI Wenzhong, WU Lixin, Li Qingquan, et al. Models and Algorithms for Three Dimensional Spatial Information System[M]. Beijing: Publishing House of Electronics Industry. 2007. (史文中,吴立新,李清泉, 等. 三维空间信息系统模型与算法[M]. 北京:电子工业出版社,2007.) |
| [8] | GUTTMAN A. R-trees: A Dynamic Index Structure for Spatial Searching[C]//Proceedings of ACM International Conference on Management of Data. Boston: ACM Press,1984: 47-57. |
| [9] | GONG Jun, ZHU Qing, ZHANG Yeting, et al. An Efficient 3D R-tree Extension Method Concerned with Levels of Detail[J]. Acta Geodaetica et Cartographica Sinic, 2011, 40(2): 249-255. (龚俊, 朱庆, 张叶廷, 等. 顾及多细节层次的三维R-索引扩展方法[J]. 测绘学报, 2011, 40(2): 249-255. ) |
| [10] | GONG Jun, ZHU Qing, ZHANG Hanwu, et al. An Adaptive Control Method of LODs for 3D Scene Based on R-tree Index[J]. Acta Geodaetica et Cartographica Sinic, 2011, 40(2): 531-534. (龚俊, 朱庆, 章汉武, 等. 基于R树索引的三维场景细节层次自适应控制方法[J]. 测绘学报, 2011, 40(4): 531-534.) |
| [11] | XIA Yu, ZHU Xinyan, LI Deren. Indexing Technology in Spatial Information Multi-grid[J]. Geospatial Information, 2006, 4(6): 4-7. (夏宇,朱欣焰,李德仁. 空间信息多级网格索引技术研究[J]. 地理空间信息,2006, 4(6): 4-7.) |
| [12] | FU Yuchen, GUO Wei, ZHOU Dongru. Hybrid-tree: An Indexing Structure of Spatial Database[J]. Computer Engineering and Applications, 2003, 39(17): 41-42. (伏玉琛, 郭薇, 周洞汝. 空间索引的混合树结构研究[J]. 计算机工程与应用, 2003, 39(17): 41-42.) |
| [13] | WANG Yanmin. A Study of Vector Spatial Data Organization for Multi-scale GIS[D]. Wuhan: Wuhan University, 2002. (王晏民. 多比例尺GIS数量空间数据组织研究[D]. 武汉:武汉大学, 2002.) |
| [14] | S RUSINKIEWICZ; M LEVOY. QSplat: A Multi-resolution Point Rendering System for Large Meshes[C]//Proceedings of SIGGRAPH 2000, New Orleans: ACM Press, 2000: 343-352. |
| [15] | ZWICKER M. Splatting Fundamentals 3D Editing and Painting[C]//Point-based Computer Graphics. Waltham: Morgan Kaufmann Publishers, 2007. |
| [16] | PFISTER H, ZWICKER M, BAAR J V, et al. Surfels: Surface Elements as Rendering Primitives. Computer Graphics[C]//Proceedings of SIGGRAPH 2000. New Orleans: ACM Press, 2000: 335-342. |
| [17] | REN L, PFISTER H, ZWICKER M. Object Space EWA Surface Splatting: A Hardware Accelerated Approach to High Quality Point Rendering[C]//Proceedings of the Eurographics 2002: 21(3). Oxford: Blackwell Publishers , 2002: 461-470. |
| [18] | MENG Fang. Studies on Interactive Rendering of Huge 3D Point-based Models[D]. Beijing: Peking University, 2005. (孟放. 大型三维点云数据的交互绘制研究[D]. 北京:北京大学, 2005.) |
| [19] | ZHANG Long, DONG Zhao, CHEN Wei, et al. High Quality Real-time Rendering of Large Scale Point Model[J]. Chinese Journal of Computer. 2005, 28(2): 241-249. (张龙, 董朝, 陈为, 等. 大规模点模型的实时高质量绘制[J]. 计算机学报, 2005, 28(2): 241-249.) |
| [20] | WAND M, BERNER A, BOKELOH M, et al. EXGRT—EXTENSIBLE GRAPHICS TOOLKIT[EB/OL]. (2006-06-11)[2011-08-10]. http://www.gris.uni-tuebingen.de/xgrt. |