Elsevier

Computers & Geosciences

Volume 35, Issue 2, February 2009, Pages 409-418
Computers & Geosciences

Method for rapid high-frequency seismogram calculation

https://doi.org/10.1016/j.cageo.2008.02.030Get rights and content

Abstract

We present a method for rapid, high-frequency seismogram calculation that makes use of an algorithm to automatically generate an exhaustive set of seismic phases with an appreciable amplitude on the seismogram. The method uses a hierarchical order of ray and seismic-phase generation, taking into account some existing constraints for ray paths and some physical constraints. To compute synthetic seismograms, the COMRAD code (from the Italian: “COdice Multifase per il RAy-tracing Dinamico”) uses as core a dynamic ray-tracing code. To validate the code, we have computed in a layered medium synthetic seismograms using both COMRAD and a code that computes the complete wave field by the discrete wave number method. The seismograms are compared according to a time–frequency misfit criteria based on the continuous wavelet transform of the signals. Although the number of phases is considerably reduced by the selection criteria, the results show that the loss in amplitude on the whole seismogram is negligible. Moreover, the time for the computing of the synthetics using the COMRAD code (truncating the ray series at the 10th generation) is 3–4-fold less than that needed for the AXITRA code (up to a frequency of 25 Hz).

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)

  • F. Hron

    Numerical methods of ray generation in multilayered media

  • K. Aki et al.

    Quantitative Seismology

    (2002)
  • P. Bernard et al.

    The Irpinia (Italy) 1980 earthquake: detailed analysis of a complex normal faulting

    Journal of Geophysical Research

    (1989)
  • M. Bouchon

    A simple method to calculate Green's functions for elastic layered media

    Bulletin of Seismological Society of America

    (1981)
  • M. Bouchon et al.

    Discrete wave-number representation of seismic-source wave fields

    Bulletin of Seismological Society of America

    (1977)
  • V. Cerveny

    Seismic Ray Theory

    (2001)
  • V. Cerveny et al.

    The ray series method and dynamic ray-tracing system for three-dimensional inhomogeneous media

    Bulletin of Seismological Society of America

    (1980)
  • T.J. Clarke

    The complete ordered ray expansion—I. Calculation of synthetic seismograms

    Geophysical Journal International

    (1993)
  • T.J. Clarke

    The complete ordered ray expansion—II. Multiphase body wave tomography

    Geophysical Journal International

    (1993)
  • F. Cotton et al.

    Dynamic stress variations due to shear faults in a plane-layered medium

    Geophysical Journal International

    (1997)
  • V. Farra et al.

    Seismic waveform modeling in heterogeneous media by ray perturbation theory

    Journal of Geophysical Research

    (1987)
  • G. Festa et al.

    The Newmark scheme as a velocity–stress time-staggering: an efficient perfectly matched layers implementation for spectral element simulations of elastodynamics

    Geophysical Journal International

    (2005)
  • 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.

    View full text