gempy

GemPy - 3D 地质建模

快速参考

import gempy as gp

# Create model
geo_model = gp.create_geomodel(
    project_name='Model',
    extent=[0, 1000, 0, 1000, 0, 500],  # [xmin, xmax, ymin, ymax, zmin, zmax]
    resolution=[50, 50, 25]
)

# Add surface points and orientations
gp.add_surface_points(geo_model, x=[100, 500, 900], y=[500, 500, 500],
                      z=[400, 350, 400], surface='TopFormation')
gp.add_orientations(geo_model, x=[500], y=[500], z=[375],
                    pole_vector=[0, 0, 1], surface='TopFormation')

# Compute and visualize
gp.set_interpolator(geo_model)
sol = gp.compute_model(geo_model)
gp.plot_2d(geo_model, cell_number=[25], direction='y')

主要类

类	用途
GeoModel	主容器——存储所有模型数据
StructuralFrame	管理地质关系
Grid	用于插值的计算网格
Solutions	模型结果（岩性、标量场）

基本操作

定义地质关系

gp.map_stack_to_surfaces(geo_model, mapping={
    'Strata1': ['TopFormation', 'BaseFormation'],
    'Basement': ['Basement']
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.ERODE

添加断层

gp.add_surface_points(geo_model, x=[300, 300, 300], y=[200, 500, 800],
                      z=[100, 250, 400], surface='Fault1')
gp.add_orientations(geo_model, x=[300], y=[500], z=[250],
                    pole_vector=[1, 0, 0.5], surface='Fault1')
gp.map_stack_to_surfaces(geo_model, mapping={
    'Fault_Series': ['Fault1'],
    'Strata1': ['Layer1', 'Layer2'],
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.FAULT

从文件加载数据

import pandas as pd

points_df = pd.read_csv('surface_points.csv')  # Columns: X, Y, Z, surface
gp.add_surface_points(geo_model, x=points_df['X'], y=points_df['Y'],
                      z=points_df['Z'], surface=points_df['surface'])

ori_df = pd.read_csv('orientations.csv')  # Columns: X, Y, Z, dip, azimuth, surface
gp.add_orientations(geo_model, x=ori_df['X'], y=ori_df['Y'], z=ori_df['Z'],
                    dip=ori_df['dip'], azimuth=ori_df['azimuth'],
                    surface=ori_df['surface'])

可视化

gp.plot_2d(geo_model, cell_number=[25], direction='y', show_data=True)
gp.plot_3d(geo_model, show_data=True, show_surfaces=True)  # Requires PyVista

获取结果

sol = gp.compute_model(geo_model)
lithology = sol.raw_arrays.lith_block
lith_3d = lithology.reshape(geo_model.grid.regular_grid.resolution)

构造关系类型

类型	描述
ERODE	较年轻的地层侵蚀较老的地层（不整合）
ONLAP	较年轻的地层超覆于较老的地层之上
FAULT	地层面为断层面
INTRUSION	侵入体

何时使用及替代方案

场景	建议
基于地表数据进行 3D 隐式地质建模	GemPy——专为此设计，原生支持 Python
使用构造框架进行复杂褶皱建模	LoopStructural——褶皱支持更好
商业级地下建模	SKUA-GOCAD——行业标准，专有软件
只需快速创建 2D 剖面	GemPy 可以胜任，但可能功能过剩；可考虑手动插值
需要根据地质数据进行重力/磁力正演建模	GemPy——内置势场支持

选择 GemPy 的情况：需要根据地表接触点和产状进行隐式 3D 地质建模，并建立断层/不整合关系，尤其是需要与 Python 工作流集成时。对于标准场景，其学习曲线比 LoopStructural 更平缓。

避免使用 GemPy 的情况：地质构造以复杂褶皱为主（使用 LoopStructural），或者需要生产级储层建模（使用商业工具）。

常见工作流

根据地表数据构建 3D 地质模型

• 加载或定义地表接触点（每个地层至少 2 个）和产状（每个地层至少 1 个）
• 根据适当的范围和分辨率创建 GeoModel
• 使用 gp.add_surface_points() 添加地表点
• 使用 gp.add_orientations() 添加产状
• 使用 gp.map_stack_to_surfaces() 映射地质序列
• 设置构造关系（ERODE、ONLAP、FAULT）
• 使用 gp.set_interpolator() 设置插值器
• 使用 gp.compute_model() 计算模型
• 使用 gp.plot_2d() 通过 2D 剖面验证
• 使用 gp.plot_3d() 可视化 3D 结果，或导出为 VTK

常见问题

问题	解决方案
模型无法计算	检查每个地层是否至少有 2 个点和 1 个产状
边缘出现伪影	将模型范围扩展到数据范围之外
断层错移错误	检查 pole_vector 的方向
内存错误	降低网格分辨率

参考资料

• 构造关系 - ERODE、ONLAP、FAULT 的详细指南
• 数据要求 - 输入数据格式和约束
• 故障排除 - 常见问题及解决方案

脚本

• scripts/validate_input.py - 验证输入数据文件
• scripts/export_model.py - 将模型导出为 VTK、numpy 或 CSV