Mathematical foundation of 3-D migration

From SEG Wiki
Jump to: navigation, search
Other languages:
English • ‎español
Seismic Data Analysis
Seismic-data-analysis.jpg
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


G.1 Implicit methods

The 2-D scalar wave equation (D-1) given in Section D.1 can be adapted to three dimensions as


(G-1

)

This equation describes propagation of a 3-D compressional zero-offset wavefield P(x, y, z, t) in a medium with constant material density and compressional wave velocity v(x, y, z), where x is the horizontal spatial axis in the inline direction, y is the horizontal spatial axis in the crossline direction, z is the depth axis (positive downward), and t is time. Given the upcoming seismic wavefield P(x, y, z = 0, t), which is recorded at the surface, we want to determine reflectivity P(x, y, z, t = 0). This requires extrapolating the surface wavefield to depth z, then collecting it at t = 0.

The solution to the 2-D scalar wave equation (D-1) is given by equation (D-7), which can be adapted to extrapolate a 3-D zero-offset wavefield in depth


(G-2

)

where


(G-3

)

with


(G-4a

)

and


(G-4b

)

In equations (G-4a,G-4b), kx and ky are the wavenumbers in the inline and crossline directions, respectively, and ω is the temporal frequency in units of radians per unit time. Note that the inline and crossline normalized wavenumbers X and Y are coupled in equation (G-3).

First, consider the constant-velocity case and rewrite the dispersion relation of equation (G-3) as for the 2-D case (equation D-49a) to get


(G-5)

where ωτ = vkz/2 is the output temporal frequency.

Suppose we perform 2-D migration on all the inlines by mapping the input frequency ω to an output frequency ω1 using the relation


(G-6)

We then sort the already migrated inlines in the crossline direction. Next, suppose that we perform 2-D migration on all the crosslines. This second-pass migration is performed by mapping ω1 of equation (G-5) to an output frequency ω2 using the relation


(G-7)

When equation (G-7) is substituted into equation (G-6), we obtain


(G-8)

which is the same as equation (G-5). This means that 3-D migration can be performed in two passes: 2-D migration in the inline direction followed by 2-D migration in the crossline direction. The above derivation of a 90-degree 3-D migration operator Levin [1] is based on the constant-velocity assumption.

Mathematically, the two-pass approach is based on the idea of full separation [2] of migration operators in the inline and crossline directions. The separation is strictly valid for the constant-velocity medium. When velocities vary spatially, splitting [2] of the inline and crossline operators is necessary (Figure 7.3-1). A rigorous mathematical treatment of separation and splitting of extrapolation operators can be found in a paper by Brown /> [3].

To implement separation or splitting, the inline and crossline wavenumbers need to be decoupled in the 3-D dispersion relation (equation G-3). By Taylor series expansion of the square root and by keeping the first three terms, we get the 15-degree dispersion relation (Section D.3)


(G-9

)

When the cross-dip component is zero (Y = 0), equation (G-9) reduces to its 2-D form (equation D-54). The inline and crossline terms in equation (G-9) are separable when velocity is slowly variable in z or independent of x and y [3].

The impulse response of the 2-D 15-degree migration operator with its dispersion relation given by equation (D-42a) is an ellipse (Figure 4.3-1). The impulse response of the 3-D 15-degree migration operator with its dispersion relation given by equation (G-9) is an ellipsoid (Figure 7.3-2). Except for the amplitudes, the vertical cross-section at center line (51) is the same as the 2-D impulse response. Within typical bandwidth of seismic data, the 3-D 15-degree impulse response departs from the ideal 3-D impulse response (a hollow hemisphere) at dips greater than approximately 30 degrees. The time slices of the impulse response for equation (G-9) are circular in shape.

Additional problems arise when more accuracy is desired; the X and Y terms in equation (G-3) no longer are decoupled. For example, the 45-degree approximation (Section D.4) to equation (G-3) introduces cross terms. Another approximation is to represent the single square root in equation (G-3) by two square roots [4] as


(G-10

)

Note that equation (G-10) has the same form as equation (D-15). By Taylor series expansion of the two square roots, equation (G-10) reduces to equation (G-9). Assuming a small cross-dip component, more accuracy in the inline direction is achieved by making the 45-degree approximation given by


(G-11)

When more accuracy is required in both directions, we can apply the 45-degree expansion to equation (G-10) in both the inline and crossline directions given by


(G-12)
Figure G-1  An algorithmic description of one-pass implicit 3-D poststack migration based on 3-D wave extrapolation and imaging in the frequency-space domain (Section G.1).

The impulse response of the 2-D 45-degree migration operator with its dispersion relation given by equation (D-58) is heart shaped (Figure 4.3-1). The impulse response of the 3-D 45-degree migration operator with its dispersion relation given by equation (G-12) is shown in Figure 7.3-3. Except for the amplitudes, the vertical cross-section at center line (51) is the same as the 2-D impulse response. Within the typical bandwidth of seismic data, the 3-D 45-degree impulse response departs from the ideal 3-D impulse response (a hollow hemisphere) at dips greater than approximately 65 degrees. The time slices of the impulse response for equation (G-12) are diamond shaped. Note that the response is accurate in both the inline and crossline directions, and that the largest error in the form of undermigration occurs in the two diagonal directions.

In summary, to implement separation or splitting, the inline and crossline terms of the 3-D migration operator need to be decoupled. This can be achieved by making rational approximations (equation G-9 or G-12) to the square root in equation (G-3). Separation works for 90-degree, constant-velocity (i.e., Stolt migration) or 15-degree, slowly varying v(z) cases. When there are large vertical velocity gradients or strong lateral velocity variations, splitting must be used. The algorithm for a one-pass implicit frequency-space 3-D poststack migration that uses the splitting method is identical to the 2-D algorithm outlined in Figure 4.1-23, except for the additional application of the diffraction term in the crossline direction (Figure G-1).


G.2 Explicit methods

Rewrite equation (G-2) for extrapolating a 3-D zero-offset wavefield in discrete depth steps of Δz


(G-13a

)

Combine equations (G-3) and (G-4) to express the vertical wavenumber kz in terms of the other transform variables as


(G-13b

)

where, v is the extrapolation velocity.

Substitute


(G-14

)

into equation (G-13b), and obtain


(G-15

)

With the vertical wavenumber kz given by equation (G-15), we now write the desired extrapolation operator D(k) given by the exponential term in equation (G-13a) for a specific ratio of 2ω/v as


(G-16

)

This operator can be, in principle, inverse transformed back to the frequency-space domain, and thus, extrapolation described by equation (G-13) can be achieved by convolution of the wavefield with the extrapolation filter in the frequency-space domain.

Here, we describe a practical scheme [5] > [6] which involves the design of a stable 2-D extrapolation filter, use of the McClellan transform to obtain a 3-D extrapolation operator without actually computing a 3-D operator, and replacing convolution for filter application by an efficient recursive scheme based on Chebychev polynomials.

Start with the extrapolation operator in the Fourier transform domain given by equation (G-16). To accommodate lateral velocity variations, we want to compute an extrapolation filter hn in the frequency-space (ω, x, y) domain with a Fourier transform that best approximates the desired transform given by equation (G-16). Since D(k) is symmetric with respect to k = 0, the complex filter coefficients hn are even — h−n = hn, and the number of filter coefficients 2N + 1 is odd. As a result, we only need to compute N coefficients. The actual Fourier transform of hn is expressed as


(G-17

)

Methods to obtain hn such that the actual Fourier transform H(k) given by equation (G-17) closely approximates the actual Fourier transform D(k) given by equation (G-16) is described by Holberg [7] and Hale [5]. The former is based on least-squares and the latter is based on a modified Taylor series. The Hale method is described in Section D.5. Since the desired amplitude spectrum is |D(k)| = 1 for all k, then we require the actual amplitude spectrum to be |H(k)| ≤ 1 for all k. This will ensure that the operator is stable — that is, amplitudes of the extrapolated wavefield will not grow from one depth level to another.

By using a special case of the recursive formula for the Chebychev polynomials [6], equation (G-17) can be rewritten in terms of cos k as


(G-18

)

so that H(k) is expressed entirely in terms of hn and cos k. As a result, convolution with the filter hn can be performed by a recursive Chebychev filter structure described by Hale [6]. This is essential for an efficient application of a 2-D explicit extrapolation filter, since direct convolution in two dimensions is very costly.

Table G-1. The McClellan filter coefficients associated with the Fourier transform given by equation (G-20).
1/8 1/4 1/8
1/4 −1/2 1/4
1/8 1/4 1/8

For the design of the 3-D explicit operator, define the transform G(k) = cos k and note from equation (G-14) that


(G-19

)

For computational efficiency, G(kx, ky) given by equation (G-19) may be approximated as [6]


(G-20

)

This approximation is exact for kx = 0 and ky = 0, best for small kx and ky, and deteriorates at large kx and ky.

G(kx, ky) of equation (G-20) is the Fourier transform of the 2-D McClellan transformation filter given by Table G-1 in the x − y domain, independent of ω and v [6].

Using the McClellan filter coefficients given by Table G-1, the 1-D extrapolation filter hn is transformed into a 2-D extrapolation filter with approximate circular symmetry [6]. Again, by way of a Chebychev filter structure, application of this filter is done more efficiently as compared with direct convolution. The impulse response of the 3-D explicit migration operator that uses the 3 × 3 McClellan transform template of Table G-1 is shown in Figure 7.3-13.

The circular symmetry of the McClellan transform filter can be improved by an alternative approximation to G(kx, ky) of equation (G-19) by [6]


(G-21

)

where c is an adjustable scalar to make the Fourier transform of the 2-D filter based on the McClellan transform match the desired Fourier transform exactly (equation G-16) at some A typical value for c = 0.0255. G(kx, ky) of equation (G-21) is the Fourier transform of the 2-D McClellan transformation filter given by Table G-2 in the x − y domain, again, independent of ω and v.

Table G-2. The McClellan filter coefficients associated with the Fourier transform given by equation (G-21).
c/8 0 c/4 0 c/8
0 1/8 1/4 1/8 0
c/4 1/4 −(1 + c)/2 1/4 c/4
0 1/8 1/4 1/8 0
c/8 0 c/4 0 c/8

The impulse response of the 3-D explicit migration operator that uses the 5 × 5 McClellan transform template of Table G-2 is shown in Figure 7.3-14.

Aside from the application of the McClellan transform, the circular symmetry of the 3-D impulse response can also be attained by applying the correction term by Li [8] to the extrapolated wavefield, intermittently, at some depth steps. The objective is to compensate for the accumulated errors caused by the approximations made to the exact 3-D extrapolation operator given by equation (G-16). Li [8] devised a correction term originally to compensate for the errors caused by the 45-degree split form of the dispersion relation given by equation (G-12). The error incurred during wave extrapolation using an implicit scheme arises from the difference between the exact form of the dispersion relation given by equation (G-13b) and an approximate form such as that given by equation (G-12).

Etgen and Nichols [9] later adapted Li’s correction to compensate for the errors caused by an approximation made to the explicit operator given by equation (G-16). The error incurred during wave extrapolation using an explicit scheme can be attributed to the difference between the exact form of the wavenumber k given by equation (G-14) and an approximate form such as that given by equation (G-19)


(G-22

)

with G(kx, ky) given by the approximate form in equation (G-21).

Write equation (G-16) for the approximate form of the wavenumber k given by equation (G-22) as


(G-23)

By using the exact (equation G-16) and the approximate (equation G-23) forms of the extrapolation operators, the error in the extrapolation accumulated after nΔz depth intervals is then given by


(G-24)

Therefore, the extrapolated wavefield P(kx, ky, nΔz, ω) given by equation (G-13a) is compensated for the errors incurred by the approximate form of G(kx, ky) given by equation (G-21) by applying the conjugate of E(kx, ky, ω) of equation (G-24) to the extrapolated wavefield


(G-25

)

Following Li’s correction adapted to the explicit schemes [9] as given by equation (G-25), the wave extrapolation is continued with the approximate form (equation G-23) until more error is accumulated.


G.3 3-D phase-shift migration

We start with the solution of the 3-D scalar wave equation to extrapolate a 3-D zero-offset wavefield P(x, y, z = 0, t) as given by equation (G-2) in the frequency-wavenumber domain. We also assume a horizontally layered earth model associated with a vertically varying velocity function v(z). By inverse Fourier transforming equation (G-2), we have


(G-26

)

where kz is defined by equation (G-13) adapted to the exploding reflectors model by replacing v with v/2.

The imaging principle t = 0 then is applied to equation (G-26) to get the migrated section P(x, y, z, t = 0) as


(G-27

)

This is the equation for the phase-shift method [10] adapted to 3-D migration. Equation (G-27) involves integration over frequency and inverse Fourier transformation along inline and crossline axes at each depth step.


G.4 3-D stolt migration

We now consider the special case of constant velocity v. As described in Section D.7, Stolt [11] devised a migration technique that involves an efficient mapping in the Fourier transform domain from temporal frequency ω to vertical wavenumber kz. We rewrite equation (G-13) to get an explicit expression for ω as


(G-28

)

By keeping the horizontal wavenumbers kx and ky unchanged and differentiating, we get


(G-29)

When equations (G-28) and (G-29) are substituted into equation (G-27), we obtain


(G-30)

This is the equation for constant-velocity 3-D Stolt migration. It involves two operations in the f − k domain. First, the temporal frequency ω is mapped onto the vertical wavenumber kz via equation (G-28). Second, the amplitudes are scaled by the quantity S:


(G-31)

which is equivalent to the obliquity factor associated with Kirchhoff migration (Section H.1).


G.5 Trace interpolation

In Section F.6 we developed the theory of the spatial prediction filter. The objective at the time was to attenuate random noise by applying the spatial prediction filter to data in the frequency-space domain [12].

Now, we want to design a complex Wiener prediction filter to interpolate data in the spatial direction. Consider a CMP-stacked data set P(x, t), where x is the CMP axis and t is the two-way zero-offset time axis. Apply Fourier transform in the t direction to decompose this 2-D data set to its frequency components P(x, ω). For each frequency component, define a complex array P : P(x, ω) in the x direction. Specifically, we want a filter F : F(x) such that, when applied to the input data array P : P(x, ω), it yields an estimate of the input array D : P(x + Δx/2, ω), at x + Δx/2, where D is the desired output array and Δx/2 is the prediction lag [13].

The spatial prediction filtering is expressed by the following convolutional relation


(G-32

)

where Y(x) represents the actual output from prediction filtering. Consider the discrete form of equation (G-32), with F(x) represented by the m-length complex array F : (F0, F1, F2, …, Fm−1), P(x, ω) represented by the n-length complex array P : (P0, P1, P2, …, Pn−1), and Y(x) represented by the (m + n − 1)-length complex array Y : (Y0, Y1, Y2, …, Ym+n−1). Equation (G-32) can then be expressed in matrix form


(G-33

)

To develop the theory for trace interpolation using spatial prediction filters, consider the special case of a three-point one-step prediction filter (F0, F1, F2) and a six-point input data array (P0, P1, P2, P3, P4, P5). We also want to specialize the output array Y(x) to be the input array P(x) one unit of distance ahead, (P1, P2, P3, P4, P5). Equation (G-33) for this special case takes the form


(G-34

)

Augment the column array on the right-hand side to the matrix on the left-hand side to obtain


(G-35

)

Write out explicitly the equations associated with rows 3, 4, and 5 which contain terms with all three filter coefficients (F0, F1, F2) as


(G-36a

)


(G-36b

)


(G-36c

)

Next, split each of these three equations into two parts — first part has terms that contain the data array elements (P0, P2, P4) and second part has terms that contain the data array elements (P1, P3, P5). Then, write the two parts in matrix form


(G-37

)

This equation suggests that given the prediction filter (F0, F1, F2), the data array (P1, P3, P5) can be calculated from the data array (P0, P2, P4). The known data array (P0, P2, P4) is associated with three consecutive traces in a CMP-stacked data set. The unknown array to be calculated (P1, P3, P5) would then be associated with the traces halfway between the traces that correspond to the known data array (P0, P2, P4). Hence, we note that equation (G-37) can be used to perform spatial interpolation.

Write equation (G-37) in matrix notation


(G-38

)

where Fd and Fu are the coefficient matrices with their elements in terms of the prediction filter coefficients (F0, F1, F2), and Pd : (P0, P2, P4) and Pu : (P1, P3, P5) are the known input data array and the interpolated data array, respectively.

Equation (G-38) is of the form d′ = Lu as for a generalized linear inverse problem (Section F.3). The least-squares solution u = (LT*L)−1LT*d is adopted for the interpolated data array Pu. In the present case, the matrix L corresponds to Fu, and the array d corresponds to FdPd, so that


(G-39

)

where T denotes matrix transpose and the asterisk denotes complex conjugate.

The question remains as to how to calculate the prediction filter from the known data (P0, P2, P4). We shall make the assumption that the input data set is made up of M traces that comprise a set of N dipping events represented by the following combination of amplitude and phase spectra in the frequency-space domain [13]:


(G-40

)

where Aj(ω) is the amplitude spectrum of the jth dipping event, Δτj is the time shift along the jth dipping event from trace to trace, and k is the trace index. Define


(G-41

)

so that equation (G-40) takes the form


(G-42

)

where the variable ω is omitted for convenience.

Now, given the data array (P0, P2, P4) and the one-step prediction filter (F0, F1, F2), we want to compute the next element of the data array P6. Refer to equation (G-36a) and replace the elements of the data array (P0, P1, P2) with the elements of the data array (P0, P2, P4). Also, replace the next element P3 of the data array (P0, P1, P2) with the next element P6 of the data array (P0, P2, P4) to get


(G-43

)

Write this equation using the form given by equation (G-42)


(G-44

)

Consider a case of three dipping events. Set N = 3 in equation (G-44) to get


(G-45a

)


(G-45b

)


(G-45c

)

which can be written in matrix form, while explicitly incorporating the variable ω [14],


(G-46

)

Equation (G-46) can be solved for the prediction filter (F0, F1, F2) for each frequency component of the input data array (P0, P2, P4). But this filter, when applied to the data array, yields the next element P6 (equation G-43). We need a filter that yields (P1, P3, P5, …), — the elements halfway between the elements of the data array (P0, P2, P4, …) equivalent to the output from interpolation.

To expedite a determination of this interpolation filter, derive the prediction filter associated with the data array that would have been created by interpolation, (P0, P1, P2). Write equation (G-43) for the case of the data array (P0, P1, P2) to compute the next element P3 using the one-step prediction filter as


(G-47

)

Write this equation using the form given by equation (G-42)


(G-48)

For the case of three dipping events, setting N = 3 in equation (G-48), we get


(G-49a

)


(G-49b

)


(G-49c

)

which can be written in matrix form, while explicitly incorporating the variable ω,


(G-50

)

Now, note from the definition expressed by equation (G-41) that


(G-51

)

With this definition, rewrite the matrix equation (G-50) as


(G-52

)

Additionally, rewrite equation (G-46) explicitly in terms of ω/2 as


(G-53

)

Finally, compare equations (G-52) and (G-53) and note that the prediction filter that we need to do trace interpolation of the data component with frequency ω actually is the prediction filter [F0(ω/2), F1(ω/2), F2(ω/2)] that is computed directly from the uninterpolated data component with half the frequency ω/2.

Figure G-2  Rayapth geometry for 3-D nonzero-offset traveltime equation (G-59) derived in Section G.6.


G.6 3-D nonzero-offset traveltime equation

Refer to the 3-D raypath geometry shown Figure G-2. For convenience, we shall consider a 3-D recording parallel to the inline direction x with zero source-receiver azimuth. Specifically, consider a midpoint M : (x, y, 0) on crossline y that is associated with a source S : (x − h, y, 0) and a receiver R : (x + h, y, 0), where h is the half-offset. The distance along the raypath SD from source S : (x − h, y, 0) to a point D(0, 0, z) in the subsurface is given by


(G-54a

)

and the distance along the raypath DR from point D : (0, 0, z) to a receiver R : (x + h, y, 0) is given by


(G-54b

)

The traveltime t associated with the raypath SDR then is given by


(G-55

)

where v is the medium velocity.

We shall show that equation (G-55) represents an ellipsoid (Figure G-2b). Square both sides to get


(G-56

)

Combine the second and third terms on the right-hand side and simplify the terms inside the square root


(G-57

)

Perform further algebraic manipulation to collect the terms in x, y, and z


(G-58

)

Finally, normalize by the terms on the right-hand side and rearrange the terms in the denominators to obtain


(G-59

)

Equation (G-59) represents an ellipsoid. The horizontal cross-section of this ellipsoid at a depth z is an ellipse. At z = 0, the ellipse has the following parameters:

  1. Semi-major axis in inline x direction: a = vt/2.
  2. Semi-minor axis in crossline y direction: .
  3. Distance from center to either focus: .

The ellipsoid of equation (G-59) in the x − y − z volume describes the kinematics of the impulse response of a 3-D nonzero-offset migration operator applied to 3-D prestack data.

When equation (G-59) is specialized to the zero-offset case, h = 0, we get


(G-60a

)

which describes a hollow hemisphere in the (x, y, z) volume for a constant t with a radius vt/2. This hemisphere represents the kinematics of the impulse response of a 3-D zero-offset migration operator applied to 3-D poststack data.

When equation (G-55) is specialized to the 3-D zero-offset case, h = 0, we get


(G-60b

)

which describes the diffraction hyperboloid of revolution in the (x, y, t) volume for a constant z.[15] [16] [17]


Additional images


References

  1. Jacubowicz and Levin (1983), Jakubowicz, H. and Levin, S., 1983, A simple exact method of 3-D migration — theory: Geophys. Prosp., 31, 34–56.
  2. 2.0 2.1 Claerbout, 1985, Claerbout, J. F., 1985, Imaging the earth’s interior: Blackwell Scientific Publications.
  3. 3.0 3.1 Brown, 1983, Brown, D. L., 1983, Applications of operator separation in reflection seismology: Geophysics, 48, 288–294.
  4. Ristow (1980), Ristow, D., 1980, 3-D downward extrapolation of seismic data in particular by finite-difference methods: PhD thesis, University of Utrecht, The Netherlands.
  5. 5.0 5.1 Hale, 1991a, Hale, D., 1991a, Stable explicit depth extrapolation of seismic wavefields: Geophysics, 56, 1770-1777.
  6. 6.0 6.1 6.2 6.3 6.4 6.5 6.6 Hale, 1991b, Hale, D., 1991b, 3-D depth migration via McClellan transforms: Geophysics, 56, 1778–1785.
  7. Holberg, 1988, Holberg, O., 1988, Towards optimum one-way wave propagation: Geophys. Prosp., 36, 99–114.
  8. 8.0 8.1 Li (1991), Li., Z, 1991, Compensating finite-difference errors in 3-D migration and modeling: Geophysics, 56, 1650–1660.
  9. 9.0 9.1 Etgen and Nichols (1999), Etgen, J. T. and Nichols, D., 1999, Application of the Li correction to explicit depth extrapolation methods: 69th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1366–1369.
  10. Gazdag, 1978, Gazdag, J., 1978, Wave-equation migration by phase shift: Geophysics, 43, 1342–1351.
  11. Stolt (1978), Stolt, R. H., 1978, Migration by Fourier transform: Geophysics, 43, 23–48.
  12. Canales, 1984, Canales, L., 1984, Random noise reduction: 54th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 525–527.
  13. 13.0 13.1 Spitz, 1991, Spitz, S., 1991, Seismic trace interpolation in the f − x domain: Geophysics, 56, 785–794.
  14. Porsani, 1997, Porsani, M. J., 1997, Seismic trace interpolation using half-step prediction filters: Preprint.
  15. Berkhout, A. J., 1985, 3-D seismic processing with an eye to the future: World Oil, 201, No. 5, 91–94.
  16. Biondi, B., Fomel, S., and Chemingui, N., 1998, Azimuth moveout for 3-D prestack imaging: Geophysics, 63, 574–588.
  17. Chemingui, N. and Biondi, B., 1997, Equalization of irregular data by iterative inversion: 67th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1115–1118.


See also


External links

find literature about
Mathematical foundation of 3-D migration
SEG button search.png Datapages button.png GeoScienceWorld button.png OnePetro button.png Schlumberger button.png Google button.png AGI button.png