基于MPI的大規(guī)模柵格影像并行瓦片化算法

劉世永，陳 犖，熊 偉，吳 燁，李 軍

LIU Shiyong,CHEN Luo,XIONG Wei,WU Ye,LI Jun

國防科學(xué)技術(shù)大學(xué) 電子科學(xué)與工程學(xué)院，長沙 410073

School of Electronic Science and Technology,National University of Defense Technology,Changsha 410073,China

1 引言

隨著衛(wèi)星傳感器技術(shù)以及無人機(jī)航拍技術(shù)的快速發(fā)展，遙感影像的空間和時(shí)間分辨率都有大幅度地提高，單幅遙感影像文件的數(shù)據(jù)量也急劇增加[1]。當(dāng)前主流GIS軟件以及互聯(lián)網(wǎng)地圖應(yīng)用在WebGIS（網(wǎng)絡(luò)地理信息系統(tǒng)）解決方案中都廣泛采用地圖切片（Tile，又稱瓦片）的發(fā)布策略[2]，這種方法通過預(yù)先將原始影像按照一定的切分規(guī)則切割成一張張大小相同的瓦片，以加載小數(shù)據(jù)量瓦片的方式達(dá)到在網(wǎng)絡(luò)帶寬受限的條件下實(shí)現(xiàn)影像的高效顯示[3-4]。但是目前的商業(yè)GIS軟件其地圖切片的生成以及發(fā)布成本較高，操作較復(fù)雜。以行業(yè)內(nèi)知名的ArcGIS系列軟件為例，要將柵格影像進(jìn)行切片不僅要安裝ArcGIS Desktop，而且還要安裝龐大的ArcGIS Server，并且操作十分復(fù)雜繁瑣，極大增加了時(shí)間和物力成本。最重要一點(diǎn)是隨著單幅遙感影像文件的分辨率以及數(shù)據(jù)量的急劇增加，其相應(yīng)的切片數(shù)量會(huì)呈現(xiàn)幾何級數(shù)式的增長，傳統(tǒng)的串行算法以及商業(yè)GIS軟件是通過單機(jī)預(yù)先切好瓦片，再對外提供，這種傳統(tǒng)切片技術(shù)計(jì)算資源利用低下，并且沒有錯(cuò)誤恢復(fù)機(jī)制，一旦切片過程中出現(xiàn)故障，整個(gè)切片工作得推倒重來，無法在原有的進(jìn)度上繼續(xù)進(jìn)行。因此在硬件性能高速發(fā)展的情況下，如何利用高性能計(jì)算技術(shù)，方便快捷并且高效快速地進(jìn)行柵格影像切片，已經(jīng)成為WebGIS地圖應(yīng)用中，快速可視化方面必須要解決的重要問題。

目前柵格影像并行切片方法主要有以下幾種思路：一種是CPU（Central Processing Unit，中央處理器）+GPU（Graphic Processing Unit，圖形處理單元）的方式，利用GPU的計(jì)算能力進(jìn)行并行加速，這種方法并行能力受限于GPU硬件的能力，并行程度有限，而且會(huì)提高系統(tǒng)架構(gòu)成本，例如劉帥等人提出了一種基于CPU+GPU的地圖切片的快速生成方法[5]；另一種是利用分布式集群系統(tǒng)，將切片任務(wù)劃分為多個(gè)子任務(wù)，各子任務(wù)在多個(gè)分布式節(jié)點(diǎn)上同時(shí)進(jìn)行，這種方法可以較方便地將原有的串行方法并行化，但是當(dāng)數(shù)據(jù)規(guī)模比較大的時(shí)候，前期的數(shù)據(jù)分布式存儲以及后期的結(jié)果合并都比較耗時(shí)[6-7]；再一種是基于多線程并行技術(shù)，各線程將任務(wù)劃分成多個(gè)子任務(wù)，每個(gè)子任務(wù)同時(shí)進(jìn)行，這種方法可以充分利用本機(jī)計(jì)算資源，但是由于線程間并行屬于細(xì)粒度并行，各線程共享父進(jìn)程的內(nèi)存空間，容易造成系統(tǒng)的不穩(wěn)定，并且方法的可擴(kuò)展性不是很好，這種方法目前是ArcGIS采用的策略；還有一種是利用并行技術(shù)進(jìn)行海量小影像的瓦片化的批處理操作，這種方法任務(wù)并行劃分簡單，但是局限性大[2]。目前基于MPI（Message Passing Interface，消息傳遞接口）的多進(jìn)程方式進(jìn)行單個(gè)大規(guī)模柵格影像的并行切片研究較少，這種方法利用共享外存的高性能集群，可擴(kuò)展性強(qiáng)，穩(wěn)定性強(qiáng)，可以充分利用多機(jī)計(jì)算資源，由于集群之間共享外存，所以數(shù)據(jù)無需提前分布式存儲，可以極大地提高影像切片的效率。

2 柵格影像瓦片化基本原理

2.1 瓦片切分模型

瓦片最早是由Google Map提出并進(jìn)行應(yīng)用的，其采用特定的切割方式對采用WebMercator投影坐標(biāo)的世界地圖柵格影像進(jìn)行切片[8-9]，由于WebMercator投影方便計(jì)算機(jī)進(jìn)行計(jì)算，隨后各大主流WebGIS和互聯(lián)網(wǎng)地圖應(yīng)用商都采用基于WebMercator投影坐標(biāo)系的方式進(jìn)行切片，例如國內(nèi)的有百度地圖、高德的圖等，國外有Google Map、Bing Map等[10]。如圖1所示，以一張世界地圖為例闡釋W(xué)ebMercator投影方式，假設(shè)地球被套在一個(gè)圓柱中，赤道與圓柱相切，然后在地球中心放一盞燈，把球面上的圖形投影到圓柱體上，再把圓柱體展開，這就形成了一幅墨卡托投影的世界地圖，其原點(diǎn)在經(jīng)緯度（0，0）處。由于理論上南北極是永遠(yuǎn)無法投影到圓柱體上，并且隨著緯度的增高其變形越大，為了方便，Web墨卡托投影忽略了墨卡托投影中南北兩級變形較大的區(qū)域，把橢圓形的地球投影成平面上邊長等于赤道周長的正方形，其大地坐標(biāo)范圍為[-180°，-85.051 128 779 8°，180°，85.051 128 779 8°]，投影坐標(biāo)范圍為[-20 037 508.342 789 2 m，-20 037 508.342 789 2 m，20 037 508.342 789 2 m，20 037 508.342 789 2 m][11]。

圖1 WebMercator投影示意圖

瓦片切分則是基于上述所說的WebMercator投影坐標(biāo)系統(tǒng)，將柵格影像進(jìn)行不同分辨率的切分，每個(gè)分辨率對應(yīng)WebGIS進(jìn)行縮放操作時(shí)相應(yīng)的層級[12]。以圖2為例，level表示瓦片的級別，tileszie為瓦片的長寬，則每一級別共有4level個(gè)tileszie×tileszie大小的瓦片，因此瓦片級別和其對應(yīng)的空間分辨率Resolution滿足公式（1）：

其中R為地球半徑，瓦片劃分采用四叉樹的方式進(jìn)行，即以赤道和本初子午線的交點(diǎn)作為中心，不斷對地圖進(jìn)行四分，直到每個(gè)格網(wǎng)大小為tilesize×tilesize為止[13]。如圖2所示，0級世界地圖由一個(gè)瓦片表示，1級世界地圖應(yīng)由4個(gè)瓦片表示，2級16個(gè)，以此類推。因此當(dāng)對一個(gè)普通柵格影像進(jìn)行切片時(shí)，首先根據(jù)柵格影像的分辨率找到與其分辨率最接近的瓦片級別，而后通過上述所述的世界地圖切分規(guī)則，計(jì)算影像所在該層世界地圖瓦片中的位置，進(jìn)行切片。

圖2 瓦片分級模型

2.2 瓦片切分模型

當(dāng)?shù)貓D被劃分為4level個(gè)瓦片后，需要對瓦片進(jìn)行編碼，存儲在文件系統(tǒng)的相應(yīng)區(qū)域以便瀏覽器前端能夠快速訪問，瓦片編碼以瓦片在瓦片坐標(biāo)系內(nèi)的行列號為基礎(chǔ)進(jìn)行編制而成。瓦片編碼是基于WebMercator投影坐標(biāo)系進(jìn)行的，如圖3所示，瓦片坐標(biāo)系的原點(diǎn)位于左上角，瓦片存儲在文件系統(tǒng)內(nèi)分為三級目錄其中第一級為瓦片級別(level)，第二級為瓦片列號(tx)，第三級為瓦片行號(ty)。以字符串“root/level/tx/ty.png”作為瓦片圖片的唯一編碼，只需要正確解析瓦片的編碼就可以獲取地圖瓦片相應(yīng)的訪問路徑，其中root為瓦片文件的根目錄。

圖3 瓦片編碼示意圖

3 并行切片任務(wù)劃分

并行瓦片化中各個(gè)進(jìn)程的任務(wù)劃分是一個(gè)非常重要的過程，任務(wù)分配的好壞直接決定著最終算法的性能。在進(jìn)行任務(wù)劃分前，需要先將原始影像投影變換到WebMercator坐標(biāo)系得到投影變換結(jié)果影像，而后在這個(gè)坐標(biāo)系下對投影變換結(jié)果影像進(jìn)行數(shù)據(jù)劃分。本文采用一種循環(huán)劃分的分配策略，具體劃分方法如圖4所示，其描述的是一個(gè)以瓦片為單位大小為5×4的柵格影像被8個(gè)進(jìn)程進(jìn)行切分的任務(wù)劃分示例。該任務(wù)劃分方法名為車輪法，其原理就是各進(jìn)程按照進(jìn)程號從小到大的順序，根據(jù)瓦片坐標(biāo)系下瓦片的分布，從左到右，從上到下依次進(jìn)行分配，當(dāng)一個(gè)周期完畢后，接著上一周期的位置重復(fù)進(jìn)行上述類似操作，直至各進(jìn)程將所有瓦片分配完畢，各進(jìn)程所分配得到的所有瓦片構(gòu)成該進(jìn)程的任務(wù)池。

可以通過如下公式計(jì)算得到每個(gè)進(jìn)程任務(wù)池內(nèi)的每個(gè)瓦片的行列號：假設(shè)進(jìn)程的進(jìn)程號為i，n表示進(jìn)程總數(shù)，[tminX,tminY,tmaxX,tmaxY]為影像所覆蓋的瓦片的行列號范圍，則其對應(yīng)瓦片行列號[tx,ty]滿足式（2）：

其中twidth和theight為x和y方向瓦片個(gè)數(shù)，tcount為瓦片總數(shù)，解算方法如式（3）所示：

圖4 并行切片任務(wù)劃分

假設(shè)在WebMercator投影坐標(biāo)系下影像地理范圍為[ominX,ominY,omaxX,omaxY]，其中ominX、omaxY為左上角點(diǎn)坐標(biāo)，omaxX、ominY為右下角點(diǎn)坐標(biāo)。那么瓦片的行列號范圍[tminX,tminY,tmaxX,tmaxY]可通過影像的地理范圍計(jì)算得出，具體解算方法如式（4）所示：

其中originShift=(2×π×6 378 137)/2，即地球周長的一半，Resolution=(2×π×6 378 137)/(tilesize×2level)，tilesize為瓦片邊長大小表示向下取整表示向上取整。

如果當(dāng)tcount%n≠0時(shí)，進(jìn)程號i＜tcount%n的部分進(jìn)程還需處理瓦片行列號[tx2,ty2]滿足以下條件的瓦片：tx2=tminX+(pos2%twidth)，ty2=tmaxY-

賦碼包括篇頭加注和詞類賦碼。其中，篇頭信息有：篇名、文本出處、作者信息等。詞類賦碼則通常使用GoTagger軟件進(jìn)行。目前，自動(dòng)詞類賦碼技術(shù)已經(jīng)成熟，英語文本的軟件自動(dòng)附碼準(zhǔn)確率可以達(dá)到96%以上。

4 數(shù)據(jù)并行讀取

當(dāng)每個(gè)進(jìn)程的任務(wù)被分配好后，各進(jìn)程依次遍歷其所屬的瓦片，根據(jù)瓦片的行列號解算其在原始柵格影像的數(shù)據(jù)讀取位置。假設(shè)rank(i)表示第i個(gè)進(jìn)程，rank(i)當(dāng)前處理的瓦片行列號為[txi,tyi]，首先進(jìn)程rank(i)在先前瓦片輸出目錄中的level文件夾內(nèi)檢查是否已經(jīng)存在名為txi的文件夾，如果不存在則在level目錄下創(chuàng)建名為txi的文件夾。然后rank(i)根據(jù)瓦片行列號以及瓦片級別反算該瓦片對應(yīng)的地理范圍[gminX,gminY,gmaxX,gmaxY]，單位為m，具體解算過程如式（5）所示：

然后求解瓦片地理范圍與投影變換結(jié)果影像地理范圍相交的區(qū)域，假設(shè)相交區(qū)域?yàn)閇gistminX,gistminY,gistmaxX,gistmaxY]，單位為m，然后計(jì)算相交區(qū)域在投影變換結(jié)果影像中的相對位置以及大小，具體解算過程如下：以m為單位假設(shè)投影變換結(jié)果影像的地理范圍為[gimgminX,gimgminY,gimgmaxX,gimgmaxY]，空間分辨率為Geores，那么投影變換結(jié)果影像中以左上角為原點(diǎn)的像素坐標(biāo)系中的讀取位置rx、ry以及讀取大小rxsize和rysize可由式（6）計(jì)算得到：

5 數(shù)據(jù)并行寫入

由于瓦片的分辨率不一定與影像分辨率一致，因此需要將各進(jìn)程讀入的數(shù)據(jù)重采樣到瓦片相應(yīng)的分辨率下，重采樣方法采用最鄰近內(nèi)插法，重采樣后數(shù)據(jù)大小wxsize、wysize計(jì)算方法如式（7）所示：

利用最鄰近內(nèi)插法將柵格數(shù)據(jù)的大小從rxsize×rysize重采樣到wxsize×wysize。

根據(jù)之前假設(shè)，rank(i)當(dāng)前處理的瓦片行列號為[txi,tyi]，首先進(jìn)程rank(i)在先前的txi文件夾目錄內(nèi)檢查是否已經(jīng)存在名為tyi的瓦片，如果存在則跳過該瓦片，回到第七步，從任務(wù)池中取出下一個(gè)瓦片的行列號，進(jìn)行下一個(gè)瓦片的處理；如果不存在則rank(i)利用GDAL（Geospatial Data Abstraction Library，地理空間數(shù)據(jù)抽象庫）[14]類庫的CreateCopy函數(shù)創(chuàng)建一個(gè)空的瓦片文件，解算重采樣數(shù)據(jù)在瓦片文件中的寫入位置(wx,wy)，計(jì)算方法如式（8）所示：

然后各進(jìn)程利用GDAL的RasterIO函數(shù)將wx、wy、wxsize、wysize填入相應(yīng)參數(shù)，將重采樣數(shù)據(jù)并行寫入瓦片文件。詳細(xì)流程請參見圖5。

圖5 算法執(zhí)行流程

6 算法執(zhí)行步驟

主進(jìn)程投影變換結(jié)果影像生成完畢后，各進(jìn)程同時(shí)打開投影變換結(jié)果影像，根據(jù)影像的空間分辨率計(jì)算其所能切出瓦片的最大級別tmaxz和最小級別tminz，判斷當(dāng)前設(shè)置的目標(biāo)瓦片級別level是否大于最大級別tmaxz或小于最小級別tminz，如果大于最大級別tmaxz則將tmaxz，重新賦值給level；如果小于最小級別tminz，則將tminz重新賦值給level。

各進(jìn)程計(jì)算當(dāng)前l(fā)evel級別下投影變換結(jié)果影像地理范圍內(nèi)所覆蓋的瓦片行列號范圍。根據(jù)瓦片行列號范圍，采用車輪法為每個(gè)進(jìn)程分配具體處理的瓦片的行列號，所有待處理的瓦片根據(jù)瓦片的行列號按照從上到下、從左到右的順序構(gòu)成該進(jìn)程的任務(wù)池。各進(jìn)程從任務(wù)池中依次取出瓦片的行列號，根據(jù)瓦片行列號以及級別在瓦片輸出路徑下創(chuàng)建如下目錄以及空瓦片文件：“root/level/tx/ty.png”，其中root為瓦片輸出根目錄，level為目標(biāo)瓦片級別，tx為瓦片列號，ty為瓦片行號，字符串“root/level/tx/ty.png”作為瓦片的唯一編碼，瀏覽器前端可以通過這個(gè)編碼對應(yīng)的URL（Uniform Resource Locator，統(tǒng)一資源定位符）直接訪問相應(yīng)的瓦片數(shù)據(jù)。

空瓦片文件創(chuàng)建完畢后，各進(jìn)程根據(jù)瓦片級別level反算該行列號的瓦片其對應(yīng)的地理坐標(biāo)范圍，根據(jù)地理范圍求解其與投影變換結(jié)果影像的相交區(qū)域，計(jì)算相交區(qū)域相對于投影變換結(jié)果影像中的像素坐標(biāo)以及范圍，而后根據(jù)像素坐標(biāo)及范圍利用GDAL類庫的RasterIO函數(shù)將該部分相交區(qū)域數(shù)據(jù)從投影變換結(jié)果影像中讀取到該進(jìn)程的內(nèi)存空間中。隨后，各進(jìn)程將相交區(qū)域數(shù)據(jù)從投影變換結(jié)果影像的分辨率重采樣到瓦片level分辨率下，最后將重采樣數(shù)據(jù)寫入之前創(chuàng)建好的空瓦片文件中，則當(dāng)前瓦片生成完畢。如果進(jìn)程的任務(wù)池中還有瓦片未曾處理，則依次對下一個(gè)瓦片進(jìn)行處理，否則該進(jìn)程瓦片切分任務(wù)完成。

7 實(shí)驗(yàn)結(jié)果分析

為了保證實(shí)驗(yàn)的準(zhǔn)確性，實(shí)驗(yàn)平臺采用兩臺硬件環(huán)境一模一樣的超微高性能計(jì)算機(jī)，其中一臺搭載Linux CentOS6.3.X86_64操作系統(tǒng)，另外一臺搭載Windows 7操作系統(tǒng)，具體的硬件環(huán)境見表1。

實(shí)驗(yàn)數(shù)據(jù)采用不同大小、分辨率以及不同投影坐標(biāo)系下的的多波段遙感影像進(jìn)行測試，實(shí)驗(yàn)數(shù)據(jù)詳細(xì)信息見表2。

表1 實(shí)驗(yàn)硬件環(huán)境

7.1 評價(jià)ParaTile并行化程度

實(shí)驗(yàn)通過采用表2所述的1.tif、2.tif、3.tif共3幅影像，切片級別采用影像最大切片級別對算法進(jìn)行測試，記錄算法在對各種規(guī)模影像進(jìn)行切片時(shí)執(zhí)行時(shí)間隨進(jìn)程數(shù)目變化情況。實(shí)驗(yàn)結(jié)果為10次實(shí)驗(yàn)的平均值，如圖7所示。

圖7 ParaTile耗時(shí)隨進(jìn)程數(shù)的變化情況

從圖7中可以看到：（1）一定范圍內(nèi)隨著進(jìn)程數(shù)的增加，ParaTile的效率逐步提升，但是隨著進(jìn)程數(shù)的持續(xù)增加，算法耗時(shí)逐步趨于某一穩(wěn)定值，如果進(jìn)程數(shù)繼續(xù)加大，算法速度在某一程度反而會(huì)出現(xiàn)下降。（2）隨著數(shù)據(jù)量的增大，算法的并行化效率愈加明顯，并且趨于算法耗時(shí)穩(wěn)定值的進(jìn)程數(shù)也相應(yīng)更大。其原因如下：隨著進(jìn)程數(shù)的增加，任務(wù)被劃分給更多進(jìn)程執(zhí)行，每個(gè)進(jìn)程對應(yīng)任務(wù)規(guī)模相應(yīng)更小，所以算法整體性能得到提升。當(dāng)進(jìn)程數(shù)目增大到一定程度后，受限于硬盤的讀寫速度，算法性能趨于穩(wěn)定，但是如果繼續(xù)增大那么可能出現(xiàn)各進(jìn)程讀寫競爭的情況，導(dǎo)致性能下降。因此針對不同規(guī)模的影像，應(yīng)采用相應(yīng)合理的進(jìn)程數(shù)才能獲得最優(yōu)的執(zhí)行效率，通過測試表明本文算法針對不同規(guī)模的影像具有良好的擴(kuò)展性。

由于測試數(shù)據(jù)1.tif數(shù)據(jù)量比較小，執(zhí)行速度比較快，所以在圖中的柱形中不易觀察到。2.tif和3.tif雖然數(shù)據(jù)量上差了3倍多，但是從實(shí)驗(yàn)結(jié)果來看，處理2.tif并沒有比3.tif快很多，這主要是因?yàn)閷?.tif進(jìn)行切片時(shí)必須先將它投影變換到WebMercator投影坐標(biāo)系下，這部分是比較耗時(shí)的，所以造成了這種現(xiàn)象。

7.2 ParaTile與ArcGIS性能對比

實(shí)驗(yàn)采用ArcGIS10.2，測試數(shù)據(jù)用上述3幅影像，切片級別都采用影像的最大切片級別，ParaTile并行切片算法進(jìn)程總數(shù)設(shè)置為16，由于ArcGIS支持多線程加速，為了保證實(shí)驗(yàn)可對比性ArcGIS切片參數(shù)中設(shè)置最大線程數(shù)16，實(shí)驗(yàn)結(jié)果為10次測試結(jié)果的平均值，如圖8所示。由于ArcGIS不支持對非WebMercator投影的影像直接切分成WebMercator投影系下的瓦片，因此對應(yīng)2.tif的ArcGIS實(shí)驗(yàn)數(shù)據(jù)無法獲得，從這也能看到當(dāng)前商業(yè)軟件切片方法的復(fù)雜與局限性。實(shí)驗(yàn)表明本文算法在面對不同大小的柵格影像都保持著穩(wěn)定高效的效率，具有良好的線性加速比，特別是隨著影像數(shù)據(jù)量以及分辨率的提高，這種優(yōu)勢相對ArcGIS來說更加明顯。

圖8 ParaTile耗時(shí)隨進(jìn)程數(shù)的變化情況

8 結(jié)束語

本文通過深入研究瓦片模型的組織方式以及瓦片化的關(guān)鍵流程，針對當(dāng)前一些柵格數(shù)據(jù)瓦片化方法存在的問題和不足，提出一種新的基于MPI的并行柵格影像瓦片化算法。算法通過重采樣以及IO兩級并行來提高效率，其性能可以隨著并行文件系統(tǒng)帶寬的擴(kuò)展得到進(jìn)一步提升。本文算法未來進(jìn)一步優(yōu)化的方向可以從集中式共享存儲架構(gòu)轉(zhuǎn)向分布式存儲架構(gòu)，這樣可以解決目前的I/O瓶頸問題。實(shí)驗(yàn)表明本文算法相比現(xiàn)有的瓦片化方法，性能得到了極大地提升，特別是隨著數(shù)據(jù)規(guī)模以及波段數(shù)的增大，優(yōu)勢越明顯。該算法已實(shí)際部署并服務(wù)于某省國土測繪部門，簡化了傳統(tǒng)基于ArcGIS的復(fù)雜切片操作流程，在相同數(shù)據(jù)量下，至少可比傳統(tǒng)工作效率提高8倍以上。

表2 實(shí)驗(yàn)數(shù)據(jù)

[1]李德仁，朱欣焰，龔健雅.從數(shù)字地圖到空間信息網(wǎng)格——空間信息多級網(wǎng)格理論思考[C]//中國地理信息系統(tǒng)協(xié)會(huì)年會(huì)，2003：642-650.

[2]劉晰，張軼，楊軍，等.利用并行技術(shù)的海量數(shù)據(jù)瓦片快速構(gòu)建[J].測繪科學(xué)，2016，41（1）：144-150.

[3]趙大龍，孫恒宇.地圖切片技術(shù)分析與簡單實(shí)現(xiàn)[J].測繪與空間地理信息，2010，33（1）：116-118.

[4]Turner N，F(xiàn)ernandez C，Lessin M.Image tile server：US，US8244770[P].2012.

[5]劉帥.一種基于CPU+GPU的地圖切片的快速生成方法：中國，CN104267940A[P].2015-01-07.

[6]Yi L.Parallel batch-building remote sensing images tile pyramid with MapReduce[J].Geomatics&Information Science of Wuhan University，2013，38（3）：278-282.

[7]杜波.基于MapReduce的柵格地圖切片系統(tǒng)[D].西安：西安電子科技大學(xué)，2014.

[8]Battersby S E，F(xiàn)inn M P，Usery E L，et al.Implications of Web Mercator and its use in online mapping[J].Cartographica the International Journal for Geographic Information&Geovisualization，2014，49（2）：85-101.

[9]Fran?ula N.Web Mercator projection[J].Geodetski List Glasilo Hrvatskoga Geodetskog Dru?tva，2014（1）.

[10]關(guān)雷，劉蕾，郭慧宇.基于瓦片技術(shù)的高分辨率遙感影像快速訪問技術(shù)在測繪生產(chǎn)中的應(yīng)用研究[J].測繪與空間地理信息，2016，39（2）：78-79.

[11]李長春，蔡伯根，上官偉，等.基于Web墨卡托投影的地圖算法研究與實(shí)現(xiàn)[J].計(jì)算機(jī)應(yīng)用研究，2012，29（12）：4793-4796.

[12]劉世永，吳秋云，陳犖，等.基于高層級地圖瓦片的低層級瓦片并行合成技術(shù)[J].地理信息世界，2015，22（6）：51-55.

[13]劉鎮(zhèn).遙感影像瓦片金字塔模型[J].科技創(chuàng)新導(dǎo)報(bào)，2008（6）：199-200.

[14]葛亮，何濤，王均輝，等.基于GDAL的瓦片切割技術(shù)研究[J].測繪與空間地理信息，2014（7）：130-132.