简介：这份资源围绕 Kriging 空间插值方法展开，提供一套可直接运行的等值线图绘制源码，面向地理信息、地质勘探及空间数据分析方向的学习者与开发者，帮助解决从散点观测数据到连续空间分布可视化的实现问题。压缩包共 84 个文件，以 h 头文件与 cpp 源文件为主体，辅以 obj 编译产物、ico 与 cur 图标资源、txt 测试数据及 exe 可执行程序，整体约 1.77MB，工程结构完整，便于直接编译调试。内容涵盖半方差函数构建、Ordinary Kriging 与 Universal Kriging 模型选择、矩阵求解与权重分配、插值预测及等值线绘制等关键环节，并附带测试数据文件，方便读者对照验证插值结果的合理性与等值线的连续性。目前已有 1214 人学习下载，适合希望深入理解 Kriging 原理并动手实践等值图绘制的读者参考。
1. Kriging 画等值线图：从散点到连续曲面的最小闭环
手里攥着几十个采样点，领导要一张能直接放进报告的等值线图，这大概是地信、水文、环境监测方向最常见的场景。直接拿点数据插值成栅格再画线，听起来简单，但真正动手就会发现：插值方法选错，等值线会扭曲成鬼画符；参数没调好，图上一片空白或者全是同心圆。Kriging 画等值线图之所以被反复提起，是因为它不只是一个插值算法，而是一套带统计假设的空间预测框架——它能在给出预测值的同时告诉你哪里可信、哪里是外推。这篇文章面向的是需要把离散观测点变成连续等值线图的从业者，无论你用的是 Python 的 PyKrige、R 的 gstat，还是 ArcGIS 里的 Geostatistical Analyst，底层的变差函数拟合、块金效应处理、搜索邻域设定逻辑是相通的。我会按“数据准备 → 变差函数建模 → 插值成栅格 → 等值线提取与出图”这条链路，把每一步的参数含义、常见翻车点和可抄作业的代码写清楚，让你拿到自己的数据就能跑通第一版。
2. 数据准备与 Kriging 的统计前提：为什么你的点不能直接扔进去
2.1 采样点分布与投影：两个最容易被忽略的预处理
Kriging 的第一条血泪经验：数据没整理好，后面全是玄学。很多人拿到经纬度坐标就直接丢进插值函数，结果等值线图在投影变形下被拉长或压扁，尤其是跨纬度较大的区域。常见做法是先把地理坐标（WGS84）投影到适合研究区的平面坐标系，比如 UTM 或高斯-克吕格，单位统一成米。这样变差函数拟合时的距离才有物理意义，搜索邻域半径也才好设定。
第二个坑是重复点和异常值。同一个位置有多个采样值时，Kriging 协方差矩阵会出现奇异，表现为插值结果全是 NaN 或者报“singular matrix”。处理方式很简单：按坐标分组取均值，或者加一个极小的随机扰动。异常值则要在变差函数拟合前识别，因为 Kriging 对离群点非常敏感，一个极端高值能把块金效应拉得很大，导致整个曲面过度平滑。
import numpy as np import pandas as pd from pyproj import Transformer # 假设原始数据有三列：lon, lat, value df = pd.read_csv("samples.csv") # 1. 投影转换：WGS84 -> UTM zone 50N（根据研究区调整） transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650", always_xy=True) x, y = transformer.transform(df["lon"].values, df["lat"].values) df["x"], df["y"] = x, y # 2. 重复点合并：按 1 米网格取均值 df["x_round"] = df["x"].round(0) df["y_round"] = df["y"].round(0) df = df.groupby(["x_round", "y_round"], as_index=False).agg( x=("x", "mean"), y=("y", "mean"), value=("value", "mean") ) # 3. 异常值粗筛：3 倍标准差之外先标记，不急着删 mu, sigma = df["value"].mean(), df["value"].std() df["is_outlier"] = np.abs(df["value"] - mu) > 3 * sigma print(f"标记异常点 {df['is_outlier'].sum()} 个，建议逐个人工确认")这段代码做了三件事：投影转换把经纬度变成平面米坐标，重复点合并避免协方差矩阵奇异，异常值标记留出人工判断余地。参数上，UTM 带号要根据研究区经度选，比如东经 117° 附近用 zone 50N；round(0)是 1 米容差，如果采样密度低可以放宽到 10 米。注意异常值不要用代码直接删，先标记出来对照原始记录，有些高值恰恰是真实污染热点，删了就丢掉了关键信息。
2.2 变差函数拟合：Kriging 的“黑匣子”其实有旋钮
变差函数（variogram）描述的是“距离越远，两点差异越大”这件事。它有三个核心参数：块金值（nugget）、基台值（sill）、变程（range）。块金值代表小于最小采样间距的随机误差，基台值是变差函数趋于平稳时的值，变程是空间自相关消失的距离。这三个参数直接决定等值线图的平滑程度和细节保留能力。
拟合变差函数时，常见做法是先算实验变差函数，再用球状模型、指数模型或高斯模型去拟合。球状模型最常用，因为它有明确的变程，超过变程后协方差为零，计算稳定。指数模型和高斯模型在变程附近是渐近的，适合连续变化的地质变量。我一般会先试球状模型，如果拟合残差大再换指数模型。
from pykrige.ok import OrdinaryKriging import matplotlib.pyplot as plt # 准备数据 data = df[~df["is_outlier"]] x = data["x"].values y = data["y"].values z = data["value"].values # 构建 Ordinary Kriging 对象，先不指定变差函数参数，让 PyKrige 自动拟合 OK = OrdinaryKriging( x, y, z, variogram_model="spherical", # 球状模型 variogram_parameters=None, # None 表示自动拟合 nlags=6, # 实验变差函数的分组数 weight=True, # 拟合时按点对数加权 verbose=False, enable_plotting=True # 画出拟合图，方便检查 ) # 打印自动拟合出的参数 print("块金值:", OK.variogram_model_parameters[0]) print("基台值:", OK.variogram_model_parameters[1]) print("变程:", OK.variogram_model_parameters[2])nlags=6是实验变差函数的分组数，点少的时候可以设 4 到 6，点多了可以到 10 到 15。weight=True让点对多的滞后距在拟合时权重更大，通常能得到更稳定的结果。enable_plotting=True会弹出一张拟合图，这张图一定要看：如果实验点散乱无规律，或者拟合曲线完全压不住点，说明数据本身空间结构性弱，强行 Kriging 只会得到一张没有信息量的平滑图。这时候要么增加采样点，要么换回 IDW 之类的确定性插值。
2.3 搜索邻域与各向异性：让等值线不出现“同心圆”
默认情况下，OrdinaryKriging 会用全部点参与每个预测点的计算。点少的时候没问题，点多了计算量爆炸，而且远处点对局部细节的干扰会让等值线出现不合理的“同心圆”或“牛眼”。常见做法是设置搜索邻域：只选距离预测点一定范围内的若干点参与计算。
# 设置搜索邻域：最多 12 个点，搜索半径 5000 米 OK = OrdinaryKriging( x, y, z, variogram_model="spherical", variogram_parameters=OK.variogram_model_parameters, # 沿用上一步拟合结果 nlags=6, weight=True, verbose=False, enable_plotting=False ) # 生成预测网格 grid_x = np.linspace(x.min(), x.max(), 200) grid_y = np.linspace(y.min(), y.max(), 200) z_pred, z_var = OK.execute( "grid", grid_x, grid_y, n_closest_points=12, # 只选最近的 12 个点 backend="C" # 用 C 后端加速 )n_closest_points=12是经验值，点密集时可以降到 8，点稀疏时升到 16。backend="C"比默认的 Python 后端快一个数量级，网格大时必开。如果研究区有明显的各向异性（比如沿河流方向污染扩散更远），需要指定anisotropy_scaling和anisotropy_angle，否则等值线会被拉成不自然的椭圆。各向异性参数一般从变差函数云图里目视判断，或者用 gstat 的fit.variogram自动估计。
3. 从插值栅格到等值线图：提取、平滑与出图参数
3.1 等值线提取：matplotlib 与 GDAL 两条路
插值完成后，你得到的是一个二维数组z_pred，对应网格grid_x和grid_y。画等值线图有两种常见路径：一是直接用 matplotlib 的contour或contourf，适合快速出图；二是写成 GeoTIFF 再用 GDAL 或 QGIS 提取等值线，适合需要矢量等值线做后续分析的场景。
import matplotlib.pyplot as plt import numpy as np # 方式一：matplotlib 直接画等值线 fig, ax = plt.subplots(figsize=(10, 8)) cs = ax.contourf(grid_x, grid_y, z_pred, levels=15, cmap="RdYlBu_r") ax.contour(grid_x, grid_y, z_pred, levels=15, colors="black", linewidths=0.5) cbar = fig.colorbar(cs, ax=ax, shrink=0.8) cbar.set_label("预测值") ax.scatter(x, y, c="black", s=10, label="采样点") ax.set_aspect("equal") ax.legend() plt.tight_layout() plt.savefig("kriging_contour.png", dpi=300)levels=15控制等值线的密集程度，太多会糊成一团，太少则丢失细节。cmap="RdYlBu_r"是红黄蓝反向色带，适合表示高低值。ax.set_aspect("equal")保证 x 和 y 轴比例一致，否则等值线会被拉伸。如果要把等值线导出成 shapefile，用 GDAL 的ContourGenerate函数，指定contourInterval和fixedLevels，注意 GDAL 的等值线是栅格像元追踪出来的，会有锯齿，需要后续在 QGIS 里平滑。
3.2 平滑与裁剪：让图能直接进报告
Kriging 插值结果本身带有预测方差，z_var数组可以单独画一张方差图，标出哪些区域外推严重。但报告里通常只需要一张干净的等值线图。常见需求是：按研究区边界裁剪、对等值线做轻微平滑、加上比例尺和指北针。
裁剪用rasterio的mask函数，传入边界矢量文件即可。平滑用scipy.ndimage.gaussian_filter，sigma取 1 到 2，太大会把真实极值抹掉。指北针和比例尺在 matplotlib 里没有原生支持，一般用matplotlib_scalebar库或者手动加箭头和文字。注意平滑要在裁剪之前做，否则边界处的插值噪声会被平滑成假异常。
from scipy.ndimage import gaussian_filter import rasterio from rasterio.mask import mask import geopandas as gpd # 平滑：sigma=1.5，轻微去噪 z_smooth = gaussian_filter(z_pred, sigma=1.5) # 裁剪：读取研究区边界 boundary = gpd.read_file("study_area.shp") # 将平滑后的栅格写成临时 GeoTIFF 再裁剪 transform = rasterio.transform.from_origin( grid_x.min(), grid_y.max(), (grid_x.max() - grid_x.min()) / (len(grid_x) - 1), (grid_y.max() - grid_y.min()) / (len(grid_y) - 1) ) with rasterio.open( "temp.tif", "w", driver="GTiff", height=z_smooth.shape[0], width=z_smooth.shape[1], count=1, dtype=z_smooth.dtype, crs="EPSG:32650", transform=transform ) as dst: dst.write(z_smooth, 1) with rasterio.open("temp.tif") as src: out_image, out_transform = mask(src, boundary.geometry, crop=True) # out_image 就是裁剪后的数组，可以继续用 contour 画sigma=1.5是经验值，数据噪声大可以到 2，数据本身平滑就设 0.5 到 1。裁剪时注意边界矢量的坐标系要和栅格一致，否则mask会报错或裁出空白。如果研究区边界复杂，裁剪后等值线在边界处会断开，这是正常的，不要强行连接。
4. Kriging 画等值线图的避坑与排查：5 个真实翻车现场
4.1 现象：等值线图全是同心圆，采样点被完全忽略
原因：变差函数拟合出的变程远大于研究区范围，导致所有点对预测点的影响权重几乎相同，Kriging 退化成全局均值插值。常见于采样点少且分布稀疏的情况。
解决：检查变差函数拟合图，如果实验变差函数一直没有趋于平稳，说明空间自相关尺度超出采样范围。此时要么缩小研究区，要么手动指定一个合理的变程（比如研究区对角线的三分之一），强制 Kriging 关注局部变化。
4.2 现象：插值结果出现负值，但原始数据全是正数
原因：Ordinary Kriging 在数据不服从正态分布时，预测值可能超出原始数据范围。这是算法本身的特性，不是 bug。
解决：如果物理上不允许负值，插值后做截断（z_pred[z_pred < 0] = 0），或者先对数据做对数变换，插值后再指数还原。对数变换还能缓解异常值影响，一举两得。
4.3 现象：等值线在采样点密集区出现“牛眼”状闭合圈
原因：搜索邻域设得太小，或者块金值拟合过低，导致每个采样点都被当成精确值，Kriging 过度拟合局部噪声。
解决：增大n_closest_points，或者手动调高块金值。块金值代表测量误差和微尺度变异，适当调高能让曲面更平滑。一般块金值占基台值的 10% 到 30% 比较合理，低于 5% 就要警惕过拟合。
4.4 现象：execute报错 “singular matrix” 或结果全为 NaN
原因：采样点中有完全重复的坐标，或者搜索邻域内的点共线。共线在 Ordinary Kriging 的矩阵求解中会导致奇异。
解决：重复点按 2.1 节的方法合并。共线问题可以通过增大搜索半径或增加n_closest_points来缓解，让邻域内包含非共线的点。如果仍然报错，改用UniversalKriging并指定漂移项，或者换用OrdinaryKriging的pseudo_inv=True参数用伪逆求解。
4.5 现象：等值线图在投影边缘严重变形
原因：用了地理坐标系（经纬度）直接插值，距离计算用的是度数，不同纬度下 1 度对应的实际距离差异巨大。
解决：回到 2.1 节，先投影再插值。如果研究区跨多个 UTM 带，用 Albers 等面积投影或 Lambert 等角圆锥投影，保证整个区域的距离变形可控。
5. 交叉验证与参数调优：让等值线图从“能看”到“可信”
Kriging 相比 IDW 的核心优势是能给出预测方差，但方差图只能告诉你相对不确定性，不能直接判断模型好坏。真正靠谱的做法是留一法交叉验证：每次剔除一个采样点，用剩余点预测该点，统计预测值与实测值的误差。PyKrige 没有内置交叉验证函数，但手动实现很简单。
from sklearn.metrics import mean_squared_error, r2_score def loo_cv(x, y, z, variogram_model="spherical", n_closest=12): preds = [] for i in range(len(x)): mask = np.ones(len(x), dtype=bool) mask[i] = False OK = OrdinaryKriging( x[mask], y[mask], z[mask], variogram_model=variogram_model, nlags=6, weight=True, verbose=False ) pred, _ = OK.execute("points", np.array([x[i]]), np.array([y[i]]), n_closest_points=n_closest) preds.append(pred[0]) preds = np.array(preds) rmse = np.sqrt(mean_squared_error(z, preds)) r2 = r2_score(z, preds) return rmse, r2, preds rmse, r2, preds = loo_cv(x, y, z) print(f"LOO-RMSE: {rmse:.3f}, R2: {r2:.3f}")n_closest=12要和最终出图时保持一致，否则交叉验证结果不能代表最终模型。RMSE 越小越好，R2 越接近 1 越好。如果 R2 低于 0.5，说明空间自相关弱，Kriging 可能不是最佳选择，可以对比 IDW 或样条插值的交叉验证结果。我一般会跑三组参数：球状模型、指数模型、高斯模型，各做一次 LOO 交叉验证，选 RMSE 最小的那个。注意交叉验证的 RMSE 不能直接和原始数据的标准差比，因为 Kriging 平滑了极值，RMSE 通常会小于标准差，这是正常的。
调参时还有一个容易忽略的点：nlags和weight的组合。nlags太少，实验变差函数不稳定；太多，计算量大且噪声大。weight=True在点对分布不均时能改善拟合，但如果点对在各滞后距上分布均匀，weight=False反而更无偏。我的习惯是先用weight=True跑一遍，看拟合图，如果曲线明显被短距离点对拉偏，再关掉权重。
最后说一个验证等值线图合理性的土办法：把等值线图叠在采样点符号图上，看高值区是否落在高值点附近，低值区是否避开高值点。如果一张等值线图的高值中心离所有高值采样点都很远，那大概率是变差函数参数出了问题，回去检查变程和块金值。这个办法虽然不严谨，但能快速筛掉明显翻车的图。希望帮到你。
本文还有配套的精品资源，点击获取
