Patent Description:
Embodiments of the subject matter disclosed herein generally relate to methods and systems and, more particularly, to mechanisms and techniques for processing seismic data sets that are recorded with different spatial sampling and different temporal bandwidths.

Reflection seismology is a method of geophysical exploration to determine the properties of a portion of a subsurface layer in the earth, information that is especially helpful in the oil and gas industry. Marine reflection seismology is based on the use of a controlled source that sends energy waves into the earth. By measuring the time it takes for the reflections to come back to plural receivers, it is possible to estimate the depth and/or composition of the features causing such reflections. These features may be associated with subterranean hydrocarbon deposits.

For marine applications, a seismic survey system <NUM>, as illustrated in <FIG>, includes a vessel <NUM> that tows plural streamers <NUM> (only one is visible in the figure) and a seismic source array <NUM>. Streamer <NUM> is attached through a lead-in cable (or other cables) <NUM> to vessel <NUM>, while source array <NUM> is attached through an umbilical <NUM> to the vessel. A head float <NUM>, which floats at the water surface <NUM>, is connected through a cable <NUM> to the head end 110A of streamer <NUM>, while a tail buoy <NUM> is connected, through a similar cable <NUM>, to the tail end 110B of streamer <NUM>. Head float <NUM> and tail buoy <NUM> are used, among other things, to maintain the streamer's depth. Seismic sensors <NUM> are distributed along the streamer and configured to record seismic data. Seismic sensors <NUM> may include a hydrophone, geophone, accelerometer or a combination thereof. Positioning devices <NUM> are attached along the streamer and controlled by a controller <NUM> for adjusting a position of the streamer according to a survey plan.

Source array <NUM> has plural source elements <NUM>, which are typically air guns. The source elements are attached to a float <NUM> to travel at desired depths below the water surface <NUM>. During operation, vessel <NUM> follows a predetermined path T while source elements (usually air guns) <NUM> emit seismic waves <NUM>. These waves bounce off the ocean bottom <NUM> and other layer interfaces below the ocean bottom <NUM> and propagate as reflected/refracted waves <NUM>, which are recorded by sensors <NUM>. The positions of both source elements <NUM> and recording sensors <NUM> may be estimated based on GPS systems <NUM> and recorded together with the seismic data in a storage device <NUM> onboard the vessel. Controller <NUM> has access to the seismic data and may be used to achieve quality control or even fully process the data. Controller <NUM> may also be connected to the vessel's navigation system and other elements of the seismic survey system, e.g., positioning devices <NUM>.

Alternatively, the seismic data may be collected in a land environment as illustrated in <FIG>. According to this scenario, a seismic acquisition system <NUM> includes plural receivers <NUM> positioned over an area <NUM> of a subsurface to be explored and in contact with, or below the surface <NUM> of, the ground. A number of dedicated seismic sources <NUM> are also placed on the surface <NUM> in an area <NUM>, in a vicinity of the area <NUM> of the receivers <NUM>. Note that a dedicated seismic source is defined as a device built by man with the main purpose of generating seismic waves to be used for a seismic survey. Alternatively, dedicated seismic sources <NUM> may be buried under surface <NUM>. A central recording device <NUM> is connected to the plurality of receivers <NUM> and placed, for example, in a station/truck <NUM>. Each dedicated seismic source <NUM> can be composed of a variable number of vibrators, typically between one and five, and can include a local controller <NUM>. A central controller <NUM> can be provided to coordinate the shooting times of sources <NUM>. A global positioning system (GPS) <NUM> can be used to time-correlate shooting of the dedicated seismic sources <NUM> and the recordings of the receivers <NUM>. With this configuration, dedicated seismic sources <NUM> are controlled to intentionally generate seismic waves, and the plurality of receivers <NUM> records waves reflected by oil and/or gas reservoirs and other structures. Patent application <CIT> discloses a method and system of seismic data processing including acquiring first and second prestack seismic datasets wherein the seismic datasets have a time-lapse relationship. Patent application <CIT> discloses source acquisition with multiple frequency components.

Patent application <CIT> discloses a method of conducting a marine seismic survey and of analysing produced data, comprising activating low-frequency sources and high frequency airguns, wherein low-frequency data are used to improve a model derived from the high-frequency data. Patent <CIT> discloses a method and a system for acquisition of seismic data, wherein high-frequency data are acquired at a dense spatial sampling and low-frequency data are acquired at a coarse spatial sampling.

Both in the land and marine seismic acquisition data, there are cases when two different sets of seismic data are acquired over the same desired subsurface. These two sets of seismic data are combined when generating an image of the subsurface. Traditionally, the data combination is possible because the two sets of seismic data are acquired to have the same spatial sampling, i.e., the data is recorded at the same locations, as schematically illustrated in <FIG>, where the first set of seismic data <NUM> is represented by circles and the second set of seismic data <NUM> is represented by crosses. Note that parameters <NUM> and <NUM> on the X and Y axes in <FIG> may be any of the x, y and z coordinates. However, it is possible that two sets of seismic data <NUM> and <NUM> have a different spatial sampling, as illustrated in <FIG>. If this is the case, existing seismic processing methods and algorithms are not capable of combining the two sets of seismic data. If the two seismic data sets further have different temporal bandwidths, in addition to different spatial sampling, then most traditional algorithms are not capable of processing this data.

Thus, there is a need to develop new algorithms and methods for processing seismic data sets having different spatial sampling and different temporal bandwidths.

According to one embodiment, there is a method for processing seismic data associated with a surveyed subsurface according to claim <NUM>. The method includes a step of receiving at least first and second seismic data sets that were recorded with different spatial sampling and different temporal bandwidths; a step of using the first seismic data to guide a processing of the second seismic data; and a step of generating an image of the subsurface based on processed second seismic data.

According to another embodiment, there is a computing seismic acquisition system configured to execute a seismic survey of a subsurface according to claim <NUM>, and including a device for processing seismic data associated with a surveyed subsurface. The computing device includes an interface configured to receive at least first and second seismic data sets that were recorded with different spatial sampling and different temporal bandwidths and a processor connected to the interface. The processor is configured to use the first seismic data to guide a processing of the second seismic data, and generate an image of the subsurface based on processed second seismic data.

According to an example, not covered by the claims, there is a method for processing seismic data associated with a surveyed subsurface. The method includes a step of receiving at least first and second seismic data sets that were recorded with different spatial sampling and different temporal bandwidths; interpolating the first seismic data to obtain interpolated first seismic data; combining the interpolated first seismic data with the second seismic data; and generating an image of the subsurface based on a combination of the interpolated first seismic data with the second seismic data.

The following description of the exemplary embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Instead, the scope of the invention is defined by the appended claims. The following embodiments are covered by the subject-matter of the claims with regard to two seismic data sets, one having a high-frequency content and the second one having a low-frequency content. However, the embodiments to be discussed next are not limited to two seismic data sets or to high- and low-frequency seismic data sets; the methods discussed herein can be applied to any number of seismic data sets and/or to other types of data sets as long as they have different spatial sampling and different temporal bandwidths, where by different temporal bandwidth it is understood that there can be some overlap in temporal bandwidth as long at the energy from the sources is separated, for example, if different pseudo-random signals are used for each source or if pseudo-orthogonal sequences are used.

Reference throughout the specification to "one embodiment" or "an embodiment" means that a particular feature, structure or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of the phrases "in one embodiment" or "in an embodiment" in various places throughout the specification is not necessarily referring to the same embodiment. Further, the particular features, structures or characteristics may be combined in any suitable manner in one or more embodiments.

According to an embodiment, a method for processing seismic data associated with a surveyed subsurface includes a step of receiving at least first and second seismic data sets that were recorded with different spatial sampling and different temporal bandwidths, a step of processing the second seismic data based on the first seismic data, and a step of generating an image of the subsurface based on the first and second seismic data.

More specifically, according to an embodiment, it is desired to process at least two seismic data sets <NUM> and <NUM> with different bandwidths (e.g., low-frequency LF, and high-frequency HF) which also have different spatial sampling. Note that data bandwidth may also be a function of the source and/or receiver group array and environment near the source/receiver, in particular, the free surface ghost in marine acquisition. In the marine context, the free surface ghost amplifies energy at some frequencies and attenuates energy at other frequencies. The difference in temporal bandwidth may be a function of frequencies being emitted but also the depth of the source. The term "different spatial sampling" means that the two seismic data sets are not co-located in at least one of the x-coordinate, y-coordinate or z-coordinate for either the source and/or the receiver. Thus, in the following embodiments, these two seismic data sets are considered for exemplifying the various methods and processes. As noted above, the methods and processes to be discussed next can be applied to more than two seismic data sets.

Examples of physical scenarios where such seismic data sets may be acquired are:.

The datasets may include a combination of sources emitting different bandwidths and receivers designed to have sensitivity to different bandwidths.

These and other examples may include land or marine vibrator acquisition with low- and/or high-frequency sources. The acquisition may use single vibrators or vibrator arrays. The sources may be stationary or moving while emitting energy. The sources may emit energy with linear sweeps, step sweeps, bandlimited pseudo-random sweeps, pseudo-orthogonal sequences or any other type of signal. In the case of simultaneous shooting, the energy from the sources may have been separated by deblending.

The following discussion may relate to land, marine (e.g., ocean bottom survey, towed streamer), or transition zone data. The receivers may be hydrophones, accelerometers, differential pressure sensors, geophones, or other sensors configured to detect seismic energy. Seismic sources may be air guns, dynamite, land vibroseis, marine vibrators, weight drop, pingers, boomers, or other seismic sources. The source may include an array of elements. Source array elements may emit energy in a synchronized way or desynchronized (e.g., "popcorn source"). A combination of source types may be used. The sources may be excited independently (without overlap of recorded energy), or simultaneously (blended acquisition). In the case of blended acquisition, the energy may have been deblended. For vibrator acquisition (marine or land), multiple sources using different pseudo-random or pseudo-orthogonal signals may be used to help source separation. The data may be before or after source designature.

<FIG> illustrates a processing method for determining an image of the surveyed subsurface based on the two seismic data sets <NUM> and <NUM>. Note that <FIG> and associated description assumes that the two data sets <NUM> and <NUM> have the same spatial sampling. The method starts with step <NUM>, in which the two seismic data sets (or more) are input and merged (optionally) with navigation information. In step <NUM>, the seismic data sets may be deblended if they were recorded with a simultaneous shooting method. In step <NUM>, a designature method may be applied, i.e., converting a known source signature used in acquisition to a desired signature. The designature method may be applied as per conventional data, either based on vertical take-off or using directionality. It may be applied in combination with deblending in the case of simultaneous shooting.

In step <NUM>, one or more noise attenuation techniques may be applied to remove incoherent (random) or coherent noise. Random noise attenuation may include swell noise attenuation. One technique may use frequency-space (FX) deconvolution. This technique works on frequency slices, and it can be applied to LF and HF data independently or on blended data. Another technique is FX projection filtering. It works on frequency slices, and it can also be applied to LF and HF data independently or on blended data. Still another technique is time-domain prediction filtering. It can be applied to LF and HF data independently or on blended data. Filters derived on the LF data may be used to constrain filters to be applied to HF data. Another technique is common reflection surface denoise (CRS). The kinematics for this technique may be derived on LF and/or HF data and used to denoise LF and/or HF data. Another technique is the impulsive denoise. Often this step consists of defining noise to be attenuated, and then reconstructing the energy from surrounding unaffected data. Normally the reconstruction works on frequency slices and may be applied to LF and HF independently. Rank reduction filtering is another technique for random noise attenuation. It normally works on frequency slices, and it may be applied to LF and HF data independently. Another random noise attenuation technique is dip filtering. A dip estimate may be derived on the LF data and used to denoise the HF data. Note that all the techniques and/or methods that are known in the art by those skilled are not explained or detailed herein.

Regarding coherent noise attenuation techniques that may be applied to seismic data processing, they include ground/mud roll attenuation. This technique may use information from non-aliased low-frequency data to help attenuate high-frequency aliased noise. Another technique is demultiple, i.e., the attenuation or removal of those waves that suffer multiple reflections before being recorded. The demultiple technique includes a class of various methods, with some of them briefly discussed now. Gapped deconvolution is one member of the class and relies on time-domain operators that are derived and applied to remove multiples in the time or tau-p domain. Gapped deconvolution may be applied on LF and HF data separately or on the combination of LF and HF data. Radon demultiple is another member of the demultiple class, and it may include high resolution Radon. Least squares Radon is often applied to frequency slices independently and may be applied to LF and HF data separately. High-resolution Radon normally calculates sparseness weights on low-frequency non-aliased data to prevent aliasing in the Radon domain at high temporal frequencies. If two seismic data sets are used, it is possible to derive sparseness weights from LF data and use the weights to prevent aliasing on the HF data. Another member of the class is surface-related multiple attenuation (SRME). This is a convolution and summation process which may be applied to LF and HF data independently. Other members of the demultiple class include wave propagation demultiple, internal multiple attenuation (IME), etc..

Another technique for coherent noise attenuation is interference noise attenuation. Amongst other strategies, this technique may consist of impulsive denoise in the time or tau-p domain. As described before, this will be effective with application on the LF and HF data independently. Still another technique for coherent noise attenuation is shear wave attenuation. Those skilled in the art would appreciate that other noise attenuation methods and techniques may exist, and those discussed above were exemplary and not intended to be a complete list.

The method then advances to step <NUM>, which deghosts the seismic data sets, if they are marine data sets. In some cases, the processes of deghosting and data regularization (to be discussed next) may be combined. Data regularization is performed in step <NUM>, and it implies regularizing the seismic data beyond aliasing. To achieve this goal, model representations of the seismic data using information from LF data propagated to HF data are commonly used. Step <NUM> is discussed later in more detail.

In step <NUM>, the velocity model is built. Velocity model building is applied after data regularization, and this step may be applied as per conventional data. For full waveform inversion (FWI), a cost function may be defined to jointly match the LF and HF data (this is a similar discussion to least squares RTM given later in the application).

In step <NUM>, the processed seismic data is used to generate an image of the surveyed subsurface. Imaging may take various shapes, e.g., Kirchhoff, beam migration, wave equation, e.g., RTM, least squares migration, etc. This step is discussed later in more detail.

In step <NUM>, post-migration processing may be applied. This processing may include interpolation (e.g., FX interpolation). FX interpolation uses operators calculated at low frequencies to interpolate high-frequency data. While after migration the LF and HF data will be co-located, this process may be applied before data regularization by deriving prediction filters on the LF data and using them for interpolation of the HF data. Another post-migration process is random noise attenuation, which is often applied on frequency slices. This process may be applied to HF and LF data separately. Still another process is Q-compensation, which is often applied on frequency slices. This process may be applied to HF and LF data separately. Another process is frequency filtering, which is applied on frequency slices. This process also may be applied to HF and LF data separately.

As discussed above, some of the steps (e.g., frequency filtering) may be applied to the low- and high-frequency data independently, but other processes use energy at one temporal frequency to process data at another temporal frequency (e.g., Radon demultiple). Many processing steps use low temporal frequency information to guide results on high frequencies. Often, these steps relate to cases where the input data is aliased. Some examples of these steps known in the literature are as follows:
Interpolation/regularization.

Traditional wisdom when dealing with any of the above steps, for two or more seismic data sets recorded with the same spatial sampling, is to (<NUM>) calculate a model domain representation of the seismic data for low frequencies and (<NUM>) to use the low-frequency model to constrain the model representation for high frequencies (e.g., for dealiasing). An alternative approach is to (<NUM>) calculate a filter at low frequencies and (<NUM>) use the filter to process the high frequencies (e.g., interpolation or denoise techniques).

However, the above-noted procedures require that the low-frequency and high-frequency data are sampled consistently in space. Thus, according to an embodiment, the high- and low-frequency data are considered to not be consistently sampled in space. If this is the case, new methods are discussed for processing seismic data sets having different spatial sampling to obtain an image of the surveyed subsurface.

In this respect, the more accurate the image, the better the chances of drilling at the correct location and of finding oil and gas reservoirs with a reduced extraction price. Thus, the embodiments to be discussed next improve the technological processes of finding and extracting oil and gas by combining two seismic data sets recorded with different spatial sampling and different temporal bandwidths.

One such embodiment is now discussed with regard to Radon demultiple, which is a component of step <NUM> of noise attenuation. As described in the art, parabolic Radon demultiple may be problematic due to spatial aliasing of data at high temporal frequencies. For this reason, it is common to use a low-frequency parabolic Radon model domain to guide the model derivation for high frequencies. In the event there is a difference in spatial positioning between low- and high-frequency data, the traditional scheme fails to produce accurate results.

To rectify this problem, consider low- and high-frequency seismic data sets relating to two seismic sources d1 and d2, respectively, as illustrated in <FIG> shows a seismic acquisition system <NUM> that includes at least a vessel <NUM> that tows a streamer spread <NUM>. Streamer spread <NUM> includes at least one streamer <NUM> having one or more types of seismic sensors <NUM> (only four are shown for simplicity), e.g., hydrophones, geophones, accelerometers, etc. Vessel <NUM> may tow the two sources d1 <NUM> and d2 <NUM>. Each source may be a source array that includes one or more source elements. As discussed above, the source element may be vibratory, impulsive, etc. Seismic data (low- and high-frequency) generated by the sources and recorded by the seismic sensors may be recorded during the same seismic survey. In one application, the two data sets may be recorded during two different seismic surveys, as is the case for a 4D survey.

Also consider that sources d1 and d2 are in relatively close proximity to each other, e.g., having a <NUM> or less cross-line separation. Source d1 may be configured to generate high-frequency signals (e.g., <NUM> to <NUM>) while source d2 may be configured to generate low-frequency signals (e.g., <NUM> to <NUM>).

It is common in data processing to sort the recorded data to 2D common midpoint gathers (CMPs), where traces are grouped based on a position located along a line between the source and receiver. <FIG> illustrates three consecutive shots <NUM>', <NUM>" and <NUM>‴ for high-frequency source <NUM>, three consecutive shots <NUM>', <NUM>" and <NUM>‴ for low-frequency source <NUM>, and three seismic receivers 507A-C distributed along streamer <NUM>. Note that the three high-frequency shots recorded by receivers 507A-C, although recorded at different instants, correspond to the same CMP gather, i.e., CMPHF. The same is true for low-frequency CMP gather CMPLF.

<FIG> shows that lateral displacement d between CMPLF and CMPHF is half the source separation D. Because low-frequency energy generally varies spatially more slowly than high-frequency energy (e.g., Fresnel zone), and assuming that sources d1 and d2 are sufficiently close together, the inventor of this application has discovered that it is possible to use the energy from low-frequency source d2 <NUM> to guide the Radon transform for high-frequency source d1 <NUM>, although the two types of energies were recorded with different spatial sampling.

To confirm this novel feature, CMP gathers are collected consistently from two sources in a configuration as shown in <FIG>. In this case, each source <NUM> and <NUM> is a source array composed of plural air guns, and each source emits broadband data, i.e., both low and high frequencies. This choice offers the possibility to compare Radon demultiple results obtained by using (i) low-frequency energy from source <NUM> to guide the high-frequency energy from source <NUM>, and (ii) low-frequency energy from source <NUM> to guide the high-frequency energy from source <NUM>, thus testing the above novel assertion.

The proposed test follows the flow illustrated in <FIG> shows that in step <NUM>, CMP data for source <NUM> is received. This data is illustrated in <FIG>. In step <NUM>, a high-resolution Radon demultiple method is applied to source <NUM>'s CMP using low-frequency data to guide the transform derivation of high-frequency data. All data is used from source <NUM>. The result of the high-resolution Radon demultiple method is illustrated in <FIG>.

In step <NUM>, the CMP data from source <NUM> is received. In step <NUM>, the CMP data from source <NUM> is low-cut at, for example, <NUM>, to generate the high-frequency data (see <FIG>), and in step <NUM>, the CMP data from source <NUM> is high-cut at, for example, <NUM>, to generate the low-frequency data (see <FIG>). A high-resolution Radon demultiple method is applied in step <NUM> to high-frequency data shown in <FIG>, where the low-frequency data shown in <FIG> is used to guide the high-frequency data of <FIG>. The results of the demultiple are shown in <FIG>. To compare the same results for the case illustrated in <FIG>, i.e., both the low- and high-frequency data received from source <NUM>, the demultiples illustrated in <FIG> are low-cut at <NUM> to generate the high-frequency result, which is illustrated in <FIG>. A difference calculated between the data illustrated in <FIG> (i.e., high-frequency data from source <NUM> guided with low-frequency data from source <NUM>, which represents the different spatial sampling case) and the data illustrated in <FIG> (i.e., high-frequency data from source <NUM> guided with low-frequency data from source <NUM>, which represents the same spatial sampling case) is calculated in step <NUM> and it is shown in <FIG> illustrates the minimal difference between the results of <FIG>, which validates the novel feature that seismic data sets of different spatial sampling may be used together for processing.

Thus, the results in <FIG> illustrate a similar effectiveness of high-frequency demultiple, whether the low frequencies for one source (<FIG>) or the other source (<FIG>) are used to guide the transform at high frequencies. The low amplitude difference illustrated in <FIG> highlights the strength of this approach. While the above method is an example of a high-resolution Radon demultiple, the use of low-frequency data to guide the processing of high-frequency data is relevant for other applications using a model space, as discussed later in this application.

Another embodiment is now discussed with regard to tau-p deconvolution, which is an implementation of gapped deconvolution. Gapped deconvolution has been discussed above with regard to <FIG>, step <NUM>, which is linked to noise attenuation. Gapped deconvolution in the tau-p domain is a robust and effective method to attenuate multiples in shallow water or land environments. The method includes the following steps: (<NUM>) forward tau-p transform the recorded seismic data, (<NUM>) apply gapped deconvolution in the tau-p domain based on the multiple lag, and (<NUM>) reverse tau-p transform the processed data to obtain the seismic data in the time-space domain. This method is proven to work for two or more seismic data sets that were recorded with same spatial sampling. However, for the case of different spatial sampling for low- and high-frequency data, the tau-p model may be designed with a linear inversion so that when the reverse tau-p transform is applied, it simultaneously satisfies both low- and high-frequency data.

In other words, for standard linear Radon, the step of deriving the tau-p model may consist of solving the following linear equation: <MAT>.

Once the tau-p domain has been calculated, deconvolution may be applied, followed by the multiple model being reverse transformed back to the offset-time domain for obtaining the low- and high-frequency seismic data sets.

Other options may include applying the process to LF and HF data independently, or combining the LF and HF data (e.g., sum) first and then applying the process to the combined data. The summing of data with different spatial positions may involve selecting the nearest LF trace to each HF trace and summing them together. An LF trace may be summed with more than one HF trace. The data may be subsequently separated.

A case of tau-p deconvolution in the shot domain with LF and HF sources as shown in <FIG> is now discussed to illustrate this method. <FIG> is a flowchart of a method for applying gapping deconvolution. In step <NUM>, two seismic data sets are received, one from source d1 and the other one from source d2. Each set of seismic data includes both LF and HF data. The data from source d1 is illustrated in <FIG> and the seismic data from source d2 is illustrated in <FIG> illustrates the HF data from source d1 and the LF data from source d2, where a <NUM> frequency was used for the low-cut and high-cut of the full data from each source.

In step <NUM>, the tau-p deconvolution based on equations (<NUM>) and (<NUM>) is applied to (a) the full data received from source d1, the results of which are illustrated in <FIG>, (b) the full data received from source d2, the results of which are illustrated in <FIG>, and (c) the combined data (i.e., HF data from source d1 and LF data from source d2 as illustrated in <FIG>, with the two data sets recorded with different spatial sampling), the results of which are illustrated in <FIG>.

In step <NUM>, the LF data (see <FIG>) is separated from the HF data (see <FIG>), for the combined tau-p deconvolved data illustrated in <FIG>. Then, to perform the comparison of the deconvolved combined data, which includes different spatial sampling, the LF data shown in <FIG> is subtracted in step <NUM> from the LF data corresponding to the deconvolved full data of source d2 (see Figure E), resulting in the residual shown in <FIG>. Similarly, the HF data shown in <FIG> is subtracted in step <NUM> from the HF data corresponding to the deconvolved full data of source d1 (see <FIG>), resulting in the residual shown in <FIG>.

The consistent demultiple results of <FIG>, along with low amplitude differences in <FIG>, confirm the equivalence of the tau-p deconvolution demultiple results of the combined data set illustrated in <FIG> compared to the independent data sets illustrated in <FIG>.

Another embodiment is now discussed with regard to data regularization. A traditional method is 3D data regularization using a guided anti-leakage Fourier transform (see for more details, <NPL>). This particular method of the Anti-leakage Fourier transform (ALFT) algorithm uses a low-frequency model of the data to define the wavenumber decomposition order for the higher temporal frequencies. However, this method was designed to work for low- and high-frequency data recorded with the same spatial sampling. In this embodiment, the modified version of the algorithm uses low-frequency data that may come (at least partially) from traces which have different spatial positioning to the high-frequency data. For example, it is possible to use the same acquisition scenario as in the embodiment illustrated in <FIG>. For the present embodiment, the data is considered in the offset-volume domain, within a spatial window in the midpoint-x and midpoint-y directions, as illustrated in <FIG>. The midpoint-x and the midpoint-y coordinates shown in <FIG> correspond to the coordinates of those points located at mid-distances between the source and the receiver. <FIG> shows low-frequency CMP <NUM> and high-frequency CMP <NUM>. Note that <FIG> indicates the different spatial sampling of the low- and high-frequency data.

In practice, it may be possible to have coarser sampling for the low-frequency data compared to the high-frequency data. The two data sets <NUM> and <NUM> both satisfactorily cover the same spatial area <NUM>. A novel regularization scheme for such data sets is now discussed with regard to <FIG>.

The method illustrated in <FIG> has a step <NUM> of receiving HF and LF data recorded in the time-space domain, i.e., HF(t,x,y) and LF(t,x,y). In step <NUM>, a Fast Fourier transform (FFT) is applied to the HF and LF data in the time-space domain, such that the HF(t,x,y) data is transformed into HF(f,x,y) in the frequency-space domain, and the LF(t,x,y) data is transformed into LF(f,x,y) in the frequency-space domain.

In step <NUM>, the irregular spatial Fourier transform for the LF data in the frequency-space domain is calculated, i.e., LF(f,x,y) is transformed to LF(f,kx,ky) in the frequency-wavenumber domain. In step <NUM>, the output LF(f,kx,ky) from step <NUM> is used to constrain the spatial Fourier transform for the HF data, i.e., to constrain the transformation of HF(f,x,y) to HF(f,kx,ky). In step <NUM>, the Fourier transform of step <NUM> may optionally be used to attenuate noise (either random or coherent) as scaling.

In step <NUM>, the spatial Fourier transform of HF and LF data may be used to generate data at new spatial positions. These spatial positions may be regularly sampled (e.g., at the center of a bin grid) or irregularly sampled (e.g., at the position of midpoints from a different data set). For example, the spatial Fourier transform may act on LF(f,kx,ky) in the frequency-wavenumber domain to obtain LF(f,x,y) in the frequency-space domain and similarly transform HF(f,kx,ky) to HF(f,x,y).

In step <NUM>, the data may be reverse Fourier transformed to the time-space domain, e.g., from LF(f,x,y) to LF(t,x,y) and from HF(f,x,y) to HF(t,x,y). If the output positions for HF and LF are consistent, then it is possible to combine in step <NUM> the data to form a data set covering the full bandwidth relating to HF and LF, e.g., F(x,y,t) = LF(x,y,t) + HF(x,y,t). In the case that there is a frequency overlap of the two data sets, a normalization weighting in the overlap zone may be required. Finally, in step <NUM>, existing processing techniques may be applied to generate an image of the surveyed subsurface. Note that this embodiment uses Fourier transforms to manipulate the data. Those skilled in the art would know that other transforms and/or combination of transforms may be used to achieve the same results.

With regard to step <NUM>, it may include a weighted summation to derive each wavenumber, as given by equation (<NUM>) <MAT> where wavenumbers in the x and y direction are given by Kj and Kk, and the input traces have positions xl and yl. The weights wI may be defined based on the Voronoi weighting and normalized to <NUM>.

With regard to step <NUM>, in one application, the data from step <NUM> is used to estimate energy amplitude at different dips for the spatio-temporal window. This may be achieved by mapping the FKxKy domain to the Fpxpy domain and summing the amplitudes in frequency. Once this estimate has been made, it may be propagated to all frequencies and used to derive the decomposition order for the anti-leakage transform.

The method described in the flowchart of <FIG> is now illustrated by using seismic data recorded with a dual-source towed streamer acquisition, similar to the embodiment illustrated in <FIG>. The data is in the common offset domain. The full bandwidth input data is illustrated in <FIG>. Using regular processing, the data is regularized as illustrated in <FIG>. This would be the reference data to be compared with the regularized data using the method illustrated in <FIG>. When applying the method of <FIG>, HF and LF data having different spatial sampling is used. This data is created as follows: the energy below <NUM> is zero for the port source to obtain the HF data having one spatial sampling (see <FIG>); the starboard source data below <NUM> is extracted to obtain the LF data having a different spatial sampling (see <FIG>). The HF data is then regularized using the LF data, i.e., the starboard source data below <NUM>, to regularize the HF data above <NUM> data for both starboard and port sources, as illustrated in <FIG>. The difference between the reference regularized data of <FIG> and the novel regularization scheme of <FIG> is shown in <FIG>. It can be seen that the difference is relatively minor.

This is true because, while in theory the results from the second regularization will not be as good as the first, in fact there is a negligible difference between the results. This is because low-frequency energy varies spatially very slowly. Therefore, even though low-frequency data is sampled more coarsely than high-frequency data, it is still sufficient as a guide for high-frequency regularization.

The previous embodiment related to 3D regularization. However, the method discussed above with regard to <FIG> is also applicable to regularization in higher dimensions (e.g., mid-x, mid-y, offset, azimuth, time, etc.). Further, although the ALFT algorithm was described in this embodiment, other schemes of deriving a model may be used, e.g., inversion, which derives a model to satisfy the input data with l2, l1, Cauchy, or other norm. Sparseness weights may be updated iteratively, e.g., iteratively re-weighted least squares. Sparseness weights from one or both data sets may be used to stabilize the inversion. In the case of inversion approaches, a model domain should be derived to simultaneously satisfy LF and HF data on their respective positions in a similar way to equation (<NUM>) given in the tau-p deconvolution section. The scheme may be applied in a frequency domain (Fourier, Laplace, Z-transform or other) or in the time domain.

Various model domains may be used for noise attenuation and/or data interpolation in the schemes previously discussed. Some examples of model domains are:.

Higher dimensional transforms using a single type or combination of types of the above-noted transforms are also possible.

While the above discussion has related to using LF to help build a dealiased model of HF, the energy used for a guide may not come exclusively from the LF data as described above. Likewise, the guide may not necessarily be used exclusively to derive the model for the HF. In fact, the guide may be derived partly from LF and partly from HF data, and may subsequently be used to derive the model for LF and/or HF data. Alternatively, the guide may be derived from LF, but used for model derivation of LF and/or HF.

As previously discussed in this application, two seismic data sets having different spatial sampling may be applied to the imaging methods. A few of these imaging methods are now discussed. In general, the traditional approach is to combine LF and HF data before migration and to migrate the combined data as per conventional data. However, according to the embodiments next discussed, the LF and HF datasets could be migrated separately and combined after migration. While the input sampling to the migration may be at different positions, the output data will be co-located and can be summed. Some of these alternative approaches are outlined below.

Kirchhoff migration (either time or depth domain), from a simplistic point of view, may be considered to include a step of travel-time calculation from a velocity estimate of the subsurface, a step of applying a migration operator including anti-alias filtering as necessary, and a step of applying the migration operator to the input data. Given an accurate velocity model, constructive and destructive interference of the migration operator results in an accurately positioned representation of the subsurface. Because the Kirchhoff operator is orthogonal in the temporal frequency sense, it is a valid approach to migrate the LF data and HF data separately and sum the migration results. A variation to the above-noted scheme would be to apply appropriate low/high-cut filtering to the LF and HF data so that when they are combined, there is no undesirable effect at the cross-over frequency between the LF and HF data.

In one embodiment, instead of migrating LF and HF data separately, it is possible, and sometimes beneficial, to regularize both data to the same sampling before migration, as discussed previously with regard to <FIG>.

If beam migration is considered, it consists of deriving slant stack representations of common-midpoint spatial window data in the common offset (or common offset-x/offset-y) volumes. While the HF and LF data may be combined prior to imaging with data regularization as discussed above with regard to <FIG>, it is also possible to combine HF and LF data during the step of forming the slant stack data representation prior to migration. This approach may consist of pre-processing LF and HF data to avoid cross-over frequency issues (e.g., scaling energy at a given frequency which has contributions from each dataset so that it is consistent with a source emission with flat amplitude spectrum), followed by the step of slant stack. The spatial sampling of LF and HF data may be different, resulting in different normalization in the slant stack procedure for LF and HF data.

Still another imaging process is based on the wave-equation migration (e.g., one-way or reverse time migration). Wave equation migration consists of a source wave field (beginning with the source signature) and a receiver wave field (the one recorded by the receiver). These wave fields are propagated in time and cross-correlated to make the output image through use of an imaging condition (e.g., cross-correlation).

<FIG> schematically illustrate the wave-equation migration process using two data sets of different spatial sampling. A one-way wave equation source wave field <NUM> (see <FIG>) is being forward-propagated into the subsurface and the receiver wave field <NUM>, <NUM> (see <FIG>) is being backward-propagated back into the subsurface. When energy correlates in the source and receiver wave fields, the cross-correlation adds energy <NUM> into the output image (see <FIG>). Note that <FIG> shows two curves <NUM> and <NUM> that correspond to detected events reflected from two subsurface layers. For this simple case, there are two horizontal reflectors at <NUM> and <NUM> depth, and the wave fields are considered at depth zero, the depth of the first reflector, and the depth of the second reflector. It can be seen that at <NUM> depth, recorded energy <NUM> relating to the first reflector in <FIG> coincides with the energy <NUM> in <FIG> of the source, giving a strong correlation <NUM> in <FIG>. At the depth of the second reflector (<FIG>), the energy <NUM> of the second arrival in <FIG> coincides with the source wave field <NUM> in <FIG>, resulting in a strong correlation <NUM> in <FIG>.

The process may be repeated for all the gathers in the data set, the result from each being summed to generate the final image. In the case that the LF and HF data relate to sources of different bandwidths, they may be migrated independently and combined post-migration as in the discussion for Kirchhoff migration.

In the case that the HF and LF data relate to different receiver types, the data may be migrated in the receiver domain.

The above scheme may be modified to cover the case where more than one source is simultaneously excited. This embodiment is illustrated in <FIG> shows two source wave fields <NUM> and <NUM> leaving the surface simultaneously and interfering at <NUM>. This arrangement generates two pairs of reflectors on the receiver side, as illustrated by fields <NUM>, <NUM>, <NUM> and <NUM> in <FIG>. Even though the receiver wave field is now more complex, there is still an exact match at the reflector depth as before, as illustrated by energies <NUM> and <NUM> in <FIG>. The only problem with this scheme relates to interference between the two sources, e.g. , where energies <NUM> and <NUM> overlap with energies <NUM> and <NUM> in <FIG>.

The analysis presented with regard to this embodiment may be extended to cover the case where the second source is excited with a delay relative to the first. This example is shown in <FIG> illustrates the time delay between shooting the first and second sources. While the timing of the events has changed, there is still consistency in timing between the source and receiver sides at the depths of the reflectors, as illustrated by energies <NUM> and <NUM>, corresponding to the two reflectors at <NUM> and <NUM> depth.

According to another embodiment, in the case of moving sources emitting energy continuously or semi-continuously (e.g., marine vibrators or de- synchronised air gun source), <FIG> may be modified to be based on continuous recording time rather than listening time. This is essentially a simultaneous shooting scenario where energy relating to earlier times overlaps with energy relating to later times, following their propagation through the subsurface. In this regard, <FIG> shows the source data and <FIG> shows the receiver data based on surface datum. In this case, the source and receiver data have been mapped on to a regular grid using spatial interpolation (e.g., sinc). The gray area <NUM> in <FIG> highlights the zone of interference based on the direct arrival of the moving source.

<FIG> illustrates a source emitting a random signal while moving and <FIG> illustrates four receivers continuously recording while moving behind a vessel. Instead of a given listening time, e.g., <NUM> seconds, a continuous recording time, e.g., <NUM> hours, is used in this embodiment while the vessel travels for the whole line of a 3D acquisition. Instead of migrating many independent shots and summing the results together, the method propagates the wave fields based on continuous recording time and derives an image for the whole sail line in one go. According to this method, a lot of interference is expected to occur, but if a pseudo-random or pseudo-orthogonal signal is used, the interference should not cause problems for the imaging condition. The above-noted method may be extended to the case where more than one source with continuous emission is used, as illustrated in <FIG>.

While the examples discussed above with regard to <FIG> relate to one-way wave equation migration, the same principle may be used for other schemes, e.g., RTM (reverse time migration). While the above example was in 2D, the same principle extends to 3D (e.g., narrow, wide, multi-azimuth, etc.).

According to another embodiment, seismic data sets having different spatial sampling may be used with a least squares migration method. One aim of a traditional least squares migration is to derive a migrated image of seismic data so that when reverse migrated, it matches the input data (e.g., in a least squares sense). This traditional principle may apply to all types of migration (e.g., Kirchhoff, beam or wave equation). An example of traditional least squares RTM may be found in <NPL>. The migration algorithm may be designed to migrate primary energy, multiples, or a combination of the two.

In the case of seismic data having different bandwidths and/or positioning, according to an embodiment, it is possible to define a novel least squares migration problem as follows: derive a subsurface image so that when reverse migrated, it satisfies the LF data for the LF bandwidth, and it also satisfies the HF data for the HF bandwidth. In practical terms, as illustrated in <FIG>, two seismic data sets are received in step <NUM>, and this seismic data is pre-processed in step <NUM> for source designature. In step <NUM>, the method reverse migrates an image of the subsurface to generate data in the time domain (sx, sy, rx, ry, t) for all recorded data positions. The image may be calculated by a previous iteration of the method as discussed later or it may be based on one of the models discussed above. In one embodiment, the image is selected based on the experience of the operator or based on previous surveys of the subsurface. In step <NUM>, the method applies a high-cut filter to the seismic data relating to the positions of the LF data, and in step <NUM> the method applies a low-cut filter to the data relating to the positions of the HF data.

In step <NUM>, a cost function (e.g., the sum of squares of the difference between modeled and recorded energy) is jointly calculated based on the difference between the LF data of step <NUM> and the HF data of step <NUM>. The cost function may have various configurations as will be appreciated by those skilled in the art. Then, in step <NUM> the method compares the result of the cost function with a given threshold. If the result of the comparison is above the threshold, the method returns in step <NUM> to update the migration image selected in step <NUM> and, in an iterative fashion, recalculates another image until the image which represents the input data as accurately as possible. When this process is performed, the method generates in step <NUM> an image of the subsurface.

In the case that the HF and LF data contain source signature effects, steps <NUM> and <NUM> may be modified to include convolution with the source wavelet. This convolution may be 1D, 2D or 3D.

In the case that a blended acquisition using HF and LF sources has been used, the modeled HF and LF data (steps <NUM> and <NUM>) after convolution with the relevant signatures may be combined (e.g., summed). The input data for this scheme may be based on continuous recording time or listening plus source signature time (listening time may be defined as the time for all useful reflected signals to have been recorded).

This approach is relevant for continuous source emission or "sweep then listen" acquisition. In the case of continuous emission, this means that an image needs to be found which, when reverse migrated based on the moving source and receiver positions, reproduces the input data in the same format as <FIG> as accurately as possible.

The above-discussed embodiments illustrated how to process two seismic data sets having different spatial sampling for various processing steps, e.g., noise removal, data regularization, imaging, etc. Based on any of the above-discussed embodiments, a method for processing seismic data associated with a surveyed subsurface, is now discussed with reference to <FIG>. The method includes a step <NUM> of receiving at least first and second seismic data sets that were recorded with different spatial sampling and different temporal bandwidths, a step <NUM> of using the first seismic data to guide a processing of the second seismic data, and a step <NUM> of generating an image of the subsurface based on processed second seismic data.

According to another embodiment, not covered by the claims, there is a method for processing seismic data associated with a surveyed subsurface. The method includes a step of receiving at least first and second seismic data sets that were recorded with different spatial sampling and different temporal bandwidths, a step of interpolating the first seismic data to positions of the second seismic data to obtain interpolated first seismic data, a step of combining the interpolated first seismic data with the second seismic data, and a step of generating an image of the subsurface based on a combination of the interpolated first seismic data with the second seismic data.

The different temporal bandwidths of the first and second seismic data sets relate to low-frequency data and high-frequency data. In one embodiment, the first and second seismic data sets relate to energy generated by a vibrating source that operates on land or water. The combination of the interpolated first seismic data and the second seismic data may be output on a regular grid or an irregular grid.

To summarize the above discussed embodiments, it is possible to use the multiple datasets as follows:.

Seismic data recorded with different spatial sampling and different temporal bandwidths as discussed above may be processed in a corresponding processing device for generating an image of the surveyed subsurface as discussed now with regard to <FIG>. For example, the seismic data collected with the characteristic noted above may be received in step <NUM> at the processing device. In step <NUM>, pre-processing methods are applied, e.g., demultiple, signature deconvolution, trace summing, motion correction, vibroseis correlation, resampling, etc. In step <NUM>, the main processing takes place, e.g., deconvolution, amplitude analysis, statics determination, common middle point gathering, velocity analysis, normal move-out correction, muting, trace equalization, stacking, noise rejection, amplitude equalization, etc., as discussed in the previous embodiments. In step <NUM>, final or post-processing methods are applied, e.g., migration, wavelet processing, seismic attribute estimation, inversion, etc., and in step <NUM> the final image of the subsurface is generated. Note that at least one process in one of steps <NUM> and/or <NUM> involves dealing with two or more seismic data sets having different spatial sampling and different temporal bandwidths. Also note that a result of this processing may result in an improvement of the image of the surveyed subsurface. The improved image quality will result in an increased likelihood of finding the oil and gas reservoirs.

An example of a representative processing (or control) device capable of carrying out operations in accordance with the embodiments discussed above is illustrated in <FIG>. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein. The processing device <NUM> of <FIG> is an exemplary computing structure that may implement any of the processes and methods discussed above or combinations of them.

The exemplary processing device <NUM> suitable for performing the activities described in the exemplary embodiments may include server <NUM>. Such a server <NUM> may include a central processor unit (CPU) <NUM> coupled to a random access memory (RAM) <NUM> and/or to a read-only memory (ROM) <NUM>. ROM <NUM> may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor <NUM> may communicate with other internal and external components through input/output (I/O) circuitry <NUM> and bussing <NUM> to provide control signals and the like. For example, processor <NUM> may communicate with appropriate valves of the source elements for controlling the air pressure inside each source element. Processor <NUM> carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

Server <NUM> may also include one or more data storage devices, including disk drives <NUM>, CD-ROM drives <NUM>, and other hardware capable of reading and/or storing information, such as a DVD, etc. In one embodiment, software for carrying out the above-discussed steps may be stored and distributed on a CD-ROM <NUM>, removable media <NUM> or other form of media capable of storing information. The storage media may be inserted into, and read by, devices such as the CD-ROM drive <NUM>, disk drive <NUM>, etc. Server <NUM> may be coupled to a display <NUM>, which may be any type of known display or presentation screen, such as LCD, plasma displays, cathode ray tubes (CRT), etc. A user input interface <NUM> is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc..

Server <NUM> may be coupled to other computing devices, such as the equipment of a vessel, via a network. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet <NUM>, which allows ultimate connection to various landline and/or mobile client/watcher devices.

As also will be appreciated by one skilled in the art, the exemplary embodiments may be embodied in a wireless communication device, a telecommunication network, as a method or in a computer program product. Accordingly, the exemplary embodiments may take the form of an entirely hardware embodiment or an embodiment combining hardware and software aspects. Further, the exemplary embodiments may take the form of a computer program product stored on a computer-readable storage medium having computer-readable instructions embodied in the medium. Any suitable computer-readable medium may be utilized, including hard disks, CD-ROMs, digital versatile discs (DVD), optical storage devices or magnetic storage devices such a floppy disk or magnetic tape. Other non-limiting examples of computer-readable media include flash-type memories or other known types of memories.

The disclosed exemplary embodiments provide a device and method for processing two seismic data sets having different spatial sampling and different temporal bandwidths. It should be understood that this description is not intended to limit the invention. On the contrary, the exemplary embodiments are intended to cover alternatives, modifications and equivalents, which are included in the scope of the invention as defined by the appended claims. Further, in the detailed description of the exemplary embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

Although the features and elements of the present exemplary embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein.

Claim 1:
A method for processing seismic data associated with a surveyed subsurface, the method comprising:
receiving (<NUM>) at least first and second seismic data sets that were recorded with different temporal bandwidths, wherein the first and second seismic data sets respectively relate to energy generated by two vibrating sources that operate on land or water, the different temporal bandwidths of the first and second seismic data sets relating respectively to low-frequency data and high-frequency data;
using (<NUM>) the first seismic low-frequency data to guide a processing of the second seismic high-frequency data; and
deriving a model which simultaneously satisfies the first and second seismic data sets; and
generating (<NUM>) an image of the subsurface based on processed second seismic data, the image being based on the model, the different temporal bandwidths relating to bandwidths overlapping by less than half of the bandwidth in octaves emitted by the vibrating source,
the method being characterized in that:
the distance between the two vibrating sources is set at <NUM> or less; and
the at least first and second seismic data sets were recorded with different spatial sampling.