Temporal Information Extraction for Afforestation in the Middle Section of the Yarlung Zangbo River Using Time-Series Landsat Images Based on Google Earth Engine

: Afforestation is one of the most efﬁcient ways to control land desertiﬁcation in the middle section of the Yarlung Zangbo River (YZR) valley. However, the lack of a quantitative way to record the planting time of artiﬁcial forest (AF) constrains further management for these forests. The long-term archived Landsat images (including the Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI)) provide a good opportunity to capture the temporal change information about AF plantations. Under the condition that there would be an abrupt increasing trend in the normalized difference vegetation index (NDVI) time-series curve after afforestation, and this characteristic can be thought of as the indicator of the AF planting time. To extract the indicator, an algorithm based on the Google Earth Engine (GEE) for detecting this trend change point (TCP) on the maximum NDVI time series within the growing season (May to September) was proposed. In this algorithm, the time-series NDVI was initially smoothed and segmented into two subspaces. Then, a trend change indicator S diff was calculated with the difference between the ﬁtting slopes of the subspaces before and after each target point. A self-adaptive method was applied to the NDVI series to ﬁnd the right year with the maximum TCP, which is recorded as the AF planting time. Based on the proposed method, the AF planting time of the middle section of the YZR valley from 1988 to 2020 was derived. The detected afforestation temporal information was validated by 222 samples collected from the ﬁeld survey, with a Pearson correlation coefﬁcient of 0.93 and a root mean squared error (RMSE) of 2.95 years. Meanwhile, the area distribution of the AF planted each year has good temporal consistency with the implementation of the eco-reconstruction project. Overall, the study provides a good way to map AF planting times that is not only helpful for sustainable management of AF areas but also provides a basis for further research on the impact of afforestation on desertiﬁcation control.


Introduction
As the social and economic center of Tibet [1], the middle section of the Yarlung Zangbo River (YZR) valley faces serious desertification due to intensive human activities and the fragile natural environment [2,3]. According to the statistics in Liu, et al. [4], there is a total of 2324.43 km 2 of aeolian sand in the middle part of the YZR basin, especially on the north bank of the wide valleys. To control desertification, the local government has launched many ecological protection and construction projects since the 1980s, including planting artificial forest (AF), artificial shrub, and artificial grassland for sand stabilization [5].
Ecological construction provides a natural laboratory to explore the impact of afforestation on natural environment improvement. To some extent, AF is of great importance in carbon sequestration [6][7][8], desertification control [9], soil improvement [10][11][12], climate regulation [13,14], and biodiversity conservation [15,16]. To investigate the detailed impact of afforestation on the local eco-environment, it is necessary to compare the environmental indicators change before and after the AF plantation [17,18]. Therefore, the AF planting time is the basic information for related researches between AF and the local environment, especially under the background of significant climate change on the Tibetan Plateau [19,20].
However, the AF plantation regions are sparsely distributed and the planting activities, including renewal and maintenance, have continued for a long period. Most of the afforestation information has been recorded manually without full spatial or temporal details. Consequently, it is difficult to obtain the exact planting time for each AF patch. Although tree-ring measurement can reveal temporal information, this point-based measurement cannot provide comprehensive AF spatial-temporal distribution recordation. Currently, the lack of this temporal information greatly constrains the assessment of the effect of ecological projects and the improvement for future ecological conservation activity. The afforestation in the YZR is conducted via transplanting the cultivated seedlings to the sandy area to improve soil conditions and avoid desertification with the increase of vegetation coverage at the same time. Therefore, the AF planting time can be reflected by the time when the land cover type transfers from sand to AF.
In recent decades, the openly available long-term remote sensing image archives, such as the Landsat satellite image archive [21,22] and the MODIS satellite image archive [23][24][25], recorded long terms of land cover changes which make the AF planting time detection possible. In the time-series remote sensing images, the land cover change caused by afforestation mentioned above is represented by the trend change of the temporal vegetation index. Taking the time-series normalized difference vegetation index (NDVI) as an example, the time-series NDVI retained a flat trend during the sandy period, while an increasing trend after AF plantation. Therefore, the trend change point (TCP) in time-series NDVI can be used to represent the AFs planting time.
In addition, the rapid development of remote sensing technology has provided efficient methods for forest monitoring, including forest mapping [26], forest change [27], and forest health assessment [28]. Kennedy, et al. [29] developed an algorithm called Land-Trendr and proposed the concept of tracking time-series images to detect forest changes. Fensholt, et al. [30] applied NDVI data from different platforms (NOAA, Terra, and SPOT) to conduct temporal-sequence analysis for long-term vegetation study. Huang, et al. [31] developed an algorithm called Vegetation Change Tracker (VCT) to monitor forest changes automatically using the Landsat time series stack (LTSS). Ochtyra, et al. [32] developed an algorithm called TVCMA and adopted a mixed-method based on deviation (threshold) and trend (regression) to detect vegetation changes. These detection methods perform well in temporal data change detection.
However, many factors, such as atmospheric noise and radiometric correction residuals distort the actual surface reflectance, causing outlier existence in the time-series NDVI [33]. These outliers usually existed with smaller duration and greater change degrees, which is easy to confuse with the AF plantation change pattern. Therefore, a trend change detection method is needed to express the change condition around each probe point while eliminating the influence of outliers.
In this study, aiming to map the AF planting times in the middle section of the YZR valley, an algorithm was proposed under the support of the Google Earth Engine (GEE) platform that tracks the temporal trend of the vegetation index acquired from long-term Landsat images. The basis of this research is that an abrupt increasing trend appears in the vegetation index curve after afforestation compared with the original sandy surface, and the AF planting time is detected by determining the time of this increase. Meanwhile, the outliers were reduced by an adaptive window smooth and subspace segmentation performed on the NDVI time series. With the developed method, a map of the AF spatial distribution and the AF planting time of the middle section of the YZR valley is acquired.

Study Area
The middle section of the YZR valley is located in the south of the Qinghai-Tibet Plateau (Figure 1a). This region is characterized by widely distributed sandy land, and land desertification is still an important environmental problem under the impact of human activities and climate change [34]. Therefore, the floodplain and sandy bank in this region have been planted with plenty of AFs to control desertification. As shown in Figure  1b, the study area surrounded by the red polygon has many green patches according to the Landsat-8 true-color image in 2020. Pictures of typical AF sites representing the AFs planted at different times are shown in Figure 1. Picture (I) shows the widely distributed AF planted on the riverbed. Picture (II) shows the mature AF with high vegetation cover. Picture (Ш) shows a newly planted AF.

NDVI Image Stack
To effectively capture the AF planting information, Landsat images acquired from the Landsat surface reflectance (Landsat-SR) product provided in the GEE platform were used in this study. In GEE, the Landsat-SR images operation was controlled by the internet-accessible application programming interactive (API) or the web-based interactive de-velopment environment (IDE), and the images were stored in the Google server and processed by the powerful parallel computing capability [35]. However, there are frequent clouds in the growing season (1 May to 30 September) because of the monsoon climate. Therefore, the cloud-free pixels were identified by the quality assessment (QA) band of the Landsat SR and used to form a composite cloud-free maximum NDVI image each year, for the maximum composite can be used to remove the effects of cloud-cover, scanangle, and solar-zenith-angle [36], and reflect the true growth condition of the vegetation [37]. The QA band was derived from the C code based on the Function of Mask (CFMASK) [38] algorithm implemented on the Landsat images to mask the score of the cloud, cloud shadow, and snow et al.
The earliest observation of the study area dated back to 1988. For this long-term NDVI series, they were derived from different Landsat platforms. To balance the systematic biases among different Landsat sensors, the inter-calibration coefficients proposed by Roy, et al. [39] were applied to each band of OLI images to reduce the effect caused by the difference between the sensors TM, ETM+, and OLI.

Validation Samples
To validate the mapping results of the planting time, a field survey was conducted from 8 to 15 August 2020, and the key AF patches (222 samples) were recorded by a handheld GPS ( Figure 2). Then, the planting time of these samples was artificially determined from the time-series annual maximum NDVI in the growing season extracted from the GEE platform. Among these samples, three AF patches are clearly indicated by a signboard around them which mentioned the exact planting year and they are noted with green color in Figure 2.

Methods
The mapping of the AF planting time was conducted via the following two steps: (1) AF extraction and (2) planting time detection. In the first step, several spectral indices related to the key land cover types in this region were calculated firstly with the cloudfree image in 2020 and imported into a random forest classifier to extract AF regions. The indices include NDVI, normalized difference sand index (NDSI), new water index (NWI), and normalized difference building index (NDBI). Then, an adaptive TCP detection method was developed to determine the planting time of the AF pixels using the yearly maximum NDVI time series. Figure 3 shows the flowchart of this research.

Artificial Forest Region Extraction
A variety of supervised classification methods have been integrated into the GEE platform, and the random forest (RF) classifier was selected here because of its wide application and reliable performance [40,41]. This method constructs multiple decision trees that depend on random vector values of independent sampling and classifies data through aggregation and guidance of decision trees [42], which is suitable for the high data dimensionality and multi-collinearity classification [43].
To perform this classification, 1840 land cover samples were collected through the high-resolution image from Google Earth and field survey, and the selection of those samples is based on the following two rules: First, the land cover type of each sample can be easily recognized in the Google Earth image; second, the pixels around the targeted sample have a similar land cover pattern. These samples include six categories: AF, cropland, sand, impervious surface, shrub, and water and the percentage of each category to the total samples is 20.93%, 16.69%, 14.62%, 19.29%, 14.67%, and 13.80%, respectively. Based on the unique features (spectral features and spectral indices) of each land cover type collected from those samples, a classification model was constructed with the RF classifier. Through times of tests, the classification accuracy reached the highest when the number of the decision tree was set to 110, and the number of the decision tree in the RF classifier was set to this value. After the classification, threefold cross-validation was implemented to evaluate the accuracy of the classification. The total training samples were divided into three parts, one of them was taken as validation samples each time and the other two parts were taken as training samples. The average value of the three times of classification accuracy was taken as the accuracy of the threefold cross-validation.

Afforestation Time Mapping
As indicated before, there is a significant increase in vegetation coverage represented by NDVI after afforestation in the sandy area, and the AF planting time can be captured by the abrupt change pattern that appears in the time-series NDVI curve. Because the NDVI images are only available since 1988, those AFs planted before 1988 are not considered in this study according to a threshold NDVI suggested by Anyamba and Tucker [44]. Therefore, when the NDVI value of any AF pixel in 1988 was higher than 0.2, it is thought that the AF was planted before 1988. After excluding the AFs planted before 1988, the planting times of other AF pixels were detected via the following two steps.

Subspace Construction
The developed algorithm for trend change detection includes six categories: threshold, difference, segmentation, trajectory classification, statistical boundary, and regression [45]. To represent the continuous trend around the probe point, a subspace segmentation was first set up to divide the NDVI trajectory into collections. The concept of the subspace in the change detection is to compare the subspaces in the trajectory between past and present and measure their dissimilarity by using a parameter to indicate the difference distance between them [18]. Figure 4 shows the sliding subspace with Y as the target year and T as the width of the subspace (take the width of three years as an example). Therefore, the time-series NDVIs from Y-T to Y made up the subspace before Y, and the NDVIs from Y to Y+T made up the subspace after Y. Probe point Y slides along the trajectory and the subspaces are fixed on both sides of the probe point and also slide along the trajectory. In the subspace segmentation, the NDVI values obtained before and after the targeted year were needed for the study period 1988-2020. However, the NDVI values before 1988 or after 2020 are not available. Therefore, the NDVI values of 1998 and 2020 were replicated to the years before or after to meet the requirements of subspace construction. Taking 1988 as an example, the subspace after 1988 was available while the subspace before 1988 was nonexistent. Therefore, the NDVI value in 1988 was replicated to be the values of the previous years to reconstruct the subspace before 1988.

TCP Indicator Construction
After subspace segmentation, a linear fit was conducted for the NDVI values in both subspaces at each probe point to get the fitting line slopes. As shown in Figure 5a, the slopes of the left subspace (Sbef) and the right subspace (Saft) were calculated for probe point Y. Then, the difference between those two slopes (Sdiff) was used as an indicator to represent the variation of the curve; the indicator was inspired by Zuo, et al. [46]. Obviously, if there was an increasing trend at this point, Sdiff should have a positive value and vice versa. Meanwhile, the absolute value of Sdiff represents the degree of the abrupt changes, the larger it is the greater the degree of the abrupt change is. To demonstrate this fact, Figure 5b shows the original NDVI time-series curve in Figure 5a and the Sdiff curve derived with the subspace width of three years as an example. There is one year with the maximum Sdiff value in 2009 which is thought of as the AF planting time.

Adaptive TCP Detection
Ideally, the detection of the TCP can be simply determined according to the maximum value of Sdiff as shown in Figure 6. However, due to some uncertainties in the NDVI data caused by drought, plant diseases, insect pests, deforestation, inundation induced by temporary flood, and also the atmospheric effect on remote sensing data, there would be some outliers distorting the actual NDVI changing pattern. To reduce these influences and enhance the detection accuracy, a sliding window smooth method was applied to the NDVI trajectory before detection to smooth the abrupt changes induced by these uncertainties. Figure 6 shows a good example of outlier smoothing. Before smoothing, the actual TCP (2004) was not accurately detected according to the maximum value of Sdiff shown in Figure 6a. However, after smoothing (Figure 6b), the interference of the uncertainties in time-series NDVI was reduced and the expected TCP was accurately found. Although the sliding window smoothing method can reduce the impacts from the uncertainties in the NDVI trajectory, how to determine the optimal width of the sliding window is important to effectively remove the disturbance from outliers. Therefore, an iteration detection scheme was implemented by changing the subspace width and the smoothing window size adaptively. Firstly, the initial subspace width was set to 2 and its maximum value was set to 5. The sliding window size was initialized at 2 with a maximum of 6. For each pixel which is hard to determine the dominated maximum Sdiff value, its subspace width was enlarged first to check the difference between the first two peak values of Sdiff. If the difference was still very small (the second biggest Sdiff is larger than 2/3 of the first biggest Sdiff), the subspace was enlarged again. Once the subspace extended to its maximum size, the smoothing window size was extended subsequently and the Sdiff series was generated again at different subspace sizes from 2 to 5. The iteration was stopped when the proportion of the second Sdiff peak value to the first peak value was no more than two-thirds. Finally, the TCP was set to the year of the maximum Sdiff.

Artificial Forest Extraction Result
With the proposed classification method, the AF distribution of the study area was mapped as shown in Figure 7a. Threefold cross-validation was conducted for the classification results with an overall accuracy of 85.1% and a Kappa coefficient of 0.821. The main errors came from the misclassification among the vegetation surfaces likes AF, cropland, and shrubland. Especially for the newly planted forests, their poor vegetation coverage makes them similar to those of shrubs or crops. Clearly, the AFs are intensively distributed on the southern shore while the north shore has relatively few AFs. This difference was mainly attributed to the implementation of ecological engineering projects on the south shore in the early stage due to the widely distributed settlements. To better present the extraction accuracy, three plots were selected to show their original image and the extracted AF at each plot. Through visual comparison between the corresponding Landsat false-color images shown in Figure 7b-d and the AF extractions in Figure 7e-g, the extracted AFs have good spatial consistency with the real distribution.

NDVI Series and the TCP Detections of Typical AF Samples
With the inspection of the NDVI trajectories of typical AFs, the NDVI series can be generally separated into three key types: ordinary AF, AF planted before 1988, and temporally inundated AF. Ordinary AFs are usually planted on barren surfaces or sandy surfaces to prevent desertification. Therefore, the NDVI curve has an obvious increase after AF plantation, and the NDVI value sticks to a relatively low level before AF planation, as shown in Figure 8a. The second type belongs to the AF patches with very good vegetation cover which were planted before 1988 (the earliest Landsat observable time). As shown in Figure 8b, their NDVI value always remains at a relatively high level and the variation is not as large as that of the ordinary AF. Therefore, the AF pixels with the mean NDVI in the initial three years bigger than a certain threshold (0.2 in this study) were considered to be planted in or before 1988. The third type was mainly distributed in the sandy islands of the riverbed, which was temporally inundated. As shown in Figure 8c To prove the reliability of representing the AF planting time by the TCP in time-series NDVI, Figure 9 shows the false-color images before and after the detected year, the NDVI trajectory, and the Google Earth image in 2019 of five AF samples. In the composite image, vegetation was shown as red color while water is dark blue and barren surface is yellow or white. The AF pixels were marked with yellow rectangles. Obviously, the AF pixels were water surfaces or sandy surfaces in the year before the detected year. After AF plantation, the pixels turn to light red, and this change can be reflected by the TCP detection shown in the NDVI trajectories. Therefore, the TCP of the time-series NDVI could represent the planting time of the AF by comparing the land cover changes in the afforestation area. Figure 9. Comparison between the Landsat false-color images before (the first column) and after (the third column) the detected year of five sample regions. The yellow squares represent the AF pixels. The NDVI curves in the second column were extracted from the averaged NDVI of the AF pixels marked with yellow squares. The fourth column shows the latest Google Earth images of these five regions. (a)-(e) are five AF samples with different planting times.

Afforestation Time Mapping Result
The proposed detection method was applied to the whole study area, and the mapping result of AF planting times is displayed in Figure 10a. The temporal-spatial distribution of AFs presents as regional clustering, which obeys the rule that AF plantation was always conducted from one region to another. Overall, the planting time of the south bank of the YZR was generally earlier than that of the north bank because prevention and control of desertification projects were mainly conducted at the south bank due to its widely distributed settlements. Take the three subregions in Figure 7 as examples, the AFs in Figure 10b were planted on the sandy island and strongly influenced by inundation. These pixels presented a patchy distribution which partly confirms the fact that AF planting was conducted at different times in different sized areas. In contrast, although there were also patches with different planting times, the planting times displayed in Figure 10c, and d present clear temporal similarity for most of the AF pixels. In Figure 10c, the AF pixels close to the riverway (at the bottom) have an earlier planting time because of better water conditions than the pixels in the north. Similarly, the AF pixels in Figure 10d, in the south were first planted to protect the road, and then the AF plantation was conducted on the river shoal.

Validation with Field Samples
The detected AF planting time result was directly validated by the collected samples first. The result is shown in Figure 11. The density scatter plot shows good agreement between the detection results and the samples, with a Pearson correlation coefficient of 0.93 and a Root Mean Square Error (RMSE) of 2.95 years. Most of these points are well distributed around the 1:1 line. Meanwhile, the three AF field samples with true planting time also perform well in the consistency with the detected results. Two of them were consistent with the recorded planting time, while only one point has one year error. Meanwhile, there are also some points with big overestimation in the planting time and this may be attributed to the variation in vegetation cover induced by water inundation or deforestation. Figure 11. The density plot between the AF afforestation time of the samples extracted by manual identification and detected by the proposed method. The three green stars represent these samples with true value in planting time.

Temporal Consistency Analysis with Implementation of Ecological Projects
With the detection results, the temporal distribution of the AF planting time over the study area from 1988 to 2020 was investigated and shown in Figure 12. The observable time series can be divided into four phases according to the temporal pattern, which is highly consistent with the implementation of key ecological construction projects during this period.
In the early years (from 1988 to 1998), the afforestation area in each year always keeps to a low level, because the afforestation was mainly conducted voluntarily by civilians around the settlements and there is small input from local or central government. It can be referred from the report from this website http://www.seac.gov.cn/seac/xwzx/200303/1010522.shtml, accessed on 25 November 2021. From 1999 to 2007, the afforestation area has slight increases because of the project of returning cultivated land to forest (grassland). However, the government mainly focused on the afforestation along the roads and railways, which occupy a small area in this region.
In the following years (from 2008 to 2014), the afforestation area had a fast increase because of the conduction of the Plan for the Protection and Construction of Ecological Safety Barriers in Tibet, and the central government invested 8.35 billion yuan to protect and restore the key ecological function zones in the Tibetan Plateau (http://www.scio.gov.cn/zfbps/32832/Document/1633895/1633895.htm, accessed on 25 November 2021). In recent years (from 2015 to 2020), the project for the protection and construction of an eco-safety barrier (around 10.7 billion investments) assists the afforestation in the basins of the Yarlung Zangbo, Nujiang, Lhasa, Nyangchu, Yalong, and Shiquan Rivers. Therefore, the afforestation area continues to be at a high level in this region. Obviously, the temporal consistency between the temporal pattern of the AF planting time and the implementation of key ecological projects in this region confirms the reliability of the detection results.

Discussion
Many algorithms have been proven to perform well in forest change mapping. However, not all of them are suitable to detect the planting time of the AFs. The afforestation generated continuous trend changes with a smaller degree than that of the outliers and the purpose of this study focuses on extracting the AF planting time from the disturbance of outliers. In this study, trend, difference, and segmentation were combined to carry out continuous trend change detection, and the adaptive scale selection was imported to reduce the influence of outliers. In relative studies, Liu, et al. [47] mapped the afforestation and deforestation using the IFZ index proposed in VCT [31] which has proven to perform well in forest disturbance detection, and the planting time was detected by the threshold of the IFZ. However, because of high spatial-temporal heterogeneity and external factors, there may be some outliers with different change degrees, distorting the planting time detection. Consequently, the threshold is different due to different afforestation ground conditions, such as sandy land, grassland, and shrubs, which may be confused by the condition of cropland. Therefore, the results using thresholds judgment may be confused by different environmental conditions in large-scale spatial-temporal distribution planting time mapping, and this error can be avoided by focusing on the trend change characteristic in the time-series vegetation index.
However, it is difficult to avoid the disturbance of the outliers during change detection. Spatial heterogeneity results in different NDVI value ranges during afforestation, and the use of threshold is always controversial, such as the TVCMA taking thresholds as change detection elements. The trend in the change detection can represent an overall tendency of a phenomenon and reduce the impact of individual outliers [48], while the use of single-point data is not enough to express the continuous trend change characteristic and is susceptible to outliers. Meanwhile, for local analysis to extract the trend features from the NDVI trajectory, the subspace was used for segmentation. Moreover, the indicator Sdiff was introduced to indicate the continuous trend change around the probe point. For the continuous construction of ecological engineering, the TCP (AF planting time) can exist at any time in the time series and produce a different mean square error (MSE) during the total least square fitting in LandTrendr, which may be confused by the MSE of outliers. Thus, the adaptive selection of the window size and the subspace size was imported to reduce the influence of outliers. Finally, aiming at the characteristics of the AF temporal vegetation index, combining the three elements minimizes the impact of outliers while expressing the trend change characteristic.
The detection results obtained in this study proved that the proposed change detection is suitable for extracting the TCP in time-series vegetation index curve which is firmly connected to the AF planting time by finding the characteristic within dense image stacks. In contrast to land cover change analysis of images from certain time points, such as the land cover mapping with images every five years or just two images at different times for analysis [49,50], the use of time-series images from each year covered the changing information with higher temporal resolution. Meanwhile, the powerful computing capability of the GEE platform greatly improved the efficiency of dense image series analysis, there are similar investigations based on GEE in the field of agriculture [51], population [52], vegetation [53], and surface water [54].
With the above basic structures and platform, the adaptive segmentation and smooth in this method shows good accuracy in finding the right TCP by automatically selecting the suitable scale of the smoothing window and segmentation subspace for each pixel, which helps remove the impact of outliers. Because changes at different degrees require different scales of window smooth and subspace segmentation, the oversize of the smoothing window will miss the exact detection of the real TCP, while the smaller size cannot remove the disturbance of outliers. Therefore, the adaptive scale selection introduced in this study alleviates this impact. Taking a sample pixel as an example ( Figure  13), the Sdiff of an outlier in 2009 was smoothed after the width of the smoothing window increased and the TCP remain unchanged. Through a wider smoothing window, the change degree at a smaller scale was reduced and the smoothed curve provided a better description of the changes in the NDVI trajectory than the initial curve. In other words, the use of smoothing windows at different scales produces different smoothing effects, which means that a single scale of the smoothing window is not suitable for each pixel to remove the disturbance. The proposed TCP detection method shows high resistance to disturbances through adaptive changes in the window width and the length of the subspace. However, for some cases when there are some significant disturbances after AF planting, such as deforestation, the NDVI trajectory does not display a continuously increasing trend but is mixed with some declining periods or a sudden increase as shown in Figure 14, which generates a larger Sdiff during the growth stage than that of the planting time. The increasing stage after this disturbance distorts the TCP detection results and introduces errors in the results. Therefore, it is necessary to combine other quantitative indicators to constrain the detected time and stick to a low biomass level to reduce the impact caused by natural or human factors during the growing stage. In addition to the uncertainty related to the disturbances, the uncertainty associated with AF extraction is another factor that should be considered. The wide distribution of planting time and different growth stages of the AF make the AF information easily confused with the spectral pattern of the cropland and shrubland. Due to the unique texture of AF planting, texture extraction [55,56] is a widely used method for extracting AF areas in remote sensing. However, the spatial resolution of the Landsat image is not sufficiently high for AF texture analysis. Therefore, there is still space to improve the AF extraction accuracy. Inclusion of the texture information from higher resolution image data, such as the Unmanned Aerial Vehicle (UAV) based Remote Sensing images, and classification at the patch level will partly assist the further improvement of AF planting time mapping.

Conclusions
In this study, an adaptive detection algorithm for identifying the planting times of AFs was proposed based on the time-series NDVI values acquired from the Google Earth Engine platform. With this method, an AF planting time map was generated for the middle section of the YZR, the AF extraction from RF classification with an overall accuracy of 85.1%, and a Kappa coefficient of 0.821. The validation of the detected planting time with samples collected in this region indicates the good performance of the proposed method with a Pearson correlation coefficient of 0.93 and an RMSE of 2.95 years. The results showed that the AFs planted in the early years were mainly distributed in the southern shore of the YZR, and the AFs of recent decades occupied a large area especially the AFs planted from 2008 to 2014, and another peak existed in 2019, which presented good consistency with the implementation time of the eco-reconstruction projects in this region.
In general, the self-adaptive process in the proposed method reduced the impacts of disturbances induced by outliers effectively. Meanwhile, the confusion with the real planting time caused by different afforestation conditions was reduced by focusing on the trend change characteristic in time-series NDVI. This study provides an important insight into the spatial and temporal distribution of AF along the YZR valley, which is very important for further assessing ecological projects and evaluating their functions in controlling widespread desertification. According to this data, the environmental changes of AFs areas with different spatial-temporal distributions can be used to quantify the role of afforestation in resisting desertification.