Automatic tracking of faults by slope decomposition

0Associated
Cases 
0Associated
Defendants 
0Accused
Products 
0Forward
Citations 
0
Petitions 
0
Assignments
First Claim
1. A computerimplemented method for automatically tracking faults in a 2D imaged seismic crosssection or a 3D imaged seismic data volume, comprising:
 (a) decomposing, with a computer, the imaged seismic data into slopes, wherein the imaged seismic data has previously undergone depth migration; and
(b) forming, with a computer, a faulthighlighted data volume or crosssection from voxels corresponding to fault discontinuities in the imaged seismic data having slopes that span a broader range of slopes than other voxels in the imaged seismic data;
(c) selecting, with a computer, one or more initial seeds for fault surfaces or fault lines within the faulthighlighted data volume or crosssection;
(d) generating, with a computer, one or more fault contours in the faulthighlighted data volume or crosssection starting from the initial seeds;
(e) displaying, with a computer, a connected, smooth fault surface or line based on the one or more fault contours and(f) exploring for hydrocarbons based at least in part upon the generated fault contours and/or the displayed smooth fault surface or line.
0 Assignments
0 Petitions
Accused Products
Abstract
Method for locating fault lines or surfaces in 2D or 3D seismic data based on the fact that fault discontinuities in the space domain span a wide range in a local slowness (slope) domain, whereas other dipping events in the space domain data, such as noise, tend to be coherent, and hence will appear focused in the slowness dimension. Therefore, the method comprises decomposing the seismic data (102) by a transformation to the local slowness domain, preferably using Gaussian slowness period packets as the local slowness or slope decomposition technique, thereby avoiding problems with the data stationary assumption. In the local slowness domain, faults may be identified (104) using the principle mentioned above, i.e. that faults are represented as a truncation in the space domain data, hence they will appear broadband in the slowness dimension.
30 Citations
No References
Method and Apparatus For Analyzing ThreeDimensional Data  
Patent #
US 20100274543A1
Filed 11/13/2008

Current Assignee
Timothy A. Chartrand, Paul A. Dunn, Kelly G. Walker

Sponsoring Entity
Timothy A. Chartrand, Paul A. Dunn, Kelly G. Walker

Systems and methods for enhancing a seismic data image  
Patent #
US 7,702,463 B2
Filed 12/12/2007

Current Assignee
Landmark Graphics Corporation

Sponsoring Entity
Landmark Graphics Corporation

Method For Geophysical and Geological Interpretation of Seismic Volumes In The Domains of Depth, Time, and Age  
Patent #
US 20100149917A1
Filed 11/20/2009

Current Assignee
Exxonmobil Upstream Research Company

Sponsoring Entity
Exxonmobil Upstream Research Company

Estimating and Using Slowness Vector Attributes in Connection with a MultiComponent Seismic Gather  
Patent #
US 20090003132A1
Filed 06/29/2007

Current Assignee
WesternGeco LLC

Sponsoring Entity
WesternGeco LLC

TimeSpace Varying Spectra for Seismic Processing  
Patent #
US 20090292475A1
Filed 05/21/2008

Current Assignee
PRIME GEOSCIENCE CORPORATION

Sponsoring Entity
PRIME GEOSCIENCE CORPORATION

Processing of seismic data using the Stransform  
Patent #
US 7,617,053 B2
Filed 05/14/2007

Current Assignee
Calgary Scientific Inc.

Sponsoring Entity
Calgary Scientific Inc.

Rapid Method for Reservoir Connectivity Analysis Using a Fast Marching Method  
Patent #
US 20080154505A1
Filed 04/10/2006

Current Assignee
ChulSung Kim, Mark Dobin

Sponsoring Entity
ChulSung Kim, Mark Dobin

System and Method for Displaying Seismic Horizons with Attributes  
Patent #
US 20080285384A1
Filed 10/19/2006

Current Assignee
Emerson Paradigam Holding LLC

Sponsoring Entity
Emerson Paradigam Holding LLC

Method for rapid fault interpretation of fault surfaces generated to fit threedimensional seismic discontinuity data  
Patent #
US 20070078604A1
Filed 09/17/2003

Current Assignee
Exxonmobil Upstream Research Company

Sponsoring Entity
Exxonmobil Upstream Research Company

Method and apparatus for seismic feature extraction  
Patent #
US 20060122780A1
Filed 11/10/2003

Current Assignee
GEOENERGY INC.

Sponsoring Entity
GEOENERGY INC.

Method and apparatus for internetworked wireless integrated network sensor (WINS) nodes  
Patent #
US 6,859,831 B1
Filed 10/04/2000

Current Assignee
Benhov GmbH LLC

Sponsoring Entity
Sensoria

Method for performing stratrigraphicallybased seed detection in a 3D seismic data volume  
Patent #
US 20040062145A1
Filed 09/03/2003

Current Assignee
Exxonmobil Upstream Research Company

Sponsoring Entity
Exxonmobil Upstream Research Company

Method for collecting data using compact internetworked wireless integrated network sensors (WINS)  
Patent #
US 6,735,630 B1
Filed 10/04/2000

Current Assignee
Benhov GmbH LLC

Sponsoring Entity
Sensoria

Apparatus for internetworked hybrid wireless integrated network sensors (WINS)  
Patent #
US 6,826,607 B1
Filed 10/04/2000

Current Assignee
Benhov GmbH LLC

Sponsoring Entity
Sensoria

Method and apparatus for distributed signal processing among internetworked wireless integrated network sensors (WINS)  
Patent #
US 6,832,251 B1
Filed 10/04/2000

Current Assignee
Benhov GmbH LLC

Sponsoring Entity
Sensoria

Method for characterization of multiscale Geometric attributes  
Patent #
US 20020035443A1
Filed 03/09/2001

Current Assignee
Exxonmobil Upstream Research Company

Sponsoring Entity
Exxonmobil Upstream Research Company

Automated seismic fault detection and picking  
Patent #
US 6,018,498 A
Filed 09/02/1998

Current Assignee
ConocoPhillips

Sponsoring Entity
Phillips Petroleum Company

Method of seismic attribute generation and seismic exploration  
Patent #
US 5,940,778 A
Filed 07/31/1997

Current Assignee
BP Amoco

Sponsoring Entity
BP Amoco

Method and apparatus for automatically identifying fault cuts in seismic data using a horizon time structure  
Patent #
US 5,999,885 A
Filed 02/06/1997

Current Assignee
Schlumberger Technology Corporation

Sponsoring Entity
Schlumberger Technology Corporation

Earthquake early warning system  
Patent #
US 5,625,138 A
Filed 02/23/1995

Current Assignee
Jack D. Elkins

Sponsoring Entity
Jack D. Elkins

System for evaluating seismic sequence lithology and property, and for evaluating risk associated with predicting potential hydrocarbon reservoir, seal, trap or source  
Patent #
US 5,475,589 A
Filed 07/08/1993

Current Assignee
SPIRAL HOLDING INC.

Sponsoring Entity
SPIRAL HOLDING INC.

Method and apparatus for finding horizons in 3D seismic data  
Patent #
US 5,251,184 A
Filed 07/23/1992

Current Assignee
Landmark Graphics Corporation

Sponsoring Entity
Landmark Graphics Corporation

Method for finding horizons in 3D seismic data  
Patent #
US 5,153,858 A
Filed 07/09/1991

Current Assignee
Landmark Graphics Corporation

Sponsoring Entity
Landmark Graphics Corporation

Timespace varying spectra for seismic processing  
Patent #
US 8,185,316 B2
Filed 05/21/2008

Current Assignee
PRIME GEOSCIENCE CORPORATION

Sponsoring Entity
PRIME GEOSCIENCE CORPORATION

METHOD AND DEVICE FOR WAVE FIELDS SEPARATION IN SEISMIC DATA  
Patent #
US 20130021873A1
Filed 07/18/2011

Current Assignee
CGG Services U.S. Inc.

Sponsoring Entity
CGG Services U.S. Inc.

Method of predicting connectivity between parts of a potential hydrocarbon reservoir and analyzing 3D data in a subsurface region  
Patent #
US 8,370,122 B2
Filed 11/13/2008

Current Assignee
Exxonmobil Upstream Research Company

Sponsoring Entity
Exxonmobil Upstream Research Company

Method Of Generating and Combining Multiple Horizons To Determine A Seismic Horizon And Its Uncertainty  
Patent #
US 20130121111A1
Filed 10/12/2012

Current Assignee
ChulSung Kim, Mark W. Dobin

Sponsoring Entity
ChulSung Kim, Mark W. Dobin

Method and System for Geophysical Modeling of Subsurface Volumes  
Patent #
US 20140278311A1
Filed 02/20/2014

Current Assignee
John S. Lien, Martin Terrell, Matthias Imhof, Pavel Dimitrov

Sponsoring Entity
John S. Lien, Martin Terrell, Matthias Imhof, Pavel Dimitrov

Automated Interpretation Error Correction  
Patent #
US 20140358445A1
Filed 05/07/2014

Current Assignee
Ethan Nowak, Matthias Imhof, Pavel Dimitrov

Sponsoring Entity
Ethan Nowak, Matthias Imhof, Pavel Dimitrov

Method of generating and combining multiple horizons to determine a seismic horizon and its uncertainty  
Patent #
US 9,239,398 B2
Filed 10/12/2012

Current Assignee
Exxonmobil Upstream Research Company

Sponsoring Entity
Exxonmobil Upstream Research Company

12 Claims
 1. A computerimplemented method for automatically tracking faults in a 2D imaged seismic crosssection or a 3D imaged seismic data volume, comprising:
(a) decomposing, with a computer, the imaged seismic data into slopes, wherein the imaged seismic data has previously undergone depth migration; and (b) forming, with a computer, a faulthighlighted data volume or crosssection from voxels corresponding to fault discontinuities in the imaged seismic data having slopes that span a broader range of slopes than other voxels in the imaged seismic data; (c) selecting, with a computer, one or more initial seeds for fault surfaces or fault lines within the faulthighlighted data volume or crosssection; (d) generating, with a computer, one or more fault contours in the faulthighlighted data volume or crosssection starting from the initial seeds; (e) displaying, with a computer, a connected, smooth fault surface or line based on the one or more fault contours and (f) exploring for hydrocarbons based at least in part upon the generated fault contours and/or the displayed smooth fault surface or line.  View Dependent Claims (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)
1 Specification
This application claims the benefit of U.S. Provisional Patent Application No. 61/898,253, filed Oct. 31, 2013, entitled AUTOMATIC TRACKING OF FAULTS BY SLOPE DECOMPOSITION, the entirety of which is incorporated by reference herein.
This application is related to U.S. application Ser. No. 13/639,802, published as US 2013/0044568, entitled Seismic Signal Processing with Gaussian SlownessPeriod Packets, which is hereby incorporated by reference in its entirety.
This disclosure relates generally to the field of geophysical prospecting for hydrocarbons and, more particularly to seismic data processing. Specifically, the disclosure concerns a computerautomated method for tracking faults in seismic data.
Many seismic signal processing techniques are applied to seismic data to enhance a migrated image, including regularization to create unrecorded traces needed by many processing algorithms, coherent noise attenuation to remove energy that does not contribute to the image, and random noise attenuation to enhance coherent events both before and after imaging. Often the underlying assumptions behind many of these signal processing techniques include an assumption of stationarity: that the events are planar in nature and that their dip or frequency content does not change with position. In reality, seismic data are nonstationary; they contain events with curvature and the frequency content changes as the recording time increases. This problem is well known, and several methods to address nonstationary data do exist. These include: breaking up the problem into overlapping spatialtemporal windows that are assumed to be locally stationary followed by processing and reassembly; the use of nonstationary filters that vary with space and time; and methods like the curvelet transform that expand the data into a compressible overrepresentation. This disclosure presents an alternative to these methods, an alternative that allows for more flexibility in handling nonstationarity in the data. See also patent application publication WO 2011/139411.
PatchBased Methods
The most common way to apply signal processing algorithms that assume a stationary input to data that are not stationary is (see
NonStationary Filtering
An alternative to the patchbased approach is to solve for filters that vary as a function of position. One example of this is the use of nonstationary predictionerror filters for either interpolation (Crawley et al., SEG Expanded Abstracts, 1999, Vol. 18, Pages 11541157) or signal/noise separation (Guitton, Geophysics, 2005, Vol. 70, Pages V97V107). This spatially variable filter is estimated on the entire dataset simultaneously by solving a large inverse problem for all of the filter coefficients. Since these filters vary with position, the number of unknown filter coefficients can be larger than the number of known data points, creating an underdetermined problem that is solved by adding regularization to create a smoothly nonstationary predictionerror filter. This has the benefit that the filter varies smoothly as a function of position and does not have the problem with visible boundaries that the patchbased approach does. However, creating a nonstationary filter is nonunique, so many of the benefits of a predictionerror filter that depend on solving an overdetermined unique problem are gone, such as a guarantee of minimum phase, among other benefits. The size of the filter also scales poorly with the number of dimensions involved, making higherdimensional filtering computationally expensive.
Another nonstationary filter often used is a local planewave destructor filter (Fomel, Geophysics 67, 19461960, (2002)) or a steering filter (Clapp, Geophysics 69, 533546 (2004)). These filters are stable and can be used to estimate local slope. This method is efficient, but when dealing with data with multiple conflicting slopes has difficulty and can alternate rapidly between the two possible slopes. Spatially aliased data also present problems as there are multiple possible solutions to the slope estimation problem.
A more recent approach to dealing with nonstationarity is from Hale (SEG Expanded Abstracts 26, 31603164 (2006)). Hale efficiently computes Gaussian windows in time and space using both the separability of multidimensional Gaussians as well as recursive filtering that he uses to solve for local crosscorrelations and autocorrelations of data, which he can then use to generate local predictionerror filters. This method works well, but assumes that the data are wellsampled. Hale uses this same efficient Gaussian smoothing to smooth the output of slope estimation from a structure tensor. Current uses of Gaussian filtering in time and space are limited to methods that depend on efficient cross and auto correlations, and for smoothing the output of other processes. Moving to a domain where prior information is more easily enforced, such as the frequencywavenumber domain, is not straightforward with this approach.
The S transform (Stockwell, Mansinha and Lowe, IEEE Signal Processing, 1996, Vol. 44, Pages 9981001; Pinnegar and Mansinha, Geophysics, 2003, Vol. 68, Pages 381385; and U.S. Pat. No. 7,617,053 to Pinnegar et al.) uses modulated Gaussian functions for time frequency analysis, but only along in a single dimension as an alternative to a spectrogram.
Gabor filtering (Daughman, IEE Trans. On Acoustics, Speech, and Signal Processing, 1988, Vol. 36, Pages 11691179), used in image analysis and edge detection, uses modulated Gaussian functions in multiple spatial dimensions. These filters are typically parameterized by a rotation and dilation, making them nonseparable. Since applications are limited to 2 or 3 spatial dimensions, separability and parallelization are less of an issue than with higherdimensional seismic data. In addition, using this type of filter on data with time and spatial axes might be confusing, as the rotation would span different frequency ranges depending on the rotation.
Curvelets
Another approach to the problem of nonstationarity is to use local basis functions. The curvelet transform (U.S. Pat. No. 7,840,625 by Candes et al) is a transform that is a partition of the frequencywavenumber domain where the data are broken up into various sized windows in scale and angle according to a parabolicdyadic scaling law, and then each window is returned to the timespace domain on a different grid by an inverse Fourier transform to produce an output in angle, scale, time, and space.
Curvelets have largely been used for compressive sensing, where the compressive qualities of curvelets on seismic data are used to eliminate noise or interpolate missing data. This method works provided that certain underlying assumptions are fulfilled. For the interpolation case, this is that the data are randomly sampled, producing a lowamplitude blur in both the frequencywavenumber and curvelet domains. For a signal/noise case, an adaptive subtraction of externally modeled noise can take place in the curvelet domain. Curvelets address the problem of nonstationarity, and work when simple operations such as thresholding or subtraction take place in the curvelet domain.
Applying operators across this curvelet domain is difficult, as the grid spacing and orientation differ for each scale and angle. Additionally, the parameterization of the curvelet space is difficult, as the parameters of scale and number of angles are not intuitive. Finally, curvelets become computationally expensive in higher dimensions, where the gridding and overrepresentation of the data increase greatly.
Many methods have been published for fault detection and tracking in seismic data, including the following: U.S. Pat. No. 6,018,498 to Neff, et al.; PCT patent application publication WO 2009/082545 by Kumaran et al.; U.S. patent application publication US 2006/0122780 by Cohen et al. (other fault detection references are listed by Cohen in his paragraph 35); U.S. patent application publication US 2002/0035443 by Matteucci et al; and U.S. Pat. No. 5,940,778 to Marfurt et al.
No single method deals with nonstationarity in a way that is stable, easily parallelizable, easily scalable to higher dimensions, and easily incorporates prior knowledge. The present invention satisfies these requirements.
One embodiment of this disclosure is a method for producing a representation of a subsurface region from multidimensional seismic data, comprising
(a) using a cascade of filtering operations to decompose the seismic data into components in a frequencywavenumber domain wherein the data are represented in terms of tiled windows;
(b) applying one or more processing operations to the decomposed data, said processing operations being designed to enhance a representation of the subsurface region from the seismic data;
(c) applying the filtering operations'"'"' adjoint operations to the decomposed data after the processing operations, then summing data components weighted by normalization factors computed to produce a flat impulse response; and
(d) using the weighted and summed data from (c) to generate a representation of the subsurface region.
The preceding method is advantageous where the processing operations require an assumption of local stationarity and the data are not stationary, that is that the dip or frequency spectra of the seismic data change with position and time. The cascade of filtering operations may be a series of 1D Gaussian filters modulated by complex exponentials, creating a series of Gaussian windows in the frequencywavenumber domain.
One embodiment of the present invention is a computerimplemented method for automatically tracking faults in a 2D or 3D seismic data volume or crosssection, comprising: (a) decomposing the seismic data into local slownesses for seismic data recorded in time, or slopes if the seismic data have undergone depth migration, using a programmed computer; and (b) identifying voxels in the seismic data which when decomposed into slownesses or slopes span a broader range of slownesses or slopes than other voxels in the seismic data, thereby forming a faulthighlighted data volume or crosssection.
The method may further comprise: (c) selecting one or more initial seeds for fault surfaces or fault lines within the faulthighlighted data volume or crosssection based on the fault highlights; (d) generating one or more fault contours in the faulthighlighted data volume or crosssection starting from the initial seeds; and (e) outputting a connected, smooth fault surface or line based on the one or more fault contours.
As with any seismic data processing method, the invention in practical applications is highly automated, i.e. is performed with the aid of a computer programmed in accordance with the disclosures herein.
The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:
The invention will be described in connection with example embodiments. However, to the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claims.
The present disclosure presents an alternative to existing methods to address nonstationary data, an alternative that allows for more flexibility in handling nonstationarity in the data, by allowing the incorporation of prior knowledge along many dimensions in the data: spatial, temporal, and their Fourier equivalents. This disclosure explains how to construct a representation that scales in an easily parallelizable way to higher dimensions, and can be used for interpolation, signal/noise separation, and decomposition of seismic data. This invention is applicable to multidimensional seismic signal processing, both before and after imaging. The invention allows for a straightforward method to process data in frequency, time, and multiple spatial axes and slopes, or any subset of these, all simultaneously.
The method of
Separable Multidimensional Modulated Gaussian Filtering
The derivation of the onedimensional Gaussian filters starts with a pair of desired multidimensional Gaussian windows in frequency (f) and wavenumber (k) with a central frequency and wavenumber f_{0 }and k_{0}, and half widths at half maxima h_{f }and h_{k }for the frequency and wavenumber axes, respectively, where
This two dimensional filter can be separated into two products of two Gaussian functions that have been convolved with shifted delta functions so that
This sequential application, or cascade, of filters can be performed in the time and space domain as the convolution of two 1D filters, each composed of a Gaussian multiplied by a complex exponential, where
If f_{0 }is larger than h_{f}, this addition of two cascaded series of filters can be approximated by taking two times the real portion of the output of one of these filters, so that
In an alternate embodiment, in the case where very low and zero frequencies are desired, the above implementation can be replaced by applying both cascades of filters in equation 3 separately to produce separate complex outputs.
The half widths at half maxima of the filters can be expressed either in time/space or frequency/wavenumber using the relation h_{f}=ln 2/πh_{t }depending on the desired input parameterization and tiling, discussed in detail next.
Implicit Tiling of the Frequency/Wavenumber Domain
Another feature of this invention is a method of tiling of the frequency wavenumber domain in a manner best suited for prestack seismic data.
Frequency sampling should be such that the impulse response of the cascaded forward and adjoint operator is flat. For the forward operator, the Gaussians should be shifted by h_{f }or less relative to each other such that the windows sum to produce a flat response. When n_{f }Gaussian windows are applied across the frequency axis with a halfwidth in time of h_{t}, the corresponding frequency halfwidth after the forward and adjoint operations is √{square root over (2)} ln 2/πh_{t}. Setting the frequency sampling equal to this value gives
as a normalization factor to rescale the frequency response to unity after the forward and inverse transforms. In equation (5), f_{max}, and f_{min }are the specified minimum and maximum central frequencies of the timedomain Gaussian windows. The tiling in wavenumber for prestack seismic data is best suited to a regular tiling in slowness, with spacing in slowness such that the tiling along wavenumber at the maximum frequency produces a flat impulse response. This maximum spacing is given by the halfwidth of the operator, either in wavenumber or in space, and is √{square root over (2)}h_{k }or √{square root over (2)} ln 2/πh_{x }and the sampling interval in p, Δp, is this value divided by the maximum frequency. (p is the slowness, which is related to the wavenumber k by k=fp.) At each output frequency, this constant sampling in p will produce a different normalization factor roughly equivalent to the ω scaling needed by a taup transform, and is
where n_{p }slownesses are sampled. This product of the present disclosure may be called Gaussian slownessperiod packets (“GaSPs”), represented by g in the equations above. A collection of these GaSPs” is shown in
Combining the two concepts of the separable Gaussian functions and the frequencyslowness tiling produces a transform that can be used in the following steps, as described in
At step 101, given desired input parameters, including a desired frequency and slowness range and either time or spatial halfwidth of operators, desired frequency and slowness discrimination, or the number of desired central frequencies and slownesses, define a tiling so that the impulse response is flat after a forward and inverse transform, which corresponds to a spacing equal to 1/√{square root over (2)} times the halfwidth of the operators in frequency and in slowness.
At step 102, precompute the necessary 1D filters required for the tiling, which is (approximately n_{f}+n_{f}n_{p }filters) as well as the n_{f}+1 normalization factors to produce a flat impulse response after forward and adjoint transforms.
At step 104, for each desired central frequency and slowness defined in step 101, apply a cascade of 1D modulated Gaussian filters to the seismic data 103 to generate a filtered version of the data for each central frequency and central slowness. It may be noted that the terms slope, dip, and slowness may be used more or less interchangeably, with slowness more applicable to an unimaged seismic data case, where the axes are time and space, whereas slope or dip are commonly used when dealing with imaged or migrated data, where the axes are often x, y, and depth.
At step 106, apply an operation or operations to the transformed datasets. These operations could include elementbyelement operations such as muting or thresholding, operators that span any of the axes of the datasets: the time axis, spatial axes, slope axes, or the frequency axis, such as applying filters along any of these axes. In general, these are data processing operations that are usually performed in the transform domain where they must assume data stationarity. Alternately, this decomposition could be done for two or more datasets and operations performed on the combination of datasets.
At step 108, apply the adjoint operation by reapplying the same filters to the decomposed data, followed by a sum weighted by the normalization factors computed in step 102.
The flowchart and examples below are for the 2D case, with time and space transforming to frequency and wavenumber. However, it can be appreciated that the separable Gaussian filters and transform tiling described here really apply to any number of dimensions, which may include as many as seven: three spatial dimensions for each of the source and the receiver, and time, although if source and receiver are located on the Earth'"'"'s surface, the number reduces to five. In fact, it is a particularly advantageous feature of this inventive method that it has beneficial scaling and efficiencies in higher dimensions that other local transforms lack.
BeyondAliasing Interpolation
A commonlyused interpolation method for prestack reflection seismic data is fk interpolation (Gülünay, Geophysics 68, 355369 (2003)) where regularlysampled data can be interpolated correctly despite the ShannonNyquist sampling criterion by incorporating prior information; in this case the prior information is a weighting function generated from the low frequencies present in the data. This is typically accomplished in the Fourier domain, and as such assumes stationarity, so the method is typically applied in overlapping spatialtemporal patches.
Gaussian slowness period packets can be used instead of these windows, and the problem can be reformulated so that the lowerfrequency GaSPs are used to constrain the higherfrequency GaSPs, and remove the aliased energy caused by the coarse sampling.
An example of this is shown in
Applying GaSPbased Gülünay factorofthree interpolation to the live traces in
Match Filtering of Multiples
Match filtering typically takes place using a filter in time and/or space to match a noise model to a dataset containing the noise as well as desired signal (Verschuur and Berkhout, Geophysics, 1992, Vol. 57, Pages 11661177). More recently, curveletbased adaptive subtraction has been used to fit a noise model to data, both using real (Herrmann et al., Geophysics, 2008, Vol. 73, Pages A17A21) and complexvalued curvelets (Neelamani et al., SEG Expanded Abstracts, 2008, Vol. 27, Pages 36503655). The matching filters can have difficulty discriminating overlapping signal from noise, but deal well with bulk shifts between the modeled and actual noise because of the length of the filter in time, while the curveletbased approaches decompose across scale and angle, but currently do not deal well with large shifts between the modeled and actual noise. By using matching filters on a GaSPdecomposed noise model and data, both overlapping slopes and frequencies and significant kinematic differences can be addressed.
The topleft panel of
Using a GaSPbased adaptive subtraction, where the input data and multiple model are both decomposed into pairs of GaSPs at the same frequencies and slownesses, followed by estimating and applying a set of matching filters with the same parameters as the standard approach for each pair of GaSPs, matching the data and noise independently at each slowness and frequency, followed by the adjoint transformation producing the matched multiples shown in the centerright panel of
These superior results can also be analyzed by viewing the results of both the GaSPbased and standard approaches in the GaSPdecomposed domain.
These results can also be examined at a single slowness and frequency at all spatial locations and times as in
Fault Detection
GaSPbased decomposition of the seismic data into a range of local slowness or slopes may also be used for fault detection, particularly adaptable to being performed in an automated manner using a computer. The inventive realization for how this could work is based on the observation that fault discontinuities in the space domain span a wide range in a local slowness (slope) domain. Imaged seismic volumes have a vertical axis of either time or depth, depending on the type of migration used. Slope or dip is the more accurate term when dealing with an image with a depth axis, whereas slowness is more accurate when the image has a time axis.
This inventive method preferably uses Gaussian slowness period packets as the local slowness or slope decomposition technique, but it is not limited to GaSPbased decomposition. Other alternatives that could be used to perform this local decomposition into slowness bands include curvelet'"'"'s, a local radon decomposition, or wave atoms. For curvelet decomposition, see, for example, Candes, et al., “Fast Discrete Curvelet Transforms,” Multiscale Modeling & Simulation 5, no. 3 (2006): 861899. (2005). For local radon decomposition, see, for example, Theune, et al., “Leastsquares local Radon Transforms for dipdependent GPR image decomposition,” Journal of Applied Geophysics 59, 224235 (2006). For wave atoms, see, for example, Demanet and Ying, “Wave Atoms and Sparsity of Oscillatory Patterns,” Appl. Comput. Harmon. Anal. 23.3, 368387 (2007). This paper defines wave atoms, details their implementation, and describes an application to sparse representation of oscillatory textures.
A 2D image may be locally decomposed into slopes to form a 3D volume of depth, space and slope. When a 3D seismic volume is locally decomposed into slopes, the 3D volume now becomes a 5D volume, where in addition to depth and two spatial coordinates, there are two axes corresponding to slopes, either parameterized as slope in the two spatial directions (x and y) or as slope and azimuth (or dip and strike in geological terms). In 2D, a truncated reflector is a point at the end of a line, while in 3D, a truncated reflector at a fault is a line at the end of a plane. For the 2D case, the truncation of a reflector at a fault will span a full range of slopes at that point in space and depth/time. In the 3D case, these truncations will appear as a line (or a blurred version of one) in either slope x/slope y or alternatively as broadband in slope centered at a single azimuth perpendicular to the fault plane.
Active contours, or snakes, are computergenerated curves that move within images to find object boundaries. Its 3D version is often known in the literature by the name of deformable models or active surfaces. They are often used in computer vision and image analysis to detect and locate objects, and to describe their shape. They were designed as a means to integrate geometric continuity information, such as connectedness or smoothness, into the image analysis. Thus, an active contour may be a deformable curve whose deformation is governed by a model component ensuring smoothness and connectedness, and a data component that lets the contour attach itself to the object during the search in the data. Therefore, they are well suited to describe objects in data with noise or artifacts. Ideally, the data component of the active contour prevails in parts of the curve where the data quality is good, while the model component prevails where the data are distorted or missing. A main application of active contours is to extract geometric objects in data with a low signal to noise ratio or in regions of missing signal.
Active contours are effective in dealing with the ambiguities arising in deciding when a fault line should be extracted from a faulthighlighted volume, because of their capability to attach themselves to the geometric structure in the data that they model. In 3D, the application of deformable models or active surfaces will enable the exploitation of geometric continuity constraints between fault lines. This could lead to a further optimization of the fault tracking procedure while retaining the ability for user intervention in cases where superior knowledge on fault interpretation is available. The optimal location of a fault can be searched in a neighborhood of an initial guess at the solution.
Published references on active contours (2D) and deformable models/active surfaces (3D) include:
 1) M. Kass, A. Witkin, and D. Terzopolous, “Snakes: Active contour models,” International Journal of Computer Vision 1, 321331 (1987).
 2) T. McInerney and D. Terzopoulos, “Deformable models in medical image analysis: a survey,” Medical Image Analysis 1, 91108 (1996).
 3) C. Xu and J. Prince, “Snakes, shapes, and gradient vector flow,” IEEE Trans. Imag. Proc. 7, 359369 (1998).
 4) J. Montagnat, H. Delingette, and N. Ayache, “A review of deformable surfaces: topology, geometry and deformation,” Image and Vision Computing 19, 10231040 (2001).
 5) J. Malik, S. Belongie, T. Leung, and J. Shi, “Contour and texture analysis for image segmentation,” International Journal of Computer Vision 43, 727 (2001).
The foregoing application is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. All such modifications and variations are intended to be within the scope of the present invention, as defined in the appended claims.