Mathematical foundation of migration

Seismic Data Analysis

Series	Investigations in Geophysics
Author	Öz Yilmaz
DOI	http://dx.doi.org/10.1190/1.9781560801580
ISBN	ISBN 978-1-56080-094-1
Store	SEG Online Store

D.1 Wavefield extrapolation and migration

A fundamental equation of reflection seismology is the double-square-root (DSR) equation. This equation describes downward continuation of both shots and receivers into the earth. It is exact for all dips and offsets. Neglecting the velocity gradient dv(z)/dz makes the DSR equation also applicable to a stratified earth. The DSR equation can be extended, with some approximation, to treat weak lateral velocity variations. A comprehensive mathematical treatise of the DSR equation is found in Claerbout .

The basic 2-D theory for wavefield extrapolation is presented here. Then, using the DSR equation, a rigorous analysis of conventional seismic data processing is made. We show that conventional implementation of the DSR equation requires zero-dip and zero-offset assumptions.

Before discussing the DSR equation, we review the basic theory for wavefield extrapolation. Once the extrapolation equations are developed, they can be used with the imaging principle to migrate 2-D or 3-D prestack and poststack data.

Start with the 2-D scalar wave equation, which describes propagation of a compressional wavefield P(x, z, t) in a medium with constant material density and compressional wave velocity v(x, z):

$\left({\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial z^{2}}}-{\frac {1}{v^{2}}}{\frac {\partial ^{2}}{\partial t^{2}}}\right)P(x,z,t)=0,$

(D-1)

where x is the horizontal spatial axis, z is the depth axis (positive downward), and t is time. Given the upcoming seismic wavefield P(x, 0, t), which is recorded at the surface, we want to determine reflectivity P(x, z, 0). This requires extrapolating the surface wavefield to depth z, then collecting it at t = 0. The process of obtaining the earth’s reflectivity P(x, z, t = 0) from the observed wavefield P(x, z = 0, t) at the surface z = 0 is called migration, and the reverse process is called modeling (Figure D-1).

It is advantageous to decompose the wavefield into monochromatic plane waves with different angles of propagation from the vertical. Therefore, we will work in the Fourier transform domain whenever possible. The wavefield can always be Fourier transformed over time t. If there is no lateral velocity variation, then the wavefield also can be Fourier transformed over the horizontal axis x. Thus,

$P(k_{x},z,\omega )=\int \int P(x,z,t)\ \exp(ik_{x}x-i\omega t)dx\,dt,$

(D-2a)

and inversely,

$P(x,z,t)=\int \int P(k_{x},z,\omega )\ \exp(-ik_{x}x+i\omega t)dk_{x}\,d\omega .$

(D-2b)

When the differential operator in equation (D-1) is applied to equation (D-2b), we get

${\frac {\partial ^{2}}{\partial z^{2}}}P(k_{x},z,\omega )+\left({\frac {\omega ^{2}}{v^{2}}}-k_{x}^{2}\right)P(k_{x},z,\omega )=0.$

(D-3)

Although v can be varied with depth z in equation (D-3), for now we assume a constant velocity case. The stratified earth case is considered later in this appendix. Equation (D-3) has two solutions, one for upcoming waves, the other for downgoing waves. The upcoming wave solution to equation (D-3) is recognized as

$P(k_{x},z,\omega )=P(k_{x},z=0,\omega )\,\ \exp \left[-i{\sqrt {{\frac {\omega ^{2}}{v^{2}}}-k_{x}^{2}z}}\,\right].$

(D-4)

Figure D-1  Relationship between migration and wavefield modeling (see Section D.1).

Equation (D-4) also is the solution to the following one-way wave equation:

${\frac {\partial }{\partial z}}P(k_{x},z,\omega )=-i{\sqrt {{\frac {\omega ^{2}}{v^{2}}}-k_{x}^{2}}}P(k_{x},z,\omega ).$

(D-5)

This solution can be verified by substituting equation (D-4) into equation (D-5).

We define the vertical wavenumber as

$k_{z}={\frac {\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{\omega }}\right)^{2}}}.$

(D-6)

Equation (D-6) often is called the dispersion relation of the one-way scalar wave equation. By using this expression, equation (D-4) takes the simple form

$P(k_{x},z,\omega )=P(k_{x},z=0,\omega )\ \exp(-ik_{z}z).$

(D-7)

To determine the reflectivity P(x, z, 0) from the wavefield recorded at the earth’s surface P(x, 0, t), proceed as follows:

1. Perform a 2-D Fourier transform over x and t to get P(k_x, 0, ω).
2. Multiply by the all-pass filter exp(−ik_zz) to obtain the wavefield P(k_x, z, ω) at depth z.
3. Perform summation over ω to obtain P(k_x, z, 0).
4. Finally, inverse Fourier transform over k_x to obtain the earth’s image P(x, z, 0) at that depth.

For the constant velocity case P(k_x, k_z, 0) can be computed by a direct mapping in the transform domain from (k_x, ω) to (k_x, k_z) using equation (D-6) .

The main objective here is to interpret equation (D-7) as a tool for downward extrapolating wavefields given at the surface. While the mathematical development of the process presented is simple, its physical basis is not obvious. In an effort to develop a physical motivation for equation (D-7), a simpler derivation follows.

Given the upcoming wavefield P(x, 0, t) recorded at the surface, we can decompose it into monochromatic plane waves, each traveling at a different angle from the vertical. We identify these plane waves by attaching each one to a unique (k_x, ω) pair. This plane-wave decomposition is equivalent to Fourier transforming the wavefield to yield P(k_x, 0, ω).

Now consider one of these plane waves as shown in Figure D-2. Imagine that this plane wave passed point P at t = 0, traveled upward, and was recorded by a receiver at surface point G at time t. For reflector mapping, we need to take the energy located at point G on the wavefront at time t, back to its position at t = 0 — to reflection point P. To return the energy to P, it makes sense to follow the same raypath used for outward propagation from P. The fact that the same path is used means that downward continuation does not alter the horizontal wavenumber k_x.

Figure D-2  Geometry for wavefield extrapolation (see Section D.1) .

Suppose that the wavefront is moved to a depth Δz = GG′ beneath the receiver at G so that the waveform at G now is at G″. If a receiver were buried at G″, it would have recorded the plane wave at t − Δt, where Δt is the traveltime between G and G″. In other words, moving the receiver at G vertically down a distance Δz to a new location G′ changes the traveltime along the raypath from G to P by −Δt.

From the geometry of Figure D-2, we have

$\Delta t={\frac {\Delta z}{v}}\cos \theta ,$

(D-8a)

where v/cos θ is the vertical phase velocity. We know the k_x and ω values for the plane wave. Suppose the distance between G and G″ is one wavelength λ. At time t − Δt, the wavefront intersects the x-axis at distance λ_x from G. From the geometric relation in Figure D-2, we have

${\frac {\lambda }{\lambda _{x}}}=\sin \theta .$

(D-8b)

By using the definitions λ = 2π/(ω/v), λ_x = 2π/k_x, and equation (D-8b), we obtain

$\sin \theta ={\frac {vk_{x}}{\omega }}$

(D-8c)

and

$\cos \theta ={\sqrt {1-\left({\frac {vk_{x}}{\omega }}\right)^{2}}},$

(D-8d)

where ω/v is the wavenumber along the raypath. By substituting equation (D-8d) into equation (D-8a), we have

$\Delta t={\frac {1}{v}}{\sqrt {1-\left({\frac {vk_{x}}{\omega }}\right)^{2}}}\Delta z.$

(D-8e)

As we move down, we do not want to change the amplitude of the plane wave. Given the change in traveltime, Δt by equation (D-8e), the corresponding phase shift is −ωΔt. At each Δz step of descent, we may propagate the plane wave with a different velocity v(z). The total phase shift to which the waveform was subjected with arrival at P is − ∫ ω dt.

To compute the wavefield at P, we use equation (D-8e) and multiply the transformed surface wavefield P(k_x, 0, ω) by

$\exp \left(-i\int _{G}^{P}\omega \,dt\right)=\exp \left[-i\int _{0}^{z}{\frac {\omega }{v(z)}}{\sqrt {1-\left({\frac {v(z)k_{x}}{\omega }}\right)^{2}}}\,dz\right].$

(D-9)

Equation (D-9) is the same operator used in equation (D-7), except that equation (D-7) was derived for constant v.

We now return to the more mathematical discussion and consider the stratified earth with a velocity v(z). Since we have not Fourier transformed P(x, z, t) over z, the one-way wave equation (D-5) also is valid for v(z):

${\frac {\partial }{\partial z}}P(k_{x},z,\omega )=-i{\sqrt {{\frac {\omega ^{2}}{v^{2}(z)}}-k_{x}^{2}}}P(k_{x},z,\omega ),$

(D-10)

in which case equation (D-6) becomes

$k_{z}(z)={\frac {\omega }{v(z)}}{\sqrt {1-\left[{\frac {v(z)k_{x}}{\omega }}\right]^{2}}}.$

(D-11)

Substitution verifies that equation (D-10) has the following solution:

$P(k_{x},z,\omega )=P(k_{x},0,\omega )\ \exp \,\left[-i\int _{0}^{z}k_{z}(z)\,dz\right].$

(D-12)

We must check whether this solution satisfies the two-way scalar wave equation (D-3). By differentiating equation (D-10) and using equation (D-11), we have

${\frac {\partial ^{2}}{\partial z^{2}}}P=-i{\frac {dk_{z}(z)}{dv}}{\frac {dv(z)}{dz}}P-ik_{z}(z){\frac {\partial }{\partial z}}P,$

(D-13a)

where P = P(k_x, z, ω). By substituting equation (D-10) for ∂P/dz, we obtain

${\frac {\partial ^{2}}{\partial z^{2}}}P=-i{\frac {dk_{z}(z)}{dv}}{\frac {dv(z)}{dz}}P-k_{z}^{2}(z)P.$

(D-13b)

If the velocity gradient dv(z)/dz is ignored, then the first term on the right side drops out. The final expression then is

${\frac {\partial ^{2}}{\partial z^{2}}}P+k_{z}^{2}(z)P=0.$

(D-13c)

When equation (D-11) is substituted into this expression, we get:

${\frac {\partial ^{2}}{\partial z^{2}}}P+\left[{\frac {\omega ^{2}}{v^{2}(z)}}-k_{x}^{2}\right]P=0,$

(D-14)

which is identical to equation (D-3) where velocity can be varied with depth z.

So far, we have shown that a 2-D wavefield recorded at the earth’s surface can be extrapolated downward using the phase-shift operator given by equation (D-9). Wave extrapolation can be done through either a constant-velocity medium (equation D-7) or a vertically varying velocity medium (equation D-12). Seismic imaging is not complete until a stopping condition is imposed on downward continuation. The process of downward continuation is terminated when the clock, which measures t − ∫ dt, reads zero traveltime.

The concepts described above can be used to downward continue a complete seismic experiment that involves many shots and receivers. The vertical wavenumber given by equation (D-6) is expressed as a sum of two square roots, one associated with downward continuation of shots and the other associated with downward continuation of receivers:

$k_{z}={\frac {\omega }{v}}DSR(G,S),$

(D-15)

where v is the medium velocity that can be varied in depth z, and DSR stands for the double square root :

$DSR(G,S)={\sqrt {1-G^{2}}}+{\sqrt {1-S^{2}}},$

(D-16)

where G and S are the normalized receiver and shot wavenumbers, k_g and k_s, respectively:

$G={\frac {vk_{g}}{\omega }}$

(D-17a)

and

$S={\frac {vk_{s}}{\omega }}.$

(D-17b)

The newly defined vertical wavenumber in equation (D-15) is inserted into the extrapolation equation (D-7) to obtain

$P(k_{g},k_{s},z,\omega )=P(k_{g},k_{s},0,\omega )\ \exp(-ik_{z}z),$

(D-18)

where P(k_g, k_s, 0, ω) is the Fourier transform of the prestack data P(s, g, z = 0, t) in shot-receiver coordinates. This new extrapolation equation then can be used to downward continue common-shot gathers, and by way of reciprocity, common-receiver gathers.

The DSR equation (D-16) is separable in terms of shot and receiver wavenumbers. This separation means that we can start with the wavefields recorded at the surface as common-shot gathers and use the first part of the DSR operator to downward continue the receivers to depth Δz. Then, we can sort the already downward-continued wavefields into common-receiver gathers and use the second part of the DSR operator to downward continue the shots to depth Δz. By alternating between common-receiver and common-shot gathers, the entire seismic experiment (whole line) can be downward continued until imaging is accomplished (Figure D-3).

Although no approximation, besides the stratified earth assumption, is made in this alternating downward-continuation scheme, it is computationally exhausting. In fact, most of today’s seismic data processing is done in midpoint-(half) offset (y, h) coordinates, rather than in shot-receiver (s, g) coordinates. Therefore, we will put DSR as defined by equation (D-16) into the (y, h) coordinates.

The following coordinate transformation is required:

$y={\frac {1}{2}}(g+s)$

(D-19a)

and

$h={\frac {1}{2}}(g-s).$

(D-19b)

After the transformation , we obtain

$G=Y+H$

(D-20a)

and

$S=Y-H,$

(D-20b)

where Y and H are the normalized midpoint and offset wavenumbers, respectively:

$Y={\frac {vk_{y}}{2\omega }}$

(D-21a)

and

$H={\frac {vk_{h}}{2\omega }}.$

(D-21b)

Figure D-3  A flow diagram of shot-geophone migration (see Section D.1).

By substituting equations (D-21a) and (D-21b) into equation (D-16), the DSR equation takes the following form in midpoint-offset coordinates:

$DSR(Y,H)={\sqrt {1-(Y+H)^{2}}}+{\sqrt {1-(Y-H)^{2}}}.$

(D-22)

The vertical wavenumber (equation D-15) now is expressed in terms of normalized midpoint-offset wavenumbers Y and H:

$k_{z}={\frac {\omega }{v}}DSR(Y,H).$

(D-23)

The newly defined vertical wavenumber in equation (D-23) is inserted into the extrapolation equation (D-18):

$P(k_{y},k_{h},z,\omega )=P(k_{y},k_{h},0,\omega )\ \exp(-ik_{z}z),$

(D-24)

where P(k_y, k_h, 0, ω) is the Fourier transform of the prestack data P(y, h, z = 0, τ) in midpoint-offset coordinates.

Figure D-4 shows the ω − k_y plane for a specific value of z and h, and the k_y − z plane for a specific value of ω and h. The radial line A corresponds to k_y = 2ω/v. The region in the ω − k_y plane below the radial line corresponds to the evanescent energy and that above the radial line corresponds to the propagating energy. The same transition between the two regions also are noted in the k_y − z plane. The zero-offset case (Figure D-4c) clearly shows the evanescent energy to the right of the point on the k_y axis labeled as k_y = 2ω/v dying off rapidly with depth. The width of the propagation region stays constant with depth. The nonzero-offset case shown in Figure D-4d, however, indicates that the width of the propagation region varies with depth — zero at the surface z = 0 and approaching rapidly to the zero-offset case immediately at shallow depths. The physical interpretation of this depth-dependency is quite intuitive — the normalized offset wavenumber H becomes increasingly less significant at greater depths, and the normalized midpoint wavenumber Y becomes the dominating wavenumber. More specifically, on a CMP gather, moveout decreases with depth which implies nearly zero H.

Figure D-5 shows the response characteristics of the dispersion relation defined by equation (D-23). Note the semi-elliptical wavefronts in the y − z plane for a single frequency ω; while in the y − t plane, note the table-top traveltime trajectories. The equations for the wavefront and traveltime trajectories are derived in Section D.2 using stationary phase approximations.

Note that in equation (D-16), the terms with different spatial wavenumbers are separable. However, we have lost the property of separation in equation (D-22) because the operators in Y and H are strongly coupled. As a result, the Taylor series expansions of the square roots in equation (D-22) yield terms that contain cross-products of the two wavenumbers. The penalty for processing in the conventional coordinate system (y, h, t) is that strong coupling in the extrapolation operator requires the entire prestack data set to be handled at the same time for each depth step.

Figure D-4  Real part of the ω − k_y plane at z = 200 m: (a) DSR (equation D-23) with h = 0 m, and (b) DSR with h = 400 m. Real part of the k_y − z plane at ω = 16 m: (c) DSR (equation D-23) with h = 0 m, and (d) DSR with h = 400 m .

Conventional processing comprises two important steps. First, the data are organized into common-midpoint (CMP) gathers, and normal-moveout (NMO) correction is applied to each gather. The time shift Δt = t(h) − t₀ associated with the NMO correction is given by

$\Delta t=t_{0}\left[{\sqrt {1+\left({\frac {2h}{vt_{0}}}\right)^{2}}}-1\right],$

(D-25)

where t(h) is the two-way traveltime for a given (half) offset h, and t₀ is the corresponding two-way zero-offset time. Here, v is the root-mean-square (rms) velocity at t₀. Equation (D-25) is based on the stratified earth (zero-dip) assumption. After NMO correction, traces of the CMP gather are stacked. This not only reduces data volume, but also enhances the signal-to-noise ratio.

Figure D-5  The response characteristics of the DSR operator (equation D-23) . (a) Real part of the y − z plane at 16 Hz and h = 400 m. Note the semielliptical wavefronts. (b) Real part of the y − z plane at t = 1024 ms, h = 400 m. Because of the wraparound in h, we observe two wavefronts, one for h = 400 m and one for h = 0. (c) Real part of the y − t plane at z = 200, 400, 600, and 800 m superimposed. These are the table-top trajectories for h = 400 m. The loci of the arrival times are determined by a stationary-phase approximation to DSR (see Section D.2) . Periodicity in y and t result from approximating Fourier integrals by sums.

Second, the CMP stack is migrated as if it were the zero-offset wavefield generated by exploding reflectors (introduction to migration). The equation used for the downward extrapolation portion of migration is the solution to the one-way wave equation (equation D-12).

To account for the one-way traveltime of the exploding reflectors model, the velocity used in extrapolation is taken as half the medium velocity. Thus, the vertical wavenumber given by equation (D-11) is expressed as

$k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{y}}{2\omega }}\right)^{2}}},$

(D-26)

where v can be varied with depth z. (Since migration is done in midpoint space, k_x has been replaced with k_y.)

Some of the current migration techniques based on wave extrapolation use certain rational approximations to equation (D-26), while some implement the exact form in the frequency-wavenumber domain (migration principles).

Conventional processing theory has an advantage over the exact theory represented by the DSR equation in midpoint-offset space (equation D-22). Unlike the approach based on the DSR equation, the conventional approach is composed of two separable operators — the NMO correction and stack applied in offset space, and migration applied in midpoint space. However, this advantage is based on zero-dip and zero-offset assumptions.

Where do we go from here? On the one hand, we have an exact theory that can handle all dips and offset angles, but is difficult to implement. On the other hand, we have a conventional approach that has the convenient property of separation, but is based on the zero-dip and zero-offset assumptions.

To examine the relationship between the two approaches, we return to the exact theory and make the same two assumptions that underlie the conventional approach. The zero-dip assumption implies that the earth model is stratified in y − t domain. The seismic energy recorded over such an earth is concentrated completely at the zero midpoint wavenumber k_y = 0. This suggests that we set the normalized wavenumber Y equal to zero in DSR as defined by equation (D-22). The resulting operator here is defined as the stacking (St) operator

$St(H)=2{\sqrt {1-H^{2}}}.$

(D-27)

As shown in Section D.2, the NMO shift given by equation (D-25) is a stationary-phase approximation to equation (D-27) . The St(H) operator condenses primary information on a CMP gather down to zero offset. After applying this operator on a CMP gather, we may keep the zero-offset trace and abandon all other offsets. Since a CMP stack can be regarded as a zero-offset wavefield, equation (D-27) is a zero-dip NMO and stack-type operator.

Application of the zero-offset (h = 0) assumption into the DSR operator is more subtle. On a CMP gather at and near h = 0, energy essentially is concentrated at zero value of the offset wavenumber k_h = 0. In fact, NMO correction tries to push the primary energy on a CMP gather toward k_h = 0. Therefore, by setting the normalized offset wavenumber H = 0 in equation (D-22), the exploding reflectors (ER) migration operator can be expressed as

$ER(Y)=2{\sqrt {1-Y^{2}}}.$

(D-28)

Equation (D-28) is an approximation, because setting H = 0 is not exactly the same as setting h = 0. By setting H = 0 in equation (D-22) and using the dispersion relation given by equation (D-23), we have the zero-offset vertical wavenumber

$k_{z}={\frac {2\omega }{v}}{\sqrt {1-Y^{2}}}.$

(D-29)

By substituting the definition for Y from equation (D-21a) into equation (D-29), we obtain

$k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{y}}{2\omega }}\right)^{2}}},$

(D-30)

which is identical to equation (D-26). We conclude that the zero-offset migration operator ER(Y) (equation D-28) derived from the DSR equation (D-22) is identical to the migration operator that is based on the exploding reflectors model of conventional processing.

Figure D-6  The response characteristics of the exploding reflectors operator ER(Y) (equation D-29) . (a) Real part of the y − z plane at 16 Hz. Note the circular wavefronts. (b) Real part of the y − z plane at t = 1024 ms. (c) Real part of the y − t plane at z = 200, 400, 600, and 800 m superimposed. These are the hyperbolic trajectories. The loci of the arrival times are determined via the stationary-phase approximation to ER(Y) (see section D.2) . Periodicity in y and t result from approximating Fourier integrals by sums.

Figure D-6 shows the response characteristics of the dispersion relation defined by equation (D-30). Note the semicircular wavefronts in the y − z plane and the hyperbolic traveltime curves in the y − t plane. Compare these with the response of the complete DSR operator for the nonzero-offset case in Figure D-5. The equations for the wavefront and traveltime trajectories are derived in Section D.2 using stationary phase approximations.

D.2 Stationary phase approximations

In this section, we shall use the method of stationary phase to derive the traveltime equation inferred by the double square root equation for nonzero-offset source-receiver separation. Consider the double square root operator defined by equation (D-22) with the wavenumbers Y and H defined by equations (D-21a,21b). We want to operate on the transformed wavefield P(k_y, k_h, z = 0, ω) with the DSR operator. Subsequent inverse Fourier transformation will yield the wavefield P(y, h, z, t):

$P(y,h,z,t)=\int \int \int P(k_{y},k_{h},z=0,\omega )\ \exp(i\Phi z)\,dk_{y}\,dk_{h}\,d\omega ,$

(D-31)

where the total phase ϕ, normalized with respect to z, is given by

$\Phi =-{\frac {\omega }{v}}DSR(Y,H)-k_{y}{\frac {y}{z}}-k_{h}{\frac {h}{z}}+\omega {\frac {t}{z}}.$

(D-32)

The main contribution to integration in equation (D-31) occurs when the phase stays nearly constant. We therefore determine the variation of the phase with respect to variables k_y, k_h, and ω

${\frac {\partial \Phi }{\partial k_{y}}}=-{\frac {\omega }{v}}{\frac {\partial DSR}{\partial Y}}{\frac {\partial Y}{\partial k_{y}}}-{\frac {y}{z}},$

(D-33a)

${\frac {\partial \Phi }{\partial k_{h}}}=-{\frac {\omega }{v}}{\frac {\partial DSR}{\partial H}}{\frac {\partial H}{\partial k_{h}}}-{\frac {h}{z}},$

(D-33b)

and

${\frac {\partial \Phi }{\partial \omega }}=-{\frac {1}{v}}DSR-{\frac {\omega }{v}}\left[{\frac {\partial DSR}{\partial H}}{\frac {\partial H}{\partial \omega }}+{\frac {\partial DSR}{\partial Y}}{\frac {\partial Y}{\partial \omega }}\right]+{\frac {t}{z}},$

(D-33c)

and set each variation to zero. Substitute equation (D-32) and carry out the differentiations in equations (D-33a,33b,33c) to obtain

${\frac {1}{2}}{\frac {G}{\sqrt {1-G^{2}}}}+{\frac {1}{2}}{\frac {S}{\sqrt {1-S^{2}}}}={\frac {y}{z}},$

(D-34a)

${\frac {1}{2}}{\frac {G}{\sqrt {1-G^{2}}}}-{\frac {1}{2}}{\frac {S}{\sqrt {1-S^{2}}}}={\frac {h}{z}},$

(D-34b)

and

${\frac {1}{\sqrt {1-G^{2}}}}+{\frac {1}{\sqrt {1-S^{2}}}}={\frac {vt}{z}},$

(D-34c)

where G and S are defined by equations (D-20a,D-20b).

Now, eliminate G and S amongst equations (D-34a,D-34b,D-34c) to get the final expression from stationary phase approximation to the double square root equation as

${\sqrt {(y+h)^{2}+z^{2}}}+{\sqrt {(y-h)^{2}+z^{2}}}=vt.$

(D-35)

This is the equation of an ellipse in the y − z plane at constant t (Section E.5). Figure D-5 shows the elliptic wavefront and the table-top traveltime trajectory described by equation (D-35).

When equation (D-35) is specialized to the zero-offset case, h = 0, we obtain

${\sqrt {y^{2}+z^{2}}}={\frac {vt}{2}},$

(D-36)

which is a circle in the y − z plane at constant t and a hyperbola in the y − t plane at constant z. Figure D-6 shows the circular wavefront and the hyperbolic traveltime trajectory described by equation (D-36).

We now consider the stacking operator defined by equation (D-27). The total phase is given by

$\Phi =-{\frac {\omega }{v}}St(H)-k_{h}{\frac {h}{z}}+\omega {\frac {t}{z}}.$

(D-37)

Differentiate equation (D-37) with respect to k_h and ω and set the results to zero to obtain

${\frac {H}{\sqrt {1-H^{2}}}}={\frac {h}{z}}$

(D-38a)

and

${\frac {1}{\sqrt {1-H^{2}}}}={\frac {vt}{2z}}.$

(D-38b)

Now, eliminate H between equations (D-38a) and (D-38b) to obtain the stationary phase approximation to the stacking operator:

${\sqrt {h^{2}+z^{2}}}={\frac {vt}{2}}.$

(D-39a)

Define the zero-offset time as t₀ = 2z/v and substitute into equation (D-39a) to get

${\sqrt {h^{2}+\left({\frac {vt_{0}}{2}}\right)^{2}}}={\frac {vt}{2}}.$

(D-39b)

Finally, rearrange to get the equation for normal moveout:

$\Delta t_{NMO}=t_{0}\left[{\sqrt {1+\left({\frac {2h}{vt_{0}}}\right)^{2}}}-1\right],$

(D-40)

where Δt_NMO = t − t₀. This is the same equation as equation (3-2b) in the main text with offset defined as x = 2h.

D.3 The parabolic approximation

Start with the dispersion relation defined by equation (D-6) recast for the exploding reflectors model for which v is replaced with v/2 to obtain

$k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}}}.$

(D-41)

Then apply Taylor series expansion and retain the first two terms:

$k_{z}={\frac {2\omega }{v}}\left[1-{\frac {1}{2}}\left({\frac {vk_{x}}{2\omega }}\right)^{2}\right].$

(D-42a)

By simplifying, we get the dispersion relation associated with the parabolic equation

$k_{z}={\frac {2\omega }{v}}-{\frac {vk_{x}^{2}}{4\omega }}.$

(D-42b)

By operating on the wavefield P(x, z, t) and replacing −ik_zP with ∂P/∂z, we write the corresponding differential equation as

${\frac {\partial P}{\partial z}}=-i\left({\frac {2\omega }{v}}-{\frac {vk_{x}^{2}}{4\omega }}\right)P.$

(D-43)

Derivation of equation (D-43) is based on the constant-velocity assumption. Nevertheless, just as we did for the 90-degree one-way wave equation (D-10), the 15-degree one-way wave-equation (D-43) can be recast using a vertically varying velocity function v(z). Going one step further, once equation (D-43) is inverse Fourier transformed from the horizontal wavenumber k_x to the horizontal axis x, we will replace v(z) with a laterally varying velocity function v(x, z). Theoretically, this may not be permissible, but in practice, its validity is widely accepted.

The effect of translation is removed by retardation (Figure 4.1-18):

$\tau =2\int _{0}^{z}{\frac {dz}{{\bar {v}}(z)}},$

(D-44)

where ${\bar {v}}(z)$ generally is chosen to be a horizontal average of v(x, z). The time shift defined by equation (D-44) is equivalent to a phase shift in the frequency domain. Therefore, the actual wavefield P is related to the time-shifted wavefield Q by

$P=Q\ \exp(-i\omega \tau ).$

(D-45)

By differentiating with respect to z, we get

${\frac {\partial P}{\partial z}}=\left({\frac {\partial }{\partial z}}-i{\frac {2\omega }{{\bar {v}}(z)}}\right)Q\ \exp(-i\omega \tau ).$

(D-46)

Finally, by substituting equations (D-45) and (D-46) into equation (D-43), we obtain

${\frac {\partial Q}{\partial z}}=i{\frac {vk_{x}^{2}}{4\omega }}Q+i2\omega \left[{\frac {1}{{\bar {v}}(z)}}-{\frac {1}{v}}\right]Q.$

(D-47a)

The first term on the right side of this equation is called the diffraction term and the second term is called the thin-lens term.

Consider the special case of $v={\bar {v}}(z).$ The thin-lens term then vanishes and we are left with

${\frac {\partial Q}{\partial z}}=i{\frac {vk_{x}^{2}}{4\omega }}Q.$

(D-47b)

After inverse Fourier transforming, we obtain the parabolic differential equation

${\frac {\partial ^{2}Q}{\partial z\partial t}}={\frac {v}{4}}{\frac {\partial ^{2}Q}{\partial x^{2}}},$

(D-48)

where, in practice, velocity can be varied horizontally as well as vertically.

In some practical implementations of equation (D-48), downward continuation is performed in τ rather than z. The migrated section then is displayed in time τ and the process is called time migration. The two variables z and τ are related by equation (D-44). The dispersion relation in equation (D-41) for the scalar wave equation in terms of ω_τ = vk_z/2, the Fourier dual of τ, takes the form

$\omega _{\tau }=\omega {\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}}}.$

(D-49a)

By squaring both sides, we get the equation of an ellipse in the ω_τ − k_x plane.

The dispersion relation in equation (D-42b) for the parabolic equation, also expressed in terms of ω_τ, takes the form

$\omega _{\tau }=\omega -{\frac {v^{2}k_{x}^{2}}{8\omega }},$

(D-49b)

which is the equation of a parabola in the ω_τ − k_x, plane. The first term is associated with a vertical time shift that can be removed by retardation. To obtain the differential equation in terms of τ, equivalent to equation (D-48), apply the second term on the right side of equation (D-49b) to the retarded wavefield Q as defined by equation (D-45) and inverse Fourier transform

${\frac {\partial ^{2}Q}{\partial \tau \partial t}}={\frac {v^{2}}{8}}{\frac {\partial ^{2}Q}{\partial x^{2}}}.$

(D-50)

This equation is the basis for the 15-degree time migration algorithms.

The dispersion relations given by equations (D-49a) and (D-49b) are plotted in Figure D-7 for constant velocity v and input frequency ω = FE. Curve 1 is associated with the 90-degree scalar wave equation and is the same as shown in Figure 4.1-25, except that the latter is in terms of k_z. Point A is mapped onto point B after migration with the 90-degree equation (D-49a). The same point A is mapped onto point C after migration with the 15-degree equation (D-49b). The dispersion curve 2 for the parabolic equation increasingly departs from the dispersion curve 1 for the exact wave equation as the dip gets steeper. (The dip before migration is measured as the angle between the vertical axis FE and radial direction FA.) Thus, the parabolic equation causes more and more undermigration as dip increases.

D.4 Frequency-space implicit schemes

Higher-order approximations to the dispersion relation (equation D-41) of the one-way equation can be found by continued fraction expansion . When equation (D-41) is rewritten, we have

$k_{z}={\frac {2\omega }{v}}R,$

(D-51a)

where

$R={\sqrt {1-X^{2}}},$

(D-51b)

with

$X={\frac {vk_{x}}{2\omega }}.$

(D-51c)

The various orders of approximations to equation (D-51b) are defined by the following recurrence relation

$R_{n+1}=1-{\frac {X^{2}}{1+R_{n}}},$

(D-52)

with the initial value R₀ = 1. By setting n = 0 in equation (D-52), we have

$R_{1}=1-{\frac {X^{2}}{2}}.$

(D-53)

Then, substitution of equation (D-53) into equation (D-51a) yields

$k_{z}={\frac {2\omega }{v}}\left(1-{\frac {X^{2}}{2}}\right),$

(D-54)

which is the same equation obtained with the parabolic approximation, equation (D-42a).

The next higher-order expansion is obtained by setting n = 1 in equation (D-52) and using equation (D-53) to obtain

$R_{2}=1-{\frac {X^{2}}{2-{\frac {X^{2}}{2}}}},$

(D-55)

This is referred to as the 45-degree approximation to equation (D-41).

Ma discovered that the recurrence relation for the continuous fraction expansion can be expressed as ratios of two polynomials for the even-ordered expansions. He also showed that the expression can be split into the following partial fractions:

$R_{2n}=1-\sum \limits _{i=1}^{n}{\frac {\alpha _{i}X^{2}}{1-\beta _{i}X^{2}}}.$

(D-56)

For example, when n =1,

$R_{2}=1-{\frac {\alpha _{1}X^{2}}{1-\beta _{1}X^{2}}},$

(D-57)

which is equivalent to the 45-degree expansion in equation (D-55) when α₁ = 0.5 and β₁ = 0.25. Note that expansions R₄, R₆, R₈, … are made up of sums of the 45-degree term, each with a different set of coefficients, α_i and β_i. Lee and Suh minimized the difference in the least-squares sense between R of equation (D-51a) and R_2n of equation (D-56) for a specified dip angle and derived optimal coefficients (α_i, β_i) for up to the 10th order (Table D-1).

Figure D-7  Dispersion relations for the 90-degree equation (D-49a), 15-degree equation (D-49b) and the 45-degree equation (D-60a) with α₁ = 0.5 and β₁ = 0.25, plotted on the ω_τ − k_x plane for a given input frequency, where ω = FE and velocity is v. Input frequency ω = FE = AP is mapped on to output frequency ω_τ which is PB, PC, and PD for the 90-degree, 15-degree, and 45-degree equations, respectively. The wavenumber at the intersection of the curves along the k_x axis is $2\omega \!/\!v,\ 2{\sqrt {2}}\omega \!/\!v,\ {\text{and}}\ \left(4\omega \right)/\left({\sqrt {3}}v\right)$ for the 90-degree, 15-degree, and 45-degree equations, respectively.

Table D-1. Coefficients of optimized, fractioned, one-way wave equations .

Order, 2n	Accuracy Deg.	αi	βi
2	45	0.5	0.25
2	65	0.478242060	0.376369527
4	80	0.040315157	0.873981642
		0.457289566	0.222691983
6	87	0.004210420	0.972926132
		0.081312882	0.744418059
		0.414236605	0.150843924
8	90-	0.000523275	0.994065088
		0.014853510	0.919432661
		0.117592008	0.614520676
		0.367013245	0.105756624
8	90	0.000153427	0.997370236
		0.004172967	0.964827992
		0.033860918	0.824918565
		0.143798076	0.483340757
		0.318013812	0.073588213

We now derive the differential equation associated with the 45-degree dispersion relation. By substituting equation (D-57) and the definition for X given by equation (51c) into equation (D-51a), we obtain

$k_{z}={\frac {2\omega }{v}}\left[1-{\frac {\alpha _{1}v^{2}k_{x}^{2}}{4\omega ^{2}}}{\frac {1}{1-{\frac {\beta _{1}v^{2}k_{x}^{2}}{4\omega ^{2}}}}}\right].$

(D-58)

The first term is a vertical shift that can be removed by retardation in the same manner as for the 15-degree approximation (equations D-42 through D-46). By simplifying the remaining terms of equation (D-58), we obtain

${\frac {\beta _{1}}{\alpha _{1}}}{\frac {v}{2\omega }}k_{z}k_{x}^{2}-k_{x}^{2}-{\frac {1}{\alpha _{1}}}{\frac {2\omega }{v}}k_{z}=0.$

(D-59a)

Finally, by operating on the retarded wavefield Q(x, ω, z) (equation D-45), we obtain

$i{\frac {\beta _{1}v}{2\alpha _{1}\omega }}{\frac {\partial ^{3}Q}{\partial z\partial x^{2}}}-{\frac {\partial ^{2}Q}{\partial x^{2}}}+i{\frac {2\omega }{\alpha _{1}v}}{\frac {\partial Q}{\partial z}}=0.$

(D-59b)

Kjartansson implemented equation (D-59b) for 45-degree modeling and migration in the frequency-space domain. Migration in the frequency-space domain (commonly known as the ω − x algorithm) involves two interleaved operations:

1. a time shift based on equation (D-45), which is velocity-independent for time migration and velocity-dependent for depth migration, and
2. focusing the diffraction energy using equation (D-59b).

Once you have a code for the basic 45-degree operator, it is easy to implement the higher-order approximations that are given by equation (D-56) with the associated coefficients in Table D-1. Note that the difference between the 45-degree and 65-degree algorithms is the values used for coefficients (α₁, β₁). Also note that the 15-degree equation (D-47b) is obtained from equation (D-59b) by setting α₁ = 0.5 and β₁ = 0.

The dispersion relation in equation (D-58) also can be expressed in terms of ω_τ = vk_z/2, the Fourier dual of the time variable τ associated with the migrated data:

$\omega _{\tau }=\omega -{\frac {\alpha _{1}\omega v^{2}k_{x}^{2}}{4\omega ^{2}-\beta _{1}v^{2}k_{x}^{2}}}.$

(D-60a)

The first term is associated with a vertical time shift that can be removed by retardation as for the 15-degree equation (D-49b).

To obtain the differential equation in terms of τ, equivalent to equation (D-59b), apply the second term on the right side of equation (D-60a) to the retarded wavefield Q as defined by equation (D-45) and inverse Fourier transform:

$i{\frac {\beta _{1}}{\alpha _{1}\omega }}{\frac {\partial ^{3}Q}{\partial \tau \partial x^{2}}}-{\frac {\partial ^{2}Q}{\partial x^{2}}}+i{\frac {4\omega }{\alpha _{1}v^{2}}}{\frac {\partial Q}{\partial \tau }}=0.$

(D-60b)

This equation is the basis for the 45-degree and related steep-dip implicit finite-difference frequency-space time migration algorithms.

The dispersion relation given by equation (D-60a) for the 45-degree equation, in which α₁ = 0.5 and β₁ = 0.25, is plotted in Figure D-7 for constant velocity v and input frequency ω = FE. Point A is mapped onto point D after migration with the 45-degree equation (D-60a). The dispersion curve 3 for the 45-degree equation lies somewhere between those of the 90-degree equation (D-49a) (curve 1) and the 15-degree equation (D-49b) (curve 2).

Figure 4.4-1 shows the impulse responses of various approximations to the one-way dispersion relation based on equation (D-56). Note that the wavefronts become increasingly closer to a semicircle as the higher-order terms are included in equation (D-56). The 15-degree approximation yields an elliptical wavefront. The 45-degree approximation yields an impulse response of a heart shape. Refer to frequency-space migration in practice for the practical aspects of 2-D frequency-space steep-dip time migration, and 3-D poststack migration for its application to 3-D migration.

D.5 Stable explicit extrapolation

The exact extrapolation filter for a specific frequency ω and velocity v is expressed in the frequency-wavenumber domain as

$D(k_{x})=\exp \left\{-i{\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}}}\Delta z\right\}.$

(D-61)

The objective is to find, again for a specific frequency ω and velocity v, an explicit filter with complex coefficients h(x) in the frequency-space domain such that, when Fourier transformed to the frequency-wavenumber domain, the difference between the actual transform H(k_x) and the desired transform D(k_x) of equation (D-61) is minimum. To ensure stability, we impose the constraint that the amplitude of H(k_x) is never greater than unity within the propagation region k_x ≤ (2ω/v).

Consider the Fourier transform H(k_x) of the actual extrapolation filter h(x), which we will write in discrete form as h_n:

$H(k_{x})=h_{0}+2\sum \limits _{n=1}^{N}h_{n}\cos(nk_{x}),$

(D-62)

where 2N + 1 is the length of the symmetric filter h_n.

To determine the filter coefficients h_n, perform Taylor series expansion of the exact D(k_x) and the actual H(k_x) extrapolators given by equations (D-61) and (D-62), respectively, and match the coefficients of the terms in each series at k_x = 0. Such a direct match of the coefficients in the Taylor series results in a response H(k_x) whose amplitude is greater than unity beyond a certain k_x, thus violating the stability constraint that the amplitude of H(k_x) is never greater than unity within the propagation region k_x ≤ (2ω/v).

A way to circumvent the unstable response of the conventional Taylor series method is to match the first M < N coefficients in the series expansion and set the remaining N − M to zero . This modified Taylor series method proceeds as follows.

Let the filter coefficients h_n be defined by the series

$h_{n}=c_{0}+2\sum \limits _{m=1}^{M}c_{m}b_{mn},$

(D-63)

where

$b_{mn}=2\cos \left({\frac {2\pi mn}{N}}\right).$

(D-64)

By way of equations (D-63) and (D-64), equation (D-62) takes the form

$H(k_{x})=c_{0}\left[1+2\sum \limits _{n}\cos(nk_{x})\right]+2\sum \limits _{m}c_{m}\left[1+2\sum \limits _{n}\cos \left({\frac {2\pi mn}{N}}\right)\cos(nk_{x})\right].$

(D-65)

Equation (D-65) can be written as a single summation

$H(k_{x})=\sum \limits _{m=0}^{M}c_{m}B_{m}(k_{x}),$

(D-66)

where

$B_{0}(k_{x})=1+2\sum \limits _{n}\cos(nk_{x}),$

(D-67a)

and

$B_{m}(k_{x})=1+2\sum \limits _{n}\cos \left({\frac {2\pi mn}{N}}\right)\cos(nk_{x}).$

(D-67b)

Now, perform the Taylor series expansion of the exact extrapolator given by equation (D-61):

$D(k_{x})=D(0)+k_{x}{\frac {\partial D(k_{x})}{\partial k_{x}}}+{\frac {k_{x}^{2}}{2}}{\frac {\partial ^{2}D(k_{x})}{\partial k_{x}^{2}}}+\ldots ,$

(D-68)

where the derivatives are evaluated at k_x = 0, and thus the terms associated with the odd derivatives vanish.

Matching the terms of the actual extrapolator H(k_x) with those of the desired extrapolator D(k_x) is equivalent to matching their even derivatives

$\sum \limits _{m=0}^{M}c_{m}B_{m}^{2j}(0)=D^{2j}(0),j=0,1,2,\ldots ,M,$

(D-69)

where the derivative terms $B_{m}^{2j}$ are

$B_{0}^{2j}(0)=(-1)^{j}\left[1+2\sum \limits _{n=1}n^{2j}\right],$

(D-70a)

and

$B_{m}^{2j}(0)=2(-1)^{j}\left[1+2\sum \limits _{n=1}\cos \left({\frac {2\pi mn}{N}}\right)n^{2j}\right].$

(D-70b)

Equation (D-69) represents a system of M linear equations that can be solved for the coefficients c_m. The extrapolation filter coefficients h_n then can be computed by substituting the solution of equation (D-69) into equation (D-63).

The modified Taylor expansion yields filter coefficients h_n with its response H(k_x) that vanishes beyond a cutoff value for k_x in the wavenumber domain that depends on the scalar M. Also, the response of the extrapolator based on the modified Taylor expansion satisfies the stability constraint H(k_x) < 1 in the passband region of the filter. The cutoff wavenumber k_x determines the maximum dip accuracy of the extrapolation filter. The larger the number of filter coefficients 2N + 1, the steeper the dip accuracy. In practice, extrapolation filter lengths 7, 11 and 25 are often associated with 30-, 50- and 70-degree dip accuracies.

An alternative method for computing the filter coefficients h_n with more accuracy at large wavenumbers k_x is based on the Remez exchange algorithm . The objective is to minimize the error function

$E(k_{x})=W(k_{x})[D(k_{x})-H(k_{x})],$

(D-71)

where W(k_x) is a weighting function, such that

$||E||_{\infty }=\max |E(k_{x})|$

(D-72)

is minimum. The criterion ||E||_∞ means that the actual response H(k_x) is equiripple. Specifically, the response H(k_x) has extrema of the same absolute value with alternating signs within the passband region specified by a cutoff wavenumber k_x. The objective of the error function is to minimize the deviation of the extrema from the desired value of unity.

D.6 Optimum depth step

The objective is to specify an optimum depth step size that yields minimum-phase errors in wave extrapolation as part of the design of finite-difference migration algorithms. First, we shall review implicit and explicit schemes. Then, we shall derive equations for optimum depth size for frequency-wavenumber implicit schemes.

Start with the one-way wave equation (D-5):

${\frac {\partial }{\partial z}}P(k_{x},z,\omega )=-ik_{z}P(k_{x},z=0,\omega ),$

(D-73)

whose solution is

$P(z+\Delta z)=P(z)\ \exp \,(-ik_{z}\Delta z),$

(D-74)

where, for simplicity, k_x and ω variables are omitted from P, and

$k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}}}.$

(D-75)

The phase of the exact extrapolation operator

$\exp \left(-ik_{z}\Delta z\right)$

(D-76a)

is

$\Phi =k_{z}\Delta z,$

(D-76b)

and its amplitude is

$A=1.$

(D-76c)

Discretize the one-way wave equation (D-73):

${\frac {\delta P}{\delta z}}=-ik_{z}P.$

(D-77)

Apply differencing approximation using the explicit scheme

${\frac {P(z+\Delta z)-P(z)}{\Delta z}}=-ik_{z}P(z),$

(D-78a)

and the implicit scheme

${\frac {P(z+\Delta z)-P(z)}{\Delta z}}=-ik_{z}{\frac {P(z+\Delta z)-P(z)}{2}}.$

(D-78b)

Rewrite the explicit equation (D-78a) as

$P(z+\Delta z)=P(z)(1-ik_{z}\Delta z).$

(D-79a)

The phase of the explicit extrapolation operator

$(1-ik_{z}\Delta z)$

(D-79b)

is

${\widehat {\Phi }}=\tan ^{-1}(k_{z}\Delta z),$

(D-79c)

and its amplitude is

${\hat {A}}={\sqrt {1+(k_{z}\Delta z)^{2}}}.$

(D-79d)

Note that the amplitude of the explicit operator is greater than unity and grows from one extrapolation step to the next. In fact, large depth steps cause the operator to yield unstable results early in the extrapolation process. In general, explicit schemes tend to be unstable unless special design considerations are made (Section D.5).

Rewrite the implicit equation (D-78b) as

$P(z+\Delta z)=P(z)\left[{\frac {1-ik_{z}\Delta z/2}{1+ik_{z}\Delta z/2}}\right].$

(D-80a)

The phase of the implicit extrapolation operator

${\frac {1-ik_{z}\Delta z/2}{1+ik_{z}\Delta z/2}}$

(D-80b)

is

${\widehat {\Phi }}=\tan ^{-1}\left[{\frac {k_{z}\Delta z}{1-(k_{z}\Delta z/2)^{2}}}\right],$

(D-80c)

and its amplitude is

${\hat {A}}=1.$

(D-80d)

Note that the amplitude of the implicit scheme is unity. This means that implicit schemes are unconditionally stable.

We now specialize the phase of the frequency-space implicit scheme given by equation (D-80b) for the 65-degree dispersion relation as in equation (D-58) and obtain

${\widehat {\Phi }}=\tan ^{-1}\left[{\frac {{\hat {k}}_{z}\Delta z}{1-({\hat {k}}_{z}\Delta z/2)^{2}}}\right],$

(D-81a)

where

${\hat {k}}_{z}={\frac {\alpha {\hat {k}}_{x}^{2}}{{\frac {2\omega }{v}}-\beta {\frac {v}{2\omega }}{\hat {k}}_{x}^{2}}},$

(D-81b)

with α and β specified as in Table D-1 for the 65-degree scheme. Equation (D-81a) is the time-retarded form of equation (D-58) and hence corresponds to the second term of the latter. The discrete forms of the transform variables are :

${\hat {k}}_{x}^{2}={\frac {4}{\Delta x^{2}}}{\frac {\sin ^{2}\left({\frac {\omega }{v}}\Delta x\sin \theta \right)}{1-4b\sin ^{2}\left({\frac {\omega }{v}}\Delta x\sin \theta \right)}},$

(D-81c)

and

${\hat {\omega }}={\frac {2}{\Delta t}}\tan \,({\frac {\omega \Delta t}{2}}).$

(D-81d)

The scalar b in equation (D-81c) is set to a value between 1/12 and 1/6. The phase error for the 65-degree implicit scheme then is given by

$\Delta \Phi ={\widehat {\Phi }}-\Phi ,$

(D-82)

where ${\widehat {\Phi }}\ {\text{and}}\ \Phi$ are given by equations (D-81a) and (D-76a), respectively. Refer to the equations (D-81a,D-81b,D-81c,D-81d) and note that the phase error depends on Δx, Δt, v, θ, ω, and the depth step Δz.

Figure D-8 shows contour plots of the phase error defined by equation (D-82) for three ω/v values as a function of dip angle and depth step size. Note the complicated behavior of the contours which indicates that the optimum depth step size is associated with a complicated interdependence of the various parameters — dip, frequency, velocity, spatial, and temporal sampling rates. Therefore, in practice, migration algorithms that require wave extrapolation at discrete depth steps usually do not incorporate an automated estimation of optimum depth steps. Instead, a constant value between one-half and one dominant period, 20 to 40 ms, depending on maximum reflector dip is specified in practice. See frequency-space migration in practice for the practical aspects of frequency-space migration.

Figure D-8  Phase error contours in degrees associated with the 65-degree implicit extrapolator (see Section D.6) for a range of dip angles θ versus depth step sizes Δz, and for specific ratios of ω/v, from top to bottom, 0.001 (low frequency or high velocity), 0.06 (medium frequency and velocity) and 0.2 (high frequency or low velocity). The parameter b in equation (D-81c) was set to 0.14 for all three cases. (Computation by Dave Nichols.)

D.7 Frequency-wavenumber migration

We start with the solution of the scalar wave equation for the zero-offset wavefield as given by equation (D-7) and assume a horizontally layered earth model associated with a vertically varying velocity function v(z). By inverse Fourier transforming equation (D-7), we have

$P(x,z,t)=\int \int P(k_{x},0,\omega )\ \exp(-ik_{z}z)\ \exp(-ik_{x}x+i\omega t)\,dk_{x}\,d\omega ,$

(D-83a)

where k_z is defined by equation (D-6) adapted to the exploding reflectors model by replacing v with v/2:

$k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}}}.$

(D-83b)

The imaging principle t = 0 then is applied to get the migrated section P(x, z, t = 0),

$P(x,z,t=0)=\int \int P(k_{x},0,\omega )\ \exp(-ik_{x}x-ik_{z}z)\,dk_{x}\,d\omega .$

(D-84)

This is the equation for the phase-shift method . Equation (D-84) involves integration over frequency and inverse Fourier transformation along midpoint axis x. Refer to Figure 4.1-29 for a flow diagram of phase-shift migration.

We now consider the special case of constant velocity v. Stolt devised a migration technique that involves an efficient mapping in the 2-D Fourier transform domain from temporal frequency ω to vertical wavenumber k_z. We rewrite equation (D-83b) to get

$\omega ={\frac {v}{2}}{\sqrt {k_{x}^{2}+k_{z}^{2}}}.$

(D-85)

By keeping the horizontal wavenumber k_x unchanged and differentiating, we obtain

$d\omega ={\frac {v}{2}}{\frac {k_{z}}{\sqrt {k_{x}^{2}+k_{z}^{2}}}}dk_{z}.$

(D-86)

When equations (D-85) and (D-86) are substituted into equation (D-84), we get

$P(x,z,t=0)=\int \int \left[{\frac {v}{2}}{\frac {k_{z}}{\sqrt {k_{x}^{2}+k_{z}^{2}}}}\right]\ P\left[k_{x},0,{\frac {v}{2}}{\sqrt {k_{x}^{2}+k_{z}^{2}}}\right]\ \exp(-ik_{x}x-ik_{z}z)\,dk_{x}dk_{z}.$

(D-87)

This is the equation for constant-velocity Stolt migration. It involves two operations in the f − k domain. First, the temporal frequency ω is mapped onto the vertical wavenumber k_z via equation (D-85). This is the same as mapping point B′ onto point B in Figure 4.1-25. Second, the amplitudes are scaled by the quantity

$S={\frac {v}{2}}{\frac {k_{z}}{\sqrt {k_{x}^{2}+k_{z}^{2}}}},$

which is equivalent to the obliquity factor associated with Kirchhoff migration (migration principles). Refer to Figure 4.1-30 for a flow diagram of constant-velocity Stolt migration.

To extend the algorithm to the variable-velocity case, yet retain efficiency, Stolt did a coordinate transformation that involves stretching the time axis to make the scalar wave equation velocity independent. A summary of the theoretical procedure is given here. Consider wavefield P(x, z, t) and the transformed wavefield P(x, d, T):

$P(x,z,t)=P(x,d,T),$

where T is the stretched time axis and d is the output variable (equivalent of z) for migration in the stretched coordinate system. Here, the x variable is identical in both coordinate systems. This coordinate transformation basically is equivalent to stretching the data using the rms velocities

$T(t)={\frac {1}{c}}\left[2\int _{0}^{t}dt'v_{rms}^{2}(t')t'\right]^{1/2},$

(D-88)

where

$v_{rms}^{2}(t)={\frac {1}{t}}\int _{0}^{t}v^{2}(t')dt',$

(D-89)

and c is an arbitrary reference velocity used to maintain the vertical axis as time after the coordinate transformation from t to T. After some tedious algebra, the time-retarded scalar wave equation takes the following form in the stretched coordinates :

${\frac {\partial ^{2}P}{\partial x^{2}}}+W{\frac {\partial ^{2}P}{\partial d^{2}}}+{\frac {2}{c}}{\frac {\partial ^{2}P}{\partial d\partial T}},$

(D-90)

where W is a complicated function of velocity and coordinate variables. In practice, it normally is set to a constant between 0 and 1. The procedure for Stolt migration with stretch follows:

• Start with the stacked section, which is assumed to be a zero-offset section P(x, z = 0, t).
• Convert this time section to stretched section P(x, d = 0, T) by the coordinate transformation (equation D-88).
• 2-D Fourier transform the stretched section P(k_x, d = 0, ω_T).
• Apply the following mapping function to perform migration:

$k_{d}=\left(1-{\frac {1}{W}}\right){\frac {\omega _{T}}{c}}-{\frac {1}{W}}{\sqrt {{\frac {\omega _{T}^{2}}{c^{2}}}-Wk_{x}^{2}}}.$

(D-91)

This equation is based on the dispersion relation of the retarded wave equation in the stretched coordinates (equation D-90). The expression for the output from migration is (omitting the heavy algebra)

$P(k_{x},k_{d},0)=\left[{\frac {c}{2-W}}(1-W+{\frac {1}{K}})\right]P\left\{k_{x},0,{\bigg [}{\frac {ck_{d}}{2-W}}(1-W+K)\right\},$

(D-92a)

where

$K={\frac {1}{\sqrt {1+(2-W){\frac {k_{x}^{2}}{k_{d}^{2}}}}}}.$

(D-92b)

• 2-D inverse Fourier transform the migrated section in the stretched coordinates P(x, d, T = 0).
• Convert back to the familiar space-time coordinates P(x, z, t = 0). This is the final migrated section.

To derive the equations for Stolt migration with the output in time τ = 2z/v, use the dispersion relation as in equation (D-49a) in lieu of equation (D-83b) and replace the vertical wavenumber k_z with the output frequency ω_τ = vk_z/2 in equations (D-84) through (D-88).

When W = 1, equation (D-91) takes the simple form

$k_{d}={\sqrt {{\frac {\omega _{T}^{2}}{c^{2}}}-k_{x}^{2}}},$

(D-93)

which makes mapping of equation (D-91) equivalent to the constant-velocity Stolt algorithm.

Note that Stolt migration with stretch tries to handle velocity variations. However, it is no substitute for depth migration; it only accommodates velocity variations that can be handled by time migration. Figure D-9 shows the flow diagram for Stolt migration with stretch. See frequency-wavenumber migration in practice for the practical aspects of frequency-wavenumber migration.

Figure D-9  Flowchart for Stolt migration with stretch.

D.8 Residual migration

The relationship between the output and input temporal frequencies for a 90-degree time migration operator (rewritten from equation D-49a) is given by

$\omega _{\tau }={\sqrt {\omega ^{2}-{\frac {v^{2}k_{x}^{2}}{4}}}},$

(D-94)

where ω_τ = vk_z/2 is the output frequency, ω is the input frequency, v is the true medium velocity, and k_x is the midpoint wavenumber. Equation (D-94) refers to single-pass migration if the migration velocity is the same as the medium velocity.

Consider a two-pass migration, first with velocity v₁, followed by a second migration with velocity v₂. The output frequency from the first pass is

$\omega _{1}={\sqrt {\omega ^{2}-{\frac {v_{1}^{2}k_{x}^{2}}{4}}}},$

(D-95a)

which is the input frequency for the second pass

$\omega _{2}={\sqrt {\omega _{1}^{2}-{\frac {v_{2}^{2}k_{x}^{2}}{4}}}},$

(D-95b)

The horizontal wavenumber k_x is fixed, since it is invariant under migration.

If the output frequency from the single-pass migration given by equation (D-94) is set equal to that from the two-pass migration given by equation (D-95b), then the necessary relationship between the residual migration velocity v₂ and the medium velocity v can be established as

$v^{2}=v_{1}^{2}+v_{2}^{2}.$

(D-96a)

One practical scheme for residual migration is implemented as follows:

1. The constant-velocity Stolt migration with a stretch factor of W = 1 is used for the first pass. Velocity v₁ in equation (D-96a) is chosen as the minimum value in the actual velocity field.
2. For the second pass, a dip-limited finite-difference migration can be used. The first-pass migration brings the dips down to within the range that the second-pass finite-difference migration can accommodate, accurately. Remember that the velocity field for the second pass is computed by the the relation

$v_{2}={\sqrt {v^{2}-v_{1}^{2}}},$

(D-96b)

where v, v₁, and v₂ are the original, first- and second-pass migration velocities, respectively. A multiple-pass application of residual migration is referred to as cascaded migration. Refer to Sections finite-difference migration in practice and frequency-wavenumber migration in practice for the practical aspects of cascaded migration and residual migration, respectively.

Other figures

• 

  Figure 4.1-18  (a) Superposition of the time sections in Figure 4.1-17; (b) removing the translational effect by retardation to place the energy at the apex of the hyperbola obtained initially along the beach line.

• 

  Figure 4.1-25  Migration in the f − k domain. (Migration in the t − x domain is illustrated in Figure 4.1-1.) (a) A dipping reflector is represented by a radial line OB in the f − k plane. (b) After migration, the radial line OB maps onto another radial line OB′, while B maps onto B′. The horizontal wavenumber is invariant under migration. For comparison, the f − k response of the dipping event before migration (a) has been superimposed on the f − k response after migration (b). Adapted from .

• 

  Figure 4.4-1  Impulse responses of the frequency-space implicit schemes with various degrees of approximations to the one-way scalar wave equation. (See Section D.4 for the theoretical basis of these responses.)

• 

  Figure 4.1-29  Flowchart for Gazdag’s phase-shift method of migration.

• 

  Figure 4.1-25  Migration in the f − k domain. (Migration in the t − x domain is illustrated in Figure 4.1-1.) (a) A dipping reflector is represented by a radial line OB in the f − k plane. (b) After migration, the radial line OB maps onto another radial line OB′, while B maps onto B′. The horizontal wavenumber is invariant under migration. For comparison, the f − k response of the dipping event before migration (a) has been superimposed on the f − k response after migration (b). Adapted from .]].

• 

  Figure 4.1-30  Flowchart for Stolt’s constant-velocity migration method in the f − k domain.

See also

• Introduction to migration
• Migration principles
• Kirchhoff migration in practice
• Finite-difference migration in practice
• Frequency-space migration in practice
• Frequency-wavenumber migration in practice
• Further aspects of migration in practice
• Exercises

External links

find literature about Mathematical foundation of migration