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: 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.