Search

Seismic datasets contain valuable information that originate from areas of interest in the subsurface; such seismic reflections are however inevitably contaminated by other events created by waves reverberating in the overburden. Multi-Dimensional Deconvolution (MDD) is a powerful technique used at various stages of the seismic processing sequence to create ideal datasets deprived of such overburden effects. Whilst the underlying forward problem is well defined for a single source, a successful inversion of the MDD equations requires availability of a large number of sources alongside prior information introduced in the form of physical preconditioners (e.g., reciprocity). In this work, we reinterpret the cost function of time-domain MDD as a finite-sum functional, and solve the associated inverse problem by means of stochastic gradient descent algorithms, where gradients are computed using a small subset of randomly selected sources. Through synthetic and field data examples, the proposed method is shown to converge more stably than the conventional approach based on full gradients. Stochastic MDD represents a novel, efficient, and robust strategy to deconvolve seismic wavefields in a multi-dimensional fashion.

Stein Variational Gradient Descent (SVGD) is a powerful algorithm for uncertainty quantification algorithm introduced by the statistics community that has recently found applications in geophysical inverse problems. It is a non-parametric, iterative method for making probabilistic predictions by sampling from a set of particles, which are updated using gradient information of the Kullback-Leibler divergence between the target posterior distribution and a user-defined proposal distribution. One of the main benefits of SVGD is its ability to handle complex, high-dimensional distributions. This is particularly useful in geophysics, where datasets can be large, and the subsurface model of interest is high-dimensional. However, the computational cost of SVGD increases with the number of particles used to characterize the posterior. Its performance also depends on the choice of the prior distribution, as it influences the SVGD ability to provide geophysical and geological realism in the posterior samples. There have been efforts to improve the efficiency of SVGD for large data sets, such as by introducing mini-batch techniques and deep learning-based prior aimed at producing high-resolution posterior samples. Along this direction, we are advancing the frontier of the SVGD algorithm by coupling it with the Plug-and-Play (PnP) framework, which allows sampling from a regularized target posterior distribution, where the target posterior distribution is regularized by a convolutional neural network (CNN) based denoiser. We showcase its ability to produce high-resolution, geologically trustworthy samples representative of the subsurface structures on a variety of geophysical problems for reservoir characterization and velocity model building (e.g., full waveform inversion).

Numerical integral operators of the convolution type form the basis of most wave-equation-based methods for the processing and imaging of seismic data. Because several of these methods require the solution of an inverse problem, multiple forward and adjoint passes of the modeling operator are generally required to converge to a satisfactory solution. We have highlighted the memory requirements and computational challenges that arise when implementing such operators on 3D seismic data sets and their use for solving large systems of integral equations. A Python framework is presented that leverages libraries for distributed storage and computing, and it provides a high-level symbolic representation of linear operators. A driving goal for our work is not only to offer a widely deployable, ready-to-use high-performance computing framework, but also to demonstrate that it enables addressing research questions that are otherwise difficult to tackle. To this end, the first example of 3D full-wavefield target-oriented imaging, which comprises two subsequent steps of seismic redatuming, is presented. The redatumed fields are estimated by means of gradient-based inversion using the full data set as well as spatially decimated versions of the data set as a way to investigate the robustness of both inverse problems to spatial aliasing in the input data set. Our numerical example shows that when one spatial direction is finely sampled, satisfactory redatuming and imaging can also be accomplished when the sampling in the other direction is coarser than a quarter of the dominant wavelength. Although aliasing introduces noise in the redatumed fields, they are less sensitive to well-known spurious artifacts compared to cheaper, adjoint-based redatuming techniques. These observations are shown to hold for a relatively simple geologic structure, and although further testing is needed for more complex scenarios, we expect them to be generally valid while possibly breaking down for extreme cases (e.g., highly heterogeneous/scattering media).

Marchenko focusing and imaging are novel methods for correctly handling multiple scattered energy while processing seismic data. However, strict requirements in the acquisition geometry, specifically the colocation of sources and receivers as well as dense and regular sampling, currently constrain their practical applicability. We have reformulated the Marchenko equations to handle the case in which there are gaps in the source geometry while the receiver sampling remains regular (or the opposite, by means of reciprocity). Using synthetic data based on a velocity model that produces strong interbed multiples, we test different solvers for the newly formulated inversion problem, and we compare these results to results obtained by applying standard Marchenko inversion to a previously reconstructed data set. When using the unreconstructed data set, the ability of the Marchenko equations to retrace multiple reflected energy deteriorates. Sparsity-promoting Marchenko inversion, while improving the appearance of focusing functions, barely decreases multiple leakage in gathers and does not visibly improve the final image when compared to standard least-squares inversion. On the other hand, reconstructing the wavefield in advance restores the proper functioning of the Marchenko methods. Further, we test a joint inversion technique designed for time-lapse data with nonrepeated geometries and originally intended to be solved using sparsity-promoting inversion. Motivated by our previous results, we compare images produced by this method to images produced by solving the same joint inversion problem without sparsity constraint. We find that the joint inversion alone hardly improves the resulting images but, when combined with the sparsity constraint, it leads to better noise and multiple suppression and produces a clean time-lapse image. Overall, none of the results from sparsity-promoting inversion techniques match the results obtained when reconstructing the wavefield in advance. We show that this can be explained by the comparatively slow convergence rate of the sparsity-promoting Marchenko inversion.

SUMMARY Full waveform inversion is a high-resolution subsurface imaging technique, in which full seismic waveforms are used to infer subsurface physical properties. We present a novel, target-enclosing, full-waveform inversion framework based on an interferometric objective function. This objective function exploits the equivalence between the convolution and correlation representation formulas, using data from a closed boundary around the target area of interest. Because such equivalence is violated when the knowledge of the enclosed medium is incorrect, we propose to minimize the mismatch between the wavefields independently reconstructed by the two representation formulas. The proposed method requires only kinematic knowledge of the subsurface model, specifically the overburden for redatuming and does not require prior knowledge of the model below the target area. In this sense it is truly local: sensitive only to the medium parameters within the chosen target, with no assumptions about the medium or scattering regime outside the target. We present the theoretical framework and derive the gradient of the new objective function via the adjoint-state method and apply it to a synthetic example with exactly redatumed wavefields. A comparison with FWI of surface data and target-oriented FWI based on the convolution representation theorem only shows the superiority of our method both in terms of the quality of target recovery and reduction in computational cost.

A model is presented to describe how low-frequency vibrations can induce changes to relative permeabilities in dual-porosity dual-permeability rocks. We show that a combination of two physical processes may cause this effect: 1) the wave-induced two-phase fluid flow between soft pores (fractures or cracks) and stiff pores (matrix); and 2) the contact angle hysteresis effect. These two effects lead to redistribution of fluid saturation between fractures and matrix, which may explain alteration of relative permeabilities by small transient strains (∼10−6) in rocks where the fluid flow occurs primarily through the fracture network, and fluid storage occurs predominantly in the porous matrix. We assume that the frequency is low enough that viscous and inertial effects are negligible and that any stress-induced increments of pore pressure are uniform within wetting and nonwetting fluid phases. We demonstrate that pulse-like vibrations are much more efficient in changing relative permeabilities than sinusoidal-shaped vibrations. Furthermore, we show that two waveforms with the same amplitude and frequency but with different polarities could have opposite effects on the alteration of relative permeabilities. Specifically, extensional pulses increase oil relative permeability and decrease water relative permeability in the water-wet rock. Contrary, compressional pulses decrease oil relative permeability and increase water relative permeability in the water-wet rock. To validate theoretical results, we developed a unique low-frequency laboratory setup and experimental methods. We performed drainage and imbibition cycles on a sandstone sample using oil and water as pore fluids. During flooding, the sample was excited by small-strain seismic pulses of a specific shape, with varying polarization direction, frequency, and amplitude. The preliminary laboratory results validate the theoretical model qualitatively. The proposed model may contribute to environmentally friendly and low-cost methods for stimulating hydrocarbon production in conventional reservoirs without invasive chemicals or fracking. In unconventional reservoirs, however, this technology can be envisioned as supplementary to standard EOR/IOR techniques. Normally, the seismic wave energy is too low to increase absolute permeability or decrease oil viscosity; only relative permeabilities can be affected.

Multidimensional deconvolution constitutes an essential operation in a variety of geophysical scenarios at different scales ranging from reservoir to crustal, as it appears in applications such as surface multiple elimination, target-oriented redatuming, and interferometric body-wave retrieval just to name a few. Depending on the use case, active, microseismic, or teleseismic signals are used to reconstruct the broadband response that would have been recorded between two observation points as if one were a virtual source. Reconstructing such a response relies on the the solution of an ill-conditioned linear inverse problem sensitive to noise and artifacts due to incomplete acquisition, limited sources, and band-limited data. Typically, this inversion is performed in the Fourier domain where the inverse problem is solved per frequency via direct or iterative solvers. While this inversion is in theory meant to remove spurious events from cross-correlation gathers and to correct amplitudes, difficulties arise in the estimation of optimal regularization parameters, which are worsened by the fact they must be estimated at each frequency independently. Here we show the benefits of formulating the problem in the time domain and introduce a number of physical constraints that naturally drive the inversion towards a reduced set of stable, meaningful solutions. By exploiting reciprocity, time causality, and frequency-wavenumber locality a set of preconditioners are included at minimal additional cost as a way to alleviate the dependency on an optimal damping parameter to stabilize the inversion. With an interferometric redatuming example, we demonstrate how our time domain implementation successfully reconstructs the overburden-free reflection response beneath a complex salt body from noise-contaminated up- and down-going transmission responses at the target level.

The dynamic time warping (DTW) distance has been used as a misfit function for wave-equation inversion to mitigate the local minima issue. However, the original DTW distance is not smooth; therefore it can yield a strong discontinuity in the adjoint source. Such a weakness does not help nonlinear inverse problems converge to a plausible minimum by any means. We therefore introduce for the first time in geophysics the smooth DTW distance, which has demonstrated its performance in time series classification, clustering, and prediction as the loss function. The fundamental idea of such a distance is to replace the $\min$ operator with its smooth relaxation. Then it becomes possible to define the analytic derivative of DTW distance. The new misfit function is entitled to the differentiable DTW distance. Moreover, considering that the warping path is an indicator of the traveltime difference between the observed and synthetic trace, a penalization term is constructed based on the warping path such that the misfit accumulation by the penalized differentiable DTW distance is weighted in favor of the traveltime difference. Numerical examples demonstrate the advantage of the penalized differentiable DTW misfit function over the conventional non-differentiable one.

The authors would like to thank Chevron for making the synthetic data available. For computer time, this research used the resources of the Supercomputing Laboratory at King Abdullah University of Science & Technology (KAUST) in Thuwal, Saudi Arabia.

The authors thank A. Krull, T.-O. Buchholz, and F. Jug foropen-sourcing their TensorFlow implementation of Noise2Void.For computer time, this research used the resources of the Su-percomputing Laboratory at King Abdullah University of Sci-ence Technology (KAUST) in Thuwal, Saudi Arabia.

With the aim of imaging subsurface discontinuities, seismic data recorded at the surface of the Earth must be numerically re-positioned inside the subsurface where reflections have originated, a process referred to as redatuming. The recently developed Marchenko method is able to handle full-wavefield data including multiple arrivals. A downside of this approach is that a multi-dimensional convolution operator must be repeatedly evaluated to solve an expensive inverse problem. As such an operator applies multiple dense matrix-vector multiplications (MVM), we identify and leverage the data sparsity structure for each frequency matrix and propose to accelerate the MVM step using tile low-rank (TLR) matrix approximations. We study the TLR impact on time-to-solution for the MVM using different accuracy thresholds whilst at the same time assessing the quality of the resulting subsurface seismic wavefields and show that TLR leads to a minimal degradation in terms of signal-to-noise ratio on a 3D synthetic dataset. We mitigate the load imbalance overhead and provide performance evaluation on two distributed-memory systems. Our MPI+OpenMP TLR-MVM implementation reaches up to 3X performance speedup against the dense MVM counterpart from NEC scientific library on 128 NEC SX-Aurora TSUBASA cards. Thanks to the second generation of high bandwidth memory technology, it further attains up to 67X performance speedup compared to the dense MVM from Intel MKL when running on 128 dual-socket 20-core Intel Cascade Lake nodes with DDR4 memory. This corresponds to 110 TB/s of aggregated sustained bandwidth for our TLR-MVM implementation, without suffering deterioration in the quality of the reconstructed seismic wavefields.

In physical acoustic laboratories, wave propagation experiments often suffer from unwanted reflections at the boundaries of the experimental setup. We propose using multidimensional deconvolution (MDD) to post-process recorded experimental data such that the scattering imprint related to the domain boundary is completely removed and only the Green's functions associated with a scattering object of interest are obtained. The application of the MDD method requires in/out wavefield separation of data recorded along a closed surface surrounding the object of interest, and we propose a decomposition method to separate such data for arbitrary curved surfaces. The MDD results consist of the Green's functions between any pair of points on the closed recording surface, fully sampling the scattered field. We apply the MDD algorithm to post-process laboratory data acquired in a two-dimensional acoustic waveguide to characterize the wavefield scattering related to a rigid steel block while removing the scattering imprint of the domain boundary. The experimental results are validated with synthetic simulations, corroborating that MDD is an effective and general method to obtain the experimentally desired Green's functions for arbitrary inhomogeneous scatterers., The Journal of the Acoustical Society of America, 149 (3), ISSN:0001-4966, ISSN:1520-8524

As a result of the world-wide interest in carbon storage and geothermal energy production, increased emphasis is nowadays placed on the development of reliable microseismic monitoring techniques for hazard monitoring related to fluid movement and reactivation of faults. In the process of developing and benchmarking these techniques, the incorporation of realistic noise into synthetic datasets is of vital importance to predict their effectiveness once deployed in the real world. Similarly, the recent widespread use of Machine Learning in seismological applications calls for the creation of synthetic seismic datasets that are indistinguishable from the field data to which they will be applied. Noise generation procedures can be split into two categories: model-based and data-driven. The distributed surface sources approach is the most common method in the first category: however, it is well-known that this fails to capture the complexity of recorded noise (Dean et al., 2015). Pearce and Barley (1977)’s convolutional approach offers a data-driven procedure that has the ability to accurately capture the frequency content of noise however imposes that noise must be stationary. Birnie et al. (2016)’s covariance-based approach removes the stationarity requirement accurately capturing spatio-temporal characterisations of noise, however, like all other data-driven approaches it is constrained to the survey geometry in which the noise data has been collected. In this work, we propose an extension of the covariance-based noise modelling workflow that aims to generate a noise model over a user-defined geometry. The extended workflow comprises of two steps: the first step is responsible for the characterisation of the recorded noise field and the generation of multiple realisations with the same statistical properties, constrained to the original acquisition geometry. Gaussian Process Regression (GPR) is subsequently applied over each time slice of the noise model transforming the model into the desired geometry.The workflow is initially validated on synthetically generated noise with a user-defined input covariance matrix. This allows us to prove that the noise statistics (i.e., covariance and variogram) can be kept almost identical between the noise extracted from the synthetic dataset and the various steps of the noise model procedure. The workflow is further applied to the openly available ToC2ME passive dataset from Alberta, Canada consisting of 69 geophones arranged in a pseudo-random pattern. The noise is modelled and transformed into a 56-sensor, gridded array, which is shown to a very close resemblance to the recorded noise field. Whilst the importance of using realistic noise in synthetic datasets for benchmarking algorithms or training ML solutions cannot be overstated, the ability to transform such noise models into arbitrary receiver geometries opens up a host of new opportunities in the area of survey design. We argue that by coupling the noise generation and monitoring algorithms, the placement of sensors can be optimized based on the expected microseismic signatures as well as the surrounding noise behaviour. This could be of particular interest for geothermal and CO2 storage sites where processing plants are likely to be in close proximity to the permanent monitoring stations.

Inverse problems lie at the core of many geophysical algorithms, from earthquake and exploration seismology, all the way to electromagnetics and gravity potential methods.In 2018, we open-sourced PyLops, a Python-based framework for large-scale inverse problems. By leveraging the concept of matrix-free linear operators – together with the efficiency of numerical libraries such as NumPy, SciPy, and Numba – PyLops solves computationally intensive inverse problems with high-level code that is highly readable and resembles the underlying mathematical formulation. While initially aimed at researchers, its parsimonious software design choices, large test suite, and thorough documentation render PyLops a reliable and scalable software package easy to run both locally and in the cloud.Since its initial release, PyLops has incorporated several advancements in scientific computing leading to the creation of an entire ecosystem of tools: operators can now run on GPUs via CuPy, scale to distributed computing through Dask, and be seamlessly integrated into PyTorch’s autograd to facilitate research in machine-learning-aided inverse problems. Moreover, PyLops contains a large variety of inverse solvers including least-squares, sparsity-promoting algorithms, and proximal solvers highly-suited to convex, possibly nonsmooth problems. PyLops also contains sparsifying transforms (e.g., wavelets, curvelets, seislets) which can be used in conjunction with the solvers. By offering a diverse set of tools for inverse problems under one unified framework, it expedites the use of state-of-the-art optimization methods and compressive sensing techniques in the geoscience domain.Beyond our initial expectations, the framework is currently used to solve problems beyond geoscience, including astrophysics and medical imaging. Likewise, it has inspired the development of the occamypy framework for nonlinear inversion in geophysics. In this talk, we share our experience in building such an ecosystem and offer further insights into the needs and interests of the EGU community to help guide future development as well as achieve wider adoption.

Structural imaging beneath complex overburdens, such as sub-salt or sub-basalt, typically characterized by high-impedance contrasts represents a major challenge for state-of-the-art seismic methods. Reconstructing complex geological structures in the vicinity of and below salt bodies is challenging not only due to uneven, single-sided illumination of the target area but also because of the imperfect removal of surface and internal multiples from the recorded data, as required by traditional migration algorithms. In such tectonic setups, most of the downgoing seismic wavefield is reflected toward the surface when interacting with the overburden's top layer. Similarly, the sub-salt upcoming energy is backscattered at the salt's base. Consequently, the actual energy illuminating the sub-salt reflectors, recorded at the surface, is around the noise level. In diapiric trap systems, conventional seismic extrapolation techniques do not guarantee sufficient quality to reduce exploration and production risks; likewise, seismic-based reservoir characterization and monitoring are also compromised. In this regard, accurate wavefield extrapolation techniques based on the Marchenko method may open up new ways to exploit seismic data.The Marchenko redatuming technique retrieves reliable full-wavefield information in the presence of geologic intrusions, which can be subsequently used to produce artefact-free images by naturally including all orders of multiples present in seismic reflection data. To achieve such a goal, the method relies on the estimation of focusing operators allowing the synthesis of virtual surveys at a given depth level. Still, current Marchenko implementations do not fully incorporate available subsurface models with sharp contrasts, due to the requirements regarding the initialization of the focusing functions. Most importantly, in complex media, even a fairly accurate estimation of a direct wave as a proxy for the required initial focusing functions may not be enough to guarantee sufficiently accurate wavefield reconstruction.In this talk, we will discuss a scattering-based Marchenko redatuming framework which improves the redatuming of seismic surface data in highly complex media when compared to other Marchenko-based schemes. This extended version is designed to accommodate for band-limited, multi-component, and possibly unevenly sampled seismic data, which contain both free-surface and internal multiples, whilst requiring minimum pre-processing steps. The performance of our scattering Marchenko method will be evaluated using a comprehensive set of numerical tests on a complex 2D subsalt model.

Summary Synthetic datasets are vital for the development and benchmarking of new processing and imaging algorithms as well as in the training of machine learning models. It is therefore important that such datasets are generated with realistic noise conditions making them resemble as much as possible their corresponding field datasets. Building on previously developed covariance-based noise modelling, we propose an extension of such an approach that aims to translate a noise model onto a user-defined geometry by means of Gaussian Process Regression. Starting from a synthetic data, we show that noise models can be generated and transformed into a desired geometry whilst keeping the same underlying statistical properties (i.e., covariance and variogram). The modelling procedure is subsequently applied to the ToC2ME passive noise dataset transforming the actual 69-sensor acquisition geometry into a gridded, 56-sensor array. The ability to generate realistic, geometry-independent noise models opens up a host of new opportunities in the area of survey design. We argue that by coupling the noise generation and monitoring algorithms, the placement of sensors could be further optimised based on the expected microseismic signatures as well as the surrounding noise behaviour.

SUMMARY Structural seismic interpretation and quantitative characterization are historically intertwined processes. The latter provides estimates of the properties of the subsurface, which can be used to aid structural interpretation alongside the original seismic data and a number of other seismic attributes. In this work, we redefine this process as an inverse problem which tries to jointly estimate subsurface properties (i.e. acoustic impedance) and a piece-wise segmented representation of the subsurface based on user-defined macroclasses. By inverting for the quantities simultaneously, the inversion is primed with prior knowledge about the regions of interest, whilst at the same time it constrains this belief with the actual seismic measurements. As the proposed functional is separable in the two quantities, these are optimized in an alternating fashion, where each subproblem is solved using the Primal-Dual algorithm. Subsequently, the final segmented model is used as input to an ad hoc workflow that extracts the perimeter of the detected shapes and produces a structural framework (i.e. seismic horizons) consistent with the estimated subsurface properties. The effectiveness of the proposed method is illustrated through numerical examples on synthetic and field data sets.

Can we image and monitor the internal ocean structure at global scales? Can we monitor in vast expanses of the Earth’s cryosphere subsurface with meter-length resolution? Can we characterize the interior structures of asteroids and comets out in space efficiently and with high confidence? At the core of these questions lies the understanding and development of wave-based imaging systems, based on seismic or radar, that rely on highly-sparse, high-quality data, but whose output image quality is comparable to that of densely sampled, wide aperture array-based data. Traditionally, exploration seismology has long relied on wide aperture, dense data sets together with high-end imaging such as reverse-time migration and full-waveform inversion to produce high resolution subsurface models. Given the recent rise of drone-like, autonomous systems, in this talk, we present approaches that can take highly-sparse data as would be recorded by autonomous platforms, into accurate high-resolution images as if they had been acquired by densely-sampled, wide aperture source and receiver arrays. We demonstrate two approaches that could achieve this goal. The first is the use of sparse multicomponent sources and receivers capable of exciting/recording fields and their spatial gradients, together with a gradient-based wavefield reconstruction approach and subsequent imaging. The second approach relies on a new deep learning architecture, the so-called Recurrent Inference Machine, designed specifically for inverse problems – showing that it can surpass the capabilities of deterministic approaches to data reconstruction and imaging. We illustrate these approaches using a numerical model for oceanic turbulence, where we show the compressive sensing potential of these acquisition, reconstruction and imaging methods for acoustic imaging of the ocean's internal structure – overcoming current limitations in data acquisition and processing for seismic oceanography. Finally, we postulate that these approaches, though still in their early days, will pave the way in enabling breakthrough imaging systems at the frontiers of geo-imaging, e.g., for oceanography at global scales, in imaging the Earth’s cryosphere or for planetary exploration.

Summary Numerical integral operators of convolution type lie at the foundations of most wave-equation-based methods for processing and imaging of seismic data. Several of such methods require the solution of an inverse problem, which in turn calls for multiple forward and adjoint passes of the modelling operator. In this abstract we provide some insights into the numerical aspects of solving such systems of integral equations and present a framework that leverages open-source libraries for distributed storage and computing as well as for high-level symbolic representation of linear operators. To validate the effectiveness of our implementation, we evaluate the scalability of the forward and adjoint operations of the well-known time-domain multi-dimensional convolution (MDC) operator with respect to increasing size of the input data and number of computational resources. Finally, we use this implementation to solve the Marchenko equations by means of least-squares inversion for a 3D synthetic dataset composed of up to 9801 sources and receivers.

Seismic methods are used in a wide variety of contexts to investigate subsurface Earth structures, and to explore and monitor resources and waste-storage reservoirs in the upper ~100km of the Earth's subsurface. Reverse-time migration (RTM) is one widely used seismic method which constructs high-frequency images of subsurface structures. Unfortunately, RTM has certain disadvantages shared with other conventional single-scattering-based methods, such as not being able to correctly migrate multiply-scattered arrivals. In principle, the recently-developed Marchenko methods can be used to migrate all orders of multiples correctly. In practice however, using Marchenko methods are costlier to compute than RTM—for a single imaging location, the cost of performing the Marchenko method is several times that of standard RTM, and performing RTM itself requires dedicated use of some of the largest computers in the world for individual datasets. A different imaging strategy is therefore required. We propose a new set of imaging methods which use so-called focusing functions to obtain images with few artifacts from multiply-scattered waves, while greatly reducing the number of points across the image at which the Marchenko method need be applied. Focusing functions are outputs of the Marchenko scheme: they are solutions of wave equations that focus in time and space at particular surface or subsurface locations. However, they are mathematical rather than physical entities, being defined only in reference media that equal to the true Earth above their focusing depths but are homogeneous below. Here, we use these focusing functions as virtual source/receiver surface seismic surveys, the upgoing focusing function being the virtual received wavefield that is created when the downgoing focusing function acts as a spatially distributed source. These source/receiver wavefields are used in three imaging schemes: one allows specific individual reflectors to be selected and imaged. The other two schemes provide either targeted or complete images with distinct advantages over current reverse-time migration methods, such as fewer artifacts and artifacts that occur in different locations. The latter property allows the recently published "combined imaging" method to remove almost all artifacts. We show several examples to demonstrate the methods: acoustic 1D and 2D synthetic examples, and a 2D line from an ocean bottom cable (OBC) field dataset. We discuss an extension to elastic media, which is illustrated by a 1.5D elastic synthetic example.

Seismic reflection data can be redatumed to a specified boundary in the subsurface by solving an inverse (or multidimensional deconvolution) problem. The redatumed data can be interpreted as an extended image of the subsurface at the redatuming boundary, depending on the subsurface offset and time. We retrieve target-enclosed extended images by using two redatuming boundaries, which are selected above and below a specified target volume. As input, we require the upgoing and downgoing wavefields at both redatuming boundaries due to impulsive sources at the earth’s surface. These wavefields can be obtained from actual measurements in the subsurface, they can be numerically modeled, or they can be retrieved by solving a multidimensional Marchenko equation. As output, we retrieved virtual reflection and transmission responses as if sources and receivers were located at the two target-enclosing boundaries. These data contain all orders of reflections inside the target volume but exclude all interactions with the part of the medium outside this volume. The retrieved reflection responses can be used to image the target volume from above or from below. We found that the images from above and below are similar (given that the Marchenko equation is used for wavefield retrieval). If a model with sharp boundaries in the target volume is available, the redatumed data can also be used for two-sided imaging, where the retrieved reflection and transmission responses are exploited. Because multiple reflections are used by this strategy, seismic resolution can be improved significantly. Because target-enclosed extended images are independent on the part of the medium outside the target volume, our methodology is also beneficial to reduce the computational burden of localized inversion, which can now be applied inside the target volume only, without suffering from interactions with other parts of the medium.

Marchenko redatuming is a revolutionary technique to estimate Green’s functions from virtual sources in the subsurface using only data measured at the earth’s surface, without having to place either sources or receivers in the subsurface. This goal is achieved by crafting special wavefields (so-called focusing functions) that can focus energy at a chosen point in the subsurface. Despite its great potential, strict requirements on the reflection response such as knowledge and accurate deconvolution of the source wavelet (including absolute scaling factor) and co-location of sources and receivers have so far challenged the application of Marchenko redatuming to real-world scenarios. I combine the coupled Marchenko equations with a one-way version of the Rayleigh integral representation to obtain a new redatuming scheme that handles internal as well as free-surface multiples using dual-sensor, band-limited seismic data (with an unknown source signature) from any acquisition system that presents arbitrarily located sources above a line of regularly sampled receivers—for example, ocean-bottom, source-over-spread streamer, and horizontal borehole seismic data. The redatuming scheme is validated using synthetic and field data, and the retrieved subsurface wavefields are used for improved structural imaging and taken as input for the computation of true-amplitude angle gathers, which can lead to more accurate amplitude-versus-angle interpretation and velocity analysis.

Resolving details in subsurface reservoir parameters from surface waveform data is a challenging problem, particularly when a reservoir is beneath a complex overburden whose properties are also poorly known. We have developed new metrics for detecting and quantifying errors in subsurface models that are well-suited for target-oriented inversion. Our metric, when combined with state-of-the art redatuming, aims at enabling waveform inversion for target volume parameters with no need to resolve model features elsewhere (e.g., in the overburden). We refer to these metrics as “interferometric objective functions” because they rely on extrapolation from reciprocity integrals commonly used in seismic interferometry. As in seismic interferometry, wavefield extrapolation retrieves the wave response between two points by combining observed data with an extrapolator that describes the response between the subsurface and the data boundary. When the source point is outside a target volume, either forward time or reverse time extrapolation produces the same field. However, because they physically rely on different components from the boundary data, the forward time and reverse time extrapolated fields are only equal when the model used is consistent with the real subsurface within the target volume. As such, we use the difference between the forward time and reverse time extrapolated fields to define subsurface-domain metrics that quantify model errors. Our approach thus provides a new metric for target-oriented nonlinear inversion in the subsurface domain, one that fundamentally differs from other subsurface domain metrics based on, e.g., focusing or image extensions.

Summary Marchenko focusing and imaging are novel methods for processing seismic data while correctly handling multiply scattered energy. Strict requirements in the acquisition geometry, specifically co-location of sources and receivers as well as dense and regular sampling, currently constrain their practical applicability. We reformulate the Marchenko equations to handle the case where there are gaps in the source geometry while receiver sampling remains regular. Comparing different solvers for the newly formulated inversion problem, we find that sparse Marchenko inversion enhances the recovery of focusing functions, filling the gaps caused by the missing sources. This does, however, not translate into a significant improvement in the subsequently produced images, as both least-squares and sparse inversion struggle to eliminate overburden effects in the presence of gaps. Further, we develop a method for co-processing two time-lapse datasets with different source geometries using sparse Marchenko inversion to image 4D effects. Sparse joint Marchenko inversion of multiple datasets results in clear time-lapse images, despite non-repeated source geometries and large fractions of missing sources.

Marchenko-based redatuming methods aim to produce full-waveform propagators from single-sided reflection data and an estimate of the subsurface velocity model. In this talk, I will first review various alternative approaches proposed in the literature and discuss how they differ from each other and from the original formulation of Wapenaar et al (2014). Using a common notation throughout, the key novelty of each method will be highlighted as well as their strengths and limitations with respect to application to field data. Finally, I will discuss the way forward, both in terms of implementation details (e.g., 3D) and theoretical advances (e.g., active role of multiples), that are the topics of current research in the community and can be seen as key enablers for a more widespread use of the Marchenko equation in geophysical applications.

Methods for wavefield injection are commonly used to extrapolate seismic data in reverse time migration (RTM). Injecting a single component of the acoustic field, for example, pressure, leads to ambiguity in the direction of propagation. Each recorded wavefront is propagated both upward and downward, and spurious (or ghost) reflectors are created alongside real reflectors in the subsurface image. Thus, wavefield separation based on the combination of pressure and particle velocity data is generally performed prior to imaging to extract only the upgoing field from multicomponent seabed or towed marine seismic recordings. By instead combining vector-acoustic (VA) data with monopole- and dipole-type propagators in the extrapolation of shot or receiver gathers, we show that wavefield separation (or deghosting) can instead be performed “on-the-fly” at limited additional cost. This strategy was successfully applied to a line of a North Sea ocean-bottom cable data set, acquired over the Volve field. We then evaluate additional advantages over standard RTM with decomposed fields such as improved handling of the directivity information contained in the acquired VA data for clearer shallow sections and better focused space-lag common image gathers, and imaging of the downgoing component without the need for additional finite-difference modeling via mirror migration. Finally, we prove the robustness of our method with respect to sparse and irregular receiver sampling.

S U M M A R Y Seismic interferometry comprises a suite of methods to redatum recorded wavefields to those that would have been recorded if different sources (so-called virtual sources) had been activated. Seismic interferometry by cross-correlation has been formulated using either two-way (for full wavefields) or one-way (for directionally decomposed wavefields) representation theorems. To obtain improved Green’s function estimates, the cross-correlation result can be deconvolved by a quantity that identifies the smearing of the virtual source in space and time, the so-called point-spread function. This type of interferometry, known as interferometry by multidimensional deconvolution (MDD), has so far been applied only to one-way directionally decomposed fields, requiring accurate wavefield decomposition from dual (e.g. pressure and velocity) recordings. Here we propose a form of interferometry by multidimensional deconvolution that uses full wavefields with two-way representations, and simultaneously invert for pressure and (normal) velocity Green’s functions, rather than only velocity responses as for its one-way counterpart. Tests on synthetic data show that two-way MDD improves on results of interferometry by cross-correlation, and generally produces estimates of similar quality to those obtained by one-way MDD, suggesting that the preliminary decomposition into upand downgoing components of the pressure field is not required if pressure and velocity data are jointly used in the deconvolution. We also show that constraints on the directionality of the Green’s functions sought can be added directly into the MDD inversion process to further improve two-way multidimensional deconvolution. Finally, as a by-product of having pressure and particle velocity measurements, we adapt oneand two-way representation theorems to convert any particle velocity receiver into its corresponding virtual dipole/gradient source by means of MDD. Thus data recorded from standard monopolar (e.g. marine) pressure sources can be converted into data from dipolar (derivative) sources at no extra acquisition cost.

Summary Marchenko redatuming is a revolutionary technique to estimate Green’s functions from virtual sources in the subsurface using data measured at the Earth’s surface, without having to place either sources or receivers in the subsurface. Such goal is achieved by crafting special wavefields (so-called focusing functions) that can focus energy at any chosen point in the subsurface. Despite its great potential, strict requirements on the reflection data have so far challenged the application of Marchenko redatuming to real scenarios. In this work, I combine the coupled Marchenko equations with a one-way version of the Rayleigh integral representation to obtain a new, more flexible redatuming scheme that handles internal as well as free-surface multiples using dual-sensor, band-limited seismic data from any acquisition system that presents arbitrarily located sources above a line of regularly sampled receivers. The proposed redatuming scheme is validated on synthetic and field data.

reflection measurements at the Earth's surface and an estimate of the direct wave from a focusing point in the subsurface to the acquisition level in order to retrieve the Green’s function between the two. However, the method suffers when an erroneous scaling of the amplitudes and/or a constant phase-shift is introduced to the reflection response, and artifacts caused by internal multiples are not fully suppressed in the retrieved subsurface wave fields. We propose a workflow to calibrate the reflection response, prior to Green's function retrieval via Marchenko focusing, using additional information in the form of a VSP dataset. First a virtual VSP dataset is estimated via Marchenko focusing, to subsequently compare its upgoing component to the upgoing part of the recorded VSP. Thereby, making it possible to correct for a constant phase shift applied to the reflection data. By identifying the minimum residual energy between the virtual VSP and the recorded VSP wave fields, an erroneous scaling of the reflection response can also be corrected. This workflow leads to a more robust Marchenko focusing, where the reflection response can be redatumed to a target zone in the subsurface and used for ghost-free imaging.

The Marchenko redatuming approach reconstructs wavefields at depth that contain not only primary reflections, but also multiply-scattered waves. While such fields in principle contain additional subsurface information, conventional imaging approaches cannot tap into the information encoded in internal multiples in a trivial manner. We discuss a new approach that uses the full information contained in Marchenko-redatumed fields, whose output are local reflection and transmission responses that fully enclose a target volume at depth, without contributions from over- or under-burden structures. To obtain the Target-Enclosing Extended Images (TEEIs) we solve a multi-dimensional deconvolution (MDD) problem that can be severely ill-posed, so we offer stable estimates to the MDD problem that rely on the physics of the Marchenko scheme. We validate our method on ocean-bottom field data from the North Sea. In our field data example, we show that the TEEIs can be used for reservoir-targeted imaging using reflection and, for the first time, local transmission responses, shown to be the direct by-product of using internal multiples in the redatuming scheme. Finally, we present local, TEEI-derived reflection and transmission images of the target volume at depth that are structurally consistent with a benchmark image from conventional migration of surface data.

A central component of imaging methods is receiver-side wavefield backpropagation or extrapolation in which the wavefield from a physical source scattered at any point in the subsurface is estimated from data recorded by receivers located near or at the Earth’s surface. Elastic reverse-time migration usually accomplishes wavefield extrapolation by simultaneous reversed-time ‘injection’ of the particle displacements (or velocities) recorded at each receiver location into a wavefield modeling code. Here, we formulate an exact integral expression based on reciprocity theory that uses a combination of velocity-stress recordings and quadrupole-dipole backpropagating sources, rather than the commonly used approximate formula involving only particle velocity data and dipole backpropagating sources. The latter approximation results in two types of nonphysical waves in the scattered wavefield estimate: First, each arrival contained in the data is injected upward and downward rather than unidirectionally as in the true time-reversed experiment; second, all injected energy emits compressional and shear propagating modes in the model simulation (e.g., if a recorded P-wave is injected, both P and S propagating waves result). These artifacts vanish if the exact wavefield extrapolation integral is used. Finally, we show that such a formula is suitable for extrapolation of ocean-bottom 4C data: Due to the fluid-solid boundary conditions at the seabed, the data recorded in standard surveys are sufficient to perform backpropagation using the exact equations. Synthetic examples provide numerical evidence of the importance of correcting such errors.