Physics-informed 2-D inversion#

pycsamt.ai.inversion.PINNInverter2D jointly optimizes layered models for all stations on a profile. It uses a differentiable 1-D MT forward recursion at each station and couples adjacent station models with lateral smoothness. The result is a pseudo-2-D resistivity section: a set of jointly regularized 1-D models, not a full finite-element or finite-difference 2-D EM forward inversion.

This distinction is central to scientific use. The workflow can provide a useful profile estimate when locally layered physics and lateral smoothness are reasonable approximations. It should not be presented as equivalent to Occam2D, MARE2DEM, ModEM, or another solver that explicitly models 2-D/3-D EM fields.

Meaning of “2-D”

“2-D” here describes the jointly optimized station–layer section and its lateral regularization. Data fit is calculated with the differentiable Wait-style MT 1-D recursion independently at each station. Strong lateral induction, off-profile structure, anisotropy, topography, and other missing physics can therefore produce a smooth but biased section.

When to use this workflow#

Consider PINN 2-D when:

  • stations form an ordered survey profile;

  • dimensionality evidence suggests locally 1-D or gently varying 2-D structure;

  • a common layered parameterization is scientifically acceptable;

  • no labelled model dataset is available or desired;

  • direct response-space optimization is preferred to supervised prediction;

  • a fast pseudo-2-D baseline or starting model is useful;

  • the result will be compared with classical inversion and independent data.

Prefer a full 2-D or 3-D classical solver when:

  • lateral induction is important to the target;

  • strong contacts, topography, coast effects, anisotropy, or 3-D structure are expected;

  • tensor components cannot be represented by per-station 1-D responses;

  • model appraisal requires the corresponding numerical forward physics;

  • regulatory or project standards require a recognized 2-D/3-D solver.

Workflow#

  1. QC and order the profile data;

  2. review strike and dimensionality evidence;

  3. choose TE, TM, or the current averaged both mode;

  4. define the common frequency grid and supported depth;

  5. choose layer count, smoothness weights, optimizer, and epochs;

  6. run baseline and sensitivity configurations;

  7. inspect convergence and response residuals;

  8. compare section, thicknesses, and reconstructed responses across scenarios;

  9. validate against classical inversion and independent evidence;

  10. report the pseudo-2-D physics and all limitations explicitly.

1. Understand the model parameterization#

For S stations and L layers, the optimizer uses:

log_rho.shape   = (S, L)
log_thick.shape = (S, L - 1)

Both arrays are optimized jointly. Resistivity is represented as \(\log_{10}(\rho/\Omega\mathrm{m})\); thickness is represented internally in log10 metres and converted back for output.

The public section methods transpose the internal station-major layout:

resistivity_section.shape = (L, S)
thickness_section.shape   = (L - 1, S)

Each station has its own interface thicknesses. Consequently, layer index is a parameter correspondence, not automatically a continuous geological horizon. Inspect cumulative interface depths before drawing boundaries between stations.

depth_max guides the layered depth parameterization and initialization; it does not establish a geophysical depth of investigation. The defensible interpretation depth must come from bandwidth, sensitivity, response fit, and scenario stability.

2. Understand the loss function#

The implemented objective combines data fit, vertical smoothing, and lateral smoothing:

\[\mathcal L = \mathcal L_{data} + \lambda_z\mathcal R_z + \lambda_x\mathcal R_x.\]

Data term#

For valid station–frequency cells, the data term combines squared differences in log10 apparent resistivity and normalized phase:

\[\mathcal L_{data} = \frac{1}{N_v}\sum_{s,f} \left[ (\log_{10}\rho_{a,sf}^{pred}-\log_{10}\rho_{a,sf}^{obs})^2 + \left(\frac{\phi_{sf}^{pred}-\phi_{sf}^{obs}}{90}\right)^2 \right].\]

N_v is the number of finite observations after common-grid interpolation. Missing cells are excluded from the data term.

The 90-degree phase normalization makes the two residual blocks numerically comparable but is not an observational covariance model. Formal errors and station-specific uncertainty are not explicit weights in this objective.

Vertical regularization#

smoothness_weight is \(\lambda_z\). It penalizes changes in log resistivity between adjacent layers at each station. Larger values favor vertically smooth models; smaller values permit sharper changes and greater instability.

Lateral regularization#

lateral_weight is \(\lambda_x\). It penalizes log-resistivity differences between adjacent stations at the same layer index. It assumes the input station order represents profile adjacency.

The lateral term does not currently scale differences by physical station spacing. Closely and widely spaced station pairs are treated as adjacent indices. Resample carefully or interpret regularization strength in light of irregular chainage.

Thickness regularization and constraints#

Resistivity smoothness is explicit in the documented objective. Layer thicknesses are optimized in log space, but the same lateral/vertical regularization description should not be assumed to apply to them unless confirmed in the backend kernel. Audit thickness stability separately.

3. Prepare and order profile data#

Load the profile canonically and verify station order before constructing the inverter:

from pycsamt.emtools._core import ensure_sites

sites = ensure_sites(
    "data/AMT/WILLY_DATA/L18",
    recursive=True,
    verbose=0,
)

The inverter uses pycsamt.ai.inversion.sites_to_obs_2d() internally and preserves the order returned by the input. It does not sort stations by chainage. Provide a Sites object already ordered along the intended profile direction.

Retain a table containing:

  • station name;

  • profile distance in metres;

  • easting, northing, and elevation;

  • original and inversion order;

  • valid frequency range;

  • TE/TM component completeness;

  • QC exclusions and corrections.

Do not infer physical x coordinates from array column number in a final plot.

4. Select polarization mode#

mode accepts "te", "tm", or "both".

te

Uses the component named by comp_te, default "xy".

tm

Uses the component named by comp_tm, default "yx". The public observation bridge stores TM phase as magnitude.

both

In the current implementation, averages TE and TM apparent resistivity and averages TE and TM phase before optimization. It does not form two independent data-misfit blocks or jointly fit both modes separately.

Warning

Averaging TE and TM can suppress genuine mode differences caused by 2-D/3-D structure. Use mode="both" only as an explicitly reviewed approximation. For scientific comparison, run TE and TM separately and inspect their disagreement.

Mode labels depend on profile orientation and strike convention. Verify these before accepting the defaults.

5. Choose the common frequency grid#

The constructor extracts valid observations and builds a log-spaced common grid with n_freqs points spanning the combined available range:

from pycsamt.ai.inversion import PINNInverter2D

inversion = PINNInverter2D(
    sites,
    n_layers=10,
    depth_max=2000.0,
    n_freqs=32,
    mode="te",
)

Each station is interpolated to this grid in log-frequency space. Values outside a station’s own range become NaN and do not contribute to the data loss.

Because the public constructor does not expose explicit freq_min and freq_max arguments, the common endpoints follow the extracted survey coverage. To enforce a project-specific range, prefilter the reviewed data or extend the API through a tested project change rather than mutating private attributes.

Inspect coverage per station. A station constrained by only a small subset of the common grid can borrow visual smoothness from neighbors without having the same depth sensitivity.

6. Choose layer count and depth#

Start with the simplest parameterization able to represent expected electrical units. Test several defensible combinations, for example:

candidates = [
    {"n_layers": 6, "depth_max": 1500.0},
    {"n_layers": 10, "depth_max": 2000.0},
    {"n_layers": 14, "depth_max": 3000.0},
]

Increasing layers can reduce response residuals while adding unstable interfaces. Increasing depth_max can allocate parameters below meaningful sensitivity. Compare:

  • reconstructed response fit;

  • interface-depth stability;

  • deep-layer variation across runs;

  • dependence on smoothing;

  • classical inversion and sensitivity depth;

  • independent borehole or geological evidence.

Do not interpret the bottom of the configured model as the investigation depth.

7. Select regularization weights#

The defaults are:

smoothness_weight = 0.01
lateral_weight    = 0.005

Treat them as starting values, not universal choices. Run a grid of plausible weights:

scenarios = [
    (0.001, 0.0005),
    (0.01, 0.005),
    (0.05, 0.02),
]

For each scenario, retain total loss, response residuals, roughness measures, section differences, and target stability. The reported convergence curve contains total loss, not separate data and regularization histories, so additional diagnostics may require calculation from outputs or backend instrumentation.

Avoid choosing weights solely for a visually smooth section. Excessive lateral smoothing can smear faults, isolated conductors, or station-specific problems. Insufficient smoothing can fit noise with alternating layers.

8. Select optimizer controls#

epochs

Maximum Adam iterations, default 300.

lr

Adam learning rate, default 1e-2.

device

PyTorch/TensorFlow device resolved through the active backend.

The current high-level fit() does not expose early stopping, validation data, learning-rate scheduling, or gradient clipping. It runs the requested optimization iterations through the backend kernel. Monitor the curve and run repeat configurations rather than assuming the final epoch is optimal.

A large learning rate can oscillate or diverge; a very small rate may appear stable without reaching an adequate fit. Compare final sections across epochs and learning rates.

9. Run a baseline inversion#

from pycsamt.ai.inversion import PINNInverter2D

inversion = PINNInverter2D(
    sites,
    n_layers=10,
    depth_max=2000.0,
    n_freqs=32,
    mode="te",
    smoothness_weight=0.01,
    lateral_weight=0.005,
    epochs=300,
    lr=1e-2,
    comp_te="xy",
    comp_tm="yx",
    device=None,
    verbose=0,
)

inversion.fit(verbose=True, log_every=25)
print(inversion)
print(inversion.stations)
print(inversion.n_sites)

The constructor performs data extraction immediately and can raise before fit() when no valid observations survive. The deep-learning backend is required when fitting.

10. Extract resistivity and thickness#

log10_rho = inversion.resistivity_section(as_log10=True)
rho_ohm_m = inversion.resistivity_section(as_log10=False)
thickness_m = inversion.thickness_section()

print(log10_rho.shape)
print(thickness_m.shape)

Build interface depths per station from cumulative thickness:

import numpy as np

interface_depth_m = np.cumsum(thickness_m, axis=0)

assert np.all(np.isfinite(rho_ohm_m))
assert np.all(rho_ohm_m > 0)
assert np.all(np.isfinite(thickness_m))
assert np.all(thickness_m > 0)

Check predictions against explicit parameter bounds used in project scenarios. The optimizer’s positivity through log representation does not guarantee geological plausibility.

11. Review convergence#

curve = inversion.convergence_curve()
print(curve.tail())

The returned pandas DataFrame contains epoch and total loss.

Review:

  • initial reduction rate;

  • late oscillation or divergence;

  • whether loss is still decreasing strongly at the final epoch;

  • sensitivity to learning rate;

  • repeatability across backend/device and repeated runs;

  • whether lower loss corresponds to better response residuals and stable structure.

There is no separate validation loss because the optimization is observation-specific rather than supervised dataset training. Overfitting can still occur by fitting noise with an overly flexible or weakly regularized model.

12. Review residuals#

residuals = inversion.residuals()
residuals["log_rho_residual"] = (
    np.log10(residuals["rho_pred"])
    - np.log10(residuals["rho_obs"])
)
residuals["phase_residual_deg"] = (
    residuals["phase_pred"] - residuals["phase_obs"]
)

summary = residuals.groupby("station").agg(
    rho_rms=("log_rho_residual",
             lambda values: float(np.sqrt(np.nanmean(values**2)))),
    phase_rms_deg=("phase_residual_deg",
                   lambda values: float(np.sqrt(np.nanmean(values**2)))),
)
print(summary)

The residual table uses each station’s original frequencies, not only the common optimization grid.

Warning

For mode="both", the optimizer fits averaged TE/TM observations, while residuals() selects TE observations for its reported observed columns. The diagnostic therefore does not directly report residuals against the averaged data used in optimization. Run TE and TM separately or construct a custom audited residual table for the averaged mode.

Also verify the shapes/components returned by the forward response in the installed version. A single aggregate RMS can hide station-, frequency-, or mode-dependent mismatch.

13. Plot the section correctly#

The section is defined by varying thicknesses, so imshow on layer index can misrepresent depth. For a quick layer-index diagnostic:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 5))
image = ax.imshow(
    log10_rho,
    aspect="auto",
    origin="upper",
    interpolation="nearest",
)
ax.set_xlabel("Station index")
ax.set_ylabel("Layer index")
fig.colorbar(image, ax=ax, label="log10 resistivity (ohm m)")
fig.savefig("review/pinn2d_layer_index.png", dpi=200,
            bbox_inches="tight")

Label this explicitly as layer index. For a physical-depth section, resample each layered model onto a common depth grid using its cumulative thicknesses, then plot against reviewed profile distances. Preserve the original layered outputs alongside any resampled image.

Do not connect layer boundaries as geological horizons without checking their meaning and stability.

14. Run TE/TM and regularization scenarios#

A defensible study includes separate modes:

results = {}
for mode in ("te", "tm"):
    inv = PINNInverter2D(
        sites,
        n_layers=10,
        depth_max=2000.0,
        n_freqs=32,
        mode=mode,
        smoothness_weight=0.01,
        lateral_weight=0.005,
        epochs=300,
        lr=1e-2,
    ).fit(verbose=False)
    results[mode] = {
        "log10_rho": inv.resistivity_section(),
        "thickness_m": inv.thickness_section(),
        "residuals": inv.residuals(),
    }

Compare:

  • structures common to TE and TM;

  • mode-specific conductive or resistive features;

  • station and frequency residual patterns;

  • boundary-depth differences;

  • sensitivity to \(\lambda_z\) and \(\lambda_x\);

  • changes when suspicious stations or bands are excluded.

Large disagreement may indicate dimensionality, distortion, component quality, or parameterization problems rather than a need to average modes.

15. Compare with 1-D and classical 2-D results#

The lateral weight should add demonstrable value over independent station models. Compare against:

  • pycsamt.ai.inversion.PINNInverter1D at each station;

  • a no/very-low lateral-weight scenario;

  • a conventional layered 1-D inversion;

  • Occam2D when a 2-D line assumption is appropriate;

  • ModEM or MARE2DEM when geometry and resources justify them;

  • boreholes and mapped geology.

Compare response fit using compatible error models. The PINN objective’s phase normalization and missing-value mask need not match the classical solver’s RMS. Do not rank methods by incomparable scalar values.

16. Use the PINN agent for orchestration#

pycsamt.agents.PINNInversionAgent wraps the workflow and returns an pycsamt.agents.AgentResult:

from pycsamt.agents import PINNInversionAgent

agent = PINNInversionAgent(
    dim=2,
    n_layers=10,
    depth_max=2000.0,
    smoothness_weight=0.01,
    lateral_weight=0.005,
    epochs=300,
    lr=1e-2,
    solver="mt1d",
)

result = agent.execute({
    "sites": sites,
    "output_dir": "outputs/pinn2d/L18",
})

if result.status == "failed":
    raise RuntimeError(result.error)

section = result["section"]
loss_table = result.get("loss_df")
residual_table = result.get("residuals_df")

The agent is convenient for figures and standardized results, but inspect its constructor mapping: the generic agent interface does not expose every PINNInverter2D option such as mode and component selection. Use the lower- level class when those controls matter.

17. Consider hybrid 2-D refinement#

pycsamt.ai.inversion.HybridInverter2D uses a fitted supervised pycsamt.ai.inversion.EMInverter2D to initialize the same joint physics refinement:

from pycsamt.ai.inversion import HybridInverter2D

hybrid = HybridInverter2D(
    sites,
    ai_inverter=fitted_ai_2d,
    n_layers=10,
    depth_max=2000.0,
    n_freqs=32,
    mode="te",
    smoothness_weight=0.005,
    lateral_weight=0.003,
    epochs=150,
    lr=5e-3,
).fit()

stage1 = hybrid.stage1_section()
stage2 = hybrid.resistivity_section()
residuals_stage1 = hybrid.residuals(stage=1)
residuals_stage2 = hybrid.residuals(stage=2)

The hybrid is justified only if Stage 2 improves response fit or scientific stability without introducing unsupported structure. Preserve both stages.

Warning

HybridInverter2D accepts an inverter object or a checkpoint path, but the current EMInverter2D class does not expose the explicit public save()/load() pair shown by some older examples/docstrings. Prefer a fitted in-memory object or a project-tested persistence mechanism; do not assume arbitrary .npz restoration works without verification.

18. Assess uncertainty#

PINNInverter2D returns one optimized solution and a total loss history. It does not itself return a posterior or calibrated interval.

Assess uncertainty with scenario ensembles:

  • TE versus TM;

  • layer count and depth_max;

  • vertical and lateral smoothness weights;

  • learning rate and epochs;

  • data masks and station exclusions;

  • alternative initialization, including hybrid Stage 1;

  • acceptable classical inversion configurations;

  • data perturbations consistent with measurement errors.

Summarize model spread on a common physical depth grid and distinguish optimizer variability from structural and data uncertainty. See AI inversion uncertainty for the wider framework.

19. Preserve a run record#

pinn2d/L18_te_r001/
├── manifest.yml
├── input/
│   ├── station_order.csv
│   ├── qc_reference.yml
│   └── frequency_coverage.csv
├── configuration/
│   └── pinn2d.yml
├── outputs/
│   ├── log10_resistivity.npy
│   ├── thickness_m.npy
│   ├── interface_depth_m.npy
│   ├── convergence.csv
│   └── residuals.csv
├── figures/
└── review/
    ├── scenario_comparison.csv
    └── scientific_review.md

Record backend, device, software version, station order, modes/components, frequency grid, layer/depth parameterization, weights, optimizer, epochs, residual definitions, scenario set, classical comparison, reviewer, and status.

Complete example#

from pathlib import Path
import json
import numpy as np

from pycsamt.ai.inversion import PINNInverter2D
from pycsamt.emtools._core import ensure_sites

output = Path("pinn2d/L18_te_r001")
(output / "outputs").mkdir(parents=True, exist_ok=True)

sites = ensure_sites(
    "data/AMT/WILLY_DATA/L18",
    recursive=True,
    verbose=0,
)

inversion = PINNInverter2D(
    sites,
    n_layers=10,
    depth_max=2000.0,
    n_freqs=32,
    mode="te",
    smoothness_weight=0.01,
    lateral_weight=0.005,
    epochs=300,
    lr=1e-2,
    comp_te="xy",
    comp_tm="yx",
    verbose=0,
).fit(verbose=True, log_every=25)

log10_rho = inversion.resistivity_section()
thickness_m = inversion.thickness_section()
interface_depth_m = np.cumsum(thickness_m, axis=0)
convergence = inversion.convergence_curve()
residuals = inversion.residuals()

np.save(output / "outputs" / "log10_resistivity.npy", log10_rho)
np.save(output / "outputs" / "thickness_m.npy", thickness_m)
np.save(output / "outputs" / "interface_depth_m.npy", interface_depth_m)
convergence.to_csv(output / "outputs" / "convergence.csv", index=False)
residuals.to_csv(output / "outputs" / "residuals.csv", index=False)

manifest = {
    "workflow": "PINNInverter2D",
    "physics": "per_station_mt1d_with_lateral_smoothing",
    "mode": "te",
    "component": "xy",
    "stations": inversion.stations,
    "n_layers": 10,
    "depth_max_m": 2000.0,
    "n_freqs": 32,
    "smoothness_weight": 0.01,
    "lateral_weight": 0.005,
    "epochs": 300,
    "learning_rate": 1e-2,
}
(output / "manifest.json").write_text(
    json.dumps(manifest, indent=2),
    encoding="utf-8",
)

Review checklist#

Check

Required evidence

Pseudo-2-D scope is explicit

Per-station MT1D physics, lateral coupling, and missing full 2-D effects.

Profile order is verified

Station names, chainage, coordinates, direction, gaps, and duplicates.

Modes are defensible

Strike convention, TE/TM components, separate results, and any averaging.

Frequency coverage is adequate

Valid cells by station, common grid, missingness, and supported depth.

Parameterization is tested

Layers, depth, thickness stability, and interface-depth scenarios.

Regularization is appraised

Vertical/lateral weight grid, roughness, response fit, and stable targets.

Optimization converges

Loss curve, learning-rate/epoch sensitivity, repeats, and no divergence.

Residuals are reviewed

Apparent resistivity and phase by station/frequency, not only one RMS.

Baselines are included

Independent 1-D, classical 2-D/3-D where appropriate, and geology.

Uncertainty is conditional

Data perturbation, configuration scenarios, initialization, and omitted physics.

Run is reproducible

Input/QC IDs, order, configuration, arrays, backend, residuals, reviewer, and status.

Common mistakes#

Avoid these errors:

  • describing the workflow as a full numerical 2-D EM inversion;

  • allowing arbitrary input order to define lateral neighbors;

  • treating layer index as physical depth;

  • calling depth_max the depth of investigation;

  • selecting both mode without recognizing TE/TM averaging;

  • assuming lateral smoothing accounts for irregular station spacing;

  • selecting weights for visual smoothness alone;

  • interpreting total loss as pure data misfit;

  • reporting the both residual table as residuals against averaged inputs;

  • comparing PINN loss directly with a differently weighted classical RMS;

  • assuming a single optimized solution provides uncertainty;

  • discarding thickness outputs when plotting resistivity;

  • relying on outdated Hybrid/EMInverter2D checkpoint examples without testing persistence.

Next steps#

Continue with: