MULTITEMPORAL ANALYSIS OF OBJECTS IN 3 D POINT CLOUDS FOR LANDSLIDE MONITORING

To date multi-temporal 3D point clouds from close-range sensing are used for landslide and erosion monitoring in an operational manner. Morphological changes are typically derived by calculating distances between points from different acquisition epochs. The identification of the underlying processes resulting in surface changes, however, is often challenging, for example due to the complex surface structures and influences from seasonal vegetation dynamics. We present an approach for object-based 3D landslide monitoring based on topographic LiDAR point cloud time series separating specific surface change types automatically. The workflow removes vegetation and relates surface changes derived from a point cloud time series directly to (i) geomorphological object classes (landslide scarp, eroded area, deposit) and (ii) to individual, spatially contiguous objects (such as parts of the landslide scarp and clods of material moving in the landslide). We apply this approach to a time series of nine point cloud epochs from a slope affected by two shallow landslides. A parameter test addresses the influence of the registration error and the associated level of detection on the magnitude of derived object changes. The results of our case study are in accordance with field observations at the test site as well as conceptual landslide models, where retrogressive erosion of the scarp and downslope movement of the sliding mass are major principles of secondary landslide development. We conclude that the presented methods are well suited to extract information on geomorphological process dynamics from the complex point clouds and aggregate it at different levels of abstraction to assist landslide and erosion assessment. * Corresponding author


INTRODUCTION
Landslides and erosion represent major challenges for natural hazard management and sustainable agriculture in mountain areas (Turner et al., 1996;Alewell et al., 2015).Being a specific type of gravitational mass movements, landslides shape the landscape and they can cause damage on humans and infrastructure (Kjekstad and Highland, 2009;Petley, 2012).Furthermore, these processes lead to a loss of soil and degrade agricultural land (Wiegand and Geitner, 2010;Alewell et al., 2015).These problems concern both the occurrence of new landslides and the reactivation and secondary erosion of existing ones.
Several case studies have quantified geomorphological surface changes with multi-temporal point clouds, e.g.via a calculation of distances between points from different acquisitions (epochs) (Lague et al., 2013;Stumpf et al., 2015;Fey and Wichmann, 2017).Usually, the detected changes are analysed based on the unstructured point clouds representing a continuous surface for each epoch, leaving the interpretation of changes to a human expert.This interpretation, i.e. the identification of the underlying process for changes, however, is often challenging, for example due to the complex surface structures and influences from vegetation dynamics.
We present an approach for object-based 3D landslide monitoring with LiDAR point clouds for separating different surface change types.The highly automated workflow removes vegetation and relates surface changes derived from a point cloud time series directly to (i) geomorphologically meaningful object classes ('landslide scarp', 'eroded area', 'deposit') and (ii) to individual objects (such as parts of the landslide scarp and clods of material moving in the landslide).

TEST SITE AND DATA
We apply the presented approach with a time series of nine TLS point cloud epochs from a slope affected by two shallow landslides (Fig. 1).This test site is located in the Schmirn valley (Tyrol, Austria) at about 1700 m a.s.l., facing south-west with a slope gradient of approximately 35°.The landslide scars, excluding the rather diffuse runout zones, are at maximum 20 x 30 m in size and have a maximum depth of approximately 2 m.Accordingly, the main (shallow) mass movements, as well as secondary erosion processes, are restricted to Quaternary deposits overlying the Bündner Schist bedrock (Wiegand et al., 2013;Mayr et al., 2014).
The lower part of the landslides' surroundings is used as a meadow and the upper part as an occasional pasture.Larch trees (Larix decidua) and shrubs are scattered over some sections of the slope.Furthermore, a few small rock cliffs in the lower part add to the geometric and semantic complexity of the scene, making it a well-suited site for testing new methods for surveying and environmental monitoring.
Historical orthophotos indicate an expansion of the eroded landslide areas since their initial development, which was triggered by a heavy precipitation event in 1965.The coarse spatial and temporal resolution of the orthophotos, however, does not reveal the detailed mechanisms resulting in this ongoing erosion and preventing a revegetation and stabilisation of the eroded areas.Striving for an enhanced process understanding of such landslides, monitoring methods based on TLS are implemented by conducting scans twice a year from 2011 to 2015.To reduce occlusions, the scans were acquired from two different positions in the flat valley bottom, at a distance of approximately 200 to 350 m from the investigated slope section.For the first five epochs an Optech Ilris-3D scanner, and for the next four epochs a Riegl VZ-6000 scanner was used.Since both TLS instruments operate at different wavelengths, intensity features cannot be used comparatively within the time series of point clouds in a straight forward manner and would require correction (e.g.Kashani et al., 2015).
Geo-referencing and registration into a common reference system was done using four spherical targets, surveyed by DGNSS (differential global navigation satellite system), and iterative closest point adjustment at stable parts of the slope (ICP; Besl and McKay, 1992).For homogenization of point density and reduction of data volume, all nine point clouds have been thinned by a 3D block filter (3 cm blocks), retaining only the point closest to each block centre.

WORKFLOW AND METHODS
The presented work builds upon an automated landslide classification pipeline for 3D point clouds as a basis for object monitoring (Mayr et al., 2017).From this as a starting point our contribution is now progressing towards the identification and characterisation of landslide changes within the entire time series of classified point clouds by 4D object-based analysis.The workflow for (i) classification and (ii) multi-temporal analysis of objects in the point cloud time series is implemented with Python scripts, combined with the free open-source geographical information system (GIS) SAGA (System for Automated Geoscientific Analysis; Conrad et al., 2015) and the proprietary SAGA add-on Laserdata LIS (Rieg et al., 2014).

Point cloud classification
The classification approach (see Mayr et al., 2017) exclusively exploits geometric point cloud features.It avoids the use of colour or LiDAR intensity since these features are not always available in sufficiently high quality and consistency (see Sect. 2).Thus, the workflow starts with calculating geometric point cloud features (such as 2D z-range, 3D/2D density ratio, slope, standard deviation from a plane, omnivariance, and geometric curvature) each at three different scales (i.e.neighbourhoods of radius 0.2, 0.4 and 1.0 m).The next step is a segmentation by seeded region growing to partition the point clouds into morphologically homogeneous subsets, aiming at oversegmentation and subsequent merging by classification.
The resulting segments are labelled with the seven classes 'scarp', 'eroded area', 'deposit', 'rock outcrop' and different classes of vegetation ('low grass', 'high grass', 'medium and high vegetation') via two sequential classification steps (Fig. 1 (bottom)).First, a supervised classification by a random forest classifier (Breiman, 2001) using the geometric features is applied to the segments from the initial oversegmentation.The classifier is trained and validated with one manually labelled point cloud epoch, using one selected landslide of the test site for training and validation, respectively.After the training phase, the classifier is applied to label the entire time series of nine point cloud epochs.Subsequently, certain errors in this classification are corrected based on approximated landslide shapes (reconstructed for each epoch) and topological rules.
The approximated reconstruction of landslide shapes i) identifies the main scarp segments and ii) considers that material from below these segments moved downward (directed by gravity).Hence, a (raster terrain model based) hydrological flow routing algorithm defines the process trajectories initiated from the main scarp segments (stopping only at the lower boundary of the area-of-interest).These trajectories roughly delimit the area affected by the landslide process.For a simple topology (inside vs. outside the landslide process zone) rules are defined.These rules reclassify and correct for instance misclassified 'erosion' segments outside the two landslides.A detailed description of the point cloud classification methods is provided in Mayr et al. (2017).

Point cloud deformation calculation
Before analysing the deformation between point cloud epochs, all points labelled as 'medium and high vegetation' (see Sect. 3.1) are removed.Subsequently, a normal vector is calculated for each point, applying a RANSAC-like (Fischler and Bolles, 1981) robust plane fitting approach for the point set within a spherical neighbourhood of 0.4 m radius around the central point.
For each time step, the local distance between the points of two subsequent point cloud epochs A and B is calculated along the surface normal direction in 3D space, using a variant of the Multiscale Model to Model Cloud Comparison (M3C2) method (Lague et al., 2013;Fey and Wichmann, 2017; see also sketches therein).In our experiment, the parameters for this distance calculation have been determined empirically, as given in the following description.For each point a in the point cloud from epoch A, the nearest points in the point cloud from epoch B both along the normal vector and along the flipped normal vector of point a are determined and the one which is closest to point a is used as point b.This search is constrained by a cylinder centred on the normal vector, with radius d/2 = 0.2 m and length L = 2 m.For both points (a and b) and their respective neighbour points within a specified fitting radius D/2 = 0.5 m least-squares fitted planes are computed, using only points with a maximum deviation of 45° from the search point's (a or b) normal vector orientation.
Finally, a signed 3D distance is measured from point a (projected on its fitted plane) along the normal (or flipped normal) of point a to the intersection with the fitted plane of point b.In case the intersection point is in the direction of the normal of point a, the surface of point cloud A is behind (or below or inside) that of point cloud B and the sign is negative.A positive sign indicates that the surface of A is in front of (or above or outside) the surface of B.

Object changes
The next step identifies geomorphologically meaningful objects by aggregating spatially connected segments of the same class (from the initial oversegmentation).This object modelling is constrained by the 3D distance, calculated in the previous step, to distinguish deforming parts of the landslide from stable parts.The resulting objects are spatially contiguous and discrete entities, each belonging to one of the seven object classes.
Finally, the mean 3D deformation is calculated (i) per object and (ii) per class, from the "real" deformations (i.e. with a magnitude larger than the level of detection (LOD); see Sect.3.4) attached to each point by the distance calculation step, including a sign (indicating the direction of surface change).This is done separately for each time step and for each landslide of the test site, i.e. split by a coarsely defined area-ofinterest polygon layer.The resulting object-based and classbased changes are attached as attributes to the point cloud epochs and exported as a table for further time series analysis and detailed visualisation of changes.

Level of detection and analysis of parameters
The applied distance calculation method (Lague et al., 2013;Fey and Wichmann, 2017) also estimates a spatially variable level of detection at the 95% confidence interval (LOD 95 ).This LOD 95 helps to differentiate between "real" changes and data noise or systematic errors, and is calculated as where A, B = point cloud epochs σ² = plane fitting variance n = number of points in the fitting neighbourhood reg = registration error.
Depending on the knowledge regarding the registration accuracy, either a constant or a spatially variable registration error reg can be used.In our experiment we calculate with a fixed registration error, thus the spatial variability of the LOD 95 depends exclusively on the plane fitting variance and the number of points in the fitting neighbourhood.These two variables in turn depend on the local surface roughness, the point density and the neighbourhood radius for fitting the plane.We use a conservative estimate for the registration error reg = 0.10 m (and reg = 0.06 m, respectively, for the deformations displayed in Fig. 5 to enhance the visibility).
Additionally, we calculate the distances within a defined range of values for the registration error (reg) starting from 0.2 m to 0.16 m with 0.2 m steps.The resulting object-based 3D deformations are aggregated and investigated exemplarily for the 'landslide scarp' class to illustrate the sensitivity of the class-based change analysis towards reg (see Sect. 4.2).

Object changes and interpretation
During the monitoring period, the sliding mass of the right landslide L2 has been reactivated once (in the fourth time step, between 2013/05/23 and 2013/10/13; Fig. 2).This is indicated by a deformation, i.e. lowering of the surface in the upper, erosional area due to erosion or downward movement (negative deformations) and accumulation of material (positive deformations) in the depositional area.This process was likely triggered by intense rainfall in early summer 2013.Contrastingly, these parts of the landslide remained stable during the rest of the monitoring period.In upcoming analyses, time series from a nearby rain gauge will be integrated to assess such relations between detailed process dynamics and potential hydro-meteorological triggers or drivers.
The scarps of both landslides have been subject to retrogressive erosion (negative deformations; Figs. 3 -5), with clods of material (turf and soil) episodically breaking-off and sliding or toppling downward and being deposited in the landslide area (positive deformations).This has resulted in an expansion of the eroded area, predominantly in uphill direction.The deposition of material is less well represented by positive 3D deformations (Figs. 2 (bottom) and 3).This corresponds with the field observation that the deposition of clods of soil is a more diffusive process, which tends to be accompanied by a disintegration of such objects and dispersal on a larger area, thus resulting in surface changes that are often beyond the level of detection (i.e.their absolute value is smaller than the LOD 95 , which in turn depends on the registration error).The ongoing erosion, transport and deposition of material within the landslides prevent the establishment of a continuous grass cover and thus soil formation, at least in the eroded area.
At landslide L1, the main type of change during the monitoring period concerns the objects of the class 'landslide scarp' (Figs. 3,5 (top)).Recurring scarp erosion is contrasted by stability of other object classes.Fig. 4 shows an oblique view and a transect of a subset from the second point cloud epoch (2012/05/02).The mean 3D surface deformation per object (compared to the previous epoch) highlights the erosion (red) of 'landslide scarp' objects and their deposition (blue) a few metres downslope within the landslide.

Level of detection -Parameter sensitivity and impact on the object-based change
The parameter sensitivity was tested using a range of registration error values for the deformation calculation between point cloud epochs.Larger registration errors increase the LOD 95 (Eq. 1) and, accordingly, result in more (per-point) deformations being excluded from class-and object-based aggregation, as they are below the local LOD 95 .
reg The test results reflect this effect, with smaller registration errors resulting in a smaller LOD 95 , hence a lower threshold for acceptance of deformations as "real" change and accordingly higher magnitudes of mean 3D deformations per class for some epochs (Tab. 1 and Fig. 5).For epoch 8, for instance, the conservative estimate reg = 0.10 m resulted in an average level of detection LOD 95mean = 0.200 m (and a standard deviation LOD 95sd = 0.002 m).If the same registration error is assumed, the spatial variability within each point cloud epoch as well as the variability between point clouds is small, at least for the investigated analysis scale, defined by the applied fitting radius D/2 and neighbourhood size for normal vector calculation (low standard deviations; Tab. 1).
Note that the per-object or per-class deformations still can be considerably smaller than the LOD 95 , as each object or class can contain points without change or change below the LOD 95 .This effect, however, is limited to some extent by the integration of the 3D distance as a criterion for object modelling (see Sect. 3.1 and 3.2).In this context, it has to be considered whether these deformations that point to "real" changes, which are of interest for geomorphological interpretation, and how much of the signal must be attributed to data noise or systematic (registration) errors.This highlights the importance of a registration error assessment to estimate the level of detection realistically.Particularly in complex natural environments the conditions for this are often non-optimal, due to a lack of welldistributed, large, planar and stable areas with different surface orientations (cf.Fey & Wichmann, 2016).Recent methodological developments aim at the automated selection of stable areas to use during registration (Wujanz et al, 2016) and the assessment of uncertainties of deformation measurements (Fey & Wichmann, 2017;Kromer et al, 2017).These are important considerations, also with prospect to the processing of longer point cloud time series, where user guided registration as well as visual inspections of (intermediate) results become less feasible.

CONCLUSIONS
This contribution is one of the first works attempting to analyse point cloud time series automatically to monitor landslides at the level of semantically meaningful objects.The presented methods extract information on geomorphological process dynamics from the point clouds and aggregate it at different levels of abstraction.The resulting data products and their visualisations demonstrate the potential for a monitoring of individual objects within a landslide system (parts of the scarp collapsing or clods of material moving downslope) as well as of more generally aggregated change trends for certain object classes (e.g.'scarp', 'eroded area', 'deposit').Reducing the spatial complexity of the point cloud can shift the focus to the temporal dimension to easily identify specific time steps where significant changes occurred (cf.Fig. 3), and then, for these time steps, back to the highly detailed spatial representation of the point cloud (cf.Fig. 4).This kind of semantic 3D scene interpretation will be particularly beneficial for larger datasets (such as time series with many epochs and large, complex point clouds) where the automated analysis workflow and the resulting visualisations assist the interpreter to grasp a quick, yet detailed overview of the most important changes, their timing and character.
The results of our analysis for the presented case study are in accordance with field observations at the test site as well as conceptual landslide models, where retrogressive erosion of the scarp and downslope movement of the sliding mass are major principles of secondary landslide development.Therefore, we conclude that object-based analysis of point cloud time series is a powerful approach to i) break-down the spatio-temporal complexity inherent to many environmental processes and ii) exploit the high level of detail and accuracy of laser scanning.By linking 3D surface deformation measurements with contextual, semantic information, the approach contributes to bridging the gap from real-world surveys to conceptual process models.In future, such methods will assist landslide experts to systematically interpret the measured changes from a geomorphological point of view.In combination with an investigation of drivers (such as rainfall thresholds, slope surface morphometry and soil structure), this will enhance landslide and erosion process understanding and assessment.

Figure 1 .
Figure 1.Terrestrial laser scanning at the test site (top) and a labelled point cloud with two landslides L1 and L2 (bottom).

Figure 2 (
Figure 2 (top).Oblique view of the point cloud coloured by mean 3D deformation per object between 2013/05/23 and 2013/10/13, showing movement of the sliding mass at landslide L2. (bottom) Cumulative mean 3D deformation of the eroded area and the deposit of landslide L2 for the point cloud time series.

Figure 4 .
Figure 4. Scarp erosion and deposition of material between epochs t 1 (2011/10/12) and t 2 (2012/05/02).(a) Oblique view of the landslide L1 coloured by mean deformation per object.(b) Transect coloured by mean deformation per object.(c) Transect overlaying the points of both epochs coloured by their class labels (classes irrelevant for the transect are greyed out in the legend).

Figure 5 .
Figure 5. Cumulative mean 3D deformation of the 'landslide scarp' class for landslides L1 (top) and L2 (bottom), with different registration errors (reg) used for calculating the LOD 95 and, accordingly, filtering "real" change.For some time steps, the class-and object-based deformations calculated with from different reg values, and hence LOD 95, vary systematically.The results obtained with larger registration error values (e.g.reg = 0.1 m) appear to represent only generalized process dynamics.In contrast, the results calculated

Table 1 .
Descriptive statistics per point cloud epoch for the LOD 95 calculated with two different registration errors (reg).