Juyter Logo

2-linear-inversion-app.ipynb

Seogi Kang
from geoscilabs.inversion.LinearInversionDirect import LinearInversionDirectApp
from ipywidgets import interact, FloatSlider, ToggleButtons, IntSlider, FloatText, IntText
app = LinearInversionDirectApp()

#Linear Inversion App

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

#Purpose

By using a simple decaying and oscillating kernel function, which emulates the physics of electromagnetic (EM) survey, we understand basic concepts of inverting data. Three items that we are going to explore are:

  • Step1: Create a model (m\mathbf{m})
  • Step2: Generate a sensitivity kernel (or matrix), G\mathbf{G}
  • Step3: Simulate data (d=Gm\mathbf{d} = \mathbf{G}\mathbf{m})
  • Step4: All three steps together
  • Step5: Invert the data, and explore inversion results

#Forward problem

Let gj(x)g_j(x) denote the kernel function for jjth datum. With a given model m(x)m(x), jjth datum can be computed by solving following integral equation:

dj=abgj(x)m(x)dxd_j = \int_a^{b} g_j(x) m(x) dx
(1)#

where

gj(x)=epjxcos(2πqjx)g_j(x) = e^{p_jx} cos (2 \pi q_jx)
(2)#

By discretizing gj(x)g_j(x) we obtain

gj(x)=epjxcos(2πqjx)\mathbf{g}_j(\mathbf{x}) = e^{p_j\mathbf{x}} cos (2 \pi q_j \mathbf{x})
(3)#

where

  • gj\mathbf{g}_j: jjth row vector for the sensitivty matrix (1×M1 \times M)
  • x\mathbf{x}: model location (1×M1 \times M)
  • pjp_j: decaying constant (<0)
  • qjq_j: oscillating constant (>0)

By stacking multiple rows of gj\mathbf{g}_j, we obtain sensitivity matrix, G\mathbf{G}:

G=[g1gN]\begin{align} \mathbf{G} = \begin{bmatrix} \mathbf{g}_1\\ \vdots\\ \mathbf{g}_{N} \end{bmatrix} \end{align}
(4)#

Here, the size of the matrix G\mathbf{G} is (N×M)(N \times M). Finally data, d\mathbf{d}, can be written as a linear equation:

d=Gm\mathbf{d} = \mathbf{G}\mathbf{m}
(5)#

where m\mathbf{m} is an inversion model; this is a column vector (M×1M \times 1).

In real measurments, there will be various noises source, and hence observation, dobs\mathbf{d}^{obs}, can be written as

dobs=Gm+noise\mathbf{d}^{obs} = \mathbf{G}\mathbf{m} + \mathbf{noise}
(6)#

#Step1: Create a model, m\mathbf{m}

The model mm is a function defined on the interval (-2,2). Here we generate a model that is the sum of a: (a) background mrefm_{ref}, (b) box car m1m_1 and (c) Gaussian m2m_2. The box car is defined by

  • m$_{background}$ : amplitude of the background
  • m1 : amplitude
  • $m1_{center}$ : center
  • $m1_{width}$ : width the Gaussian is defined by
  • m2 : amplitude
  • $m2_{center}$ : center
  • $m2_{sigma}$ : width of Gaussian (as defined by a standard deviation)
  • M: # of model parameters
Q_model = app.interact_plot_model()

#Step2: Generate a sensitivity kernel (or matrix), G\mathbf{G}

By using the following app, we explore each row vector of the sensitivity matrix, gj\mathbf{g}_j. Parameters of the apps are:

  • M: # of model parameters
  • N: # of data
  • p: decaying constant (<0)
  • q: oscillating constant (>0)
  • ymin: maximum limit for y-axis
  • ymax: minimum limit for y-axis
  • show_singular: show singualr values
Q_kernel = app.interact_plot_G()

#Step3: Simulate data

The jj-th datum is the inner product of the jj-th kernel gj(x)g_j(x) and the model m(x)m(x). In discrete form it can be written as the dot product of the vector gjg_j and the model vector mm.

#
dj=gjm d_j = \mathbf{g}_j \mathbf{m}
(7)#

If there are NN data, these data can be written as a column vector, d\mathbf{d}:

d=Gm=[d1dN]\begin{align} \mathbf{d} = \mathbf{G}\mathbf{m} = \begin{bmatrix} d_1\\ \vdots\\ d_{N} \end{bmatrix} \end{align}
(8)#

#Adding Noise

Observational data are always contaminated with noise. Here we add Gaussian noise N(0,ϵ)N(0,\epsilon) (zero mean and standard deviation σ\sigma). Here we choose

ϵ=%d+floor\epsilon = \% |d| + \text{floor}
(9)#
Q_data = app.interact_plot_data()

#Step4: All three steps together

app.interact_plot_all_three_together()

#Inverse Problem

In the inverse problem we attempt to find the model m\mathbf{m} that gave rise to the observational data dobs\mathbf{d}^{obs}. The inverse problem is formulated as an optimization problem:

minimize   ϕ(m)=ϕd(m)+βϕm(m)\text{minimize} \ \ \ \phi(\mathbf{m}) = \phi_d(\mathbf{m}) + \beta \phi_m(\mathbf{m})
(10)#

where

  • ϕd\phi_d: data misfit
  • ϕm\phi_m: model regularization
  • β\beta: trade-off (or Tikhonov) parameter 0<β<0<\beta<\infty

Data misfit is defined as

ϕd=j=1N(gjmdjobsϵj)2\phi_d = \sum_{j=1}^{N}\Big(\frac{\mathbf{g}_j\mathbf{m}-d^{obs}_j}{\epsilon_j}\Big)^2
(11)#

where ϵj\epsilon_j is an estimate of the standard deviation of the jjth datum.

The model regularization term, ϕm\phi_m, can be written as

ϕm(m)=αs(mmref)dx+αx(dmdx)dx\phi_m(\mathbf{m}) = \alpha_s \int (\mathbf{m}-\mathbf{m}_{ref}) dx + \alpha_x \int (\frac{d \mathbf{m}}{dx}) dx
(12)#

The first term is referred to as the "smallness" term. Minimizing this generates a model that is close to a reference model mrefm_{ref}. The second term penalizes roughness of the model. It is generically referred to as a "flattest" or "smoothness" term.

#Step5: Invert the data, and explore inversion 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.

#Parameters

  • mode: Run or Explore

    • Run: Each click of the app, will run n_beta times of inversion
    • Explore: Not running inversions, but explore result of the inversions
  • noise option: error contaminated or clean data

#Misfit

  • percent: percentage of the uncertainty (%)

  • floor: floor of the uncertainty (%)

  • chifact: chi factor for stopping criteria (when ϕd=N\phi_d^{\ast}=N \rightarrow chifact=1)

#Model norm

  • mref: reference model

  • alpha_s: αs\alpha_s for smallness

  • alpha_x: αx\alpha_x for smoothness

#Beta

  • beta_min: minimum β\beta

  • beta_max: maximum β\beta

  • n_beta: the number of β\beta

#Plotting options

  • data: obs & pred or normalized misfit

    • obs & pred: show observed and predicted data
    • normalized misfit: show normalized misfit
  • tikhonov: phi_d & phi_m or phi_d vs phi_m

    • phi_d & phi_m: show ϕd\phi_d and ϕm\phi_m as a function of β\beta
    • phi_d vs phi_m: show tikhonov curve
  • i_beta: i-th β\beta value

  • scale: linear or log

    • linear: linear scale for plotting the third panel
    • log: log scale for plotting the third panel
app.interact_plot_inversion()
References
  1. 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