Internal Multiples Research Articles

Marchenko imaging: Imaging with primaries, internal multiples, and free-surface multiples

Recent work on retrieving the Green’s function with the Marchenko equation shows how these functions for a virtual source in the subsurface can be obtained from reflection data. The response to the virtual source is the Green’s function from the location of the virtual source to the surface. The Green’s function is retrieved using only the reflection response of the medium and an estimate of the first arrival at the surface from the virtual source. Current techniques, however, only include primaries and internal multiples. Therefore, all surface-related multiples must be removed from the reflection response prior to Green’s function retrieval. We have extended the Marchenko equation to retrieve the Green’s function that includes primaries, internal multiples, and free-surface multiples. In other words, we have retrieved the Green’s function in the presence of a free surface. The information needed for the retrieval is the same as the current techniques, with the only difference being that the reflection response now also includes free-surface multiples. The inclusion of these multiples makes it possible to include them in the imaging operator, and it obviates the need for surface-related multiple elimination. This type of imaging with Green’s functions is called Marchenko imaging.

Elastic P- and S-wave autofocus imaging with primaries and internal multiples

Conventional seismic imaging methods rely on the single-scattering Born approximation, requiring the removal of multiply scattered events from reflection data prior to imaging. Additionally, many methods use an acoustic approximation, representing the solid earth as an acoustic (fluid) medium. We have developed imaging methods for (solid) elastic media that use primaries and internal multiples, including their PS and SP conversions, thus obviating the need for internal multiple removal and improving handling of internal conversions. Our methods rely on the elastic autofocusing method, which estimates full multicomponent elastodynamic Green’s functions from virtual sources interior to the medium to receivers placed on the surface. They require only surface seismic reflection data and estimates of the direct waves from virtual sources interior to the medium, both of which are commonly available at the imaging step of seismic processing. We tested our methods on a synthetic model with constant P and S velocities and vertical and horizontal density variations, by producing, for the first time to our knowledge, PP and SS images from elastic autofocusing that are compared to reference seismic images based on conventional methods. The effects of multiples are greatly attenuated in the images, with fewer spurious reflectors than are observed when using Born imaging.

An illustration of adaptive Marchenko imaging

In Marchenko imaging, wavefields are retrieved at specified focal points in the subsurface through an iterative scheme derived from the multidimensional Marchenko equation. The method requires seismic-reflection data at the earth's surface (after free-surface multiple elimination) and an estimate of the direct wavefield from the surface to each focal point, which can be computed, for instance, in a macrovelocity model. In the first iteration, the direct wavefield is crosscorrelated with the reflection data. This operation is identical to inverse-wavefield extrapolation as is applied commonly in various imaging schemes, for instance, in reverse time migration (RTM). At each succeeding iteration, the result of the previous iteration is truncated in time and crosscorrelated with the reflection data again. To obtain a seismic image, a multidimensional deconvolution-based imaging condition can be applied to the retrieved wavefields. By this approach, both primary reflections and internal multiples contribute to the construction of the image. Alternatively, a crosscorrelation-based imaging condition can be used in which only the primary reflections are imaged and the contributions of internal multiples are subtracted. The latter strategy offers more flexibility because the subtraction of redatumed internal multiples can be implemented adaptively. Through this approach, the artifacts from internal multiples can be removed effectively from a conventional RTM image.

The possibilities of least-squares migration of internally scattered seismic energy

Approximate images of the earth’s subsurface structures are usually obtained by migrating surface seismic data. Least-squares migration, under the single-scattering assumption, is used as an iterative linearized inversion scheme to suppress migration artifacts, deconvolve the source signature, mitigate the acquisition fingerprint, and enhance the spatial resolution of migrated images. The problem with least-squares migration of primaries, however, is that it may not be able to enhance events that are mainly illuminated by internal multiples, such as vertical and nearly vertical faults or salt flanks. To alleviate this problem, we adopted a linearized inversion framework to migrate internally scattered energy. We apply the least-squares migration of first-order internal multiples to image subsurface vertical fault planes. Tests on synthetic data demonstrated the ability of the proposed method to resolve vertical fault planes, which are poorly illuminated by the least-squares migration of primaries only. The proposed scheme is robust in the presence of white Gaussian observational noise and in the case of imaging the fault planes using inaccurate migration velocities. Our results suggested that the proposed least-squares imaging, under the double-scattering assumption, still retrieved the vertical fault planes when imaging the scattered data despite a slight defocusing of these events due to the presence of noise or velocity errors.

The electromagnetic response in a layered vertical transverse isotropic medium: A new look at an old problem

ABSTRACTWe determined that the electromagnetic vertical transverse isotropic response in a layered earth can be obtained by solving two equivalent scalar equations, which were for the vertical electric field and for the vertical magnetic field, involving only a scalar global reflection coefficient. Besides the complete derivation of the full electromagnetic response, we also developed the corresponding computer code called EMmod, which models the full electromagnetic fields including internal multiples in the frequency-wavenumber domain and obtains the frequency-space domain solutions through a Hankel transformation by computing the Hankel integral using a 61-point Gauss-Kronrod integration routine. The code is able to model the 3D electromagnetic field in a 1D earth for diffusive methods such as controlled source electromagnetics as well as for wave methods such as ground penetrating radar. The user has complete freedom to place the source and the receivers in any layer. The modeling is illustrated with ...

Imaging blended vertical seismic profiling data using full-wavefield migration in the common-receiver domain

For vertical-seismic-profiling (VSP) measurements, the use of blended acquisition, with time-overlapping shot records, can greatly reduce the downtime and, thereby, provide large cost savings. For directly imaging blended VSP measurements, we have used full-wavefield migration (FWM). FWM is an inversion-based imaging scheme that enables us to use any kind of complex source wavefield to estimate the subsurface reflectivity, using all the multiples (surface and internal) in the recorded data. The multiple scattering helps in improving the illumination as well as the vertical resolution of the image. In this scheme, active deblending is not required because the imaging process itself acts as a deblending procedure. We tested the potential of FWM to image blended VSP data, using simple and complex synthetic models. We clearly determined that using the primaries, surface multiples, and internal multiples enhanced the illumination away from the well trajectory, in which blending noise was suppressed due to the inversion scheme. We observed that some blending crosstalk noise leaked into the images with large blending factors. Such noise could be further reduced with additional constraints in the involved least-squares inversion process.

Successive application of the layer-related CFP method for internal multiple removal

Internal multiples are very difficult to remove due to their complex raypaths and poor velocity discrimination with primaries. The layer-related common focus point (CFP) method has proven to be effective for internal multiple removal, however, single application leads to multiple leakage when there are several strong reflecting boundaries in the subsurface. In order to reduce this leakage, we propose successive application of the layer-related CFP method for internal multiple removal through cascaded processing of several time levels. In this paper, we concisely reformulate the theory of the boundary- and layer-related CFP methods and compare their robustness to velocity errors. For the layer-related version in particular, we illustrate the specific steps of the method using synthetic data examples. Finally, the successive layer-related CFP method is tested on Mississippi Canyon field data.

Internal multiple adaptive subtraction using Huber norm

There are two steps in internal multiple elimination, the first step is multiple prediction (predicted multiple model by inverse scattering series (ISS) method from marine seismic data only), the second step is adaptive matching subtraction, subtract multiple from seismic data. Because of the complex generation mechanism of internal multiple, the second step is always the main challenge of internal multiple suppression. A new subtraction method using Huber norm is proposed to remove internal multiples in this paper, it can deal with the negative influence of amplitude and phase in internal multiple subtraction. In this approach, the multi-channel Huber norm function minimized by quasi-Newtonian method that has the potential for more accurate and robust multiple subtraction. Tests with synthetic and field data suggest that internal multiple adaptive subtraction using the Huber norm gives highly satisfactory result, and shows its effectiveness.

Point‐spread functions for interferometric imaging

ABSTRACTInterferometric redatuming is a data‐driven method to transform seismic responses with sources at one level and receivers at a deeper level into virtual reflection data with both sources and receivers at the deeper level. Although this method has traditionally been applied by cross‐correlation, accurate redatuming through a heterogeneous overburden requires solving a multidimensional deconvolution problem. Input data can be obtained either by direct observation (for instance in a horizontal borehole), by modelling or by a novel iterative scheme that is currently being developed. The output of interferometric redatuming can be used for imaging below the redatuming level, resulting in a so‐called interferometric image. Internal multiples from above the redatuming level are eliminated during this process. In the past, we introduced point‐spread functions for interferometric redatuming by cross‐correlation. These point‐spread functions quantify distortions in the redatumed data, caused by internal multiple reflections in the overburden. In this paper, we define point‐spread functions for interferometric imaging to quantify these distortions in the image domain. These point‐spread functions are similar to conventional resolution functions for seismic migration but they contain additional information on the internal multiples in the overburden and they are partly data‐driven. We show how these point‐spread functions can be visualized to diagnose image defocusing and artefacts. Finally, we illustrate how point‐spread functions can also be defined for interferometric imaging with passive noise sources in the subsurface or with simultaneous‐source acquisition at the surface.

Internal multiple prediction and removal using Marchenko autofocusing and seismic interferometry

Standard seismic processing steps such as velocity analysis and reverse time migration (imaging) usually assume that all reflections are primaries: Multiples represent a source of coherent noise and must be suppressed to avoid imaging artifacts. Many suppression methods are relatively ineffective for internal multiples. We show how to predict and remove internal multiples using Marchenko autofocusing and seismic interferometry. We first show how internal multiples can theoretically be reconstructed in convolutional interferometry by combining purely reflected, up- and downgoing Green’s functions from virtual sources in the subsurface. We then generate the relevant up- and downgoing wavefields at virtual sources along discrete subsurface boundaries using autofocusing. Then, we convolve purely scattered components of up- and downgoing Green’s functions to reconstruct only the internal multiple field, which is adaptively subtracted from the measured data. Crucially, this is all possible without detailed modeled information about the earth’s subsurface. The method only requires surface reflection data and estimates of direct (nonreflected) arrivals between subsurface virtual sources and the acquisition surface. The method is demostrated on a stratified synclinal model and shown to be particularly robust against errors in the reference velocity model used.

The electromagnetic response in a layered vertical transverse isotropic medium: A new look at an old problem

We determined that the electromagnetic vertical transverse isotropic response in a layered earth can be obtained by solving two equivalent scalar equations, which were for the vertical electric field and for the vertical magnetic field, involving only a scalar global reflection coefficient. Besides the complete derivation of the full electromagnetic response, we also developed the corresponding computer code called EMmod, which models the full electromagnetic fields including internal multiples in the frequency-wavenumber domain and obtains the frequency-space domain solutions through a Hankel transformation by computing the Hankel integral using a 61-point Gauss-Kronrod integration routine. The code is able to model the 3D electromagnetic field in a 1D earth for diffusive methods such as controlled source electromagnetics as well as for wave methods such as ground penetrating radar. The user has complete freedom to place the source and the receivers in any layer. The modeling is illustrated with three examples, which aim to present the different capabilities of EMmod, while assessing its correctness.

Generalized internal multiple imaging

Internal multiples deteriorate the image when the imaging procedure assumes only single scattering, especially if the velocity model does not have sharp contrasts to reproduce such scattering in the Green’s function through forward modeling. If properly imaged, internal multiples (internally scattered energy) can enhance the seismic image. Conventionally, to image internal multiples, accurate, sharp contrasts in the velocity model are required to construct a Green’s function with all the scattered energy. As an alternative, we have developed a generalized internal multiple imaging procedure that images any order internal scattering using the background Green’s function (from the surface to each image point), constructed from a smooth velocity model, usually used for conventional imaging. For the first-order internal multiples, the approach consisted of three steps, in which we first back propagated the recorded surface seismic data using the background Green’s function, then crosscorrelated the back-propagated data with the recorded data, and finally crosscorrelated the result with the original background Green’s function. This procedure images the contribution of the recorded first-order internal multiples, and it is almost free of the single-scattering recorded energy. The cost includes one additional crosscorrelation over the conventional single-scattering imaging application. We generalized this method to image internal multiples of any order separately. The resulting images can be added to the conventional single-scattering image, obtained, e.g., from Kirchhoff or reverse-time migration, to enhance the image. Application to synthetic data with reflectors illuminated by multiple scattering (double scattering) demonstrated the effectiveness of the approach.

Full‐wavefield migration of vertical seismic profiling data: using all multiples to extend the illumination area

ABSTRACTVertical seismic profiling data can provide a high‐resolution reservoir image because the receivers are located close to the reservoir and, hence, the wavefield encounters less distortion in the overburden compared to surface seismic data. However, conventional migration of vertical seismic profiling data using only upgoing primary wavefields often suffers from poor illumination, imaging artefacts and low image reliability, especially at image locations away from the well. We propose full‐wavefield migration to image vertical seismic profiling data using the primaries and all orders of surface multiples and internal multiples. The downgoing internal or surface multiples are treated as an additional source of illumination. In full‐wavefield migration, we aim to estimate a high‐resolution reflectivity image of the subsurface with this extended illumination. The algorithm is recursive in depth like conventional wavefield extrapolation‐based migration, however, it incorporates the non‐linear transmission and scattering effects at each depth level. Further, full‐wavefield migration is posed as a constrained least‐squares inversion problem that could be solved using a conjugate gradient method. We propose an iterative full‐wavefield forward modelling engine as the backbone of this inversion scheme. The parameter used in the modelling is subsurface reflectivity. Full‐wavefield modelling iteratively incorporates the non‐linearity of the wavefield due to multiple scattering, where every iteration utilizes one higher order of scattered wavefields to estimate the subsurface reflectivity. In addition, the constrained inversion helps in reducing the extrapolation artefacts and provides a high‐resolution image of the reservoir. In this paper, we discuss the concept of full‐wavefield migration and demonstrate its potential as an imaging tool for vertical seismic profiling data using synthetic examples.

Data-driven wavefield focusing and imaging with multidimensional deconvolution: Numerical examples for reflection data with internal multiples

Standard imaging techniques rely on the single scattering assumption. This requires that the recorded data do not include internal multiples, i.e., waves that have bounced multiple times between reflectors before reaching the receivers at the acquisition surface. When multiple reflections are present in the data, standard imaging algorithms incorrectly image them as ghost reflectors. These artifacts can mislead interpreters in locating potential hydrocarbon reservoirs. Recently, we introduced a new approach for retrieving the Green’s function recorded at the acquisition surface due to a virtual source located at depth. We refer to this approach as data-driven wavefield focusing. Additionally, after applying source-receiver reciprocity, this approach allowed us to decompose the Green’s function at a virtual receiver at depth in its downgoing and upgoing components. These wavefields were then used to create a ghost-free image of the medium with either crosscorrelation or multidimensional deconvolution, presenting an advantage over standard prestack migration. We tested the robustness of our approach when an erroneous background velocity model is used to estimate the first-arriving waves, which are a required input for the data-driven wavefield focusing process. We tested the new method with a numerical example based on a modification of the Amoco model.

Quantification of signal and noise in two-way transmitted wavefields to improve well ties in Cooper Basin, Australia

Permian coal beds at 2400–2900 m depth in Cooper Basin, Australia have normal-incident reflection coefficient values as large as [Formula: see text]. If internal multiples are included in synthetic seismograms, excellent correlations exist between the synthetic seismogram and seismic, even when more than 50 coal beds are present. However, neither the synthetic seismogram nor the seismic tie the well-log lithologic boundaries because the incident wavefield that strikes a lithologic boundary and returns to the surface contains a signal wavelet followed by high-amplitude noise, which are interbed multiple reflections. Because the spectra of the signal and noise coda at a given two-way time normally do not overlap, time-varying Gaussian filters applied to the near-offset stack enhance the signal and suppress the noise coda. After filtering, the apparent time delay of reflections introduced by the coal beds is removed with variable time shifts (time compression), based on the estimated time-varying signal wavelets. The two-step process of low-pass filtering and compression yields seismic events that successfully tie the lithologic boundaries in the borehole, although with limited resolution. Our preliminary tests on a seismic line indicate that the horizon event associated with the base of a 500-m-thick coal sequence is more coherently imaged with our processing than with conventional processing.

Marchenko imaging

Traditionally, the Marchenko equation forms a basis for 1D inverse scattering problems. A 3D extension of the Marchenko equation enables the retrieval of the Green’s response to a virtual source in the subsurface from reflection measurements at the earth’s surface. This constitutes an important step beyond seismic interferometry. Whereas seismic interferometry requires a receiver at the position of the virtual source, for the Marchenko scheme it suffices to have sources and receivers at the surface only. The underlying assumptions are that the medium is lossless and that an estimate of the direct arrivals of the Green’s function is available. The Green’s function retrieved with the 3D Marchenko scheme contains accurate internal multiples of the inhomogeneous subsurface. Using source-receiver reciprocity, the retrieved Green’s function can be interpreted as the response to sources at the surface, observed by a virtual receiver in the subsurface. By decomposing the 3D Marchenko equation, the response at the virtual receiver can be decomposed into a downgoing field and an upgoing field. By deconvolving the retrieved upgoing field with the downgoing field, a reflection response is obtained, with virtual sources and virtual receivers in the subsurface. This redatumed reflection response is free of spurious events related to internal multiples in the overburden. The redatumed reflection response forms the basis for obtaining an image of a target zone. An important feature is that spurious reflections in the target zone are suppressed, without the need to resolve first the reflection properties of the overburden.

Autofocus Imaging: Image reconstruction based on inverse scattering theory

Conventional imaging algorithms assume single scattering and therefore cannot image multiply scattered waves correctly. The multiply scattered events in the data are imaged at incorrect locations resulting in spurious subsurface structures and erroneous interpretation. This drawback of current migration/imaging algorithms is especially problematic for regions where illumination is poor (e.g., subsalt), in which the spurious events can mask true structure. Here we discuss an imaging technique that not only images primaries but also internal multiples accurately. Using only surface-reflection data and direct-arrivals, we generate the up- and down-going wavefields at every image point in the subsurface. An imaging condition is applied to these up- and down-going wavefields directly to generate the image. Because the above algorithm is based on inverse-scattering theory, the reconstructed wavefields are accurate and contain multiply scattered energy in addition to the primary event. As corroborated by our synthetic examples, imaging of these multiply scattered energy helps eliminate spurious reflectors in the image. Other advantages of this imaging algorithm over existing imaging algorithms include more accurate amplitudes, target-oriented imaging, and a highly parallelizable algorithm.

Seismic reflector imaging using internal multiples with Marchenko-type equations

We present an imaging method that creates a map of reflection coefficients in correct one-way time with no contamination from internal multiples using purely a filtering approach. The filter is computed from the measured reflection response and does not require a background model. We demonstrate that the filter is a focusing wavefield that focuses inside a layered medium and removes all internal multiples between the surface and the focus depth. The reflection response and the focusing wavefield can then be used for retrieving virtual vertical seismic profile data, thereby redatuming the source to the focus depth. Deconvolving the upgoing by the downgoing vertical seismic profile data redatums the receiver to the focus depth and gives the desired image. We then show that, for oblique angles of incidence in horizontally layered media, the image of the same quality as for 1D waves can be constructed. This step can be followed by a linear operation to determine velocity and density as a function of depth. Numerical simulations show the method can handle finite frequency bandwidth data and the effect of tunneling through thin layers.

Acoustic nonlinear full-waveform inversion on an outcrop-based detailed geological and petrophysical model (Book Cliffs, Utah)

Analog outcrops are commonly used to develop predictive reservoir models and provide quantitative parameters that describe the architecture and facies distribution of sedimentary deposits at a subseismic scale, all of which aids exploration and production strategies. The focus of this study is to create a detailed geological model that contains realistic reservoir parameters and to apply nonlinear acoustic full-waveform prestack seismic inversion to this model to investigate whether this information can be recovered and to examine which geological features can be resolved by this process.Outcrop data from the fluviodeltaic sequence of the Book Cliffs (Utah) are used for the geological and petrophysical two-dimensional model. Eight depositional environments are populated with average petrophysical reservoir properties adopted from a North Sea field. These units are termed lithotypes here. Synthetic acoustic prestack seismic data are then generated with the help of an algorithm that includes all internal multiples and transmission effects. A nonlinear acoustic full-waveform inversion is then applied to the synthetic data, and two media parameters, compressibility (inversely related to the square of the compressional wave velocity vP) and bulk density, ρ, are recovered at a resolution higher than the shortest wavelength in the data. This is possible because the inversion exploits the nonlinear nature of the relationship between the recorded data and the medium contrast properties. In conventional linear inversion, these details remain masked by the noise caused by the nonlinear effects in the data. Random noise added to the data is rejected by the nonlinear inversion, contributing to improved spatial resolution. The results show that the eight lithotypes can be successfully recovered at a subseismic scale and with a low degree of processing artifacts. This technique can provide a useful basis for more accurate reservoir modeling and field development planning, allowing targeting of smaller reservoir units such as distributary channels and lower shoreface sands.

Modelling and Removing WAZ OBC Interbed Multiples

Interbed (or internal) multiples, caused by ringing between strong reflectors, can play havoc with the interpretation of primary reflections deeper within a seismic section, and have been historically difficult to suppress. True wide-azimuth, relatively sparsely acquired 3D data (like OBC) add to the complexity of this problem.For this case study, we present our testing and on-going developments on interbed multiple suppression. Two methods are shown. The first is an interpretive pattern- recognition technique applied after migration. The second, applied to unmigrated data, is a new version of wavefield modelling that attempts to predict interbed multiples without any knowledge of the actual multiple generators. A well-constrained adaptive subtraction methodology is required, and the testing and development of a new 3D subtraction algorithm is reviewed. Both methods have been shown to perform well.