AIRBORNE LASER SCANNING CHANGE DETECTION FOR QUANTIFYING GEOMORPHOLOGICAL PROCESSES IN HIGH MOUNTAIN REGIONS

: Disturbance by geomorphic processes is a key factor for current and potential plant species distributions in high mountain regions. We implemented and tested an approach for the detection, quantification, and geomorphological classification of 3D topographic change using airborne laser scanning (ALS) data. This approach is applied to two point cloud epochs (from 2006 and 2017) to identify and analyse changes related to five different process categories in a study area in the European Alps: Fluvial Processes , Rock Glacier Movement & Landslides , Channel Erosion & Debris Flow , Rockfall Release & Deposition and Vegetation & Anthropogenic Change . The results are assessed through comparison with a manually produced geomorphological map. The analysis covering the eleven-year time period provides detailed information on the magnitudes and spatial distribution of change per identified geomorphological process. The study shows that the workflow is capable of providing the fundamental basis for mapping vegetation disturbance by geomorphic activity in a high alpine site.


INTRODUCTION
Disturbance by geomorphic processes (such as erosion or deposition, landsliding and rockfall) is a key factor for current and potential plant species distributions in high mountain regions (Gentili et al. 2013, Eichel et al., 2018. Thus, geomorphic process dynamics and their spatial distribution must be considered in addition to microlimatic variation and snow cover effects in order to understand the site-specific growing conditions, the vegetation succession and to estimate the rescue potential of the high-mountain flora in a warming macroclimate. High-resolution topographic data has proven its usefulness in numerous studies of Earth-surface processes at a wide range of scales (Passalacqua et al., 2015) and has the potential to provide the fundamental basis for mapping vegetation disturbance by geomorphic activity. With repeat acquisitions of airborne laser scanning (ALS) becoming increasingly available, this opens possibilities to investigate morphodynamics through topographic change detection for large areas (i.e. several km²). However, the small magnitude of many geomorphological processes makes it challenging to design a standardised and automated procedure that is reproducible and allows comparable and reliable results for different study areas . This work presented here aims to automatically detect, quantify and classify active geomorphological processes in a high-mountain area using repeat ALS point clouds and to compare the results against independently mapped geomorphological reference data.
The paper is divided into five main sections. First, in section 2 related work is presented. Section 3 presents the main methods used in our workflow and gives an overview of the datasets and the study area. Section 4 gives an insight into the accuracy of the datasets and presents the detected topographic changes and the geomorphological map. A detailed comparison of the two results is presented and the results are examined regarding possible * Corresponding author sources of error and the level of detection. Section 5 summarises the results, illustrates the practicability of the approach and identifies potential for future extensions.

Geomorphological mapping
Geomorphological mapping provides information on the spatial distribution and the type of geomorphic processes. In addition, quantitative information about the activity of processes is collected by change detection analysis comparing the state of the Earth's surface at selected points in time (Schrott and Glade, 2007). An automated way for obtaining geomorphological information is aimed at saving time, reducing workload and ensuring reproducibility of results. However, this does not exclude the need for insitu observations e.g. as validation information (Fekete et al., 2015). Otto and Smith (2013) basically distinguish two approaches in digital geomorphological mapping: (1) manual mapping and (2) automated or semiautomated mapping. Manual mapping can lead to detailed results but requires trained experts and may lack repeatability. Such a heuristic method combines visual interpretation of the terrain to identify landforms using orthophotos, terrain models and their derivatives. The second, more objective and reproducible (semi-)automated mapping approach focuses on the application of algorithms to delineate and distinguish various landforms (Otto and Smith, 2013). Morphometric features such as slope, curvature in different directions, flow path length, flow accumulation and elevation differences in the terrain are used for this purpose. The presented study primarily focuses on the latter approach and briefly compares it with results obtained with the first approach (manual mapping).

Topographic change detection and semantic labelling
Calculating elevation differences and volumetric changes from two high-resolution digital elevation models (DEMs) is straightforward and well established in geomorphology (Passalacqua et al., 2015). This DEM-of-difference (DoD) approach to topographic change detection has been applied to rasterized laser scanning data from repeat airborne (Okyay et al., 2019), terrestrial (e.g. Baldo et al., 2009) and unmanned aerial vehicle (e.g. Mayr et al. 2019) surveys. However, working directly on the 3D point clouds obtained by laser scanning, is increasingly preferred over raster-based approaches as it preserves the high geometric information content needed to reliably detect even low-magnitude topographic change (Bernard et al., 2021). Lague et al (2013) developed the Multiscale Model to Model Cloud Comparison (M3C2) method for calculating point cloud distances orthogonally to the local terrain surface as a measure for topographic change.
Beyond the detection of topographic change, a semantic labelling of processes or landforms (i.e. geomorphic objects) related to the change is a crucial step for higher-level information extraction in many applications. An automation of this task is especially challenging in complex environments and with a diverse set of processes potentially driving the change. In a raster-based approach, Anders et al. (2013) applied object based image analysis (OBIA) combining volumetric changes with a rulebased classification to analyse change per geomorphic process or landform. To distinguish erosion and deposition from tree growth or removal, Mayr et al. (2019a) combined a DoD with a machine learning classification of land cover (using laser reflectance and color orthophotos). Specifically for fluvial environments, Kasprak et al. (2017) presented two rule-based approaches to process mapping based on repeat DEMs, their differential and morphometric derivatives as well as object shape and relative orientation.
Regarding terrestrial laser scanning (TLS) in geomorphological monitoring and mapping applications, a research focus has been on semantic labelling of point clouds by machine learning. This includes the use of machine learning on geometric features (Brodu and Lague, 2012), machine learning on geometric features, laser return intensity and M3C2 distances (Weidner et al., 2019), and machine learning on geometric features combined with a rule-based topological refinement of the classification in an object-based framework (Mayr et al., 2019b). Recently, approaches for labelling point clouds or change detected therein have successfully integrated also color features (Weidner et al., 2021) and temporal features (Anders et al., 2020) extracted from photogrammetric and permanent terrestrial laser scanning point clouds, respectively. Larger-scale studies based on airborne data usually lack such features since ALS surveys are repeated rarely and with a time lag to orthophoto acquisitions. One of the few studies detecting topographic change from ALS point cloud distances is presented in Bernard et al. (2021), focusing on mapping landslide source areas and deposits in a 5-km² area. In the presented study, we apply a similar approach to detect a wider range of geomorphic process types in a 30-km² high-mountain area. We furthermore derive a geomorphological map and information for changes on a process basis from the distance calculation, using a rule-set with morphometric parameters.

Study site
The study site covers an area of ca. 30 km² around Mt. Schrankogel (3497 m a.s.l.), the second highest peak in the Stubai Alps (Tyrol, Austria). It is located in the Sulztal, a tributary valley of the Ötztal (Fig. 1). The area is characterized by metamorphic siliceous rocks and a strong geomorphic reworking by glacial, periglacial, gravitational and fluvial processes (Hoinkes et al., 2021).

Software and scripting
The processing workflow is implemented in Python 3.7.0 using the System for Automated Geoscientific Analyses (SAGA) Version 7.10.0 and the extension Laserdata LiS (Release 6432, Date 2021-04-16), which provides tools for topographic LiDAR point cloud analysis (Conrad et al., 2015;Rieg et al., 2014). The outputs are visualised with GRASS GIS for control purposes (GRASS Development Team, 2017). Figure 2 gives an overview of the implemented workflow showing the combination of the two main information sets, i.e. the ALS point cloud based topographic change and processes delineated by geomorphological mapping. On the one hand, we detected topographic changes resulting from recent process dynamics in the eleven-year time span framed by the ALS datasets from 2006 and 2017 and classify them according to the dominant processes and landforms. On the other hand, we manually produced a geomorphological map, which was used for a comparative evaluation of the semi-automatic approach. Specifically, we discuss the results with regard to their potential for assessing the geomorphic disturbance of vegetation sites (Sect. 4.2).

Figure 2.
Overall workflow and combination of the main sources of information.

Data basis
Two To restrict the performance requirements for the operating system and for an efficient data handling and data management the ALS data was divided into a graticule with 500 x 500 m tiles ( Fig. 1 and 3) as part of the data preprocessing (Fig. 4, step 1). Point cloud subsets for each tile were extracted and processed sequentially. Digital elevation models (DEMs) were derived by triangulating and gridding the Z-coordinate of laser points. In the case of multiple points in one cell, the lowest point was used for gridding. The resulting DEM 2006 and DEM 2017 had a grid resolution of 0.5 m.

Point cloud registration
An accurate co-registration of the datasets is the basis for all further processing, as it guarantees an overall data comparability. The aim is to minimize differences in the two point clouds by reducing registration errors and, thus, to improve the identification of actual changes in the terrain (Bernard et al., 2021).
For the alignment of the point clouds the iterative closest point (ICP) algorithm is applied (Besl and McKay, 1992). The point cloud from the 2017 survey with the higher point density was set as a reference during registration. The registration is carried out using selected stable areas that have neither been affected by obvious geomorphological influences, vegetation cover, nor have undergone changes due to human activity. These patches are selected to be well distributed in the area of interest, not too rough, also not overgrown with vegetation and covering a range of different orientations. A DEM of difference (DoD) based on the DEM 2006 and the DEM 2017 was calculated for determining and selecting stable areas. Moreover, the aforementioned criteria for stable areas were controlled by field observations and orthophoto interpretation. Figure 3 gives an overview of the areas used for registration and validation.
For all selected stable areas points with a 4 m buffer were created and areas with elevation differences >20 cm were excluded. Outlier points in stable area subsets were filtered by the Segmentation by Planes algorithm, which segments a point cloud by plane fitting, in this case to a minimum of three points within a neighbourhood of 2 m search radius. The maximum distance in the robust plane fitting was set 30 cm to discard outliers from the computation above this value. The ICP algorithm (Besl and McKay, 1992) was applied on the registration area subset to determine the transformation matrix, which was then applied to co-register the entire 2006 point cloud.

Detection and quantification of topographic change
After registration of the ALS 2006 to the 2017 point cloud the distance can be calculated. This detection of changes was not performed with the DoD but directly on the point clouds in order to receive more detailed results in 3D (Fig. 4, step 2). The M3C2 implementation in SAGA LIS (Fey and Wichmann, 2016) was applied to the entire dataset to create a preliminary 3D change map.
In the following, processes and landforms are mapped based on ALS point clouds, an accuracy check is carried out and outputs are created for visual control and comparison. ALS point cloud change detection handles differences data quality and data preparation by comprising detection (step 2), quantification (step 3) and characterisation of changes (step 4) (Fig. 4), as conceptually differentiated in Lague et al. (2013). This automated workflow allows to separate topographic change represented in the 3D point clouds from artefacts e.g. caused by point density differences and referencing errors. In step 3 and 4 (Fig. 4) 3D distances between the two point clouds are calculated by the M3C2 method consequently identifying a selected range of geomorphological processes. Additional validation areas were identified as explained in Section 3.5 but this time with respect to the 3D change map (Fig.  4, step 3: quantification). As an estimate of the smallest detectable change, a spatially varying level-of-detection (LOD), at the 95% confidence interval is calculated. We used the LOD calculation Eq. (1) integrated in the M3C2 implementation by Fey and Wichmann (2016).
where 2 = plane fitting variance = number of points in fitting radius = registration error = positional error , = point cloud A and B This LOD takes into account the registration error and the local plane fitting variance of both point clouds (Lague et al. 2013). To estimate the registration error between the two point clouds the M3C2 algorithm was applied on additional stable (validation) areas (Fig. 3). This preliminary registration error estimate was combined with the internal dataset quality, taking into account strip adjustment errors and the height accuracy specifications of the ALS survey. From this, a conservative estimate of the error in each point cloud was derived and used as parameter reg (registration error) to calculate the LOD. Computed distances larger than the LOD were labelled as real change.

Characterization and semi-automatic classification of topographic change
With the results of the distance calculation (Sect. 3.6) and the LOD-based thresholding of real changes, the topographic changes between 2006 and 2017 are investigated. Contiguous areas of real change were segmented by a connected components analysis and classified in combination with morphometric features (Fig. 4).

Table 1. Relevant datasets for classification and derived
information. Table 1 summarizes the datasets used for the characterisation and classification step 4 (Fig. 4). First, areas with unexpectedly high 3D distances due to interference signals and tile effects, especially on snow-and ice-covered surfaces (e.g., glacial areas), were masked out. These effects are caused by the relatively increased changes, e.g. on glacier surfaces and their volume change over the 11-year period, as well as by their specific surface reflectance-induced roughness. Subsequently, we focused on the remaining areas, which are geomorphologically relevant (investigation mask, Tab. 1).
A K-Means cluster analysis was used to create an unsupervised classification for landform elements based on morphometric features, shown in Tab. 1. As Piloyan and Konečný (2017) have shown, the geomorphometric parameters elevation, slope, profile curvature, plan curvature and flow path length are suitable for fast and efficient semi-automated landform classification. In other studies, such as Guo et al. (2021), similar parameters were found to be suitable for classifying geomorphological units. We apply K-Means cluster analysis based on iterative minimum distance (Forgy, 1965) with the maximum of 25 iterations and the input variables elevation, slope, profile curvature, plan curvature and flow path length to generate ten clusters (morphometric groups). Two of the ten clusters were deleted due to low area fractions and no meaningful assignment to a process domain, and eight clusters were reclassified as morphometric groups.
Based on the eight morphometric groups and the segmented real change layer as mask, a simple rule-set was used to identify geomorphic processes and landforms. An overview of the rulebased classification employing the morphometric parameters (Tab. 1), can be found in the supplementary material (Fig. S1). This rule-based classification allocated the change objects to the

Manual geomorphological mapping
Following the approach described by Colwell (1983) and Otto and Smith (2013), we manually produced a geomorphological map by combining field observations and photographs with visual interpretation of remote sensing data (i.e. orthophotos, shaded relief maps), and geological maps (Moser, 2011;Kreuss, 2011). The legend combines different legend concepts based on Eichel (2016) and Seijmonsbergen (2013) and is adapted to high mountain specific conditions. This manual mapping of geomorphic units was carried out independently before the digital analysis started to avoid influence of change detection results on the manual mapping process.

Registration accuracy and level of detection
At the additional stable areas used for assessment of the registration accuracy (Fig. 3), the mean point cloud distance is 0.001 m with a standard deviation of 0.023 m (Fig. 5). The typical approach would be to use this standard deviation as an estimate for the registration error across the entire study area (Fey andWichmann 2016, Bernard et al. 2021). However, a visual inspection of the calculated point cloud distances revealed distinct spatial patterns indicating local registration problems, which we attribute to internal errors in one or both point clouds (e.g. inaccurate strip adjustment where multiple strips overlap).

Figure 5.
Histogram of the 3D distances of the additional stable validation areas.
As a result, when a registration error of 0.023 m is applied, point cloud distances related to such inaccuracies are misclassified as real change in certain parts of the study area (Fig. 6). In the left map (a)), light purple areas indicate LOD values as low as 0.14 m and blue to white areas indicate high LOD values up to a maximum of 0.52 m. Yellow areas in the right map (b)) show the identified real changes. Obviously incorrectly detected real changes are highlighted and marked in the maps. An example for this is shown by the stripe-like pattern of the LOD (Fig. 6, marks  1 and 2).
For the flight campaign in 2006 a height accuracy of +/-15 cm (one standard deviation) was specified as a requirement (Federal State of Tyrol, 2011). Even if the actual point cloud accuracy is probably much better for most of the area, it seemed reasonable to apply these accuracy specifications stated by the flight operator, instead of the registration error, for testing purposes.
With the more conservative height accuracy value of 15 cm substituting the registration error (reg parameter) in the LOD calculation, previously existing problems of overestimation of the real changes were almost eliminated (Fig. 7).  With the registration error set to 15 cm a more appropriate real change result with fewer incorrectly detected areas for the whole study area can be generated on costs of change detection sensitivity.

Change detection results
Besides the location of changes, also the magnitude of change and the type of geomorphic process are of interest. Therefore, identified real changes were extracted and characterised. First, the most important changes are shown and discussed with respect to their effect in terms of vegetation disturbance (Fig. 8). This will be followed by a closer look at associated process types or landforms (Fig. 8). A detailed change map and the classified real changes can be found in the Supplementary material (Fig. S2, S4 and S5).
In Figure 8 the changes in the Schrankar on the northwest side of Schrankogel are depicted more closely. Blue areas show positive changes and red areas show negative changes in the distance measurement (Fig. 8a). Larger areas belonging to a rock glacier could be correctly classified by the rule-set and morphometric features as Rock Glacier Movement & Landslides. The rock glacier starts probably a bit higher up but, there, it is classified as Channel Erosion & Debris Flow (Fig. 8b). This is correct as such processes dominate (at least) on the surface but also demonstrates that the assignment of a single class can be problematic (depending on the application), since two or more classes may overlap. At the steep front of the rock glacier, signs of continuous erosion were observed in the field, precluding any establishment of a vegetation cover. However, topographic changes do not always relate to erosion and deposition, as it is typical for Channel Erosion & Debris Flows. The case of rock glaciers shows that geometric surface change can also result from material creeping and moving downhill or, possibly, a loss of ice content. This does not always result in disturbance of surface material and, therefore, affect the vegetation. Despite identified changes (red and blue areas in Figure 8a), vegetation patches exist on parts of the rock glacier. In the upper slope area, where changes due to channel erosion and deposition predominate, this is not the case, as the surface and the vegetation are considerably disturbed by sediment transfer. This illustrates that geometric surface change (e.g. subsidence by rock glacier degradation or movement) does not necessarily result in strong vegetation disturbance and, thus, needs a more differentiated analysis. Figure 8 also shows a problem that occurred during the detection of changes. Due to interannual variation of snow patches or different acquisition dates of the ALS data, areas of snow cover differed i.e., remain at higher altitudes or in shaded areas until late summer (see orthophoto, Fig. 8). This additionally leads to areas misclassified as real change, which further affects the correct assessment of the point cloud distances.
Changes per process are presented below as histograms (Fig. 9). The maximum surface changes were limited by the maximum search distance <2 m for all process classes. The process type Alluvial & River Area in plot 1 shows less points with negative than with positive change, a similar frequency distribution can only be found in plot 5 (Fig. 9). Note that the y-axis scale in plot 1 is 10 times lower compared to plot 2 Rock Glaciers & Landslides (Fig. 9). Rock Glaciers & Landslides have the largest share of changes identified here, followed by Channel Erosion & Debris Flow (Fig. 9, (2) and (3)). The class Vegetation & Anthropogenic Change summarizes all changes that are not caused by natural, abiotic geomorphic processes, and is included for completeness. While some of these changes are still difficult to interpret without closer inspection, the relatively sharp peak of point cloud distance frequencies around 0.35 m can be explained by tree growth, as the spatial distribution largely matches the forest extent ( Fig. 9, (5)).

Geomorphological inventory and comparison of results
Geomorphological processes were manually mapped at a scale of 1:5000 as described in section 2.2 and 3.8 (Fig. S3). Processes identified by the semi-automatic change detection and classification were verified by comparison to the geomorphological map (Fig. 10). Overlying processes could be assessed less precisely in the manual mapping. In the example of Figure 8 it has already been pointed out that Channel Erosion & Debris Flow areas can also affect rock glacier areas, i.e. the classes are not mutually exclusive. In the end, the comparison of the two results is only suitable as a simple validation possibility to check the meaningfulness of the classification. In future, not only more sophisticated approaches for a (semi-)automatic classification of topographic change according to process domains but also improved methods for validation should be developed. In cases where additional data (such as individual flight strips and trajectory information) is available for the point clouds, the LOD could be further improved (i) by a refinement of the strip adjustment, (ii) by a quantification of the positional uncertainties at point level, and (iii) by a filtering of the points with high positional uncertainties (Mayr et al. 2020). Furthermore, a correction of laser return intensity values to obtain reflectances would be interesting for derivation of a snow mask (e.g., Höfle et al. 2007).

CONCLUSIONS
We present a workflow that detects and quantifies changes between two airborne laser scanning (ALS) point clouds and, subsequently, assigns them to different geomorphological process types. Overall, the tested approach demonstrates good potential for an efficient mapping, quantification, and geomorphological interpretation of topographic change. The results indicate the spatial distribution of changes and their process affiliation, which is an important information for botanical and ecological research in high-mountain environments.
The outcomes will be further used for detecting areas where geomorphic disturbance affects the conditions for vegetation establishment and growth. A differentiated assessment of the magnitude and frequency of geomorphic disturbance from bitemporal ALS point cloud distances alone remains difficult. However, the proposed process-oriented classification of topographic changes can inform on locally dominant process types, and future work could use this for a heuristic assessment of the frequency and degree of vegetation disturbance. Thus, we plan to further improve the methods for automated process type classification, amongst others by a tighter integration of geometric and spectral features, allowing us to extend the analyses to a regional scale.

SUPPLEMENTARY MATERIALS
Supplementary Materials: The following are available online at https://doi.org/10.5281/zenodo.6402513. Figure  S1: Classification preparation and decision tree. Figure S2: 3D distance change and change offset overall map. Figure S3: Geomorphological map inventoried by manual mapping. Figure  S4: Classified real changes of the valley area between 2006 and 2017. Figure S5: Classified real changes between 2006 and 2017.