Model parameterization

Optimizing model parameterization in seismic traveltime tomography

Introduction

In seismic traveltime tomography, a set of linearized equations is solved for the unknown slowness perturbations. The matrix of this set of equations (the tomographic matrix) is usually ill-conditioned, because the model space contains more details than can be resolved using the available data. The condition of the tomographic matrix can be characterized by its singular value spectrum: the larger the normalized singular values and the greater the rank of the matrix, the smaller the null space and the better the condition of the inversion problem. The structure of the tomographic matrix depends on the source-receiver configuration and the model parameterization. Since the source-receiver geometry is often fixed (e.g. in earthquake tomography), conditioning can only be influenced through the model parameterization, i.e. the structure of the model space.

Most often in traveltime inversion, the Earth is modelled by regular grid of cells, where the local slowness is assumed to be constant. The element G_ij of the tomographic matrix G is then the length of the ith ray in the jth cell. Because of the regular parameterization not taking into account the raypath geometry, cells not hit by any rays can often be found and it is also often the case that some rows of G are linearly dependent. This means that G may be rank-deficient and has a null space. As a consequence, an infinite number of solutions can be obtained by adding any linear combination of vectors in the null space to the solution of the tomographic equations.

In order to minimize the null space of the tomographic system of equations, I have proposed a method for finding an optimal, irregular triangular cell parameterization that best suits the given raypath geometry (Wéber 2001).

Optimization method

The greater the rank and the larger the singular values of G, the more stable the inversion problem and the more independent pieces of information may be gained from the data. The singular values of G are the positive square roots of the eigenvalues of G^TG. The first eigenvalue e₁ can be calculated quickly using the power method, while the sum of the eigenvalues is equal to the trace of G^TG. Unfortunately, determining the rank p would require the whole singular value decomposition (SVD) of G, which is a computationally expensive task, given the size of the matrix. However, p can be estimated by the number of model cells in which ray coverage is greater than a predefined threshold. This threshold should be chosen according to the average ray coverage. Although high ray coverage is not a sufficient criterion to eliminate the null space, it is a necessary criterion for a well-determined inverse problem.

Now, let's define the cost function

E(G) = ne₁ ⁄ trace(G^TG) + q

where n is the number of the model cells and q = n - p denotes the estimated number of singular values in the null space. This cost function has the minimum of 1, and any reduction in E(G) results in a better conditioned tomographic matrix with smaller null space.

The best parameterization that minimizes the above cost function is expected to consist of irregularly sized cells, so I use Delaunay triangular cells to describe the model. Delaunay triangulation is a procedure to generate a unique set of triangles from arbitrarily distributed nodes in two dimensions. Since the parameterization depends on the nodal configuration, I seek the best model by searching for the optimal nodal distribution.

Numerical experiments

Now I illustrate the proposed procedure by finding the optimal triangular cell parameterization for a cross-borehole tomographic inversion problem. The ray coverage threshold used to estimate q in the cost function is chosen to be equal to one tenth of the average ray coverage (sum of ray lengths divided by the number of cells).

The cross-borehole raypath geometry with rectangular grid parameterization and with the optimized triangular cell parameterization is illustrated, respectively, in the a and b parts of the figure below. Although the gridded model consists of 40 cells, while the triangulated one contains only 39 triangles, the presented results can be readily compared.

The distribution of the triangular cells in Fig. b is correlated with the raypath distribution: the greater the ray and angular coverage, the smaller the triangles. The shape and orientation of the triangles also depend on the ray distribution. In the upper part of Fig. b, for example, where sub-horizontal raypaths dominate, the triangles have horizontally elongated shape according to the well-known fact that sub-horizontal rays produce little horizontal resolution. In other words, the optimized model with its triangles of different sizes and shapes reflects the local resolution of the inversion.

The figure above compares the normalized eigenvalue spectra of G^TG for the two cell types. The normalized eigenvalues for the optimized triangulated model are considerably larger than those for the gridded one. The grid parameterization has resulted in four zero eigenvalues with cost value of 9.28, while for the triangulated model, the tomographic matrix has full rank and the cost value is 4.47.

Further numerical experiments and their detailed discussion can be found in Wéber (2001).

Acknowledgements

The reported investigation was financially supported by the Hungarian Scientific Research Fund (No. F019277).