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

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

Tomography methods of recovering the image of medium 323

For a point source s in the homogeneous medium considered, the ﬁelds ϕ

and ϕ

out

are of the same structure, i.e.,

4π|x − x



t − t

−

|x − x



out

4π|x − x



− t −

|x − x



where x

and t

specify the location of the source s and the time of its engagement,

respectively; x

and t

describe the location of the receiver and the current time of

recording respectively. In this particular case, the tomographic functional is of the

form

(4π)

|∆x



− t −

|∆x



∂

∂t

|∆x



t − t

−

|∆x



|∆x

| = |x − x

|, |∆x

| = |x − x

Expanding |∆x

| and |∆x

| into series and retaining the ﬁrst terms of the expansions

in the denominators and the ﬁrst two terms in the arguments of the δ-function

(which corresponds to the Fraunhofer approximation), we obtain

p =

(4π)

||x

dt δ



− t −

−

· x



× δ



t − t

−

· x



(4π)

| · |x

× δ



− t

) −

(|x

| + |x

−

+ n

, x)



(4π)

||x

− n · x) ,

= (t

− t

− (|x

| + |x

|) ,

, n

, n = n

+ n

Here, the vector n corresponds to the normal to the plane relative to which the

beam ϕ

with vector n is mirror-likely reﬂected in the direction (−n

). Finally,

the scattering ﬁeld is related to the function ν(x) via the Radon transform, i.e.,

˜ϕ = hp|ν(x)i =

(4π)

||x

ν(x)δ

− n · x) dx

= const

∂

∂p

ν(x)δ(p − n · x)dx .

The Radon transform, formally represented as the integral over the points of the

plane {X : p = n · x}, provides, for diﬀerent values of the recording time (t

), the

following approximate representation of ˜ϕ:

˜ϕ =

(4π)

||x

[(t

− t

− (|∆x

| + |∆x

|)] .

324 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Here, for diﬀerent values of the recording time t

, it is implied that integration is

carried out over the corresponding ellipsoid of revolution with axes of the symmetry

passing through the points of the source (x

) and receiver (x

) (which are the focuses

of the ellipsoids), which leads to the following clear physical interpretation: the

integration is performed over sets of kinematically equivalent points in the space,

i.e., over level of the constant total phase.

We consider the generalized Radon inversion in the idealized diﬀraction tomog-

raphy on the more complicated example presented in (Beylkin, 1985). In scalar

wave equation in a reference medium with the velocity of propagation c

(x)



(x)

∂

∂t

− ∆



ϕ(x, t) = 0

can be represented in the form of the Helmholtz equation

(x) + ∆)ϕ

(x, k) = 0 ,

where k

(x) = ω

(x), ω is a cyclic frequency. Introducing the dimensionless

function µ(x) via the relations

µ(x) = −

ν(x) , ν(x) =

(x)



1 −

(x)



we write the Helmholtz equation in a medium with the velocity c(x) in the following

form:

(x) + µ(x)) + ∆]ϕ(x, k) = 0 .

The perturbation operator (an operator of multiplication) is of the form δL =

µ(x). We assume that the reference medium is such that the ﬁelds ϕ

and

out

are suﬃciently accurately approximated by the zero-order terms of their ray

expansion. Then the scattering ﬁeld is represented in the form

˜ϕ(k, x

) = k

hϕ

out

|µ(x)|ϕ

i = hp

|µ(x)i

= k

out

(x, x

ikτ (x,x

)

µ(x)A

(x, x

−kτ (x,x

)

dx .

The amplitudes A satisfy the transport equation, and τ satisﬁes the eikonal equa-

tion, so that

(

∇

∇τ,

∇

∇τ) = n

(x) ,



(∇∇τ, ∇) +

∆τ



A = 0 .

The tomographic functional p is represented by the integral kernel

p = k

out

(x) A

exp{ik[τ(x, x

) + τ (x, x

)]}

∆

= k

A(x) exp(ikτ ) .

We represent the generalized causal Radon transform associated with the rep-

resentation of the scattering ﬁeld ˜ϕ = hp|µi in the form

Rµ = hA(x)δ(q − τ (x))|µ(x)i,

Tomography methods of recovering the image of medium 325

and set R ≡ 0 for t < 0.

Note that the scattering ﬁeld ˜ϕ can be represented up to the factor as the Fourier

transform of the generalized Radon transform Rµ, i.e.,

˜ϕ = c F

Rµ .

Let the point x

of the source lies inside a closed region V with the boundary

∂V . We assume that the reference medium and the boundary ∂V are such that each

point of the boundary can be connected with the point of the source by a unique

ray. In addition, we assume that the rays do not intersect. In this case, the points

of the boundary ∂V can be parameterized by the corresponding points of the unit

sphere centered at the source point. Then the action of the conjugate operator R

∗

is determined in the following way:

∗

ψ =

∂V

ψ(q, x

q=τ (x

,x)

w(x, x

) dx

where w(x, x

) is a smooth nonnegative weight function, which is chosen in the

form

w(x, x

) = D(x, x

−1

(x, x

)h(x, x

) .

Here, h(x, x

) is a cut-oﬀ factor insuring the nonnegativity of w(x, x

), which

implies that the points of the boundary ∂V and of the unit sphere are in one-to-one

correspondence;

D(x, x

) =



∂

τ ∂

∂

x1 x

τ ∂

x2 x

τ ∂

x3 x

∂

x1 x

τ ∂

x2 x

τ ∂

x3 x



−1

(x, x

) = [A

out

(x, x

) A

(x) ]

−1

It can be shown (Beylkin, 1985), that D dx

= D

dΩ (Ω is a solid angle) and that

(x, x

) = n(x)[1 + cos θ(x)], where (cos θ(x) = e

out

·e

and e = e(x) is the unit

ray vector at the point x).

We write the generalized integral Fourier operator, deﬁned as follows:

Φµ(x) =

(4π)

∞

∂V

ikτ (x,x

)

˜a(x, x

, x

)µ(x

) dx

dk , (11.8)

where τ (x, x

, x

) = τ(x

, x

) − τ (x, x

, x

a(x, x

, x

) =

A(x

, x

)

A(x, x

)

D(x, x

) h(x, x

) .

Observing that the generalized Fourier operator (Φ) includes the Fourier operator

), where

: F

ϕ =

(2π)

3/2

∞

˜ϕ e

−ikq

dk ,

326 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

we can represent the generalized Fourier transform as the following composition of

three operators:

Φ = R

∗

In order to analyze the meaning of the generalized Fourier operator introduced

above, we write the linear term of the phase function:

τ(x, x

) ≈ ∇

τ(x, x

, x

)(x − x

) ,

∂V

µ =

(4π)

∞

∂V

ik∇

τ(x,x

)(x−x

)

× D(x, x

)µ(x

)dx

dk .

Introducing the new variable p = k∇

τ(x, x

, x

), we represent the latter integral

in the form

∂V

µ =

(4π)

Ω

(x)

ip(x−x

)

µ(x) dx

dp .

Here, Ω

(x) is the image of the domain of integration under the above change

of variable.

Examining formula (11.8), we see that the generalized Fourier operator provides

for the recovery of the spectrum of the function µ(x), which is determined by the

spectral domain Ω

(x). In its turn, the latter domain depends on the part of

the surface of the boundary ∂V

over which the scattering ﬁeld ˜ϕ is integrated,

and also on the nonsingularity of the Jacobian of the passage from the coordinates

on the surface (coordinates of recording) to the ray coordinates. Therefore, the

approximate operator equality

∂V

≈ R

∗

determines the reconstruction operator

µ ≈ R

∗

˜ϕ .

We note that in our opinion, the solution obtained is mathematically elegant.

But it does not actually concern to practical interpretation problems. Even un-

der our idealized assumptions on the relationship between the scattering ﬁeld and

the medium structure to be recovered, the linearized model, which is based on the

ray representation of the wave ﬁelds in the reference medium, has been used. The

Jacobian occurring in the generalized conjugate operator R

∗

certainly fails to be

nonsingular if the dependence of the ﬁeld ϕ

on time diﬀers from the δ-function.

This is easily seen from the generalized Radon transform, which involves the man-

ifold δ(q − τ(x)) of the dimension coinciding with that of the problem. Each value

corresponds to the integration over a three-dimensional domain. In the case of a

point source, this domain is the surface of the ellipsoid with focuses at the points

and x

if the temporal dependence is described by the δ-function and a layer of

an ellipsoid of revolution if the duration of the ﬁeld ϕ

is ﬁnite (see Fig. 10.1(b)).

Tomography methods of recovering the image of medium 327

11.3 Construction of Algorithms of Reconstruction Tomography

In order to construct a recovery algorithm, we write the model (10.9) using the

notion of a tomographic functional:

µ n

|ν

(δθ)i + ε

(11.9)

Here, u

is a sampling of the seismogram, ρ

µn

is the tomographic functional, ν

(δθ)

is an unknown function of the variation of parameters (µ is the index of the cor-

responding ﬁeld, µ = 1 ÷ M). The error ε

includes a random component, which

is generated by the ﬁelds of unidentiﬁed sources and instrument noises, as well as

the potential inadequacy of the model. The latter is caused, ﬁrst, by errors in

the representation of the apparatus function; second, by errors in the prior picture

of the medium, ignoring dissipative and nonlinear eﬀects in the description of the

propagation of sounding signals; third, by errors caused by the linearization of the

ﬁeld ϕ; fourth, by assumptions on the smoothness of the function θ(x), and so on.

Since, physically the stochastic structure of the noise ﬁeld can be regarded as the

sum of a large number of independent sources, it is natural to assume that the

random component is additive and normally distributed (Goltsman, 1971; Troyan,

1984, 1985). The statistical analysis of experimental data is generally conﬁned to

the ﬁrst and second moments of the random component. As it was mentioned in

Sec. 9.3, given the ﬁrst and second moment, the distribution of maximum entropy

is normal, which provides yet another reason for the choosing of Gaussian random

variables in the constructing of the model.

In constructing of practical algorithms, the most simple and frequently used

model (K

) is that of uncorrelated equal-accuracy errors, i.e., K

= σ

I, where

I is an identity matrix, and the only parameter σ

is the variance of the random

component. In the case of uncorrelated measurements with the unequal accuracy,

the correlation matrix is of the form



0 . . . 0

0 σ

. . . 0

. . . . . .

0 0 . . . σ



The simplest mathematical description of the correlated errors of measurements is

provided by the Markov’s model of a random process, where the correlation matrix

is of the form



1 γ γ

. . . γ

γ 1 γ . . . γ

N−1

γ 1 . . . γ

N−2

. . . . . . .

N−1

N−2

. . . 1



328 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

This matrix has the advantage that its inverse can be represented in the explicit

form, namely

−1

1 − γ



1 −γ 0 . . . 0

−γ 1 + γ

−γ . . . 0

. . . . . . . . .

0 0 0 . . . −γ

0 0 0 . . . 1



Note that the matrix K

−1

is the ﬁnite diﬀerence counterpart of the symmetric

positive deﬁned operator (−∆ + γI), which is frequently used in the Tikhonov’s

regularization (Tikhonov and Arsenin, 1977).

In what follows, we will assume that the covariance matrix of the measure-

ment errors is known. However, the simultaneous estimation of the distribution

parameters and of the desired parameters of the medium presents no great diﬃ-

culties, although additional is needed because, in this case, the estimation problem

becomes certainly nonlinear (Turchin et al., 1971).

Since the random component of ε is principally unremovable, any method of

solving inverse problems related to the interpretation of practical data must nec-

essarily deal with a problem of the statistical estimation. Note that regularization

methods make it possible to construct formally stable solutions of tomographic

problems, which represent the desired ﬁelds ν(δθ) at any point of the domain un-

der the investigation. The ﬁelds representing the medium can be considered as

multicomponent one (Tarantola, 1984). However, the solution of such problems in

interpreting data of remote sensing is useless because of the insuﬃciency of informa-

tion, i.e., practically a priori picture of the ﬁelds of parameters of the medium is not

reﬁned. The idea of recovering ﬁelds from values of a ﬁnite number of functionals

of measurement has been elaborated in geophysics (see (Backus and Gilbert, 1967,

1968)). In accordance with this approach, we try to replace the recovery of the

ﬁeld of parameters of the medium in the whole space by the recovery of values of

a number of linear functionals l of the ﬁeld ν(δθ). This is a well-posed problem

only in the case where the unknown functional l belongs to the linear hull of the

tomographic functionals {p

, n = 1 ÷ N }, i.e., to the subspace

∗

= {l : ∃α : (l − P

∗

α)ν(δθ) = 0 ∀ν}.

This regularization method has evolved spontaneously when solving practical in-

verse problems, and diﬀerent modiﬁcations of this method are commonly used,

although they are not properly justiﬁed.

If no a priori information on ν(δθ) is available, then the stability of the solution

obtained by the maximum likelihood method is determined by the Fisher’s opera-

tor, which is certainly degenerate. The degeneracy of the Fisher’s operator means

that there is a direction in the space Θ that also determines elements of the func-

tional space, the information on which is unavailable, whence the variance of the

corresponding linear functional is unbounded. On the other hand, if the functional

Tomography methods of recovering the image of medium 329

l is the linear combination of the tomographic functionals {p

} or, in other words,

if the conjugate equation P

∗

α − l = 0 is solvable, then the variance of its solution

is ﬁnite.

In practice, one frequently solves inverse problems by representing the solution

ν(δθ) in the form of a ﬁnite series:

|ν(δθ)i =

|ψ

ihψ

|ν(δθ)i ↔

|ψ

In this case, it is necessary that estimations of the coeﬃcients ν

have a ﬁnite vari-

ance, i.e., the equations P

∗

−ν

= 0 must be solvable. As the basis of {ν

} one uses

either a number of functions from a complete orthonormalized set (for example,

Fourier series, Fourier-Bessel series, Korhunen–Loeve basis and so on; the choice

of the basis is usually determined by an intuitive a priori idea of the medium),

or a certain incomplete basis (characteristic functions of nonoverlapping space do-

mains, B-splines etc. (Troyan, 1981a, 1983)), or a number of heuristic (incomplete,

nonorthogonal) basis functions of the kind of “model bodies” (for example, the ﬁeld

representations in the magnetic and gravimetric prospecting).

The error ∆ν of the recovery of the ﬁeld ν(δθ) when using the representation in

the form of ﬁnite series

|ν(δθ)i =

|ψ

is determined, ﬁrst, by errors ∆ν

of the recovery of the coeﬃcients ν

, i.e., of the

values of linear functionals hψ

|νi, and, second, by the component of the ﬁeld ν(δθ)

not involved in the above representations, i.e.,

∆ν(δθ) =

∆ν

|ψ

i + ˜ν ,

where ˜ν ∈ Θ =

I −

|ψ

ihψ

Θ .

Note that inﬁnite error related to

θ is usually ignored when analyzing solutions

of practical problems, i.e., it is intuitively assumed that the solution sought belongs

to the subspace corresponding to the basis chosen.

The problem of the recovery of an arbitrary linear functional can be made a

well-posed one if one invokes additional information by introducing either a (non-

degenerate) probability measure in the space Θ or deterministic (usually quadratic)

constraints. The introduction of a probability measure, for example, of the Gaus-

sian one with positive correlator K

, ensures that the variance σ

of an arbitrary

linear functional l becomes ﬁnite:

= l(K

∆l) = l((F + K

−1

)

−1

l) .

Here F = P

∗

−1

P is the Fisher’s operator and ∆l = l − P

∗

α. If one introduces

a deterministic constraint by specifying a closed domain in the parameter space,

330 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

for instance, the sphere ||ν(δθ)||

= 1, then, obviously, the variance of any linear

functional will be ﬁnite. But such an approach is sometimes senseless because the

error can prove to be of the same order as the diameter of a priori chosen domain,

although the equation P

∗

α = l is solvable (Morozov, 1993).

The regularization through the introducing of compacts in the space Θ, for

example, of a set of monotone bounded functions in L

(Bakushinski and Gon-

charski, 1994), is performed by applying methods of mathematical programming.

In contrast to the case of quadratic constraints, the corresponding algorithms are

multistep methods and they can not be represented by closed formulas

It should be noted that all the conclusions mentioned above are implied by the

identity

(l, ν) = l(ˆν), which means that the best estimate of a linear functional of the

ﬁeld ν coincides with the value of this functional at the best estimate of the ﬁeld. All

basic formulas and relations can be obtained in the same way as the optimal decision

operator, the construction of which was proposed in Sec. 9.3. Nevertheless, in order

to emphasize the similarities and diﬀerence in the Backus-Gilbert approach and that

suggested here, below we outline the algorithm for solving the inverse problem.

The problem is to recover the value of a given liner functional l of an unknown

ﬁeld ν

of parameters of the medium from N values u

of the sum of M linear

functional p

. For any linear functional l(ν), we seek a solution in the form of a

linear combination of measurements, i.e.,

l(ν) =

= (α, u) , (11.10)

where l(ν) ≡ hl|ui. Using the model of experimental data (11.9), in the vicinity of

a point X the linear estimation (11.10) can be represented as follows:

(δθ) =

h(α

, p

)|ν

(δθ)i

+ (α

, ε) .

Now we will reﬁne the meaning of the term “the ﬁeld in the vicinity of a point

X”. As ν(δθ

) one can accept, e.g., the mean value in the sphere of the radius ρ

centered at the point X, i.e., the value of the linear functional of the unknown ﬁeld

ν(δθ

(δθ

) = hl

|ν



πρ



−1

H (ρ − |x − X|) ,

H(ξ) =



0 if ξ < 0 ,

1 if ξ ≥ 0 .

Irrespective of the speciﬁc choice of the linear procedure, the error η

of the recovery

|ν

i can be formally represented as

= hl|ν

i −

(α, hp

)|ν

i) − (α, ε)

= hl − (α, p

)|ν

6=µ

h(α, p

)|ν

i − (α, ε) . (11.11)

Tomography methods of recovering the image of medium 331

The quantitative interpretation of this expression is as follows. The error η

depends on the unknown ﬁeld ν

(δθ) (the ﬁrst term in (11.11)) as well as on the

other ﬁelds ν

(µ

6= µ), and also on the random component ε. The less the norm

of the functional l − (α, p

), i.e., the better the approximation of the functional l

under the estimation by the linear combination of the tomographic functionals p

the less the component of the error that depends on the ﬁeld ν

. As is readily seen,

a reasonable approximation of the functional l, whose support is localized in the

vicinity of the point X can be obtained only if the set {p

µn

} includes functionals

whose supports have nonempty intersections with this vicinity. The error η

is the

sum of linear forms of the unknown ν

( µ = 1÷M) and ε, and each of these forms

depends on the coeﬃcient vector α, which speciﬁes a particular recovery procedure.

The vector α should be chosen in such a way that the recovery errors |η

| or

(η

)

be minimized. In turn, the squared error can be described by a quadratic

form, which must involve all of the linear forms occurred in the representation

(11.11) of η

We consider the solution that minimizes the squared recovery error, i.e.,

ˆα = arginf



(Λ − Q

∗

α)

∗

K (Λ − Q

∗

α)



where Q = ||P |I|| and I is the identity operator of the dimension equal to the

number of the measurement;

∗

α)

∗

= ||α

. . . α



(x)| . . . hp

µ1

(x)| . . . hp

(x)| 1 . . . 0

. . . . . . . . . . . . . . . . . . . . . 0

(x)| . . . hp

µn

(x)| . . . hp

(x)| 0 . . . 0

. . . . . . . . . . . . . . . . . . . . . 0

(x)| . . . hp

µN

(x)| . . . hp

(x)| 0 . . . 1



∗

= ||λ

∗

. . . λ

∗

. . . λ

∗

. . . λ

∗

M+1

. . . λ

∗

M+N

||,

where for i = 1, . . . , M, the λ

∗

are equal to h0| (except for λ

∗

= hl(x)|), whereas

for i = (M + 1), . . . , (M + N), the λ

∗

are zero. The operator K must be positive

deﬁnite.

In this case, the optimal estimation ˆα, which, by the linearity of the model

and the estimation of the functional, also determines the optimal estimation of the

functional l itself, is given by the expression

ˆα = (QKQ

∗

)

−1

QKΛ . (11.12)

Since the estimation error (11.11) depends on unknown ﬁelds ν

(δθ) and on the

random variable ε, the choice of a speciﬁc operator K in the criterion determines a

class of functions {ν

, µ = 1, . . . , M} and the statistical structure of the random

component ε, which is equivalent to the introduction of a priori information.

Consider the statistical interpretation of the operator K and the estimation ˆα.

If the set of ﬁelds {ν

(δθ)} accepted as admissible solutions is characterized by the

Gaussian measure in the corresponding functional space and the random component

332 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

is normally distributed, then the solution that minimizes the squared error of the

recovery E(η

) (E denotes the mean over the ﬁelds {ν

(δθ)} and realizations of ε):

K =



. . . K

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

. . . K

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

εν

. . . K

εν

. . . K

εν

. . . K

εε



Here K

is an operator of the integral type such that, for any linear functionals

ξ and η, the condition

hξ|K

|ηi = Ehξ|ν

ihν

|ηi,

where

εε

= E(εε

) , K

εν

: Eεhν

|ηi = K

εν

|ηi

is satisﬁed

The optimal estimation ˆα in (11.12) determines the structure of the recovery

algorithm with the minimal error variance. This estimation can be represented

in the following form, in which all statistical correlations of the ﬁelds ν

and ν

occur in the explicit form (under the assumption that the ﬁelds ν

and the random

component ε are not correlated, we have K

µε

= 0):

ˆα =



∗

6=µ

∗

+ P

∗

)

6=µ

∗

+ K

ε ε



−1



6=µ



|li.

In the particular case where there is no statistical correlations either between dis-

tinct ﬁeld ν

(δθ) and ν

(δθ) or between, ν

(δθ) and ε, the optimal estimation

(11.12) reduces to the form

ˆα



∗

6=µ

∗

+ K

ε ε



−1

|li,

where

∗



1µ

i . . . hp

1µ

Nµ

. . . . . . . . .

nµ

1µ

i . . . hp

nµ

Nµ

. . . . . . . . .

Nµ

1µ

i . . . hp

Nµ

