企业官网建设流程全解析

简介：这份ArcGIS教程文档面向地理信息、环境科学与资源管理领域的学习者与从业者，系统讲解协同克里金法这一空间插值技术。内容围绕主要变量Z1的自相关性及其与协变量Z2、Z3的互相关性展开，说明如何借助半变异函数、协方差与互协方差量化空间结构，并对比普通协同克里金法与泛协同克里金法、简单协同克里金法、指示协同克里金法、概率协同克里金法、析取协同克里金法等变体的适用条件。文档还涉及测量误差处理、趋势移除及多变量交互影响等实际问题，可帮助读者理解如何利用变量间关系提升Z1预测精度。资源包为1个docx文件，约43KB，结构紧凑，适合作为方法查阅与案例参考。目前已有3249人学习下载，适合希望在地理国情监测、资源评估与环境影响分析中熟练运用协同克里金法的中高级用户。

1. 协同克里金法到底在解决什么麻烦

做过空间插值的人多半有过这种经历：手头的主变量采样点稀稀拉拉，硬跑普通克里金，出来的图要么过度平滑像被熨斗烫过，要么在无数据区放飞自我。协同克里金法（Co-Kriging）就是冲着这个痛点来的——它允许你把一个采样密集、但和主变量有相关性的"辅助变量"拉进来，借它的空间结构去约束主变量的估计。在 ArcGIS 的 Geostatistical Analyst 里，这套东西落地成了 CoKriging 工具，输入一个主变量数据集加一到多个协变量数据集，输出预测栅格和误差面。

它适合谁？适合手里有两套空间数据、且两套数据存在物理或统计相关性的人：比如用高程协变量去插值气温，用遥感植被指数协变量去插值土壤有机质，用密集的电磁测量协变量去插值稀疏的化学采样。不适合谁？主辅变量相关系数低于 0.5 还硬上，或者协变量本身也是稀疏点、覆盖范围还不如主变量——那纯属给自己找麻烦。这一章先把"值不值得用"讲清楚，后面几章再拆怎么在 ArcGIS 里跑通、参数怎么调、坑在哪。

2. 协同克里金法的数学骨架与 ArcGIS 里的落地形态

2.1 从普通克里金到协同克里金：多出来的那块协方差

普通克里金的估计量是主变量自身采样点的加权和，权重由主变量的变异函数决定。协同克里金在此基础上多了一层：估计值不仅用主变量采样点，还用协变量采样点，权重矩阵里因此多出了主-协交叉变异函数（cross-variogram）和协变量自身的变异函数。

用公式说，估计量写成：

Z*(x0) = Σ λi · Z(xi) + Σ μj · Y(xj)

其中 Z 是主变量，Y 是协变量，λ 和 μ 是待求权重。求解权重的方程组里，除了主变量的变异函数 γZZ，还要塞进交叉变异函数 γZY 和协变量变异函数 γYY。交叉变异函数刻画的是"主变量在空间上变化时，协变量怎么跟着变"——它是对称的，γZY(h) = γYZ(h)。

ArcGIS 的 Geostatistical Analyst 在底层用的是线性协同克里金模型（LMC，Linear Model of Coregionalization）。它要求所有变异函数（主、协、交叉）都由同一组基函数结构叠加而成，只是各结构的贡献系数不同。这个约束不是 ArcGIS 独有的，是协同克里金能保证估计方差非负的数学前提。理解这一点很关键：你在界面里调主变量的变异函数模型时，协变量和交叉变异函数的模型结构会被联动约束，不能各调各的。

2.2 在 ArcGIS Pro 里跑通 CoKriging 的最小流程

先说数据准备。主变量和协变量必须是同一投影坐标系下的点要素，属性表里各有一个数值字段。如果协变量是栅格，得先转成点，或者用 Extract Values to Points 把栅格值提取到主变量点上——但注意，后者会让协变量退化成"只在主变量位置有值"，失去协同克里金"借密集辅助点"的核心优势。所以常见做法是：协变量本身也是点数据，且点数明显多于主变量。

下面是用 ArcPy 跑协同克里金的最小脚本。Geostatistical Analyst 的 Python 接口是arcpy.ga模块，核心类是CoKriging。

import arcpy from arcpy.ga import CoKriging # 环境设置：覆盖同名输出，指定工作空间 arcpy.env.overwriteOutput = True arcpy.env.workspace = r"C:\proj\interp\data.gdb" # 输入：主变量点、协变量点、主变量字段、协变量字段 in_primary = "soil_samples" # 主变量：稀疏的土壤有机质采样点 in_covariate = "ndvi_points" # 协变量：密集的植被指数点 z_field = "om_gkg" # 主变量字段：有机质 g/kg cov_field = "ndvi" # 协变量字段：NDVI 值 # 创建协同克里金对象 ck = CoKriging(in_primary, z_field) # 添加协变量，第二个参数是协变量字段名 ck.covariates = [[in_covariate, cov_field]] # 设置搜索邻域：主变量用 12 个点，协变量用 8 个点 ck.searchNeighborhood = "Standard" ck.searchNeighborhoodNumberOfPoints = 12 ck.covariateSearchNeighborhoodNumberOfPoints = 8 # 指定变异函数模型：主变量用球状模型，滞后步长自动 ck.modelType = "SPHERICAL" ck.lagSize = 500 # 滞后步长，单位与投影一致（米） # 输出预测栅格和预测标准误差栅格 out_pred = "om_cokriging_pred" out_err = "om_cokriging_err" ck.execute(out_pred, out_err) print("协同克里金完成，输出：", out_pred, out_err)

这段脚本的逻辑说明：CoKriging对象先绑定主变量，covariates属性接收一个列表，每个元素是[协变量数据集, 字段名]，支持多个协变量。searchNeighborhoodNumberOfPoints控制主变量参与估计的邻居数，covariateSearchNeighborhoodNumberOfPoints单独控制协变量的邻居数——这是协同克里金和普通克里金在参数上最直观的区别，两个数要分开调。

参数说明：lagSize是变异函数拟合的滞后步长，经验值是研究区最大距离的 1/10 到 1/20。设太小，变异函数点云噪声大；设太大，短程结构被抹平。modelType可选 SPHERICAL、EXPONENTIAL、GAUSSIAN 等，球状模型在土壤、气象插值里最常用，因为它有明确的变程（range），超过变程后空间自相关归零，符合多数自然变量的行为。

2.3 变异函数拟合：协同克里金真正的门槛

跑通脚本只是第一步，协同克里金能不能出好结果，八成看变异函数拟合。在 ArcGIS Pro 的 Geostatistical Wizard 里，走到"Model the semivariogram"那一步，你会看到主变量、协变量、交叉变异函数三个面板。关键操作顺序是：

先拟合主变量的变异函数，确定块金值（nugget）、基台值（sill）、变程（range）三个参数。再切到协变量面板，拟合协变量的变异函数。最后看交叉变异函数——它的形状应该介于主变量和协变量之间，如果交叉变异函数出现明显的负值或剧烈震荡，说明主辅变量的空间相关性不成立，这时候硬跑协同克里金还不如退回普通克里金。

ArcGIS 提供自动拟合（Fit）按钮，但自动拟合经常把块金值压得过低，导致预测面出现"牛眼"（bull's-eye）——每个采样点周围一圈同心圆。血泪经验是：自动拟合后手动把块金值往上抬 10% 到 20%，预测面会平滑很多，交叉验证的均方根误差通常也更稳。

3. 参数怎么设：搜索邻域、各向异性与协变量筛选

3.1 搜索邻域：主辅变量分开调的理由

协同克里金的搜索邻域有两个维度：主变量邻居数和协变量邻居数。为什么分开？因为主变量稀疏、协变量密集，如果共用一个邻居数，要么主变量邻居不够导致估计不稳，要么协变量邻居太多导致计算爆炸。

常见做法是：主变量邻居数取 8 到 16，协变量邻居数取主变量的 1.5 到 2 倍。但协变量邻居数不是越多越好——协变量对估计的贡献随距离衰减，超过变程的协变量点权重接近零，加进来只增加计算量不增加精度。我一般先用主变量变程的 1.2 倍作为搜索半径，在这个半径内数协变量点，如果超过 20 个，就限制邻居数为 20。

在 ArcPy 里，搜索邻域还可以设成扇形（Sector）或椭圆（Ellipse），对应各向异性场景。比如河谷地区的土壤变量，沿河流方向的空间相关性明显强于垂直方向，这时候用椭圆邻域，长轴沿河流走向，短轴垂直，比圆形邻域合理得多。

3.2 各向异性：什么时候必须开

各向异性分两种：几何各向异性（geometric）和带状各向异性（zonal）。几何各向异性是空间相关性在不同方向上"变程不同、基台相同"，带状各向异性是"变程和基台都不同"。ArcGIS 的变异函数拟合界面里，勾选"Anisotropy"后可以设主方向角度和长短轴比。

判断要不要开各向异性的土办法：把采样点按不同方向分组，分别算变异函数，如果各方向的变程差异超过 30%，就值得开。开了之后，主方向角度用变异函数点云最"扁"的方向，长短轴比用长轴变程除以短轴变程。

注意一个坑：协同克里金里，主变量和协变量的各向异性方向必须一致，否则交叉变异函数的拟合会崩。如果两个变量的各向异性方向差很多，说明它们的空间结构不是同一套机制驱动的，这时候要么放弃协同克里金，要么只在一个方向占主导的子区域里做。

3.3 协变量筛选：不是越多越好

协同克里金支持多个协变量，但每加一个协变量，变异函数拟合的复杂度就上一个台阶。三个以上协变量时，交叉变异函数的拟合基本靠猜，出来的结果很难解释。

我的筛选流程是：先算主变量和每个协变量的皮尔逊相关系数，只保留 |r| > 0.5 的；再对保留的协变量两两算相关，如果两个协变量之间 |r| > 0.8，只留和主变量相关性更高的那个。这样通常能把协变量压到一到两个，拟合可控，结果也好解释。

4. 避坑与排查：协同克里金翻车的五个典型场景

4.1 预测面出现大面积负值

现象：土壤有机质预测结果里出现 -5 g/kg 这种物理上不可能的值。原因：协同克里金的估计量是线性加权和，数学上不保证非负；当协变量与主变量在某些局部呈负相关，或者变异函数拟合出负的交叉协方差，权重就可能为负，估计值被拉到零以下。解决：在 CoKriging 对象里开启transformType = "LOG"做对数变换，插值完再反变换回原尺度；或者用searchNeighborhood限制邻居数，减少远距离负权重的机会。

4.2 交叉验证误差比普通克里金还大

现象：跑完协同克里金，交叉验证的 RMSE 反而比不用协变量的普通克里金高。原因：协变量和主变量的相关性在研究区内不均匀，某些子区域相关性强、某些弱，全局拟合的交叉变异函数被弱相关区域拖累。解决：分区建模，把研究区按相关性强度切成两到三块，每块单独跑协同克里金；或者退回普通克里金，别为了用而用。

4.3 变异函数拟合时交叉变异函数面板一片红

现象：Geostatistical Wizard 里交叉变异函数点云散乱，拟合线怎么调都对不上。原因：主辅变量的采样位置完全不重叠，交叉变异函数需要主辅变量在相近位置都有采样才能算，位置错开太远就算不出有效值。解决：检查两套点的空间重叠度，如果重叠率低于 30%，协同克里金不适用；考虑先把协变量插值到主变量位置，再评估相关性。

4.4 计算时间爆炸

现象：主变量几千个点、协变量几万个点，跑一次协同克里金要几个小时。原因：协同克里金的方程组规模是主变量邻居数加协变量邻居数的三次方，协变量邻居数一上去，计算量指数增长。解决：把协变量邻居数压到 10 以内；或者先用协变量做一次普通克里金，把结果栅格作为"趋势面"，再用主变量做简单克里金残差插值——这是协同克里金的近似替代，速度快一个量级，精度损失通常可接受。

4.5 输出栅格在边界外出现异常值

现象：预测栅格在研究区边界外一圈出现明显偏离的条带。原因：搜索邻域在边界处只覆盖到单侧的点，权重分布不对称，外推失控。解决：在环境设置里把输出范围裁到研究区边界内；或者用掩膜（Mask）把边界外区域设为 NoData。这个坑在普通克里金里也有，但协同克里金因为多了协变量维度，边界效应更明显。

5. 进阶：用交叉验证挑模型，以及一个我常备的替代方案

协同克里金跑完之后，别急着出图，先做交叉验证。ArcGIS 的 Geostatistical Wizard 最后一步会给交叉验证统计量，重点看三个：均方根误差（RMSE）、平均标准误差（MSE）、标准化均方根误差（RMSSE）。RMSSE 最理想是接近 1，大于 1 说明预测方差被低估（模型过于自信），小于 1 说明预测方差被高估（模型过于保守）。我一般会跑三组参数：一组自动拟合、一组手动调块金、一组换变异函数模型，比 RMSSE 谁更接近 1，再比 RMSE 谁更小。

如果协同克里金怎么调都不理想，我常备的替代方案是"回归克里金"（Regression Kriging）：先用主变量对协变量做普通最小二乘回归或广义线性回归，拿到回归残差，再对残差做普通克里金，最后把回归预测面和残差插值面相加。这个方案在 ArcGIS 里没有一键工具，但用arcpy.ga.OrdinaryKriging加arcpy.analysis里的回归工具能拼出来。它的好处是：回归部分处理协变量的全局趋势，克里金部分处理局部空间自相关，两者解耦，调参比协同克里金直观得多。代价是它不像协同克里金那样在数学上统一求解权重，理论上不如协同克里金优雅，但实操中往往更稳。

import arcpy from arcpy.ga import OrdinaryKriging # 第一步：用协变量对主变量做回归，输出残差字段 # 这里用 Ordinary Least Squares，实际可用 GWR 提升局部拟合 arcpy.stats.OrdinaryLeastSquares( in_features="soil_samples", dependent_variable="om_gkg", explanatory_variables=[["ndvi", "LINEAR"]], output_features="soil_with_resid" ) # 第二步：对回归残差做普通克里金 ok = OrdinaryKriging("soil_with_resid", "resid_field") ok.modelType = "SPHERICAL" ok.lagSize = 500 ok.searchNeighborhoodNumberOfPoints = 12 ok.execute("resid_pred", "resid_err") # 第三步：回归预测面 + 残差插值面 = 最终预测面 # 回归预测面由 OLS 工具直接输出，用栅格计算器相加 arcpy.sa.Plus("ols_pred", "resid_pred").save("final_pred")

这段脚本的关键在第三步：回归预测面是全局趋势，残差插值面是局部修正，两者相加才是最终结果。参数上，回归部分用 OLS 还是 GWR（地理加权回归）取决于趋势是否空间平稳——如果回归系数在研究区内明显变化，用 GWR；否则 OLS 够用。残差克里金的参数调法和普通克里金一样，但注意残差的变异函数通常比原始变量更"碎"，块金值占比更高，这是正常的，因为回归已经抽走了大尺度趋势。

最后说个习惯：我每次跑协同克里金之前，一定先画一张主变量和协变量的散点图，再画一张两套点的空间分布叠加图。散点图看相关性是否线性、有没有离群点；叠加图看两套点的空间重叠度。这两张图花不了十分钟，但能挡掉八成"跑完才发现数据不适合"的后悔药。协同克里金不是万能插值器，它是一个有明确适用条件的工具，条件不满足时，退回普通克里金或者回归克里金，比硬调参数明智得多。希望帮到你。

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