Constraining the inversion

Constraining the inversion#

This notebook explains some additional approaches we employ to smallness regularization. In this scenario, we aim to invert for the sediment-basement contact with some outcropping basement at the surface. We employ the following constraints:

Inverted basement should not cross the surface topography (upper_confining_layer)
Inverted basement should not go below 100 m elevation ((lower_confining_layer))
Limit changes to starting model in areas of known basement depth using a weighting grid.
For outcropping basement locations, exclude prisms from the inversion to reduce computation.

Again, we will use the same synthetic data from the past examples.

Import packages#

[1]:

# set EPSG for plotting functions
import os

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

import invert4geom

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

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

Gravity anomaly#

The input gravity anomaly should be corrected for topographic mass (i.e. a Bouguer anomaly). We will forward model the gravity effect of topography and remove it to get the Bouguer anomaly. We already calculated this to create the observed gravity, so we will just subtract it.

[8]:

grav_data["gravity_anomaly"] = grav_data.observed_gravity - grav_data.surface_grav

Create starting model#

Create a starting model from the interpolation of the 10 constraint points. Also use each grid cell of outcropping basement as an additional constraint point so the starting model equals the basement where it outcrops.

During interpolation, use upper_confining_layer and lower_confining_layer to make the starting model doesn’t ever go above the surface, or below a flat layer at 100 m elevation.

[9]:

# make all points where basement outcrops as constraint points
constraint_points = outcropping_basement.to_dataframe().reset_index().dropna()
constraint_points

[9]:

	northing	easting	upward
4	0.0	2000.0	617.558404
5	0.0	2500.0	612.739407
6	0.0	3000.0	607.952814
7	0.0	3500.0	603.200063
8	0.0	4000.0	598.482591
...	...	...	...
4936	30000.0	38000.0	540.007207
4937	30000.0	38500.0	540.369652
4938	30000.0	39000.0	540.812371
4939	30000.0	39500.0	541.335200
4940	30000.0	40000.0	541.937950

2870 rows × 3 columns

[10]:

# grid the sampled values using verde
starting_topography_kwargs = dict(
    method="splines",
    region=basement_model.region,
    spacing=basement_model.spacing,
    constraints_df=constraint_points,
    dampings=None,
    upper_confining_layer=true_surface_topography,
    lower_confining_layer=xr.full_like(true_surface_topography, 100),
)

starting_topography = invert4geom.create_topography(**starting_topography_kwargs)

_ = ptk.grid_compare(
    true_basement_topography,
    starting_topography.upward,
    grid1_name="True topography",
    grid2_name="Starting topography",
    robust=True,
    hist=True,
    inset=False,
    verbose="q",
    title="difference",
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)

../_images/how_to_constrain_the_inversion_14_0.png

Add a mask variable to the starting topography. Grid cells with a mask of NaN denote where basement depth is known, and these grid cells are excluded from the inversion and thus not altered. Additionally, since these prisms are excluded, the inversion will be faster due to fewer forward calculations, and a smaller Jacobian matrix.

[11]:

# distance to nearest constraint
min_dist = invert4geom.normalized_mindist(
    constraint_points,
    starting_topography.upward,
)
mask = xr.where(min_dist < 500, np.nan, 1)
mask.plot()

[11]:

<matplotlib.collections.QuadMesh at 0x7f41ac50f560>

../_images/how_to_constrain_the_inversion_16_1.png

[12]:

ds = starting_topography

# add mask to dataset
ds["mask"] = mask

[13]:

model = invert4geom.create_model(
    zref=true_basement_topography.to_numpy().mean(),
    density_contrast=basement_density_contrast,
    topography=ds,
    upper_confining_layer=true_surface_topography,
)
model.inv.plot_model(
    color_by="density",
    zscale=20,
)
model

Inversion#

[16]:

# setup the inversion
inv = invert4geom.Inversion(
    grav_data,
    model,
    solver_damping=0.01,
    # set stopping criteria
    max_iterations=20,
    l2_norm_tolerance=0.03,
    delta_l2_norm_tolerance=1.1,
)

[17]:

inv.invert(
    plot_dynamic_convergence=True,
    results_fname="../tmp/constraining_the_inversion",
)

../_images/how_to_constrain_the_inversion_24_0.png

[18]:

inv.stats_df

[18]:

	iteration	rmse	l2_norm	delta_l2_norm	iter_time_sec
0	0.0	1.217460	1.103386	inf	NaN
1	1.0	0.033469	0.182944	6.031264	2.057257
2	2.0	0.006045	0.077747	2.353073	1.680057
3	3.0	0.003238	0.056902	1.366331	1.451522
4	4.0	0.002237	0.047295	1.203141	1.441257
5	5.0	0.001716	0.041427	1.141637	1.424009
6	6.0	0.001395	0.037353	1.109065	1.459401
7	7.0	0.001177	0.034307	1.088795	1.470792
8	8.0	0.001019	0.031914	1.074967	1.458463

[19]:

inv.plot_inversion_results(
    iters_to_plot=2,
)

../_images/how_to_constrain_the_inversion_26_0.png

../_images/how_to_constrain_the_inversion_26_1.png

../_images/how_to_constrain_the_inversion_26_2.png

[20]:

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

../_images/how_to_constrain_the_inversion_27_0.png

[21]:

# sample the inverted topography at the constraint points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    inv.model.topography,
    "inverted_topography",
)

rmse_without_weighting = invert4geom.rmse(
    constraint_points.upward - constraint_points.inverted_topography
)
max_error_without_weighting = vd.maxabs(
    constraint_points.upward - constraint_points.inverted_topography
)

print(f"RMSE at constraints: {round(rmse_without_weighting, 1)} m")
print(f"max error at constraints: {round(max_error_without_weighting, 1)} m")

RMSE at constraints: 0.0 m
max error at constraints: 0.0 m

Constraining the inversion

Contents

Constraining the inversion#

Import packages#

Create observed gravity data#

Gravity anomaly#

Create starting model#

Gravity misfit#

Inversion#