1. 引言与背景

空间自相关分析是地理信息科学（GIScience）和空间统计学的基石。自Tobler（1970）提出“地理学第一定律”——“任何事物都与其他事物相关，但近处的事物比远处的事物更相关”——以来，如何量化空间上的依赖性与异质性便成为核心议题。传统单变量空间自相关（如Global Moran's I）用于检测单一变量是否在空间上呈现聚集、随机或离散格局。然而，现实世界中的地理现象往往不是孤立存在的：城市热岛效应与PM2.5浓度是否在空间上协同增强？新冠感染率与人口密度、医疗资源可达性是否存在空间错位？生态系统供给服务与调节服务之间是否存在权衡或协同的空间格局？这些问题均需要二元空间自相关分析来解答。

二元空间自相关系数（Bivariate Spatial Autocorrelation Coefficient）最早由Anselin等人（2002）在探索性空间数据分析（ESDA）框架下提出，作为单变量Moran's I的直接扩展。其核心思想是：对于某一空间单元i，考察变量x在该单元的观测值，与变量y在i的邻居单元的观测值之间的相关性。这种“变量x在i，变量y在邻居”的设定，使得二元空间自相关能够捕捉跨变量的空间溢出效应。近年来，随着多源地理大数据（如遥感、社交媒体、移动定位数据）的涌现，二元及多元空间自相关方法在生态学、流行病学、城市科学、犯罪地理学等领域获得了广泛关注和重要进展。

本文旨在提供一份关于二元空间自相关系数计算的全面技术指南。我们将从基础概念出发，逐步深入到理论推导、算法实现、显著性检验，并介绍最新的多尺度扩展方法。文章还特别关注了国内外研究者在方法创新上的贡献，例如基于图卷积的神经网络空间自相关、时空二元LISA、以及针对零膨胀数据的修正统计量。通过丰富的案例和可复现的代码，我们希望读者能够快速掌握并应用这一强大的分析工具。

2. 空间自相关基础概念

在深入二元情形之前，有必要回顾空间自相关的基本要素。空间自相关描述的是变量在空间上的取值与其邻近位置取值之间的依赖关系。其数学基础依赖于空间权重矩阵（Spatial Weights Matrix, W），该矩阵定义了空间单元之间的邻接或距离关系。

2.1 空间权重矩阵

空间权重矩阵W是一个n×n的对称矩阵（通常行标准化），其中元素w_ij表示单元i和j之间的空间关系。常见的构建方式包括：

• 邻接权重（Contiguity）：基于多边形共享边界（Rook邻接：共享边；Queen邻接：共享边或顶点）。适用于面状数据。
• 距离权重（Distance-based）：基于质心或点坐标的欧氏距离。常用形式包括反距离权重（w_ij = 1/d_ij^α）或阈值距离（d_ij < d_max时w_ij=1，否则为0）。
• K近邻权重：每个单元选择最近的K个邻居，保证每个单元有相同数量的邻居。

行标准化（row-standardization）是最常见的处理方式，即每行元素之和为1。标准化后的权重矩阵具有明确的解释：w_ij表示邻居j对i的平均影响权重。

2.2 单变量全局Moran's I

单变量全局Moran's I是最经典的空间自相关统计量，其定义为：

其中z_i = x_i - x̄，为变量x的离差；S₀ = Σ_i Σ_j w_ij，为所有权重的总和。I的取值范围理论上在[-1, 1]之间（但实际范围可能不同），正值表示正空间自相关（相似值聚集），负值表示负空间自相关（相异值聚集），零表示空间随机分布。

单变量局部Moran's I（Local Indicators of Spatial Association, LISA）则进一步分解全局指标，识别每个空间单元对全局自相关的贡献，并可以绘制LISA聚类图（HH、HL、LH、LL四种类型）。

3. 二元空间自相关系数理论

二元空间自相关系数将上述框架扩展至两个变量。其核心思想不再是考察同一变量的空间依赖，而是考察变量x在位置i的值与变量y在位置i的邻居j的值之间的相关性。这一设计使得我们能够识别“变量x的高值是否倾向于与变量y的高值在空间上相邻”（协同聚集），或者“变量x的高值是否与变量y的低值相邻”（权衡/错位）。

3.1 二元全局Moran's I (Bivariate Moran's I)

二元全局Moran's I的公式为：

其中z_i^x和z_j^y分别表示变量x在i的离差和变量y在j的离差。注意，求和中的下标i和j：i遍历所有空间单元，j遍历i的邻居。该统计量本质上衡量的是x的观测值与其邻居y的观测值之间的相关性。如果I_xy > 0，则表明x的高值区域周围倾向于聚集y的高值（或x的低值周围聚集y的低值），即空间协同；若I_xy < 0，则表明x的高值区域周围倾向于聚集y的低值（或反之），即空间权衡/错位。

一个重要的变体是交叉Moran's I（Cross-Moran's I），其公式为：

该统计量更接近传统的Pearson相关系数的空间版本，但分母不同。实践中，二元Moran's I更常用。

3.2 二元局部Moran's I (Bivariate LISA)

局部版本对于识别空间异质性至关重要。对于每个空间单元i，二元局部Moran's I定义为：

对所有的局部I_i^xy求和（乘以常数）即得到全局统计量。局部值可正可负，其符号和大小揭示了位置i在双变量空间关联中的角色。通过将z_i^x的符号与空间滞后项Σ_j w_ij z_j^y的符号组合，可以定义四种二元LISA类型：

• HH (High-High)：z_i^x > 0 且 Σ_j w_ij z_j^y > 0，表示x的高值区域被y的高值邻居包围（协同热点）。
• HL (High-Low)：z_i^x > 0 且 Σ_j w_ij z_j^y < 0，表示x的高值区域被y的低值邻居包围（错位/孤岛）。
• LH (Low-High)：z_i^x < 0 且 Σ_j w_ij z_j^y > 0，表示x的低值区域被y的高值邻居包围（冷点中的热点）。
• LL (Low-Low)：z_i^x < 0 且 Σ_j w_ij z_j^y < 0，表示x的低值区域被y的低值邻居包围（协同冷点）。

值得注意的是，二元LISA的解读与单变量LISA有本质区别：单变量LISA中的HH表示该单元自身高值且邻居高值，而二元LISA中的HH表示该单元的x高值且邻居的y高值。因此，二元LISA揭示的是跨变量的空间共现模式。

3.3 二元Geary's C

作为Moran's I的补充，Geary's C更侧重于度量空间上的差异而非协方差。二元Geary's C定义为：

注意，该公式中直接比较的是x_i与y_j的差异。C_xy的值在0到2之间变化（理论上），0表示完全正相关（x_i与邻居y_j完全相同），1表示无空间关联，2表示完全负相关。实践中，Geary's C对局部变异更为敏感，但不如Moran's I常用。

4. 计算方法与算法实现

计算二元空间自相关系数的核心步骤包括：数据标准化、构建空间权重矩阵、计算交叉乘积、进行显著性检验。以下给出详细的算法流程。

4.1 算法步骤

1. 数据准备：确保两个变量x和y为数值型，且无缺失值。对x和y分别进行Z-score标准化（减去均值，除以标准差），得到z^x和z^y。标准化使得不同量纲的变量可比。
2. 构建空间权重矩阵W：根据数据几何类型（点、线、面）选择合适的邻接或距离规则。对W进行行标准化。
3. 计算空间滞后：对于每个单元i，计算变量y的邻居加权平均（空间滞后）：(W·z^y)_i = Σ_j w_ij z_j^y。
4. 计算全局统计量：
   I_xy = (n / S₀) · ( (z^x)^T · W · z^y ) / ( (z^x)^Tz^x )^1/2 · ( (z^y)^Tz^y )^1/2
   注意，这里分母是两个标准差的乘积，分子是交叉空间协方差。
5. 计算局部统计量：对于每个i，I_i^xy = z_i^x · (W·z^y)_i。

4.2 伪代码与计算复杂度

以下为计算二元全局Moran's I的伪代码：

Input: vectors x, y (length n), weight matrix W (n x n, row-standardized)
Output: I_xy (global bivariate Moran's I)

1. zx = (x - mean(x)) / std(x)
2. zy = (y - mean(y)) / std(y)
3. wy_sum = W %*% zy   # spatial lag of y
4. numerator = t(zx) %*% wy_sum
5. denom = sqrt(t(zx) %*% zx) * sqrt(t(zy) %*% zy)
6. S0 = sum(W)   # or n if row-standardized
7. I_xy = (n / S0) * numerator / denom
8. return I_xy

计算复杂度主要由矩阵乘法W·z^y决定，为O(n²)。对于大规模数据（n > 10万），可使用稀疏矩阵表示W（如scipy.sparse），将复杂度降至O(kn)，其中k为平均邻居数。局部统计量的计算复杂度同样为O(kn)。

4.3 多变量扩展与广义二元指标

近年来，研究者提出了多种扩展形式。例如，多元Moran's I（Multivariate Moran's I）通过将多个变量向量化，计算向量之间的空间协方差矩阵的迹。具体地，若我们有p个变量，构成n×p矩阵X，则多元Moran's I定义为：

该指标可以看作是所有变量对之间二元Moran's I的加权平均。另一种重要的扩展是偏空间自相关（Partial Spatial Autocorrelation），用于在控制第三个变量的条件下，衡量两个变量的空间共变。

5. 多尺度与时空扩展

传统二元空间自相关假设空间关系是全局同质的（即权重矩阵不随尺度变化）。然而，地理过程往往具有多尺度特征。例如，城市热岛效应与植被覆盖的关系在街区尺度（100m）和区域尺度（10km）上可能截然不同。为此，学者们发展了多尺度二元空间自相关方法。

5.1 多尺度空间权重矩阵

通过构建不同尺度的空间权重矩阵（如不同距离阈值、不同阶数的邻接关系），可以计算尺度依赖的二元Moran's I。设W_d为距离阈值d下的权重矩阵，则尺度d下的二元Moran's I为I_xy(d)。通过绘制I_xy(d)随d变化的曲线（称为空间相关图，Spatial Correlogram），可以识别出双变量空间关联的特征尺度。

进一步地，多尺度地理加权回归（MGWR）的思想也可迁移至此：每个位置i可以具有最优的带宽，从而得到局部尺度的二元LISA。Fotheringham等（2017）提出的多尺度GWR框架，通过后向拟合算法为每个变量选择不同的带宽。类似地，我们可以为变量x和y分别选择不同的空间权重尺度，然后计算跨尺度的二元自相关。

5.2 时空二元空间自相关

当数据包含时间维度（如面板数据或时空快照）时，需要引入时空权重矩阵。设数据为n个空间单元在T个时间点的观测值，则时空权重矩阵可以表示为：

其中⊗表示Kronecker积，W_space是n×n的空间权重矩阵，W_time是T×T的时间权重矩阵（例如，只考虑相邻时间点）。然后，将nT个时空观测值视为一个整体，计算时空二元Moran's I。另一种更灵活的方法是使用时空核函数，同时考虑空间距离和时间距离，例如：

其中θ_s和θ_t分别为空间和时间带宽。这种方法的优势在于可以灵活地定义时空邻域。

5.3 非平稳二元空间自相关

空间非平稳性（Spatial Nonstationarity）指空间关系随地理位置变化。传统的全局二元Moran's I假设关系是全局平稳的，而局部LISA虽然可以揭示局部模式，但未对非平稳性进行建模。一种创新的方法是将地理加权回归（GWR）的思想引入二元自相关：对于每个位置i，使用局部加权回归估计x与y的局部空间滞后关系，然后计算局部化的二元自相关统计量。具体地，定义局部二元Moran's I为：

其中w_ij(h)是基于核函数（如高斯核）的局部权重，h是带宽。通过改变h，可以探索不同空间尺度下的局部关系。这种方法被称为局部双变量空间关联指标（Local Bivariate Spatial Association, LBSA）。

6. 显著性检验与推断

计算得到的统计量需要经过显著性检验，以判断观察到的空间模式是否显著异于随机分布。最常用的方法是置换检验（Permutation Test），也称为蒙特卡洛模拟。

6.1 置换检验原理

对于二元全局Moran's I，零假设H₀为：变量x和y在空间上是随机共变的（即观察到的I_xy不显著异于期望值）。在H₀下，变量x的观测值在空间上随机排列，而变量y的观测值固定（或反之），从而打破任何潜在的空间结构。具体步骤：

1. 计算观察值I_xy^obs。
2. 对于b = 1,..., B（通常B=999或9999）：随机打乱变量x的观测值（保持y不变），得到置换后的x^(b)；重新计算I_xy^(b)。
3. 计算p值：p = (1 + count(I_xy^(b) ≥ I_xy^obs)) / (B+1)（单侧上尾检验）。
4. 若p < 0.05，则拒绝H₀，认为存在显著的二元空间自相关。

对于局部指标，同样进行置换检验，但需要注意多重检验问题（n个局部检验同时进行）。常用方法包括Bonferroni校正（p_adj = p / n）或FDR（False Discovery Rate）控制。

6.2 近似分布检验

当样本量较大时，可以基于渐近正态分布进行检验。二元Moran's I的期望值在H₀下为：E[I_xy] = -1/(n-1)（与单变量相同）。方差可通过Cliff-Ord公式计算，但需要假设数据服从多元正态分布。实践中，置换检验更为稳健，尤其适用于非正态分布数据。

6.3 针对特殊数据类型的修正

对于含有大量零值的数据（如疾病发病率、犯罪事件计数），传统的置换检验可能低估方差。近年来，条件自回归模型（CAR）和层次贝叶斯方法被用于对零膨胀数据进行建模。例如，可以在贝叶斯框架下，将二元空间自相关参数作为模型参数进行估计，通过MCMC采样得到后验分布。这种方法虽然计算量较大，但能够更准确地刻画不确定性。

7. 国内外研究前沿

近年来，二元空间自相关领域涌现出大量创新性研究。以下从方法创新和应用拓展两个方面进行梳理。

7.1 方法创新

• 基于图神经网络的二元空间自相关：2023年，Li等人提出了一种基于图卷积网络（GCN）的二元空间自相关学习框架。该方法将空间权重矩阵作为图邻接矩阵，通过GCN学习节点嵌入，然后