The work is concentrated primarily on the numerical modelling of seismic wavefields in complex, laterally varying, isotropic and anisotropic, layered and block structures. Asymptotic methods, such as the ray method and its extensions (the paraxial ray method, the method of summation of Gaussian beams and packets, the perturbation ray method, various hybrid methods) and the methods of finite differences are mainly used. All these methods supplement suitably each other.
An inevitable prerequisite for development, testing and application of the wavefield modelling algorithms is a sophisticated computer representation of seismic models of geological structures, suitable for all applicable methods and numerical algorithms and facilitating the description of a very broad class of 3-D, isotropic or anisotropic, laterally inhomogeneous block structures. The research on the computer representation of 3-D seismic models of geological structures and on the physical meaning of the seismic models has thus been paid a considerable attention.
The seismic model describes the spatial distribution of quantities influencing the propagation of seismic waves (e.g., P and S wave velocities, density, quality factors, or even all 21 complex-valued frequency-dependent components of the stiffness tensor). Seismic model of the medium is, as a rule, specified by a finite set of values (model parameters), and by the rules or procedures expressing the dependence of spatial variations of material properties on these values.
Program package MODEL for the specification of general 3-D isotropic and anisotropic heterogeneous block models with curved interfaces has been developed. The package is designed for a joint use with all contemporary methods (e.g., ray methods and elastic finite differences).
The program package ANRAY can be applied to 3-D heterogeneous anisotropic layered models.
In complex models, the behaviour of rays becomes chaotic and the geometrical spreading, number of arrivals and density of caustic surfaces are exponentially increasing with travel time. The exponential increment may be quantitatively described in terms of the Lyapunov exponents. The quantitative estimations of the Lyapunov exponents and mean Lyapunov exponents for rays in 2-D models without interfaces should be generalized to 3-D models, to models with structural interfaces, and extensively numerically tested.
Lyapunov exponents enable to determine quantitative criteria on models, suitable for ray tracing, in terms of Sobolev scalar products. This enables to construct the optimum models of geological structures for ray tracing by inversion of seismic data. Considerable effort devoted to the construction of the optimum models for ray tracing will continue.
Introduction of ray histories and division of the ray-parameter domain into regions of equal ray histories enabled to considerably simplify and generalize 2-D two-point ray tracing algorithm (see the program package SEIS) and to propose highly accurate and reliable 3-D two-point ray tracing algorithm (see the program package CRT). The concept of ray histories facilitated to introduce calculation of travel times and other quantities on dense rectangular grids of points by interpolation within ray cells, which is of crucial importance in prestack seismic migrations and nonlinear determination of seismic hypocenters. The relevant algorithms may simultaneously be used for improving very perspective wavefront tracing methods, and to allow for the calculation of diffracted edge waves. Note that during last years, the ray histories have also become very important in 3-D computer graphics and animation, where, in contrast to seismology, the rays are straight, which considerably simplifies the determination of possible boundaries between ray histories.
A great effort has been concentrated on the algorithms and computer programs for ray tracing in isotropic 3-D heterogeneous block structures composed of heterogeneous blocks separated by curved interfaces and in anisotropic 3-D heterogeneous layered structures. The ray-tracing is supplemented with dynamic ray tracing and amplitude calculation.
The ray tracing packages have been supplemented with sophisticated two-parametric shooting algorithms for two-point ray tracing in general 3-D heterogeneous block or layered structure models, with programs to calculate ray-theory synthetic seismograms generated by sources situated in isotropic or anisotropic media, and with programs to calculate ray-theory synthetic ground motion diagrams.
The two-point ray tracing code used in the CRT program is universal, applicable not only to the point source (common-shot) initial conditions, but also to other initial conditions, e.g., zero-offset rays or diffracted rays. Moreover, the algorithm is also applicable to other initial-value ray tracing programs than CRT. Attention has also been devoted to the calculation of two-point rays diffracted from edges and corner points.
The two-point ray tracing code used in the present version of the ANRAY program is applicable only to the point source (common-shot) initial conditions.
Program package CRT was designed to compute ray diagrams, travel times, vectorial ray amplitudes, polarization vectors, ray synthetic seismograms and ray-theory Green functions for seismic body waves propagating in complex 3-D structures, specified in program package MODEL. Arbitrary type of elementary seismic body wave corresponding to the zero-order ray theory may be treated. Program package may be also applied to compute the qS waves propagating in weakly anisotropic media (using the coupling ray theory), to the Gaussian beam and Gaussian wave packet computations, to the solution of various inversion problems, etc. The aquisition schemes under consideration include surface seismics (land and marine), VSP, cross-hole, OBS, OBC and various seismological ones. Program package CRT has been continuously updated and innovated.
The program package ANRAY to compute rays, travel times, ray amplitudes, ray synthetic seismograms, ray synthetic ground motions and ray-theory Green functions in laterally varying 3-D anisotropic layered structures was designed. It also includes the computation of Green functions of qS waves in weakly anisotropic media using the quasi-isotropic approximation. The package can be applied in surface seismics (land and marine), VSP, cross-hole, OBS and OBC configurations. The package has been permanently updated and innovated.
Attention is devoted to singular directions, classification of singularities, phase shift due to caustics, applicability of Gaussian beams at the singularities, perturbation methods, Finslerian wave-propagation metric tensor, Hamiltonian formalism and Fermat's principle in anisotropic media, and to many other aspects of ray methods in anisotropic media.
The coupling ray theory, which provides continuous transition between the isotropic and anisotropic ray theories, is particularly important at degrees of anisotropy, inhomogeneity and frequencies typical in seismic exploration and structural seismology on all scales (see the program packages CRT and ANRAY). Attention will be devoted to the rederivation of the coupling ray theory directly from the elastodynamic equations and to the accuracy of various modifications of the coupling ray theory. The dynamic ray tracing corresponding to the anisotropic common S-wave ray tracing will be derived. The errors of the coupling ray theory due to the differences between the reference and exact rays will further be studied.
Weakly anisotropic media are typical in seismic exploration and structural seismology. For such media simplified formulae describing various aspects of wave propagation can be used. Simplified formulae provide a clear insight into the relations between wave phenomena and model. They have various applications in forward and inverse wave modeling.
Recently derived equations for the second-order and higher-order perturbations of travel time in heterogeneous isotropic and anisotropic media can be applied to solution of many forward and inverse problems of wave propagation, and to estimation of the accuracy of various computational methods. First applications include the estimation of the error of the coupling ray theory due to isotropic and anisotropic common S-wave ray approximations, but many other applications are anticipated. The equations may be useful in solving many problems of seismic wave propagation. Effort will be expended to generalize the equations for the second-order and higher-order perturbations to the models with structural interfaces.
First-order formulae for the phase velocity and polarization vectors have found many applications in forward and inverse problems. They are important in tomographic studies as well as in the local determination of anisotropy from observed measurements of travel times and polarization. Accuracy of the first-order formulae is studied by using higher-order terms in the perturbation series. Using the perturbation theory, approximate formulae relating group-velocity and phase-velocity vectors were found.
Relations of the reflection or, possibly, transmission coefficients to parameters of media surrounding a reflector is a basic part of the amplitude-versus-offset (AVO) or amplitude-versus-azimuth (AVA) analyses. General formulae for anisotropic media are rather complicated and non-transparent. Using the perturbation theory, simplified and transparent formulae fore the R/T coefficients were derived.
An important role in the recent ray methods and in their various extensions is played by the dynamic ray tracing. The dynamic ray tracing can be performed numerically along a known ray in several ways. Some of them are preferable for the numerical efficiency, some others are more useful in theoretical applications. Transformation relations connecting the results of dynamic ray tracing performed by various methods and in different coordinate systems have been proposed. In particular, the transformation equations between six-dimensional dynamic ray tracing in (phase-space) Cartesian coordinates and 4-dimensional dynamic ray tracing in (phase-space) ray-centred coordinates in anisotropic media were derived. The derivation concentrates on the explicit transformation equations of dynamic ray tracing between Cartesian and ray-centred coordinates. Many of the transformation equations have not been published before even for isotropic medium. Also proposed was an efficient way of reducing the number of equations being solved when numerically evaluating the paraxial-ray propagator matrices, both in Cartesian and ray-centred coordinates.
When selecting the travel-time calculation method, one should consider in the first place the kind of travel times to calculate. In particular, the first-arrival travel times (often but not always required in seismic travel-time tomography) should strictly be distinguished from ray-theory travel times (required, for instance, in the studies of seismic energy propagation and in the ray-theory based migrations).
The first-arrival travel times do not correspond to the ray-theory high-frequency asymptotic approximation, they describe the fastest possible propagation of information and are often connected with a negligible amount of energy. The first-arrival travel times are the property of the exact solution of the elastodynamic equations. In this way, they are related rather to the full-wavefield finite differences than to the ray method. First arrival travel times are single-valued and are defined in the whole model. The calculation of the first-arrival travel times is thus much easier than the calculation of the ray-theory travel times.
The ray-theory travel times ("multivalued" travel times) correspond to a selected elementary wave in the asymptotic decomposition of the wavefield. They describe the propagation of energy in the high-frequency asymptotic approximation. In dispersive media, the ray-theory travel times correspond to the velocities at the prevailing frequency. The multivalued ray-theory travel times may thus be deterministically calculated just in models which are sufficiently globally smooth, i.e., which are not only locally smooth but contain just a reasonably small number of distinct heterogeneities. In more complex models, the statistical properties of the multivalued ray-theory travel times instead of the individual values have to be estimated.
The interpolation within ray cells in 3-D models, composed of heterogeneous blocks separated by curved structural interfaces, may be incorporated into a very efficient algorithm of travel-time calculation for Kirchhoff migrations.
The following 5-step interpolation algorithm has been proposed to evaluate the migrated sections:
The algorithm of the network shortest path calculation of rays and travel times of the first arrivals has been developed. In the proposed algorithm, the travel-time error of the computations is estimated in order to achieve the best accuracy. The rough estimate of the relative travel-time error is evaluated locally at all network nodes prior to network ray tracing, and is minimized by means of a proper choice of the sizes of forward stars. In this way, the structure of the network is adjusted for a particular model and for a particular computer memory. After network ray tracing, the error estimate is refined and the absolute error bounds of the calculated travel times are evaluated. The last version of program NET also includes the accuracy improvement of the two-point rays by means of the network shortest-path ray tracing within the Fresnel volumes.
A new, very accurate computational scheme for calculating the first-arrival travel times on a rectangular grid of points has been proposed. Whereas the network shortest path ray tracing and most of the presently used "finite-difference" eikonal solvers are the methods of the first order, the new proposed method is of the second-order accuracy. It means that the relative propagation-velocity error of calculated travel time is proportional to the second power of the grid step. That is why the method should be sufficiently accurate for all applications in smooth seismic models. Another advantage is the direct calculation of the slowness vectors that increases considerably the accuracy of estimating travel-time derivatives. The equations enabling to estimate the travel-time errors of the method have been derived. The accuracy of the proposed second-order method has been numerically compared with other related techniques.
The intention of this research is to extend the procedures developed for elastic isotropic and anisotropic media to viscoelastic (dissipative) media.
Ray perturbation theory was used to derive equations for two-point complex-valued rays in anisotropic inhomogeneous weakly dissipative media without structural interfaces. The reference medium is assumed to be perfectly elastic (non-dissipative), and the reference ray connecting two fixed points to be real-valued. It is intended to extend the procedure also for the medium containing structural interfaces.
The properties of inhomogeneous plane waves propagating in viscoelastic anisotropic media are studied. General procedures to determine the complex-valued slowness vector of inhomogeneous plane waves are proposed. It was shown that the inhomogeneous waves propagating in viscoelastic anisotropic media exhibit certain phenomena not known from viscoelastic isotropic and from elastic anisotropic media. It is intended to apply these algorithms locally even in general inhomogeneous anisotropic viscoelastic media.
Reflection/transmission coefficients of inhomogeneous plane waves at a plane interface between two homogeneous isotropic dissipative halfspaces have been studied. Material dispersion connected with the dissipation is also considered. Particular attention has been devoted to the effects of inhomogeneity of the incident wave and to the dissipative properties of the model on the reflection/transmission coefficients.
The "non-ray" wave propagation effects, corresponding to the structural details which cannot be included in the models for ray tracing, can often be taken into account by a proper extension of the used ray method. The combination of the matrix methods designed for 1-D stacks of fine layers with 2-D and 3-D ray methods is a nice example.
Complex-valued rays are important in the presence of singularities in elastic media and in dissipative media. The properties of complex-valued rays in both elastic and dissipative media will be studied. Attention will be devoted to the reflection and transmission of complex-valued rays, to the Hamiltonian formulation and possibility to trace complex-valued rays, and to the paraxial and perturbative approximations of complex-valued rays.
The Gaussian-beam and Gaussian-packet summation methods may be obtained by means of the high-frequency ray theory applied to a special case of the Gabor transform of the wavefield. The solution of the respective "modified" eikonal equation need not be unique and may yield more general high-frequency approximations of the wavefield than have been derived yet. Moreover, the high-frequency ray theory applied to a general Gabor transform (or coherent-state transform) of the wavefield can yield even more general forms of the "modified" eikonal and transport equations. These generalizations of the ray theory should be paid considerable attention.
The radiation patterns of seismic point sources situated at the Earth's surface or close to it have been studied. The algorithms of the computation of synthetic seismograms for sources and receivers situated at the Earth's surface, at structural interfaces, at thin low-velocity layers, and close to them have been proposed.
Importance of incorporating the first-order additional terms of the ray method into the ray computations in isotropic media has been studied.
Similarly as in the recent years, considerable attention will be concentrated on the estimation of the accuracy and efficiency of various methods of the calculation of seismic waves and their properties, and of various methods of inversion of seismic data.
Attempts will be made to derive quantitative criteria of applicability and accuracy of ray methods.
When evaluating the Kirchhoff integral, only some regions, often concentrated close to rays, considerably contribute to the integral, whereas the contribution of rapidly oscillating parts of the integrand is usually negligible. Moreover, the numerical quadrature of rapidly oscillating parts is very expensive.
The accuracy and efficiency of the numerical quadrature of the Kirchhoff integrals were studied, and the weighting function localizing the Kirchhoff integral was derived. The weighting function minimizes the number of discrete values necessary for the numerical quadrature. The effective Fresnel zones were then derived as the minimum areas of integration with the specified accuracy. The new, explicit, and very general definition of effective Fresnel zones is purely local, independent of the reference travel times. The definition of the effective Fresnel zones is expressed in terms of the first and second derivatives of the sum of travel times along the surface of integration. The new, quantitative definition of effective Fresnel zones is conformal with other, usually qualitatively or empirically derived definitions, within the regions of their applicability.
Simultaneously, the properties of the classical Fresnel zones and volumes, defined in terms of isochrones of the Kirchhoff integral or representation theorem, were studied.
Ray-based Born approximation enables to study the information on the geological structure, carried by the wavefield, and to study its mapping on various models of the medium. We will continue in the study of the resolution of wave inversion methods, both refraction and reflection, including the resolution and accuracy of seismic migrations. Attention will be paid to the physical meaning of effective material parameters, of seismic models and of the migrated sections, and to the sensitivity of the migrated sections to the velocity model, including its anisotropy.
Whereas the exact material properties are fractal in the nature, the material properties described by the model are, as a rule, smoothed approximations of exact ones. The smoothness is thus a natural property of a seismic model. Moreover, probably all contemporary numerical methods of wavefield or travel-time calculation require in some sense smooth models.
To estimate the accuracy of the model, i.e., the difference between the model and the exact material properties, we need to estimate at least the second statistical moment, usually expressed in terms of the medium correlation function. Possibility to estimate the correlation function under assumption of a self-affine (scalable) random medium was thus studied. The power-law medium correlation function resulted in the travel-time variances proportional to a power of ray lengths. A method designed to estimate the medium correlation function using field travel-times has then been proposed.
The correlation functions describing the differences between the model and the geological structure can be estimated from travel times and other seismic data. Equations for the travel-time tomography enabling not only to invert the travel times but also to estimate the accuracy of the resulting seismic model have been derived. Travel-time inversion based on the minimization of the differences between the model and the geological structure provides automatic determination of travel-time weighting factors necessary for uneven data coverage, and yields the estimation of the correlation functions, which may be used to estimate the accuracy of synthetic travel times and other quantities calculated in the model, the accuracy of the hypocentre positions determined in the model, etc.
We shall continue the development of the theory, algorithms and programs applicable in seismic travel-time tomography and in inversion of coherency panels, with emphasis on the estimation of their accuracy.
A procedure for the determination of parameters of an arbitrary anisotropic medium at a receiver in a borehole from the qP-wave data obtained from a multi-azimuth multiple-source offset VSP experiment is proposed and tested. The data are travel times at several receivers in a borehole and polarization information at the considered receiver. Sensitivity of the results of inversion to the number and orientation of profiles, their length and density of sources along them, to the number and type of the wave is studied. Stability of individual inverted parameters of the medium is investigated.
Gabor signals in the Gabor transform have fixed envelopes, which makes the forward and inverse Gabor transform simple. Gaussian packets, optimized with respect to the elastodynamic equation, have envelopes considerably dependent on the frequency and the direction of propagation, and moderately dependent on the position and time. This makes the decomposition of a general wavefield into Gaussian packets intricate. However, the last research results suggest a possible existence of an approximate expansion into Gaussian packets in the form of a linear integral transform. The inaccuracy of the expansion is probably caused by moderate space-time dependence of the envelopes of Gaussian packets. The inaccuracy can probably be controlled by slight modifications of the envelopes.
Decomposition of a general wavefield into Gaussian packets of given, optimized envelopes is crucial for Gaussian-packet true-amplitude prestack depth migrations. It may also enable to develop general hybrid methods, combining the Gaussian packet summation method with finite differences, finite elements, and other highly accurate methods which can be applied only to small parts of large models at high frequencies. If a successful decomposition of a general wavefield into Gaussian packets were made sufficiently fast, it would allow to develop a new method using Gaussian packets locally, between two consecutive time levels. Such a method could fill the gap between finite-difference and ray methods.
The basic steps of the true-amplitude common-source prestack depth migrations based on Gaussian packets are:
The local paraxial ray approximation allows for deriving the equations for errors of finite-difference schemes due to the grid dispersion. The accuracy of "point" finite-difference schemes of the 2nd and 4th order in 2-D and 3-D regular rectangular grids was studied, and the method of designing the schemes and estimating their accuracy has been proposed. The inaccuracy of finite-difference schemes is governed, above all, by the error in the phase velocity, caused by discretization. This error was estimated for several finite-difference schemes. It is explicitly dependent on the direction of propagation and on wave polarization. The maximum phase-velocity error over all directions of propagation enables the accuracy of the individual schemes to be appreciated in order to select the best one. The proposed approach is general, and applicable to other finite-difference schemes, for example, of the 6th and higher orders.
A local application of the paraxial ray approximation to finite-difference schemes enables optimum finite-difference schemes at structural interfaces and curved earth surface to be derived. Whereas the above common, simple, and highly symmetric point schemes are applicable only in smooth parts of seismic models, inside geological blocks, the method of designing the finite-difference schemes of the 2nd and 4th order at structural interfaces is required. A local application of the paraxial ray approximation to finite-difference schemes enables optimum finite-difference schemes at structural interfaces and curved earth surface to be derived. These special schemes are linear operators acting on the wavefield gridpoint values, and are dependent on the positions and slopes of the interfaces.
SW3D
- main page of consortium Seismic Waves in Complex 3-D Structures .