CN103455624B - 一种轻量级全球多维遥感影像网络地图服务实现方法 - Google Patents

Abstract

一种轻量级全球多维遥感影像网络地图服务实现方法，包括构建全球地理空间金字塔模型，根据模型建立Tile文件编码体系；对遥感数据进行预处理，主要包括投影转换、文件格式转换等，将遥感数据切割成瓦片并进行文件编码；根据分布策略将Tile文件进行分组并生成相应的目录结构，将不同的图片组更新到不同的Web服务器上；分析客户端请求的地图可视范围计算出在此范围内的所有图片的存储路径，并将图片下载返回给客户端。本发明可以为遥感数据中心(如中国资源卫星应用中心)提供一种轻量级的在线地图服务实现方法，能够更为方便地为用户提供遥感数据服务。

Description

一种轻量级全球多维遥感影像网络地图服务实现方法

技术领域

本发明涉及一种轻量级全球多维遥感影像网络地图服务实现方法，主要用于遥感数据提供者通过网络地图方式为用户提供遥感影像数据在线服务。

背景技术

卫星遥感已经成为对地观测的主要数据来源，因为卫星能够快速地提供具有一致性数据质量的覆盖全球的遥感数据。通过卫星遥感，人类在过去几十年更为广泛地了解了关于地球以及地球是如何工作的科学知识。这些知识的应用，从天气预测到资源管理，已经证明遥感对于人类社会的巨大益处。仅仅是在过去20年的时间里，遥感数据使得五天内的天气预测的准确性从65%提高到85%。遥感数据对于灾害预防和灾后的快速响应是非常重要的，2004年的亚洲海啸和2005年的卡特里娜（Katrina）飓风就是很好的例子。遥感也已经成为自然资源开发与保护的重要工具，传统上认为自然资源问题主要是局部和区域性的问题，但过去几十年区域性问题已经变得越来越来重要，例如酸雨问题已经显示出局部的排放物是如何影响远距离地区的环境，区域性的污染源是如何影响全球的空气质量。自然资源和环境保护是经济可持续发展和地区与全球安全的基础，遥感数据对于监测和有效管理这些资源是非常重要的。商业利用遥感数据来提高它们运作的效率，例如交通运输、发电和航运等。世界经济的成熟将日益依赖遥感数据来支撑商业的发展和发展中国家经济建设。总而言之，卫星遥感已经成为人类社会繁荣和可持续发展的重要支撑技术。

随着航天技术和卫星遥感技术的不断进步，遥感影像数据的时空和光谱分辨率得到显著提高，使得遥感数据量呈几何指数增长。同时，越来越多的卫星提供了不同地理区域、不同分辨率、不同时相和不同类型（多光谱，全色和SAR）的遥感影像数据，这些多数据源、多分辨率、多时相和多数据类型的遥感影像数据统称为多维遥感影像数据。另一方面，随着Web技术的不断发展，越来越多的遥感数据通过Web方式传递给公众用户。特别是2005年Google公司推出了一种全新的遥感影像数据在线服务Google Map，Google Map彻底改变了与地图交互的方式，即传递给用户的不再是整张地图，而是一个个影像Tile（瓦片），从而大大提高了网络地图服务的响应时间。其它商业公司，如微软和雅虎等，也推出了类似的网络地图服务。但这些商业地图服务大多采用大型关系数据库系统，如微软的VirtualEarth和NASA的OnEarth分别采用SQL Server和Oracle，或者采用专有的分布式文件系统和非关系数据库，如Google Map采用GFS和Bigtable，Yahoo Map采用开源的Hadoop和HBase。这些重量级地图服务的部署和管理都是非常复杂的，底层的存储系统特别针对超大规模的集群而设计。此外，这些商业网络地图服务系统不支持遥感数据的多维特征（如多数据类型，多时间和多光谱特征等）。而很多遥感数据提供者，如中国资源卫星应用中心，需要一种轻量级的网络地图服务实现方法来向公众提供他们的遥感数据，即系统运行在通用的硬件和软件平台之上（Intel X86服务器，Windows或Linux操作系统和Apache Web服务器等），系统的规模相对较小，系统的构建成本低并且使用方便，系统对于数据访问是高效的。

发明内容

为了实现全球多维遥感数据的网络快速发布与高效交互式地图访问，本发明提出了一种轻量级的全球多维遥感影像网络地图服务实现方法。

本发明的技术方案为一种轻量级全球多维遥感影像网络地图服务实现方法，包括以下步骤：

步骤1，通过自顶向下逐级划分全球地理空间范围构建全球地理空间金字塔模型，采用全球四叉树结构表示全球地理空间金字塔模型，并建立相应Tile文件编码体系，其中Tile表示瓦片；

步骤2，对遥感影像进行预处理，包括投影变换为一致的投影方式和格式转换为一致的图像格式；

步骤3，根据全球四叉树结构对遥感影像进行切分得到Tile图片，并根据Tile文件编码体系对每个Tile图片相应存储数据的Tile文件进行编码获得文件名；

步骤4，根据数据分布策略对Tile文件进行分组，包括确定每个Tile文件存放的Web服务器，存储在相同服务器上的Tile文件为一组；

步骤5，生成每个分组上所有Tile文件的分层目录结构；

步骤6，根据步骤5所得分层目录结构将各组的Tile文件更新到相应的服务器上；

步骤7，当用户通过客户端查看地图时，计算出地图可视范围内需要下载的Tile图片相应Tile文件文件名，并根据步骤4的数据分布策略计算出每个相应Tile文件实际存储的Web服务器地址；

步骤8，根据步骤7所得Web服务器地址，按照文件名从相应Web服务器下载相应Tile文件到客户端；

步骤9，在客户端根据步骤8获取的Tile图片进行显示。

而且，全球地理空间金字塔模型每一层的遥感影像来自相同或不同的数据源。

而且，所述一致的投影方式采用球形墨卡托投影，所述一致的图像格式采用JPEG格式。

而且，所述全球四叉树结构，四叉树中的每个节点表示一个虚拟瓦片，设全球四叉树结构有n层，顶层用一个虚拟瓦片表示全球地理空间范围，第2层用4个虚拟瓦片表示全球地理空间范围，第3层用16个虚拟瓦片表示全球地理空间范围，依次类推…第n层由4的n-1次方个虚拟瓦片拼接而成。

而且，所述Tile文件编码体系为，每个Tile图片相应存储数据的Tile文件根据Tile图片的全球四叉树网格坐标<level，<row,col>>生成文件名，其中level表示Tile图片在全球四叉树结构中的层数，row表示Tile的X坐标值，col表示Tile的Y坐标值。

而且，步骤3中，对每一层遥感数据分别处理方式如下，

（1）根据遥感数据的球形墨卡托投影坐标范围叠加到全球四叉树结构中该层相应虚拟切割网格上，然后沿网格线将遥感数据切割成Tile图片；

（2）根据遥感数据的最小包围矩形的左上角XY坐标和右下角XY坐标，按照Tile文件编码体系对有效Tile图片进行编码，所述最小包围矩形是遥感数据的影像图幅在虚拟切割网格中所覆盖的以网格为单位的最小外接矩形，所述有效Tile图片是指在最小包围矩形内被遥感影像完全覆盖的Tile图片。

而且，步骤4中，根据数据分布策略对Tile文件进行分组的实现方式如下，

（1）从全球四叉树结构的第二层开始采用Hilbert空间填充曲线进行网格填充，金字塔第二层采用一阶Hilbert曲线进行填充，第三层采用二阶Hilbert曲线进行填充，依次类推…第n层采用n阶Hilbert曲线进行填充；

（2）计算出每一层每一个网格的Hilbert编码，对每个Tile文件根据Hilbert编码映射到一个Web服务器。

而且，步骤5中，生成的每个分组上所有Tile文件的分层目录结构如下，

预定义表示不同数据图层的一级目录，每个一级目录下定义n个二级目录，每一个二级目录下的目录层次结构根据遥感影像的特征设置。

而且，步骤7中，计算地图可视范围内需要下载的Tile图片相应Tile文件文件名时，根据InnerDiv相对OuterDiv的位置来计算需要下载的Tile图片的网格坐标，InnerDiv表示容纳相应层数的所有Tile图片的虚拟网格区域，OuterDiv表示地图可视范围。

而且，计算出每一层每一个网格的Hilbert编码的实现方式如下，

首先建立以下快查表，

表2子曲线方位快查表

表3象限划分过程中Hilbert编码快查表

表4Hilbert编码尾码快查表

然后执行步骤如下，

(1)对当前层中相应的Hilbert曲线，确定父曲线的方位，计算目标Tile图片所属目标象限，然后根据所得目标象限和父曲线的方位，通过查询表2获得位于目标象限的子曲线的方位；

(2)判断是否能够将目标象限划分为4个子象限，是则进行划分并进入(3)，否则根据Tile图片所处目标象限和父曲线的方位，通过查询表4执行Hilbert编码的最后处理，然后返回处理结果，结束流程；

(3)确定目标Tile图片属于哪个目标子象限，根据父象限的曲线的方位和目标子象限的位置，确定目标子象限内的曲线的方位；之后，根据子曲线所处象限和父曲线的方位，通过查询表3得到目标子象限曲线的初始顺序值；

(4)跳转到(2)，以(3)确定的目标子象限作为当前目标象限，进一步划分子象限。

本发明与现有技术相比的优点在于：

(1)本发明的全球地理空间金字塔模型GGPM不同于传统影像金字塔模型，即从已有的遥感影像数据通过逐层抽取来构建一个多分辨的金字塔结构（单数据源金字塔模型），GGPM是通过自顶向下逐级划分全球地理空间范围来构建金字塔结构，金字塔的每一层都对应于整个地球平面，金字塔每一层的遥感影像数据可以来自不同的卫星和传感器平台，GGPM是一种多数据源的金字塔模型。

(2)本发明的GGPM具有固定的金字塔层次，结构清晰简单，不同数据源在相同金字塔层次上的空间分辨率相同。如果添加新的遥感影像数据，只需根据其地理范围和分辨率特点映射到相应的金字塔层上即可。

(3)本发明的全球四叉树结构GQ能够有效地组织全球多维遥感数据，GQ提供了统一的Tiling分块规则，从而可以保证不同景的遥感影像在客户端能够无缝拼接。

(4)本发明的Tile分层目录存储结构能够有效支持遥感数据增量式更新和目标Tile快速检索。

(5)本发明的Hilbert曲线Tile分布方法HCTD能够保证金字塔每一层相邻的Tile尽可能的分布到不同的Web服务器上，从而能够并行地从Web服务器上读取Tile数据，大大减少了数据读取的时间。HCTD还能使得服务器间的负载得到均衡，即每一个服务器上的Tile图片数量基本相同。

附图说明

图1为本发明实施例的流程图；

图2为本发明实施例的全球地理空间金字塔模型GGPM示意图；

图3为本发明实施例的全球四叉树GQ示意图；

图4为本发明实施例的瓦片间相邻关系示意图

图5为本发明实施例的瓦片间父子关系示意图

图6为本发明实施例的Tile网格编码示意图；

图7为本发明实施例的Tile分层目录存储结构示意图；

图8为本发明实施例的最小包围矩形示意图；

图9为本发明实施例的Hilbert曲线Tile分布方法HCTD示意图；

图10为本发明实施例的遥感图像网格坐标定位；

图11为本发明实施例的HCTD象限划分。

图12为本发明实施例的Hilbert曲线方位图(曲线类型)，其中图12a为方位1的示意图，图12b为方位2的示意图，图12c为方位3的示意图，图12d为方位4的示意图。

图13为本发明实施例的Hilbert曲线编码顺序示意图。

图14为本发明实施例的需要显示瓦片的示意图。

图15为本发明实施例的Hilbert曲线整体走向示意图，其中图15a为整体走向是方位1的示意图，图15b为整体走向是方位2的示意图，图15c为整体走向是方位3的示意图，图15d为整体走向是方位4的示意图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

如图1所示，本发明实施例具体步骤如下：

1、构建全球地理空间金字塔模型

该方法基于全球地理空间金字塔模型（Global Geospatial Pyramid Model，GGPM）来组织全球多维遥感数据，GGPM不同于传统影像金字塔模型从已有的遥感影像数据通过逐层抽取来构建一个多分辨率的金字塔结构，GGPM是通过自顶向下逐级划分全球地理空间范围来构建金字塔结构，金字塔的每一层都对应于整个地球平面，如图2所示，从顶向下依次构建了金字塔第1层、金字塔第2层、金字塔第3层…金字塔第19层…。具体实施时本领域技术人员可根据需要指定金字塔的层数，实施例以19层为例进行说明。GGPM是一种多数据源的金字塔模型，金字塔的每一层的遥感影像数据可以来自不同的卫星和传感器平台，例如，金字塔第10层采用来自卫星A的数据，金字塔第18层采用来自卫星B的数据。GGPM采用了统一的金字塔策略，具有固定的金字塔层次，结构清晰简单，不同数据源在相同金字塔层次上的空间分辨率相同。GGPM能够实现遥感影像的无缝漫游，在从一个分幅地区漫游到另一个分幅地区的时候，只需要在当前金字塔层上检索相应的遥感影像数据，从而简化了调度算法。对于GGPM，如果添加新的遥感影像数据，只需根据其地理范围和分辨率特点映射到相应的金字塔层上即可。

实施例设计一种全球四叉树结构（Global Quadtree，GQ）来表示GGPM，如图3所示，然后根据GQ提供的统一Tiling分块原则对遥感影像数据进行Tiling分块，GQ的每一个节点对应一个256×256大小（具体大小本领域技术人员可自行设定）的正方形，即四叉树中的每个节点表示一个虚拟瓦片（Virtual Tile,VT），虚拟瓦片代表一个地理空间范围，它包含数据来源、获取日期和数据类型等属性。GQ的顶层用一个VT来表示全球地理范围，GQ的第2层用4个VT来表示全球地理范围，它们将地球平面均匀地划分为4个部分，分别为西北NW、东北NE、西南SW和东南SE，第3层用16个VT来表示全球地理范围，依次类推…第n层由4的n-1次方个虚拟瓦片拼接而成，GQ每一层VT的数量是上一级的4倍。将多维遥感影像数据按照统一的规则切割成256×256大小Tile后，再根据Tile数据分布规则将Tile数据分布到不同的Web服务器上，当客户端请求某一地理范围的遥感数据时，可以并行地从不同的Web服务器上获取Tile数据，大大减少了数据延迟，提高了网络地图服务的响应时间。

实施例还基于上述全球地理空间金字塔模型和全球四叉树建立Tile文件编码体系，这个编码体系能够体现每个Tile图片在整个GQ中的相对位置，这个相对位置主要体现在两个方面：一是兄弟关系，即在同一分辨率层次下Tile图片的相邻关系，如图4所示，瓦片（X,Y）的相邻瓦片有（X,Y-1）、（X-1,Y）、（X+1,Y）、（X,Y+1）；另一个是父子关系，即Tile图片的分辨率层次关系，如图5所示，瓦片（X,Y）作为父瓦片，子瓦片有（2X,2Y）、（2X,2Y+1）、（2X+1,2Y）、（2X+1,2Y+1）。GQ的每一层Tile构成了一个网格矩阵，左上角Tile的坐标为<0,0>，X坐标轴沿水平方向从左至右依次递增，Y坐标轴沿垂直方向由上至下依次递增，如图6所示，图中X和Y后面的数值0、1、2、3分别表示在X轴和Y轴上的坐标值(正整数)。Tile的全球四叉树网格坐标表示为<level，<row,col>>，其中level表示Tile在GQ中的层级数(为自然数)，row表示Tile的X坐标值，col表示Tile的Y坐标值，row和col为非负整数，取值范围为[0，2^{z -1}-1]。

设已知给定Tile的编码为<z，<tx,ty>>，z为缩放级别，tx为X坐标值，ty为Y坐标值，与它相邻的Tile的编码计算如下所示：

左相邻Tile编码：

右相邻Tile编码：

上相邻Tile的编码：

下相邻Tile的编码：

四个子Tile的编码：

父Tile的编码：为向下取整

由上述计算公式可以看出，在计算同一缩放级别相邻Tile编码时，只涉及到简单的加减法运算，这两种运算的计算速度是最快的。在计算相邻缩放级别的Tile编码时使用了乘除法，但都是2倍率运算，可以用位运算代替，即2x=x<<1,x/2=x>>1，而位运算的计算速度和加减法差别不大。由此可见，Tile全球四叉树网格编码能够有效地实现Tile快速的检索。对于给定的范围查询，只要计算出边界Tile的坐标就能够获得请求范围内所有Tile的坐标。

为了寻址每一个Tile文件，本发明实施例进一步提出采用文件名来表示Tile全球四叉树网格编码，这种方式的优点在于可以直接使用文件名作为Tile的索引键名，不需要为Tile建立单独的索引文件，另外一个优点是便于数据的局部增量更新，因为只需要覆盖文件名相同的Tile，而其它的Tile数据不会受到影响。与虚拟瓦片的全球四叉树网格坐标<level，<row,col>表示一致，Tile文件的命名方式为X_rowY_colZ_level.ext，row和col为Tile的网格坐标X坐标值、Y坐标值，level为Tile所在的金字塔缩放级别，ext为文件的扩展名，例如JPEG或PNG。

2、遥感影像数据预处理

对于新增的遥感影像数据都采用一致的处理方式。实施例包括以下子步骤：

（1）将遥感数据进行投影变换，即将UTM转换为SMP；因为目前大部分卫星遥感影像采用的是UTM投影，投影变换后根据遥感数据的空间分辨率查询表1来确定它所在的缩放级别，如果空间分辨率正好等于表中的某一个空间分辨率则直接确定它的缩放级别，如果介于两个分辨率之间，则需要对遥感数据进行分辨率调整。

（2）对遥感数据进行格式转换，将GeoTiff格式转换为JPEG格式。

遥感影像数据预处理主要包括投影变换、分辨率调整和图像格式转换等图像处理操作，而这些处理操作是非常费时的。

由于全球地理空间金字塔模型GGPM是以全球地理平面作为金字塔构建的基础，因而地图投影的选择是非常重要的。为了实现全球多维遥感数据的统一地图展现，需要一种能够满足全球制图的投影方式，将所有的遥感影像数据统一到同一个投影下，实现数据的自然统一。目前卫星遥感影像数据主要采用的是UTM投影（横轴等角割圆柱投影），这种投影方式使得属于不同带（UTM zone）的影像不能拼接在一起，因此需要将UTM投影转换为SMP投影。由于墨卡托投影系统能够统一和无缝的描述地球球体，没有跨投影的问题，非常适合全景式展现。球形墨卡托投影(Spherical Mercator Projection,SMP)是墨卡托投影的一种变体，它在地图投影时采用球体映射。SMP可以简化计算并且不影响可视化应用的效果，因为仅仅在Y轴方向大约有0.33%的比例失真，而当比例尺越大和地物更详细的时候，这种差别基本可以忽略。地理坐标系统采用基于WGS84基准面的经纬度坐标系统，因为WGS84非常适合全球数据集。

为了实现SMP用于显示正方形的遥感影像数据，可以对地球的纬度进行一些限制。由于赤道半径为6378137米，则赤道周长为2×pi×r=20037508.3427892，其中，pi表示圆周率，r表示半径。因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。当纬度接近两极，即90°时，Y值趋向于无穷，因此通常把Y轴的取值范围也限定在[-20037508.3427892,20037508.3427892]之间，这将克服Y坐标在极点接近无穷大的问题并且实现显示正方形遥感数据的目的。通过Y轴取值范围为[-20037508.3427892,+20037508.3427892]反计算可得到纬度范围为[-85.05112877980659，+85.05112877980659]。GGPM实际的全球地理坐标范围是<-180,-85.05112877980659,+180,+85.05112877980659>。

具体实施时，重投影可采用现有技术中GDAL库（Geospatial Data AbstractionLibrary）提供的gdalwarp命令实现。完成重投影后还需要根据遥感影像的空间分辨率确定它的金字塔层级（缩放级别），如果遥感数据空间分辨率正好等于某一个金字塔层级，则不需要进行分辨率调整，否则需要进行分辨率调整，使之等于某一个金字塔层级的空间分辨率，为了保证图像质量，一般是将分辨率降到与它最接近的层级，例如将分辨率为3米的遥感影像调整为4.8米的分辨率（实施例中金字塔第16层）。最后扩充图像的边缘，使切割后的边缘处的Tile大小都为256×256像素。

通过全球四叉树GQ可以确定金字塔模型中每一层的空间分辨率（米/像素），空间分辨率计算公式如下所示：

其中zoom为金字塔缩放级别，res_zoom表示缩放级别为zoom的空间分辨率，tilesize为Tile的像素大小（实施例中tilesize=256），R_spheroid表示地球赤道半径（R_spheroid=6378137米）。

通过上述计算公式可以计算出金字塔第一级的空间分辨率为157286.4米，第二级的空间分辨率为78643.2米，以此类推。GGPM的缩放级别理论上是可以不断扩展的，但在实际应用中还需要考虑已有遥感影像数据源的最高空间分辨率，目前全球在轨运行的遥感卫星获取的遥感数据的最高分辨率为0.6米，那么GGPM的最大缩放级别为19。表1列出了实施例整数缩放级别、空间分辨率和数据源之间的关系。具体实施时，该关系和本领域技术人员所选用的数据源有关。

表1

卫星遥感影像数据通常为GeoTIFF文件格式，而浏览器只支持JPEG、GIF、BMP和PNG这四种图像格式（网络地图服务是通过Web浏览器进行访问），其中JPEG适合色彩丰富的遥感影像，此外JPEG所需的存储空间也比较小，而存储空间对于海量遥感数据是非常关键的，例如对于一张具有30万种颜色和2088×1128像素大小的图像，JPEG格式的文件大小为147Kb，而PNG格式的大小达到2.67Mb，因而可将GeoTIFF格式转换为JPEG格式。

为了减少数据的存储空间，可选择合适的JPEG图像压缩率。压缩率越高，文件所需的存储空间就越小，但图像质量就会降低，通过实验分析发现图像质量从90到80时，存储空间有一个较大的减少，图像质量从80再往下降低时，存储空间没有明显的下降，当图像质量设置过低时图像就会非常模糊，建议具体实施时高分辨率的遥感影像使用较高的图像设置，而较低分辨率的遥感影像使用较低的图像设置。

3、Tiling遥感影像并编码

Tiling遥感影像表示切分遥感影像，即将图像切割分成瓦片。本步骤根据步骤1中所建立的全球地理空间金字塔模型相应全球四叉树结构和Tile文件编码体系实现。实施例的遥感数据Tiling和编码规则为：

（1）根据遥感数据的SMP坐标范围将它叠加到虚拟切割网格上，然后沿网格线将遥感数据切割成256×256大小的Tile；

（2）根据最小包围矩形的左上角XY坐标和右下角XY坐标对有效Tile进行编码，将XY网格坐标和缩放级别作为Tile的文件名。可以根据本发明实施例设计的编码体系对切割好的Tile进行编码，编码通过文件名来表示。

Tiling遥感影像原理为：将全球四叉树以虚拟瓦片为网格划分的某一层作为虚拟切割网格，以金字塔第16层（即第16级）为例，其切割网格X轴取值范围从0到2^16-1-1=32767，Y轴的取值范围同样从0到32767。切割时将图像按照SMP投影坐标对应的原则叠加到该级别虚拟切割网的相应位置。由于遥感影像图幅（景）的尺寸往往不是切割网格的整数倍，即不能够与切割网格线完全对齐，因此本发明提出遥感影像图幅切割网格最小包围矩形概念，它表示影像图幅在切割网格（图中是由虚线构成的）中所覆盖的以网格为单位的最小外接矩形，如图8所示，（X_start,Y_start）为最小包围矩形左上点的坐标,（X_end,Y_end）为最小包围矩形右下点的坐标。切割网格最小包围矩形内的Tile分为有效Tile（有效瓦片）和无效Tile（无效瓦片）。有效Tile是指在最小包围矩形内被遥感影像完全覆盖的Tile，无效Tile是指不能被遥感影像完全覆盖的Tile。

例如有一幅覆盖北京地区的遥感影像经过投影转换（UTM→SMP）和分辨率调整后的坐标范围为(12934501.384米,4860561.694米,12982894.984米,4814251.294米)，空间分辨率为4.8米，对应于金字塔的第16级，如图10所示。图10中有两个原点。一个是网格坐标的原点，即屏幕显示的原点，一个是Mercator投影的原点，即投影平面的中心，以投影平面的中心建立坐标系可得象限1、2、3、4。金字塔16级的Tile数量为1073741824，每个Tile的跨度为256×4.8=1228.8米，图中北京地区遥感影像的左上角的坐标(+12934501.384米,+4860561.694米)，右下角的坐标为(+12982894.984米,+4814251.294米)。根据北京地区遥感影像的左上角的坐标(12934501.384米,4860561.694米)可以计算出它在SMP坐标系统第一象限内X方向上跨越的Tile数TileNumX=Math.floor(12934501.384/1228.8)=10526，Y方向上跨越的Tile数TileNumY=Math.floor(4860561.694/1228.8)=3955，Math.floor()函数表示向下取整。因而，最小包围矩形左上角的网格坐标为Xstart=214+10526-1=26909和Ystart=214-3955-1=12428，同理右下角的网格坐标为Xend=26949和Yend=12466。通过最小包围矩形左上角和右下角的网格坐标就可以确定遥感影像所有有效Tile的网格坐标。

在确定每个Tile的网格坐标后就可以对遥感影像进行切分并为每个生成的256×256大小的Tile进行文件命名，按照本发明实施例设计的编码体系，形式为X_rowY_colZ_level.jpeg，其中row和col为网格坐标，level为缩放级别。

4、根据数据分布策略对Tile进行分组

实施例的数据分布策略为：

（1）本发明将Tile图片根据Hilbert曲线Tile分布方法（HCTD）分布到不同的Web服务器上。该方法首先从GQ的第2层开始采用Hilbert空间填充曲线进行网格填充，如图10所示，第二层采用一阶Hilbert曲线进行填充，第三层采用二阶Hilbert曲线填充，依次类推…第19层采用18阶Hilbert曲线填充。

（2）计算出每一层每一个网格的Hilbert编码，然后使用HCTD方法进行Tile的分布。

将遥感影像切分为Tile后，需要根据数据分布策略对Tile进行分组。本发明提出的瓦片分布算法称为数据分布策略HCTD，是一种基于Hilbert空间填充曲线的Tile数据分布模式，它以Tile网格矩阵的左上角Tile的中心点作为起点开始填充每一个Tile并且只填充一次，这样生成的曲线赋予了每个Tile一个唯一的Hilbert编码，表示为h(i,j),其中i和j为Tile网格坐标，i和j∈[0,2^N-1]，其中N代表缩放级别；然后根据一个映射函数f(i,j)=h(i,j)mod M将每个Tile分配到不同的服务器磁盘上，其中M为服务器的数量。

求解Hilbert编码是Tile分布的重要基础，经典的Hilbert码生成算法由Faloutsos和Roseman提出并由Warren M.Lam和Jerome M.Shapiro进行改进，由于这种算法是基于空间目标点所在的栅格格网行列数的二进制位操作，故其算法的复杂度为，其中n_i为与空间目标点i所在最终栅格格网行列数中较大者所对应的二进制位数。本发明可以利用经典的Hilbert码生成算法实施，但这种Hilbert编码算法的时间复杂度会随着曲线阶数的增加而呈几何级数增加，为此本发明提出一种基于象限空间划分的高效Hilbert编码算法，能够进一步有效降低时间复杂度。

实施例根据Hilbert曲线的重要性质设计了三个快查表用于Hilbert编码算法。表2表示子曲线的方位是由它的父曲线的方位和子曲线所处的象限确定的，表格中第一列的数字表示基本曲线的Orientation（方位，也称之为类型），四个象限的位置参见图11，第1个象限处于左上，记为象限I；第2个象限处于右上，记为象限II，第3个象限处于左下，记为象限III，第4个象限处于右下，记为象限IV。表3表示Hilbert编码是由父曲线的方位和子曲线所处象限确定的。表4表示当Hilbert曲线不能再被划分时，Hilbert编码是由基本曲线的方位和Tile所处的象限确定的。

表2子曲线方位快查表

表3象限划分过程中Hilbert编码快查表

表4Hilbert编码尾码快查表

表中，Hcoding为Hilbert(希尔伯特)编码，+=表示相加后再赋值，<<表示移位操作，n为缩放级别。各方位1、2、3、4参见图12，其中图12a为方位1（类型1）的示意图，按左上、左下、右下、右上的顺序遍历2×2的网格；图12b为方位2（类型2）的示意图，按左上、右上、右下、左下的顺序遍历2×2的网格；图12c为方位3（类型3）的示意图，按右下、左下、左上、右上的顺序遍历2×2的网格；图12d为方位4（类型4）的示意图，按右下、右上、左上、左下的顺序遍历2×2的网格。

实施例的基于象限空间划分的Hilbert编码过程，采用对某层中任一目标Tile的算法描述如下：

(1)对该层中相应的Hilbert曲线，按照大致方向（即曲线整体走向），如图15确定一个给定的Hilbert曲线的方位，作为父曲线的方位。计算目标Tile属于哪个象限（目标象限），包括根据目标瓦片的Z坐标值(缩放级别)确定这一级的瓦片矩阵，对这个瓦片矩阵进行象限划分，就可以得到每个象限的最小XY坐标值和最大XY坐标值，然后将目标瓦片的XY值与它们进行比较，落在哪个范围内就表示在哪个象限。然后所得象限和父曲线的方位，通过查询表2获得位于目标象限的子曲线的方位。

(2)判断是否能够将目标象限划分为4个子象限，是则进行划分并进入(3)。如果目标象限不能够再被进一步划分，根据Tile所处象限和父曲线的方位，通过查询表4执行Hilbert编码的最后处理，然后返回处理结果，结束流程。

(3)确定目标Tile属于哪个目标子象限。根据父象限的曲线的方位和目标子象限的位置，确定目标子象限内的曲线的方位。之后，根据子曲线所处象限和父曲线的方位，通过查询表3得到目标子象限曲线的初始顺序值。如图13中是一个三阶Hilbert填充曲线（H3），三阶Hilbert填充曲线经过的网格依次编码为0、1、2…63，将它划分为4个象限，每个象限的子曲线的初始顺序值即进入该象限的第一个格子的编码，例如象限III子曲线初始顺序值为16。

(4)跳转到(2)，以(3)确定的目标子象限作为当前目标象限，进一步划分子象限。

5、生成每个分组上Tile文件的分层目录结构

设每个服务器上的所有Tile文件属于一个分组，具体实施时，本领域技术人员可根据数据来源和缩放级数等自行设计Tile分层目录结构。一般预定义表示不同数据图层的一级目录，每个一级目录下定义n个二级目录，每一个二级目录下的目录层次结构可根据遥感影像的特征设置。合理的Tile分层目录存储结构设计能够有效支持遥感数据增量式更新和目标Tile快速检索。

本发明实施例设计的一个支持全球多维遥感数据的Tile分层目录结构，如图7所示，构建方法如下：

（1）预定义5个表示不同数据类型（图层）的一级目录，其中BaseMap表示多光谱合成彩色遥感影像数据，PanMap表示全色波段遥感影像数据，SARMap表示合成孔径雷达遥感数据，ReliefMap表示数字高程遥感数据（NASA SRTM），VectorMap表示光栅化的矢量图层数据，不同数据类型目录下的目录层次结构基本相同；

（2）定义19个二级目录，这19个二级目录代表不同的缩放级别，19级遥感影像的空间分辨率为0.6米，每一个缩放级别都对应整个地球平面，并且是由大小固定的Tile矩阵构成，每一级覆盖全球的Tile数量为4^z-1(z为缩放级别)；以BaseMap为例，在BaseMap下定义了19个二级目录，这19个二级目录代表不同的缩放级别，19级遥感影像的空间分辨率为0.6米；每一个缩放级别都对应整个地球平面，并且是由大小固定的Tile矩阵构成；每一级覆盖整个地球的Tile数量为4^z-1(z为缩放级别）；

（3）每一个缩放级别目录下的目录层次结构根据遥感影像数据的特征而有所不同，具体的目录结构见(4)～(6);

（4）对于1到7级的遥感影像，由于它们的空间分辨率较低，主要用于导航，因此将所有的Tile图片直接存储在缩放级别目录下，如图7的1到7级中某级z₁下直接存储Tile图片，例如采用X_rowY_colZ_level.ext格式的某Tile图片X_iY_jZ_k.jpg;

（5）对于8到13级的遥感影像，由于空间分辨率比较高并且Tile图片数量比较多，例如第13级的Tile数量达到16777216张，同时为了支持历史数据查询，设计了第3级目录和第4级目录，其中第3级目录用于地理空间的划分，这样使得每个路径下的目录或文件个数不至于太多，从而影响文件系统的性能，划分的原则是使得最后一级目录下的图片数量不超过4096个，这样第3级目录的数量等于4^z-7，如图7的8到13级中某级z₂下有目录1…e…4^{（z -7）}，第3级目录名称可根据网格坐标命名为XiYj，0≤i≤2^z-7-1，0≤j≤2^z-7-1，z表示缩放级别，此处用网格坐标（i，j）表示更大的地理范围，i，j为网格坐标X坐标值、Y坐标值，即第3级的一个目录包含4096个瓦片，即64×64瓦片矩阵；第4级目录用于辨识数据获取的时间，如图7中第3级目录中某目录i下有第4级目录Latest、Date1…DateN，Latest表示最新时间（即这个目录i下保存最新获取的遥感数据）第4级目录名称命名为YYYYMM，用于存储历史数据，其中YYYY表示年，MM表示月；

（6）对于14到19级的遥感影像数据，由于Tile图片数量巨大，例如14级的Tile数量达到67108864张，为了保证最后一级目录的文件数量不超过4096个，需要对地理空间进行两级划分，第一级地理空间划分的目录数量为4⁶，即14到16级的目录数量固定为4⁶，如图7的14到19级中某级z₃下有目录1…s…4⁶，第3级目录名称命名为XiYj，0≤i≤2⁶-1，0≤j≤2⁶-1，此处用网格坐标表示更大的地理范围；第二级地理空间划分的目录数为4^(z-13)，如图7的第3级目录1…i…4⁶之一下有目录1…t…4^(z-13)，第4级目录名称命名为XiYj，0≤i≤2^(z-13)-1，0≤j≤2^(z-13)-1，此处用网格坐标表示更大的地理范围；第5级目录用于表示数据获取的时间，具体表示参见（5）中的第4级目录。

实施例创建目录结构的顺序如下：

(a)根据Tile数据类型创建第1级目录(数据图层)；

(b)根据空间分辨率创建第2级目录(缩放级别)；

(c)对于8到13级的遥感影像根据Tile的文件名(网格坐标)创建第2级目录，对于14到19级的遥感影像根据Tile的文件名(网格坐标)创建第3级和第4级目录；

(d)对于8到13级的遥感影像根据获取的时间创建第4级目录，对于14到19级的遥感影像根据获取的时间创建第5级目录。

对于每个分组的Tile可根据本发明实施例设计的以上Tile分层目录结构生成相应的目录结构，对于缩放级别大于7的遥感影像数据的每一次更新都会创建一个新的时间目录并将Tile数据存储到这个目录下，同时覆盖Latest目录下的Tile数据，这样既保证了Latest目录下的数据总是最新的，同时也使得历史数据不会被覆盖。具体实施时，本领域技术人员可采用现有集群服务器技术实现服务器端，在服务器端包括管理节点和多个Web服务器，管理节点执行步骤1～5，负责构建整个模型以及处理数据，并将数据分发到不同的Web服务器上。

6、将Tile数据更新到相应的Web服务器上，可以根据步骤5所得分层目录结构将各组的Tile文件更新到相应的服务器上。

7、当用户通过客户端查看地图时，计算出地图可视范围内需要下载的Tile图片相应Tile文件文件名。客户端包括智能手机、计算机等，在用户使用的客户端请求下载地图进行显示时，可由客户端根据用户与服务端的交互来计算所需的Tile的文件名，并根据Tile分布策略计算出这些Tile所在的Web服务器的存储地址。具体实施时，可由本领域技术人员采用HTML和JavaScript设计相应客户端程序，使用时运行在浏览器中。

一般只需要下载新增Tile，客户端已经有的地图数据无需重复下载。实施例计算新增Tile的方法采用使用现有的HTML和CSS技术为：参见图14，根据InnerDiv相对OuterDiv的位置（Left，Top）来计算需要显示的Tile的网格坐标，InnerDiv表示虚拟网格区域，可以容纳相应缩放级别的所有Tile，OuterDiv表示窗口可视区域，即客户端提供的地图可视范围。其中，Left表示InnerDiv相对OuterDiv的位置在X轴上的差距，Top表示InnerDiv相对OuterDiv的位置在Y轴上的差距。

8、根据步骤7得到的所需的Tile的文件名，从相应服务器获取这些Tile。根据本发明的Tile分层目录存储结构能够实现目标Tile快速检索，具体实施时可以并行地从多个服务器获取这些Tile。

9、在客户端根据步骤8获取的Tile图片进行显示。通常是在客户端将步骤8所得新增Tile和原有Tile无缝拼接成一张完整地图，一般需要将来自不同服务器的Tile图片无缝拼接成完整的地图。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种轻量级全球多维遥感影像网络地图服务实现方法，其特征在于，包括以下步骤：

步骤1，通过自顶向下逐级划分全球地理空间范围构建全球地理空间金字塔模型，采用全球四叉树结构表示全球地理空间金字塔模型，并建立相应Tile文件编码体系，其中Tile表示瓦片；

步骤2，对遥感影像进行预处理，包括投影变换为一致的投影方式和格式转换为一致的图像格式；

所述对遥感影像进行预处理包括以下步骤，

(1)将遥感数据进行投影变换，将UTM转换为SMP，投影变换后根据遥感数据的空间分辨率查询关系表来确定它所在的缩放级别，如果空间分辨率正好等于表中的某一个空间分辨率则直接确定它的缩放级别，如果介于两个分辨率之间，则需要对遥感数据进行分辨率调整；所述关系表包括整数缩放级别、空间分辨率和数据源之间的关系；

(2)对遥感数据进行格式转换，将GeoTiff格式转换为JPEG格式；

步骤3，根据全球四叉树结构对遥感影像进行切分得到Tile图片，并根据Tile文件编码体系对每个Tile图片相应存储数据的Tile文件进行编码获得文件名；

所述Tile文件编码体系为，每个Tile图片相应存储数据的Tile文件根据Tile图片的全球四叉树网格坐标<level，<row,col>>生成文件名，其中level表示Tile图片在全球四叉树结构中的层数，row表示Tile的X坐标值，col表示Tile的Y坐标值；

对每一层遥感数据分别处理方式如下，

(1)根据遥感数据的球形墨卡托投影坐标范围叠加到全球四叉树结构中该层相应虚拟切割网格上，然后沿网格线将遥感数据切割成Tile图片；

所述球形墨卡托投影坐标范围是以距离米为单位的坐标体系而非经纬度坐标体系，将全球四叉树以虚拟瓦片为网格划分的某一层作为虚拟切割网格，虚拟切割网格表示在SMP坐标体系下的距离范围，而不是经纬度范围，切割时将图像按照SMP投影坐标对应的原则叠加到该级别虚拟切割网的相应位置；

(2)根据遥感数据的最小包围矩形的左上角XY坐标和右下角XY坐标，按照Tile文件编码体系对有效Tile图片进行编码，所述最小包围矩形是遥感数据的影像图幅在虚拟切割网格中所覆盖的以网格为单位的最小外接矩形；

切割网格最小包围矩形内的Tile分为有效Tile和无效Tile，有效Tile是指在最小包围矩形内被遥感影像完全覆盖的Tile，无效Tile是指不能被遥感影像完全覆盖的Tile；

步骤4，根据数据分布策略对Tile文件进行分组，包括确定每个Tile文件存放的Web服务器，存储在相同服务器上的Tile文件为一组；所述根据数据分布策略对Tile文件进行分组的实现方式如下，

(1)从全球四叉树结构的第二层开始采用Hilbert空间填充曲线进行网格填充，金字塔第二层采用一阶Hilbert曲线进行填充，第三层采用二阶Hilbert曲线进行填充，依次类推…第n层采用n-1阶Hilbert曲线进行填充；

(2)计算出每一层每一个网格的Hilbert编码，对每个Tile文件根据Hilbert编码映射到一个Web服务器；

步骤5，生成每个分组上所有Tile文件的分层目录结构；

步骤6，根据步骤5所得分层目录结构将各组的Tile文件更新到相应的服务器上；

步骤7，当用户通过客户端查看地图时，计算出地图可视范围内需要下载的Tile图片相应Tile文件文件名，并根据步骤4的数据分布策略计算出每个相应Tile文件实际存储的Web服务器地址；所述计算地图可视范围内需要下载的Tile图片相应Tile文件文件名时，根据InnerDiv相对OuterDiv的位置来计算需要下载的Tile图片的网格坐标，InnerDiv表示容纳相应层数的所有Tile图片的虚拟网格区域，OuterDiv表示地图可视范围；

步骤8，根据步骤7所得Web服务器地址，按照文件名从相应Web服务器下载相应Tile文件到客户端，包括基于Tile分层目录存储结构实现目标Tile快速检索，支持并行地从多个服务器获取Tile；

步骤9，在客户端根据步骤8获取的Tile图片进行显示，包括将来自不同服务器的Tile图片无缝拼接成完整的地图。

2.根据权利要求1所述的轻量级全球多维遥感影像网络地图服务实现方法，其特征在于：全球地理空间金字塔模型每一层的遥感影像来自相同或不同的数据源。

3.根据权利要求2所述的轻量级全球多维遥感影像网络地图服务实现方法，其特征在于：所述一致的投影方式采用球形墨卡托投影，所述一致的图像格式采用JPEG格式。

4.根据权利要求3所述的轻量级全球多维遥感影像网络地图服务实现方法，其特征在于：所述全球四叉树结构，四叉树中的每个节点表示一个虚拟瓦片，设全球四叉树结构有n层，顶层用一个虚拟瓦片表示全球地理空间范围，第2层用4个虚拟瓦片表示全球地理空间范围，第3层用16个虚拟瓦片表示全球地理空间范围，依次类推…第n层由4的n-1次方个虚拟瓦片拼接而成。

5.根据权利要求1所述的轻量级全球多维遥感影像网络地图服务实现方法，其特征在于：步骤5中，生成的每个分组上所有Tile文件的分层目录结构如下，

预定义表示不同数据图层的一级目录，每个一级目录下定义n个二级目录，每一个二级目录下的目录层次结构根据遥感影像的特征设置。

6.根据权利要求1所述的轻量级全球多维遥感影像网络地图服务实现方法，其特征在于：计算出每一层每一个网格的Hilbert编码的实现方式如下，

首先建立以下快查表，

表2子曲线方位快查表

表3象限划分过程中Hilbert编码快查表

表4 Hilbert编码尾码快查表

表中，Hcoding为Hilbert编码，+＝表示相加后再赋值，<<表示移位操作，n为缩放级别；然后执行步骤如下，

(1)对当前层中相应的Hilbert曲线，确定父曲线的方位，计算目标Tile图片所属目标象限，然后根据所得目标象限和父曲线的方位，通过查询表2获得位于目标象限的子曲线的方位；

(2)判断是否能够将目标象限划分为4个子象限，是则进行划分并进入步骤(3)，否则根据Tile图片所处目标象限和父曲线的方位，通过查询表4执行Hilbert编码的最后处理，然后返回处理结果，结束流程；

(3)确定目标Tile图片属于哪个目标子象限，根据父象限的曲线的方位和目标子象限的位置，确定目标子象限内的曲线的方位；之后，根据子曲线所处象限和父曲线的方位，通过查询表3得到目标子象限曲线的初始顺序值；

(4)跳转到步骤(2)，以步骤(3)确定的目标子象限作为当前目标象限，进一步划分子象限。