Variations over 20 Years in Vegetation Dynamics and Its Coupled Responses to Individual and Compound Meteorological Drivers in Sichuan Province, China

Abstract

This study presents an innovative investigation into the spatiotemporal dynamics of vegetation growth and its response to both individual and composite climatic factors. The Normalized Difference Vegetation Index (NDVI), derived from SPOT satellite remote sensing data, was employed as a proxy for vegetation growth. Multiple analytical methods, including the coefficient of variation, Mann–Kendall trend analysis, and Hurst index, were applied to characterize the spatiotemporal patterns of the NDVI in Sichuan Province from 2000 to 2020. The Standardized Precipitation Evapotranspiration Index (SPEI) was calculated using monthly precipitation and temperature data from 45 meteorological stations to examine the influence of composite climatic factors on vegetation growth, while the time lag effects between the NDVI and various climatic variables were also explored. Our findings unveil three key insights: (1) Vegetation coverage in Sichuan Province exhibited an overall increasing trend, with the highest NDVI values in summer and the lowest in winter. Significant NDVI fluctuations were observed in spring in the western Sichuan plateau and in winter in northern, eastern, and southern Sichuan. (2) A significant upward trend in the NDVI was detected across Sichuan Province, except for Chengdu Plain, where a downward trend prevailed outside the summer season. (3) On shorter time scales, the NDVI was positively correlated with precipitation, temperature, and the SPEI, with a one-month lag. The response of the NDVI to sunlight duration showed a two-month lag, with the weakest correlation and a five-month lag in western Sichuan. This research advances our understanding of the complex interactions between vegetation dynamics and climatic factors in Sichuan Province and provides valuable insights for predicting future vegetation growth trends.

1. Introduction

Vegetation, as a crucial link connecting the atmosphere, soil, and biotic elements, plays a vital role in soil and water conservation, mitigating the rise of greenhouse gas concentrations, balancing surface radiation, carbon, and nitrogen cycles, and maintaining ecosystem stability. The dynamic changes in vegetation also continuously influence the Earth’s environment and human well-being [1,2]. The relationship between vegetation and climate is complex and intertwined; climate change is a significant driving factor for vegetation change, while vegetation is essential feedback to the degree of climate change impact. Vegetation is highly sensitive to environmental changes, and over a long time series, multiple climatic factors jointly influence and constrain changes in vegetation cover [3,4,5]. Therefore, studying vegetation change trends and their inter-relationships with climatic factors not only provides a scientific basis and technical support for ecological and environmental protection, resource utilization, and human survival and development but is also key to understanding the mechanisms of dynamic vegetation changes. However, previous studies have mainly focused on the relationship between vegetation and a single meteorological factor (such as precipitation or temperature), with less consideration of the combined effects of multiple meteorological factors. In fact, vegetation growth is influenced by a combination of meteorological factors, and there may be complex interactions and lag effects between different factors. Thus, a comprehensive study of the lag effects of the Normalized Difference Vegetation Index (NDVI) with single-factor and composite meteorological factors is crucial for revealing the response mechanisms of vegetation to climate change and providing important references for predicting future vegetation dynamics under climate change.

The NDVI can accurately reflect the surface vegetation cover and is an essential indicator for characterizing vegetation growth and spatiotemporal changes [6,7,8,9]. In recent years, many scholars have used the NDVI to assess the characteristics of vegetation changes at different spatial and temporal scales and their response to climate change [10,11,12,13,14,15,16,17,18,19]. However, most of these studies have focused on the relationship between vegetation and a single climatic factor, lacking a comprehensive analysis of the impact of integrated hydrothermal climatic factors on vegetation changes in Sichuan Province.

Sichuan Province is not only an ecological barrier in the upper and middle reaches of the Yangtze River but also an ecologically fragile and climate-sensitive region [20,21]. Studying the spatiotemporal changes of the NDVI and its correlation with single hydrothermal factors and comprehensive factors is conducive to revealing the adaptability of vegetation in different regions to changes in climatic conditions. An in-depth analysis of the impact of climatic factors on vegetation growth can clarify the dominant factors affecting vegetation growth, which is not only helpful for ecological and environmental protection and climate change response in Sichuan but also has important significance for enhancing the ecosystem service functions and promoting sustainable development in Sichuan Province [22,23].

The innovation and highlights of this study lie in that it comprehensively considers the lag effects of single-factor and composite meteorological factors on the NDVI, providing a more complete understanding of the complex relationship between climate and vegetation. Hence, this study aims to analyze the spatiotemporal variation characteristics of vegetation cover in Sichuan Province using methods such as Theil–Sen median trend analysis, the Mann–Kendall statistical test, and Hurst index analysis. Furthermore, by combining single hydrothermal factors and comprehensive hydrothermal climatic factors, the study explores the influencing factors and lag effects of vegetation cover changes in Sichuan Province using correlation analysis. The findings of this study are expected to provide scientific support for ecological and environmental protection, climate change adaptation, and sustainable development in Sichuan Province and contribute to the understanding of the complex relationship between climate and vegetation.

2. Materials and Methods

2.1. Study Area

Sichuan Province in southwestern China (97°21′ E to 108°12′ E, 26°03′ N to 34°19′ N) spans 486,000 square kilometers (Figure 1). It includes five major geomorphic units: the Sichuan Basin, Qinghai–Tibet Plateau, Hengduan Mountains, Yunnan–Guizhou Plateau, and Qinling–Daba Mountains. The terrain is higher in the west and lower in the east, with an average elevation of approximately 2500 m and a maximum elevation difference exceeding 6 km [24,25,26]. The significant topographical variations and differing monsoon circulation patterns contribute to a rich and diverse climate in Sichuan Province, exhibiting marked spatial variation [27]. The western plateau and mountainous areas are characterized by a cold temperate climate, with cold winters and cool summers, abundant sunshine, an average annual temperature ranging from 4 to 12 °C, and annual precipitation between 500 and 900 mm; the eastern basin features a subtropical monsoon climate, with warm winters and hot summers, less sunshine, an average annual temperature of 16–18 °C, and annual precipitation ranging from 1000 to 1300 mm [28]. The region boasts a variety of vegetation types, primarily subtropical shrublands, evergreen broad-leaved forests, and alpine meadows on the plateau mountains. This paper divides Sichuan Province into five regions—West Sichuan (WS), North Sichuan (NS), East Sichuan (ES), Central Sichuan (CS), and South Sichuan (SS)—based on administrative divisions for the study.

2.2. Research Datasets

The vegetation types in Sichuan Province are derived from the “1:1,000,000 China Vegetation Atlas” provided by the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (https://www.resdc.cn/ accessed on 1 November 2023), with an elevation dataset also sourced from the same center, featuring a spatial resolution of 1 km. The monthly spatial distribution dataset of the Normalized Difference Vegetation Index (NDVI) (2000–2020) originates from SPOT satellite remote sensing data supplied by the Resource and Environment Science and Data Center of the Chinese Academy of Sciences. The data format is ArcGIS GRID, featuring a spatial resolution of 1 km. The NDVI datasets are calculated by taking the maximum value over a 10-day period. The Maximum Value Composite (MVC) [29] selects the highest pixel value within a specified time window to represent that period’s composite NDVI. Higher NDVI values are more likely to indicate actual vegetation status, whereas lower values can be affected by atmospheric interference from clouds or aerosols. Therefore, the NDVI data derived through the MVC more accurately reflect vegetation conditions.

The climate data are sourced from the National Meteorological Science Data Center (https://data.tpdc.ac.cn/zh-hans/ accessed on 5 November 2023), selecting daily precipitation, temperature, and sun-shine duration records from 45 meteorological stations, including Shiqu, Dege, and Ganzi, ranging from 2000 to 2020. The Collaborative Kriging interpolation method within ArcGIS [30,31,32] was utilized to interpolate the meteorological data while incorporating station elevation data to enhance interpolation accuracy, thereby producing a meteorological dataset that aligns spatially with the resolution of the NDVI raster data.

2.3. Methods

2.3.1. Coefficient of Variation Method

The coefficient of variation [33] was used to represent the fluctuation characteristics of the NDVI over time, indicating the degree of variation or dispersion among values. This metric is often applied to analyze vegetation dynamics across various temporal scales. While the CV is effective for data with differing units and means, it has limitations when the mean is close to zero or the data include extreme, negative, periodic, or trend-based values. The calculation process is shown in Equation (1).

$$\mathrm{C}\mathrm{V}={\displaystyle \frac{\sqrt{\frac{1}{n-1}\sum _{i=1}^n{({\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_i-\overline{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}})}^2}}{\overline{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}}}$$

(1)

where $\mathrm{C}\mathrm{V}$ is the coefficient of variation of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$; ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_i$ is the mean $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$ of the i-th year; $\overline{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}$ is the annual average value of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$. The value of the $\mathrm{C}\mathrm{V}$ is positively correlated with the volatility of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$ pixel values, and its magnitude represents the extent of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$ fluctuations.

2.3.2. Theil–Sen Median Trend Analysis and Mann–Kendall Statistical Test

Theil–Sen median trend analysis [34,35,36] is a nonparametric statistical method that offers high computational efficiency, does not require the test sample to follow a specific distribution, and is unaffected by a few outliers. It has significant advantages in small sample data statistics, and the calculation process is shown in Equation (2).

$$\beta =median\left({\displaystyle \frac{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_j-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_i}{j-i}}\right),2000\le i<j\le 2020$$

(2)

where $\beta$ is the trend of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$ changes, $\beta >0$ indicates an upward trend of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$, and $\beta <0$ indicates a downward trend of the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$.

Theil–Sen median trend analysis is often combined with the Mann–Kendall significance test [37]. The Mann–Kendall test, developed based on the principles proposed by Mann and Kendall, is a nonparametric rank test method suitable for testing the significance of trends in long time series data [38,39,40]. When using the Mann–Kendall method for NDVI trend testing, the time series NDVI values were treated as a set of independently distributed sample data, with the parameter Z serving as the indicator of NDVI decline at the pixel level.

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

(3)

where

$$S={\sum }_{j=i+1}^n{\sum }_{i=1}^{n-1}sgn\left(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_j-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_i\right)$$

(4)

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

(5)

$$sgn\left({\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_j-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_i\right)=\left\{\begin{array}{c}{1,\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_j>{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_i\\ 0,{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_j={\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_i\\ -1,{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_j<{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_i\end{array}\right.,j=i+1,\dots ,n$$

(6)

where $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_j$ and $\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_i$ are sample time series datasets, $n$ is the length of the time series dataset, and $sgn$ is the sign function. For the given confidence level $\alpha$, when $\left|Z\right|>Z_{1-\alpha /2}$, it is considered that the $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}$ value trend is significant at the $\alpha$ level. Otherwise, the trend is considered not significant. In this paper, a significant change is defined as $\alpha =0.1$.

2.3.3. Hurst Exponent Analysis

The Hurst exponent [41] is used to quantitatively describe the persistence of long-term time series data and is widely applied in fields such as climatology and hydrology. Accordingly, this study employs the Hurst exponent to predict the persistence characteristics of the NDVI in Sichuan across various time scales. The calculation process of the Hurst exponent is shown in Equations (7)–(11).

For the time series $\left\{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(\tau \right)\right\}$, where $\tau$ = 1, 2, …, n, the mean sequence is defined as follows:

$$\overline{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(\tau \right)}={\displaystyle \frac{1}{\tau }}{\sum }_{t=1}^{\tau }\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(\tau \right),\tau =\mathrm{1,2},\dots ,n$$

(7)

Cumulative deviation ($X\left(t,\tau \right)$) is defined as follows:

$$X\left(t,\tau \right)={\sum }_{t=1}^t\left(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(t\right)-\overline{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(\tau \right)}\right),1\le t\le \tau$$

(8)

The range sequence ($R\left(\tau \right)$) is defined as follows:

$$R\left(\tau \right)=maxX\left(t,\tau \right)-minX\left(t,\tau \right),1\le t\le \tau$$

(9)

Standard deviation ($S\left(\tau \right)$) is defined as follows:

$$S\left(\tau \right)={\left[{\displaystyle \frac{1}{\tau }}{\sum }_{t=1}^{\tau }{\left[\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(t\right)-\overline{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(\tau \right)}\right]}^2\right]}^{½},\tau =\mathrm{1,2},\dots ,n$$

(10)

The Hurst exponent (H) is defined as follows:

$${\displaystyle \frac{R\left(\tau \right)}{S\left(\tau \right)}}={\left({\displaystyle \frac{\pi \tau }{2}}\right)}^H$$

(11)

where H is the Hurst exponent, with a range of [0, 1]. When $H\in [0, 0.5]$, this indicates that the NDVI is inconsistent with past trends. When $H\in [0.5, 1]$, this indicates that the NDVI is consistent with past trends. When $H=0.5$, this indicates that the development trend of the NDVI is random.

2.3.4. Standardized Precipitation Evapotranspiration Index (SPEI)

The Standardized Precipitation Evapotranspiration Index (SPEI) quantifies the comprehensive impact of various meteorological factors on the dry and wet conditions of a study area by considering multiple elements. This approach makes the assessment of dryness and wetness more closely aligned with the actual conditions of the study area. In this paper, the Penman–Monteith formula [42] was used to calculate the SPEI. By normalizing the precipitation and potential evapotranspiration data, it is possible to compute the multi-temporal scale SPEI indices for each meteorological station. Since the NDVI data in this paper were calculated on a monthly scale, SPEI data on a one-month time scale were chosen to calculate the temporal lag characteristics of the NDVI in response to comprehensive meteorological factors. The calculation method can be found in references [43,44,45].

2.3.5. Correlation Analysis

The growth of vegetation depends on various climatic factors, which also interact with each other. To avoid this interaction, it is necessary to control one factor while analyzing the correlation between other factors and vegetation. This type of analysis is known as partial correlation analysis [46,47]. When calculating partial correlations, the correlation coefficient must first be determined using the following formula:

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

(12)

where n represents the year or quarter, while x and y denote the NDVI and temperature, precipitation, and sunshine hours, respectively. $x_i$ is the NDVI for the ith year or quarter, and $y_i$ is the average temperature, precipitation, and sunshine hours for the ith year or quarter. $r_{xy}=0$ indicates no correlation between $x$ and $y$; $r_{xy}>0$ suggests a positive correlation between $x$ and $y$; and $r_{xy}<0$ implies a negative correlation between $x$ and $y$ [48].

The calculation process of the partial correlation coefficient is shown in Equation (13):

$$r_{x,yz}={\displaystyle \frac{r_{xz-r_{xy}{r}_{yz}}}{\sqrt{(1-r_{xy}^2)(1-r_{yz}^2)}}}$$

(13)

where $r_{x,yz}$ is the partial correlation coefficient of $y$ and $z$ after controlling $x$.

2.3.6. Lag Time

The time lag effect of climate factors on vegetation was characterized by calculating the Pearson correlation coefficient. The maximum correlation coefficient, $r_i$, is considered the optimal correlation, and the corresponding month (i) is regarded as the optimal time delay.

$$R_i=corr(NDVI,Y_i) 1\le i\le 12$$

(14)

$$r_{max-lag}=max\left(r_i\right)$$

(15)

where $r_i$ represents the Pearson correlation coefficient at a lag of $i$ months, where i ranges from 1 to 12 (1 indicates no lag effect, while 2 to 11 indicates lags of 1 to 11 months). The NDVI denotes the monthly scale data for Sichuan Province from 2000 to 2020, and Y represents the monthly scale meteorological data (precipitation, temperature, sunshine hours, and the SPEI) from 1999 to 2019. $r_{max-lag}$ is the maximum value of $r_i$.

3. Results and Discussion

3.1. Changes in Single and Multiple Hydrothermal Factors in Sichuan Province

3.1.1. Single Hydrothermal Factors

Intra- and Inter-Annual Variations in Precipitation

The annual precipitation trends across different regions of Sichuan Province from 2000 to 2020 are shown in Figure 2. The average precipitation in Sichuan during this period was 798.58 mm. Precipitation showed a decreasing trend from 2000 to 2011, followed by an increase from 2012 to 2020. Between 2008 and 2011, all regions experienced varying degrees of decline in precipitation. SS had the highest multi-year average precipitation, reaching 955.49 mm. The highest annual precipitation in this region occurred in 2016 (1237.64 mm) and the lowest in 2011 (757.8 mm), with an overall increasing trend over the years. NS had the lowest average precipitation of 704.01 mm in the last 21 years, with the highest precipitation in 2018 (930.4 mm) and the lowest in 2002 (535.98 mm). The lower precipitation in NS and higher precipitation in SS are mainly due to the uplift effect of the terrain, the path of monsoonal airflows, and the combined impact of atmospheric circulation. The summer monsoon blows from a southeast direction, carrying a large amount of moisture. When the moist monsoonal airflow enters the Sichuan Basin, it first encounters elevated terrain in the south, which tends to produce precipitation, leading to the uneven distribution of precipitation within Sichuan Province.

The average seasonal precipitation in Sichuan from 2000 to 2020 is presented in

Table 1. Spatial variations in precipitation in Sichuan Province.

	WS/mm	NS/mm	ES/mm	CS/mm	SS/mm	Total/mm
Spring	148.00	160.25	251.95	177.02	240.23	172.69
Summer	415.09	343.94	334.26	349.79	407.67	394.39
Autumn	182.88	177.99	249.88	195.63	234.99	197.19
Winter	18.56	25.85	52.69	36.43	73.39	30.47
Growing season	656.85	579.42	639.81	594.78	714.37	650.23

Spring: March to May, Summer: June to August, Autumn: September to November, Winter: December to February of the following year, Growing Season: April to September.

. The growing season had the highest precipitation, reaching 650.23 mm, while winter recorded the lowest at 30.47 mm. In spring, ES had the highest rainfall, followed by SS, with WS having the lowest at 148.00 mm. In summer, rainfall was relatively high in both WS (415.09 mm) and SS (407.67 mm). In winter, the lowest precipitation was recorded in WS (18.56 mm) and the highest in SS (73.39 mm). Across regions, precipitation during the growing season is ranked as follows: SS > WS > ES > CS > NS.

Intra- and Inter-Annual Variations in Temperature

Figure 3 depicts the inter-annual temperature variations in Sichuan Province from 2000 to 2020. Over the past 21 years, the province’s average temperature was 12.64 °C, with a steady upward trend. ES recorded the highest multi-year average temperature, peaking at 18.06 °C in 2006 and reaching its lowest in 2012 (16.82 °C). The CS and NS followed, with average temperatures of 16.73 °C and 16.44 °C, respectively. WS had the lowest average (10.51 °C) but also showed a steady warming trend. The relatively warmer climate in ES, CS, and NS can be attributed to their lower latitude and subtropical location, combined with higher urbanization, dense populations, and intensive human activity, which elevate local temperatures. In contrast, WS, on the eastern edge of the Tibetan Plateau, has higher altitudes, which reduce temperatures. Additionally, the prevalence of snow, reflecting solar radiation, further lowers surface heat absorption, contributing to the region’s lower annual average temperature.

The average temperatures for each season in Sichuan Province over the past 21 years are shown in

Table 2. Spatial variations in temperature in Sichuan Province.

	WS/°C	NS/°C	ES/°C	CS/°C	SS/°C	Total/°C
Spring	11.76	18.04	18.42	18.28	16.58	13.89
Summer	18.64	26.41	27.90	26.79	24.44	21.39
Autumn	11.29	17.42	18.58	17.98	16.47	13.53
Winter	2.64	7.16	7.93	7.53	6.80	4.50
Growing season	16.58	23.72	24.86	24.11	21.99	19.09

. The average temperatures in each season in NS, ES, and CS were generally higher than those in WS and SS, and WS had the lowest temperatures in all seasons, which was due to the fact that most of the geographic location of WS is at high altitude, where the air is thin, and heat is lost quickly, resulting in lower temperatures. During the growing season, regional temperatures are ranked as follows: ES (24.86 °C) > CS (24.11 °C) > NS (23.72 °C) > SS (21.99 °C) > WS (16.58 °C).

Intra- and Inter-Annual Variations in Sunshine Duration

Figure 4 illustrates the inter-annual variation in sunshine hours across Sichuan Province from 2000 to 2020. Over this period, the province’s average sunshine duration was 1658.72 h. Between 2003 and 2005, sunshine hours in all regions showed a declining trend, followed by a gradual increase, peaking in 2013 at 1831.62 h. WS exhibited considerable variation, with an average sunshine duration of 1902.35 h. The highest annual duration in this region occurred in 2020 (2044.74 h) and the lowest in 2005 (1786.61 h). The region’s high elevation, low temperatures, limited precipitation, and minimal monsoon influence contribute to its relatively high annual sunshine duration.

As shown in

Table 3. Spatial variations in sunshine duration in Sichuan Province.

	WS/h	NS/h	ES/h	CS/h	SS/h	Total/h
Spring	533.86	396.04	380.22	346.05	348.65	476.04
Summer	452.09	422.96	522.83	406.19	402.53	452.26
Autumn	432.82	232.32	244.70	192.67	208.62	359.51
Winter	485.29	221.12	135.06	167.22	186.43	373.96
Growing season	935.90	794.24	903.73	732.56	729.88	891.15

, Sichuan’s seasonal average sunshine duration from 2000 to 2020 indicates that summer sunshine hours (452.26 h) were higher than in autumn (359.51 h) and winter (373.96 h) but slightly lower than in spring (476.04 h). The higher sunshine duration in spring compared to summer is mainly due to the clearer weather conditions and less cloud cover in spring, while summer is affected by the monsoon, resulting in more rainy and cloudy weather. In addition to having a longer sunshine duration in summer than in western Sichuan, the eastern region has the highest sunshine duration in spring, autumn, and winter. This is because western Sichuan has less precipitation and relatively fewer clouds, especially in winter and spring when sunny days are more frequent. Moreover, the higher elevation of western Sichuan allows for more direct sunlight to reach the surface, resulting in a longer sunshine duration in this region.

3.1.2. Multiple Hydrothermal Factors

Temporal Trends

When analyzing the influencing factors of vegetation, we should not only consider single factors, such as precipitation, air temperature, and sunshine hours, but also comprehensively consider the influence of multiple meteorological factors on vegetation growth. The multi-timescale SPEI value can clearly reflect regional dry and wet evolution and the status of available water resources. Figure 5 shows the multi-scale SPEI for each sub-region of Sichuan Province from 2000 to 2020. It can be seen that when the cumulative time scale is small, the SPEI fluctuates sharply, with obvious alternation between wet and dry conditions. As the cumulative time scale increases, the SPEI fluctuations weaken, showing consecutive dry and wet periods. The annual-scale SPEI can clearly indicate the switch between dry and wet periods. Sichuan Province was generally wet from 2000 to 2020, but there were dry periods in 2006–2007 and 2011–2012, with higher degrees of dryness. The wet–dry alternation in the WS region is pronounced, and it is wetter compared to other regions. The SS region experienced a relatively arid period from 2006 to 2012, with prolonged droughts. After 2015, the climatic conditions in this area gradually became wetter, with rarely occurring significant droughts. The CS region saw notable droughts in 2003, 2008, and 2020, with moderate to severe drought conditions in 2020. The ES and NS experience more frequent wet–dry alternations, but the droughts are relatively mild and have shorter durations.

Spatial Trends

After interpolating the intra- and inter-annual SPEI data for Sichuan Province, as shown in Figure 6, the spatial and temporal differences in dry and wet conditions during a year are obvious. Meteorological droughts occur mainly in the spring and winter seasons, followed by autumn. In spring, drought events primarily occur in the western part of WS, with the SPEI values varying between −0.504 and 0.393. The central and northern parts of WS, as well as the ES and SS, are relatively humid. During the summer, all regions in Sichuan Province are generally humid, with the SPEI values ranging from 0.473 to 1.378. The WS is more susceptible to the influences of monsoon and topographic precipitation during this season, resulting in relatively higher rainfall amounts and thus a higher degree of humidity compared to other regions. In winter, the drought area expands and intensifies, with the most severe drought occurring primarily in ES and the southeastern and northeastern parts of WS. The winter SPEI values range from −1.189 to −0.788. Overall, droughts mainly occur in spring and winter in the WS, which has the widest drought area and the greatest drought intensity. The main contributing factor to this phenomenon is low precipitation due to the high altitude of WS during this time period. The alternation of dryness and wetness in the ES region has an obvious seasonal characteristic, primarily due to the higher temperature levels in all time scales combined with seasonal precipitation deficits.

3.2. Spatiotemporal Variations in NDVI in Sichuan Province

3.2.1. Temporal Changes

Figure 7 illustrates the temporal dynamics of the NDVI in Sichuan Province. Overall, the province experienced a steady upward trend in the NDVI on both the intra-annual and inter-annual scales. However, from 2006 to 2013, the NDVI growth slowed or declined in some periods, particularly during the growing season and on the inter-annual scale. Over the past 21 years, the summer NDVI fluctuated between 0.684 and 0.776, with a multi-year average of 0.728. The peak value occurred in 2019, while the lowest was in 2000, showing a gradual upward trend of 0.0045 per year. The winter NDVI ranged from 0.372 to 0.456, with an average of 0.415. The highest winter value was recorded in 2016, and the lowest in 2000. The NDVI trends during the growing season and on an inter-annual scale were similar, with multi-year averages of 0.737 and 0.739, respectively. The multi-year mean values of the NDVI for all seasons in Sichuan Province, from the largest to the smallest, were as follows: summer (0.728) > autumn (0.659) > spring (0.522) > winter (0.415). In winter, due to low temperatures, limited precipitation, and short sunshine hours, the NDVI growth rate was the smallest, and vegetation cover was the lowest throughout the year. In spring, the temperature warmed up, the precipitation increased, and vegetation cover gradually increased. The NDVI growth rate in the growing season was the largest, and its multi-year average value was the highest for the year.

3.2.2. Spatial Changes

Mean Spatial Variations

The spatial distribution of the mean NDVI values at different time scales in Sichuan Province from 2000 to 2020 is depicted in Figure 8. In spring, vegetation coverage is low in the northwest of WS, gradually increasing from west to east, with higher coverage in the Sichuan Basin. In summer, significant improvements in vegetation coverage are observed throughout the province, except in the western and central regions of WS. This increase is attributed to summer being the main growing season for vegetation in Sichuan. In autumn, vegetation coverage declines compared to summer, particularly in the northwest and eastern regions of WS. Notable decreases are also observed in the southern areas of ES and NS, as well as the northern parts of CS and SS. This decline is primarily due to decreasing temperatures and precipitation, which slow down plant growth. In winter, vegetation coverage significantly decreases across the province, especially in WS and NS. The heavy snowfall in WS during winter covers the vegetation and leads to leaf drop, resulting in low coverage. Throughout the growing season, the northwest and central regions of WS and Chengdu Plain show lower vegetation coverage. Yearly observations indicate that vegetation coverage aligns with the growing season. Comparatively, WS and Chengdu Plain have lower vegetation coverage than other regions. WS, located on the eastern edge of the Tibetan Plateau, is characterized by high altitudes and mountainous terrain, with significant summer rainfall leading to natural disasters. The steep slopes and loose soil contribute to low vegetation coverage. In contrast, urbanization in the densely populated Chengdu Plain results in continuous reductions in vegetation coverage.

Spatial Fluctuations

The intra-annual and inter-annual NDVI fluctuation characteristics in Sichuan Province from 2000 to 2020 are shown in Figure 9. The NDVI in Sichuan Province mostly fluctuates at a low level, especially in the NS, ES, and SS, while the central and western regions of WS exhibit significant fluctuations. This is mainly due to the economic backwardness of WS, where economic development is dominated by agriculture, and the degradation of vegetation cover is caused by unreasonable farming methods. In spring, the northwest and central regions of WS have higher fluctuations; in summer, fluctuations in the whole province decline, and the vegetation cover is relatively stable, but some areas in the central part of WS still have higher-level fluctuations; in autumn, the fluctuation levels in the eastern part of WS, the southwestern part of NS, and the northern part of SS shift from low to lower; in winter, the fluctuation level of the whole province increases comprehensively, mainly concentrated in most areas of WS, NS, ES, and SS, with lower stability of vegetation cover; and the fluctuation changes in the growing season are consistent with those of the whole year, with a small number of areas in the western and central parts of WS still having high-level fluctuations and other areas having low-level fluctuations. The fluctuation in summer is the lowest compared with other seasons because the vegetation grows vigorously in summer, and the vegetation cover is high. The highest fluctuation in winter is due to the low temperatures and precipitation in winter, the non-vegetation growing season, withering leaves, and low vegetation cover.

Spatial Distribution

The spatial slope distribution of the intra-annual and inter-annual NDVI in Sichuan Province from 2000 to 2020 calculated by the Theil–Sen median method is shown in Figure 10. The NDVI in WS and Chengdu Plain showed a declining or unchanged trend in all periods; in summer, the growing season, and inter-annual vegetation cover, except for some areas around CS and a few areas in WS that showed unchanged or even decreasing trends, the NDVI in other areas showed an upward trend. The WS faces economic underdevelopment and adheres to traditional agricultural practices, leading to significant vegetation loss. This issue is compounded by the region’s naturally low vegetation cover, which is further exacerbated by frequent natural disasters. Chengdu Plain has experienced economic development and urban expansion, resulting in deforestation to make way for increased urban land use. Consequently, this has led to a noticeable degradation of vegetation cover in Chengdu Plain.

Based on the trend significance test criteria, the slope estimation results were combined with the Mann–Kendall test. The spatial distribution of the NDVI trend changes over different time scales was obtained according to the trend classification standards (Figure 11). From 2000 to 2020, Sichuan Province exhibited a continuous increase in vegetation cover area. The region with no change in the NDVI accounted for the largest proportion at 55.2%; areas with a highly significant increase in vegetation cover accounted for 27.3%, and areas with a decrease in vegetation cover made for 3.4%. In spring, regions with a significant decline in vegetation cover were primarily located in Chengdu Plain, while areas with significant increases in vegetation cover were mainly found in ES and the northern part of NS. During winter, vegetation cover in the northern region of WS and Chengdu Plain showed a significant downward trend, with the majority of Sichuan Province experiencing insignificant or slightly significant increases in vegetation cover. Inter-annual changes in vegetation cover aligned with those during the growing season. The areas where vegetation coverage levels have decreased are still Chengdu Plain and the central region of WS, the northern region of SS, and the southern region of ES and NS. The areas where vegetation coverage has increased are mainly the northern and southern regions of WS, the southern region of SS, and the northern region of ES. In contrast, areas with increasing vegetation cover were mainly located in the northern part of WS, the southern part of SS, and the southern parts of ES and NS. Overall, areas with reduced vegetation cover were primarily situated in the basin’s plain regions, which are urban agglomerations in Sichuan Province. The economic development of these towns has led to a continuous decline in vegetation cover. Conversely, areas with increased vegetation cover were mainly situated around the periphery of the basin—in WS, the southern part of SS, ES, and the northern part of NS—where the economy is relatively underdeveloped and vegetation destruction is minimal.

Spatial Future Change Trends

Using the Hurst exponent combined with Sen’s trend analysis, a persistence analysis was conducted on the inter-annual and intra-annual NDVI data of Sichuan Province (Figure 12). In spring, regions with continuous improvement in vegetation cover were distributed in central WS, ES, northern parts of NS, northern parts of CS, and southern parts of SS. Regions transitioning from an increase to a decrease were mainly located in the northern and southern regions of WS, southern parts of ES and NS, southern of CS, and northern parts of SS. Areas with continuous degradation were primarily found in Chengdu Plain. In summer, the area of regions with continuous improvement increased, mainly expanding in the eastern parts of WS. Regions transitioning from an increase to a decrease decreased and were mainly located in ES, southern parts of NS, and SS. Chengdu Plain continued to show a persistent decline. In autumn, the areas with continuous improvement in vegetation decreased, and the regions transitioning from an increase to a decrease increased. The areas that changed from a decrease to an increase are mainly distributed in the western part of WS. In winter, the areas with continuous improvement decreased significantly, with major reductions in WS and southern parts of SS. Regions transitioning from an increase to a decrease were distributed in SS, CS, northern and southern parts of NS, and southern parts of ES. Areas with continuous degradation were located in southern WS and Chengdu Plain. Regions transitioning from a decrease to an increase increased significantly, mainly in WS. The future trend changes during the growing season and inter-annually were almost aligned with those in summer.

3.3. Correlation Analysis and Lag Effect of NDVI on Climatic Factors in Sichuan Province

3.3.1. Correlation of NDVI with Single Climatic Factors

The spatial distribution of the partial correlation coefficients between the NDVI and temperature, precipitation, and sunshine hours in Sichuan Province from 2000 to 2020 that passed the significance test is shown in Figure 13. In spring, the NDVI is primarily positively correlated with temperature, accounting for 91.87% of the area, mainly distributed in WS and SS. The NDVI is mainly positively correlated with precipitation, covering an area of 93.6%, mainly located in WS, NS, and CS. The NDVI is primarily positively correlated with sunshine hours, accounting for 59.73% of the area, distributed in the northwestern part of WS and NS.

In summer, the NDVI was primarily positively correlated with temperature, accounting for 93.06% of the area, distributed in WS, SS, and the northern part of ES. The NDVI was mainly positively correlated with precipitation, with a positive correlation area ratio of 93.06% and a negative correlation area ratio of 24.07%, distributed in the southern part of WS and the northern part of SS. The NDVI was mostly positively correlated with sunshine duration, with a positive correlation area ratio of 78.37%, located in NS, CS, and the western part of WS and a negative correlation area ratio of 21.63%, mainly distributed in SS and the southeastern part of ES.

In autumn, the NDVI was primarily positively correlated with temperature, with the positive correlation area accounting for 69.9%, mainly distributed in the central part of WS; the negative correlation area accounted for 30.1%, mainly located in CS and SS. The negative correlation area between the NDVI and precipitation accounted for 43.86%, mainly distributed in NS and parts of WS. The NDVI was primarily negatively correlated with sunshine duration, with the negative correlation area accounting for 81.14%, distributed in NS, ES, and the eastern part of WS.

In winter, the area with a positive correlation between the NDVI and temperature accounted for 57.54%, mainly distributed in the southwestern part of WS; the negative correlation area accounted for 42.46%, located in NS, ES, and the eastern part of WS. The area with a positive correlation between the NDVI and precipitation accounted for 61.36%, mainly distributed in NS, CS, and parts of WS. The NDVI was positively correlated with sunshine duration, with the positive correlation area accounting for 84.14%, distributed in ES, SS, and the western part of WS.

3.3.2. Time Lag Effect of NDVI on Single Climatic Factors

Precipitation

The Pearson correlation coefficients of monthly NDVI and precipitation in Sichuan Province from 2000 to 2020 are shown in Figure 14. The response of the NDVI to precipitation in Sichuan Province mainly lags by one month. In various regions, when the cumulative time scale is short, the NDVI and precipitation show a positive correlation, and the lag time is shorter, indicating that vegetation growth is sensitive to precipitation. The response of the NDVI to precipitation in NS, ES, SS, and CS lags by 0 months, meaning the vegetation in these regions responds more quickly to precipitation. The correlation between the NDVI and precipitation in CS is the weakest. Due to less human activity in parts of WS, the vegetation there is more affected by natural factors, so the correlation between the NDVI and precipitation in this region is stronger compared to other areas, with a lag time of one month. The influence of precipitation on vegetation has a time lag, mainly due to the comprehensive effects of precipitation infiltration into the soil, plant physiological adjustment, growth cycle requirements, and climatic environmental factors. This leads to the effect of precipitation on vegetation not being immediate but requiring some time to become apparent.

Temperature

The Pearson correlation coefficients of the monthly NDVI and temperature in Sichuan Province from 2000 to 2020 are shown in Figure 15. The response of the NDVI to temperature in Sichuan Province mainly lags by one month. The lag time of the NDVI’s response to temperature in various regions of Sichuan Province is the same as that for precipitation, with the response in NS, ES, SS, and CS lagging by 0 months. However, the correlation between the NDVI and temperature is stronger than that with precipitation. The correlation between the NDVI and temperature in CS is relatively weak, while in WS, the correlation between the NDVI and temperature is stronger compared to other regions, with a main lag time of one month. Vegetation’s adaptation to temperature changes usually has a certain delay. After experiencing a rise in temperature, vegetation needs some time to adjust its physiological mechanisms to adapt to new climate conditions. Different types of vegetation have varying response speeds to temperature changes; some types may require more time to complete this adaptation process, showing a significant lag effect.

Sunshine Duration

The partial correlation coefficients of the monthly NDVI and sunshine duration in Sichuan Province from 2000 to 2020 are shown in Figure 16. The lag time of the NDVI’s response to sunshine duration in Sichuan Province is two months. The shorter response times of the NDVI to precipitation and temperature indicate that vegetation growth is more sensitive to these factors. The correlation between the NDVI and sunshine duration is weakest in WS, with the longest lag at five months. In NS, ES, CS, and SS, the response of the NDVI to sunshine duration mainly lags by zero months, with the strongest correlation between the NDVI and sunshine duration found in CS. The impact of sunshine duration on vegetation has a time lag primarily because plant growth and development require time to accumulate in response to light exposure. This involves energy conversion through photosynthesis, the combined effects of environmental factors, photoperiodism, regulation by plant hormones, and the response speeds of different growth organs.

3.3.3. Correlation Between NDVI and Hydrothermal Climatic Factors

The spatial distribution of correlation coefficients between the NDVI and the SPEI in Sichuan Province from 2000 to 2020, which passed the significance test, is shown in Figure 17. In spring, the NDVI and the SPEI are mainly positively correlated, accounting for 84.65% of the area, mainly distributed in the northern part of NS and parts of WS and SS. In summer, regions with a negative correlation between the NDVI and the SPEI are mainly located in the western and southern parts of WS, with a negative correlation area accounting for 35.82%. In autumn, the NDVI and the SPEI are predominantly positively correlated, covering 90.53% of the area, mainly found in the northern part of WS, as well as parts of ES and SS. In winter, the NDVI and the SPEI are mainly negatively correlated, with a negative correlation area accounting for 77.67%, primarily distributed in the central and western parts of WS. During the growing season and inter-annually, both the NDVI and the SPEI are predominantly positively correlated, with positive correlation areas accounting for 89.57% and 89.73%, respectively. The regions with strong positive correlations are similar for both, mainly located in NS, CS, and the northern part of WS. Overall, the correlation between the NDVI and the SPEI in various seasons across Sichuan Province is predominantly positive. However, in winter, due to lower temperatures and reduced precipitation, conditions reflecting drought or water shortage limit normal vegetation growth, leading to a predominantly negative correlation between the NDVI and the SPEI.

3.3.4. Time Lag Effect of NDVI on Multiple Climatic Factors

The Standardized Precipitation Evapotranspiration Index (SPEI) is an index used to measure the severity of drought, and based on previous analysis, precipitation and temperature have a stronger impact on vegetation growth compared to sunshine duration. The SPEI is based on the balance between precipitation and evapotranspiration to assess moisture conditions, so it was selected as the comprehensive climate factor to analyze its relationship with vegetation growth. The partial correlation coefficients of the monthly NDVI and the monthly SPEI in Sichuan Province from 2000 to 2020 are shown in Figure 18. The lag time of the NDVI response to the SPEI in Sichuan Province is one month. In NS, ES, and CS, the lag time is zero months. In WS and SS, the main lag time for the NDVI response to the SPEI is one month. Due to less human activity interference in parts of WS, climatic factors play a major role in vegetation growth, resulting in a higher correlation between the NDVI and the SPEI in WS compared to other regions of Sichuan Province.

4. Conclusions

This research, utilizing data from 2000 to 2020, including the NDVI, precipitation, temperature, sunshine duration, and SPEI, employs comprehensive methodologies to examine the spatiotemporal distribution, trends, and future shifts of the vegetation NDVI and its correlation with climatic factors in Sichuan Province. The study also quantifies the lag time of vegetation growth response to different climatic factors. The primary findings of the study are summarized as follows:

(1)

From 2000 to 2020, Sichuan Province saw an upward trend in precipitation and temperature, with geographical variations observed across the province. The province experienced periods of drought, most notably between 2006 and 2007 and between 2011 and 2012;

(2)

Since the implementation of the Natural Forest Protection Project in 2000, vegetation coverage in the province has consistently increased. The NDVI levels are highest during the summer and growing season, with spatial distribution showing higher values in the east and lower in the west. The area with no change in vegetation coverage accounted for the majority (55.2%) from 2000 to 2020. The future trend of the NDVI in Sichuan Province is characterized by continuous improvement, with specific regions showing shifts from falling to rising;

(3)

The NDVI in Sichuan Province generally exhibits a positive correlation with precipitation, temperature, and sunshine hours. However, the correlation with sunshine hours is primarily negative in autumn. The correlation between the NDVI and the SPEI is also primarily positive, except in winter;

(4)

The lag time of the NDVI response to precipitation, temperature, and the SPEI is one month, indicating a high sensitivity of vegetation growth to these climatic factors. The lag time for sunshine hours is two months, with the longest lag time of up to five months observed in western Sichuan.