Analyzing Ecological Vulnerability and Vegetation Phenology Response Using NDVI Time Series Data and the BFAST Algorithm

: Identifying ecologically vulnerable areas and understanding the responses of phenology to negative changes in vegetation growth are important bases for ecological restoration. However, identifying ecologically vulnerable areas is di ﬃ cult because it requires high spatial resolution and dense temporal resolution data over a long time period. In this study, a novel method is presented to identify ecologically vulnerable areas based on the normalized di ﬀ erence vegetation index (NDVI) time series from MOD09A1. Here, ecologically vulnerable areas are deﬁned as those that experienced negative changes frequently and greatly in vegetation growth after the disturbances during 2000–2018. The number and magnitude of negative changes detected by the Breaks for Additive Season and Trend (BFAST) algorithm based on the NDVI time-series data were combined to identify ecologically vulnerable areas. TIMESAT was then used to extract the phenology metrics from an NDVI time series dataset to characterize the vegetation responses due to the abrupt negative changes detected by the BFAST algorithm. Focus was given to Jilin Province, a region of China known to be ecologically vulnerable because of frequent drought. The results showed that 13.52% of the study area, mostly in Jilin Province, is ecologically vulnerable. The vulnerability of trees is the lowest, while that of sparse vegetation is the highest. The response of phenology is such that the relative amount of vegetation biomass and length of the growing period were decreased by negative changes in growth for dense vegetation types. The present research results will be useful for the protection of environments being disturbed by regional ecological restoration.


Introduction
Frequent natural disasters and the irrational exploitation of natural resources by human beings have led to a gradual increase in ecosystem vulnerability [1][2][3][4]. The definition of vulnerability by Turner et al. is the degree to which a system, subsystem, or system component is likely to experience harm due to exposure to a hazard, either a perturbation or a stress [5]. In other words, compared with ecologically stable areas, the normal functions and structures of an ecosystem in vulnerable areas are more susceptible to damage. Finally, this damage may subsequently trigger and, consequently, identify ecologically vulnerable areas. Several methods have been proposed to detect abrupt changes in vegetation, such as the Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr) [36], the Breaks for Additive Season and Trend (BFAST) algorithm [37], the Vegetation Change Tracker (VCT) [38], and the Detecting Breakpoints and Estimating Segments in Trend (DBEST) method [39]. The BFAST algorithm, developed by Versbesselt et al. in 2010, is used to identify the breakpoints of the trend components and seasonal components in NDVI time series [40]. With BFAST, high-precision data (with each pixel in the time series fitted individually) can be analyzed without setting a threshold, while detailed information, including the number, time, and magnitude of changes, are estimated. BFAST has been applied to monitor vegetation changes caused by deforestation, forest disturbances, and restoration [41,42]. Considering the phenological parameters of NDVI and trajectory characteristics, some packages, such as HANTS [43], TiSeG [44], TSPT [45], and TIMESAT [46], were developed to extract phenological information. Among the proposed packages, the TIMESAT software package is applied to parameterize the phenological and temporal behavior of NDVI on a seasonal scale. This package integrates a series of filter functions and extraction methods for extracting phenological information, which extract phenological information according to different land types. TIMESAT was developed to generate several seasonal data by fitting a smooth model function to the upper envelope of the RS time-series data. There are three different curve-fitting models: Savitzky-Golay (SG) filtering, asymmetric Gaussian (AG), and double logistic (DL). The SG model has been successfully used for the extraction of phenological parameters due to its high fitting accuracy [47]. In this paper, TIMESAT is used to extract the phenological indicators of different vegetation types.
It is necessary to develop a method for identifying vulnerable areas based on negative changes in vegetation growth and to explore the response mechanisms of the phenology to negative changes. The objectives of this study were (1) to identify areas where the number of negative changes is greater than 0 and the change direction is negative in vegetation growth from the NDVI time series trajectory extracted by the BFAST algorithm in Jilin Province, China during 2000-2018; (2) to compare the stability and sensitivity of different vegetation types; and (3) to explore the response mechanism of the phenology to negative changes by analyzing the relationship between the phenological indicators extracted by TIMESAT and the negative changes decomposed by BFAST.

Study Area
Jilin Province is located in the northeast of China (121 • 38 -131 • 19 E, 40 • 50 -46 • 19 N) and covers an area of approximately 18.74 × 10 4 km 2 . The topography in Jilin Province is rugged in the east and smooth in the west. The terrain is divided into the southeast Changbai Mountains and the northwest Songliao Plain, taking the Dahei Mountains as the boundary, as shown in the LULC map. The study area has a temperate continental monsoon climate. From the southeast to the northwest, the climate transitions from a humid climate to a semi-humid climate and then to a semi-arid climate. The annual total precipitation forms a southeast-northwest gradient, decreasing from 800 to 400 mm. The annual evapotranspiration increases from 1200 to 1800 mm, and the mean temperature decreases from 7 to 3 • C. The climate in the eastern region is humid and rainy, where the annual precipitation is 650 mm. The western part is dry, where the annual precipitation is below 500 mm. The average wind speed is 3.4-3.8 m/s. The natural vegetation of the Jilin Province has been intensively disturbed by drought, leading to a decline in growth. The dominant vegetation types are trees, herbaceous, cropland, and grass. The northwest area of the region belongs to one of the most sensitive areas to disturbances (dominated by drought) in China and plays a vital role in ecological security.

Indexes Calculating and LULC Map
MOD09A1, a MODIS reflectance data product with a spatial resolution of 500 m composited over 8 days for the period from April 2000 to September 2018, was used in this study. This product was downloaded from the NASA website (https:/search.earthdata.nasa.gov/search). The April-September period was selected because it corresponds to the vegetation growth period for the study area. The raw sinusoidal projection was converted to an Albers equal area projection with an ellipsoid WGS-84 projection by the MODIS Reprojection Tool (MRT). Pixels contaminated by clouds were masked by the MODIS quality assurance (QA) band. NDVI time-series were calculated by the Google Earth Engine (GEE) platform for band 1 (Red 620-670 nm) and band 2 (NIR 841-875 nm) of MOD09A1.
The ESA CCI-LC maps were released as the time series of annual global land cover maps at 300 m from 1992 to 2018 by merging multiple earth observation products (MERIS, SPOT-VGT, and PROBA-V) based on the GlobCover products of the ESA, with an overall accuracy of about 80% [50]. Global ESA CCI-LC maps describing the land surface with 22 classes were provided and were defined using the United Nations Food and Agriculture Organization's (UN FAO) Land Cover Classification System (LCCS). The LULC of Jilin Province in 2018 was downloaded (https://cds.climate.copernicus.eu) and projected into the Albers equal area projection to remain consistent with the NDVI data. Twelve land types are included in the study area, including cropland, herbaceous land, mosaic natural vegetation, tree cover, mosaic tree and shrub, mosaic herbaceous cover, shrub land, grassland, sparse vegetation, bare areas, urban areas, and water bodies. The mosaic herbaceous cover and mosaic tree/shrub (combined) were removed due to their small area (the sum of both was 0.6%) and the fact that they also belong to mixed land. The urban areas and water bodies were omitted, as they were not the focus of this research. Bare areas are areas with no dominant vegetation cover of at least 90%. Thus, bare land was excluded. Sparse vegetation included tree, shrub, and herbaceous land <15%. Finally, the seven remaining classes are shown in Figure 1. Cropland and sparse vegetation are distributed across the northwest area of Jilin Province. Cropland and herbaceous areas are concentrated and contiguous in the central region. The vegetation type in the south is dominated by trees, followed by shrubs. We differentiated between the reclassified land use cover maps in 2000 and 2018 and found that the area of change between the two was less than 5%. This change mainly comes from the conversion of grassland to cropland in Baicheng. Therefore, it was assumed that the land covers remained essentially unchanged during the study period. The LULC was used to statistically analyze the changes in information for each type of BFAST and as the input data of TIMESAT to extract the phenological indicators of different vegetation types.

Ecologically Vulnerable Areas and Phenological Responses Methods
Identifying vulnerable areas and exploring the relationship between negative changes and phenology were conducted in two steps. In the first step, areas with frequent and large magnitudes of negative changes in vegetation growth were screened by the number of changes and the magnitude of changes extracted using BFAST from the NDVI time series. Ecological vulnerability was deduced from the number of negative changes. In the second step, changes in phenological indicators obtained from the TIMESAT and NDVI trends decomposed by BFAST were compared by the year of negative change to explore the response of the phenology to negative changes. The flowchart is shown in Figure 2.

Ecologically Vulnerable Areas and Phenological Responses Methods
Identifying vulnerable areas and exploring the relationship between negative changes and phenology were conducted in two steps. In the first step, areas with frequent and large magnitudes of negative changes in vegetation growth were screened by the number of changes and the magnitude of changes extracted using BFAST from the NDVI time series. Ecological vulnerability was deduced from the number of negative changes. In the second step, changes in phenological indicators obtained from the TIMESAT and NDVI trends decomposed by BFAST were compared by the year of negative change to explore the response of the phenology to negative changes. The flowchart is shown in Figure 2.

Identifying Ecologically Vulnerable Areas
The abrupt changes in vegetation growth were obtained by BFAST based on the NDVI time series during 2000-2018. Breakpoints detected by BFAST indicate a sudden change of the NDVI, while the trend of the NDVI time series is different on both sides of the breakpoint. The number of changes represents the number of breakpoints in the NDVI trend decomposed by BFAST. The magnitude of change is the difference between the NDVI values before and after the breakpoint, while the time of change refers to the year of the change. The frequency of negative changes refers to the pixels whose numbers of breakpoints are greater than or equal to 1. Great negative changes mean that the magnitude of change is less than −0.05. Then, pixels with frequent and great negative changes were superimposed to identify ecologically vulnerable areas, while the vulnerability level was determined by the number of negative changes in ecologically vulnerable areas.  Figure 2. Flowchart of the identification of ecologically vulnerable areas and exploration of the phenological response mechanism for disturbances using Breaks for Additive Season and Trend (BFAST) and TIMESAT from the remote sensing (RS) time series.

Identifying Ecologically Vulnerable Areas
The abrupt changes in vegetation growth were obtained by BFAST based on the NDVI time series during 2000-2018. Breakpoints detected by BFAST indicate a sudden change of the NDVI, while the trend of the NDVI time series is different on both sides of the breakpoint. The number of changes represents the number of breakpoints in the NDVI trend decomposed by BFAST. The magnitude of change is the difference between the NDVI values before and after the breakpoint, while the time of change refers to the year of the change. The frequency of negative changes refers to the pixels whose numbers of breakpoints are greater than or equal to 1. Great negative changes mean that the magnitude of change is less than −0.05. Then, pixels with frequent and great negative changes were superimposed to identify ecologically vulnerable areas, while the vulnerability level was determined by the number of negative changes in ecologically vulnerable areas. The statistics of the BFAST results for different vegetation types according to the ESA CCI-LC and BFAST results were used to explore the vulnerability of different vegetation types.
BFAST is an iterative algorithm used to detect changes by decomposing time-series NDVI into three components: the trend, seasonality, and remainder components. The decomposition model is as follows: where is the NDVI at time t; and Tt, St, and et, respectively, represent the trend, seasonality, and remainder components. p represents the number of observations, which is 380 in this paper. Tt is fitted by piecewise linear models with specific slopes βi and intercepts αi on different segments, as in formula (2): where m represents the number of breakpoints in the trend component, which is the number of changes suffered by the area. is the time at the breakpoint i, and St is fitted as the piecewise harmonic on different n + 1 segments as shown in model (3): BFAST is an iterative algorithm used to detect changes by decomposing time-series NDVI into three components: the trend, seasonality, and remainder components. The decomposition model is as follows: where Y t is the NDVI at time t; and T t , S t , and e t , respectively, represent the trend, seasonality, and remainder components. p represents the number of observations, which is 380 in this paper. T t is fitted by piecewise linear models with specific slopes β i and intercepts α i on different segments, as in formula (2): where m represents the number of breakpoints in the trend component, which is the number of changes suffered by the area. τ i is the time at the breakpoint i, and S t is fitted as the piecewise harmonic on different n + 1 segments as shown in model (3): where n is the number of breakpoints in the seasonal component, τ j is the time at the breakpoint j, and k is the order of the harmonic function, which was set as 3 to more accurately detect complex changes, such as double or even triple-vegetation patterns, within a year. As a result, not only simple seasonal patterns, such as a yearly growth cycle, but also complex growth patterns could be accurately detected. The frequency f was set as 20 for annual observations within a year of the time-series NDVI data. The iteration process was initialized by estimating the seasonal component S t with the STL algorithm. First, the NDVI time series data were decomposed into their seasonal, trend, and remainder components with the seasonal-trend decomposition procedure based on Loess (STL). Then, an ordinary least squares moving sum (OLS-MOSUM) test was used to determine if any breakpoints were present in the time series. If the OLS-MOSUM test indicated a significant change (p < 0.05), then the numbers and locations of the breakpoints were estimated separately for the seasonal and trend components using OLS fitting. A third-order harmonic model was automatically fitted by BFAST. The result was a set of piecewise season-trend models that minimizes the errors of the whole time series. The difference between the intercept and slope terms of consecutive models was used to calculate the magnitude of change between the breakpoints. The number of changes is the count of the breakpoints. The time of the change is the time corresponding to the breakpoint location.
The BFAST results are shown in Figure 3, which provides an example of the original data (Y t ) and the seasonal (S t ), trend (T t ), and remainder (e t ) components of the NDVI time series for a pixel of grassland (Lat. 46 • 39 N, Lon. 123 • 13 E) from 2000 to 2018 that was decomposed by the BFAST algorithm. The dashed lines in Tt represent the negative changes in the NDVI trend components. One breakpoint in the NDVI trend was detected by the OSL-MOSUM algorithm. The trend was fitted by y = 0.012x + 0.3369 and y = 0.021x − 1.1761. The NDVI change between the two equations at breakpoint was −0.154. Thus, there was a breakpoint that occurred in 2008, which is, therefore, the year of change. The magnitude of change was −0.154. At the breakpoint, the multi-year trend of NDVI changes and the NDVI values experienced a sudden change, which means that the value of NDVI did not follow the change trend before the breakpoint and manifest as a sudden drop or rise. The number of changes was 1.
changes, such as double or even triple-vegetation patterns, within a year. As a result, not only simple seasonal patterns, such as a yearly growth cycle, but also complex growth patterns could be accurately detected. The frequency was set as 20 for annual observations within a year of the timeseries NDVI data. The iteration process was initialized by estimating the seasonal component St with the STL algorithm.
First, the NDVI time series data were decomposed into their seasonal, trend, and remainder components with the seasonal-trend decomposition procedure based on Loess (STL). Then, an ordinary least squares moving sum (OLS-MOSUM) test was used to determine if any breakpoints were present in the time series. If the OLS-MOSUM test indicated a significant change (p < 0.05), then the numbers and locations of the breakpoints were estimated separately for the seasonal and trend components using OLS fitting. A third-order harmonic model was automatically fitted by BFAST. The result was a set of piecewise season-trend models that minimizes the errors of the whole time series. The difference between the intercept and slope terms of consecutive models was used to calculate the magnitude of change between the breakpoints. The number of changes is the count of the breakpoints. The time of the change is the time corresponding to the breakpoint location.
The BFAST results are shown in Figure 3, which provides an example of the original data (Yt) and the seasonal (St), trend (Tt), and remainder (et) components of the NDVI time series for a pixel of grassland (Lat. 46°39′ N, Lon. 123°13′ E) from 2000 to 2018 that was decomposed by the BFAST algorithm. The dashed lines in Tt represent the negative changes in the NDVI trend components. One breakpoint in the NDVI trend was detected by the OSL-MOSUM algorithm. The trend was fitted by y = 0.012x + 0.3369 and y = 0.021x − 1.1761. The NDVI change between the two equations at breakpoint was −0.154. Thus, there was a breakpoint that occurred in 2008, which is, therefore, the year of change. The magnitude of change was −0.154. At the breakpoint, the multi-year trend of NDVI changes and the NDVI values experienced a sudden change, which means that the value of NDVI did not follow the change trend before the breakpoint and manifest as a sudden drop or rise. The number of changes was 1.

BFAST Result Verification
The negative abrupt changes in the NDVI trend detected by BFAST were used to represent the growth change of the vegetation after being disturbed. In this paper, the two drought events detected by SPEI in 2007 and 2010 are compared with the breakpoints in the NDVI trend obtained from BFAST. The drought conditions can be characterized by SPEI, which includes the drought events, the trend of drought, and the relative intensity of the drought in two events. These drought conditions were compared with the breakpoints, the time of occurrence, the trend, and the magnitude of the change

BFAST Result Verification
The negative abrupt changes in the NDVI trend detected by BFAST were used to represent the growth change of the vegetation after being disturbed. In this paper, the two drought events detected by SPEI in 2007 and 2010 are compared with the breakpoints in the NDVI trend obtained from BFAST. The drought conditions can be characterized by SPEI, which includes the drought events, the trend of drought, and the relative intensity of the drought in two events. These drought conditions were compared with the breakpoints, the time of occurrence, the trend, and the magnitude of the change in the NDVI trend. The negative changes in vegetation growth can be used to indicate the disturbed changes in vegetation growth.
According to official records, external disturbances in the northwest region of Jilin Province were dominated by drought during these years, particularly in 2007 and 2010. Therefore, the Qian'an meteorological station in the northwest was selected to monitor the drought. The precipitation at the site in 2007 and 2010 was 308.3 and 244.4 mm, respectively, which is lower than the average annual precipitation of 414.37 mm. The SPEI data over 12 months were computed with the R package "SPEI", and the SPEI values for April to September were selected during 2000-2018. The NDVI trend detected by BFAST for the pixels near the station was extracted for a comparison with SPEI. The consistency of the magnitude and time of the negative changes detected by BFAST and the degree and time of the drought detected by SPEI were compared to verify the relationship between negative vegetation changes and disturbances. The analysis of this comparison is provided in Section 3.2.

Exploring the Phenological Response Mechanism to Negative Changes
To explore the phenological responses to the negative changes of vegetation, three steps were implemented. First, based on the year of change extracted by BFAST, the year in which each type of vegetation had the most negative changes was determined, and the average NDVI time series of the negative change pixels in that year was calculated. At the same time, the 2000-2018 phenological indicators of these pixels were extracted. The average value was then calculated to obtain the phenological indicators of each type of vegetation from 2000 to 2018. Then, the overall time series similarity of the NDVI trends and the phenological indicators were compared to screen the phenological indicators related to growth. Finally, the changes in phenological indicators and NDVI trends during the year of change were compared to explore the responses of the phenology to negative changes.
TIMESAT was developed to extract phenological indicators based on time-series NDVI data derived from satellite RS images [51]. This process can be summarized in three steps.
(1) The number of seasons and their approximate years were identified. For all vegetation types in the study area, one year was taken as a growing season. Every interval of 20 data points was considered a growth cycle, and 19 growth phases were counted between 2000 and 2018.
(2) The best fitting model was selected according to the characteristics of the time-series trajectory.
In TIMESAT, there are three fitting models: double logistic, asymmetric Gaussian, and SG filtering. The SG method was chosen, with which subtle and rapid changes in the simulation of local variations can be captured [52]. The expression for SG filtering is as follows: where Y * j is the reconstructed time-series NDVI data, C i is the filtering coefficient, Y j is the original NDVI data, N is the number of data points in the sliding window (N = 2m + 1), and 2m + 1 is the filtering window width. Two parameters are required in SG filtering: the filter window width and the order of the polynomial for the smooth fitting process. The filter coefficients of the SG filter were determined by an unweighted linear least-squares regression and a second-order polynomial model. In this paper, the size of the filter window was set to 5, and the order of the fitting polynomial was taken as 2.
(3) The time series trajectory characteristics of NDVI were extracted by BFAST and TIMESAT.
The difference is that BFAST detects the inter-annual variation characteristics of the long-term time series, while TIMESAT detects the local variation characteristics of the NDVI time series (that is, the seasonal changes). The phenological indicators extracted by TIMESAT were divided into 3 categories according to their mean growth phase: start of season (SOS) and end of season (EOS). The maximum during the growth phase is taken as the peak and amplitude, and the integral during the growth phase includes the large integral (the L. integral) and small integral (the S. integral). The S. integral is equal to the L. integral minus the integral of the base line from SOS to EOS. In this study, the SOS and EOS were determined using a fixed threshold approach, with which the smoothed 8-day NDVI reached 25% of the mean amplitude for each growth phase for each pixel. The amplitude was calculated as the peak minus the minimum of the smoothed NDVI values. The minimum of the smoothed NDVI values was set to zero if the smoothed NDVI values were negative. The length of season (LOS) was calculated as the EOS minus SOS. The position of each indicator in the time-series trajectory is shown in Figure 4.
which the smoothed 8-day NDVI reached 25% of the mean amplitude for each growth phase for each pixel. The amplitude was calculated as the peak minus the minimum of the smoothed NDVI values. The minimum of the smoothed NDVI values was set to zero if the smoothed NDVI values were negative. The length of season (LOS) was calculated as the EOS minus SOS. The position of each indicator in the time-series trajectory is shown in Figure 4.

Identification of Ecologically Vulnerable Areas by BFAST
The number and magnitude of abrupt changes in vegetation indicate the instability and sensitivity of the ecosystem, respectively. The number, magnitude, and year of change were detected in the NDVI trend component decomposed using BFAST for each pixel. The spatial distribution and area percentage of the number of changes are shown in Figure 5a,g, respectively. The areas that experienced changes were mainly concentrated in the northwestern part of Jilin Province and accounted for 23.85% of the total area. This indicates that the ecosystem in the northwest is unstable against environmental changes. Using numbers from 1 to 5, as the area percentage decreased, the number increased. The areas with a score of 1 for negative changes were distributed in Siping, Changchun, and Songyuan, accounting for 55.56% of the changed area. Thus, the ecosystems in these areas are relatively stable compared to the areas with numbers greater than 1. The areas with numbers greater than 1 were mainly located in Baicheng city and Songyuan city, as shown in Figure 5d. This suggests that the system stability levels of Baicheng and Songyuan are the worst in the whole region. The larger the magnitude of change, the more sensitive the ecosystem is to external changes. The magnitude of change is shown in Figure 5b,h. The magnitude varied from −2 to 2, within which a magnitude of variation between −0.4 and 0.4 accounted for 86% of the area of change. The percentage of magnitude close to 0 was the highest, which indicates that the ecosystems in most areas are not sensitive to environmental changes. The areas that experienced a negative change, located in Baicheng and Songyuan, occupied 52% of the area of change, which is greater than the areas of positive change, which were located in Siping and Tonghua. This suggests that the areas with the worst anti-disturbance ability in the ecosystem produce a decline in vegetation growth. Figure 5c

Identification of Ecologically Vulnerable Areas by BFAST
The number and magnitude of abrupt changes in vegetation indicate the instability and sensitivity of the ecosystem, respectively. The number, magnitude, and year of change were detected in the NDVI trend component decomposed using BFAST for each pixel. The spatial distribution and area percentage of the number of changes are shown in Figure 5a,g, respectively. The areas that experienced changes were mainly concentrated in the northwestern part of Jilin Province and accounted for 23.85% of the total area. This indicates that the ecosystem in the northwest is unstable against environmental changes. Using numbers from 1 to 5, as the area percentage decreased, the number increased. The areas with a score of 1 for negative changes were distributed in Siping, Changchun, and Songyuan, accounting for 55.56% of the changed area. Thus, the ecosystems in these areas are relatively stable compared to the areas with numbers greater than 1. The areas with numbers greater than 1 were mainly located in Baicheng city and Songyuan city, as shown in Figure 5d. This suggests that the system stability levels of Baicheng and Songyuan are the worst in the whole region. The larger the magnitude of change, the more sensitive the ecosystem is to external changes. The magnitude of change is shown in Figure 5b,h. The magnitude varied from −2 to 2, within which a magnitude of variation between −0.4 and 0.4 accounted for 86% of the area of change. The percentage of magnitude close to 0 was the highest, which indicates that the ecosystems in most areas are not sensitive to environmental changes. The areas that experienced a negative change, located in Baicheng and Songyuan, occupied 52% of the area of change, which is greater than the areas of positive change, which were located in Siping and Tonghua. This suggests that the areas with the worst anti-disturbance ability in the ecosystem produce a decline in vegetation growth. Figure 5c  The areas with abrupt declines (negative changes) in vegetation growth after disturbances during the monitoring period were considered to be ecologically vulnerable areas due to their sensitivity and instability. Since the negative changes are mainly concentrated from −0.05 to −0.4, areas with a magnitude of change less than −0.05 and a number of changes greater than 0 were considered ecologically vulnerable areas ( Figure 6). Vulnerability was divided into five categories according to the number of changes: low vulnerability, slight vulnerability, moderate vulnerability, high vulnerability, and serious vulnerability, as the number of changes increased, as shown in Table 1. The areas with abrupt declines (negative changes) in vegetation growth after disturbances during the monitoring period were considered to be ecologically vulnerable areas due to their sensitivity and instability. Since the negative changes are mainly concentrated from −0.05 to −0.4, areas with a magnitude of change less than −0.05 and a number of changes greater than 0 were considered ecologically vulnerable areas ( Figure 6). Vulnerability was divided into five categories according to the number of changes: low vulnerability, slight vulnerability, moderate vulnerability, high vulnerability, and serious vulnerability, as the number of changes increased, as shown in Table 1.    Figure 6 shows that ecologically vulnerable areas accounted for 13.52% of the total region, mainly concentrated in the northwest (Baicheng, Songyuan) of Jilin Province, where disturbances are caused by drought. The vulnerability of the study area was at a low level, which indicated that resistance against drought disasters was common. Areas of low vulnerability appeared in the southwest of Baicheng and the east of Songyuan, accounting for 46.9% of the vulnerable area; these areas were dominated by cropland and herbaceous. In areas with less precipitation but a large effective irrigation area to reduce the impact of regional drought, wind erosion weakened due to the cropland's shelterbelts. Areas of slight vulnerability accounted for 37.91% of the vulnerable area and were found in the northwest of Baicheng and to the north of Songyuan, where there was sparse vegetation mixed with cropland and a mosaic of natural vegetation. The precipitation in the region was relatively low, and the irrigation conditions were poor. The landscape was fragmented, the surface structure was simple, and the wind erosion and irrigation conditions were slightly worse in this region than those in areas of low vulnerability. Areas of moderate, high, and serious vulnerability accounted for 15.19% of the vulnerable area and were located in the northeast of Songyuan, which is dominated by sparse vegetation. The wind erosion and soil salinization in this area were the most serious in the entire region. These areas usually suffered from hazardous drought and the most serious vegetation degradation. The poor ecological environment in this area led to vegetation degradation and the poorest ability to cope with drought.

BFAST Result Verification
A comparison between the NDVI trend and SPEI at the sites is illustrated in Figure 7. In the 5 months of the study (May to September), the NDVI trends corresponded well to the SPEI in terms of their trend and slope. Before April 2007, the drought at the site became serious, leading to a gradual decline in the NDVI trend. The relatively low intensity of the drought (severe drought) period from April to September 2007 (SPEI < −1.5) led to the first drop (abrupt change) of 0.04 in the NDVI trend within a month. After September 2007, the drought at the site gradually weakened as the SPEI increased, which led to an upward NDVI trend. The second highest intensity drought (extreme drought) period from May 2009, to April 2010 (SPEI < −2) led to a dramatic drop in the NDVI trend. Within two years, the average NDVI value decreased by 0.1. After May 2010, the humidity of the site gradually increased, leading to an upward trend in the NDVI. The decreasing trends were all significant (p < 0.001). In general, the negative changes of NDVI detected by BFAST reflected changes in vegetation after being disturbed.

Analysis of the Vulnerability of Different Vegetation Types
The vulnerability of different vegetation types against drought was also explored. Figure 8a shows the number of negative changes that represents the instability of each vegetation type. The

Analysis of the Vulnerability of Different Vegetation Types
The vulnerability of different vegetation types against drought was also explored. Figure 8a shows the number of negative changes that represents the instability of each vegetation type. The tree cover, which experienced a negative change of 1, accounted for 99% of the change areas during 2000−2018. Numbers from 1 to 5 were determined for sparse vegetation, while areas with 1 and 2 accounted for 40% and 34%, respectively. The area percentage with a score greater than 1 for other vegetation types was less than that of the sparse vegetation. Tree cover was relatively stable and changed little in response to drought, while sparse vegetation changed frequently and was thus the most unstable in resisting external drought. Figure 7. The 5-month (May-September) standardized precipitation evapotranspiration index (SPEI) (orange bar) and NDVI trend (line) components decomposed by BFAST at Qianan station. C, D, E, F, respectively, represent the NDVI trend for 4 pixels near the site extracted by BFAST.

Analysis of the Vulnerability of Different Vegetation Types
The vulnerability of different vegetation types against drought was also explored. Figure 8a shows the number of negative changes that represents the instability of each vegetation type. The tree cover, which experienced a negative change of 1, accounted for 99% of the change areas during 2000−2018. Numbers from 1 to 5 were determined for sparse vegetation, while areas with 1 and 2 accounted for 40% and 34%, respectively. The area percentage with a score greater than 1 for other vegetation types was less than that of the sparse vegetation. Tree cover was relatively stable and changed little in response to drought, while sparse vegetation changed frequently and was thus the most unstable in resisting external drought.  The magnitude of change was then compared to determine the sensitivity of different vegetation types. Figure 8b illustrates the mean and the 95% confidence range of the negative magnitude. For the negative changes, sparse vegetation (−0.429) had the lowest mean, while tree cover (−0.092) had the highest mean for the negative magnitude. The 95% confidence range of sparse vegetation was the highest at 0.623, while that of the tree cover was the lowest at 0.175. The results indicated that the order of sensitivity of negative changes from small to large was tree cover, cropland, herbaceous land, shrubland, mosaic natural vegetation, grassland, and sparse vegetation. In general, the strongest and the weakest abilities to resist external drought belonged to tree cover and sparse vegetation, respectively. Figure 9 shows the percentage of the negative change area for each type of vegetation compared to the total negative change area during 2000−2018. The years with negative changes ranged from 2002 to 2016. The area proportion from large to small was herbaceous land, cropland, mosaic natural vegetation, grassland, shrubland, sparse vegetation, and tree cover. The year with the largest negative change area ranged from 2010 to 2014, which indicates that the northwest area of Jilin suffered more drought during 2010 to 2014 than during the other years. Figure 9 shows the percentage of the negative change area for each type of vegetation compared to the total negative change area during 2000−2018. The years with negative changes ranged from 2002 to 2016. The area proportion from large to small was herbaceous land, cropland, mosaic natural vegetation, grassland, shrubland, sparse vegetation, and tree cover. The year with the largest negative change area ranged from 2010 to 2014, which indicates that the northwest area of Jilin suffered more drought during 2010 to 2014 than during the other years.

Vegetation Phenological Responses to Negative Changes
Finding the pattern underlying the most negative changes for each vegetation type was a prerequisite for exploring the response mechanism of phenology to negative changes. First, the year with the largest proportion was selected according to Figure 9, as shown in Table 2. Table 2. Statistics of the negative changes with the largest areas for each vegetation type.

Vegetation Phenological Responses to Negative Changes
Finding the pattern underlying the most negative changes for each vegetation type was a prerequisite for exploring the response mechanism of phenology to negative changes. First, the year with the largest proportion was selected according to Figure 9, as shown in Table 2. Table 2. Statistics of the negative changes with the largest areas for each vegetation type. Then, the correlation between the NDVI trend representing the growth detected by BFAST and the phenological indicators obtained by TIMESAT were used to screen the phenological indicators related to growth. As shown in Table 3, for the growth phase, the correlation coefficient of SOS was higher than that of EOS in the cropland, mosaic natural vegetation, and grassland. The EOS correlation was higher than that of SOS and length for the herbaceous land, tree cover, shrubland, and sparse vegetation. For all types, the correlation of the peak was higher than that of the amplitude at its maximum. The correlation of the L. integral was higher than that of the S. integral in the integral value. The EOS, SOS, peak, and L. integral were selected as the sensitivity indicators related to growth to explore the phenological responses to negative changes. The peak and L. integral were the indicators for vegetation biomass [53], while the EOS and SOS were indicators of the growth phase.

Vegetation Types Number of Changes Year of Change
The trajectory of the NDVI trend component decomposed by BFAST and the sensitive phenological indicators are shown in Figure 10. For the overall change, the NDVI trend was more similar to the L. integral than the peaks for cropland, herbaceous land, shrubland, mosaic natural vegetation, tree cover, and sparse vegetation for the vegetation biomass. For the growth phase, the EOS was consistent with the NDVI trend for all types. The results also demonstrated that the L. integral and EOS were more consistent with growth in the vegetation biomass and the growth phase. For local changes, the NDVI trend and the sensitive phenological indicators were compared by the years of negative changes. In terms of vegetation biomass, the L. integral decreased due to negative changes of the NDVI trend for cropland, herbaceous land, mosaic natural vegetation, tree cover, shrubland, and grassland, while the L. integral of sparse vegetation increased. Thus, the negative changes had a positive effect on the L. integral for dense vegetation types (Figure 10a-f). For sparse vegetation, the negative change negatively impacted the L. integral. Further, the length of the growth phase was shortened under negative changes for cropland, herbaceous land, tree cover, mosaic natural vegetation, shrubland, and grassland. However, the length for sparse vegetation was extended. For the dense vegetation types, a positive impact on the length of the growth phase was caused by the negative change. For the sparse vegetation, the negative change negatively affected the length.
For the peak changes, both EOS and SOS decreased, while the growth phase advanced. The time of the peak was earlier than that of the previous year, whereas the peak declined for mosaic natural vegetation, shrubland, grassland, and sparse vegetation, as shown in Figure 10c,e-g. For cropland and herbaceous land, EOS and SOS both increased, and the growth phase was delayed. The time of the peak was late, and an increase in the peak is shown in Figure 10a,b. When the decrease in EOS was greater than the increase in SOS, the growth phase had advanced. When the year of the peak was earlier, the peak declined for tree cover, as shown in Figure 10d. In summary, for all vegetation types, the variation of the peak depended on the peak time. Specifically, the peak time was earlier when the peak declined. Otherwise, the time of the peak was delayed as the peak rose.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 22 Then, the correlation between the NDVI trend representing the growth detected by BFAST and the phenological indicators obtained by TIMESAT were used to screen the phenological indicators related to growth. As shown in Table 3, for the growth phase, the correlation coefficient of SOS was higher than that of EOS in the cropland, mosaic natural vegetation, and grassland. The EOS correlation was higher than that of SOS and length for the herbaceous land, tree cover, shrubland, and sparse vegetation. For all types, the correlation of the peak was higher than that of the amplitude at its maximum. The correlation of the L. integral was higher than that of the S. integral in the integral value. The EOS, SOS, peak, and L. integral were selected as the sensitivity indicators related to growth to explore the phenological responses to negative changes. The peak and L. integral were the indicators for vegetation biomass [53], while the EOS and SOS were indicators of the growth phase. The trajectory of the NDVI trend component decomposed by BFAST and the sensitive phenological indicators are shown in Figure 10. For the overall change, the NDVI trend was more similar to the L. integral than the peaks for cropland, herbaceous land, shrubland, mosaic natural vegetation, tree cover, and sparse vegetation for the vegetation biomass. For the growth phase, the EOS was consistent with the NDVI trend for all types. The results also demonstrated that the L. integral and EOS were more consistent with growth in the vegetation biomass and the growth phase. For local changes, the NDVI trend and the sensitive phenological indicators were compared by the years of negative changes. In terms of vegetation biomass, the L. integral decreased due to negative changes of the NDVI trend for cropland, herbaceous land, mosaic natural vegetation, tree cover, shrubland, and grassland, while the L. integral of sparse vegetation increased. Thus, the negative changes had a positive effect on the L. integral for dense vegetation types (Figure 10a-f). For sparse vegetation, the negative change negatively impacted the L. integral. Further, the length of the growth phase was shortened under negative changes for cropland, herbaceous land, tree cover, mosaic natural vegetation, shrubland, and grassland. However, the length for sparse vegetation was

Spatial Distribution of the Ecologically Vulnerable Area
The results in Section 3.1 show that the ecologically vulnerable areas were located in the northwest area of Jilin Province. Compared with the southeastern part, which is covered by trees, the northwestern part was mainly covered by cropland and sparse vegetation, where more frequent and severe drought occurred. This result is consistent with the fact that the stability of the trees was greater than the stability of sparse vegetation and cropland, as shown in Section 3.3. Moreover, the results of this paper are consistent with previous research [54]. Consequently, the method based on abrupt change detection from RS time-series data is able to identify ecologically vulnerable areas. Thus, the regional disaster prevention capabilities were identified effectively by BFAST. The northwest area contained arid and semi-arid areas and the vulnerable ecological zone of cropland and grass in northeastern China. This study area belongs to the junction between a temperate continental climate and a temperate monsoon climate, which is the reason for the complex community structure, high degree of environmental heterogeneity due to poor anti-interference abilities, and climatic sensitivity of the ecological system, as well as the space-time fluctuations with a high frequency and weak edge effect in the area. In the agricultural belt of the northwest, as a result of long-term continuous cropping and the excessive application of pesticides and fertilizers, soil fertility has declined. The degradation and erosion of black soil have become gradually more serious under the impact of non-point-source pollution and wind erosion. The surface temperature and biological distribution structure have also been changed by human interference. However, soil salinization and grassland degradation have reduced drought resistance in the northern region, which is located in the agricultural and pastoral ecotone zone. At the same time, blind reclamation, wetland environment deterioration, function decline, and the extent of vulnerable ecosystems were much worse.

The Response of Vegetation Phenology to Negative Changes
As shown in Figure 10, parts of the NDVI trends were not similar to the phenology indicators. In this paper, the relative amount of vegetation biomass and length of the growing period were decreased by negative changes. However, the peak and SOS were less affected by negative changes. This indicates that the impact of negative change on vegetation phenology has long-term effects. When vegetation growth suddenly drops, the balanced state of the interactive process between vegetation and the environment is broken. Therefore, vegetation constantly adjusts itself to adapt to the changing environment. For example, vegetation closes its stomata to reduce water loss and fit the drought status. Vegetation needs a period of recovery before reaching its next stable period. Therefore, negative changes had a more significant effect on the cumulative values than the instantaneous values. Additionally, when a negative change occurred after the time of peak, the effect of the negative change was also independent of the peak. The peak declined as the peak time advanced, while the peak rose as the peak time was delayed. Advancement of the peak time indicated that the vegetation was not located in the correct ecological environment and was much susceptible to the stress of water and nutrients. Thus, the peak declined. Otherwise, when the peak time was delayed, there was enough available soil water, heat, and fertilizer to improve photosynthesis efficiency. Consequently, the peak rose. In general, the results of this study show that when the changes in phenology are known, changes in the surrounding environment can be inferred. When the length of the growing period and the relative amount of vegetation biomass decreased significantly, vegetation growth was under environmental stress, indicating that the environment could gradually deteriorate.

Limitations and Future Work
In this paper, the negative changes of NDVI were used to identify the ecologically vulnerable areas under drought, although this method is also applicable to other disturbances that cause declines in vegetation growth, such as insect infestations and flood disasters, because the negative changes in the NDVI reflect the sensitivity and stability of vegetation being disturbed. However, specific disturbances need to be compared with the results of BFAST to verify that the vegetation change is, in fact, caused by the disturbance.
Three major limitations of this study should be discussed further. First, vulnerable areas were identified based only on negative changes in vegetation growth. This limitation led to overlooking areas of low vegetation cover. So, areas with less than 10% vegetation cover were excluded in the vulnerable area analysis. This was a necessary trade-off because it is difficult to assess ecological vulnerability in areas of low vegetation cover. This, however, should not be a great concern as only a small proportion of the study area (1.54%) was characterized by low vegetation cover. Second, in this paper, the relationship of change between phenology and negative changes was revealed. However, the sensitivity between negative changes and phenological changes could be further explored, such as the range of negative changes that will not trigger changes in phenology and the threshold of negative changes that will affect changes in phenology. Third, when comparing the negative changes for different vegetation types, this paper assumed that the soil profile and depth were homogeneous. The influences of different soil depths and slopes on the drought resistance of vegetation will be considered in future work. Further investigations on this issue should thus be carried out in the future.

Conclusions
On the basis of the changes in vegetation growth in Jilin Province, China, during 2000-2018, a method to identify ecologically vulnerable areas was developed based on detecting abrupt changes in NDVI time series. Then, the response of phenology to negative changes was explored by analyzing the phenological indictors and the NDVI trend. The following conclusions were reached: (1) In this paper, an ecologically fragile zone identification framework based on the breakpoint detection of the BFAST NDVI time series was proposed. By identifying ecologically vulnerable areas to detect the number of negative changes and the magnitude of negative changes in vegetation growth over many years, we fully considered the long-term stability of vegetation growth and sensitivity to specific disturbances. This method can accurately reflect the long-term and short-term changes in vegetation growth and has better applicability in semi-arid regions. (2) During the past 19 years, northwest Jilin Province located in the semi-arid area was identified as ecologically vulnerable, primarily with low and slight vulnerability, where rainfall is scarce. Moreover, artificial irrigation cannot meet the needs of vegetation growth. Moisture is the main limiting factor leading to the vulnerability of the region's ecological environment. (3) Compared to other vegetation types with dense coverage, 60% of the area had a number of negative changes greater than 1, a much larger area than that of other vegetation types (less than 50%). The magnitude of change in sparse vegetation was −0.429, which is the lowest among the vegetation types. This shows that sparse vegetation is more susceptible to drought. Therefore, increasing the vegetation coverage or changing to more stable vegetation types can reduce the fragility in ecologically vulnerable areas. (4) For vegetation types with dense coverage, the impact of negative changes on vegetation phenology shows long-term effects. The negative changes in the NDVI trends of various types of vegetation led to a fluctuation range of the integral value from −0.06 to −6.9 and a phenological period length from −0.3 to −5.4; the peak value varied from −0.11 to 0.12. Negative change had a significant effect on the cumulative values of the growth phase, such as the relative amount of vegetation biomass and the length of the growing period, but less of an effect on the instantaneous value of the peak. Detecting changes in the growth phase or the integral value could be used to predict whether the vegetation growth experiences a negative change.