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

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

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

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. Create mesh
hx = np.ones(80) * 5.0
hz = np.ones(40) * 2.5
mesh = TensorMesh([hx, hz], origin='CN')

# 2. Build survey (dipole-dipole)
n_electrodes = 24
electrode_spacing = 5.0
elec_x = np.arange(n_electrodes) * 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)
)
dobs = simulation.dpred(np.log(sigma_true))
dobs += 0.02 * np.abs(dobs) * np.random.randn(len(dobs))  # Add noise

# 4. Set up inversion
obs_data = data.Data(survey, dobs=dobs,
                     standard_deviation=0.05 * np.abs(dobs))
dmis = data_misfit.L2DataMisfit(data=obs_data, simulation=simulation)
reg = regularization.WeightedLeastSquares(
    mesh, alpha_s=1e-4, alpha_x=1, alpha_z=1
)
opt = optimization.InexactGaussNewton(maxIter=20)
inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)
dir_list = [
    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