一种基于SGOG瓦块的大区域真三维地理场景自适应构建方法与流程

本发明属于地理信息科学技术领域，具体涉及一种基于SGOG瓦块的大区域真三维地理场景自适应构建方法。

背景技术：

1.研究背景

地理信息科学中，三维地理场景的构建方法一般有两种：一种是面元模型加纹理的方法，往往可以达到较好的地形和建筑可视化效果，在3DGIS领域应用广泛；另一种是基于体素模型的建模，常用来表达地质体垂向的多个界面，并方便水平与垂直剖面的制作，在三维地学模拟信息系统中应用广泛。面模型不能很好表达洞穴、拱桥和悬崖等复杂的地理实体，难以实现地形体数据表达和实时编辑，基本上不可能实现地理空间过程的模拟；而基于体素的模型可以实现虚拟地理场景的快速构建与多分辨率实时拆分，实现体素模型、表面模型的一体化表达、多维地理对象及地理过程的一体化表达，地理空间数据可定制及自适应表达，面向复杂地学分析，能够较好地表示空间粒子之间的相互作用关系，适用于表达地理空间过程的动力学机制，进而实现对多尺度地理现象与过程的模拟、分析与预测。但上述传统方法，无论是面元模型还是体素模型，都是面向局部的、小范围的欧氏投影模型。数字地球背景下，传统的投影模型在面向大尺度和全球性问题的应用时，存在着严重的局限。

在大区域或全球尺度，地球曲率不容忽略。因而，投影模型不再适用，必须基于顾及地球曲率的全球离散网格模型(Discrete Global Grid Systems，DGGs)构建三维场景,利用退化四叉树网格实现了全球多尺度地形无缝可视化。其实，以Google为代表的50余种第一代数字地球平台，都是利用球面剖分网格(瓦片或面片)构建全球DEM框架。但以上所有基于DGGs的三维模型归根结底仍属于表面模型的范畴，不涉及地表上下的空间。近些年来，在地理信息科学领域，球体剖分网格，即地球系统空间网格(Earth System Spatial Grid，ESSG)模型被提出,利用球体退化八叉树网格(Spheriod Degenerated Octree Grid,SDOG)建立了青藏高原岩石圈的三维建模。基于ESSG瓦块的地理空间真三维表达方兴未艾，亟待发展。需要说明的是，地球系统空间网格与传统的地球物理和地球系统网格有很大区别。而地球系统网络是一种科学计算的工具，基本不涉及或不强调地球场域的统一定位框架及一体化空间数据的组织、管理与应用。

2.SGOG剖分方案

球体大圆弧QTM八叉树网格(Sphere Geodesic Octree Grid，SGOG)是球体三维网格剖分理论。QTM(Quaternary Triangular Mesh)指四元三角网。简单的理解就是把一个三角形三边的中点用直线连接，就可以把一个三角形划分为四个子三角形，同理，四个子三角形又可以划分为16个孙三角形……，以此类推。设定剖分层次，就可以构成一张三角形格网。显然，四元三角网是一种面剖分方式，可以是平面或曲面。球体QTM八叉树网格剖分方法可分为两步：把地球看做一个球体，首先进行球面离散剖分，然后再进行球体径向剖分。具体为：第一步，球面剖分。首先，以0°～180°首子午圈和与之垂直的东西经90°子午圈和赤道把球面分为等积的8个球面三角形。接着，对每个三角形分别取三边弧的中点，并用大圆弧连接，形成四个二级的球面三角形；然后依据此方法依次对此后各级的球面三角形进行递归剖分，直到满足应用为止(从八分之一球体起算，递归的层次数为n，下同)。第二步，径向剖分；设球的半径为R，从网格的每个节点向球心引直线，然后以R/2ⁿ将半径等分，在各分层处以球面进行剖分(以球心为原点，以k*R/2ⁿ为半径画球，k＝1,2,…2ⁿ-1)。这样，就可以把球分割为整齐一致的球面三棱台(上下底面为球面，侧面为平面)，其中球心处为球面三棱锥。以棱台/锥的几何中心为网格参考点。把球理解为包括大气层在内的地球，可实现对整个重力场的剖分，其三维可视化效果如图1、2所示。

球体SGOG剖分还包括径向不等长八叉树剖分，由于本发明不涉及不等长剖分，所以这里不再介绍。

3.SGOG编码模型

已有研究表明：在重力场中，从地心到磁层，16倍的地球半径就可以满足包含地球系统所有圈层、基本覆盖人类活动空间的要求。考虑到网格定位及计算的方便，SGOG网格采用以下编码模型：圈层码(十六进制)+八分体标识码(八进制)+球面位置编码(四进制)+径向深度码(二进制)。这里的圈层，不是严格意义上的地球圈层，而是以整数倍地球半径表示的从地心和磁层的距离。八分体标识码表示网格所在的八分体的位置，从0°开始，沿经度增加的方向，按每90°为一个卦限，北半球分为0～3、南半球分为4～7共八个卦限。球面QTM四叉树网格已有很多种成熟的编码方案，如固定方向编码、ZOT编码和LS编码等。本发明采用固定方向编码，径向深度采用二叉树编码(靠近球心的位置取0，如图4所示)。

球面QTM四叉树固定方向编码的原理：

如图3所示，该法将八面体按等边三角形投影(Equal-Triangles Projection，ETP投影)，首先把经纬度通过ETP投影转换为x，y值；然后比较位置点(x，y)到四个三角形中心点的距离；选择距离最小的三角形，记录其地址码；依此类推，直到划分到一定层次或位置点到三角形中心的距离小于一定的值。

径向二叉树编码原理：

如前所述，球体八叉树剖分是球面四叉树和球径二叉树剖分的组合，用球面四叉树编码加上径向二叉树编码即可确定任一八叉树网格的空间位置(一个完整的八叉树网格编码还包括首位0,1,2，……7的八分体识别码)。二叉树的基本编码方法如下：把每次剖分的格按从球心到球面的次序转化为与层次相对应的二进制位编码即可。如图4所示：第一次剖分的编码为球心处0、1，第二次剖分为00、01、10、11，依此类推。

4.SGOG解码方案——编码到经纬度和空间直角坐标的相互换算

1)编码到经纬度坐标的相互换算

所谓离散网格编码到经纬度的相互转换是指：已知网格的编码，求网格体元中心点的经纬度；或者反过来，给定某点的经纬度求该点所在某层次(或满足某种精度)网格的体元编码。本发明以一倍的地球半径为例，并将网格的定位分为球面四叉树经纬度求解和径向二叉树深度(离地心的距离)求解两个步骤。

依据相关研究实现方向四叉树编码到经纬度相互转换。该法的核心是将八面体按等边三角形投影，把经纬度转换为笛卡尔直角坐标，并以其为桥梁，实现从地址码到经纬度的相互变换。其主要数学模型如式(1)和(2)所示。

式中，λ、φ为经纬度；x、y为ETP投影坐标。

径向二叉树编码的实质就是以某层次的网格棱长去度量地球的半径。径向二叉树的编解码原理：设某一点到球心的距离为d，剖分的层次为n，用d除以地球的平均半径得商e，把e转化为以2ⁿ为分母的分数，取其分子的整数部分加上1(如果为整数则不加1)，然后再将其化为n位二进制码，即得到其二叉树码；反之，已知某二叉树码(对于一个完整的八叉树码，需根据其规则截取二叉树码)，设其位数为n，把其转化为十进制数，再乘以地球的平均半径与2ⁿ的商就得到该编码所对应的点到球心的距离。由此可知，二叉树的编码没有误差，而解码的精度在一个网格范围内。变长八叉树的编解码方法与八叉树原理上是一样的，其精度也相仿，但剖分的比例不同。将变长比例代入相应的解码公式即可。

2)编码到空间直角坐标的相互换算

这里网格编码到空间直角坐标的转换是指由网格的编码求网格单元各顶点的空间直角坐标，进而可以求得SGOG瓦块内外三角形的重心坐标，将二者取平均值，可以得到SGOG瓦块的几何中心的坐标；而空间直角坐标到网格编码的转换是指已知某点的空间直角坐标和剖分层次，求其所在网格单元的编码。由于SGOG采用了大圆弧中分的剖分规则，所以整个网格体系与空间直角坐标对应十分整齐。下层网格瓦块新增顶点与上层瓦块的顶点之间存在着简单的中分关系。通过求出弦的中点，将其投影到相应的剖分球面上即可。

本发明在算法设计时，对上面的编码模型做了一些变动：将整个码用标示符A和B分为三个码段。A之前的为圈层码，由0～n个1构成。若没有1，则表示是固体地球本身，即1倍的地球半径；每多一个1，球的半径就乘以2，以此类推。A和B之间的二叉树码表示网格体元瓦块的径向位置。B之后的第一位八进制码为瓦块所在的八分体标识码，其余的四进制码表示瓦块在球面的横向位置。根据以上编码和剖分规则，给定初始条件，设计八叉树编码与空间直角坐标转换的算法思路如图5所示。图5为编码到空间直角坐标，图6为空间直角坐标到编码。

5.SGOG网格可视化绘制技术

SGOG网格体系包括球面三棱台和球面三棱锥两种瓦块。瓦块是离散网格的基本单位，也是地理空间建模的基础。因而，对瓦块结构体的设计与实现是首要的关键技术。本发明中的瓦块绘制程序基于C++平台的应用程序接口OSG开发，其结构体如下：

由结构体绘制线框模型共需3步：

(1)外层3个顶点两两之间绘制大圆弧；

(2)内层3个顶点两两之间绘制大圆弧；

(3)内外层对应顶点之间绘制直线，绘制时用多段的弦来逼近弧。

技术实现要素：

本发明所要解决的技术问题是大区域真三维地理场景构建问题，利用IDL语言处理DEM模型源数据，以OSG作为图形引擎，以VS2010为平台，以标准C++作为开发语言，构建实验系统。基于球体大圆弧QTM八叉树网格剖分理论，利用全球共享数据，建立了中国大陆地区的真三维地形框架，实现了大尺度地理场景无缝集成、自适应建模与可视化，为真三维地理空间计算打下坚实基础。

本发明所采用的技术方案为：一种基于SGOG瓦块的大区域真三维地理场景自适应构建方法，包括如下步骤：

S1，SGOG瓦块网格点坐标计算与存储；

据SGOG剖分的大圆弧中分规则，由大圆弧的两个端点坐标求其弦的中点坐标，然后将其沿直线延伸到球面上，得到该网格点的坐标；

S1a，每层网格的坐标存储三个TXT文件，分别存放X、Y、Z坐标；

S1b，按方向编码顺序存储第一层网格中N个三角形的顶点坐标；

S1c，按照第一层父三角形存储的顺序，存储第二层网格中子三角形，每四个子三角形仍按方形编码排序，依次类推，存储下一层自三角形，直至存储完成；

S2，进行DEM高程匹配与计算；

S2a，在IDL平台下读取存储的SGOG最外层网格点X、Y坐标，并输入地球半径，换算为经纬度；

S2b，遍历DEM影像数据，根据影像的最大、最小行列号，换算出影像的经纬范围，判断网格点是否落在当前影像中；

S2c，网格点落在当前影像则读入影像高程值，并把落在影像内部的网格点经纬度换算为所处当前影像的行列号，根据行列号读出像元值赋给网格点高程，如果网格点不在当前影像范围内则判断下一幅影像；

S3，绘制瓦块构建真三维地形框架；

基于VS和OSG平台，在获取最外层网格点高程的基础上，按照所有可视范围的SGOG网格点坐标，按照网格可视化绘制方法绘制网格，将高程值适当放大，增强视觉效果，进而实现大区域真三维地形框架的建模与可视化；

S4，自适应建模；

根据地形起伏度大小和建模需求，设定相应的阈值，灵活采用不同剖分层次的网格进行建模可视化。以某一层为基础层次，起伏度大于阈值的地方，进一步细分；起伏度小于阈值的地方，不再细分，直到满足精度要求为止。这样既有效减少了数据量，又保留了重要的地形特征，进而达到较好的可视化效果。

S5，进行分层设色和纹理渲染；

根据网格点的不同高程值为每个网格点赋上相应的颜色值，利用OSG渲染引擎自带的函数，随着瓦块编码的遍历而绘制瓦块的各个表面并进行颜色渲染，可得到DEM场景晕渲图。当然，也可以在网格表面粘贴遥感影像纹理，构建虚拟地理场景。

进一步，所述的阈值可以采用某层网格的起伏度平均值或其相应的比例。

进一步，所述的高程插值可以采用反距离加权法、线性内插、双线性内插或最邻近像元法。

本发明产生的有益效果是：

1、基于地表大数据的一体化无缝组织，实现了考虑地球曲率的、大区域地壳真三维可视化建模，克服了基于投影空间模型的变形和裂缝等缺陷，恢复了地理空间的自然流形性质。

2、提出了基于地表起伏度阈值的自适应可视化建模方法，既减少了数据量，又保留了重要地形特征细节，保证了建模与可视化的效果。

3、本发明的方法为建立真三维数字地球平台打下了坚实基础，并可以应用于天地一体化的空间数据的组织、管理与应用。

附图说明

图1为以八分之一球体为例的球体QTM八叉树剖分效果图；

图2为整球剖分下的球体QTM八叉树第2层剖分效果图；

图3为ETP投影算法原理与编码方案；

图4为径向二叉树编码方案；

图5为SGOG编码转换空间直角坐标的流程图；

图6为空间直角坐标转换SGOG编码的流程图；

图7为真三维地理场景构建技术路线；

图8为高程匹配方法流程图；

图9为大区域真三维地形框架的第4层瓦块和第9层瓦片；

图10为大区域真三维地形框架中第6、7、8、9层的自适应可视化；

图11为图10的侧视图；

图12为自适应真三维渲染图；

图13为图12的局部放大图；

图14为真三维地理场景效果图；

图15为图14的局部放大图。

具体实施方式

以下为本发明的一种实施例。

考虑到单机的计算与存储能力，本实施例选择中国大陆地球为研究区域，以SGOG第9层剖分瓦块为建模基础，该网格球面边长约为19.5km，径向棱长约12.4km。从共享网站下载覆盖中国大陆地区的1144幅90m分辨率的DEM数据。基本参数：投影UTM/WGS-84，GeoTIF格式，3601×3601像元，高程为相对于WGS-84椭球的大地高。

本发明利用IDL语言处理DEM模型源数据，以OSG作为图形引擎，以VS2010为平台，以标准C++作为开发语言，构建实验系统。

如图7所示，首先进行SGOG瓦块网格点坐标计算与换算，通过经纬度坐标与DEM网格匹配，接着进行DEM高程匹配与计算，然后绘制瓦块构建真三维地形框架，最后进行分层设色和纹理渲染。

对SGOG瓦块网格点坐标计算与存储

SGOG采用大圆弧中分递归剖分，网格上下层之间存在着整齐的对应关系。网格点坐标计算的方法是：首先由大圆弧的两个端点坐标，求其弦的中点坐标，然后将其投影即沿径向直线延伸到球面上即可。由于在现有的研究中已经实现网格编码、网格顶点包括几何中心点的直角坐标及其经纬度三者之间的相互换算。

不同层次的大圆弧QTM网格之间有着明确的包含与被包含关系，按照一定的顺序存储各层次的网格点坐标可以方便的查找与某一网格有着包含或被包含关系的网格。由SGOG剖分思想可知，每个三角网递归剖分1次形成4个三角形，所以可推出第N层网格的一个三角形网格包含的第N+I层网格的三角形网格数量为4(I-1)，网格点数为三角形数量的3倍。第1层三角形无父三角形，末层三角形无子三角形。这样除第1层和末层之外，每个父三角形有4个子三角形。从第1个子三角形开始，每4个子三角形有1个公共父三角形。从坐标来看，父三角形的3个顶点坐标与4个子三角形12个顶点相关联，存储时保留这种对应关系，方便后续数据操作时对父子三角形的统一处理。

最外层网格具体存储方法为：

1)每层网格的坐标存储3个TXT文件，分别存放X、Y、Z坐标；

2)按方向编码顺序存储第1层网格中N个三角形的顶点坐标；

3)按照第1层父三角形存储的顺序，存储第2层网格中子三角形，每4个子三角形仍按方向编码排序。然后依次类推，存储下一层子三角形，直至存储完成。

在进行真三维DEM可视化和自适应可视化操作时，径向各层网格的存储与操作与最外层网格一致，以保证严格的对应关系。

SGOG网格高程匹配

要建立基于ESSG的大区域地形模型，必须将传统的基于投影的DEM网格高程“映射”到离散网格上。这里借助IDL语言实现，以经纬度坐标为桥梁，将WGS-84椭球近似为球体后采用最近邻像元法匹配高程。其高程匹配方法如图8所示。

具体步骤为：

1、在IDL平台下读入存储的SGOG最外层网格点X、Y坐标，并输入地球半径，换算为经纬度；

2、遍历DEM影像数据，根据影像的最大、最小行列号，换算出影像的经纬范围，判断网格点是否落在当前影像中；

3、网格点落在当前影像则读入影像高程值，并把落在影像内部的网格点经纬度换算为所处当前影像的行列号，根据行列号读出像元值赋给网格点高程，如果网格点不在当前影像范围内则判断下一幅影像。

在高程匹配过程中，各种坐标的变换是关键环节。DEM影像中像元值按行列存储，而网格点位置是以右手坐标系下的三维坐标形式表达，所以需要把网格点的三维坐标和影像像元的行列号都换算为经纬度，才能进行网格点定位，匹配像元值。在匹配完所有DEM影像数据后，由于在OSG平台可视化时使用的是右手坐标系下的三维坐标，则需要再将所有网格点的经纬度和匹配好的高程值转换成相应的三维坐标。

真三维场景可视化

DEM可视化基于VS2010和OSG平台实现。在进行最外层网格点高程的基础上，按照所有可视范围的SGOG网格点坐标，利用网格可视化技术绘制网格，实现大区域真三维地理可视化。为突出DEM高低起伏的效果，将高程扩大30倍。由于机器性能的限制，为达到较好的显示效果，灵活设定模型的网格层次。图9是利用第9层的近7万个面网格叠加第4层体网格建立的中国大陆地区的三维地形框架。

大区域真三维地理场景自适应可视化

根据中国大陆地区的地形起伏度和建模需求，设定相应阈值，灵活采用多层网格体系，例如采用第6、7、8、9层网格进行自适应可视化。具体的方法是：选取第6层次作为基层网格，然后依次递增网格层级进行自适应剖分建模。对于第6层三角形网格而言，取每个三角形三个顶点的最大高程差代表该三角形的地势起伏，取该层所有三角形地势起伏的平均值作为该层的地势起伏度，设置为阈值。大于等于此值，则认为该三角形内地势起伏较大，则继续对其进行第7层剖分处理；小于此值，则认为地势平坦，不再对三角形进行剖分处理。对于第7、8层按照上述方法处理，以第9层为最终分层网格。但是在实验过程中，发现采用上述的阈值，得到的自适应效果并不好，网格过于稀疏，因此我们将该阈值修改为原来的三分之一，即为取该层所有三角形最大高程差的平均值的三分之一作为阈值。效果如图10、图11所示。

模型高程表达精度分析

源数据、图9和图10模型的地面高程表达的特征值统计如表1所示：

表1各种DEM模型高程表达的特征值

从表中可以看出，第9层网格和自适应网格的点数都远远小于源数据，自适应网格是第9层网格的1/4。两者的高程最小值都远大于源数据的最小值，自适应网格比第9层高28m。二者高程最大值是一样的，比源数据低2262米。从平均值来看，第9层网格比源数据低830m，而自适应网格要比源数据高530m，二者相差1383m。从最值来衡量，二者的精度相当。从均值来看，二者相差较大。90mDEM网格与SGOG第9层网格属于不同的剖分类型，二者网格点在进行“匹配”时，具有一定的偶然性，它们之间的差异主要由偶然因素引起；而从自适应网格的剖分方法可知，高差越大，剖分越细，因而它与前述两种方法相比，总是“趋高”。

场景渲染

根据网格点的不同高程值为每个网格点赋上相应的颜色值，利用OSG渲染引擎自带的函数，随着瓦块编码的遍历而绘制瓦块的各个表面并进行颜色渲染，可得到中国大陆地区的DEM场景晕渲图，如图12和图13所示。青藏高原侧面的瀑布状渲染，并不代表积雪的覆盖范围，而是由于渲染技术细节不完善造成，有待进一步优化。

选择四川盆地局部地区，以地质分层示意的SGOG第8、9、10层构建三维地壳框架，表面粘贴遥感影像纹理，得到地理场景效果如图14和图15所示。