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

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

1. 中國(guó)科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100101； 2. 中國(guó)科學(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)系查詢(xún)的異構(gòu)多核架構(gòu)并行算法

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

1. 中國(guó)科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100101； 2. 中國(guó)科學(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)系查詢(xún)中常用的Plane Sweep算法是一種串行算法，在處理海量空間數(shù)據(jù)時(shí)效率較低，而已有的并行計(jì)算方法對(duì)于普通的計(jì)算機(jī)并不適用。本文針對(duì)這個(gè)問(wèn)題，提出了一種多邊形間空間關(guān)系查詢(xún)的異構(gòu)多核架構(gòu)并行算法，該算法先利用STR樹(shù)索引過(guò)濾掉不相交的多邊形，然后將過(guò)濾后的多邊形數(shù)據(jù)集合分解為點(diǎn)集合和邊集合，并對(duì)其構(gòu)建四叉樹(shù)索引；在保證數(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)系查詢(xún)條件相比對(duì)，輸出查詢(xún)結(jié)果。最后通過(guò)試驗(yàn)驗(yàn)證了算法的準(zhǔn)確性與高效性。

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

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

由于計(jì)算機(jī)硬件系統(tǒng)的不斷升級(jí)換代，基于異構(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]。

針對(duì)傳統(tǒng)算法在處理海量空間數(shù)據(jù)時(shí)存在的不足，本文打破原有的計(jì)算方法，設(shè)計(jì)了一并行計(jì)算方法，該方法先利用CPU+GPU架構(gòu)并行計(jì)算得到兩個(gè)圖層中所有最小外包矩形相交的多邊形的線(xiàn)段相交情況，然后根據(jù)線(xiàn)段相交情況判斷多邊形環(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)系查詢(xún)的并行算法

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

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

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

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

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

本文在進(jìn)行空間關(guān)系查詢(xún)時(shí)，首先遍歷兩個(gè)多邊形數(shù)據(jù)集合的所有多邊形要素，將所有的點(diǎn)和線(xiàn)段信息轉(zhuǎn)換并統(tǒng)一整合到同一個(gè)具有拓?fù)涮匦缘墓?jié)點(diǎn)和邊集合中。同時(shí)，獲取的空間信息還包括原始數(shù)據(jù)的最大空間范圍以及節(jié)點(diǎn)集中節(jié)點(diǎn)的個(gè)數(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ì)算量，首先需要先對(duì)兩個(gè)圖層的多邊形進(jìn)行過(guò)濾，根據(jù)多邊形的最小外包矩形過(guò)濾掉一定不相交的多邊形組合。本文選擇STR樹(shù)索引進(jìn)行過(guò)濾，STR樹(shù)索引是一種使用Sort-Tile-Recursive(STR)算法創(chuàng)建的只負(fù)責(zé)查詢(xún)處理的R樹(shù)。它針對(duì)二維空間數(shù)據(jù)，十分易于操作并且可以最大化地利用內(nèi)存空間存儲(chǔ)數(shù)據(jù)[15]。主要做法是：假設(shè)有兩個(gè)多邊形集合分別為A和B，先對(duì)A和B分別構(gòu)建STR樹(shù)，將多邊形的最小外包矩形作為STR樹(shù)節(jié)點(diǎn)的空間范圍，多邊形本身作為STR樹(shù)節(jié)點(diǎn)中存儲(chǔ)的空間對(duì)象，然后遍歷B中多邊形的最小外包矩形在A的STR樹(shù)中查詢(xún)，遍歷A中多邊形的最小外包矩形在B的STR樹(shù)中查詢(xún)，最終得到兩個(gè)過(guò)濾后的多邊形集合A和B。

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

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

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

圖3　空間四叉樹(shù)索引建立過(guò)程Fig.3　The process of establishing the spatial quadtree index

1.2線(xiàn)段相交判斷

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

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

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

圖4　4種線(xiàn)段相交的類(lèi)型Fig.4　Four types of line segments intersect

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

一個(gè)多邊形是由多個(gè)環(huán)組成，一個(gè)環(huán)是由多條邊組成。欲求兩個(gè)多邊形之間的空間關(guān)系，可以先從環(huán)間的拓?fù)潢P(guān)系計(jì)算入手，而環(huán)間的拓?fù)潢P(guān)系需要根據(jù)線(xiàn)段相交類(lèi)型進(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ì)算的兩個(gè)環(huán)為Ha、Hb，則根據(jù)線(xiàn)段相交類(lèi)型判斷兩環(huán)關(guān)系的主要過(guò)程如下：①若兩環(huán)沒(méi)有相交線(xiàn)段出現(xiàn)，則關(guān)系一定為Disjoint、Inside或Contains中的一種，再通過(guò)兩環(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；④若除了普通相交，其他相交類(lèi)型均有出現(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樹(shù)索引過(guò)濾后得到兩個(gè)多邊形集合，需要依次對(duì)有可能相交的多邊形組合判斷它們之間的空間關(guān)系，為了進(jìn)一步提高效率，本文使用BOOST庫(kù)多線(xiàn)程并行判斷多邊形間的空間關(guān)系。并行判斷涉及數(shù)據(jù)劃分，由于本文計(jì)算量最大的線(xiàn)段相交的判斷已經(jīng)在最開(kāi)始由CPU和GPU協(xié)作完成，則對(duì)于每一對(duì)多邊形組合，只需調(diào)用計(jì)算好的線(xiàn)段相交類(lèi)型來(lái)判斷環(huán)間拓?fù)潢P(guān)系，進(jìn)而判斷它們的空間關(guān)系，所以在計(jì)算量上不會(huì)有很大差距。那么對(duì)多邊形數(shù)據(jù)的劃分，可以按照有可能相交的目標(biāo)多邊形按數(shù)量平均分組，以保證每個(gè)線(xiàn)程的負(fù)載基本相同，然后將每一組多邊形組合分配給一個(gè)線(xiàn)程進(jìn)行空間關(guān)系的判斷，最后合并計(jì)算的結(jié)果。

多邊形之間的空間關(guān)系由維度擴(kuò)展九交模型(DE-9IM)來(lái)表達(dá)，根據(jù)維度擴(kuò)展九交模型的參數(shù)值，比較空間關(guān)系的查詢(xún)條件，即可輸出滿(mǎn)足條件的多邊形集合。根據(jù)OGC的標(biāo)準(zhǔn)，與DE-9IM參數(shù)值對(duì)應(yīng)的7類(lèi)空間關(guān)系查詢(xú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)系分析模塊，分別得到的兩個(gè)多邊形內(nèi)、外環(huán)之間的拓?fù)潢P(guān)系，利用此拓?fù)潢P(guān)系逐步判斷兩個(gè)多邊形間內(nèi)部、邊界和外部的相交值，從而確定多邊形間的DE-9IM參數(shù)值，最后將計(jì)算出的DE-9IM參數(shù)值與7類(lèi)空間關(guān)系對(duì)應(yīng)的DE-9IM參數(shù)進(jìn)行比對(duì)，以最終確定多邊形間的空間關(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。

編程語(yǔ)言：C++、CUDA。

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

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

任意兩個(gè)帶洞多邊形(即除了一個(gè)外環(huán)還可能有任意個(gè)內(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)系都可能對(duì)應(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

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

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

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

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

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

Tab.1The test results of parallel algorithm and ArcGIS

序號(hào)拓?fù)潢P(guān)系對(duì)應(yīng)查詢(xún)條件并行算法查詢(xún)結(jié)果ArcGIS查詢(xún)結(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組件式開(kāi)發(fā)技術(shù)，制作了一套僅具有數(shù)據(jù)輸入輸出、“按位置選擇”和時(shí)間記錄3種功能的小型測(cè)試軟件。先記錄各個(gè)空間關(guān)系查詢(xún)條件下并行算法的測(cè)試時(shí)間，然后利用開(kāi)源代碼庫(kù)GEOS中的空間關(guān)系計(jì)算測(cè)試時(shí)間，再利用小軟件記錄ArcGIS在各個(gè)空間關(guān)系查詢(xún)條件下的測(cè)試時(shí)間，最后對(duì)比3個(gè)測(cè)試的計(jì)算結(jié)果和測(cè)試時(shí)間。

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

將并行算法計(jì)算出的多邊形個(gè)數(shù)，與ArcGIS中“按位置選擇”功能相對(duì)應(yīng)的空間關(guān)系查詢(xún)條件計(jì)算得到的多邊形個(gè)數(shù)相對(duì)比，試驗(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)系對(duì)應(yīng)于“按位置選擇”功能中的條件描述ArcGIS計(jì)算結(jié)果并行算法測(cè)試結(jié)果相等與源圖層要素完全相同60006000相交與源圖層要素相交62486248接觸接觸源圖層要素的邊界77包含包含源圖層要素1616被包含在源圖層要素范圍內(nèi)1111

本文對(duì)5種空間關(guān)系查詢(xún)的計(jì)算時(shí)間進(jìn)行了統(tǒng)計(jì)，包括數(shù)據(jù)預(yù)處理(構(gòu)建STR樹(shù)索引及幾何要素轉(zhuǎn)化為圖形要素)，構(gòu)建四叉樹(shù)索引，線(xiàn)段相交判斷(包括精度判斷時(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ì)對(duì)比，每個(gè)算法的測(cè)試時(shí)間對(duì)比如圖8所示。經(jīng)過(guò)對(duì)比可以發(fā)現(xiàn)，在5種關(guān)系查詢(xún)的時(shí)間測(cè)試中，并行算法的測(cè)試時(shí)間比GEOS和ArcGIS的計(jì)算時(shí)間都要短，在計(jì)算效率上具有較大優(yōu)勢(shì)。以相等關(guān)系查詢(xú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并行算法對(duì)5種空間關(guān)系查詢(xún)時(shí)間統(tǒng)計(jì)

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

s

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

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

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

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

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

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

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

3結(jié)論

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

參考文獻(xiàn)：

[1]王結(jié)臣, 王豹, 胡瑋, 等. 并行空間分析算法研究進(jìn)展及評(píng)述[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)系的基本問(wèn)題與研究進(jìn)展[J]. 測(cè)繪學(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]陳軍, 劉萬(wàn)增, 李志林, 等. 線(xiàn)目標(biāo)間拓?fù)潢P(guān)系的細(xì)化計(jì)算方法[J]. 測(cè)繪學(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]陳斐, 周曉光. 基于平面分割的線(xiàn)/線(xiàn)相接、交叉類(lèi)型判斷[J]. 測(cè)繪學(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]. 北京: 中國(guó)水利水電出版社, 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]王舒鵬, 方莉. 混合積判斷線(xiàn)段相交的方法分析[J]. 電腦開(kāi)發(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ǎn)介：第一 謝傳節(jié)(1971—)，男，博士，副研究員，研究方向?yàn)椴⑿械乩碛?jì)算。

收稿日期：2014-09-04

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

中圖分類(lèi)號(hào)：P208

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

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

引文格式：謝傳節(jié)，龍舟，馬益杭,等.多邊形間空間關(guān)系查詢(xún)的異構(gòu)多核架構(gòu)并行算法[J].測(cè)繪學(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.