Evaluation of GPR techniques for civil-engineering applications: study on a test site

Journal of Applied Geophysics 45 Ž2000. 141–156 www.elsevier.nlrlocaterjappgeo

Evaluation of GPR techniques for civil-engineering applications: study on a test site G. Grandjean ) , J.C. Gourry, A. Bitri BRGM, Direction de la Recherche, Departement Geophysique et Imagerie Geologique, AÕe. Claude Guillemin, BP 6009, F-45060 Orleans Cedex 2, France Received 29 October 1998; accepted 12 July 2000

Abstract Different ground-penetrating radar ŽGPR. techniques have been tested on the same site in order to establish the performance and reliability of this method when applied to civil-engineering problems. The Laboratoire Central des Ponts et Chaussees ´ ŽLCPC. test site at Nantes, France, was selected because it includes most of the underground heterogeneities commonly found in urban contexts, such as pipes, small voids, etc. The GPR survey consisted in recording measurements in tomographic Žsurface to horizontal borehole measurements., monostatic Ž2D surface profiling. and bistatic ŽCommon Mid Point wCMPx. modes above various buried heterogeneities. Different processing techniques were also performed, such as tomographic inversion, 2D and 3D migration, velocity analysis, as well as numerical simulations, the results of which can be summarized in three points. Ž1. Although the different filling materials of the site can be distinguished by velocity and attenuation tomography, the buried heterogeneities are more difficult to identify because of limited resolution related to angular aperture and Fresnel zone. Ž2. 2D surface profiling can detect the different shallow heterogeneities, such as pipes and voids, down to a depth of several meters. Additional processing, such as forward modeling and attenuation curve analysis, provides more quantitative information related to the medium. A comparison between 2D and 3D migrated data highlights the error introduced when the structures are considered to be perfectly cylindrical. Ž3. CMP analysis gives relatively good estimations of vertical velocity contrasts when the medium is layered. A lithologic log can be derived assuming that the velocity changes are related to material variations. q 2000 Elsevier Science B.V. All rights reserved. Keywords: Geotechnics; Test site; Ground-penetrating radar; Velocity and attenuation tomography; Modeling

1. Introduction The detection of underground heterogeneities using non-destructive methods is a crucial problem in

) Corresponding author. Tel.: q33-2-38-64-34-75; fax: q33-238-64-33-61. E-mail address: [email protected] ŽG. Grandjean..

urban environments, especially for trenchless works. For example, unknown concrete house foundations or sandstone blocks located along the path of a drilling machine can cause major equipment damage and thus a loss of time and money. Geophysical methods, particularly the ground-penetrating radar ŽGPR., can detect such superficial bodies with a relative efficiency depending on the field context, the dielectric properties of the host material and the

0926-9851r00r$ - see front matter q 2000 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 6 - 9 8 5 1 Ž 0 0 . 0 0 0 2 1 - 5

142

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

nature and size of the bodies. The objective of this work, funded by the BRGM and the National Project on Trenchless Works, is to review the capabilities of GPR in urban environments for civil-engineering applications. We approached this problem by comparing the results obtained at the same site using different survey configurations and different processing techniques. The success of such work is mainly conditioned by the features of the site where data acquisition is carried out. Firstly, the site must be representative of the urban environment, which means that the properties of the host materials and buried heterogeneities must be consistent with those commonly found beneath cities. Furthermore, the selected site must be well known and calibrated so as to be able validate the tested methods, i.e. the nature of the host material and the position of the buried heterogeneities must be precisely known. Finally, to ensure optimum data acquisition, the site must be easily accessible and as free as possible from noise sources Že.g. electrical installations., trees, etc. By respecting these criteria, a good compromise should be reached between reality and an idealized underground model. The geotechnical test site of the Laboratoire Central des Ponts et Chaussees ´ ŽLCPC. at Nantes, France ŽChazelas et al., 1997. was selected because it satisfies all these requirements. Numerous studies describe efficient GPR techniques for detecting and imaging underground pipes, voids, etc. ŽZeng and McMechan, 1997; Powers and Olhoeft, 1996; Tong, 1993; Annan et al., 1990.. However, because each technique is generally considered individually in a specific context, it is difficult to compare reliability when applied together under the same field conditions. In view of this, we present here the results of three different experiments systematically tested and compared at the same test site. The first experiment consisted in tomographic measurements from the surface to a horizontal borehole, and was dedicated to estimating the velocity and attenuation fields across the site. The second was a series of monostatic 2D surface profiling above each known buried heterogeneity and recorded using different antenna frequencies, complemented by 3D coverage over a specific area. Finally, bistatic Common Mid Point ŽCMP. measurements in the subhorizontally layered part of the site constituted the

third experiment. To guarantee the quality of interpretations, we tested different processing techniques, such as velocity and attenuation inversion, 2Dr3D migration plus forward modeling, and velocity analysis on the CMPs, which led to the definition of the most appropriate technique for imaging each object and characterizing the dielectric behavior of each host material. After a description of the LCPC test site, an inventory of the different survey configurations used and the associated results are presented, followed by a discussion concerning the contribution of each acquisition and processing technique to imaging underground heterogeneities.

2. Field measurements and processing 2.1. The LCPC test site In 1996, LCPC built a test site for geophysical measurements ŽChazelas et al., 1997. composed of a pit of dimensions 26 = 20 = 4 m, divided into five compartments filled with different host materials. Fig. 1 is a sketch of the site showing its main characteristics. Compartments 1, 3, 4 and 5 are respectively filled with silt, limestone sand, crushed gneiss with a grain size from 14 to 20 mm, and crushed gneiss with a grain size from 0.1 to 20 mm. Compartment 2 is horizontally layered with these four materials. Several types of object are buried in the different compartments: polystyrene steps Žcompartment 1., iron pipes and PVC pipes filled with water or air Žcompartments 1, 3, 4, and 5., stones and large limestone blocks Ž3, 4, and 5. and a masonry and an iron girder Ž5.. Depending on the type of host material Žsilt, limestone, gneiss., the grain size distribution and the buried object Žpipe, void, block, stone, etc.., the GPR signals acquired will vary because they are affected by lithology, granulometry and the presence of diffractors. In addition, the nature and dimension of host materials and buried objects make this site representative of true civil-engineering contexts. Furthermore, as LCPC positioned with high precision all the objects and the limits of each compartment during the construction of the site, geophysical anomalies can be correlated to them without doubt.

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

143

Fig. 1. Schematic diagram of the LCPC test site Žfrom Chazelas et al., 1997.. The five compartments, labeled 1 to 5, contain different host materials Žsee legend. and different buried objects, such as pipes ŽE1 to E9., voids ŽD., blocks ŽB., masonry ŽM. and stones ŽS..

2.2. Surface–borehole tomography As described by Leggett et al. Ž1993. for seismic imaging, and by Olsson et al. Ž1992., Valle and Zanzi Ž1997., and Vasco et al. Ž1997. for GPR,

tomographic methods consist of inverting the observed radar wave travel times or amplitudes to determine the spatial distribution of velocity or attenuation fields. Observed travel times and amplitudes are those measured after the radar wave has traveled

Fig. 2. Schematic cross-section of the LCPC test site corresponding to the tomographic plane Ža., source–receiver geometry. The straight lines connect source–receiver pairs used in the crosswell experiment Žb..

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

144

through the medium to be studied, from the transmitter to the receiver antenna. In our configuration, a 100-MHz transmitter was displaced at the surface of the site according to the horizontal borehole direction, while the 100 MHz receiver was displaced inside the borehole ŽFig. 1.. Both antenna positions were spaced 0.5 m apart, producing 490 transmitter–receiver pairs for the experiment. The tomographic plane dimensions reach 27 m long by 4.5 m height, crossing perpendicularly the different compartments ŽFig. 2.. The deformation of the pulse, i.e. time delay and amplitude decreasing, depends on the velocity and attenuation of the medium, and is a function of the path adopted by the radar wave between the transmitter and the receiver antennas. By picking the time of the earlier signal and its corresponding amplitude, one can relate each antenna position to these travel time and amplitude values. These data are then inverted according to the technique proposed by Jackson and Tweeton Ž1994., where an iterative scheme is used to recover the velocity and attenuation fields. Observed travel time t results from the sum along the ray path l of the product between the elementary portion of ray d l and the slowness p, taken as the reverse of the local velocity value. t s pŽ l . dl

Ž 1.

Hl

For a large number of observations, a matrix notation can be used to express the slowness field P

Ž y20 log A q 20 log

P s Dy1 T

As Eq. Ž4. is linear with regard to a and r, the left-hand term can be estimated for each transmis-

Ž 2.

Eq. Ž2. is resolved using an Algebraic Reconstruction Technique ŽART. to reduce computation time and improve stability ŽMason, 1981.. As mentioned by Hollender Ž1999., the interaction between the antenna and the medium is a crucial point in GPR attenuation tomography because the interactions between the antenna and the soil affect the transmitted signal in amplitude and frequency. In our approach, where only amplitude effects are studied, we will correct the observed signals for antenna radiation pattern and coupling effects. The radiation pattern correction is estimated according to the work of Arcone Ž1995; Appendix.. The antenna coupling correction is introduced in the reference amplitude A 0 as described in Olsson et al. Ž1992.. The measured amplitudes A of the transmitted signal are then inverted assuming the following relation: A s A0 D Ž u T ,f T . D Ž u R ,f R .

s a l.

l

Ž 3.

Ž 4.

sion measurement so that:

D Ž u T , f T . D Ž u E , f E . y 20 log l q 20 log A 0 . 20 log e

eya l

Where DŽ u , f . is the radiation pattern correction for a ray with azimuth f and elevation u from the vertical, with Ž f T , u T . and Ž f R , u R . being azimuths and elevations of the ray at the transmitter and the receiver, respectively. Parameter a refers to the attenuation factor and l is the travel path of the radar wave. To calculate the a parameter, Eq. Ž3. is linearized using a decimal logarithmic conversion:

D Ž u T , f T . D Ž u E , f E . 20 log l q 20 log A 0 . 20 log e

Ž y20 log A q 20 log

vs. the travel times T and the partial derivative D matrices.

s a Ž l . d l.

Hl

Ž 5.

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

Eq. Ž5. now has the same form as Eq. Ž1. and can therefore be resolved using the algorithm already described for travel time inversion. The parametrization used to discretize the velocity field consisted in 108 by 16 squared cells. Computations were per-

145

formed using MIGRATOM software ŽJackson and Tweeton, 1994. based on an ART scheme and using either straight or curved ray geometries. According to the proposition of Ivansson Ž1987. for avoiding ray bending complications, the path of the radar

Fig. 3. Schematic cross-section of the LCPC test site corresponding to the tomographic plane Ža., velocity tomogram using straight rays Žb., velocity tomogram using curved rays Žc., and attenuation tomogram Žd.. The black circles indicate the antenna positions. See Fig. 1 for legend.

146

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

wave was first approximated to straight rays ŽFig. 3b. and convergence was reached after seven iterations. We then used curved rays ŽFig. 3c. for the following three iterations to see whether convergence improved. Two more iterations were run using curved rays before the model diverged. Fig. 3c and d shows respectively the velocity and attenuation distributions on a longitudinal cross-section of the site. The different compartments can be distinguished on the velocity and attenuation tomograms. The averaged values of inverted velocity for compartments 1, 3, 4 and 5 are respectively estimated at 0.09, 0.12, 0.17 and 0.12 mrns. Similarly, the averaged attenuation values are estimated respectively at 9.5, 4.3, 2.6 and 6.9 dBrm ŽTable 1.. Compartments 3 and 5 are characterized by mean velocity Ž0.12 mrns. and attenuation Ž4.3–6.9 dBrm. values compared to compartments 1 and 4 that are filled with highly contrasted materials: high attenuation Ž9.5 dBrm. and low velocity Ž0.09 mrns. for compartment 1, and low attenuation Ž2.6 dBrm. and high velocity Ž0.17 mrns. for compartment 4. This indicates that various kinds of material can be identified from velocity and attenuation tomograms, provided that their nature and granulometry show a sufficient contrast.

Concerning the buried objects, no discernible signal can be related to a specific target. Although some visible anomalies can be related to certain objects ŽD, B1, B2, M., this is not the case for other signals Že.g. a1, a2.. This difficulty to relate velocity and attenuation anomalies to buried objects can have two origins. The first one involves artifacts generated by errors when determining time zero, picking events or computing ray path geometry ŽHollender, 1999.. The other is related to resolution limits and uniqueness in tomography. This problem, examined in Williamson Ž1991., Williamson and Worthington Ž1993. and Rector and Washbourne Ž1994., can be due to projection truncation, limited angular aperture or Fresnel zone dimensions. According to Rector and Washbourne Ž1994., projection truncation can explain the resolution drop in the left and right extremities of the tomograms. Elsewhere, a limitation of the angular aperture defined by the geometry of the tomographic device can alter the spatial resolution. These authors state that the aperture-related resolution is approximately equal to the Fresnel zone when:

lrR s 58.8 tan Ž DF .

y2

where l is the wavelength, R is the raypath length and DF is the angular aperture. In our case, where l

Table 1 Averaged wave parameters Ž V: velocity; a : attenuation., dielectric parameters Ž K : relative permittivity; Q: quality factor. and GPR performance parameters Ž P: wave penetration; r: wave resolution. estimated from the different GPR techniques A, B and C carried out in compartments 1 to 5. ŽA. Velocity and attenuation values averaged for each node of the tomographic planes belonging to a specific compartment. ŽB. Wave parameters averaged from the 12 forward models for the 900, 500 and 300 MHz frequencies. ŽC. Velocity values averaged from the 13 NMO analyses. For parameters strongly dependent on frequency, extreme values for 300 and 900 MHz frequencies are indicated Compartment

1: Silt

3: Limestone sand

4: Gneiss 14r20

5: Gneiss 0r20

(A) Surface borehole tomography V Žmrns. 0.09 a ŽdBrm. 9.5

0.12 4.3

0.17 2.6

0.12 6.9

(B) 2D profiling K Q V Žmrns. a 300 – 900 MHz ŽdBrm. P300 – 900 MHz Žm. r 300 – 900 MHz Žm .

13 7 0.07 15–45 1.5–1 0.15–0.03

6 20 0.12 6–20 4.5–2 0.14–0.03

3 30 0.17 1.5–4.5 4.5–4.5 0.23–0.06

5.5 7 0.13 9–27 2.5–1.5 0.16–0.04

(C) CMP analysis V Žmrns.

–

0.10

0.15

0.10

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

is comprised between 0.9 and 1.7 m, the Fresnel zone can be considered as the main limiting factor in terms of spatial resolution.

rand, 1999. according to the processing flow sequence:

v

2.3. 2D surface profiling This field experiment consisted in 12 surface GPR profiles located above the buried objects. They were devoted to analysing the diffracted signal in the classical minimum-offset configuration with 300, 500 and 900 MHz centered frequency antennas. The profiles were processed using the software Radar Unix developed at the BRGM ŽGrandjean and Du-

147

v

v

v

spatial re-sampling of scans along the antenna displacement axis to correct for velocity variations of the antenna during displacement; recovering amplitudes vs. time with an adaptive gain function to correct for geometrical spreading and attenuation; coherence and frequency filtering according to the observed amplitude spectrum to improve the signal to noise ratio; static corrections to remove topographic effects;

Fig. 4. Schematic plan showing the location of profile 6 in compartment 3 of the LCPC test site Ža., and radar sections along profile 6 at 300 MHz Žb., 500 MHz Žc., and 900 MHz Žd.. P s pit limit; D s void; E1–E9 s three sets of pipes ŽE1–E4–E7, E2–E5–E8, and E3–E6–E9. made respectively of iron, PVC filled with water, and PVC filled with air.

148 v

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

migration to correctly position the objects, the migration velocities being calculated from the diffraction hyperbolas curvature.

As an example, Fig. 4 shows three typical processed sections before migration corresponding to profile 6 ŽFig. 4a., and recorded respectively with 300 ŽFig. 4b., 500 ŽFig. 4c., and 900 MHz ŽFig. 4d. antennas. These sections show reflected and diffracted signals in response to the different objects located in compartment 3: the lateral limit of the pit ŽP., a void ŽD. and different kinds of pipes ŽE1–E9.. In particular, these sections illustrate the variations in the maximum depth of penetration and resolution with frequency. The radar wave penetrates down to the bottom of the pit with the 300 and 500 MHz antennas, as its reflection can be observed at about 85 ns, whereas with the 900 MHz antenna, the wave is entirely attenuated at about 50 ns. Taking the velocities calculated from the transmission measurements, the corresponding penetration depths can be estimated at around 4.5 and 2 m for the 300–500 and 900 MHz antennas, respectively. Similarly, the signal resolution, which we define here as a quarter of the ratio between the velocity and nominal frequency of the returned signal, can be estimated from these sections. For example, on profile 6 ŽFig. 4., the resolution increases from 0.03 to 0.14 m when considering 900 and 500 MHz antennas. Other GPR sections were similarly analysed in order to estimate such parameters ŽTable 1.. In addition, 53 parallel profiles, with a 10-cm spacing along the y-axis, were carried out with the 900 MHz antenna in compartment 4 so as to obtain a 3D data cube measuring 22 = 5.2 m. To reduce the volume of data without degrading the quality, data were resampled with inter-scans of 2 and 10 cm along the y- and x-axes, respectively, meaning that the final 3D data set was composed of 1101 = 53 scans. This experiment was performed in order to highlight the out-of-plane signals produced by side echoes, already described by Olhoeft Ž1994.. The study area is known to present longitudinal heterogeneities — such as pipes — along the y-axis. If these are perfectly cylindrical, the out-of-plane signals would be insignificant, and thus, the migrated and non-migrated sections in the Ž y,t . plane identical.

Fig. 5 shows two vertical sections through the data cube with no migration ŽFig. 5a., 2D migration in the Ž x,t . plane ŽFig. 4b., and 3D migration in the Ž x, y,t . volume ŽFig. 5c.. Fig. 5d indicates the location of the two sections and that of the buried objects. Analysis of these sections, particularly for the two signals enhanced by white arrows, reveals a remarkable difference depending on whether 2D or 3D migration was performed, with 3D migration offering a better focused signal than 2D migration.

Fig. 5. Schematic block-diagram of compartment 4 showing the buried objects and the vertical sections used to analyse the 3D data cube Žd.. Vertical sections through the 3D data cube without migration Ža., with 2D migration of the Ž x,t . plane Žb., and with 3D migrations of the Ž x,t . and Ž y,t . planes Žc.. The arrows highlight the signatures of a single pipe ŽO. and a set of pipes ŽE..

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

This confirms the existence of out-of-plane signals, a phenomenon that is observed every time a GPR

149

profile crosses a structure that is not perfectly cylindrical, which is most often the case.

Fig. 6. CMP gather recorded along the x-axis at a y-distance of 13 m Ža., with the corresponding semblance diagram Žb., and interval velocity curve Žc.. Interval velocity curves calculated for all the CMPs were superposed onto the interpolated velocity field Žd.. Migrated section of the compartment 2 along the x-axis and borehole information for comparison Že.. Note the high velocity layer at a depth of 1.5 m corresponding to the gneiss at B 14r20 mm.

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

150

2.4. CMP measurements CMPs were recorded in compartment 2, which is composed of seven horizontal layers of the different materials. The bistatic acquisition device consisted of two antennas that were moved symmetrically apart from a central point. The resulting data set — CMP gathers — is composed of a series of scans with the same CMP, and recorded with increasing transmitter–receiver offsets. This measurement configuration is used to estimate velocity variations with depth by applying the appropriate Normal Move Out ŽNMO. corrections to each scan ŽFisher et al., 1992; Tillard and Dubois, 1995; Greaves et al., 1996.. In our case, this operation was repeated so as to obtain 13 different CMPs distributed along the x-axis with measurements realized for each using 500 MHz antennas every 0.2 m with offsets ranging from 0.4 to 5 m. Fig. 6a shows an example of a CMP gather. NMO analyses were processed using Eq. Ž6., giving the expression of the time correction D t as a function of velocity V and offset x:

ž

Dtst 1y

x2 V 2t2

1r2

/

.

Ž 6.

The appropriate velocity law was then estimated from the semblance diagram, calculated with ŽNeidel and Taner, 1971.:

žÝ /

3. Synthesis and discussion We shall now compare and discuss the data recorded at the LCPC test site using the different GPR acquisition techniques and the results obtained from the various types of processing. After a synthesis concerning radar wave penetration, resolution, velocity and attenuation for each sounded material, we will consider the effective contribution of GPR applied to civil-engineering auscultation. 3.1. On penetration and resolution

2

ny1

maxima of semblance diagrams, and the second relating to lateral velocity variations due to the differential compaction of materials constituting each layer. In order to image appropriately the monostatic section recorded in compartment 2, a smoothed velocity field was computed from previous velocity analyses ŽFig. 6d.. A correlation between borehole information and the migrated section ŽFig. 6e. indicated that most of the reflections are related to dielectric contrasts due to lithology. Because the materials in compartment 2 are the same as those used to fill the other compartments, the CMP averaged velocities and those calculated from the tomograms or estimated from the 2D profiles are compared in Table 1.

AŽ t , j .

SŽ t . s

js0 ny1

n

Ý AŽ t , j .

.

Ž 7.

2

js0

For each CMP, the corresponding semblance diagram ŽFig. 6b. was calculated from Eq. Ž7. giving the best velocity law found from the NMO corrections. These velocities were then converted to interval velocities ŽFig. 6c. according to the Dix equation ŽDix, 1955.. Depending on the CMP gather, the number of constant velocity layers varies from five to seven, and depending on the layer considered, the interval velocity ranges from 0.07 to 0.17 mrns. Lateral velocity variations are observed within layers. Two principal causes can be proposed, the first involving the processing sequence and essentially due to errors introduced when picking velocities as

The GPR profiles in Fig. 4 show reflections from the boundaries of the test site ŽP., and some diffractions from a void ŽD. or pipes ŽE1–E9. buried in compartment 3. The frequency effect of the source is clearly illustrated in terms of penetration and resolution. In order to estimate the penetration and resolution of radar waves according to the physical characteristics of the different compartments, we tried to match observed GPR profiles with those calculated by forward modeling. The advantage of this method is the full integration of the wave propagation phenomena. The modeling principle, taken from Bitri and Grandjean Ž1998. and fully presented in the Appendix, is based on upward extrapolation of the wavefield in the frequency–wavenumber domain by using a phase-shift technique. The medium parame-

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

ters are explicitly introduced in the numerical grid as the relative permittivity: Ks

´0 ´ void

,

and the quality factor: Qs

1 s tan d

v´ X Ž v . s q v´ Y Ž v .

,

taken as the reverse of the loss tangent ŽBano, 1996.. It was demonstrated in these studies that the medium parameters are respectively associated to propagation parameters, namely velocity, attenuation and wave dispersion. As suggested by Turner and Siggins Ž1994., the approximation of a constant Q attenuation model leads to consider the Q values as valid only in a restricted frequency bandwidth. Consequently, identical models calculated with very different frequency bandwidths should not be comparable. A trial and error match of the modeled and observed sections was performed to estimate these quantities

151

characterizing the materials composing the site. The position and nature of the buried objects was known and imposed in these models. The section presented in Fig. 4c was modeled according to the dielectric models presented in Fig. 7a and b; the corresponding synthetic section is presented in Fig. 7c. The dielectric models were adjusted so that the curvature of the hyperbolas and the limit of GPR signal penetration respected Fig. 4c, indicating a good estimation of the wave velocity and attenuation. This procedure was repeated for all the GPR profiles. The values of the K and Q parameters attributed to the different media of the compartments, and the related parameters for the extreme 300 and 900 MHz frequencies are presented in Table 1. The response to the heterogeneities also merits analysis because this depends on resolution compared to size. When objects are sufficiently large compared to wavelength Že.g. blocks, voids, etc.., their interpretation is facilitated because the reflected signal contains complete information about their

Fig. 7. Relative permittivity Ža., and Q factor distribution Žb. used to model the 2D 500 MHz GPR profile recorded in compartment 3 Žcompare with Fig. 3c., and the computed synthetic section Žc.. The method for estimating velocity from hyperbola curvature and the limit of penetration are indicated. P s pit limit reflection; D s void diffraction; E1–E9 s diffraction of three sets of pipes ŽE1–E4–E7, E2–E5–E8, and E3–E6–E9. made respectively of iron, PVC filled with water, and PVC filled with air.

152

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

shape. On the contrary, when objects are too small compared to wavelength, they are difficult to identify because diffraction phenomena occur, causing complex interactions between the wave and the objects. Consequently, no major difference was observed between the GPR signature for the cables and the small pipes, except for the pipe containing water. In this case ŽFig. 4c., pipe diameter is 4 cm, with pipe E1 made of metal iron and pipes E2 and E3 made of PVC containing water and air, respectively. The GPR responses to pipes E1 and E2 show one and two hyperbolas, respectively. For E1, the observed signal represents diffraction of the wave on the metallic pipe, whereas for E2, echoes from the top and bottom of the pipe are observed, which occurs when the medium inside the pipe Žwater in this case. has a sufficiently low velocity to separate the two echoes of the returned signal. In comparison, the signal of pipe E3 containing air shows a different pattern because the velocity through air is higher, thus rendering the two echoes mixed and indiscernible. 3.2. On Õelocity and attenuation The velocity values obtained from surface–borehole tomography, 2D profiling and CMP analysis are summarized in Table 1. For each compartment of the site, the deviation between the velocity values estimated from the three considered techniques is only a few centimeters per nanosecond, which indicates coherent results and a good reliability for each method. However, local heterogeneities in the medium can produce a dispersion of the attenuation and velocity values, leading to some difficulties in recovering the main structures. For example, Fig. 6d and e shows the complexity involved in correlating the information derived from interval velocity curves with the a piori known velocity layers. The high velocity value Ž0.15–0.17 mrns. observed for compartment 4 using the different techniques could be explained by the high porosity of the gneiss at B 14r20 mm. The effective characteristics of this material, considered as a mixture of bulk rock and air, can be approached by a simple Lichtenecker’s law ŽOlhoeft, 1980.: K s K 1x 1 K 2x 2 . . . K ix i

Ž 8.

where K ) is the effective permittivity, and K i and x i the permittivity and volume fraction of the ith medium, respectively. Assuming that bulk gneiss permittivity is that of compartment 5 Ž6.9 — gneiss at B 0r20 mm and thus low porosity., that air permittivity is 1, and that K ) s 3 Ži.e. the permittivity measured in compartment 4., Eq. Ž8. gives a porosity for compartment 4 of about 0.3, which is in good agreement with the sample measurements. On the contrary, the values of attenuation factor estimated from the attenuation tomography and 2D profiling show a higher disparity, insofar as the former are systematically lower than that latter. Because the surface–borehole and 2D surface surveys were carried out respectively with low Ž100 MHz. and high Ž300, 500 and 900 MHz. antennas, this disparity could result to the approximation introduced in the constant Q attenuation model ŽTurner and Siggins, 1994.. This model assumes that Q, taken as the slope of the attenuation factor a vs. frequency, is constant for restricted frequency bandwidths, but not for a wide range in frequency, for example 100 to 900 MHz. The increasing value of Q with frequency observed here has already been described by Powers and Olhoeft Ž1994. and Hollender and Tillard Ž1998.. This shows the limitations of the constant Q model and prevents representation of a material attenuation by a single parameter because it also depends on the frequency used. In the case of the LCPC test site, the maximum depth of penetration and resolution vary respectively from 1 to 5 m and from a few centimeters to 0.25 m, depending on the considered frequency and the dielectric properties of the medium. The present study has demonstrated the potentiality of three different GPR techniques when applied to civil-engineering investigations, provided that the appropriate recommendations are respected. When boreholes are available, velocity and attenuation tomography is suitable for estimating the dielectric properties of a medium. Heterogeneities can also be imaged provided that they are larger than the wavelength. The ray coverage and angular resolution used in the reconstruction algorithm must also be adapted. 2D surface profiling is best suited to quickly imaging small diffracted objects such as pipes, voids, etc. During migration of the GPR sections from time v

v

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

to depth domain, out-of-plane signals generate artifacts when the structures are not perfectly cylindrical. For advanced studies, modeling algorithms can be used to estimate dielectric properties of the medium and heterogeneities. For sub-horizontally layered media, CMP processing by NMO analysis and semblance diagrams can be used to recover velocity variations with depth. This technique offers a resolution in relation to the frequency used. In our case, the 100 MHz signal has too large a wavelength to resolve the velocity layers, but the velocity distribution found gives a first approximation of the medium properties. Consequently, the GPR experiments carried out on the LCPC test site led us to consider the different capabilities of these techniques for civil-engineering sounding applications in urban environments. Although velocity and attenuation distributions can be efficiently recovered by surface–borehole tomography, boreholes are generally not available. In this case, velocity distribution with depth can be successfully estimated from CMP measurements, provided that the medium is sub-horizontally layered. Monostatic 2D surface profiling is well adapted to quickly locating underground heterogeneities with a relatively high reliability and, provided that control by modeling is carried out, velocity and attenuation data can be approximately estimated. v

153

tomography is suitable for estimating the dielectric properties of a medium. Heterogeneities can also be imaged depending on limitations due to aperture-related resolution and Fresnel zone dimension. CMP analysis can be used to recover the velocity field if the substratum is sub-horizontally layered. This information is then used to correctly image the monostatic sections. The performance of each technique is mainly conditioned by the material properties. In our case, depending on the kind of material sounded Žsilt, limestone sand, gneiss with different grain size. and the source frequency used, resolution and attenuation vary from a few centimeters to 0.25 m and from 2.5 to 45 dBrm, respectively. The penetration depth varies from approximately 1 to 5 m, making GPR a convenient tool for civil-engineering sounding applications.

Acknowledgements This work was carried out with the support of BRGM and the National Project on Trenchless works. We would like to thank LCPC for letting us use the test site, and for the technical support during the acquisition campaign. We also thank Rowena Stead for careful editing of the final manuscript.

Appendix A. GPR modeling method 4. Conclusion In the framework of a study dedicated to estimating the performance of GPR for civil-engineering applications, three sounding techniques were tested on the LCPC test site at Nantes, France: surface– borehole tomography, 2D surface profiling and CMP analysis. After processing according to three techniques, the results are coherent in terms of velocity and attenuation values. Monostatic 2D surface profiling is the fastest and most reliable way of imaging sub-surface structures. Forward modeling can also be used to analyse this response in order to obtain more quantitative interpretations. In reality, the choice of the sounding technique depends on the kind of problem to be resolved and the field conditions. If boreholes are available, velocity and attenuation

The starting point is a solution of Maxwell’s equation in the frequency–wavenumber domain Ž v , k .. This expression gives the electric field E generated by a source located at depth z, propagating upward through a homogeneous medium prior to being recorded at the surface ŽBitri and Grandjean, 1998.: E Ž k x ,k y , z s 0, v . s E Ž k x ,k y , z , v . eyi k z z

Ž A1.

where k x , k y , and k z are the wavenumbers in the x, y and z direction, respectively, v is the angular frequency, and i is the imaginary unit in the complex space. Eq. ŽA1. is associated to the dispersion equation relating k z to k x , k y and k, defined as the wavenumber in the propagation direction:

(

k z s k 2 y k 2x y k 2y

Ž A2.

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

154

k is generally a complex quantity k s b q i a , where the real and imaginary parts refer respectively to the phase and attenuation factors, and depend on the dielectric properties of the medium according to the classical relation ŽWard and Hohmann, 1987.:

(

k s v 2´m y i vms

Ž A3.

where ´ , m and s are respectively the complex dielectric permittivity, the magnetic permeability and the electric conductivity. To take into account the dispersion process in which dielectric parameters are dependent on frequency, a Q constant relaxation model can be used, as mentioned by Turner and Siggins Ž1994.. The dielectric permittivity becomes ŽBano, 1996.: ny 1

v

´ Ž v . s´0

ž /

Ž yi .

v0

ny 1

Ž A4.

in which ´ 0 and v 0 are two constants taken respectively as the reference permittivity and frequency. The index n is related to the quality or dissipation factor Q: ns

2

p

v as

Vp

p tan

4

v

Ž 1 y n . and b s

Vp

Ž A8.

v for Vp s V0 Ž

v0

.Ž1yn.r2 defined as the phase veloc-

ity. Because we only consider the monostatic mode, i.e. the transmitter at the same location as the receiver, we can split the electric field E into its downgoing Ey and upgoing Eq components, and state that: E Ž k x ,k y , z s 0, v . s Eq Ž k x ,k y , z , v . eyi k z z Ž A9a . provided that the phase velocity is divided by two. This assumption refers to the exploding source model ŽClaerbout, 1985.. In the case of a ID medium, where variations of dielectric parameters are only observed along the z direction, we can express the above equation as: E Ž k x ,k y , zs z y D z , v . s R Ž k x ,k y , z . Eq Ž k x ,k y , z , v . eyi k z D z

Ž A9b.

X

tany1 Ž Q . where Q s

v´ Ž v .

Ž A5.

s q v´ Y Ž v .

By combining Eqs. ŽA3. and ŽA4. k, which is now dependent on frequency:

v

Ž ny1 .r2

ž /

(

k Ž v . s v m´ 0

v0

p = cos

p

Ž 1 y n . q isin

4

4

Ž 1 y n.

Ž A6.

Assuming that we define a reference velocity V0 as: V0 s

and the expressions of a and b becomes:

v´ 0 Ž v .

1r2

Re w k x

,

we find

where D z is the discretization step along z, and RŽ k x ,k y , z . is the space Fourier transform of the medium impedance contrasts. Using a phase-shift method ŽGazdag, 1978., the upgoing wave is convolved to the R function from z s z max prior to being extrapolated upward to z–D z until z s 0. Because the extrapolation in the Fourier domain demands a homogeneous medium in Ž x, y,D z . layers, two computational steps are needed. In a first instance, the medium properties are averaged in layers Ž x, y,D z ., so that the medium becomes horizontally stratified and Eqs. ŽA9a. and ŽA9b. can be applied. For that, the wavenumber k z m is calculated from the averaged medium properties ´ m , sm and mm . In a second instance, a correction term is applied to the extrapolated field in the Ž x, y, z, v . space to take into account the lateral variations of the medium. The exact correction term is: y z x yk y . D z e i Ž 'k z myk x yk y y 'k x yk 2

v kŽ v. s

V0

v

ž / v0

Ž ny1 .r2

p 1 q tan

ž

4

Ž 1 y n.

/ Ž A7.

2

2

2

2

2

Ž A10.

where k x y z is the exact wavenumber calculated with non-averaged dielectric parameters. In view of the difficulties in calculating the square roots, Eq. ŽA10.

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

is approximated with a zero order Taylor expansion, giving the approximated correction term: e iŽ k m zyk x y z .D z .

Ž A11.

Finally, the antenna spectral signature and its radiation pattern have to be taken into account. The antenna spectral content is integrated by multiplying the extrapolated field by the tapering function:

ž ž

G Ž v . s exp y ss

v y vc s

2

//

with

Ž v N y 1 . P Bv

Ž A12.

vc

where vc and v N are respectively the centre and Nyquist angular frequencies, and Bv is the bandwidth. The term related to the 3D antenna pattern is then applied to EŽ k x ,k y , z s 0, v . according to analytical solutions of a single element steady state dipole ŽArcone, 1995.: E Ii s K

cos u sin w

Ž 1 y n2 sin2w .

E Ie s yiK

1r2

y ncos u

cos u sin w

Ž n2 sin2w y 1 .

1r2

q incos u

E Hi s K cos w 1r2

Ž 1 y n2 sin2w . q ncos u sin u cos u 1r2 n Ž 1 y n2 sin2w . y cos u 2

=

cos 2u y

Ž 1 y n2 sin2w .

1r2

y ncos u

E He s K cos w 1r2

=

Ž n2 sin2w y 1 . y incos u sin u cos u 1r2 n Ž n2 sin2w y 1 . q icos u 2

cos 2u qi

Ž n2 sin2w y 1 .

1r2

q incos u

Ž A13.

where K is a constant related to the antenna, n is the refraction index in the medium, the symbols 5 and H and define the parallel and normal planes with

155

respect to the antenna displacement, e and i subscripts refer to the external and internal field, respectively, with respect to the critical angle.

References Annan, A.P., Scaife, J.E., Giamou, P., 1990. Mapping buried barrels with magnetics and ground penetrating radar. 60th Annu. Int. Mtg., Soc. Expl. Geophys., 422–423. Arcone, S.A., 1995. Numerical studies of the radiation patterns of resistively loaded dipoles. Appl. Geophys., 33, pp. 1–3, 39–52. Bano, M., 1996. Constant dielectric losses of ground penetrating radar waves. Geophys. J. Int. 124, 279–288. Bitri, A., Grandjean, G., 1998. Frequency–wavenumber modelling and migration of 2D GPR data in moderately heterogeneous dispersive media. Geophys. Prospect. 46, 287–301. Chazelas, J.L., Leparoux, D., Hollier-Larousse, A., 1997. A test site for geophysical methods. Proc. EEGS 3rd Meet. Environ. Eng., 387–390. Claerbout, J.F., 1985. Imaging the Earth’s Interior. Blackwell, 398. Dix, C.H., 1955. Seismic velocities from surface measurements. Geophysics 20, 68–86. Fisher, E., McMechan, G.A., Annan, A.P., 1992. Acquisition and processing of wide-aperture ground penetrating radar data. Geophysics 57, 1933–1942. Gazdag, Z., 1978. Wave equation migration with the phase-shift method. Geophysics 43, 1342–1351. Grandjean, G., Durand, H., 1999. Radar unix: a complete package for GPR data processing. Comput. Geosci. 25, 141–149. Greaves, R.J., Lesmes, D.P., Lee, J.M., Toksoz, ¨ M.N., 1996. Velocity variations and water content estimated from multioffset, ground-penetrating radar data. Geophysics 61, 683–695. Hollender, F., 1999. Interpretation de la distorsion des signaux ´ georadar propages Developpement d’une tomogra´ ´ et reflechis. ´ ´ ´ phie par bande de frequence. PhD, INPG, Grenoble, France. ´ Hollender, F., Tillard, S., 1998. Modeling ground-penetrating radar wave propagation and reflection with the Jonsher parametrization. Geophysics 63, 1933–1942. Ivansson, S., 1987. Cross-hole transmission tomography. In: Nolet, G. ŽEd.., Seismic Tomography. Dordrecht, Holland. Jackson, M.J., Tweeton, D.R., 1994. MIGRATOM-Geophysical tomography using wavefront and fuzzy constraints. US Bureau of Mines, Minneapolis, RI9497, 32. Leggett, M., Goulty, N.R., Kragh, J.E., 1993. Study of traveltime and amplitude time-lapse tomography using physical model data. Geophys. Prospect. 41, 599–619. Mason, I.M., 1981. Algebraic reconstruction of a 2D velocity inhomogeneity in the High Hazels Seam of Thoresby Colliery. Geophysics 46, 298–308. Neidel, N.S., Taner, M.T., 1971. Semblance and other coherency measures for multichannel data. Geophysics 34, 482–497. Olhoeft, G.R., 1980. Electrical properties of rocks. In: Touloukian, Y.S., Judd, W.R., Roy, R.F. ŽEds.., Physical Properties of Rocks and Minerals. McGraw-Hill, New York.

156

G. Grandjean et al.r Journal of Applied Geophysics 45 (2000) 141–156

Olhoeft, G.R., 1994. Modeling out-of-plane scattering effects. Proc. 5th Int. Conf. Ground Penetrating, 133–144. Olsson, O., Falk, L., Forslund, O., Lundmark, L., Sandberg, E., 1992. Borehole radar applied to the characterization of hydraulically conductive fracture zone in crystalline rock. Geophys. Prospect. 40, 109–142. Powers, H.M., Olhoeft, G.R., 1994. Modeling dispersive ground penetrating radar data. Proc. 5th Int. Conf. Ground Penetrating, 173–184. Powers, H.M., Olhoeft, G.R., 1996. Modelling the GPR response of leaking, buried pipes. Proc. SAGEEP, Environ. Eng. Geophys., 525–534. Rector, J.W., Washbourne, J.K., 1994. Characterization of resolution and uniqueness in crosswell direct-arrival traveltaime tomography using the Fourrier projection slice theorem. Geophysics 59 Ž11., 1642–1649. Tillard, S., Dubois, J.C., 1995. Analysis of GPR data: wave propagation velocity determination. J. Appl. Geophys. 33, 77–91. Tong, L.T., 1993. Application of ground penetrating radar to locate underground pipes. Terr. Atmos. Ocean Sci. 4, 171–178.

Turner, G., Siggins, A.F., 1994. Constant Q attenuation of subsurface radar pulse. Geophysics 59, 1192–1200. Valle, S., Zanzi, L., 1997. Traveltime radar tomography for NDT on masonry and concrete structure. Eur. J. Environ. Eng. Geophys. 2, 229–246. Vasco, D.W., Peterson, J.E., Ha Lee, K., 1997. Ground-penetrating radar velocity tomography in heterogeneous and anisotropic media. Geophysics 62, 1758–1773. Ward, S.H., Hohmann, G.W., 1987. Electromagnetic theory for geophysical applications. In: Nabighian, M.N. ŽEd.., Electromagnetic Methods in Applied Geophysics — Theory vol. 1 S.E.G, Tulsa, OK. Williamson, P.R., 1991. A guide to the limits of resolution imposed by scattering in ray-tomography. Geophysics 56 Ž2., 202–207. Williamson, P.R., Worthington, M.H., 1993. Resolution limits in ray tomography due to wave behaviour: numerical experiments. Geophysics 58 Ž5., 727–735. Zeng, X., McMechan, G.A., 1997. GPR characterization of buried tanks and pipes. Geophysics 62, 797–806.