*Article* **Spatio-Temporal Evolution and Driving Factors of Ecosystem Service Value of Urban Agglomeration in Central Yunnan**

**Lei Yang and Fenglian Liu \***

Institute of Land & Resources and Sustainable Development, Yunnan University of Finance and Economics, Kunming 650221, China

**\*** Correspondence: zz2105@ynufe.edu.cn; Tel.: +86-137-2010-8008

**Abstract:** Urbanization and human activity have recently resulted in land use/cover change (LUCC), which has had a detrimental effect on the biological environment, on keeping the ecosystem's sustainable growth and on comprehending the ecosystem's quality and changes over the past 20 years in the central Yunnan urban agglomeration. The equivalent factor method and hotspot analysis were used to analyze the spatio-temporal changes in land use and ecosystem service value (ESV) in the urban agglomerations of central Yunnan province, and the effects of land use change on ESV were then examined. This study is based on the grid data of land use in 2000, 2005, 2010, 2015, and 2020. Finally, Geodetector was used to investigate the possible causes of ESV. The results showed that: (1) The urban agglomerations in central Yunnan's land-use structure and pattern clearly changed between 2000 and 2020, with continual declines in grassland, cultivated land, and woodland, and constant increases in construction land. There was significant growth in both speed and area. (2) The average ESV of the land decreased consistently, the hotspot areas shrank, and the cold-spot areas grew as the ecosystem service function declined and the total amount of ESV decreased by 1.517 billion Yuan. These events were mostly explained by an increase in construction land and a decrease in grassland, cultivated land, and woodland. (3) The synergistic effect of numerous factors is what causes the change in ESV in the urban agglomerations of central Yunnan. The key forces behind ESV change in the research area were land-use intensity, normalized difference vegetation index (NDVI), slope, and people density. The results can help decision makers establish policies for ecological conservation and land use.

**Keywords:** land-use change; ecosystem service value; hotspot analysis; Geodetector; central Yunnan urban agglomeration

### **1. Introduction**

The level of social and economic development worldwide has altered substantially as urbanization has accelerated. Humans continue to seek development from nature during this phase. Numerous issues, including climatic warming, biodiversity loss, and aggravated soil erosion, have been brought on by the growing population, unsustainable land use patterns, and excessive demand for resources [1,2]. As a result, ecosystem functions have been severely compromised, endangering the long-term welfare of people and the sustainable development of ecosystems. In this setting, it is crucial for human survival and growth to control ecosystem quality and change, identify internal root causes of these issues, and implement appropriate solutions [3]. There are numerous indicators to track the dynamic changes in ecosystems, but few can accurately capture the value and expense of ecological deterioration [4–6]. Ecosystem service value may quantitatively evaluate ecosystem service value in monetary units, indicating the ecological consequences caused by diverse activities and effectively connecting ecosystem study with management decisions [7–11]. This concept has increasingly gained popularity in geography, ecology, and other related sciences.

**Citation:** Yang, L.; Liu, F. Spatio-Temporal Evolution and Driving Factors of Ecosystem Service Value of Urban Agglomeration in Central Yunnan. *Sustainability* **2022**, *14*, 10823. https://doi.org/ 10.3390/su141710823

Academic Editors: Ronald C. Estoque, Kikuko Shoyama and Rajarshi Dasgupta

Received: 29 June 2022 Accepted: 26 August 2022 Published: 30 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Costanza et al. adopted the research method of economics, introduced "willingness to pay", and conservatively estimated the value of ecosystem services per unit area of different types of ecosystems and the economic value of global ecosystem services. Furthermore, they divided ecosystem services into 17 categories, which laid a theoretical foundation for the follow-up research [12]. The Millennium Ecosystem Assessment issued by the United Nations in 2005 formally established the concept of ecosystem services. They defined ecosystem services as the benefits that people directly or indirectly obtained from ecosystems and divided them into four categories: supply (food production, raw material supply, etc.), regulation (gas regulation, climate regulation, etc.), support (soil conservation, biodiversity maintenance, etc.) and culture (leisure and entertainment, aesthetic landscape, etc.) [13]. Since then, research on ecosystem service value has become more diversified, and research methods and framework systems have become more mature; The research objects include the forests [14,15], farmland [16], wetlands [17,18], grassland [19], river basin [20,21], city [22], island [23], coastal zone [24], etc. The research method is primarily quantitative research based on remote sensing technology [25]. The research content also extends from the detailed research of ecosystem service value to the correlation analysis with land use change, landscape pattern and landscape ecological risk, etc. [26–28] In addition, the existing research has gradually changed from evaluating and analyzing the total ecosystem service value of the research object to studying the value of a specific ecosystem or a single ecosystem service [29–31], from studying ESV itself to the comprehensive application of ESV (such as proposing an ecological compensation mechanism based on the ecosystem service value) [32].

In conclusion, the research on ESV is generally mature and productive; however, there are still two areas that require development. To begin with, the majority of prior research subjects were chosen from normal ecosystems, but there was a dearth of direction for the ecological preservation of administrative units. Especially in urban agglomerations with more complex structures and functions, the contradiction between social and economic development and ecosystem protection is more prominent. Secondly, the ecosystem is highly complex, and its changes resulting from the synergy of multiple factors. However, there are few studies on the changes in ecosystem structure and function of urban agglomerations, and the relationship between rapid urban development and ecosystem changes is still unclear. The potential driving elements of ESV have not been fully examined, especially the research on the influence degree and spatial differentiation of different types of driving factors on ESV is scarce. As a result, using the urban agglomerations in central Yunnan as the research object, this paper first examines the temporal and spatial evolution characteristics of land use and ecosystem service value, then assesses the effects of land use change on ecosystem service value, and finally, investigates the main driving forces of ecosystem service value. The research results can provide quantitative information for decision makers to grasp the ecological status of the urban agglomeration in central Yunnan to optimize the structure of land use and reduce the loss of ESV caused by unreasonable land use patterns. This will help to promote socioeconomic–ecological sustainable development in central Yunnan and urban agglomerations with similar conditions [33].

### **2. Materials and Methods**

### *2.1. Study Area*

The Central Yunnan Urban Agglomeration is one of the 19 state-level urban agglomerations located in the central part of Yunnan Province (101–104.5◦ E, 24–26.5◦ N), which is a typical plateau mountainous urban agglomeration with a topography dominated by mountains and intermountain basins (Figure 1). It is also the most economically developed, densely populated, and strongly developed region in Yunnan Province, with the highest concentration of intermountain basins. According to the Central Yunnan Urban Agglomeration Development Plan, which was published in July 2020, the region's scope includes the entire territory of Kunming, Qujing, Yuxi, and Chuxiong, as well as seven counties and cities in the north of Honghe Prefecture, totaling 49 counties (urban areas) with a combined area of 111,400 square kilometers, or 28.3% of Yunnan Province's total area. The urban agglomeration in central Yunnan has grown quickly in recent years, with a GDP of 1507.394 billion Yuan in 2020, 11.45 times that of 2000 (136.105 billion Yuan), and making up about 61.47 percent of the total GDP of Yunnan Province. The pace of urbanization has also accelerated significantly, rising from 30.62 percent in 2000 to 59.90 percent in 2020. area. The urban agglomeration in central Yunnan has grown quickly in recent years, with a GDP of 1507.394 billion Yuan in 2020, 11.45 times that of 2000 (136.105 billion Yuan), and making up about 61.47 percent of the total GDP of Yunnan Province. The pace of urbanization has also accelerated significantly, rising from 30.62 percent in 2000 to 59.90 percent in 2020.

Agglomeration Development Plan, which was published in July 2020, the region's scope includes the entire territory of Kunming, Qujing, Yuxi, and Chuxiong, as well as seven counties and cities in the north of Honghe Prefecture, totaling 49 counties (urban areas) with a combined area of 111,400 square kilometers, or 28.3% of Yunnan Province's total

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 3 of 20

**Figure 1.** Overview of the study area.

#### **Figure 1.** Overview of the study area. *2.2. Data Sources and Methods*

#### *2.2. Data Sources and Methods*  2.2.1. Data Source

2.2.1. Data Source The five periods of land use/cover remote sensing monitoring data used in this study from 2000 to 2020 were downloaded from the Center of Resources and Environmental Science and Data, Chinese Academy of Sciences (http://www.resdc.cn, accessed on 25 May 2022), among which the land use grid data in 2000, 2005 and 2010 were interpreted by Landsat TM/ETM remote sensing image data, and the land use grid data in 2015 and 2020 were obtained by interpreting Landsat 8 remote sensing image data, with a spatial resolution of 30 m and an interpretation accuracy of over 90%, including 6 primary land types (cultivated land, woodland, grassland, water area, construction land, and unused land) and 25 secondary land types [34]. DEM (digital elevation model) is a branch of digital terrain model, whose purpose is to describe the terrain surface morphology and other ground elevation information in digital form [35]. The DEM images of urban agglomeration in central Yunnan in 2020 used in this study were downloaded from the geospatial data cloud (http://www.gscloud.cn, accessed on 25 May 2022), with a spatial resolution of The five periods of land use/cover remote sensing monitoring data used in this study from 2000 to 2020 were downloaded from the Center of Resources and Environmental Science and Data, Chinese Academy of Sciences (http://www.resdc.cn, accessed on 25 May 2022), among which the land use grid data in 2000, 2005 and 2010 were interpreted by Landsat TM/ETM remote sensing image data, and the land use grid data in 2015 and 2020 were obtained by interpreting Landsat 8 remote sensing image data, with a spatial resolution of 30 m and an interpretation accuracy of over 90%, including 6 primary land types (cultivated land, woodland, grassland, water area, construction land, and unused land) and 25 secondary land types [34]. DEM (digital elevation model) is a branch of digital terrain model, whose purpose is to describe the terrain surface morphology and other ground elevation information in digital form [35]. The DEM images of urban agglomeration in central Yunnan in 2020 used in this study were downloaded from the geospatial data cloud (http://www.gscloud.cn, accessed on 25 May 2022), with a spatial resolution of 30 m, and were derived from the stitching and cropping of multiple remote sensing images. The grain yield and price data used to calculate the value of ecosystem services were obtained from the Yunnan Statistical Yearbook, the National Compilation of Costs and Benefits of Agricultural Products, and official government websites.

Considering the natural conditions and social development level of the urban agglomeration in central Yunnan, and combining it with the related research result [36], nine

drivers were selected from the natural factors and human factors, including elevation, slope, average annual temperature, average annual precipitation, NDVI, soil erosion, population density, GDP, and land-use intensity in 2020. Elevation and slope data were extracted from DEM images; temperature, precipitation, NDVI and soil erosion data were downloaded from the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on 25 May 2022); population density and GDP data were obtained from the Yunnan Statistical Yearbook in 2020; and land-use intensity data were calculated by Equation (3).

### 2.2.2. Land-Use Dynamic Degree

The dynamic degree of single land use is an essential measure of the speed and magnitude of regional land-use change, which can effectively reflect the increase or decrease in the area of various land-use types in the study area in a certain period [37], the expression is:

$$K = \frac{\mathcal{U}\_b - \mathcal{U}\_a}{\mathcal{U}\_a} \times \frac{1}{T} \times 100\% \tag{1}$$

where *K* is the dynamic degree of a specific land-use type in the *T* period, which is used to measure the degree to which the change in the land-use type is disturbed by other factors, and the larger the value of *K*, the more unstable this land-use type is in the *T* period; *U<sup>a</sup>* and *U<sup>b</sup>* indicate the area of this type of land use at the beginning and end of the study period, respectively.

### 2.2.3. Land Use Transfer Matrix

The land use transfer matrix is a quantitative research method through systematic analysis, which is mainly used in related research on land-use change. The mathematical expression is [38]:

$$\mathbf{S}\_{ij} = \begin{bmatrix} \mathbf{S}\_{11} & \mathbf{S}\_{12} & \dots & \mathbf{S}\_{1n} \\ \mathbf{S}\_{21} & \mathbf{S}\_{22} & \dots & \mathbf{S}\_{2n} \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{S}\_{n1} & \mathbf{S}\_{n2} & \dots & \mathbf{S}\_{nn} \end{bmatrix} \tag{2}$$

where *Sij* represents the area change in a certain land-use type from the beginning to the end of the research period, and *n* represents the number of land-use types involved in the study.

### 2.2.4. Land-Use Intensity Index

The land-use intensity index reflects the degree of land-use change caused by human activities, and regional land development and utilization to a certain extent, and measures the depth and breadth of regional land use. The calculation formula is [39]:

$$L = 100 \times \sum\_{i=1}^{n} \frac{A\_i P\_i}{A\_t} \tag{3}$$

where *L* is the comprehensive land-use intensity index of the study area; *n* is the number of land-use types, *A<sup>i</sup>* is the area of different land-use types, *P* represents the land-use intensity parameter, and *A<sup>t</sup>* represents the total area of the study area. According to the previous research results, the unused land, woodland, grassland, water area, arable land, and construction land are divided into four grades and assigned values of 1–4; *C<sup>i</sup>* is the area of the type *i* land-use type.

### 2.2.5. Calculation of Ecosystem Services Value

Costanza et al. proposed the equivalent factor method, which makes the assumption that the ecosystem service value per unit area of the same ecosystem type is constant; the total ecosystem service value can be obtained by multiplying this value by the area of each ecosystem type in the study area. Based on this, Chinese academics Xie et al. enhanced this method to determine the ecosystem service value equivalent per unit area acceptable for China [40,41]. Referring to the research results of Xie et al. and The Millennium Ecosystem Assessment [13], this study categorizes ecosystem services into four primary service kinds supply, support, regulation, and culture—and 11 secondary service types. Combined with the situation of urban agglomeration in central Yunnan, the equivalent scale of ecosystem service value per unit area of ecosystem in the central Yunnan urban agglomeration is shown in Table 1, and the construction land equivalent coefficient is referred to in the study of Deng Shuhong [42]. The economic value of one ecosystem service value equivalent factor can be defined as 1/7th of the market value of the average grain yield per hectare of farmland in the study area under natural conditions without external interference [43]. The calculation equation is as follows:

$$E = \frac{1}{7} \times \sum\_{a=1}^{n} \frac{S\_a P\_a Q\_a}{S} \tag{4}$$

where *E* is the ecosystem service value of a standard equivalent factor, *P<sup>a</sup>* is the average price of *a* grain in the study area during the study period, *Q<sup>a</sup>* is the grain yield per unit area of *a* grain, *S<sup>a</sup>* and *S* is the sum of the sown area of grain *a* and the sown area of three-grain crops. The value of ecosystem services for one common equivalent factor in the central Yunnan urban agglomeration was 1149.34 Yuan/hm<sup>2</sup> . The formula for calculating the value of ecosystem services is as follows:

$$ESV = \sum\_{q=1}^{n} \mathcal{S}\_q \times V \mathcal{C}\_q \tag{5}$$

$$ESV\_k = \sum\_{q=1}^{n} \mathcal{S}\_q \times V \mathcal{C}\_{qk} \tag{6}$$

$$AESV = \frac{\sum\_{q=1}^{n} \mathbb{S}\_q \times V \mathbb{C}\_q}{\sum\_{q=1}^{n} \mathbb{S}\_q} \tag{7}$$

where *ESV*, *ESV<sup>k</sup>* and *AESV* represent the total ecosystem service value, the individual ecosystem service value, and the average ecosystem service value, respectively; *q* is the land type *q*; *S<sup>q</sup>* is the area of the land type *q*; *VCq*, *VCqk* represent the ecosystem service value per unit area of land type *q* and the value coefficient of the *k*th ecosystem service, respectively [26].

**Table 1.** Equivalent scale of ecosystem service value per unit area of ecosystem in central Yunnan urban agglomeration.


### 2.2.6. Hotspot Analysis

The hotspot analysis method can judge whether there is a clustering of high and low values of ecosystem service values in the study area and then identify the spatial distribution locations of the hotspot and cold-spot areas of ecosystem service values to reveal the spatial difference of ecosystem service capacity provided by various parts of the study area [44]. The formula is as follows:

$$G\_i^\* = \frac{\frac{\sum\_{j=1}^n w\_{ij} x\_j - \overline{X} \sum\_{j=1}^n w\_{ij}}{\mathcal{S} \sqrt{\frac{\left[n \sum\_{j=1}^n w\_{ij}^2 - \left(\sum\_{j=1}^n w\_{ij}\right)^2\right]}{n-1}}} \tag{8}$$

$$\overline{X} = \frac{1}{n} \sum\_{i=1}^{n} x\_i \tag{9}$$

$$S = \sqrt{\frac{1}{n}} \sum\_{i=1}^{n} \mathfrak{x}\_i^2 - (X)^{-2} \tag{10}$$

where *G* ∗ *i* index is positive, it is a hotspot area with more ESV increase, and if the *G* ∗ *i* index is negative, it is a cold-spot area with more ESV loss; *n* is the number of grids; *x<sup>i</sup>* , *x<sup>j</sup>* is the ecosystem service values of the *i*th and *j*th grids respectively.

### 2.2.7. Sensitivity Analysis

The sensitivity index is similar to the principle of the elasticity coefficient in economics. It is used to verify the reliability of the ESV results calculated using the ecosystem service value coefficient of different land-use types and then to analyze the impact of different land-use type changes on the value of ecosystem services [45]. The formula is as follows:

$$\text{CS} = \left| \frac{(ESV\_b - ESV\_a) / ESV\_a \times 100\%}{\left(VC\_{jk} - V\mathbb{C}\_{ik}\right) / VC\_{ik}} \right| \tag{11}$$

where *CS* is the elasticity coefficient, which indicates the degree of influence of different land-use types on the ecosystem service value; *ESV<sup>b</sup>* and *ESV<sup>a</sup>* indicate the ecosystem service values at the end of the study period and the beginning of the study period; *VCik*, *VCjk* represent the ecosystem service value coefficient of the *K* land-use type before and after adjustment, respectively.

### 2.2.8. Analysis of Driving Factors

Wang J.F.'s team first proposed Geodetector. It is a statistical analysis model, which mainly includes four parts: ecology, factor, risk, and interaction detectors. Each detector has different functions. The factor detector mainly indicates the explanatory strength of the independent variable to the dependent variable and the spatial heterogeneity of the dependent variable [46], The calculation method is as follows:

$$q = 1 - \frac{\sum\_{h=1}^{L} N\_h \sigma\_h^2}{N\sigma^2} = 1 - \frac{SSW}{SST} \tag{12}$$

where *q* denotes the explanatory power of a factor on ESV, and the value is between 0 and 1. The larger the value of *q*, the stronger the explanatory power of the factor on ESV; *h* is the partition number of the independent variable; *L* is the total number of partitions *N<sup>h</sup>* and *N* is the total number of samples in each partition, and the whole region, respectively; *σ* 2 *h* and *σ* <sup>2</sup> are the variances of each partition and the variance of ESV in the whole region; *SSW* and *SST* are the sum of variance within the stratum and the total variance in the whole region, respectively.

### **3. Results** From 2000 to 2020, the most crucial land-use type in the central Yunnan urban ag-

**3. Results** 

#### *3.1. Analysis of Land Use Change* glomeration was woodland (Figure 2), accounting for about half of the total area, followed by grassland and cultivated land. During the study period, the land-use change showed

*3.1. Analysis of Land Use Change* 

whole region, respectively.

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 7 of 20

From 2000 to 2020, the most crucial land-use type in the central Yunnan urban agglomeration was woodland (Figure 2), accounting for about half of the total area, followed by grassland and cultivated land. During the study period, the land-use change showed "three increases and three decreases, that is, cultivated land, woodland, and grassland decreased, among which the grassland area decreased the most (a total decrease of 716.3 km<sup>2</sup> ); the area of construction land, water area, and unused land increased, among which the growth of construction land was more pronounced (1269.01 km<sup>2</sup> )". As can be seen from Table 2, the dynamic degree of various land-use types in 2015–2020 was significantly higher than that of the previous three research periods, indicating that the land use in the study area changed drastically during this period. Overall, the land-use dynamics of woodland, cultivated land, and grassland from 2000 to 2020 were −0.03%, −0.10%, and −0.12%, respectively, indicating that the change rates of these three land-use types accelerated sequentially during the study period. The land-use dynamics of construction land, the water area, and unused land were all positive, with the land-use dynamics of construction land having changed relatively significantly, increasing from 1.42% in 2000–2005 to 7.40% in 2015–2020. "three increases and three decreases, that is, cultivated land, woodland, and grassland decreased, among which the grassland area decreased the most (a total decrease of 716.3 km2); the area of construction land, water area, and unused land increased, among which the growth of construction land was more pronounced (1269.01 km2)". As can be seen from Table 2, the dynamic degree of various land-use types in 2015–2020 was significantly higher than that of the previous three research periods, indicating that the land use in the study area changed drastically during this period. Overall, the land-use dynamics of woodland, cultivated land, and grassland from 2000 to 2020 were −0.03%, −0.10%, and −0.12%, respectively, indicating that the change rates of these three land-use types accelerated sequentially during the study period. The land-use dynamics of construction land, the water area, and unused land were all positive, with the land-use dynamics of construction land having changed relatively significantly, increasing from 1.42% in 2000–2005 to 7.40% in 2015–2020.

the partition number of the independent variable; *L* is the total number of partitions and *N* is the total number of samples in each partition, and the whole region, respectively;

<sup>ଶ</sup> and ଶ are the variances of each partition and the variance of ESV in the whole region; *SSW* and *SST* are the sum of variance within the stratum and the total variance in the

**Figure 2.** Spatial distribution map of land use. In the legend, I = cultivated land, II = woodland, III = grassland, IV = water area, V = construction land, and VI = unused land. **Figure 2.** Spatial distribution map of land use. In the legend, I = cultivated land, II = woodland, III = grassland, IV = water area, V = construction land, and VI = unused land.


**Table 2.** Land use dynamics in central Yunnan urban agglomeration %.

The land use transfer matrix can describe the change in land-use quantity and transfer direction from dynamic and static aspects, which is of great significance for analyzing the change in the ecosystem service value caused by land-use changes. Based on the land-use data, Table 3 and Figure 3 were obtained by processing with the raster calculator in ArcGIS to systematically reflect the quantitative characteristics and spatial distribution of land-use changes. Table 3 shows that the main characteristics of land-use transfer in the central Yunnan urban agglomeration are the outflow of grassland, cultivated land, and woodland, and the inflow of construction land. The outflow of grassland goes mainly to woodland (2187.90 km<sup>2</sup> ), followed by cultivated land and construction land. The cultivated land was mainly converted into woodland (1208.45 km<sup>2</sup> ), grassland (1195.59 km<sup>2</sup> ), construction land (936.21 km<sup>2</sup> ), and water areas (105.05 km<sup>2</sup> ), which is due to the project of returning farmland to forests, lakes, and grasslands implemented in Yunnan in recent years, as well as the rapid development of the city, resulting in a sharp increase in the demand for construction land. Another reason for the increase in the water area is the comprehensive implementation of the strategy of "strengthening Yunnan with water" and the construction of a series of water conservancy projects such as the "Central Yunnan Water Diversion Project", which has increased the land and water area for water conservancy facilities. Construction land inflow and outflow are quite different, and the primary source is cultivated land (936.21 km<sup>2</sup> ). From the spatial distribution of land use transfer, the conversion of arable land and woodland to construction land mainly occurred in Anning City, Chenggong District, Guandu District, Qilin District, Mengzi City, Hongta District, and Chuxiong City.

**Land-Use Type Grassland Cultivated Land Construction Land Woodland Water Area Unused Land** Grassland - 1393.76 400.04 2187.90 98.44 12.88 Cultivated land 1195.59 - 936.21 1208.45 105.05 6.72 Construction land 23.25 142.28 - 23.67 10.25 1.39 Woodland 2086.96 1325.68 272.82 - 116.23 4.73 Water area 40.13 70.31 31.35 24.36 - 2.37 Unused land 12.30 7.00 0.14 2.55 1.03 -

**Table 3.** Land-use transfer matrix of central Yunnan urban agglomeration from 2000 to 2020 km<sup>2</sup>

.

**Figure 3.** Map of land-use change in central Yunnan urban agglomeration. **Figure 3.** Map of land-use change in central Yunnan urban agglomeration.

#### *3.2. Analysis of Ecosystem Service Value 3.2. Analysis of Ecosystem Service Value*

### 3.2.1. Time Series Change in Ecosystem Service Value 3.2.1. Time Series Change in Ecosystem Service Value

Over the past 20 years, the central Yunnan urban agglomeration's ecosystem service value (ESV) has decreased by 1.517 billion Yuan (Table 4). During the study period, with the exception of 2015 to 2020, the total ESV continued to decline, especially from 2010 to 2015, when ESV decreased rapidly, and the ecosystem service value decreased by nearly 1 × 109 Yuan. The reason for the decrease in the total amount of ESV lies in the decrease in cultivated land, woodland, and grassland area and the increase in the construction land area, resulting in the loss of ESV of 2.37 × 108 Yuan, 7.74 × 108 Yuan, 9.39 × 108 Yuan, and 19 × 108 Yuan, respectively. From 2015 to 2020, although the ESV of other land types decreased, the water area increased, and its ecosystem service value coefficient was signifi-Over the past 20 years, the central Yunnan urban agglomeration's ecosystem service value (ESV) has decreased by 1.517 billion Yuan (Table 4). During the study period, with the exception of 2015 to 2020, the total ESV continued to decline, especially from 2010 to 2015, when ESV decreased rapidly, and the ecosystem service value decreased by nearly <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>9</sup> Yuan. The reason for the decrease in the total amount of ESV lies in the decrease in cultivated land, woodland, and grassland area and the increase in the construction land area, resulting in the loss of ESV of 2.37 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, 7.74 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, 9.39 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, and 19 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, respectively. From 2015 to 2020, although the ESV of other land types decreased, the water area increased, and its ecosystem service value coefficient was significant, which made the ecosystem service value of the study area rise in 2020.

**Table 4.** ESV of different types of land from 2000 to 2020 1 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan.

cant, which made the ecosystem service value of the study area rise in 2020.


Unused land 0.04 0.00% 0.04 0.00% 0.04 0.00% 0.04 0.00% 0.04 0.00% Total 1937.92 100.00% 1935.48 100.00% 1930.36 100.00% 1920.35 100.00% 1922.75 100.00%

In terms of the value of individual ecosystem services (Table 5), from 2000 to 2020, the regulation service was dominant (67.74%), followed by the support service (22.15%), and their contribution rates to the total amount of ESV were close to 90%. Hydrological regulation and climate regulation were the most prominent types of secondary ecosystem services; followed by soil conservation, maintenance of biodiversity, gas regulation and environmental purification, at 218.58 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, 186.77 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, 182.13 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, 153.10 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan at the end of the study period, respectively. Water supply and nutrient cycle maintenance were relatively weak. From the change in individual ecosystem service value, the ecosystem service value of raw material production and climate regulation decreased by 12.54 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan and 8.44 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan, respectively. During the study period, only the ecosystem service value of hydrological regulation and nutrient cycle maintenance increased.

**Table 5.** Value of individual ecosystem services from 2000 to 2020 1 <sup>×</sup> <sup>10</sup><sup>8</sup> Yuan.


### 3.2.2. Spatial Change in Ecosystem Service Value

Compared with the total ESV in the study area, the average local ESV can better reflect the regional ecological environment quality by excluding the influence of the administrative area. With the help of the average ESV calculation formula, the average ESV of 49 counties (cities, districts) in the central Yunnan urban agglomeration and its spatial distribution (Figure 4) was obtained, and then the spatial characteristics of ESV distribution were explored. On the whole, the ESV of the central Yunnan urban agglomeration shows the spatial distribution characteristics being "high in the west and low in the east." The lowvalue areas are concentrated in the eastern region, with an expanding trend; the lower value areas are mainly distributed in the northeast; the higher value areas are distributed in the southwest and northwest of the central Yunnan urban agglomeration; the high-value areas are fewer and mainly concentrated in the central region. Specifically, high-value areas include the Xishan District, Chenggong District in Kunming, Chengjiang County, and Jiangchuan District in Yuxi, which correspond to the spatial distribution of water areas. In 2000, the low-value areas included seven counties in Kunming, Qujing, and Honghe (Wuhua, Qilin, Fuyuan, Luliang county, Luoping, Shizong county, and Luxi county), and in 2020, three counties in Kunming (Anning, Songming, and Guandu) were added. In general, the quality of the ecological environment has shown a downward trend. For example, Guandu District has changed from a high-value area to a low-value area due to the surge in construction land. However, the quality of the ecological environment has also improved in some areas, among which Eshan and Luquan have changed from middle-value areas to higher-value areas, and the value of ecosystem services has increased.

**Figure 4.** Spatial distribution of land average ESV of central Yunnan urban agglomeration. **Figure 4.** Spatial distribution of land average ESV of central Yunnan urban agglomeration.

To further explore the specific location of the spatial change in ESV, the hotspot analysis (Getis-Ord Gi\*) in the ArcGIS spatial statistical tools was used to analyze the cold and hotspots of ESV in the urban agglomeration of central Yunnan in 2000, 2005, 2010, 2015, and 2020. The results are shown in Figure 5: the hotspots of ESV in the central Yunnan urban agglomeration are distributed in the western and northeastern regions, while the cold spots are concentrated in the central region. Specifically, the hotspots in 2000 included Dongchuan District, Xuanwei City, and Shuangbai County. The sub-hotspot area included Chuxiong City, Xinping Yi, Dai Autonomous County, and Yuanjiang Hani, Yi, and Dai Autonomous County. The cold-spot areas include eight counties and districts in Kunming City (Wuhua District, Panlong District, Songming County, Xishan District, Guandu District, Chenggong District, Yiliang County, and Jinning District) and Chengjiang County in Yuxi City. The subcold-spot areas involve Fumin County and Anning City in Kunming City and Hongta District, and Jiangchuan District in Yuxi City. By 2020, the number of hotspots decreased (Xuanwei City), and the area of cold spots expanded (Fumin County and Anning City), indicating the degradation of ecosystem service functions in the study area. To further explore the specific location of the spatial change in ESV, the hotspot analysis (Getis-Ord Gi\*) in the ArcGIS spatial statistical tools was used to analyze the cold and hotspots of ESV in the urban agglomeration of central Yunnan in 2000, 2005, 2010, 2015, and 2020. The results are shown in Figure 5: the hotspots of ESV in the central Yunnan urban agglomeration are distributed in the western and northeastern regions, while the cold spots are concentrated in the central region. Specifically, the hotspots in 2000 included Dongchuan District, Xuanwei City, and Shuangbai County. The sub-hotspot area included Chuxiong City, Xinping Yi, Dai Autonomous County, and Yuanjiang Hani, Yi, and Dai Autonomous County. The cold-spot areas include eight counties and districts in Kunming City (Wuhua District, Panlong District, Songming County, Xishan District, Guandu District, Chenggong District, Yiliang County, and Jinning District) and Chengjiang County in Yuxi City. The sub-cold-spot areas involve Fumin County and Anning City in Kunming City and Hongta District, and Jiangchuan District in Yuxi City. By 2020, the number of hotspots decreased (Xuanwei City), and the area of cold spots expanded (Fumin County and Anning City), indicating the degradation of ecosystem service functions in the study area.

**Figure 5.** Analysis of ESV hotspots of central Yunnan urban agglomeration. **Figure 5.** Analysis of ESV hotspots of central Yunnan urban agglomeration.

#### 3.2.3. Sensitivity Analysis of Ecosystem Service Value 3.2.3. Sensitivity Analysis of Ecosystem Service Value

The sensitivity indices of different land-use types are different (Table 6). It is found that the sensitivity indices of each land-use type in the central Yunnan urban agglomeration are woodland, grassland, water area, cultivated land, unused land, and construction land in descending order, which is consistent with the contribution of each land-use type to the ESV. The sensitivity index of woodland was the highest, remaining above 0.66 during the study period, followed by grassland, with a sensitivity index close to 0.2. Different from the other land types, the sensitivity index of construction land was negative in different years, indicating that the increase in construction land would lead to a decrease in the ecosystem service value in the study area; the sensitivity index of unused land was shallow, indicating that the change in unused land had no significant effect on the change in ESV. In addition, the sensitivity indices of land-use types in each period were less than 1, indicating that ESV is inelastic to land-use changes. The study uses ecosystem service value coefficients of different land types suitable for the study area, and the ESV calcula-The sensitivity indices of different land-use types are different (Table 6). It is found that the sensitivity indices of each land-use type in the central Yunnan urban agglomeration are woodland, grassland, water area, cultivated land, unused land, and construction land in descending order, which is consistent with the contribution of each land-use type to the ESV. The sensitivity index of woodland was the highest, remaining above 0.66 during the study period, followed by grassland, with a sensitivity index close to 0.2. Different from the other land types, the sensitivity index of construction land was negative in different years, indicating that the increase in construction land would lead to a decrease in the ecosystem service value in the study area; the sensitivity index of unused land was shallow, indicating that the change in unused land had no significant effect on the change in ESV. In addition, the sensitivity indices of land-use types in each period were less than 1, indicating that ESV is inelastic to land-use changes. The study uses ecosystem service value coefficients of different land types suitable for the study area, and the ESV calculation results are reliable.

Unused land 0.000018 0.000018 0.000019 0.000019 0.000019



#### *3.3. The Impact of Land Use Change on Ecosystem Service Value* matrix combined with the ecosystem service value coefficients of different land-use types,

*3.3. The Impact of Land Use Change on Ecosystem Service Value* 

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 13 of 20

The ESV profit and loss matrix (Table 7) is calculated based on the land-use transfer matrix combined with the ecosystem service value coefficients of different land-use types, which can reflect the impact of different land-use changes on the quantity change in ESVs. The land use pattern that led to the loss of ESV was the change from a high ESV land-use type to a low ESV land-use type. In particular, the conversion of woodland to cultivated land and grassland, and cultivated land to construction land contributed the most to the decrease in ESV, reducing the ESV by 24.84 <sup>×</sup> <sup>10</sup><sup>6</sup> Yuan, 22.65 <sup>×</sup> <sup>10</sup><sup>6</sup> Yuan, and 17.20 <sup>×</sup> <sup>10</sup><sup>6</sup> Yuan, respectively. The land-use patterns that contribute to the increase in ESV are grassland and cultivated land converted into woodland (23.75 <sup>×</sup> <sup>10</sup><sup>6</sup> Yuan, 22.64 <sup>×</sup> <sup>10</sup><sup>6</sup> Yuan), and other land types converted into water (43.51 <sup>×</sup> <sup>10</sup><sup>6</sup> Yuan). which can reflect the impact of different land-use changes on the quantity change in ESVs. The land use pattern that led to the loss of ESV was the change from a high ESV land-use type to a low ESV land-use type. In particular, the conversion of woodland to cultivated land and grassland, and cultivated land to construction land contributed the most to the decrease in ESV, reducing the ESV by 24.84 × 106 Yuan, 22.65 × 106 Yuan, and 17.20 × 106 Yuan, respectively. The land-use patterns that contribute to the increase in ESV are grassland and cultivated land converted into woodland (23.75 × 106 Yuan, 22.64 × 106 Yuan), and other land types converted into water (43.51 × 106 Yuan). **Table 7.** ESV profit and loss matrix of central Yunnan urban agglomeration from 2000 to 2020 106 Yuan.

The ESV profit and loss matrix (Table 7) is calculated based on the land-use transfer

**Table 7.** ESV profit and loss matrix of central Yunnan urban agglomeration from 2000 to 2020 10<sup>6</sup> Yuan. **Land Use Type Grassland Cultivated Land Construction Land Woodland Water Area Unused Land**  Grassland 0.00 −10.98 −10.50 23.75 12.99 −0.16


According to previous research results, land-use change has an important impact on ecosystem service value. However, few people have explored the impact of LUCC on ESV. In this study, the Pearson coefficient was used in the Origin2021 software (Origin 2021 edition was developed by OriginLab, USA, and downloaded from https://www. originlab.com/OriginProLearning.aspx, accessed on 25 May 2022) to analyze the relationship between land-use intensity and ESV and its change, and a correlation heat map was drawn (Figure 6) to discuss the impact of land-use intensity on both. The figure shows a significant negative correlation between land-use intensity and the total amount of ESV, and the correlation coefficient was between 0.68 and 0.72. The later the year, the stronger the correlation. There is also a negative correlation between land-use intensity and ESV variation. With the gradual increase in land-use intensity from 2000 to 2020, the influence on the ESV change was also gradually increasing. ecosystem service value. However, few people have explored the impact of LUCC on ESV. In this study, the Pearson coefficient was used in the Origin2021 software (Origin 2021 edition was developed by OriginLab, USA, and downloaded from https://www.originlab.com/OriginProLearning.aspx, accessed on 25 May 2022) to analyze the relationship between land-use intensity and ESV and its change, and a correlation heat map was drawn (Figure 6) to discuss the impact of land-use intensity on both. The figure shows a significant negative correlation between land-use intensity and the total amount of ESV, and the correlation coefficient was between 0.68 and 0.72. The later the year, the stronger the correlation. There is also a negative correlation between land-use intensity and ESV variation. With the gradual increase in land-use intensity from 2000 to 2020, the influence on the ESV change was also gradually increasing.

According to previous research results, land-use change has an important impact on

**Figure 6.** Correlation heatmap. The 2000 L indicates the land-use intensity of the central Yunnan urban agglomeration in 2000; 2000 ESV represents the ecosystem service value in 2000; 2000–2005 denotes the amount of ESV change in central Yunnan urban agglomeration between 2000 and 2005.

### *3.4. Analysis of the Driving Factors of Ecosystem Service Value 3.4. Analysis of the Driving Factors of Ecosystem Service Value*  The Geodetector can explore how different factors explain the value of ecosystem

*Sustainability* **2022**, *14*, x FOR PEER REVIEW 14 of 20

The Geodetector can explore how different factors explain the value of ecosystem services, thereby identifying essential drivers of ESV changes. Kriging interpolation was used in the ArcGIS 10.3 software (Version: 10.3.0.4322; Founder: Environmental Systems Research Institute, Inc. American; Website: https://www.esri.com, accessed on 25 May 2022) to obtain a spatial distribution map of the nine drivers (Figure 7). The results of the factor detection are shown in Table 8. The explanatory strength of each driver in descending order is: land-use intensity (0.532), NDVI (0.497), soil erosion (0.314), slope (0.301), population density (0.29), temperature (0.233), GDP (0.214), elevation (0.21), and precipitation (0.134); the q value of each driving factor was more significant than 0.1, which indicates that these nine factors have an important impact on ESV. However, elevation, temperature, precipitation, soil erosion, and GDP failed the significance test of *p* < 0.05, so these five driving factors were excluded. The final result shows that the explanation of land-use intensity reached 53.2%, which is the most crucial driving factor leading to the change in ESV in the central Yunnan urban agglomeration, indicating that different landuse patterns and human disturbance will have a significant impact on ESV. The explanation of the NDVI to ESV is 49.7%, which indicates that the land cover has an important influence on the total amount and change in ESV. The reason for this is that the ESV coefficients of different vegetation types are significantly different, and the ability to provide ESV is also very different. Slope and population density are also important driving factors of ESV in the central Yunnan urban agglomeration. The greater the population density, the more intense the human activities, and the greater the impact on the ecosystem and its functions. services, thereby identifying essential drivers of ESV changes. Kriging interpolation was used in the ArcGIS 10.3 software (Version: 10.3.0.4322; Founder: Environmental Systems Research Institute, Inc. American; Website: https://www.esri.com, accessed on 25 May 2022) to obtain a spatial distribution map of the nine drivers (Figure 7). The results of the factor detection are shown in Table 8. The explanatory strength of each driver in descending order is: land-use intensity (0.532), NDVI (0.497), soil erosion (0.314), slope (0.301), population density (0.29), temperature (0.233), GDP (0.214), elevation (0.21), and precipitation (0.134); the q value of each driving factor was more significant than 0.1, which indicates that these nine factors have an important impact on ESV. However, elevation, temperature, precipitation, soil erosion, and GDP failed the significance test of *p* < 0.05, so these five driving factors were excluded. The final result shows that the explanation of land-use intensity reached 53.2%, which is the most crucial driving factor leading to the change in ESV in the central Yunnan urban agglomeration, indicating that different landuse patterns and human disturbance will have a significant impact on ESV. The explanation of the NDVI to ESV is 49.7%, which indicates that the land cover has an important influence on the total amount and change in ESV. The reason for this is that the ESV coefficients of different vegetation types are significantly different, and the ability to provide ESV is also very different. Slope and population density are also important driving factors of ESV in the central Yunnan urban agglomeration. The greater the population density, the more intense the human activities, and the greater the impact on the ecosystem and its functions.

**Figure 6.** Correlation heatmap. The 2000 L indicates the land-use intensity of the central Yunnan urban agglomeration in 2000; 2000 ESV represents the ecosystem service value in 2000; 2000–2005 denotes the amount of ESV change in central Yunnan urban agglomeration between 2000 and 2005.

**Figure 7.** Spatial distribution of the main drivers of ESV. (**a**) Elevation; (**b**) Slope; (**c**) Temperature; (**d**) Precipitation; (**e**) NDVI; (**f**) Soil erosion; (**g**) Population density; (**h**) Land-use intensity; (**i**) GDP.


**Table 8.** Factor detection results.

The results of the interaction detection (Table 9) showed that the interaction between factors explained the dependent variable significantly more strongly than individual factors, involving both two-factor enhancement and non-linear enhancement, indicating that the spatial differentiation of ESV in the urban agglomeration of central Yunnan is the result of the synergistic effect of multiple factors. Specifically, the q value of land-use intensity ∩ NDVI has the largest value (0.785), indicating that the interaction between the two has the greatest explanation for ESV. In addition, the interaction types with q value of > 0.7 include precipitation ∩ temperature (0.743), temperature ∩ slope (0.740), NDVI ∩ temperature (0.711), land-use intensity ∩ soil erosion intensity (0.711), land-use intensity ∩ precipitation (0.705), soil erosion ∩ elevation (0.704), GDP ∩ slope (0.703), and GDP ∩ NDVI (0.702). The interaction among the remaining factors also had an enhancing effect on the spatial differentiation of ESV compared with an individual driver alone, especially the factor with a higher q value coupled and coordinated with other factors would amplify the effect on ESV.


**Table 9.** Interactive probe results.

### **4. Discussion**

This research employs land-use dynamic degree, transfer matrix, and hotspot analysis to methodically investigate the spatial and temporal dynamics of land-use and ecosystem service values. In order to support the sustainable development of urban agglomerations in terms of economy, ecology, and resources, we use this foundation to investigate the influence of LUCC on ESV by combining the ESV profit and loss matrix and Pearson correlation coefficient and further exploring other potential factors affecting ESV changes.

### *4.1. Land-Use Change*

The results show that the land-use pattern of the central Yunnan urban agglomeration has changed dramatically in the last 20 years; cultivated land, woodland, and grassland has been significantly reduced, and the construction land, water area, and unused land have increased, which is consistent with the previous research results [47]. People destroyed forests and reclaimed land in the 1990s to pursue higher economic benefits, turning

woodlands and grasslands into cultivated land. This trend was effectively stopped by the policy of returning farmland to forests and grasslands, which was fully implemented in early 2002 as the nation's economic strength and emphasis on environmental protection increased. The urbanization rate of the urban agglomeration in central Yunnan nearly doubled (from 30.62% to 59.90%) between 2000 and 2020 due to the rapid development of urbanization, and the demand for construction land soared, encroaching on a significant portion of high-quality agriculture. In the future, pertinent departments should concentrate on developing policies and procedures to enhance the protection of cultivated land, ensure its quantity and quality, and increase its productivity.

### *4.2. Changes in Ecosystem Service Value*

In this paper, we used the equivalent factor method of Costanza et al. [12], referred to as the equivalent factor table of China's terrestrial ecosystem in 2015 by Xie G.D. et al. [41], combined with the data of the third land and resources survey in Yunnan Province, to calculate the value of ecosystem services in the central Yunnan urban agglomeration for the past 20 years. In the past, the construction land was not taken into account when calculating the total ESV because much research on the evaluation of ecosystem service value assumed that it served no purpose for ecosystem services [48,49]. However, this study believes that urban development has led to a substantial increase in construction land, which has a growing impact on the development pattern of land space and ecological security. On the one hand, certain tourist destinations provide excellent aesthetic and landscaping purposes (the Eiffel Tower, Forbidden City, Sydney Opera House, etc.). However, the production of raw materials, the availability of water resources, and other ecosystem services are all adversely impacted by construction land. Therefore, the ecosystem service function of construction land should not be neglected in ESV evaluation. Additionally, it is discovered that whether the ESV of construction land is considered or not and the selection of the equivalent coefficient for construction land have a significant effect on the estimation results of ESV. Even the ESV estimated using the same land use data may vary if the equivalent coefficients of the construction land used are quite different. Therefore, to enable crossregional comparison, future study should establish a set of standard construction land equivalent coefficients.

### *4.3. Analysis of Driving Factors*

There are many methods used to study the driving factors of ESV, such as the canonical correspondence analysis (CCA) model [50], decomposition analysis [29] and grey correlation analysis [51], etc. However, these methods are not thorough enough to analyze the spatial differentiation of factors. Geodetector, on the other hand, can not only quantitatively detect the interaction between different factors but also judge the significant differences in the influences of different driving factors on the research objects. According to the analysis results of driving factors, it is clear that land-use intensity is the most important driving factor for ESV change, and NDVI, slope, and population density also have significant effects on ESV, which has been verified in previous studies [47,52–55]. The type and quantity of data, the method of dispersion, and the carrier employed by the Geodetector may be factors in why elevation, temperature, precipitation, soil erosion, and GDP failed the significance test.

### *4.4. Innovation, Shortcomings, and Outlook*

In order to examine the ESV changes and driving factors in the central Yunnan urban agglomeration, this study chose it as the research object. It also took into account the spatial heterogeneity of ESV and the interaction between the driving variables. The analysis was more thorough than the study on the impact of a single component on ESV, which helped improve the ecosystem service value. However, there are some restrictions. Firstly, the changed dynamic equivalent coefficients are not employed with the static ESV accounting method [56], which may cause a little inaccuracy in the evaluation of ESV. Secondly, the options for the driving factors are too limited, and the multicollinearity test was not performed beforehand, leaving only four components that passed the significance test. Additionally, the disparities between districts (counties) and years for ESV drivers were not examined. The next phase should use a more precise technique of ESV accounting in conjunction with the study region itself, gathering as much information as possible regarding the driving causes. To suggest tailored ecological protection plans, the driving variables of Geodetector should be studied for various administrative districts in various years [57,58].

### **5. Conclusions**

This research reveals the temporal and spatial characteristics and laws of land use and ESV in the central Yunnan urban agglomeration in 2000, 2005, 2010, 2015, and 2020, then investigates the influence of LUCC on ESV, and lastly, examines the main driving factors of ESV changes using Geodetector, which can provide a valuable reference for optimizing regional territorial spatial patterns and ecological protection and restoration. The results show that woodland is the primary land use type in the urban agglomeration of central Yunnan, accounting for about half of the total area; the land use change is manifested as "three increases and three decreases," that is, cultivated land, woodland, and grassland decreased, while water, construction land, and unused land increased. Among them, the construction land increased the most and changed the fastest. From 2000 to 2020, the total ecosystem service value of the urban agglomeration in central Yunnan decreased by 1.517 billion Yuan. Among the ecosystem types, woodland, grassland, and water area provide more ecosystem services. Among the single ecosystem services, the ESV of climate and hydrological regulation was the highest. The land average ESV showed a spatial distribution characteristic of "high in the west and low in the east." The low-value areas were mainly concentrated in the eastern region and tended to increase with time, and the high-value areas were located in the central region, consistent with the spatial distribution of water. The results of hotspot analysis show that hotspots are mainly distributed in the west and northeast of the urban agglomeration in central Yunnan, while cold spots are concentrated in the central region, and the hotspots decrease, while the cold spots increase. Land-use intensity, NDVI, slope, and population density are the main driving factors of ESV changes. The interaction detection results showed that the effects of pairwise interactions between factors on ESV were significantly enhanced compared with single factors.

In recent years, under the influence of the national urban agglomeration policy, the urban agglomeration in central Yunnan has developed rapidly, and the land use structure has changed, significantly impacting the ecosystem. Therefore, preliminary suggestions are put forward according to the above research results. First of all, from 2000 to 2020, the area of construction land in the urban agglomeration of central Yunnan continued to increase, mainly encroaching on other types of land, primarily cultivated land. This trend is not conducive to food security. Therefore, the disorderly expansion of construction land should be strictly controlled and should strictly abide by the cultivated land "occupation– compensation balance." Secondly, woodland and grassland is important ecological land in the urban agglomeration of central Yunnan, with a large area. In the future, ecological restoration methods should be implemented according to local conditions to strengthen the ecosystem services of woodlands and grasslands. In addition, the unit area of the water area ESV was higher, and it is better to strengthen the comprehensive treatment of water pollution, and ecological and environmental protection and restoration. Finally, it is recommended to tap into the potential of land use, and improve the efficiency of land use by reorganizing and rehabilitating urban villages, "hollow villages," and abandoned industrial and mining land, in addition to adhering to both development and protection, and gradually improving and perfecting the system of paid use of resources and an ecological compensation mechanism. The land use data used in this study are only divided into six categories, eliminating wetlands with high ecosystem services, which will have an effect on the outcomes of ESV accounting. Additionally, it is relatively challenging to obtain data

for multiple driving factors in a long time series, so this study did not explore the driving factors of ESV in different research periods, which had some limitations. The interpretation of remote sensing images will be used in the future to obtain data on different land use types, and attempts will be made to compile a long-term series of data on driving factors for further research.

**Author Contributions:** Conceptualization, writing—review and editing, funding acquisition, F.L.; methodology, software, validation, formal analysis, resources, data curation, writing—original draft, L.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** Supported by the Scientific Research Fund Project of Yunnan Education Department (Grant no. 2021J0592), and Graduate Student Innovation Fund Project of Yunnan University of Finance and Economics (Grant no. 2022YUFEYC098).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

