企业官网建设流程全解析

简介：这是一套基于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 变异函数模型的拟合

拿到经验半方差散点图之后，需要用理论模型去拟合它，因为插值时要对任意距离求半方差，不能只靠离散的箱均值。常用的理论模型有三种：

这里的 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插值结果');

核心的变异函数计算我这里用computeVariogram：

function [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, numBins+1); hMean = zeros(numBins, 1); gammaMean = zeros(numBins, 1); nPairs = zeros(numBins, 1); for k = 1:numBins idx = (distVec >= edges(k)) & (distVec < edges(k+1)); 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(k+1)); 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时取基台值C0+C1。目标函数用点对数加权的最小二乘，点对数越多的箱子权重越大：

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]; % 参数约束：C0>=0, C1>=0, a>0 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 (n+1 × n+1) A = zeros(n+1, n+1); A(1:n, 1:n) = Gamma; A(n+1, 1:n) = 1; A(1:n, n+1) = 1; A(n+1, n+1) = 0; % 为避免矩阵奇异，加一个小扰动 A(1:n, 1:n) = A(1:n, 1:n) + 1e-10 * eye(n); % 预分解矩阵，加速多个点的插值（对同一组采样点只需分解一次） [L, U, p] = lu(A, 'vector'); b = zeros(n+1, 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(n+1, 1); sol(p) = U \ (L \ rhs); lambda = sol(1:n); mu = sol(n+1); % 插值估计 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在你脑子里就不是一个黑箱子了，后面遇到再复杂的问题也知道从哪个环节入手去改。

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