空间权重矩阵基础:邻接矩阵、距离矩阵、K最近邻、权重标准化
空间权重矩阵,说白了就是一张“谁跟谁挨着”的表格。Moran I 的整个计算都依赖它——没有权重矩阵,空间自相关就只是个统计玩具。我刚开始学 GIS 时总觉得这步很简单,不就是 1 和 0 嘛?后来踩了好几次坑才明白,权重矩阵的构造直接决定 Moran I 的显著性。
这一节我们就把四种最常用的权重矩阵掰开揉碎讲清楚,包括它们的数学定义、代码实现,以及我实际项目中遇到的坑。
1. 邻接矩阵(Contiguity Matrix)
邻接矩阵是最直观的——两个区域在地理上有没有共享边界?有就是 1,没有就是 0。根据共享边界的方式又分两种:
- Rook 邻接:只考虑共享边(线段),不考虑顶点。像国际象棋的车,只能横竖走。
- Queen 邻接:共享边或共享顶点都算。像象棋的后,八个方向都能走。
我个人习惯在大多数项目中先用 Queen,因为它更保守,不容易遗漏潜在的空间关联。不过要注意:有些不规则多边形(比如海岸线曲折的小岛)在 Queen 定义下可能和很远的多边形通过一个微小的顶点“连接”上,这算是个坑。我曾经在处理欧洲县级数据时就遇到过——两个不相邻的县因为一个极小的测绘误差被 Queen 连上了,Moran I 直接偏掉。后来我增加了最小边界长度过滤,才把这种虚假连接去掉。
用 Python 的 libpysal 库可以轻松生成:
import libpysal
# 读取 shapefile 或 GeoDataFrame
w_rook = libpysal.weights.Rook.from_shapefile('regions.shp')
w_queen = libpysal.weights.Queen.from_shapefile('regions.shp')
# 查看邻居关系
print(w_rook[0]) # 第0个区域的邻居ID列表
2. 距离矩阵(Distance Matrix)
邻接矩阵只适用于面状数据。如果你的数据是点(比如商店位置、气象站),或者区域过大导致少数共享边界无法代表空间关系,那距离矩阵就派上用场了。
计算很简单:两个点(或区域质心)之间的欧氏距离在某个阈值 d 以内,就算邻居。权重可以设为 1(二值)或是距离的倒数(距离越近权重越大)。
关键参数:阈值 d 的选择。设太小会有一堆“孤岛”(没有邻居),设太大又让所有区域都连在一起,失去空间异质性。
我建议先用 k-d 树算所有点对的距离,然后画个距离直方图,取出现频率拐点作为阈值。下面是用 libpysal 生成距离矩阵的示例:
import libpysal
import geopandas as gpd
gdf = gpd.read_file('points.shp')
# 计算质心坐标
centroids = gdf.geometry.centroid
coords = [(p.x, p.y) for p in centroids]
# 基于阈值 100 km 的二值距离矩阵
w_dist = libpysal.weights.DistanceBand.from_array(
coords, threshold=100000, binary=True
)
print(w_dist[0])
另外,如果你的数据范围很大(比如全国各城市),注意坐标系要用投影坐标,避免经纬度直接算距离的误差。我在项目里吃过这个亏——用 WGS84 经纬度算出来的“距离”非常不准,后来改成 UTM 投影才正常。
3. K 最近邻矩阵(K-Nearest Neighbors)
距离矩阵有个硬伤:阈值如果选不好,有的区域邻居很多,有的区域邻居很少甚至为零。K 最近邻解决了这个问题——不管你离多远,我都固定给你找 k 个最近的邻居。
这种做法在样本稀疏的场景特别好用。我记得有一次做城市犯罪率分析,市中心区域非常密集,郊区只有零星几个点。如果用距离矩阵,市中心每个点有一百多个邻居,郊区只有两三个,Moran I 几乎被市中心带偏。换成 K=6 最近邻之后,每个区域都平等对待,结果一下就合理了。
实现也很简单:
import libpysal
coords = [(x1,y1), (x2,y2), ...] # 所有点的坐标列表
w_knn = libpysal.weights.KNN.from_array(coords, k=6)
w_knn.transform = 'R' # 行标准化,稍后讲
print(w_knn[0])
小技巧:K 的取值没有严格规则,我通常从 4 试到 10,对比不同 K 下 Moran I 的变化趋势,选一个稳定区间的值。
4. 权重标准化
不管用什么方法生成了原始权重(0/1 或距离倒数),都需要标准化才能用于 Moran I 计算。最常见的标准化方法有两种:
| 标准化方法 | 公式 | 效果 |
|---|---|---|
| 行标准化(Row-Standardization) | wij' = wij / Σj wij | 每行和为 1 |
| 全局标准化(Global Standardization) | wij' = wij / (n-1) 或 wij' = wij / Σi Σj wij | 所有权重和为 1 或 n |
实践中,行标准化用得最多。原因很简单:每个区域的邻居权重之和固定为 1,Moran I 的分子分母在量级上不会因为某区域邻居爆炸而失衡。我做过的二十多个空间分析项目,90% 都用了行标准化。
避坑指南:有些软件(比如 GeoDa)默认使用行标准化,但 PySAL 的 libpysal 默认不做标准化。我刚开始不知道,直接用原始二值矩阵算 Moran I,结果出来个巨大的正自相关,其实是权重没归一化导致的。一定要显式设置 w.transform = 'R' 来启用行标准化。
另外,如果数据中有孤岛(没有任何邻居的区域),行标准化会出问题——分母为 0。我的处理方式是在生成权重矩阵后手动检查每个区域的邻居数,把孤岛剔除或合并到最近的区域。
# 检查并剔除孤岛
isolated = [i for i in range(len(w)) if w.cardinalities[i] == 0]
if isolated:
print(f"发现 {len(isolated)} 个孤岛,建议处理")
# 可以删除这些样本或使用 K 最近邻强制连接
一张图搞懂四种权重矩阵
下面这张图概括了三种矩阵的构建逻辑以及最终都要经过标准化步骤:
三种输入数据类型分别对应三种矩阵构造方法,最终全部需要经过标准化步骤才能得到可用的标准权重矩阵 W。别小看标准化这步,我见过不少初学者直接用原始矩阵跑 Moran I,结果汇报了一个显著但虚假的正自相关——其实只是权重量级不对。