Brief communication: Landslide motion from cross correlation of UAV-derived morphological attributes

Abstract

Unmanned aerial vehicles (UAVs) can provide observations of high spatio-temporal resolution to enable operational landslide monitoring. In this research, the construction of digital elevation models (DEMs) and orthomosaics from UAV imagery is achieved using structure-from-motion (SfM) photogrammetric procedures. The study examines the additional value that the morphological attribute of “openness”, amongst others, can provide to surface deformation analysis. Image-cross-correlation functions and DEM subtraction techniques are applied to the SfM outputs. Through the proposed integrated analysis, the automated quantification of a landslide's motion over time is demonstrated, with implications for the wider interpretation of landslide kinematics via UAV surveys.

Introduction

Landslides are a form of mass movement, which can often be complex in nature, leading to “slope” failure and the formation of characteristic surface morphological structures. Monitoring of these structures can provide a valuable insight into a landslide's subsurface dynamic failure mechanism and thereby help mitigate hazards (Gunn et al., 2013). Conventionally, in addition to geotechnical and geophysical monitoring of the subsurface, survey markers are often used to quantify surface displacement by monitoring discrete locations through periodic observations. However, such surveying can be hazardous and generally provides limited spatial resolution. The development of low-cost, mini, consumer-grade unmanned aerial vehicles (UAVs) – also known as remotely piloted aircraft systems and drones – equipped with off-the-shelf compact cameras, in combination with structure-from-motion (SfM) and multi-view stereo (MVS) algorithms, has facilitated high-spatio-temporal-resolution topographic surveys using image-based approaches. In particular, the implementation of the SfM–MVS pipeline into user-friendly commercial software packages, such as PhotoScan (PhotoScan, 2016) and Pix4D (Pix4D, 2016), has enabled the generation of high-spatio-temporal-resolution point clouds, digital elevation models (DEMs) and orthomosaics in the earth sciences (Remondino et al., 2014; James et al., 2017).

The differencing of successive co-registered DEMs constitutes a standard approach to estimate ground accumulation and depletion in monitoring applications (Daehne and Corsini, 2013; Travelletti et al., 2014). Moreover, the application of image-cross-correlation functions to optical imagery has long been successfully implemented for the quantification of surface planimetric movement in the context of landslides, glaciers, etc. (Leprince et al., 2007; Ayoub et al., 2009a; Heid and Kääb, 2012). Nevertheless, the application of image-cross-correlation functions to UAV-derived orthomosaics can increase noise due to variations in illumination conditions (Lucieer et al., 2014). Recent studies have demonstrated that the implementation of image-cross-correlation functions with DEM morphological derivatives can automatically determine the movement of surface features that preserve their structural patterns over time (Daehne and Corsini, 2013; Lucieer et al., 2014; Travelletti et al., 2014; Fey et al., 2015). Among these, Lucieer et al. (2014) and Turner et al. (2015) found the UAV-derived morphological attribute of “shaded relief”, combined with image-cross-correlation functions, to provide better surface displacement estimation of a landslide than single bands from the corresponding orthomosaic. To date, however, there has been no comprehensive evaluation of image-cross-correlation functions with various UAV-derived morphological attributes for landslide deformation monitoring.

This paper reports on the analysis of horizontal motion and elevation differences of an active landslide from multi-temporal co-registered UAV-derived outputs, including DEMs, orthomosaics and morphological attributes. Firstly, image-cross-correlation functions are evaluated through comparative analysis with synthetic datasets. Secondly, the surface deformation of a landslide is determined by integrating image-cross-correlation functions with morphological attributes and DEM differencing. The paper illustrates how to exploit a time series of UAV survey derivatives in order to quantify and interpret landslide kinematics.

Study area

The Hollin Hill study site comprises a slow-moving earth-slide, earth-flow landslide with an extent of 290 m E–W, 230 m N–S and a south-facing slope of an average 12${}^{\circ }$, located in the Lias mudrocks of North Yorkshire, UK (54${}^{\circ }$6${}^{\prime }$38.90${}^{\prime \prime }$ N, 0${}^{\circ }$57${}^{\prime }$36.84${}^{\prime \prime }$ W). The site has been monitored since 2009 by the British Geological Survey (BGS) using various methods, including terrestrial and airborne laser scanning, as well as ground-based geotechnical and geophysical investigations. BGS investigations have revealed that the landslide has a complex behaviour with seasonal surface variations and episodic failures mostly triggered by intensive rainfall and increased pore-water pressures within the constituent geological materials (Gunn et al., 2013; Uhlemann et al., 2017).

Data acquisition and processing

Image acquisition was performed using a mini fixed-wing UAV (Quest 300) equipped with a Panasonic Lumix DMC-LX5 compact camera of 5.1 mm nominal focal length and an image array of 3648 pixels $\times$ 2736 pixels. RGB UAV imagery was captured during six field campaigns in December 2014, March 2015, June 2015, September 2015, February 2016 and May 2016. The Quest 300 was flown from a nominal flying height of 90 m at 18 m s${}^{-\mathrm{1}}$, with images acquired approximately every 2 s. During every field campaign, a Global Navigation Satellite System (GNSS) base station was established over stable terrain and surveyed in static GNSS mode. Average absolute accuracies of 0.01 m in planimetry and 0.02 m in elevation were delivered. Circular targets of 0.40 m diameter (equal to 8–10 pixels), with centres easily recognisable in the imagery, were established. Between 11 and 20 targets were surveyed for each of the different campaigns using the rapid static GNSS.

A self-calibrating bundle adjustment, incorporated into the SfM–MVS pipeline, was utilised to process the UAV imagery using PhotoScan software, as described in Peppa et al. (2016). The observed coordinates of five circular targets were utilised as control in each SfM–MVS bundle adjustment, with the remainder used as independent check points. This resulted in the reconstruction of six dense point clouds, one per epoch, georeferenced in the Ordnance Survey Great Britain 1936 coordinate system. From an average 0.03 m ground sample distance, DEMs were generated at each epoch with an average 0.06 m spatial resolution. The 3-D co-registration accuracy, calculated from differences between the surveyed and observed coordinates at independent check points after the SfM–MVS bundle, was estimated as an average root mean square error (RMSE) of 0.03 m. Peppa et al. (2016) described an approach to derive the vertical sensitivity with the use of DEM standard deviations. An approximate $\pm$0.10 m sensitivity level, corresponding to the lowest detectable change, was estimated by applying error propagation (with a 95 % confidence level) to the 3-D RMSE values, calculated at check points. Both approaches resulted in a sensitivity level of the same order of magnitude.

Methodology

Four morphological attributes (shaded relief, slope, openness and curvature) were computed from each epoch's DEM. Shaded relief was created with the aid of the ambient occlusion tool in the SAGA GIS package. This applies homogenous illumination to the DEM, smoothing the shadow effect usually produced by lighting from a single direction (Fey et al., 2015). The remaining three morphological attributes were all generated using the Orientation and Processing of Airborne Laser Scanning data (OPALS) software (Pfeifer et al., 2014). In this paper (a) slope indicates the steepest slope angle of the surface; (b) openness represents the minimum angle of a cone fitted in the DEM, as viewed from above the surface (Yokoyama et al., 2002); and (c) curvature constitutes the average of minimum and maximum curvature, representing concave and convex surface features respectively. All three attributes were computed using a 3 pixels $\times$ 3 pixel window, equivalent to 0.18 m at 0.06 m pixel resolution.

Signal-to-noise ratio (SNR) maps of stable terrain derived from COSI-Corr with (a) openness and (b) shaded relief superimposed over the December 2014 orthomosaic. (c) Box plots of SNRs for stable terrain in Regions A and B, as derived from the implementation of COSI-Corr with shaded relief, slope, openness and curvature applied to synthetic datasets. The median is displayed as a red line, the mean as a red rectangle, the whiskers as black horizontal lines and the outliers as black crosses (Peppa et al., 2017).

[Figure omitted. See PDF]

An experiment was conducted with synthetic epoch pairs to evaluate the performance of the statistical normalised cross-correlation (NCC) function, implemented in the Co-registration of Optically Sensed Images and Correlation (COSI-Corr) software (Leprince et al., 2007; Ayoub et al., 2009b), as applied to these four morphological attributes. To generate the synthetic displacement, known translations of (a) 0.050 m in easting and $-$0.100 m in northing (i.e. 0.112 m total magnitude) were applied to Region A (see Fig. 1a), approximating to the $\pm$0.10 m sensitivity level; and (b) shifts of 0.455 m in easting and $-$0.544 m in northing (0.709 m total magnitude) were applied to Region B in the December 2014 DEM, simulating typical inter-epoch movement of the real landslide. Four pairs of morphological attributes were then derived from both the original December 2014 DEM and the synthetically shifted DEM. Each pair, comprising the pre- and post-event images, was imported into the COSI-Corr function. This computes the maximum absolute value of the correlation coefficient by sliding a rectangular patch from the pre-event image systematically within a window in the post-event image. The computed displacements in easting and northing, determined by the matched correlation peak between the two images, have a spatial resolution equal to a specified step parameter used for the sliding (Ayoub et al., 2009b; Lucieer et al., 2014). After a trial and error procedure, a window size of 64 pixels $\times$ 64 pixels (3.84 m) with a step of 16 pixels $\times$ 16 pixels (0.96 m) and a patch of 20 pixels $\times$ 20 pixels (1.20 m) were chosen for this research. These settings ensured that the maximum imposed shift over Region B could be detected and was therefore chosen in line with a priori knowledge of the Hollin Hill landslide movement rates (Uhlemann et al., 2017). The computed displacements in easting and northing were combined to provide 2-D motion maps across successive epochs.

Apart from the displacements in easting and northing, the COSI-Corr function also calculates a signal-to-noise ratio (SNR), which is indicative of the correlation quality. SNR values closer to unity are indicative of more reliable results. A comparative analysis of the estimated displacements and derived SNRs, obtained with the four morphological attributes, was then performed to determine which of the morphological attributes produced the optimal results. The chosen morphological attribute, together with the COSI-Corr function, was applied to successive epoch pairs of the Hollin Hill landslide to estimate 2-D motion. The COSI-Corr result was cross validated with the surface displacements calculated from 27 sample points manually measured across the orthomosaics. These points were identified on visually identifiable characteristic surface breaks and evenly distributed across the site with displacement magnitudes from cm to m level. The elevation change was derived by subtracting each DEM from the subsequent DEM on a pixel-by-pixel basis.

Having generated a time series of horizontal motions and elevation differences across the site, an additional investigation over subregions with the largest deformations was then performed. The morphological attribute of openness was then chosen due to its unique representation of discernible surface patterns within the landslide body. The NCC function was applied to openness for December 2014 and May 2016 datasets, as implemented in the Correlation Image Analysis (CIAS) package (CIAS, 2012; Heid and Kääb, 2012) using the aforementioned window and patch sizes. Unlike COSI-Corr, CIAS allows individual feature tracking. Thus, characteristic surface structures were manually located over the December 2014 openness image and the derived 2-D coordinates were used as input to the CIAS tool. The planimetric vectors of these locations, between December 2014 and May 2016, were automatically derived with the same tool. Manual cleaning to remove spurious vectors was also necessary, although this process could be automated by the application of various threshold parameters, if necessary. For instance, the sensitivity level could serve as a threshold to remove vectors of lengths lower than $\pm$0.10 m. Based on previous knowledge of the Hollin Hill landslide (Uhlemann et al., 2017), a specific azimuth range could be used as an additional threshold to exclude vectors showing, for example, backward motion due to rotational failures.

Results

Before presenting the horizontal and vertical displacements over the Hollin Hill landslide, the results of the synthetic experiment are firstly described. All four morphological attributes underestimated the imposed displacement of Region A, delivering an average displacement 0.030 $\pm$ 0.027 m in easting and 0.054 $\pm$ 0.030 m in northing. For Region B, the closest result to truth in easting was delivered by openness, with an average value of 0.435 $\pm$ 0.145 m, whereas shaded relief detected the best average displacement in northing of $-$0.528 $\pm$ 0.131 m. Figure 1a and b depict the SNR results, derived from openness and shaded relief respectively, over stable terrain outside Regions A and B. Figure 1c presents the box plots of the comparative SNR analysis. SNR values close to zero (Fig. 1b) indicated decorrelation, which is also illustrated as outliers in the box plot of shaded relief over stable terrain, whereas the other three morphological attributes were less noisy (Fig. 1c). For Regions A and B all morphological attributes with the exception of curvature produced similar box plots. The box plots reveal greater variation in SNR in Region B than in Region A (Fig. 1c), possibly due to the noise caused by the extreme local surface variations around Region B. Overall, slope and openness provided comparable displacements and noise levels. In this study, openness was finally chosen for the estimation of Hollin Hill landslide motions, as it highlights characteristic breaks in slope sliding downwards over time.

The comparison of the COSI-Corr-derived displacements with the manually observed surface movements at 27 sample points (Fig. 2) indicates the sensitivity of the NCC function to different displacement magnitudes. The scatterplot in Fig. 2 shows a general systematic overestimation of the displacement magnitude derived from COSI-Corr. Some scattered points fell within the $\pm$0.10 m 3-D sensitivity level shown in grey, especially for the March–June and June–September 2015 epoch pairs. The 33 and 38 sample points across all epoch pairs with displacement magnitude larger than $\pm$0.10 m were manually observed on orthomosaics and automatically derived with COSI-Corr respectively. Significant movement was observed mostly between the December 2014–March 2015, September 2015–February 2016 and February–May 2016 epoch pairs. Overall, the NCC function delivered results in good agreement with the manual measurements (closer to the straight line) for small displacements but miscalculated the surface movement of the last epoch pair.

The planimetric displacements across the Hollin Hill landslide between December 2014–March 2015, March 2015–February 2016 and February–May 2016 are mapped in Fig. 3a, b and c respectively. Observations from the June 2015 and September 2015 campaigns were excluded from the maps in Fig. 3 due to small displacements and additional noise caused by vegetation change. Blue hatched polygons represent areas with more reliably estimated surface displacements, as the SNR is greater than 0.7. This value is equivalent to the lowest whisker of the openness box plot (Fig. 1c), representing the outlier threshold, as derived from Eq. (1): $$\mathrm{lowest}\mathrm{whisker}=Q\mathrm{1}-\mathrm{1.5}\times (Q\mathrm{3}-Q\mathrm{1}),$$ where $Q$1 and $Q$3 are the 25 and 75 % percentiles of the data respectively.

Scatterplot of estimated surface displacements determined by COSI-Corr with openness plotted against the manual observation per epoch pair (Peppa et al., 2017).

[Figure omitted. See PDF]

Maps of surface displacements and elevation differences of (a and d) December 2014–March 2015, (b and e) March 2015–February 2016 and (c and f) February–May 2016 respectively. Manually derived planimetric vectors at sample points are also superimposed (Peppa et al., 2017).

[Figure omitted. See PDF]

There are a few erroneous displacements, mostly at the edges of the study site, around vegetated areas and outside the blue hatched polygons, as evidenced in Fig. 3a, b and c. The elevation differences between the same epoch pairs are depicted in Fig. 3d, e and f, excluding deformations within the $\pm$0.10 m sensitivity level. Part of the western lobe collapsed, creating a dramatic change of $-$0.70 m in maximum ground loss and $+$0.50 m in maximum ground accumulation within 11 months (Fig. 3e). The surface ruptured at the upper part of the slope, yielding a maximum ground subsidence of approximately $-$1.70 m and a maximum elevation increase of approximately $+$1.05 m, as seen in Fig. 3f. In addition, Fig. 3f depicts the grass growth at the foot of the slope, which in turn caused false surface movement in Fig. 3c. Also, over the regions with extreme deformations (e.g. back scarp in Fig. 3f), decorrelation created voids on the displacement map (Fig. 3c).

To further investigate these significant deformations, the May 2016 openness image was superimposed over the corresponding image from December 2014 and is presented in Fig. 4a and b. Figure 4c illustrates that narrow angles of openness can distinguish surface undulations sliding down-slope. For instance, point 1 moved 1.10 m along the profile AB towards the south. To visualise these structures a threshold of 63${}^{\circ }$ was applied to the openness images (Fig. 4a and b). Different thresholds can visualise different morphological features. This threshold was derived with the aid of visual inspection along profiles at multiple locations over active parts of the landslide. Openness also captured the surface rupture that occurred at the top of the slope between February and May 2016 (Figs. 3f and 4b). The planimetric vectors of distinctive features are plotted in Fig. 4a and 4b, as automatically determined after applying the NCC function implemented in CIAS. Spurious vectors at the edges of the back scarp, which were manually removed, were possibly generated due to rotational failures investigated by BGS (Uhlemann et al., 2017).

Detailed view of December 2014 and May 2016 openness maps over the (a) eastern lobe and (b) back scarp with elevation and openness plotted along (c) profile AB (Peppa et al., 2017).

[Figure omitted. See PDF]

Discussion

The comparative analysis of the NCC function with synthetic data was necessary to tune the function's optimal settings. If small displacements close to the UAV-derived sensitivity level do not fit within the specified window size, they cannot be precisely estimated (e.g. Region A), as was noted by Fey et al. (2015). Small step and window sizes improved the spatial resolution of the surface displacement magnitude map but increased the computational time and noise. This occurred as features with similar/repetitive patterns within the vicinity of the specified window sizes generated false displacements (Travelletti et al., 2014; Fey et al., 2015). Hence, the choice of the function's parameters is usually based on the required spatial resolution, the computational effort and the displacement magnitude (Daehne and Corsini, 2013; Travelletti et al., 2014; Fey et al., 2015).

The analysis with synthetic data also demonstrated that imagery derived from various morphological attributes can generate different displacement estimations and noise levels. Slope, openness and curvature outperformed shaded relief in terms of noise over stable terrain, even though all attributes are insensitive to illumination variations and shadows (Daehne and Corsini, 2013; Lucieer et al., 2014; Fey et al., 2015). A possible error source could be the grass cover, which is well known to affect the results of image cross-correlation (Lucieer et al., 2014; Stumpf et al., 2017).

The production of reliable surface displacements with the image-cross-correlation functions over vegetated terrain constitutes a significant challenge. As vegetation covers surface features, the NCC function generates additional noise. Conversely, grassy surfaces produce images with low texture and without distinctive surface features which can also affect the NCC function's performance (Travelletti et al., 2014), as evidenced in Figs. 1b, 3a, b and c around the eastern lobe. Hence, noisy results attributed to vegetation presence cannot be entirely removed, even with UAV surveys of high temporal resolution. The use of morphological attributes computed with larger spatial distances, thereby producing a higher level of smoothing, can potentially decrease this noise. Finally, to generate optimal NCC results with the least amount of noise possible while accounting for the vegetation variation, winter would constitute the best period to conduct UAV surveys.

Independently of the NCC function's sensitivity to displacement magnitude and vegetation presence, the presented analysis also revealed other limitations, which have already been well reported in previous studies (Daehne and Corsini, 2013; Lucieer et al., 2014; Travelletti et al., 2014; Fey et al., 2015; Stumpf et al., 2017). A priori knowledge of the displacement magnitude is required for tuning the function's settings, therefore somewhat limiting the automated fashion of the workflow. In addition, decorrelation occurs when a surface has significantly changed between two consecutive epochs. Finally, image-cross-correlation functions generate unreliable estimations over regions with rotational failures, creating spurious vectors or voids, whereas performance is much better over translational earth-flow slides. Even though threshold definition can automatically remove spurious vectors, it is not a straightforward process as it relies on a priori knowledge of the landslide. Where such information is unavailable, additional field data may be used. This demonstrates that image-cross-correlation performance is strongly related to the landslide movement type. For mixed types, such as the Hollin Hill landslide (a combination of rotational failures with earth flow, as shown in Uhlemann et al., 2017), the successful application of image cross-correlation is not entirely guaranteed.

Overall, heterogeneous horizontal motions and elevation differences were observed at the Hollin Hill landslide through the combination of multiple co-registered UAV products. The orthomosaics supported the identification of vegetated areas and cross validation of the results. The use of openness, together with the COSI-Corr tool supported the quantification of the movement over the whole site. DEM differencing was also applied to quantify the episodic surface ruptures and interpret the generated voids on displacement maps. The episodic surface ruptures generated vertical ground loss and accumulation, as seen in Fig. 3e and f. The horizontal downward motion of the front part of the eastern lobe was illustrated as positive elevation change. This motion was also identified with the image-cross-correlation analysis (Fig. 3). The CIAS tool applied to openness tracked the evolution of discernible surface patterns over the 8-month duration in a semi-automated fashion. Openness maps of different angle thresholds express surface formations in different ways and as a result can complement the investigation of landslide motion. The exploitation of available image-cross-correlation tools (COSI-Corr and CIAS) with openness decreased the intensive task of manual feature tracking. However, this task is still essential for cross validation, especially in cases where ground truth observations are lacking over the monitoring period.

Conclusions and future work

This paper has presented an investigation of UAV-derived products of DEMs and orthomosaics along with DEM morphological derivatives of openness to automatically quantify the spatio-temporal motion of an active landslide. The research has demonstrated the successful integration of image-cross-correlation functions with morphological attributes and the importance of the comparative analysis with synthetic data. The analysis has illustrated that openness implemented with image-cross-correlation functions can be used in conjunction with DEM differencing to support the comprehensive interpretation of landslide behaviour, providing a holistic overview of horizontal and vertical deformation patterns. Major limitations include the reliance on a priori knowledge of the landslide type and displacement magnitude to tune the image-cross-correlation function parameters, use of field data for cross validation, manual surface feature identification and manual cleaning, or threshold definition to remove erroneous displacement vectors. These limitations affect the performance of the resulting horizontal motions and elevation changes. Future work will assess the performance of image-cross-correlation functions with shaded relief, slope, curvature and other possible DEM derivatives computed with various pixel radial distances and implemented with real-world data. It will also apply other techniques to automatically filter spurious results. Ultimately, future research will investigate the correlation of the horizontal motions and elevation differences with rainfall observations to enhance the understanding of the landslide mechanisms.

Data supporting this publication are openly available under an “Open Data Commons Open Database License”. Additional metadata are available at 10.17634/154300-58 (Peppa et al., 2017). Please contact Newcastle Research Data Service at [email protected] for access instructions.

The authors declare that they have no conflict of interest.

This article is part of the special issue “The use of remotely piloted aircraft systems (RPAS) in monitoring applications and management of natural hazards”. It is a result of the EGU General Assembly 2016, Vienna, Austria, 17–22 April 2016.

Acknowledgements

This research was jointly funded by a Natural Environment Research Council (NERC) BGS BUFI award (S241) and an Engineering and Physical Sciences Research Council (EPSRC) DTA award (EP/L504828/1) at Newcastle University, UK. The BGS contribution to this paper is published with the permission of the Executive Director of the British Geological Survey (NERC). The authors wish to acknowledge the Hollin Hill site landowners for allowing access to their property. Finally, many thanks to Martin Robertson, Elias Berra, Magdalena Smigaj, Polpreecha Chidburee and Ben Grayson, as well as all of Newcastle University, for fieldwork assistance. Edited by: Daniele Giordan Reviewed by: two anonymous referees