Detecting Spatiotemporal Changes in Vegetation with the BFAST Model in the Qilian Mountain Region during 2000–2017

: The Qilian Mountain ecosystems play an irreplaceable role in maintaining ecological security in western China. Vegetation, as an important part of the ecosystem, has undergone considerable changes in recent decades in this area, but few studies have focused on the process of vegetation change. A long normalized difference vegetation index (NDVI) time series dataset based on remote sensing is an effective tool to investigate large-scale vegetation change dynamics. The MODerate resolution Imaging Spectroradiometer (MODIS) NDVI dataset has provided very detailed regional to global information on the state of vegetation since 2000. The aim of this study was to explore the spatial-temporal characteristics of abrupt vegetation changes and detect their potential drivers in the Qilian Mountain area using MODIS NDVI data with 1 km resolution from 2000 to 2017. The Breaks for Additive Season and Trend (BFAST) algorithm was adopted to detect vegetation breakpoint change times and magnitudes from satellite observations. Our results indicated that approximately 80.1% of vegetation areas experienced at least one abrupt change from 2000 to 2017, and most of these areas were distributed in the southern and northern parts of the study area, especially the area surrounding Qinghai Lake. The abrupt browning changes were much more widespread than the abrupt greening changes for most years of the study period. Environmental factors and anthropogenic activities mainly drove the abrupt vegetation changes. Long-term overgrazing is likely the main cause of the abrupt browning changes. In addition, our results indicate that national ecological protection policies have achieved positive effects in the study area.


Introduction
Vegetation, as a medium of material cycling, the water cycle, information transfer, and energy flow, is the most critical component of terrestrial ecosystems [1][2][3][4][5][6].Within the context of global climate change, it is important to explore vegetation dynamics spatially and temporally at global and regional scales [7].Through this information, global change scientists, researchers, natural resource managers and policy makers can provide more accurate evaluations and forecasts to inform decisions related to vegetation.Satellite remote sensing technology is a unique and useful method for monitoring vegetation dynamics and environmental changes in a repeatable manner due to its high spatial coverage and long temporal series [8,9].The normalized difference vegetation index (NDVI) is one of the most important vegetation indices, and has been widely used to analyze vegetation growth history, monitor current growth conditions, and predict future vegetation dynamics [10][11][12][13][14].
From the perspective of monitoring and mapping vegetation changes using long-term remote sensing data, interannual vegetation change dynamics can be summarized into three categories: (1) continual gradual increasing or decreasing change trends; (2) abrupt changes; and (3) no obvious change trends.The first and third situations, which include seasonal and interannual changes, are usually caused by human cultivation or a response to long-term climate variations.The second situation is relatively complicated and can be caused by human activities and natural events, such as fire, agricultural expansion, deforestation, urbanization, and extreme weather events.In response to global warming, many studies have investigated vegetation change dynamics using long time series of vegetation indices.Whereas previous studies mainly focused on the long-term trends of vegetation and climate change and the relationship between them [15][16][17][18][19][20], few studies have aimed to detect both the seasonal patterns and long-term changes in vegetation [21,22].To better analyze the driving factors underlying vegetation changes, it is important to detect where and when the vegetation changed, which is hereafter referred to as an abrupt change.
Several methods have been proposed to detect abrupt changes in vegetation, such as Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr) [23], the Breaks for Additive Season and Trend (BFAST) method [24], the vegetation change tracker (VCT) [25], and Detecting Breakpoints and Estimating Segments in Trend (DBEST) method [26].Among these methods, BFAST has been successfully used in several areas recently because it can automatically analyze each pixel individually in time series data without setting a threshold.For example, BFAST has been used for deforestation monitoring [27], the establishment of where and when a vegetation disturbance occurred [28], the detection of spatiotemporal change patterns in lakes [29], forest disturbance and regrowth monitoring [28], and urban expansion analyses [30].BFAST is an effective method for detecting multiple changes in time series because it integrates the decomposition of time series into trends, seasonal parameters, and remainder components.BFAST can be used to estimate the numbers and dates of changes and analyze the magnitude of a change.Watts and Laffan [22] indicated that BFAST is an effective algorithm to detect abrupt change trends using the number of vegetation pixels changed by known floods and fires.
Qilian Mountain ecosystems play an important and irreplaceable role in maintaining ecological security in western China and are important ecological barriers in northwestern China [31][32][33].The Qilian Mountain region is an important water conservation forest area and a snow and ice water resource area in northwest China [34,35].This area is also an important forest grassland ecosystem and wildlife protection area in northwest China.This region maintains the ecological balance and economic and social development of the Hexi Corridor oasis, slows and stops the confluence and forward movement of the Kunmingtag, Badanjilin and Tengger deserts, and establishes a system to contain sandstorms in northern China [36].Over the last few decades, the Qilian Mountains have experienced increasing rates of changes in climate and vegetation.The temperature and precipitation extremes in the Qilian Mountains have exhibited significant increasing trends, especially at high altitudes [37].Over the past 50 years, the average temperature has risen by 0.26 • C/decade in this region, which is higher than the national rate of 0.14 • C/decade [38].Several studies have focused on vegetation changes in the Qilian Mountains in recent years.Chen et al. [39] employed the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI to explore the characteristics of vegetation cover from 1982 to 2006 and showed that vegetation coverage has exhibited an increasing trend over the last 20 years.Wu et al. [40] obtained similar results for the vegetation changes in this area over the last decade using the MODIS NDVI from 2000 to 2012.Wang et al. [41] indicated that changes in vegetation over time exhibited an obvious elevational difference.However, Zhang et al. [42] used the Landsat time series NDVI data from 1986 to 2015 to investigate the land use change over the last 30 years, and the results showed that the areas of cropland, forest and grassland declined and the area of Gobi increased, while the extents of grassland fragmentation and desertification increased during the investigated decades.Therefore, there is no consistent conclusion on how the vegetation has changed in the Qilian Mountains in recent decades.In addition, previous studies focus on long-term trends of vegetation change and its response to climate, while limited studies have analyzed the time of change, the different trends exhibited by different vegetation types, and the possible driving factors.Therefore, the present study was necessary to better understand the changes in vegetation dynamics for different vegetation types and potential drivers and to assist policy makers in developing effective environmental protection and management policies.
The main objective of this study was to investigate the temporal and spatial change characteristics of the vegetation of the Qilian Mountains and their likely driving factors using the BFAST algorithm based on the newest MODIS NDVI dataset.Specifically, we aimed: (1) to detect the multiple abrupt changes in the seasonal parameters and trends for the study period; (2) to further explore the gradual and abrupt changes in different vegetation types by determining the time, magnitude, and direction of the change since 2000; and (3) to estimate the potential natural climate and anthropogenic drivers of variations in the NDVI.

Study Area
The Qilian Mountains, which are located in the center of the Eurasian continent (93 • 25 -103 • 50 E, 35 • 52 -39 • 52 N), are the largest marginal mountain range on the Qinghai-Tibet Plateau, with an elevation of 3500-5000 m.The north slope is the headwaters of three inland rivers in China, the Heihe, Shulehe and Shiyanghe Rivers.The south slope is an important water supply area for the Yellow River and Qinghai Lake.The Qilian Mountains belong to a typical continental climate and plateau climate.The natural conditions in this area are complex, with large differences in precipitation and temperature conditions, which vary greatly with changing topography.The annual average temperature is −0.8 • C, the precipitation decreases from east to west, and the annual average precipitation in the mountainous area is 150-800 mm [36].The main vegetation includes alpine meadow, grassland, desert steppe, brush, and alpine vegetation (Figure 1).trends of vegetation change and its response to climate, while limited studies have analyzed the time of change, the different trends exhibited by different vegetation types, and the possible driving factors.Therefore, the present study was necessary to better understand the changes in vegetation dynamics for different vegetation types and potential drivers and to assist policy makers in developing effective environmental protection and management policies.
The main objective of this study was to investigate the temporal and spatial change characteristics of the vegetation of the Qilian Mountains and their likely driving factors using the BFAST algorithm based on the newest MODIS NDVI dataset.Specifically, we aimed: (1) to detect the multiple abrupt changes in the seasonal parameters and trends for the study period; (2) to further explore the gradual and abrupt changes in different vegetation types by determining the time, magnitude, and direction of the change since 2000; and (3) to estimate the potential natural climate and anthropogenic drivers of variations in the NDVI.

Study Area
The Qilian Mountains, which are located in the center of the Eurasian continent (93°25′-103°50′E, 35°52′-39°52′N), are the largest marginal mountain range on the Qinghai-Tibet Plateau, with an elevation of 3500-5000 m.The north slope is the headwaters of three inland rivers in China, the Heihe, Shulehe and Shiyanghe Rivers.The south slope is an important water supply area for the Yellow River and Qinghai Lake.The Qilian Mountains belong to a typical continental climate and plateau climate.The natural conditions in this area are complex, with large differences in precipitation and temperature conditions, which vary greatly with changing topography.The annual average temperature is −0.8 °C, the precipitation decreases from east to west, and the annual average precipitation in the mountainous area is 150-800 mm [36].The main vegetation includes alpine meadow, grassland, desert steppe, brush, and alpine vegetation (Figure 1).

Datasets
The datasets used in this study include NDVI data, vegetation data, and Landsat images.The MODIS Terra time series NDVI products (MOD13A1 V006) with a 1 km spatial resolution and onemonth time step for 2000-2017 were used in this study.The data were obtained from the NASA Land Processes Distributed Active Archive Centre (LP DAAC) (https://lpdaac.usgs.gov/).The MODIS Reprojection Tool (MRT), which can be obtained from LP DAAC, was used to extract the desired bands and perform format and projection conversions.The raw HDF format images were transformed into GeoTIFF format images, and the raw sinusoidal projection was converted to WGS84/Albers Equal Area Conic projection.
The vegetation atlas of China in 2001 [43] at a spatial resolution of 1 km was used to identify the main vegetation types in the study area.The atlas is mainly derived from long-term field surveys by Chinese ecologists, and it describes the geographical distribution of 796 groups and subgroups of 54 vegetation types in 11 vegetation groups.There were eight vegetation types in our study area according to the atlas.The main natural vegetation types, which were grassland, shrub, meadow, forest, and desert steppe, were chosen for the analysis using BFAST because of the potential difficulties in using BFAST to monitor plantations [44].

Methodologies
An overview of the methods for this research is shown in Figure 2. The four major steps were as follows: (1) high-quality NDVI time series data preparation; (2) breakpoint detection using the BFAST method; (3) temporal and spatial analysis of the break point; and (4) validation of the BFAST method and discussion of the potential factors driving vegetation change.The datasets used in this study include NDVI data, vegetation data, and Landsat images.The MODIS Terra time series NDVI products (MOD13A1 V006) with a 1 km spatial resolution and onemonth time step for 2000-2017 were used in this study.The data were obtained from the NASA Land Processes Distributed Active Archive Centre (LP DAAC) (https://lpdaac.usgs.gov/).The MODIS Reprojection Tool (MRT), which can be obtained from LP DAAC, was used to extract the desired bands and perform format and projection conversions.The raw HDF format images were transformed into GeoTIFF format images, and the raw sinusoidal projection was converted to WGS84/Albers Equal Area Conic projection.
The vegetation atlas of China in 2001 [43] at a spatial resolution of 1 km was used to identify the main vegetation types in the study area.The atlas is mainly derived from long-term field surveys by Chinese ecologists, and it describes the geographical distribution of 796 groups and subgroups of 54 vegetation types in 11 vegetation groups.There were eight vegetation types in our study area according to the atlas.The main natural vegetation types, which were grassland, shrub, meadow, forest, and desert steppe, were chosen for the analysis using BFAST because of the potential difficulties in using BFAST to monitor plantations [44].

Methodologies
An overview of the methods for this research is shown in Figure 2. The four major steps were as follows: (1) high-quality NDVI time series data preparation; (2) breakpoint detection using the BFAST method; (3) temporal and spatial analysis of the break point; and (4) validation of the BFAST method and discussion of the potential factors driving vegetation change.

Preparation of High-Quality NDVI Datasets
Although the monthly maximum value composite (MVC) method has been used to decrease cloud and other atmospheric effects in the original NDVI data [45], residual noise resulting from poor atmospheric conditions, cloud cover, and unfavorable sun-sensor-surface viewing geometries remain [46][47][48].Therefore, to reduce the persistent noise and obtain high-quality NDVI time series data, the corresponding MODIS quality assurance products were used to eliminate obvious noise error, and subsequently, an effective reconstruction algorithm named the Savitzky-Golay filter was

Preparation of High-Quality NDVI Datasets
Although the monthly maximum value composite (MVC) method has been used to decrease cloud and other atmospheric effects in the original NDVI data [45], residual noise resulting from poor atmospheric conditions, cloud cover, and unfavorable sun-sensor-surface viewing geometries remain [46][47][48].Therefore, to reduce the persistent noise and obtain high-quality NDVI time series data, the corresponding MODIS quality assurance products were used to eliminate obvious noise error, and subsequently, an effective reconstruction algorithm named the Savitzky-Golay filter was used [49].The improved performance of the filtering technique on the MODIS NDVI was confirmed by Geng et al. [46], who indicated that the Savitzky-Golay filter performed best in a comparison of eight noise-reduction techniques.To eliminate uncertainty, pixels with an annual average maximum NDVI below 0.1 were considered non-vegetated regions and masked out [50].

NDVI Dynamics Changes Based on BFAST
For this study, BFAST was used to identify abrupt changes in the MODIS NDVI time series data from 2000 to 2017.BFAST is an effective algorithm to integrate the iterative decomposition of time series into seasonal, trend and remainder components, and using this method we can detect changes within a long time series [21,24].The general model can be expressed by the following equation: where Y t is the observed data at time t, Tt is the trend component, St is the seasonal component, et is the remainder component, and n is the total number of observed data.The remainder component is the remaining variation in the data beyond the seasonal and trend components [51].
It is assumed that the trend component (Tt) is a piecewise linear function with breakpoints t * 1 , . . ., t * m and defined t * 0 = 0; therefore, Tt can be expressed by the following equation [24]: where j = 1, . . ., m and t * j−1 < t ≤ t * j .α j and β j t are the intercept and slope of consecutive linear models, respectively.These variables can be used to derive the magnitude and direction of the abrupt change by calculating the difference between Tt at t * j−1 and t * j .The equation is as follows: The seasonal component (St) is also a piecewise linear seasonal model based on seasonal dummy variables [24].These variables are defined as τ * 0 = 0 and τ * m+1 = n.St represents the piecewise phenological cycle on different p + 1 (p ≥ 0) segments divided by the seasonal breakpoints, τ # 1 , . . ., τ # p (τ # 0 = 0 and τ # p+1 = n) [30], and the equation is as follows: where γ j,k and θ j,k denote the coefficients, k is the number of harmonic terms, and f is the frequency.We set k = 3 and f = 12 for the annual observation of one month of time series data for this study according to Tsutsumida et al. [30].
The iterative algorithm to detect breakpoints begins with an estimate of the seasonal component using the seasonal-trend decomposition procedure.Then, the following four steps are implemented until the number and the position of the breakpoints exhibit no changes: (1) the ordinary least squares (OLS) residuals-based moving sum (MOSUM) test is used to assess the existence of abrupt changes [52]; (2) the trend coefficients α j and β j are estimated using a robust regression based on the M-estimation; (3) the number of breakpoints is determined by minimizing the Bayesian information criterion (BIC) from the seasonally adjusted data, and the date and confidence interval of the date are estimated for each breakpoint; and (4) the seasonal coefficients γ j,k and θ j,k are estimated based on the M-estimation.
In this study, R statistical software was used to perform the BFAST analysis by adding the BFAST package.All parameters described above were determined automatically.Only one parameter, h, needed to be set, which determines the minimal segment size between potentially detected breaks.The h parameter value will affect the accuracy of the BFAST method [22,42].A high h parameter may lead to the omission of certain abrupt changes, and a low value may result in the detection of abrupt changes that do not represent the actual abrupt changes in vegetation [22].For this study, we hypothesized that the minimum interval between adjacent abrupt changes was approximately two years; therefore, h = 1/7 was used, which was also recommended by Fang et al. [44].
We also chose the "harmonic" seasonal model when executing the BFAST algorithm because it is suitable for natural vegetation while the "dummy" model is often used for croplands according to Verbesselt et al. [24].Figure 3  We also chose the "harmonic" seasonal model when executing the BFAST algorithm because it is suitable for natural vegetation while the "dummy" model is often used for croplands according to Verbesselt et al. [24].Figure 3

Linear Regression Analysis
To analyze the trends in vegetation changes during the study period, a simple linear regression was employed using the yearly maximum NDVI.The slope of the trend line in the multiyear regression equation for each pixel represents the interannual variation rate, and the function is shown below [53]: where n is the accumulative number of years in the study period, and NDVI is the maximum NDVI of the ith year.Significant interannual variability of the NDVI can be determined according to the correlation of annual time series sequences and the maximum NDVI.Specifically, a positive slope indicates an increasing trend, while a negative slope indicates a decreasing trend.

Number of Abrupt Changes
For this study area, the number of abrupt changes detected in the trend component for each pixel varied from zero to four.Approximately 80.1% of the vegetation areas experienced at least one abrupt

Linear Regression Analysis
To analyze the trends in vegetation changes during the study period, a simple linear regression was employed using the yearly maximum NDVI.The slope of the trend line in the multiyear regression equation for each pixel represents the interannual variation rate, and the function is shown below [53]: where n is the accumulative number of years in the study period, and NDVI i is the maximum NDVI of the ith year.Significant interannual variability of the NDVI can be determined according to the correlation of annual time series sequences and the maximum NDVI.Specifically, a positive slope indicates an increasing trend, while a negative slope indicates a decreasing trend.

Number of Abrupt Changes
For this study area, the number of abrupt changes detected in the trend component for each pixel varied from zero to four.Approximately 80.1% of the vegetation areas experienced at least one abrupt change from 2000 to 2017 (Table 1).In the study region, 20.7% of the area had one breakpoint, 26.6% had two breakpoints, 22.5% had three breakpoints, and 10.2% had four breakpoints.The areas of the first and fourth abrupt changes throughout the vegetation area were smaller than those of the second and third abrupt changes.Differences in the abrupt change percentages for each breakpoint number were observed among different vegetation types (Table 1).A large number of grassland, meadow, shrub, and desert steppe pixels showed one or more abrupt change.The abrupt change areas were greater than 70% for each of those four vegetation types (grassland, 83.9%; meadow, 75.0%; shrub, 76.4%; and desert steppe, 93.1%) for the entire study period, and desert steppe exhibited the largest changes area.In comparison, relatively smaller breakpoint pixels were observed for alpine vegetation (8050) and forest (1925).Alpine vegetation showed the smallest abrupt change areas (75.0%) compared with the other five vegetation types.For each abrupt change, the second abrupt change area was the largest, followed by the third, first and fourth (i.e., grassland, alpine vegetation, and desert steppe).For the meadow, shrub, and forest areas, the first breakpoint area was larger than the third.The differences in the abrupt change areas for the second and third breakpoint areas were small, whereas the differences between the fourth breakpoint and the other three breakpoints were large.The spatial distribution of the estimated number of breakpoints in the trend component is shown in Figure 4. Most of the regions with breakpoints experienced two or three abrupt changes.These regions are mainly distributed in the eastern, middle, and northern parts of the study area.A small part of the areas with four breakpoints were in the southeastern and western areas of the study region, with those in the southeast mainly surrounding the Qinghai Lake area and those in the west mainly concentrated in the desert steppe regions.The percentages of area exhibiting abrupt changes for each year in different vegetation types are listed in Table 2.For all evaluated vegetation types, significant abrupt changes were not observed in 2011, and the change in area was less than 7% during this year (Table 2).The years with large abrupt change areas were 2015 (24.5%) and 2010 (21.

Timing of the Abrupt Changes
The spatial distribution of the estimated timing of abrupt changes for each breakpoint time in the trend component is identified.The main changes are primarily distributed in the northern part of the study region, especially for the desert steppe vegetation of this area (Figure 5).However, the main change years are different for each abrupt change.The first abrupt change mainly occurred in 2002, 2003, and 2004 (Figure 5a).The second abrupt change mainly occurred in 2008, 2010, and 2015, and the changes were scattered throughout the study region (Figure 5b).The third and fourth abrupt changes were concentrated in 2014 and 2015, and these changes were also scattered throughout the study area (Figure 5c,d).The percentages of area exhibiting abrupt changes for each year in different vegetation types are listed in Table 2.For all evaluated vegetation types, significant abrupt changes were not observed in 2011, and the change in area was less than 7% during this year (Table 2).The years with large abrupt change areas were 2015 (24.5%) and 2010 (21.0%), followed by 2014 (16.0%), 2013 (15.2%), and 2012 The percentages of area exhibiting abrupt changes for each year in different vegetation types are listed in Table 2.For all evaluated vegetation types, significant abrupt changes were not observed in 2011, and the change in area was less than 7% during this year (Table 2).The years with large abrupt change areas were 2015 (24.5%) and 2010 (21.0%), followed by 2014 (16.0%), 2013 (15.2%), and 2012 (15.1%).For each vegetation type, the largest abrupt change areas occurred in 2015 and 2010 (except for desert steppe, which occurred in 2013 and 2015), and the abrupt change areas covered more than 20% of the area of each vegetation type in 2015.The abrupt change area for alpine meadow was the largest at 33.06% in 2015, followed by desert steppe and grassland at 27.69% and 25.46%, respectively.The year with the lowest abrupt change areas for all six vegetation types 2011.The abrupt change area was greater than 10% for most years of the study period (Table 2).

Magnitude of Abrupt Changes
The magnitude of the abrupt change in the NDVI for each breakpoint time in the trend component was calculated based on the pixels.Figure 6 shows the spatial distribution of the estimated magnitude of each abrupt change.The magnitudes of the four abrupt changes changed from −0.27 to 0.45 absolute NDVI units.Most of the study area presented abrupt negative changes in the trend component, including grassland, meadow, shrub, and forest areas, especially for the first abrupt change (Figure 6a,b).The negative and positive magnitudes can be used to analyze abrupt vegetation browning and greening, respectively [44].The area percentages of the different changes in the magnitude of each breakpoint for different vegetation types are listed in Table 3.The magnitude was mainly concentrated between −0.06 and −0.01 for the first abrupt change, while it was mainly concentrated between −0.01 and 0.01 for the second to fourth abrupt changes.For each of the six vegetation types, the greening area (magnitude between 0.01 and 0.45) was smaller than the browning area (magnitude between −0.27 and −0.01) for most breakpoint situations (Table 3), which means that most vegetation was undergoing some degree of degradation over certain periods in past decades.For grasslands, meadows, shrubs, and forests, the browning area was much larger than the greening area.
For the entire study period, the abrupt change area was less than 20% for most years (Table 4).The years with the largest abrupt browning areas were 2013 (13.3%) and 2014 (11.7%), followed by 2008 (9.8%) and 2015 (9.3%).The years with the largest abrupt greening areas were 2015 (15.1%) and 2010 (10.8%), followed by 2012 (9.5%) and 2009 (7.9%).The area of abrupt browning was larger than that of greening for most study years.The years with smallest abrupt browning areas were 2005 (1.8%) and 2009 (2.5%), followed by 2011 (3.8%).Small abrupt greening areas were observed in 2003, 2004 and 2013, which were all less than 2.0%.Overall, most of the study area was stable during the study period.For each type of vegetation, significant negative and positive trend changes are illustrated in  For each type of vegetation, significant negative and positive trend changes are illustrated in

Effectiveness of the BFAST Model
The Qilian Mountains represent an ecologically important area in China, and the function and structure of this ecosystem have changed considerably in recent decades [54,55].Accurate detection of vegetation change trends and change times are important for the protection and restoration of vegetation.In this study, we used NDVI to analyze long-term changes in vegetation; however, NDVI has some limitations when exploring vegetation changes, such as being less sensitive to saturation conditions and variable viewing angles [56].However, NDVI serves as a measurement of vegetation greenness and has been widely used on regional and global scales.To analyze the effectiveness of the BFAST results, a linear regression analysis was performed to calculate the slope coefficient of the trend line for each pixel in the study area.The stable area showed a relatively stable trend during the study period, and, according to the BFAST results, the stable area was between 75.6% and 93.5% (Table 4).This result is mostly consistent with the results of the NDVI change trends from 2000 to 2017 (Figure 8), which are based on the slope coefficient of the trend line.The stable area in Figure 8 covers 79.5% of the study area (change trend between −0.005 and 0.005).However, in Figure 8, the area of negative change is much larger than the area of positive change during the study period (Table 4); 20.5% of the area shows an increasing trend (change trend larger than 0.005), and less than 1% of the area shows a decreasing trend (change trend less than −0.005).These findings are seemingly different from those of the BFAST algorithm.However, negative or positive abrupt changes in one year do not mean that the same change trends will occur in the following years, and the trend of vegetation growth may be opposite the magnitude of abrupt changes.As shown in Figure 8, although three abrupt negative changes are observed in the pixel (Figure 9a), the trend of vegetation growth is positive in the following years.The same situation is observed in Figure 9b, in which the magnitude of abrupt changes is positive and the following growth trends are negative.Therefore, the slope can reflect only the overall change trend of vegetation during the study period, whereas the BFAST algorithm can display the specific changes in vegetation accurately during the entire research interval.

Effectiveness of the BFAST Model
The Qilian Mountains represent an ecologically important area in China, and the function and structure of this ecosystem have changed considerably in recent decades [54,55].Accurate detection of vegetation change trends and change times are important for the protection and restoration of vegetation.In this study, we used NDVI to analyze long-term changes in vegetation; however, NDVI has some limitations when exploring vegetation changes, such as being less sensitive to saturation conditions and variable viewing angles [56].However, NDVI serves as a measurement of vegetation greenness and has been widely used on regional and global scales.To analyze the effectiveness of the BFAST results, a linear regression analysis was performed to calculate the slope coefficient of the trend line for each pixel in the study area.The stable area showed a relatively stable trend during the study period, and, according to the BFAST results, the stable area was between 75.6% and 93.5% (Table 4).This result is mostly consistent with the results of the NDVI change trends from 2000 to 2017 (Figure 8), which are based on the slope coefficient of the trend line.The stable area in Figure 8 covers 79.5% of the study area (change trend between −0.005 and 0.005).However, in Figure 8, the area of negative change is much larger than the area of positive change during the study period (Table 4); 20.5% of the area shows an increasing trend (change trend larger than 0.005), and less than 1% of the area shows a decreasing trend (change trend less than −0.005).These findings are seemingly different from those of the BFAST algorithm.However, negative or positive abrupt changes in one year do not mean that the same change trends will occur in the following years, and the trend of vegetation growth may be opposite the magnitude of abrupt changes.As shown in Figure 8, although three abrupt negative changes are observed in the pixel (Figure 9a), the trend of vegetation growth is positive in the following years.The same situation is observed in Figure 9b, in which the magnitude of abrupt changes is positive and the following growth trends are negative.Therefore, the slope can reflect only the overall change trend of vegetation during the study period, whereas the BFAST algorithm can display the specific changes in vegetation accurately during the entire research interval.The effectiveness of the BFAST model in detecting abrupt vegetation changes and long-term change trends can be demonstrated based on known observation sites.It was difficult to obtain these known sites for pixels with a 1 km spatial resolution.In this research, a coal mine pixel with two abrupt changes and six years of negative change trends was detected accurately by the BFAST model, which confirms the credibility of the model.Figure 10 shows a coal mine pixel with abrupt humaninduced changes in the study area from 2000 to 2017.Two breaks were detected for this pixel, which were in 2007 and 2014.The first break was a negative trend change, and the second was a positive trend change (Figure 10b).Before the first abrupt change, the annual maximum NDVI (NDVImax) value of the trend component for this pixel was relatively stable and greater than 0.5.After 2007, the NDVImax value of the pixel was lower than 0.5 and decreased each year, and the NDVImax was less than 0.2 until 2014.Meanwhile, the land cover changes in the pixel from the Landsat images from 2006 to 2011 can also confirm the abrupt vegetation changes (Figure 10c).Coal mines accounted for a small area of the pixel in 2006, although their area increased sharply in 2007 and then increased every until 2011 when more than 60% of the pixel represented coal mines.Since 2014, however, the government has fully promoted the improvement to the ecological environment of mining areas, terminated exploration and mining activities, and invested hundreds of millions of dollars in surface repair and vegetation restoration in these areas.Figure 10b shows that the NDVImax value has increased rapidly since 2014.The NDVImax was only approximately 0.2 in 2014 but exceeded 0.3 in 2016 and exceeded 0.5 in 2017.The land cover changes in the pixels from Landsat images from 2013 to 2017 can also confirm the abrupt changes in vegetation (Figure 10d).

Driving Factors Underlying Abrupt Vegetation Changes in the Qilian Mountains
The factors that affect vegetation growth are very complicated and include environmental factors and anthropogenic activities.Because vegetation growth depends on the background thermal, moisture, and nutrient conditions, other factors that affect the meteorological and nutrient conditions of the growth environment can impact the growth of vegetation, such as the topography, snowpack, solar radiation, atmospheric CO2 concentration, and N deposition.In mountainous areas, topography can control local background thermal and moisture conditions over large scales [57].Snowpack is an important water reservoir for the growth of mountain area vegetation and thus can protect vegetation from damage caused by freezing and improve the temperature of the shallow soil layer.In addition, snow melt in spring enhances soil moisture and vegetation growth [58,59].The effectiveness of the BFAST model in detecting abrupt vegetation changes and long-term change trends can be demonstrated based on known observation sites.It was difficult to obtain these known sites for pixels with a 1 km spatial resolution.In this research, a coal mine pixel with two abrupt changes and six years of negative change trends was detected accurately by the BFAST model, which confirms the credibility of the model.Figure 10 shows a coal mine pixel with abrupt human-induced changes in the study area from 2000 to 2017.Two breaks were detected for this pixel, which were in 2007 and 2014.The first break was a negative trend change, and the second was a positive trend change (Figure 10b).Before the first abrupt change, the annual maximum NDVI (NDVImax) value of the trend component for this pixel was relatively stable and greater than 0.5.After 2007, the NDVImax value of the pixel was lower than 0.5 and decreased each year, and the NDVImax was less than 0.2 until 2014.Meanwhile, the land cover changes in the pixel from the Landsat images from 2006 to 2011 can also confirm the abrupt vegetation changes (Figure 10c).Coal mines accounted for a small area of the pixel in 2006, although their area increased sharply in 2007 and then increased every until 2011 when more than 60% of the pixel represented coal mines.Since 2014, however, the government has fully promoted the improvement to the ecological environment of mining areas, terminated exploration and mining activities, and invested hundreds of millions of dollars in surface repair and vegetation restoration in these areas.Figure 10b shows that the NDVImax value has increased rapidly since 2014.The NDVImax was only approximately 0.2 in 2014 but exceeded 0.3 in 2016 and exceeded 0.5 in 2017.The land cover changes in the pixels from Landsat images from 2013 to 2017 can also confirm the abrupt changes in vegetation (Figure 10d).In addition to multiple environmental factors, anthropogenic activities also affect vegetation, and these activities include grazing, afforestation, policy-driven land use conversions, ecological restoration, tourism development, and other human-induced land cover changes (e.g., mining, urban expansion, and construction of hydropower stations).The Qilian Mountains represent an important pastoral and mineral resource area in China, and meadow and grassland are the two main vegetation types, accounting for more than 60% of the total vegetation areas.The main anthropogenic activities over this region that can influence vegetation are grazing and human-induced land cover changes.The former is a gradual process that can lead to grassland and meadow degradation and even desertification.A considerable challenge was observed in disentangling the relative effects of environmental and human factors on vegetation dynamics because these factors influence vegetation growth at different spatial and temporal scales [60,61].Moreover, only six meteorological stations were present in the study region, and the climate presented large differences because of the complex topography.Therefore, accurately analyzing the contribution of changes in climate factors to sudden changes in vegetation is difficult for large study areas.
The Qilian Mountains are also an important pastoral area, and livestock is the main economic source; thus, significant changes have occurred as a result of long-term grazing.Overgrazing has led to the gradual degradation of grassland areas, which has resulted in increased soil erosion and desertification and declines in water conservation and ecological services in recent decades [55].Grassland and meadow areas account for more than 60% of the total vegetation areas in the study area, and long-term overgrazing is likely the main reason for the widespread abrupt browning

Driving Factors Underlying Abrupt Vegetation Changes in the Qilian Mountains
The factors that affect vegetation growth are very complicated and include environmental factors and anthropogenic activities.Because vegetation growth depends on the background thermal, moisture, and nutrient conditions, other factors that affect the meteorological and nutrient conditions of the growth environment can impact the growth of vegetation, such as the topography, snowpack, solar radiation, atmospheric CO 2 concentration, and N deposition.In mountainous areas, topography can control local background thermal and moisture conditions over large scales [57].Snowpack is an important water reservoir for the growth of mountain area vegetation and thus can protect vegetation from damage caused by freezing and improve the temperature of the shallow soil layer.In addition, snow melt in spring enhances soil moisture and vegetation growth [58,59].
In addition to multiple environmental factors, anthropogenic activities also affect vegetation, and these activities include grazing, afforestation, policy-driven land use conversions, ecological restoration, tourism development, and other human-induced land cover changes (e.g., mining, urban expansion, and construction of hydropower stations).The Qilian Mountains represent an important pastoral and mineral resource area in China, and meadow and grassland are the two main vegetation types, accounting for more than 60% of the total vegetation areas.The main anthropogenic activities over this region that can influence vegetation are grazing and human-induced land cover changes.The former is a gradual process that can lead to grassland and meadow degradation and even desertification.A considerable challenge was observed in disentangling the relative effects of environmental and human factors on vegetation dynamics because these factors influence vegetation growth at different spatial and temporal scales [60,61].Moreover, only six meteorological stations were present in the study region, and the climate presented large differences because of the complex topography.Therefore, accurately analyzing the contribution of changes in climate factors to sudden changes in vegetation is difficult for large study areas.
The Qilian Mountains are also an important pastoral area, and livestock is the main economic source; thus, significant changes have occurred as a result of long-term grazing.Overgrazing has led to the gradual degradation of grassland areas, which has resulted in increased soil erosion and desertification and declines in water conservation and ecological services in recent decades [55].Grassland and meadow areas account for more than 60% of the total vegetation areas in the study area, and long-term overgrazing is likely the main reason for the widespread abrupt browning changes relative to greening changes for most years of the study period.Moreover, mining and hydropower station construction programs may explain certain abrupt vegetation changes [55].Mining can cause partial vegetation damage and induce soil erosion and surface collapse, while hydropower stations in upstream areas may lead to reduced water availability in downstream reaches and thereby affect vegetation growth.
Ecological projects and policies from the national and local government may lead to positive vegetation development.During the study period, the Chinese government enforced a national conservation policy in support of the Grazing Withdrawal Program (GWP) ecological project over the north China grassland in 2003 [53].The aim of the program was to conserve grassland through the implementation of rotational grazing, banning of grazing or conversion of grazing land to cultivated pasture [62], and it has led to obvious declines in the number of livestock.In the same year, the local government of the Qilian Mountain area initiated an ecological protection project and moved 52,000 pastoral herders to the valley area over 3-5 years, which had a great impact on improving vegetation.The enforcement of all policies and projects has helped reduce the grazing pressure over the study area and led to positive changes in the vegetation.The results of this study indicate that the positive abrupt change areas have increased significantly since 2005 (Table 4), especially in 2005, 2009, 2010, 2012 and 2015, and the area percentages of positive changes were larger than those of abrupt changes for all vegetation types.These results are consistent with those of Cai et al. [63], who indicated that ecological protection and restoration projects have reversed the degradation of some areas of the central Tibetan Plateau since 2005.Although the early policy was implemented in 2003, a time lag was observed until the restoration projects became effective [64].Moreover, the effectiveness of ecological restoration projects was based on local socioeconomic, eco-environmental and restoration measures [65].A series of ecological projects have had a significant impact on the reduction in abrupt vegetation change areas.Positive abrupt change areas increased for most study years since 2005, which indicated that heavy grazing pressure has been reduced in the grassland ecosystem of the Qilian Mountains.However, overgrazing still occurs [55], and the area of negative abrupt vegetation changes is larger than that of positive changes throughout the study period (Table 4).

Conclusions
In this study, the spatial and temporal patterns of abrupt vegetation changes in the Qilian Mountains of the Northeast Qinghai-Tibet Plateau were investigated based on the BFAST algorithm using 1 km monthly MODIS NDVI data for the period of 2000-2017.The results indicate that the stable vegetation covered approximately 75.6-90.8% of the study area for each year of the study period, although approximately 80.1% of the vegetation areas experienced one or more abrupt change period from 2000 to 2017, with most of these areas distributed in the southern and northern parts of the study area, especially surrounding Qinghai Lake.Abrupt vegetation changes varied with different vegetation types.In 2015, 24.5% of the studied vegetation experienced abrupt changes, which was the highest percentage in all studied years.In 2011, only 6.5% of the studied vegetation experienced abrupt changes, which was the lowest percentage in all studied years.The positive abrupt change area increased significantly for most study years since 2005.However, for most years of the study period, the abrupt browning change area was larger than the greening area for most vegetation types.Both environmental factors and anthropogenic activities drive abrupt vegetation changes.Long-term overgrazing leads to grassland degradation and likely represents the main cause of the large area of abrupt browning changes relative to the abrupt greening area.Moreover, our results indicate that national ecological protection policies have a significant impact on reducing the area of abrupt vegetation changes, and they have achieved positive effects in the study area since 2005.The findings of this study can serve as basic knowledge for promoting regional ecological protection and improving regional environmental management.
Author Contributions: T.C. conceptualized the study; methodology, formal analysis and writing were conducted by L.G.; X.W. contributed to formal analysis; and H.W. helped to acquire some of the research data.

Figure 1 .
Figure 1.Geographic location of the Qilian Mountains and distribution of vegetation types and weather stations around the study area.

Figure 1 .
Figure 1.Geographic location of the Qilian Mountains and distribution of vegetation types and weather stations around the study area.

Figure 2 .
Figure 2. Schematic diagram representing the processing steps undertaken in this research.

Figure 2 .
Figure 2. Schematic diagram representing the processing steps undertaken in this research.
provides an example of the original data (Yt) and the seasonal (St), trend (T t ) and remainder (et) components of the MODIS NDVI time series at a meadow pixel from 2000 to 2017 that were decomposed using the BFAST procedure.The dashed lines in St and Tt represent the times of abrupt change in NDVI in the seasonal and trend components, respectively.There were three breaks, which were in 2008, 2010, and 2013, and the largest magnitude (−0.072) occurred in 2010.Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 19 provides an example of the original data (Yt) and the seasonal (St), trend (Tt) and remainder (et) components of the MODIS NDVI time series at a meadow pixel from 2000 to 2017 that were decomposed using the BFAST procedure.The dashed lines in St and Tt represent the times of abrupt change in NDVI in the seasonal and trend components, respectively.There were three breaks, which were in 2008, 2010, and 2013, and the largest magnitude (−0.072) occurred in 2010.

Figure 3 .
Figure 3.The seasonal, trend and remainder components of a MODIS NDVI time series over a grassland pixel (Lat.38°49′N, Lon.98°25′E) from 2000 to 2017 that were decomposed using the BFAST procedure.Yt is the original NDVI data, St, Tt and et stand for decomposed seasonal, trend and remainder components.Three breaks were detected for the pixel, which were in 2008, 2010, and 2013.

Figure 3 .
Figure 3.The seasonal, trend and remainder components of a MODIS NDVI time series over a grassland pixel (Lat.38 • 49 N, Lon.98 • 25 E) from 2000 to 2017 that were decomposed using the BFAST procedure.Yt is the original NDVI data, St, Tt and et stand for decomposed seasonal, trend and remainder components.Three breaks were detected for the pixel, which were in 2008, 2010, and 2013.

Figure 4 .
Figure 4. Spatial distribution of the estimated number of abrupt change points in the trend component in the Qilian Mountains in the period of 2000-2017."Other" means no vegetation and managed vegetation areas.

Figure 5 .
Figure 5. Spatial distributions of the estimated timing of abrupt changes in the trend component in the Qilian Mountains in the period of 2000-2017: (a) first abrupt change distribution; (b) second abrupt change distribution; (c) third abrupt change distribution; and (d) fourth abrupt change distribution.Other means no vegetation and managed vegetation area.
0%), followed by 2014 (16.0%), 2013 (15.2%), and 2012 (15.1%).For each vegetation type, the largest abrupt change areas occurred in 2015 and 2010 (except for desert steppe, which occurred in 2013 and 2015), and the abrupt change areas covered more than 20% of the area of each vegetation type in 2015.The abrupt change area for alpine meadow was the largest at 33.06% in 2015, followed by desert steppe and grassland at 27.69% and 25.46%, respectively.

Figure 4 .
Figure 4. Spatial distribution of the estimated number of abrupt change points in the trend component in the Qilian Mountains in the period of 2000-2017."Other" means no vegetation and managed vegetation areas.

Figure 4 .
Figure 4. Spatial distribution of the estimated number of abrupt change points in the trend component in the Qilian Mountains in the period of 2000-2017."Other" means no vegetation and managed vegetation areas.

Figure 5 .
Figure 5. Spatial distributions of the estimated timing of abrupt changes in the trend component in the Qilian Mountains in the period of 2000-2017: (a) first abrupt change distribution; (b) second abrupt change distribution; (c) third abrupt change distribution; and (d) fourth abrupt change distribution.Other means no vegetation and managed vegetation area.

Figure 5 .
Figure 5. Spatial distributions of the estimated timing of abrupt changes in the trend component in the Qilian Mountains in the period of 2000-2017: (a) first abrupt change distribution; (b) second abrupt change distribution; (c) third abrupt change distribution; and (d) fourth abrupt change distribution.Other means no vegetation and managed vegetation area.

Figure 6 .
Figure 6.Spatial distributions of the estimated magnitudes of abrupt changes in the trend component in the Qilian Mountains in the period of 2000-2017: (a) first breakpoint distribution; (b) second breakpoint distribution; (c) third breakpoint distribution; and (d) fourth breakpoint distribution.Other means no vegetation and managed vegetation areas.

Figure 6 .
Figure 6.Spatial distributions of the estimated magnitudes of abrupt changes in the trend component in the Qilian Mountains in the period of 2000-2017: (a) first breakpoint distribution; (b) second breakpoint distribution; (c) third breakpoint distribution; and (d) fourth breakpoint distribution.Other means no vegetation and managed vegetation areas.

Figure 7 .
The largest positive trend change occurred in 2015, and the largest negative change occurred in 2013 for most vegetation.Differences were observed between the negative and positive changes of each vegetation type during the study period.The negative change areas were much larger than the positive change areas for all vegetation types except desert steppe before 2005, and the positive change areas were larger than the negative change areas for all vegetation types in 2005, 2012 (except desert steppe), and 2015 (except grassland).

Figure 7 .
The largest positive trend change occurred in 2015, and the largest negative change occurred in 2013 for most vegetation.Differences were observed between the negative and positive changes of each vegetation type during the study period.The negative change areas were much larger than the positive change areas for all vegetation types except desert steppe before 2005, and the positive change areas were larger than the negative change areas for all vegetation types in 2005, 2012 (except desert steppe), and 2015 (except grassland).

Figure 7 .
Figure 7. Area percentages of negative and positive trend changes for each year of the study period (2000-2017).Negative and positive values indicate abrupt vegetation browning and abrupt vegetation greening, respectively.

Figure 7 .
Figure 7. Area percentages of negative and positive trend changes for each year of the study period (2000-2017).Negative and positive values indicate abrupt vegetation browning and abrupt vegetation greening, respectively.

Figure 8 .
Figure 8. Spatial distribution of NDVI change trends from 2000 to 2017 based on the slope coefficient of the trend line.

Figure 8 .
Figure 8. Spatial distribution of NDVI change trends from 2000 to 2017 based on the slope coefficient of the trend line.

Figure 9 .
Figure 9. BFAST decomposition results of the NDVI time series from 2000 to 2017 for a pixel with a negative magnitude (a) and a positive magnitude (b) of major abrupt changes.

Figure 9 .
Figure 9. BFAST decomposition results of the NDVI time series from 2000 to 2017 for a pixel with a negative magnitude (a) and a positive magnitude (b) of major abrupt changes.

Figure 10 .
Figure 10.Coal mine pixel with abrupt human-induced changes in the study area: (a) magnitude and location of the coalmine; (b) BFAST decomposition results between 2000 and 2017; and (c,d) Landsat images from 2006 to 2011 and from 2013 to 2017, respectively.The red, green, and blue bands of the Landsat images in the growing season with no cloud cover in the pixels were used to identify the vegetation and mining areas.

Figure 10 .
Figure 10.Coal mine pixel with abrupt human-induced changes in the study area: (a) magnitude and location of the coalmine; (b) BFAST decomposition results between 2000 and 2017; and (c,d) Landsat images from 2006 to 2011 and from 2013 to 2017, respectively.The red, green, and blue bands of the Landsat images in the growing season with no cloud cover in the pixels were used to identify the vegetation and mining areas.

Funding:
This research was supported in part by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant No. XDA19070101) and in part by the National Natural Science Foundation of China under grant Nos.41601482 and 41871250.

Table 1 .
Statistical results of the abrupt change pixels and percentage of different vegetation types for each number of breakpoints in the trend component between 2000 and 2017.

Table 2 .
Statistical results of the percentage of abrupt change area of different vegetation types for each year between 2000 and 2017.The highest and lowest area percentages of abrupt changes for each vegetation type are shown in bold.

Table 3 .
Statistical results of area percentages for different changes in the magnitude of different vegetation types between 2000 and 2017.

Table 3 .
Statistical results of area percentages for different changes in the magnitude of different vegetation types between 2000 and 2017.Vegetation Breakpoints −0.27 to −0.06 −0.06 to −0.01 −0.01 to 0.

Table 4 .
Statistical results of area percentages of negative, positive, and stable trend changes for each year of the study period.

Table 4 .
Statistical results of area percentages of negative, positive, and stable trend changes for each year of the study period.