地震行业网的电子地图本地化应用研究

郭啟倩，郑 蕊，刘珠妹，任海军

(1.中国地震局第二监测中心，陕西 西安 710054; 2.湖北省地震局，湖北 武汉 430071)

0 引 言

地理信息系统(geographic information system或geo-information system，GIS)作为多源化防震减灾信息平台，提供各类电子地图信息服务。GIS以其强大的制图功能，结合不同类型的属性数据和空间数据，对防震减灾事业中各类需求分布图进行绘制，已广泛应用于地震台站分布图、地质构造图、地震烈度图等方面。

Web GIS是GIS技术的进一步延伸，实现了网络环境下地理信息数据的访问、存储、分析。自从Google提出瓦片地图方案思想以来，瓦片技术被广泛应用于Web GIS开发中，实现了快捷电子地图在线访问。瓦片技术通过对电子地图切片、渲染处理后，存放在服务器中供用户直接访问，瓦片地图数据的缓存减少了空间数据和磁盘的访问次数，有效减轻了服务器压力，减少了网络负载，提高了响应速度。目前，瓦片地图的服务都是基于Web在线的方式，对于地震行业内涉密数据，如何有效地在离线网络环境下充分利用瓦片技术实现涉密数据可视化，结合瓦片专题地图，提供更完善的行业地图服务对于地震行业的发展有着极其重要的意义。

1 相关研究

近年来，国内外各大地图公司都对相关的地图服务进行了不断研究与改进。如Google推出的谷歌地图就提供了瓦片地图服务，同时还为各类电子地图服务爱好者提供了第三方编程接口。与之相对的强力对手就是微软的MSN(Microsoft service network)虚拟地球服务，其中包括了3D影像、多维地图服务等。国内在地图服务方面取得优秀成果的公司有百度地图、高德地图、搜狗地图等，在给用户带来不错体验效果的同时，推动着国内外研究学者对于瓦片地图技术的不间断研究。主要有以下几个主流方向：

(1)基于瓦片技术的地图服务。

王伟等[1]从整合有效资源扩大应用需求的角度，提出了一种基于数据类型分级组织管理的数据库一体化方法。苟丽美[2]等依据切片地图Web服务(open GIS Web map tile service，WMTS)规范，利用Web技术实现了一种RESTful风格类型的瓦片地图服务，从实现的角度给出了地图瓦片服务的体系架构、整体流程和具体结构。殷福忠等[3]通过研究瓦片金字塔模型，从最底层进行了瓦片金字塔地图服务发布技术的研究。

(2)服务器交互性能。

聂云峰等[4-5]从瓦片切图、空间索引及系统部署方面对地图交互访问进行改进，进一步证明了改进后的中间件方法的效率要优于传统的仓库管理系统(Web map service，WMS)。张凡等[6]提出了一种全新的地图服务实现方法，从响应速度方面通过实验验证了方法的优越性。

(3)地图动态更新。

黄祥志等[7]运用动态缓存的方法对WMS缓存技术进行改进，取得了地图更新的效果。杜清运等[8]在研究金字塔模型的基础上，结合瓦片地图服务，解决了瓦片地图在多分辨率切换时地图更新不连续的问题。

(4)瓦片切图工具的开发。

韦胜[9]则利用AE进行了瓦片地图绘制和投影。王小军等[10]提出了一种基于AE开发的切图工具对瓦片地图进行获取，并进行地图拼接。

由于地震行业系统的数据涉密性，行业内对于瓦片地图的研究在涉密数据离线可视化方向并没有进行深层次研究。针对该实际问题，文中将主流网络电子地图服务免费对外开放的应用程序接口API和部分基础地图数据迁移到本地，并通过坐标系转换算法进行地理位置纠偏，利用地图瓦片技术对离线网络环境下使用现今的电子地图服务进行研究，借鉴国内外目前瓦片技术的研究成果，进一步基于主流GIS软件开发地图切片程序，研究基础地图、专题地图和应用程序接口调用服务，有效解决业内基于电子地图的涉密数据可视化问题，并开发基于地震行业网络的电子地图示例系统。

2 离线电子地图技术

在现如今分布式网络环境下，实现空间数据高效访问的关键是合理有效地利用存储空间。传统的Web GIS通过实时发送请求获得地图数据信息，由于实时性特征，导致网络传输压力大，与服务器交互效率低。因此，按照区域请求空间数据，以瓦片为单位组织空间数据满足访问请求，使得访问更加直接、快捷，优化策略多样性明显。

瓦片地图技术是按照一定的数学规则，通过配置固定的多级显示比例尺，把连续比例的地图划分为多级离散比例的地图，在服务器端提前将地图切割成为一定规格的瓦片矩阵(支持矢量格式GIS数据和栅格格式GIS数据)，依据不同的缩放比例，将矩阵数据存储到服务器不同的目录中，进一步建立瓦片地图名称与地图坐标的映射关系。根据用户不同的请求范围，根据范围内的地图坐标找到对应的已生成大小固定的瓦片数据，返回给用户，在客户端拼接显示范围内的地图。

瓦片地图技术直接返回用户请求坐标区域的瓦片地图数据，有效缩短了地图生成和传送时间，提高了系统响应速度，同时静态图片的处理进一步降低了服务器负担，提升了地图浏览的高效性。

2.1 瓦片金字塔模型研究

假设瓦片地图金字塔模型的最底层为第0层，最上层为第N层，每一层的瓦片单元为4N-n(n=0,1,…,N,N为地图的缩放级别)。根据该规则，则第N层的瓦片个数为1，即一个瓦片单元就能看到整个地图的范围。在第N-1层就有4个瓦片单元，每一个单元显示四分之一的地图范围，以此类推，瓦片单元显示的地图范围逐渐变小，整张地图范围的分辨率逐渐增大，如图1所示。

由图1可以看出，随着比例尺逐渐变大，地图分辨率逐渐提高，唯一不变的是每一层所显示的地图范围。金字塔模型总体的构建思想为：首先确定地图的缩放级别N，将分辨率最高的地图置于金字塔模型的最底层，通过切分方法将地图切割成为相同分辨率的矩形瓦片，即第0层的瓦片集合；其次，第0层的每4个像素为一个单元合并成为1个像素，即形成了第1层的瓦片集合，以此类推，逐渐形成全部的瓦片集合。瓦片金字塔构建完成后供服务器进行不同级别的地图服务调用。

图1 瓦片地图金字塔模型

2.2 地图投影及瓦片组织方式

(1) 投影坐标系。

谷歌地图、百度地图均使用Web Mercator投影方式，其与常规的Mercator投影的区别在于：Web Mercator投影方式是将地球模拟为球体，而Mercator投影参照的是椭球模型。投影坐标系以赤道为标准纬线，本初子午线为中央经线，两者交点为坐标原点，东、北为正方向。由于整个图幅X轴、Y轴取值范围应相等，已知地球半径取值为6 378 137 m，赤道周长(2πR)为20 037 508.342 789 2，则X轴、Y轴的取值范围均为[-20 037 508.342 789 2,200 375 08.342 789 2][11]。

(2) 瓦片坐标系。

瓦片金字塔的每一个层级对应一个相应级别的瓦片坐标系，用于组织地图瓦片。谷歌地图的瓦片坐标原点在西经180°，北纬85.051 13°，即Web Mercator投影坐标系的左上角，坐标系横向往东、纵向往南为正，1级谷歌地图将全世界范围划分为4张瓦片，图2(a)是地图级别为2时的瓦片组织方式。百度地图的瓦片坐标原点以经纬度(0,0)点为基准作了偏移量为(-865,15 850)的偏移处理[12]，即地图瓦片坐标原点(0,0)对应Web Mercator投影坐标系的(-865，15 850)，坐标系横向往东、纵向往北递增[13]，图2(b)是地图级别为2时的瓦片组织方式。

2.3 地图瓦片本地化

构建基于行业网的电子地图服务，首要实现的是地图瓦片数据本地化，即在行业网服务器构建瓦片金字塔数据资源。拟提供的瓦片地图资源包括基础地理地图和专业底图，其中基础地理地图是通过Http请求实现瓦片本地化，专业底图是通过地图切片程序自定义生成本地瓦片资源。

0,01,02,03,00,11,12,13,10,21,22,23,20,31,32,33,3

(a)谷歌地图

-2,1-1,10,11,1-2,0-1,00,01,0-2,-1-1,-10,-11,-1-2,-2-1,-20,-21,-2

(b)百度地图

图2 瓦片坐标系

2.3.1 Http请求实现瓦片本地化

Http请求获取地图瓦片主要是通过瓦片坐标值参数实现的。以百度电子地图提供的基础地理地图为例，在浏览器开发者模式加载百度地图过程中查看加载的网络资源，以URL请求http://online2.map.bdimg.com/pvd/?qt=tile&x=752&y=249&z=12&styles=pl&p=0&cm=1&limit=80&v=088&udt=20170926返回的image资源则为百度地图瓦片，其中参数x、y分别为瓦片X坐标和Y坐标，z为地图级别Z。通过换算某一地理位置对应的瓦片坐标值，替换URL中对应的参数：

URL url=new URL(http://online2.map.bdimg.com/pvd/?qt=tile&x={X}&y={Y}&z={Z}&styles=pl&p=0&cm=1&limit=80&v=088&udt=20170926);

若要获取已知经纬度(lon,lat)地点的百度地图瓦片，首先结合墨卡托投影函数将经度、纬度值转换为瓦片X、Y坐标，转换算法如下：

//经度转换瓦片X坐标

X=(int)((lon+180)/360)*Math.Pow(2,Z);

//纬度转换瓦片Y坐标

doublesinLat=Math.sin(Math.PI*lat/180);

double y=0.5-Math.Log((1+sinLat)/(1-sinLat))/(4*Math.PI);

Y=(int)(y*Math.Pow(2,Z));

然后根据X、Y、Z得到获取资源的URL，使用HttpURLConnection链接资源即可获得对应瓦片：

HttpURLConnection conn=(HttpURLConnection)url.openConnection();

conn.setConnectTimeout(100);

conn.connect();

InputStream in=conn.getInputStream();

2.3.2 自定义地图切片程序

地图切片程序是在Visual Studio平台利用C#语言开发的，基于ESRI公司发布的嵌入式GIS组件库ArcGIS Engine 10.0，主要使用的程序接口为IExport、IActive View、IEnvelope等。数据源为Web Mercator投影方式的ArcGIS格式文档(.mxd)。地图切片流程为[14]：首先确定切片图幅地理坐标范围，有两种方式：自定义和利用IEnvelope接口提供的QueryCoords方法计算地图数据最小外包矩形范围。然后以瓦片坐标系原点为起始点，由比例尺、瓦片尺寸及瓦片跨度循环计算每个瓦片的地理范围并输出。

切片程序对各个层级的地图切割完成的瓦片按照瓦片地图级别、相应级别下的列值、相应列下的行值进行分级组织存储，即切图完成后的结果存储为：一级目录下，文件夹以瓦片地图级别Z命名；二级目录下，文件夹以瓦片坐标系X值命名；三级目录下存放瓦片，以瓦片坐标系Y值命名。

2.4 瓦片异步加载

文中研究的瓦片异步加载主要是指离线状态(互联网脱机状态和行业网状态)下的瓦片数据加载。需要的资源包括本地地图瓦片和本地地图API(应用程序结构)，2.3节已经介绍了地图瓦片本地化方法，本节主要介绍地图API本地化方法及本地瓦片地图加载。

2.4.1 地图API本地化

通过网络资源分析的方式获取主流电子地图程序接口，即地图API本地化。以百度地图API为例，在浏览器开发者模式加载百度地图过程中查看加载的网络资源，将Http请求返回的依赖模块API文件、图标素材获取至本地，主要包括基础API JavaScript文件apiv1.3.min.js，基础css样式文件bmap.css，基本图标素材images以及map、oppc、tile、control等基础modules文件。

将地图程序引用的网络API JavaScript地址修改为本地文件路径：

type="text/javascript" src="js/apiv1.3.min.js">，即完成地图API本地化

2.4.2 瓦片地图加载

本地瓦片地图加载就是按照异步加载方式，将按照金字塔模型组织的本地地图瓦片加载至指定位置[15]。以百度地图为例，自定义LocalMapType方法及BMap.TileLayer()对象LocalMapType，定义LocalMapType对象调用的地图瓦片png文件和瓦片坐标系tileCoord的对应关系，实现本地瓦片地图加载：

//本地瓦片加载函数

functionLocalMapType(){

LocalMapType=new BMap.TileLayer();

LocalMapType.getTilesUrl=function(tileCoord, zoom)

{

var x=tileCoord.x;

var y=tileCoord.y;

var strURL="maptile/baidumaps/";

strURL+=zoom+"/"+x+"/"+y+".png";

returnstrURL;

}

}

//定义地图对象

var map=new BMap.Map("allmap");

//加载本地瓦片地图

var localMapType=new LocalMapType();

map.addTileLayer(localMapType);

2.5 离线地图坐标纠偏

网络电子地图服务提供坐标纠偏的API接口，例如谷歌地图，就是将正确的WGS-84坐标发送至地图服务器，由地图服务器上的算法将坐标纠偏为GCJ-02(火星坐标系)坐标，再返回客户端，才能显示到和偏移处理后的地图对应的正确位置上。而这种在线方式的坐标纠偏必须在接入互联网的状态下才能进行。

文中对离线状态下WGS-84坐标转换为GCJ-02坐标的算法进行整理，对地理坐标值进行纠偏处理，解决了离线状态下真实WGS-84坐标在离线地图上显示偏移的问题，转换算法如下：

double dLat=transformLat(wgLon-105.0,wgLat-35.0);

double dLon=transformLon(wgLon-105.0,wgLat-35.0);

double magic=1-ee*Math.Sin(wgLat/180.0*pi)*Math.Sin(wgLat/180.0*pi);

dLat=(dLat*180.0)/((a*(1-ee))/(magic*Math.Sqrt(magic))*pi);

dLon=(dLon*180.0)/(a/Math.Sqrt(magic)*Math.Cos(wgLat/180.0*pi)*pi);

double mgLat=wgLat+dLat;

double mgLon=wgLon+dLon;

以上算法中，a为克拉索夫斯基椭球参数长半轴(637 824 5.0)，ee为克拉索夫斯基椭球参数第一偏心率平方(0.006 693 421 622 965 943 23)，经度差值dLon计算函数transformLat(double x,double y)和纬度差值dLat计算函数transformLon(double x,double y)如下：

double transformLat(double x,double y) {

double ret=-100.0+2.0*x+3.0*y+0.2*y*y+0.1*x*y+0.2*Math.Sqrt(Math.Abs(x));

ret+=(20.0*Math.Sin(6.0*x*pi)+20.0*Math.Sin(2.0*x*pi))*2.0/3.0;

ret+=(20.0*Math.Sin(y*pi)+40.0*Math.Sin(y/3.0*pi))*2.0/3.0;

ret+=(160.0*Math.Sin(y/12.0*pi)+320*Math.Sin(y*pi/30.0))*2.0/3.0;

return ret;

}

double transformLon(double x,double y) {

double ret=300.0+x+2.0*y+0.1*x*x+0.1*x*y+0.1*Math.Sqrt(Math.Abs(x));

ret+=(20.0*Math.Sin(6.0*x*pi)+20.0*Math.Sin(2.0*x*pi))*2.0/3.0;

ret+=(20.0*Math.Sin(x*pi)+40.0*Math.Sin(x/3.0*pi))*2.0/3.0;

ret+=(150.0*Math.Sin(x/12.0*pi)+300.0*Math.Sin(x/30.0*pi))*2.0/3.0;

return ret;

}

3 离线地图服务应用示例

目前，基于地震行业网的电子地图服务已投入使用，为部署在行业网、涉密网的网站提供方便、安全的电子地图服务，如图3所示。

图3 国家地震科学数据共享中心

国家地震科学数据共享中心对内数据共享网站http://10.5.109.26:8080/csds/index.html汇集了地震行业历史、实时和准实时数据，包括测震波形、前兆数据、地震目录、地震专题数据、空间观测数据等，站内数据总量约65 TB。网站部署在地震行业网内，为行业内人士提供便利的数据共享服务。基于地震行业网的电子地图服务为该网站提供百度地图应用接口，方便观测数据的展示、搜索。

西部形变数据分中心建立的大地形变流动观测数据库管理系统，部署于涉密网内。大地形变流动观测数据库存储了区域精密水准观测数据、流动GNSS观测数据、一二等水准点之记、GNSS点之记、重力点之记等涉密数据，而这些地球物理场流动观测数据需要地图数据辅助其展示、应用。基于地震行业网的电子地图服务为该网站提供谷歌卫星地图应用接口，方便观测数据的展示、搜索，如图4所示。

图4 区域精密水准观测数据可视化

4 结束语

针对行业内涉密数据离线可视化问题，利用地图瓦片技术着手研究在离线网络环境下使用现今的电子地图服务，进一步基于主流GIS软件开发了地图切片程序，研究基础地图、专题地图和应用程序接口调用服务，有效解决了业内基于电子地图的涉密数据可视化问题，并开发了相关的地震行业网络电子地图示例系统。下一步将对地图的更新策略及信息访问并发问题进行研究。