The fine structure of the normal incidence synthetic seismogram

Summary . The spectral function of a perfectly elastic, horizontally stratified medium has been demonstrated previously to provide an attractive formulation to describe the properties of the one-dimensional synthetic seismogram (Robinson & Treitel). Here we examine the mathematical framework of the Model in still greater detail. Knowledge of this fine structure of the synthetic seismogram leads to the solution of two particular seismic inverse problems. First, we consider a layered medium with an arbitrary surface reflection coefficient c _o, where | c _o|<1, and which contains an impulsive source immediately above the surface. Given the corresponding synthetic seismogram, we develop an inverse, or backward recursion formalism which recovers the entire series of original reflection coefficients. Second, we consider a similar problem for an impulsive source located just below the surface. Both inversion procedures constitute a continuation of the work of Goupillaud and of Sherwood & Trorey, and represent a generalization of the classical technique originally proposed by Kunetz which, however, only holds for the marine case, c_o =±1. The present approach is not so constrained and thereby becomes applicable to land seismograms as well.
If products of third or higher order in the reflection coefficients can be neglected, significant simplifications arise in the theory. In that event the usual representation of the synthetic seismogram as a ratio of two polynomials in the complex variable z becomes particularly revealing. The numerator polynomial is then approximately equal to the z transform of the reflection coefficient series, while the denominator polynomial is approximately equal to the z transform of the autocorrelation of these reflection coefficients. The resulting simplified theory affords important computational savings in the appropriate backward recursion algorithms.     

Stable Iterative Methods for the Inversion of Geophysical Data

Summary. Interpretation of earth electrical measurements can often be assisted by inversion, which is a non-linear model-fitting problem in these cases. Iterative methods are normally used, and the solution is defined by' best fit'in the sense of generalized least-squares.
The inverse problems we describe are ill-posed. That is, small changes in the data can lead to large changes in both the solution and in the iterative process that finds the solution. Through an analysis of the problem, based on local linearization, we define a class of methods that stabilize the iteration, and provide a robust solution. These methods are seen as generalizations of the well-known Singular Value Truncation and Marquardt Methods of iterative inversion.
Here, and in a companion paper, we give examples illustrating the successful application of the method to ill-posed problems relating to the resistivity of the Earth.     

Inversion of time domain three-dimensional electromagnetic data

We present a general formulation for inverting time domain electromagnetic data to recover a 3-D distribution of electrical conductivity. The forward problem is solved using finite volume methods in the spatial domain and an implicit method (Backward Euler) in the time domain. A modified Gauss–Newton strategy is employed to solve the inverse problem. The modifications include the use of a quasi-Newton method to generate a pre-conditioner for the perturbed system, and implementing an iterative Tikhonov approach in the solution to the inverse problem. In addition, we show how the size of the inverse problem can be reduced through a corrective source procedure. The same procedure can correct for discretization errors that inevidably arise. We also show how the inverse problem can be efficiently carried out even when the decay time for the conductor is significantly larger than the repetition time of the transmitter wave form. This requires a second processor to carry an additional forward modelling. Our inversion algorithm is general and is applicable for any electromagnetic field  ( E , H , d B / dt )  measured in the air, on the ground, or in boreholes, and from an arbitrary grounded or ungrounded source. Three synthetic examples illustrate the basic functionality of the algorithm, and a result from a field example shows applicability in a larger-scale field example.     

Stability considerations for one-dimensional inverse problems

Summary. We investigate the issues of stability and conditioning for the one-dimensional seismic inverse problem. We show that these issues are distinct; i.e. that numerically stable implementations of solutions to the inverse problem will not give accurate results if the problem is ill-conditioned. In addition, we identify the factors which determine the condition of the inverse problem. These are the total variation of the acoustic impedance profile being sought and the accuracy of the low-frequency content of the reflection data. We illustrate these results on implementations of two numerically stable algorithms for the inverse problem, one of which has a reputation for being unstable. A comparison shows nearly identical results for the two methods on noise-contaminated and frequency band limited reflection data. In fact, we conjecture that all of the well-known 'layer-stripping' inverse scattering methods share the same mathematical stability characteristics. On the other hand, we also show that ill-conditioning can lead to failure of such algorithms, through amplification of error due either to inaccurate data or to discretization or roundoff. Finally, we observe that appropriate smoothing of the seismic data for an ill-conditioned inverse problem (high-variation impedance profile) can cause the problem to become well-conditioned (lower-variation profile). As is typical with regularizations, the price paid for the newly-acquired ability to solve the problem is a loss of accuracy in the solution.     

A general inverse method for modelling extensional sedimentary basins

A two-dimensional inverse model for extracting the spatial and temporal variation of strain rate from extensional sedimentary basins is presented and applied. This model is a generalization of a one-dimensional algorithm which minimizes the misfit between predicted and observed patterns of basin subsidence. Our calculations include the effects of two-dimensional conduction and advection of heat as well as flexural rigidity. More importantly, we make no prior assumptions about the duration, number or intensity of rifting periods. Instead, the distribution of strain rate is permitted to vary smoothly through space and time until the subsidence misfit has been minimized. We have applied this inversion algorithm to extensional sedimentary basins in a variety of geological settings. Basin stratigraphy can be accurately fitted and the resultant spatiotemporal distributions of strain rate are corroborated by independent information about the number and duration of rifting episodes. Perhaps surprisingly, the smallest misfits are achieved with flexural rigidities close to zero. Spatiotemporal strain rate distributions will help to constrain the dynamical evolution of thinning continental lithosphere. The strain rate pattern governs the heat-flow history and so two-dimensional inversion can be used to construct accurate maturation models. Finally, our inversion algorithm is a stepping stone towards a generalized three-dimensional implementation.     

Moment tensor inversion of complex earthquakes

Summary. Moment tensor inversion methods can be applied with success in the determination of source properties of simple earthquakes. However, these methods utilize the assumption of a point source, which is inadequate for modelling many complicated, shallow earthquakes. For complex earthquakes, an inversion using finite faulting models is desirable but the number of parameters involved requires that a good starting model be found or that independent constraints be placed on some of the parameters. A method is presented for low-pass filtering both the data and Green's functions, passing only signals with wavelengths greater than the dimension of the entire fault. The filter tends to smooth complications in the waveforms and allows application of the point source moment tensor inversion. This method is applied to body waves from the 1978 Thessaloniki, Greece, earthquake, the 1971 San Fernando earthquake and to a multiple-point source synthetic model of the San Fernando event. For the Thessaloniki event, although a multiple-source mechanism has been suggested, inversion results before and after filtering were essentially identical, indicating that a point source mechanism is sufficient in modelling the long-period, teleseismic body waves. In the case of the San Fernando earthquake, the point source Green's functions were incapable of simultaneously modelling the P - and SH -waves. Inversion of P -waves alone resulted in extreme parameter resolution problems, but allowed constraint in one axis of the moment tensor and suggested an overall source time function. Inversion of a synthetic San Fernando data set yielded similar results, but allowed an investigation of the shortcomings of the method under controlled circumstances. Although the results may require substantial interpretation, the method presented represents a simple first step in the analysis of complex earthquakes.     

Automatic differentiation in geophysical inverse problems

Automatic differentiation (AD) is the technique whereby output variables of a computer code evaluating any complicated function (e.g. the solution to a differential equation) can be differentiated with respect to the input variables. Often AD tools take the form of source to source translators and produce computer code without the need for deriving and hand coding of explicit mathematical formulae by the user. The power of AD lies in the fact that it combines the generality of finite difference techniques and the accuracy and efficiency of analytical derivatives, while at the same time eliminating 'human' coding errors. It also provides the possibility of accurate, efficient derivative calculation from complex 'forward' codes where no analytical derivatives are possible and finite difference techniques are too cumbersome. AD is already having a major impact in areas such as optimization, meteorology and oceanography. Similarly it has considerable potential for use in non-linear inverse problems in geophysics where linearization is desirable, or for sensitivity analysis of large numerical simulation codes, for example, wave propagation and geodynamic modelling. At present, however, AD tools appear to be little used in the geosciences. Here we report on experiments using a state of the art AD tool to perform source to source code translation in a range of geoscience problems. These include calculating derivatives for Gibbs free energy minimization, seismic receiver function inversion, and seismic ray tracing. Issues of accuracy and efficiency are discussed.     

Elastic full-waveform inversion for earthquake source parameters

2-D full-waveform inversion of double-couple earthquake sources is implemented. Temporally and spatially extended sources are represented by superposition of double-couples. The source parameters solved for are the spatial location, origin time, amplitude and orientation of each double-couple. The velocity and density distribution and source time function are assumed to be known a priori but may be arbitrarily complicated. The non-linear inverse problem is solved by iterative linear approximation. The Jacobian matrix elements for source depth and rupture angle are computed by wavefield extrapolation forward in time, while those for origin time and amplitude are computed analytically. A smoothing technique that results in faster convergence and avoids local minima associated with cycle skipping is applied at each iteration. A spatial sampling interval, between discrete sources, of one-quarter wavelength of the dominant shear wave is optimal for inversion if high uniqueness of the result is desired. The presence of a fault is inferred from the spatial continuity of the rupture solution, rather than being imposed a priori. The method is illustrated by successful application to three synthetic source models: a single double-couple, a single extended rupture and a double extended rupture. The resolutions of the source depth and origin time are higher, and their posterior covariances are lower than those of the amplitude and rupture angle at each source point. Source depth, origin time and amplitude are primarily determined by the data; the rupture angle is more strongly influenced by the a priori information.     

Displacements and stresses due to a single force in a half-space in welded contact with another half-space

Closed-form analytical expressions for the displacements and stresses induced by a single force of arbitrary orientation located in an elastic half-space in welded contact with another elastic half-space are obtained. These expressions are valid for arbitrary values of the Poisson's ratio and for arbitrary source and observer locations. The final results are given in a form that makes numerical computation straightforward and accurate.     

A search for the ratio between gravity variation and vertical displacement due to a surface load

A data space approach to magnetotelluric (MT) inversion reduces the size of the system of equations that must be solved from M × M , as required for a model space approach, to only N × N , where M is the number of model parameter and N is the number of data. This reduction makes 3-D MT inversion on a personal computer possible for modest values of M and N . However, the need to store the N × M sensitivity matrix J remains a serious limitation. Here, we consider application of conjugate gradient (CG) methods to solve the system of data space Gauss–Newton equations. With this approach J is not explicitly formed and stored, but instead the product of J with an arbitrary vector is computed by solving one forward problem. As a test of this data space conjugate gradient (DCG) algorithm, we consider the 2-D MT inverse problem. Computational efficiency is assessed and compared to the data space Occam's (DASOCC) inversion by counting the number of forward modelling calls. Experiments with synthetic data show that although DCG requires significantly less memory, it generally requires more forward problem solutions than a scheme such as DASOCC, which is based on a full computation of J .     

P-wave traveltime and polarization tomography of VSP data

The presence of anisotropy requires that tomographic methods be generalized to account for anisotropy. This generalization allows geological structure to be correctly imaged and allows the anisotropic parameters to be estimated. Use of isotropic inversion for imaging anisotropic structures gives systematic trends in the traveltime and polarization residuals. However, due to the limited directional coverage, the traveltimes along may not be sufficient to study the anisotropic properties of the structure. Polarizations can provide independent information on the structure. Traveltime and polarization inversion are applied to synthetic examples simulating VSP experiments. Transverse isotropy and 1-D structure are assumed. Plots of traveltime and polarization residuals are an important tool to detect the anomalies due to the presence of anisotropy. For receivers located in anisotropic layers, polarization residuals display consistent anomalies of several degrees. The synthetic examples show that even the simple 1-D problem is difficult, when using direct arrivals only. Large a posteriori errors in anisotropic parameters are obtained by traveltime inversion in layers where available incidence angles are less than 45°. Resolution of the tomographic image of VSP data is greatly improved by a combination of traveltime and polarization information. In order to obtain accurate inversion results, the measurement error of polarization data should be kept to within a few degrees.     

The use of the Born approximation in seismic scattering problems

Summary. In view of recent work on seismic scattering by small-scale heterogeneities in the Earth, which has been based on single-scattering perturbation theory (that is, the Born approximation), we attempt to define the region within which this approximation may be regarded as reasonably accurate. Comparison with the exact solution for scattering by an embedded sphere shows that the inequalities obtained, governing the ranges of the parameters of the problem, are appropriate.
However, application of these constraints on the parameters imply that, in almost all cases, application of the Born approximation to upper and lower mantle scattering is probably invalid.     

Two efficient algorithms for iterative linearized inversion of seismic waveform data

We present two equivalent algorithms for iterative linearized waveform inversion for 3-D Earth structure with respect to an arbitrary 3-D starting model; one is a matrix formulation, and the second is a wavefield formulation. Both algorithms require the computation of accurate synthetic seismograms, but neither requires that any particular method be used to compute the synthetics. The matrix formulation is equivalent to our previously published algorithm (Hara, Tsuboi & Geller 1991), but requires less than 10 per cent of the CPU time of the previous algorithm. The wavefield algorithm is equivalent to that of Tarantola (1986) and Mora (1987), but appears to be substantially more efficient.     

Inversion of seismic data using tomographical reconstruction techniques for investigations of laterally inhomogeneous media

Summary. For the determination of lateral velocity or absorption inhomogeneities, methods such as the generalized matrix inversion and its damped versions, for example the stochastic inverse, are usually applied in seismology to travel-time or amplitude anomalies. These methods are not appropriate for the solution of very extensive systems of equations. Reconstruction techniques as developed for computer tomography are suitable for operations with extremely large numbers of equations and unknown parameters. In this paper solutions obtained with the BPT (Back Projection Technique), ART (Algebraic Reconstruction Technique) and SIRT (Simultaneous Iterative Reconstruction Technique) are compared with those obtained from a damped version of the generalized inverse method. Data of 2-D model-seismic experiments are presented for demonstration.     

Surface wave tomography for azimuthal anisotropy in a strongly reduced parameter space

Large scale seismic anisotropy in the Earth's mantle is likely dynamically supported by the mantle's deformation; therefore, tomographic imaging of 3-D anisotropic mantle seismic velocity structure is an important tool to understand the dynamics of the mantle. While many previous studies have focused on special cases of symmetry of the elastic properties, it would be desirable for evaluation of dynamic models to allow more general axis orientation. In this study, we derive 3-D finite-frequency surface wave sensitivity kernels based on the Born approximation using a general expression for a hexagonal medium with an arbitrarily oriented symmetry axis. This results in kernels for two isotropic elastic coefficients, three coefficients that define the strength of anisotropy, and two angles that define the symmetry axis. The particular parametrization is chosen to allow for a physically meaningful method for reducing the number of parameters considered in an inversion, while allowing for straightforward integration with existing approaches for modelling body wave splitting intensity measurements. Example kernels calculated with this method reveal physical interpretations of how surface waveforms are affected by 3-D velocity perturbations, while also demonstrating the non-linearity of the problem as a function of symmetry axis orientation. The expressions are numerically validated using the spectral element method. While challenges remain in determining the best inversion scheme to appropriately handle the non-linearity, the approach derived here has great promise in allowing large scale models with resolution of both the strength and orientation of anisotropy.     

Joint three-dimensional inversion of coupled groundwater flow and heat transfer based on automatic differentiation: sensitivity calculation, verification, and synthetic examples

Inverse methods are useful tools not only for deriving estimates of unknown parameters of the subsurface, but also for appraisal of the thus obtained models. While not being neither the most general nor the most efficient methods, Bayesian inversion based on the calculation of the Jacobian of a given forward model can be used to evaluate many quantities useful in this process. The calculation of the Jacobian, however, is computationally expensive and, if done by divided differences, prone to truncation error. Here, automatic differentiation can be used to produce derivative code by source transformation of an existing forward model. We describe this process for a coupled fluid flow and heat transport finite difference code, which is used in a Bayesian inverse scheme to estimate thermal and hydraulic properties and boundary conditions form measured hydraulic potentials and temperatures. The resulting derivative code was validated by comparison to simple analytical solutions and divided differences. Synthetic examples from different flow regimes demonstrate the use of the inverse scheme, and its behaviour in different configurations.     

Inversion of Love-wave data for velocity and anelasticity using exact kernels

Summary. The general problem of inverting Love-wave dispersion and amplitude data to obtain a velocity and Q_s structure is considered. A formulation is used which incorporates attenuation into the Haskell-Thompson matrix method in an exact manner and thus retains the inherent non-linearity in the anelasticity. The resulting exact inversion kernels allow simultaneous inversion for velocity and intrinsic attenuation parameters. The method is applied to synthetic data which allows a comparison to be made with inexact kernels. The results indicate that the use of inexact kernels may introduce spurious oscillations into the Q_s structure and that a simultaneous inversion can be more stable than inverting for velocity alone.     

Effects of topography and crustal heterogeneities on the source estimation of LP event at Kilauea volcano

The main goal of this study is to improve the modelling of the source mechanism associated with the generation of long period (LP) signals in volcanic areas. Our intent is to evaluate the effects that detailed structural features of the volcanic models play in the generation of LP signal and the consequent retrieval of LP source characteristics. In particular, effects associated with the presence of topography and crustal heterogeneities are here studied in detail. We focus our study on a LP event observed at Kilauea volcano, Hawaii, in 2001 May. A detailed analysis of this event and its source modelling is accompanied by a set of synthetic tests, which aim to evaluate the effects of topography and the presence of low velocity shallow layers in the source region. The forward problem of Green's function generation is solved numerically following a pseudo-spectral approach, assuming different 3-D models. The inversion is done in the frequency domain and the resulting source mechanism is represented by the sum of two time-dependent terms: a full moment tensor and a single force. Synthetic tests show how characteristic velocity structures, associated with shallow sources, may be partially responsible for the generation of the observed long-lasting ringing waveforms. When applying the inversion technique to Kilauea LP data set, inversions carried out for different crustal models led to very similar source geometries, indicating a subhorizontal cracks. On the other hand, the source time function and its duration are significantly different for different models. These results support the indication of a strong influence of crustal layering on the generation of the LP signal, while the assumption of homogeneous velocity model may bring to misleading results.     

One-way representations of seismic data

A general one-way representation of seismic data can be obtained by substituting a Green's one-way wavefield matrix into a reciprocity theorem of the convolution type for one-way wavefields. From this general one-way representation, several special cases can be derived.
By introducing a Green's one-way wavefield matrix for primaries , a generalized Bremmer series representation is obtained. Terminating this series after the first-order term yields a primary representation of seismic reflection data. According to this representation, primary seismic reflection data are proportional to a reflection operator, 'modified' by primary propagators for downgoing and upgoing waves. For seismic imaging, these propagators need to be inverted. Stable inverse primary propagators can easily be obtained from a one-way reciprocity theorem of the correlation type.
By introducing a Green's one-way wavefield matrix for generalized primaries , an alternative representation is obtained in which multiple scattering is organized quite differently (in comparison with the generalized Bremmer series representation). According to the generalized primary representation, full seismic reflection data are proportional to a reflection operator, 'modified' by generalized primary propagators for downgoing and upgoing waves. Internal multiple scattering is fully included in the generalized primary propagators {either via a series expansion or in a parametrized way). Stable inverse generalized primary propagators can be obtained from the one-way reciprocity theorem of the correlation type. These inverse propagators are the nucleus for seismic imaging techniques that take the angle-dependent dispersion effects due to fine-layering into account.     

Source investigation of a small event using empirical Green's functions and simulated annealing

We propose a two-step inversion of three-component seismograms that (1) recovers the far-field source time function at each station and (2) estimates the distribution of co-seismic slip on the fault plane for small earthquakes (magnitude 3 to 4). The empirical Green's function (EGF) method consists of finding a small earthquake located near the one we wish to study and then performing a deconvolution to remove the path, site, and instrumental effects from the main-event signal.
The deconvolution between the two earthquakes is an unstable procedure: we have therefore developed a simulated annealing technique to recover a stable and positive source time function (STF) in the time domain at each station with an estimation of uncertainties. Given a good azimuthal coverage, we can obtain information on the directivity effect as well as on the rupture process. We propose an inversion method by simulated annealing using the STF to recover the distribution of slip on the fault plane with a constant rupture-velocity model. This method permits estimation of physical quantities on the fault plane, as well as possible identification of the real fault plane.
We apply this two-step procedure for an event of magnitude 3 recorded in the Gulf of Corinth in August 1991. A nearby event of magnitude 2 provides us with empirical Green's functions for each station. We estimate an active fault area of 0.02 to 0.15 km² and deduce a stress-drop value of 1 to 30 bar and an average slip of 0.1 to 1.6 cm. The selected fault of the main event is in good agreement with the existence of a detachment surface inferred from the tectonics of this half-graben.