核密度估计KDE实战：从直方图到空间热力图（ArcGIS/GeoPandas）

2026/7/9 8:05:31

核密度估计KDE实战：从直方图到空间热力图（ArcGIS/GeoPandas）

当我们需要分析城市商业网点的空间分布特征时，传统点状地图往往难以直观展示密度差异。核密度估计（Kernel Density Estimation，KDE）作为空间分析利器，能够将离散的POI点数据转化为连续的热力图表面，揭示隐藏的空间分布模式。本文将手把手带您完成从基础原理到地理空间应用的完整流程。

1. 核密度估计的核心思想

想象一下城市中快餐店的分布——有些区域密集如林，有些则零星散落。直方图虽然能统计各区间的店铺数量，但存在两个明显缺陷：

边界不连续（相邻区间可能显示截然不同的值）
对区间划分敏感（不同的区间宽度会导致完全不同的图形）

核密度估计通过"滑动的概率云"解决了这些问题。其数学本质是：每个数据点都贡献一个核函数（如高斯曲线），最终密度估计是所有核函数的叠加。二维空间中的计算公式为：

$$ \hat{f}(x,y) = \frac{1}{nh^2}\sum_{i=1}^n K\left(\frac{d_i}{h}\right) $$

其中：

$K$ 为核函数（常用Epanechnikov或高斯核）
$h$ 是带宽（bandwidth），决定平滑程度
$d_i$ 是点$(x,y)$到第$i$个样本点的距离

重要提示：带宽选择直接影响结果——过小会导致噪声突出，过大会掩盖真实特征。Scott规则建议：$h = n^{-1/6} \times \sqrt{\sigma_x^2 + \sigma_y^2}/2$

2. 地理空间KDE的特殊考量

与传统一维KDE不同，地理空间分析需要额外注意：

投影选择：

在经纬度坐标（WGS84）下直接计算会导致距离失真
建议先转换为等面积投影（如Albers或UTM）

边缘校正：

靠近研究区域边界的点会损失部分核函数面积
可通过权重补偿或边界缓冲解决

空间自相关：

地理数据常呈现聚集性（如商业中心）
可能需要考虑各向异性核函数

常用核函数对比：

核类型	数学形式	优点	缺点
高斯核	$e^{-d^2/2}$	无限可微	计算量大
Epanechnikov	$(1-d^2)I(d≤1)$	效率最优	不光滑
四次核	$(1-d^2)^2I(d≤1)$	平衡性好	需自定义

3. 实战演练：河南省KFC分布分析

3.1 数据准备

使用Python获取POI数据（示例代码）：

import requests from geopy.geocoders import Nominatim def get_poi(region, keywords, api_key): url = f"https://restapi.amap.com/v3/place/text?key={api_key}&keywords={keywords}&city={region}" response = requests.get(url).json() return [(p['location'], p['name']) for p in response['pois']] # 示例：获取郑州KFC数据 kfc_zhengzhou = get_poi("郑州", "肯德基", "your_amap_key")

3.2 GeoPandas实现

完整处理流程：

import geopandas as gpd from shapely.geometry import Point import matplotlib.pyplot as plt from sklearn.neighbors import KernelDensity import numpy as np # 创建GeoDataFrame geometry = [Point(float(lon), float(lat)) for (lon,lat), _ in kfc_data] gdf = gpd.GeoDataFrame(geometry=geometry, crs="EPSG:4326") # 转换为等面积投影 gdf = gdf.to_crs("EPSG:3395") # World Mercator # 准备网格 xmin, ymin, xmax, ymax = gdf.total_bounds grid_size = 1000 # 1km网格 xi = np.linspace(xmin, xmax, int((xmax-xmin)/grid_size)) yi = np.linspace(ymin, ymax, int((ymax-ymin)/grid_size)) xx, yy = np.meshgrid(xi, yi) # 计算KDE coords = np.vstack([gdf.geometry.x.values, gdf.geometry.y.values]).T kde = KernelDensity(bandwidth=5000, kernel='epanechnikov') # 5km带宽 kde.fit(coords) zi = np.exp(kde.score_samples(np.vstack([xx.ravel(), yy.ravel()]).T)) zi = zi.reshape(xx.shape) # 可视化 fig, ax = plt.subplots(figsize=(12,10)) gdf.plot(ax=ax, color='red', markersize=5) contour = ax.contourf(xx, yy, zi, levels=20, cmap='RdYlGn_r') plt.colorbar(contour, label='核密度值')

3.3 ArcGIS Pro方案

对于企业用户，ArcGIS提供更完整的工具链：

使用Point Density工具生成原始密度栅格
通过Kernel Density工具选择不同核函数
用Reclassify工具进行密度分级
结合Spatial Analyst扩展进行热点分析

关键参数设置建议：

搜索半径（Search Radius）：根据城市规模设置（中小城市3-5km，大城市5-10km）
输出像元大小（Output Cell Size）：保持与搜索半径1:10的比例
面积单位（Area Units）：推荐使用"平方千米"

4. 进阶技巧与陷阱规避

4.1 动态带宽选择

固定带宽可能不适用于不均匀分布。可采用自适应带宽：

from sklearn.cluster import DBSCAN # 识别密集区域 coords = gdf.geometry.apply(lambda p: [p.x, p.y]).tolist() clustering = DBSCAN(eps=3000, min_samples=5).fit(coords) # 为不同集群设置不同带宽 bandwidths = [3000 if label == -1 else 5000 for label in clustering.labels_]

4.2 时空核密度

结合时间维度的STKDE（Spatio-Temporal KDE）： $$ \hat{f}(x,y,t) = \frac{1}{nh^2h_t}\sum K_s\left(\frac{d_i}{h}\right)K_t\left(\frac{\Delta t}{h_t}\right) $$ 其中$h_t$是时间带宽，$K_t$是时间核函数。