2026/9/17 12:19:01
FEATURED · 精选文章 
做空间数据分析这些年被问得最多的一句话就是我这组数据到底算不算有空间聚集。问的人手里往往已经有一张行政区划图、一列指标值Excel 里排得整整齐齐但接下来就卡住了——画个分级设色图看着好像有规律可这个东西能不能写进报告、能不能下结论心里没底。这时候该上场的就是空间自相关而莫兰指数Morans I是其中使用频率最高的那把尺子。它干的事情其实很朴素把相邻的地方数值像不像这件事压缩成一个带显著性检验的数字。数值为正说明高值挨着高值、低值挨着低值也就是聚集为负说明高值旁边净是低值也就是分散或者叫 checkerboard 式的交错接近零则说明空间上基本随机邻居是谁跟你的取值没关系。近几年用得越来越多的双变量空间自相关则是在追问一个更进一步的问题A 指标的高值区是不是恰好落在 B 指标的高值区旁边。这篇内容我打算按自己平时干活的顺序讲——先讲清楚它在问什么、公式每一项在算什么再讲权重矩阵这个最容易翻车的地方然后是 Python 和 R 两套可以直接抄的代码最后是双变量版本和我这些年踩过的坑。适合刚接触空间统计的研究生、做区域分析的规划与市场同学以及需要把空间聚集写进结论但不确定结论站不站得住的人。1. 先搞清楚空间自相关到底在问什么问题1.1 从地理学第一定律说起Tobler 那句被引用到烂的话——任何事物都与其他事物相关但近的东西比远的东西更相关——就是空间自相关全部的理论地基。它描述的是一个叫空间依赖spatial dependence的现象一个位置的观测值不是独立的它受邻居影响。跟它成对出现的另一个概念是空间异质性spatial heterogeneity指的是同一个关系在不同位置上强度不一样比如沿海城市 GDP 和港口吞吐量的关系跟内陆城市就完全是两码事。这两个概念经常被搞混但处理手段完全不同。空间依赖用莫兰指数、Gearys C、Getis-Ord G 这类全局或局部指标来度量空间异质性则要靠地理加权回归GWR、多尺度地理加权回归MGWR这类允许系数随位置变化的模型来处理。我见过不少人拿着明显存在异质性的数据硬套一个全局莫兰指数算出来接近零就宣布没有空间效应这属于典型的工具误用。正确的做法是先看莫兰散点图和局部指标确认是不是正负相消导致的假性零值。为什么非得用专门的指标不能直接算皮尔逊相关系数因为普通相关系数要求样本独立而空间数据的样本天然不独立直接算会高估显著性。更关键的是普通相关系数只关心成对数值同向变化压根不管这对数值在地图上隔多远。莫兰指数把距离通过权重矩阵塞进了公式这才是它区别于普通相关系数的本质。1.2 莫兰指数到底在算什么用一句大白话概括莫兰指数是把每个位置的取值和它邻居的平均取值做一次加权相关然后把所有位置的结果汇总起来。只不过这个相关不是简单的两两相关而是通过空间权重矩阵给不同邻居分配了不同的话语权。它的取值范围理论上落在 -1 到 1 之间-1 表示完全负相关棋盘格0 表示随机1 表示完全正相关同类抱团。这里有个细节必须提醒在权重矩阵行标准化之后莫兰指数的上下界并不是严格锁死在 ±1 的它只保证在极端构造下不超过某个由权重结构决定的值。所以看到 0.85 这种数不用惊讶看到 -1.2 也未必是算错了先去检查权重有没有问题。提示初学者最常犯的错误是把莫兰指数的绝对值大小当成聚集强度直接跨研究比较。不同研究用了不同的权重矩阵、不同的空间尺度0.6 和 0.4 之间没有可比性。要比也是在同一套权重下比不同年份、不同指标。1.3 全局和局部两套指标不要混着用全局莫兰指数Global Morans I给出的是一个整体判断整个研究区有没有空间聚集方向是正还是负显著不显著。它只有一个数回答的是全局格局。局部莫兰指数Local Morans I简称 LISA则是给每一个空间单元算一个值回答的是哪里聚集。它会输出四类结果高-高HH、低-低LL、高-低HL也叫高值被低值包围的异常点、低-高LH。前两类是同向聚集后两类是空间异常值。这两个东西的结论可以不一致而且不一致是正常的。全局指数接近零但局部可能同时存在显著的 HH 和 LL它们互相抵消了。遇到这种情况正确写法是全局未检测到显著空间自相关但局部存在高值与低值聚集区呈两极分化态势而不是简单写一句无空间自相关就完事——那等于把最有价值的信息扔掉了。2. 核心公式拆解与权重矩阵的选型逻辑2.1 公式逐项拆解别背符号要理解动作全局莫兰指数最常见的写法是这样n Σi Σj wij (xi - x̄)(xj - x̄) I --------- * -------------------------------- S0 Σi (xi - x̄)² 其中 n 空间单元个数 xi 第 i 个单元的观测值x̄ 是全部观测值的均值 wij 空间权重矩阵第 i 行第 j 列的元素 S0 所有权重之和即 Σi Σj wij不要被双重求和吓到拆成三步看就清楚了。第一步把每个观测值减去均值得到中心化后的偏离量这一步决定了后面的乘积是正是负。第二步对每一个 i把它自己和所有邻居 j 的中心化偏离量相乘再按权重 wij 加权求和——如果 i 和它的邻居都高于均值乘积为正都低于均值负负得正也是正一高一低乘积为负。第三步把所有 i 的结果加起来除以总离差平方和做归一化再乘上 n/S0 这个缩放系数。所以本质上莫兰指数的正负号是由邻居之间同向偏离均值的程度决定的数值大小则取决于这种同向程度相对于总变异的比例。理解了这一点你就能明白为什么对数据做对数变换、去趋势之后结果会变——因为中心化的基准和离差结构都变了。顺便说一个等价表达在行标准化权重下莫兰指数可以写成I (z W z) / (z z)其中 z 是中心化向量。这个形式在代码里很好用你完全可以用稀疏矩阵自己实现一遍验证库函数的结果。2.2 空间权重矩阵邻接、距离、KNN 怎么选权重矩阵是莫兰指数里唯一由分析者主观决定的部分也是结论最容易受影响的部分。常见的几类权重类型构造逻辑适用场景主要风险Rook 邻接共享边才算邻居规则格网、行政区无缝隙边缘单元邻居少方差大Queen 邻接共享边或点都算邻居最常用的面数据默认选择密集区邻居过多权重被稀释距离带欧氏距离小于阈值即为邻居点数据、连续分布现象阈值选择主观易出现孤立点KNN取最近的 k 个单元点数据、密度不均场景k 值需敏感性分析k4~8 常见反距离wij 1/dij^α强调距离衰减效应需同时设截断距离否则远距离也有影响核函数高斯/四次核加权连续衰减、平滑分析带宽参数敏感我在实际项目里的选择习惯是这样面数据行政区、格网优先用 Queen因为它对边界的细微缝隙更宽容点数据监测站、企业点位优先用 KNN但一定要做 k 的敏感性测试通常 k6 和 k8 各跑一遍看结论稳不稳。如果研究对象存在明显的距离衰减机制比如空气污染的扩散那反距离或核函数更贴合实际但必须同时设一个截断距离否则一个 200 公里外的城市还会以微小权重影响结果逻辑上说不通。注意权重矩阵的选择没有唯一正确解。审稿人或者上级真正在意的是你有没有做敏感性分析而不是你选了哪一种。至少换两套权重跑一遍结论一致才算稳。2.3 行标准化这一步为什么几乎不能省行标准化指的是把权重矩阵的每一行元素除以该行之和使每一行的权重加起来等于 1。做完之后 S0 就等于 n公式简化为I (Σi Σj wij (xi-x̄)(xj-x̄)) / Σi (xi-x̄)²。为什么几乎所有人都这么做核心原因是消除邻居数量差异带来的偏差。一个处于市中心的行政区可能有 12 个邻居边缘的只有 2 个如果不做行标准化邻居多的单元在求和里天然占更大比重结果就被这些单元主导了。行标准化之后每个单元无论有几个邻居它对整体的贡献权重都是 1这样更公平。代价也不是没有。行标准化会让权重矩阵不再对称原本对称的 W 除以不同的行和后变成非对称这在一些理论推导和双变量分析里会带来微妙影响后面讲双变量的时候我会再提。另外反距离权重做完行标准化之后距离衰减的绝对强度信息也被压缩了只剩下相对结构。所以如果你特别关心距离衰减的定量刻画可以考虑保留原始权重另做分析。2.4 显著性检验随机化假设和置换检验莫兰指数本身只是一个统计量光看数值大小没法判断是不是真的有聚集必须做检验。检验的原假设是观测值在研究区内的空间分布是完全随机的。理论上有两种方差计算方式正态假设normality assumption和随机化假设randomization assumption。前者假设数据本身服从正态分布后者只要求数据在空间上可随机置换条件更宽松所以实践中更常用。用随机化假设算出的方差通常略大一些给出的 p 值更保守。不过现在真正主流的做法是置换检验permutation test把 n 个观测值随机打乱后重新分配到 n 个位置上重新计算莫兰指数重复 999 或 9999 次得到一个经验分布。如果真实观测到的 I 落在这个分布的极端尾巴上就认为显著。它的好处是不依赖任何分布假设而且和随机化假设的解析方差高度吻合属于双保险。代码里唯一要注意的是置换次数。999 次是最低配置能得到三位有效数字想要更稳的 p 值比如 0.001 级别的显著性判断得上 9999 次。我自己的习惯是初探用 999正式出结果用 9999。这个计算量在现在的机器上完全可以接受n 在一两千以内也就几十秒。3. 手把手实操从数据到莫兰指数3.1 数据准备与坐标系统这两个坑动手之前先确认三件事。第一空间单元必须没有重叠和缝隙如果是从各种来源拼接的矢量数据边界对不齐会导致本应相邻的两个区县被判为不相邻。第二距离类权重必须用投影坐标系WGS84 经纬度直接算欧氏距离是没有物理意义的误差随纬度变化。国内常用 CGCS2000 对应的投影带或者干脆用 Web Mercator 做相对分析。第三缺失值要提前处理莫兰指数不接受 NaN而删除缺失单元又会改变整张图的拓扑结构所以要么补全要么明确说明删除后重新构建了权重。下面用一份模拟数据演示完整流程。这里刻意用模拟数据而不是真实数据是为了让代码在任何机器上都能跑通你把数据源换掉就能直接用。import numpy as np import pandas as pd import geopandas as gpd import libpysal from libpysal.weights import Queen, KNN, DistanceBand from esda.moran import Moran, Moran_Local # ---------- 第一步读数据并投影 ---------- gdf gpd.read_file(county.shp) gdf gdf.to_crs(epsg3857) # 投影坐标距离类权重必须走这一步 gdf gdf.dropna(subset[green_ratio]).reset_index(dropTrue) # ---------- 第二步构造并检查权重 ---------- w Queen.from_dataframe(gdf, use_indexTrue) w.transform r # 行标准化 print(孤立单元, w.islands) # 必须为空否则要单独处理 print(邻居数描述) print(pd.Series(w.cardinalities).describe()) # 连通性检查 print(连通分量数, len(w.component_labels))跑完这段最先要看的是w.islands。所谓孤立单元island就是在这套邻接下没有任何邻居的单元常见于岛屿、飞地、或者矢量化时边界没接上。只要它不为空莫兰指数的计算要么报错要么给出错误结论必须先解决。# ---------- 第三步全局莫兰指数 ---------- y gdf[green_ratio].values mi Moran(y, w, permutations9999) print(fMorans I {mi.I:.4f}) print(f期望值 E[I] {mi.EI:.4f}) # 理论上约等于 -1/(n-1) print(fz 值 {mi.z_sim:.4f}) print(fp 值(置换) {mi.p_sim:.4f}) print(fp 值(正态) {mi.p_norm:.4f}) # ---------- 第四步局部莫兰指数 ---------- lisa Moran_Local(y, w, permutations9999) gdf[local_I] lisa.Is gdf[p_sim] lisa.p_sim gdf[quadrant] lisa.q # 1HH 2LH 3LL 4HL # 多重检验校正FDR这一步很多人漏掉 from statsmodels.stats.multitest import multipletests reject, p_adj, _, _ multipletests(gdf[p_sim], alpha0.05, methodfdr_bh) gdf[p_adj] p_adj gdf[lisa_sig] np.where(reject, gdf[quadrant], 0) gdf.to_file(lisa_result.shp, encodingutf-8)3.2 R 实现对照spdep 的老牌流程团队里如果有人习惯 R用sfspdep的组合更顺手功能上完全对等。library(sf) library(spdep) nc - st_read(county.shp) nc - st_transform(nc, 3857) nc - nc[!is.na(nc$green_ratio), ] # 构建邻接权重queen TRUE 即 Queen 邻接 nb - poly2nb(nc, queen TRUE, snap 1) cat(孤立单元数, sum(card(nb) 0), \n) lw - nb2listw(nb, style W, zero.policy TRUE) # 全局莫兰指数解析法与置换法对照 moran.test(nc$green_ratio, lw, randomisation TRUE, zero.policy TRUE) set.seed(1234) moran.mc(nc$green_ratio, lw, nsim 9999, zero.policy TRUE) # 局部莫兰指数 locm - localmoran(nc$green_ratio, lw, zero.policy TRUE) nc$Ii - locm[, Ii] nc$pval - locm[, Pr(z ! E(Ii))] nc$padj - p.adjust(nc$pval, method BH)poly2nb里的snap参数值得说一句它按指定容差把几乎相接的边界强行判为相邻单位跟你的坐标系统一致。当权重的连通分量数大于 1也就是出现了互不相邻的孤岛时调大 snap 往往能救回来但调太大又会凭空造出邻居关系我的经验是不要超过最小区县边长的百分之一。3.3 结果怎么读散点图四个象限讲清了故事莫兰散点图的横轴是中心化后的观测值 z纵轴是它的空间滞后也就是邻居观测值的加权平均代码里就是w.sparse.dot(z)。图上每个点代表一个空间单元回归直线的斜率恰好就是莫兰指数。四个象限的含义象限观测值空间滞后含义常见成因右上 HH高高高值被高值包围核心城区、产业集聚区左下 LL低低低值被低值包围连片欠发达区、生态保护区左上 LH低高低值被高值包围大城市边缘的洼地右下 HL高低高值被低值包围资源型孤点、飞地开发区真正要重点看的是LH 和 HL 这两类异常点。它们往往是最有解释价值的部分。比如你在做商业选址分析时发现某个高消费能力的街道被一圈低消费街道包围那大概率存在行政边界效应或者商业配套断层值得实地看看。只报告 HH 和 LL 而忽略异常点等于把分析做了一半。4. 双变量空间自相关热词背后的真实需求4.1 单变量不够用的时候我们在追问什么单变量莫兰指数只能回答某个指标自己有没有抱团。但实际研究里问题往往长这样城市绿地覆盖率高的地方空气质量是不是也更好产业集聚度高的地方人均收入是不是也更高。这类问题里有两个变量而且强调空间上的匹配关系——不是简单的相关而是A 的高值区与 B 的高值区在空间上是否重叠。双变量空间自相关bivariate spatial autocorrelation就是为这个需求量身定做的。它的核心思想是仍然只给一个变量配权重因为权重必须挂在空间单元上但把另一个变量的空间滞后当作对象来比较。4.2 双变量莫兰指数的计算与三个解读陷阱标准化的双变量莫兰指数公式通常写成n Σi Σj wij (zx_i)(zy_j) I_b ----- * ------------------------- S0 sqrt(Σi zx_i²) * sqrt(Σj zy_j²) 其中 zx、zy 分别是两个变量做 z-score 标准化后的值结构上和单变量版本几乎一样差别只在于乘积的两端换成了两个不同的变量。Python 里直接用esda的Moran_BVfrom esda.moran import Moran_BV, Moran_Local_BV x gdf[green_ratio].values # 绿地覆盖率 y gdf[aqi].values # 空气质量指数 bv Moran_BV(x, y, w, permutations9999) print(f双变量 Morans I {bv.I:.4f}, p {bv.p_sim:.4f}) # 双变量局部 lisa_bv Moran_Local_BV(x, y, w, permutations9999) gdf[bv_local_I] lisa_bv.Is gdf[bv_p] lisa_bv.p_sim gdf[bv_q] lisa_bv.q解读的时候有三个陷阱必须绕开。第一双变量莫兰指数是不对称的。I(x, y)和I(y, x)通常不相等。因为行标准化之后的权重矩阵不再对称谁当被滞后的一方会影响结果。这意味着你必须明确说明谁是自变量、谁是因变量不能含混地写两者存在空间关联。第二它只反映空间共现不反映因果。双变量莫兰指数显著只能说明 x 的高值区与 y 的高值区在空间上重叠可能的原因是 x 导致 y、y 导致 x、或者第三个因素同时驱动两者也可能纯粹是共同的空间趋势造成的伪相关。要做因果推断得靠空间计量模型或者准实验设计指数本身给不了。第三象限命名容易混淆。双变量 LISA 的四个象限同样是 HH、LH、LL、HL但这里的H和L指的是两个不同变量的高低。比如第一象限表示x 高、y 的空间滞后也高。在地图上解读时一定要在标注里写清楚哪个轴对应哪个变量我见过不止一次把两组变量的象限标反的图。4.3 一个可以复现的演示流程把双变量分析串起来完整的操作顺序是先各自做单变量莫兰指数确认两个变量本身都有空间结构如果两个都接近零双变量基本不用做了再对两个变量做标准化用同一套权重矩阵计算双变量全局指数然后算双变量局部指数并做 FDR 校正最后把显著单元分类落到地图上。有一个容易被忽略的细节如果两个变量各自都有很强的空间趋势双变量指数会天然偏高因为两者的空间滞后都在系统性地同向偏离。这种情况下应该先对两个变量做趋势面回归把残差拿出来再算。我做过一次测试一个明显存在南北梯度差异的指标不去趋势时双变量指数是 0.41用一阶多项式去趋势之后掉到 0.13结论完全不同。所以只要地图上能看出明显的方向性梯度就老老实实先去趋势。5. 常见问题排查与避坑清单5.1 常见问题速查表下面这张表基本覆盖了我这些年被问过的高频问题遇到报错或结果反常可以先照着扫一遍。现象可能原因排查与解决提示存在孤立单元岛屿、飞地、矢量边界缝隙调大 snap 容差或改用距离/KNN 权重或单独剔除并说明权重连通分量数 1研究区分成互不相连的几块同上先确认是否真的地理隔离I 值异常大1 或 -1权重未行标准化、定义域不匹配检查w.transform确认权重行列顺序与数据一一对应p 值恒为 0.001置换次数太少达到下限把 permutations 提到 9999全局 I 显著但局部一个都不显著多重检验校正后损失功效属正常现象可用 FDR 而非 Bonferroni或放宽到 0.1 并说明换一套权重结论就翻转空间结构对尺度/权重敏感做敏感性分析并如实报告不要只挑好看的那个Moran 散点图点全挤在一团变量严重右偏取对数或 Box-Cox 变换后再算每次运行结果略有不同置换检验的随机性固定随机种子np.random.seed(1234)5.2 那些报告里不会写的实操心得先说权重矩阵的敏感性分析。这是最有性价比的一步。我现在的标准动作是至少跑三套Queen 邻接、KNNk6 或 8、反距离带。如果三套结果在方向和显著性上一致写结论的时候底气就足如果只有一套显著那大概率是权重构造挑出来的假象得老实说明。这套流程多花不了 10 分钟但能挡掉很多质疑。再说多重检验校正这件事。做局部莫兰指数的时候你在同一个数据集上做了 n 次检验。如果 n 是 500用 0.05 的阈值纯随机情况下也会期望出现 25 个显著单元。所以不校正直接出图等于在噪声里挑信号。我通常用 Benjamini-Hochberg 方法控制错误发现率比 Bonferroni 温和得多不会把真正有效果的单元一刀切掉。关于空间尺度的提醒。同一份数据换成市级和区县级莫兰指数可能天差地别这就是经典的可变面元问题MAUP。市级单元大内部差异被平均掉了往往显示出更高的自相关区县级单元小细节保留得多指数通常更低。所以报告结果时一定要写清楚分析单元是什么层级别让人误以为换个尺度还能复现。关于去趋势的实践。判断要不要去趋势有个很土但很有效的办法把变量的分级设色图打印出来贴在墙上退后两米看一眼。如果能一眼看出明显的南北或者东西梯度那就先做一阶多项式趋势面回归用残差的莫兰指数来度量扣除大尺度趋势之后还剩多少局部聚集。这两个结果是两件事不是谁替代谁我在报告里通常两个都放一个讲大格局一个讲局部交互。最后是关于变量预处理的。我见过不少直接把原始计数比如人口总数、GDP 总量扔进模型的例子。总量类变量跟区域面积高度相关本质上度量的是规模而不是强度算出来的莫兰指数反映的其实是大区抱团这种常识。更合理的做法是转成人均、地均或者密度类相对指标。如果一定要用总量那就把面积作为控制变量先回归掉。还有一个细节值得单独拎出来莫兰指数的期望值不是零而是 -1/(n-1)。当 n 很小的时候比如只有 20 个单元期望值是 -0.053这时候即使真实 I 等于 -0.05也不代表存在负相关因为那就是随机状态下的正常表现。所以小样本研究一定要看 z 值而不是光看 I 的绝对值n 小于 30 的时候整个检验的功效都不太够这种情况下结论要写得非常克制。libpysal的权重对象可以序列化保存w.to_file(weights.gal)下次直接读回来省得每次重新构建。不过要留意一点序列化时保存的是单元索引顺序如果中间对数据做过筛选或者排序索引就错位了权重和数值会对不上结果荒诞。我现在固定的习惯是在构建权重之后立刻把w.id_order打印出来跟数据框的索引比对一遍确认无误再往下走。这个动作成本极低但帮我挡过至少三次隐性错误——那种不报错但结果全错的坑才是最要命的。
莫兰指数实战：空间自相关、权重矩阵与双变量分析
做空间数据分析这些年被问得最多的一句话就是我这组数据到底算不算有空间聚集。问的人手里往往已经有一张行政区划图、一列指标值Excel 里排得整整齐齐但接下来就卡住了——画个分级设色图看着好像有规律可这个东西能不能写进报告、能不能下结论心里没底。这时候该上场的就是空间自相关而莫兰指数Morans I是其中使用频率最高的那把尺子。它干的事情其实很朴素把相邻的地方数值像不像这件事压缩成一个带显著性检验的数字。数值为正说明高值挨着高值、低值挨着低值也就是聚集为负说明高值旁边净是低值也就是分散或者叫 checkerboard 式的交错接近零则说明空间上基本随机邻居是谁跟你的取值没关系。近几年用得越来越多的双变量空间自相关则是在追问一个更进一步的问题A 指标的高值区是不是恰好落在 B 指标的高值区旁边。这篇内容我打算按自己平时干活的顺序讲——先讲清楚它在问什么、公式每一项在算什么再讲权重矩阵这个最容易翻车的地方然后是 Python 和 R 两套可以直接抄的代码最后是双变量版本和我这些年踩过的坑。适合刚接触空间统计的研究生、做区域分析的规划与市场同学以及需要把空间聚集写进结论但不确定结论站不站得住的人。1. 先搞清楚空间自相关到底在问什么问题1.1 从地理学第一定律说起Tobler 那句被引用到烂的话——任何事物都与其他事物相关但近的东西比远的东西更相关——就是空间自相关全部的理论地基。它描述的是一个叫空间依赖spatial dependence的现象一个位置的观测值不是独立的它受邻居影响。跟它成对出现的另一个概念是空间异质性spatial heterogeneity指的是同一个关系在不同位置上强度不一样比如沿海城市 GDP 和港口吞吐量的关系跟内陆城市就完全是两码事。这两个概念经常被搞混但处理手段完全不同。空间依赖用莫兰指数、Gearys C、Getis-Ord G 这类全局或局部指标来度量空间异质性则要靠地理加权回归GWR、多尺度地理加权回归MGWR这类允许系数随位置变化的模型来处理。我见过不少人拿着明显存在异质性的数据硬套一个全局莫兰指数算出来接近零就宣布没有空间效应这属于典型的工具误用。正确的做法是先看莫兰散点图和局部指标确认是不是正负相消导致的假性零值。为什么非得用专门的指标不能直接算皮尔逊相关系数因为普通相关系数要求样本独立而空间数据的样本天然不独立直接算会高估显著性。更关键的是普通相关系数只关心成对数值同向变化压根不管这对数值在地图上隔多远。莫兰指数把距离通过权重矩阵塞进了公式这才是它区别于普通相关系数的本质。1.2 莫兰指数到底在算什么用一句大白话概括莫兰指数是把每个位置的取值和它邻居的平均取值做一次加权相关然后把所有位置的结果汇总起来。只不过这个相关不是简单的两两相关而是通过空间权重矩阵给不同邻居分配了不同的话语权。它的取值范围理论上落在 -1 到 1 之间-1 表示完全负相关棋盘格0 表示随机1 表示完全正相关同类抱团。这里有个细节必须提醒在权重矩阵行标准化之后莫兰指数的上下界并不是严格锁死在 ±1 的它只保证在极端构造下不超过某个由权重结构决定的值。所以看到 0.85 这种数不用惊讶看到 -1.2 也未必是算错了先去检查权重有没有问题。提示初学者最常犯的错误是把莫兰指数的绝对值大小当成聚集强度直接跨研究比较。不同研究用了不同的权重矩阵、不同的空间尺度0.6 和 0.4 之间没有可比性。要比也是在同一套权重下比不同年份、不同指标。1.3 全局和局部两套指标不要混着用全局莫兰指数Global Morans I给出的是一个整体判断整个研究区有没有空间聚集方向是正还是负显著不显著。它只有一个数回答的是全局格局。局部莫兰指数Local Morans I简称 LISA则是给每一个空间单元算一个值回答的是哪里聚集。它会输出四类结果高-高HH、低-低LL、高-低HL也叫高值被低值包围的异常点、低-高LH。前两类是同向聚集后两类是空间异常值。这两个东西的结论可以不一致而且不一致是正常的。全局指数接近零但局部可能同时存在显著的 HH 和 LL它们互相抵消了。遇到这种情况正确写法是全局未检测到显著空间自相关但局部存在高值与低值聚集区呈两极分化态势而不是简单写一句无空间自相关就完事——那等于把最有价值的信息扔掉了。2. 核心公式拆解与权重矩阵的选型逻辑2.1 公式逐项拆解别背符号要理解动作全局莫兰指数最常见的写法是这样n Σi Σj wij (xi - x̄)(xj - x̄) I --------- * -------------------------------- S0 Σi (xi - x̄)² 其中 n 空间单元个数 xi 第 i 个单元的观测值x̄ 是全部观测值的均值 wij 空间权重矩阵第 i 行第 j 列的元素 S0 所有权重之和即 Σi Σj wij不要被双重求和吓到拆成三步看就清楚了。第一步把每个观测值减去均值得到中心化后的偏离量这一步决定了后面的乘积是正是负。第二步对每一个 i把它自己和所有邻居 j 的中心化偏离量相乘再按权重 wij 加权求和——如果 i 和它的邻居都高于均值乘积为正都低于均值负负得正也是正一高一低乘积为负。第三步把所有 i 的结果加起来除以总离差平方和做归一化再乘上 n/S0 这个缩放系数。所以本质上莫兰指数的正负号是由邻居之间同向偏离均值的程度决定的数值大小则取决于这种同向程度相对于总变异的比例。理解了这一点你就能明白为什么对数据做对数变换、去趋势之后结果会变——因为中心化的基准和离差结构都变了。顺便说一个等价表达在行标准化权重下莫兰指数可以写成I (z W z) / (z z)其中 z 是中心化向量。这个形式在代码里很好用你完全可以用稀疏矩阵自己实现一遍验证库函数的结果。2.2 空间权重矩阵邻接、距离、KNN 怎么选权重矩阵是莫兰指数里唯一由分析者主观决定的部分也是结论最容易受影响的部分。常见的几类权重类型构造逻辑适用场景主要风险Rook 邻接共享边才算邻居规则格网、行政区无缝隙边缘单元邻居少方差大Queen 邻接共享边或点都算邻居最常用的面数据默认选择密集区邻居过多权重被稀释距离带欧氏距离小于阈值即为邻居点数据、连续分布现象阈值选择主观易出现孤立点KNN取最近的 k 个单元点数据、密度不均场景k 值需敏感性分析k4~8 常见反距离wij 1/dij^α强调距离衰减效应需同时设截断距离否则远距离也有影响核函数高斯/四次核加权连续衰减、平滑分析带宽参数敏感我在实际项目里的选择习惯是这样面数据行政区、格网优先用 Queen因为它对边界的细微缝隙更宽容点数据监测站、企业点位优先用 KNN但一定要做 k 的敏感性测试通常 k6 和 k8 各跑一遍看结论稳不稳。如果研究对象存在明显的距离衰减机制比如空气污染的扩散那反距离或核函数更贴合实际但必须同时设一个截断距离否则一个 200 公里外的城市还会以微小权重影响结果逻辑上说不通。注意权重矩阵的选择没有唯一正确解。审稿人或者上级真正在意的是你有没有做敏感性分析而不是你选了哪一种。至少换两套权重跑一遍结论一致才算稳。2.3 行标准化这一步为什么几乎不能省行标准化指的是把权重矩阵的每一行元素除以该行之和使每一行的权重加起来等于 1。做完之后 S0 就等于 n公式简化为I (Σi Σj wij (xi-x̄)(xj-x̄)) / Σi (xi-x̄)²。为什么几乎所有人都这么做核心原因是消除邻居数量差异带来的偏差。一个处于市中心的行政区可能有 12 个邻居边缘的只有 2 个如果不做行标准化邻居多的单元在求和里天然占更大比重结果就被这些单元主导了。行标准化之后每个单元无论有几个邻居它对整体的贡献权重都是 1这样更公平。代价也不是没有。行标准化会让权重矩阵不再对称原本对称的 W 除以不同的行和后变成非对称这在一些理论推导和双变量分析里会带来微妙影响后面讲双变量的时候我会再提。另外反距离权重做完行标准化之后距离衰减的绝对强度信息也被压缩了只剩下相对结构。所以如果你特别关心距离衰减的定量刻画可以考虑保留原始权重另做分析。2.4 显著性检验随机化假设和置换检验莫兰指数本身只是一个统计量光看数值大小没法判断是不是真的有聚集必须做检验。检验的原假设是观测值在研究区内的空间分布是完全随机的。理论上有两种方差计算方式正态假设normality assumption和随机化假设randomization assumption。前者假设数据本身服从正态分布后者只要求数据在空间上可随机置换条件更宽松所以实践中更常用。用随机化假设算出的方差通常略大一些给出的 p 值更保守。不过现在真正主流的做法是置换检验permutation test把 n 个观测值随机打乱后重新分配到 n 个位置上重新计算莫兰指数重复 999 或 9999 次得到一个经验分布。如果真实观测到的 I 落在这个分布的极端尾巴上就认为显著。它的好处是不依赖任何分布假设而且和随机化假设的解析方差高度吻合属于双保险。代码里唯一要注意的是置换次数。999 次是最低配置能得到三位有效数字想要更稳的 p 值比如 0.001 级别的显著性判断得上 9999 次。我自己的习惯是初探用 999正式出结果用 9999。这个计算量在现在的机器上完全可以接受n 在一两千以内也就几十秒。3. 手把手实操从数据到莫兰指数3.1 数据准备与坐标系统这两个坑动手之前先确认三件事。第一空间单元必须没有重叠和缝隙如果是从各种来源拼接的矢量数据边界对不齐会导致本应相邻的两个区县被判为不相邻。第二距离类权重必须用投影坐标系WGS84 经纬度直接算欧氏距离是没有物理意义的误差随纬度变化。国内常用 CGCS2000 对应的投影带或者干脆用 Web Mercator 做相对分析。第三缺失值要提前处理莫兰指数不接受 NaN而删除缺失单元又会改变整张图的拓扑结构所以要么补全要么明确说明删除后重新构建了权重。下面用一份模拟数据演示完整流程。这里刻意用模拟数据而不是真实数据是为了让代码在任何机器上都能跑通你把数据源换掉就能直接用。import numpy as np import pandas as pd import geopandas as gpd import libpysal from libpysal.weights import Queen, KNN, DistanceBand from esda.moran import Moran, Moran_Local # ---------- 第一步读数据并投影 ---------- gdf gpd.read_file(county.shp) gdf gdf.to_crs(epsg3857) # 投影坐标距离类权重必须走这一步 gdf gdf.dropna(subset[green_ratio]).reset_index(dropTrue) # ---------- 第二步构造并检查权重 ---------- w Queen.from_dataframe(gdf, use_indexTrue) w.transform r # 行标准化 print(孤立单元, w.islands) # 必须为空否则要单独处理 print(邻居数描述) print(pd.Series(w.cardinalities).describe()) # 连通性检查 print(连通分量数, len(w.component_labels))跑完这段最先要看的是w.islands。所谓孤立单元island就是在这套邻接下没有任何邻居的单元常见于岛屿、飞地、或者矢量化时边界没接上。只要它不为空莫兰指数的计算要么报错要么给出错误结论必须先解决。# ---------- 第三步全局莫兰指数 ---------- y gdf[green_ratio].values mi Moran(y, w, permutations9999) print(fMorans I {mi.I:.4f}) print(f期望值 E[I] {mi.EI:.4f}) # 理论上约等于 -1/(n-1) print(fz 值 {mi.z_sim:.4f}) print(fp 值(置换) {mi.p_sim:.4f}) print(fp 值(正态) {mi.p_norm:.4f}) # ---------- 第四步局部莫兰指数 ---------- lisa Moran_Local(y, w, permutations9999) gdf[local_I] lisa.Is gdf[p_sim] lisa.p_sim gdf[quadrant] lisa.q # 1HH 2LH 3LL 4HL # 多重检验校正FDR这一步很多人漏掉 from statsmodels.stats.multitest import multipletests reject, p_adj, _, _ multipletests(gdf[p_sim], alpha0.05, methodfdr_bh) gdf[p_adj] p_adj gdf[lisa_sig] np.where(reject, gdf[quadrant], 0) gdf.to_file(lisa_result.shp, encodingutf-8)3.2 R 实现对照spdep 的老牌流程团队里如果有人习惯 R用sfspdep的组合更顺手功能上完全对等。library(sf) library(spdep) nc - st_read(county.shp) nc - st_transform(nc, 3857) nc - nc[!is.na(nc$green_ratio), ] # 构建邻接权重queen TRUE 即 Queen 邻接 nb - poly2nb(nc, queen TRUE, snap 1) cat(孤立单元数, sum(card(nb) 0), \n) lw - nb2listw(nb, style W, zero.policy TRUE) # 全局莫兰指数解析法与置换法对照 moran.test(nc$green_ratio, lw, randomisation TRUE, zero.policy TRUE) set.seed(1234) moran.mc(nc$green_ratio, lw, nsim 9999, zero.policy TRUE) # 局部莫兰指数 locm - localmoran(nc$green_ratio, lw, zero.policy TRUE) nc$Ii - locm[, Ii] nc$pval - locm[, Pr(z ! E(Ii))] nc$padj - p.adjust(nc$pval, method BH)poly2nb里的snap参数值得说一句它按指定容差把几乎相接的边界强行判为相邻单位跟你的坐标系统一致。当权重的连通分量数大于 1也就是出现了互不相邻的孤岛时调大 snap 往往能救回来但调太大又会凭空造出邻居关系我的经验是不要超过最小区县边长的百分之一。3.3 结果怎么读散点图四个象限讲清了故事莫兰散点图的横轴是中心化后的观测值 z纵轴是它的空间滞后也就是邻居观测值的加权平均代码里就是w.sparse.dot(z)。图上每个点代表一个空间单元回归直线的斜率恰好就是莫兰指数。四个象限的含义象限观测值空间滞后含义常见成因右上 HH高高高值被高值包围核心城区、产业集聚区左下 LL低低低值被低值包围连片欠发达区、生态保护区左上 LH低高低值被高值包围大城市边缘的洼地右下 HL高低高值被低值包围资源型孤点、飞地开发区真正要重点看的是LH 和 HL 这两类异常点。它们往往是最有解释价值的部分。比如你在做商业选址分析时发现某个高消费能力的街道被一圈低消费街道包围那大概率存在行政边界效应或者商业配套断层值得实地看看。只报告 HH 和 LL 而忽略异常点等于把分析做了一半。4. 双变量空间自相关热词背后的真实需求4.1 单变量不够用的时候我们在追问什么单变量莫兰指数只能回答某个指标自己有没有抱团。但实际研究里问题往往长这样城市绿地覆盖率高的地方空气质量是不是也更好产业集聚度高的地方人均收入是不是也更高。这类问题里有两个变量而且强调空间上的匹配关系——不是简单的相关而是A 的高值区与 B 的高值区在空间上是否重叠。双变量空间自相关bivariate spatial autocorrelation就是为这个需求量身定做的。它的核心思想是仍然只给一个变量配权重因为权重必须挂在空间单元上但把另一个变量的空间滞后当作对象来比较。4.2 双变量莫兰指数的计算与三个解读陷阱标准化的双变量莫兰指数公式通常写成n Σi Σj wij (zx_i)(zy_j) I_b ----- * ------------------------- S0 sqrt(Σi zx_i²) * sqrt(Σj zy_j²) 其中 zx、zy 分别是两个变量做 z-score 标准化后的值结构上和单变量版本几乎一样差别只在于乘积的两端换成了两个不同的变量。Python 里直接用esda的Moran_BVfrom esda.moran import Moran_BV, Moran_Local_BV x gdf[green_ratio].values # 绿地覆盖率 y gdf[aqi].values # 空气质量指数 bv Moran_BV(x, y, w, permutations9999) print(f双变量 Morans I {bv.I:.4f}, p {bv.p_sim:.4f}) # 双变量局部 lisa_bv Moran_Local_BV(x, y, w, permutations9999) gdf[bv_local_I] lisa_bv.Is gdf[bv_p] lisa_bv.p_sim gdf[bv_q] lisa_bv.q解读的时候有三个陷阱必须绕开。第一双变量莫兰指数是不对称的。I(x, y)和I(y, x)通常不相等。因为行标准化之后的权重矩阵不再对称谁当被滞后的一方会影响结果。这意味着你必须明确说明谁是自变量、谁是因变量不能含混地写两者存在空间关联。第二它只反映空间共现不反映因果。双变量莫兰指数显著只能说明 x 的高值区与 y 的高值区在空间上重叠可能的原因是 x 导致 y、y 导致 x、或者第三个因素同时驱动两者也可能纯粹是共同的空间趋势造成的伪相关。要做因果推断得靠空间计量模型或者准实验设计指数本身给不了。第三象限命名容易混淆。双变量 LISA 的四个象限同样是 HH、LH、LL、HL但这里的H和L指的是两个不同变量的高低。比如第一象限表示x 高、y 的空间滞后也高。在地图上解读时一定要在标注里写清楚哪个轴对应哪个变量我见过不止一次把两组变量的象限标反的图。4.3 一个可以复现的演示流程把双变量分析串起来完整的操作顺序是先各自做单变量莫兰指数确认两个变量本身都有空间结构如果两个都接近零双变量基本不用做了再对两个变量做标准化用同一套权重矩阵计算双变量全局指数然后算双变量局部指数并做 FDR 校正最后把显著单元分类落到地图上。有一个容易被忽略的细节如果两个变量各自都有很强的空间趋势双变量指数会天然偏高因为两者的空间滞后都在系统性地同向偏离。这种情况下应该先对两个变量做趋势面回归把残差拿出来再算。我做过一次测试一个明显存在南北梯度差异的指标不去趋势时双变量指数是 0.41用一阶多项式去趋势之后掉到 0.13结论完全不同。所以只要地图上能看出明显的方向性梯度就老老实实先去趋势。5. 常见问题排查与避坑清单5.1 常见问题速查表下面这张表基本覆盖了我这些年被问过的高频问题遇到报错或结果反常可以先照着扫一遍。现象可能原因排查与解决提示存在孤立单元岛屿、飞地、矢量边界缝隙调大 snap 容差或改用距离/KNN 权重或单独剔除并说明权重连通分量数 1研究区分成互不相连的几块同上先确认是否真的地理隔离I 值异常大1 或 -1权重未行标准化、定义域不匹配检查w.transform确认权重行列顺序与数据一一对应p 值恒为 0.001置换次数太少达到下限把 permutations 提到 9999全局 I 显著但局部一个都不显著多重检验校正后损失功效属正常现象可用 FDR 而非 Bonferroni或放宽到 0.1 并说明换一套权重结论就翻转空间结构对尺度/权重敏感做敏感性分析并如实报告不要只挑好看的那个Moran 散点图点全挤在一团变量严重右偏取对数或 Box-Cox 变换后再算每次运行结果略有不同置换检验的随机性固定随机种子np.random.seed(1234)5.2 那些报告里不会写的实操心得先说权重矩阵的敏感性分析。这是最有性价比的一步。我现在的标准动作是至少跑三套Queen 邻接、KNNk6 或 8、反距离带。如果三套结果在方向和显著性上一致写结论的时候底气就足如果只有一套显著那大概率是权重构造挑出来的假象得老实说明。这套流程多花不了 10 分钟但能挡掉很多质疑。再说多重检验校正这件事。做局部莫兰指数的时候你在同一个数据集上做了 n 次检验。如果 n 是 500用 0.05 的阈值纯随机情况下也会期望出现 25 个显著单元。所以不校正直接出图等于在噪声里挑信号。我通常用 Benjamini-Hochberg 方法控制错误发现率比 Bonferroni 温和得多不会把真正有效果的单元一刀切掉。关于空间尺度的提醒。同一份数据换成市级和区县级莫兰指数可能天差地别这就是经典的可变面元问题MAUP。市级单元大内部差异被平均掉了往往显示出更高的自相关区县级单元小细节保留得多指数通常更低。所以报告结果时一定要写清楚分析单元是什么层级别让人误以为换个尺度还能复现。关于去趋势的实践。判断要不要去趋势有个很土但很有效的办法把变量的分级设色图打印出来贴在墙上退后两米看一眼。如果能一眼看出明显的南北或者东西梯度那就先做一阶多项式趋势面回归用残差的莫兰指数来度量扣除大尺度趋势之后还剩多少局部聚集。这两个结果是两件事不是谁替代谁我在报告里通常两个都放一个讲大格局一个讲局部交互。最后是关于变量预处理的。我见过不少直接把原始计数比如人口总数、GDP 总量扔进模型的例子。总量类变量跟区域面积高度相关本质上度量的是规模而不是强度算出来的莫兰指数反映的其实是大区抱团这种常识。更合理的做法是转成人均、地均或者密度类相对指标。如果一定要用总量那就把面积作为控制变量先回归掉。还有一个细节值得单独拎出来莫兰指数的期望值不是零而是 -1/(n-1)。当 n 很小的时候比如只有 20 个单元期望值是 -0.053这时候即使真实 I 等于 -0.05也不代表存在负相关因为那就是随机状态下的正常表现。所以小样本研究一定要看 z 值而不是光看 I 的绝对值n 小于 30 的时候整个检验的功效都不太够这种情况下结论要写得非常克制。libpysal的权重对象可以序列化保存w.to_file(weights.gal)下次直接读回来省得每次重新构建。不过要留意一点序列化时保存的是单元索引顺序如果中间对数据做过筛选或者排序索引就错位了权重和数值会对不上结果荒诞。我现在固定的习惯是在构建权重之后立刻把w.id_order打印出来跟数据框的索引比对一遍确认无误再往下走。这个动作成本极低但帮我挡过至少三次隐性错误——那种不报错但结果全错的坑才是最要命的。 
