企业官网建设流程全解析

简介：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_{i=1}^{N(h)} [Z(x_i) - Z(x_i+h)]^2 ]

通俗解释就是：所有距离为 (h) 的点对，其属性值差异的方差的一半。(h) 叫步长（lag），也就是点对之间的距离。

实际计算时，我们会先算“实验变异函数”（experimental variogram），然后再用理论模型去拟合。常用的理论模型有三个：

其中变程 (a) 是个关键阈值：在变程范围内，空间点有相关性；超过变程，点之间就没有相关性了。这个参数直接决定每个待估点周围有多少已知点参与计算，影响非常大。

2. 普通克里金与协同克里金：怎么选、区别在哪

2.1 普通克里金的数学原理

普通克里金（Ordinary Kriging, OK）假设变量的均值是未知但恒定的常数。估值的核心就是求解权重 (\lambda_i)，使得：

[ \hat{Z}(x_0) = \sum_{i=1}^{n} \lambda_i Z(x_i) ]

在无偏约束 (\sum \lambda_i = 1) 下，通过拉格朗日乘数法求估计方差的最小值，最终得到克里金方程组：

[ \begin{cases} \sum_{j=1}^{n} \lambda_j \gamma(x_i, x_j) + \mu = \gamma(x_i, x_0), \quad i = 1,2,...,n \ \sum_{j=1}^{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 代码整体结构

我把代码拆成几个模块：

1. 计算实验变异函数
2. 拟合理论变异函数模型
3. 普通克里金预测
4. 协同克里金预测
5. 可视化对比

这里直接分享我实际在用的核心代码，已经做了简化处理，去掉了一些无关的业务逻辑，方便大家直接看懂思路。

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) 矩阵，(n=1000) 就是100万个数，占内存8MB，还好；但 (n=10000) 就是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(n+1, 1); b_vec(n+1) = 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(n+1); 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_{i=1}^{N(h)} [Z(x_i) - Z(x_i+h)][Y(x_i) - Y(x_i+h)] ]

核心代码：

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(n+1: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）；二是把邻域搜索加进去，支持更大规模的格点预测；三是把代码封装成函数库，方便批量处理不同时次的数据。沿着这个方向扩展，你就能拥有一套属于自己的、趁手的空间插值工具包。

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