Solving or resolving inadequate and noisy tomographic systems

Abstract

Tomography in seismology often leads to underdetermined and inconsistent systems of linear equations. When solving, care must be taken to keep the propagation of data errors under control. In this paper I test the applicability of three types of damped least-squares algorithms to the kind of sparse matrices encountered in seismic tomography: (1) singular value decomposition with Lanczos iteration, (2) conjugate gradient iteration with the LSQR algorithm, and (3) the Dines-Lytle method. Lanczos iteration may be applied to large sparse systems of low rank to calculate solutions by singular value decomposition but becomes impractical with problems of larger size. The Paige-Saunders algorithm (LSQR), which incorporates Lanczos' iteration into a conjugate gradient method, provides a least-squares solution to the system with acceptable filtering properties. In a synthetic tomographic experiment, it proved to be converging up to an order of magnitude faster than the Dines-Lytle algorithm, a stationary iterative process. For large tomographic systems, where restrictions in the available computer time pose limitations on the number of iterations, this indicates that conjugate directions methods are to be preferred to the more commonly applied Gauss-Seidel type of algorithms. To avoid unwarranted conclusions when the problem is severely underdetermined or undermined with large data errors, resolution analysis must be applied to the data set. An algorithm for the determination of the resolution in sparse systems is given.