Inversion Theory
(Lecture 6)
MICHAEL ZHDANOV
1
Formulation of well-posed and ill-posed problems
In the first lecture we have formulated an inverse problem as the solution of
an operator equation
d = A (m) ,
(1)
where m ∈ M, is some function (or vector) from a metric space M of the
model parameters, and d ∈ D is an element from a metric space D of data
sets. There are two important classes of inverse problems: well-posed and
ill-posed problems. We will give detailed description of these problems in
this section.
1.1
Well-posed problems
Following classical principles of regularization theory we can give the following definition of the well-posed problem.
Definition 1 The problem (1) is correctly ( or well) posed if the following
conditions take place: (i) solution m of the equation (1) exists, (ii) solution
m of the equation (1) is unique, (iii) solution m continuously depends on the
right-hand side d.
In another words the inverse operator A−1 is defined throughout the space
D and is continuous.
1
Note, that well-posed inverse problem possesses all the properties of the
”good” solution discussed in the previous lectures: the solution exists, is
unique and stable.
Definition 2 The problem (1) is ill-posed if at least one of the conditions,
listed above, fails.
We have demonstrated that the majority of geophysical inverse problems
are ill-posed, because at least one of the conditions listed above fails. However, it may happen that if we narrow the class of the models which are used
in inversion, the originally ill-posed inverse problem may become well-posed.
Mathematically it means that instead of considering m from entire model
space M, we can select m from some subspace of M, consisting of more simple and/or more suitable models for the given inverse problem. Thus, we
arrive at the idea of the correctness set and conditionally well-posed inverse
problem.
1.2
Conditionally well-posed problems
Suppose that we know a priori that the exact solution belongs to a set C
such, that the inverse operator A−1 defined on the image AC is continuous.
We call C the correctness set.
Definition 3 The problem (1) is conditionally well-posed (Tikhonov’s wellposed) if the following conditions are met: (i) we know a priori that a solution
of (1) exists and belongs to a specified set C ⊂ M , (ii) the operator A is
one-to-one mapping of C onto AC ⊂ D, (iii) the operator A−1 is continuous
on AC ⊂ D.
In contrast to the standard well-posed problems, a conditionally wellposed problem doesn’t require solvability over the whole space. Also the
requirement of continuity of A−1 over all M is substituted by the requirement
of continuity over the image of C. Thus, introducing a correctness set makes
even ill-posed problem to be well-posed, if we reduce the class of possible
solutions.
Tikhonov described the correctness set C. For example , if the model is
described by a finite number of bounded parameters, they form a correctness
set C in Euclidean space of model parameters.
2
1.3
Quasi-solution of ill-posed problem
We assume now that the problem (1) is conditionally well-posed (Tikhonov’s
well-posed). Suppose also that the right-hand side of (1) is given with some
errors:
dδ = d+δd,
(2)
where
µD (dδ , d) ≤ δ
(3)
Definition 4 A quasi-solution of the problem (1) on the correctness set C
is an element mδ ∈ C which minimizes the distance µD (Am, dδ ), i.e.:
µD (Amδ , dδ ) =
inf µD (Am, dδ ).
m∈CµD
(4)
It can be proved also, that the quasi-solution is a continuous function of
dδ .
Figure 1 illustrates the definition of a quasi-solution. Element m ∈ M is
exact solution of the inverse problem
d = A (m) .
(5)
Subset AC of the data space D is an image of the correctness set C as the
result of operator A applications. A quasi-solution, mδ , is selected from a
correctness set C under the condition that its image, A (mδ ) , is the closest
element to the observed noisy data, dδ ,from subset AC.
The idea of quasi-solution makes it possible to substitute the inverse problem solution by minimization of the distance µD (Am, dδ ) on some appropriate class of suitable models. The standard methods of functional minimization can be used to solve this problem and therefore to find the quasi-solution.
In this way, we significantly simplify the inverse problem solution.
2
2.1
The regularization method in the solution
of the inverse problem
Regularizing operators
Consider inverse geophysical problem described by the operator equation:
d = A(m),
3
(6)
m = A-1(d)
D
M
~
m
m
AC
mδ
δ
Am δ
C
d = A(m)
d
δd
d=d+
δ
m δ = A-1 (d δ )
Figure 1: A quasi-solution, mδ , is selected from the correctness set C under
the condition that its image, A (mδ ) , is the closest element to the observed
noisy data, dδ , from the subset AC : µD (Amδ , dδ ) = inf m∈C µD (Am, dδ )
4
where m presents model parameters and d is observed geophysical data.
In general case the inverse operator A−1 is not continuous and therefore
the inverse problem (6) is ill-posed. The main idea of any regularization
algorithm is to consider instead of one ill-posed inverse problem (6) a family
of the well-posed problems
d = Aα (m),
(7)
which approximate in some sense the original inverse problem. Scalar parameter α > 0 is called a regularization parameter. We require also that
mα → mt , if α → 0,
where mα = A−1
α (d) is the solution of the inverse problems (7), and mt is the
true solution of the original problem (6). Thus, we substitute the solution
of one ill-posed inverse problem by the solutions of the family of well-posed
problems, assuming that these solutions, mα , asymptotically go to the true
solution if α goes to zero.
Let us give now the more accurate definitions.
Definition 5 Operator R(d, α) (depending on a scalar parameter α) is called
the regularizing operator in some vicinity of the element dt = A(mt ) if there
is a function α(δ) such that for any > 0 it can be found a positive number
δ( ) such that if
µD (d, dt ) < δ( ),
then
µM (mα , mt ) < ,
where
mα = R(d, α(δ)).
In other words mα is a continuous function of the data and
mα = R(d, α(δ))→ mt
(8)
when α → 0.
Figure 2 illustrates the basic properties of the regularizing operators. Let
mt is exact solution for exact data dt = A(mt ). However, we can observe
only the noisy data dδ = dt + δd. If we apply some rigorous inverse operator
to noisy data, we could get a result mδ which lies far away from true solution.
5
h one can get another completely different from
By perturbing noisy data d
δ
i δ . The main advantage of the regularizing operators
the true solution result m
R is that they provide the stable solution in any situation. If we apply R
to the noisy data dδ , we will get a solution mδ = R(dδ , α) which is very
h
close to a true model mδ −mt < . Application of R to a noisy data d
δ
h
i δ = R(dδ , α), which is still close to mt . The
will result in another solution m
accuracy of the true solution approximation by the regularized one depends
on the regularization parameter α. The smaller the α the more accurate is
the approximation.
We can see that regularizing operators can be constructed by approximating ill-posed equation (6) by the system of well-posed equations (7), where
corresponding inverse operators A−1
α are continuous. These inverse operators
can be treated as regularizing operators
A−1
α (d) = R(d, α).
The only problem now is how to find the family of the regularizing operators. Tikhonov suggested the following scheme for constructing regularizing
operators. It is based on introducing special stabilizing and parametric functionals.
2.2
Stabilizing functional
Stabilizing functional (or stabilizer ) is used to select from the space M of
all possible models the subset Mc which is a correctness set.
Definition 6 A nonnegative functional s(m) on some metric space M is
called a stabilizing functional if for any real number c > 0 the subset Mc of
elements m ∈ M for which s(m) ≤ c is a correctness set.
We will give now an example of the stabilizing functionals.
Example 7 Consider the Hilbert space L2 formed by the functions integrable
on the interval [a, b]. The metric in the space L2 is determined according to
formula:
µ(m1 , m2 ) =
+]
b
a
2
,1
[m1 (x) − m2 (x)] dx
6
2
,
(9)
m'δ = A-1 (d δ )
D
M
m~' δ
m' δ
~
mδ
mt
mδ
~
~
m'δ = A-1 (d δ ) ~
d
~
~
m δ = R (d δ ,α)
dt
δ
d δ = d t + dδ
m δ = R (d δ ,α)
Figure 2: The scheme illustrating the construction of the regularizing operators. The bold point mt denotes the true solution, and the point dt ∈ D
denotes the true data. The noisy data are shown by a point dδ = dt + δd.
Application of a formal inverse operator A−1 to the noisy data generates
i δ , which are unstable with respect to a small
formal solutions, mδ and m
h . However, application of the regularizing
perturbation in the data, dδ or d
δ
h , produces stable
operators, R(d, α) to any of the observed data, dδ or d
δ
i δ are close to each other and to the
results: the inverse models mδ and m
h , are close to each other and to
true solution, if the observed data, dδ and d
δ
the true data dt .
7
It can be proved that any ball
b (m0 ,c) = {m : µ(m, m0 ) ≤ c }
is a correctness set in a Hilbert space. Therefore, we can introduce the stabilizing functional as following:
s(m) = µ(m, m0 ),
(10)
where m0 is any given model from M = L2 and a ball :
s(m) = µ(m, m0 ) ≤ c,
(11)
is a correctness set.
Let us analyze now more carefully how one can use the stabilizers to select
an appropriate class of the models. Assume that the data dδ are observed
with some noise dδ = dt +δd, where dt is the true solution of the problem. In
other words, we assume that the misfit (distance) between the observed data
and true data is less than the given level of the errors in the observed data,
equal to δ ≥ δd :
µD (dδ , dt ) ≤ δ.
(12)
In this situation it is naturally to search for an approximate solution in the
set Qδ of the models m such that
µD (A(m), dδ ) ≤ δ.
(13)
Thus Qδ ⊂ M is the set of possible solutions.
Figure 3 helps to understand this role of the stabilizing functional. The
main application of stabilizers is to select from the set of possible solutions
Qδ the solutions, which continuously depend on the data and which possesses
a specific property depending on the choice of a stabilizer. Such solutions can
be selected by the condition of the minimum of the stabilizing functional:
s(m; m ∈ Qδ ) = min.
(14)
We have introduced a stabilizing functional under the condition that it selects
a correctness subset MC from a metric space of the model parameters. Thus,
8
d = A (m)
D
M
Am δ
Qδ
δ
dδ
mδ
dt
MC
m δ = R (d
δ ,δ)
Q δ = { µD (Am, dδ)
s (m δ ; m δ
< δ}
dδ = d + δd
Qδ ) = min
Figure 3: The stabilizing functional selects from a set of the possible solutions, Qδ , a solution, mδ , which at the same time belongs to the correctness
set MC .
9
we can say that a stabilizer selects from the set of possible solutions Qδ a
solution, which at the same time belongs to a correctness set MC .
The existence of the model, minimizing (14), was demonstrated in (Tikhonov
and Arsenin, 1977). We will denote this model as mδ :
s(mδ ; mδ ∈ Qδ ) = min.
(15)
One can consider a model mδ as the result of application to the observed
data dδ an operator R(dδ , δ), depending on the parameter δ:
mδ = R(dδ , δ)
(16)
It can be proved that the operator R(dδ , δ) is the regularizing operator
for the equation (??) and mδ can be used as the approximate solution of the
inverse problem ( note that in this case α = δ, while in general case α = α(δ)).
Thus in the framework of developed approach the problem of the solution
of eq. (1) with approximate left hand part dδ can be reduced to the problem
of minimization of stabilizing functional on the set Qδ :
s(m; m ∈ Qδ ) = min,
where
Qδ = {m; µD (A(m), dδ ) ≤ δ}.
2.3
(17)
Tikhonov parametric functional
It can be proved that for a wide class of stabilizing functionals they reach
their minimum on the model mδ such that µD (A(mδ ), dδ ) = δ. Thus , we
can solve the problem of minimization (14) under the condition that
µD (A(mδ ), dδ ) = δ.
(18)
The last problem is the problem of the stabilizing functional (14) minimization when the model m is subject to the constrain (18). A common way to solve this problem is to introduce an unconstrained functional
P α (m, dδ ), m ∈ M, given by
P α (m, dδ ) = µ2D (A(m), dδ ) + αs(m),
10
(19)
and to solve the problem of minimization of this functional:
P α (m, dδ ) = min.
(20)
Here the unknown real parameter α is similar to the Lagrangian multiplier
and it is determined under the condition
ρD (A(mα ),dδ ) = δ,
(21)
where mα is the element on which P α (m, dδ ) reaches minimum. The functional P α (m, dδ ) is called Tikhonov parametric functional.
Thus for any positive number α > 0 and for any data dδ ∈ D we have
determined an operator R(dδ , α) with the values in M such that the model
mα = R(dδ , α)
(22)
gives the minimum of the Tikhonov parametric functional P α (m, dδ ).
The fundamental result of the regularization theory is that this operator
R(dδ , α) is the regularizing operator for the problem (1).Thus as an approximate solution of the inverse problem (1) we take the solution of another
problem (20) (problem of minimization of Tikhonov parametric functional
P α (m, dδ )) , closed to initial problem for the small values of the errors δ of
the data dδ .
It is important to underline that in the case when A is a linear operator, D
and M are Hilbert spaces and s(m) is a quadratic functional , the solution of
the minimization problem (??) is unique. Note that the quadratic functional
is a functional q(m) such, that
q(βm) = β 2 q(m).
2.4
Determination of the regularization parameter
One of the most critical problems of regularization method is the selection
of the optimum value of the regularization parameter α. The solution of this
problem can be based on the following consideration.
Suppose that the data dδ is observed with some noise dδ = dt + δd,
where dt is the true solution of the problem, and the level of the errors in
the observed data is equal to δ:
µD (dδ , dt ) ≤ δ.
11
(23)
Then the regularization parameter can be determined by the misfit condition
(22)
µD (A(mα ), dδ ) = δ.
(24)
To justify this approach we will examine more carefully the properties of
all three functionals involved in the regularization method: the Tikhonov
parametric functional, the stabilizing and misfit functionals.
Let us introduce the following notations:
p(α) = P α (mα , dδ );
i(α) = µD (A(mα ), dδ );
(25)
s(α) = s(mα ).
Examine some properties of the functions p(α), i(α), s(α).
Property 1
Functions p(α), i(α), s(α) are monotone functions: p(α) and i(α) are not
decreasing and s(α) is not increasing.
Property 2
It can be proved that the functions p(α), i(α), s(α) are continuous functions (if the element mα is unique).
Note also that
p(α) → 0 f or α → 0,
and
p(0) = 0.
(26)
From the fact that
i(α) + αs(α) = p(α) → 0 f or α → 0,
it follows that
i(0) = 0.
(27)
Thus one can prove the following theorem.
Theorem 8 If i(α) is one-to-one function then for any positive number δ <
δ0 = µD (A(m0 ), dδ ) (where m0 is some a priori model) there exists α(δ) such
that µD (A(mα(δ) ), dδ ) = δ
12
p, i, s
p(α)
i(α)
δ
2
s(α)
αo
α1
α
Figure 4: Illustration of the principle of optimal regularization parameter
selection
Note that i(α) is one-to-one function when element mα is unique. It
happens, for example, when A is a linear operator, D is a Hilbert space and
s(m) is a quadratic functional.
Figure 4 helps to understand the principle of optimal regularization parameters selection. One can see that due to monotone character of the function
s(α), there is only one point, α0 , where s(α0 ) = δ. Let us consider one simple
numerical method of determination of parameter α. Consider for example
the progression of numbers:
αk = α0 q k , k = 1, 2, ..., n; q > 0.
(28)
For any number αk we can find element mαk , minimizing P αk (m, dδ ) and
calculate the misfit µD (A(mαk ), dδ ). The optimal value of the parameter α
is the number αk0 , for which with the necessary accuracy we have the equality
µD (A(mαk0 ), dδ ) = δ
The equality (29) is called the misfit condition.
13
(29)
2.5
L-curve method of regularization parameter selection
L-curve analysis (Hansen, 1998) represents a simple graphical tool for qualitative selection of the quasi-optimal regularization parameter. It is based
on plotting for all possible α the curve of the misfit functional, i(α),versus
the stabilizing functional, s(α) (where we use notations (25)). The L-curve
illustrates the trade-off between the best fitting (minimizing a misfit) and
most reasonable stabilization (minimizing a stabilizer). In a case where α
is selected to be too small, the minimization of the parametric functional
P α (m) is equivalent to the minimization of the misfit functional; therefore
i(α) decreases, while s(α) increases. When α is too large, the minimization
of the parametric functional P α (m) is equivalent to the minimization of the
stabilizing functional; therefore s(α) decreases, while i(α) increases. As a
result, it turns out that the L-curve, when it is plotted in log-log scale, very
often has a characteristic L-shape appearance (Figure 5), that justifies its
name (Hansen, 1998).
The distinct corner, separating the vertical and the horizontal branches
of this curve, corresponds to the quasi-optimal value of the regularization
parameter α.
14
log s(α)
α quasi-optimal
log i(α)
Figure 5: L-curve represents a simple curve for all possible α of the misfit
functional, i(α), versus stabilizing functional, s(α), plotted in log-log scale.
The distinct corner, separating the vertical and the horizontal branches of this
curve, corresponds to the quasi-optimal value of the regularization parameter
α.
15
Inversion Theory
Assignment 4
Due Date : 10/07/2020
1. Contraction operator
An operator 𝐶 of 𝑋 into itself is called a
contraction operator if there exists a positive
real number 𝑟 < 1 with the property that
𝜇(𝐶(𝑥), 𝐶(𝑥0)) ≤ 𝑟𝜇(𝑥, 𝑥0) for all points 𝑥 and
𝑥0 in 𝑋.
Prove that the contraction operator is a
continuous operator.
2. Bounded operator
Linear operator 𝒚 = 𝐿(𝒙) is called
bounded if there exists a real number M with
the property that
𝐿(𝑥) ≤ 𝑀 𝑥
for every 𝑥 ∈ 𝑋.
Prove that the linear operator 𝐿 is
continuous if it is bounded.
(1)
3. Gram - Schmidt orthogonalization process
It is easier solving a system of linear equation if use an
orthonormal set of basis elements. The set of elements
{𝒅1, 𝒅2, 𝒅3, ⋯ , 𝒅𝑛} can be converted into a corresponding
orthonormal set 𝒆1, 𝒆2, 𝒆3, ⋯ , 𝒆𝑛 as follws.
The first step is
𝒅1
𝒆1 =
𝒅1
(2)
The 𝑒2 and 𝑒3 are calculated as follows:
𝒅2 − 𝒅2, 𝒆1 𝒆1
𝒆2 =
𝒅2 − 𝒅2, 𝒆1 𝒆1
(3)
𝒅3 − 𝒅3, 𝒆2 𝒆2 − 𝒅3, 𝒆1 𝒆1
𝒆3 =
𝒅3 − 𝒅3, 𝒆2 𝒆2 − 𝒅3, 𝒆1 𝒆1
The tasks: a) proof that 𝒆1, 𝒆2, 𝒆3 form an
orthonormal set of elements
b) present similar expressions for 𝒆4
and 𝒆5.
(4)
4. 2-D Gravity Inverse problem
Let us consider a 2-D gravity problem. The
anomalous gravity field 𝑔𝑧 is generated by some mass
distributed in rectangular domain S within the lower
half plane with the unknown excess density Δ𝜌(𝑥, 𝑧).
The gravity field observed on the earth’s surface is
given (same as in HW3).
Problem:
Solve the inverse gravity problem of determining
the excess density distribution, Δ𝜌(𝑥, 𝑧), using Riesz
representation theorem.
Directions:
The anomalous 2-D gravity field can be calculated by the
formula:
𝜌 𝑥, 𝑧 ∙ 𝑧
′
(5)
𝑔𝑧 𝑥 , 0 = 2𝛾
𝑑𝑧𝑑𝑥
𝑥 − 𝑥′ 2 + 𝑧2
𝑆
We have observations at 41 points, which can be written
as follows:
𝜌 𝑥, 𝑧 ∙ 𝑧
𝑑𝑗 = 𝑔𝑧 𝑥𝑗 , 0 = 2𝛾
𝑑𝑧𝑑𝑥 =
2
2
+𝑧
𝑆𝑖 𝑥 − 𝑥𝑗
(6)
𝑧
=
𝜌 𝑥, 𝑧 2𝛾
𝑑𝑧𝑑𝑥
2
2
𝑥 − 𝑥𝑗 + 𝑧
𝑆
According to Riesz representation theorem, there exist the
vectors 𝑙 𝑗 (the elements of the model space 𝑀), which can be used to
represent the observed data:
𝑑𝑗 = 𝜌 𝑥, 𝑧 , 𝑙 𝑗
(7)
𝑗 = 1, 2, ⋯ 𝑛; 𝑙 𝑗 ∈ 𝑀.
Vectors 𝑙 𝑗 are called “the data kernels”. Comparing (7) and (6),
we see that
𝑧
𝑗
𝑗
𝑙 = 𝑙 𝑥, 𝑧 = 2𝛾
𝑥 − 𝑥𝑗 2 + 𝑧2
(8)
The unknown distribution of the density, according to the Riesz
representation theorem, can be expressed as
𝑛
𝑛
𝑧
𝑗
𝜌 𝑥, 𝑧 =
𝛽𝑗𝑙 = 2𝛾
𝛽𝑗
(9)
𝑥 − 𝑥𝑗 2 + 𝑧2
𝑗=1
𝑗=1
where the coefficients 𝛽𝑗 satisfy the following system of linear
equations:
𝑛
𝑑𝑖 = 𝑔𝑧 𝑥𝑗, 0 =
Γ𝑖𝑗𝛽𝑗
𝑗=1
(10)
We can calculate the Gram matrix as follows:
𝑖
Γ𝑗𝑖 = 𝑙 , 𝑙
𝑗
=
𝑧2
4𝛾2
𝑆
𝑥 − 𝑥𝑖
2
+
𝑧2
𝑥 − 𝑥𝑗
𝑗 = 1, 2, ⋯ , 𝑛; 𝑖 = 1, 2, ⋯ , 𝑛.
2
+
𝑧2
𝑑𝑧𝑑𝑥
(11)
You are recommended to use Matlab for this assignment.
Assume that resulting density image consists of a number of
small enough pixels. Take 100 cells in x direction and 50 cells in z
direction. The image has to occupy a rectangular subsurface
region (0,1000) along the X axis and (25,225) along the Z axis.
Note that in this approach, we approximate one data kernel
𝑙 (𝑗) (𝑥𝑘 , 𝑧𝑙 ) by a matrix (100*50). The density image, 𝜌(𝑥𝑘 , 𝑧𝑙 ), will have the
same dimensions as the data kernel, and it will be a matrix (100*50). One
can reshape these matrices into the row vectors 𝑙 (𝑗) and 𝜌 of the length
5000. You can introduce now a matrix F, with the rows, formed by the row
vectors 𝑙 (𝑗) . The Gram matrix consists of the combination of the dot
products (a sum over products of scalar components of elements) of the
data kernels, which can be written in matrix notations as
(12)
𝜞 = 𝑭 ∙ 𝑭𝑇
The linear system of equations (10) takes the matrix form:
(13)
𝜞∙𝜷=𝒅
where 𝛽 is the vector of unknown coefficients, and d is the vector of the
observed data. Solve system (13) using MATLAB.
Now, you can find the density distribution as a linear combination
of the data kernels with the coefficients 𝛽𝑖 , which can be written in matrix
form as
(14)
𝜌 = 𝐹𝑇𝛽
where 𝜌 is the vector of the density distribution. For visualization, this
vector can be reshaped back to matrix (100*50) of 𝜌(𝑥𝑘 , 𝑧𝑙 ).
Find the gravity field produced by this density
distribution (in other words, solve the forward modeling
problem). Input the density distribution and the
predicted field into Matlab and view it using your own
visualization routine or the one provided in the previous
assignment. The data files and visualization scripts will be
e-mailed to you.
Note that, the density image does not look like the
image from the previous assignment (hw3), but produces
the same data. These two solutions (from this assignment
and assignment #3) are called the “equivalent solutions”.
Turn in the plots of the obtained density image,
predicted field values, and the program listing.
𝑭=
(1)
𝑙1
(1)
𝑙2
(2)
𝑙1
(2)
𝑙2
⋮
⋮
(1)
𝑙41
⋯
⋯
(5000)
𝑙1
(5000)
𝑙2
⋮
(2)
𝑙41
⋯
𝑙𝑖
𝑥′
(5000)
𝑙41
𝛽1
𝑑01
0
𝛽2
𝑑
0
2
,𝜷 =
,𝒅 =
⋮
⋮
𝑑041
𝛽41
𝑧
= 2𝛾
𝑥 − 𝑥′ 2 + 𝑧2
𝜞 = 𝑭 ∙ 𝑭𝑇 (41 × 41)
𝜞𝜷 = 𝒅0
𝝆 = 𝑭𝑻 ∙ 𝜷
Gz(x),mGal
Vertical component of gravity field
0.1
0.05
0
0
100
200
300
400
500
600
700
800
900
1000
700
800
900
1000
Inversion result
z, meters
50
100
150
200
0
100
200
300
400
500
600
x, meters
Inversion Theory
Assignment 4
September30,2020
Contraction operator
An operator C of X into itself is called a contraction operator if there
exists a positive real number r < 1 with the property that (C (x) ; C (x0 ))
r (x; x0 ) for all points x and x0 in X.
Prove that the contraction operator is a continous operator.
Bounded operator
Linear operator y = L (x) is called bounded if there exists a real number
M with the property that
kL (x)k
M kxk
(1)
for every x 2 X.
Prove that the linear operator L is continuous if it is bounded.
Gram - Schmidt orthogonalization process
When we have a orthognormal set of elements, it will be much simpli ed
for solving the system of equation.The set of elements fd1 ; d2 ; d3 ; ; dn g
can be convered into a corresponding orthonormal set fe1 ; e2 ; e3 ; ; en g .
The rst step is
e1 =
d1
kd1 k
1
(2)
The e2 and e3 is
e2 =
e3 =
d3
kd3
(d2 ; e1 ) e1
(d2 ; e1 ) e1 k
d2
kd2
(d3 ; e1 ) e1
(d3 ; e1 ) e1
(d3 ; e2 ) e2
(d3 ; e2 ) e2 k
(3)
(4)
Write the e4 and e5 .
2-D Gravity Inverse problem
Let us consider the 2-D gravity problem. The anomalous gravity eld gz
is generated by some mass distributed in the rectangular domain S within
the lower half plane with the unknown excess density
(x; z). The observed
on the earth's surface gravity eld is given (same as in HW3).
Problem:
Solve the inverse gravity problem of determining the excess density distribution
(x; z) ; using Riesz representation theorem.
Directions:
The anomalous 2-D gravity eld can be calculated by the formula:
gz (x0; 0) = 2
Z Z
(x; z) z
dzdx:
(x x0)2 + z 2
S
(5)
We have observations at 41 points, which can be written as
di = gz (xi ; 0) = 2
=
Z Z
Z Z
(x; z)f2
S
S
(x
(x; z) z
dzdx =
(x xi )2 + z 2
z
gdzdx:
xi ) 2 + z 2
(6)
According to Riesz representation theorem, there exist the vectors l(j) (the
elements of the model space M ), which can be used to represent the observed
data:
dj = ( (x; z); l(j) );
(7)
2
j = 1; 2; :::n; ; l(j) 2 M:
The vectors l(j) are called \the data kernels". Comparing (7) and (6), we see
that
z
l(j) = l(j) (x; z) = 2
:
(8)
(x xj )2 + z 2
The unknown distribution of the density, according to the Riesz representation theorem, can be expressed as
(x; z) =
n
X
n
X
(i)
=2
il
i=1
where the coe cients
i
i
i=1
(x
z
;
xi ) 2 + z 2
(9)
satisfy the following system of linear equations:
dj = gz (xi ; 0) =
n
X
ji i :
(10)
i=1
We can calculate the Gram matrix as follows:
ji
= (l(i) ; l(j) ) = 4
2
RR
z2
S [(x xi )2 +z 2 ][(x xj )2 +z 2 ] dzdx;
(11)
j = 1; 2; ::::n; i = 1; 2; ::::n:
You are recommended to use Matlab for this assignment. Assume that
resulting density image consists of a number of small enough pixels. Take
100 cells in x direction and 50 cells in z direction. The image has to occupy
a rectangular subsurface region (0,1000) along the X axis and (25,225) along
the Z axis.
Note that in this approach, we approximate one data kernel l(j) (xk ; zl ) by
a matrix (100*50). The density image, (xk ; zl ); will have the same dimensions as the data kernel, and it will be a matrix (100*50). One can reshape
these matrices into the row vectors l(j) and of the length 5000. You can
introduce now a matrix F, with the rows, formed by the row vectors l(j) : The
Gram matrix consists of the combination of the dot products (a sum over
products of elements) of the data kernels, which can be written in matrix
notations as
= F FT :
(12)
The linear system of equations (10) takes the matrix form:
= d:
3
(13)
where
is the vector of unknown coe cients, and d is the vector of the
observed data. Solve the system (13) using MATLAB.
Now, you can nd the density distribution as a linear combination of the
data kernels with the coe cients i , which can be written in matrix form as
= FT
;
(14)
where is the vector of the density distribution. For visualization, this vector
can be reshaped back to the matrix (100*50) of (xk ; zl ):
Find the gravity eld, produced by this density distribution (in other
words, solve the forward modeling problem). Input the density distribution
and the predicted eld into Matlab and view it using your own visualization routine or the one provided in the previous assignment. Data les and
visualization scripts will E-mail to you.
Note that the density image does not look like the image from the previous
assignment (hw3), but produces the same data. These two solutions (from
this assignment and assignment #3) are called the \equivalent solutions".
Turn in the plots of the obtained density image, predicted eld values
and the program listing.
4
Purchase answer to see full
attachment