Downward continuation

Adjointness and ordinary differential equations

It is straightforward to adapt the simple programs to depth variable velocity. As you might suspect, the two processes are adjoint to each other, and for reasons given at the end of chapter  it is desirable to code them to be so. With differential equations and their boundary conditions, the concept of adjoint is more subtle than previous examples. Thus, I postponed till here the development of adjoint code for phase-shift migration. This coding is a little strenuous, but it affords a review of many basic concepts, so we do so here. (Feel free to skip this section.) It is nice to have a high quality code for this fundamental operation.

Many situations in physics are expressed by the differential equation

$${du \over dz} - i \alpha u \quad =\quad s(z)$$ (18)

In the migration application, $u(z)$ is the up-coming wave, $\alpha=-\sqrt{\omega^2 /v^2 - k_x^2}$, $s(z)$ is the exploding-reflector source. We take the medium to be layered ($v$ constant in layers) so that $\alpha$ is constant in a layer, and we put the sources at the layer boundaries. Thus within a layer we have $du / dz - i \alpha u = 0$ which has the solution

$$u(z_k+\Delta z ) \quad =\quad u(z_k) e^{i \alpha \Delta z}$$ (19)

For convenience, we use the ``delay operator'' in the $k$-th layer $Z_k= e^{-i \alpha \Delta z}$ so the delay of upward propagation is expressed by $u(z_k) = Z_k u(z_k+\Delta z )$. (Since $\alpha$ is negative for upcoming waves, $Z_k= e^{-i \alpha \Delta z}$ has a positive exponent which represents delay.) Besides crossing layers, we must cross layer boundaries where the (reflection) sources add to the upcoming wave. Thus we have

$$u_{k-1} \quad =\quad Z_{k-1} u_k + s_{k-1}$$ (20)

Recursive use of equation (7.20) across a medium of three layers is expressed in matrix form as

$$\bold M \bold u \quad =\quad \left[ \begin{array}{cccc}... ..._1 \\ s_2 \\ s_3 \end{array} \right] \quad =\quad \bold s$$ (21)

A recursive solution begins at the bottom with $u_3=s_3$ and propagates upward.

The adjoint (complex conjugate) of the delay operator $Z$ is the time advance operator $\bar Z$. The adjoint of equation (7.21) is given by

$$\bold M' \tilde {\bold s} \quad =\quad \left[ \begin{ar... ...\\ u_2 \\ u_3 \end{array} \right] \ \quad =\quad \bold u$$ (22)

where $\tilde s(z)$ (summed over frequency) is the migrated image. The adjointness of equation (7.21) and (7.22) seems obvious, but it is not the elementary form we are familiar with because the matrix multiplies the output (instead of multiplying the usual input). To prove the adjointness, notice that equation (7.21) is equivalent to $\bold u = {\bold M}^{-1}\bold s$ whose adjoint, by definition, is $\tilde{\bold s} = ({\bold M}^{-1})' \bold u$ which is $\tilde{\bold s} = {(\bold M')}^{-1} \bold u$ (because of the basic mathematical fact that the adjoint of an inverse is the inverse of the adjoint) which gives $\bold M' \tilde\bold s = \bold u$ which is equation (7.22).

We observe the wavefield only on the surface $z = 0$, so the adjointness of equations (7.21) and (7.22) is academic because it relates the wavefield at all depths with the source at all depths. We need to truncate $\bold u$ to its first coefficient $u_0$ since the upcoming wave is known only at the surface. This truncation changes the adjoint in a curious way. We rewrite equation (7.21) using a truncation operator $\bold T$ that is the row matrix $\bold T = [1,0,0,\cdots ]$ getting $\bold u_0 = \bold T\bold u = \bold T \bold M^{-1} \bold s$. Its adjoint is $\hat {\bold s} = (\bold M^{-1})' \bold T' \bold u_0' = (\bold M')^{-1} \bold T' \bold u_0'$ or $\bold M' \hat {\bold s} = \bold T'\bold u_0$ which looks like

$$\bold M' \tilde {\bold s} \quad =\quad \left[ \begin{ar... ...in{array}{c} u_0 \\ 0 \\ 0 \\ 0 \end{array} \right] \$$ (23)

The operator 7.23 is a recursion beginning from $\tilde s_0 = u_0$ and continuing downward with

$$\tilde s_k \quad =\quad \bar Z_{k-1} \tilde s_{k-1}$$ (24)

A final feature of the migration application is that the image is formed from $\tilde \bold s$ by summing over all frequencies. Although I believe the mathematics above I ran the dot product test to be sure!

Downward continuation