Title: Lecture 24 Exemplary Inverse Problems including Vibrational Problems
1Lecture 24 Exemplary Inverse Problemsincluding
Vibrational Problems
2Syllabus
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
3Purpose of the Lecture
solve a few exemplary inverse problems
tomography vibrational problems determining
mean directions
4Part 1
tomography
5tomography 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
6discretization model function divided up into M
pixels mj
7data kernelGij length of ray i in pixel j
8data kernelGij length of ray i in pixel j
heres an easy, approximate way to calculate it
9start with G set to zero
ray i
then consider each ray in sequence
10divide each ray into segments of arc length ?s
?s
and step from segment to segment
11determine 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
12You can make this approximation indefinitely
accurate simply bydecreasing the size of
?s(albeit at the expense of increase the
computation time)
13Suppose 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
14In a typical tomographic experimentsome pixels
will be missed entirelyand some groups of
pixels will be sampled by only one ray
15In 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)
16soyou 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
17Solution 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
18Solution 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
19sources and receivers
True model
y
x
20just a few rays shown else image is black
Ray Coverage
y
x
21minor data gaps
Lesson from Radons Problem Full data coverage
need to achieve exact solution
22True model
23Estimated model
streaks due to minor data gaps they disappear if
ray density is doubled
24but what if the observational geometry is
poorso that broads swaths of rays are missing ?
25complete angular coverage
26incomplete angular coverage
27Part 2
vibrational problems
28statement of the problem
- Can you determine the structure of an object
- just knowing the
- characteristic frequencies at which it vibrates?
frequency
29the Fréchet derivativeof frequency with respect
to velocity is usually computed using
perturbation theoryhence a quick discussion of
what that is ...
30perturbation 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
31simple example
32(No Transcript)
33we know the solution to this equation x0c
34(No Transcript)
35(No Transcript)
36(No Transcript)
37Heres the actual vibrational problemacoustic
equation withspatially variable sound velocity v
38acoustic equation withspatially variable sound
velocity v
patterns of vibration or eigenfunctions or modes
frequencies of vibration or eigenfrequencies
39v(x) v(0)(x) ev(1)(x) ...
assume velocity can be written as a
perturbation around some simple structure v(0)(x)
40eigenfunctions known to obey orthonormality
relationship
41now represent eigenfrequencies and eigenfunctions
as power series in e
42now represent eigenfrequencies and eigenfunctions
as power series in e
represent first-order perturbed shapes as sum of
unperturbed shapes
43plug 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)
44result
45result for eigenfrequencies
write as standard inverse problem
46standard continuous inverse problem
47standard continuous inverse problem
perturbation in the eigenfrequencies are the data
perturbation in the velocity structure is the
model function
48standard continuous inverse problem
data kernel or Fréchet derivative
depends upon the unperturbed velocity
structure, the unperturbed eigenfrequency and the
unperturbed mode
491D organ pipe
open end, p0
unperturbed problem has constant velocity
0
perturbed problem has variable velocity
closed end dp/dx0
h
x
50p1
p2
p3
p0
0
modes
h
dp/dx0
x
x
x
x
frequencies
??
??1
??2
??3
0
51solution to unperturbed problem
52velocity structure
velocity, v
position , x
53How to discretize the model function?
our choice is very simple
- m is veloctity function evaluated at sequence of
points equally spaced in x
54the dataa list of frequencies of vibration
frequency
true, unperturbed
true, perturbed
observed true, perturbed noise
55the data kernel
56Solution 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
57Solution 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
58the solution
estimated
true
59the solution
estimated
true
60the model resolution matrix
61the model resolution matrix
what is this?
62This 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
63this behavior is very commonand is why
eigenfrequency dataare usually supplemented with
other datae.g. travel times along raysthat
are not subject to this nonuniqueness
64Part 3
determining mean directions
65statement of the problem
- you measure a bunch of directions (unit vectors)
- whats their mean?
data
mean
66whats a reasonableprobability density
functionfor directional data?
Gaussian doesnt quite work because its defined
on the wrong interval (-8, 8)
67coordinate system
central vector
?
datum
?
distribution should be symmetric in ?
68Fisher distributionsimilar in shape to a
Gaussian but on a sphere
69Fisher distributionsimilar in shape to a
Gaussian but on a sphere
?5
?1
precision parameter quantifies width of p.d.f.
70solve bydirect application ofprinciple of
maximum likelihood
71maximize joint p.d.f. of data
with respect to ? and cos(?)
72x Cartesian components of observed unit
vectorsm Cartesian components of central unit
vector must constrain m1
73likelihood function
constraint
0
C
unknowns m, ?
74Lagrange multiplier equations
75Results
valid when ?gt5
76Results
central vector is parallel to the vector that you
get by putting all the observed unit vectors
end-to-end
77Solution Possibilities
- Determine m by evaluating simple formula
- Determine ? using simple but approximate formula
- 2. Determine ? using bootstrap method
only valid when ?gt5
our choice
78Application to Subduction Zone Stresses
- Determine the mean direction of
- P-axes
- of deep (300-600 km) earthquakes
- in the Kurile-Kamchatka subduction zone
79N
E