Integral Representations in Wavefield Inversion

Chapter 15 INTEGRAL REPRESENTATIONS IN WAVEFIELD INVERSION In this chapter I will outline the general principles of seismic inversion and imaging based on integral representations of scalar and vector wavefields developed in Chapter 14. The integral equation representations and the corresponding linear and nonlinear approximations introduced in the previous chapter constitute a powerful tool for seismic inversion and for imaging seismic cross-sections. We will see that this method makes it easier to determine the sensitivity matrix, which lies at the cornerstone of any inversion algorithm. Also, the integral equation approach provides an analytical insight into the migration transformation, which is now one of the most widely used methods of seismic imaging. As in the previous chapter our discussion will be restricted mostly to the fields satisfying wave equations. However, in the last section I will outline the basic principles of nonlinear elastic inversion. 15.1

Linear inversion methods

Consider a 3-D seismic model with a background (normal) slowness distribution ^^(r) and a local inhomogeneity D with an arbitrarily varying square of slowness 5^(r) = 5^(r)-|-A5^(r). We will examine, in parallel, two cases: the propagation of the acoustic field and of the vector wavefield in this model. The acoustic field in the frequency domain, p(r,a;), satisfies the equation VMv,oj)+uh'iv)p{v,uj)

= -r{T,oj),

(15.1)

where /^(r,(x;) is the strength of an external source of energy. The vector wavefield in the model is excited by an arbitrary source with the distribution of volume force characterized by cf)^ (r, uj). In frequency domain this field satisfies the vector Helmholtz equation: V^u{r,Lu)+uj^s'^{r)u{r,uj)

= - 0 ^ (r,a;).

(15.2)

Both the acoustic and the vector wavefields satisfy the corresponding radiation conditions at infinity (see Chapter 13).

468

Integral representations in wavefield inversion

The scattered acoustic field and the vector wavefield in this model can be represented, according to formulae (14.19) and (14.80), by the following integrals p'{rj,u;)=iu^

III

G'^{vj\v',uj)/^s\v)p{v,uj)dv

= u^G^ {As^

,

(15.3)

and u^(r^.,c.;) = Lu^ fff

G^{vj\r;ij)

• A^^ (r) u{T,uj)dv = uj^G^ [As^u] ,

(15.4)

where G'^{rj\ T;UJ) and G^(rj| r;a;) stand for the Green's function and tensor defined for an unbounded medium with the background slowness s^. We assume that the scattered fields, p^{rj,uj) or u*(r.,a;), are given on some surface of observation S. The goal is to find the anomalous square slowness As^ ( r ) . 15.1.1 Born inversion of acoustic and vector wavefields The most common approach to the solution of the wavefield inverse problem is based on linearization of the integral representations (15.3) and (15.4) using Born approximations: P'{TJ,UJ)^UJ^G^

[A^V]

=^' ///

G^{rj\T;u;)As^T)p'{T,uj)dv

(15.5)

and u'{r.,uj)^uj^G^[As^u']

= uj^ fff

G^{rj\r',uj)

• As^ {r)u'{r,uj)dv.

(15.6)

Comparing formulae (15.5), (15.6) and (14.29), (14.84), we see that Born approximations can be expressed as the Frechet derivative operators (or Frechet differentials) calculated for the background square slowness 5^ and the anomalous square slowness As'^: p'(Yj,uj) ^ Lu^G^ [AsY] = Fa{slAs^), (15.7) u%r^,uj) ^ uj^G^ [As^u'] = F,{slAs^).

(15.8)

The argument in the expressions for the Frechet diff'erentials Fa(5^,As^) and F^(5^,A5^) consists of two parts. The first part, 5^, is the background square slowness distribution, from which we calculate the forward modeling operator variation; the second part, As^, is the anomalous square slowness. We will use below the following simplified notations for the Frechet differentials:

and

Linear inversion methods

469

where the upper subscript "6" denotes that the differential is calculated for the background slowness. Substituting the Born approximations (15.7) and (15.8) into formulae (15.3) and (15.4), we arrive at the linearized form of the forward modeling operators p^ = F^{As''),

u^ = Ft{As').

(15.9)

In the inverse problem solution we assume that the scattered fields, p^ or u*, and the background slowness st are given. The goal is to find the anomalous square slowness distribution, As'^. As a result, we have hnear equations (15.9) with respect to A 5 ^ Following the general principles of regularization theory (Chapter 2), we solve the linear inverse problems (15.9) by regularization methods, imposing additional conditions on the class of inverse models with corresponding stabilizing functionals. We introduce a Hilbert space Da of the acoustic field data, given on the surface of observations S with the metric {Pi Q)Da^^^

/ / P(r,cc;)g'*(r,a;)d5da;; p, q e Da, oo e ft,

(15.10)

and a Hilbert space Dy of the vector wavefield data, given on the same surface with the metric (u, v)^^ = Re / / / u(r,c<;) • V* {r,uj)dsduj] u , v G Dy, cj G f2,

(15.11)

where ds denotes a differential element of the surface S, and the integration is conducted over the variable r e S. We also consider a Hilbert space M of models m (square slowness m = 5^) with the metric {m^'\ m(2))^ = / / / m(^) (r) m^^) (r) dv; m^'\m^^^ G M.

(15.12)

According to the conventional Tikhonov regularization method, we replace the solution of the linear inverse problems (15.9) by a minimization of the corresponding parametric functionals, for example, with a minimum norm stabilizer: P: = (F^iAs^) - da, F^iAs') + a {Wa {As'-Asl^),

- da)^^

Wa (As'-Asl^))^

= min,

(15.13)

or P: =

{Fl{As')-d,,Fl{As')-d^)^

+a {W. {As'-Asl^), W, {As'-Asl^))^ = min,

(15.14)

470

Integral representations in wavefield inversion

where da = p^ and d^ = u* are observed acoustic and vector wavefield data, As^p^ is some a priori model of the anomalous square slowness distribution, and Wa, Wy are the weighting operators. The minimization problems (15.13) and (15.14) can be solved, for example, by the conjugate gradient or re-weighted conjugate gradient methods introduced in Chapter 5. 15.1.2 Wavefield imaging by the Born approximations We develop now a fast wavefield imaging technique based on the Born approximations. This technique employs formula (5.91) for an approximate regularized solution of the linear inverse problem. For acoustic or vector wavefield inverse problems (15.9), this formula takes the forms As^^k{W:WaVF',*{f), (15.15) or As^^k

{W:Wy)-^ F ^ (u^),

(15.16)

where F^* and F^* are the adjoint Frechet derivative operators, and the scalar coefficient A: is a solution of the minimization problems: \\kFl {{W:Wa)-'

Ft if))

- pf

= min,

(15.17)

||/cFt {{W^Wy)-' F ^ (u^)) - u^ll' = min.

(15.18)

or For a practical realization of the imaging algorithms (15.15) or (15.16) we have to find explicit expressions for the adjoint Frechet derivative operators F^* and F^*, according to these formulae: {FUAs'),q)^^

=

{As',Ft{q))^,aud

(Ft {As') , v)^^ = (A5^ Ft ( v ) ) ^ .

(15.19)

Using the definitions (15.10), (15.11) and (15.12) of inner products and the expressions (15.7), (15.8) for the forward operators, we can rewrite the first formula in (15.19) as {F^ {As') ,q)^^=ReJ

u;^Re

u? j j G^ [A^V] q*dsdw

HI III ^"^^"'"''^^^'' ^""^^^'''^^'^'"^* ^'^ dw

= / / / As' (r) Re [ / oj'p"{r, w) ff JJJD

[JU

G r ( r | r; w)? (r, uj) dSdu; dv

JJS

= (A.^F^(^))^.

(15.20)

Linear inversion methods

471

From the last formula we have (AS^

IF^

\

L

(q) -Re f ^Y'{T,U;) Jn

0.

ff Gr{r\ v-u:)q {T,UJ) dSdw J Js

(15.21)

Equation (15.21) should hold for any As^, for example for As^ = F^ (q) -Re f ct;V*(r, cu) [f G^*(?| r; uj)q (?, ou) dsdou. Jn J Js Therefore, 0.

Ft {q) -Re f u;V*(r, a;) / / Gr*(?l r; ^)9 (?, a;) dScLv M

From the last formula we conclude that the adjoint Frechet derivative operator F^* is given by the formula F^' {q) = Re [ UJ^P'\T, UJ) [[ G^ (?| r; cu) q (?, UJ) dSduj. Jn J Js

(15.22)

In a similar way we can obtain the formula for the adjoint Frechet derivative operator F^'*' :

F,^* (v) = Re / U'^U'\Y,UJ) • / / G^ (?! r;a;) • v (?,a;) d ? ^ . Jn J Js

(15.23)

Note that we can drop the symbol Re in the above formulae if we integrate over a symmetrical frequency interval ft (for example, from — oo to -foo) because the imaginary part of the integrand is anti-symmetric in UJ. On this understanding, and substituting expressions (15.22) and (15.23) into formulae (15.15) and (15.16), we finally arrive at the following imaging formulae: As^^k{W:Way' f u^v'\v,u) ffGr{v\Y',uj)f{v,uj)dSdw, JVL J Js

(15.24)

or As^^kiW^Wy)-^

/ u ; V * ( r , a ; ) - [[ G^ {r\r;uj) - u'{7,LU) dSduj. J JS

(15.25)

JQ

The weighting operators Wa and Wy for acoustic and vector wavefields are Unear operators of the multiplication by the functions Wa (r) and Wy (r) equal to the square root of the integrated sensitivity (see formula (3.78) in Chapter 3): Wa

sSa, and Wy = y/Sy.

(15.26)

472

Integral representations in wavefield inversion

The integrated sensitivity, Sa^ in accordance with the definition (3.75) and formula (14.33), is calculated by the formula:

\yG^{r\T;u;)p\r,u;)\\^^. Formula (15.27) acoustic field to vector r. Similarly, field to the local

(15.27)

should be treated as the integrated sensitivity of the scattered the local slowness anomaly located at a point with the position we can find the integrated sensitivity of the scattered vector waveslowness anomaly, using formula (14.88): Oy

(15.28)

u;^G^(?|r;a;)-u^(r,a;) Dr,

The weighting functions, according to (15.26), are equal to Wa (r) =

Wy{r)

^\\UJ^G-^{T\T;U)P'{T,UJ)\\J,^

c.;2G^(r|r;a;)-u^(r,a;)

(15.29) Dy

Substituting formulae (15.29) into (15.24) and (15.25), we finally find 2 ,.^ ^^ J ^ a ; V * ( r , a ; ) / / ^ G r {r\r;uj)f {v.u)dSdu As^ (r) ^k ||a;2G-(r|r;a;)pHr,^)||p.

(15.30)

and 2 ,.^ ^ ^ J ^ c ^ V * ( r , a ; ) - / / 5 G r ( ? | r ; c c ; ) - u M ? , a ; ) r f ? c L ; As' (r) ^k cj2G^(r|r;a;) •u^(r,a;)

(15.31)

Dv

We can give a clear physical interpretation of formulae (15.30) and (15.31). According to the reciprocity theorem 27 (formula (14.24)), the surface integral term in equation (15.30) can be treated as a complex conjugate acoustic field p ^ . p ^ (r,a;) = j j G^ (?| r; a;) p - (?, a;) ds,

(15.32)

generated at the point r by the fictitious acoustic source, distributed over the surface of observations S with the surface density / | , equal to the complex conjugate scattered acoustic field: rs{v,uj)=p'*{t,uj). (15.33)

Linear inversion methods

473

Note that this interpretation is similar to one given in Chapter 10 for electromagnetic Born imaging (formulae (10.33) and (10.34) ). The approximate anomalous square slowness in formula (15.30) is calculated as an integral over the frequency range Q of the product of this auxiliary field p^* and the complex conjugate incident field at the point r, normalized by the norm of the product of the incident field and the corresponding Green's function at the same point r:

As' (r) ~ kAs' (v) As {r)^kAs,(v),

and As' (v) In'^'p'*i^''^)p''* i^'^)^ a n d A 5 i ( r ) - ||^2G»(y|r.^)pi(r,a;)||^/

.15 34) ^^^•^^>

where the scalar coefficient k is the solution of the minimization problem

Formula (15.31) has a similar interpretation, except in this case we should introduce, by reciprocity, a vector wavefield u^, u ^ (r, uj) = 11 G^ (?| r; LJ) • u ^ (?, u) ds, generated at the point r by the fictitious vector wavefield source, distributed over the surface of observations S, with a surface density H?,u;)-u^*(r,u;).

(15.35)

As a result, we arrive at the following formula for anomalous square slowness As^ (r) ^ kAsl ( r ) , and Asl (r) =f'—::—-^—^ -^^—, |a;2G-(r|r;c^)-uHr,a;)

(15.36)

ID.

where the scalar coefficient k is the solutionQ ofof the the minimization mini: problem | / c F ^ A 5 2 ) - u ^ | | ' = min. Note that in the definition of the auxiliary fields p^* (r,a;) and u^* (r, uj) we use the complex conjugate Green's functions. It can be demonstrated that the transition to the complex conjugate function in frequency domain corresponds to changing a sign in front of time variable t in time domain. Thus, using the convolution theorem (Arfken and Weber, 1995, p. 863-865) and taking into account reciprocity theorem 27, we can write the inverse Fourier transform P^ (r,t') of the auxiliary scalar wavefield in the time domain as follows: 1

+00 r+oo

iujt'

du

474

Integral representations in wavefield inversion r+oo

27rJ_^ JJss

-ILL

or

(?| r; iu) p' (?, Lu) dse-'^^^'dcj

+00 _

G^ (r', t'|r, t) p' (r, t) dtds,

(15.37)

5."

where we use the notation

G^iv',t'\ r,t) = G^ir',-t'\

r,-t),

and the wave above the acoustic field P^ reflects the fact that we change a sign in front of time variable t'. The minus sign in front of time variable t in the expression for the Green's function indicates that the wave propagation occurs backward in time. Therefore, the field P^ (r,t') is a so-called "back-propagated wavefield." We will discuss the physical properties of this field in more details below in the section on seismic migration. Taking complex conjugate of both sides of equations (15.34) and (15.36), and noting that anomalous square slowness As^ (r) is a real function, we present these formulae in equivalent form: o/ X

L(jo'^pHr,(jj)p^ (r,(ju)duj

^^

||a;2G-(r|r;a;)pHr,a;)|lz)/

'

, ^

2/ ^ _ /^a;V(r,u;)-u^(r,u;)dc^ As^ (r) ^ k-p '-^-^ ' ' \. . ||a;2G^(?|r;c^).u^(r,u;)||^

^

(15.39)

According to the correlation theorem (Press et al., 1987, p. 383), the inverse Fourier transform of the product of spectrum of one function and the complex conjugate spectrum of another function is equal to correlation of these functions. Therefore, we can write the numerator in formula (15.38) as a cross correlation of the time derivatives of the back-propagated scattered field and the incident wavefield: 4-00

/

r+oo

uP'p'{Y,uj)p^{Y,uj)duj -00

[-iujp\Y,ijo)\

[-za;p^* (r,a;)] * do;

J —00

/.+00 , i

=

= I

/ J-00

P {r,t)P

1^

I

{r,t-t')dt\

r+oo

= l,^^

^i

•^R

P (r,t)P

{r,t)dt,

J-00

where the dot over the functions P* and P^ denotes their diff"erentiation with respect to time. Formula (15.39) for the vector wavefield permits a similar physical interpretation. We will show below in this chapter that the back-propagated field corresponds to a migration transformation of seismic data. Thus we see that Born imaging can be treated as an algorithm from a family of migration transformations.

Linear inversion methods

475

15.1.3 Iterative Born inversions of the wavefield The ideas of iterative Born inversion, introduced in Chapter 10 for an electromagnetic field, can be extended to wavefield inversion as well. Following the basic principles of this method, we first write the original integral equations for the acoustic or vector wavefields (14.20) and (14.81) as the domain equations for the wavefield inside the anomalous domain D: P{T',UJ)

= uj^G^ [As^ (r)p(r,a;)] +p'(r',a;), r' G D,

u(r ,0;)

= U ; 2 Q ^ [AS^

(r) u(r,

CJ)]

+ u ' ( r ' , a ; ) , r' G D,

(15.40) (15.41)

where Gyj and G^^ are the scalar and vector Green's operators determined by formulae (14.18) and (14.74). We can now use similar equations to connect the scattered acoustic [da = p^ = p — p'^ ) or vector (dv = u® = u — u*) wavefield observed on the surface S with the corresponding wavefield within the anomalous domain da{Tj,Lu) = LO^G^ [As^ (r)p(r,u;)] , r^- G 5, dv(rj,a;) =a;2G^ [As^ (r) u(r,cj)] , r^- G S.

(15.42) (15.43)

Note that each equation, (15.40), (15.41) and (15.42), (15.43), contains the product of unknown functions, As^ and p, or As^ and u. Therefore, these equations are bi-linear with respect to the corresponding unknowns. However, if we specify one of the unknowns, the equations become linear. For example, we can subsequently find As^ from equations (15.42), (15.43) for specified p or u, and then update p or u from equations (15.40) or (15.41) for predetermined As^, etc. Within the framework of this method the Green's functions G'^ and G^^ and the incident wavefields p^ or u* stay unchanged. As in the electromagnetic case, there is also another technique, the distortedBorn iterative method, which is based on updating the incident field and Green's functions after each iteration, according to the updated parameter As^ (Chew, 1990). Note in conclusion of this section that iterative Born inversion requires the application of the regularizing methods to make the solution stable. 15.1.4 Bleistein inversion In the case of a high frequency scalar wavefield it is possible to find an explicit inverse formula for Born approximation. The method of solving this problem has been systematically developed and studied in the works of Bleistein, Cohen and Stockwell (Bleistein, 1984; Bleistein et al., 2001). It is based on a high frequency WKBJ approximation for the Green's functions (13.80) and (13.83). I will outline here only the basic ideas of this technique, referring the reader to the books cited above. The Bleistein method can be more easily formulated for a 1-D wavefield inverse problem, which will be the starting point of our discussion.

476

Integral representations in wavefield inversion

One dimensional inversion Let us assume again that the wavespeed, c, the pressure field, p, and the source field, /^, are functions of the vertical coordinate z only. We also assume that the wavespeed in a 1-D model can be represented as follows:

where Cfe(z) is the background (normal) 1-D velocity distribution, and A5^(2:) is some anomalous square slowness which is non-vanishing only within some interval of positive z : A52(z) = 0 , if z < 0 . * The 1-D field p may describe, for example, plane acoustic wave propagation in the vertical direction in the medium. This field, according to (13.57), satisfies the equation — p ( z , z', a;) -h ^273^(^' ^'' ^ ) = -^{z-

z'),

(15.45)

where we use the impulse source function /^ [z) = 6 {z — z') located at the point z = z'. As one can see, the function p{z^ z' ^uo) depends on three variables: frequency a;, source position z'^ and receiver position z. We assume also that p(z,z'^uj) is bounded for all z and satisfies the 1-D radiation condition -—^(z, z', uj) =f —j-^p{z, z', cj) —^ 0, as z —> ±oo. (15.46) dz c[z) The total field in this model, as in the general 3-D case, can be represented as the sum of the incident and the scattered fields: p{z, z\ UJ) = p'{z, z', uj) + p'{z, z', u). The incident field p*(z, z', uS) satisfies the 1-D Helmholtz equation with the background velocity distribution -J{z,z',uj) + -^j-^p\z,z\uj) = -6{z- z'),

(15.47)

and the radiation condition -—p'(z, z', uo) =f —7-TP'('2^5 ^^ (x;) —> 0, as z —> ±00.

dz

Cb (z)

The scattered field p^{z, z' ^uj) also satisfies the radiation condition and the same 1-D Helmholtz equation, but with a different right-hand part, -^f{z,z',uj)^^--f{z,z\u)

= -a;2As2(z)p(z,z',a;),

(15.48)

*In the section on Bleistein inversion and later on we will direct the vertical axis z downward.

Linear inversion methods

477

where As^ (z) is the anomalous square of slowness determined by an equation similar to (14.8):

Differential equation (15.48) can be transformed into an integral equation for the scattered field, using the 1-D Green's function g'^{z\ z'\uj) for the 1-D Helmholtz equation, which is dependent on the position of the source z', the observation point z, and the frequency u. This function satisfies the 1-D Helmholtz equation

'^l^S^^

+ - i ^ , n - l -';-) = -S{z - z%

(15.49)

and the corresponding radiation condition at infinity. Evidently in our case the Green's function is equal to the incident wave: V\z,u)=p\z,z',u)=gl^{z\z'-u)

.

(15.50)

Note that Green's theorem in one dimension can be written as / '°{u{z)Hv{z) J-oo

- v{z)Hu{z)}dz

=

u{z)^^ dz

*°°

du{z)

-oo

"^

l+oo

(15.51)

;re H is the 1-D Helmholtz operator: „

(f

w'

We choose functions u{z) and v{z) to be equal to the scattered field p^{z, z\ uj) and the Green's function g'^{^z\ z'\u). Using Green's theorem, equations (15.48) and (15.49), radiation condition (15.46), and the fact that lS.s^{z) = 0 for 2 < 0 we obtain an integral equation for the scattered field p^{z^ z\uj): p+00

f{z,,z',uj)=u^

/

^s\z)[f{z,z\uo)^v\z,z\u)]g]^{z\zy,uj)dz,

(15.53)

We can also use the 1-D Green's theorem, equation (15.51), to prove the reciprocity theorem in one dimension, similar to the general reciprocity theorem 27: gl^{z'\z";u)

= gt{z"\z';iv).

(15.54)

Therefore, equation (15.53) can be written in equivalent form as r+oo

p\zj,z',Lj)=u?

/ Jo

i\s''{z)\p\z,z',uj)+p\z,z',uj)]gl^{zj\z-u)dz.

(15.55)

478

Integral representations in wavefield inversion

In the case of the Born approximation we neglect terms containing the scattered field in (15.53) and leave only the incident field: /•4-00

f{zj,z',uj)=u^

/ /^s^z)v\z,z\u)g^{z\zj',uj)dz. (15.56) Jo Taking into account that, according to equation (15.50), the incident field coincides in our case with the Green's function, we can write: p%zj,z\uj)=uj^

/ As^z)g^{z\z';uj)g^iz\zj;Lu)dz. (15.57) Jo The original Bleistein method is based on direct inversion of integral equation (15.57) with respect to anomalous square slowness As'^ (z). In a general case, this equation is ill-posed, as is any geophysical inverse problem, and its solution requires application of the corresponding regularization methods outlined in the first part of this book. However, in the Bleistein method the Green's functions are represented by their high frequency (WKBJ) approximations. This reduces equation (15.57) to a specific class of integral equations which do have stable inverse solutions. Actually, this class is formed by Fourier integral operators, or by so-called pseudo-differential operators (Bleistein et al., 2001). We will see in the next section how this reduction can be done, and what analytic methods can be used for inversion of the corresponding integral equations. Note that numerical implementation of these theoretically stable inverse solutions will still require an application of the corresponding regularization methods. 1-D inversion in a constant-background medium We will discuss in this and the following sections the solution of the inverse problem for two cases: a) a constant background wavespeed, and b) a variable background wavespeed. We consider a simple "zero-offset problem" with source and receiver at the same location at the origin: z' = Zj = 0. In this case equation (15.57) can be simplified as r+oo

p'{0,u;)=p'{0,0,uj)=u;^

/

As^ {z)[g^{z\0]u)]^

dz.

(15.58)

Jo

According to equation (13.79), the Green's function in a model with a homogeneous background, c^ (z) = CQ, takes the form iu;\z-z'\/co

iuj\z\/co

^"(z| / ; a;) = -c^—-.

= -CQ-ZICJ

, if z' = 0.

(15.59)

ZILU

Substituting the last expression in equation (15.58), we obtain p'{0,uj)=p'{0,0,u;)

= - ^

/

4 Jo

As^z)e'^'^'^''dz,

(15.60)

Linear inversion methods

479

where \z\ was replaced by z, because the variable of integration is positive. Introducing an auxiliary "wavenumber" 2a; /c = — , Co

and taking into account that As^ (z) = 0 for z < 0, we can view the last formula as a Fourier transform from the vertical coordinate z domain to the wavenumber k domain: p'{0,Lo)=p'{0,cok/2)

= —j

I

As''{z)e'^'dz.

(15.61)

Applying the inverse Fourier transform to the last formula, we arrive at the following expression 2

/'+00

1

^As^ {z) = -— / p^(0, cok/2)e-"''dk. (15.62) 4 27r J_^ From the last formula we find at once the explicit integral expression for anomalous square slowness: As^ (z) =

. / TTCQ

f{0,ij)e-^'^'^''duj.

(15.63)

J-OO

Reflectivity function and bandlimiting of data The Bleistein method helps to determine not only the anomalous slowness distribution, but also another important characteristic of a seismic model, a reflectivity function. In a 1-D case this function is a spike train with spikes located at the positions of reflectors in the seismic model, whose heights are proportional to the reflection coefficients of the corresponding reflectors. It was demonstrated by Bleistein et al. (2001) that the reflectivity function can be represented as a combination of bandlimited delta-functions. The observed seismic data are usually bandlimited, because the frequency range is limited by the natural time constant of the source process, by receiver and survey design, by seismic preprocessing required for noise removal, etc. (see Aki and Richards, 2002, and Bleistein et al., 2001 for further discussions). The bandlimiting is an important characteristic of observed seismic data which should be taken into account, especially if we use the high frequency asymptotics for data inversion. Mathematically the bandlimited data can be introduced by filtering the theoretical signals. We consider, for example, the effect of this filtering on a delta-function. There is a well-known integral representation for a delta-function: 6{z-z')

1 /'+°° = — / exp [-iio {z - z')] duj. 27r 7_oo

(15.64)

480

SB {Z

Integral representations in wavefield inversion

-

Following Bleistein et al. (2001), we define the bandlimited Z') , as SQ{Z-Z')

1

= —

/*+°°

F (CJ) exp

[-iuo {z - z')] dcu,

delta-function,

(15.65)

where F (LJ) is some symmetric and non-negative function, characterizing the bandlimited filter in frequency domain. Some examples of bandlimited delta-functions can be found in mathematical textbooks on the Fourier transform. The main difference between the original delta-function, 6 {z — z'), and its bandhmited analogue, bs {z — z'), is that instead of the infinitely narrow and infinitely high spike oi6{z — z') at the point z', the bandhmited delta-function has a spike of limited width and height, located at the same point. Using a simple 1-D model, we will illustrate below the basic principles of the reflectivity function definition based on bandlimited data. Let us assume first that the anomalous square slowness distribution As^ {z) is described by a piecewise constant function as shown in Figure 15-1. This means that the wavespeed is constant within each layer, equal, say, to CQ within the first layer, ci within the second layer, etc. We can describe, for example, the first step on this curve by the Heaviside step function As^{z) = bH{z-h),

(15.66)

where h is the depth of the first reflecting boundary and b is the size of the step in square slowness, which is determined according to equation (15.44) by the wavespeed change on the first boundary: 6 = 1 - 1

(15.67)

Note that in this case we can express the wave speed c^ as co/^l

+ bcl

(15.68)

The Heaviside step function in (15.66) is equal to zero for z < /i, and to 1 for z > h.lt \s well known that it can be expressed by the following Fourier integral: 1

+00 r /•+°°

27ry_oo

L +00

1

exp [—iuj {z — h)] duo

^^

1 1 r"^ 1 = T; — TT '^exp[—iuj{z — h)]duj. (15.69) 2 27r y_oo luj It is easy to find the exact solution of the wavefield equation within the first layer. As usual, it can be written as a superposition of the incident and scattered fields,

Linear inversion

methods

481

slowness

reflectivity

z=h

Figure 15-1 A model with piecewise constant distribution of slowness (left panel). The corresponding bandlimited reflectivity function (right panel) is represented by a combination of the bandlimited delta-functions with the positions of extremum at the reflecting boundaries, and the sizes of the peaks proportional to the corresponding reflection coefficients

p (z, uj) = p' (z, uj) H- p' (z, cj), z
and the scattered field can be calculated as p^ (z, (jj) = —A—- exp {—ioj {z — 2h) /CQ) .

(15.70)

The reflection coefficient A in the last formula is calculated based on the usual definition: Ci - C o A = (15.71) Cl + C o

We assume now that the scattered field is observed at the surface of the model at 2; = 0 : Co (15.72) p^ (0,0;) = -A^r— exp {2iu;h/co). Substituting formula (15.72) into inverse formula (15.63) of the Bleistein method, we obtain

482

Integral r e p r e s e n t a t i o n s in wavefield i n v e r s i o n

4

f^^^ Cn

A3^ (z) = —^A / TTCQ

- ^ exp [-2iuj {z - h) /CQ] du,

J-OO

(15.73)

220;

where we denote by AS^ {z) the anomalous square slowness recovered by inversion. Using the Fourier representation of the Heaviside function (15.69), we write: A?" {z) = — ^ A /

^

exp [~2iuj {z - h) /CQ] d (2a;/co)

/f(z-/i)-i

4A

(15.74)

Comparing formula (15.74) for the predicted anomalous square slowness AS^ {z) with the expression (15.66) for the original As^ {z), we find that the inverse result has the step at the right depth, but the magnitude of the step is not exact and the entire curve AS^ {z) is shifted with respect of the true curve As^ {z) by a constant (2A/co). The shift of the inversion curve actually depends on the value of the observed scattered field at zero frequency, which can never be practically determined in a real-world experiment (see Bleistein et al., 2001, for more details). The discrepancy in the step sizes can be easily resolved if we evaluate the reflection coefficient value for a small step 6, taking into account formula (15.68) ^^ ^ 1 - v T T 6 ^ ^ 1 - (1 + \bcl)+0{bcl) l + yrT6^ l + {l + | 6 c 2 ) + 0 ( 6 c 2 ) IhfP.

= -^ S °

1

+ O {bcl) = —bcl + O (6c2) , for small bcl

(15.75)

Comparing the leading order term in (15.75) with the step value in (15.74), we see that the inverse formula provides the correct estimate of a step, 6, with the accuracy determined by the leading term in the representation (15.75). However, Bleistein inversion is based on a Born approximation formula which is valid only to leadingorder in hc^. Thus the theoretical inversion result for this simple model fits well the basic assumption of the Bleistein method. The main significance of this example for the solution of our goal, outlined in the beginning of this section, is that it illustrates how one can find the reflection coefl^icient itself from the observed scattered field data. Indeed, let us consider the bandlimited data, v%{^^^)^ by including a bandlimited filter F {u) in expression (15.70): p% [z, uj) = -F [uj) A ^ exp {-iu {z - 2h) /CQ) . Applying the Bleistein inverse formula to this data, we obtain A

A 4 (Z) = : r ^ A /

r+oo

F{uj)-^

exp [-2iuj {z - h) /CQ] d (2a;/co)

(15.76)

Linear inversion methods

483 HB {Z

(15.77)

- h)

^0

where HB {Z — /i) is a bandlimited version of the Heaviside step function: 1

f^"^

HB{z-h) = - J ^ F{u;)

7T6 (UJ)

— ^ - exp [—iuj {z — h)] dw.

(15.78)

Differentiating function HB {Z — h) with respect to variable z, we arrive at the bandhmited delta-function (15.65) dHB {z-h)/dz

1 /"^^ = — F (uj) exp [-iou {z - h)] duj = 6B{z-h). 27r7-oo

(15.79)

Therefore, differentiating the function ASQ (Z) , we obtain dAs^B i^) /dz = — 2 ^ ^ ^ {z-

h).

(15.80)

Thus, the true reflection coefficient can be found from the formula ^2

"—i^m"^^'^'^""^.'^-

(15.81)

This result gives us an idea how one can modify the original inversion method to find the distribution of the reflecting boundaries and the corresponding reflection coefficients in the medium. Note also that differentiating the anomalous square slowness with respect to z in the Bleistein inversion formula (15.63) is equivalent to multiplying the scattered field in this formula by —2iuj/co. Based on these observations, Bleistein and his co-authors (Bleistein et al., 2001) introduced a bandlimited reflectivity function^ determined by the following formula: 2

4

PB (^) = - y — 4

r+oo

/

{2iuj/co)p'{0,u;)e-''-^/^du;

TTCQ J _ O O

+00

2

iujp'{0,u;)e-^'^'^''duj, (15.82) / -oo where they excluded a term 6B (0) in the denominator, assuming that the corresponding values of the reflection coefficient were scaled by the area under the filter F (UJ) . Indeed, applying formula (15.82) to the bandlimited scattered field (15.76), we obtain TTCI

PB {Z) = A6B

(z-h).

(15.83)

Formula (15.83) justifies the introduction of the reflectivity function. In the case of a multilayered model, this function represents a combination of the bandlimited delta-functions with the positions of extremum at the reflecting

Integral representations in wavefield inversion

484

boundaries, and the sizes of the peaks proportional to the corresponding reflection coefficients (see Figure 15-1). Thus, to transform the inverse formula for the anomalous square slowness into the inverse formula for corresponding reflectivity function distribution we have to multiply the observed scattered field by the factor {—2iuj/co) • (—Co/4) = iu;co/2. We will use this simple rule in subsequent sections on Bleistein inversion as well. 1-D inversion in a variable-background medium We would like to develop a similar approach to 1-D inversion in a medium with variable background wavespeed. Unfortunately, in a general case there is no simple analytical expression for the Green's function. However, when investigating a high frequency acoustic field, one can use the WKBJ approximation (13.83) for the calculation of the Green's function: 9b{^\^'^^)

= --—^—

exp

ZUJ

i:

dz

Cb{z)

(15.84)

Considering again "zero-offset data", / = Zj = 0,

we can write equation (15.57) as follows: r+oo

f{0,uj)=p'{0,0,uj)=u;'

/

As'{z)[g^{z,0,uj)]'dz.

(15.85)

Jo

In the last formula ^^^''''^^

=

^(^)

^

'JUJT{Z,0)

(15.86)

where A{z) = and T{Z.

-^Cb{z)cbiO) dz Jo Cb(z)

(15.87) (15.88)

Using notations (15.86) and (15.88), we can write (15.89) Taking into account that As^ (z) = 0 for z < 0, we can express the last formula in the form +00

/

-oo

4>{z')As'^ {z') e^^""'^'°Uz'

(15.90)

Linear inversion methods

485

where 4>iz') =

(15.91)

-\A\Z'),

and we use a new notation z' for the depth (this will be more convenient for the further algebraic transformations). Note that expression (15.90) has an appearance similar to a Fourier transform from the spatial z domain to the frequency LJ domain. Of course, this is only an apparent similarity, because we have a travel time r in the argument of the exponential function, instead of the depth, z. However, we can exploit this apparent similarity and construct an expression similar to the inverse Fourier transform (15.62): H-oo

/

p'{0,u;)e-^^^^'-°'^doj,

(15.92)

•OO

where ip (z) denotes the result of this "inverse Fourier transform." Let us substitute (15.90) into (15.92) and change the order of integration over the depth, z\ and frequency, LU: +00

/

/•+00

(l){z')As\z')

/

e'^'^^^^'^'^dwdz' ,

(15.93)

J —OO

OO

where

iz,z')= f

dz

J z'

We know that

f

g2ia,r(z,.')^ = 2TT6{2T{Z, Z')) .

(15.94)

J —c

Therefore equation (15.93) can be written as follows: +00

/

(t>{z')As\z')6{2T{z, z'))dz'. •OO

On the other hand, according to the property of the delta function, we have 6{2T{Z, z'))dz' = 5(2

r

- ^ )

dz'

By substituting equation (15.96) into (15.95), we obtain /• ++ 00 00

V'

/

(l){z')As'^{z')cb{z')6{z - z') dz' = TT4>{z)As'^{z)cb{z). -OO

(15.95)

486

Integral representations in wavefield inversion

From the last equation we have at once

^s'{z)= J['\

..

(15.97)

Substituting equations (15.91) and (15.92) into (15.97), we finally arrive at an explicit formula for the anomalous square slowness:

It can be easily shown that the same formula holds true for bandlimited wavefield data (Bleistein et al., 2001). Note also that for a model with constant background wavespeed, Cb{z) = CQ, the last formula naturally reduces to formula (15.63) derived in the previous subsection. As in the case of constant background velocity, we can obtain the inverse operator for the reflectivity function from formula (15.98) by multiplying the observed scattered field by the factor icjci, (z) /2 : (3{z) =

—--

/

iujp\{),uj)e-^'''^^'^^^du.

(15.99)

Thus while the inverse operator (15.98) determines the distribution of the anomalous square slowness within the medium, formula (15.99) solves for the reflectors and the corresponding reflection coefficients. In geophysical applications, for example in seismic methods, the reflecting boundaries are the main target of exploration. That is why inversion formula (15.99) plays an important role in the interpretation of seismic data. In the next section we will show that this method can be extended to a 3-D case. This technique provides the basis for modern methods of seismic data interpretation. T h r e e - d i m e n s i o n a l inversion Consider now a 3-D acoustic model with the arbitrary distribution of the wavespeed c (r) given by the formula c^ (r)

c^(r)

where Cb{r) is the background (normal) velocity distribution, and As'^{r) is some anomalous square slowness, which is unequal to zero only within the local anomalous zone D. The acoustic field in this model is excited by a point source located at r'. According to (14.3), this field satisfies the equation Vyr,a;) + — - p ( r , a ; ) - - ( 5 ( r - r O ,

(15.100)

Linear inversion methods

487

and the corresponding radiation conditions. The total wavefield in the model can be expressed, as usual, as the sum of the incident field p*(r,c<;) and the scattered field p*(r,cj) (see equation (14.9)), where the incident field p*(r,a;) is the solution of the following Helmholtz equation with background velocity distribution: V V ( r , Lo) + ^ p X r , a;) = - 6 (r - r').

(15.101)

Note that the incident field generated by a point source is nothing more than the corresponding Green's function for the model with background wavespeed: p\r,ij)

= p\r, r',uj) = G'"(r| r';oj).

(15.102)

The Bleistein method is based on a Born approximation of the scattered field, which has the form p'{rj,uj)=uj^

fff

G''{rj\r;uj)As\r)p\T,uj)dv.

(15.103)

Substituting (15.102) into (15.103), we have p'{rj,uj)^p'{rj,r\uj)=iu^

fff

As^ (r) G^(r,| r;cc;)G^(r| r';cj)dt;.

(15.104)

Taking into account the reciprocity theorem (14.24), we can write the last formula in the form p\rj,uj)=p\rj,r',uj)=uj'

III

As' {T)G^{v\Tf,u;)G^{r\v';uj)dv.

(15.105)

The solution of the inverse problem can be simplified for the case of "zero-offset data," when r^ = r': P'{V\LU)=P'{T\T\LU)=UJ^

fff

A52(r)[G"(r|r';u;)f dt;.

(15.106)

Thus, we arrive at integral equation (15.106) for the anomalous square slowness. The key to the solution of this equation for a general variable-background wavespeed is in using the WKBJ approximation (13.80) for the Green's function. By substituting (13.80) into (15.106), we finally find P'{T\UJ)

= u^ fff

As^ (r) A" (r, r') e'^'^^^'^'^di;.

(15.107)

This equation becomes relatively simple in a model with constant background velocity.

I n t e g r a l r e p r e s e n t a t i o n s in w a v e f i e l d i n v e r s i o n

Let us consider the case where Cb{r) = CQ. Then, according to (13.74), the Green's function is equal to G-(r,|

1 iuj\r-r'\/co 47r|r — r'l

r »

(15.108)

Therefore equation (15.106) takes the form 1

UJ"

fir^cu)

As'

(r) •

^2iuj\r-r'\/co^^^

(15.109)

Using the new notation TC = | r — r' |, we can rewrite the last equation: 1

^2iujrc/co

fir'^u)

As^ (r)

(4^)2

(jj^

—dv.

(15.110)

D

Differentiating both sides of equation (15.110) with respect the frequency a;, we have

> ( ^ »

^//X-^«

o2iu;rc/co

-dv.

(4

(15.111)

From the last equation we find at once -27rzco^OUJ

1 -y{r\cu)

III

As\v)g^{T,,u)dv,

(15.112)

UJ''

where

^2iojrc/co

g'^{x - x \ y - y^z-

(15.113)

z'.uj) = go{rc,uj) ATTTC

We assume now that the scattered field is observed on the horizontal surface z = 0 and the domain D with nonzero As^ belongs to the lower half-space. In this situation we can extend the integration volume on the right-hand side of equation (15.112) to cover the entire lower half-space 2: > 0 : -27rzco

_d_

1

duo

UJ^

V\x\y\{),u)

A5^ (x, 2/, z) g^ix -x',yy', z, u)dxdydz. (15.114) J J Jz>0 It is now clear that we have a convolution integral on the right-hand side of equation (15.114) over the horizontal variables x and y. The natural way of solving a convolution integral equation of this type is based on a Fourier transform over the spatial variables x and y.

Linear inversion methods

489

Following this idea, let us introduce a pair of forward and inverse Fourier transforms in the spatial domain: n r+oo

SfiK. ky) = II f{x, y) = ^

/(x, y)e-^'^^-^^^yy^dxdy,

ff " Sfih,

ky)e^'^^^^^^yy^dk,dky.

(15.115) (15.116)

Applying this Fourier transform to (15.112), we obtain 1 1 /*+°° —:^Sps{kx,ky,^,uj)\ = I SAs^{kx,ky,z)SgQ{kx,ky,z,uj)dz, (15.117) du /^ J Jo where S^s'^i^x^ ky, z) and SgQ{kx, ky, z,u)) are the spatial Fourier spectra of the functions As^ (r) and g^{rc,uj) respectively: —27rico

r

r+oo

SAs^kx^ ky^ z) = / /

As2(x, y,

z)e-^'^'^^^'yy^dxdy,

and r /•+00

Our next task is to find Sg^. It can be shown that this function satisfies a 1-D Helmholtz equation: Sgoikz^z) = -6{z)

^+(«.)'•

(15.118)

where kl = ujycl-kl-kl.

(15.119)

Thus SgQ is the 1-D Green's function for the medium with the constant velocity c=l/2: 1 ^2ikz\z\ Sg^ykz^z) (15.120) Substituting (15.120) into (15.117), we find -SnkzCo

d r 1 duo \-oSj,s{kx,ky,^,u)\

I r"^ =/

SAs^{kx.ky,z)e^'^^'dz,

(15.121)

where \z\ was replaced by z because the variable of integration is positive. Taking into account that As^ (x, y, z) = 0 for z < 0, we can present the right-hand side of the last formula as a Fourier transform from the z domain to the kz domain: -STTkzCo—OUJ

1

1

—Sj,s{kx,ky,0,uj)\ ^

/*"^°°

- / J

J —oo

SAs'^(kx,ky,z)e^'^^'dz.

(15.122)

Integral representations in wavefield inversion

490

Note, however, that we should be careful in selecting the proper branch of the square root function when calculating kz based on expression (15.119). The standard choice for k^ is always to select the real-valued and positive square root: kz = sign u^u'^/cl

- k^ - k^ for |a;/co| > ^A:^ -f k^,

(15.123)

and kz = i^kl

+ kl-u'^/cl

for \uj/co\ < ^Jkl + kl

(15.124)

The last formulae show that kz is real-valued only for |ct;/co| > yjk^ -f- k"^. The multipHer (sign a;) ensures that Sg^ is an "outgoing " Green's function, which satisfies the radiation condition at infinity. The imaginary values of kz in formula (15.124) correspond to so-called evanescent waves, which are usually neglected in analyses of geophysical data. At the same time, we see that the slowness distribution, As^ ( r ) , can be found from the purely real kz data. That is why we will assume that our data are free from evanescent waves. Under this assumption, we can invert formula (15.123) with respect to UJ : uj = u {kx, ky, kz) = cosign kzJkl

-f A;^ -h /c^, for \UJ/(^\ > Jkl

-h k^.

(15.125)

Formulae (15.123) and (15.125) show that the interval —oo < CJ/CQ < \/^i~-i-~^ of a; variations corresponds to the interval — oo < /c^ < 0 of/c^ variations, and, respectively, interval ^k^ + /c^ < (JO/CQ < -f-oo is translated into the interval 0 < fc^ < H-oo of kz variations. Therefore, we can use fc^ as a parameter conjugate to the variable z in a Fourier transform, and introduce another pair of forward and inverse Fourier transforms: +CXD

/

m

7r7-oo

me

•CXD

(15.126)

Sf{kz)e-'''''dkz.

Applying the inverse Fourier transform (15.126) to (15.122), we write + 00

S^s'^{kx,ky,z)

= -8co

k. /

-OO

d_

1

bps U' 2-p-

-2ikz

'dkz

(15.127)

Substituting (15.127) into (15.116), we finally arrive at the solution of our inverse problem: As^{r)=As'^{x,y,z) = +00

/»+00

'-^ICL

O

kz~^

2^p^ \kx) ky^v, UJ) e^'if^^^^f^yy^e-^^'^'dkxdkydkz.

LU

(15.128)

Linear inversion methods

491

According to the basic principles of the Bleistein method outhned above for a 1-D case, we can introduce a 3-D reflectivity function by multiplying the observed scattered field by the factor iujco/2 :

+CXD

/•4-00

O

-'^II:L?^

duj

1

e^^(f^^-+^yy)e-^'^^'dk,dkydk,.

(15.129)

UJ

Note in conclusion that a similar result can be derived for a variable-background velocity model. I refer readers to the text by Bleistein et al. (2001) for more details about this method. 15.1.5 Inversion based on the Kirchhoff approximation As we have discussed already in the previous section on Bleistein inversion, the goal of inversion in many applications is the reconstruction of the reflectivity function, which provides information about the distribution of the reflecting boundaries in the medium under investigation. This problem can also be addressed by using the Kirchhofl" approximation (14.66), which we reproduce here for convenience:

p-(r',.) = K(A) = / / ^ A ( r . ) ^ M £ - ^ > | l M £ ^ . . . ,

(15.130)

where n is upward-directed normal to 5 , and i^(A) denotes the corresponding integral operator applied to the reflection coeflHicients A ( r ^ ) . In particular, we can select as surface B an arbitrary horizontal plane, B^ {z = h). in the lower half-space. Then we have p'{r ,uj) = -

A{rB)-^-^

g^

-dsB,

(15.131)

where we take into account that d/dn= —d/dzB, assuming that the axis z is directed downward. The Kirchhoff approximation represents a linear Kirchhoff integral operator K{A) with respect to the unknown reflection coefficients A (r^): da=p'{r\Lj)=^K{A),

(15.132)

where da denotes the observed scattered field. In an inverse problem we assume that the scattered field, da = p^, and the background wavespeed distribution are given. The goal is to find the reflection coefficients A ( r ^ ) . As a result, we have a linear equation with respect to the unknown function A (r^) • As in the Born inversion method, we solve the linear inverse problem (15.132) by regularization methods, imposing additional conditions on the class of inverse models with corresponding stabilizing functionals.

492

Integral representations in wavefield inversion

We introduce again a Hilbert space Da of the acoustic field data, given on the surface of observations S with the metric determined by formula (15.10) and a Hilbert space M of reflection coefficients A, given on the surface B with the metric (A(i), A(2))^ = 1 ^ AW (ra) A(2) {rs)dss.

(15.133)

The corresponding parametric functional for Kirchhoff inversion with a minimum norm stabilizer has the form P^ = {K{A) - da, K{A) - da)o^ +a {Wa (A'-A^p,) , Wa {K^-hl^))^

- min,

(15.134)

where da = p^ are the observed acoustic wavefield data, Aapr represents some a priori values of the reflection coefficients, and Wa is the weighting operator. Note that we use the minimum norm stabilizer in this example just for convenience. In principle, any stabilizer introduced in Part I can be apphed. The minimization problem (15.134) is solved, for example, by the conjugate gradient or re-weighted conjugate gradient methods introduced in Chapter 5. One possible iterative scheme of the conjugate gradient method is given below: Rn = K{An) - da,

(a)

/ - = P ( A , ) =K^ {Rn) +aW'W{An-Aapr),

{b)

k^ is determined from the minimization problem P - (An+i) - P - ( A , - k^l^^ = min.

(e)

For practical realization of this iterative inversion scheme we have to flnd the explicit expressions for the adjoint Kirchhoff operator K* according to formula {KiA),q)j,^

= {A,K*{q))^.

(15.136)

Using the definitions (15.10) and (15.133) of inner products and expression (15.130) for the forward operator, we can rewrite formula (15.136) as {K (A) , q)j^^ =Re

f [[ K (A) q^'dSdu;

Linear inversion methods

493

= [ //Re / / A M g[p'(rB,a>)G"'(r|rB;a>)] dsBq* (r, uj) dsdu) On = ff A{rs)Re\f ff d\p'*{vB,co)G^*{T\rB;iv)] q (r, uj) dsdu lUJJs dn JJB IJnJJs =

dsf

(15.137)

{A,K*iq))^.

From the last formula we have A, K'

5[p-(r5,a;)G-*(?|rB;a;)] q (r, (jj) dsduj 9n

{q) -Re f ff Jn JIsJs

-0.

(15.138)

M

Equation (15.138) should hold for any A, for example for A = K^ (q) - Re

q{r,uj) dsdu.

dn

JnJJs/s Therefore,

K*{q)-Re f ff —-—^ JnJJss

^-^

—q (r, u) dsdu

dn

0. M

From the last formula we conclude that the adjoint Kirchhoff operator K^ is given by the formula K"^ {q) = Re Jn J Js

9[p-*(rB,a;)G"*(?kB;^)] g(r,a;) dsdu. dn

(15.139)

Assuming that we integrate over a symmetrical frequency interval f2 (for example, from — oo to -hoc), we can drop the symbol Re in the above formulae because the imaginary part of the integrand is anti-symmetric in UJ. Note also that we calculate the normal derivative in expression (15.139) with respect to the variable r^ G B, while the integration is conducted along the variable r on the surface of observation S. Therefore we can take the normal derivative operator outside of the integration symbols: K' (g) = ^

y p'^iTB.u;) 11 G-*

{VIVB; UJ)

q (?,u) d^du.

(15.140)

According to formula (15.32) we can introduce an auxiliary acoustic field p^* (r,u;) : p^* (r,a;) - / / G^* (?| r]uj)f

{r,uj) d:S.

(15.141)

Substituting equation (15.141) into (15.140), we obtain K* ( / ) = - ^ 1 p'*{vB, co)p^* {TB,CO) du.

(15.142)

494

Integral representations in wavefield inversion

Finally, if surface B coincides with the horizontal plane, Bh {z = h), we have K*{f)

= -^^Jp'*{vB,uj)p''*{rB,oj)(Lj,

(15.143)

where we take into account that d/dn= —d/dz because the axis z is directed downward. Note that, as in the case of the Born imaging method, jwe can use the convolution theorem to write the expressions for the auxiliary field P^ in the time domain as follows: i'^(r/) = ^ y

/ / . /

P^*(r,a;)e--''dc;

G"" (r', t'\T, t)p' (r, t) dtds,

(15.144)

where (5"(r,t| r ' , 0 = G " ( r , - t | r ' , - 0 . Therefore, the auxihary field P^ (r,f) is the back-propagated field, and the adjoint Kirchhoff" operator K* for a residual field can be calculated by back propagation of the residual field. We will show in the next sections that this procedure is similar to Stolt's Fourier-based migration transformation. Thus, this result demonstrates that migration is similar to applying an adjoint Kirchhoff operator to observed scattered wavefield data. In the general case of an iterative inversion scheme (15.135), we should apply the adjoint Kirchhoff operator K^ to the residual field Rn = K{An) — da computed on each iteration. Each of these applications is equivalent to the standard migration of the residual field. Thus, iterative Kirchhoff inversion can be treated as an iterative migration algorithm. 15.1.6 Traveltime inverse problem The traveltime inverse problem provides another typical example of application of the linearization method to the solution of inverse problems for wave phenomena. In the framework of the "geometric optics" approach to seismic problems, the traveltime T(r', r^) of the seismic ray can be related to the local seismic velocity c(r) by the relationship r(r',r,)=/ -^, (15.145) A(r',r,) C(r) where L{r',rj) denotes the raypath between source r' and receiver r^. The last equation can be written in another form, using slowness s(r) = l/c(r): r(r',r,) = /

s{r)dl.

(15.146)

Linear inversion methods

495

Note that this equation looks Hnear, but is actually nonlinear, because the raypath L{r\rj) also depends on the slowness 5(r). However, we can calculate the variation of the traveltime using Fermat's principle, which states that the traveltime is stationary with respect to a variation of the raypath L{r',rj) : 5r(r',r^)- /

6s{r)dL

(15.147)

Let us suppose that we know some background model ^^(r) of the slowness distribution and that the current model S{T) is obtained by a small perturbation of Sb{r): s{r) =Sb{r)-\-As{r). Denoting by r^ the traveltime along the unperturbed raypath Lt in the model 55 (r), we can use (15.147) as a linear integral equation for the difference A5(r) with respect to the difference of traveltime AT = T — r^: Ar(r',r^) = /

As{r)dL

(15.148)

Note that this equation has been obtained using Fermat's principle in which the result of a raypath perturbation is set equal to zero. Equation (15.148) can be considered the linearization of the nonlinear equation (15.146). For a whole vector of A^ travel time observations, one obtains a system of linear integral equations which after quadrature discretization gives a linear system of equations for slowness perturbations As. For practical applications a rectangular grid is generally chosen and the wave slowness is discretized by treating it as constant in the cells determined by the grid. Thus one forms a vector s of L slowness perturbations: S = (A5i,A52,

ASL)-

Usually the first arrival traveltimes r for waves propagated between sources and receivers are also considered. Differences of N travel times for theoretical models and actual observations produce the vector d of the data: d = (Ari,Ar2,....ATAr). Thus we obtain a linear system of equations for slowness perturbations: d = As.

(15.149)

Now we can use any of the algorithms described in Chapters 3 and 4 for regularized inversion of this matrix equation (15.149).

496

15.2

Integral representations in wavefield inversion

Quasi-linear inversion

In full analogy with the electromagnetic case, we can apply the quasi-linear (QL) approximations introduced in Chapter 14 for acoustic and vector wavefield inversion. We begin our discussion of the basic principles of QL inversion with the simpler scalar case of acoustic waves. 15.2.1 Quasi-linear inversion of the acoustic wavefield In the QL approximation, we assume that the scattered acoustic field p^ inside the inhomogeneous domain D is linearly proportional to the incident field p^ through some reflectivity coefficient A: P\Y,UJ)

= \{r,uj)p\r,uj).

(15.150)

At the same time, according to (14.38), p^ can be calculated as an integral involving A and the incident field: p'{Tj,uj) = u^G^ { A52 (r) [1 + A(r, UJ)]P\Y, U)] .

(15.151)

Following the basic principles of electromagnetic QL inversion, we can introduce a material property parameter m, equal to m(r, iu) - As^ (r) [1 -h A(r, u)].

(15.152)

Substituting expression (15.152) into (15.151), we arrive at a linear equation for the material property parameter m: P'{TJ,UJ)

= u^G^ [m{Y,u)p\Y,u)\

.

(15.153)

The reflectivity coefficient A can be determined from the following linear equation inside the inhomogeneous domain L>, as long as we know m : \{Y,U)P\Y,U)

= u'^Gy, [m{Y,uj)p\Y,u)\

.

(15.154)

After determining m and A it is possible to evaluate the anomalous square slowness distribution As^ (r) from equation (15.152). Note that equation (15.152) should hold for any frequency, because in a general case the reflectivity. A, and the material property parameter, m, are functions of frequency as well. In reality, of course, it holds only approximately. Therefore the square slowness As^ (r) can be found by using the method of least squares to solve equation (15.152)) with respect to A52(r): ||m(r,a;) - As^ (r) [1 + A(r,cj)]||^ = min, (15.155) where the L2 norm

||...||Q

is determined as

.-)L = y ^

l ^ ( r , ^ ) L = \\ I |m(r,a;)|^cL;. In

Quasi-linear inversion

497

In the case of high frequency asymptotics, we can select the reflectivity coefficient to be frequency independent, A = A ( r ) . In this case the material property parameter is also frequency independent, m(r) = As'^ (r) [1 + A(r)],

(15.156)

and equation (15.153) takes the form p'ivj.ij)

= cu'^G^ [m(r)p'(r,a;)] .

(15.157)

Now we can find A(r) by solving the minimization problem: \\X{r)p\rj,u;)

- uj^G^ [m{r)p'{r, u;)] \\^ = min.

The anomalous square slowness is then determined by the simple expression A,= (r) = ^

.

(15.158)

Note that the quasi-linear inversion method, outlined above, can be easily extended to the case of the vector wavefield. We leave the detailed derivation of the corresponding formulae (which look very similar to the analogous formulae for electromagnetic field QL inversion) as an exercise for the interested reader. 15.2.2 Localized quasi-linear inversion based on the Bleistein method We have noticed already in electromagnetic sections of the book that the quasilinear inversion, introduced above, cannot be used for interpretation of multi-source data, because both the reflectivity coefficient A and the material property parameter m depend on the illuminating incident wavefield. However, in many geophysical applications, for example in seismic exploration or in cross-well tomography, the data are collected using moving transmitters. In this case one can build an effective inversion scheme based on the localized quasi-linear approximation introduced in Chapter 9, which is source independent (Zhou and Liu, 2000; Zhdanov and Tartaras, 2002). Indeed, we can introduce a new function, ruL (r) = A^' (r) [1 + AL(r)],

(15.159)

which we call a localized material property parameter. Note that the function TUL (r) is independent of the source position because we select a localized reflectivity coefficient XL which does not depend on the incident field (and frequency). This is the main difference between the localized QL inversion and the original QL inversion considered in the previous section. We assume now that the scattered acoustic field p^{rj,cj) (generated by a source with one or different positions) is measured at a number of observation points.

498

Integral representations in wavefield inversion

Tj. Using LQL approximation for the observed field, p*(rj,a;), we arrive at the following equation, similar to (15.157): p'{rj,uj) = uj^G^ {rriL

{T)P\V,UJ))

,

(15.160)

which is linear with respect to the localized material property function TTIL ( r ) . The main difference between equations (15.160) and (15.157) is that in the LQL approximation, the material property function rriL (r) does not depend on an illuminating source. We can solve the linear equation (15.160) with respect to rriL ( r ) , which is source independent. Now, a scalar reflectivity coefficient \L (r) can be determined, based on condition (14.54): \\XL{V) - J'G^

K ( r , u ; ) ] | | ^ - min .

(15.161)

Knowing A^ (r) and m^ (r), we can find As^ (r) from equation (15.159):

This inversion scheme can be used for a multi-source technique because A^ and vfiL are source independent. It reduces the original nonlinear inverse problem to the same Born-type inversion we discussed in the previous sections. That is why some methods of Born-type inversion can be applied in this situation as well. In particular, the Bleistein method can be very useful in combination with localized QL inversion. For example, in the case of zero-offset data, we can use a formula similar to (15.106) to describe integral equation (15.160) v\r'M

= c^' / / / m^ (r) [G^(r| x'',uj)fdv.

(15.163)

Thus, we arrive at integral equation (15.163) for the material property function TTiL (r). Applying the Bleistein method to the solution of this equation, we use the WKBJ approximation (13.80) for the Green's function. By substituting (13.80) into (15.163), we obtain V\v\u)^ij

fff

mL{r)A^r,r')e^'^^^'^'''^dv.

(15.164)

We can now use the standard technique of Bleistein inversion to find mi ( r ) . The corresponding formulae become especially simple for a model with constant background wavespeed. In this case, according to equation (15.128), we have rriLir)

_8co fr°°+ 00 f/•+00

n' JJ-^ J-.

^ duo

=mL{x,y,z)

^«bps(A:x, r!,yi 0 , to) enk'^+''yy)e-^'''^'dk^dkydk„

(15.165)

Nonlinear inversion

499

where Sps{kx, ky, 0,uj) is the spatial Fourier transform of the observed scattered field on the ground: p r+oo

Sps{k:,,ky,0,uj) = / /

p'{x,y,0,Lj)e-'^'^''^''-^^yy^dxdy,

(15.166)

Thus, localized QL inversion is reduced to Bleistein inversion with respect to the material property function TTIL (X, y, z) and then to a simple correction of this inversion result by solving the minimization problem (15.161) and applying the algebraic transformation (15.162). Localized quasi-linear inversion increases the accuracy and efficiency of wavefield data interpretation because it is based on a much more accurate forward modeling solution than the Born approximation, used in the original Bleistein method. An example of successful apphcation of the localized QL approximation in radardiffraction tomography can be found in (Zhou and Liu, 2000). 15.3

N o n l i n e a r inversion

15.3.1 Formulation of the nonlinear wavefield inverse problem In this section, we consider general nonlinear acoustic and vector wavefield inverse problems in parallel : da = Aa{As^), (15.167) d^ =- A^ {As^) ,

(15.168)

where Aa and A^ are nonlinear forward modeling operators, given by the formulae (15.3) and (15.4): ^ . ( A ^ ^ ) = uj^G^ {As' {r)p{T,Lu)) ^cj^ fjf A, {As') =u'G^[As'u]

=u'

fff

G^{vj\r;uj)

G]S!{rj\T;uj)As'{r)p{T,u;)dv, ' As' {r)u{r,u)dv,

(15.169)

and da, d-^ represent the scattered acoustic or vector wavefield observed on the surface S. In the last formulae G'^{rj\r]Lu) and G^(rj|r;a;) stand for the Green's function and tensor respectively, defined for an unbounded medium with background slowness 56.

The nonlinearity of the operators Aa and A^ with respect to the anomalous slowness As' is related to the fact that the wavefields p(r, uj) and u(r, uj) also depend on As'. As usual, following the conventional Tikhonov regularization method, we substitute for the solution of the inverse problems (15.167) or (15.168) a minimization of the corresponding parametric functional with, for example, a minimum norm stabilizer: P - {As')

= {AaA^s')

- da,.,

AaA^s')

- <„)^^_^ +

500

Integral r e p r e s e n t a t i o n s in wavefield i n v e r s i o n

a {W (As2-A5^^) , W ( A s ^ - A s ^ ^ ) ) ^ ,

(15.170)

where W is the corresponding weighting operator, and the inner products (.., .)D and (.., .)^ are determined by formulae (15.10), (15.11) and (15.12). The minimization problem (15.170) can be solved by the conjugate gradient or re-weighted conjugate gradient methods introduced in Chapter 5. Let us describe, for example, the algorithm based on the regularized conjugate gradient method (5.92), which we reproduce here with small modifications^:

/ - = r ( A 4 ) =F*,^„ (i?„,„;„) +aW*W{/\sl-^sl^),

[b)

PZ = \\lnf I l I C i II'.

(c)

AsU,

t-l^^^Tu-x.

II = 1%.

= Asl-~k^t

(d)

(15.171)

k^ is determined from the minimization problem P - ( A 4 + : ) = P " {^sl

- k±)

= min,

(e)

where, as usual, we denote by A5^ the corresponding anomalous slowness distribution at the n-th iteration. We can see that practical implementation of this algorithm requires computing the Prechet derivative operator Fa,v;n of the corresponding forward modeling operators Aa^y on each n-th iteration. 15.3.2 Frechet derivative operators for wavefield problems In Chapter 14 we developed expressions for the Frechet differentials of the forward modeling wavefield operators (14.29) and (14.84), which we reproduce here for convenience: Fa{s\6s^) = Fair \ 6s^) = Lu^ fff

G'^{r\T;uj)6s^r)p{r,uj)dv.

(15.172)

Fy{s\8s^)=Fy{r\Ss^)=uj^

G'^i^lv.u;) - 6s\r)u{r,uj)dv.

(15.173)

and fff

^For conciseness, the scalar acoustic and vector wavefield equations have been combined into single equations in this summary, with alternative subscripts a and v to distinguish them; it is understood, of course, that d^, Ay, and Ky-n (but not the adjoint operator F*.^), would be printed in boldface in separate vector wavefield equations since they refer to 3-D vector elements in D^.

Nonlinear inversion

501

In these formulae, 5^ = 5^ + As'^ is the square slowness model for which we calculate the variation of the forward modeling operator; ^5^ is the corresponding variation of the square slowness 5^, which is obviously equal to the variation of the anomalous square slowness, 6s'^ = SAs"^; expressions G^(r|r;a;) and G ^ ( r | r ; a ; ) stand for the Green's function and tensor defined for the given square slowness 5^; and the function p{r,u) and vector u(r,a;) represent the total acoustic and vector wavefields for the given square slowness 5^. Note that in the RCG algorithm (15.171), the expression F*^.^ (Ra^y-n) denotes the result of an application of the adjoint Frechet derivative operator to the corresponding acoustic or vector residual field Ra,v;n = ^a,v(A5^) — da^y on the n-th iteration. The expressions for the adjoint Frechet derivative operators are given by formulae (15.22) and (15.23). Based on these formulae, we can write FL (Ran) = f ojY^ir, uj) ff GT (r| r; uj) Ran (r, a;) dScLj, Jn J Js

(15.174)

and

FL {Rvn) = / a;2u;(r, a;) • / / GZ* (?| r; a;) • R.„ (?, UJ) dScLj. Jo. J Js

(15.175)

where Pn(i*,^) = p'(r,a;) -f-p^(r,(x;) and Un(r,a;) = u'(r,(x;) -f u^(r,u;) are the wavefields, and G5^(r|r;a;) and G'^(T\Y,UJ) are the Green's functions and tensors, computed on the n-th iteration for the square slowness distribution, s^ (r) = s^^ (r) -|As^ ( r ) . Here the asterisk * denotes the operation of complex transposition. According to reciprocity theorem 27 (formula (14.24)), the surface integral terms in the last formulae can be treated as the complex conjugate acoustic, p ^ (r,(x;), or vector u^ (r,u;), wavefields, p^ {v,w) = J J G™ (r| r; uj) K „ (?, LJ) dS,

(15.176)

11 G - (r| r; cj) • R : „ (?, a;) ds,

(15.177)

u^ {r,oj)

generated at the point r inside the volume V, by the fictitious scalar or vector wavefield sources, distributed along the surface of observation S with the surface density determined by the complex conjugate residual fields: / ! „ ( ? , a;) =.i?:„(?,u.),

(15.178)

and (f>'snir,^) = Knir,^)-

(15.179)

502

Integral representations in wavefield inversion

Therefore, from (15.174), (15.175) and (15.176), (15.177) we find that the result of applying the adjoint Frechet operator to the residual field is equivalent to integrating over the frequency range Q of the product of the complex conjugate reciprocal fields and the corresponding wavefields p*(r,u;) and u*(r,c<;), computed at the n-th iteration: Kn (Ran) = J ujY„ir, u;)p^* {v,uj) dw, Jn

(15.180)

F:„ (R„„) = / a;2u:(r, LJ) • u f (r,u;) cko. (15.181) Jet Thus, we can conclude, that the calculation of the terms F^^ {Ran) a-nd F^^ (R^vn) requires just one additional forward modeling solution on each iteration for the auxiliary wavefields p^* (r,6t;) and u^* (T^LJ) due to the reciprocal sources. Note that, as in the case of the Born imaging method, we can use the convolution theorem to write the expressions for the auxiliary field P^ (r,t) in the time domain as follows:

-ILL

"

G^ (r, t|r, t') Ran (r, t) dtds,

J JS J -oo

(15.182)

where 1

/'"^°°

Ran (?, t) = —

Ran {v,Uj)

e-'^'diJ

^7T J~oo and

_

G"(?,^|r,0 = G"(?,-t|r,-0Therefore, the auxiliary residual field P^ (r,^') is the back-propagated residual field! Taking complex conjugate of both sides of equations (15.180) and (15.181), and noting that the adjoint operators transform residuals into the real functions of square slowness distribution, we write these formulae in equivalent form: Kn (Ran) = / uj^Pn{r,u;)p^ (v^u) dw, Jn

(15.183)

F:^ {R,n) = / uj^Un{r,uj) . u ^ (r,a;) dou. (15.184) Jn We can also write the result of applying the adjoint Frechet operator to the residual field as a correlation of the time derivatives of the back-propagated residual field and the corresponding wavefield Pni'r, t) computed at the n-th iteration: Kn

{Ran)

=

/ UJ^Pn{r, Uj)p^ (r,Cj) diJ =

Jn

[-iujpn{T,

J—oo

Uj)] [-2CJp^* (r,Cj)] * doj

Principles of wavefield migration

503

Pn{r,t)P„{r,t)dt.

(15.185)

•CX)

Thus we can see that the adjoint Frechet derivative operator for the residual field can be found by back propagation (or "back projection," see Devaney (1984)) of the residual field. At the same time, this field can be obtained by using forward modehng calculations, because the wave equation is symmetric with respect to time (Tarantola, 1987; Mora, 1987). In the next section we will see that this backpropagated field can be reconstructed using migration of the residual data. However, in the inversion algorithm we should apply this migration iteratively on each step of the inversion. 15.4

Principles of wavefield migration

Migration algorithms are widely used in seismic exploration for imaging geological cross-sections. They became first known in the middle of the twentieth century as a facility for compensation of the structural distortions occurring in seismic time cross-sections as a result of seismic drifts.^ However, their utilization in seismic data processing became really widespread only in the late 1960s and 1970s, after the petroleum industry developed the corresponding technological basis required for migration imaging. Now seismic migration is one of the most powerful tools in the interpretation of seismic data. Originally, migration and inversion were treated as two diff^erent approaches to interpretation of the wavefield data. Traditional migration solves only for the traveltime parameter, thus providing the geometrical information about the position of the reflectors. In traditional inversion one searches for the physical properties of the medium, including its elastic parameters and wavespeed distributions, which requires a careful analysis of the true amplitudes of the wavefield. It has been demonstrated, however, that traditional migration can be treated as a first iteration in the full solution of the waveform inversion problem (Tarantola, 1987; Mora, 1987). This result is similar to one formulated in this text for potential (Chapter 7) and electromagnetic (Chapter 11) fields. Thus, modern developments in theoretical geophysics have led to dissolving the diff'erence between these two approaches to interpretation of seismic data (Bleistein et al., 2001). In this section of the book I will discuss the basic ideas underlying the principles of wavefield migration, and will show how these principles are related to the general inversion technique developed in the previous sections. 15.4-1 Geometrical model of migration transformation I will illustrate the basic idea of "geometrical migration" by considering a simple example of zero-off'set scattered wavefield data generated by an impulsive source and ^See the pioneer paper by Hagedoorn (1954), which explicitly formulates the migration problem, specifies some approaches to its solution, and introduces the terminology which is now widely used in geophysics.

504

Integral representations in wavefield inversion

observed on a horizontal surface of the earth. We assume that the earth can be treated as an acoustic medium with layered structure and constant wavespeed within each layer. We assume also that we can apply the "geometrical optics" approach to our data (which means that we restrict our analysis to high frequencies in this example). The observed seismic data within this model can be presented as the so-called "time cross-section," which represents the ensemble of responses from all source-receiver pairs. In the plot of time cross-section, the horizontal axis shows the location of each source-receiver pair, while the vertical axis corresponds to traveltime. In other words, we draw the seismic oscillations ("seismic trace") along the vertical lines going just beneath the corresponding source-receiver pair. Each pulse-like signal on a seismic trace (called a seismic "event") corresponds to a specific reflector in the medium. The continuous correlated distributions of the events in this plot are usually identified with the respective reflecting boundaries. However, the corresponding time shown on the vertical axis represents two-way traveltime. This is the time required for the pulse to travel from the source to the reflector and back to the receiver position, which coincides with the source position. Therefore, the depth to the reflector can be calculated by multiplying the traveltime by half the corresponding wavespeed. Applying this procedure to each seismic trace (which is equivalent to scaling the time on the time section by half the wavespeed), we arrive at the depth section, which could provide us with geometrical information about the location of the reflecting boundaries. The question is how accurate and reliable is this information? Unfortunately, this simple imaging technique works well only for horizontal or quasi-horizontal structures, because it is based on a simple model of vertically propagating waves and completely neglects the direction of their approach to the surface of observation. In seismic cross-sections, this presumption inevitably leads to considerable distortions of the curvilinear and steep-inclined reflecting boundaries which are typical for geological structures (Figure 15-2). We call this phenomenon a "seismic drift." The flrst migration algorithms were developed speciflcally to compensate for these distortions and to reconstruct the correct conflgurations and positions of the real reflecting boundaries in the seismic sections. Let us consider these distortions in detail. We denote by to (x) the traveltime of a wave travelling along the ray normal to the reflecting boundary, where x is the horizontal coordinate of the associated source-receiver pair. However, in the time-section this traveltime is plotted along the vertical axis, whereas in reality it may correspond to a propagation of the wave at some angle to the horizontal axis. It means that the reflection point may lie on a semicircle centered on the sourcereceiver position with radius equal to to{x). In a seismic cross-section, obtained by scaling the time by half the wavespeed c, the corresponding reflection point is located on a semicircle with radius equal to to (x) c/2. Note that the actual reflector must be tangent to this semicircle, so that the ray trajectory is orthogonal to the reflector (Figure 15-2). Considering the full ensemble of the source-receiver pairs and

Principles of waveReld migration

505

Figure 15-2 An example of the distortion in seismic section caused by a "seismic drift." a) The true reflecting boundary is shown by the bold solid line, b) The corresponding synthetic zero-offset seismic section. The image of the reflecting boundary is shown by the bold solid line. In order to generate a true reflecting point one should move the observed event (shown by squares, circles, or triangles) along the corresponding semicircles centered in the source-receiver positions with the radius equal to to (x) c/2. Note that the actual reflector must be tangent to this semicircle, so that the ray trajectory is orthogonal to the reflector. The actual reflecting boundary is drawn as an envelope of the ensemble of these semicircles.

506

Integral representations in wavefield inversion

Figure 15-3 An example of graphical migration. Reconstruction of the reflecting boundary by drawing an envelope of the ensemble of semicircles centered on the sourcereceiver positions with radius equal to to (x) c/2. The dotted line shows the position of the actual reflecting boundary.

the corresponding semicircles, we arrive at the conclusion that the actual reflecting boundary may be drawn as an envelope of the ensemble of these semicircles. Thus, to restore the correct configuration of the reflecting boundary in the seismic section, it is necessary to transform the position of each reflection point by moving it along the corresponding semicircle. Figure 15-2 illustrates the simple graphical method of this transformation. Following the previous discussion, we construct an envelope of the ensemble of semicircles centered on the source-receiver positions with radius equal to to (x) c/2 (Figure 15-3). The procedure just described represents the most elementary version of seismic migration, developed in the pioneer paper of Hagedoorn (1954). However, it demonstrates the basic idea behind the migration method: the geometrical images of seismic events are moved graphically or migrated from their original false positions in seismic cross-sections to the correct positions corresponding to true reflecting boundaries. In modern migration methods this fundamental idea is implemented using the sophisticated methods of ray tracing and reverse time wavefield extrapolation. Nevertheless,

Principles of wavefield migration

507

migration can still be treated as an imaging technique that provides a geometrical image of the earth's interior. 15.4-2 Kirchhoff integral formula for reverse-time wave equation migration We see from the previous discussion that the early migration schemes completely ignored the true amplitudes of the wavefield and their variation with time (the so-called dynamic characteristics of wavefield propagation). Introducing the dynamic characteristics of the seismic wavefield in the migration procedure has become a strong catalyst for the subsequent development of seismic interpretation methods. Remarkable achievements were made in this respect back in 1970-1972 by John Claerbout, who took advantage of the transition from pure kinematic models of wavefield propagation to a description of their dynamic characteristics with the help of the wave equation, and formulated the well-defined principles of seismic imaging (Claerbout, 1970, 1976, 1992; Claerbout and Doherty, 1972). Those principles were subsequently developed by Loewenthal et al. (1976) and many other researchers. In formulating the basic ideas of wave equation migration, I will follow mostly the paper by Loewenthal et al. (1976) and the text by Zhdanov et al. (1988). Let us consider some hypothetical wave traveling upwards to the surface of observation with its front coinciding in space at a certain moment of time (say at t = 0) with one of the reflecting boundaries. Then for time t > 0 this wave will obviously travel along the same rays as the actual zero-off"set scattered wavefield, refiected from the same boundary. This happens because, in a zero-oflPset case, the incident and reflected wave rays must coincide, which is possible only if the corresponding ray trajectory is normal to the reflecting surface. Thus, the actual zero-offset scattered wavefield and the hypothetical wave introduced above (which we call the Claerbout upgoing wave) would travel along the same raypaths. However, the Claerbout wave will take half the time to reach a certain point on the observation surface as will the actual scattered wave which has to travel the same distance twice. At the same time, if we assume that the Claerbout wave velocity in the medium is twice as slow as that of the real wave, then both waves will reach any point on the surface of observation simultaneously. These considerations form the basis of Claerbout's principle of seismic imaging, according to which the wavefronts of these fictitious upgoing waves at time t = 0 may be viewed as an image of the respective reflecting boundaries. Thus, according to Claerbout, the process of migration includes two elements: 1) backward extrapolation of the scattered wavefields (i.e. continuation of the waves in the direction opposite to that of their actual propagation) and 2) synthesis of the medium image as a snapshot (at time /; = 0) of the spatial structure of the wavefield produced by backward extrapolation. These principles form the foundation of the majority of algorithms of time section migration (Berkhout, 1980, 1984; Claerbout, 1985; Kozlov, 1986; Gardner, 1985). As we see from the previous discussion, the key element of this approach is the process of backward continuation of the scattered wavefield. The diversity of available migration algorithms is related mostly to the differences in the numerical

508

Integral representations in wavefield inversion

schemes of this continuation It is known from Chapter 13 that the accurate mathematical expression of the Huygens-Presnel principle in the wavefield case is the Kirchhoff formula (13.112) which makes it possible to continue (extrapolate) the wavefield from a certain surface S into the domain of space under examination. Subsequently, we generalized this formula for unbounded domains (13.196) using the Green's function for the wave equation in infinite space, characterizing the divergent spherical wave at infinity according to radiation conditions (13.163). With this choice of Green's function, the wavefields determined by equations (13.196) or (13.197) describe waves outgoing from the surface S in both upward and downward directions. According to the basic principles of migration transformation, outlined above, the backward extrapolation of the observed scattered field requires the reconstruction of upgoing Claerbout waves, i.e. waves traveling towards the surface of observation in the lower half-space. These waves should not satisfy the traditional Sommerfeld radiation condition at infinity, but rather the radiation condition (13.162) which corresponds to convergent waves. We can repeat the derivation of the Kirchhoff formula for an unlimited domain with the convergent wave radiation condition. In this case, however, we should use as the fundamental solution of the corresponding wave equation the function adjoint to the Green's function of the wave equation, G^(r,t| r ' , 0 = G ^ ( r , - t | r ' , - ^ ,

(15.186)

to characterize the propagation of convergent spherical waves. In the case of constant wavespeed c, the adjoint Green's function takes the form

which, according to (13.158), describes a convergent spherical wave, i.e. one arriving from infinity. In this case the wavefield represented by the corresponding Kirchhoff-type formula will be a combination of waves traveling towards the surface of observation as well: f/-(r', t) = 11 j ^ \G- (r', t'|r, t) ^U{r,

-U{T,t)^G^

{T\t'\T,t) dtds,

t)

T'eV"",

(15.188)

where [/^ is the migration field obtained as a result of backward extrapolation (continuation) of the Claerbout upgoing wave, 1/+ denotes the lower half-space, and the direction n in the differentiations is the upward pointing normal to the surface S. Thus, the selection of the fundamental solution in integral formula (13.196) determines the direction of propagation of the waves specified by that relation.

Principles of wavefield migration

509

A similar result takes place in the frequency domain. Substituting the Green's function G'^ (r'|r;a;) in equation (13.197) with its complex conjugate, G^* (r'|r;a;), which satisfies the radiation condition for convergent waves, we arrive at the following migration formula in the frequency domain: ^^^(r^a;)

G'"*iv'\r-oj)-^uiT,u) //.

-u(r,a;)|^G-*(r'|r;a;) ds,

r'ev-^

(15.189)

I would like to point out that expression (15.188) may also be obtained from the original Kirchhoff formula (13.196) by replacing conventional time t by reverse time T = T — t, where T is time of observation of the^avefield on the surface of the earth. Indeed, let us introduce an auxiliary wavefield C/(r, r ) = C/(r, T — r) and write the Kirchhoff integral (13.196) for this field as follows: [/^(r^r) =

G-(r',r'|r,r)—f/(r,r) J Js J~c

-[/(r,r)—G-(r',r'|r,r)

dr (is.

v' e F+,

(15.190)

where t/^'(r^ r) - C/^(r, T -r). Transforming from reverse time r to conventional time t = T — r in equation (15.190), we obtain expression (15.188). Note that with the convolution theorem, one can show that, in the frequency domain, the formula for backward extrapolation in reverse time takes the form

U^*{V',UJ) = jj^ G^{V'\V;LO)~U*{V,LO)

-u*{v,uj)—G^{v'\v-uj)

ds,

^' ev^

(15.191)

Thus, we draw the important conclusion that we can select upgoing waves with the aid of the integral Kirchhoff formula by making the transition from direct time t to reverse time r (or by replacing the observed scattered field by its complex conjugate in the frequency domain). That is why this transformation is referred to as reverse-time wave-equation migration. Note that the same idea can be applied to electromagnetic migration as well, ^ shown in Chapter 11.

510

Integral representations in wavefield inversion

15.4-3 Rayleigh integral It can be demonstrated that the boundary values of the wavefield and its normal derivative on surface S are not independent (Zhdanov, 1988). Indeed, if we substitute any two continuous functions, / (r, ^) and ^ (r, ^), instead of the wavefield U ir^t) and its normal derivative dU (r, t) jdn in expression (15.188), we will still generate some wavefield in the lower half-space. This wavefield can be continuously extrapolated back to the surface of observation S. However, this field and its normal derivative may not be equal to the original functions / (r, ^) and g (r^t) on 5. In order for formula (15.188) to reproduce the correct wavefield values in the lower half-space and on the surface S^ it is necessary and sufficient that integral (15.188) be identically equal to zero at the points of the upper half-space, from which we have at once

J Js J-c

G^{r",t'\r,t)—U{r,t)dtds an

^ILL

U{y,t)—G'"{Y",t'\T,t)dtds, Y" i 1/+, (15.192) sJ-oo on where r" is the position vector of a point lying in the upper half-space. This result follows immediately from the Kirchhoff integral formula (13.196). According to formulae (15.189), (15.191), and (13.197), similar relationships hold for the wavefield and its normal derivative in the frequency domain: 11 G'"'

{T"\Y]UJ)

^u{Y,uj)ds

= II

U(Y,UJ)^G''*

(Y"\Y]uj)ds,

Y"

iV^,

(15.193)

U\Y,UJ)^G'"

{Y"\Y]uj)ds,

Y"

^ V^.

(15.194)

and j l G'"{Y"\Y',uj)^u\Y,uj)ds=

jj

The solution of the integral equations represented by the relationships (15.192), (15.193), and (15.194) with respect to functions du{Y^t)/dn or du{Y,uj)/dn proves to be a difficult problem in the case of an arbitrary surface S. However, if the surface 5 is a horizontal plane, these integral formulae can be simplified using a simple geometrical method. Let us position the points r" and r' in the upper and lower half-spaces respectively, as shown in Figure 15-4, i.e. symmetrically with respect to the plane S. Then the following equalities hold true for any point r located on S: G^(r^t'|r,t)

-

G-(r',t'|r,t),

(15.195)

^G-(r^t>,t)

-

-Ag-(r',t'|r,t).

(15.196)

Principles of waveHeld migration

511

V ho) Figure 15-4 The points r' and r" are located in the lower, V'^, and upper half-spaces respectively, symmetrically with respect to the horizontal plane S.

Writing equality (15.192) for horizontal surface S{z = 0), and replacing the Green's function and its normal derivative at the point r" by the Green's function at the symmetrical point r', according to equations (15.195) and (15.196), we obtain the following identity:

ILL ILL

G'"{r',t'\r,t)—U{r,t)dtds

U{r,t)—G'^{r',t'\r,t)dtds, r' e V.

(15.197)

Substituting relation (15.197) into (15.188), we arrive at an important modification of the Kirchhoff integral known as the Rayleigh integral formula (Schneider, 1978; Berkhout, 1980): U'^{r\t') = 2 II

I

U{v,t)-^G'" {Y',t'\v,t)dtds,

(15.198)

where we have noted that d/dn = ~d/dz for the horizontal surface S (where the axis z is directed downward as shown in Figure 15-4, while the normal vector n is directed upward). By replacing conventional time t by reverse time r = T — t (where T is time of observation of the wavefield on the surface of the earth) and introducing an auxiliary

512

Integral representations in wavefield inversion

wavefield f/(r,r) = U{Y^T — r ) , we can write the Rayleigh integral (15.198) for this field as follows: f7(r', T') = 2 ff

r

[7(r, ^ ) ^ G ^ (r', / | r , r ) drds.

(15.199)

Thus, at points in the lower half-space V~^{z > 0), the Rayleigh integral describes an upgoing wavefield which assumes the given values U{r,t) on the observation plane 5, i.e. it may serve as an analytical tool for utilization of the main element of migration - backward extrapolation (continuation) of the observed scattered wavefield. In a similar way we can obtain the Rayleigh integral formula in the frequency domain: ^^(r',a;) = 2 f f uir^uj)—^^* {r'\r;iu)ds. (15.200) For reverse time backward extrapolation in the frequency domain, equation (15.199) can be written as u'^*{r\uj)=2

ffu*{T,u)-^G''{r'\v;uj)ds.

(15.201)

In accordance with the basic principles of migration imaging described above, we produce the image of the refiecting boundaries by plotting the distribution of the migrated wavefield in the lower half-space at time t' = 0. For example, in a homogeneous medium with the constant wavespeed c, the corresponding imaging condition takes the form t/"^(r', 0) = 2 / / +00

•ILL

/ " C/(r, 0 ^ 5 " (r', 0|r, t) dtds

Q

[47r|r'-r| +°°

^

6[t'-t+''' \

'"'

c/2

H{2\Y'~v\/c-ty

dtds t'=0

dtds,

(15.202)

where we have used the half-wavespeed c/2 according to the principle of migration imaging discussed in the previous section. Relation (15.202) provides a theoretical basis for migration transformation in a homogeneous medium, but it is of little use for the construction of practical migration algorithms because of the presence of a singularity in the integrand. However, we can transform formula (15.202) to make it more suitable for practical implementation. Indeed, taking into account the following identity for the Green's function,

Principles of waveReld migration

513

and bringing the derivative with respect to the parameter z' outside the integral, we can rewrite formulae (15.198) - (15.199) as follows: t/-^(r', t') = - 2 ^ / / / U'^{Y\ T') = -2-^

ff

r

C/(r, t)G^ (r', t'|r, t) dtds,

(15.203)

f7(r, T ) G " (r', r ' | r , r ) drds,

(15.204)

The migration formulae in the frequency domain (15.200) - (15.201) can be transformed in a similar way: u'^ir^u)

= - 2 ^ Ij

tx^*(r',a;) = - 2 ^

u{Y,uj)G'"*{Y'\Y',uj)ds.

f f U%Y,U)G''{Y'\Y;U)

(15.205)

ds.

(15.206)

- t)dtds.

(15.207)

In particular, imaging conditions (15.202) are cast in the form [/-(r', 0) = - l - - ^ J 1^ ^

£ ^ U{r, t)6{2\v'-v\/c

The inner integral in relation (15.207) describes the convolution between the observed wavefield t/(r, t) and a delta-function. Therefore, this integral is equal to the values of the field at time t = 2|r'—r|/c. Thus, we can cast the last formula in the form

Expression (15.208) represents a simple algorithm of diffraction transformation, based on summation of the time section along the hodographs of the scattered wavefield (Timoshin, 1978). Another method of practical realization of the Kirchhoff type reverse-time migration is based on the Rayleigh integral formula in the frequency domain (15.200). Applying an inverse Fourier transform to both sides of the Rayleigh formula, we obtain for f = 0 [/^(r',0) =2

f

f I f

= 2 f f f

U(Y,UJ)^G'"''

{Y'\Y]uj)eyi^{-iujt')duj\

u{Y,Uj)^G'^^Y'\Y]Uj)dLUds.

ds

(15.209)

In a homogeneous medium this formula takes the form exp(—z2u;|r'—r|/c)

^"•<''°'=-2^//jr»<''"'

— I ^^^^' 7 7 ' \ ^''^' I duds,

(15.210)

where again we use the half-wavespeed c/2 according to the principle of migration imaging discussed in the previous section.

514

Integral representations in wavefield inversion

15.4-4 Migration in the spectral domain (StolVs method) The backward extrapolation (continuation) of the wavefield is especially easy to reaUze in the spectral domain (Stolt, 1978; Gazdag, 1978; Chun and Jacewitz, 1981). This approach is based on representation of the migration field as a sum of plane harmonic waves characterized by the most elementary law of propagation: U'^{Y\

^) ^ -^

11

/ ^"^(^a;, ky, z', uj) exp \i(kx^ + kyy' - o;^)] dujdk^dky, (15.211)

where u^{kx,ky,z,uj)

= //

u^{x,y,z,uj)exp[—i{kxX-\-kyy)]dxdy.

(15.212)

Substituting expression (15.211) into the wave equation, we obtain a simple second-order differential equation:

dz^

u^{kx,ky,z,uj)

^f^^-kl-

k^] u^{kx^ky,z,uj)

= 0,

(15.213)

which determines the variation of the amplitudes and phases of the plane harmonic waves propagating along the z axis (which is directed vertically downward). The general solution of equation (15.213) can be written as follows: u^ikx, ky, z, oj) = A{kx, ky, cj) exp{ikzz) + B{kx, ky, ou) exp{—ikzz),

(15.214)

where k, = ^a;Vc2 - /c2 -

kl,

and the unknown coefficients A{kxjky,uj) and B{kx^ky,uj) are determined from the boundary condition and the condition at infinity z —> -hoc (with the positive z directed downward). The boundary condition requires that the migration wavefield found in the lower half-space from relations (15.214) and (15.211) should be equal to the observed scattered field on the plane z = 0 , i.e. A{kx,ky,uj) + B{kx,ky,uj)

= u{kx,ky,0,uj),

(15.215)

where u{kx, ky,0,uj) is the space-temporal spectrum of the scattered wavefield measured on the surface of observation. However, condition (15.215) is insufficient for unique determination of the coefficients A{kx,ky,(jij) and B{kx,ky,cj). Following Claerbout (1985), the additional condition may be introduced by specifying that the migration field should be a combination of only upgoing waves. For a real-valued and

Principles of waveReld migration

515

positive square root Ju'^/c^ — k^ — k"^ we choose parameter kz according to formula (15.123), i.e. fc, = sign uj^uyc^

-kl-kl

for \uj/c\ > yJk[Vk^.

(15.216)

With this choice of A:^, the following expressions, u^{kx, ky,0,uj) = A{kx,ky,Lj)exp[i

{k^x + kyy -\- k^z — cut)],

(15.217)

u~{kx^ky^O,(jj) — B{kx^ky^uj)exp [i {kxX + kyy — kzZ — out)]

(15.218)

represent, for a; > 0, the downgoing, u^^ and upgoing, u~^ waves respectively.^ From this it follows that, in accordance with the principles of migration transformation, only the second term should be retained on the right-hand side of equation (15.214), when cj > 0, i.e. coefficient A{kx^ ky, UJ) should be made equal to zero. Then constant B{kx, ky.uj) is found from relation (15.215), B{kx,ky,uj)

==

u{kx,ky,0,uj),

from which we deduce for \uj/c\ > yjk\ -h k^j and a; > 0 ii!'\h,,,ky.^.^)

= u{kx, ky, 0,cj) exp{-iJuj'^/c^

- kl - klz).

(15.219)

Conversely, for a; < 0 the terms on the right-hand side of (15.214) will correspond to upgoing and downgoing waves respectively, so that for |CJ/C| > ^Jk!^ 4- k'l, u; < 0 we have u^'^kx.ky.z.uo)

= u{kx,ky,0,uj)exp{iJuj^/c^

- ^I - k^z).

(15.220)

Combining formulae (15.219) and (15.220) into one expression, we find u^{kx^ ky^ z^uj) — u[kx, ky^ 0,cc;) exp {—% kzz)

(15.221)

where kz is given by expression (15.216). Thus, we conclude that formula (15.221) determines the spectrum of the regular upgoing wavefield in the lower half-space equal to the observed scattered field u on the surface of observation. In other words, it realizes the required principles of ^As usual, we exclude evanescent waves from consideration (arising when the squre root is imaginary, i.e. kz = i\ kl + ^^ — uj'^/c^ for \UJIC\ < Jk'^ + /c^) to simplify the discussion. However, the final migration formula will be the same even if we take evanescent waves into account (Zhdanov et al., 1988).

516

Integral representations in wavefield inversion

wavefield migration. Substituting expression (15.221) into (15.211), we arrive at the time domain representation of the migration field: f/^(r',t)

-I

p r-\-oo

= —- / / STT*^ J J-OO

r+oo

/

u{kx, ky, 0, uj) exp [i{kxx' + kyy' — kzz' — uot)] dwdkxdky.

7-00

(15.222) In the spectral domain this formula can be written as U'^{Y', Lu) = —^ / /

u{kx, ky, 0, uj) exp [i{kxx' -f kyy' - k^z')] dk^dky.

(15.223)

We can produce an image of the medium under examination by reconstructing the migration field in the lower half-space at the initial time t = 0: ix-(x',2/',/,0) = ^

///

/ ^(^^' ^2/' 0'^) ^xp \-i sign u)yjAu'^/c^ - k^ - k^ z'j • • exp [i{kxx' + kyy')] dkxdkydw,

(15.224)

where we use the half-wavespeed c/2 according to the basic principles of migration imaging. Expression (15.224) represents Stolt's (1978) Fourier-based migration formula. 15.4-5 Equivalence of the spectral and integral migration algorithms In the previous sections we analyzed independently two different approaches to wavefield migration: 1) based on the Kirchhoff integral formula, and 2) Stolt's Fourierbased migration. It is important to understand that these two approaches are equivalent. Indeed, it can be demonstrated that, in z > 0, exp[-ik,z]

= -2—g^\kx,ky,z',uj), 'dz

(15.225)

where n r+oo n-\-oo

g'^{kx,ky,z]uj)

= //

G'^{x,y,z]Lu)exp[-i{kxX-\-kyy)]dxdy,

and G^(x, y, z; u) is the Green's function of the Helmholtz equation G^{x,y,z;u;)

= G^{r;u;) = - ^ exp (iJ-^] , 47r |r| \ cJ

(15.226)

Principles ofwavefield migration

517

where r is the position vector of the observation point (x, y^ z). Thus, in the spectral approach, wavefield extrapolation from one horizontal level to another is implemented by multiplying its spectrum by the complex conjugate of the vertical derivative of the spatial spectrum of the Green's function for the Helmholtz equation. Let us prove the equality (15.225). It can be shown that the function g'^{kx, ky, z; uj) satisfies the 1-D Helmholtz equation d^

.2'

g'"{k^,ky,z;u)) = -6{z).

(15.227)

Thus, g^ is the 1-D Green's function for the medium with the constant velocity CQ = 1 (see formula (13.79)): g'"{k^,ky,z;u)

= - ^ e x p ( i f c , \z\).

(15.228)

Therefore, for z > 0 —g'"\k:,,ky,z\Lu)

=

--exp{-ik^z),

which was to be proved. Substituting equation (15.225) into (15.223), we arrive at the integral formula

= - 2 — / /

u{k:r,ky,0,iu)—g'^*{kx,ky,z]Lj)exp[i{k:r:X

-\-kyy]dk:,dky.

It is well-known that the Fourier transform of the product of spectra of two functions is equal to the convolution of these functions; therefore we obtain ^i"^(r^cJ) = - 2 ^ IJu{Y,uo)G'"*{Y'\Y]u)ds,

(15.229)

which is exactly the Rayleigh integral formula (15.205) in the frequency domain. Thus we have demonstrated that the spectral extrapolation based on Stolt's method (formulae (15.222) and (15.223)) is equivalent to the reverse-time wave equation migration using the Rayleigh integral formulae (15.203) and (15.205). 15.4.6 Inversion versus migration In the section of this chapter describing the Kirchhoff inversion method and general non-linear inversion techniques, we have demonstrated that the calculation of the Kirchhoff adjoint operator (15.142) and the adjoint Frechet derivative operator

518

Integral representations in wavefield inversion

(15.180) requires back propagation of the residual field. This back propagation transformation can actually be performed by migration algorithms. Thus, conventional migration is equivalent to the calculation of the gradient direction (adjoint Frechet derivative) of inversion, except that it does not include convolution with the incident field (similar to one shown in formulae (15.142) and (15.180)). Actually, we can say that conventional migration and inversion algorithms employ diff'erent imaging conditions. This diff"erence can be easily explained by the fact that most of the inversion methods outlined above reconstruct the material property (slowness) distribution, while the migration method was designed to image the reflecting boundaries. In this sense, the closest similarities are between migration and the Bleistein and Kirchhoff inversions, because both of these methods recover the reflectivity function which helps to locate the reflection boundaries. For example, it was demonstrated by Bleistein et al., (2001) that inversion formula (15.129) can be reduced to Stolt's migration formula (15.224). Another important diff'erence between inversion and migration is that the general inversion algorithms are iterative. Without iterations, an inversion algorithm will not be able to reconstruct both geometrical and material properties of the medium. At the same time, iterative inversion, being a solution of an ill-posed problem, requires regularization. By contrast, migration is a stable procedure because it is equivalent to one forward modeling for computing the back-propagated field. However, we have demonstrated in this chapter that any inversion algorithm can be treated as an iterative migration of the residual fields obtained on each iteration. In this way, the difference between these two approaches to the interpretation of seismic data disappears. In conclusion, I should note that modern seismic migration algorithms are based mostly on finite-difference methods of the solution of back propagation problems. However, a description of these methods lies beyond the scope of this book. 15.5

Elastic field inversion

The problem of elastic field inversion is much more complicated than acoustic or vector wavefield inversion, considered in the previous sections of the book. However, the fundamental principles of elastic inversion resemble those discussed above for more simple models of seismic waves. I will present in this section a brief overview of the.basic ideas underlining the elastic field inversion. 15.5.1 Formulation of the elastic field inverse problem Consider a 3-D elastic medium. The propagation of elastic waves in the frequency domain can be described by the Lame equation (13.34) c^VV • u(r,a;) - c^V x V x M{Y,UJ)-{-iJ\i{v,uj)

= --f^

(r,cj),

(15.230)

where u(r,a;) is the elastic field, f^ ( r , ^ ) is the strength of an external source of energy located within some bounded domain Q, p is the density, and Cp (r) and Cs (r)

Elastic field inversion

519

are the Lame velocities. We also assume that the elastic field satisfies the radiation conditions (13.202) - (13.204) at infinity. We will consider a general nonlinear elastic inverse problem: d^ = AL ( m ) ,

(15.231)

m(r)=(|j^j)

(15.232)

where

is a column vector function describing the distribution of the model parameters (square velocities, Cp{r) and Cg(r)), and d^, represents the elastic field observed on the surface S. The goal is to find the unknown parameters Cp{r) and Cg(r). Operator Ai denotes the nonlinear forward modeling operator, given by the Lame equation (15.230) and radiation conditions (13.202) - (13.204). This operator can be calculated, for example, from the general integral representation of the elastic field in the frequency domain (13.97), which we write here in the form U{T\UJ) = III

- r ( r , a ; ) • G^{Y'\Y,u)dv

= A^. ( m ) ,

(15.233)

where G^'(r|r',a;) is the Green's tensor for the Lame equation in the frequency domain, which satisfies c^VV•G^^(r|r^a;)-c^V x V X&\Y\V',UO)^UJ'^&\Y\Y',uj)

= - I ( 5 ( r - r ' ) . (15.234)

Following our standard approach to regularized solution of the inverse problem, we substitute for the solution of the inverse problem (15.231) a minimization of the corresponding parametric functional with, for example, a minimum norm stabilizer: P^ (m) - (A^(m) - d^, Az.(m) - d^)^^ -\-a {W (m - mapr), W {m - mapr))^^ = min,

(15.235)

where W is the corresponding weighting operator, and the inner product (.., . ) ^ and (.., .)^^ in the Hilbert spaces D^ and M^ of data and model parameters are determined by formulae similar to (15.11) and (15.12): (u, v)^^ = R e / / / u(r,a;) - ^r* {Y,UJ) dsdw] u , v G Z^L, ^ e 1^, (15.236) Jn J Js where ds denotes a differential element of the surface S^ and integration is conducted over a variable ?, and

(mi, m2)^^ = / / / m f (r)m2(r)c/^; m i , m 2 G ML,

(15.237)

520

Integral representations in wavefield inversion

where the superscript T denotes transposition, and rrii, m^ are the column vectors of the type, shown in formula (15.232). In particular, taking into account (15.232), we can write the last expression in explicit form: (mi, m s ) ^ ^ = j j j

[cl^{r)c%{T) + CI^{Y)CI^{T)\

dv.

The minimization problem (15.235) can be solved by the conjugate gradient or reweighted conjugate gradient methods introduced in Chapter 5. Let us describe, for example, the algorithm based on the regularized conjugate gradient method (5.92), which in the case of the elastic inverse problem takes the form: ^Ln = AL(m„) - dL,

(a)

m.+i=m.-fc^I^,

{d)

(15-238)

k^ is determined from the minimization problem P« (m,+i) = P " ( m , - fc-I-) = min,

n = 0,1, 2,... (e)

where, as usual, we denote by nin the distribution of the model parameters at the n-th iteration. As is usual in the inversion method, the practical implementation of this algorithm requires computing the Frechet derivative F^.^ of the corresponding forward modeling operator A^ on each n-th iteration. 15.5.2 Frechet derivative for the elastic forward modeling operator Let us assume that the Lame velocities are known and fixed everywhere in space with the exception of some local domain D. We can find the equations for the Frechet derivative by applying the variational operator to both sides of the corresponding Lame equation (15.230) c^VV-5u(r,a;) - c^V X V X 5u(r,u) H- J^6\I{Y,UJ)

- -f^(r,cj),

(15.239)

where f^(r,a;) - 6cl (r) V V • u(r,a;) - dc] (r) V x V x u(r,a;),

(15.240)

Elastic field inversion

521

6Cp (r) and Sc^ (r) are the square velocity variations, which are unequal to zero only within domain D, and Su is the corresponding elastic field variation. For the sake of simplicity, we assume that there are no density variations in equation (15.239), i.e. 5p = 0. We also take into account that the perturbation of an external source is equal to zero: 6f^ = 0. Applying integral representation (15.233) to the elastic field variation satisfying equation (15.239), we obtain <5u(r',a;) = JJJ

f^(r,a;) • G''{v'\r,uj)dv

= A^ {4,0',) .

Therefore, the Prechet differential of the forward modeling elastic field operator is given by the following expression: FL(m,5m) = ^u [6cl (r) V V • u(r,a;) - Sc^ (r) V x V x u(r,c^))] •

IIL

Srn^ (r) Du(r, uj) • G^(r'|r, Lj)dv

G^{r'\r,uj)dv. (15.241)

'D

where 6m is the 2 x 1 matrix given by the formula ,5m=,5m(r)=(^|j^j),

(15.242)

and the symbol D is the 2 x 1 matrix differential operator, formed of the gradient divergence V V - and curlcurl V x V x differential operators, such that ° = ( ^ V XVx ) •

(1^-243)

Thus the symbol D u can be regarded as the 2 x 1 matrix of "ordinary" vectors V V • u and —V X V X u (i.e. 3-D vectors in real physical space):

\ - V XV X u Therefore, the first product 6m^ (r) Du(r,ct;) in (15.241) above is a (matrix) multiplication of a row vector and a column vector, while the second product -G^ is the tensorial inner product of this ordinary vector with the tensor G^. In particular, it is easy to show that [6ra^ (r) Du(r,a;)] • G^ = ^m^ (r) [Du(r,a;) • G^

522

Integral representations in >vavefield inversion

if we assume that

G^ = Note that vector u in expression (15.241) represents the elastic field for the given velocities Cp and Cg, and the Green's tensor is calculated for the same Lame velocities. We will use below the following simplified notation for the Frechet differential FL{m,6m) = F2^((5m). Thus, if we know the elastic field for the given velocity model, we can find the corresponding Frechet diff'erential using the integral representation (15.241). This formula can be treated as the Born approximation for the elastic field. Using this formula, one can construct a family of linear and nonlinear approximations of the elastic field similar to the corresponding approximations for electromagnetic, acoustic, and vector wavefields. However, this topic goes beyond the scope of our discussion. 15.5.3

Adjoint Frechet derivative operator and hack-propagating elastic field

Now we can address the problem of calculating the adjoint Frechet derivative operator. Note that in the RCG algorithm (15.238) expression F^^ (RLU) denotes the result of an application of the adjoint Frechet derivative operator to the corresponding residual field R^n =" A£,(inn) — AL on the n-th iteration. The explicit expression for the adjoint Frechet derivative operator F ^ can be found according to the definition: (Fr((5m),v)^^ = ( ^ m , F r ( v ) ) M , .

(15.244)

Using the definitions (15.236) and (15.237) of inner products and the expression (15.241) for the Frechet derivative operator, we can rewrite formula (15.244) as (Fr(^m),v)^^ = Re / / / [ iff Ju J Js U J

IIL

Sm'^{r)Rei

6m^ (r) Du(r, w) • G^(^\v,uj)dv • V* (r, w) dsdu) JD

/ Du*(r,w)

JjG'^*{r\r,u;)-v{T,oj)ds]

= (^m,Fr(v)U.

duoy dv (15.245)

From the last formula we have (5m, F r (v) - Re f Bu%r,uj) • \ [[ G^%r\T, cu) - v {r, u) ds du\ ==0. J ML (15.246)

Elastic field

523

inversion

Equation (15.246) should hold for any 6m, for example for 6m = F r (v) - Re / Du*(r,cj) • | / / G^*(?|r,cj) • v (?,u;) dS du. Jn [J Js Therefore, F r (v) - Re / Du*(r,u;) • [ / / G^*(?|r,a;) • v (?,a;) dS dw JQ

U Js

= 0. M

From the last formula we conclude that the adjoint Frechet derivative operator F^* is given by the formula F ^ * (v) = Re / Du*(r,a;) • | / / G^*(?|r,a;) • v (Y,U) d? du). Jn \_J Js

(15.247)

which can be displayed in the form of 2 x 1 matrix: ^[VV.u*(r,a;)]. F r (v) = Re

//^G^*(?|r,a;)-v(?,c^)dS'

duj

/ J V X V X u*(r,c^)] . [//<,G^*(?|r,a;) • v(?,a.)dS-_ dw

Based on these formulae, we can write Fr,:(RLn) in [ V V • u;(r, w)] • [/J^G^f (r|r,w) • R L „ (r,c^) £&] duj Rt

, (15.248) - ^ [V X V X < ( r , Lo)] • [//^ G^-(?|r, uj) • K^n (?, uj) dsd(jj

where u^(r,u;) and G^'(r|r,cj) are the elastic fields and Green's tensors, computed on the n-th iteration for the Lame velocities, c^^ (r) and c^^ (r) . The asterisk * denotes the operation of complex conjugate. The surface integral term in the last formula can be treated as the complex conjugate elastic field, u ^ (r,a;), u^ (r,o.) = II

G^{r\v,Lo) • R*„ {r,Lo) ds,

(15.249)

generated at the point r by fictitious sources distributed along the surface of observation S, with the surface density determined by the complex conjugate residual field: FL(?,cc')=RL(?,^), (15.250) Therefore, formula (15.248) can be cast as vm* . p

^^-n^f

In ["^^ • < ( r , w)] • < * (r>w) d-^

(15.251)

524

Integral representations in wavefield inversion

Note again that we can drop the symbol Re in the above formulae if we integrate over a symmetrical frequency interval f^ (for example, from — oo to +00) because the imaginary part of the integrand is anti-symmetric in a;,

or in equivalent from pm* (n

\ - f In [ ^ ^ • ^n{r, w)] • < (r.w) dw

\

We can get a better understanding of the physical meaning of the elastic adjoint Prechet operator if we consider the first iteration in the inversion scheme (15.238). Let us assume that the initial distribution of Lame velocities in the model is given by some background parameters Cpb (r) and Cst (r): mo = mo (r) = ( J^^j^j ") .

(15.253)

Then the first approximation of the model parameters, m i ( r ) , is mi = mo + Amo = mo - fco lo, where Amo = - ^ 0 lo and lo = F^o (RLO) - F*o (AL(mo) - dz.).

(15.254)

Application of the elastic forward modeling operator to the initial model with the background Lame velocities results in nothing more than the incident elastic field in the background model, u*(?, cj). The diff'erence between this incident field and the observed total field d/, gives us the scattered elastic field (with a minus sign): R;.o = AL(mo) - d^ = u^(?,a;) - dL{7,uj) = -u^(?,a;), r e S.

(15.255)

Substituting (15.255) into (15.254) and taking into account expression (15.252) for the adjoint Frechet operator, we find Amo = koFl^ (u^)

= ko(^^Z-''^Ml-''7Hf' ,, ) , where, by reciprocity the auxiliary field u^* {r,u) is given by the formula u^* (r,a;) = J J Gf *(r|r, oj) • u ' (?, LJ) dS,

(15.256)

Elastic field inversion

525

and G^(r|r,a;) is the elastic Green's tensor computed for the background velocity model. We may now recall that according to (15.253). the column vector ^mo represents corrections to the velocities Cpb (r) and Cgb (r) respectively, on the first iteration:

Therefore, formula (15.256) can be separated into two expressions Ac%{r) = kJ [ V V • vi»{vM] Jn

• u«* (r,a;)
= ko [ [VV-u^(r,a;)] - u^ {r,uj) du, Jn

(15.257)

and A 4 ( r ) = -/Co / [V X V X u^*(r,a;)] • u^* (r,a;) dw Jn = - ^ 0 / [V X V X u^(r, uj)] • u^ {r,uj) du. (15.258) Jn As in the case of the waveform inversion method, we note that by using the convolution theorem and a reciprocity, we can write the expressions for the inverse Fourier transform U^^ (r,^^ ^^ ^^^ auxiliary field u^* {Y.UJ) as follows: U^ {Y,t') = — / u^* (r,cj) e-'^^^'dw 27r7_oo = // /

W(y,t)-G^(i,t\v,t')dtd'S,

(15.259)

J JS J —oo

where Gj?,i| r,i')=G^(?,-t|

r,-0.

Therefore, the scattered field U ^ (r,t') is the back-propagated field. We can also write the result of an application of the adjoint Frechet operator to the scattered field as the correlation of this back-propagated scattered field and the incident field U*(r,t) acted on by the relevant difiPerential operators: +00

/

_

[ V V • U' (r,i)] • U« (r,i) dt,

(15.260)

-OO

and +00

/

_

[V X V X U ' (r,i)] • U ^ (r,i) dt. -OO

(15.261)

526

Integral representations in wavefield inversion

Thus, the corrections to the velocities Cpt (r) and Cgb (r) on the first iteration can be obtained by correlating the back-propagated scattered elastic field with the derivatives of the incident field U*(r,t). This transformation is similar to wavefield migration described previously. Thus we see that elastic field inversion can be constructed by iterative migration of the residual elastic data. Formulae (15.260) and (15.261) can be further simplified, if we assume that we have constant background velocities Cpb (r) = const and Csb (r) = const. In this case the incident field can be separated into the compressional and shear waves U'(r, t) = U;(r, t) + U:(r, t),

(15.262)

where according to (13.37) and (13.38) V X u ; - 0,

V • u; = V • u\

V X U*, - V X U \

V • U^, = 0.

(15.263) (15.264)

Substituting representation (15.262) into (15.260) and (15.261), and taking into account the properties of the compressional and shear waves (15.263) and (15.264), we arrive at a result that calculation of the anomalous compressional wavespeed requires only the compressional incident wavefield, while the anomalous shear wavespeed is determined solely by the shear incident wavefield: + CXD

_

/

[VV-U;M]-U«Mdi,

and

(15.265)

•CO +00

_

/ [V X V X U ; {r,t)] • U ^ (r,i) dt. (15.266) We have found in Chapter 13, that the compressional and shear waves in a homogeneous domain satisfy the vector wave equations (13.44) and (13.46). Using the vector identity V^U = V V - U - V x V x U •OO

and properties of the compressional and shear waves (15.263) and (15.264), we can express these equations, in a domain free of external forces, in the form

cl,v'v; --gf dV'u: - ^

= cl,vv. u;- u^= o, = -civ

X V X U^,- u ; = . 0,

(15.267) (15.268)

where dots over the vectors U^ and U'^ denotes their differentiation with respect to time.

Elastic field inversion

527

Substituting the last two formulae into equations (15.265) and (15.266), we obtain Ac^(r) = ^

/

Up (r,i) • U « (r,i) dt,

(15.269)

Ad(r) = ^

/

U , (r,i) • U « (r,^) di.

(15.270)

and

Integrating by parts, we can move one time derivative from the incident field to the back-propagated residual field. As a result, we arrive at the following expressions for anomalous velocities ACp^(r) and Ac^5(r) Acl,{r) = - 4

/

U^ (r^^) • U

(r,^) dt,

(15.271)

(r,t) dt

(15.272)

and ka

f^"^

Aclir) = -4

•^

-^

U , (r,t) • U

^s6 J —oo

Thus the anomalous square velocities can be found as the correlations of the time derivatives of the back-propagated residual field and the corresponding incident wavefields - the compressional wavefield for determining ACp^(r) and the shear wavefield for determining Ac^^(r). This result is similar to one obtained previously for acoustic and scalar wavefield inversions.

References to Chapter 15 Arfken, G., and H. J. Weber, 1995, Mathematical methods for physicists. Fourth Edition: Academic Press, San Diego, New York, Boston, London, Tokyo, Toronto, 1028 pp. Aki, K., and P. G. Richards, 2002, Quantitative seismology, 2nd Edition: Columbia University Press, 576 pp. Berkhout, A. J., 1980, Seismic migration - Imaging of acoustic energy by wave field extrapolation: Elsevier, Amsterdam, Oxford, New York, 339 pp. Berkhout, A. J., 1984, Seismic migration - Imaging of acoustic energy by wave field extrapolation, B. Practical aspects: Elsevier, Amsterdam, Oxford, New York, 274 pp. Bleistein, N., 1984, Mathematical methods for wave phenomena: Academic Press Inc., (Harcourt Brace Jovanovich Publishers), New York, 341 pp. Bleistein, N., Cohen, J. K., and J. W., Stockwell, Jr., 2001, Mathematics of multidimensional seismic imaging, migration, and inversion: Springer, New York, Berlin, London, Tokyo, 510 pp. Chew, W. C , 1990, Waves and fields in inhomogeneous media: Van Nostrand Reinhold, New York, 608 pp.

528

Integral representations in wavefield inversion

Chun, Joong, H., and Ch. A. Jacewtiz, 1981, Fundamentals of frequency domain migration: Geophysics, 46, No. 5, 717-733. Claerbout, J. F., 1970, Coarse grid calculations of waves in inhomogeneous media with application to delineation of complicated seismic structure: Geophysics, 35, 3, 407-418. Claerbout, J. F., 1976, Fundamentals ofgeophysical data processing: McGrawHill, New York, 274 pp. Claerbout, J. F, 1985, Imaging the Earth's interior: Blackwell Scientific PubUcations, Oxford, London, Edinburgh, 399 pp. Claerbout, J. F., 1992, Earth soundings analysis. Processing versus inversion: Blackwell Scientific Publications, Oxford, London, Edinburgh, 304 pp. Claerbout, J. F. and S. M. Doherty, 1972, Downward continuation of moveout corrected seismograms: Geophysics, 37, No. 5, 741-768. Devaney, A. J., 1984, Geophysical diffraction tomography: IEEE, Geoscience and Remote Sensing, 22, 3-13. Gardner, G. H. F., Editor, 1985, Migration of seismic data: Geophysical reprint series, No. 4, Society of Exploration Geophysicists, Tulsa, Oklahoma, 462 pp. Gazdag, J., 1978, Wave equation migration with the phase-shift method: Geophysics, 43, No. 7, 1342-1351. Hagedoorn, J. G., 1954, A process of seismic reflection interpretation: Geophysical prospecting, 2, 85-127. Kozlov, Ye. A., 1986, Migration transformations in seismic prospecting (in Russian): Nedra, Moscow, 247 pp. Loewenthal, D., Lu, L., Roberson, R., and J. Sherwood, 1976, The wave equation applied to migration: Geophysical Prospecting, 24, 380-399. Mora, P. R., 1987, Nonlinear two-dimensional elastic inversion of multi-offset seismic data: Geophysics, 52, No. 9, 1211-1228. Press, W. H., Flannery, B. P., Teukolsky, S. A., and W. T. Vettering, 1987, Numerical recipes. The art of scientific computing: Cambridge University Press, Cambridge, 818 pp. Schneider, W. A., 1978, Integral formulation for migration in two and three dimensions: Geophysics, 43, No. 1, 49-76. Stolt, R. M., 1978, Migration by Fourier transform: Geophysics, 4 3 , No. 1, 23-49. Tarantola, A., 1987, Inverse Problem Theory: Elsevier, Amsterdam, Oxford, New York, Tokyo, 613 pp. Timoshin, Yu. V., 1978, Pulse seismic holography (in Russian): Nedra, Moscow, 321 pp. Zhdanov, M. S., 1988, Integral transforms in geophysics: Springer, New York, Berlin, London, Tokyo, 367 pp.

Elastic field inversion

529

Zhdanov, M. S., Matusevich, V. U., and M. A. Prenkel, 1988, Seismic and electromagnetic migration (in Russian): Nauka, Moscow, 376 pp. Zhdanov, M. S., and E. Tartaras, 2002, Inversion of multi-transmitter 3-D electromagnetic data based on the localized quasi-linear approximation: Geophys. J. Int., 148, No 3. Zhou, C , and L. Liu, 2000, Radar-difTraction tomography using modified quasi-linear approximation: IEEE Transactions on Geoscience and Remote Sensing, 38, No. 1, 404-415.