A practical regularization for seismic tomography

Regularization is usually necessary in solving seismic tomographic inversion problems. In general the equation system of seismic tomography is very large, often making a suitable choice of the regularization parameter difficult. In this paper, we propose an algorithm for the practical choice of the regularization parameter in linear tomographic inversion. The algorithm is based on the types of statistical assumptions most commonly used in seismic tomography. We first transfer the system of equations into a Krylov subspace by using Lanczos bidiagonalization. In the transformed subspace, the system of equations is then changed into the form of a standard damped least squares normal equation. The solution to this normal equation can be written as an explicit function of the regularization parameter, which makes the choice of the regularization computationally convenient. Two criteria for the choice of the regularization parameter are investigated with the numerical simulations. If the dimensions of the transformed space are much less than that of the original model space, the algorithm can be very computationally efficient, which is practically useful in large seismic tomography problems.     

Sequential integrated inversion of refraction and wide-angle reflection traveltimes and gravity data for two-dimensional velocity structures

A new algorithm is presented for the integrated 2-D inversion of seismic traveltime and gravity data. The algorithm adopts the 'maximum likelihood' regularization scheme. We construct a 'probability density function' which includes three kinds of information: information derived from gravity measurements; information derived from the seismic traveltime inversion procedure applied to the model; and information on the physical correlation among the density and the velocity parameters. We assume a linear relation between density and velocity, which can be node-dependent; that is, we can choose different relationships for different parts of the velocity–density grid. In addition, our procedure allows us to consider a covariance matrix related to the error propagation in linking density to velocity. We use seismic data to estimate starting velocity values and the position of boundary nodes. Subsequently, the sequential integrated inversion (SII) optimizes the layer velocities and densities for our models. The procedure is applicable, as an additional step, to any type of seismic tomographic inversion.
We illustrate the method by comparing the velocity models recovered from a standard seismic traveltime inversion with those retrieved using our algorithm. The inversion of synthetic data calculated for a 2-D isotropic, laterally inhomogeneous model shows the stability and accuracy of this procedure, demonstrates the improvements to the recovery of true velocity anomalies, and proves that this technique can efficiently overcome some of the limitations of both gravity and seismic traveltime inversions, when they are used independently.
An interpretation of field data from the 1994 Vesuvius test experiment is also presented. At depths down to 4.5 km, the model retrieved after a SII shows a more detailed structure than the model obtained from an interpretation of seismic traveltime only, and yields additional information for a further study of the area.     

Inversion of seismic attributes for velocity and attenuation structure

We have developed an inversion formuialion for velocity and attenuation structure using seismic attributes, including envelope amplitude, instantaneous frequency and arrival times of selected seismic phases. We refer to this approach as AFT inversion for amplitude, (instantaneous) frequency and time. Complex trace analysis is used to extract the different seismic attributes. The instantaneous frequency data are converted to t * using a matching procedure that approximately removes the effects of the source spectra. To invert for structure, ray-perturbation methods are used to compute the sensitivity of the seismic attributes to variations in the model. An iterative inversion procedure is then performed from smooth to less smooth models that progressively incorporates the shorter-wavelength components of the model. To illustrate the method, seismic attributes are extracted from seismic-refraction data of the Ouachita PASSCAL experiment and used to invert for shallow crustal velocity and attenuation structure. Although amplitude data are sensitive to model roughness, the inverted velocity and attenuation models were required by the data to maintain a relatively smooth character. The amplitude and t * data were needed, along with the traveltimes, at each step of the inversion in order to fit all the seismic attributes at the final iteration.     

Wavepath traveltime tomography

The elastic-wave equation is used to construct sensitivity kernels relating perturbations in elastic parameters to traveltime deviations. Computation of the functions requires a correlation of the forward-propagating seismic wavefield with a backward propagation of the residual wavefield. The computation of the wavefields is accomplished using a finite difference algorithm and is efficiently executed on a CM-2 parallel processor. The source and receiver locations have maximum sensitivity to velocity structure. The sensitivity kernels or wavepaths are well suited for transmission traveltime inversion such as cross-borehole tomography and vertical seismic profiling. Conventional ray tomography and wavepath tomography are applied to a set of P -wave arrival times, from a cross-borehole experiment at Kesterson, California. Because the wavepaths have increased sensitivity near the source and receiver there are differences in resolution of the velocity structure. Both techniques recover the same relative variations in velocity where the coverage is adequate. The wavepath solution is more laterally continuous and the dominant variation is vertical, as is expected for the layered sediments in this region.     

P- and S-velocity images of the lithosphere–asthenosphere system in the Central Andes from local-source tomographic inversion

About 50 000 P and S arrival times and 25 000 values of t * recorded at seismic arrays operated in the Central Andes between 20°S and 25°S in the time period from 1994 to 1997 have been used for locating more than 1500 deep and crustal earthquakes and creating 3-D P , S velocity and Qp models. The study volume in the reference model is subdivided into three domains: slab, continental crust and mantle wedge. A starting velocity distribution in each domain is set from a priori information: in the crust it is based on the controlled sources seismic studies; in slab and mantle wedge it is defined using relations between P and S velocities, temperature and composition given by mineral physics. Each iteration of tomographic inversion consists of the following steps: (1) absolute location of sources in 3-D velocity model using P and S arrival times; (2) double-difference relocation of the sources and (3) simultaneous determination of P and S velocity anomalies, P and S station corrections and source parameters by inverting one matrix. Velocity parameters are computed in a mesh with the density of nodes proportional to the ray density with double-sided nodes at the domain boundaries. The next iteration is repeated with the updated velocity model and source parameters obtained at the previous step. Different tests aimed at checking the reliability of the obtained velocity models are presented. In addition, we present the results of inversion for Vp and Vp/Vs parameters, which appear to be practically equivalent to Vp and Vs inversion. A separate inversion for Qp has been performed using the ray paths and source locations in the final velocity model. The resulting Vp , Vs and Qp distributions show complicated, essentially 3-D structure in the lithosphere and asthenosphere. P and S velocities appear to be well correlated, suggesting the important role of variations of composition, temperature, water content and degree of partial melting.     

Crosshole seismic waveform tomography – I. Strategy for real data application

The frequency-domain version of waveform tomography enables the use of distinct frequency components to adequately reconstruct the subsurface velocity field, and thereby dramatically reduces the input data quantity required for the inversion process. It makes waveform tomography a computationally tractable problem for production uses, but its applicability to real seismic data particularly in the petroleum exploration and development scale needs to be examined. As real data are often band limited with missing low frequencies, a good starting model is necessary for waveform tomography, to fill in the gap of low frequencies before the inversion of available frequencies. In the inversion stage, a group of frequencies should be used simultaneously at each iteration, to suppress the effect of data noise in the frequency domain. Meanwhile, a smoothness constraint on the model must be used in the inversion, to cope the effect of data noise, the effect of non-linearity of the problem, and the effect of strong sensitivities of short wavelength model variations. In this paper we use frequency-domain waveform tomography to provide quantitative velocity images of a crosshole target between boreholes 300 m apart. Due to the complexity of the local geology the velocity variations were extreme (between 3000 and 5500 m s⁻¹), making the inversion problem highly non-linear. Nevertheless, the waveform tomography results correlate well with borehole logs, and provide realistic geological information that can be tracked between the boreholes with confidence.     

Crosshole seismic waveform tomography – II. Resolution analysis

In an accompanying paper, we used waveform tomography to obtain a velocity model between two boreholes from a real crosshole seismic experiment. As for all inversions of geophysical data, it is important to make an assessment of the final model, to determine which parts of the model are well-resolved and can confidently be used for geological interpretation. In this paper we use checkerboard tests to provide a quantitative estimate of the performance of the inversion and the reliability of the final velocity model. We use the output from the checkerboard tests to determine resolvability across the velocity model. Such tests can act as good guides for designing appropriate inversion strategies. Here we discovered that, by including both reference-model and smoothing constraints in initial inversions, and then relaxing the smoothing constraint for later inversions, an optimum velocity image was obtained. Additionally, we noticed that the performance of the inversion was dependent on a relationship between velocity perturbation and checkerboard grid-size: larger velocity perturbations were better-resolved when the grid-size was also increased. Our results suggest that model assessment is an essential step prior to interpreting features in waveform tomographic images.     

On the resolving power of tomographic images in the Aegean area

b
The imaging of upper mantle heterogeneity by seismic tomography is strongly limited by the uneven global distribution both of seismic recording stations and earthquake sources. This can result in a loss of resolution and significance in the final image, particularly when a sparse data set contains few ray paths which intersect at sufficiently high angles in the volume of interest. In order to investigate the theoretical resolving power of a previously published tomographic image of the Aegean area, synthetic tests of the inversion procedure using a ray-path matrix obtained in this previous study for local and teleseismic P -waves were carried out. The aim was to examine the extent to which the shape of a synthetic lithospheric slab penetrating to different depths is inherently distorted by the tomographic imaging procedure, and to compare the synthetic tomographic images with the results from the actual inversion. The distortion is found to take the form of an artificial stretching of the lithospheric slab. The maximum 'stretching factor', as indicated by the downdip displacement of the peak amplitude of the synthetic high-velocity anomaly, is found to be a factor of 2 or so, though the distortion is usually less than this. The peak amplitude of the tomographic image of a lithospheric slab is found from the inversion of traveltime data to be at depths at or below 400 km. This indicates that the high-velocity lithospheric slab in the Aegean penetrates deeper than the Benioff zone seismicity of about 200 km. However, no constraints of the maximum depth of penetration could be established with the data set used in the present work.     

Seismic tomographic imaging of P- and S-waves velocity perturbations in the upper mantle beneath Iran

The inverse tomography method has been used to study the P - and S -waves velocity structure of the crust and upper mantle underneath Iran. The method, based on the principle of source–receiver reciprocity, allows for tomographic studies of regions with sparse distribution of seismic stations if the region has sufficient seismicity. The arrival times of body waves from earthquakes in the study area as reported in the ISC catalogue (1964–1996) at all available epicentral distances are used for calculation of residual arrival times. Prior to inversion we have relocated hypocentres based on a 1-D spherical earth's model taking into account variable crustal thickness and surface topography. During the inversion seismic sources are further relocated simultaneously with the calculation of velocity perturbations. With a series of synthetic tests we demonstrate the power of the algorithm and the data to reconstruct introduced anomalies using the ray paths of the real data set and taking into account the measurement errors and outliers. The velocity anomalies show that the crust and upper mantle beneath the Iranian Plateau comprises a low velocity domain between the Arabian Plate and the Caspian Block. This is in agreement with global tomographic models, and also tectonic models, in which active Iranian plateau is trapped between the stable Turan plate in the north and the Arabian shield in the south. Our results show clear evidence of the mainly aseismic subduction of the oceanic crust of the Oman Sea underneath the Iranian Plateau. However, along the Zagros suture zone, the subduction pattern is more complex than at Makran where the collision of the two plates is highly seismic.     

Sensitivities of seismic traveltimes and amplitudes in reflection tomography

Seismic traveltimes and amplitudes in reflection-seismic data show different dependences on the geometry of reflection interfaces, and on the variation of interval velocities. These dependences are revealed by eigenanalysis of the Hessian matrix, defined in terms of the Fréchet matrix and its adjoint associated with different norms chosen in the model space. The eigenvectors and eigenvalues of the Hessian clearly show that for reflection tomographic inversion, traveltime and amplitude data contain complementary information. Both for reflector-geometry and for interval-velocity variations, the traveltimes are sensitive to the model components with small wavenumbers, whereas the amplitudes are more sensitive to the components with high wavenumbers. The model resolution matrices, after the rejection of eigenvectors corresponding to small eigenvalues, give us some insight into how the addition of amplitude information could potentially contribute to the recovery of physical parameters.
In order to cooperatively invert seismic traveltimes and amplitudes simultaneously, we propose an empirical definition of the data covariance matrix which balances the relative sensitivities of different types of data. We investigate the cooperative use of both data types for, separately, interface-geometry and 2-D interval-velocity variations. In both cases we find that cooperative inversions can provide better solutions than those using traveltimes alone. The potential benefit of including amplitude-data constraints in seismic-reflection traveltime tomography is therefore that it may be possible to resolve the known ambiguity between the reflector-depth uncertainty and the interval-velocity uncertainty better.     

Application of stacking and inversion techniques to three-dimensional wide-angle reflection and refraction seismic data of the Eastern Alps

We present new methods for the interpretation of 3-D seismic wide-angle reflection and refraction data with application to data acquired during the experiments CELEBRATION, 2000 and ALP 2002 in the area of the Eastern Alps and their transition to the surrounding tectonic provinces (Bohemian Massif, Carpathians, Pannonian domain, Dinarides). Data was acquired on a net of arbitrarily oriented seismic lines by simultaneous recording on all lines of seismic waves from the shots, which allows 2-D and 3-D interpretations. Much (80%) of the data set consists of crossline traces. Low signal to noise (S/N) ratio in the area of the young orogens decreases the quality of travel time picks. In these seismically heterogeneous areas it is difficult to assign clearly defined arrivals to the seismic phases, in particular on crossline record sections.
In order to enhance the S/N ratio, signal detection and stacking techniques have been applied to enhance the Pg -, Pn - and PmP phases. Further, inversion methods have been developed for the interpretation of WAR/R-data, based on automated 1-D inversion ( Pg ) and the application of the delay time concept ( Pn ). The results include a 3-D velocity model of the crust based on Pg waves, time and depth maps of the Moho and a Pn -velocity map. The models based on stacked data are robust and provide a larger coverage, than models based on travel time picks from single-fold (unstacked) traces, but have relatively low resolution, especially near the surface. They were used as the basis for constructing models with improved resolution by the inversion of picks from single-fold data. The results correlate well with geological structures and show new prominent features in the Eastern Alps area and their surrounds. The velocity distribution in the crust has strong lateral variations and the Moho in the investigation area appears to be fragmented into three parts.     

Upper mantle beneath the Eger Rift (Central Europe): plume or asthenosphere upwelling?

We present the first results of a high-resolution teleseismic traveltime tomography and seismic anisotropy study of the lithosphere–asthenosphere system beneath the western Bohemian Massif. The initial high-resolution tomography down to a depth of 250 km did not image any columnar low-velocity anomaly which could be interpreted as a mantle plume anticipated beneath the Eger Rift, similar to recent findings of small plumes beneath the French Massif Central and the Eifel in Germany. Alternatively, we interpret the broad low-velocity anomaly beneath the Eger Rift by an upwelling of the lithosphere–asthenosphere transition. We also map lateral variations of seismic anisotropy of the mantle lithosphere from spatial variations of P -wave delay times and the shear wave splitting. Three major domains characterised by different orientations of seismic anisotropy correspond to the major tectonic units—Saxothuringian, Moldanubian and the Teplá-Barrandian—and their fabrics fit to those found in our previous studies of mantle anisotropy on large European scales.     

Inversion of controlled-source seismic tomography and gravity data with the self-adaptive wavelet parametrization of velocities and interfaces

A self-adaptive automated parametrization approach is suggested for the sequential inversion of controlled-source seismic tomography and gravity data. The velocities and interfaces are parametrized by their Haar wavelet expansion coefficients. Only those coefficients that are well constrained by the data, as measured by the number of rays that cross the corresponding wavelet function support area and their angular coverage, are inverted for, others are set to zero. This approach results in a reasonable distribution of resolution throughout the model even in cases of irregular ray coverage and does overcome the trade-off between different types of model parameters. A modified sequential inversion approach is suggested to join the traveltimes and gravity anomalies inversion. An algorithm is developed that inverts for smooth velocity and density variations inside the seismic layer, the position of its bottom interface as well as for optimal values of the velocity-to-density regression coefficients. The algorithm makes use of direct (diving), reflected and head (critically refracted) wave traveltimes. The algorithm workflow is demonstrated on a synthetic data example.     

A transformational approach to geophysical inverse problems

A set of coordinate transformations is used to linearize a general geophysical inverse problem. Statistical and analytic techniques are employed to estimate the parameters of such linearization transformations. In the transformed space, techniques from linear inverse theory may be utilized. Consequently, important concepts, such as model parameter covariance, model parameter resolution and averaging kernels, may be carried over to non-linear inverse problems. I apply the approach to a set of seismic cross-borehole traveltimes gathered at the Conoco Borehole Test Facility. the seismic survey was conducted within the Fort Riley formation, a limestone with thin interbedded shales. Between the boreholes, the velocity structure of the Fort Riley formation consists of a high-velocity region overlying a section of lower velocity. It is found that model parameter resolution is poorest and spatial averaging lengths are greatest in the underlying low-velocity region.     

A crustal seismic velocity model for the UK, Ireland and surrounding seas

A regional model of the 3-D variation in seismic P -wave velocity structure in the crust of NW Europe has been compiled from wide-angle reflection/refraction profiles. Along each 2-D profile a velocity–depth function has been digitised at 5 km intervals. These 1-D velocity functions were mapped into three dimensions using ordinary kriging with weights determined to minimise the difference between digitised and interpolated values. An analysis of variograms of the digitised data suggested a radial isotropic weighting scheme was most appropriate. Horizontal dimensions of the model cells are optimised at 40 × 40 km and the vertical dimension at 1 km. The resulting model provides a higher resolution image of the 3-D variation in seismic velocity structure of the UK, Ireland and surrounding areas than existing models. The construction of the model through kriging allows the uncertainty in the velocity structure to be assessed. This uncertainty indicates the high density of data required to confidently interpolate the crustal velocity structure, and shows that for this region the velocity is poorly constrained for large areas away from the input data.     

Length scales, patterns and origin of azimuthal seismic anisotropy in the upper mantle as mapped by Rayleigh waves

We measure the degree of consistency between published models of azimuthal seismic anisotropy from surface waves, focusing on Rayleigh wave phase-velocity models. Some models agree up to wavelengths of ∼2000 km, albeit at small values of linear correlation coefficients. Others are, however, not well correlated at all, also with regard to isotropic structure. This points to differences in the underlying data sets and inversion strategies, particularly the relative 'damping' of mapped isotropic versus anisotropic anomalies. Yet, there is more agreement between published models than commonly held, encouraging further analysis. Employing a generalized spherical harmonic representation, we analyse power spectra of orientational (2Ψ) anisotropic heterogeneity from seismology. We find that the anisotropic component of some models is characterized by stronger short-wavelength power than the associated isotropic structure. This spectral signal is consistent with predictions from new geodynamic models, based on olivine texturing in mantle flow. The flow models are also successful in predicting some of the seismologically mapped patterns. We substantiate earlier findings that flow computations significantly outperform models of fast azimuths based on absolute plate velocities. Moreover, further evidence for the importance of active upwellings and downwellings as inferred from seismic tomography is presented. Deterministic estimates of expected anisotropic structure based on mantle flow computations such as ours can help guide future seismologic inversions, particularly in oceanic plate regions. We propose to consider such a priori information when addressing open questions about the averaging properties and resolution of surface and body wave based estimates of anisotropy.     

Velocity structure of upper-mantle transition zones beneath central eurasia from seismic inversion using genetic algorithms

We present velocity constraints for the upper-mantle transition zones beneath Central Siberia based on observations of the 1982 RIFT Deep Seismic Sounding (DSS) profile. The data consist of seismic recordings of a nuclear explosion in north-western Siberia along a 2600 km long seismic profile extending from the Yamal Peninsula to Lake Baikal. We invert seismic data from the mantle transition zones using a non-linear inversion scheme using a genetic algorithm for optimization and the WKBJ method to compute the synthetic seismograms. A statistical error analysis using a graph-binning technique was performed to provide uncertainty values in the velocity models.
Our best model for the upper-mantle velocity discontinuity near 410 km depth has a two-stage velocity-gradient structure, with velocities increasing from 8.70–9.25 km s⁻¹ over a depth range of 400–415 km, a gradient of 0.0433 s⁻¹, and from 9.25–9.60 km s⁻¹ over a depth range of 415–435 km, a gradient of 0.0175 s⁻¹. This derived model is consistent with other seismological observations and mineral-physics models. The model for the velocity discontinuity near 660 km depth is simple, sharp and includes velocities increasing from 10.15 km s⁻¹ at 655 km depth to 10.70 km s⁻¹ at 660 km depth, a gradient of 0.055 s⁻¹.     

Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature

We present a regional surface waveform tomography of the Pacific upper mantle, obtained using an automated multimode surface waveform inversion technique on fundamental and higher mode Rayleigh waves, to constrain the   V_SV   structure down to ∼400 km depth. We have improved on previous implementations of this technique by robustly accounting for the effects of uncertainties in earthquake source parameters in the tomographic inversion. We have furthermore improved path coverage in the South Pacific region by including Rayleigh wave observations from the French Polynesian Pacific Lithosphere and Upper Mantle Experiment deployment. This improvement has led to imaging of vertical low-velocity structures associated with hotspots within the South Pacific Super-Swell region. We have produced an age-dependent average cross-section for the Pacific Ocean lithosphere and found that the increase in   V_SV   with age is broadly compatible with a half-space cooling model of oceanic lithosphere formation. We cannot confirm evidence for a Pacific-wide reheating event. Our synthetic tests show that detailed interpretation of average   V_SV   trends across the Pacific Ocean may be misleading unless lateral resolution and amplitude recovery are uniform across the region, a condition that is difficult to achieve in such a large oceanic basin with current seismic stations.     

Large-magnitude Central American earthquakes, 1898–1994

We have collected and re-examined macroseismic information for large Central American earthquakes since the beginning of the period of instrumental recording about one hundred years ago, and combined this with a reassessment of early instrumental information to produce a catalogue of 51 events that, we believe includes ail those with magnitudes ( Ms ) greater than 7.0. We have reassessed surface-wave magnitudes by consulting station bulletins and we have derived a correction that gives an equivalent Ms for events of intermediate depth. We have also developed a regional relationship between Ms and seismic moment, which enables us to estimate the seismic slip rate across the Middle American Trench. Our best estimates give an average slip rate several times smaller than suggested convergence rates, but with the seismic slip in the central segment of the trench almost an order of magnitude smaller than that in the segments on either side. The low seismic slip rate may indicate aseismic crustal deformation     

Modelling and inversion of single-ended refraction data from the shot gathers of multifold deep seismic reflection profiling—an approach for deriving the shallow velocity structure

A multifold crustal-scale deep seismic near-vertical reflection profile generates a large number of single-ended shot gathers, which provide redundant data sets because of overlapping coverage of the shallow refractors. We present an approach for deriving the shallow velocity structure by modelling and inversion of single-ended seismic refraction first arrival traveltime data. We apply this method to a data set acquired with a 12-km long spread with 100 m spacing of shots and receivers, of the Neoproterozoic Marwar basin in the NW Indian shield. The approach is shown to be quite successful for delineating the shallow refractor depths, steep dips and velocities, even in the absence of regular reverse refraction profiles. The study reveals two-layered sedimentary formations, Malani volcanics and a complicated basement configuration of the Marwar basin, and provides a measure of resolution and uncertainty of the estimated model parameters. A seismic section of the near-trace gather is found to be qualitatively consistent with the derived structural features of the basin. The relative highs and lows, observed in the Bouguer gravity profile, further corroborate the derived velocity model. The present approach can be especially useful in offshore areas and elsewhere, where the single-ended multifold seismic profiles are the only available data sets.