Fracture thickness from GPR measurements

Rock investigation is definitely not a recent application of Ground Penetrating Radar (GPR) technique, as first studies date back to the seventies. However, only in the last decade research activities have started to address GPR characterization of rock fracture parameters, namely aperture and filling material. Rock fractures can generally be considered as thin beds, i.e., two interfaces whose distance is smaller than radar range resolution. Most of the past studies analyzed thin-bed response in the time domain, addressing time resolution, the linear relationship between bed thickness and reflected amplitude, and the derivative effect upon the incident signal. Amplitude calibration might permit to estimate fracture features for arbitrarily thin beds, but it is difficult to achieve and could be applied only to favorable cases. In this paper we explore the possibility to estimate fracture thickness and filling in the frequency domain by means of GPR. After reviewing the theoretical aspects of thin-bed response, we processed GPR data collected on ornamental marble blocks, where fractures of known aperture were simulated. We also performed numerical modelling tests to support the analysis of real datasets. Our approach consists of a 4-step procedure in which deterministic deconvolution is used to retrieve magnitude and phase thin-bed response in the selected frequency band. The procedure provided satisfactory outcomes when applied to real as well as to modelled thin-bed reflections. Results are encouraging and suggest that, under favorable circumstances, GPR could be a fast and effective tool to determine fracture parameters in non-destructive manner. Further testing is needed in order to fine-tune the processing sequence and to extend the validity of our preliminary findings to more complex case studies.


I. INTRODUCTION
Ground-penetrating radar (GPR) investigation of rocks is a well-known application and has been performed for nearly 40 years now.According to the desired trade-off between resolution and penetration depth, the full frequency range of commercial GPR systems (i.e., from tens of MHz to few GHz) has been employed in field investigations and laboratory experiments.Provided that water content is low and high conductivity minerals are not present, rocks generally offer favorable propagation conditions to electromagnetic radiation because of their relatively low absorption over the frequency band of GPR systems [1].In the last decades several GPR surveys have addressed the detection of fractures within rock bodies to characterize mining sites [2], quarries [3], rock masses [4], road cuts [5] and unstable rock slopes [6].In many instances the detection and location of fractures is obviously of great importance for safety reasons, but, for example, can also be valuable in the production of ornamental stones, both before and after rock blocks are quarried [3,7].In addition, some authors have been dealing with GPR characterization of rock fracture parameters, namely thickness and filling material [8].
Fractures can generally be envisaged as layers embedded in a homogeneous rock formation.This gives rise to two signals with opposite polarities reflected by the two sides of a fracture.If two distinct signals can be clearly identified and are well separated in time, fracture thickness can be directly determined by time difference measurements (provided velocity is known) and the fracture is referred to as a thick layer.On the contrary, the fracture is considered to be a thin layer, or a thin-bed, when the two reflections overlap in a way that just a single composite wavelet is identifiable.In such a case there may be constructive or destructive interference between reflections (tuning effect), and the information about fracture parameters is encoded in the features of the composite reflection.The difference between a thick and a thin layer is dictated by the range resolution of a radar system, i.e. its ability to distinguish two or more targets on the same bearing but at different ranges.For an impulse radar like GPR, transmitted pulse width is the primary factor in range resolution, and it is widely accepted that two pulses are said to be resolved if their envelopes are separated by half the envelope half width [1].Radar range resolution limit is derived from the Rayleigh's resolution criterion, and becomes the wellknown quarter-wavelength limit when the dominant frequency of the propagating wavelet is used to compute the wavelength.Relatively high propagation velocity of the electromagnetic energy within rock fractures coupled to the frequency range of commercial GPR antennas, yields wavelengths that are large compared to the thickness of fractures commonly encountered in surveys.For this reason, most rock fractures may be considered as thin layers from a GPR point of view.In such a case, thin-bed response is the sum of the primary reflections and multiple reflections due to the reverberation of the signal back and forth within the layer.The composite reflected event is given by the interference between reflected and transmitted signal replicas occurring at fixed time delays and scaled according to the Fresnel's reflection and transmission coefficients of the interfaces.Thin-bed response can be expressed in terms of magnitude and phase (Fig. 1) as follows [2,9] |Γ| = (2Rsin(2πd/λ)) / ((1-R 2 ) 2 +(2Rsin(2πd/λ)) 2 ) 1/2   (1) where R is the normal-incidence reflection coefficient between the medium and the thin-bed, d is the layer thickness and λ is the dominant wavelength of the signal within the bed.Thin-bed magnitude response as a function of frequency oscillates between 0 and a maximum value corresponding to destructive and constructive interference between primary reflections and multiples respectively.If layer material has some loss, the oscillations damp out as layer thickness increases.A similar trend can be observed when analyzing thin-bed response as a function of bed thickness; maxima and minima occur at thicknesses multiple of λ/4 and λ/2 respectively (Fig. 1).In fact, the fracture acts as a frequency filter whose mask is dependent upon the frequency band of the incident signal and on the properties of the layer.It is easy to show that, when the bed is very thin compared to the wavelength, the response of the thin-layer is simplified into [1] with i being the imaginary unit, v the velocity of the wave within the bed and f being the frequency.Equation ( 3) points out the linear relationship between frequency (or thickness) and reflected amplitude and the 90° phase-shift of the reflected signal (Fig. 1), i.e., the thin-bed acts as a derivative.In the following we examine thin-bed response through laboratory and synthetic GPR experiments.Analysis of the results is performed both in time and frequency domains with the aim of giving new insights for rock fracture characterization.

A. Laboratory Experiments
We carried out laboratory GPR measurements with a K2 IDS system and a dual-polarized antenna with a nominal frequency of 2GHz.Radar traces were collected setting a timetriggering mode with a stacking factor of 4, a 10ns acquisition window and a sampling frequency of about 70GHz.We employed 0.75m x 0.60m x 0.21m Carrara marble blocks, quarried in Tuscany, central Italy.We first performed radar measurements to determine the velocity of the electromagnetic signal within the marble that was estimated to be 9.67cm/ns, thus giving a relative permittivity of about 9.63.Low-loss approximation is reasonable because data collected on blocks with different thicknesses showed that spherical divergence is essentially the main factor causing wave attenuation.Afterwards we collected a reference reflection by placing a metal screen on the opposite side of the block where the antenna was placed.By comparing amplitudes of the signals reflected by the marble/air and the marble/metal interfaces we were able to compute the marble/air reflection coefficient to be 0.51, in accordance with the value of the permittivity estimated previously.To simulate an air-filled fracture, acquisitions were carried out using just two marble blocks one in front of the other.The thicknesses of the air thin bed to be tested were computed using the dominant wavelength of the signal reflected at the marble/air interface.Though this procedure is not correct theoretically, we considered it to be a good approximation assuming there is little or no frequencydependent attenuation during signal propagation in the marble block.The dominant wavelength was computed to be 21.58cm,corresponding to a dominant frequency of about 1.39GHz.Radar traces were collected for apertures ranging from λ/32 to λ/2 at λ/32 step.The most favorable polarization is the one with the dipoles having their longer axis perpendicular to the longer side of the block.This is quite obvious because their directivity patterns are less disturbed by the edges of the block.

B. Modelling Tests
In order to support the analysis of GPR data collected in the lab, we designed numerical tests in which thin-bed response to incoming GPR signals was explored.The propagation of electromagnetic pulses generated by GPR was modeled by solving Maxwell's equations using a 2D Finite-Difference Time-Domain (FDTD) code [10].According to the set-up of the laboratory GPR experiments, we considered two marble blocks at different distances from one another in order to simulate the presence of an air-filled fracture with variable thickness.Electromagnetic parameters were considered to be frequency-independent, neither conductivity nor magnetic properties are taken into account.Relative permittivity is set to 1 and 9.63 for air and marble respectively.Excitation of the model was achieved by supplying a current to a line source placed on top of one block and electromagnetic signals were collected by a receiving element adjacent to the transmitting one to have nearly zero-offset acquisition.A 2-loop Ricker wavelet with a peak frequency of 0.83GHz was fed to the transmitting element and, as a result, a mixed-phase 3-loop wavelet was collected.The peak frequency of the radiated wavelet was tuned to have the same dominant frequency of the actual radiated pulse (1.39GHz).The discretization step of the model was set to 1 mm, so as to properly resolve the thinnest fractures and to limit numerical dispersion.Data were computed in the H-plane.We also modeled the reflection from the semi-infinite air half-space as well as from the marble/metal interface to obtain reference traces for the study of the thin-bed reflection features.

III. DISCUSSION
The analysis of modeled and collected traces suggests that as fracture gets thinner thin-bed reflections tend to resemble the first derivative with respect to time of the incident wavelet, as expected.This can also be shown by means of a crosscorrelation approach [11].Considering that below the Rayleigh's limit (λ/4) amplitude information encodes thickness variations, and provided the entire amplitude variation is caused by tuning effects, amplitude calibration may permit thickness calculations for arbitrarily thin-beds.This is true for the case of fractures at constant range and with the same electromagnetic properties of the filling material.If this not the case, all factors affecting the reflected amplitude must be carefully evaluated.They may include absorption, spherical divergence, presence of additional reflectors between the radar and the analyzed fracture, differences in the radiated signal due to variations in antenna coupling and to radiation pattern, along with all the factors affecting the reflection coefficient, namely frequency, antenna polarization (this implicitly including the orientation of GPR with respect to the target) and variations between the electromagnetic properties of the host and the fracture filling materials.Provided that the fracture is sufficiently far from the acquisition plane, normal incidence condition can generally be assumed (i.e., nearly zero-offset acquisition).For all these reasons, amplitude analysis of thinbed reflection in the time domain is not straightforward and may be applied only in very favorable cases, for instance on small quarried blocks.
As far as frequency domain is concerned, thin-bed frequency response may be obtained by means of deterministic deconvolution if the propagating wavelet can be extracted from the collected radar data.By performing the ratio of thin-bed reflection to source wavelet complex spectra, both amplitude and phase deconvolutions are obtained.After deconvolution has been performed, recovered amplitude and phase curves can be fitted in order to estimate fracture parameters.Regarding amplitude spectra, the use of normalized curves (Fig. 1b) has to be preferred so as to discard factors affecting the amplitude of the collected signal, but not related to thin-bed features.Moreover, assuming low-loss materials and electromagnetic parameters not dependent upon frequency in the band of interest, the reflection coefficient at an interface can be expressed in terms of the permittivities of the two materials.Hence, the fitting procedure would provide an estimate of fracture aperture and permittivity of the filling material.It is clear that the fitting process would be significantly sensitive to thickness, while the inversion of the dielectric constant would be weakly constrained (Figs.1b and 1c).More in detail, fracture thickness controls the position along the frequency axis of minimum and maximum amplitude values, as well as of zero phase points.On the contrary, the reflection coefficient (in this case the dielectric constant) mainly controls just the dip (i.e., the first derivative) of both amplitude and phase curves.In keeping with this, we propose a 4-step assessment procedure as follows: i) identify a reference signal and its dominant frequency; ii) perform deterministic deconvolution to obtain thin-bed amplitude and phase response in the selected frequency band; iii) fit obtained curves with analytical relationships (Eqs. 1 and 2) by means of a least squares approach to determine fracture thickness in terms of wavelength (Figs.2a and 2c), i.e., regardless of the dielectric constant of the filling material (this is actually carried out with an arbitrary permittivity as reflection coefficient has very little influence on the processed curves); iv) fit deconvolved spectra considering a small range of thicknesses obtained from the previous step, and a reasonable set of possible permittivities (Figs.2b and 2d).It is obvious that, whenever the analyzed thin-bed reflection reveals a polarity inversion with respect to the radiated signal, the explored permittivities could be significantly reduced.The procedure yielded good results in terms of both thickness and dielectric constant and no significant differences were observed between modeled and real datasets (Fig. 2).Anyway, it must be stated that fitting errors were lower for synthetic data.Thickness estimate with phase curves is slightly better than the one obtained by considering amplitude ones.Instead, both curves provide the same results as far as the permittivity is concerned; the best fit is always obtained for the relative dielectric constant of air.Deconvolved amplitude and phase curves provide similar results.The accuracy of the inversion process is mainly dictated by the band over which the analysis can be performed together with the aperture of the fracture, which in turn delimit the section of the curves that are used in the inversion process.Reasonably, the wider the band, the higher the accuracy.

IV. CONCLUSIONS
In this work we proposed a methodology to estimate rock fracture thickness and filling material based on the analysis of common-offset GPR data.We address the thin-bed frequency response because we deem it to be more robust than thin-bed amplitude analysis in the time domain.The approach implies deterministic deconvolution and provided satisfactory outcomes when applied to both modeled and collected data.Dissimilarities between the two datasets could be due to small errors in setting bed thickness, noise and side-wall reflections affecting real acquisitions.In addition, the real source wavelet is likely to be more complex than a simple 3-loop Ricker wavelet and may alter the way in which the signals reflected by the bed interfere to make up the composite reflection.As a matter of fact, we observed that the real wavelet has the same dominant frequency of the synthetic one, but narrower bandwidth.As a result, the real wavelet has higher side lobes and longer duration.Also, the differences between modeling results and the analytical expressions (Eqs. 1 and 2) are due to the shape and duration of the incident waveform which are directly related to its frequency characteristics, namely peak frequency and spectral bandwidth.Processing steps such as Wiener filtering may be needed to shape the collected waveforms before running the inversion process.If deterministic deconvolution cannot be performed because of the lack of a reference signal, statistic deconvolution may also be attempted.

Fig. 1 .
Fig. 1.Thin-bed response as a function of fracture thickness and reflection coefficient.(a) Amplitude response, (b) amplitude response normalized with respect to the maximum value along the thickness axis for each reflection coefficient, and (c) phase response.

Fig. 2 .
Fig. 2. Estimation of normalized fracture thickness (a) -(c) and dielectric constant (b) -(d) of fracture filling material for synthetic (a) -(b) and real (c) -(d) datasets.Black lines are the deconvolved spectra.Color lines are theoretical expectations from eqs. 1 and 2 for different parameters according to the legends shown on the right column.