邱強(qiáng)，曹磊，方金云

（1.中國(guó)科學(xué)院大學(xué)，北京 100049；2.中國(guó)科學(xué)院計(jì)算技術(shù)研究所，北京 100190）

矢量數(shù)據(jù)疊加分析是空間分析中最基礎(chǔ)、最重要的分析方法之一。大數(shù)據(jù)集的疊加分析耗時(shí)較長(zhǎng)，耗時(shí)與數(shù)據(jù)規(guī)模具有正相關(guān)性。然而，隨著數(shù)據(jù)采集手段的不斷更新和GIS應(yīng)用的持續(xù)深入，矢量數(shù)據(jù)呈現(xiàn)出海量化和處理實(shí)時(shí)化的要求，這對(duì)疊加分析在計(jì)算性能、處理能力等方面提出了新的要求［1］。并行計(jì)算將大的數(shù)據(jù)處理與計(jì)算分解為多個(gè)可相互獨(dú)立、同時(shí)進(jìn)行的子任務(wù)，并通過(guò)這些子任務(wù)相互協(xié)調(diào)的運(yùn)行，實(shí)現(xiàn)快速、高效的問(wèn)題求解［2］。點(diǎn)面疊加分析存在數(shù)據(jù)密集、計(jì)算量大等特點(diǎn)。由于點(diǎn)數(shù)據(jù)之間不存在數(shù)據(jù)加鎖，利于劃分，并且其歸并過(guò)程不會(huì)帶來(lái)冗余計(jì)算，因此點(diǎn)面疊加分析適用于并行計(jì)算環(huán)境。本文利用MPI消息傳遞模型，基于Linux集群實(shí)現(xiàn)了高效、魯棒的并行點(diǎn)面疊加算法。

1 點(diǎn)的包含性測(cè)試

點(diǎn)的包含性測(cè)試是計(jì)算幾何的基本問(wèn)題，其含義是判斷一個(gè)點(diǎn)是否位于多邊形的內(nèi)部。隨著地理信息系統(tǒng)的發(fā)展，點(diǎn)與多邊形的包含性測(cè)試成為空間數(shù)據(jù)處理的一個(gè)基本問(wèn)題。研究人員提出了一系列點(diǎn)包含性測(cè)試算法。一些算法僅僅適用于特殊的多邊形，如文獻(xiàn)［3］中的方法僅適用于凸多邊形，文獻(xiàn)［4］中的網(wǎng)格法適合于柵格數(shù)據(jù)。另外一些算法如射線(xiàn)法、角度之和法、Winding Number法以及基于三角形的方法等［5，6］適用于任何多邊形。

這些方法對(duì)單個(gè)點(diǎn)的包含性測(cè)試是高效的，但并不適合處理批量點(diǎn)。對(duì)大量點(diǎn)進(jìn)行包含性測(cè)試最有效的方法是基于Cell的方法［4］和基于多邊形凸剖分的方法［7］?；贑ell的方法將多邊形柵格化，然后判斷不含多邊形邊的柵格在多邊形內(nèi)部還是外部。對(duì)點(diǎn)進(jìn)行包含性測(cè)試時(shí)，首先計(jì)算點(diǎn)在哪個(gè)柵格內(nèi)，若所在柵格不含任何邊，則可以在O（1）時(shí)間內(nèi)得到點(diǎn)與多邊形的關(guān)系；否則需要計(jì)算點(diǎn)與當(dāng)前柵格所包含的邊的關(guān)系。此方法對(duì)部分點(diǎn)的包含性測(cè)試時(shí)間復(fù)雜度是O（1），因此對(duì)于簡(jiǎn)單多邊形是高效的。基于多邊形凸剖分的方法是將多邊形進(jìn)行凸剖分，并用二叉空間分割（Binary Space Partition，BSP）樹(shù)管理這些凸剖分，對(duì)任何一個(gè)點(diǎn)，查詢(xún)BSP樹(shù)，計(jì)算可能包含該點(diǎn)的凸剖分，然后通過(guò)點(diǎn)與凸剖分的關(guān)系得到點(diǎn)與多邊形的關(guān)系。其包含性測(cè)試非常迅速，因此適合點(diǎn)集大的情形，該算法缺點(diǎn)是預(yù)處理時(shí)間較多，尤其對(duì)于邊數(shù)較多的復(fù)雜多邊形。

本文利用基于平均條帶劃分的點(diǎn)包含性測(cè)試方法［8］作為實(shí)現(xiàn)并行點(diǎn)面疊加算法的基礎(chǔ)。該方法首先對(duì)多邊形進(jìn)行簡(jiǎn)單預(yù)處理，然后采用射線(xiàn)法對(duì)每個(gè)點(diǎn)進(jìn)行包含性測(cè)試。預(yù)處理時(shí)，將多邊形的最小包圍矩形劃分為相同大小的豎直條帶，并計(jì)算每條邊在哪些條帶中，該預(yù)處理過(guò)程簡(jiǎn)單快速。利用預(yù)處理信息對(duì)點(diǎn)進(jìn)行包含性測(cè)試時(shí)，只需用當(dāng)前點(diǎn)構(gòu)造的射線(xiàn)和點(diǎn)所在的條帶包含的邊進(jìn)行相交測(cè)試即可。每個(gè)條帶內(nèi)包含的邊的數(shù)目有限，常常接近一個(gè)常數(shù)，遠(yuǎn)遠(yuǎn)小于多邊形的邊數(shù)，因此相交測(cè)試計(jì)算次數(shù)顯著減少，從而可以快速得到點(diǎn)與多邊形的關(guān)系。

2 基于MPI的并行算法

本文實(shí)現(xiàn)了基于靜態(tài)調(diào)度和基于動(dòng)態(tài)調(diào)度兩個(gè)版本，其差別在于數(shù)據(jù)劃分和調(diào)度策略上。主要包括輸入圖層的讀取、計(jì)算單元的劃分、子節(jié)點(diǎn)運(yùn)行點(diǎn)面疊加運(yùn)算、結(jié)果生成4個(gè)步驟［9］。

2.1 動(dòng)態(tài)調(diào)度

圖1為基于主、從結(jié)構(gòu)的點(diǎn)面疊加動(dòng)態(tài)調(diào)度體系結(jié)構(gòu)。動(dòng)態(tài)調(diào)度流程如圖2所示：首先，主節(jié)點(diǎn)從文件系統(tǒng)中讀取擬求交的點(diǎn)、面圖層（ShapeFile文件），完成文件讀入后主節(jié)點(diǎn)開(kāi)始對(duì)面要素圖層建立R－tree［10］索引，并為每個(gè)面要素建立一個(gè)潛在的相交集合，用以存儲(chǔ)與其可能相交的點(diǎn)；然后，主節(jié)點(diǎn)開(kāi)始遍歷每個(gè)點(diǎn)要素，如果點(diǎn)要素在某個(gè)面要素的MBR內(nèi)部，則將點(diǎn)加入到對(duì)應(yīng)面要素的集合中，作為一個(gè)基本的計(jì)算單元。算法1的偽代碼描述了基于R－tree的數(shù)據(jù)劃分過(guò)程。

圖1 動(dòng)態(tài)調(diào)度體系結(jié)構(gòu)Fig.1 Architecture with dynamic load strategy

預(yù)處理過(guò)程完成后，主節(jié)點(diǎn)將所有計(jì)算單元加入到任務(wù)隊(duì)列中并開(kāi)始調(diào)度，如圖3所示。開(kāi)始主節(jié)點(diǎn)順次向每個(gè)子節(jié)點(diǎn)發(fā)送一個(gè)計(jì)算單元，子節(jié)點(diǎn)接收計(jì)算單元后，運(yùn)行串行的點(diǎn)面疊加算法，完成計(jì)算后子節(jié)點(diǎn)暫時(shí)保存計(jì)算結(jié)果，并向主節(jié)點(diǎn)繼續(xù)請(qǐng)求計(jì)算單元；主節(jié)點(diǎn)接收到請(qǐng)求后，查看調(diào)度隊(duì)列中是否還有計(jì)算單元，如果有，就順次取出下一個(gè)向該子節(jié)點(diǎn)發(fā)送，反之就向子節(jié)點(diǎn)發(fā)送完成消息；子節(jié)點(diǎn)接收到完成消息后，將計(jì)算結(jié)果全部發(fā)送給主節(jié)點(diǎn)，主節(jié)點(diǎn)接收到所有子節(jié)點(diǎn)的計(jì)算結(jié)果則表明調(diào)度完成；最后主節(jié)點(diǎn)將計(jì)算結(jié)果一并寫(xiě)入文件系統(tǒng)。

圖2 動(dòng)態(tài)調(diào)度主、子節(jié)點(diǎn)流程Fig.2 Master／slave node process of dynamic load strategy

圖3 動(dòng)態(tài)調(diào)度過(guò)程Fig.3 Process of dynamic load strategy

2.2 靜態(tài)調(diào)度

相比于動(dòng)態(tài)調(diào)度，靜態(tài)調(diào)度比較直接和簡(jiǎn)單，如圖4所示。首先，主節(jié)點(diǎn)從文件系統(tǒng)中讀取要求交的點(diǎn)、面圖層；然后主節(jié)點(diǎn)根據(jù)子節(jié)點(diǎn)個(gè)數(shù)N，將面圖層根據(jù)要素ID分成N份，算法2描述了基于要素ID的數(shù)據(jù)劃分過(guò)程；之后主節(jié)點(diǎn)將相應(yīng)的面要素以及所有的點(diǎn)數(shù)據(jù)發(fā)送到相應(yīng)的子節(jié)點(diǎn)，即子節(jié)點(diǎn)i接收數(shù)據(jù)Di后，運(yùn)行串行的點(diǎn)面疊加算法，一旦完成計(jì)算，子節(jié)點(diǎn)便將結(jié)果發(fā)送給主節(jié)點(diǎn)；主節(jié)點(diǎn)接收到所有的子節(jié)點(diǎn)的數(shù)據(jù)后，便將計(jì)算結(jié)果寫(xiě)入到文件系統(tǒng)。

圖4 靜態(tài)調(diào)度主、子節(jié)點(diǎn)流程Fig.4 Master／slave node process of static load strategy

2.3 動(dòng)態(tài)調(diào)度與靜態(tài)調(diào)度對(duì)比

2.3.1 數(shù)據(jù)劃分 如算法1所示，在動(dòng)態(tài)調(diào)度中，主節(jié)點(diǎn)對(duì)面要素圖層建立R－tree索引，并為每一個(gè)面要素建立一個(gè)潛在的相交集合，用以存儲(chǔ)與其可能相交的點(diǎn)。索引建立完成后，主節(jié)點(diǎn)開(kāi)始遍歷每一個(gè)點(diǎn)要素，如果點(diǎn)要素在某一個(gè)面要素的MBR內(nèi)部，則將點(diǎn)加入到對(duì)應(yīng)面要素的集合中，作為一個(gè)基本的計(jì)算單元。

在靜態(tài)調(diào)度中，如算法2所示，主節(jié)點(diǎn)根據(jù)子節(jié)點(diǎn)個(gè)數(shù)N，將面圖層根據(jù)要素ID分成N份，然后將相應(yīng)的面要素以及所有的點(diǎn)數(shù)據(jù)作為一份發(fā)送到相應(yīng)的子節(jié)點(diǎn)。按ID劃分的方法相對(duì)粗略，不能表達(dá)出劃分對(duì)象的空間特性。當(dāng)空間數(shù)據(jù)分布不均衡時(shí)，容易造成劃分出來(lái)的各集合中空間對(duì)象分散，鄰近性差；雖然各個(gè)子任務(wù)的要素個(gè)數(shù)相等，但由于各個(gè)要素大小、復(fù)雜程度不同，容易導(dǎo)致數(shù)量?jī)A斜嚴(yán)重，降低并行運(yùn)算效率。

2.3.2 調(diào)度策略 動(dòng)態(tài)調(diào)度中，如圖3所示，主節(jié)點(diǎn)維護(hù)一個(gè)任務(wù)隊(duì)列，并動(dòng)態(tài)的根據(jù)子節(jié)點(diǎn)任務(wù)的完成情況，不斷地向空閑的子節(jié)點(diǎn)發(fā)送計(jì)算單元，直到任務(wù)隊(duì)列為空；靜態(tài)調(diào)度中，主節(jié)點(diǎn)只在一開(kāi)始向子節(jié)點(diǎn)發(fā)送數(shù)據(jù)，一旦發(fā)送完成，主節(jié)點(diǎn)就等待接收子節(jié)點(diǎn)的結(jié)果數(shù)據(jù)。

3 實(shí)驗(yàn)與分析

3.1 實(shí)驗(yàn)環(huán)境與數(shù)據(jù)

實(shí)驗(yàn)所采用的硬件環(huán)境為：CPU Intel Xeon，主頻3.4GHz，共8核；內(nèi)存4GB DDR3 1 333MHz ECC；外存500GB 7200RPM SATA硬盤(pán)；采用千兆以太網(wǎng)連接各個(gè)節(jié)點(diǎn)。軟件環(huán)境：操作系統(tǒng)為64位CentOS專(zhuān)業(yè)版；MPI采用OpenMPI版本；集成開(kāi)發(fā)環(huán)境為Eclipse＋CDT；采用GCC作為編譯器。

實(shí)驗(yàn)所用的點(diǎn)數(shù)據(jù)是全國(guó)居民點(diǎn)數(shù)據(jù)，點(diǎn)要素的個(gè)數(shù)為649 592個(gè)，對(duì)應(yīng)的Shapefile文件大小為139MB；面數(shù)據(jù)是全國(guó)縣界，面要素的個(gè)數(shù)為2 449個(gè)，對(duì)應(yīng)的Shapefile文件大小為36.8MB。

3.2 動(dòng)態(tài)調(diào)度與靜態(tài)調(diào)度對(duì)比

圖5顯示了基于動(dòng)態(tài)調(diào)度的并行點(diǎn)面疊加算法的實(shí)驗(yàn)結(jié)果，節(jié)點(diǎn)數(shù)分別是1、2、4、6、8。從總時(shí)間、I／O時(shí)間、計(jì)算時(shí)間3個(gè)維度進(jìn)行對(duì)比。其中，I／O時(shí)間包括主節(jié)點(diǎn)從文件系統(tǒng)讀入點(diǎn)、面數(shù)據(jù)以及主節(jié)點(diǎn)將最終的計(jì)算結(jié)果寫(xiě)入到文件系統(tǒng)中的時(shí)間；計(jì)算時(shí)間包括主節(jié)點(diǎn)的預(yù)處理時(shí)間、通信時(shí)間及子節(jié)點(diǎn)的運(yùn)算時(shí)間；總時(shí)間是I／O時(shí)間與計(jì)算時(shí)間之和。從圖5中可以看出，節(jié)點(diǎn)數(shù)從1增加到8，總時(shí)間下降顯著；I／O時(shí)間隨著節(jié)點(diǎn)數(shù)的增加，基本保持不變；而計(jì)算時(shí)間與總時(shí)間變化趨勢(shì)類(lèi)似，即節(jié)點(diǎn)數(shù)從1增加到8，計(jì)算時(shí)間明顯降低。

圖5 動(dòng)態(tài)調(diào)度Fig.5 Dynamic load strategy

圖6顯示了基于靜態(tài)調(diào)度的并行點(diǎn)面疊加算法的實(shí)驗(yàn)結(jié)果。總時(shí)間、I／O時(shí)間以及計(jì)算時(shí)間的變化趨勢(shì)與動(dòng)態(tài)調(diào)度的趨勢(shì)類(lèi)似。

圖6 靜態(tài)調(diào)度Fig.6 Static load strategy

圖7、圖8、圖9分別從總時(shí)間、計(jì)算時(shí)間和I／O時(shí)間對(duì)動(dòng)態(tài)調(diào)度和靜態(tài)調(diào)度做了對(duì)比。從總時(shí)間看，除了1、2個(gè)節(jié)點(diǎn)外，動(dòng)態(tài)調(diào)度都優(yōu)于靜態(tài)調(diào)度，但優(yōu)勢(shì)并不明顯。究其原因，一方面是由于動(dòng)態(tài)調(diào)度的預(yù)處理中構(gòu)建面要素的計(jì)算單元要占用一部分時(shí)間，而靜態(tài)調(diào)度幾乎不需要預(yù)處理，直接根據(jù)要素ID順次分配給子節(jié)點(diǎn)即可；另一方面，動(dòng)態(tài)調(diào)度每次給子節(jié)點(diǎn)發(fā)送一個(gè)面要素以及可能相交的點(diǎn)要素，而靜態(tài)調(diào)度一次性把子節(jié)點(diǎn)所需的所有數(shù)據(jù)發(fā)送完，因此相比于靜態(tài)調(diào)度，動(dòng)態(tài)調(diào)度增加了通信開(kāi)銷(xiāo)。由于上述兩個(gè)原因，雖然動(dòng)態(tài)調(diào)度提高了各個(gè)節(jié)點(diǎn)間的負(fù)載均衡，但相比于靜態(tài)調(diào)度，性能提升不明顯。對(duì)于I／O時(shí)間，由于要素的讀入與寫(xiě)回都是串行的，因此隨著節(jié)點(diǎn)數(shù)的增加，I／O時(shí)間并未減少，這也成為并行點(diǎn)面疊加算法性能提升的瓶頸。

圖7 總時(shí)間Fig.7 Total time

圖8 計(jì)算時(shí)間Fig.8 Computing time

圖9 I／O時(shí)間Fig.9 I／O time

4 結(jié)論與展望

本文利用MPI實(shí)現(xiàn)了基于Linux集群的并行點(diǎn)面疊加算法，并對(duì)比了動(dòng)態(tài)調(diào)度和靜態(tài)調(diào)度兩種調(diào)度策略。實(shí)驗(yàn)表明，動(dòng)態(tài)調(diào)度策略總體優(yōu)于靜態(tài)調(diào)度，但由于動(dòng)態(tài)調(diào)度增加了預(yù)處理時(shí)間和通信時(shí)間，但沒(méi)有減少I(mǎi)／O時(shí)間，因此優(yōu)勢(shì)不明顯。

由于動(dòng)態(tài)調(diào)度每次只給子節(jié)點(diǎn)發(fā)送一個(gè)面要素以及可能相交的點(diǎn)要素，因此增加了通信開(kāi)銷(xiāo)。這種開(kāi)銷(xiāo)在多機(jī)多核和集群下的增加更為明顯。下一步考慮增大數(shù)據(jù)劃分的粒度，即在動(dòng)態(tài)調(diào)度中，主節(jié)點(diǎn)一次向子節(jié)點(diǎn)傳輸多個(gè)面要素，從而減少通信開(kāi)銷(xiāo)。

由于I／O不可并行，因此隨著節(jié)點(diǎn)數(shù)的增加，I／O時(shí)間并沒(méi)有減少，這成為并行點(diǎn)面疊加算法性能提升的瓶頸。利用并行文件系統(tǒng)減少I(mǎi)／O開(kāi)銷(xiāo)將作為未來(lái)工作的重點(diǎn)。

［1］ 王結(jié)臣，王豹，胡瑋.并行空間分析算法研究進(jìn)展及評(píng)述［J］.地理與地理信息科學(xué)，2011，27（6）：1－5.

［2］ 陳國(guó)良.并行計(jì)算：結(jié)構(gòu)·算法·編程［M］.北京：高等教育出版社，2003.4－7.

［3］ HUANG C，SHIH T.On the complexity of point－in－polygon algorithms［J］.Computers ＆ Geosciences，1997，23（1）：109－118.

［4］ ZALIC B，KOLINGEROVA I.A cell－based point－in－polygon algorithm suitable for large sets of points［J］.Computers ＆ Geosciences，2001，27（10）：1135－1145.

［5］ HORMANN K，AGATHOS A.The point in polygon problem for arbitrary polyogns［J］.Computational Geometry：Theory and Applications，2001，20（3）：131－144.

［6］ FEITO F，TERRES J C，URENA A.Orientation，simplicity，and inclusion test for planar polygons［J］.Computers ＆ Graphics，1995，19（4）：595－600.

［7］ LI J，WANG W，WU E.Point－in－polygon tests by convex decomposition［J］.Computer ＆ Graphics，2007，31（4）：636－648.

［8］ 朱效民，劉焱，廖浩均.基于平均條帶劃分的點(diǎn)的包含性測(cè)試方法［J］.高技術(shù)通訊，2011，21（8）：810－816

［9］ AGARWAL D，PURI S，HE X，et al.A system for GIS polygonal overlay computation on linux cluster－an experience and performance report［A］.Parallel and Distributed Processing Symposium Workshops ＆PhD Forum （IPDPSW），2012IEEE 26th International IEEE［C］.2012.1433－1439.

［10］ GUTTMAN A.R－trees：A dynamic index structure for spatial searching［J］.ACM，1984，14（2）：47－57.