MATLAB实现克里金与协同克里金插值：从变异函数拟合到交叉验证-北京尧图网络

简介一套面向地理空间数据插值的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 i1: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); % 变程必须为正 endp0的初始化并不是随便拍的块金从最小值开始是因为理论上距离为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(n1, n1); 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(n1, 1:n) 1; A(1:n, n1) 1; A(n1, n1) 0; np size(pred_coords, 1); Z_pred zeros(np, 1); Kvar zeros(np, 1); for k 1:np b zeros(n1, 1); for i 1:n b(i) variogram_value(norm(coords(i,:) - pred_coords(k,:)), params, vg_model); end b(n1) 1; w A \ b; Z_pred(k) w(1:n) * data; Kvar(k) max(sum(w(1:n) .* b(1:n)) w(n1), 0); end end注意这里可能有人会问不是说克里金用协方差矩阵吗为什么这代码里全是变异函数两者在无偏性约束下是等价的因为C(h) 基台值 - γ(h)而无偏性要求权重和等于1常数基台值会自然抵消。直接组变异函数矩阵可以省掉一次转换代码也更好理解。克里金方差出现负数的情况实测里会碰到主要原因是变异函数模型参数不合理或者矩阵病态。我在这里用max(…,0)兜底但更重要的是后面第4章里讲怎么从根源上避免这个问题。3. 协同克里金MATLAB实现把辅助变量写进方程组3.1 协同克里金的分块矩阵结构协同克里金的本质就是把普通克里金的单变量系统扩展成多变量系统。假设有主变量Znz个采样点和辅助变量Vnv个采样点待估点x₀的估计涉及四个块的变异/交叉变异信息矩阵块含义维度Γzz主变量自身的变异函数值nz × nzΓvv辅助变量自身的变异函数值nv × nvΓzv / Γvz主变量与辅助变量的交叉变异函数值nz × nv / nv × nz拉格朗日行/列无偏性约束Σλ1Σν02 × (nznv2)方程组右侧也有对应的三个部分待估点到主变量样点的变异函数值、待估点到辅助变量样点的交叉变异函数值以及两个无偏性约束条件1和0。整个系统是一个(nznv2) × (nznv2)的分块矩阵组起来不难但索引别搞错。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(nzi, nzj) 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, nzj) rho * sqrt(gzz * gvv); A(nzj, i) A(i, nzj); end end % 无偏性约束Σλ1, Σν0 A(1:nz, nd1) 1; A(nz1:nd, nd1) 0; A(nd1, 1:nz) 1; A(nd1, nz1:nd) 0; A(1:nz, nd2) 0; A(nz1:nd, nd2) 1; A(nd2, 1:nz) 0; A(nd2, nz1:nd) 1; np size(pred_coords, 1); Z_pred zeros(np, 1); Kvar zeros(np, 1); for k 1:np b zeros(nd2, 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(nzj) rho * sqrt(gzz * gvv); end b(nd1) 1; % Z权重无偏约束 b(nd2) 0; % V权重无偏约束 w A \ b; Z_pred(k) w(1:nz) * Zdata w(nz1:nd) * Vdata; Kvar(k) max(sum(w(1:nd) .* b(1:nd)) w(nd1), 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效果会非常明显。这也是我多年用下来的取舍标准供你参考。本文还有配套的精品资源点击获取

MATLAB实现克里金与协同克里金插值：从变异函数拟合到交叉验证

相关推荐

HyperFrames constellation-hub 蓝图：节点环聚中心的镜头编排与三大收尾变体

Matlab实现PCA与BP神经网络的手写字母识别

从 Borrow One 理解 Rust 生命周期标注：绑定返回引用与参数借用的完整指南

更多资讯

SAP PS项目类型与编码方案匹配：OPSK/OPSJ配置与排查实战

GPT-Researcher CURSOR_RULES 解析：AI IDE 规则文件如何为 LLM 自主研究项目锚定结构与开发规范

多目标退火算法求解含P2X综合能源系统调度问题Matlab实现

ToolJet Map 组件深度解析：属性、事件、组件特定动作与源码级实现原理

现在主流的AI写论文工具有哪些品牌？亲测后说说真心话

有道季报解读：营收增17%净利降50%，智能硬件能否成为第二增长曲线

今日更新

CopilotKit Vue Demo 与 React Demo 的 Parity 实现指南：路由、运行时端点与交互语义的全框架对齐

Harbor RBAC 实战：db_auth 模式下非管理员用户管理项目成员的完整验证与权限模型解析

STM32平衡车纯平衡程序从拆包到参数整定的完整调试指南

本周精选

超人会飞不算本事：系统稳定依赖清晰规则与边界设计

超人VS蜘蛛侠：拆解超级IP的影响力与传播方法论

基于CNN的调制信号识别：MATLAB实现时频图分类实战

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？