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

康　順，瞿珊珊

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

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

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

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

1　Voronoi圖定義

1.1　一般Voronoi圖

一般Voronoi圖，又稱普通Voronoi圖，或常規(guī)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)

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

圖2　加權Voronoi圖

(4)

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

2　Voronoi地圖代數(shù)計算

2.1　雷利距離變換

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

圖3　不同商圈的顧客引力

圖4　城區(qū)規(guī)模引力

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

定為

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

(5)

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

2.2　地圖代數(shù)計算

將連續(xù)空間離散化的柵格數(shù)據(jù)是場模型在信息系統(tǒng)中的具體實現(xiàn)。柵格的“位”數(shù)據(jù)蘊含了豐富的關系數(shù)據(jù)，更加適用于網(wǎng)絡分析，在“數(shù)字地球”等大型地理信息工程建設的空間分析中，柵格數(shù)據(jù)展現(xiàn)了明顯的優(yōu)勢[22]。地圖代數(shù)研究的正是度量空間下的抽象點狀要素目標集，即相應空間位置上柵格值的代數(shù)計算。與代數(shù)運算所不同的是，代數(shù)運算只有矩陣的加法和減法，是相應位置上的計算。即地圖代數(shù)遵循了位置一一對應的轉換規(guī)則，是空間計算的重要工具，如圖5、圖6所示的兩個圖幅的地圖代數(shù)乘法計算。以圖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　地圖代數(shù)乘法運算

圖6　最低位置算子￡

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

圖7　地圖代數(shù)的雷利Voronoi圖生成過程

雷利Voronoi圖的地圖代數(shù)生成方法為:

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

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

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

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

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

(5) 算法結束。

3　算例與分析

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

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

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

圖8　試驗區(qū)域

圖9　市區(qū)的一般Voronoi

圖10　市區(qū)的雷利Voronoi圖

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

圖12　各市區(qū)的雷利變換距離

4　結　語

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

參考文獻：

[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的任意發(fā)生元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].中國礦業(yè)大學學報，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]涂偉，李清泉，方志祥.基于網(wǎng)絡Voronoi圖的大規(guī)模多倉庫物流配送路徑優(yōu)化[J].測繪學報，2014，43(10):1075-1082.

[18]佘冰，葉信岳，房會會，等.基于局部聚類的網(wǎng)絡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的基本理論問題—地圖代數(shù)的空間觀[J].武漢大學學報(信息科學版)，2002，27(6)：616-621.