Estimating Landslide Surface Displacement by Combining Low-Cost UAV Setup, Topographic Visualization and Computer Vision Techniques

Abstract

Many techniques are available for estimating landslide surface displacements, whether from the ground, air- or spaceborne. In recent years, Unmanned Areal Vehicles have also been applied in the domain of landslide hazards, and have been able to provide high resolution and precise datasets for better understanding and predicting landslide movements and mitigating their impacts. In this study, we propose an approach for monitoring and detecting landslide surface movements using a low-cost lightweight consumer-grade UAV setup and a Red Relief Image Map (a topographic visualization technique) to normalize the input datasets and mitigate unfavourable illumination conditions that may affect the further implementation of Lucas–Kanade optical flow for the final displacement estimation. The effectiveness of the proposed approach in this study was demonstrated by applying it to the Ruinon landslide, Northern Italy, using the products of surveys carried out in the period 2019–2021. Our results show that the combination of different techniques can accurately and effectively estimate landslide movements over time and at different magnitudes, from a few centimetres to more than several tens of meters. The method applied is shown to be very computationally efficient while yielding precise outputs. At the same time, the use of only free and open-source software allows its straightforward adaptation and modification for other case studies. The approach can potentially be used for monitoring and studying landslide behaviour in areas where no permanent monitoring solutions are present.

1. Introduction

Landslides are considered one of the most devastating geohazards, causing significant damage and loss of life. Accurate and timely monitoring of landslide motion is critical for studying and predicting their possible behavior, as well as to further define suitable mitigating measures to lessen their impact. The growing interest in using remote sensing technologies in recent years is evident in many fields, as technologies such as air- and spaceborne sensors provide invaluable and insightful datasets with constantly increasing capabilities from the spatial, temporal and spectral perspectives [1]. In the context of satellite-based sensing and the domain of landslide hazards, many new approaches and methodologies have been developed and implemented around high-to-medium resolution optical multispectral missions [2,3,4] and weather-invariant radar sensors [5,6], where the main applications are related to mapping historical events, monitoring active and potential slopes [7,8,9], as well as susceptibility mapping and hazard assessment.

Similarly, there is increasing interest in the use of Unmanned Aerial Vehicles (UAV) for landslide studies [10], as UAVs offer several advantages over traditional monitoring methods, including the ability to collect high-resolution imagery, to access difficult-to-reach areas, and to collect data frequently and at a lower cost [11,12]. Therefore, it is natural for scholars to adopt and further exploit this technology also for landslide time-series monitoring [13,14], mapping and characterization [15,16,17] as well as combining their outputs with machine learning [18] and computer vision algorithms [19,20].

For the estimation of landslide surface changes through UAV flyovers, scholars have traditionally used straightforward techniques, such as measuring the difference between Digital Surface Model (DSM) elevation products from consecutive surveys (see [12,16,21,22]) or the direct difference between point clouds (e.g., [11]), which usually results in the quantification of volumetric changes. Such approaches are usually powerful when the landslide under investigation exhibits significant and abrupt changes, however, for some smaller changes and displacements, they can be less efficient, as most of the time alignment error between two products is introduced [23,24], which makes their estimation more difficult. At the same time, scholars have mostly been utilizing high-precision positioning tools for the computation of apparent velocities (e.g., [21,24], which is not always feasible due to the specifics of a terrain. A notable issue in the literature is often related to uneven weather conditions, illumination conditions in particular [20,21,24], which mostly affects the final products from the surveys and then the analyses. Therefore, in the current work, we propose an approach that aims to tackle most of the issues affecting precise surface displacement estimation.

One of the challenges in using UAV-derived imagery for landslide monitoring lies in accurately measuring the apparent surface displacement of the landslide body over time. One solution is using computer vision techniques such as optical flow algorithms [20]. Optical flow algorithms can estimate the motion of objects or features in a sequence of images or videos by tracking the movement of pixels between images/frames and computing the apparent motion of those pixels [25]. The Lucas–Kanade [26] optical flow algorithm is a widely used method that is based on the assumption of small or constant pixel intensity changes between frames. In fact, the latter can cause reconstruction issues and further erroneous displacement estimations when applied in a time-series manner for landslide monitoring, mainly due to unfavorable illumination conditions [19,20] that may occur in one or more epochs. Such conditions are related to different sun directions, shadows and cloud cover, all factors that directly influence image brightness and contrast, i.e., that can further alter the algorithm if directly applied to orthorectified imagery.

An approach for resolving similar issues and ensuring a consistent image intensity is to input normalized products from the UAV surveys. An alternative to the use of orthomaps could be the implementation of a Red Relief Image Maps (RRIM) [27] on three-dimensional data, with photogrammetric products such as DSM derived from UAV surveys. RRIM is a topographic visualization technique, and an alternative to hillshading, which can more successfully highlight morphologic features. The visualization relies on topographic openness [28] and a red-colored slope gradient, and has already been adopted in various studies (e.g., [29,30]).

In this study, we propose a method for monitoring and estimating landslide surface displacements using the products of time-series surveys with a low-cost, consumer-grade UAV and the derived photogrammetric products, such as Digital Surface Models, to which the Red Relief Image Map visualization technique combined with the Lucas–Kanade optical flow algorithm was applied in order to estimate the displacement of pixels between images. The proposed approach allows the measurement of landslide surface motion in 2D space. The effectiveness of this method is demonstrated by applying it to the Ruinon landslide in Northern Italy using datasets from surveys carried out in the period 2019–2021. Processing is done using entirely free and open-source software, which broadens the possibility of adoption and reproducibility. Our method can potentially improve our understanding of landslide behavior and contribute to the development of effective risk mitigation actions in locations where they are lacking.

2. Materials and Methods

2.1. Case Study—Ruinon Landslide

The Ruinon landslide (Figure 1) is located in the Upper Valtellina, Northern Italy, and is one of the most active rockslides in the Italian Alps. Throughout the years it has been under continuous investigation and monitoring from scholars, experts and local authorities [31,32,33]. The active Ruinon landslide actually lies on top of a deep-seated gravitational slope deformation (DSGSD), recognized in 1998 [34], while local failures were documented from 1954, with several airborne surveys of the area having been organized as early as 1954 and 1981 by the Lombardy regional authorities [31].

While the DSGSD reaches up to 3000 m a.s.l almost on the whole hillslope, the Ruinon landslide is divided into two scarps, the lower ranging from around 1400 to 1900 m a.s.l, and the upper at around 2100 m a.s.l. In the last few years, the more active scarp has been the lower one, with planar and toppling failures which has resulted in displaced debris of various sizes (from gravel size to about 10 m). The main body of the lower scarp is bounded on the left-flank side by the Confinale river and on the right by an almost vertical rockface, and currently, the toe reaches a newly constructed rock barrier for protecting Provincial Road SP29 (Figure 2).

In fact, SP29 is the main connection to Santa Caterina di Valfurva, a small village home to famous ski resorts. Throughout the years, reactivations at Ruinon were a direct threat to people’s lives, but also to the socioeconomic state of the region. The occurrence of failures has a recurrent character, with activations predominantly occurring during the summer–autumn months, with relatively calmer periods over the winter period. Such recurrent activities have already been monitored and reported almost every year since 2006, when a permanent Ground-Based Interferometric Synthetic Aperture Radar (GBInSAR) system was installed [32,34,35,36]. Aside from the GBInSAR, the landslide has been monitored since 1997 through a manual and automatic geotechnical network [33]. In addition to the data from permanent ground stations during the last two decades, the Ruinon landslide was also monitored through spaceborne optical and radar missions [36,37], which recorded the most disastrous reactivations during the summers of 2016, 2018 and 2019. During the year 2021, no significant critical activities on the landslide body were reported [38].

2.2. Workflow

Excluding the actual UAV surveys, the implemented workflow in this work is generally divided into four main steps (Figure 3): processing the acquired datasets; post-processing their products; producing Red Relief Image Maps and applying a finer co-registration, and final displacement computation. To carry out a time-series analysis using the UAV-derived products, we used two separate datasets: six surveys carried out during the period 2019–2020 by the Regional Agency for the Protection of the Environment (ARPA) of the Lombardy Region (Section 2.3.1), and surveys carried out by us in the period 2021–2022 (Section 2.3.2). Both datasets will be further presented in the relevant sections. As there were generic differences between the datasets, the 2019–2020 set being composed only of point clouds and orthophotos, a post-processing step was needed to put them in a common framework, for example, with spatial subsampling of the point clouds, co-registration, etc. This step was followed by producing the RRIM for each epoch and a step-wise finer co-registration of the maps carried out into two sub-steps—a local co-registration of the stable regions and a global one for the unstable parts. The last step was the computation of the dense optical flow and vector magnitude. To ensure higher reproducibility, all the steps were carried out using only Free and Open Software Solutions (FOSS) and Python packages. Namely, for producing the point clouds, DSM and orthoimages were processed with OpenDroneMap [39], CloudCompare [40] was used for their co-registration, Relief Visualization Toolbox for Python [41] for producing RRIMs, and the Lucas–Kanade implementation in the scikit-image Python package [42] for displacement computation.

2.3. UAV Surveys and Their Outcomes

2.3.1. Period 2019–2020

In addition to the permanent monitoring station, during the period of 2019–2020, and due to continuous Ruinon activations, the local agency ARPA carried out UAV surveys over the area. Upon request, the resulting point clouds and orthophotos were shared with us. In

Table 1. UAV surveys carried out from ARPA Lombardy in the period 2019–2020.

	26 July 2019	4 September 2019	27 September 2019	25 October 2019	10 September 2020	19 October 2020
Total points	12,807,656	21,408,313	21,835,284	81,466,454	53,722,814	116,372,943
GCPs	X	X	X	X	X	X
Area covered [km${}^2$]	0.63	0.75	0.67	0.62	0.71	0.69
Point density [pts/m${}^2$]	20	28	33	130	75	168

, the dates of acquisition, as well as the reconstructed points are reported, while the point clouds are shown in Figure 4. As we were provided only with the output products, there is no information related to the survey or flight parameters, nor captured images per survey. Moreover, during the data acquisition or surface reconstruction steps for those specific surveys, Ground Control Points (GCPs) were not used, and therefore they could not be evaluated in terms of accuracy. In addition, by visual inspection, significant local and global misalignments among the different epochs were noted. A different amount of reconstructed points was also noticeable, where the first products were sparser compared to the latest. However, the aforementioned issues were overcome with post-processing (Section 2.4), and the point clouds were sufficient as initial input to be incorporated into the analysis.

2.3.2. Period 2021–2022

As part of our surveys, four flyovers were carried out during the period 2021–2022, over the same Ruinon landslide area. The main idea was the continuation of the surveys in terms of time and adopting novel approaches for monitoring landslide phenomena. Moreover, carrying out the surveys first-hand provides full control over all activities from the planning stage to the final outputs. As the recent reactivations were during summer periods, the surveys were planned to be before and after potentially unstable periods, i.e., spring and autumn periods were most suitable. The current work follows and extends the scope of a previous one [43] in which the main principles of UAV survey planning and processing were explained in detail, to be easily adoptable and reusable in other research domains (e.g., archaeology and glaciology [44,45]). A brief summary of the setup and steps will be presented below.

All of the surveys were made using a consumer-grade quadrotor DJI Mini 2 (249 gr.) equipped with a built-in position system able to connect to GPS, GLONNAS and GALILEO constellations. The UAV has a 12 MP CMOS sensor with an f/2.8 lens onboard. Some of the advantages of using this aircraft can be related to its size and weight. The low body weight puts the model into a low-risk category for which, according to the European regulations [46], the operator needs to do a short online course and exam. The small aircraft has wind resistance at level 5, an advantage considering the weather conditions which may occur during a landslide survey in mountainous areas. On the other hand, there are some drawbacks of the UAV, again related to body size and weight. For instance, the quadrotor lacks an obstacle avoidance system, from through which other optional features such as designing and following a survey plan are also missing. To overcome the lack of mission planning, Flight Planner [47] was used, which is a QGIS [48] plugin useful for planning airborne photogrammetric surveys. The plugin needs three main inputs—vector shape of the area of interest, digital elevation model (DEM), and target ground sampling distance (GSD), with the outputs being projection centers, flight paths and waypoints.

After data acquisition, the collected images were used to produce point clouds, digital surface models and orthophotos, using the OpenDroneMap tool [39], which is a FOSS for processing aerial data based on state-of-the-art Structure from Motion (SfM) and Multi-View Stereo (MVS) approaches. A summary of the point cloud characteristics is reported in

Table 2. UAV surveys carried out in the period 2021–2022.

	6 July 2021	29 October 2021	2 May 2022	5 July 2022
Total points	297,421,177	407,189,985	438,663,249	136,304,323
GCPs		X
GCP RMS error [m]	0.15	X	0.16	0.13
Area covered [km${}^2$]	0.45	0.52	0.45	0.45
Point density [pts/m${}^2$]	653	783	975	303

and visualized in Figure 5.

During the surveys, at least five GCPs were placed and measured using Leica AX1200 in real-time kinematic positioning (RTK) mode connected to Bormio permanent reference station (https://www.spingnss.it/, accessed on 28 October 2022). Three independent measurements were made for each marker. The GCPs were distributed in an area considered stable and far from landslide movements (see Figure 4 in [43]). It should be noted that the GCPs used are not distributed in optimal positions around the landslide body, but, considering its size, terrain and the fact that it is still an active area, field activities at higher altitudes towards the crown were not allowed.

The final products achieved the targeted GSD and all of the outputs have a spatial resolution of at least 10 cm/pix or better. All orthophotos are available and free for download from the platform for sharing aerial images OpenAerialMap (accessed on 30 October 2022) https://map.openaerialmap.org/#/10.358734130859375,46.45630839492814,10/square/120221312?_k=hzzlob. It should be noted that all surveys followed the same case-specific flight plan and processing scheme. Namely, maximum flight altitude was always 110 m above ground; when flying in a longitudinal direction 11 tracks with vertical speed maintained at 3 m/s were followed, horizontal speed at 4.5 m/s; camera inclination was at −57${}^{\circ }$, while in transversal directions; 14 paths were flown with the vertical speed at zero and horizontal at 6 m/s; the camera was set at −90${}^{\circ }$. However, as mentioned before, the UAV lacks autonomous flight, therefore all were carried out manually by an operator, maintaining the above-mentioned parameters, which can explain the small differences in the reconstructed points. Among all of them, the one with the sparsest point cloud is that of the 5 July 2022 survey. Even after a deep investigation, the reason for this lower point density could not be determined, as the detected and reconstructed features were very similar in location and count compared to the rest of the products. Nevertheless, the output was still suitable for the following processing steps.

At a first glance, it can be noted that all of the point clouds cover the whole landslide body (including 2019–2020) and the main changes on its surface are visible. However, it can also be noted that point clouds are produced under different light conditions and vegetation states. This problem will create challenges and errors when applying surface displacement analysis, an issue which we are trying to solve with the current work.

2.4. Postprocessing

The survey outputs are from different sources and used different processing techniques which, naturally, leads them to exhibit diverse characteristics and the need to be placed in a common framework, mainly to produce coherent DSMs for the 2019–2021 period. These steps were performed mostly using CloudCompare [40] with only final filtering applied in QGIS, including:

• Spatial sampling of the point clouds to achieve relatively similar density (set to a minimum distance of 10 cm).
• Co-registration of the point clouds: In most cases, Iterative Closest Point (ICP) co-registration was sufficient; when it was not, it was also performed manually using common points between the clouds. The 2019–2020 dataset had unknown accuracy, and it was decided that those epochs should be co-registered according to 2021–2022, specifically, the 6 July 2021 survey was set as a reference point. However, the changes during the 2019–2020 period were severe ,and therefore the moving primary scheme was applied in reverse order, i.e., 19 October 2020 was co-registered to 6 July 2021, 25 October 2019 to 19 October 2020, etc. On the other hand, the epochs in the 2021–2022 period were co-registered using a fixed primary scheme and all three products were transformed according to 6 July 2021. The achieved co-registration error varied among the epochs in the range of 0.20 m < RMS < 0.30 m, and most of the highest values were in the 2019–2021 dataset.
• The co-registered point clouds were rasterized into DSM with a spatial resolution of 10 cm/pix.
• Filtering the resulting DSMs: As the point density was reduced to a certain level, the interpolation of the point cloud into raster led, in some cases, to noisy and irregular terrain, which may further alter the analysis. A Gaussian filter (kernel type—square and radius of 3 pixels) was applied to obtain a more regular and smooth surface product without losing details (e.g., Figure 6).

2.5. Red Relief Image Map

Red Relief Image Map (RRIM) is a visualization technique introduced by Chiba et al. [27] as a method for highlighting topographic features in greater detail. The implementation of RRIM is an extension of the parameter openness [28], introduced as an image-processing technique which emphasizes the concavities and convexities of the surface, combined with the slope gradient. The RRIM technique is proposed as an alternative visualization approach to topographic information, based on the better perception of the color red by the human eye. In particular, the method relies on differential openness:

$$\frac{O_p-O_n}{2},$$

(1)

where O${}_p$ is the positive openness, i.e., the convexity of a surface, and O${}_n$ is the negative openness, i.e., the concave features of a surface. Once computed, the differential openness is blended with a slope gradient layer which is pseudo-colored in red. An example is shown in Figure 7, depicting the area of Upper Valtellina as separate products needed to produce RRIM using SRTM30 data, including also a satellite view and a hillshade visualization.

Black color RRIM depicts concave features, and white convex and grey colours are flat surfaces. In the red gamut, steep slopes are highlighted with a saturated red colour, with convexity and concavity in the lower range. An important parameter for the computation of both positive and negative openness is the radial limit (i.e., the search radius of the moving window), which is directly linked to the level of detail to be highlighted.

Looking at the figure, the advantages of using the RRIM visualization are evident. On the one hand, it enhances the perception of detail and depth, compared to the optical image or the hillshade visualization, on the other hand, it is not altered by the sunlight direction, shadows, and vegetation state. Those particular features are of great benefit when applying algorithms for computer vision, especially for a landslide monitoring task. The RRIM clearly depicts topographic features of various sizes and shapes, even on flat surfaces, and even if the level of the details is obviously highly dependent on the resolution and quality of the input dataset.

The tool used for producing RRIM was the Relief Visualization Toolbox for Python [41,49]. It is a set of instruments aiming to generate various visualizations of elevation models with a focus on small-scale features. Some of the implemented methods include but are not limited to, multiple-direction hillshading, a simple local relief model, sky-view factor, etc. RRIM is not directly available off the shelf but the Relief Visualization Toolbox can be easily adjusted for that purpose.

2.6. Dense Optical Flow

In this work, a traditional dense optical flow technique for determining the apparent motion of features in the image space [50] between two consecutive instances was used. In the last couple of decades these techniques have been widely used in computer vision problems [51,52], as well in the domain of Earth Observation [20,53,54]. For the current implementation, we took advantage of the Le Besnreais and Champagnat [55] implementation of the Lucas–Kanade [26] dense optical flow using iterative registration, which works with a ’local’ search window approach for brightness pattern shifts. It is very cost-efficient due to the coarse-to-fine inverse pyramid strategy for the window search. Since it is a method based on image intensity, a few drawbacks can occur during its implementation: an abrupt change in illumination between two images can create shadows which influence pixel brightness; vegetation state can also change between two surveys. In addition to intensity-related issues, it is important to have the products properly co-registered because local or global shifts would be also considered as motion. Such problems will lead to de-correlated and wrongly oriented vector fields in the final product.

In this work, the Iterative Lucas–Kanade approach was used through its implementation in the scikit-image Python package [42]. As it was also used for image registration, we implemented the dense optical flow computation in two steps: (1) finer image co-registration (Section 2.6.1) and (2) the actual landslide displacement computation (Section 2.6.2).

2.6.1. Finer Co-Registration

Since the phenomenon we are trying to model is dynamic, a direct image-to-image registration cannot be applied, as some of the areas can be wrongly wrapped to their initial position, even though it is an actual displacement. To overcome the issue and to obtain a finer co-registration of two consecutive products, a dense optical flow approach was applied at two instances for the secondary image (Figure 8). Firstly, a local co-registration was performed on stable parts outside the active deformation, followed by more rigid shift correction which does not interfere with the actual debris displacements, and finally, both co-registered segments are combined in one. Note that high and consistent shift values may be an indication of decorrelation due to insufficient overlap between two images.

2.6.2. Displacement Computation

The final step of processing is the actual displacement computation. The merged co-registered secondary image from the previous step was used directly to compute the displacements compared to a previous epoch. The output of the processing gives the displacement vectors in the u,v directions, from which the magnitude of the displacement can be easily computed:

$$\left|\right|u,v\left|\right|=\sqrt{u^2+v^2},$$

(2)

and the orientation of the vector

$$\theta =ta{n}^{-1}\left(\frac{v}{u}\right).$$

(3)

Since the UAV was manually operated during the surveys, some ambiguities in the reconstruction part can be seen. This is evident in the reconstruction of vegetation (bushes and trees). Due to this drawback, after the displacement magnitude and orientation computation, many decorrelations can be depicted. They are most notable when inspecting the vector orientations, depicting movements parallel to the general slope or even upslope. To mitigate this issue and keep only the relevant displacement, a filter according to the displacement orientation was applied. The general slope orientation of the landslide was estimated to be at ${225}^{\circ }$ (North is at $0^{\circ }$). The applied filter took into account this dip direction $\pm {45}^{\circ }$, and therefore the kept displacements were only the ones orientated between ${180}^{\circ }$ and ${270}^{\circ }$. The $\pm {45}^{\circ }$ range deviation from the actual dip direction was chosen due to the current work scale and the fact that the movement trajectory of some small debris does not need to exactly follow direct downslope shifts, since they can easily deviate from it.

3. Results

3.1. RRIM Implementation

Following the previously discussed processing workflow scheme, ten red relief image maps for the period 2019–2022 were obtained as outputs from the UAV surveys of the Ruinon landslide with a non-regular time-series. The advantages of using RRIMs compared to orthophotos for applying computer vision approaches are better highlighted at lower scale representation. To show this fact, details of the RRIM implementation in comparison to orthophotos and hillshade visualization are depicted in Figure 9, while the full-scale RRIMs for the period 2021–2022 are added as an Appendix A.1, Appendix A.2, Appendix A.3 and Appendix A.4 (available [56]).

For a full comparison with the orthophotos, the latter can be accessed as mentioned in Section 2.3.2 from the https://openaerialmap.org/ (OpenAerialMap.org, accessed on 30 October 2022) platform. The two rows of image subsets (Figure 9) depict the same subregions of the landslide body during two different periods—summer (6 July 2021) and spring (5 May 2022). Clear differences are notable in the orthophotos (Figure 9a,b), mainly associated with the lighting conditions and vegetation state. In the first example strong shadows are present which are increase the visual depth perception of the image, but so are defining areas with high contrast and a lack of detail in the underexposed areas. On the other hand, the second instance was a result of a survey carried out during an overcast day lacking direct sunlight, which created a consistent illumination of the terrain without harsh shadows but also lacks clear contrast among the objects. Examples of the different conditions are the two highlighted boulders A and B with their related displaced locations (A’,B’). Moreover, there is a strong difference in the vegetation state growing on the landslide body with seasonal variation—in July the leafage is at its maximum, while in May there are still signs of wilted grass and shrubs. From the perspective of determining displacements through intensity-based computer vision approaches, these factors are critical as they may greatly worsen the final results. As a solution, implementing RRIM can significantly aid in visualizing a normalized representation of the terrain and provide a better map on which to apply intensity-based algorithms.

3.2. Displacement Computation

In the following section, the results of the application of dense optical flow for displacement computation between two consecutive survey epochs, i.e., using the so-called moving primary strategy, are presented. The approach was applied to the RRIMs which were also finely co-registered by two independent iterations—one reducing local shifts over the stable areas outside the unstable slope, followed by another more rigid (global) co-registration of the actual landslide body. As highlighted by the results (Figure 10), the landslide was mostly active in the period 26 July 2019 to 6 July 2021, while the exact initial and final dates cannot be determined due to a lack of observations.

It is important to mention the fact that the reported displacements are already filtered according to the moving direction and refer only to two-dimensional XY displacements. Moreover, the visualized arrows depict only the movement direction, and their length is not associated with the magnitude of the displacement. Upon inspection of the maps in Figure 10, the most distinguishable movement can be seen to be amid the central landslide body and, especially during the initial four periods, the most active and severe parts are at the upper and central areas.

Interestingly, there is the appearance of a zero-displacement patch of the magnitude field starting from the first period (Figure 10a) and “moving” downhill at each following epoch (Figure 10b–d). Upon visual interpretation of the RGB orthomaps, it was determined that those areas, at the initial stages, were forested and then shattered by displaced debris at different stages of the mass movement. In contrast, the actual reason for them to appear “stable” with zero displacements was actually the decorrelation between two consecutive RRIMs in those specific regions, i.e., the abrupt landcover change from forest to debris made it impossible for the algorithm to find corresponding pixels and works only in 2D space without accounting for changes along the elevation axis. Similar decorrelation zones can be also noted in the sequential periods. After 14 September 2019, there is a significant, yet gradual, intensity decrease of the estimated surface displacements, whereas after 6 July 2021, some sporadic movements were mapped which upon visual validation are either from single displaced boulders or inconsistent reconstructions of shrubs and trees. Such false positive displacements, due to the presence of vegetation, are not rare, however, in these cases, the “detected” displacements are isolated and do not exhibit a spatio-temporal coherence among the epochs.

On the other hand, an example of a localized movement with a spatial consistency over time can be easily recognized on the left landslide flank (Figure 10e–g), exhibiting a triangular shape active during the main body reactivations. However, during the period between 6 July 2021 and 29 October 2021, this was the only main debris displacement mapped, and it was not further detected during the following observations.

The overall displacement magnitude can be inspected by epoch using box plots (Figure 11), highlighting the general decrease in the trend of surface displacements and re-occurring major reactivation during the fourth period (25 October 2019–10 September 2020), best seen in the median value (vertical line in the boxes). From the available observations, only a couple of epochs exhibit magnitudes higher than 10 m, namely the first and fourth, which does not mean that the rest are of negligible intensity. The summary box plots of magnitude (Figure 11) is further evidence of the false-positive displacements marked as outliers (denoted as dots) from the overall magnitude field of the epoch.

4. Discussion

The current work proposed an approach for obtaining landslide displacements using time-series-derived UAV datasets. The primary research problem that needed to be solved was related to inherited discrepancies in the datasets themselves. Some common ground characteristics were present among them, such as that the surveys were carried out over the same area of interest and the same type of sensor was used, namely, an optical RGB sensor. This is the extent of the commonalities and different features such as camera sensor model, airborne platform, the experience of the UAV operator, post-processing techniques, illumination conditions and vegetation state prevailed. Thus, present datasets over time can provoke difficulties in their usage, especially when the approach to be adopted, such as dense optical flow, is sensitive to image intensity [20,25]. To overcome such issues, the adoption of RRIM [27], a visualization technique directly applied to DSM products of the UAV surveys, was proposed. With the use of normalized images, the use of the resulting maps was instrumental in overcoming an issue related to different and unfavorable illumination conditions that could alter the application of dense optical flow directly on orthoimages.

The combination of different approaches resulted in a successful method for detecting movements at various scales, from highly intensive periods where the central part of the landslide collapsed, as a whole (e.g., Figure 10a–f and Figure 11), to some smaller localized instabilities such as the aforementioned zone on the right landslide flank (best visible at Figure 10g).

In fact, the latter would have been difficult to distinguish through more conventional approaches, such as point cloud comparison using a Multiscale Model to Model Cloud Comparison (M3C2) [57] approach or cross-sections [58]. As depicted in Figure 12, our approach accurately estimates the spatial extent and magnitude field for the localized movement which is not really clearly highlighted within the same intensity range in the M3C2 comparison of the point clouds (Figure 12b) for the same period. On the other hand, using DSMs it is also possible to manually inspect the elevation profiles from both surveys at a smaller scale (Figure 12c) and to compare the measurements for a particular location. The represented elevation profile actually is a fraction of the central landslide body (>500 m), and using cross-sectional elevation profiles could make the detection of such localized movement almost impossible.

By changing the spatial and intensity scale, the approach succeeded in precisely estimating velocity vectors through the epochs of the most severe reactivations of the central part of the landslide body, which obviously cannot be considered for estimating volumetric changes and more abrupt land cover changes (e.g., loss of forested areas). Some of our results derived from the 2019 surveys (Figure 10a–c) are consistent with the results obtained from another approach we implemented for estimating Ruinon reactivations through Sentinel-2 imagery [37]. Furthermore, comparable outputs for the Ruinon landslide are reported by Carla et al. [36], with the difference that their work utilized a GBInSAR monitoring technique. In the aforementioned work, some of the observations for 2019 are reported per month, in the current work, in contrast, we do not have surveys that cover exactly a full month and could be used for comparison. However, the intensity of our initial 3 epochs from 2019 is similar to their monthly cumulative displacements from August to October 2019. Moreover, the estimated extreme displacements from the GBInSAR coincide with our patches with lost coherence (already discussed in Section 3.2). Similarly, if we produce a cumulative map for the period 4 September 2019 to 25 October 2019 (Figure 13), our patch of no coherence is at a location that coincides with a zone with displacements higher than 50 m for the same period (in [36]). Moreover, the areas directly around the patch of no coherence exhibit comparable magnitude fields.

It needs to be highlighted that the optical flow computation, i.e., in this case the displacement computation, is sensitive to the search radius parameter. By applying a low value, the algorithm cannot find the corresponding location outside a certain range. Therefore, such cases as the ones mentioned before can occur, and the relevant result can exhibit zero displacement, even if this is not true. The issue can be overcome by simply increasing the radius, but this has to be done without overestimating since this approach can lead to the opposite effect—falsely “detecting” certain features or exaggerating their magnitudes. Therefore, this parameter’s value is highly dependent on the case study under investigation.

A possible limitation of the proposed approach is related to the already mentioned fact that it can only estimate displacements in the 2D-space plane, disregarding movements along the vertical axis, and further implementation of similar approaches for 3D motion estimation either on point clouds or voxels is the natural continuation of the current work. Similarly, the approach estimates displacements in the 2D-space, but an essential part of it is a DSM to derive an RRIM. Without elevation data, the topographic visualization could not be applied to obtain a normalized image map. Still, the displacement estimation could be applied directly on an orthomap, yet the results should be considered with great attention as many errors could occur if orthomaps have inherently different image intensities (i.e., great variations in brightness and contrast). Further improvements could also be related to an application of spatial consistency filters which can disregard inherited errors from the reconstruction process, similar to the previously discussed issues with vegetation. It is worth pointing out that this is not a very straightforward task as there is always the possibility that such a filter can ignore movements that are not spatially or temporally coherent for the particular time series.

The presented methodology successfully delivered more than satisfactory results, as it was set to solve challenging problems related to, on the one hand, the utilization of input field surveys from different sources, and on the other, the analysis of a significant landslide which exhibits reactivations that are non-coherent in terms of space, time and intensity. Due to the higher spatial resolution outputs from this approach, it can be further utilized for validation purposes regarding other methodologies that use lower-resolution datasets, such as satellite-derived ones. Overall, the setup can be considered low-cost and computationally inexpensive, allowing deployment for regions and mass movements where no permanent monitoring stations and tools are available, with further tuning of the time intervals according to specific needs.

5. Conclusions

The Ruinon landslide is a landslide that has been documented for more than two decades and which has occasionally reactivated, as in the last years. As such, it has served as an interesting case study for scholars to perform various investigations, whether they are field-based or remote (close-range or air-/space-borne), which further provides invaluable input and output for further developing and testing new approaches. The current work presents a method combining a low-cost UAV surveying setup with a topographic visualization technique, RRIM, and a computer-vision-based algorithm for movement estimation. The utilized consumer-grade quadrotor proved to be convenient from a regulation point of view and in delivering high-resolution images, which were further used for scene reconstruction in the form of point clouds, DSMs and orthorectified maps. Along with our own surveys, the time series dataset was further extended with other datasets provided to us. Combining data products from different sources entailed a normalization step to be included in processing. For this reason, the Red Relief Image Mapping technique was applied to the elevation models, which had a two-fold effect. On the one hand, the datasets from different sources were put into a common framework, and on the other hand, it solved the issue of different unfavorable illumination conditions during the surveys. The latter is in fact a significant issue when an intensity-sensitive algorithm for motion estimation, in this case—the Lucas–Kanade optical flow, is applied. With the results from the method we have presented in this paper, we successfully derived surface mass movements within different spatial and intensity scales, from a few centimeters to more than several tens of meters. Moreover, the resulting accurate outputs are produced in a computationally efficient manner through the use of free and open-source software which further allows for the method’s straightforward adaptation and modification in other case studies. Several follow-up research directions could potentially further contribute to the approach, namely, a transition from 2D space to 3D and an application of a filter for spatial coherency, helping to disregard false estimations. Nevertheless, even in this simplified version, our approach can at least contribute to the validation of remotely sensed lower-resolution products and to the monitoring of landslides that are not permanently monitored in situ.