430 Inverse Problems
poorer ones for components on which the seismogram depends
weakly.
This issue can be addressed in several ways. One is to exclude
small eigenvalues from the inversion. Another, which avoids
the eigenvector decomposition, is to modify the function used
to measure the misfit between the data predicted by the model
and those observed (Eqn 7.2.11) to
χ
σ
ε
2
2
2
2
2
1
=
−
⎛
⎝
⎜ ⎜
⎞
⎠
⎟ ⎟
+
∑
∑
∑ ( ) .
i
i
i j
j
j
i
j
j
d
G m
m
∆
∆
∆
(34)
This function is the sum of the net misfit and the change in
length of the model vector, weighted by ε 2 . Hence minimizing
it is a compromise between the best fit to the data and the least
change from the starting model. The resulting solution, written
with the ∆s suppressed,
m = (G T G + ε 2 I) −1 G T d,
(35)
is called the damped least squares solution. If ε is zero, we have
the best-fit solution (Eqn 7.2.17), whereas larger values of ε
reduce or damp the change in the starting model by accepting
a poorer fit to the data. The damping parameter ε is chosen
empirically to yield a solution that seems plausible, and thus of
necessity reflects our ideas about the solution sought, because
damping the poorly constrained and undesired changes in the
model also damps the better constrained and desired changes.
Another common situation is that we want some data to have
greater effect on the solution, usually because we consider them
to be better known. We thus incorporate a data-weighting
matrix W d into the solution. The simplest is to weight by W d =
(σ 2
d ) −1 , the inverse of the variance–covariance matrix of the
data, so the data with the smallest uncertainties have the greatest effect. Problem 5 shows that this weighted least squares
solution is
m = (G T W d G) −1 G T W d d.
(36)
We may also want to have the model change smoothly,
such that each element varies only slightly with respect to its
neighbors. For instance, if the model were a continuous function of one variable, we measure the smoothness, or flatness,
f, of the changes by forming
f
m ,
=
−
⋅
⋅
−
⋅
⋅
⋅
⋅
⋅
⋅ ⋅ ⋅
⋅
⋅
⋅ ⋅ ⋅
⋅
−
⎛
⎝
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎞
⎠
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⋅
⋅
⎛
⎝
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎞
⎠
⎟
⎟
⎟
⎟
⎟
⎟
⎟
=
1 1 0
0
0 1 1
0
0 0 0
0
0
0
0 0 0 0 1 1
1
2
3
m
m
m
m
F
r
(37)
where F is the flatness matrix, which is a numerical approximation to the derivative at the edges of each element. The overall
flatness of the solution is then
f T f = m T F T Fm = m T W m m,
(38)
so the matrix W m = F T F is a weighting matrix for the model. For
more complicated model geometries, F is changed appropriately.
We can combine the model and data weighting in a weighted
damped least squares inversion, which yields the solution
m = (G
T
W d G + ε
2 W m )
−1 G
T
W d d.
(39)
As noted earlier, the damping parameter ε is chosen empirically. If we do not weight the data and model, the weighting
matrices W m and W d are identity matrices, and Eqn 39 is just
the simple damped least squares solution (Eqn 35).
An example of such an inversion was shown for P-wave
velocities at the base of the mantle in Fig. 3.5-17. A grid of 660
nodes that were roughly equally spaced were used to represent
the base of the mantle. The damping factor, ε = 1.2, was a compromise between the best fit, which minimizes the prediction
error, and minimizing the undetermined part of the solution.
Because each node is surrounded by 5 or 6 nodes that are
roughly equidistant, the rows of the model flatness matrix F
were chosen with the diagonal term equal to −1 and the terms
of the nearest N neighbors equal to 1/N (with N = 5 or 6). The
data were weighted empirically so that the diagonal elements
of the W d matrix ranked the quality of the observations from 9
(excellent) through 4 (good) to 1 (poor). These choices again
bear out that we have various ways of solving inverse problems, so the solution we develop depends on choices about
the data we use and the model we seek, based on our ideas
about what seems reasonable. Hence our solutions are in part
objective and in part subjective, and different approaches yield
different solutions.
7.3.5 Examples
Studies using travel time tomography yield interesting results
for various areas. For example, Fig. 7.3-5 (top) shows the model
geometry used in a study of the upper mantle in the region
including Central Europe, the Mediterranean, and the Middle
East. The model contains nine layers, each divided into 1040 1°
by 1° blocks. The layer thickness increases with depth from
33 km at the top to 130 km at a depth of 670 km. The data consist of approximately half a million travel times from about
25,000 earthquakes, recorded at stations both within the
model region and at distances to 90°.
The data used are travel time anomalies relative to the
Jeffreys–Bullen values, which can result from earthquake
mislocations as well as variations in seismic velocities. The
location and origin time of the earthquakes were thus also
inverted for, so the number of unknowns reflects both the
number of blocks (9360) and four times the number of earthquakes used. To reduce these numbers, procedures were used
to combine data from nearby earthquakes and from stations
close to each other. The problem to be solved thus involves
approximately 300,000 equations for 20,000 unknowns.
poorer ones for components on which the seismogram depends
weakly.
This issue can be addressed in several ways. One is to exclude
small eigenvalues from the inversion. Another, which avoids
the eigenvector decomposition, is to modify the function used
to measure the misfit between the data predicted by the model
and those observed (Eqn 7.2.11) to
χ
σ
ε
2
2
2
2
2
1
=
−
⎛
⎝
⎜ ⎜
⎞
⎠
⎟ ⎟
+
∑
∑
∑ ( ) .
i
i
i j
j
j
i
j
j
d
G m
m
∆
∆
∆
(34)
This function is the sum of the net misfit and the change in
length of the model vector, weighted by ε 2 . Hence minimizing
it is a compromise between the best fit to the data and the least
change from the starting model. The resulting solution, written
with the ∆s suppressed,
m = (G T G + ε 2 I) −1 G T d,
(35)
is called the damped least squares solution. If ε is zero, we have
the best-fit solution (Eqn 7.2.17), whereas larger values of ε
reduce or damp the change in the starting model by accepting
a poorer fit to the data. The damping parameter ε is chosen
empirically to yield a solution that seems plausible, and thus of
necessity reflects our ideas about the solution sought, because
damping the poorly constrained and undesired changes in the
model also damps the better constrained and desired changes.
Another common situation is that we want some data to have
greater effect on the solution, usually because we consider them
to be better known. We thus incorporate a data-weighting
matrix W d into the solution. The simplest is to weight by W d =
(σ 2
d ) −1 , the inverse of the variance–covariance matrix of the
data, so the data with the smallest uncertainties have the greatest effect. Problem 5 shows that this weighted least squares
solution is
m = (G T W d G) −1 G T W d d.
(36)
We may also want to have the model change smoothly,
such that each element varies only slightly with respect to its
neighbors. For instance, if the model were a continuous function of one variable, we measure the smoothness, or flatness,
f, of the changes by forming
f
m ,
=
−
⋅
⋅
−
⋅
⋅
⋅
⋅
⋅
⋅ ⋅ ⋅
⋅
⋅
⋅ ⋅ ⋅
⋅
−
⎛
⎝
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎞
⎠
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⋅
⋅
⎛
⎝
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎞
⎠
⎟
⎟
⎟
⎟
⎟
⎟
⎟
=
1 1 0
0
0 1 1
0
0 0 0
0
0
0
0 0 0 0 1 1
1
2
3
m
m
m
m
F
r
(37)
where F is the flatness matrix, which is a numerical approximation to the derivative at the edges of each element. The overall
flatness of the solution is then
f T f = m T F T Fm = m T W m m,
(38)
so the matrix W m = F T F is a weighting matrix for the model. For
more complicated model geometries, F is changed appropriately.
We can combine the model and data weighting in a weighted
damped least squares inversion, which yields the solution
m = (G
T
W d G + ε
2 W m )
−1 G
T
W d d.
(39)
As noted earlier, the damping parameter ε is chosen empirically. If we do not weight the data and model, the weighting
matrices W m and W d are identity matrices, and Eqn 39 is just
the simple damped least squares solution (Eqn 35).
An example of such an inversion was shown for P-wave
velocities at the base of the mantle in Fig. 3.5-17. A grid of 660
nodes that were roughly equally spaced were used to represent
the base of the mantle. The damping factor, ε = 1.2, was a compromise between the best fit, which minimizes the prediction
error, and minimizing the undetermined part of the solution.
Because each node is surrounded by 5 or 6 nodes that are
roughly equidistant, the rows of the model flatness matrix F
were chosen with the diagonal term equal to −1 and the terms
of the nearest N neighbors equal to 1/N (with N = 5 or 6). The
data were weighted empirically so that the diagonal elements
of the W d matrix ranked the quality of the observations from 9
(excellent) through 4 (good) to 1 (poor). These choices again
bear out that we have various ways of solving inverse problems, so the solution we develop depends on choices about
the data we use and the model we seek, based on our ideas
about what seems reasonable. Hence our solutions are in part
objective and in part subjective, and different approaches yield
different solutions.
7.3.5 Examples
Studies using travel time tomography yield interesting results
for various areas. For example, Fig. 7.3-5 (top) shows the model
geometry used in a study of the upper mantle in the region
including Central Europe, the Mediterranean, and the Middle
East. The model contains nine layers, each divided into 1040 1°
by 1° blocks. The layer thickness increases with depth from
33 km at the top to 130 km at a depth of 670 km. The data consist of approximately half a million travel times from about
25,000 earthquakes, recorded at stations both within the
model region and at distances to 90°.
The data used are travel time anomalies relative to the
Jeffreys–Bullen values, which can result from earthquake
mislocations as well as variations in seismic velocities. The
location and origin time of the earthquakes were thus also
inverted for, so the number of unknowns reflects both the
number of blocks (9360) and four times the number of earthquakes used. To reduce these numbers, procedures were used
to combine data from nearby earthquakes and from stations
close to each other. The problem to be solved thus involves
approximately 300,000 equations for 20,000 unknowns.