Use of High-Resolution Multi-Temporal DEM Data for Landslide Detection

Abstract

Landslides in urban areas have been relatively well-documented in landslide inventories despite issues in accuracy and completeness, e.g., the absence of small landslides. By contrast, less attention has been paid to landslides in sparsely populated areas in terms of their occurrences and locations. This study utilizes high-resolution and LiDAR-derived digital elevation models (DEMs) at two different times for landslide detection to (1) improve the localization and detection accuracies in landslide inventories, (2) minimize human intervention in the landslide detection process, and (3) identify landslides that cannot be easily documented in the current state of the practice. To achieve this goal, multiple preprocessing steps were used to ensure the spatial alignment of the multi-temporal DEMs. Map algebra was then used to calculate the vertical displacement for each cell and create a DEM of Difference (DoD) to obtain a quantitative estimation of ground deformations. Next, the elevation changes were filtered via an appropriate Level of Detection (LoD) threshold to mark potential landslide candidates. The landslide candidates were further assessed with the aid of customized topographic maps as auxiliary data and pattern recognition to distinguish landslides (true positive changes) from construction, erosion, and deposition (false positives). The results from the proposed method were compared with existing landslide inventories and reports to evaluate its performance. The new method was also validated with temporal high-resolution Google Earth images. The results showed the successful application of the method in landslide detection and mapping. Compared with traditional methods, the proposed method provides a semi-automatic way to obtain landslide inventories with publicly available yet lowly utilized DEM data, which can be valuable in preliminary analysis for landslide detection.

1. Introduction

Landslides are among the most consequential natural hazards with devastating outcomes such as human casualties as well as damages to civil infrastructure, environment, and cultural heritage [1]. This makes landslides an ongoing concern and demands new methods to better understand the mechanisms of landslides and promote the sustainable use of ecosystems for disaster risk reduction [2]. The first step in mitigating the risk associated with landslides is compiling a landslide inventory to provide systematic information about landslides [3]. The main challenge in compiling these inventories is the landslide detection approach. There are different approaches for landslide detection, which are mainly focused on two types of data, i.e., satellite imagery and surface morphology. These approaches are utilized to detect and map landslides mostly manually via the visual interpretation of human experts [4].

In satellite imagery, high-resolution images were used to find temporal changes in land cover, structure, and lithology properties, indicating landslide occurrence [5]. Although mapping landslides based on land cover change is a well-established approach [6], the accuracy of such methods diminishes with more extended time intervals due to the vegetation growth and revegetation, which covers the scars [7]. Moreover, applying these methods to wildfire-prone areas, such as California, is complicated and could prove to be inaccurate due to wildfire-induced land cover changes [8]. However, investigators have adopted different types of satellite imagery to identify landslides via land cover change detection. For example, optical sensors have long been used for mapping large landslides through temporal change detection to identify variations in soil cover or vegetation [9]. New very high resolution (VHR) sensors [10] and synthetic-aperture radar (SAR) [11] are other examples of satellite imagery in landslide detection. The main drawback of satellite imagery is the complex manual preprocessing steps needed to prepare the data for interpretation, e.g., pansharpening, orthorectification, and radiometric correction [12]. To detect and map landslides using different types of satellite imagery, visual interpretation or semi-automated procedures were performed by a human expert to identify the above-mentioned temporal changes in the study area. The main semi-automatic approaches for using satellite imagery for landslide detection purposes are object-based and pixel-based. The more widely-used approach is pixel-based, which includes traditional maximum likelihood or minimum distance methods or more recent machine learning models such as decision trees, artificial neural networks, convolutional neural networks, and random forests [13]. These approaches have issues with high-resolution satellite imagery and its rich information content [14]. Object-based analysis, however, uses a segmentation process that groups pixels into meaningful units instead of analyzing them individually. This enables researchers to utilize spectral information, textural data, and spatial properties of data that were not possible before [15]. For example, Ghorbanzadeh et al. [16] used Sentinel 2 data and ALOS DEMs with fully convolutional networks, i.e., U-Net and ResU-Net, to detect landslides.

In addition to satellite imagery, surface morphology has been extensively used to map landslides. In surface morphology analysis, digital elevation models (DEMs) of the surface and their derivatives, such as shaded reliefs and surface roughness, allowed the mapping of complex landslides [17]. These derivatives were used to identify patterns in morphology, movement history, material types, and topographic variability to detect landslides through the interpretations of a trained geomorphologist [18]. The interpreters search for geomorphological features that indicate landslide activity, such as roughness characteristics, major and minor scarps, and cracks [19]. For example, McKean and Roering [20] utilized statistical, Laplacian, and spectral analysis of DEMs to obtain spatial patterns of roughness to identify landslides and their morphologic domains. LiDAR (Light Detection and Ranging) revolutionized these DEM-related practices and helped create high-resolution DEMs of the surface with sub-meter accuracy [21]. In addition to high accuracy, LiDAR can penetrate vegetation [22] and forest canopies and create an accurate estimation of the surface elevation [23]. The major drawback of using LiDAR data is the time-consuming and expensive data acquisition process for large areas [12]. Despite disadvantages, the employment of LiDAR-derived DEMs for obtaining an accurate map of landslides in a small area can aid in understanding the landslide process and is used in data-driven methods to identify landslides for other purposes or future applications [24].

There are some issues with the methods mentioned above. First, the number of highly-skilled and experienced human experts essential for visual interpretation is declining [12]. This fact hampers our ability to produce high-quality landslide maps, which subsequently makes it harder to validate new landslide detection methods. Second, compiling a landslide inventory for a large area is a time-consuming and expensive procedure that requires technological and economic resources to acquire and process the data [25]. Third, the accuracy of the conventional methods, which rely on changes in morphology and land cover diminishes with longer time intervals, especially in areas that experienced wildfires which causes physical weathering and vegetation removal. Fourth, landslide inventories are subjective [26], and their quality depends on the accuracy of the data, the expertise of the interpreter, geological and morphological complexity, and many other factors [27]. These issues prompted researchers to develop quantitative approaches to process multi-temporal DEM data, i.e., DEMs generated from data acquired two or more times, for landslide assessment studies. A popular type of study is to analyze landslide boundary changes to monitor the morphologic and volumetric changes in an active landslide over short periods [28,29]. Such studies were often focused on understanding the evolution and kinematics of specific deep-seated and well-known landslides [30,31] rather than identifying new landslides [32]. Another application of multi-temporal LiDAR data was to identify ground deformations and landslides in a specific study area. For example, Burns et al. [33] used temporal DEMs, site reconnaissance visits, and satellite imagery to detect landslides in Western Oregon; Mora et al. [34] employed probabilistic change detection with temporal DEMs to map surface changes and determine the probability of them being landslides. However, these studies used new LiDAR surveys, either collected by themselves or by a third party, to test their hypothesis instead of using publicly available DEMs. So far, no study has attempted to detect landslides for a large area with publicly available DEM data.

To improve the current practice of landslide detection, a new semi-automatic and nonsubjective method is presented in this paper. This semi-automatic feature means that the proposed method does not require a human expert interpreter, which translates into greater consistency and faster results. The nonsubjective feature means that the performance of the new method depends solely on the input data, and factors such as the expertise of the interpreter and the length of the time interval do not affect the outcomes. The reliance on data also means that the new method analyzes sparsely populated areas with the same procedure and accuracy as it treats urban areas, which is beneficial for the existing landslide inventories that are biased toward documenting landslides that cause socio-economic disruption and casualties. Additionally, this study aims to take advantage of the mentioned benefits of high-resolution LiDAR-derived DEMs while avoiding their time-consuming and costly acquisition processes as their main disadvantage. For this purpose, we used publicly available multi-temporal DEMs from the 3D Elevation Program (3DEP), managed by the U.S. Geological Survey (USGS). Given its low cost and the use of publicly available data, this method can be a useful tool for preliminary analysis in landslide detection projects. To prepare the DEM datasets for map algebra operations, the DEM data were preprocessed to ensure their spatial alignment. Next, map algebra was employed to detect vertical displacements within the study area and create a differential DEM called the DEM of Difference (DoD). Then, the areas with an elevation change above a certain threshold were marked as landslide candidates. The customized vector basemaps and multidirectional Hillshade maps were employed as auxiliary data to infer the causes of elevation changes to divide these landslide candidates into landslides (true positives) and “nonlandslides” (false positives). The nonlandslide category includes new construction, excavation, fluvial settlements, and erosion. The performance of the proposed method was evaluated with landslide inventories and publicly available Google Earth images to verify its effectiveness.

2. Workflow and Study Area

This study comprises three main parts (Figure 1). The first part details the theoretical foundation of this study. The DEM data are presented in two groups: 1-m DEMs from 2018 and 1/9 arc-second DEMs from 2010. Each group contains a series of topographic quadrangles, i.e., tiles, that cover the study area. The tiles in each group have the same coordinate system and cell resolution. Multiple steps were taken to preprocess the DEMs and prepare them for use, i.e., DEM treatments. The tiles in each group were merged into a single raster using geoprocessing tools to aid in later visualization and analysis. Next, the coordinate system of the 1/9 arc-second DEM dataset was altered through map projections to match that of the 1-m DEMs. After that, the DEM datasets were resampled to acquire similar resolution DEMs with aligning grid cells. In the second part, map algebra [35] was used to obtain a DoD from the resultant DEMs. The DoD was then analyzed and filtered via a Level of Detection (LoD) threshold to distinguish information (actual vertical displacement) from noise. Auxiliary maps were then used to better interpret the differences between the two raster datasets. The purpose of these maps was to help distinguish the changes in the landscape caused by landslides (true positives) from those caused by constructions, changes in the river basin, hydro-flattening effects, and erosion (false positives). The third part was intended to verify and compare the results of the proposed method using two landslide inventories and Google Earth images. The following two sections will present more technical details for the theoretical basis and method.

The study area is located in the state of California in the USA. It includes parts of San Francisco and Marin Counties that cover more than 2000 km². Figure 2 displays the location of the study area using a red rectangle. This area is marked with dense vegetation, and its elevation ranges from 104 m below sea level to 785 m above sea level, with a mean and a standard deviation of 120 m and 117 m, respectively. The local climate is Mediterranean, and precipitation falls mainly from January to May and from October to December.

3. Theoretical Basis

3.1. Data Description

The advances in remote sensing, particularly LiDAR technology, have provided means for acquiring higher-resolution imagery. The ability of LiDAR to penetrate vegetation, cover large areas, and obtain precise surface models renders it an enticing tool for landslide detection and mapping studies [36]. Moreover, the repeatability of high-resolution DEM data provides avenues to extract spatial information, such as temporal changes in the landscape. To harness this potential, a set of techniques and conventions, i.e., map algebra developed by Tomlin and Berry [35], are used to manipulate and analyze cartographic data. The cartographic data consist of perpendicular and equally spaced grids that define locations so that each location is uniquely associated with a square called a “grid cell” or a “pixel” in a raster. Therefore, the overall scheme for a raster is a rectilinear grid of columns and rows, and each cell is a discrete geometric entity representing a recorded characteristic/variable, e.g., elevation [37]. Map algebra operations use these variables to calculate and output new variables and extract meaningful information from the data.

In the original definition of map algebra, locations are expressed in a Cartesian $\left(x,y\right)$ coordinate system in a two-dimensional space. In three-dimensional space, locations are shown with cubes instead of squares and form cubic zones and layers. The rest is quite similar to the 2D space; hence the existing functions can be extended to the 3D space with $\left(x,y,z\right)$ coordinates. Moreover, the 3D map algebra can be extended to the temporal dimension by regarding the $z$ dimension as temporal in nature [38]. This assumption allows us to combine spatio-temporal values and extract information such as temporal topographic changes. To achieve this, we discretized the temporal dimension and evaluated the elevation values in two snapshot states representing two points in time. A prerequisite for performing such operations is the spatial alignment of the grids across the temporal dimension. This is noteworthy because any misalignment in temporal grids undermines the assumption that the locations are shown by cubes, consequently introducing errors in the map algebra operations. That is why quantitative analysis of changes requires special care to ensure the spatial alignment of the DEM datasets. The preprocessing steps used to address these misalignments include coordinate system transformations, cell resolution adjustments, and grid alignment adjustments. More details regarding these preprocessing steps can be found in Section 4.2. Knowing the background processes of these preprocessing steps is essential because the increasing access to publicly available DEM data and GIS software has made modeling more accessible than ever; however, using this DEM data without a proper understanding of the digital database is of little practical significance. The sensitivity of the DEM data to grid cell resolution and alignment is often overlooked due to the misconception that digital spatial data are scaleless. This results in an unrealistic sense of model accuracy that is detrimental to the modeling outcomes [39]. That is why we adopted well-established resampling methods to manipulate and prepare the DEM datasets for use.

3.2. Challenges

After ensuring spatial alignment, arithmetic operations can be utilized to quantify terrain changes and create a DoD for two-time snapshots, $t_k$ and $t_{k+1}$. In other words, the elevation values of the older DEM were subtracted from the elevation values of the newer DEM on a cell-by-cell basis. To create a DoD, spatially corresponding cells of the input grids, $z_{\left(k+1\right)}\left(i_c,j_c\right)$ and $z_{\left(k+1\right)}\left(i_c,j_c\right)$, were subtracted from one another and stored in a matching output grid, $DoD\left(i_c,j_c\right)$:

$$DoD\left(i_c,j_c\right)=z_{k+1}\left(i_c,j_c\right)-z_k\left(i_c,j_c\right)$$

(1)

These DoDs have many advantages and can be utilized to identify and estimate ground deformations and study their causes. DoDs have been used to map erosion, deposition, calculate volumetric changes, and other applications [40,41]. The use of DoDs gives rise to a set of nuanced questions: “What types of errors are involved in the DoD products and how to mitigate them?” and “What is the root cause of the ground deformations obtained by DoDs?”

Regarding the first question, several studies have tried to obtain the vertical errors of LiDAR and other high-resolution DEM data for different conditions using methods such as geodetic control networks [42]. Although there is no unanimous agreement on the exact numbers, studies suggest that the error is somewhere between 17 cm to 26 cm, depending on the vegetation cover [43]. However, issues arise when arithmetic operations are performed to derive products from DEMs such as DoD. This fact is due to the difference in acquisition tools and methods and is exacerbated when using multi-temporal DEMs spanning long periods. Preprocessing the DEM data is pivotal in such cases to ensure that issues, including the difference in coordinate systems, resolution, and misalignment of DEM grids, do not introduce unnecessary errors in the analysis. Appropriate preprocessing techniques can help minimize the DoD errors to inevitable ones, as described in the following. Mathematically, a simple error propagation theory can be used to analyze these errors for a derived variable such as DoD, provided that the two processed DEMs are independent [44]. This means that, if the errors for the two DEMs involved in DoD calculations are ${\delta }_{D1}$ and ${\delta }_{D2}$, the total estimated error is ${\delta }_{DT}:{\delta }_{DT}=\sqrt{\left({\delta }_{D1}{}^2+{\delta }_{D2}{}^2\right)}$ [45]. Any elevation change below this propagated error value will be filtered out as they are indistinguishable from noise. This is the primary reason for applying a minimum LoD threshold.

The second question is also crucial since not all ground deformations are due to landslide activities and natural processes. Human geomorphic activities, including urban development, road construction, mining, and landfills can lead to topographic changes. Furthermore, natural processes such as erosion and deposition can also contribute to such changes. This is especially important for DoDs that cover extended periods in which the accumulation of seemingly insignificant natural processes results in measurable changes in elevation. Several postprocessing steps are needed to classify these changes into landslide and nonlandslide related. Choosing an appropriate LoD threshold is the first of these steps. The LoD is a filter that aids in obtaining an abstract form of representation for the DoD layer by classifying each cell into one of the two discrete classes to distinguish noise from real changes:

$$L_{\left(i,j\right)}=\{\begin{array}{cc}1\hfill & \hfill \mathrm{if} Do{D}_{\left(i,j\right)}>\mathrm{LoD}\\ 0\hfill & \hfill \mathrm{otherwise} \end{array}$$

(2)

where $Do{D}_{\left(i,j\right)}$ is the elevation change for cell $\left(i,j\right)$ and $L_{\left(i,j\right)}$ is the resultant raster dataset.

Another critical attribute of the identified elevation changes is their relative location to built-up areas and topographical features such as roads and rivers. This relative location aids in identifying other nonlandslide-related ground deformations, including urban development, fluvial settlements, erosion, and construction. Additionally, visual pattern recognition of contiguous subsidence (downward vertical movement) in upslope and accumulation in downslope can help characterize elevation changes as landslides.

4. Method

4.1. DEM Acquisition

The DEM data were acquired from the 3D Elevation Program (3DEP) of The National Map (TNM), managed by the U.S. Geological Survey (USGS) [46]. The 3DEP provides various open-access products and services, including DEMs with different resolutions [47]. This study utilized the two highest resolution DEMs available in this program: 1-m DEMs acquired in 2018 and 1/9 arc-second DEMs acquired between April 2010 and March 2011. The 1-m DEMs were exclusively produced from LiDAR source data with the Universal Transverse Mercator (UTM) as the spatial reference with the North American Datum of 1983 (NAD83) [48]. The 1/9 arc-second DEMs have the second-best resolution (approximately 3 m). These DEMs also used NAD83, but they were distributed in a geographic coordinate system in the unit of decimal degrees. The 1/9 arc-second DEMs were developed from LiDAR, photogrammetry, and other high-resolution sources [49,50]. Elevation data in both of these products represent the hydro-flattened bare-earth topographic surface [47]. These DEM sets were both obtained from 3DEP in multiple tiles with overlapping areas. The spatial extent of each tile in the 1/9 arc-second and 1-m DEMs were 15 min and 1000 m, respectively. In terms of quality, the 1-m DEMs meet the quality level II requirements for vertical accuracy, whereas the 1/9 arc-second DEMs meet the quality level III. The specifications for each type of DEM, including the aggregate nominal pulse density (ANPD), aggregate nominal pulse spacing (ANPS), root mean square error in the vertical direction (RMSE_z), vegetated vertical accuracy (VVA), and non-vegetated vertical accuracy (NVA) are given in

Table 1. DEM specifications [52].

Quality Level	ANPD (pls/m2)	ANPS (m)	RMSEz (Non-Vegetated) (m)	NVA at the 95% Confidence Level (m)	VVA at the 95th Percentile (m)
II	≥2.0	≤0.71	≤0.100	≤0.196	≤0.3
III	≥0.5	≤1.41	≤0.200	≤0.392	≤0.6

[47,51]. This table defines the minimum requirements for accuracy and density of returns. The RMSE_z is defined as

$${\mathrm{RMSE}}_{\mathrm{z}}=\sqrt{{{\displaystyle \sum }}^{\text{}}\frac{{\left(z_n-z_n^{\prime }\right)}^2}{N}},$$

(3)

where $n$ is the identification number of each check point ranging from 1 to $N$, $z_n$ is the $n\mathrm{th}$ set of $z$ elevation values, $z_n^{\prime }$ is the corresponding check points, and $N$ is the total number of check points.

4.2. DEM Treatments

Seven 1/9 arc-second and twenty 1-m DEM tiles were merged into a single raster with continuous coverage, i.e., seamless. The average elevation value of the adjacent tiles was used to create a midpoint for each overlapping cell to determine the output cell value in the seamless DEM. Additionally, each DEM had a corresponding grayscale colormap for which the minimum and maximum values were proportionate to the minimum and maximum elevations of the area that it covered (Figure 3a). This was considered when merging the raster DEMs so that the resultant colormap could match the minimum and maximum elevations of the new raster to create a seamless DEM raster (Figure 3b). Next, the two DEM datasets needed to be processed and adjusted in multiple steps to prepare them for creating a DoD. The following paragraphs are organized to detail these steps.

The first preprocessing step was to address the discrepancy between the horizontal coordinate systems of the 1/9 arc-second DEM and the 1-m DEM. The former uses a geographical coordinate system based on a three-dimensional ellipsoidal surface; thus, angular measurements for the geodetic longitude and latitude are used to specify the location, i.e., $\left(\phi ,\lambda \right)$. The latter adopts a planar projected coordinate system and uses two-dimensional coordinates and linear distances as units, i.e., (x,y). The coordinate system of the 1/9 arc-second DEM raster was altered to match that of the 1-m DEM. For this purpose, the Universal Transverse Mercator (UTM) projection [53] was used to transform the geodetic locations of the 1/9 arc-second DEM to the planar system.

After projection, the second preprocessing step was to adjust the cell resolution so that both raster datasets had the same cell size. To achieve this goal, the resolution of the 1-m DEM raster was altered to match the coarser resolution of the 1/9 arc-second DEM (approximately 3 m). The resampled raster mapped each output cell to the input cells and assigned a value to the output based on the location of nearby inputs. Since the center of the output cells did not align with the inputs, it was necessary to use a resampling technique to rectify them. Different techniques could be employed to determine the value of the output cells, such as nearest neighbor, cubic convolution, and bilinear resampling. This study utilized the value of the four nearest input cells and bilinear resampling to obtain the weighted average of the input values based on the distance between inputs and the center of the output cells. This technique is most suitable for continuous surfaces, such as elevation, because it considers the location and distance between the known points and the unknown/output cells.

Next, the resultant DEM grid cells from the 1/9 arc-second dataset were plotted and compared with those of the 1-m DEM (Figure 4). This figure highlights the misalignment between the grids of the DEM data. This figure was constructed to show that although the cell resolutions in the $x$ and $y$ directions were the same for the two DEMs, the grids in Figure 4a,b did not align. The ArcGIS Pro software (Version 2.6) [54] was used in this study, and it was found that it does not consider this misalignment when comparing the two DEMs. Although the grid points were more than 1 m off in both directions (1.0671 m in the x-direction and 1.284 m in the y-direction), the software subtracted the elevations of the nearest points in the two raster layers to obtain the DoD. The third preprocessing step addressed this issue. In this step, the 1-m DEM was resampled using the 1/9 arc-second DEM as the snap raster. In this process, the bottom left corner of the 1-m DEM was stretched outward to the nearest cell corner of the 1/9 arc-second DEM raster. The same process was repeated for the upper right corners of the two raster layers, and bilinear resampling was used to obtain elevation values for the modified cells.

Following these preprocessing steps, the DEM datasets were prepared for comparison. A spatial analyst tool, i.e., a raster calculator in ArcGIS Pro, was utilized to subtract the 1/9 arc-second DEM from the 1-m DEM. This yielded a direct estimate of the vertical displacement during the considered period, i.e., the DoD shown in Figure 5a. This figure shows the resultant DoD layer formed from 16,352 rows and columns with a square cell size of 3.0587 m. This figure illustrates uplift and subsidence with green and red, respectively. To distinguish real surface deformation from noise in this map, an LoD threshold was used. Several thresholds were tested in a trial-and-error approach to obtain the optimum value for the LoD, starting with the propagated error based on the vegetated vertical accuracy of the DEM datasets, i.e., 0.7 m. A suitable threshold value should be large enough to minimize the noise in the new layer but not too large to omit essential landscape changes in the study area. The noise is most evident in areas with dense vegetation due to the extended time interval between the acquisition dates of the two DEMs and the difference in acquisition months relating to leaf-on and leaf-off conditions. Despite this, the trial-and-error procedure is quite fast, and optimum results can be obtained in a timely manner. It was found that the 1.5-m threshold was the optimum value to strike a balance between the above two considerations for landslide detection purposes of this study. The LoD threshold of 1.5 m was applied to Figure 5a, and the results were stored in a new raster in Figure 5b. This figure uses three colors to mark each pixel (cell) within the study area: green cells mark upward vertical displacements above the LoD threshold, red shows downward vertical displacements above the LoD threshold, and grey denote areas with vertical displacements below the LoD threshold. It is important to note that because the cell size for both DEMs is approximately 3 m, the smallest marked area is a square with a side length of 3 m.

4.3. DoD Analysis

The DoD layer represents the changes in elevation between the two DEMs with acceptable accuracy. However, not all elevation values are accurate, and RMSEz does not include all the uncertainties associated with DEMs. As a helpful indicator for accuracy, RMSEz is a measure to quantify the DEM correspondence to the source data, whereas factors such as image distortions, clouds, and steep terrains could also affect the quality of the DEM data. Furthermore, not all elevation changes were produced by landslides. Various other factors could cause elevation change, including but not limited to construction, changes in river basins, erosion, and water body treatments in DEMs. Therefore, marked areas in the filtered DoD layer, i.e., the DoD layer filtered via an LoD threshold, could be regarded as landslide candidates, and further assessment was needed to classify them into landslides (true positives) and nonlandslides (false positives). Two types of auxiliary maps were generated for this purpose: multidirectional Hillshade maps and customized vector basemaps.

Multidirectional Hillshade presents a 3D representation of the terrain with a grayscale color ramp. This technique combines six different light sources to visualize and shade the terrain. Though being qualitative—the elevation values are not displayed—it is a valuable tool for displaying more details and enabling a closer look at the terrain surface. Therefore, two multidirectional Hillshade maps were generated from the multi-temporal DEM data. These newly generated layers were used alongside the filtered DoD layer to assess the marked areas. Multidirectional Hillshades were used to visually identify patterns of contiguous elevation change in the study area and find areas with a downward movement in the upslope followed by an upward movement in the downslope. Such patterns helped distinguish landslides from other vertical displacements.

The customized vector basemaps were also important. Evaluating these maps confirmed our assumption regarding the DoD layer: many of the marked areas belonged to rivers, lakes, and the construction of new buildings and roads. For example, the elevation values for many water bodies differed between the two DEMs. Figure 6 was constructed to emphasize this point. In this figure, the blue rectangles mark some of the areas in which elevation changes were caused by the difference in elevation values of lakes and water bodies. Moreover, the customized basemaps helped distinguish real elevation changes from the ones caused by artifacts, such as remnants of buildings in the bare-earth DEM data.

It was found that the primary causes of the noises in the DoD layer (false positives) were water bodies and changes due to the construction of buildings and roads. Moreover, it is known that the risk of landslides increases near roads and rivers. Thus, another set of auxiliary maps was generated to mark roads, water bodies, and buildings. For this purpose, Vector Tile Style Editor [55], an online tool for customizing vector tile layers, was used. This tool enables editing and customizing a basemap from ArcGIS Living Atlas of the World [56] that contains many basemaps. In this study, a classic ESRI topographic map named World Topographic Map [57] was customized to only show the desired areas. The customized map uses three colors to denote water bodies, buildings, and roads in the study area. Figure 7 demonstrates an example within the study area for the customized map. In this map, water bodies, buildings, and roads are marked blue, purple, and black, respectively. All the other features in this vector basemap are transparent to aid visualization. In the next section, this map will be overlaid on the filtered DoD layer to demonstrate the geographical location of these features relative to the marked areas.

5. Results

5.1. Verification Using Landslide Inventories

To study the landslide locations, the filtered DoD and Hillshade layers were overlaid by the customized basemap. In addition, the grey color in the DoD layer was set to transparent to help visualization. This means that the vertical displacements below the LoD threshold were filtered out. Next, two landslide inventories and their associated reports for the events were used to assess the recorded landslide locations. The employed landslide inventories were the U.S. Landslide Inventory by USGS [58] and the Global Landslide Catalog (GLC) by NASA. The USGS landslide inventory is an interactive map compiled from different local, state, and federal organizations. This inventory includes links that provide access to the original inventory files for the associated reports and news articles [59]. The GLC is a database that was created to identify rainfall-triggered landslides in the world. This inventory uses media reports, scientific reports, disaster databases, and other sources to compile its dataset [60,61]. As is often the case in landslide inventories, the exact locations were unknown for many landslide records. Instead, the location of landslides was only known to be within a specific range around the marked location. The newly developed method detected these landslides and marked them in the filtered DoD layer. The marked landslide locations were then compared with the associated reports and articles in the news to verify their accuracy.

One example is the 2017 landslide(s) near Fairfax, California (Figure 8). In this figure, the blue cells represent water bodies, the purple cells show buildings, the black cells denote roads, the green cells denote areas with uplift, the red cells represent subsidence, the red dot marks the location of a landslide in the inventories, and the dotted blue boxes mark the outer boundaries of areas with uplift and subsidence. The first two images on the left belong to a multidirectional Hillshade created from the 1/9 arc-second DEMs in 2010 with (Figure 8a) and without (Figure 8b) the filtered DoD layer. Figure 8c is another Hillshade derived from the 1-m DEMs and depicts the same area nine years later in 2019. Figure 8a was intended to show the landslide location in the inventories and the detected changes in the DoD layer. The following two subfigures were constructed to compare pre-event (Figure 8b) and post-event (Figure 8c) Hillshades to evaluate the accuracy of the detections. Although DoDs are not present in these two figures, the blue dotted boxes mark the approximate location of elevation changes. The location of this landslide was known to be within 1 km from the labeled location in the inventories, i.e., the red dot. The DoD layer did not detect significant changes above the LoD threshold in the immediate vicinity of this area. Further investigations of the associated news articles and reports revealed that the red dot represents multiple landslides triggered by a winter storm in that neighborhood, which concurs with results obtained with our filtered DoD layer. Analyzing the DoD layer revealed several locations marked with elevation changes (Figure 8).

In fact, many landslide reports suggested the occurrences of several other landslides in the area were caused by the same triggering factors of the recorded event, i.e., rainfall and winter storms. However, the responsible agencies could not gather enough information and locate them, so they were omitted from the inventory. The new method can detect these landslides and mark their locations. Another example of such events is presented in Figure 9, which uses the same legends and sequence presented in the previous figure. As shown in this figure, there are many changes in the landscape of this region, with the vertical displacement values reaching the maximum of 4 m for the marked areas. Despite this substantial elevation change, the inventories did not record the marked landslides. This figure also demonstrates elevation changes due to construction: a building marked with both purple and red markers indicates elevation changes due to recent construction.

Figure 10 is constructed to highlight a case in which not only the exact location of the landslide was not known, but also only one single landslide from a series of occurred landslides was reported. The location of this landslide was known within 1 km, but the reports indicated that other landslides were reported by firefighters but not included in the inventories. Although several landslides are visible in this figure (red polygons), other elevation changes were also obtained nearby this area that supported the evidence found in the reports. Figure 11 is constructed to demonstrate another advantage of the developed method. That is, this method can detect landslides in unpopulated or sparsely populated areas, which have been paid less attention to when documenting landslides and compiling inventories. The proposed method identified many areas with an apparent change in the elevation, for which no landslide event was found in the inventories. Figure 11 is one example of such an area. As can be seen in this figure, there is a discernable landslide in this region. This landslide was not reported or documented.

5.2. Validation against Google Earth Images

High-resolution Google Earth images were used to verify the detection accuracy of the new method. For this purpose, marked areas in the filtered DoD layer were converted to polygon features and exported into a geodatabase feature class. Next, the feature polygons were imported into Google Earth Pro in two sets to show upward and downward vertical displacements marked via the filtered DoD layer. These layers were combined with the historical imageries of Google Earth Pro to view the terrains before and after the DEM data were acquired. Several factors were considered in an attempt to choose the most appropriate dates. These factors include the percentage of cloud-covered areas, distortions, and the quality of the images. Moreover, it was imperative to find the closest image acquisition months for pre-event and post-event images to minimize the effect of seasonal vegetation growth and atmospheric variations. The most suitable date for pre-event images belonged to October 2009. September 2018 and June 2019 were the best available dates for post-event images. The filtered DoD layer was then overlaid on the Google Earth images to visually assess the areas marked with elevation change. It is important to note that since Google Earth images do not provide bare earth elevations, this validation method was only performed for areas without dense canopy and vegetation. In the following, some of these cases are presented to verify and discuss the results of the proposed method.

Figure 12 is one of these cases showing an apparent change in the topography of the area. This figure uses yellow and red polygons to show upward and downward vertical movements, respectively. Although there is evidence of an exposed landslide scar that had happened prior to 2009 in the middle of Figure 12a, new crown cracks and new minor scarps are visible in the 2019 image as well (Figure 12c). Another noteworthy change is the altered shape of the old scarps and lines that were left from the previous landslide event. This indicates an active landslide that had been evolving between the two acquisition dates. Additionally, patterns of downward vertical movement in the upslope and upward movement in the downslope (landslide deposits) are visible in this figure. This observation implies that earth mass from the uphill moved downward and accumulated at lower locations.

Another example is presented in Figure 13. For the marked area on the left, there are traces of a scarp and a crack in Figure 13c. This subfigure also shows a post-event image of the landslide on the right. This landslide is easily detectable from the cracks parallel to the slope and noticeable scars in the 2019 Google Earth image. Figure 14 shows another area with patterns of accumulation of debris in the downslope and subsidence in the upslope, indicating landslide activities. A closer assessment reveals landslide scars that affected the patchy vegetation in the area in Figure 13c. Besides, signs of debris flow are visible from the top right of the image to the bottom left.

The following two figures are constructed to show examples of the destruction and removal of vegetation cover due to landslides, resulting in landslide scars. Figure 15 is presented to emphasize the effectiveness of the proposed method for landslide detection in vegetated areas. Moreover, the landslide in Figure 16 is identified similarly due to the scar left by the slope failure. However, the side length of the depicted area in Figure 16 is approximately 6 m, which demonstrates the high potential of the new method in identifying small landslides.

6. Discussion

This study employed map algebra operations to analyze high-resolution multi-temporal DEM data and to extract information regarding surface deformation and vertical displacements. Several preprocessing steps were followed to prepare the high-resolution DEMs for use, including map projections and resampling techniques. The results were verified against landslide inventories and Google Earth images. The new method can find an accurate location for all the detected landslides. This is a beneficial feature that can help improve the existing landslide reports and inventories that offer mere approximations regarding landslide locations.

Moreover, the proposed method minimizes human intervention: most preprocessing steps and procedures can be implemented via computer programs. This obviates the need for a human expert interpreter and provides new possibilities for landslide monitoring and the preliminary stage of landslide detection projects. This study accomplished the above tasks without introducing any cost for DEM acquisition by using publicly available DEM datasets of 3DEP. It is also worthwhile to mention that this performance was achieved by comparing 1/9 arc-second and 1-m DEMs. Considering the recent developments in remote sensing and GIS data, this method could be further enhanced by incorporating higher-resolution DEM data and the continuous feeding of temporal data into the model, i.e., adding DEM data from other times. It is also important to note that the accuracy of this method is only as good as the data built into it. Therefore, the use of LiDAR data in this study makes it prone to errors in areas with dense vegetation, heavy rain, and high altitudes.

The proposed method enables us to differentiate between the erosion caused by water bodies and landslides by supplementing the filtered DoD layer with customized basemaps, which is a noteworthy improvement over the existing methods. It was observed that many false positives, i.e., cases where elevation changes are falsely identified as landslides, belong to water bodies. These false positives indicate that the proposed method is unable to differentiate landslides from erosions caused by water bodies. Eliminating these erosions was especially important since they have similar patterns of erosion and deposition to landslides [32] and can be easily misclassified as landslides. However, erosion can also alter the stress distribution in the soil and disturb the static equilibrium that causes a landslide. Therefore, it is important not to label all elevation changes adjacent to water bodies and riverbanks as erosions. Site reconnaissance is perhaps the most suitable way to distinguish between erosion and landslides, but this study uses a different approach to minimize these misclassifications. A high LoD threshold ensures that most of the erosions are eliminated from the dataset. Additionally, landslides that are as small as one pixel and fall on the riverbanks are eliminated from the dataset.

Moreover, the performance of the new method was good, even though many geomorphic features were obscured by dense vegetation cover, which negatively impacted the accuracy of the DEM data due to limitations in ground detection. The dense vegetation also limited our ability to use Google Earth images to verify the detected landslides. For these areas, multidirectional Hillshade maps were generated and used to study ground deformations. Additionally, when considering DEMs that are nine years apart, revegetation can adversely affect the accuracy of the DEM data and our ability to identify elevation change. That is one of the reasons we adopted a relatively high LoD threshold, even though the low RMSE of the high-resolution DEMs suggested that the proposed method can identify smaller elevation changes. Therefore, better results are possible by acquiring the DEM data annually and preferably right after the rainy season in the study area because that is when most landslides occur and before most of the leaves appear. Annual DEM acquisition is also beneficial since the proposed method is limited by data availability. Data availability is most crucial for pre-event DEMs since these data are lost, and acquisition is impossible after the landslide event. Additionally, annual data are beneficial for detecting the rate of change and other kinematic landslide features for each landslide which is an important factor used to characterize landslides.

An area that the proposed method struggles with is identifying slow-moving landslides that creep at low rates. This is especially important to consider when applying this method to annual DEMs or when the rate of change is as low as a few millimeters or centimeters per year. In this study, the time interval is nine years which renders this issue less significant, but this is worthwhile to consider when interpreting the results. Slope-parallel and translational landslides are other examples of landslides that cannot be identified by the proposed method. These types of landslides do not change the elevation significantly, and most of the displacements are in the horizontal direction and along a planar surface of weakness.

7. Conclusions

A new method for landslide detection and mapping was proposed and tested in this study. The 3DEP products by USGS were utilized to obtain multi-temporal high-resolution DEMs and detect changes in the landscape for a study area in Northern California. Several preprocessing steps were followed to prepare the DEM datasets for use and ensure their spatial alignment. A DoD layer filtered via an appropriate LoD threshold was employed to indicate local elevation changes between the DEM datasets. It was found that the filtered DoD layer, vector basemaps, and multidirectional Hillshade layers can supplement one another to obtain a more realistic and detailed understanding of the changes in the landscape and the geomorphological evolution of landslides. Using this method, subtle ground displacements caused by small and medium-sized landslides, which were not reported previously, were successfully detected. These vertical displacements were verified with landslide inventories and Google Earth images. One of the contributions of this method is that the filtered DoD layer can be combined with DEM datasets to provide pre-event and post-event DEM data for each detected landslide event. DEM data are especially important for landslides in sparsely populated areas because they have been paid less attention due to the lack of fatalities and damage to assets. These data can be employed to extract knowledge from past landslides and study their mechanisms, features, triggering factors, and myriads of other applications.