简介面向空间统计与地学环境领域的数据处理需求这是一套MATLAB实现的克里金与协同克里金插值代码。资源覆盖半方差函数计算与模型拟合、克里金权重求解、多变量协同预测等核心环节可帮助初学者及需要快速上手的研发人员把理论方法转化为可运行的程序。压缩包共8个文件含7个.m脚本和1个.mat测试数据文件整包约44KB体量轻巧。脚本功能划分清晰普通克里金与协同克里金各自独立实现半方差函数拟合、参数寻优和测试模块均单独封装配合自带测试脚本与示例数据能直接运行观察从变异函数分析到插值输出的完整流程。目前已有2680人学习可见该主题在空间插值实践中的需求较高。通过学习这套代码读者可掌握半方差函数建模、协方差矩阵构建、辅助变量引入等关键操作并基于示例数据验证插值效果之后还能替换实际数据或调整模型参数用于科研实验与工程项目中的空间预测与分析。1. 克里金插值为什么它追求“方差最小”而不是“误差最小”克里金插值最反直觉的地方在于它并不追求把误差做到最小——因为它根本不知道真实值在哪。它做的是一件看似绕路的事在满足无偏约束的前提下把估计方差做小然后在每个待估点同时输出一个预测值和一个方差值。前者用来画等值线后者用来回答“这条线到底可信多少”。这也是克里金插值和反距离加权法的本质差别反距离加权只是拍脑袋给一个距离衰减指数而克里金从数据本身拟合空间相关性。后续所有问题——变程怎么取、块金怎么调、协同克里金怎么引入辅助变量——都从“最小化估计方差”这个出发点派生出来搞懂它才算真正看懂了这套方法的数学骨架。2. 变异函数半方差到底在拟合什么平稳性假设为什么是本命克里金插值从数据中提取空间相关性的工具是变异函数variogram。一组样点画出来横轴是采样点对之间的距离 h纵轴是半方差 γ(h)。这个散点云本身不能直接参与方程求解必须先拟合出一个理论模型。为什么需要这一道绕弯因为克里金方程组需要的不是一个样本点到另一个样本点的实际距离而是一个稳定的协方差函数能把任意距离映射成“统计意义上的相似度”。经验半方差离散、抖动大直接拿去解方程组会得到负权重甚至奇异矩阵先拟合模型就是为了把一个黑匣子变成可解释的参数——这也是为什么我们把第 2 章放在所有源码讲解之前。2.1 经验半方差的计算点云里藏着什么信息对于采样点 x_i 和 x_j距离记为 h半方差的定义是 γ(h) 1/2 * E[(Z(x_i) - Z(x_j))²]。实际操作中我们用所有距离接近 h 的样本对求平均这就是经验半方差。公式里的 1/2 不是随手写的在二阶平稳假设下这个值恰好等于 C(0) - C(h)它把“变异”和“相关”统一了起来。function [gamma_mean, lag_mean, N_pairs] exp_variogram(xy, z, maxLag, nLags) % xy: 采样点坐标N×2 % z: 采样点上的值N×1 % maxLag: 计算半方差的最大距离一般取最大采样间距的一半左右 % nLags: 把距离分成多少段常用 8~12 d2 pdist2(xy, xy); % 两两距离矩阵 dif (z - z) .^ 2; % 所有点对的平方差 g_all 0.5 * dif; % 每个点对的半方差贡献 lagEdges linspace(0, maxLag, nLags 1); gamma_mean zeros(1, nLags); lag_mean zeros(1, nLags); N_pairs zeros(1, nLags); for k 1:nLags idx find(d2 lagEdges(k) d2 lagEdges(k 1)); if isempty(idx) continue; end gamma_mean(k) mean(g_all(idx)); % 该距离区间内的平均半方差 lag_mean(k) mean(d2(idx)); % 区间内实际平均距离 N_pairs(k) length(idx); % 该区间的点对数量反映可靠性 end end这段代码里有一个细节值得停下来看pdist2 生成的距离矩阵是对称的所以每个样本对会被统计两遍。均值不受影响但 N_pairs 会翻倍。真正生产级的实现会用triu或 mask 掉对角线及下三角只统计每个点对一次否则点对数量偏大你在判断“尾部点云是否可靠”时会被误导。maxLag 的选择是第一个容易翻车的地方取得过大远端 lag 内点对太少点云尾部抖动剧烈取得过小长距离结构被截掉模型变程根本拟合不出来。通用的经验是取采样点最大间距的一半再检查每个区间 N_pairs 是否都大于 30。2.2 块金、基台、变程三个参数和空间结构的对应关系经验变异函数只是中间产物最终拿来解克里金方程的是一个理论模型最常用的是这三类参数块金nugget、基台sill、变程range。它们分别描述空间结构的三个层次几乎所有代码里的model结构体都围绕这三个字段展开。参数数学含义物理场景代码里最常见的坑块金 C0h 趋近 0 时 γ(h) 的取值测量误差 微尺度随机变异初始值给 0导致矩阵在原点处过于刚硬基台 C0C1γ(h) 在变程外稳定到的上限空间变异的总方差基台与样本方差偏离过大拟合方向基本错了变程 a空间相关性消失的距离采样方案的有效作用半径变程大于最大采样点对距离相当于全平面相关三个参数之间的关系并不独立。变程决定了“多远的点能参与预测”如果变程只有平均采样间距的一半那么在最小间距上就已经没有相关性插值结果退化为局部均值如果变程设得超过数据覆盖范围每个待估点都会被所有样本点拉住权重矩阵会变得非常稠密数值上更容易病态。基台则和样本方差对应如果拟合出来的基台只有样本方差的三分之一就该怀疑数据是否去趋势过头或者模型选错。2.3 三种理论模型怎么选球状的经典指数的顽强高斯的陷阱拟合模型的终极目标只有一个让协方差矩阵正定且接近真实空间结构。球状模型写作 γ(h) C0 C1 * (1.5 * h/a - 0.5 * (h/a)³)在 h a 时使用之后恒为 C0C1。它是有变程模型里的“默认选项”适合大多数土壤、气象、地质变量因为它在变程处平滑地到达基台。指数模型 γ(h) C0 C1 * (1 - exp(-h/a)) 在原点附近线性上升收敛更慢实际变程大约在 3a 的位置适合空间变化比较剧烈的变量。高斯模型 γ(h) C0 C1 * (1 - exp(-(h/a)²)) 在原点附近几乎为零曲线非常缓适合极端平滑的变量但也是三种模型里最危险的它在短距离内对微小距离变化不敏感组装出来的协方差矩阵行与行之间高度相似条件数偏高。我一般会这么选先看经验变异函数在原点附近的行为。如果第一个 lag 的半方差显著高于 0优先考虑带块金的球状模型如果点云在原点附近接近线性上升指数模型更稳除非你明确知道变量非常光滑否则不要一上来就用高斯。拟合时还有个血泪经验实测数据的经验变异函数在长距离处往往呈波浪状这是短距离结构之外还存在更大尺度结构的信号不要为了把尾部拟合漂亮而强行加第三个结构协同克里金阶段那个更大的尺度结构会由协变量去承接。3. 克里金方程组的源码级拆解λ权重是怎么解出来的克里金的数学本质是一个约束优化问题在权重之和为 1 的条件下最小化预测值的估计方差。用拉格朗日乘子法转成一个线性方程组后问题就变成组装矩阵、解线性方程。多数人在这一步不是被公式难住而是被“该用协方差还是半方差”绕晕。两种写法在数学上完全等价但代码里的细节完全不同这一章直接进入一个可运行的普通克里金函数。3.1 从协方差矩阵到增广方程组一个完整的普通克里金函数function [z_hat, var_hat, lambda] ordinary_kriging(xy, z, xy_target, model) % xy: 采样点坐标N×2 % z: 采样点上的实测值N×1 % xy_target: 单个待估点坐标1×2 % model: 结构体包含 .type 模型类型、.nugget 块金、.sill 基台、.range 变程 n size(xy, 1); A zeros(n 1, n 1); for i 1:n for j 1:n d norm(xy(i, :) - xy(j, :)); if i j A(i, j) model.sill; % 原点处协方差等于基台不是基台减块金 else A(i, j) model.sill - var_model(d, model); % C(h) sill - gamma(h) end end end A(1:n, n 1) 1; % 无偏约束放在最后一列 A(n 1, 1:n) 1; % 无偏约束放在最后一行 A(n 1, n 1) 0; % 增广位置固定为 0 b zeros(n 1, 1); for i 1:n d norm(xy(i, :) - xy_target); b(i) model.sill - var_model(d, model); % 待估点与样本点的协方差 end b(n 1) 1; % 权重之和为 1 的约束值 w A \ b; lambda w(1:n); % 各样本点权重 mu w(n 1); % 拉格朗日乘子 z_hat lambda * z; % 预测值 % 估计方差: sigma^2 C(0) - lambda * C_d mu var_hat model.sill - lambda * b(1:n) mu; if var_hat 0 var_hat 0; % 矩阵接近奇异时浮点误差可能带来负值先截断再看条件数 end end注意代码里特意处理了 ij 的情况点在自身位置时距离为零半方差理论值应为 0协方差应为 C(0)sill。如果不加这个判断直接调用var_model(0)大多数实现会返回块金值于是对角元变成sill - nugget。块金一旦偏大对角优势被削弱矩阵条件数急剧恶化最后解出来的权重就会乱跳。这是初版克里金代码里最常见的翻车点症状是预测结果看不出空间趋势反而像噪声放大。3.2 拉格朗日乘子无偏约束如何被“惩罚”进最小方差如果不加任何约束最小二乘意义下的最优权重之和通常不等于 1预测结果会有系统性偏差。克里金的做法是把约束塞进目标函数在 Var(ẑ - z) 上加上一项 2μ(Σλ - 1)然后对每个 λ_i 求偏导。这一项里的 μ 就是增广矩阵最后一行的解。它是有实际含义的μ 度量了这个无偏约束给方差带来的“代价”。这也解释了为什么克里金的权重里会出现负值。当采样点在某个区域过于聚集时为了满足权重和为 1矩阵会给冗余的点分配负权重来抵消过量信息这就是克里金里的屏蔽效应。出现负权重不一定是算法出错但负权重绝对值超过几倍正常权重时通常意味着矩阵条件数很高数据本身存在严重的空间共线性。在代码里请养成一个习惯解出 λ 后先打印cond(A)和sum(lambda)前者告诉你数值稳定性后者告诉你约束是否被精确满足。如果条件数超过 1e12先回头调整块金或变程而不是继续追着权重看。3.3 估计方差怎么算以及它为什么不等于“真实误差”预测点处的估计方差公式是 σ² C(0) - Σλ_i·C(x_i, x_0) μ。在源代码里lambda * b(1:n)是权重向量和样本点-待估点协方差向量的内积μ 因为有符号位置的问题务必和推导公式保持一致。这个方差值只反映模型内部的方差结构不包含“模型选错”造成的偏差。换句话说当你的变异函数拟合偏乐观时σ² 会系统性偏小置信区间漂亮得不真实。这是克里金最迷人的地方也是它最需要交叉验证兜底的原因。4. 协同克里金Co-Kriging辅助变量是怎么被“加权”进主变量的协同克里金跟普通克里金的差别不是多插值一个变量那么简单。它在同一个方程组里同时使用主变量和协变量的采样数据核心是交叉半方差。辅助数据密度更高时协同克里金能显著缩小主变量的估计方差这就是“Co-Kriging”这个名字背后真正的价值。4.1 交叉半方差衡量两个变量之间的空间互相关交叉半方差的定义是 γ₁₂(h) 0.5 * E[(Z₁(xh) - Z₁(x)) · (Z₂(xh) - Z₂(x))]和普通半方差唯一的区别是少了平方。没有平方意味着它天然可以为负表示两个变量在某个距离尺度上负相关。这个性质决定了交叉半方差不能像普通半方差那样随便套一个正定模型。function gamma12 cross_variogram(xy, z1, z2, maxLag, nLags) % 假设主变量和协变量采样在同一组坐标上 % z1: 主变量采样值 % z2: 协变量采样值 d2 pdist2(xy, xy); % 距离矩阵 g12 0.5 * (z1 - z1) .* (z2 - z2); % 交叉半方差的核心乘积差不平方 % 按 lag 分箱求平均过程与 exp_variogram 完全相同 lagEdges linspace(0, maxLag, nLags 1); for k 1:nLags idx find(d2 lagEdges(k) d2 lagEdges(k 1)); gamma12(k) mean(g12(idx)); end end实际协作中最大的坑在于主变量和协变量很少在同一组坐标上采样这时需要用两组坐标分别计算距离矩阵交叉项变成 n₁×n₂ 的矩形块。代码实现时还要保证 C₁₂ 和 C₂₁ 互为转置很多初版代码只组装了一半矩阵不对称求解结果直接发散。拟合时交叉半方差模型不能独立拟合它必须与 C₁₁、C₂₂ 一起满足联合正定条件也就是 C₁₂(h)² ≤ C₁₁(0)·C₂₂(0)。最常见的思路是给两个自变异函数和一个交叉变异函数选择同一类型模型、同一变程用线性模型近似权重由相关系数控制。4.2 协同克里金的联合方程组为什么协变量权重之和必须为 0协同克里金的增广矩阵比普通克里金复杂一个维度。假设主变量和协变量各有 n 个采样点权重向量分成两段主变量权重 λ₁ 和协变量权重 λ₂。矩阵结构是分块的。% 分块矩阵结构示意假设主变量与协变量采样点相同n1 n2 n % A 矩阵维度为 (2n 2) × (2n 2) A [C11, C12, ones(n,1), zeros(n,1); C21, C22, zeros(n,1), ones(n,1); ones(1,n), zeros(1,n), 0, 0; zeros(1,n), ones(1,n), 0, 0]; b [C1_target; C2_target; 1; 0]; w A \ b; lambda1 w(1:n); % 主变量权重 lambda2 w(n1:2*n); % 协变量权重注意最后两个约束一个的值是 1一个是 0。主变量权重之和为 1 保持无偏性这个好理解协变量权重之和为 0 才是协同克里金的辨识条件。为什么不是 1因为协变量在这里充当的是“辅助修正项”它的均值本身被吸收了。如果允许协变量权重和为任意值协变量的整体均值就会混入主变量估计而主变量本身还有一个待估计的未知常数均值两个参数无法被同时识别。把这个约束理解透你就明白协同克里金并不是简单地把两个克里金结果叠加它是在一个统一的模型里做联合推断。这个约束也决定了协同克里金的一个适用前提协变量与主变量必须有真实的交叉相关性。如果 ρ₁₂ 接近 0解出来的 λ₂ 会普遍接近 0增广矩阵却白白增加了维度数值开销和噪声都上升。跑 Co-Kriging 之前先看一眼交叉散点图没有相关性的数据硬上协同克里金结果通常还不如普通克里金。4.3 协同克里金 vs “两步走”为什么不能先插值协变量再回归很多人的朴素做法是先把协变量用普通克里金插值到整个网格然后把主变量对这个插值面做回归或者对回归残差再插值。这种两步法的问题在于第一步的插值误差被当成真实协变量值传给了第二步。协变量在稀疏区域会带着明显的插值方差这个信息在两步法里完全丢失最终的主变量预测面会过度自信。协同克里金把协变量的采样数据直接放进方程组让方差结构同时沿两条路径传播一步得到 BLUP。协变量数据的密度越高从这条额外路径获得的信息就越多主变量的预测方差降得越明显。协同克里金在代码实现上的代价是求解规模从 O(n³) 涨到 O((n₁n₂)³)。如果协变量样本有几万个直接解全矩阵会慢到不可接受。实际可行的做法是用搜索邻域限定参与计算的样本量比如只取距离待估点最近的 30~50 个主变量样本和 50~80 个协变量样本。邻域半径再结合变程来定超过变程 1.5 倍的点对预测贡献很小还容易引入尾部噪声。5. 克里金与协同克里金实战避坑从去趋势到神秘负权重这一章值得每一个拿这套代码跑出“好看”结果的人反复看。以下五类问题我都在不同数据集上见过现象、原因和对应的处理方式都按真实排查路径写出来。5.1 数据有趋势经验变异函数出现假周期性现象经验变异函数在长距离处不趋于基台反而持续上升甚至出现波浪式起伏拟合出来的变程远大于数据覆盖范围。 原因区域本身存在系统性趋势比如高程从南到北递增。数据不满足二阶平稳半方差里混入了趋势项呈现出伪周期或伪相关。 解决对主变量做趋势面回归或至少做均值中心化用残差去拟合变异函数。残差做完克里金预测后再把趋势面加回来。这个策略在代码里只多了两行一个polyfit一个polyval但效果天差地别。5.2 坐标单位不统一变程形同虚设现象变程、基台拟合结果异常插值面沿某一方向被拉伸或压缩预测结果像被揉皱的纸。 原因x 坐标用的是米y 坐标用的是公里或者一部分数据经纬度、一部分数据投影坐标。克里金对欧氏距离没有任何“单位感知”它只认数值。 解决把全部坐标统一投影到同一坐标系并做一次缩放让两个轴的数值范围相近。检查方法很简单计算所有点对的极差如果 x 方向的极差是 y 方向的 10 倍以上先归一化再跑拟合。5.3 主变量与协变量数量级差异过大交叉半方差直接崩溃现象协同克里金解出的权重里协变量权重普遍极小或极大预测方差出现负值。 原因主变量浓度是 0~100 mg/kg协变量高程是 200~1000 m交叉半方差里差值的乘积高达百万量级矩阵不同分块之间尺度失衡求解器在病态矩阵上挣扎。 解决建模前把主变量和协变量都做标准化各自减去均值除以标准差。预测完成后把结果乘回原尺度。这一步对普通克里金不是必须的但对协同克里金几乎是强制要求。5.4 矩阵条件数过大出现神秘负权重现象一部分样本点的权重为负绝对值超过其他权重数值几倍预测值超出样本数据的物理范围。 原因采样点在某个区域过度聚集加上块金设置过小、变程设置过大协方差矩阵接近奇异。屏蔽效应在极端情况下会把聚类内部的点变成强负权重去抵消冗余信息。 解决打印条件数cond(A)定位问题。先增大块金把矩阵从奇异边缘拉回来如果还不行对采样点做空间去冗余比如最小距离内只保留一个点。克里金权重里出现小负值是正常的大面积强负权重则是矩阵病态而不是算法错误。5.5 模型拟合贴着经验点云交叉验证误差却惨不忍睹现象拟合曲线肉眼看起来完美贴合但留一交叉验证 RMSE 比预测方差大一倍多。 原因初始参数给得不好。最常见的是块金初始值设为 0拟合算法为了在高拟合优度方向收敛把基台压到样本方差附近模型在原点处过度光滑等于默认“短距离内的样点完全相关”。这会让克里金对采样点附近的山峰过度信任预测方差严重低估。 解决强制块金的初始值不小于经验变异函数第一个 lag 的值并检查基台是否落在样本方差的 0.8~1.2 倍区间内。拟合完成不是终点必须用交叉验证结果做最终判定。6. 落进源码之前先问一句这个模型的“最优”会不会骗人克里金的最优性有一个隐含前提变异函数模型是对真实空间结构的正确描述。模型偏差不在“估计方差”里这也是为什么一套代码跑完图面漂亮不等于结果可靠。我每次拿到新的克里金结果都会先用留一交叉验证审视一遍而不是直接看对结果负责的路径是否一致。% 留一交叉验证每次抽掉一个样本点用其余点预测被删点 err zeros(n, 1); pred_var zeros(n, 1); for i 1:n idx setdiff(1:n, i); [pred(i), pred_var(i)] ordinary_kriging(xy(idx, :), z(idx), xy(i, :), model); err(i) pred(i) - z(i); end rmse sqrt(mean(err.^2)); avg_sigma sqrt(mean(pred_var)); ratio rmse / avg_sigma; % 经验判断1.2 以上过平滑0.8 以下过激进这个比值比任何单点的 R² 都有意义。它揭示的是克里金设计中的根本取舍块金既是测量误差和微尺度变异的描述也起到了正则化的作用。块金设得大矩阵更稳定但预测面对细节的响应变得迟钝结果退化成平滑均值块金设得小预测面波峰波谷分明方差看起来很“自信”但交叉验证会立刻拆穿它。参数搜索的方向不是“拟合优度最高”而是“交叉验证 RMSE 和平均方差尽量接近”。从那以后我每次跑完克里金都会强制走一遍 LOOCV比较这两个值如果差距超过 20%就先回头调块金和变程的初始值再重跑而不是急着出图。希望帮到你。本文还有配套的精品资源点击获取