ANN-based sub-surface monitoring technique exploiting electromagnetic features extracted by GPR signals

In this work we consider the problem of determining the dielectric characteristics of sub-surface layers by means of GPR systems. In particular, a suitable electromagnetic feature (theR0 parameter), strictly related to the geophysical parameters of the scenario, is first extracted from the GPR e.m. signal and then fed to an artificial neural network (ANN) in order to derive the dielectric permittivity of the sub-surface layer.


Introduction
The ground penetrating radar (GPR) is a probing system designed primarily for the detection of objects and/or interfaces below the earth's surface (Daniels et al., 1988), (Carin, 2001). The simplest architecture consists of two separated broadband antennas, typically "bow-tie" dipoles: a transmitter (TX) for the generation of microwave short-pulses and a receiver (RX) which measures the backscatter radiation from buried targets.
Due to its non-destructive nature and proven diagnostic effectiveness, GPR-based solutions have been developed in many application fields, from infrastructure analysis to pavement profiling, from railroad deterioration assessment to environmental monitoring. In particular, a very challenging issue is the exploitation of GPR as a remote sensing tool for detecting buried objects or characterizing the subsurface structure and properties (Vellidis et al., 1990;Mellet, 1995;Golovko, 2004;Soldovieri et al., 2007).
In our past research (Caorsi and Cevini, 2005b), embracing the ANN approach, the evaluation of layers permittivity was addressed focusing on the analysis of a numerical simulation of the electromagnetic field scattered by the buried layer, obtained as the difference between the whole signal at the receiver and the field emitted by transmitter in freespace. Unfortunately, this kind of approach -although providing good and encouraging results -is too dependent on the source properties (for example, current/voltage pulse that feeds the antenna) and, in some cases, needs complex NN architectures.
In order to reduce as much as possible the a priori information about the electromagnetic source and to pursue a strong generalization, a new approach is introduced which manages only the amplitude of the received waveforms and exploits the fixed time scale, forced by scenario and antennas geometrical configuration, to determine the signal time dynamics.
In the next section, the algorithm is shown, first describing the theoretical bases and then indicating the rationale behind the actual implementation. Detailed results and error assessments are presented in Sect. 3. Conclusions hint final remarks and future developments.  layer, another homogeneous medium of B and infinite thickness is placed.
The GPR consists of two antennas (one TX and one RX) separated by d 0 , and the whole device is positioned at distance x 0 above the ground. In general, the signal S(t) sensed at the receiver can be decomposed into the contribution of the direct link (S d (t)) and the one due to the discontinuity effect of the air-ground interface (S r (t)).
The ratio is a feature strictly related to propagation conditions and to dielectric characteristics of the layer. Measured nearby the discontinuity plane and in case of normal incidence, it would be the theoretical reflection coefficient ; therefore, given a fixed geometric configuration of the GPR system -i.e. x o and d 0 -we can avoid the dependence on propagation and oblique incidence effects (as they are constant), so that: If this parameter can be extracted from the overall signal at RX, it is possible to determine L by inverting Eq. 2. In the ideal case, S r (t) and S d (t) are completely separated, and R can be simply computed as: where t max corresponds to the maximum of the considered signal (r or d) and the absolute value has been introduced to prevent from phase inversions that can happen at the interface. Unfortunately, the typical working conditions (distances and emitted pulse length) do not allow to keep the direct and the reflected signals separated: usually the echo 0 0.5 1 1.5 2 2.5 x 10 reaches the receiver when the tail of the direct pulse is still present. It is straightforward that Eq. 3 loses its meaning and the maxima-seeking is now impracticable. Anyway, the fixed geometry enables the detection of two further instants that can be chosen for a new, but consistent, characterization of R : where t fp is the time instant corresponding to the first peak of the signal we are considering (r or d). This instant can be easily evaluated within the direct signal: assuming that the layer is non-dispersive, we expect the echo to show the same shape and time length of the original waveform. This means that, also within the reflected signal, the first peak occurs after the same time shift found for the direct signal. Such hypothesis, although strong in certain cases, is absolutely reasonable for the considered problem and in many GPR applications, since it is verified for most dielectric materials within a frequency range up to a few GHz (Pozar, 1990). The keypoint of the new formulation is that R is still univocally bound to L , since the tail of the direct signal can be simply viewed as a constant additive term that does not depend on the layer's permittivity: and consequently it can be used within the inversion process to retrieve the layer's permittivity. Actually, since this procedure is intended to face the most general cases, without taking into account the pulse shape and length, from the operational viewpoint t fpd is not chosen as the instant when the first signal peak occurs, but as the instant when the signal reaches its maximum value within Therefore, the parameter can be rewritten as: As can be easily seen, the last formalization in Eq. 6 comprises all the previous ones: depending on the considered case, t mw coincides with t max or t fp . In fact, in the former case (distinct signals), the theoretical arrival of the echo would arise after the direct wave complete attenuation and the "search" window would hold the absolute signal maximum. In the latter case (overlapping signals), depending on the geometry, the echo would start in any point comprised between the start and the end of the direct pulse: if at least one peak (a local maximum) is detectable, then our algorithm would choose the corresponding time position, otherwise the ending instant of the 'search' window.
It is interesting to notice that, within the whole framework, the role of the layer's thickness x L is marginal. In fact, in the worst case (thin layer and low permittivity) the first echo could be covered by the reflected wave generated at layerground discontinuity. Nevertheless, if this overlap does not alter S(t mw r ), the relationship between R and L continues to be univocal.
Once the parameter R has been extracted, a mapping function for recovering the layer's permittivity is needed. The proposed approach does not perform an analytical inversion process, but exploits the potentialities of an artificial neural network (Haykin, 1994). More in detail, a Multi-Layer Perceptron (MLP) has been employed, made up of three layers (input-hidden-output) of neurons, which are ac- tivated according to a bipolar sigmoid function. Furthermore, MLPs are supervised neural networks and they require a training phase during which a set of "examples" (input data and the corresponding expected outputs, the target data) is employed for the mapping process. In particular, internal weights and biases, are updated, according to a modified version of the Back Propagation (BP) algorithm, by a recursive minimization of the error between the reconstructed outputs and the expected values, computed as: where T is the number of training examples, N L the number of neurons and x and x e are the actual and the expected values of the output, respectively. The typical trade-off in supervised approaches resides in the higher accuracy with respect to the dimensions of the training set, whose creation is a very time consuming and sometimes difficult operation. For our actual case, this onerous task of collecting "on field" samples can be avoided by a numerical simulation of the electromagnetic field. We will next show that the proposed approach enables reducing this significant dependency from the training set: having only one input parameter and one output feature, the ANN architecture can be very simple, and the number of needed training samples very small. As can be seen, this also allows an easier retrieval of experimental training data, the typical bottleneck of such procedures.

Numerical results
In this paragraph the performances of the proposed method are discussed. In order to carry out a wide assessment of its potentialities, the study scenario ( Fig. 1) is simulated through a numerical computation based on the Finite Element Method (FEM). In fact, this allows us to fully control the characteristics of the system (antennas, signal, layer), and thus to vary them for a more detailed analysis. In particular, for the GPR system, the TX is modeled as a current line fed with a gaussian-modulated sinusoidal pulse (center frequency of 500 MHz) of length 3.5 ns, while the design of the RX is not taken into account, since the re- ceived signal is simply the electromagnetic field calculated at the position of the receiving antenna (at distance d 0 =30 cm).
Regarding the medium properties, the discontinuity layer is characterized by L values ranging from 2 to 20, with 0.5 step, for a total dataset consisting of 37 samples. As regards the layer's thickness, whose effect on the R parameter has been discussed in the previous section, we opted, for the sake of simplicity, for a x L of 90 cm, so that, in none of the cases, the second echo covers the first one. Actually, for the described scenario, it is an excessive value, since even from x L ≥20 cm (in the worst case of L =2) the second reflection starts arriving after the first half of the first echo, without altering the peaks of interest. From the whole dataset, a part of the "input parameters/output features" (R / L ) pairs is necessary for the training of the neural network, while the rest is used for the test phase (the R parameters will be fed to the ANN and the outputs compared with the actual permittivities).
The first part of the discussion is right intended to show the performances of the algorithm depending on the number of training samples. Then, being the simulated signals in an ideal scenario free of disturbances (except for numerical fluctuations), the robustness of the method is assessed by modelling the most common noise sources that affect GPR systems, such as additive and multiplicative noise. Additive noise is usually due to electronic devices (we assume it is a white noise with gaussian probability density function, AWGN), while the so-called fading is the multiplicative noise that characterizes propagation effects, for instance multiple reflections. Therefore, the signals' waveforms are altered by adding and/or multiplying random-generated sequences of increasing powers (SNR from 10 to 50 dB), according to: S n (t) = (1 + n m ) · S(t) + n a (8) As a final remark, the analysis is both carried out for separate and overlapping signals, to show the high degree of generalizability of the R parameter. For simulating these two different cases, it is either possible to tune the geometrical parameters, d 0 and x 0 , or the pulse duration. Nevertheless, since the GPR characteristics usually depend on system design, the only aspect users can manage is the antenna-ground distance. Consequently, we decided to keep the GPR's model unchanged (d 0 = 30 cm and pulse length = 3.5 ns) and to choose two different antenna-ground distances: x 0 = 60 cm (separate) and x 0 = 30 cm (overlapping). Fig. 2 depicts the total received signal when S r (t) and S d (t) are separated. From the fixed geometry, we can derive: · S d (t) start = 1 ns · S r (t) start = 4.12 ns ⇒ "search window" ∈ [1; 4.12] ns Therefore, in this configuration, the maximum value reached within the "search" window is the global maximum of S d (t), which corresponds to t mw d = 2.04 ns. The time shift t from the signal start is 1.04 ns and thus t mw r is 5.16 ns.

Separate signals
Once the R parameter is evaluated, it is associated to the related permittivity and then fed to the neural network for the training phase. The graph in Fig. 3 clearly shows that performance stability is reached with very small training set dimensions (even with 4 samples mean errors are below 5%).
Of course, the greater the training set, the higher the accuracy, but differences are not significant. The best trade-off in terms of accuracy and needed ground-truth references can be an 8-sample training set, whose error analysis is presented in Fig. 4. It can be noticed that the mean relative error is only 1.4%, which is exceeded in very few cases, for a maximum error of 6%. Moreover, 8 is definitely a small number, so that it is possible to easily create an experimental dataset of references. Otherwise the only solution would be the numerical simulation of the environment, which is sometimes a difficult task.
For the robustness analysis, the signal is corrupted with the additive and multiplicative terms n a and n m of Eq. 8. A simple practical device for facing disturbances consists in selecting not only a single amplitude value but performing a mean over a small neighborhood of that point.
The three curves in Fig tiplicative (A+M). As can be seen, for low SNRs the contribution of additive noise is more sensible than multiplicative one; nevertheless, for SNRs higher than 30 dB, the effect is almost identical. The joint effect of AWGN and fading is higher than the separate contributions between 25 and 40 dB; above results are comparable, as noise effect is not still significant. By the way, it must be stressed that the error absolute values are certainly promising, since the maximum mean error is only 12% (in case of a 10 dB SNR, which alters the original signals so that they are almost undistinguishable, see Fig. 6) and it rapidly decreases to 4%.

Overlapping signals
In Fig. 7 is shown the total received signal when S r (t) and S d (t) overlap. From the fixed geometry, it comes that: · S d (t) start = 1 ns · S r (t) start = 2.24 ns ⇒ "search window" ∈ [1; 2.24] ns The "search" window is of course narrower, but the maximum signal value that can be resolved is still the global maximum. It is not important that the echo is mixed up with the tail of the first signal, because we infer the position of its maximum through the a priori knowledge of the system geometry: the time shift is still 1.04 ns, thus t mw r = 3.28 ns. Following the same outline of the previous subsection, the error trend due to the number of training samples is studied. Again, with 2 and 3 training inputs, errors are very high, but as the set increases they rapidly decrease (Fig. 8). With respect to the values found for separate signals, we have lower errors, even though the presence of the direct pulse's tail. This is simply due to the fact that the distance covered by the echo is shorter: the power of the received signal is greater and the amplitude differences due to diverse L are more sensible. Anyway, the tendency remains similar, thus also in this case we decided to focus on the 8-sample training set. As reported in fig. 9, the mean error value is 0.94%, while the maximum is 4.65%; in-between there are only 3 reconstructed values around 3% and 3 values around 1-1.5%.
The robustness analysis (Fig. 10) strengthens the previous results, since we have the maximum error of 10% with a SNR of 10 dB, that steeply slopes down for SNRs higher than 20 dB. Nevertheless, differently from the other configuration, here additive noise's effect becomes similar to multiplicative a bit earlier, around 25 dB, whereas the combination of the two contributions is always greater than the single ones, at least until the SNR nullifies any interferences.

Conclusions
In this paper, a new method for the reconstruction of subsurface layers' permittivity by means of GPR signals has been presented. The problem is solved by extracting a simple, but strong electromagnetic parameter which is then associated to permittivity through a neural network. The aim of this approach is to reduce as much as possible the a priori information required for parameter extraction and inversion processes, as well as to minimize computational and time costs.
The performance assessment, carried out over a large amount of scenario conditions (acquisition geometry, layer's permittivity, noise), shows very high accuracies and generalization capabilities. Moreover, the simple neural architecture employed for the inversion enables a very fast computation, which makes the method suitable for real-time applications. Besides a more detailed analysis about the influence of layer's thickness on the R parameter, next developments will address the implementation of two in-series networks for recovering the thickness itself, as well as the reconstruction of multi-layered surfaces. Finally, the method will be tested on pulses within the range of super high frequencies, where materials become dispersive and the use of non-monochromatic source is critical.
Edited by: F. Soldovieri Reviewed by: two anonymous referees