基于核密度估計的點群密度制圖應用研究

盧 敏,楊 柳,王金茵,黃 煌,王結臣

(1.江蘇省地理信息技術重點實驗室，江蘇 南京 210023;2.南京大學 地理信息科學系，江蘇 南京 210023;3.江蘇省基礎地理信息中心，江蘇 南京 210013)

基于核密度估計的點群密度制圖應用研究

盧 敏1,2,楊 柳3,王金茵1,2,黃 煌1,2,王結臣1,2

(1.江蘇省地理信息技術重點實驗室，江蘇 南京 210023;2.南京大學 地理信息科學系，江蘇 南京 210023;3.江蘇省基礎地理信息中心，江蘇 南京 210013)

隨著信息技術的發(fā)展，點狀專題數據激增，需要突破傳統(tǒng)制圖方法的表達局限，利用海量點群數據對專題信息進行快速、直觀、綜合的制圖表達。文中以南京市銀行的ATM網點數據以及全國KTV營業(yè)點分布為制圖數據源，探討了點群數據核密度估計圖以及三維核密度曲面圖的制作方法，并對其中的關鍵技術作了深入的研究。研究表明，基于核密度估計的網點密度制圖方法，消除了邊界分隔對密度計算精度的影響，能更好反映點群數據的空間分布密度、區(qū)域差異、動態(tài)變化等特征，具備更豐富的表現力。

核密度估計；點群密度制圖；最佳帶寬計算；三維核密度曲面圖

點狀分布是專題要素空間分布的最常見狀態(tài)之一，比如銀行、學校、超市、停車場、餐廳等[1]，這些點狀分布的專題要素種類繁多、分布廣泛，點群目標的制圖表達適用范圍廣泛，通過合理運用豐富的專題地圖直觀、有效的可視化表達，可以綜合多樣地展現研究成果、直觀地傳達專題信息，并反映其空間分布的分布密度、區(qū)域差異、動態(tài)變化等特征，使讀圖者一目了然。從整體上看，點群要素最重要的特征是空間分布特征，如何直觀、快速、有效地表達點群目標在空間中的數量、質量、密度和相互關系，是制圖表達的重點[2]。隨著信息技術的發(fā)展，數據獲取能力的增強，許多領域、行業(yè)的點狀專題數據激增，如何突破傳統(tǒng)制圖方法的表達局限，利用海量點群數據對專題信息進行快速、直觀、綜合的制圖表達，成為了制圖領域亟待解決的問題。

對點群要素空間分布特征的表達，應力求展現其分布密度的區(qū)域差異，進而揭示某些潛在規(guī)律，如居民行為習慣、事件發(fā)生熱點、商圈吸引力等[3-4]。利用離散點要素進行制作密度圖時，通常有兩種思路。一種是將整個區(qū)域按指定的邊界(如行政區(qū))進行分割，統(tǒng)計各個區(qū)域內落入的點的數量，計算各區(qū)域的要素密度，從而將點要素轉換為具有指定邊界和不同值的面要素。在該思路中可采用量底法、點值法等方法制圖，但是在這類依邊界分隔的密度圖中，制圖表達的精度與區(qū)域劃分的精細程度密切相關，區(qū)域劃分越小，精度越高，區(qū)域劃分越大，精度也就越低。因此，另一種思路是不進行區(qū)域分割，根據離散的點要素數據內插計算整個區(qū)域的值，從而生成一個連續(xù)的密度表面，通過顏色、符號等在底圖上的疊加來展現專題要素的聚集特點[5]。根據內插原理的不同可分為簡單點密度制圖和核函數密度制圖，可以采用等值線法、量底法、三角形圖表法等方法制圖。點密度制圖與核密度制圖均是通過將點狀要素轉換為連續(xù)曲面，并以面上某點的值代表該點的要素密度，來消除邊界分隔對密度計算精度的影響，較好地展現了點群數據的聚集特點。

本次研究將以南京市所有銀行的ATM(Automatic Teller Machine，自動柜員機)營業(yè)網點數據，以及中國大陸地區(qū)KTV(Karaoke Television，通常指提供卡拉OK影音設備與視唱空間的場所)營業(yè)網點數據為例，探討核密度估計方法在密度制圖中的應用。

1 核密度估計

1.1 核密度估計的基本原理

事件在空間上的出現具有一定的隨意性，但這種隨意性在一定的空間過程作用下會受到影響，從而導致其在不同位置上出現的概率或高或低。如果某一區(qū)域出現更多的事件，則可以認為事件在受到某種空間作用下在此處出現的概率會更高，反之則更低。結合地理學第一定律空間自相似的規(guī)律特征，這種密度(或者概率)上的變化可以用來區(qū)分事件在空間上的分布模式[6-8]。核密度估計(Kernel Density Estimation, KDE)方法就反映了這一思想，它認為區(qū)域內任意位置都有一個可測度的事件密度(也稱強度)，該位置的事件密度可以通過其周圍單位面積區(qū)域內的事件點數量來估計[9-10]。核密度估計可以想象成一個移動的三維函數(內核)，其影響范圍內的事件點到s的距離視為事件點對s影響的權重值，進而計算出點s處的事件密度[11]。

(1)

式中：k( )表示的是核的權重函數，τ>0，稱作帶寬，即以s為原點的曲面在空間上延展的寬度，τ值的選擇會影響到分布密度估計的光滑程度；s-si是密度估值點s和si之間的距離。對于核密度估計的結果，可以采用合適的等值線算法或者以柵格的形式，將研究區(qū)域R以密度變化的連續(xù)表面來表示。

圖1 對點模式的核密度估計

(2)

式中：τ為核函數的帶寬，di為事件點到s的距離。研究表明，不同數學形式和核函數對密度估計的影響很小，對KDE的曲面特征影響較大是帶寬τ的選取。因此,帶寬的選擇比較重要，而核函數可以根據實際需求選擇即可。

1.2 最適帶寬的計算方法

在KDE方法中，帶寬的確定至關重要，在很大程度上會影響到分析的結果。實際的經驗表明，如果想要更光滑的密度表面，τ的值可以設的大一些，然而此時就有可能掩蓋想要探測的熱點區(qū)域；而如果將τ的值設的太小，所得到的密度表面變化就會變得突兀不平，即局部的特征會得到凸顯，而忽略了稍大尺度上的整體關聯性，使得熱點區(qū)域過于分散。因此，在KDE方法的實際應用中，需要采用合適的方法計算得到τ的取值，探索相應帶寬值τ下密度表面的光滑程度，以檢驗該帶寬下的KDE結果是否合理。

一些學者經過研究總結出一些最適值的計算方法，其中最簡單的一種是根據事件點數據的最小外接矩形確定，即用最小外接矩形的較短邊除以30來確定τ值，即

(3)

這種方法缺點是僅考慮了研究區(qū)的形狀特征，而沒要考慮事件點數量的影響。

一種改進的τ值確定方式是，首先確定整個事件點的平均中心，然后計算平均中心到各個事件點的距離，取這些距離的中位數Dm，并計算事件點的標準距離SD(所謂標準距離與標準差類似，是對事件點的一種概括性描述)，則τ值滿足公式[12]

(4)

式中,n為事件點數量。

Bailey和Gatrell給出τ的建議值[13]，其采用事件點數量的0.68倍的-0.2次方與研究區(qū)面積的平方根的乘積確定，即

(5)

式中:n為事件點數量，A為研究區(qū)面積。

上述幾種方法中均沒有考慮事件點之間的空間關系，因此，有學者提出一種基于事件點之間的k階最鄰近距離方法來確定τ值[14-15]，即

(6)

其中,dij表示k階最鄰近距離，即從一個事件點一直到第k個最鄰近點的距離均值。k值決定了密度表面的光滑程度，k值越大帶寬τ也就越大，生成的密度表面也就越光滑。

2 核密度估計的密度制圖

核密度制圖的思路是，給每個點要素定義一個圓形鄰域，在其上方覆蓋一個平滑曲面，依據二次核函數運算使曲面值呈正態(tài)分布，在該點所在位置處曲面的值最高，隨著與點的距離的增大曲面值逐漸減小，到鄰域邊界處減小為零。其中，曲面最高值所在點到圓形鄰域值為零處的距離，就是帶寬。每個輸出像元的密度均為疊加在該柵格像元中心的所有密度曲面值之和。

平面視角下的核密度估計，涉及的主要參數為帶寬。根據已有文獻的研究成果[16-17]，總結出幾種帶寬的計算方法。但是由于這些帶寬的計算方法對于數據集的側重點各不相同，因此分別計算所得到的帶寬值差異較大。在具體研究中，通常的做法是：首先利用這幾種方法分別計算來得到可選的帶寬，然后將其作為參數分別進行核函數估計，最后根據曲面的平滑程度，以及是否能兼顧全局趨勢信息與局部熱點信息，從而最終確立合適的核密度估計帶寬。

本次研究以南京市主城區(qū)的ATM網點數據為例，研究區(qū)包括了南京市的主要城區(qū)范圍(面積約為254.86 km2)，銀行ATM數據來源于網絡電子地圖，研究中主要選擇了南京市主城區(qū)范圍內網點數量大于50的商業(yè)銀行ATM網點，包括工商銀行、建設銀行、農業(yè)銀行等合計749家ATM網點。其他主要參數和帶寬的計算結果如表1所示。

表1 研究區(qū)帶寬計算結果

根據表1的計算結果分別就4種不同的帶寬進行密度估計計算，其結果如圖2所示。由于不同帶寬下計算所得到密度估計值的方法各不相同，為了更好地比較4種帶寬密度估計的影響，在分級渲染的過程中并不能采用統(tǒng)一的分級斷點，這里統(tǒng)一采用“自然斷點”進行分級。

圖2 不同帶寬下的核密度估計結果

由圖2可以看出，隨著帶寬的增大，密度估計表面越來越趨向于光滑，特別是當τ= 4 590 m時，整個密度曲面基本上呈現一種由中心向外圍擴散的圈層結構，雖然這在一定程度上反應了ATM網點在南京主城區(qū)范圍內分布的大體趨勢(由中心城區(qū)向外圍遞減)，但對網點分布的熱點區(qū)域探測并沒有太大的用處。對于其他3種方法而言，采用τ=643 m和τ=1 022 m這兩種帶寬的效果均較為理想，在反映網點分布的整體趨勢的同時，也能夠很好地呈現局部的熱點區(qū)域信息；k階最鄰近距離方法這里取k=3時計算所得的帶寬，即τ=345 m，可以看出其探測出的熱點區(qū)域最為分散，對網點分布熱點區(qū)域的探測而言尺度過小。但是事實上，k階最鄰近距離方法具有很強的靈活性，上述其他幾種方法計算得到的帶寬都能在k階最鄰近距離方法相應的階數下找到對應的值，因此，一種選擇合適帶寬的方式可以是：首先通過默認方法計算其帶寬，然后在k階最鄰近距離方法中找到對應的階數，并以此為基準，增加或減少階數直到找到適合的帶寬，得到較好的視覺效果。

KDE的方法消除了樣方分析中由于樣方的尺寸和形狀等對局部密度的影響，并且具有更好的可視化效果，在探索事件分布熱點、測度局部密度概括指標上具有一定的優(yōu)勢。然而KDE雖然能夠對具有一階密度或平穩(wěn)性的分布模式做出有效且實用的檢驗，但對于二階距離特性的探測還顯得不夠。

3 三維核密度曲面圖

將普通二維專題地圖進行三維顯示，通過對專題信息的立體化表達，能給讀圖者帶來更直接更具沖擊力的視覺感受。同時，三維制圖還有助于減弱專題地圖中對色彩的依賴，使地圖適應傳統(tǒng)印刷工藝對色彩的限制。在立體的現實世界中，具有高低起伏特點的景觀有很多，縱橫交錯的山脈、鱗次櫛比的樓宇，都可以作為三維制圖的效果參照。受人們閱讀習慣和設備制作工藝的限制，目前常用的圖形顯示設備依然采用平面顯示方式，三維效果通常需要依靠物體形狀、顏色、陰影、輔助線等的輔助渲染來體現[18]。

相對于平面地圖而言，三維模擬更形象、直觀，也更符合人們日常觀察事物的習慣。三維表面圖的應用較為廣泛，不僅可以直觀地展現高程、氣溫、人口密度等連續(xù)變量，還可以通過疊加其它要素參數來表達綜合專題信息。在點群要素三維制圖表達中，可根據數據特點和地圖用途的不同，選取適當的制圖表達方法，結合構圖和色彩搭配，可以形成較好的效果。較高的區(qū)域會遮擋后面較矮的區(qū)域，若要較好地掌握數據的整體狀況，需從多個角度進行觀察。對地形的三維模擬多從這兩種視角出發(fā)，或是如身臨其境般的虛擬現實場景，或是用沙盤等立體模型展現地勢起伏。為了更好地展現制圖區(qū)域內的總體情況，地形模擬更多地采用了俯瞰的視角。然而，從地圖制圖的角度來看，在現有的專題地圖制圖表達方法中，三維地圖制圖方法一般用于自然地理要素的表達，如地形的高低起伏、溫度的變化、降雨量地區(qū)差異的展示等[19]，但對于人文地理要素中社會經濟現象的空間分布特征而言，三維的制圖表達手段更加直觀、豐富多樣，但卻是少有人運用該方法來呈現此類人文地理要素的空間差異特征。劉銳、胡偉平等人[20]利用核密度估計模型與基尼系數為切入點，以我國沿海11個省份2000—2012年的數據為例，繪制了海洋經濟發(fā)展的動態(tài)演變趨勢與規(guī)律，但也依舊沒有使用三維制圖進行專題制圖展示。

本次制圖所用的KTV營業(yè)網點數據是通過網絡電子地圖整理得到的。作為一種典型的娛樂休閑行業(yè)，KTV行業(yè)在中國的發(fā)展并不均衡。經統(tǒng)計，在占大陸地區(qū)全部國土面積13.5%的東部地區(qū)12省(或自治區(qū)、直轄市)中，共有KTV營業(yè)網點44 637個，占大陸地區(qū)KTV網點總量的50.9%；在占大陸地區(qū)國土面積29.3%的中部地區(qū)9省(或自治區(qū))中，擁有KTV營業(yè)網點24 883個，占大陸地區(qū)KTV網點總量的28.4%；西部地區(qū)的10個省(或自治區(qū)、直轄市)占據大陸地區(qū)國土面積的57.2%，共有KTV營業(yè)網點18 210個，僅占大陸地區(qū)KTV網點總量的20.7%。利用前文所描述的核密度最佳帶寬確定方法，首先分別計算出4種方法對應的帶寬，然后分別用這4個帶寬作出核密度估計圖，通過視覺的對比，選出其中能兼顧全局趨勢信息與局部熱點信息的密度圖，從而確定制圖的最佳帶寬。通過該最佳帶寬做出平面圖后添加俯瞰視角的三維顯示效果，可以得到如圖3所示的效果圖。

可以看到，將平面的核密度圖拉伸成立體曲面后，密度較高的區(qū)域得到有效拉伸，輔以色彩上的區(qū)分，能夠有效突出了區(qū)域間的差異。經過三維立體化并拉伸顯示后的核密度曲面，高低值之間的差異被適當地夸大，這樣可以有效地凸顯出熱點，如圖3所示的上海、北京、廣州等，均是KTV密度最高的幾個城市。此外，在三維核密度曲面上疊加了注記、圖例、行政邊界等輔助讀圖的地圖要素之外，由于三維透視中存在顯示比例近大遠小的變化，圖中無法添加比例尺，為了彌補這一缺憾，在制圖中還添加了規(guī)則的格網，它不僅使讀圖者對空間位置的把握更準確，并且有指示方向的作用。因此在制圖的過程中，應當根據制圖數據的特點和制圖需要選取合適的角度，全面、清晰、直觀地展現空間和屬性信息，提升地圖的可讀性和可用性。

圖3 中國大陸KTV營業(yè)網點三維密度圖

4 結束語

大多數專題制圖要素在空間上的分布是不均衡的，使用傳統(tǒng)的專題地圖制圖表達手法往往受區(qū)域幾何面積、分區(qū)單元大小等的影響，不能直觀地展現局部統(tǒng)計信息的密度差異。在諸如量底法這類依邊界分隔的密度圖中，制圖表達的精度與區(qū)域劃分的精細程度相關，區(qū)域劃分越小，精度越高，就越能消除分區(qū)間的密度差異，分區(qū)間的密度分布越平滑；反之區(qū)域劃分越大，精度也就越低，分區(qū)間的密度差異越大，密度分布越粗糙。

本文將核密度估計的相關理論運用到空間點群數據的制圖上，分析了不同數據源的情況下KDE函數帶寬的變化對核密度制圖的影響。研究表明，利用核密度估計的方法能夠計算得到較為平滑的密度分布，能有效消除邊界分隔對密度計算造成的影響。而三維曲面圖用立體模型展現其凸起和凹陷，從而使專題信息的展示更加直觀多樣。不僅有效地反映了專題要素的空間形態(tài)、數量特征，而且具有方向性和對比性，表現力出色。這兩種方法是專題制圖方法的擴充與發(fā)展，為點群數據的制圖表達提供了一種新的思路。

[1] 周愛華,張寶秀,董恒年,等. 基于GIS的北京城區(qū)老字號小吃店空間分布特征研究[J].北京聯合大學學報(自然科學版),2015,29(4):29-35.

[2] 曹一冰,付洪貝,潘志超. 地圖設計中的美學研究[J]. 測繪與空間地理信息,2011,34(1): 231-233.

[3] 王勁峰. 地圖的定性和定量分析[J]. 地球信息科學學報,2009,11(2): 169-175.

[4] 禹文豪,艾廷華. 核密度估計法支持下的網絡空間POI點可視化與分析[J]. 測繪學報,2015,44(1):82-90.

[5] 祝國瑞,郭禮珍,尹貢白,等. 地圖設計與編繪[M]. 武漢：武漢大學出版社,2001.

[6] ANDERSON T K. Kernel density estimation and K-means clustering to profile road accident hotspots[J]. Accident Analysis & Prevention,2009,41(3): 359-364.

[7] 閆慶武,卞正富,趙華. 人口密度空間化的一種方法[J]. 地理與地理信息科學,2005，21(5): 45-48.

[8] 張桂銘,朱阿興,楊勝天,等. 基于核密度估計的動物生境適宜度制圖方法[J].生態(tài)學報,2013,33(23):7590-7600.

[9] XIE Z,YAN J. Kernel Density Estimation of traffic accidents in a network space[J]. Computers,Environment and Urban Systems. 2008,32(5): 396-406.

[10] XIE Z,YAN J. Detecting traffic accident clusters with network kernel density estimation and local spatial statistics: an integrated approach[J]. Journal of Transport Geography. 2013,31(0): 64-71.

[11] OKABE A,SATOH T,SUGIHARA K. A kernel density estimation method for networks,its computational method and a GIS-based tool[J]. International Journal of Geographical Information Science. 2009,23(1): 7-32.

[12] ArcGIS10.2在線幫助文檔中核密度估計函數的默認帶寬設定[EB/OL]. http://resources.arcgis.com/en/help/main/10.2/index.html#/How_Kernel_Density _works/009z00000011000000/ . (2014,10).

[13] BAILEY T C,GATRELL A C. Interactive spatial data analysis[M]. Longman Scientific & Technical Essex,1995.

[14] 王遠飛,何洪林. 空間數據分析方法[M]. 北京：科學出版社,2007.

[15] WILLIAMSON D,MCLAFFERTY S,GOLDSMITH V,et al. A better method to smooth crime incident data[J]. ESRI ArcUser Magazine. 1999: 1-5.

[16] SCOTT D W. Multivariate density estimation: theory,practice,and visualization[M]. John Wiley & Sons,2015.

[17] SIMONOFF J S. Smoothing methods in statistics[M]. Springer Science & Business Media,2012.

[18] 徐智勇,艾廷華,危擁軍，等. 三維地圖符號視覺參量研究[J]. 武漢大學學報(信息科學版),2006，31(6): 557-560.

[19] 葉思源,吳樹仁,歐陽永龍. 地質三維可視化建模與其剖面自動制圖應用研究[J]. 地質與勘探,2011,47(3):498-504.

[20] 劉銳,胡偉平,王紅亮,等. 基于核密度估計的廣佛都市區(qū)路網演變分析[J]. 地理科學,2011,31(1):81-86.

[責任編輯：劉文霞]

Applications of pointgroup density cartography based on kernel density estimation

LU Min1,2, YANG Liu3, WANG Jinyin1,2, HUANG Huang1,2, WANG Jiechen1,2

(1. Jiangsu Province Key Laboratory of Geographic Information Science and Technology, Nanjing University, Nanjing 210023, China;2. Dept.of Geographic Information Science, Nanjing University, Nanjing 210023,China;3. Jiangsu Province Surveying & Mapping Engineering Institute, Nanjing 210013,China)

With the development of information technology and the increase of thematic point data, it is necessary to break through the limitation of the traditional cartography and apply mass point group data to express thematic information more quickly, intuitively and comprehensively. This paper, taking ATM network data of banks in Nanjing and KTV network data in the whole country as given data, studies the methods to make kernel density map and 3D kernel density surface map and presents the key technology further. The study shows that the point density cartography based on the kernel density eliminates the influence of edge on the accuracy of calculating density and reflects the density of point distribution, the regional difference, dynamic change and other characteristics, so which has more abundant expression.

kernel density estimation; point group density cartography; optimal bandwidth computation; 3-D kernel density surface map

引用著錄：盧 敏,楊 柳,王金茵，等.基于核密度估計的點群密度制圖應用研究[J].測繪工程，2017，26(4)：70-74，80.

10.19349/j.cnki.issn1006-7949.2017.04.013

2016-02-25；

2016-03-19

國家自然科學基金資助項目(41571377)；國家自然科學基金資助項目(41401450)

盧 敏(1991-)，男，碩士研究生.

P283

A

1006-7949(2017)04-0070-05