改進(jìn)核密度估計(jì)的空間點(diǎn)密度算法

饒加旺，馬榮華

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

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

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

1 核密度算法

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

式中，n為點(diǎn)的個(gè)數(shù)，λ(u)可通過(guò)核密度算法進(jìn)行無(wú)參化計(jì)算得到[1]。

1.1 劃分格網(wǎng)

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

1.2 獲取密度貢獻(xiàn)值

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

1.3 輸出核密度

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

e為邊緣校正參數(shù)，可表示為λ(v)經(jīng)過(guò)核密度函數(shù)平滑后的估計(jì)值，可視化輸出的效果是將二維的離散點(diǎn)變成三維連續(xù)的光滑平面，通常具有偏差，因?yàn)榍罢呦它c(diǎn)位的離散特征、真實(shí)位置等具體信息[1]。

2 空間點(diǎn)密度算法

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

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

圖1 算法流程圖

2.1 離散點(diǎn)的分箱與降維

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

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

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

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

偽代碼如下：

2.3 維度的再膨脹

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

3 實(shí)證研究

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

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

3.1 可視化運(yùn)算實(shí)驗(yàn)結(jié)果

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

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

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

3.2.2 算法效率

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

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

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

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

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

表2 算法運(yùn)行時(shí)間對(duì)比情況s

4 總結(jié)與展望

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

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