Receiver Wavefields Research Articles

Reverse time migration using water‐bottom‐related multiples

ABSTRACTReverse time migration of multiples can be used to construct subsurface structures where primaries cannot illuminate well. However, the images generated using multiples suffer from severe artefacts due to the cross‐talks created by interference among unrelated multiples. We developed a migration approach using water‐bottom‐related multiples to reduce these cross‐talk artefacts. This approach first isolates primaries from the original data and predicts water‐column primaries. The nth‐order water‐column multiples can be obtained by auto‐convolving the water‐column primaries n times, followed by convolving the nth‐order water‐column multiples with the primaries to extract the (n+1)th‐order water‐bottom‐related multiples. The approach takes the nth‐order water‐column multiples as the secondary source and regards the (n+1)th‐order water‐bottom‐related multiples as the receiver wavefield, followed by a cross‐correlation imaging condition. Numerical examples from synthetic and field data sets reveal that our approach can provide images with substantially fewer cross‐talk artefacts than conventional reverse time migration using multiples, as well as greatly improving shallow imaging compared with reverse time migration of primaries.

Mono-component multiparameter acoustic full waveform inversion in vertically transverse isotropic media using converted vector wavefields

In multiparameter full-waveform inversion (FWI) with acoustic-approximation for vertically transverse isotropic (VTI) media, it is important to choose an appropriate modeling technique for computational efficiency and numerical stability. In addition, because the gradient is determined by the modeling algorithm used in the FWI process, we need to examine if the modeling algorithm yields a gradient similar to that obtained from the elastic scattering theory. The decomposed equation with mono-component has been proposed because it is computationally efficient and does not generate shear-wave artifacts. However, the decomposed equation has the limitation that its numerical scattering potentials for the perturbations of anisotropic parameters do not properly simulate the theoretical elastic scattering potentials. To overcome the limitation, we propose using the vector virtual sources to compute the gradient in the FWI algorithm by converting the mono-component (pressure) source and receiver wavefields into the two-component particle displacements using the equation of motion on the staggered grids. Unlike the conventional method, our approach properly simulates PP radiation patterns for anisotropic parameters obtained under the elastic assumption, while retaining computational efficiency achieved by using the mono-component. Numerical examples show that anisotropic parameters are updated in the correct gradient direction by our method in the FWI process.

Hybrid multiplicative time-reversal imaging reveals the evolution of microseismic events: Theory and field-data tests

The generation of microseismic events is often associated with induced fractures/faults during the extraction/injection of fluids. A full characterization of the spatiotemporal distribution of microseismic events provides constraints on fluid migration paths in the formations. We have developed a high-resolution source imaging method — a hybrid multiplicative time-reversal imaging (HyM-TRI) algorithm, for automatically tracking the spatiotemporal distribution of microseismic events. HyM-TRI back propagates the data traces from groups of receivers (in space and time) as receiver wavefields, multiplies receiver wavefields between all groups, and applies a causal integration over time to obtain a source evolution image. Using synthetic and field-data examples, we revealed the capability of the HyM-TRI technique to image the spatiotemporal sequence of asynchronous microseismic events, which poses a challenge to standard TRI methods. Moreover, the HyM-TRI technique is robust enough to produce a high-resolution image of the source in the presence of noise. The aperture of the 2D receiver array (azimuth coverage in 3D) with respect to the microseismic source area plays an important role on the horizontal and vertical resolution of the source image. The HyM-TRI results of the field data with 3D azimuthal coverage further verify our argument by producing a superior resolution of the source than TRI.

Least-squares reverse time migration using analytic-signal-based wavefield decomposition

An intrinsic problem during migration and imaging of seismic wavefields using the two-way wave equation is the crosstalk interference between the up/down propagation of the corresponding source and receiver wavefields. To mitigate this crosstalk, the downgoing source and upgoing receiver wavefield imaging condition (IC) is adopted at an early stage of the inversion process, improving convergence and obtaining cleaner reflection images. A wavefield decomposition methodology can also be incorporated into a least-squares reverse time migration (LSRTM) algorithm. The separation of wavefields based on the propagation direction in the early iterations of LSRTM is to reduce interference noise during the inversion process given that the IC considers only primary reflections. Wavefields decomposed with respect to the vertical direction can be easily obtained by Fourier transforms on the time and vertical axes; however, they usually require significantly higher computational effort especially for 3D applications. Vertical wavefield decomposition by a complex-valued analytic signal is an alternative method implemented by the Hilbert transform, which can be conducted by 1D Fourier transform only on the vertical axis. An LSRTM algorithm adopting this decomposition method has a disadvantage in that it requires two additional wave modelings at each iteration. However, by adapting the deprimary IC into LSRTM, only one more modeling is additionally required in the backward wavefield propagation as compared with conventional LSRTM. Our LSRTM using wavefield decomposition has the ability to produce broader band reflectivity images than conventional LSRTM. This is demonstrated with numerical examples using synthetic and real data resulting artifact-free migration results and broadband reflectivity images.

Amplitude-preserving scalar PP and PS imaging condition for elastic reverse time migration based on a wavefield decoupling method

To eliminate crosstalk within the imaging results of elastic reverse time migration (ERTM), we can separate the coupled P- and S-waves from the forward source wavefield and the backpropagated receiver wavefield. The P- and S-wave decoupling method retains the original phase, amplitude, and physical meaning in the separated wavefields. Thus, it is a vital wavefield separation method in ERTM. However, because these decomposed wavefields are vectors, we could consider how to retrieve scalar images that reveal the real reflectivity of the subsurface. For this purpose, we derive a scalar P-wave equation from the velocity-stress relationship for PP imaging. The phase and amplitude of this scalar P-wave are consistent with the scalarized P-wave. Therefore, this scalar P-wave can be exploited to perform PP imaging directly, with the imaging result retaining the amplitude characteristics. For PS imaging, it is difficult to calculate a dynamic preserved scalar S-wave. However, we have developed a scalar PS imaging method that divides the PS image into energy and sign components according to the geometric relationship between the wavefield vibration and propagation directions. The energy is calculated through the amplitude crosscorrelation of the forward P-wave and backpropagated S-wave from the receivers. The sign is obtained from the dot product of the forward P-wave vector and the backpropagated S-wave vector. These PP and PS imaging methods are suitable for 2D and 3D isotropic media and maintain the correct amplitude information while eliminating polarity-reversal phenomena. Several numerical models are used to verify the robustness and effectiveness of our method.

Reverse time adjoint migration

From a mathematical perspective, it is desirable to apply the adjoint operator for back propagating the receiver wavefield for migration of data residuals for inversion. For frequency-domain direct solvers, it is straightforward to apply the adjoint operator, whereas in most real applications of time-domain finite-difference (TDFD) stencils, the forward-propagation kernel function is reused for backward propagation for simplicity. However, when applying the exact adjoint operator, the migration result will be different, especially when strong variations exist in the velocity model. Actually, these three operators (forward, backward, and adjoint of the forward) can be expressed in matrix forms under the Born approximation for acoustic wave equations. These expressions linearly relate model perturbations to recorded seismic data. Every element in the matrices is well-defined, and all involved operations, such as the imaging condition and wavefield sampling, are included in the closed-form matrix expressions. These matrix expressions provide a platform for analyzing seismic modeling and inversion via mature linear algebra methodologies and provide clear strategies for developing computer algorithms. By analyzing the similarity of the matrix expressions, one can find that the time stepping approaches for all three operators are essentially the same. Based on this observation, a new time-marching stencil can be designed to realize the TDFD adjoint operator. Compared with traditional reverse time migration, the new method using the adjoint operator can provide better image quality, especially at sharp velocity boundaries.

Source-independent elastic least-squares reverse time migration

Elastic least-squares reverse time migration (LSRTM) has been developed recently for its high accuracy imaging ability. The theory is based on minimizing the misfit between the observed and simulated data by an iterative algorithm to refine seismic images toward the true reflectivity. We have developed a new elastic LSRTM with the same modeling equations for source and receiver wavefield extrapolations, except for their source terms. The LSRTM has a natural advantage to solve the source and receiver wavefields using the same modeling system; thus, it is easy to implement LSRTM. In practice, it is difficult to obtain an accurate source wavelet, so a convolution-based objective function is used in our source-independent elastic LSRTM. Such an objective function can relax the requirement of an accurate wavelet, and improve the robustness of the inverse problem in the presence of noise. The numerical examples indicate that our method has the ability to recover the reflectivity models with an incorrect source wavelet from noisy data.

An efficient source wavefield reconstruction scheme using single boundary layer values for the spectral element method

In the adjoint-state method, the forward-propagated source wavefield and the backward-propagated receiver wavefield must be available simultaneously either for seismic imaging in migration or for gradient calculation in inversion. A feasible way to avoid the excessive storage demand is to reconstruct the source wavefield backward in time by storing the entire history of the wavefield in perfectly matched layers. In this paper, we make full use of the elementwise global property of the Laplace operator of the spectral element method (SEM) and propose an efficient source wavefield reconstruction method at the cost of storing the wavefield history only at single boundary layer nodes. Numerical experiments indicate that the accuracy of the proposed method is identical to that of the conventional method and is independent of the order of the Lagrange polynomials, the element type, and the temporal discretization method. In contrast, the memory-saving ratios of the conventional method versus our method is at least N when using either quadrilateral or hexahedron elements, respectively, where N is the order of the Lagrange polynomials used in the SEM. A higher memory-saving ratio is achieved with triangular elements versus quadrilaterals. The new method is applied to reverse time migration by considering the Marmousi model as a benchmark. Numerical results demonstrate that the method is able to provide the same result as the conventional method but with about 1/25 times lower storage demand. With the proposed wavefield reconstruction method, the storage demand is dramatically reduced; therefore, in-core memory storage is feasible even for large-scale three-dimensional adjoint inversion problems.

Application of the continuous wavelet transform in the extraction of directional data on RTM imaging condition wavefields

Low-frequency artifacts in reverse time migration result from unwanted cross-correlation of the source and receiver wavefields at non-reflecting points along ray-paths. These artifacts can hide important details in migrated models and increase poor interpretation risk.
 Some methods have been proposed to avoid or reduce the number of these artifacts, preserving reflections, and improving model quality, implementing other strategies such as modification of the wave equation, proposing other imaging conditions, and using image filtering techniques. One of these methods uses wavefield decomposition, correlating components of the wavefields that propagate in opposite directions.
 We propose a method for extracting directional information from the RTM imaging condition wavefields to obtain characteristics allowing for better, more refined imaging. The method works by separating directional information about the wavefields based on the continuous wavelet transform (CWT), and the analysis of the main changes on the frequency content revealed within the scalogram obtained by a Gaussian wavelet family.
 Through numerical applications, we demonstrate that this method can effectively remove undesired artifacts in migrated images. In addition, we use the Laguerre-Gauss filtering to improve the results obtained with the proposed method.

Comparative analysis of 3D RTM Implementation Strategies for an Efficient Use of Memory in a Single GPU:

Reverse-Time Migration (RTM) is a two-way wave-equation based method used to generate images of the Earth’s subsurface. RTM has been successfully used in seismic imaging as it allows defining complex structural areas. However, RTM is a highly computational expensive algorithm that requires the computation of both the source and the receiver wavefields for each shot. Fortunately, numerical methods that compute the wave propagation using the wave equation are highly parallelizable, so they can take leverage on GPU features. However, the main problem of a GPU-RTM implementation is memory management. To take advantage of the GPU computing capabilities, the transfers to host RAM memory storage, or more expensive hard disk storage must be avoided. We present the analysis of three different strategies to implement RTM using only the memory available on a single GPU: (1) Stored wavefield checkpointing (2) Backpropagation of source wavefield using stored boundaries, and (3) Backpropagation of source wavefield using the two last snapshots and random boundaries, showing that the large amount of memory required in the first two strategies becomes a restriction over the model size. The last method (using random boundary conditions) is shown as a suggested solution to the memory problem of using a single GPU.

Viscoacoustic reverse time migration using a time-domain complex-valued wave equation

During seismic wave propagation, intrinsic attenuation inside the earth gives rise to amplitude loss and phase dispersion. Without appropriate correction strategies in migration, these effects degrade the amplitudes and resolution of migrated images. Based on a new time-domain viscoacoustic wave equation, we have developed a viscoacoustic reverse time migration (RTM) approach to correct attenuation-associated dispersion and dissipation effects. A time-reverse wave equation is derived to extrapolate the receiver wavefields, in which the sign of the dissipation term is reversed, whereas the dispersion term remains unchanged. The difference between the forward and time-reverse wave equations is consistent with the physical insights of attenuation compensation during wavefield backpropagation. Due to the introduction of an imaginary unit in the dispersion term, the forward and time-reverse wave equations are complex valued. They are similar to the time-dependent Schrödinger equation, whose real and imaginary parts are coupled during wavefield extrapolation. The analytic property of the extrapolated source and receiver wavefields allows us to explicitly separate up- and downgoing waves. A causal imaging condition is implemented by crosscorrelating downgoing source and upgoing receiver wavefields to remove low-wavenumber artifacts in migrated images. Numerical examples demonstrate that our viscoacoustic RTM approach is capable of producing subsurface reflectivity images with correct spatial locations as well as amplitudes.

Latest field trial confirms potential of new seismic method based on continuous source and receiver wavefields

eSeismic is a novel seismic methodology based on the emission and recording of continuous source and receiver wavefields. One of the motivations behind developing the methodology has been the increased focus on the potential environmental impact of marine seismic acquisition, which the new methodology seeks to reduce. A particular focus has been placed on the peak sound pressure levels emitted from seismic sources and their potential impact on marine mammals and fish with swim bladders. Consequently, authorities across the world have started to introduce stronger regulations concerning the use of seismic sources. The industry has responded by engaging in the development of marine vibrator systems that emit lower-amplitude transient signals and hence are expected to comply with stricter environmental regulations. Different marine vibrator systems are currently being developed or tested but none have reached full-scale commercial readiness. The methodology described in this paper has not been developed with any specific marine source technology in mind. The desired signal for the outlined methodology is that of white noise, as this enables deconvolution of the data with the total source wavefield. Indeed any type of mechanical device that produces a source signal that approaches the properties of white noise can be utilized in the eSeismic method. As described later, existing air gun equipment used on board modern seismic vessels has been used to generate the results discussed here and is equally suited for this method as are any marine vibrator systems of the future.

Making the transition from discrete shot records to continuous seismic records and source wavefields, and its potential impact on survey efficiency and environmental footprint

ABSTRACTA marine seismic method based on continuous source and receiver wavefields has been developed. The method requires continuous recording of the seismic data. The source that may consist of multiple source elements can emit signals continuously while moving. The ideal source wavefield to be used with this method should be as white as possible both in a temporal and a spatial sense to avoid deep notches in the spectrum enabling a stable multi‐dimensional deconvolution. White noise has such properties. However, equipment that can generate white noise does not exist. In order to generate a continuous source wavefield that is approaching the properties of white noise using existing equipment onboard marine seismic vessels, individual air‐guns can be triggered with short randomized time intervals in a near‐continuous fashion. The main potential benefits with the method are to reduce the environmental impact of marine seismic surveys and to improve acquisition efficiency. The peak sound pressure levels are significantly reduced by triggering one air‐gun at a time compared to conventional marine seismic sources. Sound exposure levels are also reduced in most directions. Since the method is based on continuous recording of seismic data and the air‐guns are triggered based on time and not based on position, there are less vessel speed limitations compared to conventional marine seismic data acquisition. Also, because the source wavefield is spread out in time, the wavefields emitted from source elements in different cross‐line positions can be designed such that the emitted wavefield is spatially white in this direction. This means that source elements in multiple cross‐line positions can be operated simultaneously, potentially improving the cross‐line sampling and/or the acquisition efficiency.

Full waveform inversion with angle-dependent gradient preconditioning using wavefield decomposition

The main objective of the full waveform inversion (FWI) is to deliver a velocity model which produces accurate and high resolution depth images. An integral part of FWI is the iterative optimization process aiming at minimizing an objective function describing the misfit between the observed and modelled data with respect to the velocity model parameters. However, FWI commonly suffers from local minima due to errors in data fitting and the starting velocity model. In this study, we develop a methodology to enhance convergence and mitigate the risk of having solutions falling into local minima. We employ optical flow to decompose the source and residual receiver wavefields and map them into the subsurface angle domain. Having evaluated the propagation vectors of the respective cleaner wavefields, subsurface reflection and azimuth angles can be robustly computed. Using such information we can apply angle-dependent preconditioning within the inversion algorithm. During a conventional FWI workflow, larger angles (i.e., transmission wavepaths) are mainly responsible for long wavelength velocity anomalies whereas narrow angles (i.e., reflection wavepaths) influence short wavelength velocity anomalies. In the proposed FWI methodology, we apply an optimized weighting function based on these estimated subsurface reflection angles in order to weight certain angle directions needed for gradient estimation using in parameter updating. Following a hierarchical approach, allowing for larger angles at the early stage of the FWI application, we initially solve for the smooth and longer wavelength velocity updates. Relaxing the opening angle from large to mid to narrow, the gradient direction will first estimate long wavelength velocity model updates (transmission-based wavepaths), following by short wavelength velocity model updates (reflection-based wavepaths). Moreover, angle-dependent-based gradient preconditioning can reduce cycle skips and improve convergence during FWI iterations.

Reverse time migration using phase crosscorrelation

Additional information regarding the continuity and resolution of selected seismic reflectors in reverse time migration (RTM) images can be beneficial for seismic interpretation. We have developed and evaluated new imaging conditions for RTM based on the phase coherence between the forward- and backward-propagated wavefields. These imaging conditions make use of the instantaneous phase and envelope of the analytical signals of the source and receiver wavefields, in addition to their real parts. Once the analytical wavefields are available, these imaging conditions can be calculated simultaneously with conventional conditions at little or no extra cost. The availability of these fields at each image point enables several alternative ways to define imaging conditions. We explore, in addition to pure phase crosscorrelation (PC), two approaches of amplitude-weighted PC. Our numerical experiments, imaging synthetic and field data sets, indicate that these new imaging conditions provide additional images that can highlight some weak reflectors by locally improving the resolution of RTM images. In our examples, this happens particularly in the deep portions of the seismic images. In addition, reflection events produced at discontinuities are enhanced as sharp signals, suggesting that the proposed imaging conditions can help to delineate stratigraphic and structural features that are harder to see in conventional images. These properties of the PC imaging conditions make them an interesting tool to provide additional information that can aid seismic interpretation in complex structural settings.

Reverse-time migration using multidirectional wavefield decomposition method

Reverse-time migration has attracted more and more attention owing to the advantages of high imaging accuracy, no dip restriction, and adaptation to complex velocity models. Cross-correlation imaging method is typically used in conventional reverse-time migration that produces images with strong low-frequency noise. Wavefield decomposition imaging can suppress such noise; however, some residual noise persists in the imaging results. We propose a 2D multidirectional wavefield decomposition method based on the traditional wavefield decomposition method. First, source wavefields and receiver wavefields are separated into eight subwavefields, respectively. Second, cross-correlation imaging is applied to selected subwavefields to produce subimages. Finally, the subimages are stacked to generate the final image. Numerical examples suggest that the proposed method can eliminate the low-frequency noise effectively and produce high-quality imaging profiles.

Directional‐oriented wavefield imaging: a new wave‐based subsurface illumination imaging condition for reverse time migration

ABSTRACTThe key objective of an imaging algorithm is to produce accurate and high‐resolution images of the subsurface geology. However, significant wavefield distortions occur due to wave propagation through complex structures and irregular acquisition geometries causing uneven wavefield illumination at the target. Therefore, conventional imaging conditions are unable to correctly compensate for variable illumination effects. We propose a generalised wave‐based imaging condition, which incorporates a weighting function based on energy illumination at each subsurface reflection and azimuth angles. Our proposed imaging kernel, named as the directional‐oriented wavefield imaging, compensates for illumination effects produced by possible surface obstructions during acquisition, sparse geometries employed in the field, and complex velocity models. An integral part of the directional‐oriented wavefield imaging condition is a methodology for applying down‐going/up‐going wavefield decomposition to both source and receiver extrapolated wavefields. This type of wavefield decomposition eliminates low‐frequency artefacts and scattering noise caused by the two‐way wave equation and can facilitate the robust estimation for energy fluxes of wavefields required for the seismic illumination analysis. Then, based on the estimation of the respective wavefield propagation vectors and associated directions, we evaluate the illumination energy for each subsurface location as a function of image depth point and subsurface azimuth and reflection angles. Thus, the final directional‐oriented wavefield imaging kernel is a cross‐correlation of the decomposed source and receiver wavefields weighted by the illuminated energy estimated at each depth location. The application of the directional‐oriented wavefield imaging condition can be employed during the generation of both depth‐stacked images and azimuth–reflection angle‐domain common image gathers. Numerical examples using synthetic and real data demonstrate that the new imaging condition can properly image complex wave paths and produce high‐fidelity depth sections.

Analytical reverse time migration: An innovation in imaging of infrastructures using ultrasonic shear waves

The emergence of ultrasonic dry point contact (DPC) transducers that emit horizontal shear waves has enabled efficient collection of high-quality data in the context of a nondestructive evaluation of concrete structures. This offers an opportunity to improve the quality of evaluation by adapting advanced imaging techniques. Reverse time migration (RTM) is a simulation-based reconstruction technique that offers advantages over conventional methods, such as the synthetic aperture focusing technique. RTM is capable of imaging boundaries and interfaces with steep slopes and the bottom boundaries of inclusions and defects. However, this imaging technique requires a massive amount of memory and its computation cost is high. In this study, both bottlenecks of the RTM are resolved when shear transducers are used for data acquisition. An analytical approach was developed to obtain the source and receiver wavefields needed for imaging using reverse time migration. It is shown that the proposed analytical approach not only eliminates the high memory demand, but also drastically reduces the computation time from days to minutes.

True amplitude recovery in reverse time extrapolation of plane and spherical waves

A challenging outstanding problem in reverse time extrapolation is recovering accurate amplitudes at reflectors from the receiver wavefield. Various migrations have been developed to produce accurate image locations rather than correct amplitude information because of inadequate compensation of attenuation, dispersion, and transmission losses. We have evaluated the requirements, and determined the theoretical feasibility, of true amplitude recovery of 2D acoustic and elastic seismic data by using the analytic Zoeppritz equations for plane-wave reflection and transmission coefficients. Then, we used synthetic acoustic and elastic wavefield data generated by elastodynamic finite differences to verify the recovery, in the reverse time propagation, of spherical waves and illustrated the salient differences between the incident wavefields reconstructed from reflection data only and from the combination of reflection and transmission data. These examples quantitatively verify that recovering an incident plane or a spherical wave requires the reverse time propagation of all reflections and transmissions in a model with the correct velocity and density. Accurate reconstruction of an incident wave is not possible by backward propagation of only reflections. As an application, we removed downgoing internal multiple reflections generated by upgoing waves incident at reflectors shallower than a horizontal well, in which geophones are deployed. The subtraction of the downgoing reflection involves wavefield reconstruction at depths shallower than the horizontal well and separation of upgoing and downgoing wavefields. This approach assumes that the correct acoustic (or elastic) velocity and density models are available in, and shallower than, the layer where the horizontal well is located. Incident-wave reconstruction works equally well for smooth models, as for models with sharp boundaries. Uncertainties in the model used for reconstruction, and incompleteness of the data aperture are propagated into the equivalent uncertainties, and incompleteness of the reconstruction.

3D elastic full-waveform inversion using P-wave excitation amplitude: Application to ocean bottom cable field data

We have developed an efficient elastic full-waveform inversion (FWI) based on the P-wave excitation amplitude (maximum energy arrival) approximation in the source wavefields. Because, based on the P-wave excitation approximation (ExA), the gradient direction is approximated by the crosscorrelation of source and receiver wavefields at only excitation time, it estimates the gradient direction faster than its conventional counterpart. In addition to this computational speedup, the P-wave ExA automatically ignores SP and SS correlations in the approximated gradient direction. In elastic FWI for ocean bottom cable (OBC) data, the descent direction for the S-wave velocity is often degraded by undesired long-wavelength features from the SS correlation. For this reason, the P-wave excitation approach increases the convergence rate of multiparameter FWI compared with the conventional approach. The modified 2D Marmousi model with OBC acquisition is used to verify the differences between the conventional method and ExA. Finally, the feasibility of the proposed method is demonstrated on a real OBC data from the North Sea.