Calculating resolution and covariance matrices for seismic tomography with the LSQR method

In seismic tomography, the LSQR algorithm is commonly used for solving the inverse problem. LSQR belongs to the family of conjugate gradient methods, so the generalized inverse is not solved explicitly. Consequently, neither the covariance nor the resolution matrix are provided by LSQR, which limits one’s ability to obtain estimates of uncertainty and errors in the computed models. In this paper we present a method, demonstrated by synthetic examples, for calculating the resolution and covariance matrices via the general inverse of LSQR. The extra computational effort is limited and only a few lines of computer code in the original LSQR routine are needed to produce the required output. The resolution matrices produced demonstrate that great care must be taken to ensure that a sufficient number of iterations are used when applying LSQR inversion. The relationship of the LSQR-based resolution estimates to those produced using other methods is briefly discussed.