Rill Erosion and Drainage Development in Post-Landslide Settings Using UAV–LiDAR Data

Abstract

Accurate microtopography data are an important input for characterizing small-scale rill erosion and its progression following disturbances. UAV–LiDAR systems are increasingly accessible and have successfully been used to measure microtopography data for several applications. Yet, the use of UAV–LiDAR systems for rill erosion studies in post-landslide landscapes have not been well investigated. Therefore, the purpose of this study was to implement and evaluate a UAV–LiDAR-based workflow to capture the microtopography of a post-landslide landscape, and by doing so, to help to determine best practices for UAV–LiDAR-based rill analysis. A commercial UAV–LiDAR system was used to map three post-landslide slopes and generate digital elevation models with a 1 cm-per-pixel ground resolution. Using data captured over multiple years, temporal rill development was assessed by comparing rill cross-sections and calculating changes to rill density and erosion volume. A flow-accumulation algorithm was adopted to automatically extract the rill network. We found that a flow accumulation algorithm with a threshold value of 5000 detected the rill network with overall accuracies of >88% and F1-scores of >93%. Vertical cross-sections of individual rills revealed an increase in the depth and width of rills over a one-year period. This study demonstrates that a commercial UAV–LiDAR system can effectively describe microtopography in a post-landslide landscape and facilitate analysis of small-scale rill characteristics and the progression of rill erosion.

1. Introduction

Soil erosion is an important process that affects the topography of landscapes [1], leads to loss of fertile soils [2], and in forests can affect the forest structure and aboveground biomass [3]. Soil erosion is often initiated or exasperated after disturbance events such as landslides or the removal of vegetation during construction, agriculture, or forestry activities. During the early stages of erosion, soil loss occurs primarily through the development of rills, which are small watercourses with steep slopes and widths. The dimensions of rill features differ between studies. For example, Li et al. [4] define rills as landforms between 2 and 20 cm, but rill widths exceeding 20 cm are reported by Eltner et al. [5], and widths exceeding 1 m are reported by Kou et al. [6]. Significant research attention has recently been given to understanding the formation and characteristics of rills, while methodological studies have investigated different measuring techniques and algorithms for the extraction and description of rills, a review of which is presented by Out et al. [7]. Improvements in the understanding of rill processes are needed to refine erosion models, such as EUROSEM [8] or GUEST [9], which are used to simulate erosion in different scenarios, and inform proactive mitigation strategies.

Due to the fine scale at which rill erosion occurs, studies use specialized measurement equipment that can measure these small landforms accurately. These include profilometers [10], terrestrial laser scanning (TLS) devices [11], and photogrammetry-based techniques. For photogrammetry, collections of aerial images collected under laboratory conditions [12] or from drone surveys [5] are processed to generate high-resolution digital elevation models (DEMs) [13]. These DEMs offer microtopography data at the spatial scale required for rill erosion studies. The introduction of new measurement approaches and higher resolution datasets have also lead to an increase in available methods for detecting rills automatically or semi-automatically [5,11,14]. In light of the fast development of measurement and detection workflows in rill research, the standardization of protocols, description of flight parameters, and data capture methods is required for accurate observation [7].

Unmanned aerial vehicles (UAVs) have advanced environmental monitoring, providing a versatile supplement to traditional remote sensing tools (i.e., satellite and airborne sensing) [15]. UAVs can carry different types of sensors, including LiDAR, hyperspectral imagers, and high-resolution cameras, which facilitate comprehensive 3D terrain modelling and vegetation analysis [16,17]. UAVs are often used in areas that are too small for traditional aerial surveys but too large for efficient ground-based studies, providing high-resolution orthorectified imagery of target areas and digital elevation models that can capture the microtopography of disturbed landscapes and slopes [18]. UAV–LiDAR systems, which integrate drones with light detection and ranging technology, is currently emerging as an alternative to UAV photogrammetry as an approach for capturing high-resolution (cm scale) topographic data [19]. Despite the high measurement accuracy generally associated with LiDAR systems, cases of integration on lightweight UAVs have historically been seldom, due to LiDAR sensors’ large size and weight, high power requirements, and prohibitive costs [20,21]. Their use was thus primarily constrained to satellite or aircraft (ALS), and employed by government or large contractors for high-accuracy topographic surveys over relatively large areas [22]. However, the miniaturization of LiDAR technology is leading to expanded use as a part of drone systems, allowing for collection of accurate topographical data [23,24].

UAV photogrammetry and UAV–LiDAR as well as the integration between both techniques have been successfully applied to measure micro-topography and its changes—for instance, in detecting landslides, gully erosion, and rill erosion [6,25,26,27,28,29]. Kou et al. [6] used UAV photogrammetry-based digital surface models to quantify rill erosion in the Loess Plateau of China. Several studies considered airborne LiDAR (ALS)-derived DEMs as topography data sources for geomorphology studies. However, ALS-derived DEMs are generally at a lower 1 m-per-pixel ground resolution, which precludes their direct use in rill erosion studies [29,30,31]. Previous work using UAV–LiDAR for small-scale erosion assessments did not consider rills specifically, but was rather interested in measuring volumetric changes caused by erosion using DEM differencing [27]. Their study also considered a 0.1 m-per-pixel DEM resolution, which may be too coarse still for assessing certain rill parameters. To the authors’ knowledge, the use of UAV–LiDAR systems applied to rill erosion studies on a post-landslide landscape has not been specifically evaluated, despite the potential advantages of such a system.

UAV–LiDAR-based measurements of microtopography in rill landscapes could have advantages over UAV photogrammetry-based measurement techniques, since it (1) directly measures topography, as opposed to using photogrammetric techniques to infer 3D information from 2D images, and (2) may require fewer field survey and computing resources to generate DEMs than photogrammetry approaches. In addition, (3) shadows cast into rills may cause errors in the 3D reconstruction by photogrammetry, leading to erroneous rill width or depth estimates [32]. Such errors may be avoided by LiDAR, which uses active remote sensing (laser) to measure topography, irrespective of light or shadow. However, due to the lack of existing UAV–LiDAR rill erosion studies, the specific considerations for this application, such as appropriate flight parameters, suitable DEM generation algorithms, achieved accuracies, necessary error removal steps, and the scale of rills that can effectively be identified from the workflow need still to be described.

The objective of this study is to demonstrate the application of an unmanned aerial vehicle light detection and ranging (UAV–LiDAR) system to monitor rill erosion by mapping microtopography on post-landslide slopes. The study adapts a flow accumulation-based algorithm commonly used in basin-scale hydrology studies in order to automate rill network extraction. Temporal progression of the rills is also assessed by calculating changes in microtopography data collected in different years. To quantify rills and their progression, cross-profiles, rill density, and volumetric changes are shown for the three slopes. Difficulties arising from the required fine resolution of the data, such as artifacts in multi-temporal change analysis and small vegetation regrowth, are also highlighted and discussed. This case study is expected to be used as a reference for other studies that wish to apply a UAV–LiDAR methodology to capture topography data for rill erosion studies, enabling future standardization of workflows and promoting the adoption of UAV–LiDAR in rill analysis studies.

2. Materials and Methods

2.1. Study Area

The 2018 Hokkaido Eastern Iburi earthquake, registering a magnitude of Mw 6.6, had profound geological and environmental impacts on the northern region of Japan, particularly in Hokkaido. This significant seismic event triggered an extensive series of over 7000 landslides across an area of approximately 466 square kilometers, predominantly affecting the towns of Atsuma and Abira. These landslides resulted in substantial alterations to the landscape, posing serious challenges to the affected communities and the environment [33,34,35,36].

This study was conducted in the Atsuma region (42°46′31.18″ N, 142°0′45.28″ E), located in the southern part of Hokkaido (Figure 1). Post-earthquake, numerous landslide sites became evident, with notable rill erosion on these exposed slopes. Many of the landslides were in areas with minimal human activity, allowing for the natural progression of rill erosion to be observed.

Most of the landslides were shallow and occurred in the air-fall pumice deposits released by the Tarumai volcano. Air-fall pumice is susceptible to landslides during seismic activity due to its loose and granular nature. Additionally, several deep-seated landslides were noted in the southeastern areas, indicating variations in the geological substrate and varying responses to seismic activity [37,38]. The distribution of these landslides typically occurred on gentle hills at elevations ranging from 100 to 250 m, with vertical heights of landslides measuring 15 to 45 m. The primary alignment of the landslides followed a south–southeast to north–northwest direction, reflecting both the underlying geological structures and the propagation direction of seismic waves during the earthquake [39].

2.2. UAV System

This research employs the DJI Matrice 300 RTK drone, integrated with the Zenmuse L1 LiDAR sensor. This UAV–LiDAR system offers high precision acquisition with a range of up to 450 m, an effective point rate of 240,000 points per second, and a spatial accuracy of around 3–5 cm [40,41]. Figure 2 shows the implementation of drone survey in the field to obtain aerial photographs and LiDAR data.

2.3. Field Survey & Flight Parameters

To measure the temporal development of rills in the study area, fieldwork was conducted 1 year apart, on 2 August 2023 and 2 September 2024. During the field surveys, the UAV–LIDAR was programmed to capture the different slopes with the flight parameters described in Table 1.

The initial survey flights were conducted at 120 m above ground level to generate a preliminary digital elevation model (DEM). For subsequent flights, the UAV was operated in terrain follow mode, which automatically adjusted the flying altitude to maintain a constant 30 m distance from the ground surface. This approach was deemed necessary to ensure consistent LiDAR point cloud density over the sloped terrain of the landslide.

2.4. Digital Elevation Model Generation

Following the acquisition of raw data from UAV–LiDAR, we used the DJI Terra software (version 4.2.5) to generate 3D point cloud LAS files. Furthermore, we processed the RGB data in DJI Terra alongside the point cloud data. The processing of RGB datasets produced orthomosaics for each acquisition date (2023 and 2024) in the study area. We used the orthomosaics to validate the results of rill detection through visual interpretation.

The triangulation function from the LidR package in R was used to transform the 3D point cloud into 2D digital elevation model (DEM) raster images with a resolution of 1 cm. To reduce local irregularities and noise in the elevation values of the DEM, we imported the data into MatLab and applied a low-pass mean filter to the DEMs. The kernel size of the mean filter affects the degree of smoothing in the DEM. However, there is a risk that smoothing may remove actual aspects of the microtopography. To avoid this issue, a range of kernels was considered, and their effect on the DEM was compared to determine smoothing parameters that can remove small local irregularities in the DEM.

Georeferencing errors caused elevation offsets of around 8 m between DEMs produced for different years. To address this, calibration reference points were selected at locations that were considered least affected by erosion and would thus have the same elevation in the multitemporal DEMs. Locations near exposed rock sections and points towards the upper slopes were prioritized while avoiding areas covered by vegetation (see Figure 3). These specific points were used as reference points for vertical alignment and normalization between two DEMs in a particular site, following the approach described in Figure 3. The corrected DEMs were then finalized and used as the primary dataset for further analyses. In addition to the vertical corrections described above, small horizontal shifts between the elevation models produced for different years were also corrected during post-processing. For this, we used the georeferencing toolbox in ArcGIS Pro to overlay the DEMs and drone orthophotos of each slope and visually select common reference points between the rasters. These were used to update the georeferencing of the rasters and thus improve the horizontal spatial alignment of the data.

2.5. Rill Extraction & Change Analysis

The workflow of this study is illustrated in Figure 4. To identify the rill network on the landslide slopes, we used a flow accumulation and stream delineation based approach traditionally used in hydrology studies [42]. This choice was made because the rill network resembles a larger watershed and stream network, although at a smaller spatial scale. Stream delineation is sensitive to the scale of the features and the resolution of input DEM data [43]. Therefore, it was necessary to adjust the flow accumulation threshold by trial and error until the generated rill network accurately resembled the visually observed network. For the 1 cm resolution DEM data, a flow accumulation threshold of 5000 returned the most complete rill network. Flow accumulation is quantified as the number of cells upstream of the current cell. The threshold denotes the accumulation value above which a stream is initiated. The workflow was executed using the hydrology tools available in the spatial analyst toolbox in ArcGIS Pro. The specific steps comprised the generation of fill, flow direction, and flow accumulation, followed by applying the flow accumulation threshold to generate a stream network. We produced a vector-based rill map by converting the raster image to a vector polyline using this threshold. This result was used to calculate the rill density in each slope for the two study years and analyze the changes. Rill density is defined as the total length of rills (m) per unit area (m²) [7].

To assess the accuracy of the extracted rill network, 60 validation sample plots of 27.8 m² each were randomly created across the three considered slopes and two years of data. A ground truth rill map was manually digitized for each sample plot and overlaid with the automatically extracted results. We counted the number of true positives (TPs), false negatives (FNs), and false positives (FPs), representing correctly detected rills, missed or undetected rills, and falsely detected rills, respectively. Two main accuracy metrics were calculated, consisting of the F1-score and overall accuracy. The F1-score represents the weighted average of precision and recall (Equations (1)–(3)), while overall accuracy quantifies the proportion of correctly detected rills among all rill-related predictions. Recall measures the proportion of actual rills that were correctly detected, whereas precision measures the proportion of detected rills that are actual rills.

(1)$\mathrm{Precision}={\displaystyle \frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}}}$

(2)$\mathrm{Recall}={\displaystyle \frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}}$

(3)$\mathrm{F}1\text{-score}=2\times {\displaystyle \frac{\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\times \mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}{\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}+\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}}$

Cross-section profiles of the main two rills on each slope was taken at three points along the length of the rill (near the origin of the rill, halfway down, and near the base of the slope). Profiles for both years were extracted and overlaid to visually assess the pattern of rill progression at the different section of the rill. Finally, to quantify the overall erosion pattern, temporal DEM difference images were generated for each slope and visually interpreted.

3. Results

3.1. DEM Filtering Effect

The cross-section profile in Figure 5 illustrates the effects of filtering kernel resolution on the smoothness of topography features in the DEM data. By visual inspection, a filter with a 0.2 m kernel size produced a smooth DEM that captured the main terrain features at the sub 1 m scale, but helped to reduce local irregularities at the 1–5 cm scale (red line in Figure 5), whereas larger kernel sizes negatively affected the rill pattern.

3.2. Rill Extraction Result and Validation

To quantify automatic rill extraction accuracy, 60 validation sample plots of 27.8 m² each were randomly created across the three considered slopes and two years of data. The ground truth rill map was manually digitized for each sample plot, as shown in Figure 6. In both years, high recall and precision were obtained (recall, precision > 90%, Table 2). Our method achieved an F1-score and overall accuracy of 93.70% and 88.14% in 2023, and of 94.28% and 89.17% in 2024.

The rill map for a larger area on the slope can be visually assessed by overlaying it to the site’s DEM (Figure 7a) and the orthophoto (Figure 7c) generated from the RGB imagery captured by the drone.

3.3. Rill Progression

In all three study sites, there were only small changes in the rill density between 2023 and 2024 (in the order of 0.01–0.02 m/m², Table 3). However, when considering individual rills, clear temporal progression in the rill depth and width was observed (Figure 8, Figure 9 and Figure 10). The cross-profiles show that the depth changes in pre-existing rills were larger than in newly formed rills. This can be seen by the difference in the depth changes between the main and secondary rills in cross-sections A2 and A3 in Figure 8. This may be a result of more intense channeling of water, leading to a higher flow rate and consequent erosion in already existing rills. In Site 2 (Figure 9), a spatial shift in the cross-profiles can be seen, which suggests a georectification error in the DEMs, causing a spatial shift in the horizontal direction.

To further show the spatial pattern of erosion between the two years, the DEMs of 2023 and 2024 were subtracted from each other, as shown in Figure 11, Figure 12 and Figure 13. The resulting maps show that rill erosion was more rapid than inter-rill erosion, as shown by the well-defined green lines along the rills. But inter-rill erosion can also be observed, especially along the lower portions of the slopes (e.g., Figure 11, lower-left portion of the slope). This may be due to faster overland flow farther down the slope. The georectification error observed in the cross-profiles of slope 2 (Figure 9) are also visible in the DEM difference map (Figure 12), where there is large deposition (red color) reported on the east-facing part of the slope and erosion (green color) on the west-facing portions. This deposition and erosion pattern is most likely an artefact of the slight westward shift in the 2024 DEM, caused by the georectification mismatch between the DEMs. Overall, on all three slopes, the deepest erosion occurred at the base of the slope, as shown by the dark green patches at the bottom of the slopes in Figure 11, Figure 12 and Figure 13.

During the field investigations, we found that the slope at Site 3 had the best vegetation recovery of the three sites. As a result, the development of rill erosion was relatively slow, which is also apparent from the lack of green color (indicating soil loss) in Figure 13. On the contrary, most of the site shows an elevation increase in the DEM (red color), which reflects the growth of shrubs and vegetation observed during the field work.

4. Discussion

The presented study shows the application of a UAV–LIDAR system to generate high-resolution digital elevation models for use in a rill erosion study. The work departs from related studies that generated microtopography data from terrestrial laser scanning [11], drone photogrammetry [44,45], or aircraft-borne LIDAR [46]. The latter—aircraft-borne LIDAR—generally generate DEMs with a spatial resolution of around 1 m, which makes it appropriate for studying larger erosion-induced landforms, such as gullies [46,47]. However, such data products are too coarse for studying rill landforms (Figure 8) [7]. In comparison, using the UAV–LIDAR methodology proposed in this study, the rill features could be extracted and their changes over time measured. UAV–LIDAR is thus a feasible alternative for drone photogrammetry and TLS for generating high-resolution DEMs in rill landscapes. Additionally, it does not suffer from the issue of shadow cast into rills leading to erroneous DEM reconstruction, which may occur with drone photogrammetry [32].

A further contribution of this work was to clarify the necessary field surveying parameters and processing steps required for UAV–LIDAR rill studies. It was observed that a 30 m above-terrain flying height was suitable to generate a DEM with 1 cm spatial resolution. However, due to the uneven terrain of a post-landslide slope, it was necessary to use a “terrain-following” flight profile to ensure that the density of the captured point cloud remained consistent across the study site. The study also proposed the extraction of the rill network by using flow accumulation and stream network delineation algorithms traditionally used in larger-scale hydrology studies. With a flow accumulation threshold of 5000, the rill network could be extracted with a high accuracy rate (overall accuracies > 88% and F1-scores > 93%) (Table 2). This thus poses a promising alternative to other approaches that rely on manual identification or machine learning-based classification [25,48]. However, the stream delineation-based method did not automatically define the edges of the rill and could thus not be used to extract rill parameters such as rill width and depth semi-automatically, as was reported by Eltner et al. [5], who used an edge detection algorithm to semi-automatically extract rill parameters. In future work, the rill network extraction algorithm presented here can be extended to a two-step process where, following extraction of the rill network, the rill edges can be identified using edge detection algorithms such as the canny edge detector. A similar multi-stage procedure was proposed for delineating agricultural terraces from DEMs [49]. This algorithm proposed by Dai et al. [49] may be particularly suitable for adaption to rill edge detection from LiDAR-based DEMs, since the edge detection algorithm is applied to the hill shade image of the terrain and thus does not require RGB image data.

However, during the data collection and analysis, some challenges were noted. A substantial portion of our effort in data processing was devoted to data calibration, especially the height correction of the DEM, which could potentially have affected the accuracy of the results. The terrain of our study site had steep slopes, complicating the placement and stabilization of ground control points (GCPs). The UAV–LIDAR drone used in the study is equipped with a high-precision GNSS, which has been suggested to make GCPs unnecessary for our considered UAV–LiDAR system [50]. Yet, spatial misalignments in generated point clouds and DEMs of subsequent years necessitated manual georectification during post-processing. Despite manual georectification efforts, small spatial misalignments still led to visible artefacts in the rill progression result in site 2 (Figure 9 and Figure 12). Future studies should consider installing stable control points at the base of slopes to enhance the precision of data and results, even in cases where centimeter-accuracy GNSS systems are used in the UAV–LiDAR system. Even centimeter misalignments in the spatial rasters may compromise further analysis when working with multiple data sources or multi-temporal data, as shown in Figure 9.

Further, during the initial stages of handling 3D point clouds, we attempted denoising and point cloud classification to accurately differentiate low vegetation from the underlying soil. However, the outcomes were unsatisfactory. Issues with ground classification are most pronounced on slopes with significant vegetation regrowth, like the slope in Site 3 (Figure 13), which shows elevation increases in the later DEM, reflecting the vegetation regrowth. If vegetation growth is not accounted for, soil loss will be underestimated in temporal change analysis. Current common algorithms for point cloud classification are designed to distinguish between trees and ground points. In our case, the challenge was differentiating between ground and low-lying bush and shrubs, which have distinct cloud characteristics from trees, such as sparser density and smaller scale, leading to significant discrepancies from the expected results. Landslide-affected slopes often have these characteristics due to vegetation loss and partial regrowth by pioneer species. Therefore, further research should investigate ground point cloud classification algorithms that are more suitable for landslide-affected slopes. For example, TLS point cloud-based vegetation classification algorithms [51] may be adopted for use with UAV–LiDAR systems. Alternatively, a lower flight altitude may increase point density, allowing more pulses to penetrate through the vegetation and measure the topography of the underlying slope.

The present study contributes to better understanding of best practices and standardization of workflows for UAV–LiDAR-based rill analysis. Nevertheless, this study is restricted to a particular region of Hokkaido, Japan, and might not accurately represent the features of rill erosion and landslides in other areas. To confirm the method’s universality and flexibility, future studies should take into account comparable investigations conducted in regions with various geological origins and climatic conditions.

5. Conclusions

In this work, a UAV–LiDAR system was able to capture microtopography data on a post-landslide slope with sufficient accuracy to allow for the extraction of rills and the measurement of erosion. This approach thus offers an alternative to terrestrial laser scanning, manual field measurement, and UAV photogrammetry for measuring microtopography in rill erosion studies. We found that a terrain-following flight plan 30 m above ground level allowed the LiDAR to capture a point cloud with sufficient point density to generate a 1 cm-ground-resolution DEM, despite the undulating and sloped landslide terrain. However, on slopes with significant vegetation regrowth, a lower flight height would improve the differentiation of vegetation and the ground beneath. Further, the study shows that the high-precision RTK GPS embedded in the drone–LiDAR system did not realize sufficient georeferencing accuracy to directly allow for an overlay of elevation models captured at different times. For microtopography applications, it is thus critical to install persistent ground control points over the course of the study period, although this is difficult to implement in steep post-landslide slopes. The flight parameters, software workflows, and described limitations are intended to be a reference for follow-up UAV-LiDAR rill erosion studies or related applications that rely on microtopography data.

Author Contributions

Conceptualization: X.C., A.S.L., A.P.Y. and R.A. Software, X.C., A.S.L., D.A.U. and R.A. Formal analysis, X.C., A.S.L., A.P.Y. and R.A. Investigation, X.C., S.A. and M.A.H.B. Resources, R.A. and A.P.Y. Data curation, X.C., A.S.L. and A.P.Y. Writing—original draft preparation, X.C., A.S.L., A.P.Y. and R.A. Writing—review and editing, S.A., D.A.U. and M.A.H.B. Supervision, R.A. Project administration, R.A. Funding acquisition, R.A. and S.A. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is available upon request.

Acknowledgments

The authors would like to thank the Geospatial Information Authority (GSI) of Japan for providing the open topographic data. Chen Xinyu would like to thank DX Scholarship (JST SPRING—No: JPMJSP2119), Hokkaido University. The authors are grateful for the guidance and data access support provided by Yuichi Hayakawa.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Study area location in the Atsuma region. (a) Map showing the location of Hokkaido Prefecture (blue) within Japan (green), (b) satellite image of Hokkaido with a red polygon indicating the study region, and (c) the three research sites (Sites 1–3) outlined in red in the Atsuma region (basemap source from Google Earth).

Figure 2 (a) Data acquisition using the M300 RTK drone in the study site. (b) The Zenmuse L1 LiDAR sensor is mounted on the drone. (c) The drone flight plan from the remote controller.

Figure 3 DEM elevation calibration workflow.

Figure 4 Methodology flowchart.

Figure 5 Cross-section profiles of digital elevation models (DEMs) with different kernel sizes used for noise removal in DEM elevation values (DEMs with different filters applied distinguished by color of the cross-section line).

Figure 6 Example of validation of the extracted rill network, where red and green circles delineate sample plots, each with an area of 27.8 m². (a) Orthomosaic; (b) DEM; (c) result of rill detection; (d) reference of rill features; (e) distribution of validation plots.

Figure 7 Example of a rill map: (a,c) LiDAR-based DEM and orthomosaic; (b,d) the detected rill features.

Figure 8 The temporal rill profile line in Site 1.

Figure 9 The temporal rill profile line in Site 2.

Figure 10 The temporal rill profile line in Site 3.

Figure 11 The subtraction result (in meters) between the 2024 and 2023 data at Site 1.

Figure 12 The subtraction result (in meters) between the 2024 and 2023 data at Site 2.

Figure 13 The subtraction result (in meters) between the 2024 and 2023 data at Site 3.

Table 1

UAV-LIDAR flight parameters.

Parameter	Value
LIDAR sensor	Zenmuse L1
Altitude above terrain	30 m
Flying speed	3 m/s
LiDAR data overlap	40%
Optical image overlap	70%

Table 2

Detection performance of the rill extraction result.

Year	Detection Result	Count	Metrics (%)
2023	Detected rills (TPs)	104	Precision = 92.86
	Missed rills (FNs)	6	Recall = 94.55
	False positives (FPs)	8	F1-score = 93.70
	Total actual rills	110	Overall detection accuracy = 88.14
2024	Detected rills (TPs)	107	Precision = 98.17
	Missed rills (FNs)	11	Recall = 90.68
	False positives (FPs)	2	F1-score = 94.28
	Total actual rills	118	Overall detection accuracy = 89.17

Table 3

Rill density by study site and year.

Study Site	Rill Density (m/m2)
	2023	2024
1	0.50	0.52
2	0.54	0.53
3	0.53	0.52

References

1. Wang, J.; Howarth, J.D.; McClymont, E.L.; Densmore, A.L.; Fitzsimons, S.J.; Croissant, T.; Gröcke, D.R.; West, M.D.; Harvey, E.L.; Frith, N.V. . Long-Term Patterns of Hillslope Erosion by Earthquake-Induced Landslides Shape Mountain Landscapes. Sci. Adv.; 2020; 6, eaaz6446. [DOI: https://dx.doi.org/10.1126/sciadv.aaz6446] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32548261]

2. Lal, R. Soil Erosion Impact on Agronomic Productivity and Environment Quality. Crit. Rev. Plant Sci.; 1998; 17, pp. 319-464. [DOI: https://dx.doi.org/10.1080/07352689891304249]

3. Milodowski, D.T.; Mudd, S.M.; Mitchard, E.T.A. Erosion Rates as a Potential Bottom-up Control of Forest Structural Characteristics in the Sierra Nevada Mountains. Ecology; 2015; 96, pp. 31-38. [DOI: https://dx.doi.org/10.1890/14-0649.1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26236887]

4. Li, C.; Wang, M.; Liu, K.; Xie, J. Topographic Changes and Their Driving Factors after 2008 Wenchuan Earthquake. Geomorphology; 2018; 311, pp. 27-36. [DOI: https://dx.doi.org/10.1016/j.geomorph.2018.03.019]

5. Eltner, A.; Baumgart, P.; Maas, H.-G.; Faust, D. Multi-Temporal UAV Data for Automatic Measurement of Rill and Interrill Erosion on Loess Soil. Earth Surf. Process. Landf.; 2015; 40, pp. 741-755. [DOI: https://dx.doi.org/10.1002/esp.3673]

6. Kou, P.; Xu, Q.; Yunus, A.P.; Dong, X.; Pu, C.; Zhang, X.; Jin, Z. Micro-Topographic Assessment of Rill Morphology Highlights the Shortcomings of Current Protective Measures in Loess Landscapes. Sci. Total Environ.; 2020; 737, 139721. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2020.139721]

7. Ou, X.; Hu, Y.; Li, X.; Guo, S.; Liu, B. Advancements and Challenges in Rill Formation, Morphology, Measurement and Modeling. Catena; 2021; 196, 104932. [DOI: https://dx.doi.org/10.1016/j.catena.2020.104932]

8. Morgan, R.P.C.; Quinton, J.N.; Smith, R.E.; Govers, G.; Poesen, J.W.A.; Auerswald, K.; Chisci, G.; Torri, D.; Styczen, M.E. The European Soil Erosion Model (EUROSEM): A Dynamic Approach for Predicting Sediment Transport from Fields and Small Catchments. Earth Surf. Process. Landf.; 1998; 23, pp. 527-544. [DOI: https://dx.doi.org/10.1002/(SICI)1096-9837(199806)23:6<527::AID-ESP868>3.0.CO;2-5]

9. Mahmoodabadi, M.; Ghadiri, H.; Rose, C.; Yu, B.; Rafahi, H.; Rouhipour, H. Evaluation of GUEST and WEPP with a New Approach for the Determination of Sediment Transport Capacity. J. Hydrol.; 2014; 513, pp. 413-421. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2014.03.060]

10. Kimaro, D.N.; Poesen, J.; Msanya, B.M.; Deckers, J.A. Magnitude of Soil Erosion on the Northern Slope of the Uluguru Mountains, Tanzania: Interrill and Rill Erosion. Catena; 2008; 75, pp. 38-44. [DOI: https://dx.doi.org/10.1016/j.catena.2008.04.007]

11. Vinci, A.; Brigante, R.; Todisco, F.; Mannocchi, F.; Radicioni, F. Measuring Rill Erosion by Laser Scanning. CATENA; 2015; 124, pp. 97-108. [DOI: https://dx.doi.org/10.1016/j.catena.2014.09.003]

12. Qin, C.; Wells, R.R.; Momm, H.G.; Xu, X.; Wilson, G.V.; Zheng, F. Photogrammetric Analysis Tools for Channel Widening Quantification under Laboratory Conditions. Soil Tillage Res.; 2019; 191, pp. 306-316. [DOI: https://dx.doi.org/10.1016/j.still.2019.04.002]

13. Yunus, A.P.; Xinyu, C.; Catani, F.; Subramaniam, S.S.; Fan, X.; Jie, D.; Sajinkumar, K.S.; Gupta, A.; Avtar, R. Earthquake-induced soil landslides: Volume estimates and uncertainties with the existing scaling exponents. Sci. Rep.; 2023; 13, 8151. [DOI: https://dx.doi.org/10.1038/s41598-023-35088-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37208531]

14. Di Stefano, C.; Nicosia, A.; Palmeri, V.; Pampalone, V.; Ferro, V. Rill Flow Velocity and Resistance Law: A Review. Earth-Sci. Rev.; 2022; 231, 104092. [DOI: https://dx.doi.org/10.1016/j.earscirev.2022.104092]

15. Avtar, R.; Watanabe, T. Unmanned Aerial Vehicle: Applications in Agriculture and Environment; Springer International Publishing: Cham, Switzerland, 2020.

16. Elaksher, A.; Bhandari, S.; Carreon Limones, C.A.; Lauf, R. Potential of UAV Lidar Systems for Geospatial Mapping. Proceedings of the Lidar Remote Sensing for Environmental Monitoring 2017; San Diego, CA, USA, 30 August 2017; Singh, U.N. SPIE: Pune, India, 2017.

17. Manfreda, S.; McCabe, M.F.; Miller, P.E.; Lucas, R.; Pajuelo Madrigal, V.; Mallinis, G.; Ben Dor, E.; Helman, D.; Estes, L.; Ciraolo, G. . On the Use of Unmanned Aerial Systems for Environmental Monitoring. Remote Sens.; 2018; 10, 641. [DOI: https://dx.doi.org/10.3390/rs10040641]

18. Hackney, C.; Clayton, A.I. Unmanned Aerial Vehicles (UAVs) and Their Application in Geomorphic Mapping; British Society for Geomorphology: London, UK, 2015.

19. Pilarska, M.; Ostrowski, W.; Bakuła, K.; Górski, K.; Kurczyński, Z. The Potential of Light Laser Scanners Developed for Unmanned Aerial Vehicles—The Review and Accuracy. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2016; XLII-2-W2, pp. 87-95. [DOI: https://dx.doi.org/10.5194/isprs-archives-XLII-2-W2-87-2016]

20. Jeong, N.; Hwang, H.; Matson, E.T. Evaluation of Low-Cost LiDAR Sensor for Application in Indoor UAV Navigation. Proceedings of the 2018 IEEE Sensors Applications Symposium (SAS); Seoul, Republic of Korea, 12–14 March 2018; pp. 1-5.

21. Winsemius, H.; Mather, S.; Gayton, I.; Chazua, I. OpenDroneMap360, an Affordable DIY Open-Source Hardware and Software Workflow for 3D Point Clouds and Terrain Models. Proceedings of the Copernicus Meetings; Online, 19–30 April 2021.

22. Lu, Y.; Wang, Y.; Parikh, D.; Xin, Y.; Lu, G. Extending Single Beam Lidar To Full Resolution By Fusing with Single Image Depth Estimation. Proceedings of the 2020 25th International Conference on Pattern Recognition (ICPR); Milan, Italy, 10–15 January 2021; pp. 6343-6350.

23. Choi, S.-K.; Ramirez, R.A.; Kwon, T.-H. Acquisition of High-Resolution Topographic Information in Forest Environments Using Integrated UAV-LiDAR System: System Development and Field Demonstration. Heliyon; 2023; 9, e20225. [DOI: https://dx.doi.org/10.1016/j.heliyon.2023.e20225]

24. Wallace, L.; Lucieer, A.; Watson, C.; Turner, D. Development of a UAV-LiDAR System with Application to Forest Inventory. Remote Sens.; 2012; 4, pp. 1519-1543. [DOI: https://dx.doi.org/10.3390/rs4061519]

25. Malinowski, R.; Heckrath, G.; Rybicki, M.; Eltner, A. Mapping Rill Soil Erosion in Agricultural Fields with UAVborne Remote Sensing Data. Earth Surf. Process. Landf.; 2023; 48, pp. 596-612. [DOI: https://dx.doi.org/10.1002/esp.5505]

26. Niculiță, M.; Mărgărint, M.C.; Tarolli, P. Using UAV and LiDAR Data for Gully Geomorphic Changes Monitoring. Developments in Earth Surface Processes; Elsevier: Amsterdam, The Netherlands, 2020; Volume 23, pp. 271-315. ISBN 978-0-444-64177-9

27. Li, P.; Li, D.; Hu, J.; Fassnacht, F.E.; Latifi, H.; Yao, W.; Gao, J.; Chan, F.K.S.; Dang, T.; Tang, F. Improving the Application of UAV-LiDAR for Erosion Monitoring through Accounting for Uncertainty in DEM of Difference. Catena; 2024; 234, 107534. [DOI: https://dx.doi.org/10.1016/j.catena.2023.107534]

28. Wang, H.; Pang, G.; Yang, Q.; Long, Y.; Wang, L.; Wang, C.; Hu, S.; Wang, Z.; Yang, A. Effects of Slope Shape on Soil Erosion and Deposition Patterns Based on SfM-UAV Photogrammetry. Geoderma; 2024; 451, 117076. [DOI: https://dx.doi.org/10.1016/j.geoderma.2024.117076]

29. Borrelli, L.; Conforti, M.; Mercuri, M. LiDAR and UAV System Data to Analyse Recent Morphological Changes of a Small Drainage Basin. ISPRS Int. J. Geo-Inf.; 2019; 8, 536. [DOI: https://dx.doi.org/10.3390/ijgi8120536]

30. Goodwin, N.R.; Armston, J.D.; Muir, J.; Stiller, I. Monitoring Gully Change: A Comparison of Airborne and Terrestrial Laser Scanning Using a Case Study from Aratula, Queensland. Geomorphology; 2017; 282, pp. 195-208. [DOI: https://dx.doi.org/10.1016/j.geomorph.2017.01.001]

31. Pelletier, J.D.; Orem, C.A. How Do Sediment Yields from Post-Wildfire Debris-Laden Flows Depend on Terrain Slope, Soil Burn Severity Class, and Drainage Basin Area? Insights from Airborne-LiDAR Change Detection. Earth Surf. Process. Landf.; 2014; 39, pp. 1822-1832. [DOI: https://dx.doi.org/10.1002/esp.3570]

32. Giménez, R.; Marzolff, I.; Campo, M.A.; Seeger, M.; Ries, J.B.; Casalí, J.; Álvarez-Mozos, J. Accuracy of High-Resolution Photogrammetric Measurements of Gullies with Contrasting Morphology. Earth Surf. Process. Landf.; 2009; 34, pp. 1915-1926. [DOI: https://dx.doi.org/10.1002/esp.1868]

33. Basher, L.; Betts, H.; Lynn, I.; Marden, M.; McNeill, S.; Page, M.; Rosser, B. A Preliminary Assessment of the Impact of Landslide, Earthflow, and Gully Erosion on Soil Carbon Stocks in New Zealand. Geomorphology; 2018; 307, pp. 93-106. [DOI: https://dx.doi.org/10.1016/j.geomorph.2017.10.006]

34. Cui, Y.; Bao, P.; Xu, C.; Ma, S.; Zheng, J.; Fu, G. Landslides Triggered by the 6 September 2018 Mw 6.6 Hokkaido, Japan: An Updated Inventory and Retrospective Hazard Assessment. Earth Sci Inf.; 2021; 14, pp. 247-258. [DOI: https://dx.doi.org/10.1007/s12145-020-00544-8]

35. Kawamura, S.; Kawajiri, S.; Hirose, W.; Watanabe, T. Slope Failures/Landslides over a Wide Area in the 2018 Hokkaido Eastern Iburi Earthquake. Soils Found.; 2019; 59, pp. 2376-2395. [DOI: https://dx.doi.org/10.1016/j.sandf.2019.08.009]

36. Wang, F.; Nam, K. Landslide Disasters Caused by the 2018 Eastern Iburi Earthquake in Hokkaido Japan and the Countermeasures to Completely Prevent the Similar Disasters in the Future. Progress in Landslide Research and Technology, Volume 1 Issue 1, 2022; Sassa, K.; Konagai, K.; Tiwari, B.; Arbanas, Ž.; Sassa, S. Progress in Landslide Research and Technology Springer International Publishing: Cham, Switzerland, 2023; pp. 305-313. ISBN 978-3-031-16897-0

37. Fujiwara, S.; Nakano, T.; Morishita, Y.; Kobayashi, T.; Yarai, H.; Une, H.; Hayashi, K. Detection and Interpretation of Local Surface Deformation from the 2018 Hokkaido Eastern Iburi Earthquake Using ALOS-2 SAR Data. Earth Planets Space; 2019; 71, 64. [DOI: https://dx.doi.org/10.1186/s40623-019-1046-2]

38. Yamagishi, H.; Yamazaki, F. Landslides by the 2018 Hokkaido Iburi-Tobu Earthquake on September 6. Landslides; 2018; 15, pp. 2521-2524. [DOI: https://dx.doi.org/10.1007/s10346-018-1092-z]

39. Zhao, B.; Wang, Y.; Feng, Q.; Guo, F.; Zhao, X.; Ji, F.; Liu, J.; Ming, W. Preliminary Analysis of Some Characteristics of Coseismic Landslides Induced by the Hokkaido Iburi-Tobu Earthquake (5 September 2018), Japan. Catena; 2020; 189, 104502. [DOI: https://dx.doi.org/10.1016/j.catena.2020.104502]

40. Diara, F.; Roggero, M. Quality Assessment of DJI Zenmuse L1 and P1 LiDAR and Photogrammetric Systems: Metric and Statistics Analysis with the Integration of Trimble SX10 Data. Geomatics; 2022; 2, pp. 254-281. [DOI: https://dx.doi.org/10.3390/geomatics2030015]

41. Štroner, M.; Urban, R.; Línková, L. A New Method for UAV Lidar Precision Testing Used for the Evaluation of an Affordable DJI ZENMUSE L1 Scanner. Remote Sens.; 2021; 13, 4811. [DOI: https://dx.doi.org/10.3390/rs13234811]

42. Merwade, V. Watershed and Stream Network Delineation Using ArcHydro Tools. Sch. Civ. Eng. Purdue Univ.; 2012; pp. 1-7. Available online: https://d1wqtxts1xzle7.cloudfront.net/41084682/terrain_processing-libre.pdf?1452737007=&response-content-disposition=inline%3B+filename%3DTerrain_Processing.pdf&Expires=1746083758&Signature=RS11xVgQmE1IbDk8FSEcwB1SBhavDMamjaXokccRI0YIabO92rPnvFEQZsJaCl1TW9xFuvA1KuoJsF5ueTO8LNNG1y3~cRVB34iiVcg4dbxi92TjLYQgMdsk5AZhGqHTQFeFUUT628Sg8ac3yO7sAUytP~z1R3ZLE96A5vXI~-7RH3Nj~rPKsJl~mkzWwTBNBLQpCZgK3pj6cE2sn6kAFWCllj4sITfC8PMnYb~fm4spRstGreWH4nZ-idcCD4HOC4DQEiLR~dK0S70P~3OZromnKXpmdMbvefgshSKVHsfI-f-d10BxWzCkn9yL6lPFFilxB6puJmh-4EMehUtNqw__&Key-Pair-Id=APKAJLOHF5GGSLRBV4ZA (accessed on 10 April 2025).

43. Azizian, A.; Shokoohi, A. DEM Resolution and Stream Delineation Threshold Effects on the Results of Geomorphologic-Based Rainfall Runoff Models. Turk. J. Eng. Environ. Sci.; 2014; pp. 1-15. Available online: https://d1wqtxts1xzle7.cloudfront.net/36680640/muh-38-1-5-1401-13-libre.pdf?1424288807=&response-content-disposition=inline%3B+filename%3DDEM_resolution_and_stream_delineation_th.pdf&Expires=1746083899&Signature=Qv8l2A7rWbJryv0NU2jXjD41mwqREZ~XP-tF-goDo4N7UHTR3XdZCqbqMvFWZ56kZdx7K2sHT-lt7P6O0RiKjvZsdonGvBx3uCPh3e5kQpQLUSO6y7nI1Pm3mW6mtwBN1~wet3EfIvSHMZBfKfuMb4HB0TZTKWVnWlaMVB7bGko5zWydCL4ql9P0arIU2VS~~cYGUnTwVMXGWqGNKb38KjTPkfixtzvJ6x6C9T6wKtilvKInikBQTQWTfoBvA6UcXZfeW~mCcFM0-ou1B6S90m0BQL9wodkVqPiKSszelDtaVNaNxr0gl6J1PS62LszuNkDMHWWXbg75Tdia9jTMJg__&Key-Pair-Id=APKAJLOHF5GGSLRBV4ZA (accessed on 10 April 2025). [DOI: https://dx.doi.org/10.3906/muh-1401-13]

44. Bazzoffi, P. Measurement of Rill Erosion through a New UAV-GIS Methodology. Ital. J. Agron.; 2015; 10, pp. 1-18. [DOI: https://dx.doi.org/10.4081/ija.2015.708]

45. Kou, P.; Xu, Q.; Yunus, A.P.; Ju, Y.; Guo, C.; Wang, C.; Zhao, K. Multi-Temporal UAV Data for Assessing Rapid Rill Erosion in Typical Gully Heads on the Largest Tableland of the Loess Plateau, China. Bull. Eng. Geol. Environ.; 2020; 79, pp. 1861-1877. [DOI: https://dx.doi.org/10.1007/s10064-019-01631-x]

46. Brecheisen, Z.S.; Richter, D.d.B. Gully-Erosion Estimation and Terrain Reconstruction Using Analyses of Microtopographic Roughness and LiDAR. Catena; 2021; 202, 105264. [DOI: https://dx.doi.org/10.1016/j.catena.2021.105264]

47. Perroy, R.L.; Bookhagen, B.; Asner, G.P.; Chadwick, O.A. Comparison of Gully Erosion Estimates Using Airborne and Ground-Based LiDAR on Santa Cruz Island, California. Geomorphology; 2010; 118, pp. 288-300. [DOI: https://dx.doi.org/10.1016/j.geomorph.2010.01.009]

48. Gomez, J.A.; Kamran-Pishhesari, A.; Sattarvand, J. Automated Rill Erosion Detection in Tailing Dams Using UAV Imagery and Machine Learning. Arab. J. Sci. Eng.; 2024; pp. 1-16. [DOI: https://dx.doi.org/10.1007/s13369-024-09550-0]

49. Dai, W.; Na, J.; Huang, N.; Hu, G.; Yang, X.; Tang, G.; Xiong, L.; Li, F. Integrated Edge Detection and Terrain Analysis for Agricultural Terrace Delineation from Remote Sensing Images. Int. J. Geogr. Inf. Sci.; 2020; 34, pp. 484-503. [DOI: https://dx.doi.org/10.1080/13658816.2019.1650363]

50. Sivanandam, P.; Turner, D.; Lucieer, A.; Segaran, D.R.C.; Virtue, J.; Melville, B.; McCallum, K. Drone Data Collection Protocol Using DJI Matrice 300 RTK: Imagery and Lidar. 2022; Available online: https://www.tern.org.au/wp-content/uploads/20230829_drone_data_collection.pdf (accessed on 8 May 2023).

51. Brodu, N.; Lague, D. 3D Terrestrial Lidar Data Classification of Complex Natural Scenes Using a Multi-Scale Dimensionality Criterion: Applications in Geomorphology. ISPRS J. Photogramm. Remote Sens.; 2012; 68, pp. 121-134. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2012.01.006]