A 1D Seismic Imaging Example

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.

Resources and Code Setup

These exercises are best solved using a single python script file, which I will call cp2_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). 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 cp2_sol

As 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 model

at the top of your python scripts.

Part 2: Data, the Residual and the Adjoint State Method

In this exercise, you will generate seismic data, solve the adjoint wave equation, and generate images of the reflectors.

In [1]:
# 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

# 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'])

Generating Data

A receiver (a seismometer or geophone) at spatial position \(x_r\) records the value of the true wavefield (or a function of it) at a point, and is written mathematically as

\[ d_r(t) = d(x_r,t) = \int_\Omega u(x,t) \delta(x - x_r)\textrm{d}x. \]

For 2D and 3D seismic imaging problems there are multiple receivers at different spatial positions recording data for a single 'shot' (instance of a source). This 'sampling' can be denoted with the operator \(\mathbf{S}\) and is written as

\[ \mathbf{d}(t) = \left[ \begin{array}{c} d_{r_1}(t) \\\\ d_{r_2}(t) \\\\ \vdots \\\\ d_{r_n}(t) \end{array}\right] = \left[ \begin{array}{c} \int_\Omega u(x,t) \delta(x - x_{r_1})\textrm{d}x \\\\ \int_\Omega u(x,t) \delta(x - x_{r_2})\textrm{d}x \\\\ \vdots \\\\ \int_\Omega u(x,t) \delta(x - x_{r_n})\textrm{d}x \end{array}\right] = \mathbf{S}u(x,t). \]

Your point_source function from Part 1 implements the adjoint operation of sampling, \(\mathbf{S^*}\), also known as interpolation or prolongation.

Problem 1: Write a function record_data(u, config), which takes a wavefield and a receiver position and returns the measured data, known as a trace.

In [2]:
# You will need to define the receiver position in your configuration

# Receiver position
config['x_r'] = 0.15

The Forward Operator

In Part 1 you wrote code to solve the forward problem,

\[ \mathcal{F}\left[ m \right] = u. \]

It will be useful to put the necessary steps into a function, as we will want to solve this problem many times, perhaps on different problems. Additionally, we will frequently want to solve the sampled problem,

\[ \mathbf{S}\mathcal{F}\left[ m \right] = d. \]

Problem 2: Write a Python function forward_operator(C, config) which returns a tuple containing the wavefields and the sampled data. Plot and label the trace. Use \(x_r = 0.15\). This function should utilize the functions you have written in the previous exercises.

In [3]:
from cp2_sol import forward_operator

us, d = forward_operator(C, config)
In [4]:
ts = linspace(0, config['T'], config['nt'], False)

figure()
plot(ts, d, label=r'$x_r =\,{0}$'.format(config['x_r']), linewidth=2)
xlabel(r'$t$', fontsize=18)
ylabel(r'$d(t)$', fontsize=18)
title('Trace at $x_r={0}$'.format(config['x_r']), fontsize=22)
legend();

Data Residuals

The incident wavefield \(u_0(x,t)\) is generated from the proposal model \(c_0(x)\). The difference between the observed or measured data \(d\) and the simulated data \(d_0\) due to \(u_0\) is called the residual \(r\), and

\[ r(t) = d(t) - d_0(t) = d(t) - (\mathbf{S}u)(t) = d(t) - \mathbf{S}\mathcal{F}[m_0] (t). \]

Problem 3: Generate the incident field \(u_0\) from the proposal model \(m_0 = \dfrac{1}{c_0^2}\). Compute and plot the residual trace. Be sure to use the same \(\Delta t\) as you did for computing \(u(x,t)\), otherwise you will have to use interpolation to compute the residual.

The Adjoint Wavefield

From the notes, we have that the computed model perturbation \(\delta m\), which we will call the reverse time migration image \(I_\text{RTM}\), is generated by applying the adjoint of the linearized forward model \(F\) to the residual,

\[ I_\text{RTM} = F^*r, \]

or, when sampling is taken into account,

\[ I_\text{RTM} = (\mathbf{S} F)^*r = F^{*} \mathbf{S}^{*} r = F^*r_{\text{ext}}, \]

where \(r_\text{ext}\) is the projection of the residual onto the full wavefield.

Also from the notes,

\[ F^*r_{\text{ext}} = -\int_0^T q(x,t) \partial_{tt}u_0(x,t) \textrm{d}t, \]

where \(q(x,t)\) is the solution of

\[ \begin{align*} (\frac{1}{c_0(x)^2}\partial_{tt}-\partial_{xx})^{*}q_t(x,t) & = r_\text{ext}(x,t), \\\\ (\frac{1}{c_0(0)}\partial_t-\partial_x)^{*}q(0,t) & = 0, \\\\ (\frac{1}{c_0(1)}\partial_t+\partial_x)^{*}q(1,t) & = 0, \\\\ u(x,t) & = 0 \quad\text{for}\quad t \ge T, \end{align*} \]

where the adjoint operation \(*\) is taken in time and collectively \(L_0^{*}q(x,t) = r_\text{ext}(x,t)\) (as opposed to \(Lu(x,t)=f(x,t)\) for the true wavefield). (Note that the argument \(r_\text{ext}\) can be replaced by any wavefield, \(F*\) is defined for any argument, not just the extended residual.)

Problem 4: Extend your residual to the whole space using your point_source function at the receiver location \(x_r\). Using your leap_frog function, solve for the adjoint field \(q(x,t)\). Plot the space-time diagram of \(q(x,t)\).

Hint: Be very careful, as the absorbing boundary conditions imply that \(L\ne L^*\). What does the adjoint operation mean in the time domain? How does this impact the time derivatives and your discretization of those derivatives?

The Imaging Condition

As mentioned before, the perturbation is

\[ I_\text{RTM} = F^*\mathbf{S}^*r = \int_0^T\ q(x,t) \partial_{tt}u_0(x,t) \textrm{d}t, \]

this inner product is known as an imaging condition and is the core of a gradient based optimization method. The image \(I_\text{RTM}\) is said to have been produced by Reverse Time Migration.

Problem 5: Write a Python function imaging_condition(qs, u0s, config) that computes \(F^*r_\text{ext}\), given the incident wavefield \(u_0\). You may want to write a separate function to compute the second time derivative of the incident field. You can assume that \(u_0 = 0\) for \(t > T\). Plot the resulting image of the reflection. Compare that image to \(\delta C = C - C_0\). How are the different, and why? How might these results change if the source and receiver are on opposite sides of the reflector? What if there are multiple receivers?

In [8]:
from cp2_sol import imaging_condition

I_rtm = imaging_condition(qs, u0s, config)
In [9]:
xs = np.arange(config['nx'])*config['dx']
dC = C-C0

figure()
subplot(2,1,1)
plot(xs, dC, label=r'$\delta C$')
legend()
subplot(2,1,2)
plot(xs, I_rtm, label=r'$I_\text{RTM}$' );

Problem 6: Visualizing the imaging condition. Plot the space-time diagram for \(\partial_{tt}u_0(x,t)\). Either print the space-time diagrams of \(q(x,t)\) and \(\partial_{tt}u_0(x,t)\) and overlay them or use matplotlib to overlay the two images. See where the adjoint field and the incident field intersect to find the location of the reflection.

Problem 7: Write a function adjoint_operator(C0, d, config) which implements the adjoint operation you developed in Problems 4 and 5 and returns the reverse time migration image. Consider optionally returning the adjoint field.