Quantifying the Temporal and Spatial Patterns of Ecosystem Services and Exploring the Spatial Differentiation of Driving Factors: a Case Study of Sichuan Basin, China

1 Introduction

Ecosystem services (ES) are ecological characteristics, functions, or processes that contribute directly or indirectly to human well-being (Costanza et al., 1997; Yang et al., 2021a; Yu et al., 2022b). The United Nations Millennium Ecosystem Assessment divides ES into four categories: support services, regulation services, supply services, and cultural services. There are complex relationships between different services (Tallis and Kareiva, 2005; Luo et al., 2021). The concept of ES as a “bridge” between natural ecosystems and human social systems serves as a foundation for improving human well-being and achieving regional sustainable development. The integrated assessment of ES is an effective hybrid that has become a frontier and research hotspot in geography, ecology, and environmental science, and it is also central to the sustainable development of ecosystems and the management of ecological resources (Costanza et al., 2017; Castillo-Eguskitza et al., 2018; Yu et al., 2021; Chen et al., 2022; Hu et al., 2022).

It is generally accepted that ES can be influenced by both natural and human factors (Burkhard et al., 2012; Burkhard et al., 2015; Zhang et al., 2020; Yang et al., 2021b). First, climate change can directly change the capacity to supply services such as food production, freshwater supply, and net primary productivity by changing ecosystem properties. In addition to climate change, elevation, slope, soil, and vegetation also contribute to changes in ES (Wolff et al., 2015; Yang et al., 2020). Human activities have both direct and indirect impacts on the ecosystem, with both beneficial and detrimental impacts. For example, the human activities of urban construction, forest, and grassland into an impermeable surface, resulting in the decline of the ability to provide ES. Payments for Environmental Services (PES) practices around the world, such as China’s “conversion of cropland to forest and grassland” and the South African government’s Working for Water program, have improved soil conservation (SC) and water yield (WY) in the region (Wei et al., 2017; Wu et al., 2019; Bing et al., 2021).

Therefore, we not only need to fully understand the spatiotemporal patterns of ES but also need to explore the driving factors of ES to develop scientifically sound ES management approaches (Geijzendorffer et al., 2015; Wei et al., 2017). Several recent studies have shown this concern, including Chen et al. (2020) the need to consider different drivers and spatial spillovers in the formulation of strategies for integrated ecosystem management and sustainable land use in urban agglomerations (Chen et al., 2020). Peng et al. (2020) discussed the impact of urbanization on the ES budget through decoupling analysis (Peng et al., 2020). Hu et al. (2021) took the Shanxi province of China as an example; this study explored the spatial differentiation of five ES and captured the main driving factors of ES change from a geospatial perspective (Hu et al., 2021).

Existing studies have shown that ES often has certain spatial heterogeneity (Li et al., 2013; Song and Deng, 2017; Zhou et al., 2018; Licheng et al., 2019; Yang et al., 2019), which will be more obvious in regions with high spatial heterogeneity in geography, society, and economy. However, most of these studies are based on global regression, mainly considering a single driving factor, without considering the comprehensive effect. Less attention has been paid to the spatial non-stationary relationship between the leading factors affecting ES. Commonly used methods such as correlation analysis and overlay analysis can only clarify the relationship between individual factors and ES and cannot quantitatively derive significant impact factors and their intensity, nor can it reveal the spatial heterogeneity of factor interaction and the possible synergistic or antagonistic effects among the factors (Jiang et al., 2022), especially when exploring multivariate driver analysis, the problem of multicollinearity between drivers often occurs (Toutenburg, 2006; Yu et al., 2022a). When there is a multicollinearity problem, the variance of the parameter estimator will be too large, the accuracy will be reduced, the significance test of the variable will be meaningless, and the influence of the explanatory variable on the explained variable cannot be correctly judged, resulting in the unreasonable value of the parameter estimator (Yu et al., 2009). Therefore, based on the PCA-MGWR method, this study combines the advantages of principal component analysis (PCA) and MGWR. Before MGWR analysis, the PCA of independent variables is carried out, and several uncorrelated principal components extracted are used as independent input variables for MGWR to analyze the spatial impact distribution of each principal component. PCA-MGWR can not only eliminate the multicollinearity problem caused by the correlation of independent variables but also explain the spatial effects of different ecosystem service factors.

To sum up, although natural-human factors have a significant impact on ES, the current research on their comprehensive consideration is still insufficient. Therefore, this study hopes to make up for the current gap. Sichuan Basin is the ecological barrier to the upper reaches of the Yangtze River. The basin has a good natural resource endowment and is one of the regions with the largest supply of ES in China (Pan and Wang, 2021). Fully understanding the spatio-temporal pattern of ES in the Sichuan Basin, and exploring the natural-human multiple impact factors, is of great importance to ensure the ecological security of the upper reaches of the Yangtze River and even the ecological security of the country.

At present, more and more attention is focused on exploring the driving force of the impact on the ES (Zhong et al., 2019; Sannigrahi et al., 2020). Relevant studies have focused on the spatial-temporal pattern of ES in Sichuan Province and its various driving factors from a statistical perspective (Liu et al., 2020; Niu et al., 2022). However, these studies generally focus on the analysis of change patterns, paying less attention to the spatial change characteristics of driving forces, and often ignoring the scale differences of different driving factors. Therefore, some studies have used the spatial regression model to test the spatial non-stationary relationship between the four driving forces on SC and WY in Sichuan Province (Huang et al., 2022). This study not only adds the measurement indicators and driving factors of ES in the Sichuan Basin but also uses the MGWR model to explore the nature and intensity of the impact of various factors on ES in the Sichuan Basin. In view of the limitations and deficiencies of previous studies, this study has the following specific objectives: (I) evaluate the ES and their temporal and spatial changes in the Sichuan Basin; (II) quantitatively assess the driving factors of temporal and spatial changes of ES in Sichuan Basin; and (III) analyze the comprehensive response and spatial differentiation of various driving factors of ES in Sichuan Basin.

2 Methodology

2.1 Study Area

Sichuan Basin lies at the intersection of the Belt and Road and the Yangtze Economic Belt (27° 68′∼ 32° 30′N, 101° 95′∼ 108° 95 E, Figure 1). Influenced by the subtropical monsoon climate, the annual average temperature is 13–20°C, and the annual average precipitation is 1,000 mm. The scope of the study includes 15 cities, including Chengdu in Sichuan and 29 districts and counties in Chongqing, with a total area of about 185,000 km². There are many types of landforms in the region, forming various kinds of vertical climatic zones; the climatic differentiation is remarkable, the natural ecological environment is diverse, and the ecological resources are rich. In 2019, the Chengdu-Chongqing urban agglomeration had a resident population of 96 million and a regional GDP of nearly 6.3 trillion yuan. The 2021 outline for constructing the Double-City Economic Circle proposes that the Chengdu-Chongqing urban agglomeration should play an exemplary role in promoting the green development of the Yangtze River Economic Belt and the ecological protection of western China.

FIGURE 1. Location and administrative division of the Sichuan Basin.

2.2 Data

This study refers to the methods of previous studies and the data requirements of the InVEST model to quantify the ES (Xie and Ng, 2013; Chen et al., 2020; Zhang X. et al., 2021; Huang et al., 2022). SC is an important indicator to measure the water and soil conservation of the ecosystem. Carbon storage (CS) can directly reflect the carbon absorption and storage capacity of the ecosystem. WY is the ability of the ecosystem to retain rainwater under the joint action of plants and soil. At the same time, considering that Sichuan Basin is an important water source and protected area of the Yangtze River and the Yellow River, soil erosion is an important factor in its flood disaster. Sichuan Basin is also rich in natural resources, including forests, meadows, wetlands, and snow mountains. It is a strategic area for ecological security in Western China. Therefore, combined with the natural landform and the regional characteristics of the social economy in the Sichuan Basin, these three ecosystems can be regarded as important ecosystems for the sustainable development of resources and the environment. It involves the use of multi-source data such as land use, digital elevation model, meteorological data, soil data, watershed boundary, and social and economic data. to comprehensively evaluate the ES in the Sichuan Basin. The data sources and their descriptions are shown in Table 1.

TABLE 1. Data source and its detailed description.

2.3 Quantification of Ecosystem Services

2.3.1 Carbon Storage (Carbon Storage)

Carbon sequestration services, also known as “Carbon Storage” services, are important regulatory services in ecosystems. The “carbon” module in the InVEST model is used to evaluate the CS in Sichuan Basin as the supply of carbon sequestration services. The input carbon storage data refer to previous local studies (Liu et al., 2004; Tang et al., 2018; Liang et al., 2021). The formula is as follows:

where $C_{tot}$ is the total carbon storage $\left(t\cdot h{m}^{-2}\right)$, $C_{above}$ is aboveground biological carbon $\left(t\cdot h{m}^{-2}\right)$, $C_{below}$ is underground biological carbon $\left(t\cdot h{m}^{-2}\right)$, $C_{soil}$ is soil organic carbon $\left(t\cdot h{m}^{-2}\right)$, and $C_{dead}$ is dead organic matter. The details of each coefficient are shown in Supplementary Table SA1 of the supplementary data.

2.3.2 Water Yield (Water Yield)

WY refers to the ability of ES to store water resources. The Sichuan Basin is located in the middle reaches of the Yangtze River, so it is vital to assess the WY for the storage and utilization of water resources. The “water yield” module is used to quantify the supply of WY in the Sichuan Basin. This module is defined as the part of precipitation minus evapotranspiration based on the principle of water balance (Sharp et al., 2014; Sun et al., 2018). The formula is as follows:

where $S_{wp}$ is the annual water yield (mm); $AE{T}_{\left(x\right)}$ is the actual annual evapotranspiration (mm) of the grid unit; $P_x$ is the annual precipitation (mm) of the grid unit; $PE{T}_{\left(x\right)}$ is the potential evapotranspiration of the grid unit; $K{c}_{\left(x\right)}$ is the crop evapotranspiration coefficient; $E{T}_{o\left(x\right)}$ is the reference crop scatter; $AW{C}_{\left(x\right)}$ is the plant available water content; $W_{\left(x\right)}$ is the empirical parameter; and Z is the Zhang coefficient (Zhang et al., 2004). The details of each coefficient are shown in Supplementary Table SA2 of the supplementary data.

2.3.3 Soil Conservation (Soil Conservation)

In this study, the modified universal soil loss equation (USLE) is used to estimate regional soil conservation and soil erosion. In sediment retention, the model takes into account the important hydrological process of plant-induced soil erosion and sediment interception, which makes the simulation results more scientific and reasonable. We used the “sediment delivery ratio” module of the InVEST model to evaluate SC service. This module uses the revised universal soil loss equation (RUSLE) to express SC, according to the difference between potential soil loss and actual soil loss (Renard et al., 2017; Rao and Xiao, 2018). The calculation formula is as follows:

where $RKL{S}_i$ grid I is based on geomorphology and climate calculation of potential soil erosion; $R_i$, $K_i$, $L{S}_i$, $C_i$, and $P_i$ are the rainfall erosivity factor, soil erodibility factor, slope length factor, vegetation cover, and management factor, SC measure factor of the unit grid $i$; $USL{E}_i$ was the actual soil erosion amount of grid $i$, and $SEDRE{T}_i$ was the SC amount of grid $i;$ $C_i\left(P_i\right)$ is 0.05,0.03,0.04,0,0 (1,1,1,0,0) for cultivated land, woodland, grassland, water area, and construction land, respectively. Other factors are input based on the model recommendations. The details of each coefficient are shown in Supplementary Table SA3 of the supplementary data.

2.4 Spatial Autocorrelation of Ecosystem Services

Spatial autocorrelation can reflect the degree of interdependence and aggregation between the attributes of a specific region and those of other regions (Hu et al., 2022). The calculation of the Moran’I equation in ArcGIS 10.5 can be used to examine spatial patterns of ES. The calculation formula is as follows:

where n is the total number of each district (county) in Sichuan Basin; $S=1/n{\displaystyle {\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}$; $x_i$, $x_j$ is the value of ES of each district (County); $W_{ij}$ is the weight matrix of unit space; and $x$ is the average of ES. When $Moran\prime I$ > 0, the space is positively correlated, when $Moran\prime I$ < 0, the space is negatively correlated, and when $Moran\prime I$ = 0, the space is not correlated.

2.5 Spatial Differentiation of Driving Factors

2.5.1 Factor Identification and Data Dimensionality Reduction

First, we selected the factors associated with ES, including the four major factors: anthropogenic factors, climatic factors, vegetation factors, and geomorphological factors (Costanza et al., 2014; Zhou et al., 2018; Chen et al., 2019; Schirpke et al., 2019; Zhang Z. et al., 2021; Xu et al., 2022). To explore the overall correlation of ecosystem service impact factors, we conducted Spearman correlation research between CS, WY, and SC services and potential drivers based on SPSS 25. The final anthropogenic factors include night light, GDP density, and population density; climatic factors include annual average temperature, annual average precipitation, and actual evapotranspiration; vegetation factors include NDVI; and geomorphological factors include elevation, slope, and aspect (Table 2).

TABLE 2. Selected potential factors.

Second, to avoid the influence of multicollinearity of all variables, we mainly based on PCA to synthesize multiple related variables with the principle of as little information loss as possible and turn them into a few irrelevant variables to reduce the dimension of data, simplify the data structure and ensure that the VIF of explanatory variables is not greater than 7 (Plaza et al., 2005; Sannigrahi et al., 2020), thus to eliminate the possible influence of multicollinearity between variables (Supplementary Table SB1).

2.5.2 Local Regression Analysis

Geographically weighted regression (GWR) is a spatial regression method that can explain different local relationships among variables (Sun et al., 2020). Its greatest advantage is the ability to measure spatial variability accurately (Brunsdon et al., 1998). Although GWR captures the spatial heterogeneity of influencing factors to some extent, it is implemented under the assumption that all covariates change on the same spatial scale, and it is not enough to explain the spatial heterogeneity level of different urban landscape factors. MGWR relaxes the assumption of “same spatial scale” and allows optimization of covariate-specific bandwidths. The GWR and MGWR formulas are given in the equations:

where $x_{ij}$ is the J prediction variable and $\left(u_i,v_i\right)$ is the centroid coordinate of each region (county) i. ${\beta }_{bwj}$ is the optimal bandwidth for calibrating the j condition relationship. By contrast, MGWR allows different processes to operate on different spatial scales by deriving the individual bandwidth of the conditional relationship between the response variables and the different predictive variables. MGWR uses the inverse fitting algorithm to calibrate and initializes the inverse fitting process by estimating the GWR parameters. Based on these initial values, the calibration process works iteratively, with all local parameter estimates and optimal bandwidth evaluated in each iteration. The iteration is terminated when the difference in the parameter estimates of the successive iterations converges to the specified threshold.

3 Results

3.1 Temporal and Spatial Patterns of Ecosystem Services

From 2000 to 2015, the average values of CS and SC services increased from 5846835 t/ha and 6868780 t/ha to 5884678 t/ha and 7777234 t/ha, respectively, and the average WY decreased from 184656 m³/ha to 174054 m³/ha. The results showed that urbanization had a negative effect on the occupation of arable land, forest land, and carbon storage, and the implementation of ecological and water conservation projects had a positive effect on soil and water conservation. As far as the spatial distribution of ES was concerned (Figure 2), the CS service level was higher in all regions except Chengdu and Chongqing. The spatial distribution characteristics of SC services were “low in the middle and high in the periphery”. The high-value area of WY service was located in the northeast of Sichuan Province. On the whole, the spatial distribution of various services was different, which was closely related to the regional geographical background. In terms of the spatial change of ES, most regions showed a decrease in CS services, especially in cities and towns and their surroundings. Except for some regions of Wushan County, Chongqing, Wuxi County, and Chengkou County, SC services showed a decreasing trend, the other regions showed an increasing trend, and the decrease of WY services was mainly located in the northeast of the study area.

FIGURE 2. Spatial distribution of ES from 2000 to 2015. (A) 2000 annual water yield (WY); (B) 2015 annual WY; (C) 2000 annual carbon storage (CS); (D) 2015 annual CS; (E) 2000 annual soil conservation (SC); and (F) 2015 annual SC.

3.2 Spatial Autocorrelation of Ecosystem Services

Spatial autocorrelation analysis of CS, WY, and SC ES in the Sichuan Basin was conducted from 2000 to 2015 (Figure 3), in which the 2000 CS Moran’s I index was 0.359 and the standardized test Z (I) value was 10.073; in 2015, Moran ‘ s I index was 0.360 and Z (I) was 10.124, respectively, at the significance level of 0.05(p = 0); the results showed that the spatial distribution of CS showed a significant positive autocorrelation, and the state of carbon accumulation was apparent.

FIGURE 3. LISA agglomeration distribution of ES from 2000 to 2015 (A) 2000 annual WY; (B) 2015 annual WY; (C) 2000 annual CS; (D) 2015 annual CS; (E) 2000 annual SC; and (F) 2015 annual SC.

Moran’s I Index of WY in 2000 was 0.459, Z (I) value of the standardized test was 12.762, Moran’s I Index of WY in 2015 was 0.470, Moran’s I Index of WY in 2000 was 0.459, Moran’s I value was 12.762, Moran’s I value was 0.05, Moran’s I Index of WY in 2015 was 0.470, the Z (I) value of the standardized test was 13.017, and it passed the test at the significance level of 0.05(p = 0). The results showed that the spatial distribution of WY supply presented a significant positive autocorrelation and the aggregation was obvious.

In 2000, Moran’s I index was 0.347, Z (I) of the standardized test was 9.819, and it passed the test at the significance level of 0.05(p = 0) In 2015, Moran’s I index was 0.409, and Z (I) of the standardized test was 11.399, which passed the test at the significance level of 0.05 (p = 0) The results showed that the spatial distribution of SC supply showed a significant positive autocorrelation and aggregation.

At the same time, the local spatial autocorrelation analysis of ES, such as CS, WY, and SC from 2000 to 2015, was carried out, and the distribution chart of local indicators for spatial autocorrelation (LISA) was drawn. As shown in Figure 3, from the results of local spatial autocorrelation analysis, the low-low aggregation types of ES were mainly distributed in the central part of the Sichuan Basin, including Ziyang, Suining, Neijiang, Deyang, eastern Chengdu, and western Chongqing, while high-concentration types were mainly distributed in the eastern and western border areas of Sichuan Basin, the northern part of Chongqing, the northern part of Dazhou, the southern part of Ya’an, and the southern part of Leshan, i.e., the northern part of Chengdu, the southern part of Deyang, the northern part of Ziyang and the eastern part of Meishan. The other places had not passed the local spatial autocorrelation significance test.

3.3 Influencing Factors

3.3.1 Correlation Analysis

To explore the overall correlation of the impact factors of ES in the Sichuan Basin from 2000 to 2015, the results of the Spearman correlation between services and factors (Table 3) found that, the spatial distribution of ES capacity in the Sichuan Basin was influenced by both natural and socio-economic factors, thus showing different spatial distribution characteristics. In general, terrain factors such as TEM and EVA contributed the most, while socio-economic factors such as GDP and POP contributed less, and the contribution rate of other factors was low. Among them, POP, GDP, and LIG had a significant negative correlation with these three services. NDVI, DEM, and SLO had a significant positive correlation with these three services.

TABLE 3. Statistical description of Spearman correlation between ES and various impact factors.

3.3.2 Spatial Heterogeneity

SPSS 25 PCA was used to visualize the factors after 10 kinds of standardization, and the spatial pattern of coefficients was obtained (Figure 4), the spatial heterogeneity of different factors on WY, CS, and SC was further explored.

FIGURE 4. Spatial effect of four kinds of influencing factors (A). Geomorphological factors; (B) climatic factors; (C) anthropogenic factors; and (D) vegetation factors).

The spatial distribution of geomorphological factors in the Sichuan Basin had significant ring-shaped spatial features (Figure 4A). From the spatial distribution trend of influence intensity, with Ziyang, the central area of Sichuan Basin as the core, the negative driving force of topographic factors on the Sichuan Basin ecosystem was gradually increasing outwards and presented obvious ring-type distribution characteristics. The central part of the Sichuan Basin was affected by the surrounding mountainous terrain. The whole environment was closed and blocked, and the foundation of heat dissipation was poor.

The climate factors in Sichuan Basin had a negative correlation to the regional ES (Figure 4B) and had a significant gradient in the spatial distribution. In terms of the spatial distribution of the overall driving mechanism, except for the peripheral region of the Sichuan Basin, a few mountainous areas, and the central area of the Zhongliang mountain range in Jinyun, there was a strong positive correlation between driving characteristics, the climate of Chengdu-Chongqing urban agglomeration as a whole had an obvious positive effect on regional ES.

The anthropogenic factors in Sichuan Basin showed a significant negative correlation between the driving characteristics of regional ES and the spatial distribution of small patches (Figure 4C). The main driving force behind this phenomenon was that, on the one hand, the economic and industrial structure in Sichuan Basin was still in the stage of transition from primary production to secondary production and the basic demand for space and resources by industries such as industry and manufacturing led to the horizontal expansion of cities. On the other hand, the corridor of the Double-City Economic Circle in Chengdu and Chongqing was the main area of negative correlation, the spatial population density had a negative driving effect on the urban development, while a few marginal urban areas in the north showed a weak positive driving effect. From the perspective of the spatial distribution trend of driving intensity, the influence degree of the two core cities, Chengdu and Chongqing, gradually decreased outward.

The vegetation factors in Sichuan Basin showed a significant positive correlation to the regional ES and the patch spatial distribution trend (Figure 4D). As a whole, Sichuan Basin was dominated by farmland and grassland ecosystems and had prominent supply-type services. The increase in species diversity and community complexity would directly increase the stock of ES. At present, the stock of forest land ecological assets in each district of Chengdu-Chongqing urban agglomeration was insufficient, and the regulation, support, and cultural services in the core area of the central city were weak. Therefore, the enhancement and optimization of vegetation and other ecosystems had a significant effect on the enhancement of the Sichuan Basin’s overall service functions.

3.3.3 Multi-Scale Geographically Weighted Regression

The global model, GWR model, and MGWR model under the ordinary least squares method (OLS) were shown in Table 4. The result of OLS fitting showed that the VIF of all the factors was less than 7.5, which showed that there was no variable redundancy and there were no multiple linear relationships among the factors. At the same time, Jarque-Bera’s result showed that the residuals did not obey the normal distribution, the model fitting was one-sided, and the GWR model should be combined to improve the fitting accuracy.

TABLE 4. Comparison of the fitting results of OLS, GWR, and MGWR.

However, the local spatial regression model based on MGWR could better reflect the impact of the aforementioned influencing factors on WY services, CS services, and SC services in specific spaces, according to the results shown in Figure 5, which showed the local R² spatial distribution of the fitting results of the MGWR model. Figure 5A showed the impact of various factors on WY services, showing a vertical strip distribution in space, with the east being more intense, followed by the west, while the middle was less affected; Figure 5B showed the impact of various factors on CS services, showing a circular impact pattern opposite to WY services in space, that was, west > northeast > middle; Figure 5C showed the impact of various factors on SC services, showing local patches and overall axisymmetric distribution in space, and the northwest was the most affected.

FIGURE 5. Spatial heterogeneity distribution of ES under different factors (A). WY; (B) CS; and (C) SC.

From the descriptive statistical results of the coefficients of the MGWR Model (Tables 5, 6, and 7), it could be seen that the influence of different factors and their trend of action were quite different in-service spaces. DEM had the largest effect on CS, followed by POP on WY, while LIG had the least effect on SC and the temperature had the least effect on WY. From the ratio of positive and negative effects, the same factor had different effects on different services, which indicated that there was significant spatial heterogeneity of the factors, and the MGWR model could show the characteristics of local coefficients.

TABLE 5. Statistical description of the MGWR coefficient of WY services.

TABLE 6. Statistical description of the MGWR coefficient of CS services.

TABLE 7. Statistical description of the MGWR coefficient of SC services.

4 Discussion

4.1 Spatio-Temporal Patterns of Ecosystem Services in the Sichuan Basin

In the Sichuan Basin, one of the major suppliers of ES in China (Pan and Wang, 2021), carbon sequestration services and SC services showed an upward trend, while water production services showed a downward trend. These results were mainly attributed to the great changes in land use types in the study area. The increase in vegetation coverage led to a significant increase in CS and SC, which was consistent with previous studies (Jiang et al., 2018; Luo et al., 2019). At the same time, large-scale planting led to a significant increase in water consumption and evapotranspiration, which is the main reason for the decline of WY(Feng et al., 2016). However, the p-value of the Moran index of the three ES is 0, indicating that there is a significant positive autocorrelation at the significance level of 0.05. L-L cluster areas are mainly concentrated in the central part of the Sichuan Basin, and H-H cluster areas are distributed in the northeast and the periphery of the cities. In general, the Chengdu Plain, basin hills, and mountain areas around the basin in the east are the concentrated areas for WY supply services, while the southwest Sichuan mountains, west Sichuan mountain valleys, and northwest Sichuan Plateau in the west are the concentrated areas for CS and SC supply services.

4.2 Quantitative Attribution of the Driving Factors of Ecosystem Services in the Sichuan Basin

Combined with the results of the impact of natural and socio-economic driving factors on the correlation of ecosystem services in the Sichuan Basin, it can be found that LIG, GDP, POP, and TEM have significant negative effects on ES, indicating that the rise of these indicators will make the state of ES tense to some extent. NDVI, DEM, and SLO have significant positive effects on ES. These results are consistent with other studies but slightly different in the magnitude of the correlation coefficient (Lin et al., 2020), this is mainly due to the differences in research methods and scales, which provide a reference for enhancing ES in the Sichuan Basin research area (Sannigrahi et al., 2020; Huang et al., 2022). At the same time, according to the mean absolute value of regression coefficient of different factors, artificial factors such as LIG, GDP, and POP have a greater impact on WY services, the water production services are very sensitive to POP. CS services are sensitive to DEM and LIG. LIG, POP, and NDVI negatively affected CS services. SC services are sensitive to SLO and TEM, among which GDP, TEM, and SLO have negative effects on SC services, the rest were positive. The aforementioned results provide a basis for exploring the impact of different driving factors on ES in Sichuan Basin.

4.3 Comprehensive Response and Spatial Heterogeneity of Driving Factors to Ecosystem Services in the Sichuan Basin

Four comprehensive driving factors show different characteristics of the spatial pattern (Figure 4): geomorphological factors show the spatial characteristics of “centripetal trough”, climatic factors show the spatial characteristics of “east high west low”, anthropogenic factors show the spatial characteristics of “east-west dual-core”, and vegetation factors show the spatial characteristics of “planarization”. The contribution rate of geomorphological factors, TEM, and EVA to the three natural factors is higher, the contribution rate of GDP and POP to the three social-economic factors is higher, and the contribution rate of other factors is lower. POP, GDP, LIG, and other factors have a significant negative correlation with these three services, these variables are the main driving factors of ecosystem service changes, which are also consistent with previous relevant researchers (Wang et al., 2019; Yongxiu et al., 2020). It shows that the key to achieving coordinated development of regional ecosystems is to reduce the negative interference of human activities in the ecosystem.

Because of the regional and spatial differences in geological and meteorological conditions, the service and driving factors of CS, WY, and SC in the Sichuan Basin are characterized by spatial heterogeneity; the results of MGWR show that the AIC_C value of the MGWR model is smaller than that of OLS and GWR model, the adjusted R² is larger, and the bandwidth of variables is rich, which indicates that the fitting effect of MGWR is more reliable (Oshan et al., 2019; Li and Fotheringham, 2020), and can provide reference and guidance for regional units to adopt relevant policy regulation, according to local conditions.

Specifically, there is strong spatial heterogeneity in the driving factors of WY services in the Sichuan Basin (Supplementary Figure SA1), especially in Chongqing in the east and Chengdu in the west. It further explains that the eastern and western parts of the Sichuan Basin also have higher LIG and GDP, resulting in excessive demand for water in highly urbanized areas. At the same time, the altitude of surrounding areas leads to changes in precipitation, average temperature, and vegetation coverage, thus affecting WY services in external areas (Cao et al., 2020). To sum up, the Sichuan Basin River basin, which is composed of eastern Sichuan Ridge Valley, central Sichuan hilly area, and Chengdu Plain Lake area, has an uneven water production service space due to its special topography and spatial distribution of precipitation.

There is also a strong spatial heterogeneity in Sichuan Basin’s CS services (Supplementary Figure SA2). The service capacity of CS in peripheral areas such as northern counties and western counties and cities is relatively large, the main reason is that the People’s Government of Sichuan Province has taken the lead in setting up a natural forest resources protection project leading group to implement the system of chief responsibility and management by objectives responsibility for the natural forest resources protection projects within the areas under their jurisdiction, to prohibit the cutting of natural forests, and to close down the timber trading market in the project areas, and the logging ban was quickly extended to the entire province.

The SC services in Sichuan Basin still have spatial differences (Supplementary Figure SA3). The main reason is that the hilly area in the middle of Sichuan and the Chengdu Plain of Kawanishi is the most severe area of soil erosion in the upper reaches of the Yangtze River, accounting for 21.9% of the total area of soil erosion, the area is located in the middle and lower reaches of Min River, Tuo River, Jialing River, and the Yangtze River. The ground is cut and broken, the hills undulate and the west is plain. In most places in the region, the rainfall is about 1,000 mm, and most are heavy rain. The outcropping rock is mainly composed of purple-red sandstone and mudstone, which is easy to be weathered and eroded, resulting in severe hydraulic erosion in this area. However, the SC services in the mountainous areas around the basin are better, mainly due to the low, middle, and high mountains in the region, especially the widespread distribution of carbonate rocks and the intercalation of detrital rocks among them, as well as various human actions in recent years, such as the ecological protection of the red line, less population density distribution, and other comprehensive reasons, resulting in this area SC services steadily improved.

4.4 Limitations and Future Sustainable Development Measures

In this study, we used the InVEST model to estimate the value of CS, WY, and SC. However, the Sichuan Basin ecosystem is an integrated system with a wide variety of service functions (Kang et al., 2018; You et al., 2021). In addition to the three types of services considered in this study, food production, habitat quality, pollution control, entertainment, and culture are also important to service types, which are also important directions for follow-up research. On the other hand, the supply and demand of ecosystems together constitutes the flow of ES from natural systems to social systems (Burkhard et al., 2012; Burkhard et al., 2015). However, the current research focuses on the “supply” side analysis of ES, given the lack of research on the “demand” side of ES, the synergy between supply and demand of ES should be quantified in the future.

5 Conclusion

The concept of ES has built a “bridge” between natural ecosystems and human social systems. This study evaluates the spatial and temporal patterns of ES in the Sichuan Basin from 2000 to 2015 and discusses the spatial heterogeneity of response of various driving factors to ES with the help of PCA-MGWR model. The results showed that the CS and SC services experienced an upward trend during the study period, while the WY services showed a downward trend. Especially from the perspective of spatial pattern, the carbon sequestration services in other regions except for the main urban areas of Chengdu and Chongqing were higher; the spatial distribution characteristics of SC services are “low in the middle and high in the periphery”; the high-value area of water production service is located in Northeast Sichuan. The results of the spatial regression model show that compared with the OLS model and GWR Model, the MGWR model has the best explanatory power, and there is a correlation between driving factors and ES change. In view of the differences in the driving factors on the geographical scale, this study helps to consider the spatial differences in the impact of different indicators when formulating the countermeasures for ES in the Sichuan Basin and provides more information for the refined management and zoning control of ES in the Sichuan Basin.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Author Contributions

CH: writing—original draft and data curation. ZW: data curation and software. JL: editing and supervision. HL: supervision and writing—review and editing. DS: supervision. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Natural Science Foundation of China under the Grant 41971162 and the Second Scientific Investigation of the Tibetan Plateau (Grant no. 2019QZKK0406).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2022.927818/full#supplementary-material

References