地震行業網的電子地圖本地化應用研究

郭啟倩，鄭 蕊，劉珠妹，任海軍

(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軟件開發了地圖切片程序，研究基礎地圖、專題地圖和應用程序接口調用服務，有效解決了業內基于電子地圖的涉密數據可視化問題，并開發了相關的地震行業網絡電子地圖示例系統。下一步將對地圖的更新策略及信息訪問并發問題進行研究。