Efficient internal multiple prediction methods

Abstract

Methods for processing seismic data are described. The method includes: obtaining seismic data; solving a series of partial differential wave equations, wherein a first partial differential wave equation describes propagation of a seismic wave going from a first reflector to a second reflector, wherein a second partial differential wave equation describes propagation of a seismic wave going from a second reflector to a third reflector, and wherein a third partial differential wave equation describes propagation of a seismic wave going from a third reflector to a seismic receiver, wherein outputting predicted internal multiples for further imaging or attenuation.

Claims

1. A method for seismic data processing comprising: obtaining seismic data for a subterranean structure, the seismic data including contribution from internal multiples; solving a series of partial differential wave equations each representing a reflected portion of predicted internal multiples of the seismic data, wherein a first partial differential wave equation describes propagation of a seismic wave going from a first reflector to a second reflector, wherein a second partial differential wave equation describes propagation of the seismic wave going from the second reflector to a third reflector, and wherein a third partial differential wave equation describes propagation of the seismic wave going from the third reflector to a seismic receiver, the third partial differential wave equation providing the predicted internal multiples; outputting the predicted internal multiples for attenuation; and removing the predicted internal multiples from an image of the subterranean structure.

2. The method of claim 1, wherein the series of partial differential wave equations are solved sequentially.

3. The method of claim 1, wherein at least one of the series of partial differential wave equations describes 2-D or 3-D propagation of the seismic wave.

4. The method of claim 1, wherein the first partial differential equation is $\left(\frac{}{z}+\mathrm{iq}\right){\hat{p}}_1\left(k;;z;k_s\right)={\hat{b}}_1\left(k;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_0\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}},$ wherein z is depth between z.sub.1 and z.sub.2, {circumflex over (p)}.sub.0 is Fourier transform of p.sub.0 which is Fourier conjugate to field point going from source to the reflections above depth z.sub.1, {circumflex over (p)}.sub.1 is Fourier transform of p.sub.1 which is Fourier conjugate to field point going from z.sub.1 to z.sub.2, {circumflex over (b)}.sub.1 is Fourier transform of b.sub.1 (b.sub.1 is first term of inverse scattering series), {right arrow over (k)} is wave vector directed towards observation point, and k.sub.s is Fourier conjugate of x.sub.s, is radial frequency, k.sub.s is Fourier conjugate of x.sub.s which is source coordinates, and $q=\frac{}{c}\sqrt{1-\frac{c^2{k}_x^2}{{}^2}}$ wherein c is velocity and k.sub.x is wavenumber.

5. The method of claim 1, wherein the second partial differential equation is $\left(\frac{}{z}-\mathrm{iq}\right){\hat{p}}_2\left(k;;z;k_s\right)=-{\hat{b}}_1\left(k;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_1\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}},$ wherein {circumflex over (p)}.sub.2 is Fourier transform of p.sub.2 which is Fourier conjugate to field point going from z.sub.2 to z.sub.3.

6. The method of claim 1, wherein the third partial differential equation is $\left(\frac{}{z}+{\mathrm{iq}}_g\right){\hat{p}}_3\left(k_g;;z;k_s\right)={\hat{b}}_1\left(k_g;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_2\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}},$ wherein {circumflex over (p)}.sub.2 is Fourier transform of p.sub.3 which is Fourier conjugate to field point going from z.sub.3 to the seismic receiver, q.sub.g is depth wavenumbers, and k.sub.g is Fourier conjugate of x.sub.g which is receiver coordinates.

7. The method of claim 1, wherein the seismic data is 2-D or 3-D.

8. A method for seismic data processing comprising: obtaining seismic data for a subterranean structure, the seismic data including contribution from internal multiples; solving a series of partial differential wave equations each representing a reflected portion of predicted internal multiples of the seismic data, wherein a first partial differential wave equation describes propagation of a seismic wave going from a first reflector to a second reflector, wherein a second partial differential wave equation describes propagation of the seismic wave going from the second reflector to a third reflector, and wherein a third partial differential wave equation describes propagation of the seismic wave going from the third reflector to a seismic receiver, the third partial differential wave equation providing the predicted internal multiples; and outputting the predicted internal multiples.

9. The method of claim 8, wherein the series of partial differential wave equations are solved sequentially.

10. The method of claim 8, wherein at least one of the series of partial differential wave equations describes 2-D or 3-D propagation of the seismic wave.

11. The method of claim 8, wherein the first partial differential equation is $\left(\frac{}{z}+\mathrm{iq}\right){\hat{p}}_1\left(k;;z;k_s\right)={\hat{b}}_1\left(k;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_0\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}},$ wherein z is depth between z.sub.1 and z.sub.2, {circumflex over (p)}.sub.0 is Fourier transform of p.sub.0 which is Fourier conjugate to field point going from source to reflection at depth z.sub.1, {circumflex over (p)}.sub.1 is Fourier transform of p.sub.1 which is Fourier conjugate to field point going from z.sub.1 to z.sub.2, {circumflex over (b)}.sub.1 is Fourier transform of b.sub.1 (b.sub.1 is first term of inverse scattering series), {right arrow over (k)} is wave vector directed towards observation point, and k.sub.s is Fourier conjugate of x.sub.s, is radial frequency, k.sub.s is Fourier conjugate of x.sub.s which is source coordinates, and $q=\frac{}{c}\sqrt{1-\frac{c^2{k}_x^2}{{}^2}}$ wherein c is velocity and k.sub.x is wavenumber.

12. The method of claim 8, wherein the second partial differential equation is $\left(\frac{}{z}-\mathrm{iq}\right){\hat{p}}_2\left(k;;z;k_s\right)=-{\hat{b}}_1\left(k;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_1\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}},$ wherein {circumflex over (p)}.sub.2 is Fourier transform of p.sub.2 which is Fourier conjugate to field point going from z.sub.2 to z.sub.3.

13. The method of claim 8, wherein the third partial differential equation is $\left(\frac{}{z}+{\mathrm{iq}}_g\right){\hat{p}}_3\left(k_g;;z;k_s\right)={\hat{b}}_1\left(k_g;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_2\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}},$ wherein {circumflex over (p)}.sub.3 is Fourier transform of p.sub.3 which is Fourier conjugate to field point going from z.sub.3 to the seismic receiver, q.sub.g is depth wavenumbers, and k.sub.g is Fourier conjugate of x.sub.g which is receiver coordinates.

14. The method of claim 8, wherein the seismic data is 2-D or 3-D.

15. A system for seismic data processing comprising: one or more seismic sources generating at least one seismic wave propagating into a subterranean structure; a plurality of reflectors reflecting the at least one seismic wave, the plurality of reflectors including a first reflector, a second reflector, and a third reflector; one or more seismic receivers capturing seismic data for the subterranean structure following reflection of the at least one seismic wave from the plurality of reflectors, the seismic data including contribution from internal multiples; and a computer system generating predicted internal multiples by solving a series of partial differential wave equations each representing a reflected portion of the predicted internal multiples of the seismic data, the series of partial differential wave equations including a first partial differential wave equation, a second partial differential wave equation, and a third partial differential wave equation, the first partial differential wave equation describing propagation of the at least one seismic wave from the first reflector to the second reflector, the second partial differential wave equation describing propagation of the at least one seismic wave from the second reflector to the third reflector, the third partial differential wave equation describing propagation of the at least one seismic wave from the third reflector to the one or more seismic receivers.

16. The system of claim 15, wherein the series of partial differential wave equations are solved sequentially.

17. The system of claim 15, wherein at least one of the series of partial differential wave equations describes 2-D or 3-D propagation of the at least one seismic wave.

18. The system of claim 15, wherein the seismic data is 2-D or 3-D.

19. The system of claim 15, wherein the computer system removes the predicted internal multiples from a seismic image of the subterranean structure.

20. The system of claim 15, wherein the computing system outputs the predicted internal multiples.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 illustrates a schematic representation of seismic wave propagating in the sub surface.

(2) FIG. 2A-2C illustrates an embodiment of the invention as described in the Example.

(3) FIG. 3 shows a flow diagram that summarizes an embodiment of the invention.

DETAILED DESCRIPTION

(4) Turning now to the detailed description of the preferred arrangement or arrangements of the present invention, it should be understood that the inventive features and concepts may be manifested in other arrangements and that the scope of the invention is not limited to the embodiments described or illustrated.

(5) The following examples of certain embodiments of the invention are given. Each example is provided by way of explanation of the invention, one of many embodiments of the invention, and the following examples should not be read to limit the scope of the invention.

(6) Geoscientists often process acquired seismic data to improve interpretability of the seismic data. Seismic processing includes steps or operations that are intended to remove or attenuate undesirable noise, enhance signal, and the like. Specific examples of seismic processing steps include random noise (e.g., noise caused by weather) attenuation, free surface multiple prediction and attenuation (for offshore data), internal multiple prediction and attenuation, velocity modeling, and imaging/migration.

(7) The present invention provides method and system for efficiently predicting and attenuating internal multiples by utilizing cascaded double-square-root one-way wave equations and proposing a partial differential equation (PDE) based internal multiple prediction method. This new method is mathematically equivalent to conventional inverse scattering series approach. Both the PDE-based approach and conventional ISS method are data-driven approaches that do not require a priori knowledge (e.g., location of reflectors, earth model type) of the subsurface. Computationally, the PDE-based approach is much more time efficient without sacrificing power or accuracy. For 2D application, the present invention can speed up internal multiple prediction by a factor of about 50 to 100. For 3D applications, the speed gain is even more dramatic, potentially improving computational efficiency by two orders of magnitudes or greater. Other advantages are apparent from the disclosure herein. 3D refers to three-dimensional, or any computer-generated image in which the displayed items have a length, width, and depth.

(8) Internal Multiple Attenuation Approach Using Inverse Scattering Series

(9) Inverse scattering series is largely based on scattering theory, which is a form of perturbation analysis in which perturbation in properties of a medium is related to perturbation in a wavefield that experiences the medium. The difference between actual and reference media can be characterized by a perturbation operator. The corresponding difference between actual and reference wavefields is referred to as the scattered wavefield. Forward scattering takes the reference medium, reference wavefield, and perturbation operator as inputs and outputs the actual wavefield. Inverse scattering takes the reference medium, reference wavefield, values of the actual field on the measurement surface as inputs and outputs the difference between actual and reference medium properties through the perturbation operator.

(10) The collected seismic data may be processed (e.g., deghosting, free-surface multiple elimination, etc.) prior to undergoing internal multiple attenuation. The inverse scattering series constructs the perturbation operator through the action of applying Green's function operators to the seismic data. A subseries (removal series) of the inverse scattering series can be linearly combined with the seismic data in order to produce seismic data that is largely free of internal multiples. The inverse scattering series for attenuating all internal multiples chooses the leading and most significant contribution from the removal series of each order of multiple (thus forming a series that attenuates rather than fully eliminates). This subseries can be linearly combined with the seismic data to produce seismic data that is largely free of internal multiples. A detailed description of ISS for multiples attenuation can be found in U.S. Pat. No. 5,757,723 and (Weglein et al. 1997 Geophysics), relevant contents of which are hereby incorporated by reference.

(11) Assuming 2-D wave propagation, equation (1) defines a first-order internal multiple attenuation term b.sub.3 (portion of the third term of the internal multiple attenuation series) in f-k domain:

(12) $\begin{array}{cc}\begin{array}{c}{b}_3=\frac{1}{{\left(2\right)}^2}{}_{-}^{}{}_{-}^{}{\mathrm{dk}}_1{e}^{-{\mathrm{iq}}_1\left({}_g-{}_s\right)}{\mathrm{dk}}_2{e}^{{\mathrm{iq}}_2\left({}_g-{}_s\right)}\\ {}_{-}^{}{\mathrm{dz}}_1{e}^{i\left(q_g+q_1\right)z_1}{b}_1\left(k_g,-k_1,z_1\right)\\ {}_{-}^{z_1}{\mathrm{dz}}_2{e}^{i\left(-q_1-q_2\right)z_2}{b}_1\left(k_1,-k_2,z_2\right)\\ {}_{z_2}^{}{\mathrm{dz}}_3{e}^{i\left(q_2+q_s\right)z_3}{b}_1\left(k_2,-k_s,z_3\right);z_1>z_2\mathrm{and}{z}_2<z_3\end{array}& \left(1\right)\end{array}$
where k.sub.g and k.sub.s are Fourier conjugates of receiver and source coordinates respectively, is the temporal frequency, q.sub.g and g.sub.s are depth wavenumbers, z.sub.1, z.sub.2, and z.sub.3 are the depth of three reflectors which generate the predicted internal multiples (FIG. 1), b.sub.1 is the reflectivity model generated by prestack seismic data using Stolt migration, b.sub.3 is the predicted internal multiples, k.sub.1 and k.sub.2 are the wavenumbers from either source side or the receiver side q.sub.1, q.sub.2 and q.sub.3 are the single-square-root propagation operators defined in frequency-wavenumber domain (see equation 3 below) .sub.g.sub.s are small depth perturbations. Higher-order terms corresponds to higher-order internal multiples.

(13) FIG. 1 illustrates a first-order internal multiple configuration. The seismic source is located at 1 and the seismic receiver is located 2. Both the source and the receiver are located at the free surface.

(14) PDE Based Internal Multiple Attenuation

(15) The present invention provides a PDE-based internal multiple prediction and attenuation approach. A 2-D acoustic wave equation for upcoming waves in frequency-space domain (neglecting spatial derivatives of velocity) can be written as equation (2):

(16) $\begin{array}{cc}\left(\frac{}{z}-\mathrm{iq}\right)p_0\left(k;;z;x_s\right)=-\left(z\right)e^{{\mathrm{ikx}}_s}& \left(2\right)\end{array}$
where p.sub.0(k; ; z; x.sub.s) is Fourier conjugate to field point going from source to reflection at depth z.sub.1. is radial frequency, x.sub.s is source coordinates, and q is given by equation (3):

(17) $\begin{array}{cc}q=\frac{}{c}\sqrt{1-\frac{c^2{k}_x^2}{{}^2}}& \left(3\right)\end{array}$
where c is velocity and k.sub.x is wavenumber, or Fourier conjugate of lateral propagation locations x.

(18) Equation (4) describes 2-D propagation of an acoustic wave from reflector located at z.sub.1 to reflector located at z.sub.2:

(19) $\begin{array}{cc}\left(\frac{}{z}+\mathrm{iq}\right){\hat{p}}_1\left(k;;z;k_s\right)={\hat{b}}_1\left(k;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_0\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}}& \left(4\right)\end{array}$
where {circumflex over (p)}.sub.0 is Fourier transform of p.sub.0, {circumflex over (p)}.sub.1 is Fourier transform of p.sub.1 which is Fourier conjugate to field point going from z.sub.1 to z.sub.2, {circumflex over (b)}.sub.1 is Fourier transform of b.sub.1 (b.sub.1 is first term of inverse scattering series, constructed by seismic data), {right arrow over (k)} is wave vector directed towards observation point, and k.sub.s is Fourier conjugate of x.sub.s.

(20) Equation (5) describes propagation of the acoustic wave from reflector located at z.sub.2 to reflector at z.sub.3:

(21) $\begin{array}{cc}\left(\frac{}{z}-\mathrm{iq}\right){\hat{p}}_2\left(k;;z;k_s\right)=-{\hat{b}}_1\left(k;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_1\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}}& \left(5\right)\end{array}$
where {circumflex over (p)}.sub.2 is Fourier transform of p.sub.2 which is Fourier conjugate to field point going from z.sub.2 to z.sub.3.

(22) Partial differential equation (6) describes propagation of the acoustic wave from reflector located at z.sub.3 to receiver 2:

(23) $\begin{array}{cc}\left(\frac{}{z}+{\mathrm{iq}}_g\right){\hat{p}}_3\left(k_g;;z;k_s\right)={\hat{b}}_1\left(k_g;-\stackrel{.fwdarw.}{k^{}};z\right){\hat{p}}_2\left(\stackrel{.fwdarw.}{k^{}};;z;k_s\right)d\stackrel{.fwdarw.}{k^{}}& \left(6\right)\end{array}$
where {circumflex over (p)}.sub.3 is Fourier transform of p.sub.3 which is Fourier conjugate to field point going from z.sub.3 to receiver 2, q.sub.g is depth wavenumbers.

(24) Solving for {circumflex over (p)}.sub.3 where z=0 results in equation (7),
{circumflex over (p)}.sub.3(k.sub.g;;z=0;k.sub.s)=.sub.0.sup..sub..sup.e.sup.i(q.sup.g.sup.+q.sup.2.sup.)z.sup.3.sub.0.sup.z.sup.3.sub..sup.e.sup.i(q.sup.2.sup.+q.sup.1.sup.)z.sup.2.sub.z.sub.2.sup.e.sup.i(q.sup.1.sup.+q.sup.s.sup.)z.sup.1{circumflex over (b)}.sub.1(k.sub.1;k.sub.s;z.sub.1)dz.sub.1{circumflex over (b)}.sub.1(k.sub.2;k.sub.1;z.sub.2)dk.sub.1dz.sub.2{circumflex over (b)}.sub.1(k.sub.g;k.sub.2;z.sub.3)dk.sub.2dz.sub.3(7)
Equation (7) is equivalent to b.sub.3, which is defined by equation (1). The new approach solves three wave equations (4, 5 and 6) sequentially. Solving each equation requires 4 computational loops over , z and two k.sub.s. Altogether, only at most 4 integration loops (p1, p2, p3, ) is required. Solving equation (7) at z=0 provides predicted internal multiples. The predicted internal multiples may be outputted on any suitable medium (e.g., computer screen, printout, mobile devices, etc.)

(25) FIG. 3 illustrates an embodiment of the present invention in a flow diagram. First step 301 involves collecting or obtaining seismic data. Typically, this seismic data will be a composite that includes contribution from internal multiples. This seismic data can be used to construct b.sub.1, the first term of inverse scattering series. In the next step 302, a series of PDE wave equations can be constructed. Each PDE wave equation describes a reflected portion of the predicted internal multiples. In step 303, the series of PDE wave equations are solved. Each PDE wave equations solves a term that can be utilized to solve later PDE wave equation(s). The third PDE wave equation solved at z=0 provides predicted internal multiples. In the last step 304, the predicted internal multiples are outputted. The predicted internal multiples may be outputted on any suitable medium (e.g., computer screen, printout, mobile devices, etc.)

EXAMPLE

(26) Internal multiple prediction was performed on a 2-D synthetic data using method of the present invention and conventional inverse scattering series. The synthetic dataset contained 641 shots with 12.5 m shot interval and 12.5 m receiver interval. Sampling rate was 4 ms and trace length was 2 s. FIG. 2A illustrates a single shot record as the input seismic data for the multiple prediction, which includes both primary and the internal multiples. FIG. 2B shows the predicted internal multiples of the same shot record using PDE-based approach. FIG. 2C shows the predicted internal multiples of the same shot record using the conventional inverse scattering series approach. The two methods provide almost identical internal multiple predictions. However, the PDE-based approach took 17 minutes while conventional inverse scattering series approach took 28 hours. The speedup ratio of the present invention on this dataset is about 99 times. This ratio is expected to be more dramatic for 3D datasets. The use of computers is essential to process the seismic data quickly and effectively.

(27) The foregoing patent application is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. All such modifications and variations are intended to be within the scope of the present invention, as defined in the appended claims. Persons skilled in the art will readily recognize that in practical applications of the invention, at least some of the steps in the present inventive method are performed on or with the aid of a computer, i.e. the invention is computer implemented.

(28) A suitable computer or computer system may be in communication with disk storage devices such as an external hard disk storage devices or conventional hard disk drives. These storage drives/devices can be implemented by way of a local area network or by remote access. The storage device may be used to store any and all of the program instructions, measurement data, and results as desired.

(29) In one implementation, data may be stored in disk storage device. The system computer may retrieve the appropriate data from the disk storage devices to process data according to program instructions that correspond to implementations of various techniques described herein. The program instructions may be written in a computer programming language, such as C++, Java and the like. The program instructions may be stored in a computer-readable medium, such as program disk storage device. Such computer-readable media may include computer storage media. Computer storage media may include volatile and non-volatile, and removable and non-removable media implemented in any method or technology for storage of information, such as computer-readable instructions, data structures, program modules or other data. Computer storage media may further include RAM, ROM, erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other solid state memory technology, CD-ROM, digital versatile disks (DVD), or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by the system computer. Combinations of any of the above may also be included within the scope of computer readable media. In one implementation, the system computer may present output primarily onto graphics display, or alternatively via printer.