Method for rapid high-frequency seismogram calculation☆
Introduction
The calculation of synthetic seismograms has recently become a useful tool in seismological research, and a wide variety of techniques have been developed. Forward modelling, which means in this case the generation of synthetic seismograms, represents an important part of many seismological studies, such as seismic tomography or the kinematic inversion of source parameters (Pengcheng and Archuleta, 2004). The solution of an inverse problem requires the repeated solving of the forward problem, so that speed is the most stringent condition that the method must supply.
For horizontally stratified Earth models the complete wave field computation can be performed by using the theory of normal mode (Aki and Richards, 2002; Rosenbaum, 1960) or the discrete wave number method (Bouchon and Aki, 1977). Methods based on the direct numerical solution of the elastodynamic equation by finite difference (Madariaga, 1976; Virieux, 1986), finite elements (Marfurt, 1984), spectral elements (Festa and Vilotte, 2005; Komatitsch and Vilotte, 1998), or discontinuous Galerkin method (Kaser and Dumbser, 2006) can be applied to study the seismic wave field in a complex laterally varying structure.
As an alternative it is possible to obtain the solution to the wave equation by using approximate high-frequency methods (Cerveny, 2001). The solution of the elastodynamic equation is composed of elementary body waves that correspond to the various rays connecting the source to the receiver. Although the computation of ray synthetic seismograms is only approximate and the ray method can fail in certain situations, the high-frequency methods are preferable to direct numerical methods for many applications, both for the shorter computing times and for the full interpretability of the seismograms. The problems arise when we want to use a synthetic seismogram as similar to the real one as possible in the high-frequency approximation, such that the seismogram should be relatively complete, although it is not necessary that it contains every feature of the full elastic wave field. It is becoming necessary to select from all of the rays connecting the source to the receiver only those that produce an appreciable amplitude on the seismogram. The problem of generating a comprehensive set of rays was approached for the first time by Hron, 1971, Hron, 1972 and Hron et al. (1986) for layered media, by grouping individual rays into families of kinematic equivalents. Afterwards, Clarke, 1993a, Clarke, 1993b developed a technique for computing synthetic seismograms based on a ray-generation algorithm that involved the symbolic manipulation of complete wave field expressions from the reflectivity theory, which were truncated to produce a finite ray series.
In this study, we propose a new technique for the rapid definition of an exhaustive set of rays that is based on the hierarchic generation of strings that describe the ray paths and the phase types. The string generation is subjected to physical constraints that are related to the propagation medium and the source–receiver geometry. The ray sets will represent the input of a kinematic or dynamic ray-tracing algorithm (i.e. Cerveny and Hron, 1980; Farra and Madariaga, 1987; Snieder and Spencer, 1993; Virieux, 1991). In particular, the technique developed has been implemented in the multiphase dynamic ray-tracing COMRAD code (from the Italian: “COdice Multifase per il RAy-tracing Dinamico”) that uses as its core the dynamic ray-tracing code provided by Farra and Madariaga (1987).
Section snippets
Method
Our goal is to develop an algorithm that rapidly generates an exhaustive number of seismic phases to calculate a high-frequency seismogram as complete as possible. In particular, the computing time is a critical parameter when the computation of synthetic seismograms has to be repeated many times in the inversion processes. The amplitudes, the ray paths and the travel times of the seismic phases are computed by the dynamic ray-tracing code provided by Farra and Madariaga (1987). In this section
Structure of the COMRAD code
The Comrad.f code is the Fortran77 version of multiphase code that, based on selection criteria, generates an exhaustive set of ray and phase strings allowing the computation of a seismogram as complete as possible. In order to compute synthetic seismograms, that is travel times, amplitudes, and ray paths relative to the selected phases, it uses as its core the dynamic ray-tracing code developed by Farra and Madariaga (1987), which needs a medium discretized by flat tilted layers. However, the
Validation of the method
The higher the number of generations we consider, the higher the number of phases the core must calculate. To make the best choice on the number of generations, we have to compare the calculation time with the complexity we need for the synthetics. We compute synthetic seismograms for two receivers (R1 and R2), respectively, at 1 and 30 km distance from the epicentre of an explosive source. The source function is a triangle with a duration of 0.1 s, but it is possible to choose other source
Conclusions and discussions
We have developed a method for rapid high-frequency seismogram computation based on a hierarchical order of ray and phase strings generation subjected to physical constraints. The method is used in a code (Comrad.f) to compute an exhaustive set of phases able to produce a complete body waves seismogram. It uses as core the dynamic ray-tracing code developed by Farra and Madariaga (1987).
We have numerically tested the method in two steps using a PC with an AMD-3 GHz processor and 2 Mb SD-RAM
Acknowledgements
The authors would like to thank Veronique Farra for providing us with her dynamic ray-tracing code, used as core of the COMRAD program. We are also grateful to Prof. Ray Durrheim for helpful suggestions and constructive comments on the manuscript.
References (25)
Numerical methods of ray generation in multilayered media
- et al.
Quantitative Seismology
(2002) - et al.
The Irpinia (Italy) 1980 earthquake: detailed analysis of a complex normal faulting
Journal of Geophysical Research
(1989) A simple method to calculate Green's functions for elastic layered media
Bulletin of Seismological Society of America
(1981)- et al.
Discrete wave-number representation of seismic-source wave fields
Bulletin of Seismological Society of America
(1977) Seismic Ray Theory
(2001)- et al.
The ray series method and dynamic ray-tracing system for three-dimensional inhomogeneous media
Bulletin of Seismological Society of America
(1980) The complete ordered ray expansion—I. Calculation of synthetic seismograms
Geophysical Journal International
(1993)The complete ordered ray expansion—II. Multiphase body wave tomography
Geophysical Journal International
(1993)- et al.
Dynamic stress variations due to shear faults in a plane-layered medium
Geophysical Journal International
(1997)
Seismic waveform modeling in heterogeneous media by ray perturbation theory
Journal of Geophysical Research
The Newmark scheme as a velocity–stress time-staggering: an efficient perfectly matched layers implementation for spectral element simulations of elastodynamics
Geophysical Journal International
Cited by (3)
Seismic networks layout optimization for a high-resolution monitoring of induced micro-seismicity
2020, Journal of SeismologyElectric and magnetic field changes observed during a seismic swarm in Pollino Area (Southern Italy)
2014, Bulletin of the Seismological Society of AmericaA comprehensive approach for evaluating network performance in surface and borehole seismic monitoring
2013, Geophysical Journal International
- ☆
Code available from server at http://www.iamg.org/CGEditor/index.htm.
- 1
Tel.: +39 824 323658; fax: +39 81 2420334.
- 2
Tel.: +39 81 2420315; fax: +39 81 2420334.