Spatial Pattern of Soil Erosion Drivers and the Contribution Rate of Human Activities on the Loess Plateau from 2000 to 2015: A Boundary Line from Northeast to Southwest

: The Loess Plateau is one of the most fragile areas in the world, where the problem of soil erosion is particularly prominent. The spatial and temporal variation characteristics and mechanisms of soil erosion in this region have always been hot topics for researchers. In this study, Revised Universal Soil Loss Equation (RUSLE) is used to estimate the soil erosion modulus of the Loess Plateau from 2000 to 2015, the dynamic characteristics of its temporal and spatial variations and driving mechanisms are determined, and meteorological data are combined with remote sensing data to quantitatively calculate the contribution rate of human activities. The results show that from 2000 to 2015, the soil erosion modulus of the Loess Plateau had a downward trend as a whole, with a rate of − 0.6408 t / ha / a, but the downward trend gradually slowed down. Precipitation mainly resulted in changes in the soil erosion modulus in the northwestern part of the Loess Plateau, where a signiﬁcant positive correlation was seen. Meanwhile, the Vegetation Fractional Coverage (VFC) mainly a ﬀ ected the southeastern part, where a signiﬁcant negative correlation was measured. The human-activity contribution rate was − 1.0774 on the Loess Plateau, which means human activities e ﬀ ectively reduced the soil erosion modulus while climate change promoted soil erosion combined with the result of the analysis of variance (ANOVA). “Hilly and gully regions” and “Gully region of Loess Plateau” as the main implementation areas of ecological projects, human activities had contribution rate of 0.5513 and 0.7805 toward the declining of soil erosion, respectively. Interestingly, the spatial di ﬀ erentiation characteristic of the soil erosion driving mechanisms and human contribution rates on the Loess Plateau showed the same boundary line from northeast to southwest, which was well explained by the 400-mm isohyetal line and Hu’s Line. This boundary can guide the geographical layout of the ecological management projects and urban development spaces on the Loess Plateau.


Introduction
Soil erosion is the displacement of soil from its formation place caused by external forces such as wind, raindrops, and water scouring, which occur on a point scale [1]. Soil erosion has always been one of the most important issues in ecology and agriculture, especially on the Loess Plateau of China. The "point-edge contact scaffold porous structure" of loess on the Loess Plateau is the primary cause of the loose soil, the development of vertical joints and the strong water permeability therein [2]. Located in the middle reaches of the Yellow River and the upper reaches of the Haihe River, the Loess Plateau is one of the four main plateaus in China. The Loess Plateau is more than 1,000 kilometers long from east to west and 750 kilometers wide from north to south; spanning seven provinces of China, including the vast areas west of Taihang Mountain, east of Riyue Mountain, north of Qinling Mountains and south of the Great Wall, with a total area of 624,000 km 2 . The geographical position is between 100°54′-114°33′ E and 33°43′-41°16′ N.
The Loess Plateau is located in a transitional zone consisting of a continental monsoon climate of temperate, arid and semi-arid zones, with cold winters, warm and humid summers, and synchronous rain and heat. The annual average temperature is between 3.6-14.3°C, which is low in the northwest and high in the southeast. The altitude is 10-1400 m, with an average altitude of about 240 m. The average annual precipitation is 150-750 mm, which is mainly concentrated in July, August and September, and accounts for 55-78% of the total annual precipitation. Due to the particularities of its geology, geomorphology, hydrology, soil, climate and vegetation, the Loess Plateau is an area that is very susceptible to serious soil erosion and is one of the most fragile ecological environments in the world as seen in (Figure 1).

Data and Materials
In this study, remote sensing data and meteorological data were combined, and two sets of soil erosion modulus simulation data were created to reveal the driving mechanism and contribution rate of human activities of soil erosion on the Loess Plateau. The specific technical process is shown in Figure 2.

Data and Materials
In this study, remote sensing data and meteorological data were combined, and two sets of soil erosion modulus simulation data were created to reveal the driving mechanism and contribution rate of human activities of soil erosion on the Loess Plateau. The specific technical process is shown in Figure 2. The meteorological data used in our study were obtained from the China Meteorological Science Data Sharing Service Network [22]. Daily data sets of meteorological data from 147 national meteorological stations on the Loess Plateau and its surrounding 150 km range were collected ( Figure  3), including the daily maximum temperature, minimum temperature, average temperature, precipitation, average relative humidity. The ANUSPLINE software was used to interpolate the site data to obtain the required spatial data set.  The meteorological data used in our study were obtained from the China Meteorological Science Data Sharing Service Network [22]. Daily data sets of meteorological data from 147 national meteorological stations on the Loess Plateau and its surrounding 150 km range were collected (Figure 3), including the daily maximum temperature, minimum temperature, average temperature, precipitation, average relative humidity. The ANUSPLINE software was used to interpolate the site data to obtain the required spatial data set. The meteorological data used in our study were obtained from the China Meteorological Science Data Sharing Service Network [22]. Daily data sets of meteorological data from 147 national meteorological stations on the Loess Plateau and its surrounding 150 km range were collected ( Figure  3), including the daily maximum temperature, minimum temperature, average temperature, precipitation, average relative humidity. The ANUSPLINE software was used to interpolate the site data to obtain the required spatial data set.  The soil-texture data were derived from the Harmonized World Soil Database version 1.2 (HWSD V1.2), which contains soil, sand, silt, clay and organic carbon content. The soil-types data were derived from China's 1:1 million soil databases built by the Institute of Soil Science, Chinese Academy of Sciences. The basic mapping units including 12 orders, 61 great groups and 235 subgroups. Two sets of soil data were superimposed using the ArcGIS software to obtain vector data of soil type and texture.
The topographic data were derived from the advanced space borne thermal emission and reflection radiometer-global digital elevation model (ASTER-GDEM), where the slope and length factor data were further calculated by us.
The Normalized Difference Vegetation Index (NDVI) data were derived from MODIS-NDVI, once every 16 days, and the Vegetation Fractional Coverage (VFC) data were further calculated using the ENVI software.

Soil Erosion Modulus Simulation
In 1965, Wischmeier and Smith proposed the famous universal soil loss equation (USLE) [23], and then American scientists established the revised universal soil loss equation (RUSLE) in 1997 [24]. With the development of the geographic information system (GIS) and remote-sensing (RS) technology, RUSLE has become the most widely used empirical model for the quantitative assessment of soil erosion at watershed and regional scales [25]. Its basic form is: where A is the average soil erosion modulus [t/(ha·a)], R is the rainfall erosivity factor [MJ·mm/(ha·h·a)], and K is the soil erodibility factor [t·ha·h/(ha·MJ·mm)], LS is the slope length and slope factor, C is the vegetation cover and management factor, and P is the soil and water conservation measure factor. The latter four factors are all dimensionless.
In this study, we generated 1 × 1 km 2 raster data of these factors in ArcGIS, and then obtained soil erosion modulus raster data of 1 × 1 km 2 resolution.
(1) Rainfall Erosivity Factor (R) The national daily rainfall fitting model proposed by Zhang, Xie & Liu [26] was used to estimate the half-month rainfall erosivity factor in this study. The formula is as follows: where R i is the rainfall erosivity value of the i th 16-day period of the year [MJ·mm/(ha·h·a)], D j is the erosive daily rainfall on the j th day of the 16-day period, and k = 16. The daily rainfall is required to be greater than or equal to 12 mm, otherwise it is calculated as 0 mm. The threshold value of 12 mm is consistent with the standard of erosive rainfall in China. α and β are undetermined parameters of the model, which can be calculated as: where P d12 represents the daily average rainfall of more than 12 mm (including 12 mm), and P y12 represents the annual average rainfall of more than 12 mm (including 12 mm).
(2) Soil Erodibility Factor (K) In 1990, Williams expressed erosion-productivity impact calculator (EPIC) [27]. In this study, we used the calculation method outlined in the "Ecological protection red line delineation guide" issued by the Ministry of Ecology and Environment of the People's Republic of China and the National Development and Reform Commission [28], which revised the formula for calculating soil erodibility in EPIC as: where K is the revised soil erodibility factor, K EPIC is the soil erodibility factor before revision, SAN, SIL and CLA are the content of soil sand, silt and clay (%), respectively, and C is the soil organic carbon content (%).
(3) Slope Length (L) and Slope Factor (S) Since the slope and slope length factors are closely related to each other, they are usually considered together. In this study we combined the method proposed by McCool et al. [29] and Liu et al. [30]. The calculation formulas are as follows: where L is the slope length factor (dimensionless), γ is the slope length (m), m is a dimensionless constant that depends on the degree of the slope, θ is the slope (%), and S is the slope factor (dimensionless).

(4) Vegetation Cover and Management Factor (C)
The algorithm established by Cai et al. [31], which is based on a large number of experiments, is widely respected because of its simplicity and accuracy, and it was used to simulate the C factor by VFC in this study: where C is the vegetation cover and management factor (dimensionless). VFC is calculated using the NDVI as follows [32]: where NDVI soil is the NDVI value of bare soil pixels, NDVI veg is the NDVI value of pure vegetation pixels, NDVI max and NDVI min are the maximum and minimum NDVI values in the region, respectively, and VFC max and VFC min are the maximum and minimum VFC values in the region, respectively. When VFC max = 100%, VFC min = 0%, and Equation (11) becomes: In the actual calculation process, NDVI max and NDVI min are the raster values for a cumulative rate of 95% and 5% from small to large in the area, respectively.

Statistical Analysis
In this study, the least squares method was used to linearly fit the multi-year trends of the soil erosion modulus, precipitation and vegetation coverage [33]. The calculation is as follows: whereb is the slope of the fitted line,â is the intercept of the fitted line, y consists of variables including the soil erosion modulus, precipitation and VFC in this study, and x is the corresponding year. y and x are the n-year average of y and x, y i and x i are the soil erosion modulus, precipitation or VFC of the i th year, and n is the number of years, i.e. the number of samples. Pearson's correlation coefficient (r) was used to describe any simple linear correlation that may be present between two variables (e.g., the soil erosion modulus, precipitation and vegetation coverage and the year), and the consistency of the temporal and spatial patterns measured from the time-series data (soil erosion modulus and precipitation, soil erosion modulus and vegetation coverage), which is calculated as follows [34]: The significance level is expressed by the Student's t-test, where the t-statistic is calculated as follows: where r is Pearson's correlation coefficient, X and Y are the two variables described, X and Y are the n-year averaged values of X and Y, respectively. X i and Y i represent the value of X and Y in the i th year, respectively, d f is the degrees of freedom, and n is the number of samples. According to the calculated t-statistic, the corresponding P value was obtained by combining the distribution threshold table of t values to determine the significance level between the considered variables (soil erosion modulus and precipitation, soil erosion modulus and vegetation coverage).

Soil Erosion Modulus Simulation under Meteorological Conditions
Calculation of the soil erosion modulus under different meteorological conditions was performed using RUSLE, which calculates the potential vegetation coverage factor (C) under meteorological conditions instead of the vegetation coverage factor obtained from RS data.
Calculation of the VFC was based on Equations (11)(12)(13)(14), which were employed using the NDVIs. The NDVIs were in turn calculated using the algorithm proposed by Woodward and Smith [35]. The specific formulas are as follows: where NDVI m is the potential NDVI simulated by the meteorological data, LAI is the potential leaf area index, P is the monthly precipitation (mm), D is the air humidity, t is the day length (hr), and G max is the maximum stomatal conductance (mmol·m −2 ·s −1 ), which is calculated using the empirical observations collated by Woodward and Smith [35]. The G max value of a cool temperate steppe is 287 mmol·m −2 ·s −1 , while G max of a cool temperate moist forest is 300 mmol·m −2 ·s −1 .
The daily length, t, was calculated using the formula proposed by the Food and Agriculture Organization of the United Nations (FAO) [36] as: where ω s is sunset hour angle (rad), ϕ is latitude of the meteorological station, δ is the solar decimation (rad), and J is the number of the day in the year between 1 (1 st January) and 365 or 366 (31 st December).

Determination of the Contribution Rate of Human Activities
By calculating the trend of the soil erosion modulus of the Loess Plateau from 2000 to 2015 under simulated and actual meteorological conditions, the method of obtaining contribution rate of human activities is postulated in this study as follows: where C h is the contribution rate of human activities (dimensionless), where the numerical value indicates the effect of human activities on the change of soil erosion rate and a positive or negative value represents an increase or decrease of the erosion modulus,b a is the change slope of soil erosion modulus in the Loess Plateau under actual conditions [t/(ha·a)], andb m is the slope of the soil erosion modulus on the Loess Plateau simulated by the meteorological data [t/(ha·a)].

Spatial Distribution Characteristics
From 2000 to 2015, the annual average erosion modulus of the Loess Plateau was 20.8856 t/ha/a (Table 1), and the yearly values showed similar spatial distribution patterns ( Figure 4). There was a Remote Sens. 2019, 11, 2429 9 of 20 high-value zone that extended from southwest to northeast, which was mainly located in the "hilly and gully region" and "gully region of Loess Plateau" of the comprehensive control regions in the Loess Plateau. Regionally, the "hilly and gully region" had the highest annual average soil erosion modulus, 32.6789 t/ha/y, while the "plateau gully region" of the Loess Plateau has the second highest, 28.5799 t/ha/a. "Sandy areas" constitute the smallest region, with 421.82 t/ha/a.

Spatiotemporal Variation of Soil Erosion
Over the past 16 years, the highest soil erosion modulus was 29.3461 t/ha/a in 2001, and the lowest was 13.7385 t/ha/a in 2015. From 2000 to 2015, the average soil erosion modulus of the whole Loess Plateau showed a downward trend, with a rate of −0.6408 t/ha/a (P < 0.05). However, the downward trend has gradually slowed down ( Figure 5). Regionally, the "hilly and gully region" had the highest annual average soil erosion modulus, 32.6789 t/ha/y, while the "plateau gully region" of the Loess Plateau has the second highest, 28.5799 t/ha/a. "Sandy areas" constitute the smallest region, with 421.82 t/ha/a.

Spatiotemporal Variation of Soil Erosion
Over the past 16 years, the highest soil erosion modulus was 29.3461 t/ha/a in 2001, and the lowest was 13.7385 t/ha/a in 2015. From 2000 to 2015, the average soil erosion modulus of the whole Loess Plateau showed a downward trend, with a rate of −0.6408 t/ha/a (P < 0.05). However, the downward trend has gradually slowed down ( Figure 5). Regionally, the annual soil erosion modulus in the six comprehensive control regions all showed a downward trend. Among them, the most obvious downward trends were in the "hilly and gully region" and the "gully region" of the Loess Plateau, which were −0.0107t/ha/a (P < 0.05) and −0.0099 t/ha/a (P < 0.05), respectively ( Figure 6).  Regionally, the annual soil erosion modulus in the six comprehensive control regions all showed a downward trend. Among them, the most obvious downward trends were in the "hilly and gully region" and the "gully region" of the Loess Plateau, which were −0.0107t/ha/a (P < 0.05) and −0.0099 t/ha/a (P < 0.05), respectively ( Figure 6). Regionally, the annual soil erosion modulus in the six comprehensive control regions all showed a downward trend. Among them, the most obvious downward trends were in the "hilly and gully region" and the "gully region" of the Loess Plateau, which were −0.0107t/ha/a (P < 0.05) and −0.0099 t/ha/a (P < 0.05), respectively ( Figure 6).

Accuracy Verification
Next, we compared the results of other researchers to verify the accuracy of our soil erosion modulus estimates (Table 2). Due to the different algorithms and different model parameters, the calculation period, area and resolution had some differences and deviations with the other considered studies. The results of this study were close to those of Sun et al. [37] and Qin et al. [38], which suggested that the calculations performed in this study were credible.

Response Mechanism of Soil Erosion Change to VFC and Precipitation
VFC and precipitation are the most direct and important factors affecting soil erosion change, and their annual variation ranges are much larger than other factors [37,41,42]. Understanding the driving relationship between them and soil erosion change can help deepen our understanding of the processes and mechanisms of soil erosion change.
Between 2000 and 2015, the VFC of the Loess Plateau increased as a whole, showing obvious regional differences. In the central part of the Loess Plateau, the VFC values of the "hilly and gully region" and the "plateau gully region" showed a significant upward trend. However, the VFC had a declining trend in most of the "plain areas." some regions of the "earth-rocky mountainous," "sandy areas" and "irrigated regions" (Figure 7a). From 2000 to 2015, precipitation in most areas of the Loess Plateau showed an upward trend, with the biggest increases in the "hilly and gully region" and the "Gully region of Loess Plateau". Only the southeastern and western marginal regions of the Loess Plateau showed a downward trend (Figure 7b). Next, we compared the results of other researchers to verify the accuracy of our soil erosion modulus estimates (Table 2). Due to the different algorithms and different model parameters, the calculation period, area and resolution had some differences and deviations with the other considered studies. The results of this study were close to those of Sun et al. [37] and Qin et al. [38], which suggested that the calculations performed in this study were credible.

Response Mechanism of Soil Erosion Change to VFC and Precipitation
VFC and precipitation are the most direct and important factors affecting soil erosion change, and their annual variation ranges are much larger than other factors [37,[41][42]. Understanding the driving relationship between them and soil erosion change can help deepen our understanding of the processes and mechanisms of soil erosion change.
Between 2000 and 2015, the VFC of the Loess Plateau increased as a whole, showing obvious regional differences. In the central part of the Loess Plateau, the VFC values of the "hilly and gully region" and the "plateau gully region" showed a significant upward trend. However, the VFC had a declining trend in most of the "plain areas." some regions of the "earth-rocky mountainous," "sandy areas" and "irrigated regions" (Figure 7a). From 2000 to 2015, precipitation in most areas of the Loess Plateau showed an upward trend, with the biggest increases in the "hilly and gully region" and the "Gully region of Loess Plateau". Only the southeastern and western marginal regions of the Loess Plateau showed a downward trend (Figure 7b).  Through the Student's t-test, a significant relationship between the annual variations of the VFC and precipitation with the soil erosion modulus on the Loess Plateau was obtained. According to Figure 8a, the VFC had a significant negative correlation with the soil erosion modulus in the southeastern part of the Loess Plateau (−0.05 < P < −0.01), where very significant negative correlations (−0.01 < P < −0.001) and highly significant negative correlations (−0.001 < P < 0) were seen for some areas. Hence, in these areas, the improvement of vegetation coverage restrained the occurrence of soil erosion. In the northwest, except for a small part of the region that had a significant positive correlation (0 < P < 0.05), there were no very significant correlations. and precipitation with the soil erosion modulus on the Loess Plateau was obtained. According to Figure 8a, the VFC had a significant negative correlation with the soil erosion modulus in the southeastern part of the Loess Plateau (−0.05 < P < −0.01), where very significant negative correlations (−0.01 < P < −0.001) and highly significant negative correlations (−0.001 < P < 0) were seen for some areas. Hence, in these areas, the improvement of vegetation coverage restrained the occurrence of soil erosion. In the northwest, except for a small part of the region that had a significant positive correlation (0 < P < 0.05), there were no very significant correlations. Figure 8b shows the relationships between the interannual variations of the precipitation and the soil erosion modulus. There was a significant positive correlation between them (0 < P < 0.05) in most areas of the northwestern Loess Plateau. In addition, most were very significantly positively correlated (0.001 < P < 0.01) or highly significantly positively correlated (0 < P < 0.001). This meant that precipitation controlled the amount of soil erosion modulus in the area as the soil erosion modulus will be the same whether the precipitation increases or decreases. Precipitation and soil erosion modulus did not show a very significant correlation in most areas in the southeastern part of the Loess Plateau. Based on the t-test results of the VFC, precipitation and soil erosion modulus and their interannual variations (Figures 7 and 8), we can make the following conclusions. From 2000 to 2015, the change of soil erosion modulus in the southeastern part of the Loess Plateau was mainly affected by the VFC, which showed a negative correlation. The decrease of soil erosion modulus in the "hilly and gully region" and the "plateau gully region" of the Loess Plateau were due to improved VFC. Moreover, precipitation mainly affected the soil erosion modulus in the northwest part of the Loess Plateau, which showed a positive correlation. The increase of soil erosion modulus in some areas in the "sandy areas" and "irrigated regions" is due to an increase of precipitation.

Contribution Rate of Human Activities
From 2000 to 2015, human activities on the Loess Plateau contributed negatively to soil erosion, with a contribution rate of −1.0774 (Table 3). An analysis of variance (ANOVA) was carried out to describe the impact of human activities and climate change on the variation of soil erosion rates, where 6210 points were randomly obtained in the study area with the sampling rate of 1%. The results show that human activities and climate change have significantly different effects on the soil erosion rate (P < 0.01), where the effect of the former is significantly greater than the latter on the decline of erosion rates ( Table 4).
The contribution rate of human activities is distinct in terms of its spatial differentiation. Most of the northwestern part of the Loess Plateau had a negative contribution, that is, human activities  Figure 8b shows the relationships between the interannual variations of the precipitation and the soil erosion modulus. There was a significant positive correlation between them (0 < P < 0.05) in most areas of the northwestern Loess Plateau. In addition, most were very significantly positively correlated (0.001 < P < 0.01) or highly significantly positively correlated (0 < P < 0.001). This meant that precipitation controlled the amount of soil erosion modulus in the area as the soil erosion modulus will be the same whether the precipitation increases or decreases. Precipitation and soil erosion modulus did not show a very significant correlation in most areas in the southeastern part of the Loess Plateau.
Based on the t-test results of the VFC, precipitation and soil erosion modulus and their interannual variations (Figures 7 and 8), we can make the following conclusions. From 2000 to 2015, the change of soil erosion modulus in the southeastern part of the Loess Plateau was mainly affected by the VFC, which showed a negative correlation. The decrease of soil erosion modulus in the "hilly and gully region" and the "plateau gully region" of the Loess Plateau were due to improved VFC. Moreover, precipitation mainly affected the soil erosion modulus in the northwest part of the Loess Plateau, which showed a positive correlation. The increase of soil erosion modulus in some areas in the "sandy areas" and "irrigated regions" is due to an increase of precipitation.

Contribution Rate of Human Activities
From 2000 to 2015, human activities on the Loess Plateau contributed negatively to soil erosion, with a contribution rate of −1.0774 (Table 3). An analysis of variance (ANOVA) was carried out to describe the impact of human activities and climate change on the variation of soil erosion rates, where 6210 points were randomly obtained in the study area with the sampling rate of 1%. The results show that human activities and climate change have significantly different effects on the soil erosion rate (P < 0.01), where the effect of the former is significantly greater than the latter on the decline of erosion rates (Table 4).  The contribution rate of human activities is distinct in terms of its spatial differentiation. Most of the northwestern part of the Loess Plateau had a negative contribution, that is, human activities inhibited the occurrence of soil erosion. The southeastern region was mainly characterized by a positive contribution, that is, human activities exacerbated the occurrence of soil erosion (Figure 9). inhibited the occurrence of soil erosion. The southeastern region was mainly characterized by a positive contribution, that is, human activities exacerbated the occurrence of soil erosion (Figure 9). Here, under the premise that the average soil erosion modulus of each region is reduced, C < −1 indicates that human activities have a negative contribution to the increase of soil erosion modulus while climate change promotes soil erosion. When −1 ≤ C < 0, it indicates that human activities cooperated with climate change to inhibit the increase of the soil erosion modulus. When C > 0, it indicates that human activities inhibits climate change and promotes the increase of soil erosion modulus.  Regionally, the "sandy areas" in the northwest of the Loess Plateau had the greatest negative contribution rate to human activities, with an average contribution rate of −4.9807. This indicated that Regionally, the "sandy areas" in the northwest of the Loess Plateau had the greatest negative contribution rate to human activities, with an average contribution rate of −4.9807. This indicated that climate change in this region should had promoted the increase of soil erosion modulus, but human activities made the soil erosion modulus decrease. This may have arisen due to the limited vegetation growth conditions in the "sandy areas" and human activities in the region, including tree and grass planting and even agricultural activities which made the conditions for vegetation growth far better than they would normally be in a typical desert climate.
The "irrigated regions" had the second highest negative contribution of human activities, with an average contribution rate of −1.2717. In the process of agricultural production in the "irrigated areas", agricultural measures such as irrigation and fertilization of cultivated land might result in the conditions of growing crops, fruit trees and other plants being better than those of vegetation under natural conditions in the area.
The average contribution rates of human beings in the "hilly and gully region" and the "plateau gully region" of the Loess Plateau are −0.5513 and −0.7805, respectively. Soil erosion declined most obviously from 2000 to 2015 in these two regions, which are the main implementation areas of the "Grain for Green Project" and the comprehensive control of small watersheds on the Loess Plateau. It can be concluded that the improvement of vegetation conditions by ecological engineering has played an important role in controlling the occurrence of soil erosion.
The contribution rate of human activities was positive in the "earth-rocky mountainous" and the "plain areas" in the southeastern part of the Loess Plateau, with average contribution rates of 0.2345 and 0.3929, respectively. These indicated that human activities in these two regions promoted the occurrence of soil erosion. The region with the largest positive contribution rate of human being was located in the Guanzhong plain urban agglomeration, where a dense population, developed industry and agriculture, frequent human production and living activities had destroyed surface vegetation and intensified the occurrence of soil erosion.

Spatial Differentiation of the Soil Erosion Driving Mechanism and the Contribution Rate of Human Activities
Several areas had significant correlations between the VFC or precipitation change and soil erosion modulus change on the Loess Plateau. Additionally, each area either had a positive or negative contribution rate of human activities to the soil erosion modulus, although all showed similar spatial characteristics. The spatial differentiation between the southeast and northwest parts of the Loess Plateau was obvious, and a clear demarcation line from northeast to southwest emerged.

Demarcation Line of the Driving Mechanism
Using the precipitation data of the Loess Plateau from 2000 to 2015 to calculate the annual average precipitation, it is found that the 400-mm isohyetal line roughly coincides with the demarcation line of the significant influencing areas of VFC and precipitation on soil erosion on the Loess Plateau ( Figure 10). The northwest, which has an annual average precipitation of less than 400 mm, is mainly affected by precipitation, while the southeast, which has an annual average precipitation of more than 400 mm, is mainly affected by VFC. The 400-mm isohyetal line is an important geographical boundary line in China, which distinguishes semi-humid and semi-arid areas, forest vegetation and grassland vegetation [43]. The southeast part of the 400-mm isohyetal line belongs to semi-humid areas with relatively better vegetation conditions, where precipitation is not a limiting factor for vegetation growth [44]. Therefore, rainfall hardly affects the soil erosion or vegetation growth, but rather the surface vegetation coverage is the main factor affecting vegetation growth in this region. However, the northwest part of 400-mm isohyetal line contains semi-arid areas, where grasslands and shrub vegetation are the main vegetation types. Arid climatic conditions make precipitation a limiting factor for vegetation growth, which can result in poorer vegetation conditions; hence, precipitation is a more effective inhibiting factor than VFC in terms of soil erosion in arid regions. Remote Sens. 2019, 11, x FOR PEER REVIEW 15 of 20 Figure 10. Correlation between VFC and precipitation and the erosion modulus determined with the student's t-test.

Demarcation Line of Human Activities Contribution Rate
The demarcation line dividing the positive and negative contribution rates of human activities on soil erosion was roughly consistent with the line proposed by Chinese geographer Hu Huanyong that expressed the distribution of population density ( Figure 11) [45], which ran from Aihui, Heilongjiang Province (renamed Heihe City in 1983) to Tengchong, Yunnan Province. To the southeast of Hu's Line, 36% of China's territory and 96% of the population were distributed, while the northwest contained 64% of China's territory and only 4% of the population. To the southeast of Hu's Line, where the Guanzhong plain urban agglomeration consists of Tianshui City of Gansu Province, Baoji City, Xianyang City, Xi'an City, Shangluo City and Weinan City of Shaanxi Province, Yongji City, Yuncheng City, Hejin City and Linyi City of Shanxi Province is located in, the population density was high and human activities mainly focused on production and living activities such as the exploitation and utilization of natural and ecological resources. The land use types are mainly construction land and production land, which destroyed surface vegetation and promoted the deterioration of the ecological environment. Therefore, the occurrence of soil erosion was aggravated, as demonstrated by the positive increase in the soil erosion modulus. On the contrary, the population density northwest of Hu's Line on the Loess Plateau was small, where both the living space and resources needed by human beings are small, and the ecological restoration and large-scale control projects have been implemented in this region, which include forests and grasslands, and which has resulted in a net inhibition of soil erosion.

Demarcation Line of Human Activities Contribution Rate
The demarcation line dividing the positive and negative contribution rates of human activities on soil erosion was roughly consistent with the line proposed by Chinese geographer Hu Huanyong that expressed the distribution of population density ( Figure 11) [45], which ran from Aihui, Heilongjiang Province (renamed Heihe City in 1983) to Tengchong, Yunnan Province. To the southeast of Hu's Line, 36% of China's territory and 96% of the population were distributed, while the northwest contained 64% of China's territory and only 4% of the population. To the southeast of Hu's Line, where the Guanzhong plain urban agglomeration consists of Tianshui City of Gansu Province, Baoji City, Xianyang City, Xi'an City, Shangluo City and Weinan City of Shaanxi Province, Yongji City, Yuncheng City, Hejin City and Linyi City of Shanxi Province is located in, the population density was high and human activities mainly focused on production and living activities such as the exploitation and utilization of natural and ecological resources. The land use types are mainly construction land and production land, which destroyed surface vegetation and promoted the deterioration of the ecological environment. Therefore, the occurrence of soil erosion was aggravated, as demonstrated by the positive increase in the soil erosion modulus. On the contrary, the population density northwest of Hu's Line on the Loess Plateau was small, where both the living space and resources needed by human beings are small, and the ecological restoration and large-scale control projects have been implemented in this region, which include forests and grasslands, and which has resulted in a net inhibition of soil erosion. Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 20 Figure 11. Human activity contribution rate of soil erosion modulus change and the location of Hu's Line (1935).

Regional Characteristics of Soil Erosion Control on the Loess Plateau under the Combined Influence of Human and Climate
Population distribution and the natural environment are not independent of each other, as indicated by the spatial coincidence of the 400-mm isohyetal line and Hu's Line. Climate conditions affected the natural environment of the Loess Plateau, and the natural environment determines the human population density and their way of production and lifestyle, which together had a combined effect on soil erosion on the Loess Plateau. The comprehensive influence of these processes led to the formation of the southeast-northwest spatial differentiation characteristic of the Loess Plateau in terms of soil erosion. In our study the Loess Plateau was divided into three parts: (1) Semi-arid area, where planting grass and drought-tolerant plants, building terraces and other engineering measures may better achieve soil and water conservation. This is consistent with the conclusion of Qiao's research that the Soil Moisture Content (SMC) should be kept in mind when carrying out revegetation in semi-arid regions on the Loess Plateau [46]. (2) Afforestation areas, which are key areas for bioengineering and are suitable for afforestation. (3) Urban area, i.e. the main areas of human activities, which requires attention to prevent soil erosion at the same time of urban expansion, and hence achieve sustainable development. (Figure 12) Wang et, al. [47] obtained similar zoning results through an analysis of rainfall. In addition to rainfall, our study considers a variety of factors, including human activities and vegetation, which represents a more robust scientific result. Population distribution and the natural environment are not independent of each other, as indicated by the spatial coincidence of the 400-mm isohyetal line and Hu's Line. Climate conditions affected the natural environment of the Loess Plateau, and the natural environment determines the human population density and their way of production and lifestyle, which together had a combined effect on soil erosion on the Loess Plateau. The comprehensive influence of these processes led to the formation of the southeast-northwest spatial differentiation characteristic of the Loess Plateau in terms of soil erosion. In our study the Loess Plateau was divided into three parts: (1) Semi-arid area, where planting grass and drought-tolerant plants, building terraces and other engineering measures may better achieve soil and water conservation. This is consistent with the conclusion of Qiao's research that the Soil Moisture Content (SMC) should be kept in mind when carrying out re-vegetation in semi-arid regions on the Loess Plateau [46]. (2) Afforestation areas, which are key areas for bioengineering and are suitable for afforestation. (3) Urban area, i.e. the main areas of human activities, which requires attention to prevent soil erosion at the same time of urban expansion, and hence achieve sustainable development ( Figure 12).
Wang et, al. [47] obtained similar zoning results through an analysis of rainfall. In addition to rainfall, our study considers a variety of factors, including human activities and vegetation, which represents a more robust scientific result. Remote Sens. 2019, 11, x FOR PEER REVIEW 17 of 20

Deficiencies and Prospects in Soil-Erosion Research
In this study, human activities in farming areas showed a negative contribution to the increase of soil erosion modulus; that is, agricultural production activities could inhibit the soil erosion. In addition, the fine identification of agricultural measures could improve the accuracy and validity of the obtained results. For example, in this study, we used RUSLE to estimate soil erosion, where the C factor in our model is only based on VFC, which ignored the great impact of vertical structure of vegetation on soil erosion. The arbor layer, shrub layer, herb layer and litter could effectively intercept precipitation and reduce the rain kinetic energy. Studies have shown that the litter layer above 1 cm can reduce the occurrence of splash erosion by more than 90% [48]. Splash erosion is the main cause of soil loss in semi-arid regions with low rainfall intensities [49]. Therefore, it is necessary to further identify and study the underlying conditions of farmland and orchard vegetation to determine whether they perform the same or similar roles related to soil consolidation as grasslands and forests.
In addition, it is an indisputable fact that human-induced changes in land use patterns will indirectly affect regional climates [46,50]. The impact of human activities on soil erosion in this study refers to the direct intervention of humans on the process of soil erosion, while ignoring the complex human-land-atmosphere coupling behind it. The contribution rate of human activities calculated in this research is not a minimal aspect of the destruction of the natural environment by human production and living activities or the improvement of the natural environment by the implementation of ecological protection or restoration works, but a comprehensive contribution to environmental protection, good or bad. Generally, the actual benefits of ecological projection are higher than the contribution rate of human activities simulated in this study if a certain area contains human life activities. If the land parcel and implementation of the project area can be accurately determined, we can further separate the contribution rate of ecological engineering or human production and living activities by refining pure engineering plots, which is our future research focus.

Deficiencies and Prospects in Soil-Erosion Research
In this study, human activities in farming areas showed a negative contribution to the increase of soil erosion modulus; that is, agricultural production activities could inhibit the soil erosion. In addition, the fine identification of agricultural measures could improve the accuracy and validity of the obtained results. For example, in this study, we used RUSLE to estimate soil erosion, where the C factor in our model is only based on VFC, which ignored the great impact of vertical structure of vegetation on soil erosion. The arbor layer, shrub layer, herb layer and litter could effectively intercept precipitation and reduce the rain kinetic energy. Studies have shown that the litter layer above 1 cm can reduce the occurrence of splash erosion by more than 90% [48]. Splash erosion is the main cause of soil loss in semi-arid regions with low rainfall intensities [49]. Therefore, it is necessary to further identify and study the underlying conditions of farmland and orchard vegetation to determine whether they perform the same or similar roles related to soil consolidation as grasslands and forests.
In addition, it is an indisputable fact that human-induced changes in land use patterns will indirectly affect regional climates [46,50]. The impact of human activities on soil erosion in this study refers to the direct intervention of humans on the process of soil erosion, while ignoring the complex human-land-atmosphere coupling behind it. The contribution rate of human activities calculated in this research is not a minimal aspect of the destruction of the natural environment by human production and living activities or the improvement of the natural environment by the implementation of ecological protection or restoration works, but a comprehensive contribution to environmental protection, good or bad. Generally, the actual benefits of ecological projection are higher than the contribution rate of human activities simulated in this study if a certain area contains human life activities. If the land parcel and implementation of the project area can be accurately determined, we can further separate the contribution rate of ecological engineering or human production and living activities by refining pure engineering plots, which is our future research focus.

Conclusions
From 2000 to 2015, the soil erosion modulus of the Loess Plateau showed a downward trend as a whole, with a rate of −0.6408 t/ha/a, although the downward trend gradually slowed down. Precipitation mainly affected the change of soil erosion modulus in the northwestern part of the Loess Plateau, as indicated by the significant positive correlation between them. The VFC mainly affected in the southeastern part of the Loess Plateau, which displayed a significant negative correlation. Within 16 years, human activities had a negative contribution to the change of soil erosion modulus, with a rate of −1.0774. The ANOVA result of the combined rate induced by human activities and climate change on the soil erosion rates was very significant (P < 0.01); that is, human activities effectively reduced the soil erosion modulus while climate change promoted soil erosion. Regionally, both the annual average soil erosion modulus and decline rate in the "hilly and gully region" and the "gully region of Loess Plateau" were the largest. The average human activity rate in these two regions were both negative, with rates of −0.5513 and −0.7805 respectively, which indicated that though climate change helped, human activities effectively controlled the occurrence of soil erosion by improving the surface vegetation conditions through ecological projects.
The soil erosion driving mechanisms and human contribution rates on the Loess Plateau showed a boundary line from northeast to southwest, which can be well explained by the identified 400-mm isohyetal line and the Hu's Line. Finally, the Loess Plateau was divided into three parts, semi-arid area, afforestation area and urban areas, where different areas should implement different measures that meet the regional characteristics to effectively prevent soil erosion.
Author Contributions: Q.S. and X.G. conceived and designed the experiments; X.G. analyzed the data and wrote the paper.
Funding: This research was funded by the "National Key R&D Program of China, grant number 2017YFC0506501", and "'Beautiful China' Ecological Civilization Construction Science and Technology Project, grant number XDA23100203".