On the Combination of Remote Sensing and Geophysical Methods for the Digitalization of the San Lázaro Middle Paleolithic Rock Shelter (Segovia, Central Iberia, Spain)

Abstract

This paper is focused on the Middle Paleolithic rock shelter called “Abrigo de San Lázaro”, placed in the Eresma River valley (Segovia, Spain). In this area, a multisource geomatic approach is used. On the one hand, the external envelope of the shelter has been digitalized by the means of an efficient combination between aerial photogrammetry and laser scanning (static and mobile). On the other hand, the ground penetrating radar and the electric tomography were used with the aim of evaluating the inner disposition of the shelter. The combination of both digitalization (external and internal) has allowed for improving the knowledge of the site characteristics that, in turn, will facilitate the future excavation works. The results of these studies allow archaeologists to know new data for a better understanding of the site formation (geology of the site, sedimentary potential, rock shelter dimensions, etc.) and the events that took place in it (knowing its historical evolution, especially the interaction between man and the environment). Additionally, the information obtained from these studies is very useful to plan future excavation works on the site.

1. Introduction

The integration of remote sensing and geophysical methodologies in archaeology has allowed advancements in the characterization of the external and internal parts of complex settlements. As a result, the integrations of these methodologies could provide two main benefits: (i) optimization of the fieldworks (archaeological excavations and surveys), (ii) improvement of the accuracy analysis and (iii) documentation methods. Last but not least, this multidisciplinary work also enhances the methodologies currently used in archaeological fieldworks, generating new procedures and thus new valuable information of the archaeological site.

The use of remote sensing techniques was introduced in the last decades in archaeology with the aim of improving the documentation processes [1,2]. The main remote sensing techniques applied in archaeology are image-based methods (photogrammetry or image based modelling) [3,4,5] and range-based methods (laser scanner) [6,7,8] (Table 1).

-

Image-based methods: this group includes terrestrial and aerial photogrammetry (on board airborne systems like drones, among others). These methods use the image based modelling strategy, more specifically the Structure from Motion (SfM) approach, making a 3D reconstruction of the scene by the use of digital images (2D information). Images must follow a certain shooting protocol to obtain a rigorous geometric model. Topography support (GPS, Total Station, rulers, etc.) is also required for the scale definition [9,10].

-

Range-based methods: these methods are based on active sensors that emit a laser beam, which is reflected from the object or scene to capture metric information. The final product only includes spatial information of the scene with a unique radiometric value that depends on the laser wavelength [11].

With these methods, it is possible to transform places and real objects into a digital format by means of a 3D points cloud. These clouds are constituted by points in space representing the captured geometry, with metric and textural rigor in a precise and detailed way [1,2,3,4,5,6,7,8]. By the use of aerial photogrammetry or airborne Light Detection And Ranging (LiDAR), archeological sites in open areas can be detected and characterized in a global or general way [12,13]. On the contrary, terrestrial photogrammetry or terrestrial laser scanning can be used in more local environments for both open and closed spaces [1,2,3,4,5,6,7,8]. A new collection of handheld sensors called Backpack Wearable Mobile LiDAR Systems (WMLS) are currently emerging. The resolution of these systems is worse, for Terrestrial Laser Scanner (TLS) devices, around 8 mm in the best case. However, these dynamic techniques can be applied to record scenarios where the terrestrial laser scanner (TLS) or terrestrial photogrammetry could not be efficient due to its complexity and the number of stations required [14,15,16,17]. In addition, these WMLS techniques are mainly interesting in the documentation of closed spaces with narrow paths such as caves, rock shelters or underground dolmens, etc. They can also be used as a complement for the airborne LiDAR or a photogrammetric flight in reservoirs located in forested areas, since they provide much more points density per square meter at much more affordable costs [18]. The WMLS save time and money in the capture and processing of data (up to 10 times faster than static laser systems) (17). On the contrary, as it has already been highlighted, their accuracies are worse (centimetric) and they require to perform specific protocols (i.e., planning of routes that start and end the scan in the same position). Fortunately, the development of Simultaneous Localization and Mapping (SLAM) techniques is allowing guaranteeing acceptable precisions [19].

The application of these methodologies is useful in several of the processes involved in an archaeological site study (Table 2).

As a result of the previous motivation, we can affirm that geotechnologies are necessary tools for the external geometric characterization of a heritage resource, enabling a more precise archaeological investigation. Thus, they are capable of carrying out studies and analysis with greater accuracy, making wide dissemination of the results according to the needs of the researcher, in a more graphic and visual way [7,39,40].

On the other hand, geophysics is the discipline that includes numerous techniques dedicated to the study of the subsoil elements. This discipline is commonly used in archaeology in the first phase of surveying, although it is also used during the different phases of the excavation process [41,42]. The main aim of this discipline in archaeology is the internal geometric characterization (subsoil) of an archaeological site. Within the techniques that constitute geophysics, this research is focused on two of them: (i) the Ground Penetrating Radar (GPR) [43,44] and; (ii) the Electrical Resistivity Tomography (ERT) [45]. Both techniques are non-invasive. GPR and ERT allow a first approximation of the characteristics of the elements located in the subsoil at a specific site (if the study case is favorable), allowing archaeologists to plan an excavation campaign, focusing interest on those areas where the best results are hoped [42]. The application of these techniques has meant a revolution for archaeologists, since they allow the optimization of resources (which in many cases are limited). Thus, they are able to investigate in the area they want to, digging in specific areas of a site (the site sometimes occupies several hectares) and avoiding unnecessary damages to the site (in many cases the archaeological remains are under current constructions, such as asphalt, buildings, gardens, etc.).

For the location of specific elements, such as tombs or artifacts, the GPR technique presents very precise results [46,47]. It is also useful for a more extensive location of construction elements (e.g., buildings, warehouses, etc.), engineering works (e.g., roadways, water pipes, aqueducts, etc.), and geological studies (e.g., location of cavities, stratigraphic study, etc.) [48,49,50,51,52]. Both ERT and GPR are effective techniques if the chosen method is correct. It is important to note that there are several factors to consider before choosing one or another technique; GPR allows limited prospecting depths and is specially designed for electricity resistive grounds, whereas ERT allows for much higher prospecting depths and a wider range of use. For both systems, the presence of different electric conductive layers favors the results reliability.

In this paper, a multidisciplinary methodological approach towards the external and internal geometric characterization of an archaeological site is presented based on the fusion of remote sensing and geophysical techniques. Within this context, the paper has been organized as follows: Section 2 introduces the study case and defines the materials and methods used to digitalize it; Section 3 describes the experimental results obtained, Section 4 exposes the discussion and finally some conclusions obtained from the present work are drawn within Section 5.

2. Materials and Methods 2.1. The San Lázaro Rock Shelter (Segovia, Spain)

2.1.1. Geographical and Geological Context

The San Lázaro rock shelter is a small cavity originated in the cretaceous dolomites of the left bank of the Eresma river (Segovia, Central System, Central Iberia) and about 15 m above the riverbed (Figure 1). The site, discovered in 2014, takes its name from the medieval bridge that crosses the Eresma river a few meters far from the rock shelter [53]. In turn, the San Lázaro rock shelter is located about 400 m upstream of the Eresma river, and the Molino rock shelter, a Mousterian site excavated between 2013 and 2018 by the same research team [54,55,56,57]. The site is placed in the upper part of the Eresma river valley, just passing through the city of Segovia.

The upper part of the Eresma river encompasses the area of its birth on the summits of the central range and its course on the northern slopes of the mountains and highland foothills [58,59]. This section of the valley is characterized for having carved its valley in the proterozoic-paleozoic igneous-metamorphic materials of the mountains and crystalline piedmont (gneisses, migmatites, marbles, granitoids), corresponding to the Central Iberian area of the Iberian Massif—this part where the river draws gorges and open defiles or valleys of tectonic origin (Valsaín), and coming from the incision in the sedimentary detrital-carbonate materials of the Mesozoic border (Upper Cretaceous) corresponding to the western sector of the Castilian Branch of the Iberian Basin, where it creates meandering canyons (sickles or snouts).

In this last sector of the upper Eresma, the river has constituted its valley along the Quaternary through successive stages of entanglement alternating with the formation of erosive terraces (shoulders of the ballast) or discontinuous deposits of gross detrital materials (rocky ground). There are only three levels of terraces in the valley bottom and only in the internal meanders faces [60].

Prior to the embedding of the river and coeval with it, a karstic system, essentially endokarstic, was developed in carbonate lithology (dolomites, limestones, marls and dolomitic sandstones). It was constituted by a wide network of galleries and conduits. Cavities are of very limited dimensions (hectometer length and metric height and width), showing a clear rocky-structural control, concentrated in certain carbonate levels previous to karstification or in which stabilization stages were concentrated in the nesting. Many of them, due to their small size, could be classified as balms, rock shelters and flaps, rather than caves, and their genesis is related to stabilization moments of the valley nesting and lateral erosion of the slopes by the fluvial-karstic system [53].

2.1.2. Site Description

The San Lázaro Rock Shelter was discovered in 2014 and three archaeological campaigns (2014, 2016 and 2018) have been carried out so far in the site. An excavation of 2 × 2 m was performed during the first campaign, unearthing an archaeological horizon attributed to the final Middle Paleolithic. In the following stages, only sampling works were carried out in order to get a better definition of the documented stratigraphical sequence (Figure 2), since the aim is to initiate systematic excavations from 2019 onwards.

These works were made in parallel to the excavations in the nearby site of the Abrigo del Molino, which presents an interesting archaeological sequence to study the last Neanderthals occupation in the interior of the peninsula [57,61]. The sequence of the Abrigo del Molino is contemporary to the San Lázaro Rock Shelter, although, in this last case, at least three additional archaeological levels were documented. Both sites have several Neanderthal occupation layers with faunal and lithics remains dated at the end of the Middle Paleolithic in Southern Europe [53,57]. To date, only one publication on the discovery of the site exists [53], further results still remain unpublished.

2.2. Integral Digitalization of the Rock Shelter by Means of Remote Sensing and Geophysical Approaches

2.2.1. Aerial Photogrammetry

For the elaboration of the cartography of the site and its surroundings, photogrammetric techniques were used. All of the photographs follow a protocol of parallel photographic capture (axes perpendicular to the ground and parallel to each other, creating coplanar planes with the ground in each station), with an overlap area or horizontal-longitudinal coverage of 70% and 30% between passes (Table 3). The implemented photographs were acquired from a drone, model DJI Mavic Pro (Table 4, Table 5), planning several flights. Results of this work are useful to locate and characterize the archaeological site and its environment in a global way. Flights were planned taking into account the drone autonomy (battery life), the location and the extension of the land, as well as the orography (Table 3). Coordinates were established by using a GNSS (Global Navigation Satellite System) instrument with accuracies of ±1 cm to points evenly distributed throughout the study area. These points are easily identifiable in the photographs, with the aim of resolving the external orientation and geo-referencing the work. The reference system was the European Terrestrial System 1989 (ETRS 89).

2.2.2. Terrestrial Laser Scanner

In order to obtain a detailed model of the environment close to the rock shelter, a terrestrial laser scanner was used. More specifically, the lightweight terrestrial laser scanner Faro Focus 3D was used (Figure 3). This device measures distances using the principle of phase shift within a range comprised between 0.60 and 120 m with a measurement rate of 976,000 points per second, an angular resolution of 0.009º and a beam divergence of 0.19 mrad. Concerning its accuracy, this laser allows for capturing scenes with ±2 mm of precision under normal lighting conditions.

2.2.3. Wearable Mobile Mapping System

Due to the geometrical complexity of the rock shelter, placed inside a narrow space where 90% of its area has a height lower than 0.80 m and around 50% under 0.50 m, the use of the TLS was not possible. In order to solve this limitation, a WMLS, Zeb-Revo [62], was employed. This device is commercialized by the GeoSLAM company (City, US State abbrev. if applicable, Country) and integrates a rotating 2D scanning device (Hokuyo UTM-30LX-F, City, Country) rigidly coupled to an inertial measurement unit mounted on a rotary engine. These sensors are complemented by a battery and a data storage system located in a small backpack (Figure 4, Table 6). Regarding the accuracy of the device, the manufacturer declared that its value is 1–3 cm in relative terms.

The motion of the scanning head, i.e., the rotary engine, within the movement of the operator and the information captured by an inertial measurement unit (IMU) allow for obtaining the 3D information of the environment. To this end, a 3D SLAM algorithm was used to combine the information coming from the 2D laser head and the data provided by the IMU. More specifically, the full SLAM approach was used, implemented in the robotic operative system (ROS) library [63]. In order to obtain optimal results, it was necessary to follow a close-loop path. In this sense, the operator needed to start and finish at the same position, following a circular loop rather than a “there and back” loop where the path simply doubles back on itself [17,18]. The use of this protocol allowed the compensation of error propagation obtained during the data acquisition.

2.3. Geophysical Methods

2.3.1. Electrical Resistivity Tomography (ETR)

The electrical resistivity tomography (ERT) method, in two dimensions, is a geo-electrical technique used in the subsoil characterization. The system measures the apparent resistivity (AR) from a given tetra electrode device. With that aim, it injects electricity with a known intensity into two electrodes called A and B and automatically registers the potential difference between the other two electrodes called M and N. The equipment automatically modifies the distances between the pairs of electrodes obtaining the apparent resistivity in multiple positions and levels (n). These data are later treated by means of mathematical inversion algorithms obtaining an image of resistivity and real depths of the subsoil. The interpretation of this image allows for identifying the different lithology and structure of the subsoil [45].

The contrast of resistivity obtained from the implementation of this technique allows for the differentiation of the subsoil materials according to their electrical behavior, that is, in terms of their resistivity value, since it mainly depends on the following factors:

1. Pore volume proportion versus total rock volume. Initially, the greater the porosity, the lower the resistivity if the porosity is filled (water, clay, etc.). In the case of empty pores (air filling), a resistivity increase will be produced due to the dielectric air nature.

2. Pores geometric disposition (formation factor). The greater the pore connection, the lower the resistivity, since the fluid and ion mobility is easier.

3. Proportion of pores filled with water versus empty pores. The greater the proportion of water-filled pores, the resistivity decreases as water allows the circulation of electric current, as opposed to air (dielectric).

4. Resistivity or conductivity of the water or fluid that is filling the pores. The higher the water conductivity, as in the case of seawater, the lower the resistivity of the formation containing it.

The selection and use of the ERT technique for the realization of the present works are based on the following factors:

1. It allows the obtaining of precise two-dimensional resistivity subsoil sections that enable the determination of the position of the limits and materials structure.

2. It is a non-invasive technique that allows the realization of work in any kind of grounds, environment or installation.

Equipment used in the field measurements was a Syscal Pro device from the Iris instruments company (Figure 5).

Given the conditions of study (i.e., type, objective, site depth and required resolution), an investigation was carried out using two high resolution ERT profiles (profile 1 Root Mean Square Error (RMSE) = 23.1% and profile 2 RMSE = 21.6%). Each of them was performed using 48 electrodes and an inter-electrode spacing of 1 m and direct and reverse Polo–Dipole recording method. With this configuration, a penetration higher than 15 m was achieved, establishing 46 registration levels and obtaining more than 2000 measurement points in each profile.

2.3.2. Ground Penetration Radar (GPR)

The GPR technique stands out for being based on the emission of short duration electromagnetic impulses, which are reflected and detected by the receiving antenna when an object or a discontinuity area is intercepted [43,44]. Of all the energy that arrives at the object, only a part of it is reflected (depending on the electrical properties of the object), continuing the remaining energy to the reflection in new objectives until its total cushioning.

If the conditions of the ground and the element to be detected are optimal (for example, if we have different materials on a dry surface and we want to find walls of different materials than those constituting the ground), GPR usually presents very good results.

On the contrary, if the ground humidity is high or there are several clayey levels, GPR will not make effective measurements if the electromagnetic waves cannot penetrate the ground. It should be noted that results obtained by the use of GPR (in the form of radar-gram) are fully interpretable by the operator. The profiles obtained in the field (depth–length) have shown an excellent image with the standard signal treatment. With further treatment time-depth transformation (depth slices), it has been possible to verify that the depth of the reflectors coincides with cave or karstic features Depending on the GPR’s MHz, it can deepen further, getting better or worse resolution (400 MHz—depth 4 m—resolution 10–20 cm), depending, among other factors, on the presence of clays or high humidity states. Therefore, GPR is an ideal instrument for locating elements at shallow depths at high resolution. With many GPR sections, extensive subsoil maps can be generated at different depths. In archaeology, it is the most used application since it works in plant layout.

The objective of using this technique in this work is to detect different reflectors linked to variations in lithology or stratifications, as well as gaps: caves and karst phenomena (chimneys, sinkholes, etc.)

The GPR system used was a commercial device of the IDS Company (City, US State abbrev. if applicable, Country) characterized by a central frequency antenna of 400 MHz. Acquisition was carried out with software IDS K2 and a window length of 120 nanoseconds with 512 samples was selected for the profiles shown in this work.

Results were processed by the GRESWIN 2D (Table 7) software applying to this end a bandpass and a noise removal filter, an improvement of the amplitude of the most notable reflectors and a time distance conversion (Figure 6). The usual research depth of this type of device is about 4 m with a resolution of approximately 10 to 20 cm (the greater the depth the signal dissipates and the wave amplitude increases, it is thus more difficult to detect hollow anomalies of small size). Thus, the signal “bounces” in the mentioned holes or reflectors and is collected after the so-called “double time”.

Different longitudinal and transverse profiles were made (and geo-located by GNSS) in the upper part of the rock shelter in order to know the cavity starting point, fracturing and possible presence of karstic chimneys.

In order to complement the data obtained by the ERT campaign, a total of 6 GPR profiles were carried out on the upper part of the rock shelter: (i) 3 profiles parallel to the rocky escarpment (P1 to P3) and (ii) 3 profiles orthogonal to the previous direction (P4 to P6) (Figure 6).

3. Results 3.1. Definition of the External Geometry by Means of the Combination SfM-SLAM

3.1.1. Digitalization at a Global Scale

Cartography based on photogrammetric techniques was generated, using as vertical aerial photography. Once the images were taken, they were processed with GRAPHOS open-source photogrammetric reconstruction software [64,65], obtaining the point cloud, the 3D model, the digital elevation model (DEM) and the orthophotography (average GSD of 3.7 cm, average scaling error of 3.1 cm, average photogrammetric error of 1.7 cm, average precision of 3.4 cm) (Figure 7).

3.1.2. Digitalization at Local Scale (Detailed Model)

Complementary to the digitalization carried out by means of the integration drone + SfM, an additional experimental campaign was carried out with the aim of obtaining a detailed model of the rock shelter. To this end, a LiDAR approach was used, combining the data coming from the TLS for the characterization of the environment close to the rock shelter and the data obtained by the WMLS for the analysis of those more complex and narrow spaces.

Previous to the collection of data, a network of registration spheres was placed along the scene with the aim of aligning the different point clouds generated in a common coordinate system (Figure 8). These registration spheres have a diameter of 200 mm and were placed along the scene, guaranteeing their visibility from different positions. It is worth mentioning that these registration spheres were set on topographic tripods equipped with an in-house platform that allows for georeferencing its centroid by means of a GPS device.

On the one hand, a total of three scan stations were required in order to digitalize the close environment of the rock shelter (Figure 8). The alignment of these scan stations was carried out by means of the target-based registration method. This method allows for aligning different point clouds through the use of geometrical features coming from artificial targets such as planar targets or registration spheres. For the present study case, a target-based registration approach able to use the centroid of each registration sphere was used as a control point for the alignment between point clouds. Within this framework, the centroid of each sphere was extracted by the RANSAC Shape Detector algorithm [66]. As a result, it was possible to align all the point clouds with an accuracy of ±3 mm. Complementary to this data acquisition, and, with the aim of evaluating the accuracy of the WMLS point cloud, an additional scan station and several registration spheres were placed inside the rock shelter (Figure 8).

On the other hand, the acquisition of the inner and more complex envelope of the rock shelter was performed through the use of the Zeb REVO WMLS. During this data acquisition, a close loop path was followed with the aim of compensating possible error accumulations (Figure 8). Data acquisition with the Zeb REVO device requires firstly the initialization of the IMU in order to establish the reference coordinate system. Then, the scanner head rotates, allowing for the capturing of 2D profiles. It is worth mentioning that the portability of the proposed WMLS allows for capturing the whole inner geometry of the rock shelter, spending a total of 2 mins for its acquisition.

In order to get the complete 3D model, a full SLAM approach was followed [67]. This strategy estimates the spatial position of the scan head through the use of the data captured by the IMU as well as the geometrical features extracted from each 2D profile. For more details about the workflow followed by the full SLAM approach, the reader is referred to [67]. The time spent to obtain the 3D model of the rock shelter by means of the WMLS was around 4 mins, 2 mins for the data acquisition and another 2 mins to solve the SLAM problem (Figure 9). In order to align the point cloud obtained by the WMLS with the TLS point clouds, another target-based registration phase was carried out, using to this end the spheres S2, S3 and S5 (Figure 8). The RMSE of this registration phase was 0.01 m.

To obtain a more in-depth evaluation of the potential of the WMLS solution for the mapping of narrow spaces, a point-to-point comparison was carried out. During this stage, the Multiscale Model to Model Cloud Comparison (M3C2) algorithm was used [68]. This approach allows for estimating the signed discrepancies between the WMLS and the TLS point clouds, and it is implemented in the open-source software CloudCompare v2.10 [69] In order to get reliable results, the Zeb-REVO point cloud was segmented, deleting non-common areas, i.e., those parts not visible from the scan station used to compare the results. The comparison of both point clouds provided a RMSE of 0.01 m (Figure 10), in line with the RMSE obtained during the alignment phase.

Finally, all of the point clouds (three coming from the TLS and one from the WMLS) were fused in a single point cloud. This point cloud was later decimated due to the high amount of overlapping points. During this stage, a distance filter with a threshold of 0.005 m was applied. As a result, it was possible to obtain a detailed model of the rock shelter and its inner environment made up by 23,402.999 points (77% of the original model).

As it was stated before, the registration spheres were equipped on topographic tripods with special platforms that allow for georeferencing the centroid of each sphere. Thanks to this, it was possible to geo-reference the point cloud obtained by means of a six-parameter Helmert transformation (three translations and three rotations) allowing for placing both models, the LiDAR and the SfM point clouds, in the same coordinate system.

3.2. Inner Characterization of the Rock Shelter through Geophysical Methods

With the purpose of obtaining the stratigraphy of the rock shelter as well as a possible presence of infill layers, ERT and GPR geophysical methods were used.

3.2.1. Electrical Resistivity Tomography

This section is based on the joint analysis of the two resistivity profiles carried out (Figure 11, Figure 12), it being possible to extract the following deductions:

1. Profiles show a similar structural style allowing defining the position of the materials of interest (conductive fillers). Profile-2 cuts profile-1 in an orthogonal way, clearly defining the morphology of the area occupied by filling conductive materials.

2. In these profiles, the presence of fillers above ground is clearly visible. In more detail, an alternation of sandstones-sands-sandstones is found. It can be differentiated by their resistivity contrast, being higher in the sandstones than in the sandy stratum.

3. Conductive materials consisting of fillings with possible archaeological interest can be perfectly characterized by their low resistivity value. They are presented only in the central part of profile-1 (Figure 11) and at the beginning of profile-2 (Figure 12). The thickness of the conductive fillers of interest is not constant and is around 2–3 m deep.

3.2.2. Ground Penetrating Radar

The amplitude of the signals captured by the GPR in each profile was later projected on a bidimensional graph, the so-called radargram. In this radargram, it was possible to observe the different echo obtained during the data acquisition. This echo represents anomalies inside the ground such as voids or even different layers. With the aim of improving the graphical representation of each GPR profile, several filters obtained were applied, namely: (i) background removal; (ii) vertical bandpass filter; (iii) linear gain; and (iv) smoothed gain.

It is worth mentioning that the GPR profiles were not closer to the cliff because this particular antenna has a signal cone that can reach the air, promoting the presence of false positives. In addition to this, the GPR profiles were mainly carried out with the aim of detecting a possible karstification. In all of them, the velocity conversion was checked with the hyperbola filling method.

As a result, it was possible to obtain information about the inner composition of the rock shelter, i.e., stratigraphic disposition, contacts and possible karstic areas (Figure 13).

1. P 1: this profile shows two clear reflectors that can correspond to the stratification level, the first of them is the contact between a more altered ground and a more compact rock mass underlying. These two reflectors (lines) are at a depth of 1 to 1.5 m. Abnormalities of amplitude in the signal that can represent small karstic cavities of 10 to 20 cm are detected. The ceiling of the cave could be placed in the last part of the profile at a depth of 1 m.

2. P 2: the profile shows a very marked line that probably indicates the transition between a more altered ground and healthy rock at about 0.80 m depth. Some karstic areas and buried metallic objects are detected (without necessarily having relevance). These features showed the typical multi-reflector signal with high amplitude rebounds. In the middle of the profile, the cave is crossed, special anomalies are not identified on it, perhaps the ceiling of the cave coincides with the stratification and the reflector is not marked.

3. P 3: a line between 1.0 and 0.8 m deep is recognized, indicating the transition between the ground and the rock.

4. P 4: as in all profiles, a reflector is detected indicating the horizon between the ground and the rock or upper more altered zone. A fracture or failure seems to be detected 5 m away from the beginning point.

5. P 5: In the middle of the profile, a low signal area is identified, surrounded by very markedly inclined reflectors that could indicate a karstic chimney.

6. P 6: A karstic area with an anomalous zone in the ceiling is detected. It could be a karst filling area.

It is worth mentioning that the maximum depth of the different GPR profiles extracted was about 2 m (Figure 13). This penetration depth is justified by the presence of humidity in the soil as well as the presence of clays.

4. Discussion

This study presents a combination of geophysical and remote sensing techniques for the applied study of archaeological sites, describing new developments that could be used to improve methodological approaches to archaeological fieldwork and documentation. These different techniques could provide a valuable means of contributing new data to Middle Palaeolithic and Neanderthal sites research, considering how the integration of these different methods can reveal morphological and geological data from a different perspective, even prior to any type of archaeological intervention in the area.

On one hand, 3D documentation provides a means of capturing the external geometry of a site, creating digital data in the form of 3D point clouds that have a high degree of accuracy. The possibility of combining data from structured light laser scanning and photogrammetric techniques with archaeological data, geo-referenced under the same reference system, can additionally reveal new information that might not have been possible from a purely archaeological perspective. Here, we have presented a study case that have been able to globally define the position of the site in 3D. From one perspective, the use of a novel WMLS have been able to facilitate the 3D documentation of complex areas, proving several advantages in comparison with other traditional remote sensing approach such as:

1. Static laser scanning: this method is able to digitalize narrow spaces on which it is not possible to place a static laser scanner, requiring also a small number of stations as stated di-Fillipo et al. [17]. Additionally, this method is able to general a 3D product, in form of a point cloud, around 10 times faster than the static laser scanner, being also in line with what the studied previously showed. Regarding the resolution, this method is able to provide a dense point cloud to perform dimensional analysis on a product whose expected accuracy is near 1 cm. This value could be considered enough for archaeological purposes.

2. Image-based modelling method: the geometrical characteristics of the shelter, with narrow spaces, hindering the use of standard SfM protocols. In this situation the only feasible solution was the use of the approaches proposed by Perfetti et al. [70,71], on which a 3D reconstruction is carried out by means of the combination of a SfM approach with a fish-eye lens. In this situation, the proposed approach would have a better performance due to the nature of the system (active sensor) not requiring artificial illumination. Regarding the data acquisition, the proposed approach only requires following a simple protocol, which could be considered simpler than the SfM one. With respect to the data processing, the SfM requires external measurements, by means of a topographic equipment or metric tapes, to scale the model. This issue is not necessary in the case of using the WMLS.

On the other hand, data provided by geophysical techniques, including GPR and ERT, have proven useful in characterizing the internal subfloor of the site, as well as a general characterization of the nature of the rock shelter’s roof (Figure 14). This GPR was used to highlight possible alterations and the internal stratigraphy of these different areas, while ERT was able to evaluate the sedimentary potential of the site alongside a number of additional geological characteristics; since the sediment of archaeological interest has a lower conductivity value than the filling materials, it was possible to spatially characterize the site.

Combined usage of these different techniques allow for a more developed understanding of the characteristics of a site (Figure 14), revealing data about the area’s formation, including geology, sedimentary and dimensional data, as well as the ability to evaluate the historical evolution of the area integrated with archaeological and environmental data. These contributions to the understanding of a site are useful for planning archaeological interventions, optimizing the available resources, and can save time and money as well as materials when preparing for the study and diffusion of a site.

5. Conclusions

The combination of geophysical and remote sensing techniques has been proven to be an excellent approach for the archaeological documentation as shown in the whole manuscript and discussion section. Multi- and inter-disciplinary studies are a fundamental component of any field of research. Integrating different disciplines is an efficient means of enhancing each individual’s potential, while providing a global and transversal vision of the object under study, while making the most of new techniques and methods to develop the data attainable. In archaeology, this is especially important and provides the most powerful means of interpreting a site, providing more and more data of higher quality. The final interpretation of the site refers to the precise identification of the grout of archaeological interest, allowing the detection of the site’s limits and optimizing the general resources required in the subsequent excavation works. Additional advantages of implementing team members from different disciplines is the difference in perspective that helps confront different problems and provides new means of solving them. This boosts our decision-making capacities for archaeological interventions and can aid in the overcoming of issues that arise on a daily basis. These multitude of advantages can be considered an important means of developing disciplines, such as archaeology, through providing faster, cost-worthy methods that could further our understanding of the archaeological register.

Moreover, WMLS has provided a greater resolution for the documented area of study, both from a global and local position, providing a means of analysing spatial relations from both an inter- and intra-site perspective. Alternatively, the data made available through studies of this type are also of great use for heritage management and diffusion, providing possibilities for the 3D digitalization of sites with the possibility of implicating these visual tools in web applications or museum digital displays. If implemented with virtual reconstructions of possible activities that could have been performed on the site, those interested could be provided with a more interactive, didactical and visual experience. Virtual, augmented and mixed reality applications could be an important component for future museums or educational applications.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

3D Digitalization Techniques	Image-Based Methods	Photogrammetry/Image Based Modelling	- Close Range Photogrammetry - Aerial Photogrammetry
	Range-based methods	Static laser scanner	- Terrestrial laser scanner
		Dynamic laser scanner	- Back-pack wearable mobile mapping system - Mobile LiDAR system

Archaeological surveys: phase characterised by the search of archaeological evidence and site location [12,13].	Orthophotography: (Present and historical) in order to find signs of anthropic activity (e.g., wall alignments, differential growth in crops that reveal geometric figures, etc.) Digital Elevation Model: used with the aim of detecting characteristic elements on the “z” axis (e.g., elevations, rocks alignments, that in view of a human eye are not perceptible in the terrain)
Site contextualization and documentation: phase of data collection of the site for its understanding [20,21,22,23,24,25].	Point cloud: geometric characterization of the heritage resource. Orthophotography: Location of characteristic elements of the studied resource.
Excavation processes: documentation of the elements appearing during the excavation process (it is an irreversible process, being impossible to recover the original state of the site) [26,27].	Points cloud: re-visualize what there was in each of the intervention processes of that resource, having georeferenced all the archaeological material extracted from it. Obtaining longitudinal and transversal profiles and surface and volumetric calculations. Orthophotography: obtaining a detailed planimetry of the excavation.
Preservation and analysis of the archaeological findings: processes in which different techniques are applied in order to check the deterioration suffered by the resource and preserve its conservation [28,29,30,31].	Points cloud: documentation of archaeological sites and objects with the aim of monitoring it and verifying its deterioration over time. Orthophotography: verification of the resource deterioration.
Taphonomical analyses of the archaeological findings: a process in which the materials unearthed in the site are studied to know the processes they suffered [32,33,34,35,36,37,38].	Points cloud: performing a morphometric and geometric analysis of the materials (or details of them) found in a given archaeological site (cutting marks, engravings in stone, coins, ceramics) for a better understanding of the activities and processes that affected the archaeological findings.
Dissemination: a process in which the resources (archaeological site, archaeological evidence) or what happened in them are put in value at the service of the public [1,3,4,5,6,7,8,22,25].	3D model: visualization of 3D models on web, applications in a graphic, visual and didactical way. Recreation of the site and the activities that were carried out in it. 3D models can also be used for the analysis of archaeological materials since they are reliable copies of the pieces. Orthophotos: diffusion of images from the site.Panoramic photographs: Generation of virtual visits.

Flight Planning
Ascent speed	0.5 m/s
Horizontal speed	3 m/s
Flight height	67 m
Rotation degree (respect to the north)	30° s
Horizontal-longitudinal overlap	70%
Overlap between strips	30%
Positioning system accuracy	50 cm
Surface	34 ha
Number of images	519

DJI Mavic Pro’s Camera
Autonomy	27 min
Total weight	734 gr
Maximum flight speed	65 km/h
Accuracy of the on-board GPS	±1.5 m (x,y) and ±0.5 m (z)
Maximum allowable height	5000 m over the sea

DJI Mavic Pro’s Camera
Sensor size	1/2.3″
Sensor resolution	12 MPix
Field Of View	78.8°
Focal length	5
Max. ISO	3200
Aperture	f/2.2 fixed
Shutter	Electronic

WMLS Zeb REVO
Measuring principle	Time of flight
Wavelength (nm)	905
Scanner resolution (°)	0.625 H × 1.8 V
Orientations system	MEMS IMU
Scanner dimensions (mm)	86 × 113 × 287
Total weight (kg)	4.10
Scanner weight (kg)	1.00
Dimensions (mm)	220 × 180 × 470
Working range (m)	0.60–30 m indoors
	0.60–15 m outdoors
Measurement rate	40,000 points per second

GRESWIN 2D Configuration
Xmin: −0.18
Xmax: 15.18
Move start time
Background removal
Vertical bandpass filter (td)
Linear gain
Smoothed gainTime-depth transformation: velocity = 110 mm/ns

Author Contributions

M.Á.M.-G.—Data collection and processing by using TLS and WMLS and manuscript writing; L.J.S.-A.—Data collection and processing by using TLS and WMLS and manuscript writing; C.S.B.—Data collection and manuscript writing; P.C.G.—Data collection, processing and interpretation by using ETR and manuscript writing; D.Á.-A.—Project direction, archaeological prospecting and manuscript writing; M.d.A.-H.—Project direction, archaeological prospecting and manuscript writing; J.C.G.-D.—Photogrammetric flight with drone and data processing; D.G.-A.—Data collection and processing by using TLS and WMLS.; M.H.R.—Photogrammetric flight with drone and data processing; L.J.B.—Data collection, processing and interpretation by using GPR; C.L.C.—Photogrammetric flight with drone and data processing; R.M. Data collection and processing by using TLS and WMLS.

Funding

This work has been partially funded by Junta de Castilla y León.

Acknowledgments

The authors would like to acknowledge the sponsorship and collaboration of the Junta de Castilla y León. We would also like to thank the TIDOP Group from the Department of Cartographic and Land Engineering of the Higher Polytechnics School of Avila, University of Salamanca, for the use of their tools and facilities. The authors also want to thank the Ministry of Education, Culture and Sport for providing an FPU Grant (Training of University Teachers Grant) to one of the author of this paper. Finally, the second author would like to thank the University of Salamanca for the program for human resources “Program II: Post-Doctoral Grant”.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Escarcena, J.C.; de Castro, E.M.; García, J.L.P.; Calvache, A.M.; del Castillo, T.F.; García, J.D.; Cámara, M.U.; Castillo, J.C. Integration of photogrammetric and terrestrial laser scanning techniques for heritage documentation. Virtual Archaeol. Rev. 2011, 2, 53–57.

2. El-Hakim, S.F.; Beraldin, J.A.; Picard, M.; Godin, G. Detailed 3D reconstruction of large-scale heritage sites with integrated techniques. IEEE Comput. Graph. Appl. 2004, 24, 21–29.

3. García, J.L.P.; Calvache, A.T.M.; Escarcena, F.J.C.; Arenas, A.L. Fotogrametría de bajo coste para la modelización de edificios históricos. Virtual Archaeol. Rev. 2011, 2, 121–125.

4. García-Morales, R.M.; Cortina, M.L.; Pintado, J.A.; Novella, L.R. Princeps Resurgens: Archaeological research and photogrammetric documentation in the study of a Roman thoracata statue of Los Bañales (Uncastillo, Zaragoza). Virtual Archaeol. Rev. 2015, 6, 65–71.

5. Sabina, J.Á.R.; Valle, D.G.; Ruiz, C.P.; García, J.M.M.; Laguna, A.G. Aerial Photogrammetry by drone in archaeological sites with large structures. Methodological approach and practical application in the medieval castles of Campo de Montiel. Virtual Archaeol. Rev. 2015, 6, 5–19.

6. Cano, P.; Lamolda, F.; Torres, J.C.; del Mar Villafranca, M. Uso de escáner láser 3D para el registro del estado previo a la intervención de la Fuente de los Leones de La Alhambra. Virtual Archaeol. Rev. 2010, 1, 89–94.

7. Cobo, E.P.; Vaz, P.M.; Labrado, E.O.; Zamora, Á.I.; Granero, A.B.M. Proyecto de musealización de los restos hallados en la estación de Ópera (Metro de Madrid). Reconstrucciones infográficas, escaneo laser 3D y digitalización del patrimonio arqueológico. Virtual Archaeol. Rev. 2012, 3, 88–92.

8. Fernández, F.L. Rescate documental de petroglifos y reconstrucción 3d del corredor Dolménico de Cubillejo de Lara, Burgos. Virtual Archaeol. Rev. 2016, 7, 43–52.

9. Koutsoudis, A.; Vidmar, B.; Ioannakis, G.; Arnaoutoglou, F.; Pavlidis, G.; Chamzas, C. Multi-image 3D reconstruction data evaluation. J. Cult. Herit. 2014, 15, 73–79.

10. Dhonju, H.; Xiao, W.; Sarhosis, V.; Mills, J.P.; Wilkinson, S.; Wang, Z.; Thapa, L.; Panday, U.S. Feasibility Study of Low-Cost Image-Based Heritage Documentation in Nepal. In International Archives of the Photogrammetry; Remote Sensing and Spatial Information Sciences: Göttingen, Germany, 2017; pp. 237–242.

11. Blais, F.; Rioux, M.; Beraldin, J.A. Practical considerations for a design of a high precision 3D laser scanner system. In Optomechanical and Electro-Optical Design of Industrial Systems; International Society for Optics and Photonics: Bellingham, DC, USA, 1988; pp. 225–246.

12. Cowley, D.C.; Moriarty, C.; Geddes, G.; Brown, G.L.; Wade, T.; Nichol, C.J. UAVs in Context: Archaeological Airborne Recording in a National Body of Survey and Record. Drones 2018, 2, 2.

13. O’Driscoll, J. Landscape applications of photogrammetry using unmanned aerial vehicles. J. Archaeol. Sci. Rep. 2018, 22, 32–44.

14. Campos, M.B.; Tommaselli, A.M.G.; Honkavaara, E.; Prol, F.D.S.; Kaartinen, H.; El Issaoui, A.; Hakala, T. A Backpack-Mounted Omnidirectional Camera with Off-the-Shelf Navigation Sensors for Mobile Terrestrial Mapping: Development and Forest Application. Sensors 2018, 18, 827.

15. Lagüela, S.; Dorado, I.; Gesto, M.; Arias, P.; González-Aguilera, D.; Lorenzo, H. Behavior Analysis of Novel Wearable Indoor Mapping System Based on 3D-SLAM. Sensors 2018, 18, 766.

16. Masiero, A.; Fissore, F.; Guarnieri, A.; Pirotti, F.; Visintini, D.; Vettore, A. Performance Eval-uation of Two Indoor Mapping Systems: Low-Cost UWB-Aided Photogrammetry and Backpack Laser Scanning. Appl. Sci. 2018, 8, 416.

17. Di Filippo, A.; Sánchez-Aparicio, L.; Barba, S.; Martín-Jiménez, J.; Mora, R.; González-Aguilera, D. Use of a Wearable Mobile Laser System in Seamless Indoor 3D Mapping of a Complex Historical Site. Remote Sens. 2018, 10, 1897.

18. Cabo, C.; Del Pozo, S.; Rodríguez-Gonzálvez, P.; Ordóñez, C.; González-Aguilera, D. Comparing Terrestrial Laser Scanning (TLS) and Wearable Laser Scanning (WLS) for Individual Tree Modeling at Plot Level. Remote Sens. 2018, 10, 540.

19. Nocerino, E.; Menna, F.; Remondino, F.; Toschi, I.; Rodríguez-Gonzálvez, P. Investigation of indoor and outdoor performance of two portable mobile mapping systems. In Videometrics, Range Imaging, and Applications XIV; International Society for Optics and Photonics: Bellingham, DC, USA, 2017; p. 1033201.

20. González-Aguilera, D.; Muñoz-Nieto, A.; Gómez-Lahoz, J.; Herrero-Pascual, J.; Gutierrez-Alonso, G. 3D digital surveying and modelling of cave geometry: Application to paleolithic rock art. Sensors 2009, 9, 1108–1127.

21. González-Aguilera, D.; Muñoz, A.L.; Lahoz, J.G.; Herrero, J.S.; Corchón, M.S.; García, E. Recording and modeling Paleolithic caves through laser scanning. In Proceedings of the 2009 International Conference on Advanced Geographic Information Systems & Web Services, Cancun, Mexico, 1–7 February 2009; pp. 19–26.

22. Torres-Martínez, J.A.; Seddaiu, M.; Rodríguez-Gonzálvez, P.; Hernández-López, D.; González-Aguilera, D. A multi-data source and multi-sensor approach for the 3D reconstruction and web visualization of a complex archaelogical site: The case study of “Tolmo De Minateda”. Remote Sens. 2016, 8, 550.

23. Carrero-Pazos, M.; Vilas-Estévez, B.; Vázquez-Martínez, A. Digital imaging techniques for recording and analysing prehistoric rock art panels in Galicia (NW Iberia). Digit. Appl. Archaeol. Cult. Herit. 2018, 8, 35–45.

24. Jalandoni, A.; Domingo, I.; Taçon, P.S. Testing the value of low-cost structure-from-motion (SfM) photogrammetry for metric and visual analysis of rock art. J. Archaeol. Sci. Rep. 2018, 17, 605–616.

25. Naranjo, J.M.; Parrilla, Á.; de Sanjosé, J.J. Geometric characterization and interactive 3D visualization of historical and cultural heritage in the province of Cáceres (Spain). Virtual Archaeol. Rev. 2018, 9, 1–11.

26. Organista, E.; Domínguez-Rodrigo, M.; Yravedra, J.; Uribelarrea, D.; Arriaza, M.C.; Ortega, M.C.; Mabulla, A.; Gidna, A.; Baquedano, E. Biotic and abiotic processes affecting the formation of BK Level 4c (Bed II, Olduvai Gorge) and their bearing on hominin behavior at the site. Palaeogeogr. Palaeoclimatol. Palaeoecol. 2017, 488, 59–75.

27. Siebke, I.; Campana, L.; Ramstein, M.; Furtwängler, A.; Hafner, A.; Lösch, S. The application of different 3D-scan-systems and photogrammetry at an excavation—A Neolithic dolmen from Switzerland. Digit. Appl. Archaeol. Cult. Herit. 2018, 10, e00078.

28. Del Pozo, S.; Sánchez-Aparicio, L.; Rodríguez-Gonzálvez, P.; Herrero-Pascual, J.; Muñoz-Nieto, A.; González-Aguilera, D.; Hernández-López, D. Multispectral Imaging: Fundamentals, Principles and Methods of Damage Assessment in Constructions. In Non-Destructive Techniques for the Evaluation of Structures and Infrastructure; CRC Press: Boca Raton, FL, USA, 2016; p. 139.

29. Sánchez-Aparicio, L.J.; Del Pozo, S.; Ramos, L.F.; Arce, A.; Fernandes, F.M. Heritage site preservation with combined radiometric and geometric analysis of TLS data. Autom. Constr. 2018, 85, 24–39.

30. Bautista-De Castro, Á.; Sánchez-Aparicio, L.J.; Ramos, L.F.; Sena-Cruz, J.; González-Aguilera, D. Integrating geomatic approaches, Operational Modal Analysis, advanced numerical and updating methods to evaluate the current safety conditions of the historical Bôco Bridge. Constr. Build. Mater. 2018, 158, 961–984.

31. Sánchez-Aparicio, L.J.; Bautista-De Castro, Á.; Conde, B.; Carrasco, P.; Ramos, L.F. Non-destructive means and methods for structural diagnosis of masonry arch bridges. Autom. Constr. 2019, 104, 360–382.

32. Yravedra, J.; Diez-Martin, F.; Egeland, C.P.; Maté-González, M.Á.; Palomeque-González, J.F.; Arriaza, M.C.; Aramendi, J.; García Vargas, E.; Estaca-Gómez, V.; Sanchez, P.; et al. FLK West (Lower Bed II, Olduvai Gorge, Tanzania): A new early Acheulean site with evidence for human exploitation of fauna. Boreas 2017, 46, 816–830.

33. Yravedra, J.; Maté-González, M.Á.; Palomeque-González, J.F.; Aramendi, J.; Estaca-Gómez, V.; San Juan Blazquez, M.; García Vargas, E.; Organista, E.; González-Aguilera, D.; Arriaza, M.C.; et al. A new approach to raw material use in the exploitation of animal carcasses at BK (Upper Bed II, Olduvai Gorge, Tanzania): A micro-photogrammetric and geometric morphometric analysis of fossil cut marks. Boreas 2017, 46, 860–873.

34. Hermon, S.; Polig, M.; Driessen, J.; Jans, G.; Bretschneider, J. An integrated 3D shape analysis and scientific visualization approach to the study of a Late Bronze Age unique stone object from Pyla-Kokkinokremos, Cyprus. Digit. Appl. Archaeol. Cult. Herit. 2018, 10, e00075.

35. Esteve, J.; Zhao, Y.L.; Maté-González, M.Á.; Gómez-Heras, M.; Peng, J. A new high-resolution 3-D quantitative method for analysing small morphological features: An example using a Cambrian trilobite. Sci. Rep. 2018, 8, 2868.

36. Aramendi, J.; Arriaza, M.C.; Yravedra, J.; Maté-González, M.Á.; Ortega, M.C.; Courtenay, L.A.; González-Aguilera, D.; Gidna, A.; Mabulla, A.; Baquedano, E.; et al. Who ate OH80 (Olduvai Gorge, Tanzania)? A geometric-morphometric analysis of surface bone modifications of a Paranthropus boisei skeleton. Quat. Int. 2019, 517, 118–130.

37. López-Cisneros, P.; Linares-Matás, G.; Yravedra, J.; Maté-González, M.Á.; Estaca-Gómez, V.; Mora, R.; Aramendi, J.; Rodríguez Asensio, J.A.; Barrera-Logares, J.M.; González-Aguilera, D. Applying new technologies to the taphonomic study of La Lluera (Asturias, Spain). Geometric morphometrics and the study of bone surface modifications (BSM). Quat. Int. 2019, 517, 107–117.

38. Yravedra, J.; Maté-González, M.Á.; Courtenay, L.A.; López-Cisneros, P.; Estaca-Gómez, V.; Aramendi, J.; de Andrés-Herrero, M.; Linares-Matas, G.; González-Aguilera, D.; Álvarez-Alonso, D. Approaching raw material functionality in the Upper Magdalenian of Coímbre cave (Asturias, Spain) through geometric morphometrics. Quat. Int. 2019, 517, 97–106.

39. García-Morales, R.M.; Cortina, M.L. Paleocatálogo 3D: Photogrammetry for the realization of a high quality, accessible and free 3D Virtual Catalog. Virtual Archaeol. Rev. 2015, 6, 35–40.

40. Pierdicca, R.; Malinverni, E.S.; Frontoni, E.; Colosi, F.; Orazi, R. 3D visualization tools to explore ancient architectures in South America. Virtual Archaeol. Rev. 2016, 7, 44–53.

41. El-Qady, G.; Metwaly, M.; Drahor, M.G. Geophysical Techniques Applied in Archaeology. In Archaeogeophysics; Springer: Cham, Switzerland, 2019; pp. 1–25.

42. Teixidó, M.T.; Peña, J.A. The Utility of Geophysical Models in Archaeology: Illustrative Case Studies. In Archaeogeophysics; Springer: Cham, Switzerland, 2019; pp. 169–181.

43. Conyers, L.B.; Goodman, D. Ground-Penetrating Radar: An Introduction for Archaeologist; AltaMira Press: Lanham, MD, USA, 1997; pp. 149–194.

44. Conyers, L.B. Analysis and interpretation of GPR datasets for integrated archaeological mapping. Near Surf. Geophys. 2015, 13, 645–651.

45. Martín Nieto, I.; Farfan Martín, A.; Sáez Blázquez, C.; González-Aguilera, D.; Carrasco, P.; Vasco, E.F.; García, J.C. Use of 3D electrical resistivity tomography to improve the design of low enthalpy geothermal systems. Geothermics 2019, 79, 1–13.

46. Teixidó, T.; Peña, J.A.; Moya, L.; Alemán, B.; Claros, J. Archaeologic Superficial Prospection and GPR Exploration in the Roman Site of Ciavieja (Elejido, Spain); Internal Document, Ref. AGA/2014-104; Andalusian Institute of Geophysics (Granada University): Granada, Spain, 2014.

47. Teixidó, T.; Peña, J.A.; Ruso, J.A.; Alcalá, S. Assessment of Experimental Device and the GPR Resolution Suitable to Characterize the Cork Oak Root System. Case of Study la Alcaidesa Natural Park, Cádiz (SPAIN); Internal Document, Ref. AGA-ID/2013-19; Andalusian Institute of Geophysics (Granada University): Granada, Spain, 2013.

48. Papadimitrios, K.S.; Bank, C.G.; Walker, S.J.; Chazan, M. Palaeotopography of a Palaeolithic landscape at Bestwood 1, South Africa, from ground-penetrating radar and magnetometry. S. Afr. J. Sci. 2019, 115, 1–7.

49. Papadopoulos, N.G.; Tsourlos, P.; Tsokas, G.N.; Sarris, A. Efficient ERT measuring and inversion strategies for 3D imaging of buried antiquities. Near Surf. Geophys. 2007, 5, 349–361.

50. Fernández, G.; Teixidó, T.; Peña, J.A.; Burillo, F.; Claros, J. Using shallow geophysical methods to characterise the monumental building at the Segeda I site (Spain). J. Archaeol. Sci. Rep. 2015, 2, 427–436.

51. Berge, M.A.; Drahor, M.G. Electrical resistivity tomography investigations of multi-layered archaeological settlements: Part I—Modelling. Archaeol. Prospect. 2011, 18, 159–171.

52. Martínez-Moreno, F.J.; Galindo-Zaldívar, J.; Pedrera, A.; Teixido, T.; Ruano, P.; Peña, J.A.; González-Castillo, L.; Ruiz-Constán, A.; López-Chicano, M.; Martín-Rosales, W. Integrated geophysical methods for studying the karst system of Gruta de las Maravillas (Aracena, Southwest Spain). J. Appl. Geophys. 2014, 107, 149–162.

53. Álvarez-Alonso, D.; de Andrés Herrero, M.; Herrero, A.D.; Hernández, J.A.R. Análisis geoarqueológico de las ocupaciones musterienses en el valle alto del río Eresma: El Abrigo del Molino (Segovia, España). Bol. Geol. Min. 2018, 129, 153–182.

54. Álvarez-Alonso, D.; de Andrés-Herrero, M.; Díez-Herrero, A.; Rojo Hernández, J. El abrigo del Molino (Segovia, España). Ocupaciones neandertales en el norte del Sistema Central. In VIII Reunión de Cuaternario Ibérico. Sevilla-La Rinconada; El Cuaternario Ibérico. Investigación en el siglo XXI; AEQUA: Avila, Spain, 2013; pp. 91–94.

55. Álvarez-Alonso, D.; de Andrés Herrero, M.; Herrero, A.D.; Hernández, J.A.R. El Abrigo del Molino (Segovia, España). Un nuevo yacimiento musteriense en el interior de la Península Ibérica. In Proceedings of the Investigaciones Arqueológicas en el valle del Duero, del Paleolítico a la Antigüedad Tardía: Actas de las III Hornadas de Jóvenes Investigadores del Valle del Duero, Salamanca, Spain, 20–22 November 2013.

56. Álvarez-Alonso, D.; de Andrés Herrero, M.; Herrero, A.D.; Rojo, J.; Medialdea, A.; Benito, G.; Vegas, J. Geoarqueología del yacimiento musteriense del Abrigo del Molino (Segovia). In Proceedings of the XIII Reunión Nacional de Geomorfología on Avances de la Geomorfología en España, Cáceres, Spain, 9–12 September 2014; pp. 446–449.

57. Álvarez-Alonso, D.; de Andrés-Herrero, M.; Díez-Herrero, A.; Medialdea, A.; Rojo-Hernández, J. Neanderthal settlement in central Iberia: Geo-archaeological research in the Abrigo del Molino site, MIS 3 (Segovia, Iberian Peninsula). Quat. Int. 2018, 474, 85–97.

58. Fernández García, P. Geomorfología del Sector Comprendido Entre el Sistema Central y el Macizo de Santa María La Real de Nieva (Segovia). Ph.D. Thesis, Universidad Complutense de Madrid, Madrid, Spain, 1987; p. 336.

59. Díez-Herrero, A.; Martín-Duque, J.F. Las raíces del paisaje. Condicionantes geológicos del territorio de Segovia. In Colección Hombre y Naturaleza, VII; Abella Mardones, J.A., Salinas, B., Yoldi, L., Eds.; Junta de Castilla y León: Salamanca, Spain, 2005; p. 464.

60. Fernández, P. Evolución cuaternaria y sistemas de terrazas en la subfosa terciaria de Valverde del Majano y el Macizo de Santa María la Real de Nieva (Segovia). Bol. Real Soc. Esp. Hist. Nat. (Secc. Geol.) 1988, 84, 69–83.

61. Kehl, M.; Álvarez-Alonso, D.; de Andrés-Herrero, M.; Díez-Herrero, A.; Klasen, N.; Rethemeyer, J.; Weniger, G.C. The rock shelter Abrigo del Molino (Segovia, Spain) and the timing of the late Middle Paleolithic in Central Iberia. Quat. Res. 2018, 90, 180–200.

62. Geo-SLAM Zeb-Revo. Available online: http://geoslam.com/hardware-products/zeb-revo/ (accessed on 17 January 2018).

63. Quigley, M.; Conley, K.; Gerkey, B.; Faust, J.; Foote, T.; Leibs, J.; Berger, E.; Wheeler, R.; Ng, A.Y. ROS: An open-source Robot Operating System. In Proceedings of the ICRA Workshop on Open Source Software, Kobe, Japan, 12-17 May 2009.

64. González-Aguilera, D.; López-Fernández, L.; Rodríguez-González, P.; Guerrero-Sevilla, D.; Hernández-López, D.; Menna, F.; Nocerino, E.; Toschi, I.; Remondino, F.; Ballabeni, A.; et al. Development of an All-Purpose Free Photogrammetric Tool; ISPRS International Archives: Prague, Czech Republic, 2016.

65. González-Aguilera, D.; López-Fernández, L.; Rodríguez-González, P.; Guerrero-Sevilla, D.; Hernández-López, D.; Menna, F.; Nocerino, E.; Toschi, I.; Remondino, F.; Ballabeni, A.; et al. InteGRAted Photogrammetric Suite, GRAPHOS. In Proceedings of the Congress: CATCON7-ISPRS, Prague, Czech Republic, 12–19 July 2016.

66. Schnabel, R.; Wahl, R.; Klein, R. Efficient RANSAC for point-cloud shape detection. Comput. Graph. Forum 2007, 26, 214–226.

67. Thrun, S. Simultaneous Localization and Mapping. In Robotics and Cognitive Approaches to Spatial Mapping; Springer: Berlin/Heidelberg, Germany, 2007; pp. 13–41.

68. Lague, D.; Brodu, N.; Leroux, J. Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (NZ). ISPRS J. Photogramm. Remote Sens. 2013, 82, 10–26.

69. Cloud Compare Project. Available online: https://www.danielgm.net/cc/ (accessed on 30 June 2019).

70. Perfetti, L.; Polari, C.; Fassi, F. Fisheye photogrammetry: Tests and methodologies for the survey of narrow spaces. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2017, 42, 573–580.

71. Perfetti, L.; Polari, C.; Fassi, F.; Troisi, S.; Baiocchi, V.; Del Pizzo, S.; Giannone, F.; Barazzetti, L.; Previtali, M.; Roncoroni, F. Fisheye Photogrammetry to Survey Narrow Spaces in Architecture and a Hypogea Environment. In Latest Developments in Reality-Based 3D Surveying and Modelling; MDPI: Basel, Switzerland, 2018; pp. 3–28.