Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques

Ground Penetrating Radar (GPR) is a geophysical method increasingly used in numerous shallow applications. Unfortunately, electronic or acquisition problems can cause the presence in the radargrams of coherent noise interfering with the useful signal. A commonly observed phenomenon, especially for not-shielded antennae, is the surface-scattering effect, due to reflection or diffraction from above-surface objects. These noise events appear with a characteristic hyperbolic moveout in the usual common-offset sections. Other frequent problems are related to the presence of horizontal or dipping features due to system-ringing or other non-geological causes. Several methods have been tried to overcome these problems, most of which involve time domain or Fourier domain filtering. This work presents an attempt to reduce some of these noise modes by an original adaptation of filtering techniques implemented in the Radon domain. The Radon Transform (RT), both in the linear (or t-p) and in the parabolic version (or t-q), has been widely used in seismic processing, especially for multiple removal, but is still quite unfamiliar to GPR practitioners. The results achieved by different RT based methods for coherent noise attenuation in a GPR field example, compared to those of more conventional techniques, are quite encouraging. Mailing address: Dr. Luigia Nuzzo, Osservatorio di Chimica, Fisica e Geologia Ambientali, Dipartimento di Scienza dei Materiali, Università di Lecce, Via per Arnesano, 73100 Lecce, Italy; e-mail: luigina.nuzzo@unile.it


Introduction
The ability of Ground Penetrating Radar (GPR) to provide, almost in real-time, very impressive images of the shallow subsurface makes this geophysical method increasingly used in geological, engineering, environmental and archaeological applications.Unfortunately, this method suffers from some electronic or acquisition problems that, in particular conditions, can cause the presence in the radargrams of coherent noise interfering with the useful signal.
A commonly observed phenomenon, especially for not-shielded antennae, is for example the surface-scattering effect.In fact, reflections or diffractions from above-surface objects, such as trees, fences, power lines, buildings, could be present in the radar sections and could obscure subsurface reflections, since, due to the much higher electromagnetic wave velocity in the air than in the ground, their arrival time could be in the time window of interest.
The problem of surface-scattering is not a marginal one.As pointed out by several authors (Young and Sun, 1994;Sun and Young, 1995;Bano et al., 2000), it can severely degrade the visibility of interesting reflections or even, if not recognised, can be the cause of interpretation pitfalls.The GPR antennae are generally towed directly above or very near to the ground surface.The radiation pattern of a dipole at the boundary of a dielectric halfspace has always a non-negligible lobe in the air, though generally considerably smaller than the main lobe oriented toward the ground.The back-lobe can be fairly energetic if the dipole antenna is placed above a lossy halfspace.Since shielding is more difficult to accomplish for low-frequency antennae, surface-scattering phenomena can be particularly troublesome for low-frequency GPR surveys.
Surface-scattering noise has generally a hyperbolic moveout on GPR constant-offset profiles, so that this problem closely resembles that of multiple attenuation on Common Shot or CDP gathers in seismic reflection.Sometimes, the relative position of the scattering sources with respect to the survey lines causes only the diffraction tails to be visible on the sections as (generally) steeply dipping lines.Linear noise features can be successfully removed by standard 2D Fourier transform (f-k) methods (Grasmueck, 1996) or the «domain filtering», i.e. a «local» f-k procedure (Young and Sun, 1999).Also methods based on the linear Radon Transform (RT), or t-p, or on the less known wavelet transform (Nuzzo, 2001;Nuzzo and Quarta, submitted) proved promising for linear noise feature attenuation, with performance comparable, if not superior, to that of time-domain methods.In the most frequent case of hyperbolic noise features, unless using particular tricks, approaches based on the linear RT change the problem of removing hyperbolic trajectories into that of removing elliptical trajectories.It is clear that the complexity of the problem could not be decreased in such way.More suited strategies, as in the seismic analogue, seem those based on the hyperbolic RT or, due to its computational suitability, on the parabolic RT or t-q transform.
In the following Sections, after a closer look at some case histories evidencing the seriousness of the surface-scattering problem, the early results of an attempt of coherent noise attenuation by t-q methods are presented and compared to those obtained by means of t-p and time-domain ones.

Some surface-scattering case histories
The first example is probably the most emblematic case of surface-scattering phenomena experienced by our research group.It was observed during a GPR survey carried out by a 35 MHz antenna near Nociglia (Lecce, Italy), shortly after the collapse of the roof of an unknown karstic cavity, only a few tens of meters away from a provincial road (Nuzzo, 1997).Because of the high density of collapse sinkholes and known karstic cavities or burrows in the neighbourhood, the presence of other unknown cavities in the uppermost 10-20 m deep layer was suspected, that could pose a serious hazard for human life.The early GPR survey concentrated mainly along and near the provincial road, oriented almost perpendicularly to the axis connecting the old and recent known sinkholes.At the intersection of the main road with a pathway bordered by seven secular pine-trees (dots in fig.1), a large, The suspicion that the anomaly could instead be due to surface-scattering from the nearby trees was originated by the observation that, after migration at the air velocity, the hyperbolas collapsed to two zones closely corresponding to the position of the two rows of trees, and was subsequently confirmed by the close agreement between the observed traveltimes and those estimated by a kinematic model of diffraction from the base of the foliage (fig.2b).Moreover, a different anomaly pattern was observed in the radar section relative to an almost orthogonal profile (fig.3a), whose general aspect was also confirmed by the corresponding kinematic model (fig.3b).In this situation, the surface-scattering problem affected zones of the radar sections free of interesting features, so that its removal was not needed, but only its recognition in order to avoid interpretation pitfalls.
Surface-scattering can be particularly evident also in urban environments where the presence of buildings can cause both inline and off-line reflections and diffractions.An exemplary case  occurred during a GPR survey carried out in the historical centre of Nardò (Lecce, Italy) to assess whether the cracks observed in the Cathedral apses and the damage affecting some historical buildings could be due to artificial causes or to local hydro-geologic conditions (Carrozzo et al., 1999).In the radar section presented (fig.4b), which refers to the central profile (heavy line in fig.4a) acquired in Piazza Pio XI by a 100 MHz antenna, clear diffraction hyperbolas at the air velocity are easily identifiable: the horizontal positions of their apexes match well those of the buildings at the opposite sides of the square.On the other hand, numerous events asymptotically approaching the direct air wave, and crossing the direct and reflected ground waves, are also visible on the CDP gather (fig.4c) acquired moving the two 100 MHz antenna along the same line symmetrically with respect to the position marked by a dot in fig.4a.Being not flat, these events are unlikely to be caused by inline objects; rather, being well fitted by reflection hyperbolas at the air velocity, they appear due to off-line reflections from buildings.Also in this case, it was not necessary to remove these spurious signals, since they did not interfere heavily with the interesting features.However, in other cases, surface-scattering problems can reduce the visibility of useful reflections so dramatically as to make necessary both the careful recognition of the spatially correlated noise modes and the development of proper approaches for their reduction.This was, for example, the situation during a GPR survey carried out in the archaeological site of the Roman Ships near the S. Rossore station in Pisa, Italy (Carrozzo et al., 2001(Carrozzo et al., , 2003b)).The investigation, having mainly stratigraphical purposes, was performed with a 35 MHz antenna on the plan of the excavation, located about 5 m below the ground level and protected on every side by a 6 m high iron sheet-piling.Surface-scattering from the highly-reflective iron enclosure and from other metallic objects near or around the partly recovered boats, was so severe as to require the development of specific processing strategies to help the recognition and removal of the coherent noise in order to uncover useful reflections (Nuzzo, 2001(Nuzzo, , 2002)).This experience represented an invaluable tool for understanding the seriousness of surfacescattering phenomena and the importance of filtering procedures to attenuate them.

Radon based filtering techniques for hyperbolic coherent noise removal
The parabolic RT (Appendix), or t-q, has been intensely employed in seismic processing for attenuating multiple reflections, showing, as well known, a hyperbolic moveout on Common Shot or CDP gathers.Since the t-q transform maps hyperbolas into pieces of hyperbolas, this transform can be a useful (and more efficient) approximation of the hyperbolic RT only in the small-offset approximation, where the hyperbola is fairly well approximated by a parabola, or if applied after modifying the data (e.g., after a partial Normal-MoveOut (NMO) correction) so that hyperbolas become more similar to parabolas.This has been, indeed, the way most frequently followed in seismic processing: in this way an almost parabolic trajectory is mapped into a point-like zone in the t-q domain, where it can be more easily removed.
The same argument can be valid for hyperbolic diffractions in zero-offset profiles, provided the offset is replaced by the common sourcereceiver abscissa relatively to the apex position and the velocity is replaced by an apparent velocity equal to half the true medium velocity (Appendix).A straightforward extension of this quite diffuse seismic procedure could therefore be envisaged for hyperbolic coherent noise attenuation in zero-offset (or monostatic) GPR profiles, and especially for surface-scattering removal, where the medium velocity is that of the electromagnetic waves in the air (0.30 m/ns).Unlike the most common seismic situation, a full hyperbola (instead of a half one) is usually present on GPR profiles; moreover, the function inside the summation in the discrete t-q is an even function, so that this transform cannot distinguish contributions coming from one or the other side of the hyperbola.To correctly handle the problem, it is necessary to split the data matrix into two sub-matrices with respect to the apex position, and to process them separately.If applied on the uncorrected data, the procedure is expected not to work perfectly, because of the already outlined departure of the diffraction path from the parabolic one at larger distances.For this reason, the t-q filtering can be applied after correcting the data so that hyperbolas become more similar to parabolas, such as by means of a t 2 -stretching procedure (eq.(A.6)).
An alternative approach can take advantage of the linear RT.Conceptually, the surface-scattering hyperbola after being flattened by means of a pseudo-NMO correction, maps to a point at zero slope in the t-p domain, where it can be easily muted, and therefore, after undoing the NMO, it would be removed from the radar section.To save computer time and to avoid distortions (as the stretching at large distances from the apex), an extraction procedure can substitute the pseudo-NMO correction: a sub-matrix can be extracted from the data in a hyperbolic window, tailored tightly along the air diffraction (using the air velocity and having the height of the estimated radar wavelet), transformed and filtered in the t-p domain, inverse transformed and reinserted in the original section (Nuzzo and Quarta, submitted).

Field example results
The radar section used to test the various Radon based approaches is a 510 × 460 matrix extracted from a profile carried out to assess the presence of underground cavities or fractured zones along a piece of the coastal street S. Isidoro-Torre Inserraglio (Lecce, Italy), where a karstic cavity collapsed in 1992 (fig. 5a) (Carrozzo et al., 2003a).The profile, carried out with a 200 MHz antenna, passed above the refilled collapse sinkhole (intense diffractions at the leftmost side in fig.5b) and below an aerial telephone cable (diffraction hyperbola at the air velocity centred at the abscissa 512 m).The fact that the hyperbola fades away far from its apex and the amplitude variation between the left and right branches make this real case a near ideal one for testing the various methodologies.
The radar section has already been used to test other filtering methods (time-domain, wavelet and t-p).In each case, only a hyperbolically windowed data matrix has been filtered and subsequently reinserted in the original section (Nuzzo and Quarta, submitted).As an example,  fig.6a-f outlines the filtering process using the linear RT: after extracting a sub-matrix in a hyperbolic window (fig.6a,b) and transforming it to the t-p domain (fig.6c), only the data near p = 0 ns/m have been preserved (fig.6d) and inverse transformed to model the air diffraction (fig.6e), which has been subsequently subtracted from the original to yield the residual noise-free data (fig.6f).
For the present work, using the parabolic RT, the full matrix has been transformed, either for the original or the t 2 -stretched domain.As previously outlined, after splitting the data matrix with respect to the apex position (vertical line in fig.7a), the two sub-matrices have been processed separately.For computational efficiency and to preserve the high-frequency content of the data, instead of a muting process a modelling and subtraction procedure has been developed, as previously carried out for the t-p process.In the t-q process, to reduce aliasing problems the maximum allowable sampling interval and range in the q parameter (curvature), have been selected compatibly with the empirical formulas (A.15) and (A.16) derived by Kabir and Verschuur (1995).Both quantities depend on the maximum frequency of the band-limited signal.In the present situation, the dominant band was between 0.1 and 1 GHz.After setting the above parameters in the forward parabolic RT, the t-q panels shown in fig.7b,c were obtained: not surprisingly, the hyperbola did not map to a point-like region; moreover some diagonal and horizontal features caused by the more intense anomalies disturbed the whole t-q sections.The weaker amplitude of the right branch of the hyperbola made the surfacescattering anomaly practically unrecognisable in the corresponding t-q panel.However, knowing the time location of the apex (t ) and having estimated the air diffraction velocity (v = 0.296 m/ns), the curvature at the apex could be predicted using formula (A.5), valid in the small-offset approximation.In this case, it yields a value of about 0.5 ns/m 2 .On the other hand, the filtering window choice can be guided by the slightly better visible anomaly in the left panel.Selecting the same t-q window for both panels and muting all the external zone (fig.7e,f), inverse transforming and recombining the two sub-matrices yields the result in fig.7d.Despite a minor edge effect at the junction of the two matrices, and some rather weak artifacts at greater distances from the apex, an acceptable reconstruction of the diffraction hyperbola is achieved.This is, however, allowed by the lucky circumstance that the hyperbola fades away far from its apex, so that only a small portion near the vertex, where the parabolic approximation holds, effectively contributes to the summation.
As a further test, instead of applying a (pseudo) partial NMO correction, the alternative approach proposed by Yilmaz (1989) involving a t 2 -stretching procedure was followed.By squaring the time axes (t' = t 2 and t' = t 2 ), the hyperbola in the t'-x plane becomes a parabola in the t'-x domain, so that the so far exposed filtering technique can be applied in a similar way in the time-stretched domains (t'-x and t'-q').This non-linear mapping, while stretching later events, compresses the earlier ones, so it is hazardous for filtering shallow events.Moreover, it needs to resample the data at a constant Dt' increment, generally taken as the total stretched time divided by the original sample number.A finer sampling, and thus a higher computational cost, could be useful to avoid shallow event distortion.On the other hand, in the stretched domain the curvature, q', is independent from the intercept time, t ', of the diffraction parabola apex (eq.(A.8)), and equal to around 45 (ns/m) 2 for the air velocity.After estimating (by a 1D Fourier transform) the useful (pseudo) frequency range, and after setting a useful sampling interval in q' (for instance in base of eq.(A.15)), the discrete parabolic RT can be performed over a range of curvatures generally much narrower than the maximum one allowed by eq.(A.16), provided it is sufficiently greater than the just obtained 45 (ns/m) 2 value.In fig.8a-f the results of the application of the parabolic RT based approach in the t 2 -stretched domain are shown: apart from a better visibility of the diffraction anomaly through the slant features in the t'-q' panels, an almost comparable quality in the reconstructed data as in the previous test is obtained.
This can be better appreciated from fig. 9a-f, where the results of the three Radon based techniques so far applied (linear RT in the hyperbolically windowed domain, parabolic RT Fig. 8a-f.Example 3: surface-scattering removal by means of the parabolic Radon transform (t'-q' ) in the t 2 -stretched domain.in the original and in the time-stretched domain) are compared.The t-q filtered section (fig.9c) is obtained by subtracting from the original section (fig.9a) the t-q reconstructed hyperbola (fig.9b), whereas to yield the t'-q' filtered one (fig.9f ) the reconstructed hyperbola (fig.9e), obtained after undoing the t 2 -stretching on the t'-q' reconstructed parabola (fig.8d), has been subtracted from the original section.The results of the two t-q based techniques are comparable to each other and to those previously obtained by means of the modified t-p approach (fig.9d).Moreover, they are also quite similar to those produced by more conventional time-domain methods applied in the hyperbolically windowed domain, such as those involving the subtraction of a running average computed in a horizontally sliding window, provided a suitable horizontal size (expressed in number of traces, N) is chosen (fig.10a-d).In the present case, the conventional techniques worked well for N varying between about 40 and about 80 (fig.10c), whereas they failed for smaller or larger values (fig.10d and b) or when subtracting the average trace computed in the whole (extracted) section (fig.10a).

Conclusions
The gravity of the surface-scattering problem, appearing as coherent noise with either linear or, more frequently, with hyperbolic moveout on the zero-offset radargrams, has been disclosed by some case histories derived from the These examples highlight the need for a physical and/or geometrical comprehension of the origin of the spurious signals in order to avoid interpretation pitfalls, and the need for modelling and processing strategies as an aid for their understanding and, eventually, for their removal.
The Radon transform could be a useful tool for coherent noise reduction in GPR data.Filtering procedures based on either the linear or the parabolic discrete RT can produce results comparable to other classical methods, but allow more flexibility.On the other hand, the effectiveness of Radon based methods depends on the parameters selected for the forward and the inverse transformation and the accuracy of the selection of the pass-or reject-zones.In the example proposed, the useful frequency (or pseudofrequency) band, and in particular the upper limit upon which the choice of the sampling and range of q in the Radon domain depends, was selected as narrow as possible, on the basis of the average 1D Fourier spectrum of the whole section, to limit the computational cost.Moreover, as usual, the damping factor was chosen as a compromise between stability and resolution: studying the trend of the RMS difference between the original and reconstructed section (without filtering) versus the damping factor (in the range 10 -5 -10 -1 ), a value of 10 -3 was selected, corresponding in all cases to the beginning of a steeper increment in the data misfit.More problematic is the selection of the zones to be muted (or preserved) in the Radon domain: in the present situation it was carried out heuristically by visual inspection of the filter results after various trials.Through a real data example it has been demonstrated that surface-scattering can be attenuated by means of either a modification of the t-p filtering technique or by t-q methods: the first allows a reduced data matrix to be processed, but needs an accurate estimation of both the coordinates of the apex of the hyperbola and the propagation velocity; the second (t-q) and the third (t-q with t 2 -stretching) need to split the data matrix in two sub-matrices for a proper amplitude handling of the left and right branch of the hyperbola.The last method seems the most effective from a computational point of view, since in the stretched domain the q-range can be restricted to a zone around the (fixed) curvature corresponding to the air velocity.The independence of the curvature from the time of the hyperbola apex could be exploited to filter all at once diffraction hyperbolas whose apexes have the same abscissa but different ordinate, as for the case of diffraction from buildings or vertical (metallic) fences.
For digital signals the discrete forward (DRT) and inverse Radon transform (IDRT) are defined respectively as (A.9) (A.10) where d(x,t) and m(q,t) represent data in time-space and in the Radon domain respectively, and N = 1 for the linear RT while N = 2 for the parabolic RT.Indeed, from an amplitude point of view a correction term is required in eq.(A.10).Computationally, it is easier to implement the transform in the temporal frequency domain and to obtain the forward transform as a least-squares solution of the inverse RT.In this way, for each frequency (or angular frequency w), a set of linear equations has to be solved that in matrix notation can be written as where A H is the Hermitian (complex conjugate transpose) of A and (A H A) -1 A H is the generalized inverse of A. To avoid singularities in the matrix A H A and thus numerical instability, a new matrix is built by adding a small positive number g, known as the regularization parameter or damping factor, to the diagonal elements: A H A + g I, where I is the identity matrix (Lines and Treitel, 1984).In general g is a function of w, and its choice depends on the noise content of the data.Finally, an inverse FFT is performed to obtain m (q,t).
The choice of the discretization parameters is very important to avoid (or at least to limit) aliasing problems.Whereas the sampling in t normally is the same as that in t, more crucial is the choice of the sampling in p (or q) as well as in x (Dp or Dq and Dx).
For the linear RT, in case of band-limited signals, with maximum frequency f max and x-range, x r , good rules for selecting Dp and Dx are (Turner, 1990) (A.13) and the analogous ∆x p f r ≤ ( ) 1 / max , which leads to a constraint for the maximum allowable p-range, p r , or the maximum value of p, p max , assuming a symmetrical range of slopes p r = 2p max (A.14)By analogy, Kabir and Verschuur (1995) derived similar formulas for the parabolic RT.For bandlimited signals, with maximum frequency f max and assuming that x ranges from 0 to x max with sampling interval Dx, whereas the curvature q ranges from 0 to q max with sampling interval Dq, these formulas are respectively
Fig. 4a-c.Example 2: Nardò, Lecce, Italy -Surface-scattering from buildings.a) Location map of the GPR survey; b) radar section relative to the profile evidenced by a heavy line in (a); c) CDP gather centred at the position marked by a dot in (a).

Fig
Fig. 5a,b.Example 3: S. Isidoro-Torre Inserraglio street (Lecce, Italy) -Surface-scattering from a cable.a) Waterfilled karstic cavity collapsed in 1992.b) Radar section showing a diffraction hyperbola from an aerial telephone cable crossing the street near the debris-refilled collapse sinkhole.

Fig. 9a -
Fig. 9a-f.Example 3: comparison of the results obtained for surface-scattering removal by means of different Radon based methods (linear RT in the hyperbolically windowed domain, parabolic RT in the original and in the time-stretched domain).

Fig
Fig. 10a-d.Example 3: comparison of the results obtained for surface-scattering removal by means of timedomain methods in the hyperbolically windowed domain (background removal and subtraction of a running average computed over different number of traces, N ).
d = Am (A.11) where d = d(w) and m = m(w) are the data and model vectors in the Fourier domain and A A lk