Detecting Patterns of Vegetation Gradual Changes (2001–2017) in Shiyang River Basin, Based on a Novel Framework

: A lot of timeseries satellite products have been well documented in exploring changes in ecosystems. However, algorithms allowing for measuring the directions, magnitudes, and timing of vegetation change, evaluating the major driving factors, and eventually predicting the future trends are still insu ﬃ cient. A novel framework focusing on addressing this problem was proposed in this study according to the temporal trajectory of Normalized Di ﬀ erence Vegetation Index (NDVI) timeseries of Moderate Resolution Imaging Spectroradiometer (MODIS). It divided the inter-annual changes in vegetation into four patterns: linear, exponential, logarithmic, and logistic. All the three non-linear patterns were di ﬀ erentiated automatically by ﬁtting a logistic function with prolonged NDVI timeseries. Finally, features of vegetation changes including where, when and how, were evaluated by the parameters in the logistic function. Our results showed that 87.39% of vegetation covered areas (maximum mean growing season NDVI in the 17 years not less than 0.2) in the Shiyng River basin experienced signiﬁcant changes during 2001–2017. The linear pattern, exponential pattern, logarithmic pattern, and logistic pattern accounted for 36.53%, 20.16%, 15.42%, and 15.27%, respectively. Increasing trends were dominant in all the patterns. The spatial distribution in both the patterns and the transition years at which vegetation gains / losses began or ended is of high consistency. The main years of transition for the exponential increasing pattern, the logarithmic increasing pattern, and the logarithmic increasing pattern were 2008–2011, 2003–2004, and 2009–2010, respectively. The period of 2006–2008 was the foremost period that NDVIs started to decline in Liangzhou Oasis and Minqin Oasis where almost all the decreasing patterns were concentrated. Potential disturbances of vegetation gradual changes in the basin are refer to as urbanization, expansion or reduction of agricultural oases, as well as measures in ecological projects, such as greenhouses building, a ﬀ orestation, grazing prohibition, etc.


Introduction
Satellite remote sensing has long been a technique of repeated temporal sampling on the earth's surface. Currently, a large body of sensors provides continuous observations at various spatial, temporal and spectral scales over decades. Compared with change-detection based on bi/multi-date imagery which focuses on comparing just a few scenes [1], timeseries products cover landscape change processes during the whole observation period, change-detections based on timeseries are scene independent and possess significant potential for capturing more subtle changes caused by climatic

Study Area
The study area is one of the three endorheic river basins in Gansu Province, northwestern China (36°45′-39°27′ N, 101°08′-103°50′ E) ( Figure 1). Owing to the continuous uplift of the Qinghai-Tibet Plateau since the Pliocene, the Shiyang River Basin became a typical temperate continental climate characterized by long-cold winter, rare and irregular precipitation, and high evaporation. The annual mean temperature is around 7.7 ℃. Annual precipitation shows a seasonal distribution. About 90% of the precipitation takes place in the period from May to October. Precipitation reduces gradually from the south to the north. The altitudes in the basin range from 1247 m to 4822 m above mean sea level. Both the climate and the landscapes in the basin show obvious zonality due to the altitudinal gradient.  The southern part is the Qilian Mountains, which has a colder semi-arid to semi-humid climate with annual precipitation 300-600 mm [26]. Vegetation covers alternate regularly from desert steppe (1800-2300 m), scrub (2300-2800 m), forest (2500-3200 m), meadow (3200-4200 m) to glacier summit. Picea crassifolia and Sabina przewalskii are the representative plants in the forest region [27]. The dominant species of shrub are Potentilla fruticose, Caragana jubata, and Salix gilashanica. The main herbs are Cares atrata, Polygonum viviparum, and Carex lanceolate. Vegetation in the Qilian Mountains plays a critical role in water conservation and runoff formation. However, a lot of them were reclaimed for crop production in the last decades of the 20st century. The northern part is dominant by typical temperate continental arid climate with an annual precipitation less than 300 mm while annual potential evaporation more than 2000 mm. The annual precipitation is less than 50 mm in the northernmost places that are surrounded by Badain Jaran Desert and Tengger Desert. Drought-resistant, salt-resistant shrub and perennial sand-loving herbaceous plants, such as Nitraria sibirica, Salsola passerine, Eriocaulon truncatum Ham, Reaumuria soongorica Pall. maxim, Peganum harmala, Artemisia arenaria, and Agriophyllum squarrosum are the native vegetation in this region [28]. Wetlands along rivers (or floodplain), were reclaimed for crop productions as early as Han Dynasty (2000 years ago) [29]. Currently, they are agricultural oases with high productivities and well-developed irrigation networks consisted of rivers, reservoirs, canals, and wells. Spring wheat, barley, corn, and cotton are the staple crops. The Shiyang River originates from the Qilian Mountains and flows northeast to terminal Qingtu Lake and Baiting Lake. It is the main surface water for vegetation growth in the basin. Most of the water in the river is used for agricultural irrigation [28]. June to September is the flood season of the river, which is also the season for crop growth.
Human activities had caused severe natural vegetation degradations during the period 1950-2000 [28,30]. Since the beginning of the 21st Century, the Chinese government has invested a lot to restore the ecological ecosystem. A series of water reallocation projects were implemented to alleviate the shortage of water resources or tradeoff available water among different reaches, e.g., "Jingdian Water Diversion Project", "Closing motor-Well and Reducing Cultivations", "Key Governance Planning for the Shiyang River Basin" (KGPSRB) [31]. Other measures, such as grazing prohibition, rotational grazing, afforestation and returning farmland to grassland or forest, also brought positive outcomes [32]. In addition, the developments of the economy and society have greatly altered the styles of utilizations of nature resources. Vegetation dynamics in the Shiyang River Basin were remarkable over the past years, which makes it an ideal place to carry out the study.

Data Resources
MODIS, on board earth observation system terra satellite of the National Oceanographic and Atmospheric Administration (NOAA), crosses the dayside equator at 10:30 a.m. local time [33]. The 250 m 16-day composite MOD13Q1 (Collection 6) datasets were downloaded from the online Data Pool at NASA Land Processes Distributed Active Archive Center (LPDAAC) [34]. The datasets were implemented several major improvements due to a new calibration approach and were released in February 2015 [33]. The datasets contain a variable number of "Science Data Sets" (SDS) that include 16-day NDVI, EVI, and quality assessment (QA) [35].
MOD13Q1 datasets were available for the period from April 2000 to December 2017 in our study. Mosaic, resampling, and reprojection to Universal Transverse Mercator (zone 47) were done firstly by using the MODIS Reprojection Tool (MRT). Then, all the pixels with NDVI values greater than zero (waterbody, snow, glacier or bare land) and the VI usefulness score in QA band greater than four, were masked out. We finally obtained all vegetation covered pixels with good observations. Lastly, a simple linear interpolation method was done to fulfill the missing values in the NDVI timeseries for each pixel.
Monthly NDVI values from January 2001 to December 2017 were produced by a maximum synthesis of the two-phase data. We calculated the mean value from May to September for each year and considered it as annual mean growing season NDVI which represented the overall state of vegetation growth expecting to reduce influences of anomalous events, e.g., rainstorm or flood. Table 1. Patterns of vegetation gradual changes are discussed in our study. Features of each pattern, including trend, transition years, duration and forecast of vegetation change are shown below. Options with the symbol ♦ will be excluded for discussion in our study, while the others will be considered for analyses in the study. (6) etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our dy are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively (Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our y are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively (Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our y are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the y stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively ( Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent rithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the n fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our dy are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively ( Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our dy are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively ( Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our dy are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively ( Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Table 1. Patterns of vegetation gradual changes are discussed in our study. Features of each pattern, including trend, transition years, duration and forecast of vegetation change are shown below. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the ears not less than 0.2 will be discussed in our study.

ethods
The Overview of Our Framework Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our dy are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively ( Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Table 1. Patterns of vegetation gradual changes are discussed in our study. Features of each pattern, including trend, transition years, duration and forecast of vegetation change are shown below. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified etation or desert vegetation is deeply influenced by fluctuations from precipitation in arid and i-arid regions, only pixels with maximum annual mean growing season NDVI (NDVImax) in the years not less than 0.2 will be discussed in our study.

. The Overview of Our Framework
Vegetation will keep in a stable state or slowly changing over time. Once disturbed, vegetation l change rapidly and shift to another state until a higher or lower equilibrium is reached [9]. refore, the processes of vegetation gradual changes within a certain period (2001-2017) in our dy are divided into four patterns based on the time and durations of the disturbance event. They ongoing change, stable in the early stage and keep changing later, continuously changing in the ly stage and then stable, change from equilibrium to a new equilibrium, which are named linear tern, exponential pattern, logarithmic pattern, and logistic pattern, respectively ( Table 1). All terns are further binned into two types, consisting of either positive or negative trends. Other nonnotonic patterns, such as the pattern like increase and then decrease (or decrease, then increase), not discussed in the framework. Table 1. Patterns of vegetation gradual changes are discussed in our study. Features of each pattern, including trend, transition years, duration and forecast of vegetation change are shown below. Options with the symbol ◊ will be excluded for discussion in our study, while the others will be considered for analyses in the study. Generally, each parameter in the logistical function has a ready ecological interpretation. An Se curve of logistic function could seem to be a combination of exponential curve and succedent arithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the in fitting function in our framework. In order to fit the patterns using a logistic function, modified Generally, each parameter in the logistical function has a ready ecological interpretation. An S-type curve of logistic function could seem to be a combination of exponential curve and succedent logarithmic curve. Therefore, a logistic mathematical model with four parameters is chosen as the main fitting function in our framework. In order to fit the patterns using a logistic function, modified operations

Noise Smooth by Using Moving Average
The smooth treatment is expected to reduce the influences of likely outliers in timeseries which may be caused by agriculture planting adjustments, inter-annual fluctuations of precipitation or other factors. A moving average (MA) is the commonly used smooth method, which could eliminate the seasonal and individual irregular changes, and highlight the long-term trend in timeseries. Specifically, for a timeseries x with n observations (n = 17 in our study), a single moving average in the forward and backward directions is calculated at every time point . Also, Then, a new timeseries y for each pixel is obtained.

Modeling All Non-Linear Patterns Using a Logistic Function
We model all non-linear patterns using a logistic function based on two core premises: (1) If vegetation growth is interfered, it will shift to an alternative state of rapid change. There must be a time p that the change rate in vegetation is maximum. The time p is consistent with the timepoint in the temporal profile of NDVI timeseries at which the slope value is maximum. (2) In order to present non-linear patterns, the smoothed NDVI timeseries y is prolonged by centering on p before applying to the logistic regression. The prolonged elements are filled by the first or last item of timeseries y. The concrete procedures are presented as follows: (1) Identifying the Timepoint P

Noise Smooth by Using Moving Average
The smooth treatment is expected to reduce the influences of likely outliers in timeseries which may be caused by agriculture planting adjustments, inter-annual fluctuations of precipitation or other factors. A moving average (MA) is the commonly used smooth method, which could eliminate the seasonal and individual irregular changes, and highlight the long-term trend in timeseries. Specifically, for a timeseries x with n observations (n = 17 in our study), a single moving average in the forward and backward directions is calculated at every time point i.
Then, a new timeseries y for each pixel is obtained.

Modeling All Non-Linear Patterns Using a Logistic Function
We model all non-linear patterns using a logistic function based on two core premises: (1) If vegetation growth is interfered, it will shift to an alternative state of rapid change. There must be a time p that the change rate in vegetation is maximum. The time p is consistent with the timepoint in the temporal profile of NDVI timeseries at which the slope value is maximum. (2) In order to present non-linear patterns, the smoothed NDVI timeseries y is prolonged by centering on p before applying to the logistic regression. The prolonged elements are filled by the first or last item of timeseries y. The concrete procedures are presented as follows: (1) Identifying the Timepoint P The differences of the forward and backward directions in timeseries y are calculated at every timepoint as follows: The year corresponding to the maximum absolute ∆y i is the timepoint p.
(2) Prolonging Timeseries Y The length of prolonged timeseries is symmetric about p. Therefore, the length of the prolonged timeseries (L) for each pixel is calculated as follows: where if p is equal to 9, the prolonged timeseries k is the original timeseries y. If p is greater than 9, the first 17 elements in prolonged timeseries k are original timeseries y, and others are supplemented by the last element in timeseries y. If p is less than 9, the last 17 elements in k are identical with y, and others are supplemented by the first element in timeseries y. Finally, new timeseries k which contains the whole timeseries y is obtained.
(3) Modeling the Prolonged Timeseries K Using a Logistic Function We fit the prolonged timeseries k by a logistic function. The function is shown as follows: where t are the serial numbers in timeseries k which range from 1 to L, f(t) are the values in the prolonged NDVI timeseries k, parameter a represents the change magnitude of NDVI over a period, the symbol of b denotes the direction of vegetation change, c is the location where the fitting value is equal to (a+d)/2 and parameter d reveals the initial background NDVI value. The goodness-of-fitting is implemented by the standard F statistics test. Only the goodness-of-fitting of the part corresponding to the period 2001-2017, are taken into consideration in our study.
(4) Distinguishing Non-Linear Patterns According to the Number of Transition Years The curvature (K) for the logistic function is computed first using the following equation: where z = e b(x−c) . a, b, c, are the parameters in Equation (4).
The transition year refers to the time at which vegetation gain/loss begins or ends. The beginning or stopping of vegetation change is attributed to a disturbance event or reaching a steady state. We define the transition years as the timepoints at which the change rate of curvature (K') in the logistic curve exhibits maximum or minimum. More information could refer to literature of Zhang et al. [36]. Specifically, if b is lower than zero, the locations of the two positive extreme values of K' are the transition years, if b is greater than zero, the locations of two negative extreme values of K' are the transition years. It is worth noting that we only concern the number and locations of the transition years in the period corresponding to years from 2001 to 2017. All the non-linear patterns are then distinguished by the number of transition times. If the number is 2, we conclude that vegetation change follows an S-type curve, while if the number is 1, vegetation change process is considered as exponential pattern or logarithmic pattern. Both patterns could further be differentiated by comparing the transition year to the parameter c in Equation (5). If the transition time is prior to c, the pattern is exponential, otherwise, it is logarithmic.

Identifying Linear Pattern of Vegetation Gradual Change
Linear regression of NDVI timeseries y against time is implemented to the pixels which do not pass the F statistics test in logistic fitting and pixels without transition year. The formula as below: The parameters are estimated based on OLS: Where i is the ith year in timeseries y, s, and slope are the parameters in linear regression, slope denotes the changing trend in vegetation, n is the length of timeseries y, NDVI i is the NDVI value in the ith year. The goodness-of-fitting is also implemented by adopting the standard F-test. If pixels pass F-test, we assume that vegetation growth in the pixels experience linear changes. The remaining pixels which do not belong to the aforementioned patterns are defined as no-trend.

Method for Validation
Field validation of year-to-year changes is often not straightforward due to a lack of consistent observation spanning two decades at a large spatial scale [22]. Fortunately, vegetation dynamics are closely related to landscape transformations or modifications caused by disturbance events [37,38]. Therefore, pattern of vegetation gradual changes obtained from our framework was evaluated by exploring its sensitivity to landscape transformations/modifications. Site validation and regional validation are combined to assess the performance of our framework in characterizing the temporal trajectory of vegetation gradual changes. If the shape of the curve in each site/region was consistent well with the process of land cover change disturbed by the known fact, and the observed transition years in the curve was consistent with the timing of the known fact, we assumed that the framework is effective in characterizing vegetation gradual changes. The known facts which led to significant landscape transformations or modifications were obtained from field investigations conducted from July to October in 2018 ( Figure 1). Most of them were further confirmed by the almost yearly high-resolution images on Google Earth Pro (version 7.3) and reports in newspaper or existing literature. It is worth noting that both the sites and the regions for validation were selected randomly from the areas where the same patterns concentrated. They are shown in Figure 3.

Patterns of Vegetation Gradual Changes in the Shiyang River Basin
Patterns of vegetation gradual changes in the Shiyang River Basin are shown in Figure 3. 87.39% of vegetation (NDVI max ≥ 0.2) changes were significant which had passed the F test (P < 0.05). The most widespread pattern was logistic (36.53%) ( Table 2). Linear patterns made up 20.16% of the entire vegetation covered areas, of which 18.43% were increasing while 1.73% were decreasing. Exponential and logarithmic patterns involved 15.42% and 15.27%, respectively. Most of them possessed increasing trends, while only about 1.88% and 0.79% had trends with negative slopes, respectively.
For the pixels with significantly changes, 87.57% were positive while only 12.43% were negative. Greening (Types I-IV in Figure 3) was universal in the basin and dominant in all patterns (Table 2)

Patterns of significant change in vegetation
No-trend Linear Linear exponential exponential logarithmic logarithmic logistic logistic

Patterns of Significant Change in Vegetation
No-trend Linear Linear Exponential Exponential Logarithmic Logarithmic Logistic Logistic

Transition Years and Durations of Vegetation Change
There is one transition year at which vegetation began or stopped to change in the exponential pattern and the logarithmic pattern (Table 1). There are two transition years in the logistic pattern. The transition years of the exponential pattern and the logarithmic pattern, and the first transition year of the logistic pattern are shown in      2008 was the dominant period that NDVI began to decline, which accounted for 55.22% of the logistic decreased pixels. The NDVI decreasing trends stopped in 2009 or 2010. The duration of vegetation change for each pixel was calculated through subtracting the end year from the start year ( Figure 6). More than half (51.19%) of logistic patterns possessed a continuous change interval of two years ( Figure 6. b and histogram below b). The longer vegetation changes lasted, the smaller the propouously reclaimed agricultural oases [39].

Forecasting Vegetation Condition in the Shiyang River Basin
The forecast for vegetation changes in the near future are shown in Figure 7. Our forecasts are not predictions, but rather are estimates of vegetation gradual changes which should continue on the trajectories of the past years. Therefore, we forecasted vegetation changes in near future according to the shape of the curve of NDVI timeseries in Table 1. Except for a few (12.91%) of vegetation covered areas which have no significant trends, 44.56% of vegetation will be in a stable state of high-level equilibrium, 7.25% will be in a low-level equilibrium, 31.97% of vegetation in the basin will keep improving persistently along with the current directions, and NDVIs in 3.61% of vegetation covered areas will keep declining. Vegetation in the Shiyang River Basin is generally in or toward a benign condition. The ongoing greenings in the high altitudes of the Qilian Mountains and Yongchang County are widespread. The continuously declining NDVIs mainly concentrated in the rural areas of

Forecasting Vegetation Condition in the Shiyang River Basin
The forecast for vegetation changes in the near future are shown in Figure 7. Our forecasts are not predictions, but rather are estimates of vegetation gradual changes which should continue on the trajectories of the past years. Therefore, we forecasted vegetation changes in near future according to the shape of the curve of NDVI timeseries in Table 1. Except for a few (12.91%) of vegetation covered areas which have no significant trends, 44.56% of vegetation will be in a stable state of high-level equilibrium, 7.25% will be in a low-level equilibrium, 31.97% of vegetation in the basin will keep improving persistently along with the current directions, and NDVIs in 3.61% of vegetation covered areas will keep declining. Vegetation in the Shiyang River Basin is generally in or toward a benign condition. The ongoing greenings in the high altitudes of the Qilian Mountains and Yongchang County are widespread. The continuously declining NDVIs mainly concentrated in the rural areas of Liangzhou oasis. More attention should be paid to find out reasons for the declined or declining NDVIs. Liangzhou oasis. More attention should be paid to find out reasons for the declined or declining NDVIs.

Results of Validation
The temporal profiles of the NDVI timeseries x, the smoothed timeseries y, and the prolonged timeseries k, the fitting functions and the locations of transition years for the eight selected reference sites are listed in Figure 8. The known fact of landscape transformation/modification in each site we obtained for validation is listed in Table 3. Specifically, the temporal trajectories in six sites of them were validated by subjectively comparing with the landscape transformations obtained from the images on Google Earth Pro. The Landscape modifications in the other two sites were identified by reports in the local newspaper. Our results indicate that the shape of the curve in every site coincide well with the processes of landscape changes obtained from field investigations, reports in newspaper, and higher resolution images on Google Earth Pro. Furthermore, a good agreement is also found between the detected transition years and the timing at which the landscape transformations or modifications took place in all the non-linear patterns.

Results of Validation
The temporal profiles of the NDVI timeseries x, the smoothed timeseries y, and the prolonged timeseries k, the fitting functions and the locations of transition years for the eight selected reference sites are listed in Figure 8. The known fact of landscape transformation/modification in each site we obtained for validation is listed in Table 3. Specifically, the temporal trajectories in six sites of them were validated by subjectively comparing with the landscape transformations obtained from the images on Google Earth Pro. The Landscape modifications in the other two sites were identified by reports in the local newspaper. Our results indicate that the shape of the curve in every site coincide well with the processes of landscape changes obtained from field investigations, reports in newspaper, and higher resolution images on Google Earth Pro. Furthermore, a good agreement is also found between the detected transition years and the timing at which the landscape transformations or modifications took place in all the non-linear patterns.  Figure. 3. The number of each subfigure corresponds to the number of reference sites in Figure 3, in which the pattern of each site detected by our framework is shown. The contents in each site include profiles of the original NDVI timeseries (blue solid line), the smoothed timeseries (green dotted line) and the fitting temporal profile (black line), the fitting function, and locations of the transition year (red word). Values in the horizontal axis in each subgraph represent the serial numbers in the prolonged timeseries (k) which ranges from 1 to L while that in the vertical axis are mean May to September NDVIs. The landscape transformation or modification for each site during the observation period are listed in Table 3.  Figure 3. The number of each subfigure corresponds to the number of reference sites in Figure 3, in which the pattern of each site detected by our framework is shown. The contents in each site include profiles of the original NDVI timeseries (blue solid line), the smoothed timeseries (green dotted line) and the fitting temporal profile (black line), the fitting function, and locations of the transition year (red word). Values in the horizontal axis in each subgraph represent the serial numbers in the prolonged timeseries (k) which ranges from 1 to L while that in the vertical axis are mean May to September NDVIs. The landscape transformation or modification for each site during the observation period are listed in Table 3. Table 3. Landscape transformations or modifications in the eight reference sites for validation ( Figure 8).

Site
Land Eight regions (approximately 5 km × 5 km) for validation are presented in Figure 3. Almost all the pixels in each selected region have the same pattern which is dominant in the region and the same or similar transition year except that in Region 5, 6, 7. The three regions are the places where linear decreasing pattern, exponential decreasing pattern and logarithmic decreasing pattern were concentrated. However, the exponential decreasing pattern in Region 6 is not dominant and the transition years in Region 6 vary greatly (Figure 4). Similarities are also found in Region 5 and Region 7. Therefore, only the other five regions are discussed in precise details (Figure 9). We first calculated the annual mean NDVI of each region. The pixels involved in the calculation in each region must have the same transition years as that listed in the table of Figure 9. Then, the pattern and the transition years of the region were also detected by adopting our framework.
We found that both the patterns and the transition years detected on the regional scale are coincident with that detected on the pixel scale. Therefore, the accuracy of our framework was further assessed on a regional scale. "Three-North Shelter Forest Program (Fourth and Fifth periods)", "Natural Forest Protection Project" and "Returning Farmland to Forest Project" were carried out successively in the Qilian Mountains before 2001. Tianzhu County has been committing to those projects over the past decades [42]. Logarithmic increasing patterns were found in the places where vegetation growths had been recovered completely and reached a high coverage [43]. Since the lower altitudinal limit of the forest belt in 1949 was 1900 m [44], linear increasing trends were common in most of the forest regions, like Region 1. Region 2 was reclaimed as large-scale farmland in 2014 or 2015, which we confirmed according to the yearly high-resolution images on Google Earth Pro. Agriculture oases are irrigated regularly, an increasing exponential pattern was detected in Region 2. The implementations of grassland restoration projects, such as "Natural Grassland Restoration and Construction Project" (2001)(2002)(2003)(2004)(2005) and "Returning Pasture to Nature Grassland" in Sunan County, Zhangye Municipality were earlier than other places (about in 2010) in the Qilian Mountains [45]. The grasslands recovered from overgrazing and reached a high stable state in 2003. Therefore, vegetation growth in Region 3 experienced a logarithmic increasing process with an earlier transition year (2003). Pastures in Region 4 were returned to nature grasslands in 2010. Vegetation recovered from overgrazing quickly [41] and logistic pattern was detected in Region 4. KGPSRB, invested largely by Chinese government was carried out in years of 2006-2011 [46]. In order to increase runoff and reduce groundwater extraction in the downstream area, the traditional crop plantings were replaced by protected agriculture (sunglass/plastic greening house) and large areas of croplands in Minqin Oasis were abandoned. Region

Strengths and Limitations of the Framework
Timeseries products present different temporal scales components, such as seasonal variations, long-term and short-term fluctuations [2]. Their high temporal resolution offers the opportunity to define time profiles of vegetation dynamics. Although many methods have been developed to explore vegetation dynamics at a global or regional scale, quantitative methods that could automatically address where, when and how vegetation change are still lacking [24]. Trajectory-based change detection is an inevitable choice to do this work by constructing a 'curve' or profile of full temporal records for each pixel. Nowadays, the method has been successfully applied in forestrelated change analyses [47,48], land use classification [49,50], and cropland variations [51].

Pattern
The transition year

Strengths and Limitations of the Framework
Timeseries products present different temporal scales components, such as seasonal variations, long-term and short-term fluctuations [2]. Their high temporal resolution offers the opportunity to define time profiles of vegetation dynamics. Although many methods have been developed to explore vegetation dynamics at a global or regional scale, quantitative methods that could automatically address where, when and how vegetation change are still lacking [24]. Trajectory-based change detection is an inevitable choice to do this work by constructing a 'curve' or profile of full temporal records for each pixel. Nowadays, the method has been successfully applied in forest-related change analyses [47,48], land use classification [49,50], and cropland variations [51].
The framework we proposed in our study is a trajectory-based change detection method. It has several desirable properties compared with other fashionable methods. Firstly, unlike the widely used trend analyses which hypothesize that ecosystems/vegetation always change linearly in a direction, our framework could detect the non-linear change patterns in vegetation. All patterns in the framework can be interpreted from unambiguous biophysical standpoints. Secondly, unlike DBEST and BFAST [18], our framework automatically generalizes processes of vegetation gradual changes into different patterns without setting any parameters. Furthermore, the patterns between different pixels or areas could be compared easily. Although DBEST and BFAST could also track vegetation gradual change by a piecewise linear model on the pixel scale, the number of breakpoints must be set by users when the method was extended to the regional scale. For example, de Jong [8] detected the short-term greening or browning periods within long-term timeseries by dividing the trend component in BFAST into four segments which were separated by three breakpoints. Then, the duration of the significant greening and browning segments, the magnitude of change in NDVI, and the abrupt changes could be compared between pixels subsequently. Lastly, trajectory-based change detection only functions when the process of vegetation change matches well with the hypothesized trajectory. Temporal trajectories of vegetation changes in different pixels may be different. Therefore, many mathematical models are needed to fit the trajectories. However, just a logistic model is used to simulate all the non-linear monotonic patterns in our framework.
Drawbacks of our framework are obvious. Firstly, only four monotonic patterns and eight types are presented. All the patterns are associated with only one disturbance event within the timescales of the observation period ( Figure 4). Generally, changes in vegetation are complicated and difficult to be generalized by a single function. Therefore, patterns that are composed of the general categories in our framework, or patterns that vegetation gradual changes were stalled or reverse by abrupt events, were not discussed in this study. New models that could detect both the monotonic and the non-monotonic patterns are needed to be induced in future. Secondly, the smoothed treatment implemented on the NDVI timeseries makes the framework unsuitable to detect subtle changes or abrupt changes [52]. However, it makes the framework unspecific to the MODIS timeseries products which do not suffer from the problem of data discontinuity. The data discontinuities which may be caused by changes in different sensors, orbital drift in the satellite overpass times, variations in sun-sensor-viewer geometry or differences in the atmospheric conditions (i.e., water vapour content, aerosols) [53] are common in other timeseries products, such as datasets from AVHRR, Landsat and VEGETATION.

Factors that Influence the Detected Patterns by Our Framework
Patterns of vegetation gradual changes are closely related to the regimes of the disturbances (intensity, time and duration) in our study. Our framework was designed based on the fact that the temporal trajectory of vegetation will alter after a major disturbance event. Therefore, the patterns of vegetation change largely depend on the regimes of the disturbance event. If the disturbance event priors to the observation period, vegetation growth may experience a logarithmic pattern or a linear pattern during the detected period. The exponential pattern or logistic pattern may be detected if the disturbance event is introduced in the midst of the observation period. The length of timeseries also affects the detected patterns. On one hand, a longer observation period means a more complicated pattern of vegetation dynamics, which our framework is incompetent to detect. Moreover, the longer the timeseries is, the more likely that this effect conceals actual short-term trends [8] and improves the power of the linear tread analysis [54]. On the other hand, if the observation period is short, or the duration of the disturbance event is long enough, vegetation changing condition may persist throughout the whole observation period. A linear pattern will also be detected. Finally, we found that the patterns of vegetation gradual changes are also influenced by the pre-existing land cover because the timescales over that the landscapes are affected by or recover from over exploitation are different. Compared with the forest ecosystem, grassland is more sensitive to disturbances [55]. It won't take them a long time to reach a new equilibrium. Therefore, grasslands or sparse vegetation regions in the lower altitudes of the Qilian Mountains were mostly characterized by the logistic pattern while the forest areas in the higher altitudes of the Qilian Mountains are dominant by the linear patterns (Figure 3).

Validation Method
The basic problem in timeseries-based change detection is the assessment of the method [56]. A general lack of reliable temporal field-based datasets spanning the duration of the satellite time-series limits the quantitative evaluation of the change detection methods. Most of the studies have no field data to support their findings or have no reliable techniques to assess the statistical significance of the detection methods [2]. Field investigation is still the most commonly used method for verification. Fuller et al. [10] examined the results of linear trend analysis of NDVI images in light of harvest measurements conducted in Senegal's rangelands and croplands. Qiu et al. [24] recently assessed the processes (dynamic patterns) of vegetation dynamics obtained by combination analysis of vegetation trends and temporal similarity trajectory through focusing on the primary and easily observable landscape changes. In addition, some researchers resorted to "validating" the trend analysis through the use of regional opinion or by invoking somewhat obliquely-related ancillary data sets and publications [54]. Since no independent datasets were available for validation, Kennedy et al. [22] tested the accuracy of a new trajectory-based method by comparing the detected changes with the direct interpreter delineation obtained from the timeseries images itself. Wessels [54] treated the simulation approach as an alternative avenue of field investigations to test the sensitivity of trend analyses (linear or non-parametric methods). We collected the evidence of vegetation gradual changes in the Shiyang River Basin in many ways, including field investigations, exiting publications and almost yearly higher resolution images on Google Earth Pro. The evidence of landscape transformations or modifications on both the pixel scale and the regional scale were used to evaluate the effectiveness of our framework.

Potential Driving Factors of Vegetation Gradual Changes
Our results demonstrated that vegetation in the Shiyang River Basin experienced significant changes during 2001-2017, and greening was widespread throughout the basin, which agrees with the recent research in this region [57]. Our study further indicated that human activities, especially various ecological restoration projects could explain a large part of gradual changes in the Shiyang River Basin because there is a good agreement between the detected patterns of vegetation changes and human-induced landscape transformations or modifications. Our results are consistent with that of Guan et al. [57] who also confirmed that vegetation changes in the east of the Qilian Mountains and oasis areas in the Hexi region, which is the Shiyang River Basin were mainly influenced by human activities.
Human activities refer to the expansion or reduction of agriculture oases, urbanization, migration and various measures in ecological restoration policy. The basin has stepped into a period of rapid urbanization after entering the 21st century. A large part of farmlands around the cities and towns had been transformed into urban areas, which caused significant logistic decreasing patterns during 2001-2017. Since there are obvious spatial differences in time at which agricultural oases expended [39], different increasing patterns were found in the newly reclaimed oases. The logarithmic increasing patterns were found in the north part of the transition zones between Liangzhou oasis and Tengger Desert (Site g in Figure 8), and the southern edge of Tengger Desert where were continually reclaimed as oases since 1990s. The Logistic patterns were found in the middle of Yongchang County. KGPSRB invested largely by Chinese government was carried out in years of 2006-2011 [46]. It largely changed the land-cover conditions in the whole basin through increasing vegetation coverage in the upstream areas, saving the water resource in the midstream areas, and increasing the runoff discharge in the downstream areas. As an important water-saving measure, a lot of greenhouses in Liangzhou and Minqin oases replaced the traditional crop productions in the early period of KGPSRB (2006KGPSRB ( -2008. Meanwhile, a total of 2323 pumping wells in Minqin Oasis involved in KGPSRB were closed from 2006 to 2010 to reduce groundwater exploitation [58]. Farmlands were abandoned because there is no water to irrigate. Therefore, widespread decreasing logistic patterns concentrated in the Minqin Oasis and the interior of Liangzhou Oasis. Forest projects such as "Three-North Shelter Forest Program (Fourth and Fifth periods)", "Natural Forest Protection Project" and "Returning Farmland to Forest Project" were carried out successively in the high altitudes of the Qilian Mountains since 2000. As mention in Section 4.4, most areas in the forest regions kept linear increasing due to tree growths. Policies for grassland protection such as "Natural Grassland Restoration and Construction Project" (2001)(2002)(2003)(2004)(2005) and "Returning pasture to nature grassland" were implemented prior to the observation period in Sunan County. Therefore, the logarithmic pattern is dominant in the west of the Qilian Mountains while the logistic pattern widely distributed in other places in the mountains. A stable and enlarging water surface reappeared in Qintu Lake in 2011 [59] and vegetation coverage enhanced gradually in the north of Minqin Oasis due to intensive afforestation and recovered agriculture productions [46]. Finally, people involved in KGPSRB and "Down from Hill and resettle Plain" moved out from the ecological conservation areas in the higher altitudes (above 2500 m) of the eastern Qilian Mountains, or the ecological deteriorated areas batch by batch in recently years [60]. Agricultural lands were returned to forest lands or grasslands soon after they left [61]. According to our results, widespread exponential increasing trends in vegetation with a transition year later than 2011 concentrated in these places. Generally, vegetation in the Shiyang River Basin has been experiencing the restoration from the past land degradations due to intensive human disturbances [62][63][64].
It is not surprising that landscape transformations such as deforestation or afforestation, expansion or reduction of the agricultural oases and urbanization are the processes of conversion that are affected by human activities other than inter-annual climatic variability. However, the processes of landscape modifications may be heavily influenced by interannual climatic fluctuation [65]. The non-linear change patterns in our study were closely related to landscape transformations caused by human activities. The impact of human activities on vegetation is more intensive and direct than the effects of climate variations in the arid and semi-arid environment. For example, the slight fluctuations of the inter-annual NDVIs in the grassland in the Qilian Mountains (Site g in Figure 7) may be caused by inter-annual variations in precipitation. However, the marked increase in 2010 was caused by policies of ecological restoration. Human activities may lead to significant show-term trends in timeseries. However, Climatic variations, especially precipitation has been proved to be the dominant natural factor affecting vegetation dynamics in arid and semi-arid areas [66,67]. Precipitation in northwestern China shows an increasing trend [68] and global warming also encourages vegetation greening in the high altitudes of the Qilian Mountains [69]. We suspected that the linear increasing patterns of vegetation in the Shiyang River Basin, particularly that in the upstream areas were partly attributed to the increasing humid-warm climate condition [70]. The assumption may be in line with the previous research [71] which confirmed that the pattern of NDVI is associated with climate variability on decadal and shorter time scales along with direct human-induced landscape conversions. However, the linear pattern of vegetation change caused by climate fluctuations, when coupled with that of various human activities, makes it difficult to devise effective schemes to identify the pattern of vegetation gradual change. Furthermore, their respective contributions to this linear pattern are difficult to quantify. Many works refer to the role of climate in vegetation gradual changes in the Shiyang River Basin and the discrimination between the human-induced changes and naturally-induced changes are needed to be done in future.

Conclusions
We presented a novel framework based on temporal trajectory detection method to fully understand processes of vegetation gradual changes. A case study in the Shiyang River Basin, northwestern China (the period 2001-2017) using MODIS NDVI datasets was shown.
It's a flexible method and superior to other fashionable methods. Our framework generalizes the processes of vegetation gradual changes into four patterns named linear pattern, exponential pattern, logarithmic pattern and logistic pattern. They represent the processes that vegetation is ongoing change, stable in the early stage and continuously changing later, constantly changing in the early stage and then stable, from an old equilibrium to a new equilibrium, respectively. All the non-linear patterns in the framework were fitted and differentiated automatically by using a logistic function with prolonging the original NDVI timeseries. Our framework not only makes comparison between different pixels/areas easily, but also identifies the timing at which vegetation changes took place, forecasts vegetation changes in the near future according to the shape of the curve during the observation period, and eventually contributes to explore the disturbance events. We found in our study that patterns of vegetation gradual changes largely depend on the regimes (intensity, time and duration) of the disturbances, as well as the length of timeseries and pre-change land cover types. Our framework is more suitable for characterizing the short-term trends in vegetation gradual changes which may be caused by a major disturbance.
Our results indicated that 87.39% of the vegetation covered areas (NDVI max ≥ 0.2) in the Shiyang River Basin experienced significant changes during the period 2001-2017, and a large part of the vegetation change patterns detected by our framework was non-linear. Increasing trends were dominant in all patterns. Spatially, vegetation in the east, the west, the high altitudes, and the lower altitudes of the Qilian Mountains experienced an exponential, logarithmic, linear, and logistic greening trends, respectively. The patterns of vegetation improvements in Yongchang County, the transition zones between Liangzhou Oasis and deserts, and the downstream areas of the river also showed great differences. 2008-2011, 2003-2004, and 2009-2010 were the main transition years of the exponential increasing pattern, the logarithmic increasing pattern, and the logarithmic increasing pattern, respectively. Ecological restoration projects and the expansion of agricultural oases were the main driving factors for vegetation improvements. Almost all the decreasing patterns located in Liangzhou Oasis and Minqin Oasis. Most of them were logistic patterns and NDVIs started to decline in the period of 2006-2008. Urbanization, large-scale greenhouse buildings and the reduction of agricultural oases for ecological purpose were the main reasons for the decreasing NDVIs. Generally, vegetation or ecosystem in the Shiyang River Basin is in or toward a benign condition.