Modeling angle domain common image gathers from reverse time migration

The present disclosure describes methods and systems, including computer-implemented methods, computer program products, and computer systems, for modeling angle domain common image gathers (ADCIG) from reverse time migration (RTM). One computer-implemented method includes calculating seismic source and receiver wavefields based on seismic data, calculating characteristic source and receiver wavefields from the seismic source and receiver wavefields, calculating propagation angles for the characteristic source and receiver wavefields, applying a wavefield decomposition algorithm on the characteristic source and receiver wavefields to obtain corresponding directional source and receiver wavefields, the wavefield decomposition algorithm decomposing wavefield amplitude of a wavefield in an angle interval centered on a propagation angle of the wavefield, and forming ADCIG by applying an image condition to the obtained directional source and receiver wavefields.

TECHNICAL FIELD

This disclosure relates to seismic data processing and, more specifically, to modeling angle domain common image gathers from reverse time migration.

BACKGROUND

High resolution depth images of earth subsurface layers are important for reservoir exploration, delineation, and development. Wave equation based seismic depth migration techniques, such as reverse time migration (RTM), are suitable techniques for seismic depth imaging in the oil industry, especially in complex environments such as subsalt exploration. However, RTM is a computationally intensive process that requires propagating waves in 2D or 3D models using time or frequency domain wave equation solvers.

SUMMARY

The present disclosure describes methods and systems, including computer-implemented methods, computer program products, and computer systems for modeling angle domain common image gathers (ADCIG) from reverse time migration (RTM). One computer-implemented method for modeling ADCIG from RTM includes calculating seismic source and receiver wavefields based on seismic data, calculating characteristic source and receiver wavefields from the seismic source and receiver wavefields, calculating propagation angles for the characteristic source and receiver wavefields, applying a wavefield decomposition algorithm on the characteristic source and receiver wavefields to obtain corresponding directional source and receiver wavefields, the wavefield decomposition algorithm decomposing wavefield amplitude of a wavefield in an angle interval centered on a propagation angle of the wavefield, and forming ADCIG by applying an image condition to the obtained directional source and receiver wavefields.

Other implementations of this aspect include corresponding computer systems, apparatuses, and computer programs recorded on one or more computer storage devices, each configured to perform the actions of the methods. A system of one or more computers can be configured to perform particular operations or actions by virtue of having software, firmware, hardware, or a combination of software, firmware, or hardware installed on the system that in operation causes the system to perform the actions. One or more computer programs can be configured to perform particular operations or actions by virtue of including instructions that, when executed by data processing apparatus, cause the apparatus to perform the actions.

The foregoing and other implementations can each, optionally, include one or more of the following features, alone or in combination:

A first aspect, combinable with the general implementation, wherein the seismic source and receiver wavefields are calculated using a time-domain wave equation solver.

A second aspect, combinable with any of the previous aspects, wherein the seismic source and receiver wavefields are calculated using a frequency-domain wave equation solver.

A third aspect, combinable with any of the previous aspects, wherein the ADCIG are formed using a RTM process.

A fourth aspect, combinable with any of the previous aspects, wherein the characteristic source and receiver wavefields are calculated for a plurality of predefined time steps with a first-order acoustic wave equation propagating model.

A fifth aspect, combinable with any of the previous aspects, wherein the characteristic source and receiver wavefields are calculated for a plurality of predefined time steps with a second-order acoustic wave equation propagating model.

A sixth aspect, combinable with any of the previous aspects, wherein the characteristic source wavefield has an absolute amplitude, at the propagation angle for the characteristic source wavefield, larger than an absolute amplitude at any other angles, and wherein the characteristic receiver wavefield has an absolute amplitude, at the propagation angle for the characteristic receiver wavefield, larger than an absolute amplitude at any other angles.

A seventh aspect, combinable with any of the previous aspects, wherein the angle interval is 1 degree, 5 degrees, or 15 degrees.

An eighth aspect, combinable with any of the previous aspects, wherein the image condition is a cross-correlation imaging condition.

A ninth aspect, combinable with any of the previous aspects, wherein the ADCIG are scattering angle gathers (SAG).

A tenth aspect, combinable with any of the previous aspects, wherein the ADCIG are dip angle gathers (DAG).

DETAILED DESCRIPTION

The following detailed description describes modeling angle domain common image gathers (ADCIG) from reverse time migration (RTM) and is presented to enable any person skilled in the art to make and use the disclosed subject matter in the context of one or more particular implementations. Various modifications, alterations, and permutations of the disclosed implementations can be made and will be readily apparent to those skilled in the art, and the general principles defined may be applied to other implementations and applications without departing from scope of the disclosure. Thus, the present disclosure is not intended to be limited to the described or illustrated implementations, but is to be accorded the widest scope consistent with the principles and features disclosed.

High resolution depth images of earth subsurface layers are important for successful reservoir exploration, delineation, and development. Wave equation based seismic depth migration techniques, such as RTM, are suitable techniques for seismic depth imaging in the oil industry, especially in complex environments such as subsalt exploration. However, RTM is a computationally intensive process that requires propagating waves in 2D or 3D models using time or frequency domain wave equation solvers.

Given a model of the subsurface and seismic data recorded at known receiver locations, an image of the subsurface can be formed by combining, with an imaging condition, the forward propagated source wavefield with the backward propagated receiver wavefield. When the imaging condition is a function of an additional parameter, common image gathers can be formed. ADCIG are particularly important among various common image gathers, and can be divided into two kinds of gathers, namely, the scattering angle gathers (SAG) and the dip angle gathers (DAG).

In general, two types of methods can be used to compute ADCIG from RTM, namely, local plane wave decomposition based methods and direction vector based methods. Local plane wave decomposition based methods are accurate, but computationally expensive due to the need to perform a plane wave decomposition for each grid point. Direction vector based methods are computationally efficient (e.g., the Poynting vector method (Yoon and Marfurt, 2006)), but can be unstable for complicated wavefields with intersecting events (Patrikeeva and Sava, 2013).

At a high level, the described approach is a direction vector based method. The described ADCIG computation method uses the characteristic directional wavefield of a first order wave system to compute the propagation angle (or propagation direction), decompose the full wavefields in an interval of angle centered on the propagation angle, and compute ADCIG based on the decomposed wavefields. The decomposition over the interval of angle centered on the propagation angle acts as a smoothing process in the angle domain, and can stabilize the decomposition process even for complicated wavefields with intersecting events. As a result, the described ADCIG computation method can produce accurate and computationally affordable ADCIG from RTM.

The example ADCIG computation method, described in the present disclosure, can achieve one or more advantages. First, the computation cost of the ADCIG computation method is low and comparable to the computation cost of the Poynting vector method (discussed in more detail below). Second, the proposed method is robust even in complicated wavefields, since the ADCIG computation method is fully based on quantity computed during the forward modelling. Assuming that the propagation is performed using appropriate grid spacing, the velocity and pressure wavefields are sufficiently accurate to perform the wavefield decomposition in complicated wavefields. Third, the proposed method can be applied to more complex propagation models, such as anisotropic and/or elastic propagation models. Fourth, the proposed method can be applied to Full Waveform Inversion (FWI). Since the ADCIG computation method can compute accurate ADCIG, the method can be applied to a select part of the wavefield (e.g., diving waves, reflections only) that needs to be considered in the FWI. In some applications, the example ADCIG computation method can achieve additional or different advantages.

∂u∂t+A⁢⁢∂u∂x+B⁢⁢∂u∂z=f(1)
where t denotes time, (x, z) denotes the 2-dimensional space coordinate, f is a source term, and u(x, z; t) (which is denoted as u in Equation (1) and below for brevity) stands for the wavefield time-space domain. Although this disclosure refers to the 2D isotropic first-order acoustic wave system for purposes of example, the subject matter of this document can be applied to other types of acoustic wave systems, including the 3D propagation system.

If vx(x, z; t) and vz(x, z; t) (which are denoted respectively as vxand vzbelow for brevity) denote the horizontal and vertical velocity displacement wavefields respectively, p(x, z; t) (which is denoted as p below for brevity) denotes the pressure wavefield (i.e., the wavefield, unless otherwise indicated), and c0and ρ0are the local P-wave velocity and density respectively, the following can be derived from Equation (1):

u=[vxvzp]T,⁢A=-(001ρ0000ρ0⁢c0200)⁢⁢and⁢⁢B=-(000001ρ00ρ0⁢c020)(2)
The matrices A and B are diagonalizable with the same eigenvalues (0, −c0, c0) as:

Four characteristic wavefields can be obtained by decomposing the total pressure wavefield (i.e., p) into the four directions (i.e., the direction of x axis, the opposite direction of x axis, the direction of z axis, the opposite direction of z axis):
wx+=¼(p−ρ0c0vx)  (5)
wx−=¼(p=ρ0c0vx)  (6)
wz+=¼(p−ρ0c0vx)  (7)
wz−=¼(p=ρ0c0vx)  (8)
where p=wx++wx−+wz++wz−.

FIG. 1illustrates example snapshots100of a wavefield and its four characteristic wavefields, according to some implementations. For example, a model with constant P-velocity equal to 4000 m/s (meter/second) is simulated (e.g., wave propagation in a homogenous model). A 15 Hz (hertz) Ricker source is applied in the center location of the 2-dimensional space coordinate (x, z) at t=0 ms (millisecond). The wave propagation is implemented using a staggered grid time domain finite difference solver, based on a first order wave equation.

FIG. 1illustrates, at t=250 ms, a snapshot of the total wavefield102, a snapshot of the directional wavefield wx+104, a snapshot of the directional wavefield wx−106, a snapshot of the directional wavefield wz+108, and a snapshot of the directional wavefield wz−110. The four directional wavefields wx+, wx−, wz+, and wz−are computed using Equations (5) to (8), respectively. For illustration purposes, amplitudes of the wavefields wx+, wx−, wz+, and wz−in snapshots104,106,108, and110, respectively, are multiplied by a factor of 4. As illustrated, the decomposed wavefield amplitudes are maximum or above a predefined threshold in the main direction of propagation from the center location, compared to amplitudes in other directions. For example, in the snapshot104, the main direction of propagation is the direction of the x axis, and amplitudes of the wavefield wx+are maximum or above a predefined threshold in the direction of the x axis from the center location. In the other directions, amplitudes decrease smoothly as the angle between the other direction and the main direction of propagation increases, and reach a minimum amplitude or an amplitude below a predefined threshold in the direction opposite to the main direction of propagation. For example, in the snapshot104, amplitudes of the wavefield wx+reach a minimum amplitude or an amplitude below a predefined threshold in the opposite direction of the x axis from the center location.

Similar to decomposing the wavefield (i.e., the total pressure wavefield) into four characteristic wavefields (i.e., characteristic wavefields in Equations (5) to (8)), the upward and backward wavefields propagating in any given direction θ (i.e., wθ+and wθ−, respectively) can be constructed by decomposing the wavefield into the direction θ. For example, applying a rotation of angle θ to the eigenvectors of the matrices A and B, assuming that the angle θ is measured clockwise with respect to z axis, the rotated upward and backward wavefields can be obtained as:
wθ+=¼(p−ρ0c0(vzcos(θ)+vxsin(θ))  (9)
wθ−=¼(p+ρ0c0(vzcos(θ)+vxsin(θ))  (10)

FIG. 2illustrates example snapshots200of a rotated characteristic wavefield, according to some implementations. InFIG. 2, the wavefield102inFIG. 1is decomposed using, for example, Equation (9) into the direction θ.FIG. 2illustrates, at t=250 ms, snapshots of the decomposed rotated upward wavefield (i.e., wθ+) from 0 to 360 degrees with a 15-degree sampling.

In the present disclosure, the characteristic wavefield, for example, obtained from Equation (9) is used to compute an instantaneous source propagation angle for the source wavefield (e.g., upward) and an instantaneous receiver propagation angle for the receiver wavefield (e.g., downward), decompose the source wavefield by a direction of the instantaneous source propagation angle and decompose the receiver wavefield by a direction of the instantaneous receiver propagation angle, and finally compute angle domain common image gathers (ADCIG) using the decomposed source and receiver wavefields. Although this disclosure refers to the upward and downward decomposition, for purposes of example, the subject matter of this document can be applied to 3D propagation, leading to the introduction of a dip and azimuth angle or extended to more complex propagation physics, such as anisotropic or elastic propagation. For example, techniques described in Metivier et al., 2014 can be used for application of the above decomposition technique to “smart” absorbing layers for anisotropic and elastic propagations.

According to Equation (9) and illustrated inFIG. 2, wθ+provides, at any given time and any spatial position, amplitude of the wavefield propagating in the direction θ. Amplitudes of the wavefield wθ+are maximum or above a predefined threshold in the direction θ, compared to amplitudes in directions other than θ. Waves propagating in directions other than θ are attenuated with a maximum attenuation or an attenuation above a predefined threshold in the direction −θ. Therefore, the instantaneous wavefield propagation direction can be obtained by finding an angle where wθ+has a maximum absolute amplitude or an absolute amplitude above a predefined threshold. This angle corresponds to an angle where the derivative of wθ+is equal to zero. Applying trigonometric computations to Equation (9) gives:

∂wθ+∂θ
is equal to zero for the two following angles:

θ0=tan-1⁡(vxvz)⁢⁢and⁢⁢θ1=tan-1⁡(vxvz)+π(12)
As a result, the instantaneous direction of propagation of the wavefield is given by:

Computing the propagation angle using Equations (11) to (13) is comparable to the Poynting vector method (Yoon and Marfurt, 2006). However, the method using Equations (11) to (13) does not include any time derivatives of the pressure wavefield, and hence, there is no need, in some implementations, to access multiple time steps.

FIG. 3illustrates example results300of propagation angle estimation, according to some implementations. InFIG. 3, the wavefield310(same as the wavefield102inFIG. 1) is used for propagation angle estimation, according to Equation (13). The results320are the estimated propagation angle.

Wavefield Decomposition by Propagation Angle Direction:

Several methods can be used to decompose the wavefield by propagation angle direction. A first method uses the rotated characteristic wavefield wθ+in Equation (9) directly. However, as illustrated inFIGS. 1 and 2, the decomposition using the first method may not sufficiently attenuate the waves propagating in directions other than θ. As a result, the first decomposition method can be insufficient for some wavefield decomposition applications.

A second method to decompose the wavefield by propagation angle direction consists first of computing the propagation angle θpropaat each location, and then of assigning amplitude at this location to the propagation angle θpropa. However, the second method can be unstable due to the nature of tan−1. As a result, the second decomposition method can result in instability in complex or noisy situations.

A third method (i.e., the preferred method) consists of increasing the directional discrimination of wθ+(i.e., wθ+(x, z; t) in Equation (14) and below) by decomposing the amplitude over an interval of angle centered on the propagation angle given by Equation (13). A windowing function defines the angle interval. An example of such a windowing function is the Gaussian function gθ+(x, z; t):

σ can be defined, for example, using the following formula:

σ=1-cos⁢⁢Lθln⁢⁢q(15)
where Lθ=NwΔθ is smoothing window size, Nwis the number of points defining the window, Δθ is the angle resolution, and q is a parameter controlling the amplitude attenuation beyond Δθ In some implementations, q is usually set to 0.01. The pressure wavefield in the direction θ(Vθ+(x, z; t)) can then be obtained by forming the following quantity:
Vθ+(x,z;t)=gθ+(x,z;t)p(x,z;t)  (16)

With Equation (16), the sum of Vθ+(x, z; t) over the angles may not be equal to the total pressure wavefield (p(x, z; t)). Therefore, the function gθ+(x, z; t) needs to be normalized to eliminate the inequality problem. The normalization process is performed by first numerically computing:
E(x,z;t)=∫02πgθ+(x,z;t)dθ(17)
and then forming the decomposed wavefield in the direction θ as:

Wθ+⁡(x,z;t)=gθ+⁡(x,z;t)E⁡(x,z;t)⁢p⁡(x,z;t)(18)
If σ is very small, then the entire wavefield p(x, z; t) is associated with the angle θ. For other values of σ, the wavefield p(x, z; t) is spread over the angles near the angle θ. With the smoothing parameter σ in the angle domain, the stability and robustness of the proposed wavefield decomposition by propagation angle direction (e.g., Equation (18)) can be improved.

FIG. 4illustrates example wavefield and wavefields400for value of a smoothing parameter associated with various angle intervals, according to some implementations.402shows an example of a rotated characteristic wavefield wθ+.404shows Wθ+for value of σ associated with a window length equals to the angle resolution (1 degree).406shows Wθ+for value of σ associated with a window length equals to 5 angle resolution (5 degrees).408shows wθ+for value of σ associated with window length equals to 15 angle resolution (15 degrees).

Decomposition over an interval of angle can then be performed by summing the directional wavefields over the required interval.

FIG. 5illustrates example effects500of decreasing a smoothing parameter on a decomposed wavefield, according to some implementations. As illustrated inFIG. 5, from top left to bottom right, the effect of decreasing σ on the decomposed wavefield Wθ+using the two angles 0 and 180 degrees (i.e., up and down decomposition) in Equation (18) is shown.

FIG. 6illustrates example effects600of decreasing a smoothing parameter on a propagation angle in a decomposed wavefield, according to some implementations. As illustrated inFIG. 6, from top left to bottom right, the effect of decreasing σ on the propagation angle in the decomposed wavefield Wθ+using the two angles 0 and 180 degrees (i.e., up and down decomposition) in Equation (18) is shown.

FIG. 7illustrates example snapshots700of a pressure wavefield propagating in the Marmousi model, according to some implementations.FIG. 7illustrates a snapshot of the full wavefield702, a snapshot of the downward wavefield decomposition704using Equation (18) with an average amount of smoothing σ, a snapshot of the upward wavefield decomposition706using Equation (18) with the average amount of smoothing σ, a snapshot of the downward wavefield decomposition708using Equation (18) with a very small amount of smoothing σ, and a snapshot of the upward wavefield decomposition710using Equation (18) with the very small amount of smoothing σ. As illustrated in708and710, there exist instabilities and noises associated with small values of the smoothing parameter σ. In some implementations, the value of σ is set according to the sampling in angle of the function, for example, wθ+. Assuming, for example, that wθ+is sampled every 1 degree, a small σcorresponds to a smoothing window of 1 degree, an average σ corresponds to a smoothing window of 3 to 5 degrees, and a large σ corresponds to a smoothing window of 10 or more degrees.

In the present disclosure, the wavefield decomposition presented in the previous section (e.g., Equation (18)) is used to compute ADCIG. Without loss of generality, the conventional cross-correlation imaging condition is used for wave equation migration to form an image as:
I(x,z)=Σs=1NshotsΣt=1NtWsrc(s;x,z;t)Wrcv(s;x,z;t)  (19)
where I(x, z) represents the image at (x, z) acquired through migration, Nshotsis the number of shots considered for ADCIG computation, Ntis the number of time steps, and Wsrcand Wrcvare the source and receiver wavefields, respectively. Although this disclosure refers to the conventional cross-correlation imaging condition, for purposes of example, the subject matter of this document can use other imaging conditions for ADCIG computation without departing from the scope of the disclosure. The imaging condition used in Equation (19) can be extended to compute ADCIG, including scattering angle gathers (SAG) and dip angle gathers (DAG) by forming the quantities as:

The source and receiver angle wavefields (i.e., Wsrc(s; x, z; t; θsrc) and wrcv(s; x, z; t; θrcv) in Equations (20) and (21)) are computed, for example, by application of the directional wavefield decomposition (e.g., Equation (18)) to the source and receiver wavefields, to compute the SAG and DAG at the desired angles.

In general, standard reverse time migration (RTM) implementations are based on second order pressure wave equations. For acoustic RTM, stacked image computations require two forward modellings (source and receiver) and storing the pressure snapshots of the pressure wavefield (1 grid) at the time steps where the imaging condition is computed.

In 2D systems, using an acoustic isotropic time domain finite difference modelling and a cross-correlation imaging condition, the number of operations required by a standard RTM (Nrtm) is approximately:
Nrtm=2*Ngrid*Nop*Nt+2*Ngrid*Nsnp  (22)
with Ngridbeing the number of grid points, Nopbeing the number of kernel operations per grid point, Ntbeing the number of modelling time steps, and Nsnpbeing the number of snapshots. Assuming a second order wave system with 16thorder stencil in space, the number of kernel operation per grid point is Nop=53.

To apply the ADCIG computation method in the present disclosure (e.g., Equation (19)), two forward modellings are computed using a first-order (velocity-pressure) acoustic wave system. In addition to the pressure wavefield, the variables necessary to compute the characteristic wavefield in all directions are stored. For example, in 2D systems, the particle velocities vxand vzare stored, while vx, vy, and vzare stored in 3D systems. To compute the imaging condition, the ADCIG computation method in the present disclosure requires computing the propagation angle for the source and receiver wavefields, computing the windowed angle decomposition, and forming the angle gathers. The number of operation for the ADCIG computation method in the present disclosure (Nmeth) is approximately:
Nmeth=2*Ngrid*Nop*Nt+2*(Nw+15)*Ngrid*Nsnp  (23)
with Nwbeing the number of angles defining the windowing function gθ+. The number 15 corresponds to the approximate cost of the function arcTan used to obtain the propagation angle (e.g., in Equation (12)). When Nwis equal to one, the computation cost of the ADCIG computation method is the same as the Poynting vector method (e.g., without spatial smoothing). In practice, Nwis small (e.g., 5 to 10), and the computation cost of the ADCIG computation method with small Nwis low and comparable to the computation cost of the Poynting vector method.

FIG. 9illustrates example images900on a homogenous model, according to some implementations. InFIG. 9, a homogenous model is considered. The migrated data corresponds to one shot located at the center of the model. Receivers are located all around the model. Dips vary from 0 to 360 degrees with a 10-degree interval. The corresponding image (stacked image) is shown on the top ofFIG. 9. The three dip angle gathers (DAG) are shown below the stacked image inFIG. 9, and all three DAG give the same stacked image. Illustrated inFIG. 9from top to bottom, the three DAG use increasing smoothing parameter σ.

FIGS. 10-14illustrate example images on a dipping layer model, according to some implementations. InFIGS. 10-14, a homogenous background with a layer forming a 4.3-degree angle with the horizontal and three diffraction points located above the layer are considered.

FIG. 10illustrates example images1000on the dipping layer model using an exact velocity model (e.g., a correct velocity model). The stacked image is shown on the top ofFIG. 10. Five scattering angle gathers (SAG) at locations identified by arrows in the stacked image are shown below the stacked image inFIG. 10. Similarly,FIG. 11illustrates example images1100on the dipping layer model using a fast velocity model (e.g., a wrong velocity model). The stacked image is shown on the top ofFIG. 11. Five scattering angle gathers (SAG) at locations identified by arrows in the stacked image are shown below the stacked image inFIG. 11. As illustrated inFIGS. 10-11, the SAG with the exact velocity model (FIG. 10) show flat events, while the SAG with the fast velocity model (FIG. 11) show curved events.

FIGS. 12-13illustrate DAG, at the same locations of the SAG inFIGS. 10-11, using two different values of the smoothing parameter σ.FIG. 12illustrates example images1200using the exact velocity model with a value of the smoothing parameter σ corresponding to a 2 degrees window.FIG. 13illustrates example images1300using the exact velocity model with a value of the smoothing parameter σ corresponding to a 5 degrees window. As illustrated inFIGS. 12-13, the diffraction points are represented by flat lines, and the reflection points are represented by curved lines with a maximum amplitude at the dip location. In addition, using a larger smoothing parameter σ results in improved signal to noise ratio in the DAG domain.

FIG. 14illustrates example images1400with selective stack, according to some implementations. InFIG. 14, the stacked image (same as the stacked image inFIG. 10) is shown on the top. Selective stack can be performed to remove unwanted events in the stack image. For example, a selection in the SAG domain is performed to keep only the SAG between −40 degrees and +40 degrees. The final image after the selective stack is shown on the bottom ofFIG. 14.

FIG. 15is a flowchart of an example method1500for modeling angle domain common image gathers (ADCIG) from reverse time migration (RTM), according to some implementations. For clarity of presentation, the description that follows generally describes method1500in the context of the other figures in this description. For example, method1500can be performed by a computer system described inFIG. 16. However, it will be understood that method1500may be performed, for example, by any suitable system, environment, software, and hardware, or a combination of systems, environments, software, and hardware, as appropriate. In some implementations, various steps of method1500can be run in parallel, in combination, in loops, or in any order.

The method1500starts at block1505where seismic source and receiver wavefields are calculated based on seismic data. In some implementations, the seismic data is generated by sending seismic source waves, for example, from a seismic energy source such as an air gun, into earth subsurface layers and digitally sampling the seismic waves reflected by the earth subsurface layers. The seismic data captures reflected waves as a function of time. The seismic data can include amplitudes, phases, or both, of the reflected waves. In some implementations, the seismic source and receiver wavefields are calculated using a time-domain wave equation solver. In some implementations, the seismic source and receiver wavefields are calculated using a frequency-domain wave equation solver.

At block1510, characteristic source and receiver wavefields are calculated from the seismic source and receiver wavefields (e.g., Equation (9)). In some implementations, the characteristic source and receiver wavefields are calculated for a plurality of predefined time steps with a first-order acoustic wave equation propagating model (e.g., Equation (1)). In some implementations, the characteristic source and receiver wavefields are calculated for a plurality of predefined time steps with a second-order acoustic wave equation propagating model.

At block1515, propagation angles are calculated for both the characteristic source and receiver wavefields (e.g., Equation (13)). In some implementations, the characteristic source wavefield has an absolute amplitude, at the propagation angle for the characteristic source wavefield, larger than an absolute amplitude at any other angles different from the propagation angle for the characteristic source wavefield. The characteristic receiver wavefield has an absolute amplitude, at the propagation angle for the characteristic receiver wavefield, larger than an absolute amplitude at any other angles different from the propagation angle for the characteristic receiver wavefield. In some implementations, the propagation angles are calculated by applying trigonometric computations on the characteristic source and receiver wavefields (e.g., Equations (11) to (13)).

At block1520, a wavefield decomposition algorithm is applied on the characteristic source and receiver wavefields to obtain corresponding directional source and receiver wavefields. The wavefield decomposition algorithm decomposes wavefield amplitude of a wavefield in an angle interval centered on a propagation angle of the wavefield (e.g., Equations (14) to (18)). In some implementations, the angle interval is defined by a windowing function, such as a Gaussian function. In some implementations, the angle interval can be 1 degree, 5 degrees, or 15 degrees. In some other implementations, the angle interval can be any other number of degrees.

At block1525, angle domain common image gathers (ADCIG) are formed by applying an image condition to the obtained directional source and receiver wavefields. In some implementations, the ADCIG are formed using a reverse time migration (RTM) process. In some implementations, the image condition is a conventional cross-correlation imaging condition (e.g., used in Equation (19)). In some implementations, the formed ADCIG are scattering angle gathers (SAG) (e.g., Equation (20)). In some implementations, the formed ADCIG are dip angle gathers (DAG) (e.g., Equation (21)).

The example method1500shown inFIG. 15can be modified or reconfigured to include additional, fewer, or different steps (not shown inFIG. 15), which can be performed in the order shown or in a different order. In some implementations, one or more of the steps shown inFIG. 15can be repeated or iterated, for example, until a terminating condition is reached. In some implementations, one or more of the individual steps shown inFIG. 15can be executed as multiple separate steps, or one or more subsets of the steps shown inFIG. 15can be combined and executed as a single step. In some implementations, one or more of the individual steps shown inFIG. 15may also be omitted from the example method1500.

The computer1602can serve in a role as a client, network component, a server, a database or other persistency, or any other component (or a combination of roles) of a computer system for performing the subject matter described in the instant disclosure. The illustrated computer1602is communicably coupled with a network1630. In some implementations, one or more components of the computer1602may be configured to operate within environments, including cloud-computing-based, local, global, or other environment (or a combination of environments).

The computer1602can receive requests over network1630from a client application (for example, executing on another computer) and responding to the received requests by processing the received requests using the appropriate software application(s). In addition, requests may also be sent to the computer1602from internal users (for example, from a command console or by other appropriate access method), external or third-parties, other automated applications, as well as any other appropriate entities, individuals, systems, or computers.

Each of the components of the computer1602can communicate using a system bus1603. In some implementations, any or all of the components of the computer1602, both hardware or software (or a combination of hardware and software), may interface with each other or the interface1604(or a combination of both) over the system bus1603using an application programming interface (API)1612or a service layer1613(or a combination of the API1612and service layer1613). The API1612may include specifications for routines, data structures, and object classes. The API1612may be either computer-language independent or dependent and refer to a complete interface, a single function, or even a set of APIs. The service layer1613provides software services to the computer1602or other components (whether or not illustrated) that are communicably coupled to the computer1602. The functionality of the computer1602may be accessible for all service consumers using this service layer. Software services, such as those provided by the service layer1613, provide reusable, defined functionalities through a defined interface. For example, the interface may be software written in JAVA, C++, or other suitable language providing data in extensible markup language (XML) format or other suitable format. While illustrated as an integrated component of the computer1602, alternative implementations may illustrate the API1612or the service layer1613as stand-alone components in relation to other components of the computer1602or other components (whether or not illustrated) that are communicably coupled to the computer1602. Moreover, any or all parts of the API1612or the service layer1613may be implemented as child or sub-modules of another software module, enterprise application, or hardware module without departing from the scope of this disclosure.

The computer1602includes an interface1604. Although illustrated as a single interface1604inFIG. 16, two or more interfaces1604may be used according to particular needs, desires, or particular implementations of the computer1602. The interface1604is used by the computer1602for communicating with other systems that are connected to the network1630(whether illustrated or not) in a distributed environment. Generally, the interface1604comprises logic encoded in software or hardware (or a combination of software and hardware) and is operable to communicate with the network1630. More specifically, the interface1604may comprise software supporting one or more communication protocols associated with communications such that the network1630or interface's hardware is operable to communicate physical signals within and outside of the illustrated computer1602.

The computer1602includes a processor1605. Although illustrated as a single processor1605inFIG. 16, two or more processors may be used according to particular needs, desires, or particular implementations of the computer1602. Generally, the processor1605executes instructions and manipulates data to perform the operations of the computer1602and any algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure.

The computer1602also includes a database1606that can hold data for the computer1602or other components (or a combination of both) that can be connected to the network1630(whether illustrated or not). For example, database1606can be an in-memory, conventional, or other type of database storing data consistent with this disclosure. In some implementations, database1606can be a combination of two or more different database types (for example, a hybrid in-memory and conventional database) according to particular needs, desires, or particular implementations of the computer1602and the described functionality. Although illustrated as a single database1606inFIG. 16, two or more databases (of the same or combination of types) can be used according to particular needs, desires, or particular implementations of the computer1602and the described functionality. While database1606is illustrated as an integral component of the computer1602, in alternative implementations, database1606can be external to the computer1602. As illustrated, the database1606holds seismic data1616, source and receiver wavefields1618, characteristic source and receiver wavefields1620, decomposed source and receiver wavefields1622, and the angle domain common image gathers (ADCIG)1624.

The computer1602also includes a memory1607that can hold data for the computer1602or other components (or a combination of both) that can be connected to the network1630(whether illustrated or not). For example, memory1607can be random access memory (RAM), read-only memory (ROM), optical, magnetic, and the like storing data consistent with this disclosure. In some implementations, memory1607can be a combination of two or more different types of memory (for example, a combination of RAM and magnetic storage) according to particular needs, desires, or particular implementations of the computer1602and the described functionality. Although illustrated as a single memory1607inFIG. 16, two or more memories1607(of the same or combination of types) can be used according to particular needs, desires, or particular implementations of the computer1602and the described functionality. While memory1607is illustrated as an integral component of the computer1602, in alternative implementations, memory1607can be external to the computer1602.

The application1608is an algorithmic software engine providing functionality according to particular needs, desires, or particular implementations of the computer1602, particularly with respect to functionality described in this disclosure. For example, application1608can serve as one or more components, modules, or applications. Further, although illustrated as a single application1608, the application1608may be implemented as multiple applications1608on the computer1602. In addition, although illustrated as integral to the computer1602, in alternative implementations, the application1608can be external to the computer1602.

There may be any number of computers1602associated with, or external to, a computer system containing computer1602, each computer1602communicating over network1630. Further, the term “client,” “user,” and other appropriate terminology may be used interchangeably as appropriate without departing from the scope of this disclosure. Moreover, this disclosure contemplates that many users may use one computer1602, or that one user may use multiple computers1602.

Accordingly, the previously-described example implementations do not define or constrain this disclosure. Other changes, substitutions, and alterations are also possible without departing from the spirit and scope of this disclosure.

REFERENCES