Geophysical Inversion:

which model do you want?

Steven Constable
Scripps Institution of Oceanography
 Acknowledgements

SEG and SEG Foundation

Sponsored by Statoil

Sponsored by Paradigm

Scripps Institution of Oceanography,

Seafloor Electromagnetic Methods Consortium
 SEG Membership Benefits

SEG Digital Library – full text articles

Technical Journals in Print and Online
y !
To d a Networking Opportunities
Join
Membership Discounts on
! Continuing Education Training Courses
! Publications (35% off list price)
! Workshops and Meetings

Membership materials are available today!

Join Online at [Link]/About-SEG/Membership

 Student Opportunities

❑ Sponsored Membership
❑ Student Chapter Programs
❑ SEG/Chevron Student Leadership Symposium
❑ SEG/ExxonMobil Student Education Program
❑ Challenge Bowl
❑ Scholarships
❑ Field Camp Grants
❑ Geoscientists Without Borders®
❑ Student Expos & IGSCs
❑ Honorary and Distinguished
Lecturers
❑ SEG Online

Learn more at [Link]/Education/Students-Early-Career

 Section/Associated Society
Opportunities

❑ Host DL, HL, and DISC Programs

❑ Council Representation
❑ Annual Meeting Booth Discount
❑ Best Papers presented at SEG
Annual Meetings
❑ Joint Conferences, Workshops, and
Forums in partnership with SEG

For more information, including a list of benefits, please visit:

[Link]/resources/sections-societies
 Geophysical Inversion:
which model do you want?
Steven Constable
Scripps Institution of Oceanography
There is an old joke…
The production manager asked a geologist, engineer,
and geophysicist what 2 + 2 was.
The production manager asked a geologist, engineer,
and geophysicist what 2 + 2 was.

The geologist thought for a bit and then said

“somewhere between 3 and 5”.
The production manager asked a geologist, engineer,
and geophysicist what 2 + 2 was.

The geologist thought for a bit and then said

“somewhere between 3 and 5”.

The engineer fiddled with a calculator and said

“3.9999999”.
The production manager asked a geologist, engineer,
and geophysicist what 2 + 2 was.

The geologist thought for a bit and then said

“somewhere between 3 and 5”.

The engineer fiddled with a calculator and said

“3.9999999”.

The geophysicist looked her in the eye and asked

“what answer do you want”
It is notable that when once searching
the web for this joke, I not only got the
joke page of an oil price blog, but also
the Wikipedia entry for “Inverse
problem”.

This says it all…

Today we will explore how the inverse

problem can give you whatever model
you want.
It is notable that when once searching
the web for this joke, I not only got the
joke page of an oil price blog, but also
the Wikipedia entry for “Inverse
problem”.

This says it all…

Today we will explore how the inverse

problem can give you whatever model
you want.

Well, almost.
With seismic reflection images you can often see the geology in the data.

For EM and potential field methods you need inversion to recover something that
can be interpreted as geology.

Same is true for seismic tomography and full waveform inversion.

Moore et al., 2007, Science.

 −10

0.25 Hz These are some

E−field amplitude, log(V/A m 2 )

−11 marine controlled

source electromagnetic
−12
CSEM data. You can’t
say much about
−13
geology just by looking
at them.
−10

0.75 Hz
E−field amplitude, log(V/A m 2 )

−11

−12

−13

−10

1.75 Hz
E−field amplitude, log(V/A m 2 )

−11

−12

−13

10 12 14 16 18 20 22 24 26
Transmitter position, km
Constable, Orange, and Key, 2015
 The same is true of magnetotelluric (MT) data…
site 1 site 2 site 3 site 4 site 5 site 6 site 7 site 8
Log rho 0.5

-0.5
50
Phase

40

30

0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2
log period, s

site 9 site 10 site 11 site 12 site 13 site 14 site 15 site 16

0.5
Log rho

-0.5
50
Phase

40

30

0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2
log period, s

Constable, Orange, and Key, 2015

With many modern inversion algorithms available, it is all-so-easy to input data
and turn the crank to get a model.

geology Data

Data errors
and misfit
Geophysical
inversion model
algorithm
Regularization

Priors/ Model
constraints Parameterization
With many modern inversion algorithms available, it is all-so-easy to input data
and turn the crank to get a model.
One of the main messages of my talk is that models from geophysical inversion
depend on much more than the data inputs:

geology Data

Data errors
and misfit
Geophysical
inversion model
algorithm
Regularization

Priors/ Model
constraints Parameterization
With many modern inversion algorithms available, it is all-so-easy to input data
and turn the crank to get a model.
One of the main messages of my talk is that models from geophysical inversion
depend on much more than the data inputs:

By the way… this still

Data isn’t geology - it is
geology
density, conductivity,
seismic velocity, etc.

Data errors
and misfit
Geophysical
inversion model
algorithm
Regularization

Priors/ Model
constraints Parameterization
Forward modeling:

data
model space space

d̂ = f (x, m) Some forward functional f

m = (m1 , m2 , ....., mN ) Model parameters (layers, blocks, ...)
x = (x1 , x2 , x3 , ......, xkM ) Independent variables (freqs., locations, ...)
d̂ = (dˆ1 , dˆ2 , dˆ3 , ....., dˆM ) Predicted data (gravity, magnetic, electric, ...)
Inverse modeling:

model space data

space

Given real (observed) data d = (d1 , d2 , d3 , ....., dM )

with errors = ( 1 , 2 , ..., M )
find an m
There are several approaches to inversion:

Stochastic
Monte Carlo, Markov Chains
Genetic Algorithms
Simulated annealing, etc.
(Bayesian Searches)

Deterministic
Newton Algorithms
Steepest descent
Conjugate Gradients
Quadratic (and Linear) Programming, etc.

Analytical
D+ (1D MT)
Bilayer (1D resistivity)
Ideal body theory in gravity and magnetism
Stochastic methods: “acceptable” data
space

model space

data
space

A useful approach, largely restricted to simple problems (because millions of

models required), with most of the subtlety in model generation methods.

The advantages are that (i) only forward calculations are made and (ii) some
statistics can be obtained on model parameters. Best for sparsely parameterized
models. One needs to be careful that bounds on explored model space don’t
unduly influence the outcome.
 Deterministic
Newton Algorithms
Steepest descent
Conjugate Gradients
“acceptable” data
model space space

data
space

starting model
The direction of the search is determined by how changing parts of the model
affects the fit to the data.
Analytical
e.g. D+ (1D MT) and Bilayer (DC resistivity)

model space
my data

best
fitting model
(guaranteed!)
These solutions are guaranteed best fitting but pathological.

Resistivity
MT

Bilayer
D+

MT: Delta functions of Resistivity: Arbitrarily thin surface layers

We don’t know for sure, but least squares (LS) fits to higher dimensional models
are probably also pathological.

But we are pretty sure that true LS solutions are maximally “rough”.
What we have talked about so far is model construction. For a great many
geophysicists this is what they think of when inversion is mentioned. More rigorous
approaches try to obtain bounds on model properties - something that is true of all
models. The classic example is total mass from gravity:
So what constitutes an “adequate” fit to the data?

geology Data

Data errors
and misfit choice
Geophysical
inversion σ model
algorithm
Regularization

Priors/
Model
constraints
Parameterization
For noisy data (read: all data), we need a measure of how well a given model fits.
Sum of squares is the venerable way:
M
⇤ 1 ⇥2
⇥2 = 2 di f (xi , m)
i=1 i

or 2
= ||Wd Wd̂|| 2

where W is a diagonal of reciprocal data errors

W = diag(1/ 1 , 1/ 2 , ...., 1/ M ) .
For noisy data (read: all data), we need a measure of how well a given model fits.
Sum of squares is the venerable way:
M
⇤ 1 ⇥2
⇥2 = 2 di f (xi , m)
i=1 i

or 2
= ||Wd Wd̂|| 2

where W is a diagonal of reciprocal data errors

W = diag(1/ 1 , 1/ 2 , ...., 1/ M ) .

I like to remove the dependence on data number and use RMS:

p
RMS = 2 /M .
2
The instinctive approach as this point is to minimize .

This is least squares.

For noisy data (read: all data), we need a measure of how well a given model fits.
Sum of squares is the venerable way:
M
⇤ 1 ⇥2
⇥2 = 2 di f (xi , m)
i=1 i

or 2
= ||Wd Wd̂|| 2

where W is a diagonal of reciprocal data errors

W = diag(1/ 1 , 1/ 2 , ...., 1/ M ) .

I like to remove the dependence on data number and use RMS:

p
RMS = 2 /M .
2
The instinctive approach as this point is to minimize .

This is least squares. In geophysics, this is dangerous!

Why is this dangerous? Because as you try to approach the LS solution, your
model tries to approach the maximally rough, pathological LS solutions, even if
your model space does not contain delta functions, etc.
Existence and Uniqueness: Is there a solution to the inverse problem? Is
there only one solution?

Finite noisy data for a linear problem (say, gravity)

An infinite number of solutions fit the data

Finite noisy data for a nonlinear problem

Either zero or an infinite number of solutions fit the data

Infinite exact data

A unique solution has been shown to exist for a
few cases. Probably true in general but … who cares?

Some people think that we can approach infinite exact data with LOTS of
VERY GOOD data. This is wrong. As Sven Treitel puts it, there is no such
thing as being a little bit non-unique.
Geophysical inversion is non-unique:

2
model space

misfit
space

A single misfit will map into an infinite number of models (or none at all!).
 It is also usually poorly constrained:
2

misfit
space
model space
2
A small distance in corresponds
to a large distance in m

2
(And don’t forget: the minimum
is likely outside your model
parameterization).
So what constitutes an adequate misfit?
2
For zero-mean, Gaussian, independent errors, is chi-squared distributed with M
degrees of freedom. The expectation value is just M, which corresponds to
RMS=1, and so this could be a reasonable target misfit. Or, one could look up the
95% (or other) confidence interval for chi-squared M.

2
for 14 data. For large data
RMS = 1 sets, RMS=1 and RMS95% are
RMS = 1.36 very much the same.

We could use other measures of fit, but the quadratic measure works with the
mathematics of minimization, and for Gaussian errors has nice statistical properties
(unbiased, maximum likelihood, minimum variance). But...
... sum-squared misfit measures are unforgiving of outliers:

-8
10 (marine CSEM data)
With 5% error bars this
data point has the same
-10
weight as 40,000 other data
10
Amplitude, V/Am^2

-12
10

-14
10

-16
10

-10000 -5000 0 5000 10000

Range, m

With Gaussian noise, the probability of a data point being misfit by 6 error bars is
about one in a billion.

All through any inversion process you should monitor weighted residuals
to ensure that there are no bad guys out there.
 Misfit by Period Misfit by Data Type

2 A 1.6 B It is also a good idea to

1.4 look at how the misfit is
1.5 1.2
partitioned across the

Local RMS
Local RMS

1
data:

TM log10(app. res.)
TE log10(app. res.)
1 0.8

Im(Tipper)
Re(Tipper)

TM phase
TE phase
0.6
0.5 0.4 Ideally it should be
0.2 random, but in practice
0 0
10
−1
10
0
10
1
10
2
10
3
10
4
Data Type
very rarely is.
Period (s)

3.5
Misfit by Site Example is from MT
3 C data.
2.5
Local RMS

2
1.5
1
0.5
0
M40 M35 M33 M27 M23 M18 M12 M10 M03 M01 631 633 634 667 668 636 671 672 639
M36 M34 M31 M25 M19 M17 M11 M06 M02 630 632 665 666 635 669 670 637 638 640

PhD thesis, Brent Wheelock, 2012.

Errors come from

• statistical processing errors (spectral estimation for MT; stacking for CSEM and
Phase Degrees
lots of other methods) −10 −5 0 5 10

250 300
2−4 km N = 1646 2−4 km N = 1646
• 200
systematic errors such as navigation
250 errors and instrument calibrations, and
m= 0.3
s = 2.8
m= 1.3
s = 2.4
200
150
Count

• “geological noise” (our inability

150 to parameterize fine details of geology).
100
100

In practice,
50 we only have a good50handle on processing errors - everything else is
lumped0 into a noise floor, which can
0 be pretty arbitrary at times.
−20 −10 0 10 20 −20 −10 0 10 20
Amplitude % Difference Phase % Difference

Phase Degrees
−10 −5 0 5 10
120
4−6 km N = 1491 4−6 km N = 1491
100 m= 0.6 150 m= 1.4
s = 6.3 s = 4.1
80 “errors” computed from
100 actual repeat tows for
Count

60

40
marine CSEM
50
20

0 0
−20 −10 0 10 20 −20 −10 0 10 20
Amplitude % Difference Phase % Difference
(modified from Myer et al., 2012)
 So we are often left without statistical guidance and have to use judgement
in determining an adequate fit. Some people like trade-off, or “L”-curves...

2.5 10
2.4
A 9
8.75
B
8

2.0 7
2

RMS misfit
RMS misfit

1.7
RMS 1.4 5

4
1.5 1.5

3
1.3 2.4
2 2.0
1.2 1.7 1.5
1.3 1.2 1.1
1.1
1
1
0 100 200 300 400 500 600 700 800 0 100 200 300 400 500 600 700 800
Roughness measure Roughness measure
 … but I am not one of them.

2.5 10
2.4
A 9 B
8.75
starting half-space
8

2.0 7
2

RMS misfit
RMS misfit

1.7
RMS 1.4 5
RMS 1.7
4
1.5 1.5

3
1.3 2.4
2 2.0
1.2 1.7 1.5
1.3 1.2 1.1
1.1
1
1
0 100 200 300 400 500 600 700 800 0 100 200 300 400 500 600 700 800
Roughness measure Roughness measure
 In fact you can get pretty well what you want simply by changing the range
of the plot and the scaling of the axes.

5.5 2.5 1.8

2.4
5 5.0 1.7
1.7

4.5
1.6
4 2
1.5 1.5
RMS misfit

RMS misfit

RMS misfit
3.5
1.7 1.4
3
1.3 1.3
2.5 2.4 1.5 1.5
1.2 1.2
2 1.3
1.7 1.2 1.11
1.5 1.1
1.5 1.3 1.11
1.2 1.06
1 1 1
0 50 100 150 200 250 0 50 100 150 200 250 300 350 400 450 0 50 100 150 200 250 300 350 400 450 500 550
Roughness measure Roughness measure Roughness measure

Constable, Orange, and Key, 2015

There really isn’t an objective way to choose misfit level except through a
good understanding of the data errors.
Model parameterization:

geology Data

Data errors
and misfit choice
Geophysical
model
inversion
algorithm
Regularization

Priors/
Model
constraints
Parameterization
 real
world
data
space
model
space

Best fit
Even with your best efforts, the real world is unlikely to be captured by your model
parameterization, and the best fitting model almost certainly won’t be either.
Understanding this can be important.

(And the best fitting model won’t have the same misfit as the real world does,
because of noise.)
Where in model space you are is determined by your parameterization - this also
determines where in data space you can be.

In non-linear geophysical problems, even forward modeling can involve a

challenging computational effort.

model space

2D 1D
data
space
3,4D
Sven Treitel once asked the question: “Can our mathematics ever completely
describe nature?”.

The trite answer, of course, is “No”. However, it is more useful to understand the
nature of the limitations:

Are the physics sufficient (e.g. scalar properties versus anisotropy)?

Is the forward computational machinery accurate? (e.g. finite

difference calculations don’t handle bathymetry well)

Is the dimensionality of model space large enough? (1D, 2D, 3D, 4D)

Is the discretization fine enough and the model size big enough?

One can rarely afford to blindly ensure these are all achieved, so intelligence and
understanding must be applied, perhaps by trial and error.
Model space parameterization:

d̂ = f (x, m) Some forward functional f

m = (m1 , m2 , ....., mN ) Model parameters

Most geophysical properties cannot go negative, but your inversion scheme

might well generate negative values in m. The easiest way to handle this is by
parameterizing as log(m), but there are other ways, such as NNLS.
Model space parameterization:

d̂ = f (x, m) Some forward functional f

m = (m1 , m2 , ....., mN ) Model parameters

In the real world, N (model size) is infinite (even in 1D). How we proceed from here
depends on whether N is small, moderately large, or infinite.

Small (sparse) parameterizations can be handled with parameterized inversions

(e.g. Marquardt) or stochastic inversions. The concept of least squares fitting
works because sparse models don’t have the freedom to mimic the pathological
true least squares solutions.
Model space parameterization:

d̂ = f (x, m) Some forward functional f

m = (m1 , m2 , ....., mN ) Model parameters

In the real world, N (model size) is infinite (even in 1D). How we proceed from here
depends on whether N is small, moderately large, or infinite.

Small (sparse) parameterizations can be handled with parameterized inversions

(e.g. Marquardt) or stochastic inversions. The concept of least squares fitting
works because sparse models don’t have the freedom to mimic the pathological
true least squares solutions.

Infinite N requires a real inverse theory mathematician. I am not one of them.

Model space parameterization:

d̂ = f (x, m) Some forward functional f

m = (m1 , m2 , ....., mN ) Model parameters

In the real world, N (model size) is infinite (even in 1D). How we proceed from here
depends on whether N is small, moderately large, or infinite.

Small (sparse) parameterizations can be handled with parameterized inversions

(e.g. Marquardt) or stochastic inversions. The concept of least squares fitting
works because sparse models don’t have the freedom to mimic the pathological
true least squares solutions.

Infinite N requires a real inverse theory mathematician. I am not one of them.

Most of the time geophysicists are working with moderately large N. Also, many
geophysical problems are non-linear, so we will concentrate on that approach.
To invert non-linear forward problems we often linearize around a starting model:

d̂ = f (m1 ) = f (m0 + m) f (m0 ) + J m

using a matrix of derivatives

⇥f (xi , m0 )
Jij =
⇥mj
and a model perturbation

m = m1 m0 = ( m1 , m2 , ...., mN )
2
Now our expression for is

2
⇥ ||Wd Wf (m0 ) + WJ m|| 2
For a least squares solution we solve in the usual way by differentiating and setting
to zero to get a linear system:

⇥= m
where
= (WJ) W(d T
f (m0 ))
= (WJ)T WJ .

So, given a starting model m0 we can find an update m :

m= 1
⇥
and iterate until we converge. (This is Gauss-Newton.)
Global versus local minima:
For nonlinear problems, there are no guarantees that Gauss-Newton will
converge.
There are no guarantees that if it does converge the solution is a global one.
The solution might well depend on the starting model.

global minimum
local minimum
(maybe)
2
Global versus local minima:
For nonlinear problems, there are no guarantees that Gauss-Newton will
converge.
There are no guarantees that if it does converge the solution is a global one.
The solution might well depend on the starting model.

global minimum
local minimum
(maybe)
2

Gauss-Newton only works for small N (it isn’t even defined for N > M). If N gets
too large then the solutions become unstable, oscillatory, and generally useless
(they are probably trying to converge to D+ type solutions).
Almost all inversion today incorporates some type of regularization, which
minimizes some aspect of the model as well as fit to data:

U = ||Wd Wf (m)||2 + µ||Rm||2

where Rm is some measure of the model and µ is a trade-off parameter or
Lagrange multiplier. In 1D a typical R might be:

⇥
1 1 0 0 0 ... 0 m1 -1

⇧ 0 1 1 0 0 ... 0 ⌃ m2 +1 -1
⇧ ⌃ m3 +1 -1
⇧ 0 0 1 1 0 ... 0 ⌃
⇧ ⌃ m4 +1
R1 = ⇧
-1
.. .. ⌃
⇧
⇧ . . ⌃
⌃
m5 +1 -1
m6 +1 -1
⇤ ⌅
m7 +1 -1
1 1 m8 +1

which extracts a measure of slope. This stabilizes the inversion, creates a single
solution, allows N > M, and manufactures models with useful properties.

This is easily extended to 2D and 3D modeling.

The trade-off between roughness and misfit:

U = ||Wd Wf (m)||2 + µ||Rm||2

When µ is small, model roughness is ignored and we try to fit the data. When µ is
large, we smooth the model at the expense of data fit.
The trade-off between roughness and misfit:

U = ||Wd Wf (m)||2 + µ||Rm||2

When µ is small, model roughness is ignored and we try to fit the data. When µ is
large, we smooth the model at the expense of data fit.

One approach is to choose µ and minimize U by least squares, but picking µa

priori is simply choosing how rough your model is.
The trade-off between roughness and misfit:

U = ||Wd Wf (m)||2 + µ||Rm||2

When µ is small, model roughness is ignored and we try to fit the data. When µ is
large, we smooth the model at the expense of data fit.

One approach is to choose µ and minimize U by least squares, but picking µa

priori is simply choosing how rough your model is.

We ought to have a decent idea of how well our data can be fit. This forms the
basis of the “Occam” approach, where a target data misfit 2 is chosen:

U = ||Wd Wf (m)||2 2
⇤ + µ||Rm||2
The trade-off between roughness and misfit:

U = ||Wd Wf (m)||2 + µ||Rm||2

When µ is small, model roughness is ignored and we try to fit the data. When µ is
large, we smooth the model at the expense of data fit.

One approach is to choose µ and minimize U by least squares, but picking µa

priori is simply choosing how rough your model is.

We ought to have a decent idea of how well our data can be fit. This forms the
basis of the “Occam” approach, where a target data misfit 2 is chosen:

U = ||Wd Wf (m)||2 2
⇤ + µ||Rm||2
For linearized, iterative inversion we use
⇥ ⇥
U = ||Rm1 ||2 + µ 1
||Wd W f (m0 ) + J(m1 m0 ) ||2 ⇥2⇥
After differentiation and setting to zero we get an expression for a new model:

⇥ 1
m1 = µR R + (WJ) WJ
T T
(WJ)T W(d f (m0 ) + Jm0 ) .
If the Occam algorithm does not get hung up in a local minimum, it will converge to
the smoothest model for a given misfit.
Sum-square misfit

p2
=
p1
en
wh starting model
0
s s=
e
hn
two parameter

ug
ro
smoothest
acceptable model

least squares model surface of acceptable

misfit
one parameter

It is important to run the inversion to convergence, and not stop as soon as the
target misfit is achieved.
I will step through a joint 2D Occam inversion of marine CSEM (3 frequencies, no
phase) and marine MT (Gemini salt prospect, Gulf of Mexico):
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, Gemini_joint_inv_2pt4_a.[Link]
Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 5.2768 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 3.7362 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.4814 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.3978 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.4027 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.3993 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.3982 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.3974 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.398 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.402 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
 10 70
All
CSEM 60
8 MT
50

Roughness
RMS misfit
6 40

30
4
20

2 10

0
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Iteration number Iteration number

Rho y, RMS: 2.4008 Gemini_joint_inv_2pt4_a.[Link]

Folder: 40
2
0

1.5
2

4 1
Depth (km)

log10(ohm−m)
6
0.5

8
0

10
−0.5

12

−1
8 10 12 14 16 18 20 22 24 26
3178157N 3179571N 3180985N 3182399N 3183814N 3185228N 3186642N 3188056N 3189471N 3190885N
337657E 339071E 340485E 341899E 343314E 344728E 346142E 347556E 348971E 350385E
Why do we start from a half-space? Because J depends on m.

MT: misaligned starting resistor - no harm done

Courtesy David Myer.

misaligned starting conductor - forever trapped by J

Courtesy David Myer.

 Even with well-estimated errors, choice of misfit can still be somewhat subjective.
2
0

1.5
2

4 1 2
0

log10(ohm−m)
Depth (km)

6
0.5 1.5
2
2
0
8
4 0 1
1.5
2

log10(ohm−m)
10

Depth (km)
6
−0.5 0.5

12 RMS 2.40 8
4 1

log10(ohm−m)
−1 0

Depth (km)
8 10 12 14 16 18 20 22 24 26
6
0.5
10
−0.5
8
12 RMS 1.70 0

10 −1
8 10 12 14 16 18 20 22 24 26
−0.5

0
2
12 RMS 1.30
−1
8 10 12 14 16 18 20 22 24 26
1.5
2

2
0
4 1
2
log10(ohm−m)

0 1.5
Depth (km)

2
6
0.5

1.5
4 2 1
8
0

log10(ohm−m)
Depth (km)

6 4 1
10 0.5

−0.5

Depth (km)
12 RMS 1.20 8 6
0
0.5

−1
8 10 12 14 16 18 20 22 24 26
10 8
0
−0.5

12
RMS 1.11 10
−0.5
8 10 12 14 16 18 20 22 24 26
12
RMS 1.06 −1

Constable, Orange, and Key, 2015 8 10 12 14 16 18 20 22 24 26

−1
What about anisotropy? It is quite common for physical properties of sediments
to be different in the vertical and horizontal directions. For example, horizontal
resistivity ⇢h is often smaller than vertical resistivity ⇢z .

The problem is how to weight the penalty between the two models.

⇢z ⇢h
m1 -1 m1 -1
m2 +1 -1 m2 +1 -1
m3 +1 -1 m3 +1 -1
m4 +1 -1 ? m4 +1 -1
m5 +1 -1 m5 +1 -1
m6 +1 -1 m6 +1 -1
m7 +1 -1 m7 +1 -1
m8 +1 m8 +1

0 1
⇥ 1 1 0 0 0 ... 0
1 1 0 0 0 ... 0 B 0 1 1 0 0 ... 0 C
⇧ 0 1 1 0 0 ... 0 ⌃ B C
⇧ ⌃ B 0 0 1 1 0 ... 0 C
⇧ 0 0 1 1 0 0 ⌃ B C
⇧ ... ⌃ R2 = B .. .. C
R1 = ⇧ .. .. ⌃ B
B . . C
C
⇧
⇧ . . ⌃
⌃ @ A
⇤ ⌅
1 1
1 1
 Joint Anisotropic Inversions versus penalty between rho-y and rho-z (all fitting to
RMS 1.2):
2 2 2
0 0 0

1.5 1.5 1.5

2 2 2

4 1 4 1 4 1

log10(ohm−m)
log10(ohm−m)

log10(ohm−m)
Depth (km)
Depth (km)
Depth (km)

6 6 6
0.5 0.5 0.5

8 8 8
0 0 0

10 10 10

12
⇢h −0.5

12
⇢h −0.5

12
⇢h −0.5

−1 −1 −1
8 10 12 14 16 18 20 22 24 26 8 10 12 14 16 18 20 22 24 26 8 10 12 14 16 18 20 22 24 26

2
2 2 0
0 0

1.5
1.5 1.5 2
2 2

4 1 4 1
4 1

log10(ohm−m)
log10(ohm−m)
log10(ohm−m)

Depth (km)
Depth (km)
Depth (km)

6 6
6 0.5 0.5
0.5

8 8
8
0 0
0

10 10
10

12
⇢z −0.5

12
⇢z −0.5

12
⇢z −0.5

−1 −1 −1
8 10 12 14 16 18 20 22 24 26 8 10 12 14 16 18 20 22 24 26 8 10 12 14 16 18 20 22 24 26

Low weight (0.1), Medium weight (1), High weight (10),

models are models look sensible. models are identical.
independent.
Constable, Orange, and Key, 2015
There are many ways to choose how to regularize the problem, and this
matters too.

geology Data

Data errors
and misfit choice
Geophysical
inversion model
algorithm
Regularization

Priors/
Model
constraints
Parameterization
 For a given misfit, the model depends on R

Minimum norm model First deriv smoothing -1 1 Second deriv smoothing 1 -2 1

0 0 0

-500 -500 -500

-1000 -1000 -1000

Depth, m

Depth, m

Depth, m
-1500 -1500 -1500

-2000 -2000 -2000

(1 0 ) (-1 1 ) (1 -2 1 )
-2500 -2500 -2500
-1 0 1 2 -1 0 1 2 -1 0 1 2
10 10 10 10 10 10 10 10 10 10 10 10
Resistivity, Ohm.m Resistivity, Ohm.m Resistivity, Ohm.m

Depth weighted smoothing Thickness weighted smoothing

0 0

-500 -500

-1000 -1000
Depth, m

Depth, m

-1500 -1500

-2000 -2000

-2500
(-1 1) x thick/depth -2500
(-1 1) x thickness
-1 0 1 2 -1 0 1 2
10 10 10 10 10 10 10 10
Resistivity, Ohm.m Resistivity, Ohm.m
 You can have fun with cuts (removing a row of R):
Cuts in smoothing: Top and bottom Cuts in smoothing: Too deep Cuts in smoothing: Too shallow
0 0 0

-500 -500 -500

-1000 -1000 -1000

Depth, m

Depth, m

Depth, m
-1500 -1500 -1500

-2000 -2000 -2000

-2500 -2500 -2500

-1 0 1 2 -1 0 1 2 -1 0 1 2
10 10 10 10 10 10 10 10 10 10 10 10

Resistivity, Ohm.m Resistivity, Ohm.m Resistivity, Ohm.m

Cuts in smoothing: Top only Cuts in smoothing: Bottom only Marquardt 3-layer model RMS 1.018
0 0 0

-500 -500 -500

-1000 -1000 -1000

Depth, m

Depth, m

Depth, m
-1500 -1500 -1500

-2000 -2000 -2000

-2500 -2500 -2500

-1 0 1 2 -1 0 1 2 -1 0 1 2
10 10 10 10 10 10 10 10 10 10 10 10

Resistivity, Ohm.m Resistivity, Ohm.m Resistivity, Ohm.m

 You can have fun with cuts (removing a row of R):
Cuts in smoothing: Top and bottom Cuts in smoothing: Too deep Cuts in smoothing: Too shallow
0 0 0

-500 -500 -500

-1000 -1000 -1000

Depth, m

Depth, m

Depth, m
-1500 -1500 -1500

-2000 -2000 -2000

-2500 -2500 -2500

-1 0 1 2 -1 0 1 2 -1 0 1 2
10 10 10 10 10 10 10 10 10 10 10 10

Resistivity, Ohm.m Resistivity, Ohm.m Resistivity, Ohm.m

Cuts in smoothing: Top only Cuts in smoothing: Bottom only Marquardt 3-layer model RMS 1.018
0 0 0

-500 -500 -500

-1000 -1000 -1000

Depth, m

Depth, m

Depth, m
-1500 -1500 -1500

-2000 -2000 -2000

-2500 -2500 -2500

-1 0 1 2 -1 0 1 2 -1 0 1 2
10 10 10 10 10 10 10 10 10 10 10 10

Resistivity, Ohm.m Resistivity, Ohm.m Resistivity, Ohm.m

Sparsely parameterized model does well! (Why?)

In 2D:

Synthetic Model

Inverted with isotropic smoothing

Inverted with 1st-difference

default regularization generates

Courtesy Brent Wheelock.
a conductive artifact
Always remember that regularization has an input into the model solution. Both
of these MT models fit the data equally well.

Default R:
⇥
1 1 0 0 0 ... 0
⇧ 0 1 1 0 0 ... 0 ⌃
⇧ ⌃
⇧ 0 0 1 1 0 ... 0 ⌃
⇧ ⌃
R1 = ⇧ .. .. ⌃
⇧
⇧ . . ⌃
⌃
⇤ ⌅
1 1

A nonlinear and adaptive

R, providing little
penalty for big contrasts

Courtesy David Myer.

Can we place errors or uncertainties on regularized models?

No!

For sparse parameterizations, data errors are often projected onto model
parameters through the Jacobian. This is a dangerous practice because

• it depends on the parameterization

• the Jacobian depends on the solution
Can we place errors or uncertainties on regularized models?

No!

For sparse parameterizations, data errors are often projected onto model
parameters through the Jacobian. This is a dangerous practice because

• it depends on the parameterization

• the Jacobian depends on the solution
For regularized models, however, we compound this by generating huge amounts
of covariance between the model parameters as part of the smoothness constraint.
It is much more useful to think of regularized models as extremal solutions, and
vary the regularization to ask the questions you may have.

Stochastic methods provide a useful way to assess model uncertainty, but they are
still restricted to simple, mostly 1D, models.
More on errors and regularization:
Regularized Fits to f(x) = 1.0x + 1.0
When creating synthetic data
for inversion tests, always
perturb the data with noise,
don’t just add error bars.

This is because regularized

inversion will use its misfit
budget to make the model
smaller.

You actually get better

models by adding noise.

Constable, 1991
It can even matter how you scale the data.

geology Data

Data errors
and misfit choice
Geophysical
inversion model
algorithm
Regularization

Priors/
Model
constraints
Parameterization
In EM, both MT apparent resistivities and CSEM amplitudes can vary by many orders
of magnitude. This suggests that one should use error floors that are percentages.

One might also parameterize the data as logs. For small ✏:

d0 ± 0.434✏ = log10 (d ± ✏d)

log data linear fractional

data data error

where 0.434 = 1/ ln(10)

In EM, both MT apparent resistivities and CSEM amplitudes can vary by many orders
of magnitude. This suggests that one should use error floors that are percentages.

One might also parameterize the data as logs. For small ✏:

d0 ± 0.434✏ = log10 (d ± ✏d)

log data linear fractional

data data error

where 0.434 = 1/ ln(10)

It ought not to matter how you parameterize the data (so long as the errors are
properly scaled and the appropriate chain rule is applied to the Jacobian):
⇥ 1
m1 = µR R + (WJ) WJ
T T
(WJ)T W(d f (m0 ) + Jm0 ) .

make data use the chain

log rule to convert
take the log @d/@m to @ log(d)/@m
of the forward
model
In EM, both MT apparent resistivities and CSEM amplitudes can vary by many orders
of magnitude. This suggests that one should use error floors that are percentages.

One might also parameterize the data as logs. For small ✏:

d0 ± 0.434✏ = log10 (d ± ✏d)

log data linear fractional

data data error

where 0.434 = 1/ ln(10)

It ought not to matter how you parameterize the data (so long as the errors are
properly scaled and the appropriate chain rule is applied to the Jacobian):
⇥ 1
m1 = µR R + (WJ) WJ
T T
(WJ)T W(d f (m0 ) + Jm0 ) .

But ...
... it does. Consider misfits in marine CSEM data:

in the linear domain, this misfit is 10 times

bigger than the other one

.9 x 10-10 (linear) or 1 (log)

.9 x 10-11 (linear) or 1 (log)

... it does. Consider misfits in marine CSEM data:
in the linear domain, this misfit is only 10%
bigger than the other one

.99 x 10-11 (linear) or 2 (log)

.9 x 10-11 (linear) or 1 (log)
Here we consider MT data
3
over a half-space, varying 10
True Model
only half-space resistivity R. 90 Ωm

RMS Misfit
2
10
Recall that MT impedance
(Z ) is:
Ex = ZHy 10
1

1 2
⇤= |Z| linear rho, phase
2⇥f µ 0 impedance, phase
10 log rho, phase

For small R, misfit flattens for 10

−1
10
0
10
1
10 10
2 3
10
4

linear Halfspace Resistivity

−6

−4
Air
−2

0 1
♦
Depth (km)

2 2
20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 ♦
4 ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦

8 Seafloor
10

−180 −160 −140 −120 −100 −80 −60 −40 −20 0

Horizontal Position (km)

(modified from Wheelock et al., 2015)

Here we consider MT data
3
over a half-space, varying 10
True Model
only half-space resistivity R. 90 Ωm

RMS Misfit
2
10
Recall that MT impedance
(Z ) is:
Ex = ZHy 10
1

1 2
⇤= |Z| linear rho, phase
2⇥f µ 0 impedance, phase
10 log rho, phase

For small R, misfit flattens for 10

−1
10
0
10
1
10 10
2 3
10
4

linear Halfspace Resistivity

This is because −6

−4

(d f (m)) ! const. −2
Air

0 1

as f (m) ! 0 but
♦
Depth (km)

2 2
20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 ♦
4 ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦

( log d log f (m)) 6

8 Seafloor

does not. 10

−180 −160 −140 −120 −100 −80 −60 −40 −20 0

Horizontal Position (km)

(modified from Wheelock et al., 2015)

 3 This effect can be even worse
10
True Model
90 Ωm
for marine MT data affected by
RMS Misfit

bathymetry.
2
10
Local minima develop, and
misfit flatlines at low R.
1
10

linear rho, phase

0 impedance
10 log rho, phase
−1 0 1 2 3 4
10 10 10 10 10 10
Seafloor Resistivity

(modified from Wheelock et al., 2015)

Here is a 2D example. Linear apparent resistivity and phase converged, but
log(resistivity) converges to the same model in half the iterations. MT
impedance, Z, did not converge (d.n.c) at all.

5 Ωm

Truth 100 Ωm
Z, 5%, 10 Ωm, d.n.c.

log10(ρa), ΦZ , 5%, 10 Ωm, 11 iter. ρa, ΦZ , 5%, 10 Ωm, 21 iter.

(modified from Wheelock et al., 2015)

log10(ρa), ΦZ , 10%, 10 Ωm, 14 iter. ρa, ΦZ , 10%, 1000 Ωm, 24 iter.

So I hope I have convinced you that models from geophysical inversion depend
on much more than the data alone:

geology Data

Data errors
and misfit choice
Geophysical
inversion model
algorithm
Regularization

Priors/
Model
constraints
Parameterization
So I hope I have convinced you that models from geophysical inversion depend
on much more than the data alone:

geology Data

Data errors
and misfit choice
Geophysical
inversion model
algorithm
Regularization

Priors/
Model
constraints
Parameterization

Does this mean that geophysical inversion is useless?

Not at all! There is plenty of evidence that geophysical inversion works very
well. You just have to know what you are Rho-z
doing, how your code and algorithm
work, and pay attention to the factors besides the data that determine the
result.

And, you will probably need to run more than one inversion… many more.

0.2 1 10 Ωm Folder: 11AI_Anis_Joint_1pt0_2

2

1500
1 1.5

2
2000
1

log10(ohm−m)
Depth (km)

2500 4 0.5

0
3000
6

7 −0.5
3500
meters bls

−1
8 10 12 14 16 18 20
0N 0N 0N 0N 0N 0N 0N
Connecting the world of applied geophysics