2.3.1. RUSLE Model

The RUSLE model is a modified version of the Universal Loss Equation (USLE). Due to its simple structure and ease of operation, scholars worldwide widely use it to estimate soil erosion. Its equation is as follows:

$$A = R \times K \times L \times S \times \mathbb{C} \times P \tag{1}$$

where *<sup>A</sup>* is the average annual soil erosion (t·ha−1·a−1), which means the average annual soil loss from fine gully or inter-gully erosion on slopes caused by rainfall and runoff. *<sup>R</sup>* is the rainfall erosivity factor (MJ·mm·ha−1·h−1·a−1), *<sup>K</sup>* is the soil erodibility factor (t·ha·h·MJ−1·mm−1·ha<sup>−</sup>1), *<sup>L</sup>* is the slope length factor (dimensionless), *<sup>S</sup>* is the slope gradient factor (dimensionless), *C* is the vegetation cover and management factor (dimensionless), and *P* is the soil and water conservation measures factor (dimensionless).

Gao modified the RUSLE model based on the correlation coefficients between exposed bedrock and surface sediments [23], with the following expressions:

$$A = \left(1 - 0.076^2 \times \alpha\right) R \times K \times L \times S \times \mathbb{C} \times P \tag{2}$$

where α is a correction factor whose value is determined by the average of the bedrock exposure rates for the different rocky desertification classes (Table 1). Different equations were chosen to simulate soil erosion depending on the lithology of the study area. Equation (1) was selected for non-karst regions, and Equation (2) was selected to simulate soil erosion in karst regions.

**Table 1.** α values for different karst rocky desertification classes.


The K factor is a measure of soil erosion resistance and reflects the sensitivity of the soil to erosion. The *K* factor is calculated using the erosion-productivity impact calculator (EPIC) model (Equation (3)) developed by Williams [36].

$$K = \begin{cases} 0.2 + 0.3 \exp\left[0.0256 \text{SAN} \left(1 - \frac{SIL}{100}\right)\right] \Big| \left(\frac{SIL}{\text{CLA} + \text{SIL}}\right)^{0.5} \\ 1.0 - \frac{0.25 \text{C}}{\text{C} + \exp(3.72 - 2.95)} \Big[ 1.0 - \frac{0.7 \text{SN}\_1}{\text{SN}\_1 + \exp\left(-5.51 + 22.9 \text{SN}\_1\right)} \Big] \end{cases} \tag{3}$$

where *<sup>K</sup>* is the soil erodibility factor in t·hm2·h·MJ−1·mm−1·hm−2. *SAN* is the sand content; *SIL* is the silt content; *CLA* is the clay content; *C* is the organic carbon content, *SN*<sup>1</sup> = 1 − *SAN*/100.

The *R* factor is the driving force behind soil erosion and reflects the potential for soil loss through precipitation. As there are no weather stations in the study area, the nearest weather station to the study area was selected and the *R* factor for the study area was obtained by interpolation using the rainfall erosion equation (Equation (4)) created by Arnoldus [37]. The equation is as follows:

$$R = \sum\_{i=1}^{12} 1.735 \ast 10^{\left[1.5 \ast \log\left(\frac{P\_i}{P}\right) - 0.8188\right]}\tag{4}$$

where *R* is the rainfall erosion factor, *Pi* is the monthly rainfall (mm), and *P* is the annual rainfall (mm).

The *L* factor and *S* factor are closely related to topography and accumulated flow, and *L* factor and *S* factor together reflect the influence of topographic features on soil erosion. LS represents the ratio of soil loss on a given slope length and gradient to soil loss on a typical slope in a standard runoff plot, all other things being equal, and the LS value is proportional to soil loss, playing an accelerating role in soil erosion [38]. We chose to calculate the *L* factor [39] and *S* factor [40] by applying the formulae to the soil and rocky mountainous areas of southwest China as follows:

$$L = \left(\frac{\lambda}{22.13}\right)^m \tag{5}$$

$$m = \frac{\beta}{\beta + 1} \tag{6}$$

$$\beta = (\sin \theta / 0.0896) / [3 \* (\sin \theta)^{0.8} + 0.56] \tag{7}$$

$$S = \begin{cases} 10.80 \ast \sin \theta + 0.03 \ (\theta < 5^{\circ}) \\ 16.80 \ast \sin \theta - 0.50 \ (5^{\circ} \le \theta \le 10^{\circ}) \\ 20.204 \ast \sin \theta - 1.2404 \ (10^{\circ} \le \theta > 25^{\circ}) \\ 29.585 \ast \sin \theta - 5.6079 \ (\theta \ge 25^{\circ}) \end{cases} \tag{8}$$

where *λ* is the sum of the slope lengths in the horizontal direction, *m* is the slope length factor, *β* is a factor related to the slope value, and *θ* is the slope angle extracted on the basis of the digital elevation.

*C* factor refers to the ratio of soil loss under a particular crop or vegetation cover to that of continuous recreational land after cultivation, all other factors being equal, and this factor measures the inhibitory effect of vegetation cover and management on soil erosion [41,42]. Currently, there are two forms for obtaining *C* factor values, equation calculation and assigning values on the basis of land use type. In this paper, the formula created by Durigon [43] was chosen to calculate the *C* value as follows:

$$\mathcal{C} = \frac{-NDVI + 1}{2} \tag{9}$$

where *NDV I* is the normalized vegetation index, calculated from the near-infrared band and visible red band in remote sensing imagery.

*P* factor refers to the ratio of soil erosion with soil conservation measures to soil erosion with down-slope planting, and its value is distributed between 0 and 1. When *p* = 1, it means that no water conservation measures are taken in the unit area; when *p* = 0, it means that no soil erosion will theoretically occur in the unit area, such as water bodies, construction sites, etc. The *p*-values were obtained by summarizing the results of previous authors in the karst region of southern China [23,44,45]. The results are presented in the Table 2.

**Table 2.** Land use types and *p*-value.


#### 2.3.2. Geographical Detector

The geographical detector is a statistical method that allows for the exploration of spatial dissimilarity and its drivers [46]. The principle of the geographical detector is based on the idea that if the independent variable has a significant effect on the dependent variable, then there should be some similarity in the spatial distribution of the two [47]. It consists of four models, namely, factor detector, interaction detector, risk detector, and ecological detector. In this paper, we analyzed the spatial heterogeneity with the factor detector, interaction detector, and risk detector.

The factor detector detects the extent to which the independent variable can explain the spatial divergence of the dependent variable, and the *q* value can measure the extent of explanation with the following expression:

$$q = 1 - \frac{SSW}{SST} \tag{10}$$

$$SSW = \sum\_{h=1}^{L} N\_h \sigma\_{h'}^2 \text{ SST} = N\sigma^2 \tag{11}$$

where *h* = 1,2,..., *L* is the partition of the independent and dependent variables; *Nh* and *N* are the number of cells in stratum *h* and all partitions; *σ*<sup>2</sup> *<sup>h</sup>* and *<sup>σ</sup>*<sup>2</sup> are the variances of the dependent variables in stratum *h* and all partitions. *SSW* and *SST* are the sum of the within-stratum variances and the total variance of all partitions, respectively. *q* has a value of [0,1], with larger values indicating that the independent variable explains more of the dependent variable.

The interaction detector can be used to detect the interaction between factors *XS*, i.e., the change in the degree of explanation of the dependent variable Y when *X*<sup>1</sup> and *X*<sup>2</sup> act together. The principle is to calculate the *q* values of *X*<sup>1</sup> and *X*<sup>2</sup> separately, then to calculate the new layer *q*(*X*<sup>1</sup> ∩ *X*2) values by superimposing the *X*<sup>1</sup> and *X*<sup>2</sup> layers and comparing them among the three. The comparison between the two factors can be divided into the following results: non-linearly diminished, one-factor non-linearly diminished, two-factor enhanced, independent, and non-linearly enhanced.

The risk detector can be used to count the mean values of attributes between subregions of a single factor and to determine whether the heterogeneity between sub-regions is significant. The mean attribute results are expressed as numerical data, while whether the means are significantly different is expressed as 'N/Y' binary data.

#### 2.3.3. Slope Units

The slope cell delineation method based on the principle of curvature-based watershed segmentation [48] proceeds as follows (Figure 2).

**Figure 2.** Flow chart of slope unit division.
