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 (DOI: 10.1093/gji/ggw390). See their paper for a detailed discussion.
The main idea is to re-sample the gravity data at a finer resolution and seperate 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.
Import packages#
[38]:
from __future__ import annotations
%load_ext autoreload
%autoreload 2
import logging
import numpy as np
import verde as vd
import xarray as xr
from polartoolkit import utils as polar_utils
from invert4geom import cross_validation, inversion, plotting, synthetic, utils
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Create observed gravity data#
True topography#
[39]:
# set grid parameters
spacing = 1000
region = [0, 40000, 0, 30000]
# create synthetic topography data
true_topography = synthetic.synthetic_topography_simple(
spacing,
region,
)
Prism layer#
[40]:
# the density contrast is between rock (~2670 kg/m3) and air (~1 kg/m3)
density_contrast = 2670 - 1
# prisms are created between the mean topography value and the height of the topography
zref = true_topography.values.mean()
# prisms above zref have positive density contrast and prisms below zref have negative
# density contrast
density = xr.where(true_topography >= zref, density_contrast, -density_contrast)
# create layer of prisms
prisms = utils.grids_to_prisms(
true_topography,
zref,
density=density,
)
Forward gravity of prism layer#
[41]:
# 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(
region=region,
spacing=spacing,
pixel_register=False,
extra_coords=1000, # survey elevation
)
# grid the coordinates
observations = vd.make_xarray_grid(
(coords[0], coords[1]),
data=coords[2],
data_names="upward",
dims=("northing", "easting"),
).upward
grav_df = vd.grid_to_table(observations)
grav_df["grav"] = prisms.prism_layer.gravity(
coordinates=(
grav_df.easting,
grav_df.northing,
grav_df.upward,
),
field="g_z",
progressbar=True,
)
grav_df
[41]:
| northing | easting | upward | grav | |
|---|---|---|---|---|
| 0 | 0.0 | 0.0 | 1000.0 | 9.534643 |
| 1 | 0.0 | 1000.0 | 1000.0 | 10.422834 |
| 2 | 0.0 | 2000.0 | 1000.0 | 9.949973 |
| 3 | 0.0 | 3000.0 | 1000.0 | 9.269279 |
| 4 | 0.0 | 4000.0 | 1000.0 | 8.532160 |
| ... | ... | ... | ... | ... |
| 1266 | 30000.0 | 36000.0 | 1000.0 | 3.332716 |
| 1267 | 30000.0 | 37000.0 | 1000.0 | 3.330307 |
| 1268 | 30000.0 | 38000.0 | 1000.0 | 3.335438 |
| 1269 | 30000.0 | 39000.0 | 1000.0 | 3.300721 |
| 1270 | 30000.0 | 40000.0 | 1000.0 | 2.858299 |
1271 rows × 4 columns
Resample to create testing / training data sets#
Note that it is imporant to do this resampling on the original observed gravity data before you calculate your misfit!
[42]:
# resample to half spacing
grav_df = cross_validation.resample_with_test_points(spacing, grav_df, region)
grav_df
[42]:
| northing | easting | test | upward | grav | |
|---|---|---|---|---|---|
| 0 | 0.0 | 0.0 | False | 1000.0 | 9.534643 |
| 1 | 0.0 | 500.0 | True | 1000.0 | 10.063805 |
| 2 | 0.0 | 1000.0 | False | 1000.0 | 10.422834 |
| 3 | 0.0 | 1500.0 | True | 1000.0 | 10.284459 |
| 4 | 0.0 | 2000.0 | False | 1000.0 | 9.949973 |
| ... | ... | ... | ... | ... | ... |
| 4936 | 30000.0 | 38000.0 | False | 1000.0 | 3.335438 |
| 4937 | 30000.0 | 38500.0 | True | 1000.0 | 3.346051 |
| 4938 | 30000.0 | 39000.0 | False | 1000.0 | 3.300721 |
| 4939 | 30000.0 | 39500.0 | True | 1000.0 | 3.104991 |
| 4940 | 30000.0 | 40000.0 | False | 1000.0 | 2.858299 |
4941 rows × 5 columns
[43]:
# contaminate gravity with 0.5 mGal of random noise
grav_df["observed_grav"], stddev = synthetic.contaminate(
grav_df.grav,
stddev=0.5,
percent=False,
seed=0,
)
grav_df.set_index(["northing", "easting"]).to_xarray().observed_grav.plot()
INFO:root:Standard deviation used for noise: [0.5]
[43]:
<matplotlib.collections.QuadMesh at 0x7fdf6e274740>
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.
[44]:
# create flat topography grid with a constant height equal the mean of the starting
# topography (zref)
starting_topography = xr.full_like(true_topography, zref)
# prisms above zref have positive density contrast and prisms below zref have negative
# density contrast
density = xr.where(starting_topography >= zref, density_contrast, -density_contrast)
# create layer of prisms
starting_prisms = utils.grids_to_prisms(
starting_topography,
zref,
density=density,
)
# gravity of starting model is 0 since its flat, so observed_grav = misfit
grav_df["misfit"] = grav_df["observed_grav"]
# 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 the
# full misfit is equal to the residual misfit.
# set regional misfit to 0
grav_df["reg"] = 0
# set the residual misfit to the full misfit
grav_df["res"] = grav_df.misfit
Get Cross Validation Score#
We will perform an inversion with a single damping value and calculate a score for it.
[45]:
# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.WARNING)
# set kwargs to pass to the inversion
kwargs = {
"input_grav_column": "observed_grav",
"prism_layer": starting_prisms,
"deriv_type": "annulus",
"zref": zref,
"density_contrast": density_contrast,
"solver_damping": 0.1,
# set stopping criteria
"max_iterations": 30,
"l2_norm_tolerance": 0.5,
"delta_l2_norm_tolerance": 1.005,
}
# run inversion, calculate the score, and plot the predicted and observed gravity for
# the testing dataset
score = cross_validation.grav_cv_score(
training_data=grav_df[grav_df.test == False], # noqa: E712
testing_data=grav_df[grav_df.test == True], # noqa: E712
progressbar=True,
plot=True,
**kwargs,
)
score
[45]:
0.5927306328269795
Cross Validation#
Now we can repeat this with a range of damping parameters to find the optimal (lowest) score.
[46]:
# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.WARNING)
# set which damping parameters to include
dampings = np.logspace(-3, 0, 8)
best_damping, _, _, scores = cross_validation.grav_optimal_parameter(
training_data=grav_df[grav_df.test == False], # noqa: E712
testing_data=grav_df[grav_df.test == True], # noqa: E712
param_to_test=("solver_damping", dampings),
progressbar=True,
plot_grids=True,
plot_cv=False,
verbose=True,
**kwargs,
)
INFO:root:Parameter value: 0.001 -> Score: 2.346310520871798
INFO:root:Parameter value: 0.0026826957952797246 -> Score: 2.762671261658884
INFO:root:Parameter value: 0.0071968567300115215 -> Score: 0.7640486003535428
INFO:root:Parameter value: 0.019306977288832496 -> Score: 0.6696671073768633
INFO:root:Parameter value: 0.0517947467923121 -> Score: 0.5940119946779362
INFO:root:Parameter value: 0.13894954943731375 -> Score: 0.5892038807688588
INFO:root:Parameter value: 0.3727593720314938 -> Score: 1.9643157211873323
INFO:root:Parameter value: 1.0 -> Score: 7.033086888044258
INFO:root:Best score of 0.5892038807688588 with parameter value=0.13894954943731375
[47]:
# Compare the scores and the damping values
plotting.plot_cv_scores(
scores,
dampings,
param_name="Damping",
logx=True,
logy=True,
)
Run inversion with optimal value#
[48]:
# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.INFO)
# make new kwargs without solver damping
new_kwargs = {
key: value
for key, value in kwargs.items()
if key
not in [
"solver_damping",
]
}
results = inversion.run_inversion(
input_grav=grav_df[grav_df.test == False], # noqa: E712, we can exlude the testing data now
plot_convergence=True,
solver_damping=best_damping,
**new_kwargs,
)
# collect the results
topo_results, grav_results, parameters, elapsed_time = results
INFO:root:starting inversion
INFO:root:extracted prism spacing is 1000.0
INFO:root:
####################################
iteration 1
INFO:root:Layer correction median: -3.5984 m, RMSE:20.8845 m
INFO:root:updated misfit RMSE: 5.2604
INFO:root:updated L2-norm: 2.2935, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1681, tolerance: 1.005
INFO:root:
####################################
iteration 2
INFO:root:Layer correction median: 10.1358 m, RMSE:15.0242 m
INFO:root:updated misfit RMSE: 3.9019
INFO:root:updated L2-norm: 1.9753, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1611, tolerance: 1.005
INFO:root:
####################################
iteration 3
INFO:root:Layer correction median: 7.3421 m, RMSE:10.9064 m
INFO:root:updated misfit RMSE: 2.9365
INFO:root:updated L2-norm: 1.7136, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1527, tolerance: 1.005
INFO:root:
####################################
iteration 4
INFO:root:Layer correction median: 5.3415 m, RMSE:8.0026 m
INFO:root:updated misfit RMSE: 2.2476
INFO:root:updated L2-norm: 1.4992, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.143, tolerance: 1.005
INFO:root:
####################################
iteration 5
INFO:root:Layer correction median: 3.8694 m, RMSE:5.9435 m
INFO:root:updated misfit RMSE: 1.7535
INFO:root:updated L2-norm: 1.3242, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1322, tolerance: 1.005
INFO:root:
####################################
iteration 6
INFO:root:Layer correction median: 2.8447 m, RMSE:4.4738 m
INFO:root:updated misfit RMSE: 1.3971
INFO:root:updated L2-norm: 1.182, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1203, tolerance: 1.005
INFO:root:
####################################
iteration 7
INFO:root:Layer correction median: 2.0905 m, RMSE:3.417 m
INFO:root:updated misfit RMSE: 1.1383
INFO:root:updated L2-norm: 1.0669, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1079, tolerance: 1.005
INFO:root:
####################################
iteration 8
INFO:root:Layer correction median: 1.5328 m, RMSE:2.6507 m
INFO:root:updated misfit RMSE: 0.9487
INFO:root:updated L2-norm: 0.974, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0953, tolerance: 1.005
INFO:root:
####################################
iteration 9
INFO:root:Layer correction median: 1.1087 m, RMSE:2.0899 m
INFO:root:updated misfit RMSE: 0.8085
INFO:root:updated L2-norm: 0.8991, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0833, tolerance: 1.005
INFO:root:
####################################
iteration 10
INFO:root:Layer correction median: 0.8308 m, RMSE:1.6753 m
INFO:root:updated misfit RMSE: 0.7033
INFO:root:updated L2-norm: 0.8386, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0722, tolerance: 1.005
INFO:root:
####################################
iteration 11
INFO:root:Layer correction median: 0.6105 m, RMSE:1.3652 m
INFO:root:updated misfit RMSE: 0.6233
INFO:root:updated L2-norm: 0.7895, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0622, tolerance: 1.005
INFO:root:
####################################
iteration 12
INFO:root:Layer correction median: 0.4615 m, RMSE:1.1303 m
INFO:root:updated misfit RMSE: 0.5616
INFO:root:updated L2-norm: 0.7494, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0536, tolerance: 1.005
INFO:root:
####################################
iteration 13
INFO:root:Layer correction median: 0.3513 m, RMSE:0.9499 m
INFO:root:updated misfit RMSE: 0.5131
INFO:root:updated L2-norm: 0.7163, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0462, tolerance: 1.005
INFO:root:
####################################
iteration 14
INFO:root:Layer correction median: 0.2603 m, RMSE:0.8094 m
INFO:root:updated misfit RMSE: 0.4744
INFO:root:updated L2-norm: 0.6887, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.04, tolerance: 1.005
INFO:root:
####################################
iteration 15
INFO:root:Layer correction median: 0.2037 m, RMSE:0.6986 m
INFO:root:updated misfit RMSE: 0.4429
INFO:root:updated L2-norm: 0.6655, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0348, tolerance: 1.005
INFO:root:
####################################
iteration 16
INFO:root:Layer correction median: 0.1592 m, RMSE:0.6099 m
INFO:root:updated misfit RMSE: 0.4171
INFO:root:updated L2-norm: 0.6458, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0305, tolerance: 1.005
INFO:root:
####################################
iteration 17
INFO:root:Layer correction median: 0.1263 m, RMSE:0.538 m
INFO:root:updated misfit RMSE: 0.3955
INFO:root:updated L2-norm: 0.6289, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.027, tolerance: 1.005
INFO:root:
####################################
iteration 18
INFO:root:Layer correction median: 0.0998 m, RMSE:0.4791 m
INFO:root:updated misfit RMSE: 0.3772
INFO:root:updated L2-norm: 0.6141, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.024, tolerance: 1.005
INFO:root:
####################################
iteration 19
INFO:root:Layer correction median: 0.0842 m, RMSE:0.4303 m
INFO:root:updated misfit RMSE: 0.3615
INFO:root:updated L2-norm: 0.6012, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0215, tolerance: 1.005
INFO:root:
####################################
iteration 20
INFO:root:Layer correction median: 0.0675 m, RMSE:0.3894 m
INFO:root:updated misfit RMSE: 0.3479
INFO:root:updated L2-norm: 0.5898, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0194, tolerance: 1.005
INFO:root:
####################################
iteration 21
INFO:root:Layer correction median: 0.0554 m, RMSE:0.3549 m
INFO:root:updated misfit RMSE: 0.336
INFO:root:updated L2-norm: 0.5796, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0176, tolerance: 1.005
INFO:root:
####################################
iteration 22
INFO:root:Layer correction median: 0.047 m, RMSE:0.3256 m
INFO:root:updated misfit RMSE: 0.3254
INFO:root:updated L2-norm: 0.5705, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0161, tolerance: 1.005
INFO:root:
####################################
iteration 23
INFO:root:Layer correction median: 0.0419 m, RMSE:0.3004 m
INFO:root:updated misfit RMSE: 0.316
INFO:root:updated L2-norm: 0.5622, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0148, tolerance: 1.005
INFO:root:
####################################
iteration 24
INFO:root:Layer correction median: 0.0368 m, RMSE:0.2787 m
INFO:root:updated misfit RMSE: 0.3076
INFO:root:updated L2-norm: 0.5546, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0137, tolerance: 1.005
INFO:root:
####################################
iteration 25
INFO:root:Layer correction median: 0.0333 m, RMSE:0.2598 m
INFO:root:updated misfit RMSE: 0.2999
INFO:root:updated L2-norm: 0.5476, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0127, tolerance: 1.005
INFO:root:
####################################
iteration 26
INFO:root:Layer correction median: 0.0287 m, RMSE:0.2433 m
INFO:root:updated misfit RMSE: 0.2929
INFO:root:updated L2-norm: 0.5412, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0119, tolerance: 1.005
INFO:root:
####################################
iteration 27
INFO:root:Layer correction median: 0.0244 m, RMSE:0.2288 m
INFO:root:updated misfit RMSE: 0.2864
INFO:root:updated L2-norm: 0.5352, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0112, tolerance: 1.005
INFO:root:
####################################
iteration 28
INFO:root:Layer correction median: 0.0235 m, RMSE:0.216 m
INFO:root:updated misfit RMSE: 0.2805
INFO:root:updated L2-norm: 0.5296, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0106, tolerance: 1.005
INFO:root:
####################################
iteration 29
INFO:root:Layer correction median: 0.0227 m, RMSE:0.2046 m
INFO:root:updated misfit RMSE: 0.2749
INFO:root:updated L2-norm: 0.5243, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.01, tolerance: 1.005
INFO:root:
####################################
iteration 30
INFO:root:Layer correction median: 0.0208 m, RMSE:0.1945 m
INFO:root:updated misfit RMSE: 0.2698
INFO:root:updated L2-norm: 0.5194, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0095, tolerance: 1.005
INFO:root:
Inversion terminated after 30 iterations with L2-norm=0.52 because maximum number of iterations (30) reached.
[49]:
plotting.plot_inversion_results(
grav_results,
topo_results,
parameters,
region,
iters_to_plot=4,
plot_iter_results=True,
plot_topo_results=True,
plot_grav_results=True,
)
[50]:
final_topography = topo_results.set_index(["northing", "easting"]).to_xarray().topo
_ = polar_utils.grd_compare(
true_topography,
final_topography,
# plot_type="xarray",
plot=True,
grid1_name="True topography",
grid2_name="Inverted topography",
robust=True,
hist=True,
inset=False,
verbose="q",
title="difference",
grounding_line=False,
reverse_cpt=True,
cmap="rain",
# diff_lims=(-20, 20),
)
Run inversion with poor choice of damping value#
[52]:
# set Python's logging level
logger = logging.getLogger()
logger.setLevel(logging.INFO)
# make new kwargs without solver damping
new_kwargs = {
key: value
for key, value in kwargs.items()
if key
not in [
"solver_damping",
]
}
results = inversion.run_inversion(
input_grav=grav_df[grav_df.test == False], # noqa: E712, we can exlude the testing data now
plot_convergence=True,
solver_damping=0.01,
**new_kwargs,
)
# collect the results
topo_results, grav_results, parameters, elapsed_time = results
INFO:root:starting inversion
INFO:root:extracted prism spacing is 1000.0
INFO:root:
####################################
iteration 1
INFO:root:Layer correction median: -10.7028 m, RMSE:104.1676 m
INFO:root:updated misfit RMSE: 2.0705
INFO:root:updated L2-norm: 1.4389, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.8619, tolerance: 1.005
INFO:root:
####################################
iteration 2
INFO:root:Layer correction median: -18.5924 m, RMSE:35.4002 m
INFO:root:updated misfit RMSE: 1.0033
INFO:root:updated L2-norm: 1.0017, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.4365, tolerance: 1.005
INFO:root:
####################################
iteration 3
INFO:root:Layer correction median: 4.4712 m, RMSE:19.5611 m
INFO:root:updated misfit RMSE: 0.5396
INFO:root:updated L2-norm: 0.7346, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.3636, tolerance: 1.005
INFO:root:
####################################
iteration 4
INFO:root:Layer correction median: -1.1404 m, RMSE:13.1625 m
INFO:root:updated misfit RMSE: 0.3814
INFO:root:updated L2-norm: 0.6176, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1894, tolerance: 1.005
INFO:root:
####################################
iteration 5
INFO:root:Layer correction median: 0.296 m, RMSE:10.7714 m
INFO:root:updated misfit RMSE: 0.3065
INFO:root:updated L2-norm: 0.5536, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.1156, tolerance: 1.005
INFO:root:
####################################
iteration 6
INFO:root:Layer correction median: -0.084 m, RMSE:9.5231 m
INFO:root:updated misfit RMSE: 0.2655
INFO:root:updated L2-norm: 0.5152, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.0744, tolerance: 1.005
INFO:root:
####################################
iteration 7
INFO:root:Layer correction median: 0.0353 m, RMSE:8.9668 m
INFO:root:updated misfit RMSE: 0.2527
INFO:root:updated L2-norm: 0.5027, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.025, tolerance: 1.005
INFO:root:
####################################
iteration 8
INFO:root:Layer correction median: -0.0097 m, RMSE:8.724 m
INFO:root:updated misfit RMSE: 0.24
INFO:root:updated L2-norm: 0.4899, tolerance: 0.5
INFO:root:updated delta L2-norm : 1.026, tolerance: 1.005
INFO:root:
Inversion terminated after 8 iterations because L2-norm (0.4899319229639865) was less then set tolerance: 0.5
Change parameter 'l2_norm_tolerance' if desired.
[53]:
plotting.plot_inversion_results(
grav_results,
topo_results,
parameters,
region,
iters_to_plot=4,
plot_iter_results=True,
plot_topo_results=True,
plot_grav_results=True,
)
[54]:
final_topography = topo_results.set_index(["northing", "easting"]).to_xarray().topo
_ = polar_utils.grd_compare(
true_topography,
final_topography,
# plot_type="xarray",
plot=True,
grid1_name="True topography",
grid2_name="Inverted topography",
robust=True,
hist=True,
inset=False,
verbose="q",
title="difference",
grounding_line=False,
reverse_cpt=True,
cmap="rain",
# diff_lims=(-20, 20),
)
Comparing this inversion with a damping value of 0.01 to the inversion with the optimal damping value of ~0.14 shows that a low damping value allows too much noise through to the final inverted topography.
In the past example (simple_inversion.ipynb), 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_cross_validation.ipynb, will show how to perform a similar cross validation for choosing a density contrast value.