# Patent

## Title
Detection of features in seismic images

## Abstract

A computer implemented method for the detection of features such as faults or channels in seismic images. First, edges are detected in a smoothed seismic image (106). To detect a fault line, an image intensity of the edges is projected in multiple spatial directions, for example by performing a Radon transform (118). The directions of maximum intensity are used to define a fault line (124c). To detect channels, smooth curves are identified within the detected edges (810). Sets of parallel smooth curves (812a) are then identified and used to define channels (812).

## Background
This application is the National Stage under 35 U.S.C. 371 of PCT/US2008/081917 that was published as WO 2009/082545 and was filed on 31 Oct. 2008, which claims the benefit of U.S. Provisional Application No. 61/008,758, filed 21 Dec. 2007, which is incorporated herein by reference, in its entirety, for all purposes.


BACKGROUND OF THE INVENTION
This invention relates generally to oil and gas production, and in particular to the processing of seismic data to facilitate oil and gas production.
Existing interpretation of seismic data faces a number of significant challenges. First, undiscovered oil and gas reservoirs are often contained in subtle traps at very low depth or in stratigraphic traps that are difficult to detect in a single seismic volume. Second, the amount of seismic data and attributes that can now be obtained for each seismic volume typically far exceeds the interpretation capability of human interpreters. As a result, a number of methods have been developed for automatically processing seismic data to facilitate the interpretation of the seismic data.
For example, there are a number of existing methods for applying pattern recognition and data mining techniques to seismic data. However, many of the existing methods have been limited to pre-processing seismic images to remove noise and/or utilize horizon tracking using local classification methods. In other instances, the existing methods provide automatic identification of geological features by relying upon training neural networks and other similar classifiers. Such existing methods for identifying geological features in seismic data are indirect and lack transparency.



BRIEF DESCRIPTION OF THE DRAWINGS
FIGS. 1a and 1b are flow chart representations of an exemplary embodiment of a method of determining the presence of a fault in a seismic image.
FIG. 2 is a graphical illustration of an exemplary embodiment of a seismic image.
FIG. 3 is a graphical illustration of the seismic image of FIG. 2 after smoothing the seismic image.
FIG. 4 is a graphical illustration of the seismic image of FIG. 3 after detecting the points of high radius of curvature in the smoothed seismic image.
FIG. 5 is a graphical illustration of a fault line in the image of FIG. 4.
FIG. 6 is a graphical illustration of Radon transforms performed on the image of FIG. 4.
FIG. 7 is a graphical illustration of the seismic image of FIG. 2 including a fault line.
FIG. 8 is a flow chart representation of an exemplary embodiment of a method of determining the presence of channels in a seismic image.
FIG. 9 is a graphical illustration of an exemplary embodiment of a seismic image.
FIG. 10 is a graphical illustration of the seismic image of FIG. 9 after smoothing the seismic image.
FIG. 11 is a graphical illustration of the smoothed seismic image of FIG. 10 after processing the image to determine channels in the image.



DETAILED DESCRIPTION OF THE ILLUSTRATIVE EMBODIMENTS
Referring to FIGS. 1a, 1b, 2, 3, 4, 5, 6, and 7, an exemplary embodiment of a method 100 of detecting faults in seismic images is provided in which a seismic image 102a is input in 102. In an exemplary embodiment, the input seismic image 102a is then smoothed in 104 to generate a smoothed seismic image 104a. In an exemplary embodiment, the smoothed seismic image 104a is generated in 104 using conventional methodologies for smoothing multi-dimensional data such as, for example, methodologies for the filtering of random noise from multi-dimensional data.
In an exemplary embodiment, the smoothed seismic image 104a is then processed to determine edges within the smoothed seismic image in 106. In an exemplary embodiment, the determination of edges within the smoothed seismic image 104a is provided in 106 using conventional methodologies for determining edges in multi-dimensional data such as, for example, the Canny or Sobel edge detection methodologies. In an exemplary embodiment, the radius of curvature of the determined edges are then determined in 108. In an exemplary embodiment, the radius of curvature determined in 108 is the iso-contour radius of curvature of the smoothed seismic image 104a. The iso-contour radius of curvature of the determined edge refers to the radius of curvature of the edge determined in 106.
In an exemplary embodiment, the method 100 then determines the points of high radius of curvature in the smoothed seismic image by determining, for each edge, whether or not the radius of curvature of the edge is less than or equal to a predetermined reference value R(x,y) in 110. If the radius of curvature of the edge is not less than or equal to the predetermined reference value R(x,y), then the radius of curvature for the edge is rounded off to be equal to zero in 112. If the last edge having a non-zero radius of curvature has been processed by 110 and 112, then the determination of the points of high radius of curvature is ended in 114 and the method 100 then generates an image of the points of high radius of curvature 116a in 116.
In an exemplary embodiment, the method 100 then performs a windowed Radon transform for a range of angular directions on the image of the points of high radius of curvature 116a in 118, 120, and 122. In an exemplary embodiment, the Radon transform utilized in 118 may be a conventional Radon transform. Alternatively, other similar transforms such as, for example, a Hough transform, may also be used in 118. In an exemplary embodiment, the window size for the Radon, or equivalent, transform, in 118 may help to control errors in due to potential small misalignments and/or some level of curviness in the faults determined by the method 100. In an exemplary embodiment, the window size used in the Radon, or equivalent, transform, in 118 should be set to be equal to a small multiple of the horizon spacing of the seismic image 102a. 
In an exemplary embodiment, the method 100 then determines the local maxima 124a of the results of the Radon transform results generated in 118, 120, and 122 in 124. In an exemplary embodiment, the local maxima 124a of the Radon transform results is associated with a corresponding angular direction 124b and a fault line 124c within the seismic image 102a and the image of the points of high radius of curvature 116a. 
In an exemplary embodiment, the method 100 then maps the fault line 124c onto the seismic image 102a in 126.
In an exemplary embodiment, the teachings of the method 100 may be applied to seismic images having N-dimensions.
In an exemplary embodiment, the method 100 may be used to determine one or more fault lines.
In an exemplary embodiment, in the method 100, grouping algorithms may be used in combination with Radon and/or Hough transforms in order to enhance the determination of fault lines.
In an exemplary embodiment, in the method 100, contextual geological may be used such as, for example, proximity to other geological structures and constraints on orientation based upon rock physics, may be used to enhance the determination of fault lines.
In an exemplary embodiment, the method 100 may be used to determine linear, non-linear, or curvilinear fault lines. In an exemplary embodiment, the determination of curvilinear faults lines by the method 100 may be enhanced by the use of grouping algorithms that generalize the Radon, or equivalent, transform. In an exemplary embodiment, such grouping algorithms search for smooth curves in scattered binary data.
Referring now to FIGS. 8, 9, 10 and 11, an exemplary embodiment of a method 800 of determining channels in seismic images is provided in which in 802 a seismic image 802a is input.
In an exemplary embodiment, the method 800 then smooths the seismic image 802 in 804, and detects edges in the smoothed seismic image in 806. In an exemplary embodiment, the method 800 then generates an image 808a of the edges in the smoothed seismic image in 808.
In an exemplary embodiment, the method 800 then identifies smooth curves in the edges of the smoothed seismic image 808a in 810. In an exemplary embodiment, the method 800 then identifies parallel smooth curves 812a in the smoothed seismic image 808a in 812. In an exemplary embodiment, the parallel smooth curves 812a define channels in the smoothed seismic image 808a. 
In an exemplary embodiment, the teachings of the method 800 may be applied to seismic images having N-dimensions.
In an exemplary embodiment, the method 800 may be used to determine one or more channels.
In an exemplary embodiment, in the method 800, grouping algorithms may be used in order to enhance the determination of channels.
In an exemplary embodiment, in the method 800, contextual geological may be used such as, for example, proximity to other geological structures and constraints on orientation based upon rock physics, may be used to enhance the determination of channels.
In an exemplary embodiment, the method 800 may be used to determine linear, non-linear, or curvilinear channels. In an exemplary embodiment, the determination of curvilinear channels by the method 800 may be enhanced by the use of grouping algorithms. In an exemplary embodiment, such grouping algorithms search for smooth curves in scattered binary data.
In an exemplary embodiment, the smoothing operations, 104 and 804, of the methods 100 and 800 may be implemented using one or more conventional techniques for reducing noise in an image. In an exemplary embodiment, the filtering provided in the smoothing operations, 104 and 804, of the methods 100 and 800 may, for example, include specifying a window size, specifying a filter matrix of the specified window size that defines a local averaging technique, and performing a convolution for each pixel of the image. In an exemplary embodiment, the filtering provided in the smoothing operations, 104 and 804, of the methods 100 and 800 may, for example, include anisotropic methods for smoothing the seismic images.
In an exemplary embodiment, the edge detection operations, 106 and 806, of the methods, 100 and 800, may, for example, include techniques for detecting discontinuities along spatial directions in an image. In an exemplary embodiment, such techniques for detecting discontinuities along spatial directions calculate the spatial derivatives of an image and search for sharp changes in the spatial derivatives. Examples of techniques for detecting discontinuities along spatial directions include the Sobel and Canny, the selection of which may depend upon the different continuous representations of the associated convolution operators.
In an exemplary embodiment, the determination of the radius of curvature operation 108 of the method 100 may be provided, for 2-dimensional images, using the following equation for determining the iso-contour curvature R(x,y) of a 2-dimensional function I(x,y) as follows:




  
    
      
        
          
            R
            ⁡
            
              (
              
                x
                ,
                y
              
              )
            
          
          =
          
            
            
              
                
                  (
                  
                    
                      I
                      x
                      2
                    
                    +
                    
                      I
                      y
                      2
                    
                  
                  )
                
                
                  3
                  
                    /
                    2
                  
                
              
              
                
                  
                    I
                    
                      x
                      ⁢
                      
                        
                      
                      ⁢
                      x
                    
                  
                  ⁢
                  
                    I
                    y
                    2
                  
                
                +
                
                  
                    I
                    
                      y
                      ⁢
                      
                        
                      
                      ⁢
                      y
                    
                  
                  ⁢
                  
                    I
                    x
                    2
                  
                
                -
                
                  2
                  ⁢
                  
                    I
                    x
                  
                  ⁢
                  
                    I
                    y
                  
                  ⁢
                  
                    I
                    
                      x
                      ⁢
                      
                        
                      
                      ⁢
                      y
                    
                  
                
              
            
            
          
        
        ⁢
        
          

        
        ⁢
        
          where
          ⁢
          
            :
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              I
              x
            
            ≡
            
              
                ∂
                I
              
              
                ∂
                x
              
            
          
          ;
        
        ⁢
        
          

        
        ⁢
        
          
            
              I
              y
            
            ≡
            
              
                ∂
                I
              
              
                ∂
                y
              
            
          
          ;
        
        ⁢
        
          

        
        ⁢
        
          
            
              I
              
                x
                ⁢
                
                  
                
                ⁢
                x
              
            
            ≡
            
              
                
                  ∂
                  2
                
                ⁢
                I
              
              
                
                  ∂
                  x
                
                ⁢
                
                  ∂
                  x
                
              
            
          
          ;
        
        ⁢
        
          

        
        ⁢
        
          
            
              I
              
                y
                ⁢
                
                  
                
                ⁢
                y
              
            
            ≡
            
              
                
                  ∂
                  2
                
                ⁢
                I
              
              
                
                  ∂
                  y
                
                ⁢
                
                  ∂
                  y
                
              
            
          
          ;
          
            
          
          ⁢
          and
        
        ⁢
        
          

        
        ⁢
        
          
            I
            
              x
              ⁢
              
                
              
              ⁢
              y
            
          
          ≡
          
            
              
                
                  ∂
                  2
                
                ⁢
                I
              
              
                
                  ∂
                  x
                
                ⁢
                
                  ∂
                  y
                
              
            
            .
          
        
      
    
    
      
        (
        1
        )
      
    
  




In an exemplary embodiment, the Radon transformation operation 118 of the method 100 provides a projection of the image intensity along straight lines in all possible directions and orientations. In an exemplary embodiment, other transformations that provide an indication of the image intensity in one or more spatial directions, may be used instead of, or in addition to the Radon transformation.
In an exemplary embodiment, one or more of the operations of the methods 100 and 800 include clustering and/or grouping operations in order to enhance the processing of the seismic images in order to determine non-linear and/or curvilinear faults and/or channels in the seismic images. In an exemplary embodiment, during the operation of the methods 100 and 800, edge segments are obtained by segmentation algorithms that break up edges into edge segments. In an exemplary embodiment, grouping operations are included in the implementation of the methods 100 and 800 in order to group such edge segments into meaningful geometries. In an exemplary embodiment, the edge segments may be clustered into coherent groups. In this exemplary embodiment, the edge segments that are associated with the same cluster may indicate the presence of boundaries of a geological feature.
A method of determining the presence of one or more faults within a seismic image has been described that includes inputting a seismic image; smoothing the input seismic image; detecting edges in the smoothed seismic image; projecting an intensity of the edges in a plurality of spatial directions; determining a direction of maximum intensity of the projections of the edges; and for the directions of maximum intensity defining a fault line in the seismic image. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises detecting edges having a radius of curvature greater than a predetermined value. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises detecting edges having a radius of curvature less than or equal to a predetermined value. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises detecting edges having a radius of curvature less than or equal to a predetermined value; and setting the radius of curvature of the edges having a radius of curvature less than or equal to a predetermined value equal to zero. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises detecting edges having a radius of curvature greater than a predetermined value; detecting edges having a radius of curvature less than or equal to a predetermined value; and setting the radius of curvature of the edges having a radius of curvature less than or equal to a predetermined value equal to zero. In an exemplary embodiment, projecting an intensity of the edges in a plurality of spatial directions comprises performing a Radon transform on the seismic image. In an exemplary embodiment, determining a direction of maximum intensity of the projections of the edges comprises determining a local maxima of the Radon transform. In an exemplary embodiment, the fault line comprises a linear fault line. In an exemplary embodiment, the fault line comprises a non-linear fault line. In an exemplary embodiment, the fault line comprises a curvilinear fault line. In an exemplary embodiment, the method further comprises searching for smooth curves in the seismic image. In an exemplary embodiment, the seismic image comprises an n-dimensional image, where n is greater than or equal to 2. In an exemplary embodiment, smoothing the input seismic image comprises filtering noise from the input seismic image. In an exemplary embodiment, filtering noise from the input seismic image comprises specifying a filter window size; specifying a filter matrix; and performing a convolution for each pixel of the input seismic image. In an exemplary embodiment, filtering noise from the input seismic image comprises anisotropically filtering noise from the input seismic image. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises calculating spatial derivatives in the smoothed seismic image. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises determining a radius of curvature of edges within the smoothed seismic image. In an exemplary embodiment, determining a radius of curvature of edges within the smoothed seismic image comprises determining an iso-contour curvature of edges within the smoothed seismic image. In an exemplary embodiment, determining an iso-contour curvature of edges within the smoothed seismic image comprises determining an iso-contour curvature of edges R(x,y) of a 2-dimensional function I(x,y) within the smoothed seismic image using the following equation:




  
    R
    ⁡
    
      (
      
        x
        ,
        y
      
      )
    
  
  =
  
    
      
      
        
          
            (
            
              
                I
                x
                2
              
              +
              
                I
                y
                2
              
            
            )
          
          
            3
            
              /
              2
            
          
        
        
          
            
              I
              
                x
                ⁢
                
                  
                
                ⁢
                x
              
            
            ⁢
            
              I
              y
              2
            
          
          +
          
            
              I
              
                y
                ⁢
                
                  
                
                ⁢
                y
              
            
            ⁢
            
              I
              x
              2
            
          
          -
          
            2
            ⁢
            
              I
              x
            
            ⁢
            
              I
              y
            
            ⁢
            
              I
              
                x
                ⁢
                
                  
                
                ⁢
                y
              
            
          
        
      
      
    
    .
  




A method of determining the presence of one or more faults within a seismic image has been described that includes inputting a seismic image; smoothing the input seismic image; detecting edges in the smoothed seismic image; determining edges within the smoothed seismic image having a large radius of curvature; generating an image of the large radius of curvature edges; projecting an intensity of the image of the large radius of curvature edges in a plurality of spatial directions; determining a direction of maximum intensity of the projections of the image of the large radius of curvature edges; and for the directions of maximum intensity defining a fault line in the seismic image.
A method of determining the presence of one or more channels within a seismic image has been described that includes inputting a seismic image; smoothing the input seismic image; detecting edges in the smoothed seismic image; identifying smooth curves within the detected edges of the smoothed seismic image; identifying one or more sets of parallel smooth curves in the smoothed seismic image; and for the sets of parallel smooth curves defining one or more channels in the seismic image. In an exemplary embodiment, the channels comprise linear channels. In an exemplary embodiment, the channels comprise non-linear channels. In an exemplary embodiment, the channels comprise curvilinear channels. In an exemplary embodiment, the method further includes searching for smooth curves in the seismic image. In an exemplary embodiment, the seismic image comprises an n-dimensional image, where n is greater than or equal to 2. In an exemplary embodiment, smoothing the input seismic image comprises filtering noise from the input seismic image. In an exemplary embodiment, filtering noise from the input seismic image comprises specifying a filter window size; specifying a filter matrix; and performing a convolution for each pixel of the input seismic image. In an exemplary embodiment, filtering noise from the input seismic image comprises anisotropically filtering noise from the input seismic image. In an exemplary embodiment, detecting edges in the smoothed seismic image comprises calculating spatial derivatives in the smoothed seismic image.
It is understood that variations may be made in the foregoing without departing from the scope of the invention. Furthermore, the elements and teachings of the various illustrative embodiments may be combined in whole or in part in some or all of the illustrative embodiments. In addition, one or more of the elements and teachings of the various illustrative embodiments may be omitted, at least in part, and/or combined, at least in part, with one or more of the other elements and teachings of the various illustrative embodiments.
Although illustrative embodiments of the invention have been shown and described, a wide range of modification, changes and substitution is contemplated in the foregoing disclosure. In some instances, some features of the present invention may be employed without a corresponding use of the other features. Accordingly, it is appropriate that the appended claims be construed broadly and in a manner consistent with the scope of the invention.

## Claims

1. A computer implemented method of determining the presence of one or more faults within a seismic image, comprising:
inputting a seismic image;
smoothing the input seismic image;
detecting edges in the smoothed seismic image;
projecting an intensity of the edges in a plurality of spatial directions;
determining a direction of maximum intensity of the projections of the edges; and
for the directions of maximum intensity defining a fault line in the seismic image.

2. The method of claim 1, wherein detecting edges in the smoothed seismic image comprises:
detecting edges having a radius of curvature greater than a predetermined value.

3. The method of claim 1, wherein detecting edges in the smoothed seismic image comprises:
detecting edges having a radius of curvature less than or equal to a predetermined value.

4. The method of claim 3, wherein detecting edges in the smoothed seismic image comprises:
detecting edges having a radius of curvature less than or equal to a predetermined value; and
setting the radius of curvature of the edges having a radius of curvature less than or equal to a predetermined value equal to zero.

5. The method of claim 1, wherein detecting edges in the smoothed seismic image comprises:
detecting edges having a radius of curvature greater than a predetermined value;
detecting edges having a radius of curvature less than or equal to a predetermined value; and
setting the radius of curvature of the edges having a radius of curvature less than or equal to a predetermined value equal to zero.

6. The method of claim 1, wherein projecting an intensity of the edges in a plurality of spatial directions comprises:
performing a Radon transform on the seismic image.

7. The method of claim 6, wherein determining a direction of maximum intensity of the projections of the edges comprises:
determining local maxima of the Radon transform.

8. The method of claim 1, further comprising:
searching for smooth curves in the seismic image.

9. The method of claim 1, wherein the seismic image comprises an n-dimensional image, where n is greater than or equal to 2.

10. The method of claim 1, wherein smoothing the input seismic image comprises:
filtering noise from the input seismic image.

11. The method of claim 10, wherein filtering noise from the input seismic image comprises:
specifying a filter window size;
specifying a filter matrix; and
performing a convolution for each pixel of the input seismic image.

12. The method of claim 10, wherein filtering noise from the input seismic image comprises:
anisotropically filtering noise from the input seismic image.

13. The method of claim 1, wherein detecting edges in the smoothed seismic image comprises:
calculating spatial derivatives in the smoothed seismic image.

14. The method of claim 1, wherein detecting edges in the smoothed seismic image comprises:
determining a radius of curvature of edges within the smoothed seismic image.

15. The method of claim 14, wherein determining a radius of curvature of edges within the smoothed seismic image comprises:
determining an iso-contour curvature of edges within the smoothed seismic image.

16. The method of claim 15, wherein determining an iso-contour curvature of edges within the smoothed seismic image comprises:
determining an iso-contour curvature of edges R(x,y) of a 2-dimensional function I(x,y) within the smoothed seismic image using the following equation:





  
    R
    ⁡
    
      (
      
        x
        ,
        y
      
      )
    
  
  =
  
    
      
      
        
          
            (
            
              
                I
                x
                2
              
              +
              
                I
                y
                2
              
            
            )
          
          
            3
            
              /
              2
            
          
        
        
          
            
              I
              
                x
                ⁢
                
                  
                
                ⁢
                x
              
            
            ⁢
            
              I
              y
              2
            
          
          +
          
            
              I
              
                y
                ⁢
                
                  
                
                ⁢
                y
              
            
            ⁢
            
              I
              x
              2
            
          
          -
          
            2
            ⁢
            
              I
              x
            
            ⁢
            
              I
              y
            
            ⁢
            
              I
              
                x
                ⁢
                
                  
                
                ⁢
                y
              
            
          
        
      
      
    
    .

17. A computer implemented method of determining the presence of one or more channels within a seismic image, comprising:
inputting a seismic image;
smoothing the input seismic image;
detecting edges in the smoothed seismic image;
identifying smooth curves within the detected edges of the smoothed seismic image;
identifying one or more sets of parallel smooth curves in the smoothed seismic image; and
for the sets of parallel smooth curves defining one or more channels in the seismic image.

18. The method of claim 17, further comprising:
searching for smooth curves in the seismic image.

19. The method of claim 17, wherein the seismic image comprises an n-dimensional image, where n is greater than or equal to 2.

20. The method of claim 17, wherein smoothing the input seismic image comprises:
filtering noise from the input seismic image.

21. The method of claim 20, wherein filtering noise from the input seismic image comprises:
specifying a filter window size;
specifying a filter matrix; and
performing a convolution for each pixel of the input seismic image.

22. The method of claim 20, wherein filtering noise from the input seismic image comprises:
anisotropically filtering noise from the input seismic image.

23. The method of claim 17, wherein detecting edges in the smoothed seismic image comprises:
calculating spatial derivatives in the smoothed seismic image.

24. A computer implemented method of determining the presence of one or more faults within a seismic image, comprising:
inputting a seismic image;
smoothing the input seismic image;
detecting edges in the smoothed seismic image;
determining an iso-contour curvature of edges within the smoothed seismic image; and
determining one or more fault lines from the iso-contour curvature.

25. The method of claim 24, further comprising:
co-rendering the iso-contour curvature with the seismic image; and
determining one or more fault lines from the co-rendering.