Probabilistic uncertainty estimation for 2D/3D refraction seismic traveltime tomography

The development of an objective method for the estimation of quantitative uncertainty in seismic tomography is one of the most challenging and urgent tasks of applied seismology. For one-dimensional models, probabilistic methods are successful in estimating uncertainty, but their application to higher-dimensional models suffers from extensive computational demands. However, the lack of quantitative uncertainties for higher-dimensional seismic inverse models often leads to speculations, erroneous projections and unfruitful debates that devalue seismic tomography in general. Two instruments are proposed to facilitate the application of Monte Carlo Markov Chains to higher-dimensional refraction-and-reflection-seismic models. First, irregular grids will be used that allow for adapting the node distribution to the resolving power of the data. For most realistic problems, this leads to a significant reduction of the number of inverse parameters. Secondly, linear constraints - most importantly the resolution matrix - will be used to counteract the decrease of the acceptance frequency of randomly perturbed models resulting from an increased number of inverse parameters. In damped least-squares inversion the resolution matrix indicates the linear dependency between the inverse nodes, as given by the ray geometry. By scaling any random perturbation at a certain node with the corresponding row of this matrix one can derive additional perturbations that compensate the bias of the modeled data that would result from a single perturbation. The additional computational expense for the resolution matrix is small compared to the forward solution, but the expected benefit for the acceptance rate is large without reducing the perturbations. Estimated CPU times show that this is a realistic approach, at least for small 2D models, perhaps even for small 3D models. The successful implementation of the proposed uncertainty estimation will significantly increase the value of the refraction seismic method.

In order to constrain the uniqueness and accuracy of seismic velocity models, we developed a global inverse algorithm for refraction seismic problems. Refraction seismology generally aims at reconstructing the elastic properties of the subsurface to constrain its structure and composition. The method is successfully applied to geologic, geotechnical, or hydrogeological problems, and for the prospection for resources. In contrast to local optimization methods, our approach provides not a single solution but a probability density function that describes the entire solution space and better captures the nature of the non-linear seismic problem. Probability functions help to avoid misinterpretations and they facilitate risk assessment. During the development, efficiency was a central task of this project in order to make the code applicable to large models. One approach is the use of staggered grids. The individual grids are coarse, and computations are relatively fast and can be made in parallel without increasing wall-clock time. The final result at the required resolution is then obtained from the superposition of the individual grids. The second approach is the compensation of the random perturbations that are applied during the course of the inversion. Those compensations are computed from measures of the correlation of individual model parameters. They increase the acceptance ratio and the step width, such that the relevant model space is sampled more efficiently.Our results are available as statistical test series, and in form of a computer code that can be applied to refraction seismic 2D and 3D surveys. This coded was tested with a synthetic model, and eventually with a real dataset from the Salzach valley near Zell am See.