简介这份ArcGIS教程文档面向地理信息、环境科学与资源管理领域的学习者与从业者系统讲解协同克里金法这一空间插值技术。资源包内含1个docx文件约43KB以图文结合的方式梳理方法原理与操作要点。教程从自相关性与互相关性出发说明如何借助协变量Z2、Z3提升目标变量Z1的预测精度并介绍半变异函数、协方差与互协方差等关键统计量的作用。内容还覆盖普通协同克里金法的模型假设与随机误差处理并延伸至泛协同克里金法、简单协同克里金法、指示协同克里金法、概率协同克里金法与析取协同克里金法等变体同时涉及趋势移除与测量误差应对策略。目前已有3249人学习适合希望深入掌握空间插值与预测、提升地理国情监测与资源评估分析能力的中高级用户参考。1. 协同克里金法当你的采样点不够用时怎么借力打力做空间插值的人迟早会撞上同一个尴尬手头只有二十几个土壤采样点却要出一张全场区重金属分布图。普通克里金跑出来的面要么平滑得像被熨斗烫过要么在点位稀疏处直接摆烂。这时候协同克里金法Co-Kriging就该上场了——它的核心思路很朴素既然主变量采样贵、采样少那就找一个采样便宜、采样密、且和主变量有相关性的辅助变量来“帮衬”一把。在 ArcGIS 的 Geostatistical Analyst 里这套逻辑被封装成了可交互的向导但参数背后的门道远比点几下“下一步”复杂。这篇笔记面向已经会用普通克里金、想进一步压榨有限采样点价值的人从协变量怎么选、半变异函数怎么配、到交叉验证怎么读一步步把协同克里金法在 ArcGIS 里跑通。2. 协同克里金法凭什么比普通克里金多算一层协变量与互变异函数协同克里金法不是把两个变量简单加权平均它是在统计学的框架下利用主变量和协变量之间的空间互相关性来约束预测。理解这一点才能明白为什么有时候加了协变量反而更差。2.1 主变量与协变量的相关性门槛协同克里金法的前提是主变量和协变量之间存在显著的相关性。这个相关性不是随便拿两个图层算个皮尔逊系数就完事它要求的是空间上的共变关系——在空间上距离近的地方主变量和协变量的值变化趋势要一致。我一般会先做两件事用“波段集统计”工具算主变量和协变量的相关系数。经验上相关系数绝对值低于 0.5 时协变量带来的信息增益非常有限强行加入反而增加模型复杂度。画散点图矩阵肉眼确认是否存在非线性关系。如果散点图呈现明显的曲线形态直接做协同克里金会低估协变量的贡献需要先对协变量做变换比如对数变换。注意相关系数高不代表空间共变强。两个变量可能全局相关但在局部空间上各走各的。所以下一步的互变异函数分析不能跳过。2.2 互变异函数协同克里金法真正的“黑匣子”普通克里金只需要拟合一个半变异函数协同克里金需要拟合三个主变量的自变异函数、协变量的自变异函数、以及主变量和协变量之间的互变异函数。互变异函数刻画的是当主变量在某个位置取某值时协变量在距离 h 处取值的协方差。在 ArcGIS 的 Geostatistical Analyst 向导里互变异函数不会单独弹出一个面板让你拟合它是在你选定协变量之后由“半变异函数/协方差建模”工具统一处理的。但你需要理解它背后的约束线性模型共区域化Linear Model of Coregionalization, LMCArcGIS 默认使用 LMC 来保证三个变异函数之间的正定性。简单说就是三个变异函数必须共享同一组变程和块金值只是基台值不同。如果你在拟合时发现某个变异函数的变程和其他两个差太多说明协变量可能不适合。互变异函数的形状理想情况下互变异函数应该从原点附近开始上升到达某个变程后趋于平稳。如果它出现负值或者剧烈震荡通常意味着主变量和协变量在空间上存在相位差这时候要么换协变量要么对协变量做空间平移。2.3 在 ArcGIS 里配置协同克里金的最小步骤下面是一套我常用的操作流程假设你已经有一个主变量点图层比如土壤重金属含量和一个协变量栅格比如高程或植被指数。第一步准备协变量栅格协变量必须是栅格格式且与主变量点图层在同一坐标系下。如果协变量是矢量先用“要素转栅格”或“克里金插值”把它栅格化。# 用 ArcPy 检查主变量点图层和协变量栅格的坐标系是否一致 import arcpy main_points rC:\data\soil_heavy_metal.shp covariate_raster rC:\data\elevation.tif # 获取坐标系 sr_points arcpy.Describe(main_points).spatialReference sr_raster arcpy.Describe(covariate_raster).spatialReference if sr_points.name ! sr_raster.name: print(f坐标系不一致点图层为 {sr_points.name}栅格为 {sr_raster.name}) # 实际项目中应统一投影这里仅做检查 else: print(坐标系一致可以继续)这段代码的作用是避免最常见的翻车场景点图层是地理坐标系栅格是投影坐标系插值结果会整体偏移。参数上arcpy.Describe返回的spatialReference对象包含name属性直接比较名称即可。如果名称不同不要急着用“投影”工具先确认哪个坐标系才是你分析所需的目标坐标系。第二步启动协同克里金向导在 ArcGIS Pro 中路径是分析选项卡 → 地统计向导 → 克里金法/协同克里金法 → 选择“协同克里金法”。在 ArcMap 中路径是自定义 → 扩展模块 → 勾选 Geostatistical Analyst → 地统计向导 → 克里金法 → 协同克里金法。第三步指定主变量和协变量在向导的“输入数据”页面主变量选择你的点图层和对应的字段协变量选择栅格图层。这里有一个容易忽略的选项“协变量类型”。如果协变量是连续型如高程选“连续”如果是分类型如土地利用类型选“分类”。分类协变量会生成指示变异函数计算量更大。第四步拟合变异函数点击“下一步”进入变异函数拟合页面。你会看到三个标签页主变量、协变量、互变异函数。ArcGIS 会自动计算经验变异函数点你需要手动调整模型类型球状、指数、高斯等和参数块金值、基台值、变程。我一般会这样做先看主变量的变异函数确定一个合理的变程范围。再看协变量的变异函数如果它的变程和主变量差一个数量级考虑对协变量做聚合或重采样。最后看互变异函数如果它没有明显的空间结构比如一条水平线说明协变量在这个尺度上帮不上忙。第五步交叉验证拟合完成后向导会要求你进行交叉验证。重点看三个指标均方根误差RMSE、平均标准误差MSE、均方根标准化误差RMSSE。RMSSE 越接近 1 越好如果大于 1 说明预测不确定性被低估小于 1 说明被高估。3. 协变量选不对整个模型白费三个筛选原则与两个替代方案协变量的选择直接决定协同克里金法的成败。我见过太多人随便拿一个栅格就往里塞结果交叉验证的 RMSSE 飙到 2.0 以上还以为是变异函数没拟合好。3.1 原则一协变量必须覆盖整个预测区域这是最硬性的要求。协同克里金法在预测时需要知道每个预测位置上的协变量值。如果你的协变量栅格只覆盖了采样区的一半另一半的预测值会直接报错或输出空值。常见做法是在插值前用“按掩膜提取”工具把协变量裁剪到和主变量点图层的外接矩形一致的范围。如果协变量本身有缺失值考虑用“焦点统计”做邻域填充或者换一个覆盖完整的协变量。3.2 原则二协变量的空间分辨率不宜过粗协变量的栅格分辨率如果比采样点间距还粗那它提供的信息在空间上已经被平均掉了。举个例子采样点平均间距 500 米协变量栅格分辨率 1000 米那么每个栅格像元里可能包含多个采样点协变量值无法区分这些点之间的差异。经验规则协变量栅格分辨率应小于等于采样点平均间距的一半。如果达不到可以对协变量做降尺度比如用样条插值但降尺度不会创造新信息只是让栅格看起来更细。3.3 原则三协变量与主变量的相关性要在局部成立全局相关系数高不代表局部相关性强。我一般会做“移动窗口相关分析”用一个圆形窗口在空间上滑动计算每个窗口内主变量和协变量的相关系数然后看这个相关系数的空间分布。如果某些区域相关系数很低甚至为负说明协变量在这些区域可能引入噪声。在 ArcGIS 里可以用“焦点统计”配合“波段集统计”近似实现但更直接的方法是在 R 或 Python 里用mgcv或scipy做局部相关。如果发现局部相关不成立可以考虑分区建模或者换一个在整个区域都稳定的协变量。3.4 替代方案一用回归克里金代替协同克里金如果协变量和主变量的关系是线性的而且你不想折腾互变异函数回归克里金Regression Kriging是更省事的替代方案。它的思路是先用普通最小二乘或广义最小二乘建立主变量对协变量的回归然后对回归残差做普通克里金最后把回归预测值和残差插值结果相加。在 ArcGIS 里没有直接的回归克里金工具但可以用“普通最小二乘回归”“克里金法”分两步实现。好处是回归部分可以用多个协变量而且不需要拟合互变异函数。3.5 替代方案二用经验贝叶斯克里金处理不确定性如果协变量本身也有测量误差协同克里金法会低估预测的不确定性。这时候可以考虑经验贝叶斯克里金EBK它通过子集模拟来估计变异函数的不确定性对协变量的误差也更稳健。ArcGIS Pro 的 Geostatistical Analyst 里有 EBK 工具但它不支持协变量所以这条路只适用于主变量单独建模的场景。4. 避坑指南协同克里金法在 ArcGIS 里的五个血泪教训这一章记录的是我在实际项目中踩过的坑每个都按“现象 → 原因 → 解决”写清楚。4.1 交叉验证 RMSSE 远大于 1现象交叉验证结果里均方根标准化误差RMSSE达到 1.8 甚至 2.0 以上预测面看起来也过于平滑。原因最常见的原因是变异函数的块金值设得太小。协同克里金法在拟合互变异函数时如果主变量的块金值被低估模型会过度相信采样点之间的空间自相关导致预测不确定性被低估。解决回到变异函数拟合页面手动把主变量和协变量的块金值调大重新拟合。如果调整后 RMSSE 仍然偏高检查协变量是否在局部与主变量负相关——负相关区域会拉低整体预测精度。4.2 预测结果出现条带状异常值现象输出的预测栅格上出现与协变量栅格纹理一致的条带或方块。原因协变量栅格本身存在条带噪声比如遥感影像的扫描条带协同克里金法把这个噪声当成了空间信号。解决在插值前对协变量做滤波处理。可以用“焦点统计”做均值滤波或者用“主成分分析”提取协变量的主要成分丢弃高频噪声成分。如果条带是系统性的考虑换一个数据源。4.3 向导卡在“计算互变异函数”不动现象点击“下一步”后ArcGIS 长时间无响应进度条卡在计算互变异函数。原因主变量采样点数量过多比如超过 5000 个或者协变量栅格分辨率过高导致互变异函数的计算量爆炸。解决对主变量做空间抽稀比如用“创建渔网”“空间连接”保留每个网格内的一个点。或者对协变量做重采样降低分辨率。协同克里金法本身对采样点数量不敏感但对协变量的像元数量很敏感。4.4 分类协变量导致内存溢出现象使用土地利用类型作为协变量时ArcGIS 报内存不足错误。原因分类协变量会为每个类别生成一个指示变异函数如果类别数超过 10 个计算量呈指数增长。解决先对分类协变量做重分类把类别数压缩到 5 个以内。或者改用连续型协变量比如用“欧氏距离”工具计算每个点到某类用地的距离作为连续变量输入。4.5 预测范围超出协变量覆盖区域现象预测栅格在协变量边缘区域出现明显的边界效应值突然跳变。原因协同克里金法在协变量缺失的区域无法计算互变异函数只能退化为普通克里金导致边界处预测值不连续。解决在插值前用“按掩膜提取”把预测范围严格限制在协变量覆盖区域内。如果必须预测边缘区域考虑用“外推”方法单独处理或者换一个覆盖更完整的协变量。5. 从交叉验证到成品图协同克里金法的参数微调与出图技巧跑通协同克里金法只是第一步真正决定成果能不能用的是交叉验证之后的参数微调和出图设置。这一章讲几个我常用的进阶技巧。5.1 用交叉验证的散点图判断系统性偏差ArcGIS 的交叉验证结果里有一个“预测值 vs 实测值”的散点图。如果散点整体偏离 1:1 线说明模型存在系统性偏差。具体来说散点整体在 1:1 线下方模型高估了低值区域低估了高值区域。这时候可以尝试对主变量做对数变换或者改用“普通克里金”的“对数”选项。散点整体在 1:1 线上方模型低估了低值区域高估了高值区域。这种情况比较少见通常意味着协变量和主变量的关系是非单调的。我一般会把这个散点图导出为表格在 Excel 里算一下回归斜率。斜率在 0.9 到 1.1 之间可以接受超出这个范围就需要调整变异函数模型。5.2 变异函数模型的选型对比ArcGIS 提供了球状、指数、高斯、线性等多种变异函数模型。协同克里金法对模型类型的选择比普通克里金更敏感因为三个变异函数必须共享同一组变程和块金值。我通常的做法是先用“自动拟合”让 ArcGIS 选一个模型然后手动改成球状模型再跑一次。球状模型在空间插值中是最稳健的选择它对变程附近的突变不敏感。如果球状模型的交叉验证结果比自动拟合差很多再考虑指数模型或高斯模型。下面是一个用 ArcPy 批量比较不同变异函数模型的脚本框架# 批量比较不同变异函数模型的交叉验证结果 import arcpy from arcpy.ga import * # 输入数据 in_features rC:\data\soil_heavy_metal.shp z_field Pb covariate rC:\data\elevation.tif out_folder rC:\data\ck_results # 模型列表 models [SPHERICAL, EXPONENTIAL, GAUSSIAN] for model in models: out_layer f{out_folder}\\ck_{model}.lyr try: # 执行协同克里金法指定变异函数模型 # 注意ArcPy 的协同克里金法工具需要先创建地统计图层 arcpy.ga.GACoKriging( in_features, z_field, covariate, out_layer, FULL, model, # 模型类型 , , # 其他参数省略实际使用时需补全 ) print(f{model} 模型执行成功) except Exception as e: print(f{model} 模型失败{e})这段代码的逻辑是遍历三种变异函数模型分别生成协同克里金图层。参数上GACoKriging的第四个参数是输出图层路径第五个参数是搜索方向“FULL”表示各向同性第六个参数是变异函数模型类型。实际使用时需要补全块金值、基台值、变程等参数或者先用“自动拟合”生成初始值。5.3 出图时的分级与配色协同克里金法的输出是一个连续栅格直接出图往往看不出空间格局。我一般会做两步处理分级显示用“重分类”工具把预测栅格分成 5 到 7 级分级方法选“自然间断点”或“分位数”。自然间断点适合识别高值聚集区分位数适合展示整体分布形态。配色方案主变量如果是污染物含量用“红-黄-绿”渐变色红色代表高值。协变量如果是高程用“蓝-绿-棕”地形色。出图时把协变量作为底图半透明叠加预测栅格能直观看出预测值和协变量的空间对应关系。5.4 一个我常用的验证习惯每次跑完协同克里金法我都会留出 20% 的采样点不参与插值作为独立验证集。具体做法是在插值前用“子集要素”工具把点图层分成训练集和验证集插值完成后用“提取多值至点”工具把预测值提取到验证集上然后算 RMSE 和 R²。这个习惯帮我避免了很多次“交叉验证好看但实际不能用”的翻车。交叉验证用的是所有点独立验证用的是没参与建模的点后者更能反映模型在未知区域的预测能力。如果独立验证的 RMSE 比交叉验证的 RMSE 高出一倍以上说明模型过拟合了需要简化变异函数模型或者减少协变量。希望这些参数和习惯能帮到你少走几个我当年走过的弯路。本文还有配套的精品资源点击获取