Abstract
The Romano-British site of Barcombe in East Sussex, England, has suffered heavy postdepositional attrition through reuse of the building materials for the effects of ploughing. A detailed GPR survey of the site was carried out in 2001, with results, achieved by usual radar data processing, published in 2002. The current paper reexamines the GPR data using microwave tomography approach, based on a linear inverse scattering model, and a 3D visualization that permits to improve the definition of the villa plan and reexamine the possibility of detecting earlier prehistoric remains.
1. Brief Description of the Archaeological Site
This paper deals with the reexamination of the data collected during GPR surveys collected with the aim of a prospecting of an ancient Roman Villa in Barcombe. The villa lies in a field on the perimeter of the Barcombe village in East Sussex, England (see Figure 1). The site came to the attention of the Mid Sussex Field Archaeological Team (MSFAT) and the University College London Field Archaeological Unit (UCL, subsequently replaced by the Centre for Continuing Education of the University of Sussex, CCE) because it was in danger of disappearing altogether without being adequately recorded [1]. In common with many other UK sites of the period, the site had been extensively robbed out in the centuries following its demise in order to provide building material for the adjacent village and its associated farms. This is a common problem with Romano-British sites in the UK. In addition, the site is positioned on the ridge of a field in agricultural use and has therefore been extensively ploughed out. As a result, the archaeological evidence was sparse and the remained part was being rapidly eroded.
2. GPR Survey and First Results
In April 2001, a Ground Penetrating Radar (GPR) survey was jointly carried out by Professor AM Alani of the Department of Civil Engineering of the University of Portsmouth and by Utsi Electronics Ltd, on behalf of the archaeological team put in charge of exploring the villa. The work was aimed at investigating the possibility of mapping both the villa and earlier prehistoric remains on the same ridge.
The survey also provided comparative material for other geophysical survey methods employed by the archaeologists, especially the electrical resistivity method. The results of the GPR survey, including a comparison with the evidence from the electrical resistivity field test, were published in 2002 [2].
Using a 40 m by 60 m grid laid out by the archaeological team, a Groundvue 1 system, manufactured by UTSI Electronics, was used to perform a survey along parallel transects at intervals of 50 cm (Figure 2). The GPR system was equipped with antennas of central frequency 400 MHz, linearly polarised, parallel to each other and orthogonal to the direction of the B-scan. The sampling interval along the line of survey was 5 cm, and probing was carried out to 40 ns. A Wide Angle Reflection and Refraction (WARR) calibration [3], carried out at the time of the survey, provided a value of the relative dielectric permittivity equal to 12.5, which implies a propagation velocity of the electromagnetic waves in the soil equal to about 0.085 m/ns. At this velocity, 40 ns is equivalent to an investigation depth of 1.7 m.
Two partial orthogonal surveys were carried out at a transect spacing of 1 m in order to check that the survey direction was not materially biasing the results. For the case at hand, the main limiting factor on the achieved GPR images proved to be the scarcity of remains, rather than the direction of survey, and thus no further use was made of the orthogonal survey.
The original GPR data were first processed (using the ReflexW software package [4]) by applying background removal [5], adding time-based gain, averaging over 2 traces in order to reduce noise resulting from the relative movement of the antennas across the ploughed field and finally applying a Bandpass Butterworth filter from 200 to 600 MHz (the nominal frequency range of the antennas). The 2-dimensional profiles were regrouped so to build a 3-dimensional cube and the constant time slices extracted [6]; on the basis of maximum signal return, the constant time slices at 16 ns, 25 ns, and 29 ns, about corresponding to the depth 68 cm, 106, cm and 123 cm, are here presented (Figures 3, 4, and 5). In particular, the images of Figures 3–5 correspond to the rectangle in dashed line in Figure 7.
Although the overall outlines of the villa reflected a pattern similar to those from the electrical resistivity (ER) surveys (see Figure 6, where these features are all visible together), the GPR time slices also illustrated the variation of patterning with depth. This variation along the depth appears to be primarily the result of differential stripping out of the stone for reuse elsewhere. A secondary effect was the necessity of underpinning the villa along its northern and western edges, presumably because of slippage due to the adjacent backfilled prehistoric ditch (Figure 7).
Processed data (time slices) did not produce any evidence of extended floors despite the presence of tesserae in and around the surveyed site. Subsequent excavation of the site did not reveal any evidence of any preserved floor either.
The prehistoric remains identified in the original site plan were not confirmed in the time slices produced, presumably because the electromagnetic properties of these features were too close to those of their surrounding environment.
Let us outline that, actually, some additional survey transects were completed to the south of the villa too, but this data was not used for interpretation because the results had been heavily affected by signals from small metal objects such as washers and coins buried within the first few centimetres of land surface. These metal objects were originally sown over the entire site in order to prevent pilfering of the site by local metal detector operators. By the time that the GPR survey was carried out, most of the metal had drifted southwards under the influence of rain and ploughing to the lower part of the field. Since no useful data could be extracted from this area, the data from the final 10 m to the south of the survey was excluded from interpretation and analysis.
3. Reexamination of the Data: The Microwave Tomographic Approach
Although the main limiting factor of this data set is and always will be the relative lack of material remains, re-processing GPR datasets with alternative methods has provided interesting results in other archaeological and civil engineering investigations [7, 8].
The decision was therefore taken to reexamine this dataset to see if the delineation of the villa could be improved and whether or not it was possible to detect the traces of prehistoric activity in the images.
This motivated the choice to reprocess the data by a microwave tomographic approach, which is an inverse scattering technique applied to the scenario at hand [9–14]. In particular, in this paper, we adopted a strategy aimed to provide a pseudo-3D representation of the investigated scene. This approach has been already applied with regard to the archaeological site of Ballachulish [8, 15].
The approach can be summarised as follows: first, a 2D linear scalar tomographic algorithm is applied to each B-scan (transects); afterward, the 2D tomographic reconstructions are joined to build up a 3D cube, from which horizontal or vertical slices can be extracted. In this paper, we will present the reconstruction outcomes as horizontal depth slices. The 2D tomographic reconstructions have been achieved thanks to a microwave imaging algorithm based on the Born Approximation (BA) that has been already described in detail in other papers [12, 13, 16] and falls within in a well-known specific literature [9–18]. BA is useful when penetrable objects are under investigation, as in the case of this paper where the attention is mainly focused on buried archaeological artefacts. As is well known, BA provides accurate results for the forward problem, that is, the determination of the field scattered from a known object, under the hypothesis of weak scattering targets [19, 20], that is when the dielectric characteristics of the buried objects are only slightly different from that of the host medium and when the extent of the target is not too large in terms of probing wavelength.
When an inverse scattering problem is tackled, the hypothesis of weak scattering from the target can be relaxed at the cost of renouncing to the “quantitative” reconstruction of the dielectric permittivity. In fact, the adoption of a Born model inversion scheme allows to detect, to localize, and to determine the geometry of the buried target even in the case of strong scattering objects, as has been theoretically postulated on the basis of the support of the secondary sources [19, 21] and as many numerical and experimental tests have shown [8, 22–24].
Here, we briefly present the adopted inversion approach just for consistency purposes.
The two-dimensional geometry of the problem is depicted in Figure 8. The background inhomogeneous scenario is built up of two homogeneous half-spaces separated by a planar interface at z = 0. The upper half-space is made up of air with absolute dielectric permittivity . The lower half-space (supposed homogeneous) shows a relative dielectric permittivity , possibly complex to account for losses in the propagation medium. The incident field source is a time-harmonic (time dependence exp(j2πft)) filamentary y-directed electric current, invariant along the y-axis. The data are multifrequency in the band []. A multimonostatic measurement configuration is considered where the transmitting and receiving antennas are placed in the same point. The relevant targets are invariant along the y-axis, and their cross-section is enclosed in the rectangular investigation domain (see Figure 8). We assume the relative dielectric permittivity profile (possibly complex to account for losses in the targets) inside the investigation domain D as the unknown of the problem. Accordingly, the inverse problem is recast in terms of the contrast function [13, 16, 22], defined as Under BA, the relationship between the unknown contrast function and the scattered field data is provided by the following linear integral equation in frequency domain [9, 13, 16, 22]: where is the circular frequency, is the wave-number in the soil, and is the electric scattered field along the air-soil interface collected at the abscissa (let us remember that we are considering a multimonostatic configuration, that is, the source and observation points are the same). The scattered field is given as the difference between the total field and the incident (or unperturbed) field. The incident field is essentially given by the sum of the direct talk between the antennas and the field reflected by the air-soil interface. When the antennas are moved at the air-soil interface, these two subparts are customarily superposed in time domain and also their spectra develop on the same band, so that it is hard to distinguish them from each other. The scattered field is instead the signal backscattered by the buried objects. is the external Green’s function of the problem. The expressions of the Green’s function and of the incident field, accounting for a half space 2D scalar model with the observation line put at the air soil interface, are given by [9, 13, 19]: where being the imaginary part of the square root less than or at most equal to zero. In (4), the source is assimilated to a filamentary current, where is the level of the current, and the absolute magnetic permeability of the free space, equal to H/m.
For a more detailed dealing, reader is referred to [9, 13, 17] and references therein. The linear integral relationship (2) is discretized into a linear algebraic system by means of the Method of Moments, where the unknowns are the expansion coefficients of the contrast function along the chosen functional basis. Point matching is adopted in the data space. The inversion of the resulting matrix is performed by a Truncated Singular Value Decomposition (TSVD) [25], which provides a satisfying robustness of the solution algorithm with respect to the uncertainties and the noise on the data.
The reconstruction results will be given as the spatial map of the modulus of the retrieved contrast, and therefore the regions where the modulus of the contrast function is significantly different from zero express the position and the geometry of the buried objects.
Some preprocessing is needed to pass from the total field GPR data collected in time domain to the scattered field data in frequency domain suitable for the inversion in frequency domain. This is accomplished by means of a two-step procedure. The first one is to mute (after zero timing) the first part of the time domain traces, which approximately corresponds to erase the direct coupling between the transmitting and receiving antennas and the contribution due to the air/soil interface. This step provides an estimation of the scattered field in time domain directly from the data. The second step consists in Fourier transforming the time-domain scattered field so to achieve the scattered field data in frequency domain.
4. Inversion Results
The data processing regarded 63 profiles with a length of 56 meters, so the reconstruction was concerned with an area of m2. The processed B-scans have been then exploited to achieve a pseudo-3D visualization of the investigated area [8, 15]. The processing was performed with a common PC, and it was not possible to process simultaneously all the data contained within a B-scan profile (transect). Thus, we adopted the following strategy: for each 56 m long scan (transect) profile, we performed 14 tomographic reconstructions (each of which relative to a 4 m long investigation domain). The parameters of the inversion for each single reconstruction are reported in the Table 1.
Accordingly, by considering the central frequency 325 MHz, achieved from the data (which means m), the volume of the investigated scene is .
The results are presented in terms of constant depth slices achieved through the following steps(i)The 3D volume is built up by joining side to side the 63 tomographic reconstruction of 56 m length;(ii)The constant depth slice is fixed and visualised in two modalities below reported. The first “visualization modality” is carried out by normalizing the modulus of the contrast function at the depth at hand with respect to its maximum value. A threshold at 0.2 has been applied too in order to enhance the visibility of the weakest scattering targets. Hence, the represented quantity is the normalized modulus of the retrieved contrast where the final scale of the images is (0.05–0.2).
The second “visualization modality” is achieved performing a preliminary interpolation of the modulus of the contrast function over the constant depth section by a median filtering. This filtering is two dimensional, and each output pixel contains the median value in the M-by-N (M = 3, N = 6) neighborhood around the corresponding pixel in the input image (this procedure is carried out by the routine MEDFILT2 command in MATLAB).
The tradeoff between the two visualizations is that the first one gives a cleaner image and so makes possible a better interpretation of the overall content of the buried scenario, whereas the second one enhances some weaker scatterer otherwise possibly “covered” by the stronger scattering targets.
Figures 9–13 depict the constant depth slices at depth of 0.5 m, 0.63 m, 0.75 m, 1, 1.25 m, respectively; for each figure, the panels (a) and (b) depict the reconstruction according to the first and second visualization modality, respectively.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
5. Results Interpretation, Conclusions, and Future Developments
When comparing the results of the resistivity survey, the excavation, and the GPR survey, it is important to realise that the area covered by the surveys is in any case only a part of the full excavated area. All three depictions of the villa (resistivity, classical GPR processing, microwave tomography) show a winged corridor villa. The ER data also indicate the presence of a circular anomaly at the southern side of the villa. This feature does not appear in the GPR data because it lies in an area containing many small metal objects which obscure the materials underneath. On the other hand, only the GPR survey provides high-resolution details about the depth of the anomalies. The distribution of strong signal returns in the time slices indicates the position of some flints.
At c. 69 cm, the main outlines are those of the villa wings, as outlined by the red dashed circles in Figure 3. By c. 1.08 m, although the western wing is still visible, the major feature is constituted by the corridor to the north and specifically the western end of this corridor. This is put into evidence by the dashed red ellipse in Figure 4. A similar pattern but with more details about the eastern end is visible by c. 1.23 m. This is put into evidence by the dashed red ellipse in Figure 5. The eastern boundary wall is also visible at this depth, as put into evidence by the dashed red rectangle in Figure 5.
Since the strongest signal responses from the wing foundations occur at a shallower depth than those from the footings of the northern corridor, the time slice sequence might suggest that the wings are successive (Figures 3 to 5). However, the excavation evidence shows that this is not the case [2], and a more likely explanation lies in the need to reinforce the footings in the area of the underlying prehistoric ditch along the inside of the west and north walls in order to either prevent or remedy the potential undercutting of these two walls (Figure 7). The excavation results indicated that the underlying ditch had caused structural problems which had required additional underpinning. Although the ditch is not itself visible in the time slices, the extent of the strong signal returns in the time slice is a good match for the position of this feature. In particular, it is put into evidence in Figures 4 and 5 by means of two yellow dashed lines and is outlined in its entirety in Figure 7.
Let us now turn to the time slices from the reprocessed GPR data, that as said have been taken at 0.5 m, 0.63 m, 0.75 m, 1 m, and 1.25 m.
The first of the slices at 0.5 m depth (see Figure 9) is shallower than the “traditionally processed” time slices. In addition to the two wings being visible, it indicates some building remains associated with the southern edge of the corridor that are not visible in the figure at the depth of 16 ns (Figure 3). The pattern is similar to the ER data of Figure 6.
At 0.63 m, the outline of the main building is very similar to that at 16 ns, corresponding to 69 cm (compare Figures 3 and 10). The southern wall of the villa appears more straight within the inversion result, and there are two faint parallel lines toward the northern part of the villa, put into evidence by the dashed rectangle in Figure 10. Unfortunately, this area was not excavated and it is therefore impossible to confirm that these are archaeological features. The strength of the signal in the slice, however, is similar to that of the eastern boundary wall as processed in Figures 4 and 5. This makes us think that these two faint lines might correspond to real anomalies.
At 0.75 m (Figure 11), the result is again similar to that achieved from the conventional processing at 0.69 m (Figure 3). The two faint parallel lines to the north of the villa are still visible (and again put into evidence by the superposed ellipses), and the definition of the back (north east) wall of the villa corridor is clearer than in the 0.69 m time slice (compare Figures 3 and 11). The eastern faint parallel line is close to the position of an underlying ditch visible in the excavation plan (see Figures 7 and 11), but it is not possible to state if the two features are in fact related. The reprocessed data provide the hypothesis that these two linear features may form three sides of a rectangular room immediately adjacent to the villa toward its northern part (Figures 7 and 11).
Unlike the shallower ones, there is a very striking difference in the time slices at about 1 m depth (see Figures 4 and 12). The conventionally processed data (Figure 4) highlights the north and south walls of the rooms along the northern edge of the villa and some details to the East of the building. With a coarser definition, but still clearly visible, part of the corridor to the south of these rooms can be appreciated too, as well as both the wings of the villa and the footprint of the boundary wall leading southwards from the villa. In the reprocessed data (Figure 12), the time-slice at the depth of 1 m in the first visualization modality shows only western wing and the junction of two walls. In the second visualization modality, instead, although the lines are faint, a large part of the structure of the villa walls is visible. The results at depths of about 1.25 m are more similar (see Figures 5 and 13). At this depth, most of the outline of the villa, with the notable exception of the east wing, is visible, and the definition in the conventional processing is sharper at this depth. In conclusion, some details seem to be better identified from the first processing and some from the other one, so that the interpretation can be enhanced by the comparative analysis of the results achieved under both approaches.
6. Conclusions
In this paper, we have presented a reprocessing of archaeological GPR data already studied before with a classical migration-based processing. Of course, the reprocessing of the Barcombe data through microwave tomography techniques cannot substitute for the lack of physical remains on site. Notwithstanding, it does show that it is possible to detect some fainter feature than those previously detected, notably in the 0.5 m time slice (where the data is effectively validated by both the excavation plan and the resistivity survey) and in the 0.75 m time slice (that suggests at least one additional feature). Had the processing been carried out at an earlier date, it would have been possible for the archaeological team to investigate this area. Unfortunately, this is no longer possible. The reprocessing of the data has illustrated the complementary nature of the two processing techniques and has a potential usefulness in the detection of faint archaeological remains.
Acknowledgments
The authors would like to thank the late David Coombs of the Mid Sussex Field Archaeology Team for the data from the resistivity survey; David Rudling of MSFAT and CCE for the excavation plan and information on the excavation. This work has been performed as an activity of “Active and Passive Microwaves for Security and Subsurface imaging (AMISS)” EU 7th Framework Marie Curie Actions IRSES Project (PIRSES-GA-2010-269157).