简介这是一套基于DACE工具箱的克里金Kriging插值MATLAB代码面向地质统计、地下水模拟、土壤制图以及计算机实验设计等领域的科研人员和工程师。克里金法以空间自协方差为基础进行最优插值相比普通插值方法能给出估计方差资源内部实现了常数、线性、二次回归模型与高斯、指数、球形、线性等多类相关模型并封装了拟合与预测两个核心函数便于用户针对不同数据特征灵活选择与组合。包内共十九个文件主要包括十六个脚本文件核心算法与辅助函数、一个说明文档、一个示例数据文件以及一个版本变更记录压缩包总大小约1.48MB整体结构清晰适合快速部署。目前已有4976人学习下载特别适合有一定MATLAB基础、希望深入理解克里金插值原理与实现细节的开发者。借助示例数据和说明文档读者可以快速掌握DACE工具包的函数调用方式对比不同模型参数对插值精度的影响是开展空间数据分析与建模的实用工具。1. 为什么我在MATLAB里自己写Kriging前几天一个做环境监测的朋友问我手头有一批离散的土壤重金属采样点想插值成连续分布图商用GIS软件太贵Python那边生态虽然全但团队不熟问我MATLAB能不能干这事。我直接告诉他能而且Kriging克里金插值本身就是从地质统计里走出来的方法MATLAB做矩阵运算有天然优势自己写一遍逻辑也不算复杂几十行核心代码就能跑通。Kriging解决的问题用大白话说就是你手里只有少数几个点的观测值想知道没测过的位置大概是多少。它跟反距离加权IDW、样条插值最大的区别在于——Kriging不只是一个加权平均公式它还会告诉你这个位置的估计值有多可信也就是Kriging方差。这个性质在环境监测、矿产储量估算、气象站点数据网格化、代理模型构建surrogate model这些场景里非常实用因为你不仅要一个数还要知道这个数的误差范围。这篇文章适合两类人看。一类是刚接触空间插值、需要用MATLAB完成课设或论文实验的在校生另一类是有数据但不确定怎么选插值方法的工程师。我会把整个流程拆开讲从变异函数怎么算、怎么拟合到Kriging方程组怎么解再到代码怎么组织、踩过哪些坑一步不落。2. 核心思路拆解Kriging的那套数学流程先说清楚一件事Kriging不是某一个固定公式它是一族基于变异函数variogram的最小方差无偏估计算法。普通克里金Ordinary Kriging是最常用的一种它假设区域化变量的均值未知但恒定这也是我下面代码里默认的实现方式。如果你要处理有明确趋势的数据比如高程随距离线性变化那要用泛克里金Universal Kriging如果只是想知道空间格局普通克里金基本够用。2.1 从数据到经验半方差Kriging的底层逻辑是空间上离得越近的点其观测值越相似。这个相似度随距离衰减的规律用半方差semivariance来量化。对任意两个采样点 i 和 j它们之间的半方差定义为γ(h) 0.5 × (z_i - z_j)²其中 h 是两点之间的距离z_i、z_j 是对应位置的观测值。把所有点对的距离算出来、按距离分箱bin每个箱子里取平均半方差就得到经验半方差散点图。这一步是Kriging的观察阶段——你先看数据在空间上的相关性能维持多远的距离。分箱的间距怎么定很关键。我一般把最大距离的 1/2 到 1/3 作为分箱数量参考每个箱子保证至少有 30 对点不然个别离群点对会把均值带偏。MATLAB里用 pdist2 算距离矩阵很直接然后配合 histcounts 或者手写循环分箱都行数据量在几千个点以内性能完全不是问题。2.2 变异函数模型的拟合拿到经验半方差散点图之后需要用理论模型去拟合它因为插值时要对任意距离求半方差不能只靠离散的箱均值。常用的理论模型有三种模型公式特点球状模型γ(h) C0 C1×(1.5h/a - 0.5(h/a)³)h≤a工业界最常用有明确变程指数模型γ(h) C0 C1×(1 - exp(-h/a))渐近趋向基台变程约为3a高斯模型γ(h) C0 C1×(1 - exp(-(h/a)²))曲线平滑适合连续性强变量这里的 C0 叫块金值nugget代表测量误差或微观尺度上的变异C1 是基台值sill减去块金值的部分代表空间结构性变异a 是变程range表示超过这个距离后空间相关性消失。拟合方式可以用最小二乘也可以用加权最小二乘——离得近的点对数量多、半方差更可靠所以通常按每个箱子里的点对数加权。我在MATLAB里用 fminsearch 做这个拟合目标函数就是加权残差平方和。初值选择上有个经验C0 取最小半方差的 0.8 倍左右C1 取 (箱均值最大值 - 最小半方差)a 取最大分箱距离的 1/3。这个初值给得好fminsearch 基本几十步就收敛给不好就会陷在局部最优里跑出来的模型完全不像样。2.3 普通克里金方程组拟合好变异函数后就可以对任意待插值点 x0 做估计了。核心思路是求一组权重 λ_i使得估计值 ẑ(x0) Σλ_i z_i 满足两个条件无偏权重之和等于1和估计方差最小。这是一个带等式约束的最优化问题用拉格朗日乘子法可以转化为线性方程组[[Γ, 1], [1ᵀ, 0]] × [[λ], [μ]] [[γ0], [1]]其中 Γ 是已知采样点之间的半方差矩阵Γ_ij γ(h_ij)γ0 是待插值点到已知点的半方差向量μ 是拉格朗日乘子。求解这个方程组得到权重后估计值就是权重的线性组合Kriging方差也可以用同样的矩阵元素算出来。这个方程组看着简单实际求解时有个非常容易踩的坑Γ 矩阵有时会接近奇异特别是有两个采样点距离极近或几乎重合时。后面我会专门讲这个问题怎么处理。3. 完整MATLAB实现从半方差到插值出图3.1 主程序框架我把整个流程组织成三个函数一个负责计算经验半方差一个负责拟合变异函数模型一个负责对给定网格做Kriging插值。这样职责清晰后面换模型参数、换数据集都方便。% 主脚本kriging_demo.m % 数据格式X为n×2矩阵(坐标)z为n×1向量(观测值) % 1. 生成测试数据模拟采样点 rng(42); n 60; X rand(n, 2) * 10; % 区域 [0,10]×[0,10] % 用一个已知函数生成真实场加上噪声模拟采样误差 trueField (x) 3*sin(0.8*x(:,1)) .* cos(0.6*x(:,2)) 2; z trueField(X) randn(n,1)*0.15; % 2. 计算经验半方差 [distBins, semivar, nPairs] computeVariogram(X, z, 15); % 3. 拟合变异函数模型球状模型 model fitVariogram(distBins, semivar, nPairs); % 4. 生成插值网格 [Xg, Yg] meshgrid(0:0.2:10, 0:0.2:10); Xq [Xg(:), Yg(:)]; % 5. Kriging插值 [Zq, varZq] ordinaryKriging(X, z, Xq, model); % 6. 可视化 figure(Position, [100 100 1200 500]); subplot(1,2,1); scatter(X(:,1), X(:,2), 40, z, filled); colorbar; axis equal; title(采样点观测值); subplot(1,2,2); Zgrid reshape(Zq, size(Xg)); imagesc(0:0.2:10, 0:0.2:10, Zgrid); set(gca, YDir, normal); colorbar; axis equal; title(Kriging插值结果);核心的变异函数计算我这里用computeVariogramfunction [hMean, gammaMean, nPairs] computeVariogram(X, z, numBins) % 计算经验半方差 % X: n×2坐标矩阵, z: n×1观测值, numBins: 分箱数量 D pdist2(X, X); % 两两距离矩阵 dz2 0.5 * (z - z).^2; % 两两半方差矩阵 % 只取上三角避免重复计数 [I, J] ndgrid(1:size(X,1), 1:size(X,1)); mask I J; distVec D(mask); dzVec dz2(mask); maxDist max(distVec); edges linspace(0, maxDist, numBins1); hMean zeros(numBins, 1); gammaMean zeros(numBins, 1); nPairs zeros(numBins, 1); for k 1:numBins idx (distVec edges(k)) (distVec edges(k1)); if sum(idx) 0 hMean(k) mean(distVec(idx)); gammaMean(k) mean(dzVec(idx)); nPairs(k) sum(idx); else hMean(k) 0.5 * (edges(k) edges(k1)); gammaMean(k) NaN; nPairs(k) 0; end end % 去掉空箱子 valid ~isnan(gammaMean) (nPairs 5); hMean hMean(valid); gammaMean gammaMean(valid); nPairs nPairs(valid); end这里有个细节pdist2算出来的是 n×n 矩阵内存占用是 O(n²)。n 在几千以内完全没问题但如果采样点超过一两万个就得改用分块计算或者pdist返回向量形式。我处理过三次采样点接近五万的项目直接pdist2会报内存不足那次是把区域切块分治处理的细节后面有机会再展开。3.2 变异函数拟合参数估计的关键代码接下来是拟合球状模型。球状模型的半方差公式是分段函数在h a时取基台值C0C1。目标函数用点对数加权的最小二乘点对数越多的箱子权重越大function model fitVariogram(h, gamma, nPairs) % 用加权最小二乘拟合球状模型 % 返回结构体: model.C0, model.C1, model.a % 初值估计 C0_init max(0.05, 0.8 * min(gamma)); C1_init max(gamma) - min(gamma); a_init max(h) * 0.5; % 目标函数加权残差平方和 objFun (theta) wssError(theta, h, gamma, nPairs); options optimset(Display, off, MaxIter, 5000, ... MaxFunEvals, 10000); theta0 [C0_init, C1_init, a_init]; % 参数约束C00, C10, a0 lb [0, 0, 0.01]; ub [inf, inf, max(h)*3]; [theta_opt, ~] fmincon(objFun, theta0, [], [], [], [], lb, ub, [], options); model.C0 theta_opt(1); model.C1 theta_opt(2); model.a theta_opt(3); % 画拟合效果图 figure; plot(h, gamma, ko, MarkerFaceColor, k); hold on; hFine linspace(0, max(h)*1.05, 200); plot(hFine, sphericalVariogram(hFine, model), r-, LineWidth, 1.5); xlabel(距离 h); ylabel(半方差 \gamma(h)); legend(经验半方差, 球状模型拟合, Location, northwest); title(变异函数拟合结果); end function wse wssError(theta, h, gamma, nPairs) % 加权残差平方和目标函数 C0 theta(1); C1 theta(2); a theta(3); gammaPred sphericalVariogram(h, [C0, C1, a]); wse sum(nPairs .* (gamma - gammaPred).^2) / sum(nPairs); end function g sphericalVariogram(h, theta) % 球状模型半方差 C0 theta(1); C1 theta(2); a theta(3); g zeros(size(h)); idx h a; g(idx) C0 C1 * (1.5*h(idx)/a - 0.5*(h(idx)/a).^3); g(~idx) C0 C1; end为什么这里我用fmincon而不是fminsearch因为fminsearch是无约束优化迭代过程中a可能变成负数球状模型公式直接崩掉。虽然可以在目标函数里做惩罚但用带边界约束的fmincon更干净。另外fmincon默认算法是内点法对这个只有三个参数的平滑问题收敛非常快完全不用担心性能。一个容易被忽视的点球状模型在h a处的连续性。拟合出来的a如果落在某个箱子里那段区间可能出现微小的不连续但对最终插值结果影响很小因为Kriging方程里用到的是具体的γ(h)值而不是导数不需要担心这个问题。3.3 求解Kriging权重并做插值核心的插值函数如下function [Zq, varZq] ordinaryKriging(X, z, Xq, model) % 普通克里金插值 % X: n×2 采样点坐标 % z: n×1 采样值 % Xq: m×2 待插值点坐标 % model: 变异函数模型结构体 n size(X, 1); m size(Xq, 1); % 采样点之间的半方差矩阵 Gamma (n×n) D pdist2(X, X); Gamma sphericalVariogram(D(:), [model.C0, model.C1, model.a]); Gamma reshape(Gamma, n, n); % 构建Kriging矩阵 A (n1 × n1) A zeros(n1, n1); A(1:n, 1:n) Gamma; A(n1, 1:n) 1; A(1:n, n1) 1; A(n1, n1) 0; % 为避免矩阵奇异加一个小扰动 A(1:n, 1:n) A(1:n, 1:n) 1e-10 * eye(n); % 预分解矩阵加速多个点的插值对同一组采样点只需分解一次 [L, U, p] lu(A, vector); b zeros(n1, m); Zq zeros(m, 1); varZq zeros(m, 1); for i 1:m % 待插值点到采样点的半方差向量 d0 sqrt((X(:,1) - Xq(i,1)).^2 (X(:,2) - Xq(i,2)).^2); gamma0 sphericalVariogram(d0, [model.C0, model.C1, model.a]); % 解方程组求权重和拉格朗日乘子 rhs [gamma0; 1]; sol zeros(n1, 1); sol(p) U \ (L \ rhs); lambda sol(1:n); mu sol(n1); % 插值估计 Zq(i) sum(lambda .* z); % Kriging方差 varZq(i) sum(lambda .* gamma0) mu; end end这里有个性能细节A矩阵只和采样点位置有关跟待插值点无关所以可以用lu做一次矩阵分解后面每个待插值点只需要做一次回代复杂度从 O(n³ m·n²) 降到 O(n³ m·n²)分解一次 O(n³)每次回代 O(n²)。如果 m 是几万个网格点这个优化能把运行时间从分钟级降到秒级。千万别在循环里对每个点重新\一次那样等于重复做 n 次 LU 分解效率差非常多。矩阵加1e-10 * eye(n)是防止半方差矩阵对角线全是 C0 时可能出现的奇异问题。严格来说克里金矩阵是正定的但数值上如果两个点距离非常接近矩阵行会近似线性相关加了这个小扰动可以显著提升数值稳定性。这个值不能加太大否则解出来的权重会偏离真实解。3.4 结果验证交叉验证不能省插值做完不等于能用。我每次都做留一交叉验证Leave-One-Out Cross Validation把第 i 个采样点当作未知用剩下的 n-1 个点去估计它的值然后对比估计值和真实值。这样可以量化Kriging的预测误差也能顺带检查经验半方差拟合得合不合理。% 留一交叉验证简版 n size(X, 1); predErr zeros(n, 1); for i 1:n idxTrain true(n, 1); idxTrain(i) false; [zPred, ~] ordinaryKriging(X(idxTrain,:), z(idxTrain), X(i,:), model); predErr(i) z(i) - zPred; end rmse sqrt(mean(predErr.^2)); fprintf(交叉验证RMSE: %.4f\n, rmse); % 画散点预测值 vs 真实值 figure; scatter(z, z - predErr, 30, filled); hold on; plot([min(z), max(z)], [min(z), max(z)], r--); xlabel(真实值); ylabel(预测值); title(sprintf(交叉验证结果 RMSE%.4f, rmse));如果 RMSE 明显大于观测值的标准差基本可以判断变异函数拟合有问题或者样本量太小导致空间结构估计不可靠。我见过有人插值结果图非常漂亮一交叉验证 RMSE 惨不忍睹就是因为只盯着插值图看完全没过验证这一关。写论文或者出报告之前这个数字必须记录在案。4. 常见问题与排查技巧实录4.1 矩阵奇异或条件数过大这是Kriging报错里最高频的一个现象是解出来的权重数值巨大、正负交替插值结果出现明显的牛眼或者大片异常值。原因通常有两种一是采样点里有位置几乎重合的数据二是模式参数里C0特别小导致矩阵对角线接近零。排查方法先检查数据里有没有重复坐标点有的话合并或删除。然后看变异函数拟合出来的C0如果接近 0矩阵对角线上就全靠距离近的点对来支撑主对角优势数值上很容易出问题。我一般会在构建矩阵时加1e-8到1e-12级别的正则项同时用cond(A)检查条件数条件数超过 1e12 就重点检查数据。4.2 拟合的变程超出数据范围fmincon跑了半天拟合出来的a等于上界max(h)*3说明经验半方差在整个距离范围内几乎都还在上升没有到达平台期。这个问题要么是数据本身就缺乏大尺度空间相关性那就是纯块金效应Kriging退化为均值估计要么是采样范围太小根本观测不到完整的空间结构。这时候我通常的做法是先画经验半方差散点图肉眼看趋势。如果确实是单调上升到最大距离都没平台我会直接告诉用户这个数据不适合做Kriging或者建议扩大采样范围。硬拟合一个巨大的变程插值结果跟反距离加权差不多Kriging的优势完全体现不出来。4.3 各向异性数据的坐标变换很多空间数据在不同方向上的相关性尺度不一样比如地质构造走向方向上相关性延续很远垂直方向上衰减很快。如果用各向同性变异函数去拟合会把两个方向混在一起平均拟合出来的模型既不贴近主方向也不贴近垂直方向。处理方式不复杂先算各方向的经验半方差比如0°、45°、90°、135°四个方向看有没有明显差异。如果有最省事的方法是在计算距离之前把坐标做线性变换X X * R其中R是旋转-缩放矩阵然后对变换后的坐标用各向同性变异函数。MATLAB里用rotz或手写旋转矩阵都可以关键是两个方向的变程比要确定好这个比例可以从四个方向的半方差图上大致读出来。4.4 插值结果出现负值或低于物理下限比如插值土壤含水量结果出现负的百分比这通常是Kriging权重出现负值导致的。克里金权重可以有负值这是数学上的正常现象但物理上不能接受。几种处理方法最简单的是改用对数变换对log(z)做Kriging再变换回来能保证正值或者用普通克里金求解后对权重做非负约束MATLAB里用lsqlin加不等式约束做一些处理可以这样做还有一种做法是用指示克里金Indicator Kriging专门处理这类非高斯、有物理边界的数据。不过要说明的是加约束的克里金会损失一部分最优性插值方差会变大这是代价。5. 写在最后的几个实操心得这个话题写到这儿核心的东西已经全部过了一遍。回顾我这些年在MATLAB里折腾Kriging的过程最大的体会是Kriging的门槛不在矩阵求解而在变异函数这一步。矩阵求解是确定的数学写对了就有结果变异函数拟合却需要你对数据有手感——分箱分得合不合理、模型选得对不对、初值给得好不好都直接影响后续所有结果的质量。还有一点想提醒的是很多初学者上来就追求高级的泛克里金、协克里金但普通克里金都没跑通就去加复杂度最后代码调不出来还找不到问题在哪。实际工程项目里80%的情况普通克里金完全够用把变异函数拟合做扎实、交叉验证做规范结果已经相当能打了。如果你手头也有空间采样数据建议按我上面的代码流程自己跑一遍不要直接复制粘贴就完事——把分箱数改成不同值看看结果变多少换指数模型和高斯模型对比一下拟合效果给数据加点噪声看看鲁棒性。这些实验做完Kriging在你脑子里就不是一个黑箱子了后面遇到再复杂的问题也知道从哪个环节入手去改。本文还有配套的精品资源点击获取