RTM using equal area spherical binning for generating image angle gathers

11487033 · 2022-11-01

Assignee

Inventors

Cpc classification

International classification

Abstract

Seismic exploration of an underground formation uses seismic excitations to probe the formation's properties such as reflectivity that can be imaged using reverse time migration. Using an equal area spherical binning at reflection points improves and simplifies RTM imaging together with adaptability to the data acquisition geometry, while overcoming drawbacks of conventional cylindrical binning.

Claims

1. A method for seismic exploration of an underground formation, the method comprising: obtaining seismic data representing seismic excitations emitted from a source position and detected at a receiver position after traveling through the underground formation; and processing the seismic data using a reverse time migration, RTM, method to generate image angle gathers, wherein with an equal area spherical binning technique is employed to define incident and azimuth angle bins at reflection points, and the image angle gathers are used to generate the underground formation image enabling to determine the presence of hydrocarbons and/or other minerals.

2. The method of claim 1, wherein the equal area spherical binning technique is determined by three numbers: two numbers characterize a base angular configuration, and a third number indicates a number of divisions on each side of a bin of the base angular configuration.

3. The method of claim 2, further comprising: selecting a combination of the three numbers using a subset of the seismic data, the three numbers being varied within predetermined value ranges to identify the combination yielding a best achievable angular binning for a seismic data acquisition system used to acquire the seismic data.

4. The method of claim 1, wherein the processing of the seismic data using RTM with the equal area spherical binning comprises: initializing partial images corresponding to angular bins; and for each collection of seismic traces in the seismic data from the same source position or the same receiver position: computing a forward-propagated wavefield from the source position through the underground formation and a backward-propagated wavefield from the receiver position through the underground formation, applying an imaging condition in a given domain, using the forward-propagated wavefield and the backward-propagated wavefield for each point in the given domain computing incident and azimuth reflection angles between forward and backward wavefields, identifying one of the angular bins defined by the equal area spherical binning technique with a bin center closest to the incident and azimuth reflection angles; and adding a result of the imaging condition in the given domain point to one of the partial images corresponding to the one of the angular bins.

5. The method of claim 4, wherein the incidence and azimuth reflection angles are calculated in a wavenumber domain, the forward and backward-propagated wavefields being transformed in the wavenumber domain using an anti-leakage Fourier transform, and angle domain image gathers being converted back in a space domain by applying a fast Fourier transform.

6. The method of claim 4, wherein the incidence and azimuth reflection angles are calculated using derivatives of the wavefields.

7. The method of claim 1, further comprising: interpolating the image angle gathers to generate another set of image angle gathers at positions different from bin centers of bins as defined by the equal area spherical binning technique.

8. The method of claim 1, further comprising: using the image angle gathers to update a model used to compute the forward-propagated wavefield and the backward-propagated wavefield.

9. A seismic data processing apparatus comprising: an interface configured to obtain seismic data representing seismic excitations emitted from a source position and detected at a receiver position after traveling through the underground formation; and a seismic data processing unit connected to the interface and configured to process the seismic data using a reverse time migration, RTM, method to generate image angle gathers, wherein with an equal area spherical binning technique is employed to define incident and azimuth angle bins, and the image angle gathers are used generate the underground formation images enabling to determine the presence of hydrocarbons and/or other minerals.

10. The seismic data processing apparatus of claim 9, wherein the equal area spherical binning technique is determined by three numbers: two numbers characterize a base angular configuration, and a third number indicates a number of divisions on each side of a bin of the base angular configuration.

11. The seismic data processing apparatus of claim 10, wherein the seismic data processing unit is configured to select a combination of the three numbers using a subset of the seismic data, the three numbers being varied within predetermined value ranges to identify the combination yielding a best achievable angular binning for a seismic data acquisition system used to acquire the seismic data.

12. The seismic data processing apparatus of claim 9, wherein, when processing of the seismic data using RTM with the equal area spherical binning, the seismic data processing unit performs: initializing partial images corresponding to angular bins; and for each collection of seismic traces in the seismic data from the same source position or the same receiver position: computing a forward-propagated wavefield from the source position through the underground formation and a backward-propagated wavefield from the receiver position through the underground formation, applying an imaging condition in a given domain, using the forward-propagated wavefield and the backward-propagated wavefield for each point in the given domain computing incident and azimuth reflection angles between forward and backward wavefields, identifying one of the angular bins defined by the equal area spherical binning technique with a bin center closest to the incident and azimuth reflection angles; and adding a result of the imaging condition in the given domain point to one of the partial images corresponding to the one of the angular bins.

13. The seismic data processing apparatus of claim 12, wherein the seismic data processing unit calculates the incidence and azimuth reflection angles in a wavenumber domain, transforms the forward and backward-propagated wavefields in the wavenumber domain using an anti-leakage Fourier transform, and angle domain image gathers are converted back in a space domain by applying a fast Fourier transform.

14. The seismic data processing apparatus of claim 12, wherein the seismic data processing unit calculates the incidence and azimuth reflection angles using derivatives of the wavefields.

15. The seismic data processing apparatus of claim 9, wherein the seismic data processing unit is further configured to interpolate the image angle gathers to generate another set of image angle gathers at positions different from bin centers of bins as defined by the equal area spherical binning technique.

16. The seismic data processing apparatus of claim 9, wherein the seismic data processing unit is further configured to use the image angle gathers to update a model used to compute the forward-propagated wavefield and the backward-propagated wavefield.

17. A computer readable recording medium storing executable codes, which, when executed on a computer make the computer perform a method for seismic exploration of an underground formation, the method comprising: obtaining seismic data representing seismic excitations emitted from a source position and detected at a receiver position after traveling through the underground formation; and processing the seismic data using a reverse time migration, RTM, method to generate image angle gathers, wherein with an equal area spherical binning technique is employed to define incident and azimuth angle bins at reflection points, and the image angle gathers are used generate the underground formation image enabling to determine the presence of hydrocarbons and/or other minerals.

18. The computer readable recording medium of claim 17, wherein the equal area spherical binning technique is determined by three numbers: two numbers characterize a base angular configuration, and a third number indicates a number of divisions on each side of a bin of the base angular configuration, the method computer readable recording medium further comprising selecting a combination of the three numbers using a subset of the seismic data, the three numbers being varied within predetermined value ranges to identify the combination yielding a best achievable angular binning for a seismic data acquisition system used to acquire the seismic data.

19. The computer readable recording medium of claim 17, wherein the processing of the seismic data using RTM with the equal area spherical binning comprises: initializing partial images corresponding to angular bins; and for each collection of seismic traces in the seismic data from the same source position or the same receiver position: computing a forward-propagated wavefield from the source position through the underground formation and a backward-propagated wavefield from the receiver position through the underground formation, applying an imaging condition in a given domain, using the forward-propagated wavefield and the backward-propagated wavefield for each point in the given domain computing incident and azimuth reflection angles between forward and backward wavefields, identifying one of the angular bins defined by the equal area spherical binning technique with a bin center closest to the incident and azimuth reflection angles; and adding a result of the imaging condition in the given domain point to one of the partial images corresponding to the one of the angular bins.

20. The computer readable recording medium of claim 17, wherein the method further comprises interpolating the image angle gathers to generate another set of image angle gathers at positions different from bin centers of bins as defined by the equal area spherical binning technique.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate one or more embodiments and, together with the description, explain these embodiments. In the drawings:

(2) FIG. 1 illustrates ocean bottom nodes used in acquiring seismic data;

(3) FIG. 2 illustrates incident and azimuth angle at a reflection point;

(4) FIG. 3 illustrates incident angle ranges' dependence on interface depths;

(5) FIG. 4 shows two-dimensional projections of incident azimuth angle pairs obtained using (a) a conventional cylindrical coordinates approach versus (b) a source-centric binning on an icosahedron;

(6) FIG. 5 is flowchart of a method using equal area spherical binning for RTM angle gathers according to an embodiment;

(7) FIG. 6 illustrates spheres with different numbers of bins;

(8) FIG. 7 shows two-dimensional projections of incident and azimuth angle pairs obtained using (a) a conventional cylindrical coordinates approach versus (b) an equal area spherical binning according to an embodiment;

(9) FIG. 8 is a flowchart of a method of seismic exploration according to an embodiment; and

(10) FIG. 9 is a schematic diagram of a data processing apparatus according to an embodiment.

DETAILED DESCRIPTION

(11) The following description of the exemplary embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Reference throughout the specification to “one embodiment” or “an embodiment” means that a particular feature, structure or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of the phrases “in one embodiment” or “in an embodiment” in various places is not necessarily referring to the same embodiment. Further, the particular features, structures or characteristics may be combined in any suitable manner in one or more embodiments.

(12) The embodiments described in this section process seismic data using RTM with an equal area spherical binning technique that defines reflection angle bins at any reflection point. The workflow can be summarized as follows: 1) Initiate partial images I.sub.p associated to bi-angular bins defined for reflection points using an equal area spherical binning technique; 2) Compute forward and backward propagated wavefields; 3) Apply imaging condition in any given domain. For each point in this domain i. Compute incident and azimuth angles between the forward and backward wavefields computed in (2) (and transformed in another domain if needed); ii. Identify center p of a bi-angular bin at the reflection point, p being closest to the pair of incident and azimuth angles computed in (i); iii. Add data point to the partial image I.sub.p initialized in (1) with the center at p as identified in (ii); 4) Optionally, interpolate to an arbitrary set of points on the sphere.

(13) The bin definition in 1) is independent of locations in the subsurface or seismic data but the bin definition is selected to best fit the data acquisition. The bins are bi-angular in the sense that each bin covers a range in θ (incidence angle at a reflection point) and a range in φ (azimuth angle at a reflection point). In other words, step 1) splits the RTM image of the underground formation I(x) in partial images I.sub.p(x) associated with equal area spherical bins.

(14) The steps 2) and 3) are repeated for each shot or any given subset of the seismic data used to compute the forward and backward wavefields. In the case of a streamer acquisition this subset is usually taken to be the collection of traces corresponding to a single shot position, while in the case of OBN acquisitions it is taken as the collection of traces for a single receiver position. In the OBN case, the forward wavefield is calculated from the receiver position and the backward wavefield is calculated from the source positions. This swap of source and receiver positions is possible due to the principle of reciprocity.

(15) As illustrated in FIG. 5, first, partial images I.sub.p are initialized at 510. Each of the partial images corresponds to a specific one of the equal area spherical bins (the center of the bi-angular bin being θ*, ϕ*):
I.sub.p(x):=I(x;θ*,ϕ*).  (2)

(16) An equal area spherical binning algorithm, Hierarchical Equal Area Isolatitude Pixelization, was described in the 2005 article “HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere” by Gorski et al., published in The Astrophysical Journal, 622(2), pp. 759-771 (which is incorporated herein by reference in its entirety). HEALPix is a binning procedure based on three integer parameters: N.sub.θ, N.sub.φ and N.sub.side. N.sub.θ and N.sub.φ define a base configuration. N.sub.side defines the level of refinement of each base configuration (for the base configuration N.sub.side=1), N.sub.side indicating the number of divisions along the side of a base-resolution pixel. FIG. 6 (corresponding to FIG. 3 from Gorski's article) illustrates spheres binned using HEALPix with different values of N.sub.side.

(17) In one embodiment, pixel (center of the bin) positions for N.sub.θ is given by an odd integer (in tests 3 has been used). Note that even integer is also possible in other embodiments. With odd N.sub.θ, the bins are divided in two regions, according to whether they are close to the pole or to the equator and are indexed by two integers: the first integer defines which ring of constant latitude the bin belongs to, while the second identifies the bin inside the ring. For the polar regime, the ring number i ranges between 1 and (N.sub.side−1) and for each ring the second number j ranges between 1 and N.sub.φi. A combined pixel number can be defined as p=j−1+N.sub.φi(i−1)/2. The bin center positions for the polar regime are given in terms of ϕ and z≡cos θ:

(18) z = 1 - ( 1 - z * ) i 2 N s i d e 2 ; ϕ = 2 π N φ i ( j - 1 2 ) . ( 3 )
where z*=(N.sub.θ−1)/N.sub.θ defines the separation between polar and equatorial regimes. A similar expression can be written for the equatorial regime. In the equatorial regime, the ring number i takes values from N.sub.side to . . . N.sub.side(N.sub.θ+1)/2 and the second number j is between 1 and N.sub.φN.sub.side. The combined pixel number is given by p=N.sub.φN.sub.side(N.sub.side−1)/2+N.sub.φN.sub.side(i−N.sub.side)+j−1. Bin positions are given by:

(19) z = z * ( 1 - 2 ( i - N side ) N side ( N θ - 1 ) ) and ϕ = 2 π N φ N ( j - s 2 ) with s = ( i - N side + 1 ) mod 2. ( 4 )

(20) The transition between polar and equatorial regimes is given by i=N.sub.side, in which case z=z*. The number of azimuths in the ring closest to the pole is given by N.sub.φ. This number increases in the polar region until it reaches the equatorial regime, and then it remains constant. The total number of pixels in the northern hemisphere, including the equator, is given by:
N.sub.pix.sup.H=N.sub.φ(N.sub.side(N.sub.side−1)/2+N.sub.side.sup.2(N.sub.θ−1)/2+N.sub.side).  (5)

(21) The total number of pixels on the full sphere is given by:
N.sub.pix=N.sub.φN.sub.θN.sub.side.sup.2.  (6)

(22) Similar to FIG. 4, FIG. 7 shows two-dimensional projections of incident azimuth angle pairs obtained using (a) a conventional cylindrical coordinates approach versus (b) an equal area spherical binning according to an embodiment. The equal area spherical binning provides a more equal distribution of the bin centers and is therefore better suited to describe a full azimuth seismic data acquisition such as the one using OBN.

(23) Returning now to FIG. 5, after a partial image is associated with each bi-angular bin at 510, for each reflection identified in the seismic data, forward and backward propagated wavefields (i.e., p.sub.F and p.sub.B) are calculated at 520. Then, on one hand, azimuth and reflection angles θ and ϕ are computed at 530 to then identify the spherical bin (i.e., the closest pair θ*, ϕ* to θ and ϕ) at 540. On the other hand, step 520 also allows applying an imaging condition for a given domain point at 550. The image condition may be applied in space time domain or in another domain such as in wavenumber domain (in the latter case, the wavefields being transformed to and back this other domain).

(24) At 560, the result of the imaging condition is added to the partial image I.sub.p corresponding to the identified bin at that particular data domain point.

(25) Optionally, at 570, one may interpolate any contribution to refine the angle gather outputs for any point on the sphere.

(26) The image I of the underground formation is the sum of all partial images I.sub.p:
I(x)=Σ.sub.pw.sub.pI.sub.p(x)  (10)
weights w.sub.p being constant and equal to 4π/N.sub.pix if the bins have equal areas, but having to account for the area differences if the original equal area bins have been converted to other bins at 570.

(27) The binning procedure explained relative to step 510 defines a partition of the (hemi-)sphere, which can in turn be used to approximate integrals of functions defined on it. In particular, an interpolation formula may yield the value of a function at an arbitrary point of a bin on the sphere, given the function's value for the bin's center. One version of the interpolation formula is:

(28) f ( r ^ q ) .Math. p f ( r ^ p ) t qp with ( 7 ) t qp = 1 N pix .Math. l = 0 l max ( 2 l + 1 ) P l ( r ^ q .Math. r ^ p ) . ( 8 )
where q identifies an arbitrary position on the sphere given by {circumflex over (r)}.sub.q, the set {{circumflex over (r)}.sub.p} is defined on the HEALPix pixel positions and P.sub.l(x) is the Legendre polynomial of degree l. The interpolation matrix t.sub.qp plays the role of the sinc interpolator on the sphere. It is defined as a sum of Legendre polynomials that depend only on the dot product between {circumflex over (r)}.sub.q and the set of {circumflex over (r)}.sub.p. The sum over l in the definition of t.sub.qp is bounded by a parameter l.sub.max that depends on the particular binning chosen (i.e., on the combination N.sub.θ, N.sub.φ, N.sub.side) and needs to be estimated. An estimate of l.sub.max is proposed in the 1996 article “An icosahedron-based method for pixelizing the celestial sphere” by M. Tegmark (published in The Astrophysical Journal Letters; 470 (2), L81, which is incorporated herewith by reference in its entirety):
l.sub.max≈√{square root over (3N.sub.pix)}.  (9)

(29) As an option, a percentage of this estimate may be used and tested for each binning configuration and data.

(30) In real and simulated data tests, different values sets (N.sub.θ, N.sub.φ, N.sub.side) have been used for a rectangular grid of OBNs. These parameters may be adjusted to adapt for different acquisition geometries (for instance, hexagonal receiver grid versus rectangular grid). N.sub.side defines the level of refinement starting from this base configuration. In our case, since reflection angles are defined up to 90°, in fact only half of the sphere is used, with a total number of pixels N.sub.pix.sup.H, given in equation (5). The (N.sub.θ, N.sub.φ, N.sub.side)=(3, 4, 6) as an average incident angle increment of about 7.6°, while the (N.sub.θ, N.sub.φ, N.sub.side)=(3, 4, 8) has an average incident angle increment of about 5.8°. These configurations were compared with a cylindrical coordinates binning of 5° incident and 9° azimuth sectors. The improvement in angular coverage in each bin was visible, with the configuration (3, 4, 6) showing even better continuity of events while still capturing variations in incident angle and azimuth.

(31) FIG. 8 illustrates a method of seismic exploration 800. Method 800 includes obtaining seismic data representing seismic excitations emitted from a source position and detected at a receiver position after traveling through an explored underground formation, at 810. Data acquisition geometry may be the one illustrated in FIG. 1, but it is not limited to OBN type of acquisition, method 800 being applicable to any type of land and marine data acquisition.

(32) Method 800 further includes processing the seismic data using RTM to generate image angle gathers at 820, wherein an equal area spherical binning technique is employed to define incident and azimuth angles at reflection points, and the image angle gathers. The image angle gathers are used to determine the presence of hydrocarbons and/or other minerals.

(33) If the equal area spherical binning technique is HEALPix, then it is determined by three numbers: one of the three numbers indicates a number of bins on each side of an equal area of a base configuration, and other two of the three numbers characterizing incident and azimuth angle bins, respectively. Method 800 may further include selecting a combination of the three numbers using a subset of the seismic data, the three numbers being varied within predetermined value ranges (e.g., (N.sub.θ, N.sub.φ, N.sub.side) may be within (2-10, 2-12, 4-20)) to identify the combination yielding a best achievable bi-angular binning for a seismic data acquisition system used to acquire the seismic data.

(34) In one embodiment, step 820 includes, for each series of amplitude versus time series recorded at the receiver position following one of the seismic excitations emitted at the source position performing: initializing a set of partial images, computing a forward-propagated wavefield from the source position through the underground formation and a backward-propagated wavefield from the receiver position through the underground formation, applying an imaging condition to estimate a reflection location using the forward-propagated wavefield and the backward-propagated wavefield, computing incident and azimuth reflection angles at the reflection location, and identifying one of the angular bins defined by the equal area spherical binning technique with a bin center closest to the incident and azimuth reflection angles; and adding the data point to the identified angular bin. Here, the incidence and azimuth reflection angles may be calculated in a wavenumber domain, the forward and backward-propagated wavefields being transformed in the wavenumber domain using an anti-leakage Fourier transform, and angle domain image gathers being converted back in a space domain by applying a fast Fourier transform. Alternatively, the incidence and azimuth reflection angles may be calculated using derivatives of the wavefields.

(35) Method 800 may further include interpolating the image angle gathers to generate another set of image angle gathers at positions different from bin centers of bins as defined by the equal area spherical binning technique.

(36) Method 800 may also include using migration velocity analysis to update a model used to compute the forward-propagated wavefield and the backward-propagated wavefield.

(37) The use of equal area spherical binning solves the over-sampling of the region around the pole that plagued conventional cylindrical binning. Since all bins have the same area, area equalization is no longer necessary. In HEALPix, the bins are organized in rings of constant latitude, which simplifies the computation of Legendre polynomials used in the analysis of the binned data.

(38) The above-discussed methods may be implemented in a computing device 900 as illustrated in FIG. 9. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein.

(39) Exemplary computing device 900 suitable for performing the activities described in the exemplary embodiments may include a server 901. Server 901 may include a central processor (CPU) 902 coupled to a random-access memory (RAM) 904 and to a read-only memory (ROM) 906. ROM 906 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 902 may communicate with other internal and external components through input/output (I/O) circuitry 908 and bussing 910 to provide control signals and the like. Processor 902 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

(40) Server 901 may also include one or more data storage devices, including hard drives 912, CD-ROM drives 914 and other hardware capable of reading and/or storing information, such as DVD, etc. In one embodiment, software for carrying out the above-discussed steps may be stored and distributed on a CD-ROM or DVD 916, a USB storage device 918 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as CD-ROM drive 914, disk drive 912, etc. Server 901 may be coupled to a display 920, which may be any type of known display or presentation screen, such as LCD, plasma display, cathode ray tube (CRT), etc. A user input interface 922 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.

(41) Server 901 may be coupled to other devices, such as sources, detectors, etc. The server may be part of a larger network configuration, as in a global area network (GAN) such as the Internet 928, which allows ultimate connection to various computing devices.

(42) The spherical binning may also be used with one-way wave equation migration (OWEM) instead of RTM and with Kirchhoff migration. The only modification of the methods illustrated in FIGS. 5 and 8 for OWEM is in the manner of propagating the forward and backward wavefields. That is, the OWEM forward- and backward-propagated wavefields are computed using a downward continuation method such that the image corresponds to a one-way wave equation migration. This difference between the OWEN and RTM wavefields is described in the 2004 article “A comparison between one-way and two-way wave-equation migration” by W. A. Mulder and R.-E. Plessix published in Geophysics, Vol. 69, No. 6, pp. 1491-1504, which is incorporated herein by reference in its entirety.

(43) Kirchhoff migration corresponds to an asymptotic version of RTM, for large frequencies. Kirchhoff angle gathers have been introduced in the 2001 article “Common-angle migration: A strategy for imaging complex media” by S. Xu et al. published in Geophysics, Vol. 66, No. 6, pp. 1877-1894, which is also incorporated herein by reference in its entirety. Because of the asymptotic approximation, the implementation in the Kirchhoff case is very different. Instead of propagating wavefields from receiver and source locations, ray-tracing is used to compute travel-times between source or receiver location to any location at the sub-surface.

(44) As far as binning is concerned, Kirchhoff migration does not cause major differences relative to RTM. The directions of the rays coming from source and receiver locations are available through ray-tracing, which allows computing incident or azimuth angles at any sub-surface location, angle that are then input for the spherical binning.

(45) The embodiments described in this section provide methods and apparatuses using equal area spherical binning for generating angle gathers. It should be understood this description is not intended to limit the invention. On the contrary, the exemplary embodiments are intended to cover alternatives, modifications and equivalents, which are included in the spirit and scope of the invention. Further, in the detailed description of the exemplary embodiments, numerous specific details are set forth to provide a comprehensive understanding of the invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

(46) Although the features and elements of the present exemplary embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein.

(47) This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. Other examples that occur to those skilled in the art are intended to be within the scope of the disclosed inventions.