基于Cesium的地表覆盖数据可视化

2021-06-10王宇恒李大成张德朝陈金勇

无线电通信技术 2021年3期

王宇恒，李大成*，张德朝，陈金勇，杨毅，于杰

(1.太原理工大学矿业工程学院，山西太原 030024；2.清华大学地球系统科学系，北京100084;3.中国电子科技集团公司第五十四研究所，河北石家庄 050081；4.太原理工大学物理与光电工程学院，山西太原 030024)

0 引言

可视化是将数据及其隐藏规律直观、浅显地展示给用户，以达到深入挖掘其内涵，诠释复杂数据的目的[1-2]。面向时空地理数据的多维度、高层次可视化表达及分析有重要的应用价值[2]。基于互联网的云计算技术、时空技术和三维技术的地理信息系统软件产品是对遥感地理信息相关数据进行分析、决策及展示的重要工具，将新型的WebGIS分析与展示技术应用到遥感地表覆盖数据的决策与表达上是实现大数据可视化的重要趋势。

地表覆盖(Land Cover)是指地球表面各种物质类型及其自然属性与特征的综合体，科学准确地测定各类物质的空间分布与动态变化，对于研究地球系统的能量平衡、碳循环及其他生物地球化学循环、气候变化等有着十分重要的意义[3-4]。地表覆盖数据是利用航天、航空影像中蕴含的丰富纹理、光谱及时相信息，按照地表覆盖物的物理特征、生物特征等进行分类、统计、标注得到的数据产品[5-6]。其在自然资源调查、地理国情监测、城市规划及环境检测等领域有广泛的应用，这些应用直接或间接地服务于国民经济建设，是基础性地理国情监测数据成果之一，是了解人类活动和全球变化之间复杂相互作用的关键信息来源[7-8]，合理利用并进行深度分析，最后将结果直观地展示显得尤为重要。

地理信息系统(Geographic Information System，GIS)是地理空间相关领域的重要可视化和分析工具，它是在计算机硬、软件支持下，对整个或部分地球表层(包括大气层)空间中有关地理分布的时空数据进行采集、储存、管理、运算、分析和显示的技术系统[9-11]。GIS客户端架构包括C/S(Client/Server)和B/S (Brower/Server)两种，C/S架构是常见的基于客户端软件的地理信息系统，使用不够灵活、扩展性较差、不能跨平台; B/S架构是基于浏览器的网络地理信息系统(Web Geographic Information System，WebGIS)，相比C/S平台在安装部署、扩展及跨平台等方面都有较强的优势，能更好地为大众服务[12-13]。WebGIS将地理信息数据放在服务端，用户通过浏览器发出请求，服务器对请求进行处理后，根据需求返回需要的地图图片或者数据，再通过浏览器接收返回信息并完成显示，可以大大减轻客户端的压力[14-15]。三维的地理信息系统在算法效率和计算能力极速提高的情况下，也得到飞速的发展，并受到越来越多行业的重视。三维GIS的独特之处在于可以模拟现实地球，将二维的图片、数据等在三维地球上显示，更符合人的认知，对空间数据的展示更加形象、直观[16-17]，并能将分析结果直观地展示出来，做到“深入浅出”。但是目前WebGIS平台在此方面有所不足，面对日益丰富的应用需求，亟需搭建基于WebGIS的遥感地表覆盖数据分析与展示的系统框架，开发一种适用于地表覆盖数据可视化的应用平台。

Cesium是AGI(Analytical Graphics Inc)公司开发的世界级的创建三维地球和地图的JavaScript 开源库，是一个地理空间三维制图与虚拟地球的平台。通过Cesium可以实现三维地球和二维地图的创建。Cesium通过WebGL(Web Graphics Library)技术实现图形的硬件加速，可实现跨平台、跨浏览器，并提供动态数据的可视化展现功能[18-20]。本文针对目前GIS平台在展示地表覆盖信息方面存在的问题，利用Cesium开源库设计并开发出一套三维WebGIS，可以实现地表覆盖数据的深入分析以及直观表达，为地表覆盖数据如何合理地分析、展示提供了解决办法，对实现地表覆盖数据潜在的应用价值创造了条件。

1 地表覆盖数据展示系统的框架设计

基于WebGIS的地表覆盖展示系统的开发框架包括数据生成、数据发布、网页开发和网页发布四部分。各个部分需要搭建的开发环境包括:Java运行环境安装、Geoserver地图数据发布服务器安装、Cesium开发平台安装和Node本地Web服务器安装。将开发的网页通过Node.js进行发布。GeoServer依赖于Java的运行环境，本地数据通过GeoServer发布到网络，然后Cesium通过开发接口链接到URL(Uniform Resource Locator)网络地址加载数据，并通过Node服务器实现网页发布。

本文的设计目标是开发一套主要用于地表覆盖相关数据展示与分析的三维WebGIS。设计内容如图1所示，设计过程中地表覆盖数据需进行一定的处理，所涉及技术包括数据的网络发布和信息统计。系统以Cesium开源库为开发平台，针对地表覆盖数据开发功能如图1所示，具体包括地表覆盖数据在三维虚拟地球上的加载，以及叠加数字高程模型(Digital Elevation Model，DEM)地形后的三维显示; 专题统计信息的二维、三维显示; 时间序列地表覆盖数据动态展示。此外，系统还可以实现三维模型数据的空间展示、CZML动态数据展示等功能。

地表覆盖数据可以作为图层在虚拟地球上进行加载，数据图层空间参考为WGS84地理坐标系。本文使用影像服务方式进行加载，需要将本地影像进行网络发布，通过网络接口加载地图服务。在系统中可以通过加载DEM数据来实现地形显示，实现多图层叠加展示，效果更加直观，对空间分析、土地规划有较大的帮助。地表覆盖统计数据是对地表覆盖数据进行数值化表达的结果。将数值可视化的传统方法是制作统计图表，本文依托Cesium强大的可视化性能，结合空间位置信息，不仅设计二维图表进行统计数据的展示，还将统计数据用三维柱状体的形式进行表达，展示效果生动。依靠Cesium时间轴特性，将时间序列地表覆盖数据与时间轴关联，通过调整时间轴来实现数据动态变化，对分析时间序列数据变化规律有较大的帮助。

图1 三维WebGIS设计内容Fig.1 Design content of 3D WebGIS

除专门针对地表覆盖数据所设计功能外，本平台系统还可以实现三维矢量数据、三维倾斜数据等的加载显示，Cesium语言CZML数据的动态显示。矢量数据、倾斜数据需要转换格式为JSON格式，矢量数据按照某字段值拉伸至三维，实现三维模型的显示。CZML是由Cesium提出的专注于描述动态地理场景的语言，其数据格式中最重要的是时间和位置属性，通过时间和位置的对应实现场景动态的效果。

2 地表覆盖分类测试数据

为了验证所设计开发的三维WebGIS对地表覆盖数据分析展示的效果，本文应用一系列测试数据对系统进行测试。主要测试数据为全球30 m地表覆盖数据(FROM-GLC 2017)[21]和自行分类的时间序列地表覆盖数据。地表覆盖数据需要进行网络发布，以实现WebGIS的加载调用。通过统计地表覆盖数据的相关专题信息，结合地表覆盖分类图层，实现在虚拟地球上的三维展示。

2.1 全球地表覆盖分类产品

地表覆盖数据是通过对影像数据中的纹理和光谱等信息进行分类标注得到的数据产品。本文所使用的全球地表覆盖数据是FROM-GLC 2017，不包括南极洲，分类类型包含10类:农田、森林、草地、灌丛、湿地、水体、苔原、不透水层、裸地和冰雪。使用随机森林分类方法，整体分类精度达72.35%，投影方式为WGS_1984地理坐标系统，单幅数据产品覆盖范围为10°×10°。数据生产过程[21-22]如下:

① 从Landsat-8、MOD13Q1提取特征。Landsat-8使用的是OLI(Operational Land Imager)表面反射率数据，当Landsat-8数据云量大于50%时使用Landsat-5数据代替，主要实验数据时间分布是2014—2015年。MOD13Q1数据是MODIS(Moderate Resolution Imaging Spectrometer)的时间序列植被增强指数(Enhanced Vegetation Index ，EVI)数据，数据时间分辨率为16天，生产时间为2015年。

② 采集土地覆盖类型的训练样本。训练样本采集使用Global Mapper软件，采集时需要考虑时间序列的EVI、月温度变化以及降水等因素，要在全球范围内四季节影像图上进行采集。

③ 在全球分区使用随机森林对2017年的Landsat-8数据分类得到初步结果，第一套为不透水层外的所有类别，第二套为不透水层和其他类别; 使用灯光数据过滤不透水层类别，生成的不透水层结果叠加到第一套结果上，生成最终土地覆盖图。

2.2 基于随机森林的时序地表覆盖分类数据

本文通过制作时间序列地表覆盖分类数据，来验证开发平台动态展示的效果。实验拟定4种分类类别：不透水层、水体、植被和裸地。在挑选实验区域时主要考虑以上4种地物，分类方法选择具有快速、稳定以及适合大区域分类等特点的随机森林方法。随机森林是一种有多棵决策树的集成算法[23-27]。它采用随机有放回的方式抽取训练样本，生成多棵决策树，在每棵树分类过程中随机选择特征变量，并根据基尼指数(Gini)最小原则进行分枝，分裂过程不进行剪枝，最终结果由多棵决策树投票决定。随机森林算法的特点体现在随机性上，一是训练样本的随机选择; 二是节点进行分裂时特征的随机选择。相应地带来了抗干扰能力强、抗过拟合能力强等优点，在大范围分类方面有较稳定的效果。Gini值的计算公式:

(1)

式中,k表示有k种分类类别，pi表示选中样本属于i类别的概率。Gini越大，说明样本集合的不确定性就越大，样本越不纯;反之，说明样本集的不确定性小，纯度高。

分类实验区域选择山西省太原市部分地区，地理坐标为东经112°23′～112°35′，北纬37°49′～38°01′。研究区东部和南部多为建筑物，西部为林地,北部有大块裸地，汾河纵贯南北，地物覆盖类型符合预期分类目标。分类实验的源数据选择Landsat-8的2级数据，来源于USGS。Landsat-8的2级数据是Landsat原始数据经过几何校正、辐射校正得到的，一共有10个波段。由于2018年该研究区影像数据云量多、质量较差，经过筛选，最终选择春季的2018年3月21日，夏季的2018年5月8日和2018年5月24日，共3个时间的影像进行分类实验，3幅影像植被有较明显的差异。

随机森林分类是一种监督分类，主要包括分类器训练和分类器分类两大过程。分类器训练具体有样本采集、训练样本数据生成等工作; 分类器分类具体包括测试样本数据生成和调试参数。具体分类步骤如下:

② 分类测试数据处理。将3幅Landsat-8影像都按照所选研究区域进行裁剪，得到3幅大小为512×512像素、特征变量为7(如图2(a)、(b)、(c))，波段数为10的测试样本影像，并生成测试样本文件。

③ 进行随机森林分类实验，设定分类树棵数，特征变量数为波段数目，初始为10。先选择2018年5月24日的测试数据进行试验，经过多次实验对比发现，分类树棵数为20，保留1～7波段进行训练分类的效果较好。在确定分类模式之后，对剩余时间的测试数据进行实验分类，得到不同时间的分类结果，如图2(d)、(e)、(f)所示，其中紫色代表不透水层，蓝色代表水体，绿色代表植被，黄色代表裸地。

本文仅展示不同时间地表覆盖的变化情况，分类精度不做过多考虑，实验结果有较明显的变化，因此本次分类实验可以满足实验要求。

(a) 2018.3.21

(b) 2018.5.8

(d) 2018.3.21

(e) 2018.5.8

(f) 2018.5.24

2.3 地表覆盖数据网络发布

Cesium支持加载的卫星影像分为影像服务和离线瓦片两种[28]，本文使用影像服务方式进行加载，需要将本地影像通过服务器发布到网络。全球尺度30 m分辨率地表覆盖数据体量较大，为方便加载需要将其压缩后通过GeoServer进行网络发布。GeoServer是一个符合开放地理空间信息联盟 (Open Geospatial Consortium，OGC)标准的共享地理空间数据的开源服务器，它支持发布地图数据并进行更新等操作。全球地表覆盖数据网络发布过程如下:

① 首先需要将全球地表分类图块进行拼接，考虑到全球尺度数据量很大，需将分类图块进行金字塔切片。

② 在Python环境下实现分类产品从本地到网络的批量数据上传，并在图层组拼接后进行发布，山西省范围数据如图3 所示。在GeoServer中发布的数据，可以通过Cesium中的数据接口进行加载。

图3 地表覆盖数据Fig.3 Land cover data

2.4 地表覆盖数据信息统计

地表覆盖数据蕴含着丰富的信息，从不同的角度进行加工挖掘可以得到不同的信息。通过统计不同区域的地表覆盖物附着情况，可以研究区域经济发展状况、城镇化建设以及环境绿化形式等问题。本文对中国范围、各省级行政区范围以及山西省各地级市行政区范围的地表覆盖分类数据进行统计，统计过程如下:

① 按地理坐标筛选出中国区域的原始全球30 m地表覆盖数据集数据，将中国区域地表覆盖数据进行图像拼接;

② 用中国国界矢量、各省级行政区边界矢量、山西省各地级市行政区边界矢量裁剪拼接图像，得到相应范围的地表覆盖数据。

③ 最后将中国、各省级行政区(因澳门区域矢量琐碎暂不统计)、各地级市行政区划单位(以山西省为例)的地表覆盖类别百分比进行统计，并以三维实体和二维图表的方式进行展示。

3 地表覆盖数据展示系统功能实现

针对目前GIS系统对地表覆盖相关数据产品展示的不足，本文所设计开发的三维WebGIS主要功能包括地表覆盖数据的图层展示、地表覆盖统计信息的二三维展示以及时序地表覆盖数据的动态展示等。系统主界面左上角是图层面板，可以控制图层的加载、透明度和放置顺序，图层面板下方是功能控件按钮，界面右上角分别是地址搜索键和复位键。Cesium中的虚拟地球为WGS84椭球，为保证加载的图层可以完美贴合椭球，所有图层的地理坐标系都选择WGS_1984。

3.1 地表覆盖数据的图层展示

地表覆盖数据以二维平面的方式进行展示时，经过投影变换后的影像与实际地形地物存在差异，不方便人们理解以及进行相应的分析。本文将GeoServer发布的全球地表覆盖数据在三维虚拟地球上进行加载，控制DEM数据的加载来实现叠加地形的效果。Cesium三维虚拟地球通过加载GeoServer发布的网页地图服务(Web Mapping Service，WMS)来实现数据图层的显示。在叠加DEM即开启地形后，地表覆盖图层会贴合到地形上，效果如图4所示，可以生动地展示地表覆盖情况，便于人们直观地理解与感受，对工程建设、国家自然资源政策制定等都有重要的意义。同时系统还可以调整添加图层的透明度，以及控制显示与移除等。系统开发功能的开启与关闭都是通过监听控件按钮的点击事件来实现。

图4 地表覆盖数据叠加地形局部效果Fig.4 Overlay terrain local effect of land cover data

3.2 地表覆盖统计信息展示

(1) 地表覆盖统计信息的三维展示

将2.4节处理得到的分类统计数据，包括中国、各省级行政单位和山西省各地级市行政单位，在开发展示系统中以柱状体的形式进行三维展示，三维展示的原理是将各类别统计信息与柱状体的高度关联，并且赋予相应类别的颜色，按次序将各类别统计柱状体叠加，在指定的位置放置统计柱状体，如图5和图6所示。

将地图通过鼠标滚轮缩放到全国各省视野范围，三维柱状图会自动切换到各省级行政区的统计信息，如图5所示。点击“分类类别”可以选择需要显示的单一类别统计信息。将地图视野移动到山西区域，放大地图到一定程度，就会展示出山西省各地级市的分类统计数据，如图6所示，同样可以选择显示单一类别的统计信息，点击 “全部类别”选项后会显示全部类别。

将地图不断缩小，将会由山西省各地级市统计显示返回到全国各省级行政单位统计显示，如图6所示。继续放大到中国视野，会返回到全国总的分类统计信息，不同级别切换只需要转动鼠标滚轮调整视野范围。选择单一类别后，转动滚轮调整视野会在选定类别基础上进行级别切换。点击任意柱状体都会在右上角显示具体的分类信息，包括类别和统计数值。地表覆盖数据统计信息的三维展示有直观、生动、灵活的特点，可以结合地理位置展示具体的信息，是传统统计信息展示所不能实现的。

图5 中国及各省级行政区三维统计柱状图Fig.5 3D statistical histogram of China and all provincial administrative region

图6 山西省各地级市三维统计柱状图及单一类别展示Fig.6 Three-dimensional statistical histogram and single category display of each prefecture-level city in Shanxi Province

(2) 地表覆盖统计信息的二维展示

将全国各省级行政区的分类统计数据使用二维图表的形式进行展示。图表绘制使用Plotly.js，二维展示可以与三维展示同时进行。结合地理位置的三维展示与传统二维图表展示的结合，可以实现对统计信息全方位的展示，如图7所示。

二维图表页面中有两个控件，左侧下拉菜单可以选择展示指定省级行政区的地表覆盖统计信息，点击下拉菜单选项，上方图表会自动切换显示所选省级行政区的分类图表，如图8(a)所示，点击右侧“全国统计”按钮可以切换回全国地表覆盖信息统计。鼠标在统计图上悬浮会显示具体的分类信息。点击柱状图表右侧图例，可以隐藏该类别，选择指定组合的类别进行展示，如图8(b)中只选择水体和草地两个类别进行显示。图表下方是全国地表覆盖情况统计饼状图，具体统计信息在图中进行显示。

图7 地表覆盖统计信息三维与二维结合展示Fig.7 3D and 2D display of statistical information of land cover

(a) 指定区域统计信息显示

(b) 指定类别组合显示

3.3 地表覆盖数据动态展示

因为卫星遥感对地观测具有重访性，随着遥感技术的发展，已经积累了大量的各种地表遥感时间序列数据，这些时间序列数据较为真实地反映了地表在长时间范围内的动态变化情况，是研究地表变化检测的重要途径。由遥感影像生产的地表覆盖数据同样具有时间序列特性，通过对时间序列地表覆盖数据进行分析，可以得到地表趋势时空动态变化情况，对土地利用变化、生态环境变化及气候变化等研究都有重要意义。将时间序列地表覆盖数据进行动态展示可以更容易地发现相邻时间数据的变化，直观地体验到时间序列地表覆盖物的差异。

将2.2节中使用随机森林分类方法得到的3幅地表覆盖分类图像，叠加到三维地球上，同时可以叠加地形进行显示，并可以随时间进行动态的播放，将相邻时间图像的变化直观地展示出来。影像动态变化的原理是根据时间变化每次添加两个图层，并且根据时间调整两个图层的透明度，达到图层动态播放的目的。图像播放的速度以及暂停、继续，可以通过界面下方时间轴进行控制，如图9所示，通过调整播放速度以及暂停可以进行更加细节的分析。

图9 地表覆盖数据叠加地形信息的时序动态展示Fig.9 Dynamic display of time series of overlay terrain information of land cover data

3.4 附加拓展功能展示

在生动展示地表覆盖数据的基础上，本系统还开发了三维数据和CZML数据的展示功能。CZML数据通过将位置数据与时间信息关联来实现动态变化模拟。图10是由CZML数据模拟的飞行轨迹，基本原理是通过计算得到飞行轨迹抛物线上间隔点的位置坐标，并给每一点赋予时间信息，从而实现飞行轨迹动态模拟，同时可以通过时间轴控制模拟轨迹的变化速度以及变化位置。

图11是将建筑物矢量数据按高度进行拉伸，使平面数据变成三维模型的立体效果，并且可以保留原矢量数据的字段信息，通过鼠标左键点击模型显示。通过改变限制条件实现不同的渲染，并且可以对三维模型进行提升高度、隐藏等操作。倾斜摄影数据、点云数据都可以在系统中进行展示，如图12所示。

图10 飞行轨迹模拟Fig.10 Flight path simulation

(a) 三维模型显示

(b) 按高度渲染

(d) 提升整体高度

(e) 初始状态

(f) 按高度隐藏

(a) 三维倾斜数据

(b) 三维点云数据图12 三维数据展示Fig.12 Display of 3D data

4 结束语

地表覆盖数据及其反映出的地表覆盖变化是环境变化研究、地理国情监测及可持续发展规划等不可或缺的重要的基础信息，将地表覆盖相关数据直观地表达出来显得尤为重要。本文主要针对地表覆盖相关数据的可视化进行研究，并利用Cesium开发了相应的展示系统，实现了地表覆盖相关数据的生动表达。

① 地表覆盖相关数据产品的生成。应用随机森林方法生成全球地表覆盖数据，并将中国、中国各省级行政区和山西各地级市的具体分类状况进行统计，作为测试数据为开发系统服务。

② 地表覆盖数据的图层展示。传统的地表覆盖数据只在二维平面进行展示，无法感受因地形起伏带来的变化，将地表覆盖数据叠加DEM数据进行展示，不仅可以有更加真实的感受，而且对地理空间分析有重大的帮助，甚至在工程建设、环境保护以及城市规划等领域都有重要的价值。时间序列地表覆盖数据的动态展示对研究地表覆盖时序性变化、展示地表动态变化结果有重要的意义。

③ 地表覆盖相关统计数据的二三维展示。传统二维图表只能根据地名进行表达，三维统计展示则可以加入地理位置信息，直接在相应位置进行展示，同时三维模型的展示更加直观，而且可以随意调整显示视野范围，展示更加灵活。二维图表与三维统计图相结合，使得展示内容更加完整，展示方式更加丰富，有不错的应用前景。

④ 在针对地表覆盖数据展示的基础上进行拓展开发，可以实现多种三维数据、CZML动态数据的展示，进一步丰富了本系统的功能，为后续的研究提供思路上的帮助。

目前系统主要面向地表覆盖相关数据，功能虽有特色但较为单一，预计后续拓展后端处理功能后，有望形成集数据采集、组织、展示、分析等一体化的高层级地表覆盖产品在线服务系统。