Lecture 24 Exemplary Inverse Problems including Vibrational Problems

1
Lecture 24 Exemplary Inverse Problemsincluding
Vibrational Problems
2
Syllabus
Lecture 01 Describing Inverse ProblemsLecture
02 Probability and Measurement Error, Part
1Lecture 03 Probability and Measurement Error,
Part 2 Lecture 04 The L2 Norm and Simple Least
SquaresLecture 05 A Priori Information and
Weighted Least SquaredLecture 06 Resolution and
Generalized Inverses Lecture 07 Backus-Gilbert
Inverse and the Trade Off of Resolution and
VarianceLecture 08 The Principle of Maximum
LikelihoodLecture 09 Inexact TheoriesLecture
10 Nonuniqueness and Localized AveragesLecture
11 Vector Spaces and Singular Value
Decomposition Lecture 12 Equality and Inequality
ConstraintsLecture 13 L1 , L8 Norm Problems and
Linear ProgrammingLecture 14 Nonlinear
Problems Grid and Monte Carlo Searches Lecture
15 Nonlinear Problems Newtons Method Lecture
16 Nonlinear Problems Simulated Annealing and
Bootstrap Confidence Intervals Lecture
17 Factor AnalysisLecture 18 Varimax Factors,
Empircal Orthogonal FunctionsLecture
19 Backus-Gilbert Theory for Continuous
Problems Radons ProblemLecture 20 Linear
Operators and Their AdjointsLecture 21 Fréchet
DerivativesLecture 22 Exemplary Inverse
Problems, incl. Filter DesignLecture 23
Exemplary Inverse Problems, incl. Earthquake
LocationLecture 24 Exemplary Inverse Problems,
incl. Vibrational Problems
3
Purpose of the Lecture
solve a few exemplary inverse problems
tomography vibrational problems determining
mean directions
4
Part 1
tomography
5
tomography data is line integral of model
function
y
assume ray path is known
ray i
x
di ?ray i m(x(s), y(s)) ds
6
discretization model function divided up into M
pixels mj
7
data kernelGij length of ray i in pixel j
8
data kernelGij length of ray i in pixel j
heres an easy, approximate way to calculate it
9
start with G set to zero
ray i
then consider each ray in sequence
10
divide each ray into segments of arc length ?s
?s
and step from segment to segment
11
determine the pixel index, say j, that the center
of each line segment falls within
add ?s to Gij repeat for every segment of every
ray
12
You can make this approximation indefinitely
accurate simply bydecreasing the size of
?s(albeit at the expense of increase the
computation time)
13
Suppose that there are ML2 voxelsA ray passes
through about L voxelsG has NL2 elementsNL of
which are non-zeroso the fraction of non-zero
elements is1/Lhence G is very sparse
14
In a typical tomographic experimentsome pixels
will be missed entirelyand some groups of
pixels will be sampled by only one ray
15
In a typical tomographic experimentsome pixels
will be missed entirelyand some groups of
pixels will be sampled by only one ray
the value of these pixels is completely
undetermined
only the average value of these pixels is
determined
hence the problem is mixed-determined (and
usually MgtN as well)
16
soyou must introduce some sort of a priori
information to achieve a solutionsaya priori
information that the solution is smallor a
priori information that the solution is smooth
17
Solution Possibilities

• Damped Least Squares (implements smallness)
• Matrix G is sparse and very large
• use bicg() with damped least squares function
• 2. Weighted Least Squares (implements
  smoothness)
• Matrix F consists of G plus
• second derivative smoothing
• use bicg()with weighted least squares function

18
Solution Possibilities

• Damped Least Squares
• Matrix G is sparse and very large
• use bicg() with damped least squares function
• 2. Weighted Least Squares
• Matrix F consists of G plus
• second derivative smoothing
• use bicg()with weighted least squares function

test case has very good ray coverage, so
smoothing probably unnecessary
19
sources and receivers
True model
y
x
20
just a few rays shown else image is black
Ray Coverage
y
x
21
minor data gaps
Lesson from Radons Problem Full data coverage
need to achieve exact solution
22
True model
23
Estimated model
streaks due to minor data gaps they disappear if
ray density is doubled
24
but what if the observational geometry is
poorso that broads swaths of rays are missing ?
25
complete angular coverage
26
incomplete angular coverage
27
Part 2
vibrational problems
28
statement of the problem

• Can you determine the structure of an object
• just knowing the
• characteristic frequencies at which it vibrates?

frequency
29
the Fréchet derivativeof frequency with respect
to velocity is usually computed using
perturbation theoryhence a quick discussion of
what that is ...
30
perturbation theory

• a technique for computing an approximate solution
  to a complicated problem, when
• 1. The complicated problem is related to a simple
  problem by a small perturbation
• 2. The solution of the simple problem must be
  known

31
simple example
32
(No Transcript)
33
we know the solution to this equation x0c
34
(No Transcript)
35
(No Transcript)
36
(No Transcript)
37
Heres the actual vibrational problemacoustic
equation withspatially variable sound velocity v
38
acoustic equation withspatially variable sound
velocity v
patterns of vibration or eigenfunctions or modes
frequencies of vibration or eigenfrequencies
39
v(x) v(0)(x) ev(1)(x) ...
assume velocity can be written as a
perturbation around some simple structure v(0)(x)
40
eigenfunctions known to obey orthonormality
relationship
41
now represent eigenfrequencies and eigenfunctions
as power series in e
42
now represent eigenfrequencies and eigenfunctions
as power series in e
represent first-order perturbed shapes as sum of
unperturbed shapes
43
plug series into original differential equation
group terms of equal power of e
solve for first-order perturbation in
eigenfrequencies ?n(1) and eigenfunction
coefficients bnm (use orthonormality in
process)
44
result
45
result for eigenfrequencies
write as standard inverse problem
46
standard continuous inverse problem
47
standard continuous inverse problem
perturbation in the eigenfrequencies are the data
perturbation in the velocity structure is the
model function
48
standard continuous inverse problem
data kernel or Fréchet derivative
depends upon the unperturbed velocity
structure, the unperturbed eigenfrequency and the
unperturbed mode
49
1D organ pipe
open end, p0
unperturbed problem has constant velocity
0
perturbed problem has variable velocity
closed end dp/dx0
h
x
50
p1
p2
p3
p0
0
modes
h
dp/dx0
x
x
x
x
frequencies
??
??1
??2
??3
0
51
solution to unperturbed problem
52
velocity structure
velocity, v
position , x
53
How to discretize the model function?
our choice is very simple

• m is veloctity function evaluated at sequence of
  points equally spaced in x

54
the dataa list of frequencies of vibration
frequency
true, unperturbed
true, perturbed
observed true, perturbed noise
55
the data kernel
56
Solution Possibilities

• Damped Least Squares (implements smallness)
• Matrix G is not sparse
• use bicg() with damped least squares function
• 2. Weighted Least Squares (implements
  smoothness)
• Matrix F consists of G plus
• second derivative smoothing
• use bicg()with weighted least squares function

57
Solution Possibilities

• Damped Least Squares (implements smallness)
• Matrix G is not sparse
• use bicg() with damped least squares function
• 2. Weighted Least Squares (implements
  smoothness)
• Matrix F consists of G plus
• second derivative smoothing
• use bicg()with weighted least squares function

our choice
58
the solution
estimated
true
59
the solution
estimated
true
60
the model resolution matrix
61
the model resolution matrix
what is this?
62
This problem has a type of nonuniquenessthat
arises from its symmetrya positive velocity
anomaly at one end of the organ pipetrades off
with a negative anomaly at the other end
63
this behavior is very commonand is why
eigenfrequency dataare usually supplemented with
other datae.g. travel times along raysthat
are not subject to this nonuniqueness
64
Part 3
determining mean directions
65
statement of the problem

• you measure a bunch of directions (unit vectors)
• whats their mean?

data
mean
66
whats a reasonableprobability density
functionfor directional data?
Gaussian doesnt quite work because its defined
on the wrong interval (-8, 8)
67
coordinate system
central vector
?
datum
?
distribution should be symmetric in ?
68
Fisher distributionsimilar in shape to a
Gaussian but on a sphere
69
Fisher distributionsimilar in shape to a
Gaussian but on a sphere
?5
?1
precision parameter quantifies width of p.d.f.
70
solve bydirect application ofprinciple of
maximum likelihood
71
maximize joint p.d.f. of data
with respect to ? and cos(?)
72
x Cartesian components of observed unit
vectorsm Cartesian components of central unit
vector must constrain m1
73
likelihood function
constraint
0
C
unknowns m, ?
74
Lagrange multiplier equations
75
Results
valid when ?gt5
76
Results
central vector is parallel to the vector that you
get by putting all the observed unit vectors
end-to-end
77
Solution Possibilities

• Determine m by evaluating simple formula
• Determine ? using simple but approximate formula
• 2. Determine ? using bootstrap method

only valid when ?gt5
our choice
78
Application to Subduction Zone Stresses

• Determine the mean direction of
• P-axes
• of deep (300-600 km) earthquakes
• in the Kurile-Kamchatka subduction zone

79
N
E