陳振杰，周琛，李飛雪，李滿春，任沂斌

(南京大學 江蘇省地理信息技術(shù)重點實驗室，南京 210023)

1 引 言

矢量數(shù)據(jù)和柵格數(shù)據(jù)是地理信息系統(tǒng)(GIS)兩種基本的地理數(shù)據(jù)模型[1]。柵格數(shù)據(jù)更適合進行空間分析和空間模擬，在遙感地學分析、空間決策分析等領(lǐng)域被廣泛應(yīng)用[2-3]。在大區(qū)域地理計算中，經(jīng)常需要將矢量多邊形數(shù)據(jù)轉(zhuǎn)換為柵格數(shù)據(jù)。目前，已有不少多邊形柵格化方法的研究，形成了內(nèi)部點擴散法、復數(shù)積分法、包含檢驗法、掃描線算法和邊界代數(shù)法[4-6]等算法。然而，現(xiàn)有算法大都是串行算法，柵格化效率提升面臨瓶頸，難以滿足海量矢量多邊形快速柵格化的需求[7-8]。

隨著并行計算技術(shù)的快速發(fā)展及并行硬件的門檻不斷降低，研發(fā)柵格化并行算法成為解決海量數(shù)據(jù)快速柵格化的有效途徑，一些學者開始研究單個矢量多邊形柵格化算法的并行實現(xiàn)。Healey等人設(shè)計了一種矢量數(shù)據(jù)向柵格數(shù)據(jù)轉(zhuǎn)換的并行算法，并給出矢量柵格化并行算法的選擇標準[9]；Wang設(shè)計并實現(xiàn)了一種基于掃描線法的矢量柵格化并行算法[10]。但將逐一改造現(xiàn)有串行算法不僅工作量大，且效率不高。

本文將在分析典型多邊形柵格化算法的基礎(chǔ)上，研究串行算法并行化思路，提出一種多邊形柵格化算法并行框架，以期解決矢量多邊形柵格化算法快速并行化的問題。

2 柵格化算法特征分析

多邊形柵格化一般過程可歸納為：逐個遍歷所有多邊形，判定每個多邊形內(nèi)部及邊界上的柵格單元，將多邊形屬性值賦給這些柵格單元。典型的多邊形柵格化算法有內(nèi)部點擴散法、復數(shù)積分法、包含檢驗法、掃描線算法和邊界代數(shù)法等。內(nèi)部點擴散法通過重復設(shè)定種子點，填充位于多邊形內(nèi)部及邊界上的種子點柵格，直至多邊形內(nèi)部區(qū)域被填滿(圖1(a))；復數(shù)積分法與包含檢驗法類似，都對每個柵格單元，逐個判定其是否包含在多邊形之內(nèi)，并將多邊形內(nèi)部的柵格單元進行填充(圖1(b)、圖1(c))；掃描線算法通過逐行掃描，識別多邊形內(nèi)部柵格像元條帶，并用多邊形的屬性值將其填充(圖1(d))；邊界代數(shù)法是一種基于積分思想的柵格化方法，通過簡單的加減代數(shù)運算將多邊形屬性值賦給多邊形內(nèi)部及邊界上的柵格單元，實現(xiàn)多邊形的柵格化(圖1(e))。

上述柵格化算法雖然實現(xiàn)原理不同，但具有一些共性：①判定多邊形內(nèi)部及邊界上的柵格單元是多邊形柵格化的關(guān)鍵步驟；②多邊形遍歷不依賴于具體算法；③單個多邊形柵格化的處理都在多邊形最小外包矩形內(nèi)。

圖1 矢量數(shù)據(jù)柵格化算法基本原理示意圖

3 柵格化串行算法并行化方法

3.1 串行算法并行化思路

將串行算法并行化，需要考慮采用何種并行模式、如何將問題域分解。此外，為了能便于將多種柵格化串行算法并行化，需要考慮如何將串行算法的核心功能(或多邊形柵格化基本算子)嵌入進來。在分析矢量多邊形數(shù)據(jù)特點和多邊形柵格化算法特征的基礎(chǔ)上，設(shè)計了一種通用的并行框架。該并行框架主要包括：雙層并行模式、矢量多邊形劃分方法、多邊形柵格化基本算子調(diào)用接口。其中，雙層并行模式是實現(xiàn)算法并行的基礎(chǔ)和核心，通過充分發(fā)揮進程級并行與線程級并行的各自優(yōu)勢，確保盡可能提升算法的并行效率；矢量多邊形劃分方法是問題域劃分的重點，通過對多邊形數(shù)據(jù)進行合理劃分，實現(xiàn)各進程之間的任務(wù)負載均衡；算法調(diào)用接口確保柵格化串行算法能直接嵌入并行框架中，以實現(xiàn)串行算法的快速并行化。

3.2 雙層并行模式

在多邊形柵格化過程中，利用進程級并行可實現(xiàn)算法的全局并行，但在各進程內(nèi)部仍是串行處理各個多邊形?？梢姡瑬鸥窕畛洳糠衷诶眠M程級并行后仍具有較大并行潛力。線程級并行由于其共享存儲等特點，可以彌補該不足，從而實現(xiàn)對進程內(nèi)處理過程的加速[11-12]。

本文將MPI和OpenMP相結(jié)合，設(shè)計了一種進程級和線程級雙層并行模式(圖2)。在該并行模式中，MPI進程負責讀寫數(shù)據(jù)，并對源矢量多邊形數(shù)據(jù)進行劃分；OpenMP線程在各進程內(nèi)部發(fā)揮作用，通過循環(huán)獲取各多邊形及其外包矩形區(qū)域的柵格單元，并對各多邊形調(diào)用多邊形柵格化基本算子，完成對多邊形區(qū)域的柵格化。具體過程如圖2所示，由2個計算節(jié)點中的4個MPI進程分別讀取相應(yīng)的矢量多邊形，進而在各進程內(nèi)部分別派生出2個OpenMP線程，并將多邊形分發(fā)給多線程并行處理；每個OpenMP線程同時處理1個多邊形，通過調(diào)用多邊形柵格化基本算子完成柵格化處理；當多邊形處理結(jié)束后，由各MPI進程將柵格化后的結(jié)果寫入目標柵格文件。

圖2 多邊形柵格化雙層并行模式示意圖

3.3 顧及負載均衡的矢量多邊形劃分方法

針對柵格化雙層并行模式，本文設(shè)計了一種顧及負載均衡的劃分方法，利用多邊形包含頂點數(shù)評估多邊形計算復雜度，并根據(jù)進程數(shù)目進行均勻劃分(圖3)。具體過程如下：

①多邊形復雜度排序。復雜度指標的選擇應(yīng)兼顧代表性和計算簡便。在柵格化過程中，頂點數(shù)影響柵格化過程中判定與計算的時間。因此，本文將頂點數(shù)作為衡量多邊形復雜度的指標，頂點數(shù)越多，則認為多邊形復雜度越高。通過統(tǒng)計各多邊形的頂點數(shù)，將多邊形按復雜度由大到小排序，形成待處理多邊形隊列。

②進程間多邊形劃分。根據(jù)多邊形復雜度、進程數(shù)目求得各進程處理頂點個數(shù)Vprocess，順序從待處理多邊形隊列中取出多邊形歸入各進程數(shù)據(jù)隊列中，使得各進程數(shù)據(jù)隊列中的多邊形頂點數(shù)之和不超過Vprocess，使各進程待處理的多邊形復雜度基本相當。其中，Vprocess的計算方法如下：

Vprocess=∑Vi/nProcess

(1)

其中，i為多邊形編號，Vi為第i個多邊形的頂點數(shù)，nProcess為進程數(shù)目。

③線程間多邊形動態(tài)分配。在各MPI進程內(nèi)部，動態(tài)地將待處理的多邊形分配各OpenMP線程。各進程數(shù)據(jù)每次只分配給各線程1個待處理多邊形，并標識已分配的多邊形；當某個線程處理完后，再從進程數(shù)據(jù)隊列取出1個多邊形分配給該線程，直至所有待處理多邊形被處理。這樣，即使進程數(shù)據(jù)隊列中的多邊形復雜度有較大差異，也可確保線程間的負載均衡。

圖3 矢量多邊形劃分方法示意圖

3.4 多邊形柵格化基本算子調(diào)用接口

多邊形柵格化基本算子用于實現(xiàn)單個多邊形的柵格化，在多邊形柵格化過程中被反復調(diào)用，是多邊形柵格化的核心和基本功能單元。為確保不同的多邊形柵格化基本算子能直接嵌入到并行框架中，必須規(guī)范多邊形柵格化基本算子的調(diào)用接口。多邊形柵格化基本算子內(nèi)部以串行方式執(zhí)行，通過被并行框架調(diào)用，實現(xiàn)串行算法的并行化。因此，多邊形柵格化基本算子可直接從現(xiàn)有的多邊形柵格化串行算法中提取。

通過對常見多邊形柵格化算法參數(shù)的分析，設(shè)計了多邊形柵格化填充函數(shù)接口參數(shù)。參數(shù)包括：①多邊形頂點信息：將頂點X坐標值與Y坐標值分開存儲，分別用double*padfX、double*padfY表示。②多邊形屬性值：用double*pValue表示。③目標柵格數(shù)據(jù)：用char*pCBData表示，用于臨時存儲多邊形外包矩形內(nèi)的目標柵格數(shù)據(jù)，柵格化后函數(shù)返回更新后的目標柵格數(shù)據(jù)pCBData。任何多邊形柵格化填充算法，只要滿足上述接口規(guī)范的要求，且程序開發(fā)語言相同，就可直接嵌入到并行框架中，實現(xiàn)算法的并行化。

4 測試結(jié)果與分析

4.1 測試環(huán)境與測試數(shù)據(jù)

測試環(huán)境為浪潮Linux高性能集群。集群包含4個節(jié)點，每個節(jié)點包含2顆四核八線程Intel Xeon CPU、12GB內(nèi)存、292GB硬盤、千兆以太網(wǎng)卡。軟件配置為Centos Linux 6.0操作系統(tǒng)、Openmpi 1.4.1、OpenMP 2.0、PostGIS 2.0.3+PostgreSQL 9.0.5數(shù)據(jù)庫。

測試數(shù)據(jù)為某省土地現(xiàn)狀地類圖斑數(shù)據(jù)。數(shù)據(jù)格式為PostGIS格式，數(shù)據(jù)空間參考系為1980西安坐標系，圖斑總數(shù)為1186萬個，數(shù)據(jù)量為5.5GB。數(shù)據(jù)共有4個類別，分別以分類1、分類2、分類3和分類4表示。原始試驗區(qū)的矢量多邊形如圖4所示。

圖4 矢量多邊形測試數(shù)據(jù)示意圖

本文選取掃描線算法和邊界代數(shù)法兩種常用多邊形柵格化算法作為算法用例，利用本文形成的并行框架將掃描線與邊界代數(shù)串行算法快速并行化。試驗設(shè)定輸出柵格單元大小為20m×20m，輸出的柵格數(shù)據(jù)大小為25880×24225。試驗結(jié)果表明，經(jīng)并行化后的掃描線并行算法與邊界代數(shù)并行算法均運行穩(wěn)定，且結(jié)果正確。

4.2 并行效率分析

運行時間和加速比是度量算法并行效率的重要指標。其中，算法的運行時間是從算法啟動直至最后1個進程執(zhí)行完所花費的時間；加速比是同一個任務(wù)在串行及并行環(huán)境下運行時間的比值[13]。通過測試本文實現(xiàn)的不同并行算法的運行時間和加速比，比較并行算法較于串行算法的性能提升和不同算法并行效率的差異。其中，并行算法計算單元數(shù)目用NP×NT數(shù)表示，NP表示并行算法使用的MPI進程個數(shù)，NT表示OpenMP線程個數(shù)。對于同一NP×NT數(shù)，常有不同的進程與線程組合，本文測試采用取得最佳效率的組合來表示。測試結(jié)果如圖5所示。圖5(a)、圖5(b)為不同并行算法與串行算法的對比結(jié)果。從中可看出，利用本文的并行化方法改造后的多邊形柵格化算法大大減少了柵格化時間，取得了良好的加速比，從而驗證了本文設(shè)計的并行化方法的有效性。圖5(c)、圖5(d)為兩種并行算法的效率對比結(jié)果。結(jié)果表明，邊界代數(shù)法較掃描線算法更為快速地實現(xiàn)柵格化轉(zhuǎn)換，取得更好的并行效率。原因在于邊界代數(shù)法以整個柵格單元為基本操作單位，從而實現(xiàn)對柵格矩陣的直接操作；掃描線算法則需要進行掃描線與多邊形交點的多次判別，因而效率低于邊界代數(shù)法。

4.3 不同并行模式對比

本文并行算法采用雙層并行模式實現(xiàn)，為比較其與常規(guī)方法的效率差異，本文實現(xiàn)了基于MPI并行模式的柵格化算法，并將兩種并行模式進行了效率對比，測試結(jié)果如圖6所示。從圖中可明顯地看出，對于每一種算法，雙層并行模式效率均優(yōu)于MPI并行模式。對于MPI并行模式，當進程數(shù)超過32時運行時間緩慢回升；但對于雙層并行模式，隨著NP×NT數(shù)的增大，運行時間一直降低，直至達到64。原因在于本文實驗集群共包含32個使用超線程技術(shù)的計算核心，對于MPI并行模式，可用計算單元數(shù)為32；對于雙層并行模式，可用計算單元數(shù)為64。可見本研究所選的雙層并行模式可以更好地利用系統(tǒng)資源，能獲得更好的并行效率。

圖5 并行算法效率測試結(jié)果

圖6 算法并行模式效率對比結(jié)果

5 結(jié)束語

本文在分析典型多邊形柵格化算法的基礎(chǔ)上，提出多邊形柵格化算法并行框架，從而形成了一種多邊形柵格化串行算法快速并行化的方法，為解決地理計算串行算法快速并行化問題提供有

益的借鑒。試驗表明，本文設(shè)計并實現(xiàn)的并行化方法可實現(xiàn)柵格化串行算法的快速并行化，利用此方法改造后的算法具有良好的穩(wěn)定性和加速比。在后續(xù)的研究中，將進一步探索該并行框架在其他矢量數(shù)據(jù)處理中的適用性，并拓展和深化相關(guān)研究與應(yīng)用。

參考文獻：

[1] MAGUIRE D J，GOODCHILD M F，WRHIND D.Geographical information systems：Principles and applications[M].England：Longman Scientific and Technical，1991.45-54.

[2] GOODCHILD M F.Scale in GIS：An overview[J].Geomorphology，2011，130(1-2)：5-9.

[3] CONGALTON R G.Exploring and evaluating the consequences of vector-to-raster and raster-to-vector conversion[J].Photogrammetric Engineering and Remote Sensing，1997，63(4)：425-434.

[4] 龔健雅.地理信息系統(tǒng)基礎(chǔ)[M].北京：科學出版社，2001：167-169.

[5] 吳立新，史文中.地理信息系統(tǒng)原理與算法[M].北京：科學出版社，2003：202-212.

[6] JIMéNEZ J J，F(xiàn)EITO F R，SEGURA R J.A new hierarchical triangle-based point-in-polygon data structure[J].Computers and Geosciences，2009，35(9)：1843-1853.

[7] WANG S，ARMSTRONG M P.A quadtree approach to domain decomposition for spatial interpolation in grid computing environments[J].Parallel Computing，2003，29(10)：1481-1504.

[8] GUAN Q，CLARKE K C.A general-purpose parallel raster processing programming library test application using a geographic cellular automata model[J].International Journal of Geographical Information Science，2010，24(5)：695-722.

[9] HEALEY R，DOWERS S，GITTINGS B，et al.Parallel processing algorithms for GIS[M].London：Taylor and Francis，1997.

[10] WANG Y，CHEN Z，CHENG L，et al.Parallel scanline algorithm for rapid rasterization of vector geographic data[J].Computers and Geosciences，2013，59：31-40.

[11] BOVA S W，BRESHEARS C P，CUICCHI C E，et al.Dual-level parallel analysis of harbor wave response using MPI and OpenMP[J].International Journal of High Performance Computing Applications，2000，14(1)：49-64.

[12] JIN H，JESPERSEN D，MEHROTRA P，et al.High performance computing using MPI and OpenMP on multi-core parallel systems[J].Parallel Computing，2011，37(9)：562-575.

[13] XIE J.Implementation and performance optimization of a parallel contour line generation algorithm[J].Computers and Geosciences，2012，49：21-28.