基于天地图的局部影像渲染技术研究及应用

余万一，张 达，冀 虎

(1.矿冶科技集团有限公司，北京 100160；2.中国-南非矿产资源综合利用联合研究中心，北京 102628；3.国家金属矿绿色开采国际联合研究中心，北京 102628)

目前，矿山安全监测系统广泛采用WebGIS来展示监测点的状态及空间位置信息[1-2]，并基于天地图提供的卫星影像，为用户提供更加直观的可视化分析应用[3-5]。但天地图提供的影像底图更新周期长，一般都是1～2年，甚至更长时间。而矿山的实际情况却已经发生了改变，无法及时呈现在卫星影像上，导致不能很好地满足矿山企业的需求。本文通过研究天地图的图文数据及其开放的API，通过无人机航拍矿区现场图像，拼接制图，替换天地图卫星影像局部图层，实现了矿区影像及时更新[6-8]。

1 天地图

天地图是由国家基础地理信息中心研制的国家地理信息公共服务平台，所提供的中国区域内地理信息数据资源是最全的。影像底图可以达到街道级的显示分辨率。天地图服务采用OGC WMTS标准，支持HTTP和HTTPS协议，影像底图的投影类型有经纬度投影和球面墨卡托投影。为了适应不同的业务场景，天地图提供了多种服务调用方式API，包括：地图API、网页API、WEB服务API、数据API等[9]。

本文所做的分析都是基于地图API和网页API来进行的，通过分析地图API加载卫星影像数据的原理以及网页API提供的各项接口服务说明，利用无人机航拍的矿区影像替换原有的卫星影像底图，来达到矿区现场影像实时更新的目的[10]。

2 技术研究

2.1 天地图底图

天地图底图可以选择显示地图、影像、地形，本文采用的是影像底图，而影像底图的来源主要包括：天地图官方提供的卫星影像底图、第三方供应商提供的卫星影像底图、自建卫星影像底图。三种方式的优缺点详见表1。

表1 天地图卫星影像底图来源方式比较

结合目前矿山安全监测实际应用需求及成本控制，最终采用天地图官方提供的卫星影像底图，然后通过无人机航拍来解决官方卫星影像更新慢的问题。

2.2 天地图影像底图加载原理分析

通过分析初始化地图及渲染地图的方式，可知影像底图是不同的瓦片拼接而成的，每个瓦片都是256×256像素png图片。缩放级别是从1到20，随着缩放级别改变分辨率也有相应变化，而每一级加载的瓦片数也是不同的，推算出：瓦片数为22n-1，列数为2n，行数为2n-1，n代表的是缩放级别。坐标原点是西经180度，北纬90度。

2.3 天地图影像底图加载方式

天地图默认是不加载卫星影像底图的，需要自定义卫星影像底图请求地址，拼接局部瓦片完成地图图层的渲染，具体加载方式可以分为以下几步：

1)地图服务选择影像底图图层服务，投影类型采用球面墨卡托投影，服务采用HTTP的方式请求，影像底图图层请求示例：

http：//t0.tianditu.gov.cn/img_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=img&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TILEMATRIX={z}&TILEROW={x}& TILECOL={y}&tk=密钥

其中，参数z为缩放级别，参数x为切片的行号，参数y为切片的列号，密钥在天地图开发者中心申请的，拥有密钥才能使用天地图提供的服务。

2)配置化参数，配置的内容包括：初始化缩放级别、初始化定位的中心点经纬度、是否可拖拽、是否可缩放、操作控件的加载、地图操作事件的监听等。

3)通过初始化、缩放、拖拽可触发可视区域的渲染，采用HTTP的方式请求天地图服务器获取图层瓦片，如果存在瓦片数据则进行拼接，同时用卫星影像图层渲染地图，否则用异常图片渲染地图，提示该区域无影像。加载流程见图1。

图1 地图加载流程图Fig.1 Map loading flow chart

2.4 行列号计算

要完成替换，需要找到经纬度与影像底图切片行列号的对应关系，通过分析天地图影像底图加载的原理及开放的API资源，推算出经纬度转换为行列号的计算公式。

1)行号计算公式

式中，Row为行号，Lat为纬度，Res为相应缩放级别的分辨率，缩放级别与分辨率对应关系见表2。

2)列号计算公式

式中，Col为列号，Lng为经度，Res为相应缩放级别的分辨率，缩放级别与分辨率对应关系见表2。

表2 缩放级别与分辨率对应关系

2.5 航拍图校准

航拍图校准包括两个方面，一方面需要确定航拍图起始和结束经纬度，另一方面需要确定航拍图的切片数量以及切片与行列号的对应关系。

通过分析网页API，可以采用在影像底图加载整张航拍图的方式来确定起始和结束经纬度。方法如下：将无人机航拍的图片处理成矩形，不规则部分用透明图填充。图片处理完成后，通过天地图提供的添加图片图层功能，校准图片在所在影像地图区域的位置，确定航拍图左下角和右上角的经纬度坐标，由经纬度转换行列号算法计算出左下角和右上角所在的行列号，即可得到起始和结束的行和列进而推算出切片数量，并推算出航拍图需要切割的切片和行列的对应关系。

1)切片数据量计算公式

T=[(Re-Rs)+1]×[(Ce-Cs)+1]

式中，T为切片数量，Re为末行值，Rs为首行值，Ce为末列值，Cs为首列值。

2)图片编号与行列对应关系计算公式

Col=(Idxmod (Ce-Cs+1))+Cs

式中，Row为行号，Col为列号，Idx为切片编号，Re为末行值，Rs为首行值，Ce为末列值，Cs为首列值。

2.6 重新渲染影像底图

采用专用切图软件的切片工具对航拍图进行切图，切图大小为256×256像素的PNG图片，不足256×256的部分用透明图像填充。切图排序采用原图从左至右，从上至下，从0开始编号。通过图片编号与行列转换公式计算出行、列号，以此行列号命名图片并存储到本地，命名规则：行号_列号。

通过参考网页API提供的天地图图层加载监听事件，监听加载的切片图层行列号是否在航拍图所在行列号区间内，对于在此区间内的切片图层，利用本地存储的影像进行替换，达到重新渲染图层的目的。重新渲染的流程见图2。

图2 渲染流程图Fig.2 Rendering flow chart

3 天地图局部渲染技术应用

本文所述技术已经应用于某尾矿库安全监测一张图功能上，首先通过无人机航拍整个矿区，然后根据方案设计的算法计算出要替换的行列号和航拍图切图数量及与行列号的对应关系，最后用航拍图库区的卫星影像，以达到及时呈现尾矿库现场情况的目的。下面以天地图缩放级别16级为例来说明在尾矿库安全监测系统上的应用方式及效果。

3.1 区域切片制图

航拍图使用的是该尾矿库排水斜槽区域的影像图，制图步骤如下：

1)将该影像图以图片图层的方式添加到天地图，通过人工调整影像图在天地图上的空间位置，利用天地图的坐标拾取功能，确定航拍图左下角和右上角的经纬度坐标。得出左下角经度为117.783 690，纬度为29.028 330；右上角经度为117.791 520，纬度为29.033 130。

2)通过经纬度坐标计算得出左下角所在的行号为27241，列号为54209；右上角所在的行号为27240，列号为54210。可知首行值为27240，末行值27241，首列值54209，末列值54210。

3)通过首末行列值计算得出切片数量为4个，从左至右，从上至下，从0开始编号，得到编号值0、1、2、3。

4)通过首末行列值和图片编号计算得出图片编号与行列号的对应关系，见表3。

表3 图片编号与行列号对应关系

5)利用专用切图软件提供的切片工具将该影像图切成4个256×256像素的png图片，分别命名为27240_54209.png、27240_54210.png、 27241_54209.png、27241_54210.png，并保存到本地资源文件夹下。最后影像图的切片效果见图3。

图3 尾矿库排水斜槽区域的影像图切图Fig.3 Image cutaway of tailings pond drainage chute area

3.2 区域局部渲染

天地图提供图层加载事件的监听，通过监听图层加载，判断该影像图所在行列号是否在图层加载切片的范围内，如果在此范围内则用本地资源下的影像切片替换原有卫星影像切片，不在此范围内则用原有卫星影像渲染。相关示意代码如下：

// 监听图层加载事件

map.getLayers()[0].addEventListener('tileloadstart'，function(params){

// 判断行列号是否需要用本地图片替换

if(coords.contains(params.coords)){

params.tile.src = '/resources/'+params.coords.x+'_'+params.coords.y+'.png' }})

在天地图缩放至16级，该尾矿库排水斜槽区域的拼接效果见图4。

图4 航拍影像拼接效果图Fig.4 Mosaicking effect of aerial image

4 结论

基于天地图的开放资源，通过分析地图API和网页API及其他相关资源，设计了相应的影像图层局部渲染的技术，解决了尾矿库区域卫星影像图更新慢、无法及时呈现尾矿库区域现有状况的问题，并应用于矿山安全监测一张图监测点的空间位置信息展示。通过现场应用的结果表明本技术切实可行，同时也为局部区域卫星影像实时性要求较高的业务领域提供了参考。