Simple inversion

2. Simple inversion#

Here we will walkthrough a very simple gravity inversion using synthetic data. The goal of the inversion is to recover the geometry of a layer. In this case, the layer is simply the surface of the Earth, which is represented by the density contrast between air and rock.

2.1. Import packages#

We need to import the Invert4Geom package, as well as a few others.

[1]:

# 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.

2.4. Create a starting topography grid#

Here we assume we know nothing about the topography, so we create a flat grid with a value of 500 m. Instead, you can use an existing topography model, as long as it’s an xarray Dataset with a variable upward with the height in meters and for prism models, projected coordinates named easting and northing or for tesseroid models, geographic coordinates named longitude and latitude.

[6]:

# make a flat grid at 500 m
grid_coords = vd.grid_coordinates(
    spacing=1000,
    region=(0, 40000, 0, 30000),
)
starting_topography = vd.make_xarray_grid(
    grid_coords,
    data=np.ones_like(grid_coords[0]) * 500,
    data_names="upward",
)
starting_topography

[6]:

<xarray.Dataset> Size: 11kB
Dimensions:   (northing: 31, easting: 41)
Coordinates:
  * northing  (northing) float64 248B 0.0 1e+03 2e+03 ... 2.8e+04 2.9e+04 3e+04
  * easting   (easting) float64 328B 0.0 1e+03 2e+03 ... 3.8e+04 3.9e+04 4e+04
Data variables:
    upward    (northing, easting) float64 10kB 500.0 500.0 500.0 ... 500.0 500.0

2.7. Initialize the Inversion class#

The class Inversion holds the gravity dataset as an instance attribute data and the model dataset as an instance attribute model, as well as all the parameters of the inversion, such as the stopping criteria, and options for regularization, and how to run the inversion. Some of the options you can set during initialization include:

max_iterations: set to avoid very long running inversions, this should be set high to avoid early termination.
l2_norm_tolerance: defines when the inversion should end based on the l2-norm of the residual gravity misfit. The l2-norm is the square-root of the root mean squared misfit.
delta_l2_norm_tolerance: defines when the inversion should end if it is not making any significant progress. A value of 1 means the l2-norm has not changed from the past iteration, while a value of 1.005 means the l2-norm has decreased by ~0.5%.
solver_damping: typically between ~0 and ~1 where high values gives smoother topography.
apply_weighting_grid: a flag to use the optional weighting_grid to ensure adhere to constraints`.

During the inversion, instance attributes will be added such as the current iteration, the L2-norm and the computation time. The attributes data and model will be updated during the inversion to track the gravity and topography during each iteration.

[12]:

inv = invert4geom.Inversion(
    data,
    model,
    solver_damping=0.1,
    # set stopping criteria
    max_iterations=30,
    l2_norm_tolerance=0.45,
    delta_l2_norm_tolerance=1.005,
)

2.8. Perform inversion#

Now that we have the Inversion object which holds the starting model and residual gravity misfit data we can start the inversion. As the inversion progresses, the attributes inv.data and inv.model will be updated.

[13]:

inv.invert(
    plot_dynamic_convergence=True,
)

../_images/tutorial_02_simple_inversion_25_0.png

The above plot shows how the l2-norm and delta l2-norm have been reduced during the inversion. The red line at the bottom marks the set tolerances for each of these. This show that the inversion terminated at the 12 inversion since the l2-norm was lower than the set l2-norm tolerance.

We can access the statistics of the inversion results with the attribute stats_df.

[14]:

inv.stats_df

[14]:

	iteration	rmse	l2_norm	delta_l2_norm	iter_time_sec
0	0.0	6.743327	2.596792	inf	NaN
1	1.0	3.737865	1.933356	1.343153	0.685592
2	2.0	2.178281	1.475900	1.309950	0.124484
3	3.0	1.350640	1.162170	1.269952	0.126689
4	4.0	0.896104	0.946628	1.227695	0.127047
5	5.0	0.634186	0.796358	1.188697	0.132842
6	6.0	0.474482	0.688827	1.156108	0.142375
7	7.0	0.371599	0.609589	1.129985	0.153705
8	8.0	0.302185	0.549713	1.108922	0.156976
9	9.0	0.253624	0.503611	1.091543	0.150536
10	10.0	0.218663	0.467615	1.076979	0.161218
11	11.0	0.192876	0.439176	1.064753	0.161538

We can see why the inversion terminated with termination_reason

[15]:

inv.termination_reason

[15]:

['l2-norm tolerance']

We can see the gravity data for each iteration:

[16]:

inv.data

We can see the topography results for each iteration:

[17]:

inv.model

The final inverted topography can be accessed from the topography variable of the prism layer.

[18]:

inv.model.topography.plot()

[18]:

<matplotlib.collections.QuadMesh at 0x7fe2245631a0>

../_images/tutorial_02_simple_inversion_36_1.png

The below three plots show:

a comparison of the starting topography (flat) and the final inverted topography
a comparison of the residual gravity misfit before and after the inversion
the results of 4 of the 12 iterations where the left column shows that iterations starting residual misfit, the middle column shows the current state of the topography, and the right column shows the correction that was applied to the topography during that iteration

[19]:

model

[20]:

inv.plot_inversion_results(
    iters_to_plot=4,
)

../_images/tutorial_02_simple_inversion_39_0.png

../_images/tutorial_02_simple_inversion_39_1.png

../_images/tutorial_02_simple_inversion_39_2.png

Below we compare the true topography with the results of the inversion to see how it performed. In this example, we assumed we knew the correct density contrast value and zref, and added minimal noise to the observed gravity, so we expect the inversion to accurately recover the true topography.

[21]:

_ = 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",
    frame=True,
)

../_images/tutorial_02_simple_inversion_41_0.png

As you can see, the inversion successfully recovered the true topography. The root mean square difference between the true and recovered topography is low, but this is not too surprising since we gave the inversion the true density contrast and appropriate zref value.

See the density optimization notebook and the reference level optimization notebook for examples of how to best choose those values in scenarios where you don’t know them.

The damping cross validation notebook will explain how to use cross validation to choose the optimal damping parameter value. Here, we simply chose an arbitrary value which appeared to work well. With real data with noise, this choice of the damping parameter becomes very important.

Simple inversion

Contents

2. Simple inversion#

2.1. Import packages#

2.2. Create observed gravity data#

2.3. Initialize the gravity data#

2.4. Create a starting topography grid#

2.5. Create a Model object#

2.6. Gravity misfit#

2.7. Initialize the Inversion class#

2.8. Perform inversion#