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

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

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

國防科學技術大學 電子科學與工程學院，長沙 410073

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

1 引言

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

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

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

2.1 瓦片切分模型

瓦片最早是由Google Map提出并進行應用的，其采用特定的切割方式對采用WebMercator投影坐標的世界地圖柵格影像進行切片[8-9]，由于WebMercator投影方便計算機進行計算，隨后各大主流WebGIS和互聯網地圖應用商都采用基于WebMercator投影坐標系的方式進行切片，例如國內的有百度地圖、高德的圖等，國外有Google Map、Bing Map等[10]。如圖1所示，以一張世界地圖為例闡釋WebMercator投影方式，假設地球被套在一個圓柱中，赤道與圓柱相切，然后在地球中心放一盞燈，把球面上的圖形投影到圓柱體上，再把圓柱體展開，這就形成了一幅墨卡托投影的世界地圖，其原點在經緯度（0，0）處。由于理論上南北極是永遠無法投影到圓柱體上，并且隨著緯度的增高其變形越大，為了方便，Web墨卡托投影忽略了墨卡托投影中南北兩級變形較大的區域，把橢圓形的地球投影成平面上邊長等于赤道周長的正方形，其大地坐標范圍為[-180°，-85.051 128 779 8°，180°，85.051 128 779 8°]，投影坐標范圍為[-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投影坐標系統，將柵格影像進行不同分辨率的切分，每個分辨率對應WebGIS進行縮放操作時相應的層級[12]。以圖2為例，level表示瓦片的級別，tileszie為瓦片的長寬，則每一級別共有4level個tileszie×tileszie大小的瓦片，因此瓦片級別和其對應的空間分辨率Resolution滿足公式（1）：

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

圖2 瓦片分級模型

2.2 瓦片切分模型

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

圖3 瓦片編碼示意圖

3 并行切片任務劃分

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

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

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

圖4 并行切片任務劃分

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

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

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

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

4 數據并行讀取

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

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

5 數據并行寫入

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

利用最鄰近內插法將柵格數據的大小從rxsize×rysize重采樣到wxsize×wysize。

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

然后各進程利用GDAL的RasterIO函數將wx、wy、wxsize、wysize填入相應參數，將重采樣數據并行寫入瓦片文件。詳細流程請參見圖5。

圖5 算法執行流程

6 算法執行步驟

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

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

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

7 實驗結果分析

為了保證實驗的準確性，實驗平臺采用兩臺硬件環境一模一樣的超微高性能計算機，其中一臺搭載Linux CentOS6.3.X86_64操作系統，另外一臺搭載Windows 7操作系統，具體的硬件環境見表1。

實驗數據采用不同大小、分辨率以及不同投影坐標系下的的多波段遙感影像進行測試，實驗數據詳細信息見表2。

表1 實驗硬件環境

7.1 評價ParaTile并行化程度

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

圖7 ParaTile耗時隨進程數的變化情況

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

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

7.2 ParaTile與ArcGIS性能對比

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

圖8 ParaTile耗時隨進程數的變化情況

8 結束語

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

表2 實驗數據

[1]李德仁，朱欣焰，龔健雅.從數字地圖到空間信息網格——空間信息多級網格理論思考[C]//中國地理信息系統協會年會，2003：642-650.

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

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

[4]Turner N，Fernandez 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的柵格地圖切片系統[D].西安：西安電子科技大學，2014.

[8]Battersby S E，Finn 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]關雷，劉蕾，郭慧宇.基于瓦片技術的高分辨率遙感影像快速訪問技術在測繪生產中的應用研究[J].測繪與空間地理信息，2016，39（2）：78-79.

[11]李長春，蔡伯根，上官偉，等.基于Web墨卡托投影的地圖算法研究與實現[J].計算機應用研究，2012，29（12）：4793-4796.

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

[13]劉鎮.遙感影像瓦片金字塔模型[J].科技創新導報，2008（6）：199-200.

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