Pattern and Trend of Ecosystem Service Value in the Loess Plateau of Northern Shaanxi

Abstract

The ecosystem service value (ESV) is a critical metric for assessing the construction and protection of the environment. The research into the ESV pattern and the future development trend in the Loess Plateau of Northern Shaanxi is important for the conservation of water and soil and the enhancement of the natural environment in the region. In this study, the variations and distribution patterns of the ESV in the study area from 2000 to 2020 were analyzed, the influence of various natural and social factors on the ESV was quantified, the weight of each factor was analyzed and evaluated using the entropy weighting method, and, finally, a prediction was made regarding how the ESV will develop going forward in this area. The results show that (1) the ESV showed a decreasing trend from 2000 to 2020, with the highest value for soil conservation and the lowest value for food production. Among the 25 districts and counties, Suide County had the lowest ESV per unit area, whereas Huanglong County had the highest. (2) The global positive correlation was clearly visible in the ESV. According to local spatial autocorrelation analysis, the area had a “high-high” agglomeration area in the south and a “low-low” agglomeration area in the middle and north. (3) Among the various influencing factors, population density had the highest weight and the distance from roads had the lowest weight. The impact status of the area generally showed a lighter impact in the southern region and a heavier impact in the northern region. (4) In 2030, the total ESV is predicted to be CNY 4343.6 million in the study area, CNY 39 million lower than that in 2020.

1. Introduction

Ecosystem services are the welfare and advantages that humans derive from ecosystems, such as supply, regulation, support, and cultural services [1,2,3], which directly or indirectly support human survival and development [4]. With the rapid socioeconomic growth worldwide, the impact of human activities on the ecological environment has been intensified, and the ecosystem has been repeatedly damaged [5]. Ecosystem service value (ESV) is a quantitative evaluation of ecosystem service function, and the study of ESV can measure the ecological condition of a region, which is important for the conservation, restoration and sustainable development of regional ecosystems.

Costanza et al. (1997) conducted the first study on estimating the value of global ecosystem services and natural capital [1]. It provided theoretical and methodological support for the measurement of the ESV. The current research methods for ESV calculation, which include the market-price method, the shadow-engineering method, the hedonic-price method, the hydrological- and water-quality-modeling method, and the willingness-to-pay method, have been applied to ecological resources, planning, construction and management, government policies and so on [6,7,8,9,10,11]. They are more accurate in assessing a single ecosystem service but lack comparability and applicability. Additionally, the equivalent factor method is more prominent in reflecting the change pattern and mechanism of action of ESV [12,13,14]. It has a wide range of applications and is highly applicable. The Chinese scholars Xie et al., based on the method of calculating of the unit-area-value-equivalent factor, obtained relevant expert knowledge based on the model proposed by Costanza et al. and modified the global ESV-equivalent factor by combining it with the actual ecological and environmental conditions in China to make it more consistent with the actual conditions in China [15,16].

The calculation of the ESV measures the quality of ecosystem services in a region. It is a monetized representation of ecosystem services. It can help governments to understand the general state of ecosystems [17,18]. However, it does not reflect the drivers that influence ecosystem service functions. It leads to the inability to analyze the specific causes of ecosystem destruction at the root cause. Therefore, it is also crucial to study the factors that influence the ESV. The ESV is directly or indirectly influenced by factors both internal and external to the ecosystem. Among them, natural factors have direct impacts on the ESV within the ecosystem [19]. Studies have shown that the elevation, slope, and direction of the terrain affect ecological service functions such as soil and water conservation, raw material production, and water-supply capacity [20]. Climate influences plant growth and enables changes in ecosystem services through the regulation of surface water and heat conditions [21,22]. The growth activities of organisms and the spatial characteristics of habitats influence ecosystem service functions such as food production and the maintenance of biodiversity [23,24]. Additionally, as human activities intensify, the impact of anthropogenic factors on the ESV become more prominent [25,26]. Population density, socioeconomic activities, and malicious destruction of the environment all have important effects on the ESV [27,28].

The study region’s natural surroundings in the Loess Plateau, with low annual precipitation, low vegetation cover, steep terrain, loose soils, and the constant disturbance of human activities, have caused a highly vulnerable biological environment to develop, with land desertification, soil erosion, and salinization. Zhang et al. analyzed the extent to which land-use change has an impact on the ESV as well as the coordination between ecological and environmental quality and socioeconomic development [29]. Jiang et al. used a mix of benefit transfer and proxies for land-use and land cover to study the ESV resulting from changes in the land-use and land cover on the Loess Plateau between 1990 and 2015 [30]. Fang et al. examined and evaluated the connection between the ESV and land-use patterns, as well as assessing the association between land-use change and the ESV at the town scale [31]. Therefore, in the Loess Plateau, studies on ecosystem services have mainly focused on analyzing the effects of land-use/cover changes on the ESV within the Loess Plateau region, while fewer studies have been conducted to examine the mechanisms, quantify the effects of various environmental and societal factors on the ESV, and forecast future trends in the field. In this study, based on the investigation of the ESV distribution and change in the Loess Plateau of Northern Shaanxi from 2000 to 2020, the influence status of multiple factors on the ESV was focused on and the future development trend was predicted. The objective of this study is to address the following questions: (1) What are the ESV’s temporal trends and regional distributions in the Loess Plateau of Northern Shaanxi from 2000 to 2020? (2) What is the degree of influence of various factors on the ESV in the Loess Plateau of Northern Shaanxi? (3) What is the status of the integrated constraints on ecosystem development within the Loess Plateau of Northern Shaanxi? (4) What is the future development trend of the ecological environment of the Loess Plateau in Northern Shaanxi?

2. Materials and Methods

2.1. Study Area

The study area is in the heart of the Loess Plateau. It consists of two prefecture-level cities, Yan’an and Yulin (Figure 1). The area is a paleotopography formed by Mesozoic bedrock, covered with a thick layer of loess and red soil from the Cenozoic, and then formed by cutting with running water. It is in the transition zone from semiarid to arid, the annual average temperature is approximately 8–14 °C, and the annual rainfall is approximately 350–660 mm. The spatial distribution reveals that the northwest receives less precipitation than the southeast due to an irregular seasonal distribution. The precipitation is mainly concentrated in summer and the precipitation intensity is high. It often appears in the form of heavy rain.

2.2. Data

The average temperature and annual precipitation data in 2000, 2010, and 2020 used were derived from the “Shaanxi Provincial Statistical Yearbook” and “China Environmental Bulletin” of each year. The land-use data for 2000, 2010, and 2020 were from the GlobeLand30 platform (http://globeland30.org/, accessed on 18 January 2022) and were divided into seven types: cropland, forestland, grassland, wetland, water bodies, construction land, and unused land. The slope and elevation data were from the Chinese Academy of Sciences Geospatial Data Cloud Platform (http://www.gscloud.cn/, accessed on 10 February 2022). The annual-rainfall, annual average temperature, population density, and soil-type data were from the Resource and Environmental Science and Data Center of Chinese Academy of Sciences (http://www.resdc.cn/, accessed on 10 February 2022). The main road data came from the BIGEMAP platform (http://www.bigemap.com/, accessed on 10 February 2022).

2.3. Methods

2.3.1. Miami Model

Xie et al. [32] proposed the ESV-equivalent factor per unit area of the Chinese ecosystem, which has been widely applied by many scholars in many aspects [33,34,35]. However, since it is based on the national average level, the ESV needs be adjusted to account for regional variations in a small area when it is calculated. The climatic productivity of an area is a key factor affecting its ecological environment and the ESV [36]. Therefore, this study calculated the climatic productivity of the area as well as the national average climatic productivity to obtain the regional difference coefficient and then revised the equivalent factor [37]. The ESV was calculated using the revised equivalent factor in combination with the land-use status of the study area.

Using the Miami model to calculate climate productivity, the specific formula is as follows:

$$\mathrm{np}{p}_R=30000\times \left(1-e^{-0.000664R}\right)$$

(1)

$${\mathrm{npp}}_{\mathrm{T}}=\frac{30000}{1+{\mathrm{e}}^{1.315}-0.119\mathrm{T}}$$

(2)

$${\mathrm{npp}}_{\mathrm{p}}=\mathrm{MIN}\left\{{\mathrm{npp}}_{\mathrm{T}},{\mathrm{npp}}_{\mathrm{R}}\right\}$$

(3)

where $R$ is the annual precipitation, $T$ is the annual average temperature, $np{p}_R$ is the net primary productivity of vegetation calculated using the annual precipitation, and $np{p}_T$ is the net primary productivity of vegetation calculated using the annual average temperature. According to Liebig’s law of least factor, $np{p}_p$ takes the minimum of $np{p}_R$ and $np{p}_T$ as the regional climate productivity.

$${\mathrm{S}}_{\mathrm{k}}=\frac{{\mathrm{npp}}_{\mathrm{k}}}{\mathrm{npp}}$$

(4)

where $npp$ is the national average climate productivity, $np{p}_k$ is the climate productivity of the study area, and $S_k$ is the regional difference coefficient.

The following is the formula for calculating the economic value of an equivalent factor in the study area:

$$E_a={E^{\prime }}_a\times S_k$$

(5)

where ${E^{\prime }}_a$ is the value corresponding to an equivalent factor nationwide (CNY 449.1/hm²).

2.3.2. Spatial-Autocorrelation-Analysis Method

The possible dependency of the observation data of variables in the same distribution region is referred to as spatial autocorrelation [38], which is one of the effective methods used to analyze the spatial correlation of the same variable between different observation objects [39,40]. In this study, spatial autocorrelation analysis was used to explore the spatial distribution pattern of the ESV. The global spatial autocorrelation and local spatial autocorrelation were used to investigate the overall distribution characteristics of the ESV in the study area and the correlation degree of the ESV between each district and neighboring districts and counties, respectively.

Global spatial autocorrelation:

Global spatial autocorrelation can reflect the general trend of spatial correlation within the entire study area [41]. The specific formula is as follows:

$$\mathrm{I}=\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{\mathrm{W}}_{\mathfrak{i}\mathrm{j}}·\left({\mathrm{Y}}_{\mathrm{i}}-\overline{\mathrm{Y}}\right)·\left({\mathrm{Y}}_{\mathrm{j}}-\overline{\mathrm{Y}}\right)}{{\mathrm{S}}^2·{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{\mathrm{W}}_{\mathrm{ij}}}$$

(6)

where I is the global spatial autocorrelation index, $S^2=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left(Y_i-\overline{Y}\right)$, $Y_i$ is the ESV of the county (district), $n$ is the number of counties (districts), and $W_{ij}$ is the spatial weight matrix.

Local spatial autocorrelation:

While global spatial autocorrelation cannot highlight the aggregation state of local area units, local spatial autocorrelation can accurately capture the aggregation and divergence characteristics of spatial units and neighboring units [42].

$${\mathrm{I}}_{\mathrm{i}}=\frac{\left({\mathrm{Y}}_{\mathrm{i}}-\overline{\mathrm{Y}}\right)}{{\mathrm{S}}^2}{{\displaystyle \sum }}_{\mathrm{j}}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}·\left({\mathrm{X}}_{\mathrm{i}}-\overline{\mathrm{X}}\right)$$

(7)

where $I_i$ represents the local spatial autocorrelation Moran index, and the other variables have the same meaning as in Equation (6).

2.3.3. Entropy Weight Method

The entropy weight method is an objective assignment method [43]. According to the explanation of the basic principles of information theory, entropy and information are measurements of the degree of disorder and order of a system, respectively [44]. If the index’s information entropy is lower, its weight should be higher and its role in the comprehensive evaluation should be larger because it provides more information [45]. The entropy weight method overcomes the influence of human subjective consciousness and makes the evaluation results more objective [46]. In this study, based on the quantitative analysis of the influence of each factor on the ESV, the entropy weight method was used to calculate the weights of each influencing factor, analyze the specific reasons limiting the development of the ESV, and used in the subsequent comprehensive impact zoning study. Its specific calculation process is divided into the following three steps:

Data standardization:

$$\{\begin{array}{c}{\mathrm{P}}_{\mathrm{ij}}=\frac{{\mathrm{x}}_{\mathrm{ij}}-{\mathrm{x}}_{\mathrm{ijmin}}}{{\mathrm{x}}_{\mathrm{ijmax}}-{\mathrm{x}}_{\mathrm{ijmin}}} \mathrm{Positive} \mathrm{indicator}\\ \\ {\mathrm{P}}_{\mathrm{ij}}=\frac{{\mathrm{x}}_{\mathrm{ijmax}}-{\mathrm{x}}_{\mathrm{ij}}}{{\mathrm{x}}_{\mathrm{ijmax}}-{\mathrm{x}}_{\mathrm{ijmin}}} \mathrm{Negative} \mathrm{indicator}\end{array}$$

(8)

Information entropy calculation:

$$E_j=-ln{\left(n\right)}^{-1}{{\displaystyle \sum }}_{i=1}^n{p}_{ij}ln{p}_{ij}$$

(9)

where $p_{ij}=Y_{ij}/{{\displaystyle \sum }}_{i=1}^n{Y}_{ij}$; if $p_{ij}$ = 0, then define $\underset{p_{ij}\to 0}{\lim }{p}_{ij}ln{p}_{ij}=0$.

Calculation of the weight of each indicator through information entropy:

$$W_i=\frac{1-E_i}{k-{{\displaystyle \sum }}^{}{E}_i} (\mathrm{i}=1,2,\dots ,\mathrm{k})$$

(10)

2.3.4. CA–Markov Model

The CA–Markov model combines the spatial dynamic evolution capability of cellular automata and the time dynamic advantage of the Markov model [47,48,49,50]. This model has good predictive ability in landscape pattern research evolution and land-use dynamic change [51].

Cellular automata are defined as cells with states that are discrete and finite in a cellular space. This represents a dynamic system that follows certain regional laws as it develops in a discrete time dimension. The most basic components of cellular automata include the cellular unit, cellular space, domain, and rules [52,53]. It can usually be expressed as follows:

$$S_{t+1}=f\left(S_t,N\right)$$

(11)

where $S$ is a set of discrete and finite states of the cellular, $t$ and $t+1$ represent time, $f$ is the transformation rule of the local cellular space, and $N$ is the neighborhood of the cellular unit. At present, cellular automata are widely applied in population distribution change, land-use change, forest fires, and other aspects [54].

The Markov model is a method of stochastic process theory formation proposed by Russian mathematician Markov. It calculates the probability of things in the initial period of different states and the transition probability between states and then shows the trend of transformation between different states to realize the prediction of the future state of things [55,56]. In the process of Markov model prediction, the state of things in a past moment has nothing to do with the development of the state. So, the transfer process has no aftereffect [57], which is called the Markov process.

$$\mathrm{P}={\mathrm{P}}_{\mathrm{ij}}=\left[\begin{array}{ccc}{\mathrm{P}}_{11}& \cdots & {\mathrm{P}}_{1\mathrm{n}}\\ ⋮& \ddots & ⋮\\ {\mathrm{P}}_{\mathrm{n}1}& \cdots & {\mathrm{P}}_{\mathrm{nn}}\end{array}\right]$$

(12)

where $P_{ij}$ is the transition probability matrix of the land-use change state, $0\le P_{ij}<1$, and ${{\displaystyle \sum }}_{j=1}^n{\mathrm{P}}_{\mathrm{ij}}=1$, $i, j=1, 2 \dots$ is the land type; in this paper, $n=7$.

In this study, the trend forecast of the ESV for 2030 in the study area was achieved by means of the CA–Markov model. The data of 2000 and 2010 were used to forecast the study area in 2020 before conducting the analysis, and the accuracy was compared with the actual condition in 2020, and the accuracy was verified to be qualified before conducting the subsequent forecast. Additionally, the forecast was made by adding constraints, such as elevation, slope, water body, distance from a road, and distance from a water body, to create an adaptation atlas to ensure the forecast results are more realistic.

3. Results

3.1. Land-Use Status Analysis

Through analysis, the land-use map of the study area from 2000 to 2020 was obtained (Figure 2). It was found that the area of grassland was the largest in the area, with an area of more than 45% in the three phases, followed by cropland, accounting for more than 32%. Wetland accounted for the smallest area, with a ratio of 0.18% (2000)–0.05% (2020). The sum of grassland, cropland, and forestland accounted for 96.19% (2020)–97.24% (2000) of the total size of the study area. The amount of construction land exhibited an upward trend, from 0.36% in 2000 to 1.68% in 2020, and the growth rate accelerated year-by-year. The area of wetlands in the study area decreased by 64.38% from 2000 to 2010, followed by water bodies, which decreased by 36.90%, and construction land grew the most, with an increase of 55.90%. From 2010 to 2020, the wetland in the study area decreased by 30.77% and the construction land increased by 199.55% (Figure 3).

3.2. ESV Analysis

3.2.1. Revision of the ESV

According to the annual average temperature and annual precipitation data of the area and the whole country (

Table 1. Annual average temperature and annual precipitation in Yan’an, Yulin, and nationwide.

Yan’an City		Yulin City		Nationwide
Annual Average Temperature	Annual Rainfall	Annual Average Temperature	Annual Rainfall	Annual Average Temperature	Annual Rainfall
/°C	/mm	/°C	/mm	/°C	/mm
11.0	367.3	9.5	264.9	9.1	633.2
11.1	465.8	9.3	363.9	9.5	681.0
10.3	695.2	10.0	526.4	10.3	645.5

), the climatic productivity of the study area was calculated to be 770.81, the national average climatic productivity was 838.27, the regional coefficient was 0.92, and the equivalent factor’s corresponding value in the area was CNY 412.96/hm² (

Table 2. The ecosystem service value (ESV) per unit area of different land types.

Primary Type	Secondary Type	Forestland	Grassland	Cropland	Wetland	Water Bodies	Barren
Supply services	Food production	136.28	177.57	412.96	148.67	218.87	8.26
	Raw-material production	1230.62	148.67	161.05	99.11	144.54	16.52
Regulating services	Gas regulation	1783.98	619.44	297.33	995.23	210.61	24.78
	Climate regulation	1680.74	644.22	400.57	5595.59	850.69	53.68
	Hydrological regulation	1689.00	627.70	317.98	5550.16	7751.23	28.91
	Waste disposal	710.29	545.11	574.01	5946.60	6132.43	107.37
Support services	Soil conservation	1660.09	925.03	607.05	821.79	169.31	70.20
	Maintain biodiversity	1862.44	772.23	421.22	1523.82	1416.45	165.18
Cultural services	Landscape provision	858.95	359.27	70.20	1936.78	1833.54	99.11
	Total	11,612.39	4819.23	3262.37	22,617.74	18,727.67	574.01

).

3.2.2. Changes in the ESV

Table 3. ESV of the Loess Plateau in Northern Shaanxi. Unit: CNY 100 million.

	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
2000	83.76	169.13	180.60	3.30	6.29	0.83	443.90
2010	83.82	167.04	181.38	1.18	3.97	0.86	438.24
2020	85.05	168.88	174.77	0.81	4.42	0.82	434.75

shows that the ESV of the study area in 2000, 2010, and 2020 was CNY 44.390 billion, 43.824 billion, and 43.475 billion, respectively, with an overall decrease of CNY 915 million, or 2.06%. From 2000 to 2010, the total value decreased by CNY 566 million, or 1.28%. From 2010 to 2020, the total value decreased by CNY 349 million, or 0.80%, indicating a slower decrease rate.

In terms of individual types, the total ESV of grassland and forestland was relatively high, while the value of unused land and wetland was low. Although wetland had the highest value per unit area, their area was smaller, so they provided lower values. The unused land refers to the natural covered land with vegetation coverage less than 10%, which had a low ecological service function.

3.2.3. Changes in the Individual ESV

Table 4. The individual ESV in the Loess Plateau of Northern Shaanxi. Unit: CNY 100 million.

	Food Production	Raw Material Production	Gas Regulation	Climate Regulation	Hydrological Regulation	Waste Disposal	Soil Conservation	Maintain Biodiversity	Landscape Provision
2000	19.35	27.72	57.08	60.08	59.74	48.59	74.71	67.81	28.82
2010	19.32	27.50	56.75	59.27	58.06	47.25	74.47	67.30	28.32
2020	19.26	27.55	56.28	58.73	57.69	46.88	73.68	66.70	27.99

shows that the single ESVs of food production, gas regulation, climate regulation, hydrological regulation, waste disposal, soil conservation, maintenance of biodiversity, and landscape provision in the study area from 2000 to 2020 all decreased. Among them, the hydrological regulation value decreased the most, by CNY 205 million, and the food production value decreased the least, by CNY 9.25 million. The value of soil conservation over the years was the highest and food production value was the lowest.

3.2.4. Changes in the ESV at the District and County Scales

The land-use status within the 25 districts and counties in the study area was obtained separately by district and county, and then the ESV per unit area within each district and county was calculated for 2000, 2010, and 2020. According to Figure 4, the ESV distribution in the study area was typically higher in the south and lower in the north. From 2000 to 2020, the ESV per unit area in Huanglong County was the highest among the 25 counties, with CNY 10,589/hm², CNY 10,604/hm², and CNY 10,429/hm², while the ESV per unit area in Suide County was only CNY 3656/hm², CNY 3637/hm², and CNY 3607/hm².

By analyzing the changes in the ESV per unit area of each district and county (Figure 5), it was found that, from 2000 to 2010, the ESV per unit area decreased in 19 districts and counties and increased in 6 districts. From 2010 to 2020, the number of counties with a decreasing ESV per unit area was 16 and the number of counties with increasing values was 9. Overall, the unit area value of most of the districts and counties showed a decreasing trend from 2000 to 2020, and the ecological condition of the area was deteriorating.

3.3. The Spatial Distribution Pattern of the ESV

3.3.1. Global Spatial Autocorrelation Analysis

Scatter plots of Moran’s I index for 2000, 2010, and 2020 were obtained by means of global spatial autocorrelation analysis (Figure 6). Figure 6 shows that the global Moran’s I indexes for the ESV per unit area in 2000, 2010, and 2020 for 25 districts and counties in the area were 0.666, 0.661, and 0.670, respectively, which were all greater than 0. The results show that the unit ESV of the 25 districts and counties showed a strong positive correlation in each year and the spatial distribution showed obvious clustering characteristics.

3.3.2. Local Spatial Autocorrelation Analysis

To determine whether there was a local agglomeration effect among the 25 districts and counties in the area, the local Moran’s I (LISA) value was calculated through local spatial correlation analysis and a LISA distribution map was created (Figure 7).

Figure 7 shows that local spatial autocorrelation pattern of the ESV per unit area of the 25 districts and counties showed a “high-high” agglomeration area and a “low-low” agglomeration area. There was no “high-low” agglomeration area or “low-high” agglomeration area. The “high-high” agglomeration area was distributed in the south of the study area, namely, Fuxian County, Luochuan County, Yichuan County, and so on. The unit ESV of these areas had a strong positive spatial correlation mechanism and improved the unit-area ESV of the surrounding area. The “low-low” agglomeration areas were primarily distributed in the central and northern areas of the study area, namely, Hengshan County, Jia County, Suide County, Zizhou County, etc., indicating that the ESV per unit area of these areas and the ESV per unit area of the neighboring areas had a positive correlation with each other at a low level.

3.4. Influencing Factors of the ESV

The spatial distribution maps of different influencing factors were obtained, and different levels were assigned according to the magnitude of the values, which are overlaid with the land-use maps to calculate the ESV of the corresponding level according to the land-use status within the different levels of each factor, so as to investigate the influence of each factor on the ESV status.

3.4.1. Impact of Annual Rainfall on the ESV

The annual rainfall data were reclassified into four levels: below 400 mm, 400–500 mm, 500–600 mm, and above 600 mm.

Table 5. ESV per unit area of different rainfall levels. Unit: CNY/hm².

Rainfall (mm)	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
<400	1113.06	14.86	2758.07	46.31	73.10	40.15	4045.55
400–500	1198.39	342.52	2829.68	9.98	53.50	5.35	4439.42
500–600	843.56	6009.19	1035.75	0.00	33.80	0.00	7922.31
>600	352.67	9896.70	153.45	0.00	0.00	0.00	10,402.83

shows that, with the increase in rainfall, the ESV per year with an annual rainfall of more than 600 mm was 2.57 times that of the area with an annual rainfall of less than 400 mm, indicating that rainfall had a great impact on the ESV in this area. According to comparison of land-use status, forestland accounted for 85.23% in areas with an annual rainfall above 600 mm and only 0.13% in areas with rainfall below 400 mm, indicating that vegetation grew vigorously in areas with abundant precipitation and had a high ESV.

3.4.2. Impact of Temperature on the ESV

The annual average temperature data were reclassified into four levels: below 8 °C, 8–9 °C, 9–10 °C, and above 10 °C.

Table 6. ESV per unit area of different temperature levels. Unit: CNY/hm².

Annual Average Temperature (℃)	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
<8	1557.62	126.84	2448.67	7.97	39.59	0.61	4181.30
8–9	793.24	892.25	3074.54	11.27	56.00	19.97	4847.27
9–10	866.94	2394.44	2394.17	19.11	31.65	14.11	5720.43
>10	1400.66	2967.57	1449.95	12.75	64.04	0.54	5895.51

shows that, with the increase in the annual average temperature, the ESV per unit area showed an increasing trend, with a total increase of 41.00%. According to the analysis of land-use status, cropland, forestland, and grassland accounted for 98.58% of the area above 10 °C and unused land accounted for the lowest proportion, among which 47.28% of forestland and 43.12% of water bodies were distributed in this area. This result indicated that the area with higher hydrothermal conditions had a strong ecological function and a high ESV.

3.4.3. Impact of Slope on the ESV

The slope data were reclassified into five levels: 0–2°, 2–6°, 6–15°, 15–25°, and above 25°.

Table 7. ESV per unit area of different slope levels. Unit: CNY/hm².

Slope(°)	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
0–2	986.39	294.08	2698.47	29.41	77.69	38.74	4124.79
2–6	1249.29	1182.98	2302.67	3.14	70.13	8.24	4816.45
6–15	1143.62	2614.30	1973.11	4.54	39.49	1.79	5776.86
15–25	882.70	3599.75	1973.09	2.46	18.36	0.75	6477.10
>25	660.94	2897.30	2577.95	22.57	18.69	2.86	6180.32

shows that the ESV per unit area in the 0–25° region showed an increasing trend, increasing by 57.03%. When the slope was greater than 25°, the ESV decreased, but only by 4.58%. In the 0–2° area, forestland accounted for only 2.53% and contained 47.47% of the construction land and 76.81% of the unused land. In the area above 15°, grassland and forestland accounted for 72.28%, playing a strong ecological function. Therefore, slope was positively correlated with the ESV.

3.4.4. Impact of Soil Types on the ESV

According to the traditional classification system of soil occurrence, the soil types were divided into 10 types, including semileached soil, calcium-layer soil, arid soil, primordial soil, etc.

Table 8. ESV per unit area of different soil types. Unit: CNY/hm².

Soil Type	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
Semileached soil	174.53	10,731.31	99.42	7.91	6.55	0.00	11,019.71
Calcium-layer soil	1643.42	969.95	1903.01	11.70	29.06	1.04	4558.18
Arid soil	1360.28	0.00	2656.53	79.92	0.00	10.14	4106.88
Primordial soil	1043.84	1907.76	2354.29	37.53	68.26	10.92	5422.60
Semihydrated soil	1054.50	68.42	2851.78	152.31	267.99	27.54	4422.55
Aquatic soil	808.96	118.01	3291.18	275.83	76.13	24.50	4594.61
Saline soils	1480.95	38.20	2251.09	744.00	308.02	7.55	4829.81
Artificial soil	2425.38	0.00	1039.72	65.94	163.80	3.35	3698.18
Lakes and reservoirs	120.83	143.36	2082.38	0.00	9710.64	0.00	12,057.22

shows that different soil types had different impacts on the ESV, with great differences. The unit area value of lakes and reservoirs was the highest, while the unit area value of artificial soil was the lowest. In lakes and reservoirs, the proportion of unused land and construction land was 0, as water bodies occupied the majority of the region and had a strong ecological function, and so this soil type had a high ESV per unit area. In the artificial soil, cropland accounted for 74.34%, which was due to long-term production activities. After fertilization, irrigation, farming, and other agricultural measures, the amount of cropland rose, as did its value, in line with the characteristics of artificial soil.

3.4.5. Impact of Elevation on the ESV

The elevation data were reclassified into four levels: below 1000 m, 1000–1200 m, 1200–1400 m, and above 1400 m. As shown in

Table 9. ESV per unit area of different elevation levels. Unit: yuan/hm².

Elevation (m)	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
<1000	1705.79	1041.27	1704.86	13.61	118.34	0.81	4584.68
1000–1200	1220.74	1648.12	2117.71	0.00	42.71	10.86	5040.14
1200–1400	676.91	2792.27	2412.73	18.06	60.53	18.40	5978.90
>1400	1041.44	2423.38	2228.81	12.06	26.64	2.14	5734.47

, the ESV in higher-elevation locations was frequently higher than that in regions with lower elevations. The ESV was the lowest in areas below 1000 m, and the ESV per unit area of 1200–1400 m was the highest, which was 1.25 times that of the area below 1000 m. Construction land was extensively dispersed in the low-altitude region, as per a comparison of land-use. A total of 65.10% of the construction land was distributed below 1200 m and human activity had a significant influence. Only 8.97% of the forestland was distributed lower than 1000 m, while 67.73% of the forestland was distributed in the area above 1200 m. This shows that the high-altitude area had less human intervention, a good ecological environment, and a strong ecological function.

3.4.6. Impact of Population Density on the ESV

The population density data were reclassified into four levels: less than 40 people/km², 40–80 people/km², 80–120 people/km², and more than 120 people/km².

Table 10. ESV per unit area of different population densities. Unit: yuan/hm².

Population Density (People/km2)	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
<40	664.01	5381.45	1548.09	28.57	56.57	3.22	7681.91
40–80	1059.75	1079.70	2653.85	9.36	41.65	14.88	4859.19
80–120	1192.12	1112.54	2454.96	5.71	40.21	11.85	4817.41
>120	1648.14	1345.15	1621.99	18.80	96.01	5.01	4735.11

shows that, with the increase in population density, the ESV per unit area showed a decreasing trend. The ESV per unit area in areas with a population density below 40 people/km² was 1.62 times that in areas with a population density above 120 people/km². In the area with less than 40 people/km², forestland and grassland accounted for 46.34% and 32.12%, respectively, and construction land accounted for 0.19%. Wetland was also mainly distributed here, while in areas with more than 120 people/km², cropland accounted for 50.52%, construction land accounted for 2.77%, and 44.54% of construction land was distributed here, indicating that, in areas with less population, human activity intensity was weak, vegetation coverage was high, the ecosystem was stable, and service value was great, but, in densely populated regions, people’s intense social and economic activity would undoubtedly have some negative effects on the ecological environment, which would subsequently restrict the ESV.

3.4.7. Impact of Roads on the ESV

The main roads were established to obtain regions within 2 km, 2–4 km, 4–6 km, and more than 6 km away from a road.

Table 11. ESV per unit area at different distances from roads. Unit: CNY/hm².

Distance from the Road (km)	Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
<2	1081.81	1436.47	2224.81	5.60	78.77	10.13	4837.57
2–4	1035.45	1727.92	2389.52	4.54	60.21	10.73	5228.37
4–6	1047.78	1802.36	2373.51	5.17	59.98	9.45	5298.25
>6	1068.62	2413.54	2103.04	13.35	47.65	10.31	5656.51

shows that, with the increase in the distance from the main roads, the ESV per unit area continued to increase. From areas within 2 km to areas beyond 6 km, the ESV per unit area increased by 16.93%. Among them, forestland accounted for only 12.37% in areas within 2 km from a road and 20.78% in areas more than 6 km away from a road, construction land accounted for 6.09% in areas within 2 km and 0.71% in areas more than 6 km away from a road, and the unit area wetland value in areas more than 6 km away from a road was much higher than other regions.

3.5. Weight Analysis of Influencing Factors and Comprehensive Evaluation

3.5.1. Weighting Analysis of Influencing Factors

To comprehensively reflect how the aforementioned seven natural and societal factors affect the ESV in the area, the unit area value of each grade of each factor was calculated based on the above quantification. The entropy weight method index calculation formula was used to standardize the unit area value of the different grades of each factor through the standardization formula to calculate its information entropy. Finally, information entropy was used to calculate each factor’s weight.

Figure 8 shows that the seven natural and social factors influenced the ESV in descending order of weight: population density (0.309), rainfall (0.165), soil type (0.116), slope (0.109), elevation (0.106), annual average temperature (0.100), and distance from the road (0.096). The highest weight was given to population density and the lowest weight was given to distance from a road. Therefore, human activities disturbed the ecosystem the most among the many influencing factors. In addition, due to the limited space occupied by roads, they provided less interference with the ecological environment in areas far away from roads.

3.5.2. Comprehensive Evaluation of Influencing Factors

To determine the comprehensive influence size of the seven factors on the ESV of the study area, the weights of the factors were multiplied by the corresponding ESV per unit area of each level based on the weights of the factors previously calculated. According to Jenks, four categories were used to divide the area: mild-impact zone, moderate-impact zone, severe-impact zone, and extreme-impact zone. The zoning map of the impact of the ESV in the area was generated.

Figure 9 shows that the influence of natural and social factors on the ESV showed obvious geographical distribution characteristics, which were mainly reflected in the higher degree of influence in the northern area and the lower degree of influence in the southern area. From north to south, there was a weakening tendency that was constant. Comparing the ESV distribution and land-use status, it was found that mild- and moderate-impact zones were basically distributed in zones with vigorous vegetation growth and a high ESV, indicating that the low-impact zones had high ecological environment quality. In more severe impact zones, strong disturbances from artificial and natural conditions affected the performance of ecological functions. A degree of ecological function would be impacted by excessive human interference.

3.6. Forecast of the ESV Trend

Before making a forecast, the land-use status maps in 2000 and 2010 were used to generate the forecast map of the study area in 2020 and for comparison with the real land-use map from 2020. After accuracy verification, the Kappa coefficient was 0.75, which showed that the forecast result based on the CA–Markov model in the area was relatively accurate and this method could be used to predict its future changes. In turn, the land-use status of the region in 2030 was forecasted (Figure 10).

Table 12. Forecast of ESV in the Loess Plateau of Northern Shaanxi in 2030. Unit: CNY 100 million.

Cropland	Forestland	Grassland	Wetland	Water Bodies	Unused Land	The Total Value
85.67	171.10	171.62	0.65	4.54	0.78	434.36

shows that the total ESV in 2030 was CNY 43.436 million, CNY 39 million lower than that in 2020. The value provided by grassland decreased the most, while the value of forestland increased the most. Among them, grassland, wetland, and unused land values fell by CNY 315, 16, and 4 million, respectively. The values of forestland, cropland, and water bodies increased by CNY 222, 62, and 4 million, respectively. Comparing the predicted land-use status, it was found that the construction land area in 2030 had increased by 22.30% compared to 2020. The ecological resources were quickly consumed because of the rapid growth of construction land, which reduced the ESV and weakened ecological functions. Therefore, in the future, it is necessary to further strengthen the protection and construction of the ecological environment of the area to prevent the appearance of negative results.

4. Discussion

4.1. Temporal and Spatial Change Mechanism of the ESV

In this study, we found that there were significant temporal dynamics and spatial heterogeneity of the ESV in the Loess Plateau region of Northern Shaanxi. This region has seen a large change in land-use status in the last 20 years, and this has led to changes in the ecosystem services of the region and their values. Construction land is a nonecosystem service site and therefore has a null ESV [58,59,60]. Additionally, within the region from 2000 to 2020, the ESV has been decreasing with increasing construction land and cropland. Numerous studies have found that urban expansion and the increase in cropland come at the cost of high-ESV lands such as forestland and grassland [61,62]. These lands occupy and affect other ecological lands, leading to the decrease in the ESV. This is consistent with the results of our study. Additionally, in our study, the ESV in the study area showed a trend of being high in the south and low in the north. After the analysis, it was found that Northern Shaanxi is a relatively rich region of coal resources in China, and coal resources are more concentrated in the northern part of the study area, near Shenmu City, Fugu County, Yulin City, etc. The concentrated coal mining has caused damage to the surrounding ecological environment to some extent. Qian et al. found that mining activities in the southern Qilian Mountains caused significant ESV losses to the grassland and wetland [63]. Gao et al. found that the dramatic expansion of mining land had destroyed large areas of forestland and reduced the ESV [64]. The southern part of the study area, Huanglong County, is one of the eight major forest areas in China. It is also a protective forest area of the Loess Plateau. The proportion of forestland in the area was above 85% in all years, and the vegetation grows vigorously. The study of this region is important for the conservation of regional biodiversity and the maintenance of ecological balance [65]. Therefore, the ESV was higher in the southern part of the study area. In addition, the spatial autocorrelation model was used to analyze the spatial aggregation pattern and association pattern of the ESV, which can further reveal the spatial differentiation characteristics of the ESV in the study area [66]. In this study, a high-high aggregation area in the south and a low-low aggregation area in the northcentral part were identified, respectively. This indicates that the ESV in the study area has obvious spatial aggregation characteristics.

4.2. Influencing of Various Factors on the ESV

The constant changes in the ESV are influenced by multiple natural and social factors [67,68,69]. In this study, the effect of each factor on the regional ESV was found to be negatively correlated with population density and the ESV. The effect of soil type on the ESV was determined by the nature of the soil itself. The other five factors were positively correlated with the ESV. These findings are consistent with previous studies [70,71]. Analysis of the weight of each factor on the ESV using the entropy weighting method revealed that population density, rainfall, and soil type in the study area have a greater effect on the ESV, and the distance from a road has the least effect. Human activities disturb local ecosystems, altering the ecological environment and affecting ecological functions, leading to a decrease in the ESV. Guo et al. found that population density was the most important factor influencing the ESV in the region of Funiu Mountain [72]. In addition, areas with high rainfall have good hydrothermal conditions that are more suitable for vegetation growth and have more complex ecosystems. Meanwhile, low rainfall causes the ecological function of the area to be weakened. Xie et al. found that the ESV in China gradually decreased from the rainfed southeast to the arid northwest [73]. Zhao et al. found that precipitation has an important impact on ecosystem services in the Yangtze River Delta region [74]. In this study, although roads are also the result of human production and construction, they only have a small impact on the ESV in a limited area of their proximity due to the limitation of their own functional performance. Therefore, the degree of impact is relatively small. After a comprehensive evaluation by using the weights of each factor, it was found that the ESV was high in the areas with a low impact and low in the areas with a high impact in the study area. This indicates that the ESV is higher in areas with abundant precipitation, sufficient sunshine, high elevation, a gentle slope, and low human intervention. However, excessive human intervention can affect ecosystem function [75]. Therefore, if the various factors mentioned above are adjusted in the future, a good ecological environment in the study area can be maintained.

4.3. Trend Forecast Analysis of the ESV

Forecasting of the regional ecological environment enables us to anticipate the future development trend of the ecological environment. Forecasting furthers our understanding of the regional ecological security situation in order to facilitate constructive suggestions for ecological environment improvement and protection [76]. CA models have powerful computational capabilities to solve complex system problems and have been widely used in land-use [77,78]. However, it is difficult for a single model to achieve a reasonable prediction of land-use change. Some scholars have coupled the CA model with logistic functions, artificial neural network models, and SD models to carry out research on the dynamic prediction of land-use [79,80,81]. However, the above integrated models are difficult to achieve in terms of nonlinear representation, parameter optimization, operability, and so on [82,83,84]. However, the CA–Markov model is more scientific and applicable. It can make the prediction results more accurate [85]. In this study, the CA–Markov model was used to impose the corresponding constraints to achieve the prediction of land-use in the study area and calculate the ESV. The results show that the ESV of the study area will decrease by 2030. This is mainly due to the continuous increase in construction land. Therefore, further rational planning of each land-use type is needed to promote the improvement of ecological service quality.

4.4. Policy Implications

First, after the analysis, we found that human activities have the strongest influence on the ecological environment within the Loess Plateau area in Northern Shaanxi. Therefore, in order to further control the interference of human activities, we need to establish reasonable ecological protection zones within the study area. Human activities are restricted within the protected areas to protect the ecological environment and biodiversity. Meanwhile, in regional territorial spatial planning, three control lines of the ecological protection red line, permanent basic agricultural land, and urban development boundary should be determined. The government should develop the regional economy and ecological security at the same time so as to coordinate regional development and ecological stability and prevent economic development at the expense of ecological environment [86].

Second, according to the results of the comprehensive zoning evaluation of each impact factor, the zoning is based on the degree of impact. Different environmental improvement measures were implemented in the study area for different impact zones. The extremely heavy-impact area and heavy-impact area are mainly located in the northcentral part of the study area. This area is characterized by a steep terrain, strong human interference, and relatively fragile ecological environment. Therefore, the utilization rate of construction land and vegetation coverage should be further improved in this area, and ecological protection projects should be vigorously carried out. Thus, the stability and health of the ecosystem should be maintained. For the moderately and lightly affected areas distributed in the south, we should continue to do a good job of protecting the ecological environment of the area, controlling large areas of unreasonable grazing, and protecting water resources.

Third, Northern Shaanxi is extremely rich in mineral resources. However, extensive mining activities have caused groundwater pollution, surface landscape fragmentation, ground subsidence, soil quality decline, and other environmental problems [87]. Measures should be taken to protect and restore the ecological environment of mining areas. Hu et al. proposed measures alongside mining to restore the ecological environment of the mining area [88]. Bi et al. proposed the use of microorganisms to improve the quality of the soil in the mining area [89]. Therefore, the corresponding methods for the ecological management of mining areas need to be proposed according to the mining characteristics of the study area so that the ecological condition of the entire Northern Shaanxi Loess Plateau can be improved.

4.5. Limitations and Future Work

The findings of this study offer some reference value for the future sustainable development and ecological environment of the region; however, there are some limitations to this study. The differences between different regions were affected by natural, social, and economic factors. In this study, the ESV’s equivalent factor was revised by calculating regional and national climate productivity, which is not comprehensive enough. In addition, in analyzing the factors affecting the ESV, this study selected annual rainfall, population density, roads, and other factors for analysis. Although both natural and social factors were included, the consideration of the factors still had some limitations. Therefore, there will be some differences between the established impact zone and the actual situation. Finally, when using the CA–Markov model for trend forecasting, the elevation, slope, water bodies, distance from the road, and distance from a water body were added, and the policy factors for future development in the area were not fully taken into account. Combined with the model itself being restricted, there will be some differences between the predicted results and the actual development in the future.

5. Conclusions

In this study, the temporal and spatial changes in the Loess Plateau of Northern Shaanxi based on the calculation of the ESV were analyzed, and an impact zoning map was established by quantifying various influencing factors. Finally, the future development trend of the ESV was predicted. The following conclusions can be reached: (1) The land-use types were mainly grassland, forestland, and cropland, and the combined area of the three types covered more than 95% of the area. In addition, construction land in the study area increased the fastest, with a 4.6-fold increase from 2000 to 2020. (2) From 2000 to 2020, the ESV decreased continuously, and a pattern of being high in the south and low in the north was evident in the spatial features. (3) The weighting analysis of influencing factors showed that population density (0.309) > rainfall (0.165) > soil type (0.116) > slope (0.109) > elevation (0.106) > annual average temperature (0.100) > and distance from a road (0.096). By establishing the impact zones, it was found that the degree of impact decreases from south to north in the area. (4) It is predicted that the ESV will remain reduced by 2030 compared to 2020 in the area.

The ESV is affected by a combination of factors. Natural factors create the foundation of the ecosystem, and with rapid socioeconomic development, humans will continue to impose on this foundation, making the ecosystem more and more unstable. Positive human behaviors lead to stronger ecosystem functions, while negative human behaviors disrupt the ecological balance and cause the ecosystem to evolve in an unsustainable direction. In this study, the factors that influence the ESV were considered and different impact areas were classified to understand the regional ecological quality, identify the root causes that limit the function of ecosystems, and predict the future development trend. The research results provide a scientific basis for ecological environmental protection and policy formulation in the Loess Plateau of Northern Shaanxi.