Search

This article explains the features of differential data that make them attractive, their shortcomings, and the situations for which they are best suited. The use of differential data is ubiquitous in the seismological community, in which they are used to determine earthquake locations via the double-difference method and the Earth’s velocity structure via geotomography; furthermore, they have important applications in other areas of geophysics, as well. A common assumption is that differential data are uncorrelated and have uniform variance. We show that this assumption is well justified when the original, undifferenced data covary with each other according to a two-sided exponential function. It is not well justified when they covary according to a Gaussian function. Differences of exponentially correlated data are approximately uncorrelated with uniform variance when they are regularly spaced in distance. However, when they are irregularly spaced, they are uncorrelated with a nonuniform variance that scales with the spacing of the data. When differential data are computed by taking differences of the original, undifferenced data, model parameters estimated using ordinary least squares applied to the differential data are almost exactly equal to those estimated using weighed least squares applied to the original, undifferenced data (with the weights given by the inverse covariance matrix). A better solution only results when the differential data are directly estimated and their variance is smaller than is implied by differencing the original data. Differential data may be appropriate for global seismic travel-time data because the covariance of errors in predicted travel times may have a covariance close to a two-sided exponential, on account of the upper mantle being close to a Von Karman medium with exponent κ≪12.

We review Gaussian process regression (GPR) and analyze it in the context of Inverse Theory—the collection of techniques used in geophysics (among other fields) to understand the structure of data analysis problems and the quality of their solutions. By viewing GPR as a special case of generalized least squares (least squares with prior information), we derive expressions for a variety of standard Inverse Theory quantities, including the data and model resolution matrices, the importance (influence) vector, and the gradient of the solution with respect to a parameter. We study the impulse response in the one-dimensional continuum limit and provide formulas for its area and width. Finally, we demonstrate how the importance vector can be used to design an optimum GPR experiment, through a process we call importance winnowing.

Generalized Least Squares (least squares with prior information) requires the correct assignment of two prior covariance matrices: one associated with the uncertainty of measurements; the other with the uncertainty of prior information. These assignments often are very subjective, especially when correlations among data or among prior information are believed to occur. However, in cases in which the general form of these matrices can be anticipated up to a set of poorly-known parameters, the data and prior information may be used to better-determine (or “tune”) the parameters in a manner that is faithful to the underlying Bayesian foundation of GLS. We identify an objective function, the minimization of which leads to the best-estimate of the parameters and provide explicit and computationally-efficient formula for calculating the derivatives needed to implement the minimization with a gradient descent method. Furthermore, the problem is organized so that the minimization need be performed only over the space of covariance parameters, and not over the combined space of model and covariance parameters. We show that the use of trade-off curves to select the relative weight given to observations and prior information is not a form of tuning, because it does not, in general maximize the posterior probability of the model parameters, and can lead to a different weighting than the procedure described here. We also provide several examples that demonstrate the viability, and discuss both the advantages and limitations of the method.

The non-double-couple (non-DC) components of the moment tensor provide insight into the earthquake processes and anisotropy of the near-source region. We investigate the behavior of the isotropic (ISO) and compensated linear vector dipole (CLVD) components of the moment tensor for shear faulting in a transversely ISO medium with an arbitrarily oriented symmetry axis. Analytic formulas for ISO and CLVD depend on the orientation of the fault relative to the anisotropy symmetry axis as well as three anisotropic parameters, which describe deviations of the medium from isotropy. Numerical experiments are presented for the preliminary reference Earth model. Both ISO and CLVD components are zero when the axis of symmetry is within the fault plane or the auxiliary plane. For any orientation in which the ISO component is zero, the CLVD component is also zero, but the opposite is not always true (e.g., for strong anisotropy). The relative signs of the non-DC components of neighboring earthquakes may help distinguish source processes from source-region anisotropy. We prove that an inversion of ISO and CLVD components of a set of earthquakes with different focal mechanisms can uniquely determine the orientation and strength of anisotropy. This study highlights the importance of the ISO component for constraining deep slab anisotropy and demonstrates that it cannot be neglected.

We solve an idealized version of the moment tensor—anisotropy inverse problem. The medium in which faulting occurs has general anisotropy described by 21 parameters. The data are a set of moment tensors for all possible configurations of a unit fault; that is, all combinations of strike, dip and rake (a very idealized assumption, for tectonic constraints may limit available configurations). We first ask whether the anisotropic parameters can be inferred uniquely from the moment tensors and arrive at the answer, “almost”. The data constrain only 20 linear combinations of the 21 parameters; the isotropic Lame parameter $$\lambda$$ cannot be determined. We provide an analytic solution algorithm, valid for both weak and strong anisotropy, that uses the data and a prior value of $$\lambda$$ . We then ask whether measurements of the trace of the moment tensor and its smallest eigenvalue (proxies for its explosive and compensated linear vector dipole components, respectively) can determine uniquely the anisotropy. Here our analysis is limited to the weak anisotropy case. We find that these data are “almost” sufficient. They constrain only 19 linear combinations of the 21 parameters; the isotropic Lame parameters $$\lambda$$ and $$\mu$$ cannot be determined. We provide a semi-analytic algorithm for solving the inverse problem, given prior values of $$\lambda$$ and $$\mu$$ . Numerical tests demonstrate that both solution methods work and can be used to assess the sensitivity of the solution to erroneous prior information. However, because the methods require data for all fault configurations, their application to the sparse data available for most seismic source regions is limited; they complement but do not supplant existing inversion techniques.

The goal of tomography is to reconstruct a spatially-varying image function s(x,m), where x is position and m is a finite-length vector of parameters. Many reconstruction methods minimize the total L2 error E ≡ eTe, where individual errors ei quantify misfit between predictions and observations, to quantify goodness of fit. So-called adjoint state methods allow the gradient ∂E/∂mi to be computed extremely efficiently from an adjoint field, facilitating image reconstruction by gradient-descent methods. We examine the structure of the differential equation for the adjoint field under the ray approximation and find that it has the same form as the transport equation, whose solution involves the well-known geometrical spreading function R Consequently, as R is routinely tabulated as part of a ray calculation, no extra work is needed to compute the adjoint field, permitting a rapid calculation of the gradient ∂E/∂mi.

SUMMARY We consider the case where the ‘solution’ to an inverse problem is an ensemble (e.g. drawn from the conditional probability density function $p( {{{\bf m}}|{{{\bf d}}^{obs}}} )$ of M model parameters ${{\bf m}}$ given observed data ${{{\bf d}}^{obs}}$). Here we presume that the ${{\bf m}}$s have a natural ordering, say in position x, so that ‘resolution’ means the ability of the inverse problem to distinguish physically adjacent model parameters. The trade-off curve for resolution and variance is constructed using the following steps: (1) the single solution ${{{\bf m}}^{est}}$ and its covariance ${{{\bf C}}_m}$ are estimated as the ensemble mean and covariance; (2) the eigenvalue decomposition ${{{\bf C}}_m} = {\rm{\ }}{{{\bf V} {\boldsymbol \Lambda} }}{{{\bf V}}^{\rm{T}}}$ is computed and the submatrix ${{{\boldsymbol \Lambda }}^{( N )}}$ of the N smallest eigenvalues, and submatrix ${{{\bf V}}^{( N )}}$of the N corresponding eigenvectors, are formed; (3) the equation ${{{\boldsymbol \mu }}^{( N )}} = {{{\boldsymbol \Phi }}^{( N )}}\ {{\bf m}}$ with ${{{\boldsymbol \mu }}^{( N )}} = [ {{{{\bf V}}^{( N )}}} ]{\boldsymbol{\ }}{{{\bf m}}^{est}}$and${{{\boldsymbol \Phi }}^{( N )}} = {[ {{{{\bf V}}^{( N )}}} ]^{\rm{T}}}$ is formed, as is its covariance ${{\bf C}}_\mu ^{( N )} = {{{\boldsymbol \Lambda }}^{( N )}}{\boldsymbol{\ }}$; (4) the equation is solved to yield a localized average ${\langle {{\bf m}} \rangle ^{( N )}} = \ {{{\boldsymbol \Phi }}^{ - g}}{{{\boldsymbol \mu }}^{( N )}}$, where ${{{\boldsymbol \Phi }}^{ - g}}$ is either the minimum length or Backus–Gilbert generalized inverse of ${{\boldsymbol \Phi }}$; (5) the resolution and covariance are computed as ${{{\bf R}}^{( N )}} = {{{\boldsymbol \Phi }}^{ - g}}{\boldsymbol{\ }}{{{\boldsymbol \Phi }}^{( N )}}$ and ${{\bf C}}_m^{( N )} = {{{\boldsymbol \Phi }}^{ - g}}{\boldsymbol{\ }}{{\bf C}}_\mu ^{( N )}{( {{{{\boldsymbol \Phi }}^{ - g}}} )^{\rm{T}}}$; (6) the spread ${K^{( N )}}$ of resolution and size ${J^{( N )}}\ $of covariance are computed using either the Dirichlet or Backus–Gilbert measures and (7) the process is repeated for $1 \le N \le M$ to build up the trade-off curve $K( J )$. We show that, in the Dirichlet case, ${K^{( N )}} = \ M - N$ and ${J^{( N )}} = \ {\rm{tr}}( {{{{\boldsymbol \Lambda }}^{( N )}}} )$. We also consider the case where the model parameters correspond to spline coefficients and a sequence ${y_i}( {{{\bf m}},{x_i}} )$ derived from these coefficients possesses natural ordering. Layered models are an example of such a parametrization. We construct the trade-off curve for ${{\bf y}}$ by converting each member of the ensemble from ${{\bf m}}$ to ${{\bf y}}$ and applying the above procedure to them. We demonstrate the method by applying it to several simple examples.

We consider a real-valued function on a plane of the form m(x,y,θ)=A(x,y)+Bc(x,y)cos(2θ)+Bs(x,y)sin(2θ)+Cc(x,y)cos(4θ)Cs(x,y)sin(4θ) that models anisotropic acoustic slowness (reciprocal velocity) perturbations. This “slowness function” depends on Cartesian coordinates and polar angle θ. The five anisotropic “component functions” A (x,y), Bc(x,y), Bs(x,y), Cc(x,y) and Cs(x,y) are assumed to be real-valued Schwartz functions. The “travel time” function d(u, θ) models the travel time perturbations on an indefinitely long straight-line observation path, where the line is parameterized by perpendicular distance u from the origin and polar angle θ; it is the Radon transform of m ( x, y, θ). We show that: 1) an A can always be found with the same d(u, θ) as an arbitrary (Bc,Bs) and/or an arbitrary (Cc,Cs) ; 2) a (Bc,Bs) can always be found with the same d(u, θ) as an arbitrary A, and furthermore, infinite families of them exist; 3) a (Cc,Cs) can always be found with the same d(u, θ) as an arbitrary A, and furthermore, infinite families of them exist; 4) a (Bc,Bs) can always be found with the same d(u, θ) as an arbitrary (Cc,Cs) , and vice versa; and furthermore, infinite families of them exist; and 5) given an arbitrary isotropic reference slowness function m0(x,y), “null coefficients” (Bc,Bs) can be constructed for which d(u, θ) is identically zero (and similarly for Cc,Cs ). We provide explicit methods of constructing each of these “equivalent functions”.

We prove that the problem of inverting Rayleigh wave phase velocity functions $$c\left( k \right)$$ , where $$k$$ is wavenumber, for density $$\rho \left( z \right)$$ , rigidity $$\mu \left( z \right)$$ and Lame parameter $$\lambda \left( z \right)$$ , where $$z$$ is depth, is fully non-unique, at least in the highly-idealized case where the base Earth model is an isotropic half space. The model functions completely trade off. This is one special case of a common inversion scenario in which one seeks to determine several model functions from a single data function. We explore the circumstances under which this broad class of problems is unique, starting with very simple scenarios, building up to the somewhat more complicated (and common) case where data and model functions are related by convolutions, and then finally, to scale-independent problems (which include the Rayleigh wave problem). The idealized cases that we examine analytically provide insight into the kinds of nonuniqueness that are inherent in the much more complicated problems encountered in modern geophysical imaging (though they do not necessarily provide methods for solving those problems). We also define what is meant by a Backus and Gilbert resolution kernel in this kind of inversion and show under what circumstances a unique localized average of a single model function can be constructed.

Previously-published tomographic images show low-velocity anomalies (shear velocity up to 5% slow relative to the AK135 Earth model) in the asthenosphere beneath the southeastern margin of North America, including the Northern Appalachian (New England), the Central Appalachian (Virginia), and the Northern Gulf (Texas-Louisiana coast) Anomalies. Travel time anomaly ratios and attenuation indicate that they are thermal anomalies associated with mantle upwelling. We use teleseismic P and S wave differential travel times, as determined by cross-correlation of USArray data, and ratios of compressional-to-shear wave travel time fluctuations to characterize the South Coastal Anomaly (SCA), a fourth, weaker anomaly (about half as strong as these others) that stretches from Georgia to Virginia. P-wave tomography indicates that the SCA is about 800 km long in its north-south dimension and at least 300 km wide in its east-west dimension (its eastern edge is not imaged) and it is strongest (∆VP ≈ − 0.1 km/s) in the 100–250 km depth range. Its western edge strikes north-south, parallel to the edge of the Laurentian Craton (LC). The maximum P and S wave travel time anomalies of the SCA are 1.09 ± 0.05 (95%) s and 4.20 ± 0.10 (95%) s, respectively, relative to the LC. The ratio of P-to-S wave differential travel time anomalies is 3.84 ± 0.10 (95%), a value close to the value of ~3.3 predicted for a thermal anomaly and very different from the value of ~1.8 predicted for a compositional anomaly. The compressional velocity is up to 0.41 km/s slower than the LC, which corresponds to a maximum temperature difference of about 700 °C, using a sensitivity of 1713 K-s/km. In the shallow asthenosphere, the SCA appears distinct from the Central Appalachian Anomaly, being weaker in magnitude and more tabular in shape, but the two anomalies may merge below ~250 km depth. Our findings support the notion of the vigorous edge convection along much of the eastern continental margin.

Seismic measurements are used in a detailed investigation of a region of extremely low asthenospheric seismic velocities along the US Gulf Coast, first imaged in continental-scale geotomography, which we term the Northern Gulf Anomaly (NGA). Differential P- and S-wave arrival times from teleseisms at a variety of back-azimuths, observed on EarthScope Transportable Array stations near the US Gulf Coast, demonstrate that asthenospheric seismic velocities are 8-10% lower than the neighboring craton, and define the spatial extent and character of the anomaly. Travel time anomalies are calculated relative to the AK135 earth model and corrected to account for the effect of the kilometers-thick sedimentary cover in the region. The NGA is most intense at the southernmost coast of Louisiana and East Texas (with an eastern edge at 89°W) and smoothly tapers away in a triangular wedge that extends inland as far as 300 km. It has sharper edges and a smaller areal extent (by ∼ 50 % ) than previously-published geotomography has indicated. Both the magnitude and ratio of delays indicate that the NGA has a thermal origin, which may represent past or present-day small-scale convective upwelling near the southeastern edge of the North American continent. The NGA suggests that large-scale but poorly understood asthenospheric processes are at work beneath the US Gulf Coast, notwithstanding this region's reputation as an aseismic, passively-subsiding continental margin.

We present an improved method for estimating the frequency‐dependent phase velocity of ambient noise cross spectra by waveform fitting of Aki’s formula (Aki, 1957). The key part of the procedure is the use of a grid search to generate an initial estimate of the phase velocity that matches the observed cross spectrum without any cycle slips. A generalized least‐squares procedure is then used to refine this initial estimate. The grid search has only three frequency nodes, at the left end, center, and right end of the frequency band of interest, but parameterizes velocity finely, with 40 increments at each of these nodes. Trial phase velocity curves are generated by linear interpolation between nodes, and the curve that best predicts the observed cross spectrum is chosen as the initial estimate. The generalized least‐squares refinement implements two types of prior information: the phase velocity is close to a smoothed version of the initial estimate, and the phase velocity varies smoothly with frequency. In addition to estimating the frequency‐dependent phase velocity function, the method provides its variance and resolution. The method is robust to low signal‐to‐noise levels and so may prove especially attractive in cases in which short array‐deployment times limit the quality of cross‐spectral estimates. It also provides useful estimates at short interstation distance, when the cross spectrum has few zero crossings.

We investigate whether 2D anisotropic travel‐time tomography can uniquely determine both the spatially varying isotropic and anisotropic components of the seismic velocity field. This issue was first studied by Mochizuki (1997) for the special case of Radon’s problem (tomography with infinitely long rays), who found it to be nonunique. Our analysis extends this result to all array geometries and demonstrates that all such tomographic inversions are nonunique. Any travel‐time dataset can be fit by a model that is either purely isotropic, purely anisotropic, or some combination of the two. However, a pair of purely isotropic and purely anisotropic velocity models that predict the same travel times are very different in other respects, including spatial scale. Thus, prior information can be used to select among equivalent solutions to achieve a unique solution embodying a given set of prior expectations about model properties. We extend the notion of a resolution test, used in traditional isotropic tomography, to the anisotropic case. Our equivalent heterogeneity analysis focuses on the anisotropic heterogeneity equivalent to a point isotropic heterogeneity, and vice versa. We demonstrate that it provides insights into the structure of an anisotropic tomography problem that facilitates both the selection of appropriate prior information and the interpretation of results. We recommend that it be routinely applied to all surface‐wave inversions where the presence of anisotropy is suspected, including those based on noise correlation.

Seismic anisotropy beneath the Eastern Lau Spreading Center (ELSC) is investigated using both Rayleigh waves and shear waves, using data from the 2009–2010 ELSC ocean bottom seismograph experiment. Phase velocities of Rayleigh waves determined by ambient noise cross correlation are inverted for azimuthally anisotropic phase velocity maps. Splitting of S waves from five intermediate and deep focus earthquakes was determined by waveform analysis. Taken together, Rayleigh wave and S wave data indicate that significant (~2%) anisotropy extends to at least 300 km depth. Both data sets indicate a fast direction aligned within a few degrees of the N10°E striking ELSC and somewhat oblique to the N25°E strike of the neighboring volcanic arc. We therefore describe the fast direction as spreading perpendicular, not convergence perpendicular and interpret it as due to ridge-parallel flow of the asthenosphere. However, the region arcward (east) of the ELSC has the stronger anisotropy, suggesting that the strongest flow gradients may occur in the region between the spreading center and the arc, in contrast to being centered beneath the ELSC. Fluids released from the underlying plate may produce anisotropic hydrous materials, but more importantly lower the viscosity, thus enhancing along-strike flow. Both could contribute to an along-strike fast direction signature. Seafloor spreading diminishes south of the seismic array, ceasing altogether south of latitude 25°S (500 km south of the array center), a region dominated by much slower tectonic extension, suggesting that asthenosphere is inflowing from the north to accommodate the increase in asthenospheric volume associated with the seafloor spreading.

This chapter extends the discrete inverse theory that is the main focus of this book to continuous problems. Localized averages are shown to be particularly useful in continuous problems, and the Backus-Gilbert techniques for developing generalized inverses with localized resolution are extended to the continuous case. The conversion of continuous problems to discrete problems, accomplished earlier in the book in an ad hoc fashion, is examined in a more systematic way. Attention is then turned to problems in which the data are sufficiently numerous that they can be considered a continuous function. The tomography problem is one of these, and the simplest straight-line ray version, Radon's problem, is derived and solved. Attention is then turned to generalizing the quantities and solution techniques of discrete inverse to the continuous case. Here, the vector is replaced with the function and the matrix with the linear operator. A key quantity, analogous to the transpose of a matrix, is the adjoint of the linear operator. It is at the heart of so-called adjoint methods of deriving the Frechet derivative of data with respect to the model, the continuous analog of the data kernel. Several useful Frechet derivatives are derived, including one relevant to the case where the data and model function are related indirectly, through a differential equation.

This chapter uses geometrical concepts drawn from the theory of vector spaces to solve inverse problems. The model parameters are considered one point in a space of all possible model parameters, the data are considered one point in a space of all possible data, and the theory is considered a way of mapping a point in model space into a point in data space. The coordinate axes of these spaces can be rotated or transformed, without changing any of the fundamental relationships. This process suggests methods of transforming a complicated inverse problem into one that is in some sense simpler, either with a rotation that preserved length (e.g., a Householder rotation) or with one that changes the measure of length in some specifically tailored fashion. This process is shown to clarify ideas of minimum error and uniqueness and leads to the definition of the natural solution of an inverse problem (and a corresponding natural generalized inverse). Singular-value decomposition is developed as a way to identify the null subspaces associated with solutions that do not affect the prediction error and to construct the natural solution. Vector space ideas are also useful in understanding problems with inequality constraints and can be used to derive general properties of its solution, as embodied in the Kuhn-Tucker theorem. These considerations lead to several new types of inverse problem solutions, including nonnegative least squares.

This chapter follows up on issues associated with the nonuniqueness of some inverse problem solutions. Nonuniqueness is shown to be equivalent to the existence of null solutions to the inverse problem, solutions that do not have any effect on the prediction error. Some localized average of model parameters are shown to be unique, even when the underlying model parameters are nonunique. The properties of a unique localized average are enumerated, their advantages are explored, and examples are given. Nonunique localized averages are shown to be useful, if a priori information bounds the overall range of the underlying model parameters. The linear programming algorithm is used to find the bounds on a localized average, given the bounds on the underlying model parameters.

This chapter explores the consequence of data having non-Gaussian statistics. In these cases, the method of least squares is inappropriate, and different solution techniques must be applied. The exponential p.d.f. is especially useful, because it corresponds to the frequently encountered case where data outliers are common. A method of transforming the linear problem with exponential statistics into a linear programming problem is developed. It leads to a class of inverse problem solutions that minimizes some combination of prediction error and model error, as quantified under the L 1 norm. A similar transformation method for the L ∞ norm is also developed.

This chapter shows that inverse problems can be classified on the basis of the mathematical structure of the theory connecting the data and model parameters. The simplest theory is the explicit linear theory, which has the form of a matrix equation relating a vector of model parameters to a vector of data. The most complicated is the implicit nonlinear problem, in which the data and model parameters are mixed together in nonlinear relationships. A series of six exemplary inverse problems are discussed, starting with the simple notion of fitting a straight line to data exhibiting a linear trend, and culminating with acoustic tomography and X-ray imaging. In each case, the mathematical relationship between the data and model parameters is derived and classified. The meaning of the term solution to an inverse problem is debated; for depending upon the context and one's point of view, any of several rather different quantities can be considered solutions. Solutions can include estimates of model parameters, deterministic or probabilistic bounds, probability density functions, and estimates of localized averages.

Seismic attenuation exhibits strong geographic variability in northeastern North America, with the highest values associated with the previously-recognized Northern Appalachian Anomaly (NAA) in southern New England. The shear wave quality factor at 100 km depth is 14