Linear Inversion Scheme Research Articles

Imaging Upper-Mantle Anisotropy with Transdimensional Bayesian Monte Carlo Sampling

ABSTRACT Underdetermination is a condition affecting all problems in seismic imaging. It manifests mainly in the nonuniqueness of the models inferred from the data. This condition is exacerbated if simplifying hypotheses like isotropy are discarded in favor of more realistic anisotropic models that, although supported by seismological evidence, require more free parameters. Investigating the connections between underdetermination and anisotropy requires the implementation of solvers which explore the whole family of possibilities behind nonuniqueness and allow for more informed conclusions about the interpretation of the seismic models. Because these aspects cannot be investigated using traditional iterative linearized inversion schemes with regularization constraints that collapse the infinite possible models into a unique solution, we explore the application of transdimensional Bayesian Monte Carlo sampling to address the consequences of underdetermination in anisotropic seismic imaging. We show how teleseismic waves of P and S phases can constrain upper-mantle anisotropy and the amount of additional information these data provide in terms of uncertainty and trade-offs among multiple fields.

Sequential Seismic Anisotropic Inversion for VTI Media with Simulated Annealing Algorithm Aided by Adaptive Setting of Optimization Parameters

As an important geophysical data processing technique, seismic inversion estimates subsurface rock properties with seismic observations. However, anisotropic inversion, intended for a vertical transverse isotropy (VTI) media that primarily describes shale gas/oil resources, suffers from high nonlinearity. Simulated annealing is a widely used global optimization algorithm for solving nonlinear seismic inverse problems, but it involves multiple optimization parameters (e.g., initial temperature, search limit, and perturbation range). The importance of such parameters has been proven whilst the relevant analysis is limited in seismic inversion studies. This work hereby proposes a sequential anisotropic inversion method for VTI media, wherein we combine Bayesian linear and simulated annealing nonlinear inversion schemes. The simulated annealing is featured by adaptive optimization parameters aided by the linear result. Rather than the conventional method, the adaptive setting can be implemented trace by trace for complex reservoirs, which endows the method with enhanced stability and extended applicability. Synthetic tests and practical application demonstrate the validity of the method, wherein the obtained stiffness parameters facilitate the characterization of potential shale reservoirs with an improved accuracy.

The influence of physical and algorithmic factors on simulated far-field waveforms and source–time functions of underground explosions using unsupervised machine learning

SUMMARY Characterizing explosion sources and differentiating between earthquake and underground explosions using distributed seismic networks becomes non-trivial when explosions are detonated in cavities or heterogeneous ground material. Moreover, there is little understanding of how changes in subsurface physical properties affect the far-field waveforms we record and use to infer information about the source. Simulations of underground explosions and the resultant ground motions can be a powerful tool to systematically explore how different subsurface properties affect far-field waveform features, but there are added variables that arise from how we choose to model the explosions that can confound interpretation. To assess how both subsurface properties and algorithmic choices affect the seismic wavefield and the estimated source functions, we ran a series of 2-D axisymmetric non-linear numerical explosion experiments and wave propagation simulations that explore a wide array of parameters. We then inverted the synthetic far-field waveform data using a linear inversion scheme to estimate source–time functions (STFs) for each simulation case. We applied principal component analysis (PCA), an unsupervised machine learning method, to both the far-field waveforms and STFs to identify the most important factors that control variance in the waveform data and differences between cases. For the far-field waveforms, the largest variance occurs in the shallower radial receiver channels in the 0–50 Hz frequency band. For the STFs, both peak amplitude and rise times across different frequencies contribute to the variance. We find that the ground equation of state (i.e. lithology and rheology) and the explosion emplacement conditions (i.e. tamped versus cavity) have the greatest effect on the variance of the far-field waveforms and STFs, with the ground yield strength and fracture pressure being secondary factors. Differences in the PCA results between the far-field waveforms and STFs could possibly be due to near-field non-linearities of the source that are not accounted for in the estimation of STFs and could be associated with yield strength, fracture pressure, cavity radius and cavity shape parameters. Other algorithmic parameters are found to be less important and cause less variance in both the far-field waveforms and STFs, meaning algorithmic choices in how we model explosions are less important, which is encouraging for the further use of explosion simulations to study how physical Earth properties affect seismic waveform features and estimated STFs.

Rayleigh-wave ellipticity in weakly heterogeneous layered media

SUMMARY We derive approximate expressions for the ellipticity (i.e. horizontal-to-vertical or vertical-to-horizontal ratio) of Rayleigh waves propagating in a layered medium. The approximation is based on the generalized energy equation for Rayleigh waves, which has been used previously to obtain perturbational results for ellipticity. For a medium with weakly heterogeneous layers, we obtain an approximation from the perturbational result by taking the background medium to be homogeneous. The generalized energy equation also requires an auxiliary function and we discuss how the various possible functions are related to the homogeneous Rayleigh-wave eigenfunction. The analysis reveals that, within the weak approximation, the product of ellipticity and squared phase velocity is linearly related to squared shear wave velocity in the subsurface. We show the accuracy of the approximation with a simple layer-over-half-space model and then demonstrate its utility in a linear inversion scheme for shear wave velocity.

A nonlinear multiscale inversion approach for ambient noise tomography

SUMMARYAmbient noise tomography (ANT) has been considerably used in the last decade in both academic and industrial research. In this work, we propose an innovative technique for ANT based on nonlinear multiscale inversions. Our method relies on a progressive increase in the model parametrization to reduce the nonlinearity of the inverse problem. The developed method is compared with conventional inversion schemes (linear and nonlinear), using different regularization techniques and two different network configurations. The inversion is tested on 22 different synthetic models including classical checkerboard tests. Furthermore, we performed the inversion using real data from a campaign in 2018 at Cumbre Vieja volcano (Canary Islands). The results obtained on both network configurations show an improvement compared to conventional linear and nonlinear inversion schemes, especially when the ray path density is low. This technique does not require expensive computational resources, making it convenient for small-scale industrial applications, especially in the framework of geothermal exploration.

Joint PP and PS Pre-stack Seismic Inversion for Stratified Models Based on the Propagator Matrix Forward Engine

Pre-stack seismic inversion of the P- and S-wave velocities and bulk density is important in seismic exploration for evaluating lithological units and fluid properties. Generally, this inversion is based on ray-tracing modeling, which introduces errors and requires substantial pre-processing for stratified models due to its oversimplified single-interface assumption. To overcome those problems, we propose a pre-stack inversion method, using wave-equation-based modeling as a forward engine. Most wave-equation-based pre-stack inversions are based on the reflectivity method and adopt nonlinear optimization algorithms, although accurate, but computationally expensive. Hence, we use a fast propagator matrix (PM) method valid for layered media. To improve the stability and accuracy, the PP data inversion is extended to joint PP and PS PM-based inversion (JPMI). A linear inversion scheme is adopted to reduce the computational cost, and the Frechet derivatives are computed accordingly. Moreover, to obtain an optimal model solution, the L-BFGS (Limited-memory Broyden–Fletcher–Goldfarb–Shanno) optimization algorithm and L-curve criterion, an adaptive regularization parameter acquisition method, are implemented. A posterior probability analysis shows that the method has a higher parameter sensitivity than the joint exact Zoeppritz-based inversion and gives better estimations than the single-data inversion. We discuss the effects of dataset weight, internal multi-reflections, time window setting, noise level and initial model by using model tests. Synthetic and real-data examples demonstrate that the algorithm is better than the single PP inversion in terms of stability and accuracy, especially for S-wave velocity and density estimations.

Joint PP and PS pre-stack AVA inversion for VTI medium based on the exact Graebner equation

Unconventional hydrocarbon reservoirs, such as shale rocks, are generally transversely isotropic with a vertical axis of symmetry (VTI). The identification of sweet spots in these reservoirs, which must take into account this anisotropic character, can be performed with amplitude variations with offset/angle (AVO/AVA) inversion. However, most approaches use the Rüger reflection approximate equation as forward modeling, which has low accuracy at moderate to large angles, where the seismic gathers contain most of useful information on anisotropy parameters. In this paper, we propose an exact-equation-based seismic pre-stack AVA inversion for stable estimations of anisotropy parameters in VTI media. A linear inversion scheme is adopted to save computational time and the Fréchet derivatives, the key factors of the scheme, are derived here. To improve the stability and accuracy, the PP data inversion is extended to joint PP and PS inversion, termed joint exact-Graebner-based prestack inversion (JEGI). It is applied to both synthetic and field multi-component seismic data, showing its feasibility and effectiveness.

Experimental validation of D parameter model for droplet sizing using off-axis lidar measurements.

Information about the size distribution of liquid droplets in a fog can be retrieved by measuring the backscattering lidar depolarization parameter D in circular polarization. Using a polarimetric off-axis lidar, measurements at different backscattering angles are performed on fogs made of water droplets and of mineral oil. Estimation of the effective droplet size is obtained using constrained linear inversion. Mie theory is used to calculate the variation in depolarization parameters for different effective droplet sizes. The calculation is performed for various scattering angles. These calculations provide a kernel for the constrained linear inversion scheme. It is shown that the refractive index has an effect on the retrieved droplet sizes as well as the choice of scattering angles. These measurements confirm that the circular depolarization parameter measured near the backscattering angle can be modeled as a function of the forward-scattering diffraction peak. The results of the constrained linear inversion of measurements are consistent with in situ measurement of the droplet size distribution.

Prestack Waveform Inversion by Using an Optimized Linear Inversion Scheme

Seismic waveform inversion has been one of the most studied topics in recent years, which can be implemented both in common shot gather and prestack common midpoint gather. However, there are still several intractable issues to be addressed urgently, such as high-computational complexity, nonuniqueness, and robustness. In this paper, we have developed a model-based prestack waveform inversion (PWI) with a generalized propagation matrix scheme as forward operator, where a regularized function is minimized with the limited memory-Broyden–Fletcher–Goldfarb–Shanno technique to determine a model update corresponding to an adaptively determined regularization weight in each iterative step. To avoid falling into local extrema in the process of solving the objective function, we introduce an optimal transport method into the objective function to improve its convexity and use L-curve method to acquire the optimal regularization weight adaptively. The model tests show that the proposed scheme performs better than the conventional method significantly both on convergence and on accuracy. Furthermore, we apply the PWI with the proposed inversion scheme to the well-logging data and real seismic data. The results demonstrate that the proposed inversion scheme is not only capable of obtaining an accurate description of subsurface properties but also has a good convergence and robustness.

A reversible-jump Markov chain Monte Carlo algorithm for 1D inversion of magnetotelluric data

This paper presents a new computer code developed to solve the 1D magnetotelluric (MT) inverse problem using a Bayesian trans-dimensional Markov chain Monte Carlo algorithm. MT data are sensitive to the depth-distribution of rock electric conductivity (or its reciprocal, resistivity). The solution provided is a probability distribution - the so-called posterior probability distribution (PPD) for the conductivity at depth, together with the PPD of the interface depths. The PPD is sampled via a reversible-jump Markov Chain Monte Carlo (rjMcMC) algorithm, using a modified Metropolis-Hastings (MH) rule to accept or discard candidate models along the chains. As the optimal parameterization for the inversion process is generally unknown a trans-dimensional approach is used to allow the dataset itself to indicate the most probable number of parameters needed to sample the PPD. The algorithm is tested against two simulated datasets and a set of MT data acquired in the Clare Basin (County Clare, Ireland). For the simulated datasets the correct number of conductive layers at depth and the associated electrical conductivity values is retrieved, together with reasonable estimates of the uncertainties on the investigated parameters. Results from the inversion of field measurements are compared with results obtained using a deterministic method and with well-log data from a nearby borehole. The PPD is in good agreement with the well-log data, showing as a main structure a high conductive layer associated with the Clare Shale formation.In this study, we demonstrate that our new code go beyond algorithms developend using a linear inversion scheme, as it can be used: (1) to by-pass the subjective choices in the 1D parameterizations, i.e. the number of horizontal layers in the 1D parameterization, and (2) to estimate realistic uncertainties on the retrieved parameters. The algorithm is implemented using a simple MPI approach, where independent chains run on isolated CPU, to take full advantage of parallel computer architectures. In case of a large number of data, a master/slave appoach can be used, where the master CPU samples the parameter space and the slave CPUs compute forward solutions.

A theoretical model for fluvial channel response time during time‐dependent climatic and tectonic forcing and its inverse applications

AbstractThe fluvial response time dictates the duration of fluvial channel adjustment in response to changing climatic and tectonic conditions. However, when these conditions vary continuously, the channel cannot equilibrate and the response time is not well defined. Here I develop an analytical solution to the linear stream power model of fluvial incision that predicts the channel topography as a function of time‐dependent climatic and tectonic conditions. From this solution, a general definition of the fluvial response time emerges: the duration over which the tectonic history needs to be known to evaluate channel topography. This new definition is used in linear inversion schemes for inferring climatic or tectonic histories from river long profiles. The analytic solution further reveals that high‐frequency climatic oscillations, such as Milankovitch cycles, are not expected to leave significant fingerprints on the long profiles of fluvially incised detachment‐limited rivers.

On the use of sensitivity tests in seismic tomography

Sensitivity analysis with synthetic models is widely used in seismic tomography as a means for assessing the spatial resolution of solutions produced by, in most cases, linear or iterative nonlinear inversion schemes. The most common type of synthetic reconstruction test is the so-called checkerboard resolution test in which the synthetic model comprises an alternating pattern of higher and lower wave speed (or some other seismic property such as attenuation) in 2-D or 3-D. Although originally introduced for application to large inverse problems for which formal resolution and covariance could not be computed, these tests have achieved popularity, even when resolution and covariance can be computed, by virtue of being simple to implement and providing rapid and intuitive insight into the reliability of the recovered model. However, checkerboard tests have a number of potential drawbacks, including (1) only providing indirect evidence of quantitative measures of reliability such as resolution and uncertainty, (2) giving a potentially misleading impression of the range of scale-lengths that can be resolved, and (3) not giving a true picture of the structural distortion or smearing that can be caused by the data coverage. The widespread use of synthetic reconstruction tests in seismic tomography is likely to continue for some time yet, so it is important to implement best practice where possible. The goal of this paper is to develop the underlying theory and carry out a series of numerical experiments in order to establish best practice and identify some common pitfalls. Based on our findings, we recommend (1) the use of a discrete spike test involving a sparse distribution of spikes, rather than the use of the conventional tightly spaced checkerboard; (2) using data coverage (e.g. ray-path geometry) inherited from the model constrained by the observations (i.e. the same forward operator or matrix), rather than the data coverage obtained by solving the forward problem through the synthetic model; (3) carrying out multiple tests using structures of different scale length; (4) taking special care with regard to what can be inferred when using synthetic structures that closely mimic what has been recovered in the observation-based model; (5) investigating the range of structural wavelengths that can be recovered using realistic levels of imposed data noise; and (6) where feasible, assessing the influence of model parametrization error, which arises from making a choice as to how structure is to be represented. © The Authors 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society.

The role of fluid migration and static stress transfer in searching connections between the May 2012 Emilia earthquakes through a fully 3D finite element modeling

On May 20 and 29, 2012, two earthquakes struck the Emilia Romagna region (Northern Italy), with similar mechanisms. The proximity in space and time between the two events required to investigate possible triggering effects in terms of stress transfer or fluid migration. Moreover, the debate concerning fluid extraction–injection activities brought to the appointment of an International Commission (ICHESE), for evaluating relationships between the hydrocarbon exploitation and the seismic activity near the focal zone: the ICHESE conclusive report did not exclude the possibility of a triggering effect associated with explorations, although this appears not supported by robust evidences. On the other hand, there is no agreement between published papers addressing classical triggering mechanisms and in this respect an improved analysis based on 3D numerical modeling is a valuable tool to provide new insights and outcomes. To this aim, we built up a 3D model honoring the complexities of the real Earth, to retrieve the slip patterns on both the rupture planes by implementing finite element computed Green’s functions in a linear joint inversion scheme of geodetic data. Then, we inspected possible mechanisms which could have contributed to promote the second event, such as Coulomb stress transfer and natural fluid migration. Our findings suggest that both phenomena are not likely to be responsible for the May 29 earthquake, even if fluid migration was found to have a role in facilitating the rupture but without reaching hydrofracturing condition. Nevertheless a suprahydrostatic pore pressure is inferred which could have promoted the reactivation of the thrust fault.

The possibilities of least-squares migration of internally scattered seismic energy

Approximate images of the earth’s subsurface structures are usually obtained by migrating surface seismic data. Least-squares migration, under the single-scattering assumption, is used as an iterative linearized inversion scheme to suppress migration artifacts, deconvolve the source signature, mitigate the acquisition fingerprint, and enhance the spatial resolution of migrated images. The problem with least-squares migration of primaries, however, is that it may not be able to enhance events that are mainly illuminated by internal multiples, such as vertical and nearly vertical faults or salt flanks. To alleviate this problem, we adopted a linearized inversion framework to migrate internally scattered energy. We apply the least-squares migration of first-order internal multiples to image subsurface vertical fault planes. Tests on synthetic data demonstrated the ability of the proposed method to resolve vertical fault planes, which are poorly illuminated by the least-squares migration of primaries only. The proposed scheme is robust in the presence of white Gaussian observational noise and in the case of imaging the fault planes using inaccurate migration velocities. Our results suggested that the proposed least-squares imaging, under the double-scattering assumption, still retrieved the vertical fault planes when imaging the scattered data despite a slight defocusing of these events due to the presence of noise or velocity errors.

The 2009 L’Aquila earthquake coseismic rupture: open issues and new insights from 3D finite element inversion of GPS, InSAR and strong motion data

<p>We present a Finite Element inverse analysis of the static deformation field for the M<sub>w</sub>= 6.3, 2009 L’Aquila earthquake, in order to infer the rupture slip distribution on the fault plane. An univocal solution for the rupture slip distribution has not been reached yet with negative impact for reliable hazard scenarios in a densely populated area. In this study, Finite Element computed Green’s functions were implemented in a linear joint inversion scheme of geodetic (GPS and InSAR) and seismological (strong motion) coseismic deformation data. In order to fully exploit the informative power of our dense dataset and to honor the complexities of the real Earth, we implemented an optimized source model, represented by a fault plane subdivided in variable size patches, embedded in a high-resolution realistic three-dimensional model of the Apenninic seismo-tectonic setting, accounting for topographic reliefs and rheological heterogeneities deduced from local tomography. We infer that the investigated inversion domain contains two minima configurations in the solution space, i.e. a single- and a double-patch slip distribution, which are almost equivalent, so that the available datasets and numerical models are not able to univocally discriminate between them. Nevertheless our findings suggest that a two high-slip patch pattern is slightly favoured.</p>

A hybrid approach for estimating range-dependent properties of shallow water environments

H hybrid approach based on statistical signal characterization and a linear inversion scheme for the estimation of range dependent sound speed profiles of compact support in shallow water is presented. The approach is appropriate for ocean acoustic tomography when there is a single receiver available, as the first stage of the method is based on the statistical characterization of a single reception using wavelet transform to associate the signal with a set of parameters describing the statistical features of its wavelet sub-band coefficients. A non-linear optimization algorithm is then applied to associate these features with range-dependent sound speed profile in the water column. This inversion method is restricted to cases where the range dependency is of compact support. At the second stage a linear inversion scheme based on modal arrivals identification and a first order perturbation formula to associate sound speed differences with modal travel time perturbations is applied to fine tune the results obtained by the optimization scheme. A second restriction of this stage is that mode identification is necessary. If this assumption is fulfilled the whole scheme may be applied in ocean acoustic tomography for the retrieval of three-dimensional features, combining inversion results at various slices.

Tectonics from fluvial topography using formal linear inversion: Theory and applications to the Inyo Mountains, California

AbstractTectonic activity generates topography, and the variability of tectonic forcing is responsible for topographic patterns and variability of relief in fluvial landscapes. Despite this basic relation, the inverse problem, by which features of the topography are used for inferring tectonic uplift rates, has proven challenging. Here we develop formal linear inversion schemes to infer a record of the rate of relative uplift as a function of space and time from the long profiles of rivers. The relative uplift rate is the difference between the rates of rock uplift and of the base level change. The inversion schemes are based on a closed‐form analytic solution to the transient linear stream power model, and to increase model resolution they make use of the multiplicity of information made available by multiple rivers and their tributaries. The distribution of the fluvial response time to tectonic perturbations is a key component of the inversion scheme, as this determines which tectonic events are preserved in the topography. We develop two inversion parameterizations that differ in their assumptions about the tectonic forcing: space‐invariant and time‐space variability with an assumed spatial distribution. The inversion schemes are applied to the Inyo Mountains, an uplifted block along the western boundary of the Basin and Range Province in California. Inversion results indicate that the range has been experiencing an acceleration of the relative uplift in the past ∼2‐3 Ma. We use the inversion results to constrain the paleotopography and paleo‐erosion rate along the range and to recover the throw rate history along the fault that bounds the Inyo range.

Crustal and uppermost mantle structure of SE Tibetan plateau from Rayleigh-wave group-velocity measurements

A shear-wave velocity model of the crust and uppermost mantle beneath the SE Tibetan plateau was derived by inverting Rayleigh-wave group-velocity measurements of periods between 10 and 70 s. Rayleigh-wave group-velocity dispersions along more than 3,000 interstation paths were measured based on analysis of teleseismic waveform data recorded by temporary seismic stations. These observations were then utilized to construct 2D group-velocity maps in the period range of 10–70 s. The new group-velocity maps have an enhanced resolution compared with previous global and regional group-velocity models in this region because of the denser and more uniform data coverage. The lateral resolution across the region is about 0.5° for the periods used in this study. Local dispersion curves were then inverted for a 3D shear-wave velocity model of the region by applying a linear inversion scheme. Our 3D shear-wave model confirms the presence of low-velocity zones (LVZs) in the crust beneath the northern part of this region. Our imaging shows that the upper-middle crustal LVZ beneath the Tengchong region is isolated from these LVZs beneath the eastern and northern part of this region. The upper–middle crustal LVZ may be regarded as evidence of a magma chamber in the crust beneath the Tengchong Volcanoes. Our model also reveals a slow lithospheric structure beneath Tengchong and a fast shield-like mantle beneath the stable Yangtze block.

MSATSI: A MATLAB Package for Stress Inversion Combining Solid Classic Methodology, a New Simplified User-Handling, and a Visualization Tool

Online Material: Figures of complete stress inversion results. The estimation of the stress‐field orientation from focal mechanisms of earthquakes is a relevant tool to understand crustal mechanics and the physics of earthquakes. In global seismology, Formal Stress Inversion (FSI) is a well‐established technique to study tectonic processes associated with the occurrence of large earthquakes (e.g., Hardebeck and Michael, 2004; Yoshida et al. , 2012). Information on the stress‐field orientation is also relevant for the exploitation of hydrocarbon and geothermal reservoirs. Knowledge of the orientation of the maximum horizontal stress ( σ Hmax) is crucial for reservoir development, such as drilling and leakage risk assessment (Terakawa et al. , 2012; Martinez‐Garzon et al. , 2013). Additionally, the FSI technique may be useful for understanding the physics of rupture processes at a microscale in the frame of rock deformation experiments in the laboratory. In seismological practice, the stress tensor is obtained either from inverting earthquake focal mechanisms or directly from first‐motion polarities. Most of the developed FSI methods share two main assumptions: Estimation of stress‐field orientation is a nonlinear inverse problem that can be solved either directly or linearized. When solving the nonlinear inverse problem, the best‐fitting stress tensor to the group of focal mechanisms is obtained using grid‐search algorithms (Gephart and Forsyth, 1984; Gephart, 1990; Arnold and Townend, 2007) or Monte Carlo sampling‐based optimization methods (Angelier, 1984; Xu, 2004). The linearized inversion scheme is solved by a generally less computationally extensive least‐squares technique (Michael, 1984; Hardebeck and Michael, 2006). Because of …

Improvement of local anisotropy estimation from VSP data through experimental design

We investigate a linear inversion scheme to determine the local anisotropy around a well using multi-azimuthal walkaway VSP data. The input data consist of the vertical components of the slowness vector and the polarization vector of direct and reflected P-wave measurements. The inversion assumes that weak anisotropy can be modelled by first-order perturbations around an isotropic reference medium. This inversion scheme is not restricted to symmetrical classes of anisotropy and the well local orientation. Based on numerical simulations, we show that the inversion is sensitive to the number and orientation of the acquired azimuths. Using conservative assumptions regarding the noise level, phase velocity can be well estimated for a limited angular opening around the well local orientation.