基于OpenGL三维地形构建技术研究

汪明

(重庆市勘测院，重庆 400020)

1 引言

在解决大规模地形数据三维可视化和实时绘制的问题中，国内外许多学者都进行了广泛、深入的研究，主要集中于地形和纹理数据的分页管理、纹理数据和地形数据的LOD控制、可见域的裁剪等多个方面。此外，随着图形处理器(Graphic Processing Unit，GPU)性能的提高，出现了基于硬件GPU的粗粒度的离散层次细节方法［1］。

本文拟研究基于OpenGL的大规模地形三维地形构建关键技术。由于高分辨率的DTM或DEM的数据集庞大，若直接使用原始数据进行可视化，则需要绘制的三角形或四边形将达到上百万个。因此，为了防止显示滞后，缓解内存和显卡缓存的压力，需要考虑以下几个方面:对复杂的地形做数据简化处理;三维地形数据LOD构建等。

2 三维地形构建技术算法及实现

本文基于多分辨率高程数据文件的层次细节构建地形，主要思想是对大范围的DEM数据先做预处理，即裁剪分层(见3.1)。由于文件按行列组织，因此某一块的索引则为i*m+j(i为行号，j为列号，m为列数)。我们只读取和显示在可视范围内的文件块，这解决了把DEM数据全部读入内存而实际只需要显示地形中很小一部分的矛盾。

2.1 高程数据组织方法

本文中三维地形的构建，将采用多层次文件的形式存储与读取数据，而Esri的GIRD数据却在保存上稍显复杂，为了简化存储与后期处理的方便，将采用自定义文件存储方法，对数据进行多层次组织与保存。

如图1所示，A为低层级数据L，它所存储的当前层级下的数据相对于高层级下的数据B(级别L－1)来说，其存储范围为2×2个B层级数据，但是A的分辨率只有B的1/4。A和B均为所在层级下的一个文件，而DEM数据由很多个这样的文件组成，而且不管是同一层级还是不同层级，每个文件的大小一样(除去边缘上不足块文件栅格数的文件)。本文采用每块文件栅格数为长宽均为256，文件大小为“长×宽×浮点数长度”，浮点数为4个字节，即文件存储大小为256×256×4=262144Byte=256kb。所分层级总数是按照如下公式计算得出:

每一层级(level)文件中相对于源文件重采样点的间隔数XYDiff为:

图1 分级数据示意图

在保存文件时，采用每一层级存放在单独文件夹下面，每一层级文件夹下有多个文件块，分别代表当前级别下256×256个栅格的数据。如图1中A和B中的红色框。每块数据以二进制格式写入文件，后缀名为.bk，文件命名方式为:“当前级别_所在行号_所在列号.bk”，如图1中A级别为L，并在i行j列，因此在级别为L的文件夹下，有文件名为L_i_j.bk。在L同级文件夹中包含一个cfg.cfg的二进制文件，其中包含了DEM的头文件信息，结构体如下:

如果要对所构建的地形进行补缝，则需要在生成数据时在每个块文件中插入相邻的最后一行(列)的数据，即数据需要有重叠部分，如图2中D包含了B有重叠行，而和C有重叠列。

图2 补缝数据分块示意图

2.2 数据裁剪及重采样

数据的裁剪程序使用了GDAL函数库中的重采样方法:

该函数中，nXOff、nYOff表示读取栅格的起始行列;nBufXSize、nBufYSize表示读取栅格的行列数。

通过设定这4个参数，可以实现栅格数据的剪裁和重采样。

2.3 三维地形构建算法

(1)分块分级的地形构建算法

在地形数据分级存储的基础上，确定观察点与观察目标点后，可以根据视点(Object)到观察(Target)点的距离计算出当前读入数据所需要的文件级别以及当前的视域。为保证三维地形显示与DEM大小无关，规定显示的地形最多只能是观察点及其周围8个文件块。这样当读入文件都在同一级别时，读入内存的数据最多为256×9 kb，从而在内存方面保证了实现海量数据的实时快速地形构建的可行性。

在地形显示时，为保证地形显示的效果，读取文件时，不一定只读取一个级别的文件，而更多的是通过四叉树与级别判定公式读取当前显示块、当前显示级别下的文件。如图3中，设L(i，j)为观察点块文件，而周围三个块文件则是经过一定计算会显示的文件块。

图3 地形构建实时分级

首先计算出视点到L(i，j)的空间距离Lij，然后转换成栅格个数RLij，通过公式(3)计算得到当前块的级别(其中blockSize为分块文件的栅格数)。把L(i，j)分为相同的4块，以相同的计算方法计算L(i，j)4块到视点的距离，并计算层级，如果有一个层级比L(i，j)中心点的级别更高，就需要读取更高级别的文件数据，如图3中L(i，j)继续分级后的文件块为图4所示。图5为某一次绘制时，读入的不同级别的文件块。

图4 各级数据文件示意图

图5 视图中的地形数据

为了保存每个文件块(block)文件的数据，我们设计了如下结构体:

该数据结构中，iX、iY分别表示当前块在当前级别(iLevel)下的列、行号;fData表示保存在内存中的文件块中的高层数据;pBlocks表示文件块划分的内存子块，这些子块将通过自绘完成地形的可视化;Width、Height表示内存中pBlocks所包含的内存子块的列、行数。

当需要获得某点高程数据时，只需要知道当前点P在总的DEM栅格中的行列号(x，y)，将行列号减去当前文件块的起始行列号(ox，oy)，可得到P在当前文件夹中的行列号(m，n)∶(m，n)=(x－ox，y－oy)，再由n*Width为索引得到栅格点的值。

(2)块级别评价函数

(3)式定义了数据分块级别判定函数，RLij表示视点到当前文件块或内存子块中心点的栅格个数，中心点的高程值取整个DEM的平均值，且level均为整数。我们定义最精细的层级为1，视点距离目标点的栅格数小于文件分块的栅格数的两倍。(3)式中log2(RLij)求取RLij的对数，－6表示－7+1，即为－log2(BlockSize)+1，文中的BlockSize为文件子块所划分的栅格数，即256。因此，判别公式实际上是(4)式。

(5)式表示一般的视距栅格数求取方法，由公式可以看出相对于当前文件来说的栅格视距RL只与层级level与采样间隔数XYDiff相关，在裁剪数据生成后BlcokSize为常数。由于当前视图读取的文件层级为level，并且当前文件中采样间隔数XYDiff由level确定(2)式，因此可确定RL只与BlockSize相关，且为一常数(6)式。

由于RL为一个常数BlockSize，就可以通过目标点坐标与RL计算出当前级别下视域范围(当前文件级别下的栅格矩形)，就可以在计算出目标点所在的文件块后(图6中文件块5)，根据RL就可以判定，只需要再读入文件块5周围的8个文件块，就可以覆盖整个视域，达到只读取和构建视域范围内地形的目的。

图6 视域中的文件块

(3)四叉树分块过程

3 算法性能分析与讨论

3.1 数据分块分析

测试中所用GRID数据为列12020，行12618，栅格大小 22.5 ×22.5，计算面积为 76782.10725 km2，文件大小578.57 MB，文件裁剪耗时5 min。文件裁剪后分为7个级别的数据，即包含7个文件夹，共3166个文件，总的数据大小为771 Mb，占用空间774 Mb。数据大小增加了33.25%。第一级别的数据大小为578 Mb，占用空间580 Mb，每降低一个级别，数据大小为上一级的四分之一。

3.2 三维地形显示性能分析

在内存消耗方面，相对于把所有的地形数据读入内存中，本文中提出的算法只读入少量数据，仅仅是读取有限的数十个文件块。可以对算法中数据读入大小做一定假设，每次读入的数据最多不超过40个小文件块，即不超过10 Mb。总的来说，读入内存的数据趋于常数，不会因为地形数据变大而变大。

在实时构建的地形中，同一时间所显示的顶点数在256×256=65536到256×256×40=2621400之间。在具体实际中也可以将分开的栅格数降为128，可以减少场景中绘制的顶点数。在实时显示时，场景动画的帧频为十帧以上。场景中绘制顶点数，与读入的地形数据相关，会在一定范围内变动。

3.3 地形缝隙修补

由于在地形构建中有不同级别的数据显示，在不同级别数据的临边上，细节较丰富的块上必定多出部分点，如图7中的两个圆圈。圆圈处不能保证左右两边的面在同一个面上，因此就需要对这些点和面做单独处理。

图7 地形缝隙示意图

对于这种处于边缘，并且相邻子块处于不同级别的子块来说，需要用三角形绘制出采样间隔造成的不重合的点，如图8。但是这种情况只有在相邻级别不超过1的情况下才能满足条件，因此，需要在前期判断和处理，使相邻级别不超过1。

图8 地形补缝原理示意图

在每个64×64的子块自绘时，进行对当前块与其下部块与左边块进行补缝，那就需要知道每个子块的级别以及子块所在的文件块，以便读取高程数据。我们把要显示的文件块用链表存储，并为存储在文件子块结构体中的子块提供相关的信息，比如，子块所在文件在链表中的索引，再通过索引获取结构体中的信息，这样就可以获得所需要的信息。

4 结论

本文提出的基于多级别栅格数据的多层次LOD划分方法，对原有DEM数据文件做层级划分与重采样，生成具有层次和不同分辨率的小文件块。在三维地形实时构建时，计算出视图中要显示的文件与级别，实现了部分地形数据的读取，使内存占用大大减少。在显示时，对内存中数据进一步做LOD分级，让显示的面片更少。

与已有的三维地形构建方法相比，本文中算法实时，由于只对每个块进行层级判别，因此可视化过程中的计算量少，可视化过程中几乎不受计算过程的影响，级别判断时间大概为0.0083 s左右;而且相对于其他算法把全部数据读入内存，本文算法在显示中只需读取可视范围内的DEM数据与影像数据，在可视化过程中占用很少的内存，大概在3Mb～10 Mb之间。在这种算法基础上，基于格网表达的DEM模型数据管理简单易于实现交互操作。

由于在内存使用及计算量方面减少，本文提出的算法总的来说达到了预期目标，提出了一种新的三维地形构建方案，并优化了现有的三维地形在某些方面的不足。

［1］吕希奎，周小平.实战OpenGL三维可视化系统开发与源码精解［M］.北京:电子工业出版社，2009:125－239.

［2］孙敏，薛勇，马蔼乃.基于格网划分的大数据集DEM三维可视化［J］.计算机辅助设计与图形学学报.2002，14(6):566－580.

［3］罗景馨，唐琎.基于改进四叉树分割和结点存储的LOD算法［J］.计算机工程.2009，35(20):202 －204.

［4］徐文学，江涛，姚兵.基于OpenGL的DEM三维仿真方法研究［J］.国土资源信息化.2008(1):7－11.

［5］黄超超，凌永顺，吕相银.ROAM动态地形渲染算法研究［J］.计算机仿真，2005，22(1):216 －219.

［6］周宏伟，杨平，陈琦.实时地形可视化ROAM算法的分块改进［J］.测绘通报，2003，(8):61 －63.

［7］江流长.ROAM地形渲染算法的实现［J］.农业网络信息，2006，(5):42 －44.

［8］涂超.ROAM算法原理及其应用研究［J］.辽宁工程技术大学学报，2003，22(2):176 －179.

［9］解志刚，马秋禾，赵金萍.基于二叉树结构的地形分块实时漫游的研究［J］.海洋测绘，2004，24(5):35－38.

［10］刘修国，张剑波.基于DEM库的地表模型实时简化方法［J］.小型微型计算机系统，2004，25(2):280 －282.

［11］http://www.gdal.org［OL］