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”

Abstract

The complexity of archaeological sites hinders creation of an integral model using the current Geomatic techniques (i.e., aerial, close-range photogrammetry and terrestrial laser scanner) individually. A multi-sensor approach is therefore proposed as the optimal solution to provide a 3D reconstruction and visualization of these complex sites. Sensor registration represents a riveting milestone when automation is required and when aerial and terrestrial datasets must be integrated. To this end, several problems must be solved: coordinate system definition, geo-referencing, co-registration of point clouds, geometric and radiometric homogeneity, etc. The proposed multi-data source and multi-sensor approach is applied to the study case of the “Tolmo de Minateda” archaeological site. A total extension of 9 ha is reconstructed, with an adapted level of detail, by an ultralight aerial platform (paratrike), an unmanned aerial vehicle, a terrestrial laser scanner and terrestrial photogrammetry. Finally, a mobile device (e.g., tablet or smartphone) has been used to integrate, optimize and visualize all this information, providing added value to archaeologists and heritage managers who want to use an efficient tool for their works at the site, and even for non-expert users who just want to know more about the archaeological settlement.

1. Introduction

Several techniques have been applied thus far for the 3D reconstruction and visualization of archaeological settlements based on the use of close-range photogrammetry [1,2], terrestrial laser scanner (TLS) [3,4] or unmanned aerial vehicles (UAV) [5]. However, due to the inherent complexity of these sites, several problems arise when 3D reconstruction of these sites is mandatory and just one type of geotechnology is applied. For instance, aerial photogrammetry exhibits problems reconstructing vertical planes, common in archaeological sites, whereas terrestrial laser scanners or terrestrial photogrammetry could provide good results. However, these terrestrial techniques are subject to problems with horizontal surfaces or elevated areas. Other authors have explored multi-data and multi-sensor approaches to record and reconstruct complex archaeological sites. Recently in [6], the potential of this type of hybrid approach is shown for the analysis and interpretation of 3D/4D information applied to archaeological settlements. In order to guarantee geometric and radiometric quality, a combination of TLS and terrestrial photogrammetry is used and applied, enabling the monitoring of the settlement based on a volume analysis. However, this approach has problems related to the recording of more elevated areas. In order to solve this limitation, in [7], the authors use an aerial multi-sensor approach for the 3D reconstruction and visualization of archaeological settlements which provides a very good coverage of those elevated areas. The quality and precision of the TLS and UAV registration has been outlined in [8] where authors reconstruct the interior and exterior of the Church of Santa Barbara (Italy) after the earthquake it suffered in 2012. Trying to overcome the main UAV limitations, payload and autonomy, other authors have proved that low-cost manned platforms such as the paratrike can be an efficient solution for the recording of large sites [9], including archaeological sites [10], and allowing to put on board multiple sensors such as thermographic or multispectral cameras [11]. In those cases where the archeological site is complex and subterranean, other hybrid, dynamic and terrestrial approaches could be interesting. Unfortunately, there are not many mobile laser scanners, photogrammetric or hybrid systems for subterranean sites available in the market [12,13], even less specific for the field of archaeological recording. Some authors such as Canter et al. [14] have developed indoor mapping systems for the generation of indoor cartography from accurate geospatial information. The main advantage of these systems is that they integrate high-precision GNSS with advanced inertial technology (accelerometers and gyroscopes) for the geo-referencing of the site using measures from its exterior, apart from an Applanix POS system for positioning and orientation that provides uninterrupted measurements of the true position, roll, pitch and yaw of the vehicle moving indoors. Other authors have advanced more sophisticated systems based on a wheeled mobile robot and a multi-sensor global registration approach [15]. In particular, a geometric model to derive depth information is proposed based on a registration of heterogeneous 3D data arising from eight ultrasonic sonars, one TLS and three visual sensors.

It thus becomes clear that with the advances in multi-sensor and multi-data from different sources, data integration has become as a valuable tool in archaeological applications. The main objective of multi-sensor data integration is to register sensor data from different sources—with different characteristics, resolution, and quality in order to provide more reliable, accurate, and useful information required for diverse archaeological applications. In addition, this multi-data and multi-sensor integration can be improved through 3D geographical information systems (3DGIS) together with the final presentation of the products based on 3D Web.

The goal of this study is to propose a multi-data source and multi-sensor workflow in order to obtain high quality archaeological products, contributing to the robust interpretation of the observed objects/scenes and providing the basis of effective planning and decision making, essential in archaeology. To this end, terrestrial scans and images will be registered with aerial images acquired from a paratrike and an unmanned aerial vehicle (UAV). The methodology employed for data processing has been extensively tested by several authors. In [16], the authors employed computer vision algorithms for image-based modelling from aerial imagery and point cloud generation in urban environments, which are similar to those used in this work.

To carry out studies of characterization, measurement and analysis of the surface and the elements presented in the site, the geometrical information (acquired with geomatic sensors) is combined with other available thematic data such as photographs, sketches, restoration reports, schedules, etc. in a 3DGIS, Geoweb3D^® [17]. In this way, a description of the site through time is possible. Also, the availability of the three-dimensional information of the settlement through the Web using mobile devices provides added value for archaeologists and heritage managers, simplifying the data acquisition, as well as its analysis in the field. Moreover, the centralization of the information and its external storage means it becomes available to different experts and organizations. Last but not least, the knowledge of the settlement is open to the general public based on an easy-to-use interface which integrates different 2D and 3D resources.

In this study, a specific simplification and optimization procedure was implemented to visualize the different 3D products through the Web using mobile devices and the Open Source library Cesium [18], developed by the company Analytical Graphics, Inc. (AGI, Greenbelt, MD, USA) [19]. Other authors [20] have developed similar works using their own system, SGIS3D, which uses VRML (Virtual Reality Modeling Language) format. However, one of the main limitations of this language is the mandatory use of plugins, as well as its lack of optimization through the graphical processing unit (GPU), crucial in those steps related to texture mapping. Another similar approach in the archaeological field is developed by [21], who implements a spatial data infrastructure known as QuaeryArch3D and applied it to the archaeological settlement of “Maya de Copan” in Honduras, a UNESCO World Heritage Site. In this case, they use the Open Source PostgreSQL and PostGIS for the integration and management of information in the database, whereas the 3D visualization is performed with Unity. However, one of the main limitations of Unity is dealing with huge 3D models coming from the point clouds. To overcome this problem, WebGL seems to provide an appropriate solution using similar developments to those presented in [22] or developing a new approach as the one presented in this paper. Furthermore, advancing a 3DGIS Web solution requires using different simplification and optimization processes, as well as different hierarchical (i.e., pyramidal) strategies for visualization which are not compatible with the Unity engine.

This paper is organized as follows: after this introduction, in Section 2, the different sensors employed and their characteristics are detailed. In Section 3, the proposed multi-data source and multi-sensor approach is described, as well as the simplification and optimization process. Experimental results are shown and discussed in Section 4. Conclusions and future directions are given in Section 5. Finally, two appendices are included; the first corresponding to the abbreviations used in the paper, whereas the second encompasses the explanation of the methodology employed for 3D Web visualization.

2. Materials

2.1. Paratrike

The main aerial platform employed to the documentation of the archaeological site was a paratrike (Figure 1a). It is a low-cost aircraft with more flexibility than conventional aircrafts, and more autonomy and payload capacity than the UAVs. This last characteristic allows the possibility of boarding better sensors than the UAVs, or even multiple sensors in a stabilized gimbal (MUSAS-MUltiSpectral Airborne Sensors). In particular, a tandem trike AIRGES (

Table 1. Technical specifications of the paratrike.

Motor	Rotax 503 Two-Stroke Motor
Trike	Tandem Trike AIRGES
Tandem paraglide	MAC PARA Pasha 4 Trike 39 ó 42
Emergency system	Ballistic parachutes GRS 350
Weight	110 kg
Weight capability	165–250 Kg
Air velocity range	30–60 km/h

) was used to map the whole archaeological site following a vertical flight and using a full-frame reflex camera.

The use of the gyro-stabilized camera platform, MUSAS (Figure 1b), guarantees the accurate orientation of the camera according to the flight planning by two servomotors arranged on the x and y axes, controlled by an Arduino board, which incorporates an IMU with 6 degrees of freedom: three accelerometers, a double-shaft gyroscope (for pitch and roll) and an additional gyroscope for yaw.

For the paratrike, a full-frame reflex camera, Canon 5D MkII, with a fixed focal length to achieve a GSD of 3 cm, was used. This camera was also used for the terrestrial photogrammetry and photorealistic texture mapping due to its better image quality.

In contrast to the UAV platform which encloses an integrated navigation system (GPS, IMU and barometric altimeter), the paratrike requires an external set of sensors in order to provide navigation capabilities and thus fulfill photogrammetric constraints for data acquisition. In particular, the planimetric position is provided by a GPS antenna (Trimble Bullet III), installed in the camera platform close to the optical centre of the camera, connected to a mono-frequency receiver Ublox EVK-6T-0. This system yields an absolute precision of ±9 m on the horizontal axis for 95% of the time [23]. During data acquisition the pilot follows the planned photogrammetric mission in a rugged table connected to the GPS system, where the real-time track is contrasted with the planned flight. The final component in the navigation system, the altimetry, which affects the GSD, is controlled by an altimetric barometer (DigiFly VL100) with an absolute precision of ±8 m. This solution is chosen instead of GPS receiver, since the absolute altimetric GPS precision is just ±15 m.

2.2. Unmanned Aerial Vehicle

For the aerial photogrammetric acquisition in the “El Reguerón” site, the paratrike platform was rejected due to the morphological characteristics of the terrain and the height of flight required to reach the high spatial resolution (GSD of 1 cm). The main drawback was the need to vary the flying height to keep the same scale, due to the presence of high reliefs along the strips. As an additional disadvantage, the walled constructions were occluded between walls of natural rock. Therefore, in order to complete the archaeological site documentation, a UAV was employed. Specifically, a Microdrone md4-200 (

Table 2. Technical specifications of the Microdrone md4-200 platform.

UAV Weight	900 g
Payload	up to 200 g
Size	54 cm between rotors
Flight time	10 to 20 min
Operating temperature	−10 to 50 °C
Max. height flight	500 m
Max. wind	5 m/s

) was used (Figure 1c) to map the most challenging area following vertical and oblique flights through use of a compact camera.

In spite of the manoeuvrability provided by this UAV, it has a limited payload required to employ a compact camera for the photogrammetric flight.

For the UAV, an ultra-compact camera, Canon IXUS 115 HS, was chosen allowing a GSD of 1 cm.

Basically, the multi-data obtained by UAV and paratrike is geometrical, i.e., two point clouds with metric properties and texture information (RGB) which are homogenized under a common reference system based on a network of control and check points.

2.3. Terrestrial Laser Scanner

In those complex zones, a phase shift terrestrial laser scanner, Faro Focus 3D, was employed (

Table 3. Technical specifications of the terrestrial laser scanner (TLS), Faro Focus 3D.

Model	Faro Focus 3D
Principle	Phase Shift
Wavelength	905 nm (Near infrared)
Field of view	360° H × 320° V
Range std. deviation	2 mm at 25 m
Measurement range	0.19 mrad
Beam divergence	8 mm at 50 m
Scanning speed	976,000 points/s

).

According to the archaeological settlement characteristics and the TLS technical specifications, laser stations were established in a network that guaranteed an average spatial resolution of 5 mm for the whole scenario. The mean distance acquisition was 15 m.

2.4. Geo-Referencing System

The establishment of the mapping frame in the study area is performed with two GNSS bi-frequency (L1, L2) receivers, Topcon manufacturer. The GNSS observation method was real-time kinematic (RTK) getting a relative and absolute precision of 1 cm and 3 cm, respectively.

The coordinate reference system was comprised of the official coordinate system established by Spanish law, a Compound Coordinate Reference System (CCRS) integrated by horizontal CRS referred to ETRS89 geodetic reference system and UTM Zone 30 mapping projection (EPSG: 25830), and vertical CRS with geoid´s origin defined in Alicante (Spain) (EPSG: 5782). This mapping frame was materialized by a GNSS surveyal using natural features for “El Reguerón” area, which could be clearly identified on aerial images (e.g., corners of well-defined objects, small features with excellent contrast, etc.). In this way, we avoid artificial targets appearing in those more emblematic parts of the archaeological settlement. On the contrary and due to the large extension area, artificial targets were used for the full recording of the archaeological settlement (with the paratrike) in order to establish a network of control and checkpoints.

3. Methodology

Given the complexity of the archaeological site and the archaeological documentation requirements, the aerial and terrestrial data acquisition was planned in order to provide an integral and integrated recording of the site.

All the data from the different sensors and plataforms were processed according to the workflow outlined in Figure 2 and explained in the following subsections.

3.1. Data Acquisition

Each technique (i.e., photogrammetry, laser scanning) requires different procedures and protocols for data acquisition.

In the case of aerial data acquisition, image projection centers and camera attitude must be previously defined according to the classical photogrammetric parameters. The flight planning was done using the in-house software, MFlip and PFlip, for the UAV and paratrike flights, respectively [24]. The main difference between both flights was the type of flight planning: for the paratrike, a standard stereoscopic photogrammetric flight was undertaken (Figure 3a), whereas for the case of UAV, oblique and vertical images were also considered. Therefore, data acquisition was completely planned in order to get better results. In particular, for the UAV, a script was prepared for an automatic photogrammetric flight; whereas for the paratrike, the flight axis and flight height were planned and followed by the pilot. Both flights were planned with higher overlaps in order to get better results in the dense matching process and both flights were planned considering the relief variation using a public digital terrain model (DTM) [24].

The result shown in Figure 3b corresponds to the UAV flight planning, which was complemented with oblique aerial images acquired manually. As can be seen through the image footprints, the relief effects make it difficult to maintain a constant GSD, so the data acquisition is complemented with terrestrial images. For its part, the paratrike was used for the whole recording of the settlement due to the autonomy and sensors limitations of the UAV.

For the “El Reguerón” area a more detailed flight with the UAV was designed using the ultra-compact camera.

An important issue for the aerial data acquisition was the image sharpness, which affects the final photogrammetric reconstruction. This issue is affected by the camera parameters (i.e., aperture, shutter time and sensibility), the platform performance (i.e., flight speed and efficiency of the stabilized gimbal to absorb the paratrike vibrations) and the scenario illumination conditions. The flight was executed on a cloudy day to avoid shadows being projected, and the camera sensibility and shutter time were set up for these conditions. The aperture and the focal length were fixed constant to avoid variations in the internal camera parameters.

As commented previously, the high vertical reliefs of the site and the level of detail required in some areas (e.g., those constructive elements that integrate the defensive system of the entrance, details of the walls, etc.) entailed that aerial images were not suitable; terrestrial images were thus acquired with the full-frame reflex camera. In addition, terrestrial laser scans were used in those complex areas where photogrammetry could entail problems requiring a lot of images to enclose the whole geometry or due to the presence of textureless objects or materials. A network of 13 TLS stations (Figure 4) was designed for an average distance of 15 m with an average spatial resolution of 5 mm.

From Figure 4, it can be observed that two laser scanning campaigns were planned for two different areas. The first campaign (orange colour) was designed for recording the existing constructions using six stations. The second campaign (blue colour) was designed for recording the walls of “El Reguerón” and the adjacent environment using seven stations.

3.2. Data Processing

Handling multiple sensors requires multi-data processing approaches that take the point cloud as the basic unit. However, since all data must be integrated to generate a single model, it is necessary to homogenize the information, establishing filtering, simplification and optimization algorithms. Thus, the result obtained by each method is the corresponding point cloud that will be registered, triangulated and textured for generating a single 3D model with metric and radiometric capabilities.

3.2.1. TLS Filtering and Alignment

Data acquired with TLS were processed with commercial software FARO SCENE 5.4 [25]. The raw TLS scans were filtered removing noise and undesired information which could affect the alignment process. Automatic filtering has been applied in those more conflictive areas. Firstly, each scan was filtered according to a distance threshold (20 m) in order to remove distant points. Afterwards, two specific filters (intense-based filter and outlier filter) were applied. The former applies a reflectance threshold to remove those points with the lowest intensity, whereas the latter analyses the point and its environment (3 × 3) using distances. For instance, if there is a distance variation of 1 cm between the point and its neighbourhood for more than the 50% of the neighbourhood, the point will be removed. These thresholds should be tested adaptively by the user depending on the type and geometry of the area. Finally, the more delicate areas (i.e., walls or vegetation) have been filtered manually. For instance, vegetation between the blocks of the wall was identified and removed manually.

Since artificial targets were not used, scan alignment was done by a solid rigid transformation of an iterative closest point (ICP) technique [26]. This iterative process was applied in pairwise stations, so the final a priori error was computed on the basis of the number of stations and the technical specifications of the TLS (

Table 4. Simplification and optimization results for “El Reguerón” area.

	Simplified	Simplified and Optimized
Number of points	11,267,122	2,816,853
Number of triangles	22,532,754	5,635,313
Spatial resolution * (Min, Avg, Max)	(60.5, 64.3, 69.1) mm	(69.2, 81.1, 98.9) mm

* Confidence Interval = 1σ.

) reaching 1.3 cm for the assumption of two consecutive overlapping point clouds. The reference system of the aligned TLS point cloud is defined through a local Cartesian system corresponding to the first scan station with coordinates (300,000, 4,000,000, 1000).

3.2.2. Photogrammetric Processing

The generation of the dense point cloud, from both aerial and terrestrial images, was automatized through the Photogrammetry Workbench (PW) in-house software [27], following a three step workflow: 1. Image registration; 2. Camera orientation and 3. Dense matching.

The aerial images, coming from the paratrike and UAV, were processed to generate a hybrid model, with a total of 293 images and two different cameras. The dataset was checked to assure the sharpness of the images, which could be decreased by the motion blur.

• The feature extraction has been carried out by the ASIFT (Affine Scale-Invariant Feature Transform) algorithm [28]. As its most remarkable improvement, ASIFT includes the consideration of two additional parameters that control the presence of images with different scales and rotations. In this manner, the ASIFT algorithm can cope with images displaying a high scale and rotation difference, common in oblique images. The result is an invariant algorithm that considers the scale, rotation, and movement between images. The main contribution in the adaptation of the ASIFT algorithm is its integration with robust strategies that allow us to avoid erroneous correspondences. These strategies are the Euclidean distance [29] and the Moisan-Stival ORSA (Optimized Random Sampling Algorithm) [30]. This algorithm is a variant of Random Sample Consensus (RANSAC) [31] with an adaptive criterion to filter erroneous correspondences by the employment of the epipolar geometry constraints. Once the feature points have been extracted and described, the final matching points are assessed based on their spatial distribution on the CCD. An asymmetric distribution (radial and angular) of matching points regarding the principal point will affect the correct determination of internal camera parameters and also the image orientation. Therefore, if the matching points do not cover an area more than two-thirds of the CCD format, the user will be alerted in order to modify the detector (ASIFT) and descriptor (SIFT) parameters. Through this quality control we try to minimize problems associated with the weakness and common deficiencies in the photogrammetric network geometry of both aerial flights (UAV and paratrike).

This result provides the next expression:

$${\mathbf{A}}_F=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]=H_{\lambda }{R}_1(\kappa )T_1{R}_2(\varpi )=\lambda \left[\begin{array}{cc}\cos \kappa & -\sin \kappa \\ \sin \kappa & \cos \kappa \end{array}\right]\cdot \left[\begin{array}{cc}t& 0\\ 0& 1\end{array}\right]\cdot \left[\begin{array}{cc}\cos \varpi & -\sin \varpi \\ \sin \varpi & \cos \varpi \end{array}\right]$$

(1)

where A_F is the affinity transformation that contains scale, λ, rotation, κ, around the optical axis (swing) and the perspective parameters that correspond to the inclination of the camera optical axis, φ (tilt) or the vertical angle between optical axis and the normal to the image plane; and ϖ (axis), the horizontal angle between the optical axis and the fixed vertical plan.

In order to accelerate the process, the overlapped aerial images were identified by their approximate camera orientations provided by the navigation system. In the case of terrestrial images, an all-to-all comparison was applied. This sub-step is a time-consuming process which increases exponentially with the number of images [32].

2.

The multi-image protocol acquisition will require robust orientation procedures. For this purpose, a combination between computer vision and photogrammetric strategies was used. This combination is fed by the resulting keypoints extracted previously. In a first step, an approximation of the external orientation of the cameras was calculated following a fundamental matrix approach [33]. Later, these spatial (X,Y,Z) and angular (ϖ-omega, φ-phi, and χ-kappa) positions are refined by a bundle adjustment complemented with the collinearity condition [34]. In this field, several open source tools have been developed such as Bundler [35] and Apero [36]. For the present case study, both were combined and integrated. In particular, a specific converter has been developed for reading Bundler orientation files (*.out) and computing the three rotation angles and three translation coordinates of the camera in Apero. In addition, a coordinate system transformation has been implemented for passing from the Bundler to the Apero coordinate system. It is remarkable that at the same time, thanks to the reliability of the photogrammetric procedures used, it is possible to integrate as unknowns several internal camera parameters (focal length, principal point, and radial distortions). This possibility allows the use of non-calibrated cameras and guarantees acceptable results. For the present case study, a self-calibration strategy supported by a basic calibration model which encloses five internal parameters (focal length, principal point, and two radial distortion parameters) was used [37,38]. In order to provide metric capabilities to the model, manual identification of ground control points (GCPs) in the images were accomplished. Including these as an input in the bundle adjustment, the model is oriented according to the global coordinate system.

$$\begin{array}{l}(x-x_0)+\Delta x=-f\frac{r_{11}(X-S_X)+r_{21}(Y-S_Y)+r_{31}(Z-S_Z)}{r_{13}(X-S_X)+r_{23}(Y-S_Y)+r_{33}(Z-S_Z)}\\ (y-y_0)+\Delta y=-f\frac{r_{12}(X-S_X)+r_{22}(Y-S_Y)+r_{32}(Z-S_Z)}{r_{13}(X-S_X)+r_{23}(Y-S_Y)+r_{33}(Z-S_Z)}\end{array}$$

(2)

where x and y are the known image coordinates, X_i, Y_i and Z_i are the corresponding known GCPs, r_ij are the unknown 3 × 3 rotation matrix elements, S_X, S_Y and S_Z represent the unknown camera position, f is the principal distance, x₀ and y₀ are the principal point coordinates and ∆x and ∆y are the lens distortion parameters. These internal camera parameters may be known or unknown by the user and thus are introduced as equations or unknowns (self-calibration), respectively.

3.

One of the greatest breakthroughs in recent photogrammetry has been exploiting, from a geometric point of view, the image spatial resolution (size in pixels). This has made it possible to obtain a 3D object point of each of the image pixels. Different strategies have emerged in recent years, such as the Semi-Global Matching (SGM) approach [39] that allows the 3D reconstruction of the scene, in which an object point corresponds with a pixel in the image. These strategies, fed by the external and internal orientations and complemented by the epipolar geometry, are focused on the minimization of an energy function [39]. However, besides the classical SGM algorithm based on a stereo-matching strategy, multi-view approaches are incorporated in order to increase the reliability of the 3D results and to better cope with the case of complex archaeological sites (where the images are captured with different sensors). Considering the two types of flights performed (UAV and paratrike), two different multi-view algorithms were used. For the vertical flight (paratrike), the multi-view MicMac algorithm [40] was used. Meanwhile, for the oblique flight (UAV), the multi-view SURE algorithm [41] was used, which allows a complete reconstruction of the scene. Both strategies consist of minimizing an energy function throughout the eight basic directions that a pixel can take (each 45°). This function is composed of a function of cost, M (the pixel correspondence cost), that reflects the degree of the similarity of the pixels between two images, x and x’, together with the incorporation of two restrictions, P₁ and P₂, to show the possible presence of gross errors in the process of SGM. In addition, a third constraint has been added to the process of SGM; it consists of the epipolar geometry derived from the photogrammetry, and it can enclose the search space of each pixel in order to reduce the enormous computational cost. In that case, it will generate a dense model with multiple images, obtaining more optimal processing times.

$$E(D)={\displaystyle \sum _x\left(M\left(x,D_x\right)+{\displaystyle \sum _{x\text{'}\in N_x}{P}_{1}T\left[\left|D_x-D_{x\text{'}}\right|=1\right]+{\displaystyle \sum _{x\text{'}\in N_x}{P}_{2}T\left[\left|D_x-D_{x\text{'}}\right|>1\right]}}\right)}$$

(3)

where E(D) is an energy function that must be minimized on the basis of the disparity (difference of correspondence) through the counterpart characteristics, the function M (the pixel correspondence cost) evaluates the levels of similarity between the pixel x and its counterpart x’ through its disparity Dx, while the terms P₁ and P₂ correspond with two restrictions that allow for avoiding gross errors in the dense matching process for the disparity of 1 pixel or a larger number of them, respectively.

3.2.3. Data Fusion

Data fusion has been performed homogenizing the data provided for each sensor and generating a common product, a point cloud, with metric properties’ multi-resolution and photorealistic texture. Concretely, in order to fuse data, both flights (UAV and paratrike) were solved under a combined photogrammetric bundle adjustment using common control points and two different cameras. Through this combined bundle adjustment, a better and more homogeneous aerial photogrammetric point cloud is obtained, avoiding errors that would be obtained and propagated using a solid rigid transformation. In particular, the combined bundle adjustment (UAV-uav and paratrike-pt) is solved through a least square adjustment based on collinearity condition Equation (2), as follows:

$$\begin{gathered} x=x\left({\overline{c}}_{pt}, {\overline{c}}_{uav},{\overline{eo}}_{pt i}, {\overline{eo}}_{uav j},{\overline{X}}_k\right) \\ y=y\left({\overline{c}}_{pt},{\overline{c}}_{uav},{\overline{eo}}_{pti},{\overline{eo}}_{uavj},{\overline{X}}_k\right) \end{gathered}$$

(4)

where:

• ${\overline{c}}_{pt}and{\overline{c}}_{uav}$ are the camera vectors used for paratrike and UAV, respectively, and which include the internal camera parameters (principal point and focal length) and lens distortion coefficients (radial-K, decentering-P and affinity parameters-b). A total of ten unknowns were used for each camera vector, $\overline{c}=\left(x_0,y_0,f,K_1,K_2,K_3,P_1,P_2,b_{1,},b_2\right)$.
• ${\overline{eo}}_{pti}and{\overline{eo}}_{uavj}$ correspond with the six unknowns of the external orientation for paratrike and UAV images, respectively. Being the external orientation vector, $\overline{eo}=\left(S_x,S_y,S_z,\omega ,\phi ,\chi \right)$.
• ${\overline{X}}_k$ represents the spatial coordinates vector (X,Y,Z) of the unknown object points.

Therefore, the equation system is defined as follows:

$$Ax-K=V$$

(5)

where A corresponds with the design matrix based on collinearity equations and linearized through a first-order Taylor series, K is the observations matrix, V is the residual vector and x is the unknown’s vector solved through a least squares adjustment as follows:

$$x={\left(A^{T}PA\right)}^{-1}{A}^{T}PK$$

(6)

P is the weight matrix which corresponds with the inverse cofactor matrix of the observations. The equation system is solved through a twofold step: first by computing the exterior orientation parameters, $\overline{eo},$ and then computing the object points, ${\overline{X}}_k$, that represent the point cloud.

Next, the remaining step is the registration of the aerial point clouds derived from the different sensors under a global coordinate system. To this end, a GNSS campaign of 3 h based on three permanent ERGNSS stations was performed in order to provide precise coordinates to the GNSS base station used in the archaeological settlement. GNSS observations were processed guaranteeing an absolute error of 3 cm. The RTK surveying of the GCPs allowed us to obtain the aerial point cloud (i.e., coming from UAV and paratrike) under a global coordinate system (EPSG: 25830 and height EPSG: 5782), reaching a final relative precision of 1 cm.

Finally, regarding terrestrial laser scanner (TLS), the different scans were aligned and then co-registered with the photogrammetric point cloud coming from UAV and paratrike, using matching points defined manually as initial approximations. This was carried out using each dataset of coordinates in its coordinate system, that is, the local system in the case of the aligned TLS point cloud and the absolute system for the photogrammetric point cloud, and then applying a variation of the ICP technique, Least Squares Matching (LSM). A figure (Figure 5) to visually illustrate how this fusion has been done is included. A is the function that represents the point cloud coming from the aerial photogrammetry (UAV and paratrike) and B is the function that represents the aligned TLS point cloud, the registration of both point clouds will be obtained as follows:

$$A_i\left(x,y,z\right)-e_i\left(x,y,z\right)=B_j\left(x,y,z\right)i,j=1,\dots ,n,i\ne j$$

(7)

The equation that relates both models is a solid rigid transformation with seven parameters:

$$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\end{array}\right]=R\left[\begin{array}{c}{X}_B\\ Y_B\\ Z_B\end{array}\right]+\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]$$

(8)

where R represents the rotation matrix and T the translation vector. The correspondence between both point clouds is obtained through a minimisation, based on least squares adjustment, of the Euclidean distances between both point clouds. Since the rotation matrix is composed of non-linear functions, first-order Taylor series were used for the linearization of the Equation (8) as follows:

$$\begin{gathered} -e_i\left(x,y,z\right)=B_j^0\left(x,y,z\right)+\frac{B_j^0\left(x,y,z\right)}{\partial x_j}d{x}_j+\frac{B_j^0\left(x,y,z\right)}{\partial y_j}d{y}_j+\frac{B_j^0\left(x,y,z\right)}{\partial z_j}d{z}_j-A_i^0\left(x,y,z\right) \\ -\frac{A_i^0\left(x,y,z\right)}{\partial x_i}d{x}_i-\frac{A_i^0\left(x,y,z\right)}{\partial y_i}d{y}_i-\frac{A_i^0\left(x,y,z\right)}{\partial z_i}d{z}_i \end{gathered}$$

(9)

where:

$$\begin{gathered} dx=d{t}_x+a_{10}dm+a_{11}d\omega +a_{12}d\phi +a_{13}d\kappa \\ dy=d{t}_y+a_{20}dm+a_{21}d\omega +a_{22}d\phi +a_{23}d\kappa \\ dz=d{t}_z+a_{30}dm+a_{31}d\omega +a_{32}d\phi +a_{33}d\kappa \end{gathered}$$

(10)

Through an iterative process, the linear system $-e=Ax-l$ is solved, x being the unknown’s vector encompassing the transformation parameters ($d{t}_x,d{t}_y,d{t}_z,dm,d\omega ,d\phi ,d\kappa$), whereas l is the observation vector enveloping the discrepancies of Euclidean distances between both point clouds, and e is the residual vector.

Figure 5 graphically outlines the process followed for data fusion. An example of the aligned TLS point clouds for two scans of the wall (yellow and blue point clouds) in a local coordinate system is outlined to the left side of the image. Subsequently, the aligned TLS point clouds are fused with the photogrammetric point cloud in a global coordinate system (green point cloud) by means of LSM.

3.2.4. Post-Processing

Once a metric and geo-referenced point cloud was generated, a mesh strategy was applied to generate the digital surface model (DSM). In this case, the incremental Delaunay triangulation algorithm was applied [42].

To improve the model quality, break lines were incorporated as geometric constraints. Break lines were manually restituted by the operator. Its use is relevant for the accurate representation of significant slope changes, as well as for correctly defining the defensive walls of the site.

In those cases where images come from different methodologies or acquisition time, a radiometric adjustment was necessary to improve the final model visualization and avoid abrupt radiometric changes in the texture.

Once the DSM was obtained, it was possible to generate a true orthoimage from the oriented images. Since the multi-data source and multi-sensor approach provide different DSM resolutions, it was possible to generate true orthoimages with different spatial resolutions. For instance, a true orthoimage that depicts the complete archaeological site was generated with a resolution of 10 cm; whereas a 2 cm spatial resolution was employed for the settlement entrance area in order to appreciate the construction features.

3.2.5. Simplification and Optimization

Different simplification and optimization strategies were applied to the different geomatics products (DSM and orthoimages) in order to access and analyse 2D and 3D information through a Web service using mobile devices such as tablets or smartphones. In particular, a pyramidal structure was generated for the Web visualization of the 2D orthoimages and 3D models, including “El Reguerón” as the area with the most resolution. The methodology used for simplifying the 3D models is based on the strategy “smooth looking” [43] which consists in a straightforward mesh generation procedure applying the Poisson algorithm [44] to the hybrid point cloud. The mesh resolution is determined by the octree level, chosen according to the user-defined spatial resolution. Unlike a direct mesh generation process, which usually requires mesh editing operations such as filtering and refinement [45,46], Poisson algorithm directly encloses a smoothing step and provides a continuous geometry. This process could be controlled by some computation parameters as the minimum number of sample points that should fall within a node of the resulting octree. This parameter controls the loss of detail by the smoothing process. Although a high number of sample points implies a decrease of the number of mesh vertices, its geometry could be affected, being the final mesh mildly shrunk. For this reason, and as the noise level of the hybrid point cloud is low, the threshold value was kept low in order to avoid redundant smoothing processes. Then, a mesh optimization strategy, based on reducing the final number of triangles through the “collapse” of non-relevant areas without losing significant level of detail, is performed. For this purpose, an iterative process is proposed where the mesh derived from the previous procedure (smooth-looking) is “collapsed” by 5% of the total number of triangles using the quadratic edge collapse algorithm [47]. The algorithm essentially removes edges by merging and regrouping nearby vertices. With the aim of minimising the distortion of the surface geometry, it is necessary to establish a precision threshold to stop the iterative process. So that, if the resulted mesh error against the original input point cloud remains unchanged (with respect to the previous iteration), the collapse process continues. An iterative procedure is required since the quadratic edge collapse algorithm implementation does not allow decimation at fixed spatial resolution. Although our approach is based on an edge-collapse algorithm for the mesh optimisation, others approaches could be applied, such as the remeshing with recursive resampling as shown in [48] based on the Marching Cubes algorithm [49]. However, this remeshing approach exhibits higher deviations from the original model, up to 20 times higher than the quadratic edge collapse algorithm [48], and is therefore inadequate for our final archaeological products.

Finally, a texture mapping of the simplified and optimised mesh was performed using the commercial software Agisoft Photoscan^®.

4. Experimental Results

4.1. Area of Study

The city of “Tolmo de Minateda” (Figure 6) was a strategic settlement of great importance for several centuries, largely because of its peculiar topography and geographical location. It is placed on a plateau hill of about 50 m of height, located at the junction of the route from Complutum to Carthago Nova, one of the principal Roman routes connecting the interior of the plateau with the southeast coast, and a road connecting Castulo with Saetabis.

The archaeological research from the last 30 years has highlighted the importance of this site, revealing a history from the Middle Bronze Age, through the Iberian era, the Roman period, and the Middle Ages to the twentieth century. The Middle Age period provided most of the information by an important Visigoth settlement located in the upper part of the “El Tolmo”, where an important Christian basilica was found between houses and cemeteries.

One of the most interesting areas of the archaeological site is “El Reguerón” (Figure 6), an area of natural drainage 12 m in width, with a main entrance at the top of the hill where the city was located. In the settlement, an important fortification system consisting of three walls of different chronology and architectural typologies has been documented [51]. The oldest fortification is represented by the so-called “embanked” wall and was built during the final phase of the Iberian period (4th−2nd century B.C.). Currently, only the remains of a wall 6 m high and 10 m wide at the top, embanked in the external front and built in irregular masonry work, was preserved in the soil. During the archaeological excavation of this ancient wall, an earlier phase, which dates back to the Middle Bronze Age, was discovered. The Iberian wall was successively used as a retaining wall for a new fortification building during Roman times when the “Tolmo de Minateda” probably received the title of municipality.

The last fortification found in “El Reguerón” comes from the period of peninsula occupation by the Visigoths (5th–6th century A.D.). This wall is presented as a solid, L-shaped bulwark that encloses the valley and flanks the main road access to the city. It is here at this point that a monumental gate, probably formed by two solid towers of blocks, was once located. Only partial remains have been preserved. The wall is formed by a line of blocks with inscriptions and architectural elements from older constructions (among which are examples from the Roman period).

The relevant archaeological stratification of this area with different structures and building types requires of a multi-data source and multi-sensor approach that allows us to properly record and classify archaeological surfaces, thereby establishing an integration of topographic data with documents of archaeological excavations.

In order to fulfill the archaeological documentation requirements, the aerial data gathering was performed with a paratrike which allows us to enclose the whole extension of the archaeological site using a full-frame reflex camera, assisted by a specific gyrostabilized platform (MUSAS-MUltiSpectral Airborne Sensors) with a ground sample distance (GSD) of 3 cm. According to the GSD desired for the whole archaeological site and considering the camera specifications, a maximum flight height of 224 m was established. As a result, a total of 268 images along seven strips (NW-SE direction) were required for guaranteeing a side overlap of >30% and a forward overlap of >70%.

To record with higher spatial resolution the area of interest of the “El Regueron” site (Figure 7), an UAV with an ultra-compact camera (

Table 5. Geo-referencing errors.

Check Points			Discrepancies with 3D Model Coordinates
X (m)	Y (m)	H (m)	∆X (m)	∆Y (m)	∆XY (m)	∆H (m)
621,321.452	4,259,638.120	448.672	0.005	0.001	0.005	−0.013
621,375.933	4,259,646.226	454.775	−0.009	0.004	0.009	0.032
621,416.258	4,259,667.686	463.661	0.007	−0.011	0.013	0.030
621,419.022	4,259,658.232	464.618	0.001	0.008	0.008	0.046
621,345.185	4,259,659.272	454.267	0.002	0.018	0.006	0.024
621,321.298	4,259,660.311	452.106	−0.010	−0.026	0.028	0.036

) acquired the fortified walls with a GSD of 1 cm. Finally, 25 images and a flight height of 32 m guaranteed the desired spatial resolution. The overlap parameters were the same.

In addition, to increase the level of detail in the fortification walls and thus avoid the occlusions due to the effects of terrain relief (i.e., areas occluded or without information in the model generated by aerial photogrammetry), a combination of terrestrial photogrammetry and terrestrial laser scanner (TLS) was employed for improving the final 3D hybrid model. As a result, a total of 36 terrestrial oblique images and 13 TLS point clouds were acquired in order to avoid areas without information, guaranteeing a subcentimeter resolution.

4.2. Workflow

The application of the multi-data source and multi-sensor fusion workflow was tested in the study area of “Tolmo de Minateda,” which includes a total of 9 ha.

Firstly, the aerial images were processed by the automated photogrammetric approach and the cameras were self-calibrated according to Section 3.2.2. As a result of this process, the raw reconstructed point cloud is obtained which contains 4,271,354 points, while the TLS integrated point cloud reaches up to 1,314,136 points. Secondly, the geo-referencing of the 3D model was solved by the employment of 28 control points homogeneously distributed across the site (but with higher density in the interest area), while ten were used as checkpoints. Finally, the final multi-resolution 3D model after triangulation is obtained (Figure 8), encompassing the whole study area with a spatial resolution of 3 cm. A total of 10,542,505 triangles were obtained.

For the detailed area which integrates different data sources (i.e., terrestrial laser scanner and photogrammetry), a spatial filtering of 0.5 cm was applied to avoid areas with excessive point density caused by the different overlaps.

By this multi-resolution approach, the final inspection of the archaeological site could be adapted to any area, as illustrated with “El Reguerón”, where the subcentimete resolution achieved for the walls and construction (Figure 9b) was integrated with the rest of the archaeological site (Figure 9a). This multi-resolution capability opens a range of possibilities for a spatial analysis, settlement interpretation, pattern recognition, and the establishing of relationships among the elements (sites, artefacts locations, etc.).

In order to validate the final integration of the different data sources, a series of checkpoints were used. The different error components are shown in Table 4. The average vertical error (2.6 cm) was higher than the horizontal error (1.4 cm), as was expected for the GNSS technique. An average precision of 3 cm for the 3D vector error was obtained, which is statistical compatible with the expected a priori error of 3.1 cm composed by the model and GNSS check point errors.

To manage more efficiently the available information of the archaeological settlement, true orthoimages were generated as derived products (Figure 10), since they combined the photorealistic texture with the metric properties in an easy-to-use document for non-experts.

However, to allow a more complete interpretation and to make use of the potential of the generated 3D product, an integration with a 3DGIS tool, GeoWeb3D, was performed. In particular, the 3D geometry provided by the multi-sensor approach has been integrated with 2D archaeological archives such as sketches, pictures, part details, etc., thus enhancing the subsequent analysis and providing problem solving and decision-making capabilities. Figure 11 shows the integration of different historical events in the reconstructed 3D model.

The integration of 2D information and 3D models allows us to extract intangible information that improves the analysis capabilities of the archaeological settlement. For instance, Figure 11 integrates an archaeological sketch of defensive constructions (2012) with the generated hybrid 3D model (2015). It should be noted that the high accuracy and proportion of the sketch is a perfect coincidence of the main homologous entities. However, some detached blocks belonging to the wall of the 3D model do not appear in the sketch. Analysing the position of blocks in the 3D model, it seems that they were spread along the natural drainage bed of the Tolmo, possibly suggesting that a runoff flow took place from the upper part of the Tolmo, providing this current block distribution.

Next, with the aim of showing the potential of analysing in situ the archaeological settlement, “El Reguerón” was tested using a mobile device (smartphone). The simplification and optimisation level was fixed considering different aspects: (i) the amount of 2D information which should be integrated; (ii) the minimum size of the interest elements and (iii) the level of reduction. A minimum resolution of 10 cm was established for the “El Reguerón” area.

Looking to the smallest and most emblematic elements of “El Reguerón,” it can be confirmed that these correspond to the blocks of the wall, which are bigger than the minimum resolution, meaning that they should be correctly represented. Subsequently, an iterative optimization process of the mesh was performed maintaining relevant information and stopping when the threshold of resolution (10 cm) was surpassed. As we can see in Table 5, the simplified and optimised models maintain the initial resolution (10 cm), with a level of reduction of 75% in comparison with the simplified model (without optimisation).

The simplification and optimisation procedures have removed 31,039,445 points of the 33,856,298 initial points. However, in order to assess the final error of the simplification and optimisation strategy, we have compared the simplified and optimised model with the original point cloud, obtaining an average error of 0.2 mm and a standard deviation of ±22 mm.

As we can see in Figure 12, higher discrepancies are located in peripheral areas due to the presence of vegetation with an irregular typology. Conversely, the walled area exhibits minimum discrepancies always less than ±5 mm.

Finally, results were presented through the Web based on the Open Source library Cesium. A specific template based on HTML language was prepared to show the 3D models and additional information using external geospatial services such as WMS, MapServer, Google Earth or Bing, among others. Thanks to the flexibility and portability of the mobile devices, it was possible for the archaeologist to interact directly with the platform at the field: recording data and adding new information with corresponding attributes and descriptions. An example is outlined in Figure 13, where the application is loaded in a smartphone BQ Aquaris E4.5 using Android 4.4.2 and Google Chrome 45. Results were incorporated to the archaeological information system and the spatial data infrastructure of the archaeological cultural heritage of Castilla La Mancha (ideARQ + SIA) [50].

Additionally, the optimal value of points and triangles for the proper management of 3D models in smartphones was examined. For this aim, three levels of simplification were analysed using three, two and one million(s) triangles, respectively. Afterwards, loading and operation times together with RAM were monitored in order to see the best simplification level for a conventional smartphone.

Table 6. Times and resources required with lighter models using mobile devices.

Model	Number of Points	Number of Triangles
3M	1,408,636	2,816,156
	Load	Operation	RAM consumption
	10 s	2–4 s	2.39 GB
2M	986,150	1,971,309
	Load	Operation	RAM consumption
	8 s	1–3 s	2.24 GB
1M	493,198	985,653
	Load	Operation	RAM consumption
	3 s	Instantaneous	1.79 GB

outlines the results obtained.

As we can see, with less than one million triangles, very good loading and operation times (instantaneous) are obtained, as well as an affordable RAM consumption for the majority of smartphones. Therefore, for optimal user experience and navigation for each visualization level (defined in terms of spatial resolution; see Section 3.2.5), the maximum display area is cropped accordingly to this empirical number of triangles.

5. Conclusions

This paper presents a methodology based on a combination of multiple sensors, platforms and techniques, which has been tested in a complex archaeological site. As is shown in the experimental results, the automation provided by the photogrammetric and laser techniques, along with the versatility of the aerial platforms (paratrike and Unmanned Aerial Vehicle UAV), provide the suitability of this methodology for complex archaeological sites.

The potential of ultralight aerial platforms (paratrike) is highlighted, due to its payload and flight autonomy, which overpasses the UAVs capabilities. Regarding the usual archaeological surveys, one of the main differential factors in this work has been the integration of aerial images at different resolutions, terrestrial images and terrestrial laser scans.

The generation of high resolution products based on photogrammetry requires high computational costs. To this end, the presented algorithms, such as Affine Scale Invariant Feature Transform ASIFT, have been developed to take advantage of Graphical Processing Unit GPU capabilities, reduce time operations and consequently improve workflow efficiency.

The management of the final models through a 3DGIS tool opens new analysis capabilities for the archaeologist (e.g., analysis though time, archaeological investigations, integration of historical and geometric models).

Finally, it has been demonstrated that through the simplification and optimization strategies, complex hybrid 3D modes which enclose 2D information can be flexibly shown through the Web and even embedded in smartphones. As a result, we can interact directly with hybrid 3D models using mobile devices, allowing the recording, storage and analysis of data in situ. The implemented methodology to optimize the visualization and interaction of geomatics products has proven to be effective. However, it lacks a fully complete user interaction, since the final navigation is constrained by the definition of the visualization level. To overcome this limitation, a future work line is aimed to a more efficient visualization management by the 3DGIS tool, avoiding the definition of level of visualization by the final user.

The technical case outlined in this paper could be of great interest for different stakeholders:

• Researchers for the interpretation, spatial and temporal analysis of the archeological settlement thanks to the integration capabilities and portability of the system.
• Managers through the monitoring of the archaeological settlement through time, the diffusion of the site using videos, documents, etc.
• Students who could exploit the didactical possibilities of the 3D inspection, interaction and superposition of thematic information.
• General public allowing a flexible and enjoyable accessibility to the archaeological settlement which complements and provides added value to a visit to an historical site.

The data acquisition and processing methodologies from Geomatic science broaden the possibilities of sensors, configurations and/or combination as future perspectives. Concretely, the resolution of the external orientation of the images directly, just using the integrated Global Navigation Satellites System GNSS/IMU Inertial Measurement Unit of the aerial platform, would allow us to speed up the field work and processing time. Furthermore, in order to obtain realistic textured models and orthoimages, the use of techniques such BRDF (Bidirectional Reflectance Distribution Function) to improve the radiometric parameters will be studied.