Time-domain Full Waveform Inversion Research Articles

Sound speed imaging of cortical bone and soft tissue based on frequency-domain ultrasound waveform tomography

Abstract As one of the ultrasound computed tomography (UCT) algorithms, full-waveform inversion (FWI) can recover clear and detailed internal structures of bones by iteratively minimizing the misfit between observed and numerical data. However, time domain full-waveform inversion (TDFWI) often lacks sufficient low-frequency information, and the computation of full-time data imposes a significant computational burden. Besides, there are few studies on whole-medium reconstruction that jointly invert the bone and surrounding soft tissues. In this study, an FWI-based frequency-domain UCT (FD-UCT) is utilized for simultaneously imaging the sound speed of bones and soft tissues, which facilitates the inversion efficiency while ensuring the existence of low-frequency information. Simulated studies show that the FD-UCT can both reconstruct bones with large acoustic impedance and distinguish subtle impedance differences within soft tissues to reveal lesions.

Improved Earthquake Source Parameters with 3D Wavespeed Models in California and Nevada

Abstract Seismic tomography harnesses earthquake data to explore the inaccessible structure of the Earth. Adjoint waveform tomography (AWT), a method of seismic tomography, updates the tomographic model by optimizing the fit between observed earthquake data and synthetic waveforms. The synthetic data are calculated by solving the wave equation through a given 3D model. An important requirement to calculating synthetics is the source information (location, centroid time, depth, and moment tensor). Errors in source information affect the quality of the synthetics produced, which in turn can limit how structure can be inferred in the AWT workflow. To test the effect of updating source information, we used MTTime (Chiang, 2020), a time-domain full-waveform moment tensor inversion code, to calculate the moment tensors and depths of 118 earthquakes that occurred in California and Nevada over a 20-yr period. We calculated 3D Green’s functions using a 3D seismic wavespeed model of California and Nevada (Doody et al., 2023b). We show that the inverted solutions provide better waveform fits than the Global Centroid Moment Tensor catalog and increase usable, well-correlated data by up to 7%. Therefore, we argue that recalculating source parameters should be considered in AWT workflows, particularly for smaller magnitude events (Mw<5.0).

ShaVi-1.0: An interface enabled open source 2D acoustic full waveform inversion package

Full waveform inversion (FWI) is an approach of seismic data inversion that employs data-fitting techniques to produce a detailed image of the subsurface. The computational expense and complex mathematics involved in implementing FWI make it challenging. To address this issue, we have developed an interface enabled 2D acoustic time domain FWI package: ShaVi-1.0 (“Sha” refers to “Shabda”, a Sanskrit word meaning “sound”, and “Vi” stands for “Viparita”, representing “inversion”, the package version is labeled as 1.0). That utilizes shots parallelization on the OpenMP API to deal with computational expenses and a graphical user interface (GUI) for its easy implementation. This allows the FWI to run on standalone workstations. The base code of FWI is written in “C” language with a modular design, enabling users to modify the modules as per their requirements. The GUI, which is developed in Python, includes features such as real-time plotting and process status updates, along with message and warning boxes to ensure smooth inversion flow. We add a toolbox to the GUI’s menu bar along with the help, which users can use to generate acquisition geometry and plot source-receiver locations. The toolbar integrated into the software provides multiple other functionalities to enhance the ease of implementation. We perform testing on the benchmark models such as SEG/EAGE Overthrust and Foothill model. A comparative study is also performed to assess its efficiency relative to CSIM-KAUST and SWIT-1.0 codes. The assessment is based on two key parameters: mean absolute percentage error (MAPE and computational time. The developed GUI not only simplifies the implementation process but also enables real-time analysis of results.

Local uncertainty quantification for 3-D time-domain full-waveform inversion with ensemble Kalman filters: application to a North Sea OBC data set

SUMMARY Full-waveform inversion (FWI) has emerged as the state-of-the art high resolution seismic imaging technique, both in seismology for global and regional scale imaging and in the industry for exploration purposes. While gaining in popularity, FWI, at an operational level, remains a heavy computational process involving the repeated solution of large-scale 3-D wave propagation problems. For this reason it is a common practice to focus the interpretation of the results on the final estimated model. This is forgetting FWI is an ill-posed inverse problem in a high dimensional space for which the solution is intrinsically non-unique. This is the reason why being able to qualify and quantify the uncertainty attached to a model estimated by FWI is key. To this end, we propose to extend at an operational level the concepts introduced in a previous study related to the coupling between ensemble Kalman filters (EnKFs) and FWI. These concepts had been developed for 2-D frequency-domain FWI. We extend it here to the case of 3-D time-domain FWI, relying on a source subsampling strategy to assimilate progressively the data within the Kalman filter. We apply our strategy to an ocean bottom cable field data set from the North Sea to illustrate its feasibility. We explore the convergence of the filter in terms of number of elements, and extract variance and covariance information showing which part of the model are well constrained and which are not. Analysing the variance helps to gain insight on how well the final estimated model is constrained by the whole FWI workflow. The variance maps appears as the superposition of a smooth trend related to the geometrical spreading and a high resolution trend related to reflectors. Mapping lines of the covariance (or correlation matrix) to the model space helps to gain insight on the local resolution. Through a wave propagation analysis, we are also able to relate variance peaks in the model space to variance peaks in the data space. Compared to other posterior-covariance approximation scheme, our combination between EnKF and FWI is intrinsically scalable, making it a good candidate for exploiting the recent exascale high performance computing machines.

Frequency-domain full-waveform inversion-based musculoskeletal ultrasound computed tomography.

Recently, full-waveform inversion (FWI) has become a promising tool for ultrasound computed tomography (USCT). However, as a computationally intensive technique, FWI suffers from computational burden, especially in conventional time-domain full-waveform inversion (TDFWI). On the contrary, frequency-domain full-waveform inversion (FDFWI) provides a relatively high computational efficiency as the propagation of discrete frequencies is much cheaper than full time-domain modeling. FDFWI has already been applied in soft tissue imaging, such as breast, but for the musculoskeletal model with high impedance contrast between hard and soft tissues, there is still a lack of an effective source estimation method. In this paper, a water-referenced data calibration method is proposed to address the source estimation challenge in the presence of bones, which achieves consistency between the measured and simulated data before the FDFWI procedure. To avoid the cycle-skipping local minimum effect and facilitate the algorithm convergence, a starting frequency criterion for musculoskeletal FDFWI is further proposed. The feasibility of the proposed method is demonstrated by numerical studies on retrieving the anatomies of the leg models and different musculoskeletal lesions. The study extends the advanced FDFWI method to the musculoskeletal system and provides an alternative solution for musculoskeletal USCT imaging with high computational efficiency.

Full waveform inversion guided wave tomography with a recurrent neural network

Corrosion quantitative detection of plate or plate-like structure materials is crucial in industrial Non-Destructive Testing (NDT) for determining their remaining life. For doing that, a novel ultrasonic guided wave tomography method, incorporating recurrent neural network (RNN) into full waveform inversion (FWI) called as RNN-FWI, is proposed in this paper. When the wave equation of an acoustic model is solved by a forward model with the cyclic calculation units of an RNN, it is shown that the inversion of the forward model can be obtained iteratively by minimizing a waveform misfit function of quadratic Wasserstein distance between the modeled and measured data. It is also demonstrated that the gradient of the objective function can be obtained by automatic differentiation while the parameters of the waveform velocity model are updated by the adaptive momentum estimation algorithm (Adam). The U-Net deep image prior (DIP) is used as the velocity model regularization in each iteration. The final thickness maps of the plate or plate-like structure materials shown can be archived by the dispersion characteristics of guided waves. Both the numerical simulation and experimental results show that the proposed RNN-FWI tomography method performs better than the conventional time-domain FWI in terms of convergence rate, initial model requirement, and robustness.

Three-dimensional time-domain full-waveform inversion for ring-array-based ultrasound computed tomography

Ultrasound computed tomography (USCT) is a promising imaging technique for breast cancer diagnosis. Full-waveform inversion (FWI)-based image reconstruction methods based on the acoustic wave equation can produce quantitatively accurate images of acoustic properties with higher spatial resolution compared to ray-based methods. One common USCT design employs a circular ring-array of elevation-focused ultrasonic transducers that can be translated orthogonally to the imaging plane to achieve volumetric scanning. Currently, slice-by-slice (SBS) reconstruction methods are often used in this setting. Such methods use a 2D wave physic model to reconstruct a cross-sectional image for each position of the ring-array and vertically stack them together to render the 3D volume. However, this approach neglects the 3D wave propagation physics and transducer’s focusing properties, leading to out-of-plane scattering artifacts and degraded spatial resolution. To address this, a 3D time-domain FWI method that utilizes measurement data acquired at multiple adjacent ring-array locations is proposed. The reconstruction method was investigated using virtual imaging studies of ring-array-based USCT that employed realistic 3D numerical breast phantoms. The results demonstrated improved spatial resolution (both in-plane and off-plane) and reduced image artifacts compared to the SBS approach. Additionally, the impact of the number of ring-array locations on image quality is assessed.

Waveform Inversion of Shallow Seismic data with Randomly Selected Sources

Applying waveform inversion with randomly selected sources (RS) increases the convergence rate of the optimization process within full-waveform inversion (FWI) workflows and reduces overall computational time. FWI has been shown to be a valuable addition to the existing geophysical methods for near-surface characterization. Accurate 3D modeling of the (visco) elastic wavefield allows to diminish assumptions about wave propagation and include surface- and body-wave based arrivals within the inversion workflow. This approach could result in reliable high-resolution subsurface models, but generally it comes at a high computational cost for each individual source modeling. Commonly, multiple sources are involved in the inversion process, which proportionally increases the resulting computations for a time-domain FWI, but seismic data with dense source/receiver coverage usually contain redundant information. This is especially true for seismic near-surface applications, where the number of recorded sources per wavelength of interest are normally excessive. The inversion performance was increased by randomly selecting a subset of sources at each FWI iteration. The method's effectiveness is obvious on a FWI near-surface void detection application. Synthetic 2D experiments for fixed and rolling spreads showed comparable results with fewer calculations. The best performance was achieved when a single random source was used for each inversion iteration. The effectiveness of the method was also evident on a shallow 3D field dataset collected in the Sonora Desert in western Arizona, where data were acquired over a 10-m deep void with known location, orientation, and dimensions.

Time-Domain Elastic Full Waveform Inversion With Frequency Normalization

Time-domain elastic full waveform inversion (FWI) uses seismic data to recover the high-resolution subsurface medium properties for structural imaging and lithologic identification. The approximate pulse generated by the explosion source in seismic exploration evolves into a limited bandwidth seismic wavelet through the propagation of the seismic wave in the underground medium, exhibiting strong energy near the dominant frequency and weak energy far away. Therefore, time-domain FWI using the band-limited seismic wavelet can only match the energy of the dominant frequencies present in a dataset to a large extent, resulting in insufficient low-wavenumber updates of model parameters because of the weak energy of low frequencies. To remove the effect of finite-frequency-band wavelet spectra from the time-domain FWI, we propose a time-domain elastic FWI without wavelet spectral limitation. In our method, a newly refined seismic wavelet with a normalized amplitude and accurate phase is used to propagate seismic waves in the time domain. Data residual measurement is performed based on the summation of single-frequency residuals between the normalized synthetic and observed frequencies. The adjoint-state method is used to approximate the gradients of the model parameters, and the decoupled wavefields obtained by the phase-sensitive detection method were involved in the gradient calculation. Overall, the proposed method helps FWI avoid falling into local minima by enhancing the low wavenumber reconstruction of the velocity models. The elastic FWI numerical tests and field data FWI application demonstrate that our approach can reliably recover high-precision inversion results.

Spyro: a Firedrake-based wave propagation and full-waveform-inversion finite-element solver

Abstract. In this article, we introduce spyro, a software stack to solve wave propagation in heterogeneous domains and perform full waveform inversion (FWI) employing the finite-element framework from Firedrake, a high-level Python package for the automated solution of partial differential equations using the finite-element method. The capability of the software is demonstrated by using a continuous Galerkin approach to perform FWI for seismic velocity model building, considering realistic geophysics examples. A time domain FWI approach that uses meshes composed of variably sized triangular elements to discretize the domain is detailed. To resolve both the forward and adjoint-state equations and to calculate a mesh-independent gradient associated with the FWI process, a fully explicit, variable higher-order (up to degree k=5 in 2D and k=3 in 3D) mass-lumping method is used. We show that, by adapting the triangular elements to the expected peak source frequency and properties of the wave field (e.g., local P-wave speed) and by leveraging higher-order basis functions, the number of degrees of freedom necessary to discretize the domain can be reduced. Results from wave simulations and FWIs in both 2D and 3D highlight our developments and demonstrate the benefits and challenges with using triangular meshes adapted to the material properties.

Seismic Imaging of Complex Velocity Structures by 2D Pseudo-Viscoelastic Time-Domain Full-Waveform Inversion

In the presented study, multi-parameter inversion in the presence of attenuation is used for the reconstruction of the P- and the S- wave velocities and the density models of a synthetic shallow subsurface structure that contains a dipping high-velocity layer near the surface with varying thicknesses. The problem of high-velocity layers also complicates selection of an appropriate initial velocity model. The forward problem is solved with the finite difference, and the inverse problem is solved with the preconditioned conjugate gradient. We used also the adjoint wavefield approach for computing the gradient of the misfit function without explicitly build the sensitivity matrix. The proposed method is capable of either minimizing the least-squares norm of the data misfit or use the Born approximation for estimating partial derivative wavefields. It depends on which characteristics of the recorded data—such as amplitude, phase, logarithm of the complex-valued data, envelope in the misfit, or the linearization procedure of the inverse problem—are used. It showed that by a pseudo-viscoelastic time-domain full-waveform inversion, structures below the high-velocity layer can be imaged. However, by inverting attenuation of P- and S- waves simultaneously with the velocities and mass density, better results would be obtained.

Full-Waveform Inversion in Time Domain by Increased Frequency Iteratively

With seismic imaging techniques using numerical solutions of seismic acoustic wave equations, full-waveform inversion in the time domain can be a tool for estimating high-resolution models from complex geological conditions. However, full-waveform inversion in the time domain has several drawbacks, namely if you determine the initials of the model that is far from the global minimum (the true model) will give a stacked result at the local minimum and to overcome this problem is to use the method of increasing the frequency iteratively. With this method, at least it will keep away the results from a local minimum and close to a global minimum. The result of using the iterative frequency increment method on the full-waveform inversion time domain produces an image that resembles the original synthetic model and produces an updated velocity value that is not too different from the true velocity value of the synthetic model.

Indirect joint petrophysical inversion of synthetic shallow-seismic and multi-offset ground-penetrating radar data

SUMMARYBoth full-waveform inversion (FWI) of shallow-seismic and ground-penetrating radar (GPR) surface recordings have received particular attention in the past decade since they can reconstruct seismic and electromagnetic properties at high resolution. For consistent near-surface imaging, seismic and GPR wavefields can be combined by joint petrophysical inversion (JPI) using classical time-domain FWI. In conventional JPI of wavefields, both shallow-seismic and multi-offset GPR data contribute to reconstructing the same petrophysical parameters. In this paper, we show that seismic and GPR wavefields have different sensitivities to these parameters assuming the widely established petrophysical model combining the Gassmann fluid substitution model, the complex refractive index model and Archie’s equation. Based on this observation, we propose a new coupling strategy of petrophysical parameters which we call indirect JPI. In indirect JPI, seismic data is primarily used for porosity reconstruction, while GPR data is used only for saturation reconstruction. Unlike conventional JPI, we first update the seismic and GPR parameters using non-petrophysical parametrizations and then transform the most reliable estimates to petrophysical parameters. 2-D synthetic tests show that indirect JPI can provide more accurate and consistent results than conventional JPI. In addition, indirect JPI is more robust when uncertainties exist in petrophysical a priori knowledge. More importantly, indirect JPI has the flexibility to integrate different types of seismic and electromagnetic waves and acquisition geometries depending on the target of interest, resulting in the best solution. Indirect JPI has been proven to be a promising approach for multiparameter reconstructions (seven parameters in this study).

On the time-domain full waveform inversion for time-dissipative and dispersive poroelastic media

This paper concerns the Time-Domain Full Waveform Inversion (FWI) for dispersive and dissipative poroelastic materials. The forward problem is an initial boundary value problem (IBVP) of the poroelastic equations with a memory term; the FWI is formulated as a minimization problem of a least-square misfit function with the (IBVP) as the constraint. In this paper, we derive the adjoint problem of this minimization problem, whose solution can be applied to computed the direction of steepest descent in the iterative process for minimization. The adjoint problem has a similar numerical structure as the forward problem and hence can be solved by the same numerical solver. Because the tracking of the energy evolution plays an important role in the FWI for dissipative and dispersive equations, the energy analysis of the forward system is also carried out in this paper.

Elastic wave full-waveform inversion in the time domain by the trust region method

In this paper, we investigate the elastic wave full-waveform inversion (FWI) by the trust region method. FWI is an optimization problem of minimizing the misfit between the observed data and simulated data. Usually the gradient-based search methods such as the steepest descent method or quasi-Newton's methods are applied to update the model parameters iteratively. And the L-BFGS method with approximated Hessian information is preferred for its high computational efficiency and accuracy. At each iteration, a line search method computes a search direction and then finds a suitable step length or how far to move along the direction. The trust region method is an effective alternative to solve nonlinear optimization problems. In this paper, we investigate the performance of the trust region method in elastic time-domain FWI. In the trust region method, a trial step length called the trust-region radius is forced within a certain neighborhood of the current iterate point and a trust-region subproblem requires to solve. In order to solve the trust-region subproblem, the dogleg method and the two dimensional subspace method are adopted in this paper. In the trust region method, as long as the trust-region radius is well updated, the model updating can be performed for every iteration with the fixed trial step and there is no more extra computation for forward problem like the line search method. Numerical computations for the benchmark Marmousi model are given. The inversion results and comparisons show that the trust region method is very efficient and behaves better than the line search method such as the well-known L-BFGS method.

High resolution 3D ultrasonic breast imaging by time-domain full waveform inversion

Ultrasound tomography (UST) scanners allow quantitative images of the human breast’s acoustic properties to be derived with potential applications in screening, diagnosis and therapy planning. Time domain full waveform inversion (TD-FWI) is a promising UST image formation technique that fits the parameter fields of a wave physics model by gradient-based optimization. For high resolution 3D UST, it holds three key challenges: firstly, its central building block, the computation of the gradient for a single US measurement, has a restrictively large memory footprint. Secondly, this building block needs to be computed for each of the 103–104 measurements, resulting in a massive parallel computation usually performed on large computational clusters for days. Lastly, the structure of the underlying optimization problem may result in slow progression of the solver and convergence to a local minimum. In this work, we design and evaluate a comprehensive computational strategy to overcome these challenges: firstly, we exploit a gradient computation based on time reversal that dramatically reduces the memory footprint at the expense of one additional wave simulation per source. Secondly, we break the dependence on the number of measurements by using source encoding (SE) to compute stochastic gradient estimates. Also we describe a more accurate, TD-specific SE technique with a finer variance control and use a state-of-the-art stochastic LBFGS method. Lastly, we design an efficient TD multi-grid scheme together with preconditioning to speed up the convergence while avoiding local minima. All components are evaluated in extensive numerical proof-of-concept studies simulating a bowl-shaped 3D UST breast scanner prototype. Finally, we demonstrate that their combination allows us to obtain an accurate 442 × 442 × 222 voxel image with a resolution of 0.5 mm using Matlab on a single GPU within 24 h.

Parsimonious truncated Newton method for time-domain full-waveform inversion based on the Fourier-domain full-scattered-field approximation

To exploit Hessian information in full-waveform inversion (FWI), the matrix-free truncated Newton method can be used. In such a method, Hessian-vector product computation is one of the major concerns due to the huge memory requirements and demanding computational cost. Using the adjoint-state method, the Hessian-vector product can be estimated by zero-lag crosscorrelation of the first-/second-order incident wavefields and the second-/first-order adjoint wavefields. Different from the implementation in frequency-domain FWI, Hessian-vector product construction in the time domain becomes much more challenging because it is not affordable to store all of the time-dependent wavefields. The widely used wavefield recomputation strategy leads to computationally intensive tasks. We have developed an efficient alternative approach to computing the Hessian-vector product for time-domain FWI. In our method, discrete Fourier transform is applied to extract frequency-domain components of involved wavefields, which are used to compute wavefield crosscorrelation in the frequency domain. This makes it possible to avoid reconstructing the first- and second-order incident wavefields. In addition, a full-scattered-field approximation is proposed to efficiently simplify the second-order incident and adjoint wavefield computation, which enables us to refrain from repeatedly solving the first-order incident and adjoint equations for the second-order incident and adjoint wavefields (re)computation. With our method, the computational time can be reduced by 70% and 80% in viscous media for Gauss-Newton and full-Newton Hessian-vector product construction, respectively. The effectiveness of our method is also verified in the frame of a 2D multiparameter inversion, in which our method almost reaches the same iterative convergence of the conventional time-domain implementation.

Receiver-extension strategy for time-domain full-waveform inversion using a relocalization approach

A receiver-extension strategy is presented as an alternative to recently promoted source-extension strategies, in the framework of high-resolution seismic imaging by full-waveform inversion (FWI). This receiver-extension strategy is directly applicable in time-domain FWI, and, unlike source-extension methods, it incurs negligible extra computational cost. After connections between different source-extension strategies are reviewed, the receiver-extension method is introduced and analyzed for single-arrival data. The method results in a misfit function convex with respect to the velocity model in this context. The method is then applied to three exploration-scale synthetic case studies representative of different geologic environments, based on the Marmousi model, the BP 2004 salt model, and the Valhall model. In all three cases, the receiver-extension strategy makes it possible to start FWI with crude initial models and to reconstruct meaningful subsurface velocity models. The good performance of the method, even considering the inaccurate amplitude prediction due to noise, imperfect modeling, and source wavelet estimation, bodes well for field data applications.

Detection of karst voids at pile foundation by full-waveform inversion of single borehole sonic data

Drilled boreholes are commonly used to detect rock voids and fractures in underground construction projects such as pile foundations, tunnel borings, and other near surface human activities. However, voids can be missed by direct drilling due to their limited size. Borehole sonic sounding can be used to image voids in the vicinity of the borehole. We developed a new time domain full-waveform inversion method for borehole sounding data. The FWI method employ an amplitude envelope-based objective function which captures low-frequency components to provide the initial velocity model for further inversion. Then a global correlation-based objective function is applied for FWI based on the updated initial velocity model. Numerical studies show that modeled anomalies of a single void, multiple voids and fractures are successfully delineated by the inverted Vs velocity profile. The FWI method was applied to field data from a bridge foundation project in a karst area and compared to complementary cross-hole acoustic tomography borehole profiles and drilled corehole data. The single borehole FWI inverted velocity profiles identify voids in agreement with the cross hole tomography and corehole observations. The results demonstrate that FWI is capable of imaging karst voids near drilling boreholes using borehole sonic data.

Tangential cone condition and Lipschitz stability for the full waveform forward operator in the acoustic regime

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.