geostatspy

GeostatsPy - 地质统计分析

快速参考

import geostatspy.GSLIB as GSLIB
import geostatspy.geostats as geostats
import pandas as pd

df = pd.read_csv('data.csv')
df['npor'], tvpor, tnspor = geostats.nscore(df, 'porosity')  # Transform

lag, gamma, npairs = geostats.gamv(df, 'X', 'Y', 'npor',      # Variogram
    tmin=-9999, tmax=9999, xlag=50, xltol=25, nlag=15,
    azm=0, atol=22.5, bandwh=9999, bandwd=9999)

vario = GSLIB.make_variogram(nug=0.0, nst=1, it1=1, cc1=1.0,  # Model
                              azi1=0, hmaj1=300, hmin1=300)

est, var = geostats.kb2d(df, 'X', 'Y', 'npor', ..., vario=vario)  # Krige

关键函数

类别	函数
可视化	locmap、pixelplt、hist
变异函数	gamv、vmodel
克里金	kb2d、kb3d
模拟	sgsim、sisim
变换	nscore、backtr
去聚集	declus

常见操作

1. 正态分数变换

df['npor'], tvpor, tnspor = geostats.nscore(df, 'porosity')
original = geostats.backtr(nscore_data, tvpor, tnspor, zmin=0, zmax=0.3)

2. 实验变异函数

lag, gamma, npairs = geostats.gamv(
    df, 'X', 'Y', 'npor',
    tmin=-9999, tmax=9999,    # Trimming limits
    xlag=50, xltol=25,        # Lag distance, tolerance
    nlag=15, azm=0, atol=22.5, bandwh=9999, bandwd=9999)

3. 变异函数模型

# Types: 1=spherical, 2=exponential, 3=gaussian
vario = GSLIB.make_variogram(
    nug=0.0, nst=1,            # Nugget, number of structures
    it1=1, cc1=1.0,            # Type, sill contribution
    azi1=0, hmaj1=300, hmin1=300)  # Azimuth, major/minor range

4. 克里金（kb2d）

est, var = geostats.kb2d(
    df, 'X', 'Y', 'npor', tmin=-9999, tmax=9999,
    nx=50, xmn=25, xsiz=50,    # Grid X: ncells, origin, size
    ny=50, ymn=25, ysiz=50,    # Grid Y
    nxdis=1, nydis=1, ndmin=1, ndmax=10,
    radius=500, ktype=0, skmean=0.0, vario=vario)  # ktype: 0=simple, 1=ordinary

5. 序贯高斯模拟

sim = geostats.sgsim(
    df, 'X', 'Y', 'npor', wcol=-1, scol=-1,
    tmin=-9999, tmax=9999, itrans=0,
    ismooth=0, dession=0, dmession=0,
    zmin=-4, zmax=4, ltail=1, ltpar=0, utail=1, utpar=0,
    nsim=1, nx=50, xmn=25, xsiz=50, ny=50, ymn=25, ysiz=50,
    nz=1, zmn=0, zsiz=1, seed=73073,
    ndmin=1, ndmax=10, nodmax=10, radius=500, radius1=500,
    sang1=0, sang2=0, sang3=0, mxctx=10, mxcty=10, mxctz=1,
    ktype=0, vario=vario)

6. 去聚集

wts, cell_size, ncut = geostats.declus(
    df, 'X', 'Y', 'porosity', iminmax=1, noff=10, ncell=20, cmin=10, cmax=500)
declustered_mean = np.average(df['porosity'], weights=wts)

变异函数模型

代码	模型	用途
1	球状模型	最常见，有限变程
2	指数模型	渐近达到基台值
3	高斯模型	非常平滑，在原点附近呈抛物线状
4	幂函数模型	无界，类似分形

关键参数

参数	描述
nug	块金效应（测量误差 + 微尺度变化）
sill	总方差（块金 + 结构贡献）
range	相关性变得可忽略的距离
azimuth	连续性最大的方向（相对北方的角度）
ktype	0=简单克里金（均值已知），1=普通克里金

适用场景与替代方案

使用场景	工具	原因
GSLIB 风格的工作流	GeostatsPy	将 GSLIB 程序直接移植到 Python
SGSIM / SISIM 模拟	GeostatsPy	完整的 GSLIB 模拟引擎
空间数据去聚集	GeostatsPy	内置 declus 函数
现代变异函数 API	scikit-gstat	API 更简洁，可集成 sklearn
仅进行克里金（不模拟）	pykrige	专注的 API，支持泛克里金
随机场生成	gstools	灵活的协方差模型和场生成
大规模 3D 地质建模	SGeMS / Petrel	基于 GUI，适用于工业工作流
指示变量模拟	GeostatsPy（sisim）	分类属性模拟

选择 GeostatsPy 的情况：需要在 Python 中使用兼容 GSLIB 的工作流，尤其是进行 序贯模拟（SGSIM/SISIM）、去聚集，或熟悉 GSLIB 参数约定时。最适用于储层表征工作流。

选择 scikit-gstat 的情况：偏好现代 scikit-learn 风格的变异函数分析与克里金 API， 并希望更好地集成到 Python 数据科学工作流中。

选择 pykrige 的情况：只需要克里金插值（不进行模拟），并且需要支持外部漂移或回归克里金的 泛克里金功能。

常见工作流

变异函数分析与克里金插值

• 将空间数据加载到 pandas DataFrame
• 使用 GSLIB.locmap() 和 GSLIB.hist() 探索数据
• 检查聚集情况；如有需要，使用 geostats.declus() 进行去聚集
• 使用 geostats.nscore() 应用正态分数变换
• 使用 geostats.gamv() 计算实验变异函数（先计算各向同性）
• 检查各方向的变异函数以判断各向异性（方位角 0、45、90、135）
• 使用 GSLIB.make_variogram() 拟合变异函数模型
• 将模型叠加到实验变异函数上，验证拟合效果
• 使用 geostats.kb2d() 执行克里金（普通克里金使用 ktype=1）
• 运行 SGSIM 量化不确定性（50-100 次实现）
• 使用 geostats.backtr() 对结果进行反变换
• 使用交叉验证或留出数据进行验证

常见问题

问题	解决方案
变异函数未达到基台值	增大 nlag 或 xlag，以覆盖完整变程
克里金产生负值	在克里金之后而非之前进行反变换
SGSIM 出现伪影	检查网格定义（xmn、xsiz）是否与数据范围一致
变异函数对数过少	增大 atol（角度容差）或 xltol（滞后容差）
变异函数出现孔穴效应	可能表示存在周期性；尝试使用嵌套结构

提示

1. 始终转换为正态分数——大多数方法都假设数据服从高斯分布
2. 从各向同性开始——只有在方向变差函数有依据时才添加各向异性
3. 检查变差函数配对数——确保每个滞后距都有足够的配对数
4. 进行多次实现——使用 50–100 次以上进行不确定性量化
5. 最后进行反变换——仅应用于最终结果

参考资料

• GSLIB Programs——GSLIB 程序参考
• Simulation Methods——SGSIM、SISIM 和其他方法

脚本

• scripts/kriging_example.py——完整的克里金工作流