Frequency-domain Full Waveform Inversion Research Articles

Improving full-waveform inversion based on sparse regularization for geophysical data

Abstract Full-waveform inversion (FWI) is an advanced geophysical inversion technique. FWI provides images of subsurface structures with higher resolution in fields such as oil exploration and geology. The conventional algorithm minimizes the misfit error by calculating the least squares of the wavefield solutions between observed data and simulated data, followed by gradient direction and model update increment. Since the gradient is calculated by forward and backward wavefields, the high-accuracy model update relies on accurate forward and backward wavefield modelling. However, the quality of wavefield solutions obtained in practical situations could be poor and does not meet the requirements of high-resolution FWI. Specifically, the low-frequency wavefield is easily affected by noise and downsampling, which influences data quality, whereas the high-frequency wavefield is susceptible to spatial aliasing effects that produce imaging artefacts. Therefore, we propose using an algorithm called sparse relaxation regularized regression to optimize the wavefield solution in frequency-domain FWI, which is the forward and backward wavefield obtained from the Helmholtz equation, thus improving FWI's accuracy. The sparse relaxation regularized regression algorithm combines sparsity and regularization, allowing the broadband FWI to reduce the effects of noise and outliers, which can provide data supplementation in the low-frequency band and anti-aliasing in the high-frequency band. Our numerical examples demonstrate the wavefield optimization effect of the sparse relaxation regularized regression-based algorithm in various cases. The improved algorithm's accuracy and stability are verified compared to the Tikhonov regularization algorithm.

Combined non-convex second-order total variation with overlapping group sparsity for full waveform inversion

Full waveform inversion (FWI) can provide an accurate velocity model by matching observed and simulated seismograms. Mathematically, FWI is a highly ill-posed inverse problem that the observed data are independent of the model and cannot be inverted accurately. For a stable and reasonable inversion, proper regularization methods have to be taken into account. We propose a novel composite regularization for frequency-domain FWI problem, which uses a non-convex second-order total variation (TV) term and an overlapping group sparse TV (OGS-TV) regularization term. Compared with the conventional TV regularization, our method has better accuracy and robustness, and could effectively make use of the sparsity in the velocity model. Furthermore, the alternating direction multiplier algorithm with the adaptive selection of the penalty parameters is developed to solve this composite constraint problem, which can improve the stability of the FWI process. To illustrate the superior method both visually and quantitatively, we experimentally compare the proposed method with the conventional TV regularized FWI and the second-order total generalized variation regularized FWI and overlapping group sparsity regularized FWI on the well-known geological models.

3D frequency-domain acoustic full-waveform inversion using diffraction-angle filtering

The success of seismic imaging techniques such as full-waveform inversion (FWI) depends on using a good starting velocity model or a background velocity model, which can be effectively built by long-wavelength components in the early stages of FWI. This study proposes an optimal strategy for performing 3D frequency-domain FWI with diffraction-angle filtering (DAF). This approach can perform optimized single-frequency FWI with lower computational costs and build a long-wavelength velocity structure by separating low-wavenumber components with DAF. Prior to applying this strategy, we confirm the performance of 3D frequency-domain FWI with DAF in the radiation pattern analysis. DAF effectively separates the diffraction energy of partial derivative wavefields according to the diffraction angle with a negligible increase in computational costs. We then apply our 3D FWI strategy with DAF on synthetic data generated from the 3D SEG/EAGE overthrust model. The wavenumber components are effectively separated, and a long-wavelength velocity model is estimated with DAF compared to conventional FWI. Finally, we apply our strategy to 3D ocean-bottom seismic data from the North Sea. As a result, we observe that 3D FWI with the DAF method successfully reconstructs a long-wavelength velocity structure on the real data using only a single-frequency component.

Inner product preconditioned trust-region methods for frequency-domain full waveform inversion

Full waveform inversion is a seismic imaging method which requires solving a large-scale minimization problem, typically through local optimization techniques. Most local optimization methods can basically be built up from two choices: the update direction and the strategy to control its length. In the context of full waveform inversion, this strategy is very often a line search. We here propose to use instead a trust-region method, in combination with non-standard inner products which act as preconditioners. More specifically, a line search and several trust-region variants of the steepest descent, the limited memory BFGS algorithm and the inexact Newton method are presented and compared. A strong emphasis is given to the inner product choice. For example, its link with preconditioning the update direction and its implication in the trust-region constraint are highlighted. A first numerical test is performed on a 2D synthetic model then a second configuration, containing two close reflectors, is studied. The latter configuration is known to be challenging because of multiple reflections. Based on these two case studies, the importance of an appropriate inner product choice is highlighted and the best trust-region method is selected and compared to the line search method. In particular we were able to demonstrate that using an appropriate inner product greatly improves the convergence of all the presented methods and that inexact Newton methods should be combined with trust-region methods to increase their convergence speed.

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.

Frequency domain full waveform inversion with small‐scale laboratory measurements for reconnaissance in mechanized tunneling

AbstractIn mechanized tunneling, so‐called tunnel boring machines (TBMs) drill through the ground in an automatized manner. Therefore, the drilling process is very efficient but the maintenance costs are high. Seismic exploration seems appropriate to identify changing ground conditions in front of the TBM to reduce costs from damages of the TBM and from the corresponding dwell times. Today's seismic exploration techniques are using only a small amount of the information, which is contained in the seismic records from field observations, whereas full waveform inversion (FWI) tries to use the whole content. The potential of different FWI approaches for the application in mechanized tunneling has been investigated. An analysis of the performance of FWI not only employing synthetic examples but additionally measured waveforms is essential. Since seismic surveys at the construction side are not performed with FWI in mind, the data sets are usually not appropriate for testing the developed algorithms. Nevertheless, a validation of FWI approaches by measured waveforms is possible by using waveform recordings from a small‐scale experimental setup. A small‐scale super high strength grout specimen is constructed for validating an adjoint frequency domain FWI approach. Frequency domain models compute the seismic response of a system for an infinite time interval. The attenuation of the material as well as the attenuation effects at the free surfaces are hard to quantify over an infinite time interval. Therefore, the specimen is designed in a way that four of the six borders can be modeled as absorbing boundaries. For this purpose, the displacement recordings are truncated just before the waves that are reflected at the excluded borders arrive at the measurement points. This design in combination with a notch representing a rectangular tunnel makes the measurements more similar to data from a tunnel construction side. A rectangular hole is embedded in the specimen and acts as material discontinuity, which the FWI approach is aiming to detect.

Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? The Gorgon data FWI case study

Frequency-domain full-waveform inversion (FWI) is potentially amenable to efficient processing of full-azimuth long-offset stationary-recording seabed acquisition carried out with a sparse layout of ocean-bottom nodes (OBNs) and broadband sources because the inversion can be performed with a few discrete frequencies. However, computing the solution of the forward (boundary-value) problem efficiently in the frequency domain with linear algebra solvers remains a challenge for large computational domains involving tens to hundreds of millions of parameters. We illustrate the feasibility of 3D frequency-domain FWI with a subset of the 2015/2016 Gorgon OBN data set in the North West Shelf, Australia. We solve the forward problem with the massively parallel multifrontal direct solver MUMPS, which includes four key features to reach high computational efficiency: an efficient parallelism combining message-passing interface and multithreading, block low-rank compression, mixed-precision arithmetic, and efficient processing of sparse sources. The Gorgon subdata set involves 650 OBNs that are processed as reciprocal sources and 400,000 sources. Monoparameter FWI for vertical wavespeed is performed in the viscoacoustic vertically transverse isotropic approximation with a classical frequency continuation approach proceeding from a starting frequency of 1.7 Hz to a final frequency of 13 Hz. The target covers an area ranging from 260 km2 (frequency ≥ 8.5 Hz) to 705 km2 (frequency ≤ 8.5 Hz) for a maximum depth of 8 km. Compared to the starting model, FWI dramatically improves the reconstruction of the bounding faults of the Gorgon horst at reservoir depths as well as several intrahorst faults and several horizons of the Mungaroo Formation down to a depth of 7 km. Seismic modeling reveals a good kinematic agreement between recorded and simulated data, but amplitude mismatches between the recorded and simulated reflection from the reservoir suggest elastic effects. Therefore, future works involve multiparameter reconstruction for density and attenuation before considering elastic FWI from hydrophone and geophone data.

Frequency-Domain Full-Waveform Inversion Based on Tunnel Space Seismic Data

Tunnel seismic detection methods are effective for obtaining the geological structure around the tunnel face, which is critical for safe construction and disaster mitigation in tunnel engineering. However, there is often a lack of accuracy in the acquired geological information and physical properties ahead of the tunnel face in the current tunnel seismic detection methods. Thus, we apply a frequency-domain acoustic full-waveform inversion (FWI) method to obtain high-resolution results for the tunnel structure. We discuss the influence of the frequency group selection strategy and the tunnel observation system settings regarding the inversion results and determine the structural imaging and physical property parameter inversion of abnormal geological bodies ahead of the tunnel face. Based on the conventional strategies of frequency-domain acoustic FWI, we propose a frequency group selection strategy that combines a low-frequency selection covering the vertical wavenumber and a high-frequency selection of anti-aliasing. This strategy can effectively obtain the spatial structure and physical parameters of the geology ahead of the tunnel face and improve the inversion resolution. In addition, by linearly increasing the side length of the tunnel observation system, we share the influence of the length of the two sides of the observation systems of different tunnels on the inversion results. We determined that the inversion results are best when the side length is approximately five times the width of the tunnel face, and the influence of increasing the side observation length beyond this range on the inversion results can be ignored. Finally, based on this approach, we invert for the complex multi-stratum model, and an accurate structure and physical property parameters of the complex stratum ahead of the tunnel face are obtained, which verifies the feasibility of the proposed method.

Three-dimensional finite-difference finite-element frequency-domain wave simulation with multi-level optimized additive Schwarz domain-decomposition preconditioner: A tool for FWI of sparse node datasets

Efficient frequency-domain full-waveform inversion (FWI) of long-offset node data can be designed with a few discrete frequencies, which lead to modest data volumes to be managed during the inversion process. Moreover, attenuation effects can be straightforwardly implemented in the forward problem without the computational overhead. However, 3D frequency-domain seismic modeling is challenging because it requires solving a large and sparse linear indefinite system for each frequency with multiple right-hand sides (RHSs). This linear system can be solved by direct or iterative methods. The former allows efficient processing of multiple RHSs but may suffer from limited scalability for very large problems. Iterative methods equipped with a domain-decomposition preconditioner provide a suitable alternative to process large computational domains for sparse-node acquisition. We have investigated the domain-decomposition preconditioner based on the optimized restricted additive Schwarz (ORAS) method, in which a Robin or perfectly matched layer condition is implemented at the boundaries between the subdomains. The preconditioned system is solved by a Krylov subspace method, whereas a block low-rank lower-upper decomposition of the local matrices is performed at a preprocessing stage. Multiple sources are processed in groups with a pseudoblock method. The accuracy, the computational cost, and the scalability of the ORAS solver are assessed against several realistic benchmarks. In terms of discretization, we compare a compact wavelength-adaptive 27-point finite-difference stencil on a regular Cartesian grid with a P3 finite-element method on h-adaptive tetrahedral mesh. Although both schemes have comparable accuracy, the former is more computationally efficient, the latter being beneficial to comply with known boundaries such as bathymetry. The scalability of the method, the block processing of multiple RHSs, and the straightforward implementation of attenuation, which further improves the convergence of the iterative solver, make the method a versatile forward engine for large-scale 3D FWI applications from sparse node data sets.

Accurate 3D frequency-domain seismic wave modeling with the wavelength-adaptive 27-point finite-difference stencil: A tool for full-waveform inversion

Efficient frequency-domain full-waveform inversion (FWI) of long-offset node data can be performed with a few frequencies. The seismic response of these frequencies can be computed with compact finite-difference stencils on regular Cartesian grid with direct or hybrid direct/iterative methods. Compactness, which is necessary to mitigate the fill-in induced by the lower-upper factorization, is implemented with second-order stencils, whereas accuracy is achieved by building a consistent mass matrix and a compound stiffness matrix with optimal weights so that the stencil covers the eight cells surrounding the central point, leading to a 27-point stencil. Classical approaches estimate constant weights by jointly minimizing numerical dispersion in homogeneous media for several numbers of grid points per wavelength ([Formula: see text]). Then, the impedance matrix is built by using these constant weights at each central point covering the heterogeneous subsurface model, leading to nonuniform wavefield accuracy. Instead, we estimate [Formula: see text]-dependent weights once and for all by minimizing dispersion separately for each value of [Formula: see text] in a range found in FWI applications. Then, we build the impedance matrix without computational overhead by selecting at each central point the weights corresponding to the local [Formula: see text], hence leading to a wavelength-adaptive stencil. This separate approach with adaptive weights makes wavefield accuracy uniform. We benchmark wavefields, which are computed in several 3D large-scale subsurface models with a sparse multifrontal direct solver and the nonadaptive/adaptive stencils, against analytical solutions when available and the highly accurate discretization-free convergent Born series method. Each benchmark reveals the higher accuracy of the adaptive stencil relative to the nonadaptive one. In the presence of sharp contrasts, the adaptive stencil also is more accurate than a classical [Formula: see text] finite-different time-domain staggered-grid stencil. The relevance of the adaptive stencil is finally illustrated with a 3D FWI case study in the 3.5–13 Hz frequency band.

Frequency Domain Full Waveform Inversion Based on Generalized FD Method

The implementation of an accurate and efficient forward/backward modelling approach plays a very important role in full waveform inversion (FWI), in which the finite difference method (FDM) with rectangular grid is mostly adopted. However, the FDM has the ripple-shaped divergence/scattering problem when dealing with uneven interface. On the contrary, a meshless method called GFDM is suitable for flexibility of simulation region, and can be of interest for solving the above problem. Hence, we develop a meshless frequency-domain acoustic modeling method with GFDM to improve the accuracy and eliminate the divergence. Then we apply the GFDM to inverse problem and propose a meshless FWI method. The proposed method reduces lots of the computations and memory usage of frequency-domain FWI, and successfully suppress the ripple-shaped scattering. Numerical examples demonstrate that our method can efficiently improve numerical accuracy and recover an accurate velocity model, which is helpful for high-precision seismic exploration through the complex earth.

Randomized Source Sketching for Full Waveform Inversion

Partial differential equation (PDE) constrained optimization problems such as seismic full waveform inversion (FWI) frequently arise in the geoscience and related fields. For such problems, many observations are usually gathered by multiple sources, which form the right-hand-sides of the PDE constraint. Solving the inverse problem with such massive data sets is computationally demanding, in particular when dealing with large number of model parameters. This paper proposes a novel randomized source sketching method for the efficient resolution of multisource PDE constrained optimization problems. We first formulate the different source-encoding strategies used in seismic imaging into a unified framework based on a randomized sketching. To this end, the source dimension of the problem is projected in a smaller domain by a suitably defined projection matrix that gathers the physical sources in super-sources through a weighted summation. This reduction in the number of physical sources decreases significantly the number of PDE solves while suitable sparsity-promoting regularization can efficiently mitigate the footprint of the cross-talk noise to maintain the convergence speed of the algorithm. We implement the randomized sketching method in an extended search-space formulation of frequency-domain FWI, which relies on the alternating direction method of multipliers (ADMM). Numerical examples carried out with a series of well documented 2D benchmarks demonstrate that the randomized sketching algorithm reduces the cost of large-scale problems by at least one order of magnitude compared to the original deterministic algorithm.

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.

Frequency Domain Full Waveform Inversion Method of Acquiring Rock Wave Velocity in Front of Tunnels

Ahead geological prospecting, which can estimate adverse geology ahead of the tunnel face, is necessary in the process of tunnel construction. Due to its long detection range and good recognition effect on the interface, the seismic method is widely used in tunnel ahead prospecting. However, the observation space in tunnels is quite narrow compared to ground seismic prospecting, which leads to some problems in the acquisition of wave velocity, including: the velocity of the direct wave is used to replace the wave velocity of the forward rock approximately; the arrival time information of seismic waves is the main factor in time-travel inversion or the tomography method, which is sufficient to provide a simple model rather than deal with complex geological conditions. In view of the above problems, the frequency domain full waveform inversion method in ground prospecting is introduced to tunnel seismic prospecting. In addition, the optimized difference format is given according to the particularity of the tunnel environment. In this method, the kinematics and dynamics of the seismic wavefield are fully used to obtain more accurate wave velocity results. Simultaneously, forward modeling and inversion simulations on tunnel samples with typical adverse geological bodies are given here, which verified the validity and reliability of the proposed method.

Quantitative Inspection of Complex-Shaped Parts Based on Ice-Coupled Ultrasonic Full Waveform Inversion Technology

Ultrasonic methods have been extensively developed in nondestructive testing for various materials and components. However, accurately extracting quantitative information about defects still remains challenging, especially for complex structures. Although the immersion technique is commonly used for complex-shaped parts, the large mismatch of acoustic impedance between water and metal prevents effective ultrasonic transmission and leads to a low signal-to-noise ratio(SNR). In this paper, a quantitative imaging method is proposed for complex-shaped parts based on an ice-coupled full waveform inversion (FWI) method. Numerical experiments were carried out to quantitatively inspect the various defects in a turbine blade. Firstly, the k-space pseudospectral method was applied to simulate ice-coupled ultrasonic testing for the turbine blade. The recorded full wavefields were then applied for a frequency-domain FWI based on the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method. With a carefully selected iterative number and frequency, a successive-frequency FWI can well detect half wavelength defects. Extended studies on an open notch with different orientations and multiple adjacent defects proved its capability to detect different types of defects. Finally, an uncertainty analysis was conducted with inaccurate initial velocity models with a relative error of ±2%, demonstrating its robustness even with a certain inaccuracy. This study demonstrates that the proposed method has a high potential to inspect complex-shaped structures with an excellent resolution.

On efficient frequency-domain full-waveform inversion with extended search space

Efficient frequency-domain full-waveform inversion (FD-FWI) of wide-aperture data is designed by limiting inversion to a few frequencies and by solving the Helmholtz equation with a direct solver to process multiple sources efficiently. Some variants of FD-FWI, which process the wave equation as a weak constraint, have been proposed to increase the computational efficiency or extend the search space. Among them, the contrast-source reconstruction inversion (CSRI) reparameterizes FD-FWI in terms of contrast sources (CS) and contrasts and updates them in an alternating mode. This reparameterization allows for one lower-upper (LU) decomposition of the Helmholtz operator to be performed per frequency inversion hence further improving the computational efficiency of FD-FWI. However, iteratively refined wavefield reconstruction inversion (IR-WRI) relies on the alternating-direction method of multipliers to extend the search space by matching the data from the early iterations via an aggressive relaxation of the wave equation while satisfying it at the convergence point thanks to the defect correction performed by the Lagrange multipliers. In contrast to CSRI, IR-WRI requires redoing one LU decomposition when the subsurface model is updated. In both methods, the CSs or the wavefields are computed by solving in a least-squares sense an overdetermined linear system gathering an observation equation and a wave equation. A drawback of CSRI is that CSs are estimated approximately with one iteration of a conjugate gradient method, whereas the wavefields are reconstructed exactly by IR-WRI with a Gauss-Newton method. We have combined the benefits of CSRI and IR-WRI to decrease the number of LU decomposition during IR-WRI with a fixed-point (FP) algorithm while preserving its search space extension capability. Application on the 2D complex Marmousi and the BP salt models shows that our FP-based IR-WRI manages to reconstruct these models as accurately as the classic IR-WRI while reducing the number of LU factorizations considerably. A theoretical complexity analysis and a recent application of 3D FD-FWI based upon direct solver suggest that the FP algorithm should reduce the cost of IR-WRI by a factor of approximately 2 and 10 for 3D dense ocean bottom cable and sparse ocean bottom node acquisitions, respectively.

Three-dimensional frequency-domain full waveform inversion based on the nearly-analytic discrete method

The nearly analytic discrete (NAD) method is a kind of finite difference method with advantages of high accuracy and stability. Previous studies have investigated the NAD method for simulating wave propagation in the time-domain. This study applies the NAD method to solving three-dimensional (3D) acoustic wave equations in the frequency-domain. This forward modeling approach is then used as the “engine” for implementing 3D frequency-domain full waveform inversion (FWI). In the numerical modeling experiments, synthetic examples are first given to show the superiority of the NAD method in forward modeling compared with traditional finite difference methods. Synthetic 3D frequency-domain FWI experiments are then carried out to examine the effectiveness of the proposed methods. The inversion results show that the NAD method is more suitable than traditional methods, in terms of computational cost and stability, for 3D frequency-domain FWI, and represents an effective approach for inversion of subsurface model structures.

Frequency domain full-waveform inversion in a fluid-saturated poroelastic medium

SUMMARYDetermining subsurface properties is of fundamental importance in exploration seismic imaging. Poroelasticity theory provides an opportunity to extract quantitative fluid- and attenuation-dependent properties from seismic data. Following Pratt’s frequency-domain full-waveform-inversion (FWI) procedure and extending the basic FWI equations from the elastic case to the poroelastic case, we implement poroelastic FWI (PFWI) of fluid-saturated porous media. By analysing the sensitivity kernels of poroelastic parameters, we explain the reason why some parameters are more difficult to recover than others. We also show the analytical and numerical radiation patterns based on which we predict the trade-offs among parameters. In numerical experiments, we invert two models to demonstrate the feasibility and effectiveness of the proposed PFWI and to verify our predictions about trade-offs for two-parameter PFWI. Finally, we discuss the various factors that influence the inversion results.

A gradient-based Markov chain Monte Carlo method for full-waveform inversion and uncertainty analysis

Traditional full-waveform inversion (FWI) methods only render a “best-fit” model that cannot account for uncertainties of the ill-posed inverse problem. Additionally, local optimization-based FWI methods cannot always converge to a geologically meaningful solution unless the inversion starts with an accurate background model. We seek the solution for FWI in the Bayesian inference framework to address those two issues. In Bayesian inference, the model space is directly probed by sampling methods such that we obtain a reliable uncertainty appraisal, determine optimal models, and avoid entrapment in a small local region of the model space. The solution of such a statistical inverse method is completely described by the posterior distribution, which quantifies the distributions for parameters and inversion uncertainties. To efficiently sample the posterior distribution, we introduce a sampling algorithm in which the proposal distribution is constructed by the local gradient and the diagonal approximate Hessian of the local log posterior. Our algorithm is called the gradient-based Markov chain Monte Carlo (GMCMC) method. The GMCMC FWI method can quantify inversion uncertainties with estimated posterior distribution given sufficiently long Markov chains. By directly sampling the posterior distribution, we obtain a global view of the model space. Theoretically speaking, statistical assessments do not depend on starting models. Our method is applied to the 2D Marmousi model with the frequency-domain FWI setting. Numerical results suggest that our method can be readily applied to 2D cases with affordable computational efforts.

High-resolution bone microstructure imaging based on ultrasonic frequency-domain full-waveform inversion* *Project supported by the National Natural Science Foundation of China (Grant Nos. 11827808, 11874289, and 11804056), the National Science Fund for Distinguished Young Scholars of China (Grant No. 11525416), Shanghai Municipal Science and Technology Major Project, China (Grant No. 2017SHZDZX01), Shanghai Talent Development Fund (Grant No. 2018112), State Key Laboratory

The main challenge in bone ultrasound imaging is the large acoustic impedance contrast and sound velocity differences between the bone and surrounding soft tissue. It is difficult for conventional pulse-echo modalities to give accurate ultrasound images for irregular bone boundaries and microstructures using uniform sound velocity assumption rather than getting a prior knowledge of sound speed. To overcome these limitations, this paper proposed a frequency-domain full-waveform inversion (FDFWI) algorithm for bone quantitative imaging utilizing ultrasonic computed tomography (USCT). The forward model was calculated in the frequency domain by solving the full-wave equation. The inverse problem was solved iteratively from low to high discrete frequency components via minimizing a cost function between the modeled and measured data. A quasi-Newton method called the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS) was utilized in the optimization process. Then, bone images were obtained based on the estimation of the velocity and density. The performance of the proposed method was verified by numerical examples, from tubular bone phantom to single distal fibula model, and finally with a distal tibia-fibula pair model. Compared with the high-resolution peripheral quantitative computed tomography (HR-pQCT), the proposed FDFWI can also clearly and accurately presented the wavelength scaled pores and trabeculae in bone images. The results proved that the FDFWI is capable of reconstructing high-resolution ultrasound bone images with sub-millimeter resolution. The parametric bone images may have the potential for the diagnosis of bone disease.