Spatial–temporal variations of NDVI and its response to climate in China from 2001 to 2020

ABSTRACT

Vegetation plays an important role in global or regional environmental change. In this study, the spatial–temporal variations of NDVI and its response to climate in China and its seven sub-regions were investigated based on MODIS NDVI data, ERA5-land precipitation (PRE) and temperature (TEM) data from 2001 to 2020. The inter-annual growth rate of NDVI in China was 0.0021/yr in the past 20 years. The inter-annual growth rates of NDVI in seven sub-regions had significant differences at regional or seasonal scales. The ratio of improved vegetation area to the total studied area reached about 70%. In summer, vegetation degradation was concentrated in East China and Southwest China. The vegetation in Central China and South China improved more obviously in autumn than in the other seasons. The vegetation of Northeast China had a remarkable degradation in autumn and winter, especially in winter. The influence degree of PRE (q = 0.54, P < 0.01) was greater than that of TEM (q = 0.27, P < 0.01) in the control of the spatial distribution of NDVI. The interaction influence degree q of PRE $\cap$ TEM was about 0.71 in the last 20 years. However, the PRE and TEM played different roles in vegetation growth in seven sub-regions.

KEYWORDS:

• NDVI
• spatial–temporal variations
• climate change
• MK test
• Geodetector

1. Introduction

Ecological degradation is one of the key issues of global environmental change, which greatly restricts the development of human society and affects the evolution of living things. As an important part of maintaining the normal operation of the terrestrial ecosystem, vegetation plays an important role in global or regional environmental change (Zhu et al. 2016; Levine 2016; Law et al. 2002). With the influences of environmental conditions, the spatial–temporal variations in vegetation have demonstrated significant changes (de Jong et al. 2011). The hydrothermal conditions in climatic environment determine the phenology, productivity and distribution pattern of vegetation (Wen et al. 2017). At present, the monitoring of long time series, large-scale vegetation and factor analyses causing vegetation change have become a hot research topic (Zhang and Wu 2020; Guo et al. 2021; Kong et al. 2020). The related research work can provide an important reference for regional environmental protection, and social and economic development.

Satellite observations provide a useful and powerful database for evaluating the dynamics of terrestrial vegetation at regional, continental, and global scales because of its global coverage and short revisit interval (Piao et al. 2003; Bondeau et al. 1999; Fensholt et al. 2012; Cao et al. 2014; Liu, Su, et al. 2015). In the field of remote sensing application, Vegetation indice (VI) has been widely used for qualitative and quantitative evaluations of vegetation cover on a land surface and its growth vitality over wide areas. The Normalized Difference Vegetation Index (NDVI), one of the most commonly used indices, is calculated from satellite imagery in which the satellite's spectrometer or radiometric sensor measures and stores reflectance values for both the red and near-infrared (NIR) bands on two separate channels or images (Crippen 1990). Due to its stability, NDVI is successful to be a measure of vegetation to permit meaningful comparisons of seasonal and inter-annual changes in vegetation growth and activity (Huete et al. 2002; Olafsson and Rousta 2021).

China, a country with a vast territory and the largest population, has paid much attention to the improvement of the quality of ecological environment in recent years (Li et al. 2022). China's achievements in ecological improvement are of far-reaching influence on the global. Climate diversity and intensive human activity have caused complex vegetation activity in China (Liu, Li, et al. 2015). Some researchers had investigated the variations of NDVI and their driving factors, mainly focused on localized areas (Cao et al. 2014; Qi et al. 2019; Eckert et al. 2015; Pang, Wang, and Yang 2017; Li, Wang, et al. 2021; Fensholt et al. 2009; Fensholt et al. 2012; Chuai et al. 2013; Gao et al. 2022; Hou et al. 2015; Rousta et al. 2020). Linear fit, residual trend analysis and correlation analysis methods were mostly used in these studies to qualitatively analyze the spatial–temporal changes of NDVI, of course, also used in other studies on a national or global scale (Gao, Jiao, and Wu 2019; Fang et al. 2004; Liu, Li, et al. 2015; Piao et al. 2003; Eckert et al. 2015; Fensholt and Proud 2012; Fensholt et al. 2009; Fensholt et al. 2012). Some studies focus on the trend analysis of non-stationary NDVI time series through decomposing the time series into different components (Rhif et al. 2021; Rhif et al. 2022; Abbes et al. 2018; Wen et al. 2017). At present, the successful and effective applications of the Mann-Kendall test (Kendall 1975; Mann 1945) and Geodetector (Wang et al. 2010; Wang, Zhang, and Fu 2016) in earth science provide new ways to investigate the temporal trend of NDVI and spatial heterogeneity of factors, such as the NDVI characteristics in the Hanjiang River Basin (Chen et al. 2020) and Loess Plateau (Nie et al. 2021; Li, Wang, et al. 2021). However, long-term systematic analysis of NDVI over the Chinese regions is still rare, especially in exploring the spatial heterogeneity of the response of vegetation to climate change.

This study concentrated on temporal and spatial variations of NDVI and its responses to temperature (TEM) and precipitation (PRE) based on data from Satellite observations over the past two decades in China. The purposes of this study are (1) to examine the inter-annual growth and seasonal growth of NDVI from 2001 to 2020 in China and its seven sub-regions; (2) to demonstrate the spatial distribution of NDVI trends in the past 20 years; (3) to detect the relationships between the spatial–temporal variations of NDVI and climatic factors. Figure 1 summarizes the research route applied in this study. The detailed descriptions will be elaborated on in the following sections. Our work will contribute to a better understanding of climate-driven vegetation variations, and promoting the implementation and management of ecological reconstruction programs.

Figure 1. A flow chart of the research route applied in this study.

2. Data and methods

2.1. Study area

China is located in the east of Asia and on the West Bank of the Pacific Ocean. It has a vast territory with a total land area of about 9.6 million square kilometers, which spans nearly 50 degrees from north to south. Most of Chinese territory is in the temperate zone, and a small part is in the tropical zone, without the frigid zone. As one of the developing countries, China has the largest population in the world. Affected by the monsoon climate, the winter is cold and dry, while the summer is warm and rainy. The PRE is unevenly distributed, mainly concentrated in summer and autumn in most areas and gradually decreased from the southeast coast to the northwest inland. In winter, the TEM decreases from the south to the north, and the TEM difference between the north and the south is obvious. In summer, high TEM is common in most parts of the country except the Qinghai-Tibet Plateau. The rich vegetation resources are formed under diverse climate and natural conditions. Physical geographical regionalization was used in the following analysis. In other words, Chinese region is divided into seven sub-regions (An et al. 2015), as shown in Figure 2. Some works would be performed at a regional scale in this study. The distribution of 20-year averaged NVDI is also presented in the figure and reflected a large regional difference of gradually decreasing from the southeast to the northwest.

Figure 2. Physical geographical regionalization of the study area. From north to south, east to west, the seven sub-regions are Northeast China, North China, Northwest China, East China, Central China, Southwest China, and South China. The pixel color represents the value of the 20-year average NDVI.

2.2. Data sources and data processing

The datasets used in this study consist of NDVI data, PRE and TEM data. NDVI data set series from March 2001 to February 2021 were obtained from Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra satellite from the National Aeronautics and Space Administration (NASA). We used 16-day intervals of NDVI product MOD13A2 with a spatial resolution of 1 km × 1 km (Fensholt et al. 2009). MODIS has been in orbit from March 2000 to the present time (Shi et al. 2021). Compared with the other NDVI dataset, the MODIS NDVI product is more reliable and suitable for reflecting the long-term vegetation dynamics over a large spatial scale, considering the spatial–temporal resolution, time span and on-orbit status (Zhu et al. 2019; Sun et al. 2015). The spatial resolution was bilinearly resampled from 250 m to 1 km (Shi et al. 2021). Nevertheless, the resolution of 1 km was used in this study instead of 250 m, due to the great data processing burden on the study area of 9.6 million square kilometers. There are 8,126,518 pixel values in each 16-day, 1 km × 1 km NDVI datum.

The MODIS Reprojection Tools (MRT) were used to cut and splice the NDVI bands of MOD13A2. These NDVI datasets were re-projected to the WGS-84 coordinate system using the ArcGIS 10.8 software. Climate year is adopted, that is, started in March of a given year and lasted until the next February over a continuous period from March 2001 to February 2021. The 16-day NDVI datasets were synthesized to the monthly values pixel-by-pixel using the Maximum Value Composite method (Sun et al. 2015; Nie et al. 2021; Chu et al. 2019). The S-G (Savitzky Golay Filter) algorithm was used for monthly NDVI image smoothing and de-noising, reducing the influence of rain and snow atmosphere (Savitzky and Golay 1964). The spring period is from March to May (Mar-Apr-May), summer is from June to August (Jun-Jul-Aug), autumn is from September to November (Sept-Oct-Nov) and winter is from December to February (Dec-Jan-Feb). Then seasonal and annual NDVI data were composited by monthly NDVI images at the regional scale and the pixel scale by averaging. The value of NDVI ranges from −1 to 1 (dimensionless), increasing positive NDVI values indicate increasing amounts of green vegetation. NDVI values near zero and decreasing negative values indicate non-vegetated features such as barren surfaces (rock and soil) and water, snow, ice and clouds (Schnur, Xie, and Wang 2010; Chuai et al. 2013; Olafsson and Rousta 2021). Therefore, the pixels with NDVI > 0 in vegetation areas were extracted in this study.

TEM and PRE raster datasets during 2001–2020 were obtained from ECMWF (European Centre for Medium-range Weather Forecasts) ERA5-Land (the land component of the 5th generation of European ReAnalysis) products (Joaquín et al. 2021), monthly averaged 2 m TEM and total PRE were adopted with spatial resolution of 9 km × 9 km. 2 m TEM refers to the dew point TEM 2 m high on the earth surface with the unit of Kelvin, which was converted into Celsius ($1{}^{\circ }\mathrm{C}=1\mathrm{K}-273.15$) in this study. Total PRE measured in m was converted to mm.

In addition, the NDVI and PRE/TEM data were resampled into a grid dataset of 50 km × 50 km for the subsequent detection of their relationships.

2.3. Methods

2.3.1. Linear regression analysis

The linear regression analysis method was applied to detect the temporal variations (inter-annual growth/seasonal growth) of NDVI from 2001 to 2020 at the regional scale, and the ordinary least squares method was employed to determine the parameters in the regression model (Nie et al. 2021; Hou et al. 2015). The linear regression equation is expressed as follows. (1) $NDV{I}_i=a+b\times i+{\epsilon }_{i}i=1,2,\dots ,n$(1) where $NDV{I}_i$ denotes the average NDVI of the i-th year in the year series. The time span in this study is 2001–2020, so n equals to 20. a and b are the intercept and slope of the linear regression equation, respectively; ${\epsilon }_i$ expresses the error term. b represents the inter-annual growth rate of the NDVI in the study period. The vegetation displays an increasing/ decreasing trend when the value of b is greater/lesser than zero (Gao et al. 2022; Sun et al. 2015). As one of the significance test methods for linear regression, the Pearson correlation coefficient ($r$) is derived to determine the significance of the NDVI linear trend. Since the sample size is less than 30 (n < 30), the unbiased correlation coefficient is more convincing for correlation analysis, as denoted as $r^{\ast }$. (2) $r^{\ast }=r\left[1+{\displaystyle \frac{1-r^2}{2(n-4)}}\right]$(2) where $r$ is the Pearson correlation coefficient and n is the size of the year series. When $|r^{\ast }|\ge 0.444$, the linear trend has passed the significance test with 95% confidence intervals ($P<0.05$).

2.3.2. Theil-Sen Median and Mann-Kendall test

Theil-Sen Median trend analysis and the Mann–Kendall test were combined to investigate the trend characteristic of NDVI in the study period (Cao et al. 2014). Theil-Sen Median method is a robust trend calculation method of non-parametric statistics, which is also known as Sen-slope estimation. This method is insensitive to measurement error and outlier data and has high computational efficiency, widely used in the long-time series trend analysis of meteorology and hydrology (Zhang et al. 2015; Li, Wang, et al. 2021; Ali et al. 2019). The Sen-slope is defined as follows. (3) $S_{NDVI}=\mathrm{Median}\left({\displaystyle \frac{NDV{I}_j-NDV{I}_i}{j-i}}\right)1\le i<j\le n$(3) where $S_{NDVI}$ means the trend of NDVI (increase or decrease), $\mathrm{Median}$ means the median value, $i$ and $j$ represent the years, $NDV{I}_i$ is the NDVI value of the year $i$. n is the size of the year series.

Mann-Kendall (MK for short) is a non-parametric test method, which does not require samples to follow a certain distribution compared with other parameter test methods. Sen-slope value is calculated first, and then the MK method is used to determine the trend significance. The test statistic $S$ is defined as follows. $Sgn()$ is a sign function. (4) $S={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^{n}Sgn(NDV{I}_j-NDV{I}_i)$(4) (5) $Sgn(NDV{I}_j-NDV{I}_i)=\{\begin{array}{c}1,NDV{I}_j-NDV{I}_i>0\\ 0,NDV{I}_j-NDV{I}_i=0\\ -1,NDV{I}_j-NDV{I}_i<0\end{array}$(5) The statistic $S$ is approximately normally distributed with the variance as follows. (6) $Var(S)={\displaystyle \frac{1}{18}[n(n-1)(2n+5)}$(6) The standardized test statistic Z is computed by (7) $Z=\{\begin{array}{c}{\displaystyle \frac{S-1}{\sqrt{Var(S)}}S>0}\\ 0S=0\\ {\displaystyle \frac{S+1}{\sqrt{Var(S)}}S<0}\end{array}$(7) The statistic Z follows the standard normal distribution with a mean of zero and a variance of one. The critical value is 1.96 at the significance level $\alpha =0.05$, that is, if $|Z|\ge 1.96$ the trend has passed the significance test with confidence intervals of 95%.

2.3.3. Geodetector

Geodetector (geographical detector) is a spatial variance analysis tool to assess the relative importance of different components controlling a geographic phenomenon, which is intrinsically immune to collinearity (Wang et al. 2010; Chen et al. 2022; Fang et al. 2021). This method has been successfully applied to determine controlling factors in many geoscience fields (http://www.geodetector.cn/). The Geodetector is classified into four types, namely Factor_detector, Interaction_detector, Risk_detector and Ecological_detector. In this study, Factor_detector, Interaction_detector and Risk_detector are used to reveal the response of vegetation to climate from multiple perspectives.

1. Factor_detector can detect the spatial differentiation of Y (NDVI) and the influence degree of Y changes caused by influence factor X (PRE, TEM), measured by $q$ statistic value, is determined by the local and global variances (Wang et al. 2010). The formulae are as follows. (8) $q=1-{\displaystyle \frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}=1-{\displaystyle \frac{SSW}{SST}}}$(8) (9) $SSW={\sum }_{h=1}^L{N}_h{\sigma }_h^2,SST=N{\sigma }^2$(9)

where $q$ is the influence degree of influence factor X to the change of Y; $L$ is the partition number of Y or X; $N_h$ is the number of samples in partition $h$ and $N$ is the total number in the whole area; ${\sigma }_h^2$ and ${\sigma }^2$ are the variances of partition $h$ and the whole region, respectively. The $q$ statistic is between 0 and –1. The value is larger and the influence degree of X on Y is greater. If $q=0$, it indicates that the distribution of X and Y is completely independent, whereas $1$ means factor X completely controls the spatial patterns of Y (Chen et al. 2022). A transformation of $q$ value satisfies a noncentral F distribution, and the P value acquired by geographic detector software is used to test the significance of $q$.

1. Interaction_detector was used to evaluate the interaction of different factors (PRE and TEM). The aim is to determine whether the interaction of two factors could enhance or weaken the influence on the dependent variable. Specifically, overlaying the individual factor strata $X1$ and $X2$ can generate new factor strata named $X1\cap X2$, where $\cap$ refers to the interaction between factor $X1$ and $X2$. The interactive $q$ value of $X1\cap X2$ is referred to as $q(X1\cap X2)$ (Wang et al. 2021).

2. Risk_detector is used to detect the suitable climate conditions in favor of vegetation growth. Risk_detector is based on the $t$-statistic. The calculation steps are as follows. First, calculate the average values of NDVI corresponding to different partitions of a factor (PRE or TEM). Next, determine whether there is a significant difference between every two values based on the $t$-test. Finally, find out the partition corresponding to the maximum value. If it is significantly different from other values, means that this partition is the most suitable. Otherwise, find the second largest value and compare whether it is significantly different from the values other than the maximum value. If it is true, there are two most suitable partitions. Otherwise, repeat the above step until the conditions are satisfied (Chen et al. 2020).

3. Results

3.1. Temporal variations from 2001 to 2020

3.1.1. Inter-annual variation

The inter-annual variations of NDVI were investigated by linear regression analysis in China and its seven sub-regions during 2001–2020. The slope of the fitted line was expressed as an inter-annual growth rate. The time series of the annual averaged NDVI over China is shown in Figure 3. The fluctuation range of NDVI was 0.49 ∼ 0.55 in the study period. The inter-annual growth rate was 0.0021/yr and $r^{\ast }$ is 0.924, greater than the critical value of 0.444 indicating a significantly increasing ($P<0.05$). In the last 20 years, the average NDVI was about 0.525, the minimum value was about 0.496 in 2001, and the maximum value was about 0.546 in 2020.

Figure 3. Time series of annual averaged NDVI over China. The line represents the linear fit of NDVI over the time series from 2001 to 2020. The text on the upper-left labels the regression equation and unbiased correlation coefficients.

Figure 4 shows the variations of NDVI in seven sub-regions of China. From 2001 to 2020, there was a significant difference in these time series of annual averaged NDVI (Figure 4(a)). Seven sub-regions could be divided into three clear groups according to the serial interval of NDVI. The first group of NDVI kept a high value (0.7 ∼ 0.9) during the study period (Figure 4(b)), including Northeast China, Central China, East China and South China. The second group of NDVI was moderate (0.4 ∼ 0.6), including North China and Southwest China (Figure 4(c). The third group of NDVI maintained a low scope (0.2∼0.3), located in Northwest China (Figure 4(d)). From 2001 to 2020, the NDVI trends of seven sub-regions were different but approximately similar to each other, demonstrating the consistent growth.

Figure 4. (a) Time series of annual averaged NDVI in seven sub-regions, which can be divided into three groups. (b) The first group includes Northeast China, South China, Central China and East China. (c) The second group includes Southwest China and North China. (d) The third group includes Northwest China. The line represents the linear fit of NDVI over the time series from 2001 to 2020.

As shown in Table 1, the inter-annual growth rates of Northeast China, North China, South China and Northwest China were all greater than 0.0024/yr. Although the NDVI in Northwest China remained in a low range, it still showed an obvious increase. The inter-annual growth rates of Central China and Southwest China were 0.0018/yr and 0.0013/yr, respectively, indicating a certain growth. The NDVI in East China was in a higher range, but the inter-annual growth rate was only 0.0006/yr, and $r^{\ast }$ was only 0.378, having not passed the significance test with 95% confidence intervals. Hence, the growth of NDVI in East China was inconspicuous, which may be related to the development of urbanization in East China.

Table 1. The regression equation and unbiased correlation coefficient of inter-annual NDVI evolution in seven sub-regions.

Sub-regions	Regression equation	$r^{\ast }$	P
Northeast China	y = 0.0024x + 0.8079	0.906	< 0.05
North China	y = 0.0029x + 0.4875	0.729	< 0.05
East China	y = 0.0006x + 0.7454	0.378	> 0.05
South China	y = 0.0029x + 0.7596	0.959	< 0.05
Central China	y = 0.0018x + 0.7647	0.798	< 0.05
Northwest China	y = 0.0024x + 0.2654	0.926	< 0.05
Southwest China	y = 0.0013x + 0.5305	0.928	< 0.05

3.1.2. Seasonal variation

Significant seasonal differences were presented in the vegetation dynamics. The NDVI time series of four seasons over China all showed obvious increases from 2001 to 2020 (Figure 5). In spring, the NDVI fluctuated between 0.26 and 0.33, and the inter-annual growth rate was 0.0023/yr. The mean value of 20-year NDVI was 0.295, the minimum value was 0.267 in 2001 and the maximum value was 0.322 in 2018. The mean value of summer NDVI was the largest (0.467) among the four seasons, with a fluctuation range of 0.44∼0.48, and the inter-annual growth rate was 0.0018/yr. The rapid growth of summer NDVI appeared in the period of 2009–2013, the maximum value of NDVI was 0.488 in 2013, and the minimum value was 0.444 in 2001. In autumn, the mean value of NDVI was 0.331, with a range of 0.30∼0.36 and an inter-annual growth rate of 0.0021/yr. The minimum value of autumn NDVI was 0.306 in 2002 and the maximum value was 0.353 in 2019. The time series of 20-year NDVI fluctuated sharply in winter, and the inter-annual growth rate was 0.0022/yr, with a minimum value of 0.183 in 2004 and a maximum value of 0.247 in 2020, showing a non-sustainable growth trend.

Figure 5. Time series of seasonal averaged NDVI in (a) Spring, (b) Summer, (c) Autumn and (d) Winter. The line represents the linear fit of NDVI over the time series from 2001 to 2020. The text on the upper left labels the regression equation and unbiased correlation coefficients.

Seasonal differences were also presented in seven sub-regions. Figure 6 displays that the NDVI values in seven sub-regions of China increased except the NDVI of East China in summer (Figure 6(a)). The inter-annual growth rates of Northeast China, North China and Northwest China in summer were higher than those in other seasons. However, the NDVI growths of East China, South China and Central China were more significant in winter, but the growth of Southwest China was more significant in spring, indicating the spatial heterogeneity of growth trend in seven sub-regions at seasonal scale. The summer growth was not significant in Southwest China, even with a slight decrease in East China (Figure 6(b)). In terms of the averaged inter-annual growth rate of four seasons, the growth rates of South China, Central China and East China were in the top three, while the growth rates of Northwest China and Southwest China were slightly smaller than those of North China and Northeast China. On the whole, the averaged inter-annual growth rates of seven sub-regions showed a decreasing trend from south to north and from east to west.

Figure 6. Inter-annual growth rates and inter-annual correlations of NDVI during 2001–2020 at regional and seasonal scales. (a) Inter-annual growth rates of seven sub-regions in four seasons, the black dot represented the averaged inter-annual growth rate of four seasons for each sub-region. (b) Unbiased correlation coefficients between NDVI series and Year series of seven sub-regions in four seasons were used to determine the trend significantly.

3.2. Spatial distribution of NDVI trend

Sen-slope estimation and MK test were carried out over China pixel-by-pixel during the study period, and the analysis results were used to determine the growth trend of vegetation. When $S_{NDVI}$ is greater than 0, the vegetation shows an improving trend, when $S_{NDVI}$ is less than 0, the vegetation shows a worsening trend. The Z statistic of the MK test illustrates the significance of improvement or degradation. By synthesizing $S_{NDVI}$ and Z, the trends were divided into five trend levels (Wang et al. 2022; Cao et al. 2014): Severe degradation (S_NDVI < −0.0005, Z < −1.96), Slight degradation (S_NDVI < −0.0005, −1.96 ≤ Z ≤ 1.96), Stable (−0.0005 ≤ S_NDVI ≤ 0.0005, −1.96 ≤ Z ≤ 1.96), Slight improvement (S_NDVI > 0.0005, −1.96 ≤ Z ≤ 1.96) and Significant improvement (S_NDVI > 0.0005, Z >  1.96). The spatial distribution of 20-year NDVI trend over China is shown in Figure 7. Meanwhile, Table 2 shows the area statistics of different vegetation trends in China and seven sub-regions.

Figure 7. Spatial distribution of the 20-year NDVI trend over China. Five trend levels were displayed with different colors.

Table 2. Area Statistics of five trend levels in China and seven sub-regions. A stands for the area (unit: $\mathrm{k}{\mathrm{m}}^2$) and P stands for percentage (unit: %).

Trend levels	China		Northeast China		North China		East China
	A/$\mathrm{k}{\mathrm{m}}^2$	P/%	A/$\mathrm{k}{\mathrm{m}}^2$	P/%	A/$\mathrm{k}{\mathrm{m}}^2$	P/%	A/$\mathrm{k}{\mathrm{m}}^2$	P/%
Severe degradation	254815.31	2.65	12253.15	0.13	33947.45	0.35	71219.05	0.74
Slight degradation	861698.13	8.98	40971.48	0.43	133912.35	1.39	100779.78	1.05
Stable	1803316.97	18.78	74137.55	0.77	258791.60	2.70	59568.99	0.62
Slight improvement	2846945.69	29.66	264963.64	2.76	580585.22	6.05	179411.75	1.87
Significant improvement	3833223.91	39.93	528560.77	5.51	682068.60	7.10	181102.54	1.89
Trend levels	South China		Central China		Northwest China		Southwest China
	A/$\mathrm{k}{\mathrm{m}}^2$	P/%	A/$\mathrm{k}{\mathrm{m}}^2$	P/%	A/$\mathrm{k}{\mathrm{m}}^2$	P/%	A/$\mathrm{k}{\mathrm{m}}^2$	P/%
Severe degradation	10966.81	0.11	35083.02	0.37	38197.77	0.40	53143.26	0.55
Slight degradation	26289.23	0.27	75426.29	0.79	169902.40	1.77	314401.04	3.28
Stable	26793.00	0.28	57692.72	0.60	907745.57	9.46	418582.75	4.36
Slight improvement	119391.65	1.24	210983.96	2.20	724320.43	7.55	767271.08	7.99
Significant improvement	247306.66	2.58	310451.08	3.23	1278866.50	13.32	604845.02	6.30

The vegetation in most of China has improved significantly in the past 20 years (Figure 7). Combined with Table 2, it can be found that the areas of slight improvement and significant improvement of vegetation in China took up 29.66% and 39.93%, respectively, with an improvement percentage of nearly 70%, which is much higher than the degraded area (about 11.6%), indicating a good growth trend. In terms of regional statistics, the improved vegetation areas were mainly distributed in Northwest China, Southwest China and North China, which accounted for 20.87%, 14.29% and 13.15% of the total area, respectively. The vegetation degradation (Slight degradation and Severe degradation) was obvious in Southwest China, Northwest China, East China and North China accounting for 3.83%, 2.17%, 1.79% and 1.74%, respectively. Certainly, the percentage of vegetation degradation was smaller than that of vegetation improvement within the same sub-region. Stable areas concentrated in the Taklimakan desert in Northwest China and Ala Shan Desert in North China, with a very low NDVI value (Figure 2).

Furthermore, four seasonal components of NDVI were extracted from the original monthly dataset from 2001 to 2020. Sen-slope estimation and MK test were carried out at a seasonal scale. The growth trends of NDVI during 2001∼2020 in spring, summer, autumn and winter are displayed in Figure 8. In summer, vegetation degradation areas were concentrated in East China and Southwest China. The vegetation severe degradation area in East China accounts for 2.3% of China's total area in summer, much higher than that in other seasons. The land use change caused by urbanization had a destructive effect on vegetation in East China and its surrounding areas (Yuan et al. 2019; Zheng et al. 2021; Haas and Ban 2014). The degradation percentage of Southwest China in summer was 6.6% and it was probably affected by natural disasters caused by heavy summer rains. Compared with other seasons, the vegetation in Central China and South China improved obviously in autumn. As the main grain-producing area in China, the crop-growing season of Northeast China lasts from April to October under cold weather. The vegetation in autumn and winter had an obvious degrading trend, especially in winter.

Figure 8. Spatial distribution of the 20-year seasonal NDVI trend over China. (a) Spring, (b) Summer, (c) Autumn and (d) Winter. Unvegetated areas were mostly ice-snow-covered, cloud-sheltered areas or water areas.

3.3. Association between NDVI and climatic factors

Climate plays an important role in vegetation dynamics; PRE and TEM are often taken as the key factors. Geodetector is adopted to explore the spatial–temporal differences and relationships between NDVI and PRE/TEM. The distributions of 20-year averaged PRE/TEM with spatial resolution of 9 km × 9 km are shown in Figure 9. The entire north-western inland area had sparse PRE throughout the year. The TEM distribution was very uneven due to the Geological and topographical conditions. The grid dataset of 50 km × 50 km mentioned in Section 2.2 was used for the input of the Geodetector. The points in Figure 9 (N = 427) are located at the central positions of the grids of 50 km × 50 km. The mean of all pixel values in a grid was computed as the value of the grid point. Specifically, each grid point corresponds to three values, namely NDVI, PRE and TEM, where NDVI is a dependent variable, PRE and TEM are independent variables, and the output is influence degree, represented as $q$ value. The interval of PRE/TEM was divided into nine partitions according to natural breaks (Jenks) (Slocum 1999).

Figure 9. Distribution of 20-year average PRE (a) and TEM (b) over China. The points are located at the central positions of the grids with a spatial resolution of 50 km × 50 km.

The influence degree $q$ was obtained by Factor_detector and reflected the quantitative relationship between NDVI and PRE/TEM for each year. The $q$ series of 20 years is displayed in Figure 10. The $q$ value of PRE was in the range of 0.46 ∼ 0.61, with a mean of 0.54 and passed the significance test (P < 0.01). The $q$ of TEM was in the range of 0.23 ∼ 0.32, with a mean of 0.27 and passed the significance test (P < 0.01). The $q$ value of PRE was apparently higher than the $q$ value of TEM. Therefore, it can be found that PRE was the main influencing factor and TEM was the secondary one in the control of spatial distribution of NDVI. The large fluctuation of $q$ value indicated that vegetation was affected by many other influencing factors besides climate in 20 years.

Figure 10. The influence degree $q$ of PRE/TEM acting on NDVI over China from 2001 to 2020.

Interaction_detector was used to determine the combined effects of PRE and TEM on NDVI. The interaction influence degree $q$ from 2001 to 2020 is shown in Figure 11. The mean of $q$ values of the combined effects of PRE and TEM on NDVI (denoted by PRE $\cap$ TEM) was about 0.71 in the last 20 years.

Figure 11. The interaction influenced degree $q$ of PRE $\cap$ TEM in China from 2001 to 2020. The black dotted line represents the mean value of $q$.

In the meantime, Factor_detector and Interaction_detector were carried out in each sub-region for comparison of influence differences. The influence degree $q$ of seven sub-regions is demonstrated in Figure 12, the values were the 20-year average $q$. In Northeast China, East China and Central China, TEM became the main influencing factor on NDVI change, while in North China, South China and Northwest China, PRE was the main influencing factor. The $q$ value of PRE in Northwest China (0.604) was much higher than that of TEM (0.174). The $q$ values of TEM and PRE in Southwest China slightly different from each other, that is, 0.663 and 0.599, respectively. PRE and TEM played different roles in seven sub-regions under different geographic and geomorphic conditions. Moreover, the $q$ value of PRE $\cap$ TEM in seven sub-regions was all above 0.72, indicating that the response of vegetation to climate was controlled by the interaction of PRE and TEM.

Figure 12. The 20-year averaged $q$ values of PRE, TEM and PRE $\cap$ TEM in seven sub-regions.

Risk_detector was adopted to explore the suitable PRE/TEM of vegetation growth. The detector analysis was based on the partitions of PRE and TEM, as displayed in Figure 9. A few grid points were located in partitions of 157.57 ∼ 256.49 mm and 256.49 ∼ 415.42 mm, as shown in Figure 9(a), so the two partitions were ignored and PRE was divided into seven partitions labeled as 1 ∼ 7. There were nine TEM partitions, as shown in Figure 9(b), labeled as 1 ∼ 9. The detector results of PRE and TEM are listed in Table 3. Table 3(a and c) lists the mean values of NDVI in each PRE/TEM partition. Table 3(b and d) shows the statistics of whether a significant difference of NDVI existed in two of the PRE/TEM partitions, with ‘Y’ meaing Significant, and ‘N’ meaning Non-significant. As seen from Table 3(a), sorted by NDVI from largest to smallest, the corresponding PRE partition sequence of 7 > 5 > 6 > 4 > 3 > 2 > 1 indicated that there was a weak positive correlation between NDVI and precipitation, the mean value of NDVI was the largest (0.85) in the 7th partition (102.43 ∼ 157.57 mm). At the same time, Table 3(b) shows that there was a significant difference between partition 7 and the other PRE partitions, which indicates that the most suitable PRE for vegetation growth is 102.43 ∼ 157.57 mm. Similarly, the order of TEM partition sequence was 9 > 8 > 4 > 3 > 7 > 5 > 6 > 2 > 1. The mean value of NDVI was the largest (0.79) in the 9th TEM partition (> 20°C) (Table 3(c)). Furthermore, there was a significant difference between TEM partition 9 and the other TEM partitions (Table 3(d)), indicating better vegetation growth when TEM was above 20°C.

Table 3. The Risk_detector results of PRE and TEM.

(a) PRE
P	1	2	3	4	5	6	7
M	0.13	0.46	0.69	0.72	0.77	0.75	0.85
(b) PRE
	1	2	3	4	5	6	7
1
2	Y
3	Y	Y
4	Y	Y	N
5	Y	Y	Y	N
6	Y	Y	N	N	N
7	Y	Y	Y	Y	Y	Y
(c) TEM
P	1	2	3	4	5	6	7	8	9
M	0.18	0.32	0.51	0.68	0.45	0.37	0.46	0.74	0.79
(d) TEM
	1	2	3	4	5	6	7	8	9
1
2	Y
3	Y	Y
4	Y	Y	Y
5	Y	Y	N	Y
6	Y	N	Y	Y	N
7	Y	Y	N	Y	N	N
8	Y	Y	Y	N	Y	Y	Y
9	Y	Y	Y	Y	Y	Y	Y	Y

‘P’ stands for the partition of PRE and TEM in (a) and (c).

‘M’ stands for the mean value of NDVI in each partition in (a) and (c).

(b) and (d) represents the statistical significance of the NDVI difference between two PRE/TEM partitions, where "Y" and "N" represent Significant and Non-significant, respectively.

4. Discussion

On the basis of the above quantitative analysis, the results were further discussed according to the distributions of 20-year average NDVI, PRE and TEM. Several representative areas were indicated by rectangular boxes in Figure 13 for a better graphic description based on Figure 2. The NDVI of the whole study area gradually decreased from the southeast to the northwest in the last 20 years (Figure 2), and its distribution was consistent with that of PRE (Figure 9). The quantitative analysis of Factor_detector also verified that PRE had a greater influence on vegetation growth than TEM (Figure 10).

Figure 13. Several representative areas were indicated by rectangular boxes based on Figure 2.

NDVI area was mainly located in Tarim Basin (Region A), Qaidam Basin (Region B) and Alexa Plateau (Region C). Tarim Basin belongs to a warm and dry climate, basically in Gobi and Taklimakan deserts with an annual average TEM of 9 ∼ 11°C (Abbas et al. 2021). The study in the Tarim River Basin indicated that the vegetation coverage of the growing season was positively correlated with surface TEM and rain, although the latter had a better correlation coefficient (Jiang et al. 2022). Qaidam Basin is characterized by a plateau continental climate, with dry and cold winters and arid and windy summers. The growing season is from May to October. The mean annual TEM is 4.5°C, and the mean TEM in the growing season can reach 12°C (Han et al. 2021). The Alxa Plateau is in a dry non-monsoon area, and the Qinghai-Tibet Plateau in its south blocks the humid airstream from the Indian Ocean, following unusually a dry and hot climate (Wang et al. 2022). The drought activity on account of less PRE and evapotranspiration of higher air TEM in Tarim Basin, Qaidam Basin and Alxa Plateau, results in very sparse vegetation on the surface (Yang et al. 2022; Jiang et al. 2018; Jin et al. 2021). Take the Taklamakan Desert as an example, a moderate negative correlation is displayed between the TEM and the NDVI under drought conditions (Abbas et al. 2021). However, in recent years, the climate transition from warm-dry to warm-humid with TEM increasing and more PRE in Qaidam Basin under global warming, we can also find that the vegetation in Qaidam Basin has been continuously improved, as shown in Figure 7. An increased PRE was generally followed by an increase in the vegetation coverage in dry lands (Guo et al. 2021; Shi et al. 2021; Gao et al. 2022; Song et al. 2021; Han et al. 2021). In the Qinghai-Tibet Plateau (Region D), vegetation species and coverage decreased from southeast to northwest (Zhang, Fang, and Liu 2022). The Kunlun Mountain area in the northwest of the Plateau was the typical area of low vegetation caused by drought and low TEM and showed continuous degradation (Figure 7). Existing studies had presented consistent conclusions with ours (Pang, Wang, and Yang 2017). The reduction in vegetation in recent years may be attributable to permafrost degradation caused by climate warming and increasing anthropogenic activity (Yang et al. 2010; Chen et al. 2013; Chen et al. 2014; Xu, Wang, and Zhang 2016).

Tianshan Mountain (Region E) is a special region in the northwest of the study area. There is a dense primeval forest distributed in the northern piedmont of the Tianshan Mountains that benefited from low TEM, high PRE and small evapotranspiration. While in the southern piedmont, the vegetable is mostly flat grasslands with a climate of high TEM, low PRE and high evapotranspiration (Liu, Yan, and Jiang 2021; Zheng et al. 2021). But there has been an obvious trend of degeneration in recent 20 years (Figure 7). Li et al. suggested that soil moisture deficit played a greater role in the recent browning of vegetation on Tianshan Mountain (Li, Chen, et al. 2021). The southeast valley of the Qinghai-Tibet Plateau (Region F), also known as the Brahmaputra Valley, has abundant vegetation with a climate of abundant PRE and an annual average TEM of 20°C, affected by the warm and humid southwest monsoon over the Indian Ocean. Climate change was the most crucial driver of vegetation change in this area, and human disturbance was responsible for abrupt changes in vegetation (Liu, Liu, et al. 2021). In addition, studies on Loess Plateau (Region G) vegetation indicated that drought was the primary climate factor for the vegetation decreasing in the Loess Plateau, human activities have contributed more to the improvement and degradation of vegetation (Zheng et al. 2021) and the spatial pattern of the cumulative afforestation driven by the Grain for Green Project matched quite well with the NDVI change trend (Shi et al. 2021; Zhao et al. 2019).

The PRE-induced change of NDVI is mainly determined by the water use efficiency of vegetation. Climate warming may exacerbate evapotranspiration and reduce soil moisture, limiting vegetation growth in the case of limited PRE (Rousta et al. 2022; Zhao et al. 2022). The increase of TEM can improve photosynthetic efficiency and prolong vegetation growth season, thus promoting vegetation growth (Gao et al. 2022; Dragoni et al. 2011). The above analysis reconfirmed that the response of vegetation to climate was dominated by the interaction of PRE and TEM (Figure 12).

5. Conclusion

Based on MODIS NDVI data, ERA5-land PRE and TEM data from 2001 to 2020, this study presented the spatial–temporal variations of NDVI and its response to PRE/TEM in China and its seven sub-regions. The following conclusions could be drawn based on the above analysis.

1. The inter-annual variations of NDVI were investigated by linear regression analysis at the regional scale. The 20-year averaged NDVI over China was about 0.525 and the inter-annual growth rate of NDVI was 0.0021/yr in the 20 years. Seven sub-regions could be divided into three clear groups according to the serial intervals of NDVI: the NDVIs of Northeast China, Central China, East China and South China were in the high range (0.7 ∼ 0.9), the NDVIs of North China and Southwest China were in medium levels (0.4 ∼ 0.6), and the NDVI of Northwest China was the lowest (0.2 ∼ 0.3). The inter-annual growth rates of Northeast China, North China, South China and Northwest China were all greater than 0.0024/yr. The inter-annual growth in East China was the lowest. The inter-annual growth rates of NDVI in China at a seasonal scale were Spring (0.0023/yr) > Winter (0.0022/yr) > Autumn (0.0021/yr) > Summer (0.0018/yr).

2. The spatial distribution of vegetation growth trend was analyzed with Sen-slope estimation and MK test at pixel scale. The improved area of vegetation accounted for nearly 70% of the total study area and was much higher than the degraded area (about 11.6%), showing a good growth trend. The analysis with seasonal scale indicated that vegetation degradation was concentrated in East China and Southwest China in summer. Among all seasons, the vegetation in Central China and South China improved the most obviously in autumn. The vegetation of Northeast China had a remarkable degradation in autumn and winter, especially in winter.

3. The association between NDVI and climatic factors was explored by Geodetector at the regional scale. The result of Factor_detector reflected the influence degree of PRE (q = 0.54, P < 0.01) was greater than that of TEM (q = 0.27, P < 0.01) in the control of the spatial distribution of NDVI over China. The interaction influence degree q of PRE $\cap$ TEM obtained by Interaction_detector was about 0.71 in the 20 years, indicating that the response of vegetation to climate was dominated by the interaction of PRE and TEM. PRE and TEM played different roles in seven sub-regions under different geographic and geomorphic conditions. With the estimation of Risk_detector, the most suitable PRE for vegetation growth is 102.43 ∼ 157.57 mm and the suitable TEM is above 20°C.

Acknowledgments

The authors are grateful to NASA for providing the Terra MODIS NDVI dataset (http://earthdata.nasa.gov), and the ECMWF for providing the ERA5-Land dataset (https://www.ecmwf.int/en/forecasts/datasets/browse-reanalysis-datasets).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are openly available in ‘figshare’ at https://doi.org/10.6084/m9.figshare.19556506.

Additional information

Funding

This work was funded by the National Natural Science Foundation [grant number 71971002], the Anhui Provincial Natural Science Foundation [grant number 2108085QD154], the Major Science and TechnologyProject of Anhui Province [grant number 202003a06020016], the Key R&D Project of Anhui Province [grant number 202004a07020050].