面向雷利法則的加權Voronoi圖生成方法

康　順，瞿珊珊

(1. 中國礦業大學(北京)地球科學與測繪工程學院，北京 100083； 2. 中南大學地球科學與信息物理學院，湖南 長沙 410083)

Voronoi雛形可見于笛卡爾的《哲學原理》，文中描述了在可視宇宙中太陽系是由許多漩渦(Vortices)組成的這一理念[1]。在渦旋空間被抽取為凸域集合概念之后，這一概念在不同學科中被賦予了種種名稱：生物學的中軸變換(medial axis transform)；化學、物理學的維格那-塞茨原胞(Wigner-Seitz zones)；晶體學的作用域(domains of action)；氣象學、幾何學的泰森多邊形(Thiessen polygons)等。最后，由數學家Dirichlet和Voronoi將這一概念正式命名為狄利克雷鑲嵌(Dirichlet tessellation)或沃洛諾伊圖(Voronoi diagram)[2]。作為計算幾何的重要內容，Voronoi圖是解決骨架線提取、凸包計算、Delaunay三角化，以及數據可視化等問題的有力工具[3]。因此，研究Voronoi圖的高效生成算法具有重要的現實意義。

在一般Voronoi圖與廣義Voronoi圖的生成方法上已有諸多學者進行了相關研究，主要分為矢量生成方式和柵格生成方式兩種。基于矢量的Voronoi圖生成方法主要有增量法、分治法、間接法，由于矢量數據的存儲量大、存儲結構比較復雜，對點線面復合目標構成的空間實體進行Voronoi計算時，多要求生長元是點或半線，對線或面要素則往往分解為點要素來構建Voronoi圖，該方式破壞了圖形的完整性，并增加了計算的復雜度[4-6]。由于矢量方法的局限性，基于柵格的Voronoi圖生成方法得到研究與發展。柵格Voronoi圖的生成方式一是依據傳統距離變換的柵格Voronoi生成[7-10]，在基于空間生長目標圖像轉變為距離圖像的柵格Voronoi圖生成方法中定義柵格距離是關鍵，已有的主要距離形式有歐氏距離、曼哈頓距離(或L1距離)、切比雪夫距離(或L∞距離)及其加權形式；二是基于形態學的膨脹算子生成柵格Voronoi圖[11]。在提高柵格Voronoi圖生成的效率方面，柵格Voronoi生成要計算每個柵格與各個發生元間的距離，因此當發生元所占柵格數量很多時，往往存在較高的時間計算復雜度。對此，文獻[12]利用四叉樹和區間算術，規避對背景柵格的逐項處理以提高柵格計算速度；文獻[13]通過查找鄰近柵格，計算背景柵格與過濾后的生長元之間的距離來判斷柵格歸屬提高Voronoi生成效率，以及依賴處理器物理性能的Voronoi并行計算方法[14]。

雖然基于柵格的Voronoi生成方法能夠處理復雜的空間實體，對生長元的約束寬松，但算法效率和計算精度是一直以來所研究的問題。通常，基于柵格生成的Voronoi圖精度低、計算量大、耗時長，主要歸因于在距離圖上需要進行多次計算、比較每一個點與每一個生長元的距離來確定歸屬，已有提高柵格Voronoi圖生成效率的方法主要采用鄰近方式過濾出與背景柵格鄰近的發生元柵格，從而減少計算次數、優化柵格Voronoi圖生成效率，柵格處理過程中各個柵格的計算并不依賴于其他柵格，即是相互獨立的。因此，從規避背景柵格與生長元柵格的距離計算方面，本文提出了基于地圖代數的加權Voronoi圖，即雷利Voronoi圖(Reilly Voronoi diagram，RVD)生成方法。該方法利用雷利法則的綜合權重生成每個柵格生長元的距離圖，對由每個生長元生成的距離圖進行地圖代數計算，最后實現空間場的柵格加權Voronoi圖剖分。

1　Voronoi圖定義

1.1　一般Voronoi圖

一般Voronoi圖，又稱普通Voronoi圖，或常規Voronoi圖，或1-site Voronoi圖。以空間R2點要素為例，若存在空間點狀生長元集合P={p1，p2，…，pn}(nN*)，對空間任一點q，若滿足如下關系

d(pi，q)≤d(pj，q)

(1)

Vor(P)={Vor(p1)，Vor(p2)，…，Vor(pn)}

(2)

圖1　一般Voronoi圖

1.2　廣義Voronoi圖

當生長元被賦予一定權重時，生成的一般Voronoi圖擴展為廣義Voronoi圖。以空間R2點要素為例，若存在空間生長元集合P={p1，p2，…，pn}，相應的權重W={w1，w2，…，wn}(nN*)，對空間任一點q，若滿足如下關系

d(pi，q)/wi≤d(pj，q)/wj

(3)

則稱由滿足條件的空間點要素所構成的區域為生長元pi的廣義Voronoi圖。生長元的權重越大，其廣義Voronoi圖就越大，所有生長元P則將該空間剖分為各自的加權Voronoi勢力范圍Vor(·)，示意圖如圖2所示。

圖2　加權Voronoi圖

(4)

根據生長元個數及生成的空間上下文不同，廣義Voronoi圖又衍生出集群性生長元作用下的k階Voronoi圖[15]、受道路網影響的變速城市Voronoi圖[16]，網絡Voronoi圖[17-18]、根據方向比例計算生成的維度比率Voronoi圖[19]、依據角度計算生成的角度Voronoi圖[20]等形式。

2　Voronoi地圖代數計算

2.1　雷利距離變換

雷利零售引力法則(Reilly’s law of retail gravitation)是由美國學者威廉·J·雷利根據牛頓力學的萬有引力理論，利用3年時間調查了美國150個城市對周圍地區的吸引力，總結了都市人口與零售引力的相互關系，提出了“零售引力規律”，被稱為雷利法則或雷利零售引力法則。該法則描述了一個城市的吸引力與它的規模成正比，與它們之間的距離成反比，用以解釋城市規模與建立的商品零售區之間的關系[21]，如圖3、圖4所示的不同商圈及不同的城區規模對顧客數量、空間引力范圍的描繪。

圖3　不同商圈的顧客引力

圖4　城區規模引力

城區規模不同產生的空間引力范圍不同，這正是傳統柵格加權Voronoi圖的權重局限性，即在計算上往往注重生長元的權重大小，即規模權重，而忽略了隨距離而產生的距離變換權重的影響。因此，依據雷利法則研究將Voronoi圖的距離約束重新界

定為

RF(·)=f(d(·))/w

(5)

式中，f(·)函數指值域隨著距離變量的增大而增大，空間目標距離生長源的距離呈現出歐氏距離變換的漸變規則，如f(d)=d+1、f(d)=d2等。盡管Voronoi圖的邊界復雜，難以計算，但是基于柵格Voronoi圖的計算則相對易于實現。

2.2　地圖代數計算

將連續空間離散化的柵格數據是場模型在信息系統中的具體實現。柵格的“位”數據蘊含了豐富的關系數據，更加適用于網絡分析，在“數字地球”等大型地理信息工程建設的空間分析中，柵格數據展現了明顯的優勢[22]。地圖代數研究的正是度量空間下的抽象點狀要素目標集，即相應空間位置上柵格值的代數計算。與代數運算所不同的是，代數運算只有矩陣的加法和減法，是相應位置上的計算。即地圖代數遵循了位置一一對應的轉換規則，是空間計算的重要工具，如圖5、圖6所示的兩個圖幅的地圖代數乘法計算。以圖5的柵格1與柵格2的第一行為例，1×2=2，1×2=2，3×3=9，得到輸出柵格結果。多圖幅的最低位置計算輸出結果為具有最小值的像素的位置，即第幾個柵格，如圖6的柵格1、柵格2與柵格3的第一行為例，像素值1 ￡0 ￡Null=Null，1 ￡1 ￡1=1(柵格1),0 ￡1 ￡0=1(柵格1)，0 ￡1 ￡0=1(柵格1)，0 ￡0 ￡0=1(柵格1)。

圖5　地圖代數乘法運算

圖6　最低位置算子￡

結合雷利法則，本文利用地圖代數的算術計算、關系計算，設計的地圖代數生成柵格加權Voronoi圖(RVD)的方法描述如下，生產的流程如圖7所示。

圖7　地圖代數的雷利Voronoi圖生成過程

雷利Voronoi圖的地圖代數生成方法為:

輸入：空間點要素集合PW。

輸出：PW的柵格加權Voronoi圖。

(2) 依據雷利距離變換RT(·)，地圖代數算術重計算柵格距離，生成新的柵格距離變換圖[NRDi]←RT(RDi)。

(3) 地圖代數關系計算所有柵格距離圖[NRDi]中的局部最小值←min([NRDi])。

(4) 地圖代數關系計算比較minval和[NRDi]，如果minval≥[NRDi]，則將柵格R歸屬為第一個小于minval的生長元p，Rp￡(minval，[NRDi])。

(5) 算法結束。

3　算例與分析

研究的試驗算例以遼寧省的沈陽市、撫順市和鐵嶺市3個市區為研究案例，研究區示意圖如圖8所示，研究提出的地圖代數雷利Voronoi圖生成方法利用Python 2.5.1+ArcScripting開發實現。在雷利Voronoi圖的權重界定上，本試驗采用f(d)=d2的歐氏距離重計算變換方式；規模權重設置根據全國第六次人口普查獲取的各市區人口統計數據作為其影響的重要性指標。從全國第六次人口普查(http:∥www.stats.gov.cn/ztjc/zdtjgz/zgrkpc/dlcrkpc/)公布的數據中得出沈陽8 106 171人，撫順2 138 090人，鐵嶺2 717 732人，以此作為生長元的規模權重。

為了規避在單幅圖中背景柵格與每一生長元的距離計算，試驗首先對每一生長元生成其歐氏距離變換圖，如圖11所示的撫順市、沈陽市、鐵嶺市的歐氏柵格距離變換圖。在此距離變換基礎上，利用最

低位置算子￡實現其一般Voronoi圖，如圖9所示。從各市區的一般Voronoi圖可看出，各市區的作用范圍與其行政區劃面積大體近似，而沈陽作為省會城市其影響范圍應大于周邊市區，從商圈與零售引力上，普通Voronoi與市區的實際影響力不符。

圖8　試驗區域

圖9　市區的一般Voronoi

圖10　市區的雷利Voronoi圖

在規模權重與距離變化權重的作用下，經對各市區的歐氏距離圖變換與最低位置算子￡計算(如圖12所示)，實現了各市區的雷利Voronoi圖(如圖10 所示)。對比圖9，從圖10中可以看出作為省會的沈陽市受綜合權重的影響具有更大的Voronoi作用域，遠超出了其行政區劃的作用范圍，更符合省會城市的實際效應。此外，以此方式生成的雷利Voronoi圖并未涉及背景柵格與每一個生長元的距離計算，通過簡單的地圖代數計算降低了計算的復雜性，規避了單幅圖中由多次計算、歸屬判斷所帶來的時間復雜度。

圖12　各市區的雷利變換距離

4　結　語

針對傳統柵格Voronoi圖生成方式中僅顧及生長元的權重局限性、單圖幅中計算、判斷背景柵格與每一生長元的柵格距離計算復雜性，研究基于雷利法則，從生長元的規模權重和距離變換權重出發，依據每個生長元獨立生成的柵格歐氏距離圖，利用地圖代數的方法重計算每個生長元的權重距離圖，最后，在最低位置算子作用下，對所有生長元新產生的權重距離圖剖分各自的競爭空間，即雷利Voronoi圖。試驗證明，該方法產生的雷利Voronoi圖比一般Voronoi圖更符合市區對空間的影響范圍，體現出了省會市區在競爭域上明顯大于其行政區劃范圍。該方法不僅完善了權重影響因素，而且結合每個生長元的獨立柵格距離圖，從多個距離圖幅角度，以地圖代數的形式規避了在單一圖幅中遍歷柵格，逐一處理背景柵格所屬的生長元而造成的計算效率問題；從柵格數據的精度角度，如何根據實際需求設定柵格的大小以滿足Voronoi在實際中的應用是下一步工作需要繼續探討和解決的問題。

參考文獻：

[1] DESCARTES R，MILLER V R，MILLER R P.Principles of Philosophy[M].Boston：[s.n.]，1908.

[2] VORONOI G.Nouvelles Applications des Paramètres Continusàla Théorie des Formes Quadratiques[J].Journal für die Reine und Angewandte Mathematik，1908(134)：198-287.

[3] 劉金義，劉爽.Voronoi圖應用綜述[J].圖學學報，2004，25(2)：125-132.

[4] 康順，相詩堯.基于馬氏距離的多高斯Voronoi圖生成方法[J].地理與地理信息科學，2016，32(3)：49-52.

[5] 李佳田，楊琪莉，羅富麗，等.線/面Voronoi圖的分解合并生成算法[J].武漢大學學報(信息科學版)，2015，40(11)：1545-1550.

[6] 王新生，劉紀遠，莊大方，等.基于GIS的任意發生元Voronoi圖逼近方法[J].地理科學進展，2004，23(4)：97-102.

[7] 劉寶舉，劉慧敏，鄧敏，等.一種基于柵格的加權Voronoi圖構建普適方法[J].地理與地理信息科學，2016，32(4)：17-22.

[8] 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.

[9] BAREQUET G，DICKERSON M T，SCOT R L.2-Point Site Voronoi Diagrams[J].Discrete Applied Mathematics，2002，122(1)：37-54.

[10]LI C，CHEN J，LI Z.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.

[11]李佳田，陳軍，趙仁亮，等.基于線性四叉樹結構的Voronoi圖反向膨脹生成方法[J].測繪學報，2008，37(2)：236-242.

[12]壽華好，袁子薇，繆永偉，等.一種平面點集Voronoi圖的細分算法[J].圖學學報，2013，34(2)：1-6.

[13]王新生，劉紀遠，莊大方，等.一種新的構建Voronoi圖的柵格方法[J].中國礦業大學學報，2003，32(3)：293-296.

[14]王結臣，蒲英霞，崔璨，等.一種基于點集自適應分組構建Voronoi圖的并行算法[J].圖學學報，2012，33(6)：7-13.

[15]LEE I，PERSHOUSE R，LEE K.Geospatial Cluster Tessellation Through the Complete Order-k Voronoi Diagrams[C]∥International Conference on Spatial Information Theory.[S.l.]：Springer-Verlag，2007：321-336.

[16]蘭連意，張有會，楊玉平.一般城市Voronoi圖的結晶生成[J].計算機工程與應用，2010(10)：216-219.

[17]涂偉，李清泉，方志祥.基于網絡Voronoi圖的大規模多倉庫物流配送路徑優化[J].測繪學報，2014，43(10):1075-1082.

[18]佘冰，葉信岳，房會會，等.基于局部聚類的網絡Voronoi圖生成方法研究[J].地理科學，2015，35(5):637-643.

[19]ASANO T.Aspect-ratio Voronoi Diagram with Applications[C]∥International Symposium on Voronoi Diagrams in Science and Engineering.[S.l.]：IEEE，2006：32-39.

[20]ASANO T，TAMAKI H.Angular Voronoi Diagram with Applications[C]∥International Symposium on Voronoi Diagrams in Science and Engineering.Banff，Canada：[s.n.]，2006：18-24.

[21]REILLY W J.The Law of Retail Gravitation[M].New York：Knickerbocker Press，1931.

[22]胡鵬，楊傳勇，胡海，等.GIS的基本理論問題—地圖代數的空間觀[J].武漢大學學報(信息科學版)，2002，27(6)：616-621.