### Presentation pptx - Web site of Jean Virieux

```Asymptotic wave solutions
Jean Virieux
Professeur UJF
Translucid Earth
S(t)
u( x, t )  A( x)S (t  T ( x))
u( x,  )  A( x)S ( )e
Source
Same shape !
iT ( x )
T(x)
Travel-time T(x)
Amplitude A(x)
Wavefront : T(x)=T0
Too diffracting medium : Wavefront preserved
wavefront coherence lost !
Eikonal equation
Two simple interpretation of
wavefront evolution
Orthogonal trajectories are rays
T=cte
Velocity c(x)
L
T+T
c( x) 
L T
1
1


  xT ( x) 
T L c( x)
c( x)
1
( xT ( x))  2
c ( x)
2
Direction ? : abs or square
Transport Equation
Tracing neighboring rays defines a ray tube :
variation of amplitude depends on section variation
d  A1 dS1T1  A2 dS2 T2
 
 
2
2
 A1 T1.n1dS1  A2 T2 .n2 dS2
 
 
2
2
0   A2 T2 .n2 dS2  A1 T1.n1dS1
 
 '
 
2
2
2
0   A2 T2 .n2 dS2  A1 T1.n 1dS1   Ac Tc .nc dSc


2
2
0   div( A T )d  div( A T )  0
 
A(2A.T  A 2T )  0
2
2


2A( x).T ( x)  A( x)2T ( x)  0
Ray tracing equations

T
Ray

x (s)

dx
Evolution of
ds
T=cte

 dx

 t 1
s curvilinear abscisse t 
ds
 


dx
dx
// T   c( x)T ( x)
ds
ds
evolution of x


d T but
Evolution of  T is given by
 ds
  




dT
d
 t .  c( x)T . therefore
 ds
dT c   2 c  1
 (T )  ( 2 )
ds
2
2 c
ds
 cT .(T )

dT  1
 ( )
ds
c 
evolution of  T

 1
d 1 dx
(
)  (
)
Curvature equation
known as the ray equation
ds c( x) ds
c( x)

 

We often note the slowness vector p  T (x) and the position q  x (s )
System of ray equations
dq
 cp
ds
dp
1

ds
c
dq
 c2 p
dT
dp
1
 c
dT
c
dq
p
d
dp 1 1
 
d c c
d  cds
with
cdT  ds
c 2 dT  d
which ODE to select for numerical solving ? Either T or sampling
Velocity variation v(z)
Ray equations are
dqy
dqx
dq
 px ;
 p y ; z  pz
d
d
d
dpy
dpx
dp
du( z )
 0;
 0; z  u ( z )
d
d
d
dz
The horizontal component of the
slowness vector is constant: the
dq x p x
px
trajectory is inside a plan which is


dq z p z  u 2 ( z )  p 2
called the plan of propagation. We
x
may define the frame (xoz) as this
Where px is a constante
plane.
dqx
dqz
 px ;
 pz
z
d
d
px
q
(
z
,
p
)

q
(
z
,
p
)

x
1
x1
x
0
x0
z u 2 ( z)  p 2 dz
dpx
dpz
du( z )
 0;
 u( z)
x
d
d
dz
For a ray towards the depth
1
0
Velocity variation v(z)
At a given maximum depth zp, the
slowness vector is horizontal following
the equation
z
z
p
qx ( z1 , px1 )  qx 0  
u ( z )  px
2
z0
zp
T ( z1 , px1 )  T0  
z0
p
px
u ( z )  px
z1
zp
u 2 ( z)
2
2
dz  
2
dz  
z1
zp
px
u ( z)  px
2
u 2 ( z)
u ( z)  px
2
2
2
dz
px  p 2  u 2 ( z p )
2
dz
If we consider a source at the free surface as well as the receiver, we get
a
zp
X ( p)  2 
0
zp
T ( p)  2 
0
p
u ( z)  p
2
2
u ( z)  p
rp
a
u 2 ( z)
2
dz
  2
2
dz
T  2
rp
p
dr
r 2u ( r ) 2  p 2 r
r 2u ( r ) 2
dr
r 2u ( r ) 2  p 2 r
with p  ru sin i
with p = usini
In Cartesian frame
In Spherical frame
Velocity structure in the Earth

Main discontinuous boundaries



Crustal discontinuity (Mohorovicic (moho – 30 km)
Mantle discontinuity (Gutenberg – 2900 km) D’’
Core discontinuity (Lehman – 5100 km)
Shadow zone at the Gutenberg discontinuity: thickness of few kms
Minor discontinuities



Interface at 100-200 km
Interface at 670-700 km
Interface at 15 km (discontinuity of Conrad)
These discontinuities are related to
lithospheric, mesospheric and sismogenic
structures.
They are not expected to be deployed over
the entire globe
System of ray equations


dx
 cp
ds

dp
1

ds
c

dx 2 
c p
dT
 
dp
1
 c
dT
c

dx 
p
d

dp 1  1
 
d c c
dy
 A( y)
d
Non-linear ODE !
d  cds
with
cdT  ds
c dT  d
2
which ODE to select for numerical solving ? Either T or sampling
square)
Time integration of ray equations
Runge-Kutta second-order integration
1D sampling of 2D/3D medium : FAST
source
Predictor-Corrector integration
stiffness
Initial conditions
?
Boundary conditions VERY DIFFICULT
?
source
Shooting
dp ?
Bending
dx ?
Continuing dc ?
EASY
Save p conditions
if possible !
AND FROM TIME TO TIME IT FAILS !
But we need 2 points ray tracing because we have a source and a receiver
to be connected !
Eikonal Solvers
Fast marching method (FMM)
Layered
medium
T
z=cte x
T
z z
T 2 T 2 1
) ( )  2
x
z
c
T
1
T 2
  2 ( )
z
c
x
(tracking interface/wavefront motion :
curve and surface evolution)
(
z + dz
EIKONAL SOLVER
From Sethian, 1999
Let us assume T is known at a level z=cte
Compute
T
x
along z=cte by a finite difference
T
z
and extend T estimation at depth
approximation
Deduce
z+dz
Assuming 1D plane wave propagation,
we have been able to estimate T at a
depth z+dz
Fast marching method (FMM)
2D case
(
T 2 T 2
1
) ( ) 
x
z
c ( x, z ) 2
(tracking interface/wavefront motion :
curve and surface evolution)
T
1
T 2


(
)
2
z
c ( x, z )
x
Strategy available in 2D and 3D BUT
only for the minimum time in each node
in the spatial domain (x,y,z).
Possible extension in the phase
domain (for multiphases) ?
Sharp interfaces are always difficult to
handle in this discrete formulation
From Sethian, 1999
Other techniques as fast sweeping methods
3D case
Wavefront tracking
Back-raytracing strategy
Smooth medium : simple case
Time over the grid
Once traveltime T is computed
over the grid for one source, we
of T from any point of the medium
towards the source (should be
Ray
The surface {MIN TIME} is convex
as time increases from the source :
one solution !
A VERY GOOD TOOL
FOR First Arrival Time
Tomography (FATT)
Ray tracing by wavefronts
Rays and wavefronts in an homogeneous medium.
(Lambaré et al., 1996)
Rays and wavefronts in a constant gradient of the velocity.
(Lambaré et al., 1996)
(Lambaré et al., 1996)
Rays and wave fronts in a complex medium.
(Lambaré et al., 1996)
(Lambaré et al., 1996)
(Lambaré et al., 1996)
Z in m
0
400
800
1200
1600
2000
(Lambaré et al., 1996)
Ray tracing by wavefronts
Evolution over time :
folding of the wavefront is allowed
Dynamic sampling :
undersampling of ray fans
oversampling of ray fans
Keep sampling of the medium by rays « uniform »
Heavy task 2D & 3D !
Example of wavefront
evolution in Marmousi model
Hamilton Formulation

dq 
p
d

dp 1  1
  2
d 2 c
1 2 1
H ( q, p )  ( p  2 )
2
c
sympletic structure (FUN!)
mechanics : ray tracing is a balistic problem
Information around the ray
Ray
 
dq
  pH
d


dp
  x H
d
 
q0  dq
 
p0  dp

dq  

dp  
Meaning of the neighbooring zone – Fresnel
zone for example but also anything you wish
Paraxial Ray theory


0


 

d (q0  dq )

H
d
q


d
pq
  p 0 dp H (q0  dq , p0  dp)


d
d d p    qq H 0



ddq
  p 0  p 0 H (q0 , p0 )dp   p 0  q0 H (q0 , p0 )dq
d
d
(dy)  dy
d

does not depend on
dy
 pp H 0  d q 

 qp H 0  d p 
: LINEAR PROBLEM (SIMPLE) !
Estimation of ray tube :
amplitude
Estimation of taking-off
angles : shooting strategy
…
Practical issues
 d qx   0
d q   0
d  z 

2
d d px  0.5 u, xx

 
2
0.5
u
d
p
, zx
 z 
0
0
0.5 u,2xz
0.5 u,2zz
1 0   d qx 
0 1   d qz 


0 0  d px 


0 0  d pz 
Four elementary paraxial trajectories
dy1t(0)=(1,0,0,0)
dy2t(0)=(0,1,0,0)
dy3t(0)=(0,0,1,0)
dy4t(0)=(0,0,0,1)
d y  d qx , d qz , d px , d pz 
t
NOT A RAY !
Practical issues
dH 0
Paraxial rays require other conservative quantities : the perturbation of
the Hamiltonian should be zero (or, in other words, the eikonal
perturbation is zero)
2
2
x
x
z
z
,x
x
,z
z
1
1
p d p  p d p  u ( x, z ) d q  u ( x , z ) d q  0
2
2
Point source conditions dqx(0)=dqz(0)=0 pxd px  pzd pz  0
d px (0)  a pz (0)
This is enough to verify this condition initially
d pz (0)  a px (0)
Solution
d y( )  a pz (0)d y3( )  a px (0)d y4( )
From paraxial trajectories, one can combine them for paraxial rays as long
as the perturbation of the Hamiltonian is zero.
For a point source, a small slowness perturbation (a=10-4) gives a good
approximate ray depending on the velocity variation in the medium: it is a kind of
derivative ….
Practical issues
Plane wave conditions dpx(0)=dpz(0)=0
u ( x, z), x d qx  u ( x, z), z d qz  0
2
d qx (0)  a u 2 ( x, z ), x (0)
d qz (0)  a u 2 ( x, z ), z (0)
2
This is enough to verify this condition initially
Solution d y( )  au2 ( x, z), x (0)d y1( )  au 2 ( x, z), z (0)d y2( )
Two independent paraxial rays exist in 2D (while we have four elementary
paraxial trajectories)
Point paraxial ray
Plane paraxial ray
Seismic attributes
2 PT ray tracing non-linear problem solved, any
R attribute could be computed along this line :
-Time (for tomography)
A ray
S
-Amplitude (through paraxial ODE integration
fast)
-Polarisation, anisotropy and so on
Moreover, we may bend the ray for a more accurate ray tracing
less dependent of the grid step (FMM)
Log scale in time
One ray
Keep values of p at
Travel time evolution with the
grid step : blue for FMM and
black for recomputed time
Grid step
Time error over the grid (0)
Errors through FMM times
Errors through rays deduced after FMM times
NOT THE SAME COLOR SCALE
(factor 100)
Coarser grid for computation
ODE resolution


Runge-Kutta of second order
Write a computer program for an
analytical law for the velocity: take
a gradient with a component along
x and a component along z
Home work : redo the same thing with a Runge-Kutta of fourth order (look
after its definition)
Consider a gradient of the square of slowness
Runge-Kutta integration
Second-order integration
Non-linear ray tracing
df
 A( f )
d
 0 0
f  f 
.A ( f )
2
f 1  f 0   . A1/ 2 ( f 1/ 2 )
1/ 2
0
Second-order euler integration for paraxial ray tracing
df
 A. f
d
f  f   . A. f
1
0
0
 d qx   0
d q   0
d  z 

2

d d px  0.5u, xx

 
2
0.5
u
, zx
d pz  
0
0
2
, xz
2
, zz
0.5u
0.5u
1 0   d qx 

0 1   d qz 


0 0  d px 


0 0  d pz 
Four elementary paraxial trajectories
dy1t(0)=(1,0,0,0)
dy2t(0)=(0,1,0,0)
d y  d qx , d qz , d px , d pz 
t
dy3t(0)=(0,0,1,0)
Remember !
dy4t(0)=(0,0,0,1)
3 is for point solution along x
4 is for point solution along z
Two points ray tracing: the paraxial
shooting method
Consider x the distance between ray point at the free surface and
sensor position
dd q
x 
di
x
. where the shooting angle is i .
The estimation of the derivative is through paraxial computation
d d qx
di
d d qx
di
d d qx
di
z 0
z 0
d d qx

dpix
d d qx

dpix
z 0
dpix d d qx

di
dpiz
d d qx
p 
dpiz
z
i
z 0
 d qx3 piz  d qx 4 pix
z 0
z 0
dpiz
d i
pix
z 0
We can compute an
estimated and,
therefore, a new shooting
angle
Amplitude estimation
Consider L the distance between the exit point of a ray and the
paraxial ray running point.
L 
d qx pz ( r )  d qz px ( r )
px ( r )  p z ( r )
2
2
L (d qx3 piz  d qx4 pix ) pz (r )  (d qz3 piz  d qz4 pix ) px (r )


px ( r ) 2  pz (r ) 2
Thanks to the point paraxial ray estimation dq3 and dq4, we may
estimate the geometrical spreading L/and, therefore, the
amplitude A(r)
L
A(r ) 

L
Stop and move to the
numerical exercise
How to reconstruct
the velocity structure ?

Forward problem (easy)
from a known velocity structure, it is possible to compute travel
times, emergent distance and amplitudes.

Inverse problem (difficult)
from travel times (or similarly emergent distances), it is possible
to deduce the velocity structure : this is the time tomography
even more difficult is the diffraction tomography related to the
waveform and/or the preserved/true amplitude.
Tomographic approach

Very general problem
medicine; oceanography, climatology

Difficult problem when unknown a priori medium (travel
time tomography)

Easier problem if a first medium could be constructed:
perturbation techniques can be used for improving the
reconstruction (delayed travel time tomography)
Travel Time tomography
We must « invert » the travel time or the emergent
distance for getting z(u): we select the distance.

zp
X ( p)  2 
0

p2
X ( p)
dz / du2
dz;
 
du2
2
2
2p
2
u 2 ( z)  p 2
u

p
u0
p
Abel problem (1826)
Determination of the shape of a hill from travel
times of a ball launched at the bottom of the hill
with various initial velocity and coming back at
the initial position.
ABEL PROBLEM
A point of mass ONE and initial
velocity v0 reaches a maximal height x
1 2
given by
gx  v0

x
d
2
We shall take as the zero value for the
potential energy: this gives us the
following equations and its integration:
ds
y
ds / d 
 ds 

2
g
(
x


)

;

t
(
x
)

d
 

2g(x   )
 dt 
0
2
x
We may transform it into the so-called Abel integral
x
where t(x) is known and f() is the shape
f ( )
t ( x) 
d of the hill to be found: this is an
x 
integrale equation.
0

The exact inverse solution
We multiply and we
integrate
We inverse the order
of integration
We change variable
of integration


0


0


0
We differentiate and
write it down the
final expression
d
d

x


t ( x)
dx
f ( )
dx  
d

x
x 
0  x 0
t ( x)
dx
dx   f ( )d 
x
0
   x x 
2
2

x


cos



sin

t ( x)
dx    f ( )d
x
0


0
t ( x)dx
   f ( )
x
1 d
f ( ) 
 d


0
t ( x)dx
 x
THE ABEL SOLUTION
x
By changing variable  in
a- and x in a-x, we get the
standard formulae
t ( x)  
a
f ( )
d
 x
1 d
f ( )  
 d
We must have
t(x) should be continuous,
t(0)=0
t(x) should have a finit derivative with a finite
number of discontinuities.
The most restrictive assumption is the continuity of
the function t(x).
a

t ( x)dx
x 
The solution HWB : HERGLOTZ-WIECHERT-BATEMAN
p2
X ( p)
dz / du2
 
du2
2p
u2  p2
u0 2
z (v )  
z (v )  
1

u2

u0
2
u
1
 
u0
z (v ) 
1

X (u )

0
X ( p) / 2 p
p2  u2
X ( p)
p u
2
2
dp2
dp
From the direct solution, we
can deduce the inverson
solution
After few manipulations, we can
move from the Cartesian expression
towards the Spherical expression
dX
cosh(pv)
R
1
ln(
)
r (v ) 
(r / v)

0
d
cosh(pv / r )
We find r(v) as a value of r/v
In Cartesian frame
In Spherical frame
Stratified medium



We may find the interface at a depth h
when considering all waves
We may reconstruct an infinity of
structures with only direct and refracted
waves.
We have a velocity jump when a velocity
decrease
Velocity structure with depth

Velocity profile built without any a priori
An difficulty arises when
the velocity decreases.
An initial model through the HWB method
An initial model can be built
The exact inverse formulae does not allow to introduce

F. Press in 1968 has preferred the exhaustive exploration of possible
profiles (5 millions !). The quality of the profile is appreciated using a misfit
function as the sum of the square of delayed times as well as total volumic
mass and inertial moments well constrained from celestial mechanics ...
Exploration through grid search, Monte Carlo search,
simulated annealing, genetic algorithm, tabou method,
hant search …

A simple case : small perturbation
Initiale structure of
velocity
 Search of small variation
of velocity or slowness
 Linear approach

Example:
Massif Central
GLOBAL Tomography

Velocity variation at a
depth of 200 km : good
correlation with
superficial structures.
Velocity variations at a
depth of 1325 km : good
correlation with the
Geoid.
Courtesy of W. Spakman

Delayed Travel-time tomography
station
 u( x, y, z)dl
t ( source, station) 
source
Finding the slowness u(x) from t(s,r) is a difficult
problem: only techniques for one variable !
Consider small perturbations du(x) of the slowness field u(x)

t (s, r ) 
t ( s, r ) 
u( x, y, z )dl 

u0 ( x, y, z )dl 
source
source

u0 ( x, y, z )dl 
source0


d u( x, y, z )dl
source
d u ( x, y, z )dl
source0
t ( s, r )  t0 ( s, r ) 

d u ( x, y, z )dl
source0
d t ( s, r ) 

source0
d u ( x, y, z )dl
This a LINEAR PROBLEM dt(s,r)=G(du)
DESCRIPTION OF THE VELOCITY PERTURBATION
The velocity perturbation field (or the slowness
field) du(x,y,z) can be described into a meshed
cube regularly spaced in the three directions.
For each node, we specify a value ui,j,k. The
interpolation will be performed with functions as
step funcitons. We have found shape functions
h,i,j,k=1 pour i,j,k, and zero for other indices.

du( x, y, z)  dui , j ,k hi , j ,k
cube
Discrete Model Space
du( x, y, z)  dui , j ,k hi , j ,k
Slowness perturbation description
cube
  du
Discretisation of the
medium fats the ray
dt ( s, r ) 
Sensitivity matrice is a
sparse matrice
dt ( s, r )  
d t  G0d u
Matrice of sensitivity or of
partial derivatives
h
i , j ,k i , j ,k
rayon 0 cube
dl   dui , j ,k
cube
h
i , j ,k
rayon 0
t
dui , j ,k
cube u
 t1
 dt1   u
 dt   1
 2  
 ...   

 
dt n 1  
 dt n   t n
 u
 1
t1   du 
1

u m  du 
 2 








t n  du m 1 
u m   du m 
dl
LINEAR INVERSE PROBLEM
Updating slowness perturbation values from time residuals
Formally one can write
with the forward problem
Existence,
Uniqueness,
Discretisation
Identifiability
Stability,
Small errors
propagates
d u  G0 1d t  m  G0 1d
d t  G0d u  d  G0 m
Robustness
Outliers effects
of the model
NON-UNIQUENESS & NON-STABILITY : ILL-POSED PROBLEM
REGULARISATION : ILL-POSED -> WELL-POSED
LEAST-SQUARES SOLUTIONS
•The linear system can be
recast into a least-square
system, which means a
system of normal
equations. The resolution
of this system gives the
solution.
AtDT=AtA DM
DM=(AtA)-1AtDT
•The system is both under-determined and overdetermined depending on the considered zone (and tne
number of rays going through.
LEAST SQUARES METHOD
E (m)  (d  G0 m)t (d  G0 m)
L2 norm
E (m)
 0 locates the minimum of E
m
t
t
normal equations
G0 G0 m  G0 d


1
mest  G0 G0 G0 d
t
t
if
G G 
G G  is a M by M matrice
1
t
0
0
1
t
0
G0 is a N by M matrice
0
exists
Least-squares estimation
Operator
G G 
t
0
1
0
the generalized inverse
G0t
g
0
G
Under-determination M > N
Over-determination
on data will derive a new model : this is called
N>M
Mixed-determination – seismic tomography
Maximum Likelihood method
One assume a gaussian distribution of data
Joint distribution could be written
 1

p(d )  exp (d  G0 m)Cd1 (d  G0 m)
 2

Where G0m is the data mean and Cd is the data covariance matrice: this method
is very similar to the least squares method
E(m)  (d  G0m)t (d  G0m)  E1 (m)  (d  G0m)Cd1 (d  G0m)
Even without knowing the matrice Cd, we may consider data weight Wd
through
1
E2 (m)  (d  G0m)Wd (d  G0m)
SVD analysis for stability and uniqueness
G0  ULV
SVD decomposition :
U : (N x N) orthogonal Ut=U-1
t
UtU=I and VtV=I (not the inverse !)
V : (M x M) orthogonal Vt=V-1
L : (N x M) diagonal matrice
V  [V p V0 ]
U  [U p U 0 ]
Null space for Li=0
L p
G0  U p U 0 
0


0
V p V0

0


t
Vp and V0 determine the uniqueness while Up and U0 determine the existence of
the solution
G0  U p L pV pt
G01  V p Lp1U tp
Up and Vp have now inverses !
Solution, model & data resolution
The solution is
 
mest  G01d  G01 (G0m)  Vp Lp1U tpU p L pVpt m  VV pt m  Rm
where

R  VpVpt

Model resolution matrice : if V0=0

then
R=VVt=I

dest  G0mest  U pU tp d  Nd


where N  U pU tp d
Data resolution matrice
importance matrice
Goodness of resolution
RI
2
N I
2
: if U0=0 then N=UUt=I
Good tools for quality estimation
PRIOR INFORMATION
Hard bounds
A  mi  B
Seismic velocity should be positive
Prior model
E3 (m)  (d  G0m)t (d  G0m)   (m  mp )t (m  mp )
 is the damping parameter controlling the importance of the model mp
Gaussian distribution
E4 (m)  (d  G0 m)t Cd1 (d  G0 m)  (m  m p )t Cm1 (m  m p )


1
G0g  G0t Cd1G0  Cm1 G0t
Model smoothness
E5 (m)  (d  G0m)t Wd (d  G0m)  (m  mp )t Wm (m  mp )
With Wd data weighting and Wm model weighting
Penalty approach
UNCERTAINTY ESTIMATION
m  G G  G d  G d
Least squares solution
t
0
est
1
0
t
0
g
0
covmest   G0g covd G0gt  G0g Cd G0gt
Cd   d2 I
covmest   
Uncorrelated data
2
d
G G 
t
0
1
0
1
1  E 
covmest    
2
2

m

 m mest
2
2
d
Model covariance :
uncertainty in the data
curvature of the error function
Sampling the error function around the estimated model often this has to be
done numerically
A posteriori model covariance matrice
True a posteriori distribution
If one can decompose this matrice
1
1
G0 Cd G0  Cm  USU
t
S diagonal matrice
t
eigenvalues
U orthogonal matrice eigenvectors
Error ellipsoidal could be estimated
WARNING : formal estimation related
to the gaussian distribution hypothesis
Tangent gaussian distribution
A priori & A posteriori information
What is the meaning of the « final » model we provide ?
acceptable
Flow chart
d obs
collected data
m
starting model
loop
d syn  g (m)
true ray tracing
d  d obs  d syn
data residual
G0  g m
sensitivity matrice
m  G d
model update
m  m  m
new model
1
0
Calculate
small model variation or small errors exit
2E
m 2
for formal uncertainty estimation
Sampling a posteriori distribution
Resolution estimation : spike test
Boot-Strapping
Jack-knifing
Natural Neighboring
Monte-Carlo
Sampling a posteriori distribution
Uncertainty estimation for P and S velocities using boot-strapping techniques
Steepest descent methods
E(mk 1 )  E(mk )
E (mk  tE ( mk ))  E ( mk )
d  E (mk )
Newton
 E0


dk  

 Ek   k d k 1 
d   2 E ( mk ) 1 E (mk )
Quasi-Newton
 2 E (mk )  Dk
Gauss-Newton is Quasi-Newton for L2 norm
approximation of E
Tomographic descent
1
2
Minimisation of this vector
d  C g (m)


m p    C m 
1/ 2
d
1/ 2
m
 Cd1/ 2Gk 
Ak  

1/ 2
 Cm

If one computes
then
 C

 C
1/ 2
d
1/ 2
m

1/ 2

C
( g (mk )  d )
t
t
d
Ak Ak dm  Ak  1/ 2

 Cm (m0  mk ) 

Gaussian error distribution of data and of a posteriori model
Easy implementation once Gk has been computed
Extension using Sech transformation (reducing outliers
effects while keeping L2 norm simplicity)
2
LSQR method
The LSQR method is a conjugate gradient method developped by Paige & Saunders
Good numerical behaviour for ill-conditioned matrices
When compared to an SVD exact
solution, LSQR gives main
components of the solution while SVD
requires the entire set of eigenvectors
Fast convergence and minimal norm solution
(zero components in the null space if any)
Corinth Gulf
An extension zone
where there is a deep
drilling project.
How this rift is opening?
What are the physical
mechanisms of
extension (fractures,
fluides, isostatic
equilibrium)
Work of Diana Latorre
Seismic experiment 1991 (and
one in 2001)
MEDIUM 1D : HWB AND
RANDOM SELECTION
Velocity structure image
Horizontal sections
Velocity structure image
P
S
Vertical sections
Vp/Vs ratio:
fluid existence ?
Recovered parameters
might have diferent
interpretation and the
ratio Vp/Vs has a strong
relation with the
presence of fluids or the
relation Vp*Vs may be
related to porosity
Other methods of exploration






Grid search
Monte-Carlo (ponctual or
continuous)
Genetic algorithm
Simulated annealing and co
Tabou method
Natural Neighboring method
Conclusion FATT








Selection of an enough fine grid
Selection of the a priori model information
Selection of an initial model
FMM and BRT for 2PT-RT
Time and derivatives estimation
LSQR inversion
Update the model
Uncertainty analysis (Lanzos or numerical)
THANK YOU !
Many figures have come from people I have worked
with: many thanks to them !
Kirchhoff approximation
1)
2)
3)
Representation theorem
Kirchhoff summation
Reciprocity
Born approximation
1)
2)
Single scattering approximation
Surface approximation
(Forgues et al., 1996)
(Forgues et al., 1996)
```