Analysis of Spatial and Temporal Changes in Vegetation Cover and Driving Forces on the Qinghai-Tibet Plateau

1. Introduction

In recent years, with the intensification of global climate change, the impact of climate change on terrestrial ecosystems and the corresponding feedback have become popular research topics. Vegetation, as an important component of ecosystems, links soil, the atmosphere and moisture through energy transfer and material transformation and can be indicative of global changes to a certain extent [1,2]. The normalized difference vegetation index (NDVI) is an important indicator used to characterize the growth status and quality of vegetation communities, and it is an effective indicator that can quantitatively characterize the growth status of vegetation [3]. At present, scholars around the world have performed extensive research on the characteristics of vegetation change and the corresponding response to climate change using the NDVI.

The Qinghai-Tibet Plateau is the largest plateau in the world and has the highest average altitude of any area on Earth, known as the "roof of the world" and the "third pole of the Earth"; its unique geographical location and complex natural environment make it the largest ecologically fragile area in China, as well as an area highly sensitive to climate change in Asia and even the Northern Hemisphere [4,5]. The Qinghai-Tibet Plateau, with its complex and diverse vegetation types, undulating terrain, and low and relatively concentrated disturbances by human activities, is an ideal site for studying vegetation cover changes and the corresponding response to climate change [6]. However, with the expansion of towns, population migration, livestock development and ecological construction on the Tibetan Plateau in recent years, the impact of human activities on vegetation on the plateau has been increasing [7], and large anthropogenic changes in land cover types may alter the way vegetation responds to climate variations [8]. Therefore, an increasing number of scholars have started to explore the response of vegetation to human activities and quantified and determined the relative effects of climate change and human activities on vegetation changes.

In previous studies on the relationship between the NDVI and various drivers, scholars [9,10,11] first mainly used correlation coefficients to reflect the spatial and temporal correlations between the NDVI and temperature and precipitation, and a few researchers considered the indirect effects of interactions between climate factors on vegetation changes [12,13]; in these studies, partial or complex correlation coefficients were used to analyze the direct and indirect effects of temperature and precipitation on vegetation changes, but no specific magnitude of the driving forces of vegetation change influenced by meteorological factors was obtained. Over time, multiple linear regression [14], artificial neural networks [15], support vector machine models [16], geographically weighted regression [17], residual analysis [18,19], and geographic probes [20,21] have been gradually applied to quantify the effects of temperature, precipitation, land use type, elevation, population, GDP, and other factors on the NDVI. Although previous studies identified the characteristics of vegetation cover changes on the Tibetan Plateau and the corresponding response to climate factors to a certain extent, the influence of climate factors on NDVI changes was not considered. Moreover, the influence of other factors, such as human activities and topographic factors, on vegetation has not been comprehensively studied.

In view of these research gaps, the spatial and temporal variation characteristics of the NDVI on the Tibetan Plateau from 1999 to 2019 are investigated in this paper based on SPOT/VEG NDVI data, land use data, topographic data and climate data from 1999 to 2019, supplemented by the Sen slope method and Mann–Kendall mutation test. The internal correlations between NDVI variation and precipitation and temperature were explored by time series correlation analysis, and the temporal response mechanisms between vegetation and climate factors in different ecosystems on the Tibetan Plateau were explored. Finally, the sensitivity of the annual mean NDVI to precipitation, temperature, elevation, slope, slope angle, soil type, and anthropogenic factors was measured by using the geographic probe method. The aim was to answer the following questions: (1) What were the dynamic spatial and temporal trends of vegetation cover on the Tibetan Plateau from 1999-2019? and (2) What is the degree of response of the annual mean NDVI to natural factors and human activities, and how do the interactions among different factors affect vegetation growth on the Tibetan Plateau?

2. Data Sources and Research Methods

2.1. Study Area Overview

The Qinghai-Tibet Plateau is in Southwest China, from the northern foot of the West Kunlun Mountain-Qilian Mountain Range in the north to the southern foot of the Himalayas and other mountain ranges in the south and from the Hengduan Mountains in the east to the Pamir Plateau in the west. The plateau covers an area from 26°00'12″- 39°46'50″N and 73°18'52″-104°46'59″E, with an average altitude of 4,400 m. The region spans six provinces, namely, Tibet, Xinjiang, Qinghai, Gansu, Yunnan and Sichuan, and covers a total area of approximately 2.57 million km², which is 26.8% of the total land area of the country [22]. Due to the unique geographical characteristics of this region, the climate is warm and humid in the southeast and dry and cold in the northwest, with a typical plateau climate, characterized by a dry and clean atmosphere, strong solar radiation, low temperatures, large diurnal temperature differences, small annual variations in major climate factors, low precipitation, and significant differences between seasons. The climate and topography of the Qinghai-Tibet Plateau are characterized by certain longitudinal geographical differences, and the temperature and precipitation generally decrease as the elevation gradually increases from east to west. The land use types in the region are grassland, unused land, and forestland in descending order of abundance, accounting for 48.5%, 32.0%, and 13.2% of the total area of the region, respectively, with water, arable land, and urban land mixed among them (Figure 1).

2.2. Data sources

2.2.1. NDVI Data

NDVI data were obtained from the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (http://www.resdc.cn). The annual vegetation index (NDVI) spatial distribution dataset in China was based on a continuous time series of SPOT/VEGETATION NDVI satellite remote sensing data. An annual vegetation index dataset since 1998 was generated using the maximum synthesis method [23]. The acquired NDVI data were processed to remove the effects of water, clouds, and aerosols, which effectively attenuated the influence of other factors on the data and improved the accuracy of the data. For the obtained NDVI data, simple preprocessing was performed, and the year-by-year vegetation index NDVI dataset was calculated based on mask-based cropping in ArcGIS and the vector boundary of the Tibetan Plateau using a raster calculator.

2.2.2. Meteorological Data

The temperature and precipitation data were obtained from the monthly average temperature and precipitation datasets for China from 1901.1-2022.12 provided by the National Earth System Science Data Center (http://www.geodata.cn), with a spatial resolution of approximately 1 km. This dataset was generated with the Delta spatial downscaling scheme in the Chinese region based on the global 0.5° climate data released by the CRU and the global high-resolution climate data released by WorldClim [24]. The data were validated with data from 496 independent meteorological observation points, and the validation results were credible. The dataset was preprocessed by using ArcGIS software to mask, crop, resample, and reproject the mean temperature and rainfall data based on the vector boundary of the Tibetan Plateau to obtain a raster image dataset of year-by-year mean temperature and rainfall with the same image size and coordinate system as the NDVI data. Then, MATLAB software was used to analyze the correlations among variables and significance of the raster data.

2.2.3. Resource Data

Land use resource data were provided by the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (http://www.resdc.cn), which provides a thematic database of Chinese national-scale multiperiod land use/land cover constructed based on manual visual interpretation at a resolution of 30 m using Landsat remote sensing imagery from the United States Landsat as the main information source [25]. Soil type data were obtained from the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (http://www.resdc.cn). These data were digitally generated from the 1:1 million Soil Map of the People's Republic of China prepared and published by the National Soil Census Office in 1995. DEM (digital elevation model) data were provided by the National Tibetan Plateau Science and Data Center (http://data.Tpdc.ac.cn), which generates 1 km data based on 1:250,000 contour and elevation points in China [26]. Population density and GDP data were obtained from the Data Center for Resource and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn). Notably, a kilometer-grid dataset of the spatial distributions of population and GDP in China at a spatial resolution of 1 km was obtained, and five periods of data for the years 2000, 2005, 2010, 2015, and 2019 were selected [27,28,29,30]. The above data were masked and cropped based on the vector boundary of the Tibetan Plateau using ArcGIS software and then resampled and reclassified to obtain the data required for the comprehensive geodetection analysis.

2.3. Research Methodology

2.3.1. Sen Trend Analysis

Theil-Sen median trend analysis is a nonparametric statistical method [31] with the advantages of not requiring the sample to follow a specific distribution, being computationally efficient and insensitive to outliers without a reduction in accuracy, and having the powerful ability to avoid measurement errors or errors associated with anomalous data. A comparison of different linear regression models shows that the method has significant advantages in the case of small samples. The corresponding equation is:

$$\mathsf{\beta }=\mathrm{M}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left\{\frac{{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{j}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}}}{\mathrm{j}-\mathrm{i}}\right\},\forall \mathrm{j}>\mathrm{i}$$

where β is the trend of vegetation change and NDVI_i and NDVI_j are the NDVI values at times i and j, respectively. β > 0 indicates an increasing NDVI trend, while the opposite relation indicates a decreasing NDVI trend.

2.3.2. Mann–Kendall Test

The Mann–Kendall (MK) test is a nonparametric test, which means that no prior assumptions about the statistical distribution of the data are required and can be used to verify the significance of the results of the Theil-Sen median estimation [32]. The normalized statistic Z is mainly used to test the trend and significance of time series. The Mann–Kendall test is used to determine the significance of the trend, and the related calculation is as follows:

$$\mathrm{Z}=\left\{\begin{array}{c}\frac{S-1}{\sqrt{var\left(S\right)}}, S>0\\ 0, S=0\\ \frac{S+1}{\sqrt{var\left(S\right)}}, S<0\end{array}\right.$$

$$\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}sgn({NDVI}_j-{NDVI}_i)$$

$$\mathrm{s}\mathrm{g}\mathrm{n}({\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{j}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}})\left\{\begin{array}{c}1, {\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{j}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}}>0\\ 0, {\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{j}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}}=0\\ -1, {\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{j}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}}<0\end{array}\right.$$

$$\mathrm{v}\mathrm{a}\mathrm{r}\left(\mathrm{S}\right)=\frac{n(n-1)(2n+5)}{18}$$

where n is the number of sample points in the time series. Sen+MK was combined to determine the trend of increasing or decreasing vegetation coverage during the study period. Based on the vegetation trend characterized by the β value and the results of the significance test of the vegetation trend, the statistical Z value at the significance levels of α = 0.05 and α = 0.01 was used as the threshold value, and the NDVI trend was classified into seven classes: highly significant decrease (β < 0, |Z| > 2.58), significant decrease (β < 0, 1.96 < |Z| ≤ 2.58), marginally significant decrease (β < 0, 1.65 < |Z| ≤ 1.96), no significant change (|Z| ≤ 1.65), marginally significant increase (β > 0, 1.65 < |Z| ≤ 1.96), significant increase (β > 0, 1.96 < |Z| ≤ 2.58), and highly significant increase (β > 0, |Z| > 2.58).

2.3.3. Pearson Correlation Analysis

To investigate the effects of climatic and anthropogenic factors on the NDVI, the correlation between the NDVI and trends in mean temperature and precipitation was calculated at the image scale [33] with the following equation:

$${\mathrm{A}}_{\mathrm{x}\mathrm{y}}=\frac{{\sum }_{i=1}^n(x_i - \overline{x})(y_i-\overline{y})}{\sqrt{{\sum }_{i=1}^n{(x_i - \overline{x})}^2}\sqrt{{\sum }_{i=1}^n{(y_i - \overline{y})}^2}}$$

where:

A_xy -- correlation of the trends of the two variables;

y_i -- the value of the NDVI in the i-th year;

x_i -- the values of average temperature and precipitation in the corresponding year;

?y and?x -- the mean values of the NDVI and different variables in the study period.

At test was used to test significance, and the correlation between the NDVI and climatic and anthropogenic factors was classified as highly significant (P < 0.01), significant (P < 0.05) or nonsignificant (P ≥ 0.05).

2.3.4. GeoDetector Calculation

Changes in the NDVI are governed by multiple natural factors, in addition to climate change, and are also influenced by a combination of human activities, topography and elevation [34]. To quantitatively evaluate and reveal the effects of natural and human factors on vegetation cover, a geoprobe approach is applied in this paper.

A geodetector is a statistical tool based on spatial statistics and spatial autocorrelation, and it can be used to explore spatial anisotropy, reveal the magnitude and significance of the effects of individual factors on a selected variable, detect risk zones and interaction strengths, and perform ecological detection [35]. In this study, we used a factor detection tool and an interaction detection tool to analyze the influence of NDVI drivers in the study region.

(1) Factor detection: The effects of the spatial heterogeneity of the dependent variable Y (NDVI values) and the magnitude of the spatial heterogeneity of the independent variable X (natural and socioeconomic factors) on Y, expressed as q, were calculated with the following equations:

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

$$\mathrm{S}\mathrm{S}\mathrm{W}=\sum _{h=1}^L{N}_h{\sigma }_h^2$$

$$\mathrm{S}\mathrm{S}\mathrm{W}=\mathrm{N}{\mathsf{\sigma }}^2$$

where h = 1, 2,... L, L is the classification of x or y; the value of q is in the range of [0, 1], and the larger the value of q is, the stronger the influence of x on the spatial variance of y. N_h and N denote the number of cells encompassed by the variable divided into h classes and the number of cells in the entire spatial region, respectively. σ_h² and σ² are the variances of h classes and Y in the region, respectively, and SSW and SST are the sum of the variances of L classes and the total variance of the region, respectively.

(2) Interaction detection: The interactions among different independent variables x are determined, and the result indicates whether the influence on y when the two factors act together is correlated or independent. The q value [q(X₁∩X₂)] is then output.

3. Results and Analysis

3.1. Characteristics of Temporal Changes in Vegetation Cover

Based on the continuous time series of SPOT/VEGETATION NDVI satellite remote sensing data, the maximum value synthesis method was used to generate the annual vegetation index data from 1999-2019, from which the fitted curve of the mean change in the annual NDVI was obtained. As shown in Figure 2, the annual average NDVI in the Tibetan Plateau region in the study period showed an overall fluctuating upward trend. Notably, the NDVI in the study area increased from 0.335 to 0.378 from 1999-2019, with a growth rate of 12.8%, or 0.0024-a^-1. The percentage of vegetation cover improved significantly and showed an increasing trend. The change in vegetation cover in the study area can be roughly divided into three stages. In the first stage, the NDVI increased rapidly from 0.335 to 0.345, with a growth rate of 3%, or 0.0036-a^-1, from 1999-2005. In the second stage, the NDVI increased slowly from 0.359 to 0.363 with a growth rate of 1.1%, or 0.0018-a^-1. In the third stage, from 2015-2019, the NDVI increased rapidly again from 0.363 to 0.378, with a growth rate of 4.1%, or 0.004-a^-1.

The trends from 1999 to 2019 observed in this study are consistent with the results of Mu Li et al. [36]. The changes in vegetation cover from 1999-2019 were closely related to the implementation of environmental protection policies and a general increase in people's awareness regarding environmental protection in the region. From 2015-2019, the provinces on the Qinghai-Tibet Plateau implemented a series of measures to protect the ecological environment in response to the requirements of the national ecological civilization construction policy, which led to significant changes in vegetation coverage in the region. Additionally, the temperature in the region warmed and rainfall increased during this period , which may have contributed to the NDVI increase.

Figure 2. NDVI, temperature and rainfall trends from 1999-2019 on the Tibetan Plateau.

Figure 2. NDVI, temperature and rainfall trends from 1999-2019 on the Tibetan Plateau.

3.2. Spatial Variation Characteristics of Vegetation Cover

To comprehensively evaluate the changes in spatial vegetation cover in the region, the NDVI was divided into low-density vegetation cover areas (0~0.3), medium-low-density vegetation cover areas (0.3~0.6), medium-density vegetation cover areas (0.6~0.8), and high-density vegetation cover areas (0.8~1). As shown in Figure 3, spatially, the vegetation cover in the region showed obvious regional differences, and the NDVI displayed a gradual decrease from southeast to northwest, which was consistent with the spatial distribution patterns of elevation(Figure 8c), rainfall(Figure 8b), and temperature(Figure 8a) [37] . The predicted changes in NDVI were most closely related to topographic and meteorological factors.

As shown in Table 1, the spatial distribution of vegetation cover in the region was determined every five years from 1999-2019, and the spatial and temporal changes in the spatial proportion of each vegetation coverage area in the region were analyzed. The proportion of low-density vegetation cover areas was the highest, and the change was significant, reaching 57.3% to 53.3%; moreover, the area proportion of this coverage class decreased at a rate of -0.206-a^-1. The proportions of medium-low density vegetation cover and medium-density vegetation cover were similar and decreased year by year. The change in medium-low density vegetation cover areas was the least significant, varying from 19% to 17.6%, with a rate of decrease of -0.068-a^-1. The change in the medium-density vegetation cover area was comparatively more significant, varying from 21.5% to 14.6%, with a rate of decrease of -0.384-a^-1. The high-density vegetation cover area increased from 2.1% to 14.6%, with the most significant increasing trend, reaching 0.668-a^-1.

Overall, although the proportion of the area characterized by low vegetation density decreased year by year, this classification is dominant in half of the study area. In contrast, the proportion of high-density vegetation cover area to the total study area is the lowest, but the growth trend is the most significant, indicating that although the vegetation cover in the area is not very high, the ecological conditions have had positive effects on vegetation.

Figure 3. NDVI coverage classes from 1999-2019 on the Tibetan Plateau.

Figure 3. NDVI coverage classes from 1999-2019 on the Tibetan Plateau.

Table 1. Temporal variations in the area share of NDVI coverage zones.

Table 1. Temporal variations in the area share of NDVI coverage zones.

3.3. Trend Analysis of Vegetation Cover Change

Theil-Sen median trend analysis was used to study the trend of vegetation cover changes in the region from 1999-2019. The effect of vegetation restoration in the region is remarkable, as shown in Figure 4a. The NDVI increase in the southeast and northeast areas of the region is most obvious, and the NDVI in the central and southwest areas is basically unchanged. Additionally, the NDVI in the northwest is decreasing. Using β=0 as the critical value of the vegetation trend, we reclassified the image element values to obtain Table 2, which shows that the NDVI in the whole region displays a clear increasing trend, and the areas with increasing trends and decreasing trends account for 67% and 33% of the total area, respectively.

Based on the MK test, the significance of the vegetation cover trend in the study area was assessed. As shown in Figure 4b, NDVI changes in the study area were dominated by nonsignificant changes, significant increases and highly significant increases. Most of the areas with slight significant increases and no significant changes were distributed in the central and western parts of the study area, and the areas with significant increases and highly significant increases were mainly distributed in the eastern part of the study area. Moreover, the areas with slight significant decreases, significant decreases and highly significant decreases were distributed along the northern part of the study area. The Z values of the statistics at the significance levels of α = 0.05 and α = 0.01 were used as the critical values, and the image element values were reclassified (Table 2). The results showed that the areas of very significant decreases, significant decreases, marginally significant decreases, no significant change, marginally significant increases, significant increases, and marginally significant increases accounted for 3.1%, 5.4%, 3.2%, 39.9%, 4.7%, 10.8%, and 33.0% of the total area of the region, respectively.

In general, there is obvious spatial variability in vegetation cover trends in the region, and these trends are mainly increasing. The areas that displayed a significant increase accounted for nearly 50% of the total area, which far exceeding the area of those with a significant decrease. This result indicated that the overall effect of Qinghai-Tibet Plateau ecosystem management is good.

Figure 4. NDVI trends and significance test results from 1999-2019 on the Qinghai-Tibet Plateau.

Figure 4. NDVI trends and significance test results from 1999-2019 on the Qinghai-Tibet Plateau.

Table 2. NDVI trends and significant trend proportions on the Tibetan Plateau from 1999-2019.

Table 2. NDVI trends and significant trend proportions on the Tibetan Plateau from 1999-2019.

3.4. Effect of Meteorological Factors on Vegetation Cover

According to Chen, J. H. et al[38], climate change is an important reason for the increase in the NDVI on the Tibetan Plateau, and temperature and precipitation are important indicators of climate change. Correlation analysis based on the image metric scale was used to explore the correlations between the NDVI and temperature and precipitation in space over time. As shown in Figure 5, the temperature showed a decreasing trend from 1999-2019, with the annual average temperature decreasing from -1.86°C to -2.06°C, a decrease of 10.8% or 0.0095-a^-1. Precipitation showed an increasing trend during the period, increasing from 32.46 mm to 33.33 mm, an increase of 2.7% or 0.041-a^-1. The spatial distribution of temperature in the region is shown in Figure 8a, with a gradual decrease from southeast to northwest, which is basically consistent with the spatial distribution pattern of precipitation in Figure 8b.

The results of the correlation analysis and significance test between the temperature and precipitation distributions and the NDVI are shown in Figure 6 and Figure 7. Then, the image values of the results were reclassified to obtain Table 3. An analysis of the results in Figure 6 and Table 3 showed that the percentage of the area with a positive correlation between the NDVI and mean temperature was 66.82%, mainly in the central-eastern part of the study area, which is consistent with the results of Zhao Qianqian [39]. The areas in the region in which the NDVI values showed highly significant and significant positive correlations with mean temperature, accounting for 2.73% and 7.97% of the total area of the region, respectively, were mainly concentrated in the eastern and central parts of the region, and insignificant correlations were observed in most other areas. This result indicated that there is high spatial heterogeneity in terms of the promotion of vegetation due to increasing average temperatures, and the correlation was not significant or of low significance. The effect of increased temperatures on the NDVI was most obvious in the eastern region, where the average temperature is relatively high. In Figure 7 and Table 3, the percentage of areas that displayed a positive correlation between the NDVI and mean annual precipitation was 66.29%, mainly in the northeastern and southwestern parts of the study area, which is basically consistent with the results of Wang Zhipeng [40]. The proportions of areas with highly significant and significant positive correlations were 4.87% and 9.55%, respectively, indicating that the NDVI and mean annual precipitation showed a simultaneous increase in most areas of the Tibetan Plateau during the study period. The spatial distributions of the NDVI and mean annual precipitation were highly and significantly correlated in the northeastern region, and the significance of trends in other regions was relatively weak. In general, the overall NDVI and precipitation correlation reached 14.4%, which was significantly higher than that of the NDVI and temperature correlation at 10.7%, indicating that precipitation has a greater influence on vegetation cover change on the Tibetan Plateau.

Figure 5. Temporal trends of rainfall and temperature on the Tibetan Plateau.

Figure 5. Temporal trends of rainfall and temperature on the Tibetan Plateau.

Figure 6. Spatial distribution of the NDVI and temperature correlations on the Qinghai-Tibet Plateau.

Figure 6. Spatial distribution of the NDVI and temperature correlations on the Qinghai-Tibet Plateau.

Figure 7. Spatial distribution of the NDVI and precipitation correlations on the Qinghai-Tibet Plateau.

Figure 7. Spatial distribution of the NDVI and precipitation correlations on the Qinghai-Tibet Plateau.

Table 3. Spatial distribution of NDVI correlation with temperature and precipitation in the Tibetan Plateau as a percentage of area.

Table 3. Spatial distribution of NDVI correlation with temperature and precipitation in the Tibetan Plateau as a percentage of area.

3.5. Analysis of the Influence of Each Driver on Vegetation Cover Change

The geodetector is an algorithm used for discrete data, so continuous variables (among the 8 independent variables in this paper, all are continuous variables except elevation, soil type, slope, and slope direction) are discretized. In this paper, the natural breakpoint method [41] was adopted, and all the driving factors used were divided into 9 categories, except for soil type, which was divided into 22 categories, as shown in Figure 8.

First, the magnitude of the single factor q-value was analyzed, and then the dominant driving factors were discerned. The driving factors in the analysis were series of data, such as temperature, precipitation, elevation, slope, slope angle, soil type, population density and GDP, among which temperature and precipitation were meteorological factors; elevation, slope, slope angle and soil type were topographic factors; and population density and GDP were anthropogenic factors. The results of the analysis based on the q-values of each driver in the region calculated with the geographic probe method are shown in Figure 9 (a), which indicates that the influence of each driver on the NDVI in the region ranked as follows: precipitation (0.709) > soil type (0.289) > elevation (0.252) > temperature (0.252) > slope (0.097) > slope direction (0.11) > population density (0.004) > GDP (0.003). The spatial variation in the NDVI in the region is the result of the combined effects of natural and human factors, among which precipitation, soil type, elevation, and temperature make strong contributions, with q-values greater than 0.25. Slope, slope direction, population density, and GDP make relatively small contributions. Among the factors that lead to spatial variations in the NDVI in the region, meteorological factors such as precipitation and temperature dominate, topographic factors such as soil type and elevation are the next-most important, and anthropogenic factors contribute the least to the NDVI.

To focus on the changes in vegetation cover considering the interactions between natural and anthropogenic factors, the relationships between different driving factors and the spatial variations in the NDVI were analyzed by using the geographic probe method, and the results are shown in Figure 9 b. Based on the comparison of the data in Figure 9, it can be seen that, we can see that (1) the combined effect of any two driving factors has a greater effect on the NDVI than the independent effect of a single factor, mainly showing two-factor enhancement and nonlinear enhancement trends. (2) The interaction between precipitation and other factors is significant, with q-values greater than 0.7, among which the interaction between precipitation and elevation is the largest, with a q-value of 0.797, followed by the interaction between temperature and precipitation, with a q-value of 0.793. Overall, the magnitude of the NDVI in the region is most obviously influenced by meteorological factors such as rainfall. (3) The interactions between anthropogenic factors and other factors have the smallest influence on the NDVI values in the region, except for the interaction with precipitation, which has a q-value of 0.7. The significance of the interaction with other factors is very low, and the q-values are all below 0.3. The results indicate that natural factors play a dominant role in influencing NDVI changes on the Tibetan Plateau, and anthropogenic factors have little effect, likely because of the high altitude and harsh climate in most areas of the Tibetan Plateau, the poor suitability of the living environment, the low and concentrated population density in the region, and the fact that most areas are uninhabited, resulting in the very limited influence of anthropogenic factors on the NDVI.

Figure 8. Drivers of vegetation changes on the Tibetan Plateau from 1999-2019: temperature(a), rainfall(b), elevation(c), slope(d), slope angle(e), soil type(f), population density(g), and GDP(h).

Figure 8. Drivers of vegetation changes on the Tibetan Plateau from 1999-2019: temperature(a), rainfall(b), elevation(c), slope(d), slope angle(e), soil type(f), population density(g), and GDP(h).

Figure 9. Interpretation of the driving forces of the drivers of NDVI change on the Tibetan Plateau.

Figure 9. Interpretation of the driving forces of the drivers of NDVI change on the Tibetan Plateau.

4. Conclusion and Discussion

4.1. Conclusion

(1) The NDVI on the Qinghai-Tibet Plateau rapidly increased from 1999 to 2019. Notably, 1999-2005, 2005-2015 and 2015-2019 were the three periods with the most obvious NDVI growth rates. The corresponding growth rates were 0.0036-a^-1, 0.0018-a^-1, and 0.004-a^-1, respectively; spatially, the NDVI pattern over the whole Tibetan Plateau was dominated by increases from northwest to southeast, and areas of growth accounted for 67% of the total area of the region.

(2) During the study period, the spatial variation in vegetation cover in the study area was large. Notably, the area share of low-vegetation-density cover decreased from 57.3% to 53.3%, and the rate of decrease was second only to that of medium-density-vegetation cover. However, low-vegetation-density cover remained the main type of vegetation cover in the region. In contrast, the proportion of high-density-vegetation cover in the total study area was the lowest, but it displayed the most significant growth in the study period, increasing from 2.1% to 14.6% of the area, a growth rate of 0.668-a^-1.

(3) The correlations between the NDVI and mean temperature and precipitation were mainly positive, indicating that the effects of temperature and rainfall on NDVI in the Qinghai-Tibet Plateau region are mainly positive. From the correlation test results, a correlation between the NDVI and temperature is only observed in 10.6% of the whole region, while the correlation between the NDVI and precipitation is more extensive, reaching 14.4%.

(4) The driving factors of NDVI changes in the past 21 years are precipitation>soil type>elevation>temperature>slope>slope direction>population density>GDP, and the main driving factors are precipitation, soil type, elevation and temperature, with q-values of 0.709, 0.289, 0.252 and 0.252, respectively. This finding suggests that natural factors such as precipitation, soil type, elevation and temperature dominate the spatial variation in the NDVI in the region, and anthropogenic factors such as population density and GDP have little influence on NDVI changes and environmental evolution.

(5) The analysis of the interactions among different drivers and their effects on the spatial variations in NDVI using the geodetector indicated that the combined effects of multiple drivers are greater than the effects of individual drivers, with two-factor and nonlinear enhancement trends. Specifically, the influence of meteorological factors on the NDVI values in the region is high, and the interactions between meteorological factors such as precipitation and topographic factors such as elevation have the most significant influence on the changes in NDVI values in the region, with a q value of 0.797. This value is significantly higher than that for the interactions between anthropogenic factors and other factors.

4.2. Discussion

In this paper, the spatial and temporal characteristics of vegetation cover changes are studied, and a quantitative analysis of the driving factors of vegetation on the Tibetan Plateau over the past 21 years is performed. The correlations between precipitation and temperature and the NDVI were analyzed, as were the levels of significance of these correlations, and the driving forces associated with relevant natural and anthropogenic factors were calculated using a geographic probe method. The main factors driving NDVI changes in the region were identified. Compared with existing studies of the Tibetan Plateau, this paper quantifies the driving forces of a series of factors influencing NDVI changes on the Tibetan Plateau with the help of a geo-probe method, and the driving force values are calculated for individual and multiple factors. Compared with the existing studies [42,43] of the spatial and temporal evolution trends of vegetation on the Tibetan Plateau and the corresponding influential factors, this paper quantifies the effects of topographic factors and human activities to determine specific drivers. In combination with the existing findings of Tibetan Plateau studies, a correlation analysis between the anthropogenic and natural factors that affect the NDVI on the Tibetan Plateau was performed, and it was found that precipitation displayed the strongest positive correlation with the NDVI and was the primary factor affecting vegetation growth and interannual variations in the region, which is consistent with the results of previous studies [44].

Generally, the correlation between the NDVI and drivers is studied using correlation, partial correlation and complex correlation coefficients. However, most studies that focused on the correlation between the NDVI and drivers only considered the correlation between the NDVI and meteorological factors (temperature, precipitation, etc.), and did not consider the spatial and temporal correlations between the NDVI and topographic factors (e.g., elevation and slope) or anthropogenic factors (e.g., population and GDP). In addition, the gray correlation [45] and Granger effect [46] methods, which are similar to the correlation coefficient approach, are both used to measure the relationship between the NDVI and other factors; only the measurement algorithms are different, and thus, the major limitations of these methods is the lack of quantifying the relationships between the NDVI and anthropogenic and topographic factors. The existing quantitative assessment methods also have various shortcomings; for example, the regression modeling method cannot yield the spatial distribution of or spatial differences in the contribution of each factor, and modeling methods based on biophysical processes can easily introduce large uncertainties due to the complex model structures and numerous model parameters; moreover, residual analysis cannot account for the differences in the climatic factors that affect vegetation growth in different ecosystems [47]. The geographic probe method can be used to separate the effects of natural and anthropogenic factors on vegetation change and calculate the corresponding effects on NDVI separately [48]; therefore, the magnitude of different driving forces of vegetation change can be assessed through a novel approach. However, quantitative analyses of the effects of anthropogenic and natural factors on the NDVI based on continuous, long time series are still lacking.

In general, a system for evaluating the relative roles of natural factors and human activities in the process of vegetation change considering the associated drivers is developed, but it still needs to be improved and modified for specific regions and influential factors to further enhance its accuracy and applicability.