Including a starting model

3. Including a starting model#

Here we will expand the simple_inversion.ipynb example by showing how to incorporate a non-flat starting model. A typical scenario for where this is useful is if you have a few point measurements of the elevations of the surface you are aiming to recover. These point measurements, referred to here as constraints, may be boreholes, acoustic basement from seismic surveys, or other types of measurements.

3.1. Import packages#

[1]:

from __future__ import annotations

%load_ext autoreload
%autoreload 2

import logging

import xarray as xr
from polartoolkit import utils as polar_utils

from invert4geom import inversion, plotting, regional, synthetic, utils

# set up logging to see what's going on
logging.basicConfig(level=logging.INFO)

3.2. Get a synthetic model#

Here we will load synthetic topography data, calculate its forward gravity with added noise, and sample the topography at 10 random points which represent places we have measurements of the topography (from boreholes for example). We will then interpolate these 10 points to make a starting model of topography.

[2]:

spacing = 1000
region = (0, 40000, 0, 30000)

# the density contrast is between rock (~2670 kg/m3) and air (~1 kg/m3)
density_contrast = 2670 - 1

(
    true_topography,
    starting_topography,
    constraint_points,
    grav_df,
) = synthetic.load_synthetic_model(
    spacing=spacing,
    region=region,
    number_of_constraints=10,
    density_contrast=density_contrast,
)

INFO:invert4geom:RMSE at the constraints between the starting and true topography: 19.262348539866753 m
../_images/tutorial_03_including_starting_model_4_1.png 
INFO:invert4geom:Standard deviation used for noise: [0.2]
../_images/tutorial_03_including_starting_model_4_3.png 
[3]:

constraint_points

[3]:
easting northing upward starting_topography
0 3052.331575 20376.899884 619.779099 580.163740
1 31196.751690 24112.171083 479.563412 481.652035
2 17536.369258 11428.233994 465.786893 458.230368
3 28938.607113 1978.090407 426.413502 438.512760
4 39119.580480 8644.367979 428.773305 425.210470
5 21539.834816 27287.805832 546.901160 547.287244
6 20044.818546 6401.560607 450.406755 455.069067
7 2882.045334 13563.718855 460.979806 499.138482
8 10737.559204 27936.180591 614.084318 624.948915
9 19995.300033 746.976827 470.409607 452.397262 
[4]:

grav_df

[4]:
northing easting upward gravity_anomaly
0 0.0 0.0 1000.0 9.569727
1 0.0 1000.0 1000.0 10.406525
2 0.0 2000.0 1000.0 10.088077
3 0.0 3000.0 1000.0 9.300146
4 0.0 4000.0 1000.0 8.434769
... ... ... ... ...
1266 30000.0 36000.0 1000.0 3.117968
1267 30000.0 37000.0 1000.0 3.472211
1268 30000.0 38000.0 1000.0 3.402659
1269 30000.0 39000.0 1000.0 3.120303
1270 30000.0 40000.0 1000.0 2.972834

1271 rows Ă— 4 columns

3.3. Gravity misfit#

Now we need to calculate the forward gravity of the starting topography. We then can subtract it from our observed gravity to get a starting gravity misfit.

[5]:

# set the reference level to the mean of the constraints
zref = constraint_points.upward.mean()

# prisms above zref have positive density contrast and prisms below zref have negative
# density contrast
density_grid = xr.where(
    starting_topography >= zref, density_contrast, -density_contrast
)

# create layer of prisms
starting_prisms = utils.grids_to_prisms(
    starting_topography,
    reference=zref,
    density=density_grid,
)

plotting.show_prism_layers(
    starting_prisms,
    color_by="thickness",
    log_scale=False,
    zscale=20,
    backend="static",
)

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen
../_images/tutorial_03_including_starting_model_8_1.png 
[6]:

# calculate forward gravity of starting prism layer
grav_df["starting_gravity"] = starting_prisms.prism_layer.gravity(
    coordinates=(
        grav_df.easting,
        grav_df.northing,
        grav_df.upward,
    ),
    field="g_z",
    progressbar=True,
)

# estimate regional with the mean misfit at constraints
grav_df = regional.regional_constant(
    grav_df=grav_df,
    constraints_df=constraint_points,
)

grav_df

INFO:invert4geom:using median gravity misfit of constraint points for regional field: 0.05486141840371294 mGal

[6]:
northing easting upward gravity_anomaly starting_gravity misfit reg res
0 0.0 0.0 1000.0 9.569727 -4.616509 14.186236 0.054861 14.131375
1 0.0 1000.0 1000.0 10.406525 -5.479088 15.885613 0.054861 15.830752
2 0.0 2000.0 1000.0 10.088077 -5.617448 15.705525 0.054861 15.650664
3 0.0 3000.0 1000.0 9.300146 -5.622104 14.922250 0.054861 14.867389
4 0.0 4000.0 1000.0 8.434769 -5.578906 14.013675 0.054861 13.958814
... ... ... ... ... ... ... ... ...
1266 30000.0 36000.0 1000.0 3.117968 -0.899244 4.017212 0.054861 3.962351
1267 30000.0 37000.0 1000.0 3.472211 -1.230122 4.702333 0.054861 4.647472
1268 30000.0 38000.0 1000.0 3.402659 -1.539613 4.942273 0.054861 4.887411
1269 30000.0 39000.0 1000.0 3.120303 -1.788878 4.909181 0.054861 4.854320
1270 30000.0 40000.0 1000.0 2.972834 -1.716132 4.688966 0.054861 4.634104

1271 rows Ă— 8 columns

3.4. Perform inversion#

Now that we have a starting model and residual gravity misfit data we can start the inversion.

[7]:

# set Python's logging level to suppress information about the inversion\s progress
logging.getLogger().setLevel(logging.WARN)

# run the inversion
results = inversion.run_inversion(
    grav_df=grav_df,
    prism_layer=starting_prisms,
    # display the convergence of the inversion
    plot_dynamic_convergence=True,
    solver_damping=0.06,
    # set stopping criteria
    max_iterations=30,
    l2_norm_tolerance=0.45,
    delta_l2_norm_tolerance=1.005,
)

# collect the results
topo_results, grav_results, parameters, elapsed_time = results
../_images/tutorial_03_including_starting_model_11_0.png 
[8]:

# collect the results
topo_results, grav_results, parameters, elapsed_time = results

plotting.plot_inversion_results(
    grav_results,
    topo_results,
    parameters,
    region,
    iters_to_plot=2,
    plot_iter_results=True,
    plot_topo_results=True,
    plot_grav_results=True,
    constraints_df=constraint_points,
)

final_topography = topo_results.set_index(["northing", "easting"]).to_xarray().topo

_ = polar_utils.grd_compare(
    true_topography,
    final_topography,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    title="difference",
    grounding_line=False,
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)
../_images/tutorial_03_including_starting_model_12_0.png
../_images/tutorial_03_including_starting_model_12_1.png
../_images/tutorial_03_including_starting_model_12_2.png
../_images/tutorial_03_including_starting_model_12_3.png 
[9]:

# sample the inverted topography at the constraint points
constraint_points = utils.sample_grids(
    constraint_points,
    final_topography,
    "inverted_topography",
)

rmse = utils.rmse(constraint_points.upward - constraint_points.inverted_topography)
print(f"RMSE: {rmse:.2f} m")

RMSE: 2.90 m

[10]:

constraint_points

[10]:
easting northing upward starting_topography inverted_topography
0 3052.331575 20376.899884 619.779099 580.163740 621.896282
1 31196.751690 24112.171083 479.563412 481.652035 481.298143
2 17536.369258 11428.233994 465.786893 458.230368 468.311972
3 28938.607113 1978.090407 426.413502 438.512760 430.035191
4 39119.580480 8644.367979 428.773305 425.210470 432.648317
5 21539.834816 27287.805832 546.901160 547.287244 548.322451
6 20044.818546 6401.560607 450.406755 455.069067 451.832598
7 2882.045334 13563.718855 460.979806 499.138482 459.139836
8 10737.559204 27936.180591 614.084318 624.948915 619.697290
9 19995.300033 746.976827 470.409607 452.397262 472.253349

The RMSE between the constraint’s true values and the inverted topography at the constraint’s is not 0. This shows that while the starting model helped the inversion, the actual values of the constraints is not adhered too. The next inversion (adhering_to_constraints.ipynb) will show how to help the model stick to the constraints.

Contents