Lecture 24 Exemplary Inverse Problems - PowerPoint Presentation

Slide1

Lecture 24

Exemplary Inverse Problems

including

Vibrational

ProblemsSlide2

Syllabus

Lecture 01 Describing Inverse Problems

Lecture 02 Probability and Measurement Error, Part 1

Lecture 03 Probability and Measurement Error, Part 2

Lecture 04 The L

2

Norm and Simple Least Squares

Lecture 05 A Priori Information and Weighted Least Squared

Lecture 06 Resolution and Generalized Inverses

Lecture 07 Backus-Gilbert Inverse and the Trade Off of Resolution and Variance

Lecture 08 The Principle of Maximum Likelihood

Lecture 09 Inexact Theories

Lecture 10

Nonuniqueness

and Localized Averages

Lecture 11 Vector Spaces and Singular Value Decomposition

Lecture 12 Equality and Inequality Constraints

Lecture 13 L

1

, L

∞

Norm Problems and Linear Programming

Lecture 14 Nonlinear Problems: Grid and Monte Carlo Searches

Lecture 15 Nonlinear Problems: Newton’s Method

Lecture 16 Nonlinear Problems: Simulated Annealing and Bootstrap Confidence Intervals

Lecture 17 Factor Analysis

Lecture 18

Varimax

Factors,

Empircal

Orthogonal Functions

Lecture 19 Backus-Gilbert Theory for Continuous Problems; Radon’s Problem

Lecture 20 Linear Operators and Their

Adjoints

Lecture 21

Fr

é

chet

Derivatives

Lecture 22 Exemplary Inverse Problems, incl. Filter Design

Lecture 23 Exemplary Inverse Problems, incl. Earthquake Location

Lecture 24 Exemplary Inverse Problems, incl.

Vibrational

ProblemsSlide3

Purpose of the Lecture

solve a few exemplary inverse problems

tomography

vibrational

problems

determining mean directionsSlide4

Part 1

tomographySlide5

d

i

= ∫

ray

i

m(x(s), y(s))

d

s

ray

i

tomography:

data is line integral of model function

assume ray path is known

x

ySlide6

discretization

:

model function divided up into

M

pixels

m

jSlide7

data kernel

G

ij

= length of ray

i

in pixel

jSlide8

data kernel

G

ij

= length of ray

i

in pixel

j

here’s an

easy,

approximate way to calculate itSlide9

ray

i

start with

G

set to zero

then consider each ray in sequenceSlide10

∆s

divide

each ray into segments of arc length

∆s

and step from segment to segment Slide11

determine the pixel

index, say

j

, that the

center

of each line segment falls within

add

∆s

to

G

ij

repeat for every segment of every ray

Slide12

You can make this approximation indefinitely accurate simply by

decreasing the size of

∆

s

(albeit at the expense of increase the computation time)Slide13

Suppose that there are

M=L

2

voxels

A ray passes through about L

voxels

G

has

NL

2 elementsNL

of which are non-zeroso the fraction of non-zero elements is1/L

hence

G

is very sparseSlide14

In a typical

tomographic

experiment

some pixels will be missed entirely

and some groups of pixels will be sampled by only one ray Slide15

In a typical

tomographic

experiment

some pixels will be missed entirely

and 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

M>N

as well

)Slide16

so

you must introduce some sort of a priori information to achieve a solution

say

a priori information that the solution is small

or

a priori information that the solution is smoothSlide17

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 functionSlide18

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 unnecessarySlide19

True model

x

y

sources and receiversSlide20

x

y

Ray Coverage

x

y

just a “few” rays shown

else image is blackSlide21

Data, plotted in Radon-style coordinates

angle

θ

of ray

distance

r

of ray to center of image

Lesson from Radon’s Problem:

Full data coverage need to achieve exact solution

minor data gapsSlide22

Estimated model

x

y

True modelSlide23

Estimated model

x

y

Estimated model

streaks due to minor data gaps

they disappear if ray density is doubledSlide24

but what if the observational geometry is poor

so that broads swaths of rays are missing ?Slide25

(A)

(B)

(C)

(D)

x

y

x

y

x

y

θ

r

complete angular coverageSlide26

(A)

(B)

(C)

(D)

x

y

x

y

x

y

θ

r

incomplete angular coverageSlide27

Part 2

vibrational

problemsSlide28

statement of the problem

Can you determine the structure of an object

just knowing the

characteristic frequencies at which it vibrates?

frequencySlide29

the

Fréchet

derivative

of frequency with respect to velocity

is usually computed using

perturbation theory

hence a quick discussion of what that is ...Slide30

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 knownSlide31

simple exampleSlide32
Slide33

we know the solution to this equation:

x

0

=±cSlide34
Slide35
Slide36
Slide37

Here’s the actual

vibrational

problem

acoustic equation with

spatially variable sound velocity

vSlide38

acoustic equation with

spatially variable sound velocity

v

frequencies of vibration

or

eigenfrequencies

patterns of vibration

or

eigenfunctions

or

modesSlide39

v(

x

) = v

(0)

(

x

) +

ε

v

(1)(x) + ...assume velocity

can be written as a perturbationaround some simple structurev(0)(x)

Slide40

eigenfunctions

known to obey

orthonormality

relationshipSlide41

now represent

eigenfrequencies

and

eigenfunctions

as power series in

εSlide42

represent

first-order perturbed

shapes as sum of

unperturbed shapes

now represent

eigenfrequencies

and

eigenfunctions

as power series in

εSlide43

plug series into original differential equation

group terms of equal power of

ε

solve for first-order

perturbation

in

eigenfrequencies

ω

n

(1)and eigenfunction coefficients b

nm

(use

orthonormality in process) Slide44

resultSlide45

result for

eigenfrequencies

write as standard inverse problemSlide46

standard continuous

inverse problemSlide47

standard continuous

inverse problem

perturbation in the

eigenfrequencies

are the data

perturbation in the

velocity structure is

the model functionSlide48

standard continuous

inverse problem

depends upon the

unperturbed velocity structure,

the unperturbed

eigenfrequency

and the unperturbed mode

data kernel

or

Fréchet

derivativeSlide49

1D organ pipe

unperturbed

problem has constant velocity

0

h

x

open end,

p=0

closed end

dp

/

dx

=0

perturbed

problem has variable velocitySlide50

0

h

x

p=0

dp

/

dx

=0

p

1

x

x

p

2

p

3

x

𝜔

1

modes

frequencies

𝜔

2

𝜔

3

𝜔

0Slide51

solution to unperturbed problemSlide52

position , x

velocity, v

perturbed

unperturbed

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0

2

4

6

8

10

velocity structureSlide53

How to

discretize

the model function?

m

is

veloctity

function evaluated at sequence of points equally spaced in

x

our choice is very simpleSlide54

the data

a list of frequencies of vibration

true, unperturbed

true, perturbed

observed = true, perturbed + noise

frequencySlide55

ω

i

m

j

the data kernelSlide56

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 functionSlide57

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 choiceSlide58

position , x

velocity, v

the solution

true

estimatedSlide59

position , x

velocity, v

the solution

true

estimatedSlide60

m

i

m

j

the model resolution matrixSlide61

m

i

m

j

the model resolution matrix

what is this?Slide62

This problem has a type of

nonuniqueness

that arises from its symmetry

a positive velocity anomaly at one end of the organ pipe

trades off with a negative anomaly at the other endSlide63

this behavior is very common

and is why

eigenfrequency

data

are usually supplemented with other data

e.g. travel times along rays

that are not subject to this

nonuniquenessSlide64

Part 3

determining mean directionsSlide65

statement of the problem

you measure a bunch of directions (unit vectors)

what’s their mean?

x

y

z

data

meanSlide66

what’s a reasonable

probability density function

for directional data?

Gaussian doesn’t quite work

because

its defined on the wrong interval

(-

∞, +∞)Slide67

θ

ϕ

central vector

datum

coordinate system

distribution should

be

symmetric

in

ϕSlide68

angle,

θ

π

-

π

p(

θ

)

Fisher distribution

similar in shape to a Gaussian but on a sphereSlide69

angle,

θ

π

-

π

p(

θ

)

Fisher distribution

similar in shape to a Gaussian but on a sphere

“precision parameter”

quantifies width of

p.d.f

.

κ

=5

κ

=1Slide70

solve by

direct application of

principle of maximum likelihoodSlide71

maximize joint

p.d.f

. of data

with respect to

κ and

cos

(

θ

)

Slide72

x

: Cartesian components of observed unit vectors

m

: Cartesian components of central unit vector; must constrain

|

m

|=1

Slide73

likelihood function

constraint

unknowns

m

,

κ

C =

= 0 Slide74

Lagrange multiplier equationsSlide75

Results

valid when

κ

>5Slide76

Results

central vector is parallel to the vector that you get by putting all the observed unit vectors end-to-endSlide77

Solution Possibilities

Determine

m

by evaluating simple formula

Determine

κ

using simple but approximate formula

2. Determine

κ

using bootstrap method

our choice

only valid when

κ

>5Slide78

Application to

Subduction

Zone Stresses

Determine the mean direction of

P-axes

of deep (300-600 km) earthquakes

in the Kurile-Kamchatka

subduction

zoneSlide79

N

E

data

central direction

bootstrap