- Copyright: © 2005 This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.
Surface-wave dispersion inversion is growing in popularity for geotechnical applications, due to its noninvasive character, relative straightforward field procedures and interpretation, especially when the subsurface structure is locally assumed to be one-dimensional (1D). Here, laser-Doppler physical modeling of surface-wave propagation is used to address issues of surface-wave depth penetration, the presence of dipping layers, and the associated limitations and systematic errors propagated in conventional 1D surface-wave inversion. Flat-layered models show that, with an active source and linear spread, the maximum resolvable wavelength of the Rayleigh-wave fundamental mode is on the order of 40% of the spread length. Linearised inversions confirm the rule of thumb that the depth penetration is 20–25% of the spread length, and that correct a priori layer interface depths from refraction analysis allow more accurate results. However, even under optimal conditions, failing to account for a dominant higher mode at low frequency when a stiff shallow layer is present, causes an overestimate of deeper layer shear-wave velocity. Moreover, a layer dip of only a few degrees can significantly bias the surface-wave inversion. If the incorrect a priori information from a single-shot refraction analysis is incorporated in the inverse problem, estimated interface depth depends on the shot position and deeper layer shear-wave velocity is underestimated. Even if correct a priori constraints are used, an underestimate of half-space shear-wave velocity of up to 25% remains.
Surface-wave methods are increasingly used in civil engineering applications to evaluate soil shear modulus with depth (Rix et al., 2001; O'Neill et al., 2003). This approach has advantages over invasive subsurface measurements (Xia et al., 2002), because it can be easily implemented along linear sections to obtain a two dimensional shear-wave velocity profile of shallow layers (Park et al., 1999; Lin et al., 2004; Hayashi and Suzuki, 2004) and it can be used as a tool for imaging subsurface heterogeneity (Leparoux et al., 2000; Shtivelman, 2002). The number of surface-wave profiling applications is growing, but questions about experimental and theoretical limitations are still outstanding (O'Neill, 2004a). One of these is the assumption that the probed medium is horizontally layered (1D), imposed by the popular 1D inverse problem formulation for data interpretation (Sambridge and Mosegaard, 2002).
The 1D assumption can be ‘violated’ by lateral geological or topographic variations (Gabriels et al., 1987), almost ubiquitously encountered in field surveys. In surface-wave phase velocity measurements, which are standard for engineering surveys, it is lateral variations under the recording array which corrupt the wavefield and thus contaminate the inversion (in group velocity measurements, it is those along the source-receiver path). To reduce these effects, a shorter recording spread could be employed, but this limits the maximum resolvable wavelength, and thus the depth of investigation (Forbriger, 2003; O'Neill, 2004b). Moreover, even if the 1D assumption is acceptable, experimental dispersion observations almost systematically underestimate surface-wave velocity at low frequency, and resolution is poor. The spread length is proposed as a factor of these low-frequency discrepancies, traditionally called ‘near-field effects’, even more influential than the source spectrum or geophone response (O'Neill, 2003). Increasing the spread length can reduce these effects, but then lateral discontinuity effects will become more influential, voiding the 1D assumption during inversion. Herein lies the tradeoff between lateral and vertical resolution in 1D surface-wave inversion methods.
Kuo and Nafe (1962) studied Rayleigh-wave propagation through non-horizontally layered media and pointed out the fact that, at a given point on the surface, there is a local surface-wave phase velocity independent of the direction of wave propagation. Assuming a linear variation in layer thicknesses, the surface-wave phase velocity determined from a pair of receivers is then considered as an effective surface-wave phase velocity. As the distance between the receivers approaches zero, effective surface-wave phase velocity approaches the local surface-wave phase velocity. An interface slope for instance, can be thus regarded as several stepped flat interfaces (Kuo and Thompson, 1963). This allowed Hayashi and Suzuki (2001) to reconstruct two-dimensional velocity structures and to develop a Common-Mid-Point (CMP) cross-correlation analysis of multi-shot surface-wave data (Hayashi and Suzuki, 2004). Even when such processing is used, dispersion determined under the spread may not describe the corresponding average flat medium when the lateral variations are too high. It is therefore necessary to evaluate the reliability of this assumption (Gucunski et al., 1996; Bodet et al., 2004). However, only under controlled conditions can this reliability be assessed. As it is uncommon to have perfectly controlled field-scale test-sites at one's disposal, either numerical or physical modeling seems the logical step.
Numerical simulations by Bodet et al. 2004 employing two-dimensional (2D) finite-element methods showed that the slope influence on dispersion curves cannot be summed up to a typical averaging of the media at the base of the receivers. In addition, the spread length is even a more limiting factor in the case of a sloping interface, and thus affect the obtained dispersion curves. Physical modeling methods have also been applied to elastic surface-wave propagation studies. Invariably, these employed non-contacting ultrasonic techniques (Scales and van Wijk, 1999; Meier et al., 1997; Nishizawa et al., 1997; Hayashi and Nishizawa, 2001; Campman et al., 2003; Campman et al., 2004; O'Neill, 2004a). The only modeling work which has investigated plane dipping-layered models is that of Kuo and Thompson (1963), albeit with a contacting receiver transducer.
In this study, laser-Doppler vibrometer measurements provide controlled analogues in the laboratory of field-scale surface-wave dispersion analyses. Both horizontally layered models, as well as models with an interface sloping a few degrees from horizontal illustrate the tradeoff between surface wave penetration depth and limitations of the 1D inverse problem.
We used a laser-Doppler vibrometer to record ultrasonic particle motion excited by a contacting transducer at the surface of a polymethylmethacrylate (PMMA) layer, also known as Plexiglas or Acrylic, over an aluminum half space. This kind of model and experimental setup has also been used by previous authors (Scales and van Wijk, 1999; Campman et al., 2003; Campman et al., 2004). Below are descriptions of the model (cf. Fig. 1 and Table 1), equipment and acquisition parameters.
Anticipating a dominant Rayleigh wavelength of about 6 mm, three models were constructed: Models 1 and 2 comprised PMMA plates of constant thickness (6 mm and 3 mm) over the aluminum substrate, used to verify the depth penetration ability in ‘ideal’ (flat-layered) cases. Model 3 comprised a sloping PMMA plate on the order of two percent, with average thickness comparable to models 1 and 2, to test the limitations of the 1D assumption.
It should be noted that the acoustic couplant used between PMMA and the aluminum can be up to 1.5 mm thick. For model 1, two sections, both with a thick (1.5 mm) and with a very thin coupling layer, were constructed, in order to ascertain its influence on the results. For models 2 and 3, the couplant was used throughout and was about 1 mm thick. In some cases, the coupling layer contains air bubbles, which can be seen through the PMMA layer. These, as well as some impurities in the PMMA plates, are significant in size with respect to the wavelengths used. Acquisition lines over both ‘regular’ and ‘irregular’ areas were surveyed for comparison.
To take advantage of the flexibility of a completely non-contacting measurement, it is preferable to also use a high frequency laser source in addition to the receivers (Campman et al., 2003; Scales and Malcolm, 2003; O'Neill, 2004a). This type of source is also used to improve the accuracy of surface waves dispersion in materials, for mechanical engineering applications (Ruiz and Nagy, 2004; Wu and Lui, 1999). But, as the source signal is an elastic pulse generated by thermal expansion, it is not possible to use it on PMMA without damaging the material. Instead, we used a 6 mm diameter P-Wave Panametrics® (V133) piezoelectric transducer having a 2.25 MHz resonant frequency, driven by a pulse generator. It provided us measured signals with a dominant frequency between 150 and 250 kHz.
A Polytec® laser-Doppler vibrometer measures absolute particle velocity on the surface of the sample via the Doppler shift. The output of the vibrometer-head is a beam of diameter less than 1 mm and a wavelength of 633 nm (red). Once the beam reflects off a moving target, its frequency is Doppler shifted. The beat frequency of the output plus the reflected signal is decoded in the hardware to give an absolute measurement of particle velocity, without contacting the medium. In contrast, contacting transducers are big compared to the dominant wavelength, and become part of the model, acting as scatterers.
The signal of the vibrometer is amplified with a low-noise preamplifier (SR 560 with 12 db/octave 10 kHz high-pass filter) and digitized at 14-bit resolution using a Gage digital oscilloscope card, attached to a PC. However, to ensure high signal-to-noise, reflective tape is applied to the surface of the model for better reflection of the vibrometer beam. While stacking (averaging) of multiple shots improves data quality, higher levels of noise are mostly caused by a decrease in reflection strength of the vibrometer beam. In practice, this often means that there is either an air-bubble between tape and the model, or the reflection strength of the tape is diminished by dirt. Such perturbations, as well as the heterogeneities of the coupling layer or the PMMA plate surface irregularities, can be regarded as noise, comparable to bad geophone coupling, local heterogeneity, topography and interface or surface roughness effects in the field-scale seismic data acquisition conditions.
Computer controlled scanning with the vibrometer allows us to obtain a high density data set, so that we can easily create a scaled version of classical field acquisition parameters. We typically recorded 101 seismic traces with 0.4 mm spacing in linear, single-channel walkaway fashion, over offsets from 10 mm to 50 mm. The time sample rate is 10 MHz and each trace is stacked 250 times and consists of 4,096 samples.
Refraction and Surface-Wave Dispersion Analysis
Model 1: First Layer Thickness on the Order of Dominant Rayleigh Wavelength
Two shot records for model 1 are shown in Fig. 2. With the shot position fixed, the two records (‘line 1’ and ‘line 2’) form a split spread acquisition. Despite the large absorption of the PMMA, we are able to record body and surface waves with signal-to-noise ratio greater than 25 dB up to 50 mm from the source. Coherent waves appear on this seismogram, such as direct body waves (A on Fig. 2) and multiples reflected from the bottom of the model, occurring at 77.5 microseconds (D on Fig. 2). These bottom-reflected waves arrive after the surface-wave for offsets greater than 50 mm, so we can compute Rayleigh-wave phase velocities without interference from these boundary effects. However, perturbations arise from defects on the reflective tape (cf. F on Fig. 2) or surficial undulations (cf. E on Fig. 2).
For each experiment, we first performed a refraction study since it can aid the inversion process (Jongmans et al., 1996). A basic analysis presented in Fig. 3, shows that the refracting surface coincides with the bottom of the coupling layer. Thus, we hereafter consider the PMMA-coupling combination as a homogenous layer for the surface-wave tests.
Experimental dispersion images are depicted in Fig. 4 and are overlain with the theoretical dispersion curves calculated from the actual physical parameters (Herrmann, 2002). The experimental dispersion images are determined using the slant-stack method in common shot gathers, followed by a 1D Fourier transform over the intercept time to obtain the wavefield in the phase velocity versus frequency domain (McMechan and Yedlin, 1981). In this plane the dispersion curve of a particular mode corresponds to the maximum of the wavefield (in black on Fig. 4), which is extracted with an estimate of standard error in phase velocity (Mokhtar et al., 1988; Herrmann, 2002). Higher modes may be present in the experimental images, but are not coherent enough to be picked.
The maximum wavelength observed in Fig. 4 is about 20 mm, which correlates with the long-wavelength cutoff criteria recommended by O'Neill (2004b), which specifies that using an active source and linear array, the maximum measurable wavelength from plane-wave transform is limited to 0.4 times the spread length. With this rule, the maximum resolvable wavelength for the 40 mm spread used should be 16 mm, and if the maximum depth of investigation is assumed as half of this wavelength (8 mm), then a 7.5 mm thick first layer should not be adequately detected. In Line 1 of Fig. 4, this effect can be seen as the distinct loss of resolution and lack of high phase velocities at low frequency.
Model 2: First Layer Thickness on the Order of Half the Dominant Rayleigh Wavelength
We performed a second experiment, using model 2 (a thinner, 3 mm PMMA layer) and recorded two split-spread lines, i.e., ‘line 3’ and ‘line 4’ of Fig. 5. Again, velocities determined from the refraction study correspond with the known material properties, taking the 1 mm coupling layer into account (cf. Fig. 6). Several multiples occur on the Line 4 data (E on Fig. 5, Line 4), which could be from the more than 3 mm diameter air bubbles in the coupling layer observed along this line. This feature appears to be a cause of discrepancies between the Line 3 and Line 4 dispersion images (cf. Fig. 7).
The low-frequency part of these images shows maxima up to the aluminum Rayleigh-wave phase velocity. The maximum observed wavelength is about 30 mm and the theoretical dispersion curves fit the data within estimated errors less than 10% only for wavelengths less than 20 mm. But the Rayleigh-wave dominant wavelength is almost twice that of the first layer thickness so the low-frequency portion of the dispersion curve is apparently responding to the aluminum substrate.
However, at very low frequency, a higher mode almost overlaps the fundamental on the Line 4 dispersion image (Fig. 7). It was also noted, but not as clearly, in the Line 2 dispersion (6 mm PMMA, Fig. 4). This is an interesting low-frequency complication that will have an influence on the inversion. This kind of dominant first higher mode at very low frequency was also identified in field data by Socco and Strobbia (2004). In sites with large elastic contrast at shallow depth, it was called a ‘pitfall’ for the inversion process if not taken into account.
Regarding the frequency range of our signals and our spread length, these two experiments allowed us to show that our investigation depth is unlikely to be greater than 8 mm, i.e., approximately half the maximum resolvable wavelength of 16 mm. This is similar to classical field test data (O'Neill, 2004b), where the rule of the thumb is a maximum depth penetration of about 20–25% the spread length in active source, MASW and linearised inversion of the fundamental mode.
Model 3: Dipping Interface
To evaluate the influence of the ‘spread length limit’ in the presence of a sloping interface, two lines (‘line 7’ and ‘line 8’) were recorded on a PMMA plate with substrate slope of almost 2.4%, i.e., 1.4 degrees. Shot gathers from shooting both ‘up-dip’ (Line 7) and ‘down-dip’ (Line 8) are shown in Fig. 8. The refraction analysis gave a velocity of 2,500 m/s for direct P-wave (A on Fig. 8). However, because of the slope, we obtain ‘apparent’ velocities for the refracted P-wave of 5,600 m/s for Line 8 and 6,200 m/s for Line 7 (B on Fig. 8), as well as apparent thicknesses around 5.5 mm (cf. Fig. 9).
The dispersion images (Fig. 10) correlate moderately well with the theoretical curves for the corresponding ‘flat average medium’ (or ‘central shear-wave velocity (VS) profile’, 6.15 mm in the down-dip case and 4.65 mm in the up-dip case). However, at lower frequencies, where wavelength reaches the spread length limitation (40%), the modal maxima are shifted in frequency and errors exceed 10%. Obviously, assuming a 1D subsurface structure will lead to the wrong a priori information and thus a systematic error in the inversion results.
Limitations of Fundamental Mode Inversion
The experimental dispersion curves were inverted to estimate a vertical shear velocity structure at the base of each line. An iterative linear inversion was performed, using the method of Herrmann (2002). It is important to note the large non-uniqueness of the problem and that this type of iterative linear inversion merely provides an estimated model (Mokhtar et al., 1988; Xia et al., 1999).
For each line, two different initial layered models were proposed, characterized by thicknesses, P-wave and S-wave velocities as well as densities of each layer. During the iterative process, Poisson's ratio was fixed and only the S-wave velocities were allowed to vary. Of the two different initial layer parameterizations, both ‘multi-’ and ‘simple-’ layering was used. This is similar to real field surveys, where often ‘blind’ and ‘borehole’ inversions are performed. ‘Blind’ inversions use a stack of many thin layers, since ‘true’ interfaces are unknown. ‘Borehole’ inversions employ a priori depth information, usually geological or physical contrasts identified from nearby boreholes, or from other geophysical tests. Here, refraction analysis provides the a priori constraints, namely half-space depth and P-wave velocities (for Poisson's ratio).
Flat-layers (Models 1 and 2)
Without any a priori information and using estimates of the medium S-wave velocities from the experimental dispersion curve, a ‘blind’ inversion for Line 1 (7.5 mm first layer thickness, case A in Fig. 11) yields the following results: a first layer can be guessed down to 5 mm, with a shear-wave velocity between 1,000 m/s and 1,400 m/s which increases to 2,600 m/s at about 7.5 mm. This final model thus shows an underestimate of more than 20% of the true half-space shear-wave velocity and it is impossible to determine an interface. However, if we include the information provided by the refraction study (i.e., P-wave velocity, thickness and empirical estimates of S-wave velocities and densities), it is possible to determine the half-space shear-wave velocity more accurately, i.e., within errors lower than 10% (cf. ‘borehole’ inversion, case B in Fig. 11).
When the first layer thickness is 6 mm, ‘blind’ inversion results (case A in Fig. 12) resemble the Line 1, case A. However a half-space appears at about 6.5 mm with a shear-wave velocity of 2,900 m/s (10% of the true half-space shear-wave velocity). ‘Borehole’ inversion presented in Fig. 12, case B, highly improves the layers velocity and thickness determination. Then, unlike the 7.5 mm first layer, the 6 mm layer is detected, with both multi-layer and simple-layer models. This finally ascertains the investigation depth of 8 mm previously inferred from the dispersion analysis.
Figure 13 presents possible models inverted from Line 3 (4 mm first layer thickness) dispersion and refraction data. A first attempt (Fig. 13, case A) shows ‘blind’ inversion results. The final model resembles the actual velocities, but it is impossible to determine the interface depth precisely. If ‘borehole’ inversion is performed, a good evaluation of the first layer thickness is possible and S-waves velocities are well determined (Fig. 13, case B). However, a progressive transition between the first layer and the half-space appears at about 4 mm and 2,600 m/s, which is typical of the poor resolution of the method at low frequency. This ‘smearing’ at depth was also noted by O'Neill (2004b), and highlights the importance of the spread length as a parameter of great influence in the dispersion analysis of surface waves.
Inversion results for Line 4 (again 4 mm first layer thickness) are presented as one possible solution in Fig. 14. The presence of heterogeneity (air bubbles) along Line 4 data reduces the dispersion image quality (cf. E on Fig. 5) which has a strong influence on the inverse problem solution. As previously noticed, the second mode predominates at lower frequencies up to 150 kHz (cf. Fig. 6) and may be difficult to discriminate from the fundamental. The inverse problem solutions consequently overestimates the half-space shear-wave velocity. It is representative of the misidentified mode pitfall in such inversion methods (Socco and Strobbia, 2004).
Dipping Interface (Model 3)
‘Blind’ inversion of the dispersion from both down-dip (Fig. 15, case A) and up-dip (Fig. 16, case A) shots shows a similar ‘smearing’ as in the flat-layered cases. For the down-dip shot the underestimate is over 25%, slightly more than in the ‘corresponding’ flat-layered case of Line 2 (Fig. 12, case A). For the up-dip shot, although the 40 mm spread length provides an investigation depth of 8 mm and the average layer thickness is only 4.65 mm, the dispersion discrepancy at low frequency propagates as a pronounced underestimate of the half-space shear-wave velocity, compared to the ‘blind’ inversion of the 4 mm thick layer in Line 3 (Fig. 13, case A). These results compare well with the numerical study of Bodet et al. 2004, where the low-frequency dispersion discrepancy, due to the spread length limitation, is exacerbated by a sloping interface.
To properly characterize a slope with refraction seismic techniques, it is necessary to perform both forward and reverse (or reciprocal) shots. In the case of surface-wave profiling techniques, data is often only recorded ‘one way’, i.e., without a reverse ‘check’ shot. If so, the refraction results (which can be used for a priori surface-wave inversion constraints) are also biased if lateral variations are present. The ‘a priori’ inversion of Case B on Figs. 15 and 16 show that, even with this model incompatibility, the non-uniqueness of the inverse problem means that the surface-wave data can still fit an incorrect model within errors lower than 10%. Compared to the ‘corresponding flat-average’, however, the interface is modeled as too shallow when shooting down-dip (i.e., receivers are down-dip from the shot) and too deep when shooting up-dip (i.e., receivers are up-dip of the shot).
Finally, even if the depth to the half-space is set at the average value over the line (the ‘central shear-wave velocity profile’), the actual aluminum shear-wave velocity is underestimated (Fig. 17 and Fig. 18). Similar to the inversion of the up-dip shot dispersion (both ‘blind’ and ‘a priori’ in Fig. 16), the average depth to the interface is apparently overestimated (case C, ‘borehole’, in Fig. 18). For the down-dip dispersion, there again exists an underestimate of approximately 25% of the true half-space shear wave velocity.
Laser-Doppler physical modeling of surface-wave propagation proved to be an efficient tool to study the effects of depth penetration and dipping layers, and the associated limitations and systematic errors propagated in conventional 1D surface-wave inversion. Flat-layered models show that, with an active source and linear spread, the maximum resolvable wavelength of the fundamental mode is in the order of 40% of the spread length. In the case of a dipping layer, the observed dispersion correlates well with that calculated from the flat average over the spread length, except for wavelengths which reach the spread length resolution limitation.
These resolution limitations and model incompatibilities propagate into the uncertainty of the estimated models. In the flat-layered cases, linearised inversion of the observed dispersion curves using ‘blind’ initial layered models all show a ‘smearing’ at depth and underestimate of the half-space shear-wave velocity. Thus, a rule of thumb that depth penetration is 20–25% of the spread length is justified. Models with correct a priori constraints from refraction analysis are less affected. However, these fail to account for a dominant higher mode at low frequency when a stiff shallow layer is present, causing an overestimate of shear-wave velocity.
When a dipping layer is present, ‘blind’ inversion results show a similar ‘smearing’ at depth and an underestimate of half-space shear-wave velocity over 25% from dispersion of down-dip shots. Even when shooting up-dip, if the average depth to the half-space is within the depth penetration limit, deep structure stiffness is underestimated. If the incorrect a priori constraints from a single-shot refraction analysis are used, the interface is modeled as too shallow in the down-dip shot and too deep in the up-dip shot, compared to the ‘flat-average’, along with an underestimate of half-space shear-wave velocity of up to 25%, the dispersion still fitting the data to acceptable limits (i.e., within errors lower than 10%). Even if the equivalent flat-average depth to the half-space is a priori incorporated, an underestimate of half-space shear-wave velocity of up to 25% manifests. Without other a priori information, an interface sloping only a few degrees can bias the surface-wave inverse problem significantly. This reinforces the need for correct a priori information at the site, or alternative methods such as cross-correlation analysis of multi-shot surface waves, when the 1D assumption is ‘violated’, even slightly.
L. Bodet thanks J. Scales and the group of the Physical Acoustics Laboratory at Colorado School of Mines for technical assistance and discussion in the laboratory. We thank Xander H. Campman for assistance, critical remarks and suggestions of improvement. L. Bodet acknowledges his thesis grant from the Laboratoire Central des Ponts et Chaussées and the Bureau de Recherches Géologiques et Minières, and thanks Céline Tirel for useful comments en technical assistance. We worked on collaboration with the Institut de Physique du Globe de Paris thanks to an agreement with the Laboratoire Central des Ponts et Chaussées. Constructive reviews by Adam O'Neil, associate editor, Axel Tillmann and an anonymous reviewer helped improve the final version of the manuscript.