Rock fracture characterization with GPR by means of deterministic deconvolution

Highlights

• A procedure to estimate thickness and filling of fractures is proposed.

• GPR response to thin layers is modeled and thin bed reflections are processed resorting to deterministic deconvolution.

• Processed data are analyzed in the frequency domain.

Abstract

In this work I address GPR characterization of rock fracture parameters, namely thickness and filling material. Rock fractures can generally be considered as thin beds, i.e., two interfaces whose separation is smaller than the resolution limit dictated by the Rayleigh's criterion. The analysis of the amplitude of the thin bed response in the time domain might permit to estimate fracture features for arbitrarily thin beds, but it is difficult to achieve and could be applied only to favorable cases (i.e., when all factors affecting amplitude are identified and corrected for). Here I explore the possibility to estimate fracture thickness and filling in the frequency domain by means of GPR. After introducing some theoretical aspects of thin bed response, I simulate GPR data on sandstone blocks with air- and water-filled fractures of known thickness. On the basis of some simplifying assumptions, I propose a 4-step procedure in which deterministic deconvolution is used to retrieve the magnitude and phase of the thin bed response in the selected frequency band. After deconvolved curves are obtained, fracture thickness and filling are estimated by means of a fitting process, which presents higher sensitivity to fracture thickness. Results are encouraging and suggest that GPR could be a fast and effective tool to determine fracture parameters in non-destructive manner. Further GPR experiments in the lab are needed to test the proposed processing sequence and to validate the results obtained so far.

Introduction

Ground-penetrating radar (GPR) investigation of rocks is a well-known application and has been performed for nearly 40 years now (Stewart and Unterberger, 1976). 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 (Sambuelli and Calzoni, 2010, Sassen and Everett, 2009). 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 (Annan, 2001). In the last decades several GPR surveys have addressed the detection of fractures within rock bodies to characterize mining sites (Annan et al., 1988), quarries (Grasmueck et al., 2010), rock masses (Pipan et al., 2003), road cuts (Longoni et al., 2012) and unstable rock slopes (Heincke et al., 2005). In many instances the detection and location of fractures are 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 (Arosio et al., 2012, Lualdi and Zanzi, 2004). In addition, some authors have been dealing with GPR characterization of rock fracture parameters, namely thickness and filling material (Deparis and Garambois, 2009, Grégoire and Hollender, 2004).

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 the time difference between reflections in a single measurement (provided that the velocity of the layer filling 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; Kallweit and Wood, 1982, Sheriff, 1985), and the information about fracture parameters is encoded in the features of the composite reflection. More in detail, the amplitude of the composite reflected wavelet, together with its amplitude and phase spectra, depends on the source wavelet, on the time delay between the reflections as well as on the magnitude and polarity of the reflection coefficients.

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 (Annan, 2001). For an impulse radar (i.e., an electromagnetic device that emits short-duration pulses) 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 in time by half the envelope width. The envelope width is commonly defined as the time between the points on the envelope where the envelope is greater than half of its peak amplitude (Annan, 2001). Radar range resolution is derived from the Rayleigh's resolution criterion, which was originally defined with an optical analogy and set the threshold for the resolvable limit. Rayleigh's criterion states that the limit of an optical device to distinguish separate images of objects lying close to one another occurs when the two diffraction images are separated by a distance equal to the peak-to-trough distance of the diffraction pattern of the device (Kallweit and Wood, 1982). If we define the dominant frequency as the reciprocal of the time between successive peaks or troughs of the propagating wavelet (also known as dominant period or breadth), Rayleigh's criterion converts to the well-known quarter-wavelength criterion. It is worth pointing out that this widely accepted criterion for range resolution is in actual fact a reasonable but arbitrary choice, and depends on the Signal-to-Noise Ratio (SNR), the duration and shape of the wavelet, as well as the analytical tools used to examine the collected data (Ricker, 1953). Although other researchers have tried to provide different resolvable limits (Ricker, 1953, Widess, 1973), it has been shown that Rayleigh's peak-to-trough time separation corresponds to the practical limit (Kallweit and Wood, 1982).

Relatively high propagation velocity of the electromagnetic energy within rock fractures coupled to the frequency range of commercial GPR antennas, results in wavelengths that are large compared to the thickness of fractures commonly encountered in surveys, which are generally few cm wide at most. 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 (Fig. 1a). 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 reflection coefficient as a function of frequency is related to the ratio of bed thickness to wavelength within the bed. In case of layers much thinner than the wavelength, reflection coefficient magnitude varies linearly with frequency, while, as the layer gets thicker, it oscillates between 0 and a maximum value corresponding to destructive and constructive interference between primary reflections and multiples respectively. If layer material is lossy, the oscillations damp out as layer thickness increases. A similar trend can be observed when analyzing the reflection coefficient as a function of bed thickness; maxima and minima occur at thicknesses that are multiples of λ/4 and λ/2 respectively. Fig. 1b and c shows amplitude and phase of the normal-incidence thin bed reflection coefficient as a function of frequency and thickness respectively, computed using the analytic expressions (see Eqs. (4), (5)) for air and water thin beds embedded in sandstone. Absorption is only considered for the water thin bed, where the dielectric constant is modeled according to the Debye formula (Debye, 1929, Section 18). 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.

Response of thin layers has been largely investigated by the seismic industry and, more recently, by the GPR community as well (Bradford and Deeds, 2006, Tsoflias and Hoch, 2006). From a theoretical point of view, Lord Rayleigh pioneered the discussion by exploring the response of a plate of finite thickness to an incoming plane sinusoidal wave (Rayleigh, 1945, Section 13). More in details, for a plate of thickness d and assuming normal incidence, thin bed response is given by

$$\mathrm{\Gamma }=\frac{\alpha -{\alpha }^{-1}}{\alpha -{\alpha }^{-1}-i2cot\left(2\mathit{\pi d}/\lambda \right)}$$

where i is the imaginary unit, λ is the wavelength within the bed and

$$\alpha =Z_2/Z_1$$

is the ratio between the acoustic impedances of the thin layer and the surrounding medium respectively. Annan et al. (1988) proposed a relationship when dealing with transient signals of pulse radar

$$\mathrm{\Gamma }=\frac{R_{12}\left(1-exp\left(-i4\mathit{\pi d}/\lambda \right)\right)}{1-R_{12}^{2}exp\left(-i4\mathit{\pi d}/\lambda \right)}$$

where R₁₂ is the normal-incidence electromagnetic reflection coefficient between the medium and the thin bed and λ may be regarded as the dominant wavelength of the radar pulse within the bed. Though dealing with Acoustics and Electromagnetics respectively, both approaches are based on the same phenomenon which assumes an infinite series of waves bouncing back and forth within a thin layer. In fact, the relationships above are identical, and with proper changes to the notation, they can be expressed in terms of magnitude and phase (see Appendix A)

$$\left|\mathrm{\Gamma }\right|=\frac{2R_{12}sin\left(2\mathit{\pi d}/\lambda \right)}{\sqrt{{\left(1-R_{12}^2\right)}^2+{\left(2R_{12}sin\left(2\mathit{\pi d}/\lambda \right)\right)}^2}}$$

$$\angle \mathrm{\Gamma }=\vartheta ={tan}^{-1}\left(cot\left(2\mathit{\pi d}/\lambda \right)\frac{1-R_{12}^2}{1+R_{12}^2}\right)=\frac{\pi }{2}-{tan}^{-1}\left(tan\left(2\mathit{\pi d}/\lambda \right)\frac{1+R_{12}^2}{1-R_{12}^2}\right).$$

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 (Annan, 2001, Koefoed and de Voogd, 1980)

$$\mathrm{\Gamma }\cong \left(\frac{4\pi R_{12}}{1-R_{12}^2}\frac{d}{v}\right)\mathit{if}$$

where v is the velocity of the wave within the bed and f is the frequency. This equation points out the linear relationship between frequency (or thickness) and reflected amplitude and the 90° phase-shift of the reflected signal, highlighted in the right-hand side of Eq. (5) as well.

In the scientific literature thin bed response has been analyzed according to different perspectives concerning time resolution (Kallweit and Wood, 1982, Neidell and Poggiagliolmi, 1977, Ricker, 1953, Sheriff, 1985), the linear relationship between bed thickness and reflected amplitude (Annan et al., 1988, Chung and Lawton, 1995b, Koefoed and de Voogd, 1980, Widess, 1973), the dispersion of the reflected event (O'Neill, 2000), and the derivative effect upon the incident signal (Annan, 2001, Chung and Lawton, 1995a, de Voogd and den Rooijen, 1983, Widess, 1973).

In the following I study thin bed response through synthetic GPR experiments. Analysis of the results is performed both in time and frequency domains to seek possible matching with findings of previous studies and to improve approaches to rock fracture characterization. I propose a 4-step procedure in which deterministic deconvolution is used to retrieve magnitude and phase thin bed response in the selected frequency band. After deconvolved curves are obtained, fracture thickness and filling are estimated by means of a fitting process.