Combining gradual and abrupt analysis to detect variation of vegetation greenness on the loess areas of China

Abstract

The annual peak growth and trend shift of vegetation are critical in characterizing the carbon sequestration capacity of ecosystems. As the well-known area with the fastest vegetation growth in the world, the Loess Plateau (LP) lands find an enhanced greening trend in the annual and growing-season. However, the spatiotemporal dynamics of vegetation peak growth and breakpoints characteristics on time series still needs to be explored. Here, we performed tendency analysis to characterize recent variations in annual peak vegetation growth through a satellite-derived vegetation index (NDVI_max, Maximum Normalized Difference Vegetation Index) and then applied breakpoint analysis to capture abrupt points on the time series. The results demonstrated that the vegetation peak trend had been significantly increasing, with a growth rate at 0.68×10^–2·a^–1 during 2001–2018, and most pixels (70.81%) have a positive linear greening trend over the entire LP. In addition, about 83% of the breakpoint type on the monthly NDVI time series is a monotonic increase at the pixel level, and most pixels (57%) have detected breakpoints after 2010. Our results also showed that the growth rate accelerates in the northwest and decelerates in the southeast after the breakpoint. This study indicates that combining abrupt analysis with gradual analysis can describe vegetation dynamics more effectively and comprehensively. The findings highlighted the importance of breakpoint analysis for monitor timing and shift using time series satellite data at a regional scale, which may help stakeholders to make reasonable and effective ecosystem management policies.

Keywords

vegetation greenness / gradual trend / breakpoint / BFAST algorithm / the Loess Plateau area

1 Introduction

Vegetation provides water retention and climate regulation services, progressively shapes and beautifies the landscape, and provides numerous invaluable ecosystem services which are beneficial for natural ecosystem and human society alike (Li and Pan, 2018; Wang et al., 2020). Recent studies have reported that global vegetation shows an obvious greening trend due to the improvement of hydrothermal conditions and human interference (Wang et al., 2018a; Piao et al., 2019), The greening areas are mainly concentrated in the northern hemisphere, especially in the high latitudes, China and India (Chen et al., 2019), and the central United States (Tang et al., 2015). The browning area is only concentrated in rainforests and metropolises (Pütz et al., 2014; Zhou et al., 2014; Zhong et al., 2019). Meanwhile, global arid and semi-arid regions have also experienced significant overall vegetation greening processes (Fensholt et al., 2012; Andela et al., 2013). From the spatial perspective, the areas with perceptible trend of turning ‘green’ are predominantly located in the driest areas, such as Saharan Africa and Australia (Poulter et al., 2014; Kaptué et al., 2015; Brandt et al., 2017). At present, many studies have paid attention to the spatiotemporal dynamics of vegetation greenness in the annual and growing seasons at the continental and global scale. However, there is still a lack of knowledge about local vegetation’s peak growth across arid and semi-arid areas.

The breaks for the additive seasonal and trend (BFAST) algorithm is proven to have the ability to detect the breakpoints of seasonal and trend on the satellite time series, indicating the spatial distribution of breakpoint type and breakpoint year (Watts and Laffan, 2014). In addition, traditional methods such as Mann-Kendall (MK) and Copula-Based Abrupt Variations Detection (CAVD) are usually used to determine the breakpoints of vegetation variation (Deng et al., 2018; Zhao et al., 2019). However, these methods may limit vegetation change information mining to some extent, neglecting the analysis on the vegetation dynamics at different time intervals and the clustering patterns on satellite time series. At present, BFAST has been successfully applied to study spatiotemporal breakpoints of satellite time series (He et al., 2020a). For example, the BFAST algorithm can effectively monitor vegetation cover changes in the Kurdistan province of Iran, and describe the normalized difference vegetation index (NDVI) behavior in terms of greenness, phenophase, abrupt changes for the different land covers (Gholamnia et al., 2019). BFAST has also been applied to monitor the GIMMS LAI data set of karst areas with most breakpoints during 2002–2004 and explained that the large-scale ecological restoration projects have enhanced carbon sequestration and thus mitigated climate warming (Tong et al., 2018). BFAST was also used to monitor ‘when’ and ‘where’ trend changes in China. It found that vegetation in northern China was mainly affected by land-use changes dominated by forestry projects, while in southern China were more affected by climate change (Ma et al., 2019).

The Loess Plateau is located in the central part of China, which is one of the most famous loess covered areas (Yang et al., 2019b). It is an arid and semi-arid region with forest, typical grassland and desert grassland extending from southeast to northwest. Since the 2000s, ecological restoration projects, including the Three-North Shelter Forest, the Grain for Green, and the Protection of Natural Forests, have promoted the vegetation growth in LP substantially (Cao et al., 2011). Ecological projects have changed regional processes of vegetation growth, nutrition cycle and energy exchange (Wu et al., 2019a), which has fundamentally improved the ecological environment and vegetation production (Fu et al., 2017). However, most previous studies in the LP mostly focused on the long-term gradual change of vegetation index and its response to climate; it is still unclear whether there were abrupt changes in the long term vegetation trend.

Based on the Google Earth Engine (GEE) cloud computing platform, this study preprocessed the Moderate-resolution Imaging Spectroradiometer (MODIS) NDVI data set. We analyzed the gradual trend and abrupt change of vegetation time series over the LP. The scientific questions of this study are mainly: 1) whether there is significant spatial heterogeneity in the interannual characteristics of peak greenness on the LP; 2) when, where, and what type of breakpoint occurred in vegetation time series during 2001–2018; 3) whether the BFAST algorithm can accurately and effectively detect the major breakpoints.

2 Materials and methods

2.1 Overview of the study area

The Loess Plateau is located in the center of China in spatial terms, the transition zone between semi-arid and arid areas, and it is the most concentrated and largest loess zone on the earth. The Loess Plateau comprising eastern Qinghai, southeast Gansu, Ningxia, Shaanxi, Shanxi, and part of Inner Mongolia (Fig. 1). The Loess Plateau covers an area of 640000 km², and the altitude is 200–4975 m (Wang et al., 2016). The climate can be classified as continental monsoon. The average annual precipitation is below 400 mm, with the precipitation exhibiting a pattern of gradual decrease from the south-east to the northwest part of Loess Plateau. The average annual temperature is around 9°C, and the average annual evapotranspiration ranges from 820 to 1650 mm (Wang et al., 2016; Wang et al., 2018c). The soil has high porosity and low capacity of water retention, facilitating a potential of evaporation at high rate (Zhao et al., 2017). The forests had been widely distributed in Loess Plateau throughout history. However, the deforestation process as a result of human activities during the past century has inevitably escalated the magnitude of soil erosion in this area (Guo et al., 2019).

2.2 Data introduction

MODIS NDVI has a positive relationship with vegetation coverage and biomass (He et al., 2020b), suitable as the preferred indicator for vegetation monitoring and evaluation over the LP (Liu et al., 2017). In this study, the GEE platform provided a MOD13A1 NDVI data set with a spatial resolution of 500 m.

To analyze vegetation dynamics at different time scales (month, year), the maximum value composite (MVC) was used to generate monthly and annual NDVI maximum values (Holben, 1986). The formula is as follows:

(1){\mathrm{NDVI}}_{\max }=\max \left({\mathrm{NDVI}}_i\right),

where NDVI_i refers to the NDVI value of time variable i (i = 1, 2, 3,..., 23); MOD13A1 provides 23 images a year. NDVI_max refers to the NDVI annual and monthly maximum corresponding each pixel during 2000–2018. It should be noted that the initial digital number of the MODIS NDVI used here is between –2000 and 10000 (the true NDVI needs to be multiplied by a scale factor of 0.0001).

We also collected Landsat surface reflectance images (Landsat) on the GEE platform, including TM, ETM+, and OLI image archives. The CFMask algorithm is used to detect and remove interference from clouds, shadows, snow, and ice, and then extract the greenness value of the vegetation. According to Eq. (1), we calculated the NDVI_max of MODIS and Landsat images, separately.

The MCD12Q1 land cover data set was used to generate a biomes type map on the LP (Fig. 1(b)). According to the International Geosphere-Biosphere Program (IGBP) plan, this data set was classified using a supervised decision tree classification with a spatial resolution at 500 m. The time resolution is one year. The land cover types were divided into 17 types, including 11 natural vegetation types, 3 land development and mosaic types, and 3 non-grassland types. On the basis of MCD12Q1 land classification, we obtained seven biome types over the LP.

2.3 Gradual analysis

This study analyzed NDVI_max from the provinces, and biomes and pixel scales over the LP. Trend analysis refers to the persistence of a specific indicator of the time series (such as NDVI, temperature, precipitation, etc.) over a long period of time and concerning overall increase or decrease. Linear regression is applied to calculate the inter-annual variation trend of NDVI_max under different time (year, month) and space (LP, province, and biome) scales. The slope of the linear regression equation is defined as the annual trend of the indicator. The formula for calculating slope is as follows:

(2)\mathrm{Slope}={\displaystyle \frac{n\times {{\displaystyle \Sigma }}_{i=\mathrm{1}}^n\left(i\times x_i\right)-{{\displaystyle \Sigma }}_{i=\mathrm{1}}^{n}i{{\displaystyle \Sigma }}_{i=\mathrm{1}}^n{x}_i}{n\times {{\displaystyle \Sigma }}_{i=\mathrm{1}}^n{i}^{\mathrm{2}}-\left({{\displaystyle \Sigma }}_{i=\mathrm{1}}^{n}i\right)}},

where the Slope is the linear regression equation fitted by NDVI_max, which indicates the changing trend from the time series; i is the time variable, and n is the number of years. When Slope<0, it means NDVI decreases with time; and when Slope>0, it means that NDVI increases with time. The larger the absolute value of Slope, the larger the magnitude of change of NDVI. To determine whether the linear trend of NDVI is significant or not, the t statistic was used to conduct a significance test at the significance level of 0.1, 0.05, and 0.01, respectively.

2.4 Abrupt analysis

We applied Break for the Additive Season and Trend (BFAST) analysis on per-pixel level to identify abrupt points in the vegetation time series and automatically decompose the time series into trend terms, season terms, and residual terms (Watts and Laffan, 2014). BFAST repeatedly and iteratively detects the time and number of abrupt points in time series, and describes vegetation’s characteristics according to the direction and amplitude of the breakpoints, which can be used to analyze satellite time series. Verbesselt et al. (2010a, 2020b) provided BFAST algorithm to detect jumping points in satellite time series, combined with actual events to explain the detected breakpoints. The BFAST algorithm has been used to analyze the long-term time series of the spectral index in remote sensing images, detecting seasonal changes (annual temperature, precipitation, or changes in vegetation types, etc.), revealing trend changes (annual average precipitation, land management changes, land degradation, etc.), and identifying abrupt points (deforestation, urbanization, floods, fires, etc.). It is advantageous in that there is no need for selecting a specific reference period, setting thresholds, or defining change orbits. Moreover, it can analyze long-term trend terms and seasonal terms in other scientific fields such as hydrology, climatology, econometrics, etc. (Verbesselt et al., 2010a, 2010b)

2.4.1 BFAST decomposition model

The BFAST algorithm decomposes the integrated time series into three parts: trend, seasonality, and residual. An additive decomposition model that iteratively fits a piecewise linear trend and seasonal model was used to detect and characterize the breakpoints in the trend and seasonality on time series. Assume that an additive decomposition model can iterate out and trend a piecewise linear model that matches seasonality. The algorithmic form of this model is:

(3)Y_t=T_t+S_t+e_t\mathrm{}\mathrm{}\mathrm{}t=\mathrm{1},\mathrm{2},\mathrm{3},...,n,

where Y_t is the value observed at time t; T_t is the trend component; S_t is the seasonal component; e_t is the residual component, which is a variable component independent of seasonality and trend. Assume that T_t’s change is a cent segment which is linear, and there are t₁, t₂,..., t_m.

(4)T_t={\alpha }_j+{\beta }_{j}t,

where the range of t is t _j–1<t<t _j and j= 1, 2, ..., m; {\alpha }_j and {\beta }_j is intercept and slope of the continuous linear model, respectively, which can be obtained from the size and direction of the breakpoints. The slope is the trend that gradually changes between the detection breakpoints. The magnitude of each breakpoint is the difference of T_t on t_j–1 and t_j. Seasonal trends can be expressed as:

(5)S_t={{\displaystyle \sum }}_{i=\mathrm{1}}^{s-\mathrm{1}}{\gamma }_{i,j}\left(d_{t,i}-d_{t,\mathrm{0}}\right),

where s is the seasonal period (such as the number of observations per year); {\gamma }_{i,j} is the impact factor for the ith time period. When time t is the ith observation, d_{t,i} = 1, otherwise 0, therefore, when t is the 0-th seasonal period, d_{t,i}-d_{t,\mathrm{0}} = –1. For t for other seasonal periods, d_{t,i}-d_{t,\mathrm{0}} = 1. Usually d_{t,i} is considered as a seasonal dummy variable; it has two values of 0 and 1 to account for seasonal periods in the regression model.

In the above decomposition model, parameters such as the location and number of breakpoints need to be determined. Under the interference of natural disasters, human activities, and sensor errors, it is susceptible to produce erroneous time series breakpoints. Verbesselt pointed out that the ordinary least squares residual-based moving sum (OLS-MOSUM) can tackle this problem more effectively, and Bayesian information criterion (BIC) can be a plausible approach to determine the optimal number of breakpoints (Verbesselt et al. 2010a).

2.4.2 Detecting major trend shifts in satellite time series data

To detect the NDVI major breakpoint of the time series, we applied the BFAST method by the “bfast01” function of package ‘bfast’ (available at Additive Season and Trend project website) in the R platform to detect the most important breakpoints at a per-pixel level during 2001–2018, and further determine the position and time of the breakpoint. The BFAST algorithm can detect 6 different types of breakpoints (supplement Fig. 1 according to the previous studies (de Jong et al., 2013)), namely monotonic increase, monotonic decrease, increase with a negative break, decrease with a positive break, increase to decrease, and decrease to increase. It can be intuitively understood that the monotonic increase and decrease, the breakpoints detected by the BFAST algorithm, are consistent with the trend direction before and after the breakpoint. However, the trends of increase with a negative break and decrease with a positive break are interrupted, although the trend direction is still the same, the value of NDVI interrupts the previous change trend at a certain point in time due to natural disturbances or human activities.

3 Result

3.1 Overall trends of peak greenness from MODIS

Figure 2(a) shows the linear growth trend of vegetation greenness in the LP during 2001–2018, with an average growth rate of 0.68×10^–2·a^–1. The multi-year average value of NDVI_max was 0.542, and the lowest and highest value appeared in 2001 and 2018, respectively. After 2011, the value of NDVI_max was all higher than 0.54. Figure 2(b) exhibited an increasing trend of the vegetation greenness of each province in LP. The annual greenness of SX, SHX, and GS was highly ranked with different proportion accounting for 29.68%, 22.64%, and 17.19% of the total greenness. From 2011 to 2018, the proportion of greenness decreased from 7.30% to 6.04% in QH, while the corresponding values in NMG increased from 12.41% to 14.39%. As shown in Fig. 2(c), the proportion of each biome’s greenness in total greenness remained basically unchanged. GL, CP, and FR possessed a relatively high proportion with 54.46%, 27.56%, and 8.58% respectively.

The linear trend and P value of greenness change over the LP in the spatial dimension is demonstrated (Fig. 3). Since 2001, we found that NDVI_max of 93.21% pixels over the LP displays an increasing trend, and 70.81% of pixels showed a significant increase. The average growth trend is 5.6×10^–3·a^–1, and the area with average growth rate exceeding 5×10^–3·a^–1 accounts for 62.73% of the entire region. Area with average growth rate exceeding 4·10^–3 and 3×10^–3·a^–1 accounted for 70.68% and 78.41% respectively. The pixels with a decrease rate greater than 3×10^–3·a^–1 accounted for only 2% and were concentrated in irrigation areas such as river valleys and cities under strong human disturbance. The MODIS NDVI_max of most pixels have significantly changed, and the proportions of P<0.1, P<0.05, and P<0.01 in the total pixels was 78.54%, 72.53%, and 58.84%, respectively (Fig. 3(b)).

3.2 Timing and shift in trends of MODIS NDVI time series

The BFAST method was applied on NDVI at a per-pixel level to identify the major breakpoints in satellite time series over the LP. The results showed that BFAST detected 6 major breakpoint types (Fig. 4(a)) and found that most of the breakpoints showed an increasing trend before and after each breakpoint; i.e., a monotonic increase. The monotonic increase is the dominant breakpoint type which accounts for 83.18% on the LP; and the number of pixels that have decreased accounted for only 2.29%. The proportion of increase with a negative break and decrease with a positive break consists 7.79% and 0.68%. And the pixels with reversed trends, decrease to increase and increase to decrease accounted for 2.43% and 3.60%, respectively. BFAST has disclosed that the satellite time series has abrupt changes between 2004 and 2015 (Fig. 4(b)), and the time nodes of breakpoints at a per-pixel level were moderately dispersed. About 58% of pixel breakpoints occur after 2010. The year with the largest number of breakpoints was 2011, and 2008 was the year with the smallest number of breakpoints.

The BFAST method revealed the change amplitude of the MODIS NDVI time series before and after the breakpoint (Fig. 5), and the changes were essentially consistent with the linear trend of NDVI during 2000–2018. The growth trend at pre-stage (before the breakpoint) accounted for 87.35% of all pixels over the LP and the number of pixels with a growth rate of more than 5×10^–4·mm^–1 accounts for 36.76%. The results reveal that most of the pixels showed a positive trend before the breakpoint, and only a small part of central Gansu and the Guanzhong Plain showed a negative trend. In post-stage (after the breakpoint), the increasing trend accounts for 85.59% of all the pixels, and the number of pixels increasing at a rate of more than 5×10^–4·mm^–1 accounts for 41.44%. The negative trend was randomly distributed without noticeable tendency of spatial clustering on the LP.

In addition, we calculated the difference in change amplitude before and after breakpoint (post minus pre). Figure 5(c) showed that the rate of NDVI change in south-east of the LP slowed down after the abrupt change, yet the growth rate was changing faster than before in the north-west of the LP. Moreover, we found that the linear trend can merely indicate the overall change of vegetation greenness during the study period, and the BFAST method can explicitly examine the specific breakpoint and changing intensity of vegetation in the study area.

As of NDVI change at the provincial scale, the growth rates of GS, NMG, NX, and QH before the breakpoint was higher than those after the breakpoint; while the growth rate of HN, SHX and SX decreased significantly after breakpoint (Fig. 6). At provincial scale, change amplitude can be correlated to the geographical location of the provinces. NDVI of provinces situated in the north-west part of Loess Plateau displayed an elevated growth rate after the breakpoint; meanwhile, NDVI of provinces situated in the south-east part of Loess Plateau displayed a declined growth rate after the breakpoint. For biomes, GL, which biomes type possesses the largest potential capacity for vegetation growth, had a growth rate exceeding 2×10^–4·mm^–1 after a breakpoint. However, the growth rate of FR and SS remains stable after a breakpoint, and it is not likely to maintain vegetation greening at previous rate in the future. The CP’s growth rate decreased after the breakpoint, which meant that the relative contribution of cropland to vegetation improvement decreased.

3.3 Effectiveness of the BFAST algorithm through two typical cases

Figure 7(a) showed the significance of the variation amplitude before and after the breakpoint. The significance level is defined as 0.05. It can be found that the vegetation greenness before and after breakpoint in 83.68% of pixels of the LP has changed significantly. 5.17% of pixels were significant before the breakpoint and 4.09% were significant after the breakpoint. Only 7.88% of the pixels were not significant. As illustrated in Fig. 8(a), the proportion of unstable pixels was more than 90% to the two segments (before and after the breakpoint), which implies that the majority of vegetation pixels in the LP were not sustained abidingly at a state of high resilience and can be altered with forestation effort.

The effectiveness of the BFAST method in monitoring abrupt changes of vegetation dynamics can be verified based on a known observation site. In this study, we selected the Junger coalmine area to verify the effectiveness of the BFAST algorithm. The Junger coalmine is the largest open-pit coalmine in Asia, the coal production of which accounts for about 1/6 of the total coal production in China. The coal seam is thick, buried shallower, and the geological structure is stable, suitable for open-pit mining. The pixel selected is in the south-east area of ​​the Junger coalmine, which was excavated after 2010. As showed in Fig. 7(c), the vegetation greenness was still high before 2010, but the greenness declined significantly after 2010. Meanwhile, we also downloaded the annual maximum NDVI data set of the Landsat satellite in the black frame circle (Fig. 7(b)) on the GEE cloud platform to verify the MODIS NDVI. The vegetation greenness of the sample points was relatively high in 2001; however, the open-pit coal mine area gradually expanded to the sample point in the south-east after 2010. Greenness exhibited a simultaneous decline with the mining area expanded into the sample point after 2010 (Fig. 7(d)). During this period, the vegetation changed from greening to browning.

Furthermore, the impact of human activities on the abrupt change from browning to greening was investigated. We selected the Hongsibu irrigation area of ​​the Ningxia Poverty Alleviation and Irrigation Project (NPAIP) as a research site (Fig. 8(b)). The Hongsibu, located in the middle of Ningxia, is the main battlefield of the national water conservancy project, and also is the largest ecological migration plan in China. Due to the lack of rainfall in southern Ningxia, agricultural production is low and socio-economic development is lagging. To solve the poverty in southern Ningxia, the government first raised the water level of the Yellow River to a certain elevation through pumping stations. Then the high-elevation water was self-irrigated to the Hongsibu area, to develop an artificial oasis. We found the BFAST algorithm can effectively monitor that the year node of breakpoint is 2009, and the vegetation greenness rises from 0.2 to 0.6 during 2009–2018 (Fig. 8(c)). The barren desert was originally poorly vegetated due to lack of water resources, but now it has been transformed into cropland with higher greenness. Landsat NDVI data set from 2001 to 2018 can also found that the vegetation greenness was very low before 2008. Since 2009, the greenness of the selected sample points increased notably compared with the previous (Fig. 8(d)), which is consistent with MODIS NDVI. In short, these two cases are provided to validate that the BFAST algorithm can effectively detected the breakpoints in vegetation time series.

4 Discussion

4.1 Enhance the growth of vegetation and its possible impact factors over the LP

We used a remote sensing image data set, a linear trend model, and a breakpoint method to study the spatiotemporal change in vegetation greenness. The peak growth of vegetation (NDVI_max) can accurately quantify the best growth state in a year (Yang et al., 2019a). Previous studies found that the variance in vegetation productivity is primarily restrained by the peak growth and the length of the growth season (Xia et al., 2015). Accordingly, the peak growth is critical in characterizing the capacity of terrestrial ecosystem productivity and configuration of the seasonality of atmospheric CO₂ concentration (Huang et al., 2018). This study found that the peak growth of vegetation across the LP consistently exhibited linearly increasing trends during the past two decades (2000–2018). More than 70% of the pixels had a significant increasing trend (Fig. 3). Meanwhile, we found that the fastest increase in peak vegetation greenness occurred in the central area of the LP, such as Northern Shaanxi, Southern Ningxia, Eastern Gansu, and Western Shanxi. The negative-trend pixels appeared mostly near cities and remote areas. Furthermore, we observed that the peak growth of each province and biome also had demonstrated widespread increasing trends over 2000–2018 (supplement Fig. 2). Notably, grassland and cropland had higher growth rates, which are 0.78×10^–2·a^–1 and 0.63×10^–2·a^–1, respectively. In particular, we noted that inter-monthly vegetation changes are a substantial component of enhanced peak growth. Therefore, we further analyzed the spatiotemporal dynamics of maximum vegetation greenness for every month, and found that the growth season (June to October) is the main period of vegetation restoration on the LP (supplement Figs. 3 and 4). Overall, the above results indicated that the vegetation greenness ​​had increased significantly at different spatiotemporal scales over the LP, and the vegetation has been restored remarkably and widespread greening has taken place during in the past two decades, which is consistent with previous research (Fu et al., 2017; Wang et al., 2018b; He et al., 2020b).

Breakpoints type of monotonic greening comprise the vast majority in the vegetation time series, which meant that the greenness at per-pixel level has predominantly increased both before and after the breakpoint. This is similar compatible with the result of gradual change analysis, and the catastrophe analysis also revealed that the vegetation in the LP presents a gradual recovery tendency. In particular, our research focused on the variations of vegetation before and after the breakpoint. We found that the growth rate in the southeast of the LP declined rapidly after the breakpoint, while the northwest significantly increased. This may be attributed to the fact that, after the rapid vegetation restoration in the southeast, the greenness had reached a certain saturation level (Ma et al., 2019). It could also be explained by the fact that the vegetation coverage in the north west was relatively scarce (Ning et al., 2015), the growth of pioneer species was relatively slow, and it only showed a rapid growth trend in recent years after pioneer species have established. It is worth noting that the main vegetation type in the north west of the LP is grassland. Therefore, we predict that under the scenario of continuous vegetation restoration in the north-west, grassland will continue to grow rapidly and become the biome with the most potential of greening in the LP. In conclusion, the combination of gradual analysis and abrupt analysis can comprehensively reveal that the widespread increase of LP’s vegetation greenness in the past two decades, which can be credited to vegetation restoration.

The arid and semi-arid ecosystem is fragile and complex, which is extremely vulnerable to disturbance from human activities (Li and Pan, 2018). A series of ecological restoration projects have been implemented on the LP, especially the implementation of the “Grain-to-Green Program,” featuring the administrative actions empowered by central government, the vast amount of investment and the far-reaching influence (Wu et al., 2019a; Han et al., 2020), which improved the vegetation coverage and carbon sequestration capacity to great extent (Xiao, 2014; Zhang et al., 2018; Wang et al., 2018d). The government has made overzealous legislative efforts to stop agriculture practices on steep slopes which are prone to heavy runoff, counteract desertification, and promote vegetation restoration by afforestation and grazing prohibition as per local conditions (Yang et al., 2019b). The gully area of the LP is where ecological restoration projects have the most conspicuous rehabilitation outcome (Cao et al., 2011; Lu et al., 2018; Ma et al., 2019). Vegetation growth had increased rapidly and significantly since 2001, the cumulative area of returning grain to green is in line with the greenness increasing, which indicated that the term ‘green’ might be preliminarily attributed to the impact of ecological restoration projects (Li et al., 2019; Han et al., 2020).

The vegetation greenness increased rapidly over the LP from 2001 to 2018, and all vegetation biomes showed an increasing trend, especially the greenness of cropland. In the past, the yield and aboveground biomass of most crops on LP were low, due to lack of fertilizer application and effective irrigation, which resulted in the low greenness value of cropland in the LP. However, the extraordinary progress of the Chinese economy after 2001 has enabled the government to devote more resources in agricultural production (Qiao et al., 2018). A series of policies that aim to increase crop yields have been implemented (Wu et al., 2019b). For example, the encouragement of chemical fertilizer application has contributed to increase of crop yields and fertilizer application rate. Agriculture management strategies and techniques (e.g. mulching) facilitate corn (wide leaves) to be planted, and the original low-yield crops are replaced by improved varieties with higher yield (Lu and Liao, 2017). It is worth mentioning that the construction of terraces helps reduce water runoff and soil erosion on the LP (Sun et al., 2019). The Chinese government has constructed a large number of stepped sections of farmland along the contour lines on the slopes of the LP (LaFevor, 2014), thus levelling the originally steep slopes from 5°–25°, and then modifying the ground slope and runoff coefficient, reducing soil erosion and increasing accumulation of organic carbon on the LP (Wang et al., 2018c). These measures improved the growth conditions of crops, and at the same time, promoted the increase of the cropland greenness in satellite images.

4.2 Effectiveness of the BFAST method on vegetation monitoring

Here, we used MODIS NDVI data sets from 2000 to 2018 to characterize and analyze the spatio-temporal change of vegetation by combining gradual analysis (trend change) and abrupt analysis (breakpoint monitoring). The gradual analysis is the most commonly adopted trend analysis technique, reflecting the interannual and seasonal change trend. It can evaluate the overall change direction and change amplitude of vegetation. However, gradual analysis is not capable of capturing the nonlinear and nonstationary characteristics on vegetation time series induced by climate change and human activities (Gholamnia et al., 2019). Previous studies have indicated that the BFAST algorithm can discriminate timing and shift in trends of satellite time series on vegetation time series, and amend for the defects of gradual analysis in vegetation monitoring (Tong et al., 2018; Niu et al., 2019). Therefore, the combination of gradual change analysis and abrupt change analysis can estimate vegetation greenness and carbon sink capacity more accurately.

Furthermore, we applied the pixel-based BFAST algorithm to detect the abrupt changes in vegetation dynamics on the LP. We found that BFAST detected significant space-time heterogeneity of vegetation greenness. Compared to other linear models or the MK test (Wu et al., 2019a), BFAST can distinguish ‘when’, ‘where’, and ‘what types’ of breakpoint occurs on time series at another level of accuracy (Niu et al., 2019), and help us analyze the vegetation change in key areas. Fang et al. (2018) used BFAST to detect vegetation dynamics in Quebec, Canada, and proposed that abrupt vegetation changes considerably varied with different land cover types. They also found that climate anomaly and frequent fire are the primary causes of abrupt vegetation changes (Fang et al., 2018). Geng et al. (2019) found that browning breakpoint was more widespread than greening breakpoint in Qilian Mountains using BFAST method, and suggested that long-term overgrazing might be the main cause of abrupt browning changes. Our research results showed that the pixels in most areas display greening monotonically over the LP, which meant that the growth in most areas is gradual, and the vegetation restoration was not interrupted. The monotonic decreases are only concentrated in urban areas where the rapid urbanization of China in the past ten years has exerted most pressure (Xiong et al., 2019). Meanwhile, some drastic vegetation greenness changes are spotted in some areas, such as the mining of open-pit coal mines (Fig. 7), which will interrupt the implementation of greening policy and lead to vegetation degradation. Moreover, it should also be mentioned that the relocation policy and soil improvement will promote the transformation of abandoned land into cropland (Fig. 8). As the conclusions, BFAST algorithm has the capability to monitor the major breakpoint and trend change on satellite time series. It can produce intuitive results of the changing status of vegetation, which cannot be obtained through linear and MK trend analysis.

5 Conclusions

This study used the MODIS NDVI (monthly, 500 m) remote sensing observation data set to study vegetation’s temporal and spatial dynamics of vegetation on the LP. Our results indicated that the vegetation on the LP had shown a greening trend in the past two decades. NDVI of 70.81% of the number of pixels increased significantly; and the proportion of pixels with positive changes during growing season exceeded 64%. The BFAST method has a unique advantage in satellite time series analysis, and it has high precision in fitting trends. We found that vegetation time series have abruptly changed during 2004–2015, and 2011 was the year with the largest number of breakpoint pixels. The monotonic increase was the main breakpoint type, accounting for 83.18% over the entire area. The implementation of ecological projects over the LP and cropland management measures are the main human factors that promote vegetation restoration. Generally, the LP has become one of the regions with fastest-growing vegetation cover, and can help us estimate our capacity of offsetting anthropogenic greenhouse gases by means of vegetation restoration.

References

[1] Andela N, Liu Y, van Dijk A I J M, de Jeu R A M, McVicar T R (2013). Global changes in dryland vegetation dynamics (1988–2008) assessed by satellite remote sensing: comparing a new passive microwave vegetation density record with reflective greenness data. Biogeosciences, 10(10): 6657–6676

[2] Brandt M, Rasmussen K, Peñuelas J, Tian F, Schurgers G, Verger A, Mertz O, Palmer J R B, Fensholt R (2017). Human population growth offsets climate-driven increase in woody vegetation in sub-Saharan Africa. Nat Ecol Evol, 1(4): 0081

[3] Cao S, Chen L, Shankman D, Wang C, Wang X, Zhang H (2011). Excessive reliance on afforestation in China’s arid and semi-arid regions: lessons in ecological restoration. Earth Sci Rev, 104(4): 240–245

[4] Chen C, Park T, Wang X, Piao S, Xu B, Chaturvedi R K, Fuchs R, Brovkin V, Ciais P, Fensholt R, Tømmervik H, Bala G, Zhu Z, Nemani R R, Myneni R B (2019). China and India lead in greening of the world through land-use management. Nat Sustain, 2(2): 122–129

[5] de Jong R, Verbesselt J, Zeileis A, Schaepman M(2013). Shifts in global vegetation activity trends. Remote Sens, 5: 1117–1133

[6] Deng H, Pepin N C, Liu Q, Chen Y (2018). Understanding the spatial differences in terrestrial water storage variations in the Tibetan Plateau from 2002 to 2016. Clim Change, 151(3): 379–393

[7] Fang X, Zhu Q, Ren L, Chen H, Wang K, Peng C (2018). Large-scale detection of vegetation dynamics and their potential drivers using MODIS images and BFAST: a case study in Quebec Canada. Remote Sens Environ, 206(1): 391–402

[8] Fensholt R, Langanke T, Rasmussen K, Reenberg A, Prince S D, Tucker C, Scholes R J, Le Q B, Bondeau A, Eastman R, Epstein H, Gaughan A E, Hellden U, Mbow C, Olsson L, Paruelo J, Schweitzer C, Seaquist J, Wessels K (2012). Greenness in semi-arid areas across the globe 1981–2007 — an Earth Observing Satellite based analysis of trends and drivers. Remote Sens Environ, 121(1): 144–158

[9] Fu B, Wang S, Liu Y, Liu J, Liang W, Miao C (2017). Hydrogeomorphic ecosystem responses to natural and anthropogenic changes in the Loess Plateau of China. Annu Rev Earth Planet Sci, 45(1): 223–243

[10] Geng L, Che T, Wang X, Wang H (2019). Detecting spatiotemporal changes in vegetation with the BFAST model in the Qilian Mountain Region during 2000–2017. Remote Sens, 11(2): 103

[11] Gholamnia M, Khandan R, Bonafoni S, Sadeghi A (2019). Spatiotemporal analysis of MODIS NDVI in the semif-arid region of Kurdistan (Iran). Remote Sens, 11(14): 1723

[12] Guo Y, Peng C, Zhu Q, Wang M, Wang H, Peng S, He H (2019). Modelling the impacts of climate and land use changes on soil water erosion: model applications, limitations and future challenges. J Environ Manage, 250(15): 109403

[13] Han Z, Huang S, Huang Q, Bai Q, Leng G, Wang H, Zhao J, Wei X, Zheng X (2020). Effects of vegetation restoration on groundwater drought in the Loess Plateau China. J Hydrol (Amst), 591(1): 125566

[14] He P, Sun Z, Han Z, Ma X, Zhao P, Liu Y, Ma J (2020a). Divergent trends of water storage observed via gravity satellite across distinct areas in China. Water, 12(10): 2862

[15] He Y, Yan H, Ma L, Zhang L, Qiu L, Yang S (2020b). Spatiotemporal dynamics of the vegetation in Ningxia China using MODIS imagery. Front Earth Sci., 14(1): 221–235

[16] Holben B (1986). Characteristics of maximum-value composite images from temporal AVHRR data. Int J Remote Sens, 7(11): 1417–1434

[17] Huang K, Xia J, Wang Y, Ahlström A, Chen J, Cook R B, Cui E, Fang Y, Fisher J B, Huntzinger D N, Li Z, Michalak A M, Qiao Y, Schaefer K, Schwalm C, Wang J, Wei Y, Xu X, Yan L, Bian C, Luo Y (2018). Enhanced peak growth of global vegetation and its key mechanisms. Nat Ecol Evol, 2(12): 1897–1905

[18] Kaptué A T, Prihodko L, Hanan N P (2015). On regreening and degradation in Sahelian watersheds. Proc Natl Acad Sci USA, 112(39): 12133–12138

[19] La Fevor M C (2014). Restoration of degraded agricultural terraces: rebuilding landscape structure and process. J Environ Manage, 138(1): 32–42

[20] Li G, Sun S, Han J, Yan J, Liu W, Wei Y, Lu N, Sun Y (2019). Impacts of Chinese grain for green program and climate change on vegetation in the Loess Plateau during 1982–2015. Sci Total Environ, 660(1): 177–187

[21] Li Z, Pan J (2018). Spatiotemporal changes in vegetation net primary productivity in the arid region of northwest China 2001 to 2012. Front Earth Sci, 12(1): 108–124

[22] Liu Z, Wu C, Liu Y, Wang X, Fang B, Yuan W, Ge Q (2017). Spring green-up date derived from GIMMS3g and SPOT-VGT NDVI of winter wheat cropland in the North China Plain. ISPRS J Photogramm Remote Sens, 130(1): 81–91

[23] Lu F, Hu H, Sun W, Zhu J, Liu G, Zhou W, Zhang Q, Shi P, Liu X, Wu X, Zhang L, Wei X, Dai L, Zhang K, Sun Y, Xue S, Zhang W, Xiong D, Deng L, Liu B, Zhou L, Zhang C, Zheng X, Cao J, Huang Y, He N, Zhou G, Bai Y, Xie Z, Tang Z, Wu B, Fang J, Liu G, Yu G (2018). Effects of national ecological restoration projects on carbon sequestration in China from 2001 to 2010. Proc Nat Acad Sci USA, 115(16): 4039–4044

[24] Lu X, Liao Y (2017). Effect of tillage practices on net carbon flux and economic parameters from farmland on the Loess Plateau in China. J Clean Prod, 162(20): 1617–1624

[25] Ma J, Xiao X, Miao R, Li Y, Chen B, Zhang Y, Zhao B (2019). Trends and controls of terrestrial gross primary productivity of China during 2000–2016. Environ Res Lett, 14(8): 084032

[26] Ning T, Liu W, Lin W, Song X (2015). NDVI variation and its responses to climate change on the northern Loess Plateau of China from 1998 to 2012. Adv Meteorol, 2015(1): 1–10

[27] Niu Q, Xiao X, Zhang Y, Qin Y, Dang X, Wang J, Zou Z, Doughty R B, Brandt M, Tong X, Horion S, Fensholt R, Chen C, Myneni R B, Xu W, Di G, Zhou X (2019). Ecological engineering projects increased vegetation cover production and biomass in semiarid and subhumid Northern China. Land Degrad Dev, 30(13): 1620–1631

[28] Piao S, Wang X, Park T, Chen C, Lian X, He Y, Bjerke J W, Chen A, Ciais P, Tømmervik H, Nemani R R, Myneni R B (2020). Characteristics drivers and feedbacks of global greening. Nat Rev Earth Environ, 1(1): 14–27

[29] Poulter B, Frank D, Ciais P, Myneni R B, Andela N, Bi J, Broquet G, Canadell J G, Chevallier F, Liu Y Y, Running S W, Sitch S, van der Werf G R (2014). Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle. Nature, 509(7502): 600–603

[30] Pütz S, Groeneveld J, Henle K, Knogge C, Martensen A C, Metz M, Metzger J P, Ribeiro M C, de Paula M D, Huth A (2014). Long-term carbon loss in fragmented Neotropical forests. Nat Commun, 5(5037): 5037

[31] Qiao J, Yu D, Wang Q, Liu Y (2018). Diverse effects of crop distribution and climate change on crop production in the agro-pastoral transitional zone of China. Front Earth Sci, 12(2): 408–419

[32] Sun W, Zhang Y, Mu X, Li J, Gao P, Zhao G, Dang T, Chiew F (2019). Identifying terraces in the hilly and gully regions of the Loess Plateau in China. Land Degrad Dev, 30(17): 2126–2138

[33] Tang G, Arnone J A III, Verburg P S J, Jasoni R L, Sun L (2015). Trends and climatic sensitivities of vegetation phenology in semiarid and arid ecosystems in the US Great Basin during 1982–2011. Biogeosciences, 12(23): 6985–6997

[34] Tong X, Brandt M, Yue Y, Horion S, Wang K, Keersmaecker W D, Tian F, Schurgers G, Xiao X, Luo Y, Chen C, Myneni R, Shi Z, Chen H, Fensholt R (2018). Increased vegetation growth and carbon stock in China karst via ecological engineering. Nat Sustain, 1(1): 44–50

[35] Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010a). Detecting trend and seasonal changes in satellite image time series. Remote Sens Environ, 114(1): 106–115

[36] Verbesselt J, Hyndman R, Zeileis A, Culvenor D (2010b). Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sens Environ, 114(12): 2970–2980

[37] Wang F, Haiyan Z, Dong X (2016). Quantifying changes in multiple ecosystem services during 2000–2012 on the Loess Plateau China as a result of climate variability and ecological restoration. Ecol Eng, 97(1): 258–271

[38] Wang J, Zhang D, Nan Y, Liu Z, Qi D K (2020). Spatial patterns of net primary productivity and its driving forces: a multi-scale analysis in the transnational area of the Tumen River. Front Earth Sci, 14(1): 124–139

[39] Wang Q, Zeng J, Leng S, Fan B, Tang J, Jiang C, Huang Y, Zhang Q, Qu Y, Wang W, Shui W (2018a). The effects of air temperature and precipitation on the net primary productivity in China during the early 21st century. Front Earth Sci, 12(4): 818–833

[40] Wang X, Xiao F, Feng X, Fu B, Zhou Z, Chan C (2018b). Soil conservation on the Loess Plateau and the regional effect: impact of the ‘Grain for Green’ Project. Earth Environ Sci, 109(3–4): 461–471

[41] Wang X, Wang B, Xu X, Liu T, Duan Y, Zhao Y (2018c). Spatial and temporal variations in surface soil moisture and vegetation cover in the Loess Plateau from 2000 to 2015. Ecol Indic, 95(1): 320–330

[42] Wang Y, Brandt M, Zhao M, Tong X, Xing K, Xue F, Kang M, Wang L, Jiang Y, Fensholt R (2018d). Major forest increase on the Loess Plateau China (2001–2016). Land Degrad Dev, 29(11): 4080–4091

[43] Watts L M, Laffan S W (2014). Effectiveness of the BFAST algorithm for detecting vegetation response patterns in a semi-arid region. Remote Sens Environ, 154(1): 234–245

[44] Wu Z, Wang M, Zhang H, Du Z (2019a). Vegetation and soil wind erosion dynamics of sandstorm control programs in the agro-pastoral transitional zone of northern China. Front Earth Sci, 13(2): 430–443

[45] Wu X, Wang S, Fu B, Feng X, Chen Y (2019b). Socio-ecological changes on the Loess Plateau of China after Grain to Green Program. Sci Total Environ, 678(15): 565–573

[46] Xia J, Niu S, Ciais P, Janssens I A, Chen J, Ammann C, Arain A, Blanken P D, Cescatti A, Bonal D, Buchmann N, Curtis P S, Chen S, Dong J, Flanagan L B, Frankenberg C, Georgiadis T, Gough C M, Hui D, Kiely G, Li J, Lund M, Magliulo V, Marcolla B, Merbold L, Montagnani L, Moors E J, Olesen J E, Piao S, Raschi A, Roupsard O, Suyker A E, Urbaniak M, Vaccari F P, Varlagin A, Vesala T, Wilkinson M, Weng E, Wohlfahrt G, Yan L, Luo Y (2015). Joint control of terrestrial gross primary productivity by plant phenology and physiology. Proc Nat Acad Sci USA, 112(9): 2788–2793

[47] Xiao J (2014). Satellite evidence for significant biophysical consequences of the “Grain for Green” Program on the Loess Plateau in China. J Geophys Res Biogeosci, 119(12): 2261–2275

[48] Xiong Y, Peng F, Zou B (2019). Spatiotemporal influences of land use/cover changes on the heat island effect in rapid urbanization area. Front Earth Sci, 13(3): 614–627

[49] Yang J, Dong J, Xiao X, Dai J, Wu C, Xia J, Zhao G, Zhao M, Li Z, Zhang Y, Ge Q (2019a). Divergent shifts in peak photosynthesis timing of temperate and alpine grasslands in China. Remote Sens Environ, 233(1): 111395

[50] Yang X, Na J, Tang G A, Wang T, Zhu A (2019b). Bank gully extraction from DEMs utilizing the geomorphologic features of a loess hilly area in China. Front Earth Sci, 13(1): 151–168

[51] Zhang S, Yang D, Yang Y, Piao S, Yang H, Lei H, Fu B (2018). Excessive afforestation and soil drying on China’s Loess Plateau. J Geophys Res Biogeosci, 123(3): 923–935

[52] Zhao A, Zhang A, Lu C, Wang D, Wang H, Liu H (2017). Spatiotemporal variation of vegetation coverage before and after implementation of Grain for Green Program in Loess Plateau China. Ecol Eng, 104(1): 13–22

[53] Zhao J, Huang S, Huang Q, Wang H, Leng G, Peng J, Dong H (2019). Copula-Based abrupt variations detection in the relationship of seasonal vegetation-climate in the Jing River Basin China. Remote Sens, 11(13): 1628

[54] Zhong Q, Ma J, Zhao B, Wang X, Zong J, Xiao X (2019). Assessing spatial-temporal dynamics of urban expansion vegetation greenness and photosynthesis in megacity Shanghai China during 2000–2016. Remote Sens Environ, 233(1): 111374

[55] Zhou L, Tian Y, Myneni R B, Ciais P, Saatchi S, Liu Y Y, Piao S, Chen H, Vermote E F, Song C, Hwang T (2014). Widespread decline of Congo rainforest greenness in the past decade. Nature, 509(7498): 86–90

RIGHTS & PERMISSIONS

Higher Education Press