Integrated Laser Scanner Techniques to Produce High-Resolution DTM of Vegetated Territory

: The paper presents the ﬁrst part of a research project concerning the creation of 3D terrain models useful to understand landslide movements. Thus, it illustrates the creation process of a multi-source high-resolution Digital Terrain Model (DTM) in very dense vegetated areas obtained by integrating 3D data coming from three sources, starting from long and medium-range Terrestrial Laser Scanner up to a Backpack Indoor Mobile Mapping System. The point clouds are georeferenced by means of RKT GNSS points and automatically ﬁltered using a Cloth Simulation Filter algorithm to separate points belonging to the ground. Those points are interpolated to produce the DTMs which are then mosaicked to obtain a unique multi-resolution DTM that plays a crucial role in the detection and identiﬁcation of speciﬁc geological features otherwise visible. Standard deviation of residuals of the DTM varies from 0.105 m to 0.176 m for Z coordinate, from 0.065 m to 0.300 m for X and from 0.034 m to 0.175 m for Y. The area under investigation belongs to the Municipality of Piuro (SO) and includes both the town and surrounding valley. It was affected by a dramatic landslide in 1618 that destroyed the entire village. Numerous challenges have been faced, caused both by the characteristics of the area and the processed data. The complexity of the case study turns out to be an excellent test bench for the employed technologies, providing the opportunity to precisely identify the needed direction to obtain future promising results.


The Need of a High Resolution DTM
Digital mapping of the territory and its related phenomena has become a fundamental requirement in order to keep track of all the alterations in the surveyed environment. By analyzing the surface morphology of the earth, it is possible to reconstruct the processes that shaped our planet as it is today. In the last few decades, the enormous advances in the field of Remote Sensing (RS) techniques have allowed for collecting an impressive quantity of data representing the earth's surface and its changes over time by means of imagery [1] and Laser Scanner (LS) [2,3]. This becomes crucial to detect, according to the employed instruments and achieved resolution, even the most localized event, to monitor environmental changes in time, thus building a constantly updated historical archive and establishing successful prevention and mitigation risk strategies [4][5][6][7]. Therefore, consideration should be given to the prediction, detection, and monitoring of those events most responsible for the hydrogeological risk [8,9], such as landslides [1,4,10,11], floods [6,12,13], active tectonics [14][15][16], and bank erosion phenomena [1,17].
In this framework, a complete representation of the area under investigation by means of Digital Terrain Models (DTMs) allows for simulating the events and therefore and the surrounding valley (Figure 1), covering an area of about 6 km 2 . Mera river runs from east to west and conveys the waters of numerous streams present in the valley. The typical alpine orography characterizes the area: both of the sides, the northern one with an average slope of 28 • and the southern one with 34 • average slope, are covered with very dense vegetation that includes traces of geomorphological events of glacial origin. This is due to the shaping action of the huge ice flow descending along the entire Chiavenna valley during the last glaciation.
The typical alpine orography characterizes the area: both of the sides, the northern one with an average slope of 28° and the southern one with 34° average slope, are covered with very dense vegetation that includes traces of geomorphological events of glacial origin. This is due to the shaping action of the huge ice flow descending along the entire Chiavenna valley during the last glaciation.
The tree and shrub species forming the forest include various native entities with a prevalence of chestnut, hornbeam, ash, sycamore, linden, oak, cherry, and mountain ash.
Moreover, large and still visible landslide sediments from to the 1618 landslide, which destroyed the town and surrounding areas, are present throughout the whole region of interest [37]. The dramatic event had a very strong impact on the evolution of the valley, both in terms of geological characterization and loss of human lives; historical documentation reports the entire destruction of Piuro and Scilano villages, causing the death of 1000-2000 inhabitants. It is considered one of the most catastrophic events ever recorded in the region [38].  The tree and shrub species forming the forest include various native entities with a prevalence of chestnut, hornbeam, ash, sycamore, linden, oak, cherry, and mountain ash.
Moreover, large and still visible landslide sediments from to the 1618 landslide, which destroyed the town and surrounding areas, are present throughout the whole region of interest [37]. The dramatic event had a very strong impact on the evolution of the valley, both in terms of geological characterization and loss of human lives; historical documentation reports the entire destruction of Piuro and Scilano villages, causing the death of 1000-2000 inhabitants. It is considered one of the most catastrophic events ever recorded in the region [38].
The research is part of the Interreg project A.M.AL.PI. 18 (Alpi in movimento, Movimento nelle ALpi. Piuro 1618), aimed at realizing a best practice to study historical behavior of dangerous hillsides applicable to all the case studies identified along the Alpine range, Remote Sens. 2021, 13, 2504 4 of 23 that extends from the Gottardo Pass to the Maloja Pass, both on Italian and Swiss territory. The ambitious goal is to increase the attractiveness of the territories of the Bregaglia-Chiavenna valley-Moesa-Ticino axis by focusing on an innovative fruition of natural and cultural resources. The project aims to demonstrate the historical specificity of the involved valleys, symbol of the richness of the Alpine routes, and to investigate the importance of geological knowledge in the prevention of natural disasters. Therefore, the research project consists of the creation of high-resolution 3D terrain models, useful to understand landslides movements of the past, integrating multi-scale and multi-techniques survey methods to produce a sort of reverse-engineering of the territory through a multidisciplinary approach involving the other A.M.AL.PI. 18 project partners. Landslides have historically occurred in the concerned territory and the hydrogeological risk is still significant today [39,40]. For this reason, high-resolution digital mapping helps improve its prevention and mitigation.
The paper presents the first part of the research that is the creation of a high-resolution DTM obtained by integrating 3D data (Point Clouds-PC) coming from several sources, starting from long and medium-range Terrestrial Laser Scanners (TLS) up to a Backpack Indoor Mobile Mapping System (IMMS). Since some specific sites are considered strategic for a proper and complete interpretation of the landslide event in the area under investigation, they are surveyed with LS techniques at a very high resolution (5 mm to 40 mm), in accordance with the other partners involved in the project.

Data Integration Framework
The integration of data characterized by different resolutions and coming from several sources is something already faced in literature. The chance of creating a multi-resolution model integrating data coming from different kinds of sensors presents the possibility to obtain a constantly updated 3D digital representation that deals with the appropriate level of detail wherever necessary. Considering the limits of a single-sensor approach and the fact that currently a sensor that works effectively in all conditions does not exist, this multi-source approach allows for creating a hybrid model that takes advantage of the potentialities of the employed instruments. Data integration is particularly investigated in the field of surveys for the documentation of complex architectures, where it allows for achieving great adaptability when taking into account 3D metric acquisition in different scenarios and where it constitutes a fundamental instrument for the analysis and interpretation of historical sites [41][42][43].
Things are different when referring to multi-source and multi-resolution surveys in the environmental field, particularly those aimed at providing DTMs for different purposes, from the modeling of gravitational phenomena, rockfall, avalanches, to surface run-off. The Heli-DEM project [44], for example, integrates data at the border territory between Italy and Switzerland in order to create a unique DTM for the alpine and subalpine area, which is useful to analyze the scenario when an event, such as a landslide, happens. Together with the aforementioned multi-source multi-resolution DTM, noteworthy data integrations were conducted at the national or even larger scale, combining data coming from preexisting DTMs with different resolutions and accuracies [45][46][47], vector contour lines from maps, break lines, geodetic points, and photogrammetric data [48].
The investigations and processing cited above include available products with a resolution from 1 m up to 100 m. Given the goal of the current project, however, it is not possible to rely only on the existing datasets, due to their inadequate resolution. The reference is to the 5 m and 20 m resolution DTMs made available by the institutions, and to the 1:10,000 contour lines. Traces left by the events that shaped the area as it appears today would not be visible, like metric and sub-metric morphological highs, erosion of riverbanks due to ancient river course, deposits of finer material, and small fractures with different dips along the sides. Consequently, the possibility of reconstructing precisely the history of the valley would be lost. The need of a big scale high-resolution DTM arises spontaneously. Therefore, the present research wants to play a key role in the multi-resolution 3D model framework, and it is the opportunity to establish some useful workflows to speed up acquisition and processing time in different scenarios, especially those regarding the monitoring of hydrogeological instability. The resulting DTM is 0.50 m in resolution and accuracy and reaches 0.10 m resolution in specific sites identified with the other A.M.AL.PI. 18 partners, allowing for recognizing even the single rocks dropped from the slopes and deposited on the valley floor. A processed high-resolution DTM such as this one constitutes the starting point from which to assume the evolution of the landform, recognizing the thickness and disposition of the different debris deposits. This will allow for making reliable predictions about the composition of the subsoil and to produce detailed 3D geological sections in order to be able to reconstruct the territory of Piuro back in time. By progressively removing the recognizable volumes from the DTM, a series of derived DTMs corresponding to specific historical periods will be produced, up to hypothesizing how the territory was before the 1618 landslide.

Materials and Methods
In this study, several surveys have been conducted, using the appropriate instrument depending on the characteristics of the survey and of the instrument itself. In the following section, data acquisition and processing procedures are presented and grouped on the basis of the employed instrumentation, i.e., long-range LS, medium-range LS and IMMS. A flowchart summarizing the whole adopted methodology and a table listing pros and cons of the employed instrumentation are proposed at the end of the chapter.

Long-Range TLS RIEGL VZ-4000
Data constituting the most general DTM were collected with long-range TLS RIEGL Vz-4000 between March 2017 and October 2020. The instrument offers an extremely long range of more than 4000 m and a wide field of view of 60 • vertical and 360 • horizontal. These peculiarities are ideal for the acquisition of both the slopes and the valley floor of the area under investigation, since it extends for about 2.5 km both in the E-O and N-S directions. The scanner is mounted on a tripod, leveled and connected to a laptop and an external battery. The characteristics of the scans, which are the desired resolution and the area range to be surveyed, are set up in the RiSCAN PRO software (RIEGL, Horn, Austria) at the moment of the acquisition. In addition, a built-in calibrated 5-Mpixel camera capturing images deflected by the laser mirror enables the coverage of the entire field of view with an appropriate number of high-resolution images automatically stitched together to create a high-resolution panorama image. This last one enables the creation of RGB colored PCs. In addition, the possibility of recording multiple echoes becomes crucial in order to split ground/non-ground points.
A total of 7 scans ( Figure 2) were acquired in the mentioned period, stationing between 450 and 900 masl-the valley floor is about 400 masl, making sure to cover the area as uniformly as possible. Despite surveying the area in times of the year when vegetation is sparser to have a better chance of getting the laser beams to the ground, scan n. 7 is characterized by dense vegetation in the area of the alluvial fan formed by the Valle Drana stream, thus not allowing an accurate ground detection. Furthermore, the high incidence angle of the laser beam did not help reach the ground in an effective way. The final total number of acquired points is 229,179,174 and the 7 scans' principal characteristics are summarized in Table 1.  Data Processing The 7 scans were registered using 21 homologous points detected throughout the whole area, obtaining an RMS error of 0.12 m. Subsequently, they have been georeferenced using 13 points coming from an RTK GNSS survey. This procedure allowed for georeferencing the scans all at once because not all of the RTK GNSS points were visible from the 7 different positions. The RMS error related to the RTK GNSS points is 0.36 m. The obtained errors are consistent with the resolution of the PCs and the 1:2000 scale for which the final product is designed.
In order to separate ground and non-ground points, the acquired territory was divided into three homogeneous areas: the valley floor, the north side and south side ( Figure  3A), due to common orographic characteristics. The surface-based Cloth Simulation Filter (CSF) algorithm [49], implemented in CloudCompare software [50], was then applied to the single + first echoes of each of the identified areas, each one tilted by ±30°, ±45° and  The 7 scans were registered using 21 homologous points detected throughout the whole area, obtaining an RMS error of 0.12 m. Subsequently, they have been georeferenced using 13 points coming from an RTK GNSS survey. This procedure allowed for georeferencing the scans all at once because not all of the RTK GNSS points were visible from the 7 different positions. The RMS error related to the RTK GNSS points is 0.36 m. The obtained errors are consistent with the resolution of the PCs and the 1:2000 scale for which the final product is designed.
In order to separate ground and non-ground points, the acquired territory was divided into three homogeneous areas: the valley floor, the north side and south side ( Figure 3A), due to common orographic characteristics. The surface-based Cloth Simulation Filter (CSF) algorithm [49], implemented in CloudCompare software [50], was then applied to the single + first echoes of each of the identified areas, each one tilted by ±30 • , ±45 • and ±60 • to preserve points also in the steepest slopes. To tilt the PCs, an interpolating plane was first fitted to the involved area and then made horizontal. The specified rotation matrices were then applied to rotate the PC by the desired angle. The choice of the returns to be selected in order to extract ground points is a matter of debate [51]. The last return is considered the most reliable since it has the highest chance to reach the terrain in dense vegetated areas. However, considering the characteristics of the area under study, which is dominated by dense vegetation on both the sides but made also of large crops, fields and private gardens, even the first echo was used to extract the terrain, in addition to the last echo.
The CSF algorithm plugin is simple to use thanks to the limited number of integer and Boolean parameters to be set. It is based on the cloth simulation 3D computer graphics algorithm [52] and consists of inverting the PC and applying a rigid cloth to it. By analyzing the interactions between the cloth nodes and the corresponding points of the PC, it is possible to define the final shape of the cloth. The plug-in allows for setting the type of terrain, which can be a steep slope, relief, or flat, which then determines the rigidness of the cloth. When dealing with very steep slopes, the internal constraints among particles do not allow a proper reconstruction of the ground. For this reason, there is the possibility to select an extra option to solve the problem. The second step involves the setting of cloth parameters. The resolution refers to the grid size of the cloth used to cover the terrain. The higher the resolution, the sparser the selection of ground points. The classification threshold defines the maximum distance inside which the points of the PC are then taken into consideration to be eventually added to the ground points set. Finally, it is possible to set the maximum number of iterations of the process.
The CSF algorithm was then applied to the 21 obtained PCs-the horizontal ones were also used to extract ground points. After several tests, the final ground points were obtained selecting a 0.50 m cloth resolution and a 0.10 m classification threshold. The obtained points (more than 42 millions) were interpolated in ArcGIS Pro (Esri, Redlands, California) [53], setting a 0.50 m raster cell and using the Inverse Distance Weighting (IDW) algorithm ( Figure 3C), which is considered to be a good interpolation technique thanks to a reasonable computational time and the non-addition of artefacts in the final product [29]. Considering the location of the area, the WGS84 UTM 32N projection was assigned to the output raster.
Remote Sens. 2021, 13, 2504 7 of 2 ±60° to preserve points also in the steepest slopes. To tilt the PCs, an interpolating plan was first fitted to the involved area and then made horizontal. The specified rotation ma trices were then applied to rotate the PC by the desired angle. The choice of the returns t be selected in order to extract ground points is a matter of debate [51]. The last return considered the most reliable since it has the highest chance to reach the terrain in dens vegetated areas. However, considering the characteristics of the area under study, whic is dominated by dense vegetation on both the sides but made also of large crops, field and private gardens, even the first echo was used to extract the terrain, in addition to th last echo. The CSF algorithm plugin is simple to use thanks to the limited number of intege and Boolean parameters to be set. It is based on the cloth simulation 3D computer graphic algorithm [52] and consists of inverting the PC and applying a rigid cloth to it. By analyz ing the interactions between the cloth nodes and the corresponding points of the PC, it possible to define the final shape of the cloth. The plug-in allows for setting the type o terrain, which can be a steep slope, relief, or flat, which then determines the rigidness o the cloth. When dealing with very steep slopes, the internal constraints among particle do not allow a proper reconstruction of the ground. For this reason, there is the possibilit to select an extra option to solve the problem. The second step involves the setting of clot parameters. The resolution refers to the grid size of the cloth used to cover the terrain. Th higher the resolution, the sparser the selection of ground points. The classification thresh old defines the maximum distance inside which the points of the PC are then taken int consideration to be eventually added to the ground points set. Finally, it is possible to se the maximum number of iterations of the process.
The CSF algorithm was then applied to the 21 obtained PCs-the horizontal one were also used to extract ground points. After several tests, the final ground points wer obtained selecting a 0.50 m cloth resolution and a 0.10 m classification threshold. The ob tained points (more than 42 millions) were interpolated in ArcGIS Pro (Esri, Redland California) [53], setting a 0.50 m raster cell and using the Inverse Distance Weightin (IDW) algorithm ( Figure 3C), which is considered to be a good interpolation techniqu thanks to a reasonable computational time and the non-addition of artefacts in the fina product [29]. Considering the location of the area, the WGS84 UTM 32N projection wa assigned to the output raster.

Medium-Range TLS Leica RTC360
TLS Leica RTC360 is used to collect data in the "Belfort" archaeological site (Figur 4), an area considered strategic for a proper and complete interpretation of the 1618 land slide. Located in the eastern side of the region of interest (highlighted in yellow in Figur 1), the site extends for about 1 ha and is one of the rare remaining and escaped example of the antique Piuro village before being destroyed by the 1618 landslide. The geometr is very challenging, considering the presence of subways, caves, very high vertical ston walls, and dense vegetation all around.
Two surveys were performed between September and November 2019, reconstruc ing the area by means of 122 scans. The first 104 were collected during the first surve while the last 18, which were necessary to complete the far eastern section, were collecte during the second one.
Thanks to the scanner's 360° horizontal and 300° vertical field of view, even the mo closed and narrow environments have been completely acquired without any problem The acquisition is very fast because the scanner is self-leveling and every scan only r

Medium-Range TLS Leica RTC360
TLS Leica RTC360 is used to collect data in the "Belfort" archaeological site (Figure 4), an area considered strategic for a proper and complete interpretation of the 1618 landslide. Located in the eastern side of the region of interest (highlighted in yellow in Figure 1), the site extends for about 1 ha and is one of the rare remaining and escaped examples of the antique Piuro village before being destroyed by the 1618 landslide. The geometry is very challenging, considering the presence of subways, caves, very high vertical stone walls, and dense vegetation all around. quires less than 3 minutes, setting the highest resolution of 3 mm at 10 m. The simultaneous acquisition of full spherical HDR images through the incorporated 36 Mpixel camera allows for obtaining the RGB information to color the PC. Moreover, the PCs are automatically pre-aligned directly on field without using any targets thanks to the real-time tracking of scanner movement between setups based on Visual Inertial System (VIS) [54].
At the end of the survey, the final total number of acquired points was 8,444,710,903.

Data Processing
As mentioned above, the 122 scans were automatically pre-aligned in the associated software Leica Cyclone REGISTER360 (Leica, Wetzlar, Germany) [54]. After an extensive manual cleaning operation and a coarse vegetation removal to improve the final alignment, links among scans were refined and added, obtaining a 4 mm alignment error. Considering that the survey was also conducted in two different months, September and November, the dense vegetation present in the area showed differing characteristics, adversely affecting the alignment phase. For this reason, it was removed from the project. All of the scans were then subsampled to 1 cm with a regular density reduction and subsequently merged, in order to obtain an easier to manage a unique product of 275,968,682 points. The PC was then georeferenced by means of 6 RTK GNSS points, with an RMS error of 0.036 m. Due to the change in scale, when compared to the long-range acquisition, two different approaches had been considered to separate points belonging to the ground. The first was the same CSF algorithm used with the long-range PCs, and the second was a Machine Learning (ML) approach based on the methodology developed in [55,56].
The ML approach uses a Random Forest (RF) classifier trained and evaluated on a manually annotated portion of the dataset (17% of the total number of points) composed of a total of nine semantic classes, together with computed geometric covariance features. The complexity of the area (presence of ancient ruins, caves, dense vegetation) proved a real challenge to the attribution of semantic meaning to the PC. The small amount of required labeled data and the speed of the training process make the RF classifier a good choice; nonetheless, the complex morphology of the surveyed area introduced a higher uncertainty during the selection of the correct training dataset. The 9 considered classes are: building, ground, bridge, fence, tie rods, electric wire, street retaining wall, trellis, and vegetation ( Figure 5).
Classification results proved satisfactory (94.56% overall acc.; 79.35% avg. recall; 82.43% avg. F1 score) and comparable to the application of the CSF algorithm. However, the time required to find the correct training set, manually segment it, compute geometric features for the whole PC, and spread the classification to the entire area made the CSF a better alternative for its faster application. This is because only the distinction among ground and non-ground points counts for the DTM generation. The ML RF approach would better suit an architectonic dataset that needs a finer classification. For these kinds Two surveys were performed between September and November 2019, reconstructing the area by means of 122 scans. The first 104 were collected during the first survey, while the last 18, which were necessary to complete the far eastern section, were collected during the second one.
Thanks to the scanner's 360 • horizontal and 300 • vertical field of view, even the most closed and narrow environments have been completely acquired without any problems. The acquisition is very fast because the scanner is self-leveling and every scan only requires less than 3 minutes, setting the highest resolution of 3 mm at 10 m. The simultaneous acquisition of full spherical HDR images through the incorporated 36 Mpixel camera allows for obtaining the RGB information to color the PC. Moreover, the PCs are automatically pre-aligned directly on field without using any targets thanks to the real-time tracking of scanner movement between setups based on Visual Inertial System (VIS) [54].
At the end of the survey, the final total number of acquired points was 8,444,710,903.

Data Processing
As mentioned above, the 122 scans were automatically pre-aligned in the associated software Leica Cyclone REGISTER360 (Leica, Wetzlar, Germany) [54]. After an extensive manual cleaning operation and a coarse vegetation removal to improve the final alignment, links among scans were refined and added, obtaining a 4 mm alignment error. Considering that the survey was also conducted in two different months, September and November, the dense vegetation present in the area showed differing characteristics, adversely affecting the alignment phase. For this reason, it was removed from the project. All of the scans were then subsampled to 1 cm with a regular density reduction and subsequently merged, in order to obtain an easier to manage a unique product of 275,968,682 points. The PC was then georeferenced by means of 6 RTK GNSS points, with an RMS error of 0.036 m. Due to the change in scale, when compared to the long-range acquisition, two different approaches had been considered to separate points belonging to the ground. The first was the same CSF algorithm used with the long-range PCs, and the second was a Machine Learning (ML) approach based on the methodology developed in [55,56].
The ML approach uses a Random Forest (RF) classifier trained and evaluated on a manually annotated portion of the dataset (17% of the total number of points) composed of a total of nine semantic classes, together with computed geometric covariance features. The complexity of the area (presence of ancient ruins, caves, dense vegetation) proved a real challenge to the attribution of semantic meaning to the PC. The small amount of required labeled data and the speed of the training process make the RF classifier a good choice; nonetheless, the complex morphology of the surveyed area introduced a higher uncertainty during the selection of the correct training dataset. The 9 considered classes are: building, ground, bridge, fence, tie rods, electric wire, street retaining wall, trellis, and vegetation ( Figure 5). of datasets, the geometric features distribution among the same class is more uniform and allows the algorithm to better distinguish among different elements. For instance, columns have a certain curvature, and it is remarkably different from that of a floor, thus allowing the algorithm to clearly distinguish among them. For the case study at hand, the covariance features that describe bare ground points are easily confused with those describing rocks and caves, introducing a strong degree of confusion. The CSF algorithm was tested on a 5 cm subsampled PC (16,160,360 points). The complexity of the archaeological site demanded a large number of tests before finding the optimal parameters and a significant part of manual labor to select all the correct ground points. Considering the operation of the algorithm, to improve the result, those environments like caves and similar were removed before applying it. The resulting points, 4,272,027, were obtained, setting both the cloth resolution and the classification threshold equal to 0.10 m. Subsequently, they were interpolated in ArcGIS Pro, setting a 0.10 m raster cell and using the IDW algorithm again for the reasons explained above ( Figure 6B). The WGS84 UTM 32N projection was assigned to the output raster. Classification results proved satisfactory (94.56% overall acc.; 79.35% avg. recall; 82.43% avg. F1 score) and comparable to the application of the CSF algorithm. However, the time required to find the correct training set, manually segment it, compute geometric features for the whole PC, and spread the classification to the entire area made the CSF a better alternative for its faster application. This is because only the distinction among ground and non-ground points counts for the DTM generation. The ML RF approach would better suit an architectonic dataset that needs a finer classification. For these kinds of datasets, the geometric features distribution among the same class is more uniform and allows the algorithm to better distinguish among different elements. For instance, columns have a certain curvature, and it is remarkably different from that of a floor, thus allowing the algorithm to clearly distinguish among them. For the case study at hand, the covariance features that describe bare ground points are easily confused with those describing rocks and caves, introducing a strong degree of confusion.
The CSF algorithm was tested on a 5 cm subsampled PC (16,160,360 points). The complexity of the archaeological site demanded a large number of tests before finding the optimal parameters and a significant part of manual labor to select all the correct ground points. Considering the operation of the algorithm, to improve the result, those environments like caves and similar were removed before applying it. The resulting points, 4,272,027, were obtained, setting both the cloth resolution and the classification threshold equal to 0.10 m. Subsequently, they were interpolated in ArcGIS Pro, setting a 0.10 m raster cell and using the IDW algorithm again for the reasons explained above ( Figure 6B). The WGS84 UTM 32N projection was assigned to the output raster.

Backpack LS IMMS Heron MS Twin
Moving towards the hillside, in the areas where dense vegetation prevails and the high incidence of long-range TLS rays and the morphological characteristics of the slope do not allow a proper terrain reconstruction, a Backpack IMMS campaign is done both (i) to obtain higher resolution data of a territory otherwise inaccessible and hidden using other survey techniques, and (ii) to simultaneously test the geometric accuracy that is achievable by an instrument relying on SLAM (Simultaneous Localization And Mapping) technology and strongly tested in such a challenging scenario, with continuous changes of slopes and constant presence of vegetation-trees and evergreen shrubs [57,58].
The surveyed territory is part of a natural reserve. Here, the glacial morphology, characterized by steep slopes, deep cracks, and cliffs, is particularly evident. Ancient caves called "Trone" present along the path are also surveyed. All these factors constitute an excellent test bench for an IMMS survey that is strongly tested since the typical recommendations implemented when dealing with SLAM algorithms are missing. Commonly, in order to refine the point cloud registration along the trajectory, it is recommended to set up a survey scheme made up of closed loops in order to create some rigid constraints that allow for adjusting both horizontal and vertical errors, following a principle similar to the one of the closed polygonal. This means starting and ending the survey at the same point and passing through the same portion of the surveyed environment at different times. The path is about 3 km long with an altitude difference of approximately 350 m, starting from 700 masl and ending at 360 masl (highlighted in red in Figure 1). It took two different days, one in October 2020 and the other one in March 2021, because of storage issues faced during the first day. The Backpack IMMS Heron MS Twin, produced by Gexcel [59], was used to perform the survey (Figure 7). The instrument is characterized by two LiDAR Velodyne Puck sensors, one rotating on the vertical axis and the other one on a 45° tilted axis, which provide 360° horizontal FOV (Field Of View) and 30° vertical FOV each. While walking, it is possible to check in real time the surveyed environment thanks to a portable tablet that can be carried with a dedicated harness, leaving the hands free.

Backpack LS IMMS Heron MS Twin
Moving towards the hillside, in the areas where dense vegetation prevails and the high incidence of long-range TLS rays and the morphological characteristics of the slope do not allow a proper terrain reconstruction, a Backpack IMMS campaign is done both (i) to obtain higher resolution data of a territory otherwise inaccessible and hidden using other survey techniques, and (ii) to simultaneously test the geometric accuracy that is achievable by an instrument relying on SLAM (Simultaneous Localization And Mapping) technology and strongly tested in such a challenging scenario, with continuous changes of slopes and constant presence of vegetation-trees and evergreen shrubs [57,58].
The surveyed territory is part of a natural reserve. Here, the glacial morphology, characterized by steep slopes, deep cracks, and cliffs, is particularly evident. Ancient caves called "Trone" present along the path are also surveyed. All these factors constitute an excellent test bench for an IMMS survey that is strongly tested since the typical recommendations implemented when dealing with SLAM algorithms are missing. Commonly, in order to refine the point cloud registration along the trajectory, it is recommended to set up a survey scheme made up of closed loops in order to create some rigid constraints that allow for adjusting both horizontal and vertical errors, following a principle similar to the one of the closed polygonal. This means starting and ending the survey at the same point and passing through the same portion of the surveyed environment at different times. The path is about 3 km long with an altitude difference of approximately 350 m, starting from 700 masl and ending at 360 masl (highlighted in red in Figure 1). It took two different days, one in October 2020 and the other one in March 2021, because of storage issues faced during the first day. The Backpack IMMS Heron MS Twin, produced by Gexcel [59], was used to perform the survey (Figure 7). The instrument is characterized by two LiDAR Velodyne Puck sensors, one rotating on the vertical axis and the other one on a 45 • tilted axis, which provide 360 • horizontal FOV (Field Of View) and 30 • vertical FOV each. While walking, it is possible to check in real time the surveyed environment thanks to a portable tablet that can be carried with a dedicated harness, leaving the hands free. The automatic acquisition at 15 Hz of spherical images (1920 × 1080 pixel) makes it possible to color the PC. A total of 7 trajectories were acquired, stopping the acquisition when the operator needed to rest and took about two hours to complete the survey. Their main characteristics are summed up in Table 2. As said before, the slope conformation did not allow for creating well-structured rigid constraints along the trajectory. However, closed loops were built where feasible, precisely in the central part of the trajectory, going off the path without putting the operator in trouble.
ble to color the PC. A total of 7 trajectories were acquired, stopping the acquisition when the operator needed to rest and took about two hours to complete the survey. Their main characteristics are summed up in Table 2. As said before, the slope conformation did not allow for creating well-structured rigid constraints along the trajectory. However, closed loops were built where feasible, precisely in the central part of the trajectory, going off the path without putting the operator in trouble.
Some portions of the surveyed path, reported as an example in Figure 7C, were also quite challenging to survey: the presence of dense stuck branches risked damaging the instrument, with the consequent risk of losing the calibration of the cameras on the top of the pole. In such cases, a second person present in the field helped prevent dangerous situations, alerting the operator when needed.
While walking to reach the upper part of the hillside from which to start the survey, a total of 5 RTK GNSS points were taken in those areas where the cover provided by the vegetation made it possible. Although limited, they are essential for a proper trajectory reconstruction.    Some portions of the surveyed path, reported as an example in Figure 7C, were also quite challenging to survey: the presence of dense stuck branches risked damaging the instrument, with the consequent risk of losing the calibration of the cameras on the top of the pole. In such cases, a second person present in the field helped prevent dangerous situations, alerting the operator when needed.
While walking to reach the upper part of the hillside from which to start the survey, a total of 5 RTK GNSS points were taken in those areas where the cover provided by the vegetation made it possible. Although limited, they are essential for a proper trajectory reconstruction.

Data Processing
The 7 trajectories were reconstructed and processed with Heron Desktop software (Gexcel, Brescia, Italy) [59], following the implemented four-step workflow. The numerous parameters that were set throughout the procedure gave the possibility to operate at will, according to the sensitivity of the operator. Once estimated, the trajectory by checking the accordance between the IMU (Inertial Measurement Unit) and the alignment of the scans (Figure 8A), the third part of the workflow involves the drift effect reduction by creating links all along the trajectory. In this step, it is possible to insert static scans in order to better constrain the IMMS PC. For this purpose, four scans were acquired with TLS FARO Focus 3D × 330 in correspondence with the final part of the surveyed path, paying attention to detect the most possible three-dimensional geometries in order to produce even more rigid constraints. Before being imported into Heron Desktop, they were georeferenced by means of 12 RTK GNSS points, achieving a mean registration error of 0.028 m.  In this phase of the process, it was possible to select specific points in the PC and assign known coordinates and the relative accuracy, named "Heron accuracy" in Table 3. The 5 RKT GNSS points were then added as constraints in order to adjust the trajectory. Specifications of the used points are summarized in Table 3. Table 3. RTK GNSS points used to adjust the Heron MS Twin trajectory. The mentioned optimization refers to the iterative process that concludes the trajectory reconstruction phase, in which the creation of rigid links among the IMMS PC and, if available, static scans, and the use of points with known coordinates allow a proper reconstruction.

Point
Heron Accuracy (m)

RTK GNSS Accuracy (m)
Error (m) At the end of an optimization iterative process, in which every link created to strengthen the trajectory is checked and, if needed, refined, the trajectories are exported to Reconstructor® software (Gexcel, Brescia, Italy) [59]. During this phase, it is possible to remove ghost points. This is very useful in case, for example, there are people walking near the surveyor. It is also possible to select a proper range around the sensor in order to export only the useful points. In the present case, a 1-30 m range was selected. Finally, a 0.05 m voxelization was chosen to speed up processing. The exported PC is reported in Figure 8B.
In order to separate ground/non-ground points, in this case, the CSF algorithm was also adopted. Due to the almost constant slope of the surveyed path, the PC was processed all at once. A cloth resolution of 0.20 m and a classification threshold of 0.50 m has ensured the preservation of rocks above the ground level, required for a correct representation of the territory. However, such a threshold has resulted in an incomplete trees removal: the part of the trunk closest to the ground was preserved. For this reason, the CSF algorithm was implemented again on the ground points that were just extracted. This time, a cloth resolution of 0.10 m and a classification threshold of 0.10 m were adopted, allowing for obtaining a correct result.
The obtained 15,707,444 points were interpolated in ArcGIS Pro setting a 0.10 m raster cell and using again the IDW algorithm to compute the DTM. The WGS84 UTM 32N projection was assigned to the output raster, reported in Figure 9.

Unique Multi-Resolution DTM
The last step necessary to obtain a unique multi-resolution DTM is to integrate the different terrain models processed so far. To do this, a mosaic is created using the homonymous function available in ArcGIS Pro. The function returns a unique grid whose resolution coincides with the minimum DTM cell; in this case, it was 0.10 m. The projected reference system WGS84 UTM 32N is associated with the output raster.
A table listing pros and cons of the employed instrumentation is reported in Table 4, while a flowchart summarizing the whole adopted methodology is provided in Figure 10.

Unique Multi-Resolution DTM
The last step necessary to obtain a unique multi-resolution DTM is to integrate the different terrain models processed so far. To do this, a mosaic is created using the homonymous function available in ArcGIS Pro. The function returns a unique grid whose resolution coincides with the minimum DTM cell; in this case, it was 0.10 m. The projected reference system WGS84 UTM 32N is associated with the output raster.
A table listing pros and cons of the employed instrumentation is reported in Table 4, while a flowchart summarizing the whole adopted methodology is provided in Figure 10.   Figure 10. Flowchart summarizing the entire adopted methodology to obtain the multi-resolution DTM.

Results
A quality assessment is reported for each one of the three generated DTMs. Vertical and horizontal accuracies are reported in Table 5. They are calculated by means of standard deviation of residuals, as reported in Equation (1): where SD refers to Standard Deviation, α corresponds to X, Y, Z coordinates and Δαi = αi,DTM − αi,CHECK Figure 10. Flowchart summarizing the entire adopted methodology to obtain the multi-resolution DTM.

Results
A quality assessment is reported for each one of the three generated DTMs. Vertical and horizontal accuracies are reported in Table 5. They are calculated by means of standard deviation of residuals, as reported in Equation (1): where SD refers to Standard Deviation, α corresponds to X, Y, Z coordinates and The extensive RTK GNSS campaign, resulting in the acquisition of almost 200 points over the entire valley floor of the test area, allowed for computing the SD on a representative sample of check points. ∆ Z of the three DTMs is reported in Figure 11.
In Figure 11A, check points for which ∆z exceeds the threshold value for a 1:2000 scale representation are colored in red. Threshold value is set in a conservative way at 0.50 m. Those points correspond to RTK GNSS points surveyed in the Belfort archaeological area, which has not been properly reconstructed in the 0.50 m resolution DTM. Considering that the final unique multi-resolution DTM has the accuracy attributable to the Leica RTC360 DTM in the Belfort archaeological site, if red points are removed from the SD calculation, it decreases to 0.176 m, respecting the scale representation error. In Figure 11A, check points for which Δz exceeds the threshold value for a 1:2000 scale representation are colored in red. Threshold value is set in a conservative way at 0.50 m. Those points correspond to RTK GNSS points surveyed in the Belfort archaeological area, which has not been properly reconstructed in the 0.50 m resolution DTM. Considering that the final unique multi-resolution DTM has the accuracy attributable to the Leica RTC360 DTM in the Belfort archaeological site, if red points are removed from the SD calculation, it decreases to 0.176 m, respecting the scale representation error.
Another important vertical accuracy check for the 0.50 m resolution DTM is carried out on the south side, in the red-highlighted area in Figure 1. Here, RTK GNSS points were collected to improve the Backpack IMMS PC accuracy. SD is calculated on the three points acquired while walking and is equal to 6.525 m. This kind of result is not unexpected. In this regard, some consideration must be made. Vegetation is very dense in that area and the incidence of LS is high due to the scan position. Even if RIEGL Vz-4000 registers multiple echoes, the probability of reaching the ground is almost null. Moreover, the distance from the scan position causes a low points density. For this reason, the interpolated DTM resulted in an incorrect ground reconstruction, degrading the vertical accuracy when moving towards the side.
Regarding the DTM produced with Backpack IMMS Heron MS Twin, by looking at the plot reported in Figure 11C, the accuracy seems adequate for the obtained product. Another important vertical accuracy check for the 0.50 m resolution DTM is carried out on the south side, in the red-highlighted area in Figure 1. Here, RTK GNSS points were collected to improve the Backpack IMMS PC accuracy. SD is calculated on the three points acquired while walking and is equal to 6.525 m. This kind of result is not unexpected. In this regard, some consideration must be made. Vegetation is very dense in that area and the incidence of LS is high due to the scan position. Even if RIEGL Vz-4000 registers multiple echoes, the probability of reaching the ground is almost null. Moreover, the distance from the scan position causes a low points density. For this reason, the interpolated DTM resulted in an incorrect ground reconstruction, degrading the vertical accuracy when moving towards the side.
Regarding the DTM produced with Backpack IMMS Heron MS Twin, by looking at the plot reported in Figure 11C, the accuracy seems adequate for the obtained product. However, it provides partial information since there are also those used for the improvement of the accuracy in the Heron Desktop software optimization phase among the points used for the calculation of the SD. A global vertical accuracy check is therefore conducted by subtracting the 5 m resolution DTM, made available by the institutions, to the 0.10 m resolution deriving from the Backpack IMMS. The decision to use the 5 m DTM instead of the 0.50 m one comes from the observation previously made about the incorrect reconstruction of the ground in this densely vegetated area. The histogram of the DTMs differences is shown in Figure 12.
Numerous considerations can be made, both about the processing and the obtained results. As for the use of CSF algorithm to separate ground/non-ground points, no manual operation was required to better perform the algorithm when applied to the long-range TLS RIEGL Vz-4000 PC, considering that data were acquired from higher station points with respect to the height of the ground. However, when dealing with the PCs coming from the medium-range TLS Leica RTC360 and Backpack IMMS Heron MS Twin, caves and rooms situated underground were previously manually segmented, losing some of the benefit of using an automated points extraction algorithm. In fact, the presence of underground areas makes it difficult to correctly select points belonging to the terrain since the algorithm turns the PC upside down.
used for the calculation of the SD. A global vertical accuracy check is therefore con by subtracting the 5 m resolution DTM, made available by the institutions, to the resolution deriving from the Backpack IMMS. The decision to use the 5 m DTM in the 0.50 m one comes from the observation previously made about the incorrec struction of the ground in this densely vegetated area. The histogram of the DTM ences is shown in Figure 12. Numerous considerations can be made, both about the processing and the o results. As for the use of CSF algorithm to separate ground/non-ground points, no operation was required to better perform the algorithm when applied to the lon TLS RIEGL Vz-4000 PC, considering that data were acquired from higher station with respect to the height of the ground. However, when dealing with the PCs from the medium-range TLS Leica RTC360 and Backpack IMMS Heron MS Twi and rooms situated underground were previously manually segmented, losing the benefit of using an automated points extraction algorithm. In fact, the presenc derground areas makes it difficult to correctly select points belonging to the terra the algorithm turns the PC upside down.
As for the survey conducted with Backpack IMMS Heron MS Twin, tests con to reconstruct the trajectory in a proper way have demonstrated the need for poi known coordinates, acquired with a regular frequency, to be added during the op tion phase. In the case reported, the maximum drift error is obtained at a distance 350 m from RTK GNSS points 19, 20, and 21, and it consists of about 20 m in the tion. The result is explained when checking the accuracy of the obtained DTM, in 12, and is something expected considering the challenging conditions under w survey was carried out. As stated in 2.3., the SLAM algorithms need specific pra order to perform well, and they were not guaranteed.
Regarding the high-resolution DTMs, the 0.50 m DTM obtained by processing Vz-4000 PC clearly shows its value. Figure 13 compares the same areas (A-B an visible in the 0.50 m DTM (A and C) and in the 5 m DTM (B and D) made availabl Institutions. By comparison, it is evident how the high-resolution DTM allows for ing precisely the incision rate that the Mera river achieves on the sediments (Figu giving the opportunity to conduct a quantitative analysis. Figure 13A then sho river-bed traces that are essential to study the evolution of the valley. Such de As for the survey conducted with Backpack IMMS Heron MS Twin, tests conducted to reconstruct the trajectory in a proper way have demonstrated the need for points with known coordinates, acquired with a regular frequency, to be added during the optimization phase. In the case reported, the maximum drift error is obtained at a distance of about 350 m from RTK GNSS points 19, 20, and 21, and it consists of about 20 m in the Z-direction. The result is explained when checking the accuracy of the obtained DTM, in Figure 12, and is something expected considering the challenging conditions under which the survey was carried out. As stated in 2.3., the SLAM algorithms need specific practices in order to perform well, and they were not guaranteed. Regarding the high-resolution DTMs, the 0.50 m DTM obtained by processing RIEGL Vz-4000 PC clearly shows its value. Figure 13 compares the same areas (A-B and C-D) visible in the 0.50 m DTM (A and C) and in the 5 m DTM (B and D) made available by the Institutions. By comparison, it is evident how the high-resolution DTM allows for detecting precisely the incision rate that the Mera river achieves on the sediments ( Figure 13C), giving the opportunity to conduct a quantitative analysis. Figure 13A then shows old river-bed traces that are essential to study the evolution of the valley. Such details are invisible on the 5 m DTM and make the high-resolution DTM crucial for the purpose of the A.M.AL.PI. 18 project.
The 0.10 m DTM obtained by processing Leica RTC360 PC related to the Belfort archaeological site is reported in Figure 14 and compared to the same area visible on the 0.50 m DTM. In this case as well, the added value clearly comes to light. It is sufficient to observe that, in the low resolution DTM (on the right), entire portions of the terrain are missing, thus mystifying the reconstruction of the area and creating negative volumes of rock masses (dark areas in the central part of the right image).
Regarding the high-resolution 0.10 m DTM processed starting from the Backpack IMMS Heron MS Twin, in Figure 9, two details are enhanced. They refer to outcropping rocks where the shaping action of the ice flow occupying the Chiavenna valley is particularly evident. Details, such as these highlighted, are very precious because they are not visible in any other type of survey, given their location. Moreover, considering that one of the main objectives of the A.M.AL.PI. 18 project is the innovative use of natural resources, a digital reconstruction of an area comprised in a natural reserve could become crucial to establish novel attractiveness strategies.  The 0.10 m DTM obtained by processing Leica RTC360 PC related to the Belfort archaeological site is reported in Figure 14 and compared to the same area visible on the 0.50 m DTM. In this case as well, the added value clearly comes to light. It is sufficient to observe that, in the low resolution DTM (on the right), entire portions of the terrain are missing, thus mystifying the reconstruction of the area and creating negative volumes of rock masses (dark areas in the central part of the right image).   The 0.10 m DTM obtained by processing Leica RTC360 PC related to the Belfort archaeological site is reported in Figure 14 and compared to the same area visible on the 0.50 m DTM. In this case as well, the added value clearly comes to light. It is sufficient to observe that, in the low resolution DTM (on the right), entire portions of the terrain are missing, thus mystifying the reconstruction of the area and creating negative volumes of rock masses (dark areas in the central part of the right image).  Along with these virtuous examples, however, it is also necessary to identify areas that have not been properly reconstructed. As mentioned in Section 2.1, one of the scans conducted with long-range TLS was characterized by very dense vegetation in the area of the alluvial fan formed by the Valle Drana stream. Together with the high incidence angle of the laser beam, it has resulted in an incorrect ground reconstruction, detectable in the 0.50 m DTM, as reported in Figure 15. The points classified with the CSF algorithm have a low density and, given the circumstances of the survey, could easily belong to the crowns of the trees. that have not been properly reconstructed. As mentioned in Section 2.1., one of the scans conducted with long-range TLS was characterized by very dense vegetation in the area of the alluvial fan formed by the Valle Drana stream. Together with the high incidence angle of the laser beam, it has resulted in an incorrect ground reconstruction, detectable in the 0.50 m DTM, as reported in Figure 15. The points classified with the CSF algorithm have a low density and, given the circumstances of the survey, could easily belong to the crowns of the trees.

Discussion and Conclusions
The paper presents the first part of an ambitious project aimed at realizing a best practice to study historical behavior of dangerous hillsides. Considering the goal of the A.M.AL.PI. 18 project, the presented research includes the creation of high-resolution multi-source 3D terrain models obtained by integrating 3D data coming from different types of Laser Scanning instruments.
Numerous challenges were faced, both due to the characteristics of the test area and to the processed Point Clouds (PCs). Regarding data coming from the long-range Terrestrial Laser Scanner (TLS), the valley conformation has caused the application of the ground points filtering procedure three times in order to also preserve points in the steepest regions. This resulted in tripling the processing time, which was already onerous, given the amount of data to be analyzed. In addition, the dense vegetation present during the second survey and the incidence of the laser scanner rays in a portion of the valley floor gave rise to artifacts that are hard to remove, as illustrated in Figure 15.
Processing and managing the scans acquired with the medium-range TLS was not easy. A huge calculation capacity was required, due to their number and size. Moreover, the presence of caves and other underground spaces demanded a significant manual labor when running the Cloth Simulation Filter (CSF) algorithm to obtain all the correct ground points. The additional test involving the Machine Learning (ML) Random Forest (RF) classification to extract ground points did not bring noticeable advantages, due to the environmental scale of the archaeological site with presence of vegetation. This increases the degree of confusion of the geometric covariance features among the different classes.

Discussion and Conclusions
The paper presents the first part of an ambitious project aimed at realizing a best practice to study historical behavior of dangerous hillsides. Considering the goal of the A.M.AL.PI. 18 project, the presented research includes the creation of high-resolution multi-source 3D terrain models obtained by integrating 3D data coming from different types of Laser Scanning instruments.
Numerous challenges were faced, both due to the characteristics of the test area and to the processed Point Clouds (PCs). Regarding data coming from the long-range Terrestrial Laser Scanner (TLS), the valley conformation has caused the application of the ground points filtering procedure three times in order to also preserve points in the steepest regions. This resulted in tripling the processing time, which was already onerous, given the amount of data to be analyzed. In addition, the dense vegetation present during the second survey and the incidence of the laser scanner rays in a portion of the valley floor gave rise to artifacts that are hard to remove, as illustrated in Figure 15.
Processing and managing the scans acquired with the medium-range TLS was not easy. A huge calculation capacity was required, due to their number and size. Moreover, the presence of caves and other underground spaces demanded a significant manual labor when running the Cloth Simulation Filter (CSF) algorithm to obtain all the correct ground points. The additional test involving the Machine Learning (ML) Random Forest (RF) classification to extract ground points did not bring noticeable advantages, due to the environmental scale of the archaeological site with presence of vegetation. This increases the degree of confusion of the geometric covariance features among the different classes.
The Backpack Indoor Mobile Mapping System (IMMS) survey was a very interesting test to verify how to improve the accuracy of the system, considering that the typical recommendations to be implemented when dealing with Simultaneous Localization And Mapping (SLAM) algorithm were missing. In addition, since the path was entirely reconstructed thanks to two surveys performed during different periods of the year, the matching algorithm to strengthen the trajectory was incredibly stressed. This was due to the presence-absence of leaves on the trees, resulting in an accuracy degradation.
The added value that the final multi-resolution Digital Terrain Model (DTM) brings to the research is clear and broadly demonstrated. Future works aimed at improving the result could regard an ad hoc survey of the alluvial fan formed by the Valle Drana stream. Considering the accessibility of the area thanks to numerous pathways, a Backpack IMMS survey could be the best choice.
Regarding the accuracy check, since the RTK GNSS campaign that has been conducted up until now mainly focuses the valley floor, it would be useful and proper to include both sides in a forthcoming survey.
Lastly, a photogrammetric approach could be added to the SLAM algorithms at the base of the Backpack IMMS to improve the achievable accuracy in those environments where the typical SLAM recommendations are hard to implement. Visual SLAM performs well. By combining the two technologies, promising results could be obtained.