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

: 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 signiﬁcantly changed in the last 30 years. From 1987 to 1997, cropland expansion came mainly from planted land and orchards, which transformed 68.99 km 2 and 64.93 km 2 , 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 − 2 year − 1 , 1305.17 ton km − 2 year − 1 , 1310.60 ton km − 2 year − 1 , and 1891.46 ton km − 2 year − 1 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 ﬁtting the effect of LULCC and soil erosion change, where the accuracy exceeded 0.62 in different periods.


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 nately had a gully development, forming crisscross gullies. The landforms were ma "Yuan" and slope land, followed by "Chuan". The landform features were characteri by high "Yuan", deep gully, and steep slope. This area was drained by three main r systems, namely, the Jinghe River, Heihe River and Nanhe River, which are tributarie the Yellow River. The dominant soil type was dark loessial soil, which is at "Yuan", loessial soil, which is at the valley slope and gully channel. The altitude ranged from m to 1274 m above sea level [23,28].
Changwu County belongs to a warm, temperate, semi-humid continental mons climate zone, with an average annual temperature of 9.1 °C. Furthermore, the frostperiod is 171 days, and the annual precipitation is 587.8 mm. The annual rainfall conc tration period is 7, 8 and 9 months, and rainfall can reach up to 321.4 mm, accounting 54.9% of the total annual precipitation. The study area was mainly planted forest, afforestation began in the 1950s. The main tree species are Robinia pseudoacacia, Platycla orientalis, Pinus tabuliformis, Armeniaca sibirica, and Populus spp. Among these, R. pseu acacia accounts for more than 90% of the forested area [22,23,27].
The research area has few land resources and has experienced agricultural recla tion for a long time. Agriculture occupies a major position in the county's national ec omy. The cropland and orchard are mainly concentrated on the plateau and gully bott and the main plants were wheat, apple, and persimmon.

Data Sources
The land use/cover maps of Changwu County on 28 August 1987, 23 August 1997 September 2007, and 27 June 2017 were derived from Landsat TM/ETM+/OLI image with 30 m resolution (derived from the USGS website, http://www.usgs.gov, accessed 16 July 2019), considering the vegetation growing season and image quality. The dig elevation model (DEM) was sourced from the geospatial data clo (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 image processing software to eliminate the scattering and absorption effects from the mosphere and to obtain authentic the surface reflectance character. Finally, the ima 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.

Data Sources
The land use/cover maps of Changwu County on 28 August 1987, 23 August 1997, 20 September 2007 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), 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]: 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.

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: 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.
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]: 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].

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: where: A is the average annual soil erosion rate (t km −2 year −1 ); R is the rainfall erosivity factor (MJ mm hm −2 h −1 year −1 ); K is the soil erodibility factor (thm 2 h hm −2 MJ −1 mm −1 ); LS is the topographic factor; C is the crop management factor; P is the conservation supporting practice factor.

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].
where: P i is the monthly amount of rainfall and P is the annual rainfall.

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).

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].
where: S is the slope factor and θ is the slope ( • ).
where: L i is the slope length factor of the i raster; A in is the confluence area of the grid entrance, m 2 ; 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:

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].

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].

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), 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. 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: where Y i is the dependent variable; X 1 to X n are independent variables; α is the intercept; β 1 to β n are estimated regression coefficients; ε is the error term.

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: where Y i is the dependent variable; X 1 to X n are independent variables; β 0 (µ 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.
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.
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: 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: where ϕ it is the spatially autocorrelated error term; ρ is the spatial autocorrelation coefficient; N ∑ j=1 W ij ϕ 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).

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.

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. 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 2 and 64.93 km 2 , respectively. The forest land was mainly converted into planted land and cropland, which were converted to 26.92 km 2 and 7.91 km 2 , 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 2 and 19.72 km 2 , respectively. The area converted from cropland to orchard was 81.46 km 2 , which was the largest conversion between all land types in this period.

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 −2 year −1 , 1305.17 From 2007 to 2017, the area converted from planted land to forested land was the largest among all land types at 54.46 km 2 , followed by the conversion of cropland to orchards, which was 52.94 km 2 . 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.

Soil Erosion Estimation and Analysis
We used the RUSLE to estimate the soil erosion in Changwu County. Figure 4 shows     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 2 , followed by light erosion converted to moderate 2 2 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 2 , followed by light erosion converted to moderate at 27.11 km 2 and then moderate erosion converted to light erosion at 13.73 km 2 , 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 2 , indicating that soil erosion intensity was increasing, followed by slight erosion converted to light erosion at 19.78 km 2 and moderate erosion converted to slight erosion at 18.91 km 2 . It noteworthy that not only the area converted from light to moderate erosion was the largest at 32.83 km 2 during the period of 2007-2017, but moderate erosion converted to severe erosion at 21.21 km 2 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. 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.

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 2 were shown in Table 3. 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 2 = 0.659), followed by SLM (R 2 = 0.653), GMR has lower accuracy (R 2 = 0.203), and OLS has the lowest accuracy (R 2 = 0.025). However, SLM has the highest accuracy (R 2 = 0.644) in 1997-2007, followed by SEM (R 2 = 0.639), GMR has lower accuracy (R 2 = 0.211), and OLS has the lowest accuracy (R 2 = 0.047). In 2007-2017, SEM and SLM have higher accuracies (R 2 ) at 0.639 and 0.627, respectively, and GMR and OLS have lower accuracies (R 2 ) 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.

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

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 2 were shown in Table 3. 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 2 = 0.659), followed by SLM (R 2 = 0.653), GMR has lower accuracy (R 2 = 0.203), and OLS has the lowest accuracy (R 2 = 0.025). However, SLM has the highest accuracy (R 2 = 0.644) in 1997-2007, followed by SEM (R 2 = 0.639), GMR has lower accuracy (R 2 = 0.211), and OLS has the lowest accuracy (R 2 = 0.047). In 2007-2017, SEM and SLM have higher accuracies (R 2 ) at 0.639 and 0.627, respectively, and GMR and OLS have lower accuracies (R 2 ) 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.

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 2 , 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 2 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 'fastexpansion 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 2 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.

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.