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

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

Computer exercises 413

Fig. A.10 The convolution of the signals from Figs. A.8(a), A.9(a) and inverse cepstrum transform

of the sum of their cepstrums.

A.3 Direct and Inverse Problem Solution

A.3.1 Computer simulation of gravitational attraction

A.3.1.1 Gravitational attraction of the sphere

Script f0sph1.m is realized the gravitational attraction of the uniform sphere.

Gravitational potential V of the uniform sphere is given by

V = γ

πa

, (A.12)

where γ is the gravitational constant; a is a radius of the sphere; ρ is a mass density;

r is a distance between a center of the sphere and an observation point.

Thus the gravitational attraction reads as

g = −γ

πa

n, (A.13)

where n

n is a unit vector, directed from a center of the sphere to an observation

point.

A.3.1.2 Gravitational attraction of the uniform cylinder

Script f0cyl1.m is realized the gravitational attraction of the uniform cylinder.

Gravitational potential of the uniform cylinder V is given by the formula

V = 2πa

ργ log

, (A.14)

where γ is the gravitational constant; a is a radius of the cylinder; ρ is a mass density;

r is a perpendicular distance between a center of the cylinder and an observation

point.

Thus the gravitational attraction reads as

g = −

2πa

ργ

n, (A.15)

where n is a unit vector, directed from a center of the cylinder to an observation

point.

414 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

A.3.2 Computer simulation of magnetic ﬁeld

A.3.2.1 Magnetic induction of a dipole

Script f0dip1.m is realized the magnetic induction of the dipole.

Magnetic potential V of the magnetic dipole is given by

V = C

(m · r)

, (A.16)

where C

is equal µ

/4π in SI units;

m is a dipole moment; r is a distance between

the dipole and the observation point. Thus the magnetic induction reads as

BB = C

[3(n

·nn

−nn

], (A.17)

where m = |

m|,

= m

m/|

m|, n

= r

r/|

r|.

A.3.2.2 Magnetic induction of the horizontal cylinder

Script f0mcyl1.m is realized the magnetic induction of the magnetic cylinder.

Magnetic potential of the cylinder is given by

V = 2C

(

m ·

)

, (A.18)

where C

is equal µ

/4π in SI units;

m is the dipole moment per unit length;

r = |r

r| is a distance between the center of the cylinder and the observation point;

r/|

r|.

Thus the magnetic induction reads as

B =

[2(n

· n

−nn

], (A.19)

where m = |m|, n

= mm/|m|.

A.3.3 Computer simulation of the seismic ﬁeld

In the script fd0mn1.m the numerical simulation of the propagation of the elas-

tic waves in the piecewise half-space (1-D case) is realized by the ﬁnite diﬀerence

method. The component of displacement ϕ is satisﬁed the partial equation

(λ + 2µ)

∂

∂z

∂(λ + 2µ)

∂z

∂ϕ

∂z

+ f = ρ

∂

∂t

, (A.20)

the boundary condition at the free surface (z = 0)

∂ϕ

∂z



z=0

= 0 (A.21)

and at welded interfaces (z = z

)

ϕ|

z=z

−0

= ϕ|

z=z

, (λ + 2µ)

∂ϕ

∂z



z=z

−0

= (λ + 2µ)

∂ϕ

∂z



z=z

, (A.22)

Computer exercises 415

where ρ is a mass density, λ and µ are Lame parameters. The wave ﬁeld is excited

by the point source

f = δ(z)

f(t), (A.23)

located in the point z = 0.

The piecewise homogeneous medium can be perturbed by a smooth inhomo-

geneity of E = λ + 2µ parameter, which is given by the formula

δE(z) =



δE = 0, if |z − ˜z| > ∆,

δE = E

0.5[1 + cos(π(z − ˜z)/∆)], if |z − ˜z| < ∆,

(A.24)

or the smooth inhomogeneity of ρ parameter

δρ(z) =



δρ = 0, if |z − ˜z| > ∆,

δρ = ρ

0.5[1 + cos(π(z − ˜z)/∆)], if |z − ˜z| < ∆,

(A.25)

where ˜z is the location of inhomogeneity, ∆ is its half-size, E

and ρ

are the

maximum values of perturbation of E and ρ respectively.

A.3.3.1 Exercises

To calculate the seismic wave ﬁeld for the piecewise homogeneous half-space:

• the source location is z

= 0.003 km;

• the receiver locations are z

= 0.003, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5

km;

• location of the interfaces z=0.3, 0.5, 0.7, 0.9 km;

• Lame parameters E = λ + 2µ = 1;

• mass density ρ = 1.0, 1.2, 0.8, 1.2, 1.0;

• time sample dt =0.001 s;

• as the source time dependence the Ricker wavelet of 25 Hz frequency is used.

In Fig. A.11 the model, wavelet and the ﬁrst trace are represented. The seismogram

is represented Fig. A.12 to implement a qualitative analysis of the wave ﬁeld.

A.3.4 Deconvolution by the Wiener ﬁlter

Script p0decon2.m gives an example of deconvolution using the Wiener ﬁlter (see

Sec. 12.10)

Let consider processing of the seismogram (see Sec. A.3.3, Fig. A.12) with the

additional noise by the script p0decon2.m. We will try to decrease the signal-to-

noise merit and to bring a shape of the signal to the δ-function. In Fig. A.13 a

fragment of the seismogram from Fig. A.12 with addition of the Gaussian noise is

represented (N(0, σA

). Here σ=0.1, A

is a maximum value of the seismic signal

on the seismogram). In Fig. A.14 the fragment of the seismogram from Fig. A.13

after the deconvolution is represented.

The Wiener ﬁlter can be considered as an example of the application of the method of least

squares Sec. 6.4.

416 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Fig. A.11 Graphic presentation of the model and the ﬁrst trace. Time dependence in the source

(a); 1 — ρ, 2 — λ + 2µ (b); longitudinal wave velocity (c); the wave ﬁeld in the ﬁrst observation

point (d).

Fig. A.12 The seismogram.

A.3.4.1 Exercise

To implement the deconvolution of the seismic ﬁeld:

• without the noise (sig1=sig2=0);

• for diﬀerent values of the regularizing coeﬃcient α (see formula (12.36)) (at

eps1=2.0, 1.0, 0.5);

Computer exercises 417

Fig. A.13 The fragment of a seismogram from Fig. A.12 with the addition of the uncorrelated

Gaussian noise.

Fig. A.14 The deconvolution of the seismogram from Fig. A.13.

• for diﬀerent values of the noise (sig1=sig2=0.5, 0.1).

A.3.5 Quantitative interpretation

Let us consider an application of the Newton–Le Cam method (see Secs. 3.2, 6.2,

6.3) to the point estimation of the parameters θ of the geophysical objects (script

418 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

p0mg1.m ). As an examples, the magnetic bodies with vertical magnetization are

considered:

thread of poles

(θθ) =

2Mh

+ (x

− ξ)

θ = {M, h, ξ}; (A.26)

thread of dipoles

(θθ) =

2M[h

− (x

− ξ

)

]

+ (x

− ξ)

]

, θθ = {M, h, ξ}; (A.27)

point magnetic pole

(θ) =

M(2h

− x

)

+ x

)

5/2

θ = {M, h}; (A.28)

point magnetic dipole

(θθ) =

+ x

)

3/2

, θ = {M, h}; (A.29)

where M is a magnetic moment, h is a depth, ξ is a coordinate of the function f

(θθ)

maximum.

Applying the maximum likelihood method to the additive model of the data

= f

(θθ) + ε

(u = f(θ) + ε), (A.30)

with uncorrelated Gaussian noise ε

(N(0, σ)), we can get the next procedure for

estimation of

θ = (θ

, . . . , θ

). At the initial point to choose an initial value

θθ. After that to solve a system of the linear equations (see Sec. 6.3)

C∆θ

(1)

= d, (A.31)

where



u −f

θ))

∂f(θ)

∂θ





θ=

(0)

˜c

∂f

(θ)

∂θ



θ=θ

(0)

∂f(θ)

∂θ



θ=θ

(0)

s = 1, . . . , S

and to calculate the ﬁrst

(1)

approximation of θ

(1)

(0)

+ ∆

(1)

. (A.32)

The iterative process is ﬁnished if the threshold condition is satisﬁed, for example,

|∆

(n)

| ≈ 10

−2

÷ 10

−3

. (A.33)

An example of the similar numerical simulation for the ﬁeld given by the formula

(A.26) is represented in Fig. A.15 for the next input data:

Computer exercises 419

Fig. A.15 Qualitative interpretation. The model ﬁeld without noise (a) and with noise (b), the

model ﬁeld (1) and the ﬁeld (2) computed for the initial parameters (c), the model ﬁeld and the

ﬁeld computed for the estimated parameters (d).

x=-4:0.1:4; % observation points

m=1; % intensity

h=1; % depth

xi=0.0; % location of the object

xm=0.0; % expectation of the noise

s=0.1; % standard deviation of the noise

m0=0.5; % initial value of m

h0=1.5; % initial value of h

xi=0.; % initial value of xi

eps m=1.0e-5; % threshold condition for m

eps h=1.0e-5; % threshold condition for h

The values of M and h are estimated. As the result of numerical simulation,

the Fisher matrix and its inverse matrix are calculated.

A.3.5.1 Exercises

Using the script p0mg1.m to implement the next exercises for the models (A.26)–

(A.29) (to pay attention to the Fisher matrix and its inverse matrix):

(1) To investigate restoration of M and h for the diﬀerent values of the standard

deviation of the Gaussian noise (for example, σ =0.1, 0.2, 0.3).

(2) To investigate the convergence of the iterative process in dependence on initial

values of the parameters M and h.

(3) To implement the numerical simulation for the diﬀerent sample sizes (for ex-

ample, 10, 20, 50, 100).

420 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

A.3.6 Qualitative interpretation

Using script p0mg2.m to implement the qualitative interpretation (see Sec. 3.4, sig-

nal and noise (H

) or noise (H

)) with the help of a posteriori probability ratio

(see Secs. 7.1, 7.2) for the signal models from (A.26)–(A.29). For the case of the

uncorrelated Gaussian noise this criterion reads as

α = ln

P (1)

P (0)

2σ



u − (u − f

)

(

u −

)



. (A.34)

If α ≥ 0, then hypothesis H

(“noise”). If α < 0, then hypothesis H

(“signal and

noise”).

We assume that a model function f

(

θ) diﬀers from the real function f

(θ +∆θ).

Such diﬀerence is given by the deviation ∆θθ of the “true” parameters θ + ∆θ from

their supposed values θ.

Fig. A.16 Qualitative interpretation. The model ﬁeld f

(

θ) (1) and “real” ﬁeld f

(

θ + ∆

θ) (2)

(a), the model ﬁeld with uncorrelated noise (b), dependence of α (c) and normalized α (d) as

functions of the noise value.

An example of the similar numerical simulation for the ﬁeld (A.26) is represented

in Fig. A.16 for the next input data:

x=-4:0.1:4; % observation points

m=1; % intensity

h=1; % depth

xi=0.0; % location of the object

ddm=0.0; % deviation of m

ddh=0.0; % deviation of h

ddxi=0.2; % deviation of xi

xm=0.0; % expectation of the noise

p0=0.5; % a priori probability of existence of the object

s=0.1; % initial standard deviation of the noise

s1=s:0.4:4; % the values of the standard deviation for the numerical simulation

Computer exercises 421

A.3.6.1 Exercises

Using script p0mg2.m to implement the next exercises for the models (A.26) – (A.29):

(1) to implement the numerical simulation for a wide range of a noise/signal ratio

(variable s1);

(2) to implement simulations for the diﬀerent deviations of ∆θθ between the “real”

and model ﬁelds;

(3) to implement the simulation for a diﬀerent number of the observation points

(for example, 10, 20, 50, 100).

A.3.7 Diﬀraction tomography

Script dt0mn1.m is implemented a numerical simulation to restore the medium pa-

rameters by the diﬀraction tomography method (Secs. 10, 11).

The 1-D direct problem solution on the propagation of the elastic waves in uni-

form space, which contains a smooth perturbation of the parameters comparable

in size with the wavelength, is solved by the ﬁnite diﬀerence method. The com-

ponent of the wave ﬁeld ϕ satisﬁes the partial equation (A.20) an the boundary

conditions (A.21), (A.22). The restoration of the Lame parameters λ, µ and mass

and mass density ρ is implemented by the diﬀraction tomography method using the

next algorithm.

We introduce the tomography functionals p

and p

, which are connected with

the values ρ and E = λ + 2µ correspondingly:

(z, z

, z

, t) = −

∞

out

(z, z

, t −τ)

∂

∂t

(z, z

, τ )dτ,

(z, z

, z

, t) = −

∞

∂

∂z

out

(z, z

, t −τ) ×

∂

∂z

(z, z

, τ )dτ, (A.35)

then, under the Born approximation, the scattering (or diﬀerence) ﬁeld ∆ϕ gener-

ated by the inhomogeneity, can be written as

∆ϕ ≈

δE + p

δρ)dl, (A.36)

where z

and z

are the source and receiver locations; ϕ

is the wave ﬁeld generated

by the source with its time dependence; ϕ

out

is the wave ﬁeld, generated by the

“artiﬁcial” source with time dependence as δ(t)-function, located in the observa-

tion point; L is the restoration region; δE, δρ are the desired parameters of the

inhomogeneity (perturbation of E and ρ relatively reference medium).

If, for example, a value of δρ ≈ 0 and we shall ﬁnd only δE, then the inte-

gral equation (A.36) after digitization can be presented as a system of the linear

equations

= δ

, (A.37)

422 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

where

is the desired vector and δδ

is a vector of the sampled diﬀerence ﬁeld.

The simplest regularization of the solution of (A.37) consist of the solution of

the next system of the linear equations:

P + εD

D)δδ

= P

, (A.38)

where ε is regularized parameter and D, for example, is an identity matrix.

Let us consider an example of the numerical simulation for the next model:

• z

= z

= 0.003 km;

• ρ = 1, E = 1, v =

E/ρ = 1;

the smooth (desired) inhomogeneity is given by the formula

δE(z) =



δE = 0, if |z − ˜z| > ∆,

δE = E

0.5[1 + cos(π(z − ˜z)/∆)], if |z − ˜z| < ∆,

(A.39)

where ˜z is inhomogeneity location, ∆ is its half-size, E

is a maximum value of a

perturbation of E.

˜z = 0.3 km, ∆ = 0.03 km, E

= 0.3; Ricker wavelet

F (ω) =





exp







exp{−2πiω/ω

} (A.40)

with the apparent frequency ω

= 25 Hz; the interval L of restoration is L = 0.21

– 0.39 km; regularization parameter ε = = 10

−5

Graphic presentation of the model and the results of the direct problem solution

is represented in Fig. A.17 and Fig. A.18 respectively. The result of restoration is

Fig. A.17 The data for simulation of the propagation process and the wave ﬁeld in the ﬁrst

observation point. This representation is similar to Fig. A.11.

represented in Fig. A.19.