Time Series Analysis of Landslide Dynamics Using an Unmanned Aerial Vehicle (UAV)

Abstract

In this study, we used an Unmanned Aerial Vehicle (UAV) to collect a time series of high-resolution images over four years at seven epochs to assess landslide dynamics. Structure from Motion (SfM) was applied to create Digital Surface Models (DSMs) of the landslide surface with an accuracy of 4–5 cm in the horizontal and 3–4 cm in the vertical direction. The accuracy of the co-registration of subsequent DSMs was checked and corrected based on comparing non-active areas of the landslide, which minimized alignment errors to a mean of 0.07 m. Variables such as landslide area and the leading edge slope were measured and temporal patterns were discovered. Volumetric changes of particular areas of the landslide were measured over the time series. Surface movement of the landslide was tracked and quantified with the COSI-Corr image correlation algorithm but without ground validation. Historical aerial photographs were used to create a baseline DSM, and the total displacement of the landslide was found to be approximately 6630 m³. This study has demonstrated a robust and repeatable algorithm that allows a landslide’s dynamics to be mapped and monitored with a UAV over a relatively long time series.

1. Introduction

Landslides are a worldwide phenomenon that can have dramatic economic impact and sometimes tragically result in fatalities [1]. There are multiple factors that can make an area prone to landslide formation. Human factors such as drainage and land clearing (removal of vegetation) are some that can be controlled, and may expedite but not necessarily directly cause a landslide. Other environmental factors such as topography and the shear strength of the slope material are more difficult to control. Landslides are typically triggered by heavy rainfall events or sometimes by earthquakes [2,3]. The collection of data about existing landslides in a given area is important for predicting future landslides in that region [4].

Landslide monitoring requires continued assessment of the extent, rate of displacement, surface topography, and detection of fissure structures that could be related to fracture processes [5]. Measurements of vertical and horizontal displacements improve the understanding of landslide mechanisms [6,7]. Three dimensional (3D) measurements of landslides, such as the creation of Digital Surface Models (DSMs), allow volume displacements to be calculated. These displacements are related to the distance travelled by the landslide, which can improve our understanding of the mechanisms responsible for landslide dynamics and are useful for predicting future movements [8]. Technologies such as Differential GPS (DGPS), robotic total stations, airborne Light Detection and Ranging (LiDAR), and Terrestrial Laser Scanners (TLS) have revolutionized the periodic collection of DSM data since the early 2000s [9]. However, ground surveys with DGPS are time-consuming and have sparse spatial coverage, which results in the omission of fine-scale terrain structure in the resulting DSM [8]. TLS can suffer from line-of-sight issues, and airborne LiDAR is often cost-prohibitive for individual landslide studies [9]. However, none of these methods are suitable for real-time or near-real-time landslide monitoring.

There are many studies that have assessed the ability of various remote sensing techniques and data sources to monitor landslides. Large-scale landslides, such as those found in the Indian Himalayas, were monitored with satellite data from Cartosat (spatial resolution 2.5 m/pixel) by Martha et al. [8]. Several techniques, such as vegetation filtering and the use of GCPs, were applied to improve the accuracy of the data, and analysis of multi-temporal imagery allowed Martha et al. [8] to measure the movement of 550,000 m³ of material from the upper landslide areas. An example of using LiDAR is provided by Bell et al. [4], who mapped hundreds of landslides in Austria and Germany to study their persistence. A TLS was used by Pesci et al. [2] to map the crater of Mt. Vesuvius in Italy in 2005 and 2009. To obtain accurate results, careful analysis of the errors in the TLS data and point cloud alignment techniques were required [2].

Akca [6] used four fixed cameras to monitor a slope on which an artificially generated landslide was triggered. Photogrammetric techniques were then used to monitor the movement of the landslide, aided by a network of pre-measured Ground Control Points (GCPs); it was found that 103 m³ of material moved down the slope in the initial landslide event [6].

A new method for monitoring landslides is to utilize Unmanned Aerial Vehicles (UAVs), also known as Unmanned Aircraft Systems (UAS), to collect ultra-high resolution imagery. The use of UAVs for research purposes has become more commonplace in recent times due to technological developments such as autopilot systems, quality digital cameras, miniature GPS, and advances in lightweight carbon fiber airframes [10]. Several recent studies have also demonstrated the power of Structure from Motion (SfM) algorithms for landslide mapping and monitoring [5,11,12,13].

Immerzeel et al. [14] demonstrated that UAVs and SfM also provide a powerful tool for monitoring glacier movement and change, which in many respects are similar to landslides from a measurement and monitoring perspective. SfM is a relatively new image processing technique based on computer vision algorithms that were first developed in the 1990s. It allows the reconstruction of a photographed surface without the need for GCPs or complex pre-calibration of the camera [9]. This is possible because the position and the orientation of the camera are solved in an arbitrary space simultaneously with a highly redundant, iterative bundle adjustment that is based on thousands of extracted image features [9]. The resulting point cloud created by intersecting the matched features can then be transformed into a real-world coordinate system with a small number of GCPs. The mapping of three-dimensional terrain with SfM techniques has been demonstrated by authors such as Harwin and Lucieer [15], James and Robson [16], Ragg and Fey [17], Turner and Lucieer [18], Lucieer et al. [11], Lucieer et al. [19], and Westoby et al. [9].

Some studies, such as Chou et al. [20] and Bendea et al. [21], have taken advantage of the ability of UAVs to access unsafe landslide areas for the purposes of disaster management. Chou et al. [20] produced a 5 m resolution DSM from imagery, but this required GCPs, terrestrial photography, and camera calibration. Bendea et al. [21] collected imagery with a fixed-wing UAV and processed the data with specialized software (Leica Photogrammetry Suite) to produce georeferenced data to support humanitarian aid after natural disasters such as landslides. In both studies, the most significant limitation was that the maps they produced suffered from low accuracy due to the IMU and GPS data used.

In a geomorphological context, Niethammer et al. [5] used a Quadrocopter UAV platform, which is very similar in ability to the UAV platform used in this study (an Oktokopter), to create a map of a the Super-Sauze landslide (France) based on SfM techniques. Later Niethammer et al. [13] collected a second dataset of the same landslide. This allowed surface features to be tracked such that movement vectors could be calculated and comparison of DSMs allowed vertical displacements to be described. Niethammer et al. [13] suggested that UAV-based measurements of changes in landslide structure can provide an excellent data source for landslide modelers. Walter et al. [22] used a different approach, combining seismic monitoring with UAV remote sensing, and concluded that each type of sensor or technique mostly compensated for the other’s deficiencies. This technique was able to generate an orthomosaic of the landslide area, but suffered the limitation of requiring 199 GCPs.

Current literature and this study show that landslide monitoring with a UAV is a viable method. Previously, we demonstrated that highly accurate models of landslides can be created by the use of SfM techniques [11]. Other authors such as Niethammer et al. [13] have also proven the utility of UAVs combined with SfM for monitoring landslides. In this study, we build on our previous research by increasing the temporal coverage of our study and thus demonstrate that the methodology is repeatable and robust enough to be applied to a substantial time series of UAV datasets. The aim of this study is to accurately co-register multi-temporal UAV datasets, and to determine landslide surface dynamics and calculate volumetric differences from the multi-temporal DSMs. In addition, we examine and measure details such as change in slope and area of the landslide in the time series. Finally, this study aims to apply the image correlation techniques for surface motion detection that we first described in Lucieer et al. [11] to a multi-temporal dataset of UAV imagery.

2. Methodology

2.1. Platform

Multi-rotor UAVs are becoming more commonplace and are frequently used for commercial and recreational aerial photography. For this study we used an Oktokopter (eight rotors) multi-rotor micro-UAV (see Figure 1) with a DroidWorx carbon fibre airframe, a Mikrokopter autopilot (www.mikrokopter.com), and a Photoship One camera gimbal. The Oktokopter has a payload capacity of around 2 kg, a flight duration of 5–10 min (with a typical payload), and a stabilized camera mount to maintain nadir photos during the flight. Mikrokopter flight electronic systems are used to automatically maintain level flight, control the altitude, log system data at 1 Hz (including airframe position as measured with the on-board navigation grade GPS), and autonomously fly the UAV through a series of predefined, three-dimensional GPS waypoints.

Flights were pre-planned using a Google earth image to lay out a grid of waypoints spaced such that a high overlap (60%–80%) was maintained. Using previous ground survey data, a height profile was used to set the altitude of the waypoints such that a constant altitude above the landslide surface was maintained.

Figure 1. Oktokopter in flight, fitted with Canon 550D camera.

Figure 1. Oktokopter in flight, fitted with Canon 550D camera.

2.2. Sensor

To collect visible imagery, we used a Canon 550D Digital Single Lens Reflex (DSLR) camera (see Figure 1) (18 Megapixel, 5184 × 3456 pixels, with Canon EF-S 18–55 mm F/3.5–5.6 IS lens). Image capture rate was controlled by the UAV’s flight control board, which can be programmed to emit a trigger pulse at a desired frequency. The flight control board was connected to a custom-made cable that triggers the remote shutter release of the camera. The camera was operated in shutter priority mode (a fast shutter speed is required to minimize motion blur), in which the desired shutter speed (typically 1/1250–1/1600 s) was set before flight and exposure was controlled by the camera automatically by varying its aperture. Images were captured in RAW format and stored on the memory card in the camera for subsequent download post flight.

2.3. Field Site

In 1996 a landslide formed on the western slopes of the Huon valley in southern Tasmania, 35 km southwest of the capital city of Hobart (see Figure 2). The landslide formed on a steep slope in a cleared agricultural field next to the Home Hill vineyard and is approximately 125 m long and 60 m wide at an average elevation of around 80 m above sea level. McIntosh et al. [23] described the landslide as having developed in strongly weathered, layered fine colluviums, which are the remains of underlying Permian mudstone and siltstone that are estimated to be 4–5 m deep. The upper areas of the landslide are described as a rotational earth slide, which develops into an earth flow in the lower areas [23]. We have completed seven separate aerial surveys of the site, details of which can be found in Table 1.

Figure 2. December 2011 ortho-mosaic of Home Hill landslide area highlighting the main features. (a) location map; (b) example of GCP with actual location marked with a black cross giving a typical example of absolute spatial errors (coordinate system: GDA94 UTM55S).

Figure 2. December 2011 ortho-mosaic of Home Hill landslide area highlighting the main features. (a) location map; (b) example of GCP with actual location marked with a black cross giving a typical example of absolute spatial errors (coordinate system: GDA94 UTM55S).

Table 1. Details of aerial surveys of Home Hill landslide.

Survey Name	Date	Interval (Days)	Weather Conditions
2010A	20 July 2010	-	Sunny, light winds
2011A	19 July 2011	364	Overcast, light rain and wind
2011B	10 November 2011	114	Sunny, moderate winds
2012A	27 July 2012	260	Sunny, light winds
2013A	5 April 2013	252	Sunny, moderate winds
2013B	29 July 2013	115	Sunny, moderate winds
2014A	25 July 2014	361	Sunny, no wind

Table 1. Details of aerial surveys of Home Hill landslide.

Survey Name	Date	Interval (Days)	Weather Conditions
2010A	20 July 2010	-	Sunny, light winds
2011A	19 July 2011	364	Overcast, light rain and wind
2011B	10 November 2011	114	Sunny, moderate winds
2012A	27 July 2012	260	Sunny, light winds
2013A	5 April 2013	252	Sunny, moderate winds
2013B	29 July 2013	115	Sunny, moderate winds
2014A	25 July 2014	361	Sunny, no wind

2.4. Three-Dimensional Model Generation

Images collected during each aerial campaign were processed with commercial software, Photoscan Professional [24], that uses SfM techniques to reconstruct the scene based on a large number of overlapping photos. The flight path was pre-programmed into the UAVs autopilot to fly a grid pattern over the landslide at approximately 40 m above ground level. As the landslide varies in height by approximately 40 m from top to bottom, the UAV had to descend as it moved down the landslide to maintain a constant height above the landslide surface and thus maintain a constant spatial resolution.

Prior to flight, the camera’s internal clock was synchronized with GPS time, which allows the images to be geo-tagged with approximate coordinates from the on-board navigation-grade GPS. Images were collected at a rate of 0.75 Hz which oversamples the study area and provides a high level of redundancy in the dataset. It also allows images with excessive blur or tilt to be excluded. Blurry images were detected and removed via a blur metric as developed by Crete et al. [25], and a further description of its implementation can be found in Turner et al. [26]. Using the positional information for each image, it was also possible to remove the images captured during the UAVs ascent and descent. Furthermore, the positional information was used to find images taken from spatially similar positions, i.e., within close proximity. Images with close to 100% overlap do not provide a significant amount of extra data and thus were also removed.

Detailed descriptions of the Photoscan workflow can be found in Turner et al. [26], Lucieer et al. [11], and Lucieer et al. [19]. However, in summary, Photoscan initially detects tens of thousands of features in each image, which are then matched between the images. Using the matched features it is then possible to use an iterative bundle adjustment to estimate the positions of the matched features, positions, orientations, and lens distortion parameters of the cameras. This information is used for dense multi-view reconstruction of the scene geometry from the aligned images.

These processing steps are carried out in a real-world coordinate system based on the camera positions as supplied in the EXIF header, i.e., the geo-tagged positions as recorded by the UAVs on-board data logger. These positions are only recorded with a navigation-grade GPS receiver and hence are quite inaccurate (5–10 m absolute geometric accuracy). To improve the accuracy we manually identify GCPs within the imagery. The GCPs consist of metal discs spray-painted with fluorescent orange paint, scattered around the landslide area and then measured with a dual frequency RTK DGPS, providing GCP coordinates with an absolute accuracy of 2–4 cm. Based on these GCPs, Photoscan optimizes the bundle adjustment, which is followed by a dense geometry reconstruction. The final step is to export a DSM and orthophoto based on the dense 3D geometry (see Lucieer et al. [11] for a detailed description of the Photoscan workflow).

2.5. Alignment of Digital Surface Models

To prepare the DSMs for change detection, it is necessary to check the co-registration of each DSM pair. As each DSM has been georeferenced by separate workflows there is the potential for some minor mis-alignments, the most important of which are any differences in the Z-axis (or height). There is a possibility that there are also rotational or scale differences between a pair of DSMs. To check for such errors, the dense point clouds (from which the DSMs are created) for each model were exported from Photoscan and then imported into the Cloud Compare software package [27].

Using Cloud Compare, each pair of dense point clouds were initially masked for two reasons; firstly to ensure that both clouds covered the same extent, i.e., they were fully overlapping; and secondly the active landslide area was removed such that only the points in the surrounding non active areas were considered. An Iterative Closest Point (ICP) algorithm was then run with Cloud Compare on each pair of point clouds, and transformation matrices were estimated. These matrices include rotational parameters, translation parameters, and a scale parameter. For all the point cloud pairs, there was no rotational correction required and the scale factor was 1.0, implying that there were no rotational or scale differences between the point cloud pairs and thus the DSMs. Unfortunately, the translation parameters were inconsistent, as the ICP algorithm would in some cases suggest shifts of 20–30 cm in one or more of the three axis. This is theoretically due to the ICP converging to a local minimum in which a large translation minimized the error. Hence, it was necessary to develop another method to determine if there was any remaining constant translational bias between the datasets.

Immerzeel et al. [14], Lucieer et al. [11], and Martha et al. [8] employed a technique to validate the co-registration of multi-temporal DEMs by analyzing areas of the DEMs that were not subject to change. Similarly, in this study, we selected four separate areas from outside the active landslide zone, totaling 4,578,868 pixels (approximately 1830 m²), to further validate the alignment of the DSM pairs. These areas were then compared for each pair and the Root Mean Squared Error (RMSE) and volume difference were calculated (see Table 2). As these areas should not be changing, these values should theoretically be close to zero and have a preferably narrow Gaussian distribution of uncertainty/noise inherent to the image. To identify whether there was a constant bias between the DSMs, an offset was iteratively applied, and the corresponding RMSE and volume difference recorded until the optimal offset was found (see Table 2). Factors such as differing vegetation height, i.e., the grass may have been longer in one dataset than another, along with errors in the absolute positioning of the model, can contribute to this constant bias. Measurement of volume change does not require an accurate absolute location; an accurate co-registration of the pair of DSMs is of greater importance [8]. For our calculations, we considered the first dataset of each pair to be the reference landslide and then applied the offset to the second dataset to bring it in line with the first before volume change measurements were made.

Table 2. Evaluation of relative elevation accuracy for each pair of DSMs before and after offsets were applied, and differences between the DSMs were measured over 4,578,868 pixels (1830 m²) outside the active landslide area.

DSM Pair	Prior to Offset Application		Offset Applied (m)	After Offset Application
	Volume Difference (m3)	RMSE (m)		Volume Difference (m3)	RMSE (m)
2010A–2011A	153	0.109	0.10	114	0.077
2011A–2011B	88	0.061	0.00	88	0.061
2011B–2012A	98	0.074	0.00	98	0.074
2012A–2013A	134	0.108	0.09	115	0.085
2013A–2013B	148	0.101	0.13	102	0.087
2013B–2014A	76	0.059	0.00	76	0.059

Table 2. Evaluation of relative elevation accuracy for each pair of DSMs before and after offsets were applied, and differences between the DSMs were measured over 4,578,868 pixels (1830 m²) outside the active landslide area.

DSM Pair	Prior to Offset Application		Offset Applied (m)	After Offset Application
	Volume Difference (m3)	RMSE (m)		Volume Difference (m3)	RMSE (m)
2010A–2011A	153	0.109	0.10	114	0.077
2011A–2011B	88	0.061	0.00	88	0.061
2011B–2012A	98	0.074	0.00	98	0.074
2012A–2013A	134	0.108	0.09	115	0.085
2013A–2013B	148	0.101	0.13	102	0.087
2013B–2014A	76	0.059	0.00	76	0.059

2.6. Measurement of Landslide Area and Volume Change

To measure the overall mass displacement of the landslide, two co-registered DSMs were subtracted from one another. Before these measurements were made, however, an offset was applied to correct for the mismatch in height (Section 2.5). Once the DSM difference is calculated, areas of particular interest, e.g., advances of the landslide toes or retreat of the scarp, were segmented and the volume for each area was calculated. The height difference per pixel was multiplied by the area of a pixel (0.0004 m² for the 0.02 m resolution DSMs used in this study) and summed, thus giving a total volume in cubic meters.

The total area of the landslide was measured with the use of GIS software in which the active landslide area was manually digitized by visually identifying the landslide edge within the orthophoto that was generated by the Photoscan software (see Section 2.4). The area was then calculated by multiplying the number of pixels in the polygon by the area of a pixel, which was 0.01 m² for the orthophotos generated in this study.

2.7. Tracking of Landslide Surface Movement

Lucieer et al. [11] demonstrated that image correlation techniques can be used to track surface features between two DSMs of the landslide. In Lucieer et al. [11] we used an image correlation method developed by Leprince et al. [28] and Leprince [29]. The correlation method is implemented in the ENVI image processing environment [30] and is referred to as COSI-Corr: Co-registration of Optically Sensed Imaged and Correlation [31,32]. In Lucieer et al. [11] we experimented with various forms of input images and parameters for COSI-Corr and found that a hillshaded DSM was best for tracking surface features. COSI-Corr performed better on the hillshaded DSM than the color orthophotos due to large changes in illumination and vegetation coloring and structure between survey events.

This study builds onto our previous work in Lucieer et al. [11] by processing and analyzing a time series of seven UAV acquisitions of the Home Hill landslide. ENVI was used to initially produce hillshaded images from the DSMs for each of the seven dates with an input sun elevation of 45° and azimuth of 315°. These settings emulate the lighting conditions that are typical for a summer afternoon in the area and highlight the main terrain features as the azimuth is aligned with the main orientation of the landslide. COSI-Corr has a number of parameters to be selected. In this study, we used the same settings in Lucieer et al. [11], i.e., the statistical correlator with a window size of 64 pixels, a step size of 8 pixels and a search radius of 50 pixels (5 m), as these proved superior in motion detection.

3. Results

3.1. Accuracy of DSMs and Orthophotos

To assess the accuracy of the DSMs and the orthophotos created by Photoscan we typically used around 30% of the GCPs as check points (see Table 3). This means that approximately 70% of the GCPs were used for the bundle adjustment and the transformation to an accurate real-world coordinate system. The remaining GCPs were used as check points to assess the accuracy of the model, which was done by measuring the distance between GCPs in the orthophoto (for the x,y position) and the DSM (for the z position), and its correct location (as measured by RTK DGPS). For all the check points, an RMSE was calculated, and a summary of the results is presented in Table 3. Due to limited resources, we only had 23 GCPs for the 2011B dataset and thus all of these were required to transform the model. It can be seen that typical RMSE values are around 4–5 cm in the horizontal direction (XY) and 3–4 cm in the vertical direction (Z). The largest RMSE, 0.076 m and 0.09 m in XY and Z directions respectively, occurred in the 2013B dataset, which is possibly due to poorly measured GCPs and/or excessive shadowing on the landslide surface.

Table 3. Summary of spatial errors for Home Hill landslide DSMs and orthophotos.

Name	Date	Number Photos Used in Model	GCPs	Checkpoints	XY RMSE (m)	Z RMSE (m)
2010A	20 July 2010	62	56	19	0.046	0.031
2011A	19 July 2011	116	41	20	0.045	0.042
2011B	10 November 2011	194	23	23 *	0.021	0.025
2012A	27 July 2012	170	66	17	0.047	0.039
2013A	5 April 2013	179	29	22	0.058	0.078
2013B	29 July 2013	241	23	21	0.076	0.090
2014A	25 July 2014	415	16	10	0.031	0.031

* GCPs also used as check points due to limited resources.

Table 3. Summary of spatial errors for Home Hill landslide DSMs and orthophotos.

Name	Date	Number Photos Used in Model	GCPs	Checkpoints	XY RMSE (m)	Z RMSE (m)
2010A	20 July 2010	62	56	19	0.046	0.031
2011A	19 July 2011	116	41	20	0.045	0.042
2011B	10 November 2011	194	23	23 *	0.021	0.025
2012A	27 July 2012	170	66	17	0.047	0.039
2013A	5 April 2013	179	29	22	0.058	0.078
2013B	29 July 2013	241	23	21	0.076	0.090
2014A	25 July 2014	415	16	10	0.031	0.031

* GCPs also used as check points due to limited resources.

3.2. Area and Slope Analysis

The main area of change is the advancement of the two toes at the bottom and the retreat of the scarp at the top of the landslide (see Figure 2). The changes in landslide area are presented in Table 4. There were five events of the little toe advancing (a total of 554 m² during the whole monitoring period) and three scarp retreat events. There was no noticeable change in area between the 2013A and 2013B datasets, possibly because much of the movement occurs in surges rather than a gradual creep, and/or because movement during this short time period was too small to be measured.

The slope of the leading edge of the toes was measured by calculating the slope from the DSM and taking a mean of all the pixels in the leading edge area (see Table 4). The steepness of the leading edge gradually builds up, particularly of the large toe, as material flows down from above until eventually the leading edge collapses and the toe surges forward. This surge effect can be seen graphically in Figure 3. There was little or no forward movement of the large toe between the 2010A and 2012A surveys, but internal shear stress was building up; then, sometime before the 2013A measurement, the leading edge gave way and the large toe advanced around 1–2 m (contributing to the 126 m² of area increase at this time) and the slope reduced in steepness. It can be seen in 2014 that the slope is beginning to steepen again, and it is likely the cycle will continue as material continues to flow down the slope. The trend is not quite so obvious for the small toe, most likely because this area underwent some large changes during the monitoring period (it advanced around 12 m between 2010 and 2011) and has thus not settled into the same pattern as the large toe, i.e., the deformation pattern is more chaotic. In 2014 the slope is becoming steep again and apparently internal pressure and stress are building up. Based on past dynamics, it is likely that the toe will surge forward again soon.

The other significant changes in landslide area came from the three retreats of the scarp, of which the first occurred between July 2011 and November 2011 when there was a large collapse (162 m²) of the northern part of the scarp. Much of this material seems to have flowed down the slope and contributed to the major advance of the small toe. Between July 2012 and July 2013, a further 95 m² of the northern part of the scarp collapsed, followed by a further collapse of 47 m² in 2014, making this a highly dynamic area of the landslide. There were no significant collapses in the southern area of the scarp during the monitoring period, despite the fact that it is also a near vertical face of similar height.

Figure 3. Transects through leading edge of the large toe.

Figure 3. Transects through leading edge of the large toe.

Table 4. Changes in area and leading edge slope of Home Hill landslide.

Name	Total Area (m2)	Slope of Large Toe (Deg)	Slope of Small Toe (Deg)	Area of Toe Advance (m2)	Area of Scarp Retreat (m2)
2010A	4887	31.05	36.26	-	-
2011A	5168	33.72	34.92	281	-
2011B	5435	34.37	34.06	105	162
2012A	5455	39.98	36.22	20	-
2013A	5675	34.17	34.78	126	95
2013B	5675	33.17	34.54	-	-
2014A	5744	33.87	37.63	22	47

Table 4. Changes in area and leading edge slope of Home Hill landslide.

Name	Total Area (m2)	Slope of Large Toe (Deg)	Slope of Small Toe (Deg)	Area of Toe Advance (m2)	Area of Scarp Retreat (m2)
2010A	4887	31.05	36.26	-	-
2011A	5168	33.72	34.92	281	-
2011B	5435	34.37	34.06	105	162
2012A	5455	39.98	36.22	20	-
2013A	5675	34.17	34.78	126	95
2013B	5675	33.17	34.54	-	-
2014A	5744	33.87	37.63	22	47

3.3. DSM Volumetric Changes

An example of the DSM change is presented in Figure 4, highlighting typical changes on the Home Hill landslide. There was a loss of material where the scarp has collapsed (see Figure 4 near (a)) and an accumulation of material at the leading edges of the toes (see Figure 4 near (b)). The spatial pattern of mass losses and gains on the landslides matches intuitive expectations of movements on such a landslide. The pattern is typical and can be found for most of the DSM comparisons where the scarp has collapsed and/or the toes have advanced. Subsequent comparisons were then made for the other DSMs. From these maps, areas of interest (e.g., the collapsed scarp, see Figure 4a for an example) were selected and volumetric change was calculated as described in Section 2.6. A summary of the volumetric changes for the small toe of the landslide is presented in Table 5, which also presents the bulking factor, the ratio of the volume gain to volume loss. For the first two comparisons, the bulking factors are typical of what might be found on other landslides with similar characteristics [8]. However, in the remaining comparisons, the values are higher than one would expect.

Table 5. Volumetric changes in m³ for the small toe of Home Hill landslide, including estimated errors at one sigma; also reported is the bulking factor, the ratio of accumulated material to lost material.

Name	Toe Accumulation and Estimated Error (m3)	Loss above Toe and Estimated Error (m3)	Bulking Factor
2010A → 2011A	572 ± 24	340 ± 14	1.68
2011A → 2011B	249 ± 21	124 ± 17	2.00
2011B → 2012A	88 ± 19	24 ± 10	3.66
2012A → 2013A	175 ± 21	41 ± 9	4.27
2013A → 2013B	no change	no change	-
2013B → 2014A	85 ± 17	22 ± 8	3.86

Table 5. Volumetric changes in m³ for the small toe of Home Hill landslide, including estimated errors at one sigma; also reported is the bulking factor, the ratio of accumulated material to lost material.

Name	Toe Accumulation and Estimated Error (m3)	Loss above Toe and Estimated Error (m3)	Bulking Factor
2010A → 2011A	572 ± 24	340 ± 14	1.68
2011A → 2011B	249 ± 21	124 ± 17	2.00
2011B → 2012A	88 ± 19	24 ± 10	3.66
2012A → 2013A	175 ± 21	41 ± 9	4.27
2013A → 2013B	no change	no change	-
2013B → 2014A	85 ± 17	22 ± 8	3.86

There are likely to be some inaccuracies in the volumetric measurements due to minor misalignments of the DSMs and environmental differences between datasets, such as the length of the grass. The alignment processes described in Section 2.5 minimized these alignment errors, but we wanted to quantify the error of the volumetric measurements. The average RMSE of the difference between the non-active areas of the landslide was 0.07 m (see Table 2), and the mean standard deviation of these differences was also 0.07 m. The errors expressed in Table 5 are thus expressed as ±0.07 m (one sigma).

Figure 4. 2012A to 2013A DSM change: (a) area of scarp collapse; (b) areas of toe advance (coordinate system: GDA94 UTM55S).

Figure 4. 2012A to 2013A DSM change: (a) area of scarp collapse; (b) areas of toe advance (coordinate system: GDA94 UTM55S).

3.4. Historical DSM

To model the terrain before the landslide occurred, we obtained a pair of overlapping historical aerial photographs of the Home Hill area from the archives at the Department of Primary Industries, Parks, Water and Environment (www.dpipwe.tas.gov.au). The images were collected in 1984 and had a scale of 1:15,000. The images were scanned at 2000 dots per inch, which gave a ground resolution of approximately 20 cm/pixel. The digital scans were imported into Photoscan and a 3D model was constructed (see Figure 5). Ground features that were identifiable in the 1984 images and in the latest Google Earth imagery (e.g., edges of roads, corners of roofs, and corners of fence lines) were used as GCPs to enable georeferencing of the model.

Figure 5. 3D model of the Home Hill landslide as constructed from 1984 aerial photography. The landslide area is circled in red.

Figure 5. 3D model of the Home Hill landslide as constructed from 1984 aerial photography. The landslide area is circled in red.

The georeferencing accuracy of the Google Earth imagery is insufficient for the 1984 model to be accurately co-registered with the recent UAV datasets. In a similar manner to the method used to align the DSMs (Section 2.5), we aligned the 1984 model with a recent dataset from July 2013. The alignment was carried out on the point cloud data rather that the DSMs, because both horizontal and vertical misalignments had to be corrected. Cloud Compare [27] was used for the point cloud co-registration using the ICP algorithm [33], which minimizes the distance between two point clouds. This process considered the 2013 dataset the true one and moved the 1984 point cloud as close to the 2013 point cloud as possible. Once a visual validation of the results was completed, the LAStools software suite [34] was used to interpolate the 1984 point cloud into a DSM.

The 1984 DSM is based on much lower resolution imagery than the data captured with the UAV and thus the final 1984 DSM only had a resolution of 50 cm/pixel, requiring the 2013 dataset to be resampled to match this resolution. Using the methods described in Section 2.6, volumetric calculations were made. The approximate values for total amount of material moving down the landslide are 6630 m³ of accumulation in the lower areas of the landslide and 3300 m³ of material lost from the upper areas, producing a bulking factor of 2.02, which is typical for this type of landslide [8].

3.5. Surface Movement

As described in Section 2.7 the COSI-Corr image correlation software was used to compare DSM datasets to track the movement of surface features. In Lucieer et al. [11], COSI-Corr was applied to the 2011A and 2011B datasets to create one map of surface movement. We now have three more data sets (2012A, 2013A, and 2014A) to which COSI-Corr was applied, expanding the dataset to four maps of surface movement (see Figure 6A–D).

Figure 6. Cosi-Corr surface movement maps; (a) 2011B–2012A; (b) 2012A–2013A; (c) 2013A–2013B; and (d) 2013B–2014A (coordinate system: GDA94 UTM55S).

Figure 6. Cosi-Corr surface movement maps; (a) 2011B–2012A; (b) 2012A–2013A; (c) 2013A–2013B; and (d) 2013B–2014A (coordinate system: GDA94 UTM55S).

In Figure 6, A–D illustrate the irregular behavior of the landslide movements between the different image acquisition dates even if one accounts for the variable numbers of days between image takes (Table 1). These observations may confirm the “surge behavior” of the landslide. The small toe at the north-east side of the landslide is clearly the most active part. Registered displacement ranges up to 7 m for the 2012A–2013A data taken with an interval of 252 days as shown in Figure 6B. The large toe, situated in the southern part of the landslide, is less dynamic and only shows a major displacement in Figure 6C for the 2013A to 2013B image period. Displacements range up to 4 m. Displacement registration at the scar is difficult for the COSI-Corr algorithm since the movements are nearly vertical, which makes them difficult to identify in UAV-based horizontal images. Furthermore, scar movements often show rotational movement hampering the identification of surface features required by COSI-Corr to compute the displacements. The scar shows clear signs of surface instability and movements in Figure 6B,C. Detached material at the scar is required to activate or re-activate the movements of the toes.

3.6. Comparison of Landslide Movement with Rainfall

The lower areas of the Home Hill landslide are described as an earth flow by McIntosh et al. [23], and we have hypothesized that movement of the landslide is triggered by extreme rainfall events, as precipitation is one of the main landslide triggering factors [35]. We also hypothesize that the rate of the earth flow may be partially related to the amount of rainfall that the landslide receives in the preceding days or weeks. To test this hypothesis, rainfall data from the two nearest rainfall stations, Huonville (3.8 km away to the south) and Grove (4.75 km away to the north) were analyzed. The total annual rainfall between each July landslide dataset, i.e., 2010A, 2011A, 2012A, 2013B, and 2014A was calculated and averaged across the two rainfall stations and then translated into a mean daily rainfall for that year. The annual rainfall rate was then compared with mean daily volume of movement for the small toe area (see Figure 7).

Figure 7. Comparison of annual rate of rainfall and annual rate of accumulation of material for the little toe area of the Home Hill landslide.

Figure 7. Comparison of annual rate of rainfall and annual rate of accumulation of material for the little toe area of the Home Hill landslide.

4. Discussion

To allow accurate and detailed comparisons of landslide changes and movement, it is important that the errors in the alignment of all the datasets be minimized. The creation of 3D models from UAV imagery has a limitation in that it does not create a true Digital Elevation Model (DEM) of the surface; instead it created a DSM. This is because the technique is based on a passive optical sensor that cannot penetrate the vegetation and find the true surface height, as a LiDAR system might. However, we have implemented a method to minimize the alignment error between datasets caused by factors such as the variable height of the grass.

Simple measurements, such as change in area and slope of the landslide, can be informative. These data can be used to infer patterns of movement, such as the cycle of leading edge slope change (see Section 3.2). Increasing the temporal coverage of the dataset, i.e., increasing the frequency of measurement and continuing to collect regular datasets into the future, could confirm the validity of these cycles.

From the data collected, the actual calculation of volumetric change was a trivial matter; however, it does suffer from a degree of uncertainty. This is caused by multiple factors that contribute to the absolute accuracy of the models, resulting in an error in the Z-axis. However, when areas of volumetric change are large, these inaccuracies will result in volumetric errors that only represent a small percentage of the total volume measured (see Table 5). Conversely, when measuring subtle changes, these errors will represent a high percentage of the total volume measured (see Table 5), i.e., the absolute error is small but the relative error is large.

Applying the same surface reconstruction methods to historical aerial photography can provide an approximate baseline dataset for the original landform prior to the formation of the landslide. Acquisition of further aerial photographs from the past has the potential to increase the temporal coverage of this landslide monitoring study.

Advanced image correlation techniques provide an automated method for tracking surface movements of the landslide, giving us an insight into the rate at which various areas of the landslide are moving and also which areas remain static. There is one area, denoted as “the island” (see Figure 6a), which does not move at all. Inspection of the surface material in this area revealed a conglomerate of rocks which is likely to be holding this area together while the landslide flows around it (like a nunatak in a glacier). The large toe is clearly steadily flowing down hill (see Figure 6) and seems to have settled into a pattern of the leading edge building up until it collapses and surges forward (see Figure 3). However, the little toe area is more dynamic (see Figure 6), with major forward surges. These surges seem to have reduced the material supporting the northern part of the scarp and thus it has collapsed on two separate occasions, providing more material to flow downhill on the little toe.

When measuring volumetric changes, we also calculated the bulking factor for accumulation areas of the landslide. Bulking factors up to a value of 2.0 are known to be typical for landslides of this type [8]. However, some very high values (up to 4.3) were found, for example the 2013A comparison (see Table 5), and the COSI-Corr movement maps can help explain these high values. The surface changes between 2012A and 2013A (see Figure 6c) show a great deal of movement above the small toe, so it is likely the area of loss was filled from above. This would reduce the measurement of volume loss and result in a high value for the bulking factor.

Analysis of rainfall data and comparison with landslide movement (see Figure 7) has shown a strong relationship for the first three annual comparisons. However, the relationship is not evident in the fourth comparison, and this may be due to several factors. The leading edge of the little toe in the 2014 dataset was very steep (see Table 4), implying that it is likely to surge forward again soon. The small toe has now extended to a similar part of the slope as the large toe, that is, an area that is not as steep, and it is also starting to fan out in a similar fashion to the large toe; these factors would slow the rate of movement (compared to the previous years) and thus reduce the volume of material moving downslope.

There is further evidence of rainfall being related to landslide movement when we examine the historical records from 1996, the year in which the landslide occurred. In 1996 the monthly rainfall for January was 194% of the long term mean; in February it was 268%; March, 166%; and in April the rainfall was 338% of the long term mean. It may very well be that this excessive rainfall triggered the landslide. We know from the 1984 imagery that the area had been cleared of vegetation, meaning that for at least the 12 years prior to 1996 the land was cleared. It is thus unlikely that the land clearing alone was the cause of the landslide. Further detailed ground study is required to better understand the structure of the landslide material and to assess its failure characteristics, which will be the subject of future research.

5. Conclusions

A significant achievement of this study was to demonstrate that it is possible to generate accurate Digital Surface Models (DSMs) of a landslide with an Unmanned Aerial Vehicle (UAV) and that this technique is robust and repeatable such that a substantial time series of datasets can be routinely collected. Another achievement was the development of a technique to check and correct for any co-registration errors between subsequent DSMs based on comparing non-active areas of the landslide, minimizing the alignment error to ±0.07 m on average. Measurements of variables such as landslide area, leading edge slope, and volumetric changes were made over a time series of seven DSMs spanning four years. The COSI-Corr: Co-registration of Optically Sensed Imaged and Correlation [31,32] image correlation algorithm was used to track and quantify surface movement vectors of the landslide. It should be noted that there is no ground validation for the movement vectors generated by COSI-Corr; however, we have demonstrated that COSI-Corr is robust and works over a series of datasets collected with the same methodology.

The main new finding of this study was to prove that a time-series of UAV images can be used to map landslide movements with centimeter accuracy. It also found that there can be a cyclical nature to the slope of the leading edge of the landslide, suggesting that the steepness of the slope can be used to predict the next forward surge of the leading edge. It was also discovered that historical aerial photography could be used to create a DSM as a baseline dataset for comparison to the modern DSMs. This revealed that a total of approximately 6630 m³ of material had moved downslope since the formation of the landslide. Evidence was also found that there is a relationship between the amount of rainfall that the landslide receives and the volume of material that flows downslope.

The methodology described in this paper advances knowledge in that it allows a series of measurements to be made of landslide variables. These measurements can assist experts in monitoring and understanding landslide movement dynamics. Whilst there have been other studies that have mapped surface movement of landslides and glaciers with UAV imagery [5,11,12,13,14], this study is unique as it has presented a methodology to accurately co-register DSMs and to perform a deformation analysis of a time series consisting of seven separate datasets. It also builds upon our previously published research into the use of image correlation techniques to track surface movement, demonstrating that the technique is robust and can be applied to a longer time series of images.

The way forward for this research will involve continued monitoring of the landslide to further validate the findings in this study, e.g., the relationship between the leading edge slope and movement of the leading edge, and also the relationship between rainfall and volume displacement of material on the landslide. It would be advantageous to be able to sample the landslide at a higher frequency, which will require improvements in the methodology to minimize the time required in the field. This could be achieved by implementing a direct georeferencing system, as described in Turner et al. [26], to remove the requirement for Ground Control Points (GCPs).