Analyzing the Spatiotemporal Vegetation Dynamics and Their Responses to Climate Change along the Ya’an–Linzhi Section of the Sichuan–Tibet Railway

Abstract

Vegetation dynamics and their responses to climate change are of significant spatial and temporal heterogeneity. The Sichuan–Tibet Railway (STR) is a major construction project of the 14th Five-Year Plan for Economic and Social Development of the People’s Republic of China that is of great significance to promoting the social and economic development of Sichuan–Tibet areas. The planned railway line crosses areas with a complex geological condition and fragile ecological environment, where the regional vegetation dynamics are sensitive to climate change, topographic conditions and human activities. So, analyzing the vegetation variations in the complex vertical ecosystem and exploring their responses to hydrothermal factors are critical for providing technical support for the ecological program’s implementation along the route of the planned railway line. Based on MOD13Q1 Normalized Difference Vegetation Index (NDVI) data for the growing season (May to October) during 2001–2020, a Theil-Sen trend analysis, Mann–Kendall test, Hurst exponent analysis and partial correlation analysis were used to detect the vegetation dynamics, predict the vegetation sustainability, examine the relationship between vegetation change and hydrothermal factors, regionalize the driving forces for vegetation growth and explore the interannual variation pattern of driving factors. The growing season NDVI along the Ya’an–Linzhi section of the STR showed a marked rate of increase (0.0009/year) during the past 20 years, and the vegetation’s slight improvement areas accounted for the largest proportion (47.53%). Among the three hydrothermal parameters (temperature, precipitation and radiation), the correlation between vegetation growth and the temperature was the most significant, and the vegetation response to precipitation was the most immediate. The vegetation changes were affected by the combined impact of climatic and non-climatic factors, and the proportion of hydrothermal factors’ combined driving force slightly increased during the study period. Based on the Hurst exponent, the future vegetation sustainability of the area along the Ya’an–Linzhi section of the STR faces a risk of degradation, and more effective conservations should be implemented during the railway construction period to protect the regional ecological environment.

1. Introduction

Vegetation influences carbon exchange and the hydrological and biochemical cycles through photosynthesis, water transpiration and nutrient capture. Accordingly, vegetation is regarded as the most active factor in the Earth’s ecosystem, which plays an irreplaceable role in regulating climate change, balancing the carbon cycle and improving habitat quality [1,2,3,4,5,6,7,8,9]. Climate is widely recognized as the main abiotic driving factor for vegetation change [10,11,12]. The world has experienced significant warming over the last few decades [13,14,15,16,17,18], resulting in a general vegetation greening trend that profoundly influences the structural and functional adjustment of the ecosystem [19,20,21,22].

The Qinghai–Tibet Plateau (QTP) is of great sensitivity to global climate change and vulnerable to the ecological environment. The QTP plays an important role in dominating biodiversity maintenance, water and carbon supplementation and ecosystem services’ regulation [23,24,25,26,27]. Over the past few decades, with the global climate changing and human activity intensifying, the vegetation ecosystem of the QTP has experienced complex variation, which impacts the regional, domestic and even global ecosystems [5,6,7]. Numerous studies have indicated a greening trend on the QTP, mainly driven by climate change [1,28,29,30]. Appropriate warming is beneficial to vegetation improvement, while extreme temperature stress could weaken or even stop plant photosynthesis and respiration, to inhibit plant growth [31,32]. Meanwhile, suitable precipitation could directly alleviate the net reduction in plant carbon uptake in arid areas, resulting in a vegetation greening trend [33]. In the context of global climate change, exploring the spatial and temporal changing mechanisms of vegetation dynamics and their responses to climate change on the QTP is essential for regional ecological protection. The Normalized Vegetation Index (NDVI) directly characterizes the status of vegetation cover and indicates the possible future plant development trend, and it is often applied as the primary vegetation indicator that is easily influenced by human activities in the short term [17,34], mostly controlled by climatic factors such as temperature and precipitation in a long time series [35,36,37].

The planned STR, located in the southeast of the QTP, starts from Chengdu city in Sichuan province and travels westward all the way to Lhasa city, the capital of the Tibet Autonomous Region, with a total length of 1838 km. The railway spans the three river basins, runs through Mt. Hengduan and extends to the hinterland of the QTP, which is considered another vital traffic artery for the transportation of the QTP after the Qinghai–Tibet Railway (QTR). The STR stretches across the areas are characterized by a complex geological condition, significant elevation difference, variable topography, vertical climatic distribution and fragile ecological environment [38,39,40]. Although a few studies have illustrated that the vegetation growth on the QTP is mainly controlled by climatic factors [28,29,30], the Ya’an–Linzhi section of the STR is a transition zone from the Sichuan basin to the QTP with a significant elevation difference and a fragile ecological environment. Additionally, the railway is of great importance for promoting social and economic development, which may lead to regional degradation of vegetation and habitat quality during construction. So, it is necessary to explore whether hydrothermal parameters are still the primary driving forces of the special geomorphic and environmental conditions along the railroad. At present, against a background of global climate change, some studies on vegetation trends have focused on the areas with significant climatic and topographic differences [41,42]. However, since the areas along the Ya’an–Linzhi section of the STR have unique geological and ecological environments, the changing patterns and driving mechanisms of vegetation dynamics and habitat quality under the influence of climate change, increased human activities and railway construction are still unclear, while the technical support for regional vegetation conservation is insufficient.

In this study, we used MOD13Q1 NDVI data (16-day temporal resolution and 250 m spatial resolution) and meteorological data for 2001–2020 to analyze the trends of vegetation and hydrothermal factors, detect the vegetation response patterns to climate change with different elevations and vegetation types, regionalize the driving factors of vegetation change and predict the future vegetation sustainability along the Ya’an–Linzhi section of the STR. Overall, this study’s aims were (1) investigating the spatiotemporal variations of NDVI and hydrothermal factors, (2) exploring the changing pattern and driving mechanism of vegetation growth and (3) providing references for vegetation dynamic monitoring, ecological restoration and protection during the railway construction period.

2. Materials and Methods

2.1. Study Area

The study area is the Ya’an–Linzhi section of the STR, which contains four cities (Ya’an, Kangding, Changdu and Linzhi) and 13 counties (Lushan, Tianquan, Luding, Yajiang, Litang, Batang, Baiyu, Jiangda, Gongjue, Chaya, Basu, Bomi and Milin). The railway stretches across five geomorphic units (Sichuan basin, western Sichuan alpine canyon area, western Sichuan plateau basin area, southeast Tibet Mt. Hengduan area and southern Tibet valley area) and six large rivers (Dadu River, Yalong River, Jinsha River, Lancang River, Nujiang River and Yarlung Zangbo River) from east to west that characterize the unique ecological and geological conditions. Specifically, the area has topographically diverse terrain with more than 3000 m elevation differences (Figure 1). The rock types of the study area are significantly complex, including granite, gneiss, migmatite, slate, sandstone, greywacke, siltstone, mudstone, claystone and limestone. The distribution of faults is markedly extensive (Figure 2) [38,39,40]. Moreover, the Ya’an–Linzhi section of planned railway line crosses various climatic zones including the mid-subtropical humid climate zone, plateau temperate monsoon semi-humid climate zone and semi-arid monsoon climate zone, characterized by the strong wind, heavy rain, low average temperatures, a large temperature difference and strong ultraviolet radiation. The main vegetation types in the study area are needleleaf forest (NF), broadleaf forest (BF), shrub, meadow, alpine vegetation (AV) and cultural vegetation (CV), with a typical “three-dimensional zonation” distribution pattern [43] (Figure 3).

2.2. Data Source and Preprocessing

MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid data were collected from the National Aeronautics and Space Administration (NASA, https://lad-sweb.modaps.eosdis.nasa.gov/, accessed on 9 October 2021) and used to detect the vegetation variations in the study area. A total of 920 original images for 2001–2020 were acquired from the website (two tiles, h26v05 and h26v06, cover the study area with 23 images per year). The Modis Reprojection Tool (MRT) was used for original image format conversion and projection, and Savitzky–Golay filtering was applied for data smoothing and denoising to eliminate interferences of clouds, atmosphere and the solar altitude angle on remote sensing images [44,45]. Then, the Maximum Value Composite (MVC) method was applied for generating mean annual growing season NDVI data based on the values from May to October each year. To exclude sparse and non-vegetated areas, the pixels with a growing season NDVI less than 0.1 were masked out [46].

The monthly mean temperature, total precipitation and total radiation records at 38 meteorological stations in the Ya’an–Linzhi section of the STR and its neighboring areas were downloaded from the China Meteorological Data Service Center (http://data.cma.cn/, accessed on 16 October 2021) for 2001–2020. To obtain the monthly meteorological raster data with a 250 m spatial resolution, the Universal Kriging interpolation method in ArcGIS was used, which is a geostatistical procedure that can generate an estimated surface from a scattered set of points. The Kriging method was applied in the following steps: select Kriging type, method properties, semivariogram/covariance modeling, search the neighborhood and cross-validation. A key step of Kriging is choosing a suitable empirical semivariogram model from provided models including circular, spherical, exponential, Gaussian and linear, to predict attribute values at unsampled locations. Another key point of this method is cross-validation, by which the prediction errors can be examined to evaluate the overall interpolation accuracy. After that, the monthly climatic raster data were synthesized into the growing season total precipitation, mean temperature and average radiation between 2001 and 2020 using the Raster Calculator in ArcGIS for further trend and correlation analyses. The climatic trend changes were calculated and examined with Theil–Sen and Mann–Kendall test methods.

The Digital Elevation Model (DEM) was collected from the Geospatial Data Cloud (https://www.gscloud.cn/, accessed on 25 September 2021) with a 30 m spatial resolution. The fault distribution [47] and formation lithology [48] data for the study area were downloaded from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, accessed on 6 July 2022). The vegetation type distribution data for the study area were extracted from the 1:1,000,000 China vegetation map, which was published by the Resource and Environment Science and Data Center (https://www.resdc.cn/, accessed on 16 September 2021).

2.3. Methods

Based on the NDVI and climatic data, the variation trends of vegetation and climatic factors were tested with Theil-Sen’s slope and Mann-Kendall method. Then, the coefficient of variation and Hurt exponent were calculated to explore the volatility and sustainability of vegetation dynamics. Finally, the partial correlation coefficient method was used to analyze the vegetation response to climatic parameters (Figure 4).

2.3.1. Trend Analysis

Theil–Sen slope analysis, a robust nonparametric statistical method, was applied to detect the long-term trend of NDVI and three climatic parameters at the pixel scale. As a mathematical-statistical method, the advantage of this method is that it is unaffected by the sample distribution and a few outliers [49]. The calculation formula is as follows:

$$\mathrm{Slope}=\mathrm{Median}\left(\frac{{\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}}{\mathrm{j}-\mathrm{i}}\right)\hspace{1em}0<\mathrm{i}<\mathrm{j}<\mathrm{n}$$

(1)

where ${\mathrm{x}}_{\mathrm{i}}$ and ${\mathrm{x}}_{\mathrm{j}}$ represent the values at times $\mathrm{i}$ and $\mathrm{j}$, respectively. $\mathrm{n}$ is the length of the time series. $\mathrm{Slope}>0$ indicates an upward trend, while $\mathrm{Slope}<0$ represents a downward trend.

The Mann–Kendall (MK) test, which measures the significance of a trend, is a nonparametric test method commonly used in meteorology and hydrology studies [50]. It has the advantage that samples do not need to obey a certain distribution and is free from the interference of outliers [51,52], meaning it has been applied in vegetation variation studies in recent years. The Mann–Kendall statistic $\mathrm{S}$ is calculated for the time series as follows:

$$\mathrm{S}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}-1}{\displaystyle \sum }_{\mathrm{j}=\mathrm{i}+1}^{\mathrm{n}}\mathrm{sgn}\left({\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}\right)$$

(2)

where ${\mathrm{X}}_{\mathrm{i}}$ and ${ \mathrm{X}}_{\mathrm{j}}$ present the variable values of the pixels at time $\mathrm{i}$ and $\mathrm{j}$, respectively. n represents the length of the time series, and $\mathrm{sgn}$ is the sign function, defined as

$$\mathrm{sgn}\left({\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}\right)=\left\{\begin{array}{rr}1,& {\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}>0\\ 0,& {\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}=0\\ -1,& {\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}<0\end{array}\right.$$

(3)

When $\mathrm{n}\ge 8$, $\mathrm{S}$ follows a normal distribution with a mean of zero and a variance of:

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

(4)

The standardization of $\mathrm{S}$ is as follows:

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

(5)

$\mathrm{Z}$ is the statistic normalized by the Mann–Kendall test in the range of (−$\infty$, +$\infty$). In this study, the significance of NDVI and three climatic factors were tested at a confidence level of $\mathsf{\alpha }=0.05$. The results of the MK test were classified as significant ($\left|\mathrm{Z}\right|>1.96$) or insignificant changes ($\left|\mathrm{Z}\right|<1.96$). Based on the NDVI trend slope analysis and MK test results, the NDVI trend was divided into five classes (

Table 1. NDVI trend classification.

NDVI Trend	Slope	Z
Significant degradation	≤−0.0005	>1.96
Slight degradation	≤−0.0005	−1.96–1.96
Stable	−0.0005–0.0005	−1.96–1.96
Slight improvement	≥0.0005	−1.96–1.96
Significant improvement	≥0.0005	>1.96

).

2.3.2. Coefficient of Variation

The coefficient of variation is a statistical index that describes the variation degree of the time-series values, reflecting the discrete degree of the NDVI and the interannual volatility of vegetation [53]. Specifically, a lower $\mathrm{CV}$ represents the lower fluctuation and higher stability of the vegetation variation, while a higher $\mathrm{CV}$ indicates the higher volatility of vegetation. The calculation formula is as follows:

$$\mathrm{CV}=\frac{1}{\overline{\mathrm{NDVI}}}\sqrt{\frac{1}{\mathrm{n}-1}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}(\Delta {\mathrm{NDVI}}_{\mathrm{i}}-{\overline{\mathrm{NDVI}})}^2}$$

(6)

where ${\mathrm{NDVI}}_{\mathrm{i}}$ is the NDVI value in year $\mathrm{i}$, while $\overline{\mathrm{NDVI}}$ is the mean NDVI value of the growing season from 2001 to 2020. In this study, the $\mathrm{CV}$ of the NDVI was calculated pixel by pixel to indicate the vegetation fluctuation for the research area during 2001–2020. The $\mathrm{CV}$ values were divided into five classes: low volatility (CV < 0.05); relatively low volatility (0.05 ≤ CV < 0.1); moderate fluctuation (0.1 ≤ CV < 0.15); relatively high volatility (0.15 ≤ CV < 0.2); high volatility (CV ≥ 0.2).

2.3.3. Hurst Exponent

Hurst exponent and rescaled range (R/S) analysis has been developed as an effective method to quantitatively describe the long-term dependence of time series and distinguish the sustainability of time-series data [54,55]. This method has been widely used in the fields of hydrology, climatology, economy and geology, and has also been applied in the research on vegetation variation carried out in recent years [48]. The calculation procedures are as follows [56]:

• To define the time series:

$$\left\{\mathrm{NDVI}\left(\mathrm{t}\right)\right\}\hspace{1em}\mathrm{t}=1, 2, \dots , \mathrm{n}.$$

(7)

2.

To define the mean sequence of the time series:

$$\overline{{\mathrm{NDVI}}_{\left(\mathsf{\tau }\right)}}=\frac{1}{\mathsf{\tau }} {\displaystyle \sum }_{\mathrm{i}=1}^{\mathsf{\tau }}{\mathrm{NDVI}}_{\left(\mathrm{t}\right)}\hspace{1em}\left(\mathsf{\tau }=1,2, \dots , \mathrm{n}\right)$$

(8)

3.

To calculate the accumulated deviation:

$$\mathrm{X}\left(\mathrm{t},\mathsf{\tau }\right)={{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{t}}({\mathrm{NDVI}}_{\left(\mathrm{t}\right)}-\overline{{\mathrm{NDVI}}_{\left(\mathsf{\tau }\right)}})\hspace{1em}\left(1\le \mathrm{t}\le \mathsf{\tau }\right)$$

(9)

4.

To create the range sequence:

$$\mathrm{R}\left(\mathsf{\tau }\right)=\underset{1\le \mathrm{t}\le \mathsf{\tau }}{\max }\mathrm{X}\left(\mathrm{t},\mathsf{\tau }\right)-\underset{1\le \mathrm{t}\le \mathsf{\tau }}{\min }\mathrm{X}\left(\mathrm{t}, \mathsf{\tau }\right)\hspace{1em}\left(\mathsf{\tau }=1, 2, \dots , \mathrm{n}\right)$$

(10)

5.

To create the standard deviation sequence:

$$\mathrm{S}\left(\mathsf{\tau }\right)=\sqrt{\frac{1}{\mathsf{\tau }}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}(\Delta {\mathrm{NDVI}}_{\mathrm{i}}-\overline{\Delta {\mathrm{NDVI}}_{\left(\mathsf{\tau }\right)}})}\hspace{1em}\left(\mathsf{\tau }=1, 2, \dots , \mathrm{n}\right)$$

(11)

6.

To calculate the Hurst exponent:

$$\frac{\mathrm{R}\left(\mathsf{\tau }\right)}{\mathrm{S}\left(\mathsf{\tau }\right)}={\left(\mathrm{c}\mathsf{\tau }\right)}^{\mathrm{H}}$$

(12)

The least-squares method is used by fitting the equation $\log {\left(\mathrm{R}/\mathrm{S}\right)}_{\mathrm{n}}=\mathrm{a}+\mathrm{H}\times \log \left(\mathrm{n}\right)$ to obtain the Hurst exponent value (H). H is in the range of 0 to 1 and is divided into three types [54,55,57]; when H = 0.5, the future trend is a random series with no sustainability, indicating that the future time-series trend is not related to the current trend; when H > 0.5, the NDVI time series is sustainable and the future trend is consistent with the current trend; when H < 0.5, the future trend is uncertain and opposite to the current trend [57]. Based on the NDVI trend slope analysis and Hurst exponent results, the NDVI sustainability was divided into five classes (

Table 2. NDVI sustainability classification.

NDVI Sustainability Level	Slope	H
Continuous degradation	≤−0.0005	>0.5
No continuous degradation	≤−0.0005	<0.5
Stable	−0.0005–0.0005	0.5
No continuous improvement	≥0.0005	<0.5
Continuous improvement	≥0.0005	>0.5

).

2.3.4. Partial Correlation Analysis

The partial correlation coefficient (PCC) and complex correlation coefficient (CCC) between the NDVI and three hydrothermal factors were estimated to evaluate the relationship between vegetation variations and climate change [58]. In this paper, neither the NDVI nor climatic variables were detrended for the correlation analysis. Detrending may be more acceptable for time series that contain a deterministic trend, but it may introduce errors for time-series data that contain a stochastic trend [59,60]. Yet, it is difficult to differentiate between deterministic and stochastic trends for a relatively short sample period [61,62]. The calculation formulas are as follows:

$${\mathrm{R}}_{\mathrm{xy}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)\left({\mathrm{y}}_{\mathrm{i}}-\mathrm{y}\right)}{\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)}^2}\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{y}}_{\mathrm{i}}-\mathrm{y})}^2}}$$

(13)

$${\mathrm{R}}_{\mathrm{y}1·23}=\frac{{\mathrm{R}}_{\mathrm{y}1·2}-{\mathrm{R}}_{\mathrm{y}3·2}{\mathrm{R}}_{13·2}}{\sqrt{\left(1-{\mathrm{R}}_{\mathrm{y}3·2}^2\right)\left(1-{\mathrm{R}}_{13·2}^2\right)}}$$

(14)

$${\mathrm{R}}_{\mathrm{y}·123}=\sqrt{1-\left(1-{\mathrm{R}}_{\mathrm{y}1}^2\right)\left(1-{\mathrm{R}}_{\mathrm{y}2·1}^2\right)\left(1-{\mathrm{R}}_{\mathrm{y}3·12}^2\right)}$$

(15)

where $\mathrm{n}$ is the number of samples, $\mathrm{i}$ is the time, ${\mathrm{R}}_{\mathrm{xy}}$ is the correlation coefficient between NDVI and individual climatic factors; ${\mathrm{R}}_{\mathrm{y}1·23}$ is the PCC between NDVI (y) and climatic factor 1 with fixing factors 2 and 3; ${\mathrm{R}}_{\mathrm{y}·123}$ is the CCC of NDVI (y) with climatic factors 1, 2 and 3. The T-test and F-test were used to determine the significance levels of PCC and CCC, respectively, according to whether they passed $\mathsf{\alpha }=0.05,$ to judge the significance of the correlation coefficient.

Due to the spatial heterogeneity of the correlation between NDVI and climatic factors of the study area, the driving force rules for vegetation change were regionalized and summarized according to the significance tests of PCC and CCC (T-test and F-test, respectively) [63,64,65]. To ensure the maximum consistency, regional continuity and spatial non-repeatability of each category, the pixels that passed the F-test at $\mathsf{\alpha }=0.05$ were extracted for further climate-driven regionalization, while the pixels that failed the F-test at $\mathsf{\alpha }=0.05$ were regarded as non-significant climate-driven. Moreover, the climate-driven factors were divided into eight categories according to the PCC values and T-test performances between NDVI and each climatic parameter (

Table 3. Regionalization rules of the driving force for vegetation change.

Driving Factors of Vegetation Change	Rules
	F	TP	TT	TR	PCC
Precipitation, positively driven	F > F0.05	t > t0.05			PCCP >0
Precipitation, negatively driven	F > F0.05	t > t0.05			PCCP ≤ 0
Temperature, positively driven	F > F0.05		t > t0.05		PCCT > 0
Temperature, negatively driven	F > F0.05		t > t0.05		PCCT ≤ 0
Radiation, positively driven	F > F0.05			t > t0.05	PCCR >0
Radiation, negatively driven	F > F0.05			t > t0.05	PCCR ≤ 0
Hydrothermal factors, combined	F > F0.05	t > t0.05	t > t0.05	t > t0.05
	F > F0.05	t > t0.05	t > t0.05
	F > F0.05	t > t0.05		t > t0.05
	F > F0.05		t > t0.05	t > t0.05
Weak climate-driven	F > F0.05	t <t0.05	t < 0.05	t < t0.05
Non-significant climate-driven	F ≤ F0.05

F: F-test significance for complex correlation analysis (p < 0.05); T_P: T-test significance for partial correlation analysis between NDVI and precipitation; T_T: T-test significance for partial correlation analysis between NDVI and temperature; T_R: T-test significance for partial correlation analysis between NDVI and radiation; PCC_P: partial correlation coefficient between NDVI and precipitation; PCC_T: partial correlation coefficient between NDVI and temperature, PCC_R: partial correlation coefficient between NDVI and radiation.

).

3. Results

3.1. Spatiotemporal Change of NDVI and Hydrothermal Factors

During 2001–2020, the growing season NDVI showed a fluctuating increase rate at 0.0009/year (p = 0.0033), and the overall trend was smooth with a peak value of 0.69 appearing in 2020 (Figure 5b). The growing season mean NDVI in the study area was 0.66, indicating a pattern of decreasing initially, followed by an increasing trend from east to west (Figure 5a). The high-NDVI areas were mainly located in the east of the Sichuan basin (Ya’an and Luding), where the elevation was relatively low (<2500 m) and the main vegetation types were BF, shrub and CV, with the distribution of cropland along the route of the railroad. In contrast, the low-NDVI areas were mainly concentrated in the southeast Tibet Mt. Hengduan area (Basu and Bomi), where the mean elevation was above 5000 m and the main vegetation type was AV. The NDVI spatial distribution was influenced by the variations in elevation, topography, climate conditions and vegetation types.

A significant “three-dimensional zonation” vegetation distribution pattern [43] was illustrated in the study area (Figure 5a), indicating that NDVI decreased gradually with increasing elevation. The low-elevation areas (<2500 m) were concentrated in the Sichuan basin (Ya’an and Luding), accounting for 3.71% of the total pixels, with the mean NDVI value of 0.86. The middle-low elevation areas (2500–3500 m) were mainly distributed in the river valley regions where human activities were abundant, accounting for 9.33% of the total pixels with the mean NDVI of 0.78. The middle-high elevation areas (3500–4500 m) were the main pasture zones of the study area, which were dominated by meadows and shrubs, accounting for 49.14% of the total pixels, with the mean NDVI of 0.75. The high elevation regions (4500–5500 m) were located in the western Sichuan plateau basin and the southeast Tibet Mt. Hengduan area (Litang–Bomi), which were dominated by meadows and AV, and which accounted for 37.01% of the total pixels, with the mean NDVI of 0.5. The extremely high elevation areas (>5500 m) only accounted for 0.8% of the pixels, with the mean NDVI of 0.11, and were largely covered by ice and snow, with a scattered distribution of AV. As the elevation of the entire study area increases gradually from east to west, the fertility, pH and depth of the soil decrease markedly, affecting the absorption of organic matter by plants and resulting in the gradual decrease of NDVI [14].

The main vegetation types in the study area were NF, BF, shrub, meadow, AV and CV (Figure 3), accounting for more than 90% of the total area. In the past 20 years, all kinds of vegetation have showed upward trends in the interannual growing season NDVI, with AV having the most rapid rate of increase (0.0014/year, p = 0.0076), followed by BF (0.0011/year, p = 0.0001), meadow (0.0010/year, p = 0.0115), shrub (0.0009/year, p = 0.0076) and CV (0.0008/year, p =0.0023), while NF had lowest growth rate (0.0006/year, p = 0.0147) (Figure 6a). The annual mean NDVI of BF was highest, followed successively by NF, CV, shrub, meadow and AV (Figure 6b). The interannual fluctuations of NF and BF were low, caused by the sensitivity reduction of NDVI in high-vegetation-density areas. In contrast, those of meadow and CV were high due to the influence of grazing and human activities.

During 2001–2020, the total precipitation and mean temperature indicated fluctuating upward trends at the rates of 1.1990 mm/year (p = 0.5596) and 0.0458 °C/year (p = 0.0012), respectively (Figure 7a,b). In contrast, the average radiation decreased markedly at the rate of 0.2680 MJm⁻²/year (p = 0.0007) (Figure 7c). In terms of spatial pattern, the precipitation of the low-elevation areas (<2500 m) increased most rapidly, and the rainfall slope decreased with the increasing elevation (Figure 7a). The spatial distribution of the temperature slope indicated that temperature in the middle-high elevation areas (3500–4500 m) and high-elevation regions (4500–5500 m) increased most markedly, while that in the low-elevation areas (<2500 m) decreased (Figure 7b). The radiation decreased over the past 20 years throughout the whole study area, which occurred more markedly in the area where the elevation was less than 4500 m (Figure 7c).

3.2. Characteristics of Vegetation Change

The Theil–Sen slope analysis and MK test were used to explore and effectively characterize the spatial distribution of the NDVI variation trend for the research area. According to the NDVI slope results, the pixels with an NDVI slope less than or equal to −0.0005 were classified as the degradation area, and the pixels with an NDVI slope greater than or equal to 0.0005 were classified as the improvement area, while the areas with an NDVI slope that ranged from −0.0005 to 0.0005 were classified as the stable area. In addition, based on the confidence level of 0.05, the significance test results were divided into significant ($\left|\mathrm{Z}\right|>1.96$) and insignificant changes ($\left|\mathrm{Z}\right|<1.96$). The NDVI trend was categorized into five classes by integrating the Theil–Sen slope analysis and MK test (Table 1). Figure 8a shows the spatial distribution of the NDVI trend; the two stacked bars indicate the percentages of trend classification based on elevation levels (Figure 8b) and vegetation types (Figure 8c), respectively.

In Figure 8a, the slight improvement areas account for the largest proportion (47.53%), mainly distributed in the western Sichuan plateau basin area (Batang–Changdu) and the west side of the southeast Tibet Mt. Hengduan area (Bomi), were mainly dominated by shrub and meadow. The stable vegetation areas accounted for more than one-fifth (22.74%) of the pixels distributed in the extremely high-elevation southeast Tibet Mt. Hengduan area (Basu), and they were mainly covered by AV. The slight degradation areas, which accounted for 20.99% of the pixels, were mainly distributed in the western Sichuan plateau basin area (Batang–Changdu) and southern Tibet valley area (Linzhi and Milin), where the main vegetation types were BF, shrub and meadow. The significant improvement areas (7.57%) were concentrated in the Sichuan basin (Ya’an and Luding), where hydrothermal conditions for the growth of BF and CV were sufficient. The significant degradation areas accounted for the least number of pixels (1.17%); these were distributed around the cities and dominated by CV, indicating the negative impact of human activities on vegetation growth. The percentages of trend classification based on elevation levels (Figure 8b) showed that the significant improvement areas’ proportion gradually decreased with increasing elevation, while the significant degradation areas were mainly distributed near cities in the middle and low elevation areas (<3500 m), where human activities were abundant. Figure 8c presents how a slight improvement accounted for the largest percentage for all vegetation types, except for CV, which had the largest proportion of significant improvement.

The coefficient of variation was calculated to investigate the volatility of NDVI variation, which was divided into five types: high volatility, relatively high volatility, moderate fluctuation, relatively low volatility and low volatility. Figure 9a shows the spatial pattern of NDVI volatility; the two stacked bars represent the percentages of volatility classification based on elevation levels (Figure 9b) and vegetation types (Figure 9c), respectively.

The volatility of NDVI variation in the study area was characterized by obvious high, moderate and low fluctuations, with low volatility dominating (Figure 9a). High-volatility areas (16.42%) were mainly distributed in the southeast Tibet Mt. Hengduan area (Chaya and Bomi) and southern Tibet valley area (Linzhi and Milin), where the elevations are extremely high. Relatively-high-volatility regions accounted for 7.24% of the pixels and were distributed as ring clusters around the boundary of high-volatility regions. Moderate fluctuation (14.97%) and relatively low volatility (40.13%) were widely distributed in the middle- and low-elevation areas of the western Sichuan plateau basin (Batang–Changdu). The low-volatility areas, which accounted for 21.24% of the pixels, were concentrated in the Sichuan basin (Ya’an and Luding), western Sichuan plateau basin (Batang–Changdu) and river valley areas where the elevation is low. In general, the volatility of NDVI variation was markedly controlled by the elevation gradient (Figure 9b), which increased with rising elevation in the areas where the elevation was less than 5000 m, and decreased significantly in the extremely high-elevation areas (>5000 m). In terms of different vegetation types, NF, BF and CV had low volatility, shrub and meadow had relatively low volatility, while AV had high volatility, which coincided with the vertical distribution characteristics of vegetation and elevation (Figure 9c).

Based on the Theil–Sen slope and Hurst exponent results, the future sustainability of vegetation was predicted and divided into five classes: continuous degradation, no continuous degradation, stable, no continuous improvement and continuous improvement (Table 2). Figure 10a shows the spatial distribution of NDVI sustainability; the two stacked bars indicate the percentages of sustainability classification based on elevation levels (Figure 10b) and vegetation types (Figure 10c), respectively.

Figure 10a illustrates that “no continuous improvement” areas accounted for the largest proportion (41%), with a wide distribution throughout the whole study area. The stable (19.5%) and “no continuous degradation” (16.6%) areas were mainly concentrated in the extremely high-elevation zones of the southeast Tibet Mt. Hengduan area (Chaya and Bomi) and western Sichuan plateau basin (Batang–Changdu). The “continuous improvement” areas were distributed in the low elevations, accounting for 16.37%, while the “continuous degradation” areas accounted for the least number of pixels (4.76%), with a sporadic distribution near the cities. In terms of the NDVI sustainability percentages of different elevations (Figure 10b) and vegetation types (Figure 10c), with increasing elevation, the “continuous degradation” proportions decreased initially, followed by increases, and the “continuous improvement” proportions declined gradually, while CV and BF might have a higher possibility of continuous degradation in the future.

3.3. Correlation Analysis between NDVI and Hydrothermal Factors

During 2001–2020, the PCC between NDVI and precipitation was in the range of −0.90 to 0.91, with the mean value of 0.03; positively correlated pixels (50.83%) were slightly more prevalent than negatively correlated pixels (44.06%), while “significant positive correlation” and “negative correlation” pixels accounted for only 3.36% and 1.75%, respectively (Figure 11a,

Table 4. Aerial proportions of partial correlation levels of NDVI and precipitation, temperature and radiation, respectively (%).

Partial Correlation Level	Significant Negative Correlation	Negative Correlation	Positive Correlation	Significant Positive Correlation
NDVI vs. Precipitation	1.75	44.06	50.83	3.36
NDVI vs. Temperature	1.93	38.12	58.12	1.84
NDVI vs. Radiation	1.95	52.76	43.49	1.80

). The PCC between NDVI and temperature ranged from −0.88 to 0.88, with the mean value of 0.03, and positive correlation accounted for the largest proportion (58.12%), followed by negative correlation (38.12%); the pixels that passed the 0.05 significance test were less than 4% of the total (significant positive correlation 1.84% and significant negative correlation 1.93%) (Figure 11b, Table 4). The PCC between NDVI and radiation was between −0.91 to 0.94, with the mean value of −0.03; there were 52.76% negatively correlated pixels and 43.49% positively correlated pixels, and pixels with a significant positive or negative correlation accounted for 1.8% and 1.95%, respectively (Figure 11c, Table 4). The spatial distribution patterns of the partial correlation level between NDVI and three hydrothermal factors also showed marked heterogeneity; the negatively correlated pixels were concentrated in the Sichuan basin (Ya’an and Luding) and western Sichuan plateau basin (Batang–Changdu), while positively correlated pixels were distributed in the southeast Tibet Mt. Hengduan area (Chaya and Bomi) and western Sichuan plateau basin (Batang–Changdu). Overall, the vegetation responses to climatic factors were not strong in the study area.

Figure 11d indicates that precipitation’s positive effect on the 2500–3500 elevation area was most marked, followed successively by the 3500–4500 area, >5500 area, <2500 area and 4500–5500 area, showing that the significant positive influences of precipitation on the vegetation were distributed in the medium-elevation regions. Temperature’s positive effects on different elevation areas ranked as follows: the >5500 area was most marked, followed successively by the <2500 area, 4500–5500 area, 3500–4500 area 2500–3500 area, indicating that vegetation of extremely high- and low-elevation areas was more sensitive to temperature. In contrast, the negative impact of radiation on vegetation decreased with a rising elevation. In terms of different vegetation types (Figure 11e), temperature’s positive impact on the AV was most significant, followed successively by CV, shrub, BF, meadow and NF, illustrating that temperature was the main contributor to the growth of AV, CV and shrub. The positive influence of precipitation on shrub was most marked, followed by NF, AV, BF, meadow and CV, illustrating that the growth of shrub, NF and AV was more sensitive to precipitation. Meanwhile, due to human interference, rainfall inhibited the growth of CV. The PCCs between different vegetations with radiation were all negative, with BF having the smallest value; that increased successively in the order of NF, shrub, meadow, AV and CV, indicating that radiation negatively impacted the vegetation growth, especially for BF and NF. Overall, the vegetation growth was more markedly correlated with temperature compared to precipitation and radiation in the study area.

The time-lag effect between NDVI and climatic factors has been demonstrated in previous studies, where a shorter lag time indicates a more immediate vegetation response to climate change [59,66,67]. We used a partial correlation analysis to investigate the hysteresis effect of NDVI and hydrothermal parameters at the 1-month interval, and we took the month corresponding to the maximum partial correlation coefficient as the lag time (Figure 12,

Table 5. Aerial proportions of lag months between NDVI and climatic factors (%).

Lag Months	0	1	2	3
Precipitation	67.45	24.28	5.99	2.28
Temperature	35.36	55.69	6.91	2.04
Radiation	24.25	34.63	36.99	4.13

).

The time-lag relationship between NDVI and precipitation was no-lag for 67.45% of the pixels, 1-month for 24.28% of the pixels, 2-month for 5.99% of the pixels and 3-month for 2.28% of the pixels, showing a relatively short response time (Figure 12a). Meanwhile, the lag time between NDVI and temperature was no-lag for 35.36% of the pixels, 1-month for 55.69% of the pixels, 2-month for 6.91% of the pixels and 3-month for 2.04% of the pixels, illustrating a moderate response time (Figure 12b). In contrast, the lag time between NDVI and radiation was no-lag for 24.25% of the pixels, 1-month for 34.63% of the pixels, 2-month for 36.99% of the pixels and 3-month for 4.13% of the pixels, indicating that the response time was relatively long (Figure 12c). The results represented that the lag relationships between vegetation and three hydrothermal parameters in the study area were different. In terms of spatial distribution, the areas with a long lag vegetation response to precipitation and temperature were mainly distributed in the western part of the study area and the watersheds of major rivers (Figure 12a,b), and the areas with a long lag vegetation response to radiation were concentrated in the watersheds, especially along the Yalong River and Dadu River (Figure 12c). There was no marked relationship between elevation and lag response time; the long lag vegetation response appeared more frequently in the extremely high-elevation areas (>5500 m). In addition, there was no significant variation between the lag-responses of different vegetation types, except for CV, which was mainly controlled by human activities. In general, the vegetation response to precipitation was the most immediate, followed by temperature and radiation.

3.4. Driving Factors for Vegetation Change

Based on the partial and complex correlations (Figure 11a–c and Figure 13) between NDVI and three hydrothermal factors over the past 20 years, we found there was spatial and temporal heterogeneity in the study area, and we regionalized the driving force rules for vegetation change.

The hydrothermal factor driving force had a scattered distribution pattern, where the non-significant climate-driven force accounted for the largest number of pixels (93.8%), followed by weak climate-driven force (1.78%) and combined hydrothermal factors (1.44%). In terms of the single factor, positively driven temperature accounted for 0.63% of the pixels, and positively driven precipitation, negatively driven temperature and positively driven radiation accounted for a similar number of pixels (0.58%, 0.51% and 0.53%, respectively), while negatively driven radiation and negatively driven precipitation accounted for 0.41% and 0.32% of the pixels, respectively. Since more than 90% of the pixels were non-significantly climate-driven, we excluded these pixels and only used climate-driven pixels for our further analysis (Figure 14a–c). According to Figure 14a, as the elevation rose to 5000 m, the positively driven temperature percentage increased gradually, while the positively driven precipitation percentage increased initially, followed by decreases, and the proportions of positively and negatively driven radiation declined smoothly. Similarly, the climate driving force for different vegetation types was also varied; from the perspective of single climatic driving factors, the growth of NF, shrub and meadow was mainly controlled by precipitation, AV was more markedly driven by temperature, and BF and CV were more significantly dominated by radiation (Figure 14b).

To explore the dynamic changes of the climatic driving force, we analyzed the classification of the climate-driven type from a 9-year duration (2001–2009) to a 20-year duration (2001–2020) (under-9-year durations were excluded for their large errors due to the insufficient sample). For 2001–2009 versus 2001–2015, the percentages of single hydrothermal factors showed an upward trend, and the combined hydrothermal factor trend slightly fluctuated, while the weak climate-driven factor indicated a downward trend. After 2001–2015, all trends flattened, with no marked variation (Figure 14c).

4. Discussion

During 2001–2020, the growing season NDVI of the whole study area showed an upward trend, with a significant spatial difference, indicating a pattern of overall improvement and partial degradation [68], which was consistent with previous studies on the vegetation dynamics of the QTP [69,70,71]. The “significant improvement” areas were mainly located in the eastern part of the study area (Figure 8), where the hydrothermal conditions were sufficient and the main vegetation types were CV, shrub and BF (Figure 3). In contrast, the distribution of “significant degradation” areas were more scattered, mainly around the cities that were more likely affected by intensive human activities. Due to the implementation of ecological programs, such as the Grazing Withdrawal Program (GWP), which aims at restoring and protecting the ecosystem structure and function by fencing, scarification and reseeding [72,73], and the Grain to Green Program (GTGP), which encourages the replacement of cropping and livestock grazing in fragile areas with trees and grass [74], the land-use and cover changes that positively contribute to the overall vegetation greening were determined. Additionally, the growth rate variations of different vegetation types were shown to possibly be caused by human activities.

Due to the special geographic characteristics of the study area, animal husbandry is the main local industry and the regional economic development is highly dependent on meadow resources. So, overgrazing is considered one of the main reasons for regional meadow degradation [42,75,76]. Apart from that, due to the lack of firewood, Yak dung is collected by residents for fuel [77]; the nutrients of Yak dung that should be provided to the meadow are removed, leading to regional meadow degradation to some extent [78]. Overall, our results illustrated that the NDVI gradually decreased with the increase of elevation, indicating the NDVI spatial distribution characteristics are high in the east and low in the west, which was consistent with the vertical zonality of elevation and vegetation types. It is worth noting that both the overall NDVI and NDVI values of AV, shrub and meadow in 2020 showed a convergence of significant increase. The possible reason is that warming leads to the melting of alpine snow and alleviates the water stress of high-elevation areas mainly covered by AV, shrub and meadow. In addition, the increase in annual precipitation strengthens this phenomenon, and offsets the evapotranspiration caused by warming to some extent.

The NDVI variation also showed marked spatial heterogeneity. Long-rooted plants, mainly distributed in the eastern part of the study area and watershed areas, are more resistant to extreme climate conditions and less affected by hydrothermal condition changes since the water and heat for vegetation growth are sufficient. As elevation increases, the water availability for vegetation growth decreases, along with the NDVI variation, while the ecological vulnerability increases consequently [79]. The mean value of the Hurst exponent was 0.45, indicating the uncertainty of future vegetation sustainability. The high randomness areas accounted for nearly 80% of the whole study area; the “continuous improvement” areas were distributed in the watershed regions where hydrothermal conditions were more appropriate for plant growth, while the “continuous degradation” areas were mainly concentrated in the densely populated areas. According to statistical results, 19–45 person/km² is considered as the optimal range for vegetation restoration; an intensive population enhances the possibility of vegetation destruction and places pressure on ecological protection, which negatively affects vegetation restoration and growth [80,81].

The vegetation dynamics of the study area showed a hysteresis response to hydrothermal conditions (lag time of 1–3 months), with precipitation and radiation representing the most and least immediate responses, respectively, a finding that was similar to those of the previous studies on the QTP [66]. Vegetation growth needs a lag time for accumulating soil moisture, heat and nutrients [82]. Plus, the process of converting precipitation into water available for vegetation growth is related to several factors including soil type, land cover, vegetation physiological characteristics and climatic conditions [83]. Due to differences in the soil type and topography, the water-holding capacity of soil varies; some kinds of soil can retain water for a long period after precipitation, which has relatively lasting effects on vegetation growth. Apart from that, deep-rooted plants, such as BF, may have a longer “memory” for precipitation [84,85], while shallow-rooted plants, such as meadow and AV, are widely distributed in the study area, but have a relatively short growth cycle, which may be the main reason for the short lag relationship between NDVI and precipitation. In contrast, the heat factor shows a longer hysteresis response than the water factor. The increasing temperature affects soil microbial activities and promotes the mineralization of soil organic nitrogen and organic phosphate, which influence the nutrient supplementation for vegetation growth [86]. In addition, warming affects the melting of glaciers and permafrost in high-elevation areas [87], resulting in an increase in runoff from watersheds, which has a lasting impact on vegetation growth [88,89]. Interestingly, all three hydrothermal factors had a marked vegetation hysteresis response in the watershed regions; sufficient runoff resources in watershed areas would weaken the vegetation response to climatic factors.

The world has experienced significant warming over the past few decades [90,91]; a phenomenon called elevation-dependent warming (EDW) [30,92,93] appears in the high-elevation areas, which indicates a higher warming rate in the mountain regions than in plain areas. However, the low temperature is still considered one of the main limiting factors for vegetation growth in high-elevation areas [94]. The study area has a relatively humid climatic condition with abundant water resources, resulting in the low sensitivity of vegetation to precipitation. Plus, the gradual decrease in radiation during the study period led to a reduction in vegetation photosynthesis [66], which also negatively affects vegetation growth. So, temperature is the most important and active hydrothermal driving factor for vegetation growth in the study area [38]. Since a warm temperature is a key driver for seed germination in temperate alpine habitats, in a phenomenon called “warm-cued germination” [95], early spring germination allows the first-year plants of shrubs and AV to grow large enough to withstand the cold winters on the QTP [96,97]. Most of the study area has high elevations where the temperature is relatively low; except for the Sichuan basin (Ya’an and Luding) in the east, warming may directly promote vegetation photosynthesis and gradually reduce the limitation of low temperature on vegetation growth [98,99], which results in the vegetation greening trend [14,95]. Interestingly, the growing season average temperature of the western Sichuan plateau basin area (Batang–Changdu) located in the central part of the study area is relative moderate (in the range of 13–16 °C), which has a certain inhibitory effect on the growth of meadows, which may be due to the fact that meadows have a weak ability to absorb deep soil moisture. As the precipitation decreases, warming could promote vegetation photosynthesis until reaching the optimal temperature for photosynthesis. After that, the continued increase in temperature would promote vegetation respiration, accelerate nutrient consumption and water evaporation and reduce the accumulation of organic matter content [100], leading to the inhibition of vegetation growth [101]. In addition, as the temperature increases in the western Sichuan plateau, permafrost environmental degradation occurs and the active layer thickness rises, which lead to a marked decrease in the soil organic matter content, increase in soil evaporation [72,102], reduction of water accessibility in the topsoil [103,104,105] and consequently, limitation of vegetation growth.

The results of climate driving factor regionalization and our study of their interannual variability indicated that the vegetation changes in the study area were influenced by both climatic and non-climatic factors. According to the previous studies on the QTP, the contribution of non-climatic drivers to grassland variations is nearly double that of climatic drivers [70]; non-climatic factors overwhelmingly dominate the impact on grass growth. The vegetation change over a long time series is a complex process that is vulnerable to non-climatic factors such as grazing, land-use change, urban expansion, project construction and plant disease in the study area. Artificial interference partially offsets the positive effect of warming on vegetation growth [14,106,107,108] and plays an important role in vegetation dynamics. Moreover, as shown in Figure 14c, weak climatic factors had a gradual downward trend, which may be because of the ecological program implementations since 2000.

Considering the impact of hydrothermal conditions on vegetation growth, we analyzed the changing pattern of NDVI and its response to hydrothermal factors over the past 20 years, and we concluded that climatic factors are not the most important driving force of vegetation growth along the Ya’an–Linzhi section of the STR. However, there were some limitations to this study. The interpolated results for hydrothermal factors may have had large errors in the west of the study area with sparse meteorological stations, more so than in eastern regions. Extra auxiliary data is necessary to improve the climatic parameter interpolation accuracy in further studies. In addition, the long time series of vegetation change is complex given that it is not only closely related to hydrothermal factors but also intricately correlated with various factors such as human activities, runoff, topography and evapotranspiration. Cao et al. [109] indicated that shifts in vegetation phenology induced by climate change are substantially modifying various ecosystem processes, and those changes can, in turn, affect weather and climate systems. In particular, the areas along the STR have complex topographic conditions and significant elevation differences; elevation and terrain could affect the regulation of the resources for vegetation growth such as temperature, humidity and light, which influence the vegetation change indirectly [110,111,112]. Some studies indicated that vegetation growth is markedly dominated by climatic factors at a large scale and also controlled by terrain at a relatively small scale [101]. Specifically, orientation leads to the differentiation of microenvironments by affecting the redistribution of solar radiation, evapotranspiration and precipitation [113,114]. The slope contributes significantly to variability in the soil water content by affecting infiltration, drainage, runoff and potential insolation [113]. Meanwhile, elevation determines the vertical zonality of soil and vegetation [114], and it affects soil moisture redistribution [115]. Therefore, other factors that may cause vegetation degradation such as project construction, terrain, runoff, evaporation and human activities should be comprehensively analyzed in the future. Moreover, the MOD13Q1 data provided every 16 days at a 250-m spatial resolution may be crude to use as the basis for regional studies. By applying cloud-removal approaches such as the MNSPI (modified-neighborhood similar-pixel interpolator) and ARRC (AutoRegression to Remove Clouds), high-precision Landsat time-series data could be used instead of MODIS data in further studies [116].

5. Conclusions

The growing season NDVI along the Ya’an–Linzhi section of the STR showed an increasing trend, with fluctuation at the rate of 0.0009/year (p < 0.05) during 2001–2020. The vegetation “significant improvement” areas were mainly distributed in the east of the study area, while the “significant degradation” areas were concentrated around cities. The regions with elevations less than 5000 m showed a greening trend with the increase of elevation. The vegetation response to climate change was temporally hysteretic and spatially heterogeneous, with a lag time of 1–3 months. Changes due to precipitation were the most immediate, followed by temperature and radiation. The vegetation growth was more positively correlated with hydrothermal factors in the west of the study area, while it was more negatively correlated with climatic conditions in the east of the study area. The hydrothermal factors’ combined driving force on vegetation accounted for 6.2% of the pixels during 2001–2020, and temperature-driven pixels accounted for the largest percentage in terms of single driving factors, indicating that vegetation change was influenced by climatic and non-climatic factors, with non-climatic factors playing a dominant role, and temperature contributed more than precipitation and radiation to the vegetation growth in the study area. Although the linear regression model showed that all types of vegetation had a positive greening trend, the future vegetation sustainability based on the Hurst exponent indicated a certain risk of degradation, meaning it deserves more active conservation to hedge the negative impact of railway construction and human activities on the vegetation improvement trend, and to reduce the risk of ecological degradation.