# Quantitative Attribution Analysis of the Spatial Differentiation of Gully Erosion in the Black Soil Region of Northeast China

Authors: Ranghu Wang; Nan Wang; Yingchao Fan; Huan Sun; Hanpei Fu; Jiuchun Yang
Journal: Geofluids
 (2022)

---

## Abstract

Gully erosion is the major soil erosion type in the black soil region of Northeast China. However, studies on multifactor synthesis at a large scale and on the driving mechanism of spatial differentiation are still relatively lacking for gully erosion in this region. In this study, the simulation of gully erosion and its quantitative attribution analysis have been conducted in the Sancha River catchment in Northeast China, based on high-resolution satellite imagery mapping and the geodetector method. A total of 18 indicators in 6 categories, including topography, climate and weather, soil properties, lithology, land use, have been taken into consideration. The influence of each influencing factor and its interactive influence on gully erosion were quantitatively evaluated. The results showed that at the large catchment scale, the submeter images had a strong capacity for the recognition of a permanent gully and obtained satisfactory results. According to the results of the geodetector, lithology and soil type are the main factors that affect the spatial differentiation of gully erosion in the Sancha River basin, because their interpretation power for gully density and gully intensity was close to 10%. The lithology belonged to gray–white matter rhyolite, spherulite rhyolite, and crystal clastic tuff, with the highest gully density and intensity. The interpretation power of the secondary factors, including rainfall erosivity, watershed area, elevation, soil erodibility, land use pattern, slope, and distance from the river, amounted to more than 1%. The interactions among most driving factors showed nonlinear enhancement. The influence of the interaction between lithology and soil type appeared to be the largest. In particular, the lithology of different soil types accounted for 28.7% and 32.5% of the gully density and gully intensity. The interaction of factors had a stronger influence on the spatial differentiation of gully erosion than any single factor.


---

## Body

## 1. Introduction

Gully erosion is the formation and subsequent expansion of erosional channels in the soil as a result of concentrated water flow [1]. Runoff water through the catchment area removes soil from the channels and accumulates in the narrow channels. After several iterations, a gully of considerable depth will appear. Hence, gully erosion is one of the most effective drivers of sediment removal and runoff from highland areas to valley floors [2]. Gully erosion was, is, and will continue to be one of the world’s most important environmental problems, especially in semiarid and arid areas, where soils are suffering from severe gully erosion [3]. For instance, Iran, located in Western Asia, is suffering from serious soil erosion, with approximately 2-2.5 billion tons of soil being lost every year, accounting for approximately 50% of the country’s land area. Gully initiation and development is a natural process that greatly impacts agricultural activities and environmental quality as it promotes land degradation, desertification, and ecosystem disruption [2].In the last 100 years, with the expansion of the agricultural land area, the black soil in Northeast China has suffered serious soil erosion, becoming one of the largest soil erosion regions in China [4]. Several studies have proposed that gully erosion is the major soil erosion type in this region [5, 6]. According to the findings of the first national water conservancy survey, there are nearly 300,000 gullies in the black soil region of Northeast China, causing a loss of cultivated land of approximately 4.83×105 hm2, and 3.62×109 kg of grain is lost annually due to the development of eroded gullies [7]. Due to the unique geographical environment, such as the long, gentle slopesm and large catchment area, the gully erosion process of the black soil region in Northeast China is obviously different from that of other areas. Gully erosion has mostly occurred on sloping farmland [8]. Ephemeral gullies and gullies are usually formed, which are characterized by a short length and small area. Mapping gullies in a wide range of temporal and spatial scales and identifying the formation mechanism and influencing factors are the key issues for soil conservation and land management in Northeast China.The conventional ground-based measurement methods, such as methods that utilize tape [9], microtopographic profilers [10], total stations [11], pins [12], and differential GPS [13], are time-consuming and labor-intensive when seeking high accuracy in field surveys. The ground-based monitoring methods are typically applicable at a small scale and for short-term monitoring [13, 14]. Unmanned aerial vehicles (UAV) with the photogrammetric technique of structure from motion (SFM) allow a higher level of detail and insights into the process of gully erosion [15, 16]. However, the flexible and promising method provides the user with limited continuous space coverage. At larger temporal and spatial scales, visual analysis or object-oriented analysis of high-resolution satellite remote sensing images have been used to quantify the temporal changes in various gully planar morphological parameters [17–19]. Of particular importance is the continuous mapping of gullies over large areas covering large regions. In recent years, many studies have used high-resolution satellite images to identify and extract the distribution of large-scale gullies and then to study the spatial distribution and dynamic changes in gullies. From 2010 to 2012, the National Administration of Surveying, Mapping and Geographic Information (NASG) conducted the first special census of gully erosion in the black soil region of Northeast China [20]. The survey was mainly based on the method of visual interpretation and field verification, making full use of high-resolution satellite images, airborne sensors, and ground data to obtain information on gullies. According to the field survey, the width of many gullies in the black soil region of Northeast China is less than 5 m. However, when the spatial resolution of satellite images is higher than 5 m, these gullies will be ignored, resulting in an increase in mapping errors. At this time, high spatial–temporal-resolution remote sensing images provide new possibilities for studying large-scale gully erosion [19, 21, 22].At present, researchers in China and elsewhere are carrying out relevant research on the influencing factors of gully erosion [23–26]. According to previous studies, the distribution of gullies is related to topography, soil, climate, precipitation, vegetation, land use, and human activities, etc. [1, 27]. With the rapid development of computer science, most studies use correlation, regression analysis, or machine learning models to identify the importance of gully erosion factors in order to assess gully erosion sensitivity [2, 3, 28–30]. However, these statistical methods cannot directly quantify the influence of the driving factors [31]. In addition, there is still a lack of studies that comprehensively compare factors and combinations of variables. Do these influencing factors operate independently or by interacting? The emergence of geodetectors provides new ideas and means to solve the above problems. A geodetector is a new tool for geographic research that can effectively analyze spatial differentiation in geographic phenomena and the factors that influence them [32]. It can not only quantitatively determine the dominant factors but can quantify the influence of two interacting explanatory variables on a specific target variable. Thus far, geodetectors have been used to analyze the driving forces and mechanisms of soil erosion [33, 34].It is still a major challenge to monitor and map the spatial differentiation of gully erosion at a large scale and quantitatively evaluate its driving mechanism in the black soil region of Northeast China. Thus, the purpose of this work is to (1) use the submeter satellite remote sensing image to perform the large-scale spatial mapping of gully erosion in the Sancha River basin located in the south of the typical black soil region of Northeast China by visual interpretation and field validation and (2) perform a quantitative attribution analysis of the spatial differentiation mechanism of gully erosion based on geodetection.


## 2. Materials and Methods

### 2.1. Study Area

Soil degradation caused by gully erosion has become an important problem of ecological restoration in the black soil region of Northeast China. The Sancha River basin, located in the south of Northeast China (Figure1(a)) and covering 163.83 km2, is a typical black soil region and suffers serious soil and water loss. This study area (Figures 1(b) and 1(c)), situated in the hilly area between the Changbai Mountains and Songliao Plain, has a mid-temperate continental monsoon climate with an annual temperature of 5.3°C; the highest average monthly temperature is 23.3°C and the lowest average monthly temperature is 16.3°C, and the soil types include black soil, meadow soil, albic soil, and dark brown soil. The annual precipitation ranges from 550 to 600 mm. The Sancha River basin has a mosaic of land cover comprising 66.2% farmland, 26.4% forest, 2.0% grassland, 2.9% water area, and 2.5% road and housing construction area. The land use type of this area is mainly slope farmland; the slope length of the slope farmland is mostly in the range of 300~500 m, where the main crops are soybean and corn. Due to the comprehensive influence of topographic factors, climatic factors, soil properties, land use types, and lithology, gully erosion is widely distributed in this area and has most occurred on the sloping farmland. Due to the substantial terrain undulation and long-term human activities, this region has little vegetation coverage and suffers from severe soil erosion by water. In particular, high-intensity reclamation and unreasonable land use have become important factors accelerating soil erosion in this region in the past 50 years.Figure 1
(a) The study area located in Northeast China. (b) DEM of the study area. (c) Pleiades-1B satellite image including the Sancha River catchment (27 April 2018).


### 2.2. Data Sources and Processing

To map gully erosion, Pleiades-1B images with a resolution of 0.5 m, collected from Google Earth, were drawn into a polygon according to the gully erosion area. The Pleiadia-1B image taken on 27 April 2018 shows a panchromatic band and four multi-spectral bands (blue, green, red, and near-infrared) (Figure1(c)). The data provider has made geometric, radiometric, and atmospheric corrections to the data.Gully erosion is a threshold-dependent process under the influence of a number of effective factors [35]. On the basis of referring to the relevant literature and considering the availability of data and the actual situation of the Sancha River basin, we finally selected 18 influencing factors. Elevation affects the types of vegetation in a region. Therefore, many researchers believe that altitude plays a crucial role in the study of gully erosion [36]. Slope shape and slope position can also cause the spatial differentiation of gully erosion. Slope aspect affects humidity, temperature conditions, and vegetation growth by affecting solar radiation. Slope degree is an important basic parameter to describe surface morphology and an important analysis factor in soil and water conservation. Slope length can control surface runoff velocity. The runoff and sediment yield increase with the increase in slope length. Plane curvature is an indicator of the turning of the ground, and sectional curvature can reflect the degree of terrain complexity. Runoff will be affected by slope shape [28], and plane curvature and sectional curvature become important indexes affecting gully erosion. The surface runoff and the resulting soil loss can be determined by the unit catchment area. It is the main factor causing the spatial differentiation of gully erosion. The topographic wetness index (TWI) is commonly employed as a proxy for the potential for surface and subsurface water accumulation due to runoff and lateral transmissivity [37]. TWI is considered to be an important factor affecting the development of gullies. Rainfall erosivity is positively correlated with soil and water loss, which can reflect the effect of rainfall on the triggering factors of gully erosion. The rate and pattern of gully development are largely controlled by soil type. Soil erodibility is influenced by the soil’s physical properties [38], which have an important impact on soil erosion and sediment yield. Lithological features are related to geomorphologic features and land surface characteristics [39]. The lithology of the geologic parent material influences the development of gullies [40]. In addition, unreasonable land use has become an important factor accelerating gully erosion in this region in the past 50 years. The distance from the river is directly related to runoff in the catchment area, so the distance from the river will affect gully erosion. Distance from the residential area and distance from the road represent the influence of human factors on gully erosion. Excessive utilization of land resources will promote land degradation and have a profound impact on gully erosion [41].The above factors were identified and divided into 6 categories: (1) topography factors, (2) climate and weather factors, (3) soil properties, (4) lithology, (5) land use, and (6) other factors (Figure2). Table 1 provides a more detailed overview of each variable. Topography factors, including elevation, slope shape, slope aspect, slope position, slope degree, slope length, plane curvature, sectional curvature, catchment area, and topographic wetness index, were mainly obtained by calculation or analysis in ArcGIS on the basis of a DEM of 5 m pixel size with a scale of 1 : 10,000. Climate and weather factors were represented by rainfall erosivity. Soil properties and lithology were represented by soil type, soil erodibility, and lithology. The soil type distribution map was provided by the Chinese Academy of Sciences. The lithology was obtained through 1 : 50,000 geological maps, provided by the Jilin Provincial Geological Database. Detailed algorithms and data on rainfall erosivity and soil erodibility are available from a previous study [42]. The land use data were derived from submeter-level images and ground data, as an important result of geographical condition monitoring. The other factors, including the distances from residential areas, rivers, and roads, were generated by using the buffer tool in ArcGIS. The other high-resolution satellite/aerial images, including Landsat8 (2013), Planet Labs (2014), Alos (2009), ZY3 (2013), Pleiades-1A (2013), Digital Mapping Camera (DMC, 2010), and Unmanned Air Vehicle (UAV, 2015), were provided by the Jilin Province Geomatics Center and the Chinese Academy of Sciences. The digital elevation model (DEM), land use, road, water, and residential sites were also provided by the Jilin Province Geomatics Center.Figure 2
Influencing factors on gully erosion ((a) rainfall erosivity, (b) elevation, (c) slope shape, (d) slope aspect, (e) slope position, (f) slope degree, (g) slope length, (h) distance from residential area, (i) distance from the river, (j) distance from the road, (k) plane curvature, (l) sectional curvature, (m) catchment area, (n) topographic wetness index, (o) soil type, (p) soil erodibility, (q) lithology, and (r) land use).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)Table 1
Information on factors affecting gully erosion.
FactorsDate typeScale/resolutionClimate and weather factorsRainfall erosivity (X1)Raster500 mTopography factorsElevation (X2)Raster5 mSlope shape (X3)Slope aspect (X4)Slope position (X5)Slope degree (X6)Slope length (X7)Plane curvature (X11)Sectional curvature (X12)Catchment area (X13)Topographic wetness index (X14)Soil propertiesSoil type (X15)Feature1 : 1,000,000Soil erodibility (X1616)Raster500 mLithologyLithology (X17)Feature1 : 50000Land useLand use (X18)Feature1 : 10000Other factorsDistance from residential area (X8)Raster10 mDistance from the river (X9)Distance from the road (X10)


### 2.3. Methods

#### 2.3.1. Gully Erosion Mapping Based on High-Resolution Satellite Imagery

The distribution information on gully erosion in the study area was mainly obtained through visual interpretation (Table2), combined with field investigation and verification, to judge the gully from the image and then obtain its data, including linear data and area data. The distribution of gully erosion was analyzed in terms of gully density and gully intensity. As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area. The distribution of gully density was estimated by using the line density tool in ArcGIS. The density was calculated as the ratio of the total length of the gully within the circular kernel (50 m search radius) and the total area of the circular kernel. Gully intensity refers to the total length of gully erosion per unit area, reflecting the degree of fragmentation of the surface and the degree of soil erosion.Table 2
Interpretation signs on Pleiades images.
LocationPleiades imagesGround photos126.137°E44.213°N126.261°E44.301°NAccording to the field validation result, the gully data were modified. The omitted and committed gullies were included and excluded, respectively. Then, the final gully distribution data were obtained. The gully density distribution was estimated by using the Focal Statistics tool in ArcGIS. The specific steps were as follows: a fishnet polygon containing cells of50m×50m was created, and then, the total length per area of each cell was calculated and a value was assigned to the corresponding cell; (2) the polygon was converted into a grid with a 5 m pixel size; and (3) focal Statistics was applied to calculate the average value of each input cell within a rectangular neighborhood of 100m×100m; thus, the distribution map of gully density was obtained. Similarly, gully intensity was calculated by using the measurement of total gully area per area.


#### 2.3.2. Geodetector Method

The geodetector method was established by Wang et al. and has been used extensively, mainly in factor detection, risk detection, interactive detection, and ecological detection. For more details about the geographical detector model, please see [43]. Briefly, the process is as follows:
(1)
The factor detector uses theq value to assess the impact of risk factors on the spatial pattern of gully erosion. A higher q value means that the risk factor has a stronger contribution to the spatial differentiation of gully erosion. It uses F-tests to compare whether the accumulated variance of each subregion is significantly different from the variance of the entire study region(1)q=1−∑h=1LNhσh2Nσ2,where h=1,⋯,L is the layer of independent variable X, Nh and N are the number of sample units in layer h and the total region, respectively, and σh2 and σ2 are the variance in the h layer and the variance in the region. The q value lies in [0, 1]. If factor X completely controls the soil erosion, the q value equals 1; if the factor X is completely unrelated to Y, the q value equals 0
(2)
The interaction detector compares the comprehensive contributions of two separate risk factors to the spatial differentiation of gully erosion, as well as their independent contributions. By doing so, it assesses whether the two risk factors weaken or reinforce each other, or whether they independently affect the spatial differentiation of gullies (Table3)(3)
Risk area detection can evaluate the differences in different driving factors in different areas of the study area and can be used to identify the distribution of gully erosion in different areas, which can be tested byt statistics(4)
The ecological detector evaluates whether the two risk factors are significantly different in the distribution and development of gully erosion. It also usesF-tests to compare the variance calculated in a subregion attributed to one risk factor with the variance attributed to another risk factorTable 3
Types of interaction between two covariates.
CriterionInteractionqX1∩X2<MinqX1,qX2Weakened, nonlinearMinqX1,qX2<qX1∩X2<MaxqX1,qX2Weakened, single factor nonlinearqX1∩X2>MaxqX1,qX2Enhanced, double factorsqX1∩X2=qX1+qX2IndependentqX1∩X2>qX1+qX2Enhanced, nonlinearDiscretization is aimed at transforming continuous data into discrete data. Compared with continuous data, discrete data are easier to understand, use, and explain and are closer to a knowledge-level representation [44]. Data discretization is the process whereby continuous data are divided into several intervals with selected cut points, where each interval is mapped to a qualitative symbol. A cut point is a value from the adjacent continuous data that divides them into two intervals. In actual applications, researchers always discretize continuous data with user-defined discretization and select the cut points according to their experience [45]. In this study, soil erodibility and distance from the residential area are divided into 4 grades. Soil type and slope shape are divided into 5 grades. Slope position, slope length, elevation, plane curvature, section curvature, distance from the road, and terrain wetness index are divided into 6 grades. Rainfall erosivity is divided into 8 grades. Regarding slope aspect, the catchment area is divided into 9 grades. Land use and distance from rivers are divided into 10 grades. Lithology is divided into 18 grades. In the process of division, while the soil type, slope shape, slope aspect, and land use are divided by a fixed cut point, the other factors are divided by the natural breaks (NB) method. The Fishnet tool in ArcGIS was used to extract the raster data to points. The sampling interval was set to 74 m. A total of 29876 points were extracted and used as the operating data for the geodetector.


## 2.1. Study Area

Soil degradation caused by gully erosion has become an important problem of ecological restoration in the black soil region of Northeast China. The Sancha River basin, located in the south of Northeast China (Figure1(a)) and covering 163.83 km2, is a typical black soil region and suffers serious soil and water loss. This study area (Figures 1(b) and 1(c)), situated in the hilly area between the Changbai Mountains and Songliao Plain, has a mid-temperate continental monsoon climate with an annual temperature of 5.3°C; the highest average monthly temperature is 23.3°C and the lowest average monthly temperature is 16.3°C, and the soil types include black soil, meadow soil, albic soil, and dark brown soil. The annual precipitation ranges from 550 to 600 mm. The Sancha River basin has a mosaic of land cover comprising 66.2% farmland, 26.4% forest, 2.0% grassland, 2.9% water area, and 2.5% road and housing construction area. The land use type of this area is mainly slope farmland; the slope length of the slope farmland is mostly in the range of 300~500 m, where the main crops are soybean and corn. Due to the comprehensive influence of topographic factors, climatic factors, soil properties, land use types, and lithology, gully erosion is widely distributed in this area and has most occurred on the sloping farmland. Due to the substantial terrain undulation and long-term human activities, this region has little vegetation coverage and suffers from severe soil erosion by water. In particular, high-intensity reclamation and unreasonable land use have become important factors accelerating soil erosion in this region in the past 50 years.Figure 1
(a) The study area located in Northeast China. (b) DEM of the study area. (c) Pleiades-1B satellite image including the Sancha River catchment (27 April 2018).


## 2.2. Data Sources and Processing

To map gully erosion, Pleiades-1B images with a resolution of 0.5 m, collected from Google Earth, were drawn into a polygon according to the gully erosion area. The Pleiadia-1B image taken on 27 April 2018 shows a panchromatic band and four multi-spectral bands (blue, green, red, and near-infrared) (Figure1(c)). The data provider has made geometric, radiometric, and atmospheric corrections to the data.Gully erosion is a threshold-dependent process under the influence of a number of effective factors [35]. On the basis of referring to the relevant literature and considering the availability of data and the actual situation of the Sancha River basin, we finally selected 18 influencing factors. Elevation affects the types of vegetation in a region. Therefore, many researchers believe that altitude plays a crucial role in the study of gully erosion [36]. Slope shape and slope position can also cause the spatial differentiation of gully erosion. Slope aspect affects humidity, temperature conditions, and vegetation growth by affecting solar radiation. Slope degree is an important basic parameter to describe surface morphology and an important analysis factor in soil and water conservation. Slope length can control surface runoff velocity. The runoff and sediment yield increase with the increase in slope length. Plane curvature is an indicator of the turning of the ground, and sectional curvature can reflect the degree of terrain complexity. Runoff will be affected by slope shape [28], and plane curvature and sectional curvature become important indexes affecting gully erosion. The surface runoff and the resulting soil loss can be determined by the unit catchment area. It is the main factor causing the spatial differentiation of gully erosion. The topographic wetness index (TWI) is commonly employed as a proxy for the potential for surface and subsurface water accumulation due to runoff and lateral transmissivity [37]. TWI is considered to be an important factor affecting the development of gullies. Rainfall erosivity is positively correlated with soil and water loss, which can reflect the effect of rainfall on the triggering factors of gully erosion. The rate and pattern of gully development are largely controlled by soil type. Soil erodibility is influenced by the soil’s physical properties [38], which have an important impact on soil erosion and sediment yield. Lithological features are related to geomorphologic features and land surface characteristics [39]. The lithology of the geologic parent material influences the development of gullies [40]. In addition, unreasonable land use has become an important factor accelerating gully erosion in this region in the past 50 years. The distance from the river is directly related to runoff in the catchment area, so the distance from the river will affect gully erosion. Distance from the residential area and distance from the road represent the influence of human factors on gully erosion. Excessive utilization of land resources will promote land degradation and have a profound impact on gully erosion [41].The above factors were identified and divided into 6 categories: (1) topography factors, (2) climate and weather factors, (3) soil properties, (4) lithology, (5) land use, and (6) other factors (Figure2). Table 1 provides a more detailed overview of each variable. Topography factors, including elevation, slope shape, slope aspect, slope position, slope degree, slope length, plane curvature, sectional curvature, catchment area, and topographic wetness index, were mainly obtained by calculation or analysis in ArcGIS on the basis of a DEM of 5 m pixel size with a scale of 1 : 10,000. Climate and weather factors were represented by rainfall erosivity. Soil properties and lithology were represented by soil type, soil erodibility, and lithology. The soil type distribution map was provided by the Chinese Academy of Sciences. The lithology was obtained through 1 : 50,000 geological maps, provided by the Jilin Provincial Geological Database. Detailed algorithms and data on rainfall erosivity and soil erodibility are available from a previous study [42]. The land use data were derived from submeter-level images and ground data, as an important result of geographical condition monitoring. The other factors, including the distances from residential areas, rivers, and roads, were generated by using the buffer tool in ArcGIS. The other high-resolution satellite/aerial images, including Landsat8 (2013), Planet Labs (2014), Alos (2009), ZY3 (2013), Pleiades-1A (2013), Digital Mapping Camera (DMC, 2010), and Unmanned Air Vehicle (UAV, 2015), were provided by the Jilin Province Geomatics Center and the Chinese Academy of Sciences. The digital elevation model (DEM), land use, road, water, and residential sites were also provided by the Jilin Province Geomatics Center.Figure 2
Influencing factors on gully erosion ((a) rainfall erosivity, (b) elevation, (c) slope shape, (d) slope aspect, (e) slope position, (f) slope degree, (g) slope length, (h) distance from residential area, (i) distance from the river, (j) distance from the road, (k) plane curvature, (l) sectional curvature, (m) catchment area, (n) topographic wetness index, (o) soil type, (p) soil erodibility, (q) lithology, and (r) land use).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)Table 1
Information on factors affecting gully erosion.
FactorsDate typeScale/resolutionClimate and weather factorsRainfall erosivity (X1)Raster500 mTopography factorsElevation (X2)Raster5 mSlope shape (X3)Slope aspect (X4)Slope position (X5)Slope degree (X6)Slope length (X7)Plane curvature (X11)Sectional curvature (X12)Catchment area (X13)Topographic wetness index (X14)Soil propertiesSoil type (X15)Feature1 : 1,000,000Soil erodibility (X1616)Raster500 mLithologyLithology (X17)Feature1 : 50000Land useLand use (X18)Feature1 : 10000Other factorsDistance from residential area (X8)Raster10 mDistance from the river (X9)Distance from the road (X10)


## 2.3. Methods

### 2.3.1. Gully Erosion Mapping Based on High-Resolution Satellite Imagery

The distribution information on gully erosion in the study area was mainly obtained through visual interpretation (Table2), combined with field investigation and verification, to judge the gully from the image and then obtain its data, including linear data and area data. The distribution of gully erosion was analyzed in terms of gully density and gully intensity. As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area. The distribution of gully density was estimated by using the line density tool in ArcGIS. The density was calculated as the ratio of the total length of the gully within the circular kernel (50 m search radius) and the total area of the circular kernel. Gully intensity refers to the total length of gully erosion per unit area, reflecting the degree of fragmentation of the surface and the degree of soil erosion.Table 2
Interpretation signs on Pleiades images.
LocationPleiades imagesGround photos126.137°E44.213°N126.261°E44.301°NAccording to the field validation result, the gully data were modified. The omitted and committed gullies were included and excluded, respectively. Then, the final gully distribution data were obtained. The gully density distribution was estimated by using the Focal Statistics tool in ArcGIS. The specific steps were as follows: a fishnet polygon containing cells of50m×50m was created, and then, the total length per area of each cell was calculated and a value was assigned to the corresponding cell; (2) the polygon was converted into a grid with a 5 m pixel size; and (3) focal Statistics was applied to calculate the average value of each input cell within a rectangular neighborhood of 100m×100m; thus, the distribution map of gully density was obtained. Similarly, gully intensity was calculated by using the measurement of total gully area per area.


### 2.3.2. Geodetector Method

The geodetector method was established by Wang et al. and has been used extensively, mainly in factor detection, risk detection, interactive detection, and ecological detection. For more details about the geographical detector model, please see [43]. Briefly, the process is as follows:
(1)
The factor detector uses theq value to assess the impact of risk factors on the spatial pattern of gully erosion. A higher q value means that the risk factor has a stronger contribution to the spatial differentiation of gully erosion. It uses F-tests to compare whether the accumulated variance of each subregion is significantly different from the variance of the entire study region(1)q=1−∑h=1LNhσh2Nσ2,where h=1,⋯,L is the layer of independent variable X, Nh and N are the number of sample units in layer h and the total region, respectively, and σh2 and σ2 are the variance in the h layer and the variance in the region. The q value lies in [0, 1]. If factor X completely controls the soil erosion, the q value equals 1; if the factor X is completely unrelated to Y, the q value equals 0
(2)
The interaction detector compares the comprehensive contributions of two separate risk factors to the spatial differentiation of gully erosion, as well as their independent contributions. By doing so, it assesses whether the two risk factors weaken or reinforce each other, or whether they independently affect the spatial differentiation of gullies (Table3)(3)
Risk area detection can evaluate the differences in different driving factors in different areas of the study area and can be used to identify the distribution of gully erosion in different areas, which can be tested byt statistics(4)
The ecological detector evaluates whether the two risk factors are significantly different in the distribution and development of gully erosion. It also usesF-tests to compare the variance calculated in a subregion attributed to one risk factor with the variance attributed to another risk factorTable 3
Types of interaction between two covariates.
CriterionInteractionqX1∩X2<MinqX1,qX2Weakened, nonlinearMinqX1,qX2<qX1∩X2<MaxqX1,qX2Weakened, single factor nonlinearqX1∩X2>MaxqX1,qX2Enhanced, double factorsqX1∩X2=qX1+qX2IndependentqX1∩X2>qX1+qX2Enhanced, nonlinearDiscretization is aimed at transforming continuous data into discrete data. Compared with continuous data, discrete data are easier to understand, use, and explain and are closer to a knowledge-level representation [44]. Data discretization is the process whereby continuous data are divided into several intervals with selected cut points, where each interval is mapped to a qualitative symbol. A cut point is a value from the adjacent continuous data that divides them into two intervals. In actual applications, researchers always discretize continuous data with user-defined discretization and select the cut points according to their experience [45]. In this study, soil erodibility and distance from the residential area are divided into 4 grades. Soil type and slope shape are divided into 5 grades. Slope position, slope length, elevation, plane curvature, section curvature, distance from the road, and terrain wetness index are divided into 6 grades. Rainfall erosivity is divided into 8 grades. Regarding slope aspect, the catchment area is divided into 9 grades. Land use and distance from rivers are divided into 10 grades. Lithology is divided into 18 grades. In the process of division, while the soil type, slope shape, slope aspect, and land use are divided by a fixed cut point, the other factors are divided by the natural breaks (NB) method. The Fishnet tool in ArcGIS was used to extract the raster data to points. The sampling interval was set to 74 m. A total of 29876 points were extracted and used as the operating data for the geodetector.


## 2.3.1. Gully Erosion Mapping Based on High-Resolution Satellite Imagery

The distribution information on gully erosion in the study area was mainly obtained through visual interpretation (Table2), combined with field investigation and verification, to judge the gully from the image and then obtain its data, including linear data and area data. The distribution of gully erosion was analyzed in terms of gully density and gully intensity. As the main indicator for assessing gully erosion, gully density within the study area was defined as the total gully length per area of the whole study area. The distribution of gully density was estimated by using the line density tool in ArcGIS. The density was calculated as the ratio of the total length of the gully within the circular kernel (50 m search radius) and the total area of the circular kernel. Gully intensity refers to the total length of gully erosion per unit area, reflecting the degree of fragmentation of the surface and the degree of soil erosion.Table 2
Interpretation signs on Pleiades images.
LocationPleiades imagesGround photos126.137°E44.213°N126.261°E44.301°NAccording to the field validation result, the gully data were modified. The omitted and committed gullies were included and excluded, respectively. Then, the final gully distribution data were obtained. The gully density distribution was estimated by using the Focal Statistics tool in ArcGIS. The specific steps were as follows: a fishnet polygon containing cells of50m×50m was created, and then, the total length per area of each cell was calculated and a value was assigned to the corresponding cell; (2) the polygon was converted into a grid with a 5 m pixel size; and (3) focal Statistics was applied to calculate the average value of each input cell within a rectangular neighborhood of 100m×100m; thus, the distribution map of gully density was obtained. Similarly, gully intensity was calculated by using the measurement of total gully area per area.


## 2.3.2. Geodetector Method

The geodetector method was established by Wang et al. and has been used extensively, mainly in factor detection, risk detection, interactive detection, and ecological detection. For more details about the geographical detector model, please see [43]. Briefly, the process is as follows:
(1)
The factor detector uses theq value to assess the impact of risk factors on the spatial pattern of gully erosion. A higher q value means that the risk factor has a stronger contribution to the spatial differentiation of gully erosion. It uses F-tests to compare whether the accumulated variance of each subregion is significantly different from the variance of the entire study region(1)q=1−∑h=1LNhσh2Nσ2,where h=1,⋯,L is the layer of independent variable X, Nh and N are the number of sample units in layer h and the total region, respectively, and σh2 and σ2 are the variance in the h layer and the variance in the region. The q value lies in [0, 1]. If factor X completely controls the soil erosion, the q value equals 1; if the factor X is completely unrelated to Y, the q value equals 0
(2)
The interaction detector compares the comprehensive contributions of two separate risk factors to the spatial differentiation of gully erosion, as well as their independent contributions. By doing so, it assesses whether the two risk factors weaken or reinforce each other, or whether they independently affect the spatial differentiation of gullies (Table3)(3)
Risk area detection can evaluate the differences in different driving factors in different areas of the study area and can be used to identify the distribution of gully erosion in different areas, which can be tested byt statistics(4)
The ecological detector evaluates whether the two risk factors are significantly different in the distribution and development of gully erosion. It also usesF-tests to compare the variance calculated in a subregion attributed to one risk factor with the variance attributed to another risk factorTable 3
Types of interaction between two covariates.
CriterionInteractionqX1∩X2<MinqX1,qX2Weakened, nonlinearMinqX1,qX2<qX1∩X2<MaxqX1,qX2Weakened, single factor nonlinearqX1∩X2>MaxqX1,qX2Enhanced, double factorsqX1∩X2=qX1+qX2IndependentqX1∩X2>qX1+qX2Enhanced, nonlinearDiscretization is aimed at transforming continuous data into discrete data. Compared with continuous data, discrete data are easier to understand, use, and explain and are closer to a knowledge-level representation [44]. Data discretization is the process whereby continuous data are divided into several intervals with selected cut points, where each interval is mapped to a qualitative symbol. A cut point is a value from the adjacent continuous data that divides them into two intervals. In actual applications, researchers always discretize continuous data with user-defined discretization and select the cut points according to their experience [45]. In this study, soil erodibility and distance from the residential area are divided into 4 grades. Soil type and slope shape are divided into 5 grades. Slope position, slope length, elevation, plane curvature, section curvature, distance from the road, and terrain wetness index are divided into 6 grades. Rainfall erosivity is divided into 8 grades. Regarding slope aspect, the catchment area is divided into 9 grades. Land use and distance from rivers are divided into 10 grades. Lithology is divided into 18 grades. In the process of division, while the soil type, slope shape, slope aspect, and land use are divided by a fixed cut point, the other factors are divided by the natural breaks (NB) method. The Fishnet tool in ArcGIS was used to extract the raster data to points. The sampling interval was set to 74 m. A total of 29876 points were extracted and used as the operating data for the geodetector.


## 3. Results

### 3.1. Analysis of the Spatial Variation in Gully Erosion

Remote sensing data need to be verified by a field investigation to ensure the accuracy of gully erosion mapping and correct errors. In the field validation work, 20 gullies were observed and recorded at the Jilin Continuous Operation Reference Station (JLORS), installed with Trimble equipment. The indoor study showed that 18 gullies were captured correctly, while a roadside gully was omitted and another field path was mistaken for an ephemeral gully (Table4). Nevertheless, the overall accuracy of gully interpretation reached 90%, and the very-high-resolution (submeter) satellite image showed a strong capacity for the detection of various gully types.Table 4
Examples of validation results.
LocationPleiades imagesGround photosValidation results126.155°E44.264°N√126.149°E44.357°N√126.134°E44.217°N×Based on the high-resolution images, gullies were extracted by visual interpretation and field investigation (Figure3). There were 611 gullies in 2018, with a total area of 128.31 km and 0.86 km2. According to the distribution map of gully density and intensity in 2018, it could be seen that the area with higher gully density was concentrated in the northeast. The overall distribution was more uniform, and the highest density reached 10.55 km/km2. The distribution area of gully intensity was concentrated in the northeast and south, with less in the north, and the highest reached 158457 m2/km2. The areas with high gully density and intensity in the study area were all in the northeast, which may have been due to the high soil erodibility in the northeast, where soil is more prone to erosion. The gully density of the whole area was 505.98 m/km2 and the average gully intensity was 3412.29 m2/km2.Figure 3
The distribution map of (a) gullies, (b) gully density, and (c) gully intensity in the study area.
(a)(b)(c)


### 3.2. Geodetector-Based Quantitative Attribution Analysis of Gully Erosion

#### 3.2.1. Significance Analysis of Factors Affecting Gully Erosion

The spatial differentiation of gully density and intensity in 2018 was attributed to the geodetector mode, and the detection results are shown in Table5. The ability of different factors to explain the spatial distribution of gully density and intensity was as follows: gully density: lithology (X17) > soil type (X15) > rainfall erosivity (X1) > elevation (X2) > soil erodibility (X16) > land use (X18) > slope degree (X6) > distance from the river (X9) > topographic wetness index (X14) > catchment area (X13) > distance from residential area (X8) > distance from the road (X10) > slope shape (X3) > slope aspect (X4) > slope length (X7) > slope position (X5) > plane curvature (X11) > section curvature (X12). Gully intensity: lithology (X17) > soil type (X15) > rainfall erosivity (X1) > elevation (X2) > soil erodibility (X16) > catchment area (X13) > slope degree (X6) > land use (X18) > distance from river (X9) > topographic wetness index (X14) > slope direction (X4) > slope shape (X3) > distance from road (X10) > slope position (X5) > distance from residential area (X8) > slope length (X7) > plane curvature (X11) > sectional curvature (X12). It can be seen that different factors have different explanatory power with regard to gully density and intensity, indicating that lithology and soil type are the main factors affecting gully erosion distribution. Among them, the explanatory power of lithology and soil type is the highest, the explanatory power of erosion density is 12.9% and 9.9%, respectively, and the explanatory power of erosion intensity is 12.7% and 11.1%, respectively, which is the main influencing factor. The explanatory power of rainfall erosivity, catchment area, elevation, soil erodibility, land use, slope, and distance from the road for gully density and intensity is more than 1%; thus, they are secondary influencing factors. The explanatory power of the terrain wetness index, settlement distance, distance from the road, slope aspect, slope length, slope position, slope shape, plane curvature, and profile curvature is less than 1%, and their influence is low.Table 5
Theq values of influencing factors in gully density and intensity.
Influencing factorX1X2X3X4X5X6X7X8X9X10X11X12X13X14X15X16X17X18Gully density0.0781940.0677570.0077110.0074290.0034100.0383310.0060360.0099340.0242700.0092500.0008880.0007350.0160670.0170620.0990090.0478050.1287440.044721Gully intensity0.0712880.0401100.0044320.0056180.0038610.0246710.0032040.0035090.0123290.0040220.0007060.0006290.0248980.0100670.1108130.0351740.1270370.016478


#### 3.2.2. Analysis of Interactions between Factors Affecting Gully Erosion

When most of the factors interact, their explanatory power of each one is enhanced. The main conclusion is that the explanatory power of the interaction between two factors is higher than that of a single factor. Among them, the dominant factors represented by lithology are more obvious. The following table considers the interactions among dominant factors in the spatial distribution of gully density and intensity (Tables6 and 7).Table 6
The dominant interactions between two covariates in gully density.
InteractionTwo-factorq valueThe sum of theq values of two factorsResultExplanationX17∩X10.217>0.207=X17∩X1C>A+BNonlinear enhancedX17∩X20.214>0.197=X17∩X2C>A+BNonlinear enhancedX17∩X30.131>0.129=X17∩X3C>A+BNonlinear enhancedX17∩X40.154>0.136=X17∩X4C>A+BNonlinear enhancedX17∩X50.145>0.132=X17∩X5C>A+BNonlinear enhancedX17∩X60.172>0.167=X17∩X6C>A+BNonlinear enhancedX17∩X70.136>0.134=X17∩X7C>A+BNonlinear enhancedX17∩X80.154>0.139=X17∩X8C>A+BNonlinear enhancedX17∩X90.172>0.153=X17∩X9C>A+BNonlinear enhancedX17∩X100.188>0.138=X17∩X10C>A+BNonlinear enhancedX17∩X110.134>0.130=X17∩X12C>A+BNonlinear enhancedX17∩X120.136>0.129=X17∩X12C>A+BNonlinear enhancedX17∩X130.224>0.203=X17∩X13C>A+BNonlinear enhancedX17∩X140.150>0.146=X17∩X14C>A+BNonlinear enhancedX17∩X150.287>0.228=X17∩X15C>A+BNonlinear enhancedX17∩X160.229>0.177=X17∩X16C>A+BNonlinear enhancedX17∩X180.211>0.173=X17∩X18C>A+BNonlinear enhancedTable 7
The dominant interactions between two covariates in gully intensity.
InteractionTwo-factorq valueThe sum of theq values of two factorsResultExplanationX17∩X10.203>0.198=X17∩X1C>A+BNonlinear enhancedX17∩X20.185>0.167=X17∩X2C>A+BNonlinear enhancedX17∩X30.128>0.127=X17∩X3C>A+BNonlinear enhancedX17∩X40.151>0.133=X17∩X4C>A+BNonlinear enhancedX17∩X50.146>0.131=X17∩X5C>A+BNonlinear enhancedX17∩X60.152=0.152=X17∩X6C=A+BIndependent of each otherX17∩X70.130=0.130=X17∩X7C=A+BIndependent of each otherX17∩X80.139>0.131=X17∩X8C>A+BNonlinear enhancedX17∩X90.156>0.139=X17∩X9C>A+BNonlinear enhancedX17∩X100.161>0.131=X17∩X10C>A+BNonlinear enhancedX17∩X110.130>0.128=X17∩X11C>A+BNonlinear enhancedX17∩X120.131>0.128=X17∩X12C>A+BNonlinear enhancedX17∩X130.241>0.202=X17∩X13C>A+BNonlinear enhancedX17∩X140.138>0.137=X17∩X14C>A+BNonlinear enhancedX17∩X150.325>0.238=X17∩X15C>A+BNonlinear enhancedX17∩X160.249>0.162=X17∩X16C>A+BNonlinear enhancedX17∩X180.178>0.144=X17∩X18C>A+BNonlinear enhancedThe following describes the interaction relationship between the dominant factors in the spatial distribution of gully density:X17∩X150.287>X17∩X160.229>X17∩X130.224>X17∩X10.217>X17∩X20.214>X17∩X180.211>X17∩X100.188>X17∩X90.1722>X17∩X60.1715>X17∩X40.1542>X17∩X80.1538>X17∩X140.150>X17∩X50.145>X17∩X120.1364>X17∩X70.1356>X17∩X110.134>X17∩X30.131. The following describes the interactive relationship between the dominant factors in the spatial distribution of gully intensity: X17∩X150.325>X17∩X160.249>X17∩X130.241>X17∩X10.203>X17∩X20.185>X17∩X180.178>X17∩X100.161>X17∩X90.156>X17∩X60.152>X17∩X40.151>X17∩X50.146>X17∩X80.139>X17∩X140.138X17∩X120.131>X17∩X110.1303>X17∩X70.1297>X17∩X30.128.The results show that the role of the dominant factors in the spatial distribution of gully density and intensity is roughly the same. There is a nonlinear, enhanced relationship between lithology as the dominant factor and most other factors. Regarding the spatial distribution of gully intensity, only two pairs of factors are independent: the interaction between slope and lithology, slope length, and lithology.


#### 3.2.3. Ecological Detection and Analysis of Factors

Ecological exploration can reflect whether there is a significant difference in the influence of various driving factors on the development of gully erosion. If there is a significant difference between the two factors, we mark it as “Y.” If there is no significant difference, it is marked “N.” Among the rainfall factors, there is no significant relationship between rainfall and topography and land use; there is only an interaction with soil and lithology. Among the topographic factors, most of them are significantly related to soil, lithology, and land use, but not to rainfall erosion. Among soil and lithologic factors, all factors are significantly related to other factors, among which lithology is particularly prominent. There is a significant relationship between land use and most topographic factors, but not with rainfall, soil, and lithology factors.


#### 3.2.4. Identification of High-Risk Areas of Gully Erosion

According to the risk detection in the geodetector mode, we can obtain the distribution characteristics of gullies and the high-risk areas of gullies. Moreover, we can further judge whether there are significant differences in the amount of erosion between different levels of influencing factors (Table8). Among them, rainfall erosivity, elevation, slope shape, slope aspect, slope position, slope length, slope, residential distance, distance from river, distance from road, catchment area, terrain wetness index, soil type, soil erodibility, lithology, and land use display significant differences among different levels. Among the rainfall factors, the average value of gully density and intensity is the highest when the rainfall erosivity is between 210 and 211. Among the topographic factors, when the elevation is 220-260 m and the slope is 4-8°, the average value of gully density and intensity is the highest, indicating that the risk of gully erosion is higher on flat land. The lithology belongs to gray–white matter rhyolite, spherulite rhyolite, and crystal clastic tuff, with the highest gully density and intensity.Table 8
High-risk areas of gully erosion and its mean value.
FactorsHigh-risk area of gully densityAverage gully density (m/km2)High-risk area of gully intensityAverage gully intensity (m2/km2)Rainfall erosivity210-2111281.461210-21111752.866Elevation220-260 m1197.326220-260 m7991.692Slope shapeConvex slope856.4095653Convex slope5781.986Slope directionWest slope948.096West slope6605.746Slope positionMidslope932.913Midslope6706.873Slope degree4-8°1218.9204-8°8326.412Slope length0-150 m834.0170-150 m5588.848Residential distance<0.5879.799<0.55806.133Distance from the river250 m-500 m955.356250 m-500 m6701.476Distance from the road100 m-200 m906.231200 m-300 m6452.205Plane curvature0-5857.6840-55802.145Sectional curvature-25–-5880.6405-155990.171Catchment area2-3 km21223.7242-3 km210554.147Topographic wetness index5-7.5939.9085-7.56449.666Soil typeHumus dark brown soil1923.932Humus dark brown soil17149.954Soil erodibility0.27-0.291250.7450.27-0.299139.993LithologyGray–white matter rhyolite, globular rhyolite, crystal clastic tuff3542.989Gray–white matter rhyolite, globular rhyolite, crystal clastic tuff22833.733Land useCultivated land1027.192Cultivated land6656.717According to risk detection, on the whole, the soil and lithology are the main factors affecting the distribution of gullies. Regarding land use types, the main influence on the distribution of gully erosion is cultivated land. This may be due to the shortage of woodland caused by residential farming, which further leads to soil erosion. Thus, it can be seen that returning farmland to forests may slow down the development trend of gullies and improve the erosion conditions.


## 3.1. Analysis of the Spatial Variation in Gully Erosion

Remote sensing data need to be verified by a field investigation to ensure the accuracy of gully erosion mapping and correct errors. In the field validation work, 20 gullies were observed and recorded at the Jilin Continuous Operation Reference Station (JLORS), installed with Trimble equipment. The indoor study showed that 18 gullies were captured correctly, while a roadside gully was omitted and another field path was mistaken for an ephemeral gully (Table4). Nevertheless, the overall accuracy of gully interpretation reached 90%, and the very-high-resolution (submeter) satellite image showed a strong capacity for the detection of various gully types.Table 4
Examples of validation results.
LocationPleiades imagesGround photosValidation results126.155°E44.264°N√126.149°E44.357°N√126.134°E44.217°N×Based on the high-resolution images, gullies were extracted by visual interpretation and field investigation (Figure3). There were 611 gullies in 2018, with a total area of 128.31 km and 0.86 km2. According to the distribution map of gully density and intensity in 2018, it could be seen that the area with higher gully density was concentrated in the northeast. The overall distribution was more uniform, and the highest density reached 10.55 km/km2. The distribution area of gully intensity was concentrated in the northeast and south, with less in the north, and the highest reached 158457 m2/km2. The areas with high gully density and intensity in the study area were all in the northeast, which may have been due to the high soil erodibility in the northeast, where soil is more prone to erosion. The gully density of the whole area was 505.98 m/km2 and the average gully intensity was 3412.29 m2/km2.Figure 3
The distribution map of (a) gullies, (b) gully density, and (c) gully intensity in the study area.
(a)(b)(c)


## 3.2. Geodetector-Based Quantitative Attribution Analysis of Gully Erosion

### 3.2.1. Significance Analysis of Factors Affecting Gully Erosion

The spatial differentiation of gully density and intensity in 2018 was attributed to the geodetector mode, and the detection results are shown in Table5. The ability of different factors to explain the spatial distribution of gully density and intensity was as follows: gully density: lithology (X17) > soil type (X15) > rainfall erosivity (X1) > elevation (X2) > soil erodibility (X16) > land use (X18) > slope degree (X6) > distance from the river (X9) > topographic wetness index (X14) > catchment area (X13) > distance from residential area (X8) > distance from the road (X10) > slope shape (X3) > slope aspect (X4) > slope length (X7) > slope position (X5) > plane curvature (X11) > section curvature (X12). Gully intensity: lithology (X17) > soil type (X15) > rainfall erosivity (X1) > elevation (X2) > soil erodibility (X16) > catchment area (X13) > slope degree (X6) > land use (X18) > distance from river (X9) > topographic wetness index (X14) > slope direction (X4) > slope shape (X3) > distance from road (X10) > slope position (X5) > distance from residential area (X8) > slope length (X7) > plane curvature (X11) > sectional curvature (X12). It can be seen that different factors have different explanatory power with regard to gully density and intensity, indicating that lithology and soil type are the main factors affecting gully erosion distribution. Among them, the explanatory power of lithology and soil type is the highest, the explanatory power of erosion density is 12.9% and 9.9%, respectively, and the explanatory power of erosion intensity is 12.7% and 11.1%, respectively, which is the main influencing factor. The explanatory power of rainfall erosivity, catchment area, elevation, soil erodibility, land use, slope, and distance from the road for gully density and intensity is more than 1%; thus, they are secondary influencing factors. The explanatory power of the terrain wetness index, settlement distance, distance from the road, slope aspect, slope length, slope position, slope shape, plane curvature, and profile curvature is less than 1%, and their influence is low.Table 5
Theq values of influencing factors in gully density and intensity.
Influencing factorX1X2X3X4X5X6X7X8X9X10X11X12X13X14X15X16X17X18Gully density0.0781940.0677570.0077110.0074290.0034100.0383310.0060360.0099340.0242700.0092500.0008880.0007350.0160670.0170620.0990090.0478050.1287440.044721Gully intensity0.0712880.0401100.0044320.0056180.0038610.0246710.0032040.0035090.0123290.0040220.0007060.0006290.0248980.0100670.1108130.0351740.1270370.016478


### 3.2.2. Analysis of Interactions between Factors Affecting Gully Erosion

When most of the factors interact, their explanatory power of each one is enhanced. The main conclusion is that the explanatory power of the interaction between two factors is higher than that of a single factor. Among them, the dominant factors represented by lithology are more obvious. The following table considers the interactions among dominant factors in the spatial distribution of gully density and intensity (Tables6 and 7).Table 6
The dominant interactions between two covariates in gully density.
InteractionTwo-factorq valueThe sum of theq values of two factorsResultExplanationX17∩X10.217>0.207=X17∩X1C>A+BNonlinear enhancedX17∩X20.214>0.197=X17∩X2C>A+BNonlinear enhancedX17∩X30.131>0.129=X17∩X3C>A+BNonlinear enhancedX17∩X40.154>0.136=X17∩X4C>A+BNonlinear enhancedX17∩X50.145>0.132=X17∩X5C>A+BNonlinear enhancedX17∩X60.172>0.167=X17∩X6C>A+BNonlinear enhancedX17∩X70.136>0.134=X17∩X7C>A+BNonlinear enhancedX17∩X80.154>0.139=X17∩X8C>A+BNonlinear enhancedX17∩X90.172>0.153=X17∩X9C>A+BNonlinear enhancedX17∩X100.188>0.138=X17∩X10C>A+BNonlinear enhancedX17∩X110.134>0.130=X17∩X12C>A+BNonlinear enhancedX17∩X120.136>0.129=X17∩X12C>A+BNonlinear enhancedX17∩X130.224>0.203=X17∩X13C>A+BNonlinear enhancedX17∩X140.150>0.146=X17∩X14C>A+BNonlinear enhancedX17∩X150.287>0.228=X17∩X15C>A+BNonlinear enhancedX17∩X160.229>0.177=X17∩X16C>A+BNonlinear enhancedX17∩X180.211>0.173=X17∩X18C>A+BNonlinear enhancedTable 7
The dominant interactions between two covariates in gully intensity.
InteractionTwo-factorq valueThe sum of theq values of two factorsResultExplanationX17∩X10.203>0.198=X17∩X1C>A+BNonlinear enhancedX17∩X20.185>0.167=X17∩X2C>A+BNonlinear enhancedX17∩X30.128>0.127=X17∩X3C>A+BNonlinear enhancedX17∩X40.151>0.133=X17∩X4C>A+BNonlinear enhancedX17∩X50.146>0.131=X17∩X5C>A+BNonlinear enhancedX17∩X60.152=0.152=X17∩X6C=A+BIndependent of each otherX17∩X70.130=0.130=X17∩X7C=A+BIndependent of each otherX17∩X80.139>0.131=X17∩X8C>A+BNonlinear enhancedX17∩X90.156>0.139=X17∩X9C>A+BNonlinear enhancedX17∩X100.161>0.131=X17∩X10C>A+BNonlinear enhancedX17∩X110.130>0.128=X17∩X11C>A+BNonlinear enhancedX17∩X120.131>0.128=X17∩X12C>A+BNonlinear enhancedX17∩X130.241>0.202=X17∩X13C>A+BNonlinear enhancedX17∩X140.138>0.137=X17∩X14C>A+BNonlinear enhancedX17∩X150.325>0.238=X17∩X15C>A+BNonlinear enhancedX17∩X160.249>0.162=X17∩X16C>A+BNonlinear enhancedX17∩X180.178>0.144=X17∩X18C>A+BNonlinear enhancedThe following describes the interaction relationship between the dominant factors in the spatial distribution of gully density:X17∩X150.287>X17∩X160.229>X17∩X130.224>X17∩X10.217>X17∩X20.214>X17∩X180.211>X17∩X100.188>X17∩X90.1722>X17∩X60.1715>X17∩X40.1542>X17∩X80.1538>X17∩X140.150>X17∩X50.145>X17∩X120.1364>X17∩X70.1356>X17∩X110.134>X17∩X30.131. The following describes the interactive relationship between the dominant factors in the spatial distribution of gully intensity: X17∩X150.325>X17∩X160.249>X17∩X130.241>X17∩X10.203>X17∩X20.185>X17∩X180.178>X17∩X100.161>X17∩X90.156>X17∩X60.152>X17∩X40.151>X17∩X50.146>X17∩X80.139>X17∩X140.138X17∩X120.131>X17∩X110.1303>X17∩X70.1297>X17∩X30.128.The results show that the role of the dominant factors in the spatial distribution of gully density and intensity is roughly the same. There is a nonlinear, enhanced relationship between lithology as the dominant factor and most other factors. Regarding the spatial distribution of gully intensity, only two pairs of factors are independent: the interaction between slope and lithology, slope length, and lithology.


### 3.2.3. Ecological Detection and Analysis of Factors

Ecological exploration can reflect whether there is a significant difference in the influence of various driving factors on the development of gully erosion. If there is a significant difference between the two factors, we mark it as “Y.” If there is no significant difference, it is marked “N.” Among the rainfall factors, there is no significant relationship between rainfall and topography and land use; there is only an interaction with soil and lithology. Among the topographic factors, most of them are significantly related to soil, lithology, and land use, but not to rainfall erosion. Among soil and lithologic factors, all factors are significantly related to other factors, among which lithology is particularly prominent. There is a significant relationship between land use and most topographic factors, but not with rainfall, soil, and lithology factors.


### 3.2.4. Identification of High-Risk Areas of Gully Erosion

According to the risk detection in the geodetector mode, we can obtain the distribution characteristics of gullies and the high-risk areas of gullies. Moreover, we can further judge whether there are significant differences in the amount of erosion between different levels of influencing factors (Table8). Among them, rainfall erosivity, elevation, slope shape, slope aspect, slope position, slope length, slope, residential distance, distance from river, distance from road, catchment area, terrain wetness index, soil type, soil erodibility, lithology, and land use display significant differences among different levels. Among the rainfall factors, the average value of gully density and intensity is the highest when the rainfall erosivity is between 210 and 211. Among the topographic factors, when the elevation is 220-260 m and the slope is 4-8°, the average value of gully density and intensity is the highest, indicating that the risk of gully erosion is higher on flat land. The lithology belongs to gray–white matter rhyolite, spherulite rhyolite, and crystal clastic tuff, with the highest gully density and intensity.Table 8
High-risk areas of gully erosion and its mean value.
FactorsHigh-risk area of gully densityAverage gully density (m/km2)High-risk area of gully intensityAverage gully intensity (m2/km2)Rainfall erosivity210-2111281.461210-21111752.866Elevation220-260 m1197.326220-260 m7991.692Slope shapeConvex slope856.4095653Convex slope5781.986Slope directionWest slope948.096West slope6605.746Slope positionMidslope932.913Midslope6706.873Slope degree4-8°1218.9204-8°8326.412Slope length0-150 m834.0170-150 m5588.848Residential distance<0.5879.799<0.55806.133Distance from the river250 m-500 m955.356250 m-500 m6701.476Distance from the road100 m-200 m906.231200 m-300 m6452.205Plane curvature0-5857.6840-55802.145Sectional curvature-25–-5880.6405-155990.171Catchment area2-3 km21223.7242-3 km210554.147Topographic wetness index5-7.5939.9085-7.56449.666Soil typeHumus dark brown soil1923.932Humus dark brown soil17149.954Soil erodibility0.27-0.291250.7450.27-0.299139.993LithologyGray–white matter rhyolite, globular rhyolite, crystal clastic tuff3542.989Gray–white matter rhyolite, globular rhyolite, crystal clastic tuff22833.733Land useCultivated land1027.192Cultivated land6656.717According to risk detection, on the whole, the soil and lithology are the main factors affecting the distribution of gullies. Regarding land use types, the main influence on the distribution of gully erosion is cultivated land. This may be due to the shortage of woodland caused by residential farming, which further leads to soil erosion. Thus, it can be seen that returning farmland to forests may slow down the development trend of gullies and improve the erosion conditions.


## 3.2.1. Significance Analysis of Factors Affecting Gully Erosion

The spatial differentiation of gully density and intensity in 2018 was attributed to the geodetector mode, and the detection results are shown in Table5. The ability of different factors to explain the spatial distribution of gully density and intensity was as follows: gully density: lithology (X17) > soil type (X15) > rainfall erosivity (X1) > elevation (X2) > soil erodibility (X16) > land use (X18) > slope degree (X6) > distance from the river (X9) > topographic wetness index (X14) > catchment area (X13) > distance from residential area (X8) > distance from the road (X10) > slope shape (X3) > slope aspect (X4) > slope length (X7) > slope position (X5) > plane curvature (X11) > section curvature (X12). Gully intensity: lithology (X17) > soil type (X15) > rainfall erosivity (X1) > elevation (X2) > soil erodibility (X16) > catchment area (X13) > slope degree (X6) > land use (X18) > distance from river (X9) > topographic wetness index (X14) > slope direction (X4) > slope shape (X3) > distance from road (X10) > slope position (X5) > distance from residential area (X8) > slope length (X7) > plane curvature (X11) > sectional curvature (X12). It can be seen that different factors have different explanatory power with regard to gully density and intensity, indicating that lithology and soil type are the main factors affecting gully erosion distribution. Among them, the explanatory power of lithology and soil type is the highest, the explanatory power of erosion density is 12.9% and 9.9%, respectively, and the explanatory power of erosion intensity is 12.7% and 11.1%, respectively, which is the main influencing factor. The explanatory power of rainfall erosivity, catchment area, elevation, soil erodibility, land use, slope, and distance from the road for gully density and intensity is more than 1%; thus, they are secondary influencing factors. The explanatory power of the terrain wetness index, settlement distance, distance from the road, slope aspect, slope length, slope position, slope shape, plane curvature, and profile curvature is less than 1%, and their influence is low.Table 5
Theq values of influencing factors in gully density and intensity.
Influencing factorX1X2X3X4X5X6X7X8X9X10X11X12X13X14X15X16X17X18Gully density0.0781940.0677570.0077110.0074290.0034100.0383310.0060360.0099340.0242700.0092500.0008880.0007350.0160670.0170620.0990090.0478050.1287440.044721Gully intensity0.0712880.0401100.0044320.0056180.0038610.0246710.0032040.0035090.0123290.0040220.0007060.0006290.0248980.0100670.1108130.0351740.1270370.016478


## 3.2.2. Analysis of Interactions between Factors Affecting Gully Erosion

When most of the factors interact, their explanatory power of each one is enhanced. The main conclusion is that the explanatory power of the interaction between two factors is higher than that of a single factor. Among them, the dominant factors represented by lithology are more obvious. The following table considers the interactions among dominant factors in the spatial distribution of gully density and intensity (Tables6 and 7).Table 6
The dominant interactions between two covariates in gully density.
InteractionTwo-factorq valueThe sum of theq values of two factorsResultExplanationX17∩X10.217>0.207=X17∩X1C>A+BNonlinear enhancedX17∩X20.214>0.197=X17∩X2C>A+BNonlinear enhancedX17∩X30.131>0.129=X17∩X3C>A+BNonlinear enhancedX17∩X40.154>0.136=X17∩X4C>A+BNonlinear enhancedX17∩X50.145>0.132=X17∩X5C>A+BNonlinear enhancedX17∩X60.172>0.167=X17∩X6C>A+BNonlinear enhancedX17∩X70.136>0.134=X17∩X7C>A+BNonlinear enhancedX17∩X80.154>0.139=X17∩X8C>A+BNonlinear enhancedX17∩X90.172>0.153=X17∩X9C>A+BNonlinear enhancedX17∩X100.188>0.138=X17∩X10C>A+BNonlinear enhancedX17∩X110.134>0.130=X17∩X12C>A+BNonlinear enhancedX17∩X120.136>0.129=X17∩X12C>A+BNonlinear enhancedX17∩X130.224>0.203=X17∩X13C>A+BNonlinear enhancedX17∩X140.150>0.146=X17∩X14C>A+BNonlinear enhancedX17∩X150.287>0.228=X17∩X15C>A+BNonlinear enhancedX17∩X160.229>0.177=X17∩X16C>A+BNonlinear enhancedX17∩X180.211>0.173=X17∩X18C>A+BNonlinear enhancedTable 7
The dominant interactions between two covariates in gully intensity.
InteractionTwo-factorq valueThe sum of theq values of two factorsResultExplanationX17∩X10.203>0.198=X17∩X1C>A+BNonlinear enhancedX17∩X20.185>0.167=X17∩X2C>A+BNonlinear enhancedX17∩X30.128>0.127=X17∩X3C>A+BNonlinear enhancedX17∩X40.151>0.133=X17∩X4C>A+BNonlinear enhancedX17∩X50.146>0.131=X17∩X5C>A+BNonlinear enhancedX17∩X60.152=0.152=X17∩X6C=A+BIndependent of each otherX17∩X70.130=0.130=X17∩X7C=A+BIndependent of each otherX17∩X80.139>0.131=X17∩X8C>A+BNonlinear enhancedX17∩X90.156>0.139=X17∩X9C>A+BNonlinear enhancedX17∩X100.161>0.131=X17∩X10C>A+BNonlinear enhancedX17∩X110.130>0.128=X17∩X11C>A+BNonlinear enhancedX17∩X120.131>0.128=X17∩X12C>A+BNonlinear enhancedX17∩X130.241>0.202=X17∩X13C>A+BNonlinear enhancedX17∩X140.138>0.137=X17∩X14C>A+BNonlinear enhancedX17∩X150.325>0.238=X17∩X15C>A+BNonlinear enhancedX17∩X160.249>0.162=X17∩X16C>A+BNonlinear enhancedX17∩X180.178>0.144=X17∩X18C>A+BNonlinear enhancedThe following describes the interaction relationship between the dominant factors in the spatial distribution of gully density:X17∩X150.287>X17∩X160.229>X17∩X130.224>X17∩X10.217>X17∩X20.214>X17∩X180.211>X17∩X100.188>X17∩X90.1722>X17∩X60.1715>X17∩X40.1542>X17∩X80.1538>X17∩X140.150>X17∩X50.145>X17∩X120.1364>X17∩X70.1356>X17∩X110.134>X17∩X30.131. The following describes the interactive relationship between the dominant factors in the spatial distribution of gully intensity: X17∩X150.325>X17∩X160.249>X17∩X130.241>X17∩X10.203>X17∩X20.185>X17∩X180.178>X17∩X100.161>X17∩X90.156>X17∩X60.152>X17∩X40.151>X17∩X50.146>X17∩X80.139>X17∩X140.138X17∩X120.131>X17∩X110.1303>X17∩X70.1297>X17∩X30.128.The results show that the role of the dominant factors in the spatial distribution of gully density and intensity is roughly the same. There is a nonlinear, enhanced relationship between lithology as the dominant factor and most other factors. Regarding the spatial distribution of gully intensity, only two pairs of factors are independent: the interaction between slope and lithology, slope length, and lithology.


## 3.2.3. Ecological Detection and Analysis of Factors

Ecological exploration can reflect whether there is a significant difference in the influence of various driving factors on the development of gully erosion. If there is a significant difference between the two factors, we mark it as “Y.” If there is no significant difference, it is marked “N.” Among the rainfall factors, there is no significant relationship between rainfall and topography and land use; there is only an interaction with soil and lithology. Among the topographic factors, most of them are significantly related to soil, lithology, and land use, but not to rainfall erosion. Among soil and lithologic factors, all factors are significantly related to other factors, among which lithology is particularly prominent. There is a significant relationship between land use and most topographic factors, but not with rainfall, soil, and lithology factors.


## 3.2.4. Identification of High-Risk Areas of Gully Erosion

According to the risk detection in the geodetector mode, we can obtain the distribution characteristics of gullies and the high-risk areas of gullies. Moreover, we can further judge whether there are significant differences in the amount of erosion between different levels of influencing factors (Table8). Among them, rainfall erosivity, elevation, slope shape, slope aspect, slope position, slope length, slope, residential distance, distance from river, distance from road, catchment area, terrain wetness index, soil type, soil erodibility, lithology, and land use display significant differences among different levels. Among the rainfall factors, the average value of gully density and intensity is the highest when the rainfall erosivity is between 210 and 211. Among the topographic factors, when the elevation is 220-260 m and the slope is 4-8°, the average value of gully density and intensity is the highest, indicating that the risk of gully erosion is higher on flat land. The lithology belongs to gray–white matter rhyolite, spherulite rhyolite, and crystal clastic tuff, with the highest gully density and intensity.Table 8
High-risk areas of gully erosion and its mean value.
FactorsHigh-risk area of gully densityAverage gully density (m/km2)High-risk area of gully intensityAverage gully intensity (m2/km2)Rainfall erosivity210-2111281.461210-21111752.866Elevation220-260 m1197.326220-260 m7991.692Slope shapeConvex slope856.4095653Convex slope5781.986Slope directionWest slope948.096West slope6605.746Slope positionMidslope932.913Midslope6706.873Slope degree4-8°1218.9204-8°8326.412Slope length0-150 m834.0170-150 m5588.848Residential distance<0.5879.799<0.55806.133Distance from the river250 m-500 m955.356250 m-500 m6701.476Distance from the road100 m-200 m906.231200 m-300 m6452.205Plane curvature0-5857.6840-55802.145Sectional curvature-25–-5880.6405-155990.171Catchment area2-3 km21223.7242-3 km210554.147Topographic wetness index5-7.5939.9085-7.56449.666Soil typeHumus dark brown soil1923.932Humus dark brown soil17149.954Soil erodibility0.27-0.291250.7450.27-0.299139.993LithologyGray–white matter rhyolite, globular rhyolite, crystal clastic tuff3542.989Gray–white matter rhyolite, globular rhyolite, crystal clastic tuff22833.733Land useCultivated land1027.192Cultivated land6656.717According to risk detection, on the whole, the soil and lithology are the main factors affecting the distribution of gullies. Regarding land use types, the main influence on the distribution of gully erosion is cultivated land. This may be due to the shortage of woodland caused by residential farming, which further leads to soil erosion. Thus, it can be seen that returning farmland to forests may slow down the development trend of gullies and improve the erosion conditions.


## 4. Discussion

### 4.1. Assessment of Gully Erosion Mapping Based on High-Resolution Satellite Imagery

With improvements in image resolution, the ground area represented by a single pixel becomes smaller and smaller, and the same gully is more likely to appear as pure pixels in the image, while the contour characteristics and internal details of the gullies become clearer. In the opposite case, it will appear in the form of mixed pixels, and the contour features of gully erosion will become more blurred. As shown in Figure4, the number of pixels constituting the same area of gully erosion increases in the order of “15 m-Landsat8, 3 m-Planet Labs, 2.5 m-Alos, 2 m-ZY3, 0.7 m-Pleiades, 0.5 m-DMC, and 0.04 m-UAV.” The number of mixed pixels decreases gradually, and the outline of the gully erosion becomes clearer. Especially for the ramified gully system, when the image resolution increases from 2.5 m, the contour of the gully becomes increasingly recognizable, and the land plots between the gullies can also be distinguished. The ridge planting direction of the land plots between the gully can even be seen in the Pleiades image.Figure 4
Image characteristics of branch gully on different images ((a) 15 m-Landsat8, (b) 3 m-Planet Labs, (c) 2.5 m-Alos, (d) 2 m-ZY3, (e) 0.7 m-Pleiades, (f) 0.5 m-DMC, and (g) 0.04 m-UAV).
(a)(b)(c)(d)(e)(f)(g)To alleviate the influence of the imaging time and the spectral resolution of different remote sensing images on the analysis results, the resampled Pleiades image was analyzed. As shown in Figure5, the permanent gully and ephemeral gully on the sloping farmland are visible in the 0.7 m Pleiades image. According to the actual measurement, the average width of the permanent gully marked by the red circle on the drawing is 2.3 m, and the average width of the ephemeral gully is 0.32 m. The contour of the gully is clear, and the ephemeral gully is very easy to identify. However, with the decrease in resolution, the ephemeral gully and permanent gully begin to become blurred. When the resolution is less than 1.5 m, the ephemeral gully becomes difficult to identify. When the resolution is 3 m, the ephemeral gully is “submerged” in the mixed pixels and cannot be identified in the image. Although there is still a gully on the 3 m resolution image, its contour and shape cannot be recognized, and the recognizable length of the trench becomes smaller.Figure 5
Image characteristics of the permanent gully and ephemeral gully on Pleiades images with different resolutions ((a) 0.7 m, (b) 1 m, (c) 1.5 m, (d) 2 m, (e) 2.5 m, (f) 3 M, and (g, h) photos of the gullies at ground).An improvement in the remote sensing image resolution will not only make the image contour characteristics of the gully more obvious but also make its internal structure characteristics clearer. On the scale of 1 : 1000, the internal features of submeter images have little visual difference. However, the opinion that the higher the spatial resolution of remote sensing image, the better remains to be discussed. This is because the improvement in resolution will also bring more noise, which greatly increases the likelihood of obtaining different objects with the same spectral characteristics or the same object with different spectra. The spectral characteristics of the gully mainly depend on the vegetation or soil on its surface, while the spectral values of vegetation or bare soil in the gully are usually consistent with other vegetation or bare soil in the surroundings. The field path is usually easy to misjudge as an ephemeral gully because the spectral and geometric features between them are relatively similar. The improvement in resolution not only improves the ability to identify the ephemeral gully but also improves the ability to find the field path, which is more likely to increase the possibility of misjudging the field path as an ephemeral gully, interfering with the interpretation of the gully. Therefore, the improvement in image resolution will increase the interference noise in remote sensing identification and impede the information extraction of the gully.


### 4.2. Analysis of Controlling Factors of Gully Erosion

In this study, the geodetector-based quantitative attribution results showed that lithology, soil type, rainfall erosivity, catchment area, elevation, soil erodibility, land use, slope degree, and distance from the river have a strong influence on the development of gully erosion. Other factors, such as TWI, slope shape, slope aspect, slope position, slope length, distance from the residential area, distance from the road, plane curvature, and section curvature, cannot sufficiently explain the changes and development of gully erosion. However, the interaction of these factors with lithology, soil type, and land use shows a stronger influence on the development of gully erosion.The formation and expansion of gullies is commonly influenced by particular soil characteristics and behavior. The most relevant properties are likely soil texture characteristics (e.g., percentage of sand, silt, and clay) and soil organic carbon content [46]. Moreover, the underlying lithology can play an important role in determining the occurrence and dimensions of gullies [47]. In this study, soil and lithology played a more important role. The explanatory power of lithology with regard to the spatial distribution of gullies was more than 10%. When the lithology belonged to gray–white matter rhyolite, globular rhyolite, and crystal clastic tuff, the gully erosion risk reached the highest level. The explanatory power of soil type regarding gully density was close to 10%, and the explanatory power for gully intensity was more than 10%, indicating that this is the dominant factor. When the soil type was humus dark brown soil, the risk of gully erosion was the greatest.Topographic variables play a key role in the prediction of both gully initiation and expansion. The present results showed that elevation, slope, and catchment area have a significant influence on the occurrence of gully erosion. These findings are similar to those of previous studies in the black soil region [5, 48]. Climate and weather conditions, and especially rainfall, are key drivers of gully erosion. According to the results of the geodetector model, it can be seen that the impact of rainfall erosivity on the distribution of gully erosion is obvious, showing a trend of first increasing, then decreasing, and then increasing. However, due to the lack of meteorological stations in the basin, rainfall erosivity is obtained by interpolating the surrounding meteorological stations’ data, which would increase the uncertainty of the results. Land cover/use can influence gully initiation through its effect on runoff production [49]. In this study, the type of land use had a significant effect on the distribution of gully erosion, and its explanatory power with regard to gully density and intensity was 4.47% and 1.65%, respectively. Among them, cultivated land has the greatest impact, which may be due to the frequent human activities on cultivated land. This result is consistent with some previous studies [50, 51]. For example, the results of Ali Azareh et al.’s study in 2019 showed that slope aspect, lithology, and land use were identified as the most important factors affecting gully sensitivity in Iran by using a maximum entropy model [3]. However, this contradicts our findings. We believe that the reason for the difference may be that the geographical environment and hydrological conditions of the study area were not fully considered in the selection of impact factors. The slope degree of sloping farmland in the black soil region in Northeast China is mostly 3°-10°. As a result, slope has no significant influence on gully erosion in Northeast China. Lithology and land use are the main factors affecting the distribution of gully erosion, which is consistent with our results.According to the results of the geodetector, the interaction of most factors showed nonlinear enhancement, which indicated that a gully is more likely to occur under the combined action of various factors. The most influential interaction groups were lithology with other factors, especially lithology in different soil types. In particular, the lithology of different soil types accounted for 28.7% and 32.5% of gully density and gully intensity, which also reflects the dominant role of lithology and soil. The interaction of lithology and land use also had a significant effect on gully erosion, which demonstrates the importance of returning cropland to forest in the process of gully erosion control in Northeast China. The results show that the interaction between lithology and soil, topography, rainfall, and land use factors can significantly improve the explanatory power in terms of gully erosion development, and the interactiveq values are more than 10%. Only two pairs of factors are independent: (1) slope and lithology and (2) slope length and lithology.The geodetector revealed the single and paired driving factors affecting the spatial differentiation of gully erosion in Northeast China and enriched the research content on gully erosion in this area. The results show that the interaction effect of soil type, land use type, and lithology is stronger than that of a single factor. This shows that the prevention and control of gullies should not only start from the main risk factors but also from the perspective of the whole. However, there are still some shortcomings. Interaction analysis based on a geodetector can only evaluate the interaction of two factors—it cannot analyze the interaction of more than two factors.


## 4.1. Assessment of Gully Erosion Mapping Based on High-Resolution Satellite Imagery

With improvements in image resolution, the ground area represented by a single pixel becomes smaller and smaller, and the same gully is more likely to appear as pure pixels in the image, while the contour characteristics and internal details of the gullies become clearer. In the opposite case, it will appear in the form of mixed pixels, and the contour features of gully erosion will become more blurred. As shown in Figure4, the number of pixels constituting the same area of gully erosion increases in the order of “15 m-Landsat8, 3 m-Planet Labs, 2.5 m-Alos, 2 m-ZY3, 0.7 m-Pleiades, 0.5 m-DMC, and 0.04 m-UAV.” The number of mixed pixels decreases gradually, and the outline of the gully erosion becomes clearer. Especially for the ramified gully system, when the image resolution increases from 2.5 m, the contour of the gully becomes increasingly recognizable, and the land plots between the gullies can also be distinguished. The ridge planting direction of the land plots between the gully can even be seen in the Pleiades image.Figure 4
Image characteristics of branch gully on different images ((a) 15 m-Landsat8, (b) 3 m-Planet Labs, (c) 2.5 m-Alos, (d) 2 m-ZY3, (e) 0.7 m-Pleiades, (f) 0.5 m-DMC, and (g) 0.04 m-UAV).
(a)(b)(c)(d)(e)(f)(g)To alleviate the influence of the imaging time and the spectral resolution of different remote sensing images on the analysis results, the resampled Pleiades image was analyzed. As shown in Figure5, the permanent gully and ephemeral gully on the sloping farmland are visible in the 0.7 m Pleiades image. According to the actual measurement, the average width of the permanent gully marked by the red circle on the drawing is 2.3 m, and the average width of the ephemeral gully is 0.32 m. The contour of the gully is clear, and the ephemeral gully is very easy to identify. However, with the decrease in resolution, the ephemeral gully and permanent gully begin to become blurred. When the resolution is less than 1.5 m, the ephemeral gully becomes difficult to identify. When the resolution is 3 m, the ephemeral gully is “submerged” in the mixed pixels and cannot be identified in the image. Although there is still a gully on the 3 m resolution image, its contour and shape cannot be recognized, and the recognizable length of the trench becomes smaller.Figure 5
Image characteristics of the permanent gully and ephemeral gully on Pleiades images with different resolutions ((a) 0.7 m, (b) 1 m, (c) 1.5 m, (d) 2 m, (e) 2.5 m, (f) 3 M, and (g, h) photos of the gullies at ground).An improvement in the remote sensing image resolution will not only make the image contour characteristics of the gully more obvious but also make its internal structure characteristics clearer. On the scale of 1 : 1000, the internal features of submeter images have little visual difference. However, the opinion that the higher the spatial resolution of remote sensing image, the better remains to be discussed. This is because the improvement in resolution will also bring more noise, which greatly increases the likelihood of obtaining different objects with the same spectral characteristics or the same object with different spectra. The spectral characteristics of the gully mainly depend on the vegetation or soil on its surface, while the spectral values of vegetation or bare soil in the gully are usually consistent with other vegetation or bare soil in the surroundings. The field path is usually easy to misjudge as an ephemeral gully because the spectral and geometric features between them are relatively similar. The improvement in resolution not only improves the ability to identify the ephemeral gully but also improves the ability to find the field path, which is more likely to increase the possibility of misjudging the field path as an ephemeral gully, interfering with the interpretation of the gully. Therefore, the improvement in image resolution will increase the interference noise in remote sensing identification and impede the information extraction of the gully.


## 4.2. Analysis of Controlling Factors of Gully Erosion

In this study, the geodetector-based quantitative attribution results showed that lithology, soil type, rainfall erosivity, catchment area, elevation, soil erodibility, land use, slope degree, and distance from the river have a strong influence on the development of gully erosion. Other factors, such as TWI, slope shape, slope aspect, slope position, slope length, distance from the residential area, distance from the road, plane curvature, and section curvature, cannot sufficiently explain the changes and development of gully erosion. However, the interaction of these factors with lithology, soil type, and land use shows a stronger influence on the development of gully erosion.The formation and expansion of gullies is commonly influenced by particular soil characteristics and behavior. The most relevant properties are likely soil texture characteristics (e.g., percentage of sand, silt, and clay) and soil organic carbon content [46]. Moreover, the underlying lithology can play an important role in determining the occurrence and dimensions of gullies [47]. In this study, soil and lithology played a more important role. The explanatory power of lithology with regard to the spatial distribution of gullies was more than 10%. When the lithology belonged to gray–white matter rhyolite, globular rhyolite, and crystal clastic tuff, the gully erosion risk reached the highest level. The explanatory power of soil type regarding gully density was close to 10%, and the explanatory power for gully intensity was more than 10%, indicating that this is the dominant factor. When the soil type was humus dark brown soil, the risk of gully erosion was the greatest.Topographic variables play a key role in the prediction of both gully initiation and expansion. The present results showed that elevation, slope, and catchment area have a significant influence on the occurrence of gully erosion. These findings are similar to those of previous studies in the black soil region [5, 48]. Climate and weather conditions, and especially rainfall, are key drivers of gully erosion. According to the results of the geodetector model, it can be seen that the impact of rainfall erosivity on the distribution of gully erosion is obvious, showing a trend of first increasing, then decreasing, and then increasing. However, due to the lack of meteorological stations in the basin, rainfall erosivity is obtained by interpolating the surrounding meteorological stations’ data, which would increase the uncertainty of the results. Land cover/use can influence gully initiation through its effect on runoff production [49]. In this study, the type of land use had a significant effect on the distribution of gully erosion, and its explanatory power with regard to gully density and intensity was 4.47% and 1.65%, respectively. Among them, cultivated land has the greatest impact, which may be due to the frequent human activities on cultivated land. This result is consistent with some previous studies [50, 51]. For example, the results of Ali Azareh et al.’s study in 2019 showed that slope aspect, lithology, and land use were identified as the most important factors affecting gully sensitivity in Iran by using a maximum entropy model [3]. However, this contradicts our findings. We believe that the reason for the difference may be that the geographical environment and hydrological conditions of the study area were not fully considered in the selection of impact factors. The slope degree of sloping farmland in the black soil region in Northeast China is mostly 3°-10°. As a result, slope has no significant influence on gully erosion in Northeast China. Lithology and land use are the main factors affecting the distribution of gully erosion, which is consistent with our results.According to the results of the geodetector, the interaction of most factors showed nonlinear enhancement, which indicated that a gully is more likely to occur under the combined action of various factors. The most influential interaction groups were lithology with other factors, especially lithology in different soil types. In particular, the lithology of different soil types accounted for 28.7% and 32.5% of gully density and gully intensity, which also reflects the dominant role of lithology and soil. The interaction of lithology and land use also had a significant effect on gully erosion, which demonstrates the importance of returning cropland to forest in the process of gully erosion control in Northeast China. The results show that the interaction between lithology and soil, topography, rainfall, and land use factors can significantly improve the explanatory power in terms of gully erosion development, and the interactiveq values are more than 10%. Only two pairs of factors are independent: (1) slope and lithology and (2) slope length and lithology.The geodetector revealed the single and paired driving factors affecting the spatial differentiation of gully erosion in Northeast China and enriched the research content on gully erosion in this area. The results show that the interaction effect of soil type, land use type, and lithology is stronger than that of a single factor. This shows that the prevention and control of gullies should not only start from the main risk factors but also from the perspective of the whole. However, there are still some shortcomings. Interaction analysis based on a geodetector can only evaluate the interaction of two factors—it cannot analyze the interaction of more than two factors.


## 5. Conclusion

The present work represents a contribution to the multifactor synthesis for gully erosion at a large watershed scale. The Sancha River basin, a typical area in the black soil region of Northeast China, was taken as the study area. A high-resolution satellite image was used to obtain the spatial distribution of the gullies by visual image interpretation with field verification. We analyzed the dominant factors affecting the spatial distribution of gully erosion and the degree of interaction between any two of these factors using the geodetector method, and we identified areas at high risk of gully erosion between strata of each factor. The following conclusions were obtained:(1)
At the large catchment scale, the submeter images show a strong capacity for the recognition of permanent gullies and obtained satisfactory results(2)
According to the results of the geodetector, lithology and soil type are the main factors that affect the spatial differentiation of gully erosion in the Sancha River basin. The interpretation power of gully density and gully intensity is close to 10%, and the explanatory power of gully erosion occurrence is greater. As rainfall erosivity, watershed area, elevation, soil erodibility, land use pattern, slope, and distance from the river accounted for more than 1% of gully density and gully intensity, they were identified as secondary factors(3)
The interaction detection results of the geodetector show that there is a nonlinear enhancement relationship between lithology and most factors. In particular, the lithology of different soil types accounted for 28.7% and 32.5% of gully density and gully intensity. The results show that the interaction between lithology and soil, topography, rainfall, and land use factors can significantly improve the explanatory power of gully erosion development, and the interactiveq values reached more than 10%. The results show that the interaction of several factors has a stronger influence on the spatial differentiation of gully erosion than a single factor. The prevention and control of gullies should not only start from the main risk factors but also from the perspective of the whole(4)
The lithology belongs to gray–white matter rhyolite, spherulite rhyolite, and crystal clastic tuff, with the highest gully density and intensity



---

*Source: 3831652-2022-07-20.xml*