In this sequence of computer problems you will implement a simple example of seismic imaging by reverse time migration. This problem sequence will guide you step-by-step through implementing everything you need to do 1D seismic imaging.
These exercises are best solved using a single python script file, which I will call cp3_sol.py (but you can name it whatever you wish).
Will will need access to the functions you wrote for the previous computer problems (cp1_sol.py, cp2_sol.py). You can either copy and paste them into your new file or you can use a python import statement to get access to them in your new file.
I suggest you develop your solutions inside this file, and then run the script in an IPython console. This is accomplished by starting an IPython console, from an command prompt, with the command ipython --pylab.
Then, from the IPython console, you can run your script with the command
%run cp3_solAs before,yYou will need a wave speed model for these problems. These are provided in msri_tools.py and can be accessed by downloading the file to the directory you will work in. From an ipython --pylab console you can access the model function by executing the command from msri_tools import model. You can also put the line
from msri_tools import modelat the top of your python scripts.
In this exercise you will solve the linear forward modeling equations and verify that the adjoint conditions hold for your code. The definition of the adjoint of a linear operator states that
\[ \left< F\delta m, d \right> _{\mathcal{D}} = \left< \delta m, F^{*}d \right> _{\mathcal{M}}, \]
where \(\left< \cdot, \cdot \right> _{\mathcal{D}}\) and \(\left< \cdot, \cdot \right> _{\mathcal{M}}\) are inner products in the data and model spaces, respectively. One test to see if your migration operator \(F^*\) is working correctly is to test if this relationship holds to high precision for any pair of \(d\) and \(\delta m\).
Finally, you will implement a simple gradient descent algorithm for minimizing the objective function
\[ J(m) = \frac{1}{2}\int_0^T||d_{x_r}(t) - \mathbf{S}\mathcal{F}[m](t)||_2^2 \textrm{d}t, \]
to 'solve' the seismic imaging problem.
To implement the adjoint test, you will need to solve the linear modeling equations, which are derived in the notes,
\[ \begin{align*} (\frac{1}{c_0(x)^2}\partial_{tt}-\partial_{xx})u_1(x,t) & = -\delta m(x) \partial_{tt}u_0(x,t), \\\\ (\frac{1}{c_0(x)}\partial_t-\partial_x)u_1(0,t) & = 0, \\\\ (\frac{1}{c_0(x)}\partial_t+\partial_x)u_1(1,t) & = 0, \\\\ u_1(x,t) & = 0 \quad\text{for}\quad t \le 0, \end{align*} \]
where \(u_1\) is the Born scattered field and equivalently \(F\delta m = u_1\).
Problem 1: Write a function linear_sources(dm, u0s, config) that takes an arbitrary model perturbation dm and a time sequence of incident wavefields u0s and generates the linear source wavefields (the right-hand-sides of the linear modeling equations). Functions you have previously written should be useful.
Problem 2: Use your leap_frog function to solve for the linear forward wavefield \(u_1\) due to a perturbation \(\delta m\). Use this code to write a function linear_forward_operator(C0, dm, config) which returns a tuple containing the wavefields and the sampled data.
When sampling is accounted for, the adjoint condition requires that,
\[ \left< \mathbf{S}F\delta m, d \right> _{\mathcal{D}} = \left< \delta m, F^{*}\mathbf{S}^*d \right> _{\mathcal{M}}, \]
hold to high precision.
Problem 3: Verify that the adjoint condition is satisfied by your implementation of the linear forward modeling and migration operators. Be careful to take into account the differences in the model and data inner product. Write a function adjoint_condition(C0, config) that implements this test. How accurate is the relationship? It should be accurate to machine precision for random values of the data \(d\) and the model perturbation \(\delta m\). Be sure that you define \(\delta m(0) = \delta m(1) = 0\), as it is nonphysical to have sources on an absorbing boundary.
Consider modifying your construct_matrices function to accept a config key 'bc', which allows you to toggle between absorbing and homogeneous Dirichlet boundaries. This result should still hold.
# Setup the problem from before
config = dict()
# Configure the wavelet
config['nu0'] = 10 #Hz
# Physical domain parameters
config['x_limits'] = [0.0, 1.0]
config['nx'] = 201
config['dx'] = (config['x_limits'][1] - config['x_limits'][0]) / (config['nx']-1)
# Source position
config['x_s'] = 0.1
# Receiver position
config['x_r'] = 0.15
# Load the model
from msri_tools import model
C, C0 = model(config)
# Set CFL safety constant
config['alpha'] = 1.0/6.0
# Define time step parameters
config['T'] = 3 #seconds
config['dt'] = config['alpha'] * config['dx'] / C.max()
config['nt'] = int(config['T']/config['dt'])
from wave_imaging import adjoint_condition
print "Absorbing BC"
adjoint_condition(C0, config)
print "Dirichlet BC"
config['bc'] = 'dirichlet'
adjoint_condition(C0, config)
Full waveform inversion is the geophysics terminology for minimizing the objective
\[ J(m) = \frac{1}{2}\int_0^T||d_{x_r}(t) - \mathbf{S}\mathcal{F}[m]||_2^2 \textrm{d}t. \]
The term arises because we will be using the wave equation simulation to match the seismic data to 'invert' or solve the seismic inverse problem.
From the notes, we have that the gradient of the objective at a point \(m_0\) is
\[ \nabla J[m_0] = -F^*(d_{x_r} - \mathbf{S}\mathcal{F}[m_0] ) = -F^*r. \]
The method of gradient (or steepest) descent attempts to solve the optimization problem using the iteration
\[ m_{k+1} = m_k -\beta_k \nabla J[m_k], \]
where \(\beta_k\) is a step parameter, best found via a 'line search,' that is necessary because the gradient provides no notion of scale. (Note: most texts use \(\alpha\) rather than \(\beta\), but I want to avoid confusion with the \(\alpha\) in the CFL condition.) With a sufficiently sized \(\beta_k\), this method can recover the true modem \(m\) approximately.
Problem 4: Write a function gradient_descent(C0, d, k, config) which takes a starting model C0, 'measured' data d, and a number of iterations k as arguments and returns the sequence of k estimates of the true model, as well as the values of the objective function \(J\) at each of those points. How might you select \(\beta_k\)? How does \(m_k\) compare to \(m\)? Why? How many iterations do you need before you start to recover your model? The model recovery is not perfect. Explain the errors.
Warning!! Be very careful to update the model correctly, as \(\delta m \ne c - c_0 = \delta c\)! How does \(\delta m\) relate to \(\delta c\)?
Plot the convergence of your gradient descent method using a log-scaled plot.
Problem 5: Download the latest version of msri_tools.py from the website. Use your gradient_descent function to try to recover the new Gaussian model gauss_model. How do things change? Explain the errors in the recovery.
Bonus Problem 1: It is compuationally intractible to compute the linear forward operator \(F\) directly. Why? If you wanted to explicitely compute this operator, how would you do it?
Bonus Problem 2: Gradient descent is going to converge slowly. How might you accelerate this process?