Method for controlled-amplitude prestack time migration of seismic data

[0002] This application claims the benefit of U.S. Provisional Application No. 60/107,120 filed Nov. 5, 1998 and U.S. Provisional Application No. 60/107,633 filed Nov. 9, 1998.

[0003] This invention relates generally to the field of seismic prospecting and, more particularly, to imaging of seismic data. Specifically, the invention is a method for controlled-amplitude prestack time migration of seismic data.

[0004] In the oil and gas industry, seismic prospecting techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. A seismic prospecting operation consists of three separate stages: data acquisition, data processing, and data interpretation. The success of a seismic prospecting operation is dependent on satisfactory completion of all three stages.

[0005] In the data acquisition stage, a seismic source is used to generate a physical impulse known as a “seismic signal” that propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflected signals (known as “seismic reflections”) are detected and recorded by an array of seismic receivers located at or near the surface of the earth, in an overlying body of water, or at known depths in boreholes. The seismic energy recorded by each seismic receiver is known as a “seismic data trace.”

[0006] During the data processing stage, the raw seismic data traces recorded in the data acquisition stage are refined and enhanced using a variety of procedures that depend on the nature of the geologic structure being investigated and on the characteristics of the raw data traces themselves. In general, the purpose of the data processing stage is to produce an image of the subsurface geologic structure from the recorded seismic data for use during the data interpretation stage. The image is developed using theoretical and empirical models of the manner in which the seismic signals are transmitted into the earth, attenuated by the subsurface strata, and reflected from the geologic structures. The quality of the final product of the data processing stage is heavily dependent on the accuracy of the procedures used to process the data.

[0007] The purpose of the data interpretation stage is to determine information about the subsurface geology of the earth from the processed seismic data. For example, data interpretation may be used to determine the general geologic structure of a subsurface region, or to locate potential hydrocarbon reservoirs, or to guide the development of an already discovered reservoir. Obviously, the data interpretation stage cannot be successful unless the processed seismic data provide an accurate representation of the subsurface geology.

[0008] Typically, some form of seismic migration (also known as “imaging”) must be performed during the data processing stage in order to accurately position the subsurface seismic reflectors. The need for seismic migration arises because variable seismic velocities and dipping reflectors cause seismic reflections in unmigrated seismic images to appear at incorrect locations. Seismic migration is an inversion operation in which the seismic reflections are moved or “migrated” to their true subsurface positions.

[0009] There are many different seismic migration techniques. Some of these migration techniques are applied after common-midpoint (CMP) stacking of the data traces. (As is well known, CMP stacking is a data processing procedure in which a plurality of seismic data traces having the same source-receiver midpoint but different offsets are summed to form a stacked data trace that simulates a zero-offset data trace for the midpoint in question.) Such “poststack” migration can be done, for example, by integration along diffraction curves (known as “Kirchhoff” migration), by numerical finite difference or phase-shift downward-continuation of the wavefield, or by equivalent operations in frequency-wavenumber or other domains.

[0010] Conversely, other seismic migration techniques are applied before stacking of the seismic data traces. In other words, these “prestack” migration techniques are applied to the individual nonzero-offset data traces and the migrated results are then stacked to form the final image. Prestack migration typically produces better images than poststack migration. However, prestack migration is generally much more expensive than poststack migration. Accordingly, the use of prestack migration has typically been limited to situations where poststack migration does not provide an acceptable result, e.g., where the reflectors are steeply dipping.

[0011] In some cases, reflector dip can exceed 90 degrees. As is well known in the art, it may be possible to image these “overturned” reflectors using data from seismic “turning rays.” Prestack migration techniques must be used in order to obtain an accurate image of overturned reflectors from seismic turning ray data.

[0012] There are two general types of prestack migration, prestack time migration and prestack depth migration. The present invention relates to prestack time migration which is used in situations where the subsurface seismic velocity varies in the vertical direction, but can be regarded as approximately constant laterally. Prestack time migration is widely used in the petroleum industry; it is generally considered applicable to most (but not all) prospects.

[0013] In practice, prestack time migration is usually approximated by a composite procedure involving successive steps of normal moveout correction (denoted as “NMO”), followed by dip moveout correction (denoted as “DMO”), followed by zero-offset migration (referred to as “ZOM”). This sequence of processes is performed on collections of seismic data traces having the same or nearly the same source-receiver separation distance and the same or nearly the same azimuthal orientation (referred to as “common-offset, common-azimuth gathers”). It is known from experience that the combination of these imaging steps results in a good approximation to prestack time migration in most situations. The approximation tends to be less accurate in situations where there are dipping reflectors in a medium whose velocity varies significantly with depth. In such situations, the inaccuracies manifest themselves in mispositioned images of seismic reflectors, as well as in distorted amplitudes. A principal advantage of the decomposition of prestack time migration into the NMO/DMO/ZOM sequence is that the required algorithms execute rapidly on commonly available digital computers, which results in economical processing.

[0014] Kirchhoff prestack migration is the only commonly known method of performing prestack time migration which overcomes the inaccuracies of the NMO/DMO/ZOM decomposition mentioned in the previous paragraph, and which is applicable to common offset gathers. However, Kirchhoff prestack migration, as compared to the NMO/DMO/ZOM decomposition, is generally very expensive and time consuming. For this reason, the seismic data processing industry is continuing its efforts to develop techniques that can be used to produce accurate prestack migrations in an economical manner.

[0015] In Kirchhoff prestack migration, the data in a common-offset gather are summed over all source-receiver midpoints (denoted by {right arrow over (ξ)}). During Kirchhoff summation, the data are time shifted by an amount that depends on the assumed seismic velocity structure in the earth. The time shift is denoted by the symbol t_D({right arrow over (ξ)}, {right arrow over (x)}), which depends on the midpoint location {right arrow over (ξ)} and the imaged point {right arrow over (x)}=(x,y,z). During summation, the data may be multiplied with a weight w that generally depends on {right arrow over (ξ)} and {right arrow over (x)}. As explained in Schleicher, J. et al., “3-D true-amplitude finite-offset migration,” Geophysics,volume 58, pp. 1112-1126 (1993), the weight w may be chosen in such a way that the migration processing preserves the seismic amplitudes. This is usually considered very desirable, since the amplitudes are widely used in interpreting the processed data.

[0016] The Kirchhoff prestack migration process may be expressed mathematically in the form of equation (1): $\begin{array}{cc}{D}_M\left(\overrightarrow{x}\right)={\left[\frac{\partial }{\partial t}\int d^2\overrightarrow{\xi }w\left(\overrightarrow{\eta },z\right)D\left(t+t_D,\overrightarrow{\xi }\right)\right]}_{t=0}& \left(1\right)\end{array}$

[0017] where D_M({right arrow over (x)}) is the migrated image at the point {right arrow over (x)}, D(t,{right arrow over (ξ)}) denotes the input data at time t and midpoint {right arrow over (ξ)}, and w({right arrow over (η)},z) is the weight function discussed above, expressed as a function of {right arrow over (η)} and z (depth). The variable {right arrow over (η)} in equation (1) is the displacement between the input trace midpoint and the surface coordinates of the image point, {right arrow over (η)}={right arrow over (x)}−{right arrow over (ξ)}. Since, as explained above, offset is held fixed during prestack migration, there is no reference to it in equation (1) or any of the other equations set forth herein.

[0018] A sketch of the geometry is included in FIG. 1. A seismic source 10 and a seismic receiver 12 are located on the surface of the earth 14. The midpoint between source 10 and receiver 12 is indicated by {right arrow over (ξ)}. Incident seismic ray 16 encounters subsurface reflector 18 at reflection point {right arrow over (x)}, and the resulting reflected seismic ray 20 is received by seismic receiver 12. The seismic signal traveltime from source 10 to reflection point {right arrow over (x)} is t₁. The seismic signal traveltime from reflection point {right arrow over (x)} to receiver 12 is t₂. The total seismic signal traveltime is denoted by t_D=t₁+t₂. As noted above, the horizontal displacement between midpoint {right arrow over (ξ)} and the surface (i.e., x and y) coordinates of reflection point {right arrow over (x)} is represented by {right arrow over (η)}.

[0019] In equation (1), the initial partial derivative indicates that the data should be differentiated with respect to time in order to preserve the shape of the seismic pulse during processing. All quantities in equation (1) are taken to be at fixed source-receiver offset. As more fully discussed below, a difficulty with implementation of equation (1) is that the summation over {right arrow over (ξ)} generally requires time-consuming computation.

[0020] Equation (1) operates in the space-time ({right arrow over (ξ)}−t) domain. For the special situation of zero-offset migration, it is well known that equation (1) can be implemented more efficiently by operating in the frequency-wavenumber (ω−{right arrow over (k)}) domain. This is described, for example, in Kim et al., “Recursive wavenumber-frequency migration,” Geophysics,volume 54, pp. 319-329 (1989).

[0021] In Dubrulle, “Numerical methods for the migration of constant-offset sections in homogeneous and horizontally layered media,” Geophysics,volume 48, pp. 1195-1203 (1983), the author hypothesized a possible approach to ω−{right arrow over (k)} migration at finite offset, and demonstrated his method applied to a simple example for 2-D data. However, Dubrulle did not discuss how to implement his method for 3-D data. In addition, Dubrulle did not address the issue of how to implement ω−{right arrow over (k)} migration to preserve seismic amplitudes, nor did he address the issue of how to handle seismic turning rays. The present invention overcomes these limitations.

[0022] A recent example of a prestack time migration method which operates in the frequency-wavenumber domain and has similar accuracy to Kirchhoff prestack migration is described in Etgen, “V(Z) F-K migration of common-offset common-azimuth data volumes, part I: theory,” 68th Annual International Meeting of the Society of Exploration Geophysicists, Expanded Abstracts (Sep. 13-18, 1998). According to Etgen, accurate prestack time migration can be performed using multiplication in the ω−{right arrow over (k)} domain of a migration operator with the data volume.

[0023] The method proposed by Etgen is shown pictorially in FIG. 2. The unmigrated common-offset data traces D(t,{right arrow over (ξ)}) are 3-D Fourier transformed from the space-time ({right arrow over (ξ)}−t) domain to the frequency-wavenumber (ω−{right arrow over (k)}) domain resulting in transformed data {tilde over (D)}(ω,{right arrow over (k)}). In FIG. 2, this is represented by the transformation from data cube 22 to data cube 24. The migration operator M(ω,{right arrow over (η)},z) is first computed in the ω−{right arrow over (η)} domain for every depth z₁, z₂, . . . z_N. In FIG. 2, the migration operator is represented by cubes 26, 28, and 30 for depths z₁, z₂, and z_N, respectively. Migration operator M(ω,{right arrow over (η)},z) is then 2-D Fourier transformed for every frequency ω and every depth z into the frequency-wavenumber (ω−{right arrow over (k)}) domain resulting in transformed migration operator {tilde over (M)}(ω,{right arrow over (k)},z). In FIG. 2, the transformed migration operator is represented by cubes 32, 34, and 36 for depths z₁, z₂, and z_N, respectively. The transformed migration operator {tilde over (M)}(ω,{right arrow over (k)},z) for each depth is then multiplied by the transformed data {tilde over (D)}(ω,{right arrow over (k)}). For each depth z₁, z₂, . . . z_N, the migrated data are then summed over frequency to form an image {tilde over (D)}_M({right arrow over (k)},z) in the wavenumber domain. In FIG. 2, this is represented by images 38, 40, and 42 for depths z₁, z₂, and z_N, respectively. Finally, the wavenumber domain migrated images {tilde over (D)}_M({right arrow over (k)},z) are 2-D Fourier transformed from the wavenumber domain back to the space domain resulting in migrated data D_M({right arrow over (x)},z), which is represented in FIG. 2 by migrated data cube 44.

[0024] The method described by Etgen is substantially faster than Kirchhoff prestack migration according to equation (1), but is dominated by the costs of computing and transforming the migration operator repeatedly for each depth. As further described below, in some embodiments of the present invention the unmigrated data traces are transformed to the ω−{right arrow over (k)} domain in the same manner as proposed by Etgen; however, the present invention computes and applies the migration operator in a unique way that requires significantly less computer resources than the Etgen method.

[0025] From the foregoing, it can be seen that there is a need for a method for economically achieving accurate prestack migration of seismic data. Such a method should be capable of accurately preserving the seismic amplitudes of the unmigrated data. Such a method should also be capable of imaging reflectors having very steep dips, including reflectors having dips above 90 degrees dip angle such as may be found, for example, adjacent to salt structures in the Gulf of Mexico. The present invention satisfies this need.

[0026] In one embodiment, the present invention is a method for prestack time migration of common-offset seismic data traces which comprises the steps of (1) selecting a seismic velocity function for use in the prestack time migration, (2) transforming the common-offset seismic data traces from the space-time ({right arrow over (ξ)}−t) domain to a preselected alternate data domain, (3) using the seismic velocity function to compute a time function in the {right arrow over (p)}−z domain for accurately migrating the seismic data traces, (4) computing a weighting function in the {right arrow over (p)}−z domain, which weighting function substantially preserves the seismic amplitudes of the seismic data traces, (5) using the time function to phase shift the seismic data traces in the alternate data domain, (6) multiplying the amplitudes of the phase-shifted seismic data traces by the weighting function, (7) summing the phase-shifted and amplitude-weighted data over an appropriately chosen variable so as to form an image in said alternate data domain, and (8) performing a 2-D transformation of the image from the alternate data domain back to the space ({right arrow over (x)}−z) domain.

[0027] The preselected alternate data domain may be the ω−{right arrow over (k)} domain, the {right arrow over (k)}−t domain, or the τ−{right arrow over (p)} domain.

[0028] The inventive method also permits imaging of seismic turning ray energy for non-zero source-receiver offsets.

[0029] The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

[0030] FIG. 1 is a sketch of the basic geometry used in common-offset prestack time migration according to the present invention;

[0031] FIG. 2 is a pictorial representation of the method of prestack time migration in the ω−{right arrow over (k)} domain proposed by Etgen;

[0032] FIG. 3 is a flow chart illustrating the basic steps of prestack time migration in the ω−{right arrow over (k)} domain according to one embodiment of the present invention;

[0033] FIG. 4 is a pictorial representation of prestack time migration in the ω−{right arrow over (k)} domain according to the present invention;

[0034] FIG. 5 illustrates the resampling from {right arrow over (η)} space to {right arrow over (p)} space used in calculating the migration operator according to the present invention;

[0035] FIG. 6 shows the contours of the migration time shift function T at fixed depth and offset in {right arrow over (η)} space;

[0036] FIG. 7 is a 3-D rendering of the migration time shift function T at fixed depth and offset in {right arrow over (η)} space;

[0037] FIG. 8 shows the contours of the K_3Dfactor in {right arrow over (η)} space; and

[0038] FIG. 9 is a 3-D rendering of the migration time shift function T after it has been converted from {right arrow over (η)} space to {right arrow over (p)} space.

[0039] The invention will be described in connection with its preferred embodiments. However, to the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications, and equivalents that are included within the spirit and scope of the invention, as defined by the appended claims.

[0040] The present invention is a method for efficiently achieving rapid, controlled-amplitude prestack time migration of seismic data. According to the inventive method, common-offset, common-azimuth gathers of seismic data traces (hereinafter referred to as “common-offset data gathers”) are transformed from the space-time ({right arrow over (ξ)}−t) domain to a preselected alternate data domain where they are migrated using a unique migration operator which is computed in the {right arrow over (p)}−z domain. The migrated data traces are then transformed back to the space domain. Examples of appropriate alternate data domains for use in connection with the present invention are the frequency-wavenumber (ω−{right arrow over (k)}) domain, the wavenumber-time ({right arrow over (k)}−t) domain, and the slant-stack (τ−{right arrow over (p)}) domain. In the following detailed description, the invention will be described with respect primarily to implementation in the ω−{right arrow over (k)} domain. However, it should be understood that this is for illustrative purposes only and does not constitute a limitation on the scope of the invention which also may be implemented advantageously in the {right arrow over (k)}−t or τ−{right arrow over (p)} domains.

[0041] Preferably, the inventive method is practiced using a suitably programmed digital computer or parallel processor. Persons skilled in the art of seismic data processing could easily develop computer software for implementing the inventive method based on the teachings set forth herein.

[0042] According to the present invention, controlled-amplitude prestack time migration in the ω−{right arrow over (k)} domain is implemented using the following migration equation: $\begin{array}{cc}{D}_M\left(\overrightarrow{x}\right)=\int d^2\overrightarrow{k}{e}^{-i\overrightarrow{k}·\overrightarrow{x}}\int d{\mathrm{\omega e}}^{\mathrm{i\omega }T\left(\overrightarrow{\eta },z\right)}\tilde{D}\left(\omega ,\overrightarrow{k}\right)\frac{w\left(\overrightarrow{\eta },z\right)}{K_{3D}^{1/2}\left(\overrightarrow{\eta },z\right)}& \left(2\right)\end{array}$

[0043] where D_M({right arrow over (x)}) is the migrated image at point {right arrow over (x)}, {tilde over (D)}(ω,{right arrow over (k)}) represents the input data after Fourier transformation with respect to time and midpoint so as to transform the data to the ω−{right arrow over (k)} domain, T({right arrow over (η)},z) is the migration time shift function, and w({right arrow over (η)},z)/K_3d^1/2({right arrow over (η)},z) is an amplitude-preserving weight function.

[0044] FIG. 3 is a flowchart illustrating the principal steps used in performing controlled amplitude prestack time migration in the ω−{right arrow over (k)} domain according to equation (2). At step 46, a set of common-offset seismic data traces is obtained. Initially, these data traces are in the space-time ({right arrow over (ξ)}−t) domain. Next, at step 48, a seismic velocity function for the subsurface region in question is selected. As noted above, a basic assumption of the present invention is that the subsurface seismic velocity is a function of depth or time, but is laterally invariant. Lateral velocity variation, for example due to a sloping sea floor, can be taken into account by first applying NMO with a selected local velocity and then applying a reverse NMO with the velocity used subsequently to migrate the data. Techniques for selecting an appropriate seismic velocity function for use in the present invention are well known to persons skilled in the art and, accordingly, will not be described further herein.

[0045] At step 50, the common-offset data traces obtained at step 46 are 3-D Fourier transformed from the {right arrow over (ξ)}−t domain to the ω−{right arrow over (k)} domain. In equation (2), the Fourier transformed data traces are denoted by {tilde over (D)}(ω,{right arrow over (k)}).

[0046] At step 52, the seismic velocity function selected at step 48 is used to compute a time function for use in migrating the transformed data traces in the ω−{right arrow over (k)} domain. This time function is denoted by T({right arrow over (η)}, z) in equation (2) and is defined by:

[0047] T({right arrow over (η)},z)=t_D−{right arrow over (p)}·{right arrow over (η)}  (3)

[0048] where $\begin{array}{cc}\overrightarrow{p}={\left(\frac{\partial t_D}{\partial \overrightarrow{\eta }}\right)}_z& \left(4\right)\end{array}$

[0049] is the midpoint ray parameter computed at fixed depth. As is well known to those skilled in the art, the ray parameter {right arrow over (p)} is the horizontal component of the slowness vector. As further described below, the time function T({right arrow over (η)},z) is converted to the {right arrow over (p)}−z (dip-depth) domain before it is applied to the transformed data traces.

[0050] At step 54, a weighting function for preserving the seismic amplitudes of the input data traces is computed. As noted above, the weighting function is denoted by w({right arrow over (η)},z)/K_3D^1/2({right arrow over (η)},z) in equation (2).

[0051] w({right arrow over (η)},z) is the weight function described by Schleicher et al., referred to above. According to Schleicher et al., the weight function depends on {right arrow over (ξ)} and {right arrow over (x)}, and is given by: $\begin{array}{cc}w\left(\overrightarrow{\xi },\overrightarrow{x}\right)=\frac{\sqrt{\cos {\alpha }_S\cos {\alpha }_G}}{v_S}\frac{\left|\det \left(N_{\mathrm{SR}}+N_{\mathrm{RG}}\right)\right|}{\sqrt{\left|\det N_{\mathrm{SR}}\right|}\sqrt{\left|\det N_{\mathrm{RG}}\right|}}& \left(5\right)\end{array}$

[0052] where α_sis the surface dip angle of the ray connecting the source {right arrow over (s)} and the reflection point {right arrow over (x)}, α_Gis the surface dip angle of the ray connecting the receiver {right arrow over (g)} and the reflection point {right arrow over (x)}, v_sis the surface compressional wave velocity, and the matrices are given by: $\begin{array}{cc}{N}_{\mathrm{SR}}=\left(\begin{array}{cc}\frac{{\partial }^2{t}_1}{\partial {\xi }_x\partial \rho }& \frac{{\partial }^2{t}_1}{\partial {\xi }_y\partial \rho }\\ \frac{{\partial }^2{t}_1}{\partial {\xi }_x\partial \sigma }& \frac{{\partial }^2{t}_1}{\partial {\xi }_y\partial \sigma }\end{array}\right)& \text{(6a)}\\ \mathrm{and},& \\ N_{\mathrm{RG}}=\left(\begin{array}{cc}\frac{{\partial }^2{t}_2}{\partial {\xi }_x\partial \rho }& \frac{{\partial }^2{t}_2}{\partial {\xi }_y\partial \rho }\\ \frac{{\partial }^2{t}_2}{\partial {\xi }_x\partial \sigma }& \frac{{\partial }^2{t}_2}{\partial {\xi }_y\partial \sigma }\end{array}\right)& \text{(6b)}\end{array}$

[0053] where ρ and σ are coordinates on a plane tangential to the reflector at the reflection point R. Given a knowledge (for example, from a time map) of the rays connecting the source and receiver to the reflection point, all of the above derivatives can be calculated. The source and receiver locations are {right arrow over (s)}={right arrow over (ξ)}+{right arrow over (h)} and {right arrow over (g)}={right arrow over (ξ)}−{right arrow over (h)}, where {right arrow over (h)} is the half-offset as a two-dimensional vector on the earth's surface. Since {right arrow over (η)}={right arrow over (x)}−{right arrow over (ξ)}, the weight function w({right arrow over (ξ)},{right arrow over (x)}) may be expressed as a function of {right arrow over (ξ)} and {right arrow over (x)} or as a function of {right arrow over (η)} and z. Alternatively, the weight function w({right arrow over (η)},z) can be computed using dynamic ray tracing, as described in Hanitzsch, “Comparison of weights in prestack amplitude preserving Kirchhoff depth migration,” Geophysics,volume 62, pp. 1812-1816 (1997).

[0054] The factor K_3D({right arrow over (η)},z) is a normalization factor used to adjust for differences between working in {right arrow over (x)} space and working in {right arrow over (k)} space. K_3D({right arrow over (η)},z) is defined by the following equation: $\begin{array}{cc}{K}_{3D}\left(\overrightarrow{\eta },z\right)\det \left[\frac{{\partial }^2{t}_D}{\partial {\eta }_i\partial {\eta }_j}\right]& \left(7\right)\end{array}$

[0055] where the indices i and j take on the values x and y corresponding to the two components of the variable {right arrow over (η)}.

[0056] As with the time function T({right arrow over (η)},z), the amplitude-preserving weighting function w({right arrow over (η)},z)/K_3D^1/2({right arrow over (η)},z) is converted to the {right arrow over (p)}−z (dip-depth) domain before it is applied to the transformed data traces. This conversion is described in detail below.

[0057] At step 56, the transformed data traces {tilde over (D)}(ω,{right arrow over (k)}) are phase shifted using the time function calculated at step 52.

[0058] Next, at step 58, the phase-shifted data traces are multiplied by the amplitude-preserving weighting function calculated at step 54.

[0059] Finally, at steps 60 and 61, the phase-shifted and weighted data traces are transformed from the ω−{right arrow over (k)} domain back to the {right arrow over (x)}−z domain. This is accomplished in two stages. First, at step 60, an image is formed in the wavenumber ({right arrow over (k)}) domain by summing with respect to frequency ω. Then, at step 61, the {right arrow over (k)} domain image is 2-D Fourier transformed back to the {right arrow over (x)} domain. The result of this process is {right arrow over (x)}−z domain migrated data traces in which the migration process has preserved seismic amplitudes.

[0060] The functions T and w/K_3D^1/2are together referred to as the “migration operator,” {right arrow over (M)}=(T,w/K_3D^1/2). A fundamental difficulty in the implementation of equation (2) arises because the functions T and w/K_3D^1/2are initially calculated in {right arrow over (η)}−z space while the transformed data traces are in ω−{right arrow over (k)} space. To overcome this difficulty, T and w/K_3D^1/2are first converted to the {right arrow over (p)}−z domain, as further described below. The conversion of T and w/K_3D^1/2to {right arrow over (p)}−z space permits an efficient implementation of equation (2) since points in {right arrow over (p)}−z space can be related to points in ω−{right arrow over (k)} space using the following relationship:

[0061] {right arrow over (p)}=ω/{right arrow over (k)}.  (8)

[0062] Implementation of the present invention is shown pictorially in FIG. 4. Preferably, the migration operator (i.e., the functions T and w/K_3D^1/2) is first computed in three-dimensional {right arrow over (η)}−z space. In FIG. 4, the three-dimensional migration operator {right arrow over (M)}({right arrow over (η)},z) is represented by cube 62. For every z, the two-dimensional {right arrow over (η)} coordinate is mapped (i.e., transformed) to a two-dimensional {right arrow over (p)} coordinate. This produces a three-dimensional migration operator {right arrow over ({tilde over (M)})}({right arrow over (p)},z) in {right arrow over (p)}−z space, which is represented in FIG. 4 by cube 64. The space-time ({right arrow over (ξ)}−t) domain data D({right arrow over (ξ)},t) are 3-D Fourier transformed into wavenumber-frequency ({right arrow over (k)}−ω) domain data {tilde over (D)}({right arrow over (k)},ω). In FIG. 4, the {right arrow over (ξ)}−t domain data are represented by cube 66, and the {right arrow over (k)}−ω domain data are represented by cube 68. According to the present invention, a single z-slice of the migration operator {right arrow over ({tilde over (M)})}({right arrow over (p)},z) in the {right arrow over (p)}−z domain (e.g., slice 64ior 64j) can be used to migrate all of the data in the ω−{right arrow over (k)} domain to a fixed output depth z by virtue of the relationship in equation (8). The migrated data are then summed with respect to frequency producing a single z-slice (e.g., slice 70ior 70j) of the migrated image {tilde over (D)}_M({right arrow over (k)},z), which is represented in FIG. 4 by cube 70. This is repeated for every z-slice in the migration operator. Finally, the migrated image is 2-D Fourier transformed to the {right arrow over (x)}−z domain, represented by cube 72. In effect, the present invention reduces the problem of computing and mapping the operator from a four-dimensional problem to a three-dimensional problem. (Compare Etgen's four-dimensional migration operator {tilde over (M)}(ω,{right arrow over (k)},z) with the three-dimensional migration operator {right arrow over ({tilde over (M)})}({right arrow over (p)},z) of the present invention.) Thus, the present invention greatly increases the computational efficiency of the process.

[0063] The imaging process is computationally most efficient if the functions T and w/K_3D^1/2are precomputed for a regularly sampled grid of {right arrow over (p)} values to be used in a table look-up procedure during the implementation of equation (2). The mapping of the two components of the migration operator, T and w/K_3D^1/2, from {right arrow over (η)} space to {right arrow over (p)} space greatly improves the efficiency of this invention relative to the technique described by Etgen.

[0064] Equation (4) relates {right arrow over (p)} space to {right arrow over (η)} space. In practice, the migration operator is initially computed for a regularly sampled grid in {right arrow over (η)} space at fixed depth z. Numerical differentiation of t_Daccording to equation (4) produces {right arrow over (p)}={right arrow over (p)}({right arrow over (η)}) from which T can be assembled according to equation (3). The function K_3Dcan also be constructed from numerical derivatives of t_Daccording to equation (7). The factors necessary to compute the weight function w can also be constructed from numerical derivatives of t_D, as described in Schleicher at al., or they can be computed using dynamic ray tracing, as described by Hanitzsch. For every depth z, the migration operator is then resampled from the regular grid in {right arrow over (η)} space to a regular grid in {right arrow over (p)} space using the computed samples of the function {right arrow over (p)}={right arrow over (p)}({right arrow over (η)}), as described in the following paragraph.

[0065] The preferred method of resampling from {right arrow over (η)} space to {right arrow over (p)} space is illustrated in FIG. 5 and is accomplished in the following manner. The Cartesian cells of the regular grid in {right arrow over (η)} space (e.g., the rectangular cell formed by grid points 74a,74b,74c,and 74d) are further decomposed into triangles as shown in FIG. 5, and the values of {right arrow over (p)} are computed at the corners of the triangles. Either diagonal may be used when decomposing the rectangle into two triangles. Within each triangle, the function {right arrow over (p)}={right arrow over (p)}({right arrow over (η)}) may be represented by an arbitrary interpolating function denoted {right arrow over (p)}={right arrow over (p)}({right arrow over (η)}), whose parameters are evaluated using the precomputed values of {right arrow over (p)} at the corners of the triangle. Using the interpolating function, the minimum and maximum values of {right arrow over (p)} within the triangle can be easily computed. Knowing the range of {right arrow over (p)}, one can deduce whether there potentially exist any points {right arrow over (p)}*={right arrow over ({tilde over (p)})}({right arrow over (η)}*) within the triangle (e.g., point 76) such that {right arrow over (p)}* is on the desired regular grid in {right arrow over (p)} space. For example, assume that p_x*=const₁and p_y*=const₂is a point on the desired regular grid in {right arrow over (p)} space. Contour 78 (FIG. 5) represents the points on the interpolating function for which const₁={tilde over (p)}_x({right arrow over (η)}). Similarly, contour 80 represents the points on the interpolating function for which const₂={tilde over (p)}_y({right arrow over (η)}) The intersection of contours 78 and 80 in {right arrow over (η)} space represents the point (p_x*,p_y*) in {right arrow over (p)} space. If any such points exist within the triangle, the value of {right arrow over (η)}* can be computed by inverting the functional representation {right arrow over (p)}*={right arrow over ({tilde over (p)})}({right arrow over (η)}*). A linear interpolating function is preferred, as inverting the linear function amounts to solving a system of linear equations. Higher order polynomials can be used if desired, in which case non-linear methods yield {right arrow over (η)}*. Once a point {right arrow over (η)}* has been computed, it is advisable to double-check and verify that its coordinates are within the boundary of the triangle. For each point {right arrow over (η)}* within the triangle, the migration operator can be interpolated from the precomputed regular grid in {right arrow over (η)} space at the point {right arrow over (η)}* to yield T({right arrow over (p)}*) and w({right arrow over (p)}*)/K_3D^1/2({right arrow over (p)}*), the desired migration operator at the point {right arrow over (p)}*. This process is repeated for every triangle in {right arrow over (η)} space. If desired, it is possible to use other decompositions of the original {right arrow over (η)} grid other than triangles to accomplish the remapping to {right arrow over (p)} space. For example, it is possible to use the original rectangular grid in conjunction with a bilinear interpolating function.

[0066] As described above, the present invention can be implemented on a conventional digital computer or on a parallel processor. When utilizing a parallel processor, there are a number of options for parallel implementation of equation (2). Equation (2) can be parallelized over depth, frequency, or wavenumber. The preferred method of parallelization is to first share the work of 3-D Fourier transforming the data across the processors according to temporal frequency and then to divide the work of computing and applying the operator across the processors according to depth. Each processor computes the operator for a range of depths and applies it over the same range of depths. This has the advantage that it not only distributes the work of computing and applying the operator, but also divides the memory requirement for holding the operator and image cubes across the processors. When parallelizing in this manner, each processor must access the entire transformed data volume. Because of the size of typical data volumes, it is usually not possible to hold an independent copy in the memory of each processor. To overcome this limitation, the data cube is stored on disk after it has been transformed to the frequency-wavenumber domain. One processor is assigned to asynchronously read frequency slices of the data cube from the disk and asychronously distribute the frequency slices to the other processors while the imaging is taking place. In this algorithm, the processor which has been assigned the task of reading and distributing frequency slices from the disk file (referred to as the I/O processor) is distributing the current frequency slice to all the other processors at the same time that all the other processors are migrating the previous frequency slice. Simultaneous with this, the I/O processor is reading the next frequency slice from the disk. This scheme assures that none of the processors are ever waiting for the data to come off the disk and into their memory. The processors performing the migration (referred to as the migration processors) apply the migration operator to each frequency slice as they receive it. Each migration processor applies the migration operator at every depth which has been assigned to it. The image is formed by summing the migrated frequency slice into an accumulation buffer for each depth. Once all the frequency slices from the input file have been migrated, the migration processors inverse Fourier transform the image from the wavenumber domain to the space domain and then output their imaged depth slices.

[0067] When seismic velocity increases with depth in the earth, seismic rays that initially travel downward can turn and head back toward the surface without being reflected from a subsurface reflector. These “turning rays” can then reflect off the underside of overturned reflectors and return to the surface by a path that initially heads downward and then turns upward. Turning rays can be used to image overturned reflectors. Hale et al. (“Imaging salt with turning seismic waves,” Geophysics,volume 57, pp. 1453-1462 (1992)) discuss an ω−{right arrow over (k)} migration algorithm for imaging turning rays for zero source-receiver separation. The present invention permits imaging of turning ray energy for non-zero source-receiver separation.

[0068] In general, turning rays are recorded at large distances from the overturned reflectors. The algorithm described herein can be used to image turning rays by simply computing the migration operator over a sufficiently large range of {right arrow over (η)} space to include any turning energy. When the migration operator includes turning ray energy, it becomes double-valued in {right arrow over (p)} space. The migration operator has a primary sheet corresponding to energy which does not turn and a secondary sheet corresponding to turning rays. These sheets are easily separated during the resampling process (i.e., the conversion from {right arrow over (η)} space to {right arrow over (p)} space) by taking advantage of the fact that the sign of K_3Dchanges from positive to negative when moving from the primary sheet to the secondary sheet. This is illustrated in FIGS. 6-9. FIG. 6 shows contours of the migration time shift function T at fixed depth and offset in {right arrow over (η)} space where the function is single-valued, and FIG. 7 is a 3-D rendering of the same migration time shift function T. Note that for any point (η_x,η_y) in {right arrow over (η)} space, there is only a single value of the migration time shift function T. FIG. 8 shows contours of K_3Din {right arrow over (η)} space. The zero contour 82 is elliptical in shape. Inside this ellipse, K_3Dis positive. Outside, it is negative. This fact can be used when mapping the operator from {right arrow over (η)} space to {right arrow over (p)} space. The primary sheet of the operator is in the region where K_3Dis positive. The secondary sheet of the operator corresponding to turning rays is the region where K_3Dis negative. FIG. 9 is a 3-D rendering of the migration time shift function T after it has been converted from {right arrow over (η)} space to {right arrow over (p)} space clearly showing the double-valued nature of the function as it is folded over on itself.

[0069] When imaging turning rays, equation (2) is computed twice, once for each sheet in the migration operator, and the results are summed. When migrating with the second sheet of the migration operator, the right hand side of equation (2) must be modified by an extra factor of e^−iπ/2in order to have the proper waveform phase on the image of the overturned reflectors.

[0070] The computational advantages of the migration technique described above with respect to Kirchhoff prestack migration can be assessed by counting the number of floating point operations required for implementation of equation (2) versus the number required for implementation of equation (1). In equation (2), there are three steps. First, the data traces must be 3-D Fourier transformed from the {right arrow over (ξ)}−t domain to the ω−{right arrow over (k)} domain. If the number of data samples in the x, y, and t dimensions are N_x, N_y, and N_t, respectively, then this initial forward 3-D Fourier transformation requires on the order of N_tN_xN_yIn(N_tN_xN_y) operations. In the second step, the data are imaged by applying the migration operator and performing a summation over frequency. This step must be performed for each image depth and converts the data from ω−{right arrow over (k)} space to {right arrow over (k)}−z space. If the number of imaged depths in N_z, then the total number of operations in this second step is proportional to N_zN_xN_yN_t. In the third step, the data are transformed from {right arrow over (k)}−z space to {right arrow over (x)}−z space by performing a 2-D Fourier transformation for each depth. The number of operations required for this third step is proportional to N_zN_xN_yIn(N_xN_y).

[0071] The Kirchhoff prestack migration method given in equation (1) requires a 2-D summation over midpoint for each of the image points. Thus, the total number of operations is proportional to N_x²N_y²N_z.

[0072] To compare the relative speeds of the two methods, note that for typical problems N_z≡N_x≡N_y≡N_t≡1000. Because all of the variables are the same order of magnitude, they can be replaced by the generic variable N. The time required for performing the ω−{right arrow over (k)} migration according to equation (2) is controlled by the second step which is proportional to N⁴. Kirchhoff prestack migration according to equation (1) is proportional to N⁵. Thus, prestack migration according to the present invention is a full power of N faster than Kirchhoff prestack migration according to equation (1).

[0073] As noted above, and illustrated in FIG. 2, the method of Etgen requires a 2-D Fourier transformation of the operator from the ω−{right arrow over (η)} domain to the ω−{right arrow over (k)} domain at each value of ω and z. This means that the method requires on the order of N⁴In(N²) operations to construct the migration operator. This means that the Etgen method is slower than the method presented herein by a factor of In(N²) operations. Since N is on the order of 1000, this amounts to a factor of about In(10⁶) or about an order of magnitude.

[0074] As noted above, the present invention may also be implemented in other alternate data domains. These alternate data domain implementations of the invention have the same computational advantages as discussed above with respect to implementation in the ω−{right arrow over (k)} domain.

[0075] Implementation in {right arrow over (k)}−t Space

[0076] By methods similar to those described above for ω−{right arrow over (k)} implementation, the following equations can be derived for controlled amplitude prestack time migration in the {right arrow over (k)}−t domain: $\begin{array}{cc}{D}_M\left(\overrightarrow{x},z\right)=\int \int d^2\overrightarrow{k}dt{e}^{-i\overrightarrow{k}·\overrightarrow{x}}\tilde{D}\left(t,\overrightarrow{k}\right){w\left(\frac{\partial }{\partial t}\right)}^{1/2}\frac{e^{i\overrightarrow{k}·\overrightarrow{\eta }}}{{\Delta }^{1/2}}& \left(9\right)\end{array}$

[0077] where now the input data {tilde over (D)}(t,{right arrow over (k)}) are in the {right arrow over (k)}−t domain (and are obtained by a 2-D Fourier transformation from the {right arrow over (ξ)}−t domain), w is the same weight function used in equations (1) and (2), and the denominator Δ is given by $\begin{array}{cc}\Delta \equiv p_2^2{M}_{11}-2p_1{p}_2{M}_{12}+p_1^2{M}_{22},p_i=\frac{\partial t_D}{\partial {\eta }_i}& \left(10\right)\end{array}$

[0078] where $\begin{array}{cc}{M}_{\mathrm{ij}}=\frac{{\partial }^2{t}_D}{\partial {\eta }_i\partial {\eta }_j}.& \left(11\right)\end{array}$

[0079] Preferably, the migration operator is initially calculated in the {right arrow over (η)}−z domain, and is converted to the {right arrow over (p)}−z domain before being applied to the transformed data traces {tilde over (D)}(t,{right arrow over (k)}), as with the previous embodiment.

[0080] Evaluation of equation (9) has the same computational efficiency as does equation (2) in the ω−{right arrow over (k)} domain. In this case, the image is formed in {right arrow over (k)} space by the summation over time t, and the data are subsequently transformed from {right arrow over (k)} space back to {right arrow over (x)} space. Just as in the ω−{right arrow over (k)} domain method, computation time is dominated by the imaging step and takes on the order of N⁴operations.

[0081] Implementation in τ−{right arrow over (p)} Space

[0082] The present invention may also be implemented in the slant-stack or τ−{right arrow over (p)} domain. By setting {right arrow over (k)}={right arrow over (p)}ω, equation (2) may be rewritten as follows: $\begin{array}{cc}{D}_M\left(\overrightarrow{x},z\right)=-\frac{{\partial }^2}{\partial T^2}\int d^2\overrightarrow{p}\tilde{D}\left(T-\overrightarrow{p}.\overrightarrow{x},\overrightarrow{p}\right)\frac{w}{K_{3D}^{1/2}}& \left(12\right)\end{array}$

[0083] The input data {tilde over (D)}(T−{right arrow over (p)}.{right arrow over (x)},{right arrow over (p)}) on the right hand side of equation (12) are in the τ−{right arrow over (p)} domain, with τ=T({right arrow over (p)},z)−{right arrow over (p)}.{right arrow over (x)}. w/K_3D^1/2is the same amplitude preserving weight function discussed above with respect to equation (2). Evaluation of equation (12) appears to be proportional to N⁵operations. However, as in the previous cases, the summation can be rewritten such that implementation is proportional to N⁴operations. This can be done by making a 2-D Radon transform to Radon transform coordinates (P,{circumflex over (ζ)}) to obtain the following equation: $\begin{array}{cc}{\tilde{D}}_M\left(P,\hat{\zeta },z\right)=\frac{\partial }{\partial P}\int du\tilde{D}\left(\mathrm{uP}+T\left(u\hat{\zeta }\right),u\hat{\zeta }\right)·\frac{w}{K_{3D}^{1/2}}& \left(13\right)\end{array}$

[0084] Equation (13) now represents an imaging condition where the summation over variable u (the Fourier conjugate of P) is analogous to the imaging summations over ω and t discussed above. The number of operations required to implement equation (13) is proportional to N⁴. The inverse Radon transform required to return from P−{circumflex over (ζ)} space to {right arrow over (x)}−z space can be evaluated as a product of two fast Fourier transforms and is analogous to the previously discussed final transforms to return from {right arrow over (k)} space to {right arrow over (x)} space. As with the ω−{right arrow over (k)} and {right arrow over (k)}−t embodiments, the migration operator preferably is initially computed in the {right arrow over (η)}−z domain, and is converted to the {right arrow over (p)}−z domain before being applied to the transformed data traces {tilde over (D)}(uP+T(u{circumflex over (ζ)})u{circumflex over (ζ)}). The data are obtained in the τ−{right arrow over (p)} domain by first transforming to the ω−{right arrow over (k)} domain with a 3-D Fourier transformation and then transforming to the τ−{right arrow over (p)} domain by a 1-D Fourier transformation at fixed {right arrow over (p)}={right arrow over (k)}/ω as in equation (8). Both of these steps require on the order of N³In(N) operations.

[0085] A description of τ−{right arrow over (p)} migration is given in Hubral, P., Slant Stack Processing,pp. 111-117, edited by G. Gardner and L. Lu, published by the Society of Exploration Geophysicists (1991). Hubral, however, confined his discussion to migration of poststack data, and also did not discuss either how to make the process amplitude preserving or how to optimize the computational efficiency of the method.

[0086] The foregoing description has been directed to particular embodiments of the invention for the purpose of illustrating the invention. It will be apparent to persons skilled in the art, however, that many modifications and variations to the embodiments described herein are possible. For example, in any of the foregoing embodiments, the final migrated image could easily be output with respect to vertical time (e.g., in the {right arrow over (x)}−t domain) rather than with respect to vertical depth (the {right arrow over (x)}−z domain). All such modifications and variations are intended to be within the scope of the present invention, as defined by the appended claims.