STOLT STRETCH THEORY

In order to simplify the references, I will begin with the textbook definitions of the Stolt migration method. The reader familiar with Stolt stretch theory can skip this section and go on to a new piece of theory in the next one.

Post-stack seismic migration is theoretically a two-stage process consisting of wavefield downward continuation in depth z based on the wave equation  

$${\partial^2 P \over \partial x^2} + {\partial^2 P \over \par... ...\, = \,{1 \over {v^2(x,z)}}\, {\partial^2 P \over \partial t^2}$$ (1)

and the imaging condition t=0 (here the velocity v is twice as small as the actual wave velocity). Stolt time migration performs both stages in one step, applying the frequency-domain operator  

$$\tilde{P_0}\left(k_x, \omega_0\right)= \tilde{P_v}\left(k_x,... ...omega_v \left(k,\omega_0\right)}\over{d\omega_0}}\right\vert\;,$$ (2)

where $P_0\left(x, t_0\right)$ stands for the initial zero-offset (stacked) seismic section defined on the surface z=0, $P_v\left(x, t_v\right)$ is the time-migrated section, and t_v is the vertical traveltime  

$$t_v=\int_{0}^{z}{dz \over {v(x,z)}}\;\;.$$ (3)

The function $\omega_v \left(k,\omega_0 \right)$ in (2) corresponds to the dispersion relation of the wave equation (1) and in the constant velocity case has the explicit expression  

$$\omega_v \left(k,\omega_0 \right)=\mbox{sign}\left(\omega_0\right) \sqrt{\omega_0^2 - v^2 k^2}\;\;.$$ (4)

The choice of the sign in (4) is essential to distinguish between upgoing and downgoing waves. It is the upgoing part of the wave field that is used in migration.

For the case of a varying velocity Stolt 1978 suggested the following change of the time variable (referred to in the literature as Stolt stretch):  

$$s(t)={\left({{2 \over v_0^2}\,\int_0^t\eta d \tau}\right)}^{1/2}\;,$$ (5)

where v₀ is an arbitrarily chosen constant velocity, and $\eta$ is a function defined by the parametric expressions

$$\eta(\zeta)=\int_0^{\zeta} v(x,z) \,dz \;,\; \tau(\zeta)=\int_0^{\zeta} { {dz} \over {v(x,z)}}\;\;.$$ (6)

With the stretch (5), seismic time migration can be related to the transformed wave equation  

$${\partial^2 P \over \partial x^2} + W\,{\partial^2 P \over \... ...2-W) \over v_0^2 }\, {\partial^2 P \over \partial \hat{t}^2}\;.$$ (7)

Here $\hat{z}$ and $\hat{t}$ are the transformed depth and time coordinates that possess the following property: if $\hat{z}=0$, $\hat{t}=s\left(t_0\right)$, and if $\hat{t}=0$, $\hat{z}=v_0 s\left(t_v\right)$. W is a varying coefficient defined as  

$$W=a^2+2b\,(1-a^2)\;,$$ (8)

where Stolt's idea was to replace the slowly varying parameter W with its average value. Thus equation (7) is approximated by an equation with constant coefficients, which has the dispersion relation  

$$\widehat{\omega}_v\left(k,\widehat{\omega}_0 \right)= \left(... ...\right)}\over W}\, \sqrt{\widehat{\omega}_0^2 - W v_0^2 k^2}\;.$$ (9)

Stolt's approximate method for migration in heterogeneous media consists of the following steps:

1.

stretching the time variable according to (5),

2.

interpolating the stretched time to a regular grid,

3.

double Fourier transform,

4.

f-k time migration by operator (2) with the dispersion relation (9),

5.

inverse Fourier transform,

6.

inverse stretching (shrinking) the vertical time variable on the migrated section.

The value of W must be chosen prior to migration. According to Stolt's original definition (8), the depth variable z gradually changes in the migration process from zero to $\zeta$, causing the coefficient b in (8) to change monotonically from 0 to 1. If the velocity v monotonically increases with depth, then $\eta''(z)={\partial v \over \partial z}\geq 0$, and the average value of b is  

$$\bar{b}={1 \over {\zeta \eta(\zeta)}}\, {\int_0^\zeta \eta(z... ...{\int_0^\zeta {\eta(\zeta) {z \over \zeta} }dz}= {1 \over 2}\;.$$ (10)

As follows from (8) and (10), in the case of monotonically increasing velocity, the average value of W has to be less than 1 (W equals 1 in a constant velocity case). Analogously, in the case of a monotonically decreasing velocity, W is always greater than 1. In practice, W is included in migration routines as a user-defined parameter, and its value is usually chosen to be somewhere in the range of 1/2 to 1.

In this paper I will describe a straightforward way to determine the most appropriate value of W for a given velocity distribution.

A useful tool for that purpose is Stewart Levin's formula for the traveltime curve. Levin 1985 applied the stationary phase technique to the dispersion relation (9) to obtain an explicit formula for the summation curve of the integral migration operator analogous to the Stolt stretch migration. The formula evaluates the summation path in the stretched coordinate system, as follows:  

$$s\left(t_0\right)= \left(1-{1\over W}\right) s\left(t_v\righ... ...ft(t_v\right) + {W\, {{\left(x-x_0\right)^2} \over v_0^2}}}\;.$$ (11)

Here x₀ is the midpoint location on a zero-offset seismic section, and x is the space coordinate on the migrated section. Formula (11) shows that, with the stretch of the time coordinate, the summation curve has the shape of a hyperbola with the apex at $\left\{x,s\left(t_v\right)\right\}$ and the center (the intersection of the asymptotes) at $\left\{x,{\left(1-{1\over W}\right)}\,s\left(t_v\right)\right\}$. In the case of homogeneous media, W=1, $s(t)\equiv t$, and (11) reduces to the well-known hyperbolic diffraction traveltime curve. It is interesting to note that inverting formula (11) for $s\left(t_v\right)$ determines the impulse response of the migration operator, which can be interpreted as the wavefront from a point source in the $\{x,\hat{z},\hat{t}\}$ domain of equation (7):  

$$\hat{z}-\hat{z_0}= \left({1\over Q}-1\right) R \pm {1\over Q}\, \sqrt{R^2 - {Q\,{\left(x-x_0\right)^2}}}\;,$$ (12)

where $R=v_0 \hat{t}$, and Q=2-W. According to equation (12), wavefronts from a point source in the stretched coordinates for W<2 have an elliptic shape, with the center of the ellipse at $\{x,\hat{z_0}+\left( {1 \over Q}-1\right)\, R \}$ and the semi-axes $a_x={R \over \sqrt{Q}}$ and $a_z={r \over Q}$. The ellipses stretch differently for W<1 and W>1 (Figure 1). In the upper part that corresponds to the upgoing waves, they look nearly spherical, since the radius of the front curvature at the top apex equals the distance from the source.

 

stofro Figure 1 Wavefronts from a point source in the stretched coordinate system. Left: velocity decreases with depth (W=1.5). Right: velocity increases with depth (W=0.5).