MULTIPLE ATTENUATION AND IMAGING METHODS FOR RECORDED SEISMIC DATA
TECHNICAL FIELD
Embodiments of the subject matter disclosed herein generally relate to processing seismic data acquired over land or water and, more specifically, to a method for imaging of a surveyed subsurface formation.
In the following, the term “formation” refers to any geophysical structure into which source energy is promulgated to perform seismic surveying, e.g., land or water based, such that a “formation” may include a water layer when the context is marine seismic surveying.
BACKGROUND
Hydrocarbon exploration and development uses waves (e.g., seismic waves or electromagnetic waves) to explore the structure of underground formations on land and at sea (i.e., formations in the subsurface). These waves are collected during a seismic acquisition campaign, either on land or on water. In such seismic acquisition campaigns, several seismic sources and sensors are used. Energy generated by a seismic source propagates as seismic waves downward into a geological formation, and part of the energy is reflected and/or diffracted and/or refracted to the receivers. Characteristics of the reflected/refracted/diffracted energy detected by seismic sensors are used to produce an image of the earth’s subsurface reflectivity. While this profile does not provide an accurate location for oil and gas, it suggests, to those trained in the field, the presence or absence of oil and/or gas and/or to explore the structure of underground formations on land and at sea (i.e., formations in the subsurface). Thus, providing a high-resolution image of the subsurface is an ongoing process for the exploration of natural resources, including, among others, oil and/or gas.
As schematically shown in Figure 1, a marine seismic acquisition system 100 is used to image the subsurface. This system comprises a source 110 that emits waves 112 at a known location that penetrate an unexplored formation 120 and are reflected at plural interfaces 122, 124, 126. The system comprises further a plurality of sensors or receivers 130, which are distributed along a cable 132 referred to as a streamer. Said streamer is towed by a vessel 102, detect the reflected waves 140, 150. The detected waves 140, 150 include primary reflections such as wave 140, which travels directly from a formation interface to a sensor, and multiple reflections such as wave 150, which additionally undergoes at least one more downward reflection and one upward reflection. In the case of wave 150, the downward reflection is at the water surface and the upward reflection is at reflector 122.
There are various types of multiples, e.g., surface-related multiples and interbed multiples. In the case of surface-related multiples, when the energy reflected from the subsurface reaches the water surface 104, it will be reflected back downwards into the water column and subsurface. This produces a second set of reflected energy containing spurious events. Interbed multiples are similar, but in this case the downward reflecting surface occurs at a rock interface in the subsurface.
Moreover, multiples can be characterized as belonging to different orders, e.g., first order, second order, third order, etc., based on the number of additional reflections involved. For example, a primary P has a single reflection between a source S and a receiver R as shown in Figure 2A. By way of contrast, a first order multiple M1 (shown in Figure 2B) can have two additional reflections relative to the primary P, whereas a second order multiple M2 (shown in Figure 2C) can have four additional reflections.
The use of efficient wide-tow streamer acquisitions limits the availability of shortoffset primary arrivals for shallow imaging. This results in primary imaging with strong acquisition striping and poor vertical resolution, making interpretation of the shallow section difficult. While data interpolation and regularization can help reduce striping, vertical resolution in the shallow section may still be problematic.
Acquisition approaches to alleviate this problem include reducing the sail-line spacing, placing the sources closer to (Dhelie et al., 2020) or above (Vinje et al., 2017) the streamers, and nearfield hydrophone (NFH) imaging (Davies and Tillotson, 2019). However, these approaches come with increased acquisition cost, or incomplete coverage in the case of NFH imaging for multi-streamer surveys.
The imaging of free-surface multiples offers a data-driven alternative which has been discussed in the literature for many years. While the source-side wavefield for primary imaging propagates from a single shot location, the down-going wavefield used in multiple imaging utilizes all receivers as secondary sources, thus significantly increasing illumination from the surface (Berkhout and Verschuur, 1994). This concept was developed further by (Whitmore et al. (2010)) who used up-going and down-going wavefields from a dual-sensor towed-streamer.
However, one complication arising from multiple imaging is cross-talk, where different orders of surface scattering can contaminate the image. Lu et al. (2016a) describes a practical approach to attenuate this noise based on predicting and subtracting causal and anti-causal cross-talk. The use of least-squares multiple imaging has also been documented as a way to attenuate cross-talk (for example, Zhang and Schuster, 2014, or Lu et al., 2018).
Cross-talk is a well-documented side-effect of multiple imaging, being exhibited by spurious events in the image of the subsurface with a similar appearance to multiples. Causal cross-talk occurs whenever reverse extrapolated up-going multiple energy is extrapolated beyond its actual reflection point and interacts with the extrapolated downgoing wavefield. With anti-causal cross-talk, this effect occurs before its actual reflection point. A discussion of causal and anti-causal cross-talk can be found in Lu et al. (2016a).
Some imaging conditions have been documented to behave differently to others in reducing cross-talk effects. The deconvolution imaging condition has been well documented to reduce cross-talk relative to the cross-correlation imaging condition (Whitmore et al., 2010, and Poole et al., 2021). The 2D deconvolution imaging condition of Valenciano and Biondi (2003) has been shown to further reduce cross-talk by deconvolving plane waves across several shots (Mujis et al., 2005, and Mujis et al., 2007). The multidimensional imaging condition solves a least squares problem at each depth level (Valenciano et al., 2002).
Lu et al. (2016a) describes ways of estimating causal and anti-causal cross-talk in the image domain. Causal cross-talk may be estimated by migrating a down-going source wavelet, with an up-going multiples wavefield. Anti-causal cross-talk may be estimated through migration of the down-going recorded data with up-going primaries. Simple approaches to estimate a cross-talk image corresponding to water-bottom contamination may involve vertical depth shifting the image of the subsurface down to the depth of the cross-talk, and adaptively subtracting.
Figure 3A shows multiple imaging using a cross-correlation imaging condition. The image is low resolution and contains heavily reverberating energy i.e cross-talk (highlighted by the arrows).
Figure 3B shows multiple imaging of the same data using the deconvolution imaging condition where the result is higher in resolution. Although less prevalent than in the cross-correlation result, some cross-talk is still evident in this deconvolution image, as shown with arrows.
Figure 3C shows multiple imaging of the same data using a least-squares migration approach. Figure 3C shows the least-squares multiple migration image which contains less cross-talk than either the cross-correlation or deconvolution approach. Never-the-less, cross-talk is still visible on this least-squares migration result.
It is desirable to have a method for imaging the surveyed that obviates all or part of the abovementioned drawbacks, in particular, to reduce the cross-talk effects.
SUMMARY OF THE INVENTION
Exemplary embodiments of the invention aim to satisfy this need and relate to a method implemented by a computer for imaging a surveyed subsurface formation, the method comprising:
a) Receiving a recorded seismic dataset associated with the subsurface, b) Estimating a multiple reflection dataset corresponding to a defined subsurface reflector,
c) Generating a partial demultiple dataset by subtracting the multiple reflection dataset from the recorded seismic dataset,
d) Receiving an earth model,
e) Generating a synthetic dataset based on a forward propagation through the earth model,
f) Calculating a residual dataset based on the partial demultiple dataset and the synthetic dataset,
g) Updating the earth model based on the residual dataset, and
h) Optionally, generating a final image of the surveyed subsurface based on the updated earth model.
Advantageously, the final image has reduced cross-talk levels compared to an image generated to the same output depth using the recorded dataset.
Said final image of the subsurface may be used to predict multiples.
In step b), multiple reflection dataset may be based on a multiple modelling. The multiple modelling may be for example based on a convolution or on a wave propagation through an earth model containing said defined subsurface reflector.
The residual dataset may be calculated in the time-space domain. Preferably, the residual dataset is back propagated using said earth model, and this back propagated residual dataset is used to update the earth model.
The residual dataset may be in the earth model domain. Preferably, the earth model is updated by back propagating the synthetic dataset through the earth model.
The back propagation may be one of a two-way propagation and a one-way propagation.
In step g), updating the earth model may be based at least in part on imaging or full waveform inversion to improve a misfit between synthetic data and recorded data. The imaging of full waveform inversion may be based on illumination of the subsurface using multiples.
The step of updating the earth model may be based on a cross-correlation imaging condition, a deconvolution imaging condition, an optimized migration, or a FWI.
In step e), the forward propagation may be initiated by the injection of a source wavelet into the earth model.
In step e), the forward propagation may be initiated by the injection of the recorded data into the earth model.
In step d), the earth model may be one or more of a velocity model, reflectivity model, density model, impedance model, absorption model, or an anisotropy model.
In step b), the multiple reflection dataset may comprise a first multiple dataset and/or a second multiple dataset, the first multiple dataset being generated using the recorded dataset and the second synthetic dataset being generated using the first multiple seismic dataset.
The first multiple dataset may comprise first and higher order multiples and the second multiple seismic dataset may comprises second and higher multiples.
The multiple dataset, notably the first and/or second multiple seismic dataset, may be computed using:
- wavefield extrapolation;
- model based multiple modelling, notably based on a horizon interpreted from the first image of the subsurface; or
- demigration of the first image of the surveyed subsurface.
The subtraction in step c) may be an adaptive subtraction of at least one of the following:
- the first multiple dataset ;
- the second multiples dataset.
The method may comprise generating:
- a first image of the surveyed subsurface using the recorded seismic and at least in part on a multiple imaging method, and
-a combined image by combining a shallow depth range from the first image with a deep depth range of the final image generated in step h).
The maximum depth of the deep depth range may be at least 0.5 times, preferably at least 1 time, more preferably at least 1.5 times, the maximum depth of the shallow depth range.
The invention also relates to a method implemented by a computer for multiple imaging a surveyed subsurface formation, the method comprising:
a) Providing a first recorded seismic dataset and a second recorded seismic dataset associated with the subsurface, wherein the first and second recorded datasets comprise primaries and/or multiples,
b) Generating a forwards propagated dataset by forward propagating the first recorded seismic dataset,
c) Generating a backward propagated dataset by backward propagating the second recorded seismic dataset,
d) Generating a first image of the surveyed subsurface, based on the forwards and backwards propagated datasets and at least in part on a multiple imaging method,
e) Computing a synthetic seismic dataset based on the first image of the surveyed subsurface,
f) Generating a modified second seismic dataset by subtracting the synthetic seismic dataset from the second recorded seismic dataset,
g) Generating a modified backwards propagated dataset by backward propagating the modified second seismic dataset,
h) Generating a second image of the surveyed subsurface based on the modified backwards propagated dataset.
Optionally, generating a final image of the surveyed subsurface by combining a shallow depth range from the first image of the subsurface with a deep depth range of the second image of the subsurface. Preferably, the maximum depth of the deep depth range is at least 0.5 times, preferably at least 1 time, more preferably at least 1.5 times, the maximum depth of the shallow depth range.
Advantageously, the second image will have reduced cross-talk levels compared to an image generated to the same output depth using the first and second datasets.
Hence, the invention takes advantage of multiples while having reduced cross-talk levels.
The method according to the invention may relate to multiple imaging or FWI using multiples corresponding to wave-propagation in a shot domain or in a receiver domain. Typically, shot domain applications may relate to towed streamer acquisitions, and receiver domain applications may relate to OBN (ocean bottom node) acquisitions.
The proposed method makes it possible to use any multiple imaging technique at steps d) and h). The multiple imaging method may involve for example a cross-correlation imaging condition, a deconvolution imaging condition, or an optimized migration.
In step e), the synthetic seismic data generated may comprise a first synthetic dataset and/or a second synthetic dataset.
The first synthetic dataset may be generated using the first recorded dataset and the second synthetic dataset being generated using the first synthetic seismic dataset.
In some embodiments, the synthetic seismic dataset is generated using the first recorded seismic and the first image derived from step d).
The first synthetic dataset may comprise more than one order of multiples, notably first and higher order multiples.
The second synthetic seismic dataset may comprise second and higher order multiples.
The first and/or second synthetic multiple seismic datasets may be computed based on several methods.
For example, wavefield extrapolation of the recorded seismic dataset may be used. Seismic data extrapolation relates to a method that simulates recordings of a wavefield at a position other to which it was recorded. For example, in horizontal towed streamer acquisition, one may record hydrophone measurements at a depth of 12 m. These measurements may be extrapolated to a new datum, for example, at 100 m depth.
The extrapolation may use an estimate of the subsurface properties, e.g. velocity, anisotropy, absorption, etc.
The extrapolation may be with one-way or two-way extrapolations and may be forward or reverse in direction, more information is for example given in (Biondi, 2006).
The first and second synthetic datasets may be generated using a model based multiple modelling based on a horizon interpreted from the first image of the subsurface, as for example disclosed in (Wang et al., 2011) or (Wiggins, 1988).
In some embodiments, the first and second synthetic datasets may be generated using a demigration of the first image of the surveyed subsurface to generate a primary model for SRME modelling.
In some embodiments, the subtraction in step f) is an adaptive subtraction of at least one of the following:
- the first synthetic dataset;
- the second synthetic dataset.
Preferably, step f) comprises the subtraction of both first and second synthetic datasets.
The first recorded seismic data may comprise at least one of primaries, multiples and direct arrivals.
The second recorded seismic data may comprise at least one of multiples and primaries.
The first and second datasets may be the same or different.
For example, towed streamer deghosted data may be used for the first and second datasets.
Ocean bottom node (OBN) data may be used for the first and second datasets. The invention also relates to a method implemented by a computer for imaging a surveyed subsurface formation using full-waveform inversion (FWI), the method comprising:
- receiving recorded data,
- deriving partial demultiple data from said recorded data,
- determining a structural model of the subsurface by applying a full waveform inversion, FWI, algorithm, wherein the input to the FWI modelling being based on a source signature or an areal source corresponding to the injection of recorded data, and,
- optionally, generating an image of the subsurface formation based on said structural model.
The partial demultiple data may be obtained from the recorded data by attenuating the primary reflection and/or multiples corresponding to a defined subsurface reflection.
The primary reflection may be attenuated by muting, or subtraction of modelled primary arrivals.
The multiples may be attenuated using a multiple prediction and multiple attenuation method (e.g. SRME (Berkhout and Verschuur, 1997), partial SRME (Hugonnet, 2002), SRMM (Pica et al., 2005), MWD (Wang et al., 2011), wave-equation targeted multiple modelling (Wiggins, 1988), deconvolution (Biersteker, 2001), wave-equation deconvolution (Poole, 2019)).
The invention also relates to a computing device for implementing the method of imaging a surveyed subsurface according to the invention, the computing device comprising:
- an interface configured to receive the recorded seismic dataset associated with the surface; and
- a processor connected to the interface, the processor configured to: ● Receive a recorded seismic dataset associated with the subsurface, ● Estimate a multiple reflection dataset corresponding to a defined subsurface reflector,
● Generate a partial demultiple dataset by subtracting the multiple reflection dataset from the recorded seismic dataset,
● Receive an earth model,
● Generate a synthetic dataset based on a forward propagation through the earth model,
● Calculate a residual dataset based on the partial demultiple dataset and the synthetic dataset,
● Update the earth model based on the residual, and
● Optionally, generate an image of the surveyed subsurface based on the updated earth model.
The invention also relates to a computing device for implementing the method of imaging a surveyed subsurface according to the invention, the computing device comprising:
- an interface configured to receive the first and second recorded seismic data; and
- a processor connected to the interface, the processor configured to: ● Generate a forwards propagated dataset by forward propagating the first recorded seismic dataset,
● Generate a backward propagated dataset by backward propagating the second recorded seismic dataset,
● Generate a first image of the surveyed subsurface, based on the forwards and backwards propagated datasets and at least in part on a multiple imaging method,
● Compute a synthetic seismic dataset based on the first image of the subsurface,
● Generate a modified second seismic dataset by subtracting the synthetic seismic dataset from the second recorded seismic dataset, ● Generate a modified backwards propagated dataset by backward propagating the modified second seismic dataset, ● Generate a second image of the surveyed subsurface based on the forwards propagated dataset and the modified backwards propagated dataset,
● Optionally, generating a final image of the surveyed subsurface by combining a shallow depth range from the first image of the subsurface with a deep depth range of the second image of the subsurface.
Preferably, the maximum depth of the deep depth range is at least 0.5 times, preferably at least 1 time, more preferably at least 1.5 times, the maximum depth of the shallow depth range.
BRIEF DESCRIPTION OF THE DRAWINGS
The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate presently preferred embodiments of the invention, and together with the description, serve to explain the principles of the invention. In the drawings:
Figure 1 is a schematic diagram illustrating seismic data that includes primaries and multiples;
Figure 2A, 2B and 2C show a primary, first order multiples and second order multiples, respectively;
Figure 3A illustrates multiple imaging resulting from cross-correlation imaging condition using recorded data, Figure 3B illustrates multiple imaging from deconvolution imaging condition using recorded data, and Figure 3C illustrates multiple imaging using Least-squares migration using recorded data;
Figure 4 schematically illustrates a wavefield that is forward propagated from a source, multiplied with a reflectivity at a given datum, and then being propagated to a receiver;
Figure 5 schematically illustrates the forward propagation of a source and reverse propagation of a primary wavefield recorded at a receiver to a given datum;
Figure 6 schematically illustrates the forward propagation of a primary wavefield and the reverse propagation of a multiple wavefield to a given datum, for example, the ocean bottom;
Figure 7 illustrates a flow chart of a method for multiple imaging according to the invention,
Figure 8A illustrates multiple imaging of waterbottom using recoded data, Figure 8B illustrates multiple imaging of data after removing multiples, and Figure 8C illustrates the combination of figures 8A and 8B;
Figure 9 shows a flow chart of a method for FWI imaging according to the invention, and
Figure 10 is a schematic diagram of a computing device configured to implement any of the above method.
DETAILED DESCRIPTION
In accordance with the invention, and as broadly embodied, a method according to the present invention is provided in Figure 7, aiming at generating an image surveyed subsurface formation.
First, in step 701, a first and second recorded seismic data are provided. The first and second recorded seismic dataset may be derived from the same recorded seismic dataset.
Alternatively, the first and second recorded seismic datasets are derived from different recorded seismic datasets.
The recorded data may be recorded by any sensor type, examples comprise hydrophones, geophones, particle motion sensors, particle velocity sensors, accelerometers, near-field hydrophones, near-field accelerometers or other sensor configured to detect seismic energy.
The recorded data may be from towed streamer, ocean-bottom sensor (node or cable) acquisition, land acquisition, transition zone campaign, or borehole (e.g. VSP, DAS).
The recorded data may be recorded at a constant datum, or a variable datum (e.g. variable depth streamer, or varying topography in land, e.g. floating datum). In the case of land data, the geophone recordings may be propagated forwards and backwards to form the multiple image, and the multiple image may be used to predict multiples.
In general, the terms OBN (ocean bottom node), OBC (ocean bottom cable), OBS (ocean bottom survey/sensor), and PRM (permanent reservoir monitoring) systems may be used interchangeably.
OBN receiver gather data may correspond to hydrophone, geophone (x, or, y, or, z), receiver-upgoing, receiver-downgoing, etc.
The recorded data may be pre-processed prior to imaging. The method may comprise, notably previous to step 701, a pre-processing step.
Said pre-processing step may involve noise attenuation (e.g. swell noise), guided wave attenuation, deblending, source signature compensation, source and/or receiver deghosting, wavefield separation, demulitple, data interpolation (e.g. based on Spitz, 1991) / regularisation (e.g. following Xu et al., 2010 for for down-going and/or up-going wavefields), redatuming (e.g. to free surface) and other approaches known in the art. Deghosting and/or wavefield separation may be on the source side or receiver side.
The first and second recorded seismic datasets may be derived by applying wavefield separation to the recorded data using a processor, at a locus of the depth receiver, to obtain a down-going wavefield corresponding to the first recorded seismic dataset and an up-going wavefield, constituting the second recorded seismic dataset.
The first recorded seismic data may comprise at least one of primaries, multiples and direct arrivals.
The second recorded seismic data may comprise at least one of multiples and primaries.
A selection of multiples have been schematically illustrated in Figures 1 to 2C. Primary arrivals P correspond to a signal leaving the source S, interacting once with the subsurface (e.g., reflection, refractions or diffractions), and being recorded by the receiver R. In the case there are two more reflections, that event is called a multiple M. First order free-surface multiples M1 may relate to two interactions in the subsurface and one reflection at the free-surface (water-air interface). With higher order multiples, these interactions may take place more than once. Internal multiples may include at least two upward interactions and one downward interaction in the subsurface. Raw seismic data may comprise a mixture of primaries, surface related multiples, and internal multiples.
At step 703, down-going wavefields and up-going, i.e first and second recorded datasets, are extrapolated from the locus of the receiver depth to a new datum different from the locus of the receiver depth. Various extrapolation methods may be used.
Two separate extrapolations may be performed: one for the up-going wavefield (the second seismic recorded dataset), and one for the down-going wavefield (the first recorded seismic dataset). In one embodiment, the up-going wavefield is inverseextrapolated (i.e., propagated backward in time) to the new datum and the down-going wavefield is forward-extrapolated (i.e., propagated forward in time) to the same new datum. The extrapolated up-going and down-going wavefields appear to have been recorded at this different datum instead of being recorded at a datum determined by the receiver depth.
At step 705, a first image of the surveyed subsurface is generated based on the forwards and backwards propagated seismic dataset generated at step 703.
An “image” of the surveyed subsurface may relate to a single image, as defined by the reflectivity of a plane in the subsurface, but may also relate to an “extended image”, for example, with subsurface offsets for one or more of the x-, y-, or z-directions and/or with non-zero lag arrivals (e.g., as with tau-gathers). In this case, the extended image of the subsurface may be described by: subsurface offsets, example for subsurface offsets in x, y, and z, or tau-gathers: Alternatively, it is possible to have different images for different surface offset ranges. Additionally, the image may include terms to describe a reflection as a function of reflection angle in the subsurface.
An image of the subsurface is generally understood to be a representation of the reflectivity in the earth, defined in space (x-z for 2D, and x-y-z for 3D). This may also be referred to as a migration, the migrated image or the reflectivity. Thus, these terms are used interchangeably herein. An image of the subsurface may be the result of a single-step or optimized migration.
In some embodiments, to form an image, an imaging condition may need to be applied to the data. In one application, the imaging condition is a mathematical function applied to the extrapolated down-going and up-going wavefields to form the image. The most common imaging condition is the cross-correlation imaging condition. Other options are the deconvolution imaging conditions, a variety of which are known in the field, for example, smoothing imaging condition, 2D deconvolution imaging condition or multidimensional deconvolution imaging condition.
The migration, also known as imaging, is a term used to generate an image of the subsurface from sensor recordings following the excitation of a source. Wave-equation migration is commonly applied in the shot domain or in the receiver domain. In the shot domain, a group of traces corresponding to different sensor positions are used for one shot excitation. In the receiver domain, a group of traces are used for one receiver position corresponding to different shot excitations.
There are different types of migration: e.g., one-way migration, two-way reverse time migration (RTM), Kirchhoff migration, common reflection angle migration (CRAM), Wave equation Kirchhoff migration (WEK), etc. One type of migration is an optimized migration. The migration may use an extended image. The migration may be of primaries, multiples, or a combination of primaries and multiples. Recorded data used for the migration may be pre-processed prior to the migration as mentioned herein above.
An optimized migration solves an inverse problem where the goal is to find an image r of the subsurface, which when demigrated, respects the input data D as illustrated in Figure 4. In other words, it is desired to find an image of the subsurface such that when it is used to model or simulate data in the space-time domain, it should equal the recorded data. This may be the result of a least-squares inversion, based on the following definitions:
: Forward extrapolated down-going wavefield;
: Up-going wavefield, reverse extrapolated into the subsurface; and
: Image of the subsurface.
In the following equation, the forward extrapolated wavefield, , is transformed into a reflecting wavefield through multiplication by the image of the subsurface, . The approach is initiated through the injection of the recorded data to form the forward extrapolated wavefield. Ideally, it is desired to find an image of the subsurface which results in a reflecting wavefield that equals the reverse extrapolated up-going wavefield, , for every depth step. An exact solution in practice is generally not possible, hence the following approximation is used:
This may be expressed as a minimization of the difference between both sides of the previous equation, for example, with an L2 norm (other norms may be used, e.g., L1, Cauchy, etc):
Alternatively, the optimization may be expressed as a maximization of the similarity between both sides of the equation, for example, with a zero-lag cross-correlative norm (other norms may be used, e.g., including appropriate normalization of the up- and down-going data, etc.):
It is also possible to pre-multiply both sides of the minimization equation by the adjoint of DF (i.e., the adjoint modelling operator), as shown below. This allows many stable solvers to be used, e.g., conjugate gradients, steepest descent, etc.
The left-hand side of this equation, i.e., may relate to a standard migration with cross-correlation imaging condition. In this case, a misfit may relate to the difference between the left and right hand sides of the equation and the misfit is expressed in the model (image) domain. One alternative is to use the deconvolution imaging condition when applying the adjoint modelling operator. This may not strictly be a least-squares problem that would pass the dot-product test, but in practice it may offer faster convergence. The deconvolution imaging condition would also be used for the adjoint modelling operation on the righthand side:
An alternative formulation may involve evaluating the receiver side (up-going) wavefield, for example:
In this case, the up-going wavefield, u, may be located at the recording datum, and the reflecting wavefield (DFr), for each depth in the image of the subsurface, is forward extrapolated to the surface and accumulated. An alternative way of writing this may be:
where the linear operator F extrapolates the reflecting wavefield from each depth in the image of the subsurface and accumulates at the receiver positions. Constraints may be added to the problem: e.g., image domain sparseness weights or total variation regularization. Data domain confidence or sparseness weights may also be used, for example, to respect a mute function or to respect recorded data positions e.g., following (Poole, 2019).
The optimization problem may be solved with conjugate gradients, steepest descent, an inverse Hessian approach, or another solver. The first iteration of most optimized migrations may relate to a standard migration with an imaging condition. Nonlinear optimized migration algorithms may also be used.
As discussed above, it is possible to image the subsurface using only the primaries, only the multiples, or both the primaries and the multiples. A short summary about primary and multiples imaging is now presented.
Primary imaging using wavefield extrapolation involves deriving an image of the subsurface (also known as reflectivity), r, through the application of an imaging condition between a forward extrapolated source wavefield (down-going), s, and a reverse extrapolated primary wavefield (up-going), p. The source wavefield may relate to the injection of a source wavelet. Figure 5 schematically illustrates the approach in which the image forms where the down-going source wavefield SF and the up-going primary wavefields PR constructively interfere in the subsurface at the location 412. Note that the available information is (1) information about the source, based on which the source wavefield is being determined, and (2) the information recorded at the receiver R. The source wavefield is down propagated to location 412, by a forward propagating operator, after which it is multiplied by the reflectivity r, while the recorded data is backward or reverse propagated from the receiver R to the location 412.
The image of the subsurface may be derived through an optimization scheme, for example, solving for r in the following equations:
where:
SF(t,x,y,z): is a source wavefield forward extrapolated into the subsurface. The use of capital-S relates to the fact the source wavefield will be multiplied (convolved) by the image of the subsurface;
r(x, y, z ) : is the image of the subsurface;
pR(t,x,y,z ) : is the primary wavefield reverse extrapolated into the subsurface; t is the recording time (seconds); and
(x,y,z) refer to the cartesian coordinates.
The optimized migration attempts to find an optimal image of the subsurface, r, to model the known reverse extrapolated primary (up-going) data, pR, which is illustrated in Table 1. This table indicates that the up-going primary arrivals are being generated in primary imaging.
Table 1
The primary wavefield may be raw data, after denoise (e.g., swell noise attenuation, demultiple), or after wavefield separation (e.g., PZ summation or deghosting). It may be provided to the algorithm at the surface datum, or another datum (e.g., ocean bottom survey data, or after re-datuming, for example, to free surface). If an extended image is used, it will be necessary to modify the modelling operator, SF, accordingly. The datasets and operators may alternatively be expressed in the temporal frequency domain.
The primaries may be modelled, i.e., it is possible to estimate primaries in the space-time domain using an image r of the subsurface, as illustrated by the following equation:
In this case the forward extrapolated source wavefield, SF, reflects from the image r of the subsurface to generate a reflecting wavefield, SFr. The reflecting wavefield is subsequently extrapolated forward to the primary, p, recording positions via the forward extrapolation operator (subscript F).
As with primary imaging, the aim of multiple imaging is also to derive an image of the subsurface but this time based on the multiples M and not the primaries P. In the case of multiple imaging, the down-going data utilizes the receivers R as secondary sources, thus significantly increasing the illumination from the surface compared to the primary imaging. Multiple imaging may involve a single iteration with an imaging condition or plural iterations to solve an optimized migration problem, for example, by solving the following equation (which expressed the misfit in the time-space domain) for r
where
DF(t, x,y, z ) : is forward extrapolated down-going data, before demultiple, utilizing all receivers R as secondary sources. Ideally this data should contain primaries and multiples. It may optionally include direct arrival data. The use of capital-D relates to the fact that the down-going wavefield will be multiplied (convolved) by the image of the subsurface; and
: is an up-going wavefield, reverse extrapolated into the subsurface. Ideally this wavefield should contain only free-surface multiples. In practice, this may not be available, in which case an up-going wavefield containing primaries and multiples may be used.
Figure 6 schematically illustrates how an image forms where down-going dF and up-going uR wavefields coincide in the subsurface 410 at the location 412. The rays are shown for a forward extrapolated primary arrival, dF, turning into a first order multiple, uR. Note that each receiver (for example R1) acts as a new source of illumination of the subsurface. This approach also illuminates the subsurface 410 where other multiple orders constructively interfere in the subsurface, as schematically illustrated in Table 2.
Table 2
It is also possible to use a multiple imaging approach using down-going and upgoing datasets derived from dual-sensor towed streamer recordings (Whitmore et al. (2010)). In general, this approach may use any recorded data type, which may optionally have undergone pre-processing. Up-going and down-going wavefields may correspond to data after wavefield separation. It is also possible to use an up-going wavefield, before or after redatum to the free-surface, for up-going and down-going wavefields. Alternatively, it is possible to directly use the recorded data before wavefield separation for up-going and down-going wavefields. Single iteration or optimized migrations may be used. Multiple imaging may result in causal and anti-causal cross-talk.
A third possible imaging approach is to combine the imaging of primaries and multiples. The combination of primary and multiple imaging may involve introducing a source wavefield (or direct-arrival) to the down-going wavefield of the multiple imaging equations, as follows:
In this case, both the primary and multiple arrivals will be used in the imaging, as illustrated in Table 3. In the above equation:
is the reverse extrapolated up-going primary arrivals;
is the reverse extrapolated up-going multiples; In general, + may relate to the recorded data, including both primaries and multiples;
is the forward extrapolated source-wavefield; and
is the forward extrapolated down-going wavefield including primaries and multiples.
Table 3
This may improve the signal-to-noise ratio of the resulting image of the subsurface and reduce cross-talk (especially anti-causal cross-talk, as any primaries in the up-going wavefield will now be modelled from the source). The image of the subsurface may be formed directly using an imaging condition, but more commonly an optimized migration scheme may be used. In general, an optimized migration may be designed to minimize a misfit between synthetic (modelled) data and recorded data. The misfit may relate to a difference, a time-lage cost function, optimal transport cost function, or another cost function. The misfit may be expressed in a time-space domain or in a model domain.
Figure 8A shows an example of a first image I1 of the subsurface obtained from 805. In the example of figure 8A, only the waterbottom was imaged.
In step 707, the first image I1 of the subsurface is used to model primaries and/or model multiples, e.g., Born modelling. This process may also be known as demigration. This process essentially generates the primaries and/or multiples that describe the first image. Modelling may be performed with Kirchhoff (e.g., diffraction modelling), one-way propagation, two-way propagation, etc.
In some embodiments, step 707 comprises generating a first multiple dataset D1 using the down-going wavefield (the first recorded seismic data) and the first image of the subsurface obtained in step 705, as schematically illustrated in table 4. The first multiple dataset D1 comprises first and higher order multiples.
Table 4
This first multiple dataset D1 is then used to generate a second multiple dataset D2. The second multiple dataset D2 is generated using the first multiple dataset D1 and the first image I1 of the surveyed subsurface, as schematically illustrated in Table 5.
Table 5
As illustrated, the second multiple dataset D2 comprises second and higher multiples.
The method according to the invention comprises a step 709, wherein a modified second seismic dataset or partial demultiple dataset is generated by subtracting the first and second multiple data D1 and D2 from the up-going seismic dataset (the second recorded seismic dataset). Such subtraction aims to attenuate the multiples in the input data. Either a straight subtraction or an adaptive subtraction may be used.
In step 711, the modified second seismic data are extrapolated from the locus of the receiver depth to the new datum to generate a modified backwards propagated dataset. As mentioned previously, various extrapolation methods may be used. The modified second seismic data is inverse-extrapolated, i.e., propagated backward in time to the new datum.
In step 713, a second image I2 is generated based on the forwards propagated dataset generated in step 703 and the modified backwards propagated dataset obtained in step 711. The second image I2 may be generated in a similar way as explained herein above. An example of image I2 is shown in Figure 8B. The second image has reduced crosstalk levels compared to an image generated to the same output depth using the first and second datasets (Figure 3B).
In step 715, a third image IF is generated using a shallow depth range from the first image of the subsurface with a deep depth range of the second image of the subsurface, as illustrated in Figure 8C. In this example, the image IF combines the image at the waterbottom from Figure 8A with the image below the waterbottom from Figure 8B. In this case, the shallow image (Figure 8A) from 0 – 100 m depth is combined with the deep image (Figure 8B) from 100 m to 400 m. The contribution from the second image may begin from the maximum depth of the first image. The maximum depth of the second image may be more than 1.5 times the maximum depth of the first image.
In comparison with Figures 3A-C, the image IF has a reduced cross-talk level. According to another embodiment, a method for full-waveform inversion (FWI) with reduced multiple contamination is now discussed with regards to Figure 9.
High-frequency FWI algorithms are designed to estimate subsurface model parameters (e.g. Vp, Vs, density, absorption, anisotropy) that model seismic data. The seismic data may contain primaries, free-surface ghosts and surface-related multiples. Typically, two-way wave-equation modelling engines including a reflecting free-surface may be used.
Multiple contamination can occur in FWI and the resulting FWI Image either when:
1- the subsurface model parameters do not represent the multiple generators accurately enough, or
2- the source-side propagation correlates with the receiver-side propagation at incorrect multiple orders. The second scenario bears some similarities to multiple imaging cross-talk contamination discussed herein above in Figure 7.
We propose to remove multiples from recorded data corresponding to a main multiple generator and input to FWI, where the FWI may only update subsurface model parameters in a depth interval away from the corresponding multiple generator. In one embodiment, the main multiple generator may be a waterbottom reflection.
In step 901, recorded data are provided.
The recorded data may relate to raw recordings such as hydrophone, Vz, Vy, or Vx measurements. Alternatively, the recorded data may be processed data such as an upgoing wavefield or a downgoing data. The recorded data may relate to a plurality of shots recorded by a common-receiver, or a plurality of receivers following excitation of a shotpoint.
In step 902, partial demultiple data is derived. This step may involve multiple modelling followed by subtraction of the modelled multiples. Multiples corresponding to a targeted generator, such as the waterbottom, will be removed and multiples corresponding to non-targeted multiple generators will be left in the data. Primaries corresponding to the targeted generator may optionally also be removed to form the partial demultiple data.
In step 903, the FWI algorithm is applied. The FWI algorithm will model using nonreflective (i.e. smooth) subsurface model parameters in the region corresponding to the targeted generator.
The FWI algorithm may be used to update the subsurface model parameters resulting in modelling the partial demultiple data based on a cost function (e.g. leastsquares (Virieux and Operto, 2009) or time-lag (Zhang et al., 2018)). The term misfit may be considered to be the difference (residual) between the modelled (synthetic) dataset and the recorded dataset, evaluation of a time-lag cost function, an optimal transport cost function, or another cost function.
The input to the FWI modelling may be based on a source signature or an areal source corresponding to the injection of recorded data (before partial demultiple). The FWI modelling may or may not involve a free-surface reflection. Table 6 illustrates the modelling of different arrivals in the data for the areal source case:
Table 6
As illustrated in Table 6, the up-going data used to constrain the FWI algorithm may consist of primaries and/or multiples.
The FWI method according to the invention comprises modifying the recorded data to reduce the level of resulting multiple cross-talk in the subsurface model parameters. The waterbottom may often be considered as the main multiple cross-talk generator. However, the method may be adapted for any other multiple generator.
The data may be modified by attenuating the:
- primary reflection, and/or
- multiples.
The primary reflection may be attenuated by muting, or subtraction of modelled primary arrivals. The multiples may be attenuated using a multiple prediction and multiple attenuation method (e.g. SRME (Berkhout and Verschuur, 1997), partial SRME (Hugonnet, 2002) SRMM (Pica et al., 2005), MWD (Wang et al., 2011), wave-equation targeted multiple modelling (Wiggins, 1988), deconvolution (Biersteker, 2001), wave-equation deconvolution (Poole, 2019)).
The effects corresponding to an interaction with the multiple generator may have been removed from the recorded data to generate the partial demultiple data.
In step 905, an image of the surveyed subsurface is generated. The FWI image using the partial demultiple data has reduced multiple-contamination compared to known methods.
This image is essential for many reasons to the oil and gas industry. One reason is the selection of a site for drilling a well. An accurate image of the subsurface provides clues to the location of the resource (e.g., oil and/or gas) and thus, suggests to those skilled in the art where to drill the well. Another reason is related to the monitoring of an existing production well. It is known that the production rates from wells slowdown in time and thus, it is necessary to inject various fluids to rejuvenate the well production. The images discussed above may be used to monitor the production well, and these images will offer clues about the timing for injecting the above noted fluids, and also in which well to inject them. In addition to oil and gas application, the image could be used to be indicative of geophysical features associated with a natural resource, carbon capture monitoring or wind turbine placement.
The above-discussed methods may be implemented in a computing device illustrated in Figure 10. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein.
Exemplary computing device 1100 suitable for performing the activities described in the exemplary embodiments may include a server 1101. Such a server 1101 may include a central processor (CPU) 1102 coupled to a random access memory (RAM) 1104 and to a read-only memory (ROM) 1106. ROM 1106 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 1102 may communicate with other internal and external components through input/output (I/O) circuitry 1108 and bussing 1110 to provide control signals and the like. Processor 1102 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.
Server 1101 may also include one or more data storage devices, including hard drives 1112, CD-ROM drives 1114 and other hardware capable of reading and/or storing information, such as DVD, etc. In one embodiment, software for carrying out the abovediscussed steps may be stored and distributed on a CD- ROM or DVD 1116, a USB storage device 1118 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as CD-ROM drive 1114, disk drive 1112, etc. Server 1101 may be coupled to a display 1120, which may be any type of known display or presentation screen, such as LCD, plasma display, cathode ray tube (CRT), etc. A user input interface 1122 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.
Server 1101 may be coupled to other devices, such as sources, detectors, etc. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet 1128, which allows ultimate connection to various landline and/or mobile computing devices.
Computing device or computing apparatus 1100 can be configured to implement any of the above-discussed procedures and methods, including combinations thereof.
The disclosed exemplary embodiments provide a computing device, software instructions and a method for seismic data processing. 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 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.
Reference
● Berkhout, A.J and Verschuur, D.J. [1994] Multiple technology: Part 2, migration of multiple reflections. 64<th >Annual International Meeting, SEG, Expanded Abstracts, 1497-1500.
● Berkhout, A. J. and Verschuur, D. J. [1997] Estimation of multiple scattering by iterative inversion, Part I: theoretical consideration. Geophysics, 62, 1586-1595. ● Biersteker, J. [2001] MAGIC: Shell’s surface multiple attenuation technique. 71<st >Annual International Meeting, SEG, Expanded Abstracts.
●Biondi, B.L. [2006] 3D seismic imaging. Investigations in Geophysics, SEG publication. ● Davies, D. and Tillotson, P. [2019] Near field hydrophones – Imaging the subsurface.
81<st >EAGE Conference and Exhibition, Extended Abstracts, We R49.
● Dhelie, P.E., Danielsen, V., Lie, J.E., Støen, S., Dustira, A., Ackers, M.A., Schjelderup, S., Dancer, K., Ramsay, M. and Larssen, B. [2020] 3D seismic in the Barents Sea using quadro point sources for zero-offsets and accurate AVO analysis. 82<nd >EAGE Conference and Exhibition, Extended Abstracts, We Dome34.
● Guitton, A. [2004] Amplitude and kinematic corrections of migrated images for nonunitary imaging operators. Geophysics, 69 (4), 1017–1024.
● Guitton, A., Valenciano, A., Becv, D. and Claerbout, J. [2007] Smoothing imaging condition for shot-profile migration. Geophysics V 77, S149-S154.
● Hugonnet, P. [2002] Partial surface related multiple elimination. 72<nd >Annual International Meeting, SEG, Expanded Abstracts.
● Johnston, D. [2013] Practical applications of time-lapse seismic data, Distinguised Instructor Short Course No.16, SEG publication.
● Lu, S., Whitmore, D., Valenciano, A., Chemingui, N. and Ronholt, G. [2016a] A practical crosstalk attenuation method for separated wavefield imaging. 86<th >Annual International Meeting, SEG, Expanded Abstracts, 4235-4239.
● Lu, S., Liu, F., Chemingui, N., Valenciano, A. and Long, A. [2018] Least-squares fullwavefield migration. The Leading Edge, 37(1), 46-51.
● Mujis, R., Holliger, K., and Robertsson, J. [2005] Prestack depth migration of primary and surface-related multiple reflections. 75<th >Annual International Meeting, SEG, Expanded Abstracts, 2107-2111.
● Mujis, R., Robertsson, J., and Holliger, K. [2007] Prestack depth migration of primary and surface-related multiple reflections: Part II – Identification and removal of residual multiples. Geophysics, Vol 72, No 2, S71-S76.
● Pica, A., Poulain, G., David, B., Magesan, M., Baldock, S., Weisser, T., Hugonnet, P.
and Herrmann, P. [2005] 3D surface-related multiple modeling, principles and results. 75<th >Annual International Meeting, SEG, Expanded Abstracts, 2080-2084.
● Poole, G. [2019] Shallow water surface related multiple attenuation using multisailline 3D deconvolution imaging.81<st >EAGE Conference and Exhibition, Extended Abstracts, Tu R15.
● Poole, T. Curtis, A., Robertsson, J, and Manen D-J. [2010] Deconvolution imaging conditions and cross-talk suppression, Geophysics, 75, W1-W12.
● Spitz, S. [1991] Seismic trace interpolation in the F-X domain. Geophysics, 56, 785-794.
● Valenciano, A., and Biondi, B. [2003] 2D deconvolution imaging condition for shot profile migration. 73<rd >Annual International Meeting, SEG, Expanded Abstracts, 1059–1062.
● Valenciano, A., Biondi, B., and Guitton, A. [2002] Multidimensional imaging condition for shot profile migration. SEP report 111, June 10, 2002, pages 71-81.
● Vinje, V., Lie, J.E., Danielsen, V., Dhelie, P.E., Siliqi, R., Nilsen, C-I., Hicks, E. and Camerer, A. [2017] Shooting over the streamer spread. First Break, 35, 97-104.
●Verschuur, D.J. [2006] Seismic multiple removal techniques: past, present and future.
EAGE publications.
●Virieux, J. and Operto, S. [2009] An overview of full-waveform inversion in exploration geophysics. Geophysics, 74, WCC1–WCC26.
●Wang, P., Jin, H., Xu, S. and Zhang, Y. [2011] Model-based water-layer demultiple.81<st >Annual International Meeting, SEG, Expanded Abstracts, 3551-3555.
●Wang, Y., Zheng, Y., Chang, X. and Yao, Z. [2014] Full Waveform Inversion Using Free-Surface Related Multiples as Natural Blended Sources. 76th EAGE Conference and Exhibition, Extended Abstracts, We E10616.
● Wang, P., Gomes, A., Zhang, Z. and Wang, M. [2016] Least-squares RTM: Reality and possibilities for subsalt imaging.86<th >Annual International Meeting, SEG, Expanded Abstracts, 4204-4209.
●Wiggins, W. [1988] Attenuation of complex water-bottom multiples by waveequation-based prediction and subtraction. Geophysics, 53, 1527-1539.
Whitmore, N.D., Valenciano, A.A., Sollner, W. and Lu, S. [2010] Imaging of primaries and multiples using a dual-sensor towed streamer. 80<th >Annual International Meeting, SEG, Expanded Abstracts, 3187-3192.
Xu, S., Zhang, Y., Lambare, G. [2010] Antileakage Fourier transform for seismic data regularization in higher dimensions. Geophysics, Vol 25, No.6, WB113-WB120. Zhang, D. and Schuster, G.T. [2014] Least-squares reverse time migration of multiples. Geophysics, 79, S11–S21.
Zhang, Z., Mei, J., Lin, F., Huang, R., and Wang, P. [2018] Correcting for salt misinterpretation with full-waveform inversion.88<th >Annual International Meeting, SEG, Expanded Abstracts, 1143-1147.