Variogram maps from LiDAR data as fingerprints of surface morphology on scree slopes

Herein, an aerial LiDAR topographic dataset is analysed and interpolated by means of geostatistical techniques in order to examine the morphology of a scree slope area in the Eastern Italian Alps. The LiDAR-derived digital terrain model (DTM) is analysed using variogram maps as spatial continuity indexes. This allowed for evaluation of the reproduction of spatial variability of topography and for the characterization and comparison of different morphological features occurring in the study site. The results indicate that variogram maps efficiently synthesise the spatial variability of topography in a local search window, representing suitable “fingerprints” of surface morphology.


Introduction
Elevation data derived from airborne LiDAR (light detection and ranging, also known as aerial laser scanning) must be carefully analysed and interpolated in order to obtain accurate high-resolution (cell size ≤1 m 2 ) digital terrain models (DTMs).LiDAR data, despite their high spatial density, are affected by errors.In particular, the footprint of the laser beam depends on observation geometry and acquisition distance: as a consequence, accuracy and resolution are variable in the same scene.The footprint, which can be interpreted as support of measurement in geostatistical terms (Isaaks and Srivastava, 1989), is strongly deformed on steep slopes.Moreover, multiple reflections lead to filteringdependent models.Finally, the performance of the inertial navigation system and the accuracy of ground-based GPS network strongly affect data acquisition and modelling.A geostatistical approach makes it possible to account for these Correspondence to: S. Trevisani (s.trevisani@irpi.cnr.it)issues, thereby maximising the informative content of Li-DAR data and filtering out noise.
In the present paper, variogram maps have been used to characterise surface morphology of a scree slope in the Italian Alps.

Study site and LiDAR data set
The study area (Fig. 1a) is located on the southern flank of the Sella Group (Dolomites, Eastern Italian Alps) at elevations ranging from 2550 to 3000 m.Two principal morphological units are present in the study area: subvertical rock slopes, comprised of Dolomia Principale (Norian) and entrenched by steep couloirs, and scree slopes, located at the base of the cliffs.
The present study focuses on the scree slope area.Debris cones, with apexes located at the outlet of couloirs, which cut the superjacent rock slope, are prominent features in the scree belt.Four main processes influence sediment dynamics on the scree slopes and cones, and contribute to their complex morphology, namely, gravitational accumulation of weathered rock fragments, rockfall, snow avalanches, and debris flows (Marchi et al., 2008).The main morphological features of the studied slope are debris-flow channels and debris-flow deposits (at different stages of evolution and activity), rockfall deposits, planar gravitational scree deposits, and isolated Published by Copernicus Publications on behalf of the European Geosciences Union.rocky outcrops emerging from the debris.It is noteworthy that these features are often present as mixed morphologies, for example as debris-flow deposits partly covering rockfall deposits.
The LiDAR data were acquired in October 2006 by means of an ALTM 3100 OPTECH (http://www.optech.ca/prodlatm.htm),flying at an average altitude of 1000 m above ground level under snow free conditions.During the LiDAR survey, high-resolution (0.15 m) digital aerial photos were also collected (camera Rollei H20).The spatial sampling density of the LiDAR-derived topographic dataset is high, with raw data having a mean density of 5 point/m 2 .Raw data were filtered using the Terrascan ™ software classification routines and algorithms (http://www.terrasolid.fi/),and filtered data show a mean density of 3 point/m 2 .

DTM generation and analysis
LiDAR data were analysed and interpolated following a geostatistical approach, which is composed of three steps: 1) explorative and spatial continuity analysis, 2) inference of a spatial continuity model, and 3) interpolation.
The first phase of the analysis is crucial, because it directly affects the choice of the subsequent interpolation method.The main goals of this stage of the analysis are characterisation of sampling geometry, detection of errors and artifacts, delineation of heterogeneities in spatial variability, decomposition between trend and residuals (i.e.decomposition between regional and local spatial variability components), and analysis of spatial continuity.The continuity analysis was performed by means of the experimental variogram (Goovaerts, 1997): where h is the separation vector between two point pairs, z(u α ) is the elevation at the location u α , and N(h) is the number of point pairs separated by h.According to Eq. ( 1), the variogram for a given value of the separation vector h, when applied to topographic data, is half the mean squared differences in elevation between all point pairs separated by the vector h.The analysis of this experimental function for different values of h allows for characterisation of bivariate spatial continuity (or, conversely, the spatial variability) of the studied dataset.After calculation of an experimental variogram, a variogram model must be inferred (phase 2 of the analysis) in order to perform kriging interpolation (phase 3).
A trend model was calculated by employing a local polynomial method, and locally fitting the polynomial surface by using ordinary least squares (search windows with a radius of 15 m).A first-order polynomial was adequate for describing the regional component of variability inside the search windows.The calculated trend surface describes the regional component of spatial variability or, from a signal-processing standpoint, represents the low frequency component of variability.The calculated trend surface was subtracted from the elevation data in order to obtain statistically stationary residuals (Goovaerts, 1997) and estimate an experimental variogram.From the standpoint of signal-processing, the The obtained results are the surface of residuals (Fig. 1b) and the standard deviation of estimation (Fig. 1c), which provide a measure of interpolation quality.Zones with high values of standard deviation indicate higher levels of uncertainty in elevation estimates and are directly related to a deterioration in sampling density.The overall low values of standard deviation indicate a satisfactory quality of the interpolation.Nevertheless, we observe in Fig. 1c peculiar patterns resulting in locally higher values of standard deviation, which could be due to a non-optimal filtering processes of the LiDAR data.Since the focus of this study is to describe morphologies of limited extent, the analysis was performed on the surface of residuals.The complete DTM could readily be obtained by adding the surface of residuals to the trend surface.
The surface of residuals (Fig. 1b), representing the spatial variability at local scale, outlines morphologies with limited width, such as debris-flow channels that are represented by linear patterns with negative values, and debris-flow lobes that are represented by elongated patterns with positive values.The primary morphologies occurring in the study domain are therefore readily identifiable by their specific patterns (Fig. 1b).The quantitative characterisation of morphologies, although more complex than their visual identification on the surface of residuals, is necessary for automated analysis of landforms, including the application of classification and pattern recognition routines (Carr, 1996;Drägut and Blaschke, 2006;Asselen and Seijmonsbergen, 2006;Ehsani and Quiel, 2008).In particular, the use of morphological indexes that are able to catch the specific patterns of the morphologies being investigated is needed.From this perspective, the variogram presented in equation ( 1) represents a good candidate.

Spatial continuity analysis
The variogram, which is a measure of spatial continuity, can be viewed as a multiscale directional measure of surface roughness.For a given value of vector h, defined by a modulus and direction, the variogram is an index of dissimilarity between all point pairs separated by vector h and contained within a given search window.The variogram, for a given direction of h, is plotted as function of the separation distance (i.e. the lag) (Fig. 2b and c).Variogram maps, which visualise the spatial variability along different directions and for different spatial scales in one shot (Fig. 2a), can represent suitable "fingerprints" of local morphology.
Variogram maps were calculated inside circular search windows with a radius of 10 m in correspondence with specific morphologies (Fig. 3).The radius of the moving windows was chosen in accordance with the prevailing size of the investigated morphological features.In order to outline the effect of interpolation in reproduction of spatial variability, the variogram maps were calculated both on raw data (central column of Fig. 3) and on interpolated data (right column of Fig. 3).Given the moving window size, the variogram was calculated using a lag separation of 1 m corresponding to DTM cell size up to a distance of 10 m.The calculation was performed using the open source statistical software R (http://www.r-project.org)installed on a desktop PC with an Intel ® Core ™ 2 Quad 2.4 Ghz processor and 3 Gb of RAM.
In Fig. 3a-c, the calculation is performed on three adjacent morphological features.Starting from the top, the search window covers debris-flow deposits, a scree slope not affected by debris flows, and a debris-flow channel.For the elongated debris-flow deposits (Fig. 3a shows moderate values of variability, with maximum values of approximately 0.08 m 2 .This indicates a strong anisotropy, with the direction of maximum continuity along the direction of deposition of debris flows.A periodicity along the direction of minimum continuity, related to the presence of sub-parallel lobes, is also visible.These results are consistent with the visual recognition of the morphology of debrisflow deposits; the added value of spatial continuity indexes resides in the quantitative and objective assessment of their topographic features.Comparison of the two variogram maps of the original and interpolated data suggests that the interpolated surfaces effectively reproduce the spatial variability, although moderate smoothing tends to slightly obscure the periodicity.The scree slope (Fig. 3b) is characterised by very low values of variability, with a maximum of 0.014 m 2 , an undefined continuity structure, and an absence of anisotropy.In the case of the DTM surface, the variogram map shows a defined continuity induced by the interpolation procedure.The debris-flow channel (Fig. 3c) shows high values of variability, with maximum values of 0.4 m 2 , and strong anisotropy parallel to the direction of maximum continuity along the channel direction.Additionally, there is some periodicity in the variogram along the direction of minimum continuity.The variogram map calculated on the interpolated surface shows a good reproduction of spatial variability.
The analysis was also conducted on two debris-flow deposits, with lobes of different sizes.In the case of larger lobes (Fig. 3e), the presence of periodicity in the direction of minimum continuity is evident for both the variogram computed on the raw data and the interpolated data.Conversely, in the case of the small lobes (Fig. 3d), periodicity is not readily recognisable in the variogram map of raw data (central column), and is absent in the variogram map calculated on the interpolated surface (right column).This can be ascribed to the spatial resolution of the data, which is too low for the description of these small morphological features.
In the case of rockfall deposits (Fig. 3f), the variogram map shows high values of variability, maximum values of 0.4 m 2 with a poorly defined structure of continuity, and a moderate anisotropy at short distances.There are no major differences between the calculation performed on the raw and interpolated data, apart from smoothing in the latter.

Conclusions
Spatial continuity measures, such as the variogram, synthesise the spatial variability structure of morphological features and provide informative tools for the characterisation of surface roughness.Spatial continuity indexes highlight and characterise terrain morphologic features, such as debrisflow deposits, channels, and rockfall deposits, resulting in "fingerprints" of surface morphology.Variogram maps can be used to derive and calibrate simpler morphological indexes that are capable of picking up the most distinctive characteristics of the spatial variability of different morphological features, such as the anisotropy, the presence of periodicity, and the spatial variability calculated at specific lags.These indexes could serve as the basis for subsequent automated classifications and pattern recognitions (Carr, 1996;Asselen and Seijmonsbergen, 2006;Drägut and Blaschke, 2006).
Although most of the morphologies considered in this study could be identified through visual mapping, the possibility of quantifying land surface features by means of spatial continuity indexes is of great importance, as allows for objective comparisons between morphologies.Moreover, the quantification of spatial variability could be useful to evaluate the relationships between spatial patterns and the processes involved.Finally, these indexes may enhance the mapping of morphological features not readily identifiable through traditional techniques, such as field observations or aerial photo interpretation, or other DTM-derived maps (e.g.shaded relief maps).
The present interpolation procedure, based on the block ordinary kriging algorithm, facilitated an accurate surface model of the residuals to be calculated, with a satisfactory representation of the spatial variability.In case of lower data densities or the need for higher resolution, a local kriging approach (Stroet and Snepvangers, 2005) could further improve the quality of the DTM.In particular, the correct representation of spatial variability of ground surface in presence of anisotropic morphologies should be closely considered, as represented in the study area by debris-flow deposits or channels.

1Figure 1
Figure 1.a) study site location with domain of analysis outlined in red; b) surface of 2 interpolated residuals; c) map of standard deviation of estimation.3 3

Fig. 1 .
Fig. 1.(a) study site location with domain of analysis outlined in red; (b) surface of interpolated residuals; (c) map of standard deviation of estimation.

Figure 2 .Fig. 2 .
Figure 2. Relationship between the variogram map and directional variograms.a) variogram 7 map of elevation data with profiles of directional variograms; b) directional variogram 8 calculated along the direction of maximum continuity, (dashed line in a); c) directional 9 variogram calculated along the direction of minimum continuity (continuous line in a). 10 11

Figure 3 .Fig. 3 .
Figure 3. Variogram maps for different morphological features calculated on raw data and on 2 interpolated data.The red circles represent the search window for calculations.a) debris flow 3 Fig. 3. Variogram maps for different morphological features calculated on raw data and on interpolated data.The red circles represent the search window for calculations.(a) debris-flow deposits; (b) scree slope; (c) debris-flow channel; (d) and (e) debris-flow lobes of different size; (f) rockfall deposits.