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

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

Algorithms of approximation of geophysical data 223

At Fig. 8.2 an example of smoothing of gravity ∆g(x) along a proﬁle by the cubic

splines is represented. Obtained curves can be used in subsequent interpretation

steps.

Fig. 8.2 Smoothing of gravity along proﬁle: 1 is initial data, 2 — DP = 1, 3 — DP = 10

−4

The spline approximation is widely used for the approximation of the seismic

interfaces. In the special case of the solution of a direct problem by the ray method,

the stability of the iterative processes depends on the smoothness of the interfaces

and their ﬁrst derivative. The cubic splines can provide the high smoothness and

the continuity of the ﬁrst derivative.

8.2 Periodic and Parametric Spline Functions

The interpolation and approximation of the closed and periodic curves can be im-

plemented using cubic splines with the periodic edge condition (Sp¨ath, 1973):

) = f

K−1

), f

) = f

K−1

), f

) = f

K−1

). (8.20)

Relations (8.20) are the conditions for the function ϕ(x

ϕ(x

) = ϕ(x

), ϕ

) = ϕ

), ϕ

) = ϕ

). (8.21)

Let us consider the interpolation problem. Using the second (join of the ﬁrst

derivatives) condition from (8.20), to write down the equation for the system (8.10)

(k = 2) in the form

2(∆x

K−1

+ ∆x

)ϕ

) + ∆x

K−1

)

= 6



∆u

∆x

−

∆u

K−1

∆x

K−1



. (8.22)

224 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

If k = K − 1 we obtain from the system of equations (8.10) with the use of the

condition (8.21) the following equation

∆x

K−1

) + ∆x

K−2

) + 2(∆x

K−2

+ ∆x

K−1

)

×ϕ

K−1

) = 6



∆u

K−1

∆x

K−1

−

∆u

K−2

∆x

K−2



. (8.23)

Taking into account the equations (8.22) and (8.23), we rewrite the system of equa-

tion (8.12) in the form

where

= [ϕ

), ϕ

), . . . , ϕ

K−1

)]

is a column vector;



∆u

∆x

−

∆u

K−1

∆x

K−1



, 6



∆u

∆x

−

∆u

∆x



, . . .

. . . , 6



∆u

K−1

∆x

K−1

−

∆u

K−2

∆x

K−2



is a column vector;



2(∆x

K−1

+ ∆x

) . . . ∆x

K−1

. . . . . . . . .

∆x

K−1

. . . 2(∆x

K−2

+ ∆x

K−1

)



The matrix A

is a cyclic three-diagonal matrix, and it is possible to prove its

positive deﬁniteness. The system of equations can be solved by the Gauss method.

At a solution of the interpolation problem of the closed curves we must depart

from the strict monotonicity of the data on the argument x

(see condition (8.1)).

We represents the curve by a set of arbitrary points (x

, u

), k = 1, 2, . . . , K, using

a parametric representation by two functions x = x(v) and u = u(v). At that, we

assume a strict monotonicity for the parameter v (v

< v

< ··· < ··· < v

We shall describe the curve using the cubic splines of two modes: (v

, x

) and

, u

), k = 1, 2, . . . , K. We suppose that the values of ﬁrst (second) derivatives

, ϕ

) by τ in the points v

and v

are given. In the practice, the parameter

v usually is determined in the following way:

= 0, v

= v

k−1

+ Γ

k−1

, k = 1, 2, . . . , K, (8.24)

where

∆x

+ ∆u

This method can be generalized to a three-dimensional case (x

, u

, v

). The pa-

rameter v is represented as

= 0, v

= v

k−1

∆x

+ ∆u

+ ∆v

, k = 2, . . . , K,

and for an each component the spline function is calculated. At the interpolation

of the closed curve in a plane or in a space it is necessary to specify the boundary

conditions: x

= x

, ϕ

= ϕ

and to use the parametric splines.

Algorithms of approximation of geophysical data 225

Let us consider the approximation problem of the non-monotone and closed

curves. We will approximate the corves given by the coordinates u

, v

, containing,

in the general case, the random component, with the help of the cubic splines v

, ψ

, ϕ

. The values of v

are the arguments and they calculated by the formula

(8.24). The coeﬃcient of proportionality between a jump of the third derivative

and a deviation of the observed and approximated data is chosen, in general case

for each coordinates separately:

000

) = g

− ψ(t

)), ψ

000

) = h

− ϕ(t

)).

For the approximation of the closed curves it is possible to use the periodic spline

functions for the variables u and v respectively. As the parameter we will consider

, which is determined by the formula (8.24) (instead of ∆x

we use ∆v

). Let’s

the next boundary condition is speciﬁed:

= u

, p

= p

, v

= v

, g

= g

, (8.25)

as well as

= ϕ

, ϕ

= ϕ

, ψ

= ψ

, ψ

= ψ

. (8.26)

Using the conditions (8.25) and (8.26), we obtain the system the same type as (8.19)

for the variables ϕ

, . . . , ϕ

k−1

, ϕ

, . . . , ϕ

k−1

, ψ

, . . . , ψ

k−1

, ψ

, . . . , ψ

k−1

. We ﬁnd

the solution of the system by the method described above.

Let us consider the results of the approximation of isolines of the complete

magnetic-ﬁeld vector (Fig. 8.3). From the analysis of Fig. 8.3 we can conclude that

Fig. 8.3 Smoothing of isolines of the complete magnetic-ﬁeld vector by the parametric splines at

the various weight coeﬃcients: 1 is initial data; 2 — DP = 10

−3

; 3 — DP = 10

−5

the smoothness degree is increased together with decreasing of the values of the

weight coeﬃcients.

226 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

8.3 Application of the Spline Functions for Histogram Smoothing

One of the tasks of mass processing of the results of geophysical observation is the

task of a histogram smoothing with the purpose to obtain a diﬀerentiable density

function.

Let us consider that histogram is described by the abscissas x

< x

< . . . x

and the corresponding staircase function, having in the interval [x

, x

k+1

] the or-

dinate H

, k = 1, 2, . . . , K − 1. The problem consists in the construction of the

smooth (diﬀerentiable) curve, if the integral of this curve at the interval [x

, x

k+1

]

is equal to a square of the appropriate rectangle of the histogram (Sp¨ath, 1973).

Let us suppose also, that for the each abscissa x

we specify the ordinate ϕ(x

which is crossed by the curve. To determine the coeﬃcients of the fourth order

polynomial in the interval [x

, x

k+1

(x) = a

(x − x

)

+ b

(x − x

)

+ c

(x − x

)

+ d

(x − x

) + l

(8.27)

under the following conditions

) = ϕ(x

), f

k+1

) = ϕ(x

k+1

), f

) = ϕ

), f

k+1

) = ϕ

k+1

(x)dx = H

∆x

Using these expressions to ﬁnd the polynomial coeﬃcients

∆x



∆x

(ϕ

k+1

) − ϕ

))

− 3(ϕ(x

k+1

) + ϕ

)) − 2H



= −

∆x



∆x

(ϕ

k+1

) − 3ϕ

))

− 7ϕ(x

k+1

) − 8ϕ(x

) + 15H



2∆x

[∆x

(ϕ

k+1

) − 3ϕ

) − 8ϕ(x

k+1

)

− 12ϕ(x

)) + 20H

], d

= ϕ

), e

= ϕ(x

The continuity condition for the second derivative is written down as f

k−1

) =

). Involving the representation (8.27) as well, we obtain the system of equa-

tions for determination of the ﬁrst derivatives ϕ

∆x

k−1

) + 3



∆x

k−1

∆x



)

−

∆x

k+1

= 20



∆x

−

k−1

∆x

k+1



+ 8

ϕ(x

k−1

)

∆x

k−1

+12



∆x

k−1

∆x



− 8

ϕ(x

k+1

)

∆x

. (8.28)

Algorithms of approximation of geophysical data 227

Now it is necessary to determine the edge conditions for the ﬁrst derivatives in

the points x

, x

: ϕ

), ϕ

). The system of equations (8.28) has a three-

diagonal and positive deﬁned matrix of the coeﬃcients. So, the solution of the

system possesses the unique existence properties.

For the calculation of ϕ(x

) can propose the next expedient:

ϕ(x

) = g

, ϕ(x

) = g

(

k − 1) + (1 − g

, ϕ(x

) = g

K−1

where g

(k = 1, 2, . . . , K) are the arbitrary coeﬃcients from the interval [0, 1].

Sometimes the values g

are chosen equal to 1/2 for all intervals.

Let us consider an example of smoothing of the velocity histograms, obtained

by the seismic data for two beds (Fig. 8.4). The obtained curves can be used for

Fig. 8.4 Histogram smoothing using the spline functions. (a) corresponds to the bed number one

1, g

= 0, 3, 0, 3, 0, 5, 0, 6, 0, 3, 0, 3, 0, 2; (b) corresponds to the bed number two, g

= 0, 0, 6, 0, 5,

0, 4, 0, 1, 0, 7, 0, 4, 0, 6, 0, 3, 0, 2, 0, 5, 0, 3, 0, 3.

the calculation of moments of the distribution. The histogram smoothing can be

applied to other processing stages as well.

8.4 Algorithms for Approximation of Seismic Horizon Subject to

Borehole Observations

Let the area under the study is covered by an arbitrary set of proﬁles registered

by the multifold coverage scheme. (Troyan and Sokolov, 1982). In addition, there

are a few deep boreholes with boring data. As the result of the processing of the

proﬁles we obtain the depth proﬁles. Let us consider that we have the interlocked

proﬁles, i.e. in the cross points of the proﬁles the depth deviations do not exceed

the threshold values. As the initial observed data

U = [uu

, y

)] we consider

the bed depths corresponding to the depth proﬁle coordinates: x

, y

; k are the

points of the area under investigation, k = 1, 2, . . . , K; m is a number of a bed,

m = 1, 2, . . . , M . In this case the model of observations can be written down as

UU = ψρρ + n, (8.29)

228 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

where

U =



. . .



, ρ =



. . .

ρρ



, nn =



. . .



ψ = [ψ

, ψ

, . . . , ψ

ρρ

is the parameter vector, ψ

is the structural matrix of m-th bed. For example,

in the case of the polynomial of power two, we have



1 x

. . . . . . . . . . . . . . . . . .

1 x



The random component nn

= [n

, y

)] in the points x

, y

is centralized,

and the correlation between the beds is determined by the covariance matrix



. . . R

. . . . . . . . . . . .

. . . R



. (8.30)

The values of the bed depth

V in the boreholes are assumed to be under linear

constraints, which determine the parameters of approximating planes:

aρ

ρ =

V ,

where

A = [A

, A

, . . . , A

], V =



. . .



is a set of the depth for m-th bed obtained using borehole data; A

is the

structural matrix for m-th bed, which one, for example, for the polynomial of power

to has the following shape



1 x

. . . . . . . . . . . . . . . . . .

1 x



The estimate

ρ is determined by the formula (see expressions (6.23), (6.25))

+ (

−

ρA

)(AB

−1

)

−1

, (8.31)

where

ρ = B

−1

UU, B = ψ

−1

ψ, (8.32)

Algorithms of approximation of geophysical data 229

and it is the estimate of the ρ by the least squares method without linear constraints.

The covariance matrix of the estimate of parameters

ρ can be represented as

(see formulas (6.26), (6.27))

˜ρ

= R

ˆρ

− B

−1

(AB

−1

)

−1

, R

ˆρ

= B

−1

where R

ˆρ

is the covariance matrix of the estimate of parameters

ˆρ

ρ, obtained by the

least squares method without linear constraints.

Let us consider the variants of introducing of the correlation dependence between

the random components of the model (8.29).

8.4.1 The Markovian type of correlation along the beds and no

correlation between beds

The matrix R

(8.30) in this case is a diagonal one and its inverse matrix R

−1

reads

as:

−1



−1

0 . . . 0

0 R

−1

. . . 0

. . . . . . . . . . . .

0 0 . . . R

−1



The matrices R

−1

are the inverse matrices with respect to the matrix of a type of

(3.22):

= σ



1 β

. . . β

1 β

. . . β

K−1

. . . . . . . . . . . . . . .

K−1

K−2

. . . 1



Here, the value of β

have a sense of correlation coeﬃcients between the points

with the distance of ∆r =

(∆x)

+ (∆y)

for the bed with the number m. The

inverse matrix R

−1

reads as

−1

(1 − β

)



1 −β

0 . . . 0 0

−β

1 + β

−β

. . . 0 0

0 −β

1 + β

. . . 0 0

. . . . . . . . . . . . . . . . . .

0 0 0 . . . −β



. (8.33)

Using the representation (8.33), we ﬁnd the matrix B:

(1 − β

)

(1 + β

)

mks

− β

(ψ

mks

m(k+1)s

+ ψ

m(k+1)s

mks

)

− β

(ψ

mKs

+ ψ

m1s

)

. (8.34)

230 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

8.4.2 Markovian correlation between the beds and no correlation

along bed

By repeating above mentioned arguments, we pass to the next representation of the

matrix B, which is similar to (8.34):

(1 − β

)

[(1 + β

)ψ

mks

− β

(ψ

mks

(m+1)ks

+ ψ

(m+1)ks

mks

) − β

(ψ

1ks

+ ψ

Mks

)]

8.4.3 Conformance inspection of seismic observation to borehole

data concerning bed depth

The next problem represents great practical signiﬁcance. Let using the maximum

likelihood method, the estimation of the parameters

ρ of approximated surface with

taking into account and without taking into account the borehole observation are

found. And it is necessary to solve the problem about the conformance of seismic

observations to the borehole data concerning the bed depth. This problem can be

solved by the test of the parametric hypothesis method (see Sec. 7.1). We consider

the linear model (8.29) (we put M = 1 and R

= σ

I for the simplicity sake) and

we test the hypothesis H

: Aρ =

V . For the hypothesis checking we shall use the

method of the likelihood ratio

q =

˜ρ

, ˜σ

)

ˆρ

, ˆσ

)

. (8.35)

As it was demonstrated in Sec. 6.2, the estimates

ρ obtained by the least squares

method and the estimates

ρρ

obtained by the maximum likelihood method coin-

cide in the case of the normal distribution of the random component n

To ﬁnd the estimates ρ

ρ under the condition of the validity of the hypothesis H

ρ =

V . In this case the maximum likelihood method reduces to the conditional

extremum problem, which can be solved by the Lagrange multipliers method. To

construct the function

F (ρρ) = ln L + 2λ

(Aρ −VV ), (8.36)

where λ is a column vector contained undetermined Lagrange multipliers. The

estimate

we ﬁnd by the minimization of (8.36) under the condition, that the

equality Aρ =

V is valid:

∂F/∂

ρ = 0. (8.37)

After simple transforms we obtain

+ 2σ

A(ψ

ψ)

−1

. (8.38)

To postmultiply the left hand side and right hand side of the equality (8.38) by A

and to use the equality

˜ρ

= V

, we obtain

λλ

2σ

(

−ρ

)(A(ψ

ψ)

−1

)

−1

. (8.39)

Algorithms of approximation of geophysical data 231

Substituting λ

from the formula (8.39) to the expression (8.38), we obtain

+ (VV

− ρ

)(A(ψ

ψ)

−1

)

−1

A(ψ

ψ)

−1

The obtained estimate ρ

˜ρ

is equivalent to the estimate

ρρ, given by the formula

(6.25). The estimate for the standard deviation can be written as

˜σ

(u − ψ

ρρ)

(u − ψ

ρρ).

Substituting the obtained estimates to the likelihood function we calculate its

maximum value for the hypothesis H

and its alternative:

max

(

ρ) = (2π

)

−K/2

, (8.40)

max

(

ρ) = (2π

)

−K/2

. (8.41)

Using the expressions (8.40), (8.41), to rewrite the formula (8.35) in the form

q =

max

(

˜ρ

ρ)

max

(

ρ)



ˆσ



K/2

. (8.42)

We use the criterion as follows: if q < q

, then the hypothesis H

is rejected, if

q ≥ q

, then the hypothesis H

is valid (and in our case the depth values calculated

using borehole data does not contradict to the depth values calculated using seismic

observations. To determine the threshold value of q

it is necessary to determine the

probability distribution of q, but the most simple way consists in the determination

of the constraint of the considered criterion with the Fisher information criterion.

If the hypothesis H

is correct one, then the relation

F =

[(

u − ψ

ρ)

− (uu −ψ

ˆρ

ρ)

]/Q

(uu − ψ

ˆρ

ρ)

/(K − S)

(8.43)

belongs to Fisher distribution (see Sec. 1.8.6) with Q and K −S degrees of freedom

(F ∈ Φ

Q,K−S

). By analyzing the formula (8.42) and (8.43), we can conclude that

the function F can be expressed through q:

F =

K − S

−2/K

− 1].

The likelihood ratio criterion is reduced to the Fisher criterion: if F > Φ

Q,K−S,α

(the value of Φ

Q,K−S,α

is calculated subject to given conﬁdence level α), then

the hypothesis H

is rejected and we conclude that with the given probability of

error, there is an inadequacy between the depth values obtained using borehole

data and depth values obtained using seismic data. The value Φ

Q,K−S,α

can be

determined using, for example, a tabular information on the Fisher distribution at

given Q, K − S and the probability of the error (conﬁdence level).

232 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

8.4.4 Incorporation of random nature of depth measurement using

borehole data

The above considered problem statement has been limited by an exact setting of the

bed depth using borehole data. Such assumption is justiﬁed, if the accuracy of the

seismic observation is more rough in comparison with the borehole measurement.

However, the consideration of the random nature of the depth measurements using

the borehole data has a practical signiﬁcance.

The model of observations in this case can be represented as: VV = A

ρ + ε. Here

VV , A,

ρ has the former sense,

ε = [

ε(x

, y

)] is a random component of the model,

measured at the points (x

, y

) with the covariance matrix P

The estimate of the parameter ρ we ﬁnd using minimizing of the following square

form

l(ρρ) = (uu − ψ

ρ)

−1

(u − ψρ) + (VV − A

ρ)

−1

(V −Aρρ). (8.44)

The formula (8.44) is the objective function in the parameter space ρρ, jointly for

seismic and borehole observations. The estimate of the parameter vector ρρ we ﬁnd

by minimizing of (8.44) with respect to ρ:

= (uu

−1

ψ + V

P ε

−1

A)(ψ

−1

ψ + A

P ε

−1

−1.

. (8.45)

The formula (8.45) can be transformed to the following shape

= (

ˆρ

)(E + B

)

−1

, (8.46)

where

ˆρ

ρ is determined by the formula (6.11):

= Ψ

−1

Ψ(A

−1

= V

−1

A(A

−1

The formula (8.46) can be simpliﬁed, if the following condition is valid

max

) < 1, (8.47)

where λ

max

is the maximum eigenvalue of the matrix B

In the most practical cases the accuracy of borehole data concerning the bed

depth is higher than the bed depth obtained with the use of seismic data. Analysis

of the matrix B

follows that in these cases the inequality (8.47) is valid and the

formula (8.46) can be rewritten in the approximation form:

≈ (

ˆρ

)(E − B

). (8.48)

In practice it is often to suppose the random components nn and εε are uncorre-

lated. In this case the formulas (8.46), (8.48) for the estimate ρρ become simpler:







−1

, (8.49)

ρ ∼



(Ψ

Ψ)(A

−1

ˆρ





E −

ˆρ

(Ψ

Ψ)(A

−1

