ALBERT

Notes: A common example of a large-scale non-linear inverse problem is the inversion of seismic waveforms. Techniques used to solve this type of problem usually involve finding the minimum of some misfit function between observations and theoretical predictions. As the size of the problem increases, techniques requiring the inversion of large matrices become very cumbersome. Considerable storage and computational effort are required to perform the inversion and to avoid stability problems. Consequently methods which do not require any large-scale matrix inversion have proved to be very popular. Currently, descent type algorithms are in widespread use. Usually at each iteration a descent direction is derived from the gradient of the misfit function and an improvement is made to an existing model based on this, and perhaps previous descent directions.A common feature in nearly all geophysically relevant problems is the existence of separate parameter types in the inversion, i.e. unknowns of different dimension and character. However, this fundamental difference in parameter types is not reflected in the inversion algorithms used. Usually gradient methods either mix parameter types together and take little notice of the individual character or assume some knowledge of their relative importance within the inversion process.We propose a new strategy for the non-linear inversion of multi-offset reflection data. The paper is entirely theoretical and its aim is to show how a technique which has been applied in reflection tomography and to the inversion of arrival times for 3D structure, may be used in the waveform case. Specifically we show how to extend the algorithm presented by Tarantola to incorporate the subspace scheme. The proposed strategy involves no large-scale matrix inversion but pays particular attention to different parameter types in the inversion.We use the formulae of Tarantola to state the problem as one of optimization and derive the same descent vectors. The new technique splits the descent vector so that each part depends on a different parameter type, and proceeds to minimize the misfit function within the sub-space defined by these individual descent vectors. In this way, optimal use is made of the descent vector components, i.e. one finds the combination which produces the greatest reduction in the misfit function based on a local linearization of the problem within the subspace. This is not the case with other gradient methods. By solving a linearized problem in the chosen subspace, at each iteration one need only invert a small well-conditioned matrix (the projection of the full Hessian on to the subspace). The method is a hybrid between gradient and matrix inversion methods. The proposed algorithm requires the same gradient vectors to be determined as in the algorithm of Tarantola, although its primary aim is to make better use of those calculations in minimizing the objective function.

Notes: This paper shows how the performance of a fully non-linear earthquake location scheme can be improved by taking advantage of problem-specific information in the location procedure. The genetic algorithm is best viewed as a method of parameter space sampling that can be used for optimization problems. It has been applied successfully in regional and teleseismic earthquake location when the network geometry is favourable. However, on a series of test events with unfavourable network geometries the performance of the genetic algorithm is found to be poor.We introduce a method to separate the spatial and temporal parameters in such a way that problems related to the strong trade-off between depth and origin time are avoided. Our modified algorithm has been applied to several test events. Performance over the unmodified algorithm is improved substantially and the computational cost is reduced. The algorithm is better suited to the determination of hypocentral location whether using arrival times, array information (slowness and azimuth) or a combination of both.A second type of modification is introduced which exploits the weak correlation between the epicentral parameters and depth. This algorithm also improves performance over the standard genetic algorithm search, except in circumstances where the depth and epicentre are not weakly correlated, which occurs when the azimuthal coverage is very poor, or when azimuth and slowness information are incorporated. On a shallow nuclear explosion with only teleseismic P arrivals available, the algorithm consistently converged to a depth very close to the true depth, indicating superior depth estimation for shallow earthquake locations over the unmodified algorithm.

Notes: Many techniques for the solution of seismic-wave propagation problems depend on the representation of the seismic wavefield in terms of a linear combination of basis functions, as for example in Fourier or Gaussian Beam expansions. A common formal representation encompasses such methods when a preferred coordinate is isolated to track the propagation path. Different techniques can be classified by the dependence of the basis functions on this preferred coordinate. The common representation provides useful insight into the relation between apparently disparate methods and can guide the development of computational techniques. This common framework allows the development of generalized propagator methods and a compact formulation of reflection and transmission problems. A general perturbation approach can be used either to add heterogeneity to an existing structure or to restore features, such as coupling between P and S waves, which have been ignored in an approximate development.

Notes: Over the last three years, a major international effort has been made by the Sub-Commission on Earthquake Algorithms of the International Association of Seismology and the Physics of the Earth's Interior (IASPEI) to generate new global traveltime tables for seismic phases to update the tables of Jeffreys & Bullen (1940). The new tables are specifically designed for convenient computational use, with high-accuracy interpolation in both depth and range. The new iasp91 traveltime tables are derived from a radially stratified velocity model which has been constructed so that the times for the major seismic phases are consistent with the reported times for events in the catalogue of the International Seismological Centre (ISC) for the period 1964–1987. The baseline for the P-wave traveltimes in the iasp91 model has been adjusted to provide only a small bias in origin time for well-constrained events at the main nuclear testing sites around the world.For P-waves at teleseismic distances, the new tables are about 0.7s slower than the 1968 P-tables (Herrin 1968) and on average about 1.8–1.9 s faster than the Jeffreys & Bullen (1940) tables. For S-waves the teleseismic times lie between those of the JB tables and the results of Randall (1971).Because the times for all phases are derived from the same velocity model, there is complete consistency between the traveltimes for different phases at different focal depths. The calculation scheme adopted for the new iasp91 tables is that proposed by Buland & Chapman (1983). Tables of delay time as a function of slowness are stored for each traveltime branch, and interpolated using a specially designed tau spline which takes care of square-root singularities in the derivative of the traveltime curve at certain critical slownesses. With this representation, once the source depth is specified, it is straightforward to find the traveltime explicitly for a given epicentral distance. The computational cost is no higher than a conventional look-up table, but there is increased accuracy in constructing the traveltimes for a source at arbitrary depth. A further advantage over standard tables is that exactly the same procedure can be used for each phase. For a given source depth, it is therefore possible to generate very rapidly a comprehensive list of traveltimes and associated derivatives for the main seismic phases which could be observed at a given epicentral distance.

Notes: The effect of the free surface can be removed from three-component seismic recordings to recover the incident upgoing wavefield, if the slowness and azimuth of the current wavefront are known as a function of time. For a single three-component station it is usually possible to estimate an azimuth for an event from the first arriving P-waves, but slowness estimates are less reliable when more than one wavetype is presented in the seismic wavetrain. However, the free surface correction operators are generally slowly varying functions of slowness and so some error in slowness can be tolerated.Effective approximations for the removal of the free surface effects can be made for hard rock sites to cover slowness bands for the main regional phases Pn, Pg, Sn and Lg. By applying these operators in turn over group velocity windows appropriate to the particular phases, the relative amplitude of the P, SV and SH contributions to the wavefield can be estimated. Because the free surface amplification effects have been removed, the amplitudes can be compared directly and provide useful constraints on the radiation characteristics of the source. This procedure is therefore helpful for developing discrimination measures for different classes of sources.

Notes: Laterally varying interfaces cause coupling between wavenumbers so that seismograms in two-dimensionally layered media can be synthesized by means of ‘supermatrices’, which include the coupled contributions of all the wavenumbers. We introduce reflection and transmission ‘supermatrices’ in order to eliminate numerical problems arising from loss of precision for evanescent waves in the seismogram synthesis. An interface is assumed to be such that the reflected and transmitted wavefields; on its two sides can be represented as purely upgoing and downgoing waves, i.e. the Rayleigh ansatz is imposed. The computational demands of this method can be kept to a minimum by exploiting propagation invariants in the coupled wavenumber domain.The superior performance of this ‘invariant embedding’ approach when compared to propagator or finite difference schemes is illustrated by application to the response of sedimentary basins to excitation by an incident plane wave or a line force. The results are in good general agreement with the other methods, but show greater numerical stability and computational efficiency. In the case of a single interface the ‘invariant embedding’ procedure for P-SV-waves takes 45 per cent less computation time and 29 per cent less memory than the propagator method of Koketsu (1987a, b). The gains are reduced in a multilayer case because of the level of computation required to calculate the addition rules for the large reflection and transmission supermatrices.

Notes: For quasi-stratified media in which the principal variation in seismic properties is with depth, propagation invariants can be constructed from certain combinations of the displacement and tractions elements of two seismic wavefields. These invariants are independent of depth and vanish for identical wavefields, and are constructed for anisotropic, laterally varying media in the spatial and wavenumber domains.These propagation invariants can be exploited to substantially simplify the construction of reflection and transmission processes in laterally varying media, including coupling between wavenumbers. The implementation of this approach is illustrated by application to the incidence of SH-waves on an irregular interface below a free surface. The results are in excellent agreement with those from other schemes but take about 20 per cent less computation time. Even greater improvements in calculation speed are possible in more complex models.

Notes: Coupled mode techniques for guided wave propagation are extended to 2-D stochastic heterogeneity superimposed on a stratified medium. This approach requires the variations to be smoothly varying and of modest size (less than ±2 per cent). By averaging over an ensemble of statistically similar models, coupled equations for the modal energy transport can be generated. The intermode coupling depends on the horizontal correlation functions for the heterogeneity in the crust and mantle, and the integrated effect of the vertical variations in velocity and the modal eigenfunctions.For a particular stochastic model, the attenuation of a single mode as a function of distance can be calculated as a superposition of intrinsic attenuation and scattering loss by energy transfer to other modes of propagation. These statistical estimates of attenuation can be compared with observations of regional phases travelling over a variety of paths in a single region. For Lg and Sn phases, intermode scattering may represent up to 30 per cent of the apparent loss.

Notes: Temporary array deployments of short-period seismometers in northern Australia have been used to build up composite record sections for waves interacting with the upper mantle. Stable measures of the seismic wavefield are provided by stacking the complex envelopes of all the seismic waveforms falling in a 10km distance interval away from the source.Two groups of sources (a) along the Flores Arc, Indonesia with propagation under northwestern Australia, and (b) in New Guinea with paths to the NNE of the array, have been used to construct composite record sections for both P and SV waves over the distance range 1300–2800 km. the timing and amplitude distributions for P waves from the two regions show noticeable differences. Detailed modelling of the record sections yields velocity models with significant variation in velocity for the two sets of propagation paths for which the midpoints are separated by about 1000km.The short-period SV-wave sections indicate efficient propagation of highfrequency S waves in a lithosphere extending down to 210km. Arrivals from the deeper mantle cannot be correlated with confidence because of a loss in high-frequency content revealed by broad-band observations. This requires a significant attenuation zone for S beneath 210 km.

Notes: Earthquakes in oceanic areas are normally located using traveltime tables which are representative of continental paths, since most seismic stations lie on continents. It should therefore be possible to improve such locations by employing a set of traveltimes more appropriate to paths from oceanic events to continental stations.A comparison has therefore been made between locations for a number of oceanic events using the recent iasp91 global traveltimes and the times for the pac91 model derived from observations of events in the Pacific. Although there were often significant differences in the location estimates for the two models, these were often no larger than the shifts induced by changing the misfit criterion used for determining the location.For events in purely oceanic regions such as Tonga and the Marianas with little nearby continent, the results from the pac91 model either provided a significantly better fit to the data or produced depth estimates in close accord with independent constraints (e.g. centroid moment tensor locations). In these cases the use of a specific set of ‘oceanic’ traveltimes can be recommended. However for marginal zones and island arcs, the situation is less clear and it is probably best to employ the global traveltime set with the use of additional phases to improve depth estimates.