The Effect of Land Use/Cover Change on Soil Erosion Change by Spatial Regression in Changwu County on the Loess Plateau in China

Abstract

Changwu County is a typical soil and water loss area on the Loess Plateau. Soil erosion is an important ecological process, and the impact of land use/cover change on soil erosion has received much attention. The present study used remote sensing images of the study area in 1987, 1997, 2007, and 2017 to analyze the land use/cover change (LULCC), and the RUSLE model was applied to estimate the soil erosion in different times. We exploited the Sankey diagram to visualize the spatiotemporal changes in land use/cover and soil erosion. We planned to obtain the most suitable model by comparing the application of different spatial regression models (Geographically weighted regression model, Spatial lag model, Spatial error model) and Ordinary least squares in LULCC and soil erosion changes. The results revealed that land use/cover has significantly changed in the last 30 years. From 1987 to 1997, cropland expansion came mainly from planted land and orchards, which transformed 68.99 km² and 64.93 km², respectively. In 1997–2007, the planted land increase was mainly through the conversion of cropland. In 2007–2017, the increase in orchard area came mainly from cropland. The forest land increase was mainly from the planted land. Soil erosion in Changwu County was dominated by slight erosion and light erosion, although the area of slight erosion and light erosion continued to decrease. The annual average soil erosion increased, which was estimated at 977.84 ton km⁻² year⁻¹, 1305.17 ton km⁻² year⁻¹, 1310.60 ton km⁻² year⁻¹, and 1891.46 ton km⁻² year⁻¹ in 1987, 1997, 2007, and 2017, respectively. These amounts of transformation mainly occurred when slight erosion was converted to light erosion, light erosion was converted to moderate erosion, and moderate erosion was converted to light and severe erosion. The Spatial lag model and Spatial error model have higher accuracy than the Geographically weighted regression model and Ordinary least squares when fitting the effect of LULCC and soil erosion change, where the accuracy exceeded 0.62 in different periods.

1. Introduction

The human understanding of nature has gone through a long process, from plunder to protection. Land use and land cover change (LULCC) reflects the process of humans getting along with nature [1,2,3]. Land use has been defined as the purpose for which humans exploit the land cover. This involves both the manner in which the biophysical attributes of land are manipulated and the intent underlying that manipulation, that is, the purpose for which the land is used. Land cover has been defined by the attributes of the Earth’s land surface and immediate subsurface, including biota, soil, topography, surface- and ground-water, and human (mainly built-up) structures [4]. In recent years, research on LULCC has become more extensive and in-depth, mainly focusing on the driving factors of LULCC (including human activities, population increase, ecological factors, economic development, and government policies) [5,6,7], the effect of LULCC (including on land quality, hydrological cycle, climate change, biodiversity reduction, and agricultural production) [8,9,10,11], and the prediction of future changes in LULCC [12,13].

Soil erosion leads to the degradation of soil fertility, the decline in soil productivity, the destruction of surface vegetation, the deterioration of the environment, the restriction of economic development in the region, and the massive transfer of sediment and pollutants in the erosion process, which causes the eutrophication and siltation of the rivers and lakes in downstream areas [14,15,16,17]. Soil erosion assessments have received extensive attention since soil erosion threatens environmental and socio-economic conditions worldwide [16]. The RUSLE has been widely used in estimating soil erosion in the Loess Plateau. Jin used the RUSLE to estimate soil erosion in the Loess Plateau to study its response to vegetation and climate change [18]. Zhang estimated soil erosion in the largest tableland of the Loess Plateau by RUSLE to study effect of land use changes on soil erosion [19].

Extensively covered by characteristic powdery, dry soil (called loess), and coupled with unreasonable human activities and heavy rain, the Loess Plateau in China has suffered serious soil erosion [20]. Since the founding of the People’s Republic of China, the government has implemented a series of measures to control soil erosion [20,21]. This has achieved significant results after long-term efforts by the Chinese government.

Changwu County is located in the loess gully region of the Loess Plateau [22]. As a typical area of soil erosion on the Loess Plateau, soil and water conservation began in the 1950s, and has undergone local decentralized governance, centralized continuous management, comprehensive watershed management, and county-scale management. In the late 1990s, the implementation of the Grain for Green Project played a key role in suppressing local soil erosion [23]. However, government policies and human activities all had an impact on the land use/cover pattern, and also brought great challenges or chances for soil and water conservation in Changwu County. Research on land use/cover and soil erosion in different areas of the loess plateau has attracted the continued attention of researchers [24,25,26]. However, with the intervention of long-term government policies and human activities, the spatial and temporal changes of the land use/cover pattern and soil erosion in the study area remains unclear and the effect of LULCC on soil erosion change is not established.

LULCC and soil erosion changes have obvious spatial attributes, which cannot meet the requirement of independent variables in classical statistics. So as to more accurately study the effect of LULCC on soil erosion changes, we chose different spatial regression models (Geographically weighted regression model, Spatial lag model, Spatial error model) to fit the relationship between the LULCC on soil erosion changes. We planned to obtain the most suitable model by comparing the application of different spatial regression models in LULCC and soil erosion changes. Ordinary least squares (OLS) is used as the reference of the spatial regression models. The present study aims to analyze the spatial and temporal change process of land use/cover by classifying the four remote sensing images of the study area in 1987, 1997, 2007, and 2017, and analyze the spatial and temporal change process of soil erosion through the quantitative estimation of soil erosion with the use of the Revised Universal Soil Loss Equation (RUSLE) model, and we choose suitable spatial regression models to fit the effect of LULCC on soil erosion changes.

2. Materials and Methods

2.1. Study Region

The study area is situated in Changwu County, Shaanxi Province in the south of the Loess Plateau of China, which extends from a latitude of approximately 34°59′ N–35°18′ N and a longitude of 107°38′ E–107°58′ E, occupying an area of 567.1 km² (Figure 1) [23,27]. The study area belongs to the loess gully region of the Loess Plateau. The “Yuan” (a highland formed by erosion in the Loess Plateau, steep on four sides, flat on top) and “Chuan” (a low-lying, open, and flat area formed by erosion between the Yuan and Yuan) alternately had a gully development, forming crisscross gullies. The landforms were mainly “Yuan” and slope land, followed by “Chuan”. The landform features were characterized by high “Yuan”, deep gully, and steep slope. This area was drained by three main river systems, namely, the Jinghe River, Heihe River and Nanhe River, which are tributaries of the Yellow River. The dominant soil type was dark loessial soil, which is at “Yuan”, and loessial soil, which is at the valley slope and gully channel. The altitude ranged from 847 m to 1274 m above sea level [23,28].

Changwu County belongs to a warm, temperate, semi-humid continental monsoon climate zone, with an average annual temperature of 9.1 °C. Furthermore, the frost-free period is 171 days, and the annual precipitation is 587.8 mm. The annual rainfall concentration period is 7, 8 and 9 months, and rainfall can reach up to 321.4 mm, accounting for 54.9% of the total annual precipitation. The study area was mainly planted forest, and afforestation began in the 1950s. The main tree species are Robinia pseudoacacia, Platycladus orientalis, Pinus tabuliformis, Armeniaca sibirica, and Populus spp. Among these, R. pseudoacacia accounts for more than 90% of the forested area [22,23,27].

The research area has few land resources and has experienced agricultural reclamation for a long time. Agriculture occupies a major position in the county’s national economy. The cropland and orchard are mainly concentrated on the plateau and gully bottom, and the main plants were wheat, apple, and persimmon.

2.2. Data Sources

The land use/cover maps of Changwu County on 28 August 1987, 23 August 1997, 20 September 2007, and 27 June 2017 were derived from Landsat TM/ETM+/OLI imageries with 30 m resolution (derived from the USGS website, http://www.usgs.gov, accessed on 16 July 2019), considering the vegetation growing season and image quality. The digital elevation model (DEM) was sourced from the geospatial data cloud (http://www.gscloud.cn/, accessed on 16 July 2019) with a resolution of 30 m.

Then, each image was radiometrically, atmospherically processed using the ENVI 5.3 image processing software to eliminate the scattering and absorption effects from the atmosphere and to obtain authentic the surface reflectance character. Finally, the images were geometrically corrected and registered to the Xian_1980 datum and Xian_1980_3_Degree_GK_CM_108E projections using the ArcGIS 10.5 image processing software to correct image distortion.

The object-based classification segmented the image into objects and identified both the spectral and spatial features between the objects for classification. Object-based classification reduces the within-class spectral variation, and generally eliminates salt-and-pepper effects to potentially improve the classification accuracy. This approach identifies the image objects using a multi-resolution segmentation approach, as implemented in the software package eCognition [29,30,31]. In the present study, the combination of the threshold and decision-tree approach classifications based on object-oriented classification methods were used by the eCognition Developer 9.0 remote sensing information extraction software.

We interpreted six land use/cover classes from the images (

Table 1. Description of land use/cover types identified in the study area.

Type	Description
Orchard	Agricultural land for planting fruit trees
Cropland	All agricultural land used for growing food crops
Forest land	The trees closed to canopy or have closed canopy
Planted land	Including open forest, bushes, shrubs and grass
Water body	Including river and reservoir
Construction land	Including settlements, roads and buildings

), which were orchard, cropland, forest land, planted land, water body and construction land. A total of 388 points (some inaccessible control points with high resolution images in Google earth) were randomly selected to determine the classification accuracy by computing the overall accuracy and Kappa coefficient. Overall accuracy is the sum of correctly classified values from the diagonal values divided by the total number of randomly generated reference values of the error matrix [32]. The Kappa coefficient, which reflects the consistency between classification results and reference data, was also calculated, as follows [1]:

$$K_c=\frac{N{{\displaystyle \sum }}_{i=1}^r x_{ii}-{{\displaystyle \sum }}_{i=1}^r x_i+x_{+i}}{N^2-{{\displaystyle \sum }}_{i=1}^r{x}_i+x_{+i}}$$

(1)

where: K_c is the Kappa coefficient; N is the total number of verification points; r is the number of type; x_i+ is the sum of the data in row i of the confusion matrix; x_+i is the sum of the data in column i of the confusion matrix; x_ii is the value of row i and column i in confusion matrix.

2.3. LULCC Analysis

After completing the classification and accuracy verification, we used the intersect tool in ArcGIS 10.6 to calculate the transfer matrix of the image classification results of different periods, and counted the area transferred by type.

The transfer matrix can intuitively reflect the transfer trend between different land use/cover types within a certain period, and the model can be expressed, as follows:

$$S_{ij}=\left[\begin{array}{ccc}\begin{array}{cc}{S}_{11}& S_{12}\\ S_{21}& S_{22}\end{array}& \begin{array}{cc}\dots & S_{1j}\\ \dots & S_{2j}\end{array}& \begin{array}{cc}\cdots & S_{1n}\\ \dots & S_{2n}\end{array}\\ \begin{array}{cc}\dots & \dots \\ S_{i1}& S_{i2}\end{array}& \begin{array}{cc}\dots & \dots \\ \dots & S_{ij}\end{array}& \begin{array}{cc}\dots & \dots \\ \dots & S_{in}\end{array}\\ \begin{array}{cc}\dots & \dots \\ S_{n1}& S_{n2}\end{array}& \begin{array}{cc}\dots & \dots \\ \dots & S_{nj}\end{array}& \begin{array}{cc}\dots & \dots \\ \dots & S_{nn}\end{array}\end{array}\right]$$

(2)

where: S_ij is the area of the land use/cover type i at time t converted to the land use/cover type j at time (t + 1); n is the total number of land use types.

$$P_R=\frac{S_{ij}}{{{\displaystyle \sum }}_{j=1}^n S_{ij}}$$

(3)

where: P_R is the land use/cover type transfer rate.

In order to reveal the changes in land use/cover area in different periods, the land use/cover annual area change rate (L_R) was calculated, as follows [33]:

$$L_R=\frac{L_{t+n}-L_t}{L_t}\times \frac{1}{n}\times 100\%$$

(4)

where: L_R is the land use/cover annual area change rate; L_t is the area of the land use land cover class at time t; L_t+n is the area of the same land use land cover class at time t + n; n is the time difference between t and t + n.

In order to visualize the temporal and spatial changes in land use/cover, we finally used the Sankey diagrams to illustrate the transfer situation of land use/cover classes at different times. Sankey diagrams are a type of flow diagram, in which the width of the arrows is proportional to the flow rate. These diagrams have mostly been applied to analyze energy- or material-flows [34]. The Sankey diagram has had few applications in land use/cover change [35,36]. In the Sankey diagrams, the width of the lines represents the transfer areas of land use/cover class. The stacked bars visualize the land use/cover class in 1987, 1997, 2007, and 2017, and the height of each component in the stacked bars is proportional to the amount of area occupied by land use/cover class in the study area. We used the R software package called networkD3 to map the Sankey diagrams [37].

2.4. Soil Erosion Estimation

In the present study, we used the revised universal soil loss equation (RUSLE) model to estimate the mean annual soil erosion in the study area at different times. The RUSLE model has been widely used to estimate soil erosion by considering the climate, topography, soil, and land use factor.

The RUSLE model equation was defined by Renard [38] as follows:

$$\mathrm{A} =100\times \mathrm{R} \times \mathrm{K} \times \mathrm{LS} \times \mathrm{C} \times \mathrm{P}$$

(5)

where: A is the average annual soil erosion rate (t km⁻² year⁻¹); R is the rainfall erosivity factor (MJ mm hm⁻² h⁻¹ year⁻¹); K is the soil erodibility factor (thm² h hm⁻² MJ⁻¹ mm⁻¹); LS is the topographic factor; C is the crop management factor; P is the conservation supporting practice factor.

2.4.1. R, Rainfall Erosivity Factor

Rainfall data was collected from the Loess plateau science data center, National Earth System Science Data Sharing Infrastructure, National Science and Technology Infrastructure of China (http://loess.geodata.cn, accessed on 6 July 2020). The monthly average rainfall data during the period of 1987–2017 was used to calculate the R factor. The study area has no record of rainfall intensity. As a result, the monthly rainfall data was used to calculate the R factor using the following equation developed by Wischmeier and Smith [39].

$$\mathrm{R}={\displaystyle \sum }_{i=1}^{12} (1.735\times {10}^{1.5\times lg\frac{P_i{}^2}{P}-0.08188})$$

(6)

where: P_i is the monthly amount of rainfall and P is the annual rainfall.

2.4.2. K, Soil Erodibility Factor

The K factor was measured under a standard plot, since this represents the susceptibility of soil to erosion, the amount and rate of runoff, and transportability of the sediment. The K factor was obtained from Loess Plateau Data Center, National Earth System Science Data Sharing Infrastructure, National Science and Technology Infrastructure of China (http://loess.geodata.cn, accessed on 6 August 2020).

2.4.3. LS, Slope Length and Steepness Factor

The LS factor reflects the impact of topographic features on soil erosion. We used the formula for the graded slope factor proposed by Foster and Wishcmeier [39], and the formula for the slope length factor proposed by Desmet and Govers [40].

$$\mathrm{S} =10.8 \sin \mathsf{\theta }+0.03\hspace{1em}\mathsf{\theta } <5°$$

(7)

$$\mathrm{S} =16.8 \sin \mathsf{\theta } -0.5\hspace{1em}5° \le \mathsf{\theta } <10°$$

(8)

$$\mathrm{S} =21.91 \sin \mathsf{\theta } -0.96\hspace{1em}\mathsf{\theta }\ge 10°$$

(9)

where: S is the slope factor and θ is the slope (°).

$$L_i=\frac{A_{out}{}^{m+1}-A_{in}{}^{m+1}}{{\left(\mathsf{∆}x\right)}^2·L_m^m·{22.3}^m}$$

(10)

where: L_i is the slope length factor of the i raster; A_in is the confluence area of the grid entrance, m²; A_out is the confluence area of the grid outlet; Δx is the grid resolution, m; L is the non-cumulative slope length related to the flow direction of the grid inlet and outlet; m is the slope length index, according to the research results of Liu et al. [41]. These take the values as follows:

$$\mathrm{m}=0.2\hspace{1em}\mathsf{\theta }<0.5°$$

(11)

$$\mathrm{m}=0.3\hspace{1em}0.5°\le \mathsf{\theta }<1.5°$$

(12)

$$\mathrm{m}=0.4\hspace{1em}1.5°\le \mathsf{\theta }<3°$$

(13)

$$\mathrm{m}=0.5\hspace{1em}\mathsf{\theta }\ge 3°$$

(14)

2.4.4. C, Crop Management Factor

The C factor reflects the effect of vegetation types and cover on soil erosion. Vegetation has the function of intercepting rainfall, slowing down surface runoff, maintaining soil and water, and inhibiting soil erosion [39]. Different vegetation types and cover show different effects in inhibiting soil erosion. The C factor values were assigned by combining the land use/land cover classification map of the study area and the results of previous studies on the C factor in the Loess Plateau and a drawn C factor map [42,43].

2.4.5. P, Conservation Supporting Practice Factor

The P factor expresses the impact of support practice in the reduction of soil erosion, which is defined as the ratio of soil erosion by a support practice to the corresponding loss with upslope and downslope tillage [39]. The P factor ranges from 0 to 1 (non-dimensional). The lower the value, the better the conservation practice. The P factor values were assigned by the results of previous studies on the P factor in the Loess Plateau and a drawn P factor map [42,43,44,45].

2.5. Soil Erosion Analysis

We classified the annual average soil erosion into six classes (slight, light, moderate, severe, very severe, and extremely severe) in the study area (

Table 2. The soil erosion classification index based on the standards for the classification and gradation of soil erosion SL190-2007 issued by the Ministry of Water Resources of China.

Code	Erosion Intensity Ton km−2 year−1	Erosion Intensity Class
1	<1000	slight
2	1000–2500	light
3	2500–5000	moderate
4	5000–8000	severe
5	8000–15,000	very severe
6	>15,000	extremely severe

), based on the standards for the classification and gradation of soil erosion SL190-2007 issued by the Ministry of Water Resources of China. In order to visualize the temporal and spatial changes in soil erosion, we referenced the Sankey diagrams to illustrate the transfer situation of soil erosion classes at different time and used the R software package called networkD3 to map this. Similarly, in the Sankey diagrams, the width of the lines represents the transfer areas of the soil erosion classes, and the height of each component in the stacked bars is proportional to the amount of area occupied by the soil erosion class in the study area.

2.6. Regression Models

2.6.1. Ordinary Least Squares (OLS)

OLS is mainly used for parameter estimation of linear regression, which obtains the minimum value of the sum of the squares of the difference between the actual value and the model estimate and uses the minimum value as the parameter estimation value [46]. In this study, we used OLS as the reference of the spatial regression models. The formula is as follows:

$$Y_i=\alpha +{\beta }_1{X}_1+\cdots +{\beta }_n{X}_n+\epsilon$$

(15)

where Y_i is the dependent variable; X₁ to X_n are independent variables; α is the intercept; β₁ to β_n are estimated regression coefficients; ε is the error term.

2.6.2. Spatial Regression Model

Geographically Weighted Regression Model (GWR)

The geographically weighted regression model (GWR) is a local spatial technique that can be used to examine the spatial variabilities of regression parameters and model performance [46]. Compared with the general liner models, the GWR model carries out separate regressions at each location (observation), considering only other observations within a specific distance to that location. The formula is as follows:

$$Y_i={\beta }_0\left(u_i,v_i\right)+{\displaystyle \sum }_{k=1}^m{\beta }_k\left(u_i,v_i\right)X_{ik}+{\epsilon }_i$$

(16)

where Y_i is the dependent variable; X₁ to X_n are independent variables; β₀ (μ_i, ν_i) is the intercept; μ_i and ν_i represent the coordinates of the sample i in space; β_k are estimated regression coefficients; ε is the error term.

$$\beta \left(u_i,v_i\right)={\left(A^T W \left(u_i,v_i\right) A\right)}^{-1} A{ }^{T}W \left(u_i,v_i\right) C$$

(17)

where β (μ_i, ν_i) is the regression parameter vector in location (μ_i, ν_i); X is a sampling matrix of the independent variable; C is a sampling vector of the dependent variable; W (μ_i, ν_i) is a diagonal matrix, whose diagonal elements represent the geographical weightings of each observation around point i.

$$W_i=\left[\begin{array}{ccc}\begin{array}{cc}{w}_{i1}& 0\end{array}& \cdots & 0\\ \begin{array}{c}\begin{array}{cc}0& w_{i2}\end{array}\\ \begin{array}{cc}⋮& ⋮\end{array}\end{array}& \ddots & \begin{array}{c}0\\ ⋮\end{array}\\ \begin{array}{cc}0& 0\end{array}& \cdots & w_{in}\end{array}\right]$$

(18)

where win is the weight assigned to the observation at location i.

Spatial Lag Model (SLM)

The spatial lag model (SLM) can solve the influence of adjacent independent variables on another independent variable, which was afresh regarded as an explanatory variable to be used for regression [47]. The formula is as follows:

$$Y_{it}=\delta {\displaystyle \sum }_{j=1}^N{W}_{ij}{y}_{jt}+\beta X_{it}+{\mu }_i+{\epsilon }_{it}$$

(19)

where t is an index for the time dimension (years in this study), with t = 1, 2, 3, …, T; i is individual in cross section, with i = 1, 2, 3, …, N; δ is the spatial autoregressive coefficient; W_ij is the row standardized spatial weight matrix; Y_it is the dependent variable at i and t; X_it represents the vector of the independent variables at i and t; β represents the coefficient vector; μ_i denotes a spatial-specific effect for i; ε_it is an error term for i and t.

Spatial Error Model (SEM)

The spatial error model assumes that the error terms have an impact on spatial dependence [48]. The formula is as follows:

$$Y_{it}=\beta X_{it}+{\mu }_i+{\phi }_{it}$$

(20)

$${\phi }_{it}=\rho {\displaystyle \sum }_{j=1}^N{W}_{ij}{\phi }_{jt}+{\epsilon }_{it}$$

(21)

where φ_it is the spatially autocorrelated error term; ρ is the spatial autocorrelation coefficient; ${\displaystyle \sum }_{j=1}^N{W}_{ij}{\phi }_{jt}$ denotes the interaction effects among the disturbance terms of the different units; λ refers to the spatial autocorrelation coefficient.

In this study, we used soil erosion change as the dependent variable and land use/cover area change as the independent variable to undertake regression analysis. The OLS model and the GWR model were run in ArcMap10.7 software. The SLM model and the SEM model were run in GeoDa 1.18 software (http://geodacenter.github.io/formats.html, accessed on 6 August 2020).

3. Results

3.1. LULCC Classification and Accuracy Evaluation

The LULCC classification maps of Changwu Country were generated for each of the following years: 1987, 1997, 2007, and 2017 (Figure 2). The overall classification accuracy in 1987, 1997, 2007, and 2017 was 89.1%, 90.2%, 92.4%, and 93.5%, respectively, and the respective Kappa coefficients were 0.81, 0.87, 0.90, and 0.92, respectively, which met the needs for the present study.

3.2. LULCC Analysis

The rate and trend of changes markedly varied between land use/cover in the intervals of the study period. We divided our data into the three research periods (1987–1996, 1997–2006, and 2007–2017), and the differences in LULCC were significant during each study period. Cropland expanded rapidly in the period of 1987–1997 at a rate of 3.48% per annum. However, from 1997 to 2007, the cropland area had a substantial shrinkage at a rate of 4.97% per annum. By 2007, this was only 146.11 km², accounting for 50.31% of the cropland area in 1997. The cropland area successively decreased in the period of 2007–2017, but the rate of decline slowed down, and the area decreased by 25%. The orchard area decreased in the first period of 1987–1997, and successively increased in the other two periods (1997–2007 and 2007–2017). A comparison of 1987 and 1997 revealed a reduction in orchard area percentage from 20.95% to 12.73%. However, from 1997 to 2007, the orchard area rapidly increased at a rate of 6.66% per year, with an increase of 48.15 km². From 2007 to 2017, the orchard area continued to increase at a rate of 3.34% per year, with an increase of 40.19 km².

The construction land area successively increased from 1987 to 2017, although the changes varied in different periods in the whole area. The construction land area increased by 37.85% and 23.58% from 1987 to 1997 and from 1997 to 2007, respectively. However, the rate of construction land increase quickly changed with the rapid economic development and increased by 313.50% from 2007 to 2017. This covered 34.96 km² in 2017, which was the largest in the last 30 years. The water body area drastically changed, with an increase of 79.53% from 1987 to 1997, a decrease of 37.67% from 1997 to 2007, and a slight increase of 8.51% from 2007 to 2017.

In the 1950s, the Changwu County Government began to carry out large-scale locusts by organizing people, mainly to obtain timber and firewood. In the 1990s, due to the need for economic development, a large number of forest trees were cut for the construction of houses and mine pillars. This led to the substantial shrinkage of forest land, reduced by 52.34% from 1987 to 1997. As the government’s awareness of forest land protection gradually increased, the government began to implement the Grain for Green Project in 1999, and the forest land was quickly restored. In 2007, forest land rapidly increased by 310.49%, when compared with 1997. In 2017, there was a continuous increase in forest land, with an increase of 37.55%, when compared with 2007.

However, the planted land and forest land change were inconsistent. The planted land decreased by 7.01% from 1987 to 1997, increased by 19.68% from 1997 to 2007, and decreased by 33.29% from 2007 to 2017.

The cropland, orchard, planted land and forested land were converted into large areas with other land types at each time period, and there was less conversion between construction land, water bodies, and other land types (Figure 3). From 1987 to 1997, the largest increase in area of cropland was mainly due to the conversion of planted land and orchards, which transformed 68.99 km² and 64.93 km², respectively. The forest land was mainly converted into planted land and cropland, which were converted to 26.92 km² and 7.91 km², respectively. The large-scale conversion of cropland and planted land occurred between 1997 and 2007. The conversions of cropland to planted land and forest land were 80.33 km² and 19.72 km², respectively. The area converted from cropland to orchard was 81.46 km², which was the largest conversion between all land types in this period.

From 2007 to 2017, the area converted from planted land to forested land was the largest among all land types at 54.46 km², followed by the conversion of cropland to orchards, which was 52.94 km². The rapid expansion of construction land mainly came from the conversion of cropland, orchards and planted land to cope with the development process of urbanization.

3.3. Soil Erosion Estimation and Analysis

We used the RUSLE to estimate the soil erosion in Changwu County. Figure 4 shows the spatial distribution of soil erosion intensity in 1987, 1997, 2007, and 2017. The annual average soil erosion for Changwu Country was estimated at 977.84 ton km⁻² year⁻¹, 1305.17 ton km⁻² year⁻¹, 1310.60 ton km⁻² year⁻¹, and 1891.46 ton km⁻² year⁻¹ in 1987, 1997, 2007, and 2017, respectively. The rates of annual average soil erosion increase were significantly different between the periods of 1987–1997, 1997–2007, and 2007–2017, with increases at the rate of 3.35%, 0.04%, and 4.43% per year in these three time periods, although the annual average soil erosion continued to increase over time.

Figure 5 shows the spatial distribution of six classes of annual average soil erosion in 1987, 1997, 2007, and 2017. The areas of slight erosion and light erosion continued to decrease, but the areas of moderate erosion and severe erosion continued to increase during these three periods (1987–1997, 1997–2007, and 2007–2017) (Figure 5). The areas of very severe erosion and extremely severe erosion decreased during the period of 1997–2007 but increased during the period of 1987–1997 and 2007–2017. The area of slight erosion was the largest, accounting for 70.88%, 66.30%, 65.79%, and 61.31% of the study area in every year. As the soil erosion intensity increased, the corresponding area gradually decreased.

The transformation amount of soil-erosion class area varied at different time periods. The amount of transformation mainly occurred when slight erosion was converted to light erosion, light erosion was converted to moderate erosion, and moderate erosion was converted to light and severe erosion (Figure 6). The area converted from slight erosion to light erosion was the largest at 33.30 km², followed by light erosion converted to moderate at 27.11 km² and then moderate erosion converted to light erosion at 13.73 km², during the period of 1987–1997. However, during the period of 1997–2007, the area converted from light to moderate erosion was the largest at 20.32 km², indicating that soil erosion intensity was increasing, followed by slight erosion converted to light erosion at 19.78 km² and moderate erosion converted to slight erosion at 18.91 km². It noteworthy that not only the area converted from light to moderate erosion was the largest at 32.83 km² during the period of 2007–2017, but moderate erosion converted to severe erosion at 21.21 km² also first appeared as the third in the three time periods (1987–1997, 1997–2007 and 2007–2017). This result further revealed that the soil erosion intensity was increasing. The transformation amount of slight erosion, light erosion to severe erosion, very severe erosion, and extremely severe erosion was very small, indicating that slight erosion and light erosion were relatively stable without drastic changes.

3.4. LULCC Effect on Soil Erosion Change and Model Comparison

The OLS and Spatial regression models (GWR, SLM and SEM) were applied to fit the relationship between LULCC and soil erosion change and their accuracies R² were shown in

Table 3. Accuracy R² value (p < 0.01) of different models in fitting the effect of LULCC on soil erosion change in different periods.

Period	OLS	GWR	SLM	SEM
1987–1997	0.025	0.203	0.653	0.659
1997–2007	0.047	0.211	0.644	0.639
2007–2017	0.004	0.300	0.627	0.639

. The results showed that the SEM and SLM have higher accuracy than the GWR and OLS in fitting the effect of LULCC and soil erosion change in different periods. In 1987–1997, SEM has the highest accuracy (R² = 0.659), followed by SLM (R² = 0.653), GMR has lower accuracy (R² = 0.203), and OLS has the lowest accuracy (R² = 0.025). However, SLM has the highest accuracy (R² = 0.644) in 1997–2007, followed by SEM (R² = 0.639), GMR has lower accuracy (R² = 0.211), and OLS has the lowest accuracy (R² = 0.047). In 2007–2017, SEM and SLM have higher accuracies (R²) at 0.639 and 0.627, respectively, and GMR and OLS have lower accuracies (R²) at 0.300 and 0.004, respectively. In summary, the SEM and SLM have better explanations for the effect of LULCC on soil erosion change.

4. Discussion

LULCC reflects the results of human–environment interactions [49,50]. The Loess Plateau has been disturbed by human activities for a long time, and human activities are the main factors that affect LULCC in the Loess Plateau [51,52]. Changwu County is an area dominated by agricultural production. Population expansion has led to long-term pressure on resources. The results indicate that the cropland area expanded significantly between 1987 and 1997 (Section 3.2). This was due to the population growth, backward economic development during this period, and per capita income of farmers which depended on the expansion of cropland to improve. Hence, large areas of planted land were converted into cropland. The orchard mainly comprised of persimmon orchard and apple orchard during this period. However, a large area of persimmon and apple trees fell and these areas were turned into cropland, due to economic benefits. However, the agricultural structure has changed as the orchard area continues to increase, while cropland has continued to decrease from 1997 to 2017, which plays an important role in alleviating pressure on resources. This change was correlated with the national policy on ecological environment protection, but the huge benefits brought by orchards were the direct driving force for the conversion of cropland into orchards. Wang [53] investigated the rural transformation from the perspective of regime shifts of socio-ecological systems in Changwu County and reported that the apple’s income (41,370 yuan/ha) was 138 times of that of grain in 2007. This result supports the present research results very well.

The results of the spatial and temporal changes in forest land in the study area reflect the changes in people’s attitudes towards the environment. Protecting the environment has attracted the government’s attention. The area of forested land has continuously increased since the implementation of Grain for Green Project in 1999. The implementation of the Grain for Green Project resulted in large areas of cropland being used for afforestation. Hence, cropland was mainly transformed into planted land and forested land between 1997 and 2007. During this period, the government’s awareness of forest land protection gradually increased, and the forest was quickly restored as the forest tree cover quickly grew. Some unclosed forests grew into forests. Therefore, forest land mainly comes from the conversion of planted land and cropland. In 2017, the forested area reached 123.93 km², accounting for 21.8% of the research area. This result indicates that the forest land was continuously recovered with the government’s emphasis on the ecological environment, and that the agricultural structure is continuously changing. The conversion of cropland to planted land was 18.78 km² between 2007 and 2017, which indicates that the conversion of cropland to forests is continuous.

The complex structure of forests could intercept rainfall and reduce soil erosion, thereby improving soil protection [54,55]. Hence, the increase in forest area not only benefits water and soil conservation, regulates the climate, and purifies the environment, but also brings business opportunities to the local tourism industry. Forest-based well-being tourism is becoming more and more popular due to the comfortable environment provided by the forest, which derived the development of the tertiary industry from 67.50 million yuan in 1997 to 1843 million yuan during the period. Furthermore, R. pseudoacacia forests are an important source of honey plants, and a large area of R. pseudoacacia forests have attracted bee farmers. The ecological service function provided by forests for humans guarantees the continuous increase in forest land. The continuous expansion of construction land area reflects the improvement of people’s quality of life. That is, from a GDP of 544.3 million yuan in 1987 to a GDP of 9.596 billion yuan in 2017. Especially from 2007 to 2017, Changwu Country witnessed rapid economic development and urbanization, which all determined the expansion of construction land. Xiao et al. [56] investigated the urbanization and land use change in Shijiazhuang in China, and reported that construction land, as a ‘fast-expansion stage’, was at the period of fast economic development. This coincides with the present research results.

Changwu County is an area with severe soil erosion on the Loess Plateau. The average annual soil erosion was 6500 ton/km² from 1949 to 1979 in Changwu County, according to the statistical data of soil and water conservation in Shaanxi province (derived from http://loess.geodata.cn, accessed on 6 August 2020).

However, it was found that the erosion in the study area was mainly slight and light erosion in Changwu County in the past 30 years. Severe erosion, very severe erosion, and extremely severe erosion account for only a small amount. It can be observed from the spatial and temporal variation of soil erosion that the transformation mainly occurred in low-level soil erosion. Furthermore, the amount of high-level soil erosion transformation was less. This shows that the earlier soil and water conservation projects played an important role in the water and soil management in the study area.

In 2020, Wen [21] studied the cost effectiveness of soil erosion control practices in the Chinese Loess Plateau and found that engineering erosion control techniques implemented in the Loess Plateau are effective in controlling soil erosion, which is consistent with the present research results.

As early as the 1950s, the state recognized the importance of protecting soil and water loss but did not put forward scientific treatment methods. The treatment of soil erosion was still in the exploratory stage, and farmers were called upon to build protection projects. It was not until 1981 that a scientific treatment method was developed to comprehensively treat mountains, rivers, cropland, forest, and roads, according to the scope of the drainage basin, and effectively control the local soil and water loss.

When we analyzed the temporal and spatial changes of soil erosion, it was found that the soil erosion changes in the study area had a kind of “inertia”. This “inertia” means that each erosion level prefers to convert to a neighboring level, rather than a large number of cross-level transformations. Furthermore, this “inertia” allows for the optimistic estimation that the erosion in the study area is at a relatively safe stage, and that this would not rapidly deteriorate.

However, what we must pay attention to is that the annual average soil erosion, moderate erosion, and severe erosion in the study area all increased with time. Furthermore, the very severe erosion and extreme erosion area was the largest in 2017, when compared with the other time points (1987, 1997, and 2007). Hence, there is a need to further strengthen the soil and water conservation work in Changwu County.

From the results of the LULCC, it can be observed that from 1987 to 1997, the forest land decreased, and the people’s logging, illegal logging and cropland expansion caused great damage to the local vegetation. This is also the direct cause of the increase in soil erosion intensity from 1987 to 1997. The vegetation began to recover after the state implemented the Grain for Green Project in 1999. Vegetation restoration can intercept rainfall, reduce surface runoff, and alleviate soil erosion [57,58,59]. The influence of vegetation on soil erosion can be reflected by the value of the C factor. The C factor value increased from 1987 to 1997, but the C factor value continued to decrease from 1997 to 2017, indicating that the increase of vegetation area may reduce soil erosion. However, the R factor value was determined by rainfall. From 2007 to 2017, the rapid increase in the R factor value was due to the increase in rainfall, which resulted in increased soil erosion intensity. This was the direct basis for the different effects of LULCC and rainfall on soil erosion in different periods.

This study used spatial regression model comparison for the first time to determine the impact of LULCC on soil erosion change and has significance for future studies. The results show that the relationship between soil erosion and LULCC obtained by SLM and SEM is more reliable than that derived by GWR and OLS. As spatial panel data models, SLM and SEM may consider the spatial and temporal impact on dependent variables to obtain accurate estimates. Fang [60] found that panel data models may be applied to the impact of land use changes on catchment soil erosion and sediment yield in northeastern China. This result supports our research, and our research further expands the application of panel data in LULCC and soil erosion change. This study used the selection of spatial regression models to obtain the effect of LULCC on soil erosion change, which can provide a reference for the relationship between LULCC and soil erosion change in other regions.

5. Conclusions

The LULCC in Changwu County has been significant over the last 30 years. The implementation of the Grain for Green Project, the change in agricultural structure, and good income from orchards, have controlled the expansion of cropland which was conducive to the harmony between humans and nature. The decrease to increase in forest land was the direct result of the implementation of the Grain for Green Project, reflecting the importance of people’s attachment to environmental issues. Construction land has continued to increase in the past 30 years, which is closely correlated with the local economic development and urbanization process. The soil erosion in Changwu County in the past 30 years has been dominated by slight erosion and light erosion, and other erosion levels have been less prevalent. The early soil and water conservation projects played an important role. The changes in soil erosion in the study area had both positive (the soil erosion changes showed a kind of “inertia”) and negative (the intensification trend of soil erosion intensity) aspects. The SEM and SLM have higher accuracy in explaining the effect of LULCC and soil erosion change in different periods. The explanatory power of LULCC on soil erosion exceeded 0.62 in different periods.