一種基于柵格的加權Voronoi圖構建普適方法

劉 寶 舉，劉 慧 敏，鄧 敏，樊 子 德

(中南大學地球科學與信息物理學院地理信息系，湖南 長沙 410083)

一種基于柵格的加權Voronoi圖構建普適方法

劉 寶 舉，劉 慧 敏，鄧 敏，樊 子 德

(中南大學地球科學與信息物理學院地理信息系，湖南 長沙 410083)

Voronoi圖是空間分析的一個重要工具。該文將Voronoi區域視為流域，將柵格加權距離視為高程，提出一種顧及非空間屬性的能夠與ArcGIS無縫集成的Voronoi圖生成方法。首先，根據Voronoi圖的原始定義直接計算柵格的最小生成元加權距離，并仿D8算法思想，確定每個柵格的流向。然后，提取所有只有流出沒有流入的柵格，并對柵格邊界進行去噪處理和矢量化，得到Voronoi區域公共邊，并生成附生成元屬性的加權Voronoi圖。最后，基于ArcEngine實現了任意生成元的帶有非空間屬性的加權Voronoi圖。通過對比實驗表明，該文所提出的方法能夠高精度構建包含任意生成元的加權Voronoi圖。

加權Voronoi圖；加權距離柵格；流域；生成元

0 引言

Voronoi圖是由俄國數學家M.G.Voronoi于1908年提出的一種空間分割算法，表達了在空間分布的若干要素對空間的最短距離剖分。Voronoi圖最初用于解決計算幾何的凸殼、Delaunay三角剖分、最小生成樹等問題[1,2]，由于鄰近性、鄰接性等特殊性質，其應用迅速擴展到眾多領域。而在地理信息學科，Voronoi圖作為表達空間鄰近關系與空間劃分的有力工具，應用于地圖綜合[3,4]、空間鄰近分析[5,6]、空間聚類分析[7]、空間關聯規則挖掘[8]等方面。

從算法實現的數據模型角度，Voronoi圖生成方法可以分為基于矢量和基于柵格兩大類。矢量算法的提出比較早，算法也較多，大都是通過精確的數學計算得到垂直平分線，進而求取Voronoi圖，最為典型的是增量構造法[9,10]。這類算法可以得到精準的邊界線，時間和空間效率較高，但是其計算過程和數據結構復雜，并且不易于擴展到加權Voronoi圖的生成[11]。柵格算法是基于生長元柵格按一定的距離模板向外擴張最終形成Voronoi圖，具有代表性的是基于動態距離變換的方法[12,13]。以上算法都是從計算幾何的角度出發，生成的Voronoi圖不能很好地與GIS結合，限制了在地學領域的應用。近年來，Dong[14]和田松[15]等先后以開發組件的形式探討了Voronoi圖與GIS的集成。但文獻[14]方法僅能處理單一生成元，而且Voronoi圖有雜質面情況存在，文獻[15]方法未考慮生成元的非空間屬性，不利于Voronoi的空間分析。此外，ArcGIS軟件中Cost Distance和Cost Allocation工具可以實現加權Voronoi圖，但誤差較大，不能準確體現生成元權重與Voronoi圖的關系。為此，本文以shp格式矢量數據作為輸入輸出，研究了一種顧及非空間屬性并能與GIS軟件集成的高精度Voronoi圖構建方法。

1 基于加權距離柵格的Voronoi圖生成策略

利用加權距離柵格構建Voronoi圖，首先從單個生成元入手，通過計算每個生成元的加權距離柵格，根據匯流思想可提取Voronoi柵格區域并獲得柵格的生成元距離，同時保留生成元非空間屬性特征，使之便于空間分析處理。在此基礎上，通過Voronoi區域柵格邊界鋸齒消除算法，以提高生成Voronoi圖的精度。

設S={p1,p2,…,pn}為二維歐氏空間(平面)上的點集，Voronoi圖是根據最近距離原則對平面進行的一種分割，分割得到的每個區域內部包含一個點，并且區域內所有點到該點的距離均小于到其他點的距離，可表達為：

(1)

式中：d(p,pi)為點p與點pi之間的歐氏距離；點pi(i=1,2,…,n)稱為生成元。

根據式(1)，則可得到每個要素的空間影響區域，它是對空間的一個完整劃分，亦稱為普通Voronoi圖，如圖1a所示。這種普通的Voronoi圖不考慮屬性對距離的影響。在實際應用中，每個要素對周圍區域的影響程度通常不同。例如，同樣是購物點，大商場比小商鋪供貨能力強，相應的正常服務范圍更廣。為此，需要構建加權Voronoi圖，才能更好地體現每個要素在屬性上的差異。

加權Voronoi圖與普通Voronoi圖的區別是在距離的設置上增加了一個權值(圖1b)，可表達為：

(2)

式中：d(p,pi)為點p、pi之間的歐氏距離,i=1,2,…,n;λi為對應點元pi的權重。當λ1=λ2=…=λn時，則對應普通Voronoi圖。加權Voronoi圖可以很好地模擬不同級別地物的影響范圍。

圖1 普通Voronoi圖和加權Voronoi圖

Voronoi圖作為一種具有特定結構的矢量多邊形，公共邊上的點到相鄰生成元的距離或加權距離相等。因此，根據距離或加權距離確定Voronoi公共邊是構建Voronoi的一個主要途徑。距離生成元越遠的區域，受到生成元的影響越小，且Voronoi公共邊上的點受與之最鄰近的兩個生成元的影響相等，可見，Voronoi邊是不同地物影響能力的等值線。如果對生成元的影響進行量化表達，并用柵格數據表示這種影響值，則可以很好地表示這種影響隨距離增長而減小的趨勢，這種柵格數據即是加權距離柵格。如圖2所示，柵格距離值隨著生成元的影響能力而逐漸變化，鄰近兩點間影響能力相等的柵格點形成線，即是要提取的Voronoi公共邊。進而，加權Voronoi區域生成轉化為加權距離柵格構建與Voronoi區域提取。

圖2 生成元影響能力變化

2 加權距離柵格構建

2.1 單個生成元加權距離柵格計算

對區域進行柵格劃分，并用生成元的加權距離表示柵格值，即為加權柵格距離。為構建加權距離柵格，需要從單個生成元入手，計算每一個生成元pi(i=1,2,…,n)的加權距離柵格fi(i=1,2,…,n)，然后將所有生成元的加權距離柵格依據特定關系進行疊加，獲得所有生成元的加權距離柵格。

如圖3所示，分別表示點、線、面狀生成元影響下的加權距離變化。其中，xoy為柵格平面，z軸表示加權距離。加權距離的變化速度與生成元的權值直接相關，權值越大，傾角α越小，距離增長越緩慢，反之亦然。普通Voronoi圖中所有生成元對應的α角為45°。隨著權重的增大，傾角變小，生成元恒定距離以內的范圍擴大。用一過點A且平行于z軸的平面切圖3a所示的錐體，得到其縱切面圖，如圖4所示。圖中AK的投影AP表示加權生成元A的影響區域，AK上任意一點D表示生成元A在投影位置C處柵格的加權距離為D的z坐標值。AF對應普通生成元的距離柵格，角β為45°，點E的z坐標值對應普通生成元A在柵格C處的距離值。

圖3 單一生成元加權距離柵格示意

圖4 單一點生成元加權距離投影

設d((x-x1),(y-y1))為普通距離函數，f(x,y)為加權距離函數，生成元A的權值為ωα，則有：

(3)

式中：x、y分別表示柵格位置坐標，x1、y1分別表示生成元坐標。類似地，可計算其他所有單個生成元的柵格加權距離fi(x,y),i=(1,2,…,n)。這種加權距離采用矢量歐氏距離獲得，其精度始終保持在一個像元誤差，與傳統的依據距離模板擴散生長的模型相比精度更高，避免了因為方向不同產生距離不同和因為距離模板的選擇而導致精度不統一的情形，亦保證了后續生成加權Voronoi圖的精度。

一般地，要素的幾何類型包括點狀、線狀和面狀，對于線狀要素，通常將其離散化為多個點進行處理；對于面狀要素，首先設置其內部點的加權距離為0，然后將其邊界離散化為多個點進行處理。

2.2 多生成元加權距離柵格確定

Voronoi圖是多個生成元對區域的分割，如圖5所示，多個生成元并存，在Voronoi圖對區域的劃分中，一個Voronoi區域的內部有且僅有一個對應的生成元，而Voronoi公共邊處其所指向的相鄰生成元加權距離相等。因此，在計算所有單個生成元的柵格距離值之后，需要確定最小的加權距離，即：

F(x,y)=min(fi(x,y)),i=(1,2,…,n)

(4)

圖5 多個生成元距離柵格示意

圖6 加權距離柵格示意

3 加權Voronoi區域提取

3.1 區域柵格劃分

如圖7所示，將柵格加權距離視為高程[16]，由于生成元處的柵格距離值為極小值0，因此，水流匯聚點即為生成元，而匯聚到同一生成元的柵格區域對應該生成元的流域，即Voronoi區域，而相鄰流域的分水嶺即為Voronoi區域邊界，從而，提取出流域并將其轉換為矢量面就可以獲得相應的Voronoi圖。為提取流域信息，需要確定各位置的水流方向，即每個柵格像元水流流出的方向，進而獲得水流匯聚情況。

圖7 加權距離柵格

引入D8算法思想，按八鄰域劃分方向，并按照2n(n=0,1,…，7)對八鄰域進行編碼，如圖8所示。根據每個加權距離柵格像元與周圍八鄰域像元的距離梯度，確定最陡下降方向，作為該柵格的流向，并用流向編碼中表示該方向的值對輸出像元進行編碼，由此，得到流向柵格。然后，從生成元柵格開始逆向搜索匯入柵格，歸屬為對應的流域，直至為每個柵格分配唯一流域，即可得到流域柵格，完成Voronoi區域劃分。

圖8 流向編碼

Voronoi圖是空間分析的重要工具，空間信息和非空間屬性信息是空間分析的兩個基礎，附屬非空間屬性的Voronoi圖更利于空間分析。為此，進行Voronoi區域劃分后，為每一個Voronoi面賦予對應生成元的非空間屬性，生成帶屬性值的矢量面結構。

3.2 邊界處理

在對每一個柵格確定流域的過程中，由于某些位于分水嶺上的特別是多流域交叉點處的柵格的流域歸屬無法準確確定，所以在將柵格區域轉換為矢量格式過程中，某些結點或邊界處會生成如圖9a所示的雜質面，這種雜質面只有一個像元大小，但形成了區域空洞，同時影響基于Voronoi圖的統計分析，因此，有必要在轉換為Voronoi區域之前對流域柵格邊界進行雜質面去除處理。

分析發現，在流域分區柵格中，當邊界上的中心像元與上下左右四鄰域像元值均不相同時，會產生雜質面，按照四鄰域像元值是否相同分為4種情形，即四鄰域分別取2種、3種和4種值的情況，如圖10所示。對于四鄰域對向相同的情形，如圖10e所示，表示匯水流域或Voronoi區域交叉，與實際情形不符，不予考慮。針對圖10a-圖10d所示的四鄰域模板及其旋轉形式，用鄰域柵格代替中心像元，其中：圖10b與圖10c所示情形用a代替，圖10a所示情形用a或b代替，圖10d所示情形用a、b、c、d任一值代替，使得Voronoi圖邊界存在一個像元的誤差。經過代替處理后，轉換得到的矢量Voronoi圖不包含雜質面，鋸齒效應亦得到明顯改善，如圖9b所示。

圖9 流域柵格處理前后生成的Voronoi圖

圖10 中心像元與四鄰域值不同的情形

由以上步驟，通過對研究區域的柵格劃分，生成任意生成元的加權距離柵格，根據流域思想提取Voronoi邊界，并對Voronoi圖進行邊界處理，生成加權Voronoi圖，最后將生成元專題屬性關聯到Voronoi圖中。這種帶有非空間屬性的Voronoi圖可以直接用于空間分析處理，擴展了其應用。

4 實驗與分析

4.1 方法驗證

為了驗證本文方法的有效性，首先用一組模擬的任意生成元進行實驗。如圖11所示，該組生成元包括點、線、面3類不同要素，并被賦予不同的權值。采用本文方法，構建其加權Voronoi圖(圖11b)，并輸出加權距離柵格(圖11c)。從圖中可以看出，本文方法能夠生成不同類別生成元的加權Voronoi圖，然而，劃分柵格的大小不可避免地影響Voronoi圖的精度并導致邊界的鋸齒效應。在計算效率方面，由于涉及處理過程較多無法準確給出時間復雜度，但通過實驗證明在250*250的柵格大小情況下，構建包含50個生成元的Voronoi圖所耗費時間為50 s左右(i3，2.27GHz，2GB RAM)，隨著柵格分辨率的變大，耗費時間也會有所增加。

圖11 不同生成元加權Voronoi圖

為定量分析本文方法的精度效果，設計了第二組實驗來檢驗本文方法的精度。如圖12所示，采用6組點狀分布生成元數據，每組數據包含14個隨機點，用本文方法構建其等權Voronoi圖，將其與普通Voronoi圖進行對比，選取對應的Voronoi區域的面積誤差百分比作為精度評價指標，對本文方法的精度進行檢驗。表1列出了由6組實驗數組計算得到的平均誤差與最大誤差。從表中可以看出，本文方法的平均誤差均在1%以下，具有較高的精度。本文方法在計算柵格距離時，直接采用普通距離作為基準，再進行加權掩膜，與傳統的依據距離模板擴散生長的模型相比，避免了柵格距離與生成元方向相關而導致精度不統一的問題，從而提高了最終Voronoi圖的精度。需要補充說明的是，可結合應用需要設定加權距離柵格的大小，精度要求越高，設定柵格越小，分辨率的提高還可以有效地避免多生成元共面的情形。

圖12 6組隨機數據生成的普通Voronoi圖

表1 面積誤差率

Table 1 Area error rate

第一組第二組第三組第四組第五組第六組最大誤差1.788%0.893%1.044%1.558%1.202%0.887%平均誤差0.590%0.352%0.462%0.579%0.421%0.336%

4.2 對比分析

ArcGIS中的Cost Distance和Cost Allocation工具也可以實現加權Voronoi圖，為了比較其與本文方法的差異，采用圖1所示數據進行對比實驗，生成結果如圖13所示，上標數字表示權重。由于加權距離柵格本身就體現了生成元的影響范圍，所以本文方法直接提取加權距離柵格中“分水嶺”作為Voronoi邊界，而ArcGIS方法通過統計加權距離柵格中各像元到鄰近生成元的最小加權距離以確定像元歸屬，即Voronoi邊界上像元到鄰近兩個生成元的加權距離相等，由于生成元權重和影響范圍呈正比，因此這種累加距離必然導致Voronoi邊相對靠近權值較大的生成元，弱化了相對重要生成元的影響力。所以當相鄰生成元權重相差比較大時，本文方法生成的Voronoi圖能夠更加準確地體現不同權重生成元的影響區域。

圖13 不同方法生成加權Voronoi圖

Voronoi圖與GIS軟件的結合是Voronoi圖生成方法廣泛推廣應用的重要前提，本文基于ArcEngine實現了能夠與GIS軟件完全集成的Voronoi構建過程。與文獻[14]的方法相比，本文方法可以生成任意生成元及其組合的加權Voronoi圖，與文獻[14]的方法僅能處理單一生成元的情況。為了驗證兩種方法在實際應用中的效果，采用2013年中國西部七省3.5級以上地震分布點數據作為驗證數據，震級作為生成元權重。結果如圖14所示，兩種方法生成的加權Voronoi圖大致相同，但是在細節處理上文獻[14]方法生成的Voronoi圖在部分邊界上有雜質面產生，需要采用本文邊界處理方法對其進一步優化處理方可得到平滑邊界。

圖14 中國西部2013年3.5級以上地震分布Voronoi圖

5 結論

本文結合Voronoi圖的結構特性，從其基本定義出發，利用普通距離柵格與生成元權值掩膜運算構建加權距離柵格，進而提出了一種加權Voronoi圖生成的普適方法，避免了在歐氏距離下擴散模板的選擇導致的精度不統一問題[13]，使整個Voronoi圖保持同一精度，并有效解決了許多柵格方法存在的鋸齒問題。同時，此方法可以與GIS軟件完全集成，顧及非空間屬性和伴生加權距離柵格的特性使得該方法得到的結果可直接用于GIS的空間分析處理。實驗證明，與現有方法相比，本文方法精度高，能夠處理任意種類生成元，并且與GIS軟件無縫集成，同時有效避免了矢量算法的結構復雜性。另外，該方法產生的任意中間數據和最終生成的加權距離柵格、Voronoi圖都可以直接作為軟件輸入數據進行后續處理。然而，對于各生成元的權重差異很大的情形，還不能較好地生成多孔洞形態的Voronoi圖;其次，算法的效率依然較為低下，生成效率和柵格精度難以兼顧，今后需進一步研究其效率優化問題，使其能在精度和效率方面得到共同提高。

[1] OKABE A,BOOTS B,SUGIHARA K,et al.Spatial Tessellations:Concepts and Applications of Voronoi Diagrams[M].2nd,New York:Wiley Press,2000.

[2] DU Q,FABER V,GUNZBURGER M.Centroidal Voronoi tessellations:Applications and algorithms[J].SIAM Review,1999,41(4):637-676.

[3] 閆浩文,王邦松.地圖點群綜合的加權 Voronoi 算法[J].武漢大學學報(信息科學版),2013,38(9):1088-1090.

[4] 李佳田,康順,羅富麗.利用層次 Voronoi 圖進行點群綜合[J].測繪學報,2014,43(12):1300-1306.

[5] LI Z L,HUANG P Z.Quantitative measures for spatial information of maps[J].International Journal of Geographical Information Science,2002,16(7):699-709.

[6] 劉慧敏,鄧敏,樊子德,等.地圖上居民地空間信息的特征度量法[J].測繪學報,2014,43 (10):1092-1098.

[7] YAN H W,WEIBEL R.An algorithm for point cluster generalization based on the Voronoi diagram[J].Computers & Geosciences,2008,34(8):939-954.

[8] 李光強,鄧敏,朱建軍.基于 Voronoi 圖的空間關聯規則挖掘方法研究[J].武漢大學學報(信息科學版),2009,33(12):1242-1245.

[9] 張有會.加權 Voronoi 圖畫法的研究[J].計算機科學,2001,28(6):126-128.

[10] GONG Y,LI G,TIAN Y,et al.A vector-based algorithm to generate and update multiplicatively weighted Voronoi diagrams for points,polylines,and polygons[J].Computers & Geosciences,2012,42:118-125.

[11] LETSCHER D.Vector Weighted Voronoi Diagrams and Delaunay Triangulations[C].Proceedings of the 19th Annual Canadian Conference on Computational Geometry.2007.

[12] 李成名,陳軍.Voronoi 圖生成的柵格算法[J].武漢測繪科技大學學報,1998,23(3):208-210.

[13] CHEN J.A raster-based method for computing Voronoi diagrams of spatial objects using dynamic distance transformation[J].International Journal of Geographical Information Science,1999,13(3):209-225.

[14] DONG P.Generating and updating multiplicatively weighted Voronoi diagrams for point,line and polygon features in GIS[J].Computers & Geosciences,2008,34(4):411-421.

[15] 田松,崔希民,孫云華,等.顧及障礙物的一般圖形 Voronoi 圖及其加權圖的 ArcGIS 柵格實現[J].地理與地理信息科學,2014,30(4):42-45.

[16] 艾廷華,禹文豪.水流擴展思想的網絡空間Voronoi圖生成[J].測繪學報,2013,42(5):760-766.

A General Raster-Based Approach for Generating Weighted Voronoi Diagrams

LIU Bao-ju,LIU Hui-min,DENG Min,FAN Zi-de

(DepartmentofGeographicInformatics,CentralSouthUniversity,Changsha410083,China)

Voronoi diagram plays a very important role in spatial analysis.Existing methods of generating Voronoi diagrams seldom consider the effect of non-spatial attributes,in this case the partition of geographical space is inconsistent with human cognition.To solve this problem,a unified raster-based algorithm is presented for generating weighted Voronoi diagrams in this paper.First,the region is divided into a set of small regular grids,and the weighted distance of each grid generator is calculated based on the original definition of Voronoi diagram.Second,flow direction of each grid is calculated by using an eight direction (abbreviated as D8) flow model,where the weighted distance is regarded as the elevation and the Voronoi region as being a basin.In this situation,the collection of grids which has no any inflow consists of the extract common edges of adjacent Voronoi regions.The weighted Voronoi diagrams are further built with a modification of the common edges by considering non-spatial attributes,which is implemented by means of ArcEngine software.Finally,practical examples are provided to illustrate the advantages of the proposed method in this paper,compared with current methods.

weighted Voronoi diagram;weighted distance grid;watershed;generator

2015-11-23；

2016-04-05

國家自然科學基金項目(41301515)；地理空間信息工程國家測繪地理信息局重點實驗室經費資助項目(201413)；中國博士后科學基金(2014M562134)

劉寶舉(1989-)，男，碩士研究生，主要從事GIS空間分析研究。*通訊作者 E-mail：lhmgis@163.com

10.3969/j.issn.1672-0504.2016.04.004

P208

A

1672-0504(2016)04-0017-06