简介：这份资源是面向地质、环境、气象及遥感领域科研人员与工程师的克里金插值工具包，基于DACE工具箱实现简单克里金、普通克里金与泛克里金等算法，用于在观测数据不完整时估计连续变量的空间分布。压缩包共28个文件，约1.5MB，以m脚本文件为主，辅以asv备份、mat示例数据、pdf说明文档、changelog更新日志和txt数据文件，覆盖变异函数建模、回归多项式、相关函数与预测器等核心模块，并附协同克里金计算示例。已有1890人学习下载。读者可借助已注释的源码理解最佳线性无偏估计、变程与nugget效应等参数含义，结合示例数据完成数据导入、变异函数拟合、插值计算与结果可视化，也可在此基础上进行二次开发与算法扩展，适合作为空间插值入门与工程应用的参考工具。

1. 克里金插值 dace 工具箱：从一组散点到一张可信曲面

做地质建模、代理模型优化或者实验设计的人，迟早会撞上同一个问题：手头只有十几二十个采样点，却要画出一张能支撑后续决策的连续曲面。线性插值太糙，多项式拟合又容易在点稀疏的地方甩出离谱的振荡。这时候克里金插值就该上场了——它不把未知点当成孤立的坐标，而是用变异函数描述空间相关性，再通过加权求和给出预测，顺带还能吐出每个预测点的方差。dace 工具箱就是干这件事的 MATLAB 老牌实现，全称 Design and Analysis of Computer Experiments，核心是dacefit和predictor两个函数，配合regpoly与corrgauss这类回归和相关的组合，能把一组散点拟合成带置信区间的响应面。它适合谁？做仿真代理模型、地质储量估算、水文参数反演，以及任何需要“少样本、带不确定性”插值的场景。下面我按自己拆包复现的顺序，把这份资源从安装到调参到排错讲透。

2. dace 工具箱的目录结构与核心函数拆解

拿到一份 dace 资源包，第一件事不是急着跑 demo，而是先看清它到底给了什么。dace 的代码组织方式很朴素，没有花哨的工程封装，所有函数平铺在根目录下，靠文件名区分职责。理解这套命名逻辑，后面调参和改核函数才不会抓瞎。

2.1 文件清单与职责划分

一个完整的 dace 包通常包含以下文件，我按调用频率排了序：

这里要强调一点：dace 的回归部分和自相关部分是解耦的。regpoly决定的是全局趋势，corr决定的是局部偏差。很多人一上来就换核函数，却忽略了趋势项选错才是预测跑偏的主因。

2.2 拟合与预测的调用链路

dace 的使用逻辑是“先拟合出模型结构体，再拿结构体去预测”，两步分离。下面这段是我常用的最小可运行骨架：

% 假设已有样本点 S（n×d）和响应 Y（n×1） % S 每行是一个 d 维样本，Y 是对应输出 % 第一步：拟合模型 % regpoly0 表示常数趋势，corrgauss 表示高斯核 % theta 是核宽度初值，lob 和 upb 是优化下上界 [model, perf] = dacefit(S, Y, @regpoly0, @corrgauss, ... [1.0], [1e-3], [1e3]); % 第二步：在待测网格上预测 % X 是 m×d 的待测点矩阵 [Yp, MSE] = predictor(X, model); % Yp 是预测均值，MSE 是预测方差 % 画曲面时常用 Yp，画置信带时用 sqrt(MSE)

逻辑说明：dacefit内部做的是极大似然估计，它会在lob和upb之间搜索最优的theta，让相关矩阵的行列式与残差平方和的组合最优。perf返回的是优化过程信息，包含最终似然值，调参时盯这个值。predictor拿到模型后，对每个待测点计算它与所有样本点的相关向量，再套克里金公式给出均值和方差。

参数说明：theta初值给 1.0 是保守做法，如果样本点坐标跨度大（比如经纬度量级），初值要给到 0.01 甚至更小，否则优化器第一步就跨出边界。lob和upb别设太窄，我一般用[1e-3, 1e3]，让优化器自己找。regpoly0适合响应整体平稳的情况，如果响应有明显倾斜趋势，换regpoly1。

2.3 相关函数的选择依据

核函数不是随便挑的，它直接决定曲面的光滑程度。高斯核产生的曲面无限可微，适合物理场连续变化的场景；指数核在原点不可微，曲面会有折痕感，适合存在突变界面的地质层；球状核介于两者之间，在变差函数拟合里最常见。判断方法很简单：把样本点按距离排序，算半变差，看散点图的形状。如果半变差在短距离内快速上升然后平缓，高斯核；如果线性上升很久才平缓，指数核；如果先抛物线后平缓，球状核。选错核函数不会报错，但预测曲面会在样本点之间出现不自然的鼓包或塌陷，这是最隐蔽的坑。

3. 从散点到曲面：完整建模流程与参数调优

这一章把上一章的骨架填成能跑通的完整流程。我按“数据准备 → 趋势项选择 → 核参数优化 → 预测与验证”四步走，每步都给可抄的代码和判断标准。

3.1 样本点预处理与重复点合并

原始数据里经常有坐标完全相同的点，或者距离近到相关矩阵接近奇异。dace 的相关矩阵求逆对重复点极其敏感，一旦出现，dacefit会报Matrix is singular或者给出 NaN。处理办法是先用dsmerge合并：

% S 是原始样本坐标，Y 是响应 % dsmerge 会把距离小于 tol 的点合并，响应取均值 tol = 1e-6; % 合并阈值，按坐标量级调整 [S_clean, Y_clean] = dsmerge(S, Y, tol); % 检查合并后样本量 fprintf('原始点数 %d，合并后 %d\n', size(S,1), size(S_clean,1));

逻辑说明：dsmerge内部用距离判据找近邻对，把它们的响应做算术平均。tol的取值要看坐标单位，如果坐标是米级，1e-6足够；如果是经纬度，1e-6可能太小，要放到1e-4。合并后如果样本量掉得太多（超过 20%），说明原始采样有问题，得回头查数据采集环节。

参数说明：合并阈值不是越小越好。设太小起不到作用，设太大会把有意义的近邻点也吞掉，导致局部特征丢失。我的经验是取样本点平均最近邻距离的十分之一作为tol。

3.2 趋势项与核函数的组合策略

趋势项和核函数的组合决定了模型的表达能力。dace 允许任意regpoly配任意corr，但实际可用的组合就那么几种。下面这张表是我踩过坑之后总结的对应关系：

选组合的方法：先画响应随各维度的散点图，看有没有明显趋势。有趋势就升regpoly阶数，没趋势就用regpoly0。然后看样本点之间的响应变化是平滑还是跳跃，平滑用高斯，跳跃用指数。别一上来就上regpoly2配高斯，参数多、优化慢，还容易过拟合。

3.3 theta 初值与优化边界设置

theta是 dace 里最玄学的参数，它控制相关函数的作用范围。theta大，相关衰减慢，曲面光滑；theta小，相关衰减快，曲面贴近样本点。dacefit会用优化器搜theta，但初值和边界设不好，优化器会陷在局部最优里。

% 多维度时 theta 是向量，每个维度一个值 d = size(S_clean, 2); % 维度 % 初值策略：按坐标跨度的倒数给 span = max(S_clean) - min(S_clean); theta0 = 1 ./ span; % 每个维度初值不同 % 边界：下界给初值的 1e-2，上界给 1e2 lob = theta0 * 1e-2; upb = theta0 * 1e2; [model, perf] = dacefit(S_clean, Y_clean, @regpoly0, @corrgauss, ... theta0, lob, upb);

逻辑说明：theta0按坐标跨度倒数给，是因为相关函数里theta乘的是距离，跨度大的维度需要小的theta才能让相关在合理距离内衰减。如果所有维度都给同一个初值，跨度大的维度会主导优化，跨度小的维度被忽略。lob和upb给初值的两个数量级范围，既给优化器空间，又不至于跑到无意义的区域。

参数说明：perf里的似然值可以用来判断优化是否收敛。如果perf返回的迭代次数等于最大迭代数，说明没收敛，得放宽边界或换初值。我一般会跑两三次，每次换不同初值，取似然最大的那个模型。

3.4 预测结果的不确定性解读

dace 的predictor除了给预测均值，还给均方误差。这个 MSE 是克里金插值的核心价值——它告诉你每个预测点有多可信。样本点附近的 MSE 接近零，远离样本点的区域 MSE 上升。画图时把sqrt(MSE)作为置信带宽度叠上去，能直观看出哪些区域的预测是“猜”的。

% 生成待测网格 [X1, X2] = meshgrid(linspace(0, 10, 50), linspace(0, 10, 50)); X = [X1(:), X2(:)]; % 预测 [Yp, MSE] = predictor(X, model); % 重塑成网格形状 Yp_grid = reshape(Yp, 50, 50); MSE_grid = reshape(MSE, 50, 50); % 画曲面和置信带 figure; surf(X1, X2, Yp_grid); hold on; % 置信带用上下界表示 surf(X1, X2, Yp_grid + 1.96*sqrt(MSE_grid), 'FaceAlpha', 0.3); surf(X1, X2, Yp_grid - 1.96*sqrt(MSE_grid), 'FaceAlpha', 0.3);

逻辑说明：1.96是 95% 置信水平的正态分位数。FaceAlpha调透明度，让置信带不遮挡主曲面。如果 MSE 在某个区域特别大，说明那里样本点太稀，要么补采样，要么在报告里明确标注该区域不可靠。

参数说明：meshgrid的分辨率按需求定，50×50 适合快速查看，200×200 适合出图。分辨率太高会让predictor变慢，因为每个点都要算与所有样本点的相关向量。

4. 避坑与排查：dace 使用中的五个血泪教训

dace 是个老工具箱，代码风格停留在早期 MATLAB 时代，很多错误不会给出友好提示，直接甩一个矩阵奇异或者 NaN。下面五条是我实际踩过的坑，按“现象 → 原因 → 解决”整理。

4.1 相关矩阵奇异导致拟合失败

现象：dacefit报Warning: Matrix is close to singular，或者返回的model里theta是 NaN。

原因：样本点里有重复坐标，或者两个点距离近到相关函数值接近 1，相关矩阵的行列式趋近零，求逆失败。

解决：先用dsmerge合并重复点，合并阈值按坐标量级取。如果合并后仍报奇异，检查是不是有维度坐标完全不变（比如所有点的某一维都是同一个值），这种情况该维度对相关矩阵没有贡献，要么删掉该维度，要么给该维度加极小扰动。

4.2 预测值在样本点处不闭合

现象：用predictor在样本点位置预测，得到的值跟原始响应对不上，差得还挺多。

原因：dacefit的优化没收敛，theta停在边界上，导致相关矩阵近似单位阵，模型退化成纯回归。这时候预测值只反映趋势项，不反映局部偏差。

解决：检查perf的迭代信息，如果theta顶在lob或upb上，说明边界设窄了。放宽边界，或者换初值重跑。另一个可能是趋势项阶数太高，把局部偏差吸收了，换成regpoly0试试。

4.3 高维输入时优化极慢

现象：输入维度超过 5 维后，dacefit跑十几分钟不出结果，CPU 单核跑满。

原因：dace 的似然优化是数值优化，每步都要算相关矩阵的行列式和逆，复杂度是 O(n³)，n 是样本量。维度高时theta向量长，优化空间大，迭代次数暴涨。

解决：降维。先用主成分分析或相关性分析筛掉不重要的维度，把维度压到 3 到 5 维再喂给 dace。如果必须用高维，把theta的优化拆成逐维优化，或者用fmincon替代默认优化器，设最大迭代数上限。

4.4 外推区域预测值爆炸

现象：在样本点覆盖范围之外预测，Yp出现极大或极小的离谱值。

原因：克里金插值本质是加权平均，外推时权重不再受样本点约束，趋势项主导，如果趋势项是regpoly1或regpoly2，外推会沿趋势无限延伸。

解决：外推区域不要用高阶趋势项，换regpoly0。如果必须外推，把预测范围限制在样本凸包附近，超出部分标注为不可信。更稳妥的做法是在外推区域补采样，而不是硬预测。

4.5 不同核函数结果差异巨大

现象：同一组数据，换核函数后预测曲面形状完全变了，不知道信哪个。

原因：核函数决定了相关矩阵的结构，不同核对应的似然曲面不同，优化到的theta也不同。高斯核倾向光滑解，指数核倾向粗糙解，两者在样本稀疏区差异最大。

解决：用交叉验证选核。留一法：每次去掉一个样本点，用其余点拟合，预测被去掉的点，算预测误差。换不同核函数各跑一遍，选误差最小的。别凭感觉选，数据说话。

5. 进阶技巧：用交叉验证和变差函数手动校准模型

dace 自带的优化器在样本量小于 30 时经常给出不稳定的theta，这时候手动校准比自动优化靠谱。我一般走两条路：一是留一交叉验证，二是手动拟合变差函数。

留一交叉验证的代码骨架：

n = size(S_clean, 1); err = zeros(n, 1); for i = 1:n % 留出第 i 个点 idx = true(n, 1); idx(i) = false; S_train = S_clean(idx, :); Y_train = Y_clean(idx); % 拟合 model_i = dacefit(S_train, Y_train, @regpoly0, @corrgauss, ... theta0, lob, upb); % 预测被留出的点 Y_pred = predictor(S_clean(i, :), model_i); err(i) = Y_clean(i) - Y_pred; end % 计算均方误差 cv_mse = mean(err.^2); fprintf('留一交叉验证 MSE: %.6f\n', cv_mse);

逻辑说明：每次留一个点，用剩余点建模，预测留出点，累积误差。cv_mse越小，模型泛化越好。这个流程可以套在不同核函数和不同theta初值上，做网格搜索。缺点是慢，n 个点要拟合 n 次，但样本量小于 50 时完全可接受。

参数说明：theta0、lob、upb在循环外定义一次即可，循环内复用。如果某次拟合报奇异，跳过该点并在最后统计有效次数。

变差函数手动校准的思路：把样本点按距离分组，算每组的半变差，画散点图，然后用球状模型或指数模型拟合这条散点曲线，从拟合参数反推theta。dace 的corrspherical和correxp里的theta与变差函数的变程有对应关系，具体换算要看核函数的定义式。这一步需要点耐心，但校准出来的theta比自动优化稳得多，尤其在地质数据上。

从那以后我每次拿到新数据，都强制走一遍留一交叉验证，再决定用哪个核函数和theta。自动优化给的是起点，交叉验证给的是终点，两者结合才敢把曲面拿去做决策。希望帮到你。

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