geostatistics-interpolation

地统计学与插值

目的：进行插值，并报告未知之处。专业产品与漂亮栅格的区别在于不确定性表面和诚实的交叉验证——两者在这里都不可或缺。

方法选择

情况	方法
现象密集且平滑，只需快速查看	IDW（报告幂参数；测试 1-3）
存在空间结构的物理现象，且需要不确定性	普通克里金（默认的专业选择）
存在明显趋势（高程梯度、海岸效应）	对协变量进行泛克里金或回归克里金
有明显的协变量可用（DEM、土地覆盖、距离）	回归克里金 / 随机森林残差克里金
目标为类别	指示克里金
需要无蜂窝状伪影的镶嵌，且不希望外推	自然邻域

IDW 是合理的基线方法，但没有误差模型，且会在极值周围产生靶心效应；仅交付 IDW 产品时，请说明这一点。

探索阶段（任何插值之前）

• 绘制带有数值的点；查找坐标完全相同的重复点（对其取平均或进行偏移，否则克里金矩阵会变成奇异矩阵）。
• 直方图 + 偏度：偏度明显的变量（降雨量、浓度）通常需要进行对数变换 / 正态分数变换；正确地对预测结果进行反变换（对数正态克里金需进行偏差校正）。
• 趋势检查：对数值与 x、y 及候选协变量进行回归；如果趋势明显，则采用泛克里金 / 回归克里金路径。
• 如果采样具有偏向性（高值区域采样密集），则进行去聚类。

变异函数规范

变异函数是一个建模决策，而不是自动拟合的结果：

import gstools as gs

bin_center, gamma = gs.vario_estimate((x, y), values, max_dist=dmax)  # dmax ≈ half extent
model = gs.Exponential(dim=2)
model.fit_variogram(bin_center, gamma, nugget=True)
print(model)   # report: nugget, sill, range — plus the fitted plot

• 最大滞后距离约为区域直径的一半；每个分箱至少包含 30 对点。
• 使用方向变异函数（0/45/90/135°）检查各向异性；地质和气象场通常具有各向异性——应拟合各向异性模型，而不是忽略它。
• 用文字解读并报告参数：块金（测量误差 + 微尺度方差）、变程（相关距离）、基台。块金接近基台，意味着数据几乎无法支持插值——请如实说明。
• 绝不可对超出最近样本变异函数变程的区域进行有意义的插值；应屏蔽或标记这些单元格。

克里金执行

from pykrige.ok import OrdinaryKriging

ok = OrdinaryKriging(x, y, values, variogram_model="exponential",
                     variogram_parameters={"sill": s, "range": r, "nugget": n},
                     coordinates_type="euclidean")
z, ss = ok.execute("grid", gridx, gridy)   # ss = kriging VARIANCE — keep it!

• 使用投影坐标参考系（除非明确使用支持地理坐标的模型，否则以度为单位的克里金距离是错误的）。
• 交付两张栅格：预测值和克里金标准差（np.sqrt(ss)）。不确定性图可用于确定下一步的采样位置，以及哪些区域不应信任地图结果。
• 搜索邻域：通常使用最近的 16-32 个点；记录并说明所用设置。

验证——留一法及其他方法

• 小样本量时使用 LOOCV：报告 ME（偏差约为 0？）、RMSE 和标准化 RMSE（如果克里金方差可信，该值应约为 1）。
• 对于聚类样本，LOOCV 会给出过于乐观的结果；还应运行空间分块交叉验证（标准方案：ml-experiment-standards → references/spatial-cv-protocol.md），并报告两者。
• 绘制观测值与预测值的散点图，并添加 1:1 线；绘制交叉验证残差图——残差在空间上呈聚类分布，意味着遗漏了趋势 / 协变量。
• 与简单基线（全局均值、IDW）进行比较：克里金必须证明其复杂性物有所值。

报告模板

## Surface: <variable>
- n = <>, extent, CRS, transform applied: <log/none>
- Variogram: <model>, nugget/sill/range = <...>, anisotropy: <...>
- Method: <OK/UK/RK + covariates>
- LOOCV: ME <>, RMSE <>, RMSE_std <>; block-CV RMSE <>
- Uncertainty: kriging SD raster delivered; area beyond reliable range masked

常见问题检查清单

• 未查看图表就接受自动拟合的变异函数。
• 在大范围区域使用 EPSG:4326 进行克里金。
• 交付预测栅格，却未提供对应的方差 / 标准差栅格。
• 使用简单的 exp() 对数克里金结果进行反变换（会产生偏差！）。
• 将远超数据凸包的外推区域呈现为具有同等置信度。
• 重复的站点坐标导致拟合崩溃，或在未提示的情况下扭曲拟合结果。

执行约定

• 工作流： 检查采样和变换需求；建模空间依赖性；拟合候选曲面；根据简单基线进行验证；量化不确定性；屏蔽不受样本支持的外推区域。
• 决策规则： 选择 IDW 或趋势方法作为透明的基线；存在可靠的变异函数时使用克里金法；协变量能改善分块验证结果时使用回归克里金法。
• 验证协议： 检查变异函数诊断、留一交叉验证（LOOCV）和空间分块交叉验证（空间 block-CV）的残差、相对于基线的提升、不确定性校准，以及超出样本凸包后的表现。
• 失败模式： 当采样过于稀疏、缺乏依赖性、重复数据占主导、距离单位无效或不确定性未经校准时，不应给出置信度高的曲面。
• 交付内容： 预测曲面、不确定性曲面、变异函数及参数、验证表、采样支持掩码和方法局限。
• 源信息时效性： 依赖软件包 API 或默认值前，请查阅权威来源注册表，并记录查阅日期。