Nuzzo.pdf 533 ANNALS OF GEOPHYSICS, VOL. 46, N. 3, June 2003 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques Luigia Nuzzo Osservatorio di Chimica, Fisica e Geologia Ambientali, Dipartimento di Scienza dei Materiali, Università di Lecce, Italy Abstract 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 Key words GPR – coherent noise – Radon Transform 1. Introduction The ability of Ground Penetrating Radar (GPR) to provide, almost in real-time, very im- pressive 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 ac- quisition 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 tow- ed 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 534 Luigia Nuzzo 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 diffi - cult 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 tra- jectories. 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. 2. Some surface-scattering case histories The first example is probably the most emblematic case of surface-scattering pheno- mena experienced by our research group. It was observed during a GPR survey carried out by a 35 MHz antenna near Nociglia (Lecce, Italy), shortly Fig. 1. Example 1: Nociglia, Lecce, Italy – Surface-scattering from trees. Location map of the radar profiles (13 and 37-38) and of the pine-trees (dots). 535 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques 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, Fig. 2a,b. Example 1: Nociglia, Lecce, Italy – Surface-scattering from trees. a) Radar section relative to profile 13 (fig. 1); b) the corresponding kinematic model. a b 536 Luigia Nuzzo complicated, intense anomaly was noted on the radar section (fig. 2a) that at first sight could suggest the presence of a large cavity. 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 a b Fig. 3a,b. Example 1: Nociglia, Lecce, Italy – Surface-scattering from trees. a) Radar section relative to profile 37-38 (fig. 1); b) the corresponding kinematic model. 537 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques 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 0 20 m 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). a b c 538 Luigia Nuzzo 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 re- flections 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, 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, 2002). This experience represented an invaluable tool for understanding the seriousness of surface- scattering phenomena and the importance of filtering procedures to attenuate them. 3. 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 hyper- bolic diffractions in zero-offset profiles, provid- ed the offset is replaced by the common source- receiver 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-scat- tering 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 539 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques 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 reinsert- ed in the original section (Nuzzo and Quarta, submitted). 4. 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. Isi- doro-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. 5a,b. Example 3: S. Isidoro-Torre Inserraglio street (Lecce, Italy) – Surface-scattering from a cable. a) Water- filled 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. a b 540 Luigia Nuzzo Fig. 6a-f. Example 3: surface-scattering removal by means of the linear Radon transform (t-p transform) on a hyperbolically windowed data matrix. Fig. 7a-f. Example 3: surface-scattering removal by means of the parabolic Radon transform (t-q). a e b c d f a e b c d f 541 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques 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 ma- trix with respect to the apex position (vertical line in fig. 7a), the two sub-matrices have been processed separately. For computational effi - ciency and to preserve the high-frequency con- tent 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 surface- scattering anomaly practically unrecognis- able 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/m2. 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-stretch- ing 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 542 Luigia Nuzzo Fig. 8a-f. Example 3: surface-scattering removal by means of the parabolic Radon transform (t'-q' ) in the t 2-stretched domain. 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). a e b c d f a e b c d f 543 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques 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). 5. 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 disclo- sed by some case histories derived from the Fig. 10a-d. Example 3: comparison of the results obtained for surface-scattering removal by means of time- domain methods in the hyperbolically windowed domain (background removal and subtraction of a running average computed over different number of traces, N ). a b c d 544 Luigia Nuzzo recent experience of GPR investigation for en- vironmental and archaeological applications. These examples highlight the need for a physi- cal and/or geometrical comprehension of the origin of the spurious signals in order to avoid interpretation pitfalls, and the need for mod- elling 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 pseudo- frequency) 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. Acknowledgements I wish to thank my Ph.D. tutors, Prof. Maria Teresa Carrozzo and Prof. Tatiana Quarta, for valuable suggestions and critical remarks during the doctoral activity where I developed the present work. I am grateful to my colleagues G. Leucci and S. Negri, and to the technicians S. Corriero, G. Fortuzzi and M. Luggeri for their valuable collaboration during data acquisi- tion. I also thank Prof. M.D. Sacchi from the University of Alberta, Canada, for making the Matlab codes on the Radon transform (for_taup.m and inv_taup.m) available on the Internet (at the website http://rubble.phys.ualberta.ca/%7Esacchi/ SEISMIC_LAB/) and the reviewers for the helpful suggestions. Appendix. The Radon transform. The Radon Transform (RT), originally introduced by Radon in 1917, is an integration procedure that maps a curve in the space (or time-space) domain to a point in the Radon domain. The RT and in particular the linear RT, has been widely used in geophysics, and especially in reflection seismic, for various purposes: plane wave decomposition, velocity analysis, multiple suppression and directional filtering (Stoffa et al., 1981; Schultz, 1982; Treitel et al., 1982; Yilmaz, 1987; Zhou and Greenhalgh, 1994; Pawlowski, 1997; Sykes and Das, 2000). The linear RT (t-p transform or slant-stack) is a particular 545 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques case of the generalized Radon transform (Beylkin, 1987), where the integration path is a line A straight line (A.1) having intercept t and slope p in the t-x plane maps to a point in the t-p domain. If t represents time and x space, p has the dimensions of the inverse of velocity and is usually named slowness (or ray- parameter). In the parabolic RT the integration path is described by (A.2) where q is the curvature of the parabola centred at zero offset. The parabolic expression (A.2) is the Taylor expansion, truncated at the second order, of the hyperbolic traveltime curve (A.3) Equation (A.3) could represent either a reflection in CMP or common-shot gathers, for a homogeneous medium with velocity v, or a diffraction in zero-offset gathers, providing v is an apparent velocity equals to half the true medium velocity. In other words, for a reflection (A.4) while for a diffraction . (A.5) Usually the parabolic RT has been applied on partially normal move-out corrected gathers, where the reflection events become approximately parabolic (centred at zero-offset). An alternative approach has been proposed by Yilmaz (1989), involving a t 2-stretching procedure. By setting t' = t 2 and t' = t 2, eq. (A.3) takes the form (A.6) which represent a parabola in the t'-q' plane. Using the same notations as before, it follows that for a reflection: (A.7) while for a diffraction . (A.8) In the t 2-stretched domain the curvature is independent from the intercept time. u p d x t px dx, , .τ τ( ) = = +( ) −∞ ∞ ∫ t px= +τ t qx= +τ 2 t x x= + ≅ +τ ν τ ν τ 2 2 2 2 21 2 . q = 1 2 2ν τ q = = 1 2 2 2 2ν τ ν τapp t x' ' + 1 ' + q' x 2 2= =τ ν τ2 q' = 2 1 ν q' = 1 appν ν 2 2 4 = 546 Luigia Nuzzo 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 d = Am (A.11) where d = d(w) and m = m(w) are the data and model vectors in the Fourier domain and A Alk= [ ] is a matrix of elements A i x qlk l N k= −( )exp .ω The (unconstrained) least-squares solution of (A.11) is (A.12) where AH is the Hermitian (complex conjugate transpose) of A and (AH A)-1 AH is the generalized inverse of A. To avoid singularities in the matrix AHA 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: AH 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 fmax and x-range, xr, good rules for selecting Dp and Dx are (Turner, 1990) (A.13) and the analogous ∆ x p fr≤ ( )1 / max , which leads to a constraint for the maximum allowable p-range, pr, or the maximum value of p, pmax, assuming a symmetrical range of slopes pr = 2 pmax (A.14) By analogy, Kabir and Verschuur (1995) derived similar formulas for the parabolic RT. For band- limited signals, with maximum frequency fmax and assuming that x ranges from 0 to xmax with sampling interval Dx, whereas the curvature q ranges from 0 to qmax with sampling interval Dq, these formulas are respectively (A.15) . (A.16) m q d x t qx N x , ,τ τ( ) = = +( )∑ d x t m q t qx N q , ,( ) = = −( )∑ τ m A A A d= ( )−H H1 ∆q x f≤ ( )1 2/ max max q x x fmax max max/≤ ( )1 2 ∆ ∆p x fr< ( )1 / max p f xmax max/ .≤ ( )1 2 ∆ 547 Coherent noise attenuation in GPR data by linear and parabolic Radon Transform techniques REFERENCES BANO, M., G. MARQUIS, B. NIVIÈR, J.C. MAURIN and M. CUSHING (2000): Investigating alluvial and tectonic features with ground-penetrating radar and analyzing diffraction patterns, J. Appl. Geophys., 43, 33-41. BEYLKIN, G. (1987): Discrete Radon Transform, IEEE Transactions on Acoustic, Speech and Signal Processing, 35, 162-172. CARROZZO M.T., V. BASILE, A. FEDERICO, S. NEGRI and L. NUZZO (1999): Indagini geofisiche nel centro del comune di Nardò, in Thalassia Salentina, Atti del I Incontro di Studi: Il Carsismo nell'Area Mediterra- nea, Castro Marina 1-2 Settembre 1997, 23, 219-230. CARROZZO, M.T., G. LEUCCI, S. NEGRI and L. NUZZO (2001): Preliminary GPR survey at Roman Ships archaeological site (Pisa, Italy), in Atti del XX Convegno Nazionale GNGTS, Roma 6-8 November 2001 (in press). CARROZZO, M.T., M. DELLE ROSE, A. FEDERICO, G. LEUCCI, V. MARRAS, S. NEGRI and L. NUZZO (2003a): Osservazioni geologiche e indagini geofisiche sul carsis- mo della costa Neretina, in Thalassia Salentina, Atti del II Incontro di Studi: Il Carsismo nell’Area Mediterranea, Castro Marina 14-16 Settembre 2001, 26 (in press). CARROZZO, M.T., G. LEUCCI, S. NEGRI and L. NUZZO (2003b): GPR survey to understand the stratigraphy of the Roman Ships archaeological site (Pisa, Italy), Archaeological Prospection, 10, 57-72. GRASMUECK, M. (1996): 3D ground-penetrating radar ap- plied to fracture imaging in gneiss, Geophysics, 61, 1050-1064. KABIR, M.M.N. and D.J. VERSCHUUR (1995): Restoration of missing offsets by parabolic Radon transform, Geophys. Prospect., 43, 347-368. LINES, L.R. and S. TREITEL (1984): Tutorial a review of least-squares inversion and its application to geophysical problems, Geophys. Prospect., 32, 159-186. NUZZO, L. (1997): Indagini georadar in zona interessata da dissesti idrogeologici (Nociglia), University of Lecce, M.Sc. Thesis (unpublished). NUZZO, L. (2001): Acquisition, processing and visualiza- tion of Ground Penetrating Radar data by means of classical and novel techniques – Some applications to archaeological and environmental studies, University of Messina, Ph.D. Thesis. NUZZO, L. (2002): Identification and removal of above-ground spurious signals in GPR archaeological prospecting, in European Geophysical Society Scientific Program, Nice, April 2002, EGS02-A-02468, CD-ROM. NUZZO, L and T. QUARTA: Improvement in GPR coherent noise attenuation using t-p and wavelet trasforms, Geophysics (submitted). PAWLOWSKI, R.S. (1997): Use of slant stack for geologic or geophysical map lineament analysis, Geophysics, 62, 1774-1778. SCHULTZ, P.S. (1982): A method for direct estimation of interval velocities, Geophysics, 47, 1657-1671. STOFFA, P.L., P. BUHL, J.B. DIEBOLD and F. WENZEL (1981): Direct mapping of seismic data to the domain of intercept time and ray parameter – A plane-wave decomposition, Geophysics, 46, 255-267. SUN, J. and R.A. YOUNG (1995): Recognizing surface scat- tering in ground-penetrating radar data, Geophysics, 60, 1378-1385. SYKES, M.P. and U.C. DAS (2000): Directional filtering for linear feature enhancement in geophysical maps, Geophysics, 65, 1758-1768. TREITEL, S., P.R. GUTOWSKI and D.E. WAGNER (1982): Plane-wave decomposition of seismograms, Geophysics, 47, 1375-1401. TURNER, G. (1990): Aliasing in the t-p transform and the removal of spatially aliased coherent noise, Geophysics, 55, 1496-1503. YILMAZ, O. (1987): Seismic Data Processing, edited by S.M. DOHERTY (SEG, Tulsa), pp. 526. YILMAZ, O. (1989): Velocity-stack processing, Geophys. Prospect., 37, 357-382. YOUNG, R.A. and J. SUN (1994): Recognition and removal of surface scattering in GPR data, in Proceedings of the Fifth International Conference on Ground Penetrating Radar, Kitchener, Ontario, 735-746. YOUNG, R.A. and J. SUN (1999): Revealing stratigraphy in ground-penetrating radar data using domain filtering, Geophysics, 64, 435-442. ZHOU, B. and S.A. GREENHALGH (1994): Linear and parabolic t-p transform revised, Geophysics, 59, 1133-1149.