WO1999064896A1

WO1999064896A1 - Seismic data interpretation method

Info

Publication number: WO1999064896A1
Application number: PCT/IB1999/001040
Authority: WO
Inventors: Trygve Randen
Original assignee: Geco As; Schlumberger Canada Limited; Services Petroliers Schlumberger
Priority date: 1998-06-09
Filing date: 1999-06-07
Publication date: 1999-12-16
Also published as: GB0028282D0; GB2353358B; AU3950599A; GB2353358A; NO20006224L; GB9904101D0; NO20006224D0; CA2334011A1

Abstract

A method of identifying a subsurface geologic feature using seismic data including the steps of: subdividing a subsurface area into a plurality of discrete cells; associating the cells with seismic data that image subsurface regions represented by the cells; determining structural characteristics for a plurality of the cells; and grouping together nearby cells based on the structural characteristics to identify the subsurface geologic feature. The grouped cells may also be suitably clustered and/or vectorized. The method is capable of identifying seismic reflector surfaces, fault surfaces, and other types of subsurface geologic features. The method is particularly suited for identifying hydrocarbon reservoir boundaries for use in three-dimensional 'shared earth' reservoir modeling.

Description

SEISMIC DATA INTERPRETATION METHOD

BACKGROUND OF THE INVENTION

This invention relates generally to seismic data interpretation methods and more particularly to a method of identifying subsurface geologic features using seismic data.

The oil industry today recovers only about 35% of the oil it finds. To make the next quantum leap in productivity, where ultimate recovery rates routinely exceed 50%, the process by which oil and gas reservoirs are managed must change. A key tool for these types of improved reservoir management activities is the shared earth model, where seismic data is integrated with all other available data, such as drilling data, well log data, well test data, and production data, to describe the static and dynamic properties of the reservoir.

Most reservoir data is obtained in the immediate vicinity of the well bore and measure conditions that are directly associated with only a small fraction of the total reservoir volume. Seismic imaging is .typically the only method available for identifying the boundary surfaces of hydrocarbon reservoirs, including seismic reflector surfaces and faults. This reservoir structure information is often crucial for effective field development and production planning activities.

Hydrocarbons typically migrate upward from their source through porous subterranean strata until the route is blocked by a layer of impermeable rock, and they accumulate beneath this sealing structure. Geologists divide traps into two types, structural and stratigraphic. Structural traps are formed by tectonic forces after sedimentary rocks have been deposited and generally fall into two categories, anticlines and faults. Stratigraphic traps are most often formed at the time the sediments are deposited. They fall into three categories: pinchouts, most common in stream environments where a channel through a flood plain has been filled with permeable sand that was then surrounded by less permeable clays or silts when the channel moved; unconformities, where a permeable reservoir rock has been truncated and covered by an impermeable layer following a depositional period or a time of erosion; and carbonate reefs, fossilized coral structures and associated deposits that arose from ancient ocean shelves and were overlain by layers of both permeable and impermeable rock.

In both structural and stratigraphic traps, the containing surfaces of hydrocarbon reservoirs invariably consist of interfaces between relatively permiable materials and relatively impermiable materials. The bottom surfaces of hydrocarbon reservoirs typically comprise hydrocarbon/water interfaces located within relatively permiable materials. In many cases, the acoustic impedances of materials on opposing sides of these surfaces will be significantly different and these boundary surfaces will act as effective reflectors of seismic energy. These differences in acoustic impedences provide those involved in hydrocarbon exploration and production activities with the opportunity to remotely identify hydrocarbon reservoir boundary surfaces using reflection seismology.

Faults are formed when tectonic forces break sedimentary rock, and displacements occur along the breaks. A fault may generate a fluid channel through a sealing layer, or seal a permeable layer.

A conventional approach to identifying both seismic reflectors and fault surfaces using seismic data is to view two dimensional cross sections of the seismic data and to manually identify either points that appear to lie on a common seismic reflector or points at which the primary geologic layers appear to be displaced (i.e. a fault). While it may be a relatively straightforward problem to manually interpret seismic reflectors or faults on a single two dimensional seismic section, it becomes an extremely tedious and difficult task to manually identify these geologic features as three dimensional (i.e. planar) surfaces because the number of points that must be identified goes up exponentially.

This is particularly true when identifying seismic reflector surfaces in stratigraphic traps. Exploration seismologists find structural traps far easier to identify than stratigraphic traps in both 2D and 3D seismic data because structural traps are seen as highly dipping reflectors and discontinuities in otherwise smooth reflections. For years, seismic data acquisition and processing techniques have been tailored to accentuate these features, allowing interpreters to concentrate their efforts on faults and anticlines rather than their more subtle stratigraphic counterparts.

As a majority of the major hydrocarbon production areas have already been explored for easily identified structural traps, an increasing emphasis is being placed on locating and evaluating complicated structural and stratigraphic traps. A conventional approach to identifying stratigraphic traps is to perform a manual interpretation of all the primary layers, analyze their termination points, and then to manually classify the type of trap based on the interpreter's training and experience. This is often an extremely tedious and difficult task which requires substantial geological experience to allow prospective stratigraphic traps to be properly identified and analyzed.

The complexity of identifying these types of traps using seismic data is compounded by changes in the type of seismic data being acquired. The increasing use of three-dimensional ("3D") acquisition techniques and high- resolution surveys have dramatically increased the quantity of the seismic data available to be interpreted. At the same time, stable to falling crude oil prices and competitive pressures have required companies to more efficiently utilise their exploration and production expenditures and to reduce the time it takes to reach decisions regarding the development of new fields and prospects.

Fortunately, modern computer systems offer particularly powerful tools for processing and aiding in the interpretation of seismic data. While the processing capabilities of computer systems have been increasing exponentially, the cost of these systems have simultaneously fallen dramatically. Modern computer hardware and software systems are being used successfully to process seismic data and to produce clearer and more detailed maps of the subsurface, which often eases the process of manual seismic data interpretation.

The continuing reliance on the manual interpretation of seismic data, however, both adds significant costs to and introduces substantial uncertainty into the process of seismic data interpretation. Even the most highly regarded seismic data interpreter is capable of making errors in judgement or suffering from lapses in concentration.

An even more pervasive problem with manual seismic data interpretation is that a human being is efficient at working on only one problem at a time and is able to only work for limited periods of time before the quality of the work performed begins to rapidly degrade. If an automated computer-based seismic data interpretation system is used, however, large quantities of seismic data can be subdivided and the task of interpreting each portion of the seismic data can be delegated to a separate computer system which can be run 24 hours a day, 7 days a week until the process is complete. This not only greatly speeds up the seismic data interpretation process, but it also eliminates the variation in results that would be inherent if more than one person was involved with interpreting the seismic data. Prior art techniques for identifying subsurface geologic features using seismic data have typically been limited to identifying features in two- dimensional cross-sections. Some of these systems have required digitized horizons or operator selected points on the subsurface feature, such as U.S. Patent No. 5,229,976, issued to Boyd et al. on July 20, 1993. Other systems have used positive/negative amplitude cross-over points and edge following techniques to track horizons, such as U.S. Patent No. 5,148,494 issued to Keskes on September 15, 1992.

These types of prior art methods have suffered from three primary problems. First, these methods have only identified a linear cross-sectional view of the subsurface feature (instead of providing the planar surface information desired for shared earth reservoir modelling). Second, these identification techniques only utilised two-dimensional input data (limiting the accuracy and effectiveness of these prior art methods). Third, the input seismic data was required to exhibit strong contrasts or the method would fail to properly identify the desired geologic feature.

It is therefore desirable to implement an improved method of interpreting seismic data that overcomes one or more of the problems exhibited by these prior art methods.

An object of the present invention is to provide an improved method of interpreting seismic data.

An advantage of the present invention is that three-dimensional geologic features may be identified and the identification process may utilize seismic data from three-dimensional neighborhoods of cells to identify the geologic feature. Another advantage of the present invention is that it is capable of identifying subtle geologic features.

An additional advantage of the present invention is that the selected cells may be clustered to distinguish between different geologic features.

A further advantage of the present invention is that the identified geologic feature may be represented as a plane or other type of polynomial surface using a vectorizing procedure.

Another advantage of the present invention is that the geologic features identified may comprise hydrocarbon reservoir boundaries for use in three-dimensional shared earth reservoir modeling.

An additional advantage of the present invention is that the method may be substantially automated and may be used to identify subsurface geologic features with little or no manual operator input.

SUMMARY OF THE INVENTION

The present invention involves a method of identifying a subsurface geologic feature using seismic data. The method includes the steps of: subdividing a subsurface area into a plurality of discrete cells; associating the cells with seismic data that image subsurface regions represented by the cells; determining structural characteristics for a plurality of the cells; and grouping together nearby cells based on the structural characteristics to identify the subsurface geologic feature. The grouped cells may be suitably clustered and/or vectorized. The interpretation results may be suitably represented and/or displayed. The invention and its benefits will be better understood with reference to the detailed description below and the accompanying figures.

BRIEF DESCRIPTION OF THE DRAWINGS

Figure 1 is a process flow chart showing steps associated with the inventive method;

Figure 2 is a schematic plan view of a seismic survey vessel and associated equipment acquiring seismic data;

Figure 3 is a schematic diagram of a seismic data interpretation workstation;

Figure 4 is a cross-sectional view through a seismic data volume with highlighted lines showing faults intersecting the section; and

Figure 5 is a three dimensional perspective view of faults shown in Figure 4.

Figure 6 is a cross-sectional view through a seismic data volume with highlighted lines showing seismic reflectors intersecting the section;

Figure 7 is a three dimensional perspective view of seismic reflectors shown in Figure 6;

Figure 8 is a cross-sectional view through a termination attribute cube in accordance with an embodiment of the present invention; and Figure 9 is a cross-sectional view through an attribute cube that highlights chaotic reflector responses in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

Figure 1 is a process flow chart showing various steps associated with the present method of identifying subsurface geologic features using seismic data. In the Acquire Seismic Data step 10, seismic data is obtained from a subsurface area. This seismic data is then processed in the Process Seismic Data step 12 to improve the signal to noise ratio of the data and to produce a more accurate representation of the subsurface geology. The subsurface area of interest is subdivided into a plurality of discrete cells in the Subdivide Area step 14, with each cell representing a different region within the subsurface area. The cells are then associated (i.e. "loaded") with seismic data that image the subsurface regions represented by the cells in the Load Seismic Data step 16.

In the Determine Structural Characteristics step 18, structural characteristics are determined for a number of the cells. Nearby cells identifying the subsurface geologic feature of interest are grouped together in the Group Cells step 20. If desired, the grouped cells can then be clustered, represented mathematically, such as by a planar surface, and/or displayed in the Represent/Display Results step 22. The individual steps associated with the present method will now be discussed in detail.

The activities typically performed in the Acquire Seismic Data step 10 are discussed with reference to the ship and associated equipment shown in Figure 2. In Figure 2, a seismic survey vessel 30 is shown towing three seismic sources 32, such as airguns or arrays of airguns. The survey vessel 30 is also shown towing a pair of seismic streamers 34 that contain large numbers of seismic sensors, such as hydrophones. Acoustic pulses produced by the seismic sources 32 are transmitted through the water and into the geologic subsurface. When the spreading acoustic pulses reach subsurface interfaces where the acoustic impedances of materials on opposing sides of the interface change, a portion of the acoustic energy is reflected and returned toward the surface (i.e. the subsurface interfaces act as seismic reflectors). Seismic sensors within the seismic streamers 34 receive this reflected acoustic energy. The acoustic energy received by the seismic sensors is measured (typically at a regular time-based sampling interval, such as every 2 milliseconds), digitized and transmitted to the survey vessel 30 where the measurement values are recorded. These types of sensed, digitized, and recorded acoustic energy-based measurement values typically comprise the type of seismic data used in connection with the present method.

Preferably, the seismic data is acquired using a seismic data acquisition system that individually digitizes and processes the acoustic energy received by each seismic sensor, such as the type of system described in PCT International Application No. PCT/GB97/02544, International Publication No. WO 98/14800. This type of "single sensor" seismic data acquisition equipment is capable of acquiring seismic data that provides significantly clearer images of geologic features than conventional coarse spatial sampling "hardwired array" types of seismic data acquisition systems.

The seismic data will also preferably be acquired using a high resolution, high fold seismic data acquisition geometry. Using such an acquisition geometry, the seismic data may be acquired at 12.5m by 12.5m, or even 6.25m by 6.25m, bin sizes and the number of independent measurements obtained at each bin (the "fold" of the seismic survey) may be 48 or higher. Seismic data acquired under these conditions will typically allow more consistent and more accurate identifications of geologic features to be made using the present method.

Time lapsed or "4D" seismic data may also be used in connection with the inventive method. Acquiring seismic images of hydrocarbon reservoirs at different times can provide important information for reservoir characterization and monitoring activities. When used with time lapsed seismic data, the present method can be used to identify fluid/fluid interface surfaces in the hydrocarbon reservoirs (oil/water and gas/oil interface surfaces) and to monitor the movements of these surfaces due to hydrocarbon production activities. This type of monitoring can provide invaluable information to a reservoir engineer contemplating alternative production enhancement and remedial well activities.

Other types of seismic data (in addition to conventional pressure/pressure transmission mode seismic data) may also be used with the present method to identify subsurface geologic features, such as pressure/shear transmission mode seismic data, shear/shear transmission mode seismic data and multi-component seismic data.

The vessel and equipment configuration shown in Figure 2 merely illustrates one of a vast number of possible seismic data acquisition systems that may be used in connection with the inventive method in marine, transition zone, or land environments. The applicability and scope of the present method is in no way limited or restricted to the use of one particular type of seismic data or seismic data acquired by one particular type of seismic data acquisition system.

After being acquired, the seismic data is processed in the Process Seismic Data step 12. This processing typically improves the signal to noise ratio of the data and transforms the seismic data into a more accurate representation of the subsurface geology. Typical processing steps may include migration (such as pre-stack depth migration), stacking, and/or filtering (such as K-F or tau-p domain filtering). The use of decomposed and reconstructed seismic data, such as the type of decomposed and reconstructed seismic data disclosed in PCT International Application Number PCT/IB98/00209, International Publication Number WO 98/37437, incorporated herein by reference, may be particularly useful when identifying extremely subtle subsurface geologic features. The seismic data may also be converted from the time domain (i.e. where the sampling interval of the seismic data is a time interval, such as 2 milliseconds) to the depth domain (i.e. where the sampling interval of the seismic data is a distance interval, such as 1 meter). Interpreting seismic data after it has been converted to the depth domain allows the identified features to be more readily integrated into a shared earth reservoir model.

Figure 3 depicts components of a seismic data interpretation workstation computer, including one or more portable storage devices 40, a hard storage device 42, a CPU 44, one or more operator input devices (such as a keyboard 46), and one or more output devices (such as a display 48). The workstation may comprise, for instance, a general purpose Silicon Graphics lndigo2 workstation computer. The portable storage device 40 may comprise magnetic or optical recording media, such as floppy disks, CD- ROMs, or magnetic tapes. The portable storage devices 40 may contain seismic data or computer software instructions that are copied to the hard storage device 42 and that allow the workstation to perform one or more of the steps or processes associated with the present method.

The method steps that follow the Process Seismic Data step 12 are preferably performed using the GeoFrame software package available from GeoQuest, a division of Schlumberger Technology Corporation, Houston, Texas, USA. This type of software package is particularly suited for subdividing the subsurface area of interest into a plurality of discrete cells (the Subdivide Area step 14) and associating the cells with portions of the seismic data that image subsurface regions represented by the cells (the Load Seismic Data step 16).

In contrast to some prior art methods, the cells in the present method preferably contain seismic data associated with a single subsurface point (i.e. the seismic data images a subsurface area nearer the centerpoint of this particular cell than the centerpoint of any other particular cell), not seismic data traces (seismic data associated with a suite or group of spaced apart subsurface points). In many cases, the single seismic data value in the cell will be a zero-offset pressure/pressure transmission mode seismic amplitude value associated with the cell centerpoint. In other cases, more than one seismic data value may be associated with a cell, such as a pressure/shear transmission mode value in addition to the pressure/pressure transmission mode value discussed above, or one or more non-zero offset value associated with the imaging point.

Next, structural characteristics are determined for at least some of the cells in the Determine Structural Characteristics step 18. A choice must typically be made at this point whether to perform a true three-dimensional ("3D") identification of the subsurface geologic features or to identify the features using more simplistic two-dimensional ("2D") methods. Subsurface feature identification using 3D methods provides better results, but is more computationally intensive. A faster way to accomplish similar results is to identify the subsurface features in 2D slices of the data and then to aggregate the results. This substantially reduces the computational requirements of the method, but the output of the process may not be as accurate or reliable. If three-dimensional calculations are made, the structural characteristics determined may comprise or include local seismic reflector dip and azimuth values. If two-dimensional calculations are made, the structural characteristics may comprise or include seismic reflector local plunge estimates (i.e. estimates of the dip angles of the seismic reflectors from the particular cross-sectional reference viewpoint of the section). The confidence or likelihood that the values of these determined structural characteristics are accurate can also be quantified during this process.

3D Structural Characteristics Determination Methods

Various three-dimensional methods for determining the local dominating dip and azimuth of the seismic reflector will now be described. In the following, if a signal, x(t t₂,tj, has a fixed dip and azimuth, then the loci of the signal's iso-value surfaces, i.e. the points where the signal has constant value, will be parallel planes. Furthermore, the gradient of the signal, V , will be perpendicular to these parallel planes.

In a real world dataset, the loci of the iso-value surfaces are not likely to be exact planes. Hence, we attempt to estimate the dominating dip and azimuth. Given a set of 3-dimensional spatial gradient vectors,

d*⁽', - ₃ ⁾

_x( ._{2 s}t ) =

^χ(t, ,t₂ ,t₃)

the dominating dip and azimuth may be determined by principal component analysis. In principal component analysis, the covariance matrix,

of the gradient vectors is constructed, and the dominating dip and azimuth selected to be given by the principal eigenvector of this matrix. The C,'s of the covariance matrix are defined by

where

Having only one dip and azimuth estimate for an entire signal may be of limited practical utility. To increase the usefulness of this estimate, the dominating local dip and azimuth is needed. A localized estimate is obtained by replacing the global covariance estimate by a windowed local estimate,

The window function, w(s) = w(t₁,t₂,t₃), will typically be a low-pass filter. An implication of the windowing function is a smoothing of the dip and azimuth estimate.

In summary, the local dip and azimuth estimation is done by: 1. Orientation selective filtering, typically by gradient estimation;

2. Windowed covariance estimation of the gradient vectors; and

3. Principal component analysis.

There are several possible implementations of the gradient estimation and window function. Below is an exemplification of particularly useful choices of these operators.

Orientation selective filtering: For discrete data, the gradient is a discrete estimate. By nature, gradient estimation is rather noise sensitive. In order to allow some noise removal, the gradient estimation is done by filtering with a derivative of Gaussian ("DoG") filter. If discretization effects are ignored, this is equivalent to smoothing the data by a Gaussian low-pass filter and then differentiation. The smoothing removes noise, and the size of the Gaussian low-pass filter determines the degree of noise removal.

The unit pulse response of a multi-dimensional DoG filter is separable, and has the equation

in the differentiated dimension and

in non-differentiated dimensions, where C_κ and C_λ are scalar factors. A 3D DoG operation consists of applying h_κ(k) once and h_λ(l) twice, to different dimensions. For example applying h_κ(k) vertically and h_λ (l) to the two horizontal dimensions gives the vertical gradient component. The result of applying all DoG filters is a gradient vector, Vx(t, ,t₂,t₃) , with one partial derivative component for each dimension. Other types of orientation selective filters may be used to produce orientation estimates that may be used instead of or in addition to the gradient estimates described above, including bandpass or high-pass filters having orientation selective properties.

Windowed covariance estimation: At strong edges in the data, the dip and azimuth of the gradient, Vx(t,,t₂,t₃) , will be consistent, but in regions of the data with weaker signal, the gradient will have less consistent dip and azimuth. This problem is primarily dealt with by the window function in the covariance estimate, which essentially smoothes the dip and azimuth estimates.

A practical implementation of this operation is to estimate the volume

and then smooth it with a Gaussian low-pass filter,

The samples of the smoothed volume are then used as the estimates of the C,_;'s. This implies the computation of 3x3=9 volumes. Note that the Gaussian filter is separable, which speeds up the computation significantly. Principal component analysis: Now we typically have one covariance matrix C for each position in the data volume. The dominating local dip and azimuth is estimated by calculating the principal eigenvector of C. Detailed information regarding suitable principal component analysis implementation methods may be found in R.A. Johnson and D.W. Wichern. Applied Multivariate Statistical Analysis. 3rd ed. Prentice-Hall International, Inc., Englewood Cliffs, NJ, 1992, incorporated herein by reference.

Estimate reliability measure: There will be three eigenvectors, v, , of the matrix C, each of them associated with one eigenvalue, λ_t . The larger v,. , the better v, describes the dip and azimuth. The larger the difference between the dominating λ,. and the two other λ, 's, the more reliable the dip and azimuth estimate is. Assuming without loss of generality that λ ≥ λ₂ ≥ λ_z , one type of calculated value that quantifies the estimates' reliability is:

Other types of estimate reliability measures may also be utilised in connection with the inventive method.

Parameter tuning: Empirically, it has been determined that useful values for σ, may be in the range of 0.5 to 3.0 and for <τ₂ may be in the range of 1.5 to 6.0. Too large coefficients may blur the dip and azimuth estimates too much, while too small values may make the estimates too noisy. The parameter values of σ, = 0.5 and for σ₂ = 3.0 have been empirically determined to yield good results in some applications for several data sets, including the data set shown in Figure 8. Other applications may require different values. 2D Structural Characteristics Determination Methods

The dip and azimuth estimation approach described thus far for three- dimensional data may also be specialised for two-dimensional signals. It can be shown that the two-dimensional specialisation of the above procedure simplifies to the approach presented below.

The gradient vector is estimated similarly as in the three-dimensional case, i.e. by orientation selective filtering, such as by DoG filtering. The angle of this vector will have a good correspondence with the angle of plunge in the image at strong edges. However, in regions with low energy (like in between seismic reflectors), the amplitude of Vx will be low and the angle estimate less reliable. Hence, some smoothing of the gradient angle, weighted by the gradient magnitude, is desired. Note that this smoothing corresponds to the windowing in 3D.

However, adjacent vectors at angles separated by approximately π (i.e. pointing in opposite directions) with similar magnitudes will tend to cancel each other with such a smoothing. To alleviate this problem, the vectors are transformed to a space where vectors separated by an angle of π are identical. This may be done by representing Vx by complex numbers, with e.g. the horizontal component corresponding to the real part and the vertical to the imaginary. Then the arguments of the squared complex numbers are smoothed, weighted by their amplitudes. To illustrate the effect of this process, consider the two complex numbers

0 5 z, = 2e and

Squaring yields

2 ^ 2 /2 0.50 z l = 2 e

and

_z ² _{2 = 2} ² e^JΑ0-^5+π) = 2² e^{j +2π)} = 2² e^Λ20-⁵⁰⁾ = z² _i .

Hence,

and z²2 are equal.

Finally, the estimate of the angle of plunge is obtained as the angle of the square root of the elements from the process above. Note that with seismic data, it is meaningless to consider angles separated by π as different (e.g. a horizontal reflector has no direction left-to-right or right-to-left).

As in 3D, several possible realizations of the smoothing and derivative filters exist. Sets of particularly good candidates are Gaussian and DoG filters. As we have seen, these filters are separable, lending themselves to particularly efficient implementations.

With these choices, the 2D angle of plunge estimator may be summarized as:

1. Orientation selective filtering, typically gradient estimation by filtering with a vertical and a horizontal partially differentiated Gaussian filter. 2. Combine these two partial derivatives into a vector.

3. Transform the resulting vector elements by squaring the complex representation.

4. Weighted Gaussian smoothing (accomplished by smoothing the complex image).

5. Square root of the complex result.

6. The angles of these complex numbers correspond to the local angle of plunge.

This approach corresponds to the three-dimensional approach in that steps 1 and 2 correspond to the gradient estimation process and steps 3, 4, 5, and 6 correspond to the windowed covariance estimation and principal component selection processes.

Other methods for determining similar structural characteristics are known and may be used in connection with the present method, including those described in T.P.H. Steeghs, Local Power Spectra and Seismic Interpretation, Ph.D. dissertation, Delft University of Technology, Delft, the Netherlands, Sept. 1997 and M.S. Bahorich and S. Farmer, "3-D seismic for faults and stratigraphic features: The coherence cube", The Leading Edge, vol. 14, pp. 1053-1058, Oct. 1995, both incorporated herein by reference.

Fault Structural Characteristics Determination Methods

The types of structural characteristics described thus far have been associated with seismic reflectors, typically the interface surfaces between different geologic sedimentary layers. If the subsurface geologic feature to be identified by the present method is a fault, these types of structural characteristics can be used to calculate structural characteristics associated with subsurface fault features. These fault structural characteristics determination methods are based on the dip and azimuth estimation technique described above.

If a fault intersects a strong seismic reflector, the reflector typically becomes discontinuous and is characterised by an abrupt change in the seismic signal. Derivative filters will enhance abrupt changes, while suppressing smooth regions. However, the reflection layers are also intrinsically abrupt signal changes, and will also be enhanced. Hence it is necessary to perform the differentiation of the seismic signal along the reflection layers. In general, it can not be assumed that the layers are horizontal. The dip and azimuth of the layers therefore need to be estimated first.

A separable differentiation may be performed, where the components of the derivative are combined according to the dip and azimuth estimates. For example, if the angle of plunge of the seismic at a particular position in a 2D section is θ and the vertical and horizontal partial derivatives are

dx(n_χ ,n₂) dn.

and

dx(n ,n₂) dn.

respectively, a structural characteristic that quantifies the likelihood that a fault feature passes through a particular cell is dx(n_λ ,n₂) dx(n_x ,n₂) y = cos(#) + sin(0) dn. dn.

This type of structural characteristic is hereinafter referred to as a "fault characteristic".

In 3D, the local dip and an azimuth value represent a plane. Consequently, the fault characteristic may be interpreted as the projection of the vector with the three partial derivatives,

dx(n ,n₂,n ) dn_x dx(n ,n₂,n^)

Vx = dn₂ dx(n_x,n₂,n₃) dn.

onto the dip/azimuth plane. In smooth regions, the projected vector will have a small magnitude, whereas it will be larger near abrupt signal changes (i.e. near a fault feature).

Fault Cell Grouping/Thinning

The fault structural characteristics values will have high values when it appears that a fault passes through a particular cell and low values when it appears that a fault does not pass through a cell. In its most simplistic form, fault cell grouping may comprise grouping or associating cells having relatively high fault structural characteristic values with nearby cells also having relatively high fault structural characteristic values. This step is shown generally in Figure 1 as the Group Cells step 20. Three primary problems, however, may complicate the fault feature identification process: 1. The faults may have discontinuous structural characteristic values. That is, the fault structural characteristic values may occur as a series of blobs instead of as a continuous region. This appears to be primarily due to the lack of appropriate seismic signals in areas between seismic reflectors;

2. The noise behind the signal may distort the determined seismic characteristics; and

3. The seismic response normal to the fault may be too weak for easy fault feature identification.

An experienced human interpreter may be able to pick the centre-line of the determined fault structural characteristic values, but for a computer to make a robust pick, a "thinning" or reduction procedure should typically be applied.

Smoothing the fault attribute will reduce problems associated with noise and lack of continuity. On the other hand, unless the smoothing filter is well aligned with the dip and azimuth of the fault, it would increase the problem with blurred responses. However, by smoothing, it can be expected that the peak of the smoothed attribute will be close to the center of the fault response.

To utilize this, the neighborhood of the fault is traced perpendicular or nearly perpendicular to the fault. Assume that we have dip and azimuth estimates in the vicinity of a fault that are perpendicular or nearly perpendicular to the fault and "flow lines" are lines following plunge angles or local dip and azimuth estimates along its path. For 3D attributes, a new dip and azimuth estimate is therefore necessary, estimating the orientation perpendicular to the fault. One way to make this estimate is to use the dip and azimuth of the gradient vector after projection onto the orientation plane. This vector will generally be semi- orthogonal to the faults. In order to have reliable dip and azimuth estimates, these estimates should also be smoothed.

Consider, for instance, a trace of fault characteristic values extending n samples to each side of the current position. During the thinning process, if the value at the current position is the largest value in the group, it is retained and otherwise it is discarded. As a result, the fault characteristic response is thinned and only cells with the largest amplitudes are retained. This type of process may also be used to "thin" or "depopulate" other types of structural characteristics discussed below.

In summary, the process of computing a 3D set of fault structural characteristics using this embodiment of the present method may be described as follows:

1. Read a 3D seismic sub-volume from a cellular seismic database.

2. At each cell, estimate the local dip and azimuth of the seismic reflectors, corresponding to a local plane. This plane is represented by its normal vector. The dip and azimuth estimate includes the selection of two Gaussian filter parameters.

3. Compute the gradient (more precisely the derivative of Gaussian or "DoG") in the seismic volume. The width of the Gaussian is a parameter of this step. To have a sharp response, this parameter should be quite low. 4. Project this gradient vector onto the local plane. This projection is now an initial fault structural characteristic, where the projection removes the derivative component caused by the reflector. The only components left are the components due to faults and other discontinuities in the signal.

5. Smooth the dip and azimuth of the projected derivative as described. It is not very critical to have a dip and azimuth estimate that provides a very accurate fault plane orientation. The use of this parameter is for tracing a flow line across the fault, and deviations from the "normal" or "perpendicular" direction can be tolerated. It is, however, important that the dip and azimuth is consistently non-parallel to the fault. Consequently, the smoothing filter should be quite large.

6. Compute a refined fault structural characteristic by Gaussian smoothing of the magnitude of the result from step 4. The smoothing removes noise and "connects" discontinuous fault responses. The size of this smoothing filter, σ₅ , is probably the most critical parameter of this approach and the parameter can be regarded a scale-space parameter. Consequently, a small σ₅ corresponds to a fine scale, and subtle fault responses and fine detail is enhanced. However, with a small σ₅ , the noise suppression is low, and a lot of minor "artificial" faults are identified. A large σ₅ , on the other hand, suppresses more noise and only enhances the faults with major discontinuity responses.

7. Compute the displacement corresponding to the dip and azimuth above. 8. Trace flow lines beginning at each cell. If the beginning cell of the flow line has the maximum structural characteristic value, this value is retained. Otherwise, it is marked as blank.

9. Write the resulting sub-volume to the seismic database.

Other types of subsurface geologic features may also be identified by the inventive method, such as seismic reflectors.

Common Reflector Cell Selection/Thinning

In its most simplistic form, cells associated with a common reflector may be selected merely by selecting cells having certain seismic reflector inclination values or estimate reliability measures. Groups of nearby cells in a time to depth converted seismic data set having zero or nearly-zero seismic reflector inclinations (i.e. lying on a horizontal or nearly horizontal surface) may, for instance, be selected. These types of seismic reflectors may comprise hydrocarbon/water interfaces and may be direct predictors of hydrocarbon deposits. These types of "flat spots" may be readily extracted from seismic data using the inventive method under certain combinations of subsurface geology and fluid fill conditions.

Nearby cells having relatively high estimate reliability measures can also be selected. Cells that have relatively high estimate reliability measures are typically located on the edges of locally dominant seismic reflectors.

In some cases, the selected seismic reflector cells will need to be thinned. An experienced human interpreter may be able to pick the "center of gravity" of a group of selected cells, but for a computer to make a robust pick, a "thinning" or reduction procedure may need to be applied. The reflector dip and azimuth estimates may be used, for instance, to calculate a vector largest value in the group, it is retained and otherwise it is discarded. Alternatively, the center-most cell in the group may be retained. This type of process may also be used to "thin" or "depopulate" selected termination cells, as discussed below.

Seismic Reflector Extraction by Flow Surface Generation

As described above, in some cases seismic reflector positions may be extracted from seismic data directly on the basis of particular seismic reflector inclination values or estimate reliability measures. In other cases, the Group Cells step 20 will involve the identification of flow surfaces. Starting from a seed position, a "flow surface" is the surface through that position having the angle of the local dip and azimuth estimates (perpendicular to the gradient angle) in its path. Flow surfaces may be 2D surfaces in a 3D volume, 1D lines in a 3D volume, or 1D lines on a 2D cross section. In the case of lines, the flow surface may be denoted a flow line. Flow lines on a cross section may be generated by using a 2D angle of plunge estimate, or by projecting the dip and azimuth onto the section. A flow line is a straight-forward specialization of a flow surface and will not be treated separately. A flow surface may be represented as a time (or depth) interpretation, i.e. as an interpretation grid, or as a set of surface patches.

The inventive method is particularly suited for extracting subtle seismic reflectors that are difficult or impossible to extract using conventional methods. Because the inventive method calculates a local seismic reflector inclination for each of the cells examined, a flow surface associated with each cell in the database may be generated, regardless of seed cell and the relative magnitude of the seismic response from that cell. Many prior art methods will fail to identify any surface whatsoever if the magnitude of the contrast in the seismic data is small. Seed Cell Selection

In some cases, seed cells may be manually input because they represent known points on a subsurface geologic feature, such as known reservoir boundary intersection points. Cells associated with beginning and end points of an oil column could be determined from a well log and input in the present method as seed cells to produce a detailed image of the reservoir cap rock and the oil/water interface, for instance.

In other cases, it may be desirable to identify a large number of seismic reflectors within a data set. Two seed cell layout schemes have proven useful in this type of application of the inventive method.

1. Seed cells are located on a set of traces with several seeds per trace.

2. Having some scanning order of the samples, additional seed cells are located at positions where the distance to the closest previously identified geologic feature is above some threshold.

It should be remembered that the selection of one or more seed cells merely allows the flow surface cell selection process a place to start. The seed cells can therefore also be thought of as beginning locations or areas in the cell grouping process. Momentum

The concept of "momentum" can also be used during the flow surface generation process. When using this momentum component, the local dip and azimuth estimates are smoothed over neighboring estimates on the same surface. Hence, if the angular estimates would make a flow surface bend using the previous approach, now the flow surface is forced to continue straighter.

This flow surface generation process may be terminated when the flow surface meets a previously generated flow surface, when an abrupt change in flow direction is indicated by the structural characteristics values of the latest grouped cell, etc.

Termination Identification

The points at which the geologic strata terminate or become relatively close to each other may also be identified. As discussed above, identifying termination points is an important activity associated with the location and evaluation of stratigraphic traps. If the dip and azimuth of a particular geologic stratum is traced and we detect where this stratum terminates against another stratum, an onlap, downlap, toplap, or erosional truncation is likely and these features are associated with potential stratigraphic hydrocarbon traps.

If two flow surfaces intersect, this may correspond to a geologic strata termination interface. However, due to interference in the seismic signal and due to the smoothing in the dip/azimuth estimate, flow surfaces tend to converge rather than intersect at the termination interfaces. One way to alleviate this is to assign a width or window to each flow surface. If a flow surface intersects an area covered by the window around another flow surface, the intersection is interpreted as a termination. Terminations may be detected practically by producing a binary volume with the flow surfaces while the flow surfaces are generated. If a flow surface approaches a cell marked as visited by an earlier flow surface, this cell is marked as a geologic stratum termination cell.

It is also noted that terminations have proven to be good indicators of other geologic features and seismic signal characteristics than terminations, e.g. faults. The identified termination cells may therefore be used for purposes other than solely for stratigraphic interpretation.

Stratigraphic Fades Identification

The structural characteristics values may also be used in the identification of the stratigraphic fades, i.e. the internal structure or texture of the strata in the subsurface area. Several alternative processes may be used to compute and represent this texture information. Some will yield an attribute value for a confined sub-volume, while others will yield one attribute value for each seismic sample. The dip/azimuth estimation technique and the flow surface descriptions discussed above provide powerful primitives for generating fades texture attributes.

Attributes can be created, for instance, that highlight chaotic regions of the seismic signal. One such method uses the smoothed magnitude of the gradient in the dip and azimuth estimates, normalized with respect to the smoothed magnitude of the seismic. If the dip/azimuth estimate is reliable, this measure will yield a high value. Chaotic regions will typically not yield reliable estimates. An alternative is the reliability measure of the dip and azimuth estimate discussed above. Similar attributes can be calculated by measuring the regularity of the flow lines. For each sample position, two flow lines can be generated starting from a distance above and below this sample position. The attribute value at this sample position is then set to be the standard deviation of the change in distance between the flow lines. In regions with chaotic signals, the distance between the flow lines will change a great deal, and the standard deviation will be high.

Attributes can also be developed that highlight other types of statigraphic fades patterns, such as divergent, straight, and wavy parallel patterns. One such type of attribute can be calculated by initiating flow lines of a certain length above and below each sample of the seismic data. Then the distance between the flow lines' respective left and right end points, d, and d_r can be calculated. Attribute values can then be calculated as \d,

.

Divergent signals will lead to flow lines with consistent differences in the differences in the distance between the flow line end points, whereas parallel signal configurations will lead to fairly constant distances. A similar method computes the attribute value as the mean value of the flow lines' separation gradient normalized by the standard deviation of the separation gradient. Divergent flow lines will have a somewhat constant gradient with nonzero mean and small standard deviation.

It is easy to see how these techniques can yield a good separation of the parallel from the divergent signal characteristics. Used in combination with each other or with the chaos attributes its is then easy to discriminate between divergent, parallel and chaotic stratigraphic fades.

Represent/Display Results

The cells grouped in the Group Cells step 20 collectively provide a granular representation of the identified geologic feature. While this type of representation will be sufficient for many purposes, the identified subsurface geologic feature may show undesired characteristics, such as holes or inaccurately grouped cells, and several different subsurface features may be connected into a single composite feature in the data volume. This latter problem is particularly likely if the feature being identified is a fault and the seismic data being interpreted has two or more faults that cross each other.

These issues may be addressed in the Represent/Display Results step 22. In this step, the grouped cells may be clustered or subdivided and the grouped cells or these clustered portions of the grouped cells may be represented as vectorized surfaces in three dimensions, such as a plane:

an₁ + bn₂ + cn₃ = ,

where n_1t n_2t and n₃ represent the three axis coordinates, or a higher order polynomial surface, such as

an² ₁ + bn₁ + cn₂ + dn³ ₃ + en₃ = f.

If the surface shape is too complex for such a model, each surface may be subdivided into multiple connected patches of polynomial models. The parameters of this surface representation may be estimated using least squared error approximation or other parameter estimation techniques known to those skilled in the art.

If less than all of the grouped cells are to be represented, the cells must be suitably clustered, i.e. which samples correspond to which surfaces must be determined. In some cases, the clustering and vectorizing processes may be performed as a sequential process. In other cases, better results may be obtained using an iterative method where the results of the vectorizing process are used to refine the results of the clustering process.

Several alternative methods for clustering the grouped cells may be used in the inventive method. In one method, the grouped cells are subdivided and connected components are extracted from these subdivided units. Connected components are sets of cells that are connected within a geometrical neighborhood. A parametric surface representation of these connected components is made and the degree of fit of this parametric surface representation with respect to all of the samples in the larger volume is measured. Next, the cells with best degree of fit are used to recalculate the surface parameters, then new degrees of fit are computed, etc. The last part of this process is iteratively repeated until a sufficient degree of fit is obtained.

Other clustering schemes may also be employed including letting an interpreter select (portions) of the connected components to start from; estimating parameters, and iteratively fitting computations and reestimations as above. Detailed information on general clustering techniques may be found in Pattern Recognition: Statistical, Structural and Neural Approaches by Robert Schalkoff, John Wiley, 1992, incorporated herein by reference.

The grouped, clustered, and/or vectorized cells may then be suitably displayed as shown in Figures 4 through 9. Figures 4 and 5 display the types of subsurface fault features that may be identified by the present method. In Figure 4, section 50 is a sectional view through a seismic data cube with fault features 52 showing the types of faults that may be identified by the present method. Figure 5 shows a perspective view of the entire seismic data cube 60 as well as perspective view of faults 62 identified by the present method. In the processing used to create Figures 4 and 5, the parameters used were σ, = 3.0 , σ₂ = 3.0 , σ₃ = 1.0 , σ₄ = 6.0 and σ₅ = 5.0.

In this example, the cells were grouped into blocks of size 40x40x40 when producing the vectorized fault representations shown in Figures 4 and 5 and the parametric surface representation used was:

p_xn i + p n i + ₃π, + p n₂ + p_sn ₃ + p₆n^~ + p_ηn₃ = 1 where the p/s are the surface parameters. These parameters were estimated by least squares approximation.

Figures 6 and 7 show the types of seismic reflector features that may be identified by the present method. In Figure 6, section 70 is a sectional view through a seismic data cube (a collection of seismic data obtained from a particular subsurface area). First seismic reflector feature 72 and second seismic reflector feature 74 represent a pair of seismic reflectors of the type that may be identified by the present method. Derrick 75 shows a hypothetical well location and first wellbore location 76 and second wellbore location 78 represent a pair of locations within the well that could have acted as seed points for the location of first seismic reflector feature 72 and second seismic reflector feature 74. Figure 7 shows the types of three-dimensional seismic reflector features that may be identified by the present method. The first reflector surface 80 and second reflector surface 82 shown in Figure 7 may represent a more complete 3D view of the seismic reflectors shown in cross section in Figure 6.

The type of seismic reflector termination information that may be determined using the inventive method is shown in the 2D cross section in Figure 8. In the processing used to create Figure 8, flow lines were initiated at every row in every 10th column, using flow lines having a width of three cell or sample positions. If a flow line intersected the window around another flow line which was initiated at least five rows away on the same column, a termination density map was updated by adding one to that cell. The momentum concept was also used, with the last 31 angular estimates along the flow line being averaged.

Figure 9 is a cross-sectional view through a facies texture attribute cube prepared in accordance with the inventive method that highlights chaotic signal patterns as discussed above. Section 90 identifies a particular subsurface area 92 where the seismic data signal is highly chaotic.

The present method is particularly suited for autonomous or semi- autonomous seismic data interpretation, where steps 18, 20, 22, and 24 are performed with little or no manual operator input.

It will be readily understood that the steps and processes associated with the disclosed embodiment of the present method are capable of a wide variety of alternative implementation methods. The entire seismic data cube may be examined or only a particular area within the cube. Many of the steps may be performed iteratively. 2D methods may be used to identify 2D portions of the subsurface geologic feature and these 2D portions may then be aggregated and interpolated to produce a 3D representation of the subsurface geologic feature. The present method is in no way limited or restricted to the particular order of steps described in the preferred embodiment above.

Claims

CLAIMS:

1. A method of identifying a subsurface geologic feature using seismic data obtained from a subsurface area, said method comprising the steps of:

subdividing said subsurface area into a plurality of discrete cells; each cell representing a different region within said subsurface area;

associating said cells with portions of said seismic data that image subsurface regions represented by said cells;

determining structural characteristics for a plurality of said cells using said seismic data; and

grouping together nearby said cells based on said structural characteristics to identify said subsurface geologic feature.

2. A method according to Claim 1 , wherein said subsurface geologic feature comprises a fault.

3. A method according to Claim 1 , wherein said subsurface geologic feature comprises a seismic reflector.

4. A method according to Claim 1 , wherein said subsurface geologic feature comprises a hydrocarbon reservoir boundary.

5. A method according to Claim 4, wherein said hydrocarbon reservoir boundary comprises a hydrocarbon/water interface.

6. A method according to Claim 1 , wherein said grouping step comprises selecting nearby cells having essentially horizontal local seismic reflector inclinations.

7. A method according to Claim 1 , wherein said determined structural characteristics comprise an estimate reliability measure.

8. A method according to Claim 7, wherein said grouping step comprises grouping together nearby cells having relatively high estimate reliability measures.

9. A method according to Claim 7, wherein said grouping step comprises grouping together nearby cells having relatively low estimate reliability measures.

10. A method according to Claim 1 , wherein said grouping step comprises the step of selecting a seed cell from said plurality of cells, said seed cell representing a point on a first seismic reflector.

11. A method according to Claim 10, wherein said grouping step further comprises the step of selecting additional cells at increasing distances from said seed cell using said local seismic reflector inclinations, said additional cells representing additional points on said first seismic reflector.

12. A method according to Claim 11 , wherein said selected cells when aggregated with neighboring selected cells produce localized seismic reflector areas generally in alignment with said local seismic reflector inclination estimates of said cells through which said seismic reflector passes.

13. A method according to Claim 10, further including the step of extracting a flow surface using said seed cell and said structural characteristics.

14. A method according to Claim 13, wherein said flow surface is extracted using momentum.

15. A method according to Claim 13, wherein a plurality of said flow surfaces are extracted.

16. A method according to Claim 15, wherein cells where said extracted flow surfaces become relatively close together are identified as seismic reflector terminations.

17. A method according to Claim 1 , wherein said seismic data has been time domain to depth domain converted prior to said determining step.

18. A method according to Claim 1 , wherein said seismic data has been decomposed and reconstructed prior to said determining step.

19. A method according to Claim 1 , wherein said cells are thinned in said grouping step.

20. A method according to Claim 1, wherein said grouped together cells are subdivided using a clustering process.

21. A method according to Claim 1 , wherein a plurality of said grouped cells are represented as a vectorized surface.

22. A method according to Claim 1 , further comprising the step of displaying said grouped cells.

23. A method according to Claim 1 , wherein said seismic data associated with said cells are associated with a single subsurface imaging point.

24. A method according to Claim 1, wherein said determining step comprises orientation selective filtering said seismic data.

25. A method according to Claim 24, wherein said orientation selective filtering comprises estimating gradients of local seismic reflectors.

26. A method according to Claim 25, wherein said gradient estimation is performed by filtering said seismic data with a derivative of Gaussian filter.

27. A method according to Claim 26, wherein said gradient estimates are smoothed.

28. A method according to Claim 27, wherein said smoothing is performed by filtering said gradient estimates with a Gaussian low-pass filter.

29. A method according to Claim 1 , wherein said structural characteristics comprise a fault characteristic that quantifies a likelihood that a fault passes through said cell.

30. A method according to Claim 29, wherein said grouping step comprises grouping together nearby cells having relatively large fault characteristic values.

31. A method according to Claim 1 , wherein said cells are three- dimensional volume elements.

32. A method according to Claim 31 , wherein said structural characteristics comprise dip and azimuth values.

33. A method according to Claim 32, wherein said dip and azimuth values are determined by principal component analysis of smoothed gradient estimates.

34. A method according to Claim 1 , wherein said cells are two-dimensional area elements.

35. A method according to Claim 34, wherein said structural characteristics comprise plunge angles.

36. A method according to Claim 1, wherein said structural characteristics comprise stratigraphic facies attributes.

37. A general purpose computer programmed to implement the method of any one of Claims 1 to 36.

38. A computer program on a carrier containing instructions for causing a general purpose computer to implement the method of any one of claims 1 to 36.