Multiplatform-SfM and TLS Data Fusion for Monitoring Agricultural Terraces in Complex Topographic and Landcover Conditions

Abstract

Agricultural terraced landscapes, which are important historical heritage sites (e.g., UNESCO or Globally Important Agricultural Heritage Systems (GIAHS) sites) are under threat from increased soil degradation due to climate change and land abandonment. Remote sensing can assist in the assessment and monitoring of such cultural ecosystem services. However, due to the limitations imposed by rugged topography and the occurrence of vegetation, the application of a single high-resolution topography (HRT) technique is challenging in these particular agricultural environments. Therefore, data fusion of HRT techniques (terrestrial laser scanning (TLS) and aerial/terrestrial structure from motion (SfM)) was tested for the first time in this context (terraces), to the best of our knowledge, to overcome specific detection problems such as the complex topographic and landcover conditions of the terrace systems. SfM–TLS data fusion methodology was trialed in order to produce very high-resolution digital terrain models (DTMs) of two agricultural terrace areas, both characterized by the presence of vegetation that covers parts of the subvertical surfaces, complex morphology, and inaccessible areas. In the unreachable areas, it was necessary to find effective solutions to carry out HRT surveys; therefore, we tested the direct georeferencing (DG) method, exploiting onboard multifrequency GNSS receivers for unmanned aerial vehicles (UAVs) and postprocessing kinematic (PPK) data. The results showed that the fusion of data based on different methods and acquisition platforms is required to obtain accurate DTMs that reflect the real surface roughness of terrace systems without gaps in data. Moreover, in inaccessible or hazardous terrains, a combination of direct and indirect georeferencing was a useful solution to reduce the substantial inconvenience and cost of ground control point (GCP) placement. We show that in order to obtain a precise data fusion in these complex conditions, it is essential to utilize a complete and specific workflow. This workflow must incorporate all data merging issues and landcover condition problems, encompassing the survey planning step, the coregistration process, and the error analysis of the outputs. The high-resolution DTMs realized can provide a starting point for land degradation process assessment of these agriculture environments and supplies useful information to stakeholders for better management and protection of such important heritage landscapes.

1. Introduction

Over the last two decades, our ability to monitor and characterize landscapes [1,2] through submetric digital terrain models (DTMs) has greatly improved due to important developments in high-resolution topography (HRT) techniques, methods, sensors, and platforms [3,4]. The explosion in the availability of HRT data is revolutionizing the way in which we can analyze features in the landscape and quantify many processes at the fine spatial resolutions at which they occur [5].

The choice of the most appropriate platform for HRT surveys must first consider the features present in the analyzed area, the required resolution, the cost and flexibility of the available and required technologies, and the spatial extent of the study area [6]. Each monitoring platform has advantages and disadvantages. They must be assessed during the survey planning phase, considering the aims of the study and the characteristics of the system under investigation. Among the most widely used technologies, airborne laser scanning (ALS) is the most suitable system for the acquisition of HRT data over remote and vegetated areas over short time periods across large areas [7,8,9,10,11]. However, this approach is costly for small area surveys. Instead, ground-based approaches might be more suitable on some sites with steep slopes and near-vertical surfaces, such as agricultural terraces, river banks, and landslides [12]. For example, ground-based platforms such as terrestrial laser scanning (TLS) are used when a higher resolution and more flexibility in the scanning angle is needed [13]. Other traditional techniques such as Global Navigation Satellite System (GNSS) and total station survey are not ideal for large areas and can thus prove costly [14]. Instead, recent developments in digital photogrammetry (e.g., structure from motion (SfM)) paired with multiview stereo (MVS) algorithms (hereafter together referred to as SfM) have increased the possibility of low-cost and high-resolution surveys, both ground- and aerial-based [15,16], opening new opportunities for landscape investigation and management [17].

The sole use of any of the above techniques can be a significant problem in rugged topography, steep slopes, inaccessible areas, zones with high altitude differences, and in vegetated environments, because each HRT method is constrained by several factors. Therefore, data fusion from different acquisition platforms is necessary for the creation of detailed DTMs in complex situations. The integration of different types of data (in terms of source, resolution, accuracy) acquired by different techniques allows the capture of more detailed information and valuable feature representation relating to the surveyed areas [18,19]. Among all the morphologically complex environments, in which data fusion could be an effective method, it is possible to identify those where the anthropic component has modified the natural landscape by introducing very different features compared to the surrounding context. An excellent example of such heterogenic systems is agricultural areas, where significant natural landscape changes have taken place, especially where there are steep slopes. Indeed, hilly and mountainous areas have been widely transformed for cultivation, either through the maintenance of natural slopes or the construction of terraces (using dry-stone walls or earth banks to reduce slope gradient, mitigate soil erosion, collect rainwater run-off, and support agriculture [20]). In this particular environment, some factors that constrain HRT surveys must be considered: (i) the existence of subvertical surfaces (i.e., terrace risers or walls) that are not easy to survey with all HRT techniques; (ii) the presence of vegetation (common in agricultural areas or in ancient, often abandoned, terraces) that covers surfaces and occludes the instrument’s point of view; (iii) the presence of steep and rugged slopes that, while not being more suitable for ground-based approaches, also restrict accessibility (e.g., for large instruments); (iv) the presence of anthropic infrastructures close to, or inside, the area of interest that impede the execution of a survey (e.g., no-fly zones, inaccessible areas, critical zones or parts with specific restrictions by flight authorities [21]); (v) the existence of cultivated areas which require surveys at a landscape scale, such as high change in elevation. The monitoring and characterization of these environments is important both in order to learn from the past and to study future landscape evolution. Agricultural terraced landscapes worldwide can be considered as part of our historical heritage (e.g., UNESCO or Globally Important Agricultural Heritage Systems (GIAHS) sites) and thus provide cultural ecosystem services, in addition to providing food security and income with significant benefits to the economy (e.g., terraced cultivation of wine grapes, fruits, and olives). Nowadays, with climate change leading to enhanced precipitation intensity [22,23], combined with land abandonment, lack of maintenance, and in some cases, unsuitable agronomic practices, such areas are under threat. Thus, terrace systems becoming one of the most erosion-prone agricultural landscapes [24,25]. Indeed, an increasing occurrence of slope failure (e.g., landslide, debris flow, and dry-stone wall collapse) has been documented, with direct consequences for people when these processes are triggered in densely populated areas [26]. Observing the evolution of agricultural landscapes cultivated in hilly and mountainous areas, often with terracing practiced, through HRT gives us important information on how to protect these sensitive environments [25,27,28]. HRT can provide suitable information through the extraction of 3D models, profiles, sections, scaled plans, and orthomosaics, simplifying and speeding-up the field geomorphological analysis [29].

The objective of this work is to test data fusion of HRT techniques (TLS and aerial/terrestrial SfM) in a particular agricultural context where, to the best of our knowledge, it has never been applied. The aim is to assess whether the exploitation of the SfM–TLS data fusion methodology applied in these anthropogenic landscapes can overcome specific detection problems such as complex topographic and landcover conditions. The terrace environments represent a new challenge for the data fusion methodology, which could be a key solution for a detailed reconstruction of terraced systems even when covered by vegetation. Such in-depth knowledge will provide a solid basis for land degradation process assessment and allow for a better management of relevant heritage terraced landscapes. In order to apply this methodology, the design of an appropriate data-fusion workflow was developed to plan surveys in difficult conditions (i.e., presence of vegetation, rugged topography, inaccessible areas, and the requirement to adequately locate topographic SfM and TLS targets) and to carry out an accurate merge of data in two agricultural terrace systems with diverse survey issues.

2. SfM–TLS Data Fusion

The data fusion concept identifies a set of methods and tools for using data derived from various sources with different characteristics [18]. In HRT, this is the process of merging data derived from different sensors but representing the same real-world object in order to produce a consistent, accurate, and useful representation and thereby overcome the limitations of single data sources [19,30]. There are different levels of data integration: low, intermediate, or high. These levels are dependent on the processing stage at which fusion takes place [18,31]. In this paper we used the low level of data fusion (i.e., the combination of several sources of raw data to produce new raw information that is more representative than the original raw inputs). We integrated raw point clouds of TLS and unmanned aerial vehicle (UAV)/terrestrial SfM to exploit the intrinsic advantages and overcome the weaknesses of each dataset. A data fusion approach, namely SfM–TLS merge, has been applied in cultural heritage environment and archaeological research [19,32,33,34,35,36] and in environmental monitoring [14,30,37,38,39]. This suggests that SfM–TLS fusion achieves better performance in various applications compared to the use of a single data type [40]. In previous research, data fusion was applied in specific complex landscapes comprising built and natural heritage within demanding areas [32,34]. In terrace systems, the complex topography and vegetation cover certainly makes it more difficult to apply the usual data fusion. For this reason, steep agricultural terrace systems (Figure 1) are a challenging environment for the testing of the data fusion approach, which to the best of our knowledge has yet to be applied in this particular context (single-mode HRT has been applied thus far to survey terrace environments [20,27,41,42,43]). The SfM–TLS merge is typically done to overcome weak points of each technique, e.g., lack of texture, gaps due to occlusions, and noncollaborative material/surfaces [19], derived from the intrinsic nature of the HRT survey method. Therefore, it was important to identify these particular aspects within the survey project.

TLS produces surveys that can reach millimeter accuracy and is considered to more reliably filter vegetation spots [44], which may be problematic in SfM data due to difficulties encountered from image matching over plant cover [37]. Indeed, TLS produces more reliable information on intercanopy surfaces and the ground as it exploits different points of view at close range to the studied object, where vegetation density is low. TLS can use the ability to provide additional information such as the intensity classifications of the surveyed surfaces (ground or canopy [45]). However, TLS is expensive and requires many hours of fieldwork (e.g., scanning time, or TLS target and tripod moving for multiple scans [46]). A TLS survey usually needs multiple scans to avoid occlusions due to the areas obscured in the sensor’s field of view [39,46,47]. Therefore, the collected point cloud can often be incomplete without accurate survey planning [48]. TLS surveys typically capture information on vertical surfaces (e.g., terrace walls, steep slopes, landslides, building facades) but with sampling limits across the upper zones within the field of view (e.g., top of vegetation canopy, building roofs [32]). Aerial SfM can overcome these data gaps through UAV systems, which have evolved significantly in the last decade. UAVs have revolutionized large-scale low-altitude imaging and geospatial information gathering [49]. The impact of UAVs has been noticeable in remote sensing [50,51], supported by the increasing use of SfM and linked to the development of user-friendly SfM software [52]. This low-cost technique with nadir UAV images minimizes the shadows of the TLS ground-based view [53] and allows the surveying of large areas which are not, or are less, accessible with the scanner in a short time. Vertical surfaces pose a problem for the UAV’s nadir viewpoint, but SfM allows the choice of a wide range of other acquisition platforms [54], including ground-based photos. They can also capture the fine-grained and otherwise hidden details of a vertical surface. Terrestrial SfM can also be a good solution when flight permissions are not granted by local legislation at certain sites or when weather conditions may not be suitable for drone flights (e.g., high wind speeds [12]). However, the ground elevation under vegetation is usually not represented in SfM, limiting the quality of the resulting point cloud [55]. Another constraint for SfM surveys, both in terms of time and cost [56], especially in steep-slope and rugged environments, is the ground control point (GCP) requirement. This is of fundamental importance for indirect georeferencing or indirect sensor orientation (InSO) [57] and the coregistration process [58]. If the surveyed zone is inaccessible, densely vegetated, hazardous, or subject to restrictions, locating and measuring GCPs with GNSS becomes a problem [59,60]. Today, this problem can be solved through the direct georeferencing (DG) method [61]. In DG or direct sensor orientation (DSO) [57], the photogrammetric solution is determined with precise and accurate knowledge of the camera positions/orientation at the time of image acquisition [62]. The exterior orientation parameters of aerial images, at the exposure time, are computed by processing the onboard GNSS and high-quality inertial measurement unit (IMU) data [21]. Nowadays, onboard multifrequency GNSS offers a real-time kinematic (RTK) option that allows a high level of image geotagging with on-site data quality checks. However, recording raw GNSS data is also desirable, as a backup in the case of radio link (i.e., between a master station and the UAV GNSS) problems, in order to refine the camera station positions (i.e., precise ephemeris data of GNSS satellites are available during postprocessing), and to consider the UAV trajectory and atmospheric error sources. Consequently, postprocessing kinematic (PPK) solution can often provide a more accurate solution than RTK [57,63,64,65].

Therefore, the SfM and TLS techniques have been chosen for their intrinsic characteristics and their ability to complement each other in terms of overcoming problems (Figure 1) [66]. In addition, research has shown how TLS and SfM surveys can lead to comparable results in terms of accuracy [46,67,68,69,70]. For this reason, SfM–TLS data fusion is a logical option in challenging contexts (e.g., rugged vegetated terrace complexes) in order to increase the reliability of the information and generate an accurate and complete 3D model of the entire landscape.

In the data fusion process, one fundamental step is coregistration. Even if several studies highlight how coregistration was essential in order to merge HRT data correctly, it is often not addressed in postprocessing workflows [71]. In addition, in topographically complex environments such as terraced areas, coregistration requires specific measures to overcome misalignments due to the presence of vegetation. The SfM–TLS data fusion needs the coregistration of TLS multiple scans (a common task in TLS surveys [72]) because the overlapping point clouds have to be located in the same reference system, computing a transformation (composed of rotation, translation, and possibly scale) that defines the optimal alignment among all of them. Furthermore, the alignment of partially overlapping TLS and SfM point clouds is also necessary; these must be coregistered by computing a similarity transformation that minimizes the distance between corresponding points.

3. Data and Methods

An integrated approach employing ground-based and UAV SfM with TLS data permitted an accurate survey of the highly vegetated terrace systems in order to produce detailed DTMs. Effective data fusion has to consider several aspects, starting from a careful survey planning, where the distribution of the SfM and TLS targets is crucial for coregistering data in the final merge. For this reason, it is important to develop a precise workflow (Figure 2). The workflow must address all the possible issues, such as the coregistration problems and the errors, that the proposed methods could generate in terms of DTM accuracy. In the following sections, we give full details on HRT surveys (i. and ii. in Figure 2), data processing (iii. to v. in Figure 2), data coregistration (vi. in Figure 2), data fusion (vii. in Figure 2), DTM generation (viii. in Figure 2), and error analysis (ix. in Figure 2).

3.1. Study Areas

We tested this approach on two agricultural terrace sites (Figure 3a): Soave traditional vineyards (a Globally Important Agricultural Heritage Systems (GIAHS) site) in the Veneto region of northeastern Italy and Martelberg in the Saint-Martens-Voeren area of eastern Belgium. These are case study sites for the TerrACE geoarchaeological research project (ERC-2017-ADG: 787790, 2019–2023; https://www.terrace.no/) that aims to create a methodological step-change in the understanding of agricultural terraces. This research project is applying a new scientific integrated methodology to agricultural terraces across Europe, bringing together landscape archaeology, geomorphology, and paleoecology. We also chose these sites because they both present complex topography and landcover, including the presence of vegetation that covers parts of the subvertical surfaces (e.g., vertical walls of terraces), slopes, and large survey areas.

3.1.1. Belgian Study Area

The Belgian site (18 ha, elevation range between 148 and 183 m a.s.l., and average slope of 20%; see Figure 3b) is located in a hilly agropastoral area close to the border with the Netherlands. The site is characterized by six braided lynchets (terraces without stone walls) covered with trees and shrubs with dense canopy, interconnected by sloping pastures that are periodically plowed. The study area spans on two different hills divided by a small road which crosses them in the lower part, forming a small valley. The studied terrace complex is located on the west side of the valley. The case presented here is particularly interesting from a geoarchaeological point of view. Indeed, geomorphological researchers had applied a tillage-erosion model to the site, which had yielded an estimate of its age—but this had not been tested [73]. It also represents and a challenging survey in terms of HRT data fusion.

3.1.2. Italian Study Area

The Italian site (3.5 ha, elevation range between 33 and 70 m a.s.l., and average slope of 35%; see Figure 3c) is located in a hilly zone which is famous for the production of excellent wines. The presence of steep slopes led local people to invest time and effort in the construction of contour terraces—a system that is required for the creation and maintenance of vineyards. The survey comprises several abandoned dry-stone terraces in the lower part. Some of these terraces have collapsed, and shrubby vegetation has invaded the flat terrace areas (Figure 3c). This obscures the view of the terraces, thus complicating the HRT survey. In the higher part of the study area, where the slope is less steep, vineyards are still under cultivation.

3.2. Data Acquisition

3.2.1. SfM Surveys

An integrated approach employing terrestrial and UAV SfM was carried out on both sites. Ground-based photography was used to survey particular elements, as the vertical parts of the terrace risers were obscured or collapsed. Meanwhile, UAV images captured flat terrace zones at a large spatial scale (Figure 1). Since the study areas were large and included considerable variation in slope morphology, they were divided into homogeneous elevation zones that were surveyed through planned and manual UAV flights. Before image acquisition, GCPs employed in SfM were distributed throughout the study areas (Figure 4a,b; Section i. in Figure 2) in flat zones but also on vertical surfaces so that GCPs were visible from different points of view (i.e., nadir and oblique images) and easily distinguishable from TLS positions. These were employed as corresponding points in the SfM–TLS point cloud coregistration. In the survey planning, the number, location, and distribution of GCPs are crucial [56,58,74] because a reliable SfM adjustment during 3D reconstruction demands a robust GCP network [75]. The GCP network was based on the features and spatial scale of the studied area and the desired resolution for the survey. The SfM targets comprised square orange cloth with a black rhombus in the middle, an old compact disk placed in the center (i.e., CD-ROMs; see Figure 4e). SfM arrangements are summarized in

Table 1. Main characteristics of SfM acquisitions (ground-based and UAV).

Study Area	Area Covered (ha)	Number of Images	Flight Height (m)	GSD (m)	Number of GCPs (CPs) 1	GNSS Positional Accuracy (X, Y—Z) (m)
Belgium	18	1219	20–45	0.005–0.015	40 (13)	<0.05
Italy	3.5	632	25–35	0.006–0.008	56 (18)	0.03–0.04

¹ Following [76], we applied a bootstrapping resampling technique to randomly select one-third of the GCPs and use them as control points (CPs) (see Section 3.6) to assess the uncertainty of each point.

.

Belgian site: In October 2019, a custom-built Quadcopter equipped with a DJI A3 flight controller and a Sony ILCE-6000 camera (24 Mpixels, focal length 16 mm, sensor size APS-C 23.5 mm × 15.6 mm) was used for aerial SfM survey (Figure 4c; Section ii. in Figure 2). The UAV was fitted with a compact multi-GNSS RTK receiver (Emlid Reach M+, Emlid Ltd., Hong Kong; Tallysman TW2708 antenna model) with RTK–PPK capability to log the raw data as UBX format using global positioning system (GPS) and GLONASS satellites. The offset shift between the antenna phase center (placed on an aluminum plate, 22.5 cm right above the camera lens center) and camera projection center was considered during the postprocessing, assuming a constant vertical offset (more details can be found in [65], where the same system was used). Moreover, an electronic system was used to integrate and synchronize the GNSS with the action camera for image geotagging. The images contained only positioning information without attitude parameters because a link between the UAV-IMU and camera was not available. During the UAV flights, a Reach RS (Emlid Ltd.) base station was mounted on a tripod located in the middle of the survey area (the maximal distance between the UAV and the base station was 400 m) to provide positioning correction input. The receiver of the base was configured to log the raw data in a RINEX file at 5 Hz. Flight missions (i.e., three flights, namely the west slope, east slope, and central part; see Figure 4a) were planned using the Autopilot app (Hangar Technology, 2018) that defined the frontal (90%) and side overlap (80%) based on the speed of the UAV and the camera trigger interval, which was set at 3 s. For the ground-based survey, the same camera with the same focal length as the UAV survey was used to minimize the integration problems of aerial and terrestrial images in the SfM models [77]. GCPs were distributed only along the western side of the study area because the eastern side was inaccessible (Figure 4a). For this reason, the direct georeferencing technique was used to reconstruct the morphology of this specific part of the study area. The GCPs were surveyed using a Reach RS (Emlid Ltd.; RTK solution) with the EUREF-IP network in NRTK (network of permanent GNSS stations) mode. The point coordinates were tied to the WGS84/UTM zone 31N (EPSG: 32631) reference system.

Italian site: In December 2019, nadir and oblique UAV images were collected with a DJI Zenmuse X4S camera (20Mpixels, focal length 8.8 mm, 1-inch CMOS Sensor; see Figure 4f) mounted on a professional quadcopter (DJI Matrice210v2). The survey area was divided into two zones (i.e., the lower part with abandoned terrace and the higher cultivated ones; see Figure 4b) of uniform elevation. The UAV flight control unit (coupled to a GNSS) was used to plan two UAV flights, using software that adjusts the height and speed of flight accordingly, and the image overlap (optimal overlap of 80% in flight direction and a flight strip overlap of 70%). In areas where a specific detail was needed, the manual flight mode was used with a time-lapse function of the camera that allowed the capturing of an image (nadir and oblique) at 3 s intervals. This was sufficient to guarantee the overlap in sequential photographs, which is essential for the image matching algorithms used in SfM [37]. In addition, terrestrial SfM surveys were carried out using a Sony NEX-5R camera (16 megapixels, focal length 16 mm, sensor size 23.5 × 15.6 mm). Ground-based photos were taken in front of some specific wall collapses (i.e., hidden by vegetation), using an adequate average depth distance from the studied objects. A mean baseline of 0.30 m was used between adjacent camera positions to avoid large jumps in scale. A GeoMax Zenith 40 GNSS allowed us to survey all target types in NRTK mode. All the points coordinates were referenced to the RDN2008/UTM zone 32N (EPSG: 7791) reference system.

3.2.2. TLS Surveys

In order to enhance the speed and efficiency of data collection, the scanner positions were chosen so as to minimize their number. These choices considered the areas not reached by SfM (Figure 1). The data acquisition started with the planning of target and TLS positions. We used four tripod-mounted Leica targets (i.e., planar HDS 4.5-inch circular black/white targets; see Figure 4e) distributed within the survey area (Figure 4a,b) in order to make these targets visible at a certain distance from all points of view of the TLS positions (Section i. in Figure 2). For each scan, four targets were positioned, three of which remained in the same position for the subsequent scan, while one was moved from time to time in order to detect the adjacent area. This was important to maintain an overlap among the scans and an adequate number of targets for multiscan coregistration. An overlap of 30–50% between adjacent scans ensured the generation of a very accurate 3D laser model [78]. The riser surfaces of the terraces were surveyed with a Leica P50 TLS (Figure 4d). This scanning system allows a larger field of view (360°H × 279°V), high measurement accuracy (e.g., 1.6 mm @ 10 m at >1 km maximum range mode), a range of up to 1 km, and a scan rate of up to 1,000,000 points/s. Furthermore, besides intensity (e.g., in Figure 4d) of the reflected beam, the Leica P50 acquired a panoramic photo of the surveyed objects through an integrated camera (4 megapixels) after each scan and made a higher density specific scan of the targets, manually identified and labeled by the user on the P50 screen preview. The position of the targets was surveyed using the same GNSS of SfM surveys for the two study areas (Section ii. in Figure 2).

Table 2. Main characteristics of TLS acquisitions. The chosen parameters were an acceptable compromise between the level of detail of the final 3D model and computing resources needed for data processing.

Study Area	Number of Scans	Mean Distance between Scans (m)	Accuracy (m)	Maximum Range Mode (m)	Scan Rate (Points/s)	Image Resolution (Mpixels)	Number of Targets
Belgium	7	70	1.6 mm@10	120	1,000,000	4	20
Italy	6	65	1.6 mm@10	270–570 1	250,000–500,000 1	4	16

¹ Different parameter for the two different lines of TLS scans (Figure 4b).

synthesizes the characteristics of TLS surveys.

Belgian site: In order to survey the vertical vegetated risers of the terrace, several scans were required. This varied depending on the complexity and the length of the topographic surface: three different scan positions were used for the upper terrace (TLS 1, 2, 3 in Figure 4a), and two scan positions each were used for the middle (TLS 4, 5 in Figure 4a) and lower (TLS 6, 7 in Figure 4a) terraces. The scanner was located in front of the braided lynchets and moved along the terrace line, keeping approximately the same elevation.

Italian site: The terrace walls covered by shrubby vegetation were surveyed through six TLS scans located at two different lines of distance: three at a distance of 150 m (TLS 1, 2, 3 in Figure 4b) and three at a distance of 35 m (TLS 4, 5, 6 in Figure 4a) from the terrace system.

3.3. Data Processing

3.3.1. SfM Processing

For the Belgian site, raw GNSS data from the UAV-mounted cameras and the base station were extracted and corrected by postprocessing using RTKLib, an open-source software package for differential positioning computing. Then, the camera position estimates were extracted through PPK GNSS solutions and were available to test direct georeferencing in the SfM workflow.

The image datasets were processed with a 2× Intel Xeon Bronze 3106 CPU @ 1.70 Ghz with 256 GB RAM and 2× NVIDIA GeForce RTX 2080 Ti, through Agisoft Photoscan Pro v 1.4.5. The computer vision routines of SfM and MVS algorithms require powerful processers to extract the 3D point clouds from the images, and, additionally, orthomosaics. The first step in SfM processing workflow was masking unwanted objects (e.g., vegetation and clouds in ground-based images) in the photos uploaded in the software. This was followed by a camera precalibration using Agisoft Lens, an automatic lens calibration routine which uses an LCD screen as a calibration target (i.e., checker-board pattern). This tool estimated camera model parameters and lens distortion coefficients, ready to import into Photoscan (Section iii. in Figure 2). This precalibration step was a starting point for the parameter refinement in the next process i.e., the SfM step (Section iv. in Figure 2), where ground-based and UAV (nadir and oblique) photos were processed together in Photoscan to avoid subsequent data fusion problems at the point-cloud level [77]. During the SfM step, common features in the set of images were identified and matched, the internal camera parameters and relative orientation of the camera were estimated, and construction of the image network took place [16,79] to realize the sparse point cloud. Several papers have highlighted how a robust self-calibration procedure during the alignment stage is useful, due to the inherent instability of consumer cameras used in UAVs, and leads to better results with the inclusion of oblique imaging [21,57,80]. This first alignment allowed the removal of unwanted or outlying data (i.e., points that are clearly located off the surface or have anomalously large image residuals) and the deletion of the photos that the software did not align for different reasons. In the next step, georeferencing of the 3D sparse point cloud was carried out, testing three different solutions for the Belgian site:

• GCP solution: The traditional solution of the GCP coordinates (evaluating the level of GCP uncertainty before including these data to avoid adversely affecting data accuracy [58]) was employed to scale and georeference the SfM-derived point cloud. This result was obtained by using a seven-parameter linear similarity transformation, through locating and manually marking GCPs on at least two photographs. In this case, some of the GCPs (one-third) were used as control points (CPs) to provide an independent measure of accuracy.
• PPK solution: The camera positions from PPK data encoded in the images were used. Here, no information on GCPs was exploited.
• PPK + 1GCP solution: Besides the camera positions, as in PPK, the coordinates of one GCP were also used. This last configuration has been tested in several projects [57,65].

For the Italian study, instead, only the GCP solution was adopted to georeference the SfM point cloud, because the UAV used in this study could not exploit the other solutions. In light of this georeferencing data, the SfM step was improved [16,54] through a bundle adjustment step: a fundamental phase that refined the camera and tie-point locations (homologous points that link different images) and the camera calibration parameters of each image, through the bundle adjustment algorithm (least-squares network optimization [81]). This improved the values during the camera alignment step by incorporating georeferencing data and removing obvious outliers and incorrect matches from the sparse point cloud. Moreover, the optimization process was done through appropriate weighting of tie and control point image observations in the bundle adjustment to enhance a real error characterization [58]. This was followed by processing of high-density point clouds and orthomosaics, which involved the use of an MVS image-matching algorithm that increased the point density by several orders of magnitude [82]. This operates at the individual pixel scale to build dense clouds [83] and orthomosaics in agreement with the resolution of the photos.

The georeferenced SfM point clouds were imported in the CloudCompare software (Omnia Version 2.10.2; http://www.danielgm.net) to be filtered through different steps: manual filtering, a “distance filter procedure”, the cloth simulation filter (CSF) [84], and the “SOR filter tool” (Section v. in Figure 2). The SOR filter was used to remove outliers through the computation of the average distance of each point to its neighbors (it rejects the points that are farther than the average distance plus a defined number of times the standard deviation). The manual filter was used to delete unwanted objects in the point cloud (e.g., trees and shrubs on the lynchet risers and treads/flat part of terrace) in both SfM point clouds. The removal of vegetation in the SfM point cloud led to a gap of data at ground level (Figure 5a) that the TLS cloud filled. The CSF filter extracted the ground points only in vegetated lynchets of the Belgian site because this tool gives the best results on slopes that are not overly steep [85], such as the Belgian braided lynchets. For the rugged slope of the Italian terraces, a “distance filter procedure” was used instead to remove shrubs on the flat terrace areas (Figure 5b). This method consists of generating a mesh of the surface using the minimum elevation points within a defined value of a grid cell (the size is established by the user based on the point cloud density and the needed resolution) through the CloudCompare Rasterize tool. The mesh of the minimum elevations (considered the ground elevation) was used to calculate the cloud to mesh distance (C2M) [86] for the data points in the original point cloud. Then the C2M values (i.e., a scalar field in CloudCompare) were used to select points in the original cloud that were at a specific threshold distance from the ground (e.g., points recognized as vegetation). Therefore, the selected points were deleted, filtering/differentiating the terrace from the vegetation.

3.3.2. TLS Processing

Leica Cyclone 9.4 software was used to postprocess the TLS data. At first, each scan was manually filtered in order to remove unwanted objects, such as unnecessary elements of the landscape (Section iii. in Figure 2). The coregistration of TLS multiple scans was performed according to the following steps: (i) The software was automatically able to recognize the center of the black/white TLS targets (identified and labeled in the field) and paired to equivalent tie-points or overlapping features in multiple scans, thereby generating groups of prealigned scans (e.g., one group of scans for each of the three terraces surveyed in Belgian site). (ii) The prealignment was further refined with the built-in Cyclone variant of the well-known iterative closest point (ICP) [87,88] automatic algorithm. With this approach, the targets enabled us to add further constraints into the registration step and thus strengthen the alignment among the scans [89]. iii) The ICP optimization, implemented in Leica Cyclone, was then reapplied to merge all the grouped scans in a global georeferenced point cloud at the landscape scale. In order to georeference the TLS point cloud, the central points of the targets placed within the scanned scenes were assigned with the coordinates measured by GNSS in the field through a registration procedure based on automatic target recognition (Section iv. in Figure 2). The entire georeferenced TLS point clouds were processed employing the CloudCompare software through the same tools used to filter the point clouds for SfM (Figure 5d,e; Section v. in Figure 2).

3.4. Coregistration and Data Fusion

The filtered TLS and SfM point clouds must be coregistered to minimize residual inaccuracies of the georeferencing process. Among the different methods proposed to perform the coregistration task [90], we chose the ICP algorithm, which is the most common fine registration method for the alignment of point clouds published thus far [91,92,93]. This is implemented in several software packages, including CloudCompare. The ICP procedure permits the more uniform distribution of the residual registration error across the scans with respect to a simple pairwise approach [89], which can only be useful in particular cases. Before applying the ICP tool, it was important to check if the georeferencing process on the individual SfM and TLS datasets was effective and placed the different point clouds quite close together, in order to improve the efficiency of fine registration. If the georeferencing process is not optimal, it is preferable to implement a manual pairwise alignment through the identification of a set of corresponding points, mainly recognizable on overlapping areas. This coarse registration, based on at least four manually identified points, can also be performed in CloudCompare (with “Point Pair Picking” tool) to provide initial values for the following fine registration [39]. Then, the automatic ICP algorithm was applied and iteratively estimated pair points, computing a rigid body or a similarity transformation at each step (Section vi. in Figure 2). The SfM point cloud was used as the “reference” for the alignment process, while the TLS point clouds were chosen as “data to align”. This choice was made considering the georeferencing errors of each cloud and the surface covered by each technique. Initially, the ICP was executed using only subsets of overlapping point clouds (i.e., areas easily recognizable and objects detected by both techniques, e.g., natural and/or anthropogenic features). TLS and SfM targets were also used as homologous surfaces that were identifiable in the different point clouds by means of the associated RGB information. Finally, the estimated transformations were applied to the entirety of the original point clouds using the parameters derived by the ICP tool [71].

The coregistered TLS and SfM point clouds were used to extract the useful parts from each survey. Only terrace walls and lynchet risers were extracted from the TLS point clouds from the Italian and Belgian surveys, respectively, while only flat terrace zones and the wider landscape were derived from the SfM survey. The suitable HRT data from each survey were finally merged in CloudCompare (Figure 5c; Section vii. in Figure 2) generating data fusion point clouds: three for the Martelberg site (we tested three different SfM solutions in the georeferencing process, see Section 3.3.2) and one for the Soave study area.

3.5. DTM Generation

The data fusion point clouds were decimated in order to reduce processing constraints and extremely high data density. For this, we used the geostatistical Topography Point Cloud Analysis Toolkit (ToPCAT) [94], which has successfully been used to decimate point clouds in several studies (e.g., [92,95,96]). This tool allows intelligent decimation by decomposing the point cloud into a set of non-overlapping grid-cells and calculates statistics for the observations in each grid. Following the work by [97], the minimum elevation within each grid cell was taken to be the ground elevation, and a grid cell of 0.10 m was selected to regularize the data set. The point clouds obtained by ToPCAT were used to calculate a Triangular Irregular Network (TIN) that was then converted to raster through natural neighbors interpolator, obtaining four DTMs (Section viii. in Figure 2).

3.6. Data Analysis

In order to identify bias and estimate accuracy and precision of the obtained point clouds and DTMs, different analyses were carried out at different steps (Section ix. in Figure 2). Firstly, following [76], a bootstrapping resampling technique was implemented within Photoscan. This produced an estimation of the quality of the point cloud, considering the three different solutions used to produce the SfM models (see Section 3.3.2). In the case of the GCP solution, the bootstrapping resampling technique randomly selected one-third of the GCPs and used them as CPs (see Table 1) to provide an independent measure of uncertainty of each point (i.e., the residuals, or the difference between the real coordinates of this point and the modeled values). For the PPK + 1GCP solution, all GCPs minus one were randomly selected as CPs; in the case of the PPK solution, all the GCPs were used as CPs. This random selection was run iteratively 1000 times. Then, after all of the iterations, during which the bundle adjustment step was reset, the accuracy and precision were obtained for each point when used as GCP or CP. The mean of the residuals provides an indication of the accuracy of the registration process and the point cloud when the GCP and CP residuals are used, respectively. The standard deviation of the residuals yields an indication of the precision. This exercise provided an opportunity to check for potentially biased points.

To assess the effectiveness of the SfM–TLS coregistration process, the Multiscale Model to Model Cloud Comparison tool (M3C2) of CloudCompare [86] was used to calculate the cloud-to-cloud distance between the SfM and TLS point clouds located in overlapping areas, before and after the coregistration process [71]. The standard deviation of the distance between point clouds was used as an indicator of the measurement precision, while the mean absolute distance was considered as the accuracy of the point clouds [77].

Finally, the transformation of the point cloud into a continuous elevation surface (i.e., data interpolation), and the subsequent gridding, introduced several uncertainties or artifacts. They could influence the accuracy/quality of the surface produced, especially in complex scenes with vertical features, steep slopes, and rough surfaces [98,99]. The accuracy and precision evaluation of the geospatial products can be done by using independent checkpoints, reference surfaces, or length measurements [17]. Therefore, these errors, in particular in the vertical component, were evaluated through a statistical comparison between the Z values of CPs and the equivalent Z measures of DTMs. For the GCP solution, CPs were located in the most stable areas, following a homogeneous distribution. Instead, for the PPK + 1GCP solution, a single GCP was excluded, placed approximately in the middle of the area. For the PPK solution, all the CPs were considered. A number of CPs were selected for each DTM to calculate errors usually caused by filtering (for classifying into ground and off-terrain points) and interpolation (i.e., for filling gaps). Some studies [51,100,101] have highlighted that DTMs derived by laser scanning or digital photogrammetry seldom show a normal distribution of errors. Therefore, the metrics employed, such as the root-mean-square error (RMSE), standard deviation (SDE), and mean error (ME) had to assume that outliers existed and errors were not normally distributed. For this reason, following the approach of [100], the outliers were removed by applying a threshold selected from an initial calculation of the error measures. In this case, the threshold for removing outliers was selected as 2 times the RMSE. This means that an error was classified as an outlier if the absolute elevation difference between the CP and the corresponding point in the DTM was higher than 2 × RMSE. The value of 2 × RMSE was chosen considering the distribution of errors in the different surveys, keeping an approach that is as precautionary as possible. Moreover, to also test a robust accuracy assessment metric, preferably not influenced by outliers, we calculated the median (a robust estimator for a systematic shift of the DTM, less sensitive to outliers) and the normalized median absolute deviation (NMAD) [100,101]. They can be considered as estimates for the SDE that are more resilient to outliers in the dataset.

4. Results and Discussion

4.1. SfM Outputs

The assessment of errors in the SfM surveys is an important aspect of the workflow and was estimated through the point quality assessment of the GCPs and CPs in terms of precision, accuracy, and registration error.

Table 3. Registration errors, precision, and accuracy of all point clouds on ground control point (GCP) and CP residuals. The registration errors evaluated the GCPs and camera residuals after the point cloud georeferencing. The SDE and the mean absolute error (MAE) of the CP errors assessed the precision and accuracy of the point clouds, respectively. These errors are provided by the bootstrapping technique applied in Photoscan Pro software v 1.4.5 (see Section 3.6).

Solution Type	Accuracy				Precision			Registration
	MAE (m)			RMSE3D 1 (m)	SDE (m)			RMSE3D 1 (m)	Cameras 2 (m)
	X	Y	Z	CPs	X	Y	Z	GCPs
Italy
GCP	0.013	0.014	0.019	0.034	0.010	0.011	0.014	0.031	-
Belgium
GCP	0.014	0.015	0.052	0.070	0.013	0.011	0.025	0.068	-
PPK	0.016	0.017	0.060	0.072	0.012	0.012	0.026	-	0.088
PPK + 1GCP	0.016	0.016	0.058	0.071	0.012	0.012	0.025	0.072	0.088

¹ RMSE_3D (3D root-mean-square error) of GCPs and CPs computed along the x, y, and z directions, following [102]. ² Measures provided by Photoscan software.

summarizes all these error aspects for each of the SfM process carried out.

The bootstrap resampling technique allowed a more detailed assessment of errors based on the different SfM solutions in terms of the georeferencing process (see Section 3.3.2). The errors of the GCPs, CPs, and camera position (e.g., the distance between the input and estimated positions of the camera in Photoscan), which determine overall SfM survey accuracy, were all in the order of magnitude of centimeters, adequate for investigating topographic features of the terrace complex. As for the test of the different solutions for georeferencing the SfM point cloud in the Belgian site, Table 3 shows that the best results were those of the GCP solution. However, the values for GCP solution were very close to the estimated errors of the PPK solution, both in terms of accuracy, precision, and RMSE_3D. These results highlighted that direct georeferencing, with very accurate positioning, was able to match the conventional GCP method. This was not possible in the past, when the direct georeferencing method could achieve only decimeter to meter accuracies, as indicated in many studies [103,104,105,106]. This was caused by the low accuracy of the navigation-grade onboard GPS units used for most current consumer UAVs [49,56,103]. Recently the DG methods have gradually improved and can deliver centimeter-level accuracy [51,107], exploiting multifrequency GNSS receivers (with devices light enough to be carried in <25 kg UAV platforms) and kinematic data postprocessing [64,108,109]. Our empirical result agrees with the simulations of [110], where georeferencing with GCPs or DG was found to yield comparable ground point vertical precisions. In addition, [60,63,64] provided high-quality camera position measurements and enabled direct georeferencing to centimeter-level accuracy through the PPK solution. Moreover, using the PPK + 1GCP solution, the performance was slightly enhanced in terms of accuracy (in particular for the Z component) and RMSE_3D (Table 3). These results were in line with previous studies [21,106,111,112] which highlighted that direct georeferencing allowed the acquisition of robust centimetric HRT data, similar to the GCP solution [57,107,113]. This is the case if biases in the elevation component are controlled using at least one GCP to strengthen the self-calibration. Indeed, Reference [65] showed that the application of PPK in direct georeferencing could provide the same centimeter-level accuracy and precision as the GCP traditional approach. Nevertheless, some UAV surveys were characterized by a vertical shift that could be mitigated using a single GCP. Using a single GCP provided a robust way to detect perturbations of the GNSS signal, given that it was difficult to assess the quality of the PPK solution without independent observation. The addition of at least one GCP to the camera stations, as specified by [59], successfully removes most of the bias (especially in elevations); though less precise than GCPs, GNSS observations of the camera stations spread the control over the block fairly homogeneously. The addition of limited but high-accuracy GCPs has had a dominant effect on the georeferencing quality and a sparse network of GCPs could be supported by a dense network of camera positions in a DG workflow, as stated by [56]. The location of the additional point might be worth investigating, since it is well known that UAV-SfM accuracy depends on the location and number of GCPs introduced in the optimization step [74]. However, the used bootstrapping technique allowed the assessment of the different possible localizations of the added GCP, among those available, choosing a different GCP at each interaction. This allows us to verify the effectiveness of adding a GCP, regardless of its location within the study area, in a more thorough and accurate way. Additionally, the bootstrapping technique also permits us to quickly check the additional benefit of using more (e.g., three or one-third of the total GCPs) well-distributed points, when the conditions on site allowed for it. These tests showed that with the addition of three or one-third of the GCPs to the PPK solution, there was no improvement in results, and the accuracy and precision values remained stationary and very similar to the PPK + 1GCP solution. This result was confirmed by [106], where the use of many more GCPs compared to just one indicated no further improvement of the final 3D model. From the above remarks, it is clear that few GCPs are still necessary even in direct georeferencing, but a clearer future picture of this technique is emerging; the PPK + 1GCP solution can certainly help to survey inaccessible areas (e.g., western side of the Belgian site; see Figure 4a), where GCP placement is impossible. Indeed, GCP location is a time-consuming process involving a considerable manual effort in field positioning, especially when the survey area includes inaccessible or hazardous terrains [64]. Therefore, one of the best solutions, as in the Belgian site, is a combination of direct and indirect georeferencing, i.e., a technique called integrated sensor orientation (ISO) or a special case of ISO without IMU data—GNSS-assisted aerial triangulation (GNSS-AT) [57,59,109]. These techniques are very useful in cases of very steep vegetated areas, motion blur and image noise, or anywhere the SfM-based approaches may not succeed. Indeed, a previous study on GNSS-AT by [113] underlined how this approach allowed a reduction in the minimum the number of GCPs required, decreasing the substantial inconvenience and cost of GCP establishment, especially in areas with difficult access.

Another aspect to be considered in order to improve the results in slope UAV surveying is the use of a tilted camera. UAV oblique images increased the robustness of the geometrical model, also providing a possible strategy for reducing the total number of GCPs adopted over a given area [50]. This helped to enhance the elevation accuracy in the Belgian UAV survey, especially in the inaccessible areas. Moreover, the accuracy of the obtained SfM results highlighted that the acquisition of images from two different observation directions (i.e., oblique and nadir) and platforms (i.e., UAV and ground-based) led to denser point clouds, an optimal camera network geometry (i.e., great image overlap and high angle of convergence), and fewer deformation errors or large-area distortions. Several pieces of research confirmed that the addition of oblique photographs in a UAV survey considerably strengthens the network geometry [37,58,68,114,115]. This also reduced error in estimated camera parameters and the likelihood of detectable systematic DTM error, such as the “doming effect”. Moreover, the direct fusion of ground-based and UAV photos (rarely processed together) in the SfM process in Photoscan avoided subsequent merging problems at point clouds level [77]. Many authors [17,54,116] recently highlighted how the development of a framework or workflow that allows data to be correctly processed and analyzed in accordance with the objective of the work is a critical part of SfM. Moreover, moving beyond such a ‘black-box’ SfM approach, in this work, an understanding of the SfM sources of error allowed us to identify the best strategy to operate in a topographically complex environment, such as vegetated terraces.

4.2. Coregistration Process

The first alignment process was applied to coregister the multistation TLS data in Cyclone software (see Section 3.3.2). The limited number of TLS targets over the study area caused a minor tilting of the point clouds. Therefore, the ICP algorithm implemented in Cyclone was applied to enhance the registration of the point clouds, as in [117]. The Cyclone output errors for the final georeferencing process were RMS values of 0.053 and 0.079 m for Italian and Belgian sites, respectively. This result, together with those shown in Table 3, demonstrate how the georeferencing process was accurate enough to place the TLS and SfM point clouds close together and to apply the ICP algorithm without a prealignment step (see Section 3.4) using coarse registration. The ICP output errors in CloudCompare for the SfM–TLS coregistration process were RMS values of 0.045 and 0.068 m for Italian and Belgian sites, respectively. In other case studies where the point clouds were not already close enough, the coarse registration was used only as a first prealignment step to place the TLS and SfM closer together, improving the quality of the following ICP registration [39]. This procedure was also applied in [32], which highlighted how the result of coarse registration was generally poor but necessary for a successful ICP if the point clouds extracted from the SfM and TLS presented differences in terms of sensing resolution, coverage, and accuracy.

To carry out the TLS and SfM alignment process, it was necessary to identify overlapping natural and/or anthropogenic surfaces easily recognizable in both corresponding point clouds. The identification of the corresponding surface in the alignment process may be prone to errors [118], especially in rugged and vegetated environments [71], due to (i) the difficulty for the operator to select the same areas in clouds generated from different methods (e.g., different point density and detected objects); (ii) the uneven point distribution on different data sets, which often do not perfectly cover the same areas [93]; and (iii) the roughness of surfaces, which reduces the accuracy of cloud matching techniques (e.g., [86]). Therefore, the use of common targets and the effective distribution of targets within the study area, so that they were visible from both the TLS and SfM surveys, helped to identify homologous surfaces for the alignment. This was also underlined in [19], which confirmed how a set of common references clearly identified in each dataset was indispensable for fusing different data, due to the facts that every data source has its own characteristics and the identification of common natural references is sometimes a hard task. Therefore, the survey planning was an essential step in which the positioning of targets had to be carefully considered so that they could be easily recognized in each dataset to enable the merging of the data. The choice to start the ICP process from the overlapping stable areas and then extend the coregistration process to the whole area certainly improved the quality of the alignment, as demonstrated in several studies [71,91]. Other previous studies [35,36,37,39,90] in which the ICP method was used to coregister TLS and SfM data showed good results.

The effectiveness of the coregistration method was numerically confirmed by the results provided by the cloud-to-cloud M3C2 distance tool applied in overlapping stable areas. As shown in

Table 4. M3C2 cloud-to-cloud distance values in the overlapping stable zones, obtained using the coregistration method explained in Section 3.4.

M3C2 Distance	Non-Coregistered Point Clouds			Coregistered Point Clouds
	Accuracy		Precision	Accuracy		Precision
Solution Type	MAE (m)	ME (m)	SDE (m)	MAE (m)	ME (m)	SDE (m)
Italy
GCP	0.169	0.169	0.006	0.048	−0.027	0.069
Belgium
GCP	0.209	0.218	0.120	0.080	0.018	0.102
PPK	0.246	−0.246	0.155	0.106	0.019	0.132
PPK + 1GCP	0.219	−0.208	0.177	0.069	0.004	0.089

, the accuracy and precision of the alignment in overlapping stable areas considerably increased after the coregistration in terms of MAE, ME, and SDE for each survey and SfM solution type. The SfM survey that shows the best results in terms of M3C2 distance from TLS point cloud were obtained for the GCP solution used at the Belgian site. However, the M3C2 values for the PPK + 1GCP solution were very close to GCP solution also in this context, further highlighting the similarity between these two solutions.

4.3. Data Fusion

Data processing was carried out for four SfM–TLS point clouds, one for Soave site (104,614,919 points, with a density of 2990 points/m²) and three for the Belgian study area (mean of 298,739,953 points, with a mean density of 1660 points/m²). The obtained data fusion clouds (Figure 6a,b) show that an integrated approach of SfM images and TLS acquisition was needed to survey the complex topography of these terraces. Indeed, TLS data provided a more accurate representation of subvertical surfaces covered by vegetation (e.g., the vertical walls or risers of terrace) where the SfM surveys were limited. Figure 5b shows a clear example of the lack of data in the Soave SfM point cloud due to the presence of shrubs in front of the terrace walls that, hidden by vegetation, cannot be detected by the SfM technique. The removal of vegetation also led to a lack of data at ground level (Figure 5a). This data gap was filled by TLS surveys that allowed us to obtain more topographic information relating to the terrace surface (e.g., wall and risers). Indeed, the upper part of Figure 6a highlights how, using only SfM data, the terrace risers were reconstructed by interpolation, creating an unrealistic surface with no roughness. Instead, the addition of the ground-based TLS cloud allowed a 3D model reconstruction (e.g., lower part of Figure 6a) with the real roughness of the terrace risers covered by vegetation. This works in areas where the cultivated vegetation does not have a dense canopy which can be penetrated by a simple TLS pulse. When the vegetation canopy is very dense, it is necessary to place the TLS station very close to the surveying object to better detect ground data. Otherwise, a full-waveform TLS is required. At the same time, ground-based TLS was unreliable over large areas due to a very oblique perspective with a very low incidence angle. Therefore, UAV SfM data integrated the survey, covering large areas on a relatively flat zone in small amounts of time, and allowed us to create an accurate HRT survey of the whole landscape (e.g., in Figure 6b). Moreover, terrestrial images provided a more accurate survey of particular objects such as wall collapses at a detailed scale, permitting the study of the ongoing processes of the terrace systems.

The different perspectives and resolutions greatly increased the point cloud density, the individual point precision, the robustness of topographic mapping, and the high-resolution detail. The union of ground-based (i.e., TLS data and terrestrial images) and aerial data allowed us to overcome many of the disadvantages of a single method, such as the line-of-sight obscuration that occurred due to vegetation and other complex objects. The effectiveness of similar data fusion methods in a complex topographic context has been demonstrated in several pieces of research [19,33] but not in agricultural terraces, which have particular geometrical properties.

In regions of complex topography with multiple local horizons, restricted lines of sight significantly hinder the use of tripod-based instruments by requiring multiple setups and taking more time to achieve full coverage of the area. Therefore, a further step forward in data fusion could be to test the use of the hand-held mobile laser scanning (MLS) technique. This offers particular promise for topographic surveys of complex environments, combining the reliability of laser techniques with the flexibility of on-foot surveying and delivering a data density typical of scanning systems [119]. This technology has emerged and considerably improved over the years [120]; it can certainly be exploited more in combination with the SfM technique for environmental monitoring. Another possible improvement in terrace system surveys could be the use of small and fully integrated laser scanning sensors for UAVs. This solution would allow the use of laser technology at a larger scale and from all possible points of view, making SfM survey no longer necessary. However, the high cost of such systems makes them currently inaccessible for most research projects. For these reasons, the use of such technologies in environmental monitoring is still greatly limited [121].

4.4. DTM Error Assesment

The SfM–TLS data fusion allowed the generation of four DTMs at 0.10 m (Figure 7), three (i.e., one for each SfM solution; see Section 3.3.2) for the whole study area in Belgium (Figure 7a) and one for the Italian site (Figure 7b). Given that the main final survey products of the HRT surveys were DTMs, the assessment of their error was necessary to understand the quality of the data fusion process. This was also useful to evaluate the different accuracies of the geospatial products based on the methods used for the SfM georeferencing. DTM error assessment was done by using independent CPs, frequently used in the literature to assess the difference in elevation between experimental evidence from real data and DTMs [16,68,122,123].

Table 5. Digital terrain model (DTM) error statistics: the difference between the Z values of the CPs and the Z measures of the same points on the DTMs.

Solution Type	MAE (m)	ME (m)	SDE (m)	RMSE (m)	Median 1 (m)	NMAD 2 (m)
Italy
GCP	0.024	−0.017	0.026	0.040	−0.022	0.035
Belgium
GCP	0.063	−0.033	0.073	0.075	−0.038	0.069
PPK	0.095	−0.095	0.035	0.101	−0.079	0.099
PPK+1GCP	0.063	−0.062	0.034	0.080	−0.054	0.078

¹ 50% quantile. ² Normalized median absolute deviation: proportional to the median of the absolute differences between errors and the median error [100].

summarizes all the statistics regarding the residuals between DTMs and validation points.

The statistics reported in Table 5 show how the interpolation process slightly increased the errors in DTM data compared to the errors at the point cloud level. It is well known that gridding processes may induce a loss of resolution and increased errors with respect to the original data [55,98,99,124,125], especially in complex topography [77]. However, the quality of the DTMs obtained was certainly high (centimeter-level), considering the values in Table 5, which demonstrated how the data fusion process proved effective in generating HRT data. The wide range of calculated metrics allowed a more robust error estimation, especially after the outlier removal that made the statistics more reliable [100]. Moreover, the use of metrics such as median and NDAM further reinforced the data obtained from other statistics as being very similar. Several authors have stated the importance of using robust accuracy assessment methodology and metrics, preferably not influenced by outliers or by a skew in the distribution of the errors [51,101,126]. As in the previous analysis for the Belgian site, more minor errors were found for the DTM generated from the GCP solution than those generated by the PPK and PPK + 1GCP solutions in terms of accuracy and precision, although they were very close. The DTM produced with the PPK method presented a lower accuracy than that resulting from the PPK + 1GCP solution, which showed values very close to the traditional procedure of the GCP solution. Similar outcomes were obtained by [59], where eight 10-cm-resolution digital surface models (DSMs) generated a mean RMSE in elevation very close (i.e., approximately 5.4 cm), using either GCPs or all of the camera stations (i.e., direct georeferencing) and one GCP. This underlines once again that the combination of direct and indirect georeferencing can lead to the generation of very accurate DTMs with a quality similar to that of DTMs obtained with traditional methods.

The DTM validations using CPs demonstrated the high quality of the DTMs obtained from the data fusion, and Figure 7a,b showed the very high resolution of these DTMs. The results highlighted the effectiveness of the developed workflow, which had proven to be valid and fundamental in two different sites with diverse survey issues. The SfM–TLS data fusion allowed us to obtain accurate DTMs with few gaps, despite the presence of considerable vegetation, and to reconstruct the complex morphology of these terraced areas both at the landscape scale and at a detailed scale. From these high-resolution DTMs, it is possible to extract 3D metric information about the terrace complex which can be very useful for geomorphological analyses. In addition, very detailed DTMs can be used as benchmarks for numerical and physical modeling or simulations of erosion process and soil formation in land degradation analysis of agricultural environments. Moreover, in a multidisciplinary context, such as that of the TerrACE project (see Section 3.1), high-resolution topographic data can be integrated with subsurface information (e.g., chronostratigraphic models of geoarchaeological trenches). This will allow the dating of terrace systems, the determination of their form and construction, and the understanding of their original and later purposes and use. Therefore, we will be able to understand how human societies have been reshaping the geomorphology of landscapes for thousands of years [127]. Data integration from different studies can provide helpful information to find the best ways of preserving agricultural terraces.

5. Conclusions

This research has highlighted how the SfM–TLS data fusion can be used to create highly accurate DTMs (centimeter-level) in the context of complex topography and vegetation cover, such as in the case of agricultural terrace systems. The data fusion of these HRT techniques made it possible to overcome the weaknesses, while still exploiting the advantages, of each surveying method. This SfM–TLS merging method, applied for the first time in vegetated abandoned terrace systems, has proved to be an effective solution for reconstructing these agricultural structures in a more accurate, complete, and realistic way, even with vegetation cover. The presence of unreachable areas, where it was impossible to place the GCPs, forced us to test alternative solutions for the SfM survey such as direct georeferencing, exploiting the benefits of the PPK method. The results showed how the addition of at least one GCP to the camera stations reduced direct georeferencing biases in the elevation component and ensured the quality of GCP solution (both at the point-cloud and the DTM level). Therefore, in inaccessible or hazardous terrains typical of rugged vegetated areas, a combination of direct and indirect georeferencing can be a useful solution, allowing the substantial inconvenience and cost of GCP placement to be reduced. Without the use of these combinations of methods and different acquisition platforms, it was impossible to achieve accurate DTMs that reflected the real surface roughness without a loss of data. To obtain an accurate data fusion, it was essential to utilize a specific workflow that considered all data merging issues, the complex topography, and the landcover condition problems in terrace systems. The most important phase in this workflow was the coregistration process, which played a key role for effective data fusion of the point clouds; this process must be considered during the planning of the survey in the field.

The high-resolution DTMs realized through SfM–TLS data fusion are the starting point for generating valuable metric information about ancient agricultural terrace complexes. These data, when integrated with subsurface soil information, can be extremely useful in studying agricultural terrace environments through a multidisciplinary approach. Information merged from diverse surveying disciplines can supply valuable and comprehensive landscape information for environmental authorities, facilitating the improvement of guidelines for the management and preservation of agricultural terrace heritage sites.