A computationally feasible approximate resolution matrix for seismic inverse problems

Geophys. J. Int. (1996) 126, 345-359 A computationally feasible approximate resolution matrix for seismic inverse problems Susan E. Minkoff* The Rice lnversion Project, Department ofComputationa1 and Applied Mathematics, Rice Unioersity, PO Box 1892, Houston T X 7725 1-1892, U S A Accepted 1996 February 29. Received 1996 February 13; in original form 1995 October 7 Key words: inversion, seismic resolution, source time functions, viscoelasticity. INTRODUCTION In a classic paper, Backus & Gilbert (1968) gave a description of the general linear inverse problem and defined model resolution. They explained that solutions of seismic inverse problems are generally not pointwise unique but may still have unique average behaviour. Resolution quantifies how close the model estimates are to the true model as a function of spatial location (such as depth). It indicates the shortest length-scale which the given data can resolve. [References include papers by Franklin (1970), Jackson (1979), Kennett & Nolet ( 1978), and the book by Menke (1989).] Wiggins (1972) defines this matrix in terms of the singular value decomposition (SVD) and analyses resolution for a very small inversion problem (approximately 10 unknowns) using observed data from surface waves and free oscillations. He comments that, ‘computation of parameter and information resolution is such a simple extension of any inversion procedure based on perturbation parameters that such inversion studies are incomplete without considering resolution.’ Unfortunately, forming the resolution *Now at: Texas Institute for Computational and Applied Mathematics, Taylor Hall 2.316, University of Texas, Austin, TX 78712, USA. 0 1996 RAS matrix via the SVD has not become common practice because it is computationally prohibitive for large problems. There are a number of papers in which SVD resolution is computed for small inverse problems. For example, using waveform inversion and a viscoelastic simulator, Martinez & McMechan (1991) examine model resolution for a simple geometry experiment of a target layer located between two homogeneous half-spaces. For field data experiments they do not compute resolution. Bishop et al. (1985) analyse the resolution of components of a simplified tomography model analytically, and substantiate their conclusions with a small numerical experiment performed on real data. There have also been attempts to compute the full SVD for realistic (large) inverse problems [see, for example, Assous & Collin0 (199O)l. Ory & Pratt (1995) examine the effect of different regularizing operators on resolution. They consider a synthetic experiment of determining three 1-D anisotropic velocity parameters from traveltime data, and state (p. 420): One practical problem that arises is the cost of computing the model resolution matrix, either because of too large a number of data, or too large a number of parameters ... or both. It took about 3 hours CPU time ... to compute the matrix in our 345 SUMMARY Seismic inversion produces model estimates which are at most unique in an average sense. The model resolution matrix quantifies the spatial extent over which the estimate averages the true model. Although the resolution matrix has traditionally been defined in terms of the singular value decomposition of the discretized forward problem, this computation is prohibitive for inverse problems of realistic size. Inversion requires one to solve a large normal matrix system which is best tackled by an iterative technique such as the conjugate gradient method. The close connection between the conjugate gradient and Lanczos algorithms allows us to construct an extremely inexpensive approximation to the model resolution matrix. Synthetic experiments indicate the data dependence of this particular approximation. The approximation is very good in the vicinity of large events in the data. Two large linear viscoelastic inversion experiments on p-t marine data from the Gulf of Mexico provide estimates of the elastic parameter reflectivities corresponding to two different seismic sources. Traditionally, one evaluates the accuracy of the two reflectivity estimates by comparing them with measured well logs. The approximate model resolution matrices agree with the well-log ranking of the two models and provide us with a way to compare different model estimates when, for example, such well-log measurements are not available. 346 S. E . Minkof one-dimensional anisotropic parameter problem ( 576 parameters and 1420 data). For larger problems, this computation becomes prohibitive. We wish to quantify our ability to resolve the individual model parameters from the forward model for the seismogram G . It is necessary to devise a generalized inverse or estimator G - g which will act on the data and return an estimate of the model parameters. The simplest way to define the model resolution matrix is to assume that there is a set of model parameters mtrueE W”which satisfies the equation Gmtrue=dabs, for the The forward modeller G E W”x n and observed data dobs E 9”. resolution matrix (Menke 1989) measures the averaging of the model mtruewhich enters into the model estimate mest in the inversion process: mest - ~ - g d o b s= ~ - g ~ ~ t = r Rmtrue u e (1) E W”x’2is the model resolution matrix. If R = I then - lyltrue, and the model parameters are perfectly resolved in the inversion. In general, R f I, and then the model estimates mest are weighted averages of the true parameters where R mest = aTmtrue . Note that this averaging vector a corresponds to a column of the resolution matrix R . In the experiments I will describe, the model parameters depend on one spatial coordinate only, namely depth ( z ) ,and the averaging vectors correspond to averages over depths. The singular value decomposition allows us to decompose the forward operator into the matrix product G = U Z V T (Golub & Van Loan 1989, page 71). The columns of U are eigenvectors of GGT (the left singular vectors of G ) , and the columns of V are the eigenvectors of G T G (the right singular vectors of G ) . Finally, the p singular values on the diagonal of C = diag(al, ... , a), are the square roots of the non-zero eigenvalues of GGT and G T G (the singular values). Of importance to us in this paper is the fact that V s columns are eigenvectors which span the model space M . Only the p non-zero singular values contribute to the model resolution matrix, so we work instead with the truncated SVD defined by G = U p C pV,’ where U p and V,’ consist of the first p columns of U and V respectively, and C, is a p x p diagonal matrix. We note that, although VE V, = I , in general, since the p right singular vectors do not span the whole space, V, V,’ # I . A version of the estimator G-g can now be expressed in terms of the truncated singular value decomposition, G-8 = V,C-’ , UT,, and the model resolution matrix [from expression ( I ) ] is simply mfst In my implementation, I do not e (...truncated)