Erosion Gully Networks Extraction Based on InSAR Refined Digital Elevation Model and Relative Elevation Algorithm—A Case Study in Huangfuchuan Basin, Northern Loess Plateau, China

Abstract

The time-effective mapping of erosion gullies is crucial for monitoring and early detection of developing erosional progression. However, current methods face challenges in obtaining large-scale erosion gully networks rapidly due to limitations in data availability and computational complexity. This study developed a rapid method for extracting erosion gully networks by integrating interferometric synthetic aperture radar (InSAR) and the relative elevation algorithm (REA) within the Huangfuchuan Basin, a case basin in the northern Loess Plateau, China. Validation in the study area demonstrated that the proposed method achieved an F1 score of 81.94%, representing a 9.77% improvement over that of the reference ASTER GDEM. The method successfully detected small reliefs of erosion gullies using the InSAR-refined DEM. The accuracy of extraction varied depending on the characteristics of the gullies in different locations. The F1 score showed a positive correlation with gully depth (R² = 0.62), while the fragmented gully heads presented a higher potential of being missed due to the resolution effect. The extraction results provided insights into the erosion gully networks in the case study area. A total of approximately 28,000 gullies were identified, exhibiting pinnate and trellis patterns. Most of the gullies had notable intersecting angles exceeding 60°. The basin’s average depth was 64 m, with the deepest gully being 140 m deep. Surface fragmentation indicated moderate erosive activity, with the southeastern loess region showing more severe erosion than the Pisha sandstone-dominated central and northwestern regions. The method described in this study offers a rapid approach to map gullies, streamlining the workflow of erosion gully extraction and enabling efficiently targeted interventions for erosion control efforts. Its practical applicability and potential to leverage open-source data make it accessible for broader application in similar regions facing erosion challenges.

1. Introduction

Gully erosion is a widespread land degradation form that is occurring around the globe [1]. It is a large source of sediment yield, usually accounting for 20~90% depending on the landscape [2,3]. Consequently, the nutrient content, like the soil organic carbon, in the topsoil decreases due to the transportation of soil and sediment [4]. The development of erosion gullies also shrinks the area of usable land by fragmenting the morphology of a surface that was once flat, forming a rugged terrain consisting of numerous gullies that link into networks [5,6]. Gully erosion has become a key concern in many regions around the world over the past century [7]. Therefore, monitoring the development of erosion gullies via extracting and mapping has become vital for enacting prevention strategies.

Efforts to monitor erosion gully have evolved over time. Traditional methods include field surveys and manual extraction of erosion gullies from satellite images, which are labor-intensive and time-consuming methods [6]. Zhao et al. [8] developed a random point sampling method for the rapid estimation of the presence of erosion gullies. However, detailed information on gully characteristics such as morphology, depth, and width is crucial for water and soil conservation studies. Therefore, the extraction of entire gullies within specific regions is still necessary. Technological advancements, such as LiDAR and unmanned aerial vehicles (UAVs), have provided high-resolution data, enabling object-based gully extraction [9,10]. Existing efforts that have been dedicated to erosion gullies extraction can be classified into two categories.

Category I is mainly supported by computer vision and deep learning. UAVs and open archives of high-resolution optical satellites (e.g., Google Earth Imagery, Esri World Imagery, etc.) provide orthographic observation of erosion gullies. By applying deep learning models with different architectures (e.g., SegNet, U-net, R2-Unet, etc.), computers obtain the features of a gully and classify the gully pixels on images to accomplish erosion gully extraction tasks [11,12]. Although these methods can generate fine boundaries, the labeling and training process is time-consuming and sensitive to the image color uniformity, contrast, and quality. For tasks covering large areas, small UAVs are unsuitable due to economic decline, and the use of high-resolution satellite images greatly increases the cost. Additionally, three-dimensional attributes such as the depth of a gully cannot be obtained from images. Xue et al. suggested a method to obtain three-dimensional river networks with a relatively novel effect. However, this workflow is mainly designed for larger river systems and may not be able to obtain information from small erosion gully networks [13].

Category II takes advantage of the digital elevation models (DEMs), a rasterized digital expression of the Earth’s surface terrain, to detect abrupt topographic points as the identifier of gullies. These methods are elaborately designed mathematically, such as the window sliding algorithm developed on finite morphologic element theory [14], the active contours method using the snake model [15], edge detection based on the LOG operator [16], and the Marr–Hildreth operator [17]. However, there is a tradeoff between accuracy and computation cost when applied to a large area. Moreover, the rapid development of erosion gullies in highly eroded areas may outdate the available public DEM datasets, hindering the quasi-real-time extraction and monitoring of erosion gullies at the regional scale.

Notably, the mentioned methods in these two categories are currently limited to smaller areas due to time and cost constraints. There is a pressing need to develop extraction methodologies that consider gully networks as an interconnected system [18]. Such approaches would enable a more comprehensive understanding of the advancing progress of gully erosion. Additionally, existing methods encounter challenges in identifying gully type (e.g., gully bank, gully floor, etc.) because a physically meaningful classification scheme and a corresponding feasible method are still lacking [6]. Currently, the definition of gullies that can distinguish their boundary from fragmented landscapes remains ambiguous, hindering the intactness of gullies’ extraction from available data. While the lower bound of the critical cross-sectional area of a gully is usually defined as the size that conventional tillage practices cannot erase [19], the upper bound remains unclear in the existing literature [6]. In the Chinese Loess Plateau, the shoulder line is commonly treated as the boundary of erosion gullies [20], as the slope below this line increases dramatically (usually > 35°).

Interferometric synthetic aperture radar (InSAR) can detect deformations of the Earth’s surface by comparing the phase difference of two different SAR images with slightly different antenna positions [21]. With its all-day functioning, short-interval revisiting, cloud-penetrating ability, and high-resolution detection, SAR provides considerable potential for time sequence DEM reconstructing [22] and in long-term gully erosion studies. Moreover, leveraging the phase information offered by InSAR enables the measurement of slight topographic variations and the generation of DEM with enhanced relief details [23]. These additional details facilitate more accurate erosion gully extraction. Additionally, SAR satellites provide cost-effective data covering a wide range of areas [24,25], allowing for the rapid acquisition of DEMs and improvement in DEM quality, particularly in capturing surface gully features. This approach partially overcomes the limitations of traditional DEM extraction methods. The existing literature has proven that InSAR-based methods are efficient in detecting surface deformation and thus reconstructing DEM with relatively high accuracy. Berardino et al. [26] proposed a small orbital separation baseline (SBAS) to limit spatial decorrelation and detect ground deformation in southern Italy on a small map scale. Zhao et al. [27] examined a large-map-scale loess landslide in Shaanxi, China, using improved SBAS by checking DEM errors with a coherence threshold and determining interferogram quality in advance to suppress imprecise elevation correction. The launch of Sentinel-1 (S1) made the open access of spaceborne radar imagery possible. Even though this project was not primarily designed for DEM reconstruction, the available imagery still meets the basic requirements for DEM reconstruction with a careful selection of image pairs [23].

The relative elevation algorithm (REA) removes the downward trend in floor elevation by respecting the baselines along the bottom of gullies [28]. It is generally known that gully floors always have a downward trend (the decreasing trend of floor elevation) extending from upper gully heads to lower confluences or outlets because of the close relationship between gully erosion and hydrologic processes. This limits small-map-scale gully extraction, as the descending elevation along the longitudinal gully-flow profile cannot be ignored, and the contours of the gully bank are always convergent. By applying REA, a floor-equal elevation model enhancing the cross-sectional characteristics of gullies can be derived. This algorithm is often applied to fluvial geomorphology as a cartographic expression for visual enhancement, but it also shows great potential in gully extraction. With the removal of the downward trend, the elevations at different terrain locations are numerically clustered, facilitating the classification of gully types and overcoming the limitation of the traditional methods that only extract gully boundaries without acquiring three-dimensional information from the inner gully area.

The Loess Plateau (LP) in north-central China is characterized by extensive gully formation, primarily due to the erodible loess soil and prolonged inappropriate land management practices [29]. The Huangfuchuan Basin, recognized as one of the most eroded areas within the Loess Plateau, is significant as it is designated as one of China’s eight national key management areas [30]. These erosional activities have significantly increased sediment transport in the Yellow River, contributing to the increased silting of river beds and heightened risks of ecological destruction [31,32]. Consequently, the development of extensive gully networks in this basin has established it as a prime location for studying erosion gully extraction, given its typical gully landforms and practical significance (Figure 1). In this study, we developed a workflow that integrates InSAR-reconstructed DEMs with REA for the rapid extraction of gully networks in a selected typical gully erosion watershed in the LP. The objectives of this study were to reconstruct digital elevation models (DEMs) and extract gully networks in the study area, and to analyze the spatial characteristics of the gully networks, as well as the performance of the method.

2. Materials and Data

2.1. Study Area

The Huangfuchuan River (HFC) basin (Figure 2) was selected as the study area to represent a typical loess gully basin suffering from intensive gully erosion. The HFC originates in northern Jungar Banner, Inner Mongolia Autonomous Region, and flows 125.1 km into the Yellow River in Fugu County, Shaanxi Province. The HFC drainage basin (39°12.7′~39°59.6′N, 110°20.3′~111°15.0′E) covers 3240 km² in the northern semiarid Loess Plateau region, with elevations ranging from 820 to 1440 m. This region has been experiencing uplift since the Quaternary period, resulting in diverse forms of erosion gullies [33]. Gully depths average 40–70 m, with some exceeding 100 m. The composition of the soil within the HFC basin predominantly consists of Pisha sandstone and loess. Pisha sandstone is a local soil texture that mainly covers the northern Loess Plateau. The feldspar in the Pisha sandstone, primarily weathered by kaolinite, exhibits low erosion resistance and serves as the primary mineral contributing to the sandstone’s susceptibility to erosion [34,35]. Loess features a unique structure where clay-size particles coat silt grains and gather at quartz contacts, resulting in high porosity and macroscopic pores. Its soil properties are influenced by microstructure, making it prone to collapse and sliding when wet [36,37]. Both these soil types exhibit weak bonding and low compressive strength, which contribute to the basin’s susceptibility to significant soil erosion [38]. The annual erosion modulus has been estimated as between 3000 and 4000 tons per square kilometer [39]. The northern region of the basin is characterized by a prevalence of Pisha sandstone, while the southern region is dominated by collapsible loess. Vegetation cover, primarily Caragana korshinskii, coniferous Prunus armeniaca (apricot), and Pinus massoniana (masson pine), decreases from southeast to northwest. Crop fields are distributed near the main HFC channel and intergully tableland area, having a gentle slope. The predominant crop types are Zea mays (maize) and Panicum miliaceum (millet). The average annual temperature is 9.1 °C, and precipitation is 350~450 mm, with extreme rainfall in July (mean 92 mm) occurring due to the continental monsoon. This results in an average of 61.2 million tons of sediment yield annually, with heightened summer fluvial erosion. Extensive water and soil conservation projects, chiefly check dams (sediment-trapping dams) and terraced fields, have been implemented across the basin. These distinct characteristics of the HFC make it an ideal area for researching hydraulic erosion gullies in the Loess Plateau.

2.2. Data

2.2.1. SAR Data

Sentinel-1 interferometric wide (IW) swath data distributed in Level-1 single look complex (L1 SLC) products were used in this study. Original SAR data were acquired by the Terrain Observation with Progressive Scans SAR (TOPSAR) on the Sentinel-1. The spatial resolution was 5 m × 20 m in a single look. Data from the L1 SLC products were 2.3 m × 14.1 m gridded for universal agreement in range and azimuth while preserving phase information. As suggested by Alessandro et al. [40], suitable temporal and perpendicular baseline (150–300 m) are two vital requirements that need to be met when selecting suitable image pairs, because improper selection may lead to failure in DEM reconstruction due to the signal being out of phase. Moreover, low water vapor in the atmosphere and low vegetation coverage on land are also suggested to avoid phase delay and temporal decorrelation. Following these principles, we selected and downloaded available data (

Table 1. SAR data details.

ID	Acquired Time	Path	Frame	Beam Mode	Polarization	Baseline
						Perpendicular (m)	Temporal (d)
1	24 November 2022 10:30:34	113	126	IW	VV + VH	127	24
2	18 November 2022 10:30:33	113	126	IW	VV + VH

) from the Data Search Vertex of Alaska Satellite Facility (https://search.asf.alaska.edu, accessed on 10 December 2023) [41].

2.2.2. Auxiliary Data

Reference ASTER GDEM V3 was downloaded from www.earthdata.nasa.gov (accessed on 10 December 2023) [42]. Version 3 of ASTER GDEM was generated using stereo data and released in 2019. According to its official user guide, the standard deviation of the elevation error has decreased to 12.1 m, an improvement compared to that of the previous versions [43]. Several studies have evaluated the accuracy of ASTER GDEM, reporting a general vertical root mean square error (RMSE) ranging from 8.52 to 30.11 m, depending on the complexity of the topography [44,45,46]. In comparison to other available DEMs such as SRTM DEM, AW3D DEM, and ALOS DEM, ASTER DEM is more suitable for assisting the phase unwrapping process in the InSAR workflow, particularly in mountainous areas [47]. Although the recently released Copernicus DEM was reported to have high accuracy [48], the central and northern regions of China have not been sufficiently observed at present [49]. Therefore, ASTER DEM was selected as a reference in this study to assist in the reconstruction of the InSAR DEM.

Landsat-9 pan-sharpened imagery was selected for manual extraction of the ground truth data of the gullies because it has a similar spatial resolution to the refined DEM. To simplify the processing workflow and focus on the main tasks of this study, we used the ©ESRI ArcGIS Online service for instant access to the preprocessed Landsat-9 imagery data (https://www.arcgis.com/home/item.html?id=a7412d0c33be4de698ad981c8ba471e6, accessed on 10 December 2023). This online platform offers a comprehensive Landsat-9 imagery collection that has been pan-sharpened and processed in a multitemporal manner. This imagery is refreshed daily using newly acquired images directly obtained from the USGS Landsat collection. It is dynamically generated in natural color with dynamic range adjustment (DRA) to facilitate visualization and analytical purposes [50]. By leveraging this online service, we efficiently accessed the prepared Landsat imagery, streamlining the complexities of traditional image processing workflows.

3. Methods

This study’s workflow comprised two stages. The initial stage involved DEM reconstruction utilizing SAR data (shown in yellow boxes in Figure 3). The ASTER DEM served as the basis for georeferencing. Selected pairs of Sentinel-1 SAR images underwent a pixel-by-pixel multiplication operation to generate an interferogram for surface deformation detection. Subsequently, a coherence image was derived from the SAR image pairs to evaluate the interferogram’s quality. Pixels with coherence values less than 0.15 were excluded. The filtered image was then subjected to phrase unwrapping and geocoding processes sequentially. Finally, the refined DEM was generated. The second stage entailed the extraction of erosion gully networks from the refined DEM (shown in green boxes in Figure 3). Hydrologic analysis was initially conducted to generate line features of streams. These features facilitated the detection of the gully floor elevation, which was used to establish the reference plane. The relative elevation algorithm was then employed to produce a relative elevation model (REM). Ultimately, the erosion gully network masks were extracted from the REM.

3.1. DEM Reconstruction Using SAR Data

Two different radar images with known perpendicular baselines can be used to estimate terrain information from the Earth’s surface by calculating the different path lengths of the two in-phase radar echo signals received by the antennas [51]. The varied delays of the radar echo signal received by the antennas contain the terrain change information. For example, a longer traveling time of a signal indicates it touches a concave-down surface $\left(\mathsf{\Delta }d\right)$ rather than a flat one before returning to the antenna (Figure 4). Respecting the sinusoidal nature of the radar signal, the signal ($S$) associated with time ($t$) can be expressed as in Equation (1) [52]:

$$S\left(t\right)=A\cos \left[2\pi f_0\right(t-\tau \left)\right]=A\cos \left[2\pi f_{0}t-\left(-\frac{2\pi }{\lambda }2R\right)\right]$$

(1)

where $A$ is the amplitude, $\tau =\frac{2R}{c}$ is the time required for signals to travel back to the antenna, $\lambda =\frac{c}{f_0}$ is the wavelength ($c$ is the speed of an electromagnetic wave), and $-\frac{2\pi }{\lambda }2R$ is known as the phase delay ($\phi$).

Above is the original theory of InSAR, assuming that initial radar signals are sent by the antennas at the same time. However, Sentinel-1 was not designed for DEM reconstructing tasks and can only work in monostatic mode [23]. Given that the Earth’s surface does not change within a short time, repeat observations still follow the principle of InSAR so can be used to measure surface elevation [23,53]. The main steps of DEM reconstruction using Sentinel-1 data in this study are briefly explained below. Note that SAR data were processed via the SARscape 5.6 framework in ENVI 5.6 software for DEM reconstruction. For details about the processing process and intermediate results, please see Appendix A for reference.

3.1.1. Retrieving Interferogram

As mentioned above, gridded S1 images preserve the phase information about the last fraction of a complete radar wave period. To retrieve the phase information, the selected images for DEM reconstruction were cross-multiplied pixel by pixel to create an interferogram. It was a representation combining the phase variation of topography, deformation, Earth’s curvature, atmospheric effect, and noise [23]. Due to the pixel-based calculation, spatial registration must be precisely conducted in advance to obtain phase differences from the pixel pairs. In this study, we accomplished spatial registration by referencing an existing DEM for homologous image point matching.

3.1.2. Calculating Coherence Image

The coherence image between the two image was calculated to inspect the quality of the interferogram. The coherence ($\gamma$) is defined by its correlation coefficient as in Equation (2) [54]:

$$\gamma =\frac{E\left\{s_1{s}_2^{*}\right\}}{\sqrt{E\left\{{\left|s_1\right|}^2\right\}E\left\{{\left|s_2\right|}^2\right\}}}$$

(2)

where $s_1$ and $s_2$ represent the SAR signals captured and rasterized on the two images pair, * is the complex conjugate operator, and $E$ means mathematical expectation. The value range of $\gamma$ is $[0,1]$. Bare land and urban areas usually have a high coherence value, while rapidly changing surfaces like water and vegetation usually have a low coherence value [23]. Signals with low coherence are unusable for DEM reconstruction. A threshold between 0.1 and 0.4 is widely accepted to filter out unusable information. Considering the number of valid pixels and the low vegetation coverage in our study area, we used 0.15 as the threshold.

3.1.3. Unwrapping Phase

The unwrapping phase can convert SAR signals from cyclic wave measurement to continuous height measurement [55]. The pixel values on an unwrapped interferogram represent topographic variation or relative height between the image pair. For proper unwrapping, we applied a Goldstein filter to the interferogram in advance to enhance the signal–noise ratio [56]. Although applying the filter may introduce errors by suppressing potential true but minor phase variations [57], the fast Fourier transformation (FFT)-based filter remains stable and functional within the usable signal range. Given the low vegetation coverage and arid climate in the study area, the benefits of this filter outweighed its drawbacks.

3.1.4. Geocoding Unwrapped Phase

Geocoding was performed to translate the unwrapped phase to metric measurement representing the true Earth’s surface height along the line between the pixel and antenna [40]. In this study, we used ASTER DEMs as the reference to rectify geometry disorders during imaging and finally retrieved reconstructed a DEM georeferenced to WGS 1984/EGM 1996.

3.2. Gullies Extraction

The erosion gully refers to the inner area between shoulder lines [15], as illustrated in Figure 5. Given the close relationship between erosion gullies and streams, we introduced the following terminology from hydrology to describe a gully network: gully heads are the farthest upstream locations of gullies, junctions are the locations where two gullies join to form a downstream one, the floor is the lower part of the inner area, the bank is the part between the shoulder line and the floor, and depth is the vertical distance from the lowest point on the floor to the shoulder line.

The reconstructed DEMs were projected onto the WGS 1984 coordinate system and EGM 1996 geoid, as previously mentioned. Each DEM pixel represented the surface elevation relative to the geoid. Due to the downward-inclined nature of an erosion gully, the elevation of both the gully floor and shoulder exhibit a downward trend along the center line of erosion gullies. Thus, applying constant elevation thresholds across the landscape to extract gullies was impractical, as local elevation variations within gullies may be obscured by global terrain undulations. Drawing on the concept of the relative elevation model in cartography, we developed an REA method to eliminate the effect of global terrain. Specifically, gully floors were extended as a reference plane, and DEM elevations were vertically aligned to the reference plane. Pixel values were assigned to represent gully depths rather than absolute elevation. Consequently, the histogram distribution of these pixel values should match the topographic relief of the gullies and could help delineate inner gullies.

3.2.1. Detecting Floor Elevation of Erosion Gullies

The floor elevations of erosion gullies can be detected using hydrological analysis methods, as gullies are primarily shaped by hydraulic erosion. SAR-derived DEMs may contain sinks created by errors. A sink refers to a cell with an undefined drainage direction. In this study, sinks were first filled to avoid unreasonable flows. The flow direction and accumulation raster were then derived from the filled DEM using the multi-flow direction algorithm [58]. Subsequently, a stream network was generated by selecting a threshold on the flow accumulation raster. The turning point where river network density ceases to substantially increase with the threshold is often used as an optimal threshold [59]. However, to better detect small gullies, a relatively smaller threshold could be selected. As such, sampling points (orange points shown in Figure 6a) were set at equal intervals along the generated river network to extract gully floor elevations from the DEM.

3.2.2. Creating a Reference Plane

The reference plane was created by interpolating floor elevations across the study area using the ANUDEM 5.3 program built in ArcGIS Pro 3.0. ANUDEM was specifically developed to improve the interpolation of DEMs from vector inputs by adaptively managing local topographic features. It incorporates a drainage enhancement algorithm that preserves hydrological features such as streamlines and drainage structures. Additionally, it removes unnecessary kinks and generalizes the data to better manage feature conflicts [60]. In comparison to the traditional rasterizing method that utilizes the Triangular Irregular Network (TIN), significant improvements in the interpolated DEM were observed, particularly in the LP region [61]. In this study, perpendicular lines to flow direction were created at each sampling point. The length of these lines was determined by the width of the erosion gully. These lines acted as the contour lines of the reference plane, with elevation values inherited from the sampling points. Both sampling points and perpendicular lines were used as inputs of the ANUDEM program to generate the reference plane (Figure 6b).

3.2.3. Extracting Mask of Erosion Gullies Using REA

The REA is the subtraction of the elevation ($E$) of the original DEM from that of the reference plane (Equation (3)):

$$E_{\mathrm{R}\mathrm{E}\mathrm{M}}=E_{\mathrm{R}\mathrm{e}\mathrm{f}\mathrm{P}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{e}}-E_{\mathrm{D}\mathrm{E}\mathrm{M}}$$

(3)

In the REM (Figure 6c), a value of 0 corresponds to the floor of the erosion gullies, and negative values correspond to the descending height of the gully banks. After conducting the REA, the base elevation of the gully and its banks on both sides were adjusted to a consistent elevation on the reference plane. This operation standardized the base elevation of the gully, thus eliminating the downward trend. As the downward trend was removed, the undulation of the gullies was reinforced. Thus, the gullies could be extracted from the terrain using the natural breaks method. The natural breaks method was first developed for unsupervised classification in choroplethic maps [62,63]. It is a data-clustering method that minimizes elevational discrepancies within categories while maximizing those between categories. It is reported that data can be properly classified using this method, with minimal information loss [64]. The method and its variants have been successfully applied in the field of remote sensing and GIS in classification and object extraction tasks based on quantitative data [65,66,67]. After REA processing, DEM data in gully areas exhibited aggregation patterns in the gully floor and inter-gully regions. This resulted from REA normalizing the relief to a uniform base elevation, causing these two areas to share similar elevations. By applying the natural breaks method, these two categories representing the gully floor and intergully areas could be distinguished. Considering the cross-sectional characteristics of the gully, the remaining data between these two categories represented the gully bank. Converting these three categories into numeric sequences in ascending order of elevation, they, respectively, represented the gully floor, gully bank, and intergully area. Specifically, the first and second categories together represented the inner-gully area. (Figure 6d).

3.3. Validation

The ground truth mask of the gully was extracted manually from Landsat-9 pan-sharpened images at random locations within the basin. Pixels of the REA-extracted gully located in the ground truth mask were considered to be correct. Here, we applied the F1 score (Equation (4)) for validation.

$$\mathrm{F}1 \mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\frac{2\times \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}}$$

(4)

where $\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}$ indicates the proportion of true positive pixels among all the positive pixels generated from the REA extraction, while $\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}$ indicates the proportion of true positive pixels identified using the REA extraction out of all the actual positive pixels in the ground truth. The calculation of these two indicators is given in Equations (5) and (6), respectively.

$$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}=\frac{{Count}_{\mathrm{T}\mathrm{P}}}{{Count}_{\mathrm{T}\mathrm{P}}+{Count}_{\mathrm{F}\mathrm{P}}}$$

(5)

$$\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}=\frac{{Count}_{\mathrm{T}\mathrm{P}}}{{Count}_{\mathrm{T}\mathrm{P}}+{Count}_{\mathrm{F}\mathrm{N}}}$$

(6)

where ${Count}_{\mathrm{T}\mathrm{P}}$ represents the count number of true positive pixels, ${Count}_{\mathrm{F}\mathrm{P}}$ represents the count number of false positive pixels, and ${Count}_{\mathrm{F}\mathrm{N}}$ represents the count number of false negative pixels. The accuracy was also computed as the ratio of correctly extracted pixels to the overall pixel count for purposes of comparative analysis. This calculation is presented in Equation (7).

$$\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{Coun{t}_{\mathrm{T}\mathrm{P}}+Coun{t}_{\mathrm{T}\mathrm{N}}}{{Count}_{\mathrm{T}\mathrm{P}}+{Count}_{\mathrm{F}\mathrm{P}}+{Count}_{\mathrm{T}\mathrm{N}}+{Count}_{\mathrm{F}\mathrm{N}}}$$

(7)

where ${Count}_{\mathrm{T}\mathrm{N}}$ represents the count number of true negative pixels. Other parameters have the same meaning as in Equations (5) and (6). To assess the effectiveness of the method described above, it was applied again to the ASTER DEM. The validation metrics for the extracted erosion gully networks were then calculated using the same validation method to establish a reference result, which was subsequently utilized for comparative analysis.

4. Results

4.1. Extraction of Erosion Gullies

To demonstrate the above procedure of erosion gully extraction, we selected the HFC basin as the case region for testing. A DEM with finer spatial resolution (Figure 7c) was first reconstructed using InSAR techniques. This DEM showed enhanced terrain details with less noise compared to the reference ASTER DEM (Figure 7b). As mentioned in Section 3.2.1, a threshold of flow accumulation was determined as 100 to better detect gully heads. Tests showed that a smaller threshold could lead to the identification of stream features with lengths shorter than the sampling intervals (50 m) of floor elevation (Figure 7d). The REA was then used, and its result was classified using the natural breaks method into three parts, namely, the gully floors, the gully banks, and the intergully tablelands (Figure 7e). The initial gully mask was created by assigning null values to the pixels classified as intergully tablelands. Continuous flat areas over 1 km² (e.g., large river channels, large tableland, etc.) were also removed using a sliding circular window with a radius of 38 pixels. The erosion gullies could finally be extracted with this optimized mask (Figure 7f).

4.2. Characteristics of Extracted Erosion Gullies

A total of 28,000 erosion gullies were identified and extracted from the HFC basin. These gullies exhibited pinnate patterns at the basin scale and trellis patterns on localized slopes (Figure 8). Using the Strahler method, the gullies were classified into five hierarchical orders. Note that the fourth- and fifth-order gullies were excluded in the following analysis because they were large river channels and were not part of this study. The results revealed that the intersecting angles at the junctions of the first- and second-, as well as second- and third-order gullies, typically exceeded 60°. Some intersecting angles in the northwestern basin reached 90°. Transverse small-watershed systems due to headward erosion were observed more commonly in this area. REA revealed that the average depth of the basin was 64 m, with 95% of the gullies ranging in depth from 50 to 78 m, and the deepest measuring 140 m.

To examine the spatial characteristics of the gullies, we calculated surface fragmentation by determining the proportion of gully planar area within a 1 km² rectangular region. The average fragmentation was found to be 54%, indicating moderate erosive activities in the basin. The southeastern part of the basin, composed primarily of loess material, exhibited more severe surface fragmentation (over 60%) compared to the central and northwestern regions, where Pisha sandstone was predominant (Figure 9).

4.3. Validation of the Proposed Extraction Method

Twenty-seven square regions (1 × 1 km) within the study area were randomly selected for validation. Ground truth gully masks were manually digitalized from Landsat-9 pan-sharpened images. These images were chosen because their spatial resolution was close to that of the DEM. The accuracy, precision, recall, and F1 score of each validation region were calculated and plotted (Figure 10). The results showed that the method achieved overall good performance. It accurately identified the majority of the gullies, with an average accuracy of 80.88%. Additionally, it efficiently distinguished erosion gullies from the inner gully areas, as evidenced by the relatively high average recall of 87.76%. However, the average precision was lower than the other metrics, reaching a minimum of 54.27%. This indicates that some gullies were missed, possibly due to the coarse spatial resolution of the DEMs, which might have hindered the detection of small-scale topographic relief near the gully heads. Nonetheless, the average F1 score of 81.94% suggested that local errors did not significantly impact the overall results. The method achieved a balanced performance by effectively identifying true positives (real inner gully area) and avoiding false positives.

5. Discussion

The spatial resolution of DEMs may have significant impacts on the accuracy of gully extraction due to the resolution effect [68]. The resolution effect refers to the generalization of terrain details when using a coarser grid (lower spatial resolution) to represent the true topographic relief. It may result in a smoothed representation of the surface [69] and limits the efficiency of gully extraction. The suitable spatial resolution proposed in previous studies for gully information extraction varies from 2 m to 15 m [68,70]. For studies covering large areas, the acquisition of these high-resolution DEMs remains challenging due to their limited availability. The use of drones to obtain DEMs with high spatial resolution can also become less feasible due to time and cost issues. Nonetheless, the extraction of gullies does not necessitate a DEM with precise absolute elevation values. Because the relative elevation difference between the inner- and intergully areas is more crucial for gully extraction purposes. Thus, this study developed a method that integrates InSAR and REA techniques to facilitate rapid gully extraction over large areas. Comparative results showed that the reconstructed DEM with a finer resolution (15 m) improved gully extraction performance, with an F1 score of 81.94% for gullies’ extraction. In contrast, the reference ASTER DEM with a coarser resolution (30 m) achieved an F1 score of 72.17%. Thus, our study confirmed that it is feasible to extract gullies from lower-resolution DEMs by employing suitable information enhancement methods, as suggested by Maerker et al. [71].

Gullies in different locations may have different characteristics due to different erosion magnitudes. This can impact the accuracy of gully extraction [1]. Gully depth usually indicates the erosion rate or development age of a gully [72]. Regression analysis (Figure 11a) revealed a positive correlation between the F1 score and gully depth (R² = 0.62), indicating that older gullies with significant elevation differences were easier to extract. Gully heads presented a higher likelihood of being missed during the extraction process, particularly when fewer than 50 continuous pixels (approximately 1 hm²) were involved (Figure 11b). Gully heads in the loess region exhibited greater fragmentation than those in the Pisha sandstone area. Indoor experiments revealed that sandy soils showed higher soil loss compared to silty soils [73]. Therefore, we could infer that the sandy-soil-dominated Pisha sandstone region was less prone to multidirectional competition at gully heads due to its high erodibility. Consequently, the degree of gully head undercutting was expected to be greater, rendering complete extraction of gullies in this region more feasible.

Future works should focus on enhancing the accuracy of gully head extraction. This can be achieved by refining the spatial resolution of the DEM or by assigning appropriate weights to the gully head area during the generation of the reference plane. This approach would emphasize the topographic gradient at the gully heads, thereby improving the accuracy of their extraction. Moreover, while this study was conducted within the context of the Loess Plateau’s erosion gullies, its applicability to other regions without distinct shoulder lines requires additional validation. Nonetheless, the method proposed here demonstrates the potential for leveraging open-source data for the rapid and extensive mapping of gully networks. The extracted gully networks provided a way to investigate gully erosion advances and geomorphology, from a systematic viewpoint.

6. Conclusions

This study integrated InSAR and REA techniques to develop a rapid method for extracting erosion gully networks. The InSAR technology reconstructed up-to-date DEMs with finer spatial resolution. And, the REA method extracted the gullies by enhancing their vertical elevation difference. The validation of the results indicated that the method achieved an F1 score of 81.94%, representing a 9.77% improvement over the reference ASTER DEM. This suggests that the application of appropriate information enhancement methods to low-resolution DEMs can be effective for gully extraction. The accuracy of the extraction was affected by the characteristics of the erosion gullies at different locations. The F1 score was positively correlated to the depth of the gullies (R² = 0.62), while the fragmented gully heads presented a higher potential of being missed due to the resolution effect. Further insights into the erosion gully networks within the HFC basin were revealed by the extraction results. The identified 28,000 gullies in the HFC basin showed pinnate and trellis patterns. Most of the gullies had notable intersecting angles exceeding 60°. The basin’s average depth was 64 m, with the deepest gully being 140 m. Surface fragmentation indicated moderate erosive activity, with the southeastern loess region showing more severe erosion than the Pisha sandstone-dominated central and northwestern regions. In summary, the method presented in this paper offers a streamlined workflow and relatively good performance for the rapid extraction and mapping of erosion gully networks. It enables quasi-real-time monitoring of gully network development and demonstrates the potential for leveraging open-source data to rapidly map erosion gully networks at the regional scale.