矢量多邊形并行柵格化數據劃分方法*

周琛，李滿春，陳振杰，姜朋輝，陳東(南京大學地理與海洋科學學院，江蘇南京210023)

矢量多邊形并行柵格化數據劃分方法*

周琛，李滿春，陳振杰，姜朋輝，陳東
(南京大學地理與海洋科學學院，江蘇南京210023)

針對多邊形并行柵格化中的負載不均衡問題提出一種新的數據劃分方法，主要包括:迭代計算劃分線的位置，在每次迭代中保證分塊間的計算量大致均衡，完成數據劃分、實現負載均衡；提出基于二叉樹的劃分結果融合策略，以解決跨邊界多邊形的融合問題。在多核CPU環境下實現并行算法，選用多個典型土地利用現狀數據集進行測試。結果表明:針對不同類型多邊形數據集，所提方法較傳統方法可獲得更高的并行加速比和更好的負載均衡；針對大數據量數據集，以多邊形節點數為度量標準可更精確地估算分塊計算量，從而更好地實現負載均衡。

地理信息系統；并行計算；多邊形柵格化；數據劃分；負載均衡

矢量數據和柵格數據是地理信息系統(Geographic Information System，GIS)中的基本數據類型［1］。柵格數據更適合進行空間分析和空間模擬，能夠高效地處理空間尺度問題，因此經常需要進行矢量多邊形數據的柵格化處理［2］。近年來，隨著對地觀測技術的快速發展，利用并行計算技術實現對大規模多邊形數據的快速、實時處理顯得十分迫切和必要［3-7］。在多邊形并行柵格化中，數據劃分方法的優劣將極大地影響各劃分分塊計算量的均衡性，進而影響并行計算效率［8-9］。同時，多邊形具有數據量大、形態各異和復雜度差異大的特點［10］，這對研究負載性良好的數據劃分方法提出了挑戰。

傳統的數據劃分方法包括基于多邊形ID順序和基于空間位置的劃分方法。基于ID順序的劃分方法根據多邊形ID的存儲順序均勻劃分成多個分塊［11-13］；該方法易于實現，但劃分較為粗略，忽略了多邊形復雜程度不同對并行效率的影響，因而效率不高。基于空間位置的劃分方法根據數據集的空間位置進行規則劃分，包括行劃分、列劃分、格網劃分和四叉樹劃分等［14-18］。在此基礎上，Lee等［19］提出了一種啟發式劃分方法，可將給定的空間范圍劃分成任意個面積相等的格網。該方法實現速度快，且能在一定程度上保證多邊形的空間聚集性，因而應用廣泛。然而，該方法采用面積相等作為劃分的標準，忽略了多邊形的大小、形狀等特征對并行效率的影響，因而很難實現負載均衡。范俊甫等［13］針對不同多邊形圖層疊置分析的并行處理提出了一種分組間關聯最小化的劃分方法，通過將所有相交多邊形分組實現對多邊形數據的劃分。該方法劃分規則明確，但僅適用于多邊形疊置分析，不具有通用性，且極易造成不同分組間的計算量失衡，導致數據傾斜。為此，提出一種負載性較好的多邊形數據劃分方法。

本文提出一種改進的基于啟發式劃分的數據劃分方法，主要包括:①基于傳統啟發式劃分方法，將多邊形節點數或多邊形數作為劃分的度量標準，通過迭代計算劃分線的位置，進而完成數據劃分；②提出一種基于二叉樹的劃分結果融合策略，將空間上相鄰的分塊依次進行兩兩融合，以解決跨邊界多邊形的融合問題。在多核CPU環境下實現并行柵格化算法，選用多個典型的中國土地利用現狀數據集進行測試，并從運行時間、并行加速比和負載均衡三個方面對數據劃分方法的有效性和穩定性進行評價。

1 算法并行性分析

典型的多邊形柵格化算法包括內部點擴散法、復數積分法、掃描線算法和邊界代數法等［20-21］。內部點擴散法通過重復設定種子點，填充位于多邊形內部及邊界上的種子點柵格，直至多邊形內部區域被填滿；復數積分法對每個柵格單元逐個判定其是否包含在多邊形之內，并將多邊形內部的柵格單元進行填充；掃描線算法通過逐行掃描，識別多邊形內部柵格像元條帶，并用多邊形的屬性值將其填充。邊界代數法通過加減代數運算將屬性值賦給多邊形內部及邊界上的柵格單元。其中，邊界代數法實現簡便、運算速度快，本文主要實現邊界代數法的并行柵格化。

多邊形柵格化算法實現原理不同，但具有相同的并行性，表現為:主要過程為判定多邊形內部及邊界上的柵格單元；對單個多邊形的處理都在最小外包矩形內，不涉及其他多邊形；不依賴于具體的柵格化填充算法；對多邊形節點數敏感性極強。上述分析表明，多邊形柵格化屬于數據密集型的局部計算類型，具有良好的可并行性。

2 改進的啟發式數據劃分方法

2.1 改進的啟發式劃分過程

本文基于傳統的啟發式劃分方法，將分塊計算量相等作為劃分的標準，通過重復迭代計算劃分線位置使得劃分后各分塊計算量大致相當，從而實現負載均衡。在估算計算量時可采用多邊形節點數或多邊形數作為度量計算量的標準。改進后的啟發式劃分方法的基本過程如下(見圖1):

圖1 改進的啟發式劃分方法示意圖Fig.1 Improved heuristic decomposition method

步驟1：確定分塊數p、最大迭代次數n、迭代中劃分線每次移動的距離d和分塊間節點數(或多邊形數)相差閾值s；

步驟2：若p=1則停止迭代計算；

步驟3：若p為偶數，則將待劃分區域沿著寬度較長的邊劃分成面積相等的分塊A和B；通過空間查詢分別獲取A和B的節點數(或多邊形數)NA和NB，若NA-NB≤s，則滿足要求，進行下一輪迭代。否則，需要對劃分線位置進行調整，具體調整過程如下:①當NA-NB＞s時，將劃分線的位置向使得分塊A面積減少的方向移動距離d；②當NB-NA＞s時，將劃分線的位置向使得分塊B面積減少的方向移動距離d；③重復過程①、過程②直至滿足NA-NB≤s或達到最大迭代次數n，在迭代過程中若當前劃分線移動方向與上一次移動方向相反，則改變d的值，使得d=d/2。

步驟4：若p為奇數，則將待劃分區域沿著寬度較長的邊劃分成分塊A和B，使得兩者面積比為［p/2］:［p/2］+1；通過空間查詢分別獲取A和B的節點數(或多邊形數)NA和NB，若NA-(［p/2］:［p/2］+1)NB≤s，則滿足要求，進行下一輪迭代。否則，需要對劃分線位置進行調整，具體調整過程如下:①當NA-(［p/2］:［p/2］+1)NB＞s時，將劃分線的位置向使得分塊A面積減少的方向移動距離d；②當(［p/2］:［p/2］+ 1)NB-NA＞s時，將劃分線的位置向使得分塊B面積減少的方向移動距離d；③重復過程①、過程②直至滿足NA-(［p/2］:［p/2］+1)NB≤s或達到最大迭代次數n，在迭代過程中若當前劃分線移動方向與上一次移動方向相反，則改變d的值，使得d= d/2。

步驟5：對子分塊A和B重復步驟1至步驟4，直至劃分完畢。

2.2 分塊處理結果融合策略

在上述啟發式劃分后，各任務分塊的邊界處存在大量跨邊界多邊形。若該類型多邊形不經處理，則會導致該類型多邊形的重復處理，從而降低并行效率。本文提出一種結果融合策略，對并行處理后的跨邊界多邊形進行快速融合，主要包括兩個步驟:基于二叉樹結構的啟發式劃分結果構建及分塊跨邊界多邊形的迭代處理。

首先，根據啟發式劃分迭代劃分空間位置的特性逐級構建二叉樹，每次劃分當前空間形成的兩個分塊分別為當前層級的左右結點(如圖2(a)所示)。這樣，當指定分塊數為p時，形成的二叉樹最大層級為［log2p］。當上述二叉樹構建完畢后按照二叉樹層級從底端逐層向上進行迭代計算，迭代次數為［log2p］。在每次迭代中，參與計算的分塊為當前層級包含的分塊(如圖2(b)所示)。主要過程如下:①不屬于該層級的分塊直接進入下一次迭代；②在當前二叉樹層級中，將擁有相同根結點的右結點分塊中的跨邊界多邊形傳遞給處理左結點分塊的進程，并由該進程負責兩個分塊中跨邊界多邊形的融合；③左結點進程剔除兩分塊中的相同多邊形，并將僅與該分塊有交集的跨邊界多邊形寫入目標文件，將跨多個分塊的多邊形保留，進入下一次迭代，在下一次迭代中，該進程作為根結點的虛擬處理進程，其處理的兩個分塊整體作為根結點的虛擬分塊；④重復步驟①～③，直至迭代結束。該策略可保證每次參與融合的兩個分塊在空間上鄰近，且融合次數最少。

2.3 并行算法實現流程

多邊形并行柵格化過程可分為預處理、并行執行和后處理過程。預處理過程包括多邊形數據劃分和任務分發；并行執行過程，即對各多邊形進行并行柵格化填充計算；后處理過程針對跨邊界多邊形進行結果融合處理。并行算法采用標準C++編程語言在Linux開發平臺下開發，并在消息傳遞接口(Message Passing Interface，MPI)并行環境下實現，矢量數據的讀寫操作通過地理數據處理類庫(Geospatial Data Abstraction Library，GDAL)實現。具體并行實現流程總體上分為以下步驟(如圖3所示):

步驟1：并行環境初始化，接收數據劃分參數，包括計算量度量標準(多邊形節點數或多邊形數)、分塊數(進程數)p、最大迭代次數n、迭代中劃分線移動距離d和分塊間節點數(或多邊形數)相差閾值s。

步驟2：讀取源矢量多邊形數據，并獲取最小范圍內的空間位置。

步驟3：對多邊形數據集進行改進的啟發式劃分，完成數據劃分，并將劃分結果傳遞給各并行進程。各并行進程分別處理1個任務分塊。

步驟4：各并行進程讀取各自任務分塊中的多邊形，并調用多邊形柵格化算法進行并行計算。

步驟5：并行執行過程結束后構建基于二叉樹的啟發式劃分結果，并根據二叉樹結構進行迭代計算，將空間上鄰近的分塊兩兩進行融合，以解決跨邊界多邊形的融合問題。

圖2 劃分結果融合策略示意圖Fig.2 Fusion strategy for the decomposed results

步驟6：輸出最終結果并退出并行環境。

圖3 并行算法實現流程圖Fig.3 Implementation flow of the parallel polygon calculation algorithm

3 實驗結果與分析

3.1 性能評價方法

本文將多邊形數據集分為三類:不同數據量、不同空間分布和不同數據復雜度。其中，數據量用數據集的內存占用量和多邊形數表示；空間分布用數據集實際面積與其最小外接矩形面積的比值表示；數據復雜度用平均多邊形節點數表示。本文分別采用運行時間、加速比和負載均衡指數來評價算法的并行性能。其中，運行時間是從算法啟動直到最后一個進程執行完所花費的時間。加速比是同一個任務在串行環境下和并行環境下運行時間的比值，如式(1)所示。

其中，Tsequential為串行時間，Tparallel為并行時間。負載均衡指數等于并行進程執行的最長時間與最短時間的比值，該指數越接近于1，表明進程間運行時間越接近、算法負載性能越好，如式(2)所示。

本文將所提數據劃分方法與傳統方法進行對比，分別測試三類多邊形數據集，并從并行運行時間、加速比和負載均衡性能三個方面對算法并行性能進行評價。

3.2 并行環境與實驗數據

程序運行選擇IBM并行集群，包含5個計算節點，每個節點的硬件配置為:CPU 2顆，規格為Intel(R)Xeon(R)CPU E5-2620(主頻2.00GHz，6核12線程)；內存為16GB(4根4GB內存條，規格為DDR3 RDIMM1600MHz)；硬盤為2TB，網絡為集成的雙口千兆以太網。軟件配置:操作系統為Centos Linux 6.3，文件系統為lustre系統，MPI的實現產品選擇OpenMPI 1.4.1。

測試數據為中國不同區域的土地利用現狀數據(見表1)。其中，數據1和2代表不同數據量的數據集，其數據量分別為5.5 GB和1.6 GB，多邊形數分別為12 126 100和2 300 723；數據3和4代表不同空間分布的數據集，實際面積與其最小外接矩形(Minimum Bounding Rectangle，MBR)面積比值分別為54.9%和27.7%；數據5和數據6代表不同復雜度的數據集，其數據量相近、但平均節點數相差很大，分別為35.05和643.15。

表1 測試數據集基本參數Tab.1 Datasets used in the experiments

3.3 不同數據劃分方法性能對比

實驗的目的為比較不同劃分方法的性能，將本文方法與傳統的基于多邊形ID順序的劃分方法和啟發式劃分方法分別應用于并行算法中。在本文方法中，采用多邊形節點數作為度量計算量的標準；測試從數據1至數據6，分別計算運行時間、加速比和負載均衡指數(見圖4)。

圖4(a)描述了對不同數據量數據集的測試結果。不同數據劃分方法表現出相同的變化趨勢:運行時間隨著進程數的增加逐漸減少；加速比逐漸上升，當達到并行環境的最大核數時達到最優；負載均衡指數逐漸降低。這表明上述三種方法均能降低算法的運行時間；比較而言，本文方法能更有效地降低算法運行時間、獲得更高的并行加速比。對數據1，串行時間為1805.34 s，三種方法的最少運行時間分別為126.34 s，122.40 s和101.48 s，最大加速比為14.29，14.75和17.79；對數據2，串行運行時間為756.89 s，最少運行時間分別為54.65 s，53.38 s和45.62 s，最大加速比為13.85，14.18和16.59。同時，對不同數據量的數據集，傳統方法負載均衡指數較大，這表明各進程間計算量均衡性較差；而本文方法負載均衡指數低于其他兩種方法，表明采用本文方法的并行算法各進程間計算量較為均衡。

圖4(b)描述了對不同空間分布數據集的測試結果。對于空間分布較為均勻的數據3，不同方法均表現出較好的性能。對于空間分布不均的數據4，傳統的啟發式劃分方法受空間分布影響較大，僅獲得9.28的加速比；同時，其負載均衡指數較高，高于3.16。原因在于傳統啟發式劃分方法僅能保證各分塊面積相等，而不能保證各分塊計算量相等，從而不適用于空間分布不均的數據集。而本文方法基本不受空間分布的影響，對不同數據集均能表現出穩定的加速效果。

圖4不同數據劃分法性能對比結果圖Fig.4 Experimental results of differentmethods on execution time，speedup ratio and load balancing

圖4 (c)描述了對不同數據復雜度數據集的測試結果。對于復雜度較低的數據5，不同方法表現出較好的性能。對于復雜度較高的數據6，傳統的兩種方法均不能適用，具體表現在:對于相同的串行時間468.95 s，傳統方法的最少運行時間分別為64.59 s和54.27 s，最大加速比為7.26和8.64；兩者的負載均衡指數均較高，最小值為3.34。主要原因在于數據6多邊形節點數相差巨大，存在多個包含大量節點數的大多邊形，上述兩種劃分方法均無法有效分配大多邊形，導致數據嚴重傾斜，從而產生并行處理過程中的負載不均衡問題。比較而言，本文方法可以較好地處理復雜度較高的數據集，獲得良好的加速比(16.15)。

總結來說，本文提出的數據劃分方法較傳統方法能更穩定、更有效地處理不同類型的矢量多邊形數據集，獲得更高的并行加速比、更好的負載均衡性能。

3.4 不同度量標準對運行時間的影響

本文提出的數據劃分方法中可采用節點數或多邊形數作為度量分塊計算量的標準。并行算法采用的度量標準不同，其并行執行效率也不同。實驗將并行算法執行總時間分為預處理、并行執行、后處理和I/O時間；分別采用多邊形節點數和多邊形數為度量標準，執行并行多邊形柵格化算法計算數據1和數據2，并統計其運行時間(見表2)。

在總時間上，以多邊形節點數為度量標準的并行算法運行時間少于以多邊形數為度量標準的并行算法，其中以大數據量的數據1更為明顯:兩者執行數據1的最少時間分別為104.54s和124.38s；執行數據2的最少時間分別為50.26s和49.49s。這表明針對大數據量的矢量多邊形數據集，以多邊形節點數為度量標準比以多邊形數為度量標準的并行算法可獲得更好的加速效果。

在各部分執行時間上，兩者的后處理時間和I/O時間相差不大，而在預處理時間和并行執行時間上相差較大，表現在:以多邊形節點數為度量標準的并行算法預處理時間高于以多邊形數為標準的算法，而其并行執行時間明顯較少。原因在于:以多邊形節點數為度量標準的并行算法在數據劃分中需要重復統計分塊中的多邊形節點數，因而耗時長；但其時間收益大于以多邊形數為標準的并行算法，主要表現在其并行執行時間明顯降低較快。這表明針對大數據量的多邊形數據集，以多邊形節點數為標準可更精確地度量計算量，從而使得各劃分分塊計算量更加均衡，獲得更高的并行執行效率。

表2 不同度量標準對并行效率影響的測試結果Tab.2 Results of the influence of differentmetrics on parallel efficiency s

4 結論

提出一種新的數據劃分方法，主要包括:將多邊形節點數或多邊形數作為度量計算量的標準；迭代計算劃分線的位置，在每次迭代中保證分塊間的計算量大致均衡，完成數據劃分；提出了基于二叉樹的劃分結果融合策略，解決了跨邊界多邊形的快速融合問題。在多核CPU環境下實現并行算法；選用多個土地利用現狀數據集進行測試。結果表明:本文提出的數據劃分方法較傳統方法針對不同類型多邊形數據集可獲得更高的并行效率、更好的負載均衡性能；針對大數據量的多邊形數據集，以多邊形節點數為度量標準可更精確地估算分塊計算量，從而更好地實現負載均衡。

References)

［1］Goodchild M F.Scale in GIS:an overview［J］.Geomorphology，2011，130(1-2):5-9.

［2］吳立新，史文中.地理信息系統原理與算法［M］.北京:科學出版社，2003.WU Lixin，SHIWenzhong.Geographic information systems principles and algorithms［M］.Beijing:Science Press，2003.(in Chinese)

［3］Mineter M J.A software framework to create vector-topology in parallel GIS operations［J］.International Journal of Geographical Information Science，2003，17(3):203-222.

［4］劉軍志，朱阿興，劉永波，等.基于柵格分層的逐柵格匯流算法并行化研究［J］.國防科技大學學報，2013，35(1):123-129.LIU Junzhi，ZHU Axing，LIU Yongbo，etal.Parallelization of a grid-to-grid routing algorithm based on grids layering［J］.Journal of National University of Defense Technology，2013，35(1):123-129.(in Chinese)

［5］Yang CW，Goodchild M F，Huang Q Y，et al.Spatial cloud computing:how can the geospatial sciences use and help shape cloud computing?［J］.International Journal of Digital Earth，2011，4(4):305-329.

［6］吳立新，楊宜舟，秦承志，等.面向新型硬件構架的新一代GIS基礎并行算法研究［J］.地理與地理信息科學，2013，29(4):1-8.WU Lixin，YANG Yizhou，QIN Chengzhi，et al.On basic geographic parallel algorithms of new generation GIS for new hardware architectures［J］.Geography and Geo-Information Science，2013，29(4):1-8.(in Chinese)

［7］程果，景寧，陳犖，等.柵格數據處理中鄰域型算法的并行優化方法［J］.國防科技大學學報，2012，34(4): 114-119.CHENG Guo，JING Ning，CHEN Luo，et al.Parallel optimization methods for raster data processing algorithms of neighborhood-scope［J］.Journal of National University of Defense Technology，2012，34(4):114-119.(in Chinese)

［8］Shekhar S，Ravada S，Chubb D，et al.Declustering and load-balancing methods for parallelizing geographic information systems［J］.IEEE Transactions on Knowledge and Data Engineering，1998，10:632-655.

［9］Ferhatosmanoglu H，Tosun A S，Canahuate G，et al.Efficient parallel processing of range queries through replicated declustering［J］.Distributed and Parallel Databases，2006，20(2):117-147.

［10］Meng L K，Huang C Q，Zhao C Y，et al.An improved hilbert curve for parallel spatial data partitioning［J］.Geospatial Information Science，2007，10(4):282-286.

［11］Hawick K A，Coddington P D，James H A.Distributed frameworks and parallel algorithms for processing large-scale geographic data［J］.Parallel Computing，2003，29(10): 1297-1333.

［12］Ye J Y，Chen B，Chen J，et al.A spatial data partition algorithm based on statistical cluster［C］//Proceedings of the 19th International Conference on Geoinformatics，NY:IEEE，2011:1-6.

［13］范俊甫，馬廷，季民，等.GIS中8種圖層級多核并行多邊形疊置分析工具的實現及優化方法［J］.地理科學進展，2013，32(12):1835-1844.FAN Junfu，MA Ting，JIMin，et al.Implementation and optimization of eight parallel polygon overlapping tools with OpenMP at the feature layer level in GIS［J］.Progress in Geography，2013，32(12):1835-1844.(in Chinese)

［14］Wang SW，Cowles M K，Armstrong M P.Grid computing of spatial statistics:using the TeraGrid for Gi*(d)analysis［J］.Concurrency and Computation:Practice and Experience，2008，20(14):1697-1720.

［15］Armstrong M P，Pavlik CE，Marciano R.Experiments in the measurement of spatial association using a parallel supercomputer［J］.Geographical Systems，1994，1(4): 267-288.

［16］Mineter M J，Dowers S.Parallel processing for geographical applications:a layered approach［J］.Journal of Geographical Systems，1999，1(1):61-74.

［17］Agarwal D，Puri S，He X，et al.A system for GIS polygonal overlay computation on Linux cluster—an experience and performance report［C］//Proceedings of the 26th International Parallel and Distributed Processing Symposium Workshops＆PhD Forum，NY:IEEE，2012:1433-1439.

［18］Wang Y F，Chen Z J，Cheng L，et al.Parallel scanline algorithm for rapid rasterization of vector geographic data［J］.Computers＆Geosciences，2013，59:31-40.

［19］Lee CK，HamdiM.Parallel image processing applicationson a network of workstations［J］.Parallel Computing，1995，21(1):137-160.

［20］Liao S B，Bai Y.A new grid-cell-based method for error evaluation of vector-to-raster conversion［J］.Computational Geosciences，2010，14(4):539-549.

［21］Zhou C H，Ou Y，Yang L，et al.An equal area conversion model for rasterization of vector polygons［J］.Science in China Series D:Earth Science，50(1):169-175.

A novel data decom position method for rapid parallel processing of vector polygon rasterization

ZHOU Chen，LIManchun，CHEN Zhenjie，JIANG Penghui，CHEN Dong
(School of Geographic and Oceanographic Sciences，Nanjing University，Nanjing 210023，China)

According to the load balance problem of large-scale parallel vector polygon rasterization，a novel data decomposition method was proposed.Firstly，the number of polygon nodes or the number of polygons was employed to evaluate the amount of calculations of a subset.The spatial locations of decomposed lineswere computed iteratively and the balanced calculations between decomposed subsetswere guaranteed，so as to realize data decomposition and load balancing.Secondly，a binary-tree based fusion strategy was put forth to merge the polygons acrossmultiple subsets.The proposed parallel algorithm was implemented under amulti-core CPU-based environment and multiple China land use datasets were employed.Experimental results show that the presentedmethod can outperform conventionalmethods for different datasets and can achieve a higher speed-up ratio and good load balancing.Moreover，when dealing with a large-scale vector dataset，the number of polygonal nodes is more appropriate to be themetric to evaluate the calculation of a subset precisely.

geographical information system；parallel computing；vector polygon rasterization；data decomposition；load balancing

TP751

A

1001-2486(2015)05-021-08

10.11887/j.cn.201505004

http://journal.nudt.edu.cn

2015-06-16

國家863計劃資助項目(2011AA120301)

周琛(1990—)，男，江蘇宿遷人，博士研究生，E-mail:njuzhouc@gmail.com；李滿春(通信作者)，男，教授，博士，博士生導師，E-mail:limanchun@yahoo.com