基于四叉树索引的三维数字地球空间数据组织渲染方法

摘要

本发明提出了一种基于四叉树索引的三维数字地球空间数据组织渲染方法，属于地图学与地理信息系统和虚拟现实技术领域。该方法通过将多比例尺、多投影变换下的常用空间数据格式和多光谱、多时相、高分辨率的遥感卫星影像、航空影像、不同比例尺的数字专题图统一到同一个坐标系统内，对每个元素的属性和四叉树瓦片参数调节进行操作，以四叉树瓦片的形式输出，对三维景观地图数据进行四叉树切割，将空间数据导入关系型数据库中，统一管理。通过使用该方法可将常用的矢量数据、栅格数据、高程数据、三维地图数据进行有机融合，并发布到三维数字地球原型中，缩短了数据预处理的时间并提高执行效率，为三维数字地球基础地理数据发布提供了新的集成方法。

说明书

技术领域

本发明属于地理信息系统、地图制图学和虚拟现实技术领域，具体涉及基于四叉树索引的三维数字地球空间数据组织渲染方法。

背景技术

数字地球是以计算机技术、多媒体技术和大规模存储技术为基础，以宽带网络为纽带运用海量地球信息对地球进行多分辨率、多尺度、多时空和多种类的三维描述。

数字地球主要是由空间数据、文本数据、操作平台、应用模型组成的。空间数据主要包括全球性的大、中、小比例尺的空间数据以及地球的各类多光谱、多时相、高分辨率的遥感卫星影像、航空影像、不同比例尺的各类数字专题图。数字地球的一个基本思想是将高分辨率的遥感影像覆盖全球，并通过建立覆盖全球的数字高程模型及矢量数据，建立地球景观模型，形成虚拟地球。

空间数据主要是通过各种不同类型的空间数据相结合的方式来实现的。地球的表面积接近(1.5×10¹²)平方公里，因此三维数字地球的数据量特别巨大而且数据类型复杂，数据量随着分辨率的提高呈几何级数递增，往往达到几百个G，甚至几十个TB字节。在目前的网络环境下或PC机上，如此庞大的数据难以达到实时传输和显示的要求，故海量地理数据的传输和三维场景的实时绘制成为数字地球研究中的重要内容。计算机存储设备的更新和价格的下降为大容量的数据存储提供了支持，对基于网络的系统来说，合理组织空间数据，提高海量数据的访问速度并尽可能的减少不必要的数据在网络上的传输，对提高全球三维场景的显示速率具有重要意义。

全球地形数据量非常巨大，例如，全球1∶25万的地形数据总的数据量达20多个G，再加上影像，其数据量可想而知。为了降低计算机内存的消耗，加快计算机处理时间和节省网络传输时间，能够最快的访问多分辨率超大规模的地理数据，必须设计一个合理的数据库存储结构用来灵活的管理从粗到细的数据，这样的适应性数据库结构必须使用空间索引来实现。通常将数据分层组织，每一层的数据又分割成小的数据块，这种组织数据的方法通常称为四叉树结构或金字塔数据结构。图1说明了四叉树的结构层次关系。

参考文献

[1]承继成，李琦，易善桢.国家空间信息基础设施与数字地球[M].北京：清华大学出版社，1999.

[2]龚键雅，杜道生，李清泉，等.当代地理信息技术[M].北京：科学出版社，2004：87～88.

[3]方裕，周成虎，景贵飞，等.第四代GIS软件研究.中国图像图形学报，2001b，6A(9)：817～823.

[4]宋关福，钟耳顺.组件式地理信息系统研究与开发，图像图形学报，1998，4：314.31.

[5]童晓冲，贲进，张永生.全球多分辨率数据模型的构建与快速显示.测绘科学.2006.1.

[6]白建军，赵国胜，陈军基于椭球面三角网的全球地形可视化，武汉大学学报，2005.05.

[7]李德仁.信息高速公路、空间数据基础设施与数字地球.测绘学报，1999.

[8]李德仁.论世纪遥感与的发展.武汉大学学报，2003.02.

[9]吴艳兰.DEM的可视化研究及若干应用.武汉测绘科技大学学位论文，1998.

[10]彭辉.基于MicroStation的数字高程模型的建立及其应用研究，武汉测绘科技大学学位论文，1998.

[11]薛勇，王剑秦，郭华东.数字地球网格计算雏议[J].遥感学报，2004.8(I).

[12]杨永崇，郭达志.基于地理坐标的数字地图量算模型[J]，西南交通大学学报，2005年6月，Vol 40，No3.

发明内容

本发明的目的在于提出了一种基于四叉树索引的三维数字地球空间数据组织渲染方法。该方法可以将全球性的多比例尺、多投影变化下的常用空间数据格式以及地球的各类多光谱、多时相、高分辨率的遥感卫星影像、航空影像、不同比例尺的各类数字专题图统一到同一个坐标系统内，并可进行对每个元素的属性调节、四叉树瓦片参数调节等操作，最终将可视化数据以四叉树金字塔瓦片形式输出，将三维景观地图数据进行四叉树切割，将空间数据的几何数据与属性数据导入关系型数据库中，统一管理。本发明为最终数据提供了一套检测框架，确保了数据可以稳定、正确的提供给数字地球平台。

本发明的技术方案，基于四叉树索引的三维数字地球空间数据组织渲染方法，包括如下设计思路(或者如下基本步骤)：

1)采用数字地球数据分层分块机制

2)生成瓦片并对其进行渲染

3)对多源数据进行融合处理

4)采用三维数字地球原型数据检测工具

本发明的有益效果在于：通过使用该方法可以将常用的矢量数据、栅格数据、高程数据、三维地图数据等进行有机融合，并以四叉树金字塔的切分存储方式发布到三维数字地球原型中，显著缩短了数据预处理的时间并提高执行效率，为三维数字地球基础地理数据发布提供了新的集成方法。

附图说明

图1四叉树的结构层次关系

图2笛卡尔坐标系

图3瓦片的宽和高

图4 36度划分

图5不同级别划分示意图

图6各级别四叉树瓦片分块表示

图7第n层四叉树到第n+1层四叉树的分割规则

图8瓦片预载入与过渡显示处理

图9地图数据来源

图10在3Ds Max中三维景观地图渲染过程

图11单幅地图切割流程

图12超大地图多幅关联切割流程

图13线性矢量图层与DEM和影像的层次结构关系

图14支持大地坐标到兰伯特投影坐标

图15支持Jpeg2000栅格数据

图16支持DEM数据

图17OpenGIS的简单要素数据模型

图18瓦片渲染方法工具原型截图

图19瓦片渲染方法工具原型导出数据到三维数字地球原型

图20多层栅格与矢量数据相融合

具体实施方式

下面结合说明书附图介绍本发明的设计思路。

1  数字地球数据分层分块机制

1.1地表瓦片划分机制

采取笛卡尔坐标，原点(X＝0，Y＝0.)在投影坐标左下方，即南极点(-90，-180)(以经纬度表示)，见图2和图3。

如图3，使用“Level Zero Tile Size”来决定每一个瓦片宽和高的大小(所有的瓦片都是正方形。标准的level zero tile size还没有制定，但是它必须满足能被180整除。

以地球360度经度和180度纬度为标准，第一层以36度划分，见图4。

共分为(360/36)*(180/36)＝50层，第二层以18度，第三层以9度以次类推，见图5和图6。

1.2瓦片定位机制

定义：

A.地图等级为从1开始、行列数为从0开始数；

B.地图范围的极值用XMin，XMax，YMin，YMax表示；

C.用xΔ表示横坐标差，yΔ表示纵坐标差；

D.四叉树地图瓦片的命名采用包含地图等级、行列数的方式，如“zoomx_row_column”，其中x表示地图等级，row和column分析表示当前瓦片所处的行数和列数；

1.2.1地图数量计算

设第一个等级地图的为m行n列，则第level级的地图行列数如下：

行数：m×2^level-1，level≥1

列数：n×2^level-1，level≥1

1.2.2地理坐标与瓦片位置的换算

1.2.2.1由瓦片位置计算地理范围

输入：行列对(i，j)

输出：Extent范围(currentXMin，currentXMax，currentYMin，currentYMax)

第level级下第i行j列位置切片的坐标范围由如下公式进行计算：

$\mathrm{currentXMin}=\mathrm{XMin}+j\times \frac{\mathrm{\Delta x}}{m\times 2^{\mathrm{level}-1}}$

$\mathrm{currentXMax}=\mathrm{XMin}+(j+1)\times \frac{\mathrm{\Delta x}}{m\times 2^{\mathrm{level}-1}}$

$\mathrm{currentYMin}=\mathrm{YMax}-(i+1)\times \frac{\mathrm{\Delta y}}{n\times 2^{\mathrm{level}-1}}$

$\mathrm{currentYMax}=\mathrm{YMax}-i\times \frac{\mathrm{\Delta y}}{n\times 2^{\mathrm{level}-1}}$

1.2.2.2由坐标计算瓦片位置

输入：(x，y)坐标

输出：i，j(图片的行列数)

其中，

$j=\sqrt{(x-X_{\mathrm{Min}})/\mathrm{\Delta x}\times m\times 2^{\mathrm{level}-1}}$

$i=\sqrt{(Y_{\mathrm{Max}}-y)/\mathrm{\Delta y}\times n\times 2^{\mathrm{level}-1}}$

采用四叉树的数据结构，参照图4，Level 0把地球按36度划分成50个区域，第0层的每一块分裂为四个子树，这一层为Level 1层，依次推算如图7。

每个块都有四个子树，包括本身信息，数据定义为：

      double West；     //四个方向的平面笛卡尔坐标值

      double East；

      double North；

      double South；

      Angle CenterLatitude；//瓦片中央经纬度

      Angle CenterLongitude；

      double LatitudeSpan；

      double LongitudeSpan；

      int Level；//层次Level

      int Row；      //行，列

      int Col；//子树

      QuadTile northWestChild；

      QuadTile southWestChild；

      QuadTile northEastChild；

      QuadTile southEastChild；

每一块对应当前Level的一个512×512的纹理，每一块的每条边分为40等份，生成41×41或43×43个顶点，再结合高度信息实时绘制。

由于分层分块的结构，对地形数据的加载和实时渲染时，只需要从服务器下载需要层次，区域的纹理图片、高程信息，还可以对有些图片进行预处理成三维果的图片直接贴图，不需要高程信息。

过渡处理：

Level(n+1)级未下载完未能处理的图片，仍然显示Level n级的效果，见如图8。

2  生成瓦片

2.1三维瓦片数据

地图的预生成是将指定范围的地图按照指定尺寸和指定格式(如TIFF、JPEG、PNG等)切成若干行及列的正方形图片，切图所获得的地图切片叫做瓦片(Tile)，瓦片的获取流程如图9所示。渲染瓦片只能得到二维矢量数据、影像数据或者数字高程数据对应的栅格图，三维景观地图的切图需要借助AutoCAD、3DMax等工具进行三维建模，并通过VRay插件进行调节后再渲染整幅三维地图。

三维电子地图中的三维实际是一种观察面丰富的2.5维，是由三维模型生成的固定视角静态渲染图片。制作前要根据区域房屋朝向特点，确定摄像机角度，在三维模型制作软件设定相机、灯光效果、制作底图，形成标准场景，供制作人员统一使用。对那些被遮挡、看不见或位于次要位置的物体，可以省略、简化。三维场景中的点、线、面城市实体信息图层，在模型制作时分层制作，三层模型叠加形成城市三维景观。图10所示为在3DSMax中三维模型渲染成景观图片的实现过程。

2.2瓦片切割

由于计算机硬件的限制，在技术实现层面.Net类库对于图像体积具有限制，对于载入的单幅的地图长宽幅度限制在15000*15000以内，因此将地图的切割分为单幅地图切割与超大地图多幅关联切割两类：

1、单幅地图切割流程：

步骤1：对将切割范围的原始值进行修正，使其长高值均为瓦片的倍数；

步骤2：从左上角(西北)开始对此固定范围进行切图，后续的不同等级地图之间采用四叉树数据结构，第level级上的一张瓦片(Tile)到第level+1级将裂变成4张，这种结构有助于切图和显示，但得出的地图没有固定的比例尺，比例尺随地理纵坐标变化，因此进行地理量算时不是根据比例尺而是根据地理坐标直接计算。地图命名规则为第N级地图图片命名为zoomN.png。

第x级chunkrow行chunkcol列的瓦片命名规则为zoomx_chunkrow_chunkcol.png，切割后瓦片的命名规则为zoomx_chunkrow_chunkcol\x_chunkrow_chunkcol。

其中各变量定义：

X：当前缩放级别；

chunkrow：当前瓦片在地图中的行数；

chunkcol：当前瓦片在地图中的列数；

2、超大地图多幅关联切割流程：

步骤1：将超大地图分割为相互关联的若干小地图，小地图命名规则为uperchunk_row_col\zoomN.png，其中row为当前小地图在所有切割后的小地图中所占行位置，col为当前小地图在所有切割后的小地图中所占列位置。

步骤2：对小地图采用单幅地图切割流程进行切割

切割后瓦片的命名规则为superchunk_row_col|zoomx_chunkrow_chunkcol\x_chunkrow_chunkcol。

其中各变量定义：

Row：当前小地图在所有切割后的小地图中所占行位置；

Col：当前小地图在所有切割后的小地图中所占列位置；

X：当前缩放级别；

chunkrow：当前瓦片在地图中的行数；

chunkcol：当前瓦片在地图中的列数；

单幅地图切割流程/超大地图多幅关联切割流程，如图11和图12。

3  多源数据融合

3.1融合数据图层

随着视点远近的不同，三维地形的清晰度实时发生变化，这时系统自动调入不同分辨率的DEM和影像数据来构造地表模型。当地形从一个细节层次向另一个细节层次过渡时，系统同时调用多幅地形数据显示，地形面片的分割将发生变化。如果地物没有与地形相融合，就会造成空间矢量实体出现“下陷”和“抬高”现象。为了加快显示的速度和图形质量，矢量图层的数据与DEM和影像一样要分级管理，如图13。

3.2支持投影变换

通过读取.shp文件相对应的.prj文件，并以WKT格式标准进行解析，如图14。

参数包括：

坐标系统是按OpenGIS WKT(Well known Text)的方式来描述的。它包括如下信息：

A.一个总的坐标系的名称。

B.一个地理坐标系统的名称。

C.大地参照系。

D.参考椭球体，椭球半长轴和。椭球扁率的导数(即a/(a-b))。

E.一个0度经线的名车以及它于本初子午线的偏离程度。

F.一个投影类型。

G.一个投影参数的列表。

H.一个单位名称以及转化为米或者弧度的转化因子。

I.轴的名称和排序。

3.3多种矢量和栅格数据类型支持

以OpenGIS标准作为统一接口，可以同台支持不同数据，如

■矢量数据：Shape数据、mapinfo数据、gml数据、geojson数据

■栅格数据：Geotiff数据、ecw数据、jp2数据(见图15)

■数字高程数据：DEM数据(见图16)

■GPS数据：.GPX格式

使用GDAL库来实现了对Shape数据、mapinfo数据、gml数据、geojson数据等矢量数据的支持。

矢量数据是基于直角坐标系统，用点、线、多边形描述地理要素得数据模型或数据结构。每个地理要素有一系列有顺序的坐标(x，y)坐标描述，这些要素与属性相结合。我们设计了严格遵循OGC OpenGIS的简单要素数据规范的Geometry模型，Geometry中定义的是矢量图形中的基本几何形状，以及一些简单的空间关联关系分析功能。最重要的几何形状是点、线、面，以及由此派生出来的多点、多线、多面、曲线、直线、环线、多边形等，见图17。

用两种方式描述栅格位置(象素坐标或线坐标)和地理参考坐标的关系。第一种也是最普遍应用的一种是仿射变换。

仿射变换由六个系数构成，它用如下的关系式将象素/线坐标映射到地理参考空间：

Xgeo＝GT(0)+Xpixel*GT(1)+Yline*GT(2)

Ygeo＝GT(3)+Xpixel*GT(4)+Yline*GT(5)

对于向下排列的图象，GT(2)和GT(4)都是零，GT(1)是象素宽，GT(5)是象素高，(GT(0)，GT(3))表示栅格左上角象素的左上角坐标。

注意：上面的象素和线坐标都是从左上角象素的左上角坐标(0.0，0.0)到右下角象素的右下角坐标(象素宽，象素高)。左上角象素/线的中心位置则是(0.5，0.5).

定义：波段

一个波段的光栅数据对应GDAL中的GDALRasterBand类。它描述了单个波段的band/channel/layer。它并不是一次描述整个影象。例如，一个24bit的RGB影象在一个dataset中一般被描述为三个波段，分别对应red/green/blue三中颜色。

每个波段解释具体为：

GCI_Undefined：默认，未知信息。

GCI_GrayIndex：为一个独立的灰度影象。

GCI_PaletteIndex：this raster acts as an index into a color table。

GCI_RedBand：RGB或RGBA影象的红色波段。

GCI_GreenBand：RGB或RGBA影象的绿色波段。

GCI_BlueBand：RGB或RGBA影象的蓝色波段。

GCI_AlphaBand：RGBA影象的alpha通道。

GCI_HueBand：HLS影象的色调。

GCI_SaturationBand：HLS影象的饱和度。

GCI_LightnessBand：HLS影象的光强。

GCI_CyanBand：CMY或CMYK影象的cyan波段。

GCI_MagentaBand：CMY或CMYK影象的magenta波段。

GCI_YellowBand：CMY或CMYK影象的yellow波段。

GCI_BlackBand：CMYK影象的black波段。

支持数据属性调节：支持数据属性调节

3.4四叉树瓦片参数调节

将瓦片按四叉树参数调节后导出，并导出参数配置文件，格式如：

<MapParam>

  <tilesize>128</tilesize>

  <ColumnNum>4</ColumnNum>

  <RowNum>2</RowNum>

  <centerX>-243383.913445131</centerX>

  <centerY>4376213.31793891</centerY>

  <ZoomParam>2</ZoomParam>

  <MaxLevel>4</MaxLevel>

  <Xmin>-4663419.25215648</Xmin>

  <Xmax>4176651.42526622</Xmax>

  <Ymin>2166195.64858323</Ymin>

  <Ymax>6586230.98729458</Ymax>

  </MapParam>

3.5空间数据库的建立

A.数据库连接

包含：①Sql方式验证；②Windows身份验证。

B.新建数据库，包括数据库文件(.MDF)和日志文件(.LDF)。

C.空间数据库部署：在指定数据库中生成字段“BLUE.GEOMETRY_COLUMNS”和“BLUE.SPATIAL_REF_SYS”，分别用于存放几何图层与空间WKT信息。

D.矢量数据导入

将shape数据导入到数据库中。shp文件生成固定字段oid，the_geom，zID，the_geom_Envelope_MinX，the_geom_Envelope_MinY，the_geom_Envelope_MaxX，the_geom_Envelope_MaxY，分别存放流水号、几何特征数据、图层号、包围盒横坐标最小值、包围盒纵坐标最小值、包围盒横坐标最大值、包围盒纵坐标最大值。由.dbf文件生成属性信息字段，例如：“城市”、“人口”。

3.6地理数据批量搜索

基于Geocoding的经纬度搜索，通过输入地理位置搜索到对应经纬度。发送格式为“http://ditu.google.com/maps/geo？output＝xml&q＝地理位置”的请求，对返回的xml格式文件进行解析，其中的point.coordinates数据即为经纬度值。通过使用数据集与数据适配器，基于sqlserver2005进行批量的数据搜索。

4  由本方法对地图进行编辑并导出数据至数字地球的演示，见图18、图19和图20。