LiDAR And Hyperspectral Data Integration For Landslide Monitoring: The Test Case Of Valoria Landslide

In the framework of the WISELAND project, funded by MIUR, we tested the integration between LiDAR and hyperspectral methodologies in the Valoria landslide (Modena province, Italy), a high risk area with vulnerable elements, subjected to periodic and abrupt reactivations. Multitemporal LiDAR Digital Terrain Models (DTMs) allowed the calculation of a differential surface, highlighting absolute height variations, recognizing the main landslide components and identifying depletion and accumulation zones. Hyperspectral data helped in the landslide terrain roughness characterization, performing the Principal Component Analysis (PCA) and correlating the results with Flatness and Organization geomorphometric parameters derived from LiDAR DTM.


Introduction
Landslide occurrence is related to a variety of factors such as underlying geology, mechanical properties of soil and rocks, degree of weathering, groundwater conditions, and the presence (or absence) of geological structures such as joints, faults, and shear zones [Fell et al., 2000]. Because of this complexity, landslide monitoring is commonly adopted both in the early detection of risk factors and as an effective tool for landslide hazard management and analysis [Sassa e Canuti, 2008]. This paper's aim is to demonstrate the possibility to successful apply high resolution LiDAR and hyperspectral airborne remote sensing techniques to landslide monitoring in the special case of an active, large earthflow characterised by rapid to moderate rate of movement (the Valoria landslide, Northern Apennines, Italy). The airborne platform is preferable instead of the terrestrial one because it allows the acquisition of remote sensing data in large areas like the Valoria landslide. Other works using these methodologies to study landslides can be found on the web, but there are no examples of the integration between differential LiDAR and hyperspectral data in literature. The activities described in this paper are part of the research project WISELAND (Integrated Airborne and Wireless Sensor Network Systems for Landslide Monitoring) funded by the Italian Government (financial years 2007-2009). The Valoria landslide is a large, active earthflow which mostly involves low-plasticity scaly clays [Manzi et al., 2004;Corsini et al., 2006]. It has been completely reactivated in 2001, and since then it has been intermittently active with displacements that in one season could be in the order of hundreds of meters. This recent evolution has caused a significant modification in the slope morphology, with quite distinct depletion and accumulation zones.

Laser scanning and hyperspectral imagery
The possibility to directly acquire a high density and accurate 3D point cloud has made LiDAR (Light Detection and Ranging) the preferred technology for topographic data collection; high-resolution DSMs (Digital Surface Models) and DTMs (Digital Terrain Models), in forestry areas are some example of the potentiality of this methodology [Wehr et al., 1999;Holmgren, 2004;Coren et al., 2006]. A typical LiDAR system consists of a laser ranging and scanning unit, together with a POS (Position and Orientation System), which encompasses an integrated Differential Global Positioning System (DGPS) and an Inertial Navigation System (INS) [Cramer, 1999]. The laser ranging unit measures the distances from the sensor to the mapped surface, while the onboard GPS/INS component provides the position and orientation of the platform. LiDAR data collection is carried out in a strip-wise fashion and the ground coordinates of the laser footprints are derived [Baltsavias, 1999]. The LiDAR we used is an Optech ALTM3100; it is a small footprint LiDAR system that is able to acquire data up to 100 kHz frequency. In spite of very dense and precise spatial data, these systems are rather poor in spectral sensitivity [Coren et al., 2006]. In order to overcome this problem, a hyperspectral dataset has been acquired. Hyperspectral remote sensors, on the other hand, collect image data simultaneously in dozens or hundreds of narrow, adjacent spectral bands. These measurements allow to derive a continuous spectrum for each image cell. After adjustments for sensor, atmospheric, and terrain effects are applied, these image spectra can be compared with field or laboratory reflectance spectra, in order to recognize and map surface materials such as particular types of vegetation or land soils. The hyperspectral system we used is an AISA Eagle system [Hyvärinen, 2003]. It is a hyperspectral sensor allowing the acquisition of a maximum of 252 bands, ranging from visible bands to near-infrared ones. This sensor is the most appropriate to precisely detect many different terrain features. AISA Eagle is a complete, pushbroom system, consisting of a hyperspectral sensor head, a miniature GPS/INS sensor, a data acquisition unit in a rugged PC with display unit and power supply.

Data acquisition
In this study the LiDAR Optech ALTM 3100 sensor was used to investigate the Valoria landslide (Modena Province, Italy). LiDAR datasets have been acquired in 2006, 2007 and 2009. The study area was surveyed from an altitude of 1500m above ground level (agl), with a mean sampling density of about 4 points/m2; the radiometric resolution of LiDAR data, scan frequency and scan width were 12bits, 70Hz and ±25° respectively. The last LiDAR survey was performed on 30th March 2009 with the same sensor and the same acquisition parameters. Hyperspectral data were acquired on 16th June 2009, using the AISA Eagle system. The flight was performed at 3000 m of altitude (agl), acquiring 252 bands and setting a ground resolution of 2 m.

Data pre-processing
All LiDAR datasets were processed using Applanix PosPac software for the trajectory computation. The final point cloud was obtained using Optech DashMap software, while TerraScan software (produced by Terrasolid Corporation) was used for data classification, in order to produce a good high resolution digital terrain model of Valoria landslide [Axelsson, 1999]. Typical vertical accuracy is lower than 0.10 m while horizontal accuracy is in the order of 0.5 m or lower, as reported by Optech Incorporated, [www.optech.ca]. The heights are ellipsoidal. The LiDAR data have been classified generating the ground class by applying the Axelsson algorithm [Axelsson, 1999] embedded in Terrascan software. DTM data have been gridded and interpolated on a regular 2m x 2m grid. This operation was necessary because the geomorphometrical algorithms can't be applied on point cloud and needs gridded data. An error have certainly been introduced, but it has been considered of no importance in this study. The final geocoded hyperspectral dataset was obtained using a self-made software called HSP (Hyper-Spectral Processor), developed by Istituto Nazionale di Oceanografia e di Geofisica Sperimentale -OGS. All the remote sensing datasets are in the following projection: WGS84 ellipsoid, UTM32 North projection.

Geomorphometric analysis
The morphometric analysis of LiDAR derived Digital Terrain Model (DTM) for Valoria orphometric analysis of LiDAR derived Digital Terrain Model (DTM) for Valoria landslide has been performed using MicroDEM software [Guth, 2008]. In this paragraph we focus on some theoretical concepts about morphometric parameters and their computation.
Slope and aspect are calculated in correspondence of every DTM grid point; the vector normal to ground is so defined by applying the steepest neighbor algorithm [Chapman, 1952]. The direction cosines of this normal vector are then calculated. A 3m x 3m matrix is computed, containing the sum of cross products values. Eigenvalues and eigenvectors are extracted from this matrix, normalizing the eigenvectors. Eigenvalues are indicated as S1, S2 e S3; usually S1> S2> S3. The morphometric terrain analysis is usually performed considering these indexes, [Guth, 2003]: 1) Flatness: defined as: Large values indicate flat terrain, low values rugged terrain. It correlates strongly and negatively with slope or relief.
2) Organization: defined as: Large values indicate a dominant linear fabric to the terrain, low values isotropic topography.
3) Orientation: trend of S3. It indicates the dominant trend to the terrain fabric; its direction is between 0 and 180°. It is used in eigenvector analysis of SSO diagrams, [Guth, 2008]. 4) Strength: defined as: Large values indicate flat terrain, low values rugged terrain. It looks very similar to the flatness parameter. It correlates strongly and negatively with slope or relief. 5) Shape: defined as: Large values indicate a dominant linear fabric to the terrain, low values isotropic topography. It correlates moderately with terrain organization. In this study we principally considered Flatness and Terrain organization; the other parameters are supposed to be analyzed in future further investigations.

Principal Component Analysis (PCA)
When dealing with hyperspectral images, as for Valoria landslide analysis, where a large number of useful bands has been acquired, a fundamental task is to perform the so-called Principal Components Analysis (PCA) [Jolliffe, 2002;Coren et al., 2005] generally used to reduce the amount of data to a smaller but significant dataset. In fact, in such images it is very likely that subsets of spectral bands are highly correlated one to each other. If this is the case, you will discover that the accuracy and reliability of a final classification image will suffer if you include highly correlated variables. As a general principle, PCA is a mathematical procedure, often applied in geodesy, transforming a number of (possibly) correlated variables into a (smaller) number of uncorrelated variables, called "principal components". Referring to hyperspectral image processing, the objective of PCA is to reduce the number of bands of the dataset but contemporary to retain most of the original variability in the hyperspectral data. The first principal component accounts for as much of the variability in the data as possible, and each succeeding component accounts for as much of the remaining variability as possible. A PCA is concerned with explaining the variance covariance structure of a high dimensional random vector through a few linear combinations of the original component variables. Considering a p-dimensional random vector: , , ..., X X X X 1 p 1 2 = " 6 6 @ @ The k principal components of X  are the k (univariate) random variables.
This means that the principal components are the linear combinations of the original variables which maximize the variance of the linear combination and which have zero covariance (and hence zero correlation) with the previous principal components. The numerical computation involving a PCA analysis is quite complicated for hyperspectral data, and only some specific software can truly accomplish it; ENVI software is one of the most extensively used for this purpose.

Differential LiDAR DTMs
In the past some photogrammetric digital elevation models have been computed and analysed. For instance, the differential analysis of a DEM of 1973 and of a DTM of 2003 resulted in a clear enough identification of major depletion and accumulation zones occurred after the 2001 reactivation event [Corsini et al. 2007]. However it was impossible to compute volumes precisely, because of an inconstant misalignment between elevation values even in stable zones, due to the stochastic errors relative to accuracies of acquisition techniques.
In this study, LiDAR data from 2006, 2007 and 2009 have been processed, and the relative DTMs have been differentiated. This technique allowed a rather precise quantification of depletion in the source area and of accumulation along the slope and at the landslide toe. LiDAR data bracket in time a quite significant acceleration event occurred in winter 2008-2009. Therefore, a significant picture of slope modification in given by the differential analysis of 2007 and 2009 DTMs (Fig. 1). More specifically, a depletion of about 460.000 m 3 has been estimated for the landslide's head zone. At the same time, the landslide toe has shown a marked bulging, associated to downslope sliding. This has been the result of movements that, on the basis of topographic total station monitoring data, have exceeded 200 m in some slope sectors.

Geomorphometry and PCA integration
The application of morphometric algorithms appeared as a powerful methodology to monitor the Valoria landslide. Flatness and Terrain organization of Valoria landslide DTM were calculated by applying the formulas previously mentioned: see Figure 2 and Figure  3 In Figure 2 we can observe low Flatness values (from 1 to 2) corresponding to zones with a higher differential displacement (see highlighted zones 1 to 5); low Flatness values are associated to rough terrains, and roughness appears to be proportional to stress on landslide surface. Zones highlighted in this figure correspond to high displacement zones in Figure 1.   In Figure 3 we can observe that Terrain Organization values are, instead, quite high in four zones (1 to 4). It's an interesting phenomenon, associated to the presence of a dominant linear fabric; topography generates very clear lineaments on landslide flux directions. Zones highlighted in Figure 3 are not characterized by high slope values (see Fig. 4); high slope values don't correspond therefore to high Terrain Organization values. It means that in Valoria the zones where slope is relatively high are not moving very fast. The PCA calculation has been performed for the Valoria landslide (see Figure 5) considering 210 bands: 42 bands were not taken into account because strongly affected by noise, especially in the Near Infrared field. It produced a new image, totally unlinked to the original one from a spectral point of view; each pixel contains the radiance information of each band, so it is proportional to the original information. In this way, a better discrimination of different terrain surfaces properties can be done, due to the consequent possibility to remove noisy bands. In Figure 5 it is possible to remark a very low PCA values in zones 1-4. It means that the band decorrelations of the hyperspectral data are very strong, probably due to the terrain roughness. Interesting results come from the analysis of the hyperspectral image after applying the PCA algorithm (see Fig. 5). In the four zones previous mentioned (1 to 4), PCA values are very low; it means that terrain roughness strongly affects hyperspectral bands decorrelation. This is a very interesting result and demonstrates that hyperspectral images can find a direct application to landslides monitoring, even if this is to be improved in further studies. Also between Figure 2 and 5 we see a correspondence among zones 1, 2, 3 and 4. Zones characterized by a Flatness value equal to 2 or less than 2 correspond to zones affected by a higher bands decorrelation in PCA image (low PCA values). The same zones are characterized by a high Terrain organization value; it means that the predominant Fabric alignment is clearly marked.
Extending the analysis to Figure 4, we can observe that zones 1, 4 and 5 are also affected by a relevant slope; zones 2, 3 and 6 are instead characterized by a low slope value, even if Flatness and Terrain organization are sensibly relevant. Zone 4 is an accumulation zone and a high slope value is expected; zone 2 and 3 are instead depletion zones, so this result needs to be further investigated. Zone 6 can be differently interpreted depending on the dataset analyzed; a correlation between results coming up from all the datasets (especially from PCA and DTM analysis) doesn't seem to exist. As overall conclusion we can state that morphometric analysis performed jointly with the use of PCA algorithms seems a promising methodology for landslides monitoring. Analysis described in this paper open the access to new research fields. Especially hyperspectral methods are worthwhile to be applied to landslide monitoring. PCA algorithms help identify some key structure in landslide dynamic. Further hyperspectral analysis may try to refine existing geological maps and to identify the spatial distribution of previously unmapped or unknown faults and shear zones through the detection of minerals alteration. Although existing spectral-map libraries can be used to identify minerals; the spectra of a particular mineral can vary depending on the specific host rock; collection of field spectral data will be necessary to ground-truth the remotely sensed data. This analysis hasn't been included in this paper because the work is still in progress and some field measurements are still to be completed; soon we'll have some preliminary results.

Conclusions
The study demonstrates the capabilities of remote sensing techniques to recognise the essential features of an active, rapid earthflow. Using a differential high resolution DTM approach a displacement can be easily detect and the zones with major displacement identified. LiDAR acquisitions in different periods need so to be performed and considered. Areas subjected to a strong landslide activity have been identified by direct DTM analysis. Zones subjected to high relative differential displacements are associated low Flatness values in the DTM analysis; these indicate rough terrains. The same zones are associated a dominant linear fabric; it doesn't seem to be correlated to the local slope obtained by DTM analysis.
Hyperspectral data revealed themselves to be very useful in roughness estimate. PCA analysis is a very useful and powerful methodology to characterize the surface landslide features because of its sensitivity to surface roughness. Future studies will focus on terrain classification by supervised algorithms applications, in order to better identify the landslide lithology, and on the integration of Wireless Sensor Network technologies, in order to complement the existing remote sensing techniques and increase their accuracy [Rosi et al. 2010].