先抛一个做区域数据分析的人迟早会撞上的现象每个区县的数据明明是独立采集的画到地图上却总是一片一片地扎堆——高值区域挨着高值区域低值区域挨着低值区域。你把这件事讲给只学过基础统计的同事听对方十有八九会反问一句这不就是相关性吗跑个皮尔逊系数不就完了问题恰恰出在这儿。皮尔逊系数假设样本之间彼此独立可地理数据天生就不独立。空间自相关Spatial Autocorrelation就是用来度量位置这个维度到底有没有在数据里起作用的一套方法而莫兰指数Morans I是其中用得最多、也最经得起推敲的一个统计量。它回答的是一个很朴素的问题一个区域的值和它邻居的值到底是同向变化、反向变化还是毫无关系。这篇内容适合做城市规划、房地产评估、公共卫生、林业遥感、经济地理、门店选址、环境监测的人看。只需要一点点基础统计概念加一点 Python 或 GeoDa 的操作能力就能跟上。我会从它到底在算什么讲起把权重矩阵这个最容易被忽略的环节掰开揉碎然后完整跑一遍代码再说说实际项目里那些让莫兰指数翻车的坑。1. 从相邻区域总是相似说起莫兰指数要解决的到底是什么问题1.1 独立同分布假设在空间数据上为什么会崩任何一门统计课开头都会反复强调一件事样本必须独立。问卷里第 3 个人和第 4 个人的答案互不影响实验组和对照组的观测互不干扰这是所有推断的前提。但地理数据从采集那一刻起就在破坏这个前提。相邻的两个小区共享同一条主干道、同一所学校、同一片商圈房价怎么可能互不影响相邻的两个县共享一片流域、同一条气流通道PM2.5 浓度怎么可能是独立的。这种近的比远的更像的现象在地理学里有个很直白的表述任何事物都与其他事物相关但近的事物比远的事物关联更强。这个事实带来的后果非常实际。如果你把一份明显存在空间聚集的数据丢进普通最小二乘回归模型给出的标准误会系统性偏小p 值会虚低你会得到一堆显著的系数其中相当一部分是假的。这不是模型写错了而是前提被违反了。所以在做任何空间数据的建模之前先跑一次空间自相关检验应该是一个肌肉记忆级别的动作。1.2 把位置关系写成数学莫兰指数的构造思路莫兰指数的聪明之处在于它把相邻这件事抽象成了一个矩阵然后把传统的相关系数结构套了上去。定义一个空间权重矩阵 W里面每个元素 w_ij 表示区域 i 和区域 j 之间的空间关系强度。相邻就是 1不相邻就是 0或者用距离衰减函数给出一个连续值。有了这个矩阵莫兰指数的公式长这样n Σi Σj w_ij (x_i - x̄)(x_j - x̄) I -------- * ---------------------------------- S_0 Σi (x_i - x̄)^2 其中 S_0 Σi Σj w_ijn 是区域总数拆开看其实很好理解。分子里那一坨(x_i - x̄)(x_j - x̄)是区域 i 和它的邻居 j 各自偏离均值的乘积。如果 i 是高值离差为正邻居 j 也是高值离差为正乘积就是正的大数如果 i 是高值而邻居是低值乘积就是负数。把这个乘积按照 w_ij 加权求和实际上是在算所有相邻配对之间同向偏离的程度总共有多大。分母Σ(x_i - x̄)^2是大家熟悉的总离差平方和起归一化作用把量纲消掉。前面的n / S_0是缩放因子让指数落在便于比较的范围里。一句话概括莫兰指数就是空间滞后项和原始离差之间的相关系数。你把每个区域的邻居加权平均值也就是空间滞后算出来然后拿它跟区域自己的值求相关系数得到的结果基本就是莫兰指数。这个直观理解在后面画莫兰散点图时特别有用。1.3 期望值不是 0解读莫兰指数最容易踩的认知偏差绝大多数人第一次拿到莫兰指数结果时会下意识地把 0 当成没有空间自相关的基准线。这是个需要纠正的习惯。在完全随机没有空间结构的零假设下莫兰指数的期望值不是 0而是E[I] -1 / (n - 1)n 是区域数量。如果你分析的是 100 个区县期望值是 -0.0101如果是 20 个区县期望值就是 -0.0526。样本越小这条基准线偏离 0 越远越不能忽略。正确的判断方式是把算出来的 I 值和这个期望值比大于期望值说明存在正空间自相关相似值聚集小于期望值说明存在负空间自相关相异值相邻类似棋盘格那种交替分布。同时还要看置换检验给出的伪 p 值而不是靠 I 的绝对值大小拍脑袋。实际算出来的 I 值通常在 -1 到 1 之间正自相关的常见取值在 0.2 到 0.8 之间负自相关在现实数据里比较少见能到 -0.3 已经算很强的离散格局了。如果拿到一个 I 接近 0.9 的结果先别高兴很可能是权重矩阵设置出了问题或者变量里混进了强趋势这个后面第 5 节会专门说。判断标准可以整理成一张表情形I 与 E[I] 的关系伪 p 值结论高值扎堆、低值扎堆I 明显大于 E[I] 0.05显著正自相关高低值交替出现I 明显小于 E[I] 0.05显著负自相关分布接近随机I 接近 E[I] 0.05无显著空间结构2. 空间权重矩阵莫兰指数里最容易被随手一带、却最影响结论的环节2.1 邻接矩阵与距离矩阵选哪个不是审美问题很多人写论文或者做报告时权重矩阵那一栏就一句话带过采用 Queen 邻接矩阵。但从我实际做过的情况看权重矩阵的选择对莫兰指数的影响往往比变量本身的选择还大。同一份数据换个矩阵定义结论从显著聚集变成不显著是常有的事。常见的几类定义方式和它们的适用场景权重类型定义方式适合场景主要风险Rook 邻接共享边才算相邻规则格网、行政区划规整边界处邻居过少Queen 邻接共享边或点都算相邻行政区、不规则多边形角落区域邻居数偏多KNN取最近的 k 个邻居点数据、分布不均的样本k 值需要试错距离带阈值内为邻居超出为 0有明确影响半径的现象阈值难定易出现孤岛核函数按距离连续衰减需要平滑效应时参数解释成本高Rook 和 Queen 的区别就像国际象棋里的车和王后车只能走直线王后能走斜线。落在行政区划数据上Queen 通常更合适因为现实中的两个区县只要挨着不管接的是边还是一个角交流都存在。但如果你的研究关心的是跨界面的相邻关系比如流域、道路网络Rook 反而更贴近实际。KNN 是我在点数据比如门店、监测站、病例点上最常用的起手式原因很简单它保证每个样本都至少有一个邻居不会出现孤岛。2.2 行标准化到底做了什么为什么几乎所有人都开它构建完权重后几乎所有的教程都会加一行w.transform r。这行代码在做行标准化把矩阵的每一行除以该行的和让每个区域的邻居权重加起来等于 1。不标准化的权重矩阵有个隐患。假如区域 A 有 3 个邻居区域 B 有 12 个邻居那么在求和时B 的贡献会被放大 4 倍仅仅因为它邻居多。这显然不是我们想要的——邻居多不代表影响力大。行标准化之后每个区域对最终统计量的投票权被拉平了。行标准化还有个附带好处S_0 直接等于 n前面那个n / S_0的缩放因子变成 1公式简化成Σi Σj w_ij (x_i - x̄)(x_j - x̄) I ------------------------------------ Σi (x_i - x̄)^2需要注意的是行标准化之后莫兰指数的取值范围不再严格是 [-1, 1]可能会略微超出。这不影响显著性判断但如果你要拿不同研究之间的 I 值做横向比较权重设置最好保持一致否则比出来的东西没有意义。2.3 距离阈值怎么定一个能直接抄的经验做法用距离带矩阵时阈值定多少是绕不过去的问题。定太小一堆区域没有邻居变成孤岛定太大全图的人都成了邻居空间结构被抹平。我自己的做法是分两步。先跑一遍最小距离诊断找出让每个区域至少有一个邻居的最小阈值这是理论下限。然后在最小距离到最大距离之间取一个让平均邻居数落在 4 到 8 之间的阈值。为什么是 4 到 8这是个经验区间。太少的话局部结构统计不出来太多的话等于没设置空间约束。当然如果你的现象本身有公认的影响半径比如某种传染病的有效传播距离、某个服务设施的服务半径应该优先用业务意义上的距离而不是纯粹从统计角度调参。代码上可以先做一个邻居数量分布诊断import numpy as np from libpysal.weights import DistanceBand, KNN thresholds np.linspace(w_min, w_max, 30) for t in thresholds: w DistanceBand.from_dataframe(gdf, thresholdt, binaryTrue, silence_warningsTrue) counts np.array([len(v) for v in w.neighbors.values()]) print(f阈值 {t:.0f}m 平均邻居数 {counts.mean():.2f} 孤岛数 {sum(counts 0)})跑一遍看哪一行孤岛数为 0 且平均邻居数落在合适区间就是它了。这个诊断过程花不了五分钟但能让后面所有的结论都站得住脚。3. 用 Python 把全局莫兰指数完整跑一遍从投影到莫兰散点图3.1 数据准备阶段最容易被跳过的一步投影拿一份经纬度坐标的 GeoJSON 直接算空间距离是新手最常见的错误。经纬度是角度单位不是长度单位同一个 1 度在赤道和在高纬度代表的距离差了一倍多。任何跟距离有关的计算都必须先投影到平面坐标系。import geopandas as gpd import numpy as np import libpysal from libpysal.weights import Queen, KNN from esda.moran import Moran, Moran_Local # 读入数据字段里要有待分析的变量 gdf gpd.read_file(districts.shp) # 投影根据研究区所处纬度带选择合适的等距投影 # 中国东部常用 EPSG:3857 或自定义 Albers gdf gdf.to_crs(epsg3857) # 检查待分析变量有没有缺失值 print(gdf[price].isna().sum()) gdf gdf.dropna(subset[price]).reset_index(dropTrue) y gdf[price].values缺失值一定要处理干净。莫兰指数对缺失值非常敏感一个 NaN 会让整个权重矩阵的邻居关系错位算出来的结果完全是垃圾。我一般会直接删掉有缺失的单元同时在报告里说明删了几个、占比多少。3.2 构建权重矩阵并做孤岛检查# 方式一Queen 邻接适合行政区划 w Queen.from_dataframe(gdf, use_indexFalse) # 检查孤岛没有邻居的区域 print(孤岛数量, len(w.islands)) print(孤岛编号, w.islands) # 方式二如果存在孤岛改用 KNN 兜底 w KNN.from_dataframe(gdf, k6) # 行标准化 w.transform r孤岛是必须处理的。它的存在意味着某些区域对空间结构完全没有贡献而且会让某些局部统计量失效。如果孤岛只有一两个可以考虑手工指定邻居或者换成 KNN。如果孤岛很多说明你的研究单元划分本身就存在问题可能需要重新考虑尺度。3.3 计算全局莫兰指数与置换检验莫兰指数的显著性检验不用正态近似而是用置换检验permutation test。做法是把变量值在所有区域之间随机打乱重新算一次 I重复几百到几千次看真实 I 值在随机分布里排在第几百分位。mi Moran(y, w, permutations999) print(fMorans I {mi.I:.4f}) print(f期望值 {mi.EI:.4f}) print(f方差 {mi.VI_sim:.6f}) print(fz 分数 {mi.z_sim:.4f}) print(f伪 p 值 {mi.p_sim:.4f}) print(f正态 p 值 {mi.p_norm:.4f})置换次数 999 是个折中值能得到 0.001 的分辨率。如果你的结论卡在 0.05 附近建议把置换次数提到 9999 再看一次。伪 p 值和正态 p 值如果差异很大说明数据的分布偏离正态比较严重这时候应该更信任伪 p 值。3.4 画莫兰散点图把统计量变成看得懂的图形莫兰散点图是我最推荐的呈现方式因为它能同时告诉你三件事整体自相关的强度、哪些区域是异常点、异常出现在哪个方向。横轴是标准化的变量值 z纵轴是它的空间滞后 Wz邻居的加权平均。拟合直线的斜率就是莫兰指数。import matplotlib.pyplot as plt from splot.esda import moran_scatterplot fig, ax moran_scatterplot(mi, aspect_equalTrue) plt.xlabel(标准化后的价格) plt.ylabel(空间滞后) plt.show()四个象限的含义要记牢右上HH自己高邻居也高典型的高值聚集核心区左下LL自己低邻居也低典型的低值塌陷区左上LH自己低邻居高被高值包围的低值洼地右下HL自己高邻居低高值孤岛拟合直线的斜率就是莫兰指数。斜率越陡空间聚集越强。散点越分散在四个角落说明空间结构越清晰如果散点大致绕成一团斜率接近 0那就是随机分布。提示莫兰散点图里的象限划分和局部莫兰指数的四分类是同一套逻辑只是散点图展示的是全体的分布形态局部指数给出的是每个区域的具体归属和显著性。4. 局部莫兰指数 LISA把一句整体有聚集拆成每个区域的结论4.1 全局显著不代表每个地方都显著全局莫兰指数只能给你一句话整个研究区存在或不存在空间自相关。但这句话对决策的价值很有限。城市管理者想知道的是具体哪几个街道是热点哪几个是冷点钱该往哪里投。更要命的是全局指数还有掩盖局部差异的风险。一个区域可能东部是强聚集、西部是强离散两者在全局计算中相互抵消最后得出无显著自相关的结论。这种情况在做城市研究时并不罕见。局部莫兰指数Local Morans I常简称 LISA就是来解决这个问题的。它把全局公式拆到每个区域计算该区域与它邻居之间的局部关联强度并为每个区域单独做显著性检验。4.2 四种局部关联类型怎么解读局部莫兰指数输出的核心信息是两个这个区域属于四个象限中的哪一个以及这个归属是否显著。类型含义典型场景HH高值被高值包围核心商圈、房价高地LL低值被低值包围连片的老旧小区、经济洼地LH低值被高值包围被豪宅区包围的老破小HL高值被低值包围孤立的商业中心、飞地HH 和 LL 叫做空间聚集说明这一片区域有共同的结构性因素在起作用。HL 和 LH 叫做空间异常它们往往是值得单独研究的对象——为什么在一片低值区里冒出一个高值可能是某个特殊设施带动也可能是数据质量问题。lisa Moran_Local(y, w, permutations999) # 提取显著区域 sig lisa.p_sim 0.05 gdf[lisa_type] np.where(sig, lisa.q, 0) # 注意 esda 的象限编码1HH, 2LH, 3LL, 4HL label_map {0: 不显著, 1: HH, 2: LH, 3: LL, 4: HL} gdf[lisa_label] gdf[lisa_type].map(label_map) gdf.plot(columnlisa_label, categoricalTrue, legendTrue, cmapcoolwarm)4.3 多重比较校正不做这一步的 LISA 图基本不可信这是我在实际评审中最常看到的一个疏漏。LISA 对每个区域都做了一次假设检验如果你的研究区有 300 个区县就做了 300 次检验。在 0.05 的显著性水平下即使数据完全随机也平均会有 15 个区域被判定为显著。这些全是假阳性。处理办法是做多重比较校正常用的有 Bonferroni 和 FDRBenjamini-Hochberg。Bonferroni 太保守容易把真的信号也压掉FDR 控制在空间分析里更常用它在控制假发现比例的同时保留了更多真实信号。from statsmodels.stats.multitest import multipletests reject, p_adj, _, _ multipletests(lisa.p_sim, alpha0.05, methodfdr_bh) gdf[lisa_type] np.where(reject, lisa.q, 0) print(f校正前显著区域数{(lisa.p_sim 0.05).sum()}) print(fFDR 校正后显著区域数{reject.sum()})两个数字一对比往往会发现校正后显著区域少了一大截。这时候你交出去的图才是能经得起追问的。校正之后如果还是显著的聚集区那才是真正值得写进结论的部分。5. 真实项目里让莫兰指数翻车的五种情况5.1 权重矩阵换一个定义结论直接反转我做过一个关于城市绿地覆盖率的研究同一份数据用 Queen 邻接跑出来 I 0.31p 0.002用 KNNk 4跑出来 I 0.18p 0.073。一个是显著正相关一个是不显著。差别仅仅来自邻居的定义。这不是数据的问题也不是方法的问题而是空间关系的定义本身就是研究假设的一部分。稳妥的做法是主分析用一种权重稳健性检验至少再跑两种把结果放在一起对比。如果三种权重下结论一致那结论就很扎实如果结论随权重变化那就要诚实地报告这个敏感性而不是挑一个好看的数报上去。5.2 尺度问题换个划分方式故事完全不一样可修改面积单元问题MAUP是空间分析里的老大难。同样一块地按街道划分和按社区划分算出来的莫兰指数可能相差很大。原因有两层一是不同尺度的空间过程本来就不一样街道层面的聚集机制和社区层面不同二是即使过程相同单元划分的边界也会影响权重结构。我的建议是至少做两个尺度。比如先按区县做一遍再按街道做一遍。如果两个尺度的方向和量级大致吻合说明这个空间结构是稳健的如果不吻合那反而是有价值的发现——它说明这个现象有尺度依赖性本身就是一个值得写进讨论的点。另外要注意边界效应。位于研究区边缘的区域邻居只算了一半它们的空间滞后值会被系统性低估。如果显著区域大量集中在研究区边界先别急着解释很可能是边界效应在作祟。常见的应对是加一圈缓冲带或者在做局部统计时对边界区域的结果保持审慎。5.3 变量本身有强趋势时莫兰指数会虚高假如你分析的是坡度数据或者一个从市中心向外单调递减的人口密度场算出来的莫兰指数会非常高甚至超过 0.9。但这不代表存在真正的空间聚集它只是反映了一个大尺度的梯度趋势。趋势和聚集是两回事。趋势是全局性的、方向一致的聚集是局部性的、围绕均值的波动。莫兰指数的公式是围绕全局均值算离差的一旦变量存在强趋势均值本身就不能代表基准水平结果自然失真。处理办法是做去趋势detrending先用一个一阶或二阶趋势面拟合变量取残差再做莫兰分析。或者直接换用能分离趋势的空间统计量。判断是否存在强趋势很简单把变量按某个坐标轴排序画个散点图如果呈现明显的单调形态就该考虑去趋势了。5.4 置换次数不足p 值在临界点反复横跳置换检验是基于随机模拟的结果本身有随机性。默认 999 次时最小可分辨的伪 p 值是 0.001但如果你的真实 p 值大概在 0.04 到 0.07 之间每次跑出来的结果可能都不一样一会儿显著一会儿不显著。这种做法在正式结论里是危险的。我的习惯是初步探索用 999 次最终结论用 9999 次并且在代码开头设一个固定的随机种子保证结果可复现。import numpy as np np.random.seed(42) mi Moran(y, w, permutations9999)5.5 样本量太小期望值偏移不能忽略前面提过期望值是 -1/(n-1)。如果你研究的是 15 个地级市期望值是 -0.0714这个偏移就不能当没看见。算出来 I -0.05看起来是负数好像有负自相关但跟期望值一比其实是大于期望值的方向是正自相关。小样本还有个附带问题置换检验的分辨率受限。15 个区域的排列组合有限999 次置换里其实存在大量重复实际有效信息量远低于 999。这时候要么扩大研究单元数量要么在报告里明确说明样本量限制不要给出过于肯定的结论。6. 莫兰指数之外还有哪些空间统计工具值得放进工具箱6.1 Getis-Ord Gi*专门用来找热点莫兰指数能告诉我们存在聚集但它对高值聚集和低值聚集的区别在全局层面是混在一起的。Getis-Ord Gi* 则是直接为热点识别设计的它计算每个区域及其邻居的加权和然后跟随机分布比较输出一个 z 分数。z 分数大于 1.96 就是统计意义上的热点高值聚集小于 -1.96 就是冷点。这个输出比 LISA 的四个象限更直观在给非技术背景的决策者做汇报时特别好用。缺点是它只能区分高低不能像 LISA 那样识别出 HL、LH 这类空间异常。我的做法通常是两个一起跑用 LISA 找异常点那些值得单独调查的地方用 Gi* 找热点区那些需要整体投入的区域。6.2 把空间自相关塞进回归模型空间滞后与空间误差如果检验发现残差存在显著空间自相关普通回归的结果就不可信了。这时候有两个主流方向空间滞后模型SAR适合邻居的值直接影响我的值的情形。比如相邻区域的房价会相互参考这是一种实质性的空间溢出。模型里会显式加上一项 Wy。空间误差模型SEM适合存在某些遗漏的、本身有空间结构的变量的情形。比如某个没被纳入模型的政策因素它在空间上聚集从而让误差项也聚集。模型通过误差项的空间结构来吸收这部分影响。怎么选跑一次拉格朗日乘子检验LM test看 LM-Lag 和 LM-Error 哪个更显著或者用稳健版本的检验。两个都显著的话再结合业务逻辑判断。6.3 时空莫兰指数把时间维度叠上去单一时间截面的莫兰指数只反映一个瞬间。如果你的数据有时间维度比如连续十年的监测数据可以算时空莫兰指数把空间权重矩阵和时间权重矩阵做克罗内克积得到一个时空权重矩阵。这样做能回答一些单截面回答不了的问题这个聚集区是在扩大还是缩小某个区域的类型有没有发生转变从 HH 变成 HL这类分析在城市扩张、疾病传播、环境污染研究中价值很高。代价是计算量。时空权重矩阵的维度是 n×T 的平方几万个单元就跑不动了。实际中一般会选择几个关键时间点比如每隔三年取一个或者用滑动窗口的方式来降维。最后分享一点个人体会。我在刚接触莫兰指数的时候花了太多时间纠结公式推导却低估了权重矩阵和尺度选择的重要性。后来做得多了才明白莫兰指数本身的计算并不难难的是让空间关系这个抽象概念真正贴合你研究的现象。邻居到底是按行政边界算还是按实际通勤联系算距离阈值取 3 公里还是 5 公里这些看似技术性的决定背后都是对现象的假设。每次动手之前先花十分钟想清楚我研究的东西影响是怎么传播的比后面反复调参数有用得多。