import numpy as np
from geoscilabs.inversion.TomographyInversion import TomographyInversionApp
import matplotlib.pyplot as plt
from ipywidgets import interact, IntSlider
From the Linear inversion notebook (1D), we learned important apsects about the linear inversion. However, real world geophysical inverse problem are not often 1D, but multidimensional (2D or 3D), and this extension of dimension allows us to put more apriori (or geologic) information through the regularization term. In this notebook, we explore these multidimensional aspects of the linear inversion by using 2D traveltime croswell tomography example.
This notebook includes four steps:
Here we set up a velocity model using a following app. Controlling parameters of the app are:
set mesh
: use active only when you want to change the 2D mesh
add block
: use active when you want to add block (if not stay inactive)
model type
: background or block
show grid?
: show grid of the mesh
v0
: velocity of the background
v1
: velocity of the block
xc
: x center of the block
yc
: y center of the block
dx
: width of the block
dy
: thickness of the block
nx
: # of cells in x-direction (this is only active when set_mesh=active
)
ny
: # of cells in y-direction (this is only active when set_mesh=active
)
Related parameters for this task are: set mesh
, nx
, ny
.
Size of the 2D domain are fixed to 200m $\times$ 400m, but the number of cells in each direction can be changed such that you can alter size of the cells in each direction. When you change either nx
or ny
make sure you choose set mesh=active
otherwise set mesh=inactive
. Note that once mesh setup is changed, velocity model is reset to a background value (v0
).
Although you can change the location, size, and velocity of the block there are few rules that you need to follow to do so.
If you want to change the parameter of the block: first set add block=active
then change parameters of the block (v1
, xc
, zc
, dx
, dy
)
Once you changed the parameters, make sure first choose model type=background
then change that to model type=block
You can also add multiple blocks. To add a block follow below steps:
Set add block=inactive
, then change the parameter of the new block using: v1
, xc
, zc
, dx
, dy
. Velocity model will not change, but you can see the white lines which illustrate boundary of the new block.
Once you are happy with the new block, set add block=active
, then velocity model will be updated with the new block that you have set.
Repeat 1 and 2 if you want to add more blocks.
app = TomographyInversionApp()
app.interact_plot_model()
Within this app, by using the velocity model set up, we simulate traveltime tomography data and add Gaussian noise. This syntehtic data set will be used in the following inversion module. Controlling parameters are:
percent (%)
: percentage of the Gaussian noisefloor (s)
: floor of the Gaussain noiserandom seed
: seed to generate random variables having normal distributiontx_rx_plane
or histogram
: choice of the plotting dataupdate
: this buttion is for the interactin between the first app. If the velocity model is changed by altering the first app, you can simply click update
to run simulation again with the updated velocity model.app.interact_data()
maxIter
: maximum number of iterationm0
: initial modelmref
: reference modelpercentage
: percent standard deviation for the uncertaintyfloor
: floor value for the uncertaintychifact
: chifactor for the target misfitbeta0_ratio
: ratio to set the initial betacoolingFactor
: cooling factor to cool betan_iter_per_beta
: # of interation for each beta valuealpha_s
: $\alpha_s$alpha_x
: $\alpha_x$alpha_z
: $\alpha_z$use_target
: use target misfit as a stopping criteria or notmodel, pred, save = app.run_inversion(
maxIter=20,
m0=1./1000.,
mref=1./1250.,
percentage=0,
floor=0.005,
chifact=1,
beta0_ratio=1e2,
coolingFactor=2,
n_iter_per_beta=1,
alpha_s=1/(app.mesh_prop.hx.min())**2 * 1e4,
alpha_x=1,
alpha_z=1,
use_target=True
)
app.interact_model_inversion(model, clim=None)
app.interact_data_inversion(pred)
save.plot_misfit_curves()
app.plot_tikhonov_curves(save);