田 松，崔 希 民，孫 云 華，宮 宇

(1.中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院，北京 100083；2.國家海洋局南海分局，廣東 廣州 510310)

?

等級Voronoi圖及加權(quán)Voronoi圖的ArcGIS矢量生成算法

田 松1，崔 希 民1，孫 云 華1，宮 宇2

(1.中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院，北京 100083；2.國家海洋局南海分局，廣東 廣州 510310)

等級Voronoi圖及加權(quán)Voronoi圖是以樹狀形式表達(dá)事物層次關(guān)系的方法，可作為一種空間數(shù)據(jù)模型應(yīng)用于城鎮(zhèn)等級體系、組織結(jié)構(gòu)關(guān)系表達(dá)等地學(xué)領(lǐng)域。該文首次提出一種等級Voronoi圖及加權(quán)Voronoi圖的ArcGIS矢量生成算法，以ArcGIS Engine為開發(fā)工具，利用增量法思想，通過區(qū)域分割和合并方法實現(xiàn)了對空間的多級劃分，為研究和發(fā)展GIS空間數(shù)據(jù)模型提供了重要的方法手段。

等級Voronoi圖；等級加權(quán)Voronoi圖；ArcGIS；區(qū)域分割；區(qū)域合并

0 引言

等級Voronoi圖(Hierarchical Voronoi Diagram)亦稱作分層Voronoi圖，簡稱等級V圖或分層V圖,是對空間區(qū)域多級劃分，以樹狀形式表示區(qū)域?qū)哟侮P(guān)系的方法，是Voronoi圖理論和應(yīng)用的擴展。等級V圖是一種數(shù)據(jù)結(jié)構(gòu)，用于提高數(shù)據(jù)組織、檢索及可視化效率，亦是一種空間數(shù)據(jù)模型，用于表達(dá)空間事物間的等級層次關(guān)系。Farin等[1]利用等級V圖進行地形建模;Gold等[2]將等級V圖用于地形數(shù)據(jù)的組織與LOD顯示;Wallgrün[3]利用等級V圖解決了機器人路徑搜索問題；趙威、陳笑筑、薛瑩等[4-6]分別利用等級V圖表達(dá)了開封市、畢節(jié)市、蘇浙滬城市群的城市等級體系。

等級加權(quán)Voronoi圖(等級加權(quán)V圖)是等級Voronoi圖的擴展，是在考慮母點權(quán)重的情況下，劃分多級空間的方法。Balzer等[7]研究了等級V圖及加法加權(quán)等級V圖的生成算法，并用于軟件度量可視化；Sud等[8]利用CPU加速技術(shù)實現(xiàn)了加法加權(quán)的等級V圖；Kang等[9]利用等級加權(quán)V圖進行空間數(shù)據(jù)的組織和管理。目前，尚無結(jié)合GIS軟件或工具實現(xiàn)等級V圖及加權(quán)圖的方法。本文在文獻(xiàn)[10]方法的基礎(chǔ)上，利用ArcGIS Engine(簡稱AE)實現(xiàn)了等級V圖及加權(quán)V圖的矢量生成方法，旨在建立一種適用于地學(xué)領(lǐng)域的空間數(shù)據(jù)模型，為表達(dá)空間事物的多級層次關(guān)系提供一種良好的可視化手段。

1 定義

1.1 等級V圖定義[9]

定義1：Vh(P)={Vh(p1),Vh(p2),…,Vh(pn)}是以pi(i=1,2,…,n)為母點的Voronoi圖，h表示等級，h=1時等級最高，Vh(pi)表示母點pi的Voronoi區(qū)域(簡稱V區(qū)域)，令點集Qi={qj|qj∈Vh(pi),j≥1,1≤i≤n}，若Voronoi圖Vh+1(Qi)和Vh(pi)存有如下關(guān)系：

Vh+1(Qi)={Vh+1(qj)|Vh+1(qj)?Vh(pi),

qj∈Vh(pi),j≥1,1≤i≤n}

(1)

且有∑Vh+1(qj)≡Vh(pi)，則稱Vh(P)為等級Voronoi圖。公式為：

Vh(P)={Vh+1(Q1),Vh+1(Q2),…,Vh+1(Qn)}

(2)

1.2 等級加權(quán)V圖定義

定義2：Vh(P,λ)={Vh(p1,λ1),Vh(p2,λ2),…,Vh(pn,λn)}是以pi(i=1,2,…,n)為母點、λi為權(quán)重的加權(quán)Voronoi圖，Vh(pi,λi)表示母點pi的加權(quán)Voronoi區(qū)域(簡稱加權(quán)V區(qū)域)，令點集Qi={qj|qj∈Vh(pi,λi),j≥1,1≤i≤n}，若加權(quán)Voronoi圖Vh+1(Qi,λ)和Vh(pi,λi)存在如下關(guān)系：

Vh+1(Qi,λ)={Vh+1(qj,λj)|Vh+1(qj,λj)?Vh(pi,λi),

qj∈Vh(pi,λi),j≥1,1≤i≤n}

(3)

且有∑Vh+1(qj,λj)≡Vh(pi,λi)，則稱Vh(P,λ)為等級加權(quán)Voronoi圖。公式為：Vh(P,λ)={Vh+1(Q1,λ),Vh+1(Q2,λ),…,Vh+1(Qn,λ)}

(4)

為方便后續(xù)描述，均用Vh表示第h級V圖或加權(quán)V圖，Vh(i)表示母點pi的V區(qū)域或加權(quán)V區(qū)域。

2 等級V圖、加權(quán)V圖ArcGIS矢量生成算法

2.1 算法基本思想

以等級加權(quán)V圖為例，設(shè)P={p1,p2,…,pn}是二維歐氏空間上的n個母點，λ為母點權(quán)重，h為母點等級，相同h值的母點處于同一等級。算法基本思想：生成點集最小外包矩形，用h=1的母點集合劃分最小外包矩形，生成加權(quán)V圖V1；判斷h=2的點集歸屬，并利用該點集劃分V1，生成第二層的加權(quán)V圖V2，以此類推，直至算法結(jié)束。

加權(quán)V圖生成方法為：利用增量法思想，每插入一個點，對原有加權(quán)區(qū)域進行重新劃分，分割出屬于新插入點的加權(quán)V區(qū)域，當(dāng)所有點插入完成后，即得到最終的加權(quán)V圖。等級V圖是等級加權(quán)V圖權(quán)重相同時的特例，二者算法實現(xiàn)思想相同。

2.2 算法實現(xiàn)步驟

Input:P={p1,p2,…,pn}，pi∈P，包含屬性：ID，pi的唯一標(biāo)識(對應(yīng)圖1中的Belong字段)；Weight，pi權(quán)重(對應(yīng)圖1中的Weight字段)，與變量λ意義相同；Level，pi等級(對應(yīng)圖1中的Level字段)，與變量h意義相同。

Output:等級V圖、等級加權(quán)V圖。

圖1 ArcGIS點集示例

Fig.1 A set of points saved in ArcGIS attribute table

算法具體步驟如下：

(1)數(shù)據(jù)獲取及預(yù)處理。獲取母點集合P，并按其Level值劃分為不同等級的母點集合pL={pL1,pL2,…,pLm}，pLh∈pL，∑pLh≡P，pL1表示最高層次點集，pLm為最低層次點集。

(2)獲取P的最小外包矩形MBR。初始化加權(quán)V圖結(jié)構(gòu)體V0，并令V0=MBR。

(3)遞歸循環(huán)pLh集合，生成等級加權(quán)V圖，偽代碼如下：

for(inth=1;h≤m;h++){

//獲取當(dāng)前等級點集

pLh=GetLevelPoints(P,h)

//將pLh歸屬不同加權(quán)V區(qū)域

Vh-1=SetPnts2BelongRegion(pLh,Vh-1)

//初始化加權(quán)V圖結(jié)構(gòu)體Vh,_Vh

for(intj=1;j≤VorCount;j++){

//VorCount：當(dāng)前加權(quán)V區(qū)域個數(shù)

Vh←Vh-1(j)

for(intt=2;t≤VorPntsCount;t++){

//VorPntsCount：加權(quán)V區(qū)域內(nèi)包含的點個數(shù)

for(intk=1;k

if(pt.Weight=pk.Weight)

//垂直平分線劃分區(qū)域

Vh(tk)=DivideRegionWithPB(pt,pk,Vh(k))

else//圓弧劃分區(qū)域

Vh(tk)=DivideRegionWithArc(pt,pk,Vh(k)) }

//區(qū)域合并

Vh(tK)={Vh(t1),Vh(t2),…,Vh(tk),…,Vh(tt-1)}

Vh(t)=UnionRegion(Vh(tK))}

Vh←Vh

_Vh=null}

Vh=_Vh

//輸出h級加權(quán)V圖

ExportLevelVorRegion(Vh)}

等級V圖偽代碼和加權(quán)V圖類似，由于不需要判斷權(quán)重，所以只需保留垂直平分線劃分區(qū)域的方法，即代碼改變?nèi)缦拢?/p>

for(intk=1;k

Vh(tk)=DivideRegionWithPB(pt,pk,Vh(k))}

其他代碼不變。

2.3 算法說明

2.3.1 點集的區(qū)域歸屬判定 第一次循環(huán)，SetPnts2BelongRegion()方法認(rèn)定MBR為點集pL1的V區(qū)域或加權(quán)V區(qū)域V0，利用pL1劃分V0生成V1，自第二次循環(huán)起，將pLh中的點歸屬于Vh-1(j)(j=1,2,…,VorCount)。點集歸屬判別方法：如果pLh中的一個點p∩Vh-1(j)≠Ф，則認(rèn)為點p屬于Vh-1(j)。

2.3.2 區(qū)域分割和合并 假定pi(xi,yi)、pj(xj,yj)為空間中兩個不重疊的點，λi、λj是兩點權(quán)重。根據(jù)文獻(xiàn)[10,11]給出的V圖及加權(quán)V圖定義可知，兩點間Voronoi邊為一直線，其公式如下：

(5)

兩點間加權(quán)Voronoi邊為一圓弧，其公式如下：

(6)

圓心為：

(7)

半徑為：

(8)

式中:d(pi,pj)表示pi和pj間的歐氏距離。

以此數(shù)學(xué)理論為基礎(chǔ)，便可利用AE接口實現(xiàn)點集的加權(quán)V區(qū)域劃分，具體包括4個要點：直線生成、圓弧生成、V區(qū)域及加權(quán)V區(qū)域分割、V區(qū)域及加權(quán)V區(qū)域合并。

(1)直線生成。AE提供ILine接口，包含F(xiàn)romPoint和ToPoint屬性用于生成直線，表示通過設(shè)定起始點和終止點坐標(biāo)(通過式(5)和MBR邊界求交獲得)，生成一條兩點間的線段。該方法的具體實現(xiàn)包含在DivideRegionWithPB()方法中。

(2)弧段生成。IConstructCircularArc.ConstructCircle(IPointCenterPoint,doubleRadius)方法可用于生成弧段，其中圓心參數(shù)CenterPoint由式(7)獲得，半徑參數(shù)Radius由式(8)獲得。該方法的具體實現(xiàn)包含在DivideRegionWithArc()方法中。

(3)區(qū)域分割。如圖2所示，假設(shè)Vh(k)是點pk的加權(quán)V區(qū)域，pt是插入點，Lkt是pk和pt間的加權(quán)Voronoi邊，利用AE接口ITopologicalOperator.Cut(IPolylinecutter,refIGeometryleftGeom,refIGeometryrightGeom)方法，可將區(qū)域Vh(k)分割為leftGeom和rightGeom兩部分，其中，cutter即為Lkt，rightGeom和leftGeom對應(yīng)Vh(k) (分割后屬于pk的加權(quán)V區(qū)域)和Vh(tk) (從Vh(k)分割出來屬于pt的區(qū)域)。在AE提供的拓?fù)浞指罱涌谥校瑓?shù)Cutter即為由ILine和IConstructCircularArc接口生成的直線及弧段。該方法的具體實現(xiàn)包含在DivideRegionWithPB()和DivideRegionWithArc()中。

圖2 區(qū)域分割

Fig.2Regiondivisionbyanarc

加權(quán)V區(qū)域分割時，由于考慮權(quán)重因素，需判斷權(quán)重相等及不相等兩種情況，從而確定直線分割區(qū)域和弧段分割區(qū)域兩種方法。由于V區(qū)域劃分不涉及權(quán)重因素，所以只需考慮直線分割區(qū)域的情況。

(4)區(qū)域合并。UnionRegion()方法指定V區(qū)域和加權(quán)V區(qū)域合并方式。如圖3，設(shè)pk是h級加權(quán)V區(qū)域Vh(k) 所包含的點，k=1,2,…,c,c=VorPntsCount，若新插入一點pt(t=c+1)，其合并操作步驟如下：計算獲得pt與pk的加權(quán)Voronoi邊{L1t,L2t,…,Lkt,…,Lct}，利用Lkt對應(yīng)分割{Vh(1),Vh(2),…,Vh(c)}，將分割后屬于pt的區(qū)域{Vh(t1),Vh(t2),…,Vh(tc)}合并，生成pt的加權(quán)V區(qū)域Vh(t)。

圖3 區(qū)域合并

Fig.3 Regions union

2.3.3 等級加權(quán)V圖輸出 ExportLevelVorRegion()方法指定等級加權(quán)V圖的輸出方式：對于每一層級點集pLh，輸出對應(yīng)的加權(quán)V圖Vh，并以矢量圖形式存儲。將輸出圖層疊加，即可顯示最終的等級加權(quán)V圖效果。

3 應(yīng)用實例

圖4、圖5為等級V圖及加權(quán)V圖示例，其中母點集合為{pL1,pL2,pL3,pL4}，分別對應(yīng)于圖例中大圓、正方形、三角形、小圓所代表的點集，其中pL1包含4個母點，pL2包含10個母點，pL3包含11個母點，pL4包含125個母點。結(jié)果為四級等級V圖及加權(quán)圖{V1,V2,V3,V4}，邊界粗細(xì)程度代表等級程度差異，邊界越粗，等級越高。生成圖4所示的等級V圖平均用時5.21 s，生成圖5所示的等級加權(quán)V圖平均用時33.08 s。造成時間差異的主要原因是生成加權(quán)V圖時拓?fù)浞指詈秃喜⒉僮饕壬蒝圖復(fù)雜，所以用時較長。

4 結(jié)論

圖4 等級V圖 圖5 加權(quán)等級V圖

Fig.4 HVD Fig.5 HMWVD

本文在文獻(xiàn)[10]的基礎(chǔ)上首次實現(xiàn)了等級V圖及加權(quán)V圖的ArcGIS矢量生成算法，該算法有如下特點：1)提供了一種適用于地學(xué)領(lǐng)域的空間數(shù)據(jù)模型，應(yīng)用該算法可以方便地建立等級V圖及加權(quán)V圖，用以表達(dá)城鎮(zhèn)等級體系、機構(gòu)布局設(shè)置、區(qū)域?qū)哟尾町惖葍?nèi)容，豐富了Voronoi圖在地學(xué)領(lǐng)域的應(yīng)用。2)利用ArcGIS提供接口，以拓?fù)浞指詈喜⒎绞綄崿F(xiàn)算法，精度高，實用性強。生成的等級V圖及加權(quán)V圖和原始生成元文件可以組合存儲于GIS空間數(shù)據(jù)庫中，便于管理，同時可以靈活地更改點、面樣式，豐富輸出結(jié)果的表達(dá)方式。等級V圖及加權(quán)V圖以矢量方式生成，以方便與其他矢量或柵格數(shù)據(jù)進行疊加分析，挖掘更深層次的隱含信息。

[1] FARIN G,HAGEN H,HAMANN B.Hierarchical and Geometrical Methods in Scientific Visualization[M].Springer,2001.

[2] GOLD C,ANGEL P.Voronoi hierarchies[A].Proceedings of Geographic Information Science-4th International Conference,Geoscience[C].2006.99-111.

[3] WALLGRüN J O.Hierarchical Voronoi Graphs-Spatial Representation and Reasoning for Mobile Robots[D].Germany:University of Bremen,2008.

[4] 趙威，張仲元，秦耀辰,等.基于Voronoi圖的開封市域城鎮(zhèn)影像范圍和等級體系研究[J].河南大學(xué)學(xué)報(自然科學(xué)版)，2008，38(8)：271-275.

[5] 陳笑筑，王博，胡偉.基于Voronoi圖的鄉(xiāng)鎮(zhèn)等級體系劃分及影響范圍研究——以畢節(jié)市為例[J].安徽農(nóng)業(yè)科學(xué)，2011，39(35)：21947-21949.

[6] 薛瑩，俞路，韓貴峰.基于Voronoi圖剖分的區(qū)域旅游中心城市體系——以蘇浙滬地區(qū)為例[J].地域研究與開發(fā)，2008，27(4)：81-84.

[7] BALZER M,DEUSSEN O.Voronoi treemaps[A].Proceedings of the 2005 IEEE Symposium on Information Visualization[C].IEEE Press,2005.

[8] SUD A,FISHER D,LEE H P.Fast dynamic voronoi treemaps[A].Proceedings of the 2010 International Symposium on Voronoi Diagrams in Science and Engineering (ISVD′10)[C].2010.85-94.

[9] KANG S,LI J T,ZHANG F Q,et al.Hierarchical weighted Voronoi tessellations on visual application of spatial data[J].Journal of Computational Information Systems,2013,9(18):7341-7348.

[10] 范熙偉.加權(quán)Voronoi圖矢量生成算法研究及其實現(xiàn)[D].西安：西北大學(xué)，2011.21-35.

[11] 張有會.加權(quán)Voronoi圖的畫法研究[J].計算機科學(xué)，2001，28(6)：126-130.

Vector-Based Realization of Hierarchical Voronoi Diagram and Multiplicatively Weighted Voronoi Diagram with ArcGIS

TIAN Song1,CUI Xi-min1,SUN Yun-hua1,GONG Yu2

(1.CollegeofGeoscienceandSurveyingEngineering,ChinaUniversityofMiningandTechnology,Beijing100083; 2.SouthChinaSeaBranch,StateOceanicAdministration,Guangzhou510310,China)

Hierarchical Voronoi diagram (HVD) and hierarchical multiplicatively weighted Voronoi diagram (HMWVD) are methods of expressing hierarchical relationships in the form of trees,which are widely applied in the field of geography,such as urban hierarchies exhibition,organization structures expression,etc.A vector-based approach is first proposed to generate HVD and HMWVD through the methods of regions division and regions union with ArcGIS Engine based on the thought of incremental method.It provides an important way to research and develop GIS spatial data model.

hierarchical Voronoi diagram;hierarchical multiplicatively weighted Voronoi diagram;ArcGIS;region division;regions union

2014-05-13；

2014-07-08

國家自然科學(xué)基金項目(71150001)

田松(1985-)，男，博士研究生，研究方向為GIS算法研究。E-mail：tsisatc@163.com

10.3969/j.issn.1672-0504.2015.02.008

P208

A

1672-0504(2015)02-0034-04