AEM Inversion & Forward Modelling

Documents the inversion and forward modelling processes for airborne electromagnetics (TEM), powered by a custom fork of SimPEG.

Source Files

System Description Architecture

The SimPEG/electromagnetics/utils/static_instrument/ directory defines a hierarchy of instrument descriptions for moving EM platforms (AEM, TTEM). The key design principle is that the geometric configuration of transmitter(s) and receiver(s) is independent of the data — only their absolute positions change per-sounding.

XYZSystem (base class)

Source: base.py

Base class for all AEM system descriptions. Core responsibilities:

  1. System instantiation: Takes a libaarhusxyz.XYZ object and optional keyword parameters (which override class-level defaults via __getattribute__)
  2. Filtering: Applies sounding_filter (row-wise) and gate_filter (column-wise) via xyzfilter.FilteredXYZ
  3. Survey construction: make_survey() iterates over all soundings, calling make_system() to create SimPEG sources for each
  4. Inversion assembly: Chain of make_thicknesses() → make_misfit() → make_regularization() → make_optimizer() → make_inversion()
  5. Model <-> XYZ conversion: inverted_model_to_xyz() and forward_data_to_xyz() convert between SimPEG model vectors and libaarhusxyz format

SingleMomentTEMXYZSystem

Source: single.py

A simple single-moment system suitable for synthetic data:

DualMomentTEMXYZSystem

Source: dual.py

Dual-moment system describing SkyTEM instruments. Cannot be instantiated directly — use DualMomentTEMXYZSystem.load_gex(gex) to generate a subclass from a GEX instrument file.

Key features: - Two sources: Low Moment (Channel 1) and High Moment (Channel 2), each as a MagDipole with separate PiecewiseLinearWaveform from the GEX file - Receiver offset: Rx position relative to Tx from gex.General['RxCoilPosition'] - Tilt correction: correct_tilt_pitch_for1Dinv() — same (cos·cos)² factor as the pipeline - Gate factor scaling: Data multiplied by gex.ChannelN['GateFactor'] (Aarhus Workbench convention) - Automatic sounding filter: Excludes soundings where no usable gates remain in either channel - Gate filtering: Configurable gate ranges via gate_filter__start_lm/end_lm and gate_filter__start_hm/end_hm

Inversion Formulation

Forward Model

Uses SimPEG.electromagnetics.time_domain.Simulation1DLayeredStitched — a 1D layered earth EM solution computed independently per sounding, stitched into a pseudo-2D model. For each sounding, the solver computes the vertical magnetic field response of a layered half-space to a magnetic dipole source with an arbitrary waveform.

The 1D solution uses the Hankel transform for the frequency-domain solution (via digital filter coefficients) and the Fourier transform (via cosine/sine transforms) to convert to time domain. Layer conductivities are transformed using maps.ExpMap (σ = exp(m)).

Solver options: - Pardiso (direct sparse, multi-threaded) — requires pymatsolver.PardisoSolver - Default spLU (SuperLU) — fallback

Parallelism: Configurable n_cpu threads for sounding-parallel computation.

Model Discretization

Layer thicknesses: Three options via startmodel__thicknesses_type: 1. Log-spaced (default): 30 layers, 1 m → 400 m depth range. build_log_spaced_layer_thick() iteratively adjusts the last layer to match the target depth exactly. 2. Geometric: Constant geometric factor (default 1.15309) between successive layers. 3. Time-based: Derived from gate times and background conductivity using get_vertical_discretization_time().

The model vector has n_param = n_layers × n_soundings parameters (plus one half-space). Resistivity is log-transformed: σ = exp(m), initialized to 1/startmodel__res (default 100 Ω·m).

Data Misfit

Weighted L2 norm with weights = 1 / uncertainty:

φ_d = ||W_d · (d_obs - d_pred)||²

where the per-datum uncertainty combines: - STD from stacking (STD_Ch## values, capped to minimum uncertainties__std_data = 3%) - Time-dependent noise floor: noise_level_1ms × (t × 1000)^(noise_exponent) / moment - W_d = diag(1 / (std × |data| + noise))

Regularization

LaterallyConstrainedInversion (LCI) using two-mesh approach:

Reference model: m_ref = 1/100 S/m (uniform 100 Ω·m half-space).

Reference: Auken et al. (2002) for laterally constrained inversion theory.

Inversion Algorithm

  1. Beta estimation: BetaEstimate_ByEig — computes initial regularization weight from eigenvalue spectrum of JᵀJ
  2. Beta schedule: BetaSchedule(coolingFactor=2, coolingRate=1) — halves beta every iteration
  3. Target misfit: Stopping criterion when φ_d ≈ N_data
  4. Optimization: InexactGaussNewton — inexact Newton with conjugate-gradient inner iterations (max 20 CG per step, max 40 GN iterations total)

Sparse Model (IRLS)

Optional iterative reweighted least squares (IRLS) produces a "focused" model with sharper boundaries:

IRLS iteration: update weights W_R^(k+1) = diag(1 / (|m_j - m_ref|^(2-p/2) + ε))

where p ≈ 1 gives blocky/sparse solutions. Runs after L2 convergence.

References: Last & Kubik (1983), Farquharson & Oldenburg (1998).

Forward Modelling

Source: forward_process.py

Takes a resistivity model (libaarhusxyz format with resistivity, dep_top, dep_bot) and computes synthetic dB/dt responses using the same system description and Simulation1DLayeredStitched.

The model is passed directly as log-conductivity: m = ln(1/resistivity). The output is a libaarhusxyz dataset with: - Synthetic Gate_Ch## data - InUse_Ch## flags set to 1 where data is finite - STD_Ch## set to GEX UniformDataSTD - Current_Ch## filled from GEX approximate current - Tilt columns defaulted to zero

Result Conversion

inverted_model_to_xyz()

Converts the inversion result back to libaarhusxyz format:

resistivity = 1 / exp(model.reshape(n_soundings, n_layers))

forward_data_to_xyz()

Converts predicted data to libaarhusxyz format, computing per-sounding residuals and RMS misfit:

resdata[i] = sqrt(mean(derr[i, :]²))  where derr = (d_obs - d_pred) × W_d

Custom Directives

Source: directives.py

References