謝傳節(jié)，龍　舟，2，馬益杭，2，由志杰，2

1. 中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100101； 2. 中國科學(xué)院大學(xué)，北京 100049

Parallel Algorithm of Spatial Relationship Query between Polygons Based on Heterogeneous Multi-core Architecture

XIE Chuanjie1,LONG Zhou1，2,MA Yihang1，2,YOU Zhijie1，2

1. State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing 100101,China; 2.University of Chinese Academy of Sciences, Beijing 100049

?

多邊形間空間關(guān)系查詢的異構(gòu)多核架構(gòu)并行算法

謝傳節(jié)1，龍舟1，2，馬益杭1，2，由志杰1，2

1. 中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100101； 2. 中國科學(xué)院大學(xué)，北京 100049

Foundation support： The National High-tech Research and Development Program of China (863 Program) (Nos.2011AA120302; 2011AA120306; 2012AA12A401); The Marine Public Welfare Project (No.201105033-6)

摘要：目前在空間關(guān)系查詢中常用的Plane Sweep算法是一種串行算法，在處理海量空間數(shù)據(jù)時(shí)效率較低，而已有的并行計(jì)算方法對于普通的計(jì)算機(jī)并不適用。本文針對這個問題，提出了一種多邊形間空間關(guān)系查詢的異構(gòu)多核架構(gòu)并行算法，該算法先利用STR樹索引過濾掉不相交的多邊形，然后將過濾后的多邊形數(shù)據(jù)集合分解為點(diǎn)集合和邊集合，并對其構(gòu)建四叉樹索引；在保證數(shù)據(jù)浮點(diǎn)運(yùn)算精度符合要求的情況下，利用GPU強(qiáng)大的批量運(yùn)算能力快速處理邊與邊的相交情況并據(jù)此逐步計(jì)算得到環(huán)間的拓?fù)潢P(guān)系，再根據(jù)環(huán)間拓?fù)潢P(guān)系計(jì)算得到多邊形間的維度擴(kuò)展九交模型(DE-9IM)參數(shù)值；根據(jù)DE-9IM參數(shù)值與空間關(guān)系查詢條件相比對，輸出查詢結(jié)果。最后通過試驗(yàn)驗(yàn)證了算法的準(zhǔn)確性與高效性。

關(guān)鍵詞：異構(gòu)多核；并行計(jì)算；拓?fù)潢P(guān)系；空間關(guān)系查詢

空間關(guān)系查詢指的是從空間數(shù)據(jù)集中查找滿足某種空間關(guān)系的空間目標(biāo)的過程，它廣泛應(yīng)用于海量數(shù)據(jù)集或海量數(shù)據(jù)庫操作中，是地學(xué)計(jì)算中必不可少且十分常用的一項(xiàng)數(shù)據(jù)操作。隨著空間信息獲取技術(shù)日趨成熟，所獲取的空間數(shù)據(jù)量急速增加，如何將這些海量、超海量空間數(shù)據(jù)快速地處理并應(yīng)用于地學(xué)計(jì)算，從而獲得更有價(jià)值的信息，已經(jīng)成為現(xiàn)今地理信息系統(tǒng)技術(shù)創(chuàng)新的熱點(diǎn)[1]。

空間關(guān)系一般包括拓?fù)潢P(guān)系、度量關(guān)系、方位關(guān)系等，而空間拓?fù)潢P(guān)系是空間關(guān)系的主要研究內(nèi)容[2]。如今，在諸多開源的GIS代碼庫中，大部分判斷空間拓?fù)潢P(guān)系的算法都是基于Plane Sweep算法及其衍生算法編寫的。文獻(xiàn)[3]最早提出了Plane Sweep算法，該算法只能檢測到是否存在相交的線段，但不能查詢到所有相交的線段，文獻(xiàn)[4]在此之后對算法進(jìn)行了完善，解決了這個問題。文獻(xiàn)[5—7]改進(jìn)了Plane Sweep算法，通過降低算法的時(shí)間復(fù)雜度來提高計(jì)算效率；文獻(xiàn)[8]在算法中加入了方位角判斷和坐標(biāo)比較，解決了算法對共線線段及3條以上線段交于一點(diǎn)不適用的情況。但因?yàn)樵撍惴ū旧硎且环N串行算法，所以在針對海量數(shù)據(jù)的計(jì)算時(shí)效率會非常低；文獻(xiàn)[9]使用矢量叉積的原理判斷兩條線段的相交類型，這種方法計(jì)算準(zhǔn)確，但需要對線段兩兩逐個比較求交，效率過低；文獻(xiàn)[10—12]基于PRAM并行計(jì)算模型，使用不同的方法對一組不相交的線段集并行構(gòu)建二叉樹，然后通過遍歷和搜索樹，得到最終相交線段的結(jié)果。但由于空間關(guān)系查詢是一項(xiàng)比較常用的操作，而這種并行計(jì)算模型在普通計(jì)算機(jī)上并不適用，所以這種并行模式難以滿足需要。

由于計(jì)算機(jī)硬件系統(tǒng)的不斷升級換代，基于異構(gòu)多核架構(gòu)的異構(gòu)計(jì)算已經(jīng)成為各高性能計(jì)算領(lǐng)域的主流技術(shù)。圖形處理單元(GPU)正廣泛應(yīng)用于大型數(shù)據(jù)集的并行處理中。與傳統(tǒng)的CPU多核技術(shù)相比較，GPU所擁有的密集型并行構(gòu)架具有更為強(qiáng)大的并行潛力以及更加顯著的加速效果[13-14]。

針對傳統(tǒng)算法在處理海量空間數(shù)據(jù)時(shí)存在的不足，本文打破原有的計(jì)算方法，設(shè)計(jì)了一并行計(jì)算方法，該方法先利用CPU+GPU架構(gòu)并行計(jì)算得到兩個圖層中所有最小外包矩形相交的多邊形的線段相交情況，然后根據(jù)線段相交情況判斷多邊形環(huán)間拓?fù)潢P(guān)系，再根據(jù)環(huán)間的拓?fù)潢P(guān)系統(tǒng)計(jì)多邊形間的維度擴(kuò)展九交模型(DE-9IM)參數(shù)值，從而得到多邊形間的空間關(guān)系。

1多邊形間空間關(guān)系查詢的并行算法

本文的空間關(guān)系查詢對象，是兩個海量的多邊形數(shù)據(jù)集合。兩個空間數(shù)據(jù)集分別作為目標(biāo)數(shù)據(jù)集合和源數(shù)據(jù)集合，其中，需要被統(tǒng)計(jì)的圖層作為目標(biāo)圖層，以空間關(guān)系條件進(jìn)行比對的圖層作為源圖層。該空間關(guān)系查詢計(jì)算流程(圖1)主要由線段相交分析、環(huán)間拓?fù)潢P(guān)系分析和多邊形間DE-9IM計(jì)算3個模塊組成。

圖1　多邊形間空間查詢計(jì)算方法Fig.1　Spatial query calculation method between polygons

1.1數(shù)據(jù)預(yù)處理

在對線段相交判斷之前，需要對兩個多邊形數(shù)據(jù)集進(jìn)行預(yù)處理，預(yù)處理過程分為以下兩個步驟。

第1步：幾何要素與圖形要素的轉(zhuǎn)換。將空間數(shù)據(jù)中的點(diǎn)(Point)轉(zhuǎn)化為圖形節(jié)點(diǎn)(Node)信息，空間數(shù)據(jù)中的線段(Segment)轉(zhuǎn)換為圖形邊(Edge)信息。邊中會記錄構(gòu)建該邊的兩個端點(diǎn)的節(jié)點(diǎn)信息，節(jié)點(diǎn)中也會記錄以該節(jié)點(diǎn)為某端點(diǎn)的所有邊信息。圖形要素中攜帶更多的標(biāo)識信息，其中很多信息是空間幾何要素?zé)o法表達(dá)和記錄的，這些信息便于位置關(guān)系的構(gòu)建和計(jì)算。

本文在進(jìn)行空間關(guān)系查詢時(shí)，首先遍歷兩個多邊形數(shù)據(jù)集合的所有多邊形要素，將所有的點(diǎn)和線段信息轉(zhuǎn)換并統(tǒng)一整合到同一個具有拓?fù)涮匦缘墓?jié)點(diǎn)和邊集合中。同時(shí)，獲取的空間信息還包括原始數(shù)據(jù)的最大空間范圍以及節(jié)點(diǎn)集中節(jié)點(diǎn)的個數(shù)。參數(shù)信息如圖2所示。

圖2　節(jié)點(diǎn)和邊的數(shù)據(jù)結(jié)構(gòu)信息Fig.2　The data structure information of the nodes and edges

第2步：構(gòu)建索引。為了減小計(jì)算量，首先需要先對兩個圖層的多邊形進(jìn)行過濾，根據(jù)多邊形的最小外包矩形過濾掉一定不相交的多邊形組合。本文選擇STR樹索引進(jìn)行過濾，STR樹索引是一種使用Sort-Tile-Recursive(STR)算法創(chuàng)建的只負(fù)責(zé)查詢處理的R樹。它針對二維空間數(shù)據(jù)，十分易于操作并且可以最大化地利用內(nèi)存空間存儲數(shù)據(jù)[15]。主要做法是：假設(shè)有兩個多邊形集合分別為A和B，先對A和B分別構(gòu)建STR樹，將多邊形的最小外包矩形作為STR樹節(jié)點(diǎn)的空間范圍，多邊形本身作為STR樹節(jié)點(diǎn)中存儲的空間對象，然后遍歷B中多邊形的最小外包矩形在A的STR樹中查詢，遍歷A中多邊形的最小外包矩形在B的STR樹中查詢，最終得到兩個過濾后的多邊形集合A和B。

對于過濾后的多邊形集合A和B，利用OpenMP并行技術(shù)[16]對它們的節(jié)點(diǎn)集合和邊集合構(gòu)建空間四叉樹索引[17-18]，并置于計(jì)算機(jī)內(nèi)存中，以便在樹節(jié)點(diǎn)中快速獲取邊邊組合，為下一步邊邊相交計(jì)算做準(zhǔn)備。如圖3(a)所示，假設(shè)一個多邊形為待構(gòu)建四叉樹索引的空間幾何體，完成空間四叉樹索引的步驟可分為以下兩步。

步驟a：提取該多邊形的所有節(jié)點(diǎn)構(gòu)成節(jié)點(diǎn)集，構(gòu)建點(diǎn)四叉樹。如圖3(b)所示，最大空間范圍決定了四叉樹根節(jié)點(diǎn)的空間范圍，每個葉子節(jié)點(diǎn)所容納結(jié)點(diǎn)的最小個數(shù)為1。由于各節(jié)點(diǎn)的分裂是完全獨(dú)立的，故可以利用OpenMP并行技術(shù)進(jìn)行新節(jié)點(diǎn)的分裂，直至各進(jìn)程內(nèi)的四叉樹都分裂完全。

步驟b：完善空間索引的邊信息。依照已構(gòu)建的空間點(diǎn)四叉樹索引，對于每一個葉子節(jié)點(diǎn)，提取其代表的空間范圍信息，依次與所有邊信息的最小外包矩形相比較，若兩個空間范圍有交集，則將此邊加入候選集合(先對它們的邊界進(jìn)行相交判斷是為了節(jié)省計(jì)算成本，進(jìn)行一次粗過濾)，然后再精確地判斷此邊與葉子節(jié)點(diǎn)所代表的空間范圍是否相交，如果相交，則被認(rèn)為屬于該葉子節(jié)點(diǎn)，最終完成邊信息的插入，構(gòu)建成一個包含了點(diǎn)信息和邊信息的空間四叉樹索引。如圖3(c)所示，這樣構(gòu)成的四叉樹索引可能會將一條邊存儲在多個葉子節(jié)點(diǎn)中，但這是為了在后續(xù)的計(jì)算中能保證每個葉子節(jié)點(diǎn)內(nèi)所有的相交線段都能被統(tǒng)計(jì)到，以保證最終計(jì)算結(jié)果的準(zhǔn)確性。

圖3　空間四叉樹索引建立過程Fig.3　The process of establishing the spatial quadtree index

1.2線段相交判斷

邊邊組合準(zhǔn)備好后，需要對其相交類型進(jìn)行判斷，在此之前，需要對邊邊組合進(jìn)行精度判斷。雖然GPU已經(jīng)支持雙精度浮點(diǎn)運(yùn)算，但與其單精度浮點(diǎn)運(yùn)算能力相比較，其雙精度浮點(diǎn)運(yùn)算能力還是偏弱，不如CPU的處理效率高。本文提取空間四叉樹索引中每一個葉子節(jié)點(diǎn)內(nèi)所有的邊邊組合(每一個邊邊組合的兩條邊來源于不同的數(shù)據(jù)集)，計(jì)算每一對邊邊組合中兩條邊的4個端點(diǎn)的相對坐標(biāo)，以此方法，大幅度降低端點(diǎn)坐標(biāo)的位數(shù)，使其有可能在計(jì)算誤差允許范圍內(nèi)參與單精度浮點(diǎn)運(yùn)算。根據(jù)單精度浮點(diǎn)計(jì)算的數(shù)值范圍，判斷相對坐標(biāo)是否滿足單精度浮點(diǎn)運(yùn)算所需的幾何計(jì)算精度要求[19-20]。若4個端點(diǎn)都滿足精度要求，則將該組存儲于滿足單精度浮點(diǎn)運(yùn)算的邊邊組合候選集中，將每一邊邊組合放入GPU的每一個線程中進(jìn)行相交判斷的并行計(jì)算。若4個端點(diǎn)中有任一端點(diǎn)不滿足，則將該邊邊組合直接置于CPU中進(jìn)行相交判斷和相交情況的計(jì)算。

本文使用CUDA編程模型進(jìn)行線段相交的并行計(jì)算，運(yùn)行在GPU上的程序稱為內(nèi)核函數(shù)(kernal)，CPU在調(diào)用內(nèi)核函數(shù)時(shí)，先通過邊邊組合的數(shù)量聲明內(nèi)核函數(shù)中的需要開啟的線程數(shù)量，每一個線程都有自己的block ID和thread ID用于與其他線程相區(qū)分。然后，CPU將數(shù)據(jù)從主存轉(zhuǎn)移到設(shè)備內(nèi)存，GPU將每一組數(shù)據(jù)分配給一個線程，由于每個線程只負(fù)責(zé)處理一組邊邊組合的相交判斷，所以各線程間無須通信，計(jì)算完成后直接將計(jì)算結(jié)果返回即可。

在進(jìn)行邊邊相交類型判斷時(shí)，利用向量混合積判斷兩條邊是否相交[21]，若判斷得到邊邊相交，如圖4所示，依據(jù)交點(diǎn)特征可將相交情況可分為普通相交、丁字相交、連接相交和疊置相交四大類：若交點(diǎn)不是相關(guān)兩條邊的4個端點(diǎn)中的任何一個，則為普通相交；若交點(diǎn)是相關(guān)兩條邊的4個端點(diǎn)中的一個，則為丁字相交；若交點(diǎn)是相關(guān)兩條邊的4個端點(diǎn)中的兩個，則為連接相交；若交點(diǎn)的相關(guān)兩條邊有公共重合部分，則為疊置相交。

圖4　4種線段相交的類型Fig.4　Four types of line segments intersect

1.3環(huán)與環(huán)之間拓?fù)潢P(guān)系判斷方法

一個多邊形是由多個環(huán)組成，一個環(huán)是由多條邊組成。欲求兩個多邊形之間的空間關(guān)系，可以先從環(huán)間的拓?fù)潢P(guān)系計(jì)算入手，而環(huán)間的拓?fù)潢P(guān)系需要根據(jù)線段相交類型進(jìn)行判斷。環(huán)間的拓?fù)潢P(guān)系分為以下8種：Disjoint、Meet、Equal、Overlap、Covers、Covered By、Contains和Inside，如圖5所示?；诰沤荒Ｐ屠碚?，環(huán)與環(huán)之間一定發(fā)生且只發(fā)生8種關(guān)系中的一種[22-23]。

假設(shè)待計(jì)算的兩個環(huán)為Ha、Hb，則根據(jù)線段相交類型判斷兩環(huán)關(guān)系的主要過程如下：①若兩環(huán)沒有相交線段出現(xiàn)，則關(guān)系一定為Disjoint、Inside或Contains中的一種，再通過兩環(huán)上兩點(diǎn)p、q與環(huán)間的位置關(guān)系具體確定；②若有普通相交出現(xiàn)，則兩環(huán)拓?fù)潢P(guān)系屬于Overlap；③若兩環(huán)的所有邊都發(fā)生了疊置相交，則兩環(huán)的拓?fù)潢P(guān)系為Equal；④若除了普通相交，其他相交類型均有出現(xiàn)，則關(guān)系一定為Covers、Covered By或Meet中的一種，再根據(jù)兩環(huán)位置關(guān)系具體判定；⑤由于九交模型只有8種情況，所以其他情況兩環(huán)的拓?fù)潢P(guān)系為Overlap。

圖5　環(huán)間的8種拓?fù)潢P(guān)系Fig.5　Eight kinds of topological relations between the rings

1.4多邊形間空間關(guān)系的判斷方法

在1.1節(jié)中由STR樹索引過濾后得到兩個多邊形集合，需要依次對有可能相交的多邊形組合判斷它們之間的空間關(guān)系，為了進(jìn)一步提高效率，本文使用BOOST庫多線程并行判斷多邊形間的空間關(guān)系。并行判斷涉及數(shù)據(jù)劃分，由于本文計(jì)算量最大的線段相交的判斷已經(jīng)在最開始由CPU和GPU協(xié)作完成，則對于每一對多邊形組合，只需調(diào)用計(jì)算好的線段相交類型來判斷環(huán)間拓?fù)潢P(guān)系，進(jìn)而判斷它們的空間關(guān)系，所以在計(jì)算量上不會有很大差距。那么對多邊形數(shù)據(jù)的劃分，可以按照有可能相交的目標(biāo)多邊形按數(shù)量平均分組，以保證每個線程的負(fù)載基本相同，然后將每一組多邊形組合分配給一個線程進(jìn)行空間關(guān)系的判斷，最后合并計(jì)算的結(jié)果。

多邊形之間的空間關(guān)系由維度擴(kuò)展九交模型(DE-9IM)來表達(dá)，根據(jù)維度擴(kuò)展九交模型的參數(shù)值，比較空間關(guān)系的查詢條件，即可輸出滿足條件的多邊形集合。根據(jù)OGC的標(biāo)準(zhǔn)，與DE-9IM參數(shù)值對應(yīng)的7類空間關(guān)系查詢條件為：分離關(guān)系、重疊關(guān)系、相等關(guān)系、接觸關(guān)系、包含關(guān)系、包含于關(guān)系以及相交關(guān)系(不屬于分離關(guān)系的所有空間關(guān)系都為相交關(guān)系)[24]。根據(jù)環(huán)間拓?fù)潢P(guān)系分析模塊，分別得到的兩個多邊形內(nèi)、外環(huán)之間的拓?fù)潢P(guān)系，利用此拓?fù)潢P(guān)系逐步判斷兩個多邊形間內(nèi)部、邊界和外部的相交值，從而確定多邊形間的DE-9IM參數(shù)值，最后將計(jì)算出的DE-9IM參數(shù)值與7類空間關(guān)系對應(yīng)的DE-9IM參數(shù)進(jìn)行比對，以最終確定多邊形間的空間關(guān)系。

2試驗(yàn)設(shè)計(jì)與結(jié)果分析

2.1試驗(yàn)環(huán)境

2.1.1硬件環(huán)境

CPU：Intel(R)Core(TM)i7-3770。

內(nèi)存：8.00 GB。

GPU：獨(dú)立顯卡NVIDIA GTX 660(Device2-Graphics)。

2.1.2軟件環(huán)境

操作系統(tǒng)：Win7 64位操作系統(tǒng)。

編程環(huán)境：Visual Studio 2012。

編程語言：C++、CUDA。

開源庫：GDAL庫(版本為1.1.0)、GEOS庫(版本為3.4.2)、BOOST庫(版本為1.53.0)、CUDA庫(版本為5.0)、OpenMP庫(版本為2.0)。

2.2驗(yàn)證并行算法正確性與穩(wěn)健性

任意兩個帶洞多邊形(即除了一個外環(huán)還可能有任意個內(nèi)環(huán)的多邊形)之間一定發(fā)生且只能發(fā)生18種拓?fù)潢P(guān)系中的一種[25]，這18種拓?fù)潢P(guān)系分別為：disjoint、fillingHole、holeFIlled、meet、equal、inside、holeCoveredBy、coveredBy、contains、holeCovers、covers、insideOverfill、containsOverfill、insideContains、coversCoveredBy、coversOverfill、CoveredByOverfill、overlap。

根據(jù)18種拓?fù)潢P(guān)系，利用ArcGIS軟件逐一繪制相應(yīng)的多邊形計(jì)算數(shù)據(jù)，如圖6所示。每一種拓?fù)潢P(guān)系都可能對應(yīng)于多種不同的多邊形模型，多邊形數(shù)據(jù)中紅色多邊形為目標(biāo)數(shù)據(jù)，黃色多邊形為源數(shù)據(jù)。這些模型中的多邊形樣式或形態(tài)不同，但其拓?fù)潢P(guān)系相同且唯一，保證了試驗(yàn)數(shù)據(jù)的穩(wěn)健性。

圖6　根據(jù)帶洞多邊形間18種拓?fù)潢P(guān)系繪制的多邊形數(shù)據(jù)Fig.6　Polygon data according to 18 kinds topological relations between the polygons with holes

對于不同的試驗(yàn)數(shù)據(jù)，測試的空間關(guān)系查詢條件也不相同，其主要目的是驗(yàn)證并行算法是否準(zhǔn)確，比對并行程序輸出的多邊形ID與ArcGIS中“按位置選擇”功能選擇出的多邊形ID是否相一致。

經(jīng)過并行算法與ArcGIS中“按位置選擇”功能對18種代表不同的拓?fù)潢P(guān)系的多邊形數(shù)據(jù)與其對應(yīng)的空間關(guān)系查詢條件進(jìn)行測試，分別得到的多邊形ID如表1所示。從表中可以看出，對于每一種多邊形數(shù)據(jù)，按照對應(yīng)的空間關(guān)系查詢條件進(jìn)行查詢，并行算法查詢到的多邊形ID與ArcGIS查詢到的多邊形ID都相同。

2.3驗(yàn)證并行算法加速效果

如圖7所示，本文試驗(yàn)數(shù)據(jù)為從1990年山東省土地利用圖中提取出的50 000個多邊形，將所有多邊形統(tǒng)一構(gòu)建成為一個海量多邊形圖層，并作為目標(biāo)圖層數(shù)據(jù)，同時(shí)對上述海量多邊形圖層的副本進(jìn)行整體或局部的修改，作為源圖層數(shù)據(jù)。

表1并行算法與ArcGIS的測試結(jié)果

Tab.1The test results of parallel algorithm and ArcGIS

序號拓?fù)潢P(guān)系對應(yīng)查詢條件并行算法查詢結(jié)果ArcGIS查詢結(jié)果1disjoint分離002fillingHole接觸003holeFIlled接觸004meet接觸01234012345equal相等006inside被包含01017holeCoveredBy被包含008coveredBy被包含01019contains包含0010holeCovers包含0011covers包含0012insideOverfill相交0013containsOverfill相交0014insideContains相交0015coversCoveredBy相交0016coversOverfill相交0017CoveredByOverfill相交0018overlap相交00

注：試驗(yàn)結(jié)果的正確與否參考ArcGIS中“按位置選擇”功能的計(jì)算結(jié)果。

圖7　試驗(yàn)數(shù)據(jù)Fig.7　The experiment data

為了驗(yàn)證并行算法的計(jì)算效率，本文利用ArcGIS組件式開發(fā)技術(shù)，制作了一套僅具有數(shù)據(jù)輸入輸出、“按位置選擇”和時(shí)間記錄3種功能的小型測試軟件。先記錄各個空間關(guān)系查詢條件下并行算法的測試時(shí)間，然后利用開源代碼庫GEOS中的空間關(guān)系計(jì)算測試時(shí)間，再利用小軟件記錄ArcGIS在各個空間關(guān)系查詢條件下的測試時(shí)間，最后對比3個測試的計(jì)算結(jié)果和測試時(shí)間。

本文把相等關(guān)系、相交關(guān)系、接觸關(guān)系、包含關(guān)系和被包含關(guān)系作為5個空間關(guān)系查詢條件進(jìn)行測試(由于相交即不分離，故不需要對分離關(guān)系再做測試)。對于相等關(guān)系的測試，本文保留原始數(shù)據(jù)的6000個多邊形位置不變，將其他多邊形進(jìn)行平移作為源圖層。對于其余4個關(guān)系的測試，本文將原始數(shù)據(jù)的副本做了全局性的平移，使其在位置上不同于原始數(shù)據(jù)，將新的數(shù)據(jù)作為源圖層。

將并行算法計(jì)算出的多邊形個數(shù)，與ArcGIS中“按位置選擇”功能相對應(yīng)的空間關(guān)系查詢條件計(jì)算得到的多邊形個數(shù)相對比，試驗(yàn)結(jié)果統(tǒng)計(jì)如表2所示，可以發(fā)現(xiàn)，并行算法的計(jì)算結(jié)果都是正確的。

表2并行算法計(jì)算與ArcGIS計(jì)算的試驗(yàn)結(jié)果統(tǒng)計(jì)

Tab.2Statistic the experiment results of parallel algorithm and ArcGIS

空間關(guān)系對應(yīng)于“按位置選擇”功能中的條件描述ArcGIS計(jì)算結(jié)果并行算法測試結(jié)果相等與源圖層要素完全相同60006000相交與源圖層要素相交62486248接觸接觸源圖層要素的邊界77包含包含源圖層要素1616被包含在源圖層要素范圍內(nèi)1111

本文對5種空間關(guān)系查詢的計(jì)算時(shí)間進(jìn)行了統(tǒng)計(jì)，包括數(shù)據(jù)預(yù)處理(構(gòu)建STR樹索引及幾何要素轉(zhuǎn)化為圖形要素)，構(gòu)建四叉樹索引，線段相交判斷(包括精度判斷時(shí)間)以及最后計(jì)算DE-9IM(包括環(huán)間拓?fù)潢P(guān)系計(jì)算)的時(shí)間，統(tǒng)計(jì)結(jié)果如表3所示。同時(shí)，本文也將并行算法的時(shí)間與ArcGIS的計(jì)算時(shí)間以及GEOS的計(jì)算時(shí)間進(jìn)行了統(tǒng)計(jì)對比，每個算法的測試時(shí)間對比如圖8所示。經(jīng)過對比可以發(fā)現(xiàn)，在5種關(guān)系查詢的時(shí)間測試中，并行算法的測試時(shí)間比GEOS和ArcGIS的計(jì)算時(shí)間都要短，在計(jì)算效率上具有較大優(yōu)勢。以相等關(guān)系查詢?yōu)槔篏EOS計(jì)算時(shí)間大于3500 s，ArcGIS計(jì)算時(shí)間為1 039.873 s，而并行算法總時(shí)間僅為401.847 s，比ArcGIS用時(shí)節(jié)省了638.026 s，計(jì)算效率是ArcGIS的2.59倍。

表3并行算法對5種空間關(guān)系查詢時(shí)間統(tǒng)計(jì)

Tab.3The query time of parallel algorithm for five kinds of spatial relations

s

圖8　針對5種空間關(guān)系查詢3種算法的時(shí)間對比Fig 8　Time comparison if three kinds of algorithm for five kinds of spatial relations query

由于本文是綜合利用CPU+GPU架構(gòu)的算法來進(jìn)行空間關(guān)系的查詢，為了驗(yàn)證本算法對線段相交判斷處理的優(yōu)勢，筆者以相交關(guān)系的判斷為例，統(tǒng)計(jì)所有線段相交均由CPU處理和均由GPU處理的時(shí)間進(jìn)行對比，對比結(jié)果如圖9所示。可以發(fā)現(xiàn)，綜合利用CPU和GPU對線段相交處理的速度比單獨(dú)用CPU或GPU對線段相交處理的速度都要快，說明綜合利用CPU和GPU進(jìn)行線段相交的判斷是優(yōu)于其他兩種方法的。

圖9　不同硬件處理線段相交的時(shí)間對比Fig.9　Time comparison of different hardware processingline segment intersection

2.4試驗(yàn)結(jié)果分析

針對以上試驗(yàn)結(jié)果，從正確性和高效性兩方面，對本文提出的空間關(guān)系并行算法的實(shí)際應(yīng)用進(jìn)行評價(jià)。

(1) 正確性：在小規(guī)模的多邊形數(shù)據(jù)情形下，依據(jù)18種帶洞多邊形的拓?fù)潢P(guān)系，構(gòu)建各式多邊形模型數(shù)據(jù)，保證了所有的多邊形之間的空間關(guān)系可以被涵蓋在測試數(shù)據(jù)中。18種多邊形數(shù)據(jù)經(jīng)過不同的查詢條件計(jì)算，與ArcGIS的“按位置選擇”功能的計(jì)算結(jié)果作比對，結(jié)果是完全正確的。當(dāng)數(shù)據(jù)擴(kuò)展到50 000個多邊形數(shù)據(jù)時(shí)，計(jì)算結(jié)果依然與ArcGIS的計(jì)算結(jié)果相同，證明本文提出的空間關(guān)系并行查詢算法，在對于各種數(shù)據(jù)規(guī)模時(shí)，計(jì)算結(jié)果都是正確的。

(2) 高效性：對于規(guī)模較大的數(shù)據(jù)，進(jìn)行空間關(guān)系查詢是非常費(fèi)時(shí)的過程。本算法綜合利用CPU和GPU并行判斷線段相交情況，突破了Plane Sweep算法逐一掃描判斷的特性；在計(jì)算多邊形間空間關(guān)系時(shí)，利用CPU多線程并行判斷，因此大幅縮減了計(jì)算時(shí)間，經(jīng)過試驗(yàn)驗(yàn)證，本文提出的空間關(guān)系查詢并行算法的計(jì)算速度是明顯快于GEOS中的空間關(guān)系算法和ArcGIS中按位置選擇功能，所以本算法具有比較顯著的高效性。

3結(jié)論

本文算法充分利用CPU和GPU的并行計(jì)算能力，提高了對海量數(shù)據(jù)空間關(guān)系查詢的處理效率。本文算法能夠保證在查詢結(jié)果準(zhǔn)確的情況下高效地處理海量數(shù)據(jù)的空間關(guān)系查詢問題，在目標(biāo)數(shù)據(jù)和源數(shù)據(jù)都為50 000個多邊形的情況下，以相等關(guān)系為例，并行算法比ArcGIS計(jì)算時(shí)間節(jié)省了638.026 s，計(jì)算效率約為ArcGIS的2.6倍，通過其他關(guān)系的查詢時(shí)間對比也可以看出，并行算法能有效地提升查詢效率。由于本算法是根據(jù)多核CPU和GPU的硬件架構(gòu)而設(shè)計(jì)的通用算法，因此可以預(yù)測在其他擁有多核CPU和GPU的硬件環(huán)境下，本算法也可以同樣提升計(jì)算效率。

參考文獻(xiàn)：

[1]王結(jié)臣, 王豹, 胡瑋, 等. 并行空間分析算法研究進(jìn)展及評述[J]. 地理與地理信息科學(xué), 2011, 27(6): 1-5.

WANG Jiechen, WANG Bao, HU wei, et al. Review on Parallel Spatial Analysis Algorithms[J]. Geography and Geo-Information Science, 2011, 27(6): 1-5.

[2]陳軍, 趙仁亮. GIS空間關(guān)系的基本問題與研究進(jìn)展[J]. 測繪學(xué)報(bào), 1999, 28(2): 95-102.

CHEN Jun, ZHAO Renliang. Spatial Relations in GIS: A Survey on Its Key Issues and Research Progress[J]. Acta Geodaetica et Cartographica Sinica, 1999, 28(2): 95-102.

[3]SHAMOS M I, HOEY D. Geometric Intersection Problems[C]∥17th Annual Symposium on Foundations of Computer Science, 1976. Houston, TX: IEEE, 1976: 208-215.

[4]BENTLEY J L, OTTMANN T A. Algorithms for Reporting and Counting Geometric Intersections[J]. IEEE Transactions on Computers, 1979, C-28(9): 643-647.

[5]BARTUSCHKA U, MEHLHORN K, NHER S. A Robust and Efficient Implementation of a Sweep Line Algorithm for the Straight Line Segment Intersection Problem[C]∥Proceedings of Workshop on Algorithm Engineering. 1997.

[6]CHAZELLE B, EDELSBRUNNER H. An Optimal Algorithm for Intersecting Line Segments in the Plane[J]. Journal of the ACM, 1992, 39(1): 1-54.

[7]BALABAN I J. An Optimal Algorithm for Finding Segments Intersections[C]∥Proceedings of the 11th Symposium on Computational Geometry. New York: ACM, 1995: 211-219.

[8]陳軍, 劉萬增, 李志林, 等. 線目標(biāo)間拓?fù)潢P(guān)系的細(xì)化計(jì)算方法[J]. 測繪學(xué)報(bào), 2006, 35(3): 255-260.

CHEN Jun, LIU Wanzeng, LI Zhilin, et al. The Refined Calculation Method of Topological Relationships between Line Objects[J]. Acta Geodaetica et Cartographica Sinica, 2006, 35(3): 255-260.

[9]陳斐, 周曉光. 基于平面分割的線/線相接、交叉類型判斷[J]. 測繪學(xué)報(bào), 2014, 43(2): 186-192.

CHEN Fei, ZHOU Xiaoguang. An Algorithm for Determining the Touching and Crossing Relations between a Pair of Lines Based on One Line Splitting Plane to Two Parts[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(2): 186-192.

[10]GOODRICH M T, GHOUSE M R, BRIGHT J. Sweep Methods for Parallel Computational Geometry[J]. Algorithmica, 1996, 15(2): 126-153.

[11]ATALLAH M J, GOODRICH M T. Efficient Plane Sweeping in Parallel[C]∥Proceedings of the Second Annual Symposium on Computational Geometry. New York: ACM, 1986: 216-225.

[12]AGGARWAL A, CHAZELLE B, GUIBAS L, et al. Parallel Computational Geometry[J]. Algorithmica, 1988, 3(1-4): 293-327.

[13]OWENS J D, LUEBKE D, GOVINDARAJU N, et al. A Survey of General-Purpose Computation on Graphics Hardware[J]. Computer Graphics Forum, 2007, 26(1): 80-113.

[14]張舒, 褚艷利. GPU高性能運(yùn)算之CUDA[M]. 北京: 中國水利水電出版社, 2009.

ZHANG Shu, CHU Yanli. CUDA GPU High-performance Computing[M]. Beijing: China Water Power Press, 2009.

[15]LEUTENEGGER S T, LOPEZ M A, EDGINGTON J. STR: A Simple and Efficient Algorithm for R-Tree Packing[C]∥Proceedings of the 13th International Conference on Data Engineering. Birmingham: IEEE, 1997: 497-506.

[16]周偉明. 多核計(jì)算與程序設(shè)計(jì)[M]. 武漢: 華中科技大學(xué)出版社, 2009.

ZHOU Weiming. Multi-core Computing and Programming[M]. Wuhan: Huazhong University of Science and Technology Press, 2009.

[17]SAMET H, WEBBER R E. Storing a Collection of Polygons Using Quadtrees[J]. ACM Transactions on Graphics (TOG), 1985, 4(3): 182-222.

[18]SAMET H, WEBBER R E. On Encoding Boundaries with Quadtrees[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1984, 6(3): 365-369.

[19]SHEWCHUK J R. Adaptive Precision Floating-point Arithmetic and Fast Robust Geometric Predicates[J]. Discrete & Computational Geometry, 1997, 18(3): 305-363.

[20]MULLER J M, BRISEBARRE N, DE Dinechin F, et al. Handbook of Floating-point Arithmetic[M]. Boston: Birkh?user, 2009.

[21]王舒鵬, 方莉. 混合積判斷線段相交的方法分析[J]. 電腦開發(fā)與應(yīng)用, 2006, 19(10): 34-35.

WANG Shupeng,FANG Li.An Analysis of Two Segments Intersection Judgment with Mixed Product[J]. Computer Development & Applications, 2006, 19(10): 34-35.

[22]EGENHOFER M J, HERRING J R. Categorizing Binary Topological Relations Between Regions, Lines and Points in Geographic Databases[R]. Orono: University of Maine 1991.

[23]虞強(qiáng)源, 劉大有, 謝琦. 空間區(qū)域拓?fù)潢P(guān)系分析方法綜述[J]. 軟件學(xué)報(bào), 2003, 14(4): 777-782.

YU Qiangyuan, LIU Dayou, XIE Qi. A Survey of Analysis Methods of Topological Relations between Spatial Regions [J]. Journal of Software, 2003, 14(4): 777-782.

[24]HERRING J R. OpenGIS?Implementation Standard for Geographic Information——Simple Feature Access-Part 1: Common architecture[S]. [S.l.]: OGC Document, 2011.

[25]MCKENNEY M, PAULY A, PRAING R, et al. Local Topological Relationships for Complex Regions[M]∥Advances in Spatial and Temporal Databases. Heidelberg: Springer, 2007: 203-220.

(責(zé)任編輯：張艷玲)

修回日期： 2015-05-04

First author： XIE Chuanjie(1971—), male, PhD, associate professor, majors in parallel geographic calculation.

E-mail： xiecj@leris.ac.cn

Parallel Algorithm of Spatial Relationship Query between Polygons Based on Heterogeneous Multi-core Architecture

XIE Chuanjie1,LONG Zhou1，2,MA Yihang1，2,YOU Zhijie1，2

1. State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing 100101,China; 2.University of Chinese Academy of Sciences, Beijing 100049

Abstract：The commonly used Plane Sweep algorithm in the spatial relationship query is a serial algorithm, when dealing with huge amounts of spatial data, the efficiency is very low，and the existing parallel computing method is not applicable to normal computer. Aiming at this problem, this paper proposes a parallel algorithm of spatial relationship query between polygons based on heterogeneous multi-core architecture. The algorithm uses STR-tree index to filter out non-intersection polygons first, then decomposes the filtered polygons data into points set and edges set, and constructs a quadtree index for them; under the premise of data accuracy meets the requirement of floating point calculation, uses GPU’s powerful batch computing power quickly processing the intersection between edges and calculates the topology relationship between rings, then uses the topology relationship between rings to calculate the DE-9IM model between polygons; compares the DE-9IM model with spatial relationship query condition and outputs the query results. At last, the efficiency and accuracy of the algorithm is verified by experiment.

Key words：heterogeneous multi-core; parallel computing; topological relations; querying spatial relationship

作者簡介：第一 謝傳節(jié)(1971—)，男，博士，副研究員，研究方向?yàn)椴⑿械乩碛?jì)算。

收稿日期：2014-09-04

基金項(xiàng)目：國家863計(jì)劃(2011AA120302;2011AA120306;2012AA12A401)；海洋公益性項(xiàng)目(201105033-6)

中圖分類號：P208

文獻(xiàn)標(biāo)識碼：A

文章編號：1001-1595(2016)01-0119-08

引文格式：謝傳節(jié)，龍舟，馬益杭,等.多邊形間空間關(guān)系查詢的異構(gòu)多核架構(gòu)并行算法[J].測繪學(xué)報(bào)，2016,45(1)：119-126.DOI:10.11947/j.AGCS.2016.20140463.

XIE Chuanjie,LONG Zhou,MA Yihang,et al.Parallel Algorithm of Spatial Relationship Query between Polygons Based on Heterogeneous Multi-core Architecture[J]. Acta Geodaetica et Cartographica Sinica,2016,45(1):119-126.DOI:10.11947/j.AGCS.2016.20140463.