Adhering to point measurements

Adhering to point measurements#

This notebook explains one approach to smallness regularization, with the goal of adhering to prior knowledge of topography we’re aiming to recover. In the scenario of a sediment-basement contact inversion, if we know the basement depth at a few points from drill holes or seismic data, we want the inverted results to adhere to these points, which we refer to as constraints.

We offer 2 methods to adhere to theses points:

  1. Using the Constraint Point Minimization technique for regional misfit estimation.

    • this defines the regional field so that the residual field is 0 and constraints, resulting in minimal changing of the topography at constraints during the inversion.

  2. Using a weighting grid, which smoothly tapers each iterations topography correction to be 0 at constraint points.

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

Import packages#

[1]:

# set EPSG for plotting functions
import os

import pandas as pd
import polartoolkit as ptk
import verde as vd

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.

Create observed gravity data#

To run the inversion, we need to have observed gravity data. In this simple example, we will first create a synthetic topography, which represents the true Earth topography which we hope to recover during the inversion. From this topography, we will create a layer of vertical right-rectangular prisms, which allows us to calculated the gravity effect of the topography. This will act as our observed gravity data.

Sample the starting topography at 10 random locations and regrid with those sampled values. This simulates only knowing the depth to this topography at 10 boreholes.

[2]:

true_topography, _, _, _ = invert4geom.load_synthetic_model(
    spacing=500,
    region=(0, 40000, 0, 30000),
)

# create lower synthetic topography data
lower_topography = invert4geom.synthetic_topography_regional(
    500,
    (0, 40000, 0, 30000),
    scale=1,
    yoffset=-1000,
)

[3]:

upper_model = invert4geom.create_model(
    zref=true_topography.to_numpy().mean(),
    density_contrast=2670 - 1,
    topography=true_topography.to_dataset(name="upward"),
)
lower_model = invert4geom.create_model(
    zref=lower_topography.to_numpy().mean(),
    density_contrast=3100 - 2670,
    topography=lower_topography.to_dataset(name="upward"),
)

invert4geom.plot_prism_layers(
    [upper_model, lower_model],
    color_by="thickness",
    zscale=20,
)
../_images/how_to_adhering_to_point_measurements_5_0.png 
[4]:

# make pandas dataframe of locations to calculate gravity
# this represents the station locations of a gravity survey
# create lists of coordinates
coords = vd.grid_coordinates(
    spacing=500,
    region=(0, 40000, 0, 30000),
    pixel_register=False,
    extra_coords=1001,  # survey elevation
)

# grid the coordinates
observations = vd.make_xarray_grid(
    (coords[0], coords[1]),
    data=coords[2],
    data_names="upward",
    dims=("northing", "easting"),
)

grav_data = invert4geom.create_data(observations)

# forward gravity of upper and lower prisms
grav_data.inv.forward_gravity(upper_model, "upper_surface_grav")
grav_data.inv.forward_gravity(lower_model, "lower_surface_grav")

grav_data["gravity_anomaly"] = (
    grav_data.upper_surface_grav + grav_data.lower_surface_grav
)
grav_data

[4]:

<xarray.Dataset> Size: 159kB
Dimensions:             (northing: 61, easting: 81)
Coordinates:
  * northing            (northing) float64 488B 0.0 500.0 ... 2.95e+04 3e+04
  * easting             (easting) float64 648B 0.0 500.0 ... 3.95e+04 4e+04
Data variables:
    upward              (northing, easting) float64 40kB 1.001e+03 ... 1.001e+03
    upper_surface_grav  (northing, easting) float64 40kB 6.966 8.536 ... 2.06
    lower_surface_grav  (northing, easting) float64 40kB -0.04507 ... 0.385
    gravity_anomaly     (northing, easting) float64 40kB 6.921 8.486 ... 2.445
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       500.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')
[5]:

fig = ptk.plot_grid(
    grav_data.lower_surface_grav,
    fig_height=10,
    title="lower prisms gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSWe", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.upper_surface_grav,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="upper prisms gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.gravity_anomaly,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Observed gravity",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig.show()
../_images/how_to_adhering_to_point_measurements_7_0.png

Create “a-priori” topography measurements#

[6]:

# create 25 random point within the outcropping basement region
num_constraints = 25
coords = vd.scatter_points(
    region=vd.pad_region(grav_data.region, -500),
    size=num_constraints,
    random_state=0,
)
constraint_points = pd.DataFrame(data={"easting": coords[0], "northing": coords[1]})

# sample true topography at these points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    true_topography,
    "upward",
)
constraint_points.head()

[6]:
easting northing upward
0 21903.726653 19057.709618 472.843511
1 28392.385289 4657.245335 417.654600
2 24007.771667 27895.398594 537.989827
3 21750.444137 15633.601331 485.088965
4 17022.537174 12525.196260 471.304408

Create starting model#

[7]:

# grid the sampled values using verde
starting_topography_kwargs = dict(
    method="splines",
    region=grav_data.region,
    spacing=grav_data.spacing,
    constraints_df=constraint_points,
    # dampings=np.logspace(-40, 0, 100),
    dampings=None,
)

starting_topography = invert4geom.create_topography(**starting_topography_kwargs)

_ = ptk.grid_compare(
    true_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_adhering_to_point_measurements_11_0.png 
[8]:

# sample the inverted topography at the constraint points
constraint_points = invert4geom.sample_grids(
    constraint_points,
    starting_topography.upward,
    "starting_topography",
)

rmse_starting = invert4geom.rmse(
    constraint_points.upward - constraint_points.starting_topography
)
max_error_starting = vd.maxabs(
    constraint_points.upward - constraint_points.starting_topography
)

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

RMSE at constraints: 0.1 m
max error at constraints: 0.5 m

[9]:

model = invert4geom.create_model(
    zref=true_topography.to_numpy().mean(),
    density_contrast=2670 - 1,
    topography=starting_topography,
)
model.inv.plot_model(
    color_by="density",
    zscale=20,
)
../_images/how_to_adhering_to_point_measurements_13_0.png

Method 1: Weighting grid#

To force the inversion to adhere to the starting model we need to supply a weighting grid. At each iteration, the correction grid is multiplied by this weighting grid to alter the iteration’s correction. Therefore, this weighting grid should be ~0 at the constraints, so that they aren’t altered from the starting model. These values should increase to ~1 at a distance to allow the inversion to be un-affected at locations far from constraints.

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. For this, the regional component will just be constant, the mean of the misfit at the constraint points.

[10]:

grav_data.inv.forward_gravity(model)

grav_data.inv.regional_separation(
    method="constant",
    constraints_df=constraint_points,
)
grav_data

[10]:

<xarray.Dataset> Size: 475kB
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:
    upward                    (northing, easting) float64 40kB 1.001e+03 ... ...
    upper_surface_grav        (northing, easting) float64 40kB 6.966 ... 2.06
    lower_surface_grav        (northing, easting) float64 40kB -0.04507 ... 0...
    gravity_anomaly           (northing, easting) float64 40kB 6.921 ... 2.445
    forward_gravity           (northing, easting) float64 40kB -2.37 ... 1.78
    misfit                    (northing, easting) float64 40kB 9.291 ... 0.6649
    reg                       (northing, easting) float64 40kB -0.617 ... -0.617
    res                       (northing, easting) float64 40kB 9.908 ... 1.282
    starting_forward_gravity  (northing, easting) float64 40kB -2.37 ... 1.78
    starting_misfit           (northing, easting) float64 40kB 9.291 ... 0.6649
    starting_reg              (northing, easting) float64 40kB -0.617 ... -0.617
    starting_res              (northing, easting) float64 40kB 9.908 ... 1.282
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       500.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')
[11]:

fig = ptk.plot_grid(
    grav_data.gravity_anomaly,
    fig_height=10,
    title="Topo free disturbance",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSWe", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.misfit,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Misfit",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.reg,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Regional misfit",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.res,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Residual misfit",
    cmap="balance+h0",
    cpt_lims=[-vd.maxabs(grav_data.res), vd.maxabs(grav_data.res)],
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
    points=constraint_points,
    points_style="x.15c",
)
fig.show()

makecpt [ERROR]: Option T: min >= max
supplied min value is greater or equal to max value
../_images/how_to_adhering_to_point_measurements_17_1.png

Create weighting grid#

[12]:

# calculate the distance between each grid cell and the nearest constraint, then
# normalize those values between 0 and 1
weighting_grid = invert4geom.normalized_mindist(
    constraint_points,
    model.starting_topography,
    low=0,
    high=1,
)
weighting_grid.plot()

[12]:

<matplotlib.collections.QuadMesh at 0x7f5ad4fd8d10>
../_images/how_to_adhering_to_point_measurements_19_1.png

Perform inversion#

Now we can perform the inversion, supplying the argument apply_weighting_grid=True and supplying the grid to argument weighting_grid.

[13]:

# setup the inversion
inv = invert4geom.Inversion(
    grav_data,
    model,
    solver_damping=0.05,
    apply_weighting_grid=True,
    weighting_grid=weighting_grid,
    # set stopping criteria
    max_iterations=300,
    l2_norm_tolerance=0.3,
    delta_l2_norm_tolerance=1.005,
)

[14]:

inv.invert(
    plot_dynamic_convergence=True,
    results_fname="../tmp/weighted_inversion",
)
../_images/how_to_adhering_to_point_measurements_22_0.png 
[15]:

inv.stats_df

[15]:
iteration rmse l2_norm delta_l2_norm iter_time_sec
0 0.0 5.796239 2.407538 inf NaN
1 1.0 3.921412 1.980256 1.215771 2.143775
2 2.0 2.819363 1.679096 1.179358 1.585192
3 3.0 2.123618 1.457264 1.152225 1.663122
4 4.0 1.658120 1.287680 1.131697 1.584170
5 5.0 1.331713 1.153999 1.115842 1.667037
6 6.0 1.093920 1.045906 1.103348 1.648916
7 7.0 0.915157 0.956638 1.093314 1.607757
8 8.0 0.777230 0.881606 1.085108 1.631489
9 9.0 0.668475 0.817603 1.078281 1.646033
10 10.0 0.581142 0.762327 1.072510 1.682685
11 11.0 0.509918 0.714085 1.067557 1.704136
12 12.0 0.451059 0.671609 1.063245 1.785189
13 13.0 0.401863 0.633927 1.059443 1.707829
14 14.0 0.360338 0.600281 1.056050 1.686785
15 15.0 0.324982 0.570072 1.052992 1.760374
16 16.0 0.294650 0.542817 1.050212 1.690624
17 17.0 0.268448 0.518120 1.047666 1.752501
18 18.0 0.245676 0.495657 1.045320 1.759439
19 19.0 0.225772 0.475155 1.043149 1.707620
20 20.0 0.208286 0.456383 1.041131 1.695318
21 21.0 0.192850 0.439147 1.039249 1.758697
22 22.0 0.179165 0.423279 1.037489 1.775427
23 23.0 0.166981 0.408633 1.035841 1.749154
24 24.0 0.156091 0.395084 1.034295 1.775392
25 25.0 0.146322 0.382521 1.032842 1.739660
26 26.0 0.137529 0.370849 1.031475 1.773850
27 27.0 0.129587 0.359982 1.030187 1.763101
28 28.0 0.122392 0.349845 1.028974 1.771846
29 29.0 0.115854 0.340373 1.027830 1.825444
30 30.0 0.109895 0.331504 1.026751 1.745752
31 31.0 0.104450 0.323188 1.025733 1.763574
32 32.0 0.099462 0.315376 1.024771 1.752474
33 33.0 0.094880 0.308026 1.023862 1.793197
34 34.0 0.090661 0.301099 1.023003 1.762726
35 35.0 0.086767 0.294563 1.022191 1.798803 
[16]:

inv.plot_inversion_results(
    iters_to_plot=2,
)
../_images/how_to_adhering_to_point_measurements_24_0.png
../_images/how_to_adhering_to_point_measurements_24_1.png
../_images/how_to_adhering_to_point_measurements_24_2.png 
[17]:

_ = ptk.grid_compare(
    true_topography,
    inv.model.topography,
    region=grav_data.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",
    points=constraint_points,
    points_style="x.3c",
)
../_images/how_to_adhering_to_point_measurements_25_0.png 
[18]:

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

rmse_with_weighting = invert4geom.rmse(
    constraint_points.upward - constraint_points.inverted_topography
)
max_error_with_weighting = vd.maxabs(
    constraint_points.upward - constraint_points.inverted_topography
)

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

RMSE at constraints: 2.1 m
max error at constraints: 4.3 m

Method 2: Constraint Point Minimization#

We can define the regional gravity field so that the residual field, which is the input to the inversion, is already ~0 mGal at constraint points. This means the inversion will minimally alter the topography at theses locations. To do this, we estimate the regional gravity misfit by sampling the values at constraints, and interpolating over the entire region. This defines the regional misfit, which when subtracted from the full misfit, gives the residual misfit.

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. From this, we estimate the regional field with the values at constraint points.

[19]:

grav_data.inv.forward_gravity(model)

grav_data.inv.regional_separation(
    method="constraints",
    constraints_df=constraint_points,
)
grav_data

[19]:

<xarray.Dataset> Size: 475kB
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:
    upward                    (northing, easting) float64 40kB 1.001e+03 ... ...
    upper_surface_grav        (northing, easting) float64 40kB 6.966 ... 2.06
    lower_surface_grav        (northing, easting) float64 40kB -0.04507 ... 0...
    gravity_anomaly           (northing, easting) float64 40kB 6.921 ... 2.445
    forward_gravity           (northing, easting) float64 40kB -2.37 ... 1.78
    misfit                    (northing, easting) float64 40kB 9.291 ... 0.6649
    reg                       (northing, easting) float64 40kB -2.798 ... 0.9961
    res                       (northing, easting) float64 40kB 12.09 ... -0.3312
    starting_forward_gravity  (northing, easting) float64 40kB -2.37 ... 1.78
    starting_misfit           (northing, easting) float64 40kB 9.291 ... 0.6649
    starting_reg              (northing, easting) float64 40kB -2.798 ... 0.9961
    starting_res              (northing, easting) float64 40kB 12.09 ... -0.3312
Attributes:
    region:        (0.0, 40000.0, 0.0, 30000.0)
    spacing:       500.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')
[20]:

fig = ptk.plot_grid(
    grav_data.gravity_anomaly,
    fig_height=10,
    title="Topo free disturbance",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSWe", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.misfit,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Misfit",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.reg,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Regional misfit",
    cmap="viridis",
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
)

fig = ptk.plot_grid(
    grav_data.res,
    fig=fig,
    origin_shift="x",
    fig_height=10,
    title="Residual misfit",
    cmap="balance+h0",
    cpt_lims=[-vd.maxabs(grav_data.res), vd.maxabs(grav_data.res)],
    hist=True,
    cbar_label="mGal",
    frame=["nSwE", "xaf10000", "yaf10000"],
    points=constraint_points,
    points_style="x.15c",
)
fig.show()
../_images/how_to_adhering_to_point_measurements_30_0.png

Inversion without weighting grid#

[21]:

# setup the inversion
inv = invert4geom.Inversion(
    grav_data,
    model,
    solver_damping=0.05,
    # set stopping criteria
    max_iterations=300,
    l2_norm_tolerance=0.3,
    delta_l2_norm_tolerance=1.005,
)

[22]:

inv.invert(
    plot_dynamic_convergence=True,
    results_fname="../tmp/non_weighted_inversion",
)
../_images/how_to_adhering_to_point_measurements_33_0.png 
[23]:

inv.stats_df

[23]:
iteration rmse l2_norm delta_l2_norm iter_time_sec
0 0.0 5.277361 2.297251 inf NaN
1 1.0 1.197874 1.094474 2.098954 1.524768
2 2.0 0.359743 0.599786 1.824775 1.625472
3 3.0 0.140696 0.375094 1.599028 1.676509
4 4.0 0.067563 0.259928 1.443068 1.660308 
[24]:

inv.plot_inversion_results(
    iters_to_plot=2,
)
../_images/how_to_adhering_to_point_measurements_35_0.png
../_images/how_to_adhering_to_point_measurements_35_1.png
../_images/how_to_adhering_to_point_measurements_35_2.png 
[25]:

_ = ptk.grid_compare(
    true_topography,
    inv.model.topography,
    region=grav_data.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",
    points=constraint_points,
    points_style="x.3c",
)
../_images/how_to_adhering_to_point_measurements_36_0.png 
[26]:

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

rmse_constraint_point_minimization = invert4geom.rmse(
    constraint_points.upward
    - constraint_points.inverted_topography_constraint_point_minimization
)
max_error_constraint_point_minimization = vd.maxabs(
    constraint_points.upward
    - constraint_points.inverted_topography_constraint_point_minimization
)

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

RMSE at constraints: 5.6 m
max error at constraints: 22.4 m