Groom-Bailey Galvanic Distortion#

pycsamt.emtools.gb estimates and optionally removes frequency-independent galvanic distortion from MT/AMT/CSAMT impedance tensors. It is designed as an auditable preprocessing step before 2-D interpretation and inversion.

The fitted model is:

\[Z_{obs}(f) \approx D Z_{2D}(f)\]

where:

  • Z_obs is the observed impedance tensor;

  • D is a real, frequency-independent 2 x 2 distortion matrix;

  • Z_2D is the best anti-diagonal regional tensor at each frequency.

Full callable signatures live in the API reference. This page explains how to fit the table, read the distortion parameters, apply the correction, and record the result in a pre-2D workflow.

When To Use Groom-Bailey#

Use this workflow when the data appear close enough to 2-D for galvanic distortion correction to be meaningful, but the impedance tensor has diagonal leakage or station-dependent distortion that should be documented before inversion.

Good use cases include:

  • preparing a 2-D inversion input after dimensionality and strike checks;

  • testing whether diagonal tensor leakage is reduced after correction;

  • documenting twist, shear, and anisotropy-style distortion parameters;

  • comparing corrected and uncorrected impedance curves at the same station.

Poor use cases include:

  • strongly 3-D data with no stable strike or 2-D period band;

  • too few valid frequencies in the selected band;

  • using the fitted gain as a unique static-shift solution;

  • applying correction without saving the fit diagnostics.

Core Assumptions#

The implementation fits a real distortion matrix that is constant over the selected period band. That is the galvanic assumption: the distortion is local and frequency-independent, while the regional tensor varies with frequency.

The regional tensor is forced to be anti-diagonal:

\[\begin{split}Z_{2D}(f) = \begin{bmatrix} 0 & u(f) \\ v(f) & 0 \end{bmatrix}\end{split}\]

The observed tensor is then approximated by multiplying this 2-D tensor by D. The fit alternates between estimating the anti-diagonal regional tensor and solving for the rows of the real distortion matrix.

The fitted matrix is normalized by determinant scale, then summarized as gain, twist, shear, and anisotropy-style parameters.

Fit A Distortion Table#

Start by estimating parameters without changing the data.

 1from pycsamt.emtools.gb import groom_bailey_table
 2
 3survey = "data/AMT/WILLY_DATA/L18PLT"
 4
 5table = groom_bailey_table(
 6    survey,
 7    band=(1e-3, 10.0),
 8    rotate_deg=None,
 9    min_freq=4,
10    max_iter=30,
11    tol=1e-6,
12    robust=True,
13)
14
15print(
16    table[
17        [
18            "station",
19            "status",
20            "n_freq",
21            "twist_deg",
22            "shear",
23            "anisotropy",
24            "rms_fit",
25            "diagonal_ratio_before",
26            "diagonal_ratio_after",
27        ]
28    ].head()
29)

   station status  ...  diagonal_ratio_before  diagonal_ratio_after
0  18-001A     ok  ...               0.614153              0.273213
1  18-002U     ok  ...               0.434659              0.261954
2  18-003A     ok  ...               0.373588              0.379570
3  18-004A     ok  ...               0.454007              0.323360
4  18-005U     ok  ...               0.496585              0.305364

[5 rows x 9 columns]

The band argument is in period seconds, not hertz. Choose a band that is justified by dimensionality, strike stability, and data quality.

Table Columns#

Successful rows have status == "ok" and include:

Column

Meaning

station

Station name.

n_freq

Number of valid frequencies used in the fit.

period_min_s / period_max_s

Period range actually used after band selection.

rotate_deg

Rotation angle applied before fitting, or NaN.

distortion_xx … distortion_yy

Entries of the fitted real 2 x 2 distortion matrix.

gain

Matrix scale from the determinant normalization.

twist_deg

Twist angle inferred from the normalized matrix.

shear

Dimensionless shear-style parameter, clipped to [-0.99, 0.99].

shear_angle_deg

atan(shear) in degrees.

anisotropy

Dimensionless anisotropy-style parameter, clipped to [-0.99, 0.99].

rms_fit

Relative fit residual.

diagonal_ratio_before

Median diagonal/off-diagonal tensor ratio before correction.

diagonal_ratio_after

Median diagonal/off-diagonal tensor ratio after applying the fitted inverse matrix to the fitted band.

robust

Whether robust residual weighting was used.

method

Current method label, gb_real_distortion_2d.

Rows with too few valid frequencies have status == "insufficient_frequencies" and include the available n_freq. Increase the band, lower min_freq only with care, or exclude that station from correction.

Reading The Parameters#

The most useful diagnostic columns are usually:

  • rms_fit: lower values indicate that the fitted model describes the selected band better.

  • diagonal_ratio_before and diagonal_ratio_after: correction is behaving sensibly when the after value is lower.

  • twist_deg: large twist can imply strong galvanic distortion or a poor 2-D assumption.

  • shear and anisotropy: large absolute values deserve station inspection.

  • n_freq: low values make the fit less stable.

Do not over-interpret gain as a unique static-shift estimate. The scalar gain ambiguity is not uniquely resolved by this decomposition.

Rank Stations For Review#

Use the table to find stations with poor fits or strong residual diagonal leakage.

 1from pycsamt.emtools.gb import groom_bailey_table
 2
 3table = groom_bailey_table(
 4    "data/AMT/WILLY_DATA/L18PLT",
 5    band=(1e-3, 10.0),
 6    robust=True,
 7)
 8
 9ok = table.loc[table["status"] == "ok"].copy()
10ok["diag_reduction"] = (
11    ok["diagonal_ratio_before"] - ok["diagonal_ratio_after"]
12)
13
14ranked = ok.sort_values(
15    ["rms_fit", "diagonal_ratio_after"],
16    ascending=[False, False],
17)
18
19print(
20    ranked[
21        [
22            "station",
23            "n_freq",
24            "rms_fit",
25            "diag_reduction",
26            "twist_deg",
27            "shear",
28            "anisotropy",
29        ]
30    ].head(10)
31)

    station  n_freq   rms_fit  diag_reduction  twist_deg     shear  anisotropy
22  18-022U      39  0.552447        0.322411  22.141866 -0.117145   -0.008538
21  18-021U      39  0.537890       -0.419458 -63.216543  0.990000    0.404853
20  18-021B      39  0.430286       -0.604539 -37.911526  0.866936    0.211520
17  18-018A      39  0.427220       -0.564491  56.099182  0.990000   -0.234017
24  18-023A      39  0.412376        0.298014  19.376516 -0.123178    0.012153
27  18-025A      39  0.390159       -0.384625 -62.091367  0.504772    0.012596
18  18-019U      39  0.363720       -0.016977  17.861709 -0.015653    0.008229
8   18-009A      39  0.335742       -0.014381  17.904448  0.118619   -0.006990
19  18-020A      39  0.326943       -0.391908 -59.746000  0.870488    0.525820
12  18-013U      39  0.322409        0.166075   8.996415 -0.388653    0.079927

Stations with high rms_fit or little diagonal reduction should be reviewed before applying correction automatically.

Use A Strike Rotation#

If you have selected a strike angle, pass it as rotate_deg before fitting.

 1from pycsamt.emtools.gb import groom_bailey_table
 2
 3strike_deg = 35.0
 4
 5table = groom_bailey_table(
 6    "data/AMT/WILLY_DATA/L18PLT",
 7    band=(1e-3, 10.0),
 8    rotate_deg=strike_deg,
 9    robust=True,
10)

The rotation is applied to the tensor before fitting the distortion matrix. Use a strike that has been justified by the strike and dimensionality workflows, not one chosen to improve the GB fit alone.

Apply A Precomputed Table#

Use apply_groom_bailey when you have already inspected and accepted a table.

 1from pycsamt.emtools.gb import apply_groom_bailey, groom_bailey_table
 2
 3survey = "data/AMT/WILLY_DATA/L18PLT"
 4
 5table = groom_bailey_table(
 6    survey,
 7    band=(1e-3, 10.0),
 8    robust=True,
 9)
10
11accepted = table.loc[
12    (table["status"] == "ok")
13    & (table["rms_fit"] < 0.25)
14    & (table["diagonal_ratio_after"] < table["diagonal_ratio_before"])
15].copy()
16
17corrected = apply_groom_bailey(
18    survey,
19    table=accepted,
20    inplace=False,
21)

Only stations present in the accepted table are corrected. Stations missing from the table or with invalid matrices are left unchanged.

Estimate And Apply In One Step#

Use groom_bailey_decomposition when you want a result container with the fitted table and optionally corrected sites.

 1from pycsamt.emtools.gb import groom_bailey_decomposition
 2
 3result = groom_bailey_decomposition(
 4    "data/AMT/WILLY_DATA/L18PLT",
 5    apply=True,
 6    band=(1e-3, 10.0),
 7    rotate_deg=None,
 8    robust=True,
 9    inplace=False,
10)
11
12print(result.summary())
13corrected_sites = result.sites
14gb_table = result.table

GroomBaileyResult(stations=28, applied=True, median_rms=0.2797)

The result container records:

  • sites: corrected sites when apply=True, otherwise loaded sites.

  • table: fitted parameter table.

  • applied: whether correction was applied.

  • method: method label.

  • n_station: number of fitted table rows.

Compare Robust And Non-Robust Fits#

Robust weighting downweights high-residual frequencies during fitting. Compare both modes when outliers are suspected.

 1from pycsamt.emtools.gb import groom_bailey_table
 2
 3survey = "data/AMT/WILLY_DATA/L18PLT"
 4
 5robust = groom_bailey_table(
 6    survey,
 7    band=(1e-3, 10.0),
 8    robust=True,
 9    api=False,
10)
11plain = groom_bailey_table(
12    survey,
13    band=(1e-3, 10.0),
14    robust=False,
15    api=False,
16)
17
18compare = robust.merge(
19    plain,
20    on="station",
21    suffixes=("_robust", "_plain"),
22)
23
24print(
25    compare[
26        [
27            "station",
28            "rms_fit_robust",
29            "rms_fit_plain",
30            "twist_deg_robust",
31            "twist_deg_plain",
32        ]
33    ].head()
34)

   station  rms_fit_robust  rms_fit_plain  twist_deg_robust  twist_deg_plain
0  18-001A        0.138687       0.138478         10.728457        10.602753
1  18-002U        0.212166       0.211317          8.600826         9.526479
2  18-003A        0.278134       0.277152          4.295350         4.190446
3  18-004A        0.276268       0.274535         17.547552        16.011004
4  18-005U        0.249210       0.248697          4.692963         3.924469

If robust and non-robust parameters differ strongly, inspect the station for outlier frequencies, poor dimensionality, or unstable strike.

Synthetic Sanity Check#

For development and training, it is useful to test the decomposition on a known distorted 2-D tensor. This example constructs a small synthetic site-like object with a known distortion matrix and checks whether the correction reduces diagonal leakage.

 1import numpy as np
 2
 3from pycsamt.emtools.gb import groom_bailey_table
 4
 5class ZBlock:
 6    def __init__(self, z, freq):
 7        self.z = z
 8        self.freq = freq
 9        self.z_err = None
10
11class Site:
12    station = "SYN001"
13
14    def __init__(self, z, freq):
15        self.Z = ZBlock(z, freq)
16
17freq = np.logspace(0, 3, 12)
18regional = np.zeros((freq.size, 2, 2), dtype=complex)
19regional[:, 0, 1] = 1.0 + 0.2j
20regional[:, 1, 0] = -0.8 + 0.1j
21
22D = np.array([[1.0, 0.25], [-0.15, 1.1]])
23observed = D[None, :, :] @ regional
24site = Site(observed, freq)
25
26table = groom_bailey_table([site], robust=False)
27
28print(table[["station", "rms_fit", "diagonal_ratio_before", "diagonal_ratio_after"]])

  station       rms_fit  diagonal_ratio_before  diagonal_ratio_after
0  SYN001  3.678917e-16               0.187225          4.388355e-17

This pattern is useful when you need to verify behavior after changing preprocessing code. Real surveys should still be assessed with their own dimensionality and strike diagnostics.

Integrate With Pre-2D Assessment#

The dimensionality guide includes pre2d_inversion_assessment. After running Groom-Bailey, record whether it was attempted and applied.

 1from pycsamt.emtools.dimensionality import pre2d_inversion_assessment
 2from pycsamt.emtools.gb import groom_bailey_decomposition
 3
 4survey = "data/AMT/WILLY_DATA/L18PLT"
 5band = (1e-3, 10.0)
 6
 7gb = groom_bailey_decomposition(
 8    survey,
 9    apply=True,
10    band=band,
11    robust=True,
12)
13
14assessment = pre2d_inversion_assessment(
15    gb.sites,
16    band=band,
17    rotation_applied=False,
18    groom_bailey_attempted=True,
19    groom_bailey_applied=gb.applied,
20    groom_bailey_reason="Applied pycsamt.emtools.gb real 2-D distortion fit.",
21)
22
23print(assessment[["station", "frac_3d", "groom_bailey_applied", "recommendation"]].head())

   station   frac_3d  groom_bailey_applied               recommendation
0  18-001A  0.974359                  True  review_3d_effects_before_2d
1  18-002U  0.974359                  True  review_3d_effects_before_2d
2  18-003A  0.974359                  True  review_3d_effects_before_2d
3  18-004A  0.974359                  True  review_3d_effects_before_2d
4  18-005U  0.923077                  True  review_3d_effects_before_2d

This makes the correction auditable in reports and manuscripts.

Reading The Results#

Use this interpretation order:

  1. Confirm that dimensionality and strike are acceptable in the selected period band.

  2. Fit groom_bailey_table without applying correction.

  3. Inspect status, n_freq, rms_fit, and diagonal ratios.

  4. Compare robust and non-robust fits if outliers are likely.

  5. Apply correction only to stations with acceptable fits.

  6. Save the table and pre-2D assessment with the inversion inputs.

Common Failure Modes#

Insufficient frequencies

The selected period band has fewer than min_freq valid tensor rows. Widen the band or skip the station.

High residual fit

The station may not be well described by a frequency-independent real distortion matrix times a 2-D regional tensor.

Diagonal ratio does not improve

Correction may not be useful for that station. Review strike, dimensionality, and the period band.

Very large twist, shear, or anisotropy

Large parameters may indicate strong galvanic distortion, but they can also indicate a poor model assumption.

Treating gain as static shift

The scalar gain ambiguity is not uniquely solved here. Use static shift workflows and independent constraints when gain matters.

Applying correction globally

Do not apply every fitted row blindly. Filter by status and quality diagnostics first.

Saving A Reproducible Bundle#

Save the fitted table, accepted subset, and pre-2D assessment.

 1from pathlib import Path
 2
 3from pycsamt.emtools.dimensionality import pre2d_inversion_assessment
 4from pycsamt.emtools.gb import apply_groom_bailey, groom_bailey_table
 5
 6survey = "data/AMT/WILLY_DATA/L18PLT"
 7band = (1e-3, 10.0)
 8out = Path("outputs/gb_l18plt")
 9out.mkdir(parents=True, exist_ok=True)
10
11table = groom_bailey_table(survey, band=band, robust=True)
12accepted = table.loc[
13    (table["status"] == "ok")
14    & (table["rms_fit"] < 0.25)
15    & (table["diagonal_ratio_after"] < table["diagonal_ratio_before"])
16].copy()
17
18corrected = apply_groom_bailey(survey, table=accepted, inplace=False)
19assessment = pre2d_inversion_assessment(
20    corrected,
21    band=band,
22    groom_bailey_attempted=True,
23    groom_bailey_applied=True,
24    groom_bailey_reason="Applied accepted Groom-Bailey station fits.",
25)
26
27table.to_csv(out / "groom_bailey_table.csv", index=False)
28accepted.to_csv(out / "groom_bailey_accepted.csv", index=False)
29assessment.to_csv(out / "pre2d_assessment_after_gb.csv", index=False)

Worked Workflow#

There is currently no dedicated Sphinx-Gallery plot_gb.py example in docs/examples/emtools. Until one is added, use the examples in this page as the worked workflow:

  • estimate the table without applying correction;

  • filter stations by fit quality;

  • apply correction to accepted stations;

  • record the result in the pre-2D assessment.