简介这份MATLAB代码资源面向地理信息、环境科学与地球科学领域的学习者与研究人员用于实现克里金及协同克里金空间插值解决由已知观测点估计未知区域变量值的问题。压缩包共8个文件以7个m脚本和1个mat数据文件为主脚本涵盖半方差函数拟合、变差函数计算与协同克里金主流程mat文件提供测试数据整体约44KB结构紧凑便于快速上手。资源围绕半方差建模、协方差矩阵构建、插值计算与精度评估等环节展开协同克里金部分引入辅助变量以提升多因素影响下的插值精度读者可借助测试数据复现完整流程并对照RMSE、R²等指标检验结果进而调整模型参数适配自身场景。目前已有2680人学习适合希望从原理到代码实践掌握克里金插值的中高级用户参考。1. 从一组散点到一个连续曲面克里金插值到底能解决什么问题手里有一堆采样点每个点带一个测量值可能是土壤重金属含量、地下水位标高、降雨量、矿石品位。想把它们变成一张连续分布的曲面图或者估算任意未采样位置的数值这就是空间插值要干的事。克里金跟反距离加权、样条函数的区别在于它不只看距离远近还看数据本身的空间自相关结构——半变异函数。换句话说它先“学”一遍你的数据在空间上怎么变化再拿这个规律去插值所以在地学、环境、矿业这些领域克里金几乎是默认选项。这份资源是一套 MATLAB 代码覆盖普通克里金和协同克里金两条线。协同克里金多了一个“辅助变量”比如用高程辅助插值气温、用电阻率辅助插值含水率主变量采样稀疏时靠协变量把精度拉上来。代码结构清晰适合拿来做课程设计、论文实验或者直接嵌入自己的数据处理流程里改。如果你手头正好有一批空间采样数据要出图、要交叉验证这套代码能省掉从零推导公式再翻译成 MATLAB 的时间。2. 代码架构与核心模块半变异函数、权重求解、协同扩展2.1 整体目录与模块划分拿到源码包先别急着跑主脚本。花两分钟把目录结构看清楚后面改起来会顺很多。这套代码的典型组织方式大致如下Kriging_CoKriging/ ├── main_kriging.m % 普通克里金主入口 ├── main_cokriging.m % 协同克里金主入口 ├── variogram/ │ ├── calc_variogram.m % 实验半变异函数计算 │ └── fit_variogram.m % 理论模型拟合球状/指数/高斯 ├── kriging/ │ ├── build_kriging_matrix.m % 构建克里金方程组矩阵 │ ├── solve_weights.m % 求解权重 │ └── predict.m % 逐点预测 ├── cokriging/ │ ├── build_cokriging_matrix.m% 协同克里金矩阵含交叉协方差 │ └── predict_cokriging.m % 协同预测 ├── utils/ │ ├── dist_matrix.m % 距离矩阵计算 │ └── cross_validate.m % 留一交叉验证 └── data/ └── sample_data.mat % 示例数据这个划分的好处是半变异函数、矩阵构建、求解、预测各管各的你想换一个变异函数模型只动fit_variogram.m想换求解方式只动solve_weights.m。常见做法是把参数配置集中在一个结构体里往下传避免全局变量满天飞。2.2 实验半变异函数计算克里金的第一步永远是从数据里估计空间结构。核心公式是γ(h) (1/2N(h)) * Σ [Z(xi) - Z(xj)]²其中 h 是滞后距N(h) 是距离落在 h 附近的数据对数量。MATLAB 实现如下function [h, gamma] calc_variogram(coords, values, n_lags, max_dist) % coords: N×2 坐标矩阵 [x, y] % values: N×1 观测值 % n_lags: 滞后距分组数 % max_dist: 最大滞后距一般取研究区最大距离的一半 N size(coords, 1); h zeros(n_lags, 1); gamma zeros(n_lags, 1); lag_width max_dist / n_lags; % 预计算所有点对距离 D dist_matrix(coords); % N×N 距离矩阵 for k 1:n_lags lower (k-1) * lag_width; upper k * lag_width; pair_sum 0; pair_count 0; for i 1:N for j i1:N if D(i,j) lower D(i,j) upper pair_sum pair_sum (values(i) - values(j))^2; pair_count pair_count 1; end end end if pair_count 0 h(k) (lower upper) / 2; % 取组中值作为代表距离 gamma(k) pair_sum / (2 * pair_count); else h(k) NaN; gamma(k) NaN; end end % 去掉空组 valid ~isnan(h); h h(valid); gamma gamma(valid); end这段代码有几个参数值得注意。n_lags一般取 10 到 20太少曲线粗糙太多每组点数不够、噪声大。max_dist取最大距离的一半是经验做法超过这个范围点对太少半变异函数值不可靠。dist_matrix用向量化算比双重循环快得多数据量大时建议替换成pdist2。2.3 理论变异函数拟合实验半变异函数是一堆散点克里金方程需要的是解析形式。常用模型有三种模型表达式适用场景球状γ(h) C0 C·(1.5h/a - 0.5(h/a)³), h≤a空间自相关随距离线性衰减后截断指数γ(h) C0 C·(1 - exp(-3h/a))自相关缓慢衰减无明确变程高斯γ(h) C0 C·(1 - exp(-3h²/a²))空间连续性强、变化平滑拟合用最小二乘MATLAB 的fit或lsqcurvefit都能做function [model_params, fitted_gamma] fit_variogram(h, gamma, model_type) % model_type: spherical | exponential | gaussian % 初始猜测块金值0基台值max(gamma)变程max(h)/2 p0 [0, max(gamma), max(h)/2]; switch model_type case spherical fit_func (p, x) p(1) p(2) * (1.5*x/p(3) - 0.5*(x/p(3)).^3) ... .* (x p(3)) (p(1)p(2)) * (x p(3)); case exponential fit_func (p, x) p(1) p(2) * (1 - exp(-3*x/p(3))); case gaussian fit_func (p, x) p(1) p(2) * (1 - exp(-3*x.^2/p(3)^2)); end % 非线性最小二乘 opts optimoptions(lsqcurvefit, Display, off); model_params lsqcurvefit(fit_func, p0, h, gamma, ... [0, 0, 0], [inf, inf, max(h)*2], opts); fitted_gamma fit_func(model_params, h); endmodel_params返回三个值块金值 C0、基台值 C0C、变程 a。块金值代表随机误差和微观变异如果拟合出来块金值接近基台值说明数据空间自相关很弱克里金插值效果不会比反距离加权好多少——这是个重要的诊断信号。2.4 普通克里金方程组构建与求解有了变异函数模型对每个待估点 x0构建方程组[ γ(x1,x1) γ(x1,x2) ... γ(x1,xn) 1 ] [ λ1 ] [ γ(x1,x0) ] [ γ(x2,x1) γ(x2,x2) ... γ(x2,xn) 1 ] [ λ2 ] [ γ(x2,x0) ] [ ... ] [ .. ] [ ... ] [ γ(xn,x1) γ(xn,x2) ... γ(xn,xn) 1 ] [ λn ] [ γ(xn,x0) ] [ 1 1 ... 1 0 ] [ μ ] [ 1 ]λ 是权重μ 是拉格朗日乘子。解出权重后预测值就是 Z(x0) Σ λi·Z(xi)。function [z_pred, kriging_var] predict(coords, values, target, model_params, model_type) % target: 1×2 待估点坐标 % 返回预测值和克里金方差 n size(coords, 1); % 构建左侧矩阵 A A zeros(n1, n1); for i 1:n for j 1:n A(i,j) eval_variogram(coords(i,:), coords(j,:), model_params, model_type); end A(i, n1) 1; A(n1, i) 1; end % A(n1, n1) 0 已由 zeros 初始化 % 构建右侧向量 b b zeros(n1, 1); for i 1:n b(i) eval_variogram(coords(i,:), target, model_params, model_type); end b(n1) 1; % 求解 weights A \ b; % 预测 z_pred weights(1:n) * values; % 克里金方差 kriging_var weights * b; endeval_variogram根据模型类型和参数计算两点间的半变异函数值。注意A矩阵是 (n1)×(n1)最后一行一列是拉格朗日约束。求解用 MATLAB 的左除\内部走 LU 分解比手动求逆稳定。如果A接近奇异检查是否有重复点或者变异函数参数是否合理。2.5 协同克里金的扩展逻辑协同克里金的核心区别在于除了主变量的自协方差还要引入主变量与协变量之间的交叉协方差以及协变量自身的自协方差。方程组规模变成 (n1n21)×(n1n21)其中 n1 是主变量点数n2 是协变量点数。function [z_pred, var_pred] predict_cokriging(coords1, vals1, coords2, vals2, ... target, params1, params2, params_cross) % coords1/vals1: 主变量 % coords2/vals2: 协变量 % params_cross: 交叉变异函数参数 n1 size(coords1, 1); n2 size(coords2, 1); n n1 n2; A zeros(n1, n1); % 主变量自协方差块 for i 1:n1 for j 1:n1 A(i,j) eval_variogram(coords1(i,:), coords1(j,:), params1, spherical); end end % 主-协交叉协方差块 for i 1:n1 for j 1:n2 A(i, n1j) eval_variogram(coords1(i,:), coords2(j,:), params_cross, spherical); A(n1j, i) A(i, n1j); end end % 协变量自协方差块 for i 1:n2 for j 1:n2 A(n1i, n1j) eval_variogram(coords2(i,:), coords2(j,:), params2, spherical); end end % 拉格朗日约束 A(:, n1) 1; A(n1, :) 1; A(n1, n1) 0; % 右侧向量 b zeros(n1, 1); for i 1:n1 b(i) eval_variogram(coords1(i,:), target, params1, spherical); end for i 1:n2 b(n1i) eval_variogram(coords2(i,:), target, params_cross, spherical); end b(n1) 1; weights A \ b; z_pred weights(1:n1) * vals1 weights(n11:n1n2) * vals2; var_pred weights * b; end协同克里金能不能比普通克里金好关键看主变量和协变量的相关性。相关系数低于 0.5 的话协变量带来的信息增益有限反而增加计算量和矩阵病态风险。我一般会先算一下两者的皮尔逊相关系数再决定要不要上协同。3. 跑通示例数据从加载数据到出图的完整流程3.1 环境准备与数据格式这套代码对 MATLAB 版本要求不高R2018b 及以上都能跑。不需要额外工具箱核心函数都是基础语法。如果你的数据是 Excel 或 CSV用readmatrix或xlsread读进来整理成两列坐标加一列值的格式。示例数据sample_data.mat里包含三个变量coordsN×2 坐标矩阵、valuesN×1 观测值、grid_x、grid_y待插值网格。加载后先看一眼数据范围和分布load(data/sample_data.mat); fprintf(样本点数: %d\n, size(coords, 1)); fprintf(X 范围: [%.2f, %.2f]\n, min(coords(:,1)), max(coords(:,1))); fprintf(Y 范围: [%.2f, %.2f]\n, min(coords(:,2)), max(coords(:,2))); fprintf(值范围: [%.2f, %.2f]\n, min(values), max(values)); % 快速散点图检查 figure; scatter(coords(:,1), coords(:,2), 40, values, filled); colorbar; title(原始采样点分布);这一步别跳过。我踩过的坑是数据里混了重复点或者坐标单位不统一一个用度一个用米后面变异函数拟合出来完全是乱的排查半天才发现是数据问题。3.2 普通克里金完整运行主脚本main_kriging.m的流程是算实验变异函数 → 拟合理论模型 → 构建网格 → 逐点预测 → 出图。%% main_kriging.m clear; clc; close all; % 1. 加载数据 load(data/sample_data.mat); % 2. 计算实验半变异函数 n_lags 15; max_dist max(pdist(coords)) / 2; [h, gamma] calc_variogram(coords, values, n_lags, max_dist); % 3. 拟合理论模型 [params, fitted] fit_variogram(h, gamma, spherical); fprintf(块金值: %.4f, 基台值: %.4f, 变程: %.4f\n, ... params(1), params(1)params(2), params(3)); % 4. 绘制变异函数拟合图 figure; plot(h, gamma, ko, MarkerFaceColor, k); hold on; plot(h, fitted, r-, LineWidth, 2); xlabel(滞后距 h); ylabel(半变异函数 \gamma(h)); legend(实验值, 球状模型拟合); title(半变异函数拟合); % 5. 构建插值网格 [Xq, Yq] meshgrid(linspace(min(coords(:,1)), max(coords(:,1)), 80), ... linspace(min(coords(:,2)), max(coords(:,2)), 80)); Zq zeros(size(Xq)); Vq zeros(size(Xq)); % 克里金方差 % 6. 逐点预测 for i 1:numel(Xq) [Zq(i), Vq(i)] predict(coords, values, [Xq(i), Yq(i)], params, spherical); end % 7. 出图 figure; contourf(Xq, Yq, Zq, 20, LineColor, none); colorbar; title(普通克里金插值结果); hold on; scatter(coords(:,1), coords(:,2), 15, k, filled); figure; contourf(Xq, Yq, Vq, 20, LineColor, none); colorbar; title(克里金方差分布);跑完这两张图一张是插值曲面一张是方差分布。方差大的地方说明该区域采样点稀疏或者变异函数在该距离上不确定性高这对后续布设采样点有直接参考价值。3.3 协同克里金运行与对比协同克里金的运行流程类似区别在于要准备两套数据。示例里假设coords2和vals2是协变量通常协变量采样密度更高。%% main_cokriging.m clear; clc; close all; load(data/sample_data.mat); % 主变量coords, values % 协变量coords2, vals2假设数据里有 % 检查相关性 r corr(values, interp1(coords2, vals2, coords, linear, extrap)); fprintf(主-协变量相关系数: %.3f\n, r); % 分别拟合变异函数 [h1, g1] calc_variogram(coords, values, 15, max(pdist(coords))/2); [params1, ~] fit_variogram(h1, g1, spherical); [h2, g2] calc_variogram(coords2, vals2, 15, max(pdist(coords2))/2); [params2, ~] fit_variogram(h2, g2, spherical); % 交叉变异函数简化处理用主变量参数近似 params_cross params1; % 网格预测 [Xq, Yq] meshgrid(linspace(min(coords(:,1)), max(coords(:,1)), 60), ... linspace(min(coords(:,2)), max(coords(:,2)), 60)); Zq_co zeros(size(Xq)); for i 1:numel(Xq) [Zq_co(i), ~] predict_cokriging(coords, values, coords2, vals2, ... [Xq(i), Yq(i)], params1, params2, params_cross); end figure; contourf(Xq, Yq, Zq_co, 20, LineColor, none); colorbar; title(协同克里金插值结果);跑完之后把普通克里金和协同克里金的结果放一起对比用cross_validate.m做留一交叉验证看 RMSE 和 MAE 哪个更低。如果协同克里金没有明显优势检查相关系数是不是太低或者协变量本身的空间结构是不是太差。3.4 交叉验证与精度评估cross_validate.m的逻辑是每次留一个点用其余点预测它最后统计误差function [rmse, mae, r2] cross_validate(coords, values, model_params, model_type) n size(coords, 1); pred zeros(n, 1); for i 1:n idx true(n, 1); idx(i) false; [pred(i), ~] predict(coords(idx,:), values(idx), coords(i,:), ... model_params, model_type); end residuals values - pred; rmse sqrt(mean(residuals.^2)); mae mean(abs(residuals)); ss_res sum(residuals.^2); ss_tot sum((values - mean(values)).^2); r2 1 - ss_res / ss_tot; fprintf(RMSE: %.4f, MAE: %.4f, R²: %.4f\n, rmse, mae, r2); endR² 低于 0.5 的话要么变异函数模型选得不对要么数据本身空间结构性太弱。这时候可以试试换指数模型或高斯模型或者检查有没有趋势项需要先剔除泛克里金的思路。4. 避坑与排查跑克里金代码时最容易翻车的几个地方4.1 变异函数拟合出来是条直线现象实验半变异函数散点图看起来没有明显的上升-平稳形态拟合出来的变程接近 0 或者块金值等于基台值。原因数据空间自相关很弱或者坐标单位太大导致滞后距分组不合理。比如坐标用经纬度1 度约 111 公里而你的采样点间距只有几百米所有点对都落在第一个滞后距组里。解决先把坐标投影到平面坐标系UTM 或高斯克吕格确保单位是米。然后调整n_lags和max_dist让每个滞后距组里至少有 30 对点。如果还是不行说明数据本身不适合克里金考虑换成反距离加权或者回归克里金。4.2 矩阵奇异或接近奇异现象A \ b报警告“矩阵接近奇异值”或者结果出现 NaN。原因有重复采样点坐标完全相同导致 A 矩阵中有两行完全相同或者变异函数参数中变程太小导致 A 矩阵条件数爆炸。解决先去重把坐标相同的点合并取均值。然后在求解前检查条件数if cond(A) 1e12 warning(矩阵条件数过大: %.2e检查变异函数参数, cond(A)); end常见做法是在对角线加一个极小的白噪声项如 1e-8这叫“正则化”能显著改善数值稳定性。4.3 协同克里金结果反而不如普通克里金现象辛辛苦苦准备了协变量交叉验证 RMSE 却比普通克里金还高。原因主变量和协变量相关性太弱或者交叉变异函数参数估计不准。协同克里金对方程组的数值条件更敏感参数稍微偏一点权重就乱飞。解决先算相关系数低于 0.5 就别硬上协同了。如果相关性够交叉变异函数不要偷懒直接用主变量的参数用calc_variogram对主-协变量对的差值算一遍。另外协变量点数不要比主变量多太多一般控制在 3 倍以内否则矩阵规模膨胀太快。4.4 插值结果出现明显条带或方块现象插值曲面图上有规则条纹或者方块状伪影。原因搜索邻域设置有问题。如果每次预测只用最近的几个点而点分布不均匀就会出现这种伪影。代码里如果没显式设置搜索半径默认用了全部点那问题可能出在网格分辨率太低。解决把网格分辨率从 60×60 提到 100×100 以上。如果还有点状伪影检查predict函数里是不是用了固定数量的最近邻。常见做法是设置一个搜索半径如变程的 1.5 倍半径内的点才参与预测这样既保证精度又避免远点干扰。4.5 运行速度慢到无法接受现象网格 100×100 要跑十几分钟。原因双重循环逐点预测每次都要重新构建 A 矩阵。实际上 A 矩阵只跟采样点有关跟待估点无关完全可以预计算一次。解决把 A 矩阵的构建提到循环外面循环里只更新右侧向量 b% 预计算 A 矩阵 A build_kriging_matrix(coords, model_params, model_type); A_inv inv(A); % 或者用分解 for i 1:numel(Xq) b build_rhs(coords, [Xq(i), Yq(i)], model_params, model_type); weights A_inv * b; Zq(i) weights(1:end-1) * values; end这样速度能快一个数量级。如果还嫌慢把网格分辨率降下来或者用parfor并行。5. 二次开发与改造换模型、加趋势、接自己的数据5.1 替换变异函数模型代码里默认用了球状模型但实际数据可能更适合指数或高斯。改造方式很简单在fit_variogram.m里加一个分支然后在主脚本里改model_type参数。更灵活的做法是把模型函数句柄作为参数传进去function gamma eval_variogram(p1, p2, params, model_type) h norm(p1 - p2); C0 params(1); C params(2); a params(3); switch model_type case spherical if h a gamma C0 C * (1.5*h/a - 0.5*(h/a)^3); else gamma C0 C; end case exponential gamma C0 C * (1 - exp(-3*h/a)); case gaussian gamma C0 C * (1 - exp(-3*h^2/a^2)); case linear gamma C0 C * h / a; % 简单线性模型 end end加一个linear分支就是几行的事。如果你有自己偏好的模型比如幂函数模型照着格式加就行。5.2 加入趋势项从普通克里金到泛克里金普通克里金假设数据没有全局趋势但实际数据往往有。比如地下水位从北到南逐渐降低这就是趋势。泛克里金的思路是在克里金方程里加入趋势项通常用多项式表示% 在 build_kriging_matrix 里扩展 % 趋势项f(x,y) b0 b1*x b2*y % 方程组增加三行三列 n size(coords, 1); A zeros(n3, n3); % ... 填充协方差部分 ... for i 1:n A(i, n1) 1; A(i, n2) coords(i,1); A(i, n3) coords(i,2); A(n1, i) 1; A(n2, i) coords(i,1); A(n3, i) coords(i,2); end这样解出来的权重会自动适应趋势。改造完之后用交叉验证对比一下R² 提升明显就说明趋势项加对了。5.3 接入自己的数据把data/sample_data.mat替换成自己的数据就行但要注意几点坐标必须是投影坐标单位米不能是经纬度数据里不能有 NaN有的话先插补或剔除如果数据量超过 5000 个点建议先做空间抽稀否则矩阵求解会很慢。我一般会写一个数据检查脚本每次换数据先跑一遍function validate_data(coords, values) assert(size(coords,2) 2, 坐标必须是 N×2 矩阵); assert(size(coords,1) length(values), 坐标和值数量不一致); assert(~any(isnan(coords(:))), 坐标中存在 NaN); assert(~any(isnan(values)), 值中存在 NaN); [~, ia] unique(coords, rows); if length(ia) size(coords,1) warning(存在 %d 个重复点建议去重, size(coords,1) - length(ia)); end fprintf(数据检查通过: %d 个点\n, size(coords,1)); end这个习惯是从一次翻车经历来的当时跑了一晚上结果全是 NaN第二天才发现数据里有个空值。从那以后我每次换数据都强制走一遍validate_data几秒钟的事省一晚上。5.4 出图与导出最后一步是把结果导出成能用的格式。MATLAB 的contourf出图好看但不好编辑建议同时导出原始网格数据% 导出为 CSV output [Xq(:), Yq(:), Zq(:), Vq(:)]; writematrix(output, kriging_result.csv); % 导出为 GeoTIFF如果有 Mapping Toolbox R georefcells([min(Yq(:)), max(Yq(:))], [min(Xq(:)), max(Xq(:))], size(Zq)); geotiffwrite(kriging_result.tif, Zq, R);CSV 可以丢进 GIS 软件或者 Python 里继续处理GeoTIFF 直接就是带坐标信息的栅格数据。如果不需要地理参考导出 CSV 就够了。改造这套代码的核心思路就一条把变异函数、矩阵构建、求解、预测这四个环节解耦每个环节都可以单独替换。你想换模型就动eval_variogram想加趋势就动build_kriging_matrix想换求解器就动solve_weights。改完之后用交叉验证的 RMSE 和 R² 做回归测试确保改动没有让精度掉下来。希望这套代码能帮到你手头的空间插值任务。本文还有配套的精品资源点击获取