Introduction to Seismology - Model Uniqueness

TRAVELTIME TOMOGRAPHY (CONT) MODEL UNIQUENESS The forward model d i Data vector = A m ik k Sensitivity matrix Model vector ⇔ d = Am (1) states the linearized relationship between data and model via a mapping or sensitivity matrix, where i = 1, 2, … , N (data) and k = 1, 2, … , M (model parameters). In the case of GRT traveltime tomography, A can be thought of as a sparse sampling matrix whose nonzero elements are the lengths of segments of the raypath within each discretized voxel of the Earth. ⎡ l11 L l1M ⎤ ← ray 1 ⎢ M A≡ M O M ⎥ ⎢ ⎥ ⎢⎣l N 1 L l NM ⎦⎥ ← ray N • Normally, N >> M – the number of data (typically ~106 – 108) is far greater than the number of model parameters (typically ~105 – 106) – making (1) nominally overdetermined. • Consequently, A is normally not square. We can invert for m by first left-multiplying (1) by AT A T Am = A T d and then, since ATA is square, we complete the inversion to yield the classic least squares solution, as follows: −1 ) m = (A T A ) A T d . (2) The penalty function that is satisfied by (2) is ε ≡ m − (A T A ) A T d −1 2 (3) ) where m (plugged in for m in (3)) minimizes ε in the least squares sense. Aside: penalty metrics other than the squared norm in (3) can be used; the absolute norm, for example. As stated above, N >> M Æ the system is nominally overdetermined. But A is often rank-deficient – mathematically, rank(A) << M – and therefore the system behaves as though it is underdetermined. In other words, many of the N observations are linearly dependent and so do not provide complementary information on the subsurface. To understand this, consider the worst-case scenario where we have 106 observations that travel along the exact same raypath like the figure shown at right. Only cells traversed by the ray are actually sampled, meaning that unsampled cells are completely unconstrained. From an inverse point of view, the unconstrained parameter cells make the problem underdetermined – there are more parameters than 1 S R there are complementary pieces of information to constrain the system. When an inverse problem is underdetermined, all the following statements are true and/or synonymous: {1} the problem is ill-posed; {2} ATA is singular/non-invertible/rank-deficient; {3} one or more rows or columns of A and ATA are linearly dependent; {4} one or more parameters are unconstrained; {5} A and ATA both have nontrivial null spaces. All parts of a discretized model that are unconstrained – or, equivalently, that are not mapped by A from model space into data space – reside in the null space of A. This means that ε is independent of unconstrained parameters; an unconstrained parameter can take on any value without affecting ε . So the least squares solution is not unique if one or more parameters are unconstrained! In fact, there is an infinite set of solutions that equally minimize ε . Namely, any prospective solution that contains the true solution (which we don’t know) linearly combined with any of the vectors in null (A) will equally minimize ε . This is termed nonuniqueness. REGULARIZATION AND DAMPING Model regularization essentially adds additional constraints to the inverse problem to address the solution nonuniqueness that arises as a consequence of ill-posedness. For example, the expression a+b=2 is ill-posed if we are asked to find values for both a and b. There are an infinite number of combinations of a and b that equally satisfy this expression, so the inverse problem is nonunique. But if we add an extra condition or constraint that specifies λ1(a – b) = 0, where λ1 is a user-specified weight (a.k.a. Lagrange multiplier), we now have a well-posed system of two equations in two unknowns (provided λ1 ≠ 0). The unique solution to this system is a = b = 1. This second expression is a form of model regularization. It assumes that a functional relationship exists between the model parameters (in this case, it basically assumes that a = b) that can be used to add information and thereby constrain the original problem. This idea can be applied to (3) by augmenting it with an additional regularization condition: ε = m − (A T A ) A T d + λ1 Lm −1 2 2 (4) where, for example, L can defined as the gradient or first-difference operator ⎤ ⎡1 − 1 0 L ⎥ ⎢0 1 − 1 0 L ⎥ L≡⎢ 1 − 1 0 L⎥ ⎢0 0 ⎥ ⎢ O O O O⎦ ⎣M that basically reduces the difference between adjacent model parameters in m. This type of constraint requires that neighboring model parameters be of similar magnitude. For this reason, the second term on the RHS of (4) is often referred to as a smoothness constraint. The strength 2 of the weight λ 1 specifies how important parameter smoothness should be in relation to the main objective (the first term on the RHS), which is to minimize the residual in the least squares sense. Regularization can be implemented by augmenting the system of equations. In block matrix form, the regularized forward problem is written ⎡A ⎤ ⎡d⎤ ⎢ L ⎥ m = ⎢L⎥ , ⎢ ⎥ ⎢ ⎥ ⎢⎣ 0 ⎥⎦ ⎢⎣λ1L⎥⎦ (5) ) and m can be found in the usual way by (2). Another constraint customarily employed in inversion is norm damping. For unconstrained parameters, it is often desirable to at least constrain them to be physically reasonable. In fact, if we don’t have any data that constrains particular parameters, it’s just as well to make them zero, or at least something small. Going back to the simple two-parameter problem at the beginning of this section, we might additionally constrain the problem by a second constraint: λ2 a = 0 λ2 b = 0 again, where λ2 is a weight (Lagrange multiplier) specifying the importance of this objective. Applying this idea to (4), we create a new penalty function, ε = m − (A T A ) A T d + λ1 Lm + λ2 m . 2 −1 Squared error 2 2 Smoothness constraint Norm Constraint (6) Damping applies a norm constraint to the problem. That is to say, it constrains the parameters to be small; large values in m would make its norm large and would therefore make ε big. In the limit as λ2 Æ ∞, all the parameters would be zero. The augmented forward problem is expressed as ⎡A ⎤ ⎡d⎤ ⎢L⎥ ⎢L⎥ ⎥ ⎢ ⎢ ⎥ ⎢λ1L⎥ m = ⎢ 0 ⎥ , ⎥ ⎢ ⎢ ⎥ ⎢L⎥ ⎢L⎥ ⎢⎣ λ2 I ⎥⎦ ⎢⎣ 0 ⎥⎦ (7) ) and, as with (5), m can be found by employing (2). A RESOLUTION Equation (1) establishes a forward mapping from model to data space (cartooned by the top arrow in the figure at right) and Equation (2) establishes the inverse mapping from data to model space 3 Data Space Model Space (A A ) T −1 AT (cartooned by the bottom arrow). We would especially like to know how well we resolve model parameters via the inverse mapping. Combining (1) and (2), we get −1 ) m ≅ (A T A ) (A T A )m , (8) and we define the resolution matrix as R ≡ (A T A ) (A T A ) −1 (9) so that ) m ≅ Rm . Clearly, if R = I then all the parameters are equally and perfectly resolved. In practice, R will often be diagonally dominant, but the main diagonals will not all be one and the off-diagonals will often be different from zero. By analyzing the diagonals of R, we gather information on how well resolved the parameters are (resolution increases as the diagonal values approach unity). The off-diagonal terms tell us about ‘smearing’ between parameters. However, in seismic tomography, A frequently contains billions or even trillions of elements, so formally computing R is impractical (impossible?). There are a number of ways to circumnavigate this difficulty. The two methods described here are the checkerboard resolution test and the hypothesis test. The checkerboard resolution test basically utilizes a synthetic Earth model wherein seismic velocities alternate between slow and fast, in a checkerboard pattern: msyn ≡ checkerboard. msyn is mapped to synthetic data by dsyn = Amsyn, and then dsyn is inverted to create the inversion model −1 ) m syn = (A T A ) A T d syn . We can then consider the resolution metric ) m syn − m syn 2 ) or we can simply visually inspect m syn to see how well it tallies with msyn. The following observations can be made about checkerboard resolution tests: • The greatest model resolution, not surprisingly, occurs over continental landmasses • Many regions over open oceans apparently cannot be resolved at all. • Deep mantle model parameters are generally better resolved than ones nearer the surface because there is better data coverage provided by deep diving waves. • Model regularizations used in the inversions tend to provide poorer amplitude resolution than spatial resolution. • These test types assume A is known. But A is an approximation of physical reality. 4 The hypothesis testing method is less arbitrary than the checkerboard method. Rather than test the resolution of an unrealistic Earth model, which is only marginally informative at best, this method tests a null hypothesis. The null hypothesis can be any Earth model that is similar to the current inversion model, but that differs in some way that would be important from an interpretive viewpoint. Let the current and null hypothesis models be denoted, respectively, m0 and mnull. The null hypothesis Earth model is mapped to data space via A and then inverted back to model space via (ATA)-1AT – so we’re considering (ATA)-1(ATA)mnull. If mnull is similar to (ATA)-1(ATA)mnull at those places where mnull was purposely made to differ from m0, then there is sufficient resolution in these regions and the null hypothesis is rejected. If (ATA)-1(ATA)mnull looks like m0 then the null hypothesis is accepted – the data do not constrain the parameters enough to distinguish between m0 and mnull. So, when the null hypothesis is rejected, we have resolution in regions of interest; but if it is accepted, then both the null hypothesis and the original inversion model can equally explain the data, which presents a problem as far as interpretation is concerned. Hypothesis testing is easily visualized using the block diagram shown in Figure 1. Tomographic model with anomaly (m0) Null hypothesis model (mnull) with dashed outline of m0 for comparison Hypothesis test Figure 1 Block diagram cartooning hypothesis testing for resolution. The top left panel depicts a tomographic image with one anomaly. We would like to know, for example, whether the data provides enough resolution to unequivocally interpret this anomaly as a contiguous body or perhaps as two separate bodies (as seen in the top right null hypothesis panel). To do T -1 T this, apply (A A) A A to mnull and check to see whether the two-anomaly hypothesis is preserved. If it is, the hypothesis is rejected; there is enough resolution to distinguish between the two models. If not, the hypothesis is accepted; we do not have enough resolution to reasonably distinguish between the null hypothesis and the original image. (A A ) (A A )m T −1 T null or Reject – Resolution is sufficient to distinguish between models. 5 Accept – Resolution is insufficient to distinguish between models. FRESNEL ZONES AND MULTIPATH/SCATTERING We have ∫ Gδµ dV d= Volume but we only need to worry about those portions of Gδµ that are nonzero. Geometric Ray Theory employs a high-frequency approximation that allows us to assume that Gδµ is only nonzero along the ray path. So the preceding integral can be reduced to δt = ∫ ∆s dl . ray This simplification breaks down when we consider finite frequencies. For finite frequencies, we employ the analytical tools of Fresnel zones and multipathing/scattering. These concepts provide more realistic 3D kernels applicable to Earth models. Fresnel Zones L’ S R L Rays traveling from S to R will follow not only the GRT raypath, but also any number of excursive paths. If the lengths of these excursive paths, L’, do not differ from L (the GRT raypath length) by more λ/4, they will constructively interfere. In other words: L′ − L < λ 4 ⎫ ⎬ constructive interference T ′ − T < T 4⎭ L′ − L = s 2 + w2 + (L − s )2 + w2 − L ↓ w S R 1/ 2 ⎛ λL ⎞ w≅⎜ ⎟ ⎝ 8 ⎠ S L-S w denotes the radius of a Fresnel Zone. It specifies a cross-sectional area through which constructively interfering rays travel from S to R. So the waveform observed at R is no longer dependent only on the GRT raypath, but the volume described by the Fresnel Zone. The arrival at R has now sampled a volume, as opposed to a line. L ~ 1000km Æ w ~ 30km L ~10,000km Æ w ~ 100km The excursive raypaths arrive later than the GRT raypath, creating constructive interference and the arrival coda. 6 The Fresnel Zone is shaped like a banana and is frequently referred to as the banana donut: In cross-section, the banana donut is represented by a series of concentric contours that correspond to the spatial sensitivity of the banana donut at each point. Interestingly, for finite frequencies, the sensitivity is identically zero at the center of the cross-section, which provides the counterintuitive insight that the GRT raypath is a path of zero sensitivity! Banana Donut Sensitivity Kernel Zero sensitivity The banana donut kernel (BDK) replaces the elements in A, as described earlier in these notes. Because the BDK is a volume and not a raypath, A is no longer nearly as sparse as it is in the case of GRT. Moreover, because sensitivity is now distributed in a volume, the BDK A is less ill-posed than the GRT A. Hence, banana donut kernels will tend to increase resolution, as compared to GRT kernels. Multipathing The preceding ideas also naturally give rise to multipathing, which acknowledges that more than one raypath can satisfy Fermat’s Principle. For example, a “scattered” ray might travel through a faster region and therefore arrive at the same time as the GRT ray. Consequently, the integral d= fast ∫ Gδµ dV Volume can be simplified to the sum of two integrals δt = ∫ ∆s dl + ∫ ∆s dl , ray1 and this is the central idea of multipathing. 7 ray 2