Assessment of the Evolution of a Landslide Using Digital Photogrammetry and LiDAR Techniques in the Alpujarras Region ( Granada , Southeastern Spain )

In this work a detailed analysis of the temporal evolution of the Almegíjar landslide is presented. It is a rock slide located in the Alpujarras region (Granada, Spain) that has developed over the last 30 years. Six datasets and photogrammetric flights corresponding to the years 1956, 1984, 1992, 2001, 2008, and 2010 were surveyed. The more recent flight of 2010 combined an aerial digital camera and a LiDAR sensor and was oriented by means of in-flight data and tie points. This 2010 flight allowed for the generation of a reliable and high-precision Digital Terrain Model (DTM). The other flights were oriented using second-order ground control points transferred from the 2010 flight, and the corresponding DTMs were prepared by automatic matching and subsequent editing from the stereoscopic models. After comparing the DTMs of different dates, it has been observed that the landslide was triggered after 1984 and since then has evolved in an irregular pattern with periods of variable activity. On average, the ground surface dropped more than 8 m in depleted zones and rose nearly 4 m in the accumulation zones, with a velocity catalogued as very slow (about 15–30 cm/year) over a time span corresponding to a degree VIII of diachroneity. The total volume of the mobilized mass of this large contemporary slide was about 300 × 103 m3.


Introduction
This article is an extended and improved version of the conference paper titled "Digital photogrammetry and LiDAR techniques to study the evolution of a landslide" [1], presented at the Eighth International Conference on Geo-information for Disaster Management (Gi4DM), celebrated in Enschede (The Netherlands) in December 2012.
The application of remote sensing techniques to natural hazards and landslide research has been steadily expanding in the last two decades [2,3], with different approaches from the optical spectrum to Synthetic Aperture Radar (SAR) techniques [4,5].Among these techniques, those based on photogrammetry and Light Detection and Ranging (LiDAR) are adequate for middle-to high-resolution studies.
Most studies based on photogrammetric techniques follow a similar methodology, starting with the orientation of the digital aerial images using photo-triangulation techniques by means of bundle block adjustment [1,7,10,11,[13][14][15][16][19][20][21][22][23]26,27].Then it is possible to use a reduced number of ground control points (GCPs), which are generally established by means of GNSS techniques, discarding in each survey any control points that are not adequately identified in the corresponding images.
The quality or positional accuracy of these analyses is easily assessed in most of the aforementioned studies by means of the root mean square error (RMS) at the residuals of the control and/or check points [19,22,29,32,33,37,39,40], including the error propagation to DEMs and their differentials [7,23,38].
The methodology of this study starts from a combined camera and LiDAR flight oriented with in-flight information (onboard GNSS and inertial data), which provides a common reference system for both the photogrammetric and the high-resolution LiDAR digital models.Taking this flight as a reference, the other photogrammetric flights were oriented by means of second-order ground control points that were transferred from the reference flight.Once these historical flights were oriented, the DTMs corresponding to each survey were obtained by outdating the reference DSM (from LiDAR data) by means of stereoscopic editing in a photogrammetric workstation.Then differential DTMs and subsequent profiles and volumetric changes were computed.It is a rather efficient methodology that considerably reduces the number of field-surveyed ground control points and allows for the geo-referencing of all data in a common reference system.

Geographical and Geological Setting
The study zone is located in the Alpujarras region, a mountainous area in the Sierra Nevada range (Granada province, Southern Spain), near the Mediterranean Sea (Figure 1).The region is widely affected by slope instability processes, giving rise to abundant rock falls, slides, and debris flows resulting from a combination of abrupt relief and metamorphic geological units of the Internal Zone of the Betic Cordillera, rather susceptible to landslides [53][54][55][56][57][58][59].The main triggering factor in the region is the combination of torrential rainfalls, the morphological slope, and a river channel evolution controlled by semiarid environmental conditions which gives place to a "rambla" segment in the Guadalfeo River.This channel has a planar bottom in which ephemeral streams alternate with variable flooding controlled by extraordinary intense storms or the quick melting of Sierra Nevada snow during spring rising of temperatures.
The landslide studied is a good example of the many landslides resulting from the slope evolution on the intensely eroded northern margin of the Guadalfeo River [59].It is settled near the village of Almegíjar (Granada, Spain), but clearly outside the urban area and affecting only terrain free of any buildings or infrastructure.Given the fresh morphological features and the sliding material, it has been classified as a translational rock slide, with a clear upper scarp and a planar crown of 20 m height.The landslide mass is 250 m long, 300 m wide, and 140 m high (Figure 2).Also, some secondary minor scarps are visible, besides rock fall and debris flow in the frontal part of the mass [1,60].The main triggering factor in the region is the combination of torrential rainfalls, the morphological slope, and a river channel evolution controlled by semiarid environmental conditions which gives place to a "rambla" segment in the Guadalfeo River.This channel has a planar bottom in which ephemeral streams alternate with variable flooding controlled by extraordinary intense storms or the quick melting of Sierra Nevada snow during spring rising of temperatures.
The landslide studied is a good example of the many landslides resulting from the slope evolution on the intensely eroded northern margin of the Guadalfeo River [59].It is settled near the village of Almegíjar (Granada, Spain), but clearly outside the urban area and affecting only terrain free of any buildings or infrastructure.Given the fresh morphological features and the sliding material, it has been classified as a translational rock slide, with a clear upper scarp and a planar crown of 20 m height.The landslide mass is 250 m long, 300 m wide, and 140 m high (Figure 2).Also, some secondary minor scarps are visible, besides rock fall and debris flow in the frontal part of the mass [1,60].The main triggering factor in the region is the combination of torrential rainfalls, the morphological slope, and a river channel evolution controlled by semiarid environmental conditions which gives place to a "rambla" segment in the Guadalfeo River.This channel has a planar bottom in which ephemeral streams alternate with variable flooding controlled by extraordinary intense storms or the quick melting of Sierra Nevada snow during spring rising of temperatures.
The landslide studied is a good example of the many landslides resulting from the slope evolution on the intensely eroded northern margin of the Guadalfeo River [59].It is settled near the village of Almegíjar (Granada, Spain), but clearly outside the urban area and affecting only terrain free of any buildings or infrastructure.Given the fresh morphological features and the sliding material, it has been classified as a translational rock slide, with a clear upper scarp and a planar crown of 20 m height.The landslide mass is 250 m long, 300 m wide, and 140 m high (Figure 2).Also, some secondary minor scarps are visible, besides rock fall and debris flow in the frontal part of the mass [1,60].The landslide is considered to be active, although alternating periods of very slow speed with more active periods during ephemeral floods of the "rambla" [53,58,59].Also, according to [60], the slide reactivation in the period 2009-2010 showed a bulging in the lower sector of the displaced mass, following a mechanism described by [61].
The mass of the Almegíjar slide is composed of an Alpujarride unit of quartzites and quartzphyllites of Permian to Triassic age [62].Nevertheless, the slide is part of an unstable but dormant slope showing an upward extension of the instability expressed by the appearance of large rock blocks and open cracks affecting the overlying marble unit, with water flows in its bottom infiltrating the mass of phyllite rock.

Materials and Methods
Digital photogrammetry and LiDAR techniques represent a significant advance in the study of ground surface evolution derived from landslides and other processes.The results using photogrammetric workstations have offered better and more accurate interpretations of changes on the terrain surface over the past decade [8,[13][14][15][16].
Thus, landslide inventories and their geomorphic features (scarps, crown, lateral flanks, cracks, etc.) may be prepared by means of 3D digital stereoplotting.In this sense, not only single landslides but also the identification of unstable slopes in a region may be analyzed with these techniques, starting from high-resolution DEMs and 3D inventories [22,44,47].In addition, if data from different surveys are available, it becomes possible to study landslide evolution.The overlaying of landslide inventories in Geographical Information Systems (GIS) allows the study of landslide evolution in a given region and the identification of reactivated landslide areas as well as the growth or generation of new landslides.Even landslide multitemporal landslides can be obtained from accurate DEMs of different epochs [44,47].Thus it is possible to determine landslide activity, differentiating between active, dormant, and relict landslides.The concept of landslide diachroneity [63] distinguishes landslides considering the time span of the active periods and the initial triggering data.
The starting point is a dataset made up by photogrammetric flights and LiDAR data.The methodology, established by the research team based on previous results [1,19], includes a set of tasks summarized in the following steps (Figure 3): Dataset compilation and digitization of images from historical flights.2. Definition of reference system and orientation of reference data.
Building and editing DTMs.

5.
Comparison between models and calculations.

Dataset Compilation and Digitization of Images from Historical Flights
First, a photogrammetric flight corresponding to the year 2010 was available, combined with a digital camera (Z/I DMC) and a LiDAR sensor (Leica ALS50-II), as well as attached GNSS/IMU systems for direct orientation.The Z/I digital camera offers a resolution of 0.20 m and has four spectral bands, three in the visible spectrum (RGB) and one in the near-infrared (NIR).The LiDAR instrument produces a ground resolution of about 1-1.5 points per m2 .A second photogrammetric flight of 2008 with the same properties of the first one was also available but without LiDAR dataset.
Besides the highly accurate and recent flights, historical images from aerial film cameras were also available, corresponding to four photogrammetric flights of the years 1956, 1984, 1992, and 2001.The first one was the so-called "American flight", as it was undertaken by the U.S. government in 1956, and it was panchromatic at a middle scale (1:33,000).The second one (1984) was a panchromatic flight made by the National Geographic Institute of Spain (IGN) at scale of 1:30,000.The other two (1992 and 2001) were also panchromatic flights at a scale of 1:20,000, made by the Andalusian Cartographic and Statistical Institute (IECA).The main flight features are indicated in Table 1 and their distributions are shown in Figure 4.

Dataset Compilation and Digitization of Images from Historical Flights
First, a photogrammetric flight corresponding to the year 2010 was available, combined with a digital camera (Z/I DMC) and a LiDAR sensor (Leica ALS50-II), as well as attached GNSS/IMU systems for direct orientation.The Z/I digital camera offers a resolution of 0.20 m and has four spectral bands, three in the visible spectrum (RGB) and one in the near-infrared (NIR).The LiDAR instrument produces a ground resolution of about 1-1.5 points per m 2 .A second photogrammetric flight of 2008 with the same properties of the first one was also available but without LiDAR dataset.
Besides the highly accurate and recent flights, historical images from aerial film cameras were also available, corresponding to four photogrammetric flights of the years 1956, 1984, 1992, and 2001.The first one was the so-called "American flight", as it was undertaken by the U.S. government in 1956, and it was panchromatic at a middle scale (1:33,000).The second one (1984) was a panchromatic flight made by the National Geographic Institute of Spain (IGN) at scale of 1:30,000.The other two (1992 and 2001) were also panchromatic flights at a scale of 1:20,000, made by the Andalusian Cartographic and Statistical Institute (IECA).The main flight features are indicated in Table 1 and their distributions are shown in Figure 4.The period considered (c.50 years) is long enough to study short-term changes in the relief and, more precisely, to analyze the landslide evolution in the study zone.

Definition of Reference System and Orientation of Reference Data
As a starting point for the entire process, the 2010 flight (camera and LiDAR) was used.This flight was oriented by means of in-flight orientation systems (direct orientation).As this dataset was considered to be of the highest quality, both in image resolution and point cloud accuracy, it has been proposed as the reference data and the previous flights have been registered with respect to it.
Thus, it is important to ensure the adjustment and minimal data discrepancies between the photogrammetric flight and the LiDAR point cloud of this 2010 reference flight.The discrepancies between them are originated as a consequence of the independent orientation systems of these sensors-despite the simultaneous recording of data-and thus its different positioning, which may prompt a lack of geometrical coincidence in the data compiled.This raises the need for an optimal adjustment between the two types of data [64].
In accordance with the above, the reference data are compiled and refined by three steps: 1. Photogrammetric orientation using in-flight GNSS and inertial data to ensure the geometrical quality of the whole photogrammetric block.This procedure was performed by means of a photogrammetric workstation and Socet Set 5.2 software (BAE Systems Plc., London, UK) [65] by measuring about 200 tie points and readjusting the orientation parameters.The RMS errors of the resulting model, based on the residual of 20 check points surveyed with GNSS (GNSSbased check points), are 0.258 m in XY and 0.100 m in Z (Table 2).2. LiDAR data orientation, by readjustment of different strips, which is employed to obtain a single homogeneous point cloud, even in overlapping zones.This task was carried out with the TerraMatch module of TerraSolid software (Terrasolid Ltd., Helsinki, Finland) [66].The period considered (c.50 years) is long enough to study short-term changes in the relief and, more precisely, to analyze the landslide evolution in the study zone.

Definition of Reference System and Orientation of Reference Data
As a starting point for the entire process, the 2010 flight (camera and LiDAR) was used.This flight was oriented by means of in-flight orientation systems (direct orientation).As this dataset was considered to be of the highest quality, both in image resolution and point cloud accuracy, it has been proposed as the reference data and the previous flights have been registered with respect to it.
Thus, it is important to ensure the adjustment and minimal data discrepancies between the photogrammetric flight and the LiDAR point cloud of this 2010 reference flight.The discrepancies between them are originated as a consequence of the independent orientation systems of these sensors-despite the simultaneous recording of data-and thus its different positioning, which may prompt a lack of geometrical coincidence in the data compiled.This raises the need for an optimal adjustment between the two types of data [64].
In accordance with the above, the reference data are compiled and refined by three steps: 1. Photogrammetric orientation using in-flight GNSS and inertial data to ensure the geometrical quality of the whole photogrammetric block.This procedure was performed by means of a photogrammetric workstation and Socet Set 5.2 software (BAE Systems Plc., London, UK) [65] by measuring about 200 tie points and readjusting the orientation parameters.The RMS errors of the resulting model, based on the residual of 20 check points surveyed with GNSS (GNSS-based check points), are 0.258 m in XY and 0.100 m in Z (Table 2).

2.
LiDAR data orientation, by readjustment of different strips, which is employed to obtain a single homogeneous point cloud, even in overlapping zones.This task was carried out with the TerraMatch module of TerraSolid software (Terrasolid Ltd., Helsinki, Finland) [66].

3.
Analysis of height coincidence of both datasets.For this analysis an internal control process of the LiDAR point cloud was made using check points (model-based check points) taken from the refined photogrammetric model of the step 1.These points were selected taking into account the following: • Points on stable and flat surfaces, avoiding high slopes or roughness zones, were selected.

•
The points had to be well distributed throughout the study area and the number high enough to get an appropriate redundancy and significant results.Taking into account the good adjustment of the image orientation process, with mean values of about 0.10 m in GNSS-based check points, what really was tested in Step 3 were the errors of the LiDAR data.The vertical accuracy of aerial LiDAR data is well documented in many works [47][48][49]51,52] although the values present a great dispersion-from 0.16 to 3.26 m [47]-depending on different factors (flight altitude, field of view, beam angle respect to the vertical, coverage reflectance, vegetation, etc.).Furthermore, according to [67], the vertical error in LIDAR increases by 0.1 m per 1000 m of flight altitude, so it reaches values of about 0.10-0.20 m in conventional flights.In this sense, some works, in which the point density is usually higher than 1 point/m 2 , establish a vertical accuracy of 0.10-0.20 m [48, 49,52,68]; but in other works, with point densities lower than 1 point/m 2 , the accuracy is about 0.50 m [47,50].Thus, the higher the flight is, the lower the point accuracy and density, and therefore the vertical errors and uncertainties usually increase with the flight altitude.Meanwhile, in [47] the vertical discrepancies-expressed as RMS-between two LIDAR surveys in a close time interval were of 0.28 m, measured in stable zones (out of a mask drawn with the landslide inventory); however, the calculated uncertainties from comparison of LiDAR data with GNSS-measured points ranges from 0.50 m in leaf-off conditions and 0.75 m in leaf-on conditions.Bearing in mind these previous data, a reference value of accuracy of 0.20 m has been established for the present work, according the flight altitude and LiDAR configuration.
Nevertheless, a set of 59 model-based check points was extracted in a wide area surrounding the studied landslide to test the adjustment of LIDAR data and the photogrammetric model.The distribution of these points, all of them located outside of the landslide limits, is shown in Figure 5.The discrepancies or errors between the datasets range from 0 to near 1 m, with a mean positive value of 0.529 m, as can be observed in Table 3 and Figure 5a,c.Considering the good adjustment of the photogrammetric orientation process, this analysis allows us to infer a systematic error upwards in LiDAR data, corresponding to a translation of about 0.5 m, due to weak absolute georeferencing of the whole LiDAR block.This discrepancy could be eliminated by a simple translation downwards in height of the LiDAR data equal to the mean error.Meanwhile, the standard deviation (0.260 m) was of the same order of expected accuracy of LiDAR [48, 49,52,67,68] and remained constant after correction.The discrepancies in the points after correction are also shown in Table 5 and Figure 5b,c.Thus, in the new situation, a refined adjustment has been reached with a zero mean error and a standard deviation around the considered accuracy of LiDAR data.Thus, in the new situation, a refined adjustment has been reached with a zero mean error and a standard deviation around the considered accuracy of LiDAR data.

Orientation of Historical Photogrammetric Flights
Once the 2010 reference flight was re-oriented and the agreement between the photogrammetric and LiDAR data was ensured, the historical flights were oriented with respect to the 2010 reference system.The processing of several historical flights corresponding to a broad time span is rather difficult as the field ground control points (GCPs) must be stable, accessible, and unequivocally recognizable, not only in images from the time period considered but also in the field, where the points have to be measured with differential GNSS observations.Thus, this paper proposes the use of GCPs transferred from the reference photogrammetric flight (second-order GCPs) recognizable in flights corresponding to each survey.As the GCPs are directly located on the photogrammetric flights (the one to be oriented and the reference flight), it will be possible to certify in real time whether the point is observable and unequivocally recognizable on both flights.In addition, with a zenithal perspective the observer would have access to the entire study area, without limitation in accessibility

Orientation of Historical Photogrammetric Flights
Once the 2010 reference flight was re-oriented and the agreement between the photogrammetric and LiDAR data was ensured, the historical flights were oriented with respect to the 2010 reference system.The processing of several historical flights corresponding to a broad time span is rather difficult as the field ground control points (GCPs) must be stable, accessible, and unequivocally recognizable, not only in images from the time period considered but also in the field, where the points have to be measured with differential GNSS observations.Thus, this paper proposes the use of GCPs transferred from the reference photogrammetric flight (second-order GCPs) recognizable in flights corresponding to each survey.As the GCPs are directly located on the photogrammetric flights (the one to be oriented and the reference flight), it will be possible to certify in real time whether the point is observable and unequivocally recognizable on both flights.In addition, with a zenithal perspective the observer would have access to the entire study area, without limitation in accessibility or observation.Furthermore, with this procedure the number of required GCPs measured in the field is optimized, as well as the field work and processing time.
The main constraint of this methodology is a theoretical reduction in the accuracy of the coordinates in the second-order GCPs with regard to the accuracy of the first-order GCPs due to error propagation.Nevertheless, as the 2010 photogrammetric flight with higher resolution and a RMS lower than the pixel size was always taken as the reference, the accuracy in the GCPs transferred in this way is proved to be better than the resolution of the historical flights.
The process of transferring points from the reference flight (more recent) to the previous flights to be oriented is made iteratively following the next steps (Figure 3): Approximate pre-orientation of the flight to be oriented in the reference system.2.
Once the two photogrammetric flights are pre-oriented in the same reference system, a set of several common points in both flights are located and uploaded simultaneously in the digital photogrammetric workstation.

3.
All these points located and measured in the reference flight are considered second-order GCPs for the flight to be oriented (Figure 6).This process includes full GCPs (known X,Y,Z), horizontal GCPs (known X,Y), and height GCPs (known Z).Table 2 shows the number and type of the second-order GCPs used in the orientation of the historical flights.

4.
The flight to be oriented is adjusted and the process is repeated from step 2 until a final proper orientation is attained, that is when the difference between coordinates of GCPs in both datasets was lower than the image resolution of the flight to be oriented.
Geosciences 2017, 7, 32 9 of 23 or observation.Furthermore, with this procedure the number of required GCPs measured in the field is optimized, as well as the field work and processing time.
The main constraint of this methodology is a theoretical reduction in the accuracy of the coordinates in the second-order GCPs with regard to the accuracy of the first-order GCPs due to error propagation.Nevertheless, as the 2010 photogrammetric flight with higher resolution and a RMS lower than the pixel size was always taken as the reference, the accuracy in the GCPs transferred in this way is proved to be better than the resolution of the historical flights.
The process of transferring points from the reference flight (more recent) to the previous flights to be oriented is made iteratively following the next steps (Figure 3): 1. Approximate pre-orientation of the flight to be oriented in the reference system.2. Once the two photogrammetric flights are pre-oriented in the same reference system, a set of several common points in both flights are located and uploaded simultaneously in the digital photogrammetric workstation.3.All these points located and measured in the reference flight are considered second-order GCPs for the flight to be oriented (Figure 6).This process includes full GCPs (known X,Y,Z), horizontal GCPs (known X,Y), and height GCPs (known Z).Table 2 shows the number and type of the second-order GCPs used in the orientation of the historical flights.4. The flight to be oriented is adjusted and the process is repeated from step 2 until a final proper orientation is attained, that is when the difference between coordinates of GCPs in both datasets was lower than the image resolution of the flight to be oriented.After applying this procedure, the spatial correspondence between the two compared flights is checked out under the premise of finding a total correspondence in stable zones where no modifications were detected between the two surveys.The comparison was made by overlaying the reference 2010 LiDAR data and the stereomodels of the historical flights with the stereoscopic viewing tools of the Socet Set and the photogrammetric workstation.If the orientation was correct, then DTM contours and points derived from LiDAR data had to be adjusted to the ground in the stable areas of the study zone.If significant discrepancies were found, then an orientation re-adjustment would be necessary.In the case studies, the stereoviewing display comparison showed a good coincidence between elevation and stereoscopic models in all the flights, and therefore no readjustments of the flight orientation were necessary.
For this comparison, different types of color compositions were applied in recent flights: true color (Blue-Green-Red) and false color (Green-Red-NIR), which was possible in recent flights (2008 and 2010), since they incorporated the near infrared (NIR) band.With these compositions the different ground features were more visible and identifiable in addition to the element over them, such as vegetation, which is very sensible to detection with near infrared.

Building and Editing DTMs
A Digital Surface Model (DSM) was obtained from the 2010 LiDAR data.Then a classification was made between ground-points and non-ground-points (vegetation, buildings, artifacts, etc.).This classification process was performed with the software TerraSolid (module TerraScan).Once this classification was finished, the ground surface corresponding to the reference survey was defined by those filtered ground-points, thus obtaining the Digital Terrain Model (DTM).
From the initial 4,240,499 points in the 2010 LiDAR point cloud, 1,943,545 were classified as ground-points while 2,296,954 were non-ground points.Then, the results of this classification were displayed in the stereo photogrammetric workstation together with the reference photogrammetric flight, and edited using the stereoscopic tools supplied by Socet Set.This editing process consisted basically on correcting wrong classified points and stereoplotting breaklines and some typical ground features, such as scarps, water streams, slope changes, etc., which allowed for a more realistic and precise DTM generation.
After preparing a reference DTM from LiDAR data classification, the models corresponding to the historical flights are obtained by editing and outdating the 2010 reference model.Other authors [7,20,21,23,24,26] have generated DEMs corresponding to each survey by a process of image correlation independent for each flight.However, in this work another approach has been adopted for the following two reasons.First, the accuracy of models is consistent with the resolution and the radiometric quality of the aerial images, which in the older flights has led to a lower accuracy in the models.A second reason is the need to build the model again by matching, which in addition requires the editing of all the surveys in the study zone, no matter if there have been any changes in the terrain or not, thereby giving rise to inevitable false ground modifications in the comparison of models corresponding to different surveys and noise in the detection of changes.
For these reasons, an outdating of the reference model-of a higher quality and geometric resolution than those models that would be obtained by a matching technique in each case-is proposed.This outdating process consists of uploading the reference model overlaid to the historical flights in the photogrammetric workstation.Then, by means of the stereoscopic viewing tools is possible to analyze the lack of coincidence between the two datasets.When that lack of coincidence occurs, the 2010 DTM is edited in that zone to adapt the terrain to the ground surface in the stereomodel of the historical flight.In the case of changes affecting large surfaces, it is necessary to launch a matching process exclusively for that area.In ground areas where no change is detected, the model is not modified, thus reducing the editing process time as well as avoiding noise in the analysis of changes and maintaining a better overall quality of the model provided by the higher positional accuracy of the 2010 data.
After finishing the outdating of the reference model to the immediately previous historical survey under analysis, this resulting model was used as a reference model in the preceding survey, repeating the process in all cases until the 1956 DTM was obtained.In Figure 7, the DTMs generated in the unstable area are shown, allowing a visual observation of significant changes between these six dates considered.repeating the process in all cases until the 1956 DTM was obtained.In Figure 7, the DTMs generated in the unstable area are shown, allowing a visual observation of significant changes between these six dates considered.

Comparison between Models and Calculations
The comparison between models was carried out from different perspectives: 1. Vertical displacements by subtracting DTMs (differential DTMs) (Figure 7).They may be negative or positive, depending on whether the model compared with the reference lies below or above, allowing the identification of areas of mass depletion or accumulation, respectively.2. Volumetric calculations between models.As in the previous case, the volumes can be negative (depletion areas) or positive (accumulation areas).3. Analysis of DTMs profiles through the landslides permitting a better observation of terrain displacements (Figure 8).With enough profiles, it is possible to attain an overview of the landslide kinematics, for instance if it is planar or there are rotations around vertical or horizontal axis indicating rotational sliding.

Comparison between Models and Calculations
The comparison between models was carried out from different perspectives: 1. Vertical displacements by subtracting DTMs (differential DTMs) (Figure 7).They may be negative or positive, depending on whether the model compared with the reference lies below or above, allowing the identification of areas of mass depletion or accumulation, respectively.2.
Volumetric calculations between models.As in the previous case, the volumes can be negative (depletion areas) or positive (accumulation areas).

3.
Analysis of DTMs profiles through the landslides permitting a better observation of terrain displacements (Figure 8).With enough profiles, it is possible to attain an overview of the landslide kinematics, for instance if it is planar or there are rotations around vertical or horizontal axis indicating rotational sliding.These calculations and analyses have been made by means of the Maptek I-Site Studio 6.0 software (Maptek Pty. Ltd., Adelaide, Australia) [69] and ArcGIS 10.2 (Esri, Redlands CA, USA) [70].The errors and uncertainties-usually expressed as root mean square (RMS) and standard deviation (SD), respectively-of these procedures have to be estimated in order to validate the results of these calculations.In the case of the 2010 flight, the value of standard deviation (0.26 m) is assumed to be vertical uncertainty or the minimum level of detection [52,71,72].Meanwhile, from previous works [7,23,32] 4).These calculations and analyses have been made by means of the Maptek I-Site Studio 6.0 software (Maptek Pty. Ltd., Adelaide, Australia) [69] and ArcGIS 10.2 (Esri, Redlands CA, USA) [70].The errors and uncertainties-usually expressed as root mean square (RMS) and standard deviation (SD), respectively-of these procedures have to be estimated in order to validate the results of these calculations.In the case of the 2010 flight, the value of standard deviation (0.26 m) is assumed to be vertical uncertainty or the minimum level of detection [52,71,72].Meanwhile, from previous works [7,23,32] 4).Uncert.YEAR1-YEAR2 = (Uncert.YEAR12 + Uncert.YEAR2 2 ) 0.5 (2) In this case, the uncertainties are lower than 0.40 m for the 2008-2010 period, 0.65-0.70m for the 2001-2008, 1992-2001, and 1984-1992 periods, and close to 1 m for the 1956-1984 period (Table 4).
In addition, the uncertainties of the volumetric calculations are estimated by multiplying the vertical uncertainties by the areas of depletion and accumulation [52].The uncertainty of the waste material (difference between depletion and accumulation volumes) is calculated as a propagated error (square root of the sum of squares of depletion and accumulation volumetric uncertainties).These last ones range from 21 × 10 3 m 3 for the 2008-2010 period, 31-33 × 10 3 m 3 for the 2001-2008, 1992-2001, and 1984-1992 periods, to about 40 × 10 3 m 3 for the 1956-1984 period (Table 4).

Results
The analysis of orthophotographs and models since 1956 of the area located in the Guadalfeo River, near Almegíjar village, has allowed for the identification of a landslide as well as its characteristics and evolution.The landslide does not appear in the aerial photographs of the 1956 and 1984 flights, so it was triggered after 1984, although located in a broadly unstable zone, with evidence of older activity upslope.Thus in the stereoscopic models of later dates (1992, 2001, 2008, and 2010), the landslide area presents some clear features that include a clean scarp, steps and terraces, minor scarps, and an accumulation zone invading the river channel.The involved mass reaches about 300 × 10 3 m 3 according to the volumetric calculations, so the landslide can be defined as moderately large.Other features identified through aerial photographs, but also field work, are a large extent of debris flows and the development of pervasive erosion in rills and gullies on the weathered rock-mass surface, even before 1984 and at the start of the main landslide.
The differential models between the 1956 and 2010 surveys (the whole period considered) show two clearly differentiated areas, the depletion area at the upper part of the mass and the accumulation area in the lower part as shown in Figures 7 and 8.As the differential models are calculated by subtracting the older models to the recent ones, negative differences in elevation and volume can be observed in the depleted area, and positive differences are recorded in the accumulation zone.In the Table 5, the proportion of areas of depletion, accumulation, and below the uncertainty is shown.If uncertainty is not taken into account, the areas of depletion and accumulation are practically compensated for (near 50%), although in some periods (1992-2001 and 2008-2010) depletion areas predominate slightly.When the uncertainty areas are introduced, some periods show a high proportion of areas below the uncertainty, such as 1956-1992 and 2001-2008 (less than 10%) or even 2008-2010 (16%).In these cases, the depleted and accumulated areas are compensated for except in the 2008-2010 period where the depleted area predominates.Meanwhile, other periods such as 1984-1992, 1992-2001, and the global periods show a higher proportion of area depleted and accumulated (over 70%), usually compensated for, except in the 1992-2001 period, when depleted area also predominates.Significant vertical displacements are recorded in the landslide area, as shown in Table 6, with mean values higher than an order of magnitude with respect to the estimated uncertainties (Table 4).Thus, in the depleted area, there is an average descent of 8.5 m of the ground surface, with maximum values around 28 m.In addition, in the accumulation area, vertical displacements related to the elevation of the ground surface show an average value of 3.9 m, with maximum peaks of almost 28 m.Considering that the landslide was triggered after 1984, the annual displacement rates reach at least 0.32 m/year and 0.15 m/year in the depleted and accumulation areas, respectively.Consequently, the landslide velocity can be estimated on the order of some few cm/year so the slide can be catalogued as "very slow" if it were considered continuous and uniform in time [73,74].
The analysis of the results of all five time periods considered (Table 6), between the six available flights, reveals that the vertical mean displacements are practically null in the first period of 1956-1984 (0.09 m and 0.08 of terrain elevation in the depletion and accumulation areas), higher in the second period  The volume assessment in Table 7 shows a depletion of more than 290 × 10 3 m 3 of mass in the upper part of the slide, but the accumulated mass in the lower part near the river channel is about 120 × 10 3 m 3 .This implies mass losses of nearly 175 × 10 3 m 3 , probably evacuated by the river flow.These values are slightly modified to 280 × 10 3 m 3 , over 110 × 10 3 m 3 and 167 × 10 3 m 3 if a starting date of 1984 is considered.By periods, the values are very low in the first period of 1956-1984 (24 × 10 3 m 3 of depleted material, 17 × 10 3 m 3 of accumulated material and 8 × 10 3 m 3 of waste material), clearly higher in the second period of 1984-1992 (169 × 10 3 , 53 × 10 3 and 116 × 10 3 m 3 , respectively), lower in the third period of 1992-2001 (86 × 10 3 , 46 × 10 3 and 40 × 10 3 m 3 ), and far lower in the last two periods of 2001-2008 (9,2 × 10 3 , 8,7 × 10 3 and 0,5 × 10 3 m 3 ) and 2008-2010 (18 × 10 3 , 7 × 10 3 and 11 × 10 3 m 3 ).The values corresponding to the periods 1984-1992 and 1992-2001 are clearly higher than the uncertainties in Table 4, but the values of the periods 1956-1984 and 2001-2008 are lower than the uncertainties.The values of the period 2008-2010 are on the same order as the uncertainties, higher in depleted material and lower in accumulated and waste materials.

Discussion
As the landslide does not appear in aerial photographs from the 1956 and 1984 flights, the origin of the landslide had to be after 1984, probably in 1989, the rainiest year in this period [75].However, due to the uncertainty about this date, the calculations presented in Tables 6 and 7 are made assuming the landslide origin in 1984, so the rate values should be considered as minimum values for the period 1984-1992.
The development of recent features, such as a clean scarp, steps and terraces, minor scarps, and an accumulation zone invading the river channel, can be clearly observed in the orthophotographs and models from 1992 to 2010, which evidences significant activity in the last 30 years.As mentioned before, other features identified in the orthoimages of 1956 and 1984 such as debris flows and the development of rills and gullies allow the zone to be described as unstable, with some evidence of older activity upslope.
Furthermore, landslides and these other features evidence a relevant geomorphic activity in the study area, confirmed by large mass losses associated with seasonal "rambla" flash-flooding events.Thus, the recorded negative displacements or volume changes (ground decrease or mass erosion) are larger than the positive changes (ground elevation or mass accumulation), so it can be deduced that the resulting mass losses are due to their evacuation by river flows.The great volume of mass waste is consistent with the significant geomorphic activity and the high indexes of tectonic activity along the southern flank of the Sierra Nevada, where the Almegíjar landslide is set [58].This currently visible activity is related to the overall tectonic evolution of the region since Late Neogene and through the Holocene [76].
The parameters related to the landslide and determined from the calculation and comparison of models, longitudinal sections, and other evidence (Tables 5-7 and Figures 7 and 8) allow for the assessment of its dimensions and morphology, as well as the identification of trends concerning age, velocity, activity, diachroneity, and intensity following published classifications [63,73,74,[77][78][79][80][81].
Thus the affected area is about 75 × 10 3 m2 , with a balanced distribution of depletion and accumulation areas (near 50%, although with a certain predominance of the depletion areas) that leads to a movement of rockslide type without development of flow mechanism, in which the material can expand downslope.In this sense, according to [60], the landslide, reactivated in the period 2009-2010, shows a strain shortening along the slope with an extension perpendicular to this axis, resulting in a bulging in the lower sector of the displaced mass, following a mechanism described by [61].The morphology of the landslide scarp and the profiles suggest a planar sliding surface and a mass thickness ranging from 40 to 70 m, partially eroded and affected by debris flow of the weathered rock-mass surface.According to it, the movement can be classified as a rock translational or planar slide [74,77] without evidence of a rotational sliding component in the models and sections.The mobilized (depleted) mass has a volume of about 300 × 10 3 m 3 , which leads us to categorize the landslide as "moderately large" [63,78,80].Part of the mobilized mass is now displaced downslope, while another part invaded the river channel and was evacuated by the water flow, which has also limited the development of the accumulation areas.
The age is catalogued as "contemporary" [63] and its velocity-in average terms-can be considered as "very slow" [73,74], with vertical displacements of 0.32 and 0.15 m/year in the depleted and in the accumulated area, respectively, over the whole period of 26 years (1984-2010).However, as discussed below, this velocity is not uniform along this period and the landslide has shown periods of different rates of movement, with inactive and active phases, where the velocity ranges from "extremely slow" to "very slow."Thus, more difficult is the assessment of the landslide activity and diachroneity, because the results do not provide a detailed pattern of active periods, which is particularly difficult in these regions where the triggering factors, mainly the rainfall, do not follow a regular or seasonal pattern.However, the analysis of the differential changes between DTMs for the five time periods considered, and some additional data from terrestrial laser scanner (TLS) techniques [60], make it possible to assess certain activity features of interest.
The style of activity [81] is "single," corresponding to the displacement of a well-delimited mass in this single planar slide.On the other hand, the distribution of the activity may be considered "in advance" and also "confined," as the failure plane-shown only at the slide crown of the main scarp-is below the mass at slide toe.Regarding the state of activity [81], the movement can be considered active, although alternating periods of low activity with more active periods during ephemeral floods of the "rambla" [53,58,59].
Furthermore, the results described in the previous section show some periods of activity with displacements rates-measured in the vertical direction-on the order of decimeters by year.Regarding them, the accuracy of the technique has to be taken into account.The accuracy has been established in this work from errors of the flight orientation at check points that allow the estimation of the uncertainty of DTMs and differential DTMs.This uncertainty ranges from 0.36 m in the most favorable case (2008-2010) to 0.95 m in the least favorable , with the other cases between 0.60-0.70m.Thus the vertical displacements measured in the periods 1984-1992 and 1992-2001, clearly higher than uncertainty, can be considered significant, but those measured in the periods 1956-1984 and 2001-2008 (next to zero values) are considered non-significant.The displacements calculated in the short period of 2008-2010-slightly higher than uncertainty-are at the limit of accuracy of the analysis, although additional data confirm the activity in this period, as will be discussed below.The volumetric accuracy and calculations agree with this, and therefore the depleted, accumulated, and waste volumes in the periods 1984-1992 and 1992-2001 are significant (much larger than the uncertainty), while the other periods (1956-1984, 2001-2008, and 2008-2010) are not significant (lower or similar to the uncertainty).
With this premise, the movement started after 1984 with the maximum rates observed (0.46 m/year of descent in the head or depleted area and 0.39 m/year in the foot or accumulation area), and slowed in the next period of 1992-2001 when lower displacements are observed (0.32 and 0.20 m/year, respectively), until the landslide seemed to stop in the period of 2001-2008, when the measured displacements are not significant.After that, in the last period of 2008-2010 a reactivation of the landslide occurred with displacement rates of about 0.21 and 0.12 m/year.The volumetric calculations confirm this evolution, with volume rates ranging between 10 to 20 × 10 3 m 3 /year in the active periods and values near to 0 in the inactive periods.
At this point, we must mention the complementary results found in this landslide by the combined use of Terrestrial Laser Scanner (TLS) and Global Navigation Satellite Systems (GNSS) between 2008 and 2010 [60], which showed a reactivation in the period March 2009 to June 2010-coinciding with a year of intense rainfall-after an inactive period between September 2008 and March 2009.In the reactivation period a downward movement of 1.20 m at the top was measured, which corresponds to a descent rate of 0.96 m/year; meanwhile, below the middle part of the mass, 1.30 m of advance was established with a displacement rate of 1.04 m/year.These values of the annual rates about 1 m/year express periods of greater activity with a velocity near to that catalogued as "slow" (higher than 1 m/year) between periods of much lower activity in which the velocity is "very slow" to "extremely slow".These data suggest that other vertical displacements and volumes measured between the flights used in this study could have taken place in shorter intervals of time than the periods considered in this work.
Thus, in a more speculative way, if the date of 1989 were considered as the starting point of the landslide-triggered by a period of heavy rainfall [75]-the estimated rates of displacements calculated for the period of 1984-1992 (decimeters by year) would increase to 1.7 m/year in the depleted area and almost 1 m/year in the accumulation area (velocity catalogued as "slow"), on the same order of magnitude as those calculated with TLS measurements [60].If the period in which the displacement occurred were even shorter, the displacement rates could be even higher, which expresses a large initial impulse for the landslide.The same situation could have occurred in other active periods such as 1992-2001, in relation to the heavy rainfall in the years 1995 to 1998, particularly in the winter of 1996-1997.Again, the displacement rates could increase to values near 1 m/year or higher.
For these landslides of discontinuous activity, an interesting concept is the activity duration or diachroneity [63].According to the scale established in these works, the degree of diachroneity assigned to the studied landslide is VIII (10 to 100 years).Nevertheless, if only the periods of activity are considered, the degree of diachroneity would be VI (1 month to 1 year).
The triggering mechanism of the landslide was the rainfall, as in other regions near to it [55] and in general in Mediterranean countries [82].In these areas the movements are triggered occasionally by "rambla" flash floods cutting down through the confined planar slide mass and dragging out the mass accumulation [75,83].These flood events occur at the beginning of autumn with heavy rain (cold drop phenomenon) or at the beginning of spring after a rainy winter, in which heavy rains coincide with snow melting from the Sierra Nevada range.
Thus, in an approximate way, a relationship between the landslide activity and the rainfalls can be established.The landslide probably was triggered in relation to a rainy event of about 800 mm between November 1989 and January 1990 (more than 400% of the average rainfall for this period).Other equivalent periods occurred between November 1996 and January 1997 with more than 600 mm (more than 300% of the average) [53,55,59,75,83], and in the early months of 2010 with 500 mm between December 2009 and March 2010 (more than 200% of the average rainfall).These events and the general rainfalls of the corresponding years (1989-1990, 1995-1998, and 2009-2011) are related to the higher activity of the periods 1984-1992, 1992-2001, and 2008-2010, respectively.Meanwhile, the lesser activity in the period 2001-2008, when minimum displacements and volumetric changes were recorded, is associated with the absence of relevant precipitation.
This irregular pattern of landslides and rainfall coincides with other evidence such as field observations and landslide inventories along the Betic Cordillera since 1979 [54-57], particularly in the vicinity of this landslide [53,58,59,75,83].The analysis of the available regional climatic data, reliable in the region since 1940, shows some intense rainy periods similar to those recorded in 1989, 1997, and 2010 [54], triggering slope failures in the study zone.These rainfall patterns, spanning periods of two or three years, register intense peaks in some months.First, we can see five periods exceeding 400 mm in three months, when the mean annual rainfall is 620 mm/y, and also another four periods of accumulated rainfall of more than 300 mm.Thus, the return period for rainfall triggering landslides in the region is around 15 years for the more intense events and eight years if less intense events are taken into account.Other analyses, based on average annual rainfall, show return periods of 18 years for annual precipitation higher than 950 mm, and five years for precipitation higher than 750 mm, these reaching the threshold of new landslides in the region [84].Finally, recent studies based on the analysis of the anomalies in antecedent cumulative rainfall establish return periods of 12.4 years for these anomalies [85], a value similar to the previous studies.
Comparing these data with climatic studies on the origin of the precipitation in southern Spain, the intense rainfall in the region appears to be related to a negative value of the North Atlantic Oscillation index (NAO index, [86]).Positive values, which are more frequent, are related to more intense precipitation in Northern Europe and drier periods with rainfall in the Mediterranean area.As the periodicity of winter values of this NAO index is eight years [87], this can be seen as coincident with the return periods identified.

Conclusions
In conclusion, the real usefulness of LiDAR and digital photogrammetry techniques, as well as the use of differential DEMs, is demonstrated in landslide studies, particularly concerning the determination of the temporal evolution of landslides.The accuracy and consistency of these techniques enables the detection of very detailed slope features and their changes.In this work a new methodology has been developed, based on the use of a stereo photogrammetric workstation for outdating DTMs of more recent flights by editing them only in unstable zones (with clear changes between DTMs), over the stereoscopic models of previous flights.This permits us to reduce the time spent editing DTMs and especially to avoid noise in the differential model calculations.
Concerning the results, the methodology has enabled the detection and quantification of the extent of depletion and accumulation areas of the landslide, with average vertical displacement of about 8.5 m in the depletion area and 4 m in the accumulation area, attaining maximum displacements of almost 30 m in both areas.These displacements can be considered significant because they are an order of magnitude higher than the estimated uncertainty.Also, significant volumetric changes of about 300 × 10 3 m 3 in depletion area and 120 × 10 3 m 3 in accumulation areas were quantified, with some 175 × 10 3 m 3 of materials evacuated afterwards by the river flow.
The displacement and volumetric change rates are not constant over time, attaining maximum values in the initial period 1984-1992 (displacements of about 40 cm/year) and more moderate in the periods 1992-2001 and 2008-2010 (displacements of 15-30 cm/year), during seasons of heavy rainfalls, particularly in 1997 and 2010.The lowest rates were recorded in the period 2001-2008, when the displacement hardly reached 4 cm/year in a long dry period and they can be considered non-significant regarding the accuracy of the technique used.Other observations and data from TLS lead to estimates of displacement rates of higher than 1 m/year in some recent reactivation (2009-2010), similar to that found during the landslide-triggering phase.This difference in slide activity features diachronic slope movements alternating periods of variable activity depending on the influence of triggering factors.From these data, an "extremely to very slow and slow" slide velocity is calculated.
Nevertheless, current research could provide more detailed results on the slide displacements and evolution by using more flights, mainly available in the last decade, and incorporating TLS and UAV photogrammetric data, as well as displacement wireless sensor networks (WSN) and local rainfall data, for easier forecasts of future slide displacements.On the other hand, other techniques such as the determination of displacement vectors in significant slope points or the calculation of absolute distances provide more precise information on the slide kinematics.Finally, regional quantitative studies of other landslides are also needed to draw better quality landslide-hazard maps.

Figure 3 .
Figure 3. Flow chart of the methodology.

Figure 3 .
Figure 3. Flow chart of the methodology.

Figure 4 .
Figure 4. Distribution of available flights.Coordinates are in ETRS89/UTM30.These historical flights were digitized with a precision photogrammetric scanner Vexcel Ultrascan 5000 at a resolution of 15 microns, which, in turn, implies a ground sample distance (GSD) of 0.60 m in the 1956 flight and 0.30 m in the 1992 and 2001 flights.The 1984 flight was already supplied in digital format after a scanning at 20 microns that provided a GSD of 0.80 m.The period considered (c.50 years) is long enough to study short-term changes in the relief and, more precisely, to analyze the landslide evolution in the study zone.

Figure 5 .
Figure 5. Distribution of errors or discrepancies in height between 2010 LiDAR and photogrammetric data: (a) before the adjustment; (b) after the adjustment; (c) histograms.

Figure 5 .
Figure 5. Distribution of errors or discrepancies in height between 2010 LiDAR and photogrammetric data: (a) before the adjustment; (b) after the adjustment; (c) histograms.
it has been observed that the vertical errors and uncertainties in DEMs obtained by means of matching from photogrammetric oriented flights are about two or three times higher than the orientation residual errors in Z. Therefore the vertical uncertainties of DTMs obtained in this work range from 0.27 to 0.30 m for 2008 and 1992 flights, about 0.60 m for 2001 and 1984 flights, to over 0.70 m for the 1956 flight (Table
it has been observed that the vertical errors and uncertainties in DEMs obtained by means of matching from photogrammetric oriented flights are about two or three times higher than the orientation residual errors in Z. Therefore the vertical uncertainties of DTMs obtained in this work range from 0.27 to 0.30 m for 2008 and 1992 flights, about 0.60 m for 2001 and 1984 flights, to over 0.70 m for the 1956 flight (Table of 1984-1992 (3.70 m of decrease in the depletion area and 3.11 m of elevation in the accumulation area), lower in the third period of 1992-2001 (2.87 m and 1.78 m, respectively), and far lower in the last two periods of 2001-2008 (0.30 and 0.27 m) and 2008-2010 (0.43 and 0.23 m).These partial values are significantly higher than the uncertainties of Table 4, except for the periods of 1956-1984 and 2001-2008 (clearly lower than uncertainties) and the period 2008-2010 (similar to uncertainties).More informative is the comparison of annual rates, with values clearly above average in the period 1984-1992 (0.46 m/year of depletion and 0.39 m/year of accumulation); values near average in the period 1992-2001 (0.32 and 0.20 m/year respectively), and slightly lower in 2008-2010 (0.22 and 0.12 m/year); finally, rate values are clearly lower in the period 2001-2008 (0.04 and 0.04 m/year) and practically null in the period 1956-1984.

Table 1 .
Properties of the image datasets.

Table 1 .
Properties of the image datasets.

Table 3 .
Statistics of errors between LiDAR and photogrammetric data of the 2010 flight.
1Errors are in meters.

Table 3 .
Statistics of errors between LiDAR and photogrammetric data of the 2010 flight.
1Errors are in meters.

Table 4 .
Estimated uncertainties of DTMs, differential models, and volumetric calculations.

Table 4 .
Estimated uncertainties of DTMs, differential models, and volumetric calculations.

Table 5 .
Areas resulting in from the comparisons between DTMs.