LANDSLIDE EVOLUTION PATTERN REVEALED BY MULTI-TEMPORAL DSMS OBTAINED FROM HISTORICAL AERIAL IMAGES

Landslides are a widespread natural hazard that cause damages to people and to the built up environment. Accurate knowledge of landslide distribution is crucial to develop planning strategies, prevention and resilient communities worldwide. One of the most diffuse way of reporting landslides distribution in a territory is by preparing landslide inventory maps. Such a task is mostly accomplished by expert photo-interpretation of historical aerial photographs, which are an invaluable source of information because they portray the landscape when the anthropic pressure was lower than the present day, providing an observation of the landscape closer to the natural conditions. Despite such a common use of aerial photographs, they are poorly exploited to obtain quantitative measures to support landslide mapping activities. In this paper we present a comparison of two photogrammetric approaches to measure elevation changes in a 50-years period for an area densely affected by landslides in Southern Italy. The obtained results allowed to revisit the original expert mapping proving that such a method is a useful tool to support geomorphological mapping and to improve the overall accuracy of landslide inventories.


INTRODUCTION
Landslides represent a major natural hazard. They are widespread phenomena that occur in many parts of the World, under different morphological, geological, climatic conditions, often causing major damages to structures, infrastructures and, most importantly, direct damage to the population (deaths, injuries and missing persons), (Alvioli et al., 2021;Donnini et al., 2017;Santangelo et al., 2020). Landslides can involve different types of material and the types of movement can be very different (Hungr et al., 2014). They can be very fast, such as rockfall or debris-flows, or slow-moving, such as earthflows. Depending on the type of failure, landslides leave discernible signs on the topographic surface and on land use, which are exploited by geomorphologists to prepare landslide inventory maps. Such features are referred to as morphologic and radiometric signature (Fiorucci et al., 2018;Guzzetti et al., 2012;Santangelo et al., 2022).
Studying the spatial and temporal pattern of landslides is a complex task usually carried out by geomorphologists by preparing landslide inventories, the simplest geographical representation of landslides distribution in a territory . Different types of inventories exist, which portray different characteristics of landslides that can be detected and mapped. The simplest form of inventory is an event landslide inventory map Ardizzone et al., 2012;Guzzetti et al., 2012). It reports evidences of landslides occurred after a specific triggering event (e.g. earthquakes, heavy rainfall, rapid snowmelt, volcanic eruptions). Geomorphological landslide inventories (Bucci et al., 2016;Guzzetti et al., 2012;Santangelo et al., 2014) report the spatial distribution of landslides occurred in the last tens to tens of thousands of years in an area, and can be considered as the results of multiple events (Malamud et al., 2004). The most complex type of inventory is a multi-temporal inventory, which reports all the evidences about * Corresponding author event landslides occurred in the last tens of years as recorded by the available aerial images (Galli et al., 2008;Guzzetti et al., 2012).
Even though inventories provide information about landslides distribution in time and space, little is done to fully exploit the quantitative information available from aerial photographs (Seccaroni et al., 2018). The increasing availability of historical aerial images and of photogrammetric techniques (Structure from Motion, SfM, and Multi-View Stereo, MVS) implemented in open source software (Rupnik et al., 2017) offer an unprecedented opportunity for landslide scientists to corroborate heuristic interpretation by photogrammetric data similarly to what already done for earthquakes (Ayoub et al., 2009).
In this work, we have adopted two photogrammetric processing pipelines, a standard manual procedure and an automatic procedure developed for multi-epoch historical images (i.e. images acquired in different years) (Zhang et al., 2021) to obtain Digital Surface Models (DSMs) of a test area in South Italy where a slow moving slide-earthflow has been detected since the 1950s. We have generated two DSMs for each available set of aerial photographs to compare their performance in measuring morphological changes. We also compare the results of the DSMs differences through time to landslide maps prepared independently by expert image interpreters. Such maps of the morphological changes of landslides over time could help geomorphologists prepare more accurate landslide maps and also be a precious tool to study landslide dynamics and evaluate hazard and risk.

DATA AND METHODS
Available images: We have used two sets of black and white stereoscopic aerial photographs ( Figure 1A), taken in September 1954 by the Gruppo Aereo Italiano, GAI, ( Figure 1B) and in May 2003 by the Istituto Geografico Militare Italiano, IGMI ( Figure  1C). The main characteristics of the two sets of aerial photographs are summarised in Table 1. Aerial photographs were provided as 800 dpi scanned images by the IGMI service. The images provided in digital version were scanned from poorly preserved images, which present some scratches and dust, a quite typical feature for images that were not preserved for photogrammetric purposes. Also, the scanning technique adopted is non-photogrammetric, i.e. using a commercial flat scanner. Unfortunately, no information about the scanner model is available.  Table 1. Main characteristics of aerial photographs used.

Available landslide data:
Landslide data were collected by expert photo-interpreters using both 1954 and 2003 aerial photographs following the approach described by Santangelo et al., (2015Santangelo et al., ( , 2014. Landslides were classified according to the type of movement and relative age according to the morphological and radiometric evidences of the aerial photographs examined. In particular, four generations of landslides were mapped in the study area (Figure 2), the first three being antecedent to the first aerial photographs acquisition (pre-1954), and the latest belonging to a landslide event occurred a few months before the acquisition of the 2003 aerial image, presumably during the winter 2002 -early spring 2003.
Photogrammetric pipelines: The images were processed using two different multi-epoch procedures. The first one, henceforth referred to as manual pipeline or manual procedure, includes the following steps: 1. Intra-epoch processing: For images from the same epoch, we (i) apply an affine transformation of the original images to the geometry of the fiducial marks, (ii) apply feature matching based on SIFT, (iii) compute relative orientations based on the sequential SfM; 2. Inter-epoch processing: We (i) perform a manual identification of a denser network of common feature correspondences then acting as virtual GCPs in the next step; (ii) compute a rough co-registration of all epochs using Ground Control Points and a 7-parameter transformation; 3. Combined processing: Based on all GCPs we perform a joint bundle adjustment to refine all the image orientations and camera calibrations; 4. DSM computation: We computed the DSMs, for each epoch individually.
The second procedure, henceforth referred to as automatic pipeline or automatic procedure, was presented by Zhang et al., (2021)  get DSMs in their arbitrary coordinate frames. Steps (i) to (iii) are the same for both the manual and automatic pipelines; 2. Inter-epoch processing: We (i) roughly co-register the DSMs and image orientations from different epochs and (ii) perform precise matching under the guidance of the roughly co-registered results to reduce the matching ambiguity; 3. Combined processing: Based on the intra-epoch and interepoch feature correspondences, we perform a joint bundle adjustment to refine all the image orientations and camera calibrations; 4. Geo-referencing and DSM computation: Using a minimum of 3 Ground Control Points (GCPs) and a 7-parameter transformation we transformed the result from an arbitrary to an absolute frame. Finally, we computed the DSMs, for each epoch individually.
Note that steps 1-3 are fully automated and seamlessly embedded within the "TiePHistoP" command in MicMac.
We then compared the Difference of DSMs (DoD) obtained after application of the two pipelines and analyse the elevation absolute differences using GRASS GIS 7.8 (GRASS Development Team, 2020). Best results are then used to estimate the landslide evolution in the test site.

STUDY AREA
The study area is located in Southern Italy, near the village of Alberona, in Puglia region ( Figure 1A). In the area elevation ranges between 346 and 1000 m a.s.l. (µ=559 m, σ=144 m), whereas slope spans between 0 and 70° (µ=11.6°, σ=7.7°). This portion of the Italian Apennines is characterized by a typical smooth topography, due to the diffuse presence of clay rich lithologies (Bucci et al., 2022, preprint), which also is a predisposing factor for widespread slow moving landslides, mainly of the slide and slide-earthflow type. In the study area, the land use is typical of rural areas (Zumpano et al., 2020), sparsely inhabited and mainly agricultural and wooded. The vegetation cover has experienced moderate wood growth over the years covered by our analysis.
The Alberona test case: The selected test site was affected by a landslide occurred during the winter/spring season 2002-2003 as the result of heavy rainfalls that hit this part of northern Puglia Region. The geomorphological inventory map prepared by expert photo-interpretation highlighted that the landslide occurred on a slope that had already been affected by large slideearthflows in its past. In particular, three generations of landslides were recognised and mapped as shown in Figure 2. The largest landslide that was recognised is also the oldest, and can be thought as the first generation failure. It is a slideearthflow 2,880 m long, and 571 m wide in its widest point ( lowest difference values in the centre of the image (i.e. dome effect), even though the procedure included a joint bundle adjustment optimisation of images belonging to the two epochs. Such an effect is absent in the DoD computed using the automatic pipeline. Furthermore, it is worth noting that in Figure 3B elevation differences change the sign along some E-W oriented borders of the aerial photographs. Such artefacts are due to the rather poor overlap that characterised the acquisitions of aerial photographs and are difficult to remove.
Despite the presence of some artefacts in the DoD computed by the automatic procedure, it is more credible than the DoD computed by the manual approach, which produced a DoD characterised by a dome effect and by similar artefacts as the automatic procedure.

Accuracy evaluation of the automatic pipeline:
The automatic pipeline provides satisfactory internal and external camera calibrations (Table 2), expressed in an arbitrary, but common, reference frame for all epochs. As it can be seen in Figure 3, the dome effect that appeared in the result of manual procedure is effectively alleviated in the DoD of automated pipeline, because the camera parameters are better estimated thanks to the numerous and precise inter-epoch feature correspondences. However, we still observe residual systematic errors as shown in Table 2, which are difficult to remove, because (i) the images are poorly preserved and scanned with nonphotogrammetric scanners, and (ii) limited number of images leads to a lack of redundant observations.
Year Bundle adjustment σ0 [pix] The evaluation of the accuracy of the DoD obtained by applying the automatic pipeline is not an easy task due to the absence of a quantitative independent measure to compare against. We addressed this problem by measuring the values of the DoD in areas where they are supposed to be null or negligible, i.e. in areas (i) where the terrain is supposed to be stable, (ii) outside of wooded and built up areas, and (iii) where no relevant landuse change have occurred. Note that the working assumption of limiting to stable areas is not necessarily true, since landslides may have occurred during the 50-years interval considered in this work and their evidences may have been cancelled by, e.g. anthropic activity and erosion. To prevent similar issues we selected areas characterised by a gentle slope, better if flat.   and that a Gaussian distribution approximates the error values, the 95% confidence interval corresponds to ± 2.92 m.

Comparison with landslide mapping:
The DoD obtained by the automatic procedure was then compared to each generation of landslides mapped in the geomorphological inventory. Results of the comparison are shown in Figure 5. Inspection of the Figure reveals that: 1. The first generation failure (thick black polygons in Figure  5A) does not show evidences of a complete remobilization in the observation period, as shown by small values of elevation difference in the scarp and deposit areas.
2. The second generation failure (thick blue polygons in Figure  5B), well evident in the photographs, shows a relevant deformation occurred in the deposit area during the observation period. Whereas the changes measured in the scarp and transit area are due to the 2003 event, the metric deformation measured in the deposit zone cannot be completely attributed to the single event but to a combined effect of a cumulate deformation during the decades and the effect of the event itself. This interpretation is based on the observation that the radiometric signature of this portion of the landslide would have shown clear evidences if it was entirely caused by the 2003 event. It is worth noting that on the 2003 image, the geomorphologists did not report any changes in this part of the landslide, so the changes were hardly detectable even to expert photo-interpreters. 3. Based on the features shown in the aerial photographs of 1954, the third generation landslide (thick yellow polygons in Figure 5C) occurred likely a few seasons (years) before the acquisition of the 1954 photograph. Measures of change obtained by the automatic pipeline cluster at the toe, and the arc-shaped pattern is compatible with the type of movement (i.e., flow) of the event landslide occurred in 2003. 4. The 2003 event landslide (thick red polygons in Figure 5D) is characterised by an alternation of positive and negative differences that is common in clay-rich earthflows (Giordan et al., 2013). Based on evidences in the DoD, the interpreters The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France considered that a more consistent delineation of the 2003 event landslide should consider the arc-shaped deformations measured at the toe of the third generation landslide (dotted red line in Figure 5D). Considering the amended (red dotted) polygon, the total area affected by the landslide is 256,918 m 2 , which represents a 57.4% increase compared to the original delineation.
Finally, Figure 6 shows a longitudinal profile of the two DSMs approximately along the central line of the 2003 landslide. Inspection of the longitudinal profile reveals that the landslide scarp was drawn correctly by the interpreters as the landslide head. Furthermore, the measures show that the escarpment area of the 2003 landslide produced a drop of about 20 metres in the slope profile and that most of the deformation is located in the upper part of the slope, whereas smaller differences are located along the transit and deposit zones, as also evident in Figure 5D.

Discussion:
The task of preparing landslide inventory maps is usually accomplished by applying photointerpretation criteria to remote sensing images such as stereoscopic aerial photographs, satellite images or even Digital Elevation Models (Fiorucci et al., 2011;Niculiţă et al., 2016;Santangelo et al., 2015).
Analyses carried out by photo-interpreters are heuristic, based on the third dimension perception of the interpreters and according to a complex multivariate analysis based on experience (Mondini et al., 2019;Santangelo et al., 2022). Given the heuristic nature of photointerpretation, several sources of error are implicit within the mapping activity Santangelo et al., 2015). A way of reducing inconsistencies in the application of mapping criteria and the systematic use of reproducible criteria and use of ancillary information. The experiment presented here introduces a further instrument that could be exploited by geomorphologists to support their mapping activity through quantitative analyses. Such data could provide information where geomorphological and radiometric signature of landslides is subtle and it can also lead to a better estimation of landslide area, which in turn is precious information for quantitative risk assessment. In our experiment, the automatic pipeline (Zhang et al., 2021) proved more effective than the manual procedure in providing a good co-registration between DSMs belonging to two different epochs, which resulted in a 1.1 m resolution DoD with an error of ± 2.95 m at 95% confidence interval. Despite it may seem a poor value, it must be considered that the images used were scanned from poorly preserved printed copies at 800 dpi using non photogrammetric commercial scanners. To this extent, the benefit from using a photogrammetric scanner, scanning the original film instead of printed copies and using a higher resolution remains unexplored and further research is needed. We emphasize here that the method used for quantifying the error associated with our measures is valid but is not as accurate and reliable as a quantitative independent measure to use as benchmark. Apart from the assumptions already explained in Section §4.1.2, three main aspect deserve specific mention: (i) the statistical representativeness of the selected points for the overall distribution of the area where the expected value of the DoD is 0, (ii) the representativeness of the measured error for areas where the values of the DoD are not expected to be null, and (iii) the effect of Ground Control Points on the values of errors in their surroundings. These three topics need to be specifically addressed to provide more accurate estimations of the errors associated with the methods proposed. Despite these issues related to the estimation of the error, we maintain that our results represent the dynamics of the studied landslide with unprecedented detail. Our results provided expert geomorphologists with new data to support and locally revisit their map. The new delineation of the event landslide includes an area 57% larger than the original map, which has impact on the estimation of landslide volume and magnitude . The data obtained in this experiment may also provide useful information for the preparation of multi-temporal landslide inventories and may be integrated in finite element models. Furthermore, accurate estimations of landslide volumes mobilized in specific time intervals may help studies on landslide mobilization rates and defining the impact of landslides on structure and infrastructures. Finally, analyses like the one presented here may be the basis for a long time baseline monitoring of landslides that may be coupled with shorter time monitiring systems, e.g. based on SAR interferometry.

CONCLUSIONS
In this paper we presented a comparison between two photogrammetric pipelines applied to an area densely affected by slow-moving landslides in Puglia Region, in Southern Italy. The standard manual procedure based on a two-phases bundle block adjustment aided by a manual selection of common points between epochs as GCPs was compared to the new automatic pipeline specifically developed for historical aerial photographs.
Our results show that the new method provides more reliable results with an acceptable error. Despite several research questions remain open, our experiment proved that photogrammetry applied to historical aerial photographs is a precious source of data that may help not only produce more accurate landslide inventory maps but also leveraging other relevant analyses within the landslide science community that span from the landscape dynamics to landslide monitoring and quantitative risk assessment.