Soil Conservation Service Spatiotemporal Variability and Its Driving Mechanism on the Guizhou Plateau, China

Abstract

The Guizhou Plateau has an extremely fragile ecological environment with prominent soil and water losses. Since 2000, conservation policies and ecological restoration projects, e.g., the Grain for Green Project (GGP), have been implemented on the Guizhou Plateau to control soil/water losses which have achieved notable accomplishments. Using the Revised Universal Soil Loss Equation (RUSLE) to estimate the soil conservation service (SCS) on the Guizhou Plateau, this study analyzed the dynamic characteristics of its spatiotemporal variation based on multiyear (2000–2018) meteorological and remote sensing data to determine its driving mechanisms. Residual analysis of the meteorological and remote sensing data was used to evaluate the effect of anthropogenic activities. Results showed a clear upward trend (1.39 t ha⁻¹ yr⁻¹) of SCS on the Guizhou Plateau during 2000–2018, and areas with a highly improved positive effect on SCS were distributed primarily in karst landform regions. Precipitation and vegetation fractional coverage (VFC) were found to be positively correlated with SCS on the Guizhou Plateau. Specifically, the highest proportion of significant positive correlation between precipitation and SCS was related to the Wildlife Conservation Nature Reserve (WCNR), and the highest proportion of significant positive correlation between VFC and SCS was related to the GGP, i.e., 76.59% and 53.02%, respectively. Residual analysis revealed a significant positive role of anthropogenic activity on SCS improvement via ecological engineering in areas with a poor ecological background, e.g., the GGP in western areas where the ecological environment is fragile and the problem of water/soil loss is serious. In areas with a more robust ecological background, e.g., the engineering area of the WCNR, the effect of anthropogenic activity has had a largely negative effect on SCS. The findings of this study could make an important contribution to the development of ecological management projects and the work to control soil/water losses on the Guizhou Plateau.

1. Introduction

The soil conservation service (SCS), which is an important regulating service of ecosystem services [1], plays a critical role in both the conservation of water and soil and the improvement of soil structure through ground cover and litter [2]. SCS provides not only an important foundation for human agricultural production, but also protection against the risk of flooding and deterioration of the ecological environment [3,4]. Aggravation of soil erosion can lead to serious soil and water losses, that result in ecosystem fragmentation and decline in both soil fertility and land quality [5], which could be detrimental to the ecological environment and the development of human society. Therefore, SCS evaluation is an important aspect in ecosystem service analysis.

The difference between the potential soil erosion and actual soil erosion can be used as a quantitative indicator of SCS [6]. In recent years, this method has been used widely to evaluate SCS on watersheds [7], regional [8], and other scales. This model is also particularly useful for simulating the dynamic changes of SCS under the effects of climatic and land use changes [9,10]. Many models have been proposed for the calculation of soil erosion. The United States Department of Agriculture proposed the Water Erosion Prediction Project model, which is a physical process model that can dynamically monitor the slope erosion process [11]; Morgan et al. [12] proposed the European soil erosion model, which divided the soil erosion instead of erosion between fine trenches and fine trenches is applied to the prediction of secondary rainfall erosion; De Roo et al. [13] proposed the Dutch soil erosion prediction model, which considers various processes in simulation of the process of soil erosion. The Revised Universal Soil Loss Equation (RUSLE) [11], issued by the United States Department of Agriculture, can be used for the assessment of soil erosion loss. In comparison with many other models, the RUSLE model has a simpler structure, a clearer physical meaning of parameters, and stronger practicable and comprehensive capabilities [14]. Ganasri and Ramesh [15] used the RUSLE model to estimate the amount of soil erosion in the Nethravathi Basin of southern India. Wu et al. [16] assessed the magnitude of soil erosion and its spatial distribution characteristics in the Mekong Basin. Wang et al. [17] studied soil erosion in the water source area of Shaanxi (China) using the RULSE model. The increased use of geographic information systems (GIS) and remote sensing (RS) technology has allowed the evaluation of ecological services to be supported by long-term dynamic information [18]. Li and Zheng [19] studied soil erosion changes during 2001–2010 in the Yanhe watershed on the Loess Plateau (China) with GIS and RS support. Wang et al. [20] investigated the degree of soil erosion in the area of the Three Gorges Reservoir (China) during 2000–2010 using the RUSLE model in conjunction with a Moderate-resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) map and digital elevation model (DEM) data. Olorunfemi et al. [21] assessed potential soil erosion in Ekiti State (Nigeria) with a GIS platform and monthly rainfall data interpolated from climate point data using the inverse distance weighting method. The common data sources of RUSLE-based models are meteorological and RS data [22].

Topography is an important factor that affects the area and distribution of soil erosion, and slope has important influence on the mechanical composition of soil and runoff [23,24]. Slope-based runoff experimental methods [25] for measuring soil erosion tend to be elaborate, so cannot be used for large-scale simulation studies. When using prediction methods such as the RUSLE model, the accuracy of the topographical factor is extremely important. Although a DEM can represent an effective source of GIS-based topographic information, the topographical factor extracted from DEMs at different scales can affect the accuracy of the final soil erosion calculation results [26]. Consequently, the scale of DEM data should be selected depending on the scale and objectives of any particular study [27]. For large-scale study areas, the use of low-resolution DEM data would erase certain details of the surface terrain, thereby affecting the accuracy of the topographical factor. Conversely, the use of high-resolution data would exacerbate problems associated with the storage and processing of large volumes of data. In addition, for small-scale soil erosion studies and hydrologic-level analyses, it can be necessary to build very high resolution DEMs, even to runoff plot or hydrological observation scale [28]. Currently, satellite-derived high-resolution RS imagery (e.g., from the Sentinel 2A satellite) and DEMs acquired using unmanned aerial vehicle technology can be incorporated in quantitative studies of soil erosion [29].

Climate change and anthropogenic activities have altered the structure and function of ecosystems. Related studies have shown that global soil erosion is at increased risk under the influence of climate change and unrestrained anthropogenic activities [30]. Anthropogenic activity is the main driving factor affecting soil conservation that can change soil erosion [31]. Residual analysis is often used to study the driving mechanisms or to evaluate the impact of anthropogenic activities. For example, Wu et al. [32] undertook a quantitative assessment of the impact of weather conditions on air quality in Hohhot (China) using residual analysis. Luo et al. [33] also used residual analysis with MODIS-NDVI data to analyze the effect of climate change and anthropogenic activities on vegetation on the island of Hainan (China).

China is among the countries within which the problem of soil erosion is most serious. By 2011, the total area of land within China affected by soil erosion had reached 294.92 × 104 km² [34]. However, strengthening ecological awareness has led China to implement ecological engineering projects for comprehensive control of soil and water losses, e.g., the Shelterbelt Forest Project (SFP), Grain for Green Project (GGP), and Comprehensive Control Project of Rocky Desertification (CCPRD) [35]. Previous research has shown that ecological projects such as GGP could control soil erosion [36]. However, previous studies have not examined the spatiotemporal response of SCS to ecological engineering activities in regions with different landforms, especially in the karst areas. Moreover, few quantitative studies have investigated the roles of climate change (e.g., in relation to rainfall) and anthropogenic activities in driving SCS under the background of ecological engineering. These subjects are of great importance with regard to further improvement of the layout and implementation of ecological engineering constructions and projects for soil erosion control in specific areas.

The Guizhou Plateau, which is one of the regions of China with the most serious soil and water losses, has an extremely fragile ecological environment. It is a region characterized by its typical karst landforms, rocky desertification, and serious soil erosion; the area of soil/water losses accounts for 43.5% of the total area [37]. The Guizhou Plateau, which is located at the confluence of the Yangtze and Pearl rivers, is an important ecological defense in the upper reaches of both river basins. Soil erosion within this region is related directly to the ecological security of the river basins and has considerable impact on the socioeconomic development of downstream areas. To improve the increasingly prominent problems of soil erosion and rocky desertification on the Guizhou Plateau, ecological restoration projects such as the CCPRD have been undertaken on a large scale since 2004 [38]. Zhang et al. [39] quantified the contribution rate of anthropogenic activities to vegetation coverage on the Guizhou Plateau under the effects of ecological engineering construction. It was found that the positive effects of anthropogenic activities outweighed the negative effects. Based on studies in the Shaanxi region of China, Liu et al. [40] concluded that GGP implementation has reduced soil erosion and enhanced the soil protection effect. The implementation of large-scale ecological projects raises certain important scientific concerns. For example, it is imperative to ascertain the changes that have taken place in terms of SCS on the Guizhou Plateau, determine the characteristics of SCS on the Guizhou Plateau under different ecological engineering constructions, examine the impact of ecological engineering on SCS, establish the characteristics of SCS in regions with different landforms regions, and investigate whether ecological engineering construction has been adapted sufficiently to local conditions.

The Guizhou Plateau is located in the largest karst area in the world, where population pressure is greatest, land degradation serious, and rocky desertification prominent. Consequently, the ecological restoration of land degradation characterized by soil erosion is a difficult problem [41]. The fragile ecological environment of the Guizhou Plateau makes it highly susceptible to the effects of anthropogenic activities. China has implemented many ecological projects in this region. However, few studies have considered the spatiotemporal variation and driving factors of SCS in relation to ecological engineering construction on the Guizhou Plateau.

This study analyzed the spatiotemporal variation and driving factors of SCS on the Guizhou Plateau during 2000–2018 using GIS in combination with long-term meteorological and RS data. First, the RUSLE model was used to calculate SCS. Then, the comparative influence of natural factors (e.g., precipitation and vegetation fractional coverage (VFC)) and anthropogenic activities on SCS was determined through regression analysis and residual analysis, following which the evolution trend of SCS and the effect of different ecological engineering projects on SCS were assessed. The objectives of this study were to evaluate the variability of SCS on the Guizhou Plateau during 2000–2018, elucidate its driving mechanisms, clarify the effects of anthropogenic activities under ecological restoration projects in regions with different landforms, and propose suggestions for reducing soil erosion and rational planning of ecological engineering areas. The findings will make an important contribution to both the work of controlling soil and water losses on the Guizhou Plateau area and the development of an ecologically aware civilization, as well as provide scientific reference regarding ecological security problems in the basin of the Yangtze and Pearl rivers.

2. Materials and Methods

2.1. Study Area

The Guizhou Plateau is in southern China (24°37′–29°13′N, 103°36′–109°35′E). It has high elevation in the west and low elevation in the east (average elevation: 1100 m), and the landforms are predominantly typical karst-type landforms. According to the “Techniques standard for comprehensive control of soil erosion and water loss in karst region” (SL 461-2009) [42] and related studies [43,44], based on landform type and soil erosion, the study area can be divided into six landform regions: karst canyon, peak-cluster depression, fault depression basin, karst plateau, karst trough valley, and non-karst landform (Figure 1). The Guizhou Plateau has various soil types, and the zonal soil that includes red soil, yellow soil, and yellow-brown soil, accounts for 60.29% of the total area. Vegetation coverage on the Guizhou Plateau is high and the range of vegetation types diverse. The zonal vegetation is subtropical evergreen broadleaved forest. The area has a subtropical humid monsoon climate with little temperature change. The annual average temperature is in the range of 5–22 °C, and the annual average precipitation is in the range of 687–1657 mm. Owing to the characteristics of the landforms, soil types, climate, and hydrology, the Guizhou Plateau is highly susceptible to soil and water losses.

With the objective of improving the ecological environment, the Chinese government has implemented a series of ecological engineering projects in the Guizhou Plateau region. By 2015, the main ecological engineering projects within the region comprised the Natural Forest Protection Project (NFPP), Wildlife Conservation Nature Reserve (WCNR), Key Public Welfare Forest Project (KPWFP), Shelterbelt Forest Project (SFP), CCPRD, GGP, and some other forestry projects (OFP). The combined areas (39,700 km²) of these projects, of which the largest two are the NFPP followed by the GGP, account for approximately 22.56% of the total area of the plateau.

2.2. Data and Processing

This study used meteorological data and RS data to evaluate SCS. The main methodological workflow is shown in Figure 2.

The meteorological data used in this study included precipitation, temperature, and other meteorological variables (1990–2018) recorded at 40 national meteorological stations located in the study area and within the surrounding 150 km (Figure 3). The dataset was obtained from the National Meteorological Information Center of China [45].

The ANUSPLIN professional meteorological interpolation software package, which is designed to target climate data, was used to interpolate the temperature and precipitation data. Although GIS software includes multiple spatial interpolation methods, ANUSPLIN specifically considers the accuracy, convenience, and time series length of the data and it has the advantage of being based on the theory of the thin plate spline function that produces high interpolation accuracy [46]. We used longitude and latitude which projected into the Krasovsky−1940 coordinate system as independent variables with which to interpolate precipitation to a spatial interpolation resolution of 250 m. The interpolated air temperature and precipitation raster data passed the generalized cross-validation procedure and the verification accuracy met the data analysis requirements. For consistency with the spatial resolution of the RS data, we processed both the precipitation and the temperature data to 250-m resolution to match the other datasets used in this study (Figure 4).

Soil data, based on the basic mapping unit for soil classes, were obtained from the China Soil Database [47], which included 12 soil classes, 61 soil types, and 227 soil subtypes. The attributes of soil sand content, soil clay content, soil silt content, and soil organic carbon content were used in the RUSLE model calculation.

The Digital Elevation Database of the Shuttle Radar Topography Mission (SRTM3) v4.1 was used to derive topographic data with 90-m resolution. These data were provided by the Geospatial Data Cloud site of the Computer Network Information Center (Chinese Academy of Sciences) [48]. The SRTM terrain data were obtained by the International Center for Tropical Agriculture using a new interpolation algorithm. SRTM terrain data can be divided into SRTM1 and SRTM3 datasets that have a corresponding resolution accuracy of 30 and 90 m, respectively. Currently, the published data have 90-m resolution (SRTM3), and these data were used in this study to calculate both the slope length and the slope of the terrain using an algorithm based on Zhang et al. [49].

The vegetation data used in this study comprised NDVI data for 2000–2018 derived from the MODIS NDVI dataset, which has 16-d temporal resolution and 250-m spatial resolution. We cropped the NDVI data according to match the study area and removed invalid fill values. The Maximum Value Composites method was used to extract the maximum value of the two mid-month NDVI values as the monthly NDVI value. This method, which can eliminate some of the effects of the atmosphere, satellite orbit drift, and other disturbances, has been proven to produce reasonable results for mountainous areas with complex terrain [39]. Then, the average of the monthly NDVI values is taken as the annual NDVI value. This method, which can characterize the annual average vegetation cover of an area, can be used to produce basic data for research and analysis. The spatial and temporal resolutions of these data are suitable for regional-scale analysis of the Guizhou Plateau. The VFC data were calculated using ENVI software based on the pixel dichotomy model. The pixel dichotomy method has been used in many studies to construct a model for vegetation coverage extraction based on satellite RS technology in combination with multispectral vegetation indices [50].

Land use data were derived from the China Land Use/Cover Dataset of 2000–2015 that includes 6 land grades and 25 land classes such as cultivated land, forest, grassland, construction land, and unused land. This dataset, which was created by Liu et al. [51], was obtained through considerable human–computer interactive interpretation based on China’s land use classification system and Landsat RS satellite data; the comprehensive evaluation accuracy is >94%.

Details of the area of ecological engineering and other related forest data related to the study area were obtained from the Fourth Survey of Forest Resources Planning and Design in Guizhou Province, provided by the Department of Forestry of Guizhou Province. These data, which cover the entire plateau, comprise 3,225,627 class records that include more than 100 property fields such as land use, altitude, slope direction, slope position, engineering category, type of land degradation, and desertification type. We mainly used the vector format of the boundary of ecological engineering areas for zonal statistics. The types and sources of the data used in this study are shown in

Table 1. Data type, resolution, and source.

Data Type	Resolution or Spatial Distribution	Data Source
Meteorological data (precipitation and temperature)	40 points	National Meteorological Information Center
Soil data	1:1,000,000	Resource and Environment Data Cloud Platform
SRTM-DEM	90 m	Computer Network Information Center, CAS
MODIS-NDVI	250 m	Computer Network Information Center, CAS
Land use dataset	1 km	Institute of Geographic Sciences and Natural Resources Research, CAS
The boundary of Eco-engineering	shapefile	Department of Forestry of Guizhou Province

Note: CAS—Chinese Academy of Sciences.

. The main software to process these data was ArcGIS.

2.3. Methods

2.3.1. Soil Conservation Service Simulation

In this study, the SCS capacity was evaluated based on the soil erosion modulus. We used the SCS per unit area of soil conservation amount [52] to evaluate the SCS capacity in the study area. In this study, SCS was estimated using the RUSLE model [53]. The SCS represents the difference between the potential soil erosion under no vegetation cover and the actual soil erosion under the pertaining land use or land cover (Equations (1)–(3)).

$$\text{Potential Soil Erosion}:A_P=R\times K\times LS,$$

(1)

$$\text{Actual Soil Erosion}:A_a=R\times K\times LS\times C\times P,$$

(2)

$$\mathrm{SCS}:A_c=A_p-A_a,$$

(3)

where A_p is the potential soil erosion modulus [t·ha⁻¹·yr⁻¹], A_a is the actual soil erosion modulus [t·ha⁻¹·yr⁻¹], A_c is the SCS [t·ha⁻¹·yr⁻¹], R is the rainfall erosivity factor [MJ·mm·ha⁻¹·h⁻¹·yr⁻¹], K is the soil erodibility factor [t·h·MJ⁻¹·mm⁻¹], LS is the slope length and slope steepness factor [dimensionless], C is the vegetation cover factor (dimensionless), and P is the soil and water conservation measure factor (dimensionless). All factors were calculated using ArcGIS and processed into 250-m-resolution raster data.

(1) Rainfall erosivity factor (R)

The rainfall erosivity factor is used to indicate the magnitude of rainfall erosivity. It is difficult to obtain data of rainfall energy and rainfall intensity; however, Wang et al. [54] found that the ratio of the sum of the squares of monthly precipitation and yearly precipitation had good correlation with R in karst areas. Hence, we used the model of Arnoldus et al. [55] (Equation (4)) to simulate rainfall erosivity in our study area:

$$R={\displaystyle \sum _{i=1}^{12}}1.735\times {10}^{\left[1.5\log \left(\frac{P_i{}^2}{P}\right)-0.8188\right]}$$

(4)

where R is the rainfall erosivity factor [MJ·mm·ha⁻¹·h⁻¹·yr⁻¹], P_i is the monthly rainfall [mm] (i = 1, 2, 3, …, 12), and P is the annual rainfall.

(2) Soil erodibility factor (K)

The soil erodibility factor can indicate the sensitivity of soil properties to erosion. In this study, we adopted the erosion-productivity impact calculator model of Williams et al. [56] (Equation (5)), and then we revised the soil erodibility calculated using this method following the study of Zhang et al. [57] (Equation (6)):

$$K\left(EPIC\right)=\left\{0.2+0.3\exp \left[-0.025SAN\left(1-\frac{SIL}{100}\right)\times {\left(\frac{SIL}{CLA+SIL}\right)}^{0.3}\right]\times \left[1.0-\frac{0.25C}{C+\exp \left(3.72-2.95C\right)}\right]\times \left[1.0-\frac{0.7SNI}{SNI+\exp \left(-5.51+22.9SNI\right)}\right]\right\}$$

(5)

$$K=\left(-0.0183+0.51575K_{EPIC}\right)\times 0.1317$$

(6)

where K is the revised soil erodibility factor, SAN is the content of soil sand (particle size: 0.1–2.0 mm, %), SIL is the content of soil silt (particle size: 0.002–0.1 mm, %), CLA is the content of soil clay (particle size < 0.002 mm, %), C is the content of soil organic carbon (%), and SNI = 1 − SAN/100.

(3) Topographical factor (LS)

The topographical factor, which is the combination of the two factors of slope length and slope steepness, represents the influence of topography on soil loss. In this study, the calculation method proposed by Wischmeier et al. [30] and McCool et al. [34] was used to calculate the topographical factor (Equations (7)–(10)):

$$L={\left(\frac{\gamma }{22.13}\right)}^{\alpha }$$

(7)

$$\alpha =\frac{\beta }{\beta +1}$$

(8)

$$\beta =\frac{\sin \theta /0.0896}{3.0\times {\left(\sin \theta \right)}^{0.8}+0.56}$$

(9)

$$S=\left\{\begin{array}{ccc}10.8\times \sin \theta +0.03& & \theta <5\%\\ 16.8\times \sin \theta -0.50& & 5\%\le \theta \le 10\%\\ 21.91\times \sin \theta -0.96& & \theta >10\%\end{array}\right.$$

(10)

where L is the slope length factor, S is the slope factor, $\gamma$ is the slope length (m), $\theta$ is the slope (%), $\alpha$ is the variable slope length index, and $\beta$ is the coefficient that varies with slope.

(4) Vegetation cover factor (C)

This factor is directly related to vegetation coverage and its value is in the range of 0–1. In this study, we adopted the method for the calculation of the vegetation cover factor proposed by Cai [58] (Equations (11) and (12)):

$$C=\left\{\begin{array}{ccc}1& & VFC=0\\ 0.6508-0.3436\lg f& & 0<VFC\le 78.3\%\\ 0& & VFC>78.3\%\end{array}\right.$$

(11)

$$VFC=\frac{\left(NDVI-NDV{I}_{soil}\right)}{\left(NDV{I}_{\max }-NDV{I}_{soil}\right)}$$

(12)

where NDVI_soil is the NDVI value of bare soil pixels and NDVI_max is the NDVI value of pure vegetation pixels. Based on the obtained vegetation index, we took the confidence intervals of NDVI of 5% and 95% as pure bare soil pixels and pure vegetation pixels, respectively, and we replaced NDVI values lower (higher) than the pure bare soil (vegetation) pixel value with the minimum (maximum) value.

(5) Soil and water conservation measure factor (P)

The soil and water conservation measure factor (P) refers to the ratio of the soil loss after adopting certain soil and water conservation measures to the soil loss without taking any measures. Following a study in Guizhou Province by Wang et al. [59], we assigned certain values to different land use types (

Table 2. P-factor values assigned for different land use types on the Guizhou Plateau.

Land Use	Paddyland	Dryland	Woodland	Grassland	Waters	Construction Land	Unused Land
P	0.15	0.5	1	1	0	0	0

).

2.3.2. Correlational Analyses

Correlation analysis is a statistical method for studying the correlation between two or more variables. We established a linear regression model to fit the multiyear trends (2000–2018) of precipitation, vegetation coverage, and SCS, respectively. The degree of correlation between pairs of variables was represented by the Pearson correlation coefficient [60], and the significance level was expressed by the t-test statistic (Equations (13) and (14)):

$$R=\frac{{\displaystyle \sum _{i=1}^n}\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{\displaystyle \sum _{i=1}^n}{\left(X_i-\overline{X}\right)}^2{\displaystyle \sum _{i=1}^n}{\left(Y_i-\overline{Y}\right)}^2}}$$

(13)

$$t=\frac{R\sqrt{n-2}}{\sqrt{1-R^2}}$$

(14)

where R is the Pearson correlation coefficient, X and Y are the two variables, $\overline{X}$ and $\overline{Y}$ are the averages of the two variables, respectively, and i is the i-th year. With the calculated t-test statistic, the corresponding p-value was used to determine the significance level between the variables. In this study, the NDVI, annual rainfall, and SCS on the Guizhou Plateau during 2000–2018 were selected as variables for the correlation analysis, which was programmed using the MATLAB package.

2.3.3. Residual Analysis

In this study, we considered that ecological engineering would change the vegetation coverage. Therefore, residual analysis was used to eliminate the climatic factors of temperature and precipitation in the long-term time series variation of the NDVI. As a series of large-scale ecological engineering projects commenced in 2004, we chose 2004 as the time node. The NDVI, temperature, and precipitation during 2000–2004 were introduced into the regression model to establish a per-pixel NDVI–climate regression model. Then, temperature and precipitation data for 2000–2018 were inserted into the model to obtain predicted NDVI values under the prevailing meteorological conditions without the implementation of ecological engineering projects [39] (Equation (15)):

$${{\displaystyle NDVI}}_{(i,t)}{}^{=}{}^a{{\displaystyle T}}_{(i,t)}{}^{+}{}^b{{\displaystyle P}}_{(i,t)}{}^{+}{}^c$$

(15)

where NDVI_(i,t), T_(i,t), and P_(i,t) represent the NDVI value, temperature, and precipitation of the i-th pixel in the t-th year, respectively; a and b are the regression coefficients of NDVI_(i,t) to T_(i,t) and P_(i,t), respectively; and c is a constant.

The SCS, NDVI, and precipitation for 2000–2004 were introduced into the regression model to establish a per-pixel SCS regression model. Then, the predicted values of the NDVI and precipitation for 2000–2018 were inserted into the model to obtain the predicted values of SCS under the prevailing meteorological conditions without the implementation of ecological engineering projects (Equation (16)):

$${{\displaystyle SCS}}_{(i,t)}{}^{=}{}^{a^{\prime }}{{\displaystyle NDVI}}_{(i,t)}{}^{+}{}^{b^{\prime }}{{\displaystyle P}}_{(i,t)}{}^{+}{}^{c^{\prime }}$$

(16)

where SCS_(i,t), NDVI_(i,t), and P_(i,t) are the SCS, temperature, and precipitation of the i-th pixel in the t-th year, respectively; a and b are the regression coefficients of SCS_(i,t) to NDVI_(i,t) and P_(i,t), respectively; and c is a constant.

The predicted SCS under the prevailing meteorological conditions without the implementation of ecological engineering projects was obtained based on the above two formulas. The actual SCS obtained using the RUSLE model was subtracted from the predicted value to obtain the residual sequence representing the effect of anthropogenic activities [61]:

$${{\displaystyle \epsilon }}_{(i,t)}={{\displaystyle SCS}}_{real(i,t)}-{{\displaystyle SCS}}_{predict(i,t)}$$

(17)

where $\epsilon$_(i,t), SCS_real(i,t), and SCS_predict(i,t) represent the residual sequence of anthropogenic activities [t·ha⁻¹·yr⁻¹], the SCS calculated using the RUSLE model [t·ha⁻¹·yr⁻¹], and the predicted SCS [t·ha⁻¹·yr⁻¹], respectively. A value of $\epsilon$ of >0 (<0) was considered to reflect a positive (negative) effect of ecological engineering on SCS; a value of $\epsilon$ ≈ 0 was considered to reflect weak effect. The spatial distribution of the intensity of the effect of ecological engineering projects on SCS was obtained based on the average value of the residual sequence of 2000–2018, following which the Mann–Kendall test was used to determine the significance of the trend. The z value of the statistic is a standard normal distribution, generally taken as α = 0.05; when |z| > 1. 96, the confidence level of the study sequence α is <0.05, and when |z| < 1.96, the confidence level α is >0.05 [62]. The Mann–Kendall test [63] was used to determine the temporal trend of the residual sequence value and to obtain the significant distribution of the change of the SCS function under the influence of ecological engineering; thus, making the results more credible scientifically.

3. Results

3.1. Spatial and Temporal Characteristics of SCS

3.1.1. Spatiotemporal Variation of SCS

The R, K, LS, C, and P factors used in the RUSLE model were simulated using GIS software. The R, C, and P factors change annually, whereas the K and LS factors are fixed (Figure 5). The average SCS on the Guizhou Plateau during 2000–2018 was 50.67 t ha⁻¹ yr⁻¹; the lowest (32.09 t·ha⁻¹·yr⁻¹) and highest (87.06 t·ha⁻¹·yr⁻¹) values occurred in 2013 and 2017, respectively. As shown in Figure 6, SCS on the Guizhou Plateau showed an overall upward trend during 2000–2018 at a rate of 1.39 t·ha⁻¹·yr⁻¹, significant at the p = 0.05 level, which indicates that the SCS function improved.

Regionally, based on the annual SCS, the area of the Guizhou Plateau with the strongest SCS was in the southeast, where non-karst landforms and peak-cluster depressions form the primary landform types. This area comprises a forest ecosystem with high vegetation coverage, low anthropogenic activity, and an excellent ecological environment. The areas with the weakest SCS function were in central and western parts of the Guizhou Plateau, characterized by karst canyons and karst plateaus with serious soil and water losses. Central parts of Guizhou Province are characterized by many cities with high levels of socioeconomic development and widespread anthropogenic activities. The areas of greatest SCS improvement (>50 t·ha⁻¹·yr⁻¹) within the study area were primarily in the southwest, whereas the regions with the greatest deterioration were mainly in the north (Figure 7).

3.1.2. Spatiotemporal Variation of SCS in Different Landform Regions

In terms of landform region, the highest SCS was 73.12 t·ha⁻¹·yr⁻¹ in the peak-cluster depression region located in the south of the study area, followed by the non-karst landform regions (59.69 t·ha⁻¹·yr⁻¹); the lowest SCS (40.22 t·ha⁻¹·yr⁻¹) was in the central karst plateau region (

Table 3. Average SCS in different landform regions on the Guizhou Plateau.

SCS t·ha−1·yr−1	Karst Canyon	Peak-Cluster Depression	Fault Depression Basin	Karst Plateau	Karst Trough Valley	Non-Karst Landform
Value	49.59	73.12	59.19	40.22	49.38	59.69

). The SCS of all landform regions showed an upward trend; the most rapid improvement (2.34 t·ha⁻¹·yr⁻¹) was in the karst canyon region and the slowest rate of improvement (0.65 t·ha⁻¹·yr⁻¹) was in the non-karst landform regions. The karst canyon and fault depression basin regions in the high-elevation area of the western Guizhou Plateau have long had fragile ecological environments with serious soil and water losses. Conversely, the non-karst region in the east of the study area has lower elevation, relatively flat terrain, and high forest coverage, and its historical ecological background is better than that of the other landform regions (Figure 8).

Among the different landform regions, comparison of the 2000 and 2018 data reveals that the geomorphic areas in which SCS has increased included the karst canyon, fault depression basin, and peak-cluster depression regions. Conversely, SCS has reduced in the non-karst landform and karst trough valley regions and remained largely stable in the karst plateau region. In terms of the total amount, the same conclusion can be reached (Figure 9), i.e., the SCS capacity of most karst areas increased, whereas that of non-karst areas declined.

3.1.3. Spatiotemporal Variation of SCS in Different Ecological Engineering Areas

In the different ecological engineering areas, the highest proportion of the total area is associated with the NFPP (7.20%), followed by the GGP (3.01%); the smallest area is associated with the KPWFP. Among the different ecological engineering areas, the WCNR has the highest SCS (77.44 t·ha⁻¹·yr⁻¹) and its implementation area accounts for 0.71% of the total area. The NFPP has the lowest SCS (42.66 t·ha⁻¹·yr⁻¹) and its implementation area accounts for 7.2% of the total (

Table 4. Area and SCS of different ecological engineering projects on the Guizhou Plateau.

Eco-Engineering	Area km2	Proportion of the Total Area %	SCS t·ha−1·yr−1	Total Amount of SCS 104 t
NFPP	12,688.50	7.20	42.66	5412.82
WCNR	1246.31	0.71	77.44	965.17
KPWFP	932.00	0.53	61.74	575.37
SFP	2634.94	1.50	53.7	1414.87
CCPRD	3964.75	2.25	52.65	2087.34
GGP	5296.38	3.01	45.37	2402.94

). The SCS of the WCNR is much higher than that of the other ecological engineering areas. The project area is focused primarily on the nature reserve, which has a reasonable ecological environment and abundant vegetation layers; however, there are few nature reserves on the Guizhou Plateau. The NFPP has the highest total amount of SCS (much higher than other ecological engineering projects) and its area of implementation is largest, followed by the GGP (Figure 10).

The SCS of the different ecological engineering projects in the different landform regions (2000–2018) was analyzed; however, there was no WCNR in the fault depression basin region. The NFPP, KPWFP, SFP, CCPRD, and GGP in the karst canyon region all showed the highest upward trend of SCS: 3.27, 2.94, 2.62, 2.23, and 2.48 t·ha⁻¹·yr⁻¹, respectively. The SCS of the WCNR improved rapidly in the karst plateau region, which is a major landform of the Guizhou Plateau, characterized by plain terrain and deep soil layers that make it a suitable habitat for the survival of wildlife. The CCPRD improved the SCS function considerably in the various karst landform areas. This reflects the various control measures implemented to improve the damaged natural ecosystem in areas of rocky desertification, areas, which have had obvious effect on the improvement of soil and water conservation in karst areas. The SCS of each ecological engineering area improved slowly in non-karst landform areas. Overall, the construction of ecological engineering projects has had marked effect on SCS in karst landform areas but not so in non-karst landform areas (Figure 11).

3.2. Analysis of Driving Mechanism of SCS

3.2.1. Driving Effect of Precipitation and VFC on SCS

SCS is affected by various factors that can be divided into natural factors and anthropogenic activities. Among the former, precipitation and VFC are the most important factors affecting soil erosion [64]. The Pearson method was used to calculate the degree of correlation of both precipitation and VFC with SCS during 2000–2018, and the t-test statistic was calculated to assess the significance of the results.

During 2000–2018, areas with significant positive correlation between the annual variations of precipitation and SCS, which accounted for 59.43% of the total area, were located in northern and northwestern parts of the Guizhou Plateau (p < 0.05) (Figure 12a). In addition, areas in which VFC also had significant positive correlation with SCS (accounting for 47.26% of the total area) were distributed widely throughout the entire study area, but concentrated primarily in western parts (p < 0.05) (Figure 12b). Neither precipitation nor VFC showed significant negative correlation with SCS on the Guizhou Plateau. Areas of significant positive correlation of SCS with precipitation were mainly in northern parts of the karst trough valley region as well as some parts of the karst canyon and karst plateau regions. Moreover, VFC mainly affected SCS in the southwest of the study area, which includes the karst canyon, fault depression basin, and peak-cluster depression regions.

In the different ecological engineering areas, there was only positive correlation between precipitation, VFC, and SCS. The highest proportion of positive correlation between precipitation and SCS was in association with the WCNR (76.59%) and the lowest was in association with the KPWFP (41.37%). For VFC, the GGP and the WCNR had the highest and lowest proportion of positive correlation with SCS, respectively, i.e., 53.02% and 12.81%, respectively (Figure 13). A probable explanation is that the WCNR is mainly a natural reserve with high vegetation coverage, in which excessive rainfall could lead to intensified potential soil erosion that might increase SCS. Moreover, the GGP directly converts large areas of cultivated land with severe soil erosion and rocky desertification into forest, restoring the vegetation, directly controlling soil erosion, and improving the SCS function. In each ecological engineering area, both precipitation and VFC were positively correlated or uncorrelated with SCS, i.e., neither exhibited any negative correlation.

3.2.2. Driving Effect of Anthropogenic Activities on SCS

We predicted the NDVI based on temperature and precipitation data without the effect of ecological engineering implementation, and then used it to simulate the predicted value of SCS. This could exclude the impact of vegetation change caused by ecological engineering implementation. The effect of anthropogenic activities was obtained through residual sequence analysis. The average R² of the modeled SCS was 0.61, and 16.1% of pixels passed the 0.05 significance level test. In addition, SCS predicted by the regression models and the real SCS exhibited high correlation (0.80) in terms of the spatial distribution, reaching significant correlation at the p < 0.05 level. Given that the time series used for the regression was short and that the study area has very complex terrain and climatic conditions, it was considered suitable for achieving the purposes of this study.

In this study, the average annual residual value of the study area was 5.71 during 2000–2018. According to the spatial distribution of the average annual residual value during 2000–2018 (Figure 14a), it can be seen that anthropogenic activities have had a primarily positive effect on SCS. The high-value regions of this positive effect were located in northeastern and western parts of the study area, i.e., primarily in karst trough valley, karst canyon, fault depression basin, and peak-cluster depression regions. Conversely, the areas of negative effect were distributed mostly in southeastern and northern parts of the study area, i.e., mainly non-karst landform regions. According to the Mann–Kendall significance trend test (Figure 14b), the areas in which the trend of the effect of anthropogenic activities was significantly positive (negative) accounted for 26.43% (3.97%) of the entire study area (p < 0.05). Overall, the positive effect of anthropogenic activities on SCS has been stronger than the negative effect. It can be concluded that anthropogenic activities such as ecological engineering construction have had positive effects in karst landform areas and negative effects in non-karst landform areas.

The distribution of ecological engineering projects in the study area (Figure 15) shows that the NFP is distributed widely in various ecological engineering areas, the GGP is distributed mainly in karst landform areas, and the WCNR is distributed primarily in non-karst landform areas. It can be seen from the spatial distribution of ecological engineering areas that the areas in which the effect of anthropogenic activity is positive (i.e., where the residual sequence is >0) were generally distributed in areas in which ecological engineering projects are implemented, especially in the karst landform regions. Conversely, in non-karst landform regions, even if certain ecological engineering projects were implemented (e.g., the WCNR), the effect of anthropogenic activities was negative, whereas in areas without implementation of ecological engineering projects, the effect of anthropogenic activities had a significant negative trend. It can be concluded that particular attention should be paid to soil and water losses as well as to the implementation of ecological engineering projects in non-karst landform regions.

The average residuals of the different ecological engineering projects during 2000–2018 are listed in

Table 5. Average residuals of the different ecological engineering projects during 2000–2018 on the Guizhou Plateau.

Eco-Engineering	NFPP	WCNR	KPWFP	SFP	CCPRD	GGP	OFP
Average annual residual	3.4	−3.11	4.46	10.21	6.05	9.45	7.04

. The SFP has the highest value (10.21) because this ecological engineering project has the greatest positive effect on the SCS function. The WCNR is the only ecological engineering project to have a negative residual value (−3.11). It can be concluded that ecological engineering construction has generally had positive effect on the SCS function, although certain anthropogenic activities in the WCNR could have had negative impact on the ecosystem that might have led to soil and water losses.

For the different ecological engineering projects, the areas of positive and negative effects of anthropogenic activities have different proportions. The ecological engineering project with the highest proportion of positive effects associated with anthropogenic activities is the SFP (82.06%). Moreover, the areas with positive effects associated with anthropogenic activities accounted for >60% of the total area for the GGP, CCPRD, and NFPP. Conversely, the area of the WCNR that exhibited negative effects associated with anthropogenic activities accounted for 68.44% of the total area. The WCNR was the only ecological engineering project to have a negative area that accounted for >50% of its total area (Figure 16).

Comparison of the change in VFC during 2000–2018 reveals that vegetation has generally improved, especially in karst canyon, fault depression basin, and karst trough valley regions (Figure 17). In considering the impact of anthropogenic activities on SCS, the most important factor is the coupling relationship between vegetation change and anthropogenic activity. Ecological engineering construction changes the vegetation coverage and affects SCS indirectly. Following the implementation of the various ecological engineering projects in the study area, vegetation coverage increased considerably in all karst landform regions, e.g., the VFC in the karst canyon region was 0.72 in 2000 and 0.80 in 2018 (

Table 6. The VFC associated with different landform regions and different ecological engineering projects in 2000 and 2018.

Landform	Karst Plateau	Karst Canyon	Fault Depression Basin	Peak-Cluster Depression	Karst Trough Valley	Non-Karst Landform
2000	0.75	0.72	0.74	0.77	0.77	0.81
2018	0.80	0.79	0.80	0.82	0.81	0.84
Eco-Engineering	NFPP	WCNR	GGP	SFP	CCPRD	KPWFP
2000	0.78	0.84	0.76	0.78	0.78	0.79
2018	0.83	0.83	0.84	0.83	0.83	0.82

). Conversely, vegetation coverage increased only slightly in the non-karst landform regions (0.81 in 2000; 0.84 in 2018). In addition, it is evident that VFC increased markedly following the implementation of all the ecological engineering projects except the WCNR. The VFC associated with the GGP was 0.76 in 2000 and 0.84 in 2018, which indicates that the implementation of GGP is more appropriate in specific areas. In contrast, VFC in the WCNR was 0.84 in 2000 and 0.83 in 2018, i.e., a slight reduction.

4. Discussion

4.1. Model Simulations of Soil Erosion Modules and SCS

In this study, the R, K, LS, C, and P factors used in the RUSLE model were simulated using local rainfall data, RS data, and regional ecological engineering data. All factors were processed into raster data with 250-m resolution and used to simulate the soil erosion modulus as well as SCS through parameter localization. The simulation of the R factor in this study was based on interpolation of rainfall data using ANUSPLIN (spatial resolution: 250 m). According to Yue et al. [65], the average absolute error of the interpolation results in each grid is changed only slightly when the size of the grid is increased from 50 to 1000 m. This is mainly related to the scale of the spatial distribution of rainfall, which usually varies little within the range of 1 km². The factor values of L and S obtained from DEMs with different spatial resolution are notably different [66], Generally, topography tends to be flattened in small-scale DEM data and the slope decreased, which means the intensity of soil erosion could be underestimated [67], whereas the slope value of a real surface cannot be reflected directly by large-scale DEM data [68]. Given the difficulty in obtaining high-resolution DEM data and recognizing the complex topography of our study area, we chose 90-m resolution DEM data. Such data, which are considered suitable for large-scale spatial analysis, made the calculation of the topographical factor more accurate and were appropriate for the karst landform and terrain fragmentation in the study area. The value of the P factor reflects the characteristics of different soils and water conservation measures, e.g., crop rotation, ground cover, and contour tillage. Generally, values of the P factor are obtained based on RS imagery and assigned according to land use type on the large scale [69]. We used land use data of 2000–2015 to assign the values of the P factor, which could well reflect the impact of land use and land cover changes on soil erosion. The original 1-km resolution data were changed to 250-m resolution by resampling because of the difficulty associated with obtaining high-resolution land use data. To improve the accuracy of the results, water bodies and construction land were directly assigned a value of zero. Despite the many uncertainties involved in the acquisition and determination of the R, K, L, S, C, and P factors, and the difficulty in assessing their accuracy [68], we consider that the selected data had sufficient precision to ensure the validity of the final result.

During our study, the annual average actual soil erosion modulus on the Guizhou Plateau during 2000–2018 was 10.51 t·ha⁻¹·yr⁻¹, and the soil erosion modulus calculated using the RUSLE model in 2000 was 18.84 t·ha⁻¹·yr⁻¹. This is in reasonable agreement with the average annual soil modulus of the Proclamation of Soil and Water Loss in Guizhou Province in 2000 (14.32 t·ha⁻¹·yr⁻¹ issued by the Water Resources Department of Guizhou, China [70]. In addition, using the RUSLE model, both Gao et al. [71] and Zeng et al. [72] calculated the soil erosion modulus of certain areas of the Guizhou Plateau as 15.08 and 14.40 t·ha⁻¹·yr⁻¹, respectively. Thus, given this level of accord, we believe that our evaluation of SCS using the RUSLE model and the data described in Section 2.2 had satisfactory accuracy.

In this study, we used the RUSLE model to investigate the sensitivity of SCS on the annual timescale. However, soil erosion has strong seasonal variation owing to the influence of natural conditions [73]. Based on a study in Danjiangkou (China), Feng et al. [74] reported that soil erosion in summer and autumn was higher than in spring and winter, the main reason for which was the uneven distribution of rainfall throughout the year. The plant canopy is another important factor influencing soil and water conservation. Wu et al. [75] found that crops with marked seasonal changes of characteristics throughout the year had notable influence on soil loss on the Loess Plateau (China). If monthly scale data could be obtained, the relative impacts of vegetation and precipitation could be reflected more clearly.

4.2. Characteristics of SCS

In the Guizhou Plateau region, SCS has generally improved during 2000–2018 with an upward trend of 1.39 t·ha⁻¹·yr⁻¹. Analysis of the change of VFC during 2000–2018 revealed an upward trend of VFC in the Guizhou Plateau region. The average annual VFC during 2000–2004 (2014–2018) was 0.77 (0.82). Moreover, the average VFC was >0.8 in each year during 2014–2018. According to Gao et al. [71], when vegetation coverage is >0.5–0.6, soil erosion decreases with a further increase in vegetation coverage. It is widely believed that vegetation restoration could improve the SCS capacity substantially [76]. The findings of this study corroborate the results of previous related research. In Guizhou Province, the experimental work of the GGP and NFPP began in 2000. In 2001, the SFP, WCNR, and KPWFP were initiated in the Yangtze River Basin. Since 2004, large-scale implementation of the CCPRD has been undertaken, and various other ecological engineering projects have been performed. The ecological engineering measures such as the GGP have been shown to increase vegetation coverage and enhance the soil and water conservation capacity [77].

On the Guizhou Plateau, SCS was found to show high spatial heterogeneity. Geomorphology is known to affect the soil erosion rate [78]. Areas with high SCS were found mainly in southeastern parts of the Guizhou Plateau, i.e., in non-karst landform and peak-cluster depression regions. Areas with lower SCS were found primarily in central and western areas, i.e., the karst plateau and karst canyon regions. Comparison of non-karst landform and karst landform areas revealed that the underlying surface of the former is more suitable for vegetation growth [44] and more conducive to vegetation recovery. In contrast, the decline in the rate of soil erosion is faster in non-karst landform regions than in karst regions that have terrain fragmentation and a fragile ecological environment [72]. Moreover, the central karst plateau region is mainly an urban ecosystem with notable urban construction activities, high population density, and considerable anthropogenic interference [79]. In the western karst region of the Guizhou Plateau, soil erosion has a complex relationship with topography, lithology, and the degree of rocky desertification [80,81].

The areas with increased SCS were found mainly in the western karst canyon and the southern peak-cluster depression regions, of which the latter exhibit the most obvious improvement. In comparison with other areas of the Guizhou Plateau, the hydrothermal conditions of southern parts are more advantageous for vegetation growth and recovery [39]. Following the implementation of many ecological engineering projects, especially the GGP that has returned areas of steep farmland to forest, SCS on the Guizhou Plateau has improved significantly. The karst canyon region in the west of the study area has a history of intense tectonic activity that has led to widespread distribution of highly eroded carbonate rocks. The regional topography is complex and the ecological environment fragile [82]. Following the implementation of ecological engineering projects such as the KPWFP and NFPP in these areas, which expanded the area of forest to protect the existing vegetation, the SCS function was improved. However, improvement of the SCS function in non-karst landform regions was not obvious; in some cases, it actually exhibited a downward trend. Historically, the ecological background of non-karst landform regions is better than that of other areas and the forest coverage rate is always high. Forest coverage can improve the soil structure and regulate runoff, which reduces the influence of rainfall [83]. Thus, the better underlying ecological conditions of non-karst landform areas meant that there was little scope for improvement in the SCS function by the implementation of ecological engineering projects. However, with continued socioeconomic development, it will become increasingly difficult to maintain the ecological environment of such areas.

4.3. Driving Mechanism of SCS Under the Effects of Ecological Engineering

The mechanism of soil erosion, which is highly complex and influenced by many factors such as rainfall, soil type, terrain, and vegetation, results in different erosion characteristics in different regions [84]. Often, SCS is influenced by the coupling of multiple factors that are difficult to distinguish. In this study, the highest (lowest) value of SCS during 2000–2018 was in 2017 (2013). Calculation using the RUSLE model showed that the actual soil erosion modulus was similar in 2013 and 2017, i.e., 6.32 and 7.04 t·ha⁻¹·yr⁻¹, respectively. However, the average annual rainfall in 2013 was less than in the other years, and the value of the rainfall erosivity factor was only 2828.87 MJ·mm·ha⁻¹·h⁻¹·yr⁻¹, resulting in less potential soil erosion that led to the lower SCS. In contrast, the average annual rainfall in 2017 was the highest of the 19 years (1248 mm) and the value of the rainfall erosivity factor was as high as 7266.64 MJ·mm·ha⁻¹·h⁻¹·yr⁻¹. These results showed that rainfall also has considerable impact on SCS. However, the vegetation coverage in 2017 (average VFC: 0.83) was also better than in 2013 (average VFC: 0.79). Thus, although the substantial rainfall increased the potential soil erosion, the higher vegetation coverage resulted in the higher SCS. These results are similar to those reported by Xiao et al. [52] in relation to the Three Gorges region (China).

In this study, precipitation and VFC were found mainly correlated positively with SCS on the Guizhou Plateau. The average VFC on the Guizhou Plateau was 0.76 in 2000 and 0.81 in 2018. Vegetation plays an important role in SCS in the Guizhou Plateau area. Of the different ecological engineering projects, the WCNR exhibited the highest proportion of significant positive correlation between precipitation and SCS. Although rainfall can increase the potential soil erosion, adequate precipitation also promotes vegetation growth, leading to strengthening of the SCS function. Ecological policies and projects such as the GGP, which aim to regulate anthropogenic activity and influence land use, will ultimately affect the condition of regional vegetation [85].

Our study showed that anthropogenic activities on the Guizhou Plateau have generally had positive effect on SCS, and that the areas in which the positive effect has been greatest are distributed primarily in northeastern and western parts of the study area. The original ecological environmental quality is poor in karst landform areas that have complex physiognomy such as peak-cluster depressions, karst canyons, and fault depression basins, and almost all these areas reflect the positive effect of anthropogenic activities. Generally, it is in such areas that many ecological engineering projects have been implemented, e.g., the SFP and GGP. Thus, the anthropogenic activities associated with the implementation of ecological engineering constructions in karst areas could account for the positive effect on SCS. In addition, following the implementation of ecological engineering projects, vegetation coverage in all karst landform areas has been increased markedly. However, in southeastern parts of the Guizhou Plateau, which have a better historical ecological background and less regional ecological engineering construction, anthropogenic activities showed significant negative effects on SCS. As shown in Table 6, vegetation coverage has increased only slightly in non-karst landform areas, which might reflect overexploitation in terms of tourism and forestry resources in areas with a robust ecological environment [39].

The area exhibiting positive effect of anthropogenic activities in relation to ecological engineering projects was highest for the SFP. In the Guizhou Plateau region, the SFP comprises mainly the Yangtze River protection forest system and the Pearl River protection forest system, the aim of which is to maintain the core area of the mountains around the basins, protect the main natural forest resource, and improve soil and water conservation. This project, which can improve the hydrological condition of a watershed by slowing surface runoff and preventing soil erosion, relies on anthropogenic activity to improve the SCS function. It can be seen from Table 6 that vegetation increased markedly following the implementation of the SFP. Conversely, in some areas with a better historical ecological environment (e.g., the WCNR), vegetation was not maintained or was even destroyed by anthropogenic activity. Therefore, it would be better to give priority to natural development and recovery in such areas and to reduce anthropogenic interference to a minimum. These results will be useful in relation to controlling soil loss and increasing awareness of ecological environmental issues.

4.4. Deficiencies and Prospects

The objective of this study was to investigate the impact of ecological engineering on the SCS function of the Guizhou Plateau during 2000–2018. Although the RUSLE model is a mature model used widely in studies on soil erosion, it tends to overestimate the magnitude of soil erosion in karst areas [86]. There has been insufficient study of the topography of karst regions and therefore the complexity of the internal characteristics of each karst landform region should be investigated in detail [87]. Moreover, the model used requires further improvement in terms of its accuracy. In our future work, different equation coefficients could be modified to reflect different landform regions to ensure scientific validity and to improve precision.

Precipitation on the Guizhou Plateau showed no significant trend of change during 2000–2018. However, some studies have shown that precipitation could affect soil and water losses [88]. The length of the time series data used in this study might be too short to reflect the effects of climate change, and the span of the time series could make it difficult to estimate the driving effects of precipitation. Furthermore, owing to the limitations of data acquisition, only five year-length datasets were used for the climate simulation prior to the implementation of the ecological engineering projects. In future research, data from before 2000 should be incorporated to increase the length of the time series and to make the conclusions more persuasive. Thus, future study on the impacts of climate change and vegetation restoration on SCS could be improved. In addition, the use of the two regression equations could have increased the error of the simulation, and the short time series could have caused the model to become unstable. We did not consider that the NDVI might present a certain temporal lag with respect to precipitation, temperature, and soil moisture [89]. Using annual precipitation as a variable with which to analyze the relationship between precipitation and the NDVI meant that the statistical results would be affected adversely by invalid precipitation data.

Vegetation contributes greatly to the ecosystem services of the Guizhou Plateau, and anthropogenic activities such as ecological engineering play an important role in SCS. Hence, obvious changes in vegetation restoration in certain areas could be attributed to the effects of ecological engineering constructions. However, in areas with insignificant changes in vegetation, especially over long periods, the role of anthropogenic activities in relation to SCS requires further investigation. In future research, the contributions of the effects of climate change and anthropogenic activities should be examined further.

5. Conclusions

This study evaluated SCS using the RUSLE model and analyzed the spatial and temporal variations of SCS in the Guizhou Plateau region (2000–2018) based on multiyear meteorological and RS data. During 2000–2018, the average annual SCS of the Guizhou Plateau showed an obvious upward trend at a rate of 1.39 t·ha⁻¹·yr⁻¹. Areas with high (low) SCS values were located mainly in southeastern non-karst landform and southern peak-cluster depression (central karst plateau and western karst canyon) regions. Moreover, the area with the most obvious improvement of SCS (2.34 t·ha⁻¹·yr⁻¹) was in the karst canyon region, whereas the areas with the worst improvement (0.64 t·ha⁻¹·yr⁻¹) were in the non-karst landform regions. Ecological engineering construction has had evident positive effect on SCS, especially in the karst area, whereas the effect is not obvious in non-karst landform regions.

This study also found that precipitation and VFC were largely correlated positively with SCS on the Guizhou Plateau, and that the areas with significant positive correlation accounted for 59.43% and 47.26%, respectively; areas with negative correlation had almost no distribution. The WCNR had the highest proportion of significant positive correlation between precipitation and SCS, i.e., up to 76.59%, and the GGP had the highest proportion of significant positive correlation between VFC and SCS, i.e., up to 53.02%. According to the Mann–Kendall test, the areas in which the trend of the effect of anthropogenic activities was significantly positive (negative) accounted for 26.43% (3.97%) of the entire study area. As for the different ecological engineering projects, the greatest area of the positive effect of anthropogenic activities was associated with the SFP, i.e., 82.06% of the total construction area. Moreover, the proportion of the positive effect of anthropogenic activities of the GGP, CCPRD, and NFPP were all >60%. Overall, the implementation of ecological engineering projects has markedly improved the SCS function on the Guizhou Plateau. However, the negative effect of anthropogenic activities in the WCNR accounts for >50% of the total construction area. For nature reserves with a historically better ecological background, the adverse effects of anthropogenic activities should be identified and then corrected and improved. In such areas, to prevent soil and water losses, the focus should be on natural evolution and restoration, reduced tourism development, and minimal anthropogenic interference.