MATLAB手写Kriging算法：从变异函数到带方差的空间插值

简介：克里金（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_{i=1}^{N(h)} (z(x_i) - z(x_i+h))^2 ]

其中(h)是两个点之间的距离，(N(h))是距离落在(h)附近的所有点对数量，(z(x_i))是位置(x_i)处的观测值。除以2是有讲究的，这个"半"字就是这么来的。

这个式子其实很好理解，你可以把(z(x_i) - z(x_i+h))想象成"两个邻居的数值差异"，差异越大，说明这个距离上数值变化越剧烈，变异函数值就越大。把不同距离上的变异函数值画出来，就会看到一条通常单调递增的曲线：距离越远，差异越大，直到一个平台——超过这个距离之后，两个点的数值基本就互不相干了。

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_{i=1}^{n} \lambda_i z(x_i) ]

关键问题是权重(\lambda_i)怎么确定。Kriging的思路是：让预测方差最小，同时保证无偏。

经过拉格朗日乘子法推导，最终得到的是一个( (n+1) \times (n+1) )的线性方程组：

[ \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时最常用的是这三个：

公式里的三个参数请你务必记牢，它们是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是元素级别的幂运算，如果写成^ 2，MATLAB会试图做矩阵乘法，结果完全不一样。二是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); end

3.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 邻域选择：全用还是只用最近的？

教科书上的克里金公式默认使用所有已知点。但实际数据量一大（超过几百个点），用所有点会带来两个问题：一是矩阵规模大，计算慢；二是距离很远的点对预测的贡献其实很小，却可能引入数值噪声。

工程上常用的做法是邻域克里金：预测每个点时，只选择距离最近的30～50个点参与求解。你可以在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代码，配合这套验证流程，在绝大多数空间插值场景里都够用了。

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