Complex analysis of GPR signals for the delineation of subsurface subtle features

In this paper, complex signal analyses of ground penetrating radar (GPR) field data over an area of farmland in Krakow were interpreted alongside the basic filtered field data. The farmland was simulated with varying degrees of soil compaction induced by tractor movement. The focus of the study was the delineation of inherent characteristics of media through which the electromagnetic energy travelled. Fourteen GPR profiles were acquired from the area. The field data were subjected to preand post-processing prior to its the presentation and interpretation. Advance analysis operations on the field data which resorted in different attributes reveal more about the effects of the compaction on the soil than indicated by the basic filtered field data. Better resolution of subsurface layers boundary and lateral variation in the physical properties of the traversing media were well elucidated. The results have demonstrated that an advanced signal processing such as used in the study has ability to depict subtle characteristics of the propagating media.


INTRODUCTION
The ability to translate electromagnetic (EM) field measurements into useful engineering or scientific information is critically important to the widespread use of the technique in Ground Penetrating Radar (GPR). However, random noise influence may in some cases obliterate subtle, yet very important, features that may have been encountered by the EM pulse energy. This may then mean that visual analysis and interpretation of the signal response may not extract all of the detailed information of the media through which the signal traverses.
The uppermost part of the Earth's crust that is biologically active and porous constitutes what is called soil. It is one of the principal substrata of life on Earth, serving as a reservoir of water and nutrients, as a medium for the filtration and breakdown of injurious wastes, and as a participant in the cycling of carbon and other elements through the global ecosystem. It has evolved through weathering processes driven by biological, climatic, geological and topographical influences. Soil can be said to be a horizon within the subsurface or near surface layer that plays vital roles for both plants and animals dwelling on it, including humans. Adequate knowledge of soil constitutive parameters and how they are related to their optimum utilization cannot be overemphasized. Among the various tools used for field studies of the near surface layer, in which soil is a part, is Ground Penetrating Radar (GPR) method.
GPR is a geophysical method that uses radar pulses to image the subsurface. It is a non-destructive method that uses electromagnetic radiation in the radio spectrum to detect the reflected signals from the subsurface.
A GPR transmitter emits electromagnetic energy into the ground. When the energy encounters a buried object or a boundary between materials having different dielectric permittivity, it may be reflected or refracted or scattered back to the surface. A receiving antenna can then record the variations in the return signal (Daniels 2004, Jol 2009, Cerquera 2017. As applicable in all the major geophysical survey methods, geoscientists are mainly interested in geophysical anomalies, i.e., local variations in the measured parameter relative to some normal background value. Such variations are attributable to a localized subsurface zone of distinctive physical properties and possible geological importance. However, a geophysical interpretation may be ambiguous. Moreover, the impact of attenuation and random noise affecting the geophysical data may limit the interpretation of it. This paper focuses on the advance interpretation of ground penetrating radar which may be used to limit the problems of attenuation and non-uniqueness inherent in the mere visual assessment of GPR field data. Qiao et al. (2015) employed a Multi-resolution Monogenic Signal Analysis (MMSA) system in GPR images to locate metal objects. The use of a wavelet to reduce noise and a fuzzy cluster approach has been used by Delbo et al. (2000) to detect buried objects. Delineation of light non-aqueous phase liquids (LNAPL) floating on water table was made possible via the analysis of field data signals by Orlando (2002). A wavelet transform approach for signal transformation (filtering) as opposed to the Hilbert transform has also been proposed by Liu & Oristaglio (1998). The depth to compacted layer estimation using GPR data has been recorded in the literature by Muñiz et al. (2016). Analyses of spectral features such as power spectral density (PSD), short-time Fourier transform (STFT) and Wigner Ville distribution (WVD) on simulated data for detection of features in GPR data was carried out by Vinicius et al. (2013). Jonard et al. (2012) have reported the evaluation of the effects of different tillage practices on soil physical properties such as soil water content (SWC) by using geophysical methods which include GPR. A detailed description of use of seismic and GPR attributes can be found in Barnes (2016) and Tomecka-Suchoń (2019).
In this article, an attempt has been made to correlate the results of advance signal processing with basic filtered field data with a view to probing deeper into the subtle information, which the analysis may provide. This is so vital because some inherent information may not be detectable from the raw field data due to the attenuation of signal energy with depth and the influence of noise.

Principles of Ground Penetrating Radar
The basic principle behind the GPR is the principle of the scattering of electromagnetic waves. A short pulse of an ultra-high frequency electromagnetic (EM) wave within the range from 10 MHz to 2.6 GHz is propagated through the Earth. As the electromagnetic wave propagates through the ground, it encounters different Earth materials of varying dielectric contrasts. Part of the wave energy is reflected and part transmitted through the material due to the bulk changes in the material's electrical properties (e.g. relative permittivity (ε r ), magnetic permeability and electrical conductivity). The relative permittivity is a material property that controls the velocity of the EM wave propagating through material (Jol 2009).
The propagation velocity (v) is related to the speed of EM wave through vacuum (c) and ε r by formula: From equation (1), deductions can be made such that changes in subsurface material properties will cause a contrast in ε r which will affect the index of refraction by producing a reflected energy at the boundary between two materials. The relative permittivity is mainly controlled by the soil water content in low loss media. An increase in water saturation in a given soil will cause an increase in ε r thereby increasing the energy of the reflected EM wave (Huisman et al. 2003, Marcak et al. 2018). The propagation velocity v can also be calculated by: where d is the depth to a target object in the propagating media and t is the two-way travel time to the reflector (Manu et al. 2014).
The physical basis of GPR measurements lies in electromagnetic (EM) theory. Maxwell's equations mathematically describe the physics of EM fields, while constitutive relationships quantify material properties. Combining the two provides the foundation for quantitatively describing GPR signals. Details of this can be found in Daniels (2004).
Other important elements of the theory of GPR are electrical conductivity (σ) and magnetic permeability (μ). Electrical conductivity characterizes free charge movement (creating electric current) when an electric field is present. Resistance to charge flow leads to energy dissipation. Magnetic permeability describes how intrinsic atomic and molecular magnetic moments respond to a magnetic field (Jol 2009).

METHODS OF STUDY Field data measurement
The GPR field data were obtained from a farmland soil, composed of two parts: loamy topsoil and sandy material, within University of Agriculture campus in Balicka street, Krakow, Poland. Compaction was induced by the movement of tractors simulating farming activities. The tractor's tracks caused the compaction of the soil. The farmland is the practical agriculture demonstration plot of the University. The GPR equipment used for the data collection was the MALA Pro-Ex GPR system manufactured by ABEM MALA Inc., Sweden with a shielded antenna of central frequency 800 MHz (Fig. 1). The choice of antenna was based on target depth and resolution. Data collections were made in the common offset mode, with acquisition parameters as shown in Table 1. A total of 14 profiles were acquired. Three of the profiles have length of 70 meters while the remaining profiles, which are about 10 meters long, are perpendicular to the long three profiles (Fig. 1A, B).

A B
The GPR systems were setup and mounted, and then data were collected by pulling the wheeled antenna along the profile lines at walking speed. The system generates radar pulses at a given central frequency which in this situation 800 MHz, and send them into the earth through a transmitting antenna. The pulses are scattered back at electromagnetic discontinuities of the subsurface, mainly due to dielectric contrasts between horizons of different materials within the soil. The back scattered pulses were collected by a receiving antenna and the data presented as signal amplitudes versus travel time in form of image (radargrams).

Data processing
To improve the signal to noise ratio (SNR), the measured data were first scrutinized. This was made possible using the software REFLEXW (Sandmeier 2012). Subsequently low frequency components were removed from the data through the application of a "dewow" filter. To resolve all traces to a common zero point, a time zero correction was applied on all the GPR traces to bring them to a fixed starting time, below the ground wave. The background removal filter was further activated to remove coherent horizontal noise from the processed data. In order to enhance the signals received from deeper depths, the gain tool was applied to enhance the drastic fall in energy of the wave before getting to the receiver. Velocity for time-depth conversion of 0.094 m/ns was used base on WARR (Wide Angle Reflection and Refraction) measurements analysis.

Advance signal analysis
The concept of analysis of electromagnetic pulse energy (signals) may include the instantaneous parameters such as Instantaneous Amplitude (Envelope), Instantaneous Phase and Instantaneous Frequency. Matheney & Nowack (1995), Taner (2001) and Gao et al. (1997) have reported the concept of instantaneous parameters as widely used in electric engineering and geophysics. Instantaneous parameters such as the instantaneous amplitude, instantaneous phase, and instantaneous frequency are directly related to geometry and physical property variations of the medium through which the radar signal propagates. The classic approach for estimating instantaneous parameters relies on using the Hilbert Transform (HT). This approach extends a real signal to an analytical signal, by doing the HT for the real signal to get its imaginary counterpart, and extracts the instantaneous parameters by comparing the imaginary part and the real part of the analytical signal (Liu & Oristaglio 1998). Numerical manipulation of the raw data through the Hilbert Transform would not only enhance the delineated features in the signals but also reveal the inherent ones that may have been obliterated by weak energy due to the attenuation and noise from various sources. The following complex signal analysis operations were performed on the field data after the basic processing.

Instantaneous attributes (amplitude, phase and frequency)
Instantaneous attributes were computed sample by sample and represent instantaneous variations of various parameters. Instantaneous values of attributes such as trace envelope, its derivatives, frequency and phase may be determined from complex traces. Instantaneous amplitude outputs the envelope of the selected data volume at the sample location, it enhances, among others, lateral variations within events (Tomecka-Suchoń 2019). The instantaneous phase calculates the instantaneous phase at the sample location; it emphasizes the spatial continuity/discontinuity of reflections by providing a way for weak and strong events to appear with equal strength. Furthermore, instantaneous frequency provides outputs for the instantaneous frequency at the sample location. The instantaneous frequency attribute responds to both wave propagation effects and depositional characteristics, hence it is a physical attribute, which can be used as an effective discriminator (Taner 2001). Barnes (1991) obtained relations between instantaneous frequency and seismic signal attenuation. A given geophysical signal (e.g., seismic or GPR time traces) s(t) can be represented by its envelope a(t) and phase ϑ(t): and the quadrature (imaginary) trace is: and the complex (analytical) trace z(t) is then given by: The quadrature (imaginary) trace is the Hilbert transform of the real trace. It is obtained by phase shifting of the recorded trace by 90 degrees. The complex trace comprises of the real trace (recorded) and the imaginary trace (Barnes 2007). Once the quadrature trace is found, the instantaneous amplitude a(t) and the instantaneous phase ϑ(t) are: and respectively: The instantaneous frequency f(t) is the time derivative of the instantaneous phase: The instantaneous frequency f(t) can have large amplitude positive and negative spikes.
Large spikes in f(t) occur when the denominator of Eq. (8), which is the square of the instantaneous amplitude a(t), approaches zero more rapidly than the numerator. To reduce the large spikes in the instantaneous frequency, two more steps, damping and weighting, can be taken to get a much smoother f(t) (Barnes 1993 In addition to damping, the instantaneous frequency can also be weighted as: where the weight w(t′) can be taken as the squared instantaneous amplitude. The weighted instantaneous frequency approaches the average Fourier spectral frequency as the time interval T becomes large (Barnes 1993). Besides damping, weighting the instantaneous frequency will further stabilize the calculations. Tomecka-Suchoń & Marcak (2015) used computed instantaneous attributes of GPR data to identify the cause of sinkhole formation.

Other computed attributes
The following attributes were also calculated for the GPR data: Derivative of envelope. Time derivative of the instantaneous amplitude is the time rate of change of the envelope. It shows the variation of the energy of the reflected events. It is used to detect sharp interfaces and discontinuities. Mathematically, it is given by: where: a(t) -envelope, t -time.

Energy.
Energy output is a response attribute that returns the energy of a trace segment. This attribute calculates the squared sum of the sample values in the specified time-gate divided by the number of samples in the gate. The energy is a measure of reflectivity in the specified time-gate. The signal energy in the signal x(t) is given by: The higher is the energy, the higher is amplitude. This attribute enhances, among others, lateral variations within electromagnetic events and is, therefore, useful for object detection.

Event frequency.
It is an attribute that quantifies an events shape or distance relative to a next event and returns frequency properties.

Pseudo relief.
This attribute is applied on wave data in order to create a more consistent image for the easier interpretation of faults and horizons. A general formula for pseudo relief can be defined as: where: m -index of energy RMS function, k -shift.

Q-factor.
This attribute is computed from the instantaneous frequency divided by the bandwidth. The bandwidth is the absolute value of the envelope time derivative. So general formula for Q-factor can be defined as: The GPR amplitude envelope is a positive value and is a measure of the total signal energy. Changes in this parameter indicate changes in the distribution of dielectric permittivity and distribution of reflection coefficient which is related to dielectric permittivity (Akinsunmade et al. 2019a(Akinsunmade et al. , 2019b. Phase continuity is one of the main parameters used in the interpretation of GPR records. The correlation of phases of the signal between traces in GPR records can be used to eliminate natural amplitude decrease with depth. Instantaneous frequency is a parameter which indicates changes in mean signal frequency, possibly interpretable with changes in the lithology of the investigated strata (Tomecka-Suchoń & Marcak 2015).
Subsequent to the field data processing, the advance signal analyses were performed on the processed data using the complex trace analysis/spectral analysis extension of the REFLEXW software developed by Sandmeier (2012) and OpendTect software Earth sciences Inc. (dGBBeheer 2014).

RESULTS AND DISCUSSION
In this section the results of the different computed attributes of the processed GPR field data (Fig. 2) are presented and discussed. It is obvious that little information is only extractable from the processed field data as most salient information may be hidden. From among the analyzed profiles where measurements were made, the article presents the analysis results for one representative of the profiles (profile 7), which conveys the same responses in other profiles where measurements were made. The inverted arrows on the radargrams are the marker points of the tractor's tracks.
Analysis of the first vertical derivative attributes (Fig. 3A) shows a near surface layer (horizon) discontinuity at some parts of GPR section (1-9). These are as coincided with zones of impact of tractor movements at different passes. Moreover, two horizons are easily discernable at the near surface at depths of approximately 0.1 meter, and which are not connected with air or ground wave, thanks to their removal during processing. Also, variations in the subsurface layer properties are observable from the attribute of the signals.

Fig. 2. GPR process data for profile 7
The energy attribute highlights the material properties of the subsurface media. Zones with weak energy suggest a zone of low reflectivity and may be ascribed to the presence of conductive media such as clay or a highly saturated layer. Of particular interest is the delineation of the effect of the tractor passes within the underlying soil layer. This attribute has shown clearly depth extent (white circle zones - Fig. 3B) of the compaction zone, created by the impact of the tractor movements. The depth response to the compaction effects by tractor passes from distance 10 to 35 meters is approximately 0.15 meters while the depth from distance 40 to 60 meters is approximately 0.28 meters (vertical extents of zones of extremely low values of energy attribute beneath tractor passes tracks depict the compaction effect). This depth determination has been made possible by the energy signal attributes whereas, this is not discernable in the raw field data. The event frequency attribute (Fig. 3C) also corroborated the appropriate location of the zones of impact of compaction due to the tractor movement. This is shown in the low values of considered attribute (blue zone - Fig. 3C) at the near surface.
Essentially it emphasizes the extent of the topmost compacted zone of reduced porosity due to the compaction. Hence, penetration resistance (PR, which represents soil's strength) zone are easily detected within the soil layers. Figure 4A is the instantaneous amplitude attribute computed from the GPR signal and has revealed a discontinuity in the near surface horizons particularly the point of impact of tractor movement. The points are properly delineated and depth extents of the effect of the vertical stress is observable at depth range of 0.15 meters and 0.28 meters as shown in points A and B (Fig. 4A), due to changes in soil structures along profile. Horizons that are not easily discernable in the raw field data are better shown on this attribute. The distinct in the horizons are well shown in the instantaneous frequency attributes (Fig. 4B). This may be connected with the dielectric properties of the various subsurface materials. A continuous strong reflection at approximately 0.05-meter depth (shown with arrow Fig. 4B) which is extensive from the beginning of the profile to the end may suggest interface of highly saturated horizon. It is interesting to note however that this continuity was interrupted at the points of impact of tractor tracts as shown with inverted triangle in figures.
The computed instantaneous phase for profile 7 is as displayed in Figure 4C. This attribute gives an insight into the variations in horizon materials properties. This variation is revealed at zones designated with white circles (Fig. 4C). The encircled ones are thought to be made of soil materials different from the surrounding soils or represent disturbances within soil layer, which are less clearly observable in the field raw data. Distortions in the continuity of the horizons are more pronounced in the pseudo-relief attribute. Subsurface responses to the vertical impact of the tractor passes are better shown in this analysis supported by the use of an appropriate color scale (Fig. 5A).
The attribute Q factor is computed from the instantaneous frequency. It has shown clearly in Figure 5B that a fuzzy zone (a white outline polygon) that is thought to coincide with the attenuation of the electromagnetic signal energy which may be suggestive of clay zone, which is hard to distinguish from the processed data (Fig. 2).

CONCLUSIONS
The correlation of the computed signal attributes and basic processed images has revealed that more information about the subsurface is better discernable from the computed attributes than basic processed data. As a result, distinguishing layers of contrasting properties is more visible. Moreover, fuzzy zones which are not clearly observed in the processed radargram are properly depicted in the computed attributes. It is suspected that this region may be point of signal attenuation which may be due to presence of clay material. It can be inferred that complex analysis of the signal may have pointed to the inherent properties of the traversing media which may be vital information for further studies of the media.
Furthermore, the effects of the compacted zones due to the tractor tracks (movement) that were hazy in the processed field data are better enhanced and more clearly discernable in the resulted attributes of the data. The appropriate locations of these zones and possible characterization of such zone can be possible with diagnostic operations via attribute computation of the GPR electromagnetic signal.
The attempts made in this study have shown the possibility of using the complex signal analysis via attribute computation to delineate subtle characteristics of the subsurface media that are inherent but hidden in the GPR field data. Thus, the spatial variation of the properties of the signal traversing media is revealed, along with obliterated features. This may be due to the weak signal response (as it travels deeper in the subsurface) and, as such, are better enhanced by considered method. This study also confirmed that the analysis of GPR data is more effective when transformed into the frequency domain.
This type of analysis in the frequency domain can reduce the ambiguity of interpretation of Ground Penetrating Radar data.