Factors Affecting Long-Term Trends in Global NDVI

Abstract

Diagnosing the evolution trends of vegetation and its drivers is necessary for ecological conservation and restoration. However, it remains unclear what the underlying distribution pattern of these trends and its correlation with some drivers at large spatial-temporal scales. Here we use the normalized difference vegetation index (NDVI) to quantify the activity of vegetation by Theil–Sen median trend analysis and the Mann–Kendall test, Pearson correlation analysis and Boosted regression trees (BRT) model. Results show that about 34% of the global continent area has experienced greening in the grid annual NDVI from 1982 to 2015. The major greening areas were observed in the Sahel, European, India and south China. Only 10% of the global continent land areas were browning, and these were observed in Canada, South America, central Africa and Central Asia. BRT model shows that rainfall is the most important factor affecting vegetation evolution (63.1%), followed by temperature (15%), land cover change (8.6%), population (6.5%), elevation (6.4%) and nightlight (0.4%). It’s about 21% of the world’s continent were affected by rainfall, mainly in arid regions such as central Asia and Australia. The main temperature-affected areas accounted for 36%, located near the equator or in high latitudes.

1. Introduction

Vegetation is the main component of terrestrial ecosystems on Earth and is crucial in the water, energy exchange, and biogeochemical cycles of terrestrial surfaces [1,2,3]. Vegetation hugely impacts the terrestrial carbon cycles, ecosystem service value, surface energy exchange and water balance from regional to global [4,5,6,7,8,9]. The NDVI (normalized difference vegetation index) is sensitive to phenology [10] and is often used as the best indicator of vegetation growth and coverage change [11]. The NDVI has also been widely used in the study of vegetation changes at different scales [12,13]. The climate exerts dominant control over the spatial distribution of major vegetation types on a global scale [14]. Vegetation dynamics and the responses to climate changes have been recognised as a core issue of global change in terrestrial ecosystems [1].

The trends and divers of Earth’s vegetation have been in the spotlight of science [15]. Many scholars have extensively analysed the temporal and spatial evolution characteristics and driving factors of the NDVI. Zhu et al. [16] use three long-term satellite leaf area index (LAI) records and 10 global ecosystem models to investigate four key drivers of LAI trends during 1982 to 2009. Brandt et al. [17] found a major increase in growing season vegetation cover from 69% in 1999 to 81% in 2017. Piao et al. [18] use three different satellite-derived LAI datasets for detection as well as five different process-based ecosystem models for attribution. Rising atmospheric CO₂ concentration and nitrogen deposition are identified as the most likely causes of the greening trend in China. Song et al. [19] studied the global land change from 1982 to 2016 and found that tree cover increased by 2.24 million km² during that period. Jong et al. [20] found that the net global figure showed greening from 1981 to 2006.

To understand the long-term NDVI trends (during 1982 to 2015) on a large spatial scale (global) and their influence factors, this study analysed the trends of NDVI by Theil–Sen median trend analysis and the Mann–Kendall test. Climate and socioeconomic factors were included to screen the driving forces for vegetation change using Pearson correlation analysis and BRT model.

2. Materials and Methods

2.1. Data

2.1.1. Global Inventory Modeling and Mapping Studies (GIMMS) NDVI 3g Data

GIMMS NDVI data represent a vegetation product obtained by the Advanced Very High Resolution Radiometer (AVHRR) sensor mounted on the National Oceanic and Atmospheric Adm (NOAA) satellite. It is considered the best dataset available for long-term NDVI trend analysis [16]. Such data are widely used to study the large-scale changes in vegetation phenology at the global and regional levels [21,22]. The data cover the period of 1982–2015 and has a temporal resolution of 16 days and a spatial resolution of 0.072 (8 km) [13,22,23]. In this study, we applied the maximum value composite (MVC) technique to remove any biases caused by atmospheric conditions [24].

2.1.2. Data Sources of Impact Factors

The Climatic Research Unit (CRU) high-resolution gridded dataset is a globally used near-surface climate dataset with a horizontal resolution of 0.5° × 0.5° longitude. This study used the TS 4.00 dataset published by the University of East Anglia CRU (http://www.cru.uea.ac.uk/). The time scale of 1982–2015 was adopted. The global nighttime lights data obtained by the Operational Linescan System (OLS) supported by the Defense Meteorological Satellite Program (DMSP) provide a new tool for large-scale urbanization research. These data exclude the effects of bright clouds and water bodies with a spatial resolution of 1 × 1 km [25]. The Shuttle Radar Topography Mission (SRTM) obtained elevation data on a near-global scale to generate the most complete high-resolution digital topographic database of the Earth. The SRTM is an international project spearheaded by the National Geospatial-Intelligence Agency, the National Aeronautics and Space Administration (NASA) and the Italian Space Agency and the German Aerospace Center. This study uses the 1 km resolution Digital Elevation Model (DEM) data produced by SRTM (http://glcf.umd.edu/data/srtm/). The climate change initiatiue land cover project (CCI-LC) delivers consistent global LC maps at 300 m spatial resolution on an annual basis from 1992 to 2015.

2.2. Methods

2.2.1. Mann–Kendall Mutation Test

The Mann–Kendall mutation test is a method often used in meteorology and climatology [26]. The formula is as follows:

Set the original time series as x₁, x₂, …, x_n, m_i is the cumulative number that x_i is greater than x_j in the i_th sample, define the statistic as:

$$d_k={\sum }_{i=1}^k{m}_i,\left(2\le k\le n\right),$$

(1)

The mean and variance of d_k are respectively:

$$\mathrm{E}\left(d_k\right)=k\left(k-0.25\right),$$

(2)

$$\mathrm{Var}\left(d_k\right)=\frac{k\left(k-1\right)\left(2k+5\right)}{72},$$

(3)

UF statistics is obtained through standardization,

$${\mathrm{UF}}_{\left(k\right)}=\frac{d_k-E_{\left(dk\right)}}{\sqrt{Va{r}_{\left(dk\right)}}},$$

(4)

By referring this method to the inverse sequence and calculating another curve, UB, the intersection of the two curves in the confidence interval is determined as the mutation point (And the intersection point is within the critical line of y = ±1.96 (α = 0.05)).

2.2.2. Theil–Sen Median Trend Analysis and Mann–Kendall Test

Theil–Sen median trend analysis is a robust non-parametric statistical method for trend calculation [27,28,29].

$$S_{NDVI}=Median\left(\frac{NDV{I}_j-NDV{I}_i}{j-i}\right)1982\le i\le j\le 2015,$$

(5)

when S_NDVI > 0, it reflects that NDVI presents a growing trend, and vice versa.

Mann–Kendall is a non-parametric statistical test method used to determine the significance of a trend and are not affected by abnormal values [30,31]. When the Mann–Kendall test is applied to analyze the trend of the sequence, the test statistic is different from the above mentioned (2.2.1). The formulas are as follows:

Set the {NDVI_i}, i = 1982, 1983, …, 2015;

Define the Z-statistic as:

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

(6)

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

(7)

$$sgn\left(NDV{I}_j-NDV{I}_i\right)=\{\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},$$

(8)

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

(9)

where, n is the number of the data in the sequence, NDVI_i and NDVI_j represent the NDVI values in the i_th year and j_th year, respectively. Sgn is a sign function. At the given significance level α, when $\left|\mathrm{Z}\right|>u_{1-\mathsf{\alpha }/2}$, it indicates the significant change of time series data in the study at the level of α. This study judged the significance of NDVI time-series change trend at the confidence level of a = 0.05.

2.2.3. Pearson Correlation Analysis

For more than a decade, scholars have used NDVI data from different regions, time series and resolutions to study the relationship between vegetation and climate [32,33]. Determining the relationship between vegetation and its driving factors is mainly accomplished by calculating and verifying the correlation coefficients [34]. The formula is presented as follows:

$$R_{xy}=\frac{{\sum }_{i=1}^n\left[\left(x_i-\overline{X}\right)\left(y_i-\overline{Y}\right)\right]}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\overline{X}\right)}^2{\sum }_{i=1}^n{\left(y_i-\overline{Y}\right)}^2}},$$

(10)

where n is the number of samples; $\overline{X}$ and $\overline{Y}$ are the means of variables x and y, respectively; and R_xy is the correlation coefficient between variables x and y. if |R| ≤ 0.296; the correlation between NDVI and its’ divers are insignificant correlation; if 0.296 ≤ |R| ≤ 0.3494, the correlation between NDVI and its’ divers correlate positively or negatively at the p = 0.01 level. If |R| ≥ 0.3494, the correlation coefficients are taken as statistically significant at the p = 0.05 level.

2.2.4. Boosted Regression Trees Analysis (BRT)

The boosted regression tree (BRT) is a quantifiable factor contribution rate and prediction model established by J. Elith et al. [35], and is one of several techniques that aim to improve the performance of a single model by fitting many models and combining them for prediction. BRT uses two algorithms: regression trees (from the classification and regression tree) and boosting builds and combines a collection of models. It is widely used in medicine [36,37].

3. Results

3.1. Global Spatial Distribution of Vegetation

According to the characteristics of global NDVI distribution, we divide the NDVI into four categories, namely, non-vegetation coverage area (NDVI ≤ 0.2), low vegetation coverage area (0.2 < NDVI ≤ 0.5), medium vegetation coverage area (0.5 < NDVI ≤ 0.8) and high vegetation coverage area (NDVI > 0.8) [38]. We found that the medium vegetation coverage area accounted for the largest proportion, reaching 42%, accounting for about half of the global continental area. In addition, the proportion of non-/very low vegetation coverage is the smallest, accounting for only 10%, followed by low vegetation coverage, accounting for 22%, and high vegetation coverage accounting for 26% (

Table 1. Global vegetation coverage distribution of all levels. NDVI: (normalized difference vegetation index).

Degree	NDVI Value	Area (km2)	Percent
Non-vegetation coverage area	NDVI ≤ 0.2	18,361,920	10%
Low vegetation coverage area	0.2<NDVI ≤ 0.5	39,134,272	22%
Medium vegetation coverage area	0.5 < NDVI ≤ 0.8	76,844,352	42%
High vegetation coverage area	NDVI > 0.8	47,376,128	26%

).

We selected several warp and weft lines (100° W, 20° E, 100° E, 60° N, 30° N, 10° S) with a large distribution of continental area to analyze the spatial distribution characteristics of NDVI values (Figure 1). 100° W meridian line passes through the North America, NDVI values are higher in Canada (with some regions exceeding 0.8 and rising) and lower in the western United States [39]. The 20° E meridian line passes through the higher NDVI value areas such as central Europe and the Africa region, as well as the extremely low NDVI value areas of the Sahara. 100° E meridian line passes through the Asian continent, NDVI value is higher in Asia except for Central Asia. In the north temperate zone where the 60° N parallel passes, the NDVI values are relatively high. In the north temperate zone where the 30° N parallel passes, the distribution of NDVI is very diverse. In the tropics where the 10° S parallel passes, the value of NDVI is higher.

3.2. Time Evolution Trend of Global Vegetation

The average annual NDVI from 1982 to 2015 was 0.38. According to the Mann–Kendall mutation test, NDVI mutation occurred in 1989 (Figure 2a), with a significant increasing trend from 1982 to 1989 at a rate of 1.9 × 10⁻³ per year, and a slowly increasing trend from 1990 to 2015. Overall, the global average annual NDVI had been rising at a rate of 4 × 10⁻⁴ per year for 34 years. (Figure 2b)

In this study, Theil–Sen median trend analysis and the Mann–Kendall test are combined to reveal the global NDVI change trends from 1982 to 2015 and classify S_NDVI. The area where S_NDVI ≥ 0.0005 is the increasing area, the area of −0.0005 < S_NDVI < 0.0005 is the constant area, and the area of S_NDVI < −0.0005 is the decreasing area. The confidence level of the Mann–Kendall test is chosen to be 0.05, and the result is divided into significant changes (Z ≥ 1.96 or Z ≤ −1.96) and the nonsignificant (−1.96 ≤ Z ≤ 1.96). The result shows the two more widely distributed NDVI evolution trend types is “Significant increase” (34%) and “constant” (45%), followed by “Significant decrease” (10%) and “Slight increase” (14%), “Slight increase” (6%), “Slight decrease” (5%). (Figure 3,

Table 2. The statistics of NDVI trends.

SNDVI	Zvalue	Trend of NDVI	Percentage
S ≥ 0.0005	Z ≥ 1.96	Significant increasing	34%
S ≥ 0.0005	−1.96 < Z < 1.96	Slight increasing	6%
−0.0005 < S < 0.0005	−1.96 < Z < 1.96	Constant	45%
S ≤ −0.0005	−1.96 < Z < 1.96	Slight decreasing	5%
S ≤ −0.005	Z ≤ −1.96	Significant decreasing	10%

)

Analyses of the changes in observed GIMMS NDVI 3g data show a widespread greening trend (34% of the continent lands) (Table 2), these are mainly observed in the Sahel (Figure 4a), European (Figure 4b), India (Figure 4c) and south China (Figure 4d) [6]. The data also show decreasing NDVI trend (browning) (10% of the continent lands), and these are observed near Canada (Figure 4e), South America (Figure 4f), central Africa (Figure 4g) and Central Asia (Figure 4h).

3.3. Correlation Analysis between NDVI and Factors

3.3.1. Temperature

About 36% of the vegetation has experienced a positive and negative correlation with temperature, and among which about 27% is statistically significant. We performed correlation analysis at the temperature and NDVI, and they showed significant negative correlation in the relatively dry areas, such as North America (especially in the southern part of the US), southern Argentina, Sahara Desert, southern Africa, the Middle East, the Arabian Peninsula, northwestern China, northern Russia and Australia; which indicates that the increase in temperature is not conducive to the growth of vegetation in the region. The negative correlation may be caused by an increase in the evaporation of soil moisture resulting from the rise in temperature, which then leads to the decline of soil moisture; in this case, the ‘physiological drought’ of plants occurs and limits the photosynthesis and growth rate of plants [40]. On the equator or low latitude areas, such as the Sino–American isthmus, near the Amazonian estuary, Coast of the gulf of guinea, India and Yarlung Zangbo Grand Canyon [41], the distribution of the NDVI is significant positively correlated with temperature. These indicate that the growth of vegetation in these areas is greatly affected by temperature and increases with the increase of temperature. (Figure 5a)

3.3.2. Rainfall

About 21% of the vegetation has experienced a positive and negative correlation with rainfall, and among which, about 13% is statistically significant. Figure 5b shows that the NDVI showed significant positive correlation with rainfall in the Midwest and Southwest regions of the US, Central Argentina, the northern and southern margins of the Sahara Desert, southern Africa, the Persian Gulf coast, western India, Mongolia and central and western Australia. These areas exhibit a sparse vegetation coverage, which indicates that water is the most important factor that limits the growth of vegetation in the area. The significant negative correlation between the NDVI and rainfall is mainly distributed in the Amazon plain, western Brazilian Plateau, southern Sahara, part of Russia and southern China. These places are usually moist; thus, the growth of vegetation benefits from relatively abundant rainfall. Such excessive rainfall leads to increased cloud cover, reduced incident radiation and increased soil moisture. These conditions lead to the increase of the ground latent heat of evaporation, which in turn results in decreased temperature, reduced photosynthesis ability and unfavourable influence on vegetation growth.

3.3.3. Nighttime Lights Index

Other studies have shown a negative correlation between nighttime lights index and the NDVI [25,42]. In our study, the NDVI and nighttime lights index are significantly negatively correlated in Eastern Europe, the Middle East, western part of North America and southern Africa. However, in north Europe, southeast America, the east China coast, India, Japan and other densely populated areas, a significant positive correlation is demonstrated. The GIMMS NDVI and nighttime lights index data are also low-resolution data whose unit pixels occupy areas that are so large that we cannot show their correlation in an in-depth manner.

From a low-resolution perspective and from an overall perspective, the regions mentioned above, such as northern Europe, show high vegetation coverage (Figure 6), dense population and high economic development [43]; these regions also include eastern China [44,45,46]. In the Middle East, the dry climate is not conducive to vegetation growth, hence the NDVI value of less than 0.2 per year. A high nighttime lights index exists in developed cities where the NDVI and nighttime lights index show a negative correlation. Moreover, the enlarged view shows that the nighttime lights index and NDVI share a significantly negative correlation, with the correlation coefficient being less than −0.3 in the areas of London, Ganges River in India and other large cities such as the urban agglomerations along the Yangtze River.

3.3.4. Digital Elevation Model (DEM)

In this study, vegetation coverage and DEM were divided into 4 grades and 9 grades respectively to analyze vegetation distribution regulation at altitude (Figure 7). The results show that, above 100 m, vegetation coverage increases with elevation, and the area of medium vegetation coverage does not change much at each altitude level, which reveal that the growth of vegetation is affected by altitude to a certain extent.

3.4. Contribution of Impact Factors

We use the Generalized Boosted Regression Models (GBM) package written by Elith in R language. Set the tree complexity = 5, learning rate = 0.005, and bag fraction = 0.5. In each calculation 50% of the data was extracted, and 5 cross validations was performed. A vegetation influencing factor model consisting of 3450 decision trees was constructed. Figure 8 shows the contribution rate of an independent variable to the influence of the dependent variable, as well as the relationship between the independent variable and the dependent variable when other independent variables remain unchanged or the mean value is taken. The results show that rainfall is the most important factor affecting vegetation evolution (63.1%), followed by temperature (15%), LUCC (8.6%), population (6.5%), DEM/elevation (6.4%) and nightlight (0.4%) (Figure 8).

Figure 9 indicate the fitted values in relation to the six predictors used in the model. Areas with high vegetation coverage are concentrated in areas with temperatures greater than 0 °C, rainfall less than 3000 mm, and elevations less than 2000 m. Nightlight and population have less impact on NDVI.

4. Discussion

4.1. The Impact of Land-Use Management

In our study, one typical country was selected from each continent to analyze the evolution trend of vegetation in the country, including the United States, Brazil, Germany, Sudan, China and Australia. We found that NDVI value of the countries on the upward trend were Brazil, China and Germany, while the countries on the downward trend were the United States, Sudan and Australia (Figure 10).

In the United States, the value of NDVI was on the upward trend before 2005, but on the downward trend after 2005, which was consistent with relevant studies [47]; the value of NDVI was showing an upward trend from 1982 to 2005 in North America, due to the severe drought and subsequent recovery, irrigated agriculture expansion, record and subsequent regeneration after insect outbreaks, and forest fire subsequent regeneration.

For Germany, we can see that NDVI showed an upward trend before 2000; in 1995 to 2010, value of NDVI showing a large “W” type fluctuation trend; and showing a small fluctuation trend after 2010. In general, NDVI showed an upward trend in Germany from 1982 to 2015. In the study of Julien [48], we can see that NDVI difference images in Europe are progressively more saturated from 1983 to 1999, and we can see in the picture in the artic that the NDVI of Germany is on the increase condition too.

There’s the world’s largest tropical rainforest in Brazil, the Amazon rainforest, which rainforest is characterized by abundant vegetation, high plant height and high density, so the NDVI value is always higher than other areas. NDVI was showing a trend of increasing slowly from 1982 to 2015, but there was a sharp decline in the vicinity of 2000, which was directly related to deforestation and land reclamation [49,50,51].

Australia is the country with the most significant decline trend among the six countries. Australia’s arid and semi-arid areas reach 75%, mainly distributed in the western and central regions. The desert area is 2.69 million square kilometers, accounting for 35% of the country’s land area. In addition, 5700 thousand hectares of Australia are affected by salinization, mainly in the southeast and southwest corners. Due to the subtropical high pressure and the southeast trade wind, it is easy to accumulate water and cause salinization as Australia has less rainfall and the terrain is flat. In addition, the soil is rich in material and easily degraded [52]. The phenomenon of excessive grazing, mining, and improper introduction of creatures leading to biological invasion [53] also led to serious desertification and low vegetation coverage in Australia [54,55]. A wide range of pests such as aphids are also part of the cause of reduced vegetation coverage [56].

In China, NDVI has experienced a slow upward trend in the past 34 years, especially in recent years, China’s vegetation coverage has experienced the strongest human activities, such as the “Returning Farmland to Forests and Grassland Project”, the “Three-north Forest Protection Project”. Liu [38] analyzed the cumulative afforestation area and NDVI changes in China’s returning farmland to forest project over the years, and found that since the implementation of the project of returning farmland to forests and grassland in 1999, China’s afforestation area has continued to rise, and NDVI has also shown an upward trend. It shows that the implementation of ecological restoration projects in recent years has greatly promoted the increase of vegetation cover in China and achieved good results. This is consistent with related research, Ke Liu and Xianfeng Liu studied the trend of NDVI in China from 1982 to 2012, Li Fei studied the trend of NDVI in China from 1982 to 2006, both indicating that NDVI in China has shown an upward trend over the years [14,57]. Here, we believe that land-use management is one of the important factors affecting vegetation evolution.

4.2. Comparison with Other Studies

The conclusions of this study are consistent with other studies. We use Theil-Sen median trend analysis to analyze the vegetation trends in the three stages of 1982–1989, 1990–2000 and 2000–2015, all showing greening (Figure 11), which is consistent with the conclusion of Chen et al. (greening in 2000–2017) [6] and Zhu et al. (greening in 1982–2009) [16], which proves that the trend analysis method used in this paper is reasonable.

Meanwhile, this study is different from other studies. Firstly, the periods examined by other studies are relatively short, while this study revealed the vegetation change trend from 1982 to 2015 (all time periods of GIMMS NDVI data), making a new contribution to the vegetation trend analysis of long time series. Secondly, some studies used LAI [6,16,18,58], Fractional Vegetation Cover (FVC) [59], Enhanced Vegetation Index (EVI) [60], and even growth season data [16] to characterize vegetation. In this study, annul of NDVI was used to reveal the evolution trend of vegetation. Finally, this study used the BRT model to analyze the contribution rate of vegetation factor for the first time, and although it has been widely used in the field of medicine [36,37] it is rarely used in vegetation research.

Limitations also exist in this study, such as the lack of analysis on the impact of climate change, carbon dioxide regulation and nitrogen deposition on vegetation [6,16,61]. It is necessary for the future study to considerate more impact factors, and factors with high contribution rate could be screened for the prediction of vegetation change trends, which is more important for the development of vegetation restoration.

5. Conclusions

We analysed the trends of the NDVI and its divers based on the global GIMMS NDVI 3g data from 1982 to 2015. We also investigated the correlation between NDVI distribution and temperature, rainfall, nighttime lights index, DEM, LUCC and population. We clarified the driving factors of the changes and development of the NDVI and then derived the following conclusions:

• Over the past 34 years, the global annual NDVI has increased by 4 × 10⁻⁴/year. About 34% of the global continental area has experienced greening in the grid annual NDVI, the greening areas were observed in the Sahel, European, India and south China. Only 10% of global vegetation land was decreasing, and this was observed in Canada, South America, central Africa and Central Asia.
• BRT model shows that rainfall is the most important factor affecting vegetation evolution (63.1%), followed by temperature (15%), LUCC (8.6%), population (6.5%), DEM (6.4%) and nightlight (0.4%) (Figure 8). About 21% of the world’s continental area were affected by rainfall, mainly in arid regions such as central Asia and Australia. The main temperature-affected areas accounted for 36%, located near the equator or in high latitudes.
• The use of the nighttime lights index to characterise human activities shows that a positive correlation exists between the nighttime lights index and the NDVI in areas with high population density, such as the Ganges and Yangtze River Basins. The use of low-resolution data analysis offers low feasibility, and the improvement of data resolution should help to accurately analyse correlations.