key: cord-0059056-n2mgp5uq
authors: Lazzari, Maurizio
title: High-Resolution LiDAR-Derived DEMs in Hydrografic Network Extraction and Short-Time Landscape Changes
date: 2020-08-19
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58802-1_52
sha: 6517c486cf2e2db2205a1d68c8f0c1af60596f43
doc_id: 59056
cord_uid: n2mgp5uq

In this paper an automatic methodology to extract the channel network from high-resolution LiDAR-derived DTMs and a semi-quantitative methodology to assess the short-time landscape evolution of a test-area, located in southern Italy, have been applied. In particular, the technique used is based on a local nonlinear filter together with the global geodesic optimization for channel head and drainage network extraction. Further, the two Lidar acquisition for the year 2012 and 2013 have been used to detect hydrographic network changes and slope evolution in terms of erosion and deposition pattern and then compare them with the slope processes (landslides and linear erosion).

During last twenty years, new methods of acquisition of high-resolution, high-quality topographic data have provided, in unprecedented manner, the opportunity to investigate landscapes and landforms at the scale of geomorphological process [1] . One of the most significant advances is airborne laser swath mapping (ALSM) or LiDAR that nowadays gives the way to a new generation of high-resolution (0.5 to 2 m) digital terrain models (DTMs) offering new possibilities to use detailed representations of surface features [2, 3] .

Thanks to the availability of these non-conventional datasets, researchers have the ability to analyze landforms at a scale ranging from few meters to decimetres which is crucial to investigate slope and fluvial processes.

This gives new insights for the study of the "geomorphological signatures" of surface processes and their significance in terms of landscape evolution.

Recent application of high resolution topography concentrated on identification and extraction of geomorphological features [4] [5] [6] , landslides morphology and activity detection [7] [8] [9] [10] , river channel morphology [11] , tectonic geomorphology [12] , geoarcheology [13] , landscape change [14] .

High resolution topography also favourites the development of new and advanced methods of terrain analysis in particular for channel network extraction. New methods based on landform curvature have been developed [2, [15] [16] [17] [18] .

The proposed methods detect thresholds in topographic curvature for channel network identification and hillslope channel transition, thus providing an alternative to classical methodologies for channel network extraction such has, for example, a constant threshold contributing area [19] [20] [21] , a slope dependent threshold contributing area [22, 23] , and also a threshold on local curvature [24, 25] .

Few studies have explored the potentiality of Lidar application to landform change detection and landscape evolution [14] , mainly due to the high costs for data acquisition. Further the recent availability of the technology, such as UAV survey, constrains often a too short timescale to detect changes in the landscape. On other hand multitemporal acquisition of Lidar data could be very useful to define spatial patterns of erosion and deposition and to estimate sediment flux within the landscape and the consequent landscape evolution. This constitutes a future challenge in geomorphology to develop new geomorphic transport law [26] .

In this paper a case study located in Basilicata region (southern Italy), whose peculiarity is to be a rapid evolving landscape, is proposed. The climate regime, the sparse vegetation and the high erodible lithology result in a deep dissected landscape with the presence of badlands and landslides on the steepest slopes [27] .

A semi-quantitative methodology to assess landscape evolution in this area on the base of multitemporal Lidar data analyses has been applied. In particular, it has been considered the technique proposed by [16, 17] , based on a local nonlinear filter together with the global geodesic optimization for channel head and drainage network extraction. Further two Lidar acquisition for the year 2012 and 2013 have been used to detect landscape changes in term of erosion and deposition pattern and then compare them with the actual hillslope channel configuration.

The study was carried out in a test site located in the central sector of the Basilicata region, southern Italy ( Fig. 1) , about 2 km east of the Calciano town on the right side of Basento river.

The area is located along the outer thrust-front of the Lucanian Apennine Chain, where the turbiditic arenaceous-marly Miocene formations of Serra Palazzo outcrop. Above the Miocene substratum there is the regressive succession of Bradanic Foredeep, in stratigraphic contact with angular discordance, represented by Argille Supappennine Fm. (lower Pleistocene), characterized by silty-clayey and sandy-clayey deposits [28] , laminates and with bioturbation, reportable to a shallow water sea of Sicilian in age (small Gephyrocapsa Zone; [29, 30] ). Toward the top this formation gradually passes, in continuity of sedimentation, to the regressive sandy-conglomerate, reportable to depositional delta, coastal systems and continental systems (plain alluvial).

The interaction of the regional uplift with other local driving factors, such as the relative changes of the sea level, sinsedimentary and post-depositional tectonics, action of the erosion, has contributed to define the actual geological and geomorphologic setting of the study area.

The dominant clayey lithological nature and the north-east exposure of the study area favour conditions the development of mudflow, shallow landslides triggering and deep retrogressive channel erosion (gullies), that determinates also high sedimentation rates at the base of the slopes.

The acquisition of laser scanning data was carried out by CNR and Aerosigma srl on May 2012 and May 2013 (Fig. 2) , using a full-waveform scanner, RIEGL Q680i on board twin-engine plane P68 victor B -Partenavia to obtain a higher spatial resolution, with a Ground Sampling Distance (GSD) of about 5 cm with overlap of 70% and sidelap of 30% ( Table 1) .

The processing of laser point cloud was made using the commercial Terrasolid software, one of the most complete and efficient software for laser points processing. In particular, TerraScan, TerraPhoto and TerraModel by Terrasolid package has been used [31] .

The DTMs of different years (2012 and 2013) of the study area were first standardized (precision, accuracy, reference systems) and then analyzed, area by area, to make differences among multitemporal DTMs (Fig. 3) .

To determine high-resolution DTM from laser scanning point clouds filtering and classification processing was performed. The extraction of unclassified.las data was followed by the post-processing phase of the point cloud, i.e. the processing chain aimed at their classification. Every peculiar feature of the scene, represented by points that satisfy a certain property, is inserted in special classes (ground, vegetation, buildings, electricity lines, roads), which include points with similar characteristics.

In the case study, the classification procedures were carried out with the help of the Terrasolid modules, specifically the Terrascan module, the most common, as well as the most robust and performing commercial classification application.

Initially, an analysis of the entire territory was made to be processed in order to apply the correct "ground" classification parameters.

The cleaning of the point cloud is a fundamental step in the classification phase as it allows you to "move" spurious points (the presence of which is inevitable) in a class that is not taken into consideration in subsequent elaborations. This discrimination was carried out through a terrascan macro carried out in two steps:

• Elimination of low spurious points, or better known as Low Points • Elimination of high points The low points have been classified by inserting a reference altitude which is far lower with the minimum altitude of the acquisition scene (>0.5 m).

The ground points were interpolated with a Delaunay triangulation, the most popular TIN interpolation. TIN is generally able to represent abrupt changes in the land surface, because the density of the triangle can be varied easily and so it is able to incorporate discontinuities. TIN-grids have been interpolated to derive the cell-based DEMs. Each raster cell contains a size of height determined before resembling, and 

LiDAR derived DTM permits direct detection of channels. New algorithms have been proposed for the automatic extraction of channel network from high resolution DTM obtained from airborne LiDAR data, such as [15] and [16] . Both of them are based on the local slope and landform curvature concepts, whilst traditional methodologies are based on drainage area and local slope thresholds. Tarolli and Dalla Fontana [2] showed that the curvature derived from high-resolution topography using a Lidarderived DEM, seems to improve channel network extraction because it permits to recognize in detail concave-convergent terrain associated with fluvial erosion processes.

According to these considerations, in order to extract the channel network from high-resolution LiDAR derived DTM of Calciano test-area, an automatic methodology has been applied. Between the two algorithms proposed by [15] and [16] it has been used the second one, i.e. [15] according the results discussed in detail in [17] . It compares the capability of [15] and [16] to extract the channel network, capturing Passalacqua et al. [16] script proposed to the community a Matlab toolbox named GeoNet, where the accumulation area was computed outside the code and the curvature threshold was visually determined from quantile-quantile plot. The choice of the parameters is automatic and the accumulation area is computed inside the code.

Before to run DTM image processing GeoNet considers DTM (as a geotiff format): to each option a tiff tag is assigned. These tags are used to store geotiff information of any type.

Passalacqua et al. [16] starts by applying a smoothing filter on the original data-set for removing noise and irregularities which can reduce the detection of the edges of features at the scale of interest and identify features as entities that persist over a range of scale.

GeoNet applies the non linear filter performed by [32] . From a numerical point of view it exactly used the Perona and Malik filter proposed by [33] .

The choice of using nonlinear filter as smoothing filter is especially motivated by the space-scale paradigm that the Perona and Malik [32] formulated to address the degradation of spatial localization of natural boundaries, especially at larger scales of smoothing, due to the isotropic criteria which characterize other smoothing filters (see for example Gaussian filtering used by [15] to extract channel network).

The isotropic nature of many smoothing filters does not allow to respect the natural boundaries of the features, which may represent important discontinuities in the landscape.

In particular, the Perona and Malik [32] filter method consider these three criteria:

1) causality: a scale-space representation should have the property that no spurious feature should be generated passing from finer to coarser resolutions, as stated by [34] and [35] ; 2) immediate localization: the region boundaries should be sharp and meaningful at each resolution; 3) piecewise smoothing: at all resolutions the smoothing should be preferentially intraregion that interregion.

If h 0 (x,y) = h(x,y;t) represents the original high resolution digital elevation model the nonlinear filtering operation can be expressed as: @ t hðx; y; tÞ ¼ r Á cðx; y; tÞrh ½ ¼cðx; y; tÞr 2 h þ rc Á rh ð1Þ

where t indicates that the derivative is taken in time, c is the diffusion coefficient dependent of the space, r is the gradient operator. The best situation to choose the diffusion coefficient c is in the way that the nonlinear filter achieves noise reduction and emphasizes edges, without smoothing across boundaries, i.e. the situation in which the channel boundaries location is known and so c = 0 at the boundary and c = 1 everywhere else.

The channel location is not known in advance but only an estimate of it or the knowledge of some geometric characteristics may be computed. According [32] a first estimate of the edgesẼðx; y; tÞ location within the non linear space-scale paradigm can be given by the gradient of the elevation h(x,y,t) at the location (x,y) and time t: Eðx; y; tÞ ¼ rhðx; y; tÞ ð 2Þ

being the diffusion coefficient c a function of the magnitude of the gradient c ¼ pð rh j jÞ where pð rh j jÞ is the edge stopping function able to avoid diffusion across boundaries and computed using the following formula: pð rh j jÞ ¼ 1

where k is a constant estimated from the 90% quantile of the probability distribution function of the absolute values of the gradient throughout the digital elevation model image [32] . Such Eq. (1) can be written as:

which is the nonlinear filter equation used in [16] . GeoNet allows to select as a parameter of the algorithm the N number of nonlinear filtering iterations, the number of steps needed to achieve noise reduction and discontinuities enhancement before proceeding with channel extraction. This value must be defined on the base of the scale of the objects it wants removing from the digital elevation model data. In particular, it has been tested as the skeleton appear for different value of N in the original data and finally the best number filter iterations for our case study is N = 50.

GeoNet by [16] proceeds creating a skeleton of the likely channelized pixels using the quantile-quantile plot of the curvature (as in the work by [15] ) and then to further narrow the skeleton by introducing a threshold of the contributing area, the so called area threshold.

Skeleton is a binary matrix where pixels that satisfy the threshold criteria (area threshold and geometric curvature) have a value of 1, while pixels that do not satisfy them have a value of 0.

For the quantile-quantile plot of the curvature [16] used the definition of geometric curvature of the isoheight contours, defined as

and computed by standard finite differences.

As concerns the quantile-quantile plot of the curvature, it could be possible to use also the Laplacian curvature (as employed in [15] ). A discussion about the advantages to obtain the skeleton of the pixels using geometric curvature is in [16] . It is a value which depends on the landscape in analysis; its selection must be small enough to not interfere with channel initiation and large enough to effectively reduce the presence of isolated convergent areas which appear as noise in the skeleton [16, 17] .

In this study the value of area threshold is of 4000 m 2 and threshold quantilequantile curvature = 1.

After the identification of the likely channelized pixels (i.e. the creation of the skeleton using the curvature and the contributing area), GeoNet by [16] performed the extraction of the channels using the geodesic curves. It is defined as curves of minimal cost among all the possible curves C connecting point a (starting point) and point b (ending point):

where WðCÞ : X !R þ represents the cost of traveling on a curve C 2 X:

The cost function W includes surface curvature and flow accumulation because they are the topographic attributes that distinguish channels from the rest of the landscape. In particular, curves with the positive curvature above threshold defined from the quantile-quantile plot and large value of flow accumulation are preferred to all the possible curves connecting two points.

Formally the geodesic curve is computed by gradient descendent on the distance function, backtracking from the downstream point b. So, the geodesic is the integral curve of rd starting at point b and the gradient is computed on the surface.

When the Skeleton obtained by the thresholding curvature and contributing area is computed and the channel extracted, the algorithms proceeds with the detection of the end points. They are identified as the points at which the skeleton end.

Because the skeleton is wider than one pixel, the point taken as end point is located at the end of the skeleton and at minimum geodesic distance from the outlet. So to define the end points a cost function is chosen to penalize the paths along which the drainage area does not have large flow accumulation and along which the curvature is not large compared to the surrounding points. The cost function is:

where, A is the contributing area, K is the geometric curvature and a e d are two constants. The contributing area was computed using the D inf algorithm [36] . The choice of the a e d constants of the cost function is based on the optimal computation of the geodesic distance. The purpose of these constants is to take care of the dimensionality of w (as A is measured in m 2 , while k in 1/m) and of the difference in the order of magnitude between the quantities employed (A varies between 1 and 5 Â 10 5 m 2 , while k has been normalized and thus varies between 0 and 1).

Other parameters used are: DEM smoothing quantile = 0.9; Point search box size = 30.

Following the previously illustrated methodology a processing of the data of the two DEMs was carried out obtaining the curvature, using Perona-Malik filtered data, and skeleton obtained by thresholding curvature (Fig. 4) .

From the skeleton of Fig. 4 , it can detect the river network outlet, as the point with the maximum flow accumulation area, computed, for example, using the Dinf algorithm [36] . After the outlet of the network has been identified, it can proceed with the detection of the end points applying the cost function. These are identified as the points at which the branches end. Since the channels are wider than one pixel, the actual point taken as end point is the one which belongs to the minimum geodesic distance path. The river network has been extracted and the geostructure for channel heads and drainage paths created (Fig. 5) .

Among the advantages that the analysis of multidemporal LIDARs by DEMs offers there is certainly the possibility of verifying the evolution of the hydrographic network in a short time-span, through the new positions reached by the channels head points.

The comparison of the two scenarios of 2012 and 2013 (Fig. 6) shows, in fact, that the minor hydrographic network is rapidly evolving with retrogressive channels, which tend to achieve an equilibrium profile. Nevertheless, along this north-facing slope, it is quickly altered and changed by the presence of several shallow landslides.

Among other applications, the subtraction of DEMs (Fig. 7) with spatial analysis using Map algebra functions in ArcGIS or QGIS, allows to identify the fast-changing areas of the slope for the combined action of erosion and deposition, exercised, respectively, by the hydrographic network and landslides [37] . 

In this paper an automatic methodology, based on Passalacqua et al. script [16] , to extract the channel network from high-resolution LiDAR-derived DTMs and a semiquantitative methodology to assess the short-time landscape evolution of a test-area, located in Basilicata (southern Italy), have been applied.

In particular, the technique used is based on a local nonlinear filter together with the global geodesic optimization for channel head and drainage network extraction. Further, the two Lidar acquisition for the year 2012 and 2013 have been used to detect hydrographic network changes and slope evolution in terms of erosion and deposition pattern.

High-resolution DEMs offer new opportunities for extracting detailed features from landscapes (e.g., channels, erosional/depositional areas, channel heads), but also challenges in hydrographic network developing.

Availability of this data provides new opportunities to study the spatial organization of landscapes and channel network features and the comparison of high-resolution DEMs for the investigation of short-term evolution of complex areas largely affected by widespread fluvial and slope processes, increasing the accuracy of environmental transport models, and inform decisions for targeting conservation practices.

Automatic extraction of detailed and localized geomorphic features of interest, such as channel networks, is very important for accurate estimation of sediment sources and flux transport, especially along slopes, where morphoevolutive mechanisms are complex and difficult to interpret using only the geomorphological field survey.

The search for a topographic signature of life

Hillslope to valley transition morphology: new opportunities from high resolution DTMs

Understanding earth surface processes from remotely sensed digital terrain models

Geomorphological mapping of glacial landforms from remotely sensed data: an evaluation of the principal data sources and an assessment of their quality

Characterizing arid region alluvial fan surface roughness with airborne laser swath mapping digital topographic data

Tectonic geomorphology of the

Objective landslide detection and surface morphology mapping using high-resolution airborne laser altimetry

Analysis of LiDAR-derived topography information for characterizing and differentiating landslide morphology and activity

Identification and mapping of recent rainfall-induced landslides using elevation data collected by airborne LiDAR

Automated landslide mapping using spectral analysis and high resolution topographic data

The effectiveness of airborne LiDAR data in the recognition of channel bed morphology

Geomorphic response to uplift along the Dragon's zone from high resolution topography: an example from the Cholame segment, Geomorphology; special issue on high resolution topography

Remote sensing and LiDAR applications in the alluvial geoarchaeology of NE Italy

Conversion of steep Hawaiian hillslope transport from soil creep to overland flow accelerates erosion rates over 100-fold

Channel network extraction from high resolution topography using wavelets

A geometric framework for channel network extraction from LiDAR: nonlinear diffusion and geodesic paths

Testing space-scale methodologies for automatic geomorphic feature extraction from LiDAR in a complex mountainous landscape

Automatic geomorphic feature extraction from lidar in flat and engineered landscapes

The extraction of drainage networks from digital elevation models

The analysis of river basins and channel networks using digital terrain data

On the extraction of channel networks from digital elevation data

Channel initiation and the problem of landscape scale

Analysis of erosion thresholds, channel networks and landscape morphologyusing a digital terrain model

Fractal River Basins. Chance and Self-Organization, 528 p

Development and comparison of approaches for automated mapping of stream channel networks

Geomorphic transport laws for predicting landscape form and dynamics

Landslide inventory of the Basilicata region (Southern Italy)

Note illustrative della Carta Geologica d'Italia: F°188

Modello stratigrafico-deposizionale della successione regressiva infrapleistocenica della Fossa Bradanica nell'area compresa tra Lavello

Il comportamento tettonico e sedimentario del bacino d'avanfossa Bradanica durante il Pleistocene inferiore. Volume in memoria di ALFREDO JACOBACCI "Evoluzione delle conoscenze geologiche dell'Appennino Apulo-Campano e Tosco-Umbro-Marchigiano

TerraScan User's Guide. Terrasolid

Scale-space and edge detection using anisotropic diffusion

Image selective smoothing and edge detection by non linear diffusion

Scale-space filtering

The structure of images

A new method for the determination of flow directions and contributing areas in grid digital elevation models

Tracking landslide displacements by multi-temporal DTMs: a combined aerial stereophotogrammetric and LIDAR approach in western Belgium