Inversion Module

Core Inversion Classes

Single-time ERT inversion functionality.

PyHydroGeophysX.inversion.ert_inversion.ERROR_SOURCES = ('file', 'estimate', 'max', 'stack', 'reciprocal')

How err is chosen. file trusts the err column the instrument wrote and only estimates where it is missing; estimate always recomputes from relative_error/absolute_error; max takes the larger of the two per datum, which is the conservative reading when the file’s errors look optimistic. stack adds each reading’s stacking spread - the scatter of the potential’s samples within the reading, as a fraction of it, which Subsurface Insights records and the stack token carries - to the estimate in quadrature. That spread is not a reading’s uncertainty and usually overstates it, so it is only ever asked for. reciprocal gives each reading the error of a reciprocal error model dR = 10**b * |R|**m fitted to the survey’s reciprocal pairs (error_model: m, b, floor; see ert_io.reciprocal_model_errors). file is the default because silently discarding measured errors makes chi2 report on an error model the data never had.

class PyHydroGeophysX.inversion.ert_inversion.ERTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs)

Bases: InversionBase

Single-time ERT inversion class.

run(initial_model: ndarray | None = None, reference_model: ndarray | None = None) → InversionResult

Run ERT inversion.

Parameters:

• initial_model – Initial model parameters (if None, a homogeneous model is used)

• reference_model – Model the regularization pulls toward. Defaults to the initial model, which is the usual smoothness-from-homogeneous setup. Pass it explicitly when continuing an interrupted run. Leaving it to default makes the restart point the new reference, so the penalty becomes roughness of the change since the restart rather than roughness of the model; accumulated structure then goes unpenalized and the run drifts steadily under-regularized.

Returns:

InversionResult with inversion results

setup()

Set up ERT inversion (create operators, matrices, etc.)

zone_prior()

The a-priori zones on this inversion’s parameter mesh, or None.

Worked out once the forward operator has built the parameter mesh, whose cell order is the model vector’s. See ert_zones.

class PyHydroGeophysX.inversion.ert_inversion.ERTRun(lam: float, chi2: float, iterations: int, stop: str, convergence: ~typing.List[float], model: ~numpy.ndarray, response: ~numpy.ndarray, mesh: ~typing.Any, coverage: ~numpy.ndarray | None = None, manager: ~typing.Any = None, metrics: ~typing.Dict[str, ~typing.Any] = <factory>, volume: ~typing.Any = None)

Bases: object

One inversion at one lambda, described the same way by either engine.

chi2: float

convergence: List[float]

coverage: ndarray | None = None

iterations: int

lam: float

manager: Any = None

mesh: Any

metrics: Dict[str, Any]

model: ndarray

response: ndarray

property result

A manager-like handle for the viewer and the VTK export.

stop: str

volume: Any = None

(mesh, model) of the whole volume, when a profile was inverted in 3-D and mesh/model are the section along the line.

PyHydroGeophysX.inversion.ert_inversion.GEOMETRIC_FACTOR_POLICIES = ('off', 'check', 'fix')

What to do about geometric factors before anything else touches the data. fix validates and, if the check fails, recomputes k numerically on the inversion mesh and rebuilds rhoa from the measured transfer resistance.

PyHydroGeophysX.inversion.ert_inversion.ensure_geometric_factors(container, mesh, *, policy: str = 'fix', tolerance: float = 0.05, rho0: float = 100.0, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Validate k before the data are used, and repair it when asked.

A repair is kept only if it actually makes the check pass; otherwise the original factors are restored, because a failed repair is worse than a reported problem. Either way the caller is told what happened.

PyHydroGeophysX.inversion.ert_inversion.profile_section(container, volume, para) → Tuple[Any, ndarray]

(section, cells): the section along a profile laid into volume (_section_mesh()) and, for each of its cells, the cell of para (the volume’s parameter mesh) whose model it shows.

PyHydroGeophysX.inversion.ert_inversion.repair_geometric_factors(container, mesh, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Recompute k numerically on the inversion mesh and rebuild rhoa from R.

The transfer resistance R = rhoa/k is what the instrument actually measured; k and rhoa are both derived from it. So the repair keeps R fixed, takes k from a forward run on the mesh the inversion will use, and rebuilds rhoa = R * k. Numerical factors also carry the topography that the analytic half-space formula cannot.

R has to be recovered with the factors that formed rhoa, which is the file’s own k_file when the loader kept it, and only otherwise the container’s k. Dividing by the wrong one would bake the discrepancy into R and leave the section scaled after the repair rather than before it.

Modifies container in place and returns what changed.

PyHydroGeophysX.inversion.ert_inversion.run_ert_manager_inversion(data_path: str | ~pathlib.Path, output_dir: str | ~pathlib.Path, *, relative_error: float = 0.03, absolute_error: float = 0.0, error_source: str = 'file', error_floor: float = 0.005, error_model: ~typing.Dict[str, ~typing.Any] | None = None, mesh_quality: float = 34.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, para_boundary: float = 2.0, surface_nodes: int = 1, outer_width: float = 0.0, outer_max_cell_size: float = 0.0, mesh_file: str = '', zones: ~typing.Sequence[~typing.Dict[str, ~typing.Any]] | None = None, conform_to_zones: bool = False, decouple_zones: bool = False, lam: float = 50.0, max_iterations: int = 20, plateau_tolerance: float = 0.005, max_total_iterations: int = 60, engine: str = 'pyhydro', model_constraints: ~typing.Tuple[float, float] = (0.01, 100000.0), solver: str = 'cgls', geometric_factor_policy: str = 'fix', geometric_factor_tolerance: float = 0.05, instrument: str | None = None, reject_outliers: bool = False, outlier_threshold: float = 3.0, outlier_passes: int = 2, min_data_fraction: float = 0.5, auto_lambda: bool = False, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_lambda_trials: int = 6, lambda_bounds: ~typing.Tuple[float, float] = (0.001, 100000.0), lambda_warm_start: bool = True, lambda_cold_retry_chi2: float = 15.0, e4d: ~typing.Dict[str, ~typing.Any] | None = None, r2: ~typing.Dict[str, ~typing.Any] | None = None, survey_3d: bool = False, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert one ERT dataset, in the order that actually lowers chi2.

The stages run back to back, each optional:

0. Geometric factors (geometric_factor_policy). A homogeneous forward run must return the model resistivity; if it does not, k disagrees with the geometry being modelled and the whole section is scaled by that factor with no trace in chi2. fix recomputes k numerically on the inversion mesh and rebuilds rhoa from the measured transfer resistance.

1. Error model. error_source decides whether the file’s own err column is trusted, recomputed from relative_error/absolute_error, combined with that estimate, built from each reading’s stacking spread, or given by a reciprocal error model (error_model; see ERROR_SOURCES). Overwriting a measured error with an assumed one makes chi2 report on an error model the data never had.

2. Inversion at the requested lambda, iterated to a plateau. A run that exhausts max_iterations is continued (up to max_total_iterations) rather than being judged where it stopped. This run is always kept, under fixed_lambda.

3. Outlier rejection (reject_outliers). Measurements the converged model cannot explain are dropped and the inversion repeated, at most outlier_passes times and never below min_data_fraction of the data.

4. Lambda search (auto_lambda). Only now, with a plateaued misfit on the cleaned data, is lambda allowed to move; every trial is itself iterated to a plateau before its chi2 counts. With lambda_warm_start each trial continues from the nearest lambda already solved rather than restarting from a homogeneous model, which is why lam should start on the smooth side: the sweep is then a relaxation from an over-regularized model down to the data, the direction in which continuation is stable.

engine selects the solver: "pyhydro" for the in-house Gauss-Newton inversion (ERTInversion), "pygimli" for ert.ERTManager, or "adtlert" for the optional differentiable ADTLERT backend when a CUDA GPU, cuDSS and NVRTC are available: 2.5D on a profile mesh, 3D on a tetrahedral mesh (mesh_file). Without them, "adtlert" falls back to "pyhydro".

ADTLERT matches the published real-data branch: fast normal sensitivity, no Robin-boundary derivative, line search, a maximum log step of one, and GPU CGLS when CUDA is available.

"e4d" runs PNNL’s E4D, an external 3-D program, with e4d holding its settings (E4DSettings: the launcher, the program, the MPI processes). E4D runs on Linux, on Windows inside WSL 2, and on macOS when built from source; elsewhere its run folder is written under output_dir/e4d and the run stops with E4DNotRun saying how to run it. It never falls back to another engine. A profile is inverted on a 3-D mesh built around the line and read back as the section along it.

"r2" and "r3t" run Andrew Binley’s R2 (2-D profiles) and R3t (3-D surveys, and profiles on a 3-D mesh around the line), with r2 holding their settings (R2Settings: the program, the launcher). They are Windows programs, run natively there and through Wine on Linux and macOS, found in an installed ResIPy unless named; elsewhere their run folder is written under output_dir/r2 (or r3t) and the run stops with R2NotRun. Both search their smoothing weight at every iteration themselves, so lam and auto_lambda do not apply to them; max_iterations and target_chi2 (as an RMS of sqrt(target_chi2)) do, and they hold fixed zones fixed.

zones carries a-priori resistivity zones (see ert_zones): polygons whose value the inversion starts from and regularizes toward, and which the in-house engine can hold fixed. The other engines use the values as the starting model only, and the log says so. With conform_to_zones a generated mesh follows the zone outlines, so no cell straddles one; with decouple_zones the smoothness constraint is dropped across them, so the model may jump at a zone’s edge (every engine honours it).

The mesh options (mesh_quality to outer_max_cell_size) are those of build_inversion_mesh().

survey_3d: the electrodes cover an area and sit at (x, y, z), so the survey is inverted in 3-D - by the in-house and PyGIMLi engines on a generated prism mesh, by any engine on an imported 3-D mesh - and the model is the volume, exported as resistivity_3d.vtk too.

PyHydroGeophysX.inversion.ert_inversion.validate_geometric_factors(container, mesh, *, rho0: float = 100.0, tolerance: float = 0.05, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Check k by forward-modelling a homogeneous half space.

A homogeneous model of resistivity rho0 must return rho0 as the apparent resistivity of every configuration. The forward response is (U/I) * k, so the returned ratio is exactly k / k_true for the geometry actually being modelled, and it is independent of the field data.

This catches what chi2 cannot. A geometric factor that is uniformly wrong by a factor X forces the inversion to scale the model by 1/X to fit the same apparent resistivities, so the section is wrong by X while the fit looks perfect. Topography, wrong electrode spacing, and half-space versus full-space convention errors all show up here.

Returns the ratio statistics, an ok flag, and a message naming the likely cause and the consequence.

The mesh an ERT inversion runs on: built, imported, cut at zones, previewed.

One module for every place that needs the inversion mesh - the single-survey and the time-lapse inversions, and the ERT page’s mesh preview - so what the page shows before a run is the mesh the run then inverts on, cell for cell:

• build_inversion_mesh() builds PyGIMLi’s parameter mesh of a survey, along the outlines of a-priori zones (ert_zones) when asked, a prism volume under a 3-D survey’s electrodes (survey_volume_mesh()), or imports one built elsewhere;

• load_inversion_mesh() reads such a mesh - PyGIMLi, Gmsh, VTK, E4D - and checks that it can hold the survey;

• mark_zone_interfaces() stops the smoothness at the zone outlines;

• mesh_preview() splits the mesh the way the inversion does, for the page to draw.

These were written in ert_inversion, which still exports them.

PyHydroGeophysX.inversion.ert_mesh.MAX_MESH_QUALITY = 34.0

Triangle stops terminating a little above a 34-degree smallest angle (36 already never finishes on a plain profile), so a higher mesh quality is refused up front instead of hanging the run.

PyHydroGeophysX.inversion.ert_mesh.MESH_SUFFIXES = ('.bms', '.msh', '.vtk', '.vtu', '.poly', '.node', '.ele', '.cfg')

Mesh formats a user can hand the inversion. .bms is PyGIMLi’s own, .msh is Gmsh (the usual route for a complex 3D domain), .node and .ele an E4D (TetGen) mesh, .cfg an E4D mesh configuration, which is meshed on import; the rest are what PyGIMLi’s loader recognises.

PyHydroGeophysX.inversion.ert_mesh.ZONE_INTERFACE_MARKER = 12

one between cells of two different zones, or of a zone and the ground around it.

Type:

Marker of an edge the smoothness constraint does not cross

PyHydroGeophysX.inversion.ert_mesh.build_inversion_mesh(data, *, mesh_quality: float = 34.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, para_boundary: float = 2.0, surface_nodes: int = 1, outer_width: float = 0.0, outer_max_cell_size: float = 0.0, conform_zones: ~typing.Sequence[~typing.Dict[str, ~typing.Any]] | None = None, mesh_file: str = '', log: ~typing.Callable[[str], None] = <function noop>, report: ~typing.Dict[str, ~typing.Any] | None = None, survey_3d: bool = False)

The mesh an ERT inversion of data runs on: imported, or built.

One function for every place that needs it - the single-survey and the time-lapse inversions, and the ERT page’s mesh preview - so what the page shows before a run is the mesh the run then inverts on, cell for cell.

A generated mesh is PyGIMLi’s parameter mesh: an inverted region under the electrodes and a coarse outer region around it that carries the boundary condition far away. Every sizing option maps onto one of PyGIMLi’s, and a value left at its default is not passed, so the mesh is the one PyGIMLi builds by itself. They are the controls an E4D mesh configuration sets too: the inverted region is its fine zone, the outer region its outer boundary.

Parameters:

• data (pygimli.DataContainerERT) – The survey; its sensor positions define the generated mesh.

• mesh_quality (float) – Smallest angle, in degrees, a triangle of a generated mesh may have: higher gives better-shaped and more cells. At most MAX_MESH_QUALITY.

• para_depth (float) – Depth of the inverted region and its largest cell area (m²); 0 lets PyGIMLi choose from the array.

• para_max_cell_size (float) – Depth of the inverted region and its largest cell area (m²); 0 lets PyGIMLi choose from the array.

• para_boundary (float) – How far the inverted region reaches past the first and the last electrode, in electrode spacings (paraBoundary).

• surface_nodes (int) – Nodes placed on the surface between two neighbouring electrodes (addNodes); more makes the cells near the surface smaller.

• outer_width (float) – How far the outer region reaches beyond the inverted region, sideways and down, in lengths of the electrode spread (boundary); 0 keeps PyGIMLi’s four.

• outer_max_cell_size (float) – Largest cell area in the outer region (m², boundaryMaxCellSize); 0 leaves it unlimited.

• conform_zones (list of dict, optional) – A-priori zones (see ert_zones) whose outlines a generated mesh follows: they become cell edges, so no cell straddles a zone’s edge.

• mesh_file (str) – A mesh built elsewhere (.bms, .msh, .vtk, .vtu, .poly, E4D). Its own domain applies, and the options above are ignored.

• survey_3d (bool) – The electrodes cover an area (ert_io.survey_is_3d) and sit at (x, y, z): the generated mesh is survey_volume_mesh(), and an imported one must be 3-D.

• report (dict, optional) – Filled with what was built: source (“generated” or “imported”), zone_outline_edges, the edges the zone outlines added (0 when the mesh does not follow them), and for an imported mesh profile_in_volume: a profile laid into a 3-D mesh, whose electrodes in data now carry the volume’s coordinates.

Returns:

The full forward mesh: parameter domain (marker 2 and above) and the outer region that carries the boundary far from the array (marker 1).

Return type:

pygimli.Mesh

PyHydroGeophysX.inversion.ert_mesh.load_inversion_mesh(mesh_path: str | ~pathlib.Path, data=None, log: ~typing.Callable[[str], None] = <function noop>, *, place_profile: bool = False, survey_3d: bool = False)

Load a user-supplied inversion mesh and check it can hold this survey.

Building a mesh from the electrode line is fine for a 2D profile and hopeless for a 3D domain with topography, boreholes or known structure, so those are meshed externally (usually in Gmsh or E4D) and brought in here. An E4D mesh is read with its .trn translation undone and its zones as cell markers, smoothed across its internal boundaries within a zone as E4D does; an E4D .cfg is meshed first (see core.e4d_mesh).

An imported mesh fails in ways a generated one cannot: electrodes outside the domain, or every cell marked background so nothing is inverted. Both surface deep inside the forward solver as errors that name nothing useful, so they are checked here where the message can say what is wrong.

A 2-D profile fits a 3-D mesh when its electrodes, laid along x at y = 0 (_volume_positions()), lie inside it. place_profile then moves the electrodes of data there, for the run that inverts on the mesh; otherwise data is left as it is. A 3-D survey (survey_3d) needs a 3-D mesh, and its electrodes are where they are.

PyHydroGeophysX.inversion.ert_mesh.mark_zone_interfaces(mesh, zones, *, marker: int = 12)

A copy of mesh on which the smoothness stops at the zone outlines.

Every edge between two cells of the parameter domain that belong to different zones - or to a zone and the ground around it - is given marker. PyGIMLi’s region manager drops the smoothness constraint across a marked edge, so the in-house engines and PyGIMLi’s own manager let the model jump there; the ADTLERT engine reads the same marks (_interface_regions()). A cell belongs to the zone its centre lies in, as it does everywhere else, so on a mesh that follows the outlines the cut runs along them exactly, and on any other mesh along the cell edges closest to them.

Returns (mesh, edges): the marked copy and how many edges were cut; the mesh itself, and 0, for no zones or a 3-D mesh, on which zones do not apply.

PyHydroGeophysX.inversion.ert_mesh.mesh_preview(mesh, data=None) → Dict[str, Any]

The inversion mesh as the ERT page’s mesh tab draws it.

The mesh is split the way the inversion splits it - through the forward operator’s region manager - into the parameter domain, the cells that are inverted, in the order of the model vector, and the outer region, which is never inverted and only carries the boundary condition far from the array: pyGIMLi’s background region, the “infinite” part of the mesh.

Parameters:

• mesh (pygimli.Mesh) – From build_inversion_mesh().

• data (pygimli.DataContainerERT, optional) – The survey, for the electrode positions.

Returns:

dim, the cell counts and the electrode positions always (x and the elevation on a 2-D mesh, x, y, z in a volume). For a 2-D mesh also the cell polygons of both parts as (cells, corners, 2) arrays, para_centers (what zone_prior() takes) and both extents as (xmin, xmax, zmin, zmax).

Return type:

dict

PyHydroGeophysX.inversion.ert_mesh.profile_in_volume(data, mesh) → bool

Whether data is a line of electrodes along x at y = 0 in a 3-D mesh: a profile inverted in 3-D, read back as the section along it.

PyHydroGeophysX.inversion.ert_mesh.survey_volume_mesh(data, *, para_depth: float = 0.0, para_max_cell_size: float = 0.0, para_boundary: float = 2.0, outer_width: float = 0.0, quality: float = 34.0, log: ~typing.Callable[[str], None] = <function noop>)

A prism mesh under a 3-D survey’s electrodes (x, y, z).

The settings mean what they mean for a profile’s mesh. The inverted region (marker 2) reaches para_boundary electrode spacings past the outermost electrodes and para_depth below the ground (0: 0.4 of the array’s larger extent, as PyGIMLi sizes a profile’s); its plan triangles are at most para_max_cell_size m² (0: a spacing squared) and every electrode is one of their corners. Its layers start a third of a spacing thick and grow by 20 % each. The defaults keep a 50-electrode grid near 5000 inverted cells: a 3-D forward model refines the mesh eightfold. The outer region (marker 1) reaches outer_width spreads further sideways and down (0: four, as on a profile’s mesh), in coarse cells. The ground passes through the electrodes, flat beyond them.

PyHydroGeophysX.inversion.ert_mesh.tetrahedral_mesh(mesh)

mesh with every prism split into three tetrahedra, markers kept.

Each quadrilateral face is cut along the diagonal from its lowest-numbered node, so two prisms sharing a face cut it alike and the tetrahedra fit together (Dompierre et al., 1999, Int. J. Numer. Meth. Eng.). A mesh of tetrahedra comes back as it is; ADTLERT’s 3-D forward takes nothing else.

Seismic Refraction Tomography (SRT) inversion functionality.

Uses PyGIMLi’s TravelTimeManager for forward modeling and Jacobian computation. Provides custom Gauss-Newton inversion with the same architecture as ERTInversion but for travel-time data.

class PyHydroGeophysX.inversion.srt_inversion.SRTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Seismic Refraction Tomography inversion class.

Inverts travel-time data for subsurface velocity via log-slowness, and uses a Gauss-Newton optimization loop.

run(initial_model: ndarray | None = None) → InversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

PyHydroGeophysX.inversion.srt_inversion.build_srt_mesh(data, *, mesh_quality: float = 32.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, log: ~typing.Callable[[str], None] = <function noop>)

Build the travel-time inversion mesh the way the ERT pipeline does.

PyGIMLi sizes the parameter domain from the array length when paraDepth is left at 0. For refraction that reaches well past where any ray turns, so the deep cells are unconstrained and only slow the run down; capping the depth removes them. Cell size and quality trade resolution against cost the same way they do for ERT.

PyHydroGeophysX.inversion.srt_inversion.run_srt_manager_inversion(travel_time_path: str | ~pathlib.Path, output_dir: str | ~pathlib.Path, *, engine: str = 'pygimli', lam: float = 50.0, max_iterations: int = 20, plateau_tolerance: float = 0.005, max_total_iterations: int = 60, mesh_quality: float = 32.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, secondary_nodes: int = 3, auto_lambda: bool = False, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_lambda_trials: int = 6, lambda_warm_start: bool = True, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert one travel-time dataset for velocity.

With auto_lambda the same machinery as the ERT pipeline applies: each lambda is iterated to a plateau before its chi2 counts, the sweep relaxes from the requested lambda downward, and each trial continues from the nearest lambda already solved. engine="pyhydro" uses the in-house Gauss-Newton solver, which is what the search can drive; "pygimli" runs TravelTimeManager once and is the historical default.

The mesh is built once here and handed to whichever engine runs, so the two invert the same domain. para_depth and para_max_cell_size take 0 to mean “let PyGIMLi size it from the array”; secondary_nodes refines the ray tracing without adding unknowns.

Time-lapse ERT inversion functionality.

class PyHydroGeophysX.inversion.time_lapse.TimeLapseERTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs)

Bases: InversionBase

Time-lapse ERT inversion class.

run(initial_model: ndarray | None = None) → TimeLapseInversionResult

Run time-lapse ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

TimeLapseInversionResult with inversion results

setup()

Set up time-lapse ERT inversion (load data, create operators, matrices, etc.)

zone_prior()

The a-priori zones on this inversion’s parameter mesh, or None.

Every survey shares the mesh, so one prior serves them all. See ert_zones.

Time-lapse Seismic Refraction Tomography (SRT) inversion functionality.

Jointly inverts multiple travel-time datasets with spatial and temporal regularization using log-slowness parameterization.

class PyHydroGeophysX.inversion.srt_time_lapse.TimeLapseSRTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Time-lapse Seismic Refraction Tomography inversion.

Follows the structure of TimeLapseERTInversion, using TravelTimeManager and log-slowness model parameters.

run(initial_model: ndarray | None = None) → TimeLapseInversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

Windowed time-lapse ERT inversion for handling large temporal datasets.

class PyHydroGeophysX.inversion.windowed.WindowedTimeLapseERTInversion(data_dir: str, ert_files: List[str], measurement_times: List[float], window_size: int = 3, mesh: pygimli.Mesh | str | None = None, engine: str = 'pyhydro', log: Callable[[str], None] | None = None, **kwargs)

Bases: object

Class for windowed time-lapse ERT inversion to handle large temporal datasets.

run(window_parallel: bool = False, max_window_workers: int | None = None) → TimeLapseInversionResult

Run windowed time-lapse ERT inversion.

Parameters:

• window_parallel – Whether to process windows in parallel

• max_window_workers – Maximum number of parallel workers (None for auto)

Returns:

TimeLapseInversionResult with stitched results

EM Inversion

Inversion utilities for Time-Domain Electromagnetic (TDEM) data.

This module provides classes for 1D TDEM inversion using SimPEG, with support for both L2 and sparse (IRLS) regularization.

class PyHydroGeophysX.inversion.tdem_inversion.TDEMInversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, source_location: ndarray | None = None, receiver_location: ndarray | None = None, n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)

Bases: object

Class for 1D TDEM inversion using SimPEG.

This class provides functionality for inverting TDEM sounding data to recover 1D layered Earth conductivity models.

Example (dobs and uncertainties are your sounding’s data, one value per time gate):

>>> times = np.logspace(-5, -2, 31)
>>> inv = TDEMInversion(
...     times=times,
...     dobs=dobs,
...     uncertainties=uncertainties,
...     source_radius=10.0
... )
>>> result = inv.run()

plot_result(result: TDEMInversionResult, true_model: Tuple[ndarray, ndarray] | None = None, save_path: str | None = None) → None

Plot inversion results.

Parameters:

• result – TDEMInversionResult from run()

• true_model – Tuple of (thicknesses, conductivity) for true model

• save_path – Path to save figure (optional)

run(starting_model: ndarray | None = None) → TDEMInversionResult

Run TDEM inversion.

Parameters:

starting_model – Initial log-conductivity model (optional)

Returns:

TDEMInversionResult with inversion results

setup() → None

Set up inversion components (survey, mesh, simulation).

class PyHydroGeophysX.inversion.tdem_inversion.TDEMInversionResult(recovered_model: ndarray = None, recovered_conductivity: ndarray = None, l2_model: ndarray = None, l2_conductivity: ndarray = None, predicted_data: ndarray = None, mesh: discretize.TensorMesh = None, thicknesses: ndarray = None, chi2: float = None, iterations: int = 0, convergence_history: List[float] = <factory>)

Bases: object

Container for TDEM inversion results.

recovered_model

Final recovered model (log-conductivity)

Type:

numpy.ndarray

recovered_conductivity

Recovered conductivity (S/m)

Type:

numpy.ndarray

l2_model

L2 model before IRLS (log-conductivity)

Type:

numpy.ndarray

l2_conductivity

L2 conductivity (S/m)

Type:

numpy.ndarray

predicted_data

Predicted data from recovered model

Type:

numpy.ndarray

mesh

TensorMesh used for inversion

Type:

discretize.TensorMesh

thicknesses

Layer thicknesses used in inversion

Type:

numpy.ndarray

chi2

Final chi-squared misfit

Type:

float

iterations

Number of iterations

Type:

int

convergence_history

History of data misfit per iteration

Type:

List[float]

PyHydroGeophysX.inversion.tdem_inversion.run_tdem_inversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, n_layers: int = 25, use_irls: bool = True, verbose: bool = True, **kwargs) → TDEMInversionResult

Convenience function to run TDEM inversion.

Parameters:

• times – Time channels (s)

• dobs – Observed data (T)

• uncertainties – Data uncertainties (T)

• source_radius – Loop radius (m)

• n_layers – Number of inversion layers

• use_irls – Use IRLS for sparse inversion

• verbose – Print progress

• **kwargs – Additional parameters for TDEMInversion

Returns:

TDEMInversionResult with inversion results

Inversion utilities for Frequency-Domain Electromagnetic (FDEM) data.

Provides 1D layered-earth FDEM inversion using SimPEG.

class PyHydroGeophysX.inversion.fdem_inversion.FDEMInversion(frequencies: ndarray, dobs: ndarray, uncertainties: ndarray, thicknesses: ndarray | None = None, source_location: ndarray | None = None, source_radius: float = 10.0, receiver_location: ndarray | None = None, receiver_orientation: str = 'z', receiver_component: str = 'secondary', waveform_type: str = 'dipole', n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)

Bases: object

1D FDEM inversion using SimPEG.

Follows the same pattern used in TDEMInversion. Source, receiver and waveform_type are checked, and an omitted receiver placed, as in FDEMSurveyConfig: the origin, unless that is directly below or above a dipole source (then 10 m along x from it) or away from a loop’s centre (then the centre on z = 0). Explicit zero horizontal dipole offsets are rejected rather than producing a non-finite prediction. Receivers at (n, 3) locations take complex dobs ordered as FDEMForwardModeling returns them: the n receivers of each frequency together.

run(starting_model: ndarray | None = None) → FDEMInversionResult

setup() → None

class PyHydroGeophysX.inversion.fdem_inversion.FDEMInversionResult(recovered_model: ndarray = None, recovered_conductivity: ndarray = None, l2_model: ndarray = None, l2_conductivity: ndarray = None, predicted_data: ndarray = None, mesh: Any = None, thicknesses: ndarray = None, chi2: float = None, frequencies: ndarray = None)

Bases: object

Container for FDEM inversion results.

chi2: float = None

frequencies: ndarray = None

l2_conductivity: ndarray = None

l2_model: ndarray = None

mesh: Any = None

predicted_data: ndarray = None

recovered_conductivity: ndarray = None

recovered_model: ndarray = None

thicknesses: ndarray = None

Joint and Multi-Method Inversion

Joint ERT-SRT inversion with structural and geostatistical constraints.

This module integrates the legacy alternating inversion strategy from cg_joint_inversion_hang into the package architecture.

class PyHydroGeophysX.inversion.joint_ert_srt.JointERTSRTInversion(ert_data: str | PathLike | pygimli.DataContainer, srt_data: str | PathLike | pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Joint ERT-SRT inversion using alternating Gauss-Newton updates.

• ERT model parameter: log-resistivity

• SRT model parameter: log-slowness

• Coupling: linearized cross-gradient terms (B1, B2)

• Regularization: smoothness or geostatistical

run() → JointERTSRTResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.joint_ert_srt.JointERTSRTResult(ert_log_resistivity: ndarray | None = None, srt_log_slowness: ndarray | None = None, ert_resistivity: ndarray | None = None, srt_velocity: ndarray | None = None, ert_predicted: ndarray | None = None, srt_predicted: ndarray | None = None, ert_coverage: ndarray | None = None, srt_coverage: ndarray | None = None, chi2_ert: float | None = None, chi2_srt: float | None = None, iteration_history: list = <factory>, mesh: Any = None, meta: Dict[str, ~typing.Any]=<factory>)

Bases: object

Container for joint ERT-SRT inversion outputs.

chi2_ert: float | None = None

chi2_srt: float | None = None

ert_coverage: ndarray | None = None

ert_log_resistivity: ndarray | None = None

ert_predicted: ndarray | None = None

ert_resistivity: ndarray | None = None

iteration_history: list

mesh: Any = None

meta: Dict[str, Any]

srt_coverage: ndarray | None = None

srt_log_slowness: ndarray | None = None

srt_predicted: ndarray | None = None

srt_velocity: ndarray | None = None

Unified multi-method geophysical inversion interface.

class PyHydroGeophysX.inversion.multi_method.GeophysicalInversion(method: str, **kwargs)

Bases: object

Unified factory for multi-method geophysical inversion.

Dispatches to the correct inversion engine.

SUPPORTED = {'ert', 'fdem', 'joint', 'joint_ert_srt', 'srt', 'tdem'}

property engine

run(**kwargs)

setup(**kwargs)

Cross-Constraint Utilities

Cross-method constraint utilities for joint/cooperative inversion.

These utilities enable information sharing between geophysical methods and provide hydrology-to-geophysics coupling helpers.

class PyHydroGeophysX.inversion.cross_constraints.PetrophysicalCoupling

Bases: object

Coupling helpers from hydrological state to multi-method geophysics.

static compare_inversions_to_hydro(hydro_wc, ert_result, srt_result=None, em_result=None, petro_params=None)

Compare inversion products back to hydrological water content.

Returns misfit statistics for available methods. SRT conversion uses the same velocity_model selection as water_content_to_all_geophysics. Linear, Wyllie and Raymer have built-in inverses. For DEM/Hertz-Mindlin (including the default), supply petro_params['velocity_inverse'], a calibrated callable (velocity, porosity) -> water_content; these models are not silently replaced with a linear mixing rule.

An SRT model that cannot be converted does not cost the ERT and EM comparisons: its entry is {"error": reason} in place of water_content and stats, a UserWarning gives the reason, and summary['skipped'] names it. SRT cells whose velocity has no water content under the chosen model (outside its saturation range, or at the forward clipping limits) are NaN in water_content, left out of the statistics and counted in stats['n_masked'].

static water_content_to_all_geophysics(water_content, porosity, **params)

Convert water content to resistivity, velocity, and conductivity.

Uses existing petrophysical functions in the package.

class PyHydroGeophysX.inversion.cross_constraints.StructuralConstraint

Bases: object

Build structural constraint terms from one method for another.

static apply_structural_weights_to_Wm(Wm, boundary_weights)

Modify an existing smoothness matrix Wm to respect boundaries.

static build_cross_gradient_operator(mesh, model_a, model_b)

Build a cross-gradient operator surrogate.

Returns a diagonal sparse matrix weighted by local cross-gradient magnitude, where small values indicate better structural agreement.

static build_linearized_cross_gradient_blocks(RCM: Any, X: ndarray, model_a: ndarray, model_b: ndarray, mode: str = 'direct') → Tuple[scipy.sparse.csr_matrix, scipy.sparse.csr_matrix]

Build linearized cross-gradient blocks B1 and B2.

B1 multiplies model_a and B2 multiplies model_b. The resulting penalty terms are ||B1 m_a||^2 and ||B2 m_b||^2.

Row i is the weighted least-squares gradient estimate over the cells j with RCM[i, j] != 0, so B1 and B2 have the sparsity of RCM and are returned as CSR matrices. RCM may be dense or sparse. Each row is assembled from its own neighbourhood: forming X^T W X with W = diag(RCM[i]) from an n x n W for every row would be O(n^3) per call, on dense n x n blocks.

static build_local_design_matrix(mesh) → ndarray

Build the local design matrix X=[x, z, 1] used by cross-gradient.

The column of ones is the constant term of the local linear fit.

static build_neighborhood_matrix(mesh, Wm: Any | None = None, source: str = 'smoothness', correlation_lengths: Sequence[float] = (4.0, 4.0), threshold: float = 0.01, binarize: bool = True) → scipy.sparse.csr_matrix

Build the neighborhood/correlation matrix used by cross-gradient.

Parameters:

• mesh – Mesh object.

• Wm – Smoothness matrix. Required when source='smoothness'.

• source – One of 'smoothness' (legacy direct mode) or 'covariance' (legacy spatial mode).

• correlation_lengths – Correlation lengths passed to pg.utils.covarianceMatrix when source='covariance'.

• threshold – Entries with absolute value below this threshold are zeroed.

• binarize – If True, convert nonzero entries to 1.

Returns:

The n x n matrix, sparse. Wm^T Wm of a smoothness operator has a handful of entries per row, so it is kept sparse; a dense array would take 0.2 GB at 5,000 cells and grow as n^2.

Return type:

scipy.sparse.csr_matrix

static from_conductivity_model(conductivity, mesh, gradient_threshold: float = 0.3)

Extract structural boundaries from an EM conductivity model.

static from_velocity_model(velocity_model, mesh, gradient_threshold: float = 0.3)

Extract structural boundaries from an SRT velocity model.

Returns a pairwise-boundary weight map that can reduce smoothing across detected structural interfaces.

Base Classes

Base classes for geophysical inversion frameworks.

This module defines:

• InversionResult: A base class to store and manage common results from an inversion process, including final model, predicted data, convergence history, and plotting utilities.

• TimeLapseInversionResult: A specialized version of InversionResult for time-lapse studies, handling multiple models over time and providing time-slice plotting and animation.

• InversionBase: An abstract base class outlining the common structure and interface for various geophysical inversion methods (e.g., ERT, SRT). It handles data, mesh, and basic parameter management.

class PyHydroGeophysX.inversion.base.InversionBase(data: pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: object

Abstract base class for geophysical inversion methods.

This class provides a foundational structure for specific inversion techniques (e.g., ERT, SRT). It manages observed data, mesh, common inversion parameters, and an InversionResult object. Subclasses must implement methods for setting up the inversion, running the inversion loop, computing Jacobians, and evaluating the objective function.

compute_jacobian(model: ndarray) → ndarray

Abstract method to compute the Jacobian matrix (sensitivity matrix).

The Jacobian J_ij = ∂d_i / ∂m_j relates changes in model parameters (m_j) to changes in observed data (d_i).

Parameters:

model (np.ndarray) – The current model parameter vector for which to compute the Jacobian.

Returns:

The computed Jacobian matrix, typically of shape (n_data, n_model_params).

Return type:

np.ndarray

objective_function(model: ndarray, data_to_fit: ndarray | None = None) → float

Abstract method to compute the value of the objective function.

The objective function typically includes a data misfit term and one or more regularization terms: Φ(m) = Φ_d(m) + λ * Φ_m(m).

Parameters:

• model (np.ndarray) – The current model parameter vector.

• data_to_fit (Optional[np.ndarray], optional) – The observed data to fit. If None, self.data (from DataContainer) is typically used. Defaults to None.

Returns:

The calculated value of the objective function.

Return type:

float

run() → InversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.base.InversionResult

Bases: object

Base class to store, save, load, and plot results from a geophysical inversion.

final_model

The final inverted model parameters (e.g., resistivity, velocity). Typically a 1D array corresponding to mesh cells.

Type:

Optional[np.ndarray]

predicted_data

The data predicted by the forward model using the final_model.

Type:

Optional[np.ndarray]

coverage

Coverage or sensitivity values for the model parameters, often derived from the Jacobian or resolution matrix.

Type:

Optional[np.ndarray]

mesh

The PyGIMLi mesh object used in the inversion.

Type:

Optional[pg.Mesh]

iteration_models

A list storing the model parameters at each iteration of the inversion.

Type:

List[np.ndarray]

iteration_data_errors

A list storing the data misfit (e.g., residuals) at each iteration.

Type:

List[np.ndarray]

iteration_chi2

A list storing the chi-squared (χ²) value or a similar misfit metric at each iteration.

Type:

List[float]

meta

A dictionary to store any additional metadata about the inversion run (e.g., inversion parameters, timings, comments).

Type:

Dict[str, Any]

classmethod load(filename: str) → InversionResult

Load inversion results from a file previously saved by the save method.

Only load trusted pickle files. Missing or unreadable mesh sidecars produce warnings and leave mesh unset.

Parameters:

filename (str) – The base path to the saved results file. If saved as name.pkl and name.bms, provide name.pkl or name.

Returns:

An instance of InversionResult (or a subclass if called on one)

populated with the loaded data.

Return type:

InversionResult

Raises:

• FileNotFoundError – If the main data file is not found.

• IOError – If there’s an error during file reading.

• pickle.UnpicklingError – If the file cannot be unpickled.

plot_convergence(ax: Axes | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot the convergence curve (chi-squared misfit vs. iteration number).

Parameters:

• ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

• **kwargs (Any) – Additional keyword arguments passed directly to ax.plot (e.g., color, marker, linestyle).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.iteration_chi2 is empty (no convergence data).

plot_model(ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot the final inverted model on its associated mesh.

Parameters:

• ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

• cmap (str, optional) – The colormap to use for visualizing model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – If provided, cells with coverage values below this threshold will be masked (made transparent or semi-transparent) in the plot. Requires self.coverage to be populated. Defaults to None (no masking).

• **kwargs (Any) – Additional keyword arguments passed directly to pygimli.show (e.g., cMin, cMax, orientation, logScale).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.final_model or self.mesh is None.

save(filename: str) → None

Save the inversion results to a file using Python’s pickle format. If a mesh is present, it is saved separately as a PyGIMLi binary mesh file (.bms).

Parameters:

filename (str) – The base path (including filename without extension) to save the results. The main data will be saved as filename.pkl (or just filename if user includes .pkl). A trailing .pkl is removed before deriving the .bms mesh filename; other suffixes are retained.

Raises:

IOError – If writing or pickling the main data fails. Mesh-save failures are printed as warnings and do not prevent saving data.

class PyHydroGeophysX.inversion.base.TimeLapseInversionResult

Bases: InversionResult

Specialized class to store and manage results from time-lapse inversions. Inherits from InversionResult and adds attributes specific to time-lapse data.

final_models

A 2D NumPy array where each column represents the inverted model for a specific timestep (shape: num_cells x num_timesteps).

Type:

Optional[np.ndarray]

timesteps

An array or list of time values (e.g., hours, days) corresponding to each model slice in final_models.

Type:

Optional[np.ndarray]

all_coverage

A list where each element is the coverage array for the corresponding timestep’s model.

Type:

List[np.ndarray]

all_chi2

A list where each element is a list of chi-squared values per iteration for the inversion of that specific timestep or window.

Type:

List[Any]

create_time_lapse_animation(output_filename: str, cmap: str = 'viridis', coverage_threshold: float | None = None, dpi: int = 100, fps: int = 2, **kwargs: Any) → None

Create and save an animation (e.g., MP4 video) of the time-lapse inversion results.

Requires ffmpeg or another Matplotlib-supported animation writer to be installed.

Parameters:

• output_filename (str) – The filename for the output animation (e.g., ‘timelapse_animation.mp4’).

• cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

• dpi (int, optional) – Dots Per Inch for the output animation. Defaults to 100.

• fps (int, optional) – Frames Per Second for the animation. Defaults to 2.

• **kwargs (Any) – Additional keyword arguments passed to plot_time_slice for each frame.

Raises:

• ValueError – If final_models or mesh is missing.

• ImportError – If matplotlib.animation cannot be imported.

classmethod load(filename: str) → TimeLapseInversionResult

Load time-lapse inversion results. Overrides base class load.

plot_time_slice(timestep_idx: int, ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot a single time slice (inverted model at a specific timestep) from the results.

Parameters:

• timestep_idx (int) – The zero-based index of the timestep to plot.

• ax (Optional[plt.Axes], optional) – Matplotlib Axes to plot on. If None, creates new.

• cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

• **kwargs (Any) – Additional arguments passed to pg.show.

Returns:

The matplotlib Figure and Axes objects.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If models or mesh are missing, or if timestep_idx is out of range.

save(filename: str) → None

Save time-lapse inversion results. Overrides base class save.