企业官网建设流程全解析

简介：一套面向地理空间数据插值的MATLAB实现，完整覆盖普通克里金（Kriging）与协同克里金（Co-Kriging）两种方法，适合地质统计、环境科学、气象遥感等领域的研究者，也可作为相关课程教学与毕业设计的参考代码。压缩包共8个文件，以7个MATLAB脚本（.m）为主，涵盖普通克里金、协同克里金、变差函数计算与拟合、带边界参数优化等核心模块，并配有1个.mat格式测试数据，整体仅44KB，便于快速下载和本地运行。代码结构清晰，通过测试脚本即可观察插值结果，也可自行替换数据并调整变差函数模型、搜索半径等参数，深入理解空间相关性建模与权重求解的完整流程。已有2677人浏览学习，说明其实用性受到一定认可。借助这套代码，读者既能掌握两类克里金方法在MATLAB中的落地实现，也能为后续开展更为复杂的空间插值研究打下基础。 做空间插值的人，应该都绕不开这两个词：克里金和协同克里金。前几年我做土壤属性空间制图，用反距离加权插出来一片“马赛克”，被审稿人追着问为什么不用克里金。后来老老实实把变异函数和克里金方程组啃了一遍，又用MATLAB把普通克里金和协同克里金的代码一行行写出来，才真正明白这套方法到底强在哪。这篇东西就是那段时间的代码和踩坑记录，给需要在项目里或论文里用MATLAB做克里金及协同克里金插值的工程师和研究生参考，不用额外工具箱，主程序和函数一起用，可以直接跑通。

1. 克里金不是玄学：它到底在优化什么

1.1 从“脑子一热给权重”到“让数据自己说话”

反距离加权（IDW）之所以让人不放心，是因为它的权重全靠人拍脑袋定：距离平方反比、距离三次方反比，凭什么指数是2而不是2.7？点多了会不会互相打架？这些在IDW里都没有答案。克里金的核心变化是——权重不是手动指定的，而是通过空间相关性结构算出来的。它要求你先回答一个问题：距离多远的时候，两个点基本上就不相关了？这个“距离多远才不相关”的规律，就是变异函数（variogram）干的事。

克里金本质上是一个线性加权估计：

Ẑ(x₀) = ΣλᵢZ(xᵢ)

但权重λᵢ不再是人定的，而是满足两个条件：无偏性（权重之和等于1）和估计方差最小。也就是说，克里金在数学上明确定义了“最优”——在所有无偏线性估计里，它的均方误差最小。这就是“最优线性无偏估计”（BLUE）这个说法的来源。

1.2 协同克里金：主变量不够，请辅助变量来帮忙

实际项目中经常遇到一种尴尬：目标变量采样点特别少，比如土壤重金属铜，打钻采样费用高，全区可能只有二三十个点；但另一个变量采样密度很高，比如地形高程、土壤有机质，几千上万个点都有数据。铜和有机质如果确实存在空间相关性，为什么不把这些高密度的辅助信息用起来？

协同克里金（Cokriging）要解决的就是这个问题。它不只对主变量Z做加权，还把辅助变量V也加进来：

Ẑ(x₀) = ΣλᵢZ(xᵢ) + ΣνⱼV(xⱼ)

无偏条件变成：主变量的权重之和为1，辅助变量的权重之和为0。辅助变量的作用不是直接“填”缺测点，而是通过它和主变量的空间交叉相关性，修正主变量相邻点之间的估计权重。这个思路在气温插值里特别常见：气温站点稀疏，但DEM（数字高程模型）有完整连续的面，把高程作为辅助变量做协同克里金，山区插值结果通常比只用站点做普通克里金平滑得多。

2. 普通克里金MATLAB实现：先拟合变异函数，再解方程组

2.1 第一步：用实测数据计算经验变异函数

变异函数是克里金的灵魂，没拟合好，后面方程组解得再漂亮都没意义。现实采集的样点不是规则的，我们先把所有点对找出来，按距离分箱（lag），每个箱子里计算半方差：

γ(h) = (1/2N(h)) · Σ [Z(xᵢ) - Z(xⱼ)]²

这个函数就是经验变异函数，也是后续拟合的基础。MATLAB实现就按这个公式直接写：

function [lag, gamma_emp] = empirical_variogram(coords, data, n_lags) % 计算经验变异函数 % coords: 样点坐标，n x 2 % data: 样点观测值，n x 1 % n_lags: 距离分箱数量，默认20 if nargin < 3, n_lags = 20; end n = size(coords, 1); pairs = zeros(n*(n-1)/2, 2); k = 0; for i = 1:n-1 for j = i+1:n h = norm(coords(i,:) - coords(j,:)); g = 0.5 * (data(i) - data(j))^2; k = k + 1; pairs(k, :) = [h, g]; end end hmax = max(pairs(:,1)); lag_step = hmax / n_lags; lag = zeros(n_lags, 1); gamma_emp = zeros(n_lags, 1); for i = 1:n_lags lo = (i-1) * lag_step; hi = i * lag_step; idx = pairs(:,1) >= lo & pairs(:,1) < hi; lag(i) = (i - 0.5) * lag_step; if sum(idx) > 0 gamma_emp(i) = mean(pairs(idx, 2)); else gamma_emp(i) = nan; end end end

这里有个技巧：n_lags不要设太多，模型点数多了反而没有足够点对支撑。我一般取15到25个箱子，每个箱子里至少有几十个点对，否则拟合出来全是跳动的噪声。判断分箱是否合理，就看后面画的散点图是不是有一个“先上升后趋于平稳”的形态。

2.2 第二步：用fminsearch拟合变异函数模型，不依赖优化工具箱

经验变异函数只是一堆散点，克里金方程需要连续函数。常用模型有三个：球状模型、指数模型、高斯模型。它们表达的空间相关性曲线略有区别：

• 球状模型：在距离小于变程a时上升，超过a后稳定在基台值，线性到曲线过渡，最常用。
• 指数模型：渐进逼近基台值，比球状模型平滑。
• 高斯模型：在原点附近非常平缓，适合连续性特别强的变量。

模型的一般形式都是块金值（nugget）+ 偏基台值（partial sill）× 形状函数。写成MATLAB函数：

function g = variogram_value(h, p, model) % 变异函数理论模型 % p = [c0, sill, a]，分别是块金、基台、变程 c0 = p(1); sill = p(2); a = p(3); if a == 0, a = eps; end switch lower(model) case 'spherical' g = zeros(size(h)); idx = h <= a; g(idx) = c0 + (sill - c0) * (1.5*(h(idx)/a) - 0.5*(h(idx)/a).^3); g(h > a) = sill; case 'exponential' g = c0 + (sill - c0) * (1 - exp(-3*h/a)); case 'gaussian' g = c0 + (sill - c0) * (1 - exp(-3*(h/a).^2)); otherwise error('未知变异函数模型'); end end

拟合的过程就是让理论曲线尽量贴近经验散点，我这里用fminsearch做最小二乘拟合，这个函数在MATLAB基础版里就有，不需要买优化工具箱：

function p = fit_variogram(lag, gamma_emp, model) % 最小二乘拟合变异函数参数 valid = ~isnan(gamma_emp) & gamma_emp >= 0; lag = lag(valid); gamma_emp = gamma_emp(valid); p0 = [min(gamma_emp), max(gamma_emp), 0.6*max(lag)]; if p0(2) <= p0(1), p0(2) = p0(1) + 1; end obj = @(p) sum((variogram_value(lag, p, model) - gamma_emp).^2); options = optimset('Display', 'off', 'TolX', 1e-6, 'TolFun', 1e-6, 'MaxIter', 1000); p = fminsearch(obj, p0, options); p(1) = max(p(1), 0); % 块金不能为负 p(2) = max(p(2), p(1) + 1e-6); % 基台必须大于块金 p(3) = max(p(3), eps); % 变程必须为正 end

p0的初始化并不是随便拍的：块金从最小值开始，是因为理论上距离为0时变异函数接近0；基台取经验变异函数的最大值，反映空间方差的上限；变程取最大距离的0.6倍，保证初始曲线不会太离谱。fminsearch容易陷入局部最优，所以初值要给得“有道理”，而不是从0开始。

2.3 第三步：构建普通克里金方程组，完成空间插值

有了变异函数模型，就可以组方程组。普通克里金对每个待估点x₀求解：

Γ · w = b

其中Γ是样点之间的变异函数矩阵，最后加一行1和一列1作为无偏性约束，b是待估点与各样点之间的变异函数值，加一个1对应拉格朗日乘子。解出来的权重w后n个值就是各样点的克里金权重，最后一个值是拉格朗日乘子，用来算克里金方差。MATLAB写出来非常直接：

function [Z_pred, Kvar] = ordinary_kriging(pred_coords, coords, data, vg_model, params) % 普通克里金插值 n = length(data); A = zeros(n+1, n+1); for i = 1:n for j = i:n h = norm(coords(i,:) - coords(j,:)); v = variogram_value(h, params, vg_model); A(i,j) = v; A(j,i) = v; end end A(n+1, 1:n) = 1; A(1:n, n+1) = 1; A(n+1, n+1) = 0; np = size(pred_coords, 1); Z_pred = zeros(np, 1); Kvar = zeros(np, 1); for k = 1:np b = zeros(n+1, 1); for i = 1:n b(i) = variogram_value(norm(coords(i,:) - pred_coords(k,:)), params, vg_model); end b(n+1) = 1; w = A \ b; Z_pred(k) = w(1:n)' * data; Kvar(k) = max(sum(w(1:n) .* b(1:n)) + w(n+1), 0); end end

注意这里可能有人会问：不是说克里金用协方差矩阵吗，为什么这代码里全是变异函数？两者在无偏性约束下是等价的，因为C(h) = 基台值 - γ(h)，而无偏性要求权重和等于1，常数基台值会自然抵消。直接组变异函数矩阵可以省掉一次转换，代码也更好理解。

克里金方差出现负数的情况，实测里会碰到，主要原因是变异函数模型参数不合理或者矩阵病态。我在这里用max(…,0)兜底，但更重要的是后面第4章里讲怎么从根源上避免这个问题。

3. 协同克里金MATLAB实现：把辅助变量写进方程组

3.1 协同克里金的分块矩阵结构

协同克里金的本质就是把普通克里金的单变量系统扩展成多变量系统。假设有主变量Z（nz个采样点）和辅助变量V（nv个采样点），待估点x₀的估计涉及四个块的变异/交叉变异信息：

方程组右侧也有对应的三个部分：待估点到主变量样点的变异函数值、待估点到辅助变量样点的交叉变异函数值，以及两个无偏性约束条件（1和0）。整个系统是一个(nz+nv+2) × (nz+nv+2)的分块矩阵，组起来不难，但索引别搞错。

3.2 交叉变异函数：工程上最省事的处理方式

交叉变异函数是协同克里金里真正麻烦的地方。严格做法需要对Z和V的所有交叉点对统计互变差函数γzv(h)，再单独拟合。但这里有个坑：如果Z和V不是在同一个站点观测的，直接统计会出现大量错位配对，拟合出来的交叉变异函数可能不满足正定性。

我在实际项目里常用的折中方案是内嵌相关模型（intrinsic correlation model）的思想：假设主辅变量的交叉相关性由一个全局相关系数ρ控制，交叉变异函数近似表达为：

γzv(h) ≈ ρ · √(γzz(h) · γvv(h))

其中ρ就是Z和V在采样点上的皮尔逊相关系数。这个近似的好处是：只要γzz和γvv都是有效的变异函数模型，γzv就能自动保持正定性，不用单独费劲拟合交叉变异函数。缺点是假设了Z和V在整个研究区有稳定的相关性，如果两者相关性随空间位置变化特别大，这个近似就会失真。所以做之前务必先算一下Z和V的相关系数，低于0.5就别硬上协同克里金，不如把辅助变量做回归残差克里金。

3.3 协同克里金预测代码

按上面思路，核心函数这样写：

function [Z_pred, Kvar] = cokriging_predict(pred_coords, Zcoords, Zdata, Vcoords, Vdata, vg_model, pzz, pvv, rho) % 协同克里金插值 % Zcoords/Zdata: 主变量采样点和观测值 % Vcoords/Vdata: 辅助变量采样点和观测值 % pzz/pvv: 主变量和辅助变量各自拟合出的变异函数参数 [c0, sill, a] nz = length(Zdata); nv = length(Vdata); nd = nz + nv; A = zeros(nd + 2, nd + 2); % 主变量自身变异函数块 for i = 1:nz for j = 1:nz A(i,j) = variogram_value(norm(Zcoords(i,:) - Zcoords(j,:)), pzz, vg_model); end end % 辅助变量自身变异函数块 for i = 1:nv for j = 1:nv A(nz+i, nz+j) = variogram_value(norm(Vcoords(i,:) - Vcoords(j,:)), pvv, vg_model); end end % 交叉变异函数块：rho * sqrt(gzz*gvv) for i = 1:nz for j = 1:nv h = norm(Zcoords(i,:) - Vcoords(j,:)); gzz = variogram_value(h, pzz, vg_model); gvv = variogram_value(h, pvv, vg_model); A(i, nz+j) = rho * sqrt(gzz * gvv); A(nz+j, i) = A(i, nz+j); end end % 无偏性约束：Σλ=1, Σν=0 A(1:nz, nd+1) = 1; A(nz+1:nd, nd+1) = 0; A(nd+1, 1:nz) = 1; A(nd+1, nz+1:nd) = 0; A(1:nz, nd+2) = 0; A(nz+1:nd, nd+2) = 1; A(nd+2, 1:nz) = 0; A(nd+2, nz+1:nd) = 1; np = size(pred_coords, 1); Z_pred = zeros(np, 1); Kvar = zeros(np, 1); for k = 1:np b = zeros(nd+2, 1); for i = 1:nz b(i) = variogram_value(norm(Zcoords(i,:) - pred_coords(k,:)), pzz, vg_model); end for j = 1:nv h = norm(Vcoords(j,:) - pred_coords(k,:)); gzz = variogram_value(h, pzz, vg_model); gvv = variogram_value(h, pvv, vg_model); b(nz+j) = rho * sqrt(gzz * gvv); end b(nd+1) = 1; % Z权重无偏约束 b(nd+2) = 0; % V权重无偏约束 w = A \ b; Z_pred(k) = w(1:nz)' * Zdata + w(nz+1:nd)' * Vdata; Kvar(k) = max(sum(w(1:nd) .* b(1:nd)) + w(nd+1), 0); end end

这段代码里最需要注意的地方是块索引。MATLAB里向量索引从1开始，nz和nv混在一起很容易写错位置。我踩过最贵的一次坑就是把拉格朗日行填错了列，结果求出来的权重全部乱套，画出来的插值图直接出现了波浪形伪影。建议在组矩阵前先把维度写注释，像上面这样标清楚“主变量块”“辅助变量块”“拉格朗日行”，或者写个小的索引测试用例，用两个点一组算一遍手推答案来验证。

4. 实战调参经验：块金、变程、搜索邻域与交叉验证

4.1 块金效应和变程：两个最容易误读的参数

块金值（nugget）在理论上描述的是“距离为0时的方差”，包含测量误差和微观尺度变异。很多新手拟合时喜欢把块金压到0，觉得这样拟合更漂亮，实际上这是个危险动作。块金为0意味着变异函数在原点处从0开始连续上升，天然假设了“距离很近的地方属性几乎一样”，但真实采样数据往往不满足这个条件。更麻烦的是，块金过小会导致克里金矩阵条件数恶化，解出来的权重波动大，插值图出现明显的“斑点”。

变程（range）的意义更直观：超过这个距离，两个点之间就没有空间相关性了，变异函数稳定在基台值。变程估得太小，插值结果几乎等于样点值的局部映射，外推区域迅速回落到均值；变程估得太大，每个预测点都会受到远处不相关点的影响，局部细节被过度平滑。我判断变程合不合适的方法是：把拟合曲线和经验散点画在同一张图里，看变程附近是不是刚好是散点开始变平的位置。如果模型在散点还在上升时就到顶了，说明变程低估；反之如果散点早就平了而模型还在长，那就高估了。

4.2 搜索邻域：别让远处的不相关点进来捣乱

理论上克里金可以用全部样点来解方程组，但全量求解有两个坏处：一是n大了矩阵求解慢，二是变程范围之外的点贡献很小却会稀释局部权重。实践中更普遍的做法是每个待估点只取周围固定数量的最近样点参与克里金求解，比如取最近的24个、48个点，或者限制搜索半径等于1.5倍的变程。这就是搜索邻域（search neighborhood）。

在MATLAB里，实现搜索邻域的办法就是在循环里对每个预测点算一次到所有样点的距离，排序后截取前K个点及对应的坐标、观测值，再调用ordinary_kriging或cokriging_predict。K选多少没有硬性标准，样点总量少就全用，样点总量大建议K取50到100之间。要注意的是，如果变程比较小，K取太多会把变程外的点也拉进来，得不偿失。

4.3 用留一交叉验证对比普通克里金和协同克里金

两种插值方法谁更好，不能靠肉眼，要靠交叉验证说话。留一法（LOO）每次去掉一个样点，用剩下的点预测它的值，最后统计所有样点的预测误差。我最常用的指标有三个：平均绝对误差（MAE）、均方根误差（RMSE）、预测值与实测值的决定系数R²。MAE和RMSE越小越好，R²越接近1越好。

这里有一个容易忽略的点：做交叉验证时，变异函数参数要在每次去掉一个点后重新估计，而不是用全量数据拟合好的固定参数。严格来说这样才公平，因为被去掉的点的信息已经不在变异函数里了。如果只是为了快速对比，可以只在第一次拟合后固定参数，但论文投稿最好用前者。我的经验是：两种方式算出来的RMSE差异通常不太大，但审稿人问起来时，能说清楚“重新拟合了变异函数”会让文章严谨不少。

5. 完整脚本：从模拟数据到插值成图

最后给一个可以直接复制运行的完整示例，包括了随机生成样点、拟合变异函数、普通克里金、协同克里金、留一交叉验证比对和画图。跑一遍这个脚本，就能直观看到普通克里金和协同克里金在结果上的差别，也能顺着代码检查自己的数据哪里出了问题。

clear; clc; rng(42); % 生成模拟采样点：在[0,100]x[0,100]区域内取30个随机站点 site = rand(30, 2) * 100; Z = 20*sin(site(:,1)/30) + 15*cos(site(:,2)/25) + 5; V = 2*Z + randn(30, 1)*4; % 辅助变量：和主变量强相关但有噪声 rho = corr(Z, V); fprintf('主辅变量相关系数 rho = %.3f\n', rho); % 待预测网格：10到90，步长10 [gx, gy] = meshgrid(10:10:90); pred_coords = [gx(:), gy(:)]; % 拟合主变量和辅助变量的变异函数 [lag1, gam1] = empirical_variogram(site, Z, 18); pzz = fit_variogram(lag1, gam1, 'spherical'); [lag2, gam2] = empirical_variogram(site, V, 18); pvv = fit_variogram(lag2, gam2, 'spherical'); % 普通克里金 [Z_ok, var_ok] = ordinary_kriging(pred_coords, site, Z, 'spherical', pzz); % 协同克里金 [Z_ck, var_ck] = cokriging_predict(pred_coords, site, Z, site, V, 'spherical', pzz, pvv, rho); % 对比两套结果的统计量 fprintf('普通克里金: 均值 %.2f, 方差范围 [%.2f, %.2f]\n', mean(Z_ok), min(var_ok), max(var_ok)); fprintf('协同克里金: 均值 %.2f, 方差范围 [%.2f, %.2f]\n', mean(Z_ck), min(var_ck), max(var_ck)); % 绘图 figure; subplot(1,3,1); scatter(site(:,1), site(:,2), 40, Z, 'filled'); colorbar; title('样点观测值'); axis equal; xlim([0 100]); ylim([0 100]); subplot(1,3,2); surf(gx, gy, reshape(Z_ok, size(gx)), 'EdgeColor', 'none'); colorbar; view(2); title('普通克里金插值'); axis equal; xlim([0 100]); ylim([0 100]); subplot(1,3,3); surf(gx, gy, reshape(Z_ck, size(gx)), 'EdgeColor', 'none'); colorbar; view(2); title('协同克里金插值'); axis equal; xlim([0 100]); ylim([0 100]); % 留一交叉验证 mae_ok = 0; rmse_ok = 0; mae_ck = 0; rmse_ck = 0; n = size(site, 1); for i = 1:n tr_idx = true(n, 1); tr_idx(i) = false; % 普通克里金 ok_i = ordinary_kriging(site(i,:), site(tr_idx,:), Z(tr_idx), 'spherical', pzz); % 协同克里金 ck_i = cokriging_predict(site(i,:), site(tr_idx,:), Z(tr_idx), site(tr_idx,:), V(tr_idx), 'spherical', pzz, pvv, rho); mae_ok = mae_ok + abs(ok_i - Z(i)); rmse_ok = rmse_ok + (ok_i - Z(i))^2; mae_ck = mae_ck + abs(ck_i - Z(i)); rmse_ck = rmse_ck + (ck_i - Z(i))^2; end mae_ok = mae_ok/n; rmse_ok = sqrt(rmse_ok/n); mae_ck = mae_ck/n; rmse_ck = sqrt(rmse_ck/n); fprintf('普通克里金 LOO-MAE=%.3f LOO-RMSE=%.3f\n', mae_ok, rmse_ok); fprintf('协同克里金 LOO-MAE=%.3f LOO-RMSE=%.3f\n', mae_ck, rmse_ck);

这段代码里的模拟数据是故意造出来的：主变量有一个大尺度的趋势，辅助变量和主变量相关系数很高但不是完全线性。跑完后通常能看到协同克里金的RMSE比普通克里金低一截，同时插值图的表面更连贯，这是因为辅助变量提供了额外的空间结构约束。如果你的真实数据跑出来协同克里金反而不如普通克里金，大概率是主辅变量相关性太弱，或者交叉变异函数近似不成立。这时别硬调参，回到数据本身去重新审视变量选取更靠谱。

在实际项目里，我通常先跑普通克里金作为基线，再用辅助变量做协同克里金，最后留一交叉验证对比。如果辅助变量的相关系数低于0.5，协同克里金一般提不了多少精度；如果高于0.8，效果会非常明显。这也是我多年用下来的取舍标准，供你参考。

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