Search

Concerns raised by Okazaki & Ueda (2022) on the paper by Sambridge et al. (2022) are addressed. Two issues are discussed and some new numerical results presented. The first concerns whether the properties of the Wasserstein time-series misfit introduced in our earlier paper will translate to model space non-uniqueness in a seismic waveform inversion setting. It is argued that this is unlikely, given the special conditions, which must exist between all observed/predicted seismic waveform pairs for non-uniqueness to result. The second issue discussed is the efficacy of using the Sliced Wasserstein algorithm of Bonneel et al. (2015) as an alternate to the marginal Wasserstein algorithm, as proposed by Okazaki & Ueda (2022). It is argued that for optimization-based waveform fitting, the Sliced Wasserstein algorithm is a viable alternate provided care is taken to ensure that conditions arise which do invalidate analytical derivative expressions of the resulting Wasserstein misfit. In practice, this would likely mean recasting the 2D Optimal Transport problem posed in our earlier paper onto unstructured grids. [ABSTRACT FROM AUTHOR]

We propose a new approach to measuring the agreement between two oscillatory time-series, such as seismic waveforms, and demonstrate that it can be used effectively in inverse problems. Our approach is based on Optimal Transport theory and the Wasserstein distance, with a novel transformation of the time-series to ensure that necessary normalization and positivity conditions are met. Our measure is differentiable, and can readily be used within an optimization framework. We demonstrate performance with a variety of synthetic examples, including seismic source inversion, and observe substantially better convergence properties than achieved with conventional L 2 misfits. We also briefly discuss the relationship between Optimal Transport and Bayesian inference. [ABSTRACT FROM AUTHOR]

Monte Carlo methods are widespread in geophysics and have proved to be powerful in non-linear inverse problems. However, they are associated with significant practical challenges, including long calculation times, large output ensembles of Earth models, and difficulties in the appraisal of the results. This paper addresses some of these challenges using generative models, a family of tools that have recently attracted much attention in the machine learning literature. Generative models can, in principle, learn a probability distribution from a set of given samples and also provide a means for rapid generation of new samples which follow that approximated distribution. These two features make them well suited for application to the outputs of Monte Carlo algorithms. In particular, training a generative model on the posterior distribution of a Bayesian inference problem provides two main possibilities. First, the number of parameters in the generative model is much smaller than the number of values stored in the ensemble, leading to large compression rates. Secondly, once trained, the generative model can be used to draw any number of samples, thereby eliminating the dependence on an often large and unwieldy ensemble. These advantages pave new pathways for the use of Monte Carlo ensembles, including improved storage and communication of the results, enhanced calculation of numerical integrals, and the potential for convergence assessment of the Monte Carlo procedure. Here, these concepts are initially demonstrated using a simple synthetic example that scales into higher dimensions. They are then applied to a large ensemble of shear wave velocity models of the core–mantle boundary, recently produced in a Monte Carlo study. These examples demonstrate the effectiveness of using generative models to approximate posterior ensembles, and indicate directions to address various challenges in Monte Carlo inversion. [ABSTRACT FROM AUTHOR]

We develop a theoretical framework for framing and solving probabilistic linear(ized) inverse problems in function spaces. This is built on the statistical theory of Gaussian Processes, and allows results to be obtained independent of any basis, avoiding any difficulties associated with the fidelity of representation that can be achieved. We show that the results of Backus–Gilbert theory can be fully understood within our framework, although there is not an exact equivalence due to fundamental differences of philosophy between the two approaches. Nevertheless, our work can be seen to unify several strands of linear inverse theory, and connects it to a large body of work in machine learning. We illustrate the application of our theory using a simple example, involving determination of Earth's radial density structure. [ABSTRACT FROM AUTHOR]

By starting from a general framework for probabilistic continuous inversion (developed in Part I) and introducing discrete basis functions, we obtain the well-known algorithms for probabilistic least-squares inversion set out by Tarantola & Valette. In doing so, we establish a direct equivalence between the spatial covariance function that must be specified in continuous inversion, and the combination of basis functions and prior covariance matrix that must be chosen for discretized inversion. We show that the common choice of Tikhonov regularization (⁠|$\mathbf {C_m^{-1}} = \sigma ^2\mathbf {I}$|⁠) arises from a delta-function spatial covariance, and that this lies behind many of the artefacts commonly associated with discretized inversion. We show that other choices of spatial covariance function can be used to generate regularization matrices yielding substantially better results, and permitting localization of features even if global basis functions are used. We are also able to offer a straightforward explanation for the spectral leakage problem identified by Trampert & Snieder. [ABSTRACT FROM AUTHOR]

Obtaining slip distributions for earthquakes results in an ill-posed inverse problem. While this implies that only limited and uncertain information can be recovered from the data, inferences are typically made based only on a single regularized model. Here, we develop an inversion approach that can quantify uncertainties in a Bayesian probabilistic framework for the finite fault inversion (FFI) problem. The approach is suitably efficient for rapid source characterization and includes positivity constraints for model parameters, a common practice in FFI, via coordinate transformation to logarithmic space. The resulting inverse problem is nonlinear and the most probable solution can be obtained by iterative linearization. In addition, model uncertainties are quantified by approximating the posterior probability distribution by a Gaussian distribution in logarithmic space. This procedure is straightforward since an analytic expression for the Hessian of the objective function is obtained. In addition to positivity, we apply smoothness regularization to the model in logarithmic space. Simulations based on surface wave data show that smoothing in logarithmic space penalizes abrupt slip changes less than smoothing in linear space. Even so, the main slip features of models that are smooth in linear space are recovered well with logarithmic smoothing. Our synthetic experiments also show that, for the data set we consider, uncertainty is low at the shallow portion of the fault and increases with depth. In addition, a simulation with a large station azimuthal gap of 180° significantly increases the slip uncertainties. Further, the marginal posterior probabilities obtained from our approximate method are compared with numerical Markov Chain Monte Carlo sampling. We conclude that the Gaussian approximation is reasonable and meaningful inferences can be obtained from it. Finally, we apply the new approach to observed surface wave records from the great Illapel earthquake (Chile, 2015, M w = 8.3). The location and amplitude of our inferred peak slip is consistent with other published solutions but the spatial slip distribution is more compact, likely because of the logarithmic regularization. We also find a minor slip patch downdip, mainly in an oblique direction, which is poorly resolved compared to the main slip patch and may be an artefact. We conclude that quantifying uncertainties of finite slip models is crucial for their meaningful interpretation, and therefore rapid uncertainty quantification can be critical if such models are to be used for emergency response. [ABSTRACT FROM AUTHOR]

Most linear inverse problems require regularization to ensure that robust and meaningful solutions can be found. Typically, Tikhonov-style regularization is used, whereby a preference is expressed for models that are somehow 'small' and/or 'smooth'. The strength of such preferences is expressed through one or more damping parameters, which control the character of the solution, and which must be set by the user. However, identifying appropriate values is often regarded as a matter of art, guided by various heuristics. As a result, such choices have often been the source of controversy and concern. By treating these as hyperparameters within a hierarchical Bayesian framework, we are able to obtain solutions that encompass the range of permissible regularization parameters. Furthermore, we show that these solutions are often well-approximated by those obtained via standard analysis using certain regularization choices which are—in a certain sense—optimal. We obtain algorithms for determining these optimal values in various cases of common interest, and show that they generate solutions with a number of attractive properties. A reference implementation of these algorithms, written in Python, accompanies this paper. [ABSTRACT FROM AUTHOR]

In geophysical inversion, inferences of Earth's properties from sparse data involve a trade-off between model complexity and the spatial resolving power. A recent Markov chain Monte Carlo (McMC) technique formalized by Green, the so-called trans-dimensional samplers, allows us to sample between these trade-offs and to parsimoniously arbitrate between the varying complexity of candidate models. Here we present a novel framework using transdimensional sampling over tree structures. This new class of McMC sampler can be applied to 1-D, 2-D and 3-D Cartesian and spherical geometries. In addition, the basis functions used by the algorithm are flexible and can include more advanced parametrizations such as wavelets, both in Cartesian and Spherical geometries, to permit Bayesian multiscale analysis. This new framework offers greater flexibility, performance and efficiency for geophysical imaging problems than previous sampling algorithms. Thereby increasing the range of applications and in particular allowing extension to trans-dimensional imaging in 3-D. Examples are presented of its application to 2-D seismic and 3-D teleseismic tomography including estimation of uncertainty. [ABSTRACT FROM AUTHOR]

This paper develops a probabilistic Bayesian approach to the problem of inferring the spatiotemporal evolution of earthquake rupture on a fault surface from seismic data with rigorous uncertainty estimation. To date, uncertainties of rupture parameters are poorly understood, and the effect of choices such as fault discretization on uncertainties has not been studied. We show that model choice is fundamentally linked to uncertainty estimation and can have profound effects on results. The approach developed here is based on a trans-dimensional self-parametrization of the fault, avoids regularization constraints and provides rigorous uncertainty estimation that accounts for model-selection ambiguity associated with the fault discretization. In particular, the fault is parametrized using self-adapting irregular grids which intrinsically match the local resolving power of the data and provide parsimonious solutions requiring few parameters to capture complex rupture characteristics. Rupture causality is ensured by parametrizing rupture-onset time by a rupture-velocity field and obtaining first rupture times from the eikonal equation. The Bayesian sampling of the parameter space is implemented on a computer cluster with a highly efficient parallel tempering algorithm.The inversion is applied to simulated and observed W-phase waveforms from the 2010 Maule (Chile) earthquake. Simulation results show that our approach avoids both over- and underparametrization to ensure unbiased inversion results with uncertainty estimates that are consistent with data information. The simulation results also show the ability of W-phase data to resolve the spatial variability of slip magnitude and rake angles. In addition, sensitivity to spatially dependent rupture velocities exists for kinematic slip models.Application to the observed data indicates that residual errors are highly correlated and likely dominated by theory error, necessitating the iterative estimation of a non-stationary data covariance matrix. The moment magnitude for the Maule earthquake is estimated to be ∼8.9, with slip concentrated in two zones updip of and north and south of the hypocentre, respectively. While this aspect of the slip distribution is similar to previous studies, we show that the slip maximum in the southern zone is poorly resolved compared to the northern zone. Both slip maxima are higher than reported in previous studies, which we speculate may be due to the lack of bias caused by the regularization used in other studies. [ABSTRACT FROM AUTHOR]

Non-linear inverse problems in the geosciences often involve probabilistic sampling of multimodal density functions or global optimization and sometimes both. Efficient algorithmic tools for carrying out sampling or optimization in challenging cases are of major interest. Here results are presented of some numerical experiments with a technique, known as Parallel Tempering, which originated in the field of computational statistics but is finding increasing numbers of applications in fields ranging from Chemical Physics to Astronomy. To date, experience in use of Parallel Tempering within earth sciences problems is very limited. In this paper, we describe Parallel Tempering and compare it to related methods of Simulated Annealing and Simulated Tempering for optimization and sampling, respectively. A key feature of Parallel Tempering is that it satisfies the detailed balance condition required for convergence of Markov chain Monte Carlo (McMC) algorithms while improving the efficiency of probabilistic sampling. Numerical results are presented on use of Parallel Tempering for trans-dimensional inversion of synthetic seismic receiver functions and also the simultaneous fitting of multiple receiver functions using global optimization. These suggest that its use can significantly accelerate sampling algorithms and improve exploration of parameter space in optimization. Parallel Tempering is a meta-algorithm which may be used together with many existing McMC sampling and direct search optimization techniques. It's generality and demonstrated performance suggests that there is significant potential for applications to both sampling and optimization problems in the geosciences. [ABSTRACT FROM AUTHOR]

SUMMARY A meaningful interpretation of seismic measurements requires a rigorous quantification of the uncertainty. In an inverse problem, the data noise determines how accurately observations should be fit, and ultimately the level of detail contained in the recovered model. A common problem in seismic tomography is the difficulty in quantifying data uncertainties, and thus the required level of data fit. Traditionally, the complexity of the solution model (defined by both the number of basis functions and the regularization) is defined arbitrarily by the user prior to inversion with only limited use of data errors. In the context of multiscale problems, dealing with multiple data sets that are characterized by different noise variances and that span the Earth at different scales is a major challenge. Practitioners are usually required to arbitrarily weigh the contribution of each data type into the final solution. Furthermore, the basis functions are usually spatially uniform across the velocity field and regularization procedures are global, which prevents the solution model from accounting for the uneven spatial distribution of information. In this work we propose to address these issues with a Hierarchical Bayesian inversion. The new algorithm represents an extension of the transdimensional tomography to account for uncertainties in data noise. This approach has the advantage of treating the level of noise in each data set, as well as the number of model parameters, as unknowns in the inversion. It provides a parsimonious solution that fully represents the degree of knowledge one has about seismic structure (i.e. constraints, resolution and trade-offs). Rather than being forced to make decisions on parametrization, level of data fit and weights between data types in advance, as is often the case in an optimization framework, these choices are relaxed and instead constrained by the data themselves. The new methodology is presented in a synthetic example where both the data density and the underlying structure contain multiple length scales. Three ambient noise data sets that span the Australian continent at different scales are then simultaneously inverted to infer a multiscale tomographic image of Rayleigh wave group velocity for the Australian continent. The procedure turns out to be particularly useful when dealing with multiple data types with different unknown levels of noise as the algorithm is able to naturally adjust the fit to the different data sets and provide a velocity map with a spatial resolution adapted to the spatially variable information present in the data. [ABSTRACT FROM AUTHOR]

There is significant seismic activity in the region around Australia, largely due to the plate boundaries to the north and to the east of the mainland. This activity results in serious seismic and tsunami hazard in the coastal areas of Australia. Hence seismicity is and will be monitored in real time by Geoscience Australia (GA), which uses a network of permanent broadband seismometers. Seismic moment tensor (MT) solutions are currently determined using 1-D, radially symmetric models of Earth and this requires augmentation by recording stations located outside of Australia. A 3-D model of the Australian continent developed recently using full waveform tomography now offers the opportunity to significantly improve the determination of MT solutions of earthquakes from tectonically active regions. A complete-waveform, time-domain MT inversion method has been developed using a point-source approximation. A series of synthetic tests using first a 1-D and then a 3-D structural model has been performed. The feasibility of deploying 3-D versus 1-D Earth structure for the inversion of seismic data has been studied and the advantages of using the 3-D structural model were illustrated with examples. The 3-D model is superior to the 1-D model, as a number of sensitivity tests show. The ultimate goal of this work is an automated MT inversion system in Australia relying on GA and other international stations, although more work remains to be done before the full implementation of such a scheme in real time. [ABSTRACT FROM AUTHOR]

We present a globally distributed data set of ∼400 000 frequency-dependent SH-wave traveltimes. An automated technique is used to measure teleseismic S, ScS and SS traveltimes at several periods ranging from 10 to 51 s. The targeted seismic phases are first extracted from the observed and synthetic seismograms using an automated time window algorithm. Traveltimes are then measured at several periods, by cross-correlation between the selected observed and synthetic filtered waveforms. Frequency-dependent effects due to crustal reverberations beneath each receiver are handled by incorporating crustal phases into WKBJ synthetic waveforms. After correction for physical dispersion due to intrinsic anelastic processes, we observe a residual traveltime dispersion on the order of 1–2 s in the period range of analysis. This dispersion occurs differently for S, ScS and SS, which is presumably related to their differing paths through the Earth. We find that: (1) Wavefront-healing phenomenon is observed for S and to a lesser extent SS waves having passed through very low velocity anomalies. (2) A preferred sampling of high velocity scatterers located at the CMB may explain our observation that ScS waves travel faster at low-frequency than at high-frequency. (3) A frequency-dependent attenuation , with , is compatible with the globally averaged dispersion observed for S waves. [ABSTRACT FROM AUTHOR]

The reversible jump algorithm is a statistical method for Bayesian inference with a variable number of unknowns. Here, we apply this method to the seismic tomography problem. The approach lets us consider the issue of model parametrization (i.e. the way of discretizing the velocity field) as part of the inversion process. The model is parametrized using Voronoi cells with mobile geometry and number. The size, position and shape of the cells defining the velocity model are directly determined by the data. The inverse problem is tackled within a Bayesian framework and explicit regularization of model parameters is not required. The mobile position and number of cells means that global damping procedures, controlled by an optimal regularization parameter, are avoided. Many velocity models with variable numbers of cells are generated via a transdimensional Markov chain and information is extracted from the ensemble as a whole. As an aid to interpretation we visualize the expected earth model that is obtained via Monte Carlo integration in a straightforward manner. The procedure is particularly adept at imaging rapid changes or discontinuities in wave speed. While each velocity model in the final ensemble consists of many discontinuities at cell boundaries, these are smoothed out in the averaged ensemble solution while those required by the data are reinforced. The ensemble of models can also be used to produce uncertainty estimates and experiments with synthetic data suggest that they represent actual uncertainty surprisingly well. We use the fast marching method in order to iteratively update the ray geometry and account for the non-linearity of the problem. The method is tested here with synthetic data in a 2-D application and compared with a subspace method that is a more standard matrix-based inversion scheme. Preliminary results illustrate the advantages of the reversible jump algorithm. A real data example is also shown where a tomographic image of Rayleigh wave group velocity for the Australian continent is constructed together with uncertainty estimates. [ABSTRACT FROM AUTHOR]

Focuses on the importance of inversion techniques in the Earth Sciences. Provision of indirect observation on the surface; Existence of the geometrical concept of algorithm annealing; Efficacy of inversion techniques.