How do we understand and visualize uncertainty?

Geophysicists are often concerned with reconstructing subsurface properties using observations collected at or near the surface. For example, in seismic migration, we attempt to reconstruct subsurface geometry from surface seismic recordings, and in potential field inversion, observations are used to map electrical conductivity or density variations in geologic layers. The procedure of inferring information from indirect observations is called an inverse problem by mathematicians, and such problems are common in many areas of the physical sciences. The inverse problem of inferring the subsurface using surface observations has a corresponding forward problem, which consists of determining the data that would be recorded for a given subsurface configuration. In the seismic case, forward modeling involves a method for calculating a synthetic seismogram, for gravity data it consists of a computer code to compute gravity fields from an assumed subsurface density model. Note that forward modeling often involves assumptions about the appropriate physical relationship between unknowns (at depth) and observations on the surface, and all attempts to solve the problem at hand are limited by the accuracy of those assumptions. In the broadest sense then, exploration geophysicists have been engaged in inversion since the dawn of the profession and indeed algorithms often applied in processing centers can all be viewed as procedures to invert geophysical data

Geophysical imaging using trans-dimensional trees.

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 trans-dimensional 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

pyprop8: A lightweight code to simulate seismic observables in a layered half-space

The package pyprop8 enables calculation of the response of a 1-D layered halfspace to a seismic source, and also derivatives (‘sensitivity kernels’) of the wavefield with respect to source parameters. Seismograms, seismic spectra, and measures of static displacement (e.g. GPS, InSAR and field observations) may all be simulated. The method is based on a ThompsonHaskell propagator matrix algorithm, described in O’Toole & Woodhouse (2011) and O’Toole et al. (2012). The package is entirely written in Python, dependent only on the mainstream libraries numpy (Harris et al., 2020) and scipy (Virtanen et al., 2020). As such, it is lightweight and easy to deploy across a variety of platforms, making it particularly suited to use for teaching and outreach purposes

Overcomplete tomography: a novel approach to imaging

Regularized least-squares tomography offers a straightforward and efficient imaging method and has seen extensive application across various fields. However, it has a few drawbacks, such as (i) the regularization imposed during the inversion tends to give a smooth solution, which will fail to reconstruct a multi-scale model well or detect sharp discontinuities, (ii) it requires finding optimum control parameters, and (iii) it does not produce a sparse solution. This paper introduces ‘overcomplete tomography’, a novel imaging framework that allows high-resolution recovery with relatively few data points. We express our image in terms of an overcomplete basis, allowing the representation of a wide range of features and characteristics. Following the insight of ‘compressive sensing’, we regularize our inversion by imposing a penalty on the L1 norm of the recovered model, obtaining an image that is sparse relative to the overcomplete basis. We demonstrate our method with a synthetic and a real X-ray tomography example. Our experiments indicate that we can reconstruct a multi-scale model from only a few observations. The approach may also assist interpretation, allowing images to be decomposed into (for example) ‘global’ and ‘local’ structures. The framework presented here can find application across a wide range of fields, including engineering, medical and geophysical tomography

Bayesian noise-reduction in Arabia/Somalia and Nubia/Arabia finite rotations since ~20 Ma: Implications for Nubia/Somalia relative motion

Knowledge of Nubia/Somalia relative motion since the Early Neogene is of particular importance in the Earth Sciences, because it (i) impacts on inferences on African dynamic topography; and (ii) allows us to link plate kinematics within the Indian realm with those within the Atlantic basin. The contemporary Nubia/Somalia motion is well known from geodetic observations. Precise estimates of the past-3.2-Myr average motion are also available from paleo-magnetic observations. However, little is known of the Nubia/Somalia motion prior to ∼3.2 Ma, chiefly because the Southwest Indian Ridge spread slowly, posing a challenge to precisely identify magnetic lineations. This also makes the few observations available particularly prone to noise. Here we reconstruct Nubia/Somalia relative motions since ∼20 Ma from the alternative plate-circuit Nubia-Arabia-Somalia. We resort to trans-dimensional hierarchical Bayesian Inference, which has proved effective in reducing finite-rotation noise, to unravel the Arabia/Somalia and Arabia/Nubia motions. We combine the resulting kinematics to reconstruct the Nubia/Somalia relative motion since ∼20 Ma. We verify the validity of the approach by comparing our reconstruction with the available record for the past ∼3.2 Myr, obtained through Antarctica. Results indicate that prior to ∼11 Ma the total motion between Nubia and Somalia was faster than today. Furthermore, it featured a significant strike-slip component along the Nubia/Somalia boundary. It is only since ∼11 Ma that Nubia diverges away from Somalia at slower rates, comparable to the present-day one. Kinematic changes of some 20% might have occurred in the period leading to the present-day, but plate-motion steadiness is also warranted within the uncertainties. Key Points Bayesian inference to reduce finite-rotation noise Reconstruction of Nubia/Somalia motion since ∼20 Ma Nubia/Somalia relative motion changed at ∼11 M

A self-parametrizing partition model approach to tomographic inverse problems

Partition modelling is a statistical method for nonlinear regression and classification, and is particularly suited to dealing with spatially variable parameters. Previous applications include disease mapping in medical statistics. Here we extend this method to the seismic tomography problem. The procedure involves a dynamic parametrization for the model which is able to adapt to an uneven spatial distribution of the information on the model parameters contained in the observed data. The approach provides a stable solution with no need for explicit regularization, i.e. there is neither user supplied damping term nor tuning of trade-off parameters. The method is an ensemble inference approach within a Bayesian framework. Many potential solutions are generated, and information is extracted from the ensemble as a whole. In terms of choosing a single model, it is straightforward to perform Monte Carlo integration to produce the expected Earth model. The inherent model averaging process naturally smooths out unwarranted structure in the Earth model, but maintains local discontinuities if well constrained by the data. Calculation of uncertainty estimates is also possible using the ensemble of models, and experiments with synthetic data suggest that they are good representations of the true uncertainty

Trans-dimensional inverse problems, model comparison and the evidence

In most geophysical inverse problems the properties of interest are parametrized using a fixed number of unknowns. In some cases arguments can be used to bound the maximum number of parameters that need to be considered. In others the number of unknowns is set at some arbitrary value and regularization is used to encourage simple, non-extravagant models. In recent times variable or self-adaptive parametrizations have gained in popularity. Rarely, however, is the number of unknowns itself directly treated as an unknown. This situation leads to a transdimensional inverse problem, that is, one where the dimension of the parameter space is a variable to be solved for. This paper discusses trans-dimensional inverse problems from the Bayesian viewpoint. A particular type of Markov chain Monte Carlo (MCMC) sampling algorithm is highlighted which allows probabilistic sampling in variable dimension spaces. A quantity termed the evidence or marginal likelihood plays a key role in this type of problem. It is shown that once evidence calculations are performed, the results of complex variable dimension sampling algorithms can be replicated with simple and more familiar fixed dimensional MCMC sampling techniques. Numerical examples are used to illustrate the main points. The evidence can be difficult to calculate, especially in high-dimensional non-linear inverse problems. Nevertheless some general strategies are discussed and analytical expressions given for certain linear problem