Introduction to Seismology - Traveltime Tomography

Seismic tomography We will study travel time tomography (body waves), which can be assimilated to imaging. But one has to be aware that imaging can be interpreted as reflection seismic, velocity reconstruction by the industry (SEG). One can also look at surface wave tomography but these two types of tomography are not part of the same theory. There is no easy way to combine them; the assumptions made are different in each case. The steps: - We have observations, represented by dots on the figure below (picking of the different waves: P, S, PP…) Image removed due to copyright considerations. Please see: Stein, S., and M. Wysession. Figure 3.5-3 in An Introduction to Seismology, Earthquakes, and Earth Structure. 1st ed. Malden, MA: Blackwell, September 2002, p. 163. ISBN: 0865420785. - People generate linear regression to give the best fit of these observations - We obtain models (ex: IASP91, PREM…) [Wiechert and Herglotz (GER) found a 1D classical model from travel time curves.] - We use an inversion technique to get the P & S wave velocity as a function of the depth (see the figure below). 1 Image removed due to copyright considerations. Please see: Stein, S., and M. Wysession. Figure 3.5-1 in An Introduction to Seismology, Earthquakes, and Earth Structure. 1st ed. Malden, MA: Blackwell, September 2002, p. 162. ISBN: 0865420785. - From these models, we can produce travel time curves again. Go back to the observations again and look at the deviation from the model → the scattering around it. δt = Tobs - Tref T . .. δt. .. . . . ∆ We look at the difference between the data and the model and NOT at the travel time itself! OBJECTIVE: find a model (3D) that minimizes δt. Tobs= 3D structure + error in source location + noise (instrument, measurement) Hypocenter: to,x,y,z (or θ,ϕ,r) ⇒ δt = δt3D+δtmislocation+δtnoise But can we isolate δt3D? What people do is combine the model (δt3D) and the source location (δtmislocation). The noise is assumed to be white, with a Gaussian distribution. 2 R S ⇒ Travel time dl (to,x,y,z) If we know the 1D model, we can apply Snell’s law to estimate the geometry: 1 Tobs = ∫ dl v(r) 3D ray path But we know that: 1 ∇T = k̂ c If c changes, the ray path changes. We end up with a non-linear problem. The ray path depends on the answer of the problem. What do we need to do? ⇒ linearize the inversion To do so, we use the Fermat’s principle: T ← Example of a minimum ray path Changes in ray path If we change a little bit the ray path around the optimum, we’ll end up with a small change in travel time. ∆v Reference ray We have two kinds of “deviations” from the reference ray: 1. First contribution: effect of changes in velocity ∆v 2. Second contribution: effect of change in the ray path Fermat’s principle says that the 2nd effect is of 2nd order: we can ignore it! 3 Linearization of the travel time δt = Tobs - Tref = T – To= 1 1 ∫ v(r) dl − ∫ v (r) dl ' 3D ray (unknown) o Reference ray (known), model = ∫ s.dl − ∫ so .dl ' Because the deviation of the ray is ignored, we only look at the effect of the heterogeneities. It means that in linearizing the travel time, we obtain: nd ∫ s.dl ≅ ∫ so .dl ' + 2 order effect 3D ray (unknown) Reference ray Reference ray (known) deviation of the ray path Heterogeneities δ t = ∫ (s − so ).dl ' = ∫ ∆s.dl ' Reference ray (known) In linearizing the problem, we get rid of the unknown ray. We can do our calculation in a reference earth model. [Rayleigh ‘s principle is similar] The travel time tomography is an iterative process: - create 1D model - ray tracing and get new rays in the model - update ray geometry - get the reference ray related to the 3D δt = Tobs - Tref (3D) (The reference model does not have to be a 1D model.) Linearization of the hypocenter mislocation ∂T ∂T ∂T ∂T δ to + δx + δy + δz ∂t o ∂x ∂y ∂z ∂⊗ δ⊗ = δ⊕ ∂⊕ with to the origin time, (x,y,z) the location of the earthquake. δ t misloc = 4 ∆s δto ,δx, δy, δz What we want to know is: From the travel time tomography, we get the wavespeed model and the estimation of the mislocation of the source. • If we want to work with the velocity, not the slowness: ∆v − ∫ 2 dl ' for small ∆v = v - vo Fermat’s principle vo ∆s ≅ − • ∆v vo2 conversion Inversion problem: linear relationship between δt and the model. observations δt = T – To model perturbation (m-mo) ∫ di = Gi (r)µ (r)dr volume sensitivity kernel/ Frechet derivatives volume d⊗ = ∂⊗ δ⊕ ∂⊕ Partial derivatives / Frechet derivatives ƒ The first thing to do next is to discretize the problem: do a parametrization (not work with a continuous integral). M µ (r ) = ∑ γ k hk (r ) k =1 γk: weight hk: basis function For example, 1. Fourier theory, harmonic functions: hk (r) = cos(kθ ) Just add a weight and we obtain a signal. 5 2. Plane wave summation: φ = ∫∫∫ Φ(x, y, z, w)ei ( kx − wt )dk.dw weight ƒ basis function We inject µ into the equation: M di = ∑ k =1 {∫ G (r)h (r)dr}γ i k k kernel projected on the basis function h M di = ∑ Aik γ k for i=1,2,…, N data k =1 sensitivity matrix kernel projected on the basis functions projector d = Am m = (γ 1 , γ 2 ,..., γ M ) model (includes the values of the weight) data ƒ Green functions: solution of a point source Gi is a Green function for a point perturbation. We need to do a convolution with a point perturbation in order to get the observations. ƒ M: number of model parameters N: number of observations (data) in general, N≠M. As a consequence, the matrix A is not square in general. A= N rows M columns m = A-1d → we cannot do that! 6 We have to multiply by AT: AT d = (AT A)m - ATA is a square matrix (M x M) → m̂ = (AT A)−1 AT d Generalized Least Squares inversion || m̂ − ( AT A) −1 AT d ||= ε → Goal is to minimize ε. ƒ Usually, we work with 10 million x 10 million matrix. So, we need to use numerical techniques to solve the problem: LSQR, SIRT… d = Am ƒ δ t = ∫ ∆s.dl → parametrize ∆s M ∆s = ∑ γ k hk k =1 One way is to take hk as a series of cells/blocks. hk will have a value for r inside the cell k and will be zero for r outside the cell k. Another way: - hk as splines, wavelets - hk: spherical harmonics (in global seismology) hklm = fk (r)Yl m (θ , ϕ ) surface spherical harmonic radial function model perturbation: µ (r)= µ (r,θ ,ϕ )∑∑∑ k α lm cos(mϕ ) + k βlm sin(mϕ ) fk (r)P ml (cosθ ) ( k l ) m for gravity (1/r)l ∑ γ × h = µ(r) ∑ weight × basisgfunction = mod el k k 7 Legendre polynomial - δ ti = ∫ ∆s.dl = path length Mcell ∑ ∆s (dl) k ik k =1 i: event-station pair in a block k dl.s = dl/v T - d=Am with A the sensitivity matrix containing the path length l1 l2……….lm Am = ∆s1 ∆s2 … ∆sm model = δt1 δt2 … δtm data vector each ray gives a row We have an average wavespeed along the ray. In order to construct a model vector, we need to get data from different rays crossing each other. A is a sparse matrix. If we look at one ray: M~300,000 ~ 20 layers the ray samples ~ 100 cells A will have only ~ 100 elements non-zero. The good thing about sparse matrix is that ATA is ~ diagonal. The problem is that there are many singularities, which makes the inversion unstable (in that case, we need to add a damping factor or regularize the problem). 8 One possibility is to not use cells of the same size. Consequently, it reduces the number of cells; the inverse matrix is less singular. Nevertheless, the computation time increases. 9