Method for Fracture Detection Using Multicomponent Seismic Data

The present invention relates to the field of geophysical exploration and more specifically to the detection of sub-surface fractures with surface seismic data. The detection of fractures is of paramount importance in so-called "tight" petroleum reservoirs in which the primary determinant of well producibility is the presence of a connected network of fractures to convey fluid into the borehole. As is well known in the prior art, seismic shear waves are among the most sensitive tools available for detecting fractures in hydrocarbon reservoirs.

Previous shear-wave seismic acquisition and processing techniques have often relied on measurements of shear-wave splitting ("birefringence") to detect sub-surface fractures. To elaborate, aligned fractures induce horizontally-transverse anisotropy in the subsurface such that a vertically-incident shear wave splits into a fast mode polarized along the fracture direction and a slow mode polarized perpendicularly to the fractures. Shear waves are seismic vibrations that are polarized perpendicularly (or nearly perpendicularly in the case of an isotropic material) to their propagation direction, i.e., waves for which the vibrations occur perpendicularly to the direction of the wave's propagation. By measuring the shear-wave splitting (which is proportional to the velocity anisotropy), the location and density (the number of fractures per unit volume) of subsurface fracturing may be determined, because higher fracture densities cause a greater amount of shear-wave splitting. The use of four-component seismic acquisition (two orthogonal sources and two orthogonal receivers active during acquisition by both sources) is described in U.S. Pat. No. 4,803,666 to Alford, which is hereby incorporated by reference herein. Alford's technique entails the acquisition of a four-component data matrix to determine the symmetry planes of the medium. Alford describes the use of rotation algorithms to transform the seismic data into a symmetry-plane coordinate system which provides useful information about the orientation of the symmetry planes of the medium while also improving data quality. Knowledge of the orientation of the medium symmetry planes provided, e.g., by Alford rotation, is useful because the direction of open fracturing and that of the symmetry planes usually coincide. Measuring the symmetry plane orientation is therefore usually tantamount to measuring the fracture orientation, a parameter frequently difficult to determine from geological data and useful in reservoir management decisions.

A second consequence of horizontally transverse anisotropy caused by aligned vertical fracturing is a variation in seismic reflection amplitude at boundaries (such as a reservoir) as a function of profile azimuth caused by changes in the intensity of fracturing at the reflecting interface. This physical phenomenon enables fracture intensity to be ascertained with a relatively high vertical resolution by comparing the reflection amplitudes of the fast and slow shear-wave seismic sections. Both of these methods for fracture detection using shear-wave seismic data are described in detail in U.S. Pat. No. 4,817,061 to Alford et al., which is hereby incorporated by reference herein. Alford et al. describes the use of at least one source polarization along each source-receiver azimuth and receivers having matched polarizations.

However, these techniques may produce unsatisfactory results if the fracturing in the reservoir is relatively weak or if more than a single direction of open fracturing is present. Consequently, a method of seismic exploration which is able to characterize fracture intensity in the subsurface in even subtly fractured reservoirs would be advantageous.

Andreas Ruger, in his Ph.D. thesis, "Reflection Coefficients and Azimuthal AVO Analysis in Anisotropic Media", extended the theoretical treatment of split shear waves in the symmetry planes of a horizontally transversely isotropic ("HTI") medium or an orthorhombically anisotropic medium to the case of non-vertical incidence, deriving equations for the plane-wave reflection coefficients as a function of the waves' incident phase angle. This thesis is hereby incorporated by reference herein. AVO stands for amplitude variation with offset. At non-vertical incidence, i.e., by using offset seismic sources and receivers, it is possible to measure two additional amplitude attributes (reflection amplitude refers to the strength of the reflected signal observed at the receivers) in a seismic waveform, its reflection amplitude intercept and slope. The intercept is the projected zero-angle amplitude of the signal and provides information about the change in acoustic impedance and fracture density across the reflecting interface for shear waves. The reflection slope is the slope of a line fitted through the observed amplitudes as a function of the incidence angle of the waves and gives important information about the change in shear wave velocity and fracture intensity at the reflecting interface. This analytical insight into shear-wave reflection at non-normal incidence offers the possibility of obtaining more information about the elastic parameters (which are related to the fracture density) of the subsurface than previously possible if the numerous practical problems associated with the acquisition, processing, and interpretation of non-vertically-incident shear waves can be overcome.

The present invention is a novel method for acquiring, processing, and interpreting shear-wave seismic data for the purpose of detecting fractures and porosity in the subsurface of the earth's crust. The invention uses the change in shear-wave reflection amplitude as a function of the source-receiver distance (the "source-receiver offset" or "offset") to detect the presence and measure the intensity of fractures. The method includes the steps of acquiring shear-wave seismic data in symmetry planes of the subsurface and estimating the intercept and slope of the reflection amplitude with offset of each acquired data set for each time sample. During the aforementioned estimation step, variances (squared standard deviations) are computed for each slope and intercept at each time sample.

These variances are used in the following step to form weights for each slope and intercept of each shear-wave seismic reflection amplitude data set at each time sample. The weights and each slope and intercept value are then optimally combined at each time sample in a weighted least squares solution for a fracture density or intensity term .DELTA..gamma., a shear-wave velocity contrast term .DELTA..beta./.beta., and a density contrast term .DELTA..rho./.rho.. Contrast terms are the fractional change in a quantity at a reflection interface. The resulting seismic sections of each of these parameters are then interpreted to indicate areas of anomalously high fracture density in the direction of primary fracture orientation and other directions. The density contrast term .DELTA..rho./.rho. is optionally interpreted to determine areas of higher porosity in the subsurface. Changes in the shear-wave velocity contrast term may indicate fractures that are open in a different direction from the predominant orientation. Consequently, this parameter is useful to ascertain whether or not two sets of open fractures are present.

FIG. 1 is a plan view of an azimuthally-anisotropic subsurface model, indicating the direction of assumed vertical fracture orientation and diagrammatically illustrating a shear-wave seismic survey layout designed for sensing the intensity of vertical fractures.

FIG. 2 is a plan view of an alternate embodiment of the seismic survey layout illustrated in FIG. 1, indicating the direction of fracture orientation and diagrammatically illustrating a second shear-wave seismic survey layout for sensing the intensity of vertical fractures in the subsurface.

FIG. 3 is a block diagram illustrating the steps utilized in processing the shear-wave seismic data for detecting fractures according to the technique of the present invention.

FIG. 4 is a sample seismic section of the shear-wave splitting parameter .DELTA..gamma. generated by the methods of the present invention.

FIG. 5 is a map view of the amplitude of a .DELTA..gamma. stack generated by the methods of the present invention over a known, fractured petroleum reservoir indicating the intensity of open fractures in the direction of maximum principal horizontal compressive stress.

FIG. 6 is a map view of the amplitude of the attribute .DELTA..beta./.beta. generated by the technique of the present invention over a known, fractured petroleum reservoir indicating the intensity of fractures open in directions other than the direction of maximum principal horizontal compressive stress.

1 . Brief theoretical background for subsurface fracture detection using shear-wave seismic data.

The presence and intensity of subsurface fracturing often exerts a strong control on the productivity of oil and gas wells in a variety of settings. Consequently, the detection of such fractures by remote means is of considerable industrial importance. As disclosed in U.S. Pat. No. 4,817,061, it has long been known that seismic shear waves (i.e. seismic energy polarized at angles which are perpendicular or nearly perpendicular to the direction of energy propagation) are very sensitive to the orientation and intensity of fractures in the subsurface of the earth. To interpret the response of these shear waves to fracturing, however, it is necessary to have a theoretical model of how the fractures influence the propagation of shear waves. Aligned fractures cause the host rock to behave in an anisotropic fashion, in which the velocities of both compressional and shear waves are a function of the angle of propagation of energy through the medium. The simplest and most common anisotropic model of fracturing is that of horizontal transverse isotropy (HTI), which can be caused, inter alia, by a single set of aligned, penny-shaped vertical fractures embedded in an isotropic rock matrix.

As shown in FIG. 1, a horizontal slice of an azimuthally-anisotropic subsurface model, this anisotropic system imposes two symmetry planes on the subsurface. The first of these is the so-called "isotropy plane" or fracture-strike plane 2, which is a vertical plane along the direction of the (aligned) subsurface fracture set. Seismic waves that have a direction of propagation confined to this plane behave isotropically, i.e., with no variation in velocity as a function of propagation angle. The second plane, perpendicular to the fracture strike plane 2, is the so-called "symmetry plane" 4. It has long been known that information about the intensity of fracturing can be had by conducting shear-wave seismic surveys along the symmetry planes of the medium. The orientation of the symmetry planes may be determined by processing four-component seismograms by the method of Alford rotation described in U.S. Pat. No. 4,803,666 or by using other information such as borehole breakout or fracture orientation measurements from well-log data or azimuthal variations in the reflection signature of P-wave data, as is well known to those skilled in the art. The presence of aligned fractures gives rise to the well-known phenomenon of shear-wave splitting, in which two shear-waves propagate in the medium, a fast shear wave polarized along the direction of the fractures and a slow shear wave polarized perpendicularly to the fractures. The fractional difference in velocity between the two shear waves .gamma. is roughly proportional to the intensity of fracturing in an HTI medium. The acquisition, processing, and interpretation of shear-wave seismic surveys in such HTI media are discussed in Multicomponent Seismology in Petroleum Exploration, which is hereby incorporated by reference herein. The theoretical aspects of seismic wave propagation are discussed in Andreas Ruger's Ph.D. thesis, "Reflection Coefficients and Azimuthal AVO Analysis in Anisotropic Media", which is also hereby incorporated by reference herein. Additionally, as discussed in U.S. Pat. No. 4,817,061, the normal-incidence shear-wave reflection coefficients of the two modes differ by the change in the shear-wave splitting parameter .DELTA..gamma. across the reflecting interface. Consequently, fracture detection may be performed over an HTI medium by examining either the difference in arrival time of the fast and slow shear-wave modes or by comparing the reflection amplitudes of the fast and slow shear-wave seismic sections at the level of interest. The use of the latter technique for fracture detection is described in "Prediction of lateral variability in fracture intensity using multicomponent shear-wave surface seismic as a precursor to horizontal drilling in the Austin Chalk", by M. C. Mueller, which is hereby incorporated by reference herein. It is important to understand that deviations from the assumptions used in such fracture detection algorithms, such as two or more aligned crack sets, dipping rather than vertical fracture sets, or an anisotropic rock matrix containing the cracks, will cause such techniques to produce erroneous results if these complications are not taken into account in processing the seismic data.

Shear-wave seismic surveys in the symmetry and isotropy planes of an azimuthally-anisotropic medium are described in U.S. Pat. No. 4,817,061. In the method described in that patent, a shear-wave survey line 6 aligned along the symmetry plane of the medium is laid out with paired in line and crossline sources 8 and 10 respectively. These are denoted as S.sub.1.sup.sym and S.sub.2.sup.sym in FIG. 1. These sources are matched with in line and crossline receivers 12 and 14 respectively. Another seismic line 16 may be situated along the fracture strike direction with in line and crossline sources 18 and 20 and matched receivers 22 and 24 as well. The seismic acquisition is conducted by activating the in line source 8 and recording on in line receiver 12. The matched crossline source-receiver pair 12 and 14 may then be activated, and the matched in line source-receiver pair 18 and 22 and crossline source-receiver pair 20 and 24 may also be activated if a strike-line survey is also being conducted. These steps are repeated at a plurality of locations along the seismic lines to create multichannel common-midpoint lines of seismic data containing a plurality of source-receiver offsets. Four-component data acquisition (e.g. recording on both receiver pairs 12 and 14 and 22 and 24 during activation of each source 8, 10, and 28 and 20) may also be done for quality control purposes or to more accurately establish the orientation of the symmetry and isotropy planes 4 and 2. The resulting seismograms are than processed in a conventional way well understood by those skilled in the art by, for example, applying amplitude recovery corrections, static time corrections in both shot and receiver domain, normal-move out corrections, and any amplitude-preserving noise suppression to emphasize the primary seismic reflections which are of interest while removing undesirable noise from the seismograms. The seismograms are usually then stacked and the reflection amplitudes and travel times of the fast and slow shear-wave seismic sections are then compared, as described in detail in U.S. Pat. No. 4,817,061, which is incorporated by reference herein.

The method for fracture detection disclosed in U.S. Pat. No. 4,817,061 encompasses near-normal-incidence illumination of the subsurface, e.g., with closely-spaced sources and receivers. While it is possible to extract much useful information regarding the fracture properties of the subsurface with such seismic surveys, it is advantageous to utilize variations in the shear-wave seismic reflectivity of the subsurface as a function of offset (or, equivalently, incidence angle) to obtain more information about fracturing. Andreas Ruger, in his Ph.D. thesis, "Reflection Coefficients and Azimuthal AVO Analysis in Anisotropic Media", derived linearized plane-wave reflection coefficient equations for the S.sub.1.sup.sym, S.sub.2.sup.sym, S.sub.1.sup.strike, and S.sub.2.sup.strike shear waves in an HTI medium. These equations are valid if the contrasts in velocities and anisotropies between the incident and reflecting media are small and if the incident angle is small or moderate (less than the critical angle and generally less than 30-40 degrees). The plane-wave reflection coefficients are: ##EQU1## It can be seen from these equations that much useful information is contained in the non-normal-incidence or slope reflection terms (the terms which are a function of the incidence phase angle .theta.). The other terms, i.e., the terms that are not a function of the incidence phase angle .theta., are the intercepts. Techniques for estimating the normal-incidence reflection coefficient and slope term (the term which is a function of the incident phase angle .theta.) from seismic data are described in detail in Offset-dependent reflectivity: Theory and Practice of AVO analysis, which is hereby fully incorporated by reference herein. Heretofore, seismic AVO analysis has focused on estimating the AVO slope and intercept on conventional, P-wave seismic data for detecting hydrocarbons and lithology changes, rather than from shear waves for fracture detection. However, estimating shear-wave AVO intercepts and gradients involves the same procedure of line fitting as will be described in the following section. The present invention puts these computed shear-wave AVO intercepts and slopes to the novel use of weighted least squares fracture detection. It is apparent from Equation I that if all AVO intercepts and slopes are known, the resulting system of equations contains more equations than unknowns. Because the system is redundant, it is advantageous to use a weighted least squares solution technique described below to estimate the elastic constants .DELTA..beta./.beta., .DELTA..rho./.rho., and the anisotropy parameter .DELTA..gamma..

2. Seismic fracture detection using weighted least squares AVO analysis.

The method of the present invention for detecting fractures using shear-wave AVO analysis may be performed using conventional azimuthal shear-wave seismic surveys illustrated in FIG. 1, or, in an alternate embodiment of the invention illustrated in FIG. 2, by laying out sources 30 and 32 and receivers 34 and 36 in lines 38 and 40 along the symmetry and isotropy planes respectively. A single source (30 or 32) is oriented along or perpendicular (the configuration illustrated in FIG. 2) to the line and the receivers may be oriented at an angle (respectively .theta..sub.1 for receiver 34 and .theta..sub.2 for receiver 36) to the source. The receivers may also be aligned with the sources. If the receivers are oriented at an angle .theta..sub.1 to the source, the component of the reflection amplitude in the direction of the source polarization is

In the embodiment illustrated in FIG. 2, both the symmetry plane and isotropy plane surveys utilize SH (i.e., transverse to the direction of the line) source polarizations. If the azimuthal orientation of the symmetry planes is known from other information, it is not necessary to conduct either four-component acquisition (as described in U.S. Pat. No. 4,803,666) or to match the polarization of sources and receivers (as described in U.S. Pat. No. 4,817,061). As described in detail below, the use of transverse polarization along both directions suffices to both determine the intensity of subsurface fracturing oriented along the fracture strike direction 4 and provide some information about potential fracturing in other directions as well. As is well understood by those skilled in the art, the source and receiver lines may be laid out in 3-D seismic acquisition geometries to maximize coverage of the subsurface in all three coordinates.

A schematic diagram of the seismic data processing flow utilized in the preferred embodiment of the present invention is illustrated in FIG. 3. The data acquisition step 50 described above includes acquiring shear-wave reflection data in the symmetry planes of the medium, which may either be known a priori or determined by Alford-style rotation analysis. As described above, the acquisition of a single component (e.g. the SH mode) of shear-wave polarizations in each of the symmetry planes is sufficient for the fracture detection method of the present invention to succeed. However, if added redundancy in the form of additional measurements is desired, it is also possible to acquire both source and receiver polarizations (e.g. 18, 20, 22, and 24 in FIG. 1) concurrently.

Subsequently to acquisition of the seismic data, the measurements are preprocessed for AVO analysis in a manner well known to those skilled in the art. Such processing may include, for example, the steps of true amplitude recovery, surface-consistent statics and deconvolution, velocity analysis and normal move out correction, Q compensation, and prestack time migration. After this preprocessing, the data are sorted into the common midpoint domain and their AVO intercepts and gradients (for each shear-wave mode acquired) estimated by linear regression well understood by those of ordinary skill in the art, as illustrated in step 52. In the preferred embodiments of the present invention, the incident phase angle .theta. is computed by either raytracing through an interval velocity model of the subsurface or assuming a straight raypath connecting source and receiver. Concurrently with estimation of each AVO intercept and gradient, variances for each intercept and gradient are computed. The inverse of these variances are assigned as weights (step 54) to the subsequent weighted least squares solution of the plane-wave reflection coefficients given in Equation 1.

Because there are typically more equations than unknowns in the plane-wave reflection coefficient equations (Equation 1) for the seismic survey acquisition techniques illustrated in FIGS. 1 and 2, a weighted least squares solution of these equations may advantageously be made to both utilize all available information on the offset-dependent reflection coefficient acquired during seismic surveying and to optimally estimate the quantities of interest, particularly the fracture density term .DELTA..gamma.. Weighted least squares solution techniques differ from conventional least squares estimation in that the problem to be solved, which may be conveniently written as

(where the vector Y is observed data (here, each shear-wave AVO intercept and gradient), X is the matrix relating Y to the variables of interest B, and e is the vector of noise or uncertainty in the measurements Y), has known values of the noise e. Because the values of e are computed in step 54 when the measurement variances are computed, the weighted least squares solution to the problem is

where B* is the estimate of the desired quantities and V is a diagonal matrix, the entries of which are precisely the inverse weights computed in step 54. A thorough mathematical description of weighted least squares regression analysis is given in Seber, Linear Regression Analysis, which is hereby incorporated by reference herein.

The technique of weighted least squares analysis is particularized to the problem of multi-component AVO analysis as follows. If the AVO intercepts and gradients are written with A.sub.S1.sup.sym being .sup.1 A, A.sub.S2.sup.sym being .sup.2 A, A.sub.S2.sup.strike as .sup.3 A, and A.sub.S1.sup.strike as .sup.4 A, and the same convention is adopted for each gradient, i.e., B.sub.s1.sup.sym being .sup.1 B and B.sub.S2.sup.strike being .sup.3 B, the system of Equation 3 becomes ##EQU2## where .DELTA..beta./.beta. is the fractional change in the fast shear-wave velocity at the interface of interest, .DELTA..rho./.rho. is the fractional change in density, and .DELTA..gamma. is the change in the shear-wave splitting parameter. In Equation 5, the SV-mode AVO gradients, i.e., the terms that are a function of sin.sup.2 .theta., are not employed due to anisotropy-related overburden focusing effects which tend to render the SV gradient measurement unreliable. If the subscripts for each of the weights w.sub.i (computed in step 54 by forming the inverse variance of each AVO intercept and gradient) are designated in the order in which the corresponding AVO intercepts and gradients appear on the left hand side of Equation 5, the solution for the fracture density term .DELTA..gamma. using the weighted least squares Equation 4 is then ##EQU3##

F=w.sub.1 w.sub.2 w.sub.3 +w.sub.1 w.sub.2 w.sub.4 +w.sub.1 w.sub.3 w.sub.5 +w.sub.1 w.sub.4 w.sub.5 +w.sub.2 w.sub.3 w.sub.6 +w.sub.2 w.sub.4 w.sub.6 +w.sub.3 w.sub.5 w.sub.6 +w.sub.4 w.sub.5 w.sub.6

and

Missing AVO attributes, e.g. from unacquired modes or components, are assigned a weight of zero in Equations 6 and 7. In the preferred embodiment, this solution (and the corresponding solutions for the density and fast shear-wave velocity contrast terms) is computed from the AVO intercepts and gradients and their corresponding weights computed at step 54 in a solution step 56 for each time sample. The result is preferably a set of 2-D or 3-D seismic data volumes of the contrast parameters which contain information about the fracture intensity, shear-wave velocity change, and density change. At step 56, the standard deviations for each quantity of interest are also preferably computed. The standard deviation of the fracture intensity .DELTA..gamma. term is the square root of its variance, ##EQU4## Following the computation step 56, displays of each of the contrast parameters .DELTA..beta./.beta., .DELTA..rho./.rho., and .DELTA..gamma. are then made (step 58) and interpreted for the presence of fracturing and/or porosity.

FIG. 4 illustrates one such display, a display 60 of the computed change in shear-wave splitting parameter .DELTA..gamma.. As is well known to those of ordinary skill in the art, events corresponding to reflections from, e.g., hydrocarbon reservoirs, may readily be identified by their time values as interpreted from conventional seismic data. Mapping the amplitude of reservoir events on fracture intensity displays such as the display 60 computed from multicomponent weighted least squares AVO analysis allows lateral changes in the amount of fracturing in the reservoir to be detected. One such map, made from a Permian Basin San Andres carbonate oil reservoir, is shown in FIG. 5, as display 62. An area of anomalously high fracture intensity is indicated at location 64 near a prominent fault bend. Another location of anomalously intense fracturing is indicated to the south of a CO.sub.2 injection well at location 66.

An advantageous feature of the present invention is its ability to distinguish fractures open in the dominant direction from fractures open in other directions. This is accomplished by comparison of the change in splitting parameter .DELTA..gamma. (which detects changes in fracturing open in the dominant direction) with the change in shear-wave velocity .DELTA..beta./.beta.. As is well known to those skilled in the art, the presence of fractures aligned in directions other than the dominant direction causes a decrease in the fast shear wave velocity .beta. with respect to the unfractured rock. It is this decrease that is detected in the contrast term .DELTA..beta./.beta. by the weighted least squares AVO analysis methods of the present invention. Consequently, the presence of significant open fracturing in the dominant direction can be detected by examining displays such as 60 and maps such as 62 of the change in the splitting parameter. The presence of fracturing open in all other directions is indicated by anomalies in the shear-wave velocity contrast term. Both displays may advantageously be combined on a single image, such as the display 68 illustrated in FIG. 6. The display 68 includes a section 70 of the intensity of the change in the fast shear-wave velocity .DELTA..beta./.beta., and is overlain by contours 72 of the splitting parameter term .DELTA..gamma.. Areas where both quantities are relatively large compared to their surroundings often denote areas where fractures are open in both the dominant and other directions. One such area 74 is shown to the south of the injection well CVU-97, where a time-lapse anomaly related to carbon dioxide injection at that well was observed on conventional shear-wave seismic data.

Changes in the density contrast term .DELTA..rho./.rho. computed by the methods described above may be mapped within hydrocarbon reservoirs according to the method of the present invention to detect variations in reservoir porosity in a fashion analogous to that used to detect areas of fracturing in the reservoir.

While Equations 5-8 assume an HTI symmetry of anisotropy in the subsurface, it should be understood that analogous equations may be derived in a straightforward fashion for orthorhombic symmetry as well, using anisotropy parameters developed by Tsvankin in his paper "Anisotropic parameters and P-wave velocity for orthorhombic media", which is hereby incorporated by reference herein. Reflection coefficient equations for such a medium are given in Ruger, "Reflection Coefficients and Azimuthal AVO Analysis for Anisotropic Media". The resulting seismic sections are interpreted for the presence of fracturing in a completely analogous fashion.

Having disclosed the fundamental concept of the present invention, it will be obvious to those skilled in the art that variations or modifications thereto can be obtained without departing from the scope of the claimed invention.