Increasing Negative Impacts of Climatic Change and Anthropogenic Activities on Vegetation Variation on the Qinghai–Tibet Plateau during 1982–2019

Abstract

Climate change, combined with ever-increasing anthropogenic activities, has had significant impacts on the vegetation of the Qinghai–Tibet Plateau (QTP). This study quantitatively analyzed the impacts of climate change and human activities on vegetation variation on the QTP from 1982 to 2019 based on AVHRR NDVI data and the residual trend method. The main results were as follows: (1) From 1982 to 2000, the vegetation of the QTP had an obvious restoration process, whereby 67.8% of vegetation coverage areas experienced an increasing trend, while it had a large range of degradation during 2001–2019, especially in the central QTP. (2) The positive effect of climate change on the vegetation of the QTP decreased, and the negative impact increased. The area of positive impact decreased from 68.54% in 1982–2000 to 47.13% in 2001–2019, while the negative-impact area increased from 31.46% to 52.87%. (3) The area negatively affected by human activities increased from 57.68% in 1982–2000 to 79.46% in 2001–2019 and was mainly concentrated in the grassland of the central QTP. The findings of this study provide a scientific basis for vegetation restoration and management in the QTP region.

1. Introduction

Vegetation is one of the most important components of terrestrial ecosystems [1,2]. The Qinghai–Tibet Plateau (QTP), due to its unique climatic conditions and topography, has formed a variety of climatic zones and nurtured diverse vegetation types [3]. A large body of evidence shows that the vegetation of the QTP experienced substantial change due to global climate change [4,5,6]. Specifically, the vegetation of the QTP changed from restoration to stagnation or even degradation prior to and following 2000 [5,7]. Given the role of vegetation on the QTP in regulating regional ecological security and even global climate change, quantifying vegetation change on the QTP has elicited considerable interest from scientists and policymakers [8].

In quantifying the spatial–temporal dynamic changes in vegetation, the Normalized Difference Vegetation Index (NDVI) was developed as the most effective indicator for the identification of vegetation activity dynamics due to their sensibility to vegetation and short revisit intervals [9,10]. A large number of studies have explored the relationship between vegetation change and its influencing factors on the QTP based on NDVI data. For instance, Immerzeel et al. (2005) and Zhong et al. (2019) identified the relationships between vegetation variations and precipitation and found that the major significant greening was mainly caused by climatic factors in the QTP based on the System Pour l’Observation de la Terre (SPOT) NDVI [11,12]. Different studies have shown that the vegetation change on the QTP shows great differences in time and space. For example, Shen et al. (2015) demonstrated that the vegetation showed an increasing trend in the northeast and a decreasing trend in the southwestern QTP from 1982 to 2011 using three different NDVI datasets (Global Inventory Monitoring and Modeling System (GIMMS) NDVI, SPOT NDVI, and Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI) [13]. Wei et al. (2021) found that vegetation restoration was mainly located in the central and southwestern QTP and that degradation was distributed in the northern and eastern QTP from 1982 to 2015 based on GIMMS NDVI [14]. However, Zhou et al. (2007) concluded that an increasing trend appeared to the west, south, north, and southeast, whereas a decreasing trend presented along the southern plateau boundary and the central–eastern region of the QTP using the Global Inventory Monitoring and Modeling System (GIMMS) NDVI data from 1982 to 2002 [3]. The data at different spatial and temporal resolutions are likely to be the main reason for these different conclusions. These studies indicated that the vegetation of the QTP has undergone great changes, and the influencing factors of the changes have become a topic of heated discussion among scholars.

Climate change and anthropogenic activities are critical factors that induce changes in vegetation [15,16]. Climate change is primarily reflected through (increasingly dramatic) changes in temperatures and precipitation, which affect the photosynthesis and respiration of plants, thus influencing plant growth [17]. Anthropogenic activities have negative and positive effects on vegetation through urbanization, deforestation, and overgrazing, as well as agricultural management, afforestation, and ecological projects [18,19,20,21]. Currently, there are several methods used to separate the contribution of human activities and climate change to vegetation [22]. The regression analysis method, based on a simple linear relationship between vegetation and impact factors, is simple and easy to implement, but it is difficult to explain this relationship to deal with vegetation and the factors that affect vegetation change [2,14,23]. The vegetation-based model can accurately reveal the relationship between vegetation and driving factors, but there are uncertainties due to different data and study areas [24]. Yet, these limitations are being overcome by applying process-based ecosystem models driven by observed historical environmental variables [16,22]. The Lund–Potsdam–Jena Dynamic Global Vegetation Model (LPJ) is a process-based model that can simulate photosynthesis, plant growth, maintenance and regeneration losses, and vegetation structure based on a combination of plant physiological relations, and generalized empirically established functions and plant trait parameters [25]. In addition, the residual trend method is simple and effective, and it has been widely used for the automated detection of vegetation change in remotely sensed vegetation and climate datasets [26]. Combining the LPJ model and the residual trend method, Kong et al. (2020) quantified the relative contributions of human activities and climate change to forest change in the Yangtze River basin.

Therefore, based on climate and long-term series AVHRR NDVI datasets, etc., along with the LPJ model and the residual trend analyses method, this study aimed to: (1) analyze the spatiotemporal dynamics of vegetation variation in different climate zones of the QTP during the past 40 years; (2) reveal the main reasons for vegetation changes, especially vegetation degradation after the year 2000; (3) quantify the impacts of climate change and anthropogenic activities to vegetation variation on the QTP. This study can improve our understanding of the driving mechanisms behind vegetation changes in plateau regions while providing useful references for vegetation restoration on the QTP.

2. Materials and Methods

2.1. Study Area

The Qinghai–Tibet Plateau resides in Southwest China with a territorial area of 2.57 million km² [27]. The plateau has an average elevation of more than 4000 m above sea level. Furthermore, the Qinghai–Tibet Plateau is characterized by a typically cold and dry alpine climate, with a mean annual air temperature ranging between −3.1 and 4.4 °C and mean annual precipitation ranging from 103 to 694 mm [28]. As a result, the Tibetan Plateau has been called “The Third Pole” of the Earth, with 61% of the plateau area being covered by alpine grasslands (Figure 1). The division of the Qinghai-Tibet Plateau climatic zone is mainly based on the results of Zheng Du et al. [29] (

Table 1. Eco-climate zones on the QTP.

Zone	Arid/Humid Region	Eco-Climate Zone
H1: Plateau Sub-frigid Zone	B: Semi-humid region	H1B1: Guoluo–Nagqu alpine shrub meadow zone
	C: Semi-arid region	H1C1: South Qinghai alpine steppe zone
		H1C2: Qutang alpine steppe zone
	D: Arid region	H1D1: Kunlun Mountains alpine desert zone
H2: Plateau Temperate Zone	A/B: Humid/Semi-humid region	H2AB1: West Sichuan alpine coniferous forest
	C: Semi-arid region	H2C1: East Tibet–Qilian alpine coniferous/steppe
		H2C2: South Tibet alpine shrub steppe zone
	D: Arid region	H2D1: Qaidam basin desert zone
		H2D2: Kunlun Mountains north desert zone
		H2D3: Ngari Mountains desert zone
H3: Subtropical Broadleaf Evergreen Forest Zone

).

2.2. Data Sources

The CDR AVHRR NDVI data used in this study were obtained from the National Science Data Center (http://www.geodata.cn/, accessed on 14 July 2022). The original CDR AVHRR NDVI dataset with spatial and temporal resolutions of 0.05 degrees and daily scale, respectively, was acquired from the NOAA AVHRR Surface Reflectance product (1982 to 2020). Based on the CDR AVHRR NDVI product, the rgee and terra package in R was used to call Google Earth Engine for monthly maximum synthesis and cropping [10]. The data were stored as TIFF format raster, with a total of 468 bands, each corresponding to one month’s data from 1982 to 2020. This study further used the maximum synthesis method to extract the annual NDVI values for the QTP from 1982 to 2019 [30]. Further, pixels with an annual mean NDVI of fewer than 0.1 were defined as bare and sparsely vegetated regions, which were excluded from this study [6]. The GIMMS NDVI-3g dataset with spatial and temporal resolutions of 1/12 degrees and 15 days, respectively, and a total of 69 NetCDF-4 format files were acquired from the Global Inventory Monitoring and Modeling Studies group and derived from the AVHRR sensor (1 July 1982 to December 2015). The maximum value was extracted from the semi-monthly scale data using the maximum synthesis method as the NDVI value of the month. The monthly SPOT NDVI dataset with 1 km spatial resolution (1998–2019) was obtained from the Resource and Environment Science and Data Center (https://www.resdc.cn/, accessed on 12 January 2022). The monthly MODIS NDVI dataset with 1 km spatial resolution (2001–2019) was obtained from the National Earth System Science Data Center (http://www.geodata.cn/, accessed on 18 July 2022).

Precipitation and temperature data were obtained from the monthly surface precipitation dataset of 0.5° × 0.5° grid points (V2.0), whereas the monthly surface temperature dataset of 0.5° × 0.5° grid points (V2.0) was provided by the National Meteorological Information Center (http://data.cma.cn, accessed on 31 January 2020) for the 1982-to-2019 timeframe. The annual average temperature and annual precipitation were obtained by calculating the monthly temperature and precipitation data pixel by pixel. For this study, ArcGIS software was used to resample the annual average temperature and annual accumulated precipitation to the same resolution as NDVI for subsequent calculations.

The vegetation type data came from the 1:1 million Chinese vegetation classification system data of the Resources and Environment Science and Data Center (http://www.resdc.cn/, accessed on 12 September 2022). Population density data for 2000 and 2015 were obtained from the WorldPop dataset (https://www.worldpop.org, accessed on 10 May 2022) at a 1 km resolution.

2.3. Methodology

The methodology of this work is shown in Figure 2. Firstly, AVHRR NDVI was compared with the other three sets of NDVI datasets to verify its reliability for the QTP from 1982 to 2019. Theil–Sen and Mann–Kendall tests were used for the analysis of the spatiotemporal dynamics of vegetation. Secondly, the Pearson correlation coefficient method was applied to analyze the relationship between climate change and vegetation. Thirdly, the potential NDVI affected by climate was obtained using the LPJ model driven by meteorological data, and the contribution of climate change and human activities to vegetation change was quantified using the trend residual method. Finally, the impacts of climate change and human activities on vegetation change on the QTP were explored.

For long-term data analysis, the Theil–Sen trend analysis method was used to calculate the changing trends pixel by pixel. The Theil–Sen trend analysis method is a robust non-parametric estimation algorithm that can effectively eliminate the influence of outliers for a long-time-series analysis and is an improvement on the least square linear regression method [31]. The calculation formula is as follows:

$$\mathrm{slope}=\mathrm{median}\left(\frac{x_j-x_i}{j-i}\right),i<j\le n$$

(1)

where slope refers to the change rate of the factor; $i$ and $j$ refer to the corresponding years; $n$ refers to the length of the time series; and $x_i$ and $x_j$ refer to the values of the corresponding years. Generally, a slope equal to 0 means that the variable has not changed, but in the calculation, a slope equal to 0 may not exist, and the value with a very small slope means that the variable in this area has not changed, so we set the slope value of threshold interval [−0.0005–0.0005] for no change. A positive value of the slope (>0.0005) represents an increasing trend, while a negative value (<−0.0005) represents a decreasing trend.

The Theil–Sen trend analysis is typically combined with the Mann–Kendall test (a nonparametric statistical test method) to judge long-term series data trends. It was initially proposed by Mann in 1945 and further improved by Kendall and Sneyers. Its advantages are that the measured values do not need to follow a normal distribution, nor does the trend need to be linear, and it is not affected by missing values and outliers. It has been widely used in trend significance tests of long-time-series data [32,33,34]. Test statistic S was calculated using the following formula:

$$\mathrm{S}={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^{n}sign(x_j-x_i)$$

(2)

$$\mathrm{sign}\left(v\right)=\left\{\begin{array}{c}1 for v>0\\ 0 for v=0\\ -1 for v<0\end{array}\right.$$

(3)

The variance of S is:

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

(4)

Statistics Z is defined as:

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

(5)

where $\left|\mathrm{Z}\right|$ ≥ 1.96 refers to when the trend passes the 95% significance test.

The Pearson correlation coefficient method was used to analyze the relationship between the NDVI, and the temperature and precipitation. A correlation coefficient greater than 0 represents a positive correlation, and one less than 0 represents a negative correlation. The value of the correlation coefficient ranges from −1 to 1 [35]. The closer the absolute value is to 1, the closer the correlation between the two variables is. The p-value was used to test the significance of the correlation coefficient, and p < 0.05 represented a significant correlation coefficient at a 95% confidence level [16].

The cumulative curve analysis method was used to distinguish the variability in time-series data. The Cumulative Sum of Departures of Modulus Coefficient (CSDMC) was used in this paper to detect the NDVI variability in different climate zones. The CSDMC is expressed by the following formula:

$$R_i=Q_i/\overline{Q}\left(i=1, 2, 3, \dots ,n\right)$$

(6)

$$K_p={\sum }_{i=1}^p\left(R_i-1\right) \left(p=1, 2, 3, \dots ,n\right)$$

(7)

where $i$ refers to a time series of N years, $K_p$ is the cumulative sum of departures of CSDMC from 1 to p years and $Q_i$ is the annual NDVI. The periods with downward trends of $K_p$ (negative slope) represent intervals of an NDVI lower than average, while upward trends (positive slope) represent intervals of an NDVI higher than average [36]. A sudden change in the direction of the CUSUM indicates a sudden shift in the average. A period where the CUSUM chart follows a relatively straight path indicates a period where the average is relatively stable [16].

The Lund–Potsdam–Jena Dynamic Global Vegetation Model (LPJ) combines process-based, large-scale representations of terrestrial vegetation dynamics and land–atmosphere carbon and water exchanges in a modular framework [25]. This model can simulate photosynthesis, plant growth, maintenance and regeneration losses, fire disturbance, soil moisture, runoff, evapotranspiration, irrigation, and vegetation structure based on a combination of plant physiological relations and generalized empirically established functions and plant trait parameters. Vegetation in a grid cell is described in terms of the fractional coverage of populations of different plant functional types (PFTs). Ten plant functional types (PFTs) are differentiated by physiological, morphological, phenological, bioclimatic, and fire-response attributes. The potential NDVI was simulated at a 0.5° spatial resolution with the LPJ model. It was run from 1961 to 1980 at a pseudo-daily time step, preceded by a 1000-year spin-up period that brought the carbon pools and vegetation cover in equilibrium with climate. The input data of the LPJ model included monthly mean temperature, precipitation, cloud cover, wet days, annual CO₂ concentration, and soil texture. The output data included the vegetation carbon, NDVI, vegetation net primary productivity (NPP), leaf area index (LAI), evapotranspiration, etc. [16]. Since the data inputs of the model were all climate factors, the simulation results could be regarded as the NDVI only being affected by climate change [16].

The effects of human activities and climate change on NDVI changes and their relative contributions were calculated via a trend residual analysis [37,38]. The satellite AVHRR NDVI was defined as ${\mathrm{NDVI}}_{\mathrm{obs}}$, which represented the combined effects of human activities and climate change on NDVI changes. The impact of climate change on NDVI changes was defined as ${\mathrm{NDVI}}_{\mathrm{pre}}$, which was predicted based on the LPJ model. The effect of anthropogenic activity on NDVI changes was defined as ${\mathrm{NDVI}}_{\mathrm{ha}}$. ${\mathrm{NDVI}}_{\mathrm{ha}}$ was obtained by calculating the residual difference between ${\mathrm{NDVI}}_{\mathrm{obs}}$ and ${\mathrm{NDVI}}_{\mathrm{pre}}$ using the following formula [2,15]:

$${\mathrm{NDVI}}_{\mathrm{ha}}={\mathrm{NDVI}}_{\mathrm{obs}}-{\mathrm{NDVI}}_{\mathrm{pre}}$$

(8)

The ${\mathrm{NDVI}}_{\mathrm{pre}}$ and ${\mathrm{NDVI}}_{\mathrm{ha}}$ Sen slopes represented the NDVI variation trends under the influence of climate change and human activities, respectively. The positive slope value indicated that climate change or human activities could increase the NDVI, which could promote the restoration of vegetation. In the opposite case, it could lead to a decrease in the NDVI, which could contribute to vegetation degradation [2]. Details regarding the identification and contribution of driving factors of NDVI changes are shown in

Table 2. Identification and contribution of driving factors of NDVI changes.

Vegetation Change	Driving Factor	Identification		Contribution (%)
		Slope (NDVIpre 3)	Slope (NDVIha 5)	CC	HA
Degradation	CC 1 and HA 2	<0	<0	$\frac{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{pre}}\right)}{\mathrm{Slope}\left({{\mathrm{NDVI}}_{\mathrm{obs}}}^{ 4}\right)}$	$\frac{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{ha}}\right)}{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{obs}}\right)}$
	HA	>0	<0	0	100
	CC	<0	>0	100	0
Restoration	CC and HA	>0	>0	$\frac{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{pre}}\right)}{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{obs}}\right)}$	$\frac{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{ha}}\right)}{\mathrm{Slope}\left({\mathrm{NDVI}}_{\mathrm{obs}}\right)}$
	HA	<0	>0	0	100
	CC	>0	<0	100	0

¹ Climate change; ² Human activities; ³ Predicted NDVI values; ⁴ Observed NDVI values; ⁵ NDVI affected by human activities.

[2,23].

3. Results

3.1. Spatiotemporal Changes in NDVI in Different QTP Climate Zones

The comparison of the monthly scale time series of AVHRR NDVI, GIMMS NDVI, MODIS NDVI, and SPOT NDVI from 1982 to 2019 for the QTP is shown in Figure 3a. Except for the outliers and missing values in October, November, and December 1989 and 1994, the AVHRR NDVI of other years had the same trend as the monthly data of the other three sets of NDVI products. The monthly values of AVHRR NDVI and GIMMS NDVI were close to each other and slightly higher than those of SPOT NDVI and MODIS NDVI. Before 1998, the annual means of AVHRR NDVI and GIMMS NDVI were close, and during 1999–2005, AVHRR NDVI was significantly higher than the other three datasets. After 2006, the annual means of AVHRR NDVI and MODIS NDVI were in good agreement and were higher than those of GIMMS NDVI and SPOT NDVI (Figure 3b). The scatter plot shows that AVHRR NDVI was significantly correlated with the other three NDVI datasets. The correlation coefficient between AVHRR NDVI, and GIMMS NDVI and SPOT NDVI was 0.94, and the correlation coefficient between AVHRR NDVI and MODIS NDVI was 0.89 (Figure 3c–e). These results indicated that AVHRR NDVI had good consistency with GIMMS NDVI, MODIS NDVI, and SPOT NDVI in the QTP, and the time series was the longest among several datasets, which could be used for long-term vegetation change monitoring in the QTP.

The spatial distribution of the mean value of GIMMS NDVI from 1982 to 2019 (Figure 4a) indicated that there were significant spatial differences in the distribution of the vegetation NDVI on the QTP. The annual average NDVI was 0.36, which increased from the northwest (0.1) to the southeast (0.84). The NDVI of the QTP showed an increasing trend from 1982 to 2019 (0.0013 per annum). However, the NDVI slope decreased from 0.0012/a to 0.0002/a before and after 2000. From 1982 to 2019, the QTP area with a significant increase (54.8%) in the NDVI was larger than the area with a significant decrease (4%), with the significantly increased areas being mainly distributed across the central, eastern, and northern regions of the QTP, including H1B1, H1C1, H1C2, H2C1, H2C2, H2AB1, and H3, while the significantly decreased areas were primarily concentrated in arid regions, such as H1D1, H2D1, H2D2, and H2D3 (Figure 4b).

Before 2000, the NDVI of the QTP mainly experienced an increasing trend (67.8%), and the area with a significant increase accounted for 16.3% and was mainly distributed in the central QTP (H1C2, H1D1, H2C2). The decreasing trend of the NDVI was mainly distributed in the south of H2C1 of the QTP (14.7%) (Figure 4c). However, the NDVI trend of the QTP after 2000 underwent great changes compared with prior to 2000. The NDVI trend in more regions of the QTP showed a downward trend (decrease of 32.5%; significant decrease of 9%) after 2000 and was mainly distributed in the H1B1, H1C2, H2C1, H2C2, H2AB1, and H3 regions of the QTP. The increasing trend of the NDVI was mainly distributed in the north of the QTP and the south of H2AB1 (increase of 36.1%; significant increase of 5.5%) (Figure 4d).

3.2. Relationships between Vegetation and Climatic Factors in Different QTP Climate Zones

From 1982 to 2019, the NDVI of the QTP was mainly positively correlated with the temperature, accounting for 83.1% of the total area, and 43.3% of the total area was significantly positively correlated and mainly distributed in the eastern, central, and southern regions (H2C1, H1B1, H1C1, H2AB1, H3, and H1C2) of the QTP. The areas with a negative correlation between the NDVI and the temperature were mainly distributed in the arid region, accounting for 15.2%, and the significant negative correlation only accounted for 1.8% (Figure 5a). The NDVI was mainly positively correlated with precipitation, but its spatial distribution differed from that of the temperature. About 71.7% of the NDVI was positively correlated with precipitation, of which 19.6% was significantly positively correlated and mainly distributed in H1C2, H2C2, H1C1, H1B1, H2C1, and H2D1. The negative correlation between the NDVI and precipitation was mainly distributed in H3 and H2AB1, accounting for 26.3% (Figure 5b).

To assess the deviations of the impacts of climatic and other factors on the NDVI during 1982–2019, the cumulative curve of climate factors and the NDVI were calculated in six representative regions with significant changes in the NDVI (Figure 6). It can be seen that except for the H3 region, the cumulative temperature curves of the other five climate regions deviated from the NDVI after 1982, and the deviation was the largest around 2000 (Figure 6a,c,e,g,i,k). The cumulative curve of precipitation and NDVI showed consistent growth changes before 2000, while the cumulative curves showed great deviation after 2000 in the six regions (Figure 6b,d,f,h,j,l). The difference between them indicated that the relationship between the NDVI and climate factors has changed since 2000.

3.3. Quantitative Attribution Analysis of Vegetation Changes

Figure 7a depicts the spatial distribution of the mean value of the predicted NDVI of the QTP from 1982 to 2019. The predicted NDVI was underestimated compared with the remotely sensed AVHRR NDVI, but it could well reflect the spatial distribution pattern of the NDVI of the QTP. The predicted spatial distribution of the NDVI gradually increased from the northwest to the southeast, and the maximum NDVI in the H3 region was 0.67. Compared with the underestimation in other regions, there was an obvious overestimation in the H1C2 region (Figure 7a). The scatter-plot comparison of AVHRR NDVI with the simulated NDVI indicates that the simulated NDVI was in good agreement with AVHRR NDVI (R = 0.65, p < 0.05) (Figure 7b). Therefore, the potential NDVI obtained through the LPJ model can be applied in this study.

The driving factors behind vegetation NDVI changes were analyzed (Figure 8), and it was found that the effects of climate change were responsible for the restoration of vegetation in most regions of the QTP from 1982 to 2019 (Figure 8a). The QTP overall experienced vegetation restoration was dominated by climate change (45.2%; H1B1, H1C1, H1C2, and H2C2) from 1982 to 2019. A total of 36.4% of the combined effects areas led to the restoration of vegetation (H1C1, H2C1, H2AB1, and H2C2), whereas 0.1% led to its degradation. Human activities mainly caused vegetation degradation (H1D1, H2D1, H2D2, and H2D3) in the arid areas in the west and north of the QTP and promoted vegetation restoration in the H3 region in the south of the plateau (Figure 8a).

Before 2000, the vegetation gradually recovered in most areas of the QTP, of which the combined effects of climate change and human activities accounted for 39.8%, where the areas affected by climate change and human activities accounted for 32% and 6.3%, respectively. Approximately 21.9% of the QTP exhibited vegetation degradation, with 18.38% being impacted by human activities (Figure 8b). In general, anthropogenic activities had a narrow range of impacts on the vegetation of the QTP, which for the most part promote vegetation restoration in the western and southern areas (H2D3, H2AB1, and H3), but induced the degradation of vegetation in the central and northeastern portions of the QTP (H2C1, H1C1, and H1B1).

Following 2000, the changes in the vegetation of the QTP were primarily characterized by restoration (3.5%) in part of the northeastern and southern areas (H2C1, H2D3, and H3) and degradation (45.7%) in the eastern, central, and western portions of the QTP (H1B1, H2AB1, H1C1, H1C2, and H2C2), dominated by human activities (Figure 8c). In addition to the anthropogenic impacts, the vegetation restoration (30%) areas on the QTP affected by climate changes were larger than those driven by human activities.

The contributions of climate change and human activities to vegetation changes on the QTP are shown in Figure 9 and Figure 10. From 1982 to 2019, the average contribution rate of climate change to vegetation restoration in each climate zone exceeded 20%, and in H1B1, H1C1, H1C2, H2C1, H2C3, and H2AB1, it exceeded 50%; the area with a contribution rate greater than 80% accounted for 49.14%. The contribution rate of climate change to vegetation degradation (<0%) was mainly concentrated in the range of −20–0%, accounting for 19.14% of the area, mainly distributed in H1D1, H2D1, H2D3, and H3 (Figure 9a,g and Figure 10). In contrast, the positive contributions of human activities were mainly distributed across the eastern, central, and northern regions of the QTP, in which the contribution rates of 20–40%, 40–60%, 60–80%, and >80% were 6.82%, 8.02%, and 8.1%, respectively, where contribution rates of >80% accounted for 9.03%. The anthropogenic influences (<0%) on the degradation of the vegetation of the QTP were mainly distributed across H1C1, H1C1, H1C2, and H2D1, which accounted for 62.06% (Figure 9d,g and Figure 10).

Compared with 1982–2019, the spatial distribution of the positive contribution of climate change to QTP vegetation before 2000 was reduced, and the average contribution rate also decreased. The positive contributions of climate change to QTP vegetation were primarily concentrated from 0 to 20% and >80%, which accounted for 16.42% and 34.49%, respectively, and were mainly distributed across H1C2, H1C1, H1B1, H2C2, H2AB1, and H2C1 (Figure 9b,h and Figure 10). The negative contribution of human activities to vegetation change was mainly distributed in the H1C2, H1C1, H1B1, H2AB1, and H2C1 regions, and the negative contribution was mainly −20–0% and <−40%, accounting for 38.14% and 19.36%, respectively. The positive contribution of human activities (42.32%) was mainly distributed across H3, H2D3, H1D1, H2D2, H1C1, and H1C2, whereas >80% accounted for 20.43%. Judging from the average value of human-activity contribution to vegetation change for each climate zone, human activities gave a positive contribution to each climate zone, except for H1C2 and H2D1, of which H1D1 had the highest average contribution rate of 61.64% (Figure 9e,h and Figure 10).

The positive contribution range of climate change to QTP vegetation shrank, and the negative contribution range expanded after 2000 (Figure 9c,i and Figure 10). Human activities played a major role in the degradation of QTP vegetation after 2000, especially H1C2, H1B1, H1C1, H2C1, and H2AB1, where the area with a contribution rate of <−40% accounted for 46.19% (Figure 9f,i and Figure 10).

In general, the positive effect of climate change on the vegetation of the QTP decreased, and the negative effect increased. Prior to and following 2000, the area proportion of the positive contribution of climate change to the QTP decreased from 68.54% to 47.13%, while the negative contribution increased from 31.46% to 52.87%. The negative effect of human activities on QTP vegetation change increased. The proportion of areas with positive contributions from human activities decreased from 42.32% to 20.54%, while the area with negative contributions increased from 57.68% to 79.46%, and especially, the area with a contribution rate < −40% increased from 19.36% to 46.19% and was mainly located in H1B1, H1C1, H1C2, H2C2, H2AB1, and H3.

The average contribution rate of climate change to forest and grassland changes in different regions showed a decreasing trend, especially in the H1C2 (forests, 81.48% to 32.10%; grassland, 64.73% to 37.58%) and H2C2 regions (forests, 55.82% to 22.39%; grassland, 61.77% to 37.73%) following 2000 (Figure 11a). The negative contribution rate of human activities to forest and grassland changes in different regions significantly expanded, especially in the H1C2 and H1B1 regions, where the contribution rate changed by more than −50% (Figure 11b). Combined with the relatively large grassland areas in H1C2, H2C, and H1B1, the negative impact of human activities on the grasslands in these areas was more serious.

From the perspective of the spatial distribution of human population density on the Tibetan Plateau, it is mainly concentrated across the H2C2, H2C1, H1C1, and H2AB1 regions. In contrast to 2000, the population density in 2015 showed a mainly increasing trend in H1B1 and H2AB1, while in H2C1, it decreased (Figure 12a–c). In general, after 2000, the intensity of human activities significantly increased in the H2C1, H1B1, H2AB1, and H2C2 regions. The regions of increased population density on the QTP were very consistent with the area of vegetation degradation caused by human activities.

4. Discussion

Long-time-series NDVI datasets have been extensively employed to monitor global vegetation changes since 1981 [8]. This study initially verified the spatial and temporal distribution of the vegetation of the QTP based on AVHRR NDVI from 1982 to 2019. The results revealed that the annual average NDVI was 0.36 and that it increased from the northwest (0.1) to the southeast (0.84) of the QTP, which was consistent with earlier studies [39,40]. Most regions of the world experienced a greening (increasing)–browning (decreasing)–greening trend over the 1981–2004, 1995–2004, and 2005–2012 time periods, respectively [41]. Vegetation changes on the QTP also showed a similar trend in this study, but with spatial differences. The QTP NDVI change rates were 0.0013/a, 0.0012/a, and 0.0002/a over the 1982–2019, 1982–2000, and 2001–2019 time periods, respectively, which aligned with Ding et al.’s results [42]. In terms of spatial changes from 1982 to 2019, 54.8% of the QTP underwent a significant increasing trend (H1B1, H1C1, H1C2, H2C1, H2C2, H2AB1, and H3); however, the arid regions, such as H1D1, H2D1, H2D2, and H2D3, of the QTP underwent significant decreasing trends. Furthermore, following 2000, the NDVI of the QTP considerably changed compared to before 2000 (Figure 4d). This was particularly the case for the H1B1, H1C2, H2C1, H2C2, H2AB1, and H3 regions, which decreased significantly. These results were similar to the findings obtained by Shen and Zhang et al. using AVHRR and MODIS NDVI [13,43].

Climate change/global warming is recognized as one of the most potent drivers behind the modification of vegetation, where a warming climate can be beneficial for vegetation growth, within certain parameters [44]. Over 28% of global greenness is attributed to climate change, such as anthropogenic warming and regional trends in precipitation; however, the contributions of climate change to vegetation greenness, or drought, can vary dramatically between regions worldwide [45]. Such variability was also evident on the QTP due to its unique orographic and topographic features, and various climate zones [4]. According to correlations between the NDVI, temperature, and precipitation (Figure 5), it could be inferred that vegetation changes were primarily impacted by temperature, whereas precipitation played a relatively minor role. Warming had a positive effect on improving the growth of vegetation in the semi-humid and semi-arid regions (H1C1, H1C2, H1B1, H2C1, and H2AB1), which also confirmed that the vegetation in this area underwent an overall growth trend from 1982–2019. These conclusions were consistent with the results obtained by Zhou, Piao, et al. [3,5]. It is possible that climate warming led to enhanced plant metabolism, advanced spring phenology, a prolonged growing season, and accelerated photosynthesis, all of which promoted vegetation growth [46,47,48].

It is worth noting that vegetation changes in the H1C2 region showed a decreasing trend after 2000, which may have been due to the QTP undergoing a warming slowdown during a so-called warming hiatus since 1999 [12]. Further, a regional drought caused by continuous warming weakened the NDVI responses to temperature changes [49]. Meanwhile, a hotter climate and increased solar radiation greatly increased vegetation respiration and evapotranspiration, thereby exacerbating soil moisture deficits [1,50,51]. In addition, a warmer climate significantly modified the permafrost of this region, which resulted in the expansion of the permafrost thawing area, the thickening of the seasonal thaw layer, and even the complete disappearance of the underlying permafrost layer. Permafrost degradation reduces soil moisture in the vegetation root layer, dries the topsoil, and reduces the availability of soil nutrients, which causes vegetation to degrade [52,53]. The H1C2 vegetation type is grassland, which is extremely sensitive to climate change. The warm and dry climate and the reduction in topsoil water were undoubtedly among the reasons for grassland degradation [54]. Changes in temperature and precipitation induced glacier mass loss, particularly after 2000, which led to the continuous expansion of the lakes in the H1C2 region, where they are widely distributed [55,56]. Although the increased lake area indirectly affected vegetation growth, it also directly caused the surrounding grassland to be submerged [39,57]. This also meant that the more expansive lake area had a stronger impact on the H1C2 region where the grassland type is the mainstay. Cheng et al. demonstrated that many grasslands and resident facilities were submerged in the Qiangtang Lake basin due to the expansion of lakes in the inner flow area of the QTP [54]. Another interesting study demonstrated that warmer temperatures (especially warm in winter) alone with higher moisture levels in the soil trigger vegetation to sprout too early when conditions are not ideal for healthy growth, which ultimately makes it fail and leads to a decline in the overall species diversity [58]. Increases in the tree and forest lines driven by elevated temperatures also contributed to the decrease in the NDVI in the H2AB1 and H3 regions [59]. In addition to climate change impacts, population and livestock increases and soil disturbances by small mammals are also potential causes of vegetation degradation in the region [60,61].

Vegetation is affected by both climate and anthropogenic factors with spatially variable effects [38,41]. The results of this study demonstrated that climate change mainly promoted vegetation change on the QTP, but prior to and following 2000, the area proportion of the positive contribution of climate change to the QTP decreased from 68.54% to 47.13%. The proportion of areas with positive contributions from human activities decreased from 42.32% to 20.54%, while the area with negative contributions increased from 57.68% to 79.46%, and especially, the area with a contribution rate < −40% increased from 19.36 to 46.19%. Vegetation changes induced by human activities were greater than those brought about by climate for different periods, which was consistent with previous findings [62]. This result was akin to the contributions of climate change to global vegetation shifts, where >28% of vegetation areas were mainly affected by climate change, albeit the contribution was only 8% [8,45]. Therefore, the anthropogenic effects on vegetation changes might be either positive or negative, contingent on the specific human activity [41]. The decrease in the positive impact of climate change and the increase in the negative impact of human activities resulted in the continuous degradation of grasslands in the central QTP.

The results of this study indicated that the impacts of human activities on vegetation changes had large disparities before and after 2000, mainly in the central–southern regions (H1B1, H1C1, H2AB1, H3, and H2C2) of the QTP. Prior to 2000, the human population distribution in these regions (H1C1, H1C2, and H3), at an altitude of over 4000 m, was relatively small, and the alpine meadow grew well under a suitable climate. The H2C1, H1B1, H2AB1, and H3 regions with relatively low altitudes were more susceptible to human activities [63]. However, in recent years, due to a continuously warming climate, the development of Western China, the opening of the Qinghai–Tibet Railway, the large-scale development of tourism, and the implementation of national development policies, the immigration and floating population in this region gradually increased, resulting in continuous population growth in the central and southern portions of the QTP (Figure 12) [64,65]. With population growth, the demand for dairy, meat, fur, and other animal products increased, with the income of herdsmen being directly related to the size of the livestock population. The stimulation of income demands led to the large-scale development of animal husbandry, which eventually contributed to the degradation of alpine grassland meadows in this region [61,65,66]. The research results obtained by Wei et al. showed that the population increased from 8.13 million in 1971 to 15.535 million in 2010, with an average annual growth of 185,100 individuals, and the total population increased nearly two-fold. The size of the livestock population on the QTP increased by 1.34 times, and the faster growth was primarily concentrated in the central and eastern regions [14,40]. Although the Grazing forbidden project had been implemented in Tibet, the impacts of the project have partially deteriorated due to activities such as farmland reclamation and abandonment [67,68]. Since 2000, the positive contributions of human activities to vegetation in the northeast portion of the QTP have been primarily due to these regions being covered by four ecological projects in China (Natural Forest Resources Protection Project, Three-North Shelterbelt Project, Grain for Green Project, and Grazing forbidden project) [69].

5. Conclusions

This study quantitatively analyzed the impacts of climate change and human activities on the vegetation of the QTP from 1982 to 2019 based on climate and remote sensing data, and residual trend analysis methods. The main results were as follows:

(1)

There was an obvious process of vegetation degradation on the QTP from 2001 to 2019 compared with 1982–2000. About 67.8% of vegetation coverage area experienced an increasing trend (significant increase of 16.3%), mainly distributed around the central part of the QTP prior to 2000. However, after 2000, the central region turned to a downward trend (decrease of 32.5%; significant decrease of 9%), and the area with a significant increase in vegetation decreased to 5.5%;

(2)

The positive effect of climate change on the vegetation of the QTP decreased, and the negative impact increased. The area of positive impact decreased from 68.54% in 1982–2000 to 47.13% in 2001–2019, and the decreasing area was mainly located in the central plateau. The negative-impact area increased from 31.46% to 52.87% and was mainly located in the H1B1, H1C1, H1C2, and H2C2 regions;

(3)

The decrease in the positive impact of climate change and the increase in the negative impact of human activities resulted in the continuous degradation of grasslands in the H1C2, H1B1, and H2C2 regions. The area negatively affected by human activities increased from 57.68% in 1982–2000 to 79.46% in 2001–2019, and especially, the area with a contribution rate < −40% increased from 19.36% to 46.19% and was mainly located in H1B1, H1C1, H1C2, H2C2, H2AB1, and H3. In particular, the contribution rate of human activities in H1C2 and H1B1 decreased by more than −50%. The findings of this study provide a scientific basis for vegetation restoration and management in the QTP region.