江錦成，吳立新，2，＊，孫文彬，楊宜舟

（1.北京師范大學(xué)民政部／教育部減災(zāi)與應(yīng)急管理研究院，北京 100875；2.中國(guó)礦業(yè)大學(xué)物聯(lián)網(wǎng)（感知礦山）國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，江蘇 徐州 221008；3.東北大學(xué)測(cè)繪遙感與數(shù)字礦山研究所，遼寧 沈陽(yáng) 110819；4.中國(guó)礦業(yè)大學(xué)（北京）地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083）

0 引言

空間鄰近分析是空間信息科學(xué)領(lǐng)域中一個(gè)極其重要的研究?jī)?nèi)容，在空間查詢與分析［1］、空間層次聚類［2］、DEM 內(nèi)插［3］、目標(biāo)預(yù)警［4］等領(lǐng)域得到廣泛應(yīng)用。目前，多采用Voronoi圖（簡(jiǎn)稱V圖）距離度量空間對(duì)象之間的鄰近程度［5，6］?；赩圖距離的K階鄰近［7］（V圖－K階鄰近）不僅可以描述相接目標(biāo)間的拓?fù)溧徑?，還能區(qū)分出相離目標(biāo)間的鄰近程度，具有廣泛的應(yīng)用價(jià)值。

盡管國(guó)內(nèi)外眾多學(xué)者已對(duì)V圖－K階鄰近串行算法的優(yōu)化進(jìn)行了大量的研究，但隨著實(shí)際應(yīng)用數(shù)據(jù)規(guī)模的不斷增加，導(dǎo)致V圖－K階鄰近算法計(jì)算時(shí)間急劇增長(zhǎng)，成為約束算法應(yīng)用的瓶頸。近年來(lái)，計(jì)算機(jī)并行技術(shù)得到了快速發(fā)展，已成為解決串行算法計(jì)算效率低下的有效手段之一［8］。國(guó)內(nèi)外已有學(xué)者對(duì)V圖－K 階鄰近算法進(jìn)行了并行化［9，10］，但大多是基于共享存儲(chǔ)的并行模型，易受存儲(chǔ)容量的限制，且可移植性差。而基于消息傳遞的并行模型，支持進(jìn)程間的分布式存儲(chǔ)，且具有更好的可移植性和擴(kuò)展性［11］。

本文提出一種基于消息傳遞模型的V圖－K階鄰近并行搜索算法——PVKN（Parallel Voronoi K－order Neighbors）算法，擬采用分治法構(gòu)建點(diǎn)集V圖，應(yīng)用波浪法搜索K階鄰近目標(biāo)，利用MPI（Mes－sage Passing Interface）并行技術(shù)分別對(duì)它們進(jìn)行并行化，以提高V圖－K階鄰近的搜索效率。

1 V圖－K階鄰近的描述

V圖距離可用來(lái)描述不同空間實(shí)體間的鄰近程度，其定義為：設(shè)P是二維空間目標(biāo)P1，P2，…，Pn的集合，Pi，Pj∈P（i，j＝1，2，…，n）所屬的 V圖單元記為PVi和PVj，則將Pi、Pj之間最少的V圖單元數(shù)K，稱為Pi與Pj的V圖距離，記為νd（Pi，Pj）。若PVi、PVj存在νd（Pi，Pj）＝K，則稱Pi與Pj之間存在V圖－K階鄰近關(guān)系。

如圖1所示，圓點(diǎn)為源目標(biāo)點(diǎn)，即0階鄰近目標(biāo)；正方形和三角形點(diǎn)集分別位于源目標(biāo)的V圖－1階和2階鄰近單元內(nèi)。

圖1 V圖－2階鄰近示意Fig.1 An example of Voronoi 2－order neighbors

2 V圖－K鄰近算法并行化分析

點(diǎn)集V圖－K階鄰近的求解過程主要包括：V圖構(gòu)建和K階鄰近搜索。目前，可用于V圖構(gòu)建和K階鄰近搜索的算法較多，但各算法間的時(shí)間復(fù)雜度（效率）或功能存在著很大差異。

2.1 構(gòu)建V圖

構(gòu)建點(diǎn)集V圖的算法主要有：逐點(diǎn)插入法、間接法及分治法［12］，時(shí)間復(fù)雜度分別為 O（n2）、O（n2）和O（n×log（n））。從時(shí)間復(fù)雜度來(lái)看，分治法的效率最高，且采用“分而治之”的思想［13］，即將數(shù)據(jù)集按區(qū)域鄰近的原則劃分成若干個(gè)子區(qū)域，先分別生成各區(qū)域的V圖，再合并成整個(gè)V圖。其中，各子區(qū)域V圖的構(gòu)建是相互獨(dú)立的，易于并行計(jì)算。

2.2 搜索K階鄰近

搜索K階鄰近的主要方法包括：對(duì)向法、穿越法和波浪法［6］。對(duì)向法僅為獲取兩給定目標(biāo)間的鄰近階數(shù)；穿越法主要用于計(jì)算兩目標(biāo)在特定方向上的鄰近目標(biāo)集，也可計(jì)算兩目標(biāo)間的鄰近階數(shù)，但所求階數(shù)是近似解；波浪法不僅可求解源目標(biāo)的K階鄰近目標(biāo)集，還可用于計(jì)算兩目標(biāo)間的鄰近階數(shù)，且所求階數(shù)為精確解。故本文選用波浪法進(jìn)行K階鄰近搜索。

波浪法搜索單源目標(biāo)的第k＋1階鄰近時(shí)，先搜索出第k階鄰近，再搜索第k階鄰近的直接鄰近目標(biāo)，并從這些目標(biāo)中剔除第k－1階目標(biāo)，即可得到第k＋1階鄰近目標(biāo)。由于各階鄰近之間具有很強(qiáng)的關(guān)聯(lián)性，無(wú)法應(yīng)用簡(jiǎn)單的靜態(tài)任務(wù)分配來(lái)實(shí)現(xiàn)并行化，需采用動(dòng)態(tài)調(diào)度策略。而當(dāng)對(duì)多個(gè)源目標(biāo)進(jìn)行搜索時(shí)，各源目標(biāo)間的計(jì)算相互獨(dú)立，將多個(gè)源目標(biāo)分配至不同進(jìn)程后，各進(jìn)程獨(dú)立搜索K階鄰近即可實(shí)現(xiàn)并行計(jì)算。

3 V圖－K階鄰近并行算法設(shè)計(jì)

利用MPI并行技術(shù)分別對(duì)V圖構(gòu)建和K階鄰近搜索過程進(jìn)行并行化，以提升時(shí)間效率。PVKN算法流程如圖2所示（設(shè)共p個(gè)進(jìn)程參與計(jì)算），主要包括以下步驟：1）讀取點(diǎn)數(shù)據(jù)集；2）按區(qū)域鄰近的原則對(duì)點(diǎn)數(shù)據(jù)集進(jìn)行劃分，并將分割后的局部數(shù)據(jù)分配至各進(jìn)程；3）各進(jìn)程先獨(dú)立構(gòu)建局部V圖，再通過迭代合并，生成全局V圖；4）各進(jìn)程在V圖的基礎(chǔ)上，并行搜索源目標(biāo)（集）的K階鄰近目標(biāo)。步驟3和步驟4是PVKN算法的主要耗時(shí)部分，對(duì)其加速效果直接決定整個(gè)算法的效率。

圖2 PVKN算法流程Fig.2 Flow chart of PVKN algorithm

3.1 并行構(gòu)建V圖

為實(shí)現(xiàn)原始數(shù)據(jù)集的快速分割，需先對(duì)點(diǎn)數(shù)據(jù)集按坐標(biāo)值排序?；诳焖倥判蚍ǖ牟⑿姓齽t采樣排序——PSRS（Parallel Sorting by Regular Sampling）算法［14－16］不僅效率高，還能保證各進(jìn)程間的負(fù)載基本均衡，故本文采用PSRS算法對(duì)點(diǎn)數(shù)據(jù)集進(jìn)行并行排序。

排序完成后，可應(yīng)用并行分治法構(gòu)建V圖，其基本思想是先構(gòu)建局部V圖，再通過兩兩迭代合并得到全局V圖，主要步驟如下：1）各進(jìn)程使用分治法分別構(gòu)建局部V圖；2）初始化合并次數(shù)t＝1；3）進(jìn)程2t×k（k＜）接收進(jìn)程2t×k＋2t－1發(fā)送的子圖消息，并在進(jìn)程2t×k中用分治法合并兩子圖；4）更新t＝t＋1，若2t－1＜p，則轉(zhuǎn)至步驟3，否則結(jié)束計(jì)算。

為使所有（或多個(gè)）進(jìn)程都參與后續(xù)的K階鄰近搜索，應(yīng)將合并到0號(hào)進(jìn)程的全局V圖廣播至所有（或多個(gè)）進(jìn)程。

3.2 并行搜索K階鄰近

波浪法采用逐層擴(kuò)張的思想，從源目標(biāo)開始向周圍逐步搜索鄰近目標(biāo)，先尋找第k（k＜K）階鄰近目標(biāo)的所有直接鄰近，再剔除掉第k－1階鄰近目標(biāo)，即可得到第k＋1階鄰近。如此可求出源目標(biāo)的1至K階或第K階鄰近目標(biāo)。

在實(shí)際應(yīng)用中，存在著單源和多源目標(biāo)鄰近搜索兩種情況。單源目標(biāo)時(shí)，只有搜索完所有第k階鄰近后，方可搜索第k＋1階鄰近，故只能在搜索第k階鄰近目標(biāo)的直接鄰近時(shí)進(jìn)行并行化；而對(duì)多源目標(biāo)進(jìn)行搜索時(shí)，不同源目標(biāo)間的計(jì)算相互獨(dú)立，易于并行計(jì)算。針對(duì)以上兩種情況，分別采用以下策略進(jìn)行并行化：

（1）單源目標(biāo)K階鄰近：1）主進(jìn)程將源目標(biāo)v廣播至所有進(jìn)程；2）進(jìn)程i（0≤i≤p－1）初始化當(dāng)前搜索階數(shù)k＝0、第k階鄰近Sk（i）?ν、第k－1階鄰近Sk－1（i）＝?；3）進(jìn)程i搜索Sk（i）的直接鄰近S′k＋1（i），得到第k＋1階鄰近Sk＋1（i）＝S′k＋1（i）－Sk－1（i）；4）主（），進(jìn)程m先匯總所有進(jìn)程的Sk＋1i即Sk＋1（i），并剔除Sk＋1（m）中重復(fù)的鄰近目標(biāo)，再剔除第k階和k－1階鄰近，即：Sk＋1（m）＝Sk＋1（m）－Sk（m）－Sk－1（m）；5）主進(jìn)程m 將Sk＋1（m）均勻地分成p等份，并分發(fā)至各進(jìn)程，進(jìn)程i更新Sk－1（i）＝Sk（i），Sk（i）為接收到的鄰近目標(biāo)；6）更新k＝k＋1，若k＞K，則結(jié)束搜索，否則轉(zhuǎn)至步驟3。最后主進(jìn)程m計(jì)算得到的Sk（m）即為源目標(biāo)的K階鄰近。

（2）多源目標(biāo)K階鄰近：1）主進(jìn)程將多個(gè)源目標(biāo)平均分成p等份S0（i）（0≤i≤p－1），并分發(fā)至各進(jìn)程；2）進(jìn)程i分別對(duì)S0（i）中的每個(gè)源目標(biāo)都使用串行波浪法求解K階鄰近；3）各進(jìn)程分別存儲(chǔ)每個(gè)源目標(biāo)的K階鄰近。

4 實(shí)驗(yàn)

為測(cè)試算法的正確性和時(shí)間效率，在圖形工作站（CPU：32核，2.0GHz；內(nèi)存：48GB）上，采用MPICH2并行庫(kù)進(jìn)行消息傳遞，應(yīng)用1∶500萬(wàn)數(shù)字地理底圖居民地（點(diǎn)數(shù)為5 042）數(shù)據(jù)對(duì)本文PVKN算法進(jìn)行相關(guān)實(shí)驗(yàn)。在不考慮I／O的情況下，測(cè)試結(jié)果如圖3－圖8所示。

圖3 點(diǎn)集V圖－3鄰近Fig.3 The 3－order neighbors of PVKN algorithm

圖3是求解單源目標(biāo)3階鄰近的結(jié)果?？梢钥闯觯琍VKN算法可正確構(gòu)建出點(diǎn)集V圖，并能搜索出源目標(biāo)點(diǎn)（三角形）的V圖－3階鄰近點(diǎn)（圓點(diǎn)），驗(yàn)證了算法的正確性。

PVKN算法主要由并行構(gòu)建V圖和并行搜索K階鄰近組成，本文測(cè)試并對(duì)比了其運(yùn)行時(shí)間。

（1）并行構(gòu)建V圖：圖4、圖5分別為并行構(gòu)建V圖的運(yùn)行時(shí)間和加速比。由圖可知，隨著進(jìn)程數(shù)的增加，構(gòu)建V圖的時(shí)間大幅度減少；加速比呈近似線性增長(zhǎng)，8進(jìn)程時(shí)，可達(dá)5倍以上。

圖4 并行構(gòu)建V圖的時(shí)間Fig.4 Time efficiency test of parallel constructing Voronoi diagrams

圖5 并行構(gòu)建V圖的加速比Fig.5 Speedup of parallel constructing Voronoi diagrams

（2）并行搜索K階鄰近：分別在不同K階數(shù)時(shí)，對(duì)單源目標(biāo)和多源目標(biāo)的K階鄰近搜索時(shí)間進(jìn)行測(cè)試，結(jié)果如圖6、圖7所示，其中圖6、圖7a的圖例表示K階數(shù)，圖7b的圖例表示源目標(biāo)數(shù)。

圖6為搜索單源目標(biāo)K階鄰近的耗時(shí)情況。由圖可知：在單源目標(biāo)情況下，隨著K階數(shù)的增加，搜索時(shí)間呈增長(zhǎng)趨勢(shì)；但由于搜索過程中存在著頻繁的消息傳遞與數(shù)據(jù)剔重過程，代價(jià)較高，導(dǎo)致加速效果不佳。因此，在單源目標(biāo)情況下，搜索K階鄰近過程不適宜并行計(jì)算。

圖7a顯示的是源目標(biāo)數(shù)固定為600，K階數(shù)變化時(shí)搜索時(shí)間與進(jìn)程數(shù)的關(guān)系。由圖可知，搜索時(shí)間隨著K階數(shù)的增加而增長(zhǎng)；對(duì)同一K階數(shù)，隨著進(jìn)程數(shù)的增多，搜索時(shí)間不斷減少，加速效果明顯，8進(jìn)程時(shí)加速比可達(dá)到5倍以上。圖7b顯示的是K階數(shù)固定為12，源目標(biāo)數(shù)變化時(shí)搜索時(shí)間與進(jìn)程數(shù)之間的關(guān)系。由圖可知，隨著源目標(biāo)數(shù)的增多，計(jì)算耗時(shí)呈增長(zhǎng)趨勢(shì)；對(duì)相同數(shù)量的源目標(biāo)，隨著進(jìn)程數(shù)的增加，計(jì)算時(shí)間快速下降，8進(jìn)程時(shí)，加速比亦可達(dá)到5倍以上。

圖6 單源目標(biāo)K階鄰近的并行搜索時(shí)間Fig.6 Time efficiency test of parallel searching neighbors for single source point

圖7 多源目標(biāo)K階鄰近的并行搜索時(shí)間Fig.7 Time efficiency test of parallel searching neighbors for multi－source points

為分析V圖構(gòu)建與K階鄰近搜索對(duì)PVKN算法時(shí)間效率的影響，本文對(duì)比了兩個(gè)過程的運(yùn)行時(shí)間（如圖8，其中，橫坐標(biāo)為源目標(biāo)數(shù)，縱坐標(biāo)表示V圖構(gòu)建與K階鄰近搜索的時(shí)間之比，圖例表示K階數(shù)）。

圖8 V圖構(gòu)建與K階鄰近搜索時(shí)間對(duì)比Fig.8 The time ratio of constructing Voronoi diagrams and searching K－order neighbors

由圖8可知，求解單源目標(biāo)（源目標(biāo)數(shù)為1）的K階鄰近時(shí)，構(gòu)建V圖時(shí)間遠(yuǎn)大于搜索K階鄰近的時(shí)間（約1 000倍）；而隨著源目標(biāo)數(shù)的增加，構(gòu)建V圖與搜索K階鄰近的時(shí)間之比迅速降低，且K越大，下降速度越快。可見：當(dāng)源目標(biāo)數(shù)或K階數(shù)較小時(shí)，構(gòu)建V圖耗時(shí)較長(zhǎng)，其加速效果決定PVKN算法的效率；而隨著K階數(shù)和源目標(biāo)數(shù)的增加，K階鄰近搜索時(shí)間快速增長(zhǎng)，成為影響PVKN算法的主要因素。

5 結(jié)論

本文使用并行分治法構(gòu)建點(diǎn)集V圖，應(yīng)用波浪法在不同的并行策略下，分別對(duì)單源目標(biāo)和多源目標(biāo)的K階鄰近搜索過程進(jìn)行了并行化。實(shí)驗(yàn)結(jié)果表明：?jiǎn)卧茨繕?biāo)時(shí)，構(gòu)建V圖占據(jù)主要時(shí)間，其良好的加速效果有效地提升了PVKN算法的運(yùn)行效率；多源目標(biāo)時(shí)，并行構(gòu)建V圖和搜索K階鄰近都有明顯的加速效果。兩種情況下，PVKN算法的加速比在8進(jìn)程時(shí)，都能達(dá)到5倍以上。與傳統(tǒng)串行算法相比，PVKN算法的運(yùn)行效率更高，為快速構(gòu)建大量點(diǎn)集V圖及基于V圖搜索K階鄰近提供了新方法。

［1］ CHRISTOPHER B J，MARK W J.Nearest neighbor search for linear and polygonal objects with contrained triangulations［A］.Proc.8th International Symposium on Spatial Data Handing，1998.13－21.

［2］ 宋曉眉，程昌秀，周成虎，利用K階空間鄰近圖的空間層次聚類方法［J］.武漢大學(xué)學(xué)報(bào)，2010，35（12）：1496－1499.

［3］ 張俊前，陳浩，趙伶俐.基于Voronoi K階鄰近的DEM內(nèi)插方法研究［J］.湖南科技大學(xué)學(xué)報(bào)，2010，25（2）：85－88.

［4］ 司海棠，秦小麟，郝學(xué)峰.基于Voronoi K階鄰近的目標(biāo)預(yù)警預(yù)報(bào)［J］.計(jì)算機(jī)應(yīng)用，2009，29（2）：598－601.

［5］ PREPARATA F P，SHAMOS M I.Computational Geometry：An Introduction［M］.New York，1985.

［6］ 趙仁亮，基于Voronoi圖的空間關(guān)系計(jì)算研究［D］.中南大學(xué)，2002.

［7］ CHEN J，ZHAO R L，LI Z L.Voronoi－based K－order neighbor relations for spatial analysis［J］.Advanced Techniques for Analysis of Geo－spatial Data，2004，59（2）：60－72.

［8］ BARNEY B.Introduction to parallel computing［EB／OL］.Available at：https：／／computing.llnl.gov／tut－orials／parallel comp／，2010－03－31.

［9］ 余婧，曹菡，靳朋飛.Voronoi圖K階鄰近并行矩陣迭代算法［J／OL］.計(jì)算機(jī)工程與應(yīng)用，2012.

［10］ LIANG S S，LIU Y，WANG C，et al.Design and evaluation of aparallel K－nearest neighbor algorithm on CUDA－enabled GPU［A］.Proceeding 2010IEEE 2ed Symposium on Web Society［C］.2010.

［11］ 張林波，遲學(xué)斌.并行計(jì)算導(dǎo)論［M］.北京：清華大學(xué)出版社，2006.

［12］ 吳立新，史文中.地理信息系統(tǒng)原理與算法［M］.北京：科學(xué)出版社，2003.

［13］ SHAMOS M I，HOEY D.Closest－point problems［A］.Proceeding of the 16th Annual IEEE Symposium on Foundation of Computer Science［C］.1975.151－162.

［14］ 王永剛.排序算法綜述［J］.電腦知識(shí)與技術(shù)，2006，29（7）：1－6.

［15］ HOARE C A R.Quicksort［J］.Computer，1962，5（1）：10－16.

［16］ LI X B，LU P，SCHAEFFER J，et al.On the versatility of parallel sorting by regular sampling［J］.Parallel Computing，1993，19（10）：1079－1103.