Report

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 ! iT ( x ) T(x) Travel-time T(x) Amplitude A(x) Receiver 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 Grad(T) orthogonal to wavefront 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 dS1T1 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(2A.T A 2T ) 0 2 2 2A( 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 dT d t . c( x)T . therefore ds dT c 2 c 1 (T ) ( 2 ) ds 2 2 c ds cT .(T ) dT 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 Many analytical solutions (gradient of velocity; gradient of slowness square) 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 Radial Structure 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 Many analytical solutions (gradient of velocity; gradient of slowness square) Time integration of ray equations Runge-Kutta second-order integration 1D sampling of 2D/3D medium : FAST source Predictor-Corrector integration stiffness receiver Initial conditions receiver ? 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 may backtrace using the gradient of T from any point of the medium towards the source (should be applied from each receiver) 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 q0 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 source and receiver ! 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 di x . where the shooting angle is i . The estimation of the derivative is through paraxial computation d d qx di d d qx di d d qx di z 0 z 0 d d qx dpix d d qx dpix z 0 dpix d d qx di 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 additional information, 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 … The symmetrical radial EARTH 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) receiver t (s, r ) t ( s, r ) receiver u( x, y, z )dl receiver u0 ( x, y, z )dl source source receiver0 receiver0 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 ) receiver0 d u ( x, y, z )dl source0 d t ( s, r ) receiver0 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)Cd1 (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)Cd1 (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 G01 V p Lp1U tp Up and Vp have now inverses ! Solution, model & data resolution The solution is mest G01d G01 (G0m) Vp Lp1U 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 SPREAD(R)= SPREAD(N)= RI 2 N I 2 Spreading functions : 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 Cd1 (d G0 m) (m m p )t Cm1 (m m p ) 1 G0g G0t Cd1G0 Cm1 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 add additional relations between model parameters (new lines) 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 ) Gradient method E (mk tE ( 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 Conjugate gradient Gauss-Newton is Quasi-Newton for L2 norm quadratic 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 Cd1/ 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 and of Vadim Monteiller 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)