Exploiting Lateral Resolution of Near-Surface Seismic Refraction Methods

INTRODUCTION

Ubiquity of Non-uniqueness

In recent years, refraction tomography (Lanz et al., 1998; Zhang and Toksöz, 1998; Stefani, 1995; Zhu et al., 1992), also known as tomographic inversion, has been widely employed to process near-surface seismic refraction data. The advantages of these methods when compared with traditional approaches include greater speed of processing, minimal requirements with technical proficiency, and attractive methods for presenting the results.

Refraction tomography is one example of modelbased inversion, in which an initial starting model is systematically updated by iteratively comparing the modeled response with the field data. The role of inversion is to provide information about the unknown numerical parameters that go into the model, not to provide the model itself (Menke, 1989). If features such as lateral or azimuthal velocity variations in the main refractor, velocity reversals, or vertical seismic anisotropy in the overburden are not included in the model, then the inversion process will not necessarily provide useful information about those features, even though they may be appropriate. Accordingly, close agreement between the modeled response and the data simply demonstrates that the model is consistent with the data: it does not "prove" that the model is "correct" or even the most probable.

Frequently, it is possible to obtain consistency between the data and a wide range of quite disparate starting models. It is now generally recognized that non-uniqueness is a fundamental reality with the inversion of virtually all sets of geophysical data (Treitel and Lines, 1988; Oldenburg, 1984). Although theoretical studies demonstrate the ubiquity of nonuniqueness with the inversion of near-surface seismic refraction data (Ivanov et al., 2005a, b), nevertheless, there can still be considerable reluctance to address the issue, largely because of the implications for litigation with geotechnical applications.

The most common source of non-uniqueness is the choice of the inversion algorithm used to generate the starting model (Oldenburg, 1984). Non-uniqueness can occur with the 1D parameterization of the layering in the vertical dimension and with the 2D lateral resolution of individual layers in the horizontal dimension.

Inversion Algorithms for Vertical and Lateral Resolutions

The 1D τ-p algorithm (Barton and Barker, 2003) is the default in virtually all refraction tomography programs, because it is commonly considered that this algorithm can achieve automatic parameterization of the traveltime graphs into the various layers. The parameterization of the traveltime graphs is often considered to be a major source of uncertainty in the processing of seismic refraction data.

The τ-p algorithm emphasizes the vertical resolution of many layers, which in most refraction tomography programs are parameterized with vertical linear velocity gradients. However, virtually any standard velocity function can be fitted to most traveltime data with acceptable accuracy (Palmer, 1992, 1986; Hagedoorn, 1955), because the arrivals used to derive the seismic velocities in each layer usually propagate within only a relatively thin section of the upper part of each layer. With reasonable vertical velocity gradients in which the increases are less than a few meters per second per meter, less than 30% of the upper section of each layer is usually sampled (Palmer, 1986).

This incomplete sampling is the major reason for the well-known undetected layer problem, whereby layers, which are below a minimum thickness in relation to the seismic velocities of the surrounding layers, are not detected as first arrivals. Using the analysis of Merrick et al. (1978) for constant velocity layering, Palmer (1986) demonstrated that the number of undetected layers is more significant than the seismic velocities within those layers.

The hyperbolic velocity function (Aki and Richards, 2002; Berry, 1971; Healy, 1963; Slichter, 1932) represents the limiting case in which the maximum vertical velocity gradient, or equivalently, a multiplicity of undetected layers can exist in the region above the main refractor. As a result, traveltime graphs consisting of linear segments can represent either constant velocity layering or vertical hyperbolic velocity gradients. The most common result of model parameterization with vertical velocity gradients is that the depth to the main refractor is increased, often by approximately 50%.

Figure 1 illustrates the ambiguity. Although the ray paths suggest complete sampling of the first layer, the reality is that the traveltime graph represents signals that only travel along the upper region of the first layer, and therefore, have minimal sampling of that layer.

Figure  1.  Schematic traveltime graphs and ray paths for a hyperbolic vertical velocity gradient in the upper layer.

DownLoad: Full-Size Img PowerPoint

Frequently, ray path coverage diagrams are presented as a validation of the outputs of automatic refraction tomographic inversion. However, such ray path trajectories are an inevitable consequence of the vertical velocity function selected to represent each layer, rather than the true seismic velocities within the various layers. Minimal penetration occurs where constant velocity layering is employed. Alternatively, complete sampling of the layer can occur where vertical velocity gradients are employed, as is shown in Fig. 1. Because the fundamental task of determining a unique vertical velocity function in each layer is not realizable with the traveltime data alone, it follows that it is not possible to derive a true or a unique ray path coverage for each layer. Accordingly, it can be concluded that, being circulus in probando, ray path coverage diagrams can be of very limited usefulness in validating the results in the majority of cases.

Another major shortcoming with the τ-p algorithm is its failure to detect large lateral changes with low seismic velocities that can be representative of shear zones. This study demonstrates this deficiency with a major shear zone that is 50 m wide. However, even shear zones in excess of 100 m wide, such as that at Welcome Reef (Palmer, 1980), are not detected with the τ-p algorithm (Palmer, 2008a).

Alternative inversion algorithms include the 1.5D algorithms of Hagiwara's method (Hagiwara and Omote, 1939), the plus-minus method (Hagedoorn, 1959), the ABC method (Nettleton, 1940), and the conventional reciprocal method (Hawkins, 1961), which are all mathematically identical, and the 2D algorithms of the generalized reciprocal method (GRM; Palmer, 2006, 1986, 1981, 1980; Palmer et al., 2005). While 1.5D and 2D methods generate much the same time structure model of the refractor (see Fig. 12), the use of 2D methods is indicated where detailed lateral variations in the seismic velocities in the refractor are required. Commonly, 1.5D methods generate artifacts in the seismic velocities in the refractor where there are irregular (2D) refractors (Palmer et al., 2005; Palmer, 1986).

The major objective of this study is to demonstrate that the algorithms for processing seismic refraction data that emphasize lateral resolution generate more useful results than the algorithms that emphasize vertical resolution. This approach is consistent with the observation that most refracted seismic energy propagates in a predominantly horizontal direction, rather than a predominantly vertical direction as is the case with reflected seismic energy. Furthermore, this approach is supported by the fact that the unique parameterization of the layering in terms of constant seismic velocities, vertical velocity gradients, or velocity reversals is not possible using the traveltimes from those layers, as is shown in Fig. 1.

In addition to the use of traveltime data, this study also employs the head wave amplitudes and the refraction convolution section (RCS; Palmer and Shadlow, 2008; de Franco, 2005; Palmer and Jones, 2005; Palmer, 2001a, b) in order to extract additional useful information from the field data. The derivation of additional parameters with which to characterize the regolith can promote more effective application of near-surface seismic refraction operations, many of which have exhibited only quite modest innovation in more than half a century (Palmer, 2008b).

The Mt. Bulga Case Study

Mt. Bulga is located near Orange in southeastern Australia, and is the site of a small massive sulphide ore body in steeply dipping to vertical Silurian metasediments. The site underwent extensive investigation several decades ago, and much of the exploration program is summarized in Whiteley (1986), who also recorded numerous seismic refraction profiles at the site.

The 2D seismic profile described in this study (Palmer, 2001a, b) was recorded by the author as part of an undergraduate geophysics students' field tutorial at the site in the 1980s. The profile is located some distance away from the ore body on relatively flat ground in the vicinity of the now-demolished core shed, and corresponds with in-line 19 in Fig. 22. Note that the in-line station numbers correspond with the cross-line line numbers, and conversely, the in-line line numbers correspond with the cross-line station numbers. Nine shot points with a separation of 60 m were recorded.

Single digital grade geophones with a natural frequency of 40 Hz were used. The response below the natural frequency is approximately 12 dB/octave. Figure 21 shows that the signal below approximately 20 Hz is strongly attenuated. The anti-aliasing filters for the one millisecond sampling interval limited the high frequency response to approximately 250 Hz.

The traveltime data, which are presented in Fig. 2, show that three layers can be recognized. They are a surface soil layer with a seismic velocity of 300–500 m/s, a weathered layer with a seismic velocity of ~1 000–1 200 m/s for which the upward concave graphs suggest cycle skipping, and in turn, the occurrence of a velocity reversal, and a main refractor at the base of the weathering with a seismic velocity of 2 500–6 000 m/s. The traveltime graphs for shot points 1, 13, 85, and 97 represent arrivals from the base of the weathering, and since they are very irregular, they suggest that 2D rather than 1D inversion algorithms would be more useful.

Figure  2.  Traveltime graphs for the Mt. Bulga 2D traverse.

DownLoad: Full-Size Img PowerPoint

Approximately one year later, a 3D refraction survey was recorded at the same location using the same 48-channel seismic recorder, in order to study azimuthal anisotropy associated with the shear zone (Palmer, 2001c). The 3D survey consisted of two separate receiver setups.

The first used two lines of 24 receivers at a 5-m spacing, which were parallel to, and located approximately 10 m either side of the earlier 2D survey and correspond with in-lines 17 and 21 in Fig. 22. The in-line station numbers of these two in-lines correspond with those of the earlier 2D profile. Five shot points, 60 m apart, were located along each line. Another four oblique shot points offset 60 m from the end of each line of receivers in the in-line direction and offset 60 m in the cross-line direction were also recorded, making a total of fourteen shot points.

The second setup used a series of seven parallel cross-lines that were 20 m apart, and each of which consisted of 12 receivers at a 5-m spacing. The crossline numbers are 45, 49, 53, 57, 61, 65, and 69 in Fig. 22. The center of each cross-line coincided with the earlier 2D profile. There were four shot points nominally 60 m apart on each cross-line. These lines were recorded in groups of four by simply rolling through from one end to the other. A total of 27 shot points were recorded in the cross-line direction.

Tomographic Inversion with the τ-p Algorithm

Figure 3 presents the starting model generated with the τ-p algorithm and the tomogram obtained with wave eikonal traveltime (WET) tomography (Schuster and Quintus-Bosz, 1993) through its implementation with the refraction tomography program RAYFRACT. The close similarity between the starting model and the final WET tomogram emphasizes the significance of the inversion algorithm used to generate the starting model. The errors are presented in Table 1.

Figure  3.  The starting model and the WET tomogram obtained with the 1D τ-p inversion algorithm.

DownLoad: Full-Size Img PowerPoint

Table  1.  WET tomography errors

 | Show Table

DownLoad: CSV

The τ-p tomogram does not demonstrate the existence of any low-velocity shear zones. The seismic velocities in the lowest layers are greater than 4 000 m/s for the entire traverse and there are no regions with seismic velocities that would be indicative of shear zones. The region between stations 40 and 50 is interpreted as a local increase in the depth of weathering and not as a shear zone because the seismic velocity in the lower part of the tomogram is still 4 000 m/s or higher.

Tomographic Inversion with GRM Algorithms

GRM inversion algorithms

The GRM employs two inversion algorithms.

The first, the refractor velocity analysis function t_V (Palmer, 1980, equation (2); Palmer, 1981, equation (1); Palmer, 1986, equation (8.1)), is shown in equation (1)

where t_FY is the traveltime between the forward shot point at F and the receiver at Y; t_RX is the traveltime between the reverse shot point at R and the receiver at X; and t_FR is the reciprocal time between the two shot points.

Equation (1) is illustrated schematically in Fig. 4, in which traveltimes that are added are shown as solid lines, those that are subtracted are shown as open lines, and both are shown as broken lines where they cancel. Where detailed lateral resolution of seismic velocities is required, equation (1) is evaluated for a range of XY separations between the reverse and forward receiver positions, and is referenced to the position G, where G is midway between X and Y, that is XG=GY.

Figure  4.  Schematic ray paths used with the GRM refractor velocity analysis algorithm.

DownLoad: Full-Size Img PowerPoint

From Fig. 4, it can be seen that equation (1) reduces to approximately the traveltime from the shot point at F to a point in the refractor below the reference point G. Therefore, the gradient of the line(s) fitted to the points computed with equation (1), using equation (2), is a measure of the seismic velocity in the refractor, viz

The second is the time model algorithm tG (Palmer, 1980, equation (10); Palmer, 1981, equation (2); Palmer, 1986, equation (8.6)), which provides a measure of the depth to the refracting interface in units of time, and which is shown in equation (3)

The term XY/V_n represents the additional traveltime in the refractor between the stations at X and Y. This term is approximated with the difference in the refractor velocity analysis function in equation (1), which has been averaged over a range of XY values, i.e.

Equation (3) is illustrated in Fig. 5. It shows that the GRM time model is the average of the reverse delay time at X and the forward delay time at Y.

Figure  5.  Schematic ray paths used with the GRM time model algorithm.

DownLoad: Full-Size Img PowerPoint

The time model is related to the thicknesses Z_i, and the seismic velocity V_i, of each layer in the weathering with equation (5)

In many cases, however, it is not possible to define all layers within the weathering because of large station spacings, as frequently occurs with reflection surveys for deep targets, or velocity reversals, as frequently occurs in arid and semi-arid regions. In such circumstances, the multiplicity of layers in equation (5) can often be replaced with a single layer of total thickness Z_G and an average seismic velocity V, as shown in equation (6), viz

A unique feature of the GRM is the computation of an average vertical seismic velocity in the weathering where an optimum XY value can be recovered from the refractor velocity analysis function, as shown in equation (7) (Palmer, 1980, equation (27); Palmer, 1981, equation (10); Palmer, 1986, equation (10.4)). The optimum XY value is selected where the seismic velocity model is the simplest, that is, the extreme lateral variations are minimized

An objective methodology

In the standard application of the GRM, a single velocity model, the so-called optimum, is derived firstly by visually determining the simplest of the refractor velocity analysis functions computed with equation (1) for a range of XY separations and then with manual curve fitting for the evaluation of equation (2) over discrete intervals (Palmer, 1991). However, Sjögren's (2000) critique of the GRM demonstrates the necessity of employing objective approaches for determining seismic velocities. Despite the fact that the mean-minus-T method is mathematically equivalent to the GRM (Sjögren, 2000), and that Sjögren's Fig. 4 is identical to Fig. 17 in Palmer (1991), apart from some minor smoothing, nevertheless Sjögren achieved a different model of the seismic velocities in the refractor. Sjögren's (2000) critique demonstrates the significance of the subjectivity in simple curve fitting, rather than his conclusion that the GRM is somehow defective.

In this study, the seismic velocities are obtained with a simple algorithm that computes a reciprocal of the gradient with equation (2) at each station. The aim is to employ an objective methodology for automatically computing a multiplicity of detailed seismic velocity models in the refractor. If an evaluation of the significance of non-uniqueness is an objective, then each of these models can be used as starting models for tomographic inversion (Palmer, 2008a). The algorithm employed in this study to determine the seismic velocity at each station is shown in equation (8)

Corrections for variations in soil composition

The traveltimes include minor variations that originate in the surface soil layers and from minor picking errors. These minor variations, which can be correlated with the short wavelength or residual static corrections in the processing of seismic reflection data (Palmer and Shadlow, 2008), must be removed prior to the application of equation (8). Otherwise, they can adversely affect the detailed determination of seismic velocities in the base of the weathering. This study employs a smoothing method based on the GRM (Palmer, 2006). The corrections for the traveltimes with shot points at stations 1 and 97 are shown in Fig. 6 and average 0.02 ms with a standard deviation of ±0.88 ms.

Figure  6.  The GRM smoothing time and amplitude corrections for variations in the surface soil layers for SP1 and SP97.

DownLoad: Full-Size Img PowerPoint

Although the application of the GRM smoothing corrections has significantly reduced most of the "noise" from the surface soil layers, it has not completely removed it. Although not shown here, the application of equation (8) results in seismic velocities that are still erratic and unrealistically high, especially between stations 25 and 49. In that region, the seismic velocities in the refractor are approximately 6 000 m/s, and the increments in the refractor velocity analysis function are less than one millisecond per station interval. As a result, even small errors of a tenth of a millisecond in these times can result in large variations in the seismic velocities computed with equation (8). Although another application of the GRM smoothing corrections would reduce these errors, it would also introduce additional smoothing and a possible loss of resolution.

Figure 7 shows the refractor velocity analysis function computed with equation (1) for XY values which range from -20 to 25 m in increments of 5 m [-20 m ≤ XY (5 m) ≤ 25 m] after one application of the GRM smoothing corrections to the traveltimes from shot points at stations 1 and 97. The major effect of the application of the smoothing corrections is that the refractor velocity analysis function is more symmetrical about the optimum XY value, which an inspection indicates is between 0 and 5 m. This study will use a value of 2.5 m, which illustrates the resolution of the optimum XY value to half the trace spacing. Previous studies (Palmer, 2003, 2001a, b) have used a 5 m value.

Figure  7.  The GRM refractor velocity analysis function for -20 m ≤ XY (5 m) ≤ 25 m together with the average.

DownLoad: Full-Size Img PowerPoint

Averaging refractor velocity function

The approach employed in this study is based on the empirical observation that the refractor velocity analysis function is essentially symmetrical about the optimum value (Palmer, 1980, Fig. 26), as is shown in Fig. 7. Therefore, averaging the refractor velocity analysis function over a range of XY values centered on the optimum preserves the shape of the graph, while also reducing any small errors through the averaging process.

Figure 8 presents the seismic velocities computed with equation (8) with ∆x=5 m, averaged over various ranges of XY values that are centered on an XY value of 2.5 m. For small ranges of XY values, the seismic velocities oscillate about an average, which is better defined where the averaging is applied over larger ranges of XY values. Although the averaging over 10 values, which represents XY values from -20 to 25 m, has been selected as providing optimum smoothing, nevertheless, some values are as large as 8 000 m/s which are clearly not geologically reasonable.

Figure  8.  Seismic velocities obtained after averaging the refractor velocity analysis function for various ranges of XY values which are symmetrical about 2.5 m.

DownLoad: Full-Size Img PowerPoint

Figure 9 shows the seismic velocities obtained with the application of equation (8) to the refractor velocity analysis function that has been averaged over the same range of XY values from -20 to 25 m. However, the window ∆x, over which the seismic velocity is computed, is not constant, but instead ranges from 5 to 50 m. Furthermore, the two times used in equation (8) are averaged over similar windows, with the center stations also being the same distance apart as the length of the window. A 30 m value has been selected as providing optimum smoothing without undue loss of resolution. For this value, the two times used in equation (8) have been averaged over seven stations or 30 m intervals, the center stations of which are 30 m apart.

Figure  9.  Seismic velocities obtained by increasing the length of the window over which the refractor velocities are computed from the refractor velocity analysis function averaged over the ranges -15 m ≤ XY (5 m) ≤ 20 m.

DownLoad: Full-Size Img PowerPoint

The unusually high values in the vicinity of station 49 correspond with low head wave amplitudes in Fig. 15, low amplitude products in Fig. 17, and in turn, high inverted amplitude products in Fig. 18. Furthermore, the 3D results in Fig. 22 suggest that the region with the high seismic velocity on cross-line 49 is offset from in-line 19, which is the location of the 2D traverse. Therefore, it can be concluded that the high seismic velocities in the vicinity of station 49 may be the result of out-of-plane side swipes and/or diffractions.

It is recognized that the extensive averaging in Figs. 8 and 9 can result in an effective reduction in resolution. However, it is considered that the need for an objective methodology is more important. Furthermore, it will be noted that the GRM WET tomogram in Fig. 13 has recovered much of the short spatial wavelength detail that is not present in the considerably smoothed starting model.

An alternative approach is the novel application of the Hilbert transform described in Chopra and Marfurt (2007). This algorithm computes an average of the seismic velocities derived over a range of distances ∆x, such as 5, 10, 15, 20, and 25, which are centered on the reference station. The results obtained with this range of distances ∆x, for the range -20 m ≤ XY (5 m) ≤ 25 m are comparable with those shown in Fig. 8 over a similar range of XY values. Figure 11, which is described below, presents the results computed with 5 m ≤ ∆x (5 m) ≤ 30 m, for a range of XY values from -2.5 to 7.5 m, using averaging windows of nine or ten XY values.

The development of simple reciprocal gradient operators is the subject of ongoing research. Experience with other sets of data has shown that the averaging achieved with that shown in Fig. 8 with equation (8) is generally adequate (e.g., Palmer, 2008a, Fig. 8). However, further averaging, such as that shown in Fig. 7 or the Hilbert transform, can be employed where the data indicate that it can be useful.

Determination of the optimum XY value

The averaging operation can be readily extended to other, non-optimum XY values with the same aim of preserving the shape of the graph, while simultaneously minimizing any noise. Figure 10 presents the seismic velocities in the refractor computed for a range of XY values from -2.5 to 7.5 m, using nine or ten XY values and an averaging window of 30 m with equation (8). The seismic velocities are smoothly varying and as a result, the determination of the most probable model of the seismic velocities and the optimum XY value(s) are more convenient than is possible with the graphs in Fig. 7.

Figure  10.  Seismic velocities computed over a range of windows of nine or ten XY values and centered on XY values for the range -2.5 m ≤ XY (2.5 m) ≤ 7.5 m.

DownLoad: Full-Size Img PowerPoint

Where plane layering is appropriate, equation (1) essentially determines the traveltime from the shot point to a position below the reference point, G in equation (1) (Palmer, 1980, equation (5); Palmer, 1986, equation (8.2)), as is shown in Fig. 4. In Fig. 7, the regions at either end of the traverse and between stations 50 and 60 show much the same values, indicating that plane layering is a reasonable approximation and that the seismic velocities can be determined reasonably accurately. However, that approximation is not valid elsewhere, such as in the vicinity of station 40 where there is a large change in depth as shown in Fig. 12, and as a result, there can be some uncertainty in the determination of the true seismic velocities.

Nevertheless, the total traveltime through the refractor in the vicinity of station 40 remains the same. Therefore, if a low seismic velocity is assigned to one narrow region, then a high value must be assigned to an adjacent region, in order to preserve the total traveltime through the combined regions. For example, if a low seismic velocity is assigned to the region between stations 32 and 40, then a high value must be assigned to the region between stations 40 and 49, as shown with XY=-2.5 m in Fig. 10.

The optimum XY value is selected where the seismic velocity model is the simplest, that is, the extreme lateral variations are minimized. Figure 10 indicates an optimum XY value of 2.5 m, and it shows that the seismic velocities range from 2 700 m/s to 6 200 m/s.

As described above, the Hilbert transform has also been trialed. Figure 11 presents the results computed with 5 m ≤ ∆x (5 m) ≤ 30 m, for an effective range of XY values from -2.5 to 7.5 m, using averaging windows of nine or ten XY values. While the determination of the optimum XY value requires more care than is the case with Fig. 10, nevertheless, a value is 2.5 m can still be recovered. Also, there is more variability in the seismic velocities, which addresses any fears for a loss in resolution, although there are minor concerns that the seismic velocities may be unrealistically high in the vicinity of stations 40 and 49.

Figure  11.  The seismic velocities computed over a range of windows of nine or ten XY values and centered on XY values for the range -2.5 m≤ XY(2.5 m)≤7.5 m, using the Hilbert transform with 5 m ≤ ∆x (5 m) ≤ 30 m.

DownLoad: Full-Size Img PowerPoint

Generating GRM tomogram

Figure 12 presents the GRM time model computed with equation (3) for a range of XY values, using the traveltimes from SP1 and SP97 that have been corrected for minor variations in the surface soil layers and picking errors. As stated previously, there is minimal variation in the time models for a wide range of XY values.

Figure  12.  The GRM time models for the range -10 m≤ XY (5 m) ≤ 15 m.

DownLoad: Full-Size Img PowerPoint

The time structure is converted to a depth model using equation (5). The depth cross-section was firstly gridded in SURFER with a kriging algorithm, and then used as a starting model for WET inversion with RAYFRACT, and is shown in Fig. 13. In this study, the default settings were used: no steps were taken to remove any obvious artifacts in the gridded starting model caused by either the gridding or the inversion algorithms. The aim was to assess whether tomography is able to recognize and accommodate artifacts of various types.

Figure  13.  The starting model and WET tomogram obtained with the GRM inversion algorithms. No account has been taken of the velocity reversal in the weathering.

DownLoad: Full-Size Img PowerPoint

The seismic velocities for the starting models for the two intervals between stations 51 and 61, and between stations 69 and 73, where the low seismic velocities occur, are much the same for 0 and 5 m XY values as is shown in Fig. 10. Accordingly, the tomogram shown in Fig. 13 is similar to that derived with a starting model obtained with 0 and 2.5 m as well as a 5 m XY value.

The GRM WET tomogram in Fig. 13 shows that two regions in the refractor with low seismic velocities can be recognized between stations 53 and 61, and between stations 68 and 73. These seismic velocities, which are less than 2 500 m/s and 40% lower than the lowest seismic velocity in the refractor indicated in the τ-p tomogram in Fig. 3, provide very strong evidence for one or more shear zones. Furthermore, the 3D cross-line velocities, shown in Fig. 22 below, support the existence of the major shear zone between stations 53 and 61. By contrast, the τ-p tomogram shows the reverse situation between stations 53 and 61, with an increase in the seismic velocity in the deepest layer together with a decrease in the depth to it.

Inversion with a Velocity Reversal

The likely occurrence of velocity reversals is indicated by the upward concave traveltime graphs in Fig. 2 and the head wave amplitudes in Fig. 16. It is supported by a comparison of the seismic velocities in the weathering of approximately 1 000 m/s computed from the direct arrivals, with the seismic velocities of ~450 m/s computed with the average vertical velocity of the GRM with equation (7), using an optimum XY value of 2.5 m.

Figure 14 shows a tomogram generated by replacing the direct traveltimes with values representative of an average velocity of 500 m/s, in order to accommodate a velocity reversal. This tomogram still supports the occurrence of a major shear zone between stations 51 and 61. Furthermore, it also indicates the possibility of another narrower low-velocity region in the vicinity of station 37. However, the isolated high-velocity fragment below station 47, which initially appears to be an artifact, is consistent with this segment being an out-of-plane side swipe. The orthogonal profiles in Fig. 22 support the likelihood of out-of-plane arrivals.

Figure  14.  The WET tomogram obtained by replacing the measured direct traveltimes with values representative of a velocity reversal with an average seismic velocity of 500 m/s in the weathering.

DownLoad: Full-Size Img PowerPoint

Figure 14 illustrates some of the benefits of explicitly addressing velocity reversals with lower average vertical velocities. The most obvious is that significantly different geological models, in this case the possibility of another shear zone, can result. However, perhaps a more important conclusion is that the use of unrealistically high seismic velocities with vertical gradients can result in a loss of resolution, and in some cases, even a failure to recognize out-of-plane events. Furthermore, it demonstrates that 2D refraction tomograms, irrespective of how closely they may fit the traveltime data, do not necessarily generate meaningful results where there are genuine 3D targets.

Velocity reversals are not well accommodated with the majority of refraction tomography programs. The replacement of the observed traveltimes in the overburden with values representative of a velocity reversal, in order to obtain the tomogram in Fig. 14, constitutes an elementary, but not necessarily an ideal solution. Although Fig. 14 generally illustrates the importance of determining an accurate parameterization of the overburden velocities in order to obtain representative depth estimates, detailed conclusions, such as the occurrence of an additional low-velocity zone in the refractor in the vicinity of station 37, can be less certain. Often, the 2D analysis of shot amplitudes can resolve these ambiguities.

Usefulness of Error Analyses

Table 1 shows that for comparable settings with RAYFRACT, in this case the default settings of 10 iterations, the GRM results provide more accurate inversion. However, the significance of these errors is questionable. Both the τ-p and the GRM tomograms may contain artifacts, such as the possible out-ofplane events, while the errors in the 1D τ-p tomogram can be greatly reduced by simply increasing the number of iterations, without any significant changes in the detail of the final tomogram.

A fundamental tenet of model-based inversion is that the response of the final model should be consistent with the original data. Even though it is of questionable validity, because no account has been taken of lateral variations in refractor velocities, or velocity reversals in the overburden layers, nevertheless, the 1D τ-p WET tomogram generated for the Mt. Bulga 2D case studies is consistent with the traveltime data. Furthermore, the tomogram that seeks to accommodate the likely velocity reversals has significantly greater errors than the alternatives.

These results demonstrate that simplistic assessments of errors do not necessarily indicate geologically appropriate tomograms, especially where there are velocity reversals. Instead, these results demonstrate that the selection of the inversion algorithm used to generate the starting model, together with appropriate model parameterization are equally important for effective refraction inversion.

In most cases, the final WET tomogram is quite similar to the starting model. In fact, it can be concluded that the starting models can be considerably more useful than the ray path coverage, because the latter depends largely on the former.

Head Wave Amplitudes

1D analysis of amplitudes

Head wave amplitudes constitute half of the volume of seismic refraction data (the other half being the traveltime data), which are too often overlooked. In many instances, there can be a reluctance to analyze head wave amplitudes because it is considered that there can be significant variations in the coupling of individual receivers, and therefore, significant variations in the outputs between individual receivers. However, Drijkoningen (2000) demonstrated that variations in coupling are generally minimal while Palmer (2006) demonstrated that any short wavelength variation in output is related to variation in the surface soil layers. Irrespective of the source of any variation in output, the application of the GRM smoothing routine (Palmer, 2006) largely removes their effect. The amplitude corrections are shown in Fig. 3.

Frequently, a 1D analysis of the gross shot amplitudes can often indicate whether uniform seismic velocities, vertical velocity gradients, or even velocity reversals are appropriate. Figure 15 presents the head wave amplitudes for the offset shot points located at stations 1 and 97. These amplitudes decrease with the source-to-receiver distance and therefore indicate that vertical velocity gradients are not significant in the main refractor. If there were vertical velocity gradients, then the head wave amplitudes would in fact increase with the source-to-receiver distance (Červený and Ravindra, 1971). The short wavelength departures from the geometric spreading, which in this figure are approximated with the reciprocal of the distance cubed, are attributed to local variations in the petrophysical properties of the refractor, and they are examined in a 2D analysis of the amplitudes below in Fig. 17.

Figure  15.  Uncorrected shot record amplitudes for SP1 and SP97 together with a reciprocal distance cubed geometrical spreading function.

DownLoad: Full-Size Img PowerPoint

Head wave amplitudes for selected near shot points located within the spread of receivers are presented in Fig. 16. They show initially over-scaling of the seismic recorder which is indicative of very high amplitudes, followed by very rapid attenuation, and then the establishment of a second mode of more "normal" amplitude decay with distance. This effect is usually recognized as cycle skipping, also known as shingling, on the shot records, and it is characteristic of velocity reversals (Whiteley and Greenhalgh, 1979; Domzalski, 1956).

Figure  16.  Uncorrected shot record amplitudes for selected shot points within the receiver spread.

DownLoad: Full-Size Img PowerPoint

This analysis of the shot record amplitudes does not support the default use of vertical velocity gradients shown in Fig. 3.

2D analysis of amplitudes

In the near field, where most surveys for geotechnical investigations take place, geometrical spreading dominates the observed head wave amplitudes. Palmer (2001a) demonstrated that a 2D analysis of shot amplitudes with the multiplication of the forward and reverse amplitudes effectively compensates for geometrical spreading and that the resulting amplitude products are essentially proportional to the square of the head coefficient, which is the refraction analogue of the Zoeppritz transmission coefficient in reflection seismology, as is shown in Fig. 17.

Figure  17.  The product of the forward and reverse corrected shot amplitudes.

DownLoad: Full-Size Img PowerPoint

Palmer (2001b) also demonstrated that the head coefficient is approximately proportional to the ratio of the specific acoustic impedance in the overburden to that in the refractor. Therefore, low seismic velocities in the refractor exhibit high amplitude products, while conversely, high seismic velocities in the refractor produce low amplitude products, provided of course that the seismic velocities in the overburden do not exhibit significant lateral variations.

For strong contrasts in seismic velocities, the head coefficient, K, is given by

where ρ and V_P are the density and compressional wave velocity in each layer.

Since the amplitude product, which has been corrected for near-surface irregularities with the GRM smoothing method (Palmer, 2006), is approximately proportional to the square of the head coefficient, it follows that

Accordingly, amplitude products inverted with equation (10) can provide an additional useful attribute which is an approximate measure of seismic velocities.

However, for the detailed attribute analysis of head wave amplitudes, it is necessary to more accurately compensate for the residual geometrical component, in addition to the use of the product of the forward and reverse values. Although the geometric spreading is widely assumed to vary with the distance squared, that function is only applicable after the signal has travelled several wavelengths and is in the far field. In the near field, the geometric component can be considerably larger. Since the great majority of seismic refraction surveys for geotechnical investigations take place in the near field, it is not uncommon to obtain geometric spreading functions that are quite different to the far field approximations, as is shown in Fig. 15.

Figure 18 presents the square root of the reciprocal of the scaled amplitudes products that have been corrected with the GRM smoothing method (see Fig. 6) and have had the geometric component applied using a range of exponents. As a result, the vertical scale correlates approximately with the seismic velocities in the refractor. In this study, an exponent of 6 has been judged to have effectively eliminated the residual geometric component.

Figure  18.  Inverted amplitude products corrected for geometric spreading with a range of exponents.

DownLoad: Full-Size Img PowerPoint

Figure 19 shows the seismic velocities computed at each station with equation (8) and with the modified Hilbert transform, and the reciprocal of the square root of amplitude product, which has been corrected for the residual geometric spreading, using a distance to the sixth power function. The correlation between the seismic velocities and the processed amplitudes is generally very good. Figure 19 does not support the occurrence of a major shear zone with a low seismic velocity zone in the refractor in the vicinity of station 37. It might be noted that the very large inverted amplitude at station 49 correlates with probable out-ofplane arrivals, as is suggested by Figs. 14 and 22.

Figure  19.  The seismic velocity determined for XY=2.5 m with a 30 m window and the Hilbert transform, and the inverted amplitude product shown in Fig. 18 for a reciprocal distance to the sixth power geometric correction function.

DownLoad: Full-Size Img PowerPoint

The determination of a representative exponent for the geometric spreading is necessary only with the analysis of single-fold data and single receiver spreads. With multi-fold data, the inversion routine automatically compensates for both variations in the source parameters and the geometric spreading (Palmer, 2009).

Refraction Convolution Section

The refraction convolution section (RCS) is generated by the convolution of pairs of individual forward and reverse traces. The convolution operation adds traveltimes and multiply amplitudes, and hence generates the same time structural model of the refractor as that obtained with the GRM algorithms. Furthermore, the multiplication of the amplitudes largely compensates for the geometrical spreading component, with the result that the RCS amplitudes are approximately proportional to the square of the head coefficient.

Figure 20 shows the RCS generated from the shot records for the offset shot points located at stations 1 and 97 for the Mt. Bulga profile, together with the GRM time model superimposed. Figure 20 demonstrates that the time model generated with the RCS is consistent with the time model generated with GRM, and that the amplitudes of the RCS first arrivals are consistent with the results shown in Fig. 17.

Figure  20.  The single-fold RCS generated with SP1 and SP97 and the GRM time model.

DownLoad: Full-Size Img PowerPoint

A major advantage of the RCS is that it facilitates the application of the suite of techniques generally known as attribute analysis (Chopra and Marfurt, 2007) to near-surface seismic refraction data. The spectral analysis of the RCS in Fig. 21 shows that the region in the sub-weathering between stations 51 and 61 with the low seismic velocity has significantly attenuated the high-frequency components. This reduction in the high-frequency response is consistent with increased fracturing and a reduction in mechanical rigidity in the shear zone. Furthermore, both the spectral analysis in Fig. 21 and the GRM WET tomogram, which seeks to accommodate the velocity reversal in the overburden in Fig. 14, suggest that the shear zone is wider and extends from station 49 to station 61.

Figure  21.  Spectral analysis of the RCS generated with SP1 and SP97.

DownLoad: Full-Size Img PowerPoint

There is also a smaller reduction in the highfrequency content in the vicinity of station 37. This region also corresponds with lower seismic velocities in Figs. 10 and 11, much lower seismic velocities in Fig. 14, and an increase in depths of weathering in Fig. 13. This region of the sub-weathering may represent low to moderate levels of fracturing but not necessarily to the extent characteristic of a shear zone. However, the combined attributes of increased depths of weathering, lower seismic velocities, and reduced high frequency response indicate that this region may be anomalous, and that it might possibly warrant further exploration, depending upon the objectives of any site investigation program. By contrast, there is no indication of any comparable variations in the subweathering in the τ-p tomogram in Fig. 3.

The major focus of this study has been on the detection and delineation of the major shear zone between stations 51 and 61. Such features are not common at most sites and more often it is necessary to define less extreme variations, such as those in the vicinity of station 37. This study demonstrates that the application of detailed attributes derived with the GRM and the RCS can significantly improve the resolution and reliability of routine geotechnical site characterization with near-surface seismic refraction methods.

Extensions to 3D

The 3D refraction survey was carried out some time after the 2D traverse described above (Palmer, 2001c). The 2D traverse corresponds with in-line 19, while the in-line station numbers are essentially the same for both surveys. Figure 22 shows the seismic velocities derived in the base of the weathering from the seven cross lines, and confirms the existence of the low-velocity shear zone along virtually all of crosslines 53, 57, and 61. Figure 22 adds two additional items of considerable geological significance.

Figure  22.  The seismic velocities computed in the sub-weathering with the GRM for the seven cross-lines 45, 49, 53, 57, 61, 65, and 69. Depths in meters are shown as white contours.

DownLoad: Full-Size Img PowerPoint

First, although the survey was oriented so that the in-lines would be orthogonal to the known shear zone and the dominant strike direction of the local geological structure, which was inferred to be N-S, the 3D results show that the true strike is closer to NW-SE. A common observation with numerous petroleum exploration and production case studies has been that 3D results can show a rotation of 45°, and even up to 90°, in the dominant strike direction of the faulting, as well as a change from a few faults with large strike lengths to many faults with smaller strike lengths, when compared with coincident 2D results (Ruijtenberg et al., 1992; Fig. 2).

Second, there are several significant ENEstriking cross-cutting structures, which are indicated by the lateral changes in the seismic velocities in the refractor. These ENE structures are not readily detectable on any single profile with any orientation, because they do not exhibit the same distinctive low seismic velocity signature as the major NW-SE shear zone. These ENE structures emphasize the fact that all geological features, whether nominally 2D such as at Mt. Bulga or elsewhere, are 3D in reality.

Figure 23 presents the cross-line corrected amplitude products processed to reflect seismic velocities with equation (10). This figure is a little more complex, probably because the seismic velocities in the weathered layer vary laterally. Nevertheless, the NW-SE strike of the shear zone can still be recognized, as can the ENE structures. Furthermore, the high amplitude effect previously recognized at station 49 of the 2D traverse in Fig. 19 can also be seen on crossline 49 near the intersection with in-line 19.

Figure  23.  The inverted amplitude products for the seven cross-lines 45, 49, 53, 57, 61, 65 and 69.

DownLoad: Full-Size Img PowerPoint

DISCUSSION

The Mt. Bulga and other case studies are compelling demonstrations of the fundamental nonuniqueness of all near-surface refraction inversion. Unless specific measures are taken to explicitly address non-uniqueness, the production of a single refraction tomogram that fits the traveltime data to sufficient accuracy does not necessarily demonstrate that the result is correct, or even the most probable.

The τ-p inversion algorithm is currently used as the default for generating starting models with most refraction tomography programs. However, the Mt. Bulga case study demonstrates that the τ-p algorithm is unable to detect even major low-velocity shear zones in excess of 10 stations in width, and that 2D rather than 1D inversion algorithms are required where the detailed lateral resolution of seismic velocities in refractors is a primary objective.

In those cases where the reality of nonuniqueness is recognized, the most common approach is to employ a priori information to limit the range of some model parameters. Normally, either borehole data or other geophysical datasets are considered as suitable a priori information. Unfortunately, such a priori data are commonly not available for geotechnical and environmental investigations in reconnaissance surveys.

However, head wave amplitudes, which constitute a major proportion of the data, provide an immediately accessible source of a priori data, which are too often overlooked. The 1D analysis of head wave amplitudes can usefully resolve many ambiguities with the parameterization of the inversion model in terms of constant velocities, vertical velocity gradients, or velocity reversals.

Most refraction tomography programs parameterize the model with vertical velocity gradients, which can result in a loss of resolution in the definition of irregular (2D) refractors as well as a failure to recognize and accommodate out-of-plane (3D) refractors. However, it is usually not possible to uniquely determine from the traveltimes from those layers whether constant velocities, vertical velocity gradients, or velocity reversals are appropriate. In this study, the combination of the head wave amplitudes and the GRM average vertical velocity indicate that velocity reversals are more appropriate at the Mt. Bulga site. Furthermore, the 2D analysis of head wave amplitudes can often confirm lateral variations of seismic velocities within individual refractors, such as the subweathering layer.

The major advantage of using head wave amplitudes is that the amplitudes and traveltimes are different components of the same seismic signal, and therefore relate to the same layers and interfaces, and the same petrophysical property of seismic velocity. Consistency between amplitudes and traveltimes constitutes internal consistency within the seismic data, and relates to the reliability of the inversion process.

There is a widespread expectation that refraction tomography can significantly improve the resolution of the final result, that is, refraction tomography can reveal features that are not apparent in either the traveltime data or the starting model. However, the close similarity between both starting models used in this study and the final WET tomograms demonstrates the importance of the inversion algorithm used to generate the starting model. The τ-p tomogram demonstrates that if the gross lateral variations in the seismic velocity in the main refractor are not included in the starting model, then it is unlikely that they will be generated by the tomographic inversion process. Furthermore, the spatial resolution achieved with the final GRM WET tomogram in Fig. 13 in the main refractor is consistent with the resolution in Fig. 8 for an XY value of 2.5 m, and it is significantly better than the heavily smoothed starting model. Therefore, it can be concluded that the gross or long wavelength spatial resolution of the final tomogram is usually comparable to that of the starting model, but that refraction tomography can improve the detailed or short wavelength resolution of the final tomogram.

The emphasis on lateral resolution can be readily extended to full trace processing with the RCS. A major advantage of the RCS is that it conveniently incorporates amplitudes into the time structure model of the refractor. Furthermore, the RCS is ideally suited to the application of the wide range of techniques for post-processing seismic reflection data, generally known as attribute analysis (Chopra and Marfurt, 2007). In fact, the methods for deriving detailed lateral variations in refractor velocities described in this study can be viewed as elementary applications of attribute analysis. Even standard operations, such as the spectral analyses of the RCS, can usefully contribute to the generation of more detailed geotechnical models.

The derivation of more detailed geotechnical models of the regolith with innovative methods of data analysis is more effective with data acquired with multi-fold roll-through operations rather than with static receiver spreads. At the present time, most near-surface seismic refraction surveys for environmental or geotechnical applications are carried out by in-house personnel, often with in-house equipment. Frequently, these operations are under-capitalized with inadequate channel counts, and inefficient source procedures (Palmer, 2008b). The data acquired with these operations are not well suited to detailed or innovative processing, such as the generation of stacked RCS, and commonly, automatic processing with refraction tomography that requires only minimal technical proficiency is employed. These operations represent an emphasis on minimal or low-cost data acquisition that produces low-resolution conventional results, rather than on specialist geotechnical services. The development of detailed and innovative methods of data processing and interpretation will require many geotechnical organizations to determine whether their core business is data acquisition and conventional processing or whether it is specialist geotechnical services.

CONCLUSIONS

This study provides a compelling demonstration that 2D inversion algorithms provide significantly greater lateral resolution than 1D algorithms. The shear zone at Mt. Bulga is a major geological feature. It is 50 m wide and the seismic velocity of 2 500 m/s provides a large contrast with the surrounding metasediments in which the seismic velocity is approximately 6 000 m/s. Nevertheless, the shear zone is neither resolved nor even detected with the τ-p inversion algorithm. By contrast, the 2D GRM inversion algorithm detected the major shear zone in the 2D profile, the existence of which is corroborated both by the 2D analysis of the head wave amplitudes and numerous orthogonal profiles recorded as part of a later 3D survey. Furthermore, the GRM analysis of the 3D set of data facilitated the recognition of several ENE cross-cutting structures that did not exhibit comparable "signatures" in the seismic velocities. It can be concluded that the τ-p inversion algorithm generates low-resolution tomograms that are unsuited to the great majority of geotechnical and environmental investigations.

The major conclusion which can be reached with this study is that the acquisition and processing of near-surface refraction seismology should aim to maximize the superior lateral resolution of the refraction method. It is acknowledged that all seismic refraction operations should aim to provide as accurate depth estimates as is practical. However, the ongoing increasing use of high-resolution 3D seismic reflection and airborne geophysical methods is driven as much by the generation of more detailed geological models facilitated by the greater spatial resolution, as by any improvement in depth inversion. Nevertheless, one somewhat paradoxical outcome of the focus on lateral resolution is the determination of average vertical velocities above irregular refractors with the GRM.

This study has demonstrated the accurate definition and orientation of major structures with lateral velocity signatures, such as the NW-SE shear zone, as well as structures without lateral velocity expressions, such as the ENE faults. Furthermore, the more accurate determination of seismic velocity can facilitate the computation of azimuthal anisotropy or rock fabric, and in turn, fracture porosity. Accordingly, it can also be concluded that 3D refraction methods can significantly improve geotechnical and environmental site characterization and that they should be adopted as a matter of some importance.