EN中文

地球信息科学学报 >

2021 , Vol. 23 >Issue 10: 1756 - 1766

DOI: https://doi.org/10.12082/dqxxkx.2021.210029

地球信息科学理论与方法

一种与地图服务结合的栅格瓦片计算模型

  • 胡毅荣 , 1, 2 ,
  • 王超 1, 2 ,
  • 杜震洪 , 1, 2, * ,
  • 张丰 1, 2 ,
  • 刘仁义 1, 2
展开
  • 1.浙江大学地球科学学院,地理信息科学研究所,杭州 310027
  • 2.浙江省资源与环境信息系统重点实验室,杭州 310028
* 杜震洪（1981— ）,男,浙江衢州人,博士,教授,主要从事遥感与地理信息系统,时空大数据与人工智能研究。E-mail:

胡毅荣（1994— ）,男,浙江嘉兴人,硕士生,主要从事地理大数据存储与计算研究。E-mail:

收稿日期: 2021-01-18

  要求修回日期: 2021-03-13

  网络出版日期: 2021-12-25

基金资助

国家重点研发计划项目(2018YFB0505000)

国家自然科学基金项目(41922043)

国家自然科学基金项目(41871287)

国家自然科学基金项目(42001323)

版权

版权所有,未经授权,不得转载、摘编本刊文章,不得使用本刊的版式设计。
收起

A Raster Tile Calculation Model Combined with Map Service

  • HU Yirong , 1, 2 ,
  • WANG Chao 1, 2 ,
  • DU Zhenhong , 1, 2, * ,
  • ZHANG Feng 1, 2 ,
  • LIU Renyi 1, 2
Expand
  • 1. Department of Geographic Information Science, School of Earth Sciences, Zhejiang University, Hangzhou 310027, China
  • 2. Zhejiang Provincial Key Lab of GIS, Hangzhou 310028, China
* DU Zhenhong, E-mail:

Received date: 2021-01-18

  Request revised date: 2021-03-13

  Online published: 2021-12-25

Supported by

National Key Research and Development Program of China(2018YFB0505000)

National Natural Science Foundation of China(41922043)

National Natural Science Foundation of China(41871287)

National Natural Science Foundation of China(42001323)

Copyright

Copyright reserved © 2021
Fold

摘要

随着遥感影像数据的快速增长,对于栅格数据高效的信息处理和价值挖掘方式提出了更大的挑战,传统地图服务聚焦于内容的共享与可视化,缺乏对影像实时分析处理功能。本研究以地图服务的形式实现了对栅格瓦片数据实时分析处理能力,将云优化的GeoTIFF（Cloud Optimized GeoTIFF,COG）作为数据组织方式,设计了分布式协同预取策略,实现了栅格瓦片数据的冷热加载,优化了从云端读取影像数据的效率。在栅格瓦片数据高效加载的基础下,提出了一种基于表达式的栅格瓦片处理模型,通过对表达式转换建模为计算工作流,在地图服务的请求中实现对栅格瓦片的实时处理,对存储在云端的海量遥感数据进行快速分析,实现原始数据到数据产品的直接可视化转换。针对全量数据参与的场景,使用合适的重采样数据进行简化计算,以满足地图服务的实时性。使用了NDVI、地物分类、植被覆盖度三类不同复杂度模型,在地图服务中对Landsat 8影像进行了实时计算分析。实验结果表明,该处理模型能对栅格瓦片进行有效分析,且能进行分布式扩展,在高并发场景下能够提供稳定的地图服务能力,适应各层级尺度的计算,对未来地图服务的发展提供了一种新思路。

关键词： 栅格瓦片; 地图服务; Web GIS; 遥感处理; 工作流; 表达式计算; 瓦片预取; 分布式

本文引用格式

胡毅荣 , 王超 , 杜震洪 , 张丰 , 刘仁义 . 一种与地图服务结合的栅格瓦片计算模型[J]. 地球信息科学学报, 2021 , 23(10) : 1756 -1766 . DOI: 10.12082/dqxxkx.2021.210029

Abstract

With the rapid growth of remote sensing data, greater challenges arise in raster data efficient processing and value mining. Traditional map services focus on content sharing and visualization, but lacking real-time image analysis and processing functions. In this study, the real-time analysis and processing capabilities of raster tile data are realized in the form of map service. The cloud optimized GeoTIFF (Cloud Optimized GeoTIFF, COG) is used as the data organization method. The distributed collaborative prefetching strategy is designed to realize the raster tile loading in a cold or hot way, which optimizes the efficiency of reading image data from the cloud. Based on the efficient raster tile data loading, an expression-based raster tile processing model is proposed. By converting the expression into a calculation workflow, the raster tile is processed in the request of the map service in real time. The massive remote sensing data stored in the cloud is quickly analyzed to realize the direct visual conversion from raw data to products. For scenarios where full data are involved, use appropriate resampling data to simplify calculations to meet the real-time performance of map services. Three types of different complexity models, NDVI, ground object classification, and fractional vegetation cover, are used to perform real-time calculation and analysis on Landsat 8 images in the map service. Experimental results show that the processing model can effectively analyze raster tiles, and can be extended in a distributed manner. It can provide stable map service capabilities in high-concurrency scenarios, adapt to calculations at various levels and scales, and contribute a new idea to the future development of map service.

Key words： raster tiles; map service; Web GIS; remote sensing processing; workflow; expression calculation; tile prefetching; distributed

1 引言

在地理信息系统（GIS）中面向服务的体系架构使用变得越来越广泛[1]。对公开数据集实现标准的WEB接口有助于隐藏数据集的技术细节,从而潜在地让更广泛的人员可以使用这些数据集。开放地理空间联盟（OGC）发布了一套用于互联网数据交换的开放标准[2,3]。其中Web地图服务（Web Map Service, WMS）规范描述了客户端（如Web浏览器或桌面应用程序）如何从服务器请求自定义生成的、地理定位的数据集图像,以及相关的元数据[4]。OGC地图服务标准聚焦于数据内容的共享与可视化,使得来自不同网络地图服务的图像可以精确地叠加在地理信息系统中,允许可视化和相互比较[5],但是未对地图服务内容的计算处理进行规范。
目前遥感传感器生成数据的能力远远超过了处理这些数据的能力,无论是“人工”解译图像还是对图像进行数字分析,流程上较复杂,且具有一定的滞后性。然而,现实生活中存在许多对时间要求严格的应用场景,例如应急响应、救援救灾和基于位置的服务等,都需要即时访问各种遥感影像数据或者影像数据的衍生产品信息,以快速做出决策和及时采取行动,这些应用需求进一步对综合型地图服务提出了更加高效、实时和动态的要求[6]。地理空间处理常用的接口标准是OGC Web Processing Service（WPS）,通过XML的形式对已部署服务的输入和输出规则以及服务的请求方法进行标准化,满足了很多对于地理数据的处理操作需求,但不适用于处理实时工作流的任务场景[7]。因此,如何以地图服务为载体,对数据内容进行便捷、高效、实时分析,并对结果进行发布与共享,对影像数据的价值提取具有重要意义,是当前遥感信息化领域的重点与难点。
遥感影像由于数据量大,在地图服务中通常表现为栅格瓦片的形式。按照影像栅格瓦片的生成原理,本文梳理了当前较常见的3种栅格瓦片地图服务形式。① 静态瓦片地图服务,即将地理空间范围内的影像数据预先生成静态瓦片图像[8],然后将这些瓦片图像作为地图服务,这种方式降低了地图服务的灵活性,减少服务器负载的同时提高了可伸缩性[9],使得该服务形式广泛用于互联网底图,例如谷歌地图和必应地图等。王浩等[10]基于静态瓦片图像研究了瓦片缓存置换算法;Ricardo等[11]提出了使用神经网络的方法实现对瓦片地图的智能预取,这些研究都进一步提高了静态瓦片地图的服务能力。② 通过GIS服务器（典型的有GeoServer、ArcGIS Server等）关联栅格数据源,并发布成地图服 务[12],该类方式能够较快提供地图服务可视化的功能,且遵循一系列OGC规范[13],但是面对管理和发布海量影像数据时,其横向扩展的能力具有一定的局限性。③ 实时地图服务,将影像数据以Cloud Optimized GeoTIFF（COG）格式存储在云端,能够获取较小的兴趣区数据并进行渲染,避免了下载和提取单个波段到本地磁盘的需要,作为一种新兴的可视化方式,具备基于像素的分析能力[14]。Steve Kopp[15]提出了实时图像交互式数据探索分析的方法,其核心是通过Web应用程序探索多维数据集并可视化结果。Hnatushenko 基于COG开发一个农业监测的Web地图应用程序,将地图瓦片服务与用于分布式计算的地理空间处理服务结合在一起,帮助种植者优化土壤处理和产量,并提供作物歉收或饥荒的早期预警[16]。表1对3种地图服务瓦片的特点进行了总结。
表1 3种地图服务瓦片特点

Tab. 1 Three map service tiles feature

地图服务瓦片类型 瓦片生成方式 数据读取方式 具备分析性
静态瓦片 预生成静态图片 直接读取静态图片,性能高 否
GIS服务器生成的瓦片 实时读取本地影像 条带状读取,性能低 是
基于COG生成的瓦片 实时读取云端影像 瓦片式读取,性能高 是
上述3种方式中,后面2个都是实时从影像数据中获取特定空间范围内的栅格瓦片,而当前主流的地图服务只考虑了将栅格瓦片通过web服务呈现给客户端,未将栅格瓦片的计算与地图服务相结合。本文结合COG数据组织形式,提出了一种与地图服务结合的栅格瓦片计算模型,使得地图服务具备计算分析的能力,从而实现对遥感数据信息的快速提取。该计算模型丰富了地图服务的内涵,为地图服务的发展提供了新的方向。

2 研究方法

2.1 Cloud Optimized GeoTIFF数据组织

GeoTIFF是建立在TIFF ISO标准之上的OGC标准,存储为二维栅格数据层,大量与遥感、地理信息系统、制图和测绘相关公司和组织以GeoTIFF作为地理栅格数据的存储标准[17]。Cloud Optimized GeoTIFF（COG）是支持云上更高效读取数据的组织形式,同时也是一个托管在HTTP文件服务器上的标准GeoTIFF文件,该文件生成方式依赖于开源的地理空间数据转换库GDAL[18],常见的GIS软件平台开始对COG格式的数据逐步进行了支持,如QGIS、Google Earth Engine等[19]。其内部数据组织主要依赖2个特性[20]：① 切片（Tile）,将条带式数据的组织方式更改为瓦片式,图像被分组成有规律大小的连续块（默认是512像素×512像素）,这些连续块有规律地存储在文件中;② 对原始数据经过重采样后生成一系列概览层（Overview）,以适应不同层级尺度获取数据的需要。COG影像组织方式如 图1所示,构成金字塔结构,该组织方式适合在网络环境中传输,结合HTTP GET范围请求可对COG内部的数据进行部分获取。通过使用GDAL Virtual Format机制可实现对影像元数据、缩略图、指定范围的空间子集数据（单/多光谱波段）的获取[21],当访问一定空间范围但不同分辨率尺寸的数据瓦片时,该机制会选择适宜的Overview层级数据进行读取,本质和读取瓦片金字塔模型类似,使得服务器读取数据的压力变小,提高了系统响应速度[22]。
图1 COG组织方式

Fig. 1 Organization of COG

2.2 分布式栅格瓦片预取策略设计与实现

在地图瓦片服务中,每个请求都由三元组（缩放层级、行号、列号）坐标表示,以此来获取指定空间区域的数据。若每个请求都从云存储中实时获取数据,则大量时间将消耗在数据获取中,降低地图服务的效率。地图瓦片服务的请求呈现规律性[23],即单张瓦片被请求时,其邻域瓦片将来也有可能会被访问到。依据该规律,本文设计并实现了一种分布式环境下服务结点间栅格瓦片的协同预取策略。一方面,实现冷热数据的混合读取,保证瓦片数据的高效加载;另一方面,当多用户同时或者经常访问相同的瓦片数据时,热数据可为服务器端和客户端之间提供了高速的数据缓存,势必会大大减少服务端不必要的处理开销,降低客户端的等待时间,提高用户体验[24]。
如图2所示,当请求瓦片T1和T2时,按照负载均衡策略,由不同服务结点进行处理。服务结点除了获取当前栅格瓦片之外,还会计算当前请求的邻域瓦片的三元组坐标,并将当前请求的影像和波段等参数封装成消息（T3和T4）发送到消息队列中,分布式消息投递方式使得消息中间件能够快速汇聚来自不同服务结点投递的瓦片预取信息。使用消息队列是为了在不阻塞当前的服务请求的同时可以异步获取邻域的栅格瓦片,监听消息队列的服务进程收到消息后,根据消息内容执行从云存储中获取栅格瓦片的操作,将得到的瓦片经过序列化压缩后存储到集中式缓存Redis中,并为其初始化生存时间（Time to Live,TTL）。消息队列的任务执行者,由多个消费者结点承担,每个消费者从云存储中读取栅格瓦片数据并缓存到Redis,分布式任务执行使得短时间内将大量云端冷数据瓦片转化为Redis中的热数据。当缓存中的瓦片数据被再次请求时,根据当前存活的时间梯度更新生存时间：
f t = T 1 + e - 0.5 t < 0.5 T T 1 + e - 0.5 - l n T t - 1 0.5 T ≤ t
图2 分布式栅格瓦片预取流程

Fig. 2 Distributed raster tile prefetching process

式中： t为当前数据瓦片的生存时间, T为设置的瓦片最大生存时间, f t为更新后的生存时间。当 t<0.5 T时,统一设置为固定生存时间;当0.5 T≤t时,对瓦片生存时间使用sigmoid函数进行非线性更新,使得TTL值低的瓦片得到更多的时间增量,而TTL接近阈值的瓦片则增加更少的存活时间。当缓存中瓦片的TTL降低到0时,则会释放掉缓存中的数据。

2.3 栅格瓦片计算模型

为了实现对瓦片数据的处理分析,本文提出一种栅格瓦片计算模型,包括输入输出的标准和瓦片数据的处理方式2个方面内容。首先,对计算模型的输入输出进行了规范定义,输入是在地图服务的请求中增加运算表达式参数,而在栅格瓦片的渲染中呈现计算输出结果,达到所见即所得的效果。接着,提出了一种以工作流的方式对运算表达式进行迭代计算的方法,在工作流执行的过程中,调用基础函数和模型函数进行运算。
栅格瓦片的运算可以通过表达式表示,如 ArcGIS和ENVI等软件都具备使用表达式对波段进行运算的功能。本文基于标准的栅格瓦片地图服务,在请求中增加一个运算表达式参数来传递栅格瓦片的计算需求,使得地图服务具备实时计算的能力。运算表达式中包含的运算符和函数表示对数据处理的流程,其中波段数据表示的方式由字母B和数字组成,例如：B1代表该影像数据的第一个波段数据,其他波段数据以此方法类推。
工作流在应用程序中是一系列数据操作和计算的步骤,可以清晰反映步骤之间的流程[25]。根据表达式中包含的运算优先级将其建模成为一个工作流。图3是1个运算表达式转换为2个计算结点的工作流的示意图,按照预设的函数自动传递参数、数据或执行任务。该工作流体现了计算结点间的执行步骤和依赖关系,构成了有向无环图。
图3 栅格瓦片计算工作流

Fig. 3 Raster tile calculation workflow

完成对计算工作流的建模,首先需要将表达式解析成栅格瓦片数据元信息和计算结点,栅格瓦片数据元信息是一个映射表,包含了获取该栅格瓦片的数据存储路径、空间三元组坐标、波段等信息。计算结点具有函数名和计算参数属性,负责执行完整的工作流程,包括参数替换、数据获取、函数计算和递归等。计算结点还拥有触发计算的方法compute。若要触发整个计算工作流的运行,调用最后一个计算结点的compute方法即可,工作流将从后往前递归地执行结点运算或瓦片数据获取。结合图3说明计算工作流的具体执行步骤：
（1）获取最后一个计算结点Node3,执行compute方法,循环迭代计算参数Node1与Node2,判断计算参数类型为结点（Node）类型,需要进行参数替换,执行步骤（2）。
（2）执行Node1的compute方法,循环迭代计算参数MetaData1和MetaData2,判断计算参数类型为数据元信息（MetaData）类型,则根据数据元信息获取栅格瓦片,返回结果赋予当前结点的函数参数,完成Node1函数的参数替换,执行步骤（3）。
（3）此时Node1的函数参数为栅格瓦片数据,根据函数名称调用sub函数进行运算,将计算结果进行返回,以完成Node3的函数参数替换,此时Node3中计算参数Node1已替换为Node1的执行结果,执行步骤（4）。
（4）同样地,执行步骤（2）和步骤（3）完成Node3中计算参数Node2的参数替换,当所有计算参数替换完毕后,执行Node3中的div函数,最终得到工作流的计算结果。
为了满足对栅格瓦片广泛的处理需求,计算结点的处理函数分为基础函数和模型函数。基础函数是对栅格瓦片进行简单的数值运算;模型函数是对栅格瓦片进行一系列复杂的处理操作,其内部可以封装模型算法或参与计算的额外数据。基础函数见表2所示。
表2 用于栅格瓦片处理的基础函数

Tab. 2 Basic function for raster tile processing

函数类型 操作函数名
基本运算 加（add）、减（sub）、乘（mul）、除（sub）
三角函数 正弦（sin）、余弦（cos）、正切（tan）反正弦（asin）、反余弦（acos）、反正切（atan）
其他数学函数 平方根（sqrt）、绝对值（abs）、自然指数（exp）、最大值（Max）、最小值（Min）
全局函数 全局最大值（GMax）、全局最小值（GMin）
指数函数 归一化植被指数（NDVI）、归一化水指数（NDWI）

2.4 全局运算效率优化

很多计算场景,需要影像中的最小、最大值作为参数参与运算,例如运算表达式： B 2 / GMax ( NDVI ( B 5 , B 4 ) )。当栅格影像整个波段的数据都需要参与运算时,将会导致数据获取和计算的压力,难以满足实时的地图服务需求。针对这种情况,本文提出使用重采样算法得到的Overview数据来代替完整影像。常见的重采样算法有最邻近内插法、双线性内插法和三次卷积法内插法等[26]。Overview不仅在直方图上较好地保留了原始影像的全局特征,在相同位置区域上也保持了近似的特征。该方式可以大幅度减少运算数据量,优化地图服务的实时运算效率。下面以计算植被覆盖度（Fractional Vegetation Cover,FVC）为例,讨论如何使用Overview进行实现。
植被覆盖度是指植被植株冠层或叶面在地面的垂直投影面积占植被区总面积的比例[27],是刻画地表植被覆盖的一个重要参数,也是指示生态环境变化的重要指标之一,计算公式如下：
FVC = ( NDVI - NDV I soil ) ( NDV I veg - NDV I soil )
式中： NDV I soil为完全是裸土或无植被覆盖区域的NDVI值,相当于 NDV I min, NDV I veg则代表完全被植被所覆盖的像元的NDVI值,即纯植被像元的NDVI值,相当于 NDV I max。获取 NDV I min和 NDV I max需要计算当前影像内全局NDVI的最小、最大值,此处使用Overview数据进行计算。具体实现步骤如下：
（1）设计模型函数,输入参数为当前栅格瓦片的可见红光波段、近红外波段和Overview数据。
（2）计算Overview数据的NDVI值,为了消除异常值的干扰,根据直方图统计取置信区间5%到95%的端点值作为全局的 NDV I min和 NDV I max。
（3）计算当前栅格瓦片的NDVI值,代入式（2）求得FVC。

3 实例与结果分析

3.1 实验数据与实验环境

本文实验所用的实验数据为从地理空间数据云下载的Landsat 8影像,数据区域涵盖以中国东部的多个省份,时间范围选取2020—2021年,影像标识编号如表3所示。将每景影像的前5个波段合成到一个GeoTIFF文件中,并将文件转化为COG。实验环境使用的操作系统为Centos7系统,开发语言为Python3.7,Redis版本为5.0.8,CPU为i7-8750H,虚拟机集群为三结点,每个结点分配4G的系统内存。
表3 Landsat 8数据列表

Tab. 3 Landsat 8 data list

标识编号 云量%
LC81200362020115LGN00 0.10
LC81190362020092LGN00 0.15
LC81190382020252LGN00 2.12
LC81190372020236LGN00 6.98
LC80960432021013LGN00 7.97
LC81200372020291LGN00 15.71

3.2 GeoTIFF与COG栅格瓦片读取性能测试

为测试从普通GeoTIFF和云优化的COG中读取栅格瓦片的效率,分别将这两种组织格式的影像存储到本地固态硬盘和HDFS集群中,然后按照不同缩放层级顺序读取一定量的栅格瓦片,并计算每张瓦片的平均读取时间,实验结果如表4所示。
表4 影像瓦片读取性能

Tab. 4 Image tile reading performance

存储位置 影像组织格式 每张瓦片平均读取时间消耗(Zoom)/ms
8 10 12 14 16
本地固态硬盘 GeoTIFF 4864 402 82 48 39
COG 84 153 177 42 34
HDFS GeoTIFF 9677 7725 2097 102 68
COG 104 180 207 53 45
分析表4可知,可得出如下结论：
（1）无论影像数据是存储在本地固态硬盘中还是存储在HDFS中,读取COG中栅格瓦片的效率要比GeoTIFF影像数据要高。如图4所示,这是因为普通GeoTIFF中数据是按照条带式组织,当读取某一兴趣区的瓦片数据时需要先获取整行的像元数据,再对符合兴趣区的数据进行裁切,对大量不必要的数据也进行了加载;而COG影像是按照瓦片式进行存储的,首先获取兴趣区相交的切片数据,再对符合兴趣区的数据进行裁切,减少了不必要数据的加载。因此在兴趣区读取的场景下,COG格式组织的影像更具优势,适用于地图服务中栅格瓦片的读取。
图4 2种TIFF数据组织模型

Fig. 4 Two TIFF data organization models

（2）虽然从本地固态硬盘中读取COG组织的影像效率比从HDFS上读取更高,但基于COG的切片式数据组织方式,读取栅格瓦片时,将只需要读取少量数据,在数据量较小的情况下,本地固态硬盘的读取性能优势较不明显。在实际生产环境中,固态硬盘价格相对较高,而且每个结点都需要存储一份数据,无法适用于海量影像的存储分析,而HDFS集群可由若干结点组成,可对外提供统一的集中式网络存储服务,方便数据管理的同时也可以大量降低存储成本。
（3）在COG中重采样可得到Overview的数据,即栅格瓦片金字塔模型。当Zoom层级较小时,可直接读取Overview层级的数据,使得各Zoom层级栅格瓦片获取的时间相对稳定。在读取普通的GeoTIFF时,随着Zoom层级越小,对应的空间范围越大,读取的数据量也变多,使得瓦片最终读取的效率也急剧下降。
综上,在分布式环境下,使用云优化COG格式组织的影像,对于栅格瓦片的读取效率具有较明显的优势。

3.3 计算性能和地图服务性能评估

为了验证在地图服务中使用不同复杂度模型的实时计算性能,在请求地图服务时,增加运算表达式参数,统计在工作流上对计算结点执行的函数时间。实验中所读取的影像存储在HDFS集群上,栅格瓦片Zoom层级为10—16层,进行测试的栅格瓦片数量为26万余个,表5是进行实验的瓦片空间范围及数量。图5是在不同并发请求数下,分别使用NDVI、基于Landsat 8影像的监督分类 [28]和植被覆盖度FVC共3类模型对栅格瓦片计算的平均耗时结果。对图5分析可得出,模型计算耗时和模型的复杂度呈紧密相关,随计算量增大而增加。NDVI属于像元的叠置分析,计算量较少;地物监督分类模型使用最短距离法实现,计算复杂度和分类数量和样本数据相关,属于像元级别的处理,计算量适中;FVC涉及Overview数据的计算,计算过程中要统计直方图并进行合适的置信度截取,其计算量比前二者要大。
表5 影像瓦片测试数量

Tab. 5 Image tile test quantity

Zoom层级 WGS84瓦片空间编号范围 瓦片
数量/个
XMin YMin XMax YMax
10 849 417 856 424 64
11 1699 835 1713 849 225
12 3399 1671 3427 1699 841
13 6799 3343 6854 3399 3192
14 13 599 6687 13 709 6799 12 543
15 27 199 13 375 27 418 13 599 49 500
16 54 398 26 750 54 837 27 199 198 000
图5 不同模型平均计算时间

Fig. 5 Average calculation time for different models

在地图服务请求中,消耗的时间主要集中在栅格瓦片的获取、模型的计算和计算结果的渲染3个步骤中,表6是不同并发数下各类模型在各个步骤的平均耗时,通过分析对比可知,获取栅格瓦片的时间与模型计算所参与的数据量成正相关。NDVI和FVC都需要2个波段的瓦片数据,FVC模型还需要获取该影像最顶层的Overview数据（471像素×481像素）。因为原始影像的大小为7531像素×7681像素,当使用最顶层的Overview数据替代原始数据进行全局计算时,可减少255倍的数据量参与计算。因此,NDVI和FVC模型总体获取的数据量相对较少,地物监督分类模型由于选取了5个波段的数据参与计算所以获取的数据量较多。此外,计算结果进行渲染的时间在不同计算模型下都相对稳定。在当前实验环境下,随着请求并发度上升,单个瓦片的平均获取时间线性增加,而计算和渲染的平均消耗时间在较小并发度下保持线性增加,在并发度较大时增加幅度不大。即在高并发场景下,瓦片获取消耗的时间会成为实时地图服务计算的瓶颈。本文提出的栅格瓦片计算模型具备分布式扩展能力,因此可通过扩展服务结点可使得服务性能进一步提升。总体而言,在并发场景下,结合分布式瓦片预取缓存策略和多服务节点计算的能力,能对栅格瓦片进行实时、高效的获取与计算,符合对地图服务响应需求。
表6 不同并发数下各阶段效率对比

Tab. 6 Comparison of the efficiency of each stage under different concurrency

并发度 计算模型 单个瓦片平均消耗时间/ms
瓦片获取 计算 渲染
10 NDVI 21.0 1.8 17.9
地物分类 57.2 6.9 18.5
FVC 26.1 87.1 20.1
50 NDVI 69.0 21.5 55.9
地物分类 203.1 48.2 49.6
FVC 76.5 274.8 49.3
100 NDVI 120.3 36.0 73.8
地物分类 273.8 47.9 54.5
FVC 108.1 357.3 69.7
200 NDVI 109.2 48.0 87.7
地物分类 259.4 67.1 76.9
FVC 156.3 369.7 91.2
500 NDVI 172.2 59.1 100.6
地物分类 328.9 60.4 99.1
FVC 226.6 418.9 107.4
1000 NDVI 294.5 64.5 113.1
地物分类 532.4 69.3 105.2
FVC 407.3 450.0 119.8

3.4 全局运算的准确性评估

本文提出的栅格瓦片计算模型对瓦片的叠置分析或基于像元的计算结果是完全准确的,对于全局数据的计算存在一定的偏差。为了验证使用重采样数据来替代全局数据的计算准确性,以FVC模型为例进行准确性验证。先使用原始影像计算得到真实的FVC数据,再分别使用3种重采样得到的Overview数据计算近似的FVC数据,随机选取相同地区的100个3像素×3像素并计算其均值,将FVC近似值和真实值进行对比,并计算平均误差率,表7是其中一景影像的误差结果。
表7 FVC计算数据

Tab. 7 VFC calculation data

全局数据 5%置信度值 95%置信度值 平均误差率/%
原始影像 -0.0397 0.4603 -
最邻近内插法overview -0.0401 0.4602 0.11
双线性内插法overview 0.0153 0.4205 7.40
三次卷积内插法overview -0.0060 0.4267 6.41
通过分析表7得出,使用最邻近内插法得到的影像Overview数据计算结果的平均误差率最小,和真实值最为接近。在3种重采样算法中,三次卷积在空间纹理信息相关性方面最接近原始影像,但最邻近内插法的误差最小,这是因为三次卷积通过增加参与内插计算的邻近像元的数目增强图像边缘,使得空间纹理信息的清晰度和均衡化效果较好,以达到最佳的重采样效果,但是会改变原来的栅格 值[29],导致重采样Overview数据计算得到的NDVI结果直方图和原始影像计算得到的NDVI直方图值分布规律发生了较大变化。在求解FVC时,关注的是全局直方图数值分布百分比,而不是空间纹理信息相关性。图6是原始影像和3种重采样Overview数据运算得到的NDVI值分布的直方图,最邻近内插法得到的直方图数值分布规律和原始数据最相似,在置信区间5%到95%的选取上和原始影像计算的得到的NDVI值接近程度高,最终导致FVC的计算结果也和原始数据相差最小。为防止影像数据的干扰,在其他编号的影像数据上也进行相同实验,得到的结果规律一致。
图6 全局NDVI的直方图结果对比

Fig. 6 Histogram comparison of global NDVI

图7是使用3种重采样方法和真实值的对比散点图。对比发现,最邻近内插法计算的到近似值和真实值的相关性较高;3次卷积法和双线性内插法计算结果与真实值的相关程度也较为显著。图8是真实值和使用3种重采样方法计算得到地图服务可视化结果,使用不同重采样的Overview对于地图服务结果可视化存在一定的误差,但正确反映了计算结果的总体趋势,在地图服务目视解译的误差接受范围内,满足栅格瓦片地图服务的实时计算需求,实现了对原始影像的分析。
图7 不同重采样方法计算结果和真实值散点图

Fig. 7 Scatter plots with calculated results of different resampling methods and true value

图8 FVC计算结果的可视化

Fig. 8 Visualization of FVC calculation results

综上所述,在FVC指数的计算中,使用重采样的Overview数据来代替全局数据进行计算时,计算结果不仅在整体误差上较少,而且在空间上也与真实值保持良好的相关性。针对具体的模型,选择合适的重采样方式,可以极大程度上减少误差,得到和真实结果较一致的结果。

4 结论

本文提出了一种栅格瓦片计算模型,将其应用在了地图服务中。通过对地图服务请求中的计算表达式进行解析,将表达式建模构成计算工作流,使用基础函数和模型函数实现对栅格瓦片实时分析处理。在地图服务中可直接处理得到的数据产品可视化结果,实现了对云端海量遥感影像信息的灵活提取,丰富了地图服务的内涵。最后对本文所做的研究内容进行如下总结：
（1）使用分布式瓦片预取策略,优化从云端读取影像数据的效率;通过设置sigmoid函数更新栅格瓦片的生存时间,实现了栅格瓦片的冷热交替,保证栅格瓦片的高效加载。
（2）使用该计算模型提供的基础函数可完成对栅格瓦片的简单计算,基于模型函数可扩展对栅格瓦片更为复杂的分析。本实验中,基于NDVI、地物分类、FVC等不同复杂度的模型,实现了Landsat 8数据产品信息的实时提取,证明该计算模型具有一定的可扩展能力。
（3）当涉及影像的全局运算时,选择重采样的Overview数据优化计算效率,为栅格瓦片进行复杂的实时分析计算提供了可能。实验结果表明,在Landsat 8数据集上,使用最邻近内插法得到的Overview数据计算得到的FVC结果与真实值相比,平均误差率小于1%;在地图服务可视化的效果上,与真实值效果对比无明显区别,适用于栅格瓦片地图服务中。
（4）该栅格瓦片处理模型具备分布式扩展能力,对于栅格瓦片的请求可分发到不同服务结点进行独立处理。在不同并发场景测试下,能够提供稳定的地图服务能力,在500并发数下,总体响应时间优于1 s。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
Brovelli M A, Minghini M, Zamboni G. Public Participation GIS: A FOSS architecture enabling field-data collection[J]. International Journal of Digital Earth, 2015, 8(5):345-363.

DOI

[2]
Haller A, Janowicz K, Cox S J, et al. The modular SSN ontology: A joint W3C and OGC standard specifying the semantics of sensors, observations, sampling, and actuation[J]. Semantic Web, 2019, 10(1):9-32.

DOI

[3]
Qiao X, Li Z, Ames D P, et al. Simplifying the deployment of OGC Web Processing Services (WPS) for environmental modelling-Introducing Tethys WPS Server[J]. Environmental Modelling & Software, 2019, 115:38-50.

[4]
Li W W, Yang C W, Yang C J. An active crawler for discovering geospatial Web services and their distribution pattern - A case study of OGC Web Map Service[J]. International Journal of Geographical Information Science, 2010, 24(8):1127-1147.

DOI

[5]
Wu H Y, Li Z L, Zhang H W, et al. Monitoring and evaluating the quality of Web Map Service resources for optimizing map composition over the Internet to support decision making[J]. Computers and Geosciences, 2011, 37(4):485-494.

DOI

[6]
Li R, Dong G, Jiang J, et al. Self-adaptive load-balancing strategy based on a time series pattern for concurrent user access on Web map service[J]. Computers & Geosciences, 2019, 131:60-69.

DOI

[7]
Herle S, Blankenbach J. Enhancing the OGC WPS interface with GeoPipes support for real-time geoprocessing[J]. International Journal of Digital Earth, 2018, 11(1):48-63.

DOI

[8]
Guan X, Cheng B, Song A, et al. Modeling Users' Behavior for Testing the Performance of a Web Map Tile Service[J]. Transactions in GIS, 2014, 18:109-125.

DOI

[9]
Stefanakis E. Web Mercator and raster tile maps: two cornerstones of online map service providers[J]. Geomatica, 2017, 71(2):100-109.

DOI

[10]
王浩, 喻占武, 曾武, 等. 网络地理信息服务中的空间数据缓存算法研究[J]. 测绘学报, 2009, 38(4):348-355.

[ Wang H, Yu Z W, Zeng W, et al. The research on the algorithm of spatial data cache in network geographic information service[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(4):348-355. ]

[11]
García R, Verdú E, Regueras L M, et al. A neural network based intelligent system for tile prefetching in web map services[J]. Expert Systems With Applications, 2013, 40(10):4096-4105.

DOI

[12]
Chen D, Shams S, Carmona-Moreno C, et al. Assessment of open source GIS software for water resources management in developing countries[J]. Journal of Hydro-environment Research, 2010, 4(3):253-264.

DOI

[13]
Růžička J. Comparing speed of Web Map Service with GeoServer on ESRI Shapefile and PostGIS[J]. Geoinformatics FCE CTU, 2016, 15(1):3-9.

DOI

[14]
Picoli M C, Simoes R, Chaves M, et al. Cbers Data Cube: a Powerful Technology for Mapping and Monitoring Brazilian Biomes[J]. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2020, 3:533-539.

[15]
Kopp S, Becker P, Doshi A, et al. Achieving the full vision of earth observation data cubes[J]. Data, 2019, 4(3):94.

DOI

[16]
Hnatushenko V V, Sierikova K Y, Sierikov I Y. Development of a Cloud-Based Web Geospatial Information System for Agricultural Monitoring Using Sentinel-2 Data[C]// 2018 IEEE 13th International Scientific and Technical Conference on Computer Sciences and Information Technologies (CSIT). IEEE, 2018, 1:270-273.

[17]
Mahammad S S, Ramakrishnan R. GeoTIFF-A standard image file format for GIS applications[J]. Map India, 2003:28-31.

[18]
van Hoek M, Zhou J E, Jia L, et al. A prototype web-based analysis platform for drought monitoring and early warning[J]. International Journal of Digital Earth, 2020, 13(7):817-831.

DOI

[19]
Panidi E, Rykin I, Kikin P, et al. Cloud-desktop remote sensing data management to ensure time series analysis, integration of QGIS and google earth engine[J]. The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 2020, 43:553-557.

[20]
Rouault E. Cloud optimized GeoTIFF[EB/OL]. https://trac.osgeo.org/gdal/wiki/CloudOptimizedGeoTIFF 2018.

[21]
Lehto L, Kähkönen J, Oksanen J, et al. Supporting wide user-base in rasterR analysis-geocubes finland[J]. International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences, 2018, 42(4):329-334.

[22]
张瀚文. 大规模影像数据集瓦片金字塔快速构建技术[D]. 湖南:国防科技大学, 2017.

[ Zhang H W. Research on Rapid Tile Pyramid Construction Technology of Large-scale Image Data Set[D]. Hunan: National University of Defense Technology, 2017. ]

[23]
陆晔, 张伟, 李飞, 等. 一种基于主题时空价值的服务器端瓦片缓存算法[J]. 浙江大学学报(理学版), 2020, 47(1):12-19.

[ Lu Y, Zhang W, Li F, et al. A server-side tile caching algorithm based on theme temporal and spatial value[J]. Journal of Zhejiang University (Science Edition), 2020, 47(1):12-19. ]

[24]
许虎, 聂云峰, 舒坚. 基于中间件的瓦片地图服务设计与实现[J]. 地球信息科学学报, 2010, 12(4):562-567.

[ Xu H, Nie Y F, Shu J A. Design and implementation of a tile map service based on the middleware[J]. Geo-Information Science, 2010, 12(4):562-567. ]

DOI

[25]
Wu H, Chen X, Song X, et al. Cost minimization of scheduling scientific workflow applications on clouds[J]. Concurrency and Computation: Practice and Experience, 2020, 32(5):e5503.

[26]
Goyal M K, Panchariya V K, Sharma A, et al. Comparative assessment of SWAT model performance in two distinct catchments under various DEM scenarios of varying resolution, sources and resampling methods[J]. Water Resources Management, 2018, 32(2):805-825.

DOI

[27]
李苗苗. 植被覆盖度的遥感估算方法研究[D]. 中国科学院研究生院(遥感应用研究所), 2003.

[ Li M M. The Method of Vegetation Fraction Estimation by Remote Sensing[D]. Graduate School of Chinese Academy of Sciences (Institute of Remote Sensing Applications), 2003. ]

[28]
Lillesand T, Kiefer R W, Chipman J. 遥感与图像解译[M]. 北京: 电子工业出版社, 2016:322-326.

[ Lillesand T, Kiefer R W, Chipman J. Remote Sensing and Image Interpretation[M]. Beijing: Electronic Industry Press, 2016:322-326. ]

[29]
鲍文东, 杨春德, 邵周岳, 等. 几何精校正中三种重采样内插方法的定量比较[J]. 测绘通报, 2009(3):71-72.

[ Bao W D, Yang C D, Shao Z Y, et al. Quantitative comparison of three resampling and interpolation methods in geometric precision correction. Bulletin of Surveying and Mapping[J]. Bulletin of Surveying and Mapping, 2009(3):71-72. ]

Options
文章导航

/

〈 〉