Method and apparatus for seismic feature extraction

• US 20060122780A1
• Filed: 11/10/2003
• Published: 06/08/2006
• Est. Priority Date: 11/09/2002
• Status: Abandoned Application

• Current Assignee: Geoenergy Incorporated
• Original Assignee: Geoenergy Incorporated
• Inventors: Vassiliou, Anthony, Cohen, Israel, Coult, Nicholas
• Application Number:
  US10/534,149
• Publication Number:
  US 20060122780A1
• Time in Patent Office: Days
• US Class Current: 702/14
• CPC Class Codes: G01V 1/288 Event detection in seismic signals, e.g. microseismics G01V1/36 takes precedence
• CPC Class Codes: G01V 1/30 Analysis G01V1/50 takes precedence
• CPC Class Codes: G01V 1/301 for determining seismic cross-sections or geostructures
• CPC Class Codes: G01V 1/364 Seismic filtering G01V1/37 takes precedence

Abstract

A method and apparatus for seismic image processing is disclosed. A preferred embodiment aids in the identification of subterranean faults, which are significant in hydrocarbon exploration. The method includes steps of: a) reading a three dimensional seismic data volume; b) computing the three-dimensional orientation of the subsurface; c) subdividing the original volume into small data volumes that are rotated at a predetermined set of dips and azimuths related to those of the subsurface orientation; volumes formed in step c; e) performing a 3-D contrast enhancement operation in each of the small volumes; f) filtering the result of the contrast enhancement with selected 3-D filters at the predetermined set of dips and azimuths; g) skeletonizing the results of the filtering operation; h) separating the individual fault surfaces, and i) labelling the individual fault surfaces for further interpretation and exploration.

Claims

1. A method for seismic exploration, comprising actions of:

(a) obtaining a set of seismic data traces representing a three-dimensional volume of seismic data samples;

(b) dividing the three-dimensional volume into a plurality of smaller subvolumes;

(c) selecting a discrete set of dip values and azimuth values;

(d) dividing the three-dimensional volume into a plurality of parallelepipeds, each of the parallelpipeds being tilted by one of the selected dip values and rotated by one of the selected azimuth values;

(e) halving each parallilepiped to obtain two half-parallelpipeds, wherein the two half-parallelpipeds have contain an equal number of samples and there exists a one-to-one relationship between corresponding samples in the two half-parallelpipeds;

(f) enumerating the samples in each of the two half-parallelpipeds so as to obtain two vectors, such that corresponding samples in the two half-parallelpipeds have corresponding indices in the two vectors;

(g) calculating a three-dimensional edge detection measure from the two vectors;

(h) associating the computed edge detection measure to the parallilepiped center point to obtain a first subresult;

(i) applying a constrast enhancement measure to the first subresult to obtain a second subresult;

(j) filtering the second subresult to obtain a set of third subresults by convolving the second subresult with a directional filter kernel that is tilted and rotated in accordance with a set of dip and azimuth values that correspond to the selected dip values and azimuth values of the computational analysis parallilepiped associated with the second subresult;

(k) selecting a maximum filtered value from the set of third subresults;

(l) applying a three-dimensional skeletonization algorithm to the maximum filtered value to generate a skeleton representing a fault surface;

(m) executing action (l) with respect to the maximum filtered value for each of the plurality of parallelpipeds to obtain a plurality of distinct skeletons; and

(n) labeling each of the plurality of distinct skeletons as a separate geologic feature.

2. The method of claim 1, wherein the method comprises an additional action of:

(b1) computing, for at least one of the subvolumes, a local seismic data orientation and an associated statistical uncertainty, and further, wherein each of the plurality of parallelpipeds is associated with a subvolume from the plurality of subvolumes and each of the plurality of parallelpipeds has a top surface and a bottom surface that are parallel to the local seismic data orientation for the subvolume associated with that parallelpiped.

3. The method of claim 2, where action (b1) further includes:

(b1a) applying a smoothing function to a subvolume;

(b1b) computing, in three dimensions, a local orientation of seismic data within the subvolume;

(b1c) averaging, over a surrounding three-dimensional space, the local orientation of seismic data within the subvolume to obtain a stable orientation estimate; and

(b1d) computing a measure of statistical uncertainty of the orientation estimate.

4. The method of claim 2, where action (b1) further includes:

(b1a) transforming the seismic data of a subvolume so as to map Cartesian coordinates of the subvolume into complex values in a complex plane;

(b1b) computing a three-dimensional gradient of the subvolume;

(b1c) convolving the three-dimensional gradient with a three-dimensional Gaussian filter to obtain a convolution result;

(b1d) extracting a local orientation estimate of the subvolume as the argument (in the sense of a complex number) of the convolution result; and

(b1e) computing an uncertainty measure (ρ

) of the local orientation estimate as

$\mathrm{\rho <br><br>}=\frac{{\mathrm{\lambda <br><br>}}_1-{\mathrm{\lambda <br><br>}}_0}{{\mathrm{\lambda <br><br>}}_1+{\mathrm{\lambda <br><br>}}_0},$

where λ

₀ and λ

₁ are eigenvalues of a scatter matrix for the subvolume'"'"'s local spectrum.

5. The method of claim 1, wherein the selected dip values are selected from an interval extending from about −

45 degrees to about 45 degrees and with a dip increment selected from an interval extending from about 5 degree to about 10 degrees.

6. The method of claim 1, wherein the selected azimuth values are selected from an interval extending from about 0 degrees to about 180 degrees and wherein the selected azimuth values are selected using an azimuth increment selected from an interval extending from about 5 degrees to about 30 degrees.

7. The method of claim 1, wherein at least one of the parallelpipeds is tilted by a selected dip value and rotated by a selected azimuth angle to form a parallelpiped having integer dimensions of the form L₁×

(2L₂+1)×

N, with N representing a number of samples measured with respect to a vertical dimension of the parallelpiped, with the selected dip value measured with respect to the vertical direction, and the azimuth measured with respect to a direction of north in the seismic data.

8. The method of claim 7, wherein each parallelpiped has a minimum dimension L₁ of 6, a minimum dimension 2L₂₊₁ of 5 samples and a minimum vertical dimension N of 41 samples.

9. The method of claim 7, wherein the set of selected azimuth values is exclusively dependent on L₁.

10. The method of claim 9, where the set of selected azimuth values has a cardinality of at least 8.

11. The method of claim 1, wherein the edge detection measure for a given parallelpiped having a central analysis point λ

is at least a function of

∥

v_1,λ

(γ

,φ

)−

v_2,λ

(γ

,φ

)∥

_p
and
∥

v_1,λ

(γ

,φ

)∥

_p+∥

v_2,λ

(γ

,φ

)∥

_p

where v_1,λ

and v_2,λ

are the two vectors generated from halving the parallelpiped, γ and

φ

are the parallelpiped'"'"'s dip, and azimuth, respectively, and p denotes a choice of vector space norm for a vector space in which v₁ and v₂ are defined.

12. The method of claim 11, wherein the three dimensional edge detection measure is a function of

$\frac{{}_{}}{}$

v 1 , λ

⁡

( γ

, ϕ

) - v 2 , λ

⁡

( γ

, ϕ

) 

p 

v 1 , λ

⁡

( γ

, ϕ

) 

p + 

v 2 , λ

⁡

( γ

, ϕ

) 

p .

13. The method of claim 12, wherein values of

∥

v_1,λ

(γ

,φ

)−

v_2,λ

(γ

,φ

)∥

_p
and
∥

v_1,λ

(γ

,φ

)∥

_p+∥

v_2,λ

(γ

,φ

)∥

_p

are computed only once per each pair of vectors v_1,λ

and v_2,λ

and dip/azimuth combination, are subsequently stored in memory so as to memoize the computed values such that the computed values need not be recomputed for that pair of vectors and dip/azimuth combination.

14. The method of claim 13, wherein the memoized computed values from a window of preceding edge detection measure computations are utilized, through dynamic programming, to compute subsequent edge detection measures.

15. The method of claim 1, wherein the contrast enhancement measure is a function of a filter that is applied by convolving the first subresult with a rotated form of a “

Mexican hat function”

of the form

f(n)=C(1−

n²)e<sup>−

n2/2</sup>,

(where e is the base of the natural logarithm function) for values of n taken from an interval extending from about −

4.5 to 4.5, the filter containing an odd number of filter coefficients and having its main direction perpendicular to slabs of the first subresult'"'"'s corresponding parallelpiped, wherein convolving the first subresult with the rotated form of the “

Mexican hat function”

results in a “

Mexican hat function convolution result.”

16. The method of claim 15, wherein the contrast-enhanced three-dimensional edge detection measure is computed as the maximum of the “

Mexican hat function convolution result” and

zero.

17. The method of claim 1, where portions of fault surfaces that are approximately aligned with a parallelpiped are extracted.

18. The method of claim 17, wherein the directional filter kernel represents a three-dimensional finite impulse response filter having a set of coefficients characterized by a primary axis, a secondary axis, and a tertiary axis, wherein the secondary and tertiary axes are significantly shorter in length than the primary axis.

19. The method of claim 18, wherein the three-dimensional filter is rotated by an angle γ

+α

with respect to a time axis and rotating it by an angle φ

around an in-line axis, wherein α

represents a being a tolerance of dip increment.

20. The method of claim 19, wherein the three-dimensional filter has a number of coefficients in the primary direction that is at least equal to a number of samples in a corresponding direction of the parallilepiped.

21. The method of claim 20, wherein the three-dimensional filter is a long three-dimensional Hanning window with very small number of coefficients along the secondary and tertiary axes in relation to the number of coefficients along the primary axis.

22. The method of claim 21, wherein obtaining the set of third subresults includes applying, to a result of the three-dimensional filter, a threshold function defined by a user-determined threshold.

23. The method of claim 22, further comprising filtering the set of third subresults to produce a three-dimensional volume output in which portions of fault surfaces parallel to the selected dip and azimuth values are distinguished.

24. The method of claim 23, further comprising selecting the maximum directionally filtered output value over the filtered set of third subresults.

25. The method of claim 24, wherein the three-dimensional skeletonization algorithm is applied to the maximum directionally filtered output value.

26. The method of claim 24, wherein action (m) further includes actions of:

(m1) computing a local orientation volume for the set of third results, the computation of the local orientation volume being accomplished by computing local gradient vectors at each volume point using a 2-point stencil in each direction and finding a principal component of the local gradient vectors using singular value decomposition;

(m2) finding the location x of a largest-valued point in an output of action (I) that has not already been marked as masked point and ending execution of action(m) if largest-valued point'"'"'s value is smaller than a threshold fault value;

(m3) selecting points that are less than a predetermined distance of N₁ from x in a horizontal plane containing x, with N₁ and that have a local orientation vector that is within a preselected orientation difference threshold;

(m4) selecting points that lie in a plane perpendicular to the local orientation vector at the point x, that are less than a predetermined distance of N₂ from x, and that have a local orientation vector that is within a preselected orientation difference threshold;

(m5) repeating actions (m3) and (m4), considering one iteration for all of the points selected in actions (m3) and (m4) until the total number of iterations exceeds a preselected number of iterations;

(m6) marking output volume points around points of maximum fault volume value as masked points in response to a determination that a total number of points selected is less than a preselected number of minimum points;

(m7) skeletonizing fault surfaces composed of the selected points from actions (m3) and (m4) in response to a determination that the total number of selected points is greater than a preselected number of points;

(m8) zeroing selected fault output volume points in response to a determination that a number of selected points is smaller than a preselected number of points;

(m9) in response to a determination that a number of points selected during a immediately preceding iteration of skeletonization is greater than a preselected number of points, labelling the skeletonized fault surfaces with labels that are unique to each fault surface.

27. The method of claim 26, wherein skeletonization includes:

selecting a point “

γ

”

in accordance with either action (m3) or action (m4);

searching along a line that passes through γ and

is parallel to the local orientation vector and finding a location of a maximum fault output volume along that line that is within N₂ grid points of γ

,

marking in a new output volume the location of the maximum fault output volume value with a non-zero value; and

removing γ

from the set of points selected in actions (m3) and (m4)

28. The method of claim 25, wherein the skeletonized filtered volume is stored in magnetic media for subsequent use.

29. The method of claim 29, wherein descriptions of individual fault surfaces are stored in magnetic media for subsequent use.

30. An apparatus for processing and analyzing seismic trace data, comprising means for:

(a) obtaining a set of seismic data traces representing a three-dimensional volume of seismic data samples;

(b) dividing the three-dimensional volume into a plurality of smaller subvolumes;

(c) selecting a discrete set of dip values and azimuth values;

(d) dividing the three-dimensional volume into a plurality of parallelipipeds, each of the parallelpipeds being tilted by one of the selected dip values and rotated by one of the selected azimuth values;

(e) halving each parallilepiped to obtain two half-parallelpipeds, wherein the two half-parallelpipeds have contain an equal number of samples and there exists a one-to-one relationship between corresponding samples in the two half-parallelpipeds;

(f) enumerating the samples in each of the two half-parallelpipeds so as to obtain two vectors, such that corresponding samples in the two half-parallelpipeds have corresponding indices in the two vectors;

(g) calculating a three-dimensional edge detection measure from the two vectors;

(h) associating the computed edge detection measure to the parallilepiped center point to obtain a first subresult;

(i) applying a constrast enhancement measure to the first subresult to obtain a second subresult;

(j) filtering the second subresult to obtain a set of third subresults by convolving the second subresult with a directional filter kernel that is tilted and rotated in accordance with a set of dip and azimuth values that correspond to the selected dip values and azimuth values of the computational analysis parallilepiped associated with the second subresult;

(k) selecting a maximum filtered value from the set of third subresults;

(l) applying a three-dimensional skeletonization algorithm to the maximum filtered value to generate a skeleton representing a fault surface;

(m) executing action (I) with respect to the maximum filtered value for each of the plurality of parallelpipeds to obtain a plurality of distinct skeletons; and

(n) labeling each of the plurality of distinct skeletons as a separate geologic feature.