简介MATLAB编写的克里金与协同克里金插值实现代码面向地学、环境、气象等领域的科研人员、研究生及空间统计初学者用于解决单变量和多变量空间插值建模问题。压缩包共8个文件包含7个.m脚本和1个.mat测试数据集整体约44KB代码覆盖克里金主程序、协同克里金程序、变差函数计算与拟合、参数寻优工具及测试脚本结构紧凑便于对照理解完整流程。目前已有2677人学习使用。运行Test.m可快速走通从半方差函数分析、协方差矩阵构建到权重计算的链路Co_kriging0.m与Co_variogram_sq.m则演示如何引入辅助变量提升插值精度配套的Test_data.mat可直接用于复现节省数据准备时间。对于希望从理论向代码落地、或在工程项目中尝试协同克里金方法的读者这套代码提供了可直接修改与复用的参考实现也可作为后续扩展克里金系列算法的起点。 做空间插值这几年我踩过不少坑从最初无脑用IDW反距离加权到后来逐步上手克里金再到为了处理“站点少、但辅助变量多”的难题逼着自己啃协同克里金整个过程最大的体会是克里金这套方法理论门槛看着高但真把原理吃透之后用MATLAB落地并没有想象中那么难。这篇文章就拿我自己在风场数据插值和土壤属性制图中实际跑通的代码为例把克里金和协同克里金的MATLAB实现一次讲透希望看完你也能直接上手。这篇内容适合这几类人刚接触地统计插值、被“变异函数、基台值、块金值”这些术语劝退的初学者已经会用griddata但觉得插值结果不够合理、想做更科学空间估计的研究生以及工作中需要把离散站点数据转成格点场比如气象风场、空气质量监测、环境污染物浓度的工程师和科研人员。1. 克里金插值的核心思想与适用场景1.1 从“反距离加权”到“最优无偏估计”很多人刚开始做空间插值时最顺手的就是IDW距离近的点权重大距离远的点权重小。它逻辑简单、计算快但有个致命问题——权重参数是人为定的完全没考虑数据本身的空间相关性而且对异常值极其敏感。克里金不一样。它名字听着唬人核心思想其实是“用已知点的加权平均去估未知点”但它有两个硬指标无偏性和估计方差最小。这意味着克里金的权重不是靠距离公式拍脑袋算的而是通过半变异函数semivariogram刻画空间自相关结构之后解一个线性方程组求出来的。举个实际例子。有次我做某区域气象站风场数据插值60多个站点用IDW插出来的结果在站点稀疏区域出现明显的“牛眼”现象明明地形平缓的区域被插出一个孤立高值点。换成普通克里金之后因为变异函数对空间连续性的刻画更合理结果平滑自然很多和实际地形分布也更符合。1.2 核心概念变异函数与空间自相关变异函数是克里金方法的灵魂。它的定义是[ \gamma(h) \frac{1}{2N(h)} \sum_{i1}^{N(h)} [Z(x_i) - Z(x_ih)]^2 ]通俗解释就是所有距离为 (h) 的点对其属性值差异的方差的一半。(h) 叫步长lag也就是点对之间的距离。实际计算时我们会先算“实验变异函数”experimental variogram然后再用理论模型去拟合。常用的理论模型有三个模型公式块金值 (c_0)、基台值 (c)、变程 (a)特点球状模型(\gamma(h) c_0 c(1.5h/a - 0.5(h/a)^3))(h \le a)最常用空间相关性随距离先增后平稳指数模型(\gamma(h) c_0 c(1 - e^{-h/a}))相关性衰减快渐近逼近基台值高斯模型(\gamma(h) c_0 c(1 - e^{-(h/a)^2}))曲线平滑适合连续性很强的变量其中变程 (a) 是个关键阈值在变程范围内空间点有相关性超过变程点之间就没有相关性了。这个参数直接决定每个待估点周围有多少已知点参与计算影响非常大。2. 普通克里金与协同克里金怎么选、区别在哪2.1 普通克里金的数学原理普通克里金Ordinary Kriging, OK假设变量的均值是未知但恒定的常数。估值的核心就是求解权重 (\lambda_i)使得[ \hat{Z}(x_0) \sum_{i1}^{n} \lambda_i Z(x_i) ]在无偏约束 (\sum \lambda_i 1) 下通过拉格朗日乘数法求估计方差的最小值最终得到克里金方程组[ \begin{cases} \sum_{j1}^{n} \lambda_j \gamma(x_i, x_j) \mu \gamma(x_i, x_0), \quad i 1,2,...,n \ \sum_{j1}^{n} \lambda_j 1 \end{cases} ]这个方程组里 (\gamma(x_i, x_j)) 就是点与点之间的变异函数值(\mu) 是拉格朗日乘子。MATLAB里解这个方程组非常简单本质上就是矩阵求解A \ b。2.2 协同克里金是如何引入辅助变量的协同克里金Cokriging, CK解决的是“主变量采样点少、辅助变量采样点多”的场景。比如研究土壤重金属含量主变量是重金属浓度可能只有20个采样点辅助变量是土壤有机质含量或地形因子可能有200个采样点。辅助变量和主变量之间有相关性就能把辅助变量的信息“借用”过来提升主变量的插值精度。协同克里金和普通克里金最大的区别是它需要同时考虑主变量自身的变异函数 ( \gamma_{ZZ}(h) )、辅助变量自身的变异函数 ( \gamma_{YY}(h) )以及两者的交叉变异函数 ( \gamma_{ZY}(h) )。求解的方程组从单变量变成了分块矩阵形式复杂度和计算量都明显上升。2.3 适用范围与选择建议我自己使用的经验是这样的主变量站点充足比如100个以上、分布均匀——直接上普通克里金不需要协同克里金。主变量站点少20~50个、但有高密度辅助数据且相关系数明显比如和地形、遥感反演变量相关性超过0.5——果断用协同克里金。主变量和辅助变量相关性太弱低于0.3——别用协同克里金交叉变异函数拟合不出来结果反而更差。还有个实操经验协同克里金提升精度的前提是辅助变量能覆盖整个研究区域如果辅助数据只覆盖一部分区域插值结果会在覆盖边缘出现莫名其妙的“断层”。3. MATLAB代码实现3.1 代码整体结构我把代码拆成几个模块计算实验变异函数拟合理论变异函数模型普通克里金预测协同克里金预测可视化对比这里直接分享我实际在用的核心代码已经做了简化处理去掉了一些无关的业务逻辑方便大家直接看懂思路。3.2 计算实验变异函数以下是计算实验变异函数的MATLAB代码function [h_exp, gamma_exp] compute_variogram(x, y, z, max_dist, n_lags) % 计算实验变异函数 % x, y, z: 坐标和变量值 % max_dist: 最大距离一般取研究区最大距离的一半 % n_lags: 距离分组数 dists pdist2([x y], [x y]); % 所有点对距离矩阵 n length(z); gamma_exp zeros(n_lags, 1); h_exp zeros(n_lags, 1); lag_width max_dist / n_lags; for k 1:n_lags h_min (k-1) * lag_width; h_max k * lag_width; % 筛选距离在区间内的点对 [I, J] find(dists h_min dists h_max); if isempty(I) h_exp(k) NaN; gamma_exp(k) NaN; continue; end % 半方差: 值差异平方的一半再平均 gamma_exp(k) 0.5 * mean((z(I) - z(J)).^2); h_exp(k) (h_min h_max) / 2; end % 去掉空的区间 valid ~isnan(h_exp); h_exp h_exp(valid); gamma_exp gamma_exp(valid); end这块有个细节容易被忽略计算距离矩阵时数据量一大pdist2生成的是 (n \times n) 矩阵(n1000) 就是100万个数占内存8MB还好但 (n10000) 就是8亿个数6.4GB电脑直接卡死。所以站点超过3000个的时候建议用循环分段计算距离而不是一次性生成全距离矩阵。3.3 拟合理论变异函数模型计算完实验变异函数后需要拟合理论模型。这里我用lsqcurvefit做拟合MATLAB自带的优化工具箱就能实现function [params, model] fit_variogram(h_exp, gamma_exp, model_type) % 拟合理论变异函数模型 % model_type: spherical, exponential, gaussian % 初始值: 块金值取最小半方差的50%基台值取最大半方差的90%变程取最大距离的1/3 c0_init 0.5 * min(gamma_exp); c_init 0.9 * max(gamma_exp) - c0_init; a_init max(h_exp) / 3; params0 [c0_init, c_init, a_init]; % 模型函数句柄 switch model_type case spherical model (p, h) p(1) p(2) * (1.5*h./max(p(3),eps) - 0.5*(h./max(p(3),eps)).^3) .* (h p(3)) p(2) .* (h p(3)); case exponential model (p, h) p(1) p(2) * (1 - exp(-h ./ max(p(3),eps))); case gaussian model (p, h) p(1) p(2) * (1 - exp(-(h.^2) ./ max(p(3)^2, eps))); end % 非线性最小二乘拟合 lb [0, 0, 0]; % 参数非负 ub [inf, inf, max(h_exp)*2]; options optimoptions(lsqcurvefit, Display, off, MaxIterations, 500); params lsqcurvefit(model, params0, h_exp, gamma_exp, lb, ub, options); end拟合并不能保证总是一次成功。我用的经验是如果拟合出来的变程非常小比如只有最大距离的1/10说明空间相关性很弱这时候即便插值出来了可信度也有限如果块金值占比太高超过基台值的50%说明数据噪音大或采样尺度不合适最好回头检查数据质量。3.4 普通克里金预测函数这是所有代码里最核心的一段其实就是解克里金方程组function [Zhat, Var_est] kriging_ok(x_obs, y_obs, z_obs, x_pred, y_pred, params, model) % 普通克里金插值预测 % params [c0, c, a] n length(z_obs); m length(x_pred); Zhat zeros(m, 1); Var_est zeros(m, 1); % 构建左侧矩阵 A点与点之间的变异函数) dist_obs pdist2([x_obs y_obs], [x_obs y_obs]); A model(params, dist_obs); A A eye(n) * 1e-10; % 增加微小扰动防止矩阵奇异 A [A, ones(n,1); ones(1,n), 0]; % 加入拉格朗日约束行/列 b_vec zeros(n1, 1); b_vec(n1) 1; % 对待插值点循环求解 for i 1:m dist_pred sqrt((x_obs - x_pred(i)).^2 (y_obs - y_pred(i)).^2); b_vec(1:n) model(params, dist_pred); weights A \ b_vec; % 求解克里金权重 Zhat(i) sum(weights(1:n) .* z_obs); Var_est(i) sum(weights(1:n) .* b_vec(1:n)) weights(n1); end end这个函数里有一个每次都要提的关键点矩阵 (A) 只和观测点的相对位置有关对于不同的待插值点只需要更新方程右端项 (b) 就行。如果有10000个格点要预测把 (A) 的逆一次算好然后反复调用速度会快很多而且能避免反复构造矩阵带来的NaN风险。3.5 协同克里金的实现协同克里金的核心是构建分块方程组。我们用 (Z) 表示主变量(Y) 表示辅助变量方程组的左侧变成[ A \begin{bmatrix} \gamma_{ZZ}(h_{ij}) \gamma_{ZY}(h_{ij}) 1 0 \ \gamma_{YZ}(h_{ij}) \gamma_{YY}(h_{ij}) 0 1 \ 1 0 0 0 \ 0 1 0 0 \end{bmatrix} ]其中交叉变异函数 (\gamma_{ZY}(h)) 的计算公式是[ \gamma_{ZY}(h) \frac{1}{2N(h)} \sum_{i1}^{N(h)} [Z(x_i) - Z(x_ih)][Y(x_i) - Y(x_ih)] ]核心代码function [Zhat] cokriging_xy(x_obs, y_obs, z_obs, x_aux, y_aux, z_aux, ... x_pred, y_pred, model_Z, model_Y, model_ZY, ... params_Z, params_Y, params_ZY) % 协同克里金插值 % 注意: 需要保证辅助变量和主变量在同一位置都有观测值或者至少距离足够近 % 实际项目中建议先对辅助变量做普通克里金统一到主变量站点坐标上 n length(z_obs); % 主变量站点数 % 这里简化处理: 假设在主变量站点位置上辅助变量值已知 % 如果辅助变量不在主变量站点上需要先插值获取 z_aux_at_obs z_aux(1:n); % 实际使用中通过插值得到 m length(x_pred); Zhat zeros(m, 1); % 构建分块矩阵 dist_obs pdist2([x_obs y_obs], [x_obs y_obs]); A11 model_Z(params_Z, dist_obs); % 主变量变异函数矩阵 A12 model_ZY(params_ZY, dist_obs); % 交叉变异函数矩阵 A21 A12; A22 model_Y(params_Y, dist_obs); % 辅助变量变异函数矩阵 % 组装大矩阵 A [A11, A12, ones(n,1), zeros(n,1); A21, A22, zeros(n,1), ones(n,1); ones(1,n), zeros(1,n), 0, 0; zeros(1,n), ones(1,n), 0, 0]; A A eye(size(A)) * 1e-10; for i 1:m dist_pred sqrt((x_obs - x_pred(i)).^2 (y_obs - y_pred(i)).^2); r1 model_Z(params_Z, dist_pred); % 主变量变异函数 r2 model_ZY(params_ZY, dist_pred); % 交叉变异函数 b [r1; r2; 1; 0]; weights A \ b; % 主变量权重和辅助变量权重分别取前n和后n Zhat(i) sum(weights(1:n) .* z_obs) sum(weights(n1:2*n) .* z_aux_at_obs); end end协同克里金调用前有个前置工作特别重要必须保证在主变量站点位置上辅助变量的值是已知的。如果辅助变量站点和主变量站点不重合大部分时候都不重合需要先把辅助变量插值到主变量站点位置否则交叉变异函数算不出来。3.6 完整调用流程示例% 模拟数据生成 load(wind_station_data.mat); % 假设有 x_sta, y_sta, wind_speed % 生成目标格点 [xg, yg] meshgrid(linspace(min(x_sta), max(x_sta), 100), ... linspace(min(y_sta), max(y_sta), 100)); % 1. 计算实验变异函数 max_dist max(pdist2([x_sta y_sta], [x_sta y_sta]), [], all) * 0.5; [h_exp, g_exp] compute_variogram(x_sta, y_sta, wind_speed, max_dist, 20); % 2. 拟合模型 [params, model] fit_variogram(h_exp, g_exp, spherical); fprintf(拟合结果: 块金值%.2f, 基台值%.2f, 变程%.2f\n, params(1), params(2), params(3)); % 3. 克里金预测 [Z_grid, Var_grid] kriging_ok(x_sta, y_sta, wind_speed, xg(:), yg(:), params, model); Z_grid reshape(Z_grid, size(xg)); Var_grid reshape(Var_grid, size(xg)); % 4. 可视化 figure; subplot(1,2,1); contourf(xg, yg, Z_grid, 20, LineColor, none); colorbar; title(克里金插值结果); hold on; scatter(x_sta, y_sta, 20, wind_speed, filled, k); subplot(1,2,2); contourf(xg, yg, sqrt(Var_grid), 20, LineColor, none); colorbar; title(克里金方差(标准差));4. 实操经验与参数调优4.1 数据预处理的三个关键步骤克里金对数据质量极其敏感。以下三个预处理步骤每次都要做第一检查数据分布。克里金在理论上并不要求数据严格正态但强烈建议偏态严重的数据做对数变换或Box-Cox变换。我遇到过土壤重金属数据偏度系数达到3.5的情况直接插值结果全是“离群值拉偏”对数变换后插值再反变换回原尺度效果好了很多。第二剔除明显异常值。用3倍标准差或局部空间异常检测比如和周围8个点的均值差超过3倍标准差把异常点筛出来。这些异常点会让变异函数在短距离上出现很大的半方差直接拉高块金值导致插值结果平滑过头。第三数据去趋势。如果变量在空间上有明显的大尺度趋势比如温度随纬度线性降低先用多项式拟合去掉趋势项对残差做克里金最后再把趋势加回去。这个操作在气象领域叫“回归克里金”的简化版能大幅提升插值精度。4.2 变异函数模型的参数初始值选择用lsqcurvefit拟合时初始值选的好坏影响很大。我总结了一个经验法则块金值初始值取最小半方差的一半左右基台值初始值取最大半方差的90%减去块金值变程初始值取最大距离的1/3如果拟合结果出现块金值为0的情况别高兴太早这通常意味着数据在极小距离上仍然高度相关是好事但也可能说明采样尺度太密、存在重复采样。如果变程拟合出非常离谱的值比如大于最大距离的两倍被上限约束到还是硬顶着这时候要考虑换模型或者重新审视数据的空间平稳性。4.3 邻域搜索别把所有点都塞进方程组初学者最容易犯的错误就是把所有观测点都放进克里金方程组。这样做有三个问题矩阵阶数太大导致计算慢离待估点很远的点权重接近0纯属浪费更严重的是远距离点上如果有异常值会把整个估计带偏。工程上建议用“邻域搜索”对待插值点只取半径 (R) 内最近的 (k) 个点参与计算。(R) 一般取变程的1~1.5倍(k) 建议16~32。我实测过对于300个站点、10000个格点的场景全局求解需要几秒钟邻域搜索后只需要不到0.5秒精度几乎不变。在MATLAB里搜索邻域点用knnsearch或rangesearch都行效率远高于自己写循环遍历。5. 常见问题与排查技巧5.1 结果出现大面积NaN或异常大值这种情况绝大多数是矩阵奇异导致的。克里金方程组矩阵在两种情况下最容易奇异一是多个观测点距离太近重合或几乎重合导致矩阵行近似线性相关二是数据点数量太少少于5个矩阵本身就不可靠。排查方法很简单先检查是否有重复坐标点合并掉再在矩阵对角线上加一个微小扰动比如1e-10的量级我上面的代码里已经加了。如果还有问题把扰动给大一些到1e-6但别太大否则插值结果会失真。5.2 批量克里金插值速度太慢做批量克里金插值比如多个时间步长的风场数据时一个常见问题是每个时间步都重新拟合一次变异函数并重建矩阵。实际上如果站点坐标不变变异函数应该相差不大。我的做法是用第一个时间步的数据拟合变异函数后面的时间步直接用第一组参数只更新右端项求解。这样速度能提升5到10倍。如果数据差异较大再考虑每隔几个时间步重新拟合一次。5.3 协同克里金结果比普通克里金还差这是很多人用协同克里金后最挫败的时刻。我排查过多次主要原因几乎都是辅助变量和主变量的相关性不够强或者交叉变异函数拟合得太差。交叉变异函数的计算本身就容易受到两变量噪声影响如果相关系数低于0.5交叉变异函数经常乱成一团。另一个容易被忽略的点是量纲问题。主变量是风速单位m/s辅助变量是气压单位hPa数值范围差异巨大直接算交叉变异函数量级小的变量会被“淹没”。解决办法是对两个变量做标准化处理减均值、除以标准差后再进行协同克里金。5.4 克里金方差和预想不一致克里金方差反映的是配置优劣和信息量不是变量本身的真实方差。如果某片区域站点密集克里金方差会很小如果站点稀疏方差就会变大。这是正常的。但如果你发现克里金方差在已知站点位置上不等于0理论上在已知点上方差应该为0那基本上是因为加了正则化扰动项或者坐标精度问题导致程序没有识别出已知点。6. 写在最后的个人体会我在实际项目里用这套代码处理过气象风场、空气质量监测和土壤采样数据最大的体会是克里金的精度提升并不在于模型选得多花哨而在于变异函数拟合得是否符合物理实际。每次插值前我都会把实验变异函数和拟合模型画出来看一眼如果拟合曲线和散点明显偏离再好的插值结果也是“垃圾进、垃圾出”。如果后续你想继续深挖建议在现有代码基础上增加三个功能一是支持带趋势项的泛克里金Universal Kriging二是把邻域搜索加进去支持更大规模的格点预测三是把代码封装成函数库方便批量处理不同时次的数据。沿着这个方向扩展你就能拥有一套属于自己的、趁手的空间插值工具包。本文还有配套的精品资源点击获取 
长河君 · 2026/9/9 21:56:52 · 阅读 1.2k 克里金与协同克里金MATLAB实现：从原理到代码实战
简介MATLAB编写的克里金与协同克里金插值实现代码面向地学、环境、气象等领域的科研人员、研究生及空间统计初学者用于解决单变量和多变量空间插值建模问题。压缩包共8个文件包含7个.m脚本和1个.mat测试数据集整体约44KB代码覆盖克里金主程序、协同克里金程序、变差函数计算与拟合、参数寻优工具及测试脚本结构紧凑便于对照理解完整流程。目前已有2677人学习使用。运行Test.m可快速走通从半方差函数分析、协方差矩阵构建到权重计算的链路Co_kriging0.m与Co_variogram_sq.m则演示如何引入辅助变量提升插值精度配套的Test_data.mat可直接用于复现节省数据准备时间。对于希望从理论向代码落地、或在工程项目中尝试协同克里金方法的读者这套代码提供了可直接修改与复用的参考实现也可作为后续扩展克里金系列算法的起点。 做空间插值这几年我踩过不少坑从最初无脑用IDW反距离加权到后来逐步上手克里金再到为了处理“站点少、但辅助变量多”的难题逼着自己啃协同克里金整个过程最大的体会是克里金这套方法理论门槛看着高但真把原理吃透之后用MATLAB落地并没有想象中那么难。这篇文章就拿我自己在风场数据插值和土壤属性制图中实际跑通的代码为例把克里金和协同克里金的MATLAB实现一次讲透希望看完你也能直接上手。这篇内容适合这几类人刚接触地统计插值、被“变异函数、基台值、块金值”这些术语劝退的初学者已经会用griddata但觉得插值结果不够合理、想做更科学空间估计的研究生以及工作中需要把离散站点数据转成格点场比如气象风场、空气质量监测、环境污染物浓度的工程师和科研人员。1. 克里金插值的核心思想与适用场景1.1 从“反距离加权”到“最优无偏估计”很多人刚开始做空间插值时最顺手的就是IDW距离近的点权重大距离远的点权重小。它逻辑简单、计算快但有个致命问题——权重参数是人为定的完全没考虑数据本身的空间相关性而且对异常值极其敏感。克里金不一样。它名字听着唬人核心思想其实是“用已知点的加权平均去估未知点”但它有两个硬指标无偏性和估计方差最小。这意味着克里金的权重不是靠距离公式拍脑袋算的而是通过半变异函数semivariogram刻画空间自相关结构之后解一个线性方程组求出来的。举个实际例子。有次我做某区域气象站风场数据插值60多个站点用IDW插出来的结果在站点稀疏区域出现明显的“牛眼”现象明明地形平缓的区域被插出一个孤立高值点。换成普通克里金之后因为变异函数对空间连续性的刻画更合理结果平滑自然很多和实际地形分布也更符合。1.2 核心概念变异函数与空间自相关变异函数是克里金方法的灵魂。它的定义是[ \gamma(h) \frac{1}{2N(h)} \sum_{i1}^{N(h)} [Z(x_i) - Z(x_ih)]^2 ]通俗解释就是所有距离为 (h) 的点对其属性值差异的方差的一半。(h) 叫步长lag也就是点对之间的距离。实际计算时我们会先算“实验变异函数”experimental variogram然后再用理论模型去拟合。常用的理论模型有三个模型公式块金值 (c_0)、基台值 (c)、变程 (a)特点球状模型(\gamma(h) c_0 c(1.5h/a - 0.5(h/a)^3))(h \le a)最常用空间相关性随距离先增后平稳指数模型(\gamma(h) c_0 c(1 - e^{-h/a}))相关性衰减快渐近逼近基台值高斯模型(\gamma(h) c_0 c(1 - e^{-(h/a)^2}))曲线平滑适合连续性很强的变量其中变程 (a) 是个关键阈值在变程范围内空间点有相关性超过变程点之间就没有相关性了。这个参数直接决定每个待估点周围有多少已知点参与计算影响非常大。2. 普通克里金与协同克里金怎么选、区别在哪2.1 普通克里金的数学原理普通克里金Ordinary Kriging, OK假设变量的均值是未知但恒定的常数。估值的核心就是求解权重 (\lambda_i)使得[ \hat{Z}(x_0) \sum_{i1}^{n} \lambda_i Z(x_i) ]在无偏约束 (\sum \lambda_i 1) 下通过拉格朗日乘数法求估计方差的最小值最终得到克里金方程组[ \begin{cases} \sum_{j1}^{n} \lambda_j \gamma(x_i, x_j) \mu \gamma(x_i, x_0), \quad i 1,2,...,n \ \sum_{j1}^{n} \lambda_j 1 \end{cases} ]这个方程组里 (\gamma(x_i, x_j)) 就是点与点之间的变异函数值(\mu) 是拉格朗日乘子。MATLAB里解这个方程组非常简单本质上就是矩阵求解A \ b。2.2 协同克里金是如何引入辅助变量的协同克里金Cokriging, CK解决的是“主变量采样点少、辅助变量采样点多”的场景。比如研究土壤重金属含量主变量是重金属浓度可能只有20个采样点辅助变量是土壤有机质含量或地形因子可能有200个采样点。辅助变量和主变量之间有相关性就能把辅助变量的信息“借用”过来提升主变量的插值精度。协同克里金和普通克里金最大的区别是它需要同时考虑主变量自身的变异函数 ( \gamma_{ZZ}(h) )、辅助变量自身的变异函数 ( \gamma_{YY}(h) )以及两者的交叉变异函数 ( \gamma_{ZY}(h) )。求解的方程组从单变量变成了分块矩阵形式复杂度和计算量都明显上升。2.3 适用范围与选择建议我自己使用的经验是这样的主变量站点充足比如100个以上、分布均匀——直接上普通克里金不需要协同克里金。主变量站点少20~50个、但有高密度辅助数据且相关系数明显比如和地形、遥感反演变量相关性超过0.5——果断用协同克里金。主变量和辅助变量相关性太弱低于0.3——别用协同克里金交叉变异函数拟合不出来结果反而更差。还有个实操经验协同克里金提升精度的前提是辅助变量能覆盖整个研究区域如果辅助数据只覆盖一部分区域插值结果会在覆盖边缘出现莫名其妙的“断层”。3. MATLAB代码实现3.1 代码整体结构我把代码拆成几个模块计算实验变异函数拟合理论变异函数模型普通克里金预测协同克里金预测可视化对比这里直接分享我实际在用的核心代码已经做了简化处理去掉了一些无关的业务逻辑方便大家直接看懂思路。3.2 计算实验变异函数以下是计算实验变异函数的MATLAB代码function [h_exp, gamma_exp] compute_variogram(x, y, z, max_dist, n_lags) % 计算实验变异函数 % x, y, z: 坐标和变量值 % max_dist: 最大距离一般取研究区最大距离的一半 % n_lags: 距离分组数 dists pdist2([x y], [x y]); % 所有点对距离矩阵 n length(z); gamma_exp zeros(n_lags, 1); h_exp zeros(n_lags, 1); lag_width max_dist / n_lags; for k 1:n_lags h_min (k-1) * lag_width; h_max k * lag_width; % 筛选距离在区间内的点对 [I, J] find(dists h_min dists h_max); if isempty(I) h_exp(k) NaN; gamma_exp(k) NaN; continue; end % 半方差: 值差异平方的一半再平均 gamma_exp(k) 0.5 * mean((z(I) - z(J)).^2); h_exp(k) (h_min h_max) / 2; end % 去掉空的区间 valid ~isnan(h_exp); h_exp h_exp(valid); gamma_exp gamma_exp(valid); end这块有个细节容易被忽略计算距离矩阵时数据量一大pdist2生成的是 (n \times n) 矩阵(n1000) 就是100万个数占内存8MB还好但 (n10000) 就是8亿个数6.4GB电脑直接卡死。所以站点超过3000个的时候建议用循环分段计算距离而不是一次性生成全距离矩阵。3.3 拟合理论变异函数模型计算完实验变异函数后需要拟合理论模型。这里我用lsqcurvefit做拟合MATLAB自带的优化工具箱就能实现function [params, model] fit_variogram(h_exp, gamma_exp, model_type) % 拟合理论变异函数模型 % model_type: spherical, exponential, gaussian % 初始值: 块金值取最小半方差的50%基台值取最大半方差的90%变程取最大距离的1/3 c0_init 0.5 * min(gamma_exp); c_init 0.9 * max(gamma_exp) - c0_init; a_init max(h_exp) / 3; params0 [c0_init, c_init, a_init]; % 模型函数句柄 switch model_type case spherical model (p, h) p(1) p(2) * (1.5*h./max(p(3),eps) - 0.5*(h./max(p(3),eps)).^3) .* (h p(3)) p(2) .* (h p(3)); case exponential model (p, h) p(1) p(2) * (1 - exp(-h ./ max(p(3),eps))); case gaussian model (p, h) p(1) p(2) * (1 - exp(-(h.^2) ./ max(p(3)^2, eps))); end % 非线性最小二乘拟合 lb [0, 0, 0]; % 参数非负 ub [inf, inf, max(h_exp)*2]; options optimoptions(lsqcurvefit, Display, off, MaxIterations, 500); params lsqcurvefit(model, params0, h_exp, gamma_exp, lb, ub, options); end拟合并不能保证总是一次成功。我用的经验是如果拟合出来的变程非常小比如只有最大距离的1/10说明空间相关性很弱这时候即便插值出来了可信度也有限如果块金值占比太高超过基台值的50%说明数据噪音大或采样尺度不合适最好回头检查数据质量。3.4 普通克里金预测函数这是所有代码里最核心的一段其实就是解克里金方程组function [Zhat, Var_est] kriging_ok(x_obs, y_obs, z_obs, x_pred, y_pred, params, model) % 普通克里金插值预测 % params [c0, c, a] n length(z_obs); m length(x_pred); Zhat zeros(m, 1); Var_est zeros(m, 1); % 构建左侧矩阵 A点与点之间的变异函数) dist_obs pdist2([x_obs y_obs], [x_obs y_obs]); A model(params, dist_obs); A A eye(n) * 1e-10; % 增加微小扰动防止矩阵奇异 A [A, ones(n,1); ones(1,n), 0]; % 加入拉格朗日约束行/列 b_vec zeros(n1, 1); b_vec(n1) 1; % 对待插值点循环求解 for i 1:m dist_pred sqrt((x_obs - x_pred(i)).^2 (y_obs - y_pred(i)).^2); b_vec(1:n) model(params, dist_pred); weights A \ b_vec; % 求解克里金权重 Zhat(i) sum(weights(1:n) .* z_obs); Var_est(i) sum(weights(1:n) .* b_vec(1:n)) weights(n1); end end这个函数里有一个每次都要提的关键点矩阵 (A) 只和观测点的相对位置有关对于不同的待插值点只需要更新方程右端项 (b) 就行。如果有10000个格点要预测把 (A) 的逆一次算好然后反复调用速度会快很多而且能避免反复构造矩阵带来的NaN风险。3.5 协同克里金的实现协同克里金的核心是构建分块方程组。我们用 (Z) 表示主变量(Y) 表示辅助变量方程组的左侧变成[ A \begin{bmatrix} \gamma_{ZZ}(h_{ij}) \gamma_{ZY}(h_{ij}) 1 0 \ \gamma_{YZ}(h_{ij}) \gamma_{YY}(h_{ij}) 0 1 \ 1 0 0 0 \ 0 1 0 0 \end{bmatrix} ]其中交叉变异函数 (\gamma_{ZY}(h)) 的计算公式是[ \gamma_{ZY}(h) \frac{1}{2N(h)} \sum_{i1}^{N(h)} [Z(x_i) - Z(x_ih)][Y(x_i) - Y(x_ih)] ]核心代码function [Zhat] cokriging_xy(x_obs, y_obs, z_obs, x_aux, y_aux, z_aux, ... x_pred, y_pred, model_Z, model_Y, model_ZY, ... params_Z, params_Y, params_ZY) % 协同克里金插值 % 注意: 需要保证辅助变量和主变量在同一位置都有观测值或者至少距离足够近 % 实际项目中建议先对辅助变量做普通克里金统一到主变量站点坐标上 n length(z_obs); % 主变量站点数 % 这里简化处理: 假设在主变量站点位置上辅助变量值已知 % 如果辅助变量不在主变量站点上需要先插值获取 z_aux_at_obs z_aux(1:n); % 实际使用中通过插值得到 m length(x_pred); Zhat zeros(m, 1); % 构建分块矩阵 dist_obs pdist2([x_obs y_obs], [x_obs y_obs]); A11 model_Z(params_Z, dist_obs); % 主变量变异函数矩阵 A12 model_ZY(params_ZY, dist_obs); % 交叉变异函数矩阵 A21 A12; A22 model_Y(params_Y, dist_obs); % 辅助变量变异函数矩阵 % 组装大矩阵 A [A11, A12, ones(n,1), zeros(n,1); A21, A22, zeros(n,1), ones(n,1); ones(1,n), zeros(1,n), 0, 0; zeros(1,n), ones(1,n), 0, 0]; A A eye(size(A)) * 1e-10; for i 1:m dist_pred sqrt((x_obs - x_pred(i)).^2 (y_obs - y_pred(i)).^2); r1 model_Z(params_Z, dist_pred); % 主变量变异函数 r2 model_ZY(params_ZY, dist_pred); % 交叉变异函数 b [r1; r2; 1; 0]; weights A \ b; % 主变量权重和辅助变量权重分别取前n和后n Zhat(i) sum(weights(1:n) .* z_obs) sum(weights(n1:2*n) .* z_aux_at_obs); end end协同克里金调用前有个前置工作特别重要必须保证在主变量站点位置上辅助变量的值是已知的。如果辅助变量站点和主变量站点不重合大部分时候都不重合需要先把辅助变量插值到主变量站点位置否则交叉变异函数算不出来。3.6 完整调用流程示例% 模拟数据生成 load(wind_station_data.mat); % 假设有 x_sta, y_sta, wind_speed % 生成目标格点 [xg, yg] meshgrid(linspace(min(x_sta), max(x_sta), 100), ... linspace(min(y_sta), max(y_sta), 100)); % 1. 计算实验变异函数 max_dist max(pdist2([x_sta y_sta], [x_sta y_sta]), [], all) * 0.5; [h_exp, g_exp] compute_variogram(x_sta, y_sta, wind_speed, max_dist, 20); % 2. 拟合模型 [params, model] fit_variogram(h_exp, g_exp, spherical); fprintf(拟合结果: 块金值%.2f, 基台值%.2f, 变程%.2f\n, params(1), params(2), params(3)); % 3. 克里金预测 [Z_grid, Var_grid] kriging_ok(x_sta, y_sta, wind_speed, xg(:), yg(:), params, model); Z_grid reshape(Z_grid, size(xg)); Var_grid reshape(Var_grid, size(xg)); % 4. 可视化 figure; subplot(1,2,1); contourf(xg, yg, Z_grid, 20, LineColor, none); colorbar; title(克里金插值结果); hold on; scatter(x_sta, y_sta, 20, wind_speed, filled, k); subplot(1,2,2); contourf(xg, yg, sqrt(Var_grid), 20, LineColor, none); colorbar; title(克里金方差(标准差));4. 实操经验与参数调优4.1 数据预处理的三个关键步骤克里金对数据质量极其敏感。以下三个预处理步骤每次都要做第一检查数据分布。克里金在理论上并不要求数据严格正态但强烈建议偏态严重的数据做对数变换或Box-Cox变换。我遇到过土壤重金属数据偏度系数达到3.5的情况直接插值结果全是“离群值拉偏”对数变换后插值再反变换回原尺度效果好了很多。第二剔除明显异常值。用3倍标准差或局部空间异常检测比如和周围8个点的均值差超过3倍标准差把异常点筛出来。这些异常点会让变异函数在短距离上出现很大的半方差直接拉高块金值导致插值结果平滑过头。第三数据去趋势。如果变量在空间上有明显的大尺度趋势比如温度随纬度线性降低先用多项式拟合去掉趋势项对残差做克里金最后再把趋势加回去。这个操作在气象领域叫“回归克里金”的简化版能大幅提升插值精度。4.2 变异函数模型的参数初始值选择用lsqcurvefit拟合时初始值选的好坏影响很大。我总结了一个经验法则块金值初始值取最小半方差的一半左右基台值初始值取最大半方差的90%减去块金值变程初始值取最大距离的1/3如果拟合结果出现块金值为0的情况别高兴太早这通常意味着数据在极小距离上仍然高度相关是好事但也可能说明采样尺度太密、存在重复采样。如果变程拟合出非常离谱的值比如大于最大距离的两倍被上限约束到还是硬顶着这时候要考虑换模型或者重新审视数据的空间平稳性。4.3 邻域搜索别把所有点都塞进方程组初学者最容易犯的错误就是把所有观测点都放进克里金方程组。这样做有三个问题矩阵阶数太大导致计算慢离待估点很远的点权重接近0纯属浪费更严重的是远距离点上如果有异常值会把整个估计带偏。工程上建议用“邻域搜索”对待插值点只取半径 (R) 内最近的 (k) 个点参与计算。(R) 一般取变程的1~1.5倍(k) 建议16~32。我实测过对于300个站点、10000个格点的场景全局求解需要几秒钟邻域搜索后只需要不到0.5秒精度几乎不变。在MATLAB里搜索邻域点用knnsearch或rangesearch都行效率远高于自己写循环遍历。5. 常见问题与排查技巧5.1 结果出现大面积NaN或异常大值这种情况绝大多数是矩阵奇异导致的。克里金方程组矩阵在两种情况下最容易奇异一是多个观测点距离太近重合或几乎重合导致矩阵行近似线性相关二是数据点数量太少少于5个矩阵本身就不可靠。排查方法很简单先检查是否有重复坐标点合并掉再在矩阵对角线上加一个微小扰动比如1e-10的量级我上面的代码里已经加了。如果还有问题把扰动给大一些到1e-6但别太大否则插值结果会失真。5.2 批量克里金插值速度太慢做批量克里金插值比如多个时间步长的风场数据时一个常见问题是每个时间步都重新拟合一次变异函数并重建矩阵。实际上如果站点坐标不变变异函数应该相差不大。我的做法是用第一个时间步的数据拟合变异函数后面的时间步直接用第一组参数只更新右端项求解。这样速度能提升5到10倍。如果数据差异较大再考虑每隔几个时间步重新拟合一次。5.3 协同克里金结果比普通克里金还差这是很多人用协同克里金后最挫败的时刻。我排查过多次主要原因几乎都是辅助变量和主变量的相关性不够强或者交叉变异函数拟合得太差。交叉变异函数的计算本身就容易受到两变量噪声影响如果相关系数低于0.5交叉变异函数经常乱成一团。另一个容易被忽略的点是量纲问题。主变量是风速单位m/s辅助变量是气压单位hPa数值范围差异巨大直接算交叉变异函数量级小的变量会被“淹没”。解决办法是对两个变量做标准化处理减均值、除以标准差后再进行协同克里金。5.4 克里金方差和预想不一致克里金方差反映的是配置优劣和信息量不是变量本身的真实方差。如果某片区域站点密集克里金方差会很小如果站点稀疏方差就会变大。这是正常的。但如果你发现克里金方差在已知站点位置上不等于0理论上在已知点上方差应该为0那基本上是因为加了正则化扰动项或者坐标精度问题导致程序没有识别出已知点。6. 写在最后的个人体会我在实际项目里用这套代码处理过气象风场、空气质量监测和土壤采样数据最大的体会是克里金的精度提升并不在于模型选得多花哨而在于变异函数拟合得是否符合物理实际。每次插值前我都会把实验变异函数和拟合模型画出来看一眼如果拟合曲线和散点明显偏离再好的插值结果也是“垃圾进、垃圾出”。如果后续你想继续深挖建议在现有代码基础上增加三个功能一是支持带趋势项的泛克里金Universal Kriging二是把邻域搜索加进去支持更大规模的格点预测三是把代码封装成函数库方便批量处理不同时次的数据。沿着这个方向扩展你就能拥有一套属于自己的、趁手的空间插值工具包。本文还有配套的精品资源点击获取