Adjointness and ordinary differential equations

It is straightforward to adapt the simple programs phasemig() and phasemod() 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 (20) across a medium of three layers is expressed in matrix form as  

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

A recursive solution begins at the bottom with u₃=s₃ and propagates upward.

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

 

$$\bold M' \, \tilde {\bold s} \quad =\quad \left[ \begin{arr... ...u_1 \\ 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 (21) and (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 (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 (22).

We observe the wavefield only on the surface z=0, so the adjointness of equations (21) and (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₀ since the upcoming wave is known only at the surface. This truncation changes the adjoint in a curious way. We rewrite equation (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{arr... ...begin{array} {c} u_0 \\ 0 \\ 0 \\ 0 \end{array} \right] \\ end{displaymath}$$ (23)

The operator 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 and the code in subroutine gazadj() , I ran the dot product test to be sure!  

# Phase shift modeling and migration.   (Warning: destroys its input!)
#
subroutine gazadj( adj, dt,dx, v,nt,nx, modl,         data )
integer            adj,          nt,nx,                  iw, ikx, iz,nz
complex eiktau, cup,                    modl(nt,nx),  data(nt,nx)
real                      dt,dx, v(nt),              pi, w,w0,dw, kx,kx0,dkx,qi
call adjnull(     adj,  0,              modl,nt*nx*2, data,nt*nx*2 )
pi = 4.*atan(1.);       w0  = -pi/dt;   dw = 2.*pi/(nt*dt);    qi=.5/(nt*dt)
nz = nt;                kx0 = -pi/dx;   dkx= 2.*pi/(nx*dx)
if( adj == 0)   call ft2axis( 0, -1., nz, nx, modl)
else    {       call ft2axis( 0, -1., nt, nx, data)
                call ft1axis( 0,  1., nt, nx, data)
                }
do ikx  = 2, nx     {       kx = kx0 + (ikx-1) * dkx
do iw   = 2, 1+nt/2 {       w  = w0  + (iw -1) * dw
        if( adj== 0) { data(iw,ikx) = modl(nz,ikx)
                do iz = nz-1, 1, -1 
                        data(iw,ikx) = data(iw,ikx) * eiktau(dt,w,v(iz)*kx,qi) +
                                       modl(iz,ikx)
                }
        else {          cup          = data(iw,ikx)
                do iz = 1, nz {
                        modl(iz,ikx) = modl(iz,ikx)  +  cup
                        cup          = cup *  conjg(  eiktau(dt,w,v(iz)*kx,qi))
                        }
                }
        }}
if( adj == 0)  {  call ft1axis( 1,  1., nt, nx, data)
                  call ft2axis( 1, -1., nt, nx, data) }
else            { call ft2axis( 1, -1., nz, nx, modl) }
return; end

Finally, a few small details about the code. The loop on spatial frequency ikx begins at ikx=2. The reason for the 2, instead of a 1, is to omit the Nyquist frequency. If the Nyquist frequency were to be included, it should be divided into one half at positive Nyquist and one half at negative Nyquist, which would clutter the code without adding practical value. Another small detail is that the loop on temporal frequency iw begins at iw=1+nt/2 which effectly omits negative frequencies. This is purely an economy measure. Including the negative frequencies would assure that the final image be real, no imaginary part. Omitting negative frequencies simply gives an imaginary part that can be thrown away, and gives the same real image, scaled by a half. The factor of two speed up makes these tiny compromises well worthwhile.