退化四叉树格网的全球多分辨率DEM无缝表达

赵学胜，范德芹，王娇娇，王 磊

1.中国矿业大学（北京）地球与测绘工程学院，北京 100083；2.北京师范大学 资源学院，北京 100875

1 引 言

随着空间数据采集技术的飞速发展和全球经济一体化的不断深入，许多应用领域如全球环境变化监测、灾害的预报预警、资源可持续开发、大型工程设计、国防安全乃至战争、“数字地球”等，越来越频繁地使用大范围（甚至全球）高分辨率地形数据进行分析决策。但是，由于受当前的计算机硬件及网络的限制，为了提高显示效率并实现全球DEM数据的无缝绘制和渲染，就需要在保证地形精度的前提下进行DEM格网简化，即构建全球多分辨率DEM表达模型。这样就不可避免地在相邻不同分辨率DEM格网之间产生裂缝，因而消除邻近格网间的裂缝成为全球地形多分辨率连续表达的关键问题之一［1］。

目前传统的裂缝消除方法主要有：垂直边缘法、渐变法、调整高程值法等。垂直边缘法（vertical skirt）［2－4］即在块的边界上建立一个由地表到水平面的垂直外包体，当有块间裂缝存在时，在视觉上裂缝将被“垂直裙”挡住，但并未从本质上消除。此法只适用于不同分辨率分层加载格网的情况，不适用于消除同一层次不同分辨率格网简化时产生的裂缝。渐变法［5－11］采用限制性四叉树（即控制邻近格网的剖分层次差），再通过平滑数据、增减节点或网格线等方式实现裂缝消除。该方法要求相邻地块的剖分层次差不能超过1，并需要时刻检测边界，计算量大［12－13］，若用于全球会产生大量冗余三角形。调整高程值法［14－15］通过调整裂缝处节点的高程值实现无缝拼接，会导致T型节及地形失真，也会带来绘制时的光照不连续现象，对于有些显卡也可能导致一些空洞小点。其他方法还有：自适应网格法［16］、簇依赖（cluster dependencies）法［17］及跳点法［18］等。

上述研究大都针对局部地形进行可视化操作，若应用于全球，将可能大大增加数据量，尤其在南北两极处，将会出现大量的数据冗余，造成不必要的计算资源消耗，降低显示效率［19］。针对上述问题，本文拟采用全球退化四叉树（degenerate quadtree grid）格网作为实现全球多分辨率DEM格网无缝表达的构模框架。全球DQG是一种类似经纬度格网的全球离散格网系统，不同的是涉及极点的格网退化为三角形，而这种退化是规则的和自适应的，既可以直接利用以经纬度格网为参考系的各种新旧数据源，又避免了经纬度格网的非均匀性和极点奇异性问题，且易于构建空间邻近关系和检索机制［20］。本文重点探讨全球退化四叉树的层次分块构建方法，设计了由于格网简化所产生的各类裂缝自适应消除算法。最后，通过属性渲染，实现了全球多分辨率DEM格网的无缝可视化表达。

2 基于退化四叉树的全球分块结构及格网简化

2.1 球面退化四叉树格网剖分原理

球面退化四叉树格网的初始剖分和QTM（quaternary triangular mesh）一样，选取球内接正八面体作为球面格网划分的基础，其顶点占据球面主要点（包括两极），而边的投影则与赤道、主子午线和90°、180°、270°子午线重合，首次剖分将球面划分成8个等正球面三角形（亦称八分体）。在进一步对每一个初始八分体（三角形）进行细分时，首先对初始三角形3个顶点的经纬度进行两两平分，得到3个新点（这3个点位于球面上），将三角形两腰上的两个新点彼此连成一条纬线，再将该纬线的中点与另一新点彼此连成一条经线，这样就形成了一个新的球面三角形和两个四边形，如图1（a）所示；在第二层次，对于三角形部分按第一次剖分的方法进行剖分，对于四边形部分，将四边形4个顶点的经纬度进行两两平分，得到4个新四边形，如图1（b）所示；在第三层次依此类推，如图1（c）所示；如此递归进行，直到满足一定的分辨率要求为止，详细请参考文献［17］。

图1 球面退化四叉树层次剖分Fig.1 Hierarchical partition of spherical degenerate quadtree grid

2.2 DQG球面分块四叉树结构模型

下面在球面DQG剖分的基础上对八分体进行分块，构建球面DQG分块四叉树结构。当剖分层次为4时，如图2所示，将北半球的一个八分体（已划分为三角形）剖分成10个部分。由于八分体结构的对称性，在用四叉树方法生成格网时，只需考虑八分体的一半即可，另一半与其结构相同。以八分体左侧部分为例，将其划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ5个地块。块Ⅰ为三角形，对应为极点附近三角形块；块Ⅱ为四边形块，对应为非四叉树四边形块，其上侧邻近格网为极点三角形块Ⅰ，下侧邻近格网为四叉树块Ⅲ；块Ⅲ的上侧邻近格网为非四叉树四边形块Ⅱ，下侧邻近格网为四叉树块Ⅳ；块Ⅳ上侧邻近格网为四叉树块Ⅲ，下侧邻近格网为四叉树块Ⅴ；四叉树块Ⅴ上侧邻近格网为四叉树块Ⅳ，下侧邻近格网为南半球八分体中与其层次相同的四叉树块。在此八分体中，极点三角形块Ⅰ和Ⅹ、非四叉树四边形块Ⅱ和Ⅸ不需要简化，因其已属最简形式，只有在去除所有顶点高程近似相等的格网或去除海面格网（各顶点高程均为0）时才考虑将其去除；对于四叉树块Ⅲ、Ⅷ、Ⅳ、Ⅴ、Ⅵ 和Ⅶ，则需考虑用四叉树格网细分简化方法进行简化。

球面分块四叉树模型的构建方法如图2所示。在一个八分体中，四叉树块Ⅲ、Ⅳ、Ⅴ、Ⅵ、Ⅶ和Ⅷ 均视为球面八分体的一个地块。将每一地块的中心点作为根节点，从这个根节点出发，检查根节点是否满足某种分割条件，即根据格网4个顶点与其各边中点高差是否在阈值内，若在阈值内则不再对格网进行细分，并将其作为叶节点保存；否则把根节点递归地不断分割成相等的4个节点区域，直到不能再分割为止，若一直不在阈值内就细分到设定的最高层次。对于八分体中的非四叉树块Ⅱ、三角形块Ⅰ，因其已经为最简形式，所以不需对其继续分割。在DQG剖分方法中，全球是由8个八分体构成的，因而当剖分层次为4时，将全球剖分成16个极点三角形块，16个极点附近四边形块，48个四叉树地块。

图2 北半球的一个八分体（剖分层次为4）Fig.2 An octahedron of northern hemisphere（partition level is 4）

该方法可以避免生成多余三角形，最大限度地实现格网简化目标。但这样在节点拼接处会产生格网分辨率差大于1的复杂裂缝。此外，由于是对全球进行分块退化四叉树分割，在各个分块间进行拼接时也将产生不同类型的裂缝，下节重点讨论裂缝的类型及消除方法。

3 多分辨率格网裂缝的自适应消除方法

对于节点分辨率层次不同的情况，消除不同类型的裂缝通常有两种方法：在拼接处增加一条边，或去掉一条边。相对来说，第1种方法更复杂，但是也更全面，适用于拼接处两个节点的分辨率相差任意大的情况。第2种方法则更加简单，但它要求拼接处的两个节点的层次差距最多不超过1。本文综合采用了这两种消除裂缝的方法，并充分利用四叉树索引结构的特点，对格网节点进行搜索，分别对四叉树块内、四叉树块间及四叉树块与非四叉树块间的裂缝进行了消除。

3.1 四叉树块内裂缝的消除方法

根据四叉树块内简化格网与邻近格网相差的层次不同，裂缝的消除分为两种情况：一种是与邻近格网相差一个层次，另一种是相差两个及以上层次。

3.1.1 邻近格网相差一个层次

此时去掉1条边，将相邻的叶节点的两个子三角形进行合并，就可以消除裂缝，具体如图3所示。

四叉树块内简化格网abcd与未简化格网adgf相差一个细分层次，即简化格网abcd的宽度为未简化格网adgf的两倍。根据设定的四叉树节点可知，简化格网abcd的中心点o为叶节点，未简化格网adgf的中心点e为根节点，根据设定的四叉树细分条件，eo之间距离为格网adgf宽度的两倍，此时若o点为根节点则将三角形ade细分，连接ep点；若o为叶节点则不再将三角形ade细分，不连接ep点。类似的，对于简化格网abcd与其下侧、左侧、右侧邻近的未简化格网层次相差1的情况，均采用此法消除裂缝。

3.1.2 邻近格网相差两个及以上层次

当相邻节点细分层次超过1时，消除裂缝的原理是：先根据3.1.1的方法消除相邻节点间相差一个层次的裂缝，然后根据四叉树节点之间的关系，搜索到已简化的格网，再分别按层次搜索与其邻近的上侧、下侧、左侧、右侧相同宽度（即相同细分层次）的格网节点的标识。若相邻格网宽度相同，则不需要消除裂缝；若相邻格网宽度大于简化的格网，亦不需要消除裂缝；若相邻格网宽度小于简化的格网宽度的1／2，则根据节点标识记录裂缝点，并同时存储裂缝三角形的坐标，最后绘制三角形时读出裂缝三角形的坐标文件，实现裂缝的消除绘制与显示。

如图4（a）所示，以简化的四边形abcd为例，只考虑四边形abcd格网上侧相邻的三角形A、B。由于P点高程插值后不一定与直线ad内插中点的高程相等，从而格网简化后可能产生裂缝三角形apd，若要对其进行消除，则需要补充绘制三角形apd。图4（b）所示为消除格网abcd简化引起的裂缝后的效果。采用同样方法进行下侧邻近、左侧邻近、右侧邻近的四叉树格网内部裂缝的消除。

图4 块内简化格网与邻近格网相差两个及以上层次的裂缝消除Fig.4 Level difference between adjacent grids is two or more

3.2 块间裂缝的消除方法

根据裂缝沿纬度方向还是经度方向不同，块间裂缝分为上下块间裂缝和左右块间裂缝，下面首先给出上下裂缝的不同类型及相应的消除方法。如图2所示，上下块间裂缝分为以下几种不同类型。

类型1：块边界上侧为极点三角形块，下侧为非四叉树块，如Ⅰ与Ⅱ。

类型2：块边界上侧为非四叉树块，下侧为四叉树块，如Ⅱ与Ⅲ。

类型3：块边界上侧为四叉树块，下侧也为四叉树块，如Ⅲ与Ⅳ、Ⅳ与Ⅴ。

对于块边界裂缝类型1，只考虑八分体的左侧部分，如图5（a）所示，可见极点三角形块Ⅰ与其底边邻近的非四叉树四边形块Ⅱ间不存在裂缝。

对于块间裂缝类型2（如图2中的Ⅱ部分和Ⅲ部分）：① 若下侧四叉树块Ⅲ已简化，如图5（b）所示，此时上下块Ⅱ、Ⅲ间不存在裂缝；② 若下侧四叉树块Ⅲ未简化（如图5（c）），由于非四叉树格网的层次比其块边界下侧四叉树格网层次大1，将可能在ab中点p处产生裂缝，可通过添加三角形apb消除裂缝，如图5（d）所示。

图5 类型1和类型2Fig.5 Type 1and type 2

对于块边界裂缝类型3，又可根据块边界上下侧简化情况不同细分为3种情况：① 块边界上侧邻近格网未简化，下侧亦未简化或简化后格网宽度小于等于块边界上侧格网，此种情况记为类型3－1；② 块边界上侧格网简化，下侧未简化或简化后宽度小于等于边界上侧格网，此种情况记为类型3－2；③ 块边界下侧格网简化后宽度大于块边界上侧邻近格网未简化或简化后的宽度，此种情况记为类型3－3。

对于块边界类型3－1，由于上下块间在进行DQG剖分时已相差一个剖分层次，而不同块的四叉树结构是相对独立的，因此在块间边界处将出现裂缝。如图6所示，块边界上侧四叉树格网iack未简化，其宽度为ab，块边界下侧格网adeb未简化，其宽度为ap＝ab／2，块边界下侧格网befc已简化，宽度为bc＝ab。可见在添加高程时上下块边界p点处出现裂缝。

图6 类型3－1裂缝消除前后Fig.6 Before and after type 3－1cracks elimination

消除裂缝方法为：

（1）根据块边界上侧四叉树格网的四叉树结构节点标识，搜索出边界上侧邻近的未简化格网的中心点j的坐标及格网宽度ab。

（2）根据j点坐标及格网宽度ab确定块边界下侧邻近格网搜索的最大宽度为ab。

（3）根据块边界下侧四叉树的标识，确定下侧邻近格网的宽度，由此可知未简化格网adeb的宽度小于上侧格网iack的宽度，因此在上下块边界拼接处产生裂缝，并确定裂缝点为p。简化格网bcef的宽度等于上侧格网iack的宽度，因此块边界上侧三角形jbc与块边界下侧三角形bec之间不存在裂缝。

（4）根据p点和已知的a、b点，存储三角形apb的坐标。将四叉树中三角形格网坐标转换为经纬度坐标，再将经纬度坐标转换为三维空间坐标，存储在三维坐标文件中，用于最终三维图形的绘制与显示。

对于块边界裂缝类型3－2、3－3，虽然块边界上下侧简化情况不同，但裂缝消除原理同第1种情况。如图7至8所示，裂缝点分别为p1和p2，消除裂缝需添加的三角形分别为ap1b和p1p2b、ap1d和p1p2d。

图7 类型3－2裂缝消除前后Fig.7 Before and after type 3－2cracks elimination

图8 类型3－3裂缝消除前后Fig.8 Before and after type 3－3cracks elimination

左右块间简化只存在一种情况：块边界左右侧均为四叉树块Ⅴ和Ⅵ（如图2），它们允许细分的最小宽度相同。当块边界两侧简化格网宽度与未简化格网宽度不同时，将会在边界处出现裂缝，消除裂缝的原理及方法与四叉树块内左右侧消除裂缝方法类似。

4 试验结果与分析

试验采用美国地质测量局（USGS）发布的地形数据GTOPO30为数据源，根据双线性插值方法获取格网点高程。应用VC＋＋6.0语言和OpenGL三维工具，设计开发了基于球面退化四叉树的全球多分辨率DEM无缝可视化原型系统。根据GTOPO30数据精度，在本次试验中格网剖分层次最高为12层。试验结果得到裂缝消除前后对照图（如图9a）以及相应的属性渲染图（如图9b）。红色面片及圈出部分表示消除裂缝所需添加的三角形。由于是在球面显示，有些裂缝三角形由于球体旋转的角度不同，在屏幕上捕捉图时看不见，而局部放大图裂缝则清晰可见（如图10）。

图9 全球DQG格网裂缝消除前后对照及相应的属性渲染图Fig.9 Before and after cracks elimination of global DQG and the corresponding render

图10 局部放大图Fig.10 Enlarged figures in local area

格网简化、裂缝消除前后的格网数目及对比分析结果见表1。其中：

（1）裂缝消除代价——定义为以消除裂缝增加的三角形数与格网简化减少的三角形数之比（以%表示）；

（2）简化效率——表示为简化前总格网数（qz）减去格网简化及裂缝消除后总格网数qj除以简化前总格网数qz（以%表示）。

可见，DQG格网简化数目随着剖分层次的递增迅速增加，而消除裂缝增加的三角形数随着剖分层次的递增仅有小幅增加，这使得最终的简化效果十分明显。其中，裂缝消除代价随着剖分层次的递增迅速下降（如图11）；而简化效率随剖分层次的增加而增加（如图12），剖分层次越高，简化效果越明显，12层时的简化效率已接近67%，优于传统的经纬度格网（其简化效率在7层之后趋近于常值45%，如表1和图12）。

5 结 论

本文提出一种基于退化四叉树格网的全球多分辨率DEM无缝可视化表达方法。该方法在不限制相邻节点间剖分层次的前提下，根据地形粗糙度对格网进行充分简化，并对简化过程中产生的各种块内、块间裂缝进行自适应消除。试验结果表明：采用本文方法进行格网简化并消除裂缝后，格网数目大大减少，呈现简化效率随格网剖分层次递增而提高的规律，当格网剖分层次为12层时，简化效率为66.8%，明显优于传统的经纬度格网。初步实现了全球多分辨率DEM的无缝可视化表达，基本满足了全球多分辨DQG格网模型绘制过程中对绘制速度和逼真度的要求。尽管如此，该方法还不是很完善，下一步的工作包括：DQG格网简化和裂缝消除的高效性、交互性设计与空间分析及其大区域地形应用模式等。

表1 格网简化数目对比分析Tab.1 Comparison and analysis of the simplified grid numbers

图11 DQG格网的裂缝消除代价Fig.11 Costs of crack elimination in DQG

图12 DQG和经纬度格网简化效率对比Fig.12 The simplification efficiency of DQG and Long／Lag grid

［1］ XING Wei，SUN Yankui，TANG Zesheng.View－dependent and Multi－resolution Simplification Algorithm of a Terrain Model［J］.Journal of Tsinghua University：Science and Technology，2004，44（1）：29－33.（邢伟，孙延奎，唐泽圣.与视点相关的多分辨率地表模型简化算法［J］清华大学学报：自然科学版，2004，44（1）：29－33.）

［2］ HU Jinxing，MA Zhaoting，WU Huanping，et al.3D Visualization of Massive Terrain Data Based on Grid Partition［J］.Journal of Computer Aided Design ＆ Computer Graphics，2004，16（8）：1164－1168.（胡金星，马照亭，吴焕萍，等.基于格网划分的海量地形数据三维可视化［J］.计算机辅助设计与图形学报，2004，16（8）：1164－1168.）

［3］ MA Zhaoting，PAN Mao，HU Jinxing，et al.A Fast Walkthrough Method for Massive Terrain Based on Data Block Partition［J］.Acta Scicentiarum Naturalum Universitis Pekinesis，2004，40（4）：619－625.（马照亭，潘懋，胡金星，等.一种基于数据分块的海量地形快速漫游方法［J］.北 京 大 学 学 报：自 然 科 学 版，2004，40（4）：619－625.）

［4］ HU Aihua，HE Zongyi，MA Xiaoping.Rendering Technique for Large Scale Terrain Based on LOD［J］.Bulletin of Surveying and Mapping，2009（12）：23－26.（胡爱华，何宗宜，马晓萍.基于LOD的大规模地形实时绘制方法［J］.测绘通报，2009（12）：23－26.）

［5］ LIU Ding，XU Huiping，CHEN Huagen.A LOD Algorithm for Large Scale Ocean Based on OpenGL Index Vertex Array［J］.Journal of Tongji University：Natural Science，2009，37（3）：414－418.（刘丁，许惠平，陈华根.基于OpenGL索引顶点数组的大尺度海面LOD算法［J］.同济大学学报：自然科学版，2009，37（3）：414－418.）

［6］ XU Miaozhong.Research on Real－time Rendering for Large Scale Terrain［J］.Geomatics and Information Science of Wuhan University，2005，30（5）：392－395.（许妙忠.大规模地形实时绘制的算法研究［J］.武汉大学学报：信息科学版，2005，30（5）：392－395.）

［7］ ROTTGER S，HEIDRICH W，SLUSALLEK P，et al.Realtime Generation of Continuous Levels of Detail for Height Fields［C］∥Proceedings of the 6th International Conference in Central European Computer Graphics and Visualization.Plzen：［s.n.］，1998：315－322.

［8］ ZENG Wei，HAN Zhanxiao，ZHU Xuefang.Research of Application of LOD Algorithm to 3DTerrain Simulation［J］.Journal of System Simulation，2009，21（1）：292－294.（曾维，韩占校，朱学芳.LOD算法在3D地表模拟中的应用研究［J］.系统仿真学报，2009，21（1）：292－294.）

［9］ LINDSTROM P，PASCUCCI V.Terrain Simplification Simplified：A General Framework for View－dependent Out－of－core Visualization ［J］.IEEE Transactions on Visualization and Computer Graphics，2002，8（3）：239－254.

［10］ WANG Yuan，LIU Jianyong，JIANG Nan，et al.A Dynamic Triangulation Algorithm for the View－dependent and Real time LOD Model of Terrain［J］.Acta Geodaetica et Cartographica Sinica，2003，32（1）：47－52.（王源，刘建永，江南，等.视点相关实时LOD地形模型动态构网算法［J］.测绘学报，2003，32（1）：47－52.）

［11］ ZHAO Xuesheng，BAI Jianjun，WANG Zhipeng.An Adaptive Visualized Model of the Global Terrain Based on QTM［J］.Acta Geodaetica et Cartographica Sinica，2007，36（3）：316－320.（赵学胜，白建军，王志鹏.基于QTM的全球地形自适应可视化模型［J］.测绘学报，2007，36（3）：316－320.）

［12］ CIGNONI P，GANOVELLI F，GOBBETTI E，et al.Planet－sized Batched Dynamic Adaptive Meshes （PBDAM）［J］.Proceedings of the 14th IEEE Visualization 2003（VIS＇03）.Washington DC：IEEE Computer Society Press，2003：147－154.

［13］ HWA L M，DUCHAINEAU M A，JOY K I.Real－time Optimal Adaptation for Planetary Geometry and Texture：4－8Tile Hierarchies［J］.IEEE Transactions on Visualization and Computer Graphics，2005，11（4）：355－368

［14］ ZHAO Youbing，SHI Jiaoying，ZHOU Ji，et al.，A Fast Algorithm for Large Scale Terrain Walkthrough［J］.Journal of Computer Aided Design ＆ Computer Graphics，2002，14（7）：624－628.（赵友兵，石教英，周骥，等.一种大规模地形的快速漫游算法［J］.计算机辅助设计与图形学学报，2002，14（7）：624－628.）

［15］ RUI Xiaoping，ZHANG Yanmin.An Improved Real－time Continuous LOD Algorithm ［J］.Journal of System Simulation，2004，16（11）：2628－2630.（芮小平，张彦敏.一种实时连续LOD技术的改进算法［J］.系统仿真学报，2004，16（11）：2628－2630.）

［16］ ZHANG Xiaohu，SHAO Yongshe，YE Qin.New LOD Method Based on Adaptive Quad－tree for Terrain Visualization［J］.Journal of Computer Applications，2009，29（9）：2596－2598.（张小虎，邵永社，叶勤.基于自适应四叉树的地形LOD算法［J］.计算机应用，2009，29（9）：2596－2598.）

［17］ YOON S E，SALOMON B，GAYLE R.Quick－VDR：Out－of－core View－dependent Rendering of Gigantic Models［J］.IEEE Transactions on Visualization and Computer Graphics，2005，11（4）：369－382.

［18］ LIU Yang，GONG Adu，LI Jing.A Model for Massive 3DTerrain Simplification Based on Data Block Partition and Quad－tree［J］.Acta Geodaetica et Cartographica Sinica，2010，39（4）：410－415.（刘扬，宫阿都，李京.基于数据分层分块的海量三维地形四叉树简化模型［J］.测绘学报，2010，39（4）：410－415.）

［19］ LI Yachen，JIANG Hongliu，XIONG Hailin，et al.3D Earth Modeling in Visual Simulation［J］.Computer Engineering，2007，33（12）：225－227.（李亚臣，蒋红柳，熊海林，等.视景仿真中三维地球的建模［J］.计算机工程，2007，33（12）：225－227.）

［20］ ZHAO Xuesheng，CUI Majun，LI Ang，et al.An Adjacent Searching Algorithm of Degenerate Quadtree Grid on Spherical Facet［J］.Geomatics and Information Science of Wuhan University，2009，34（4）：479－482.（赵学胜，崔马军，李昂，等.球面退化四叉树格网单元的邻近搜索算法［J］.武汉大学学报：信息科学版，2009，34（4）：479－482.）