Multitemporal Landslide Inventory and Activity Analysis by Means of Aerial Photogrammetry and LiDAR Techniques in an Area of Southern Spain

: This paper deals with the use of aerial photogrammetry and LiDAR techniques to analyze landslide activity over a long time span—just over 32 years. The data correspond to several aerial surveys (1984, 1996, 2001, 2005, 2009, 2010, 2011, 2013 and 2016) covering an area of about 50 km 2 along highway A-44, near Ja é n (Southern Spain). An ad hoc combined photogrammetric and LiDAR aerial survey of 2010 was established as the reference ﬂight. This ﬂight was processed by means of direct orientation methods and iterative adjustments between both data sets. Meanwhile, historical ﬂights available in public geographical data servers were oriented by transferring ground control points from the reference ﬂight. Then, digital surface models (DSMs) and orthophotographs were generated, as well as the corresponding differential models (DoDs), which, after the application of ﬁlters and taking into account the estimated uncertainty of ± 1 m, allowed us to identify true changes on the ground surface. This analysis, complemented by photointerpretation, led us to obtain a landslide multitemporal inventory in the study area that was analyzed in order to characterize the landslide type, morphology and activity. Three basic typologies were identiﬁed: rock falls–collapses, slides and ﬂows. These types present different morphometric properties (area, perimeter and height interval) and are associated with different conditions (height, slope, orientation and lithology). Moreover, a set of monitoring areas, common for the different ﬂights, was also used to analyze the activity throughout the study period. Thus, some more active periods were identiﬁed (2009–2010, 2010–2011, 2011–2013 and 1996–2001) among other less active ones (1984–1996, 2001–2005, 2005–2009 and 2013–2016), which are related to rainy events and dry years, respectively.


Introduction
The estimation of landslide activity and its dating is crucial in risk analyses, since it allows the determination of hazards in terms of spatiotemporal probability, one of the components of the risk equation [1]. Landslide dating is difficult in most cases, so it is usually reduced either to the direct dating of a reduced sample of individual landslides or to the dating of the triggering factors in order to determine the activity in an indirect manner. Direct dating is performed with complex and sometimes invasive techniques, which restricts its use generally to individual cases. For wider areas, geomorphological analysis is a primary tool for obtaining multitemporal inventories [2,3].
This geomorphological analysis is usually based on photointerpretation, which is a time-consuming technique and introduces a certain subjectivity into the analysis [3,4]. However, in recent years, geomatic techniques have opened up an interesting research line  In the area, diverse materials appear, belonging to different tectonic [40] units of a complex and widely discussed structure ( Figure 1):  Triassic materials, consisting of clays, shales, carbonates and evaporites;  Middle Subbetic, formed by thick carbonated series (Jurassic) overlayed by marls and marly limestones (Cretaceous);  Intermediate Units, where the previous series are repeated with variations;  Prebetic, in which marls, marly limestones, limestones and dolomites (Cretaceous) outcrop;  Guadalquivir Units, where different materials predominate: evaporite and shale sections (Triassic), loams and clays (Cretaceous-Paleogene) and marly clayey sediments (Lower Miocene) [41];  Pliocene conglomerates and sands;  Quaternary fillings (colluvial and alluvial materials, debris, terraces and soils). In the area, diverse materials appear, belonging to different tectonic [40] units of a complex and widely discussed structure ( Figure 1):

Materials
Our methodology focuses on the combined use of aerial photogrammetry and LiDAR techniques, with data captured from conventional aerial platforms with decametric resolution. First, an ad hoc combined flight, with a digital camera (Z/I DMC) and LiDAR sensor (Leica ALS50-II), was flown in 2010 with GNSS/inertial measurement unit (IMU) systems for direct georeferencing and flight guidance. The aerial images have four radiometric bands (RGB and near infrared (NIR)), and the ground sample distance (GSD) of the images is 0.20 m (Table 1). Meanwhile, the LiDAR resolution is about 1-1.5 points/m 2 . Initially, this flight was planned to study the instability processes along the A-44 highway cuts, so several strips were flown along the length of the highway (Figure 1).   (Table 1). All data are available in different image formats (TIFF/JPG/ECW) from public spatial data infrastructures (SDIs) or data servers, such as the photo library of the IGN [42]

Methodology
The methodology consists of several phases partially described in previous studies [10,16,21,24,44,45], which can be summarized in the following steps ( Figure 3):

1.
Orientation of the 2010 reference flight (photogrammetric and LiDAR).

2.
Orientation of the remaining flights in the same coordinate reference system as the reference flight. 3.
Generation of the DSMs and orthophotographs. 4.
Elaboration of the multitemporal landslide inventory.

Orientation of the Reference Flight
First, the reference 2010 flight was processed by means of direct orientation tec niques based on the flight parameters (GNSS and IMU data) in the digital photogramme ric workstation (DPW) and the software Socet Set 5.6 [46]. The reference coordinate syste was ETRS89 UTM 30N (EPSG: 25830).
After the first orientation step, some residual misalignment between the photogram metric and LiDAR data persisted, despite the simultaneity of the acquisition process both data sets. In this case, the photogrammetric flight was finely aligned, using a set 25 height ground control points (Z GCPs) extracted from the LiDAR point cloud followin the methodology developed in [44,45]. Table 2 shows the results of this fine adjustmen Thus, the errors of the photogrammetric orientation of the 2010 flight, expressed as t root mean square (RMS), were 0.030 and 0.093 m for XY and Z, respectively.

Orientation of the Reference Flight
First, the reference 2010 flight was processed by means of direct orientation techniques based on the flight parameters (GNSS and IMU data) in the digital photogrammetric workstation (DPW) and the software Socet Set 5.6 [46]. The reference coordinate system was ETRS89 UTM 30N (EPSG: 25830).
After the first orientation step, some residual misalignment between the photogrammetric and LiDAR data persisted, despite the simultaneity of the acquisition process of both data sets. In this case, the photogrammetric flight was finely aligned, using a set of 25 height ground control points (Z GCPs) extracted from the LiDAR point cloud following the methodology developed in [44,45]. Table 2 shows the results of this fine adjustment. Thus, the errors of the photogrammetric orientation of the 2010 flight, expressed as the root mean square (RMS), were 0.030 and 0.093 m for XY and Z, respectively.

Orientation of the Remaining Flights
Next, the previous historical flights (1984, 1996, 2001, 2005 and 2009) and subsequent flights (2011, 2013 and 2016) were oriented according to a similar methodology used by authors in previous studies [10,16,44,45], where second-order GCPs located in stable areas were transferred from the 2010 reference flight. This procedure ensures that all flights are oriented in the same coordinate reference system as the 2010 reference flight. Furthermore, this approach implies that there is no requirement of field-surveyed GCP. Table 2 shows the number and properties of the GCPs for each flight.
As the second-order GCPs are directly located on the photogrammetric flights on the DPW, it is possible to ensure in real time whether the point is observable and unequivocally recognizable on both flights. Moreover, transferring GCPs to historical flights facilitates the location of these GCPs when the landscape has dramatically changed, and it is sometimes almost impossible to identify suitable present points in the field from old aerial images. The theoretical reduction in accuracy with the use of the second-order GCP is considered negligible, as a photogrammetric flight with higher resolution and a low RMS error (lower than the GSD of all flights) was used as a reference [16]. Thus, the accuracy estimated in the transferred GCPs was proven to be better than the resolution of the historical flights. Table 2 also shows the results of this process, expressing the RMS of the orientation residual errors. Thus, horizontal RMS XY errors ranged from values equal to or lower than 0. However, as the historical flights were oriented using transferred GCPs obtained from the reference 2010 flight, the precision of any measurement in those flights is influenced by the error propagation from the reference flight. Therefore, the propagated error (RMS Prop) of these flights can be calculated according to the following expression: Thus, RMS XY errors were very similar to the previous ones, since the error of the 2010 flight was very small, and thus they ranged from around 0.05 m to 0.27. Meanwhile, RMS Z errors ranged between 0.10 and 0.22 m.

Generation of the DSMs and Orthophotographs
After the orientation of all the flights, automatic image correlation (dense matching) and differential rectification techniques were used to generate the digital surface models (DSMs) and orthophotographs for every epoch, using the Next Generation Automated Terrain Extraction (NGATE) module of Socet Set 5.6 [46,47]. A DSM resolution of 5 times (2.5 m) the value of the original imagery GSD (0.5 m) was considered, following the conventional usage in photogrammetry of applying a multiple value of the GSD for DSM generation, which is gradually being reduced thanks to the improvement of algorithms for the automatic measurement of points. The GSD considered for the orthoimage generation is equal to the original imagery GSD (0.5 m). Finally, both the DSMs and orthophotographs were exported to a raster format (TIFF) in order to be used in the GIS software.
As shown in previous studies [11,16,45], the vertical uncertainties in DEMs obtained by means of image-matching in photogrammetric flights were about two to three times the RMS propagated errors in Z. Therefore, the vertical uncertainties of DSMs obtained in this study ranged from 0.23 m for the 2010 reference flight to 0.25-0.40 m for most flights and 0.55 m for the 2009 flight (Table 3).  [48]. These models allow the objective characterization of areas that show significant vertical changes in the ground surface or height differences between consecutive models. These height differences may be negative or positive, depending on whether each model lies below or above the preceding one, which would allow the identification of landslide areas with ground ascent (mass depletion) or descent (mass accumulation), respectively.
Unlike other studies [16,31], in this paper, DSMs were used instead of DTMs, because the models are of photogrammetric origin without stereo edition. Therefore, the processes of automatic classification and filtering of point clouds do not ensure accurate results in areas where there is a dense vegetation cover or other elements. Moreover, the stereo edition of models with a DPW to obtain DTMs in a wide area would be extremely timeconsuming. However, this produces a distortion in the landslide detection because some height differences can be caused by vegetation (growing, clearing, etc.), construction or other elements included in the DSMs, which requires some filtering operations, which are discussed in the following sections.
Regarding the vertical uncertainties of the DoDs, they were estimated as follows [11,27,32,49]: Thus, the uncertainty of the differential DSMs ranged from 0.35 to 0.70 m, being generally lower in the comparison of more recent flights and higher in the comparison of older flights (Table 3).

Filtering and Masking of DoDs
Once the differential models were obtained, it was necessary to identify the unstable areas or landslides, since not all the vertical changes measured in the DoDs were due to landslides. In fact, these were also due to vegetation changes, new constructions, maintenance work, etc., as mentioned previously. In addition, there was noise due to residual image misalignments and inaccuracies in DSM generation.
Thus, as the manual stereo edition of DSMs was discarded, other solutions were addressed based on the application of filters and masks in order to identify landslides from Remote Sens. 2021, 13, 2110 9 of 32 the DoDs. First, vegetation indexes, such as the normalized difference vegetation index (NDVI), were calculated from orthophotographs available in Spanish and Andalusian spatial data infrastructure [50]. These indexes allowed the exclusion of areas with dense vegetation by establishing threshold values. The NDVI [51] was calculated from the following expression:  [52], which could be calculated with the following expression: where R, G and B are the DN in the corresponding RGB bands. Finally, the older images corresponding to the 1984, 1996 and 2001 aerial surveys are panchromatic, so it was not possible to obtain vegetation indexes for these.
The index images were then classified by means of different thresholds, and several masks were obtained and applied to the DoDs. For the relative NDVI estimated, the threshold was fixed at 0.4, analyzing the images and their histograms in order to discriminate between vegetated and non-vegetated areas. This value allowed us to separate trees (pines, oaks, riparian, olives, etc.) from soil, water and urbanized areas, but also from brush and grass. Since all the flights were performed during the summertime, which is usually very dry in Southern Spain, these last covers present lower NDVI values. Meanwhile, for GLI, the threshold was fixed at 0.05, in accordance with the strategy described previously, which also separated trees from the other land covers. Finally, the panchromatic images did not allow us to establish a threshold to distinguish vegetated areas from other land cover, such as soils, urban and water. These masks, resulting from the binarization of NDVI and GLI images, were then filtered by means of a mode filter in order to obtain continuous masks. The window size of this filter was large enough (21 × 21 pixels, corresponding to 10.5 × 10.5 m) to filter the olive trees, which were integrated in the soil, but it allowed the exclusion of densely vegetated areas, especially those made up of trees (pines, oaks and riparian). Since the vegetated areas were changing from year to year and, in addition, there were no masks for some years (panchromatic images), a unique vegetation mask was generated. The final mask was obtained with the global extension of all the individual masks, and then filtered again with the same filter (21 × 21 mode). Afterwards the filter mask was resampled to a pixel size of 2.5 m, the same as the DoDs. The process of this vegetation filter is shown in Figure 4.

Elaboration of the Multitemporal Landslide Inventory
The final step in elaborating the landslide inventory and database was the photointerpretation and digitalization on orthophotographs and DoDs for each date and period. The most common method for identifying and delineating landslides is photointerpretation, both in 2D from orthophotographs and in 3D using analog stereoscopes or digital photogrammetric stations [3,4]. However, the process is relatively cumbersome, especially when a multitemporal approach is required, and subjective; thus, some movements frequently remain undetected to inexpert photointerpreters [4].
In this case, once the DoD calculation and filtering had been carried out, the results were checked using photointerpretation techniques in orthoimages and stereo pairs. This photointerpretation was performed considering the areas in which the estimated differences were higher than the estimated propagated uncertainty (1 m in absolute terms). This allowed for discrimination between landslides and other objects (DSM artefacts and misalignments), as well as true changes in the terrain itself, such as vegetation, constructions, excavations, fillings and quarries, not detected in the previous analysis. A useful criterion for identifying movements was the detection of adjacent areas of descent (or depletion) and ascent (or accumulation), the first area being topographically higher than the latter one [31].
Once the landslides had been identified, they were digitized in the GIS, obtaining different vector layers for each date. Then, different analyses were performed to characterize the movements and to complete the database with landslide typology, morphometric parameters, lithology, activity, etc. Landslide typology was assigned based on the height differences calculated in the DoDs and the length of the estimated movements, according to the well-known landslide classifications after Varnes [38,39], but it was also based on photointerpretation. Morphometric parameters, such as landslide areas (A), perimeters (P), height interval or range (H), length (L) and H/L ratio [54], were calculated by means of different tools (geometric, arithmetic, zonal, etc.) of QGIS software, based on the Then, a second mask with urban areas (populations, industrial, isolated buildings, roads, constructions, etc.) and water bodies (river and wetlands) was applied. These new masks were obtained from vector layers (shape files) downloaded from the Andalusian SDI and then rasterized to a resolution of 2.5 m and applied to the DoDs. Thus, urban and water areas, the same as the vegetated areas, were discarded as potential areas of landslides. This filter is also shown in Figure 4.
Third, a new convolution mean filter in a 5 × 5 pixel window (12.5 × 12.5 m) was applied again to the DoDs in an iterative manner in order to eliminate or reduce the noise caused by local residual misalignments in the orientation process, propagated to DSMs and DoDs, as well as the remains of the vegetation and other elements.
Finally, based on the uncertainty of the DoDs, another filter was applied in order to discard those areas with height differences with an absolute value lower than 1 m and also to reduce noise. In this sense, since, in this study, DoDs were only used for detecting unstable areas and estimating their activity but not for quantifying volumes, more complete approaches were not considered [49,53].

Elaboration of the Multitemporal Landslide Inventory
The final step in elaborating the landslide inventory and database was the photointerpretation and digitalization on orthophotographs and DoDs for each date and period. The most common method for identifying and delineating landslides is photointerpretation, both in 2D from orthophotographs and in 3D using analog stereoscopes or digital photogrammetric stations [3,4]. However, the process is relatively cumbersome, especially when a multitemporal approach is required, and subjective; thus, some movements frequently remain undetected to inexpert photointerpreters [4].
In this case, once the DoD calculation and filtering had been carried out, the results were checked using photointerpretation techniques in orthoimages and stereo pairs. This photointerpretation was performed considering the areas in which the estimated differ-ences were higher than the estimated propagated uncertainty (1 m in absolute terms). This allowed for discrimination between landslides and other objects (DSM artefacts and misalignments), as well as true changes in the terrain itself, such as vegetation, constructions, excavations, fillings and quarries, not detected in the previous analysis. A useful criterion for identifying movements was the detection of adjacent areas of descent (or depletion) and ascent (or accumulation), the first area being topographically higher than the latter one [31].
Once the landslides had been identified, they were digitized in the GIS, obtaining different vector layers for each date. Then, different analyses were performed to characterize the movements and to complete the database with landslide typology, morphometric parameters, lithology, activity, etc. Landslide typology was assigned based on the height differences calculated in the DoDs and the length of the estimated movements, according to the well-known landslide classifications after Varnes [38,39], but it was also based on photointerpretation. Morphometric parameters, such as landslide areas (A), perimeters (P), height interval or range (H), length (L) and H/L ratio [54], were calculated by means of different tools (geometric, arithmetic, zonal, etc.) of QGIS software, based on the inventory vector layers and the DSMs. In the same way, landslide conditions (average height, slope and orientation, as well as modal lithology) were also obtained in the GIS from the inventory, DSMs, the slope and orientation maps, as well as the lithological units derived from the geological map [40].
Regarding the activity, this could be estimated from the zonal analysis of the inventory on the DoDs. However, since the inventory is more oriented towards landslide characterization, another approach was used which is the analysis of monitoring areas. Thus, the activity of some representative examples of landslide monitoring areas with different characteristics (geological, topographical and landslide typology) was evaluated by means of the zonal analysis of filtered DoDs. This approach allowed the estimation of average height differences in the monitoring areas for the different periods, resulting in negative or positive balances that corresponded to general descents or ascents of the ground surface in these areas. Considering only the sectors of filtered DoDs with negative values or positive values, the average descents or ascents of the ground surface were also calculated. Dividing the height differences by the period length, the average vertical velocity could be estimated, although more accurate calculations of the velocity (vertical and horizontal) would require the measurement of displacements of points sampled in the unstable area [21,24]. Nevertheless, this estimated vertical velocity allowed the characterization of the activity of the landslide monitoring areas, according to published classifications [55,56] throughout the different periods, which could be related to the rainfall regime as the main triggering factor of landslides.

Results
The filtered DoD maps were the first results of this study from which a multitemporal landslide inventory was obtained. In this section, this inventory is presented in order to characterize and analyze the landslides of the study area. Moreover, a set of monitoring areas was also analyzed in order to determine the landslide activity of this area. Figure 5 shows the inventory differenced by periods, and Table 5 displays the corresponding analysis. Three basic landslide typologies were identified: falls-collapses, slides and flows [38,39]. In general terms, landslides cover 0.98 km 2 , representing 1.90% of the study area. In the first class, rock falls and collapses, only one rock fall in the strict sense was identified. Moreover, collapses in natural and engineering slopes were distinguished. Natural collapses extended for 0.46 km 2 (46.63% of the landslide area), collapses in engineering slopes extended for 0.16 km 2 (16.57%), slides extended for 0.24 km 2 (25.00%) and flows extended for 0.11 km 2 (11.78%). Average areas ranged from about 1000 m 2 of collapses (both natural and engineering slopes), to nearly 4500 m 2 of slides and 10,800 m 2   Regarding the average conditions in which landslides appear, the height was around 600 m in most cases and near 1000 m in the rock fall. The average slope was around 18 • in flows, while values of between 25 and 30 • occurred in slides and collapses (both in natural and engineering slopes), and 60 • in rock falls. Orientation is mostly towards the south and southeast, except collapses, which are mainly oriented to the west, and flows, which are oriented to the NE. Modal lithologies are Betic carbonates in the rock fall, Miocene marls in natural collapses, Betic marls from the Cretaceous age in slides and flows, and quaternary glacis and colluvial in collapses in engineering slopes. Finally, the average height differences between DSMs (DoDs) ranged from −0. All these values were negative, which implies a general descent of the ground surface.

Analysis of Monitoring Areas
As well as the multitemporal inventory, a set of common landslide areas throughout the whole period (1984-2016) was considered in order to monitor the landslide activity in the study area ( Figure 6). Table 6 shows the results of this analysis, which was mainly focused on the rate of height differences between DSMs, both the average of the whole areas and the average of sectors with negative and positive values.

Analysis of Monitoring Areas
As well as the multitemporal inventory, a set of common landslide areas throughout the whole period (1984-2016) was considered in order to monitor the landslide activity in the study area ( Figure 6). Table 5 shows the results of this analysis, which was mainly focused on the rate of height differences between DSMs, both the average of the whole areas and the average of sectors with negative and positive values.

Accuracy and Uncertainties
As previously explained, the uncertainty of the DSMs was estimated from the RMS Z (vertical) errors calculated at the GCPs used in the photogrammetric orientation. In previous studies [11,16,45], this uncertainty has been considered as two or three times the residual RMS Z errors at these points. Moreover, as the different photogrammetric flights had been oriented from second-order GCPs extracted from the reference flight (2010), there was an error propagation. Finally, the uncertainties of DoDs were estimated from the uncertainties of the pairs of flights compared. Thus, the obtained RMS Z errors ranged from 0.09 to 0.22 m with the DSM uncertainties between 0.25 and 0.55 m, and the DoD uncertainties between 0.35 and 0.70 m (Table 3). Thus, in a conservative manner, a minimum level of detection (minLoD) [32,49,53] was established at ± 1 m.
These values are comparable to those obtained in previous studies of the research group in which similar methodologies were applied in order to analyze the evolution of an active landslide [16] and a gully system [45], although with some differences depending on the available data. In both studies, photographs of the same dates and properties (national and regional flights) were used. In addition, they are of same order of magnitude as other studies by different authors [11,14].

Height Differences
The calculation of the height differences of the ground surface between successive dates allowed the observation of some general processes. First, there were some sectors where extensive changes occurred related to the growth or decrease in vegetation due to natural, agricultural or forestry activities. Thus, pine and holm oak forests could grow or be cleared, and scrubs or shrubs could suffer changes, although in the last case, their effects were less significant. In the same way, olive groves-the dominant crop in the study areaunderwent more or less significant modifications throughout the study period, although these are difficult to detect as they are single isolated trees separated by several meters from each other. In any case, most of these changes were discarded by means of the application of an NDVI filter to remove extensive areas and, in some cases, with the help of mode or mean filters to remove the effect of isolated trees (Figure 4).
In addition to vegetation, urban areas (populations, industrial, isolated buildings, roads, constructions, etc.) and water bodies (river and wetlands) were removed by applying masks from vector layers (Figure 4), since in these areas changes not related to the ground surface could occur. However, some changes in the ground surface in the surroundings of urban or construction areas, such as human-made excavations, explorations or infills, had remained. Among them, the earthworks related to the construction of the A-44 highway in the mid-1990s could be observed as well as other minor works of cutting, grading and fillings related to civil engineering and agricultural activities. These areas could have been removed by means of the application of corridors (buffers) along these constructions, but in this case, we preferred to keep them because some landslides originated in the cut slopes after the construction works. The effect of roads and construction was considered and then discarded for the landslide inventory by means of photointerpretation. Thus, the highway roadworks appeared clearly in the first period of 1984-1996 and to a lesser extent in the second period of 1996-2001, while other minor works could be observed in these and other intervals.
Once the effects of previous elements (vegetation, urban areas, water bodies and earthworks) had been removed, the height differences observed in the maps (Figures 7-9) should correspond to changes in the ground surface, although some isolated vegetated and urban areas could remain, as well as some artifacts and misalignments. These last elements were also discarded by photointerpretation. Thus, natural changes in the ground surface, such as landslides, erosion processes and subsidence, could be detected. The interpretation of these changes as landslides was undertaken from orthoimages but also from DoD analysis. For example, the observation of adjacent depletion and deposition areas from DoDs could be related to landslides (Figures 7-9). When discarding general subsidence processes in the area, gully erosion or small karstic subsidence processes in relation to gypsum bodies could occur. Erosion processes were well identified [45], and thus the depletion and accumulation areas in the bottom and lateral walls of the gullies were not considered as landslides, although some small landslides (mainly collapses) in the lateral walls were inventoried. Small karstic subsidence cavities were also well identified in the orthoimages, and thus they were discarded as landslides.

Landslide Inventory and Factor Analysis
As previously mentioned, three basic typologies were distinguished in the area: fallscollapses, slides and flows. Practically no falls in the strict sense were identified in the study period , but there were many examples of collapses that occurred in the lateral walls of riverbeds and gullies, as well as in engineering slopes, which were considered in this study as separate classes. Meanwhile, the slides were mostly earth slides, but they were always shallow movements of a few meters in depth. Flows were mostly earth flows, although in some cases, they may have been in transition to mud flows; nevertheless, they were also movements with a depth of a few meters and a length of hectometric magnitude. Deep slides or large flows that appear in this area [35] and similar ones [57][58][59] have not presented any evidence of activity in recent years. The same occurs with debris flows that also appear in the surrounding ranges but with no activity in recent years in the study area.
The percentage of mobilized area, near to 2%, was higher than other values found for recent landslides (active or suspended, according to [56]) in previous studies in a wider area that includes this one (1.15 % in [35]) and in other areas with similar conditions (0.56% in [58,59]). This may be due to the fact that this is a more active area with natural and artificial conditions prone to landslides. Such conditions are basically high slopes in an active geological front with marly and clayey lithologies, the excavation of a river and the construction of a highway. However, it could also be due to the elaboration of a more precise landslide inventory by means of DoD calculation and photointerpretation. Nevertheless, these values were lower than those calculated in landslide inventories, which also include ancient (dormant and relict) movements, in similar environments to the Betic Cordilleras [35,[57][58][59], where percentages between 5 and 10% have been reached. Analysis by typologies showed an area distribution ( Figure 10) with a higher proportion of collapses (63%) and slides (25%) than flows (12%). In the extended area, which includes ancient landslides, the percentage of rock falls reached 10%, and flows increased to 55%, while slides decreased to 35% [35]. If we consider the number of landslides, the proportion of collapses in natural and engineering slopes increased by nearly 85%, while slides and flows decreased to 12 and 3%, respectively, due to the small individual area of collapses (about 1000 m 2 ), with respect to slides and flows (4500 and 10,800 m 2 , respectively) ( Figure 10). The analysis of determinant factors and the conditions associated with landslides is usually an initial step in susceptibility modelling [59,60]. These conditions are different for each typology since they correspond to different mechanisms [38,39,59]. Thus, the rock fall appeared in a steep slope (average of 60°) of the higher part of the area, where limestones of the Betic ranges outcropped. Collapses were concentrated in a narrow range of height between 600 and 630 m, with average slopes of around 30°, excavated mostly in Miocene and Cretaceous marls and also in Pliocene-Quaternary deposits (conglomerates and sands). Slides were associated with slopes of between 20 and 30° formed mostly in Cretaceous and Neogene marls and Triassic lutites. Flows originated in slopes of between 10 and 25°, oriented mainly to N and NE, and formed mostly in Cretaceous-Neogene marls and clays. These conditions are very similar to those in which the different landslides occur in other areas of the Betic Cordilleras [35,[57][58][59], although the slopes are somewhat steeper in this area.
Finally, the average height differences (DoDs) corresponded to average descents ranging from very low values (−0.13 m) in flows to values near −1 m in slides, −1.2 to −1.6 m in collapses and −2.2 m in rock falls. This seems consistent with the kinematics and evolution of landslides observed in other studies [19,27,31]. Usually, greater descents of the ground surface are observed in the landslide head or scarp, although more concentrated in space, while lower ascents are observed at the foot area, although these can be more widely distributed [61]. This causes a higher proportion of the ascents to be below the uncertainty threshold with respect to the descents, which leads to an average surface descent in the balance. Moreover, a proportion of the deposit masses accumulated at the foot of the slope can be eroded and evacuated by the drainage network [16,21]. This effect also depends on the landslide typology. Thus, the movements with steep scarps, such as collapses, show significant descents of the ground surface in the landslide head and much lower ascents (in many cases under the uncertainty threshold) in the foot or accumulation zone. Meanwhile, flows undergo small descents in the head and even smaller ascents in the foot, since the horizontal displacements predominate over the vertical ones. Regarding the morphometry, the height interval was around 40 m in the rock fall that produced a high H/L ratio of 2.23, typical of a small movement with a high potential activity [54]. In collapses, the height interval had an average of 14-19 m for the natural and engineering slopes, and the H/D ratios were in the range 0.54-0.60, respectively. These relatively high values were also related to small movements with a moderate potential activity. Meanwhile, the height interval associated with slides and flows was about 35-45 m, with larger areas and lower H/L ratios (0.44-0.19), suggesting a lower level of potential activity.
The analysis of determinant factors and the conditions associated with landslides is usually an initial step in susceptibility modelling [59,60]. These conditions are different for each typology since they correspond to different mechanisms [38,39,59]. Thus, the rock fall appeared in a steep slope (average of 60 • ) of the higher part of the area, where limestones of the Betic ranges outcropped. Collapses were concentrated in a narrow range of height between 600 and 630 m, with average slopes of around 30 • , excavated mostly in Miocene and Cretaceous marls and also in Pliocene-Quaternary deposits (conglomerates and sands). Slides were associated with slopes of between 20 and 30 • formed mostly in Cretaceous and Neogene marls and Triassic lutites. Flows originated in slopes of between 10 and 25 • , oriented mainly to N and NE, and formed mostly in Cretaceous-Neogene marls and clays. These conditions are very similar to those in which the different landslides occur in other areas of the Betic Cordilleras [35,[57][58][59], although the slopes are somewhat steeper in this area.
Finally, the average height differences (DoDs) corresponded to average descents ranging from very low values (−0.13 m) in flows to values near −1 m in slides, −1.2 to −1.6 m in collapses and −2.2 m in rock falls. This seems consistent with the kinematics and evolution of landslides observed in other studies [19,27,31]. Usually, greater descents of the ground surface are observed in the landslide head or scarp, although more concentrated in space, while lower ascents are observed at the foot area, although these can be more widely distributed [61]. This causes a higher proportion of the ascents to be below the uncertainty threshold with respect to the descents, which leads to an average surface descent in the balance. Moreover, a proportion of the deposit masses accumulated at the foot of the slope can be eroded and evacuated by the drainage network [16,21]. This effect also depends on the landslide typology. Thus, the movements with steep scarps, such as collapses, show significant descents of the ground surface in the landslide head and much lower ascents (in many cases under the uncertainty threshold) in the foot or accumulation zone. Meanwhile, flows undergo small descents in the head and even smaller ascents in the foot, since the horizontal displacements predominate over the vertical ones.

Multitemporal Inventory Analysis
The analysis of multitemporal inventories is suitable for landslide hazard modelling [3,62]. In this case, it led us to observe changes in activity (Figure 11

Multitemporal Inventory Analysis
The analysis of multitemporal inventories is suitable for landslide hazard modelling [3,62]. In this case, it led us to observe changes in activity (Figure 11 However, if we consider the number of landslides, the proportion was more uniform between the different periods, about 12% for slides, 3% for flows and more than 80% for collapses. Regarding individual areas and morphometric parameters, collapses maintained relatively uniform areas and height intervals around the aforementioned average values (1000 m 2 and 15-20 m, respectively). Slides showed some differences in the areas (2000-8000 m 2 ) and height intervals (25-60 m) but without a clear pattern between more or less active periods. Meanwhile, flows showed a certain tendency between less active periods (areas lower than 8000 m 2 and height intervals lower than 35 m) and the more active periods (larger than 10000 m 2 and greater than 40 m, respectively). However, the largest dimensions were found in 2013-2016 (25000 m 2 and 66 m). This is due to the fact that some flows of larger dimensions, originated or reactivated in the previous periods (2009-2013), were still active in the first months of this period before their stabilization [21]. Meanwhile, the H/L ratio was also uniform throughout the different periods for all types of landslides.
Most conditions in which the different typologies of landslides occurred, such as height, slope and geology, did not change significantly over the different periods, since these conditions are more related to the typology than the activity. Nevertheless, when we analyze all the movements, we can find some differences related to the different proportions of landslide typologies throughout the different periods. Thus, less active periods with a higher proportion of collapses presented greater slopes (31-33°) than the most active (about 28-29°). In the same way, the predominant lithology in periods with a high proportion of collapses in engineering slopes (1996-2001 and 2001-2005) were the Quaternary deposits, while in the remaining periods, the Miocene marls predominated.
Finally, the average height differences between DSMs were always negative (ground descent), but they showed greater absolute values in the periods of high activity in landslide number and area. Thus, in the less active periods (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(2001)(2002)(2003)(2004)(2005) and 2013- Regarding individual areas and morphometric parameters, collapses maintained relatively uniform areas and height intervals around the aforementioned average values (1000 m 2 and 15-20 m, respectively). Slides showed some differences in the areas (2000-8000 m 2 ) and height intervals (25-60 m) but without a clear pattern between more or less active periods. Meanwhile, flows showed a certain tendency between less active periods (areas lower than 8000 m 2 and height intervals lower than 35 m) and the more active periods (larger than 10,000 m 2 and greater than 40 m, respectively). However, the largest dimensions were found in 2013-2016 (25,000 m 2 and 66 m). This is due to the fact that some flows of larger dimensions, originated or reactivated in the previous periods (2009-2013), were still active in the first months of this period before their stabilization [21]. Meanwhile, the H/L ratio was also uniform throughout the different periods for all types of landslides.
Most conditions in which the different typologies of landslides occurred, such as height, slope and geology, did not change significantly over the different periods, since these conditions are more related to the typology than the activity. Nevertheless, when we analyze all the movements, we can find some differences related to the different proportions of landslide typologies throughout the different periods. Thus, less active periods with a higher proportion of collapses presented greater slopes (31)(32)(33) 28-29 • ). In the same way, the predominant lithology in periods with a high proportion of collapses in engineering slopes (1996-2001 and 2001-2005) were the Quaternary deposits, while in the remaining periods, the Miocene marls predominated.
Finally, the average height differences between DSMs were always negative (ground descent), but they showed greater absolute values in the periods of high activity in landslide number and area. Thus, in the less active periods (1984-1996, 2001-2005 and 2013-2016), the height differences were between −0.5 and −0.7 m, while in the more active periods (1996-2001, 2009-2010, 2010-2011 and 2011-2013), the height differences were about −1 m. By typologies, collapses remained relatively uniform at values of between −1 and −2 m. However, slides presented greater descents in the more active periods (more than 1.0 in absolute value for 1996-2001 and 2009-2010) than in the less active ones (under 0.5 m for 1984-1996, 2001-2005, 2005-2009 and 2013-2016), which agrees with observations of individual slides studied previously [16]. Flows presented very low and even null average values, but they reached significant average descents (between −0.2 and −0.4 m) in those periods of more activity (1996-2001, 2009-2010 and 2011-2013).
In summary, in the periods of higher activity, the movements were more abundant; the area affected was also larger; some typologies, such as flows, had greater dimensions; and, in particular, they presented a greater average descent of the ground surface. Meanwhile, in the periods with lower activity, fewer movements were originated, the affected areas were smaller and the average descents of the ground surface were lower and even null. The exception is the period 2013-2016, which was of general low activity, in which some large flows and slides originating in the previous periods continued their activity [21].

Analysis of the Monitoring Areas
This analysis allowed the monitoring of a set of common landslide areas throughout the whole period, which can be used in further landslide hazard modelling [62]. Regarding the morphometry, in areas 1, 6, 7, 8, 9 and 10 (Figures 7 and 9), where collapses predominate, the dimensions are around 2800 m 2 (area) and 30-50 m (height interval), somewhat greater than the average of collapses in the inventories, since in these areas there are also slides and flows. In areas 2 and 5, where slides predominate, the dimensions are around 9800 m 2 and 30-55 m, respectively, which are similar to those found in the inventories. The changes observed in some areas in the period 1984-1996, when the relatively higher descents were compensated with ascents, were mainly related to earthworks of the A44 highway. Nevertheless, these rates allowed us to catalogue the landslides as very slow (0.016-1.6 m/year) in average terms throughout the periods considered, although they could probably reach the category of slow (1.6 m/year-13 m/month) or even moderate (13 m/month-1.8 m/hour), according to the classification of WP/WLI [39,55].  1984-1996 and 1996-2001, related to the road works of the A4 highway. In all these periods, the descent rates were greater than 1 m/year, and the ascent rates were lower than 1 m/year, in absolute terms. The average velocities for the collapses were in the boundary of very slow to slow, although this average velocity is not significant since this type of movement usually reaches rapid to very rapid velocity (more than 1.8 m/hour).
Areas 2 and 5, in which slides predominate, showed significant descent rates in the period 2009-2010 (about 1.6 m/year) and in 2009-2013 and 1996-2001 (lower than 0.5 m/year), with some minor ascent rates in some cases that produced a general ground descent in these periods. With these average rates, the movements were catalogued as slow, although they could reach moderate and even rapid velocity. Meanwhile, areas 3 and 4, with two significant flows studied previously by means of UAS surveys [21,36], presented sectors with significant descent and ascent rates in several years. The main changes took place in 2009-2010, when significant descent rates (1-3 m/year) and ascent rates (0.5-1.5 m/year) were estimated to result in a general negative balance, although this also occurred in other active periods, such as 2010-2011 and 2011-2013. The average velocities were in the boundary of very slow to slow, but more precise measurements, in both the vertical and the horizontal component, allowed us to catalogue the flows as slow in the more active episodes. Moreover, one of these flows (area 4; Figure 8) maintained its downslope movement until the first year of the period 2013-2016 [21]. Meanwhile, the other flow (area 3; Figure 8) had a previous activity in the active period 1996-2001, with significant descents in the upper part of the slope, after an ascent in the same area in the previous period (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), which could be due to a fill of material derived from nearby quarries and works.
In summary, some periods were more active, such as the consecutive 2009-2010, 2010-2011 and 2011-2013 periods, especially the first, and previously the period 1996- In all these periods, the descent rates were greater than 1 m/year, and the ascent rates were lower than 1 m/year, in absolute terms. The average velocities for the collapses were in the boundary of very slow to slow, although this average velocity is not significant since this type of movement usually reaches rapid to very rapid velocity (more than 1.8 m/hour).
Areas 2 and 5, in which slides predominate, showed significant descent rates in the period 2009-2010 (about 1.6 m/year) and in 2009-2013 and 1996-2001 (lower than 0.5 m/year), with some minor ascent rates in some cases that produced a general ground descent in these periods. With these average rates, the movements were catalogued as slow, although they could reach moderate and even rapid velocity. Meanwhile, areas 3 and 4, with two significant flows studied previously by means of UAS surveys [21,36], presented sectors with significant descent and ascent rates in several years. The main changes took place in 2009-2010, when significant descent rates (1-3 m/year) and ascent rates (0.5-1.5 m/year) were estimated to result in a general negative balance, although this also occurred in other active periods, such as 2010-2011 and 2011-2013. The average velocities were in the boundary of very slow to slow, but more precise measurements, in both the vertical and the horizontal component, allowed us to catalogue the flows as slow in the more active episodes. Moreover, one of these flows (area 4; Figure 8) maintained its downslope movement until the first year of the period 2013-2016 [21]. Meanwhile, the other flow (area 3; Figure 8) had a previous activity in the active period 1996-2001, with significant descents in the upper part of the slope, after an ascent in the same area in the previous period (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), which could be due to a fill of material derived from nearby quarries and works.
In summary, some periods were more active, such as

Relationships with Rainfall
The relationships between rainfall and landslides have been well established in different environments and regions. In fact, rainfalls are the main triggering mechanism of landslides all over the world [63,64], but also in Mediterranean countries [65] and particularly in the regions near the study area [37,59].
Some of the movements, usually those with lower dimensions, such as collapses and smaller slides, presented an episodic activity related to intense rainy events, while others, such as larger slides and flows, presented an intermittent activity with periods of lower and higher rates of movement or velocity. This is typical of landslides in Mediterranean or arid climates [65], where the rainfall regime shows an irregular distribution alternating between dry and wet years [66]. Figure 13. Relationships between landslide activity and rainfalls: (a) 1-day rainfalls and landslide number; (b) 1-week rainfalls and total area (km 2 ); (c) 1-month rainfalls and average height differences. Figure 13. Relationships between landslide activity and rainfalls: (a) 1-day rainfalls and landslide number; (b) 1-week rainfalls and total area (km 2 ); (c) 1-month rainfalls and average height differences.

Conclusions
The techniques of photogrammetry and LiDAR have been revealed as very useful tools for landslide identification and the characterization of their activity. LiDAR and aerial photogrammetry allow us to work with individual movements of a certain magnitude (hectometer to kilometer), but also with areas of up to several kilometers, as with the studied case (about 50 km 2 ), in order to elaborate multitemporal inventories.
From historical and recent flights oriented through direct orientation procedures and iterative adjustment between models, we were able to obtain the corresponding DSMs and orthophotographs, and especially the differential models (DoDs), which showed potential unstable areas with changes in the ground surface. After a filtering process, noise was partially eliminated as well as vegetation, construction areas and water bodies. Furthermore, a threshold based on the uncertainty of DoDs of about ±1 m was considered. Then, by means of photointerpretation in addition to DoDs analysis, landslides were identified and digitized for the different periods. GIS tools allowed the inventory and factor analysis to characterize the landslides dimensions (area, perimeter, height interval and H/L ratio) and their properties (height, slope, orientation and lithology). Moreover, the zonal analysis of DoDs in the different periods considered led us to estimate vertical velocity as an expression of landslide activity. Nevertheless, this analysis was complemented with other analyses focused on a common set of landslide monitoring areas, in which the activity was monitored in a more specific way.
The inventory first shows a more extensive and detailed inventory than that obtained only with photointerpretation for recent movements in other areas with similar characteristics to the Betic Cordilleras, although it is less extensive than those inventories that include ancient movements, as can be also observed in other areas of the region. Thus, the proportion of small landslides, such as collapses and slides, in this area is higher than large slides and flows, in terms of affected area but especially in landslide number. However, the conditions in which landslides occur are practically the same, that is, great slopes (30 •  Both the multitemporal inventory elaborated in this study and the consequent factor and temporal analyses constitute the starting point of susceptibility and hazard modelling. Currently, a great variety and quantity of landslide hazard approaches are available [60,62], but they must be based on accurate multitemporal inventories [2,3].
Future research will focus on adding more flights, both those corresponding to official plans, as they become available in public SDI, and those carried out ad hoc, such as UAS flights, in order to perform a more complete monitoring of the activity of the area. LiDAR data should be also incorporated into the analyses since they allow us to obtain more accurate DSMs and DTMs. Meanwhile, filtering techniques must be improved to perform better discrimination between areas with changes in the DSMs and areas where landslides have actually occurred. Machine and deep learning techniques should also be implemented for the automatic detection of landslides and other processes, such as erosion and subsidence.