Spatial–Temporal Changes and Driving Force Analysis of Ecosystems in the Loess Plateau Ecological Screen

Abstract

The ecological degradation caused by unreasonable development and prolonged utilization threatens economic development. In response to the development crisis triggered by ecological degradation, the Chinese government launched the National Barrier Zone (NBZ) Construction Program in 2006. However, few in-depth studies on the Loess Plateau Ecological Screen (LPES) have been conducted since the implementation of that program. To address this omission, based on the remote sensing image as the primary data, combined with meteorological, soil, hydrological, social, and economic data, and using GIS spatial analysis technology, this paper analyzes the change characteristics of the ecosystem pattern, quality, and dominant services of the ecosystem in the LPES from 2005 to 2015. The results show that from 2005 to 2015, the ecosystem structure in the study area was relatively stable, and the area of each ecosystem fluctuated slightly. However, the evaluation results based on FVC, LAI, and NPP showed that the quality of the ecosystem improved. The vegetation coverage (FVC) increased significantly at a rate of 0.91% per year, and the net primary productivity (NPP) had increased significantly at a rate of 6.94 gC/(m²∙a) per year. The leaf area index (LAI) in more than 66% of the regions improved, but there were still about 8% of the local regions that were degraded. During these 10 years, the soil erosion situation in LPES improved overall, and the amount of soil conservation (ASC) of the ecosystem in the LPES increased by about 0.18 billion tons. Grassland and forest played important roles in soil conservation in this area. Pearson correlation analysis and redundancy analysis showed that the soil conservation services (SCS) in the LPES were mainly affected by climate change, economic development, and urban construction. The precipitation (P), total solar radiation (SOL), and temperature (T) can explain 52%, 30.1%, and 17% of the change trends of SCS, respectively. Construction land and primary industry were negatively correlated with SCS, accounting for 22% and 8% of the change trends, respectively. Overall, from 2005 to 2015, the ecological environment of LPES showed a gradual improvement trend, but the phenomenon of destroying grass and forests and reclaiming wasteland still existed.

1. Introduction

The healthy development of human society depends on the sustainable and stable provision of ecosystem services (ES) [1]. Many signs indicate that human demand for various ES will continue to increase substantially in the coming decades [2]. According to current estimates, by 2050, the world’s population will increase by more than 3 billion and the world economy will quadruple, resulting in a sharply increase in human demand for and consumption of biological and material resources [3]. Meanwhile, the impact of unreasonable human activities such as deforestation, grassland reclamation, overgrazing, and excessive digging on the ecosystem has accelerated the ecosystems to degenerate and reduce the ES. Over the past quarter-century, China’s economy has expanded 10-fold [3], which is considered to be at the cost of eco-environmental damage. In the face of increasingly acute ecological problems, ecological protection and green development have become the consensus of the whole society in China: lucid waters and lush mountains are invaluable assets [3]. Protecting and restoring ecosystems to maintain and enhance ES significantly has become a critical work content of all Chinese government levels. In the 21st century, the Chinese government is implementing more stringent environmental protection policies and larger-scale ecological construction/restoration projects [4].

In 2006, the Chinese government started implementing an essential national ecological security strategic pattern construction plan—the national barrier zone (NBZ), which aims to achieve national ecological security by protecting and constructing some key areas [5]. It is composed of the Tibet Plateau ecological barrier, the Loess Plateau–Sichuan–Yunnan ecological barrier, the northeast forest belt (NFB), the northern sand-stabilization belt (NSB), and the southern hilly land zone. The NBZ includes 21 provinces (cities), including Fujian Province, Gansu Province, Guangdong Province, Hebei Province, and Heilongjiang Province, etc., reaching a total area of 3,135,200 km². This plan aims to improve the quality and services of different ecosystems in many key ecological regions in China through targeted protection and restoration. It will be able to better play the role of ecological barriers as “filtering”, “buffering”, and “shelter”, which is of far-reaching significance for China to achieve sustainable socio-economic development and high-quality. Among them, the Loess Plateau Ecological Screen (LPES), as a vital component of the NFB, is not only an energy and mineral concentration area in China but also a substantial barrier to curb desertification, control soil erosion, and protect the agricultural areas in north China [6,7].

Many studies have been conducted on the Loess Plateau and achieved fruitful research results [8,9,10,11]. However, there are few studies related to the LPES. The LPES has been reported systematically only in the “National Ecological Environment Ten-Year Change Survey, and Evaluation” project jointly carried out by the Ministry of Environmental Protection (now the Ministry of Ecology and Environment) and the Chinese Academy of Sciences in 2012 [3]. Using remote sensing survey and sample plot survey, the project systematically and comprehensively studied and analyzed the eco-environment of China for the first time at the national scale, i.e., the pattern, quality, ES characteristics and temporal and spatial changes of the ecosystems of the NBZ. Nevertheless, the threat of global climate change to the eco-environment is increasing [12,13,14]. The original research results are insufficient to support new decision-making forms, and a new in-depth analysis is needed to scientifically evaluate the status of the ecosystems in the NBZ, especially in the LPES, after ten years of construction.

Based on remote sensing as the primary data, combined with meteorological, soil, hydrology, social, and economic data, this study aims to use GIS spatial analysis technology to systematically analyze the evolution characteristics of the pattern, quality, and ES of LPES ecosystem from 2005 to 2015. This study can provide decision support and policy reference for the construction of the LPES.

2. Materials and Methods

2.1. Study Area

The LPES is located in the range of $105.1°–112.21° \mathrm{E}$, $34.01°–38.13° \mathrm{N}$, covering an area of about 12.11 $\times {10}^4{ \mathrm{km}}^2$, ranging from 1200 to 1600 m above sea level, accounting for 19.5% of the total area of the Loess Plateau (Figure 1). It involves about 60 counties Shanxi, Shaanxi, Gansu, and Ningxia Hui Autonomous Region. The region has a semi-arid continental climate with an average annual rainfall of 400–776 mm [15]. The rainfall is concentrated mainly between July and September and varies significantly from year to year [9]. The soil in the study area is mainly derived from loess, and it is fine silt to silt in texture and weakly resistant to erosion [10]. The LPES is not only a sensitive area for climate change but also a key area for soil erosion control in the middle and upper reaches of the Yellow River [16,17,18,19], and an essential area for the Chinese government to implement the project of returning farmland to forest and grassland [20].

2.2. Data

Many data types are involved in this study, including meteorological data, digital elevation model, soil data, remote sensing data, etc. The main data required for the study are shown in

Table 1. List of data.

Data Name	Description	Source
Meteorological data	Temperature, precipitation, relative humidity, sunshine hours, etc.	National Meteorological Information Center
Harmonized World Soil Database_China_Subset	Soil texture, including sand, silt, clay, and organic carbon	National Tibetan Plateau/Third Pole Environment Data Center
DEM	90 m	Resource and Environment Science and Data Center, Chinese Academy of Sciences
Land cover	250 m	NASA
MOD13Q1 NDVI	250 m	NASA
MOD11A12 LST	250 m	NASA
MOD15A2 LAI	250 m	NASA

. Meteorological data are derived from the National Meteorological Information Center (http://data.cma.cn, last accessed on 12 September 2021), including temperature, relative humidity, and sunshine hours at 52 stations. Before use, spatial interpolation of meteorological data generates raster data with a spatial resolution of 250 m. Topographical factors are extracted using 30 m SRTM-DEM data products from the Resource and Environment Science and Data Center, Chinese Academy of Sciences.

The land cover, LAI, LST, and NDVI Products of Moderate Resolution Imaging Spectroradiometer (MODIS) are the essential data to evaluate the ecosystem quality and service in the study area, derived from the U.S. National Aeronautics and Space Administration (NASA). The MODIS Land Cover Product includes 11 natural vegetation types, 3 land development and mosaic types and 3 non-grassland types. After preprocessing, such as projection transformation and mosaic, the ecosystem types in the study area were extracted according to the IGBP land cover classification system [20]. The land cover classification system of IGBP is shown in

Table 2. Classification system of first-class of LPES.

Ecosystem	Description
Forest	deciduous broad-leaved forest, evergreen broad-leaved forest, deciduous needle-leaf forest, evergreen coniferous forest, mixed broadleaf-conifer forest, open forest
Grassland	grassland, meadow, hassock
Waterbody	river course, lake, bog
Cropland	cultivated land, garden land
Built-up land	residence, road, urban green space
Desert	desert
Others	bare land, sandy soil, saline-alkali land,

. MODIS Reprojection Tools (MRT) were used for mosaics, format conversion, and projection of MODIS NDVI products in each ten-day data. Additionally, then, the maximum NDVI of the whole year was calculated by the maximum value synthesis method, which can reflect the best growth condition of vegetation in one year, and further eliminate the influence of factors such as atmosphere and solar elevation angle [21,22]. Near-infrared and short-wave infrared data were extracted from the MODIS LST Product as parameters for estimating NPP. The MODIS LAI product incorporated the MODIS Land Cover Product and MODIS Surface Reflectance Product into the radiative transfer model to derive LAI.

2.3. Methodology

2.3.1. Analysis of the Evolution of Ecosystem Pattern

Transition Matrix

The ecosystem transfer matrix can analyze the structural characteristics of regional ecosystem changes and the direction of various changes [23]. It can directly reflect the classification structure of ecosystems in the initial and final stages of the research, but it can also reflect the transformation of various ecosystem types during the research period. The conversion relationship between the ecosystem types in the study area at the beginning and the end of a particular period is calculated using the following formula [24]:

$$S_{ij}=\left[\begin{array}{ccccc}{S}_{11}& S_{12}& S_{13}& \cdot \cdot \cdot & S_{1n}\\ S_{21}& S_{22}& S_{23}& \cdot \cdot \cdot & S_{2n}\\ S_{31}& S_{32}& S_{33}& \cdot \cdot \cdot & S_{3n}\\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ S_{n1}& S_{n2}& S_{n3}& \cdot \cdot \cdot & S_{nn}\end{array}\right]$$

(1)

where: n is the number of ecosystems, i and j are the ecosystems at the beginning and end of the research period, $S_{ij}$ is the area where the i-th ecosystem is converted to the j-th ecosystem during the study period.

Land Cover Change Index

The land cover change index ($LCC{I}_{ij}$) quantifies the transformation intensity of ecosystem types [25,26]. When the value of $LCC{I}_{ij}$ is positive, it means that the overall ecosystem type of the study area has improved, while the value of $LCC{I}_{ij}$ is negative, indicating that the overall ecosystem type in this study area has deteriorated [19]. The calculation formula is:

$$LCC{I}_{ij}=\frac{\sum \left|A_{ij}\times \left(D_a-D_b\right)\right|}{A_{ij}}\times 100\%$$

(2)

where: i is the research area, j is the ecosystem type, $j=1,\cdots ,n$, $A_{ij}$ is the area of a specific ecosystem type converted once, $D_a$ is the level before the conversion and $D_b$ is the conversion post-class level. The classification standard of ecosystems is shown in

Table 3. Classification standard of various ecosystems.

Degree		Ecosystem	Level
High		Waterbody	level 1
		Forest	level 2
		Grassland	level 3
		Cropland	level 4
		Built-up land	level 5
		Desert	level 6
Low		Others	level 7

[27].

PNTIL Model

The forward and reverse transformation index (PNTIL) models can analyze the evolution trend of the ecosystem patterns [4,27]. The conversion rules are shown in

Table 4. Rules of forwarding/reverse ecosystem transformation.

Ecosystem	Converted Type	Transformation Direction
Forest (I)	III	+
	I, II, IV, V, VI, VII	-
Grassland (II)	I, III	+
	IV, V, VI, VII	-
Water body (III)	-	+
	I, II, IV, V, VI, VII	-
Cropland (IV)	I, II, III	+
	V, VI, VII	-
Built-up land (V)	I, II, III, IV	+
	VI, VII	-
Desert (VI)	I, II, III, IV, V, VII	+
	-	-
Others (VII)	I, II, III, IV, V,	+
	VI	-

. The higher the forward conversion rate, the more stable the ecosystem, while the higher the reverse conversion rate, the more vulnerable the ecosystem [5,28].

$$LCT{R}_{i,k}=\left(\Delta S_{i,k}/S_i\right)\times 1/t\times 100$$

(3)

where: $LCT{R}_{i,k}$ is the transformation rate of ecosystem types in the i-th area to k directions, $S_i$ is the total area of all ecosystems in the i-th area, $\Delta S_{i,k}$ is the total area transformed to type k, k = 1. The positive transformation rate of the ecosystem in the region, k = 2 indicates the reverse transformation rate of the ecosystem in the region, $t$ is the time variable.

2.3.2. Analysis of the Ecosystem Quality

Vegetation Coverage

Vegetation Coverage (FVC) is the most commonly used and most important indicator to measure the quality of ecosystems and the status of surface vegetation [7,29]. With the development of remote sensing technology, remote sensing images are often used to calculate FVC at a regional scale. Usually, the FVC of a specific area is calculated by the Normalized Difference Vegetation Index (NDVI) [30]. The calculation formulas of NDVI and FVC are as follows:

$$NDVI=\left({\rho }_{NIP}-{\rho }_R\right)/\left({\rho }_{NIP}+{\rho }_R\right)$$

(4)

where: ${\rho }_{NIP}$ and ${\rho }_R$ are, respectively, the reflectance of the near-infrared band and red band.

$$FVC=\frac{NDVI-NDV{I}_{soil}}{NDV{I}_{veg}-NDV{I}_{soil}}$$

(5)

where:$FVC$ is the vegetation coverage,$NDV{I}_{veg}$ is the$NDVI$ value of pure vegetation pixels, and $NDV{I}_{soil}$ is the$NDVI$ value of no vegetation pixels, which is the soil background value.

Net Primary Productivity

Net primary productivity (NPP) refers to the total dry organic matter produced by plants per unit time and unit area [31,32]. NPP can directly characterize the productivity of vegetation and reflect the quality status of the ecosystem. At present, there are many models for calculating NPP. Due to the convenience of remote sensing data acquisition, ecological models based on remote sensing data are applied at regional and global scales. The CASA (Carnegie–Ames Standford Approach) model is the most widely used in the world, and it is a process model based on light use efficiency (LUE) [33,34]. Therefore, this paper uses the CASA model proposed by Potter et al. (1993) to estimate the NPP of LPES [35]. The calculation formula is:

$$NPP=APAR\left(t\right)\times \epsilon \left(t\right)$$

(6)

$$APAR=FPAR\times PAR$$

(7)

where: $APAR$ is the photosynthetic active radiation absorbed by the plant; $\epsilon$ is the efficiency of converting $APAR$ into organic carbon; $FPAR$ is the effective absorption ratio of $APAR$ by the plant; $PAR$ is the driving energy for photosynthesis of the vegetation. The calculation method of each parameter is as follows:

$$PAR=0.48\times K_{24}^{\downarrow }\left(t\right)$$

(8)

$$K_{24}^{\downarrow }\left(t\right)=\left\{0.25+\frac{0.50n\left(t\right)}{N\left(t\right)}\right\}{K}_{24}^{\downarrow exo}\left(t\right)$$

(9)

$$K_{24}^{\downarrow exo}\left(t\right)=\frac{24\times 60}{\pi }{G}_{sc}{d}_r\left[{\omega }_s\sin \left(\phi \right)\sin \left(\delta \right)+\cos \left(\phi \right)\cos \left(\delta \right)\sin \left({\omega }_s\right)\right]$$

(10)

$$\delta =0.409\sin \left(\frac{2\pi }{365}J-1.39\right)$$

(11)

$$d_r=1+0.033\cos \left(\frac{2\pi }{365}J\right)$$

(12)

$${\omega }_s=\arccos \left[-\tan \left(\phi \right)\tan \left(\delta \right)\right]$$

(13)

where: $K_{24}^{\downarrow }\left(t\right)$ is the total solar radiation, $K_{24}^{\downarrow exo}\left(t\right)$ is the extraterrestrial solar radiation $\left(MJ/m^2·t\right)$, $G_{sc}$ is the solar constant, taken as 0.0820 $MJ/m^2·min$ (approximately 1366.67 $W/m^2$), $d_r$ is the relative distance between the sun and the earth, ${\omega }_s$ is the sunset angle (rad), $\phi$ is the Latitude (rad), $\delta$ is the angle (rad) between the earth’s equatorial plane and the line connecting the sun and center of the earth (rad), $J$ is the Julian day, $N\left(t\right)$ is the potential sunshine hours, $n\left(t\right)$ is the actual sunshine hours [36].

$$FPAR=\frac{\left(SR-S{R}_{min}\right)\times \left(FPA{R}_{max}-FPA{R}_{min}\right)}{S{R}_{max}-S{R}_{min}}+FPA{R}_{min}$$

(14)

$$SR=\frac{NIR}{RED}=\frac{1+NDVI}{1-NDVI}$$

(15)

where: the related studies have shown that the maximum and initial values of $FPAR$ have nothing to do with the vegetation type. By calculating the NDVI value of the corresponding vegetation type, the lower percentiles of 5% and 95% are obtained. Among them, $FPA{R}_{min}$ is 0.001 and $FPA{R}_{max}$ is 0.95; $SR$ is the ratio vegetation index, $NIR$ is the reflectivity in the near-infrared band, and $RED$ is the reflectivity in the red band.

$$\epsilon \left(t\right)={\epsilon }^{*}\times T_1\left(t\right)\times T_2\left(t\right)\times W\left(t\right)$$

(16)

$$T_1=0.0005{\left(T_{opt}-20\right)}^2+1$$

(17)

$$T_2=\frac{1}{1+\exp {}^{\left\{0.2\left(T_{opt}-10-T_{mon}\right)\right\}}}\times \frac{1}{1+\exp {}^{\left\{0.3\left(-T_{opt}-10+T_{mon}\right)\right\}}}$$

(18)

$$W\left(t\right)=\frac{EET\left(t\right)}{PET\left(t\right)}$$

(19)

where: ${\epsilon }^{*}$ is the maximum lutilization rate, $g/MJ$, T₁ and T₂ are the effects of temperature and vegetation photosynthesis at different degrees, respectively, W is the stress coefficient of the influence of water on vegetation growth, $T_{opt}$ is the time when the NDVI value of the vegetation reaches the highest value in the growth stage. The average temperature of the current month, °C, $T_{mon}$ is the average monthly temperature, °C, EET is the actual evapotranspiration, mm, and PET is the potential evapotranspiration mm. Among them, T₁, T₂, and W are dimensionless parameters.

Leaf Area Index

The leaf area index (LAI) refers to the total area of plant leaves from a specific land area to the land area [37,38,39]. It can reflect the ecological status of vegetation from the vertical structure and is a critical structural parameter in ecosystem research. There are many LAI measurement methods, but from the perspective of large-scale regional research, remote sensing technology can realize large-scale, dynamic real-time monitoring of regional LAI. The vegetation canopy radiation transmission model was commonly used to obtain LAI using remote sensing technology. This model could establish the relationship between the vegetation canopy structure, leaf optical parameters, and vegetation canopy reflectivity from the optical and structural characteristics [40]. The radiation transmission process is described as follows:

$$-\mu \frac{\partial L\left(z,\Omega \right)}{\partial \tau }+G\left(\tau ,\Omega \right)L\left(z,\Omega \right)=\frac{\omega }{4\pi }{{\displaystyle \int }}_{4\pi }^{}p\left({\Omega }^{\prime }\to \Omega \right)G\left({\Omega }^{\prime }\right)L\left(z,{\Omega }^{\prime }\right)d{\Omega }^{\prime }$$

(20)

where: L is the brightness,$\mu$ is the cosine of the zenith angle in the light radiation transmission direction,$G$ represents the leaf inclination angle distribution function,$\omega$ is the leaf albedo, and $p\left({\Omega }^{\prime }\to \Omega \right)$ is the vegetation canopy phase function.

2.3.3. Soil Conservation Services

The soil conservation service (SCS) of the ecosystem is characterized by the amount of soil conservation (ASC) of vegetation [41,42]. The soil conservation of vegetation is equal to the difference between the potential soil conservation and the actual soil conservation. Among them, the potential soil conservation amount refers to soil erosion under the assumption that there is no vegetation. The actual ASC refers to soil erosion that occurs when vegetation exists [43,44]. In this study, the universal soil loss equation (USLE) was selected to evaluate the soil conservation of the ecosystem in the LPES [45]. The calculation formula is as follows:

$$SC=S{E}_p-S{E}_a=R·K·LS·\left(1-COG\right)$$

(21)

where: SC is the amount of soil conservation,$t/m^2·\mathrm{a}$. $S{E}_p$ and $S{E}_a$ are the potential and actual soil erosion, respectively, $t/(m^2·\mathrm{a})$. R is the rainfall erosivity factor, $MJ·\mathrm{mm}/\left({\mathrm{hm}}^2·\mathrm{h}·\mathrm{a}\right)$. K is the soil erodibility factor, $t·{\mathrm{hm}}^2·\mathrm{h}/\left({\mathrm{hm}}^2·\mathrm{MJ}·\mathrm{mm}\right)$. LS and COG are the topographic factors, and the vegetation coverage factor, respectively, and all of them are dimensionless parameters. The calculation methods of each factor are as follows.

(1)

Rainfall erosivity factor (R)

The R of 52 meteorological stations in and around LEPS was calculated using rainfall data, and then the continuous grid surface of R in the study area was obtained by Kriging interpolation. The formulas of R are as follows.

$$\overline{R}={\displaystyle \sum }_{k=1}^{24}\overline{R_k}$$

(22)

$$\overline{R_k}=\alpha {\displaystyle \sum }_{j=1}^m{P}_j^{\beta }$$

(23)

$$\alpha =21.239{\beta }^{-7.3967}$$

(24)

$$\beta =0.6243+\frac{27.346}{{\overline{P}}_{d12}}$$

(25)

$${\overline{P}}_{d12}=\frac{1}{n}{\displaystyle \sum }_{l=1}^n{P}_{dl}$$

(26)

where: $\overline{R}$ is the R of mean annual (MJ⋅mm/(hm²⋅h⋅a)). $\overline{R_k}$ is the R in the k-th half month (MJ⋅mm/(hm²⋅h)). $P_j$ is the daily rainfall greater than or equal to 12 mm on the j-th day of the k-th half month. α, β are parameters.${\overline{P}}_{d12}$ is the mean daily rainfall greater than or equal to 12 mm. $P_{dl}$ is greater than or equal to 12 mm on the l-th day in the study period. k represents the order of 24 half months for a year (k = 1, 2, …, 24). j represents the number of days with daily rainfall greater than or equal to 12 mm in the kth half month (j = 1, 2, …, m), and l represents the number of days with daily rainfall greater than or equal to 12 mm during the study period (l = 1, 2, …, n).

(2)

Soil erodibility factor (K)

In this paper, the Williams model was used to estimate the soil erodibility factor, and the Zhang Keli model was used to correct the deviation [46,47]. The formulas of K are as follows.

$$\begin{gathered} K_0=\left\{0.2+0.3\exp \left[-0.0256SAN\left(1-\frac{SIL}{100}\right)\right]\right\}\times {\left(\frac{SIL}{CLA+SIL}\right)}^3\times \\ \left[1-\frac{0.25C}{C+\exp \left(3.27-2.98C\right)}\right]\times \left[1-\frac{0.7SNI}{SNI+\exp \left(-5.51+22.9SNI\right)}\right]\times 0.1317 \end{gathered}$$

(27)

$$K=\left(-0.01383+0.51575K_0\right)\times 0.1317$$

(28)

where: K is the soil erodibility factor $(t·{\mathrm{hm}}^2·\mathrm{h}/\left({\mathrm{hm}}^2·\mathrm{MJ}·\mathrm{mm}\right)$). SAN, SIL, CLA, and C are the percentage contents of sand, silt, clay, and organic carbon, respectively. SIN equates to1-SAN/100. The value 0.1317 represents the coefficient of conversion from United States units to international units.

(3)

Topographical factor (LS)

In this paper, the slope factor (S) was calculated using the formula proposed by Liu and the slope length factor (L) was calculated using the formula proposed by Desmet [48,49]. The formulas of LS are as follows.

$$S=\{\begin{array}{c}10.8\sin \theta +0.03 \theta <5°\\ 16.8\sin \theta -0.50 5°\le \theta <10°\\ 21.9\sin \theta -0.96 \theta \ge 10°\end{array}$$

(29)

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

(30)

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

(31)

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

(32)

where: θ is the slope (°). λ is the slope length (m).

(4)

Vegetation cover factor (COG)

According to the experts’ suggestions and references, we assign vegetation coverage factors according to different ecosystem types and vegetation coverage [50,51]. The ecosystem types involved include forests, grasslands, water bodies, farmland, towns, deserts, etc. For waterbody, built-up, desert, and others, the parameters in N-SPECT were assigned 0, 0.01, 0.7 and 0, respectively [52]. Forest and grassland are assigned according to different vegetation coverage (

Table 5. COG of different land cover types.

Ecosystem	Vegetation Coverage (%)
	<10	10~30	30~50	50~70	70~90	>90
Forest	0.1	0.08	0.06	0.02	0.004	0.001
Grass land	0.4	0.22	0.14	0.085	0.040	0.011

Based on the grid layer of each factor, the soil retention (t/(hm²⋅a)) was calculated by ArcGIS.

).

3. Results

3.1. Evolution of Ecosystem Pattern

3.1.1. Ecosystem Structure Characteristics and Its Changes

The results of ecosystem classification in LPES are shown in Figure 2. It can be seen from the figure that grassland and farmland are the main ecosystems. Moreover, the statistical results show that the area of these two ecosystems is about 9.38 $\times {10}^4{ \mathrm{km}}^2$, accounting for more than 74%. Forests are mainly distributed in the east and middle of LPES. Water bodies, deserts, and others are scattered in LPES.

From 2005 to 2015, the structure of the LPES ecosystems was relatively stable, but there were frequent conversions between different types. Figure 3 shows the ecosystem transfer in the LPES from 2005 to 2015.

From 2005 to 2010, farmland was the only ecosystem type whose area was reduced. Cropland was mainly converted to grassland, about 168.04 ${\mathrm{km}}^2$, followed by forest, about 68.16 ${\mathrm{km}}^2$. However, some grasslands and forests are still converted into cropland, and the conversion area is about 71.43 ${\mathrm{km}}^2$. This shows that the Grain for Green Project in the region from 2005 to 2010 has achieved good results. However, while implementing ecological construction, there is still the phenomenon of destroying forests and opening up wastelands. In 5 years, about 63.96 ${\mathrm{km}}^2$ of grasslands were reclaimed as croplands.

From 2010 to 2015, a total area of about 653.88 ${\mathrm{km}}^2$ cropland was mainly converted into built-up land, grasslands, and forests, ranking the primary type of ecosystem transferred out. During this period, the transfer of farmland was the primary driver of urban area increment, and about 284.45 ${\mathrm{km}}^2$ of cropland was converted into built-up land, accounting for 76.16% of the transferred area. Nevertheless, 357.73 ${\mathrm{km}}^2$ grassland and forest were still converted into cropland, with an increase of 286.30 ${\mathrm{km}}^2$ from 2005 to 2010, which was more than 400%. During 2010–2015, the urban area increased steadily, and the transformation rate was higher than 2005–2010.

Between 2005 and 2015, the three types of forest and grassland, grassland and cropland, and built-up land and cropland had more mutual transformations. Generally, the project of reforesting the cultivated land has achieved good results, the phenomenon of destroying forests and opening up wasteland existed in some areas, and this phenomenon increased slightly. Furthermore, we should be alert to the phenomenon of blind expansion of built-up land and the occupation of cropland.

3.1.2. Ecosystem Transformation Intensity

The ecosystem transformation intensity is shown in

Table 6. Land Cover Change Index of the Loess Plateau Ecological Screen (%).

Year	2005~2010	2010~2015
$LCC{I}_{ij}$	7.32	−1.47

. From 2005 to 2010, the $LCC{I}_{ij}$ index is positive, which indicates that the ecosystem has been improving in the past five years. Additionally, the index of 2005–2010 was significantly higher than that of 2010–2015, which indicated that the degree of ecosystem improvement from 2005 to 2010 was higher than that of 2010–2015. The $LCC{I}_{ij}$ index showed a negative number from 2010 to 2015, mainly because of the acceleration of urbanization during that period, and a large amount of forest, grassland, and cropland were converted into built-up land. The results of ecosystem transformation intensity analysis show that from 2005 to 2015, especially from 2010 to 2015, the eco-environment faces certain pressure and degradation to a certain extent.

3.1.3. LULC and Its Changes

Table 7. Forward/reverse transformation rate in the Loess Plateau Ecological Screen in 2005~2010, 2010~2015.

Year	Transformation Direction	Grassland	Built-Up Land	Cropland	Others	Desert	Forest	Waterbody
2005–2010	TFR	0.0020	0.0049	0.0054	0.0000	0.0000	0.0000	0.0000
	RFR	0.0014	0.0000	0.0007	0.0000	0.0000	0.0013	0.0124
2010–2015	TFR	0.0027	0.0267	0.0076	0.0000	0.0000	0.0001	0.0000
	RFR	0.0071	0.0000	0.0069	0.0000	0.0000	0.0082	0.0235

Note: TFR is the forward transformation rate, and RTR is the reverse transformation rate.

shows the forward/reverse conversion rate (TFR and RFR) of the LPES from 2005 to 2010 and from 2010 to 2015. In the two periods, the RFR of water and forest is greater than the TFR, indicating that the two ecosystems were more converted to lower-grade ecosystems. Moreover, The RFR in 2010–2015 was more than twice that in 2005–2010, indicating that the environmental threats faced by water bodies and forests are intensifying. Only the TFR of built-up land was greater than the RFR in two periods in all ecosystems. Moreover, in the case of rapid expansion of built-up land, the TFR in 2010–2015 is about five times that in 2005–2010, which shows that, in urban expansion, the implementation of policies such as reforestation and limiting the scope of economic development activities will have a positive impact on ecosystem change. In the two periods, both TFR and RFR of farmland and grassland show an increasing trend, but RFR increases more due to the transformation of more farmland and grassland into built-up land. From the perspective of TFR and RFR, the RFR of each ecosystem in LPES was generally greater than TFR and showed an increasing trend during 2010–2015, which indicated that there were ecological projects such as De-farming and Reafforestation in this region, but the downward trend of its eco-environment is worrying.

It can be seen from the distribution of transformation rate in Figure 4. From 2005 to 2015, the degree of positive/reverse transformation of the ecosystem in all counties was generally low. The areas with a high degree of positive transformation in the first 5 years are mainly distributed in the middle of the LPES. The Ziwu Mountains and Huang long Mountain are key construction areas for De-farming and Reafforestation, and soil conservation. From 2010 to 2015, the degree of positive transformation in most areas was low, and the environmental improvement was poor. However, in the Baota District of Yan’an City in the middle of LPES, Gujiao City, Jiaocheng County and Fenxi County in the east end have a higher positive transformation.

3.2. Ecosystem Quality

3.2.1. Dynamic Characteristics of FVC

Statistical analysis shows that the average FVC vegetation coverage of the LPES is about 56%. It can be seen from Figure 5 that FVC in this area is gradually decreasing from southeast to northwest. FVC is low in parts of the west side—the Luliang Mountains—the east side of Huanglong Mountain and the west side of the Ziwu Mountains. Areas with high FVC are mainly concentrated in the Lvliang Mountains, Ziwu Mountains, Huanglong Mountains, and the northern edge of the Qinling Mountains, especially in the Ziwu Mountains and Huanglong Mountains. FVC was divided into low, middle, and high levels by the natural breakpoint method. In 2015, the areas covered by vegetation at each level were 3.5 $\times {10}^4{ \mathrm{km}}^2$, 4.3 $\times {10}^4 {\mathrm{km}}^2$, and 4.33 $\times {10}^4{ \mathrm{km}}^2$, accounting for 28.64%, 35.45%, and 35.71% of the total area, respectively.

Statistical analysis shows that the area of FVC at all levels changed slightly in 2005, 2010, and 2015 (Figure 6). In the past ten years, FVC in the LPES has increased at an average annual rate of 0.91%, and FVC increased in the region of about 7.94 $\times {10}^4{ \mathrm{km}}^2$ (Figure 7). The area with reduced vegetation coverage is about 8.20 $\times {10}^4{ \mathrm{km}}^2$, accounting for about 6.76% of the total area. It is mainly in the Jinzhong plain at the eastern foot of the Luliang mountains in Shanxi Province and scattered in the central and western parts of the barrier area, such as Heshui County, Huachi County, Qingcheng County, etc.

3.2.2. Dynamic Characteristics of NPP

In 2015, the NPP per unit area of the LPES was 428.92 $gC/\left({\mathrm{m}}^2·a\right)$, the NPP per unit area of forest was the largest at 500.05$gC/\left({\mathrm{m}}^2·a\right)$, the grassland was 427.23$gC/\left({\mathrm{m}}^2·a\right)$, and the cropland was 397.31$gC/\left({\mathrm{m}}^2·a\right)$. It can be seen from Figure 8 that high-grade areas are mainly distributed in most areas of the LPES, and most concentrated in the middle of the LPES. The areas with low NPP are mainly distributed in the districts and counties on the west side of the Lvliang Mountains, and the east side of Huanglong Mountain in the LPES, FVC in this area is low throughout the year, so the NPP value is also low. NPP was also divided into low, middle, and high levels by the natural breakpoint method, and the classification results are shown in Figure 9. NPP showed a better trend in these 10 years from the classification results. The range of low grades continues to decrease. The middle-level shows a state of fluctuation, but this change is positive because many middle-level regions are transformed into high-level regions.

During 2005 to 2015, the NPP in the LPES showed an overall improvement, increasing at an average annual rate of 6.94$gC/\left({\mathrm{m}}^2·a\right)$, degraded only in some parts (Figure 10). Statistical analysis shows that stable and improved NPP areas in the LPES are 4.36 $\times {10}^4{ \mathrm{km}}^2$ and 4.39 $\times {10}^4{ \mathrm{km}}^2,$ respectively, accounting for 36.02% and 36.35%, respectively. NPP enhancement mainly occurred in the areas with FVC less than 50%, mainly distributed in the eastern side of Huanglong and the counties to the north of Ziwuling and Qinling Mountains, such as Yanchuan County, Yanchang County, Ansai County, Zhidan County, Jingchuan County, Changwu County, Lingtai County, etc. While the degraded area of NPP is about 3.34 $\times {10}^4{ \mathrm{km}}^2$, accounting for 27.62%, mainly distributed in Huanglong Mountain, Ziwu Mountain, and Luliang Mountain.

3.2.3. Dynamic Characteristics of LAI

It can be seen from Figure 11 and Figure 12, in 2015, the LAI of the LPES was mainly 2–4, covering an area of about 7.98 $\times {10}^4 {\mathrm{km}}^2$, accounting for 65.86%, which mainly distributed in western Shanxi, central and southwestern parts of Shaanxi Province, southeastern Gansu Province. The areas with LAI greater than four were mainly distributed in the western part of Taiyuan, the eastern part of Luliang, and the central and southern part of Shaanxi Province, with an area of about 2.39 $\times {10}^4{ \mathrm{km}}^2$, accounting for 19.66%. The area with LAI less than one is mainly located in Shaanxi and Shanxi Province, with an area of about 1.32 $\times {10}^4 {\mathrm{km}}^2$. The LAI of the LPES decreases from southeast to northwest.

Comparing the LAI of different levels (Figure 12), we found that, from 2005 to 2015, the LAI of the LPES has an increasing trend in the area ratios of 2–3 and 4–5, and the area of LAI in 2–3 increases faster than that in 4–5, the area with LAI from 0 to 0.5, 0.5 to 1, and 3 to 4 showed a downward trend, and the area with LAI of 0.5 to 1 had a greater area reduction rate than the other two intervals. The LAI increased first and then decreased in 1–2 and 5–6.

From 2005 to 2015, the LAI in the LPES has improved significantly (Figure 13), with an area of about 8.09 $\times {10}^4 {\mathrm{km}}^2$ improved, accounting for 66.70%, especially in the central and western regions of Shanxi Province, the central region of Shaanxi Province, and the southeast of Gansu Province. The analysis from the perspective of LAI shows that the quality of the ecosystem has improved significantly in the past ten years, indicating that returning farmland to forest and grassland in this area has achieved remarkable results. The area of LAI degradation is mainly in the area with a more immense LAI value in the previous period, and these areas are similar to FVC and NPP degradation areas, but the degradation area is far less than that of FVC and NPP, about 1.05 $\times {10}^4{ \mathrm{km}}^2$.

3.3. Soil Conservation Services

In 2005, 2010, and 2015, the total amount of soil conservation (ASC) in the LPES was 2.35 $\times {10}^8 \mathrm{t}$, 2.39 $\times {10}^8 \mathrm{t},$ and 2.53 $\times {10}^8 \mathrm{t}$, respectively. In the past ten years, the ASC has increased by 0.18 $\times {10}^8 \mathrm{t}$. Nevertheless, from the classification results of ASC (Figure 14 and Figure 15), from 2005 to 2015, low-level soil conservation areas were the central area, accounting for about 87%, mainly distributed in Suide County and Zizhou County in Shaanxi, as well as Xiaoyi City and Fenxi County in Shanxi. The high-level soil conservation area is only 0.15$\times {10}^4 {\mathrm{km}}^2$, mainly concentrated in Huangling County, Huanglong county, and Yichuan County.

Among all ecosystem types (

Table 8. Soil conservation services of various ecosystem types.

Ecosystem	2005		2010		2015
	$\mathbf{ASC}/{10}^8\mathbf{t}$	Ratio/%	$\mathbf{ASC}/{10}^8\mathbf{t}$	Ratio/%	$\mathbf{ASC}/{10}^8\mathbf{t}$	Ratio/%
Cropland	0.58	24.70	0.60	25.16	0.63	24.73
Forest	0.89	38.06	0.89	37.22	0.96	37.98
Grassland	0.85	36.40	0.88	36.78	0.92	36.26
Water body	0.02	0.76	0.00	0.00	0.00	0.02
Built-up land	0.00	0.00	0.02	0.76	0.02	0.91
Others	0.00	0.08	0.00	0.08	0.00	0.11

), the forest has the highest soil conservation service, with an average of (3414 t/${\mathrm{km}}^2)$. The ASC per unit area of grassland (1935.38 t/${\mathrm{km}}^2$) and cultivated land (1476.88 t/${\mathrm{km}}^2$) were also significantly higher than other ecosystem types. The ASC of water body (1317.92 t/${\mathrm{km}}^2)$ and Built-up land (1220.09 t/${\mathrm{km}}^2)$ were close. In terms of annual growth rate, the growth rate of soil conservation services of the forest is remarkably rapid, with the rate of 22.92 t/$\left({\mathrm{km}}^2·a\right)$, followed by that of grassland with the rate of 13.60 t/$\left({\mathrm{km}}^2·a\right)$. In the LPES, forests and grasslands have played a critical role in soil conservation. The ASC of these two ecosystems accounts for more than 74%. The study also found that cropland has a soil conservation service that cannot be ignored. In the three periods, the ASC of farmland was about 25% of the total.

From 2005 to 2015, the soil conservation services (SCS) in the LPES increased significantly (Figure 16). The SCS in the area of about 5.45 $\times {10}^4 {\mathrm{km}}^2$ showed an increasing trend. From 2005 to 2015, the ASC per unit area increased by 150.81 t/$\left({\mathrm{km}}^2·a\right)$. The areas with a significant increase in ASC are mainly distributed in Xiaoyi City, Fenxi County, Shanxi Province. The ASC per unit area in these areas increases from 1170.64 t/$\left({\mathrm{km}}^2·a\right)$ to 1503.02 t/$\left({\mathrm{km}}^2·a\right)$. About 4.48 $\times {10}^4 {\mathrm{km}}^2$, accounting for 37.25% of the barrier area, the soil conservation function is unchanged. The area of soil conservation decreased is only 3/8 of the increased area, about 2.09 $\times {10}^4{ \mathrm{km}}^2$, which is mainly distributed in the west of Luliang Mountain to the east of Huanglong Mountain, and the west of Ziwuling mountain. In the past ten years, although there is still local degradation in the LPES, the overall soil erosion situation has improved, and ecological protection policies such as reforesting the cultivated land have played a positive role.

3.4. Analysis of Driving Forces of Soil Conservation Services Changes

Based on the research and analysis on the pattern, quality, and SCS of the ecosystem in this paper, it can be seen that: the assessment based on landscape pattern can only reflect the quantity and structural changes of each ecosystem, while the assessment based on quality (FVC, NPP, and LAI) and SCS can provide more information about the actual situation of the ecosystem. Additionally, there are many factors involved in SCS assessment, which can better characterize the state of the ecosystem. Concomitantly, LPSE is the main soil erosion control area in China. Studying the driving mechanism of SCS change can provide more direct ecological information for decision-making. As natural and social stresses both affect ES [35]. To analyze the spatial pattern of SCS and the driving mechanism of their spatial changes in the LPES, we used Pearson correlation analysis and redundancy analysis to explore the effects of natural stress and social activity stress factors on the spatial distribution SCS. Natural stress is mainly affected by climate change. In this study, total solar radiation (SOL), precipitation (P), and temperature (T) are selected to represent the impact of climate change. Social stress is studied by selecting the primary industry (PI) and built-up land (BL). As can be seen from

Table 9. The linear regression correlation coefficient of soil conservation services in the Loess Plateau Ecological Screen.

	FVC	P	T	SOL	PI	BL
Standardized Coefficient	0.18	0.88 **	−0.06	0.61 **	−0.18	−0.43 **

**. p < 0.01

, from 2005 to 2015, SCS is positively correlated with climate change (P, SOL) and vegetation cover (FVC), and there is a significant positive correlation with P and SOL (R = 0.88, p < 0.001; R = 0.61, p < 0.001). The average annual precipitation contribution rate is the largest at 88%. It is negatively correlated with economic development (PI), BL, and T, and BL is significantly negatively correlated (R = 0.56, p < 0.01).

We use redundancy analysis to further study the contribution rate of different influencing factors to changes in ES. The results are shown in Figure 17. Climate factors are the main factors affecting the change of ES in the LPES. Precipitation, total solar radiation (SOL), and temperature can explain 52%, 30.1%, and 17% of the change trends of SCS, respectively. Construction land and primary industry negatively correlated with SCS, explaining 22% and 8% of the change trends, respectively. Therefore, in climate change, optimizing the industrial structure and rational urbanization is of great significance to the protection of ES in the LPES.

4. Discussion

In this paper, multi-source data were used to study the ecosystem of LPES by various methods. This study can help us to understand its dynamic changes and provide suggestions for ecosystem management in this region. Our research results showed that grassland and farmland were the main ecosystems in LPES. From 2005 to 2015, the ecosystem’s structure was relatively stable, and FVC, LAI, and NPP were significantly improved. The average FVC of LPES was about 56 % and is growing at an average annual rate of 0.91 %. The high annual FVC value was the maximum value synthesis method, and the MODIS NDVI data, including the growing season (June–September), was selected in the calculation. The annual maximum value of NDVI can reflect the best growth condition of vegetation in a year and further eliminate the influence of factors such as the atmosphere and solar altitude. This method has been widely used in the study of the ecosystem in large regions. The evaluation results showed that NPP was generally improved with an average annual growth rate of 6.94 gC/(m²), consistent with the study of Wu et al. [53]. However, the promotion effects of NPP in our study was more significant than those of Wu et al., which may be related to the different scope of the study area. As is a part of the broad Loess Plateau, the hydrothermal conditions in our study area are better than other regions of the Loess Plateau. Our results showed that LAI in LPES was significantly improved which was consistent with the research results of Kou et al., Wang et al. [54,55]. LPSE is a vital soil conservation district in China, so we also evaluated the changes of primary ES (SCS) in the study area. ASC continued to increase by 0.18 $\times {10}^8$t. From 2005 to 2015, and we found that forests and grasslands played a crucial role in soil conservation. Due to the apparent improvement of SCS in the region, it is particularly vital to explore the driving factors behind it. The annual average precipitation contributed largest to the variance of SCS in natural factors, reaching 88%. The results of this study are consistent with the conclusions of Zhang Siyu (2019) and Liu Xiaona (2018) [56,57]. Correspondingly, increased precipitation may lead to increased ASC. However, Zhang Jinxi claimed that rainfall will make soil erosion more serious. Additionally, according to Formula (21), rainfall will also increase the potential erosion. The implementation of ecological engineering, such as returning farmland to forest, could improves the vegetation ecosystem, increasing the difference between potential soil erosion and actual soil erosion. That is, the amount of soil conservation will increase. The pressure of human society, especially urban expansion, negatively affected the amount of SCS. Although farmland has a higher amount of SCS in this study, which is negatively correlated with SCS. The expansion of farmland will worsen soil erosion, as the expansion of farmland is bound to lead to shrinking grasslands and forests, and the ASC of these ecosystems is higher.

This study also showed that degradation still exists in some areas of the LPES. In the study area, land reclamation and construction land occupation of farmland continue to exist simultaneously, showing an increasing trend. For example, from 2010 to 2015, about 284.45 ${\mathrm{km}}^2$ of farmland was converted to construction land. In 2010–2015, about 357.73 ${\mathrm{km}}^2$ of grassland and forest were converted to cultivated land, about five times that of the first five years. Degradation can be observed through the research results of NPP, LAI, FVC, and SCS, and the results may be more reliable because these results are simulated through complex process models, which describe the inner state of the ecosystem in detail [58,59,60,61]. In the degradation assessment, the area of vegetation cover reduction was the largest, reaching 8.20 $\times {10}^4 {\mathrm{km}}^2$, and the area of ASC reduction was about 2.09 $\times {10}^4 {\mathrm{km}}^2$. However, the area of LAI degradation was smaller than that of FVC, about 1.05 $\times {10}^4 {\mathrm{km}}^2$. The above-degraded areas are mainly located in the west of Luliang Mountain, the east of Huanglong Mountain, and the west of Ziwuling mountain. One of the reasons is that many residential areas in these areas are greatly affected by interference of human activities. Therefore, in the future, it is still necessary to reduce unreasonable urban expansion behavior and continue to implement ecological projects such as returning farmland to forests to promote the continuous improvement of the ecological environment.

This paper also has some shortcomings and the future need for further in-depth research. This research analyzed the first-level ecosystem pattern of the LPES. Although this can reflect the changes in the ecosystem from 2005 to 2015 to a certain extent, a detailed exploration was lacking. In the next step, we will analyze the pattern changes of secondary and even tertiary ecosystems, study the evolution of the ecosystem pattern of LPES more precisely, and reveal the quantitative relationship between the horizontal pattern, vertical structure, quality, and function of the ecosystem. Additionally, we will determine the mutation point and its spatial range of ecosystem change, carry out ecosystem simulations in different scenarios, and put forward countermeasures and suggestions. In addition, among the multiple ES, only the SCS of the LPES in this study was evaluated, and the synergy and trade-off between multiple ES will be reported in the following work. In research, we also found that the differences in the degradation areas was determined by each index when using NPP, LAI, FVC, and SCS to evaluate the ecosystem. Using one indicator alone can lead to wrong conclusions while multiple indicators can lead to complex conclusions. We will report in detail on this issue in another study and build a composite index to assess the state of ecosystems. Next, we will explore the occurrence regularity, variation characteristics, and driving factors of these synergy or trade-off relationships over space and time scale and determine the conditions, scope, and threshold of the occurrence of the synergy mechanism of the main ES items. In the sustainable management of the LPES, it is necessary to combine regional development goals and take into account other ES to maximize ecological benefits, to help avoid the decision of sacrificing other ES in pursuit of specific ES.

5. Conclusions

This paper uses multi-source data and multiple indicators to evaluate the changes in the LPES ecosystem from 2005 to 2015. The main conclusions are as follows. (1) Farmland and grassland are the main ecosystems. During the past 10 years, the cultivated land area continued to decrease, and the forest and grassland area continued to expand. However, the most rapid expansion during this period was the construction of land. (2) The eco-environment showed an overall improvement, FVC, LAI, NPP, and SCS increased significantly in the past 10 years, but local degradation still occurred. (3) As an important ecosystem service in the study area, SCS was positively correlated with climate change (P, SOL) and FVC, and negatively correlated with PI, BL, and T. Despite the results highlighted the vital ecosystem service of farmland, the increased risk of soil erosion triggered by deforestation was observed. Therefore, in the context of climate change, the study area is expected to continue to implement the transition of farmland to forests, further curb the behavior of deforestation and reclamation, and also prevent the blind expansion of cities, which is very necessary to promote the sustainable improvement of the eco-environment of LPES.