GPR Investigation at the Archaeological Site of Le Cesine, Lecce, Italy

Abstract

In this contribution, we present some results achieved in the archaeological site of Le Cesine, close to Lecce, in southern Italy. The investigations have been performed in a site close to the Adriatic Sea, only slightly explored up to now, and where the presence of an ancient Roman harbour is alleged on the basis of remains visible above all under the current sea level. This measurement campaign has been performed in the framework of a short-term scientific mission (STSM) performed in the framework of the European Cost Action 17131 (acronym SAGA), and has been aimed to identify possible points where future localized excavation might and hopefully will be performed in the next few years. Both a traditional elaboration and an innovative data processing based on a linear inverse scattering model have been performed on the data.

1. Introduction

Ground Penetrating Radar (GPR) is a useful tool to study and document cultural heritage sites [1,2,3]. Archaeology is one of the classic fields of application of GPR prospecting [4,5], especially when localized or even extensive excavations are scheduled. However, the areas of potential interest are customarily too large for an excavation of the entire zone, which would result to be slow and expensive and might not provide the expected results, and so a preventive geophysical investigation can reduce the extensions of the areas where targets of possible archaeological interest are expected, even if within the natural (technological but, above all, physical) limits of the technique (in particular, the not-infinite available spatial precision and resolution). In most cases, data processing is performed by means of a commercial code, and there are several codes for this, among which the most popular nowadays are the Reflexw (https://www.sandmeier-geo.de/reflexw.html (accessed on 4 October 2021)) and the GPR-SLICE (https://www.gpr-survey.com/ (accessed on 4 October 2021)). Alternatively, or as a complementary tool, linear Microwave Tomographic (MWT) approaches, which face the imaging as an inverse scattering problem (by exploiting the Born Approximation to model the wave-materials interaction underling the GPR survey) can be exploited too [6,7]. To these pros, a valuable example dealing with a large set of field data gathered at the Archaeological Park of Paestum and Velia is provided in [8], where a computational effective implementation of the MWT approach has been used.

The first attempts to apply an inverse scattering approach [9,10] to an electrically large domain have been based on the joining of adjacent reduced investigation domains, related to overlying corresponding adjacent subsections of the comprehensive observation line [11]. The fractioning of the investigation/observation domains, however, poses a relevant problem of reduced view angle [12] for all the targets close to any bound between two adjacent investigation domains. In [13], a method for mitigating such a problem has been introduced, called shifting zoom (or “shift and zoom”), and in [14] it has been applied to the creation of depth slices based on a linear inverse scattering algorithm. In this contribution, thanks to the shifting zoom, we will apply an inverse scattering algorithm to three electrically large sets of data and we will create depth slices based on this procedure. We will do this in relation to the archaeological site of Le Cesine, close to the Adriatic Sea and on the outskirts of Lecce, in southern Italy. We have checked the penetration time-depth of the signal on the basis of the results, and it was quite satisfying (about 50 ns) in spite of the closeness of the sea, probably because the soil is quite rocky and hinders the infiltration of the salty sea water beyond the shore. In particular, three areas have been investigated close to archaeological remains, visible above all under the sea level close to the shore. To these pros, it is also important to say that it is estimated that, at the Roman Age, the sea level in this area was about two meters lower than the current level. Indeed, these remains have been well known to archaeologists for years, but the hypothesis of a Roman harbour had not been formulated until recent studies, also due to the fact that, at 1.5 km from the area considered in this paper, there is another well-known dock of Roman Age [15]. Consequently, these recent archaeological investigations have opened up an intriguing problem regarding the reason why two Roman harbours would have been built so close to each other. This is not the focus of the present paper, but witnesses the cultural relevance of the performed investigations.

In the next section, a presentation of the site is proposed. In Section 3, the GPR investigations are described together with the main results achieved through a traditional data processing. In Section 4, we will show some results achieved with an inverse scattering algorithm applied to the same data. Conclusions will follow.

2. The Archaeological Site of Le Cesine

The waterfront of the northern coastal stretch of “Le Cesine”, a State Nature Reserve on the Adriatic coast near Lecce, is an interesting archaeological zone where imposing ancient remains have been found and mapped also thanks to the use of new detection technologies (Figure 1). This work was carried out in September–October 2020 and in May–June 2021 as part of the UnderwaterMuse Project (Italy-Croatia 2014–2020 Cooperation Programme).

The coast and the shallow waters just off the Reserve have yielded interesting submerged, semi-submerged and terrestrial structures that contribute to a more complete diachronic and contextual understanding of the area.

Recent researches have allowed the identification of further archaeological evidences, which enriches the previous knowledge of the area. In addition to the identification of a wall from the Bronze Age in the “Specchiuddhri” site, it is worth noting the discovery of a submerged port complex presumably of the Roman period, in the locality of Posto San Giovanni, which offers a much more articulated view of the maritime façade of Lupiae (ancient Lecce). Among the evidences already investigated and documented in the 1990s, some more were recognized in the latest photographs and video footage obtained by drone.

Worthy of note, here there is a structure known as the “submerged church” located in the sea off the channel, 150 m from the shore, at a depth up to 5 m, measuring 33 × 15 m. A rocky platform elevated up to 3 m from the bedrock (here, depth 5 m) appears regularly cut and divided into three rooms (the central one is complete, the other ones are only partially preserved due to the erosion of the rock) enclosed in concrete walls preserved up to 1.5 m of depth. Some squared scattered blocks lie close to the structure, and probably penetrate it (Figure 1 and Figure 2, UT 3).

At about 100 m south of the modern drainage channel and 125 m from the shore, at an average depth of 3 m, there is a rectangular structure, which is a foundation already known and documented, composed of blocks of local calcarenite of varying size, arranged parallel and perpendicular to each other, that develops in a rectangular area of 24 × 30 m; the complex is indeed quite larger and conspicuous, as some rows of blocks continuing under the sand demonstrate (Figure 1, UT 4).

In terms of both position and typological and technical characteristics, this structure seems to be connected to the imposing evidence identified during the 2020 surveys. In particular, the latter is a pier’s foundation built using the typical technique of landing stages in the Adriatic and other areas of the Mediterranean Basin (especially eastern Mediterranean), with “caissons”, namely, with two parallel walls built with large parallelepiped blocks of local stone, and with the space between them filled with irregular stones, sometimes reinforced with internal partitions for the distribution of the forces [16,17,18,19].

This structure begins about 15 m away from the current shore and 65 m southeast of the outlet channel; it is bordered by two external sides of large blocks and stretches east/northeast (84° N). It is about 8 m wide and its overall length, in this preliminary investigation phase, seems to be 100 m (83 m visible). The two parallel side alignments are composed of large parallelepipedal blocks (1.50 × 0.65 cm) put along the entire length of the pier (Figure 1 and Figure 2, UT 5). In some places, two or more overlapping rows are preserved, but on both sides many collapsed blocks are scattered across a wide area.

In the terminal stretch, in a secondary position, sections of a pipe dug out of long blocks of limestone, sometimes overturned and no longer aligned (probably also removed by the force of the waves), can be seen. The depth of the blocks at the (assumed) beginning of the pier is smaller than one meter, while at the end it reaches 4.0 m. It is therefore an important marker of sea-level variation, considering a relative rise since the pier’s construction of about 2 m, the draft of ancient ships being compatible with the remaining 2 m. This structure shows an affinity of type and construction technique with the Hadrian’s pier north of the wide bay of San Cataldo, and the two structures are of similarly impressive dimensions and monumental character [15,18,19].

The three submerged structures of San Giovanni were probably interrelated and belonged to an important port complex, whose overall geometry has yet to be determined by means of a targeted study. Certainly, dating to the Augustan age, these port structures form a very appealing hypothesis. In particular, sources mention Octavian’s (the future emperor Augustus, 63 B.C.–14 A.D.) landing at the port of Lupiae on his journey from Apollonia to Rome. Consequently, in the late Republican and early Imperial period, some infrastructures already existed in this area. Only later the port would be moved further north, with the construction of a new large pier wanted by the emperor Hadrian (76–138 A.D.).

Clearly, this imposing complex had to be equipped with service facilities on the seaside; actually, some wall foundations were visible up to the end of the last century, before the construction of some breakwater series built up for protecting the coast. They were immediately to the north of the modern drainage channel and appeared to be perpendicular to each other, defining some rooms (Figure 1 and Figure 2, UT 1).

Finally, visible starting from the southern side of the modern channel, along the seaside and filled with sand, in the bedrock there are numerous elongated excavations, aligned and 30–40 cm deep, long between 2 and 4 m, probably referable to saltpans (small basins for the salt production) (Figure 1, UT 2).

In Figure 2, photographs of some of the most interesting points exposed in Figure 1 are shown, as well as the geographic localization of the archaeological area of Le Cesine.

The present GPR survey aimed at identifying terrestrial evidences now disappeared, in order to redraw the general plan and development of this coastal complex.

3. GPR Surveying and Results

As can be appreciated from Figure 2, three rectangular areas have been prospected, labelled as Area 1, Area 2 and Area 3. They correspond to three rectangles whose size is 15 × 8 m, 15 × 7.5 m and 16 × 12 m, respectively. In each of them, an orthogonal grid of Bscans has been gathered with transect equal to 50 cm. The system exploited for this was a RIS Hi-mode pulsed GPR manufactured by IDSGeoradar equipped with a dual antenna with nominal central frequencies at 200 and 600 MHz.

To avoid pedantic repetitions, here we will present only the results that appeared to us heuristically more interesting, and, in particular, we will show results at 600 MHz with regard to Area 1 and Area 3, whereas we will show results at 200 MHz with respect to Area 2. The “classic” data processing was performed with the Reflexw code and consisted of zero timing, background removal, linear and exponential gain vs. depth, 1D filtering (so as to remove undesired low-frequency peaks in the spectrum of the signal related to the gain vs. depth) and Kirchhoff migration [20]. The propagation velocity exploited for the focusing of the targets and for the time-depth conversion has been evaluated on the basis of the visible diffraction hyperbolas. In particular, we have chosen 11 cm/ns in Area 1 and Area 2, and 12 cm/ns in Area 3. On the basis of the processed Bscans, we evaluated in all the areas a penetration time-depth of about (roughly) 50 ns. Time slices [21,22] have been retrieved so to achieve flat horizontal maps at different depth levels. In Figure 3, the most meaningful results in Area 1 (at 600 MHz) are shown; in Figure 4 the depth levels of largest interest in Area 2 (at 200 MHz) are shown and, finally, in Figure 5, the most interesting results in Area 3 (at 600 MHz) are shown.

With regard to the interpretation of the data, it is to be premixed that we are in a context where the structures (if any) are certainly fragmentary because of the continuous anthropic frequentation and the action of the sea. Hence, we cannot expect precise maps of building following a rectangular pattern, as often happens in ancient Roman sites [23,24]. Moreover, the site does not necessarily correspond to a Roman town but maybe just to a Roman harbour. So, insulated anomalies are outlined in Figure 3, where future localized excavations might and hopefully be performed to test the underground scenario. In Figure 3, we have put into evidence with circles the anomalies showing also some continuation along the depth. Figure 3 also shows a large reflection all over the lower side of the rectangle. As can be appreciated from Figure 2, on that side, the prospected area lies close to a concrete platform, and indeed this platform contains a welded steel mesh (probably on more levels). So, this large anomaly is certainly at least partially related to the reflections from this welded mesh. However, it is also to be outlined that the anomaly appears as large as two meters, which would be unusual if due to the only concrete platform. So, we cannot exclude the presence of some structure beyond the effect of the platform. It is also noticeable, in the upper panel of Figure 3 (i.e., at shallow depth levels), in the central part of the investigated Area, a zone with low reflectivity in the upper panel, i.e., at shallow depth levels, which might also deserve a localized excavation in order to be better understood. In Figure 4 (Area n. 2), also is visible an anomaly on the lower side of the prospected rectangle, and is still related to the same welded mesh as easily understood from Figure 2, but also in this case it is unusually large. The anomaly outlined in Figure 4 resembles instead a possible structure, even if we cannot be sure of this hypothesis. The structure is slightly deeper than those that have been outlined in Figure 3 and in Figure 5, and this is the reason why the data at 200 MHz provide a clearer image with respect to those at 600 MHz (not shown for Area 2). With regard to Area 3 (Figure 5), we have outlined a large anomaly on the left hand side (Area 3 is not close to any modern structure) and an insulated anomaly on the opposite side at the depth level of about 90 cm.

4. Inverse Scattering Results

In this Section, we show homologous slices to those shown in Figure 3, Figure 4 and Figure 5 achieved by means of the processing chain detailed in [7,13]. This processing chain involves time domain filtering procedures and a Born approximation-based MWT approach, whose implementation is made effective by means of the Shift and Zoom (SZ) procedure implemented, as described in [13,14].

The MWT approach provides a focused image of the investigated region by computing the regularized solution of a linear inverse scattering problem. Specifically, the relationship between the GPR data and the unknown targets is expressed in the frequency domain by means of a linear scattering operator [25], which is discretized into a matrix. The targets are expressed by means of the unknown variations of the dielectric permittivity with respect to the assumed background value that describes the investigated medium, and they are accounted for by means of the dielectric contrast function χ. The regularized solution is achieved by means of the truncated singular value decomposition (TSVD) of the relevant operator [26]. In this paper, the useful frequency range is from 100 MHz to 450 MHz for the data collected by the antenna working at the nominal central frequency of 200 MHz, and it is from 250 MHz to 850 MHz for the data collected at the nominal central frequency of 600 MHz. Moreover, the truncation threshold is set all the times in such a way as to filter out all those singular values of the scattering matrix whose amplitude is lower than −30 dB with respect to the maximum one.

The linear inverse-scattering approach is no novelty in itself, but rarely has the problem of the comprehensive size of a real (or realistic) investigation domain for GPR applications been considered in the literature on inverse scattering. Of course, in real GPR case histories, the size of the investigation domain is huge with respect to the central internal wavelength, and coherently also the observation lines (i.e., the Bscans) are much longer than the central wavelength. This makes it computationally impossible to invert the entire physical investigation domain laying under the entire observation line. A first strategy (adopted, e.g., in [11]) can be that of joining side by side several reconstructions according with the schematic of Figure 6 (upper panel). However, in this way, the targets close to the edge between any two adjacent investigation domains will be seen under a reduced-view angle, and so their reconstruction will be necessarily poorer than that of the targets centred with respect to the relative investigation domain and overlying observation line. Moreover, the different average intensity of the reflections between adjacent domains in general produces seam-effects at the vertical bounds between the two domains. These drawbacks are overcome by means of the SZ procedure, where two consecutive investigation domains are not adjacent, but rather superposed to each other. The main steps of the adopted implantation of the SZ procedure are:

• The measurement line O and the investigated domain D are divided into K, partially overlapping subdomains O_k and D_k, k = 1, 2, …, K;
• For each subdomain D_k, the tomographic reconstruction of the contrast function χ_k is obtained by exploiting the TSVD inversion procedure;
• The tomographic image of the overall surveyed domain D is obtained by joining only the central parts of the K reconstructions χ_k achieved for each subdomain D_k, as shown in Figure 6 (lower panel).
• Besides all the central parts of all the sub-reconstructions, also the left hand side of the first sub-reconstruction and the right hand side of the last one are retained and added at the beginning and at the end of the comprehensive image, respectively. This because it is of primary importance not to miss any information retrievable from the data.

A schematic description of the SZ implementation is provided in Figure 6 (lower panel). In particular, in this figure, the green vertical strip is the part retained from the reconstruction of the (k − 1)th investigation domain, the red strip is the homologous quantity relative to the kth investigation domain and the blue strip is the homologous quantity relative to the (k + 1)th investigation domain. It is worth pointing out that, if we adopt the SZ strategy, the investigation domains to be considered are many more than those needed by just joining side by side the sub-reconstructions. However, the computationally time-consuming part of the linear MWT approach is the calculation of the scattering operator, namely, of the elements of the corresponding matrix, and the subsequent numerical evaluation of its SVD. This is done once and for all on one matrix and subsequently applied to all the subdomains except, of course, in the rare case of strong variations of the wave propagation velocity within the length of a single Bscan, in which case different matrixes and SVDs should be calculated. Therefore, customarily, the extra computational burden related to the SZ is not cumbersome and the comprehensive computation burden essentially depends on the length of O_k and the horizontal and vertical size of D_k. Herein, the length of O_k and the horizontal size of D_k are both equal to 2 m, for both 200 MHz and 600 MHz data. Instead, the vertical size of D_k is fixed according to the expected penetration depth and it is 2 m for the 600 MHz data and 4 m for the 200 MHz data. The SZ procedure has been implemented choosing the shifting step equal to 20 cm for all the processed data sets.

Since the classical processing with migration shown in Figure 3, Figure 4 and Figure 5 and inverse scattering algorithms have been applied independently from each other, the two values of the propagation velocity do not correspond exactly (but the two quantities are quite close to each other). Indeed, the MWT approach has been applied by setting the relative permittivity of the investigated medium equal to 7, which means a wave propagation velocity of 11.34 cm/ns, for all the surveyed areas. Moreover, the pre-focusing steps aimed to enhance the signal vs. depth and to filter out disturbing effects were not the same. Therefore, the estimated depth and shape of the main anomalies are not exactly the same. That said, the main anomalies are reproduced by both algorithms and this makes us more confident that they correspond to physical targets, even if we cannot be sure about their nature. In Figure 7, Figure 8 and Figure 9, some slices are shown, where the main anomalies put into evidence are essentially the same as those shown in Figure 3, Figure 4 and Figure 5. The depth of the slices in Figure 7, Figure 8 and Figure 9 are provided directly in cm because, as said, the inversion algorithm works in the frequency domain and its output comes out directly into the spatial domain. Only an excavation will be able to say what reconstructions are closer to the ground truth.

5. Conclusions

In this contribution, we have presented the results of a GPR measurement campaign performed at the archaeological site of Le Cesine, where several remains (part of which are placed underwater) suggest the presence of an ancient Roman harbour. We have investigated three rectangular areas and have provided suggestions for possible future localized test excavations. The data processing has been performed both with a commercial code (focusing the data thanks to a migration algorithm) and with a linear inverse scattering algorithm accompanied by a shift-and-zoom procedure for the inversion of large investigation domains. There is a good correspondence between the traditional processing and the inverse scattering approach (which is a mathematically more refined but also less assessed algorithm), and this guarantees the coherence between the two procedures.

The future developments that we hope to perform are, on the archaeological side, the implementation of at least localized excavations in order to investigate the nature of the anomalies put into evidence by the GPR results.

Moreover, on the side of the GPR data processing, future development might regard the enhancement of the inverse scattering approach to large-scale data, in particular with regard to all those amplification and filtering operations preliminary to the focusing.

Finally, in the future, the strategy of the shifting zoom might be applied also to non-linear approaches, making it possible to expand the comprehensive investigation domain also in this case.