简介克里金Kriging插值算法的MATLAB实现面向需要开展空间插值或地质统计建模的研究人员与工程师可解决克里金插值、变差函数拟合与空间预测的编程问题。代码源自地质统计学中常用的DACE工具箱思路适用于地下水模拟、土壤制图、环境监测等领域的网格化空间数据处理场景。压缩包共19个文件以16个M脚本为主并附有PDF说明文档、MAT数据文件及更新日志整体仅1.48MB轻量易用。M文件覆盖dacefit、predictor、gridsamp、corrgauss、correxp、corrlin等功能模块同时包含线性与二次回归模型及多种相关函数从模型训练到空间预测均有对应实现data1.mat提供实测样例用于快速验证PDF文档对算法原理和调用方式做了系统梳理便于使用者对照排错。目前已有4976人浏览学习可作为课程设计、科研项目或工程应用中的参考实现适合具备一定MATLAB基础并希望将克里金方法落地到实际项目中的中高级使用者。 在MATLAB里做空间插值大多数人第一反应是scatteredInterpolant或者interp2。但如果你需要的不只是一张漂亮的等值线图还需要每个预测点上的不确定性——也就是预测方差——那常规插值就完全无能为力了。这就是我研究Kriging算法的起点。Kriging这个名字最早来自南非矿业工程师Danie Krige的地质统计实践后来在气象、土壤、环境监测、实验设计响应面分析等一堆领域都成了标配方法它最核心的价值是在所有线性无偏估计器里Kriging给出的预测方差最小所以它有一个很硬核的称号——最优线性无偏预测BLUP。这篇文章我把我自己在MATLAB里从零手写Kriging算法普通克里金的完整思路、代码、验证过程、踩坑心得都整理了出来。适合这几类人看不想用现成工具箱、想彻底搞懂Kriging内部原理的人做空间插值但发现常规方法给不出误差估计的人以及手里有一堆离散点数据、需要生成平滑网格场的学生和工程师。1. Kriging最值钱的地方不只是插值是带方差输出的插值1.1 从一个数据插值场景说起假设你在某个区域布了40个采样点测了土壤重金属含量现在想画一张整个区域的污染分布图。用scatteredInterpolant三行代码就能得到一张图看起来也没问题。但紧接着一个问题就把我噎住了插值结果在每个位置到底有多可信采样点密集的地方肯定比稀疏的地方可信但是常规插值函数一概不告诉你这些。Kriging不一样它天生自带一套机制能在输出预测值的同时输出克里金方差。这个方差不是拍脑袋给的而是基于数据自身的空间相关性推出来的直接反映了这个位置预测得有多准。还有一个更实际的问题常规插值方法比如反距离加权IDW给的权重是凭经验定的距离越近权重越大但多大是人为指定的。Kriging不这么干它通过变异函数把数据点之间的空间关系定量地挖出来然后用这个空间关系自动求解每个已知点对未知点的最优权重。所以它不叫插值叫预测语义上也是有讲究的。1.2 半变异函数Kriging的地基Kriging的一切都是建立在半变异函数semivariogram之上的。半变异函数描述的是两个点之间的距离越近它们的值越相似距离越远差异越大。这种相似程度随距离衰减的关系就是空间自相关性。它的数学定义是[ \gamma(h) \frac{1}{2N(h)} \sum_{i1}^{N(h)} (z(x_i) - z(x_ih))^2 ]其中(h)是两个点之间的距离(N(h))是距离落在(h)附近的所有点对数量(z(x_i))是位置(x_i)处的观测值。除以2是有讲究的这个半字就是这么来的。这个式子其实很好理解你可以把(z(x_i) - z(x_ih))想象成两个邻居的数值差异差异越大说明这个距离上数值变化越剧烈变异函数值就越大。把不同距离上的变异函数值画出来就会看到一条通常单调递增的曲线距离越远差异越大直到一个平台——超过这个距离之后两个点的数值基本就互不相干了。1.3 简单克里金、普通克里金、泛克里金怎么选Kriging有一整个家族MATLAB里最常碰到的是这三个简单克里金Simple Kriging假设研究区域内的均值是已知的常数。这个在现实中很难满足因为我们做插值恰恰是因为不知道空间分布怎么可能提前知道全局均值所以实际用得少。普通克里金Ordinary Kriging假设均值是未知的常数在求解权重时加一个约束条件让所有权重之和等于1。这是应用最广的版本也是我这篇文章的默认版本。泛克里金Universal Kriging假设均值不是常数而是随坐标变化的趋势项比如线性趋势、二次趋势。适用于数据有明显空间漂移的场景比如地形高程从南到北逐渐抬升。对大多数第一次接触Kriging的人来说先把普通克里金吃透就够了。泛克里金的内核跟普通克里金几乎一样只是在方程组里额外多了一些趋势项的基函数和约束条件理解了普通克里金再看泛克里金就是顺水推舟的事。2. MATLAB手写Kriging前先把这个数学映射理顺2.1 从变异函数到权重矩阵的数学映射Kriging的推导过程在各种教材里能写十几页但落到代码层面核心就一件事解一个线性方程组。假设有(n)个观测点已知它们的坐标和值现在要预测任意一个位置(x_0)的值。Kriging假设预测值是所有观测值的线性组合[ \hat{z}(x_0) \sum_{i1}^{n} \lambda_i z(x_i) ]关键问题是权重(\lambda_i)怎么确定。Kriging的思路是让预测方差最小同时保证无偏。经过拉格朗日乘子法推导最终得到的是一个( (n1) \times (n1) )的线性方程组[ \begin{bmatrix} \gamma(x_1-x_1) \gamma(x_1-x_2) \cdots \gamma(x_1-x_n) 1 \ \gamma(x_2-x_1) \gamma(x_2-x_2) \cdots \gamma(x_2-x_n) 1 \ \vdots \vdots \ddots \vdots \vdots \ \gamma(x_n-x_1) \gamma(x_n-x_2) \cdots \gamma(x_n-x_n) 1 \ 1 1 \cdots 1 0 \end{bmatrix} \begin{bmatrix} \lambda_1 \ \lambda_2 \ \vdots \ \lambda_n \ \mu \end{bmatrix}\begin{bmatrix} \gamma(x_0-x_1) \ \gamma(x_0-x_2) \ \vdots \ \gamma(x_0-x_n) \ 1 \end{bmatrix} ]简单说就是左边的大矩阵描述的是已知点与已知点之间的空间关系右边的向量描述的是待预测点与已知点之间的空间关系解出来的权重就是让两种关系最协调的权重。那个额外的(\mu)是拉格朗日乘子它没有实际物理意义但在求解过程中保证了无偏性约束预测方差的计算里还会用到它。2.2 理论变异函数的三种常用模型实验算出来的变异函数是一堆离散的点不能直接拿来解方程组。你需要用一个数学函数去拟合它这个函数就是理论变异函数模型。MATLAB手写Kriging时最常用的是这三个模型公式特点球形模型(\gamma(h) C_0 C(1.5\frac{h}{a} - 0.5(\frac{h}{a})^3))(h \le a)有明确的变程到达变程后平稳指数模型(\gamma(h) C_0 C(1 - e^{-3h/a}))渐近到达基台值没有硬边界高斯模型(\gamma(h) C_0 C(1 - e^{-3(h/a)^2}))原点处行为像抛物线适合平滑连续场公式里的三个参数请你务必记牢它们是Kriging调参的核心块金值nugget(C_0)距离为0时的变异函数值。理论上距离为0时变异函数应该也为0但由于测量误差和微观尺度变异实际数据里往往有一个小的跳跃。块金值越大说明数据噪声越大。基台值sill(C_0 C)变异函数达到平台时的值反映数据的整体方差。变程range(a)变异函数达到基台值时的距离。超过这个距离点与点之间就没有空间相关性了。2.3 把方程组写成矩阵形式解这个方程组在MATLAB里非常直接就是一次矩阵左除lambda A \ b;MATLAB的\运算符会自动选择合适的求解方法一般是LU分解矩阵对称时可能走Cholesky比inv(A)*b快得多数值稳定性也更好。我自己在写代码的时候几乎从来不用inv去解线性方程组这是MATLAB使用的一个基本素养。这里要提醒一下很多人第一次接触Kriging会花很多时间纠结数学公式推导反而忽略了代码实现。我的经验是先写代码跑通再回头补数学。当你看到矩阵方程和代码里的一一对应关系时那些公式自然就理解了。3. 完整可用的MATLAB代码可直接复制3.1 实验半变异函数的MATLAB实现写Kriging的第一步不是预测而是算实验变异函数看看数据到底有没有空间相关性。下面是完整的函数代码function [h_mean, gamma_hat, n_pairs] experimental_variogram(x, y, z, n_bins) % 计算实验半变异函数 % 输入 % x, y - 观测点坐标列向量 % z - 观测值列向量 % n_bins - 距离分箱数量 % 输出 % h_mean - 每个箱子的平均距离 % gamma_hat- 每个箱子的平均半变异值 % n_pairs - 每个箱子的点对数 n length(z); if n 4 error(样本点太少无法计算变异函数); end % 计算所有点对之间的距离矩阵 D pdist2([x(:) y(:)], [x(:) y(:)]); % 计算所有点对之间的半变异值注意除以2 Z_diff2 (z(:) - z(:)) .^ 2 / 2; % 只取上三角元素去掉对角线和重复点对 tri_idx triu(true(n), 1); distances D(tri_idx); semivars Z_diff2(tri_idx); max_d max(distances); h_mean zeros(1, n_bins); gamma_hat zeros(1, n_bins); n_pairs zeros(1, n_bins); for k 1:n_bins lower (k - 1) * max_d / n_bins; upper k * max_d / n_bins; mask distances lower distances upper; n_pairs(k) sum(mask); if n_pairs(k) 0 h_mean(k) mean(distances(mask)); gamma_hat(k) mean(semivars(mask)); end end % 去掉没有点对的空箱子 keep n_pairs 0; h_mean h_mean(keep); gamma_hat gamma_hat(keep); n_pairs n_pairs(keep); end这段代码里有两个细节值得说一是Z_diff2 (z(:) - z(:)) .^ 2 / 2这里用了MATLAB的隐式展开一个列向量减一个行向量得到一个矩阵。这里用到了点乘和直接乘的区别.^ 2是元素级别的幂运算如果写成^ 2MATLAB会试图做矩阵乘法结果完全不一样。二是triu(true(n), 1)这个掩膜技巧一次性把所有的点对距离和半变异值取出来避免了写双层for循环数据量几百个点时性能毫无压力。3.2 变异函数拟合用不依赖工具箱的最小实现拿到实验变异函数的离散点之后就需要拟合理论模型。很多教材推荐用lsqcurvefit拟合但这个方法需要Optimization Toolbox。考虑到不少用户用的MATLAB版本未必装了全量工具箱我写了一个基于fminsearch的最小实现不依赖任何工具箱function gamma_model_val variogram_model(h, model, params) % 计算理论变异函数值 % model: spherical, exponential, gaussian % params: [nugget, sill, range] % 注意sill 在这里指基台值实际拟合时基台值 nugget partial sill nugget params(1); sill params(2); % 这里的 sill 是完整基台值 rnge params(3); psill sill - nugget; % 偏基台值 switch lower(model) case spherical gamma_model_val zeros(size(h)); idx h rnge; gamma_model_val(idx) nugget psill * ... (1.5 * h(idx) / rnge - 0.5 * (h(idx) / rnge).^3); gamma_model_val(~idx) sill; case exponential gamma_model_val nugget psill * (1 - exp(-3 * h / rnge)); case gaussian gamma_model_val nugget psill * (1 - exp(-3 * (h / rnge).^2)); otherwise error(不支持的变异函数模型%s, model); end end然后是拟合函数。直接对params做fminsearch有个问题参数可能被优化成负值而块金值、变程在物理上必须非负。我的做法是把参数变换到对数域再优化这样无论优化器怎么迭代指数变换回来始终是正数function [nugget, sill, rng] fit_variogram(h, gamma_hat, model) % 拟合变异函数参数不依赖优化工具箱 % 参数变换到对数域保证优化过程中始终为正 nugget_init min(gamma_hat) * 0.6; % 初始块金值 sill_init max(gamma_hat); % 初始基台值 rng_init max(h) / 3; % 初始变程 % 目标函数预测值与实测值的均方误差 obj_fun (p) sum((gamma_hat - variogram_model(h, model, exp(p))).^2); p0 log([nugget_init, sill_init, rng_init]); p_opt fminsearch(obj_fun, p0, optimset(Display, off)); p_opt exp(p_opt); nugget p_opt(1); sill p_opt(2); rng p_opt(3); end3.3 普通克里金主函数与网格插值有了变异函数参数预测就水到渠成了。下面是普通克里金的主函数function [Z_pred, var_pred] ordinary_kriging(x_obs, y_obs, z_obs, x_pred, y_pred, model, params) % 普通克里金预测 % 输入 % x_obs, y_obs, z_obs - 观测点坐标与观测值 % x_pred, y_pred - 待预测点坐标可以是网格化后的向量 % model - 变异函数模型 % params - [nugget, sill, range] % 输出 % Z_pred - 预测值 % var_pred - 克里金方差 n length(z_obs); obs_coords [x_obs(:) y_obs(:)]; % 观测点之间的变异函数矩阵 D_obs pdist2(obs_coords, obs_coords); Gamma variogram_model(D_obs, model, params); % 构造克里金矩阵 A A [Gamma, ones(n, 1); ones(1, n), 0]; % 预分配输出 m length(x_pred(:)); Z_pred zeros(m, 1); var_pred zeros(m, 1); % 对待预测点逐个求解 for k 1:m % 待预测点到所有观测点的距离 d0 sqrt((x_obs(:) - x_pred(k)).^2 (y_obs(:) - y_pred(k)).^2); % 待预测点与观测点之间的变异函数向量 gamma0 variogram_model(d0, model, params); % 右端向量 b [gamma0; 1]; % 求解权重 lambda A \ b; % 预测值 Z_pred(k) sum(lambda(1:n) .* z_obs(:)); % 克里金方差 var_pred(k) sum(lambda(1:n) .* gamma0) lambda(end); end % 保持输出维度与输入一致 Z_pred reshape(Z_pred, size(x_pred)); var_pred reshape(var_pred, size(x_pred)); end主脚本调用示例% 生成一组模拟数据 rng(42); n 50; x_obs rand(n, 1) * 10; y_obs rand(n, 1) * 10; % 模拟一个真实的空间场正弦余弦噪声 z_obs sin(x_obs) cos(y_obs) 0.3 * randn(n, 1); % 第一步计算实验变异函数 [h, gamma_hat, n_pairs] experimental_variogram(x_obs, y_obs, z_obs, 15); % 第二步拟合理论模型 model spherical; [nugget, sill, rng] fit_variogram(h, gamma_hat, model); fprintf(拟合结果块金值%.3f基台值%.3f变程%.3f\n, nugget, sill, rng); % 第三步生成网格 [xx, yy] meshgrid(0:0.2:10, 0:0.2:10); [Z_pred, var_pred] ordinary_kriging(x_obs, y_obs, z_obs, xx, yy, model, [nugget, sill, rng]); % 可视化 subplot(1, 2, 1); contourf(xx, yy, Z_pred, 30); colorbar; title(Kriging预测值); hold on; plot(x_obs, y_obs, ko, MarkerSize, 4); subplot(1, 2, 2); contourf(xx, yy, var_pred, 30); colorbar; title(克里金方差);4. 实测样本点回代、与IDW对比、误差验证4.1 验证一预测值在样本点处是否等于实测值写完之后测试第一件事是回代验证用同一个数据集预测样本点本身的位置。普通克里金有一个数学上可以证明的性质在观测点处预测值等于实测值克里金方差等于0在无块金效应时。但实际上用含噪声的数据测试时会发现因为有块金值的存在方差不会刚好是0而是趋近于块金值。这是合理的不是代码bug恰恰说明噪声被正确量化了。我自己第一次跑回代时看到方差不是0还以为是矩阵求逆出了问题排查了半天才发现这本来就是对的。这个心得写在前面可以帮你少花一天时间。4.2 验证二与反距离加权插值法对比光看Kriging自己跑得顺不顺还不够得找个参照物对比。我拿Kriging和IDW反距离加权插值在同一组数据上做了对比。IDW的实现很简单function Z_idw idw_interp(x_obs, y_obs, z_obs, x_pred, y_pred, power) % 反距离加权插值 Z_idw zeros(size(x_pred)); for k 1:numel(x_pred) d sqrt((x_obs - x_pred(k)).^2 (y_obs - y_pred(k)).^2); d(d 0) 1e-10; w 1 ./ (d.^power); Z_idw(k) sum(w .* z_obs) / sum(w); end end对比结果中有一件事特别有意思当数据分布均匀、噪声比较小时IDW和Kriging的预测结果非常接近肉眼很难看出差别。但一旦数据有块金效应噪声大或者采样点分布不均IDW就开始露馅了——它会画出很多以采样点为中心的牛眼状斑块而Kriging的等值线明显更平滑也更符合真实的空间连续性。更深层的差异在于IDW的权重是人为设定的距离幂次通常取2它没有从数据本身学习任何东西而Kriging的权重是从变异函数里解出来的它知道每个方向、每个距离上数据到底有多相关。4.3 预测方差的现实意义有了方差输出你可以做一件常规插值永远做不到的事画置信区间图。假设Kriging预测误差服从正态分布在无偏高斯过程假设下成立那么(Z_pred \pm 1.96 \times \sqrt{var_pred})就是95%置信区间。实际应用中的价值非常大。比如环境监测布点时你可以先对已有数据做一次Kriging插值看方差分布图方差大的区域就是采样点稀疏的区域也就是应该追加采样的位置。这就是所谓的自适应采样策略。5. 调参与避坑这些细节才是Kriging能否落地的关键5.1 变异函数拟合别过度追求好看用fit_variogram拟合变异函数时我第一次拿到结果非常兴奋——拟合曲线把点连得几乎完美。但用这个完美的变异函数去做预测反而出现了很奇怪的预测结果远处出现莫名其妙的波动方差图上也出现了不规则的条带。后来我才反应过来这是过拟合了。实验变异函数的末尾几个点往往只靠少数几个点对支撑统计意义很弱强行让拟合曲线经过这些点会扭曲整个预测。正确的做法是计算实验变异函数时只保留点对数大于某个阈值的箱子比如至少20对。这样尾部那些不稳定的点就直接不参与拟合了。手动检查拟合曲线如果尾部明显被一两个异常点带偏可以考虑手动调整变程和基台值不要一味相信优化器。变异函数的三个参数是有物理解释的块金值不能超过整体方差的50%如果拟合出来块金值比数据方差还大基本说明数据本身就没有空间自相关性或者你的坐标轴搞错了单位。5.2 矩阵奇异与数值稳定性的处理克里金矩阵(A)的数值稳定性是个大问题尤其是当样本点之间存在非常近的距离时。两个点几乎重合变异函数矩阵对应的两行几乎相同矩阵就接近奇异了。A \ b不报错但会给出一个非常离谱的结果。我处理这个问题的经验有两条一是在计算变异函数矩阵时给对角线加一个很小的抖动值。具体做法是Gamma variogram_model(D_obs, model, params) 1e-10 * eye(n)。这个抖动在数值上不明显影响预测结果但能让矩阵从病态变为良态。二是用torst判断矩阵条件数如果cond(A) 1e10说明矩阵接近奇异需要检查数据。有时候是重复点太多有时候是变程设置得太小导致距离比较远的点对之间的协方差几乎为0。5.3 邻域选择全用还是只用最近的教科书上的克里金公式默认使用所有已知点。但实际数据量一大超过几百个点用所有点会带来两个问题一是矩阵规模大计算慢二是距离很远的点对预测的贡献其实很小却可能引入数值噪声。工程上常用的做法是邻域克里金预测每个点时只选择距离最近的3050个点参与求解。你可以在ordinary_kriging函数里加一个max_neighbors参数先按距离排序取最近的K个点构成矩阵。我实测下来的经验是数据比较平滑时30个邻居和全部点几百个的结果差异几乎看不出来但计算速度能快一个数量级。数据噪声大时邻域克里金甚至比全量克里金更稳。5.4 MATLAB环境里常见的几个翻车现场写代码过程中我遇到过几个特别典型的MATLAB问题写在这里给你避雷第一点乘和直接乘的教训。我在最早版本里写gamma_model_val nugget psill * (1.5 * h / rnge - 0.5 * (h / rnge).^3)时漏掉了最后那个.^3里的点号。结果h是向量时(h/rnge)^3尝试做矩阵乘法直接报错或者给出错误的结果。在MATLAB里凡是涉及数组元素级运算的一律用.*、./、.^。第二pdist2的输入必须是列向量或矩阵。我从数据表里读出来的坐标常常是行向量直接扔进pdist2会得到维度不匹配的报错。养成习惯进入函数前先x_obs x_obs(:)强制转成列向量。第三脚本里出现变量名gamma会覆盖MATLAB自带的gamma函数。我一开始顺手把变异函数值存成变量gamma后面再调用gamma(n)计算伽马函数时结果全乱了。这种变量遮蔽内置函数的问题在长时间运行的脚本里非常隐蔽建议用gamma_val或者gamma_hat这样的名字。第四注意variogram_model函数里的向量化写法。我用zeros(size(h))预分配了输出然后对h rnge的索引做赋值这样无论h是标量还是向量都能正确处理。h rnge的边界情况也要注意球形模型在h rnge时应该刚好等于基台值我用了idx h rnge把等号包含进去了这样保证连续性。第五如果你用的是比较老的MATLAB版本可能没有pdist2。有个不依赖统计工具箱的替代写法D sqrt((x - x).^2 (y - y).^2)这个写法用隐式展开实现任何版本都能跑数据量几千个点以内没有任何问题。6. 什么时候别用Kriging最后聊点实在的。Kriging不是万能的在我自己的项目里这些情况我会直接换方法样本量太少少于20个变异函数本身的估计就不可靠拟合出来的参数更是随缘不如老老实实用IDW。数据没有空间相关性如果拟合出来的变程很小、块金值占比很大说明数据本身几乎就是随机场Kriging和取平均值没什么区别。数据量极大几十万点全量克里金的矩阵是(n \times n)的几十万个点根本算不动这种情况应该用局部邻域克里金或者换高斯过程回归的近似算法。需要快速出图、精度要求不高Kriging的搭建成本比scatteredInterpolant高得多如果只是看一眼分布趋势没必要杀鸡用牛刀。多说我个人的体会Kriging算法在MATLAB里其实有现成的工具箱比如mGstat、DACE但自己从零写一遍的真正价值不在于省那几个工具钱而在于你写完一遍之后对空间插值的理解深度完全不一样了。以前看到空间自相关变程块金值这些词只觉得玄乎亲手用代码实现一遍之后它们全都变成了脑子里看得见摸得着的概念。最后再分享一个小技巧在跑网格插值之前先用随机抽样的方式做一次交叉验证——把样本随机分成两份一份训练一份测试反复做几次统计预测误差的均值和标准差。这个操作能让你在正式出图之前就发现参数设置的问题比等图出来之后发现不合理再回头改参数效率高太多。自己手写的Kriging代码配合这套验证流程在绝大多数空间插值场景里都够用了。本文还有配套的精品资源点击获取