-
随着空间数据采集技术的飞速发展和全球经济一体化的不断深入,许多应用领域如全球环境变化监测、灾害预报预警、资源可持续开发、大型工程设计、国防安全乃至战争、“数字地球”等,越来越频繁地使用大范围(甚至全球)高分辨率地形数据进行分析决策。但是,由于受当前的计算机硬件及网络的限制,为了提高显示效率并实现全球数字高程模型(digital elevation model, DEM)数据的无缝绘制和渲染,就需要在保证地形精度的前提下进行DEM格网简化,即构建全球多分辨率DEM表达模型[1]。
目前, 全球多分辨率DEM表达多以四叉树结构为基础,如球面QTM格网[2-3]、球面DQG格网[4-5]、球面菱形格网[6-7]。但是,四叉树结构在实现全球多分辨率DEM表达时,不可避免地在相邻的不同分辨率DEM格网之间产生裂缝。尽管裂缝消除的方法有很多,但每种方法都有各自的局限性,如顶点拉伸法、“限定四叉树”法、针对球面DQG格网的块内块间裂缝消除方法和裙边法等。顶点拉伸法[8]是对裂缝处的顶点高程进行调整,可能会导致T型节及地形失真;“限定四叉树”法[9]要求相邻格网剖分层差不超过1,导致数据存在冗余,并且需要反复检测相邻格网的剖分层差,计算量大;块内块间裂缝消除方法[1]需要根据格网位置和相邻格网剖分层差的不同,采用不同的裂缝消除方法,方法复杂;裙边法[4]的思想是在块的边界上建立一个由地表到水平面的垂直外包体,当有块间裂缝存在时,在视觉上裂缝将被“垂直裙”遮挡,但并未从本质上消除。
分析发现,裂缝产生的根本原因是:四叉树将地形抽象为格网面的集合,四叉树剖分过程操作的唯一拓扑实体是面,未涉及格网边的细分;另外,现有的裂缝消除方法都是裂缝产生后的“修补”操作,属于“滞后”处理。为此,本文引入多分辨率组合映射理论[10],将格网抽象为dart的集合,以多分辨率半边结构存储和管理多分辨率格网,构建一种无缝的全球多分辨率DEM建模方法,有望从格网剖分结构上消除全球多分辨率DEM建模中的“裂缝”问题。
-
二维组合映射(2-maps)模型是一个三元组M =(D, α0, α1)。其中,D表示集合,集合中的元素被称为dart,α0是对合运算,α1是排列运算;2-maps是半边数据结构的理论模型,即2-maps特指半边数据结构时,dart具体化为半边,D表示半边集合,α0和α1分别具体化为反向半边指针和后继半边指针[11-13];多分辨率二维组合映射(多分辨率2-maps)是2-maps在多分辨率表达方面的扩展[14-16]。基于以上理论,可以实现多种全球离散格网的统一表达。如图 1所示,第i层次的格网均表达为Mi =(Di, α0i, α1i),其中,α0i和α1i分别表示第i层半边的反向半边指针和后继半边指针,Di表示第i层的格网半边集合,满足表达式[17]:
图 1 基于多分辨率组合映射的3种全球离散格网统一表达
Figure 1. Unified Expression of Three Global Discrete Grids Based on Multi-resolution Combinatorial Maps
$$ {D^0} \subset {D^1} \subset \cdots {D^i} \subset \cdots, {D^i} = \bigcup\limits_{j = 0}^i {{L^j}}, D = \bigcup\limits_{i \ge 0} {{L^i}} $$ -
基于多分辨率组合映射理论,格网面的表达式为orbit < α0∘α1 > ,即多分辨率半边对格网面的管理通过复合运算α1(α0(e))实现。
如图 2所示,假设面A内部的半边的排列顺序为逆时针,提取面A的所有顶点,可以从面A的任意一个半边开始,比如e1,先对其求对合运算α0,得到反向半边e2,再对e2求排列运算α1,得到面A的下一条半边e3,反复执行α1(α0(e)),直到起始半边e1,最终得到e1e3e5e7e9e11,半边的起点集合就是面的顶点,该顶点用于地形格网的渲染,实现了半边对格网面的管理。面A此时已不再是四边形,而是由6个半边首尾相连组成的六边形,多分辨率半边结构能够自适应地管理多边形。多分辨率半边结构通过半边间接管理格网面,能够实现格网面与边的同步剖分,因而可以从前期的格网剖分结构上,避免不同层次相邻格网边界处出现裂缝。
-
多分辨率半边编码的目的是建立多分辨率半边结构。球面菱形格网半边编码由5个元素组成,分别是四分体码、剖分层次、行号、列号、半边方向码,即(Q+k+i+j+D)。其中,四分体码与常规方法相同[1],值为0、1、2、3,在编码中设置为1位;剖分层次是半边起点的引入层次,设置为2位;在如图 3所示的坐标系中,行号和列号是半边起点所在的行列号,设置为4位;半边方向码如图 3黑点处所示,在编码中设置为1位。以图 3黑点处方向码为2的半边为例,假如其位于第0个四分体内,其编码为002000300022。
-
坐标转换是指将多分辨率半边编码转换为空间直角坐标,其目的是用于可视化、度量关系分析和高程点的获取等。坐标转换分为两个过程:一是半边编码转经纬度坐标;二是经纬度坐标转化为空间直角坐标。第二个过程与常规方法相同,不再赘述;关于第一个过程,如图 3所示,处于同一水平高度的点的行列号(i, j)满足图右侧表达式,通过这种规律,稍加推导,便可以得到半边编码转经纬度的公式:
$$ {N_\theta } = Q \times 90^\circ + \frac{{90^\circ }}{{{2^k}-\left( {i-j} \right)}} \times j $$ (1) $$ {N_\varphi } = 90^\circ-\frac{{90^\circ }}{{{2^k}}} \times \left( {i-j} \right) $$ (2) $$ {S_\theta } = {S_\varphi } = Q \times 90^\circ + \frac{{90^\circ }}{{{2^k} + \left( {i-j} \right)}} \times i $$ (3) 式中,Nθ、Nφ分别是北半球的经度和90°-纬度;Sθ、Sφ分别是南半球的经度和90°+纬度。
-
基于多分辨率半边的全球多分辨率DEM无缝表达算法分为3个阶段。
1) 对每个四分体进行单层次多分辨率菱形格网细分,分3个小步骤:
(1) 对初始格网建立半边结构,并加入节点哈希表中。
(2) 初始格网建立后,遍历节点哈希表中上一层次剖分的节点(原节点)。如果原节点位于边界处,则直接对3个半边建立半边结构;如果节点为非边界点,那么求原节点相邻的节点(邻节点)。然后,判断邻节点是否存在,如果不存在,则对3个半边建立半边结构;如果存在,继续判断邻节点是否剖分。如果没剖分,则建立3个半边的半边结构;如剖分但未建立半边结构,则建立4个半边的半边结构。
(3) 递归执行第(1)、(2)步,直到遍历完上一层次所有节点。
2) 对每个四分体进行多层次多分辨率菱形格网细分,分2个小步骤:
(1) 计算本层节点的细分准则,计算4个顶点和中心点的最大高差。
(2) 如果超出阈值则继续剖分,将下一层节点添加到节点哈希表,标注该层节点已剖分,求得新加入节点哈希表中的节点数n,并将判断层次加1。递归执行(1)过程,直到没有新的节点加入,即n等于0时,建模过程结束。
3) 全球菱形格网面的提取及可视化,分2个小步骤:
(1) 遍历半边哈希表,如果半边没有被绘制,执行复合运算α1(α0(he)),实现单个格网面的提取,并将查找到的半边标注为已绘制,直到半边哈希表的结尾。
(2) 应用DirectX工具,分别绘制线带(LINESTRIP)和三角形扇(TRIANGLEFAN),得到未三角化和三角化后的球面菱形格网。
-
本文实验数据为美国地形测量局于1996年发布的GTOPO30数据,数据分辨率为30 s(约为1 km),球面菱形格网点高程值的获取采用双线性插值。初期实验发现,虽然四分体内部不产生裂缝,但是四分体间是独立的,如果采用统一的细分标准,在四分体的边界处和极点处仍会出现少量的裂缝现象,所以边界处和极点处结点的细分准则需要特殊处理,如图 4所示。其中,图 4(a)表示边界结点,当A、B两结点至少有一个满足细分条件时,两个都要细分;图 4(b)表示极点处结点,C、D、E、F这4个结点至少有一个满足细分条件时,4个都要细分。因为多分辨率半边结构可以保证四分体内部无裂缝,所以边界和极点处结点可以采用这种简单的方式解决裂缝问题。
本文实验采用C++语言,在VS2010集成开发环境下完成,实验结果的三维显示由Direct3D实现,如图 5所示。其中,图 5(a)、5(b)是在未消除裂缝的情况下,分别以线框和面的形式显示多分辨率球面菱形格网,通过局部放大,对应图 5(e)和图 5(f),可以看出,在不同层次相邻格网边界处存在大量的裂缝问题,因可视化角度问题,有些裂缝未能看清楚,因此只对能够看清的裂缝以椭圆标出。图 5(c)和图 5(d)是基于多分辨率半边结构,根据本文算法生成多分辨率全球菱形格网,并分别以线框和面的形式显示,通过局部放大,对应图 5(g)和图 5(h),可以看出,不同层次相邻格网边界处不存在裂缝问题。
表 1详细列出了裂缝消除前后所需渲染三角形的个数,可以看出,随着剖分层次的增加,简化效率逐渐增大。以第10层为例,地形未简化时需要渲染三角形个数为16 777 216个,按照阈值500 m简化后,在不消除裂缝的情况下,需要渲染三角形个数为580 852个。采用本文提出的算法消除裂缝后,需要多渲染三角形27 614个,比例为4.75%。本文算法对全球地形的简化效率为96.37%,渲染三角形个数大幅降低。通过定义裂缝消除代价指数,发现裂缝消除的代价随剖分层次的增加逐渐降低,当剖分层次为10时,该值仅为0.17,裂缝消除代价很小。
表 1 裂缝消除前后对比分析
Table 1. Comparison and Analysis of the Spherical Diamond Grids Before and After Crack Elimination
层次
(A)需渲染三角形个数 裂缝消除前后对比 简化效率(%)
(G=1-D/B)裂缝消除代价(%)
(F=E/(B-D))未简化
(B)简化但未消除裂缝(C) 简化并消除裂缝(D) 多渲染三角形个数(E=D-C) 多渲染比例(%)
(F=E/C)4 4 096 2 140 2 396 256 11.96 41.50 15.06 5 16 384 4 900 5 470 570 11.63 66.61 5.22 6 65 536 12 208 13 200 992 8.13 79.86 1.90 7 262 144 31 864 33 868 2 004 6.29 87.08 0.88 8 1 048 576 83 560 87 928 4 368 5.23 91.61 0.45 9 4 194 304 221 116 231 402 10 286 4.65 94.48 0.26 10 16 777 216 580 852 608 466 27 614 4.75 96.37 0.17 -
分别以限定四叉树和多分辨率半边作为全球多分辨率DEM无缝表达的数据结构,实现多分辨率全球菱形格网无缝表达,如图 6所示。图 6(b)和图 6(d)分别是图 6(a)和图 6(c)的局部放大图。从椭圆形区域能清晰看到,限制性四叉树算法需要限制相邻格网的剖分层差不大于1,导致一些不满足剖分条件的格网被剖分;而基于多分辨率半边结构,相邻格网的剖分层差没有限制,可实现完全依据地形起伏的地形简化,无数据冗余。与此同时,得到基于以上两种结构简化后的三角形个数,如表 2所示。由于限定四叉树算法[18]需要保证相邻格网的剖分层差不大于1,需要反复检测相邻格网的剖分层次,计算量大,当阈值为10 m时,本文实验未能得到第九层和第十层的实验数据;多分辨率半边结构比限定四叉树结构的简化程度更高,当剖分层次为10,阈值为500 m时,限定四叉树比多分辨率半边多渲染三角形91 922个,比例为15.11%。
图 6 基于限制性四叉树和多分辨率半边的全球多尺度DEM无缝可视化表达
Figure 6. Seamless Visualization of Global Multi-scale DEM Based on Restrict Quad-Tree and Multi-resolution Half-Edges
表 2 格网简化对比分析
Table 2. Comparison and Analysis of the Simplified Grid Numbers
层次 简化前
(A)限定四叉树简化后三角形个数 多分辨半边简化后三角形个数 简化程度对比 ε=500
(B)ε=10
(C)ε=500
(D)ε=10
(E)ε=500
(F=B-D)ε=500
(G=B/D-1)(%)4 4 096 2 396 2 924 2 396 2 924 0 0.00 5 16 384 5 927 8 514 5 470 8 198 457 8.35 6 65 536 15 143 27 629 13 200 25 930 1 943 14.72 7 262 144 39 506 94 473 33 868 88 912 5 638 16.65 8 1 048 576 101 530 335 991 87 928 321 672 13 602 15.47 9 4 194 304 266 279 - 231 402 1 201 072 34 877 15.07 10 16 777 216 700 388 - 608 466 4 623 802 91 922 15.11 -
本文基于多分辨率半边结构提出了一种全球多分辨率DEM无缝表达方法。该方法从格网剖分结构上解决了多分辨率DEM表达中的裂缝问题,且不需要限制相邻格网的剖分层差,可完全根据地形起伏对格网进行彻底简化。实验结果表明, 以多分辨率半边结构为基础,通过边界和极点处结点细分准则的特殊处理,四分体内部、边界和极点处均不会产生裂缝。当剖分层次为10时,裂缝消除代价很小,仅为0.17;当剖分层次为10,阈值为500 m时,“限定四叉树”算法比本文算法多渲染三角形91 922个,比例为15.11%。本文实现了全球多分辨率DEM的无缝可视化表达,基本满足了全球多分辨DEM建模对绘制速度和逼真度的要求。
本文算法在高效建立多分辨率半边结构中存在不足,下一步将继续研究多分辨率半边的存储方式,并结合并行技术提高建模算法效率。
Seamless Expression of Global Multi-resolution DEM Based on Multi-resolution Half-Edges Structure
-
摘要: 应用传统四叉树结构进行全球多分辨率数字高程模型(digital elevation model,DEM)表达时,不同层次相邻格网间会产生裂缝问题,尽管目前有许多裂缝消除的方法,但大都是属于"后期处理",且存在诸多限制。为此,引入多分辨率组合映射理论,提出了一种基于多分辨率半边结构的全球多分辨率DEM无缝表达方法。首先给出了基于多分辨率组合映射的多种全球离散格网统一表达方法和裂缝消除原理,并以球面菱形格网为例,提出多分辨率半边编码方法和坐标转换方法;然后,设计并实现了一种全球多分辨率DEM无缝表达算法;最后,应用C++语言和DirectX工具,开发了相应的可视化实验系统。实验结果表明,利用多分辨率半边结构,通过边界结点细分方法的特殊处理,球面菱形格网的四分体内部、边界和极点处均不会产生裂缝。与传统方法相比,该方法从前期的格网剖分结构上解决了多分辨率DEM格网的裂缝问题。Abstract: When performing the multi-resolution DEM expression based on Quad-tree structure, topological crack generates at the meeting of areas of different resolution levels. Holes appear in the mesh between neighboring faces of different refinement depths. Although there are many methods of eliminating the crack, but mostly belong to the "post processing", and there are many restrictions. A new theory of multi-resolution combination maps is introduced, and a seamless expression of global multi-resolution DEM is discussed based on multi-resolution half-edges. The main contents are as follows. Firstly, the unified expression and principle of crack elimination of Global discrete grids are studied based on multi-resolution combination maps, and some methods are put forward based on spherical diamond grids, including encoding and coordinate transformation. Then, seamless expression algorithm is designed and achieved. Finally, by using C++ language and DirectX tools, an experimental system is developed, and the results show that seamless expression of the global multi-resolution DEMs can be carried out with a method of making a tiny adjustment of node subdivision. Compared with the traditional methods, the method in the paper can solve the crack problem from the underlying data structure.
-
表 1 裂缝消除前后对比分析
Table 1. Comparison and Analysis of the Spherical Diamond Grids Before and After Crack Elimination
层次
(A)需渲染三角形个数 裂缝消除前后对比 简化效率(%)
(G=1-D/B)裂缝消除代价(%)
(F=E/(B-D))未简化
(B)简化但未消除裂缝(C) 简化并消除裂缝(D) 多渲染三角形个数(E=D-C) 多渲染比例(%)
(F=E/C)4 4 096 2 140 2 396 256 11.96 41.50 15.06 5 16 384 4 900 5 470 570 11.63 66.61 5.22 6 65 536 12 208 13 200 992 8.13 79.86 1.90 7 262 144 31 864 33 868 2 004 6.29 87.08 0.88 8 1 048 576 83 560 87 928 4 368 5.23 91.61 0.45 9 4 194 304 221 116 231 402 10 286 4.65 94.48 0.26 10 16 777 216 580 852 608 466 27 614 4.75 96.37 0.17 表 2 格网简化对比分析
Table 2. Comparison and Analysis of the Simplified Grid Numbers
层次 简化前
(A)限定四叉树简化后三角形个数 多分辨半边简化后三角形个数 简化程度对比 ε=500
(B)ε=10
(C)ε=500
(D)ε=10
(E)ε=500
(F=B-D)ε=500
(G=B/D-1)(%)4 4 096 2 396 2 924 2 396 2 924 0 0.00 5 16 384 5 927 8 514 5 470 8 198 457 8.35 6 65 536 15 143 27 629 13 200 25 930 1 943 14.72 7 262 144 39 506 94 473 33 868 88 912 5 638 16.65 8 1 048 576 101 530 335 991 87 928 321 672 13 602 15.47 9 4 194 304 266 279 - 231 402 1 201 072 34 877 15.07 10 16 777 216 700 388 - 608 466 4 623 802 91 922 15.11 -
[1] 赵学胜, 范德芹, 王娇娇, 等.退化四叉树格网的全球多分辨率DEM无缝表达[J].测绘学报, 2012, 41(6):918-925 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=chxb201206023&dbname=CJFD&dbcode=CJFQ Zhao Xuesheng, Fan Deqin, Wang Jiaojiao, et al. Seamless Expression of the Global Multi-resolution DEMs Based on Degenerate Quadtree Grids[J].Acta Geodaetica et Cartographica Sinica, 2012, 41(6):918-925 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=chxb201206023&dbname=CJFD&dbcode=CJFQ [2] 赵学胜, 白建军, 王志鹏.基于QTM的全球地形自适应可视化模型[J].测绘学报, 2007, 36(3):316-320 http://www.wenkuxiazai.com/doc/3882640949649b6648d747e4.html Zhao Xuesheng, Bai Jianjun, Wang Zhipeng. An Adaptive Visualized Model of the Global Terrain Based on QTM[J].Acta Geodaetica et Cartogr-aphica Sinica, 2007, 36(3):316-320 http://www.wenkuxiazai.com/doc/3882640949649b6648d747e4.html [3] 白建军, 赵学胜, 陈军.基于线性四叉树的全球离散格网索引[J].武汉大学学报·信息科学版, 2005, 30(9):805-808 http://ch.whu.edu.cn/CN/abstract/abstract2276.shtml Bai Jianjun, Zhao Xuesheng, Chen Jun.Indexing of Discrete Global Grids Using Linear Quadtree[J]. Geomatics and Information Science of Wuhan University, 2005, 30(9):805-808 http://ch.whu.edu.cn/CN/abstract/abstract2276.shtml [4] 王磊, 赵学胜, 殷楠.基于DQG的全球地形实时可视化系统研究[J].系统仿真学报, 2014, 26(9):1975-1979 http://d.old.wanfangdata.com.cn/Periodical/xtfzxb201409018 Wang Lei, Zhao Xuesheng, Yin Nan.Research on DQG-Based Global Terrain Real-Time Visualization System[J]. Journal of System Simulation, 2014, 26(9):1975-1979 http://d.old.wanfangdata.com.cn/Periodical/xtfzxb201409018 [5] 崔马军, 赵学胜, 孙乾坤, 等.基于球面DQG的数字高程建模算法与可视化表达[J].地理与地理信息科学, 2008, 24(5):33-35 http://d.old.wanfangdata.com.cn/Periodical/dlxygtyj200805008 Cui Majun, Zhao Xuesheng, Sun Qiankun, et al. Algorithm and Visualization of Digital Elevation Modeling Based on Spherical DQG[J]. Geography and Geo-Information Science, 2008, 24(5):33-35 http://d.old.wanfangdata.com.cn/Periodical/dlxygtyj200805008 [6] 张玉梅, 陈维华, 聂洪山, 等.球面菱形网格递归剖分方法研究[J].地理与地理信息科学, 2010, 26(6):34-37 https://www.cnki.com.cn/lunwen-1013047708.html Zhang Yumei, Chen Weihua, Nei Hongshan, et al. Study on Sphere Rhombus Grid Recursive Subdivision[J]. Geography and Geo-Information Science, 2010, 26(6):34-37 https://www.cnki.com.cn/lunwen-1013047708.html [7] 赵学胜, 白建军.基于菱形块的全球离散格网层次建模[J].中国矿业大学学报, 2007, 36(3):397-401 http://www.doc88.com/p-6397751632772.html Zhao Xuesheng, Bai Jianjun.Hierarchical Model of Global Discrete Grids Based on Diamonds[J]. Journal of China University of Mining & Technology, 2007, 36(3):397-401 http://www.doc88.com/p-6397751632772.html [8] 孙文彬, 赵学胜.基于QTM格网的空间数据无缝层次建模[J].中国矿业大学学报, 2008, 37(5):675-679 http://d.old.wanfangdata.com.cn/Periodical/zgkydxxb200805018 Sun Wenbin, Zhao Xuesheng.A Hierachical Seamless Model of Spatial Data Based on QTM[J]. Journal of China University of Mining & Technology, 2008, 37(5):675-679 http://d.old.wanfangdata.com.cn/Periodical/zgkydxxb200805018 [9] 王洪斌. 基于Morse复形的地表形态建模及应用研究[D]. 北京: 中国矿业大学(北京), 2014 http://cdmd.cnki.com.cn/Article/CDMD-11413-1014371349.htm Wang Hongbin. Research on Modeling Terrain Morphology and its Application based on Morse Complex[D]. Beijing: China University of Mining & Technology (Beijing), 2014 http://cdmd.cnki.com.cn/Article/CDMD-11413-1014371349.htm [10] Kraemer P, Untereiner L, Jund T, et al. CGoGN: N-dimensional Meshes with Combinatorial Maps[C]. The 22nd International Meshing Roundtable, New York, 201 [11] Damiand G, Dupas A. Combinatorial Maps for 2D and 3D Image Segmentation[M]. Netherlands:Springer, 2012 [12] Dupas A, Damiand G. First Results for 3D Image Segmentation with Topological Map[C]. The 14th IAPR International Conference, Lyon, France, 2008 [13] Damiand G. Topological Model for 3D Image Representation:Definition and Incremental Extraction Algorithm[J]. Computer Vision and Image Understanding, 2008, 109(3):260-289 doi: 10.1016/j.cviu.2007.09.007 [14] Kraemer P, Cazier D, Bechmann D. A General and Efficient Representation for Multiresolution Meshes: Application to Quad/Triangle Subdivision[C]. Canadian Conference on Computational Geometry, Ottawa, Ontario, 2007 http://hal.cirad.fr/INSA-STRASBOURG/hal-01208563 [15] Mousa M H, Hussein M K. Multi-resolution Surface Representation Using Combinatorial Maps[J]. Computers, 2012, 2(6):103-110 [16] Untereiner L, Cazier D, Bechmann D. N-dimen-sional Multiresolution Representation of Subdivision Meshes with Arbitrary Topology[J]. Graphical Models, 2013, 75(5):231-246 doi: 10.1016/j.gmod.2013.03.003 [17] Kraemer P, Cazier D, Bechmann D. Extension ofHalf-Edges for the Representation of Multiresolution Subdivision Surfaces[J]. The Visual Computer, 2009, 25(2):149-163 doi: 10.1007/s00371-008-0211-6 [18] Hill D. An Efficient, Hardware-Accelerated, Level-of Detail Rendering Technique for Large Terrains[D]. Toronto: University of Toronto, 2002 http://www.researchgate.net/publication/249886465_AN_EFFICIENT_HARDWARE-ACCELERATED_LEVEL_OF-DETAIL_RENDERING_TECHNIQUE_FOR_LARGE_TERRAINS -