Time domain viscoelastic full waveform inversion

Summary Viscous attenuation can have a strong impact on seismic wave propagation, but is rarely taken into account in full waveform inversion (FWI). In the time domain, when viscoelasticity is considered, the displacement formulation of the wave equation is usually used. However, the adjoint state equations are quite different for the velocity-stress formulation than for the displacement formulation. In this paper, we derive the adjoint state equations for the viscoelastic wave equation based on the velocity-stress formulation. Using a modified definition of the memory variables commonly found in the literature, we define a Lagrangian from which the adjoint state equations and the misfit gradient are derived. The resulting expressions are similar to the displacement formulation, but differ by the source term and by the wavefield cross-correlations giving the misfit gradient. To validate our results, the misfit gradient obtained by the adjoint state method is compared to the misfit gradient calculated by finite-difference for a simplified problem, giving an excellent agreement. In short, this work gives the right adjoint state equations for the velocity-stress formulation, which is commonly used for time-domain viscoelastic modeling. Further studies are required to evaluate the performance of this approach in real FWI viscoelastic experiments.

In this paper, we study 3-D multiparameter full waveform inversion (FWI) in viscoelastic media based on the generalized Maxwell/Zener body including arbitrary number of attenuation mechanisms. We present a frequency-domain energy analysis to establish the stability condition of a full anisotropic viscoelastic system, according to zero-valued boundary condition and the elastic–viscoelastic correspondence principle: the real-valued stiffness matrix becomes a complex-valued one in Fourier domain when seismic attenuation is taken into account. We develop a least-squares optimization approach to linearly relate the quality factor with the anelastic coefficients by estimating a set of constants which are independent of the spatial coordinates, which supplies an explicit incorporation of the parameter <it>Q</it> in the general viscoelastic wave equation. By introducing the Lagrangian multipliers into the matrix expression of the wave equation with implicit time integration, we build a systematic formulation of multiparameter FWI for full anisotropic viscoelastic wave equation, while the equivalent form of the state and adjoint equation with explicit time integration is available to be resolved efficiently. In particular, this formulation lays the foundation for the inversion of the parameter <it>Q</it> in the time domain with full anisotropic viscoelastic properties. In the 3-D isotropic viscoelastic settings, the anelastic coefficients and the quality factors using bulk and shear moduli parametrization can be related to the counterparts using <it>P</it> and <it>S</it> velocity. Gradients with respect to any other parameter of interest can be found by chain rule. Pioneering numerical validations as well as the real applications of this most generic framework will be carried out to disclose the potential of viscoelastic FWI when adequate high-performance computing resources and the field data are available.

Time-domain full waveform inversion (FWI) in the acoustic regime comprises a parameter identification problem for the acoustic wave equation: pressure waves are initiated by sources, get scattered by the earth’s inner structure, and their reflected parts are picked up by receivers located on the surface. From these reflected wave fields the two parameters, density and sound speed, have to be reconstructed. Mathematically, FWI reduces to the solution of a nonlinear and ill-posed operator equation involving the parameter-to-solution map of the wave equation. Newton-like iterative regularization schemes are well suited and well analyzed to tackle this inverse problem. Their convergence results are often based on an assumption about the nonlinear map known as tangential cone condition. In this paper we verify this assumption for a semi-discrete version of FWI where the searched-for parameters are restricted to a finite dimensional space. As a byproduct we establish that the semi-discrete seismic inverse problem is locally Lipschitz stable, in particular, it is conditionally well-posed.

This study aims at designing optimal crosshole seismic approach for CO2EOR (enhance oil recovery by CO2 injection) monitoring in Abu Dhabi. It is not easy to perform successful reservoir monitoring for carbonate reservoirs in the Middle East because rock properties may not change much associated with fluid replacement. We investigate advantage of newly available elastic full waveform inversion (FWI) for detection of CO2 flooded zones. To investigate advantage, applicability and future problem of time-lapse FWI for CO2EOR monitoring, we conducted rock physics study, elastic seismic simulation, and test application of FWI to synthetic data. We first made realistic velocity models based on Gassmann's equation and actual well logs obtained before and after CO2 injection at carbonate reservoirs in Abu Dhabi. The in-situ light oil at the target reservoirs was replaced by CO2 (Secondary mode CO2EOR). We prepared models with different width and thickness of CO2 flooded zones and analyzed detectability of seismic waveform change. We also conducted test application of elastic full waveform inversion to synthetic data for delineation of CO2 flooded zones. Results of the rock physics study show 1.6% decrease in P-wave velocity, little change in S-wave velocity and density associated with secondary mode CO2EOR. The small property changes are caused by the rigid rock frame and small density difference in light oil and supercritical CO2 at the reservoir condition. Next we performed synthetic seismic analysis based on the rock property models, available crosshole seismic source and realistic well distances. Because of the small property change, seismic traveltime difference and waveform change are quite small when well distance is 116ft. When well distance is 700 ft or greater, guided waves, or trapped energy inside low velocity layers can be observed; relatively large waveform change can be observed in the guided waves due to CO2 injection. Although guided waves have not been used for reservoir monitoring, they are sensitive to small Vp and thus can be an indicator of the CO2 flooding. The elastic FWI result for a model with well distance of 1160 ft demonstrated good representation of the differential Vp with vertical resolution of 20-30 ft and some horizontal resolution. Although we need further investigation on repeatability and signal transmittability using actual field data, our synthetic study indicates possibility of identifying layers with CO2 using guided waves and elastic FWI.

Time-domain multi-scale full waveform inversion (FWI) is an important method for seismic imaging. Using different frequency contents of seismic data, this method mitigates the effect of the local-minima problem by sequential inversion. One of the important factors determining the convergence rate of the nonlinear inversion is the precondition of the gradient. Virtual-source precondition has been successfully applied to Laplace-domain and Laplace-Fourier-domain FWI. In this paper, we will examine the application of the virtual-source precondition to time-domain multi-scale FWI. The research results reveal that, compared to the gradient method without precondition, the virtual-source-based gradient method achieves much greater convergence rate while its additional computational cost is very low. Therefore, virtual-source precondition is very effective for time-domain multi-scale FWI. keywords—virtual-source; precondition; full-waveform inversion

In acoustic logging, the dipole-mode wave is a type of guided wave that propagates along the borehole, and its dispersion characteristics are typically used to invert the shear-wave velocity of the formation around the borehole. However, traditional dispersion-based inversions rely on the layered model assumption, limiting their applicability to formations that are either uniformly distributed along the borehole axis or exhibit gradual variations. As a method that directly fits observed waveforms, full waveform inversion (FWI) can be applied to various formation models without being constrained by the assumption of two-dimensional homogeneity. Therefore, this technique is introduced into the processing of array acoustic logging data to reconstruct the shear-wave velocity structure. Finite-difference methods are used for forward modeling and the adjoint equations and gradient expressions of the elastic wave equation in cylindrical coordinates for dipole sources are derived. To address the challenges posed by the strong gradients generated by guided waves near the borehole wall, the gradients are preprocessed with an energy-weighted method. Additionally, an initial model-building method based on the slowness-time-coherence method and dispersion inversion are proposed, providing a high-accuracy initial velocity model. Numerical tests show that the FWI results match the true model, accurately reconstructing the shear-wave velocity structure.

An accurate earth model is key to any successful depth imaging effort. Full waveform inversion is an advanced velocity model building process which uses the full two way wave equation. Existing methods use a ray tracing approach to distribute velocity errors within the model. In this presentation, we will show examples of seismic data processed with the latest technology, including earth model building with full waveform inversion. Introduction The industry has moved to using two-way wave equation migrations(commonly known as reverse time migration), especially in areas of complex geology such as the salt bodies in the Gulf of Mexico and offshore Brazil. The velocity model, including velocity anisotropy, is key to any depth migration effort. The natural next step is to use the two-way wave equation for velocity model building. One of the most advanced tools for velocity model building using the two-way wave equation is Full Waveform Inversion(FWI). Full-waveform inversion uses computer intensive forward modeling of the seismic measurement combined with residual wavefield back propagation to iterate to a final velocity model, which can provide greater detail than tomographic ray tracing approaches. Velocity inversion using full waveform inversion has been used on a number of projects around the world including the GOM, Brazil and the North Sea.

Time-domain full waveform inversion (FWI) in the acoustic regime comprises a parameter identification problem for the acoustic wave equation: Pressure waves are initiated by sources, get scattered by the earth’s inner structure, and their reflected parts are picked up by receivers located on the surface. From these reflected wave fields the two parameters, density and sound speed, have to be reconstructed. Mathematically, FWI reduces to the solution of a nonlinear and ill-posed operator equation involving the parameter-to-solution map of the wave equation. Newton-like iterative regularization schemes are well suited and well analyzed to tackle this inverse problem. Their convergence results are often based on an assumption about the nonlinear map known as tangential cone condition. In this paper we verify this assumption for a semi-discrete version of FWI where the searched-for parameters are restricted to a finite dimensional space. As a byproduct we establish that the semi-discrete seismic inverse problem is Lipschitz stable, in particular, it is conditionally well-posed.

Velocity models presenting sharp interfaces are highly relevant in seismic imaging, e.g., for imaging the subsurface of the Earth in the presence of salt bodies. In order to mitigate the oversmoothing of classical regularization strategies such as the Tikhonov regularization, we propose a shape optimization approach for the sharp-interface reconstruction in time-domain acoustic full waveform inversion. Our main result includes the shape differentiability of the cost functional measuring the misfit between observed and predicted data. In particular, it reveals the expression of the distributed shape derivative in tensor form, built on a Lagrangian-type approach and regularity results for the wave equation with discontinuous coefficients. Based on the achieved distributed shape derivative and the level set method, we propose a numerical approach and present several numerical tests supporting our approach.

Velocity macromodel building is an essential step of the seismic imaging workflow. Indeed, obtaining acceptable results through migration or full waveform inversion is highly dependent on the kinematic accuracy of the background/initial velocity model. Two decades ago, stereotomography was proposed as an alternative to reflection traveltime tomography, the first relying on semi-automatic picking of locally coherent events associated with small reflection or diffraction segments tied to scatterers in depth by a pair of rays, while the latter on interpretive picking of laterally continuous reflections. The flexibility of stereotomography paved the way for many developments that have shown the efficiency of the method whilst emphasizing on the complementary information carried out by traveltimes and slopes of locally coherent events. A recent formulation recast stereotomography under a matrix-free formulation based on eikonal solvers and the adjoint-state method. In the latter, like in the previous works, the scatterer positions and the velocity field are updated jointly to tackle the ill-famed velocity–position coupling in reflection tomography. Following on from this adjoint-state formulation, we propose a new parsimonious formulation of slope tomography that offers the chance to restrain the problem to minimizing the residuals of a single data class being a slope, in search of a sole parameter class being the subsurface velocity field. This parsimonious formulation results from a variable projection, which is implemented by enforcing a consistency between the scatterer coordinates and the velocity macromodel through migration of kinematic attributes. We explain why the resulting reduced-parametrization inversion is more suitable for tomographic problems than the most common joint inversion strategy. We benchmark our method against the complex Marmousi model along with a validation through time domain full waveform inversion and then present the results of a field data case study.

Glaciers, which store nearly 70% of the Earth's freshwater, are undergoing significant changes due to accelerating melting caused by climate change. A better understanding of their behavior and mechanisms is therefore crucial for the years to come. To study these processes, a dense seismic array was deployed on the Argentière Glacier in the Mont Blanc massif (French Alps). The sensor array consists in 98 3-component seismic stations which continuously recorded surface displacements over 35 days in early spring 2018. This period coincided with rising temperatures and rapid glacier evolution, so the sensors captured thousands of events, mostly icequakes.The aim of this study is to reconstruct the glacier's structure and study icequake mechanisms using elastic Full Waveform Inversion (FWI) on the 3-component data. As the data come from a passive seismic experiment, we have no information about the sources. Before reconstructing the structure it is therefore necessary to work on source parameters estimation. These parameters include spatial localization and mechanism. First, we detect and localize the icequakes inside the glacier using a beamforming method called Matched Field Processing (MFP). The detected icequakes are observed to be located mostly at the positions of crevasses on the glacier's surface. Then, we decompose the icequake mechanism into a moment tensor and a time signature wavelet. To estimate these two parameters, we develop a joint inversion method based on waveform analysis using an iterative alternating minimization algorithm. The type of mechanism and the source orientation are then interpreted through the eigenvalue and eigenvector decomposition of the estimated moment tensor. The Fundamental Lune representation is employed to statistically study the distribution of more than 14,000 icequake mechanisms within the glacier, revealing a significant proportion of opening and closing mechanisms associated with crevasses. In certain areas, Double-Couple (DC) mechanisms can also be observed, potentially corresponding to crevasse fault slip events.Using the estimated source parameters, FWI can be applied to reconstruct the glacier structure. A 3D synthetic crevasse model was created to mimic reality, incorporating the three observed crevasse clusters on the glacier, to evaluate the effectiveness of FWI in reconstructing the model in a given frequency-band. The parameters used in the model include S-wave and P-wave velocities, as well as density. The inversion results reveal several key findings: first, multi-parameter inversion with both S-wave velocity and density yields better results. Second, crevasses can be accurately reconstructed within the considered frequency band, provided the source parameters are well-estimated. Finally, the accuracy of source mechanism estimation significantly impacts the quality of crevasse reconstruction. Importantly, iterating between mechanism estimation and structure reconstruction yields improved results, providing promising insights.

To improve the understanding of flow and transport processes in the critical zone, high-resolution and accurate estimation of the small-scale heterogeneity is essential. Preferential flow paths related to high-porosity layers and clay lenses in gravel aquifers greatly affect flow and transport processes in the subsurface, and their high electrical contrast to their surrounding matrix and limited extent can act as low-velocity electromagnetic waveguides. In the past decade, time-domain full-waveform inversion (FWI) of crosshole ground-penetrating radar (GPR) data has shown to provide 2D decimeter-scale resolution images of relative permittivity and electrical conductivity of the subsurface, which can be related to porosity and soil texture. Most studies using crosshole GPR FWI resolved high-porosity zones that were identified by an amplitude analysis approach. But clay lenses or zones with higher electrical conductivity that act as low-velocity waveguides are hard to distinguish in the measured data and amplitude analysis because of the absence of characteristic wave-propagation features. We have investigated a set of nine crosshole GPR data sets from a test site in Hermalle-sous-Argenteau near the Meuse River in Belgium to characterize the aquifer within a decimeter-scale resolution and to improve the understanding of a previously performed heat tracer experiment. Thereby, we extend the amplitude analysis to identify two different types of low-velocity waveguides either caused by an increased porosity or a higher electrical conductivity (and higher porosity). Combining the GPR amplitude analysis for low-velocity waveguide zones with the standard FWI results provided information on waveguide zones, which modified the starting models and further improved the FWI results. Moreover, an updated effective source wavelet is estimated based on the updated permittivity starting models. In comparison with the traditional FWI results, the updated FWI results present smaller gradient of the medium properties and smaller root-mean-squared error values in the final inversion results. The nine crosshole sections are used to generate a 3D image of the aquifer and allowed a detailed analysis of the porosity distribution along the different sections. Consistent structures of the permittivity and electrical conductivity show the robustness of the updated FWI results. The aquifer structures obtained by the FWI results agree with those results of the heat tracer experiment.

Full-waveform inversion (FWI) includes migration and tomography modes. The tomographic component of the gradient from reflection data is usually much weaker than the migration component. To use the tomography mode to fix background velocity errors, it is necessary to extract the tomographic component from the gradient. Otherwise, the inversion will be dominated by the migration mode. We have developed a method based on nonstationary smoothing to extract the tomographic component from the raw gradient. By analyzing the characteristics of the scattering angle filtering, the wavenumber of the tomographic component at a given frequency is seen to be smaller than that of the migration component. Therefore, low-wavenumber-pass filtering can be applied to extract the tomographic component. The low-wavenumber-pass smoothing filters are designed with Gaussian filters that are determined by the frequency of inversion, the model velocity, and the minimum scattering angle. Thus, this filtering is nonstationary smoothing in the space domain. Because this filtering is carried out frequency by frequency, it works naturally and efficiently for FWI based on frequency-domain modeling. Furthermore, because the maximum opening angle of the reflections in a typical acquisition geometry is much smaller than the minimum scattering angle for the tomographic component, which is generally set at 160°, there is a relatively large gap between the wavenumbers of the tomographic and migration components. In other words, the nonstationary smoothing can be applied once to a group of frequencies for time-domain FWI without leaking the migration component into the tomographic component. Analyses and numerical tests indicate that two frequency groups are generally sufficient to extract the tomographic component for the typical frequency range of time-domain FWI. The numerical tests also demonstrate that the nonstationary smoothing method is effective and efficient at extracting the tomographic component for reflection waveform inversion.

Full-waveform inversion (FWI) is firmly established within our industry as a powerful velocity model building tool. FWI carries significant theoretical advantages over conventional velocity model building methods such as refraction and reflection tomography. Specifically, by solving a nonlinear inverse problem through the wave equation, FWI is able to recover a broadband velocity model containing both high and low spatial wavenumbers, thus extending the approximation of residual moveout correction inherent in traditional velocity model building approaches. Moreover, FWI is capable of inverting information from the entire wavefield (i.e., early arrivals, reflections, refractions, and multiple energy) rather than from a subset as in conventional approaches (i.e., first break and primary reflections), thereby availing itself of more information to better constrain its model estimate. However, these theoretical benefits cannot be realized easily in practice because various complexities of real seismic data often conspire to violate algorithmic assumptions, leading to unsatisfactory results. Dynamic matching FWI (DMFWI) is a newly developed algorithm that solves an inversion problem that maximizes the cross correlation of two dynamically matched data sets — one recorded and the other synthetic. Dynamic matching of the two data sets de-emphasizes the amplitude impact, which allows the algorithm to focus on minimizing their kinematic differences rather than amplitude in the data-fitting process. The multichannel correlation makes the algorithm robust for data with low signal-to-noise ratio. Applications of DMFWI across different types of acquisition and geologic settings demonstrate that this novel FWI approach can resolve complex velocity errors and provide high-quality migrated images that exhibit a high degree of geologic plausibility. Additionally, reflectivity images can be obtained in a straightforward manner as natural byproducts through computation of the directional derivative of the inverted FWI velocity models.

In order to decrease the fluid drag on an underwater robot manipulator, an optimal trajectory method based on the variational method is presented. By introducing the adjoint variables, which are Lagrange multipliers, we formulate a Lagrange function under certain constraints related to the target angle, target angular velocity, and dynamic equation of the robot manipulator. The state equation (the partial differentiation of the Lagrange function with respect to the state variables), adjoint equation (the partial differentiation of the Lagrange function with respect to the adjoint variables), and sensitivity equation (the partial differentiation of the Lagrange function with respect to torques) can be derived from the stationary conditions of the Lagrange function. Using the state equation, we can calculate the state variables (angles, angular velocities, and angular acceleration) at every time step in the forward time direction. These state variables are stored as data at every time step. Next, by using the adjoint equation, we can calculate the adjoint variables by using these state variables at every time step in the backward time direction. These adjoint variables are stored as data at every time step. Third, the sensitivity equation is calculated by using both the state variables and the adjoint variables. Finally, the optimal trajectory of the manipulator is obtained using the sensitivities. The proposed method is applied to the problem of two-link manipulators. It can obtain the optimal drag reduction trajectory of the manipulator under the constraints mentioned above.