# Part 2: The Adjoint State Method¶

In the previous section, you implemented the nonlinear forward operator

In the notes, the linearized version of this operator is derived and denoted by the operator \(F\). In this section, you will implement the adjoint of this operator, \(F^*\), and see how the adjoint state method can be used to 1) produce ‘images’ from the imaging condition (Reverse Time Migration) and 2) compute the gradient of the full waveform inversion nonlinear least squares objective function.

The linear forward operator \(F\) is addressed in full in the next section.

## 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

Problem 2.1

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.

## Adjoint Wave Solver¶

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,

or, when sampling is taken into account,

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

Also from the notes,

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

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 2.2

Extend your residual to the whole space using your `point_source`

function at the receiver location \(x_r\), to create a list of
*adjoint source wavefields* . 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?

## Imaging Condition¶

As mentioned before, the perturbation is

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 2.3

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\).

Hint

Don’t forget to account for the measure in your integrals.

Plot the resulting image of the reflection. Compare that image to \(\delta C = C - C_0\). Think about how and why they are different? How might these results change if the source and receiver are on opposite sides of the reflector? What if there are multiple receivers?

```
# Compute the image
I_rtm = imaging_condition(qs, u0s, config)
# Plot the comparison
xs = np.arange(config['nx'])*config['dx']
dC = C-C0
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(xs, dC, label=r'$\delta C$')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(xs, I_rtm, label=r'$I_\text{RTM}$')
plt.legend()
```

Visualizing the imaging condition as a space-time diagram is a very good way to understand how the intersection of the incident and adjoint waves locates the reflectors.

Problem 2.4

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. Try to understand why the intersections appear where they do.

## Adjoint Operator¶

The adjoint operator is a linear operator and it will be useful to have a way
to apply it as a black box, much like we have a function `forward_operator`

for applying the nonlinear forward operator as a black box.

Problem 2.5

Write a function `adjoint_operator(C0, d, config)`

which implements the
adjoint operation you developed in Problems 2.2 and 2.3 and returns the
reverse time migration image. Consider optionally returning the adjoint
field.

Hint

Remember that the argument `d`

can be anything that ‘looks like’
data. For example, the residual.