企业官网建设流程全解析

简介：面向空间统计与地学环境领域的数据处理需求，这是一套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结构体都围绕这三个字段展开。

三个参数之间的关系并不独立。变程决定了“多远的点能参与预测”，如果变程只有平均采样间距的一半，那么在最小间距上就已经没有相关性，插值结果退化为局部均值；如果变程设得超过数据覆盖范围，每个待估点都会被所有样本点拉住，权重矩阵会变得非常稠密，数值上更容易病态。基台则和样本方差对应，如果拟合出来的基台只有样本方差的三分之一，就该怀疑数据是否去趋势过头，或者模型选错。

2.3 三种理论模型怎么选：球状的经典，指数的顽强，高斯的陷阱

拟合模型的终极目标只有一个：让协方差矩阵正定且接近真实空间结构。球状模型写作 γ(h) = C0 + C1 * (1.5 * h/a - 0.5 * (h/a)³)，在 h < a 时使用，之后恒为 C0+C1。它是有变程模型里的“默认选项”，适合大多数土壤、气象、地质变量，因为它在变程处平滑地到达基台。指数模型 γ(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

注意代码里特意处理了 i==j 的情况：点在自身位置时距离为零，半方差理论值应为 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₁(x+h) - Z₁(x)) · (Z₂(x+h) - 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(n+1: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%，就先回头调块金和变程的初始值再重跑，而不是急着出图。希望帮到你。

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