Driving Factor Analysis of Ecosystem Service Balance for Watershed Management in the Lancang River Valley, Southwest China

: Revealing the spatio-temporal change of the supply, demand and balance of ecosystem services (ESs) associated with human activities and land-use changes is of great signiﬁcance for watershed ecosystem management. Taking the Lancang river valley as a case, we explicitly studied the ES spatial characteristics, using the land use/land cover (LULC) matrix model, Optimized Hot Spot Analysis and landscape pattern analysis. Furthermore, we screened out the dominant explanatory variables that had signiﬁcant inﬂuence on the ES supply, demand and balance by means of the Geographical Weighted Regression (GWR) method at pixel scale. The results showed that the ES demand intensity varied little throughout the watershed, while the downstream ES supply capacity and balance values were greater than upstream ones. Meanwhile, the hotspots of ES supply and demand were mainly distributed in the south part with coldspots in the north part. Human activity factors integrating landscape pattern variables were veriﬁed to have a negative impact on the ES balance in general. Among them, the Largest Patch Index (LPI) had a negative inﬂuence on the majority of pixels, while the Gross Domestic Product (GDP), cultivated land ratio and Area Weighted Average Patch Fractal Dimension (AWAPFD) had positive effects on a few pixels. This study will provide scientiﬁc support for regional ecosystem service trade-off and regulation at multiple scales.


Introduction
With the research development of ecosystem services (ESs), the relationship between ESs and human well-being have gradually become the focus of ES research [1][2][3]. James et al. (2007) defined ESs as components of nature, direct enjoyment, consumption or human well-being realization from the perspective of economic accounting [4]. However, the human well-being realization relying on the consumption of ESs will result in changes in the state of the ecosystem and ESs, thus affecting the continuous ES supply [5]. Therefore, regional ES consumption should be constrained within their supply range [6]. ES supply is the natural resource and service provided by natural components that can be directly utilized by humans [7]. Thus, some studies have emphasized that human beings are the direct beneficiaries [7,8]. In this context, quantitative assessment and scientific management of ES supply and demand have become the current research hotspots [9].

Study Area
The Lancang river valley is located in Southwest China (93 • 48 -101 • 48 E, 21 • 07 -33 • 48 N) with a total area of 165,383 km 2 flowing through Qinghai, Tibet, Sichuan and Yunnan Province [34]. The valley is in the shape of a narrow strip in the north-south direction, wide upstream and downstream, and narrow midstream. The landform types are complex and diverse, and the upstream originates from the Qinghai-Tibet Plateau (QTP) with many alpine valleys. The midstream is hilly with a rugged terrain and multiple terraces, and the downstream is of alluvial plains with a flat terrain and frequent floods. The terrain is undulant with great spatial heterogeneity. The elevation in the upstream exceeds 4500 m and the relative vertical distance in the midstream is about 2000 m with a deep valley, while the terrain in the downstream is gentle and the minimum elevation is only 486 m.
The vertical difference in climate is remarkable. The upstream belongs to the alpine climate zone of the QTP with low temperature and low rainfall [35]. The mean annual temperature varies between −3 • C and 3 • C, and the mean annual precipitation ranges from 400 to 800 mm [35]. The midstream belongs to the subtropical climate of northwestern Yunnan, and the vertical change in temperature is obvious, which increases from north to south [36]. The mean annual temperature ranges between 12 • C and 15 • C, and the annual precipitation varies between 1000 and 2500 mm [36]. The downstream belongs to the tropical and subtropical climate of the southwestern Yunnan with temperatures increasing from north to south [37]. The mean annual temperature ranges between 15 • C and 22 • C, and the annual precipitation varies between 1000 and 3000 mm [37]. Overall, the whole basin belongs to the southwest monsoon climate with distinct dry and wet seasons [38]. The valley is very rich in water resources, and the mean annual runoff is approximately 4750 × 10 8 m 3 , which is about 9-fold greater than the annual runoff of the Yellow River [39]. The complex terrain and climate features form rich vegetation types with obvious latitude zonality and vertical zonality ( Figure 1). Vegetation types mainly include alpine vegetation, meadows, shrubs, coniferous forests, broad-leaf and mixed forests, broad-leaf forests, cultivated vegetation and swamps, while meadows and shrubs are the main vegetation types of upstream [40]. The main land use types in the valley include agricultural land, forest land, grassland, urban land, water body, bare land and wetland. The proportion of cultivated land is relatively low, accounting for only 8.43% of the total area of the valley distributed in the downstream ( Figure 1). Apart from that, there are large areas of bare rock and bare soil in the north part of the valley, accounting for 7.04% of the total area.

Data Sources and Processing
(1) Land use and land cover.
In this study, land cover data with a 100 m × 100 m resolution in 2015 were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC, Beijing) (http://www.resdc.cn, 1 September 2020). Based on the land cover data in 2015, we calculated the proportion of cultivated land and landscape pattern indexes, including the Largest Patch Index (LPI), Area Weighted Fractal Dimension of Patches (AWMPED) and Contagion Index (CONTAG) by Fragstats 4.2. Among them, LPI represents shape, AWMPED represents area and CONTAG represents aggregation aspects.
In this study, GDP spatial distribution grid data with a 1 km × 1 km resolution in 2015 were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC, Beijing) (http://www.resdc.cn, 5 September 2020). Township location point data were obtained from the GeoHey (https://geohey.com/, 10 September 2020). The road distribution data were obtained from the National Agriculture Science Data Sharing Center (http://www.agridata.cn, 8 September 2020).

Data Sources and Processing
(1) Land use and land cover.
In this study, land cover data with a 100 m × 100 m resolution in 2015 were prov by the Data Center for Resources and Environmental Sciences, Chinese Academy o ences (RESDC, Beijing) (http://www.resdc.cn, 1 September 2020). Based on the land c data in 2015, we calculated the proportion of cultivated land and landscape patter dexes, including the Largest Patch Index (LPI), Area Weighted Fractal Dimensio Patches (AWMPED) and Contagion Index (CONTAG) by Fragstats 4.2. Among them represents shape, AWMPED represents area and CONTAG represents aggregatio pects.
In this study, GDP spatial distribution grid data with a 1 km × 1 km resolution in were provided by the Data Center for Resources and Environmental Sciences, Ch Academy of Sciences (RESDC, Beijing) (http://www.resdc.cn, 5 September 2020).
Township location point data were obtained from the GeoHey (https://geohey.c 10 September 2020). The road distribution data were obtained from the National A culture Science Data Sharing Center (http://www.agridata.cn, 8 September 2020).

Quantification and Spatial Patterns of ES Supply, Demand and Balance
According to the LULC matrix model constructed by Burkhard

Quantification and Spatial Patterns of ES Supply, Demand and Balance
According to the LULC matrix model constructed by Burkhard et al., (2015) [8], we obtained the ES supply matrix and demand matrix by merging 44 land use types and 22 ES types of the original matrix into 25 land use types and 22 ES types ( Figure 2). We assigned different values to the supply capacity and demand intensity of various ES types, including provisioning service (11 types), regulating service (9 types), and cultural service (2 types) for each land use type. The ES supply matrix was used to define the supply capacity of different land use types for various ESs (Figure 2a), while the ES demand matrix was used to reflect the actual demand intensity of humans for ESs for each land use type (Figure 2b). In these two matrixes, the evaluation levels of ES supply capacity and demand intensity vary from 0 to 5. It was reported that this method is a propagable and convenient method for ES quantification [3]. By subtracting the ES demand intensity by the ES supply capacity, we calculated the ES balance to identify the surplus and deficit of ES ( Figure 2c). The evaluation scale ranges between −5 and 5. Negative values indicate that the demand for ESs exceeds the supply in a deficit state, and positive values indicate that ESs are in a surplus state. Based on the derived assessment matrixes of land use, the spatial quantitative analysis of ES supply, demand and balance of various ESs was carried out in the Lancang river valley. Then, the ES supply, demand, balance and total values were calculated by multiplying the area of 22 land use types by the quantitative values of ES supply, demand and balance to reflect the spatial variation of the ESs in the Lancang river valley. The formulas used are as follows: where V ESi is the quantitative value of the ES category for land cover types i (i = 1, 2, 3 . . . 22); S i is the area of land cover i; ES j is the ES value for ES categories j. In order to further explore the local performance of ES changes in space, Optimized Hot Spot analysis/Getis-Ord G* index in ArcGIS 10.5 was used to analyze the local correlation between coldspots and hotspots. The spatial pattern of the total ES supply and demand for each pixel (3 km × 3 km) was recognized to identify spatial clusters at a 99% significance level in the Lancang river valley in 2015. High-value clustered pixels were counted as hotspots; low-value clustered pixels were counted as coldspots; the high-value To verify the reliability of the method, the coefficient of sensitivity (CS) was further introduced to evaluate the response of total ES value to the change of the quantitative value. The formula used is as follows: where i and j are the initial and adjusted values, respectively; adjustment quantity = ±50%; k is the land use type; VC is the quantitative value; and ESV is the total ES value. If CS ≥ 1, then ESV is flexible to VC, and the assessment method has poor accuracy and low credibility. If CS < 1, then ESV is inelastic to VC, and the method is credible. In order to further explore the local performance of ES changes in space, Optimized Hot Spot analysis/Getis-Ord G* index in ArcGIS 10.5 was used to analyze the local correlation between coldspots and hotspots. The spatial pattern of the total ES supply and demand for each pixel (3 km × 3 km) was recognized to identify spatial clusters at a 99% significance level in the Lancang river valley in 2015. High-value clustered pixels were counted as hotspots; low-value clustered pixels were counted as coldspots; the high-value and low-value interval distributions were insignificant pixels. The formulas used are as follows: where W ij is the spatial weight; X i and X j are the observation values of the ith ESV evaluation unit and the jth ESV evaluation unit, respectively; and n is the sample size.

Driving Factor Analysis of ES Balance Based on Geographically Weighted Regression Model
The Geographical Weighted Regression (GWR) model is a kind of regression model that can detect the changes in the variable structure and variable relationship caused by geographical location changes, which is called spatial non-stationarity [41]. Compared with the general linear regression model, the GWR model effectively improves the goodnessof-fit, which has been widely used in spatial regression analysis [42][43][44]. The GWR model was calculated according to Equation (7) in ArcGIS 10.5.
where (u i , v i ) is the spatial position of the ith sample; x ik is the independent variable; is the weight in the regression equation calculated according to the attenuation function, which is used to modify the parameter x ik ; p is the number of samples and ε i is the error correction term. The analysis process of the driving factor based on the GWR model is described as follows.
In this study, 7 factors were initially selected to reflect human activity intensity as explanatory variables, including 4 direct reflection factors (GDP, cultivated land ratio, distance to road and distance to township point) and 3 indirect reflection factors describing the degree of landscape fragmentation (LPI, AWMPED and CONTAG). The balance values of four ES types were selected as the dependent variables, including regulating service balance, provisioning service balance, cultural service balance and total ES balance. In the GWR model, all variables were standardized, using the deviation standardization method.
We divided the Lancang river valley into 18,325 grids (3 km × 3 km) and extracted the values of 4 dependent variables and 7 explanatory variables of each grid. Based on the above data, GWR analysis was performed in the 3 following steps: Significant explanatory variables screening: The Ordinary Least Squares (OLS) method was used to carry out the regression analyses on 4 dependent variables and 7 explanatory variables to eliminate the insignificant variables and redundant variables, respectively. The screening criteria were that p-value was less than 0.05 and the Variance Inflation Factor (VIF) was less than 7.5. Then, OLS was performed again on the screened explanatory variables, and the obtained R 2 was the goodness-of-fit of the global regression model. According to the OLS analysis, 4 corresponding explanatory variables were selected for provisioning service balance, including LPI, GDP, cultivated land ratio and distance to township points. For cultural services balance, 5 of the 7 corresponding explanatory variables (GDP, distance to road, distance to township point, LPI and AWMPED) were selected, except the CONTAG and cultivated land ratio. For regulating service balance and total ES balance, all of the 7 initial explanatory variables were screened out as corresponding explanatory variables.
GWR model performing: Using the Spatial Statistics/GWR tool in ArcGIS10.5, GWR model was respectively performed on 4 dependent variables and their corresponding explanatory variables. The weight function calibration method was adaptive, and the bandwidth determination method was AICc.
Moran's I verification: The Spatial Autocorrelation (Moran's I) tool was used to calculate the global Moran's I of standardized residual in GWR model to test whether it was in a random distribution state.
(1) Spatial stationarity analysis of GWR regression coefficients, and spatial differentiation analysis of dominant explanatory variables.
The coefficient of variation (C v ) was calculated based on each regression coefficient in the GWR model to describe the spatial stationarity as follows: where C v was the coefficient of variation; σ was the standard deviation; x was the average value. The influence degree of various explanatory variables on ES balance value can be reflected by the regression coefficients of the pixel. In this study, we screened the maximum regression coefficients at pixel scale and performed the spatial differentiation analysis, which reflected the dominant explanatory variables of each pixel and their influence directions.
The research framework and corresponding methods were shown in Figure 3. lected for provisioning service balance, including LPI, GDP, cultivated land ratio and distance to township points. For cultural services balance, 5 of the 7 corresponding explanatory variables (GDP, distance to road, distance to township point, LPI and AWMPED) were selected, except the CONTAG and cultivated land ratio. For regulating service balance and total ES balance, all of the 7 initial explanatory variables were screened out as corresponding explanatory variables. GWR model performing: Using the Spatial Statistics/GWR tool in ArcGIS10.5, GWR model was respectively performed on 4 dependent variables and their corresponding explanatory variables. The weight function calibration method was adaptive, and the bandwidth determination method was AICc.
Moran's I verification: The Spatial Autocorrelation (Moran's I) tool was used to calculate the global Moran's I of standardized residual in GWR model to test whether it was in a random distribution state.
(1) Spatial stationarity analysis of GWR regression coefficients, and spatial differentiation analysis of dominant explanatory variables.
The coefficient of variation ( ) was calculated based on each regression coefficient in the GWR model to describe the spatial stationarity as follows: where was the coefficient of variation; σ was the standard deviation; ̅ was the average value.
The influence degree of various explanatory variables on ES balance value can be reflected by the regression coefficients of the pixel. In this study, we screened the maximum regression coefficients at pixel scale and performed the spatial differentiation analysis, which reflected the dominant explanatory variables of each pixel and their influence directions.
The research framework and corresponding methods were shown in Figure 3.

Spatial Distributions of ES Supply, Demand and Balance
As shown in Figure 4, the supply capacity of the provisioning service in the study area was low in the upstream and high in the downstream, similar to the spatial differentiation of elevation and climate [13,45]. The demand intensity of the provisioning service was generally low, and the high-value areas were scattered in the downstream, mainly in the Dali and Xishuangbanna cities. The deficit areas of provisioning service showed a continuous distribution in the upstream and a scattered distribution in the downstream. For regulating service, the supply capacity was overall high, and the low-value areas were mainly distributed in the Yushu, Puer and Lincang cities. The demand intensity of the regulating service was overall low with higher values in the downstream than in the upstream. Dali was the region where the high value of demand intensity was relatively concentrated. For cultural service, the supply capacity was high, especially in the southeastern part of the downstream, including Xishuangbanna and Puer. The demand of cultural service was low scattered in the whole study area. For total ES, the spatial characteristics of supply capacity was high in the downstream and low in the upstream. Although the demand intensity of total ES was also high in the downstream, the surplus of the downstream was much higher than that of the upstream, with only a small number of deficit areas scattered. In the upstream, deficit areas of total ES had a continuous distribution, but the degree of deficit was relatively low.

Spatial Patterns of Total ES Supply and Demand
The supply and demand of total ES in the study area were in a local, spatial aggregation pattern, as shown in Figure 5. For the total ES supply, the area of hotspots accounted for 35.87% of the Lancang river valley, concentrated in the eastern part of the downstream and scattered in the west. The area of coldspots accounted for 34.79% of the study area, concentrated in the upstream and scattered in the downstream. Compared with the total ES supply, the hotspots and coldspots of total ES demand represented a more dispersed distribution, being mainly distributed in Dali city and Nujiang Lisu Autonomous Prefecture with a small amount distributed in the western part of the downstream. The area of demand hotspots accounted for 18.84% of the total area. The number of coldspots was small, accounting for 13.71% of the study area, scattered in the Yushu and Changdu cities located in the upstream. When comparing the spatial pattern of the total ES supply with that of demand, the distribution of supply coldspots was more concentrated and northerly than that of demand coldspots. Both supply hotspots and demand hotspots were distributed in the downstream, representing mosaic distribution in the east-west direction, which explained the large spatial differences of total ES balance in the southern part of the study area.

OLS Analysis and Corresponding Explanatory Variable Selection
Comparing the OLS results of four dependent variables revealed that LPI was the major impact factor of all dependent variables ( Table 1). The correlation coefficients were all greater than 0.7 and all were negatively impacted.

Spatial Patterns of Total ES Supply and Demand
The supply and demand of total ES in the study area were in a local, spatial aggregation pattern, as shown in Figure 5. For the total ES supply, the area of hotspots accounted for 35.87% of the Lancang river valley, concentrated in the eastern part of the downstream and scattered in the west. The area of coldspots accounted for 34.79% of the study area, concentrated in the upstream and scattered in the downstream. Compared with the total ES supply, the hotspots and coldspots of total ES demand represented a more dispersed distribution, being mainly distributed in Dali city and Nujiang Lisu Autonomous Prefecture with a small amount distributed in the western part of the downstream. The area of demand hotspots accounted for 18.84% of the total area. The number of coldspots was small, accounting for 13.71% of the study area, scattered in the Yushu and Changdu cities located in the upstream. When comparing the spatial pattern of the total ES supply with that of demand, the distribution of supply coldspots was more concentrated and northerly

Validation of GWR Model and Spatial Stationarity of Explanatory Variables
Comparing the adjusted R 2 of GWR and OLS demonstrated that the goodness-of-fit values of GWR models were higher than those of the corresponding OLS models in various ESs, indicating that the regression effect of the GWR model was better than that of the OLS model (Table 1). Meanwhile, based on the global Moran's I, the z-scores of all four dependent variables were less than 1.65 and p-values were greater than 0.1, indicating that the standardized residuals showed random distribution and the GWR models passed the validation of the global Moran's I ( Table 1). than that of demand coldspots. Both supply hotspots and demand hotspots were distributed in the downstream, representing mosaic distribution in the east-west direction, which explained the large spatial differences of total ES balance in the southern part of the study area.

OLS Analysis and Corresponding Explanatory Variable Selection
Comparing the OLS results of four dependent variables revealed that LPI was the major impact factor of all dependent variables ( Table 1). The correlation coefficients were all greater than 0.7 and all were negatively impacted.

Validation of GWR Model and Spatial Stationarity of Explanatory Variables
Comparing the adjusted R 2 of GWR and OLS demonstrated that the goodness-of-fit values of GWR models were higher than those of the corresponding OLS models in various ESs, indicating that the regression effect of the GWR model was better than that of the For the explanatory variables in the GWR model, the smaller the C v of the regression coefficient was, the larger the spatial stationarity and the more stable impact on ES balance were. In contrast, the larger C v was, the greater the fluctuation in space was and the more unstable the impact was. For different dependent variables, the spatial stationarity of different explanatory variables was as follows: (i) for the provisioning services balance, the explanatory variable with the largest spatial variability was GDP and the most stable explanatory variable was LPI; (ii) for the regulating services balance, the spatial stationarity of different explanatory variables was not much different, but the ratio of cultivated land showing the highest stability, and the explanatory variable with the largest spatial variability was GDP; (iii) for the cultural services balance, the spatial variabilities of the distance to road and distance to township points were the largest; and (iv) for the total ES balance, the explanatory variables with the largest spatial variabilities were GDP and AWMPED (Table 1).

Spatial Differentiation Analysis of Dominant Explanatory Variables
The screening results of dominant explanatory variables at pixel scale indicated that four kinds of ES balance were affected by seven types dominant explanatory variables with negative or positive directions ( Figure 6). For provisioning service balance, 93.81% of the pixels were negatively affected by LPI, and 4.95% of the pixels were positively affected by GDP. The remaining 1.24% of the pixels were affected by the other three variables. For regulating service balance, 74.70% of the pixels were negatively affected by the cultivated land ratio and 21.52% of the pixels were negatively affected by LPI. In addition, 2.58% of the pixels were negatively affected by GDP, concentrated in Yushu city. Cultural service balance only had two dominant explanatory variables. Specifically, 92.07% of the pixels were negatively affected by LPI, and the remaining pixels were negatively affected by GDP. The dominant explanation variables of total ES balance included six types, and 92.88% of the pixels were negatively affected by LPI. Overall, the impact of human activities on all kinds of ES balance was basically negative. Only the GDP, cultivated land ratio and AWMPED showed a positive influence on a few pixels.

Spatial Variability of the ESs in the Lancang River Valley
The spatial patterns of the ESs indicated that the supply, demand and balance of various ESs presented low values in the upstream and high values in the downstream, which showed similar trends to the topographic and climatic distribution characteristics in study area [13,45]. The fundamental reason was that the spatial distribution of land use types is restricted by both natural conditions and human activities [46]. The restriction of natural conditions brought by the topographic and climatic features, and the contradictions between human and land jointly affected the land use status, thereby influencing the ES supply and demand [47].
Previous studies focusing on the values of ESs have proved that the ES values presented a regular distribution on the topographic and climatic gradients, which were affected by both elevation, slope, temperature and precipitation [45,48]. Taking the upper stream of the Min river as the study area, Zhu et al. (2017) analyzed the impact of the terrain niche on the ES values from the perspective of the selectivity of various land use types on topography [49]. Their results showed that, driven by the topography limitation and the human-land contradiction, the woodland and cultivated land with better topographic conditions were occupied by construction land, grassland and water, and the ES values varied in an inverted "V" shape with the rising topography [49]. Fang et al. (2020) quantified ES variations in elevation-associated vegetation types in a tropical region of China and found that ESs showed spatial differences along with different elevation levels in the Xishuangbanna area [47]. Furthermore, Jiang et al. (2016) assessed the ecosystem changes in the Three-River Headwater Region, and found that the increasing precipitation-intensified water erosion, and elevating temperature-induced glacier melting and permafrost degradation pose a threat to the ESs [50]. Lin et al. (2015) also concluded that climate, topography, land-use patterns and vegetation cover, as well as the interactions among them were the main factors that affected the spatial patterns of soil retention ESs in the Three Parallel Rivers Region [30].
The elevation is high in the upstream of the study area, and the midstream has a large elevation drop and steep slope, and the elevation in the downstream is constant. Thus,

Spatial Variability of the ESs in the Lancang River Valley
The spatial patterns of the ESs indicated that the supply, demand and balance of various ESs presented low values in the upstream and high values in the downstream, which showed similar trends to the topographic and climatic distribution characteristics in study area [13,45]. The fundamental reason was that the spatial distribution of land use types is restricted by both natural conditions and human activities [46]. The restriction of natural conditions brought by the topographic and climatic features, and the contradictions between human and land jointly affected the land use status, thereby influencing the ES supply and demand [47].
Previous studies focusing on the values of ESs have proved that the ES values presented a regular distribution on the topographic and climatic gradients, which were affected by both elevation, slope, temperature and precipitation [45,48]. Taking the upper stream of the Min river as the study area, Zhu et al., (2017) analyzed the impact of the terrain niche on the ES values from the perspective of the selectivity of various land use types on topography [49]. Their results showed that, driven by the topography limitation and the human-land contradiction, the woodland and cultivated land with better topographic conditions were occupied by construction land, grassland and water, and the ES values varied in an inverted "V" shape with the rising topography [49]. Fang et al., (2020) quantified ES variations in elevation-associated vegetation types in a tropical region of China and found that ESs showed spatial differences along with different elevation levels in the Xishuangbanna area [47]. Furthermore, Jiang et al., (2016) assessed the ecosystem changes in the Three-River Headwater Region, and found that the increasing precipitation-intensified water erosion, and elevating temperature-induced glacier melting and permafrost degradation pose a threat to the ESs [50]. Lin et al., (2015) also concluded that climate, topography, land-use patterns and vegetation cover, as well as the interactions among them were the main factors that affected the spatial patterns of soil retention ESs in the Three Parallel Rivers Region [30].
The elevation is high in the upstream of the study area, and the midstream has a large elevation drop and steep slope, and the elevation in the downstream is constant. Thus, due to the undulating terrains, the land use types vary greatly in different parts of the Lancang river valley. Integrated with the driving effects of human activities, the study area forms the following land use characteristics: the downstream is mainly covered by woodland, construction land and cultivated land, and the upstream is mostly covered by grassland, bare soil and bare rock. Our study revealed the spatial variability of ESs at the watershed scale and found that the land use characteristics affect the spatial distribution of ES supply, demand and balance, similar to the distribution of topographic features [47].

Links between Landscape Pattern Indices on ES Balance
Landscape patterns can quantify the landscape structural composition and spatial configuration to reflect landscape dynamics and the ecological process [51]. Landscape composition and configuration impact the ability of ecosystems to provide goods and services, and the improvement or degradation in ESs can be detected by landscape pattern changes [52,53]. Landscape pattern indexes may be a valid tool to explore the effects of land-use changes and human activities on ESs [54,55]. The landscape pattern changes will affect the variation in landscape indexes and, thus, affect the change in ES supply and demand [56]. Therefore, combining landscape pattern indexes and ESs will contribute to exploring how landscape patterns influence the ES balance.
The results of our study showed that the ES balance was significantly correlated with the landscape pattern indexes, including LPI, AWMPED and CONTAG (Table 1). LPI was the dominant explanatory variable of ES balance for the majority pixels, characterized by strong influence and high spatial stability ( Figure 6). Previous studies proved that the changes of landscape patterns had a regional impact on ES values [57]. Wang et al. found that there were significant positive correlations between ES values and landscape richness and aggregation, which showed a strong negative correlation between ES values and the complexity of patch shape [58]. The reliable links between landscape indexes and ESs in terms of land pattern and their implications should still be determined. In addition, the landscape patterns play a key role in the spatial variation of ES balance (Figure 6), consistent with the results of previous studies [52,59]. Yushanjiang et al. also quantified the spatial correlations between landscape patterns and ES values in Ebinur Lake Basin, Xinjiang, and their results showed that landscape metrics had direct spatial correlations with different ES types [52]. LPI had negative correlation with various ES balances, which can influence landscape fragmentation and connectivity [60,61]. It has been proved that social policies can effectively drive changes in landscape patterns, thus affecting the ES supply and demand [62]. Results of this research may provide scientific support for the formulation of regional ecosystem management policies from landscape perspectives.

Application and Advantages of GWR Model
Recently, several studies have proved that human activities and landscape patterns are important factors influencing ES supply, demand and balance [51,55]. Our results contributed to providing effective information and guidance for ES management. Previous studies on the driving factors of ESs focused on traditional analyses, such as correlation analysis, cluster analysis (CA) and principal component analysis (PCA) [63]. These methods could infer a preliminary assessment result about which human activity factors will have significant influence on ES balance. However, these studies lacked exploring the spatial distribution features of human activity factors' effects on ES balance on a pixel scale, so they could not screen out the specific sites and types of human activity on spatial scale.
In this study, we used the OLS and GWR models to evaluate the driving influence of different human activity factors and landscape pattern indexes on ES balance, which could supply more detailed information for policy decision-making. The results of the OLS model showed that LPI had the greatest negative influence on various service types of ES balance, which has been proved by the front section. In agreement with the OLS model analysis, the GWR model also revealed that LPI showed the strongest negative effect on ES balance for the majority of the pixels ( Figure 6). However, the difference between these two methods was that the GWR model could determine the weight of each sample point in the regression equation by establishing the distance attenuation function, whereas the OLS analysis could not. The constructed local regression equation conformed to the first law of geography and showed a higher goodness-of-fit than that of the OLS model (Table 1), which was consistent with the results of previous studies [42,44]. For example, Wu (2020) explored the spatio-temporal variations of the ecological footprint and its influencing factors in China's provinces using the GWR model and concluded that the GWR model was superior to OLS considering the regression analysis of the goodness-of-fit, variance comparisons and the spatial auto-correlation of residual [42]. In addition, compared with global regression, the GWR model can obtain the influence degree of different human activity factors on the ES balance at pixel scale. All these findings imply the advantages of the GWR model over the OLS model in exploring the spatial distribution and driving force of ES balance at regional and spatial scales. By screening the dominant explanatory variables, the results can provide more detailed and effective guidance for the regulation of regional ESs.

Limitations and Implications for Ecological Management at the Watershed Scale
In this study, we applied the LULC matrix model to quantify and map ES supply, demand and balance in the Lancang river valley. According to the previous studies, due to the simplicity, the LUCC matrix method has been widely applied to quantify ES supply, demand and balance in many areas [2,8,64]. However, there are some uncertainties and criticisms for this method. This method was an expert-based evaluation method with high subjective limitations and lacked conformity [65]. Only land-use data were used in the ES assessment, but other influencing factors, such as population and GDP, were ignored, which may bring estimation mistakes. These above uncertainties would bring uncertain results for ES balance estimation with great discrepancies [7]. Thus, we further calculated the CS to verify the reliability of the LULC assessment method. The results showed that the CS values for various land use types were less than 1, which indicated the rationality and reliability of our method and results [2,17]. Therefore, in future studies, all possible uncertainties should be further considered at watershed scale to provide more valid assessment and support for ES management at the watershed scale.
In addition, many human activity factors and landscape indexes, namely GDP, cultivated land ratio, distance to road, distance to township point, LPI, AWMPED and CONTAG, have been used to explore the driving forces of ES balance in the Lancang river valley. These factors are some representative factors that we selected to better reflect the human activity intensity in the study. In previous studies, Zhao et al., (2015) selected population density, number of villages and road length to evaluate the human influence intensity on the QTP [66]; Li et al., (2018) introduced the electricity infrastructure, population density, land use intensity, road and railways, and grazing intensity to reflect human activity [27] but these human activity factors are insufficient. In addition, GDP, cultivation activity and landscape pattern indexes may be also important factors influencing human activity [3,52]. Therefore, compared to previous studies, we had a comprehensive consideration for the influence of human activities on ESs, which provided more reliable results. However, some other human activity and landscape factors may also influence the ES balance, such as population density, traffic volume, night-time light and tourism belonging to human activity factors, and patch density (PD), landscape division index (DIVISION) and landscape shape index (LSI) belonging to landscape pattern indexes [27,52,57]. The inadequate consideration of these impact factors in the Lancang river valley may underestimate the assessment results. Many studies showed that tourism and traffic volume would bring a significant and negative influence on ESs [67,68], but these data are usually difficult to obtain at watershed scale; thus, we excluded the data. Many potential human activity factors may threaten the ecosystem, thus more attention should be paid to consider all the possible factors of human activity intensity at the watershed scale in future studies.

Conclusions
Regional land use status is an important factor driving the change in ES supply and demand. In this study, based on the LULC matrix model, four kinds of ES supply, demand and balance in the Lancang river valley were quantified, and the spatial distribution maps of various ESs were obtained. On this basis, we analyzed the spatial pattern of total ES supply and demand and obtained the quantity and spatial characteristics of the hotspots and coldspots of total ES supply and demand, which performed as hotspots in the south and coldspots in the north. Furthermore, using the OLS and GWR models, we conducted a driving factor analysis on various ES balances, and screened human activity factors with significant influence. According to the regression coefficients of each explanatory variable, the dominant explanatory variables of each pixel were obtained. The results showed that the provisioning service and cultural service of the ES balance were significantly affected by four and five human activity factors, respectively, and regulating service of the ES balance and total ES balance were significantly affected by seven human activity factors. Among the factors, LPI was proved to be the dominant explanatory variable of the vast majority pixels and its influence direction was negative. The research also indicated that the impacts of human activities on various service types of ES balances were basically negative, and only GDP, cultivated land ratio and AWMPED had a positive influence on a small quantity of pixels. The research may provide scientific support for regional ES management and regulation.