SUPPRESSING NEAR-SURFACE SCATTERED SURFACE WAVES

Abstract

Seismic data of a subterranean region is accessed and an image orientation of each data sample of the seismic data is computed. The seismic data is divided, in a spatial domain, into a plurality of sections based on spatial attributes of the seismic data, each section including respective seismic data samples having the spatial attributes within the respective section. For each section of the plurality of sections, a dominant orientation of the section is determined and a noise component of the seismic data is predicted based on the dominant orientation of the section.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of priority to U.S. Provisional Application Ser. No. 62/155,053, filed on Apr. 30, 2015, and U.S. Provisional Application Ser. No. 62/155,207, filed on Apr. 30, 2015. This application is also related to our reference, Attorney Docket Number 38136-0071001, U.S. patent application Ser. No. ______, filed on Monday, May 2, 2016. The entire contents of U.S. Provisional Application Ser. No. 62/155,053, U.S. Provisional Application Ser. No. 62/155,207, and Attorney Docket Number 38136-0071001, U.S. patent application Ser. No. ______ are each hereby incorporated by reference in their entirety.

TECHNICAL FIELD

This disclosure relates to exploration seismology and, more specifically, to seismic data processing.

BACKGROUND

Certain types of seismic waves, such as direct and scattered surface (Rayleigh) waves, contain a significant amount of seismic energy. They present great challenges in land seismic data acquisition and processing, especially in arid regions with complex near-surface heterogeneities (for example, dry river beds, wadis/large escarpments, caves, carbonate rocks, karst features, and any abrupt change in the elastic properties of the earth material (for example, velocity and density) not only due to lithology changes but also can be due to other material property variations such as pressure and fluid saturations).

The near-surface scattered body-to-surface waves have comparable amplitudes to reflections and can mask the seismic reflections. Combining the near-surface scattered body-to-surface wave with large amplitude direct and back-scattered surface (Rayleigh) waves, creates a significant reduction in a signal-to-noise ratio and degrades image quality of a generated final subsurface image that is contaminated with near-surface scattered surface waves.

SUMMARY

This disclosure relates to predicting and reducing noise from seismic data, for example, for reservoir analysis of a subterranean region.

In an implementation, seismic data of a subterranean region is accessed and an image orientation of each data sample of the seismic data is computed. The seismic data is divided, in a spatial domain, into a plurality of sections based on spatial attributes of the seismic data, each section including respective seismic data samples having the spatial attributes within the respective section. For each section of the plurality of sections, a dominant orientation of the section is determined and a noise component of the seismic data is predicted based on the dominant orientation of the section

Particular implementations of described methods and systems can include corresponding computer systems, apparatuses, or computer programs (or a combination of computer systems, apparatuses, and computer program) recorded on one or more computer storage devices, each configured to perform the actions of the methods. A system of one or more computers can be configured to perform particular operations or actions by virtue of having software, firmware, hardware, or a combination of software, firmware, or hardware installed on the system that, in operation, causes the system to perform the actions. One or more computer programs can be configured to perform particular operations or actions by virtue of including instructions that, when executed by data processing apparatus, cause the apparatus to perform the actions.

The above-described implementation is implementable using a computer-implemented method; a non-transitory, computer-readable medium storing computer-readable instructions to perform the computer-implemented method; and a computer system comprising a computer memory interoperably coupled with a hardware processor configured to perform the computer-implemented method/the instructions stored on the non-transitory, computer-readable medium.

The subject matter described in this specification can be implemented in particular implementations so as to realize one or more of the following advantages. First, conventional filtering algorithms for locally linear coherent noise removal are either global (for example, remove a range of slopes from the entire data set) or local (for example, remove a range of slopes within local time-space windows) and cannot precisely predict and separate noise components in seismic data traces based on space-variant slopes (for example, near-surface scattered surface wave noise). Using the described subject matter, the separation and removal of locally linear coherent noise features in seismic data can be improved, leading to superior subsurface images and interpretation results. The proposed method improves noise attenuation and avoids unnecessary signal distortions resulting from using imprecise conventional global velocity filtering approaches. Second, the new capabilities can be used not only to predict and suppress noise features from seismic data, but also to estimate a shear wave velocity of upper-most near-surface layers, which can be utilized as an initial model for velocity estimation schemes. Third, the method can also be applied to ocean-bottom-cable (OBC) marine data, in which the sea bottom acts in a similar way to the free surface, and therefore, propagating waves can scatter to Scholte waves by irregularities or heterogeneities near the sea bottom. Fourth, by predicting a true slope of noise at each offset location, the proposed filter will avoid attenuating part of the reflected signal (for example, reflections from dipping layers) with similar trace-to-trace slope as the noise component. Fifth, applying a directional filter (for example, a non-linear median filter) in the time domain will minimize signal distortions and smearing due to leakage in the transform domain. Although the directional filter is applied in the time domain, the approach can also be generalized to directional filters in other domains, using the predicted dominant slope of the noise as a guide in suppressing the scattered noise. Other advantages will be apparent to those of ordinary skill in the art.

The details of one or more implementations of the subject matter of this specification are set forth in the accompanying drawings and the description below. Other features, aspects, and advantages of the subject matter will become apparent from the description, the drawings, and the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Patent and Trademark Office upon request and payment of the necessary fee.

FIG. 1A is a diagram showing an example earth model where seismic energy is trapped and scattered in the near-surface layer; FIG. 1B is a diagram showing an example ideal model after removing the effects of surface-wave scattering in the example earth model, according to an implementation.

FIG. 2A is a plot showing an example basis filter, an oriented filter with orientation θ=0°; FIG. 2B is a plot showing an example of another basis filter, an oriented filter with orientation θ=90°; and FIG. 2C is a plot showing an example oriented filter with orientation θ=80°, according to an implementation.

FIG. 3A is a plot showing an example seismic image filtered at orientation θ=0°; FIG. 3B is a plot showing an example seismic image filtered at orientation θ=90°; and FIG. 3C is a plot showing an example seismic image filtered at orientation θ=80°, according to an implementation.

FIG. 4 is a diagram showing an example synthetic 2-dimensional (2D) earth model with multiple dipping layers and five scatterers embedded in the uppermost layer, according to an implementation.

FIGS. 5A-5F show example scattering effects for the earth model in FIG. 4 based on finite difference simulations, according to an implementation.

FIG. 6 is a flowchart showing an example process for predicting and reducing noise from seismic data by suppressing near-surface scattered body-to-surface waves, according to an implementation.

FIG. 7 is a diagram showing a frequency-wavenumber domain, according to an implementation.

FIG. 8A shows histograms of local slopes calculated using steerable filters for three patches of offsets from one shot gather shown in FIG. 8B, according to an implementation.

FIG. 9A is a plot showing the input image; FIG. 9B is a plot showing the filtered image; and FIG. 9C is a plot showing the removed noise, according to an implementation.

FIG. 10A is a plot showing an example filtered image using a directional filter, according to an implementation.

FIG. 10B is a plot showing an example filtered image using the f-k filter, according to an implementation.

FIG. 11A is a plot showing an example difference between the input data (with noise) and the de-noised results using a directional filter; FIG. 11B is a plot showing an example difference using f-k filter, according to an implementation.

FIGS. 12A-C show example results of applying the example steered median filter approach to field data, according to an implementation.

FIG. 13A shows histograms of local slopes calculated using steerable filters, corresponding to patches shown in FIG. 13B, respectively, according to an implementation.

FIGS. 14A-C show input data, filtered data, and the residual of the application of the steered median filter approach to field data after a Normal Move Out (NMO), running average filter, and inverse NMO to enhance the reflections, according to an implementation.

FIGS. 15A-C illustrates input field data, filtered field data, and a residual (the difference or noise removed), respectively, of the application of a f-k filter to field data after NMO, running average filter, and inverse NMO to enhance reflections as illustrated in FIGS. 15A-B, according to an implementation.

FIGS. 16A-C are plots showing three-dimensional (3D) finite difference simulation results when the scatterers are in-line with the receivers, according to an implementation.

FIGS. 17A-C are plots showing 3D finite difference simulation results when the scatterers are in the cross-line direction relative to the receivers, according to an implementation.

FIGS. 18A-D are plots showing aspects of a 3D irregular (Gaussian) bedrock interface model, according to an implementation.

FIGS. 19A-C are plots showing aspects of finite difference results for the irregular (Gaussian) bedrock interface model simulated using 2D finite difference techniques, according to an implementation.

FIGS. 19D-F are plots showing aspects of finite difference results for the irregular (Gaussian) bedrock interface model simulated using 3D finite difference techniques, according to an implementation.

FIG. 20 is a plot showing an example 3D earth model with multiple dipping layers and near-surface scatterers.

FIGS. 21A-B are plots showing aspects of modeling and filtering of scattered body-to-surface waves in 3D, according to an implementation.

FIGS. 22A-B are plots showing aspects of application of a 3D f-k filter to spatially dense sampled 3D simulated data, according to an implementation.

FIG. 23 is a block diagram illustrating an example computer used to provide computational functionalities associated with described algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure, according to an implementation.

Like reference numbers and designations in the various drawings indicate like elements.

DETAILED DESCRIPTION

This disclosure describes computer-implemented methods, software, and systems for suppressing near-surface scattered surface waves based on local slope estimation and separation of scattered surface waves from reflected body waves in seismic data.

For purposes of this disclosure, “near-surface” refers to the shallowest one wavelength depth from the earth's surface. For example, in case of 30 Hz dominant frequency and 1800 m/s wave speed, the shallowest one wavelength is equivalent to the top 60 m depth from the Earth's surface. A “P-wave” denotes a primary or compressional wave, an “S-wave” denotes a secondary or shear wave, and “surface waves” refer to Rayleigh waves that are confined with propagation along the earth's free surface.

The ultimate goal in exploration seismology is to obtain an optimum image of the subsurface to map reservoir boundaries (for example, structural mapping) and to characterize reservoir properties (for example, porosity, saturation, permeability, and fracture density and orientation) for field evaluation and development. Extracting and monitoring reservoir properties are essential steps for reservoir modeling, flow simulation studies and for horizontal drilling to accomplish optimal recovery. These processes can be achieved by inverting formation elastic properties (for example, P-impedance, S-impedance, and density) from seismic data and by using rock physics information derived from core and log data to link between the rock and seismic elastic properties. However, in the real world, many factors affect the quality of seismic data. Surface topography and near-surface heterogeneities can significantly complicate the data with large static travel time variations, signal attenuation, and wave scattering. These factors can deteriorate the image quality and increase the uncertainty of subsurface images and formation properties.

Near-surface complexity and surface topography (for example, in arid environments) present major challenges to land seismic data acquisition and processing. These challenges can severely affect data quality and introduce uncertainty into reservoir imaging and characterization as the near-surface complexity and surface topography can severely contaminate the seismic data with complicated wave phenomena that cannot be accounted for by single layered models and conventional processing methods. Both the direct surface waves and the up-going reflections are scattered by the heterogeneities. The scattering takes place from body waves to surface waves, and from surface waves to body waves. Further, the near-surface scattered body-to-surface waves, which have comparable amplitudes to reflections, can mask the seismic reflections. These difficulties, added to large amplitude direct and back-scattered surface (Rayleigh) waves, create a major reduction in signal-to-noise ratio and degrade the final subsurface image quality. Also, strong contrast in velocity (and therefore impedance) between the near-surface and the substrata (for example, carbonate layers) may add to the complexity of the energy penetration, as most of the seismic energy is reflected and trapped near the surface, generating both internal and surface related multiples. In interpreting noise-contaminated data, the main challenges are caused by complex topographic and near-surface features such as sand dunes, wadis/large escarpments, karsted carbonates, caves, dry river beds, and any abrupt change in the elastic properties of earth material (for example, velocity and density) not only due to lithology changes but also can be due to other material property variations such as pressure and fluid saturations. In such environments, it is essential to understand the various types of near-surface effects such as direct waves, mode-conversions, scattering, and attenuation.

The scattered elastic waveforms have strong effects on both the phase and amplitude of the recorded signal and are neglected in most conventional imaging and interpretation schemes. A scattering model (for example, a Born approximation) assumes a weak single scattering and an incident wave propagating in a homogeneous or smooth background media. The Born approximation fails in the case of strong scatterers and multiple scattering waves. More widely used, the acoustic approximation assumes only pressure waves are involved and neglects the effects of P to S mode-conversion. The acoustic approximation fails when upcoming body waves scatter from near-surface heterogeneities and generate scattered surface waves. These assumptions can lead to significant deviations in amplitude and phase of the scattered wavefield, due to strong impedance and velocity anomalies. Re-datuming techniques, based on wave-equation migration, common-focus-point technology, and interferometry, can handle dynamic time shift (that is, statics) of the P wave reflected energy (acoustic response) due to delays caused by near-surface velocity variations. However, these techniques cannot handle multiple scattering waves in fully elastic media even when an exact sub-surface velocity model is available. Therefore, the Born and acoustic assumptions can severely affect the interpretation and imaging techniques that require high quality seismic attribute maps and are sensitive to the level of signal-to-noise ratio, such as velocity estimation and migration, elastic and anisotropic parameter estimations for Amplitude Variations with Offset (AVO) and Amplitude Variations with Offset and Azimuth (AVOz) studies, and 4D time-lapse for reservoir monitoring.

Understanding and alleviating the effects of near-surface complexities on the seismic records can improve the signal-to-noise ratio, mainly reducing the noise components due to scattered body-to-surface waves. A noise component is modeled using numerical analysis in order to explain observations made in field data and to understand noise component characteristics to develop an effective algorithm for its subsequent removal. Building on the analysis of the numerical modeling results, a technique is described to reduce scattered body-to-surface waves based on spatially varying slope prediction and separation, using steerable and non-linear median filters. Suppression of these scattered surface waves can be difficult using conventional filtering methods, such as an f-k filter, without distorting the reflected signal.

Modeling and studying the contribution of near-surface heterogeneities on seismic wavefield scattering can provide better understanding of land seismic data, provide an effective approach to filter out the scattered noise from the seismic records to enhance the signal-to-noise ratio, and to accurately image and locate the near-surface heterogeneities. This can be accomplished by techniques to: 1) spatially vary local-slope estimation, aimed at alleviating the effects of scattered surface waves and enhancing the quality of the seismic signal and 2) simulating the effects of seismic wave scattering from buried, shallow, subsurface heterogeneities through finite difference numerical forward modeling.

In spatially vary local-slope estimation, the technique is based on predicting a spatially varying slope of the noise, using steerable filters, and separating the signal and noise components by applying a directional non-linear filter oriented toward the noise direction to predict the noise and then subtract it from the data. The slope estimation step using steerable filters is very efficient, as it requires only a linear combination of a set of basis filters at fixed orientation to synthesize an image filtered at an arbitrary orientation. We apply our filtering approach to simulated data as well as to seismic data recorded in the field to suppress the scattered surface waves from reflected body-waves, and demonstrate its superiority over conventional f-k techniques in signal preservation and noise suppression.

While this disclosure emphasizes two-dimensional (2D) earth models and presents three-dimensional (3D) examples when it is essential to demonstrate the differences between the two modes (as 2D models are generally more feasible than 3D earth models in terms of the computational efficiency), this disclosure is not limited to only 2D or 3D models. As will be appreciated by those of ordinary skill in the art, the described techniques are applicable with appropriate modifications to other types of earth models, other dimensionality, and similar Seismic energy recorded in the field represents a true earth response that includes all complicated wave types (for example, multiple scattering, elastic mode conversion, and viscoelastic attenuation). In land seismic data, near-surface complexity not only causes time shifts as assumed by static corrections, but also causes amplitude distortions to the recorded signal. A significant fraction of the seismic energy is trapped and scattered in the near-surface layers (in the form of coherent noise) and masks the body-wave reflections from deeper structures. It is desirable to separate the scattered surface (Rayleigh) wave energy and remove the effects of near-surface heterogeneities.

In some instances, upcoming body-wave reflections (for example, primaries, multiples, and mode conversions) are scattered close to the free surface due to near-surface heterogeneities, and hence excite scattered surface waves. In contrast to the direct Rayleigh wave, the scattered surface-wave noise can cover the whole time-space domain with locally linear (also known as piece-wise linear) features. For example, the scattered surface-wave noise can include plane waves propagating along the earth's free surface that propagate in both positive and negative (that is, forward and backward) directions. The scattered surface waves are more sensitive to the uppermost earth layer properties. Therefore, their slope (that is, orientation) can vary due to lateral velocity variations (for example, mostly shear-wave velocity). The slope or image orientation of a seismic wave, also referred to as the wave apparent velocity, is defined at each data sample in the time-offset domain as the rate of change of recorded offset or distance over time (for example, gradient) calculated from immediate neighboring traces (that is, receiver locations). At each receiver location, however, direct and scattered surface waves can have a constant slope depending on the shear wave velocity, independent of the source depth and location. In some implementations, the slope of the surface Rayleigh wave noise can be assumed to be locally linear and correspond to a velocity of about 0.92 times the shear velocity. The local slope can be determined for each data sample (also referred to as a data point or data record).

Because the local slope of the noise only varies spatially with different receiver location (or patch of receivers), unlike global velocity filters, only a narrow range of slopes may be filtered at each offset. The offset is a spatial attribute and refers to the distance from the seismic source to the seismic receiver in surface seismic acquisition. As such, the seismic data can be divided into multiple sections or windows in the spatial domain, for example, as a function of offset. In some implementations, the example techniques can suppress locally linear scattered surface waves by estimating the dominant local slope in each spatial window, predicting the noise component based on the dominant local slope, and removing the noise component from the seismic data.

In some implementations, the image orientation of each data sample can be predicted, for example, using steerable filters or other techniques. For example, a steerable filter oriented at a particular orientation can be implemented by a linear combination of a set of basis filters at fixed orientation. A bank of orientation filters that span all possible seismic data orientations can be constructed for the estimation of image orientation at each data sample. Image orientations at each data sample can be predicted using the bank of orientation filters spanning all possible supplementary orientations (for example, 0°/180°, 1°/179°, 2°/178°, etc.), by finding the orientation with minimum energy at each data sample.

The spatially variant slopes can be determined by finding the orientation that maximizes the mass density function of image orientations within laterally sectioned windows. The orientation that maximizes the mass density function can be determined, for example, using histograms of predicted image orientations, and assigned to the laterally sectioned windows. Such an orientation can be referred to as the dominant image orientation or slope and be regarded as the image orientation for the noise component (also referred to as noise orientation) for the respective window.

In other words, the bank of orientations can be used to predict the image orientation at each data sample and the image orientations of all data samples within each lateral section can be used to predict the dominant slope within each respective section.

Based on the determined dominant local slopes for the multiple spatial windows, the noise component can be predicted or enhanced, for example, by using a directional filter (for example, a non-linear median filter) or other techniques. In some implementations, a directional filter (for example, based on dip, slope, or velocity) oriented towards the noise orientation can be applied to predict the noise component. Then the noise component can be subtracted from the data. The example techniques can exploit the multidimensionality of the seismic data (for example, temporal, spatial, and directional variables) to obtain more reliable signal and noise separation.

In some implementations, the directional filter is narrow in the f-k domain (for example, tackling a specific slope instead of a range of slopes), which ensures minimal signal distortion and a focused subsurface image. The application of the filter operator on the total data will predict only the noise component; the signal component can be obtained by subtracting the predicted noise from the input data to produce the filtered record.

The advantages of the example techniques include improvement of the separation and removal of locally linear coherent noise components in seismic data, better signal preservation and noise suppression, increase of the signal-to-noise ratio of the seismic data, enhancement of the quality and accuracy of seismic data processing, and superior subsurface images and interpretation results. The example techniques can improve noise attenuation and avoid unnecessary signal distortions resulting from using imprecise conventional global velocity-filtering approaches. In some applications, the example techniques can achieve additional or different advantages.

The example techniques can be used not only to predict and suppress noise features from seismic data, but also to estimate the shear wave velocity of the uppermost near-surface layers (as the spatially varying slopes can also relate to the shear wave velocity of the uppermost near-surface layer), which can be utilized as an initial model for velocity estimation schemes. The example techniques can be applied to simulated data, seismic data recorded in the field, or other data to suppress the scattered surface waves from reflected body waves. The example techniques can also be applied to ocean-bottom-cable (OBC) marine data, in which the sea bottom acts in a similar way to the free surface; therefore, propagating waves can scatter to Scholte waves by irregularities or heterogeneities near the sea bottom.

FIG. 1A is a diagram showing an example earth model 100 where seismic energy is trapped and scattered in the near-surface layer. The example earth model 100 includes a seismic source 110 (for example, a seismic air-gun or vibroseis) that can emit energy (for example, air-gun energy) beneath or on the free surface 115. Acoustic waves (for example, body waves 140 generated by the source 110) can travel through the earth 116, be reflected by seismic reflectors (for example, layer boundaries 125, 135), and be recorded by borehole receivers 120 (for example, geophones) located above free surface 115 or any other place. In some implementations, besides the reflected body waves 145, there can be coherent noise 155, for example, induced by the “scattered reflections” 160 caused by scatterers 165. Specifically, the earth model 100 shows scattered surface (Rayleigh) waves 155 excited by the direct surface (Rayleigh) waves 130, direct body waves 140, and upcoming body wave reflections 145. FIG. 1B is a diagram showing an example ideal model 105 after removing the effects of surface-wave scattering in the example earth model 100, according to an implementation.

In some implementations, a computer system 150 can be used to acquire seismic data. The computer system 150 can include one or more computing devices or systems. The computer system 150 or any of its components can be located apart from the other components shown in FIG. 1A. For example, the computer system 150 can be located at a data processing center, a computing facility, or another suitable location. The earth model 100 can include additional or different features, and the features of the earth model 100 can be arranged as shown in FIG. 1A or in another configuration.

The acquired seismic data can include both signals and noises. In the case of a surface seismic source, the main components of surface-related noise include: (1) direct-surface waves; (2) forward- and back-scattered surface waves; and (3) body-to-surface scattered waves. The upcoming body-wave reflections (for example, P and S waves, including primaries, multiples and mode conversions) impinge on the near-surface heterogeneities and scatter to weak P and S waves, and also to Rayleigh waves since the heterogeneities act as secondary elastic sources. The scattered surface-to-body and body-to-body wave components, however, pose less challenges as they are much weaker in amplitude and attenuate much faster with propagation distance than surface waves.

Surface waves, in general, are characterized by low frequency, linear move-out, large amplitudes, and slower amplitude decay with distance. The direct-surface wave is confined within a fan-shaped window in the time-space domain and is much larger in amplitude and lower in frequency than body-wave reflections. These characteristics of direct-surface waves make it effective to apply filtering techniques within a narrow fan-shaped window. In contrast, the scattered body-to-surface waves have comparable amplitudes to reflections (especially in the cases of large, high-contrast and shallow scatterers), and they could contaminate the entire dataset.

In some implementations, source-generated noise can be filtered and attenuated based on its characteristics (for example, frequency, velocity, amplitude, and polarization). In some implementations, example filtering techniques can exploit the multidimensionality of the seismic data (for example, temporal, spatial, and directional variables) to obtain more reliable signal and noise separation. The example filtering techniques can separate locally linear scattered surface waves by first estimating the dominant local slope as a function of offset using steerable filters and then applying a non-linear median filter oriented in the noise direction. The example techniques can suppress and reduce the noise by subtracting the predicted noise from the data to obtain the signal.

Steerable Filters

In some implementations, oriented filters can be used in many vision and image processing tasks such as edge detection, texture analysis, segmentation, motion analysis, and image enhancement and compression. The steerable filter is one of a class of filters in which a linear combination of a set of basis filters at a fixed orientation is used to synthesize a filter at an arbitrary orientation. It can be used to examine filter outputs (images) as a function of phase and orientation, efficiently, without the need to apply many versions of the same filter rotated at different angles. An example illustration of the steerable filters is discussed in the following using the partial derivative of a 2D symmetric Gaussian function in Cartesian coordinates x and y:

G(x,y)=e^−(x2+y2).  (1)

The first x derivative of a Gaussian is:

G^0°=−2xe^−(x2+y2)  (2)

and they derivative is rotated 90° as:

G^90°=−2ye^−(x2+y2).  (3)

Therefore, a filter at an arbitrary orientation θ can be synthesized by taking a linear combination of x and y derivatives (basis filters):

G^θ°=cos(θ)G^0°+sin(θ)G^90°.  (4)

FIGS. 2A-2C show example derivatives of Gaussian filters, according to an implementation. FIG. 2A is a plot showing an example basis filter 201, an oriented filter with orientation θ=0°; FIG. 2B is a plot showing an example of another basis filter 202, an oriented filter with orientation θ=90°; and FIG. 2C is a plot showing an example oriented filter 203 with orientation θ=80°, which is obtained by linearly combining the basis filters 201 and 202.

In some implementations, a seismic image filtered at an arbitrary orientation R (x, y, θ) can be synthesized by convolving the input image I(x, y) with the orientation filter G^θ°:

R(x,y,θ)=cos(θ)(G^0°*I(x,y))+sin(θ)(G^90°*I(x,y)),  (5)

where,

• • * denotes convolution.

FIGS. 3A-3C show example convolution of the input image with different orientation filters, according to an implementation. FIG. 3A is a plot showing an example seismic image 301 filtered at orientation θ=0° (the convolution of an input seismic image with the basis filter oriented at 0°); FIG. 3B is a plot showing an example seismic image 302 filtered at orientation θ=90° (the convolution of the input seismic image with the basis filter oriented at 90°); and FIG. 3C is a plot showing an example seismic image 303 filtered at orientation θ=80°, which is obtained by linearly combining the convolution of the input seismic image with the basis filters, that is, linearly combing the filtered seismic images 301 and 302.

Slope estimation of local plane-waves

In some implementations, the local plane-wave partial differential equation can be represented by

$\begin{array}{cc}\left(\frac{\partial }{\mathrm{dx}}+s\left(t,x\right)\frac{\partial }{\mathrm{dt}}\right)u\left(t,x\right)=0,& \left(6\right)\end{array}$

where,

• • u denotes the seismic wavefield and s denotes the local slope, which can vary in both time t and space x. The time-space derivative of the plane-wave equation is equivalent to the convolution operator G(θ) applied to the data u (seismic wavefields):

G(θ)u=0.  (7)

The orientation θ in equation (7) is related to the local slope s in equation (6). The convolution operator can be efficiently constructed using steerable filters. When the steerable filter orientation is parallel to the feature orientation in the image, the output amplitude is maxima. In contrast, when the steerable filter orientation is perpendicular to the feature orientation in the image, the output is zero. In some implementations, two oriented filters with supplementary orientations (for example, θ and 180°−θ, where θ=1, 2, . . . , 90°) can be applied in sequence to the image. The local slope (that is, local-image orientation) can be determined at each data sample by spanning all possible orientations of supplementary angles (for example, 0°/180°, 1°/179°, 2°/178°, etc.) to find the orientation corresponding to the minimum amplitude of the steerable filter outputs, as described by minimizing the left hand side of equation (7). Alternatively, the local slope can be determined by identifying the maximum amplitude but with rotating the orientation/angle by 90°. The range of orientations corresponds to different apparent velocities of both the signal and noise. Supplementary orientations can be used to take into account that surface waves propagate at the same velocity (that is, image orientation) not only in the forward (that is, positive) but also the backward (that is, negative) direction. An instantaneous slowness (local slope) plot can be constructed.

The Gaussian filter in the steerable filter definition plays an important role as a regularizing (smoothing) term to avoid unwanted, oscillatory, instantaneous slowness estimates. The smoothing of this operator is determined by the variance of the Gaussian operator.

In some implementations, an instantaneous slowness (local slope) table can be constructed using steerable filters, for example, by computing the local image orientations at each time-space data sample. The slowness table includes the local slope of some or all seismic data records, for example, calculated according to the techniques described with respect to equation (7) or another method. In some implementations, the local slope of some or all seismic data records can be stored in another data structure.

In some implementations, the table can be sectioned laterally in the spatial domain across the data records as a function of offset. A modified instantaneous slowness table can be obtained by assignment of each section with a single constant slope or image orientation value (to each data sample in the section) that has maximum probability within each time-offset range:

$\begin{array}{cc}{\hat{\theta }}_n\left(x\right)=\arg \underset{\theta }{\max }P\left(\theta ,x\right),\theta =0,1,\dots ,90°,& \left(8\right)\end{array}$

where,

• • {circumflex over (θ)}_n(x) is the dominant local slope as a function of offset x; θ refers to a range of supplementary orientations (that is, θ and 180°−θ, where θ=1, 2, . . . , 90°) corresponding to different apparent velocities, and P(θ,x) is the probability as a function of orientation θ and offset x. The {circumflex over (θ)}_n(x) can be regarded as the image orientation of the noise component.

Signal and Noise Separation

The input data d can include unknown signal component d_s and noise component d_n:

d=d_s+d_n.  (9)

The modified instantaneous slowness table with a constant slope value for each spatial window can be used to steer the directional filter (for example, based on dip, slope, or velocity) or the non-linear median filter toward the noise direction. For example, for each temporal data point within a spatial window, a non-linear local median filter steered towards the predicted noise direction {circumflex over (θ)}_n can be applied. The application of the non-linear local median filter can enhance (predict) the data component (that is, noise) that is aligned with the orientation {circumflex over (θ)}_n and attenuate all other components (that is, signal) with different orientations:

F({circumflex over (θ)}_n)ds=0 and

F({circumflex over (θ)}_n)d_n=d_n,  (10)

where,

• • F ({circumflex over (θ)}_n) is the filter operator steered toward the predicted noise orientation {circumflex over (θ)}_n, and d_s and d_n are the signal and noise components, respectively.

As an example, a 2M-point steered non-linear median filter can be applied to the data:

d[m,n]=  (11)

where,

• • −M≦i≦i+M; −M·tan({circumflex over (θ)}_n)≦j≦j+m·tan({circumflex over (θ)}_n); d[m,n] is the data point at the m trace and n time sample; i and j are the sample indices in the spatial (x) and temporal (t) directions, respectively; and {circumflex over (θ)}_n is the dominant local slope (that is, orientation) of the noise as a function of offset. The center and median values of the time-space window can be very similar and represent the amplitude value of the estimated noise component when there is no reflection (signal) information (F({circumflex over (θ)}_n)d_n=d_n). Subsequently, if there are reflections, the value of the center sample will be replaced by the median value of its neighboring points (F({circumflex over (θ)}_n)d_s=0).

Therefore, the application of the median filter on the total data will predict only the noise component:

F({circumflex over (θ)}_n)d=F({circumflex over (θ)}_n)(d_s+d_n)=F({circumflex over (θ)}_n)d_n+F({circumflex over (θ)}_n)d_s=d_n.  (12)

The signal component can be obtained by subtracting the predicted noise (d_n) from the input data (d) to produce the filtered record (d_s):

d_s=d−d_n=d−F({circumflex over (θ)}_n)d.  (13)

In some implementations, additional or different filters or techniques can be used to predict the noise component based on the predicted dominant orientation for each spatial window, and to suppress the noise from the data.

Synthetic Example

The example filter techniques of noise separation based on spatially varying slopes are discussed with a numerical example. In some implementations, the modeling and understanding of the complicated scattered wave features play an important part in the success of their subsequent removal. In some instances, seismic forward modeling based on finite difference methods can handle complex lateral variations in material properties and produce a complete and accurate solution to the elastic wave equation, with all direct, converted, scattered, and guided waves.

Finite Difference Modeling

To fully model elastic waves in the presence of heterogeneity, an accurate implementation of the standard staggered-grid finite difference scheme can be used with a convolution-matched, layer-absorbing boundary condition. The staggered-grid scheme is fourth-order accurate in space (including the free surface boundary) and second-order accurate in time. The internal interfaces can be represented by the so-called effective medium parameters to avoid spurious numerical diffractions caused by material discontinuities and discretization. The density is calculated by arithmetic average, and Lame parameters are calculated by harmonic average.

FIG. 4 is a plot showing an example synthetic 2D earth model 400, according to an implementation. Multiple dipping layers 410, 420, 430, 440, and 450 with five circular scatterers 460a-e (circles near the free surface) are embedded in the shallow layer 410. The scatterers 460a-e are located at 15 m depth below the free surface, and each has a 10 m diameter and an impedance contrast of 0.36. The source (not shown) is located at (x, z)=(150 m, 0 m). The receivers (not shown) are located on the surface with 50 m near-offset and 5 m space intervals. The color scale (on the right) and associated numbers refer to material properties given in Table 1:

TABLE 1
Material properties (P wave velocity, S wave velocity,
and density) of the model shown in FIG. 4.
	Material Index	Vp (m/s)	Vs (m/s)	Density (kg/m3)
	1 - Dark blue	1800	1000	1750
	2 - Blue	2200	1200	1900
	3 - Green	2500	1300	2000
	4 - Orange	2700	1400	2100
	5 - Brown	3000	1500	2250

The domain has N_x=1001 and N_z=501 grid points with 1 m grid spacing (that is, Δx and Δz), that is, 500 m depth (along the z-axis) and 1000 m distance (along the x-axis). The grid size is small enough to capture the shape of the scatterers. The time step is 0.2 ms. A vertical source is used with a Ricker wavelet and 30 Hz central frequency (˜75 Hz maximum frequency). In this example, only the vertical component (v_z) of the particle velocity field is considered. The scatterers 460a-e are treated in the numerical scheme as density and velocity perturbations. Calculated waveforms for the scattering model with and without the direct-surface waves are shown in FIGS. 5A-F.

FIGS. 5A-5F show example scattering effects for the earth model 400 in FIG. 4 based on finite difference simulations (v_z-component), according to an implementation. FIGS. 5A-C are plots 510, 520, and 530 showing free surface and Rayleigh waves generated by the source, and FIGS. 5D-F are plots 540, 550, and 560 showing the reflections without the source-generated Rayleigh waves, respectively. FIGS. 5A and 5D show the total wavefield simulated using the earth model 400 with scatterers; FIGS. 5B and 5E show the reference wavefield simulated using the earth model 400 without scatterers; and FIGS. 5C and 5F show the scattered wavefield (that is, the difference between the total and reference wavefields).

The direct and back-scattered surface waves are removed as they are much larger in amplitude than the scattered body-to-surface waves. This can be achieved by computing the wavefield for a homogeneous full-space, with and without the near-surface scatterers, and then subtracting the direct-surface waves from the incident and total wavefields, respectively. The effective medium properties of the heterogeneities are wavelength-dependent. Therefore, the near-surface scatterers, though relatively small compared to the incident wavelength, produce high-frequency features. In this example, scatterers with a size of ⅙ of the wavelength (10 m) produce significant scattering.

FIG. 6 is a flowchart showing example process 600 for predicting and reducing noise from seismic data by suppressing near-surface scattered surface waves, according to an implementation. The process 600 can be implemented, for example, as computer instructions stored on computer-readable media and executable by data processing apparatus (for example, one or more processor(s) of the computer system 150 in FIG. 1A). In some implementations, some or all of the operations of process 600 can be distributed to be executed by a cluster of computing nodes, in sequence or in parallel, to improve efficiency. The example process 600, individual operations of the process 600, or groups of operations may be iterated or performed simultaneously (for example, using multiple threads). In some cases, the example process 600 may include the same, additional, fewer, or different operations performed in the same or a different order.

At 610, seismic data of a subterranean region can be received or otherwise accessed by the data processing apparatus (for example, one or more processor(s) of the computer system 150 in FIG. 1A). The seismic data can be seismic data collected, for example, by receivers (for example, the receivers 120 in FIGS. 1A and 1B) and transmitted to the data processing apparatus. The seismic data could be sorted as common shot gathers (that is, seismic traces corresponding to one source but multiple receiver locations) or common receiver gathers (that is, seismic traces corresponding to one receiver but multiple source locations). In some implementations, the seismic data can be stored in a computer-readable media (for example, memory) and the processing apparatus can load the seismic data from the computer-readable media.

The seismic data can include a number of data samples or points. In some implementations, the seismic data can include, for example, input data traces D (S, R, T) with a source location S, receiver location R, and the travel time T for multiple data points. In some implementations, the data points can be represented, for example, as seismic wavefields in a time-space domain such as u (t, x) in equation (6). The wavefield can represent the amplitude of the seismic signal as a function of time (t) and offset (x). The data points can be represented in other domains (for example, time-frequency, frequency-distance, or other domains) or formats. The data points can be converted from one domain to another by appropriate transformations. The seismic data can include signal and noise components. The seismic data can be field data, synthetic data, or a combination of these and other types of data. The seismic data can include 2D, 3D, or higher-dimensional data.

At 620, dip decomposition can be applied to the seismic data, for example, by using a velocity filter (for example, in the f-k domain) to reduce the interference between signal and noise (for example, mainly by reducing the signal component) to narrow the search space of image orientations such that the slope estimation yields a more reliable result and can narrow the search space of local-image orientations. As such, only a narrow range of orientations is searched (for example, passed by the dip decomposition). For example, dip decomposition using a frequency-wavenumber fan filter reduces the signal and noise interference for enhanced slope estimations of the noise component. The process can be performed by passing dip components in the data that are within a close range of the true noise dip or slope, and rejecting (that is, filtering out) all other dip components as guided by the fan-filter:

d(x,t)=F⁻²[V(k,w)F²[d(x,t)]],  (14)

where,

• • d denotes the dip decomposed part of the data, V corresponds to the fan filter in the spectral domain, and F² and F⁻² denote the forward and inverse 2D Fourier transforms, respectively.

At 630, an image orientation of each data sample of the seismic data is computed. For example, local image orientations (that is, local slopes) can be computed and predicted for one or more shot gathers (for example, source locations such as the source 110 in FIGS. 1A and 1B). The image orientation of each data sample can be computed using a number of oriented filters with supplementary orientations (for example, θ and 180°−θ, where θ=1, 2, . . . , 90°). In some implementations, the image orientation of each data sample can be computed according to the example techniques described with respect to equation (7) or another method. For example, steerable filters G(θ) and G(180°−θ) can be applied successively to a data sample u and its immediate neighboring samples (as the operator of the steerable filter G can be of different size (for example, 3 by 3 matrix)). For example, the data samples can include raw data samples or dip-decomposed data sample as a result of operation 620 according to the directional derivatives equation (7). By searching over all image orientations of the steerable filters G(θ) and G(180°−θ), the image orientation that minimizes the left-hand side of equation (7) (that is, returns minimum amplitude of the steerable filter outputs) can be identified. Such an image orientation can be regarded as the image orientation of the data sample. In some implementations, an instantaneous slowness table can be constructed to record the image orientation and the minimum amplitude of the steerable filter output for each data sample.

In some implementations, in addition to or as an alternative to the steerable filters for computing image orientations, other slope estimation techniques such as slant-stacking or plane-wave destruction filters can be used. In some implementations, steerable filters can be more computationally efficient than others as they can be constructed by a linear combination of a set of basis filters at fixed orientations.

FIG. 7 is a diagram 700 showing a frequency (f)-wavenumber (k) domain, according to an implementation. The solid lines 710a-b show the range of wavenumbers constrained by the f-k filter (for example, as a result of operation 620), and the dashed colored lines 720-780 show the frequency-wavenumbers corresponding to different steerable filter orientations (for example, the directional derivatives computed based on equation (7)).

Referring back to FIG. 6, at 640, the seismic data can be divided, in a spatial domain, into a number of sections based on spatial attributes (for example, offset values) of the seismic data. Each section includes respective seismic data samples having the spatial attributes within the respective section. For example, each spatial section can be represented by a spatial range having a lower endpoint and an upper endpoint as a function of offset (for example, measured in meters, feet, etc.). The sections can be overlapping; or they can be adjacent or disjointed, non-overlapping or with each other. The sections can have the same widths. The number of sections can be a default or configurable value. For example, it can depend on how the shear wave velocity changes laterally in the upper layer. Strong lateral variations may require shorter windows and hence more sections to resolve the variations. FIG. 8B shows that seismic data are divided into three example spatial sections 815, 825, and 835 as a function of offset in meters, according to an implementation. FIG. 13B shows that seismic data are divided into four example spatial sections 1315, 1325, 1335 and 1345 as a function of offset (for example, in meters), according to an implementation.

In some instances, the scattered body-to-surface wave noise can cover the whole time-space domain with locally linear noise features. The local slope of the noise only varies spatially with different receiver location (or patch of receivers) whether the scattered surface waves are generated by upcoming primary, multiple, or mode converted waves, and regardless of the source location. Therefore, the scattered noise from one or more common shot gathers makes a repetitive pattern of constant slopes at each offset, which can be analyzed by computing local image orientations and identifying dominant slopes from statistical histogram distributions. As a result, only a narrow range of slopes may be filtered at each offset instead of filtering-out a fixed range of slopes (for example, global f-k filter) throughout the entire data, giving rise to the advantage that dividing the sections in the spatial domain.

At 650, for each of the plurality of sections, a dominant orientation of the section is determined. The dominant orientation of the section can be determined according to the example techniques described with respect to equation (8) or another method. For example, the dominant orientation of the section can be determined based on probability distributions of respective image orientations of data samples within the section from one or more shot gathers. The dominant orientation can be identified as the orientation that maximizes the mass density function of the orientations of data samples in the section. The dominant orientation can be regarded as the dominant local orientation of the noise component of the section. The predicted dominant orientations or slopes can also represent the shear-wave speed at the upper-most earth layer.

For example, FIG. 8A shows histograms 810, 820 and 830 of local slopes calculated using steerable filters for three patches of receiver offsets from one shot gather, according to an implementation. The top, middle, and bottom histograms 810, 820, and 830 show the orientation distributions of the data samples in sections 815, 825, and 835 shown in FIG. 8B, respectively. The lines 840, 850, and 860 are the true orientation of the forward and backward scattering of the data samples in the sections 815, 825, and 835, respectively. As illustrated, the bins 845, 855, and 865 having the highest values can be identified as the dominant orientations for the three sections 815, 825, and 835, respectively, and they match the true slopes of the noise components 840, 850, and 860 in all cases.

In this case, the histograms 810, 820, and 830 look identical (because f-k values are about the same). The offset sections allow for spatially varying slope estimation, which can be useful in the case of lateral shear wave velocity variations near the earth's surface. Each orientation represents the supplementary angles (θ) and (180°−θ). Similarly, the scattered surface waves propagating at the forward and backward directions are recorded at two supplementary angles (that is, slopes) at each offset location.

Referring back to FIG. 6, at 660, the noise component can be predicted based on the dominant orientation of the section. For example, the noise component can be predicted by applying, to each of the number of sections of the seismic data in 610, a directional filter (for example, a non-linear median filter) guided by the spatially-varying dominant slope. In some implementations, predicting the noise components can also be achieved by not only using a non-linear median filter but also by using directional filters such as local velocity, slope, or dip filters not only in the time domain, but also in other domains (for example, f-k or tau-p). In some instances, non-linear median filters show less noise and artifacts compared to other directional filters.

The operator of the above-described directional filter guided by the spatially-varying dominant slope returns the amplitude of the surface wave when there is no reflection interference, and it replaces the amplitude of the sample with its median of neighboring samples when reflections are present. According to an implementation, the results of separating the signal and noise components are shown in FIGS. 9A-C. Specifically, FIG. 9A is a plot showing the input image 910; FIG. 9B is a plot showing the filtered image 920, filtered according to the example process 600 using the above-described spatially-varying directional filter; and FIG. 9C is a plot showing the difference (that is, the removed noise) 930 between the input image 910 and the filter image 920. The filtered signal (that is, the signal after removing noise) is free of the scattered body-to-surface wave noise, and only very weak scattered P waves with apparent velocity (that is, slope) similar to the reflected signal are passed by the filter.

In some implementations, the data can be divided in spatial sections and transformed into other slope/orientation domains (for example, tau-p slowness-time intercept domain). The dominant slope for each section can be estimated using statistical histogram distribution in the respective transformed domain.

Referring back to FIG. 6, at 670, the noise predicted based on the predicted noise image orientations can be subtracted from the received seismic data, according to the example techniques described with respect to equations (11)-(13) or another method.

At 680, the filtered seismic data can be output, for example, by transmitting over a communication interface, storing in a data store, displaying on a user interface (for example, as shown in FIGS. 9B and 10A), etc. The filtered seismic data can have a greater signal to noise ratio than the received seismic data. The filtered seismic data can be used for other processing, for example, for reservoir analysis.

FIG. 10A is a plot showing an example image 1010 filtered according to the example process 600 using the above-described spatially-varying directional filter, according to an implementation. FIG. 10B is a plot showing an example image 1020 filtered using an f-k filter, according to an implementation. FIG. 10B shows the edge effects 1030 and smearing of the reflected signal caused by the f-k filter due to leakage in the transform domain. Even in the case of local f-k filtering, there will still be some artifacts, which are window-dependent.

FIG. 11A is a plot showing an example difference 1110 between the input data (with noise) and the de-noised results according to the example process 600 using a spatially varying directional filter, according to an implementation. FIG. 11B is a plot showing an example difference 1120 using the f-k filter with edge effects 1130. Note that the f-k filter removed part of the reflected signal, which directional filtering has preserved, according to an implementation. The spatially varying slope filter has a narrow reject band in the f-k domain that minimizes signal smearing and distortion.

Field Data Example

FIGS. 12A-C show example results of applying the example spatially-varying directional filter approach (for example, the spatially-varying directional filter of FIG. 6, 660) to example field data, according to an implementation. FIG. 12A shows the input data 1210; FIG. 12B shows the filtered data 1220 according to the example process 600; and FIG. 12C shows the residual 1230 (the difference or noise removed) between the input data 1210 and the filtered data 1220.

The field data used in FIGS. 12A-C has a receiver group interval of 6.25 m and a maximum offset of 1000 m. Each receiver station includes 12 bunched geophones arranged randomly within a 25 cm radius from the center of the trace. This “single-sensor” acquisition can provide true high-resolution sampling, both spatially (to avoid aliasing) and temporally (no array stacking applied), of both the signal and noise for optimal implementation of noise prediction and suppression.

The main ground roll mode (dominant direct and back-scattered ground roll) has a velocity of about 1050 m/s and frequencies up to 40 Hz. The minimum wavelength is about 25 m. This ground roll is back-scattered and is visible on the time-offset domain, mainly at negative wave numbers. The high-amplitude, slow-scattered noise is due to wave propagation through near-surface heterogeneities, as shown in FIGS. 12A-C. The example spatially varying directional filter is applied to the field data. First, the interference between signal and noise is reduced using an f-k filter, and then the steerable filters are computed for the initial noise component. Then the most probable orientations are estimated.

In this example, the seismic data is divided into four offset patches. FIG. 13A shows histograms 1310, 1320, 1330 and 1340 of local slopes calculated using steerable filters, corresponding to patches or sections 1315, 1325, 1335 and 1345 shown in FIG. 13B, respectively, according to an implementation. The most probable orientation within four offset patches are estimated. In this example, two orientations, measured counter-clockwise from the positive y-axis and perpendicular to the scattered noise, have been predicted by the histograms 1310-1340 in FIG. 13A. The dashed lines 1314 and 1324 show the dominant orientation of 67.5° for patch one and two, and the dashed lines 1332 and 1342 show the dominant orientation of 70° for patch three and four, respectively. The change in slope of the direct and scattered surface waves depends on the shear wave velocity of uppermost layers. This could be due to either an increase in the shear wave velocity of the rock at far offsets, or due to attenuation of high frequencies traveling within slow, shallow-layer velocities and dominance of the low-frequency component that travels with faster, deeper layer velocities, assuming velocity increases with depth.

Based on the predicted slopes, a non-linear median filter steered toward the noise directions is applied, and the filtered image is shown in FIG. 12B. The residual image FIG. 12C shows that only locally linear noise features have been removed by the filter.

In some implementations, Normal Move Out (NMO) and running average filter can be applied to enhance the reflected signal. FIGS. 14A-C show input data 1410, filtered data 1420, and the residual (the difference or noise removed) 1430 of the application of the steered-median filter approach to field data after NMO, running average filter, and inverse NMO to enhance the reflections. Reflected P waves modeled using ray-tracing are shown in dashed lines 1440 and 1450.

FIGS. 15A-C illustrate input field data 1510, filtered field data 1520, and a residual (the difference or noise removed) 1530, respectively, of the application of a f-k filter to field data after NMO, running average filter, and inverse NMO to enhance reflections as illustrated in FIGS. 15A-B, according to an implementation. Compared to the example steered median filter approach in FIGS. 14A-C, the f-k filter not only caused smoothing and edge effects, but also removed part of the signal due to leakage in the transform domain as shown in FIGS. 15A-C.

The wavelength of the scattered waves is directly proportional to the size of the scatterers, and, as demonstrated by the difference plot 1530, the wavelength of the back-scattered surface waves is much shorter than direct-surface waves. This suggests that most of the near-surface scattering is due to individual or sharp interface scatterers.

Effects of 3D Heterogeneities on 2D Data Acquisition and Processing

The example techniques can be extended to 3D data as well. In the 2D earth models, receivers are in-line with the scatterers. In reality, the Earth is a 3D structure, and scatterers can exist in a cross-line direction. The scattered noise due to heterogeneities in the cross-line direction appears as non-linear events in 2D sections; however, the scattered waves appear as cones in 3D data. Therefore, sampling densely in both the in-line and cross-line directions makes it possible to capture and filter out the noise using 3D spatially varying directional filters.

The effects of 3D heterogeneities on 2D acquisition and processing are shown by simulating seismic elastic waves in 3D. The waves in 2D are modeled as cylindrical waves while they are modeled as spherical waves in 3D. Similarly, a circular scatterer and a point source in 2D can be represented by a cylindrical heterogeneity and a line of sources in 3D, respectively. These differences between the 2D and 3D modes not only cause alterations to the phase, but can also alter the amplitude decay due to geometrical spreading (that is, attenuation with propagating distance). Body waves attenuate as r⁻¹ in energy or r^−1/2 in amplitude in 2D, whereas in 3D they attenuate as r⁻² in energy or r⁻¹ in amplitude. On the other hand, surface waves do not attenuate with propagation distance in 2D, whereas in 3D they attenuate as r⁻¹ in energy or r^−1/2 in amplitude.

FIGS. 16A-C are plots 1610, 1620 and 1630 showing 3D finite difference simulation results when the scatterers are in-line with the receivers, according to an implementation. The results are similar to the 2D cases shown in FIGS. 5A-5F. Specifically, FIG. 16A shows the reference wavefield simulated using an earth model without scatterers; FIG. 16B shows the total wavefield simulated using the earth model with scatterers; and FIG. 16C shows the difference between the total and reference wavefields.

FIGS. 17A-C are plots 1710, 1720 and 1730 showing 3D finite difference simulation results when the scatterers are in the cross-line direction, according to an implementation. Similarly, FIG. 17A shows the reference wavefield simulated using an earth model without scatterers; FIG. 17B shows the total wavefield simulated using the earth model with scatterers; and FIG. 17C shows the difference between the total and reference wavefields. As illustrated, the apex of the scattered surface waves appear as non-linear arrivals (for example, hyperbolic as shown in 1720 and 1730) as opposed to the conventional linear direct surface waves. These features may rise in case of surface and body-to-surface wave scattering from scatterers in the cross-line direction. These features cannot be easily attenuated in 2D data using conventional linear noise removal techniques such as f-k filter.

In the case of direct surface wave generated by sources in the cross-line direction, it may be possible to correct for the non-linear move-out given that the source location and surface wave velocity are known. However, this is not trivial for the case of scattered body-to-surface waves excited by scatterers in the cross-line direction, as the locations of the sources, scatterers in this case, are not known.

In addition, the effects of 2D and 3D wave scattering from an irregular (Gaussian) interface are studied. The Gaussian surface can be constructed by generating a random Gaussian surface that is correlated (smoothed) with a Gaussian operator as shown in FIGS. 19A-D.

FIGS. 18A-D are plots 1810, 1820, 1830 and 1840 showing aspects of a 3D irregular (Gaussian) bedrock interface model, according to an implementation. Specifically, FIG. 18B shows a Gaussian surface profile with 70 m standard deviation. FIG. 18B shows a Gaussian smoothing operator with 5 m correlation length. FIG. 18C shows a smoothed surface generated by convolving the Gaussian surface with the smoothing operator. FIG. 18D shows an earth model with near-surface irregular bedrock interface at 15 m depth below the free surface and deeper flat reflector at 200 m depth.

FIGS. 19A-C are plots 1910, 1920, and 1930 showing aspects of finite difference results for the irregular (Gaussian) bedrock interface model simulated using 2D finite difference techniques. FIGS. 19D-F are plots 1940, 1950, and 1960 showing aspects of finite difference results for the irregular (Gaussian) bedrock interface model simulated using 3D finite difference techniques. FIGS. 19A and 19D show the total wavefield simulated using the model with Gaussian shallow interface. FIGS. 19B and 19E show the incident wavefield simulated using the model with plane shallow interface. FIGS. 19C and 19F show scattered wavefields (that is, the difference between the total and incident wavefields).

The 3D simulations as shown in FIGS. 19D-F include scattering phases coming from the cross-line direction. The 2D velocity model was taken as a slice from the 3D model, at the same location of the receiver line. The results for the model with a planar shallow interface look similar for the 2D and 3D cases, except the weak amplitude due to spherical divergence in the 3D case. On the other hand, the simulated data for the model with Gaussian bedrock interface are different in 2D than in 3D. The 2D case includes scattering from only the in-line direction, and thus the scattered surface waves are locally linear. The 3D case, however, includes scattering from both the in-line and cross-line directions, and thus the recorded scattered surface waves are no longer linear.

FIG. 20 is a plot showing an example 3D earth model 2000 with multiple dipping layers 2010, 2020 and 2030 and near-surface scatterers 2040, according to an implementation. The multiple scatterers 2040 are randomly distributed in 3D.

FIGS. 21A and 21B are plots 2110 and 2120 showing aspects of modeling and filtering of scattered body-to-surface waves in 3D, according to an implementation. Specifically, FIG. 21A shows simulated 3D finite difference results of input data based on the 3D earth model 2000 in FIG. 20. As shown in FIG. 21A, the scattered waves exhibit mixed phases between linear arrivals due to the in-line scatterers and non-linear phases with different curvatures due to the randomly distributed scatterers in the cross-line direction. These noise features have frequency content comparable to the size of the scatterers. An amplitude comparable to the reflected signal can mask the entire data in the time-space domain. Therefore, filters based on frequency, amplitude contrast, or windowing are not feasible in this case. However, the scattered noise appears as cones in 3D, and therefore sampling densely in both the in-line and cross-line directions makes the problem similar to the 2D case, in which it is possible to capture and filter out the cone.

FIG. 21B shows estimated (filtered) signal after application of a 3D f-k filter to a spatially dense sampled dataset. The scattered waves have a cone shape in 3D; showing circular phases in the time slices (see FIG. 21A), while they appear as non-linear surface waves in the time-offset domain. Similar to the 2D case, the filtered signal is obtained according to the example process 600 except that, in 660, a directional filter (for example, f-k filter) is applied in 3D instead of a non-linear median filter as shown for the 2D case. The receivers are divided into patches in the spatial domain based on offset. The dominant local-slopes are predicted by computing the 3D steerable filters (see FIG. 22B) to guide a 3D f-k filter.

FIGS. 22A and 22B are plots 2210 and 2220 showing aspects of application of a 3D f-k filter to spatially dense sampled 3D simulated data, according to an implementation. Specifically, FIG. 22A shows the difference between the input data and filtered data (filtered noise) shown in FIGS. 21A and 21B, respectively. FIG. 22B shows a histogram 2230 of local-slopes calculated using 3D steerable filters showing the dominant slope 2240 for the receiver patch. The histogram distribution 2230 in FIG. 22B corresponds to the receiver patch highlighted in red in FIG. 22A.

Example techniques are described for estimating dominant local slopes in seismic data as a function of offset (or patch of offsets) using steerable filters and to separate the signal and noise components using a directional filter based on dip, slope, or velocity. The slope estimation step using steerable filters requires only a linear combination of a set of basis filters at fixed orientations to synthesize an image filtered at an arbitrary orientation. This makes the process very efficient. The example techniques can handle locally linear coherent noise that has small amplitude contrast, and varying slope with offset. The example techniques have been successfully implemented on synthetic and field data for the separation of scattered surface waves from reflected body waves. The results show that the example techniques are superior to conventional f-k techniques.

Although only surface waves excited by near surface body wave scattering are discussed, the example techniques can also handle surface-to-surface wave scattering. The method can also be applied to ocean-bottom-cable (OBC) marine data, in which the sea bottom acts in a similar way to the free surface; therefore, propagating waves can scatter to Scholte waves by irregularities or heterogeneities near the sea bottom.

FIG. 23 is a block diagram 2300 illustrating an example computer (for example, computer system 150 of FIG. 1A) used to provide computational functionalities associated with described algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure, according to an implementation. The example computer system 150 can be located at or near one or more well survey systems or at a remote location.

The illustrated computer 2302 is intended to encompass any computing device such as a server, desktop computer, laptop/notebook computer, wireless data port, smart phone, personal data assistant (PDA), tablet computing device, one or more processors within these devices, or any other suitable processing device, including both physical or virtual instances (or both) of the computing device. Additionally, the computer 2302 may comprise a computer that includes an input device, such as a keypad, keyboard, touch screen, or other device that can accept user information, and an output device that conveys information associated with the operation of the computer 2302, including digital data, visual, or audio information (or a combination of information), or a GUI.

The computer 2302 can serve in a role as a client, network component, a server, a database or other persistency, or any other component (or a combination of roles) of a computer system for performing the subject matter described in the instant disclosure. The illustrated computer 2302 is communicably coupled with a network 2330. In some implementations, one or more components of the computer 2302 may be configured to operate within environments, including cloud-computing-based, local, global, or other environment (or a combination of environments).

At a high level, the computer 2302 is an electronic computing device operable to receive, transmit, process, store, or manage data and information associated with the described subject matter. According to some implementations, the computer 2302 may also include or be communicably coupled with an application server, e-mail server, web server, caching server, streaming data server, business intelligence (BI) server, or other server (or a combination of servers).

The computer 2302 can receive requests over network 2330 from a client application (for example, executing on another computer 2302) and responding to the received requests by processing the said requests in an appropriate software application. In addition, requests may also be sent to the computer 2302 from internal users (for example, from a command console or by other appropriate access method), external or third-parties, other automated applications, as well as any other appropriate entities, individuals, systems, or computers.

Each of the components of the computer 2302 can communicate using a system bus 2303. In some implementations, any or all of the components of the computer 2302, both hardware or software (or a combination of hardware and software), may interface with each other or the interface 2304 (or a combination of both) over the system bus 2303 using an application programming interface (API) 2312 or a service layer 2313 (or a combination of the API 2312 and service layer 2313). The API 2312 may include specifications for routines, data structures, and object classes. The API 2312 may be either computer-language independent or dependent and refer to a complete interface, a single function, or even a set of APIs. The service layer 2313 provides software services to the computer 2302 or other components (whether or not illustrated) that are communicably coupled to the computer 2302. The functionality of the computer 2302 may be accessible for all service consumers using this service layer. Software services, such as those provided by the service layer 2313, provide reusable, defined business functionalities through a defined interface. For example, the interface may be software written in JAVA, C++, or other suitable language providing data in extensible markup language (XML) format or other suitable format. While illustrated as an integrated component of the computer 2302, alternative implementations may illustrate the API 2312 or the service layer 2313 as stand-alone components in relation to other components of the computer 2302 or other components (whether or not illustrated) that are communicably coupled to the computer 2302. Moreover, any or all parts of the API 2312 or the service layer 2313 may be implemented as child or sub-modules of another software module, enterprise application, or hardware module without departing from the scope of this disclosure.

The computer 2302 includes an interface 2304. Although illustrated as a single interface 2304 in FIG. 23, two or more interfaces 2304 may be used according to particular needs, desires, or particular implementations of the computer 2302. The interface 2304 is used by the computer 2302 for communicating with other systems in a distributed environment that are connected to the network 2330 (whether illustrated or not). Generally, the interface 2304 comprises logic encoded in software or hardware (or a combination of software and hardware) and operable to communicate with the network 2330. More specifically, the interface 2304 may comprise software supporting one or more communication protocols associated with communications such that the network 2330 or interface's hardware is operable to communicate physical signals within and outside of the illustrated computer 2302.

The computer 2302 includes a processor 2305. Although illustrated as a single processor 2305 in FIG. 23, two or more processors may be used according to particular needs, desires, or particular implementations of the computer 2302. Generally, the processor 2305 executes instructions and manipulates data to perform the operations of the computer 2302 and any algorithms, methods, functions, processes, flows, and procedures as described in the instant disclosure.

The computer 2302 also includes a memory 2306 that holds data for the computer 2302 or other components (or a combination of both) that can be connected to the network 2330 (whether illustrated or not). For example, memory 2306 can be a database storing data consistent with this disclosure. Although illustrated as a single memory 2306 in FIG. 23, two or more memories may be used according to particular needs, desires, or particular implementations of the computer 2302 and the described functionality. While memory 2306 is illustrated as an integral component of the computer 2302, in alternative implementations, memory 2306 can be external to the computer 2302.

The application 2307 is an algorithmic software engine providing functionality according to particular needs, desires, or particular implementations of the computer 2302, particularly with respect to functionality described in this disclosure. For example, application 2307 can serve as one or more components, modules, applications, etc. Further, although illustrated as a single application 2307, the application 2307 may be implemented as multiple applications 2307 on the computer 2302. In addition, although illustrated as integral to the computer 2302, in alternative implementations, the application 2307 can be external to the computer 2302.

There may be any number of computers 2302 associated with, or external to, a computer system containing computer 2302, each computer 2302 communicating over network 2330. Further, the term “client,” “user,” and other appropriate terminology may be used interchangeably as appropriate without departing from the scope of this disclosure. Moreover, this disclosure contemplates that many users may use one computer 2302, or that one user may use multiple computers 2302.

Described implementations of the subject matter can include one or more features, alone or in combination.

For example, in a first implementation, a computer-implemented method of processing seismic data for reservoir analysis, comprising: accessing, by operation of one or more computers, seismic data of a subterranean region; computing an image orientation of each data sample of the seismic data dividing, in a spatial domain, the seismic data into a plurality of sections based on spatial attributes of the seismic data, each section comprising respective seismic data samples having the spatial attributes within the respective section; for each of the plurality of sections: determining a dominant orientation of the section; and predicting a noise component of the seismic data based on the dominant orientation of the section.

The foregoing and other described implementations can each optionally include one or more of the following features:

A first feature, combinable with any of the following features, further comprising removing the noise component from the seismic data based on an image orientation of the predicted noise component.

A second feature, combinable with any of the previous or following features, further comprising applying dip decomposition to reduce interference prior to computing the image orientation of each data sample of the seismic data.

A third feature, combinable with any of the previous or following features, wherein computing an image orientation of each data sample of the seismic data comprises using a steerable filter, a slant-stacking filter, or a plane-wave destruction filter.

A fourth feature, combinable with any of the previous or following features, wherein the steerable filter comprises a linear combination of a set of basis filters at fixed orientations.

A fifth feature, combinable with any of the previous or following features, wherein computing an image orientation of each data sample comprises using a plurality of oriented filters with supplementary orientations to identify the image orientation of the data sample.

A sixth feature, combinable with any of the previous or following features, wherein determining a dominant orientation of the section comprises determining the dominant orientation of the section based on probability distributions of respective image orientations of data samples within the section.

A seventh feature, combinable with any of the previous or following features, wherein predicting a noise component comprises applying, to each section of the plurality of sections, a directional filter or a non-linear media filter aligned with the dominant orientation of the section.

An eighth feature, combinable with any of the previous or following features, wherein the directional filter comprises a local velocity, slope, or dip filter applicable in a time-space domain or in an f-k or tau-p domain.

A ninth feature, combinable with any of the previous or following features, wherein the spatial attributes of the seismic data comprises offset values of the seismic data, the offset values representing distance between a seismic source and a seismic receiver of the seismic data.

In a second implementation, a non-transitory, computer-readable medium storing one or more instructions executable by a computer system to perform operations for imaging and locating near-surface heterogeneities, comprising: accessing seismic data of a subterranean region; computing an image orientation of each data sample of the seismic data; dividing, in a spatial domain, the seismic data into a plurality of sections based on spatial attributes of the seismic data, each section comprising respective seismic data samples having the spatial attributes within the respective section; for each of the plurality of sections: determining a dominant orientation of the section; and predicting a noise component of the seismic data based on the dominant orientation of the section.

The foregoing and other described implementations can each optionally include one or more of the following features:

A first feature, combinable with any of the following features, further comprising one or more instructions to remove the noise component from the seismic data based on an image orientation of the predicted noise component.

A second feature, combinable with any of the previous or following features, further comprising one or more instructions to apply dip decomposition to reduce interference prior to computing the image orientation of each data sample of the seismic data.

A third feature, combinable with any of the previous or following features, wherein computing an image orientation of each data sample of the seismic data comprises using a steerable filter, a slant-stacking filter, or a plane-wave destruction filter.

A fourth feature, combinable with any of the previous or following features, wherein the steerable filter comprises a linear combination of a set of basis filters at fixed orientations.

A fifth feature, combinable with any of the previous or following features, wherein computing an image orientation of each data sample comprises using a plurality of oriented filters with supplementary orientations to identify the image orientation of the data sample.

A sixth feature, combinable with any of the previous or following features, wherein determining a dominant orientation of the section comprises determining the dominant orientation of the section based on probability distributions of respective image orientations of data samples within the section.

A seventh feature, combinable with any of the previous or following features, wherein predicting a noise component comprises applying, to each section of the plurality of sections, a directional filter or a non-linear media filter aligned with the dominant orientation of the section.

An eighth feature, combinable with any of the previous or following features, computer-readable medium of claim 18, wherein the directional filter comprises a local velocity, slope, or dip filter applicable in a time-space domain or in an f-k or tau-p domain.

A ninth feature, combinable with any of the previous or following features, wherein the spatial attributes of the seismic data comprises offset values of the seismic data, the offset values representing distance between a seismic source and a seismic receiver of the seismic data.

In a third implementation, a computer-implemented system to perform operations for imaging and locating near-surface heterogeneities, comprising: a computer memory operable to store recorded near-surface scattered wave data; and a data processing apparatus interoperably coupled with the computer memory and configured to: a computer memory operable to store seismic data of a subterranean region; and a data processing apparatus operable to: access seismic data of a subterranean region; compute an image orientation of each data sample of the seismic data; divide, in a spatial domain, the seismic data into a plurality of sections based on spatial attributes of the seismic data, each section comprising respective seismic data samples having the spatial attributes within the respective section; and for each of the plurality of sections: determine a dominant orientation of the section; and predict a noise component of the seismic data based on the dominant orientation of the section.

The foregoing and other described implementations can each optionally include one or more of the following features:

A first feature, combinable with any of the following features, further operable to remove the noise component from the seismic data based on an image orientation of the predicted noise component.

A second feature, combinable with any of the previous or following features, further operable to apply dip decomposition to reduce interference prior to computing the image orientation of each data sample of the seismic data.

A third feature, combinable with any of the previous or following features, wherein computing an image orientation of each data sample of the seismic data comprises using a steerable filter, a slant-stacking filter, or a plane-wave destruction filter.

A fourth feature, combinable with any of the previous or following features, wherein the steerable filter comprises a linear combination of a set of basis filters at fixed orientations.

A fifth feature, combinable with any of the previous or following features, wherein computing an image orientation of each data sample comprises using a plurality of oriented filters with supplementary orientations to identify the image orientation of the data sample.

A sixth feature, combinable with any of the previous or following features, wherein determining a dominant orientation of the section comprises determining the dominant orientation of the section based on probability distributions of respective image orientations of data samples within the section.

A seventh feature, combinable with any of the previous or following features, wherein predicting a noise component comprises applying, to each section of the plurality of sections, a directional filter or a non-linear media filter aligned with the dominant orientation of the section.

An eighth feature, combinable with any of the previous or following features, wherein the directional filter comprises a local velocity, slope, or dip filter applicable in a time-space domain or in an f-k or tau-p domain.

A ninth feature, combinable with any of the previous or following features, wherein the spatial attributes of the seismic data comprises offset values of the seismic data, the offset values representing distance between a seismic source and a seismic receiver of the seismic data.

Implementations of the subject matter and the functional operations described in this specification can be implemented in digital electronic circuitry, in tangibly embodied computer software or firmware, in computer hardware, including the structures disclosed in this specification and their structural equivalents, or in combinations of one or more of them. Implementations of the subject matter described in this specification can be implemented as one or more computer programs, that is, one or more modules of computer program instructions encoded on a tangible, non-transitory computer-storage medium for execution by, or to control the operation of, data processing apparatus. Alternatively or in addition, the program instructions can be encoded on an artificially generated propagated signal, for example, a machine-generated electrical, optical, or electromagnetic signal that is generated to encode information for transmission to suitable receiver apparatus for execution by a data processing apparatus. The computer-storage medium can be a machine-readable storage device, a machine-readable storage substrate, a random or serial access memory device, or a combination of computer-storage mediums.

The terms “data processing apparatus,” “computer,” or “electronic computer device” (or equivalent as understood by one of ordinary skill in the art) refer to data processing hardware and encompass all kinds of apparatus, devices, and machines for processing data, including by way of example, a programmable processor, a computer, or multiple processors or computers. The apparatus can also be or further include special purpose logic circuitry, for example, a central processing unit (CPU), an FPGA (field programmable gate array), or an ASIC (application-specific integrated circuit). In some implementations, the data processing apparatus or special purpose logic circuitry (or a combination of the data processing apparatus or special purpose logic circuitry) may be hardware- or software-based (or a combination of both hardware- and software-based). The apparatus can optionally include code that creates an execution environment for computer programs, for example, code that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of execution environments. The present disclosure contemplates the use of data processing apparatuses with or without conventional operating systems, for example LINUX, UNIX, WINDOWS, MAC OS, ANDROID, IOS or any other suitable conventional operating system.

A computer program, which may also be referred to or described as a program, software, a software application, a module, a software module, a script, or code, can be written in any form of programming language, including compiled or interpreted languages, or declarative or procedural languages, and it can be deployed in any form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program may, but need not, correspond to a file in a file system. A program can be stored in a portion of a file that holds other programs or data, for example, one or more scripts stored in a markup language document, in a single file dedicated to the program in question, or in multiple coordinated files, for example, files that store one or more modules, sub-programs, or portions of code. A computer program can be deployed to be executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a communication network. While portions of the programs illustrated in the various figures are shown as individual modules that implement the various features and functionality through various objects, methods, or other processes, the programs may instead include a number of sub-modules, third-party services, components, libraries, and such, as appropriate. Conversely, the features and functionality of various components can be combined into single components as appropriate.

The processes and logic flows described in this specification can be performed by one or more programmable computers executing one or more computer programs to perform functions by operating on input data and generating output. The processes and logic flows can also be performed by, and apparatus can also be implemented as, special purpose logic circuitry, for example, a CPU, an FPGA, or an ASIC.

Computers suitable for the execution of a computer program can be based on general or special purpose microprocessors, both, or any other kind of CPU. Generally, a CPU will receive instructions and data from a read-only memory (ROM) or a random access memory (RAM) or both. The essential elements of a computer are a CPU for performing or executing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to, receive data from or transfer data to, or both, one or more mass storage devices for storing data, for example, magnetic, magneto-optical disks, or optical disks. However, a computer need not have such devices. Moreover, a computer can be embedded in another device, for example, a mobile telephone, a personal digital assistant (PDA), a mobile audio or video player, a game console, a global positioning system (GPS) receiver, or a portable storage device, for example, a universal serial bus (USB) flash drive, to name just a few.

Computer-readable media (transitory or non-transitory, as appropriate) suitable for storing computer program instructions and data include all forms of non-volatile memory, media and memory devices, including by way of example semiconductor memory devices, for example, erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), and flash memory devices; magnetic disks, for example, internal hard disks or removable disks; magneto-optical disks; and CD-ROM, DVD+/−R, DVD-RAM, and DVD-ROM disks. The memory may store various objects or data, including caches, classes, frameworks, applications, backup data, jobs, web pages, web page templates, database tables, repositories storing dynamic information, and any other appropriate information including any parameters, variables, algorithms, instructions, rules, constraints, or references thereto. Additionally, the memory may include any other appropriate data, such as logs, policies, security or access data, reporting files, as well as others. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry.

To provide for interaction with a user, implementations of the subject matter described in this specification can be implemented on a computer having a display device, for example, a CRT (cathode ray tube), LCD (liquid crystal display), LED (Light Emitting Diode), or plasma monitor, for displaying information to the user and a keyboard and a pointing device, for example, a mouse, trackball, or trackpad by which the user can provide input to the computer. Input may also be provided to the computer using a touchscreen, such as a tablet computer surface with pressure sensitivity, a multi-touch screen using capacitive or electric sensing, or other type of touchscreen. Other kinds of devices can be used to provide for interaction with a user as well; for example, feedback provided to the user can be any form of sensory feedback, for example, visual feedback, auditory feedback, or tactile feedback; and input from the user can be received in any form, including acoustic, speech, or tactile input. In addition, a computer can interact with a user by sending documents to and receiving documents from a device that is used by the user; for example, by sending web pages to a web browser on a user's client device in response to requests received from the web browser.

The term “graphical user interface,” or “GUI,” may be used in the singular or the plural to describe one or more graphical user interfaces and each of the displays of a particular graphical user interface. Therefore, a GUI may represent any graphical user interface, including but not limited to, a web browser, a touch screen, or a command line interface (CLI) that processes information and efficiently presents the information results to the user. In general, a GUI may include a plurality of user interface (UI) elements, some or all associated with a web browser, such as interactive fields, pull-down lists, and buttons operable by the business suite user. These and other UI elements may be related to or represent the functions of the web browser.

Implementations of the subject matter described in this specification can be implemented in a computing system that includes a back-end component, for example, as a data server, or that includes a middleware component, for example, an application server, or that includes a front-end component, for example, a client computer having a graphical user interface or a Web browser through which a user can interact with an implementation of the subject matter described in this specification, or any combination of one or more such back-end, middleware, or front-end components. The components of the system can be interconnected by any form or medium of wireline or wireless digital data communication (or a combination of data communication), for example, a communication network. Examples of communication networks include a local area network (LAN), a radio access network (RAN), a metropolitan area network (MAN), a wide area network (WAN), Worldwide Interoperability for Microwave Access (WIMAX), a wireless local area network (WLAN) using, for example, 802.11 a/b/g/n or 802.20 (or a combination of 802.11x and 802.20 or other protocols consistent with this disclosure), all or a portion of the Internet, or any other communication system or systems at one or more locations (or a combination of communication networks). The network may communicate with, for example, Internet Protocol (IP) packets, Frame Relay frames, Asynchronous Transfer Mode (ATM) cells, voice, video, data, or other suitable information (or a combination of communication types) between network addresses.

The computing system can include clients and servers. A client and server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other.

In some implementations, any or all of the components of the computing system, both hardware or software (or a combination of hardware and software), may interface with each other or the interface using an application programming interface (API) or a service layer (or a combination of API and service layer). The API may include specifications for routines, data structures, and object classes. The API may be either computer language independent or dependent and refer to a complete interface, a single function, or even a set of APIs. The service layer provides software services to the computing system. The functionality of the various components of the computing system may be accessible for all service consumers using this service layer. Software services provide reusable, defined business functionalities through a defined interface. For example, the interface may be software written in JAVA, C++, or other suitable language providing data in extensible markup language (XML) format or other suitable format. The API or service layer (or a combination of the API and the service layer) may be an integral or a stand-alone component in relation to other components of the computing system. Moreover, any or all parts of the service layer may be implemented as child or sub-modules of another software module, enterprise application, or hardware module without departing from the scope of this disclosure.

While this specification contains many specific implementation details, these should not be construed as limitations on the scope of any invention or on the scope of what may be claimed, but rather as descriptions of features that may be specific to particular implementations of particular inventions. Certain features that are described in this specification in the context of separate implementations can also be implemented in combination in a single implementation. Conversely, various features that are described in the context of a single implementation can also be implemented in multiple implementations separately or in any suitable sub-combination. Moreover, although features may be described above as acting in certain combinations and even initially claimed as such, one or more features from a claimed combination can in some cases be excised from the combination, and the claimed combination may be directed to a sub-combination or variation of a sub-combination.

Particular implementations of the subject matter have been described. Other implementations, alterations, and permutations of the described implementations are within the scope of the following claims as will be apparent to those skilled in the art. While operations are depicted in the drawings or claims in a particular order, this should not be understood as requiring that such operations be performed in the particular order shown or in sequential order, or that all illustrated operations be performed (some operations may be considered optional), to achieve desirable results. In certain circumstances, multitasking or parallel processing (or a combination of multitasking and parallel processing) may be advantageous and performed as deemed appropriate.

Moreover, the separation or integration of various system modules and components in the implementations described above should not be understood as requiring such separation or integration in all implementations, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products.

Accordingly, the above description of example implementations does not define or constrain this disclosure. Other changes, substitutions, and alterations are also possible without departing from the spirit and scope of this disclosure.

Furthermore, any claimed implementation below is considered to be applicable to at least a computer-implemented method; a non-transitory, computer-readable medium storing computer-readable instructions to perform the computer-implemented method; and a computer system comprising a computer memory interoperably coupled with a hardware processor configured to perform the computer-implemented method or the instructions stored on the non-transitory, computer-readable medium.

Claims

1. A computer-implemented method of processing seismic data for reservoir analysis, the method comprising:

accessing, by operation of one or more computers, seismic data of a subterranean region;

computing an image orientation of each data sample of the seismic data;

dividing, in a spatial domain, the seismic data into a plurality of sections based on spatial attributes of the seismic data, each section comprising respective seismic data samples having the spatial attributes within the respective section;

for each of the plurality of sections: determining a dominant orientation of the section; and predicting a noise component of the seismic data based on the dominant orientation of the section.

2. The computer-implemented method of claim 1, further comprising removing the noise component from the seismic data based on an image orientation of the predicted noise component.

3. The computer-implemented method of claim 1, further comprising applying dip decomposition to reduce interference prior to computing the image orientation of each data sample of the seismic data.

4. The computer-implemented method of claim 1, wherein computing an image orientation of each data sample of the seismic data comprises using a steerable filter, a slant-stacking filter, or a plane-wave destruction filter.

5. The computer-implemented method of claim 4, wherein the steerable filter comprises a linear combination of a set of basis filters at fixed orientations.

6. The computer-implemented method of claim 4, wherein computing an image orientation of each data sample comprises using a plurality of oriented filters with supplementary orientations to identify the image orientation of the data sample.

7. The computer-implemented method of claim 1, wherein determining a dominant orientation of the section comprises determining the dominant orientation of the section based on probability distributions of respective image orientations of data samples within the section.

8. The computer-implemented method of claim 1, wherein predicting a noise component comprises applying, to each section of the plurality of sections, a directional filter or a non-linear media filter aligned with the dominant orientation of the section.

9. The computer-implemented method of claim 8, wherein the directional filter comprises a local velocity, slope, or dip filter applicable in a time-space domain or in an f-k or tau-p domain.

10. The computer-implemented method of claim 1, wherein the spatial attributes of the seismic data comprises offset values of the seismic data, the offset values representing distance between a seismic source and a seismic receiver of the seismic data.

11. A non-transitory, computer-readable medium storing one or more instructions executable by a computer system to perform operations comprising:

accessing seismic data of a subterranean region;

computing an image orientation of each data sample of the seismic data;

dividing, in a spatial domain, the seismic data into a plurality of sections based on spatial attributes of the seismic data, each section comprising respective seismic data samples having the spatial attributes within the respective section;

for each of the plurality of sections: determining a dominant orientation of the section; and predicting a noise component of the seismic data based on the dominant orientation of the section.

12. The non-transitory, computer-readable medium of claim 11, further comprising one or more instructions to remove the noise component from the seismic data based on an image orientation of the predicted noise component.

13. The non-transitory, computer-readable medium of claim 11, further comprising one or more instructions to apply dip decomposition to reduce interference prior to computing the image orientation of each data sample of the seismic data.

14. The non-transitory, computer-readable medium of claim 11, wherein computing an image orientation of each data sample of the seismic data comprises using a steerable filter, a slant-stacking filter, or a plane-wave destruction filter.

15. The non-transitory, computer-readable medium of claim 14, wherein the steerable filter comprises a linear combination of a set of basis filters at fixed orientations.

16. The non-transitory, computer-readable medium of claim 14, wherein computing an image orientation of each data sample comprises using a plurality of oriented filters with supplementary orientations to identify the image orientation of the data sample.

17. The non-transitory, computer-readable medium of claim 11, wherein determining a dominant orientation of the section comprises determining the dominant orientation of the section based on probability distributions of respective image orientations of data samples within the section.

18. The non-transitory, computer-readable medium of claim 11, wherein predicting a noise component comprises applying, to each section of the plurality of sections, a directional filter or a non-linear media filter aligned with the dominant orientation of the section.

19. The non-transitory, computer-readable medium of claim 18, wherein the directional filter comprises a local velocity, slope, or dip filter applicable in a time-space domain or in an f-k or tau-p domain.

20. The computer-implemented method of claim 11, wherein the spatial attributes of the seismic data comprises offset values of the seismic data, the offset values representing distance between a seismic source and a seismic receiver of the seismic data.

21. A computer-implemented system comprising one or more computers that include:

a computer memory operable to store seismic data of a subterranean region; and

a data processing apparatus operable to: access seismic data of a subterranean region; compute an image orientation of each data sample of the seismic data; divide, in a spatial domain, the seismic data into a plurality of sections based on spatial attributes of the seismic data, each section comprising respective seismic data samples having the spatial attributes within the respective section; and for each of the plurality of sections: determine a dominant orientation of the section; and predict a noise component of the seismic data based on the dominant orientation of the section.

22. The computer-implemented system of claim 21, further operable to remove the noise component from the seismic data based on an image orientation of the predicted noise component.

23. The computer-implemented system of claim 21, further operable to apply dip decomposition to reduce interference prior to computing the image orientation of each data sample of the seismic data.

24. The computer-implemented system of claim 21, wherein computing an image orientation of each data sample of the seismic data comprises using a steerable filter, a slant-stacking filter, or a plane-wave destruction filter.

25. The computer-implemented system of claim 24, wherein the steerable filter comprises a linear combination of a set of basis filters at fixed orientations.

26. The computer-implemented system of claim 24, wherein computing an image orientation of each data sample comprises using a plurality of oriented filters with supplementary orientations to identify the image orientation of the data sample.

27. The computer-implemented system of claim 21, wherein determining a dominant orientation of the section comprises determining the dominant orientation of the section based on probability distributions of respective image orientations of data samples within the section.

28. The computer-implemented system of claim 21, wherein predicting a noise component comprises applying, to each section of the plurality of sections, a directional filter or a non-linear media filter aligned with the dominant orientation of the section.

29. The computer-implemented system of claim 28, wherein the directional filter comprises a local velocity, slope, or dip filter applicable in a time-space domain or in an f-k or tau-p domain.

30. The computer-implemented system of claim 21, wherein the spatial attributes of the seismic data comprises offset values of the seismic data, the offset values representing distance between a seismic source and a seismic receiver of the seismic data.