Linear Tikhonov Inversion¶
Note This Jupyter notebook is a linked copy of the LinearTikhonovInversion_App and contains additional cells that reproduce the specific examples presented throughout 2 Linear Tikhonov Inversion.
from geoscilabs.inversion.LinearInversionDirect import LinearInversionDirectApp from ipywidgets import interact, FloatSlider, ToggleButtons, IntSlider, FloatText, IntText import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams['font.size'] = 14
app = LinearInversionDirectApp()
This app is based upon the inversion tutorial: "INVERSION FOR APPLIED GEOPHYSICS" by Oldenburg and Li (2005).
Douglas W. Oldenburg and Yaoguo Li (2005) 5. Inversion for Applied Geophysics: A Tutorial. Near-Surface Geophysics: pp. 89-150. eISBN: 978-1-56080-171-9 print ISBN: 978-1-56080-130-6 https://doi.org/10.1190/1.9781560801719.ch5
We illustrate how a generic linear inverse problem can be solved using a Tikhonov approach. The default parameters provided for the Forward and Inverse problems below generate a reasonable example for illustrating the inversion but the learning comes when these parameters are changed and outcomes are observed.
The app is divided into two sections:
Mathematical Background for the Forward Problem
Step 1: Create a model, .
Step 2: Generate a sensitivity matrix .
Step 3: Simulate data () and add noise.
These steps are explored individually but additional text is given in 2 Linear Tikhonov Inversion. For convenience, the widgets used to carry out all three steps are consolidated at the end of the section. A brief mathematical description is also provided.
Mathematical Background for the Inverse Problem
Step 4: Invert the data, and explore the results
Here we provide widgets to adjust the parameters for the inverse problem. Some basic information is provided but details about the parameters are provided in the text 2 Linear Tikhonov Inversion.
Mathematical Background for the Forward Problem¶
Let denote the kernel function for th datum. With a given model , the th datum can be computed by solving following integral equation:
is the kernel function. By integrating over cells of width and using the midpoint rule cell we obtain the sensitivities
: th row vector for the sensitivty matrix ()
: model location ()
: decaying constant (<0)
: oscillating constant (>0)
By stacking multiple rows of , we obtain sensitivity matrix, :
Here, the size of the matrix is . Finally data, , can be written as a linear equation:
where is an inversion model; this is a column vector ().
In real measurments, there will be various noise sources, and hence observation, , can be written as
Step 1: Create a model, ¶
The model is a function defined on the interval [0,1] and discretized into equal intervals. It is the sum of a: (a) background , (b) box car and (c) Gaussian .
m_background: background value
The box car is defined by
The Gaussian is defined by
m2_sigma: width of Gaussian (as defined by a standard deviation )
M: number of model parameters
Q_model = app.interact_plot_model()
app.return_axis = True ax = app.plot_model_only( m_background = 0., m1 = 1, m1_center = 0.2, dm1 = 0.2, m2 = 2, m2_center = 0.75, sigma_2 = 0.07, M=100 ) ax.set_xlabel("x") app.return_axis = False
Step 2: Generate a sensitivity matrix, ¶
By using the following app, we explore each row vector, , of the kernel or sensitivity matrix , . Parameters of the apps are:
M: number of model parameters
N: number of data
pmax: minimum and maximum of the -length range of decaying constant values (<0)
qmax: minimum and maximum of the -length range of oscillating constant values (>0)
ymax: minimum and maximum of the y-axis
Q_kernel = app.interact_plot_G()
#plot for 3 kernels app.return_axis = True axs = app.plot_G( N=3, M=100, pmin=0, pmax=-2, qmin=1, qmax=3, scale='log', fixed=False, ymin=-0.005, ymax=0.011, ) axs.set_title("kernel functions") app.return_axis = False;
Step 3: Simulate data, , and add noise¶
The -th datum is the inner product of the -th kernel and the model . In discrete form it can be written as the dot product of the vector and the model vector .
If there are data, these data can be written as a column vector, :
Observational data are always contaminated with noise. Here we add Gaussian noise (zero mean and standard deviation ). Here we choose
Q_data = app.interact_plot_data()
#plot accurate data app.return_axis = True ax = app.plot_data_only( add_noise=True, percentage=0, floor=0.03, ) app.return_axis = False
Composite Widget for Forward Modelling¶
# Default parameters: accurate data app.return_axis = True axs = app.plot_model( m_background = 0, m1 = 1, m2 = 2, m1_center = 0.2, dm1 = 0.2, m2_center = 0.75, sigma_2 = 0.07, percentage = 0, floor = 0.0, pmin=-0.25, pmax=-3, qmin=0., qmax=5, ) axs.set_title("Model") app.return_axis = False
# Default parameters: noisey data app.return_axis = True axs = app.plot_model( m_background = 0, m1 = 1, m2 = 2, m1_center = 0.2, dm1 = 0.2, m2_center = 0.75, sigma_2 = 0.07, percentage = 0, floor = 0.03, pmin=-0.25, pmax=-3, qmin=0., qmax=5, ) axs.set_title("Noisy data") app.return_axis = False
Mathematical Background for the Inverse Problem¶
In the inverse problem we attempt to find the model that gave rise to the observational data . The inverse problem is formulated as an optimization problem:
- : data misfit
- : model regularization
- : trade-off (Tikhonov) parameter
Data misfit is defined as
where is an estimate of the standard deviation of the th datum.
The model regularization term, , can be written as
The first term is referred to as the "smallness" term. Minimizing this generates a model that is close to a reference model . The second term penalizes roughness of the model. It is generically referred to as a "flattest" or "smoothness" term.
Step 4: Invert the data, and explore the results¶
In the inverse problem we define parameters needed to evaluate the data misfit and the model regularization terms. We then deal with parameters associated with the inversion.
Run: Each click of the app, will run
Explore: Not running inversions, but explore result of the previously run inversions
percent: estiamte uncertainty as a percentage of the data (%)
floor: estimate uncertainty floor
chifact: chi factor for stopping criteria (when
mref: reference model
alpha_s: weight for smallness term
alpha_x: weight for smoothness term
n_beta: the number of
obs & predor
obs & pred: show observed and predicted data
normalized misfit: show normalized misfit
phi_d & phi_mor
phi_d vs phi_m
phi_d & phi_m: show and as a function of
phi_d vs phi_m: show tikhonov curve
i_beta: i-th value
linear: linear scale for plotting the third panel
log: log scale for plotting the third panel
app.return_axis = True axs = app.plot_inversion( mode="Run", #"Explore" mref=0.0, percentage=app.percentage, floor=app.floor, beta_min=1e-3, beta_max=1e5, n_beta=81, alpha_s=1, alpha_x=0, tikhonov="phi_d & phi_m", data_option="obs & pred", scale="log", i_beta=0, chifact=1, ) app.return_axis = False
- Oldenburg, D. W., & Li, Y. (2005). 5. Inversion for Applied Geophysics: A Tutorial. In Near-Surface Geophysics (pp. 89–150). Society of Exploration Geophysicists. 10.1190/1.9781560801719.ch5