THE STOLT METHOD

On most computers the Stolt [1978] method of migration is the fastest one--by a wide margin. For many applications, this will be its most important attribute. For a constant-velocity earth it incorporates the Huygens wave source exactly correctly. Like the other methods, this migration method can be reversed and made into a modeling program. One drawback, a matter of principle, is that the Stolt method does not handle depth variation in velocity. This drawback is largely offset in practice by an approximate correction that uses an axis-stretching procedure. A practical problem is the periodicity of all the Fourier transforms. In principle this is no problem at all, since it can be solved by adequately surrounding the data by zeroes.

A single line sketch of the Stolt method is this:

$$p(x,t) \ \ \rightarrow \ \ P(k_x , \omega ) \ \ \rightarrow ... ...sqrt{ \omega^2 / v^2 \ -\ k_x^2 } ) \ \ \rightarrow \ \ p'(x,z)$$

To see why this works, begin with the input-output relation for downward extrapolation of wavefields:  

$$P( \omega , k_x , z ) \ \ =\ \ e^{ i\,k_z z } \ P( \omega , k_x , z = 0 )$$ (22)

Perform a two-dimensional inverse Fourier transform:

$$p(t,x,z) \ \ =\ \ \int \int \ e^{ i\,k_x x - i \omega t + i\,k_z z } \ P( \omega , k_x , 0 ) \ d \omega \ dk_x$$ (23)

Apply the idea that the image at (x,z) is the exploding-reflector wave at time t = 0:  

$$\hbox{Image} (x,z) \ \ =\ \ \int \int \ e^{ i\,k_x x } \ e... ...\omega , k_x )\, z } \ P( \omega , k_x , 0 ) \ d \omega \ dk_x$$ (24)

Equation (24) gives the final image, but it is in an unattractive form, since it implies that a two-dimensional integration must be done for each and every z-level. The Stolt procedure converts the three-dimensional calculation thus implied by (24) to a single two-dimensional Fourier transform.

So far nothing has been done to specify an upcoming wave instead of a downgoing wave. The direction of the wave is defined by the relationship of z and t that is required to keep the phase constant in the expression exp$(-i \omega t \ +\ ik_z z )$.If $\omega$ were always positive, then +k_z would always refer to a downgoing wave and -k_z to an upcoming wave. Negative frequencies $\omega$ as well as positive frequencies are needed to describe waves that have real (not complex) values. So the proper description for a downgoing wave is that the signs of $\omega$ and k_z must be the same. The proper description for an upcoming wave is the reverse. With this clarification the integration variable in (24) will be changed from $\omega$ to k_z.

$$\begin{eqnarray} \omega \ \ & =& \ \ -\,\sgn (k_z ) \ v \sqrt { k_x^2 \ +\ k_z^... ... = &\ \ {-\,v \vert k_z \vert \over \sqrt { k_x^2 \ +\ k_z^2 } }\end{eqnarray}$$ (25) (26) (27)

Put (25), (26) and (27) into (24), and include also a minus sign so that the integration on k_z goes from minus infinity to plus infinity as was the integration on $\omega$. 

$$\hbox{Image}(x,z) \quad =\quad \int \int e^{ ik_z z + ik_x x... ...ert \over \sqrt { k_x^2 + k_z^2 } } \, \right\} \, dk_z \, dk_x$$ (28)

Equation (28) states the result as a two-dimensional inverse Fourier transform. The Stolt migration method is a direct implementation of (28). The steps of the algorithm are

• Double Fourier transform data from p(t,x,0) to $P( \omega , k_x , 0)$.
• Interpolate P onto a new mesh so that it is a function of k_x and k_z. Multiply P by the scale factor (which has the interpretation $\cos \, \theta$).
• Inverse Fourier transform to (x,z)-space.

Samples of Stolt migration of impulses are shown in Figure 9.

 

stolt
Figure 9 Response of Stolt method to data with impulses. Semicircles are seen, along with computation artifacts.

You can see the expected semicircular smiles. You can also see a semicircular frown hanging from the bottom of each semicircle. The worst frown is on the deepest spike. The semicircular mirrors have centers not only at the earth's surface z=0 but also at the bottom of the model $z=z_{\rm max}$.It is known that these frowns can be suppressed by interpolating more carefully. (Interpolation is the way you convert from a uniform mesh in $\omega$,to a uniform mesh in k_z). Interpolate with say a sinc function instead of a linear interpolator. A simpler alternative is to stay away from the bottom of the model, i.e. pad with lots of zeroes.

 

stolt4 Figure 10 Periodicity of the output of the Stolt migration program.

It seems that an extraordinary amount of zero padding is required on the time axis. To keep memory requirements reasonable, the algorithm can be reorganized as described in an exercise. Naturally, the periodicity in x also requires padding the x-axis with zeroes.

Below is a subroutine for Stolt migration. Its input is assumed to be data after 2-D Fourier transformation to $( \omega , k_x )$-space, and its output is the image in to (k_z, k_x)-space, so it needs a final transformation to (z,x)-space.  

subroutine stoltmig( nt,nz,nx, dt,dz,dx, vel, cftdata, cftimage)
integer it,iz,ix, nt,nz,nx, iw,ikz,ikx, nw,nkz,nkx
real    dt,dz,dx, w0,kz0,kx0, dw,dkz,dkx, pi, vel,vhalf, w,kz,kx, signum
complex ckzkx, cftdata(nt,nx), cftimage(nz,nx)
pi = 3.1415926
nw  = nt;       w0  = -pi/dt;   dw  = 2.*pi/(nt*dt);
nkz = nz;       kz0 = -pi/dz;   dkz = 2.*pi/(nz*dz);
nkx = nx;       kx0 = -pi/dx;   dkx = 2.*pi/(nx*dx);
do ikz= 1, nkz {        cftimage( ikz,   1) = (0.,0.) } # zero the Nyquist.
do ikx= 1, nkx {        cftimage(   1, ikx) = (0.,0.) } # zero the Nyquist.
vhalf = 0.5 * vel
do ikx = 2, nkx {       kx = kx0 + dkx * (ikx-1)
do ikz = 2, nkz {       kz = kz0 + dkz * (ikz-1)
        w = - signum( kz) * vhalf * sqrt( kx*kx + kz*kz)
        call cinterp1( w, nw,w0,dw, cftdata(1,ikx), ckzkx)
        if( w != 0.)
                cftimage(ikz,ikx) = ckzkx * abs( kz / w)
        else
                cftimage(ikz,ikx) = 0.
        }}
return;         end

stoltmig() uses the linear interpolation routine cinterp1() and the signum function $\sgn(x) = x/\vert x\vert$ done by routine signum() . In stoltmig() I have omitted the theoretically required $\cos\theta$ multiplication. This angle-dependent scaling can be left as an exercise. The effect of such an angle-dependent multiplier is to scale the signal amplitude along the hyperbola according to the local slope of the hyperbola.  

# Nearest neighbor interpolation would do this:    cbx = cb( 1.5 + (x-x0)/dx)
# Linear interpolation with the same definitions of variables is:
#
subroutine cinterp1( x, nx,x0,dx, cb, cbx)
integer ix,ixc,         nx
real    x, xc, x0, dx, fraction
complex cb(nx), cbx
xc      = (x-x0) / dx
ixc     =  xc
fraction = xc - ixc
ix      = 1 + ixc
if( ix < 1 )
                cbx = cb( 1)
else if( ix+1 > nx)
                cbx = cb( nx)
else
                cbx = (1.-fraction) * cb(ix)  +  fraction * cb(ix+1)
return; end

 

real function signum(x)
real x
        if      (x>0)  { signum =  1 }
        else if (x<0)  { signum = -1 }
        else           { signum =  0 }
return;end