Reconstructing the Roman Site “ Aquis Querquennis ” ( Bande , Spain ) from GPR , T-LiDAR and IRT Data Fusion

This work presents the three-dimensional (3D) reconstruction of one of the most important archaeological sites in Galicia: “Aquis Querquennis” (Bande, Spain) using in-situ non-invasive ground-penetrating radar (GPR) and Terrestrial Light Detection and Ranging (T-LiDAR) techniques, complemented with infrared thermography. T-LiDAR is used for the recording of the 3D surface of this particular case and provides high resolution 3D digital models. GPR data processing is performed through the novel software tool “toGPRi”, developed by the authors, which allows the creation of a 3D model of the sub-surface and the subsequent XY images or time-slices at different depths. All these products are georeferenced, in such a way that the GPR orthoimages can be combined with the orthoimages from the T-LiDAR for a complete interpretation of the site. In this way, the GPR technique allows for the detection of the structures of the barracks that are buried, and their distribution is completed with the structure measured by the T-LiDAR on the surface. In addition, the detection of buried elements made possible the identification and labelling of the structures of the surface and their uses. These structures are additionally inspected with infrared thermography (IRT) to determine their conservation condition and distinguish between original and subsequent constructions.


Introduction
Spain is very rich in cultural heritage, with castles, walls, churches, amphitheaters and monumental buildings of outstanding workmanship spread all over the country.These structures are vulnerable to changing weather patterns and other environmental hazards, thus requiring special protection against these events for the sake of conservation [1].Their protection and management include material characterization, structural stability analyses and archaeological prospecting.Regarding the latter, the full characterization of archaeological sites using conventional techniques is a long process, while coring and excavations, which are most frequently applied in archaeological assessment, are ground disturbing.The use of non-destructive testing (NDT) is therefore recommended, as these techniques can provide useful information about the conditions of conservation without intrusive intervention [2].Geophysics, aerial archaeology and other remote sensing techniques can be used to enhance the identification and understanding of the hidden archaeological remains [3].
The application of three-dimensional (3D) technologies to reconstruct archaeological sites has represented a great benefit within the heritage community [4].In this regard, the use of Terrestrial Light Detection and Ranging (T-LiDAR) for digital recording and documentation of cultural heritage items has increased significantly in recent years [5][6][7].This technique generates point clouds with 3D Cartesian coordinates and possibly color information that are useful for the visualization of the as-built realities.Furthermore, due to its capacity to provide fast, dense and accurate measurements, T-LiDAR presents several applications for structural monitoring of cultural heritage structures [8] and civil infrastructure facilities, such as tunnels [9,10], bridges [11], breakwaters [12] and roads [13].
However, buried structures cannot be assessed by visual or optical inspection and consequently, data integration and fusion techniques are needed for the correlation of sub-surface data with its corresponding superficial element.Alternative methods are therefore used for sub-surface inspection, including geophysics, such as magnetic, electrical and electromagnetic methods [14].Regarding the matter of archaeogeophysics, one of the most noteworthy methods has been ground-penetrating radar (GPR).GPR is a fast data acquisition technique that has been commonly applied for high-resolution imaging in many archaeological applications [15].However, the understanding of the GPR data has been a long-term challenge among the non-geophysical community.The analysis of the radar signal and the interpretation of the measured data (2D-GPR images) are not straightforward.The use of 3D-GPR data processing and visualization has advanced the acceptance of the use of GPR for archaeological prospecting [16].The use of 3D imaging techniques and processing software produces more realistic images of buried archaeological remains [17][18][19][20][21][22][23], which allows not only the discovery, but also the obtainment of 3D reconstructions of buried structures [24,25].Furthermore, all the data produced can be combined in a Geographic Information System (GIS) to achieve a more comprehensive archaeological interpretation [26][27][28].
Infrared thermography (IRT) is a technique for measuring characteristics related to the thermal state of the materials, allowing for the detection of pathologies mainly on their surface [29], but also in the sub-superficial areas if they present a high influence on the thermal state of the object [30,31].
The aim of this paper is to show the results of a morpho-geophysical survey at the archaeological site of "Aquis Querquennis" (Bande, Spain) by means of T-LiDAR and GPR, complemented with IRT.The first techniques are powerful tools for the identification and mapping of targeted remains, while IRT is a key technique for the diagnosis and determination of the conservation state of these remains.Thus, the three of them may result in very useful support for archaeologists when dealing with the collection, exploration and preservation of our cultural heritage.GPR data processing is performed using the in-house tool "toGPRi", which has been developed for GPR signal processing that includes the generation of both 2D and 3D imaging.With this objective, the tool implements different filters for signal amplification and noise removal as well as additional topographic correction.In particular, the product generated is a 3D model of the sub-surface and its subsequent XY images or time-slices at different depths, in addition to overlaid images connecting reflections at different levels that provide a better understanding of the occupied underground space.All of these products are georeferenced and can be subsequently imported into a GIS environment and merged with other sources of information for integral data treatment.Thus, the hybrid outcomes derived from the methods of GPR and T-LiDAR, and IRT and T-LiDAR, can be merged into a single image, allowing the reconstruction of the visible reality of the archaeological remains while integrating the unexcavated structures detected in the sub-surface.

The Roman Site "Aquis Querquennis"
"Aquis Querquennis" was a military camp from Roman times, built in Bande (Galicia, Spain), on the banks of the Limia river.Its occupation dates from the last quarter of the first century to the middle of the second and, according to the findings of the conquest of the site by the legio septimagemina detachment, it was a mixed unit of cavalry and infantry.Findings and historical proofs indicate the establishment of the legion at the site during the reign of Vespasian (69-79 AD).The encampment was probably built to monitor the construction of the roads that communicated Bracara Augusta (Braga, Portugal) and Asturica Augusta (Astorga, Spain).The first archaeological excavations were made in the 1920s and new studies were authorized in 1975, centered on the Northwestern area [32].
The settlement presents a classic layout with a rectangular shape and two main orthogonal paths (Figure 1), occupying 3 hectares.The border of the site is surrounded by a wall with softened corners, 3.20 m in height, with semi-cylindrical battlements, which are separated from the interior constructions by 11 m.There is a V-shaped moat outside the wall, 4 m deep and 3 m wide.The encampment had four main doors, corresponding to the ends of the two aforementioned paths, but only two have been excavated.Three barracks (strigia) have been identified, each one dedicated to housing two centuria with their respective commands.The barracks are composed of two rectangular wings faced around a courtyard (compluvium) with a cistern to collect rainwater.Each centuria was formed by ten contubernium, which is a space for eight soldiers, distributed as six contubernium in one wing and the guard quarter and four contubernium in the opposite wing, with the suite of rooms for the centurion at the entrance.Each contubernium was divided into a small front room (arma) for the equipment of the soldiers, and a rear room (papilio) where the soldiers slept.The gates and hollows of the barracks are oriented to the south.
Moreover, there were two rectangular granaries raised on stone pillars and delimited by thick walls with buttresses, so the roofs are believed to have been vaulted.A building of almost square shape and rooms arranged around an impluvium (a sunken area designed to carry away the rainwater) has been found and is considered to be designated to health care (valetudinarium).
The central building of the site was probably used as headquarters, due to its greater dimensions of 34.8 × 32.1 m, the presence of a vestibule flanked by ambulatories, a large central courtyard, a basilica with three doors and a sacral-administrative area or official temple.were made in the 1920s and new studies were authorized in 1975, centered on the Northwestern area [32].
The settlement presents a classic layout with a rectangular shape and two main orthogonal paths (Figure 1), occupying 3 hectares.The border of the site is surrounded by a wall with softened corners, 3.20 m in height, with semi-cylindrical battlements, which are separated from the interior constructions by 11 m.There is a V-shaped moat outside the wall, 4 m deep and 3 m wide.The encampment had four main doors, corresponding to the ends of the two aforementioned paths, but only two have been excavated.Three barracks (strigia) have been identified, each one dedicated to housing two centuria with their respective commands.The barracks are composed of two rectangular wings faced around a courtyard (compluvium) with a cistern to collect rainwater.Each centuria was formed by ten contubernium, which is a space for eight soldiers, distributed as six contubernium in one wing and the guard quarter and four contubernium in the opposite wing, with the suite of rooms for the centurion at the entrance.Each contubernium was divided into a small front room (arma) for the equipment of the soldiers, and a rear room (papilio) where the soldiers slept.The gates and hollows of the barracks are oriented to the south.
Moreover, there were two rectangular granaries raised on stone pillars and delimited by thick walls with buttresses, so the roofs are believed to have been vaulted.A building of almost square shape and rooms arranged around an impluvium (a sunken area designed to carry away the rainwater) has been found and is considered to be designated to health care (valetudinarium).
The central building of the site was probably used as headquarters, due to its greater dimensions of 34.8 × 32.1 m, the presence of a vestibule flanked by ambulatories, a large central courtyard, a basilica with three doors and a sacral-administrative area or official temple.

Scan Data Acquisition and Processing
A terrestrial laser scanner (T-LiDAR) generates 3D coordinates of an object point by measuring the distance between the scanner's center and the object and both the horizontal and vertical angles.

Scan Data Acquisition and Processing
A terrestrial laser scanner (T-LiDAR) generates 3D coordinates of an object point by measuring the distance between the scanner's center and the object and both the horizontal and vertical angles.Depending on the number of scans, the distance from the object to the scanner and the minimum angle increment of the system, a very dense point cloud can be achieved.Additionally, the reflectance of the surface can be measured by recording the intensity of the reflected laser beam, and for those scanners with an RGB camera attached, the RGB-value of the reflection is also captured.
A Faro Focus 3D X330 scanner (Valladolid, Spain), was used for this work.This is a phase-shift scanner with an effective range of 330 m at 90% reflectivity.Specifications indicate that the ranging noise of a single measurement is 0.3 mm at ranges up to 25 m (90% reflectivity) [33].The scanner has an integrated 8Mpix HDR camera to provide color to the point clouds and texture maps for the triangulated point cloud data.A topographic tripod was also used for its placement.
The location of the scan positions was planned beforehand in order to minimize the fieldwork and number of scans, while avoiding occlusions, to ensure the full coverage of the site.Actually, a total of 48 scan positions were performed with a base distance between 10 and 15 m, covering the whole monitoring area and obtaining their corresponding point clouds with a spatial resolution of 6-7 mm/10 m distance.
As the locations of the scan positions were not the same in different epochs, the coordinates of identical points sampled in different epochs were not expected to be equal either.Consequently, registration was required to combine multiple data into a single set of range data.The Faro Scene software provides an automatic registration method, using fixed plane targets in common surfaces as control points and is supported by the Iterative Closest Point (ICP) algorithm.Although registration may introduce small misalignment errors, these were in the order of a few millimeters with no significant effect on the final quality of the resulting point cloud [34,35].Finally, the raw laser point clouds were filtered and segmented in order to remove the points that were not part of the Roman site "Aquis Querquennis".

GPR Survey
Ground-penetrating radar (GPR) is a geophysical method based on the propagation of very short electromagnetic pulses (1-20 ns) in the frequency band of 10 MHz-2.5 GHz.A transmitting antenna emits the electromagnetic signal into the ground, which is partly reflected when it encounters media with different dielectric properties and is partly transmitted into deeper layers.Then, the reflections produced are recorded upon arrival at the receiving antenna; which is either in a separate antenna box, or in the same antenna box as the transmitter.The strength or intensity of the reflection is provided in terms of amplitude value.The amplitude is higher as the dielectric contrast between two different media is greater.While the antenna is moved along the ground surface, a two-dimensional image (radargram) is obtained, which is an XZ graphic representation of the detected reflections.The x axis represents the antenna's displacement along the survey line, and the z axis represents the two-way travel time of the pulse emitted (in terms of nanoseconds).If the time required for the electromagnetic pulse to go from the transmitting antenna to the reflector into the ground and return to the receiving antenna is measured, and the velocity of this pulse in the subsurface medium is known, then the position in depth of the reflector can be determined.

GPR Data Acquisition
The GPR survey was conducted using a RAMAC system from MALÅ Geosciences (Malå, Sweden).A 500-MHz antenna was used, because this frequency provides proper vertical resolution as well as sufficient penetration.The data acquisition was carried out using the common-offset mode with the antenna polarization perpendicular to the direction of data collection, and the acquisition parameters Remote Sens. 2018, 10, 379 5 of 16 were a 2 cm trace-distance interval and a total time window of 75 ns, composed of 512 samples.To complete the trace-distance interval calculation while measuring the profile lengths, the GPR antenna was mounted on a survey cart with an encoder (odometer wheel).
Three-dimensional (3D) GPR methodologies were performed in this work.As illustrated in Figure 1, two grids of approximately 156 m 2 and 787 m 2 were designed to prospect different unexcavated areas of the Roman site.To cover the entire grids, parallel 2D profiles were registered at regular intervals of 20 cm spacing, resulting in a total of 41 and 123 profiles, respectively.All the profile lines were collected in the same direction, starting from the lower left corner along the x axis.
3.2.2.GPR Data Processing: "toGPRi" Tool Data processing was performed using an in-house customized software tool called "toGPRi" [36], which was programmed using GNU Octave programming language [37].The developed software is compatible with the GNU General Public License v.3 and runs on GNU/Linux.It has been especially created to process data from MALÅ GPR systems in order to generate 3D cubes of data and subsequent georeferenced raster images at different levels of depth representing the underground elements.
As described in Figure 2, the required workflow consists of the following steps: (1) importing the data into the software; (2) checking the data configuration and geometry of the profiles (number of traces, frequency, number of samples, etc.); (3) displaying the 2D images (also called XZ images or radargrams) produced, including pre-processing of the data; (4) applying different filters for data processing, namely to remove noise and to amplify the signal as well as for topographic corrections; (5) exporting the processed 2D data as a series of ".xyz" profile files; (6) creating an Octave 3D matrix for the positioning of the profiles; (7) generating a 3D point cloud as well as XY images (or time-slices) at different levels of depth.
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 16 of 20 cm spacing, resulting in a total of 41 and 123 profiles, respectively.All the profile lines were collected in the same direction, starting from the lower left corner along the x axis.

GPR Data Processing: "toGPRi" Tool
Data processing was performed using an in-house customized software tool called "toGPRi" [36], which was programmed using GNU Octave programming language [37].The developed software is compatible with the GNU General Public License v.3 and runs on GNU/Linux.It has been especially created to process data from MALÅ GPR systems in order to generate 3D cubes of data and subsequent georeferenced raster images at different levels of depth representing the underground elements.
As described in Figure 2, the required workflow consists of the following steps: (1) importing the data into the software; (2) checking the data configuration and geometry of the profiles (number of traces, frequency, number of samples, etc.); (3) displaying the 2D images (also called XZ images or radargrams) produced, including pre-processing of the data; (4) applying different filters for data processing, namely to remove noise and to amplify the signal as well as for topographic corrections; (5) exporting the processed 2D data as a series of ".xyz" profile files; (6) creating an Octave 3D matrix for the positioning of the profiles; (7) generating a 3D point cloud as well as XY images (or time-slices) at different levels of depth.The "toGPRi" tool can open and read ".rd3" and ".rad" data files from MALÅ GPR systems, and the radargram is built and displayed using different color maps and contrast levels to optimize its visualization towards interpretation.
Regarding the data processing step, different filters were implemented to reduce clutter or any unwanted noise in the raw-data and to amplify the signal received.The purpose was to enhance the extraction of information from the signal received and to produce better data presentation, to make the data interpretation easier.The filters implemented were time-zero correction, temporal filters (subtract mean trace and vertical smooth), spatial filters (background removal and horizontal smooth), gain application (gain function and temporal constant gain or manual gain), and topographic corrections.Table 1 includes the sequence of filters and parameters considered for the data processing in the particular case study presented in this work.The "toGPRi" tool can open and read ".rd3" and ".rad" data files from MALÅ GPR systems, and the radargram is built and displayed using different color maps and contrast levels to optimize its visualization towards interpretation.
Regarding the data processing step, different filters were implemented to reduce clutter or any unwanted noise in the raw-data and to amplify the signal received.The purpose was to enhance the extraction of information from the signal received and to produce better data presentation, to make the data interpretation easier.The filters implemented were time-zero correction, temporal filters (subtract mean trace and vertical smooth), spatial filters (background removal and horizontal smooth), gain application (gain function and temporal constant gain or manual gain), and topographic corrections.Table 1 includes the sequence of filters and parameters considered for the data processing in the particular case study presented in this work.
Table 1.Filters and parameters used for data processing in the case study of the Roman Site "Aquis Querquennis".

Time-zero correction
Traces require adjustment to a common time-zero position, since thermal drift, electronic, instability, cable length differences or variations in the antenna air gap can cause "jumps" in the air/ground wavelength's first arrival time.This is usually achieved using some particular criteria such as the air wave first break point or the maximum amplitude peak of the trace.

Subtract mean trace
On some occasions, a continuous or very low frequency component (DC component) appears in the traces recorded by the radar, since the average level of the signal is moved from zero amplitude to a different value.This filter computes the corresponding sample mean value in each pixel of the trace to reduce this component.The filter is applied by selecting a percentage of traces on each side of the target trace to obtain the mean value.The subtract mean trace filter was applied in this work with a percentage of 50%.

Gain function
Signals at greater depths have very low energy due to signal attenuation and geometrical spreading.Gaining consists of amplifying the received signal by multiplying the data using a mathematical function or manually entering gain values.This filter applies a linear and exponential function to each sample.The function is defined as p * (t/s) * n + eq * (t/s) * n, where p is the linear factor, q is the exponential factor, t is the time window, s is the number of samples per trace and n is the number of samples filtered.
In this work, a gain function was applied from the first sample (once the time-zero was corrected) with a linear factor of 250 and an exponential factor of 10.

Temporal constant gain
This filter allows for the amplification of the signal from a certain depth by multiplying a selected number of samples by a certain factor.It affects all the traces composing the radargram at this temporal distance.This manual gain was applied in this work from sample 1 to sample 158 with a direct factor of 3.

Vertical smooth
Softens data vertically, in a specific horizontal section and time window, to remove high-frequency noise.This filter is applied by selecting a percentage of the mean of the number of samples.Vertical smooth was applied in this work from sample 1 to sample 150 and from trace 100 to trace 380 with a percentage of 35%.

Background removal
The main objective of these filters is to remove potential low-frequency noise which appears in the form of continuous horizontal bands along the recorded traces or only in some parts of them.This noise usually originates by bad coupling between the antenna and the medium and to eliminate ringing.This filter estimates an average of all traces in a window and removes it from every single trace.The main effect in the data is to suppress flat-lying reflectors, emphasizing smaller reflections.For every trace, this filter operates by calculating the mean value of the pixels (positive or negative), inferior to a pre-established percentage.The background removal filter was applied in this work with a percentage of 25%.

Horizontal smooth
Softens data horizontally, in a specific horizontal section and time window, to suppress horizontal reflections coming from the upper and lower interfaces or heterogeneous medium.This filter is applied by selecting a percentage of the mean of the number of traces.Horizontal smooth was applied in this work from sample 75 to sample 130 and from trace 1 to trace 949 with a percentage of 50%.

Topographic corrections
Compensating for topography is often important in improving the accuracy of imaging subsurface features.Features that are not directly underneath the antenna are recorded as if they actually were.The topographic data is loaded through a ".txt" file containing three columns with the XYZ coordinates for several points of the surface line, although only the XZ coordinates are used for correction.
The selection of the filters to be applied and their parameters are configured by the user.The results produced after each filter or after a sequence of different filters can be saved as an Octave binary matrix and the corresponding 8-bit image file (TIFF format).Figure 3 illustrates the results of the data processing described in Table 1.

correction.
The selection of the filters to be applied and their parameters are configured by the user.The results produced after each filter or after a sequence of different filters can be saved as an Octave binary matrix and the corresponding 8-bit image file (TIFF format).Figure 3 illustrates the results of the data processing described in Table 1.For the generation of the 3D matrix, a new pattern (.gpr) had to be established regarding the relative position of the profiles.The parameter of distance interval between traces is contained in the ".rad" file for each profile line.The distance between the parallel profile lines needs to be manually introduced.If the distance between profile lines is not regular along the entire grid prospected, the profile lines will be introduced in different phases.Other parameters to be introduced were the total time window, the original maximum number of samples (before the time-zero cut) and the average velocity of propagation of the GPR signals.Amplitude interpolation was used to fill the spaces between consecutive profile lines.The 3D cube obtained was finally saved as a variable named AY3D in a GNU Octave binary file with the extension "3D.b".Once this 3D matrix is built, the next step is to introduce the ".xyz" profile files in order to generate the 3D point cloud in-depth (Figure 4a).In addition, raster images can be generated through horizontal time-slices of the 3D cube and their georeferencing.Regarding the generation of the raster images, these time-slices are directly obtained with the "toGPRi" tool.These images can be obtained for different depth values depending on the length of the time window selected and the slice thickness.Moreover, the user-selectable option "overlaid raster" is also implemented for the computation of the maximum, minimum and mean pixel values of the fusion of layers (Figure 4b).This overlay analysis allows the strongest reflectors at different depths on a single image to be obtained.All the 3D point clouds generated can be exported to be managed in additional software like Cloud Compare (Paris, France) or MeshLab (ISTI-CNR, Pisa, Italia).Furthermore, the raster images produced can be saved in a TIFF image format after the configuration of an alpha channel and a range of transparent pixels (from 0 to 255, according to the 8-bit format) for the establishment of the threshold between black and white.However, if no alpha channel is selected, it will be an indexed color image.The georeferencing of images is also possible For the generation of the 3D matrix, a new pattern (.gpr) had to be established regarding the relative position of the profiles.The parameter of distance interval between traces is contained in the ".rad" file for each profile line.The distance between the parallel profile lines needs to be manually introduced.If the distance between profile lines is not regular along the entire grid prospected, the profile lines will be introduced in different phases.Other parameters to be introduced were the total time window, the original maximum number of samples (before the time-zero cut) and the average velocity of propagation of the GPR signals.Amplitude interpolation was used to fill the spaces between consecutive profile lines.The 3D cube obtained was finally saved as a variable named AY3D in a GNU Octave binary file with the extension "3D.b".Once this 3D matrix is built, the next step is to introduce the ".xyz" profile files in order to generate the 3D point cloud in-depth (Figure 4a).In addition, raster images can be generated through horizontal time-slices of the 3D cube and their georeferencing.Regarding the generation of the raster images, these time-slices are directly obtained with the "toGPRi" tool.These images can be obtained for different depth values depending on the length of the time window selected and the slice thickness.Moreover, the user-selectable option "overlaid raster" is also implemented for the computation of the maximum, minimum and mean pixel values of the fusion of layers (Figure 4b).This overlay analysis allows the strongest reflectors at different depths on a single image to be obtained.All the 3D point clouds generated can be exported to be managed in additional software like Cloud Compare (Paris, France) or MeshLab (ISTI-CNR, Pisa, Italia).Furthermore, the raster images produced can be saved in a TIFF image format after the configuration of an alpha channel and a range of transparent pixels (from 0 to 255, according to the 8-bit format) for the establishment of the threshold between black and white.However, if no alpha channel is selected, it will be an indexed color image.The georeferencing of images is also possible given that the "toGPRi" tool includes the geospatial GDAL library tools (libgdal and gdal-bin v2.1.1)[38], from OSGEO (USA), to generate GeoTIFF images and compressed KML files in a WGS84 coordinate system with UTM projection, whose zone can be selected by the user.As an example, Figure 4c shows the superposition of the 3D overlaid image in Figure 4b with the satellite image within Google Earth.
given that the "toGPRi" tool includes the geospatial GDAL library tools (libgdal and gdal-bin v2.1.1)[38], from OSGEO (USA), to generate GeoTIFF images and compressed KML files in a WGS84 coordinate system with UTM projection, whose zone can be selected by the user.As an example, Figure 4c shows the superposition of the 3D overlaid image in Figure 4b with the satellite image within Google Earth.

IRT Inspection
Infrared thermography (IRT) is a technique that measures radiation emitted by the bodies in the thermal infrared band of the spectrum.This radiation is a function of the thermal state of the bodies (temperature), and consequently, the surface temperature can be computed from it.Provided that the bodies tend towards thermal equilibrium if no heating is applied, anomalies in the temperature distribution of the bodies are associated with pathologies in their composition or state.Following this fact, a thermographic inspection is performed for the walls of the site, with the aim of providing information on the different stages of construction and determining which part of the remains are original and which belong to posterior reconstructions.This thermographic inspection is performed following the approach of passive thermography; that is, no artificial heating is applied, and the acquisition is performed without direct incidence of the sun on the surfaces.The emissivity of the materials is set to 1, so that the temperatures measured are apparent temperatures, and the differences in temperature in this case are due to differences on the superficial state.This way, temperature differences are associated to (1) different materials and (2) pathologies on the surface, such as cracks, mould or presence of vegetation (Figure 5).In addition, thermographic images are acquired with a 10% overlap between consecutive images and from a position of 10° with

IRT Inspection
Infrared thermography (IRT) is a technique that measures radiation emitted by the bodies in the thermal infrared band of the spectrum.This radiation is a function of the thermal state of the bodies (temperature), and consequently, the surface temperature can be computed from it.Provided that the bodies tend towards thermal equilibrium if no heating is applied, anomalies in the temperature distribution of the bodies are associated with pathologies in their composition or state.Following this fact, a thermographic inspection is performed for the walls of the site, with the aim of providing information on the different stages of construction and determining which part of the remains are original and which belong to posterior reconstructions.This thermographic inspection is performed following the approach of passive thermography; that is, no artificial heating is applied, and the acquisition is performed without direct incidence of the sun on the surfaces.The emissivity of the materials is set to 1, so that the temperatures measured are apparent temperatures, and the differences in temperature in this case are due to differences on the superficial state.This way, temperature differences are associated to (1) different materials and (2) pathologies on the surface, such as cracks, mould or presence of vegetation (Figure 5).In addition, thermographic images are acquired with a 10% overlap between consecutive images and from a position of 10 • with perpendicular to the walls.This angle avoids the measurement of radiation reflected from the operator, and the overlap minimizes missing information from any part of the walls.IRT also allows for the detection of sub-superficial anomalies if their effects are high enough to influence the state of the surface.In archaeology, buried elements have been found in [39], where the presence of air chambers under the surface produced temperature differences on the surface from which the presence of buried graves was identified with IRT, and in [40] where the soil coverage was freed from vegetation in order to detect small objects buried next to the surface.In the case of the "Aquis Querquennis" site, thermographic tests were performed from a zenithal point of view, in order to cover bigger extensions of the terrain from an angle of approximately 45°.The objective of the tests was to check whether buried walls could be located with IRT.

Results of the T-LiDAR Survey
The overall point cloud acquired by the Faro Focus 3D X330 system provides a model of the site (Figure 6a) that allows for detailed examination of individual barracks, including the qualitative and quantitative identification of geometric anomalies of structural elements.The walls and openings were all able to be captured, adding to the utility of the model for cultural heritage interpretation.Although there are a few areas inside the three barracks not accessible, this will not affect the overall documentation process.
Previous surveying of the site involved manual measurements, which required several days to measure and draw in CAD software packages.The resulting model was only accurate at the points where the measurements were taken, therefore leading to an approximate drawing.To overcome this setback, building dimensions were obtained from T-LiDAR orthoimages (Figure 6b).They were generated by applying an affine transformation to each one of the planes of the triangular mesh obtained from the point clouds.The same mathematical model directly relates the pixel coordinates of the registered RGB values from the 8Mpix HDR camera to those of the orthoimages.
The average accuracy of the resulting barrack lengths is typically in the order of one centimeter, which is within the tolerance expected given the minor misalignment errors and the accuracy of single measurements.In addition, thanks to the information about the portion of energy reflected by the surface of the barracks, which depends on its reflectance, it is possible to detect pathologies, measure the size and orientation of rock fissures or compute the damaged area.Lastly, it is noteworthy to mention that T-LiDAR derived measurements can be used as the starting points for advanced studies in the structural analysis of monuments and historical constructions.This fact IRT also allows for the detection of sub-superficial anomalies if their effects are high enough to influence the state of the surface.In archaeology, buried elements have been found in [39], where the presence of air chambers under the surface produced temperature differences on the surface from which the presence of buried graves was identified with IRT, and in [40] where the soil coverage was freed from vegetation in order to detect small objects buried next to the surface.In the case of the "Aquis Querquennis" site, thermographic tests were performed from a zenithal point of view, in order to cover bigger extensions of the terrain from an angle of approximately 45 • .The objective of the tests was to check whether buried walls could be located with IRT.

Results of the T-LiDAR Survey
The overall point cloud acquired by the Faro Focus 3D X330 system provides a model of the site (Figure 6a) that allows for detailed examination of individual barracks, including the qualitative and quantitative identification of geometric anomalies of structural elements.The walls and openings were all able to be captured, adding to the utility of the model for cultural heritage interpretation.Although there are a few areas inside the three barracks not accessible, this will not affect the overall documentation process.

Results of the GPR Prospection
Figure 7a represents the 3D overlaid GPR image obtained in grid 1 (see Figure 1), which corresponds to the unexcavated line of contubernia at the right wing of the barrack located in the Southern sector of the Roman camp.This overlaid image was built considering the time-slices (XY Previous surveying of the site involved manual measurements, which required several days to measure and draw in CAD software packages.The resulting model was only accurate at the points where the measurements were taken, therefore leading to an approximate drawing.To overcome this setback, building dimensions were obtained from T-LiDAR orthoimages (Figure 6b).They were generated by applying an affine transformation to each one of the planes of the triangular mesh obtained from the point clouds.The same mathematical model directly relates the pixel coordinates of the registered RGB values from the 8Mpix HDR camera to those of the orthoimages.
The average accuracy of the resulting barrack lengths is typically in the order of one centimeter, which is within the tolerance expected given the minor misalignment errors and the accuracy of single measurements.In addition, thanks to the information about the portion of energy reflected by the surface of the barracks, which depends on its reflectance, it is possible to detect pathologies, measure the size and orientation of rock fissures or compute the damaged area.Lastly, it is noteworthy to mention that T-LiDAR derived measurements can be used as the starting points for advanced studies in the structural analysis of monuments and historical constructions.This fact might help managers, archaeologists, and conservators to design and plan future interventions on this rapidly deteriorating Roman military camp.

Results of the GPR Prospection
Figure 7a represents the 3D overlaid GPR image obtained in grid 1 (see Figure 1), which corresponds to the unexcavated line of contubernia at the right wing of the barrack located in the Southern sector of the Roman camp.This overlaid image was built considering the time-slices (XY images) from 105 cm to 140 cm in depth (assuming a pre-calibrated velocity of propagation of 11.6 cm/ns).Given this image, it was possible to corroborate that, as usual, that sector of the barrack was composed of six contubernia of 6 × 3 m, divided into two zones (papilio and arma) and presenting a guard quarter (8 × 3 m size) at its entrance, its structure being similar to the barracks in the Northwestern sector. Figure 7b

Results of the GPR Prospection
Figure 7a represents the 3D overlaid GPR image obtained in grid 1 (see Figure 1), which corresponds to the unexcavated line of contubernia at the right wing of the barrack located in the Southern sector of the Roman camp.This overlaid image was built considering the time-slices (XY images) from 105 cm to 140 cm in depth (assuming a pre-calibrated velocity of propagation of 11.6 cm/ns).Given this image, it was possible to corroborate that, as usual, that sector of the barrack was composed of six contubernia of 6 × 3 m, divided into two zones (papilio and arma) and presenting a guard quarter (8 × 3 m size) at its entrance, its structure being similar to the barracks in the Northwestern sector. Figure 7b allows the validation of the GPR results after a recent excavation conducted in a 2016 campaign.Figure 8 presents the results produced for the GPR prospection acquired in grid 2 (see Figure 1).Although the wall reflections begin to appear from 50 cm in depth, the overlaid image in Figure 8 was built considering the time-slices from 70 cm to 150 cm.The purpose was to avoid unwanted reflections from wall falls in order to produce a clear image of the sub-surface structures, which makes the data interpretation easier.To enhance the wall reflections, the image was generated by selecting the maximum amplitude (or intensity) values and these were colored into a black/white palette with an alpha channel.Thus, the white pixels in the imaging represent the maximum amplitude values corresponding to wall reflections-either wanted reflections from sub-surface structures or undesired reflections from fallen walls.Observing the 3D image, it is also important to mention that some of the structure alignments highlighted are interrupted or missed due to the absence of walls by collapse.
The results exhibit the partial structure of a barrack, with two rectangular wings faced around a courtyard.The lower wing presents five contubernia (numbers 1-5 in Figure 8) and the guard quarter (number 9) at the entrance of the courtyard, while its counterpart has three contubernia (numbers 6-8) and the centurion's suite of rooms (number 10) at the entrance of the courtyard.All contubernia have similar dimensions and characteristics, as described for Figure 4. Unfortunately, it was not possible to prospect the complete structure, and the two opposite contubernia at the back of the courtyard and central cistern were missed in the GPR imaging.It is important to highlight that although the other excavated barracks present the centurion's suite at the left wing of the entrance (see Figure 1), the barrack under this grid has the centurion's suite at the right wing of the entrance.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 16 Figure 8 presents the results produced for the GPR prospection acquired in grid 2 (see Figure 1).Although the wall reflections begin to appear from 50 cm in depth, the overlaid image in Figure 8 was built considering the time-slices from 70 cm to 150 cm.The purpose was to avoid unwanted reflections from wall falls in order to produce a clear image of the sub-surface structures, which makes the data interpretation easier.To enhance the wall reflections, the image was generated by selecting the maximum amplitude (or intensity) values and these were colored into a black/white palette with an alpha channel.Thus, the white pixels in the imaging represent the maximum amplitude values corresponding to wall reflections-either wanted reflections from sub-surface structures or undesired reflections from fallen walls.Observing the 3D image, it is also important to mention that some of the structure alignments highlighted are interrupted or missed due to the absence of walls by collapse.
The results exhibit the partial structure of a barrack, with two rectangular wings faced around a courtyard.The lower wing presents five contubernia (numbers 1-5 in Figure 8) and the guard quarter (number 9) at the entrance of the courtyard, while its counterpart has three contubernia (numbers 6-8) and the centurion's suite of rooms (number 10) at the entrance of the courtyard.All contubernia have similar dimensions and characteristics, as described for Figure 4. Unfortunately, it was not possible to prospect the complete structure, and the two opposite contubernia at the back of the courtyard and central cistern were missed in the GPR imaging.It is important to highlight that although the other excavated barracks present the centurion's suite at the left wing of the entrance (see Figure 1), the barrack under this grid has the centurion's suite at the right wing of the entrance.

Results of IRT
The thermographic images acquired from the walls were registered with the T-LiDAR data through six control points in each image and their projective transformation.This way, the thermographic image presents temperature and geometric information, and the visualization of a complete wall in one image is possible in such a way that pathologies appearing in several images can be more clearly identified and interpreted.Figure 9 shows the thermographic orthoimages generated, where the following pathologies are highlighted: stones losing material, lack of union material and presence of vegetation.In addition, the presence of stones from two different times (before and after restoration activities) can be determined from their different temperatures, and because the temperature distribution in the newer stones is homogeneous.

Results of IRT
The thermographic images acquired from the walls were registered with the T-LiDAR data through six control points in each image and their projective transformation.This way, the thermographic image presents temperature and geometric information, and the visualization of a complete wall in one image is possible in such a way that pathologies appearing in several images can be more clearly identified and interpreted.Figure 9 shows the thermographic orthoimages generated, where the following pathologies are highlighted: stones losing material, lack of union material and presence of vegetation.In addition, the presence of stones from two different times (before and after restoration activities) can be determined from their different temperatures, and because the temperature distribution in the newer stones is homogeneous.For the test of detection of buried elements, although some images showed a high-temperature grid pattern on the terrain (Figure 10a), closer acquisitions revealed that high temperatures belong to areas with a lack of vegetation cover (Figure 10b).In addition, a partial excavation in the site showed that the current depth of the top of the walls is 21 cm (Figure 10c), while their width is 60 cm.In addition, the walls are made of sedimentary rock (sandstone), with very similar thermal properties (thermal resistance and diffusivity) as the burying sediments.This fact confirms the lack of capability of IRT to detect the buried parts of the wall, since it has been determined that the dimensions of the anomaly should be greater than its depth and that the materials should present different characteristics, in order to be detected with this technique [41].In the "Aquis Querquennis" site, although the width of the walls is greater than its depth, their geology and that of the burying sediments are the same (solid and disaggregated sandstone), not allowing for the existence of different heat fluxes and the consequent appearance of different temperatures.In order to confirm this incapability for detection, the inspection was performed after direct radiation of the sun on the area, trying to maximize the heat excitation of the material and forcing the appearance of a thermal flux in cases where it was possible.The non-appearance of a heat flux as a consequence of the geological equilibrium of the wall and burying materials, acts as a demonstration of the hypothesis.For the test of detection of buried elements, although some images showed a high-temperature grid pattern on the terrain (Figure 10a), closer acquisitions revealed that high temperatures belong to areas with a lack of vegetation cover (Figure 10b).In addition, a partial excavation in the site showed that the current depth of the top of the walls is 21 cm (Figure 10c), while their width is 60 cm.In addition, the walls are made of sedimentary rock (sandstone), with very similar thermal properties (thermal resistance and diffusivity) as the burying sediments.This fact confirms the lack of capability of IRT to detect the buried parts of the wall, since it has been determined that the dimensions of the anomaly should be greater than its depth and that the materials should present different characteristics, in order to be detected with this technique [41].In the "Aquis Querquennis" site, although the width of the walls is greater than its depth, their geology and that of the burying sediments are the same (solid and disaggregated sandstone), not allowing for the existence of different heat fluxes and the consequent appearance of different temperatures.In order to confirm this incapability for detection, the inspection was performed after direct radiation of the sun on the area, trying to maximize the heat excitation of the material and forcing the appearance of a thermal flux in cases where it was possible.The non-appearance of a heat flux as a consequence of the geological equilibrium of the wall and burying materials, acts as a demonstration of the hypothesis.For the test of detection of buried elements, although some images showed a high-temperature grid pattern on the terrain (Figure 10a), closer acquisitions revealed that high temperatures belong to areas with a lack of vegetation cover (Figure 10b).In addition, a partial excavation in the site showed that the current depth of the top of the walls is 21 cm (Figure 10c), while their width is 60 cm.In addition, the walls are made of sedimentary rock (sandstone), with very similar thermal properties (thermal resistance and diffusivity) as the burying sediments.This fact confirms the lack of capability of IRT to detect the buried parts of the wall, since it has been determined that the dimensions of the anomaly should be greater than its depth and that the materials should present different characteristics, in order to be detected with this technique [41].In the "Aquis Querquennis" site, although the width of the walls is greater than its depth, their geology and that of the burying sediments are the same (solid and disaggregated sandstone), not allowing for the existence of different heat fluxes and the consequent appearance of different temperatures.In order to confirm this incapability for detection, the inspection was performed after direct radiation of the sun on the area, trying to maximize the heat excitation of the material and forcing the appearance of a thermal flux in cases where it was possible.The non-appearance of a heat flux as a consequence of the geological equilibrium of the wall and burying materials, acts as a demonstration of the hypothesis.

Data Fusion and Visualization of "Aquis Querquennis" T-LiDAR and GPR Results
This final step involves the registration of T-LIDAR data with data from the GPR sensor in order to provide end-users with merged orthoimages of the archaeological site and its surrounding environment.In particular, the data integration process consists of finding corresponding control points in both datasets, comparable to the LiDAR registration process explained in Section 3.1.For this task, eight ground control points were located on the site in each of the four corners of both rectangular grids surveyed with the GPR technique (see Figure 1).
After that, both GPR grids were accurately registered in relation to the 3D geometric model defined by the TLS survey (see Figure 11).This multi-data source approach gives end-users the opportunity to generate orthoimages with different spatial resolutions.In addition, if the point cloud was already georeferenced, the Cartesian representation would enable the user to directly retrieve the coordinates of the points in an absolute coordinate system.In fact, those local coordinates could be transformed into global ones using, for example, a GNSS (Global Navigation Satellite System) base station.

Data Fusion and Visualization of "Aquis Querquennis" T-LiDAR and GPR Results
This final step involves the registration of T-LIDAR data with data from the GPR sensor in order to provide end-users with merged orthoimages of the archaeological site and its surrounding environment.In particular, the data integration process consists of finding corresponding control points in both datasets, comparable to the LiDAR registration process explained in Section 3.1.For this task, eight ground control points were located on the site in each of the four corners of both rectangular grids surveyed with the GPR technique (see Figure 1).
After that, both GPR grids were accurately registered in relation to the 3D geometric model defined by the TLS survey (see Figure 11).This multi-data source approach gives end-users the opportunity to generate orthoimages with different spatial resolutions.In addition, if the point cloud was already georeferenced, the Cartesian representation would enable the user to directly retrieve the coordinates of the points in an absolute coordinate system.In fact, those local coordinates could be transformed into global ones using, for example, a GNSS (Global Navigation Satellite System) base station.

Conclusions
The application of geomatic techniques at the Roman site "Aquis Querquennis", recognized as a cultural heritage resource, has widely contributed to its documentation and 3D reconstruction.The integration of terrestrial LiDAR and ground-penetrating radar was used for the first time in the encampment in order to support and complement the visual observations and petrophysical characterization carried out by previous restorers.The T-LiDAR method has shown its capabilities for measuring the 3D geometry of the overall settlement, letting users obtain, for example, the

Conclusions
The application of geomatic techniques at the Roman site "Aquis Querquennis", recognized as a cultural heritage resource, has widely contributed to its documentation and 3D reconstruction.The integration of terrestrial LiDAR and ground-penetrating radar was used for the first time in the encampment in order to support and complement the visual observations and petrophysical characterization carried out by previous restorers.The T-LiDAR method has shown its capabilities for measuring the 3D geometry of the overall settlement, letting users obtain, for example, the lengths, heights and volumes of every building, or the width of a particular stone pillar, without making direct contact with them.In particular, it provides highly accurate and more detailed 3D virtual reconstructions, incomparable with other large-volume measurement systems and their derived products, like the freely-available Google Earth or satellite imagery.Other interesting T-LiDAR derived products, such as 2D plans, sections and orthoimages, can also play important roles for heritage conservation and preservation purposes.
In addition, the use of GPR results revealed hidden constructions or structural details, such as the right wing of one barrack in grid 1 and the partial structure of another barrack in grid 2, which cannot be observed during visual exploration of the surface.With GPR data, the tool "toGPRi" presented in the paper allowed the generation of precise multi-source orthoimages, using the 3D-GPR imaging and the orthoimages from TLS, for a deeper analysis of both visible and underground space.At this point, the tool only supports 3D matrices for the representation of the buried 3D reality and 2D raster layers, but it is expected that for future work, vector layers will be implemented in order to compute parameters such as the volume of buried structures.In the aforementioned case study, the volumetric reconstruction would facilitate the dimensioning of the buildings while providing an intuitive and easily understood layout of the underground space distribution.In particular, the future goal of the "toGPRi" tool is to directly connect the processed results from GPR with LiDAR (both terrestrial and mobile) and DEM (Digital Elevation Model) surface data in geospatial databases with 2.5D and 3D geometries as SQLite/SpatiaLite/RasterLite2 and PostGIS.This way, the tool will be a single application based in a data infrastructure light model.
The inspection was complemented with IRT for the evaluation of the conservation state of the unburied elements.The registration of the thermographic images with T-LiDAR data allows for the direct location of the pathologies on their corresponding positions on the walls, as well as for their geometrical measurement in terms of area, number of stones and depth affected.The fact that buried elements and burying materials are the same (sandstone) has not allowed the detection of the first with IRT under the thermal excitation of the sun, regardless of the depth and dimensions of the buried element to detect.
Regarding novel integrations of techniques, future surveys at the site will deepen the application of infrared thermography with two objectives: first, the thermal characterization of the materials will be performed differentially in order to allow for quantitative diagnosis of the severity of the pathologies; second, information provided by GPR about the elements buried near the surface and their presence and location will be corroborated, through the application of specific thermal excitations and with knowledge of the thermal properties of the buried materials.

Figure 1 .
Figure 1.Situation of the Roman site "Aquis Querquennis" in Galicia (Northwest of Spain) and grids surveyed with the ground-penetrating radar (GPR) technique.

Figure 1 .
Figure 1.Situation of the Roman site "Aquis Querquennis" in Galicia (Northwest of Spain) and grids surveyed with the ground-penetrating radar (GPR) technique.

Figure 2 .
Figure 2. Workflow of the "toGPRi" tool developed for GPR data processing.

Figure 2 .
Figure 2. Workflow of the "toGPRi" tool developed for GPR data processing.

Figure 4 .
Figure 4. Different GPR results obtained with the "toGPRi" tool: (a) visualization of a horizontal slice of the 3D GPR cube in a point cloud viewer; (b) raster image generated through the overlay of several time-slices and (c) superposition of the 3D overlaid image with the aerial photograph of the site in Google Earth viewer.

Figure 4 .
Figure 4. Different GPR results obtained with the "toGPRi" tool: (a) visualization of a horizontal slice of the 3D GPR cube in a point cloud viewer; (b) raster image generated through the overlay of several time-slices and (c) superposition of the 3D overlaid image with the aerial photograph of the site in Google Earth viewer.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 16 perpendicular to the walls.This angle avoids the measurement of radiation reflected from the operator, and the overlap minimizes missing information from any part of the walls.

Figure 5 .
Figure 5. Examples of thermographic images acquired during the inspection of the excavated walls.(Left) completely original wall.(Right) partially reconstructed wall.

Figure 5 .
Figure 5. Examples of thermographic images acquired during the inspection of the excavated walls.(Left) completely original wall.(Right) partially reconstructed wall.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 16 might help managers, archaeologists, and conservators to design and plan future interventions on this rapidly deteriorating Roman military camp.

Figure 6 .
Figure 6.Results of the T-LiDAR survey: (a) Detailed view of a colored 3D point cloud of a stone wall and (b) orthoimage of the barrack in the Southern area of the Roman site

Figure 6 .
Figure 6.Results of the T-LiDAR survey: (a) Detailed view of a colored 3D point cloud of a stone wall and (b) orthoimage of the barrack in the Southern area of the Roman site.
Figure7arepresents the 3D overlaid GPR image obtained in grid 1 (see Figure1), which corresponds to the unexcavated line of contubernia at the right wing of the barrack located in the Southern sector of the Roman camp.This overlaid image was built considering the time-slices (XY images) from 105 cm to 140 cm in depth (assuming a pre-calibrated velocity of propagation of 11.6 cm/ns).Given this image, it was possible to corroborate that, as usual, that sector of the barrack was composed of six contubernia of 6 × 3 m, divided into two zones (papilio and arma) and presenting a guard quarter (8 × 3 m size) at its entrance, its structure being similar to the barracks in the Northwestern sector.Figure7ballows the validation of the GPR results after a recent excavation conducted in a 2016 campaign.

Figure 6 .
Figure 6.Results of the T-LiDAR survey: (a) Detailed view of a colored 3D point cloud of a stone wall and (b) orthoimage of the barrack in the Southern area of the Roman site

Figure 7 .
Figure 7. Results of grid 1 in the Southern sector of the camp: (a) 3D overlaid image and interpretation, and (b) aerial picture of a recent excavation.Figure 7. Results of grid 1 in the Southern sector of the camp: (a) 3D overlaid image and interpretation, and (b) aerial picture of a recent excavation.

Figure 7 .
Figure 7. Results of grid 1 in the Southern sector of the camp: (a) 3D overlaid image and interpretation, and (b) aerial picture of a recent excavation.Figure 7. Results of grid 1 in the Southern sector of the camp: (a) 3D overlaid image and interpretation, and (b) aerial picture of a recent excavation.

Figure 8 .
Figure 8. Results of grid 2 in the Northern sector of the camp: 3D overlaid image showing the configuration and interpretation of the contubernia identified (1-8), as well as the guard quarter (9) and the centurion's suite(10).

Figure 8 .
Figure 8. Results of grid 2 in the Northern sector of the camp: 3D overlaid image showing the configuration and interpretation of the contubernia identified (1-8), as well as the guard quarter (9) and the centurion's suite(10).

Figure 9 .
Figure 9. Thermographic image of the walls, registered on the Terrestrial Light Detection and Ranging (T-LiDAR) data for the visualization of complete walls within one product.The following pathologies are highlighted: stones losing material (red circle), lack of union material (red square) and presence of vegetation (green circle).Stones from before and after restoration activities are also marked (orange and purple curves, respectively).

Figure 10 .
Figure 10.(a) Thermographic image of the terrain, showing a possible grid pattern of high temperatures; (b) thermographic image of a detail of the terrain, where high temperatures are associated with a lack

Figure 9 .
Figure 9. Thermographic image of the walls, registered on the Terrestrial Light Detection and Ranging (T-LiDAR) data for the visualization of complete walls within one product.The following pathologies are highlighted: stones losing material (red circle), lack of union material (red square) and presence of vegetation (green circle).Stones from before and after restoration activities are also marked (orange and purple curves, respectively).

16 Figure 9 .
Figure 9. Thermographic image of the walls, registered on the Terrestrial Light Detection and Ranging (T-LiDAR) data for the visualization of complete walls within one product.The following pathologies are highlighted: stones losing material (red circle), lack of union material (red square) and presence of vegetation (green circle).Stones from before and after restoration activities are also marked (orange and purple curves, respectively).

Figure 10 .Figure 10 .
Figure 10.(a) Thermographic image of the terrain, showing a possible grid pattern of high temperatures; (b) thermographic image of a detail of the terrain, where high temperatures are associated with a lack Figure 10.(a) Thermographic image of the terrain, showing a possible grid pattern of high temperatures; (b) thermographic image of a detail of the terrain, where high temperatures are associated with a lack of vegetation; (c) visual demonstration of the depth of the top of the walls, from a wall partially excavated.

Figure 11 .
Figure 11.GPR results produced in grids 1 and 2, fused on the T-LiDAR based orthoimage of the Roman site.

Figure 11 .
Figure 11.GPR results produced in grids 1 and 2, fused on the T-LiDAR based orthoimage of the Roman site.