Forward Module

ERT and SRT Forward Modeling

Forward modeling utilities for Electrical Resistivity Tomography (ERT).

class PyHydroGeophysX.forward.ert_forward.ERTForwardModeling(mesh: pygimli.Mesh, data: pygimli.DataContainer | None = None)

Bases: object

Class for forward modeling of Electrical Resistivity Tomography (ERT) data.

classmethod create_synthetic_data(xpos: ndarray, ypos: ndarray | None = None, mesh: pygimli.Mesh | None = None, res_models: ndarray | None = None, schemeName: str = 'wa', noise_level: float = 0.05, absolute_error: float = 0.0, relative_error: float = 0.05, save_path: str | None = None, show_data: bool = False, seed: int | None = None, xbound: float = 100, ybound: float = 100) → Tuple[pygimli.DataContainer, pygimli.Mesh]

Create synthetic ERT data using forward modeling.

This method simulates an ERT survey by placing electrodes, creating a measurement scheme, performing forward modeling to generate synthetic data, and adding noise.

Parameters:

• xpos – X-coordinates of electrodes

• ypos – Y-coordinates of electrodes (if None, uses flat surface)

• mesh – Mesh for forward modeling

• res_models – Resistivity model values

• schemeName – Name of measurement scheme (‘wa’, ‘dd’, etc.)

• noise_level – Level of Gaussian noise to add

• absolute_error – Absolute error for data estimation

• relative_error – Relative error for data estimation

• save_path – Path to save synthetic data (if None, does not save)

• show_data – Whether to display data after creation

• seed – Random seed for noise generation

• xbound – X boundary extension for mesh

• ybound – Y boundary extension for mesh

Returns:

Tuple of (synthetic ERT data container, simulation mesh)

forward(resistivity_model: ndarray, log_transform: bool = True) → ndarray

Compute forward response for a given resistivity model.

Parameters:

• resistivity_model – Resistivity model values

• log_transform – Whether resistivity_model contains natural-log resistivity; when True the response is also returned in natural-log space.

Returns:

Apparent resistivity (ohm-m), or its natural logarithm when log_transform=True.

forward_and_jacobian(resistivity_model: ndarray, log_transform: bool = True) → Tuple[ndarray, ndarray]

Compute forward response and Jacobian matrix.

Parameters:

• resistivity_model – Resistivity model values

• log_transform – Whether resistivity_model is log-transformed

Returns:

Tuple of (response, Jacobian). With log_transform=True these are ln(apparent resistivity) and its derivative with respect to ln(model resistivity); otherwise both use linear resistivity.

get_coverage(resistivity_model: ndarray, log_transform: bool = True) → ndarray

Compute a sensitivity-based coverage proxy for each model cell.

Parameters:

• resistivity_model – Resistivity model values

• log_transform – Whether resistivity_model is log-transformed

Returns:

Coverage values for each cell

set_data(data: pygimli.DataContainer) → None

Set ERT data for forward modeling.

Parameters:

data – ERT data container

set_mesh(mesh: pygimli.Mesh) → None

Set mesh for forward modeling.

Parameters:

mesh – PyGIMLI mesh

PyHydroGeophysX.forward.ert_forward.check_resistivity_model(model: Any, *, log: bool = False, name: str = 'resistivity model', hint: str = 'Look for masked or inactive hydrologic cells, cells outside the interpolated data, or a failed conversion.') → ndarray

Refuse a model the ERT forward operator would model as something else.

pyGIMLi’s ERTModelling does not refuse a NaN cell: the NaN never reaches the stiffness matrix, and the apparent resistivities come back finite and plausible (fielddataline2 with one NaN cell, pgcore 1.5.5.post1). An infinite cell is accepted with only a printed warning, the rows it empties forced to Dirichlet. Zones refuse both and the petrophysical models clip to a finite range, so either one is a hole in the model - a masked or inactive hydrologic cell, a cell outside the interpolated data - and is refused before the forward call, as zero and negative resistivity are.

Parameters:

• model – Resistivity of each cell (ohm-m), or its natural log with log.

• log – The model holds ln(resistivity), so any finite value is valid.

• name – What the model is, for the error message.

• hint – Where a NaN or infinite cell may have come from, for the message.

Returns:

The model as a new flat float array.

Raises:

ValueError – naming how many cells are NaN, infinite or not positive, and the first of them.

PyHydroGeophysX.forward.ert_forward.ertforandjac(fob: Any, rhomodel: Any, xr: Any) → Any

Forward model and Jacobian for ERT.

Parameters:

• fob (pygimli.ERTModelling) – ERT forward operator.

• rhomodel (pg.RVector) – Resistivity model.

• xr (np.ndarray) – Log-transformed model parameter.

Returns:

Log-transformed forward response. J (np.ndarray): Jacobian matrix.

Return type:

dr (np.ndarray)

PyHydroGeophysX.forward.ert_forward.ertforandjac2(fob: Any, xr: Any, mesh: Any, *, response: Any = None, with_response: bool = False) → Any

Alternative ERT forward model and Jacobian using log-resistivity values.

Parameters:

• fob (pygimli.ERTModelling) – ERT forward operator.

• xr (np.ndarray) – Log-transformed model parameter.

• mesh (pg.Mesh) – Mesh for the forward model.

• response – The linear response of fob’s last forward solve, when that solve was for this very model and fob has solved nothing since - as ertforward2(..., with_response=True) just returned it. The solve is then not repeated. pyGIMLi’s createJacobian reads the potentials of the operator’s last response() call and does not check which model they belong to: given another model’s, it returns a Jacobian of neither. So pass it only on that guarantee.

• with_response – Also return the linear response.

Returns:

Log-transformed forward response. J (np.ndarray): Jacobian matrix. With with_response, the linear response as a third value.

Return type:

dr (np.ndarray)

PyHydroGeophysX.forward.ert_forward.ertforward(fob: Any, mesh: Any, rhomodel: Any, xr: Any) → Any

Forward model for ERT, with only the cells of marker 2 taken from xr.

Parameters:

• fob (pygimli.ERTModelling) – ERT forward operator.

• mesh (pg.Mesh) – Mesh whose cell markers pick the cells to update.

• rhomodel (pg.RVector) – Resistivity of every cell (ohm-m).

• xr (np.ndarray) – Natural-log resistivity of the cells of marker 2.

Returns:

Log-transformed forward response. rhomodel (pg.RVector): Updated resistivity model (ohm-m).

Return type:

dr (np.ndarray)

PyHydroGeophysX.forward.ert_forward.ertforward2(fob: Any, xr: Any, mesh: Any, *, with_response: bool = False) → Any

Simplified ERT forward model.

Parameters:

• fob (pygimli.ERTModelling) – ERT forward operator.

• xr (np.ndarray) – Log-transformed model parameter.

• mesh (pg.Mesh) – Mesh for the forward model.

• with_response – Also return the linear response the solve gave, which ertforandjac2(..., response=) can build the Jacobian from.

Returns:

Log-transformed forward response, and with with_response the linear response as a second value.

Return type:

dr (np.ndarray)

Forward modeling utilities for Seismic Refraction Tomography (SRT).

class PyHydroGeophysX.forward.srt_forward.SeismicForwardModeling(mesh: pygimli.Mesh, scheme: pygimli.DataContainer | None = None)

Bases: object

Class for forward modeling of Seismic Refraction Tomography (SRT) data.

classmethod create_synthetic_data(sensor_x: ndarray, surface_points: ndarray | None = None, mesh: pygimli.Mesh = None, velocity_model: ndarray | None = None, slowness: bool = False, shot_distance: float = 5, noise_level: float = 0.05, noise_abs: float = 1e-05, save_path: str | None = None, show_data: bool = False, verbose: bool = False, seed: int | None = None) → Tuple[pygimli.DataContainer, pygimli.Mesh]

Create synthetic seismic data using forward modeling.

This method simulates a seismic survey by placing geophones along a surface, creating a measurement scheme, and performing forward modeling to generate synthetic travel time data.

Parameters:

• sensor_x – X-coordinates of geophones

• surface_points – Surface coordinates for placing geophones [[x,y],…] If None, geophones will be placed on flat surface

• mesh – Mesh for forward modeling

• velocity_model – Velocity model values

• slowness – Whether velocity_model is slowness (1/v)

• shot_distance – Distance between shots

• noise_level – Level of relative noise to add

• noise_abs – Level of absolute noise to add

• save_path – Path to save synthetic data (if None, does not save)

• show_data – Whether to display data after creation

• verbose – Whether to show verbose output

• seed – Random seed for noise generation

Returns:

Tuple of (synthetic seismic data container, simulation mesh)

static draw_first_picks(ax, data, tt=None, plotva=False, **kwargs)

Plot first arrivals as lines.

Parameters:

• ax (matplotlib.axes) – axis to draw the lines in

• data (GIMLI::DataContainer) – data containing shots (“s”), geophones (“g”) and traveltimes (“t”)

• tt (array, optional) – traveltimes to use instead of data(“t”)

• plotva (bool, optional) – plot apparent velocity instead of traveltimes

Returns:

ax – the modified axis

Return type:

matplotlib.axes

forward(velocity_model: ndarray, slowness: bool = True) → ndarray

Compute forward response for a given velocity model.

Parameters:

• velocity_model – Velocity model values (or slowness if slowness=True)

• slowness – Whether velocity_model is slowness (1/v)

Returns:

Forward response (travel times)

set_mesh(mesh: pygimli.Mesh) → None

Set mesh for forward modeling.

Parameters:

mesh – PyGIMLI mesh

set_scheme(scheme: pygimli.DataContainer) → None

Set seismic data scheme for forward modeling.

Parameters:

scheme – Seismic data scheme

EM Forward Modeling

Forward modelling for time-domain electromagnetic (TDEM) soundings.

The layered-earth response comes from SimPEG. On top of it this module applies the instrument model that turns a modelled decay into the numbers a receiver records: convolution with the transmitter’s turn-off waveform, a cascade of first-order receiver low-pass stages, superposition of the bipolar transmitter train, and integration of each gate over its recorded window. Each of those is a standard step, and their parameters are read from the acquisition files rather than assumed; see PyHydroGeophysX.data_processing.em1d.

The whole instrument model is linear in the modelled response, so it is built once as a matrix and applied to the forward and the Jacobian alike. A dataset that does not describe an instrument keeps SimPEG’s direct receiver path.

PyHydroGeophysX.forward.tdem_forward.COMPILED_ENV = 'PHGX_TDEM_COMPILED'

Set to 0 to compute every TDEM response and Jacobian through SimPEG’s own dpred and getJ.

PyHydroGeophysX.forward.tdem_forward.GATE_SHAPE_NAMES = {1: 'tukey', 2: 'square'}

GateShape values, as the acquisition files record them.

Only the shapes this module implements appear. A file asking for one that does not is refused rather than quietly given the nearest thing, because a wrong window in the forward has nothing downstream to catch it. See PyHydroGeophysX.forward.em1d._gate_window_name().

PyHydroGeophysX.forward.tdem_forward.GATE_WINDOWS: Tuple[str, ...] = ('tukey', 'square', 'centre')

How a modelled decay becomes one gate value.

A recorded gate is an integral of the transient over a window, so the model of it is the same integral. The acquisition files name the window: GateShape selects its shape and GateShapePar1 gives the fraction of the width that is cosine-tapered, so the usual 0.667 leaves a flat top over the middle third.

tukey is that window. The integral runs in linear time, and the response between the samples of its own grid is interpolated with the local cubic Hermite rule described in _local_log_hermite_matrix(). Reading GateShapePar1 as the tapered fraction rather than as the flat top changes the result by 1.5 to 4 percent; the tapered-fraction reading is the one that reproduces responses stored in survey files.

Two things separate it from a flat average of the window, and they pull the same way. The taper draws weight in from the edges, and the window is symmetric in linear time while the gate centre a file records is the geometric mean of open and close, so it sits slightly late relative to it. Neither rule is trying to reproduce the centre value.

square is the same integral with no taper, for a survey whose files ask for the plain window. centre reads the response at the gate centre and is the fallback for data that records no gate windows at all.

class PyHydroGeophysX.forward.tdem_forward.TDEMForwardModeling(thicknesses: ndarray, survey_config: TDEMSurveyConfig | None = None, survey: simpeg.electromagnetics.time_domain.Survey | None = None)

Bases: object

Class for forward modeling of Time-Domain Electromagnetic (TDEM) data.

This class provides functionality for 1D layered Earth TDEM forward modeling using SimPEG’s time_domain module.

Example

>>> # Define layer model
>>> thicknesses = np.array([10.0, 30.0])
>>> conductivity = np.array([0.01, 0.1, 0.001])  # S/m
>>>
>>> # Create forward modeler
>>> fwd = TDEMForwardModeling(thicknesses=thicknesses)
>>>
>>> # Compute response
>>> response = fwd.forward(conductivity)

forward(conductivity: ndarray, log_input: bool = False) → ndarray

Compute forward response for a given conductivity model.

Parameters:

• conductivity – Conductivity values for each layer (S/m)

• log_input – If True, conductivity is log-transformed

Returns:

Forward response (magnetic flux density, T)

forward_with_noise(conductivity: ndarray, noise_level: float = 0.05, seed: int | None = None, log_input: bool = False) → Tuple[ndarray, ndarray, ndarray]

Compute forward response with added Gaussian noise.

Parameters:

• conductivity – Conductivity values for each layer (S/m)

• noise_level – Relative noise level (default 5%)

• seed – Random seed for reproducibility

• log_input – If True, conductivity is log-transformed

Returns:

Tuple of (noisy_data, clean_data, uncertainties)

get_times() → ndarray

Get the time channels from the survey.

property n_data: int

Number of data points.

sensitivity(conductivity: ndarray) → ndarray

Analytic d(response)/d(conductivity), averaged over the gate windows.

The same reduction the forward applies has to be applied to the Jacobian, or the two describe different data. The conductivity term comes from _compiled_sigma_jacobian() where it applies, and from SimPEG’s getJ otherwise.

class PyHydroGeophysX.forward.tdem_forward.TDEMSurveyConfig(source_location: ndarray = None, source_radius: float = 10.0, source_current: float = 1.0, receiver_location: ndarray = None, receiver_orientation: str = 'z', times: ndarray = None, waveform_type: str = 'step_off', source_turns: int = 1, source_moment: float | None = None, waveform_times: ndarray | None = None, waveform_currents: ndarray | None = None, gate_open: ndarray | None = None, gate_close: ndarray | None = None, gate_window: str = 'centre', gate_window_par: float = 0.667, waveform_period: float | None = None, waveform_repetitions: int = 3, analog_points_per_decade: int = 150, analog_model_points_per_decade: int = 40, analog_lowpass: dict | None = None, instrument_points_per_decade: int = 10, instrument_model_points_per_decade: int = 10, gate_quadrature_order: int = 8, receiver_type: str = 'b')

Bases: object

Configuration for TDEM survey geometry.

source_location

[x, y, z] location of source center (m)

Type:

numpy.ndarray

source_radius

Radius of circular loop source (m)

Type:

float

source_current

Peak current amplitude (A)

Type:

float

receiver_location

[x, y, z] location of receiver (m)

Type:

numpy.ndarray

receiver_orientation

Component to measure (‘x’, ‘y’, or ‘z’)

Type:

str

times

Time channels for measurement (s)

Type:

numpy.ndarray

waveform_type

Type of waveform (‘step_off’, ‘ramp_off’, ‘custom’)

Type:

str

analog_lowpass: dict | None = None

Analog receiver electronics parsed from a GEX file. The supported fields are receiver_damping / receiver_cutoff_hz for the receiver-coil two-pole filter and tib_order / tib_cutoff_hz for the transmitter-interface-board low-pass filter.

analog_model_points_per_decade: int = 40

Resolution of the grid SimPEG is asked to model, in samples per decade. Kept separate because SimPEG’s setup cost grows with it while the filter’s accuracy does not depend on it.

analog_points_per_decade: int = 150

Resolution of the internal grid the analog filter is integrated on, in samples per decade of time. See _analog_sampling().

gate_close: ndarray | None = None

gate_open: ndarray | None = None

Gate windows, for averaging the response over each one instead of reading it at the gate centre.

gate_quadrature_order: int = 8

Gauss points per smooth panel of the gate window. Eight integrates it to a few parts per million, negligible beside one SimPEG call.

gate_window: str = 'centre'

How a modelled decay becomes a gate value. See _gate_sampling().

gate_window_par: float = 0.667

Total cosine-taper fraction of the gate window, as stored in the project’s GateShapePar1. 0.667 leaves a 0.333 flat top.

instrument_model_points_per_decade: int = 10

SimPEG step-response density, in samples per decade. Ten matches the response grid and avoids asking SimPEG for redundant samples; raise it independently for convergence checks.

instrument_points_per_decade: int = 10

Response grid density, in samples per decade. Ten matches the grid the instrument’s own processing works on, with local interpolation between those nodes.

receiver_type: str = 'b'

“b”, or “dbdt” for its time derivative.

Type:

The recorded field

source_moment: float | None = None

Transmitter moment in A m^2, overriding current * turns * area. Set it to 1.0 for data already normalized by the transmitter moment, which is how TEMcompany instruments report dB/dt (V/A/m^4): dividing the measurement by the moment and then modelling with that moment counts it twice.

source_turns: int = 1

waveform_currents: ndarray | None = None

waveform_period: float | None = None

Half-period of the bipolar transmitter cycle, in seconds. Set it to model the earlier pulses of the train as well as the last one.

waveform_repetitions: int = 3

How many earlier half-cycles the repetition sums.

waveform_times: ndarray | None = None

Turn-off waveform as (times, currents) nodes. A real ramp is not a step, and the earliest gates of a ground system sit only microseconds after it ends, which is exactly where the difference bites.

PyHydroGeophysX.forward.tdem_forward.create_tdem_survey(times: ndarray, source_radius: float = 10.0, source_current: float = 1.0, source_location: ndarray | None = None, receiver_location: ndarray | None = None, receiver_orientation: str = 'z', waveform_type: str = 'step_off') → simpeg.electromagnetics.time_domain.Survey

Create a TDEM survey for 1D sounding.

Parameters:

• times – Time channels (s)

• source_radius – Loop radius (m)

• source_current – Peak current (A)

• source_location – Source center [x, y, z] (m)

• receiver_location – Receiver position [x, y, z] (m)

• receiver_orientation – Measurement component (‘x’, ‘y’, ‘z’)

• waveform_type – Waveform type (‘step_off’, ‘ramp_off’)

Returns:

SimPEG TDEM Survey object

PyHydroGeophysX.forward.tdem_forward.hydro_to_tdem(*args, **kwargs)

Deprecated alias for simulate_tdem_sounding_from_hydro().

PyHydroGeophysX.forward.tdem_forward.simulate_tdem_sounding_from_hydro(water_content: ndarray, porosity: ndarray, layer_thicknesses: ndarray, sigma_w: float | ndarray = 0.05, m: float | ndarray = 1.5, n: float | ndarray = 2.0, sigma_s: float | ndarray = 0.0, times: ndarray | None = None, source_radius: float = 10.0, noise_level: float = 0.05, seed: int | None = None, verbose: bool = False) → Tuple[ndarray, ndarray, ndarray, ndarray]

Convert hydrological properties to TDEM response.

This function takes water content and porosity from hydrological models and computes the expected TDEM response using petrophysical relationships.

Parameters:

• water_content – Water content for each layer (-)

• porosity – Porosity for each layer (-)

• layer_thicknesses – Thickness of each layer except bottom (m)

• sigma_w – Pore water conductivity (S/m)

• m – Cementation exponent

• n – Saturation exponent

• sigma_s – Surface conductivity (S/m)

• times – Time channels (s), default is logspace(-5, -2, 31)

• source_radius – Loop radius (m)

• noise_level – Relative noise level for synthetic data

• seed – Random seed for reproducibility

• verbose – Print progress information

Returns:

Tuple of (noisy_data, clean_data, uncertainties, conductivity)

Forward modeling utilities for Frequency-Domain Electromagnetic (FDEM) data.

Uses SimPEG’s frequency-domain module for 1D layered-earth simulations.

class PyHydroGeophysX.forward.fdem_forward.FDEMForwardModeling(thicknesses: ndarray, survey_config: FDEMSurveyConfig | None = None, survey: simpeg.electromagnetics.frequency_domain.Survey | None = None)

Bases: object

Forward modeling of Frequency-Domain EM data using SimPEG.

Supports 1D layered-earth conductivity models.

forward(conductivity: ndarray) → ndarray

Compute FDEM response for a given conductivity model.

One complex value per frequency, field and receiver, in that order: with receivers at (n, 3) locations the n values of a frequency (and field, for receiver_component="both") are consecutive.

forward_with_noise(conductivity: ndarray, noise_level: float = 0.05, seed: int | None = None) → Tuple[ndarray, ndarray, ndarray]

Compute noisy and clean FDEM responses with data uncertainties.

static hydro_to_fdem(water_content: ndarray, porosity: ndarray, layer_thicknesses: ndarray, **petro_params)

Convert hydrological properties to FDEM response via petrophysics.

class PyHydroGeophysX.forward.fdem_forward.FDEMSurveyConfig(source_location: ndarray = None, source_radius: float = 10.0, receiver_location: ndarray = None, receiver_orientation: str = 'z', receiver_component: str = 'secondary', frequencies: ndarray = None, waveform_type: str = 'dipole')

Bases: object

Configuration for FDEM survey geometry.

source_location is one XYZ point, shape (3,) or (1, 3), and is stored as (3,). receiver_location is one point, (3,), or several, (n, 3); it keeps the shape it was given. Several receivers share every frequency.

An omitted receiver sits at the origin, except where the origin cannot be modelled: below or above a dipole source (zero horizontal offset, a non-finite response) it goes 10 m along x from the source, and for a loop away from the origin it goes to the loop centre on z = 0, the only receiver position the layered-earth loop kernel supports. An explicit zero horizontal dipole offset is rejected before simulation.

frequencies: ndarray = None

receiver_component: str = 'secondary'

receiver_location: ndarray = None

receiver_orientation: str = 'z'

source_location: ndarray = None

source_radius: float = 10.0

waveform_type: str = 'dipole'