Phenology-based cropland retirement remote sensing model: a case study in Yan’an, Loess Plateau, China

ABSTRACT Cropland retirement is a widespread phenomenon across the world. The conversion of inefficient cropland to forest or grassland is a policy aimed at restoring ecology, improving the environment, and promoting economic development. However, in most developing countries, the results of cropland retirement and land restoration are characterized by spatial fragmentation, and there are significant temporal differences as a result of poor agricultural intensification, human interference, and regional environmental differences. This substantially increases the difficulty of information extraction and reduces the extraction accuracy of remote sensing methods. In this paper, we developed a new phenology-based cropland retirement remote sensing (PCRRS) model to detect the extent and timing of cropland retirement. Considering the characteristic growth of crops, the normalized difference vegetation index (NDVI), at the start, middle, and end of the growth cycle, is the phenological metric to distinguish cropland from other vegetation types. In addition, the interannual variation of phenological metrics are significant after cropland retirement, which is the key to effectively identify retired cropland. High-resolution Google Earth images were used to verify the accuracy of the algorithm. The results suggested that the overall accuracy of our algorithm exceeded 85%, and was more suitable for sloping cropland. In comparison with other cropland retirement extraction methods, the PCRRS model had high sensitivity and stability. We found it was common for sloping cropland to be retired earlier, and we also identified the existing inter-planting phenomena between crops and shrubs in areas with gentle slopes. Overall, this study provided a basis for understanding the drivers of cropland retirement and evaluating their environmental effects.


Introduction
Cropland retirement is a major form of agricultural land-use change and affect water resources, carbon storage and biodiversity (Isbell et al. 2019). In the last century, large-scale cropland retirement has occurred in North America (Brown et al. 2005;Klooster 2003), Europe (Estel et al. 2015;Griffiths et al. 2013;Prishchepov et al. 2012b), and Central Asia (Dara et al. 2018;Kraemer et al. 2015). However, much less is known about the spatial and temporal distribution results of cropland retirement.
In recent decades, cropland retirement and reforestation have become widespread in developing countries (Díaz et al. 2011;Li et al. 2017), especially China . The ecological restoration "Grain for Green" (GFG) program enacted in 1999 covers vast tracts of China. Several benefits in the project implementation region have been reported, including increased vegetation coverage, control of water and soil erosion, and reduced sediment loads in water bodies Chansheng, He, and He 2016;Lü et al. 2017). In the conversion area, environmental consequences (Xu et al. 2019) such as ecological restoration, carbon storage (Deng, Liu, and Shangguan 2015;Scott et al. 2015), and water resources (Benayas 2007) depend on the time that has elapsed since cropland retirement. Few studies have focused on the spatiotemporal patterns of cropland retirement. The lack of spatially explicit data complicates evaluations of the drivers and environmental effects of cropland retirement.
Analysis of satellite imagery is the most convenient and efficient way to monitor changes in agricultural fields, such as cropland retirement (Duveiller and Defourny 2010;Zheng et al. 2015). Because agricultural land abandonment is defined as cropland going unused for a minimum of 2-5 years (Yin et al. 2020), previous studies that relied on satellite imagery at specific times (Bergen et al. 2012;Falcucci, Maiorano, and Boitani 2007) resulted in retired cropland being easily confused with fallow areas (Prishchepov et al. 2013). Land use datasets are as source of data about cropland-forest conversions Yin et al. 2018b); however, the process of dataset interpretation is time-consuming and labor-intensive, and mapping error has transitivity (Cao, Chen, and Yu 2009;Liu, Liu, and Kuang 2020). Recently, several timeseries methods that map the temporal trajectories of spectral indices to land-cover variation (Dara et al. 2018;Nguyen et al. 2018;Yin et al. 2018a), such as the normalized difference vegetation index (NDVI), have been used to monitor changes within specific land cover categories (Devries et al. 2015;Yan and Roy 2015). However, it is difficult for these methods to identify agricultural land conversions accurately due to insufficient spectral reflectance, which also varies greatly after changes in crop types and planting patterns on cultivated land (Pazúr et al. 2020;Yin et al. 2020). Clouds and haze present a challenge for data availability in optical time-series algorithms, and high-quality training data for identifying cropland retirement maps are also difficult to obtain (Gomez, White, and Wulder 2016;Xu et al. 2018).
Differences in phenology are key to separating retired cropland from cropland still in use (Alcantara et al. 2012). Multiple images covering different periods of the growing season, especially in the early and late stages, are typically used to obtain the highest classification accuracy in any land cover category (Homer et al. 2004;Pax-Lenney and Woodcock 1997). Because of the large particles investigated and low complexity of land use types, these studies have not been affected by the mixed pixel problem (Yang et al. 2019). Cultivated land and retired cropland are spatially dispersed and easily disturbed by mixed pixels due to the low agricultural productivity and scarce per capita cultivated land area in most developing countries. Therefore, coarse resolution remote sensing images may not be able to achieve the extraction accuracy required to identify retired cropland (Yin et al. 2018b). In addition, studies based on phenological changes have significant limitations, because phenology changes with the climate (Beurs and Henebry 2004).
Considering the complexity of cropland retirement recognition and inapplicability of existing mainstream methods, this study developed an efficient method to identify and extract the spatiotemporal distribution of cropland retirement, by combining the phenological metrics of the start, middle, and end of the crop growth cycle. With the rapid growth of crops during the growing period and rapid decline of the NDVI value in the harvesting period ), a reasonable threshold of the temporal change rate of the NDVI value can effectively eliminate the interference of the mixed-pixel problem on model accuracy. More importantly, differences in phenological metrics caused by climate change are slight due to the low vegetation coverage post-abandonment (Ma, Churkina, and Trusilova 2012). We simulated the distribution of cropland retirement using the relationship between the interannual variation of crop phenological metrics and the proportion of cropland conversion area. Our method can effectively avoid interference by the phenological changes caused by climate change. With an increase in the time density of mapping, our model can effectively prevent any impact of fallow fields on the mapping accuracy of cropland retirement.
We selected a typical semi-arid region in China as the study area. The area has experienced large-scale cropland retirement in the past 20 years. The low agricultural intensification, single economic source, and human disturbance has led to irregularities in the size and shape of retired cropland areas. The study area of Yan'an is located on the Loess Plateau of China and is a focal region for the implementation of the GFG program. Specifically, we developed a high-precision remote-sensing extraction model for cropland retirement in areas with severely mixed pixels and revealed the temporal and spatial patterns of cropland retirement on the Loess Plateau of China.

Study area
The Yan'an region is located in the center of the Loess Plateau in China, and a total land area of 37,000 km 2 ( Figure 1). The terrain is largely hills and gullies. The semi-arid climate has a low average annual precipitation of 250 mm, and there is significant intra-year variability most rainfall occurs between June and September. Due to the steep terrain and irregular rainfall, long-term and extensive human activities have resulted in serious soil erosion in Yan'an area . In order to promote ecological restoration, Yan'an was the initial site adopted by the Chinese government to implement the GFG program in 1999 (Chen and Chen 2006). According to current statistics, approximately 7,183.3 km 2 of cropland has been afforested (Liu et al. 2018).
Maize and wheat are the main food crops in Yan'an area. There are significant spatial differences in the planting patterns of maize and wheat. The terrain in the northern part of Yan'an is dominated by hills and ravines, and the agricultural land is mainly grown with maize. The terrain in the south is dominated by plains and ravines, and the agricultural land adopts the rotation method of summer maize and winter wheat. The growing period for the summer maize typically spans from mid-May to early-October, and that for winter wheat spans from mid-October to early-May in the following year. We mainly estimated the results of maize crop retirement. Generally, the results of summer maize crop retirement could almost reflect the overall cropland retirement results of the cropland adopts the rotation method in this study area.

Data
The high temporal resolution of daily coverage provided by Moderate Resolution Imaging Spectroradiometer (MODIS) satellites allows detecting seasonal vegetation phenology dynamics. MOD13Q1 16-day synthetic data can effectively reduce the interference of cloud and rain (Ganguly et al. 2010). We used a time-series of NDVI data from MOD13Q1 (Collection 6) in 2001 to construct intra-annual timeseries NDVI curves for various surface vegetation types. Then, the phenology metrics (such as the start, middle and end of the growth period) of different surface vegetation types were extrapolated based on the NDVI curves, which were used to indicate the correct Landsat data to collect.
We collected four Landsat footprints (WRS-2 paths/ rows 128/34, 127/34, 127/35, and 126/35) covering the Yan'an area for all growing seasons between 2000 and 2020 from the United States Geological Survey (USGS) and obtained reflectance data. We applied a cloud cover threshold of < 20% and obtained terrain-corrected image. Unfortunately, Landsat images showing different land-use stages and phenological metrics were not always available as a result of the influences of clouds. Therefore, based on the modified neighborhood similar pixel interpolator theory (Yang et al. 2020), we used healthy pixels to substitute any pixels with cloud pollution. These replaced pixels were acquired on the same or adjacent dates in nearby years (within a 3-year window). Considering the differences in spectral reflectance of the sensors, we applied previously reported parameters to normalize Operational Land Imager (Landsat8) reflectance to Landsat Thematic Mapper (Landsat5) and Enhanced Thematic Mapper Plus (Landsat7) values .
We calculated NDVI based on the red band and near-infrared band. The ALOS Global Digital Surface Model "ALOS World 3D -30 m" (AW3D30) in 2015 was resampled to obtain the slope data of the study area. Here, 0-15° is defined as a gentle slope, and above 15° is defined as a steep slope. We used the global land cover (GLC) data (obtained from http://data.ess.tsinghua. edu.cn) to identify impervious areas and water bodies. The GLC and MCD12Q1 data were used to provide a contrast experiment with our method. Additionally, the support vector machines classifier (SVM) was used to derive land cover maps to capture the transfer of agricultural land use types.
According to their respective results of retired cropland, confusion matrices were estimated based on Google Earth historical images samples. Figure 2 shows the detailed research process. In the first step, we created an annual Landsat NDVI dataset based on crop phenology metrics obtained from the timeseries curves of MODIS NDVI. Second, we created a cropland extraction index (CEI) in Section 3.2.1. The construction of the NDVI dataset depends on the type of crops in the study area. This paper only calculated the CEI of the maize crop, we used the results of maize crop retirement to approximate reflected the overall retired results of cropland using the rotation method in this study area. And section 3.2.2 explains the calculation of CEI as a variable of the phenology-based cropland retirement remote sensing model (PCRRS) for deriving the extent of cropland retirement. Considering the effect of climate change on differences in phenological metrics, we therefore evaluated the consistency of annual CEI values by visual interpretation and analytical samples obtained from Google Earth high-resolution imagery. Analytical samples were stable in terms of their land cover type and range; i.e. they were always cropland or forest from 2000 to 2020. Third, the thresholds of CEI and PCRRS were analyzed using Google Earth highresolution imagery due to their critical role in improving  model accuracy. Finally, we distinguished retired cropland and stable cropland based on the PCRRS model and evaluated the mapping accuracy. And the cropland retirement results obtained by GLC product, SVM model and PCRRS model werecross-validated at 30 m spatial resolution scale.

Selection of phenological metrics based on a phenological analysis
We used time-series Google Earth imagery to monitor the dynamic processes of land-cover types and gradual cropland retirement in the Yan'an area. For the NDVI curves of crops in this study area, as shown in Figure 3a, there was a significant discrepancy with other vegetation types during the key phenological stages (from DOY161 to DOY289) (Yang et al. 2020). In the NDVI time-series curves, the crop NDVI began to increase on DOY 161 and peaked between DOY 225 and 241, after which crop NDVI values began to decrease, reaching their lowest levels around DOY 273. In our analysis, the NDVI curve for cropland clearly differed from the curves for forests and grasslands. As shown in Figure 3b, we found that the NDVI time-series curve of the cropland retirement sample was relatively steep from DOY161 to DOY289 when more than 95% crops are covered. Relative to the crop coverage rate lower than 10% (Figure 3d), the NDVI values of this cropland retirement sample at the start and end of the crop growth period were significantly higher, and the NDVI time-series curve was smooth (Qin et al. 2016). The start and end of the crop growth period are the most important phenological parameters that reflect cropland retirement. Therefore, we determined the NDVI values of crops at the start (DOY161), middle (DOY225) and end (DOY273) as phenological metrics, and the length of the periods described above, as shown in Table 1. We used Landsat data from the same or nearby dates as a phenological metric for the analysis of cropland retirement.

Cropland extraction index based on phenological characteristics
Based on the growth characteristics of vegetation during the growing season, multi-time or time series satellite images have been widely used to extract land cover information (Huang et al. 2017). During the growing period of a crop, differences in phenology between crops and other surface vegetation categories, such as grass or forest, are huge as a result of the differing sowing and harvesting times among crops (as shown in Figure 3a). We constructed the CEI based on phenological metrics. The parameters were configured as follows: (1) where NDVI grow max , NDVI grow begin , and NDVI matu end represent the NDVI at the middle, start, and end of the growth period, respectively, with Δt grow maxÀ grow begin and Δt grow maxÀ matu end representing the corresponding time intervals. The formula NDVI grow max À NDVI grow begin À � =Δt grow maxÀ grow begin is defined based on the difference in growth fluctuations between the target crop and background data, and we selected two periods with a large difference for comparison. To facilitate the calculation and threshold extraction, we multiplied the result by 10,000.
The setting of a reasonable CEI threshold is a key step (Maxwell and Sylvester 2012) in distinguishing between cropland and other features. Extracting pure crop pixels is a relatively simple process. But choosing a reasonable CEI threshold for extracting crop areas from mixed pixels is difficult. To simplify the interference of mixed pixels, a Landsat pixel is defined as "cropland" if more than 50% of the pixel is covered by crops (Yang et al. 2019). The percentage of crops within the pixel was verified by visual interpretation using high-resolution Google Earth images.

Recognition of the extent of cropland retirement based on the PCRRS model
Cropland surface phenology can change dramatically due to changes in crop cover proportions or the alteration of land use practices. We constructed a CEI model based on the NDVI values of crop phenology characteristic parameters. Consequently, we explored using the CEI value as a variable in the PCRRS model to map cropland retirement.
The six main steps of our approach are described below. First, we calculated the CEI index of crops (such as corn) in the start year M and end year N (M and N refer to the two years before and after the cropland retirement phenomenon occurred). Then, the CEI thresholds were used to separate cropland from other land-use classes. Second, we used the inter-annual CEI value as an input parameter for the PCRRS model to calculate the PCRRS value Eq. (2). Third, a relationship between PCRRS values and the ratio of cropland retirement area in Landsat pixels was derived by visually interpreting a part of analytical samples (126 samples in total) from Google Earth high-resolution imagery. Then, within the cropland area present in year M, we were able to determine the extent of cropland retirement from M to N based on the appropriate PCRRS threshold. Fourth, it is assumed that other crops (such as wheat) should not be grown in the same area after corn is retirement. The interannual variation in phenological parameters caused by alternating cultivation between crop types can cause the misidentification of the model in this paper. Growing other type of crops in the retired area is defined as a shift in crop type. Repeat steps first to fourth to obtain active and abandoned cropland for other crops. Fifth, we distinguished cropland retirement from stable cropland based on the area of cropland in N(The phenomenon of cropland retirement occurs within the scope of farmland and ends in the scope of non-farmland), and eliminated interference from water bodies and built-up areas using land-use data for year N. Finally, in the sixth step, farmland fallow is defined as the early result of cropland abandonment which was marked as active farmland in the next stage, so we can eliminate the impact of fallow land.
where i represents various crops. In steps 1-3, we characterized the changes in phenological parameters caused by changes in crop coverage, extracted the extent and timing of cropland retirement based on an appropriate threshold, as shown in Figure 4. In step 4-5, because alterations in the crop type and land-use category may result in false judgments, we highlighted these changes and provided a solution. In the final step, we removed the interference from fallow cropland. Thus, the extent and timing of cropland retirement could be obtained with high precision.

Land-cover classification based on support vector machines classifier
SVM is a computer learning method for supervised analysis and classification of high-dimensional datasets. Equivalent to other learning methods such as neural networks and decision trees, SVM employ nonparametric optimization algorithms to locate optimal boundaries between classes. (Prishchepov et al. 2013). We used SVM and three NDVI images, such as at the start, middle and end of the crop growth period, to obtain annual coverage maps of crops (maize and wheat) in the study area (Alcantara et al. 2012). Then, we identified active farmland and abandoned farmland based on annual cover transfer results for crops. Training data were extracted from Global Land Cover (GLC) 2017, and the land cover categories include forest, grassland, water, impervious surface and crops. Among them, the crop types are divided into corn and wheat through Google historical images combined with crop phenological parameters.

Validation of the cropland retirement map
Confusion matrices were used to examine the accuracy of the retired cropland maps by calculated the producer's accuracy (PA), user's accuracy (UA), and overall accuracy (OA). Taking into account the sampling bias caused by the sparsity of samples (Card 1982), we adopted disproportionate stratified sampling for verification (Olofsson et al. 2014). We stochastically selected 100 pixels each from the stable farmland (500 samples in total) and stable nonfarmland classes, and 150 pixels per year from the retired cropland class (600 samples in total). The sample size is 50 m*50 m, covering one Landsat pixel to reduce the recognition error caused by registration. The sample distribution is shown in Figure 4. These validation points, which record crop percentage and cropland abandonment percentage, was derived from high-resolution Google Earth imagery by visual interpretation.

Parameter extraction of the PCRRS model
We initially evaluated the impact of the Landsat image date selected on the range of CEI values, and the results showed good consistency. All available Landsat images in this paper were showed in Figure 5a. Deviations in the imaging dates may cause fluctuations in the CEI values obtained for vegetation types. The CEI values of forest and cropland were analyzed (Figure 5b). Cropland had higher CEI values than forest, with values exceeding 140. Pure forest had low CEI values, ranging from 20 to 40. The mean CEI value of crops was highest in 2005 and lowest in 2018. For the same vegetation type, the distribution range of CEI values may differ due to differences in the imaging date. For example, in 2018, an image collected on DOY302 was used to construct the CEI. Because the NDVI value of forests also dropped to its lowest level at that time, forests had larger CEI values than during the other four periods, while the corresponding CEI values for cropland were smaller.
Based on the distribution of CEI values described above, relationships between CEI values and the proportion of cropland within mixed pixels were obtained for both 2005 and 2018 (Figure 6a). The results showed that CEI values increased with increasing the proportion of cropland. Moreover, during different periods, the same CEI value corresponded to  similar crop proportions. Based on the analysis above, when the proportion of crop area exceeded 50%, we used the same CEI threshold of >70 for all periods.
The PCRRS method by applying the difference among CEI values in different years to each percentage of cropland retirement, obtaining thresholds that relate directly to the percentage of cropland retirement within a pixel. A high-resolution Google Earth image taken on 07/01/2019 of an area of cropland retirement is shown in Figure 6b. For this image, the percentage of cropland retirement was obtained by applying the PCRRS to CEIs for 2001 and 2018. The results showed that the PCRRS method enabled a calculation of the approximate percentage of cropland retirement within a Landsat pixel.
The relationship between a set of PCRRS value and the proportion of cropland retirement is shown in Figure 6c. Land cover changes as small as 25% of a pixel can be effectively monitored (Vargas, Montalban, and Leon 2019). However, estimating the percentage of cropland retirement within a pixel in areas where cropland is highly dispersed is often difficult (Ozdogan and Woodcock 2006). Therefore, we marked the pixels with the proportion of cropland retirement more than 25% as the pixel of retired cropland. Considering the consistency of the constructed datasets and the same calculation method, we used the same PCRRS threshold across periods, with a corresponding PCRRS value of > 20.

Cropland retirement mapping and accuracy assessment
The OA of active cropland maps throughout the study ranged from 86.8% (2005)  were obtained using the PCRRS model. The retired cropland class also had a high accuracy, with the PA ranging from 84. 2% (2005-2009) to 95.2% (2001-2005) and the UA ranging from 83. 9% (2014-2018) to 92% (2001-2005) (Table 2). In Yan'an, land-use policies have played a decisive role in cropland retirement (Lu et al. 2012). The GFG program launched by the Chinese government aimed to restore ecosystems by converting cropland in mountain and barren areas to forests or grasslands.
All cropland retirement samples were used to cross-validate the accuracy of the model. Croplands in mountainous areas are especially likely to be retired because they tend to be unprofitable (Zimmermann 2007), The overall accuracy of the model in steep terrain is as high as 93.4%. The recognition accuracy of the PCRRS model in flat terrain indicated that the impact of crop type changes on OA cannot be ignored. In addition, based on the land use type results obtained by SVM with multi-date images, we obtained the recognition results of the conversion of cropland to woodland and grassland. And the accuracy of the recognition results was still lower than with our method (Table 3), although the result of cropland retirement also had good accuracy.
We identified widespread cropland retirement in the Yan'an region, with the timing and extent of cropland retirement varying greatly from 2001 to 2018. The results showed that of the ~722,800 ha of cropland in 2001, roughly 82.86% had been retired by 2018. We found high rates of cropland retirement after the beginning of the GFG program and a substantial contraction of cropland (Figure 7

Temporal and spatial patterns of cropland retirement
We observed strong differences in cropland retirement among counties in our study area. In particular, the temporal patterns of retirement differed among counties. Six counties had substantial retirement rates before 2005: Zichang, Zhidan, Wuqi, Ansai, Yanchuan, and Yanan in the north, among which Yanchuan had the highest retirement rates (up to 66.8%). On the other hand, southern counties had relatively low retirement rates before 2005 and higher rates after 2005. In general, our analyses revealed that the implementation time of land-use policy differed for each county. All counties in Yan'an showed a slowing trend of cropland retirement, although the extent of the change varied regionally (Figure 8a). Cropland retirement was also related to terrain (Figure 8b). Before 2005, only 21% and 34% of cropland was converted to forest or grassland on slopes of 0-6° and 6-15°, respectively, while nearly half of cropland on slopes of 15-25° was converted to forests or grassland. The greatest reforestation occurred on slopes > 25°. We found that, compared with flat terrain, cropland retirement occurred earlier and increased more on steeper slopes, reflecting the success of initiatives to restore the ecological environment of the mountainous area. This outcome was in line with the design of the GFG program, which encourages the conversion of cropland to forest on steep slopes. Inter-cropping of shrubs and crops is also common in retired cropland areas with flatter terrain, the phenomenon of gradually cropland retirement is widespread in this study area (Figure 9). We found that 87% and 88% of cropland was converted to forest or grassland on slopes of 0-6° and 6-15°, respectively, from 2001 to 2018. In all topography classes there was a trend for the pace of cropland retirement to slow.  Retired cropland may lead to either a abrupt or gradual cessation of cultivation under the combined influence of policies, topography and environmental (Estel et al. 2015). Figure 9 shows the gradual retirement process of cropland in this study area. The land cover type within an approximate MODIS pixel was dominated by maize crops from 2008 to 2014 (Figure 9a

Discussion
To accurately estimate the area of cropland retirement, a CEI threshold of crop was used to exclude the interference from vegetation such as woodland and shrubs. Then, how to detect the extent of cropland retirement in the mixed pixels is particularly important. We developed a PCRRS model for identifying and extracting information regarding cropland retirement. Our model eliminated the influence of mixed-pixels problem and allows mapping of retired cropland with high precision for a complex underlying surface.

Importance of phenological interannual variation for the recognition of retired cropland
Based on the specific phenological cycle of vegetation types, only three images per year are needed to classify land uses, with classification accuracies of up to 80% achieved (Prishchepov et al. 2012b). Multiple important phenological periods have been widely used in land use classification (Foerster et al. 2012). Considering that the mixed pixels of cropland also have NDVI fluctuations during the growth cycle, we simulated the proportion of cropland according to the rate of change in NDVI values in the growing and harvesting periods. Based on change reasonable to the value of threshold, the interference caused by other vegetation types could be effectively eliminated. When monitoring changes of agricultural landuses based on inter-annual variation of key phenological "nodes" in farmland, a change of crop type will generate a large change in NDVI values . Therefore, we conducted a cooperative analysis of multiple crop types to eliminate the identification error caused by conversion between crop types. Previous studies have shown that the phenological characteristics change with climate (Chmielewski and Rtzer 2001;Khanduri, Sharma, and Singh 2008). Changes in phenological parameters (start, mid and end of the growing season) due to climate can cause misidentification by our model. Considering that the early stage of cropland retirement was characterized by low vegetation coverage, which prevented the phenology of these areas from being significantly advanced or delayed by climate change (Penuelas and Rutishauser 2009). Therefore, a model based on the characteristics of phenological interannual variation can effectively extract the distribution of cropland retirement.
Trees are common in cropland areas (Xu et al. 2018), and the majority of Landsat pixels are a mixture of various land-cover types. Therefore, a mixed pixel analysis is preferred for analyzing changes in regions where cropland is highly fragmented (Jain et al. 2013). The difference in spatial resolution between MCD12Q1 and GLC 30 m may have been the main factor for the large difference between the retired cropland results (Yin et al. 2020) (Figure 10), while the temporal characteristics of gradual cropland retirement should be ignored due to the difficulty of identifying agricultural land retirement via machine learning. This is especially true for the analysis of gradual changes over time, as defined classification rules tend to oversimplify the analysis of the process (Prishchepov et al. 2012b), leading to a lag in the temporal results for retired cropland. The rapid cropland loss before 2005 was in line with the initial regulations of an ecological restoration project (Liu et al. 2008). However, following a period of forest expansion at the cost of diminishing grain production, the target for increased forest cover was reduced after 2005 (Xu et al. 2006).

Uniqueness analysis of the PCRRS model
The cropland retirement in China is unique. In contrast to the natural succession of abandoned agricultural land in other countries (Blair, Shackleton, and Mograbi 2018), Chinese managers have developed afforestation plans on abandoned cropland (Zhou, Rompaey, and Wang 2009). Similar human interventions have resulted in retired cropland being covered with woody vegetation in a relatively short period of time. Compared with the agricultural abandonment caused by social (Prishchepov et al. 2012a) and armed conflicts 2019, the agricultural land retirement in China is based on the purpose of improving the ecological environment (Gao et al. 2018). On the premise of ensuring food security (Chen et al. 2015), cropland with low efficiency is planned and purposefully abandoned (Bakker et al. 2011). Obtaining a valid map of cropland retirement from a social, economic, and environmental perspective helps to assess the social, economic, and environmental outcomes of agricultural land abandonment. Provides a basis for managers to develop sustainable agricultural abandonment plans.
MODIS with high temporal resolution provides a reliable basis for us to obtain phenological parameters of a single crop type. Based on the dates of MODIS phenological parameters (start, middle and end of crop growth period), we collected Landsat data for the corresponding dates. Considering that crops are at the beginning and harvesting stage of the growing season, low vegetation cover corresponds to low NDVI values. Compared with the 16day synthetic MODIS data, the Landsat image data more realistically reflects the surface vegetation coverage of the corresponding date; in the middle of the growing season, high vegetation coverage corresponds to a high NDIV value. Considering the saturation property of NDVI, the Landsat NDVI image in the middle of crop growth period can approximately reflect the lush growth state of crops. Therefore, it seems reasonable to choose Landsat data from the same date as the MODIS phenological parameters for the study.
More importantly, although the annual maximum NDVI increased in gradually retired cropland areas, the trend of the annual maximum NDVI increase may be less significant in the later stages of retirement as a result of the saturation of NDVI signal (Beck et al. 2006). On the contrary, the vegetation growing season in regions of cropland retirement was significantly longer owing to NDVI signal is most sensitive to moderate vegetation coverage (Figure 9) (Debeurs and Henebry 2004). Several studies have found trends of an earlier greening onset and longer growing season through satellite and station observations of phenology (Karlsen et al. 2007). Longer growing periods of surface vegetation will increase drought stress and evapotranspiration with global warming (Barber, Juday, and Finney 2000;Zhang et al. 2009), and the extent of carbon sequestration also increased (White, Running, and Thornton 1999). Therefore, observations of the distribution of cropland retirement might be more effective than maximum NDVI values (or annual mean, growing season integral, etc.) as an indicator of stress or soil restoration on the Loess Plateau (Jong et al. 2011).
The environmental and climatic factors in relatively arid areas are such that the cultivated land is mainly planted with single-season crops. At the same time, the short rainy season makes it possible to capture images of key phenological periods. In addition, poor quality pixels obtained in the key phenological period are also highly substitutable because the NDVI values of single-season crops are relatively stable before sowing and after harvesting. Therefore, in arid and semi-arid regions, our model is able to achieve a high recognition accuracy. However, in wetter regions, our model may not be applicable due to the complexity of farming practices, perennial crops and influence of cloud and rain.

Uncertainty analysis of the PCRRS model
While its relatively high mapping accuracy highlights the potential of our new approach, there are a few uncertainties. First, although high accuracy was achieved by the PCRRS model in the study area, the model may not be able to accurately simulate areas with complex planting structures in humid regions. Fine transitions between crop types in complex Figure 10. Estimates of the cropland retirement area for the study region based on different models. farming areas will be easily overlooked as a result of the comparatively coarse temporal resolution of Landsat imagery, but with the enrichment of data sources, such as Sentinel series, these interference factors will be greatly reduced. Second, cropland pixels with phenological changes will be marked as cropland conversion, which will impact on the mapping accuracy of cropland retirement. For example, the urban sprawl surrounding the cities in Yan'an was sometimes incorrectly labeled as vegetation restoration caused by an ecological restoration program, despite the use of GLC30 data for masking urban areas. Third, the assumption that phenology and reflectance remain similar over many years may have led to incorrect classifications, so we cannot completely rule out overestimation of the accuracy of the results. At the same time, the use of the same phenological parameters in different regions will cause the model to be disturbed by regional heterogeneity, which will lead to incorrect assessments of the results. By averaging regional phenological parameters, the interference of regional heterogeneity is alleviated to a certain extent. Overall, our analysis demonstrated that cropland retirement can be mapped with a relatively high accuracy at the scale of Landsat imagery using the PCRRS model.

Conclusion
The PCRRS model was developed for determining the spatiotemporal pattern of retired cropland in this paper. We eliminated interference from other land-cover types by combining phenological metrics for crops with multi-temporal NDVI images. Then, the phenological change characteristics were used to identify the extent of cropland retirement. The results suggested that our algorithm has a high accuracy of over 85%, and our model was most suitable for sloping cropland areas. Gradual cropland retirement was a common phenomenon in the study area, and our approach successfully mapped its extent and timing. This new approach for satellite-based imagery highlighted the possibility of obtaining more accurate maps of cropland retirement change in areas with small-scale farms and high landscape heterogeneity. Future studies will focus on combining the PCRRS model with MODIS time-series images, applying them at the regional scale to reveal the relationships among interference by human activity, policy implementation, regional environmental differences, and agricultural land-use change.