THE DIFFERENTIAL EQUATION

Many situations in physics are expressed by the differential equation

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

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, and $\tilde s(z)$ (summed over frequency) is the derived image. 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}$$ (2)

For convenience, I introduce 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 )$.Besides crossing layers, we must cross layer boundaries where the (reflection) sources add in to the upcoming wave. Thus we have  

$$u_{k-1} \quad = \quad Z_k u_k + s_{k-1}$$ (3)

Recursive use of equation (3) 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$$ (4)

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

The conjugate of the delay operator Z is the time advance operator $\bar Z$.The conjugate operation to that in equation (4) is given by

 

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

The conjugacy of equation (4) and (5) seems obvious, but it is not quite the elementary form we are so familiar with because the matrix multiplies the output (instead of the customary matrix multiplying the input). To prove the conjugacy, notice that equation (4) is equivalent to $\bold u = {\bold M}^{-1}\bold s$ whose conjugate, 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 conjugate of an inverse is the inverse of the conjugate) which gives $\bold M' \tilde\bold s = \bold u$which is equation (5).

We observe the wavefield only on the surface z=0, so the conjugacy of equations (4) and (5) is academic because it relates the wavefield at all depths with the source at all depths. Thus we need to truncate $\bold u$ to its first coefficient u₀. That will change the conjugate relation. Remember that the conjugate to truncation is zero padding. Thus the conjugate to solving equation (4) and then abandoning the (unobservable) waves beneath the surface is the solution to

 

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

which amounts to 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}$$ (7)

A final feature of the migration application is that the image is formed from $\tilde \bold s$by summing over all frequencies. The proposed subroutine is

# Phase shift modeling and migration.   (Warning: destroys its input!)
#
subroutine gazconj( conj, dt,dx, v,nt,nx, modl,         data )
integer	            conj,          nt,nx,			 iw, ikx, iz,nz
complex	iktau,  	          	  modl(nt,nx),  data(nt,nx)
real	                  dt,dx, v(nt),                  pi, w,w0,dw, kx,kx0,dkx
temporary complex               cs(nt)
call conjnull(     conj,  0,              modl,2*nt*nx,  data,2*nt*nx )
pi = 4.*atan(1.);	w0  = -pi/dt;	dw = 2.*pi/(nt*dt)
nz = nt;		kx0 = -pi/dx;	dkx= 2.*pi/(nx*dx)
if( conj == 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, nt {       w  = w0  + (iw -1) * dw
	if( conj== 0) { data(iw,ikx) = modl(nz,ikx)
		do iz = nz-1, 1, -1 
			data(iw,ikx) = data(iw,ikx) * iktau(dt,w,v(iz)*kx,dw) +
				       modl(iz,ikx)
		}
	else {		cs(   1) =     data(iw,ikx)
		do iz = 1, nz-1
			cs(iz+1) = cs(iz) *   conjg(  iktau(dt,w,v(iz)*kx,dw) )
		do iz = 1, nz
			modl(iz,ikx) = modl(iz,ikx)  + cs(iz)
		}
	}}
if( conj == 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

complex function iktau( dt, w, vkx, nyep)
real                    dt, w, vkx, nyep 	# nyep = nyquist * epsilon
iktau = cexp( - dt * csqrt( cmplx( nyep, w) ** 2   +   vkx * vkx /4. ) )
return; end