4. Damping value cross validation

The choice of the damping value used in the inversion directly affects the results and therefore needs to be carefully chosen. We provide some tools to perform a cross-validation with a range of damping values in order to chose the optimal value.

This cross validation routine is adapted from Uieda & Barbosa 2017. See their paper for a detailed discussion.

The main idea is to re-sample the gravity data at a finer resolution and separate the data into a testing and training set. The inversion is performed with a range of damping values, each time using just the training set. Each resulting inverted topo from each of the damping values is then forward modelled onto the locations of the testing data, which were not included in the inversion. The optimal damping value is the one which results in the smallest difference between the forward modeled data and the testing data.

4.1. Import packages

# set EPSG for plotting functions
import os

import numpy as np
import polartoolkit as ptk
import verde as vd

import invert4geom

os.environ["POLARTOOLKIT_EPSG"] = "3031"

/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/UQpy/__init__.py:6: UserWarning:

pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.

4.2. Create synthetic data

Below we will create some synthetic topography data to act as the real topography, and it will forward model its gravity effect with some added noise. It will forward model it at 2x the grid spacing to create testing and training datasets to use for the damping parameter cross validation.

# the density contrast is between rock (~2670 kg/m3) and air (~1 kg/m3)
true_density_contrast = 2670 - 1
true_topography, _, _, observed_gravity = invert4geom.load_synthetic_model(
    spacing=1000,
    region=(0, 40000, 0, 30000),
    density_contrast=true_density_contrast,
    gravity_noise=0.2,
    plot_topography=True,
)
true_zref = true_topography.values.mean()
print(f"mean of true topography: {true_zref} m ")

mean of true topography: 492.2704164812973 m

4.3. Gravity misfit

As in the last notebook, for simplicity here we assume that we know the true density contrast and the reference level of the true topography, and use these values to create our starting model. Note that in a real world scenario, these would be unknowns which would need to be carefully chosen, as explained in the following notebooks.

data = invert4geom.create_data(observed_gravity)
data

<xarray.Dataset> Size: 21kB
Dimensions:          (northing: 31, easting: 41)
Coordinates:
  * northing         (northing) float64 248B 0.0 1e+03 2e+03 ... 2.9e+04 3e+04
  * easting          (easting) float64 328B 0.0 1e+03 2e+03 ... 3.9e+04 4e+04
Data variables:
    upward           (northing, easting) float64 10kB 1e+03 1e+03 ... 1e+03
    gravity_anomaly  (northing, easting) float64 10kB 9.57 10.41 ... 3.12 2.973
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       1000.0
    buffer_width:  3000.0
    inner_region:  (3000.0, 37000.0, 3000.0, 27000.0)
    dataset_type:  data
    model_type:    prisms
    coord_names:   ('easting', 'northing')

xarray.Dataset

• Dimensions:
  • northing: 31
  • easting: 41
• Coordinates: (2)
  • northing
    (northing)
    float64
    0.0 1e+03 2e+03 ... 2.9e+04 3e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000.])

  • easting
    (easting)
    float64
    0.0 1e+03 2e+03 ... 3.9e+04 4e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000., 31000., 32000., 33000., 34000., 35000.,
           36000., 37000., 38000., 39000., 40000.])

• Data variables: (2)
  • upward
    (northing, easting)
    float64
    1e+03 1e+03 1e+03 ... 1e+03 1e+03

    array([[1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           ...,
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.]], shape=(31, 41))

  • gravity_anomaly
    (northing, easting)
    float64
    9.57 10.41 10.09 ... 3.12 2.973

    array([[ 9.56972748, 10.40652495, 10.08807729, ..., -4.99960775,
            -4.59899176, -4.28599967],
           [11.11673226, 12.23445744, 11.61235841, ..., -6.19476666,
            -6.21913632, -5.24379889],
           [10.72913615, 11.77128303, 11.23716478, ..., -7.14516969,
            -6.74908688, -5.63212927],
           ...,
           [17.39173674, 19.48062958, 19.07342115, ...,  2.23956654,
             2.26864547,  2.2228519 ],
           [17.37848222, 19.9306603 , 19.07296993, ...,  3.37340862,
             3.16044899,  2.94565266],
           [15.60962256, 17.71225442, 17.39604885, ...,  3.40265911,
             3.12030258,  2.97283375]], shape=(31, 41))

• Attributes: (7)

  region :

  (0.0, 40000.0, 0.0, 30000.0)

  spacing :

  1000.0

  buffer_width :

  3000.0

  inner_region :

  (3000.0, 37000.0, 3000.0, 27000.0)

  dataset_type :

  data

  model_type :

  prisms

  coord_names :

  ('easting', 'northing')

# make a flat grid at the height of the zref
grid_coords = vd.grid_coordinates(region=data.region, spacing=data.spacing)
starting_topography = vd.make_xarray_grid(
    grid_coords, data=np.ones_like(grid_coords[0]) * true_zref, data_names="upward"
)

model = invert4geom.create_model(
    zref=true_zref,
    density_contrast=true_density_contrast,
    topography=starting_topography,
)

data.inv.forward_gravity(
    model,
    progressbar=True,
)

In many cases, we want to remove a regional signal from the misfit to isolate the residual signal. In this simple case, we assume there is no regional misfit and set it to 0. The function below will subtract the column forward_gravity from gravity_anomaly to create misfit. The reg column (chosen here to be 0) is then subtracted from misfit to get res, the residual gravity misfit.

data.inv.regional_separation(
    method="constant",
    constant=0,
)
data

<xarray.Dataset> Size: 102kB
Dimensions:                   (northing: 31, easting: 41)
Coordinates:
  * northing                  (northing) float64 248B 0.0 1e+03 ... 3e+04
  * easting                   (easting) float64 328B 0.0 1e+03 ... 3.9e+04 4e+04
Data variables:
    upward                    (northing, easting) float64 10kB 1e+03 ... 1e+03
    gravity_anomaly           (northing, easting) float64 10kB 9.57 ... 2.973
    forward_gravity           (northing, easting) float64 10kB -0.0 ... -0.0
    misfit                    (northing, easting) float64 10kB 9.57 ... 2.973
    reg                       (northing, easting) float64 10kB 0.0 0.0 ... 0.0
    res                       (northing, easting) float64 10kB 9.57 ... 2.973
    starting_forward_gravity  (northing, easting) float64 10kB -0.0 ... -0.0
    starting_misfit           (northing, easting) float64 10kB 9.57 ... 2.973
    starting_reg              (northing, easting) float64 10kB 0.0 0.0 ... 0.0
    starting_res              (northing, easting) float64 10kB 9.57 ... 2.973
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       1000.0
    buffer_width:  3000.0
    inner_region:  (3000.0, 37000.0, 3000.0, 27000.0)
    dataset_type:  data
    model_type:    prisms
    coord_names:   ('easting', 'northing')

xarray.Dataset

• Dimensions:
  • northing: 31
  • easting: 41
• Coordinates: (2)
  • northing
    (northing)
    float64
    0.0 1e+03 2e+03 ... 2.9e+04 3e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000.])

  • easting
    (easting)
    float64
    0.0 1e+03 2e+03 ... 3.9e+04 4e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000., 31000., 32000., 33000., 34000., 35000.,
           36000., 37000., 38000., 39000., 40000.])

• Data variables: (10)
  • upward
    (northing, easting)
    float64
    1e+03 1e+03 1e+03 ... 1e+03 1e+03

    array([[1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           ...,
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.]], shape=(31, 41))

  • gravity_anomaly
    (northing, easting)
    float64
    9.57 10.41 10.09 ... 3.12 2.973

    array([[ 9.56972748, 10.40652495, 10.08807729, ..., -4.99960775,
            -4.59899176, -4.28599967],
           [11.11673226, 12.23445744, 11.61235841, ..., -6.19476666,
            -6.21913632, -5.24379889],
           [10.72913615, 11.77128303, 11.23716478, ..., -7.14516969,
            -6.74908688, -5.63212927],
           ...,
           [17.39173674, 19.48062958, 19.07342115, ...,  2.23956654,
             2.26864547,  2.2228519 ],
           [17.37848222, 19.9306603 , 19.07296993, ...,  3.37340862,
             3.16044899,  2.94565266],
           [15.60962256, 17.71225442, 17.39604885, ...,  3.40265911,
             3.12030258,  2.97283375]], shape=(31, 41))

  • forward_gravity
    (northing, easting)
    float64
    -0.0 -0.0 -0.0 ... -0.0 -0.0 -0.0

    array([[-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.],
           ...,
           [-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.]], shape=(31, 41))

  • misfit
    (northing, easting)
    float64
    9.57 10.41 10.09 ... 3.12 2.973

    array([[ 9.56972748, 10.40652495, 10.08807729, ..., -4.99960775,
            -4.59899176, -4.28599967],
           [11.11673226, 12.23445744, 11.61235841, ..., -6.19476666,
            -6.21913632, -5.24379889],
           [10.72913615, 11.77128303, 11.23716478, ..., -7.14516969,
            -6.74908688, -5.63212927],
           ...,
           [17.39173674, 19.48062958, 19.07342115, ...,  2.23956654,
             2.26864547,  2.2228519 ],
           [17.37848222, 19.9306603 , 19.07296993, ...,  3.37340862,
             3.16044899,  2.94565266],
           [15.60962256, 17.71225442, 17.39604885, ...,  3.40265911,
             3.12030258,  2.97283375]], shape=(31, 41))

  • reg
    (northing, easting)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(31, 41))

  • res
    (northing, easting)
    float64
    9.57 10.41 10.09 ... 3.12 2.973

    array([[ 9.56972748, 10.40652495, 10.08807729, ..., -4.99960775,
            -4.59899176, -4.28599967],
           [11.11673226, 12.23445744, 11.61235841, ..., -6.19476666,
            -6.21913632, -5.24379889],
           [10.72913615, 11.77128303, 11.23716478, ..., -7.14516969,
            -6.74908688, -5.63212927],
           ...,
           [17.39173674, 19.48062958, 19.07342115, ...,  2.23956654,
             2.26864547,  2.2228519 ],
           [17.37848222, 19.9306603 , 19.07296993, ...,  3.37340862,
             3.16044899,  2.94565266],
           [15.60962256, 17.71225442, 17.39604885, ...,  3.40265911,
             3.12030258,  2.97283375]], shape=(31, 41))

  • starting_forward_gravity
    (northing, easting)
    float64
    -0.0 -0.0 -0.0 ... -0.0 -0.0 -0.0

    array([[-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.],
           ...,
           [-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.],
           [-0., -0., -0., ..., -0., -0., -0.]], shape=(31, 41))

  • starting_misfit
    (northing, easting)
    float64
    9.57 10.41 10.09 ... 3.12 2.973

    array([[ 9.56972748, 10.40652495, 10.08807729, ..., -4.99960775,
            -4.59899176, -4.28599967],
           [11.11673226, 12.23445744, 11.61235841, ..., -6.19476666,
            -6.21913632, -5.24379889],
           [10.72913615, 11.77128303, 11.23716478, ..., -7.14516969,
            -6.74908688, -5.63212927],
           ...,
           [17.39173674, 19.48062958, 19.07342115, ...,  2.23956654,
             2.26864547,  2.2228519 ],
           [17.37848222, 19.9306603 , 19.07296993, ...,  3.37340862,
             3.16044899,  2.94565266],
           [15.60962256, 17.71225442, 17.39604885, ...,  3.40265911,
             3.12030258,  2.97283375]], shape=(31, 41))

  • starting_reg
    (northing, easting)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(31, 41))

  • starting_res
    (northing, easting)
    float64
    9.57 10.41 10.09 ... 3.12 2.973

    array([[ 9.56972748, 10.40652495, 10.08807729, ..., -4.99960775,
            -4.59899176, -4.28599967],
           [11.11673226, 12.23445744, 11.61235841, ..., -6.19476666,
            -6.21913632, -5.24379889],
           [10.72913615, 11.77128303, 11.23716478, ..., -7.14516969,
            -6.74908688, -5.63212927],
           ...,
           [17.39173674, 19.48062958, 19.07342115, ...,  2.23956654,
             2.26864547,  2.2228519 ],
           [17.37848222, 19.9306603 , 19.07296993, ...,  3.37340862,
             3.16044899,  2.94565266],
           [15.60962256, 17.71225442, 17.39604885, ...,  3.40265911,
             3.12030258,  2.97283375]], shape=(31, 41))

• Attributes: (7)

  region :

  (0.0, 40000.0, 0.0, 30000.0)

  spacing :

  1000.0

  buffer_width :

  3000.0

  inner_region :

  (3000.0, 37000.0, 3000.0, 27000.0)

  dataset_type :

  data

  model_type :

  prisms

  coord_names :

  ('easting', 'northing')

data.inv.plot_anomalies()

makecpt [ERROR]: Option T: min >= max

4.4. Get Cross Validation Score

We will perform an inversion with a single damping value and calculate a score for it. The grid shows how to score is calculated, as the RMS difference between the predicted and observed gravity data.

# setup the inversion
inv = invert4geom.Inversion(
    data,
    model,
    solver_damping=0.1,
    # set stopping criteria
    max_iterations=30,
    l2_norm_tolerance=0.2**2,  # gravity error is .2 mGal or L2-norm of ~0.45
    delta_l2_norm_tolerance=1.01,
)

# resample data at 1/2 spacing to include test points for cross-validation
inv.data = invert4geom.add_test_points(inv.data)
inv.data

<xarray.Dataset> Size: 401kB
Dimensions:                   (northing: 61, easting: 81)
Coordinates:
  * northing                  (northing) float64 488B 0.0 500.0 ... 3e+04
  * easting                   (easting) float64 648B 0.0 500.0 ... 4e+04
Data variables:
    test                      (northing, easting) bool 5kB False True ... False
    upward                    (northing, easting) float64 40kB 1e+03 ... 1e+03
    gravity_anomaly           (northing, easting) float64 40kB 9.57 ... 2.973
    forward_gravity           (northing, easting) float64 40kB 0.0 0.0 ... 0.0
    misfit                    (northing, easting) float64 40kB 9.57 ... 2.973
    reg                       (northing, easting) float64 40kB 0.0 0.0 ... 0.0
    res                       (northing, easting) float64 40kB 9.57 ... 2.973
    starting_forward_gravity  (northing, easting) float64 40kB 0.0 0.0 ... 0.0
    starting_misfit           (northing, easting) float64 40kB 9.57 ... 2.973
    starting_reg              (northing, easting) float64 40kB 0.0 0.0 ... 0.0
    starting_res              (northing, easting) float64 40kB 9.57 ... 2.973
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       1000.0
    buffer_width:  3000.0
    inner_region:  (3000.0, 37000.0, 3000.0, 27000.0)
    dataset_type:  data
    model_type:    prisms
    coord_names:   ('easting', 'northing')

xarray.Dataset

• Dimensions:
  • northing: 61
  • easting: 81
• Coordinates: (2)
  • northing
    (northing)
    float64
    0.0 500.0 1e+03 ... 2.95e+04 3e+04

    array([    0.,   500.,  1000.,  1500.,  2000.,  2500.,  3000.,  3500.,  4000.,
            4500.,  5000.,  5500.,  6000.,  6500.,  7000.,  7500.,  8000.,  8500.,
            9000.,  9500., 10000., 10500., 11000., 11500., 12000., 12500., 13000.,
           13500., 14000., 14500., 15000., 15500., 16000., 16500., 17000., 17500.,
           18000., 18500., 19000., 19500., 20000., 20500., 21000., 21500., 22000.,
           22500., 23000., 23500., 24000., 24500., 25000., 25500., 26000., 26500.,
           27000., 27500., 28000., 28500., 29000., 29500., 30000.])

  • easting
    (easting)
    float64
    0.0 500.0 1e+03 ... 3.95e+04 4e+04

    array([    0.,   500.,  1000.,  1500.,  2000.,  2500.,  3000.,  3500.,  4000.,
            4500.,  5000.,  5500.,  6000.,  6500.,  7000.,  7500.,  8000.,  8500.,
            9000.,  9500., 10000., 10500., 11000., 11500., 12000., 12500., 13000.,
           13500., 14000., 14500., 15000., 15500., 16000., 16500., 17000., 17500.,
           18000., 18500., 19000., 19500., 20000., 20500., 21000., 21500., 22000.,
           22500., 23000., 23500., 24000., 24500., 25000., 25500., 26000., 26500.,
           27000., 27500., 28000., 28500., 29000., 29500., 30000., 30500., 31000.,
           31500., 32000., 32500., 33000., 33500., 34000., 34500., 35000., 35500.,
           36000., 36500., 37000., 37500., 38000., 38500., 39000., 39500., 40000.])

• Data variables: (11)
  • test
    (northing, easting)
    bool
    False True False ... True False

    array([[False,  True, False, ..., False,  True, False],
           [ True,  True,  True, ...,  True,  True,  True],
           [False,  True, False, ..., False,  True, False],
           ...,
           [False,  True, False, ..., False,  True, False],
           [ True,  True,  True, ...,  True,  True,  True],
           [False,  True, False, ..., False,  True, False]], shape=(61, 81))

  • upward
    (northing, easting)
    float64
    1e+03 1e+03 1e+03 ... 1e+03 1e+03

    array([[1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           ...,
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.],
           [1000., 1000., 1000., ..., 1000., 1000., 1000.]], shape=(61, 81))

  • gravity_anomaly
    (northing, easting)
    float64
    9.57 10.06 10.41 ... 3.038 2.973

    array([[ 9.5697279 , 10.06032902, 10.40652466, ..., -4.59899187,
            -4.43701935, -4.28599977],
           [10.4641428 , 10.94956768, 11.39148229, ..., -5.48267764,
            -5.15328769, -4.80049103],
           [11.1167326 , 11.66342121, 12.23445702, ..., -6.21913624,
            -5.75835735, -5.24379873],
           ...,
           [17.37848282, 18.75796282, 19.93066025, ...,  3.16044903,
             3.00968942,  2.94565272],
           [16.60377824, 17.84360372, 18.83705711, ...,  3.2070532 ,
             3.07967874,  3.00271942],
           [15.60962296, 16.81211543, 17.71225357, ...,  3.12030268,
             3.03813769,  2.97283363]], shape=(61, 81))

  • forward_gravity
    (northing, easting)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(61, 81))

  • misfit
    (northing, easting)
    float64
    9.57 10.06 10.41 ... 3.038 2.973

    array([[ 9.5697279 , 10.06032902, 10.40652466, ..., -4.59899187,
            -4.43701935, -4.28599977],
           [10.4641428 , 10.94956768, 11.39148229, ..., -5.48267764,
            -5.15328769, -4.80049103],
           [11.1167326 , 11.66342121, 12.23445702, ..., -6.21913624,
            -5.75835735, -5.24379873],
           ...,
           [17.37848282, 18.75796282, 19.93066025, ...,  3.16044903,
             3.00968942,  2.94565272],
           [16.60377824, 17.84360372, 18.83705711, ...,  3.2070532 ,
             3.07967874,  3.00271942],
           [15.60962296, 16.81211543, 17.71225357, ...,  3.12030268,
             3.03813769,  2.97283363]], shape=(61, 81))

  • reg
    (northing, easting)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(61, 81))

  • res
    (northing, easting)
    float64
    9.57 10.06 10.41 ... 3.038 2.973

    array([[ 9.5697279 , 10.06032902, 10.40652466, ..., -4.59899187,
            -4.43701935, -4.28599977],
           [10.4641428 , 10.94956768, 11.39148229, ..., -5.48267764,
            -5.15328769, -4.80049103],
           [11.1167326 , 11.66342121, 12.23445702, ..., -6.21913624,
            -5.75835735, -5.24379873],
           ...,
           [17.37848282, 18.75796282, 19.93066025, ...,  3.16044903,
             3.00968942,  2.94565272],
           [16.60377824, 17.84360372, 18.83705711, ...,  3.2070532 ,
             3.07967874,  3.00271942],
           [15.60962296, 16.81211543, 17.71225357, ...,  3.12030268,
             3.03813769,  2.97283363]], shape=(61, 81))

  • starting_forward_gravity
    (northing, easting)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(61, 81))

  • starting_misfit
    (northing, easting)
    float64
    9.57 10.06 10.41 ... 3.038 2.973

    array([[ 9.5697279 , 10.06032902, 10.40652466, ..., -4.59899187,
            -4.43701935, -4.28599977],
           [10.4641428 , 10.94956768, 11.39148229, ..., -5.48267764,
            -5.15328769, -4.80049103],
           [11.1167326 , 11.66342121, 12.23445702, ..., -6.21913624,
            -5.75835735, -5.24379873],
           ...,
           [17.37848282, 18.75796282, 19.93066025, ...,  3.16044903,
             3.00968942,  2.94565272],
           [16.60377824, 17.84360372, 18.83705711, ...,  3.2070532 ,
             3.07967874,  3.00271942],
           [15.60962296, 16.81211543, 17.71225357, ...,  3.12030268,
             3.03813769,  2.97283363]], shape=(61, 81))

  • starting_reg
    (northing, easting)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(61, 81))

  • starting_res
    (northing, easting)
    float64
    9.57 10.06 10.41 ... 3.038 2.973

    array([[ 9.5697279 , 10.06032902, 10.40652466, ..., -4.59899187,
            -4.43701935, -4.28599977],
           [10.4641428 , 10.94956768, 11.39148229, ..., -5.48267764,
            -5.15328769, -4.80049103],
           [11.1167326 , 11.66342121, 12.23445702, ..., -6.21913624,
            -5.75835735, -5.24379873],
           ...,
           [17.37848282, 18.75796282, 19.93066025, ...,  3.16044903,
             3.00968942,  2.94565272],
           [16.60377824, 17.84360372, 18.83705711, ...,  3.2070532 ,
             3.07967874,  3.00271942],
           [15.60962296, 16.81211543, 17.71225357, ...,  3.12030268,
             3.03813769,  2.97283363]], shape=(61, 81))

• Attributes: (7)

  region :

  (0.0, 40000.0, 0.0, 30000.0)

  spacing :

  1000.0

  buffer_width :

  3000.0

  inner_region :

  (3000.0, 37000.0, 3000.0, 27000.0)

  dataset_type :

  data

  model_type :

  prisms

  coord_names :

  ('easting', 'northing')

# run inversion, calculate the score, and plot the predicted and observed gravity for
# the testing dataset
grav_cv_object = inv.gravity_score(
    plot=True,
)
inv.gravity_best_score

Inversion terminated due to max_iterations limit. Consider increasing this limit.

np.float64(0.05102558192068222)

4.5. Cross Validation

Now we can repeat this with a range of damping parameters to find the optimal (lowest) score.

This is a classic optimization problem, where we vary the damping parameter to find which value produces the best score. For this we use the Python package Optuna.

For each trial a .pickle file containing the inversion results will be created. Once all trials are complete, all but the best trial’s pickle files will be deleted. Additionally, another .pickle file will be created containing the Optuna study object with all the trial information and results.

By default, optimize_inversion_damping will use a sampler to choose damping values between the provided damping_limits and repeat this n_trials times. This may skip large portions of the parameter space. If instead you would like to systematically test damping values at a constant interval, you can set grid_search to True, and this will test all damping values as defined by:

np.logspace(<<exponent of lower damping limit>>, <<exponent of upper damping limit>>, n_trials)

damping_cv_obj = inv.optimize_inversion_damping(
    damping_limits=(0.001, 1),
    n_trials=8,
    grid_search=True,
    plot_scores=False,
    fname="../tmp/04_damping_CV",
)
inv.solver_damping

Inversion terminated due to max_iterations limit. Consider increasing this limit.

np.float64(0.019306977288832496)

damping_cv_obj.study.trials_dataframe()

	number	value	datetime_start	datetime_complete	duration	params_damping	user_attrs_fname	system_attrs_grid_id	system_attrs_search_space	state
0	0	1.730399	2026-02-17 02:17:36.651284	2026-02-17 02:17:43.101102	0 days 00:00:06.449818	0.372759	../tmp/04_damping_CV_trial_0	0	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
1	1	0.119877	2026-02-17 02:17:43.102307	2026-02-17 02:17:44.636571	0 days 00:00:01.534264	0.007197	../tmp/04_damping_CV_trial_1	1	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
2	2	0.809352	2026-02-17 02:17:44.637739	2026-02-17 02:17:45.716687	0 days 00:00:01.078948	0.002683	../tmp/04_damping_CV_trial_2	2	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
3	3	6.319568	2026-02-17 02:17:45.717765	2026-02-17 02:17:46.259139	0 days 00:00:00.541374	1.000000	../tmp/04_damping_CV_trial_3	3	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
4	4	0.023230	2026-02-17 02:17:46.260044	2026-02-17 02:17:52.302433	0 days 00:00:06.042389	0.019307	../tmp/04_damping_CV_trial_4	4	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
5	5	0.886706	2026-02-17 02:17:52.303473	2026-02-17 02:17:53.242700	0 days 00:00:00.939227	0.001000	../tmp/04_damping_CV_trial_5	5	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
6	6	0.097773	2026-02-17 02:17:53.243739	2026-02-17 02:17:59.959784	0 days 00:00:06.716045	0.138950	../tmp/04_damping_CV_trial_6	6	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE
7	7	0.026461	2026-02-17 02:17:59.960753	2026-02-17 02:18:06.618580	0 days 00:00:06.657827	0.051795	../tmp/04_damping_CV_trial_7	7	{'damping': [0.001, 0.0026826957952797246, 0.0...	COMPLETE

fig = invert4geom.plot_scores(
    damping_cv_obj.study.trials_dataframe().value.values,
    damping_cv_obj.study.trials_dataframe().params_damping.values,
    param_name="Damping",
    logx=True,
)

4.6. Plot inversion results using optimal value

inv.plot_convergence()

inv.plot_inversion_results(iters_to_plot=2)

_ = ptk.grid_compare(
    true_topography,
    inv.model.topography,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    title="difference",
    reverse_cpt=True,
    cmap="rain",
)

4.7. Run inversion with poor choice of damping value

# setup the inversion
inv = invert4geom.Inversion(
    data,
    model,
    solver_damping=0.005,
    # set stopping criteria
    l2_norm_tolerance=0.2**2,  # gravity error is .2 mGal or L2-norm of ~0.45
    delta_l2_norm_tolerance=1.01,
)

inv.invert()

_ = ptk.grid_compare(
    true_topography,
    inv.model.topography,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    title="difference",
    reverse_cpt=True,
    cmap="rain",
    # diff_lims=(-20, 20),
)

Comparing this inversion with a damping value of 0.005 to the inversion with the optimal damping value of ~0.1 shows that a low damping value allows too much noise through to the final inverted topography.

In the past example (simple inversion), we arbitrarily chose a damping value. Here, we show how by finding the optimal damping value using a cross-validation, we have improved the inversion’s performance.

The next notebook, density optimization, will show how to perform a similar procedure for choosing a density contrast value.