Quantiﬁcation of Natural and Anthropogenic Driving Forces of Vegetation Changes in the Three-River Headwater Region during 1982–2015 Based on Geographical Detector Model

: The three-river headwater region (TRHR) supplies the Yangtze, Yellow, and Lantsang rivers, and its ecological environment is fragile, hence it is important to study the surface vegetation cover status of the TRHR to facilitate its ecological conservation. The normalized difference vegetation index (NDVI) can reﬂect the cover status of surface vegetation. The aims of this study are to quantify the spatial heterogeneity of the NDVI, identify the main driving factors inﬂuencing the NDVI, and explore the interaction between these factors. To this end, we used the global inventory modeling and mapping studies (GIMMS)-NDVI data from the TRHR from 1982 to 2015 and included eight natural factors (namely slope, aspect, elevation, soil type, vegetation type, landform type, annual mean temperature, and annual precipitation) and three anthropogenic factors (gross domestic product (GDP), population density, and land use type), which we subjected to linear regression analysis, the Mann-Kendall statistical test, and moving t -test to analyze the spatial and temporal variability of the NDVI in the TRHR over 34 years, using a geographical detector model. Our results showed that the NDVI distribution of the TRHR was high in the southeast and low in the northwest. The change pattern exhibited an increasing trend in the west and north and a decreasing trend in the center and south; overall, the mean NDVI value from 1982 to 2015 has increased. Annual precipitation was the most important factor inﬂuencing the NDVI changes in the TRHR, and factors, such as annual mean temperature, vegetation type, and elevation, also explained the vegetation coverage status well. The inﬂuence of natural factors was generally stronger than that of anthropogenic factors. The NDVI factors had a synergistic effect, exhibiting mutual enhancement and nonlinear enhancement relationships. The results of this study provide insights into the ecological conservation of the TRHR and the ecological security and development of the middle and lower reaches.


Introduction
The vegetation cover is an important component of surface ecosystems that connects the atmosphere, hydrosphere, pedosphere, and areas inhabited by humans [1]; thus, the study of regional vegetation cover is essential to regional ecological conservation. The normalized difference vegetation index (NDVI) can accurately reflect the status of surface vegetation cover, which is the best indicator of vegetation coverage and the most effective indicator for monitoring regional vegetation change and the ecological environment [2,3]. Regional vegetation coverage changes and their drivers have been studied at different and alpine grasslands. The Qinghai-Tibet Plateau is a vast semi-natural area with relatively little artificial influence [23]. The TRHR is located in the central part of the Qinghai-Tibet Plateau, with high altitude and sparse population. The TRHR is an important ecological barrier in China with a fragile ecological environment; therefore, its ecological conservation is crucial for the sustainable development of a vast area in China.

Data and Processing
The data included in this study included natural factors, such as the NDVI, digital elevation model (DEM), climate data, landform type, soil type, and vegetation type, as well as anthropogenic factors, such as land use type, population density, and GDP in the TRHR. NDVI data were obtained from the Big Earth Data Platform for Three Poles, using GIMMS NDVI3g data with a spatial resolution of 8 km and a temporal resolution of 16 days [24,25], the NDVI images for each year from 1982 to 2015 were obtained by maximum value composite (MVC) [26]. Vegetation coverage was divided into five classes according to NDVI values: low coverage (≤0.2), medium-low coverage (0.2-0.4), medium coverage (0.4-0.6), medium-high coverage (0.6-0.8), and high coverage (> 0.8). The annual mean temperature and annual precipitation data were obtained from daily standard meteorological data of 26 meteorological stations in and around the TRHR from 1982 to 2015 using the inverse distance weighting method. DEM data were GDEMV2 30 m resolution digital elevation data from the Geospatial Data Cloud of the Chinese Academy of Sciences (http://www.gscloud.cn/), and the elevation, slope, and aspect data were obtained from the DEM data. Other data were obtained from the Resource and Environmental Sciences Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/). All the above data were extracted according to the vector boundary of the TRHR [27] and were resampled to make them consistent with the 8-km NDVI data image. Using ArcGIS to create a fishnet, 5,853 random sampling points were generated based on 8 × 8 km grids according to the scope of the TRHR, and the spatial attribute values were extracted.
We selected NDVI as the dependent variable and six categories of topography (slope, aspect, elevation), soil (soil type), vegetation (vegetation type), landform (landform type), climate (annual mean temperature, annual precipitation), and human activity (GDP, population density, land use type), a total of 11 representative and easily quantifiable factors, as independent variables. Precipitation and temperature are important factors affecting vegetation changes [21], elevation, slope, aspect, and landform type affect vegetation distribution by changing water and heat conditions [28]. Soil and vegetation types are important environmental elements for vegetation growth [29]; economic development affects ecological environment, land use type, GDP, and population density are indicators that can quantify changes in socioeconomic development [30]. The independent variables in the geographical detector model must use discrete quantities, therefore we have to classify the factors. According to the actual situation of the TRHR, slope was divided into 7 categories according to the Technical Regulations for Land Use Status Survey; aspect was divided into 9 categories according to slope orientation; Soil type was divided into 10 categories according to the traditional "Soil Occurrence Classification" system; vegetation type was divided into 9 categories according to the 1:1,000,000 Chinese Vegetation Atlas; landform type was divided into 6 categories according to the 1:1,000,000 Landform Atlas of the People's Republic of China; land use type was divided into 6 categories according to the 1:1,000,000 Land Use Map of China; the elevation, annual mean temperature and annual precipitation were divided into 9 categories according to the natural breakpoint method [31], and the GDP and population density were divided into 7 categories according to the natural breakpoint method [31] ( Figure 2).

Linear Regression Analysis
Linear regression analysis can analyze the trend of each raster in an image [32]. The raster calculator of ArcGIS was used to analyze the NDVI trend of each image element in the TRHR from 1982 to 2015, and categorized the NDVI change trend into seven classes according to the natural breakpoint method [31]: significant degradation, moderate degradation, slight degradation, basically unchanged, slight improvement, moderate improvement, and significant improvement. The slope can be calculated through Equation (1) [32]: In Equation (1): is the total number of the year series ( = 34 in this study), ranges from 1 to n, is the NDVI value of the th year, and is the variation trend of the NDVI; if > 0, the vegetation coverage shows an increasing trend; if

Data and Processing
The data included in this study included natural factors, such as the NDVI, digital elevation model (DEM), climate data, landform type, soil type, and vegetation type, as well as anthropogenic factors, such as land use type, population density, and GDP in the TRHR. NDVI data were obtained from the Big Earth Data Platform for Three Poles, using GIMMS NDVI3g data with a spatial resolution of 8 km and a temporal resolution of 16 days [24,25], the NDVI images for each year from 1982 to 2015 were obtained by maximum value composite (MVC) [26]. Vegetation coverage was divided into five classes according to NDVI values: low coverage (≤0.2), medium-low coverage (0.2-0.4), medium coverage (0.4-0.6), medium-high coverage (0.6-0.8), and high coverage (>0.8). The annual mean temperature and annual precipitation data were obtained from daily standard meteorological data of 26 meteorological stations in and around the TRHR from 1982 to 2015 using the inverse distance weighting method. DEM data were GDEMV2 30 m resolution digital elevation data from the Geospatial Data Cloud of the Chinese Academy of Sciences (http://www.gscloud.cn/ (accessed on 3 August 2021)), and the elevation, slope, and aspect data were obtained from the DEM data. Other data were obtained from Remote Sens. 2021, 13, 4175 5 of 24 the Resource and Environmental Sciences Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/ (accessed on 3 August 2021)). All the above data were extracted according to the vector boundary of the TRHR [27] and were resampled to make them consistent with the 8-km NDVI data image. Using ArcGIS to create a fishnet, 5853 random sampling points were generated based on 8 × 8 km grids according to the scope of the TRHR, and the spatial attribute values were extracted.
We selected NDVI as the dependent variable and six categories of topography (slope, aspect, elevation), soil (soil type), vegetation (vegetation type), landform (landform type), climate (annual mean temperature, annual precipitation), and human activity (GDP, population density, land use type), a total of 11 representative and easily quantifiable factors, as independent variables. Precipitation and temperature are important factors affecting vegetation changes [21], elevation, slope, aspect, and landform type affect vegetation distribution by changing water and heat conditions [28]. Soil and vegetation types are important environmental elements for vegetation growth [29]; economic development affects ecological environment, land use type, GDP, and population density are indicators that can quantify changes in socioeconomic development [30]. The independent variables in the geographical detector model must use discrete quantities, therefore we have to classify the factors. According to the actual situation of the TRHR, slope was divided into 7 categories according to the Technical Regulations for Land Use Status Survey; aspect was divided into 9 categories according to slope orientation; Soil type was divided into 10 categories according to the traditional "Soil Occurrence Classification" system; vegetation type was divided into 9 categories according to the 1:1,000,000 Chinese Vegetation Atlas; landform type was divided into 6 categories according to the 1:1,000,000 Landform Atlas of the People's Republic of China; land use type was divided into 6 categories according to the 1:1,000,000 Land Use Map of China; the elevation, annual mean temperature and annual precipitation were divided into 9 categories according to the natural breakpoint method [31], and the GDP and population density were divided into 7 categories according to the natural breakpoint method [31] (Figure 2).

Linear Regression Analysis
Linear regression analysis can analyze the trend of each raster in an image [32]. The raster calculator of ArcGIS was used to analyze the NDVI trend of each image element in the TRHR from 1982 to 2015, and categorized the NDVI change trend into seven classes according to the natural breakpoint method [31]: significant degradation, moderate degradation, slight degradation, basically unchanged, slight improvement, moderate improvement, and significant improvement. The slope can be calculated through Equation (1) [32]: In Equation (1): n is the total number of the year series (n = 34 in this study), i ranges from 1 to n, NDV I i is the NDVI value of the ith year, and Slope is the variation trend of the NDVI; if Slope > 0, the vegetation coverage shows an increasing trend; if Slope < 0, the vegetation coverage shows a decreasing trend; if Slope = 0, there is no significant change in the vegetation coverage.

Mann-Kendall Test
The Mann-Kendall method is a nonparametric statistical test used to determine the significance of trends [33]. The change trend was significant when |Z| > Z 0.05 . In this study, the Mann-Kendall statistical test was used to test the mutation points of the NDVI. The significance level was set at 0.05. The intersection of UF and UB is the mutation point; if there is more than one intersection, it is not certain whether it is the mutation point, and further testing is needed [34].

Moving t-Test
The moving t-test was used to test for mutations by examining whether the difference between the means of the two sample groups was significant [35]. If the difference between the mean values of the two subsequences exceeded the significance level of p = 0.05, the mutation was considered to be present; otherwise, no mutation was considered to be present.

Geographical Detector
The geographical detector is a statistical method used to detect spatial heterogeneity and its driving factors [11]. We used a geographical detector to compare the spatial distribution of NDVI vegetation with the spatial distribution characteristics of the detection factors; if a factor drives the NDVI variation, then the spatial distribution of the NDVI will be similar to the spatial distribution of that factor. This method has been successfully used to study the drivers of NDVI change [12][13][14][15][16][17].
(1) Factor detector. The factor detector q-statistic measures the degree of spatial stratified heterogeneity of a variable Y; and the determinant power of an explanatory variable X of Y. A factor detector is used to detect the spatial heterogeneity of the NDVI and the explanatory power of the independent variable X on the dependent variable Y, expressed by the q value [11]: In Equations (2) and (3): q is the explanatory power of the independent variable X on the dependent variable Y, with a value range of [0, 1]; the larger the q value, the more obvious the spatial heterogeneity and the stronger the explanatory power of X on Y. The study area is divided into h = 1, 2 . . . , L regions; N h and N are the number of units in layer h and the whole region, respectively; σ h 2 and σ 2 are the variances of the Y values of layer h and region, respectively; SSW and SST are the sum of variance within layer and total variance of region, respectively.
In this study, the independent variable X represents the detection factor Xs (s = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11), as is shown in Table 1, and the dependent variable Y is the NDVI. (2) Interaction detector. The interaction detector reveals whether the factors X1 and X2 (and more X) have an interactive influence on a response variable Y. Because the factors in nature do not exist independently, there are interactions between the factors, and the interactions between the factors need to be analyzed in the study. An interaction detector Remote Sens. 2021, 13, 4175 7 of 24 was used to detect the interaction between NDVI detection factors. It can detect any relationship between factors as long as they exist. The assessment methods used are presented in Table 2. Table 2. Types of interaction.

Foundation
Interaction Bivariate enhancement q (X1∩X2) = q (X1) + q (X2) Independent q (X1∩X2) > q (X1) + q (X2) Nonlinear enhancement (3) Risk detector. A risk detector was used to compare whether there was a significant difference between the mean values of the dependent variables in the two regions. This study was used to detect the appropriate range or types of the driving NDVI factors. The t-statistic used was the following [11]: In Equation (4): − Y h denotes the attribute mean within subregion h, n h is the number of samples within subregion h, and Var denotes the variance. According to the null hypothesis H 0 : , if H 0 is rejected at confidence level α, it is considered that there is a significant difference in the attributed means between the two subregions.
(4) Ecological detector. An ecological detector was used to detect whether there was a significant difference in the influence of different factors on NDVI changes, as measured by the F-statistic [11]: In Equations (5) and (6): N X1 and N X2 denote the sample sizes of the two factors X1 and X2, respectively, SSW X1 and SSW X2 denote the sum of within-layer variances of the strata formed by X1 and X2, respectively, and L1 and L2 denote the number of strata of the variables X1 and X2, respectively. According to the null hypothesis H 0 : SSW X1 = SSW X2 , if H 0 is rejected at the significance level of α, this indicates that there is a significant difference in the effect of the two factors X1 and X2 on the spatial distribution of attribute Y.

Spatial and Temporal Variation Characteristics of the NDVI in the TRHR
The NDVI values in the TRHR were high in the southeast and low in the northwest ( Figure 3). Regions with low vegetation coverage were mainly distributed in the northwest, with most being low coverage grassland; regions with high vegetation coverage were mainly in the southeast, where the hydrothermal conditions were better, the elevation was relatively low, and the vegetation was mainly high coverage grassland and forest. were mainly in the southeast, where the hydrothermal conditions were better, the elevation was relatively low, and the vegetation was mainly high coverage grassland and forest. In this study, the annual mean NDVI values from 1982 to 2015 were selected to represent the annual vegetation coverage status in the TRHR. The change observed exhibited an increasing trend, which is consistent with the findings of Zhai et al. [36]. With an increase rate of 0.002/10 a, the mean NDVI value increased from 0.454 in 1982 to 0.458 in 2015, and the maximum (0.493) and minimum (0.430) NDVI values occurred in 2010 and 1995, respectively ( Figure 4). These results indicate that the vegetation coverage of the TRHR has been improving, but with small changes from 1982 to 2015. Due to overgrazing, the ecological degradation of the TRHR as serious and the vegetation coverage was low. After 2005, the vegetation coverage gradually increased due to the increase in artificial precipitation and the implementation of ecological projects, such as the return of grazing to grass.

Trend Analysis of NDVI Changes in the TRHR
In 1982 and 2015, high and medium-high vegetation coverage areas accounted for more than 33% and 31%, respectively, of the TRHR area, while low and medium-low areas accounted for approximately 42% and 41%, respectively, of the total TRHR area. From 1982 to 2015 the low and high vegetation coverage areas decreased, with the former type decreasing the most (by 3.08%). During the same period, the medium vegetation coverage area increased the most (by 2.78%; Table 3). In this study, the annual mean NDVI values from 1982 to 2015 were selected to represent the annual vegetation coverage status in the TRHR. The change observed exhibited an increasing trend, which is consistent with the findings of Zhai et al. [36]. With an increase rate of 0.002/10 a, the mean NDVI value increased from 0.454 in 1982 to 0.458 in 2015, and the maximum (0.493) and minimum (0.430) NDVI values occurred in 2010 and 1995, respectively ( Figure 4). These results indicate that the vegetation coverage of the TRHR has been improving, but with small changes from 1982 to 2015. Due to overgrazing, the ecological degradation of the TRHR as serious and the vegetation coverage was low. After 2005, the vegetation coverage gradually increased due to the increase in artificial precipitation and the implementation of ecological projects, such as the return of grazing to grass. were mainly in the southeast, where the hydrothermal conditions were better, the elevation was relatively low, and the vegetation was mainly high coverage grassland and forest. In this study, the annual mean NDVI values from 1982 to 2015 were selected to represent the annual vegetation coverage status in the TRHR. The change observed exhibited an increasing trend, which is consistent with the findings of Zhai et al. [36]. With an increase rate of 0.002/10 a, the mean NDVI value increased from 0.454 in 1982 to 0.458 in 2015, and the maximum (0.493) and minimum (0.430) NDVI values occurred in 2010 and 1995, respectively ( Figure 4). These results indicate that the vegetation coverage of the TRHR has been improving, but with small changes from 1982 to 2015. Due to overgrazing, the ecological degradation of the TRHR as serious and the vegetation coverage was low. After 2005, the vegetation coverage gradually increased due to the increase in artificial precipitation and the implementation of ecological projects, such as the return of grazing to grass.

Trend Analysis of NDVI Changes in the TRHR
In 1982 and 2015, high and medium-high vegetation coverage areas accounted for more than 33% and 31%, respectively, of the TRHR area, while low and medium-low areas accounted for approximately 42% and 41%, respectively, of the total TRHR area. From 1982 to 2015 the low and high vegetation coverage areas decreased, with the former type decreasing the most (by 3.08%). During the same period, the medium vegetation coverage area increased the most (by 2.78%; Table 3).

Trend Analysis of NDVI Changes in the TRHR
In 1982 and 2015, high and medium-high vegetation coverage areas accounted for more than 33% and 31%, respectively, of the TRHR area, while low and medium-low areas accounted for approximately 42% and 41%, respectively, of the total TRHR area. From 1982 to 2015 the low and high vegetation coverage areas decreased, with the former type decreasing the most (by 3.08%). During the same period, the medium vegetation coverage area increased the most (by 2.78%; Table 3). The NDVI transfer matrix of the TRHR showed that there was a transformation in the NDVI at all levels from 1982 to 2015 ( Table 3). The shifted-out areas were mainly medium-high vegetation coverage, which shifted mainly to medium vegetation coverage, and the shifted-in areas were mainly medium-low, medium, and medium-high vegetation coverage, with a significant increase in medium vegetation coverage and a substantial decrease in high vegetation coverage.
Although the trend of the NDVI value of the TRHR was increasing, it was still dominated by low, medium, and medium-high vegetation coverage, which all accounted for more than 25% of the area, while the high vegetation coverage area accounted for the smallest proportion and decreased significantly. Previously, the ecological environment was severely damaged, and the restoration was difficult and slow. Land use is still dominated by low-coverage grassland, thus, the status of the vegetation coverage of the TRHR was still not optimistic.
From the linear regression analysis, it was concluded that the vegetation coverage of the TRHR showed an increasing trend from 1982 to 2015 ( Figure 5), thereby indicating that the vegetation coverage of the TRHR gradually recovered. The area with the largest increase in vegetation coverage was mainly distributed in the west and north, covering a total of 14.5 × 10 4 km 2 and accounting for 37.86% of the total area; this area was mainly dominated by grassland, meadow, and alpine vegetation. The area with the largest decrease in vegetation coverage was mainly concentrated in the center and the south, covering a total of 12.6 × 10 4 km 2 and accounting for 32.87% of the total area. Areas with unchanged vegetation were distributed throughout the region ( Table 4). The NDVI change trend in the TRHR increased in the north and west and decreased in the south and center. The desert in the northeast of the TRHR has gradually transformed into grassland and meadow vegetation types [37]. The unused land in the Sanjiangyuan Ecological Protection Project area has been transformed into low-coverage grassland, and the area of high-coverage grassland has increased significantly (Table A1), therefore the implementation of ecological projects has significantly improved the vegetation coverage of Zhiduo, Qumalai and Mado counties in the north and northwest of the TRHR. The decrease in vegetation coverage in Yushu, Jiuzhi, and Banma counties in the south may be due to the decrease in precipitation.
The M-K test showed that none of the intersection points of the UF and UB exceeded the critical value. Significance test indicated that |Z| = 0.048 < Z 0.05 = 0.236, thereby indicating that the trend of the NDVI change in the TRHR was not significant, but had multiple intersection points ( Figure 6). Therefore, the mutation points needed to be further examined using the moving t-test. Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 23 The M-K test showed that none of the intersection points of the UF and UB exceeded the critical value. Significance test indicated that | |=0.048 < Z0.05=0.236, thereby indicating that the trend of the NDVI change in the TRHR was not significant, but had multiple intersection points ( Figure 6). Therefore, the mutation points needed to be further examined using the moving t-test. The moving t-test showed that 2008 was the mutation point of the NDVI (Figure 7), which experienced a decreasing trend before 2008 and a significantly increasing trend in 2008 according to the cumulative departure method ( Figure 8). Therefore, the combination of the M-K test, moving t-test, and cumulative departure method led to the conclusion that the NDVI of the TRHR was mutated in 2008.   The M-K test showed that none of the intersection points of the UF and UB exceeded the critical value. Significance test indicated that | |=0.048 < Z0.05=0.236, thereby indicating that the trend of the NDVI change in the TRHR was not significant, but had multiple intersection points ( Figure 6). Therefore, the mutation points needed to be further examined using the moving t-test. The moving t-test showed that 2008 was the mutation point of the NDVI (Figure 7), which experienced a decreasing trend before 2008 and a significantly increasing trend in 2008 according to the cumulative departure method (Figure 8). Therefore, the combination of the M-K test, moving t-test, and cumulative departure method led to the conclusion that the NDVI of the TRHR was mutated in 2008. The moving t-test showed that 2008 was the mutation point of the NDVI (Figure 7), which experienced a decreasing trend before 2008 and a significantly increasing trend in 2008 according to the cumulative departure method (Figure 8). Therefore, the combination of the M-K test, moving t-test, and cumulative departure method led to the conclusion that the NDVI of the TRHR was mutated in 2008.

Factor Detection
According to the q values of each factor obtained from the factor detector (Table 5) The q value of annual precipitation was the largest, with an explanatory power of 55%, which was much more influential than other factors; therefore, annual precipitation was the main driving factor of vegetation change in the TRHR, followed by annual mean temperature, vegetation type, and elevation, with an explanatory power of more than 30%; land use type, landform type, population density, soil type, and slope had an explanatory power of more than 10%; GDP and aspect had little explanatory power.
The q value of annual precipitation was the largest, with an explanatory power of 55%, which was much more influential than other factors; therefore, annual precipitation was the main driving factor of vegetation change in the TRHR, followed by annual mean temperature, vegetation type, and elevation, with an explanatory power of more than 30%; land use type, landform type, population density, soil type, and slope had an explanatory power of more than 10%; GDP and aspect had little explanatory power.

Factor Detection
According to the q values of each factor obtained from the factor detector (Table 5), the magnitude of the influence of each factor on the NDVI of the TRHR was as follows: annual precipitation (0.550) > annual mean temperature (0.463) > vegetation type (0.409) > elevation (0.350) > land use type (0.244) > landform type (0.216) > population density (0.204) > soil type (0.147) > slope (0.141) > GDP (0.088) > aspect (0.055). The q value of annual precipitation was the largest, with an explanatory power of 55%, which was much more influential than other factors; therefore, annual precipitation was the main driving factor of vegetation change in the TRHR, followed by annual mean temperature, vegetation type, and elevation, with an explanatory power of more than 30%; land use type, landform type, population density, soil type, and slope had an explanatory power of more than 10%; GDP and aspect had little explanatory power.

Ecological Detection
Ecological detection reflects whether there is a significant difference in the influence of each detection factor on the NDVI in the TRHR. The results showed that there were significant differences in the effects of all factors on NDVI, except for the effects of soil and slope, population, and landform on NDVI ( Table 6). The effects of annual precipitation on the NDVI in the TRHR were significantly different from those of the other factors. The factor detection showed that annual precipitation was the dominant driver of NDVI changes in the TRHR, and the results of ecological detection further proved that the effects of annual precipitation were stronger than those of other factors. The non-significant differences between the effects of soil and slope, population density, and landform on the NDVI indicated that both had little influence on vegetation. Table 6. Ecological detection of factors.
Note: Y indicates a significant difference in the influence of two factors on vegetation NDVI, while N indicates no significant difference (confidence level is 95%).

Interaction Detection
A single variable could not explain the spatial variation in the NDVI, and the synergistic effects of multiple natural and anthropogenic factors needed to be considered. The geographical detector can reveal the interactions among the factors and their effect on NDVI changes. The results showed that all factor interactions enhanced the influence of a single factor on the NDVI, showing a bivariate and non-linear enhancement effects. Among them, the interactions of aspect with elevation, annual mean temperature, GDP, and population density, the interactions of soil type with GDP and population density, and the interactions of GDP with population density and land use type showed a non-linear enhancement effect, and the interactions of other factors showed a bivariate enhancement effect. Among them, the q value of the interaction between annual precipitation and other factors was high, with the explanatory power reaching more than 58%. This was higher than the explanatory power of the single factor of annual precipitation on vegetation, whose q value of interacting with elevation, annual mean temperature, and vegetation type was the largest, reaching approximately 68% ( Table 7). The annual precipitation was the dominant factor influencing the NDVI changes in the TRHR, and the interaction between the annual precipitation and other factors could further increase its influence on the NDVI in the TRHR. Among other factors, the q value of the interaction among vegetation type, elevation, and annual precipitation was larger, reaching approximately 60%; although the influence of GDP and aspect on the NDVI was small, their interaction with other factors greatly increased their explanatory power of the NDVI.

Risk Detection
We used the risk detector to determine the range or types of factors suitable for vegetation growth (Table 8); the suitable range or types of factors is very important for vegetation growth, the larger the NDVI value, the better the vegetation growth. The results of the risk area detection can be applied to the ecological protection project of the TRHR. The suitable range or types of different factors can be combined with the spatial distribution of temperature, precipitation, and population density to increase the vegetation coverage.

Annual Precipitation
The spatial distribution of vegetation coverage in the TRHR was consistent with the spatial distribution pattern of annual precipitation. The annual precipitation was divided into nine subzones. The mean NDVI value generally increased with the increase in annual precipitation and peaked in the 578 to 708 mm range, thereby indicating that this range promoted vegetation growth. The results showed that the annual precipitation subzone 9 was significantly different from the other subzones, so that the vegetation coverage was best in the 578 to 708 mm annual precipitation range in the TRHR ( Table 9). The interaction detector showed that interaction with other factors can further enhance the influence of annual precipitation on the NDVI.

Annual Mean Temperature
The factor detector showed that the annual mean temperature also had an important influence on the NDVI in the TRHR. The annual mean temperature was divided into nine subzones. The mean NDVI value increased and then decreased with the increase in the annual mean temperature, and peaked in the 1.65 • C to 3.82 • C range. There were significant differences between the mean NDVI values in subzone 6 and other subzones ( Table 10). The interaction of annual mean temperature with other factors enhanced the effect of the former on the NDVI. Temperature changes can cause changes in other factors in the region, and within the temperature range suitable for vegetation growth, the higher the temperature, the better the vegetation coverage, beyond which vegetation growth will be inhibited.

Vegetation Type
The vegetation type had an important influence on the NDVI of the TRHR, and the interaction with other factors further enhanced its influence on the NDVI. Vegetation types were divided into nine subzones. The mean NDVI values peaked in the coniferous forest vegetation type. There was no significant difference among the mean NDVI values in vegetation type subzones 2, 3, and 4. There were significant differences between the coniferous forest vegetation type and other vegetation type subzones; the coniferous forest, broadleaf forest, and scrub vegetation covers were better (Table 11).

Elevation
Elevation affects the spatial distribution of natural elements and human activity. The elevation was divided into nine subzones. The mean NDVI value increased and then decreased with the elevation of the TRHR, and it was better in the 3446 to 3851 m range. There were significant differences between this elevation range and other elevation subzones (Table 12). At elevations higher than 3851 m, the NDVI decreased as the elevation increased.

Land Use Type
The land use types were divided into six subzones. The NDVI value peaked in construction land, with no significant difference from the value obtained in forest land, and with significant differences from other land use types; therefore, construction land and forest land had the best vegetation coverage. The main land use type in the TRHR was grassland, accounting for 68%, of which low-coverage grassland accounts for 38%, followed by unused land, water area and forest land, which accounted for 23%, 5%, and 4%, respectively, of the total area. The cropland, forest land, middle-coverage grassland, and low-coverage grassland areas in the TRHR decreased from 1980 to 2015, while the high-coverage grassland, water area, construction land, and unused land areas increased, with the low-coverage grassland area decreasing the most and the high-coverage grassland and unused land area increasing the most (Table A1). Both the forest and construction lands were small, but both were distributed in the middle-high and high vegetation coverage areas east and south of the study area, with better hydrothermal conditions; the construction land was affected by human activities and had more green vegetation, thus the NDVI values were higher there.

Synergistic Effects of Other Factors
The factor detector demonstrated that the single factors of landform type, soil type, slope, aspect, GDP, and population density had small effects on NDVI changes in the TRHR, but the interactions of these factors with others could enhance the effects on NDVI changes.
The landform types of the TRHR were diverse and affected the distribution of vegetation. The landform types were divided into six subzones. The mean NDVI values peaked in the medium-undulating mountains; there were significant differences between this and other landform types, thereby indicating that the vegetation coverage in the mediumundulating mountains was the best. The soil types were divided into 10 subzones. The mean NDVI value peaked in semi-leached soil; there were significant differences between the mean NDVI value in this soil type and other soil types. Therefore, semi-leached soil had the best vegetation coverage. Different slopes and aspects led to differences in climatic elements, and suitable slopes and aspects were conducive to vegetation growth. The slope was divided into seven subzones. The mean NDVI value increased and then decreased with slope increases, and peaked in the 35 • to 45 • range. The vegetation of this slope consisted mainly of scrubs and alpine meadows. There were no significant differences between the mean NDVI value of this slope and those of slope subzones 6, 5, and 7, while there were significant differences with other subzones. Therefore, the vegetation growth conditions were better in the slope range of >25 • . As shown by the q value (Table 5), aspect had a minimal effect on the NDVI. The aspect was divided into nine subzones. The mean NDVI value fluctuated little with aspect changes. The NDVI value of the eastern slope was the largest, with no significant NDVI differences between this and aspect subzones 2, 3, 8, and 9, and significant differences with the other aspect subzones. Therefore, the vegetation coverage of the northern, northeastern, eastern, western, and northwestern aspects was the best.
Among the anthropogenic factors, both GDP and population density had little influence. The GDP was divided into seven subzones. The NDVI value peaked at a GDP of 12 × 10 4 -37 × 10 4 yuan/km 2 and was not significantly different from that of the area with a GDP of 1.04 × 10 4 -2.42 × 10 4 yuan/km 2 ; therefore, the vegetation growth was good in both of these areas. The population density was divided into seven subzones. The largest NDVI value was observed in the area with a population density of 74.95-94.31 people/km 2 , with no significant differences with the area with a population density of 8.43-19.37 people/km 2 ; therefore, the vegetation coverage was optimal in both areas. The area with a population density in the 74.95-94.31 people/km 2 range was very small, accounting for only 0.02% of the total area, which may have led to inaccurate results. If this area is not considered in the analysis, the NDVI will increase and then decrease with increasing population density, with larger values in the range of 8.43 to 19.37 people/km 2 , this result would be more accurate.

Discussion
Global warming over the last few decades has led to changes in the regional environment. Under the influence of climate change and human activities, vegetation green has generally increased in China [38]; the NDVI has shown an increasing trend in northern China over the past 40 years [39]; the Qinghai-Tibet Plateau tends to become warm and wet, and the vegetation status has gradually improved [40]. This study showed that the NDVI of vegetation in the TRHR also showed an increasing trend from 1982 to 2015, which is consistent with the trend of the NDVI change in China and Qinghai-Tibet Plateau during this period.
In this study, four geographical detectors were used to quantify the main drivers of the NDVI in the TRHR and the interaction of the factors. In the following sections, we will discuss the effects of natural and anthropogenic factors separately.

Effects of Natural Factors
The Qinghai-Tibet Plateau is a sensitive area for climate change in China [41]. This study indicated that climate factors were the main drivers of the NDVI changes in the TRHR, which is consistent with the findings of Chen [42]. The factor detector showed that the q value of annual precipitation was the largest and was the dominant factor influencing NDVI changes in the TRHR, which is consistent with the findings of Zheng [43] and Xiong [44]. In contrast, Xu [45] and Zhu [46] considered temperature as the dominant factor influencing NDVI variation in the TRHR; the differences in the results may be attributed to the different time scales of the study or the different spatial resolutions of the NDVI used. The warming trend in the TRHR was greater than the Chinese, as well as global average during 1982-2015, and precipitation was lower compared to the global [20]. Extreme temperature increases, and extreme precipitation is relatively stable. The rapid increase in temperature and slow increase in precipitation in the TRHR has led to regional warming and drought [47], while studies have shown that precipitation is the main factor affecting changes in vegetation NDVI in arid and semi-arid alpine meadow and alpine grassland regions [48]. The M-K test showed that the annual precipitation in the TRHR changed abruptly in 2004 and 2006 ( Figure A1), and extreme drought events occurred frequently. In 2006, the TRHR suffered an extreme drought, and the growth of forage grasses was disrupted and the grassland ecosystem was damaged [49], resulting in a decrease in NDVI. Precipitation increased abruptly around 2007 [50], since the NDVI has a lag effect on precipitation [51], the NDVI increased abruptly from 2008 onwards. The influence of extreme precipitation events on NDVI in the Qinghai-Tibet Plateau region is more pronounced than that of extreme temperature events, indicating that vegetation is more sensitive to changes in precipitation. Extreme wetness would offset the negative effects caused by extreme drought, and extreme high temperature events occurring in May would stimulate vegetation growth, while extreme low temperatures would inhibit vegetation growth [52]. The effects of extreme climatic events on the vegetation of the TRHR need further study. The influence of temperature on the NDVI gradually decreased, while precipitation occupied a more dominant position [42]. The annual precipitation and annual mean temperature of the TRHR decreased from southeast to northwest. The increasing trend of temperature in the TRHR was significantly greater than that of precipitation, can lead to the warming and drying of the TRHR, which will inhibit vegetation growth [53]. Seasonally, precipitation in the TRHR increases in spring and winter, and in summer when the temperature rises, precipitation also increases [54]. The growing season of vegetation in the TRHR is from May to September, and the climate is conducive for vegetation growth. Water resources are closely related to vegetation, and vegetation changes interact with hydrological processes [55]. Changes in temperature and precipitation lead to changes in vegetation patterns, which can alter surface hydrological characteristics, which, in turn, can affect changes in vegetation coverage [56]. The artificial rainfall implemented by the ecological project of the TRHR has restored the vegetation coverage and the increase in precipitation is beneficial to the growth of vegetation, but the excessive precipitation may cause soil erosion [57], which will instead damage vegetation, therefore the artificial rainfall project should be implemented scientifically and consistently to promote the growth of vegetation in the TRHR.
The vegetation types of the TRHR were mainly alpine meadows and alpine grasslands. During 1982-2015, part of the desert vegetation was converted to grassland and meadow vegetation types, increasing the vegetation coverage. Coniferous forests were mainly distributed in the elevation range of 3446 to 3851 m, and natural environmental conditions were more suitable for vegetation growth. The medium-undulating mountains were mainly dominated by meadow and scrub vegetation types, which were distributed in the southern part of the TRHR, with sufficient hydrothermal conditions and relatively suitable elevation, which are favorable for vegetation growth. The soil is the basis for vegetation growth. The fertilizer retention capacity of semi-leached soil is high, and the semi-leached soil of the TRHR is mainly distributed in mountainous areas, which are favorable for vegetation growth. In this study, the influence of soil type on vegetation change was small, but the interaction with other factors could enhance this influence; for example, the interaction of soil type with temperature and precipitation had a higher influence on the NDVI than did soil by itself. Soil temperature has an important effect on vegetation growth [58].
Topography affects vegetation distribution by changing water and heat conditions [28]. According to Chen [59], the 3500 to 3800 m elevation range is relatively low and precipitation and temperature conditions are good, thus the NDVI value is the largest in this elevation range. In elevations higher than 3800 m, the natural conditions become worse as the elevation increases, thus the NDVI value decreases as well. Slope affects vegetation growth by changing surface runoff, and vegetation coverage generally decreases with increasing slope. However, in this study, the gentle slope was more influenced by human activities; the vegetation coverage was low, while, with increasing slope, human influence decreased and vegetation coverage was relatively high. Aspect affects light intensity, which, in turn, changes the hydrothermal conditions for vegetation growth. The sunny slope has strong light, less soil water content, less nutrient accumulation, and lower vegetation coverage, while the shady, semi-shady, and semi-sunny slopes have sufficient soil water and high nutrient content [60], which are suitable for vegetation growth.

Effects of Anthropogenic Factors
According to the results of factor detection, anthropogenic factors had little influence on the NDVI. However, the combination of anthropogenic with natural factors can increase the impact. The population density in the TRHR was relatively small, and economic development was slow. Land use type had the greatest influence on the NDVI among anthropogenic factors. Low-coverage grassland is mainly located in the northwest, where water resources are scarce and the altitude is high, while high-coverage grassland is mainly located in the southeast where water and heat conditions are better. From 1980 to 2015, the conversion area between unused land and grasslands is large, and most unused land is converted into low-cover grasslands, but overall the increase in the area of unused land is greater than the decrease. Due to increase in population, land for construction has expanded. Ecological protection projects have increased the area of waters and lakes and improved the condition of wetlands. Before 2000, overgrazing led to the degradation of grassland; therefore, although the grassland area was large in the TRHR, the NDVI value was low. The implementation of the Sanjiangyuan Ecological Project in 2005 resulted in the slight recovery of the grassland, but the effects were short-term [61]. The areas with the highest NDVI values under the influence of GDP and population density were all located in the northeastern part of the TRHR, which is relatively densely populated, vegetation is affected by human activities, and the population is usually distributed in areas with better vegetation coverage [62] which have good survival conditions. Such natural conditions are also suitable for the growth of vegetation, but the increase in population will also cause some damage to vegetation, and the NDVI of vegetation will decrease beyond a certain range of population numbers.
The main conclusion of this study is that, compared with natural factors, anthropogenic factors had less influence on the NDVI of the TRHR. Natural factors, especially climatic factors, dominated the changes in the NDVI in the vegetation of the TRHR. In the context of global climate change, climatic factors have a strong association with vegetation change [56]. This is also verified by this study. The TRHR is at a high altitude, the population is sparse and the area of cultivated land is small. The impact of human farming activities is small, and although there is a certain degree of grazing, the impact is minimal relative to the climate, so the study area in this paper is basically equivalent to an undisturbed area. Therefore, the impact of human activities on the vegetation of the TRHR is very limited. As the impact of anthropogenic factors is short-lived, ecological engineering needs to be implemented continuously. Effective interventions for the restoration of vegetation in the TRHR can be based on the appropriate range or types of factors or a combination of factors. Separating natural factors from anthropogenic factors and quantitatively studying the influence of factors on vegetation is important for the ecological protection and sustainable development of the TRHR, as well as the middle and lower reaches of the region.

Effectiveness, Limitations, and Future Directions
To the best of our knowledge, this study is the first to use a geographical detector to quantify the effects of natural and anthropogenic factors on vegetation activity and effectively distinguish between the effects of natural and anthropogenic factors on the NDVI in the TRHR. The natural environment of the TRHR is complex, diverse, and spatially heterogeneous. Previous studies on vegetation drivers have used correlation analysis, which assumes a linear relationship between the NDVI and drivers, whereas correlation studies have shown a nonlinear relationship [63]; in contrast to traditional methods, the geographical detector can quantify the non-linear effects of factors and their interactions on vegetation change, making it well suited for this study. We also made the selection of factors with reference to existing studies, and the factors selected in this paper have been shown to be effective many times [13,14,29,64,65], so that the factors selected can be non-independent and the geographical detector method selected for this study allows the analysis of interactions between factors that have been neglected by traditional methods. However, the independent variable input to the geographical detector consists of type quantities, thus the numerical quantities must be classified. This study was based on the natural break method of classifying independent variables, which has been applied before and proven to be effective [29]; different methods of classification can affect the results. To ensure the length and completeness of the time series, NDVI data with a spatial resolution of 8 km were used in this study, which may have had some influence on the results owing to the low data resolution. Although NDVI is currently considered to be the most effective indicator for detecting vegetation change, it has shortcomings, such as the NDVI can reach saturation in dense vegetation canopies, which may lead to inaccurate trends in areas of dense biomass, and the effect on soil background in low vegetation coverage areas is not addressed [66], which were did not consider in this paper. Additionally, the different time ranges of the selected data may lead to some differences in the results, for example, if the growing season data are selected for analysis, the spatial and temporal distribution pattern of the NDVI in the growing season is basically the same as that of the whole year, the influence results are opposite to the annual data, and the influence of temperature (0.458) is slightly greater than the influence of precipitation (0.448). Although there are some differences in the results, the influence of climate factor is still the largest and is the dominant factor of vegetation coverage change in the TRHR, and this main result is unchanged. Therefore, for further research, data resolution should be further improved, while classification methods also need further improvement. To obtain a more accurate result, future studies could use the improved enhanced vegetation index (EVI) for comparison. The effect of growing season climate change on vegetation NDVI also needs further study.

Conclusions
In this study, which was based on GIMMS-NDVI data from 1982 to 2015 and 11 detection factors from the same period, we analyzed the spatial and temporal variation characteristics of the NDVI in the TRHR using linear regression analysis, the Mann-Kendall test, and the moving t-test. We also analyzed its spatial heterogeneity and driving factors using a geographical detector, and determined the appropriate range or types of factors suitable for vegetation growth. The main conclusions of the study are as follows: (1) The NDVI distribution of the TRHR was high in the southeast and low in the northwest; the change had an increasing trend in the west and north and a decreasing trend in the center and south. The annual mean value of the NDVI from 1982 to 2015 generally followed a slow increasing trend with a growth rate of 0.002/10 a; regions with low and high vegetation coverage decreased, while other regions increased. The NDVI increased abruptly in 2008. Overall condition of the TRHR has been improving, but vegetation coverage remains poor.
(2) The magnitude of the influence of each factor on the NDVI was as follows: annual precipitation > annual mean temperature > vegetation type > elevation > land use type > landform type > population density > soil type > GDP > aspect. Among them, annual precipitation had an explanatory power of more than 50% and was the dominant factor influencing NDVI changes in the TRHR. The annual mean temperature, vegetation type, and elevation had an explanatory power of more than 30% and also explained NDVI changes well. Land use type, landform type, and population density had an explanatory power of more than 20%, while other factors had less explanatory power. Compared with the natural factors, the influence of anthropogenic factors on the NDVI of vegetation in the TRHR was smaller. Climatic factors were the main drivers of NDVI changes in the TRHR.
(3) Interactions of bivariate and non-linear enhancements among the NDVI factors were observed, and there were no factors with weakening and independent effects. The interactions of annual precipitation, elevation, mean annual temperature, and vegetation type enhanced the influence of the factors to the greatest extent. Although factors such as the GDP and aspect had small influence on the NDVI, their interaction with other factors greatly increased their explanatory power on the NDVI.
(4) We analyzed the NDVI changes in the TRHR from 1982 to 2015, revealed the natural and anthropogenic factors driving NDVI changes, and determined the appropriate range or types of factors, which is important for ecological conservation and the sustainable development of the TRHR.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.