改進核密度估計的空間點密度算法

饒加旺，馬榮華

1.江蘇省測繪工程院 空間信息技術(shù)研究中心，南京210013 2.中國科學(xué)院 南京地理與湖泊研究所，南京210008

點密度是一定范圍內(nèi)點數(shù)量的統(tǒng)計值，是地理空間分析的重要任務(wù)[1-2]，而核密度分析是點密度分析常用的重要方法，它是根據(jù)核密度估計函數(shù)（Kernel Density Estimator）將平面的二維離散點生成連續(xù)的三維表面，計算事件點在設(shè)定的周圍鄰近空間的密度的過程，直觀地反映點群的聚集或離散分布特征[1]。相比較于其他空間分析方法，核密度分析需要參數(shù)較少，受主觀因素影響較小，因而成為了地理空間分析中應(yīng)用最廣泛的方法之一[2]，被廣泛應(yīng)用至地物空間分布[3-4]、空間區(qū)域格局分析[5]、疫情監(jiān)測與分析[6]、地質(zhì)災(zāi)害與自然資源環(huán)境監(jiān)測[7-8]、路徑優(yōu)化與分析[9]等諸多領(lǐng)域，從空間上獲取事件的宏觀分布特征。研究者相繼開發(fā)出適應(yīng)R、Python、ArcGIS等多種編程語言和軟件環(huán)境下的算法。Kern-Smooth是R語言環(huán)境下用于計算核密度的功能包，按照核函數(shù)的不同可分為bkde、bkde2D、bkfe等函數(shù)，其中kde2d是R語言環(huán)境下應(yīng)用最為廣泛的密度函數(shù)[10-11]，該函數(shù)主要使用高斯核密度估計實現(xiàn)函數(shù)的平滑，以三維連續(xù)表面的形式輸出點的概率密度[12]。KernelDensity為Python語言環(huán)境下計算核密度常用的功能包，實現(xiàn)方式、功能與R語言中的KernSmooth類似；ArcGIS軟件在Spatial Analysis工具箱中集成了核密度分析功能，用于衡量單位面積內(nèi)點要素的量值[13]，在可視化輸出方面常借助核密度函數(shù)（或平滑函數(shù)），將離散點在空間上擬合為連續(xù)的光滑錐形表面，再以不同的顏色分級表示密集程度[14]。

相關(guān)研究發(fā)現(xiàn)，核密度分析對于描述連續(xù)現(xiàn)象的應(yīng)用場景效果較好[1]，而對于離散點，雖能表示空間格局、事件的分布特征及趨勢變化，但密度空間分布模式上存在的“距離衰減效應(yīng)”[15]，在利用ArcGIS軟件或kde2d等算法實現(xiàn)核密度分析時常存在如下問題：（1）分析的結(jié)果通常會把本身并不活躍的區(qū)域納入計算范圍，不僅改變了點原有的空間位置，也改變點的離散性質(zhì)，存在過度擬合的情況；（2）運算時需預(yù)先設(shè)置計算空間范圍[10，13]，常導(dǎo)致額外的資源開銷，當(dāng)數(shù)據(jù)量較大時，運算耗時較長、效率較低；（3）可視化輸出時，不同的顏色分級表示的密度分布難以反映離散點真實的密度情況；（4）主觀影響較強，設(shè)置帶寬、搜索半徑與輸出單元格等參數(shù)通常依靠經(jīng)驗，較大的搜索半徑與輸出單元格常導(dǎo)致難以準確地反映事件分布情況，整體空間形狀變形較大，而較小的搜索半徑與輸出單元格又會對核密度分析帶來巨大的計算量與計算資源的消耗。在不改變離散點位置的基礎(chǔ)上獲取點的真實密度值在疫情監(jiān)測[16]、軍事領(lǐng)域、災(zāi)害監(jiān)測與統(tǒng)計[17-18]、地圖制圖[19]、社會治安與應(yīng)急管理[20]等領(lǐng)域中具有重要應(yīng)用。為此本文在核密度分析的基礎(chǔ)上提出了空間點密度算法，利用幾何投影與哈希數(shù)據(jù)結(jié)構(gòu)，在不改變離散點空間位置以及離散特性的基礎(chǔ)上，降低主觀性對結(jié)果的影響，快速度量點真實數(shù)量的方法，以提高點密度的真實分布與高效的可視化輸出。

1 核密度算法

對于非均勻分布的離散點而言，假設(shè)用λ(u)表示空間位置u的密度，將任意空間范圍B劃分為無限多個子區(qū)域ai，范圍B內(nèi)點密度可表示為λ(u)ai，因此B范圍內(nèi)點密度E[B]可表示為所有子區(qū)域點密度之和，則：

式中，n為點的個數(shù)，λ(u)可通過核密度算法進行無參化計算得到[1]。

1.1 劃分格網(wǎng)

將研究區(qū)域劃分為多個單元格網(wǎng)，定義搜索窗口W，逐一遍歷單元格網(wǎng)，對于特定位置v，當(dāng)其在W內(nèi)時f(u)=κ(v-u)（f為空間位置u的函數(shù)，κ為核函數(shù)），反之，則f(u)=0。

1.2 獲取密度貢獻值

獲取W內(nèi)各個格網(wǎng)對整體的密度貢獻值，其值為W內(nèi)所有點xi數(shù)量之和，用坎貝爾方程來表示[1]：

1.3 輸出核密度

計算并輸出每個單元格網(wǎng)的密度值，該值為W內(nèi)所有點對格網(wǎng)密度貢獻值的和，根據(jù)公式（2），則：

e為邊緣校正參數(shù)，可表示為λ(v)經(jīng)過核密度函數(shù)平滑后的估計值，可視化輸出的效果是將二維的離散點變成三維連續(xù)的光滑平面，通常具有偏差，因為前者消除了點位的離散特征、真實位置等具體信息[1]。

2 空間點密度算法

本文提出的空間點密度算法是在保持點的離散屬性與初始空間位置的基礎(chǔ)上，統(tǒng)計并計算每個搜索鄰域內(nèi)點的單位真實數(shù)量，快速輸出離散點的位置與點密度值。假設(shè)圓心O代表了網(wǎng)格點最近的一個事件，研究區(qū)域內(nèi)有n個點，其緯度、經(jīng)度分別由lati、loni表示，g代表格網(wǎng)的大小，r表示用于搜索的鄰域半徑，tlati和tloni為最鄰近格網(wǎng)點的緯度和經(jīng)度；為了使構(gòu)建點密度算法具有可解析性以及保證算法運行的高效性，定義兩個哈希表bin_density_hash（）與density_hash（），其key值分別為鄰近格網(wǎng)點與鄰域內(nèi)格網(wǎng)點的經(jīng)緯度坐標值。

算法分為離散點的分箱與降維操作、計算格網(wǎng)內(nèi)點的密度、維度的再膨脹三個步驟。算法的流程圖如圖1所示。

圖1 算法流程圖

2.1 離散點的分箱與降維

執(zhí)行時初始化參數(shù)r與g，將研究區(qū)域劃分為等距的格網(wǎng)，使鄰域中心點與格網(wǎng)相疊加，對鄰域內(nèi)的離散點執(zhí)行分箱操作，分箱的規(guī)則是通過近似的方式獲取每個點在其最近格網(wǎng)點的坐標。當(dāng)執(zhí)行完成時，數(shù)據(jù)的維度由n降至m，bin_density_hash（）的地址值為裝箱操作后點的數(shù)量值，該值為改變點的位置所得到的，為相對值，作為后續(xù)運算重要的數(shù)據(jù)基礎(chǔ)，偽代碼如下：

2.2 計算格網(wǎng)點內(nèi)的密度

該步驟是本文算法的核心，示意圖如圖2所示。首先遍歷m個格網(wǎng)點，獲取鄰域內(nèi)格網(wǎng)點F的坐標(lat1,lon1)；其次沿著中央經(jīng)線方向自下而上迭代運行，依次從F點遍歷至G點，通過上下移動的距離l，以及鄰域半徑r計算容忍距t；然后以g為步長遍歷經(jīng)度方向的點E1至E4，并將其經(jīng)緯度坐標存儲至density_hash（）表中，統(tǒng)計被記錄次數(shù)，依此類推直至遍歷至點G完成鄰域內(nèi)所有格網(wǎng)點的循環(huán)。容忍距t是識別鄰域內(nèi)點經(jīng)度的重要依據(jù)，確保了只有在鄰域范圍內(nèi)的密度值是增加，避免了算法額外的循環(huán)，是本文算法高效運行的關(guān)鍵。當(dāng)遍歷結(jié)束后，density_hash(latt,lont)的key值為經(jīng)過分箱操作后點的經(jīng)緯度坐標值，地址值為離散點真實的密度值，此時數(shù)據(jù)的維度仍為m。

圖2 計算格網(wǎng)內(nèi)點的密度示意圖

偽代碼如下：

2.3 維度的再膨脹

遍歷n個點，依次從density_hash(latt,lont)的地址值中讀取點的密度值；將離散點初始的經(jīng)緯度坐標與實際的點密度值相關(guān)聯(lián)，若未完成遍歷，重新執(zhí)行步驟2.1～2.2，直至完成整個研究區(qū)域的遍歷，最終輸出結(jié)果為n個離散點的初始坐標與對應(yīng)的密度值。偽代碼如下：

3 實證研究

美國地質(zhì)調(diào)查局（United States Geological Survey，USGS）的國家空間數(shù)據(jù)基礎(chǔ)設(shè)施節(jié)點（National Spatial Data Infrastructure Node）提供了美國50個州的水資源數(shù)據(jù)信息[21]，共計約190多萬條，包括地表水、地下水等類型的站點編號、站點名稱、經(jīng)緯度坐標等字段，該數(shù)據(jù)作為開放的數(shù)據(jù)集，具有數(shù)據(jù)格式化程度高、區(qū)域跨度和數(shù)據(jù)量均較大等特點。本文基于R語言編寫網(wǎng)絡(luò)爬蟲程序，抓取了截止至2014年12月31日的美國大陸（不包括阿拉斯加、夏威夷、關(guān)島等地區(qū)）所有的地下水資源的數(shù)據(jù)信息，共由1 233 186個離散點組成的數(shù)據(jù)集，作為本文算法的驗證數(shù)據(jù)，如圖3所示。

圖3 驗證數(shù)據(jù)的情況

3.1 可視化運算實驗結(jié)果

在R語言環(huán)境下，讀取美國大陸地下水資源數(shù)據(jù)集。兼顧可視化效果與運算次數(shù)，將格網(wǎng)大小g與鄰域半徑r分別設(shè)置為1、5，借助ggplot2功能包[22]，將點密度執(zhí)行的結(jié)果可視化，返回結(jié)果包括離散點的初始坐標、密度值，如圖4所示。

圖4 空間點密度可視化結(jié)果

3.2 與核密度算法的比較

3.2.1 可視化效果

在ArcGIS 10.4.1軟件環(huán)境下，執(zhí)行Kernel Density操作，分別將搜索半徑與輸出單元格的大小設(shè)置為0.1、10，得出如圖5（a）所示的結(jié)果。由圖可知，可視化輸出的結(jié)果為連續(xù)的，圖4中空白區(qū)域也被擬合成連續(xù)的區(qū)域，整體形狀發(fā)生了改變，圖中點的位置發(fā)生了變化，存在過度擬合的情況，當(dāng)搜索半徑與輸出單元格的值均較大時，整體的形狀改變越劇烈；而過于密集的網(wǎng)格對核密度分析帶來巨大的計算量和計算機軟硬件資源的消耗，此外，由于整體形狀發(fā)生了改變，此時的密度值為空間加權(quán)平均值，導(dǎo)致密度值的真實性與可識別性均較低；在R語言環(huán)境中執(zhí)行kde2d算法[14]，分別設(shè)置帶寬與輸出單元個數(shù)為0.1、3 000，輸出的核密度結(jié)果為碎片化的面狀形式如圖5（b）所示，存在離散點位置發(fā)生變化而被過度擬合的情況，密度的可識別性較差，同時存在點缺失的問題。當(dāng)設(shè)置較大的帶寬與輸出單元時，存在更多區(qū)域被過度擬合的情況。本文算法最大程度降低了主觀性對于結(jié)果的影響，通過設(shè)置不同的格網(wǎng)與鄰域半徑，如圖4與圖5（c）（格網(wǎng)大小與鄰域半徑分別設(shè)置為10、100）所示，改變的僅是點的相對密度值，而不會改變離散點整體的形狀，輸出的結(jié)果仍為設(shè)置的格網(wǎng)與鄰域半徑下真實的點密度值。因此本文算法的可視化效果與點密度的可識別性均優(yōu)于ArcGIS 10.4.1軟件與kde2d算法。

圖5 與ArcGIS核密度、kde2d的可視化效果比較

3.2.2 算法效率

評價算法的效率，通常通過統(tǒng)計算法中基本操作重復(fù)執(zhí)行的次數(shù)（時間復(fù)雜度）來衡量，使用O(n)表示[23]。kde2d算法的時間復(fù)雜度為O(nm)[11]，本文算法的時間復(fù)雜度為O(n(2r/g)2)。為驗證不同樣本數(shù)量下兩種算法的時間復(fù)雜度對比情況，依次從數(shù)據(jù)集中隨機取5 000、10 000、50 000、100 000、500 000、1 000 000個點共六組樣本，樣本分布如圖6（a）～（f）所示，依次通過本文點密度算法，得到的結(jié)果如圖7（a）～（f）所示，兩種算法的時間復(fù)雜度的比值情況如表1所示。

本文算法與kde2d算法的時間復(fù)雜度均與樣本數(shù)量n有關(guān)，前者還與設(shè)置的格網(wǎng)大小g與鄰域半徑r相關(guān)，后者與離散點分箱操作后格網(wǎng)點的數(shù)量m相關(guān)。將兩者的樣本數(shù)量設(shè)置為相同，如表1所示，當(dāng)數(shù)據(jù)量從小到大時，由于m?(2r/g)2本文的空間點密度算法的時間復(fù)雜度均低于kde2d算法；當(dāng)數(shù)據(jù)量逐漸增大時，本文的算法效率優(yōu)勢較明顯，表明本文算法更適合大數(shù)據(jù)環(huán)境下計算離散點的密度值。因Environmental Systems Research Institute（ESRI）官方尚未公布ArcGIS軟件Kernel Density算法的源碼[13]，因此以運行時間為指標將本文空間點密度算法與ArcGIS軟件Kernel Density算法進行對比，計算機運行環(huán)境為CPU：Intel I7-6700HQ 2.60 GHz，內(nèi)存：16 GB，操作系統(tǒng)：Windows 7專業(yè)版64位，軟件系統(tǒng)：R 3.5.0、ArcGIS 10.4.1版本，執(zhí)行時將Kernel Density搜索半徑與輸出單元格設(shè)置為1、5，將本文算法的格網(wǎng)大小與鄰域半徑設(shè)置為1、5，兩種算法各執(zhí)行20次，取平均所耗費的時間如表2所示。

總體而言本文算法的運行時間比ArcGIS的Kernel Density較長，尤其是當(dāng)數(shù)據(jù)量較大時。究其原因，首先ArcGIS軟件的Kernel Density函數(shù)封裝與模塊化應(yīng)用較好，可較大程度地利用計算機資源，程序讀寫與運算的效率較高；其次在算法機理上，Kernel Density是對二維離散點表面進行內(nèi)插，計算格網(wǎng)對整個區(qū)域密度的貢獻值[13]，因此遍歷的次數(shù)相對較少，而本文算法為保持離散點的初始位置與獲取真實的點密度，在執(zhí)行分箱操作、計算格網(wǎng)內(nèi)點的密度、維度的再膨脹過程時遍歷次數(shù)較多，導(dǎo)致耗時較長。

圖7 不同樣本數(shù)量下的空間點密度計算結(jié)果

表1 不同樣本數(shù)量下算法效率對比情況

表2 算法運行時間對比情況s

4 總結(jié)與展望

針對目前核密度分析計算點密度存在著改變離散點位置與屬性、密度值可識別性較低、運算效率較低以及主觀因素對輸出結(jié)果影響較大等局限，提出了空間點密度算法，該算法在保持點的離散屬性與初始空間位置的基礎(chǔ)上，通過設(shè)定分箱規(guī)則，獲取離散點最近的格網(wǎng)點坐標；依次遍歷、統(tǒng)計并計算搜索鄰域內(nèi)點的數(shù)量，最終輸出離散點初始坐標與點密度值，算法共分為離散點的分箱與降維、計算格網(wǎng)內(nèi)點密度、維度的再膨脹三個步驟。

引入USGS網(wǎng)站發(fā)布的美國大陸123多萬條地下水資源信息作為驗證數(shù)據(jù)，通過驗證與分析，該算法最大程度地降低了主觀因素對于計算結(jié)果的影響，不改變離散點的屬性與整體形狀，在可視化效果與點密度識別性上均優(yōu)于ArcGIS 10.4.1版本的核密度法與kde2d算法，運算效率優(yōu)于常用的kde2d算法。由此可視化效果與算法時間復(fù)雜度驗證均證明了本文空間點密度算法的有效性。但實驗結(jié)果同時表明，當(dāng)數(shù)據(jù)量為百萬級別以上時，將格網(wǎng)點與鄰域設(shè)置較小時，存在消耗計算機較大內(nèi)存、運算時間較長的情況。為此，后期工作重點將進一步優(yōu)化算法，減小內(nèi)存消耗、縮短運算時間。