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

(南京大學(xué) 江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，南京 210023)

1 引 言

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

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

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

2 柵格化算法特征分析

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

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

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

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

3.1 串行算法并行化思路

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

3.2 雙層并行模式

在多邊形柵格化過程中，利用進(jìn)程級并行可實(shí)現(xiàn)算法的全局并行，但在各進(jìn)程內(nèi)部仍是串行處理各個多邊形。可見，柵格化填充部分在利用進(jìn)程級并行后仍具有較大并行潛力。線程級并行由于其共享存儲等特點(diǎn)，可以彌補(bǔ)該不足，從而實(shí)現(xiàn)對進(jìn)程內(nèi)處理過程的加速[11-12]。

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

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

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

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

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

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

Vprocess=∑Vi/nProcess

(1)

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

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

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

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

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

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

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

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

測試環(huán)境為浪潮Linux高性能集群。集群包含4個節(jié)點(diǎn)，每個節(jié)點(diǎn)包含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西安坐標(biāo)系，圖斑總數(shù)為1186萬個，數(shù)據(jù)量為5.5GB。數(shù)據(jù)共有4個類別，分別以分類1、分類2、分類3和分類4表示。原始試驗(yàn)區(qū)的矢量多邊形如圖4所示。

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

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

4.2 并行效率分析

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

4.3 不同并行模式對比

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

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

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

5 結(jié)束語

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

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

參考文獻(xiàn)：

[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].北京：科學(xué)出版社，2001：167-169.

[5] 吳立新，史文中.地理信息系統(tǒng)原理與算法[M].北京：科學(xué)出版社，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.