Spatial-Temporal Dynamic Evaluation of Ecosystem Service Value and Its Driving Mechanisms in China

: Understanding the spatial differentiation and driving mechanisms of ecosystem service value (ESV) is helpful for the protection and sustainable development of the ecological environment. Despite the fact that various studies on ESV have been conducted in various regions, few studies have discussed the spatial differentiation characteristics of ESV in a long time series at a national scale, and even fewer studies have thoroughly examined the driving mechanism of the spatial differentiation of ESV from the perspective of different regions. On the basis of China’s land use data from 1990 to 2018, this paper used the methods of land use dynamics, the ESV evaluation model, hot spot analysis, the barycenter model, and the geographical detector model to study the temporal and spatial differentiation characteristics of land use and ESV in the study area. Moreover, it analyzes the driving mechanisms of the spatial differentiation of ESV at the national scale and in different regions of China. Our results showed the following: (1) Other land types have increased overall, with the exception of grassland. Obvious differences were observed in the single land use dynamics of each land type, especially the construction land, where farmland was the primary source of construction land. With the passage of time, the dynamic degree of comprehensive land use increased. (2) During the study period, ESV generally showed a decreasing trend, with distinct characteristics in high and low ESV areas. The center of gravity of ESV was constantly in Dingxi County and Pingliang City, Shaanxi Province, and its trajectory was generally “S”-shaped. (3) From the perspective of national scale and different regions, the dominant factors affecting the spatial differentiation of ESV were different, and the interaction among multiple factors was signiﬁcantly stronger than that of a single factor. The ﬁndings of the study can provide more scientiﬁc decision-making services for China in order to promote regional environmental protection and develop ecological civilization.


Introduction
Ecosystem services (ES) refer to the life-supporting products and services provided to human beings directly or indirectly by ecosystems under good ecological conditions [1,2]. Ecosystem service value (ESV) is an effective quantitative assessment method which highlights the importance of natural assets for human welfare and translates ecological services into practical applications [3,4]. Numerous factors, including natural and social factors, influence the spatial pattern and evolution of ESV. We can reflect the state of the ecosystem in a deeper level by studying changes in the external spatial pattern and internal influence mechanism of ESV, which is highly important for regional planning and effective protection of the ecosystem [5]. Xinjiang Uygur Autonomous Region was missing, making some data impossible to collect. To interpolate these data, we used data from adjacent years' statistics yearbooks and data from surrounding cities and counties in the same province. spatial differentiations of ESV in 2000 and 2018. The digital elevation model (DEM) [27] and NDVI data were obtained from the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on 1 January 2018). DEM data could be used to extract slope and elevation distribution characteristics within the study area. Meteorological data were acquired from the National Geospatial Data Cloud platform of China (http://www.geodata.cn/data/, accessed on 1 January 2018), which could reflect the regional changes in precipitation and temperature. The population data were obtained from the World Pop Data Center. GDP was obtained from the China Urban Statistical Yearbook. Construction land intensity was calculated with land use data. Road data were provided by OpenStreetMap (which accessed on 1 January 2018, http://www.openstreetmap.org/).

Land Use Dynamic Degree
Land use dynamic degree is an important index for analyzing the change dynamics of land use, including single and comprehensive dynamic degrees [9,28]. To explore the socioeconomic and natural factors driving the change of ESV, we selected nine driving factors, namely, slope, elevation, temperature, precipitation, normalized difference vegetation index (NDVI) [26], population density, gross domestic product (GDP), construction land intensity, and distance from road. Due to the lack of data on some driving factors, this paper only studied the driving mechanism of the temporal and spatial differentiations of ESV in 2000 and 2018. The digital elevation model (DEM) [27] and NDVI data were obtained from the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on 1 January 2022). DEM data could be used to extract slope and elevation distribution characteristics within the study area. Meteorological data were acquired from the National Geospatial Data Cloud platform of China (http://www.geodata.cn/data/, accessed on 1 January 2022), which could reflect the regional changes in precipitation and temperature. The population data were obtained from the World Pop Data Center. GDP was obtained from the China Urban Statistical Yearbook. Construction land intensity was calculated with land use data. Road data were provided by OpenStreetMap (which accessed on 1 January 2022, http://www.openstreetmap.org/).

Land Use Dynamic Degree
Land use dynamic degree is an important index for analyzing the change dynamics of land use, including single and comprehensive dynamic degrees [9,28]. (1) Single dynamic degree (K). Single dynamic degree refers to the change of a certain land use type in a certain period. The calculation formula is as follows: where U a and U b refer to the areas of the specific land use type at the start date and end date (km 2 ), respectively; and T is the study period.
(2) Comprehensive dynamic degree (LC). Comprehensive dynamic degree refers to the change rate of all land use types in a certain period in the study area. The calculation formula is as follows: where U ai and U bi represent the areas of the specific LULC type at the start date and end date, respectively; T is the study period; and n is the number of the land use type.

Construction of Ecosystem Service Value Evaluation Model
Costanza (1997) [6] used the utility and equilibrium value theories to evaluate the global ESV for the first time. Xie (2008) [11] established the ESV equivalent per unit area of the ecosystem in China by combining the structural characteristics of China's ecosystem with the global ESV evaluated by Costanza (1997) [6]. The ESV of cities at the prefecture level and above in China was measured by Xie's (2008) [11] "Equivalent scale of ESV per unit area of terrestrial ecosystem in China". The economic value of the ESV equivalent was about 1/7 of the national average grain yield value. According to previous studies [29], when evaluating the ESV of construction land, various factors must be considered, and the calculation may be difficult, so this study would not estimate it. While estimating the ESV of unused land, the equivalent factor parameters were estimated based on the average of grassland and desert equivalent factors [30].
The economic value of grain output per unit area was calculated by actual grain yields per unit area in China's provinces, autonomous regions, and municipalities directly under the central government. The calculation formula is as follows: where E a represents the economic value of grain yield per unit area in the study area (yuan/hm 2 ), and P refers to the average value of national grain in 2018 (yuan/kg). The average price of grain in 2018 was 2.403 yuan/kg. This price was based on the average value of the national minimum purchase price of early, middle, late indica rice and the national minimum purchase price of wheat from the national development and reform commission's price monitoring center in 2018. Q is the average grain yield per unit area of each study area from 1990 to 2018 (kg/hm 2 ), estimated according to the grain yield data from 31 provinces, autonomous regions, and municipalities directly under the central government in China's statistical yearbooks in 1990, 2000, 2010, and 2018 (Table 1). However, as society develops, the degree of social development in the study area changes with time, and then ESV will also change. This study used 2018 as the base year, and it introduced the revised willingness to pay the coefficient related to social development (D t ) to improve the accuracy of the ESV assessment model [31]. D t is calculated as follows: where l is the social development stage coefficient in relation to actual willingness to pay; E nt refers to the Engel coefficient of each region at year t; E tc and E tr refer to the Engel coefficient of urban and rural regions at year t, respectively; P tc and P tr refer to the proportion of urban and rural population at year t, respectively; k is the stage coefficient of social development; l m refers to the stage coefficient of social development of different cities, and l n refers to the stage coefficient of the social development of China. Combined with the revised willingness to pay the coefficient related to social development (D t ), we constructed an ESV evaluation model suitable for cities at the prefecture level and above in China [31]. The calculation formula is as follows: where ESV is the total value of ecosystem services in the study area (10 10 yuan); A i is the area of class i land use type (km 2 ); M ij is the table value of the equivalent factor corresponding to j ecosystem service function of class i land use type ( Table 2); m is the number of ecosystem types, and n is the number of ecosystem service function items; E a is the economic value of grain output per unit area (yuan/hm 2 ); D t is the coefficient of the social development stage after the introduction of the Engel coefficient and modification by the logistic regression model.

Hot-Spot Analysis
The Getis-Ord G i * can be used to describe the spatial agglomeration degree of ESV. The hot and cold spots represent statistically significant high-value and low-value spatial agglomeration of ESV and its changes [32]. The Getis-Ord G i * index is calculated as [33].
where x j is the ESV of unit j, X is the average value of ESV, W ij refers to the spatial weight coefficient between geospatial units, and n is the total number of units. G i * is represented by the Z-value. For positive z-scores, with statistically significant positive values, the higher the Z-score is, the more intense the clustering of high-values (hot spots) will be. Similarly, for statistically significant negative z-scores, the lower the z-score is, the more intense is the cluster of low-values (cold spots) will be.

Barycenter Model
The barycenter model is an effective tool for analyzing how spatial variables change over time [34]. We introduced the model to reveal the spatial evolution characteristics of ESV. The center of gravity of ESV in China has changed over the study period, and its movement represents the spatial trajectory of ESV change in China. The calculation formula is as follows: where X and Y are the longitude and latitude of the barycenter, respectively; x i and y i are the longitude and latitude of unit i, respectively; w it is the ESV of unit i at year t; and D is the moving distance of the center of gravity.

Geographical Detector Model
Wang (2016) [35] proposed the Geographicol detector model, a statistical tool for detecting spatial differentiation caused by geographical elements and revealing the driving force behind it. Without excessive assumptions, this method can overcome the limitations of traditional statistical analysis methods. This method, which has been widely used in the fields of social economy and environment, can not only detect the spatial differentiation of a single factor but can also assess the interaction of multiple factors [18,36,37]. To detect the driving factors influencing the spatial difference of ESV in the study area, this study primarily used the two modules of factor and interaction detection in the geographic detector model. The calculation formula is as follows: where q is the detection index of the influencing factors of the spatial differentiation of ESV in China, and the interval is [0, 1]. The greater the value of q is, the stronger the heterogeneity of the spatial stratification will be. Conversely, the lesser the value of q is, the randomness of the spatial distribution will be stronger. The value of q is 0, indicating that the factor has no influence on ESV; N and σ 2 are the sample size and variance, respectively; N h and σ 2 h are the secondary sample size and variance, respectively; and h = 1, 2, . . . , L, which is the stratification of variable Y or factor X.

Change in Land Use from 1990 to 2018
On the basis of the land use data of four periods in 1990, 2000, 2010, and 2018, the land use transfer matrix in China was calculated, and the transfer chord diagram between different land use types was created ( Figure 2). During the study periods of 1990-2000, 2000-2010, 2010-2018, and 1990-2018, the evolution characteristics of various land use types in the study area were analyzed. The results showed that throughout each time, the area of each land type in the study area changed to varying degrees. Woodland and grassland were always the two primary land types. Specifically, the change in grassland net area always showed a decreasing trend. The net area declined by 12.71% across the study period. Especially from 2010 to 2018, the net area of grassland dropped by 11.4% and mostly shifted to farmland, woodland, and unused land; Over the study period, the change in woodland was relatively moderate, with a steady decreasing trend from 1990 to 2000, followed by a slowly increasing trend from 2000 to 2018; From 1990 to 2000, the net area of farmland increased by 1.61%, whereas from 2000 to 2018, it began to slowly decrease, with most of the land being moved to woodland, grassland, and construction land; From 1990 to 2010, the unused land and water body changed steadily, whereas the increase from 2010 to 2018 showed an increasing trend, with increases of 11.39% and 6.07%, respectively; The change in construction land was the most significant, and the change rate continued to increase throughout the three study periods. The change rate of construction land reached 34.95%, particularly from 2010 to 2018.
According to the land use data of each year, Formulas (1) and (2) Tables 3 and 4. From the perspective of single dynamic degree, the land use dynamic degree of grassland was always negative, with the greatest decline range (1.17%) between 2010 and 2018. Farmland and woodland change rates were relatively slow, but overall showed an increasing trend, with increases of 1.50 × 10 4 km 2 and 2.33 × 10 4 km 2 , respectively. The land use dynamic degree of farmland showed an increasing trend from 1990 to 2000, and then decreased by a tiny margin in 2000-2010 and 2010-2018, and the decreasing range continued to increase. The land use dynamic degree of woodland was was always less than 0.1% in each period, and the change range was not obcious; The water body kept increasing throughout the study period. Especially from 2000 to 2010, The land use dynamic degree of water body reached 0.53%; The unused land changed in a slowly decreasing trend from 1990 to 2010, and then it increased slightly from 2010 to 2018; During the study period, the area of construction land continued to grow rapidly. From 1990 to 2018, the dynamic degree of land use in the whole study period from 1990 to 2018 increased by 2.50%, with a total increase of 10.72 × 10 4 km 2 . This discovery was mostly due to population growth, urban development, and other factors, which resulted in a progressive rise in the amount of construction land. The comprehensive land use dynamics were always less than 0.05% from 1990 to 2010, but this changed greatly from 2010 to 2018, reaching 0.36%.  Tables 3 and 4. From the perspective of single dynamic degree, the land use dynamic degree of grassland was always negative, with the greatest decline range (1.17%) between 2010 and 2018. Farmland and woodland change rates were relatively slow, but overall showed an increasing trend, with increases of 1.50 × 10 4 km 2 and 2.33 × 10 4 km 2 , respectively. The land use dynamic degree of farmland showed an increasing trend from 1990 to 2000, and then decreased by a tiny margin in 2000-2010 and 2010-2018, and the decreasing range continued to increase. The land use dynamic degree of woodland was was always less than 0.1% in each period, and the   During the study period, the total amount of ESV decreased by 6.67%, with an average yearly decline rate of 0.24%. Figure 3 shows how the ESV produced by various land use types and ecosystem subservice structure alterations varied. From the perspective of various land use types, the ESV produced by woodland was the largest, followed by grassland, farmland, unused land, and water body. The ESV generated by farmland and unused land increased by 3.50% and 7.41%, respectively, whereas the ESV of other land types decreased, especially grassland, which decreased by a total of 23.80%. From the perspective of various ecosystem subservice structures, the ESV generated by the four ecosystem service structures was in the following order: regulating services > supporting services > provisioning services > cultural services. The change trend was consistent across the three research periods, showing a trend of initially decreasing, then increasing, and finally decreasing. During the whole study period, supporting services saw the greatest decrease in ESV, with a total decrease of 7.35%, whereas provisioning services experienced the least decrease, with a total decrease of 4.79%.  ESV was divided into five levels ( Figure 4) using the natural breakpoint method in ArcGIS 10.2: extremely low ESV (<3.48 × 10 10 yuan), low ESV (3.48 × 10 10 -8.27 × 10 10 yuan), medium ESV (8.47 × 10 10 -16.68 × 10 10 yuan), high ESV (16.68 × 10 10 -45.77 × 10 10 yuan), and extremely high ESV (>45.77 × 10 10 yuan). In the four periods, the spatial distribution of ESV was largely similar. The high and extremely high ESV were primarily distributed in Tibet, Xinjiang, and Neimenggu, due to the huge distribution of grassland in the western region and the large ESV produced by grassland. The medium, low, and extremely low ESV were primarily distributed in the central, southern, and eastern coastal regions due to the large amount of farmland in these areas which generated less ESV than woodland and grassland.

Evolution Characteristics of Cold and Hot Spots of Ecosystem Service Value
With the help of the hot spot analysis tool in ArcGIS 10.2, this study conducted a hot ESV was divided into five levels ( Figure 4) using the natural breakpoint method in ArcGIS 10.2: extremely low ESV (<3.48 × 10 10 yuan), low ESV (3.48 × 10 10 -8.27 × 10 10 yuan), medium ESV (8.47 × 10 10 -16.68 × 10 10 yuan), high ESV (16.68 × 10 10 -45.77 × 10 10 yuan), and extremely high ESV (>45.77 × 10 10 yuan). In the four periods, the spatial distribution of ESV was largely similar. The high and extremely high ESV were primarily distributed in Tibet, Xinjiang, and Neimenggu, due to the huge distribution of grassland in the western region and the large ESV produced by grassland. The medium, low, and extremely low ESV were primarily distributed in the central, southern, and eastern coastal regions due to the large amount of farmland in these areas which generated less ESV than woodland and grassland.
Land 2022, 11, x FOR PEER REVIEW 11 of 20 decreased, with the majority of them concentrated in Mudanjiang City. The spatial distribution of hot spots was concentrated mostly in the central and southwest areas. The number of sub hot spots slightly increased compared with the previous stage, and was primarily distributed in the central area.

Evolution Characteristics of Cold and Hot Spots of Ecosystem Service Value
With the help of the hot spot analysis tool in ArcGIS 10.2, this study conducted a hot spot analysis on the ESV in 1990, 2000, 2010, and 2018 ( Figure 5). In the four periods, the distribution of cold spots and hot spots in the study area was relatively consistent, and both were relatively concentrated. The hot spots were primarily distributed in the whole region of Tibet, parts of Xinjiang, Heihe City and Tahe County in northeast Heilongjiang Province, and Hailar City in Neimenggu. The cold spots were primarily distributed in the eastern region. The sub-hot and sub-cold spots were sporadically scattered. The variations in cold and hot spots were less noticeable with time. The number of sub-cold spots has increased over time, but they were still mostly distributed around the cold spots. The distribution of sub-hot spots varied with time, and the distribution areas were different.  To further explore the temporal and spatial distribution characteristics of ESV changes in China from 1990 to 2018, the cold and hot spot map in the change of ESV in the study area was described again with the use of a hot spot analysis tool ( Figure 6). From 1990 to 2000, the cold spots were primarily distributed in Tibet, the sub-cold spots were distributed in Yuxi City, the hot spots were primarily distributed in north China, and the sub-hot spots were widely distributed in north China, the eastern coast, and the southwest; From 2000 to 2010, the cold spots were primarily distributed in the south, whereas the sub-cold spots increased slightly, primarily distributed around the cold spot area. In comparison with the stage from 1990 to 2000, the number of hot spots increased significantly, and the degree of aggregation also increased, mostly distributed in Qinghai, Sichuan, and Heilongjiang Province. The number of sub-hot spots decreased significantly and scattered around the hot spots; From 2010 to 2018, the regional distribution of cold and hot spots changed significantly. Cold spots were mainly distributed in Tibet and some areas of Heilongjiang and Inner Mongolia in the northeast. The number of sub cold spots decreased, with the majority of them concentrated in Mudanjiang City. The spatial distribution of hot spots was concentrated mostly in the central and southwest areas. The number of sub hot spots slightly increased compared with the previous stage, and was primarily distributed in the central area.

Barycenter Evolution of Ecosystem Service Value
According to the barycenter model in Formulas (6) and (7), the change in China's ESV center of gravity trajectory was calculated. As shown in Figure 7, the center of gravity of ESV was primarily distributed in Pingliang City and Dingxi County of Shaanxi Province from 1990 to 2018, and the overall change difference was obvious, with an "S"-shaped trend.

Barycenter Evolution of Ecosystem Service Value
According to the barycenter model in Formulas (6) and (7), the change in China's E center of gravity trajectory was calculated. As shown in Figure 7, the center of gravit ESV was primarily distributed in Pingliang City and Dingxi County of Shaanxi Provi from 1990 to 2018, and the overall change difference was obvious, with an "S"-sha trend.

Driving Factors of Regional Differences in Ecosystem Service Value
Taking 2000 and 2018 as examples, a total of nine driving factors were selected fr the two aspects of physical geography and social economy. The ESV was used as the pendent variable of the geographical detector. The independent variables were the sl (X1), elevation (X2), average annual temperature (X3), annual precipitation (X4), ND

Driving Factors of Regional Differences in Ecosystem Service Value
Taking 2000 and 2018 as examples, a total of nine driving factors were selected from the two aspects of physical geography and social economy. The ESV was used as the dependent variable of the geographical detector. The independent variables were the slope (X1), elevation (X2), average annual temperature (X3), annual precipitation (X4), NDVI (X5), population density (X6), GDP (X7), construction land intensity (X8), and distance from road (X9). The influence degree of each driving factor on the spatial differentiation characteristics of ESV in the study area, as well as the interaction characteristics between the driving factors, were described using the geographical detector's factor detection and interactive detection functions. China spanned a huge area from east to west and north to south, and different natural and socioeconomic factors influenced different regions in various ways. This paper divided the study area into six regions (i.e., northeastern region, north China, east China, south central region, northwest region, and southwest region) for further analysis to better detect the effects of various driving factors on the temporal and spatial differentiation of ESV in the study area. As shown in Table 5, the influence of each driving factor on the spatial differentiation characteristics of ESV in the study area was mostly explained using the Q statistic. From a national scale, the p value of all the driving factors, which were all significant characteristics, was less than 0.01 during the study period. The main natural factors that influenced the spatial differentiation of ESV were elevation (X2) and average annual temperature (X3). The main social factors that influenced the spatial differentiation of ESV were always population density (X6) and distance from road (X9). Simultaneously, the influence of construction land intensity (X8) on the spatial differentiation of ESV changed significantly from 0.19 to 0.33 between 2000 and 2018. The effects of other driving factors on the spatial differentiation of ESV did not change significantly over time. This result was mostly due to the wide latitude and longitude spans of the various regions in China, as well as the fast social and economic growth with the passage of time. Moreover, construction land encroached on other ecological lands, making construction land intensity (X8) one of the main driving factors of ESV spatial differentiation. From the regional scale, the main driving factors influencing the spatial differentiation characteristics of ESV were different. The most important socioeconomic factors for the spatial differentiation of ESV in the northeast region and north China were population density (X6), construction land intensity (X8), and distance from road (X9). The annual average temperature (X3) had a significant influence on the spatial differentiation of ESV in the northeast region in 2000, but not in 2018. The Q statistic of the annual average temperature (X3) on ESV In north China between 2000 and 2018 was 0.74, indicating significant characteristics; In east China and the south central region, slope (X1), population density (X6), and construction land intensity (X8) had the strongest explanations. The influence of NDVI (X5) increased over time, whereas that of distance from road (X9) on the spatial differentiation of ESV in the two regions was insignificant; In the northwest region, the main driving factors were NDVI (X5) and distance from road (X9), with the annual precipitation (X4) and construction land intensity (X8) gradually increasing their explanatory strength for ESV over time; Natural factors such as elevation (X2), average yearly temperature (X3), and NDVI (X5) exhibited significant influence on the ESV of the southwest region, in addition to distance from road (X9), which explained the largest degree of spatial differentiation of ESV.
The interactive detection results showed that bi-enhanced or nonlinear enhanced effects were primarily interaction types between natural and socioeconomic factors. This finding indicates that the spatial differentiation of ESV in the study area was caused by the interaction of multiple driving factors, rather than being purely influenced by one factor. From a national scale (Figure 8), the interaction degree between slope (X1) and distance from road (X9) was the largest, with an interaction intensity of 0.74. The influence of the interaction between distance from road (X9) and other factors on the spatial differentiation of ESV was greater than 0.60, and the interaction between distance from road(X9) and natural factors was significantly greater than the interaction between distance from road Land 2022, 11, 1000 14 of 20 (X9) and socioeconomic factors. From the regional scale, we listed the top three dominant interactions of influence of different driving factors on the spatial differentiation of ESV in each region (Table 6). Overall, construction land intensity (X8) and distance from road (X9) had the most significant interactions with other factors, especially with natural factors. Thus, the interactive effect between natural and socioeconomic factors further enhanced the effect on the spatial differentiation of ESV. Specifically, in the northeast, northwest, and southwest regions, the main interactions came from the interaction between natural factors and distance from road (X9) in 2000 and 2018. The interaction between elevation (X2) and distance from road (X9) in the northeast region was always the largest; The interaction between distance from road (X9) and altitude (X2) and the interaction between distance from road (X9) and annual average temperature (X3) in Northwest China always ranked in the top two. In the southwest region, the interactions between distance from road (X9) and natural factors were above 0.90, and the interaction intensity was relatively large. The primary interactions in north China, east China, and the south central region were produced by the interaction of natural factors and the intensity of development land (X8). During the research period, the top three interactions in north China were primarily the interaction intensity among construction land intensity (X8), with elevation (X2), annual average temperature (X3), and NDVI (X5), all of which had interaction strengths over 0.92; In east China and the south central region, the interaction between slope (X1) and construction land intensity (X8) was the strongest, and the interaction between them kept growing as the years passed. of the interaction between distance from road (X9) and other factors on the spatial differentiation of ESV was greater than 0.60, and the interaction between distance from road(X9) and natural factors was significantly greater than the interaction between distance from road (X9) and socioeconomic factors. From the regional scale, we listed the top three dominant interactions of influence of different driving factors on the spatial differentiation of ESV in each region (Table 6). Overall, construction land intensity (X8) and distance from road (X9) had the most significant interactions with other factors, especially with natural factors. Thus, the interactive effect between natural and socioeconomic factors further enhanced the effect on the spatial differentiation of ESV. Specifically, in the northeast, northwest, and southwest regions, the main interactions came from the interaction between natural factors and distance from road (X9) in 2000 and 2018. The interaction between elevation (X2) and distance from road (X9) in the northeast region was always the largest; The interaction between distance from road (X9) and altitude (X2) and the interaction between distance from road (X9) and annual average temperature (X3) in Northwest China always ranked in the top two. In the southwest region, the interactions between distance from road (X9) and natural factors were above 0.90, and the interaction intensity was relatively large. The primary interactions in north China, east China, and the south central region were produced by the interaction of natural factors and the intensity of development land (X8). During the research period, the top three interactions in north China were primarily the interaction intensity among construction land intensity (X8), with elevation (X2), annual average temperature (X3), and NDVI (X5), all of which had interaction strengths over 0.92; In east China and the south central region, the interaction between slope (X1) and construction land intensity (X8) was the strongest, and the interaction between them kept growing as the years passed.

ESV Dynamic Change in Response to LUCC in China
The ESV of each Chinese region was calculated using Xie's (2008) [11] ESV equivalent per unit area of the ecosystem. Simultaneously, the sensitivity index was used to reflect the change of ESV caused by the change in ESV coefficient by 1% (See References [38,39] for details). The results showed that the ESV has decreased by 1.57 trillion yuan, with an average annual decline rate of 0.24%. During the study period, the sensitivity indexes of different land use types have been in a stable state and were less than 1 (Figure 9), indicating that the ESV coefficient was reasonable in the study area. In addition, the sensitivity index of different land use types in the study area was in the following order from high to low: woodland > grassland > farmland > unused land > water body. Woodland and grassland contributed the most ESV, whereas unused land and water body contributed the least. Therefore, to ensure the stability and coordination of ESV, we should improve the protection of woodland and grassland. From the perspective of the spatial distribution of ESV, we found an obvious phenomenon of high and low value aggregation of ESV. The high and extremely high ESV was mainly distributed in the northwest region with more grassland. The low ESV was mainly distributed in the central region and the eastern coastal region; this part of the region had rapid economic development and serious occupation of ecological resources. The center of gravity of ESV in the study area was always distributed in Shaanxi Province, and continued to move north by east with the passage of time, generally in an "S"-shaped trend. The ESV was affected by changes in land use. To maintain the security and stability of ecosystem service capacity, relevant departments should properly adjust the land use structure and strengthen the ecological land occupation and compensation mechanism in each region.

Impact Factors on ESV Distribution
The spatial differentiation of ESV was affected to some extent by the interaction between human activities and the ecological environment [40]. Exploring the impact mechanism of the spatial differentiation of ESV, which was of great significance to China's ecosystem management and ecological pattern optimization, could help in further explaining the generating mechanism of ecological problems.
This study quantitatively analyzed the relative importance of the driving factors to the ESV and the degree of interaction among driving factors in the national and regional scales. The results showed that the primary driving factors affecting the spatial differentiation of ESV were different in the whole nation and different regions. From a national scale, elevation (X2), average annual temperature (X3), population density (X6), and distance from road (X9), were the primary driving factors affecting the spatial differentiation of ESV, especially distance from road (X9). From the different regional scales, the driving factors affecting the spatial differentiation of ESV in the northeast, north, east, and south central regions were mainly socioeconomic factors, such as population density (X6) and distance from road (X9). In northwest and southwest regions, the primary driving factor was not only the distance from road (X9), but also the annual average temperature (X3), NDVI (X5), and other natural factors. At the same time, the influence of construction land intensity (X8) on the spatial differentiation of ESV also increased over the years, and it has recently become one of the major driving factors affecting the spatial differentiation of ESV in the nation scale and different regional scales. The interaction between driving factors had a greater influence on ESV than that of a single factor, and the interaction between socioeconomic and natural factors had a greater influence on the spatial differentiation of ESV. In particular, the interaction among construction land intensity (X8), distance from road (X9), and other natural factors was always stronger than that with socioeconomic factors. This finding was made because, as the social economy develops, the intensity of human activities has increased, the speed of urban expansion also increases, and some irrational reclamation and construction has occupied green vegetation land, such as farmland, woodland, and grassland, thereby putting great pressure on the ecosystem's stability [41]. Therefore, China should properly control the rate of urbanization and optimize the land use structure. While focusing on the coordinated development of ecological and economic construction, more ecological protection and restoration projects must be implemented to increase vegetation coverage and mitigate water and soil loss, thereby ensuring the stability of the ecosystem [42,43].
The ESV of each Chinese region was calculated using Xie's (2008) [11] ESV per unit area of the ecosystem. Simultaneously, the sensitivity index was us the change of ESV caused by the change in ESV coefficient by 1% (See Refere for details). The results showed that the ESV has decreased by 1.57 trillion yu average annual decline rate of 0.24%. During the study period, the sensitivit different land use types have been in a stable state and were less than 1 (Fig  cating that the ESV coefficient was reasonable in the study area. In addition, th index of different land use types in the study area was in the following orde to low: woodland > grassland > farmland > unused land > water body. Wo grassland contributed the most ESV, whereas unused land and water body the least. Therefore, to ensure the stability and coordination of ESV, we shou the protection of woodland and grassland. From the perspective of the spatial of ESV, we found an obvious phenomenon of high and low value aggregation high and extremely high ESV was mainly distributed in the northwest region grassland. The low ESV was mainly distributed in the central region and coastal region; this part of the region had rapid economic development and s pation of ecological resources. The center of gravity of ESV in the study area distributed in Shaanxi Province, and continued to move north by east with th time, generally in an "S"-shaped trend. The ESV was affected by changes in l maintain the security and stability of ecosystem service capacity, relevant d should properly adjust the land use structure and strengthen the ecological l tion and compensation mechanism in each region.

Limitations and Improvement
Xie (2008) [11] created the ESV equivalent per unit area of the ecosystem in China, which has been widely used by Chinese scholars in the assessment of ESV [44,45]. The revised willingness to pay coefficient related to social development (D t ) was introduced in our paper based on this table to modify the evaluation model of ESV. However, due to the large span from north to south and east to west, significant differences were observed in the geographical environments of the various regions in the study area. This paper did not formulate the ESV equivalent per unit area of terrestrial ecosystem suitable for each region according to the actual situation of each region, resulting in a certain deviation in calculation. Furthermore, while using the geographical detector model to analyze the driving mechanism of ESV spatial difference, this study did not consider the relevant policy factors. At the same time, due to the limitations of the model's own characteristics, it could not accurately describe the spatial characteristics of each driving factor's influence on ESV. In the future, we will select appropriate policy factors and use a geographically and temporally weighted regression model to advance our arguments.

Conclusions
(1) During the study period, grassland always showed a decreasing trend, with a total decline of 38.74 × 10 4 km 2 , due to the large quantity of grassland transferred to farmland, woodland, and unused land. Farmland increased from 1990 to 2000, then slightly decreased from 2000 to 2018, but it kept growing throughout the study period, with a net area increase of 1.5 × 10 4 km 2 . The change in woodland in each research stage was not obvious. The land use dynamic degree of the water body was always positive, increasing by 2.06 × 10 4 km 2 throughout the study period. From 1990 to 2010, the amount of unused land decreased, but increased from 2010 to 2018, with an overall increase. Construction land changed the most throughout the study period, and the dynamic degree of construction land was always positive, with an increase of 10.72 × 10 4 km 2 , and farmland was the largest source of growth of construction land. The comprehensive dynamic degree of land use has also increased over time, especially from 2010 to 2018.
(2) From 1990 to 2018, the high ESV was mostly in the western region. The total amount of ESV decreased by 6.67% over the last three decades, with an average yearly decline rate of 0.24%. From the perspective of different land use types, the ESV generated by woodland and grassland was the main part of the total ESV. During the study period, farmland and unused land showed an increasing trend, whereas other land use showed a decreasing trend, especially grassland. From the perspective of different ecosystem subservice structures, regulating services produced the most ESV, whereas cultural services produced the least. Supporting services decreased the most during the study period, with a total decrease of 7.35%. The distribution of ESV cold and hot spots in the study area was relatively consistent in 1990, 2000, 2010, and 2018. The hot spots were mostly concentrated in the northeast and northwest regions, whereas the cold spots were mostly concentrated in east China. During the study period, the numbers of both remained relatively consistent. The number of sub-cold spots increased with time, mainly in the vicinity of cold spots. The distribution of sub-hot spots changed with time. The center of gravity model showed that ESV's center of gravity was always distributed in Shaanxi Province, and it continues to migrate northeast of China in an "S"-shaped trend.
(3) The results of the Geographical detector model revealed that the spatial difference of ESV in China was driven by a combination of natural and socioeconomic factors, and that the influence of different driving factors on ESV in the study area was significantly different. From a national scale, the ESV is mostly affected by the distance from the road (X9). From the regional scale, the primary factors affecting the spatial differentiation of ESV were different. Socioeconomic factors dominated in the northeast area, north China, east China, and the central south region. Distance from road (X9) was the most important factor in the northwest and southwest areas (X9). In southwest China, natural factors such as annual average temperature (X3) and NDVI (X5) were also primary factors. At the same time, the spatial differentiation of ESV was the result of the interaction of multiple driving factors rather than a single factor. In particular, the interaction among construction land intensity (X8), distance from road (X9), and other driving factors was significantly enhanced.
(4) By analyzing the spatiotemporal differentiation of China's ESV and exploring the driving mechanism of ESV at national and regional scales, this paper proposed some policy implementation suggestions for environmental conservation in China. First, we should pay more attention to the areas of low ESV, especially in the central and eastern coastal regions, where fast economic growth and high degree of urbanization have transformed a huge quantity of ecological land into construction land. For these areas, some ecological land for specific applications should be earmarked, and protective measures should be taken to avoid being occupied by construction land. Simultaneously, land consolidation for residential areas and related infrastructure in rural areas should be promoted. Second, some areas with high ESV, especially in northwest China, have more woodlands and grasslands, which produce more ESV. However, the effects of human activities on these areas continues to increase over time. Indiscriminate land reclamation and overgrazing have led to the destruction of woodland and the degradation of grasslands. For these areas, to limit unreasonable human activities, various ecological protection measures should be implemented. For example, policy subsidies should be offered to relevant employees to encourage them to return farmland to forest and grassland.