简介克里金(Kriging)插值法又称空间自协方差最佳插值法是以矿业工程师D.G.Krige命名的一种最优内插方法广泛用于地下水模拟、土壤制图等地质统计领域。这套MATLAB实现共包含19个文件以16个m脚本为主涵盖dacefit回归拟合、predictor预测、corrgauss相关模型、gridsamp网格生成等核心函数并配有changelog更新说明、PDF理论文档和示例数据mat文件压缩包仅1.48MB轻量易用。目前已有4976人学习下载。借助这套代码读者可以快速搭建克里金插值流程理解DACE工具箱的建模与预测逻辑通过PDF文档梳理算法原理结合mat数据直接运行验证m函数模块化设计也便于按需修改与二次开发适合从事空间插值研究或工程应用的地学、测绘及数据科学人员参考学习。 第一次在MATLAB里实现Kriging插值我差点被一个矩阵奇异问题劝退。Kriging这个算法在空间统计学里已经存在了几十年核心思路可以概括为一句话用已知采样点的加权平均去估计未知位置的值。但跟普通插值方法不同的是它的权重不完全由距离决定而是由一个从数据中估计出来的空间变异函数决定所以它能给出“最优无偏估计”还带着一个可用的预测方差这是IDW这类方法做不到的。正因如此Kriging在气象、水文、地质、土壤调查、环境监测这些需要空间预测的领域基本是标配算法。而MATLAB因为矩阵运算方便、可视化工具完整成了实现Kriging最顺手的语言之一。这篇文章我会从Kriging的数学原理讲起然后给出一套完整的普通克里金MATLAB代码包括实验变异函数计算、理论模型拟合、预测主函数、网格化出图和交叉验证最后把我在调参和踩坑过程中遇到的高频问题一并整理出来。文中代码我尽量只依赖MATLAB基础函数不用额外工具箱也能跑通。1. Kriging到底在算什么从空间自相关说起1.1 空间自相关的数学刻画变异函数做空间插值的人都知道地理学第一定律距离越近的事物越相似。IDW反距离加权的做法很简单直接按距离的倒数分配权重距离越近权重越大。但IDW有个问题它没有考虑样本点之间的相互位置关系如果两个采样点离得很近它们提供的信息其实是重复的IDW却会把它们当两条独立信息用。Kriging的思路聪明在它先通过变异函数也叫半变异函数把数据的“空间结构”提取出来。变异函数定义长这样γ(h) 0.5 * E[(Z(xh) - Z(x))²]简单理解就是把所有距离为h的点对找出来求它们观测值差值的平方平均再除以2。这个值越小说明这个尺度上空间连续性越好h越大γ(h)通常也越大说明距离远了相似性就下降。这个曲线就是Kriging用来分配权重的依据它决定了预测结果在不同方向、不同尺度上的平滑程度。实际计算时不会真的去求数学期望而是把点对按距离分组比如0到0.05一个箱子、0.05到0.10一个箱子这样分下来每个箱子里的点对半方差取平均得到的散点就是“实验变异函数”。后续要做的是用一个理论模型去拟合这些散点。1.2 普通克里金、简单克里金、泛克里金选哪个Kriging不是一个单一算法而是一个家族。简单克里金假设研究区域的均值已知且恒定比如区域平均气温已知那插值就围绕这个已知均值展开普通克里金假设均值未知但恒定这是绝大多数实际场景下的合理假定泛克里金则处理均值本身存在趋势的情况比如污染物浓度从污染源向外明显衰减时需要先把趋势项剥离出来再插值。三者里普通克里金应用最广原因是它不需要额外输入太多先验知识只需要采样点位置和观测值就能跑。实际项目里我默认先上普通克里金效果不够再考虑泛克里金。如果你直接把简单克里金的均值填错了结果往往还不如普通克里金稳。所以这篇博文的所有代码我全部围绕普通克里金展开。1.3 为什么用MATLAB实现现在实现Kriging的语言选择很多Python有pykrige、GSTools这些现成库R有gstat但MATLAB仍然有一批忠实用户。原因很直接克里金方程组本质上是求解一个线性系统Aλ b这是MATLAB最擅长的运算几条代码就能搞定。再加上MATLAB的surf、contourf、scatter3一套可视化组合拳非常适合快速验证插值效果。不过有一点要提醒网上很多MATLAB版Kriging代码用了pdist2、knnsearch这类函数它们是Statistics and Machine Learning Toolbox里的不是基础MATLAB自带。如果只有基础环境跑起来会直接报错。我这篇文章里的代码尽量避开工具箱函数距离矩阵用手写矩阵运算完成保证绝大多数读者能直接跑。2. 动手编码前先吃透三个数学概念2.1 变异函数的三个关键参数块金值、基台值、变程拟合变异函数时你一定会碰到三个参数块金值、基台值、变程。把它们理解清楚代码才能调得动。块金值nugget描述的是h趋近于0时变异函数仍然存在的跃变值。它来源于测量误差、采样定位误差或者微小尺度上的随机变异。用生活场景类比你问邻居房价两个房子几乎挨着但报价就是会有差异这个差异就是“块金”。基台值sill是变异函数达到的平稳上限可以理解为空间方差的总水平。变程range是从原点开始到变异函数不再明显增加的距离超过变程之后采样点之间已经不存在空间相关性Kriging权重会趋于均匀。拟合球状模型时代码里要估计三个参数块金值(c0)、偏基台值(c也就是基台值减去块金值)、变程(a)。初值的选取很关键我一般先画出实验变异函数散点图目测一下块金值、基台值和变程的大致范围再丢给优化器去精调直接给随机初值很容易陷入局部最优。2.2 球状、指数、高斯变异函数模型怎么选理论变异函数模型有很多种最常用的三个是球状模型、指数模型和高斯模型。它们的表达式和适用场景差别很大我整理了一个对照表模型数学形式h为距离特点适用场景球状γ(h) c0 c(1.5h/a - 0.5(h/a)³)h≤aγ(h)c0cha在变程处严格到达基台曲线有明确拐点地质储量、土壤属性最常用最稳妥指数γ(h) c0 c(1 - exp(-h/a))渐近趋于基台有效变程约为3a污染物扩散、气象要素过渡平滑高斯γ(h) c0 c(1 - exp(-(h/a)²))原点附近非常平滑连续性强高程、水位等高连续性场景但易数值病态实际编程时球状模型最容易处理因为它在h超过a之后直接等于基台值行为很清晰。我大多数项目里第一选择就是球状模型如果交叉验证误差不满意再换指数模型对比。高斯模型虽然听起来高级但它的变异函数在原点附近太平滑导致克里金矩阵严重病态新手容易在这里翻车后面我会专门讲这个问题。2.3 克里金方程组无偏与最优怎么变成线性代数普通克里金的目标是找一组权重λ_i使得预测值等于Σλ_i * Z(x_i)同时满足无偏和预测方差最小两个约束。无偏条件很简单权重之和必须等于1这样才能保证预测值和真实值在均值水平上一致。最小方差条件则需要用拉格朗日乘子法引入一个中间变量μ。最终得到的线性方程组是对每个已知点i都有Σλ_j * γ(x_i - x_j) μ γ(x_i - x_0)再加上Σλ_i 1这行约束。写成矩阵形式就是Aλ bA矩阵是n1阶方阵左上角n×n块是已知点两两之间的变异函数值第n1列和n1行是1右下角是0。b向量前n项是预测点到每个已知点的变异函数值最后一项是1。λ向量前n项就是权重最后一项是拉格朗日乘子。在MATLAB里这个方程就一句话lambda A \ b。不是lambda inv(A)*b直接左除就可以了数值稳定性更好这是MATLAB最常用也是我强烈推荐的做法。3. Kriging的MATLAB完整实现3.1 构造模拟数据用随机采样点验证整个流程要验证Kriging代码最好先有一份已知答案的数据。我在这里构造一个二维正弦场作为“真实分布”然后随机取30个采样点并叠加噪声模拟现实中的观测误差。% 设定随机种子保证结果可复现 rng(42); % 生成30个随机采样点坐标范围 [0,1] x [0,1] n 30; dataXY rand(n, 2); x dataXY(:, 1); y dataXY(:, 2); % 用二维正弦场叠加噪声构造观测值 z sin(2 * pi * x) .* cos(2 * pi * y) 0.3 * randn(n, 1);这里有一个MATLAB新手常踩的坑sin(2*pi*x)里面用的.*是点乘因为x是列向量如果写成sin(2*pi*x) * cos(2*pi*y)MATLAB会试图做矩阵乘法直接报维度不匹配。这类逐元素运算一定要用点乘后面计算距离矩阵时也一样。3.2 求实验变异函数并拟合球状模型有了采样点之后第一步是计算点对之间的距离矩阵和半方差矩阵。% 距离矩阵第(i,j)个元素是第i个点和第j个点的欧氏距离 distMat sqrt((x - x).^2 (y - y).^2); % 半方差矩阵0.5 * (z_i - z_j)^2 gammaMat 0.5 * (z - z).^2; % 按距离分组计算实验变异函数 maxLag 0.6; % 最大滞后距离约为坐标范围的一半 nBins 12; % 分组数量 binEdges linspace(0, maxLag, nBins 1); hExp zeros(nBins, 1); gammaExp zeros(nBins, 1); for k 1:nBins % 去掉自己也去掉两端只取(h,k]区间的点对 idx find(distMat binEdges(k) distMat binEdges(k 1)); if ~isempty(idx) hExp(k) mean(distMat(idx)); gammaExp(k) mean(gammaMat(idx)); else hExp(k) NaN; gammaExp(k) NaN; end end % 删除没有点对的距离区间 valid ~isnan(hExp); hExp hExp(valid); gammaExp gammaExp(valid);% 球状模型函数h0时赋00ha用公式ha取基台值 sphVariogram (h, c0, c, a) ... c0 c * (1.5 * h / a - 0.5 * (h / a).^3) .* (h 0 h a) ... (c0 c) * (h a); % 用fminsearch做最小二乘拟合初值通过目测散点估计 objFun (p) sum((sphVariogram(hExp, p(1), p(2), p(3)) - gammaExp).^2); p0 [0.05, 1, 0.4]; pOpt fminsearch(objFun, p0); c0 pOpt(1); c pOpt(2); a pOpt(3); fprintf(拟合结果块金%.4f, 偏基台值%.4f, 变程%.4f\n, c0, c, a);这里最需要留意的是球状模型在h0处的定义。理论球状模型在原点处取0但公式本身会算出c0所以我在匿名函数里强制做了h 0的判断把原点位置单独处理成0。如果这里不加判断之后构建克里金矩阵时对角线上会出现一个c0直接污染整个求解结果。这个细节代码里体现得特别清楚。3.3 普通克里金预测核心函数预测函数是整段代码的心脏。它接收已知点坐标、观测值、待预测点坐标和变异函数参数返回预测值。function zPred ordinaryKriging(dataXY, dataZ, p0, c0, c, a) n length(dataZ); % 已知点两两之间的距离矩阵 distMat sqrt((dataXY(:,1) - dataXY(:,1)).^2 ... (dataXY(:,2) - dataXY(:,2)).^2); % 构建A矩阵左上角块 A11 c0 c * (1.5 * distMat / a - 0.5 * (distMat / a).^3) .* (distMat 0 distMat a) ... (c0 c) * (distMat a); A11(distMat 0) 0; % 对角线归零 % 拼装完整克里金矩阵 A [A11, ones(n, 1); ones(1, n), 0]; % 预测点与已知点之间的距离向量 distPred sqrt(sum((dataXY - p0).^2, 2)); % 构建b向量 bVec c0 c * (1.5 * distPred / a - 0.5 * (distPred / a).^3) .* (distPred 0 distPred a) ... (c0 c) * (distPred a); bVec [bVec; 1]; % 求解克里金方程组 lambda A \ bVec; % 预测值前n个权重与观测值加权求和最后一个元素是拉格朗日乘子不参与预测 zPred lambda(1:n) * dataZ; end有几个细节值得讲。第一A11(distMat 0) 0这行很关键因为对角线上是点和自己的距离半方差理应为0而我们的匿名函数对h0返回的是c0必须人工覆盖。第二b向量的最后一行是1这与A矩阵最后一行的约束对应用于求解过程中强制权重和为1。第三最后预测时只取lambda的前n个元素最后的拉格朗日乘子虽然参与了求解但不参与加权求和当初我第一次写代码时直接拿lambda全部分量去乘dataZ结果预测值偏得离谱。3.4 网格化批量预测与结果可视化现在把预测函数应用到网格上生成连续的插值曲面。% 生成60x60的预测网格 [Xg, Yg] meshgrid(linspace(0, 1, 60), linspace(0, 1, 60)); Zg zeros(size(Xg)); % 逐点预测 for i 1:numel(Xg) Zg(i) ordinaryKriging(dataXY, z, [Xg(i), Yg(i)], c0, c, a); end % 可视化采样点和预测面叠在一起看 figure; surf(Xg, Yg, Zg, EdgeColor, none); hold on; scatter3(x, y, z, 60, z, filled); colormap(turbo); % 比默认jet的颜色过渡更自然 colorbar; xlabel(X); ylabel(Y); zlabel(Z); title(Kriging插值结果);跑完这段代码你会看到预测曲面在采样点附近自然弯曲、穿过观测值在数据稀疏区域逐渐回归到整个区域的均值水平。这是Kriging的正常行为也解释了为什么Kriging在很多场景下比IDW稳健——它不会在远离数据的地方给出一个突兀的极端值。如果你用scatter3把采样点叠上去还能直观检查预测面是否在采样点处“穿过了”真实观测值。3.5 留一交叉验证评估插值质量插值做完只是第一步评估插值好坏的常规手段是留一交叉验证。做法很简单每次拿掉一个采样点用剩下的n-1个点在该位置做预测再和真实观测值对比最后统计RMSE和MAE。n length(z); predCV zeros(n, 1); for i 1:n idx true(n, 1); idx(i) false; xv dataXY(idx, :); zv z(idx); % 交叉验证时为了速度直接用全样本拟合好的变异函数参数 predCV(i) ordinaryKriging(xv, zv, dataXY(i, :), c0, c, a); end RMSE sqrt(mean((z - predCV).^2)); MAE mean(abs(z - predCV)); fprintf(RMSE%.4f, MAE%.4f\n, RMSE, MAE);交叉验证的严格做法是每一折都重新拟合一次变异函数参数因为理论上去掉一个点后空间结构会发生变化。但在点数量不大的情况下全样本拟合的参数已经比较稳定重新拟合性价比不高。如果你想严谨可以在循环体里把fminsearch步骤加进去只是运行时间会长不少。4. 调试经验与性能优化速查4.1 矩阵奇异、NaN结果是最大拦路虎我在文章开头提到的矩阵奇异问题是新手最容易遇到的第一道坎。常见原因有三个采样点之间存在重复坐标或距离极近的点导致A矩阵两行几乎线性相关。变程a估计过小使得大距离范围内的半方差都落到基台值附近矩阵对角线外的元素差异极小接近秩亏。使用了高斯模型它在原点附近太平滑构造出的矩阵条件数可以轻松到达1e17甚至更高。对策也很直接。第一数据预处理阶段做去重对距离小于某个阈值的点合并或剔除。第二在A矩阵对角线上加一个极小正则化量比如1e-10相当于给求解过程加了一个数值阻尼不影响最终权重但能显著改善条件数。第三如果使用高斯模型导致异常换回球状模型或者指数模型多数情况下能直接解决。A A 1e-10 * eye(size(A));这行代码建议加在求解之前是成本最低的保险。我当时第一次遇到NaN结果时排查了很久最后发现就是高斯模型在作祟换成球状模型后问题消失。4.2 权重为负、预测值出现跳变怎么处理Kriging的一个特性是某些权重可能是负的。这不是bug是数学上允许的结果相当于用周围点的振荡去校正局部趋势。但如果负权重过多或者绝对值过大预测值会莫名跳出样本值范围出现不合理的负浓度、负高度等。遇到这种情况第一件事是检查打印出来的lambda向量。如果某个权重绝对值超过0.5甚至到了1以上说明变异函数模型和参数有明显问题。常见修复路径包括适当增大块金值c0让随机噪声承担更多方差权重会变得更均匀。换用指数模型它的曲线更平滑不容易产生大幅度负权重。检查变程是否过小如果a只有坐标范围的十分之一权重很容易振荡。另外一个容易被忽略的点是预测点位置。Kriging在数据范围内部表现正常但一旦跑到采样点凸包之外外推区域容易出现明显的振荡和偏差。如果有外推需求务必缩小网格范围或者明确告知读者结果仅用于可视化参考。4.3 数据量变大后如何优化性能基础版代码在数据点少几十个的情况下完全够用。但如果采样点超过1000预测网格到几百×几百循环逐点预测的耗时就会猛涨。原因在于每个预测点都要构建一个n1阶线性系统并求解复杂度是O(n³)。优化方向有两个。第一改成局部克里金每个预测点只选择最近的m个采样点参与求解m通常取20到50。这样每次求解的系统都很小同时也符合空间相关性的实际规律——变程以外的点本来就没有多少信息贡献。第二预计算距离矩阵和变异函数矩阵避免重复计算。如果采样点布局固定可以在循环开始前把A矩阵中所有点对半方差算好每个预测点只需更新b向量能省掉大量冗余运算。如果需要处理几万点以上的数据或者网格特别密那就别硬扛MATLAB循环了考虑降维思路先做规则网格的粗插值再用细网格的局部克里金只对局部残差做修正这样可以把计算量降几个数量级。对于大多数项目局部克里金这一个优化已经足够。我个人现在写Kriging代码已经形成固定习惯拿到数据先画实验变异函数散点图目测参数范围再拟合跑完预测必做交叉验证求解前无条件加正则化项模型首选球状再看误差指标决定是否换指数。这套流程在几十到几百个点的项目里非常稳很少再被数值问题卡住。如果你也被Kriging的代码调参折磨过希望这份完整的MATLAB实现能帮你省下至少一个下午的排查时间。本文还有配套的精品资源点击获取 MATLAB实现普通克里金插值：变异函数与完整代码解析
简介克里金(Kriging)插值法又称空间自协方差最佳插值法是以矿业工程师D.G.Krige命名的一种最优内插方法广泛用于地下水模拟、土壤制图等地质统计领域。这套MATLAB实现共包含19个文件以16个m脚本为主涵盖dacefit回归拟合、predictor预测、corrgauss相关模型、gridsamp网格生成等核心函数并配有changelog更新说明、PDF理论文档和示例数据mat文件压缩包仅1.48MB轻量易用。目前已有4976人学习下载。借助这套代码读者可以快速搭建克里金插值流程理解DACE工具箱的建模与预测逻辑通过PDF文档梳理算法原理结合mat数据直接运行验证m函数模块化设计也便于按需修改与二次开发适合从事空间插值研究或工程应用的地学、测绘及数据科学人员参考学习。 第一次在MATLAB里实现Kriging插值我差点被一个矩阵奇异问题劝退。Kriging这个算法在空间统计学里已经存在了几十年核心思路可以概括为一句话用已知采样点的加权平均去估计未知位置的值。但跟普通插值方法不同的是它的权重不完全由距离决定而是由一个从数据中估计出来的空间变异函数决定所以它能给出“最优无偏估计”还带着一个可用的预测方差这是IDW这类方法做不到的。正因如此Kriging在气象、水文、地质、土壤调查、环境监测这些需要空间预测的领域基本是标配算法。而MATLAB因为矩阵运算方便、可视化工具完整成了实现Kriging最顺手的语言之一。这篇文章我会从Kriging的数学原理讲起然后给出一套完整的普通克里金MATLAB代码包括实验变异函数计算、理论模型拟合、预测主函数、网格化出图和交叉验证最后把我在调参和踩坑过程中遇到的高频问题一并整理出来。文中代码我尽量只依赖MATLAB基础函数不用额外工具箱也能跑通。1. Kriging到底在算什么从空间自相关说起1.1 空间自相关的数学刻画变异函数做空间插值的人都知道地理学第一定律距离越近的事物越相似。IDW反距离加权的做法很简单直接按距离的倒数分配权重距离越近权重越大。但IDW有个问题它没有考虑样本点之间的相互位置关系如果两个采样点离得很近它们提供的信息其实是重复的IDW却会把它们当两条独立信息用。Kriging的思路聪明在它先通过变异函数也叫半变异函数把数据的“空间结构”提取出来。变异函数定义长这样γ(h) 0.5 * E[(Z(xh) - Z(x))²]简单理解就是把所有距离为h的点对找出来求它们观测值差值的平方平均再除以2。这个值越小说明这个尺度上空间连续性越好h越大γ(h)通常也越大说明距离远了相似性就下降。这个曲线就是Kriging用来分配权重的依据它决定了预测结果在不同方向、不同尺度上的平滑程度。实际计算时不会真的去求数学期望而是把点对按距离分组比如0到0.05一个箱子、0.05到0.10一个箱子这样分下来每个箱子里的点对半方差取平均得到的散点就是“实验变异函数”。后续要做的是用一个理论模型去拟合这些散点。1.2 普通克里金、简单克里金、泛克里金选哪个Kriging不是一个单一算法而是一个家族。简单克里金假设研究区域的均值已知且恒定比如区域平均气温已知那插值就围绕这个已知均值展开普通克里金假设均值未知但恒定这是绝大多数实际场景下的合理假定泛克里金则处理均值本身存在趋势的情况比如污染物浓度从污染源向外明显衰减时需要先把趋势项剥离出来再插值。三者里普通克里金应用最广原因是它不需要额外输入太多先验知识只需要采样点位置和观测值就能跑。实际项目里我默认先上普通克里金效果不够再考虑泛克里金。如果你直接把简单克里金的均值填错了结果往往还不如普通克里金稳。所以这篇博文的所有代码我全部围绕普通克里金展开。1.3 为什么用MATLAB实现现在实现Kriging的语言选择很多Python有pykrige、GSTools这些现成库R有gstat但MATLAB仍然有一批忠实用户。原因很直接克里金方程组本质上是求解一个线性系统Aλ b这是MATLAB最擅长的运算几条代码就能搞定。再加上MATLAB的surf、contourf、scatter3一套可视化组合拳非常适合快速验证插值效果。不过有一点要提醒网上很多MATLAB版Kriging代码用了pdist2、knnsearch这类函数它们是Statistics and Machine Learning Toolbox里的不是基础MATLAB自带。如果只有基础环境跑起来会直接报错。我这篇文章里的代码尽量避开工具箱函数距离矩阵用手写矩阵运算完成保证绝大多数读者能直接跑。2. 动手编码前先吃透三个数学概念2.1 变异函数的三个关键参数块金值、基台值、变程拟合变异函数时你一定会碰到三个参数块金值、基台值、变程。把它们理解清楚代码才能调得动。块金值nugget描述的是h趋近于0时变异函数仍然存在的跃变值。它来源于测量误差、采样定位误差或者微小尺度上的随机变异。用生活场景类比你问邻居房价两个房子几乎挨着但报价就是会有差异这个差异就是“块金”。基台值sill是变异函数达到的平稳上限可以理解为空间方差的总水平。变程range是从原点开始到变异函数不再明显增加的距离超过变程之后采样点之间已经不存在空间相关性Kriging权重会趋于均匀。拟合球状模型时代码里要估计三个参数块金值(c0)、偏基台值(c也就是基台值减去块金值)、变程(a)。初值的选取很关键我一般先画出实验变异函数散点图目测一下块金值、基台值和变程的大致范围再丢给优化器去精调直接给随机初值很容易陷入局部最优。2.2 球状、指数、高斯变异函数模型怎么选理论变异函数模型有很多种最常用的三个是球状模型、指数模型和高斯模型。它们的表达式和适用场景差别很大我整理了一个对照表模型数学形式h为距离特点适用场景球状γ(h) c0 c(1.5h/a - 0.5(h/a)³)h≤aγ(h)c0cha在变程处严格到达基台曲线有明确拐点地质储量、土壤属性最常用最稳妥指数γ(h) c0 c(1 - exp(-h/a))渐近趋于基台有效变程约为3a污染物扩散、气象要素过渡平滑高斯γ(h) c0 c(1 - exp(-(h/a)²))原点附近非常平滑连续性强高程、水位等高连续性场景但易数值病态实际编程时球状模型最容易处理因为它在h超过a之后直接等于基台值行为很清晰。我大多数项目里第一选择就是球状模型如果交叉验证误差不满意再换指数模型对比。高斯模型虽然听起来高级但它的变异函数在原点附近太平滑导致克里金矩阵严重病态新手容易在这里翻车后面我会专门讲这个问题。2.3 克里金方程组无偏与最优怎么变成线性代数普通克里金的目标是找一组权重λ_i使得预测值等于Σλ_i * Z(x_i)同时满足无偏和预测方差最小两个约束。无偏条件很简单权重之和必须等于1这样才能保证预测值和真实值在均值水平上一致。最小方差条件则需要用拉格朗日乘子法引入一个中间变量μ。最终得到的线性方程组是对每个已知点i都有Σλ_j * γ(x_i - x_j) μ γ(x_i - x_0)再加上Σλ_i 1这行约束。写成矩阵形式就是Aλ bA矩阵是n1阶方阵左上角n×n块是已知点两两之间的变异函数值第n1列和n1行是1右下角是0。b向量前n项是预测点到每个已知点的变异函数值最后一项是1。λ向量前n项就是权重最后一项是拉格朗日乘子。在MATLAB里这个方程就一句话lambda A \ b。不是lambda inv(A)*b直接左除就可以了数值稳定性更好这是MATLAB最常用也是我强烈推荐的做法。3. Kriging的MATLAB完整实现3.1 构造模拟数据用随机采样点验证整个流程要验证Kriging代码最好先有一份已知答案的数据。我在这里构造一个二维正弦场作为“真实分布”然后随机取30个采样点并叠加噪声模拟现实中的观测误差。% 设定随机种子保证结果可复现 rng(42); % 生成30个随机采样点坐标范围 [0,1] x [0,1] n 30; dataXY rand(n, 2); x dataXY(:, 1); y dataXY(:, 2); % 用二维正弦场叠加噪声构造观测值 z sin(2 * pi * x) .* cos(2 * pi * y) 0.3 * randn(n, 1);这里有一个MATLAB新手常踩的坑sin(2*pi*x)里面用的.*是点乘因为x是列向量如果写成sin(2*pi*x) * cos(2*pi*y)MATLAB会试图做矩阵乘法直接报维度不匹配。这类逐元素运算一定要用点乘后面计算距离矩阵时也一样。3.2 求实验变异函数并拟合球状模型有了采样点之后第一步是计算点对之间的距离矩阵和半方差矩阵。% 距离矩阵第(i,j)个元素是第i个点和第j个点的欧氏距离 distMat sqrt((x - x).^2 (y - y).^2); % 半方差矩阵0.5 * (z_i - z_j)^2 gammaMat 0.5 * (z - z).^2; % 按距离分组计算实验变异函数 maxLag 0.6; % 最大滞后距离约为坐标范围的一半 nBins 12; % 分组数量 binEdges linspace(0, maxLag, nBins 1); hExp zeros(nBins, 1); gammaExp zeros(nBins, 1); for k 1:nBins % 去掉自己也去掉两端只取(h,k]区间的点对 idx find(distMat binEdges(k) distMat binEdges(k 1)); if ~isempty(idx) hExp(k) mean(distMat(idx)); gammaExp(k) mean(gammaMat(idx)); else hExp(k) NaN; gammaExp(k) NaN; end end % 删除没有点对的距离区间 valid ~isnan(hExp); hExp hExp(valid); gammaExp gammaExp(valid);% 球状模型函数h0时赋00ha用公式ha取基台值 sphVariogram (h, c0, c, a) ... c0 c * (1.5 * h / a - 0.5 * (h / a).^3) .* (h 0 h a) ... (c0 c) * (h a); % 用fminsearch做最小二乘拟合初值通过目测散点估计 objFun (p) sum((sphVariogram(hExp, p(1), p(2), p(3)) - gammaExp).^2); p0 [0.05, 1, 0.4]; pOpt fminsearch(objFun, p0); c0 pOpt(1); c pOpt(2); a pOpt(3); fprintf(拟合结果块金%.4f, 偏基台值%.4f, 变程%.4f\n, c0, c, a);这里最需要留意的是球状模型在h0处的定义。理论球状模型在原点处取0但公式本身会算出c0所以我在匿名函数里强制做了h 0的判断把原点位置单独处理成0。如果这里不加判断之后构建克里金矩阵时对角线上会出现一个c0直接污染整个求解结果。这个细节代码里体现得特别清楚。3.3 普通克里金预测核心函数预测函数是整段代码的心脏。它接收已知点坐标、观测值、待预测点坐标和变异函数参数返回预测值。function zPred ordinaryKriging(dataXY, dataZ, p0, c0, c, a) n length(dataZ); % 已知点两两之间的距离矩阵 distMat sqrt((dataXY(:,1) - dataXY(:,1)).^2 ... (dataXY(:,2) - dataXY(:,2)).^2); % 构建A矩阵左上角块 A11 c0 c * (1.5 * distMat / a - 0.5 * (distMat / a).^3) .* (distMat 0 distMat a) ... (c0 c) * (distMat a); A11(distMat 0) 0; % 对角线归零 % 拼装完整克里金矩阵 A [A11, ones(n, 1); ones(1, n), 0]; % 预测点与已知点之间的距离向量 distPred sqrt(sum((dataXY - p0).^2, 2)); % 构建b向量 bVec c0 c * (1.5 * distPred / a - 0.5 * (distPred / a).^3) .* (distPred 0 distPred a) ... (c0 c) * (distPred a); bVec [bVec; 1]; % 求解克里金方程组 lambda A \ bVec; % 预测值前n个权重与观测值加权求和最后一个元素是拉格朗日乘子不参与预测 zPred lambda(1:n) * dataZ; end有几个细节值得讲。第一A11(distMat 0) 0这行很关键因为对角线上是点和自己的距离半方差理应为0而我们的匿名函数对h0返回的是c0必须人工覆盖。第二b向量的最后一行是1这与A矩阵最后一行的约束对应用于求解过程中强制权重和为1。第三最后预测时只取lambda的前n个元素最后的拉格朗日乘子虽然参与了求解但不参与加权求和当初我第一次写代码时直接拿lambda全部分量去乘dataZ结果预测值偏得离谱。3.4 网格化批量预测与结果可视化现在把预测函数应用到网格上生成连续的插值曲面。% 生成60x60的预测网格 [Xg, Yg] meshgrid(linspace(0, 1, 60), linspace(0, 1, 60)); Zg zeros(size(Xg)); % 逐点预测 for i 1:numel(Xg) Zg(i) ordinaryKriging(dataXY, z, [Xg(i), Yg(i)], c0, c, a); end % 可视化采样点和预测面叠在一起看 figure; surf(Xg, Yg, Zg, EdgeColor, none); hold on; scatter3(x, y, z, 60, z, filled); colormap(turbo); % 比默认jet的颜色过渡更自然 colorbar; xlabel(X); ylabel(Y); zlabel(Z); title(Kriging插值结果);跑完这段代码你会看到预测曲面在采样点附近自然弯曲、穿过观测值在数据稀疏区域逐渐回归到整个区域的均值水平。这是Kriging的正常行为也解释了为什么Kriging在很多场景下比IDW稳健——它不会在远离数据的地方给出一个突兀的极端值。如果你用scatter3把采样点叠上去还能直观检查预测面是否在采样点处“穿过了”真实观测值。3.5 留一交叉验证评估插值质量插值做完只是第一步评估插值好坏的常规手段是留一交叉验证。做法很简单每次拿掉一个采样点用剩下的n-1个点在该位置做预测再和真实观测值对比最后统计RMSE和MAE。n length(z); predCV zeros(n, 1); for i 1:n idx true(n, 1); idx(i) false; xv dataXY(idx, :); zv z(idx); % 交叉验证时为了速度直接用全样本拟合好的变异函数参数 predCV(i) ordinaryKriging(xv, zv, dataXY(i, :), c0, c, a); end RMSE sqrt(mean((z - predCV).^2)); MAE mean(abs(z - predCV)); fprintf(RMSE%.4f, MAE%.4f\n, RMSE, MAE);交叉验证的严格做法是每一折都重新拟合一次变异函数参数因为理论上去掉一个点后空间结构会发生变化。但在点数量不大的情况下全样本拟合的参数已经比较稳定重新拟合性价比不高。如果你想严谨可以在循环体里把fminsearch步骤加进去只是运行时间会长不少。4. 调试经验与性能优化速查4.1 矩阵奇异、NaN结果是最大拦路虎我在文章开头提到的矩阵奇异问题是新手最容易遇到的第一道坎。常见原因有三个采样点之间存在重复坐标或距离极近的点导致A矩阵两行几乎线性相关。变程a估计过小使得大距离范围内的半方差都落到基台值附近矩阵对角线外的元素差异极小接近秩亏。使用了高斯模型它在原点附近太平滑构造出的矩阵条件数可以轻松到达1e17甚至更高。对策也很直接。第一数据预处理阶段做去重对距离小于某个阈值的点合并或剔除。第二在A矩阵对角线上加一个极小正则化量比如1e-10相当于给求解过程加了一个数值阻尼不影响最终权重但能显著改善条件数。第三如果使用高斯模型导致异常换回球状模型或者指数模型多数情况下能直接解决。A A 1e-10 * eye(size(A));这行代码建议加在求解之前是成本最低的保险。我当时第一次遇到NaN结果时排查了很久最后发现就是高斯模型在作祟换成球状模型后问题消失。4.2 权重为负、预测值出现跳变怎么处理Kriging的一个特性是某些权重可能是负的。这不是bug是数学上允许的结果相当于用周围点的振荡去校正局部趋势。但如果负权重过多或者绝对值过大预测值会莫名跳出样本值范围出现不合理的负浓度、负高度等。遇到这种情况第一件事是检查打印出来的lambda向量。如果某个权重绝对值超过0.5甚至到了1以上说明变异函数模型和参数有明显问题。常见修复路径包括适当增大块金值c0让随机噪声承担更多方差权重会变得更均匀。换用指数模型它的曲线更平滑不容易产生大幅度负权重。检查变程是否过小如果a只有坐标范围的十分之一权重很容易振荡。另外一个容易被忽略的点是预测点位置。Kriging在数据范围内部表现正常但一旦跑到采样点凸包之外外推区域容易出现明显的振荡和偏差。如果有外推需求务必缩小网格范围或者明确告知读者结果仅用于可视化参考。4.3 数据量变大后如何优化性能基础版代码在数据点少几十个的情况下完全够用。但如果采样点超过1000预测网格到几百×几百循环逐点预测的耗时就会猛涨。原因在于每个预测点都要构建一个n1阶线性系统并求解复杂度是O(n³)。优化方向有两个。第一改成局部克里金每个预测点只选择最近的m个采样点参与求解m通常取20到50。这样每次求解的系统都很小同时也符合空间相关性的实际规律——变程以外的点本来就没有多少信息贡献。第二预计算距离矩阵和变异函数矩阵避免重复计算。如果采样点布局固定可以在循环开始前把A矩阵中所有点对半方差算好每个预测点只需更新b向量能省掉大量冗余运算。如果需要处理几万点以上的数据或者网格特别密那就别硬扛MATLAB循环了考虑降维思路先做规则网格的粗插值再用细网格的局部克里金只对局部残差做修正这样可以把计算量降几个数量级。对于大多数项目局部克里金这一个优化已经足够。我个人现在写Kriging代码已经形成固定习惯拿到数据先画实验变异函数散点图目测参数范围再拟合跑完预测必做交叉验证求解前无条件加正则化项模型首选球状再看误差指标决定是否换指数。这套流程在几十到几百个点的项目里非常稳很少再被数值问题卡住。如果你也被Kriging的代码调参折磨过希望这份完整的MATLAB实现能帮你省下至少一个下午的排查时间。本文还有配套的精品资源点击获取