Irregular interface parametrization in 3-D wide-angle seismic traveltime tomography

Summary

1 Introduction

A basic goal in seismic tomography is to maximize the extraction of structural information from data. In traveltime tomography, the data consist of source–receiver traveltimes and the structural information is in the form of variations in wave speed and/or reflector depth and shape. For most realistic problems, ray path coverage will not be uniform, so the minimum scalelength of structure that the data are able to resolve will vary across the model. This suggests that the data should play an important role in deciding how the subsurface is represented. However, this has rarely been the case in tomographic studies.

The traditional way of representing structure in seismic tomography is to adopt a regular parametrization. In global and regional traveltime tomography, layers of constant velocity cells (or cells in which the velocity varies linearly) which have a uniform size in latitude, longitude and depth are often used (Clayton 1983; Inoue et al. 1990; van der Hilst et al. 1997; Gorbatov et al. 2000). At the local scale, regular cells or grids of nodes with specified interpolation functions are common. For example, a parametrization which has been widely used in local earthquake tomography (e.g. Thurber 1983; Eberhart-Phillips 1986; Haslinger et al. 1999) and teleseismic tomography (e.g. Zhao et al. 1994; Steck et al. 1998) is a rectangular grid of velocity nodes with a trilinear interpolation function. Cardinal spline functions have also been used in teleseismic tomography by Thomson & Gubbins (1982), and by Sambridge (1990) in a combined local earthquake and controlled source study. Cubic splines have been employed in 2-D wide-angle tomography by Farra & Madariaga (1988), Lutter & Nowack (1990) and McCaughey & Singh (1997) and splines under tension have been used by Neele et al. (1993), VanDecar et al. (1995) and Ritsema et al. (1998) in teleseismic tomography. An alternative to a block or grid approach to parametrization is to discretize velocity in the wavenumber domain rather than the spatial domain. An example is a spectral parametrization that uses some form of truncated Fourier series (e.g. Wang & Pratt 1997). The minimum wavelength of structure permitted in the model can be controlled by the number of harmonic terms in the series. Hildebrand et al. (1989), Hammer et al. (1994) and Wiggins et al. (1996) have used this type of spectral parametrization in wide-angle tomography. In global traveltime tomography, truncated spherical harmonic expansions are sometimes used to represent structure (e.g. Dziewonski et al. 1977; Morelli & Dziewonski 1987), although they are more commonly used in global surface waveform tomography (e.g. Dziewonski & Woodhouse 1987; Trampert & Woodhouse 1995).

The drawback of a regular parametrization, whether it be cells, interpolated nodes or a truncated Fourier series, is that the minimum wavelength of structure permitted is essentially invariant throughout the model. This type of formalism is therefore not consistent with the inhomogeneous distribution of information that characterizes most data sets. In a solution, this inconsistency manifests as regions of limited resolution corresponding to cells or nodes that are poorly constrained by the data, and other regions where the parametrization is not flexible enough to respond to information contained in the data. The presence of poorly constrained parameters leads to an underdetermined inverse problem. In such cases, uniform spatial smoothing and/or damping are usually introduced, but regularizing the problem in this way can lead to a further loss of information in the well-constrained parts of the model due to oversmoothing or overdamping.

A better approach would be to have parameters placed only where they are required by the data so as to produce a solution that is well constrained everywhere yet satisfies the data to an acceptable level. However, there have been very few implementations of such schemes. Reasons for this include: (1) the difficulty of constructing a tomographic scheme that can handle irregular parametrization; (2) choosing node/cell positions so that the resulting model is well constrained and satisfies the data; and (3) the difficulty of interpreting a model described by an irregular distribution of parameters. Ray tracing and inversion become more problematic when a parametrization is irregular, and are likely to require additional computational effort for the same number of unknowns. However, this trades off against the fewer nodes that will be required to satisfy the data. In large problems, positioning nodes manually according to ray path concentration or some other measure of resolving power can be extremely time consuming. A model composed of basis functions with variable length-scales imposes an additional burden in the interpretation stage because the minimum wavelength of permitted structure will vary according to the distribution of parameters.

In the past, irregular parametrizations have generally only been used in problems involving relatively few unknowns within a static framework. That is, the parameters are positioned by the user a priori and are not repositioned or refined by the inversion step. The commonly used 2-D wide-angle inversion method of Zelt & Smith (1992) is an example of this. They employ a layered, variable block size representation with velocity variation and interface depth controlled by nodes that may be positioned irregularly by the user to improve the trade-off between data fit and resolution. Several global tomography studies have also used static irregular parametrizations. Bijwaard et al. (1998), Bijwaard & Spakman (2000) and Spakman & Bijwaard (2001) use a spatially variable cell size parametrization based on ray sampling. An underlying regular grid is used to construct a mosaic of non-overlapping irregular cells. Sambridge & Gudmundsson (1998) propose a more sophisticated scheme for irregular parametrization based on Delaunay and Voronoi cells and illustrate how it may be used in whole-earth tomography.

Recently, several studies have investigated data adaptive parametrizations for tomographic inversion. This class of scheme does not necessarily impose an irregular grid a priori, but allows the nodes or blocks to adapt to the structural information extracted from the data by the inversion. Michelini (1995) formulated an adaptive-grid method for cross-hole tomography using cubic B-splines in parametric form to describe a 2-D velocity field. From an initially regular grid of nodes, the scheme progresses by simultaneously inverting for both the velocity and the position of the nodes. The aim of the approach is to locate the optimum grid configuration for the model. Curtis & Snieder (1997) opt for a slightly different approach to adaptive cross-hole tomography; they use a genetic algorithm to search for the configuration of Delaunay triangles that results in the best conditioned inverse problem (as measured by the condition number of the matrix to be inverted).

Vesnaver et al. (2000) use an adaptive scheme in 3-D reflection tomography with structure represented by subhorizontal layers. The velocity field within a layer is described in terms of Voronoi polyhedra, which are adjusted during the inversion. This adjustment involves increasing the grid density where the velocity gradient and ray density is significant, and decreasing it where velocity gradients are negligible or ray coverage is poor. Böhm et al. (2000) use a similar scheme to investigate the use of Delaunay triangles and Voronoi polygons in 3-D adaptive tomography (see also Böhm et al. 1996). A dual criterion is used for adding and removing cells from the model based on the magnitude of the velocity gradient and the null space energy, as computed by a singular-value decomposition of the tomographic matrix. The scheme is demonstrated on the SEG/EAEG salt model using reflection and cross-well data in a joint inversion.

In the context of global tomography, Sambridge & Faletic (2003) use an initial model composed of a coarse grid of Delaunay tetrahedra, which are then locally refined after each inversion step according to some measure of the structural signal extracted from the data. Chiao & Kuo (2001) implement a multiscale tomographic parametrization based on 2-D spherical wavelets with permitted bandwidths that are locally controlled by the resolving power of the data. They apply the method to S–SKS traveltimes to image lateral shear wave structure in the D′′ layer.

Both static and adaptive irregular grid tomography show promise, but further effort is required towards understanding more precisely what their advantages and disadvantages are compared with regular grid tomography. In this paper, we describe the benefits and pitfalls of exploiting a static irregular interface parametrization in 3-D wide-angle (i.e. refraction and reflection) traveltime tomography. The scheme we adopt describes an interface by a mosaic of cubic B-spline surface patches. Refracted and reflected rays are traced from source to receiver using a shooting method, and a subspace inversion method is used to invert the traveltimes for interface structure. A series of tests are performed to compare the reconstruction of synthetic models using regular and irregular parametrizations in the presence of highly variable ray coverage. The effectiveness of resolution estimates from linear theory in examining these two classes of solution is also investigated. An application of irregular parametrization tomography to 3-D wide-angle data collected in Tasmania is then presented to help verify the conclusions drawn from the synthetic examples.

2 Method

The method of model parametrization, ray tracing and inversion we use is described in Rawlinson et al. (2001a), so we only provide a brief summary below. However, we will discuss the characteristics of the irregular parametrization and our approach to analysing solution robustness in some detail, as they play an important role in our tests with synthetic and observational data.

2.1 Model Parametrization

Model structure is represented by a set of subhorizontal layers. Within a layer, the velocity varies linearly with depth. Each layer interface is described by a mosaic of cubic B-spline surface patches, the shape and size of which are controlled by a grid of depth nodes with a rectangular association. Thus, for a given set of control vertices p_i,j= (x_i,j, y_i,j, z_i,j), where i= 1, …, m and j= 1, …, n, the (i, j)th surface patch is

(1)

where {b_i} are the uniform cubic B-spline basis functions (Bartels et al. 1987). Phantom vertices are used to extend the surface out to the boundary nodes (Rawlinson et al. 2001a). Any point on a surface patch is a function of the two independent variables u (0 ≤u≤ 1) and v (0 ≤v≤ 1). A surface defined by eq. (1) for any given set of nodes will be everywhere continuous in curvature, and will not necessarily interpolate the control vertices. The surface may also be multivalued, but this feature is not utilized by the inversion method.

The property of eq. (1) that we wish to exploit in this paper is its ability to deal with an irregular distribution of nodes—a smooth (i.e. C₂ continuous) surface will always be constructed. However, the scheme has two relatively minor drawbacks. First, the grid of nodes that control the shapes of the surface patches have a rectangular association. This means that some node configurations may result in surface patches with undesirable shapes. Secondly, there must be m×n nodes, which may lead to some redundancy. Fig. 1 shows several examples of how the node distribution influences the shape of surface patches. In Fig. 1(a), nodes are placed an equal distance apart on a 6 × 7 rectangular grid. Surface patch corners coincide with horizontal node positions in this case. Fig. 1(b) compresses the inner grid of 4 × 5 nodes so that they lie within a 20 × 20 km² area. The nodes are more randomly positioned in Fig. 1(c), but their rectangular association is still evident. Fig. 1(d) has the same node configuration as Fig. 1(a) except that the nodes along the line y= 0 have been placed along three sides of a square at y > 100 km. Similarly, the nodes along y= 120 km have been placed along three sides of a square at y < 20 km. The surface now has highly variable patch sizes and self-intersects; this example demonstrates why nodes cannot be completely randomly positioned.

1

Examples of how node distribution influences surface patch shape using the B-spline formulation of eq. (1). Nodes (grey dots) and surface patches (with boundaries denoted by solid lines) are projected on to the xy plane. Refer to the text for a description of each case.

Open in new tabDownload slide

2.2 Inversion Scheme

Rays are traced through the model by a shooting method—initial-value ray tracing is rapid since ray paths consist of piecewise circular arc segments for which analytic expressions are calculated. A multistage fan shooting scheme is used to target the receiver region and the final two-point path is found using a Gauss–Newton-type iterative method. Multiple two-point paths for a single phase may be found, in which case the path with minimum traveltime is selected.

The inverse problem can be formulated as one of minimizing an objective function consisting of a data residual term and one or more regularization terms. In our case, the objective function we choose is

(2)

where d_obs are the observed data and g(m) are the model predictions; the data residual term is weighted by the data covariance matrix C_d. The regularization term is a function of the difference between the current model (m) and the initial model (m₀) and is weighted by the a priori model covariance matrix C_m. ε is a damping factor that governs the trade-off between how well the data are fitted and how near the final model is to the initial model.

The inverse problem is solved iteratively using a subspace inversion method that only requires the solution of a small system of linear equations and naturally deals with different parameter classes (Kennett et al.1988). To address the non-linear nature of the inverse problem, rays are retraced after each iteration. All refraction and reflection data can be simultaneously inverted for all three classes of model parameter (interface node depth, layer velocity and layer velocity gradient).

2.3 Analysis Of Solution Robustness

A common way of investigating solution non-uniqueness in seismic tomography is to examine the formal estimates of resolution given by the resolution matrix R (e.g. Aki et al. 1977; Menke 1989; Benz et al. 1992; Graeber & Asch 1999). For the objective function of eq. (2), the resolution matrix is defined by

(3)

where G=∂g/∂m represents the Fréchet matrix of partial derivatives calculated at the solution point (see Tarantola 1987 for a derivation). The diagonal elements of R range between 0 and 1; ideally, a value of unity indicates that the model parameter has been uniquely determined by the data; a value of zero indicates that the model parameter is completely unconstrained by the data; and a value between zero and one indicates that the data partially constrains the model parameter. Since geometrical ray theory is used in the calculation of G, eq. (3) does not account for the finite-frequency effects of seismic waves, which impose a limit on the length-scale of objects that may be resolved, regardless of data coverage. In the examples we present, the node separations on the interface grids are always larger than the seismic wavelength, which means that finite-frequency effects are not significant.

The resolution matrix is a function of both data coverage and the model parametrization that has been chosen. For example, if a model were defined by a single parameter (e.g. the average velocity of a volume or the average interface depth), then it is likely that the resolution matrix will yield a value in the vicinity of unity, because the parameter will be sampled by all the data, and hence its value will be extremely well constrained. On the other hand, the spatial resolution of the model will be poor as it permits no structural variation. As the number of parameters is increased, model variance becomes poorer (uncertainties associated with model parameters will increase) but spatial resolution improves since finer features can be detected. This describes the well-known trade-off between resolution and variance (see Backus & Gilbert 1968; Menke 1989). Tomographic schemes that use an irregular parametrization exploit this trade-off differently to schemes that use a regular parametrization.

With a regular parametrization, spatial resolution is kept at the same level throughout the model, which means that model variance will change according to data coverage. In contrast, an irregular parametrization can be used to vary the spatial resolution so that model variance does not change significantly across the model. By quantitatively measuring the level of constraint imposed on each model parameter by the data, the resolution matrix is potentially a valuable tool for assessing tomographic models described by both regular and irregular grids.

It is important to remember that the resolution matrix does not measure: (1) the accuracy of the model—the model parameters may be well constrained by the data, but may not be able to represent the true structural variations or (2) whether the model satisfies the data. The last point is important, because it means that the resolution matrix cannot be used exclusively to analyse solution quality. For our examples, the ability of the model to satisfy the data is measured by the rms traveltime misfit (the rms difference between all observed and predicted traveltimes). The optimum model we seek is one that satisfies the data and is well constrained throughout.

To help analyse the results of the synthetic tests, the similarity between the recovered model and the true model is examined using two measures of misfit: the correlation coefficient and rms vertical separation. The correlation coefficient varies between −1 and 1; a correlation of 1 means that the two surfaces are identical; a correlation of 0 means that they have no similarity; and a correlation of −1 means that the two surfaces are inverted copies of each other. Rms vertical separation is a measure of the average vertical separation distance between the two surfaces. Differences in horizontal extent between surfaces are accounted for in the calculation of both misfit measures, which allows a direct comparison of these values for different solutions.

3 Synthetic Tests

3.1 Results

A series of six synthetic tests using a single interface and fixed layer velocities are performed. An array of five sources and 142 receivers was used to generate the synthetic data (see Fig. 2a). While this particular array geometry is unlikely to be encountered in any observational study, it does produce a highly heterogeneous ray coverage in a way that makes it clear which parts of the model are sampled by data and which parts are not; this makes the results easier to interpret.

2

Interface structure and source–receiver array used in the synthetic tests. (a) Interface nodes and surface patch boundaries are indicated by dots and lines, respectively. Sources are denoted by stars and receivers by triangles. (b) Projection of the surface defined by the nodes in (a). Sources are denoted by large dots and receivers by small dots. Figs 2–9 may be viewed in colour in the on-line version of the journal ().

Open in new tabDownload slide

The model chosen for the synthetic tests is defined by a rectangular grid of 17 × 21 = 357 nodes spaced 10 km apart (see Fig. 2a). The B-spline surface described by these nodes varies between 3 and 12 km in depth (see Fig. 2b); above the interface, the velocity gradient is 0.05 s⁻¹ and below the interface, the velocity gradient is 0.06 s^{− 1}. The average velocity contrast across the interface is 1.1 km s⁻¹ and the average velocity in the top layer along the interface is 4.6 km s⁻¹. A total of 537 refracted and 260 reflected traveltimes comprise the synthetic data set; the MPEG animation referenced in (electronic copy only) shows these rays traced through the interface surface and clearly illustrates the inhomogeneous nature of the data distribution. To simulate the noise content of real data, Gaussian noise with a standard deviation of 45 ms was added to the synthetic traveltimes.

All inversions were carried out using six iterations of a 4-D subspace scheme with a damping factor of ε= 1.0 and an initial model consisting of a horizontal planar interface at 7 km depth. The damping factor controls the trade-off between model perturbation and data fit. Tests with different values of the damping factor show that a value of ε= 1.0 tends to produce the optimum trade-off (see Rawlinson et al. 2001a, for more details on damping).

The first test we perform (see Fig. 3 and Table 1) attempts to recover structure using a grid with the same number and spacing of nodes as was used to describe the synthetic model. The rms data misfit curve (Fig. 3a) decreases monotonically from an initial value of 267 ms to a final value of 44 ms, which is approximately equal to the standard deviation of the added Gaussian noise. This behaviour is characteristic of a stable inversion. The horizontal node distribution and ray-interface hit points of the solution model are shown in Fig. 3(b)—it is evident that a large number of surface patches are not intersected by any rays. The resolution diagram (Fig. 3c), which plots the diagonal elements of the resolution matrix at the solution point, strongly correlates with the ray path distribution. A discontinuity in the colour palette is introduced at the resolution value of 0.5 to help distinguish between well-constrained regions (≥0.5) and poorly constrained regions (<0.5). Of course, such a separation is rather arbitrary and may vary depending on the value chosen for the damping parameter. However, as a general guide, we have found that this cut-off level is a reasonable indicator of which parts of the interface are well constrained by the data and which parts are not.

3

Results from the first test. Structure is recovered only in regions of good data coverage. (a) Rms data misfit versus iteration number. (b) Node distribution and ray-interface hit points. Nodes are denoted by dots and surface patch boundaries are denoted by lines. Refracted ray hit points are denoted by circles and reflected ray hits by asterisks. (c) Diagonal elements of the resolution matrix. (d) True model interface structure. (e) Recovered model interface structure. (f) Misfit between the true model and the recovered model. Triangles represent receiver locations and stars represent source locations in all plots.

Open in new tabDownload slide

Table 1.

Data and model misfit information for the six synthetic tests. The i= 0 misfit refers to the rms data misfit for the initial model; i= 6 misfit refers to the rms data misfit for the solution model. Note that Figs 6 and 7 use a different synthetic data set to the other four tests.

Open in new tabDownload slide

Fig. 3(d)–(f) compare the true model, recovered model and model misfit (i.e. the difference between the true model and the recovered model), respectively. The resolution plot (Fig. 3c) appears to be an excellent indicator of which parts of the model have been recovered. In regions of poor resolution (<0.5), the solution model tends not to deviate from the initial model, but in regions of good resolution (≥0.5), structure is accurately recovered. In particular, Fig. 3(f) shows that there is virtually no leftover structure within the inner ring of receivers where resolution is especially good.

The patchy nature of the solution model (Fig. 3e), caused by poorly constrained parts of the interface, is symptomatic of schemes that map information from inhomogeneously distributed data into models with regular parametrizations. A common way of treating this effect is to apply smoothing (e.g. Sambridge 1990). Second-derivative smoothing is equivalent to increasing the node spacing with our parametrization. The second test (see Fig. 4) is similar to the first except that it uses only 120 nodes on a regular grid, a 66 per cent reduction in the number of unknowns. As one might expect, the resultant solution (Fig. 4e) essentially looks like a smoothed version of the true model in the well-constrained regions.

4

Results from the second test. Only broad scale features are recovered in this case. See the caption of Fig. 3 for an explanation of each plot.

Open in new tabDownload slide

A comparison of the resolution plots of Figs 3(c) and 4(c) shows that the interface is well constrained in a larger region of the Fig. 4 solution—this is to be expected since there are much fewer nodes but the same volume of data. However, this does not mean that the second test model is superior. In fact, the traveltime residual at iteration six is now 74 ms compared with 44 ms for the first test, and Fig. 4(f) shows that there is significant short-wavelength structure that has not been recovered. This is also reflected in the poorer correlation coefficient and rms model misfit (Table 1) for this solution compared with the Fig. 3 solution. Nevertheless, the improved area of resolution is meaningful provided it is understood that structure is recovered in these regions only at the scalelengths permitted by the parametrization. As an example of this, we see that the shallow anomaly at (x < 60, y= 80) km, which is oriented in the x-direction, is both better recovered and better constrained in the second test (Fig. 4) compared with the first (Fig. 3). Another advantage of the second test is that it only required 49 per cent of the CPU time of the first test.

The third test (Fig. 5) uses the same number of nodes as the second but spaces them irregularly. Nodes are placed in the vicinity of ray-interface hit points (Fig. 5b), with the greatest concentration of nodes occurring within the inner ring of receivers where ray paths are most numerous. The resolution plot (Fig. 5c) indicates that this configuration of nodes results in a solution model that is well constrained within the entire region spanned by the rays. The fine-scale structure in the vicinity of the sources is accurately retrieved, and where data coverage is sparse, the broader-scale features of the true model have been extracted. Both measures of model misfit (Table 1) suggest a superior reconstruction compared with the first and second tests. The data fit of 52 ms is a significant improvement over the solution of Fig. 4, and the CPU time is only slightly more at 54 per cent of the first test.

5

Results from the third test. The irregular parametrization is adapted to the data coverage. See the caption of Fig. 3 for an explanation of each plot.

Open in new tabDownload slide

The fourth (Fig. 6) and fifth (Fig. 7) tests are similar to the first (Fig. 3) and third (Fig. 5) tests except that the synthetic structure we wish to recover is now defined on the irregular grid of Fig. 5(b). The fourth test attempts to recover the structure using the same irregular grid, while the fifth uses the 357 node regular grid of Fig. 3(b). The reconstruction of Fig. 6 is clearly superior to that of Fig. 7, even though they both have identical data misfits of 45 ms. The correlation coefficients and rms vertical separations for these two tests also support this assertion.

6

Results from the fourth test. The true model is now described by an irregular grid. The recovered model is based on this same irregular grid. See the caption of Fig. 3 for an explanation of each plot.

Open in new tabDownload slide

7

Results from the fifth test. As in the fourth test, the true model is described by an irregular grid, but the recovered model is now based on a regular grid. See the caption of Fig. 3 for an explanation of each plot.

Open in new tabDownload slide

The sixth and final synthetic test (Fig. 8), is a caveat on the use of irregular parametrizations. The structure we attempt to recover is once again the regular grid model (Fig. 3d), and a 120-node irregular grid (Fig. 8b) is used to satisfy the data. Nodes are placed only in the vicinity of ray interface hits but in a different configuration to that of Fig. 5(b). The inversion solution in this case is much poorer than that of Fig. 5; in fact, as also indicated by the model misfit measures, it is even worse than the regular parametrization of the second test (Fig. 4) despite having a slightly better data misfit (67 ms compared with 74 ms). The resolution plot (Fig. 8c) indicates that this model is better constrained by the data than any of the previous models, despite the poor recovery.

8

Results from the sixth test. An highly irregular grid is used to try and capture structural information described on a regular grid. See the caption of Fig. 3 for an explanation of each plot.

Open in new tabDownload slide

3.2 Discussion

The above synthetic tests succinctly reveal the drawbacks of using a regular parametrization in the presence of a highly heterogeneous data distribution. These include the presence of unnecessary nodes, which do not affect data fitting, increased computing time and a less well-determined inverse problem, as reflected in the resolution plots. Fig. 3 results show that a fine-scale regular grid accurately reconstructs the interface in regions of good data constraints, but elsewhere, the interface tends not to deviate from its initial value. Using a coarser regular grid, as in Fig. 4, results in a larger area of good data constraint, but a poorer data fit. This is to be expected because only the longer-wavelength features are recovered, and therefore the ray sampling of the interface need not be as intense to achieve good constraints on the model parameter values. The poorer data and model fit (Table 1) are due to the finer-scale structure not being recovered (Fig. 4f).

The results of the first (Fig. 3) and second (Fig. 4) tests suggest that a better outcome may be obtained by using a finer-scale parametrization in regions of good data coverage and a broader-scale parametrization in regions of poor data coverage—this is demonstrated in Fig. 5. The results illustrated in Figs 6 and 7 complement those of Figs 3 and 5 and provide another demonstration that the irregular parametrization is capable of recovering structure in regions of poor data coverage where the fine-scale regular parametrization does not deviate from the initial model. The regular parametrization solution is conservative because structure is only introduced where there are constraints on the minimum wavelength permitted by the parametrization; elsewhere, structure tends not to be introduced even though the data are capable of resolving broader-scale features. Put another way, there will always be some wavelength of structure that the data will constrain, no matter how inhomogeneously it is distributed. Regular parametrizations are not consistent with this reality.

Although we have highlighted the advantages of an irregular parametrization, the final synthetic test (Fig. 8) is clear evidence that care must be taken in choosing an appropriate node distribution. In this case, nodes are placed only where the concentration of ray-interface hit points is high; consequently, the solution is extremely well constrained by the data within the region bounded by the source–receiver array. One therefore may be tempted to conclude that all the structure contained in this model is ‘real’. However, comparison of Figs 8(d)–(f) shows that this is not the case. The other two indicators of something being amiss with the solution are the data fit and the nature of the recovered anomalies.

The rms data misfit of the final test solution is much worse than that of the Fig. 3 test and the Fig. 5 test and only marginally better than the second test (Fig. 4). The latter result seems to suggest that Fig. 8(e) is a more desirable model than Fig. 4(e), which is not actually the case. The marginal improvement in misfit of the irregular model can be attributed to the improved reconstruction within the inner circuit of receivers, which contains many more parameters than the coarse regular grid. Therefore, the data misfit for paths that travel to the inner circuit of receivers is much smaller for the irregular grid solution (55 ms) than it is for the regular grid solution (77 ms). For paths to the outer array of receivers, the data misfit is 75 ms for the irregular case and 72 ms for the regular case. This is significant because the paths to the outer receiver array are largely responsible for constraining structure that lies outside the inner circuit of receivers; for the irregular solution, this is the region where reconstruction is poor (Fig. 8f).

The shapes of the anomalies in Fig. 8(e) tend to bear a strong resemblance to the shapes of the underlying surface patches (Fig. 8b). When this occurs, it suggests that the data contain information on smaller-scale structure which is being filtered out by the parametrization choice. In some ways, this is not a bad thing—in fact, something similar occurs when coarse regular parametrizations are used. However, with this irregular node distribution, the minimum wavelength of structure that is permitted varies strongly both spatially and directionally. For example, the set of patches to the left and right of the inner ring of receivers are oriented at approximately 45° to the x and y axes. Structure in these regions (Fig. 8e) is smeared along the long axis of these patches, but perpendicular to this direction, the pattern of anomalies is quite well recovered.

From the plots of data coverage for our synthetic tests (e.g. Fig. 5b), it is clear that there are large areas of the interface that do not have any ray penetration. Yet we apparently claim, on the basis of data misfit and resolution estimates, that in the general case, structural information in these regions can somehow be extracted from the data set. For example, one may argue that in regions of ‘no’ data coverage, a wide variety of structures could be introduced without affecting any ray path or traveltime. Essentially, this is true, but it must be remembered that the role of any parametrization is to regularize (e.g. by smoothing) the solution everywhere between infinitely thin data (under the high-frequency assumption). In other words, the inverse problem is inherently underdetermined and strictly speaking, structural information is really inferred rather than extracted. Consequently, the accuracy of a solution (i.e. how close it is to the truth) will also depend on our underlying assumptions concerning the type of structures that can exist, in addition to considerations of model robustness and data fit.

In our examples, we have assumed that layer boundaries vary smoothly, and that structure can be decomposed into features of different length-scales. While this assumption is reasonable, it will not always be true; however, this could also be said for any parametrization that one might care to choose. From this point of view, all tomographic solutions must be treated with some caution, whether they are based on regular or irregular parametrizations. The approach we adopt, which varies the length-scale of structure that is recovered according to data distribution, is consistent with the underlying assumptions of our parametrization. As such, we could claim that it produces a more robust and desirable solution compared with a regular parametrization, but it is difficult to guarantee that it will always be more accurate.

4 Application To The Tasmanian Data Set

4.1 Results

In this section, we examine the use of an irregular interface parametrization in the inversion of wide-angle data collected in Tasmania, SE Australia in 1995. To create this data set, the Australian Geological Survey Organization (AGSO) research vessel Rig Seismic fired ∼36 000 air-gun shots at a 50 m spacing during a circumnavigation of Tasmania. An array of 44 digital and analogue recorders were deployed throughout Tasmania to record the data. Rawlinson et al. (2001b) inverted 2148 P_mP and 442 P_n traveltimes from this data set for the structure of the Moho beneath Tasmania using the scheme of Rawlinson et al. (2001a) with a regular interface parametrization. We aim to repeat this experiment with an irregular grid of interface nodes, and then compare the results with those of Rawlinson et al. (2001b). The procedure we follow is identical to that of Rawlinson et al. (2001b) so we will not describe it in detail here.

The irregular grid used to parametrize the Moho surface is composed of 240 nodes preferentially distributed to regions with good ray coverage (Fig. 9a). A number of nodes in our scheme are essentially redundant (e.g. above line 4 and in the vicinity of line 12)— this is hard to avoid due to the nature of the data distribution and the rectangular association of surface patches. Despite this, we have more than halved the number of nodes compared with the regular parametrization of Rawlinson et al. (2001b), which uses 600 nodes (Fig. 9b). Four iterations of a 14-D subspace inversion scheme are used to produce the solution model from an initially flat interface at 30 km depth. The layer velocities and the layer velocity gradients are permitted to vary during the inversion in addition to the interface depth.

9

Comparison of inversion results using regular and irregular grids and wide-angle data from Tasmania. In all diagrams, stars indicate recorder positions and small triangles indicate shot points from which data were picked—contiguous triangles form shot lines. (a) Ray interface hit points and surface patches for the irregular grid. Refracted ray hits are denoted by circles and reflected ray hits are denoted by crosses. Interface nodes are represented by dots and thin lines represent surface patch boundaries. (b) Same as (a) but for the regular grid. (c) Diagonal elements of the resolution matrix for the irregular grid solution. (d) Same as (c) but for the regular grid. (e) Solution model for the irregular grid. (f) Solution model for the regular grid.

Open in new tabDownload slide

The solution model (Fig. 9e) and resolution plot (Fig. 9c) produced by the inversion with an irregular grid show that significant structure is recovered, but regions of the Moho, particularly beneath central Tasmania, are still poorly resolved. However, if we compare these results with the regular parametrization solution of Rawlinson et al. (2001b), as shown in Figs 9(d) and (f), we see that the area of good solution robustness has been substantially increased.

The initial model has an rms data misfit of 371 ms; this is reduced to 191 ms for the irregular grid solution and 176 ms for the regular grid solution. Thus, there is only a small trade-off in data misfit for the large saving in nodes (60 per cent) and computing time (the irregular solution required only 60 per cent of the computing time of the regular solution). An MPEG animation of the irregular grid solution model and constraining ray paths is referenced in (electronic copy only).

4.2 Discussion

The use of an irregular interface parametrization in the inversion of 3-D wide-angle data from Tasmania demonstrates that we can also improve solution robustness without significantly degrading the data fit. A comparison of Figs 9(c) and (d) shows that the robustness of the irregular parameter solution is generally superior to that of the regular parametrization solution of Rawlinson et al. (2001b). In fact, towards the centre of Tasmania, the resolution has become non-zero. This indicates that the model structure in this region is now influenced by the data to some extent. In comparing Figs 9(e) and (f), we see that the crust is thicker near central Tasmania in the irregular model compared with the regular model—this is consistent with the increased surface elevation in this region. However, the low resolution values associated with central Tasmania (Fig. 9c) mean that we should treat this result with some caution.

Another noticeable difference between the two Moho solution models is that Fig. 9(e) has a less patchy appearance than Fig. 9(f), but at the same time does not appear to have features influenced by the shape of the underlying surface patches. For example, the two shallow anomalies in the vicinity of station 17 (Fig. 9f) have become a single shallow anomaly in the irregular interface solution (Fig. 9e). Ultimately, though, we are not inclined to modify the overall structural interpretation given by Rawlinson et al. (2001b) based on the new irregular model, as they share most of the significant features. However, the irregular solution (Fig. 9e) does reveal a broad scale shallowing in the NE of Tasmania that is both less distinct and less well constrained in the regular solution (Fig. 9f). The location of this feature coincides with that of the Northeast Tasmania Element (Brown et al. 1998), a tectonic element with a distinct surface geology that was juxtaposed with western Tasmania during the Mid-Devonian Tabberabberan Orogeny (Elliot et al. 1993). Evidence from coincident reflection data (Barton 1999) and a previous refraction survey (Richardson 1980) also indicates a thinning of the crust beneath the Northeast Tasmania Element.

5 Conclusions

In this study, we have investigated the use of irregular interface parametrization in 3-D wide-angle seismic traveltime tomography. Our results from synthetic tests indicate that an irregular approach has many potential benefits when compared with a regular approach if data are inhomogeneously distributed. These include faster computation time and the ability to satisfy the data with a better constrained model. We cannot guarantee that such a model will be more accurate, but if the underlying structural assumptions of the parametrization hold true in reality, then the likelihood of this happening will increase.

The synthetic tests clearly show that broader-scale structure can be recovered in regions where a finer regular parametrization is unable to recover structure. The inversion of 3-D wide-angle data from Tasmania for Moho depth also supports these findings: several new broad scale features, which were absent from a previous regular grid inversion, were revealed by the irregular grid solution. At the same time, the well-constrained shorter-scalelength features of the regular grid solution were preserved.

Although we have established that regular parametrizations have undesirable properties in the presence of inhomogeneously distributed data, care must be taken when adopting an irregular scheme because the minimum wavelength of structure permitted in the model will vary both spatially and directionally. One consequence of this is that there is a possibility of obtaining an extremely well-constrained solution that poorly represents the true structure. This type of solution should be revealed by a poor data fit, but another indicator is the presence of structure that conforms to the shape of the underlying basis functions.

For relatively small tomographic inversions such as those carried out in this paper, manual positioning of nodes is feasible. However, for large tomographic problems, such as those which image 3-D velocity variations, data adaptive schemes would be more efficient and desirable. One automated way of finding an optimal parametrization, which is consistent with the findings of this study, is to begin with a coarse parametrization that is constrained everywhere by the data (as measured by the diagonal elements of the resolution matrix). Parameters may then be added to regions where resolution is greatest until a suitable trade-off between model robustness and data fit is found. This remains a direction for further study.

Acknowledgments

We thank A. Vesnaver and an anonymous reviewer for constructive reviews of the original manuscript. Fig. 2(b) and the MPEG animations referenced in the appendices were created using Geomview, a freeware interactive 3-D viewing program (see ). All remaining figures were created using the PGPLOT graphics package, which is also freely available over the internet at. This work was supported by Australian Research Council Discovery Project DP0208039.

Appendix

Appendix A: Visualization Of Synthetic Model

An MPEG movie has been created to help visualize the distribution of ray paths through the synthetic model and the geometry of the interface surface. The movie can be downloaded from. The ray paths shown correspond to those which produced the synthetic traveltimes which were inverted in Figs 3–5 and 8, and the interface surface corresponds to the ‘true model’ that we attempt to reconstruct in each of these four tests.

A frame from the MPEG movie is shown in Fig. A1; note that the variations in depth of the interface have been colour contoured using the same colour scheme as in Figs 3–8. A 3-D interactive visualization program called Geomview is used to visualize the model, which is represented as an OOGL (Object Oriented Graphics Library) object. The actual animation sequence is coded as a TCL script, which is executed by a Geomview module called Stagetools to create the MPEG movie. The Geomview package, which includes the Stagetools module, is freely available over the internet at.

A1

Snapshot from the MPEG movie of the synthetic model. Refracted rays are shown in red; reflected rays are shown in green. Pink dots denote receivers and red dots denote sources. Cf. Fig. 3(d).

Open in new tabDownload slide

Appendix

Appendix B: Visualization Of Tasmanian Moho Model

An MPEG movie of the Tasmanian Moho solution model (Fig. 9e) and its constraining ray paths (Fig. 9a) has been constructed. It can be downloaded from the same web site as the movie described in. A frame from the MPEG movie is shown in Fig. A2; note that variations in interface depth have been colour contoured using the same colour scheme as in Fig. 9(e). The MPEG movie was constructed using the method described in.

A2

Snapshot from the MPEG movie of the Tasmanian Moho. Refracted rays are shown in red; reflected rays are shown in green. Pink dots denote sources and red dots denote receivers. The Tasmanian coastline is shown in purple. Cf. Figs 9(a) and (e).

Open in new tabDownload slide

References