PART 3

Intro to inversion and sparse

regularization
Mauricio Sacchi
University of Alberta, Edmonton, AB, Canada

April 2014

Outline

• Geophysics & Sensing

• Forward Problems, Inverse Problems & Regularization Methods

• Quadratic and non-quadratic regularization techniques

• Examples
– Deconvolution
– AVO

• Discussions

Part 3 Intro to inversion 2

1
Geophysics & Sensing

What you would like What you can

to know measure

MODEL DATA

Part 3 Intro to inversion 3

Geophysics & Sensing

What you would like What you can

to know measure

Equations of
Math Physics
MODEL DATA

Part 3 Intro to inversion 4

2
Geophysics & Sensing

Unknown Physical Observations (x,y,z=0)

Parameters in the Earth’s What you can measure !
interior (x,y,z)
What you would like to know !

Density Gravity

Elastic Parameters Seismic Wavefiels

Susceptibility Magnetic Fields

Conductivity Electric Potential

Part 3 Intro to inversion 5

Geophysics & Sensing

What you would like What you can

to know measure

Equations of
Math Physics
MODEL DATA

Part 3 Intro to inversion 6

3
Geophysics & Sensing

What you would like What you can

to know measure

Equations of
Math Physics
MODEL DATA
Inverse Theory

Part 3 Intro to inversion 7

Inverse Theory
Forward Problem

MODEL DATA
Inverse Problem

Remarks:

§ Model:
Properties are continuously distributed on a 3-D volume
§ Data:
Observations are sparsely distributed on the surface of the Earth

Part 3 Intro to inversion 8

4
Inverse Theory
Forward Problem

MODEL DATA
Inverse Problem


F(m) = d


m = m(x, y,z) d = [d1,d2 ,d3 ....,dN ]T

F:€Operator that maps Model Parameters to Observations

Finding F −1 is non trivial; we are dealing with ill-posed problems

Part 3 Intro to inversion 9

A Well-Posed Problem

J. Hadamard (1923)

• For all data, a solution exists

• For all data, the solution is unique
• The solution depends continuously on the data

Otherwise: ill-posed

Part 3 Intro to inversion 10

5
 Ill-posed problems

Why a problem is ill-posed?

• Insufficient data
• Noise
• Inadequate sampling, aperture
• Inadequate band-width
• Poor illumination
• Math-Phys equations (averaging property of forward kernels)

Proper description of physical phenomena should lead to well-

posed problems….. (not of much help in practical situations)

Part 3 Intro to inversion 11

Averaging Property of Kernels

g(x)
Forward problem - smoothes the edge

Smooth gravity profile

€
Inversion - unsmoothing

ρ1 ρ1 €
Sharp change in density

€ g(x) =€∫ K(x', x) ρ(x')dx'

Part 3 Intro to inversion 12

€

6
• Forward problem produce data that often looks like a “smooth” version of material
properties
• Smoothing is a stable process
• The process of undoing smoothing is unstable (Inversion is about controlling
unstability)

Smoothing is a low pass operation

Un-smoothing is a high pass operator

Assignment: Consider travel-time modeling via Dix formula and

look into the inverse problem.

Part 3 Intro to inversion 13

Regularization Methods to Solve Inverse Problems

• Looking for models that explain observations (Tikhonov, 1936)

Examples:
• Tomography
– Medical
– Seismic
• Inverse Scattering in Physics
• Inverse Heat Conduction Problems
• Deconvolution
– Geophysics
– Astronomy
– image processing
• Design of semi-conductors
• Non-destructive & non-invasive testing
• Wireless communication
• Data Analysis
– Density Kernel Estimators
– Interpolation/Data Reconstruction
– Harmonic Analysis
– Machine Learning

Part 3 Intro to inversion 14

7
 Regularization Methods to Solve Inverse Problems

• Replacing an ill-posed problem by a parameter dependent family of well-posed

problems (.. Let’s find a solution)

• Additional information is required about the model

– Positivity
– Compact support
– Smoothness
– Sparseness
– Common sense [talk to geologists and engineers]

Proper description of the problem and a priori information about the

unknown model should lead to proper regularization strategies…

Part 3 Intro to inversion 15

An introduction to regularization with a

minimum amount of math

- A trivial problem (toy problem)

Part 3 Intro to inversion 16

8
 Lab experiment

v1 v2

€ L1 € L2
L1 L2
+ ≈t
v1 v 2
0 t
€ €
L1m1 + L2 m2 ≈ d

m 
(L1 L2 ) 1  ≈ d or Lm ≈ d
 m2 

Part 3 Intro to inversion 17

Model Space
m2
L1 m1 + L2 m2 ≈ d

€
€

m1

Part 3 Intro to inversion 18

9
 Model Space
m2
L1 m1 + L2 m2 ≈ d

€
€

m1
Solution is anywhere close to the line
Part 3 Intro to inversion 19

Smallest Solution
m2 m1 m2
Mimimize R = m12 + m22
subject to L1m1 + L2 m2 ≈ d
€

€
€

m1

Part 3 Intro to inversion 20

10
 Smallest Solution
m2 m1 m2
2 2
Mimimize R = m + m
1 2

subject to L1m1 + L2 m2 ≈ d
€

€
€

m1

Part 3 Intro to inversion 21

Smallest Solution
m2 m1 m2
Mimimize R = m12 + m22
subject to L1m1 + L2 m2 ≈ d
€ These solutions were avoided :
m1 
→ ∞
€ m2 
→ −∞
€ and
m1 
→ −∞
m2 
→ ∞

m1
Smallest Solution = Minimum Norm Solution

Part 3 Intro to inversion 22

11
 Smallest Perturbation
m2 m1 m2
2 2
Mimimize R = (m 1− mi1 ) + (m 2− mi2 )
subject to L1m1 + L2 m2 ≈ d
€

€ (mi1,mi2 )
€

m1

Part 3 Intro to inversion 23

Smallest Perturbation
m2 m1 m2
Mimimize R = (m 1− mi1 ) 2 + (m 2− mi2 ) 2
subject to L1m1 + L2 m2 ≈ d
€

€ (mi1,mi2 )
€

m1

Part 3 Intro to inversion 24

12
 Flattest Solution
m2 m1 m2
2
Mimimize R = (m2 − m1 )
subject to L1m1 + L2 m2 ≈ d
€

€
€

m1

Part 3 Intro to inversion 25

Flattest Solution
m2 m1 m2
2
Mimimize R = (m2 − m1 )
subject to L1m1 + L2 m2 ≈ d
€

€
€

m1

Part 3 Intro to inversion 26

13
 Sparse Solution*
m2 m1 m2
Mimimize R =| m1 | + | m2 |
subject to L1m1 + L2 m2 ≈ d
€

€
€

m1

Part 3 Intro to inversion 27

Sparse Solution*
m2 m1 m2
Mimimize R =| m1 | + | m2 |
subject to L1m1 + L2 m2 ≈ d
€

€
€

m1
* Parsimonious or minimum entropy solution

Part 3 Intro to inversion 28

14
 Problems in exploration seismology
where regularization methods are
important

• Data Reconstruction
• Focusing with Radon Transforms
• Deconvolution
• Probing a reflector with seismic waves (AVO)
• Imaging (regularized LS imaging/Migration)

These problems can all be tackled with regularization methods

Part 3 Intro to inversion 29

Deconvolution
- debluring
- equalization

Increase resolution and our ability

to “see” thin layers, stars,
imperfections in materials, etc etc

Part 3 Intro to inversion 30

15
 Seismic Deconvolution

Earth versus depth Earth versus time

s r

Time
Depth

3
4

r W d0 d
d(t) = ∫ w(τ - t)r(τ )dτ + n(t)
d = d0 + n
= Wr + n
€ € € €
Part 3 Intro to inversion 31

d(t) = ∫ w(τ - t)r(τ )dτ + n(t)

d = Wr + n

Ill-posed ?
This term can be extremely
D(ω ) = D0 (ω ) + N(ω ) large when the source
amplitude spectrum is
= W (ω )R(ω ) + N(ω ) small

D(ω ) D0 (ω ) + N(ω ) N(ω )

Rˆ (ω ) = = = R(ω ) +
W (ω ) W (ω ) W (ω )

€
Part 3 Intro to inversion 32

16
Seismic Deconvolution
The very Naïve Solution (Just try to fit the data…)

d =W r+n

J =||Wr − d ||22

∇J = W T W r − W ' d = 0

rˆ = (W T W ) −1W T d

Small perturbations on d will produce large perturbations on r

€
Part 3 Intro to inversion 33

Seismic Deconvolution
d =W r+n
Low Resolution Solution:

J =||Wr − d ||22

r = (W T W ) −1W T d

c I ≈ W TW

r˜ ≈ c −1W T d
Some sort of “backprojection”
Matching filter (Good to suppress noise)

Part 3 Intro to inversion 34

€

17
 Example
Time

r W d0 d r˜

Right phase but pulse has not been compressed

Part 3 Intro to inversion 35

€ € € € €

Seismic Deconvolution
Quadratic Regularization (add a quadratic penalty term)

d =W r+n Quadratic
penalty term
J =||Wr − d ||22 + µ || Lr ||22

Misfit ∇J = W T W r − W T d + µLT L r = 0

r = (W T W + µLT L) −1W T d

L=I⇒

Final solution: estimator

of the reflectivity r = (W 'W + µI) −1W ' d

Part 3 Intro to inversion € 36

18
 Example
Time

r W d0 d r˜ rq

Right phase, pulse has been compressed, ringing due to minimum norm requirement

Part 3 Intro to inversion 37

€ € € € € €

Quadratic Regularization
Quadratic Regularization (general case)
Interpretation
Model Norm: Measure of
“bad” features

Minimize J =||Wr − d ||22 + µ || Lr ||22

Misfit: Measure of Data Fidelity

€
Trade-off parameter

Part 3 Intro to inversion 38

19
 Quadratic Regularization
J =||Wr − d ||22 + µ || Lr ||22

€4
µ = 10
2 Low resolution / underfitting
||Wr − d || 2
∗
€ ∗
Unstable / overfitting
€ ∗
∗ µ = 10−6
€ ∗
€ || Lr ||22
Trade-off curve / L-curve €
€ Morozov's Discrepancy Principle for Tikhonov regularization

€
Part 3 Intro to inversion
€
39

Quadratic Regularization
Quadratic Regularization
Typical regularizers
J =||Wr − d ||22 + µ || Lr ||22

L=I zero order quadratic regularization (smallest model)

€
L = D1 First order quadratic regularization (flattest model)

L = D2 Second order quadratic regularization (smoothest model)

D1, D2 Are high pass operators (First and Second order Derivatives)

Part 3 Intro to inversion 40

€

20
 "1 0 0 0 0 0 %
Quadratic Regularization $
$0 1 0 0 0 0
'
'
$0 0 1 0 0 0 '
L = I =$ '
Quadratic Regularization $0 0 0 1 0 0 '
$0 0 0 0 1 0 '
$ '
#0 0 0 0 0 1 &

"1 −1 0 0 0 0 0 %
$ '
$0 1 −1 0 0 0 0 '
$0 0 1 −1 0 0 0 '
L = D1 = $ '
$0 0 0 1 −1 0 0 '
$0 0 0 0 1 −1 0 '
$ '
#0 0 0 0 0 1 −1&

"1 −2 1 0 0 0 0%
$ '
$ 0 1 −2 1 0 0 0'
$0 0 1 −2 1 0 0'
L = D2 = $ '
$0 0 0 1 −2 1 0'
$0 0 0 0 1 −2 1 '
$ '
#0 0 0 0 0 1 −2&

Part 3 Intro to inversion € 41

Seismic Deconvolution:
Trying to increase resolution/focusing via mathematical
tricks….

When in Doubt, Smooth

Sir Harold Jeffreys (Quoted by Moritz, 1980; Taken from Tarantola, 1987)

• Or… define a regularization term that can switch

smoothing on/off in order to focus/resolve signals/
models

• Let’s analyze high-resolution estimates obtained

via non-quadratic regularization

Part 3 Intro to inversion 42

21
High-resolution via non-quadratic regularization
(GOAL: Suppress ringing & Increase BW)

J =||Wr − d ||22 + µ R(r)

R(r)= ∑ | ri | L1 norm Non Quadratic Norm

i

2
R(r)= ∑ | ri | L2 norm Quadratic Norm
i

ri2
R(r) = ∑ ln(1+ ) Cauchy norm Non Quadratic Norm
i σ c2

Let’s study the Cauchy norm inversion

Part 3 Intro to inversion 43

High-resolution via non-quadratic regularization

Non-quadratic norms are used to estimate SPARSE models

J =||Wr − d ||22 + µ R(r)

R(r)= ∑ | ri | L1 norm Non Quadratic Norm

i

2
R(r)= ∑ | ri | L2 norm Quadratic Norm
i

ri2
R(r) = ∑ ln(1+ ) Cauchy norm Non Quadratic Norm
i σ c2

Let’s study the Cauchy norm inversion

Part 3 Intro to inversion 44

22
 Seismic Deconvolution: Non-quadratic norms and
sparse deconvolution

q Non-quadratic norms can be used to retrieve sparse models - solution with high frequency
content
q Connected to the retrieval of non-Gaussian models/signals
q Sparse models are associated to broad-band solutions
D. W. Oldenburg et al, 1982, Recovery of the acoustic impedance from reflection
seismograms, Geophysics, vol. 48, No. 10, pp. 1318-1337

q New interest on these methods because of the importance for reconstruction & sparse
coding algorithms:
q Data reconstruction and spectral analysis from short records
Sacchi, M.D., Ulrych, T.J, and Walker, C., 1998, Interpolation and extrapolation using
a high resolution discrete Fourier transform: IEEE Trans. on Signal Processing, 46,
No. 1, 31-38.

q Modeling simple-cell receptive primary visual cortex in mammals

Olshausen, B. A., & Field, D. J. (1996). Emergence of simple-cell receptive field
properties by learning a sparse code for natural images. Nature, 381, 607-609

q Compressed sensing and data recovery

Baraniuk, Compressive Sensing, IEEE Signal Processing Magazine, July 2007)

Part 3 Intro to inversion 45

Cauchy Norm Deconvolution

High-resolution via non-quadratic regularization

J =||W r − d ||22 +µ S(r)

ri2
S(r) = ∑ ln(1+ )
i σ2 Non-Linear system of
equations - can be solved
with IRLS
∇J = W 'W r − W ' d + µ Q(r) r = 0
Iterative Reweighted least
squares
r = (W 'W + µQ(r))−1W ' d

Part 3 Intro to inversion 46

€

23
 Example - High Res. Result with Cauchy regularization
Time

r W d0 d r˜r˜ rq rnq
r˜ : Matched filtering (cross-correlation of wavelet with trace)
rq : Quadratic regularization (l2 norm on reflectivity)
rnq : Non-quadratic regularization (Cauchy norm on the reflectivity)
47
€ €Part 3 Intro€to inversion€ €€ € €
€
€
€
Cauchy Norm Deconvolution

High-resolution via non-quadratic regularization

Assignment:
Derive Q for the Cauchy norm, ell-1 norm, and ell1-ell2 norm and the ell2 norm

Part 3 Intro to inversion 48

24
 Cauchy Norm Deconvolution

Super-resolution via non-quadratic regularization

J =||W r − d ||22 +µ S(r)

ri2
S(r) = ∑ ln(1+ )
i σ2

∇J = W T W r − W T d + µ Q(r) r = 0

Non-Linear system of equations can be

r = (W T W + µQ(r)) −1W T d solved with IRLS (Iterative Reweighted
Least Squares)

Part 3 Intro to inversion 49

€

Cauchy Norm Deconvolution

r = (W T W + µQ(r)) −1W T d
W 'W : Autocorrelation matrix of the wavelet (Toeplitz)

# µQ(r1 ) &
% (
% µQ(r1 ) (
µQ(r) = % ! (
% (
% µ Q(rM −1 ) (
% (
$ µQ(rM )'

r = (W T W + µQ(r)) −1W T d

1
Qii = Q(ri ) = , we now have an adaptive damping
σ 2 (1+ ri2 /σ c 2 )

ri →small, Q(ri ) →1/σ 2 →large if σ is small

ri →large, Q(ri ) →small →if ri >> σ

Part 3 Intro to inversion 50

€

25
 Cauchy Norm Deconvolution

The problem is non-linear, the reflectivity

depends on the reflectivity… IRLS (Iterative Re-weighted Least Squares)
T −1 T iter = 0
r = (W W + µQ(r)) W d
r iter = 0

For iter = 1: iter_max

1
Qii =
σ (1+ ri2 /σ 2 )
2
1
Qiter ii =
σc 2 (1+ (riiter ) 2 /σ 2 )

r iter+1 = (W T W + µQiter ) −1W T d

€
End

€
Part 3 Intro to inversion 51

I used IRLS…

But but . there are 1010 methods algorithms to

solve the l1-l2 problem..

IRLS
ISTA
FISTA = Fast Iterative Soft Threshold Algorithm
Bregmann Iterations
Proximity Methods
ADMM

Part 3 Intro to inversion 52

26
 Cauchy Norm Deconvolution

For large inverse problems, IRLS requires

the following modification: the inversion IRLS (Iterative Re-weighted Least Squares)
stage is done with a semi-iterative solver :
iter = 1
• Conjugate Gradients (CG)
r iter = 0
• Gauss-Seidel
• Pre-conditioned CG For iter = 1: iter_max

• LSQR 1
Qiter ii =
Iterative solvers can compute an σ 2 (1+ (riiter ) 2 /σ 2 )
approximation to the solution of the linear
system r iter+1 = (W 'W + µQiter )−1W ' d
• HR Radon uses the above trick
End

€
Part 3 Intro to inversion 53

R(r)= ∑ | ri |
i

2
R(r)= ∑ | ri |
i

ri2
R(r) = ∑ ln(1+ )
i σ2

R(r) = ∑ F(ri )
i

General Form of the

€ norm

Part 3 Intro to inversion 54

27
 R(r) = ∑ F(ri )
i

dF(x)
: Influence Function
dx

Part 3 Intro to inversion 55

• L2 norm: all samples are weighted the

same

r = (W 'W + µQ(r))−1W ' d • Cauchy norm: the weight is model

dependent. This is what allows the
solution to become sparse
1
Qii =
σ 2 (1+ ri2 /σ 2 )
−1
Qii

µ Qii : Model Dependent

Pre - whitening €

Qii−1 : Can be intepreted as

model - dependent variance

Part 3 Intro to inversion 56

28
Cauchy Norm Deconvolution:
Algorithm [IRLS]
True reflectivity Trace
r = zeros(N,1);!
sc=0.01; !
mu =0.01!
iter_max = 10;!
!
R = W'*W;!
iter_max;!
!
for k=1:iter_max;!
Estimated reflectivity
!
Q = diag(1./(sc^2+r.^2));!
Matrix = R + mu*Q;!
r = inv(Matrix) * W'*s;!
!
end

Part 3 Intro to inversion 57

Cauchy Norm Deconvolution:

Algorithm [IRLS]
True reflectivity Trace
r = zeros(N,1);!
sc=0.01;!
mu = .01!
iter_max = 10;!
R = W'*W;!
!
for k=1:iter_max;!
!
Q =diag(1./(sc^2+r.^2));! Estimated reflectivity
Matrix = R + mu*Q;!
r = inv(Matrix) * W'*s;!
!
sc = 0.01*max(abs(r));!
!
Adaptive
end;!
version
!
!

Part 3 Intro to inversion 58

29
Cauchy Norm Deconvolution:
Algorithm for real data … (HFRestoration or Cauchy
Norm Inversion)

Similar to MATLAB prototype but in f77

• Pre-conditioning is used to avoid specifying the trade-off parameter
(mu)
• The number of iteration of the external loop is used as trade-off
• Linear inversion is done in the flight with CG (the internal loop)
• Pre-conditioning makes it robust to parameter selection (you do not
have to play to much with trade-off parameters!)
• We will briefly discuss the pre-conditioning problem when dealing with
de-migration/migration operators

Part 3 Intro to inversion 59

High frequency imaging methods….

You can use any norm that drives the algorithm to

sparse solutions (Ell-1, Ell-p with p close to 1, etc)

DEBEYE, H. W. J. & RIEL, P. (1990) Lp-NORM DECONVOLUTION. Geophysical

Prospecting 38 (4), 381-403.

Other methods exits - they all attempt to retrieve a sparse reflectivity sequence:
Atomic Decompostion/Matching Pursuit / Basic Pursuit (like an L1)

Chopra, S., J. Castagna, and O. Portniaguine, 2006, Seismic Resolution and Thin-Bed Reflectivity
Inversion: Recorder, 31, 19-25. ([Link])

Portniaguine, O., and J. P. Castagna, 2004, Inverse spectral decomposition: 74th Annual
International Meeting, SEG, Expanded Abstracts, 1786-1789.

Part 3 Intro to inversion 60

30
 3D Seismic Volume (Venezuela)
≈ 20 km

t = 2z /v
t ∝z

t = 3s
Time

< v >= 2000m /s

z = 3km

Conventional Deconv High Res. Deconv

Time: Time that a wave takes to go to and come back from a layer
Red: Positive polarity Blue: negative polarity - The image portrays some sort of reflectivity

Part 3 Intro to inversion 61

Conventional Deconv High Res. Deconv

Part 3 Intro to inversion 62

31
 Amplitude spectrum

Can you trust the new bandwidth? I am not saying this is the true answer but a possible
answer if you believe in the sparse reflectivity assumption.

Part 3 Intro to inversion 63

Cauchy Norm Deconvolution:

with impedance constraints

Ik +1 − Ik
rk = Reflectivity as a function of P-impedance
Ik +1 + Ik

k
1 Log Approximation (Logarithm impedance and
ξ k = log(Ik /I0 ) ≈∑ rj
2 j= 0
reflectivity are related by a linear operator)

J =||W r − d ||22 +µ S(r) + β || Cr − ξ ||22

Fit the data Fit impedance constraints

Solution must be sparse (High freq)

€
Part 3 Intro to inversion 64

32
 Cauchy Norm Deconvolution:
with impedance constraints
True reflectivity Trace

β =0
Constraints are
not honored

Inverted impedance &

true impedance (red)

Part 3 Intro to inversion 65

Cauchy Norm Deconvolution:

with impedance constraints
True reflectivity Trace

β = 0.25

Constraints are
€ honored

Inverted impedance &

true impedance (red)

Part 3 Intro to inversion 66

33
Cauchy Norm Deconvolution:
with impedance constraints - Algorithm [IRLS + soft
bounds]

beta=0.25!
iter_max=10!
R = W'*W+beta*C'*C;!
r = zeros(N,1);!
!
for k=1:iter_max;!
Q = diag(1./(sc^2+r.^2));!
Matrix = R + mu*Q;!
r = inv(Matrix) * (W'*s+beta*C'*psi)!
sc = 0.01*max(abs(r));!
end;!
!

Part 3 Intro to inversion 67

Angle dependent reflectivity

AVO: Amplitude versus offset

AVA: Amplitude versus angle
AVP: Amplitude versus Ray Parameter

Estimation of rock
properties & HCI

Part 3 Intro to inversion 68

34
 AVO as a Multi-channel Seismic Deconvolution
Problem
Offset (h)
w(t) d(t,α , x)

α α
€ € Vp1,Vs1, ρ1

€ € Vp2 ,Vs2 , ρ 2
€
x
€
d(t,α , x) = ∫ r(τ ,α, x) w(τ − t)dτ

€
r(t,α , x) = f (Vp1,Vs1, ρ1,Vp2 ,Vs2 , ρ 2 ) = A(t, x) + B(t, x) × F(α )
Shuey’s approximation

Part 3 Intro to inversion 69

€

AVO as a Multi-channel Seismic Deconvolution

Problem

d(t,α , x) = ∫ r(τ ,α, x) w(τ − t)dτ

r(t,α , x) = f (Vp1,Vs1, ρ1,Vp2 ,Vs2 , ρ 2 ) = A(t, x) + B(t, x) × F(α )

A 
d ≈ Lm, m= 
 B

J =|| d − Lm ||22 + µR(m)

€ Part 3 Intro to inversion 70

35
Real data example - comparison of quadratic vs non-
quadratic regularization

• Quadratic - damped least-squares

• Non-quadratic - Cauchy regularization + Internal iteration with Conjugate Gradients

• To stress: The two solutions that I will show fit the data at the same level

Part 3 Intro to inversion 71

AVO as a Multi-channel Seismic Deconvolution

Gather

Offset Gather

Wavelets
Offset

Far Offset Wavelet

Near Offset Wavelet
Spectra

Wavelets estimated from first gather

Data courtesy Xin-Gong Li
(Techseis)

Part 3 Intro to inversion 72

36
 A(t,x) B(t,x)

Quadratic regularization solution

Part 3 Intro to inversion 73

A(t,x) B(t,x)
Nonquadratic regularization solution
• More details: Dey, A.K., M.D. Sacchi, and A. Gisolf, 2006, High-resolution reservoir rock properties via joint prestack seismic amplitude inversion:
76th Ann. Internat. Mtg.: SEG, Expanded Abstracts, 2156-2160.

Part 3 Intro to inversion 74

37
Data fitting for first gather
First Gather / Fit

Input
Prediction
Residuals

Quadratic regularization solution

Part 3 Intro to inversion 75

Data fitting for first gather

First Gather / Fit

Input
Prediction
Residuals

Nonquadratic regularization solution

Part 3 Intro to inversion 76

38
Remarks and some new directions

• Tools for reflectivity inversion, impedance inversion, AVO etc etc are quite similar

• For extensions to 2D-3D one should add constraints for lateral continuity
– Juefu Wang et al, 2007, High Res Deconvolution with sparsity and lateral
continuity constraints ([Link]
[Link])

• Simultaneous AVO deconvolution and wavelet estimation:

– Inversion for source wavelet and AVA parameters form prestack data, Routh,
Anno, Baumel and Chavarra, SEG 2003, Abstracts

• Including correlation between sparse estimators of intercept and gradient

– See Jon Downton PhD dissertation, CREWES-UofC.

Part 3 Intro to inversion 77

Juefu Wang et al., CSEG 2007

WCSB

Part 3 Intro to inversion 78

39
 Large Scale Problems

Inverse problems -
Angle dependent reflectivity and
LS migration

Part 3 Intro to inversion 79

Large inverse problems: Pre-stack depth

Imaging

Imaging Angle Dependent Reflectivity

Wang J., Kuehl H. and Sacchi M.D., 2005, High-resolution wave-equation AVA
imaging: Algorithm and tests with a data set from the Western Canadian Sedimentary
Basin: Geophysics, 70, 891-899

Sacchi M.D., Wang J., and Kuehl H., 2006, Regularized migration/inversion: new
generation of imaging algorithms, CSEG Recorder, Volume 31, Special Edition, 54-59.

Wang J. and Sacchi M.D., 2007, High-resolution wave equation AVP imaging with
sparseness constraints, Geophysics, 72 (1), S11-S18.

Part 3 Intro to inversion 80

40
 Imaging - Large scale problems

m
h x

α €
m(α, x,z)

€ €

m: midpoint
h: offset
Wang, Kuehl and Sacchi, Geophysics 2005
Part 3 Intro to inversion 81

Imaging in Operator Form..

Lm≈d

sin(α )
m = m(x,z, p), p=
v (x,z)

d = d(g,s,t)

d: is now pre-processed seismic data

€
Part 3 Intro to inversion 82

41
 Migration/Inversion Hierarchy..

Lm = d Modeling (De-migration operator)

mmig = LT d Migration

€
minv = L−1d Inversion

Part 3 Intro to inversion 83

Regularized Imaging

J=|| W(Lm − d) ||22 + µ R(m)

Remarks:
€
• We don’t have L, we only have a code that knows how to apply L to m and,
another, that knows how to apply LT to d
• L and L’ are built to pass the dot product test
• Importance of sampling Matrix W for reducing footprints
Kuehl H. and Sacchi M.D., 2003, Least-squares wave-equation migration
for AVP/AVA inversion: Geophysics, 68, 262-273

Part 3 Intro to inversion 84

42
Regularized Imaging
m Function L d
Dot Product Test d Function LT m’
(only for programmers..)

1) m1, d1 = L m1

2) d2, m2 = LT d2

3) if d1T d2 = m1T m2 ,

then, LT is the adjoint of L

Goal: Guarantee adjointness of L and LT (numerically)

Part 3 Intro to inversion 85
€

Regularized Imaging with Conjugate Gradients

J=|| W(Lm − d) ||22 + µ F(m)

m = m(x,z, p)

F : Smoothness in x, y, p ,
Sparseness in z

€[Link], PhD 2005, UofA, [Link]

Part 3 Intro to inversion 86

43
Regularized Imaging with Conjugate Gradients
+ Preconditioning

Reduce computational cost

J=||W (Lm − d) ||22 +µ || Rm ||22

J' =||W (L Pν − d) ||22 + µ|| ν ||22

m =Pν , P ≈ R −1

Part 3 Intro to inversion 87

Regularized Imaging with Conjugate Gradients

+ Preconditioning + Regularization by Iteration

J' =||W (L Pν − d) ||22 + µ|| ν ||22

k : CG Iteration

J'' =||W (LPν k − d) ||22 < tol

Hansen, 1987, Rank-Deficient and Discrete Ill-Posed Problems:

€ Aspects of Linear Inversion (SIAM Monographs on
Numerical
Mathematical Modeling and Computation)

Part 3 Intro to inversion 88

44
 Examples:
• Regularized Imaging to attenuate the influence of missing data
• Increase resolution of seismic images

Part 3 Intro to inversion 89

Marmousi model (P-wave velocity and density)

Horizontal distance (km)

3 4 5 6 7 8
0

5300
Depth (km)

4800 1
4300
m/s

3800
3300
Test AVA Target
2800
2300 2 Horizontal distance (km)
1800 3 4 5 6 7 8
1300 0

3
1.0
1
Depth (km)
kg/m^3

0.8

0.6
2
0.4
x10 4

Part 3 Intro to inversion 90

45
 CMP position [m] Ray parameter [mu s/m]
7000 7500 8000 0 200 400 600 800
0 0
Migration
(Complete data)
0.5 0.5

1.0 1.0
Depth [km]

Depth [km]
1.5 1.5

2.0 2.0

2.5 2.5

Part 3 Intro to inversion 91

CMP position [m] Ray parameter [mu s/m]

7000 7500 8000 0 200 400 600 800
0 0

Least Squares
Migration with
Smoothing
regularization along 0.5 0.5
ray parameter p
(Complete data)

1.0 1.0
De pth [ km]

De pth [k m]

1.5 1.5

2.0 2.0

2.5 2.5

Part 3 Intro to inversion 92

46
 Ray parameter [mu s/m] Ray parameter [mu s/m]
0 200 400 600 800 0 200 400 600 800
0.5 0.5

0.7 0.7

D ept h [km]
Dept h [k m]

0.9 0.9

1.1 1.1

0.2 0.2
Reflection coefficient

0.18 0.18

Reflection coefficient
0.16 0.16
0.14 0.14
0.12 0.12
0.1 0.1
0.08 0.08
0.06 0.06
0.04 Migration 0.04 LSM
0.02 0.02
00 10 20 30 40 50 00 10 20 30 40 50
Angle of incidence (deg) Angle of incidence (deg)

Part 3 Intro to inversion 93

Offset [m] Offset [m]

200 400 600 800 1000 1200 200 400 600 800 1000 1200
0 0

Complete Incomplete
0.5 CMP 0.5 CMP (30 %)

1.0 1.0
Ti m e [sec ]

Ti m e [ sec ]

1.5 1.5

2.0 2.0

2.5 2.5

All the data (240 shots) were decimated

Part 3 Intro to inversion 94

47
 Ray parameter [mu s/m] Ray parameter [mu s/m]
0 200 400 600 800 0 200 400 600 800
0.5 0.5

0.7 0.7

D ept h [km]
Dept h [k m]

0.9 0.9

1.1 1.1

0.2 0.2
0.18
Reflection coefficient

0.18

Reflection coefficient
0.16 0.16
0.14 0.14
0.12 0.12
0.1 0.1
0.08 0.08
0.06 0.06
0.04 0.04
0.02
Migration 0.02 LSM
00 10 20 30 40 50 00 10 20 30 40 50
Angle of incidence (deg) Angle of incidence (deg)

Part 3 Intro to inversion AVP from Incomplete data 95

Field 3D
ERSKINE data
Dataexample
/ WCSB /Alberta

ERSKINE orthogonal 3-D sparse land data set

dx=50.29m

dy=33.5m

Part 3 Intro to inversion 96

48
 CIG at x-line #10, in-line #71

Migration CG 4 CG 11 PCG 4
Comparison of LSM solution with regularization in the ray parameter direction. CG is the solution using
conjugate gradients, PCG is the solution with Pre-conditioned CG.

Part 3 Intro to inversion 97

Migration

RLSM (11)

PLSM (4)

Position of the CIG gather CIG / AVP Gather

Part 3 Intro to inversion 98

49
 Wang J., Kuehl H. and Sacchi M.D., 2005, High-
resolution wave-equation AVA imaging: Algorithm and
tests with a data set from the Western Canadian
Sedimentary Basin: Geophysics, 70, 891-899.

VP = (VS −1360) /1.16 m /s (Castagna, 1985)

The analysis is restricted to the Ellerslie (sandstone)
and Banff (shales) formations.
The economical target, the Leduc Reef, is
excluded from the AVA analysis due the
difficulty of estimating carbonate shear velocities.

Part 3 Intro to inversion 99

Synthetic trace

Part 3 Intro to inversion 100

50
 CIG (synthetic vs. inverted) AVA (synthetic vs. inverted)

Red: synthetic
Blue: inverted

Synthetic CIG is computed

using Aki & Richards’ ap-
proximation

A: synthetic, B: inverted

Part 3 Intro to inversion 101

Large scale Imaging - Including Resolution

Enhancement (using Cauchy norm R )

Cost function for SLSM

J(m) =|| W(Lm − d) ||22 + µ R(SH(m))

H : High Pass operator along p - gathers

S : Stacking on p − gathers

m = Pz

Preconditioning implementation
2
J'(m) =|| W(LPz − d) ||2 + µ R(Sz)

P is chosen to behave like H -1

Part 3 Intro
€to inversion 102

51
 Marmousi data / Common image gathers

Part 3 Intro to inversion 103

Erskine data /Common image gathers

LS Migration +
Migration LS Migration Sparseness

Part 3 Intro to inversion 104

52
 Imaging & High Resolution Imaging

Imaging with Adjoint Operators [Current technology]

Imaging + Quadratic Regularization * [Not in production]

Imaging + Non-Quadratic Regularization ** [??]

Resolution

Part 3 Intro to inversion 105

Addendum

The Bayesian Framework

Part 3 Intro to inversion 106

53
A bit about Bayes

Thomas Bayes was born in London in 1702 into a religious atmosphere. His, father, the
Rev. Joshua Bayes, was one of the first six Nonconformist ministers to be ordained in
England. Like his father, Thomas was ordained a Nonconformist minister and assisted
his father until the late 1720s when he become a Presbyterian minister.
Bayes theory appeared posthumously in “Essay Towards Solving a Problem in the
Doctrine of Chances” published in the Philosophical Transactions of the Royal Society of
London in 1764. Thomas Bayes died in England in 1761.

From: A Bayes tour of Inversion: A tutorial. Ulrych, Sacchi and Woodbury, Geophysics, 2001, 66-1, p.55-64.

Part 3 Intro to inversion 107

Cost functions derived from Bayes rule

• Bayes rule
• MAP solution
• Cost functions induced by a Posterior
• Examples

Part 3 Intro to inversion 108

54
Goal of this section

• Make the connection between probabilities and regularization

• Laplace distribution à l1 norm
• Gauss distribution à l2 norm
• Cauchy distribution à Cauchy criterion

Part 3 Intro to inversion 109

Bayes theorem

p(A | B) p(B) = P(A) p(B | A)

p(A) p(B | A)
p(A | B) =
p(B)

p(A) : Prior
p(B | A) : Probability of B given A
p(A | B) : Posterior

Update rule for

probabilities = state of information

Part 3 Intro to inversion 110

€

55
 Relation to Inverse problems

• Consider a linear problem d = Lm + n

p(A) p(B | A)
• Apply Bayes rule p(A | B) =
p(B)
€
p(A) = p(m) Prior
p(B | A) = p(d | m) Likelihood
p(A | B) = p(m | d) Posterior

p(m) p(d | m)
p(m | d) =
p(d)

Part 3 Intro to inversion 111

€

Relation to Inverse Problems

∫ p(m | d)dm = 1 ⇒ p(d) = ∫ p(m) p(d | m)dm

p(m) p(d | m)
p(m | d) =
p(d)

p(m | d) ∝ p(m) p(d | m)

Posterior ∝ Prior × Likelihood

Part 3 Intro to inversion 112

€

56
Example: Gaussian Noise and “Gaussian Model”
1
p(d | m) ∝ exp[− (d − Lm)T Cd −1 (d − Lm)]
2

1 −1
p(m) ∝ exp[− (m − m0 )T Cm (m − m0 )]
2

m0 = 0, Cm = σ m2 I, Cd = σ d2 I

p(m | d) ∝ p(m) p(d | m)

1 1
p(m | d) ∝ exp[−( 2
mT m + 2 (d − Lm)T (d − Lm))]
2σ m 2σ d

€ 3 Intro to inversion
Part 113

Example: Gaussian Noise and “Gaussian Model”

Maximze
1 1
p(m | d) ∝ exp[− mT m − 2 (d − Lm)T (d − Lm)] = exp[−J]
2σ m2 2σ d

Minimize

1 1
J= mT m + 2 (d − Lm)T (d − Lm)
2σ m2 2σ d
or
σ d2 T
Cost = m m + (d − Lm)T (d − Lm)
σ m2
= µmT m + (d − Lm)T (d − Lm)

€
Part 3 Intro to inversion 114

57
 Statistical Interpretation of the trade-off parameter

σ d2
µ= 2
σm

Noise then µ
Noise then µ
€
€
Part 3 Intro to inversion 115

MAP Solution (Maximum a Posteriori)

Minimize
σ d2 T
Cost = 2 m m + (d − Lm)T (d − Lm)
σm
= µmT m + (d − Lm)T (d − Lm)
⇒
m = (LT L + µI)−1 LT d

€
Part 3 Intro to inversion 116

58
Remarks

• Today I’ve used the idea of regularization methods to construct a cost function
• Now I will derive the cost function for quadratic and non-quadratic regularization and
quadratic misfit using Bayes theorem:

• Probability of data given model is the probability of the noise and we assume
Gaussian noise. The assuption will lead to a Quadratic Misfit

• If we model our state of information about the model is encoded in a Gaussian

pdf then the regularization term is Quadratic

• Other priors are possible, they all lead to non-quadratic regularization

methods

Part 3 Intro to inversion 117

More on Bayes

Assume that the parameters to estimate (model) are non-Gaussian (say the follow a
Laplace or a Cauchy distribution)

p(x) ∝ exp[− λ | x |p ]

p(x) ∝ exp[− λ | x |] Laplace p=1

€
1
p(x) ∝ Cauchy

x2
1+ 2
σ

Part 3 Intro to inversion 118

59
 Tails

Part 3 Intro to inversion 119

For instance use Cauchy

Maximize
1 1
p(m | d) ∝ ∏ × exp[− 2 (d − Lm)T (d − Lm))]
i 1+ mi2 /σm2 2σd

Minimize
Cost = µ∑ ln(1+ mi2 /σm2 ) + (d − Lm)T (d − Lm)
i

= µ∑ ln(1+ mi2 /σm2 )+ | d − Lm |22

i

Cauchy “norm”

€

Part 3 Intro to inversion 120

60
For instance, let’s use the Laplace prior

Maximize
1
p(m | d) ∝ ∏ exp(− λ | x i |) × exp[− (d − Lm)T (d − Lm))]
i 2σ d2

Minimize
Cost = µ∑ | mi | + (d − Lm)T (d − Lm)
i

= µ | m |1 + | d − Lm |22

Ell-1 norm

€

Part 3 Intro to inversion 121

Example (Non-Gaussian reflectivity) retrieval with

sparse inversion

• We assume a non-Gaussian reflectivity and known wavelet

SPARSENESS

Part 3 Intro to inversion 122

61
 Mixing Parameter p=0.8 (very sparse)
Cauchy inversion

Data

True Impulse
Response

Predicted
Data

Estimate Impulse
Response

Part 3 Intro to inversion 123

Mixing Parameter p=0.2

Cauchy inversion

Data

True Impulse
Response

The solution is too Predicted

simple (sparse).
Data

Estimate Impulse
Response

Part 3 Intro to inversion 124

62
Kurtosis

K=3

Error = difference between true and estimated impulse response (reflectivity)

Results obtained after 20 realizations

Part 3 Intro to inversion 125

Remark

• Key feature for proper recovery of the impulse response with sparse inversion..
• Sparsity (of course)….. But something else is needed

Part 3 Intro to inversion 126

63
Remark

• Key feature for proper recovery of the impulse response with sparse inversion..
• Sparsity (of course)….. But something else is needed

MEASURABLE INFORMATION…. Which is controlled by the

source band-width

Part 3 Intro to inversion 127

Source BW vs Estimation error for a sparse reflectivity

(p=0.5)

Error = difference between true and estimated impulse response

Part 3 Intro to inversion 128

64
Source functions used in the simulation

Part 3 Intro to inversion 129

• BW is not only important for de-phasing (e.g. maximum Kurtosis de-phasing) but also
for the inversion of the reflectivity
• Nothing new: more BW means that more information about the reflectivity is preserved
in the seismic trace.

– In the absence of sufficient BW, the prior will dominate the solution and dictate
how the solution looks like. The geologist/interpreter (with access to borehole-
derived reflectivity) will feel deeply discontented with the simplicity of your model.

Part 3 Intro to inversion 130

65