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

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

Methods for parameter estimation 193

The variable t

belongs to the Student’s distribution with k = n − S degrees of

freedom. For given β and (a number of degrees of freedom) k to determine (for

example, using tabulated data) γ such that

P (|t

n−S

| ≤ γ) = β.

The conﬁdence interval has the form

± γˆσ

[ψ(ψ

ψ)

−1

]

1/2

After calculation of I

for i = 1, 2, . . . , n, we obtain the conﬁdence band, which

covers the function

f(θ) with the probability β (Fig. 6.5).

Fig. 6.5 Graphic representation of conﬁdence band I

Let us note, that the asymptotically eﬃcient estimates and the eﬃcient estimates

under suﬃciently general conditions have the lowest conﬁdence intervals.

6.16 The Method of Backus and Gilbert of the Linear Inverse

Problem Solution

At the solution of the inverse geophysical problems the model of observed data u

(i = 1, 2, . . . , n) can be represented in the form

(t)θ(t)dt + ε

, (6.72)

where ε

∈ N(0, R

), θ(t) is unknown desired function of the parameters, ϕ

(t) is

a given function. In a seismic case u

is a sample of the seismogram, ϕ

(t) is a

signal shape, θ(t) is a function, which characterized the reﬂectivity properties of a

medium, [0, T ] is an observation interval. The desired function θ(t) in point t

will describe by the mean value

hθ(t

)i =

C(t

, t)θ(t)dt, (6.73)

where C(t

, t) is an averaging function. If the averaging function is the delta func-

tion, then C(t

, t) = δ(t − t

) and the mean value hθ(t

)i is equal to the true value

194 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

θ(t

). Hence, under the practical calculation it is necessary to ﬁt the averaging

function close to the delta function, and to fulﬁl the normalization condition:

C(t

, t)dt = 1.

In the method of Backus and Gilbert the averaging function (Backus and Gilbert,

1967) appears as a linear combination of the signals ϕ

(t) with unknown coeﬃcients

C(t

, t) =

i=1

(t)ρ

). (6.74)

Substituting the expression (6.74) to the right hand side of the equality (6.73), we

obtain

hθ(t

)i =

i=1

, (6.75)

where

(t)θ(t)dt.

As it follows from the formula (6.75), the solution is represented as a linear

combinations of the functionals f

. Though, the random component ε

brings a

diﬀerence between the observed data u

and the model function f

. Therefore, as

the estimate, it is taken a linear combination of the observed data u

θ(t

)i =

i=1

. (6.76)

The problem is to ﬁnd ρ

) in accordance with the stated criterion. To construct

a decision rule we introduce the proximity function s(t

), which is a measure of

deviation between the averaging function C(t

, t) and δ-function:

s(t

) = 12

+L/2

−L/2

(t − t

)

, t)dt, (6.77)

where coeﬃcient 12 is chosen in accordance with the condition: if C(t

, t) is a

function with a squared shape (Fig. 6.6)

C(t

, t) =



1/L, if t ∈ [t

− L/2 ≤ t ≤ t

+ L/2],

0, if t 6∈ [t

− L/2 ≤ t ≤ t

+ L/2],

then s(t

) = L. Substituting the formulas (6.74) to the right hand side of the

equality (6.77), we obtain the proximity function

s(t

) = ρ

)S(t

)ρ(t

Methods for parameter estimation 195

Fig. 6.6 Graphic representation of C(t

, t).

where

ρ(t

) = [ρ

), ρ

), . . . , ρ

)],

S(t

) = kS

µν

µ,ν

µν

) = 12

+L/2

−L/2

(t)ϕ

(t)(t − t

)

dt.

The second component of the decision rule is the variance

) = h(hθ(t

)i − h

θ(t

)i)

i = h(

i=1

)(u

− f

))

= ρ

ρ(t

). (6.78)

At the development of the formula (6.78) the representations (6.75) and (6.76) are

used. The decision function looks usually as the following

g(ρ

ρ(t

)) =

ρ(t

), G = S cos θ + βR

sin θ,

where θ is a parameter belonging to the interval 0 ≤ θ ≤ π/2, β is a coeﬃcient

to ﬁll in the order of the deviation and the variance. If θ = 0, we get the decision

function, which is entirely determined by the deviation s(t

) and the standard de-

viation σ

). The deviation and standard deviation characterize correspondingly

the resolution and accuracy of the solution. We ﬁnd the desired vector using the

condition

ρ = arg min g(ρρ(t

)),

if BB

ρ(t

) = 1,

B = kB

ν=1

, B

(t)dt. (6.79)

Taking into account the theorem about minimum of the quadratic form with the

linear constraints, we can obtain

ˆρ

ρ(t

) =

−1

. (6.80)

A linear combination of the obtained estimate

ρρ(t

) with the values of the seismic

record allows to get the reﬂectivity property in the given point (t

196 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

6.17 Genetic Algorithm

Genetic algorithms (GA), ﬁrst proposed by John Holland (Holland, 1992) are based

on analogies with the processes of biological evaluation. GA are very important

methods for the solution of non-linear problems. The basic steps in GA are: coding,

selection, crossover, mutation.

6.17.1 Coding

Common to any basic GA is the digitization of the list of model parameters using

a binary coding scheme.

Consider, for example, the coding of the pressure wave velocity. The low velocity

limit of interest may be 1500 m/s and the upper limit 1810 m/s. Assume that the

desired resolution is 10 m/s. For this coding ﬁve bits is required (see Fig. 6.7). The

algorithm must now determine the ﬁtness of the individual models. This means

that the binary information is decoded into the physical model parameters and the

forward problem is solved. The resulting synthetic data is estimated, then compared

0 0 0 0 0

min

= 1500 m/s

0 0 0 0 1 v = 1510 m/s

0 0 0 1 0

v = 1520 m/s

0 0 0 1 1

v = 1530 m/s

··············· ···············

1 1 1 1 1

max

= 1810 m/s

Fig. 6.7 The binary coding of the values of velocity.

with the actual observed data using the speciﬁc ﬁtness criteria. Depending on the

problem, the deﬁnition of the ﬁtness will vary. For example, we represent the

normalized correlation function, recently used in the geophysical inversion,

F (θ) =

⊗ u

(θ)

⊗ u

)

1/2

⊗ u

(θ))

1/2

and the squared error function

F (

θ) =

i=1

− u

(θθ))

where ⊗ is a sign of the correlation, u

and u

(θθ) correspond to the observed data

and the model data for the parameter θθ.

Methods for parameter estimation 197

6.17.2 Selection

The selection pairs of the individual models for the reproduction is based on their

ﬁtness values. Models with the higher ﬁtness values are more likely to get the selec-

tion than models with low ﬁtness values. Let us consider the ﬁtness proportionate

selection.

The most basic selection method uses the ratio of each model’s ﬁtness function

to the sum of the ﬁtness of all models in the population to deﬁne its probability of

the solution, i.e.

(θ

) =

F (θ

)

j=1

F (θ

)

where n is the number of models in the population. The selection, based on this

probability, proceeds until a subset of the original models have been paired.

In a basic GA, if the population originally contained 100 models, ﬁfty pairs are

selected based on their ﬁtness values. We shall take the models, which satisﬁed to

the inequality

(θ

) ≥ δ,

where δ is a given threshold value. Let we obtain L models, which form L/2 pairs.

Each pair of the models will now produce two oﬀsprings using the genetic operation

of a crossover and a mutation. This will result in a completely new population of

individuals.

6.17.3 Crossover

A crossover is the mechanism that allows the genetic information between the paired

models to be shared. In the terminology of the geophysical inversion, the crossover

caused the exchange of some information between the paired models thereby gener-

ating new models. The crossover can perform in two modes: single and multi-point.

In the single point crossover, one bit positioning the binary string is selected at ran-

dom from the uniform distribution. All of the bits to the down of this bit are now

exchanged between the two models, generating two new models (see Fig. 6.8(a) In

multi-point crossover, this operation is carried out independently for each model

parameter in the string (Fig. 6.8(b)).

Example 6.1. The pair velocity model parameters v

and v

are the extreme mem-

bers for the coding scheme from 1500 m/s to 1810 m/s (Fig. 6.9 and Table 6.1).

6.17.4 Mutation

The last genetic operation is a mutation. The mutation is a random alternation of

a bit. It can be carried out during the crossover process. The mutation rate is also

speciﬁed by a probability determined by the algorithm designer.

198 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

(a)

←→

θθ

(b)

←→

Fig. 6.8 An example of the single-point crossover (a) and the multi-point (b) crossover (CP is a

crossover point).

0 0 0 0 0

= 1500 m/s

5 4 3 2 1 0

1 1 1 1 1 v

= 1810 m/s

Fig. 6.9 The pair of the extreme positions for the coding scheme.

Table 6.1 An example of the crossover po-

sition.

Crossover position v

, m/s v

, m/s

0 1500 1810

1 1510 1800

2 1530 1780

3 1570 1740

4 1650 1660

5 1810 1500

Methods for parameter estimation 199

Example 6.2. The process of the mutation for an initial velocity value of 1760 m/s

(Fig. 6.10).

5 4 3 2 1

1 1 0 1 0

= 1760 m/s

= 1740 m/s

Fig. 6.10 An illustration of the mutation procedure.

If the random value ε, obtained from the random generator, less the threshold value

, then the mutation procedure is produced, if ε > P

, then the mutation is not

produced.

6.17.5 Choice

We choose from each of the L/2 pairs a model, which has the larger ﬁtness function

(Fig. 6.11). Thus we obtain L/2 models, which form randomly L/4 pairs. Then we

1 2

F (θ

) F (θ

)

3 4

F (θ

) F (θ

)

···

L − 1 L

F (

L−1

) F (θ

)

If F (θ

) > F (

)

If F (

) < F (θθ

)

···

If F (θθ

L−1

) < F (θθ

)

| {z }

L/2

Fig. 6.11 A reduction of the models during GA.

200 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

produce the procedures: the crossover, the mutation and the choice. This procedure

is continued until we obtain the optimal model. The main advantage of GA is that

makes possible to ﬁnd the global maximum much faster than, for example, the

Monte-Carlo method. We note that the local linearization methods can give only

local maximum and they do not always lead to the global maximum.

Example 6.3. Test of GA on the synthetic seismogram.

In the speciﬁc case of seismic signals the ﬁtness function has the form

F (A

, . . . , A

, τ

, . . . , τ

) =

i=1

−

µ=1

ϕ(t

− τ

)

where the shape of the seismic wave is

ϕ = exp{−β|t|}cos ωt,

and A

is the amplitude, τ

is the arrival time of the seismic wave with the number µ.

Possible bound for each of the unknown parameters (A

, τ

) and the discrete interval

is determined a priori. The parameters of the attenuation β and the frequency ω

are ﬁxed also. The estimates of parameters to be ﬁnd by minimization of the ﬁtness

(or objective) function F (A

, . . . , A

, τ

, . . . , τ

During applying of GA to the optimization problem we choose the type of the bit

string encoding, the size of the working population L and the values of probabilities

of the crossover and the mutation.

Example 6.4. The size of the population, probabilities of the crossover and muta-

tions are respectively:

L = 50, P

= 0.6, P

= 0.001.

Algorithm is tested under the variety population size, the mutation and crossover

rates.

In a case of two interfering seismic waves the threshold of the correct wave

separation is about 1/4 wave period and signal-to-noise ratio must be equal or

greater than two. For the case of three interfering seismic waves, if |τ

− τ

| =

|τ

−τ

| ≥ 0.016 s. Three waves are separated with correct values of the estimated

parameters, which correspond to the value of the wave displacement between 1/4

and 1/3 wave period. The signal-to-noise ratio must be also greater or equal to two.

Chapter 7

Statistical criteria for choice of model

Before the processing of geophysical data the interpreter faces with a choice of

the suitable model of observed data. For some of these problems, for example, the

extraction of a seismic signal, the resolution of the seismic signals, the determination

of a number of interfering signals, the determination of the degree of the polynomial

etc., in mathematical statistics a number of the methods for the solution of such

problems are developed. Widely it is used the test of the parametric hypothesis with

the subjective setting of the signiﬁcance level. The information criterion, based on

the properties of the maximum likelihood estimates and the Fisher information, is

free from this imperfection. Let us start our consideration with the method of the

testing of parametric hypothesis.

7.1 Test of Parametric Hypothesis

The problem statement in the case of the test of parametric hypothesis usually is

the next: as a null hypothesis H

is regarded to the statement, that the desired

parameter vector θθ is equal to θθ

; the ﬁrst hypothesis (alternative hypothesis) H

consists in the statement, that

θ 6= θθ

. The estimate of the parameter vector can be

determined by the maximum likelihood method and it is tested: the consilience of

the hypothesis H

with the given error probability. Let on a basis of the model of the

medium with the parameter vector, the model ﬁeld uu

is calculated. As the result

of the observation we obtain the measured ﬁeld uu, which is described by the model

(3.2). To ﬁnd the estimate

θ using the maximum likelihood method. Practical

signiﬁcance has a problem of adjustment of the model parameter vector θ

with

the obtained estimate

θ. This problem is a typical one for the test of parametric

hypothesis at that, the null hypothesis H

consists in the statement

θ =

As a decision criterion we shall use the method of likelihood ratio, which has at

an enough great number of observations possesses the same optimal properties as

the maximum likelihood method. Moreover, for enough great size of the measured

data uu, the criterion of likelihood ratio possesses an important asymptotic property:

if θ is the S-dimensional vector and the likelihood function is regular in a sense of

201

202 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

the ﬁrst and the second derivatives on

θ, then in the case of the validity of the null

hypothesis and at enough great sample size, the distribution of the likelihood ratio

λ = −2 ln

L(u, θ

)

L(u,

θ)

tends to χ

distribution with S degrees of freedom.

The application of the criterion of the likelihood ratio corresponds in this case

to the next procedure:



if λ < χ

α,S

, then H

does not contradict u,

if λ > χ

α,S

, then H

is rejected,

where χ

α,S

is the value of the function χ

with a signiﬁcance level α and a number

of degrees of freedom S.

For the model

u =



either

f(θ

) + εε

: H

∈ N(0, σ

I),

f(θ

) + εε

: H

∈ N(0, σ

I),

(7.1)

where ε

, ε

are the random vectors with the normally distributed and independent

components with unknown variances, σ

and σ

. The criterion of the likelihood

ratio can be written as

λ = −2 ln



L(u,

, ˆσ

)

L(u,

, ˆσ

)



. (7.2)

The estimates ˆσ

and ˆσ

we ﬁnd by the maximization of the appropriate likeli-

hood functions L(

, σ

) and L(

, σ

ˆσ

(uu −

))

, ˆσ

(

u −

f(θ

))

Substituting the estimates σ

, σ

in the expression (7.2), and after a simple

transformation we obtain

λ = n ln

ˆσ

= n ln



1 +

ˆσ

− ˆσ

ˆσ



. (7.3)

The criterion χ

for the function λ reduces to the following: if λ < χ

α,S

, then

obtained estimate

θ does not contradict with the signiﬁcance level α to the hy-

pothesis H

θ =

, if λ > χ

α,S

, then the model θ

θ =

does not ﬁll in obtained

estimates and the null hypothesis is rejected. The value of χ

α,S

is determined using

tabulated quantities of χ

with a number of degrees of freedom S and a signiﬁcance

level α.

7.2 Criterion of a Posteriori Probability Ratio

Along with considered above solution of the problem of test of statistical hypothesis,

using the properties of the criterion of likelihood ratio, at a solution of the problems