geophysical-inversion

Geophysical Inversion Workflow

End-to-end pipeline for inverting geophysical data, from survey data loading through mesh creation, forward modelling, inversion, gridding, and 3D visualization of recovered physical property models.

Skill Chain

simpeg / pygimli      verde            pyvista
[Mesh + Inversion] --> [Gridding]    --> [3D Visualization]
  |                     |                |
  Survey geometry       Interpolate      Volume render
  Forward model         Grid to raster   Slice views
  Misfit + reg          Trend removal    Overlay data
  Recover model         Cross-validate   Export mesh

Decision Points: SimPEG vs pyGIMLi

Criterion	SimPEG	pyGIMLi
DC resistivity / ERT	Yes	Yes (simpler API)
Magnetics	Yes	Limited
Gravity	Yes	Limited
Electromagnetics (TDEM, FDEM)	Yes	No
Seismic refraction (SRT)	No	Yes
Induced polarization	Yes	Yes
Built-in electrode arrays	Manual setup	Built-in (Wenner, Schlumberger, etc.)
Mesh types	TensorMesh, TreeMesh, CurvilinearMesh	Triangular, tetrahedral, structured
Joint inversion	Yes (Wires maps)	Limited
API complexity	More boilerplate, more flexible	Less boilerplate, opinionated

Rule of thumb: Use pyGIMLi for standard near-surface ERT/SRT surveys with conventional arrays. Use SimPEG for multi-physics, EM methods, potential fields, or research-grade custom inversions.

Step-by-Step Orchestration

Choose one inversion branch. Each branch explicitly produces cell_centers_xz_m, resistivity_ohm_m and model_vtk for the shared stages. X is distance along the profile in metres, Z is elevation in metres (negative below the surface); resistivity is in ohm metres. These short synthetic runs check API and data flow, not survey resolution or field-model reliability.

Stage 1a: Inversion with SimPEG (DC Resistivity Example)

# example: simpeg-inversion
import numpy as np
from discretize import TensorMesh
from simpeg.electromagnetics.static import resistivity as dc
from simpeg import maps, data, data_misfit, regularization
from simpeg import optimization, inverse_problem, inversion, directives

# 1. Small synthetic 2D mesh. SimPEG's second coordinate represents elevation.
hx = np.full(24, 5.0)
hz = np.full(12, 5.0)
mesh = TensorMesh([hx, hz], origin='CN')

# 2. Build survey (dipole-dipole)
n_electrodes = 8
electrode_spacing = 5.0
elec_x = (np.arange(n_electrodes) - (n_electrodes - 1) / 2) * electrode_spacing
elec_locs = np.c_[elec_x, np.zeros(n_electrodes)]

source_list = []
for i in range(n_electrodes - 3):
    rx = dc.receivers.Dipole(elec_locs[[i+2]], elec_locs[[i+3]])
    src = dc.sources.Dipole([rx], elec_locs[i], elec_locs[i+1])
    source_list.append(src)
survey = dc.Survey(source_list)

# 3. Forward model (for synthetic test)
sigma_true = np.ones(mesh.nC) * 0.01  # 100 ohm-m background
sigma_true[mesh.cell_centers[:, 1] < -20] = 0.1  # Conductive layer
simulation = dc.Simulation2DNodal(
    mesh, survey=survey, sigmaMap=maps.ExpMap(mesh)
)
clean_data = simulation.dpred(np.log(sigma_true))
# Voltage / unit-current uncertainties: fractional term plus nonzero floor.
standard_deviation = 0.02 * np.abs(clean_data) + 1e-6
rng = np.random.default_rng(42)
dobs = clean_data + standard_deviation * rng.standard_normal(survey.nD)

# 4. Set up inversion
obs_data = data.Data(survey, dobs=dobs, standard_deviation=standard_deviation)
dmis = data_misfit.L2DataMisfit(data=obs_data, simulation=simulation)
reg = regularization.WeightedLeastSquares(
    mesh, alpha_s=1e-4, alpha_x=1, alpha_y=1
)
opt = optimization.InexactGaussNewton(maxIter=3)
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)
dir_list = [
    directives.BetaEstimate_ByEig(beta0_ratio=1.0, random_seed=42),
    directives.BetaSchedule(coolingFactor=2),
    directives.TargetMisfit()
]
inv = inversion.BaseInversion(inv_prob, directiveList=dir_list)

# 5. Run inversion
m0 = np.log(np.ones(mesh.nC) * 0.01)  # Starting model
mrec = inv.run(m0)
sigma_rec = np.exp(mrec)  # Recovered conductivity