Troyan V., Kiselev Y. Statistical Methods of Geophysical Data Processing
Подождите немного. Документ загружается.
Construction of tomographic functionals 313
We obtain following representation for ˜ϕ:
˜ϕ
1
=
∂
t
˜ϕ
(3)
˜σ
(1)
˜σ
(2)
, ˜ϕ
2
=
˜σ
(3)
∂
t
˜ϕ
(1)
∂
t
˜ϕ
(2)
,
where ∂
t
ϕ
(i)
is the velocity of the particles in the direction of i-th ort; σ
i
= σ
3i
is
the component of the stress tensor; the tilde denotes the Fourier transform.
In this case
˜
L
1
=
1
λ+2µ
λS
1
λ+2µ
λS
2
λ+2µ
λS
1
λ+2µ
a −
µ(3λ+2µ)S
1
S
2
λ+2µ
λS
2
λ+2µ
−
µ(3λ+2µ)S
1
S
2
λ+2µ
b
,
where a = ρ −
4µ(λ+µ)S
2
1
λ+2µ
− µS
2
2
; b = ρ − µS
1
2
−
4µ(λ+µ)S
2
2
λ+2µ
;
˜
L
2
=
ρ S
1
S
2
S
1
µ
−1
0
S
2
0 µ
−1
.
Here S
i
= p
i
/ω (i = 1; 2) are the components of the slowness vector; β = −iω.
The equation (10.37) describes the plane P–SV-wave with the normal vector
n = e
1
, for the initial condition ˜σ
2
= ∂
t
˜ϕ
2
= 0, and the equation (10.37) describes
SH-wave, for the initial condition S
1
= 0, ˜σ
1
= ˜σ
3
= ∂
t
˜ϕ
1
= ∂
t
˜ϕ
3
= 0. Finally, if
a normal vector to the wave front is e
3
and initial conditions are ˜σ
1
= ˜σ
2
= ∂
t
˜ϕ
1
=
∂
t
˜ϕ
2
= 0, then the equation (10.37) describes P plane wave.
Let us consider an electromagnetic case. The Maxwell equations for an isotropic
medium reads as
∇ × E = −µ
∂H
∂t
,
∇
∇ × H = σE + ε
∂H
∂t
.
For the Fourier transform, by analogy with the relation (10.38), we obtain the
following representation
˜
L
1
,
˜
L
2
, ˜ϕ
1
and ˜ϕ
2
from the equation (10.37):
β = −iω , ˜ϕ
1
=
E
1
E
2
, ˜ϕ
2
=
−H
2
H
1
,
˜
L
1
=
µ −
S
2
1
γ
−
S
1
S
2
γ
−
S
1
S
2
γ
µ −
S
2
2
γ
,
˜
L
2
=
γ −
S
2
2
µ
S
1
S
2
µ
S
1
S
2
µ
γ −
S
2
1
µ
.
Where S
i
= p
i
/ω , i = 1; 2 and γ = σ/(ωi) + ε.
According to the boundary conditions the tangential components E and H must
be continuous at the crossing of the boundary. In the case of horizontally layered
medium it leads to continuity the vectors
˜
ϕ
1
and
˜
ϕ
2
. The components of the vector
˜
ϕϕ provide the following representation of the energy flow in x
3
direction
R
3
= −
1
4
( ˜ϕ
∗
1
˜ϕ
2
+ ˜ϕ
∗
2
˜ϕ
1
) = −
1
4
˜ϕ
∗
M ˜ϕ ,
314 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
M =
0 I
I 0
(I is identity matrix with dimension 2 × 2).
In the case of the acoustic medium, the vectors
˜
ϕ
1
,
˜
ϕϕ
2
and matrices
˜
L
1
,
˜
L
2
(from the equation (10.37)) are represented as
β = −iω , ˜ϕ
1
= ||∂
t
˜ϕ
3
||, ˜ϕ
2
= ||˜σ||,
˜
L
1
=
1
λ
−
S
2
ρ
,
˜
L
2
= ||ρ||, S
2
=
S
2
1
+ S
2
2
ω
2
.
The solution of the evolutionary equation (10.37) in the general form can be
represented with the help of the transition matrix P (x
3
, x
0
3
), which satisfies the
equation
∂
∂x
3
P (x
3
, x
0
3
) = β
˜
L(x
3
) P (x
3
, x
0
3
) (10.39)
and initial conditions P (x
0
3
, x
0
3
) = I.
Taking into account that the equation (10.39) can be represented in the equiv-
alent form:
P (x
3
, x
0
3
) = I + β
x
3
Z
x
0
3
dx
0
3
L(x
0
3
) P (x
0
3
, x
0
3
) ,
we can write the transition matrix in an iterative form:
P = I +
x
3
Z
x
0
3
˜
L(x
0
3
) dx
0
3
+
x
3
Z
x
0
3
dx
0
3
˜
L(x
0
3
)
x
0
3
Z
x
0
3
dx
00
3
˜
L(x
00
3
)
+
x
3
Z
x
0
3
dx
0
3
˜
L(x
0
3
)
x
0
3
Z
x
0
3
dx
00
3
˜
L(x
00
3
)
x
00
3
Z
x
0
3
dx
000
3
˜
L(x
000
3
) + . . .
Let note, that the transition matrix for the adjoint equation has an opposite
sign of β.
Chapter 11
Tomography methods of recovering the
image of medium
The interpretation of the reconstruction tomography as a problem of recovering
the image of an object from its ray projections has gradually transformed to its
interpretation as a problem of integral geometry, i.e. as the problem of recovering
a function θ(x) from the values of a given function
u(x) =
Z
Σ(x)
θ(x
0
) w(x, x
0
)dσ ,
where x ∈ R
n
; Σ(x) is a set of k-dimensional manifolds in R
k
, k < n; w(x, x
0
) is a
weight function; dσ is an element of a measure on Σ(x).
As we have seen in Chap. 10, where the notion of the tomographic functional was
introduced, any problem of recovering distributed parameters of a system (in partic-
ular, an arbitrary problem of remote sensing) can be considered as a problem of the
reconstruction tomography. Integrals of parametric functions defined on manifolds
of the dimension coinciding with that of the definition domain are treated as the
initial data of the corresponding tomographic problem. These manifolds are fixed
in the linear approximation; in the general case, they are implicitly transformed at
each iteration.
11.1 Elements of Linear Tomography
Mathematically, classical tomography is based on the Radon transform (Helgason,
1999; Herman, 1980):
u(p, n) =
ˆ
Rϕ =
Z
X
ϕ(x) delta(p − x · n)dx
∆
= hδ(p − x · n)|ϕ(x)i,
where u(p, n) is the Radon transform of a function ϕ(x), defined in an N -
dimensional space; δ(p −x ·n) is the singular kernel of the Radon transform, which
cuts out an (N −1)-dimensional hyperplane (a straight line in the two-dimensional
case), that is defined by the normal vector n and the aiming-distance parameter p,
i.e., the distance from the origin of coordinates to the hyperplane:
x · n = p.
315
316 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
The Radon transform arose in the classical tomography as a mathematical model
for the description multiforeshortening radiographies by X-rays of biological “soft”
objects, on which the scattering of X-rays may be assumed small. The transport
equation for a point (x
0
), monoenergetic (E
0
), collimation (ν
0
) source of X-rays
radiation is a reduced transport equation of the form
(νν ·∇∇)I(x) + ϕ(x)I(x) = cδ(x −x
0
)δ(ν −ν
0
)δ(E −E
0
). (11.1)
Here, I(x) is the beam of primary X-ray photons; ν
ν
ν is the direction of collimation of
the source; ϕ(x, E) is a linear coefficient; E is the energy of photons. It is assumed
that the vectors
ν
ν are coplanar, i.e. the problem considered is two-dimensional.
The solution of (11.1) can be written, upon performing the change of variables
x −→ (p, n), as
I(p, n, E
0
) = c
0
exp{−
Z
L(p,n)
ϕ(x) dx}, (11.2)
where L(p, n) is the line which radiation propagates. Taking the logarithm of the
relation (11.2), we obtain the Radon transform for the two-dimensional problem
−ln
I(p, n, E
0
)
c
∆
= u(p, n) =
Z
L
ϕ(x) dx =
Z
X
ϕ(x)δ(p − x ·n) dx .
Here, in the rightmost expression, the integration is carried out over the entire
domain of the definition of ϕ(x).
The model of a planar cross section of an object arises in describing plasma
sounding by a planar laser ray (“light knife”). After its transition through the
volume of plasma under investigation, the planar beam is focused on a photodiode.
As it is not difficult to see, in both cases tomographic experiments are based on the
adequate use of the ray description of the propagation of the sounding signal (the
wavelength is small in comparison with the characteristic size of an inhomogeneity
with a strongly stretched forward scattering indicatrix).
Before passing to the basic properties of the direct Radon transform, we recall
some properties of linear continuous functionals, i.e., of generalized functions
l(ϕ)
∆
= (l, ϕ) = (l(x), ϕ(x)),
where x ∈ R
n
.
11.1.1 Change of variables
Let x = Ly + x
0
be a nonsingular mapping, then
[l(Ly + x
0
), ϕ(y)] = |L|
−1
[l(x), ϕ(L
−1
(x − x
0
))] ,
whence for the delta-function we obtain
hδ(Ly + x
0
)|ϕ(y)i
∆
= hδ(x)||L|
−1
ϕ(L
−1
(x − x
0
))i.
In particular,
Tomography methods of recovering the image of medium 317
(a) hδ(x − x
0
)|ϕi = ϕ(x
0
) ;
(b) hδ(ay)|ϕ(y)i = hδ(x)||a|
−n
ϕ(a
−1
x)i;
(c) if L : L
∗
L = I and x
0
= 0, then
hδ(Ly)|ϕ(y)i = hδ(x)|ϕ(L
∗
x)i;
(d) hδ(L(x))ϕ(x)i =
P
k
hδ(y)|| L
0
|
x
k
|
−1
ϕ(L
0
|
−1
x
k
y + x
k
)i
=
P
k
| L
0
|
x
k
|
−1
ϕ(x
k
),
where x
k
: L(x
k
) = 0.
For example:
δ(x
2
− c
2
t
2
) =
1
2ct
[δ(x − ct) + δ(x + ct)] ;
δ(sin x) =
∞
X
k=−∞
δ(x − kπ) ;
δ(x − x
0
) =
Y
i
|h
i
|
−1
Y
i
δ(q
i
− q
0
i
) .
Here, q
i
= q
i
(x), and h
i
=
P
k
(∂x
k
/∂q
i
)
2
, which corresponds to the passage from
Cartesian coordinates to curvilinear orthogonal coordinates q
i
(x).
11.1.2 Differentiation of generalized function
Let l ∈ C
n
and let |α| ≤ n. Then
(D
α
l, ϕ) = (−1)
α
(l, D
α
ϕ) .
In particular,
hD
α
δ|ϕi
∆
= (−1)
|α|
D
α
ϕ(0) .
Below we list the basic properties of the Radon transform (Pikalov and Preo-
brazhenski, 1983; Helgason, 1999).
1. Its linearity
ˆ
R[α
1
ϕ
1
+ α
2
ϕ
2
] = α
1
ˆ
R[ϕ
1
] + α
2
ˆ
R[ϕ
2
]
follows from the definition.
2. In view of (11.2), the rotation of the object leads to the same rotation of the
vector n in the argument of the Radon transform, i.e.,
ˆ
R[ϕ(Ox)] = u(p, On); (11.3)
the similarity transformation with coefficient α 6= 0 increases the aiming distance
by the factor α and decreases the Radon transform by the factor α
1−n
, i.e.,
R[ϕ(αx)] = α
−n
u
p,
n
α
= α
1−n
u(αp, n) .
318 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
3. The shift of the object by x
0
amounts to the shift of the aiming distance by
n · x
0
, i.e.,
ˆ
R[ϕ(x − x
0
)] = u(p − n · x
0
, n) .
4. The Radon transform of the derivative of the function ϕ, which describes the
object, is given by
ˆ
R[(k · ∇)ϕ] = k · n
∂u(p, n)
∂p
,
for the mixed second derivative we have
ˆ
R[(k ·∇∇)(l · ∇)ϕ] = (k · n)(l · n)
∂
2
u(p, n)
∂p
2
,
whence it follows that
ˆ
R[(∇ ·∇∇)ϕ(x)] = |n|
2
∂
2
u
∂p
2
=
∂
2
u(p, n)
∂p
2
. (11.4)
Finally, for an arbitrary linear differential operator
ˆ
L = L(∂/∂x
1
, ∂/∂x
2
, . . . ,
∂/∂x
k
) with constant coefficients, we have
ˆ
R[
ˆ
Lϕ(x)] =
ˆ
L
n
1
∂
∂p
, . . . , n
k
∂
∂p
u(p, n) .
An example of using property (11.4) is provided by the application of the direct
Radon transform to the 3-dimensional wave equation
∇
2
−
1
c
2
∂
2
∂t
2
ϕ(x, t) = 0 ,
ˆ
R
∇
2
−
1
c
2
∂
2
∂t
2
ϕ
=
∂
2
∂p
2
−
1
c
2
∂
2
∂t
2
u(p, n) ,
which makes possible to reduce the 3-dimensional wave equation to the correspond-
ing 1-dimensional equation.
5. Differentiation the Radon transform
(k · ∇)
ˆ
R[ϕ(x)] = −
∂
∂p
ˆ
R[(k · x)ϕ(x)]
for two vectors k and l and the derivatives with respect to components of the unit
vector n, we arrive at the equality
(k · ∇)(l · ∇)u(p, n) =
∂
2
∂p
2
ˆ
R[(k · x)(l · x)ϕ(x)] .
6. The Radon transform of the convolution of functions is as follows:
ˆ
R
Z
g(x − x
0
)s(x
0
)dx
0
=
Z
s(x
0
)dx
0
Z
g(x − x
0
)δ(p − n · x)dx.
Tomography methods of recovering the image of medium 319
Performing the change of variables x − x
0
= y, we derive
u(p, n) =
Z
s(x
0
) dx
0
Z
g(y)δ(p − n · x
0
− n · y) dy
=
Z
s(x
0
)
ˆ
R[g(p − n · y, n)] dx
0
=
Z
ˆ
R[g(p − p
0
, n)] dp
0
Z
s(x
0
)δ[p
0
− n · x
0
] dx
0
,
i.e.
ˆ
R[g ∗ s] =
ˆ
R[g] ∗
ˆ
R[s].
Note that the Radon transform of a multidimensional convolution is equal to the
one-dimensional convolution of the Radon transforms with respect to the aiming
distance p
∗
.
7. The relation between the Radon and Fourier transforms.
We write the Fourier transform with respect to the aiming distance of the Radon
transform u(p, n):
˜u(p, n) = (2π)
−1/2
Z
e
−ipν
dp
Z
ϕ(x) δ(p − x · n) dx
= (2π)
−1/2
Z
dx ϕ(x)
Z
e
−ipν
δ(p − x · n) dp
= (2π)
−1/2
Z
ϕ(x) e
−iν(x·n)
dx,
or, in operator form,
F
p
[u(p, n)]
∆
= F
p
[
ˆ
R[ϕ(x)]] = (2π)
(n−1)/2
F
x
[ϕ(x)] . (11.5)
The one-dimensional Fourier of the Radon transform is equivalent to the multi-
dimensional Fourier transform of the function ϕ(x), which characterizes the object.
This relation between the Fourier transform of the Radon transform and multidi-
mensional Fourier transform of ϕ is formulated as the following generalized theorem
on the central section: for a fixed n, the one-dimensional Fourier transform of the
Radon transform yields the central section of the multidimensional Fourier spectrum
along the ray nν, passing through the origin of coordinates.
Based on the relation (11.5), the inverse Radon transform is derived in the
following way:
ϕ(x) = (2π)
(−n+1)/2
F
−1
x
F
p
u(p, n)
= C
n
Z
d(νn) e
iνn·x
Z
dp e
−iνp
u(p, n)
= C
n
Z
d(νn)
Z
dp e
iν(q−p)
u(p, n),
where q = n · x. Taking into account that d(νn) = dnν
n−1
dν, we represent ϕ(x)
in the form
ϕ(x) = C
n
Z
|n|=1
dn
∞
Z
0
dν ν
n−1
Z
dp e
iν(q−p)
u(p, n) = C
n
Z
|n|=1
dn J(n, q) .
320 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
We extend the limits of integration with respect to ν to −∞ and +∞, either by
introducing an even function J
e
(n, q),
J
e
(n, q) =
1
2
[J(n, q) + J(−n, −q)]
or by using the following explicit representation of the integral J(n, q):
J
e
(n, q) =
1
2
∞
Z
0
ν
n−1
dν
∞
Z
−∞
u(p, n) e
iν(q−p)
dp
+
1
2
∞
Z
0
ν
n−1
dν
∞
Z
−∞
u(p, −n) e
−iν(q−p)
dp
=
1
2
∞
Z
−∞
|ν|
n−1
dν
∞
Z
−∞
u(p, n) e
iν(q−p)
dp .
The cases of even and odd values of n are considered separately.
1. Let n be odd. Integrating J
o
(n, q) by parts, n −1 times, we obtain
J
o
(n, q) =
1
2
(i)
1−n
∞
Z
−∞
dν
∞
Z
−∞
dp
∂
∂p
n−1
u(p, n) e
iν(q−p)
=
1
2
(i)
1−n
∞
Z
−∞
∂
∂p
n−1
u(p, n) δ(q − p) dp
=
1
2
(i)
1−n
∂
∂p
n−1
p=q=n·x
u(p, n) .
Thus, for odd values of n, the calculation of the integral J
o
(n, q) is reduced
to the computation of the derivative of order (n − 1) with respect to the aiming
distance p of the Radon transform u(p, n) at the point p = g = n · x. The ultimate
formula for Radon inversion in this case is
ϕ(x) = C
n
Z
|n|=1
∂
∂p
n−1
p=n·x
u(p, n) dn , (11.6)
where C
n
= (2πi)
1−n
/2.
2. Let n be even. Similarly to the case of add values of n, we integrate J
e
(n, g)
by parts n − 1 times, which yields
J
e
(n, q) =
1
2
(i)
1−n
∞
Z
−∞
sgn(ν)dν
Z
dpe
−iν(p−q)
∂
∂p
n−1
u(p, n)
=
1
2
(i)
1−n
∞
Z
−∞
dp
∂
∂p
n−1
u(p, n)
∞
Z
−∞
dνsgn(ν)e
−iν(p−q)
=
(i)
−n
2π
Z
(∂/∂p)
n−1
u(p, n)
p − q
dp .
Tomography methods of recovering the image of medium 321
Here we have used the relation
F [sgn(ν)] =
i
π
P
1
p
,
where P is the symbol of integration in the sense of Cauchy principal value. De-
noting the Hilbert transform by H, i.e., setting
H(u) = −
1
π
P
1
p
∗ u(p)
,
we can write the following representation of the Radon inversion in an even-
dimensional space:
ϕ(x) = −iC
n
Z
|n|=1
dn
H
∂
∂p
n−1
u(p, n)
(n · x, n)
= −
iC
n
π
Z
|n|=1
dn P
∞
Z
−∞
(∂/∂p)
n−1
u(p, n)
p − n · x
dp . (11.7)
Thus, in the case of an even value of n, in order to recover the function ϕ(x) one
needs to compute the derivative of the order (n − 1) with respect to the aiming
distance p of the Radon transform u(p, n), thereafter perform the Hilbert transform
with respect to the aiming distance p, and, finally, carry out integration over the
unit sphere.
Analyzing representations (11.6) and (11.7) of the Radon inversion for even- and
odd-dimensional spaces, we readily see that the Radon inversion is “local” for odd
values of n. In this case, Radon surfaces on hyperplanes passing in the vicinity of
the space point x of the recovery are used. In even-dimensional spaces, the Radon
inversion is essentially nonlocal, which is explicitly demonstrated by the integral
Hilbert transform. The manifestation of the locality or nonlocality of the Radon
inversion is caused by a close relation between the Radon inversion and harmonic
functions. In particular, this relation manifests itself in the nonlocality of wave
fronts in even-dimensional spaces, where a wave front is followed by a diffusion
train, which is in contrast to sharply pronounced wave fronts in odd-dimensional
spaces.
We want to represent the Radon transform in the operator form. To this end,
for an arbitrary function ψ(p, n) (where p = n ·x), such that ψ(p, n) = ψ(−p, −n),
we introduce the conjugate Radon operator R
∗
:
R
∗
ψ =
Z
|n|=1
ψ(n · x, n)dn .
Then, for odd values of n, the Radon inversion can be represented in the operator
form as
R
−
= R
∗
Γ
e
,
322 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
where
Γ
e
= C
n
∂
∂p
n−1
p=n·x
,
whereas for even values of n, it can be represented in the form
R
−
= R
∗
Γ
a
, Γ
a
= −iC
n
H
∂
∂p
n−1
.
The operators Γ
o
and Γ
e
are counterparts of the operator (RR
∗
)
−
, which occurs
in the following representation of the solution of the integral equation Rϕ = u
provided by the least squares method:
ϕ = R
∗
(RR
∗
)
−
u.
The subscript in the formula for the operator R
−
of the Radon inversion is used
because Radon inversion is an unstable procedure, and thus a regularization must
be used (Pikalov and Preobrazhenski, 1983; Tikhonov et al., 1987).
11.2 Connection of Radon Inversion with Diffraction Tomography
The linearized model of measurements in problems of sounding an object by a signal
ϕ with the equation of propagation Lϕ = s is represented by the relation (10.9):
˜u
n
= h(L
−1
0
)
∗
h
n
|δL
θ
ϕ
0
i + ˜ε
n
∆
= hϕ
n (out)
|δL
θ
|ϕ
in
i + ˜ε
n
.
The idealized analog of this model is obtained by representing the apparatus
function (h
n
) by a point receiver with the infinite band-pass, i.e., by letting
h
n
⇒ δ(x − x
n
)δ(t − t
n
) and by neglecting noises: ˜ε
n
⇒ 0. Then ˜u
n
will be
represented as the scattering field at the measurement point x
n
, i.e.,
˜u
n
= ˜ϕ|
x
n
∆
= (ϕ − ϕ
in
)|
x
n
.
We write this model in the form ˜ϕ = hϕ
out
|δL
θ
|ϕ
in
i.
As an example, consider the scalar wave equation (c
−2
∂
2
/∂t
2
−∆)ϕ = 0 in the
case of a homogeneous reference medium, i.e., in the case c
0
(x) = c
0
= const. As it
was shown in Sec. 10.3, the tomographic functional p for the scalar wave equation
is of the form
p : ˜ϕ = hp|ν(x)i,
where ν(x) = c
−2
0
(1 − c
2
0
/c
2
(x)), and the integral kernel of the tomographic func-
tional is given by the expression
p = hϕ
out
∂
2
∂t
2
ϕ
in
i
T
.