ArcGIS协同克里金法实战：从互协方差建模到交叉验证避坑指南

简介：这份ArcGIS教程文档面向地理信息、环境科学与资源管理领域的学习者与从业者，系统讲解协同克里金法这一空间插值技术。内容围绕如何利用目标变量Z1的自相关性及其与协变量Z2、Z3之间的互相关性提升预测精度展开，涵盖普通协同克里金法的模型假设、半变异函数与协方差的量化方式，并延伸至泛协同克里金法、简单协同克里金法、指示协同克里金法、概率协同克里金法与析取协同克里金法等变体，同时涉及测量误差处理、趋势移除与互协方差应用等实操要点。资源包为1个docx文档，大小约43KB，结构紧凑，便于快速查阅与对照学习。目前已有3249人学习下载，适合希望深入理解空间插值原理、提升地理国情监测与资源评估分析能力的中高级用户参考。

1. 协同克里金法在 ArcGIS 里到底怎么落地：从 Z1 和 Z2 的互相关说起

如果你手头有一份主变量采样点稀疏、但旁边还躺着一份采样密集的协变量数据，普通克里金插出来的面往往在局部会出现明显的“条带感”或者高值低估。协同克里金法就是冲着这个场景来的：它不只吃 Z1 的自相关，还把 Z1 与 Z2 之间的互相关一起塞进方程组，用协变量的信息去补主变量的窟窿。在 ArcGIS 里，这套方法落在 Geostatistical Analyst 扩展模块的“协同克里金”工具里，支持普通、简单、泛、指示、概率、析取六种变体。适合做地理国情监测、土壤属性制图、地下水位预测这类“主变量难测、协变量易得”的活。但别急着上，它比普通克里金多估一层互协方差，参数翻倍，翻车概率也翻倍。

2. 协同克里金法的数学底子：为什么多一个变量反而更麻烦

2.1 自相关、互相关与协方差函数的关系

协同克里金的核心假设是：主变量 Z1 和协变量 Z2 各自有均值 µ1、µ2，随机误差 ε1(s)、ε2(s) 既自相关又互相关。普通克里金只估 Z1 的自协方差 C11(h)，协同克里金还要估 Z2 的自协方差 C22(h) 以及两者的互协方差 C12(h)。这三个函数共同决定权重矩阵。用半变异函数表达时，互变异函数 γ12(h) 通常由 C12(h) 转换而来，ArcGIS 内部用协方差更稳，因为半变异函数在互相关为负时会出现非正定问题。

常见做法是先用“半变异函数/协方差建模”工具分别拟合 Z1、Z2 的自变异函数，再拟合互变异函数。注意：互变异函数不能随便套球状模型，它可能取负值，ArcGIS 允许用“协方差函数”直接建模，或者用“线性模型”处理负相关。如果 Z1 低时 Z2 高，互协方差就是负的，这时候硬套指数模型会得到一个正定但完全错误的权重。

2.2 普通协同克里金与泛协同克里金的选型边界

普通协同克里金假设均值未知但恒定，适合趋势不明显的数据。泛协同克里金则把均值当成坐标的函数（一次或二次趋势面），适合有区域趋势的变量，比如降水随海拔升高。简单协同克里金要求均值已知，实际中很少用，除非你做过大量先验调查。指示协同克里金针对分类数据，比如把 Z1 按阈值切成 0/1，再预测某个阈值以上的概率，它对异常值不敏感，但需要每个阈值单独建模。

选型时先看数据：如果 Z1 和 Z2 的散点图呈明显线性且无趋势，普通协同克里金够用；如果 Z1 有南北梯度，泛协同克里金更合适；如果 Z1 是重金属超标与否这类二值变量，直接上指示协同克里金。别为了“更精确”硬上泛协同，多估的趋势参数会吃掉样本量，小数据集上反而更差。

2.3 互协方差估计的代价：参数翻倍与不确定性

协同克里金需要估计的参数量大致是普通克里金的 2 到 3 倍：Z1 自相关（块金、基台、变程）、Z2 自相关、互相关（块金、基台、变程），如果加趋势还要多几个系数。每多估一个参数，估计方差就多一份。理论上即使 C12(h)=0，协同克里金也能退化成普通克里金，但实际拟合时互协方差模型可能“假装”存在相关，导致权重矩阵病态。血泪经验是：先做互相关检验，Pearson 相关系数绝对值低于 0.3 就别硬上协同，直接用普通克里金加回归残差更稳。

3. 在 ArcGIS 里跑通普通协同克里金：从数据准备到交叉验证

3.1 数据准备与坐标系检查

ArcGIS 的协同克里金工具要求输入点要素，主变量和协变量可以在同一图层，也可以分两个图层。如果分图层，两个图层的点位置不必重合，这是协同克里金的优势之一。但坐标系必须一致，且投影坐标系优先，因为变程是基于地图单位的。地理坐标系下变程单位是度，拟合出来的各向异性会失真。

常见做法是先把 Z1 和 Z2 分别做成点要素，用“投影”工具统一到同一投影坐标系，比如 UTM 或 Albers。然后检查 Z1 和 Z2 的直方图，如果严重偏态，先做对数或 Box-Cox 变换。ArcGIS 的协同克里金工具内置了变换选项，但建议在预处理阶段做，方便对比。

# ArcPy 预处理：投影 + 提取分析 import arcpy from arcpy import env env.workspace = r"C:\gis\ck_data.gdb" env.overwriteOutput = True # 投影到 UTM 50N arcpy.Project_management("Z1_raw", "Z1_utm", arcpy.SpatialReference(32650)) arcpy.Project_management("Z2_raw", "Z2_utm", arcpy.SpatialReference(32650)) # 检查点数量与范围 desc = arcpy.Describe("Z1_utm") print("Z1 要素数:", arcpy.GetCount_management("Z1_utm")) print("Z1 范围:", desc.extent)

这段代码做两件事：统一投影、输出要素数和范围。参数arcpy.SpatialReference(32650)是 WGS84 UTM 50N，国内东部常用；如果你在西部，换成对应的 UTM 带号。GetCount_management返回的是字符串，实际用的时候转 int。范围检查是为了确认 Z1 和 Z2 的覆盖区域有重叠，如果 Z2 只在 Z1 外围，互相关估计会不可靠。

3.2 用 Geostatistical Analyst 向导拟合互协方差

ArcGIS 的协同克里金向导分五步：数据选择、趋势移除、半变异函数/协方差建模、搜索邻域、交叉验证。关键在第三步。向导里会分别列出 Z1 的自协方差、Z2 的自协方差、Z1-Z2 的互协方差。互协方差图通常点很散，别追求 R² 高，先看块金效应。如果互协方差块金接近总基台，说明短距离互相关弱，这时候可以把互协方差模型设成纯块金，相当于放弃互相关，退化成普通克里金。

参数设置上，自协方差用球状或指数模型，互协方差用线性模型或高斯模型。变程一般参考自协方差的变程，互相关变程不会比自相关大太多。如果互协方差变程远大于自协方差，多半是趋势没移除干净，回去做一阶趋势移除。

# 调用协同克里金工具（普通协同克里金） arcpy.GACok_ga( in_features="Z1_utm", # 主变量点图层 z_field="Z1_value", # 主变量字段 covariate="Z2_utm", # 协变量图层 covariate_field="Z2_value", # 协变量字段 out_ga_layer="Cok_layer", # 输出地统计图层 out_raster="Cok_raster", # 输出栅格 cell_size=100, # 栅格分辨率 search_radius=5000, # 搜索邻域半径 min_neighbors=5, # 最少邻点数 max_neighbors=15 # 最多邻点数 )

这是 ArcPy 里GACok_ga的简化调用。cell_size根据你的采样密度定，一般取平均点间距的 1/5 到 1/10。search_radius要大于变程，否则邻域内点太少，权重不稳。min_neighbors设 5 是底线，低于 5 会出现无解。max_neighbors设 15 是防止计算量爆炸，协同克里金的矩阵是普通克里金的 2 倍，邻点数超过 20 会明显卡顿。

3.3 交叉验证与精度对比

跑完插值必须做交叉验证。ArcGIS 的“交叉验证”工具会输出预测误差的均值、均方根、平均标准误差等。协同克里金的均方根误差应该比普通克里金低 10% 以上才值得用，否则多估的参数没带来实际收益。注意看“标准化均方根”是否接近 1，如果大于 1.2，说明互协方差模型高估了相关性，预测过于自信。

对比时把普通克里金和协同克里金的交叉验证结果放一张表：

如果协同克里金的均方根只降了 3%，而标准化均方根飙到 1.5，说明互协方差是硬凑的，果断退回普通克里金。

4. 避坑与排查：协同克里金在 ArcGIS 里的五个翻车现场

4.1 互协方差图一片散点，拟合出来全是块金

现象：互协方差云图没有明显上升趋势，拟合模型块金等于基台。原因：Z1 和 Z2 的采样位置完全不重叠，或者重叠极少，互协方差只能靠距离很近的异位置点对估计，噪声大。解决：检查两个图层的点距离，如果最近点对超过变程的 1/3，互相关无法估计。要么补充重叠采样点，要么放弃协同克里金，改用回归克里金（先用 Z2 回归 Z1，再对残差做普通克里金）。

4.2 交叉验证报“矩阵奇异”或“无法求逆”

现象：工具运行到一半报错，提示协方差矩阵非正定。原因：互协方差模型和自协方差模型不兼容，比如自协方差用球状模型，互协方差用高斯模型，导致协方差矩阵出现负特征值。解决：统一模型类型，自协方差和互协方差都用球状或都用指数。如果必须用不同模型，在“半变异函数/协方差建模”里勾选“允许非正定”，但预测结果可能不稳定，只用于探索。

4.3 预测面出现明显“牛眼”或条带

现象：插值结果在某个点周围出现同心圆高值，或者沿某方向条带。原因：搜索邻域设得太小，或者各向异性没设对。协同克里金对邻域更敏感，因为协变量的点会参与权重。解决：把搜索邻域半径调到变程的 1.5 倍，邻点数上限提到 20。如果数据有各向异性，在“搜索邻域”里设方向角和长短轴比，别用默认各向同性。

4.4 协变量比主变量还多，工具直接卡死

现象：选了 Z2、Z3、Z4 三个协变量，点击运行后 ArcGIS 无响应。原因：协同克里金每加一个协变量，互协方差矩阵维度翻倍，三个协变量就是 4×4 的块矩阵，计算量指数上升。解决：ArcGIS 的协同克里金工具一次只支持一个协变量。多个协变量要用“泛协同克里金”或者分步：先选相关性最高的一个，跑完看残差，再对残差用第二个协变量。别指望一次塞三个。

4.5 变换后的数据忘记反变换，预测值量纲错误

现象：预测栅格的值在 0 到 1 之间，但原始 Z1 是 0 到 100。原因：预处理时做了对数变换，插值后没有反变换。ArcGIS 的协同克里金工具如果在内置变换里选了“对数”，输出会自动反变换；但如果你在 ArcPy 里手动arcpy.Log变换了字段，工具不知道，输出就是对数值。解决：手动变换的，插值后自己用arcpy.Exp反变换，或者直接用工具内置的变换选项。

5. 进阶技巧：用指示协同克里金处理阈值预测与异常值

指示协同克里金是协同克里金里最实用的变体之一，尤其当 Z1 有异常值或者你只关心“是否超过某个阈值”。做法是把 Z1 按阈值切成 0/1 指示变量，然后对每个阈值分别做协同克里金，最后得到超过阈值的概率面。它对异常值不敏感，因为异常值只影响它所在的那一个阈值，不会像普通协同克里金那样把整个变异函数拉偏。

具体操作：在 Geostatistical Analyst 向导里选“指示协同克里金”，设置阈值，比如 Z1 的 75% 分位数。协变量 Z2 保持连续。工具会为每个阈值生成一个指示变异函数，互协方差也变成指示变量与 Z2 的互相关。注意：阈值不要设太多，3 到 5 个足够，每多一个阈值就多一套模型，交叉验证工作量翻倍。

# 指示协同克里金：阈值 75% 分位数 import numpy as np z1_values = arcpy.da.FeatureClassToNumPyArray("Z1_utm", "Z1_value") threshold = np.percentile(z1_values["Z1_value"], 75) print("阈值:", threshold) arcpy.GAICok_ga( in_features="Z1_utm", z_field="Z1_value", covariate="Z2_utm", covariate_field="Z2_value", out_ga_layer="ICok_layer", out_raster="ICok_raster", cell_size=100, threshold=threshold, # 指示阈值 search_radius=5000 )

这段代码先算 75% 分位数作为阈值，再调用GAICok_ga。threshold参数是指示协同克里金特有的，普通协同克里金没有。输出栅格的值是概率，0 到 1 之间。验证时看概率面的校准度：如果预测概率 0.8 的区域实际超标频率只有 0.5，说明模型过于自信，需要回到互协方差模型调块金。

我一般会先跑普通协同克里金，如果交叉验证的标准化均方根超过 1.3，或者数据里有明显离群点，就切到指示协同克里金。指示协同克里金的概率面不能直接当连续值用，但做风险区划足够。从那以后我每次做协同克里金之前，都强制先跑一遍互相关散点图和普通克里金基线，只有协同克里金的均方根降低超过 10% 才继续往下走。希望帮到你。

本文还有配套的精品资源，点击获取