Troyan V., Kiselev Y. Statistical Methods of Geophysical Data Processing

Подождите немного. Документ загружается.

Tomography methods of recovering the image of medium 353

Using (11.36) we can obtain a coordinate form of the Green function

(x, t; x

, t

) =

cos

(t − t

)ϕ

(x)ϕ

∗

) ,

(x, t; x

, t

) =

sin

√

(t − t

)

√

(x)ϕ

∗

) .

If we have the system of eigenfunctions for the operator H, then we can write

the particular representation for the Green functions of the wave equation.

Let us consider, for example, the wave equation without the source (homogeneous

wave equation) for the case of the uniform space:

∆ϕ −

∂

∂t

ϕ = 0.

The eigenfunctions of the operator, which in this case is equal to (−∆) (x ∈ R

are determined as

(2π)

3/2

ikx

, x ∈ R

Here the normalized coeﬃcient corresponds to the condition (ϕ

, ϕ

) = δ(k −k

Thus, the coordinate form of

reads as

(x, t; x

, t

) =

(2π)

∞

−∞

sin c|k|(t − t

)

c|k|

ik(x−x

)

4πc|x − x

{δ(|x − x

| + c(t − t

)) − δ(|x − x

| − c(t − t

))}.

From this representation it is followed that under condition t > t

just only the

second part of the Green function “works”, i.e.

(x, t; x

, t

) =

4πc



−

δ(|x − x

| − c(t − t

))

|x − x



∆

= G

(−)

(11.37)

for t > t

, that corresponds the signal propagation in the direct time. In the case

of the inverse time t < t

the Green function is equal to

(x, t; x

, t

) =

4πc



−

δ(|x − x

| + c(t − t

))

|x − x



∆

= G

(+)

for t < t

11.7.2 Green function for “Poisson equation”

Let us write the Poisson equation in the form

Hϕ = s ,

where H > 0, s = s(x) is the source function. The Green function in the operator

form can be written as G = H

−1

, and in the coordinate form it reads as

G(x, x

) =

6=0

−1

(x)ϕ

∗

) . (11.38)

354 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Here λ

are eigenvalues of the operator H. In the expression (11.38) the summation

is implemented only if λ

6= 0, that is equivalent to the representation of the solution

in the form

ϕ = H

−1

s + kϕ

where ϕ

is an eigenfunction of the operator H that corresponds to λ

= 0, that

leads to the equality

Hϕ = s = Hϕ + kHϕ

for an arbitrary k. If the operator H is equal to −∆ (x ∈ R

) (we consider the

Poisson equation in the inﬁnite space), then the Green function in the coordinate

form can be written as

G(x, x

) =

(2π)

exp[ik(x − x

)]

2π

|x − x

∞

sin k|x − x

dk =

4π

|x − x

. (11.39)

The solution of the Poisson equation can be written down as follows

ϕ(x) =

4π

s(x

) dx

|x − x

+ kϕ

(x),

where ϕ

: ∆ϕ

= 0; k is an arbitrary number.

11.7.3 Green function for Lame equation in uniform isotropic

inﬁnite medium

The vector Lame equation Lϕ

ϕ = s has an unique solution, which can be represented

with the help of the Green tensor

G such that its integral kernel is the tensor function

ˆg(x, t; x

, t

) which satisﬁes the equation

Lˆg = δ(x − x

) δ(t − t

)I.

In the inﬁnite uniform space the retarded Green function g

−

(x, x

; t, t

) satisﬁes the

condition g

−

≡ 0 where t < t

. To construct the fundamental solution of the Lame

equation in the explicit form, we use, taking into account the Helmhotz theorem,

the representation of the vector ﬁeld by scalar and vector potentials (4.20):

ϕϕ = −∇∇Φ +

∇

∇ × A ,

s = −∇∇Φ

+ ∇∇ × A

which satisfy the scalar and vector wave equation respectively

Φ − v

∆Φ =

, (11.40)

A − v

∆A =

. (11.41)

Tomography methods of recovering the image of medium 355

The representation of the vector s as a sum of the vectors of potential and solenoidal

ﬁelds allow us to introduce the vector ﬁeld H such as

−∇

H = s = −∇Φ

+ ∇ × A

Taking into account the identity law −∇

H = −∇(∇ · H) + ∇ × (

∇

∇ × H), it is

possible to ﬁnd that

= ∇

∇

∇ · H , (11.42)

= ∇∇ × H . (11.43)

Let us write a solution for the vector Poisson equation with the help of the

solution for the scalar Poisson equation for the case of the vector point source,

located in the origin of coordinates and oriented in the direction of the ﬁrst ort

of the coordinates, i.e. s(x, t) = e

s(t)δ(x). The solution of the vector Poisson

equation ∇

H = −s, using Green function for the scalar wave equation (11.39)

reads as

H(x, t) =

4π

s(t)δ(x

)

|x − x

4π

s(t)

|x|

. (11.44)

Using the relations for the scalar (11.42), vector (11.43) potential and the represen-

tation for H in the form (11.44), we obtain

= ∇·



4π

s(t)

|x|



s(t)

4π

∂

|x|

−1

, (11.45)

∇

∇ ×



4π

s(t)

|x|



s(t)

4π



∂

|x|

−1

∂

|x|

−1



. (11.46)

Substituting the source functions in the forms (11.45) and (11.46), we rewrite the

wave equations (11.40) and (11.41):

Φ − v

∆Φ =

s(t)

4πρ

∂

|x|

−1

, (11.47)

A − v

∆A =

s(t)

4πρ



∂

|x|

−1

−∂

|x|

−1



. (11.48)

We write retarded Green functions for the wave equation (11.47), (11.48), using

the Green function, obtained in the case of L = (∆ −

∂

) (see (11.37)):

(x, x

; t −t

) =

4πv

δ(|x − x

| − v

(t − t

))

|x − x

4πv

δ((t − t

) − |x − x

|/v

)

|x − x

, (11.49)

356 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

(x, x

; t −t

) = I

4πv

δ((t − t

) − |x − x

|/v

)

|x − x

. (11.50)

To represent the solutions of the wave equations (11.47) and (11.48) in the form

Φ(x, t) =

(x, x

; t − t

)



s(t

)

4πρ

∂

−1



, (11.51)

A(x, t) =

(x, x

; t −t

)

s(t

)

4πρ



∂

|x|

−1

−∂

|x|

−1



. (11.52)

Then, using the expressions for the Green function (11.49) and (11.50) the solutions

(11.51), (11.52) we write down as follows

Φ(x, t) =

(4π)

ρv

δ((t − t

) − |x − x

|/v

)

|x − x

s(t

) ∂

−1

, (11.53)

A(x, t) =

(4π)

ρv

δ((t − t

) − |x − x

|/v

)

|x − x

× s(t

)



∂

−1

−∂

−1



. (11.54)

After the integration over t

the solution (11.53) and after the change of variables

|x − x

|/v

= τ, we obtain the representation for the scalar potential

Φ(x, t) =

(4π)

ρv

∞

dτ

s(t − τ)

|x−x

|=v

∂

−1

dσ . (11.55)

To bear in mind the calculation the surface integral from the right hand side of

(11.55), we consider the following integral

∂

|x−x

|=r

dσ

, (11.56)

where r = v

τ, dσ = 2πr

sin θ dθ, x

= x − r,

, x

) = (x, x) −2(x, r) + (r, r)

(Fig. 11.2) or x

= x

− 2xr cos θ + r

. To write down an exact diﬀerential of the

last expression and to obtain

= xr sin θ dθ ,

and from it follows

dσ =

2πr|x

|dx

To calculate the integral (11.56):

Tomography methods of recovering the image of medium 357

Fig. 11.2 Illustration of the calculation of the integral from right hand side of the equality (11.55).

∂

|x−x

|=r

dσ

= 2πr∂

|x|

−1

|x|+r

|x|−r



4πr

∂

|x|, if r < |x|,

0, if r > |x|.

(11.57)

Using representation (11.57) for the surface integral, we rewrite the expression

(11.55) for the potential Φ(x, t):

Φ(x, t) =

4πρ

(∂

|x|

−1

)

|x|/v

τs(t − τ)dτ. (11.58)

Taking into account the representation (11.57) and after change of variable |x −

|/v

= τ in the solution (11.54) we can calculate the vector potential:

A =

4πρ

|x|/v

τs(t − τ)dτ



∂

|x|

−1

−∂

|x|

−1



. (11.59)

Using obtained solutions for scalar and vector potentials (11.58) and (11.59) we

restore the potential and solenoidal components of the ﬁelds ϕϕ

and ϕ

= −

∇

∇Φ = −

4πρ



∂

|x|

−1

∂

|x|

−1

∂

|x|

−1



|x|/v

dτs(t − τ)τ

−

|x|s(t − |x|/v

)

4π ρv



∂

|x|

−1

∂

|x|

∂

|x|

−1

∂

|x|

∂

|x|

−1

∂

|x|



= ∇

∇

∇ × A =

4πρ



−∂

|x|

−1

− ∂

|x|

−1

∂

|x|

−1

∂

|x|

−1



|x|/v

τs(t − τ)dτ

|x|s(t − |x|/v

)

4π ρv



−∂

|x|

−1

∂

|x| − ∂

|x|

−1

∂

|x|

∂

|x|

−1

∂

|x|

∂

|x|

−1

∂

|x|



358 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Introducing the notations γ

= ∂

|x| and under the assumption that the source

direction is e

, i.e. after change 1 on j, we obtain, after some mathematical manip-

ulations,

ϕϕ = ϕ

+ ϕϕ

4π ρv

|x|



t −

|x|



+ (

I −

P )

4π ρv

|x|



t −

|x|



+ (3

P −

4π ρ|x|

|x|/v

τs(t − τ)dτ.

The retarded Green function after the integration over τ can be represented in

the form

g(x, x

; t −t

4πρ



δ(t − t

− r/v

)

(

I −

P )

δ(t − t

− r/v

)

+(3

P −

(t − t

)[θ(t − t

− r/v

) − θ(t − t

− r/v

)]



where θ(t) is the rump function,

θ(t) =



1, if t ≥ 0 ,

0, if t < 0 ;

r = x − x

, r = |x − x

|; e

|r|

, P = (e

, e

11.7.4 Green function for diﬀusion equation

Let us to write a “diﬀusion equation” in the form



∂

∂t

+ H



ϕ = 0 , (11.60)

where H > 0. The operator H acts on space variables only. To ﬁnd the Green func-

tion of the equation (11.60) we represent its solution as expansion on eigenfunctions

of the operator H:

ϕ(x, t) =

∞

n=0

(t)ϕ

(x). (11.61)

Substituting the expansion (11.61) to the equation (11.60), we obtain the equation

for the coeﬃcients A

(t):

∞

n=0

(

+ A

)ϕ

= 0 ,

i. e.

(t) = A

(0) e

−λ

. (11.62)

Tomography methods of recovering the image of medium 359

In the expansion (11.61) the coeﬃcient A

(t) is equal to projection of the func-

tion ϕ(x, t) on the ort of the proper basis of the operator H ϕ

, i.e.

(t) = (ϕ

, ϕ) .

Denoting ϕ(x, t = 0) using ϕ

, we write down:

(0) = (ϕ

, ϕ

) . (11.63)

Substituting the expression (11.63) into the representation (11.62) and, further,

substituting the representation (11.62) into the expansion (11.61), we ﬁnd a link

between ϕ(x, t) and ϕ(x, t = 0):

ϕ(x, t) =

∞

n=0

−λ

(x)(ϕ

, ϕ

) . (11.64)

From the formula (11.64) we can see, that the Green function in the operator form

can be written as

G = e

−

, in the coordinate representation we have

G(x, x

; t, t

) =

∞

n=0

−λ

(t−t

)

(x)ϕ

∗

) .

Let the operator is H = −a

∇

(x ∈ R

), then the Green function for the

diﬀusion equation (t > t

) reads as

G(x, t; x

, t

) =

a(2

π(t − t

))

exp



−

(x − x

)

(t − t

)



11.7.5 Green function for operator equation of the second genus

Let us write the equation under the consideration in the form

ϕ + Hϕ = s , (11.65)

where H > 0. If we represent the equation (11.65) in the form (I + H)ϕ = s, then

we write the coordinate representation of the Green operator G = (I + H)

−1

follows

G(x, x

) =

∞

n=0

(1 + λ

)

−1

(x)ϕ

∗

) .

As an example to ﬁnd the Green function of the equation

(−∇

+ m)ϕ = s.

The coordinate representation of the Green function can be written as follows

G(x, x

) =

8π

exp[ik(x − x

)]

+ m

4π

i|x − x

∞

−∞

exp[ik|x − x

+ m

k dk =

4π

exp[−m|x − x

|/2]

|x − x

This expression is known as Yukawa’s potential.

360 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

11.8 Examples of the Recovery of the Local Inhomogeneity

Parameters by the Diﬀraction Tomography Method

Let us consider the results of the numerical simulation on the recovery of the param-

eters of local inhomogeneities by the diﬀraction tomography method. The resolution

of the diﬀraction tomography is demonstrated by the recovery results of a bulk of

inhomogeneities with sizes much smaller of the wavelength, which are placed in the

3-D uniform reference medium (Troyan and Ryzhikov, 1994, 1995). The results of

the recovery of the parameters of local inhomogeneities of a complex shape and

the sizes comparable with the wavelength (Kiselev and Troyan, 1998) allows us to

estimate the errors sequent from the Born approximation, which is used here to

linearize the inverse problem solution.

11.8.1 An estimation of the resolution

Let us consider the numerical experiments (Troyan and Ryzhikov, 1994, 1995), for

the of the resolution of the diﬀraction tomography method under the simply ordered

array of three-component seismometers. The reference medium is assumed to be

inﬁnite and uniform. The incident wave ﬁeld is a plane p-wave with a deﬁnite time

shape and a given normal vector to the wave front. The resolution is tested on a

regular set of the point diﬀractors. The synthetic seismograms (direct problem solu-

tion) are solved under the Born approximation. The geometry of such experiments

is represented at Fig. 11.3. All results are normalized on a dominant wavelength λ

of the incident p-wave in the reference medium (v

=4000 m/s, v

=2000 m/s).

The velocity of the p-wave inside diﬀractors is equal to 3300 m/s. The recovery of

the set of diﬀractors is implemented by the processing of the several experiments

with the various normal to the wave front and the various signal shape. The recov-

ery consists in the estimation of the parameter w = (v

−1). In the experiment

the following normal vectors are used: n

, n

and n

. The next equalities are valid

, e

) = 1, (n

, e

) =

√

2/2, (n

, e

) =

√

2/2, where e

and e

are unit vectors,

directed along z and x axis respectively. For the simulation of the real seismic

data the white noise with the zero mathematical expectation and the variance σ

is added.

The resolution of the diﬀraction tomography method is investigated depend-

ing on the distance between receivers, the type of the incident wave, noise

level.

At Fig. 11.4 the results of the recovery for distances between receivers d = 0, 5

λ, d = 3 λ and d = 5 λ are represented. The noise level is 5 %. Two incident plane

waves are used for the calculation of the synthetic seismogram. The best recovery

is obtained for the case of d = 3 λ , Fig. 11.4(b). The better recovery is observed

in the case of two plane incident waves in comparison with one plane wave. The

applied algorithm can be used for the case of the noise level 20 %.

Tomography methods of recovering the image of medium 361

Fig. 11.3 The geometry of the numerical experiment: A is the plane (z = 0) contained the

geophones, B is the plane with diﬀractors (y = const), intervals between diﬀractors along z axis

are equal to 1.5 λ and intervals between diﬀractors along x axis are equal to 1.0 λ).

11.8.2 An estimation of the recovery accuracy of inhomogeneities

parameters

Let us consider the results of the numerical simulation on the estimation of the

parameters of local elastic inhomogeneities and inhomogeneities of an electrical

conductivity by the diﬀraction tomography method. In the cases of both types of

inhomogeneities, the direct problem on the elastic wave propagation (the solution of

the Lame equation (4.17), (10.11)) and the direct problem on the electromagnetic

wave propagation (the solution of the Maxwell equation (4.61)–(4.64)) is solved

by the ﬁnite diﬀerence method. It allows us to take into account the diﬀraction

phenomenon on the target object with the complex geometry and estimate the

accuracy of the diﬀraction tomography method depending on the contrast of the

target object relatively to the reference medium.

Let us consider the restoration of the elastic parameters. Using notion of the

tomography functional (see Secs. 10.2, 10.3), the vector of the diﬀerence ﬁeld (which

is the ﬁeld, scattered by the desired inhomogeneity — see Sec. 10.1) can be written

as follows

u(x

, x

, t) =

[δλ(x)p

(x, x

, x

, t) + δµ(x)p

(x, x

, x

, t)

+ δρ(x)p

(x, x

, x

, t)]dx, (11.66)

362 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Fig. 11.4 The results of the recovery for two incident plane waves with the normals n

and n

to the wave front (w = v

− 1): (a) — d = 0.5 λ; (b) — d = 3 λ; (c) — d = 5 λ.