S-wave anisotropy estimate by automated image registration

Abstract

The present disclosure provides a system and method for estimating fracture density within a subsurface formation from S-wave seismic data. In one embodiment, the S-wave seismic data is separated into fast (“S.sub.1”) and slow (“S.sub.2”) data. A computer is used to compute local similarity of the S.sub.1 and S.sub.2 data and to compute a cumulative time-difference by which the S.sub.2 data lags the S.sub.1 data from the local similarity. Based on the computed cumulative time-difference, the fracture density of a subsurface formation is estimated.

Claims

1. A method, comprising: obtaining S-wave seismic data of a subsurface formation, wherein the S-wave seismic data is acquired from the subsurface formation using multi-component receivers adapted to measure a plurality of motion vector components including two horizontal components, the S-wave seismic data comprising a plurality of seismic time samples; separating the S-wave seismic data into fast (“S.sub.1”) and slow (“S.sub.2”) data; using a computer to compute local similarity of the S.sub.1 and S.sub.2 data at each time sample; using a computer to compute at each time sample a cumulative time-difference by which the S.sub.2 data lags the S.sub.1 data from the local similarity, wherein computing the cumulative time-difference includes, performing a local similarity attribute scan between the S.sub.1 and S.sub.2 data; generating a warping function, w(t), by detecting areas of similarity in the local similarity attribute scan, and determining the cumulative time-difference, ΔT.sub.S.sup.cum(t), based upon w(t); calculating an instantaneous S-wave time-difference, δt.sub.S.sup.ins(t), at each seismic time sample from the cumulative S-wave time difference; estimating fracture density within the subsurface formation from the cumulative time-difference, ΔT.sub.S.sup.cum(t); and drilling a well into the subsurface formation based at least in part on the fracture density.

2. The method of claim 1, wherein the cumulative time-difference, ΔT.sub.S.sup.cum(t), is determined using a formula that can be expressed as
ΔT.sub.S.sup.cum(t)=(w(t)−1)t, where t is reflection travel time.

3. The method of claim 1, further comprising: generating time-compensated S.sub.2 data based upon the .fwdarw.T.sub.S.sup.cum(t) to time shift the S.sub.2 data.

4. The method of claim 1, wherein the subsurface formation has a primary fracture direction, the S.sub.1 and S.sub.2 data is generated by vector-rotating the horizontal components to directions parallel and normal to the primary fracture direction.

5. The method of claim 1, wherein the S.sub.1 and S.sub.2 data are generated by performing a layer stripping process on the S-wave seismic data to generate the S.sub.1 data and a registered S.sub.2 data which has a time-difference applied and removing the time-difference applied to the registered S.sub.2 data to generate the S.sub.2 data.

6. The method of claim 1, wherein the instantaneous S-wave time-difference is calculated using a formula that can be expressed as $\delta t_S^{\mathrm{ins}}\left(t\right)=\frac{d\Delta T_S^{\mathrm{cum}}\left(t\right)}{dt}.$

7. The method of claim 6, further comprising: estimating the fracture density based on δ.sub.S.sup.ins(t).

8. The method of claim 1, further comprising: producing hydrocarbons from the well.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the U.S. Patent and Trademark Office upon request and payment of the necessary fee.

(2) The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

(3) FIG. 1 is a flow chart showing basic steps for the fracture characterization according to one embodiment of the present disclosure;

(4) FIG. 2A depicts an exemplary similarity scan and warping function for a seismic pair;

(5) FIG. 2B depicts the cumulative S-wave time difference for the seismic trace information depicted in FIG. 2A;

(6) FIG. 3 illustrates S.sub.1 and uncompensated S.sub.2 traces;

(7) FIG. 4 illustrates S.sub.1 and time compensated S.sub.2 traces;

(8) FIG. 5 illustrates a cumulative S-wave time-difference volume based on the seismic traces shown in FIGS. 3 and 4;

(9) FIG. 6 illustrates an instantaneous S-wave time-difference volume based on the cumulative S-wave time-difference volume shown in FIG. 5.

(10) The invention will be described in connection with example embodiments. However, to the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claims. Persons skilled in the technical field will readily recognize that in practical applications of the present inventive method, it must be performed on a computer, typically a suitably programmed digital computer.

DESCRIPTION OF THE SELECTED EMBODIMENTS

(11) The disclosed methodology describes a process to generate a time-continuous time-difference volume between the fast S-wave (S.sub.1) and slow S-wave (S.sub.2) modes. The time-difference between S.sub.1 and S.sub.2 modes is an indicator of fracture intensity within a subsurface formation. As a result, a larger time-difference indicates a larger fracture intensity or density. This new, additional information will help to locate highly anisotropic or fractured zones: critical information for well planning in unconventional reservoirs. As used throughout this document, a time-continuous volume is understood mean an attribute volume in which there is a defined value at each time sample.

(12) The flowchart of FIG. 1 will be referred to in describing one embodiment of the present disclosure. The depicted process (100) first requires acquisition of multi-component, multi-azimuth data (101). Azimuth is defined for a particular source-receiver combination. The direction (relative to true North or some other reference direction) of the line connecting a source-receiver pair is called the azimuth of that particular source-receiver pair and associated seismic data. Traditionally, only the vertical component of the seismic wavefield, which is dominated by the P-wave energy, is acquired. For certain applications, all three vector components of the wavefield are also acquired (using a motion-detector type of seismic receiver).

(13) In this type of acquisition, a P-wave source is used which may be, but is not limited to, dynamite, a vertical vibrator, or air gun. The vertical component of the data mostly contains P-wave energy and the two horizontal components carry converted-wave PS energy. The PS energy is defined as the P-wave energy reflected back from a reflector as S-wave energy; i.e., P-wave goes down and some of that energy is reflect back up in an S-wave mode.

(14) The process continues as the acquired PS energy is rotated or resolved into radial (direction along the line connecting source to receiver) and transverse (direction perpendicular to the line connecting source to receiver) components (103). Free-surface related seismic noise, such as surface-waves and free-surface multiples, is then removed from the radial and transverse components (105). After noise correction, normal moveout (NMO) correction is applied (107) on the data to flatten the reflections. In other embodiments, the reflections may be flattened by pre-stack time migration. Though not depicted in FIG. 1, some embodiments may include the additional step of stacking the flattened data. By stacking the data, the signal-to-noise ratio may be improved. As appreciated by those skilled in the relevant art, steps 101-107 are standard processing steps and are routinely applied in seismic data processing.

(15) The two horizontal components of the recorded wavefields, which primarily have S-wave reflected energies, carry a mixture of fast and slow S-waves. The S-wave data recorded over a fractured subsurface is seldom separated in pure S.sub.1 and S.sub.2 modes. In order to do so, the horizontal components (radial and transverse components) first need to be rotated to fast and slow S-wave directions before S.sub.1 and S.sub.2 images can be time-aligned. The way this is done is dependent upon whether there is a single dominant fracture direction in the subsurface. If there is only a single dominate fracture direction in the subsurface (109), the horizontal components can be simply vector-rotated to fracture strike and normal directions to produce S.sub.1 and uncompensated S.sub.2 seismic sections (111).

(16) However, if the fracture orientation varies with depth (113), it is not possible to apply a single vector rotation on the whole section. In this case, layer-stripping is performed to produce an S.sub.1 seismic section and a registered (or time-compensated) S.sub.2 section (115). A number of methods have been published on layer stripping from surface seismic data. To name a few, Gaiser [9, 14] published a method called “3-D converted shear wave rotation with layer stripping”. Another method was published by Thomsen et al., [11] called “coarse layer stripping of vertically variable azimuthal anisotropy from shear-wave data”. Granger et al., [12] developed a method to find the fast S-wave direction which corresponds to the fracture orientation. Haacke et al. [17] proposed a method of layer-stripping in marine data. Crampin [13] gave a detailed description of S-wave propagation in fractured media which led to development of the layer-stripping technique. Bansal et al. [18] developed a method to perform true-amplitude layer-stripping using surface seismic PS data. Bansal's method performsf a scan of offset and azimuths of the surface seismic data to obtain the optimum dataset to perform layer stripping.

(17) Conventional layer-stripping techniques typically consist of four primary steps: (1) division of data in several layers (or windows), (2) rotation of the radial and transverse components to fracture strike and normal in first window to generate fast and slow S-waves, (3) cross-correlation of the fast and slow S-waves to estimate the time-difference between the two in the first window and (4) shifting the slow S-wave by the estimated time-difference to align the fast and slow S-wave events or reflections in first window. This process is performed multiple times in a top-down fashion revealing fracture orientation and S-wave time-difference of each subsequent fractured layer (or window). This process produces pure S.sub.1- and registered S.sub.2-mode seismic sections, with the registered S.sub.2-mode seismic sections resulting from the time shifting noted above in step 4. This time shifting is necessary so that the fast and slow waves can be accurately determined in the next layer (or window). In one non-limiting embodiment, the method proposed by Bansal et al. [18] is used as the layer stripping technique. However, it must be noted that the embodiments of the present disclosure are not dependent on the type of layer stripping method used.

(18) After the S.sub.1- and S.sub.2-modes have been generated, the process continues by determining the time-difference between the two S-waves. Fomel and Jin [19] previously developed a method to register time-lapse seismic images using the local similarity attribute which provides a smooth continuous measure of similarity between two images. Their solution produces a point-by-point aligned time images, time-continuous time-difference between the two images and a warping function used to align the images. There are two main advantages of this approach over traditional cross-correlation or cross-equalization approaches: (1) stability of the method and (2) generation of time-continuous time-lags and warping functions, a direct indicator of velocity change in the reservoir. However, once the S.sub.1 and S.sub.2 data are generated by the present inventive method, the cumulative S-wave time difference may alternatively be calculated by cross-correlation, cross-equalization, or any similar comparison method, but using the similarity attribute is preferred because of, at least, stability and robustness.

(19) One embodiment of the present disclosure uses a similar approach to register S.sub.1 and S.sub.2 images and to produce a time-continuous S-wave time-difference. As appreciated by those of ordinary skill, Fomel teaches an image registration technique to remove artifact time-lapse differences caused by velocity changes. Where Fomel aligns P-wave data at different times, certain embodiments of the present disclosure align simultaneous S-wave data after it has been separated into its S.sub.1 and S.sub.2 modes. Because of the time-continuous nature of the time-difference, the inventive methodology allows for the computation of the S-wave time-difference between any two reflection events. Such information can lead surveyors to highly fractured zone in the subsurface.

(20) In order to perform a similarity-based S-wave time-difference calculation, the window-based S-wave time-difference applied in step 115 on the registered slow S-wave volume is removed in order to generate an uncompensated slow S-wave (S.sub.2) volume (117). A local similarity scan between S.sub.1 and uncompensated S.sub.2 volumes is then performed (119) to find the optimum warping function, which is the name given to a function of time that compensates the S.sub.2 volume so that it is in phase with the corresponding S.sub.1 volume. The warping function also indicates the extent to which the cumulative S.sub.2 time shift departs from a linear function of time for a particular fracture system. The correct warping function (or warping function trend), w(t), is generated by detecting the areas of strong similarity (121) and is calculated to minimize the difference between the S.sub.1 and uncompensated S.sub.2 volumes when applied to the uncompensated S.sub.2 volume. In some embodiments, a smoothness constrain is applied on the warping function while it is being calculated. Any kind of smoothness constrain me be applied on the warping function in time (or vertical) and/or lateral (X and Y) directions. In one embodiment, the user is prompted to provide the number of points in each direction to use in order to smooth out the warping function.

(21) FIG. 2A is an illustration depicting a similarity scan (201) as the background and the chosen warping function (203) for a seismic trace pair in S.sub.1 and uncompensated S.sub.2 seismic sections.

(22) Referring again to FIG. 1, based upon the generated warping function, a cumulative S-wave time-difference is then calculated (123). In one embodiment, the cumulative S-wave time-difference, ΔT.sub.S.sup.cum(t), is determined using the following equation, which follows from the definition of the warping function:
ΔT.sub.S.sup.cum(t)=(w(t)−1)t,  (1)
where t is the reflection travel time. FIG. 2B illustrates the calculated cumulative S-wave time-difference (205) as a function of time based on the warping function depicted in FIG. 2A.

(23) Based upon the cumulative S-wave time-difference, an instantaneous S-wave time difference is calculated (125). The fracture density of a subsurface formation is directly proportional to the instantaneous S-wave time difference. In one embodiment, the instantaneous S-wave time-difference, δt.sub.S.sup.ins(t), is determined at each seismic data time sample using the following equation:

(24) $\begin{array}{cc}\delta t_S^{\mathrm{ins}}\left(t\right)=\frac{d\Delta T_S^{\mathrm{cum}}\left(t\right)}{dt}& \left(2\right)\end{array}$

(25) The described process may then be performed for every trace pair in the volume, thereby resulting in a cumulative S-wave time-difference volume and an instantaneous S-wave time-difference volume being generated. In another aspect of the present disclosure, a new compensated S.sub.2 volume may be generated by using ΔT.sub.S.sup.cum(t) volume to time-shift uncompensated S.sub.2 volume (127).

(26) It is important to note that the steps depicted in FIG. 1 are provided for illustrative purposes only and a particular step may not be required to perform the inventive methodology. As a non-limiting example, a survey may decide not to calculate an instantaneous S-wave time difference and also forego the generation of a new compensated S.sub.2 volume based upon survey objectives and/or computer processing speed or memory. Again, the claims, and only the claims, define the inventive system and methodology.

EXAMPLE

(27) To demonstrate applicability of the disclosed methodology, the present inventors applied the method on synthetic seismic modeling data. FIG. 3 shows S.sub.1 traces (301) in black and uncompensated S.sub.2 traces (303) in red. FIG. 3 illustrates that events in S.sub.1 and uncompensated S.sub.2 volumes are not aligned in time.

(28) FIG. 4 illustrates the S.sub.1 traces (401) in black and the new compensated S.sub.2 traces (403) in red. As shown, the events in S.sub.1 and compensated S.sub.2 volumes are now aligned. FIGS. 5 and 6 depict the Δ.sub.S.sup.cum(t) and δt.sub.S.sup.ins(t) volumes, respectively. This information allows an operator or surveyor to compute an S-wave time lag and/or instantaneous S-wave time difference between any two reflection events. This new, additional information will help to locate highly anisotropic or fractured zones: critical information for well planning in unconventional reservoirs.

(29) A review of FIGS. 5 and 6 shows that the fracture intensity of the subsurface formation appears to be large in an area corresponding to a travel time between 3 and 3.5 seconds. More specifically, the area of greatest fracture intensity is identified at (501) and (601) in FIGS. 5 and 6, respectively.

(30) While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiments have been shown and described and that all changes and modifications that come within the wording of the claims are desired to be protected. It is also contemplated that structures and features embodied in the present examples can be altered, rearranged, substituted, deleted, duplicated, combined, or added to each other. The articles “the”, “a” and “an” are not necessarily limited to mean only one, but rather are inclusive and open ended so as to include, optionally, multiple such elements.

REFERENCES

(31) 1. Aguilera, R., “Naturally fracture reservoirs”, PennWell Book, Tulsa, 1995. 2. Nelson, R. A., “Geologic analysis of naturally fractured reservoirs”, Gulf Publishing Company, Houston, 2001. 3. Sinha, B. K., Norris, A. N. and Chang, S., “Borehole flexural modes in anisotropic formations”, Geophysics 59, 1037-1052 (1994) 4. Schoenberg, M. and Douma, J., “Elastic wave propagation in media with parallel fractures and aligned cracks”, Geophysical Prospecting 36, 571-590 (1988). 5. Ruger, A. and Tsvankin, I., “Using AVO for fracture detection: Analytic basis and practical solutions”, The Leading Edge 16, 1429 (1997). 6. MacBeth, C. and Crampin, S., 1991, “Comparison of signal processing techniques for estimating the effects of anisotropy”, Geophysical Prospecting 39, 357-386. 7. Alford, R. M., “Multisource multireceiver method and system for geophysical exploration”, U.S. Pat. No. 5,343,441 (1994). 8. Winterstein, D. F. and Meadows, M. A., “Shear-wave polarization and subsurface stress directions at Lost Hills field”, Geophysics 56, 1331-1348 (1991). 9. Gaiser, J. E., “3-D converted shear wave rotation with layer stripping”, U.S. Pat. No. 5,610,875 (1997). 10. Tsvankin, I., “Seismic signatures and analysis of reflection data in anisotropic media”, Pergamon, New York (2001). 11. Thomsen, L., Tsvankin, I. and Mueller, M. C., “Coarse-layer stripping of vertically variable azimuthal anisotropy from shear-wave data”, Geophysics, 64, 1126-1138 (1999). 12. Granger, P. Y., Bonnot, J. M., Gresillaud, A. and Rollet, A., “C-wave resolution enhancement through birefringence compensation at the Valhall field”, Society of Exploration Geophysicist Annual Conference, 2001. 13. Crampin, S., “Evaluation of anisotropy by shear-wave splitting”, Geophysics, 50, 142-152 (1985). 14. Gaiser, J. E., “Application for vector coordinate systems of 3-D converted-wave data”, The Leading Edge, 18, 1290-1300 (1999). 15. Alford, R. M., “Shear data in the presence of azimuthal anisotropy: Dilley, Tex.”, SEG Expanded Abstracts, 5, 476-479 (1986). 16. Thomsen, L., “Converted-wave reflection seismology over inhomogeneous media”, Geophysics, 64, 678-690 (1999). 17. Haacke, R. R., Westbrook, G. K. and Peacock, S., “Layer stripping of shear-wave splitting in marine PS waves”, Geophysical Journal International, 176, 782-804 (2009). 18. Bansal, R., Matheney M. and Liu, E., “True—amplitude Layer-stripping in fractured media”, U.S. Patent Application No. 61/484,949. 19. Fomel, S. and Jin, L., 2009, “Time-lapse image registration using the local similarity attribute”, Geophysics, 74, A7-A11. 20. Bakulin, A., Grechka, V., and Tsvankin, I., “Estimation of fracture parameters from reflection seismic data—Part I: HTI model due to a single fracture set”, Geophysics 65, 1788 (2000)