An Efficient Approach to the Seismogram Synthesis for a ... - CiteSeerX

Bulletin of the Seismological Society of America, Vol. 86, No. 2, pp. 379-388, April 1996

An Efficient Approach to the Seismogram Synthesis for a Basin Structure Using Propagation Invariants b y H. Takenaka, M. Ohori, K. Koketsu, and B. L. N. Kennett

Abstract

The A N - L a m e r method is one of the cheapest methods for synthetic seismograms in irregularly layered media. In this article, we propose a new approach for a two-dimensional S H problem, solved originally by Aki and Lamer (1970). This new approach is not only based on the Rayleigh ansatz used in the original AkiLamer method but also uses further information on wave fields, i.e., the propagation invariants. We reduce two coupled integral equations formulated in the original AkiLarner method to a single integral equation. Applying the trapezoidal rule for numerical integration and collocation matching, this integral equation is discretized to yield a set of simultaneous linear equations. Throughout the derivation of these linear equations, we do not assume the periodicity of the interface, unlike the original AkiLarner method. But the final solution in the space domain implicitly includes it due to use of the same discretization of the horizontal wavenumber as the discrete wavenumber technique for the inverse Fourier transform from the wavenumber domain to the space domain. The scheme presented in this article is more efficient than the original Aki-Larner method. The computation time and memory required for our scheme are nearly half and one-fourth of those for the original Aki-Lamer method. We demonstrate that the band-reduction technique, approximation by considering only coupling between nearby wavenumbers, can accelerate the efficiency of our scheme, although it may degrade the accuracy.

Introduction we further develop the Aki-Larner method by using the propagation invariants that were first introduced by Kennett (1984). A general form of the propagation invariants for time harmonic elastic fields in irregularly layered media was shown by Kennett (1984), which is closely related to Betti's theorem. The propagation invariants depend on the interaction of two distinct wave fields and vanish if both fields are identical. Haines (1988) extended the concept of propagation invariants to the spatial wavenumber domain for acoustic wave fields in laterally varying media. Subsequently, Kennett et al. (1990) used the propagation invariants to describe elastic wave fields in anisotropic, laterally varying media. Koketsu et al. (1991) exploited the propagation invariants in the coupled wavenumber domain to simplify the representation of the reflection and transmission of elastic waves in irregularly layered media. Recently, Takenaka et al. (1993) derived a number of new forms of propagation invariants for elastic waves in irregularly layered media, which are based on the transformation properties of integral operators represented in terms of displacement and traction components. Here we use a form of the propagation invariants presented by them.

The Aki-Larner method (Aki and Lamer, 1970), based on the Rayleigh hypothesis (e.g., Millar, 1973) and the discrete wavenumber representation, is one of the widely used methods for synthetic seismograms in irregularly layered media. In the AN-Lamer method, the wave field in each layer is expressed as a superposition of plane harmonic waves, including inhomogeneous plane waves, and the boundary condition is matched in the horizontal wavenumber domain by taking advantage of the fast Fourier transform (FFr) in space. The method has been extended to the time domain by Bouchon (1973) and Bard and Bouchon (1980), to the case of multiple layers by Kohketsu (1987a) and Horike (1987), and to the case of vertically inhomogeneous layers by Bard and Gariel (1986). Recently, the Aki-Larner method has been extended to solve three-dimensional (3D) problems with two-dimensionally irregular interfaces by Horike et al. (1990), Ohori et al. (1990), and Uebayashi et al. (1992). Since the AN-Lamer method has the advantage that it requires less memory and less computation than most of the alternative numerical methods (e.g., Axilrod and Ferguson, 1990), it may be one of the cheapest methods. In this article, 379

380

H. Takenaka, M. Ohofi, K. Koketsu, and B.L.N. Kennett

In this article, we consider a two-dimensional S H problem for the simplest configuration of irregularly layered media, which consists of a half-space and only one sedimentary layer over it. Although this problem, solved originally by Aki and Lamer (1970), may now be too simple for seismologists to investigate seismic wave propagation in realistic media, this problem is the best to illustrate the development of the A N - L a m e r method based on the propagation invafiants and to demonstrate the feasibility of our approach. Of course, the present method can be extended to the P - S V system, to the 3D system, and to multi-layered systems, which are also outlined in this article. Formulation of Problem In this section, we show a formulation of a two-dimensional S H problem for a basin structure that consists of a half-space and a sedimentary layer over it. The formulation described in this section is essentially the same as that of Aki and Lamer (1970). We take a Cartesian coordinate system [x,y,z] with zaxis taken positive downward and assume that the free surface is horizontal and the interface varies in the x-z plane. Therefore, the wave field does not depend on the coordinate y. We assume a time dependence exp( + icot) but will not normally represent the time variation explicitly. When we take the Fourier transform with respect to x of some functionf(x,z), as follows:

f(k,z) = f~_ f(x,z) e x p ( - i k x ) dx,

(1)

in a homogeneous and isotropic layer having the S-wave velocity fl and the density p, the y-component of the displacement in the k (horizontal wavenumber) domain can be written as a harmonic solution with weighting factors H v and He:

9(k,z) = HUgH e x p ( + iv~z) + HDgH exp(-- iv~z),

(2)

z=O

A

_ VDA

VUA

VUB

VDB

B z--H

Figure 1. Configuration of calculation. Cartesian coordinate system [x, y, z], with x- and z-axes taken positive right and downward, respectively, is employed. A homogeneous layer (basin) A and a halfspace (bedrock) B are separated by an irregular interface Sx. The free surface is horizontal. The source denoted by a star is located below St and generates the upgoing wave field fsu at a level z = H between Sz and the source. which fluctuates around the reference level z = d (Fig. 1). We assume that the free surface is horizontal (z = 0), that a source is located below the deepest level of St, and that there is no source in layer A. When we define the upgoing and downgoing waves

Vv = H u e x p ( + ivlfl),

Vo = Ho e x p ( - iv~d),

(5)

respectively, the displacement in layer A can be represented in the spatial domain as the sum of upgoing and downgoing waves of all k's:

V(X,Z) = ~

1;

-1- ~

VuA(k) exp[ + iV~A(Z -- d)] exp( + ikx) dk

YDA(k) exp[-- iV~A(Z -- d)] exp(+ ikx) dk.

(6)

where Yfl = ifo2/fl 2 -- k2) 1/2,

Im(D) _--- - h(x)],

(4)

This is the same as equation (1) in Aki and Lamer (1970), except for the use of different notations, and valid over all of region A, even near the interface. The displacement on the interface is then expressed as

v[x,z = d + h(x)] = (~Y[UA~UA) iX) -I- i3VfDAYDA) ix),

(7)

where the integral operators Mu, o have been defined as

if

(:34u,DVv.D) (X) =- ~

exp[ + iv~h(x) + ikx]Vu,D(k) dk.

(8)

The displacement in the half-space B can be represented in the same way as equation (6), which is also a complete rep-

381

An Efficient Approach to the Seismogram Synthesis for a Basin Structure Using Propagation lnvariants

resentation. The boundary condition of continuity of displacement on the interface can then be written as a continuous function of x: (.J~UAYUA) (X) -Iv (,.~DAVDA) (X) = (.~UBVUB) (X) + (:;~40BVD~) (X).

(9)

Taking the normal [nx, 0, n z] to the interface to lie in the positive z-direction, i.e.,

that must lie below the deepest part of St and above the source depth. Note that this is where the so-called "Rayleigh ansatz" (see equation 13.100 and Figure 13.13 in Aki and Richards, 1980) is invoked, which ignores upgoing waves scattered from low points on the interface. Applying (O/Oz) v(x, z = 0) to equation (6), the tractionfree condition for the horizontal free surface yields the following relation between Vva(k) and VDA(k): exp(--iv¢ad)VvA(k ) = exp( + iv¢ad)VoA(k ).

nx(X) = -h'(x)/~/1 + [h'(x)] 2,

(10)

(15)

Then,

nz(X ) = 1/~/1 + [h'(x)] 2,

where h'(x) = dh(x)/dx, the boundary condition of continuity of traction on the interface can be written as OV[x,z = d + h(x) - 0] #a nx 0~

ov

+ nz ~ Ix, z = d + h(x) -

VDA(k ) =

exp(-- i2vl~Ad)VvA(k)

= ~ exp(-i2V~Ad)~(k'

- k)vua(k')dk'

(rRflU VUA) (k).

)

O] ,/1 + [h'(x)f

(11)

= /z~ n x-~x [x,z = d + h(x) + 0] + nz-~z Ix, z = d + h(x) + 0]

1 + [h'(x)]2,

where/.tA and/t 8 are the rigidity of layer A and half-space B, respectively, and we have multiplied the traction by ]1 + [h'(x)] 2 because its inverse is common to the traction forms in both layer A and half-space B (see equation 10). Substituting equation (6), its counterpart for half-space B, and equation (10) into equation (11), we get (NuAVuA) (X) "}- (,~DAYDA) (X) = (~uBVuB) (X) + (NDsVo~) (X), (12)

(16)

Using equations (14) and (16), the integral equations (9) and (12) can be considered as simultaneous equations having two unknown functions VuA(k) and VoB(k). These integral equations correspond to two coupled integral equations (equation 10) in Aki and Lamer (1970). They reduced the integral equations to the infinite-sum equations by assuming a periodicity in the shape of the interface and approximated them by the finite-sum equations. Instead of solving these equations directly in x, they took the Fourier transform of these equations by using FFT. Then they solved the resultant simultaneous linear equations. In this article, we propose a more efficient approach than that of Aki and Lamer (1970) to solve the coupled integral equations for VvA(k) and Vow(k), i.e., equations (9) and (12). Our approach reduces them to a single integral equation only for VVA by employing the propagation invariants. In the next section, we derive the propagation invariants for our problem by following Takenaka et al. (1993).

where the integral operators Nv, o have been defined as Propagation Invariant (Nv,gVv,D) (x) =-- ~

[-- h'(x)" k +- v~] • exp( +-ivph(x) + ikx)vv, o(k) dk.

(13)

The Rayleigh hypothesis is that the reflected and transmitted fields due to the irregular surface consist only of waves moving away from the surface. Applying the Rayleigh hypothesis for our problem, the upgoing wave in the basement VVB is taken as the upgoing part of the wave field radiated by the source, which in the wavenumber domain can be represented as

Kennett et al. (1990) established a spatial propagation invariant f [ulT(x, z)t2(x,z) -- tT(x,z)U2(X,Z)] dS

for two elastic displacement fields ul and u2 and their associated traction fields tl and t2. For any two surfaces $1 and $2 spanning the whole x-y plane, s~ [ur(x'z)t2(x'z)

vvB(k) = exp[--iv~B(H -- d)]~U(k),

(17)

- t~(x'z)u2(x'z)] dS

(14)

where 9IN(k) is the field produced by the source at a level H

: fs2 [ur(x'z)t2(x'z) - t~(x,z)u2(x,z)] dS

(18)

382

H. Takenaka, M. Ohori, K. Koketsu, and B•L.N. Kennett

[see Kennett (1984) or Kennett et al. (1990) for details]. In two-dimensional SH cases, the propagation invariant (18) is reduced to

i

MaOBNVB -- NAOBS~fu8= ~ 22B,

(26)

where 22 is an operator that acts on a wave, as follows:

s, [v~(x,z)t~2(x,z) - t:(x,z)v2(x,z)] dS (22v) (k) =-- f 2ttvl~6(k' - k)v(k') dk'

= fs2 [Vl(X'Z)ty2(X'Z) - tyl(X'Z)V2(x'z)] dS,

(19)

= 21~v~. v(k). where V1,2(X,Z) are y-components of the displacement fields and tym(X,Z) are y-components of their associated traction fields. Here we consider a homogeneous full-space having the same medium parameters as the half-space B of our original problem and set the horizontal surface z = d and the irregular surface z = d + h(x) in this full-space. When the invariant (19) is adopted for these two surfaces,

f

If the scaling parameter were taken e/~ = (2/tvp)- 1/2 instead of eH = 1, the operator 22 would reduce to the identity operator f,,as shown in Takenaka et al. (1993)• In the next section, by embedding the propagation invariants derived above, i.e., equations (24) and (26), in the coupled integral equations (9) and (12), we obtain a single integral equation• A Single Integral Equation

[vl(x,a)ty2(x,d) - trl(x,d)v2(x,d)] dS (20)

= f [va(x,d + h(x))ty2(x,d + h(x)) - tyl(x,d + h(x))v2(x,d + h(x))] ~/1 + [h'(x)]2 dx.

If we take Vl = v2 = Mo~VDB in (20), some manipulations (see Takenaka et al. (1993) for details) reduce (20) to

0 = f VDB(--k) [(MADBNDB-- NABMDB)VDB] (k) dk,

(21)

where the A transform for an integral operator

(234.) (x) = f M(x,k) . dk

(27)

Now we return to the integral equations (9) and (12). In this section, we construct a single integral equation only for VVA from these integral equations. First, we eliminate voB(k) by using the propagation invariants (24) and (26). The action of NaoB on (9) and of 5~IAB on (12) followed by subtraction and use of the propagation invariants gives i

(KDB.vAVvA) (k) + (KDB,DAVDA) (k) = ~ 2]IBVflBVuB(k),

(28)

where (22)

~DB,UA -- ~/~ADB~:UA

:A

(29)

KDB,DA ~- OV[~)B~[DA -- O~ADBMDA.

is defined as

(LMA. ) (k) = f M(x, - k). dx.

(23)

Equation (21) can be cast as the following relation between the integral operators:

By some manipulations, the kernels of these integral operators are found as follows: k_______~' -

KDB,uA(k,k') = ~

(/'taV~B+/dAPtflA) + Q2Bk -}- 11Ak') FflB -- "V~A

•J_~ exp[ - i(v/~ - V'~A)h(x)]exp( -- i(k - k')x) dx, M ~ B N ~ -- N ~ 8 ~ B

= O.

(24) (30)

If we take vl ulations give

~/[DBVDB and v2

=

=

ffV(UBUUB, some manip-

KDB,DA(k,k') = ~ i

fVDB(--k)'~-~2flBVflBVuB(k)dk

i [

• f+~

= fVD13(--k)"

, k-k'l (ltBVpB-- ~lAFtfla) -}- (tzak + ,UAk ) ~ J

e x p [ -- i(v/~B +

V'~A)h(x)] e x p (

-- i ( k -

k')x) dx. (31)

[(SMABNVB - :vAaJq4v~)Vw] (k) dk.

(25)

This yields another relation between the integral operators:

Substituting (14) and (16) into equation (28) gives the following single integral equation for VUA:

An Efficient Approach to the Seismogram Synthesis for a Basin Structure Using Propagation Invariants N

[( ff~DB,UA -1- ff~DB,DA~RflU)llUA] ( k ) i = 2---~2/tBv~ exp( - iv~n(H - d))dN(k).

(32)

P(k,k') = 2TrY(k-k') + ~] .xl - xt-I [exp(ipt) - exp(ipl_ 1)] / = 1 lpl -- ipl--I 1 , - i(£ - k ) {exp[ - i(k - k')b] - exp[ - i(k - k')a] } (b - a) (39)

-

We set

A(k,k')

=-

KDB,uA(k,k' ) q- KDs,DA(k,k')'exp(--i2v'~ad )

383

(for k ~ k') (for k = k')'

(33)

and

where i

b(k) -- ~ 2/.tsv~ exp( - ivps(H - d))f/N(k),

(34)

where we have used (16). Then, equation (32) can be written as

f ) ? A(k,k')VuA(k' dk' = b(k).

(35)

Note that the function b(k) is a known source quantity and that the unknown function in this equation is only VuA(k). This single integral equation is one of the most important results of this article. If the integral equation (35) can be solved, the displacement at the free surface Vo(X) can be obtained from (6) and (16) as follows:

Pl

-'= V h l

(1 = 0 , ' " ", N), a = x 0 < x I < . . . < x~v-1