Application of 3 D visualization techniques in the analysis of GPR data for archaeology

In this work, some results of a GPR survey carried out in a 10 000 m large archaeological site, located in Lecce (Italy) near to a necropolis dating from the Messapian to the Roman imperial age, are reported. After a preliminary survey, performed on the entire area along parallel 1 m spaced profiles using a 200 MHz and a 500 MHz antenna in single-fold continuous mode, some smaller areas were selected, where the survey was repeated decreasing the profile spacing down to 0.50 m for the lower frequency antenna and to 0.25 m for the higher one. For two selected zones (D and B) the processed data were visualized in 3D space not only by the standard time slice technique, but also by two recently proposed approaches, namely by iso-amplitude surfaces of the complex trace amplitude and by 3D projection of energy and envelope stacks. The immediacy in revealing the spatial positioning of highly reflecting bodies, such as the anomaly interpreted as an old refilled cistern in zone D, makes 3D visualization techniques very attractive in archaeological applications of GPR. Their sensitivity to the signal/noise ratio is, on the other hand, highlighted by the quite poor performance in zone B, where the only reliable result provided by all the techniques was the soil/bedrock reflection, whereas none of them could effectively enhance the visibility of weak dipping reflections noted on 2D sections and probably related to fractures or bedding planes in the calcarenitic basement. The performance of the various techniques in these two different situations allowed insights into their main advantages and drawbacks to be gained.


Introduction
Ground-Penetrating Radar (GPR) is a fast and cost-effective electromagnetic (EM) method which, in favourable conditions, i.e. mainly resistive non-magnetic environments, can provide valuable information on the shallow subsurface.Since it is based on the propagation and reflection of EM waves, it is sensitive to variations of the EM parameters in the subsoil, especially the dielectric constant and the electric conductivity (Davis and Annan, 1989).Despite its relatively low penetration depth (especially with high-frequency antennae and in moderately conductive environments), the GPR resolution capability (also depending on frequency and soil properties), by far greater than that obtained by other geophysical methods, makes this technique suitable for high-resolution shallow studies like archaeological applications and shallow stratigraphy mapping.
The increasing need for detailed 3D imaging of the shallow subsurface makes 3D GPR one of the most important current topics.Although the significant advantages of 3D georadar surveying are well documented, especially for mapping geological features (Grasmueck, 1996;Grandjean and Gourry, 1996;Sigurdsson and Overgaard, 1998), the higher horizontal and vertical resolution required in archaeological applications makes 3D GPR more expensive for large-area surveys.Because of the small size of common archaeological targets, for a proper 3D acquisition a submeter line separation is generally needed, but in most cases 1 m spacing is used.Subsequently data are processed by means of 1D and 2D techniques, since the crossline spacing is still too large to obtain remarkable improvements from 3D processing techniques.Furthermore, the latter are very time-consuming and rarely available on common GPR processing software.So far, although widely appreciated, 3D acquisition techniques, using submeter spacing, have been applied only in test-areas of limited extensions or, using a coarser spacing, also in larger areas, but to image only large-scale archaeological features (Goodman et al., 1994;Malagodi et al., 1996;Pipan et al., 1999;Basile et al., 2000).
The use of 3D visualization techniques is of primary importance in archaeological applications in order to display complex data in an easily understandable fashion, thus improving the quality and efficiency of the archaeological interpretation.The most widespread way to display 3D radar data is in «time slice» (or depth slice) maps (Conyers and Goodman, 1997).Horizontal slices may not be the more suitable visualization technique in the case of great subsurface complexity since, for example, false amplitude anomalies can occur when the slicing planes cross dipping or undulating reflectors.However, time slices still remain the easiest and most rapid means to provide a plan synthetical view of the anomaly pattern, especially for large areas.For small zones a more complete understanding of the subsurface can be achieved by means of various 3D data presentations, including 3D cubes, chair views and slices parallel to the axes or along arbitrary directions.Two interesting visualization approaches, involving the extraction and 3D visualization of the most promising signal attributes, have been proposed in recent papers (Zanzi and Valle, 1999;Valle et al., 2000), where they have been successfully applied in imaging three dimensional bodies in mine detection and Non-Destructive Testing applications.
In the present work, a trial of application of the two above-mentioned 3D visualization techniques, along with classical time slice representation, in archaeological prospecting has been made.The results obtained from two different situations encountered in an archaeological case history («Buon Pastore» -Lecce, Italy) are presented, which allowed a comparison of the performance of the different visualization techniques, outlining the main advantages and drawbacks of each one.

Site description
The study area (fig. 1) is located in Lecce (Italy), about 500 m NW of «Porta Napoli», near the present cemetery and the Norman monastery of «SS.Niccolò e Cataldo» (12th century).The site is just outside the north-western border of the ancient city.It has remained almost unchanged from the Messapian to the modern age, as testified by the local coincidence of the quite well-preserved Aragonese walls with the Roman and the Messapian ones.The potential archaeological interest of the site arise from the documented recovery of both Messapian (4th-2nd century B.C.) and Roman tombs in the neighbourhood.The typology and sizes of the tombs recovered in Lecce are quite variable.The burials can be classified in three main types: -«Ipogei» (hypogeum type: dimensions of the order of 3 m × 1.5 m × 1.5 m).
-«Tombe a fossa» (grave pit type: average dimensions of 1.8 m × 0.5 m × 0.5 m; some of which covered by stone slabs).
The second type seems to be the most common.Also the depth of burial is highly variable from a minimum of 1 m to a maximum of about 4 m from the present ground level (Giardino, 1994).As well as for the Classical age, the site is of archaeological interest also for the Middle Ages, since in past times it could have been belonged to the «SS.Niccolò e Cataldo» monastery and, potentially, remnants of structures related to the monastic property could be buried there.
At present the site belongs to the University of Lecce, which intends to enlarge the existing building, named «Buon Pastore», by adding further edifices.The risk of destruction of potential archaeological structures during the building works motivated archaeological investigations.The relatively large dimension of the area (almost 10 000 m 2 ), together with time and cost constraints, made necessary the recourse to geophysical investigations as a faster means to ascertain the presence of relevant archaeological features and to delineate the geological

N
stratigraphy, with particular regard to the calcarenite bedrock.Since the most important targets (tombs or wall remnants) were expected to be located at the soil-bedrock boundary, its depth was an important parameter for planning archaeological tests.

Data acquisition and preliminary analysis
The interpretation of a previous four-lines seismic refraction survey (S1, ..., S4 in fig.2) led to a simple two-layer model for the entire area: an upper soil or weathered layer (V p 500-600 m/s), extending from ground level to a relatively constant depth of about 1 m (locally increasing to less than 2 m), and a deeper layer, interpreted as the calcarenite bedrock (V p 2000 m/s).Taking into account the two-fold objective (stratigraphical and archaeological) and the seismic indications, a GPR survey, using a GSSI SIR-System2 equipped with 200 MHz and 500 MHz antennae, was performed in autumn 1999 and spring 2000.Because of the variety of sizes of the expected archaeological targets, a 1 m line spacing was used for the reconnaissance survey, carried out in continuous mode on all the accessible parts of the site.Apart from rare exceptions, due to the presence of trees, shrubs and other obstacles, almost all profiles were surveyed with both antennae.Subsequently, on the basis of the preliminary survey indications or proximity to zones being excavated in the meantime (test 1 to 7), some areas were selected (zones A, B and C), where the crossline spacing was lowered to 0.50 m and 0.25 m for the 200 MHz and 500 MHz surveys, respectively.Moreover some CMP gathers (dots in fig.2) were carried out for velocity estimations.Despite some small differences, they allowed us to estimate an average velocity of 0.09 m/ns for the entire area, also confirmed by the velocity analysis made on diffraction hyperbolas in the radar sections.This average value was used for depth conversion and migration.
The radar sections shown in figs.3a-d and 4a-c, exemplify the signal behaviour in different parts of the site, probably related to different rock and soil conditions and to variable moisture content.Figure 3a-d shows the sections relative to two profiles acquired in the southern part of the site (L2 and L1 in fig.2).On the 200 MHz section relative to profile L2 (fig.3a) at time ranging from about 25 ns to about 35 ns (1.20 m ÷ 1.60 m in depth), a reflection event, slightly undulating and sometimes interrupted by diffraction hyperbolas, is easily identifiable.Because of its high amplitude, denoting a strong electromagnetic contrast, this event was interpreted as due to the soil-bedrock interface.The same reflection, on the other hand, is barely recognisable on the 500 MHz section (fig.3b), due to the proximity of the antenna frequency to the «clutter frequency» (Annan and Cosway, 1994), which causes an enhancement of the soil heterogeneity diffraction effect.Apart from the bedrock reflection, a very strong anomaly (as that denoted by C in fig.3c,d) was noted on the radar sections relative to L1 and to some profiles close to it.
In fig.4a,b which refers to a profile acquired in the northern part of the site (L4 in fig.2), the bedrock reflection appears shallower (almost flat at about 20 ns) and more continuous.Moreover, the overlying soil appears more homogeneous, while at times greater than 20 ns, although disturbed by noise, a slightly dipping event can be identified, probably related to fractures or bedding planes in the calcarenite bedrock.In such more favourable conditions the bedrock reflection is quite well distinguishable also on the 500 MHz section, while the more visible diffraction hyperbolas better mark its interruptions.As can be seen from the section (fig.4c) relative to a profile from the middle part of the site (L3 in fig.2), probably because of more gradual changes, the bedrock reflection is less evident even on the 200 MHz antenna (and sometimes is overlain by a shallower reflection).Moreover, the signal behaviour in the deeper part of the section could be due to a lower degree of bedrock integrity (fractures) in this zone.
From the interpretation of all the radar sections, the estimated bedrock depth varied between 1.2 and 2.0 m in the southern part and between 0.8 and 1.2 m in the northern part of the area, in agreement with the results of the seismic survey and of the archaeological excavations.Moreover, the bedrock reflection was generally characterized by marked roughness (although less or more evident in different parts of the site) and in the overlying soil layer many anomalies of various types were noted.At a first sight, both the reflector roughness and the soil anomalies could be due to natural or anthropic causes, so recourse to visualization schemes which reveal geometrical (and in particular horizontal) relationships between the anomalies appeared unavoidable to resolve this ambiguity.As a first trial, two zones of different radar signal complexity were selected to test the methodology, namely zone D (fig.2), where a very strong anomaly of limited size was noted (fig.3c,d), and zone B (fig. 2) characterized by weaker but more complex GPR signal (fig.4c).

Data processing and visualization
Radar data were previously processed using standard 1D and 2D techniques, mainly consisting of: trace editing and horizontal normalization (0.1 m inline spacing), spectral analysis and band-pass filtering, horizontal highpass filtering (background removal), gain control and 2D Kirchhoff migration using the average velocity value of 0.09 m/ns.3D Kirchhoff migration was tested on the 500 MHz data from zone A, but the minor improvements obtained in comparison to 2D migration, confirmed that even in this case the spacing (0.25 m) was too large for 3D migration to be effective.2D migration appeared the best-suited method for the largerspaced 200 MHz data.The processed data were then visualized with three different techniques whose main characteristics are briefly outlined.t.Sometimes a particular complex-trace attribute, the instantaneous amplitude or envelope (modulus of the Hilbert transform), is used instead.Being a measure for reflectivity strength, it helps to evidence high amplitude anomalies.Previous spatial averaging is also useful to reduce small-scale heterogeneity noise.Finally data are interpolated and gridded on a regular mesh (Conyers and Goodman, 1997).Selecting the various parameters involved (Basile et al., 2000) and in particular the width of the slice, t, is crucial.Typically t must be of the order of the dominant period, but different widths can be used to enhance particular features.In common practice, non-overlapping time windows are chosen, although sliding windows could be used instead, with the advantage of greater resolution but higher computational costs.
A two-fold approach for visualizing 3D radar data has been proposed by Zanzi and Valle (1999) for automatic mine detection.In the first case, after an appropriate processing of radar data, a 3D image of the sought diffracting or reflecting object could be easily obtained by: 1. Extraction of a particular complex signal attribute (trace envelope).
Whereas this was effective in the case of a laboratory experiment, the low signal-noise ratio observed in a real case induced the authors to propose an alternative approach consisting of: 1. Extraction of the most promising complex signal attributes (trace energy and envelope).
2. Three stacks separately performed along each coordinate axis, providing separate 2D results: stacking of the energy along the depth or Z axis, in order to obtain a plan view of the high-energy suspected zones; stacking of the trace envelope along X; stacking of the envelope along Y.
3. Thresholding.4. 3D rendering of the presumed target by projection in 3D space of the automatically selected thresholded data.
As pointed out by the authors, in both cases the threshold calibration is a very delicate task.
In the following figures two examples of application of the three visualization procedures in an archaeological context are given.In both cases the 200 MHz GPR data were used because they are less noisy than the 500 MHz ones.For every figure the viewpoint is from the south corner.

Zone D (crossline spacing: 1 m)
In fig.5a,b two time slice representations, using the absolute amplitude, are shown for zone D, obtained using the same grid cell sizes, but different time intervals: 8 ns (or approximate-ly 0.36 m depth windows) and 16 ns (about 0.72 m).A clear rectangular anomaly, 6 m long and 4 m wide, is visible from the 8 ns (0.36 m depth) to the 32 ns (1.44 m) slices.Moreover, although weak and thin, a NE-SW linear anomaly is visible on 24 ns and 32 ns slices ending at the base of the main anomaly.It corresponds to the hyperbola seen on radar sections (as, for example in fig.3a) slightly deeper than the soil-bedrock reflection.The higher mean amplitude on slice 24 ns (1.08 m) is related to the soil-bedrock reflection, whose roughness is evident on this slice.The same main lineaments are also evident in fig.5b, but the coarser slicing in t axis could lead to a less correct depth estimation, since for the non-overlapping window choice the depth error is at least half the depth window.As expected, amongst the other parameters, the time window selection appeared the most critical.Apart from minor differences, using the square amplitude or the modulus of the Hilbert transform led to analogous conclusions.
In fig.6a-e the same data set is displayed with iso-amplitude surfaces using two threshold values: 50% (a) and 40% (b) of the maximum complex trace amplitude.Obviously, lowering the threshold value, increases the visibility of the main anomaly and smaller objects (the linear anomaly denoted by an arrow in fig.6b), but also heterogeneity noise.To better evidence the main characteristics, three different views, obtained by simple rotations of the 3D volume (for the 40% threshold) are displayed: plan (c), front (d) and lateral view (e).In particular, the latter figures clearly show that the small-size anomalies are located at depth corresponding to the soilbedrock boundary (about 1.2 m), while the main body has a nearly flat top at about 0.5 m and a thickness of about 1 m.
Following the third approach, the trace energy and envelope were extracted from the processed data, stacked along the depth (energy) and along X and Y (envelope) axes, normalized to the respective maximum and visualized as the three coordinate planes of a 3D volume (fig.7a); after a careful selection of the threshold value (fig.7b), data were projected in 3D space to produce the final display (fig.7c).Although a relatively strong continuous reflection is visible on the stacked and thresholded volumes (between 20 and 30 ns), the backprojection procedure images only the most energetic 3D anomalies.Using a different attribute for the top plane, e.g., the complex trace amplitude, no major differences were noted in this case, provided a careful selection of the threshold value is made, which appeared in both cases the most delicate parameter.

Zone B (crossline spacing: 0.5 m)
The same visualization techniques were employed on the zone B processed data set.In the time slice technique the lower signal/ noise ratio than in zone D is clearly apparent from the more chaotic anomaly pattern on nearly every slice (fig.8a).A general character separates the two areas.Moreover, some localized strong anomalies are visible on different slices at the same location, probably indicating near vertical objects, although we could not rule out the possibility that the vertical distribution of energy could also be due to reverberations occurring at the top of the targets, because of local antenna/soil characteristics.The west-east trend of the amplitude anomalies with depth and the central low-amplitude zone may be more apparent using the 16 ns time window (fig.8b), but this implies less precision in the vertical direction.
An even lower threshold value (42% or 37%, as in fig.9a,b), than that used for zone D, is needed in this case for the iso-amplitude visualization, because of the absence of highly reflecting bodies with respect to the host environment.The already outlined west-east differences are easily visible through this technique both in the perspective (b) and in the plan views (c), whereas the lateral views (d and e) better evidence the shallower location (less than 1 m) of the soil-bedrock anomalies and, to a much lesser extent, some dipping features (lines).
Within the third technique, by stacking the envelope along X and Y coordinates, the visibility of the weak dipping features was slightly increased with respect to the host diffraction energy (fig.10a).Nevertheless, the thresholding and 3D rendering procedure inherently enhanced the contribution of threedimensional bodies in comparison to that of 2D (surfaces) or 1D (linear) features, so that the final image (fig.10c) is quite similar to that obtained by the previous technique (fig.9b).In this case, besides the threshold value, the technique is also rather sensitive to the complex attribute (energy or envelope) used for the top plane (in this case the envelope was used also for the top plane instead of the energy).
With respect to the results obtained in zone D, a general consideration is that, in this lessfavourable case, every technique suffered from the lower signal/noise ratio.

Data interpretation
In zone B the anomaly concentration on the western corner at the soil-bedrock boundary, indicating a local higher roughness of the interface, could be compared with the results of the nearest archaeological excavation (test 6 in The strong squared anomaly found in zone D seemed unlikely to be an hypogean tomb, both in base of GPR results and oral tradition.Its shape and size, clearly evidenced by all the visualization techniques adopted, and the connected linear anomalies, visible on particular slices using proper time windows, suggested that it could have been an old refilled cistern with connected water canals probably excavated in the rock. On the other hand, no evident grave pit anomaly was found by the analysis of the detailed GPR survey data, neither indication of grave pit presence was derived by the reconnaissance survey, although in the last case the large crossline spacing (1 m) could not completely exclude this possibility.However, no grave pit was found in any of the archaeological tests.

Discussion and conclusions
The primary importance of 3D visualization comes from the fact that it can provide a powerful and intuitive means of communicating complex information to non-geophysicists.Aside from the archaeological meaning of the features imaged, the present study allowed a comparison of the performance of different visualisation techniques of GPR data in relatively complex environments.
-Time slice maps, used mainly to enhance the horizontal relationships between amplitude anomalies at the almost flat soil-bedrock interface, evidenced some alignments along directions consistent with those of the artificial cuts observed in the nearby excavations and allowed distinguishing zones probably more heavily quarried from other less disturbed.Obviously this method failed in imaging dipping events.They are a more objective representation of the results than the other methods explored, but do not furnish an instantaneous view of the entire volume with the same immediacy.
-The 3D contouring and isoanomaly plotting suffered from a certain degree of subjectivity in selecting the threshold value, but furnished very impressive pictures of 3D (and to a lesser degree also 2D or 1D) reflecting bodies, provided the data had a relatively high S/N ratio.
-The alternative visualisation approach proposed, which involves a procedure of stacking, thresholding and 3D rendering of selected data attributes, has a higher degree of subjectivity (selecting the trace attribute and the threshold value), but is the most selective and effective in imaging only true 3D bodies.
In the archaeological case study proposed, none of the visualisation methods was able to resolve the singular bedrock cuts due to intrinsic limits: -Low thickness and width in comparison to the vertical and spatial resolution allowed by the antennae and acquisition geometry used.
-Relatively high depth that prevented the use of higher frequency antennae.
-Proximity of the archaeological target sizes to the natural bedrock undulation wavelengths and to the dimensions of the host heterogeneity.
In reality, 3D visualisation techniques, which enhance particular characteristics present in the data, can give very impressive images for understanding the 3D spatial relationship and, hence are very helpful for archaeological or geological interpretation.Nevertheless the fundamental rule remains the same: features must be detected before  they can be enhanced, so that particular care is needed both during the acquisition and the processing stages.

Fig. 2 .
Fig. 2. Location map of the geophysical surveys: seismic refraction (dashed lines S1,..., S4) and GPR (continuous lines); heavy lines refer to profiles whose sections are shown in the next figures.A, B and C denote zones of detailed 3D GPR surveys.The dots mark the centres of the CMP gathers carried out for velocity estimation.Tests 1 to 7 denote archaeological excavations.

Fig
Fig. 3a-d.Examples of radar sections from profiles relative to the southern part of the site: Profile L2, 200 MHz (a) and 500 MHz (b); Profile L1, 200 MHz (c) and 500 MHz (d).I denotes the soil-bedrock reflection; H a hyperbolic diffraction probably due to a ditch in the bedrock; C a strong anomaly probably due to an old refilled cistern.

Fig
Fig. 4a-c.Examples of radar sections from profiles relative to the northern part of the site: Profile L4, 200 MHz (a) and 500 MHz (b); Profile L3, 200 MHz (c).I denotes the soil-bedrock reflection.

Fig
Fig. 5a,b.Zone D. Time slices obtained averaging the absolute amplitude of the data within different time windows: 8 ns (a) and 16 ns (b).The strong rectangular anomaly C is probably related to an old refilled cistern; arrows denote a weak and thin linear anomaly probably due to a ditch in the bedrock.a b

Fig
Fig. 8a,b.Zone B. Time slices obtained averaging the absolute amplitude of the data within different time windows: 8 ns (a) and 16 ns (b).

Fig. 11 .
Fig. 11.Photo of the excavation results from Test 6: artificial cuts in the bedrock interpreted as ancient (maybe Messapian) quarry works.

Fig
Fig. 12a,b.Photo of the founding pit wall close to zone B revealing the presence of colluvial deposits and fractures in the calcarenite bedrock dipping to the south (a) and east (b).