A method for the computation of finite frequency body wave synthetic seismograms in laterally varying media

Summary. Body wave synthetic siesmograms for laterally varying media are computed by means of a slowness implementation of the extended WKBJ (EWKBJ) theory of Frazer & Phinney. An EWKBJ seismogram is computed by first tracing rays through a particular model to obtain conventional ray information (travel time, ray end point, ray slowness) and then using these data in the finite frequency integral expression for the EWKBJ seismogram. The EWKBJ seismograms compare favourably to geometrical ray theory (GRT) seismograms but are significantly better because of the finite frequency nature of the EWKBJ calculation. More realistic behaviour is obtained with EWKBJ seismograms at normal seismic frequencies near caustics, where the GRT amplitude is infinite, and within geometrical shadow zones where GRT predicts zero amplitudes. In addition the EWKBJ calculation is more sensitive than GRT to focuses and defocuses in the ray field. The major disadvantage of the EWKBJ calculation is the additional computer time over that of GRT, necessary to calculate one seismogram although an EWKBJ seismogram costs much less to compute than a reflectivity seismogram. Another disadvantage of EWKBJ theory is the generation of spurious, non-geometrical phases that are associated with rapidly varying lateral inhomogeneities. Fortunately the amplitudes of these spurious phases are usually much lower than that of neighbouring geometrical phases so that the spurious phases can usually be ignored. When this observation is combined with the moderately increased computational time of the EWKBJ calculation then the gain in finite frequency character significantly outweighs any disadvantages.     

A method for calculating synthetic seismograms in laterally varying media

Summary An effective algorithm for computing synthetic seismograms in laterally inhomogeneous media has been developed. The method, based on zero-order asymptotic ray theory, is primarily intended for use in refraction and reflection studies and provides an economical means of seismic modelling.
A given smoothed velocity-depth-distance model is divided into small squares with constant seismic parameters and first-order interfaces are represented by an arbitrary number of dipping linear segments. The computation of ray propagation and amplitudes through such a model does not involve complicated analytic expressions and therefore minimizes computer time.
Amplitudes are determined by geometrical spreading of spherical wave-fronts and energy partitioning at interfaces. Synthetic seismograms calculated for laterally homogeneous models are in good agreement with those obtained by the Reflectivity Method.     

Synthetic body wave seismograms for three-dimensional laterally varying media

Summary. Two methods of computing body wave synthetic seismograms in three-dimensional laterally varying media are discussed. Both these methods are based on the summation of Gaussian beams. In the first, the initial beam parameters are chosen at the source, in the second at the beam endpoints. Both these variants eliminate the ray method singularities. The expansion of the wavefield into plane waves may be considered as the limiting case of the first approach and the Chapman–Maslov method as the limiting case of the second approach. Computer algorithms are briefly described and numerical examples presented. In the first numerical example, the comparisons of the two approaches, based on summing Gaussian beams, with the reflectivity method indicate that the computed synthetic seismograms are satisfactorily accurate even in the caustic region. The next example suggests that the two methods discussed can be simply and effectively applied to 3-D laterally inhomogeneous structures.     

The theory of finite frequency body wave synthetic seismograms in inhomogeneous elastic media

b
A new method is presented by means of which one can compute finite frequency synthetic seismograms for media whose velocity and density are continuous functions of two or three spatial variables. Basically, the method is a generalization of the familiar phase integral method, to which it reduces in a stratified medium. For a given source location the travel-time and distance functions needed to compute synthetics are obtained by numerically tracing rays through the model. This information is then used to evaluate a double integral over frequency and take-off angle at the source. The solution obtained reduces to the geometrical optics solution wherever that is valid but it also works in shadows and at caustics without knowing explicitly where these may be located. The method can be used as a spectral method, in which the integral over take-off angle is evaluated first, or as a slowness method, in which the frequency integral is evaluated first.     

Synthetic body wave seismograms for laterally varying layered structures by the Gaussian beam method

Summary. Several approaches to computing body wave seismograms in 2–D and 3–D laterally inhomogeneous layered structures are suggested. They are based on the Gaussian beam method, which has been recently applied to the evaluation of time-harmonic high-frequency wavefields in inhomogeneous media. Three variants are discussed in some detail: the spectral method, the convolutory method and the wave-packet method. The most promising seems to be the wave-packet approach. In this approach, the wavefield, generated by a source, is expanded into a system of wave packets, which propagate along rays from the source in all directions. The wave packets change their properties due to diffusion, spreading, reflections/transmissions, etc. The resulting seismogram at any point of the medium is then obtained as a superposition of those packets which propagate close to the point. The final expressions in all the three methods are regular even in regions, in which the ray method fails, e.g. in the vicinity of caustics, in the critical region, at boundaries between shadow and illuminated regions, etc. Moreover, they are not as sensitive to the minor details of the medium as the ray method and, what is more, they remove the time-consuming two-point ray tracing from computations. Numerical examples of synthetic seismograms computed by the wave-packet approach are presented.     

First-order P-wave ray synthetic seismograms in inhomogeneous weakly anisotropic media

We propose approximate equations for P -wave ray theory Green's function for smooth inhomogeneous weakly anisotropic media. Equations are based on perturbation theory, in which deviations of anisotropy from isotropy are considered to be the first-order quantities. For evaluation of the approximate Green's function, earlier derived first-order ray tracing equations and in this paper derived first-order dynamic ray tracing equations are used.
The first-order ray theory P -wave Green's function for inhomogeneous, weakly anisotropic media of arbitrary symmetry depends, at most, on 15 weak-anisotropy parameters. For anisotropic media of higher-symmetry than monoclinic, all equations involved differ only slightly from the corresponding equations for isotropic media. For vanishing anisotropy, the equations reduce to equations for computation of standard ray theory Green's function for isotropic media. These properties make the proposed approximate Green's function an easy and natural substitute of traditional Green's function for isotropic media.
Numerical tests for configuration and models used in seismic prospecting indicate negligible dependence of accuracy of the approximate Green's function on inhomogeneity of the medium. Accuracy depends more strongly on strength of anisotropy in general and on angular variation of phase velocity due to anisotropy in particular. For example, for anisotropy of about 8 per cent, considered in the examples presented, the relative errors of the geometrical spreading are usually under 1 per cent; for anisotropy of about 20 per cent, however, they may locally reach as much as 20 per cent.     

Seismic body waves in anisotropic media: synthetic seismograms

Summary. Synthetic seismograms and particle motion diagrams are computed for simple, layered Earth models containing an anisotropic layer. The presence of anisotropy couples the P, SV and SH wave motion so that P waves incident on the anisotropic layer from below produce P, SV and small-amplitude SH waves at the surface both the P velocity and the amplitudes of the converted phases vary with azimuth. Significant SH amplitudes may be generated even when the wavelength of the P wave is much greater than the thickness of the anisotropic layer. Incident SV or SH waves may each generate large amplitudes of both SV and SH motion. This strong coupling is largely independent of the degree of velocity anisotropy of the medium. The arrivals from short-period S waves exhibit S-wave splitting, but arrivals from longer period S waves superpose into a modified waveform. This strong coupling does not allow the arrival of separate phases with pure SV and SH polarization except along directions of symmetry where the motion decouples.     

Ray amplitudes of seismic body waves in laterally inhomogeneous media

Summary . The most complicated part in the computation of ray amplitudes of seismic body waves in laterally inhomogeneous media with curved interfaces lies in the evaluation of the geometrical spreading. Geometrical spreading can be simply expressed in terms of the Jacobian J of the transformation from the Cartesian into ray coordinates. Several systems of ordinary differential equations to compute the function J are suggested. For general three-dimensional media, in which the velocity changes with all the three spatial coordinates, a system of three non-linear ordinary differential equations of the first order is derived. If the velocity does not depend on one coordinate, the system of equations reduces to only one non-linear differential equation. The initial conditions for these differential equations at point (or line) source and at points of intersection of the ray with curved interfaces are presented.     

A WKBJ spectral method for computation of SV synthetic seismograms in a cylindrically symmetric medium

Summary. A method of synthetic seismogram computation for teleseismic SV -waves is developed in order to treat quantitatively SV -waves in problems of body wave source inversion and source—receiver structure studies. The method employs WKBJ theory for a generalized ray in a vertically inhomogeneous half-space and the propagator matrix technique for waves in near-surface homogeneous layers. Wavenumber integration is done along the real axis of the wavenumber plane and anelasticity is included by using complex velocity in all regions of the earth model. The near-surface source structure is taken into account in the computation for the case of the shallow source by allowing a point source to be located in the homogeneous layers. Source and receiver area structures are also allowed to differ. A general moment tensor point source is considered.     

A new method for computing synthetic seismograms

Summary. The computation of theoretical seismograms for models in which the elastic parameters and density vary only with depth (in a plane, cylindrical or spherical geometry) reduces to the solution of an ordinary differential equation plus the evaluation of inverse transformations. In principle, the problem is straightforward. In practice, many techniques and approximations can be used at each stage and many combinations and variants are possible. In this paper, we discuss a new method of evaluating the inverse transforms. Any method can be used to solve the differential equation and we only discuss a few analytic approximations to illustrate the new method. The inverse transformations are a frequency and wavenumber integral. Essentially four techniques can be used to evaluate these depending on the order of integration and whether the wavenumber integral is distorted from the real axis. Three of these have been widely used, but the technique of evaluating the frequency integral first and keeping the wavenumber real is new. In this paper, we discuss some of the advantages of this combination.     

Efficient calculation of Gaussian-beam seismograms for two-dimensional inhomogeneous media

Summary. This paper discusses several aspects of the calculation of theoretical seismograms for two-dimensional inhomogeneous media with the method of Gaussian beams. The most important steps of this method, kinematic and dynamic ray tracing, can be performed very efficiently, if the model cross-section is subdivided into triangles with linear velocity laws. Each Gaussian beam is characterized by a complex beam constant ε which determines its width and phase-front curvature. Various possibilities to choose ε are discussed, including cases where beam properties at the beam endpoint (and not at the beginning) are prescribed; for instance, the beam width at the endpoint can be specified. In such cases the beam constant is a function of the radiation angle at the source, and the decomposition of a cylindrical wave into beams has to take this into account by weighting the beams differently, at least in principle. The exact weight function is derived and shown to be reasonably well approximated by the weight function, corresponding to angle-independent ε Theoretical seismograms are presented for a laterally heterogeneous model of the crust–mantle transition which is characterized by complications in the reflection from the transition and in the refraction from below. These complications are modelled by and large with success. The seismograms, however, depend to a certain extent on the choice of the beam constant. Moreover, according to the reciprocity principle calculations with source and receiver interchanged should have the same results as calculations for the original configuration. In practice this is not so, and the difference increases with the strength of lateral heterogeneities. Hence, for a successful application of Gaussian beams the model should not vary too strongly in lateral direction.     

Computation of complete synthetic seismograms for laterally heterogeneous models using the Direct Solution Method

We use the Direct Solution Method (DSM) together with the modified operators derived by Geller & Takeuchi (1995) and Takeuchi, Geller & Cummins (1996) to compute complete synthetic seismograms and their partial derivatives for laterally heterogeneous models in spherical coordinates. The methods presented in this paper are well suited to conducting waveform inversion for 3-D Earth structure. No assumptions of weak perturbation are necessary, although such approximations greatly improve computational efficiency when their use is appropriate.
An example calculation is presented in which the toroidal wavefield is calculated for an axisymmetric model for which velocity is dependent on depth and latitude but not longitude. The wavefield calculated using the DSM agrees well with wavefronts calculated by tracing rays. To demonstrate that our algorithm is not limited to weak, aspherical perturbations to a spherically symmetric structure, we consider a model for which the latitude-dependent part of the velocity structure is very strong.     

Solution of the inverse problem of seismology for laterally inhomogeneous media

Summary. The Backus-Gilbert method has been extended to the estimation of the seismic wave velocity distribution in 2-D or 3-D inhomogeneous media from a finite set of travel-time data. The method may be applied to the inversion of body wave as well as surface wave data. The problem of determining a local average of the unknown velocity corrections may be reduced to a choice of a suitable δ-ness criterion for the averaging kernel. For 2-D and 3-D inhomogeneous media the simplest criterion is to minimize a sum of 'spreads' over all the coordinates. The use of this criterion requires the solution (the averaged velocity corrections) to be represented as a sum of functions, each of which depends only on one coordinate. This is a basic restriction of the method. In practice it is possible to achieve good agreement between the solution and a real velocity distribution by a reasonable choice of the coordinate system.
Numerical tests demonstrate the efficiency of the method. Some examples of the application of the method to the inversion of real seismological data for body and surface waves are given.     

Efficient global matrix approach to the computation of synthetic seismograms

Summary. A numerically efficient global matrix approach to the solution of the wave equation in horizontally stratified environments is presented. The field in each layer is expressed as a superposition of the field produced by the sources within the layer and an unknown field satisfying the homogeneous wave equations, both expressed as integral representations in the horizontal wavenumber. The boundary conditions to be satisfied at each interface then yield a linear system of equations in the unknown wavefield amplitudes, to be satisfied at each horizontal wavenumber. As an alternative to the traditional propagator matrix approaches, the solution technique presented here yields both improved efficiency and versatility. Its global nature makes it well suited to problems involving many receivers in range as well as depth and to calculations of both stresses and particle velocities. The global solution technique is developed in close analogy to the finite element method, thereby reducing the number of arithmetic operations to a minimum and making the resulting computer code very efficient in terms of computation time. These features are illustrated by a number of numerical examples from both crustal and exploration seismology.