简介这份资源围绕Kriging空间插值方法展开提供一套可直接运行的等值线图绘制源码面向地理信息系统、地质勘探及空间数据分析方向的学习者与开发者帮助解决从离散观测点生成连续空间分布图的实际问题。压缩包共84个文件以h头文件与cpp源文件为主体辅以obj编译中间文件、ico与cur图标资源、txt测试数据及rc资源脚本等整体约1.77MB工程结构完整便于在VC环境下编译调试。目前已有1214人学习下载具备一定参考热度。源码覆盖数据预处理、半方差函数构建、Kriging类型选择、插值模型求解与等值线绘制等关键环节并包含矩阵运算、逆距离插值等辅助模块读者可据此理解空间插值的实现逻辑调整参数观察插值效果并在此基础上扩展自己的等值线可视化应用。1. 从散点到连续曲面Kriging 画等值线图到底在解决什么手头有几十个采样点每个点带一个数值——土壤重金属含量、地下水位埋深、气温观测值。把这些点标在地图上一眼看过去全是孤立的数字看不出任何空间分布规律。真正想要的是那张能直接放进报告里的等值线图哪里浓度高、哪里低、梯度往哪个方向走一目了然。Kriging 就是干这个的——它不只是把点连成线而是先根据采样点的空间自相关性建一个变异函数模型再用这个模型去预测未采样位置的数值最后在预测网格上画等值线。和反距离加权插值相比Kriging 的优势在于它能给出预测方差也就是说你不仅知道某处预测值是多少还知道这个预测有多可靠。适合做地质、环境、水文、气象领域空间插值的人尤其是需要出正式图件、对插值精度有交代的场景。下面从工具选型一路讲到出图参数和踩坑记录。2. 工具链选型与变异函数拟合从 pykrige 到 gstools 怎么选2.1 为什么我最终选了 pykrige 而不是 gstoolsPython 生态里做 Kriging 插值主流就两个库pykrige 和 gstools。两个都能做普通克里金、泛克里金都支持各向异性但实际用下来差别不小。pykrige 的 API 更直白OrdinaryKriging类初始化时把坐标、值、变异函数模型、各向异性参数一次性传进去然后调execute方法指定网格分辨率就能拿到插值结果和方差。文档虽然不算特别详细但示例代码够用出图流程短。gstools 底层更灵活支持自定义变异函数、三维插值、条件模拟但学习曲线陡一些参数命名偏学术。我一般做二维等值线图就用 pykrige代码量少、调试快。如果要做三维插值或者需要条件模拟才切到 gstools。下面以 pykrige 为主线讲完整流程gstools 的差异在关键处会提。安装很简单pip install pykrige matplotlib numpypykrige 依赖 numpy 和 scipymatplotlib 用来出图。版本方面pykrige 1.7 以上对 numpy 2.x 兼容性好了很多如果遇到np.float报错大概率是 numpy 版本太新而 pykrige 太旧升级 pykrige 即可。2.2 变异函数模型怎么选球状、指数、高斯的使用边界Kriging 的核心是变异函数它描述的是“距离多远之后两个点的值就不再相关了”。pykrige 内置了三种常用模型模型适用场景关键参数球状spherical空间自相关随距离线性衰减到变程后归零变程、基台值、块金指数exponential自相关衰减较缓渐近归零变程、基台值、块金高斯gaussian自相关在短距离内变化平缓适合光滑曲面变程、基台值、块金选哪个模型不是拍脑袋。常见做法是先算实验变异函数看散点图的形状如果散点在小距离处快速上升然后很快走平球状模型拟合好如果上升缓慢、拖尾长指数模型更合适如果散点非常光滑、短距离内几乎不下降高斯模型可以考虑但高斯模型容易导致插值结果过度光滑在数据稀疏区域产生不真实的波动。pykrige 里指定模型用variogram_model参数from pykrige.ok import OrdinaryKriging import numpy as np # 假设 x, y 是经纬度或投影坐标z 是观测值 OK OrdinaryKriging( x, y, z, variogram_modelspherical, # 可选 spherical / exponential / gaussian variogram_parametersNone, # None 表示自动拟合 nlags6, # 实验变异函数的分组数 weightTrue # 拟合时按滞后距分组点数加权 )variogram_parameters传None时 pykrige 会自动拟合变程、基台值和块金。自动拟合在数据量足够一般大于 30 个点时效果尚可但数据少的时候拟合结果可能很离谱。我一般会先让它自动拟合把参数打印出来看看是否合理——变程不应该超过研究区最大距离的一半块金不应该大于基台值。如果自动拟合结果不合理就手动指定OK OrdinaryKriging( x, y, z, variogram_modelspherical, variogram_parameters{sill: 1.2, range: 5000, nugget: 0.1}, nlags6 )sill是基台值约等于数据方差range是变程单位跟坐标一致nugget是块金值反映测量误差和微尺度变异。手动调参时先把nugget设小一点比如方差的 5%然后调range看交叉验证误差。2.3 各向异性参数怎么设坐标旋转与各向异性比现实中的空间现象很少是各向同性的。比如河流沿岸的污染物扩散沿水流方向的相关距离可能比垂直方向大好几倍。pykrige 支持各向异性通过anisotropy_scaling和anisotropy_angle两个参数控制。anisotropy_scaling是各向异性比即长轴方向变程与短轴方向变程的比值。anisotropy_angle是长轴相对于正北方向的旋转角度度。设置方法OK OrdinaryKriging( x, y, z, variogram_modelexponential, anisotropy_scaling2.5, # 长轴变程是短轴的 2.5 倍 anisotropy_angle45, # 长轴方向为北偏东 45 度 nlags8 )这两个参数怎么定常见做法是先在各向同性条件下拟合变异函数然后在不同方向上分别计算实验变异函数看哪个方向的变程最长。那个方向就是长轴角度用np.arctan2算。各向异性比就是长轴变程除以短轴变程。如果拿不准先设anisotropy_scaling1.0跑一遍看看插值结果有没有明显的方向性拉伸再决定要不要调。3. 网格插值与等值线绘制从 execute 到 contourf 的完整链路3.1 生成插值网格分辨率与范围怎么定拿到变异函数模型之后下一步是定义要预测的网格。execute方法支持两种模式grid和masked。grid模式接收网格的 x 坐标数组和 y 坐标数组输出完整的二维矩阵masked模式额外接收一个掩膜只对掩膜内的区域插值适合研究区边界不规则的情况。import numpy as np # 定义网格范围一般比采样点范围略大一点 grid_x np.arange(x.min() - 500, x.max() 500, 100) # 步长 100 米 grid_y np.arange(y.min() - 500, y.max() 500, 100) z_pred, z_var OK.execute(grid, grid_x, grid_y)grid_x和grid_y的步长决定了等值线的光滑程度。步长太大等值线会有棱角步长太小计算量上去但视觉提升有限。我一般让网格步长约等于采样点平均间距的三分之一到五分之一。比如采样点平均间距 500 米步长就取 100 到 150 米。execute返回两个数组z_pred是预测值矩阵z_var是预测方差矩阵。z_pred直接用来画等值线z_var可以用来画不确定性图或者用来判断哪些区域的插值结果不可信。3.2 用 matplotlib 画等值线图contour 与 contourf 的取舍拿到z_pred之后画图本身不复杂但有几个参数直接影响出图质量。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 8)) # 填充等值线 cf ax.contourf(grid_x, grid_y, z_pred, levels15, cmapRdYlBu_r) # 叠加等值线 cs ax.contour(grid_x, grid_y, z_pred, levels15, colorsblack, linewidths0.5) # 标注等值线数值 ax.clabel(cs, inlineTrue, fontsize8, fmt%.1f) # 画采样点 ax.scatter(x, y, cblack, s20, markero, label采样点) # 颜色条 cbar fig.colorbar(cf, axax, shrink0.8) cbar.set_label(观测值) ax.set_xlabel(X (m)) ax.set_ylabel(Y (m)) ax.set_title(Kriging 插值等值线图) ax.legend() plt.tight_layout() plt.savefig(kriging_contour.png, dpi300) plt.show()levels15表示把数值范围分成 15 个区间。区间太少等值线稀疏细节丢失区间太多颜色条太密读图困难。我一般取 10 到 20 之间看数据范围决定。cmap选RdYlBu_r是因为它从蓝到红过渡自然适合表示浓度、温度这类有方向性的变量。如果数据有正负用RdBu_r以零为中心对称。clabel的inlineTrue会把等值线数值嵌在线条里比旁边放图例清爽。fmt%.1f控制小数位数根据数据量级调整。3.3 掩膜与边界裁剪只对研究区插值如果研究区不是矩形比如是一个流域边界或者行政区划直接画矩形网格会把区域外的插值结果也画出来看起来不专业。pykrige 的masked模式可以解决这个问题但需要先准备一个掩膜数组。from matplotlib.path import Path # 假设 boundary 是一个 Nx2 的数组表示研究区边界多边形 boundary np.array([[x1, y1], [x2, y2], ...]) # 生成网格点 grid_x np.arange(x.min() - 500, x.max() 500, 100) grid_y np.arange(y.min() - 500, y.max() 500, 100) xx, yy np.meshgrid(grid_x, grid_y) # 判断每个网格点是否在多边形内 points np.vstack((xx.ravel(), yy.ravel())).T path Path(boundary) mask path.contains_points(points).reshape(xx.shape) # 用 masked 模式插值 z_pred, z_var OK.execute(masked, grid_x, grid_y, maskmask)mask是一个布尔矩阵True表示该点参与插值False表示跳过。execute返回的z_pred在False位置是np.nan画图时 matplotlib 会自动留白。掩膜边界如果很复杂contains_points可能比较慢但一般几千个网格点也就几秒钟的事。如果边界特别复杂可以考虑用shapely做预处理把边界简化一下再传进去。4. 避坑与排查Kriging 插值翻车实录4.1 插值结果全是 NaN 或者异常值现象execute返回的z_pred里大量np.nan或者出现远超数据范围的极端值。原因最常见的是坐标数组里有重复点或者变异函数拟合失败导致参数为nan。另外如果variogram_parameters手动指定时range设得过大协方差矩阵可能接近奇异求解时数值不稳定。解决先检查x、y是否有重复坐标用np.unique去重。然后打印OK.variogram_model_parameters看拟合结果是否合理。如果自动拟合失败改用手动参数把range设为研究区最大距离的三分之一左右nugget设为方差的 5% 到 10%。4.2 等值线图出现“牛眼”状同心圆现象每个采样点周围都有一圈圈密集的同心等值线看起来像牛眼。原因变异函数的nugget设得太小或者range设得太小导致插值结果过度依赖最近邻点远处点的权重衰减太快。解决增大nugget或range。常见做法是先把nugget调到方差的 10% 到 20%然后逐步增大range直到牛眼消失。如果数据本身噪声大nugget大一点反而更合理。4.3 各向异性参数设反导致方向性错误现象设置了anisotropy_scaling和anisotropy_angle之后插值结果的方向性跟预期相反。原因anisotropy_angle的定义是长轴相对于正北的顺时针角度但很多人会把它当成相对于正东的角度或者把长短轴搞反。解决先用anisotropy_scaling1.0跑一遍确认各向同性结果。然后分别计算 0°、45°、90°、135° 四个方向的实验变异函数看哪个方向变程最长。那个方向就是长轴。anisotropy_angle用np.degrees(np.arctan2(dy, dx))算注意 pykrige 的角度定义是北偏东为正。4.4 网格步长过粗导致等值线锯齿现象等值线图看起来一格一格的不光滑。原因grid_x和grid_y的步长太大插值网格分辨率不够。解决把步长缩小到采样点平均间距的三分之一以下。但步长太小会导致计算量指数增长一般 100 到 200 米步长对区域尺度研究够用了。如果还嫌不够光滑可以在画图时用scipy.ndimage.zoom对z_pred做插值放大但注意这只是视觉平滑不增加信息量。4.5 交叉验证误差大但不知道哪里出了问题现象留一交叉验证的均方根误差很大但不知道是变异函数模型选错了还是数据本身有问题。解决pykrige 提供了OK.cross_validate()方法返回预测值、方差和误差。先看误差的空间分布——如果误差集中在某个区域可能是该区域采样点太少或者存在趋势面。如果误差跟观测值大小相关可能需要先做数据变换比如对数变换再插值。常见做法是对偏态分布的数据先取对数插值完再指数还原。5. 进阶技巧用预测方差图判断插值可信度5.1 预测方差图怎么读execute返回的z_var是预测方差它反映的是插值结果的不确定性。在采样点附近方差接近零远离采样点方差逐渐增大。把z_var画成等值线图可以直观看到哪些区域的插值结果可信、哪些区域是“猜”出来的。fig, axes plt.subplots(1, 2, figsize(16, 6)) # 左图预测值等值线 cf1 axes[0].contourf(grid_x, grid_y, z_pred, levels15, cmapRdYlBu_r) axes[0].scatter(x, y, cblack, s20) axes[0].set_title(Kriging 预测值) fig.colorbar(cf1, axaxes[0], shrink0.8) # 右图预测方差等值线 cf2 axes[1].contourf(grid_x, grid_y, z_var, levels15, cmapGreys) axes[1].scatter(x, y, cred, s20) axes[1].set_title(Kriging 预测方差) fig.colorbar(cf2, axaxes[1], shrink0.8) plt.tight_layout() plt.savefig(kriging_with_variance.png, dpi300)左图看趋势右图看可信度。如果右图某区域方差特别大而左图那个区域又有明显的等值线密集区那就要小心了——那个地方的“规律”可能是插值算法编出来的不是数据支持的。5.2 用方差图指导补充采样预测方差图还有一个实用场景决定下一步去哪里补采样。方差大的区域就是信息最缺乏的区域优先去那里采样能最大程度降低整体不确定性。我一般会把方差图叠加到采样点分布图上圈出方差最大的三个区域作为下一轮采样的候选位置。5.3 从教训到习惯早期做 Kriging 插值我只盯着预测值图看觉得等值线画出来漂亮就行。直到有一次评审会上被问“这片高值区有多少采样点支撑”翻出方差图一看那块区域方差大得离谱等值线完全是外推出来的。从那以后我每次出 Kriging 图都强制走一遍“预测值 方差”双图流程方差大的区域要么补采样要么在报告里明确标注不确定性。这个习惯帮我避开了好几次潜在的误判。希望帮到你。本文还有配套的精品资源点击获取