SAR Sentinel 1 Imaging and Detection of Palaeo-Landscape Features in the Mediterranean Area

Abstract

The use of satellite radar in landscape archaeology offers great potential for manifold applications, such as the detection of ancient landscape features and anthropogenic transformations. Compared to optical data, the use and interpretation of radar imaging for archaeological investigations is more complex, due to many reasons including that: (i) ancient landscape features and anthropogenic transformations provide subtle signals, which are (ii) often covered by noise; and, (iii) only detectable in specific soil characteristics, moisture content, vegetation phenomenology, and meteorological parameters. In this paper, we assessed the capability of SAR Sentinel 1 in the imaging and detection of palaeo-landscape features in the Mediterranean area of Tavoliere delle Puglie. For the purpose of our investigations, a significant test site (larger than 200 km²) was selected in the Foggia Province (South of Italy) as this area has been characterized for millennia by human frequentation starting from (at least) the Neolithic. The results from the Sentinel 1 (S-1) data were successfully compared with independent data sets, and the comparison clearly showed an excellent match between the S-1 based outputs and ancient anthropogenic transformations and landscape features.

1. Introduction

For most of the 20th century, aerial photography was the primary remote sensing tool adopted for landscape archaeology and for detecting buried archaeological structures through visual interpretation based on archaeological proxy indicators [1]. The most common archaeological proxy indicators are generally known as crop, soil, shadow, and damp marks and are caused by the presence of buried remains and traces of ancient environments still fossilized in the modern landscape. These features induce spatial anomalies (in vegetation growth and/or status, surface moisture content, and micro-reliefs) that are generally not visible in situ but only evident from above [2]. In recent decades, significant improvements were obtained in the identification of archaeological proxy indicators from very high resolution (VHR) satellite optical and synthetic aperture radar (SAR) data [3].

Starting from this heritage, new applications and developments are expected [4,5,6,7,8,9,10,11], particularly from the use of open data, such as satellite Sentinel 1 and 2 (S-1 and S-2) that are part of Europe’s ambitious Copernicus program [11]. All the Sentinel missions were released under an open data policy to foster knowledge, innovative applications, and advanced developments. In particular, S-1 data were used for the monitoring of cultural heritage (CH) [12], and recent pioneering studies have highlighted that S-1 can be useful to extract information regarding the contemporary and past landscape (the detection of buried or emerging archaeological remains), to infer changes in the current and former environment, and support the discovery of buried archaeological sites. In more detail, early investigations addressed the underwater palaeo-landscape (bathymetric assessment) of the North Sea [13], and land-based archaeological landscape in Egypt to discover potential buried archaeological sites in the agriculture environment [14].

Despite these very promising, early results, the use of S-1 for archaeology is still an open issue and presents current and future challenges. In this context, compared to optical data, the use and interpretation of radar data for archaeological investigations is much more complex for many reasons, in particular as archaeological investigations are focused on the detection of subtle signals (often covered by SAR noise), and only detectable in specific conditions [15]. For these reasons, over the years, SAR based archaeology has been historically limited by the low spatial resolution of the early sensors, even if important archaeological discoveries were made in vast deserted areas, as in the case of the Saharan radar-river [16].

Even if the use of satellite radar in archaeology is still today in its experimental stage, it undoubtedly offers great potential for manifold applications, such as the reconstruction of the palaeo-landscape [16], the detection of buried sites [17], along with the documentation and monitoring of cultural heritage as performed in recent studies conducted using very high-resolution satellite radar data such as TerraSAR-X and COSMO-SkyMed [18,19]. SAR based archaeology poses pressing critical issues, linked with the lack of investigations, particularly those conducted in diverse environments (desert, semiarid, Mediterranean, etc.), different land use/cover types, along with the needs to improve the interpretation and modeling approaches, which may be adjusted or developed ad hoc for different environments.

This study provides a contribution in this context, focusing on the use of S-1 data in the Mediterranean environment to detected ancient landscape features and anthropogenic transformations. To this aim, a significant test area (larger than 200 km²) was selected in the Foggia Province (South of Italy), as this region is characterized by long human frequentation (from the Neolithic time until now) and, consequently, for millennia has been continuously involved in and affected by anthropogenic human transformations. Across and around the study area, large Neolithic villages, pre-Roman routes, and Roman villas have been discovered during the last century.

2. Materials and Methods

2.1. Study Area

The area of interest, known as Tavoliere delle Puglie, lies between the city of Foggia and Lucera in the Apulia Region (Figure 1a). The study area is situated within geographical coordinates between longitude 15°20′30′′ and 15°37′ E, and latitude 41°28′30′′ and 41°35′ N, and covered approximately 224 km² (Figure 1b). The Apulian palaeo-environment coast has been reconstructed by the multilayers of sediment cores drilled in various coastal areas and characterized by different landscapes, in addition to its various ancient civilization sites (e.g., Neolithic, Greek/Roman, and medieval settlements) [20].

From the geological point of view, the investigated area belongs to the central part of the domain of the Apennine foredeep. The deposits formed during several middle and late Quaternary depositional phases were strictly linked with the interplay between the regional uplift and sea-level fluctuations. The oldest terrains, cropping out mainly in the western area, are represented by the argille subappennine unit that form the hills of Lucera (very close to the study area) and consists of a poorly bedded alternation of clays and silty clays of marine origin, lower Pleistocene in age. All the deposits of this unit can be referred to several systems of alluvial fans in seven systems. Terraced deposits are located at different elevations on present-day river beds. Our study area is included in a system, dated to the Upper Pleistocene, named Motta del Lupo, which derives its toponymal from the site of Motta de Lupo, archaeologically known for a Neolithic settlement. The Motta del Lupo system consists of alluvial terraces composed of brow fine sands interbedded with thin layers of pelites. This unit lays in erosive unconformity on the argille subappennine unit and/or on the deposits of older systems. The observed thickness varies from a few meters up to a maximum of 10 m [21].

The Tavoliere is characterized by multi-layered settlements generally referable to Neolithic, Roman, and medieval ages, which make it one of the most important archaeological areas in Southern Italy [22]. The first ancient phase of human settlement is the Neolithic, herein dated between 6500 and 6200 BC (e.g., Puglia region) [23,24], or earlier [25,26]. The Tavoliere includes the area of interest of this study, which remained unknown until the 1940s, when John Bradford recognized Neolithic features from the vegetation marks visible from the RAF and USAAF aerial photographs [27]. Up to 2014, approximately 570 georeferenced sites were found in the areas covered by the original RAF and USAAF air surveys (less than 50% of the total area of the plain), to which we can add 206 known sites that fall outside these mapped areas [28]. Large Neolithic villages, pre-Roman routes, and Roman villas have been discovered during the last century across and around the study area located on the Northeast of Lucera and on the Northwest of Foggia.

The landscape is characterized by the presence of several palaeorivers and palaeochannels, and this is of interest considering that, in the Mediterranean Basin regions, river management generally [29] started in the Neolithic societies during the humid phase of the Holocene, characterized by a humid and rainy climate from approximately 7500 to 4000 BC [21,22,23,24,25,26,27,28,29,30].

Land management in the early stages of Neolithic sedentary agriculture and its impact on the landscape and the organization of settlements is an open issue for archaeologists. In particular, the Neolithic culture of Puglia and southern Italy was one of the first that developed in Europe. The Neolithic culture in Apulia was complex, with settlements built according to their functions: (i) mountain or sub-mountain settlements for pasture (e.g., Murgia area); (ii) near the extraction points and caves (e.g., the coast to the north-east of Foggia); and (iii) in the Tavoliere, along the rivers, for agriculture [31,32,33,34,35,36].

The practice of agriculture at the end of the glaciation period is well documented in the area by archaeologists and archeobotanists and was closely linked to the use of hydrological resources [37,38] (Figure 2). In particular, in [36,39] (Figure 2), the relationship between settlements and rivers is well represented.

Therefore, the detection of palaeo-hydrography is fundamental for the study of the first human settlements, where our knowledge starting from the Neolithic, or earlier, depends on the understanding of the local river dynamics and hydrological variability [29].

2.2. Rationale, Data Set, and Processing

The recognition of archaeological and palaeo-environmental proxy indicators from radar imaging is more complex compared to optical indicators due to a greater number of parameters that characterize SAR data, including (i) surface roughness; (ii) moisture content, linked to dielectrical properties of the target; (iii) the radar system in terms of the operating frequency, polarization, angles, and viewing geometry (ascending or descending); (iv) the characteristics of the sensed surface in terms of the land use and cover types, topography, and relief; along with (v) the morphological characteristics. It is crucial to consider that these features are generally not permanent signals but only “visible” in specific observation conditions depending mainly on the crop phenology, soil type, and moisture content, and in turn, the weather conditions. Therefore, in order to assess the potentiality of Sentinel-1 in the detection of ancient anthropogenic transformations and palaeo-hydrography, a time series for 2014–2019 of Sentinel-1 GRD (IW) was used for the study area to capture both the inter- and intra-year variability expected in the “visibility” of the archaeological proxy indicators. Considering that SAR data are influenced by the radar system itself, in terms of angles, viewing geometry (ascending or descending), and polarization, the data set was selected including both:

✓

ascending and descending modes and;

✓

VV and VH polarizations.

The inter- and intra-year analysis of the data collection, downloaded for free from the ESA [41] and Alaska Satellite Facility (ASF) data [42] websites (

Table 1. Data collection of the study, including the Sentinel 1 sensor type, acquisition date and time, the pass direction, and the source.

Sensor	Acquisition Date	Acquisition Time Sensing Start: Sensing Stop:	Pass Direction	Source
Sentinel 1A	15 October 2014	T16:48:49.539Z T16:49:14.538Z	ASCENDING
Sentinel 1A	21 October 2015	T05:03:12.848Z T05:03:37.848Z	DESCENDING
Sentinel 1A	28 October 2016	T16:48:57.066Z T16:49:22.064Z	ASCENDING	ASF data
Sentinel 1B	17 October 2017	T16:48:27.644Z T16:48:52.643Z	ASCENDING
Sentinel 1A	17 October 2018	T05:03:36.225Z T05:04:01.224Z	DESCENDING
Sentinel 1B	17 March 2019	T16:48:31.045Z T16:48:56.044Z	ASCENDING
Sentinel 1B	29 March 2019	T16:48:31.421Z T16:48:56.419Z	ASCENDING
Sentinel 1B	9 April 2019	T05:02:45.340Z T05:03:10.338Z	DESCENDING
Sentinel 1B	16 April 2019	T16:49:12.556Z T16:49:37.554Z	ASCENDING
Sentinel 1B	15 May 2019	T05:02:46.831Z T05:03:11.830Z	DESCENDING
Sentinel 1A	22 May 2019	T16:49:14.289Z T16:49:39.287Z	ASCENDING
Sentinel 1B	9 June 2019	T16:48:34.549Z T16:48:59.547Z	ASCENDING
Sentinel 1B	21 June 2019	T16:48:35.301Z T16:49:00.299Z	ASCENDING
Sentinel 1A	21 July 2019	T16:49:17.773Z T16:49:42.772Z	ASCENDING	ESA data
Sentinel 1B	26 July 2019	T05:02:57.215Z T05:03:22.214Z	DESCENDING
Sentinel 1A	13 August 2019	T05:03:40.261Z T05:04:05.260Z	DESCENDING
Sentinel 1B	20 August 2019	T16:48:38.937Z T16:49:03.934Z	ASCENDING
Sentinel 1B	12 September 2019	T05:02:59.931Z T05:03:24.930Z	DESCENDING
Sentinel 1A	19 September 2019	T16:49:20.893Z T16:49:45.891Z	ASCENDING
Sentinel 1A	13 October 2019	T16:48:40.700Z T16:49:05.697Z	ASCENDING
Sentinel 1B	19 October 2019	T16:49:21.404Z T16:49:46.402Z	ASCENDING
Sentinel 1A	25 October 2019	T16:48:40.819Z T16:49:05.819Z	ASCENDING
Sentinel 1B	12 November 2019	T16:49:21.367Z T16:49:46.365Z	ASCENDING
Sentinel 1B	24 November 2019	T16:48:40.406Z T16:49:05.405Z	ASCENDING
Sentinel 1A	24 December 2019	T16:49:19.980Z T16:49:44.978Z	ASCENDING
Sentinel 1B	29 December 2019	T05:02:52.669Z T05:03:17.668Z	DESCENDING

), enabled us to identify, for some representative known archaeological buried remains, both the best period of the year and the best SAR based parameters to capture the archaeological proxy indicators (Figure 3).

Radar S-1A and B satellite images were downloaded from the Sentinel Scientific Data Hub and Alaska Satellite Facility (ASF) data. The pre-processing and processing were done using SNAP 6.0 (Italian National Research Council, Tito Scalo, Potenza, Italy) coupled with tools available in ArcMap 10.4.1 (National Authority for Remote Sensing and Space Sciences, Cairo, Egypt) and ENVI 5.1 software (National Authority for Remote Sensing and Space Sciences, Cairo 1564, Egypt).

2.3. Data Processing

In the current study, we tested the Sentinel-1 GRD (IW) acquisition mode’s observations for measuring its suitability in detecting the buried settlements in the study area. In more detail, we used the Apply Orbit File, Thermal Noise Removal, Calibration, Speckle filtering, Range-Doppler Terrain Correction, Conversion to dB, and Filtered band tools to analyze the Sentinel 1 data. Both VV and VH polarizations investigated. The polarisation of SAR imagery is commonly denoted by two letters, the first is the transmitted, and the second is the received [43]. In this study, the VV (vertical transmission; vertical reception) and VH (vertical transmission; horizontal reception) amplitudes were used and projected to ground range using an Earth ellipsoid model, WGS84. The acquired radar data for the study were from Sentinel 1A/B (C-band) interferometric wide (IW) swath mode with the mode products at SAR Level-1 Ground Range, Multi-look, detected (GRD) on 15 October 2014, 21 October 2015, 28 October 2016, 17 October 2017, 17 October 2018, 19 October 2019, and 25 October 2019.

The data processing steps were: Apply Orbit File, Thermal Noise Removal, Calibration, Speckle Filtering, Range Doppler Terrain Correction, and Conversion to dB, along with the Filtered band as follows:

2.3.1. Apply Orbit File

Generally, orbit state vectors, contained within the metadata information of SAR products, are not accurate. The precise orbits of SAR satellites are determined after many days and became available after the generation of the final product by days-to-weeks. SNAP software tools allow the SAR data users to update the orbit state vectors for each scene in its product metadata, provide velocity information, and give an accurate position for the satellite according to the radar tool [44]. Sentinel-specific processing facilities also host the relevant Level-1 instrument processor components [45]. We considered the expected accuracy of “resituated orbits” sufficient for computing interferograms without artifacts [46].

2.3.2. Thermal Noise Removal

In addition to the speckle noise, SAR images suffer from additive thermal noise, especially when the backscattered power is low. Thermal noise removal was used to reduce the noise effects in the inter-sub-swath texture, in particular, to normalize the signal of the backscatter within the entire Sentinel-1 radar image and resulted in reduced discontinuities between sub-swaths for every scene in the multi-swath acquisition modes. The ESA provided thermal noise information for each image in XML formatted files [47]. For removing the thermal noise for the level-1 product, the SNAP provided operator for Sentinel-1 radar data was used to re-introduce the noise signal and allow for re-application of the correction [48].

2.3.3. Calibration

This step includes two main tools; radiometric correction and the multilooking technique. Radiometric correction of an image was conducted so that the pixel values truly represent the backscattering of the reflecting surface. After data pre-processing, the calibrated data were filtered to reduce the inherent SAR speckle noise and clipped in range and azimuth to remove the image border noise. The result of the pre-processing is a re-projected, radiometrically calibrated, and rescaled normalized radar cross-section (NRCS) image [49]. For a particular azimuth time, the thermal noise was estimated in the slant range coordinates by calculating the range spreading loss vector and the elevation beam pattern vector, and applying the scalar contributing factors. According to the Sentinel-1 GRD products, in subtracting the noise from the power-detected image, the thermal noise vectors were converted to ground range coordinates and applied to the data. When calibrating the product to σ, γ, the noise vector must be scaled by the corresponding calibration look-up table (LUT) ($\beta ,\mathsf{\sigma }$, γ, or $DN$, respectively) as in Equations (1)–(4) [49]:

$$noise\left(i\right)=\frac{{\mathsf{\eta }}_i}{A_i^2}$$

(1)

where depending on the LUT selected to calibrate the image data:

$$noise\left(i\right) =\mathrm{calibrated}\mathrm{noise}\mathrm{profile}\mathrm{for}\mathrm{one}\mathrm{of}{\beta }_i^o,{\sigma }_i^o\mathrm{or} {\gamma }_i\mathrm{or}\mathrm{original} {DN}_i$$

(2)

$${\mathsf{\eta }}_i=\mathrm{noiseLut}\left(\mathrm{i}\right)$$

(3)

$$A_{i }=\mathrm{one}\mathrm{of}\mathrm{betaNought}\left(\mathrm{i}\right),\mathrm{sigmaNought}\left(\mathrm{i}\right),\mathrm{gamma}\left(\mathrm{i}\right),\mathrm{dn}\left(\mathrm{i}\right)$$

(4)

After obtaining the calibrated noise profile, removing the noise from the GRD data became available by subtraction. During the TOPSAR sub-swath merging, the noise vectors for any pixel ί that fell between points in the LUT value were also merged into one vector by merging two adjacent ground range images [50].

The calibrated noise profile can be applied to remove the noise by subtraction. Application of the radiometric calibration LUT and the calibrated de-noise LUT can be applied in one step as in Equation (5) [51]:

$$value\left(i\right)=\frac{\left(D{N}_i^2-{\mathsf{\eta }}_i\right)}{A_i^2}$$

(5)

After the radiometric calibration, the spatial resolution is degraded but the image noise is reduced and approximate square pixel spacing is achieved. Therefore, the Multilooking technique as a tool in the SNAP software is required to reduce the speckle noise effect, reaching a spatial resolution of 20 m [52].

2.3.4. Speckle Filter

Speckle filtering is needed to suppress the noise and to remove observations that are not affected by noise and contain valuable land surface information. To remove the speckle, the refined Lee filter was applied as it maintains the detail of the standing boundary [53].

2.3.5. Range-Doppler Terrain Correction

Geometric correction using the “Range-Doppler Terrain Correction (RDTC)” tool in the SNAP software was used for converting the Sentinel-1 GRD data from ground-range geometry into a map coordinate system [54]. In more detail, terrain correction was applied to transfer the single 2D of raster radar geometry to accurately geocode the images by correcting SAR geometric distortions using digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) (that allowed us to consider the local elevation variations) for geocoding Sentinel 1 imagery [55]. Step one was choosing the imported SAR data as an input in orthorectification and the directory output to save the orthorectified image. Both the input and output data were saved in one projected file. The second step will define the parameters and bands amplitude VV and VH that will be processed using Shuttle Radar Topography Mission (SRTM 3 Sec) as inputs to the digital elevation model (DEM) data. This was automatically downloaded, as both of DEM and Image Resampling use the bilinear interpolation method. We defined pixel spacing to 30 m, the corrected image was resampled to 30 m from 10 m, and thus this new size fitted both the used map projection, the orthorectified SAR data size in the GLS-2000, UTM, and also the datum, WGS 1984. The areas without elevation were masked based on DEM data. Then, the process was running automatically after defining the parameters [54].

2.3.6. Conversion to dB

This step was applied to change the “Sigma0_VH_db” and “Sigma0_VV_db” bands from the virtual bands to files. In this conversion, the unitless backscatter coefficient was converted to dB using a logarithmic transformation [52].

2.3.7. Filtered Band

Filtration tools are included in SNAP software that can be used in the processing steps for the improvement of the analyzed data. These include detect lines (e.g., horizontal edges, vertical edges, compass edge detector, etc.), detect gradient (e.g., emboss), smooth and blur, sharpen, enhance discontinuities, and nonlinear filters. There are also morphological filter tools (e.g., erosion, dilation, opening, and closing) that can be used in the morphological field.

2.3.8. Time Series Analysis for Extracting the Backscatter Values

The SAR data were chosen to make the comparison between the backscatter values for three classes; Neolithic, palaeo-hydrography (palaeoriver and palaeochannel), and agriculture (non-archaeological sites). The data were clipped for all the dates of October in 2014, 2015, 2016, 2017, 2018, and 2019 to test the best dates for obtaining the potential archaeological and palaeo-hydrographic features according to the backscatter values of the chosen points. The time series analysis was tested for ten months in 2019 to cover the four seasons in March, April, May, June, July, August, September, October, November, and December. One batch was built to apply on all tested data with calibrate, terrain correction, speckle filter, and liner to dB. Two samples were measured for every class, and the median was calculated cross the VV and VH polarizations for both pass directions (ascending and descending) according to Equation (6):

$$AC\left(x\right)={\left(x^{\prime }-\mu \right)}^T{c}^{-1} \left(x^{\prime }-\mu \right)$$

(6)

where $x$ is the image pixel, $x’$ is the vector that formed by the band values of the image, the $x$, $\mu$ is a vector composed of the mean value of the background pixels in each image band, and $C$ is the covariance matrix of the image bands (that computed cross the background pixels) [56,57].

2.3.9. Computation of Polarimetric Indicators

The ratio of σ° VV/VH, which is the ratio between σ° VH and σ° VV.

3. Results

In this section, we will discuss the results obtained from the analysis of the whole time series of 2014–2019 made to characterize the (i) intra- and inter-year variability in the visibility of palaeo-hydrographic features and archaeological proxy indicators (as explained in Section 2.2).

In order to make the interpretation of S-1 data easier, the SAR data set related to 2019 were compared with the 2019 multi-date S-2-data set elaborated as detailed in [58] (Figure 4).

It is important to consider that the whole area under investigation is classified as arable land in the Corine land cover [59] and devoted to agricultural activities, mainly cereal crops. So that, in the study area, the October month is a period before snow, and therefore, expected to be mainly characterized by the presence of bare soil and spontaneous herbaceous cover (see Figure 4b) as evident from the S-2 RGB false color where the green areas are related to vegetation (mainly herbaceous cover) whereas the brown areas are mainly related to bare soil. These surface conditions can be considered ideal for the detection of soil moisture anomalies linked to archaeological proxy indicators, such as soil and damp marks. For SAR sensors operating at the C band as S-1, we have to consider the interfering effects that are introduced by surface roughness.

Sentinel-1 for the Identification of Archaeological Proxy Indicators and Palaeo-Landescape Elements

SAR analysis conducted in the study area using Sentinel-1 data allowed us to assess the potential of this tool, thanks to the integration of data obtained with multispectral/multitemporal images obtained from Sentinel-2 [58] and ancillary data from bibliography and other sources.

SAR radar backscatter measurements were influenced by both the terrain structure and surface roughness, and it is expected that the more roughness, the greater the backscatter (resulting in a bright feature). The diverse polarizations, in our case, VV and VH helped in discriminating and estimating the different contributions due to (i) the moisture content and (ii) roughness.

In the final analyzed S-1-A and B images (October 2014, 2015, 2016, 2017, 2018, and 2019), some black ditches and white complex areas were identified with archaeological potential for buried archaeological remains in the study area (Figure 5a–f).

We compared the SAR based anomalies of potential archaeological interest with the results from (i) previous investigations and ancillary data as the aerial photographs acquired by RAF and USAAF air surveys photographs and (ii) published researches [60,61,62,63]. The results from these comparisons, along with the analyses of the shapes and size of the discovered archaeological features suggested that they were related to the Neolithic period. Due to the noise present in the Sentinel-1 data, the first activity was to identify known archaeological complexes so as to understand the type of result provided by the S-1 SAR sensor (Figure 5). The features in the S-1 analyzed images (Figure 5c,f) had the typical shape of the Neolithic settlements of the Tavoliere as circular with concentric ditches that enclosed smaller ditches and huts [36] (Figure 6).

These detected shapes are supposed to date back to the Neolithic Era according to the findings of previous works [30,62,63,64]. The analysis allowed us to identify SAR data sites already known in the bibliography, such as those east of Lucera as: (i) Masseria Sarcone and (ii) Masseria Villano (Figure 7). Within the S-1 SAR data, the outer ditches of the settlements and the smaller inner ditches are visible, according to Neolithic settlement shape (Figure 6).

The results obtained from the diverse acquisitions of S-1 showed how the visibility of features of archaeological interest and other elements useful for the reconstruction of the ancient landscape are influenced by factors, such as: (i) the weather; (ii) soil moisture; (iii) temperature; and (iv) soil roughness (see Discussion). In the present study, the images that gave the best results were the ones produced by the acquisition of 19 October 2019.

In addition, the ability to identify macro-elements of the ancient landscape was used on the comparison of known data present in previous archaeological and remote sensing studies applied to the province of Foggia [36,37,38,39,40,66,67,68,69] and S-1 data. The SAR data were used with S-2 data to trace the network of palaeo-rivers, communication routes, and Neolithic settlements, using as a test-area the area north-west of Foggia, known for the presence of ancient remains [66] (Figure 8).

Comparing S-1 and S-2 data, it is possible to clearly observe the palaeo-hydrographic features from S-1 (Figure 8a), due to changes in the moisture content, and the same palaeo-hydrographic features from the soil map (Figure 8b) obtained from S-2 [58] due to changes in moisture content and also probably due to soil organic nutrients.

The importance of an integrated approach of S-1 and S-2 for archaeology has been seen in the comparison of Sentinel data with high-resolution data on Google Earth (Figure 9).

For the same scene, a very high-resolution image from Google Earth (Figure 9d) acquired in the same month (October), but different year (2014), did not exhibit any indicators of palaeo-environmental and archaeological features.

The comparison between the S-1 and S-2 scenes, evidence generally in correspondence of dark marks from S-2, corresponded with clear tone marks from S-1 (Figure 9I–III). However, in some cases, the tone of the features were inverted, appearing with light grey from S-2, and dark grey from S-1 (Figure 9IV).

In other cases (Figure 10a,b), the brightest signal was not due to the presence of palaeo-rivers but likely due to the roughness of the soil. However, the integration of the Sentinel’s data with Google Earth data (Figure 10c) has allowed us to discover that the curvilinear palaeo-channel (characterized by two branches) goes around a Neolithic site, which is characterized by the circular features typical of the Neolithic settlements in this area.

The behavior of the SAR signal in the identification of the palaeo-channel has proven to change depending on diverse factors. In addition to those already mentioned, the penetration capability and the depth of the buried palaeo-channel played a fundamental role. Likely due to the depth that is greater than the expected penetration capability, that the C-band waves were completely absorbed and the backscatter was very low: similar to the “radar river” behavior observed in deserted areas [67].

This behavior was observed for the area of Motta la Regina (Figure 8III and Figure 11). Contrary to what is generally observed for other sites (Figure 8I,II, Figure 9 and Figure 10), in these images, the palaeo-channels were characterized by a darker grey tone.

The great anthropogenic impact of the Neolithic trenched sites on the territory left visible traces, due to the abundant storage of moisture inside trenches and ditches, which make it possible to hypothesize the existence of new archaeological sites, also in accordance with the bibliography about the relationship between Neolithic sites, communication routes, and water resources (Figure 8).

The area represented in Figure 12 is located north of Foggia and is very close to the ancient city of Arpi. According to previous observations on already known archaeological sites (Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10), the data suggest the presence of (i) Neolithic settlements in diverse locations, characterized by large circular or subcircular ditches; and (ii) roads and canals [23,24,25,26,27,28,29]. These elements are visible in the comparison of Figure 12b–d. The PCA from [58] mainly highlights the features left by roads and palaeo-channels. The same features are visible in the S-1 images more than in the RGB of 20 October 2019. The capability of SAR data to highlight elements of the ancient landscape, such as roads and rivers, is certainly useful in the identification of archaeological sites. In fact, settlements are often located along communication routes or water-courses for easy access to resources. As shown in Figure 12b–d, III, in the center of the image, ditches and streets are clearly visible in both the S-1-2 and S-1 images. In this case, there is a very articulated road system, which develops along two directions: a south-north road and an east-west road that crosses a trenched site, visible from the dark outline.

Features in Figure 12, evident in both the VH and VV SAR polarizations, are due to the combined effects of moisture and penetration capability of S-1 that are present and conditioned by the soil grain size and, above all, soil moisture content. In particular, concerning specifically the soil moisture, results from studies [64] conducted in several geographic areas, showed that VV enables better estimation of soil moisture compared with VH.

4. Discussion

The analysis of S-1 SAR data, according to what is described in the archaeological bibliography, aerial archaeology research, and remote sensing studies in the Foggia [36,37,38,39,40,66,67,68,69] area, proved useful for the identification of features of archaeological interest and for the identification of elements of the ancient landscape. However, as shown in Figure 6 and the following, the results are not always optimal for archaeological research and are closely related to several factors including (i) the type of buried evidence; (ii) sensor acquisition mode; and (iii) weather and soil conditions.

Buried remains or elements with a high impact on the landscape (rivers, roads, ditches, etc.) are more visible when they tend to retain soil moisture and are large enough to be appreciated with the use of S-1 SAR and its spatial resolution.

However, other factors deserve further discussion. In Section 4.1 and Section 4.2, the results are discussed in relation to the acquisition mode (ascending or descending orbit of the satellite) and the meteorological data that influenced the identification of the elements of archaeological interest within the multitemporal series of images.

4.1. Identification the Impact of Ascending/Descending Acquisition Modes

To better understand the intra-year variations in the visibility of the archaeological proxy indicators, the SAR data were first analyzed for the whole 2019 year (

Table A2. Backscatter values for the Sentinel 1 2019 Sigma0-VH ascending pass direction.

Date	Agriculture 1	Agriculture 2	Median	Palaeo-river 1	Palaeo-river 2	Median	Neolithic 1	Neolithic 2	Median
17 March	−18.0	−17.0	−17.5	−17.0	−20.0	−18.5	−19.0	−22.0	−20.5
29 March	−18.0	−18.0	−18.0	−19.0	−19.0	−19.0	−20.0	−24.0	−22.0
16 April	−15.0	−19.0	−17.0	−13.0	−19.0	−16.0	−14.0	−18.0	−16.0
22 May	−17.0	18.0	−17.5	−11.0	−15.0	−13.0	16.0	15.0	−15.5
9 June	−15.0	−19.0	−17.0	−13.0	−17.0.0	−15.0	−20.0	−16.0	−18.0
21 June	−15.0	−23.0	−19.0	−18.0	−23.0	−20.5	−23.0	−16.0	−19.5
21 July	−17.0	−20.0	−18.5	−17.0	−21.0	−19.0	−19.0	−17.0	−18.5
20 August	−19.0	−24.0	−21.5	−17.0	−23.0	−20.0	−18.0	−14.0	−16.0
19 September	−17.0	−20.0	−18.5	−19.0	−21.0	−20.0	−20.0	−20.0	−20.0
19 October	Nan	Nan	Nan	Nan	Nan	Nan	Nan	Nan	Nan
25 October	−19.0	−20.0	−19.5	−17.0	−18.0	−17.5	−20.0	−21.0	−20.5
12 November	−11.0	−15.0	−13.0	−15.0	−14.0	−14.5	−19.0	−15.0	−17.0
24 November	−13.0	−13.0	−13.0	−16.0	−13.0	−14.5	−18.0	−14.0	−16.0
24 December	−18.0	−21.0	−19.5	−18.0	−19.0	−18.5	−18.0	−22.0	−20.0

,

Table A3. Backscatter values for the Sentinel 1 2019 Sigma0-VV descending pass direction.

Date	Agriculture 1	Agriculture 2	Median	Palaeo-river 1	Palaeo-river 2	Median	Neolithic 1	Neolithic 2	Median
9 April	−8.0	−16.0	−12.0	−10.0	−15.0	−12.5	−14.0	−9.0	−11.5
15 May	−10.0	−13.0	−11.5	−12.0	−13.0	−12.5	−14.0	−10.0	−12.0
26 July	−13.0	−16.0	−14.5	−12.0	−15.0	−13.5	−11.0	−10.0	−10.5
13 August	−15.0	−14.0	−14.5	−13.0	−13.0	−13.0	−13.0	−11.0	−12.0
12 September	−10.0	−15.0	−12.5	−13.0	−16.0	−14.5	−12.0	−12.0	−12.0
29 December	−11.0	−13.0	−12.0	−7.0	−8.0	−7.5	−8.0	−9.0	−8.5

,

Table A4. Backscatter values for the Sentinel 1 2019 Sigma0-VH descending pass direction.

Date	Agriculture 1	Agriculture 2	Median	Palaeo-river 1	Palaeo-river 2	Median	Neolithic 1	Neolithic 2	Median
9 April	−17.0	−19.0	−18.0	−17.0	−21.0	−19.0	−15.0	−22.0	−18.5
15 May	−17.0	−18.0	−17.5	−17.0	−19.0	−18.0	−17.0	−18.0	−17.5
26 July	−20.0	−23.0	−21.5	−19.0	−25.0	−22.0	−21.0	−14.0	−18.0
13 August	−22.0	−24.0	−23.0	−19.0	−22.0	−20.5	−20.0	−15.0	−17.5
12 September	−19.0	−24.0	−21.5	−18.0	−24.0	−21.0	−23.0	−22.0	−22.5
29 December	−19.0	−22.0	−20.5	−18.0	−19.0	−18.5	−18.0	−21.0	−19.5

and

Table A5. Backscatter values for the Sentinel 1 2014, 2015, 2016, 2017, 2018, and 2019 Sigma0-VV polarization.

Date	Agriculture 1	Agriculture 2	Median	Palaeo-river 1	Palaeo-river 2	Median	Neolithic 1	Neolithic 2	Median
15 October 2014	−7.0	−10.0	−8.5	−10.0	−8.0	−9.0	−9.0	−11.0	−10.0
21 October 2015	−5.0	−9.0	−7.0	−10.0	−7.0	−8.5	−6.0	−7.0	−6.5
22 October 2016	−7.0	−9.0	−8.0	−10.0	−10.0	−10.0	−11.0	−8.0	−9.5
28 October 2016	−7.0	−10.0	−8.5	−9.0	−7.0	−8.0	−8.0	−7.0	−7.5
17 October 2017	−9.0	−10.0	−9.5	−14.0	−11.0	−12.5	−12.0	−9.0	−10.5
17 October 2018	−9.0	−10.0	−9.5	−9.0	−9.0	−9.0	−10.0	−11.0	−10.5
19 October 2019	Nan	Nan	Nan	Nan	Nan	Nan	Nan	Nan	Nan
25 October 2019	−6.0	−8.0	−7.0	−10.0	−9.0	−9.5	−17.0	−15.0	−16.0

), considering, along with the VV/VH polarizations, the impact of ascending/descending acquisition modes over features related to known areas, such as agricultural areas, palaeo-river, and buried Neolithic settlements, as shown in Figure 13a–d. If we focalized on the single date acquisition, it was clear that according to the period of the year, both ascending and descending acquisitions enabled the discrimination of the diverse targets, and, if we analyzed one by one the single date acquisitions, the ascending tracks could be preferred because they are characterized by the higher backscattered coefficients that are due to the lowest incidence angle. As a whole, the best discrimination of the diverse considered targets is April and October, as evident from (Figure 13a–f) and also from the comparison made between the 2019 multi-dates S-1 and S-2 data set.

In the same context, SAR data for 2014, 2015, 2016, 2017, 2019, and 2019 were analyzed along with the VV/VH polarizations for the month of October (see Table in Appendix A and Figure 13f,g), which is considered the best time for extracting some buried shapes (Figure 13a,e,f). As a whole, the analysis performed herein suggested that:

✓

First, the best periods for the visibility of the archaeological proxy indicator were summer and autumn, as also suggested from other previous studies [68,69,70]. Autumn, in particular for SAR imaging, enabled a better capture of the soil moisture in bare soil;

✓

Secondly, considering that radar imaging outputs are influenced by the radar system itself, in terms of angles, viewing geometry (ascending or descending), and polarization,

○

both ascending and descending modes appeared acceptable; however, we selected the ascending mode acquisition, as it is characterized by higher backscattering values due to the lower view angle compared to the descending pass, and;

○

both VV and VH polarizations were suitable to enhance soil moisture and, in turn, archaeological marks; however, VV appeared more suitable.

4.2. Sentinel 1 Multiyear (2014–2019) Data Analyses

Outputs from the multiyear (2014–2019) data analysis highlighted that, in the study area, the highest capability of SAR Sentinel 1 in capturing archaeological proxy indicators was observed during summer and autumn. This is linked to the characteristics of the area in terms of the vegetation types and crop phenology, along with the impact of the meteorological conditions on the soil moisture content. To consider the contribution of soil moisture, we analyzed the SAR data coupled with the data on rain precipitation acquired from the meteorological stations in the areas and available from the Protezione Civile of the Apulia Region [71] (Figure 14a,b).

A comparative analysis of Figure 13 and Figure 14 suggests that the reason why the most appropriate months (July to November) obtained good results in S-1 was linked with the long duration of drying (with little of rains) and the higher SAR backscatter related to the agriculture land. In the study area, mainly cultivated with cereals (grain), the season of summer and autumn were characterized by a low presence of vegetation, which made the visibility of the soil and damp marks easier; thus, the best period to capture these archaeological marks is the month of October.

5. Conclusions

This study aimed to test the capability of Sentinel-1 Radar data in the identification and imaging of the subtle features linked to the palaeo-landscape and ancient human-induced transformation of the landscape in the Mediterranean environment.

The SAR data were difficult to understand and certainly less immediate than the S-2 optical sensor, especially in the detection of archaeological sites. In order to understand the SAR data of S-1, the results of the analysis were compared with information on the archaeological sites already known in the literature [36,37,38,39,40,66,67,68,69]. This comparison allowed a process of learning and understanding of the S-1 SAR data and the simultaneous use of S-2 and ancillary data to support and validate the diverse hypotheses. Of course, this approach was made possible by the rich archaeological documentation already produced for the area of study from research institutions and universities (see Section 2.1 and Section 3). The identification of sites from validated data, as in the case of the linear and circular features visible in Figure 7 [65] or Figure 11 [66], allowed us to understand which features were related to Neolithic settlements and to hypothesize, based on this experience, the existence of new archaeological sites as in Figure 12.

However, the visibility of the proxy indicators related to the palaeo-environmental and archaeological features, as described, depends on a number of factors, including (i) the presence of micro-topographical features; (ii) changes in the moisture content; and (iii) variation of crop growth and soil nutrients. These factors vary from one season to another and are dependent on the meteorological conditions, land use, and pedological and geological characteristics [2,21].

For these reasons, investigations based on time series are important to analyze and characterize the intra- and inter-year variability in the visibility of subtle features linked to the palaeo-landscape. To this aim, a significant test area (larger than 200 km²) was selected in the Foggia Province (South of Italy) and investigated using a time series of 2014–2019 for S-1 data to analyze and characterize the (i) intra- and inter-year variability in the visibility of features linked to the palaeo-landscape, along with the (ii) impact of different acquisitions (ascending and descending) modes and polarization VV/VH.

As a whole, the results from these investigations indicated that (i) the best period for the visibility of the archaeological proxy indicator was summer and autumn as already determined for optical data [68]; (ii) both ascending and descending modes well captured the subtle features of interest; however, the ascending mode was preferable as it was characterized by higher backscattering values. In more detail, on the basis of our analysis, as expected being that archaeological features are not permanent signals, several specific acquisitions from the S-1A, B data (captured on 15 October 2014, 21 October 2015, 28 October 2016, 17 October 2017, 17 October 2018, 19 October 2019, and 25 October 2019) provided the highest visibility of the subtle palaeo-landscape features. On the basis of these S-1 acquisitions, several new anomalies of potential archaeological interest were identified. In particular, palaeo-rivers and palaeo-channels along with ditches in shape and size “compatible” with the Neolithic era were detected and successfully compared with independent studies and data (such as the air photographs acquired by the RAF during the Second World War).

Insights into the strong connection between Neolithic settlements, roads, and palaeo-rivers/palaeo-channels were evident and highlighted by our S-1 based findings (see Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12). This is very important, not only for the study areas but also because this opens new potential S-1 applications and research lines. In fact, the results from our investigations clearly highlighted the high potentiality of S-1 in the identification of the very subtle features linked to palaeo-landscape features and ancient anthropogenic transformation in vegetated areas, such as in the Mediterranean ecosystems. This is a very important result, considering that, since the early assessment of the potentiality of SAR in archaeology, satellite L-band SAR sensors SAR data were traditionally used for subsurface imaging in arid regions to reconstruct the palaeo-landscape, given the ability of longer wavelengths to penetrate more deeply into sand [72]. This study demonstrated the high capabilities of shorter-wavelengths (in particular the S-1 C-band sensor) in subsurface imaging in vegetated areas of the Mediterranean ecosystems, and these results offer important new insights in the application of SAR for palaeo-environmental investigations:

✓

First, the capability of S-1 to detect traces of past environment and landscape with particular reference to ancient surface water rivers highlighted that can S-1 data can provide a major contribution in archaeological investigations considering that rivers are and were crucial to past human activity and are, therefore, considered important targets of archaeological prospection;

✓

Secondly, the successful results of this investigation can be replicated in different geographic areas considering the free worldwide availability of S-1 data along with data processing tools;

✓

Thirdly, considering that the main critical, challenging aspect of the use of SAR in archaeology is a lack of correspondence between the great amount of spaceborne SAR data (as in the case of S-1) and effective methods to extract information linked to past human activity, this paper contributes in providing a methodological approach to exploit S-1 data in archaeological investigations.

As a whole, the results from our investigations highlighted that S-1 data can provide a major contribution in archaeological investigations and also overcome the limits of passive optical data: active sensors are able, to some extent, to “penetrate” vegetation and soil and unveil important targets of archaeological prospection, such as palaeo-rivers and palaeo-channels, that are currently not evident in the modern landscape and from the optical data. The detection of palaeo-hydrography is fundamental for the study of human history and the first human settlements, in the Neolithic, or earlier, which depends on the understanding of the local river dynamics and hydrological variability [36,37,38,39,40,66,67,68,69,73,74,75,76]. Palaeo-hydrography can inform us regarding the ancient site’s occupation in the prehistoric period as water is, and was, the most immediate requisite element for human survival. Therefore, the identification of palaeo-rivers and palaeo-channels, along with ancient roads, can suitably provide information regarding buried lost settlements, which generally require the data of large areas to investigate.

Once proven that S-1 can be used for the identification of remains of archaeological interest and palaeo-environmental elements, there are many possible future developments, some of which have been already successfully tested with other types of data (e.g., TerraSAR-X, COSMO-SkyMed, or Gaofen-1), such as the identification, the extraction, or the classification of linear and sub-linear features [18,19,77,78,79].

The approach proposed herein provides a great advantage to easily identify large features (such as palaeo-channels and ditches) or linear routes, and to narrow down surface surveys to be further investigated using very high-resolution data and/or geophysical prospections, thus reducing the time of analysis and, potentially, the cost for excavations. This type of approach provides a tool, in addition to those already used in the Earth Observation techniques and in combination with some of them (e.g., S-2, Landsat TM, or PRISMA Satellite data), at the service of landscape archaeology, regardless of the area of application. On the other hand, knowledge of the territory, anywhere in the world, opens up opportunities for the protection, discovery, and preservation of cultural heritage.