張樹清，張策，楊典華，張俊巖，潘欣，姜春雷

（1.中國(guó)科學(xué)院東北地理與農(nóng)業(yè)生態(tài)研究所，吉林 長(zhǎng)春 130102；2.首都師范大學(xué)三維信息獲取與應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048）

0 引言

當(dāng)前，地理信息系統(tǒng)（GIS）正面臨著密集計(jì)算和海量數(shù)據(jù)兩大挑戰(zhàn)。傳統(tǒng)地理信息系統(tǒng)軟件中的串行算法由于在處理器、內(nèi)存以及數(shù)據(jù)傳輸方面的瓶頸，難以有效應(yīng)用，從而極大地限制了地理信息系統(tǒng)的應(yīng)用。并行計(jì)算是高性能計(jì)算的重要技術(shù)［1］，它通過(guò)將計(jì)算和數(shù)據(jù)重新組合形成新的任務(wù)集并分配到多個(gè)計(jì)算單元中，在計(jì)算單元并行執(zhí)行，極大地提高了計(jì)算效率，對(duì)科學(xué)的發(fā)展產(chǎn)生了深遠(yuǎn)的影響。這為推動(dòng)具有計(jì)算與數(shù)據(jù)密集特征的地理信息科學(xué)的創(chuàng)新與發(fā)展帶來(lái)了機(jī)遇［2，3］。并行計(jì)算的引入不僅極大地提高了GIS計(jì)算速度，同時(shí)可以重新思考GIS新的表達(dá)、管理和分析方法［4］。如何有效地利用并行計(jì)算技術(shù)已成為當(dāng)前地理信息系統(tǒng)迫切需要解決的科學(xué)問(wèn)題。

矢量疊加分析（Overlay）是對(duì)用戶感興趣區(qū)域不同數(shù)據(jù)層及其所含屬性進(jìn)行空間疊置分析［5］，Overlay屬于GIS集合操作，為GIS重要空間分析功能之一［1］。研究GIS矢量疊加數(shù)據(jù)高效并行算法，解決諸如海量復(fù)雜空間圖層或動(dòng)態(tài)數(shù)據(jù)快速計(jì)算問(wèn)題，提高GIS實(shí)時(shí)分析與動(dòng)態(tài)決策能力，將是十分必要的。然而，目前Overlay的并行，大多基于計(jì)算效率低、內(nèi)存資源消耗大的傳統(tǒng)掃描線法。因此有必要以簡(jiǎn)單要素模型的Overlay為基礎(chǔ)，從多邊形對(duì)象分類特征角度，研究高效的空間劃分策略與并行算法。盡管多邊形Overlay包含交集、并集和差集等復(fù)雜計(jì)算，但一般可進(jìn)一步劃分為兩類最基本操作：1）點(diǎn)包含，一個(gè)對(duì)象的一個(gè)（或者多個(gè)）節(jié)點(diǎn)包含于另外一個(gè)對(duì)象中；2）線穿越（線相交），一個(gè)對(duì)象的一條邊界與另外對(duì)象的邊界相交［6］。謝忠等［5］給出了詳細(xì)的線穿越（線相交）計(jì)算公式以及交點(diǎn)方向判別方法，有關(guān)“點(diǎn)包含”檢測(cè)方法也有較深入研究，常見的方法為射線法［7］。然而，這兩種類型只能處理拓?fù)漕愋蚾verlapping（部分重疊）及contain（完全包含）或contained by（完全被包含）情形，卻無(wú)法計(jì)算共線相鄰（adjacent）或共線包含（cover）或共線被包含（covered by）拓?fù)漕愋?。事?shí)上，多邊形Overlay的另一基本操作類型：“共線”（一個(gè)主多邊形和疊加多邊形有部分相同邊界的情形）需要深入研究。為此，本文以多邊形交集計(jì)算為例，深入研究了共線計(jì)算方法，并給出多重?cái)?shù)據(jù)包圍盒（Box）的對(duì)象空間劃分、篩選與并行計(jì)算策略。

1 交集邊界共線情形處理

如圖1，設(shè)交點(diǎn)為P，線段PiPi＋1相對(duì)方程為：

其中：t為點(diǎn)在線段P0Pi＋1中的相對(duì)位置，P0為多邊形第一個(gè)節(jié)點(diǎn)：

其中：t為雙精度數(shù)值，整數(shù)部分是節(jié)點(diǎn)編號(hào)，小數(shù)部分是交點(diǎn)在線段P0Pi＋1中的相對(duì)位置：t∈［0，1）。

圖1 交點(diǎn)P在節(jié)點(diǎn)Pi和Pi＋1中的相對(duì)位置Fig.1 Relative position of intersection Pbetween node Piand Pi＋1

共線數(shù)據(jù)鏈兩端點(diǎn)間交點(diǎn)數(shù)量計(jì)算：1）若2個(gè)交點(diǎn)都滿足t≠0，兩交點(diǎn)間存儲(chǔ)的節(jié)點(diǎn)數(shù)為兩交點(diǎn)的相對(duì)位置，先取整，再相減；2）若2個(gè)交點(diǎn)中1個(gè)位于節(jié)點(diǎn)上，兩交點(diǎn)間存儲(chǔ)的節(jié)點(diǎn)數(shù)為兩交點(diǎn)的相對(duì)位置，先相減，再取整；3）若2個(gè)交點(diǎn)都位于節(jié)點(diǎn)上，兩交點(diǎn)的相對(duì)位置相減，再減1。

共線性質(zhì)判斷定理：兩多邊形某一區(qū)間共線（不分內(nèi)、外環(huán)），交點(diǎn)在A中后一交點(diǎn)相對(duì)位置與前一交點(diǎn)相對(duì)位置差（DA）乘以交點(diǎn)在B中后一交點(diǎn)相對(duì)位置與前一交點(diǎn)相對(duì)位置差（DB），若乘積為正，共線是“交集面積區(qū)域”邊界線；若乘積為負(fù)，共線不是“交集面積區(qū)域”邊界線，只是兩多邊形相鄰邊界線。若共線為“交集面積區(qū)域”邊界線情形，需進(jìn)一步鏈接產(chǎn)生新的交集面積區(qū)域多邊形，而僅僅為兩多邊形相鄰邊界線的情形，無(wú)需鏈接產(chǎn)生新的交集面積區(qū)域多邊形。

圖2顯示了兩個(gè)對(duì)象A（灰色，外環(huán)節(jié)點(diǎn)為Pi，內(nèi)環(huán)節(jié)點(diǎn)為Si）和B（白色，節(jié)點(diǎn)為Qi）相交的兩類情形：1）圖2a，共線邊界不是交集面積區(qū)域邊界；2）圖2b、圖2c，共線邊界為“交集面積區(qū)域”（斜線填充區(qū)域）邊界線，其中圖2c中B對(duì)象含有島區(qū)（島區(qū)節(jié)點(diǎn)分別為S0、S1、S2、S3、S4）。為簡(jiǎn)便起見，圖中的交點(diǎn)假定與節(jié)點(diǎn)重合，因此可直接用節(jié)點(diǎn)參與計(jì)算相對(duì)位置。下面根據(jù)共線性質(zhì)定理判定：1）圖2a：DA＝2－3＝－1；DB＝3－2＝1；DA＊DB＝（－1）＊1＝－1＜0：共線為A、B相鄰邊界線，不是“交集面積區(qū)域”邊界線。2）圖2b：DA＝5－4＝1；DB＝3－2＝1；DA＊DB＝1＊1＝1＞0：是“交集面積區(qū)域”邊界線。3）圖2c：DA＝2－3＝－1；DB＝1－2＝－1；DA＊DB＝（－1）＊（－1）＝1＞0：“交集面積區(qū)域”邊界線。

圖2 交集區(qū)域共線示意Fig.2 Illustration of intersection region collinear

2 多邊形對(duì)象并行相交邏輯框架

2.1 多層數(shù)據(jù)包圍盒（Box）篩選結(jié)構(gòu)

在數(shù)據(jù)結(jié)構(gòu)中，無(wú)論是二維數(shù)據(jù)，還是三維數(shù)據(jù)，都可內(nèi)部設(shè)置一些Box，以減少查詢、檢索、空間分析時(shí)間，而不用進(jìn)行預(yù)處理。既然所有對(duì)象數(shù)據(jù)都是以對(duì)象為單位建模的，也就適用于以對(duì)象／子對(duì)象為單位進(jìn)行并行計(jì)算?，F(xiàn)有二維數(shù)據(jù)，除對(duì)象外，無(wú)Box?？赏ㄟ^(guò)不改變?cè)袛?shù)據(jù)，附加一文件，提取構(gòu)造Box，以降低查詢時(shí)間。

圖3給出了多層Box遴選示意圖，其中對(duì)象A其矩形包圍盒為BOX（A），由兩個(gè)環(huán)構(gòu)成A1和A2，用圓表示，其包圍盒分別為BOX（A1）和BOX（A2）；而環(huán)A1的邊界又可能分為n條鏈（或段），相應(yīng)鏈的包圍盒分別為 BOX（A11），BOX（A12），…，BOX（A1n）。同樣，對(duì)象B包圍盒為BOX（B），也由兩個(gè)環(huán)構(gòu)成B1、B2，包圍盒分別為BOX（B1）和 BOX（B2），設(shè)B1邊界可能分為m條鏈，其包圍盒分別為BOX（B11），BOX（B12），…，BOX（B1m）。這樣通過(guò)多層Box篩選，可以迅速檢索真正相交部分，如與B1相交的A1邊界位于BOX（A1i）和BOX（A1j）內(nèi)。這里“鏈”（或“段”）指：兩交點(diǎn)間（不含其它交點(diǎn)）的多邊形連續(xù)邊界，當(dāng)無(wú)交點(diǎn)存在時(shí)，多邊形的邊界構(gòu)成一條“閉合鏈”。

圖3 多層數(shù)據(jù)包圍盒空間對(duì)象相交示意Fig.3 Illustration of spatial object intersection with multi－level bounding boxes

求交以對(duì)象為單位做并行運(yùn)算。為降低求交時(shí)間，劃分方法如圖4：以吉黑兩省土地利用數(shù)據(jù)（JHLUCC）和土壤數(shù)據(jù)（JHSOIL）為例，根據(jù)報(bào)告類型簡(jiǎn)單對(duì)象多的為Ai（i∈［1，n］），不分段；對(duì)象少、分為多鏈的為Bij（i∈［1，u］，j∈［1，v］）（i為對(duì)象編號(hào)，v為該對(duì)象分段數(shù)，u為該類對(duì)象總數(shù)）。據(jù)此可建立復(fù)雜對(duì)象多線程或多進(jìn)程劃分與并行計(jì)算模式，故可將u個(gè)對(duì)象Bij分配在不同計(jì)算節(jié)點(diǎn)上，各節(jié)點(diǎn)單獨(dú)進(jìn)行Bij與所有A類對(duì)象（Ai（i∈［1，n］））（一對(duì)多）的疊加運(yùn)算（圖4）。

圖4 復(fù)雜對(duì)象Overlay多線程或多進(jìn)程任務(wù)劃分示意Fig.4 Illustration of multi－thread or multi－processor task partition for complex object Overlay

2.2 多邊形分類及其Overlay計(jì)算處理流程框架

在簡(jiǎn)單要素模型中，一般一個(gè)多邊形對(duì)象可包含如下類型：?jiǎn)苇h(huán) （Single Polygon，SP）、多環(huán)（Multi－polygon，MP）；而根據(jù)節(jié)點(diǎn)數(shù)量的多少，每個(gè)環(huán)又劃分為單鏈（Single Line，SL）或者多鏈（Multi－Lines，ML）。針對(duì)這些不同類型，可分別設(shè)置分支處理函數(shù)。通過(guò)讀取主圖層和疊加圖層數(shù)據(jù)，可分別統(tǒng)計(jì)主圖層（A）和疊加圖層（B）多邊形數(shù)據(jù)分類情形，并將主要類型設(shè)為主體處理函數(shù)。這里設(shè)定A（單環(huán)單鏈SPSL）對(duì)B（單環(huán)單鏈SPSL）為主函數(shù)，進(jìn)入第一層函數(shù)，若對(duì)象BOX相交，則繼續(xù)進(jìn)行“線穿越（線相交）”、“點(diǎn)包含”或“共線”計(jì)算，判定多邊形是否為邊界相交、共線或整體包含。第二層共有四分支函數(shù)中：1）當(dāng)對(duì)象A為單環(huán)單鏈（SPSL），對(duì)象B為單環(huán)多鏈（SPML），轉(zhuǎn)入該層分支函數(shù)1；2）當(dāng)對(duì)象A仍為單環(huán)單鏈（SPSL），對(duì)象B為多環(huán)（MP），轉(zhuǎn)入該層分支函數(shù)2，該函數(shù)與主函數(shù)的差異在于：對(duì)象B中多于一定節(jié)點(diǎn)的環(huán)要采用BOX篩選；3）當(dāng)對(duì)象A為多鏈（ML），對(duì)象B為一般對(duì)象時(shí)，轉(zhuǎn)入該層分支函數(shù)3，同樣需要對(duì)多鏈進(jìn)行BOX篩選，如此類推；4）當(dāng)對(duì)象A為多環(huán)（MP），對(duì)象B為一般對(duì)象時(shí)，轉(zhuǎn)入該層分支函數(shù)4。第三層分支函數(shù)用于處理第二層分支函數(shù)中函數(shù)3與4的情況，若對(duì)象A為單環(huán)多鏈，對(duì)象B為單環(huán)或多環(huán)時(shí)，分別轉(zhuǎn)入該層分支函數(shù)1和2；當(dāng)對(duì)象A為多環(huán)，對(duì)象B為單環(huán)時(shí)，轉(zhuǎn)入該層分支函數(shù)3；當(dāng)對(duì)象A為多環(huán)多鏈，對(duì)象B為多環(huán)時(shí)，轉(zhuǎn)入該層分支函數(shù)4（圖5），該類處理需要多重BOX篩選（圖3）。

圖5 不同類型多邊形Overlay分支處理流程Fig.5 Illustration of branch operation for different types of polygon Overlay

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)開發(fā)環(huán)境及數(shù)據(jù)

處理器：Intel Core 2DuoE7500＠1.6GHz；內(nèi)存：Ramaxel Technology PC2－6400DDR2 800MHz，2×2GB，800MHz；系統(tǒng)：Microsoft XP，Professional 2002，Service Pack 2；軟件：Microsoft Visual Studio 2005C?。目前Overlay程序是在單線程下實(shí)現(xiàn)，有關(guān)多線程或多進(jìn)程并行程序正在開發(fā)中。

采用的數(shù)據(jù)為：吉黑兩省土地利用數(shù)據(jù)（JHLUCC）和吉黑兩省土壤數(shù)據(jù)（JHSOIL），其數(shù)據(jù)格式：ESRI Shape。數(shù)據(jù)量：JHSOIL為47.7MB，JHLUCC為355.0MB；數(shù)據(jù)類型：JHLUCC中總計(jì)對(duì)象數(shù)量為179 111，其中單環(huán)對(duì)象數(shù)量為172 054，JHSOIL總計(jì)對(duì)象數(shù)量為16 728，單環(huán)對(duì)象數(shù)量為15 439，二者單環(huán)數(shù)據(jù)占絕大多數(shù)（大于90%），故可將單環(huán)對(duì)象求交設(shè)為主函數(shù)：（單環(huán)單鏈SPSL）對(duì)疊加圖層（單環(huán)單鏈SPSL）。

3.2 實(shí)驗(yàn)結(jié)果

以共線情形為例，本文給出了處理結(jié)果并驗(yàn)證其準(zhǔn)確性。圖6是本文軟件下兩對(duì)象共線相交的實(shí)例，JHLUCC為主層、JHSOIL為疊加層。其中，圖6a為JHLUCC中編號(hào)為7667多邊形與JHSOIL中編號(hào)為394多邊形的疊加結(jié)果，這里的共線邊界為“交集面積區(qū)域”邊界線；圖6b為JHLUCC中編號(hào)為93181與JHSOIL中編號(hào)為5429的疊加結(jié)果，其中的共線邊界是兩個(gè)多邊形的鄰域邊界線，而非“交集面積區(qū)域”邊界線，說(shuō)明所給出的算法能夠準(zhǔn)確處理實(shí)際情形。在本文軟件下，無(wú)論以JHLUCC（對(duì)象數(shù)量為179 111）為主圖層、JHSOIL（對(duì)象數(shù)量為16 728）為疊加圖層，還是以JHSOIL為主圖層、JHLUCC為疊加圖層，其多邊形疊加時(shí)間均大約為350s左右。

圖6 共線情形實(shí)際測(cè)試結(jié)果Fig.6 Real test result of collinear situation

4 討論與結(jié)論

多邊形Overlay是GIS重要分析功能，商業(yè)化軟件如ArcGIS等其Overlay版本不斷升級(jí)。然而在學(xué)術(shù)界對(duì)多邊形對(duì)象疊加的研究卻不夠系統(tǒng)：首先，已有研究只考慮了“線穿越（線相交）”和“點(diǎn)包含”基本類型［5，6］，對(duì)“共線”類型的研究涉足較少，無(wú)法處理共線相鄰（adjacent）或共線包含（cover）或共線被包含（covered by）拓?fù)漕愋?，為此，本文提出了共線節(jié)點(diǎn)數(shù)量的計(jì)算方法，以及共線是否為“交集區(qū)域邊界”的詳細(xì)判斷定理，從而可以有效處理各類空間拓?fù)潢P(guān)系。其次，多邊形對(duì)象類型劃分及其要素篩選策略的研究不夠深入，ESRI shape文件中僅包含對(duì)象Box，這仍不能有效地遴選復(fù)雜對(duì)象（如多環(huán)構(gòu)成的對(duì)象）的相交成分。事實(shí)上，在構(gòu)建復(fù)雜空間對(duì)象過(guò)程中應(yīng)顧及對(duì)象與子對(duì)象層次關(guān)系，而對(duì)于由大量節(jié)點(diǎn)構(gòu)成的子對(duì)象應(yīng)考慮子對(duì)象數(shù)據(jù)包圍盒，這里的子對(duì)象可以是多邊形構(gòu)成的閉合環(huán)（如島），也可以是由兩個(gè)交點(diǎn)間多邊形邊界構(gòu)成的“獨(dú)立幾何單元”，從而提出了面向?qū)ο髥苇h(huán)、多環(huán)、單鏈、多鏈特征的多重Box數(shù)據(jù)篩選及對(duì)象空間劃分與高效并行運(yùn)算策略。初步試驗(yàn)表明：該方法具有時(shí)間復(fù)雜度低、編程簡(jiǎn)單、內(nèi)存等資源消耗少等優(yōu)點(diǎn)。同時(shí)，通過(guò)有關(guān)新增數(shù)據(jù)包圍盒等，可保證與當(dāng)前數(shù)據(jù)的兼容性。以吉黑兩省土地利用（JHLUCC）和土壤（JHSOIL）ESRI Shape文件數(shù)據(jù)為例，對(duì)以上假定進(jìn)行初步實(shí)驗(yàn)，結(jié)果表明，所提出的方法可有效處理任何復(fù)雜的共線情形。但限于篇幅，本文僅給出了串行計(jì)算實(shí)驗(yàn)的初步結(jié)果，有關(guān)多線程、多進(jìn)程并行計(jì)算以及相關(guān)計(jì)算效率的比較等有待進(jìn)一步討論。

［1］ DING Y M，DENSHAM P J.Spatial strategies for parallel spatial modeling［J］.International Journal of Geographical Information Science，1996，10（6）：669－698.

［2］ WANG S，ARMSTRONG M P.A theoretical approach to the use of cyberinfrastructure in geographical analysis［J］.International Journal of Geographical Information Science，2009，23（2）：169－193.

［3］ WANG S，LIU Y.TeraGrid GIScience gateway：Bridging cyberinfrastructure and GIScience［J］.International Journal of Geographical Information Science，2009，23（5）：631－656.

［4］ DANGERMONJD J，MOREHOUSES S.Trends in hardware for geographic information systems［A］.Proceedings of Auto－Carto 8（Bethesda：American Congress for Surveying and Mapping）［C］.1987.380－385.

［5］ 謝忠，葉梓，吳亮.簡(jiǎn)單要素模型下多邊形疊置分析算法［J］.地理與地理信息科學(xué)，2007，23（3）：19－23.

［6］ ZHANG S Q，ZHANG J Y.Theoretical analytics of stereographic projection on 3Dobjects′intersection predicate［J］.International Journal of Geographical Information Science，2010，24（1）：25－46.

［7］ YEHUDA E K.Determining the spatial containment of a point in a general polyhedron［J］.Computer Graphics and Image Processing，1982，19（4）：303－334.