scikit-gstat

SciKit-GStat - 地统计学

快速参考

import skgstat as skg
import numpy as np

# Create variogram
V = skg.Variogram(coordinates=coords, values=values, n_lags=15)

# Fit model
V.model = 'spherical'
print(f"Range: {V.parameters[0]:.2f}, Sill: {V.parameters[1]:.2f}")

# Kriging interpolation
ok = skg.OrdinaryKriging(V)
predictions = ok.transform(grid_coords)

关键类

类	用途
Variogram	经验变差函数和理论变差函数
OrdinaryKriging	基于空间相关性的插值
DirectionalVariogram	各向异性变差函数
SpaceTimeVariogram	时空分析

基本操作

创建并拟合变差函数

import skgstat as skg

V = skg.Variogram(
    coordinates=coords,      # (n, 2) array of x, y
    values=values,           # (n,) array of measurements
    n_lags=15,
    maxlag='median'          # or specific distance
)

# Fit model: 'spherical', 'exponential', 'gaussian', 'matern', 'stable'
V.model = 'spherical'

# Get parameters
print(f"Range: {V.parameters[0]:.2f}")
print(f"Sill: {V.parameters[1]:.2f}")
print(f"Nugget: {V.parameters[2]:.2f}")
print(f"RMSE: {V.rmse:.4f}")

普通克里金

import skgstat as skg
import numpy as np

V = skg.Variogram(coords, values, model='spherical')
ok = skg.OrdinaryKriging(V)

# Create prediction grid
x = np.linspace(0, 100, 50)
y = np.linspace(0, 100, 50)
xx, yy = np.meshgrid(x, y)
grid_coords = np.column_stack([xx.ravel(), yy.ravel()])

# Predict
predictions = ok.transform(grid_coords)
Z = predictions.reshape(xx.shape)

# Get variance
ok.return_variance = True
predictions, variance = ok.transform(grid_coords)

方向变差函数

import skgstat as skg

DV = skg.DirectionalVariogram(
    coordinates=coords,
    values=values,
    azimuth=45,          # Direction in degrees
    tolerance=22.5,      # Angular tolerance
    bandwidth='q33'      # Perpendicular bandwidth
)

# Check anisotropy
for az in [0, 45, 90, 135]:
    DV.azimuth = az
    print(f"Azimuth {az}: Range = {DV.parameters[0]:.2f}")

交叉验证

import skgstat as skg
from sklearn.model_selection import cross_val_score

V = skg.Variogram(coords, values, model='spherical')
ok = skg.OrdinaryKriging(V)

scores = cross_val_score(ok, coords, values, cv=5, scoring='neg_mean_squared_error')
print(f"CV RMSE: {np.sqrt(-scores.mean()):.4f}")

稳健估计量

import skgstat as skg

# Use robust estimator for noisy data
V = skg.Variogram(
    coords, values,
    estimator='cressie'  # 'matheron', 'cressie', 'dowd', 'genton'
)

模型快速参考

模型	特性
spherical	最常见，原点附近呈线性
exponential	不会达到基台，逐渐趋近
gaussian	原点附近呈抛物线状，平滑
matern	可灵活控制平滑度

与替代方案的适用场景对比

使用场景	工具	原因
变差函数分析和克里金	scikit-gstat	现代 API，兼容 sklearn
GSLIB 风格模拟（SGSIM）	GeostatsPy	完整的 GSLIB 模拟引擎
带趋势/漂移的克里金	pykrige	泛克里金、回归克里金
随机场生成	gstools	灵活的协方差、SRF 生成
时空变差函数	scikit-gstat	内置 SpaceTimeVariogram
生产级地质建模	SGeMS / Petrel	GUI、大规模 3D 模型
稳健变差函数估计	scikit-gstat	Cressie、Dowd、Genton 估计量
集成到 ML 流水线	scikit-gstat	sklearn fit/transform 接口

在以下情况下选择 scikit-gstat：你需要用于变差函数拟合和克里金的 Pythonic、兼容 scikit-learn 的 API。最适合结合交叉验证并集成到 ML 流水线中的探索性地统计分析。

在以下情况下选择 GeostatsPy：你需要兼容 GSLIB 的模拟工作流（SGSIM、SISIM），或正在使用传统地统计约定。

在以下情况下选择 pykrige：你需要使用外部漂移变量进行泛克里金，或需要结合地统计学与机器学习预测的回归克里金。

常见工作流

变差函数拟合和普通克里金

• 将空间数据作为 numpy 数组加载（坐标和值）
• 使用适当的 n_lags 和 maxlag 创建 Variogram 对象
• 测试估计量：matheron（默认）与 cressie（稳健，适用于含噪数据）
• 拟合多个模型（spherical、exponential、gaussian）并比较 RMSE
• 使用 0、45、90、135 度的 DirectionalVariogram 检查各向异性
• 根据 RMSE 和视觉拟合效果选择最佳模型
• 根据拟合后的变差函数创建 OrdinaryKriging 对象
• 定义预测网格并运行 ok.transform(grid_coords)
• 设置 ok.return_variance = True 以获取克里金方差
• 使用 cross_val_score() 进行交叉验证，评估预测质量
• 绘制预测值和克里金方差图

常见问题

问题	解决方案
变差函数平坦或不规则	调整 n_lags 和 maxlag（尝试 maxlag='median'）
模型拟合不佳（RMSE 较高）	尝试不同的模型类型或嵌套结构
克里金速度太慢	减少条件点数量或降低网格分辨率
块金效应过大	可能表示测量误差；尝试稳健估计量
各向异性不明确	在 DirectionalVariogram 中使用较小的角度容差

提示

1. Maxlag 应约为研究区域对角线的 50%
2. 对含噪数据使用稳健估计量（cressie、dowd）
3. 测试多个模型并比较 RMSE
4. 在克里金之前，使用方向变差函数检查各向异性

参考资料

• 变差函数模型 - 模型方程和参数
• 克里金方法 - 克里金类型和配置

脚本

• scripts/variogram_analysis.py - 完整的变差函数分析工作流