• Nie Znaleziono Wyników

Propagation of transient elastic waves in stratified anisotropic media

N/A
N/A
Protected

Academic year: 2021

Share "Propagation of transient elastic waves in stratified anisotropic media"

Copied!
328
0
0

Pełen tekst

(1)

P R O P A G A T I O N OF

T R A N S I E N T ELASTIC WAVES

IN S T R A T I F I E D

A N I S O T R O P I C M E D I A

Joseph H . M . T . van der Hijden

TR diss

1535

(2)

*-"\

Mil

'"t*

PROPAGATION OF

TRANSIENT ELASTIC WAVES

IN STRATIFIED

(3)

PROPAGATION OF TRANSIENT ELASTIC WAVES

IN STRATIFIED ANISOTROPIC MEDIA

-.--o'.-..'./. .'■'-yjif'-lzn,-.

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus, Prof. Dr. J. M. Dirken, in het openbaar te verdedigen ten overstaan van een commissie aangewezen door

. het College van Dekanen op 16 april 1987 te 14.00 uur

door

Joseph Henricus Maria Titus van der Hij den,

geboren te Maastricht, Elektrotechnisch Ingenieur.

(4)
(5)

Aan Inez

(6)

The author is employed by Schlumberger-Doll Research, Ridgefield, Con­ necticut, U.S.A., and gratefully acknowledges the support for research re­ ported here and the permission to write and publish this monograph.

(7)

Contents

I N T R O D U C T I O N 1 1.1 Statement of the problem 1

1.2 The method of solution employed 3

1.3 Numerical considerations 6 B A S I C R E L A T I O N S F O R E L A S T I C WAVES IN S T R A T ­

I F I E D , P I E C E W I S E H O M O G E N E O U S , A N I S O T R O P I C

M E D I A 9 2.1 Introduction 10

2.2 Description of the configuration and formulation of the prob­

lem 10 2.3 Basic equations for the elastic wave motion 14

2.4 The transform-domain equations 17 2.5 The motion-stress vector in a homogeneous subdomain . . . 19

2.6 The transform-domain wave vector in a source-free domain 22 2.7 The transform-domain wave vector and the motion-stress

vector generated by a localized source in a homogeneous

sub-domain 26 2.8 The transform-domain wave vector and the motion-stress

vector generated by a localized source in a stratified medium 29 2.9 The transform-domain generalized-ray wave constituents . . 34 2.10 Transformation of the solution back to the space-time do­

(8)

B A S I C RELATIONS F O R ELASTIC WAVES I N S T R A T ­ IFIED, P I E C E W I S E H O M O G E N E O U S , I S O T R O P I C

M E D I A 4 1 3.1 Introduction . 42

3.2 Description of the configuration and formulation of the prob­

lem 42 3.3 Basic equations for the elastic wave motion 44

3.4 The transform-domain equations 47 3.5 The motion-stress vector in a homogeneous subdomain . . . 47

3.6 The transform-domain wave vector in a source-free domain 50 3.7 The transform-domain wave vector and the motion-stress

vector generated by a localized source in a homogeneous

sub-domain 58 3.7.1 Explosive source 59

3.7.2 Force source 60 3.8 The transform-domain wave vector and the motion-stress

vector generated by a localized source in a stratified medium 61 3.9 The transform-domain generalized-ray wave constituents . . 63 3.10 Transformation of the solution back to the space-time do­

main 64 3.11 Basic relations for acoustic waves in a fluid 64

R A D I A T I O N F R O M A N I M P U L S I V E S O U R C E I N A N U N B O U N D E D H O M O G E N E O U S I S O T R O P I C SOLID 67

4.1 Introduction 68 4.2 Transformation of the solution back to the space-time do­

main 68 4.3 The behavior of s3' in the complex s plane 70

4.4 Cagniard-de Hoop contours in the complex s plane . . . 71

4.5 Space-time expression for the motion-stress vector 74

4.5.1 Full solution in special directions 76 4.6 Alternative implementation of the Cagniard-de Hoop

method 77 4.6.1 The behavior of sf'5 in the complex s plane 80

4.6.2 Cagniard-de Hoop contours in the complex 5 plane . 81 4.6.3 Space-time expression for the motion-stress vector . 84

(9)

4.6.4 Full solution in special direction 91 4.7 Approximations and derived results 93

4.7.1 The Cagniard-de Hoop method for large horizontal

offset 93 4.7.2 The Cagniard-de Hoop method for large vertical off­

set 97 4.7.3 The two-dimensional problem 99

4.7.4 Far-field approximation independent of the propaga­

tion direction 101 4.8 Numerical results 102 4.9 Conclusion I l l 5 R A D I A T I O N F R O M A N I M P U L S I V E S O U R C E I N A S T R A T I F I E D I S O T R O P I C M E D I U M 115 5.1 Introduction 116 5.2 Transformation of the solution back to the space-time do­

main 116 5.3 The two-dimensional problem 118

5.3.1 The behavior of s3!m in the complex s plane 118 5.3.2 Cagniard-de Hoop contours in the complex 5 plane

(two-dimensional case) 119 5.3.3 Space-time expression for the motion-stress vector

(two-dimensional case) 127 5.4 The three-dimensional problem 128

5.4.1 The behavior of s£'^ in the complex 5 plane 129 5.4.2 Cagniard-de Hoop contours in the complex 5 plane

(three-dimensional case) 130 5.4.3 Space-time expression for the motion-stress vector

(three-dimensional case) 138 5.5 Alternative implementation of the Cagniard-de Hoop

method 143 5.5.1 The behavior of s3.'m in the complex s plane 146

5.5.2 Cagniard-de Hoop contours in the complex 5 plane . 146 5.5.3 Space-time expression for the motion-stress vector . 154

(10)

5.6.1 The Cagniard-de Hoop method for large horizontal

offset 158 5.6.2 The Cagniard-de Hoop method for large vertical off­

set 163 5.6.3 Far-field approximation independent of the propaga­

tion direction 165 5.6.4 Relation between the two-dimensional result and the

large horizontal offset approximation 166

5.7 Numerical results 166 5.8 Conclusion 172 6 R A D I A T I O N F R O M A N I M P U L S I V E S O U R C E I N A N U N B O U N D E D H O M O G E N E O U S A N I S O T R O P I C SOLID 179 6.1 Introduction 180

6.2 Transformation of the solution back to the space-time do­

main 181 6.3 The behavior of sfn in the complex s plane 185

6.3.1 Illustration of the behavior of s*n in the complex s

plane using the slowness surface 191 6.4 Cagniard-de Hoop contours in the complex s plane 196

6.5 Space-time expression for the motion-stress vector 208

6.5.1 Full solution in special direction 215 6.6 Approximations and derived results 216

6.6.1 The Cagniard-de Hoop method for large horizontal

offset 216 6.6.2 The Cagniard-de Hoop method for large vertical off­

set 221 6.6.3 The two-dimensional problem 224

6.6.4 Relation between the two-dimensional result and the

large horizontal offset approximation 226 6.7 Numerical aspects of the computation of Cagniard-de Hoop

contours in anisotropic media 227

6.8 Numerical results 231

(11)

7 R A D I A T I O N F R O M A N I M P U L S I V E S O U R C E I N A

S T R A T I F I E D A N I S O T R O P I C M E D I U M 255

7.1 Introduction 256 7.2 Transformation of the solution back to the space-time do­

main 256 7.3 The behavior of sf.^ in the complex 5 plane 259

7.4 Cagniard-de Hoop contours in the complex 5 plane 260 7.5 Space-time expression for the motion-stress vector 266 7.6 Approximations and derived results . 269

7.6.1 The Cagniard-de Hoop method for large horizontal

offset 269 7.6.2 The Cagniard-de Hoop method for large vertical off­

set 271 7.6.3 The two-dimensional problem 273

7.6.4 Relation between the two-dimensional result and the

large horizontal offset approximation 275

7.7 Numerical results 275

7.8 Conclusion 284 A P P E N D I C E S 285 A The slowness surface of a homogeneous anisotropic medium 285

B Up- and down-going waves in a homogeneous anisotropic

medium 291

C The values of afn for large \a\ 294

D The transversely isotropic medium and its stiffness tensor . 296 E Simplified formulation for unbounded medium in the

two-dimensional problem 298

(12)

Chapter 1

INTRODUCTION

S u m m a r y

Seismic waves are one of the standard diagnostic tools that are used to determine the distribution of the mechanical parameters (volume density of mass, compressibility, elastic stiffness) in the interior of the earth and of the geometry of subsurface structures. There is increasing evidence that in the interpretation of the seismic data the influence of anisotropy has to be taken into account, especially when shear-wave data interpretation is aimed at. With this application in mind, we present a method to compute the seismic waves that are generated by an impulsive source in a stratified anisotropic medium.

1.1 S t a t e m e n t of t h e problem

Although the present monograph has been written with the seismic ap­ plications in mind, the method that is developed is not limited to solid-earth geophysics. The classical example of waves in anisotropic media are the elastic waves in crystals (Musgrave 1970). Further, many ultrasonic devices are constructed by using layers of piezoelectric (and, hence, an­ isotropic) materials; often these devices, too, are appropriately described by the idealized geometries discussed in this book. Electromagnetic waves in the ionosphere (Budden 1961, Felsen and Marcuvitz 1973, p.740) and in integrated optical devices (crystal optics) (Born and Wolf 1980, p.665,

(13)

Yariv and Yeh 1984) are further examples where anisotropy of the medium has to be included in the description. In short, the methods discussed in this monograph are applicable wherever waves do propagate in stratified, anisotropic media.

Seismic anisotropy is very common in the earth. There are several causes of anisotropy: the interleaving of thin sedimentary beds, the pres­ ence of preferentially oriented cracks, the occurrence of stress-induced ef­ fects, and the alignment of crystals or grains. Until recently, seismic data were usually interpreted on the basis of compressional-wave analysis. In the latter, anisotropy is expected to be of less importance; so anisotropy could be ignored. With the increasing resolution of seismic observations, the influence of anisotropy is often noticeable, and of particular importance if shear waves are analyzed. In fact, one can safely state that in the real earth there exists no isotropic rock. Therefore, the question to ask is: how strong is the anisotropy and how strong does it manifest itself in elastic wave propagation? In order to answer that quantitative question we need a good understanding of how the observed seismic wave fields change due to anisotropy. Such an understanding can be gained by building appropriate computer models and calculating in some, idealized, model geometries the acoustic wave fields, and their accompanying seismic records. The present monograph describes one of those models, viz. the stratified geometry. In this model, the influence of anisotropy in each layer separately can be stud­ ied.

There are some fundamental differences between acoustic wave prop­ agation in isotropic and in anisotropic elastic media. These differences are already manifest when considering the propagation of uniform plane waves in a homogeneous region with a horizontal boundary. In an isotropic medium one can then distinguish between the compressional (P) and the, vertically or horizontally polarized, shear waves (SV and SH, respectively). The decomposition into these three eigenwaves is carried out on the basis of the polarization of their particle displacement, or their particle velocity, with respect to the horizontal plane, so that "V" means "vertical," and "H" means "horizontal." The particle velocity of the P waves is curl-free, while the particle velocities of the SV and the SH waves are divergence-free.

In a weakly anisotropic medium (weak in the sense that the medium is characterized by constitutive coefficients that are only small

(14)

perturba-tions of the ones pertaining to some isotropic medium) the plane waves in a certain direction of propagation can still be labeled as (with q for quasi) a qP wave, with approximately longitudinal particle velocity polar­ ization, and qSV and qSH waves, with approximately transverse particle velocity polarization (Keith and Crampin 1977b, Aki and Richards 1980, p.188). The decomposition into these three modes is again carried out with respect to the horizontal plane so that "V" means "mostly vertical," and

"H" means "mostly horizontal." In a strongly anisotropic medium there are three plane waves in every direction of propagation with mutually orthog­ onal particle velocity polarizations and with wave speeds that are different and vary with direction. Identification of the waves according to specific dominant particle velocity polarizations is now meaningless. Except in cer­ tain specific symmetry directions, none of the three waves is either curl-free or divergence-free in its particle velocity.

In this monograph we investigate the features of the wave field radiated by a concentrated source in the case where the medium is horizontally strat­ ified and each layer is a homogeneous, and arbitrarily anisotropic, solid. As such, the configuration serves as a canonical problem. Two types of sources which are of seismic interest are considered in detail, viz. a point source of expansion (model for an explosive source) and a point force (model for a mechanical vibrator). The theory is also applicable to point sources of deformation rate, which are considered to be adequate models for earth­ quake generation. The results can serve as a guide to the interpretation of experimental data acquired in the more complicated situations met in practice, where the layers may not be all parallel.

1.2 T h e m e t h o d of solution employed

The standard approach to the problem stated in Section 1.1 is to employ a Fourier transformation with respect to time, and Fourier transformations with respect to the horizontal spatial coordinates. To obtain numerical results, the relevant inverse transforms have to be evaluated numerically, with the possible use of asymptotic evaluations in certain regions of space-time. This set-up of the problem is similar to the frequency-wavenumber integration method used by Booth and Crampin (1983), and by Fryer and Frazer (1984).

(15)

We solve the problem by a different method, viz. by applying the Cagniard-de Hoop method (de Hoop 1960, 1961, see also Miklowitz 1978, p.302, and Aki and Richards 1980, p.224), which, in general, leads to con­ siderably less computational effort. Since by this method the computational results can relatively easily be obtained with any degree of accuracy, they can also be used as a check on the accuracy of the numerical procedures that are employed to evaluate the inversion integrals in the standard treat­ ment of the problem, i.e. the frequency-wavenumber integration method. The latter numerical technique seems to be the only available procedure that can be used in those cases where the materials have an arbitrary loss mechanism, and/or in the case where the geometries involve curved sur­ faces.

We shall address the advantages and disadvantages of our method as compared to various alternative methods at appropriate places in this monograph. The methods that can deal with anisotropy in layered media, besides the Cagniard-de Hoop method, are: the frequency-wavenumber in­ tegration method mentioned above, ray-tracing methods for asymptotic high-frequency results and first-arrival analysis, and finite-element and finite-difference methods, where the latter can deal with arbitrarily inhomo-geneous media, but require enormous computational effort. The Cagniard-de Hoop method and the frequency-wavenumber integration method are both integral-transformation methods and are restricted to a special geom­ etry, such as horizontal stratification and time invariance of the configura­ tion.

The monograph consists of two parts: acoustic waves in isotropic media and acoustic waves in anisotropic media. The separate formulation for isotropic media has been included for didactical purposes. Although the general method is the same, the resulting expressions are simpler than the ones in anisotropic media. By offering the opportunity to look back at the analogous expression for isotropic media, the physical interpretation of the results pertaining to anisotropic media becomes much easier.

However, there is still another reason why we write down explicitly, in Chapters 3 to 5, the analysis of acoustic wave propagation in strati­ fied isotropic media with the aid of the Cagniard-de Hoop method. The reason is that nowhere in the literature the solution of this problem has completely been written up in its simplest form. Many publications on

(16)

this problem can be found, but we feel there is always something miss­ ing. First, Cagniard's (1939) book and the translation by Flinn and Dix

(1962) contain the original, intricate, version of the transformation back to the time domain with the intermediate complex time variable. The latter difficulty was circumvented by de Hoop (1960), who managed to simplify the transformation scheme such that the time variable is kept real all the way through, but de Hoop's (1960) paper only deals with the case of a source in infinite space. Meanwhile, many authors have used the method to solve a variety of specific problems. We mention the solution of Lamb's problem, where the elastic half space is considered (de Hoop 1961, Gaken-heimer 1969). Some authors only used the two-dimensional version of the method (Achenbach 1973, p.298). Others studied multilayered structures, but either missed de Hoop's modification (Pao and Gajewski 1977), or only solved the scalar acoustic case (Aki and Richards 1980), or introduced ap­ proximations to the method (Wiggins and Helmberger 1974). Of course, there are also many positive aspects to the previously mentioned papers; in fact, we have used them as much as possible; we mention great educational clarity (Aki and Richards 1980), elegant notation (Pao and Gajewski 1977) and a very clear and detailed statement of the results (Johnson 1974, Wig­ gins and Helmberger 1974). In particular, we shall present the solution to the problem in its simplest form and state its concise solution for the all of these encompassing case of an arbitrarily layered isotropic solid.

Still another reason to include the analysis for isotropic media is that we want to present a formalism for wave propagation in stratified media that is generally applicable, i.e. in anisotropic media as well. Therefore, we have avoided concepts that are only advantageous in isotropic media, such as circularly cylindrical coordinates and the wave equations for the scalar and the one-component vector potentials. In the context of this book it becomes clear that Cartesian coordinates and the wave equations for the particle velocities and stresses are the more general ingredients since they are the keystones for the analysis of wave propagation in the anisotropic case. In view of the absence of directional independence in anisotropic media, it is also in general no longer advantageous to introduce circularly cylindrical coordinates (and the corresponding separation of variables) in the analysis of the acoustic waves radiated by a concentrated source in a stratified medium. Since, further, the properties of anisotropic media

(17)

are most easily expressed through Cartesian tensors, a Cartesian reference frame is also the simplest one to express the wave phenomena in. This does not withstand the fact that when studying wave propagation generated by a source in a geometry with circularly cylindrical interfaces, it may still be necessary to employ circularly cylindrical coordinates to satisfy the boundary conditions at the interfaces; the latter type of problem is beyond the scope of our present analysis.

For the same reason of general applicability, we have chosen to use in Chapters 2, 6, and 7 a formalism that applies to arbitrarily anisotropic media. Several authors who have discussed the influence of anisotropy have limited themselves to special symmetry cases like transversely isotropic media (Payton 1983), or to weakly anisotropic media (Booth and Crampin 1983). The range of anisotropy in geophysical applications, however, is not limited to these special cases.

1.3 Numerical considerations

In the numerical treatment of the problem the following steps can be dis­ tinguished: (1) Selection of the generalized rays that have to be included in the calculation, (2) Calculation of the Cagniard-de Hoop contour for each generalized ray, (3) Its use to construct, by inspection, the time-domain Green's function, (4) Convolution of the Green's function with the source pulse to arrive at the complete waveform at the receiver position. It turns out that the numerical methods for this consist of simple algorithms: an eigenvalue procedure to obtain the wave speeds (in anisotropic media), an iterative root-finding procedure to get the Cagniard-de Hoop contour, and the evaluation of a finite-range convolution integral.

The standard argument that is always used to denounce the generalized-ray/Cagniard-de Hoop method is: "In a many-layered model there are far too many rays for efficient computation" (Chin, Hedstrom and Thigpen e 1984). This argument can be refuted by using an appropriate energy-based

criterion in the ray selection procedure.

As a matter of fact, the selection 'procedure of the generalized rays is crucial to the success of the method. First of all, we observe that all gen­ eralized rays are causal functions of time and hence they arrive, one after another, at an observation point. In practice, one is only interested in a

(18)

synthetic seismogram within a finite time interval. This clearly puts, in a particular situation, an upper limit to the number of generalized rays that contribute. But this number may be very large. Secondly, as has been ob­ served by Hron (1972), many of the generalized rays that contribute within the given time window, can be kinematic or dynamic analogs; these rays must be gathered into a single ray with the proper multiplicity. In case the time window is large, still a large number of rays can be involved even after the gathering of analogs. Many of these, however, have already a negligibly small amplitude due to the successive reflections and/or transmissions they have undergone. Therefore, we have, in addition, applied an energy-based criterion, in accordance with which we only take into account those rays that do contribute to the synthetic seismogram on the scale on which the amplitudes can actually be observed. The relevant selection takes place on the basis of a user-defined accuracy parameter. We have found that our energy-based criterion is extremely efficient in the selection procedure of the generalized rays; even in the cases with many layers, the computation time then stays within reasonable bounds.

The whole method as described in this monograph has been imple­ mented in Fortran 77 and runs on a VAX/8600. Numerical results for typical cases are shown. For layered isotropic media they include vertical-seismic-profiling (VSP) geometries. For anisotropic media they include, apart from again VSP geometries, the canonical cases of an unbounded medium and a single layer over a half space. In one case we include pictures of the full wave field as generated by a source, in the form of a "snapshot" in time.

Since the method does employ the generalized ray method, none of the difficulties associated with the frequency-wavenumber integration method (accuracy of incorporating evanescent wave fields) occur. The additional advantages of using the generalized-ray/Cagniard-de Hoop method are: (l) Causality is satisfied automatically, (2) The separate calculation of each causal generalized ray allows an immediate identification of each arrival in the synthetic seismogram, and (3) For simple models with only a few layers, the computation times are extremely short.

Upon comparing the two methods for the computation of acoustic wave phenomena in layered anisotropic media, we observe the following major dif­ ference: a three-fold infinite integral has to be evaluated for the

(19)

frequency-wavenumber integration method, whereas apart from the convolution with the source pulse, a single finite integral and a contour in the complex ray-parameter plane have to be computed for the Cagniard-de Hoop method. In the latter, the time-domain Green's function for a point source then is expressed as a finite integral. Further, for large horizontal or verti­ cal source-receiver separations very efficient approximations exist (Wiggins and Helmberger 1974). These approximations will be discussed as well in our monograph.

(20)

Chapter 2

BASIC RELATIONS FOR

ELASTIC WAVES IN

STRATIFIED, PIECEWISE

HOMOGENEOUS,

ANISOTROPIC MEDIA

Summary

The basic relations for pulsed elastic waves, radiated by sources in strat­ ified, piecewise homogeneous, anisotropic media are presented. First, the decomposition of the wave field in the interior of each homogeneous layer into up- and downgoing wave constituents is studied. After that, the cou­ pling between the wave constituents at the interfaces is accounted for by invoking the appropriate boundary conditions.

(21)

2.1 Introduction

In this chapter we present the general aspects of the radiation of an acous­ tic wave field that is generated by a concentrated source in a stratified medium, each (horizontal) layer of which is a homogeneous and arbitrar­ ily anisotropic solid. The formalism employed is one that enables us to apply later on the Cagniard-de Hoop method by which method space-time-domain expressions for the wave field constituents shall be obtained. In accordance with this method, a Laplace transformation with respect to time, and Fourier transformations with respect to the horizontal Cartesian coordinates are employed.

2.2 Description of t h e configuration and for­

mulation of the problem

We investigate theoretically the pulsed elastic wave motion in a piecewise homogeneous, perfectly elastic, anisotropic solid. To specify position in the configuration, we employ the coordinates {x\,X2,Xz} with respect to a Cartesian reference frame with origin O and three mutually perpendicular base vectors {ii,i2,i3} of unit length each. In the indicated order, the base vectors form a right-handed system. The subscript notation for Cartesian vectors and tensors is used, and the summation convention applies to re­ peated lowercase Latin subscripts; the latter range over the values 1, 2, and 3. The time coordinate is denoted by t. Whenever appropriate, the position is also specified by the position vector x = xp\p. The configura­ tion consists of a number of homogeneous subdomains: the layers, with mutually parallel boundaries: the interfaces. The X3 axis is taken to be perpendicular to the interfaces. This direction is denoted as the vertical direction, and in accordance with geophysical convention the X3 coordinate increases downwardly. The interfaces are then horizontal and parallel to the Xi,x2 plane. Obviously, the configuration is shift invariant in i i and

xi\ it is taken to be time invariant as well. The further nomenclature of

the configuration is presented in Table 2-1 (see also Fig. 2-1). The upper and lower half-spaces are treated in the same way as layers.

(22)

Table 2-1. The layered configuration

Domain Vertical range Layer

£>x - o o < x3 < x3 ; 1 1

D2 x3 ; 1 < x3 < x3;2 2

Dpf-i X3;j\r_2 < X3 < X3;N..i N — 1

(23)

o

1

—^i1

^ 2 © —

D1 - 3

D

2

Layer 2

> > > > > > > > > > > ^ ^ x- = x..

? • • •

< < ( < ( ( < ( ( ( ( ( < < ( < < ( ( ( < ( ( ( ( ( ( ( ( ( ( ( ( < ^

X

3 -

X 3 I N - 2 DN - 1

> > > > > > > > > > > > > > ^ ^ ^ ^ ^

X

3 = * » N - 1

D

N

^Interface N-1

(24)

concurrently the stiffness c,JP, and the compliance 5,; p ?, while in the volume density of mass p,-y we have included a possible anisotropy. The tensors c,;p(7,

Sijpq and pjk must be symmetric and positive definite for a deformation-energy density and a kinetic-deformation-energy density associated with the wave motion in a passive medium to exist. The elastic wave motion will be characterized by the stress Tij and the particle velocity v,-. Further, the action of sources is represented by the introduction of the volume source densities of force ƒ, and of strain rate /i,y. Si-units will be used throughout the analysis.

The source that generates the wave motion is assumed to be located in the interior of one of the homogeneous subdomains. As long as the re­ flections from the interfaces bounding this layer have not arrived at some point of observation in the layer, the wave motion is identical to the one that would be generated as if the source were located in an unbounded medium. Starting with the latter, it is found that the wave field that is causally related to the action of the source can be regarded as a super­ position of six waves, three of which are present in the domain above the source, and these propagate upwardly, the remaining three being present in the domain below the source, propagating downwardly. The amplitudes of the different wave-field constituents are determined by the nature of the source. Once the wave field radiated by the source has been determined, the influence of the interfaces where the configuration changes discontinu-ously in its mechanical properties is incorporated by taking the appropriate reflections at and transmissions across the different interfaces into account. The solution to the wave propagation problem proceeds as follows: (l) First, we present the fundamental equations that govern the wave motion in elastic media. (2) Next, we transform these equations from the space-time domain to the transform domain, using a Laplace transformation with respect to time and. Fourier transformations with respect to the horizontal space coordinates. (3) Taking the vertical traction and the particle velocity as the state quantities characterizing the acoustic wave field, we rewrite, in the transform domain, the resulting equations by introducing a 6 by 6 ma­ trix formulation that is inspired by the propagator-matrix method (see, e.g., Aki and Richards 1980, p.273). (4) Through a linear transformation, we solve for the independent field components using a wave-vector approach. (5) The reflection at and transmission across interfaces is accounted for by the boundary conditions. The latter are the ones that apply to a rigid

(25)

con-tact between two adjacent solids, i.e. the traction and the particle velocity are continuous across the interfaces. (6) We conclude by indicating how the subsequent transformation of the expressions for the field components back to the space-time domain, using the Cagniard-de Hoop method, is carried out. The latter method will be discussed in full detail in Chapters 4 and 5 for isotropic media and in Chapters 6 and 7 for anisotropic media.

2.3 Basic equations for the elastic w a v e m o ­

tion

As the starting point for the description of the elastic wave motion in a ho­ mogeneous subdomain of the configuration we take the linearized equation of motion

diT<3 ~ PijdtVj - - ƒ,-, (2.1)

and the linearized deformation rate equation

\{diVj + djVi) - sijpqdtTpq = hij, (2.2) where the stress rtJ- and the volume source density of strain rate /i,y are

symmetric tensors. The names of the quantities and their Si-units are listed in Table 2-2. Equations (2.1) and (2.2) both have a source term on the right-hand side and only first-order derivatives of the particle velocity and the stress occur in them. The particle velocity and the stress are taken as the two fundamental quantities for describing the elastic wave motion because of the fact that the opposite of their inner product represents the elastodynamic power-flow density.

Since in what follows we repeatedly need the equations from which the stress has been eliminated, we also introduce the stiffness Cijpq as the inverse of the compliance s,-yM, i.e.

CijpqSpqrt = &ijpqcpqra = = 2 V *r J* "■ Oi»0jr)i V*'*v where £,> denotes the symmetric unit tensor of rank two, or Kronecker

tensor. The constitutive tensors satisfy the symmetry relations p,y = pj,,

Sijpq = Spgij, and ctJP, = cp,tJ-; further, they are positive definite. In addition, the stiffness satisfies the symmetry relations

(26)

Table 2-2. Acoustic wave-field quantities and their Si-units symbol[SI-unit] quantity a, [1] cp [ ms x ] cs [ m s- 1 ] cijpq [ Pa ] dii [ s-1 ] «y [i] /i [ Nm-3 ] / o [ N ] hii [ s-1 ] ka [ sm x ] P [ s -1 ] P° [ P a ] Si [ s m- 1 ] sp [ s m- 1 ] ss { s m '1 ] sijpq [ Pa"1 ] Ui [ m ] v,- [ m s- 1 ] V0 [ m V1 ] Xi [m]

unit vector in the direction of point force compressional wave speed

shear wave speed stiffness

deformation rate strain

volume density of body force normalized strength of body force

(we always normalize /o = 1) source strain rate

spatial horizontal Fourier transformation parameter time Laplace transformation parameter

acoustic pressure slowness

compressional wave slowness shear wave slowness

compliance

particle displacement particle velocity

normalized strength of injected volume rate (we always normalize V0 = 1)

(27)

Table 2-2. (Continued)

symbol[SI-unit] quantity

<f>(t) [l] source pulse shape

^>(p) [s] Laplace transformed source pulse shape

6(xi,x2,X3) [ m- 3 ] three-dimensional spatial unit pulse (Dirac function)

8ir [l] Kronecker tensor

K [Pa-1] acoustic compressibility of an ideal fluid A, \i [Pa] Lame coefficients of an isotropic elastic solid

p [ kg m- 3 ] scalar volume density of mass

pij [ kg m~3 ] tensor volume density of mass

rap [ P a ] stress

di [ m- 1 ] partial differentiation with respect to Xi

dt [ s_ 1 ] partial differentiation with respect to t

~ tilde denotes transform-domain variable. Dimension is [m2s] times the dimension in space-time

(28)

similar relations hold for the compliance. With the aid of Eq. (2.3), Eq. (2.2) can be rewritten as

dtTij - Cijpq\{dpVq + dqVp) = -Cijpghpq. (2.5)

In its source-free version, Eq. (2.5) is equivalent to the constitutive relation

Tij = Cijpqepq , where the linearized strain epq =■ \{dpuq + dqup), ut- being the particle displacement, is related to the strain rate dpq = \{dpvq + dqvp) via dpq = dtepq. By using the symmetry relations for the stiffness (cf. Eq.

(2.4)), we can also write Eq. (2.5) as

dtTij - Cijpqdqvp = —Cijpqhpq. (2.6)

It is assumed that />t;-, Sijpq, and c,JM are independent of t; this property reflects the time invariance of the configuration.

At the interfaces, /),-,- and/or s,; p, (c,;p,) jump by finite amounts. Hence, across these discontinuities not all components of v,- and r,;- remain con­ tinuous and therefore these quantities are no longer differentiable. Across the interface the connection between the quantities takes place via bound­ ary conditions of the continuity type. Assuming that the solid layers are in rigid contact, all components of the particle velocity v,- and the three components of the vertical traction r<3 must be continuous.

Only causal solutions of the differential equations (2.1) and (2.2) are acceptable from a physical point of view. Assuming that the sources start to act at the instant t = 0, the causality of the wave motion in the time domain is ensured by putting the values of u,- and r,y equal to zero when

t < 0 .

2.4 The transform-domain equations

The causality condition and the invariance properties of the configuration are most easily taken into account by carrying out appropriate integral transformations. First, we take a one-sided Laplace transformation with respect to time with real, positive, transformation parameter p. To show our notation, we write down the expression for the particle velocity:

f°°

t),(x;p) = / exp(-pt)vi(x;t)dt. (2.7) Jo

(29)

The usefulness of the transformation (2.7) is directly related to the time invariance of the configuration. The lower bound at the integration in the right-hand side of Eq. (2.7) accounts for the vanishing of the disturbance in the configuration prior to t = 0, while the boundedness of the transforms for real, positive, p ensures causality. The latter property is based on Lerch's theorem, which states that upon considering Eq. (2.7) as an integral equation to be solved for u,- if i>i is known for real, positive p, the solution of this integral equation (which is also known as Carson's integral equation) is unique, while it vanishes for negative values of t (Widder 1946, p.63). Now, the condition of causality is replaced by a condition of boundedness everywhere in space, of the transformed quantities that characterize the wave motion.

Further, to take advantage of the shift invariance of the configuration in the horizontal direction, we introduce the Fourier transformations with respect to the horizontal coordinates x\ and x2 (with transformation pa­ rameters pki and pk2, respectively, where kt and k2 are real). For the particle velocity we have

/

oo roo

/ ex-p(ipk(rxa)vi(-x.;p)dxidx2, (2.8) -oo J — oo

with the inverse transformation

/

OO /"OO

/ exp(-tpkaxa)vi(iki,ik2,x3;p)dk1dk2, (2.9) -oo J—oo

where Greek subscripts are used to indicate the "horizontal" parts of vectors and tensors, i.e. they are to be assigned the values 1 and 2 only.

For later convenience we also introduce the variables sa, which will be denoted as the "horizontal slownesses," with sa purely imaginary, defined through ika = sa. In terms of these, Eq. (2.8) changes into

/

oo roo

/ exv(psaxa)vi(x;p)dxidx2, (2.10) -oo J—oo

and Eq. (2.9) into

/

too rioo

/ exp(-psaxa)vi(si,S2,x3;p)dsids2. (2.11) -too J— «oo

(30)

The relation between sa and the real-valued slownesses that characterize the propagation of uniform plane waves in the medium, will become clear later (cf. Appendix A).

Applying the rule that dt^+p and da—*■ — psa, the application of the transformations indicated in Eqs. (2.7) and (2.10) to Eqs. (2.1) and (2.6), leads to

d3Ti3 = pptjVj + psaTia - ƒ,-, (2.12) and

CijpzdzVp = pCijpi/SvVp + pfij + Cijpqhpq. (2-l3) From these equations it is clear that only the differentiation with respect

to z3 is left. Hence, we can, in the interior of each homogeneous layer, use standard methods from the theory of ordinary differential equations with constant coefficients to solve for the wave field in the transform domain. This is done by introducing an appropriate matrix formalism.

The boundary conditions at interfaces are replaced by their transform-domain counterparts, viz. the continuity of t/,- and f,3. Further, the transform-domain quantities must remain bounded as £3—► — 00 and as x3—>oo.

2.5 T h e motion-stress vector in a homoge­

neous subdomain

Equations (2.12) and (2.13) will be rewritten in a 6 by 6 matrix formalism in which the motion-stress vector occurs as the quantity that characterizes the transform-domain elastic wave motion. A similar method is used in the propagator-matrix method that is commonly used in geophysics (see, e.g., Aki and Richards 1980, p.273, Kennett 1983, p.40). Our aim is, however, different: we prefer to write the total wave motion as a superposition of propagating waves, and in accordance with this, the transform-domain so­ lution will be cast in a form where the generalized-ray picture is manifest (Spencer 1960). The motion-stress vector is chosen as the 6 by 1 matrix

(31)

where v,- = {vi,ï>2,v3)T and f,3 = (ri3,f23»^33)T- Here T means transpose. Note that the field components in bj are the ones that are continuous across a horizontal interface between two elastic media in rigid contact. In our matrix formalism uppercase Latin subscripts take on the values 1 to 6, and the summation convention applies to them over this range. Our goal will be to write Eqs. (2.12) and (2.13) as a system of first-order ordinary differential equations of the form

d3bj = -pAubj + FIt (2.15)

where Au is the 6 by 6 elastodynamic system's matrix and Fi is the 6 by 1 source matrix or source vector. Hence, the quantities raii must be eliminated from the equations. Since the matrix formalism will only be used in the transform domain, we have, for notational simplicity, omitted the tilde over the relevant matrices.

We start the elimination process by decomposing Eq. (2.13) into the two equations

Ci3p3d3vp = pci3pvsuvp + pri3 + ci3pqhpq, (2-16)

Ciopq hpq •

(2.17) Solving Eq. (2.16) for d3vp we get

Ö3Vk = {c.3.3)k?{pcr3pt,s,,vp + pfr3 + cr3pqhpq), (2.18) where (c.3.3)^.1 denotes the 3 by 3 matrix that is the inverse of the 3 by 3

matrix cr3p3, i.e.

[c.3.3)krCr3P3 = Skp. (2.19)

It is important to note that (c.3.3)^.1 differs from Sfc3r3. Substituting Eq. (2.18) in Eq. (2.17) we obtain

ha = -QiopvSvVp + Ci(,k3(c.3.3r fr 3 - p~ Qiapqhpq, (2.20)

where

Qiapq = Ciopq ~ C,<7A:3(c.3.3)jtr Cr3pq. ( 2 - 2 l )

Now substituting Eq. (2.20) in Eq. (2.12), we arrive at

(32)

Equations (2.18) and (2.22) have the structure that is required in Eq. (2.15) with the system's matrix given by

(2.23)

Sip Tir

in which

and

and the source matrix

in which and

Tkp = -«„(c.3.3) krCr3pu, (2.24)

Ckr = (c.3.3)^1 = Crk, (2.25)

Sip = pip — sasvQiopt, = Spi, (2.26) T? = -saci<7k3(c.3.3)ïr\ (2.27)

Fi = ( FJr ) , (2-28)

** = (c.3.3)^cr3p,fcM, (2.29)

F- = fi + saQiapqhpq. (2.30)

It is clear that Tkp,T?,Sip, and Ckr are 3 by 3 partitions of Au. In view of the symmetry properties of pi}-, s,yp,, and c,ypg, the matrices 5,p and Ckr are symmetric.

The structure of AJJ reflects the operations that were used to derive it. Observe that the fundamental equations (2.1) and (2.2) only show a cross-coupling between the field quantities v,- and rtJ-. This cross-coupling is still present in the anti-diagonal 3 by 3 parts of the matrix Au. The diagonal parts of Au are due to the process of elimination of the stress components Tag. Finally, the elimination process has also complicated the anti-diagonal part 5,-p with the term containing Qiapv, and it has changed the compliances from Eq. (2.2) into the expression for Ckr in Eq. (2.25).

(33)

Equations similar to Eqs. (2.14), (2.15), and (2.23) - (2.27) have also been derived by Fryer and Frazer (1984). Since these authors started out with a Fourier transformation from the time domain to the frequency do­ main, they have dt-+ — iw, whereas we have dt—*p . Another difference is that our motion-stress vector bj is composed of the particle velocity and the stress, while Fryer and Frazer (1984) use a motion-stress vector composed of the particle displacement and the stress. The advantage of our choice over the one of Fryer and Frazer is that there is no factor of "p" difference between the first three and last three components in motion-stress vector

bj; this will turn out to be an advantage in the transformation back to the

time domain.

The system's matrix Au is a function of the horizontal slownesses si and 52, the volume density of mass pij, and the stiffness c,;pg only. It is independent of the Laplace transformation parameter p and, since the coefficients in the differential equations (2.1) and (2.2) are real-valued, it is real-valued for real values of si and 62. (At the moment, Si and s2 are still imaginary.)

2.6 T h e transform-domain w a v e vector in a

source-free domain

To elucidate the structure of the wavelike solutions of Eq. (2.15), we carry out a linear transformation on the motion-stress vector bj and, through it, want to arrive at a wave-vector formalism in which the presence of up- and downgoing waves is manifest. Let WN be the wave vector that is related to the motion-stress vector via the linear transformation

bj = DJNwN, (2.31)

in which DJN is subject to a convenient choice. On the assumption that

Djff is non-singular, its inverse Dj}j exists and the relation inverse to

Eq. (2.31) is

wM = Dtttbj. (2.32)

In the homogeneous medium we assume DJN to be independent of X3 and the substitution of Eq. (2.31) into Eq. (2.15), followed by the

(34)

premul-tiplication by D^i yields

d3wM = -PAMNWN + Dj£jFIt (2.33) where the 6 by 6 matrix AMN is given by

AMN = D^AUDJN. (2.34)

The desired structure of Eq. (2.33) is now arrived at if AMN is diagonal. From the theory of matrices it is known that this is accomplished by taking

DJX to be the eigencolumn matrix of Au; in this case, A^^r is the diagonal

matrix of the corresponding eigenvalues of AJJ. The relevant eigenvalues follow from the determinantal equation

det{Au - \{N)Iu) = 0, (2.35)

where IJJ is the 6 by 6 identity matrix, and A ^ is the TV th diagonal element of AMJV. Equation (2.34) is equivalent to AIJDJN — -D/MAA/JV and to DMJAU = AMND~^J. Hence, we can write either

AlJbW = \Wbf\ (2.36)

or

g^Au = \MgM, (2.37)

where b^ is the N th column of the eigencolumn matrix DJN of AJJ , and

g\ ' is the M th row of the eigenrow matrix D^ of Ajj. (The vector g\M' is also denoted as a reciprocal, or left, eigenvector of Au-) In view of the fact that

DMIDIN — IMN, (2.38)

the normalization of b\ ' and g\ ' must be taken such that

rf">6<"> = 6MN. (2.39) The simple way to 'find D^j from DJN without resorting to numerical

inversion, is described extensively by Fryer and Frazer (1984). We only summarize the procedure here. Each column of DJN is an eigenvector of the form

(35)

The algorithm for computing the eigenrow matrix Df}1 is simply: (1) Con­ struct each of the six columns of the raw eigencolumn matrix somehow to form DJN, (2) interchange the traction and velocity 3-vectors ( —T}3 ' and v] ' ) in each column, (3) transpose the result, and normalize each of the six rows g\ ' of the resulting matrix in accordance with Eq. (2.39). Then, the latter matrix is the eigenrow matrix D^}T .

Since the equations in Eq. (2.33) are now mutually uncoupled, their solutions in a source-free domain are the six linearly independent functions

WN = constant X e x p ( — p \ ^ x3) , (2.41) each of which has the structure of a wave, propagating in the vertical di­

rection. For the imaginary values of si and s2 that are considered thus far, the eigenvalues can be shown to separate into two sets of three, viz. { • s j1, ^2, ^3} , whose real parts are negative: they represent the "vertical slownesses" of the upgoing waves (see Fig. 2-2), and { s j1, ^2, ^3} , whose real parts are positive: they represent the "vertical slownesses" of the down-going waves (for the proof, see Appendix B). Choosing a particular ordering we may write

AMJV = d i a g f o1, ^2, ^3, ^1, ^2, ^3) , (2.42) or, alternatively,

KMN = diag(A<^A<^ A<^ A«,AM, A « ) . (2.43) At this point we want to motivate why we use two different symbols to

denote the six eigenvalues of AJJ, and under what circumstances each of the notations will be used. When we focus on the 6 by 6 matrix formulation and the properties that are derived from linear algebra, we shall write A W where JV € {l,..., 6}, in the ordering indicated in Eq. (2.43). However, when we focus on the eigenvalues as the vertical slownesses for three upgoing and three downgoing waves, we shall write s^n, where n G {1,2,3}, in the ordering indicated in Eq. (2.42). The correspondence between the two sequences is inferred from Eqs. (2.42) and (2.43).

Whenever advantageous, the superscripts {—1,-2,-3} will henceforth generally be used to denote wave fields propagating in the direction of decreasing x3 (upgoing waves), while the superscripts { + l , + 2 , + 3 } will

(36)

•a"

1

S - 2

3

downward propagating ) i s - 3

waves

3

* * 3

s

3

upward propagating

+ 1 I ^ waves

s + 2

3

S+3

3

Fig. 2-2. The six admissible vertical slownesses 53fcn. Positive superscripts denote downgoing waves, negative superscripts denote upgoing waves.

(37)

generally denote wave fields propagating in the direction of increasing X3 (downgoing waves). The elements of w^ will be ordered accordingly, i.e.

wN = {w~,w+)T, (2.44)

where,

w- = {w-\w-2,w-3)T, and u;+ = {w+\w+2,w+3)T. (2.45) After ordering AMN according to Eq. (2.42), the ordering of the columns

and rows of the matrices DJN and Dj£j is fixed.

For isotropic media, and for certain classes of anisotropic media as well (the latter if the X3 direction is perpendicular to a symmetry plane of the medium), we can pair the vertical slownesses such that s$n = — s^", but for arbitrarily anisotropic media no such simple relationship exists. Therefore, the identification as to which waves are upgoing and which are downgoing is more complicated in anisotropic media. For a detailed discussion on how we sort the six eigenvalues A ^ into the six ordered vertical slownesses sfn we refer to Section 6.3.

Concluding this section we see that the matrix Au has two convenient properties: its eigenvalues are the six admissible vertical slownesses sfn for a given {51,52} pair (51 and 52 imaginary) in a homogeneous subdomain, and there is a simple relationship between the associated eigencolumn and eigenrow vectors, which allows the calculation of the inverse Dj^j of the eigencolumn matrix DJS, without using numerical inversion.

2.7 T h e transform-domain w a v e vector and

t h e motion-stress vector generated b y a

localized source in a homogeneous

sub-d o m a i n

We now proceed to determine the transform-domain wave motion in a ho­ mogeneous domain in the presence of sources. We assume that the source vector Fj in Eq. (2.33) contains concentrated sources that act at a single level only. Distributed sources can be included in the analysis by consider­ ing them as a continuous superposition of concentrated sources. Without

(38)

loss of generality, we can assume that the source is located at the point (0,0,X3;j!). The source vector is then of the form

FI = 4>{p)XI6(x3-x3.t>), (2.46) where S(x3 — x3;,) is the one-dimensional unit pulse (Dirac delta function)

acting at the source level x3 = x3.t, <j>(p) is the Laplace transform of the source pulse shape (signature of the source), and Xj is a vector that only depends on the nature of the source (e.g. an explosion source or a point force; examples of these will be discussed later on).

The wave vector WM of the wave field generated by the source must then be determined from (cf. Eq. (2.33))

d3wM = -P^MNWN + UP)WMS[X3 - x3 ;,), (2.47) where the amplitude vector WM , defined by

WM = (W-,W+)T = DM\Xj, (2.48)

contains the "excitation coefficient" for each type of wave. The latter ex­ presses how strongly each of the six waves is coupled to the source. Causal­ ity entails that above the source WM must be bounded as x3—► — 00 and hence only the waves with Re(5^ri) < 0 (i.e. {s3n }) are admissible, whereas below the source WM must be bounded as x3—►oo and hence only the waves with R e ( s ^ ) > 0 (i.e. { 4 " }) are admissible. In accordance with this, we have arranged WM as indicated in Eqs. (2.44) and (2.45). Now, we construct the "one-sided" modified amplitude vector WM = {W~,W+)T as

W-{x

3

) = {

and ^:(x3) = \ - e x p [ - p s3" ( z 3 - x3.tt)]W- if - 0 0 < x3 < x3;„ 0 if x3;, < x3 < 00, (2.49) 0 if - c o < x3 < x3i,, , e x p [ - p 4n( x3 - x3.a)\W+ if x3 ; j < x3 < 00, (2.50)

in which the propagation factors that occur in Eq. (2.41) axe included. Further, we know that in order to satisfy Eq. (2.47), the first derivative

(39)

of WM with respect to x3 must contain a Dirac delta pulse at the source level x3 = x3.t,. Hence, WM itself must, at the source level x3 = z3 ; s, have a positive step of magnitude

lim WM — lim WM = 4>{P)WM- (2-51)

From Eqs. (2.49), (2.50), and (2.51) we conclude that the transform-domain wave vector, generated by a localized source in a homogeneous medium is given by

WM{X3) = UP)WM(X3). (2.52)

By combining Eq. (2.31) with Eq. (2.52) we arrive at the final expression for the motion-stress vector in the transform domain as

M5i > 52, S3; P) = ^>{p)DjM{su S2)WM{X3). (2.53)

The right-hand side of Eq. (2.53) is a superposition of six waves (of which, in agreement with Eqs. (2.49) and (2.50), either the three upgoing waves or the three downgoing waves are zero, depending on the position of the point of observation relative to the source position). In Eq. (2.53), the propagation factor that occurs in WM can be written as exp[—PA.MN{^3 — x3;3)}; for our later method of transforming back to the time domain we need the properties of the different vertical slownesses explicitly, and therefore we shall use s3n instead of AMN>

To illustrate the transformation of the right-hand side of Eq. (2.53) back to the time domain, we shall discuss the details of this transformation for a single wave constituent. The expression for the motion-stress vector of an individual wave constituent has the general shape

6j(5i,s2,x3;p) = <j>{p)Bj(s1,s2)exp[-p(x3 - x3;a)s3(si,s2)}, (2.54) where Bj is independent of p and where, for notational simplicity, we have written 53(^1,s2) = s3n(si,s2). The six different expressions of Bfn, one for each of the wave constituents, follow from Eq. (2.54) and the definitions of DJS and WM as

B±n = ±b±ngfnXI, (2.55)

where bjn and gfn are the right eigenvector and the left eigenvector of AJJ, respectively, corresponding to the eigenvalue sf". The method of trans­ forming the right-hand side of Eq. (2.54) back to the space-time domain is indicated in Section 2.10.

(40)

2.8 The transform-domain w a v e vector and

the motion-stress vector generated by a

localized source in a stratified m e d i u m

In this section we shall derive an expression for the transform-domain wave vector and the motion-stress vector generated by a localized source in a stratified medium. The final result will have a structure similar to Eq. (2.53), but of course more complicated, because at each interface between the layers the boundary conditions have to be satisfied. The added com­ plexity due to the stratification is reflected in the extra indices for layers and interfaces in the expressions. The symbols can get quite awkward, but fortunately they are only used in this section and the next one. Once it has been established that the total solution can rigorously be expressed as a sum of generalized-ray wave constituents (which will be defined later in this section), and once it is clear that there can be no confusion about reflec­ tion and transmission coefficients and the trajectory of each generalized-ray wave constituent, the extra indices will be omitted.

In the interior of each homogeneous layer we separate the wave motion into the up- and downgoing waves. In each layer, the relation between the motion-stress vector and the wave vector is, in agreement with Eq. (2.31), given by

*i? = M W , (2-56)

where the superscript t applies to the relevant layer (see Fig. 2-3). In this layer, with domain x3;,_i < x$ < x^i, the wave vector (cf. Eq. (2.44)) is written as

« # ( * 3 ) = [Wn]~ expf+pS^Zs;.- ~ X3)], U>£)+ e x p [ - j w j ? ( x3 - *3;i-l)]) ,

(2.57) where we have used the convention that in every layer of the medium the waves inside that layer have zero phase at the interface at which they orig­ inate.

The boundary conditions at interface *', located at x3 = x3;t-, require the continuity of the six elastic field components in the motion-stress vector, i.e.

(41)

N3 ; i - 1

layer

w'

(0+

4 w

( ö -

n( i ) UJ N ^i — x3',i "" x3;i

-i + 1

w

0 +

1)-W (1 + 1 ) -D

(i

+ D

UJ N

x

3;j

interface i

h,

+

.

x3;i + 1

Fig. 2-3. Notation for layers and interfaces for a stack of homogeneous layers.

(42)

This can be rewritten, by using Eq. (2.56) and by suppressing the trans­ formation parameters p, Si, and S2 in the argument, as

lim D%w^(x3) = lim D ^ w ^ (x3). (2.59)

In Eq. (2.59) the wave vector values are given by (cf. Eq. (2.57))

lim wjp(xs) = (u»W-,u/W+ e x p ( - p 4;^ , ) )T , (2-60) and

lim i 4+ 1 ,( *3) = («tf+1)- exp(+P537+ 1/i,+ 1), utf+ 1>+)T , (2.61) in which only the thicknesses of each layer, given by hi = x3-,i — a^ji-i, occur.

Next, we apply the scattering formalism at interface i. This means that we express the amplitudes of the waves propagating away from the interface in terms of the amplitudes of the waves propagating towards the interface (see Figs. 2-3 and 2-4). We write

«ff" = * f f i > ?, + « p ( - j < f c ) + Tj^wt1]- exp(+ps;»+1hi+1), (2.62) and

™£+1)+ = T j ^ + e x p t - p ^ (2.63)

where R$* and Tj^ are the 3 by 3 reflection and transmission coefficient matrices that characterize the interactions at interface i. The definition of the reflection and transmission matrices is shown schematically in Fig. 2-4. For example, R^x is the amplitude of the upward wave in layer i with vertical slowness s^?{, generated by reflection from interface i of a downward wave of unit amplitude in layer i with vertical slowness s^\. The values of the reflection and transmission coefficients are found by substituting Eqs. (2.60)-(2.63) in Eq. (2.59), and requiring that the resulting equations hold for arbitrary values of w£+1)~ and w$+. The result is written in the form of the scattering matrix S$N, containing all reflection and transmission

(43)

>(i) +

< i )

-interface i

interface i

,(0-Fig. 2-4. Definition of reflection and transmission matrices at interface i.

R(%)+ and T^+ describe the reflection and transmission of a wave initially travelling downward; flW- and T®' of a wave initially travelling upward.

(44)

coefficients for scattering a t t h e interface i:

c(») _

DMN —

(2.64) where Q^n are the 3 by 3 partitions of t h e 6 by 6 m a t r i x QMN t h a t is given

by / Q1* Q12 *v mn mn QMN = = {DMyylD%. (2.65) \ Q21 Q22 \ T» mn ™ mr

For imaginary values of s^, the matrices ( % j )_ 1 a n <* (Qmn)"1 prove t o be nonsingular, which implies that the inversion can b e carried out.

At t h e level x$ = x^&, in the interior of t h e layer Ds, a concentrated

source is located. For the wave vector of the wave motion generated by this source, t h e expression from Eq. (2.52) holds. It will be convenient to position an "artificial" interface at t h i s source level. A t this interface 5 there is a perfect transmission and n o reflection and hence t h e reflection and transmission matrices, defined in E q s . (2.62) and (2.63) reduce t o

J2W± = 0 , « n d T ^ = Imn, (2.66)

where Imn is t h e unit m a t r i x of rank t h r e e . Incorporating the source terms

in accordance with Eqs. (2.49), (2.50), and (2.52) we have

u/W- = Imnwlt+V- exp(+ps^,+1h,+1) + HP)W-(X3), (2.67)

a n d

^+ 1 ) + = U ' 1 + exp(-ps£hs) + HP)W+(X3). (2.68) Combining these results with the p r o p e r t y t h a t in t h e t w o half-spaces

extending to X3—► — 00 and £3—►00, only waves propagating in t h e direction of decreasing a n d increasing x3, respectively, are present, we c a n determine all transform-domain elements of t h e motion-stress vector i n each layer (cf. E q s . (2.62), (2.63), (2.67), (2.68), and (2.56)). I n the transformation back t o the space-time domain, as well as for the interpretation of the results, it will be advantageous to introduce t h e concept of generalized ray and to write t h e total wave motion as t h e superposition of generalized-ray wave constituents. This is discussed in the next section.

(45)

2.9 The transform-domain generalized-ray

wave constituents

Following Cisternas, Betancourt, and Leiva (1973), it will be shown that the full wave field excited by a localized source in our stratified, anisotropic medium can be expressed as the sum of generalized-ray wave constituents. To this end, we first define the modified reflection and transmission coeffi­ cient matrices R$* and T&}*, by including in them the propagation factors that occur in Eqs. (2.62) and (2.63), which are rewritten as

«"ff" = A2k

+

«tf

,+

+ 2ffi-«i'

+ 1 )

-, (2.69)

and

« £

+ 1 , +

= Jï&

+

«tf

,+

+ i & U i

, + 1 )

- . (2.70)

Equations (2.67) and (2.68) are rewritten in a similar way. Next, we arrange in one large column vector X the wave amplitudes of the different waves in all layers. Then, Eqs. (2.69), (2.70), (2.67), and (2.68) can be rewritten as ( • \ ( ■ \ I ■ \ m V

J

0 ■"-ran T(i)+ 0 0 0 0 0 0 0 R(i+1)+ 0 ■'■'■mn 0 V w

y v

while at the source level we have

f : \ f -w(«+i)+ tu. (•+i)-/ 0 0 run 0 0 0 0 0 0 0 0 0 0 *mn 0 0 J

+

(2.71) 0 o m + V In a more compact way, Eqs. (2.71) and (2.72) are replaced by

x = nx + s,

(2.72)

(46)

where Ï1 contains the modified reflection and transmission matrices, and S is a column vector containing the source. Equation (2.73) can be rewritten as

X = ( I - n )_ 1S , (2.74)

which leads to the power-series expansion

oo

X = £ n « S . (2-75) It remains to be shown that this power series converges.

Now, from Eq. (2.73) it follows that

Q Q Q

£n«s = £n«x - £n«

+1

x

= x-n«

+1

x,

or Q X = £ n ' S + n g + 1 X- (2.77) ?=o

Further, in practice we are only interested in the wave motion in a certain finite time window. However, from some finite value of Q onward, the factor

ftQ+1 leads to such a delay (note that none of the exponential functions occurring in Jl has a vanishing argument) that fI<'+ 1X yields in the time domain only contributions that start after the time window of interest has elapsed.

A nice property of this formulation is its physical interpretation. We have seen that S has only non-zero elements at the level of wfi~ and wj*+1)+ in Eq. (2.72). Hence, the term corresponding to q = 0 gives fi° = I. This means that the first term of the series gives X = S, which, in agreement with the transform-domain wave vector in a homogeneous medium from Eq. (2.52), gives the direct waves in the form

(2.78) « £+ 1 , + (*3) = HP)W^(X3),

(47)

i.e., the downgoing waves in layer 5 + 1 and the upgoing waves in layer s. The term corresponding to q — 1 gives Q1 = fï, which means that this term gives the once reflected and transmitted arrivals as shown in Fig. 2-5, i.e., it gives the complete solution in case multiples are neglected. The term corresponding to q = 2 generates the arrivals that have two interactions at interfaces. For a receiver located in layer s + 3 (see Fig. 2-5) the first arrival will be due to the term corresponding to q = 2. For a receiver located in layer 5 + 1, we see that the the term corresponding to q = 2 gives rise to two arrivals with different paths already.

Equation (2.75) expresses that in each layer the wave amplitude is a sum of a number of terms, each of which (due to the summation over q and the matrix multiplications fi') consists of a product of modified re­ flection/transmission coefficients, in a unique order, multiplied by a source term. Each term in this summation can be seen as a wave constituent with a unique trajectory, determined by interactions at interfaces and propagation through layers. Such a term is called a generalized-ray wave constituent or, in short, a generalized ray (Spencer 1960). The concept of ray is usually connected with high-frequency approximations or ray asymptotics (Felsen and Marcuvitz 1973, p.123); therefore, it is important to note that the sum of generalized rays in Eq. (2.75) gives the full exact solution to the wave equation, and is just one of the representations of the full solution.

Since in each anisotropic layer there are three upgoing and three down-going waves, the number of generalized rays that might contribute to the wave field at a particular point of observation can be very large. For the example in Fig. 2-5, the term corresponding to q = 0 contains 3 generalized rays in layer s, the term corresponding to q = 1 contains 9 generalized rays in layer s, and the term corresponding to q = 2 contains 54 generalized rays in layer s. In Chapter 7 we s'hall discuss how these contributions can be efficiently combined.

As a consequence of the generalized-ray decomposition, we see from Eq. (2.75) that the wave amplitude in the receiver layer can be written symbolically as

„(receiver)* = ^ £(J J ££)*) ^(source)* ? (2.79)

(48)

q = O s + 1 q = s - 1 s + 1 q = 2 s - 2 s - 1 s + 1 s + 2

Fig. 2-5. The generalized rays that are included in the terms with q = 0,1,2 for a structure with a source located at interface a.

(49)

or

^(receiver)* = fa) ^(R R$±) eXpfaZ 4^)^°^> (2-80)

q i wi,n

where I l i i E ^ is the product of q interface reflection/transmission coeffi­ cients, J2m,n ^m^m" contains q terms, and /i*n is the total vertical path length of wave ± n in layer m. By comparing Eq. (2.79) with Eq. (2.52), we see that they have a similar structure, where the latter is just the special case of the generalized rays that are associated with the direct arrivals in an unbounded anisotropic medium.

Upon combining Eq. (2.56) with Eq. (2.80), we arrive at the final transform-domain motion-stress vector in the receiver layer of a stratified medium

b^"\sus2;p) = D{J^Wei)w^eWet). (2.81) The total wave motion represented by the motion-stress vector bj can be

split into the contributions from t h e individual generalized rays by just taking one individual term in the right-hand side of Eq. (2.80). Note that a particular value of q gives rise t o a number of different terms (gener­ alized rays). The expression for t h e motion-stress vector due to a single generalized ray is written symbolically as

bj(si,s2;p) = <j>{p)Bj(s1,s2)exp[-pYï,S3;m{si,S2)hm], (2.82)

m

where Bj is independent of p, and where, for notational simplicity we have written

Es*AsuS2)hm = ±"£s^m(sus2)hin. (2.83) The method of transforming the right-hand side of Eq. (2.82) back to the

space-time domain is schematically indicated in the next section.

2.10 Transformation of the solution back t o

the space-time domain

The transformation back to the space-time domain is carried out by us­ ing a particular version of the Cagniard-de Hoop method. With the aid

Cytaty

Powiązane dokumenty

Therefore, Theorem 4.3 may be generalized to all line graphs of multigraphs which possess maximal matchable subsets of vertices – for example, the line graphs of multigraphs

In particular, compact convex sub- sets of R n with nonempty interior, fat subanalytic subsets of R n and sets in Goetgheluck’s paper [G] (where a first example of Markov’s

(Given a Steiner tree for a set of k vertices, one possible closed walk through those vertices would trace each edge of the Steiner tree twice.) The k-Steiner distance plus one

By Lemma 3.2 this maximum is attained at a point all of whose coordinates, with possibly one exception, lie on the unit circle... The assertion of our theorem

(This is trivial for s = 1 and was shown for s = 2 in [11].) In this connection (especially for numerical integration) the notion of good lattice points plays an outstanding

It is also remarked there that this fact is a consequence of a lemma of [11] which in turn is proved via Kloosterman sums and Kuznetsov’s trace formulas.. We shall prove Lemma 3

If X is a real Hilbert space condition (d) can be replaced by “F ( · , x) has a strongly measurable selection” and the values of F need only be closed convex.. This is Theorem 10.5

More- over, our results and methods used in the proof suggest that in the class of bounded pseudoconvex complete Reinhardt domains the symmetry of the Green function is equivalent