Retrieval of eucalyptus planting history and stand age using random localization segmentation and continuous land-cover classification based on Landsat time-series data

ABSTRACT Obtaining robust change-detection results and reconstructing planting history are important bases for conducting forest resource monitoring and management. The existence of multiple change points in a very short period can lead to a global segmentation method incorrectly locate the change points, because they could impact each other during model initialization. This is especially true for monitoring plantations such as eucalyptus, which has a unique growth cycle with short rotation periods and frequent disturbances. In this study, we proposed a method to find critical change points in a normalized difference vegetation index (NDVI) time series by combining random localization segmentation and the Chow test. Features of the NDVI time series calculated on the divided segments and change points were used to train a Random Forest classifier for continuous land-cover classification. The proposed method was successfully applied to a eucalyptus plantation for identifying the management history, including harvest time, generation, rotation cycle, and stand age. The results show that our method is robust for different lengths of NDVI time series, and can detect short-interval (cut and stability) change points more accurately than the global segmentation method. The overall accuracy of identification was 80.5%, and successive generations in 2021 were mainly first- and second-generation, accounting for 69.0% and 27.9% of the total eucalyptus area, respectively. The rotation cycle of eucalyptus plantation is usually 5–8 years for 66.9% of the total area. The eucalyptus age was accurately estimated with an R2 value of 0.91 and RMSE of 13.3 months. One-year-old eucalyptus plantations accounted for the highest percentage of 14.5%, followed by seven-year-old plantations (12.9%). This study provides an important research basis for accurately monitoring the rotation processes of short-period plantations, assessing their timber yield and conducting carbon- and water-cycle research.


Introduction
Eucalyptus is an important plantation for national timber security, national economic development, and ecological protection (Liu and Xie 2020;Wang et al. 2020). Short-period and pure forest succession are adopted worldwide in eucalyptus management. Accurate reconstruction of the planting and growth history is vital to better evaluate the quality of the eucalyptus ecosystem and achieve multi-objective sustainable management. Eucalyptus plantation has a very unique disturbance-growth cycle: the trees are harvested every several (4-7) years and replanted for the next rotation a few months later, and the canopy is closed with saturated reflectance signal in 1-2 years (Le Maire et al. 2014;Qiao et al. 2016). The growth trajectories are irregular and various depending on the soil fertility and management measures. There are several "change points" in very short periods produced by clear-cut, planting, and canopy closure. It is important to correctly locate these points for reconstructing the forest planting history and management.
Remote sensing, with over 40 years of data available, has become the most popular data for conducting retrospective analysis and characterization of long-term changes in forest resources (Gomez, White, and Wulder 2016;Wulder et al. 2019). Various algorithms in satellite-based, time-series analysis have been developed for forest change detection and characterization (Verbesselt et al. 2010;Zhu and Woodcock 2014;Zhu, Woodcock, and Olofsson 2012). These algorithms can be categorized into offline and online methods according to their monitoring strategies (Bullock, Woodcock, and Holden 2020;Zhu 2017). Offline algorithms consider the entire timeseries dataset at once (thus called global segmentation), and look back in time to recognize where changes occurred. They include Breaks For Additive Season and Trend (BFAST) (Verbesselt et al. 2010), BFAST Lite (Masiliunas et al. 2021), Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) , Cao et al. (2020), and Deng et al. (2020). BFAST and LandTrendr are the two most commonly used offline change-detection methods in forest change detection. Both receive the whole time series at once for model initialization and do not need to specify a stable period. BFAST estimates break magnitude using the models fitted before and after the preset breakpoints (Verbesselt et al. 2010). It calculates the corresponding minimum residual sum of squares (RSS) for different numbers of breakpoints and determines the optimal segmentation scheme of the model based on minimization Bayesian Information Criterion (BIC). LandTrendr adopts a residual-error criterion to identify vertices based on regression vertex identification . The point with the largest absolute deviation from the fitted regression line is designated as the breakpoint and divides the time series into two segments. Then, new regressions are fitted for the new segments and identify the next potential breakpoint until it reaches the maximal segment count and control parameters. This kind of change-detection process is also called Binary Segmentation (BS). A recent study showed that BFAST performed relatively consistently across a range of different change-severity levels for the break/trend set than other change-detection methods (Saxena et al. 2018). However, this test is based on only one breakpoint in the time series, and it is unclear whether good results for multiple breakpoints detection can be achieved.
Global segmentation is a commonly used method for time-series multiple change-point detections. It searches for the first change point from the entire dataset. After the first breakpoint is detected, the time-series observations are then divided into two subsegments (hence named "binary") based on the detected breakpoint. Then, a similar process is executed on each subsegment. This BS strategy was widely used in time-series remote sensing data for change detection in different ecosystems (Cao et al. 2020;Deng et al. 2020;Kennedy, Yang, and Cohen 2010). However, the global segmentation method suffers some inherent issues because it uses all the data for model fits (Fryzlewicz 2014). If the timeseries data have multiple change points, the model initialization will possibly fit the wrong model. This may lead to catastrophic impacts as the changedetection method can incorrectly locate all change points, because the change points could impact each other during the model initialization. The global segmentation method can be consistent when the minimum spacing between any two adjacent change points is large enough; that is, if the true breakpoints occur close to each other, which is the case in eucalyptus plantations, the result can be spurious change-point detection (Fryzlewicz 2014). In addition, the inability to effectively detect breaks at the beginning and end of the time series is a commonly known limitation of this type of method (Li et al. 2019b;Robbins et al. 2011).
The amount of data taken into the model initialization is a trade-off between observations for model fit and change detection. Taking too much data into the model initialization increases the probability of fitting the model with incorrect change points, and brings large omission errors because fewer observations are available for change detection. Taking too little data into the training period might fit an insensible model and detect spurious change points when compared to the predictions of the monitoring process. The choices of observations used to model initialization can be treated as determination of an optimal window size. Incorporation of variable window sizes and periods for change detection may provide a good solution to the trade-offs between small and large (whole) windows. Fryzlewicz (2014) proposed a random segmentation strategy that solves some of these problems. The approach -Wild Binary Segmentation (WBS) -is based on maximizing a test statistic over randomly sampled subsequences of time-series observations. By considering subsequences, WBS is able to detect change points that only affect local parts of the sequence and avoids the mutual impacts from multiple change points (Fan et al. 2022). However, this strategy to detect frequent forest change points in time-series remote sensing data has not been thoroughly tested. The impacts of seasonal change, noises, and variant data lengths in these data also bring many obstacles to its application.
Wide concerns have been raised about the negative environmental impacts -consumption of soil nutrients, biodiversity reduction, depletion of soil water, and production of allelopathic chemicals by eucalyptus plantings (Diego et al. 2020;Li et al. 2019a;Nouvellon et al. 2012). Some studies have tried to map the distribution of eucalyptus and stand age at regional scale (Deng et al. 2020;Le Maire et al. 2014). However, it is not known what exactly the global or regional rotation of eucalyptus is. This study proposed a new method to detect frequent (adjacent) change points of eucalyptus by combining a random localization segmentation approach and the Chow test based on Landsat timeseries data, and comparing results with the changedetection results from the global segmentation method. Specifically, we want to answer the following three questions: (1) Can the proposed method accurately detect the frequent change points of eucalyptus during the period of 1986-2021? (2) Is the random localization segmentation more robust than the global segmentation method in a short-rotation eucalyptus plantation? (3) Can the proposed framework successfully provide historical information on eucalyptus planting, such as harvest time, generation, rotation cycle, and stand age?

Study area
Yunxiao County, located in southeastern Fujian Province, China, was taken as our study area ( Figure 1). The county has a subtropical marine monsoon climate with annual average temperature of 21.3°C and annual precipitation of 1730.6 mm, and is frost free more than 340 days of the year. It has 583 km 2 of forestland, accounting for 40.1% of the county's area. According to our survey and China's statistical yearbook, Yunxiao County focused on the development of the fruit industry, and a large amount of farmland was converted to orchards in the 1980s. With the decrease in economic benefits of the fruit industry and the gradual development of the eucalyptus industry, many forests and orchards were transformed to eucalyptus plantations in the early 2000s. In the process of eucalyptus management, the first round adopts seedling cultivation and the second round generally adopts coppice regeneration. Seedlings are usually planted from March to June to ensure their survival rate. The height of seedlings is about 20-30 cm, and they take 1-2 months to establish root systems. The trees can grow to heights of 6-7 m in the first year. The coppice regeneration method can be divided into two types -using slashand-burn or not. The eucalyptus plants sprout profusely from stumps right after the trees are cut, and can grow to heights of more than 2 m in 3-4 months. The harvest, planting (regeneration), and stability in such a short period leads to several change points in the growth cycle.

Datasets and preprocessing
Landsat NDVI data with cloud cover of less than 80% from August 1986 to January 2021 in Yunxiao County were selected from USGS (593 total images, including TM 259, ETM + SLC-off 203, ETM + SLC-on 32, and OLI 99) (Figure 2). There were 17 images on average per year with the most images (29) from 2018, and the fewest images (3) from 1986. There were 49 images on average in each month, and the most (68) and the fewest (31) were in October and March, respectively. Clouds, shadows, and strips in each image were removed based on the quality band.
We randomly collected 4,172 land-cover samples through field survey records and visual interpretations on Google Earth in July 2021. Land cover was classified into six types: eucalyptus (851), non-eucalyptus forest (601), farmland (633), orchard (670), non-vegetation (1288), and logging slash (129). Seventy percent (70%) of the samples were used to train a random forest classifier model, and the remaining 30% samples were used to test the model. We randomly selected 154 eucalyptus samples to collect additional information, including generations, times of planting and cutting during 1986-2021, rotation cycles, stand age (months) in 2021, planting type (coppice regeneration or seedling cultivation), management practices, etc. Further manual verification was carried out in combination with time-series Landsat data to ensure there were no obvious mistakes in the field survey questionnaire responses. Based on data from the field survey and Landsat NDVI time series, we divided the change processes of eucalyptus into the following types: (1) an obvious cutting-rapid-growth-stable process; (2) no cutting and rapid-growth process with the transformation of non-eucalyptus (or eucalyptus) to eucalyptus; and (3) transformation from eucalyptus (or non-eucalyptus) to eucalyptus with negative abrupt change (cutting). See section 3.1.3 for more detailed information on these processes. Figure 3 illustrates the general framework we created, including the following: (1) Using NDVI time series as the input data after invalid data such as cloud, shadow, and outliers were eliminated. (2) The NDVI time series of each pixel was randomly divided into M segments; the maximum cumulative sum (CUSUM) of each segment and its corresponding breakpoint was calculated and determined by strengthened Schwarz information criterion (sSIC). (3) The harmonic model and Chow test were used to judge if the determined breakpoints were spurious change points based on F statistics. (4) The number, location, and order of the detected breakpoints were recorded. (5) Based on the features of change points and each segment, the spatialtemporal distribution of eucalyptus during 1986-2021 was classified using a random forest classifier. (6) Combined with the breakpoint characteristics and the eucalyptus growth process, the times of cutting and recovery to stability, generations, rotation cycles, and forest age in 2021 were determined based on the information gleaned from processes (1)-(5).  number of change points (N) and locate places or times of change points that divide the time-series data into several subsequences. Each subsequence can be referred to as a state of the time series, or period of time when the parameters governing the process do not change. Two consecutively dissimilar states are divided by a breakpoint. In this study, the breakpoints were defined as the abrupt changes of NDVI time series driven by land-cover conversion, forest cutting, and significant change in growth state (from rapid growth period to no obvious trend period). The strategy of a random localization mechanism used to detect change can be summarized as follows (Fryzlewicz 2014;Korkas and Fryzlewicz 2017): First, the CUSUM statistics are calculated over M subsequences [s, e]. The starting s and ending e points are randomly drawn from a uniform distribution U (0, T − 1). The assumption is that some of the subsequences will only contain single change points for a very large M. These subsequences are especially appropriate for CUSUM-based change detection. The location that has achieved the largest maximum CUSUM over all subsequences is labeled as the first breakpoint candidate, and is tested for significance. After detecting a breakpoint, the changedetection process is resolved into two subproblems that find further change points for each segment in the same way. The method proceeds until a certain criterion is satisfied. A program for reducing the computational complexity is proposed by calculating the CUSUM statistics once at the beginning of the method for the M random subsequences. As the random localization segmentation algorithm proceeds at a generic subsequence [s, e], the calculated CUSUM statistics can be reused by ensuring the starting and end points fall within [s, e].

Random localization segmentation-based
We randomly drew a number of subsequences (M), that is, vectors (NDVI s , NDVI s+1 , . . ., NDVI e ), where s and e were the start and end observations for a subsequence, and computed the CUSUM statistic for each subsequence (Eq.1). Then, we maximized each CUSUM, chose the largest among the M CUSUMs (Eq.2), and took it as the first breakpoint candidate to be tested against information criteria or a predefined threshold. Once the breakpoint was determined to be significant, the same procedure was then performed to the left and right time series of it. The detection of a change point directly from CUSUM statistics was based on the following considerations: (1) Evergreen trees covered most of the study area. Although there were some seasonal changes, their impact was relatively small compared to disturbance. (2) It was necessary to reduce the limiting conditions of change point judgment to capture all potential breakpoints and reduce the problem of high omission in breakpoint monitoring.
(3) The detected breakpoints were recorded in descending order -each breakpoint is less important than the breakpoint before it. In addition, a method based on the Chow test was proposed to eliminate spurious change points.
g NDVI b s m ;e m ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ð Þ is the subsegment with the maximum g NDVI b s m ;e m and its corresponding change point. By "localizing" the CUSUM statistic in a randomized manner, it can solve the mutual impacts of multiple breakpoints which are confronted in the "global" CUSUM algorithms. It can also dramatically alleviate the allowed minimum intervals between neighboring breakpoints compared to current methods, and the permitted change magnitudes. By creating subsequences with different lengths, we do not need to fix the span or window size, which is present in some prevailing methods to localize the CUSUM statistic. Furthermore, the algorithm returns estimated breakpoints in a descending order by the maxima

Determination of candidate change points.
The sSIC was used to determine the locations and number of breakpoints. The attraction of the sSIC approach lies in the fact that the choice of the parameter is easier than the thresholding approach. It is difficult to determine whether change occurs or not by selecting the optimal threshold value. These values may be application dependent and they may change over time. The only parameter of the sSIC is the constant α, which determines the degree of penalty.
where α = 1 represents the standard Schwarz Information Criterion (SIC) penalty. Here we chose α = 1.01 to ensure the results remained close to those obtained by SIC; K is the number of candidate breakpoints; σ k 2 is the corresponding maximum likelihood estimator of the residual variance; and N is the number of observations in the time-series data.

Removal of spurious change points.
A relatively loose condition based on CUSUM was used to identify the change points, which might result in many more breakpoints than expected due to its sensitivity to local changes and noises in the NDVI time series. Although the sSIC was performed to test the significance of detected breakpoints, it was necessary to further eliminate the spurious change points. The Chow test is a statistical test developed by economist Gregory Chow that is used to test whether the coefficients in two different regression models on different datasets are equal (Chow 1960). The Chow test compared the effectiveness of two separately adjacent models with one single model that spans the entire time period (two segments). We used two harmonic terms (i.e. j = 2) to robustly detect phenological changes within the NDVI time series to fit the observations before and after the potential change points (Eq. 4), and computed the error sum of squares (ESS). Another ordinary least squares (OLS) model for all observations of two adjacent segments (Eq. 5) with a restricted sum of squares (RSS) was calculated.
where i is the ith segment divided by the detected breakpoint, and τ � k is the kth breakpoint; b ρ i; t ð Þ NDVI is the fitted NDVI value by the harmonic model with four coefficients including the intercept β, harmonic components a and b, and slope c at time t; T equals 365 days per year; b ρ t ð Þ NDVI is the fitted NDVI for two adjacent subsequences; F is the Chow test statistic; RSS is the sum of squared residuals (SSR) of the two segments before and after breakpoint; ESS = SSR of the model before breakpoint + SSR of the model after breakpoint; k is the number of model coefficients; and N s is the number of observations in the two subsequences before and after breakpoint. If the statistic F determines that the coefficients are not equal between the regression lines, there is significant evidence that a structural break exists in the data. The Chow test determines whether the equality of coefficients is significant or not. If the p-value associated with the Chow test statistic is less than a certain significance level (0.05 in this study), we reject the null hypothesis and conclude that there is a structural break point in the data. Otherwise, the two segments should be merged and the change point removed from the candidate breakpoints.

Classification of segments using random forest
We classified the land-cover types (i.e. eucalyptus, non-eucalyptus forest, farmland, orchard, nonvegetation, and logging slash) of each temporal segment based on the features of subsequence NDVI and terrain. The features include harmonic terms, mean NDVI of segment, change trend, length of segment, magnitude of change, elevation, slope, and aspect. Combined with the field survey data of land-cover types gathered in 2021, the random forest classification model of the latest temporal segment was established based on 70% of the samples (Breiman 2001). The parameters of mtry and ntree were set as default values with square root of the number of variables and 500, respectively. The established model was used to classify all segments (Li et al. 2019;Zhu and Woodcock 2014). The land-cover type during the period of each segment was not changed, so the continuous land cover can be obtained as each segment can be classified by the trained random forest model.

Identification of eucalyptus planting history and stand age
Once all the critical change points of the NDVI time series were accurately located, the forest growth cycle could be depicted clearly. The attribution of each eucalyptus-related change point was determined by the land-use and -cover types before and after the change point, and its change characteristics. Cutting and stability are two main change types of the NDVI time series for eucalyptus. The stability point is the turning from the rapid-growth phase to a stable phase without disturbance. The rules to recognize them were established as follows: (1) If the land cover was converted from logging slash to eucalyptus, the start time and end time of the segment were labeled as cut time and stable time; the two segments were grouped as one rotation cycle. (2) The change point from non-eucalyptus to eucalyptus was labeled as cut time/stable time (the cut time was the same as stable time); the second segment was labeled as the first cultivation period. (3) The change point was labeled as both cut time and stable time if the land cover converted from eucalyptus to eucalyptus, and the generation was increased by one. (4) The forest age increased by six months if the NDVI time series only had a stable point and the previous planting was non-eucalyptus. (5) We also postprocessed the classification result: when the segments were classified as logging slash and the length of segment lasted more than three years, they were relabeled as eucalyptus. The segment was labeled as logging slash if the duration was less than three years and it was not classified due to lack of feature data.

Accuracy assessment
The classification accuracy was evaluated by overall accuracy (OA), producer's accuracy (PA), and user's accuracy (UA) based on 30% the land-cover samples. The accuracies of generations, rotation cycle and forest age were assessed based on the surveyed 154 eucalyptus samples. OA, PA, and UA were calculated for generations, and the R 2 and RMSE were adopted for rotation cycle and forest age. Only eucalyptus samples that had complete cultivation periods were counted for rotation cycle evaluation. To validate the change-detection accuracy, we further randomly selected half of the 154 eucalyptus samples and manually interpreted the cut and stability points through NDVI trajectories and TimeSync , and calculated the PA, UA, and OA to assess the accuracy of the abrupt change results.

Comparison with global segmentation change detection
The global segmentation (e.g. BFAST) algorithm has excellent performance in many comparison experiments with change-detection algorithms (Cohen et al. 2017;Saxena et al. 2018). This article compares the differences in results between random localization segmentation and global segmentation in calculating the same NDVI time series. We further assessed their robustness at different lengths of NDVI time series in change detection. Here we briefly introduce the basic theory of global segmentation in BFAST (Verbesselt et al. 2010), which was adopted from "strucchange" and proposed by Bai and Perron (2003). The basic idea is that the total number of possible segments is at most S(S + 1)/2 with a sample size S. The global SSR for any m-partition (S 1 , . . ., S m ) and for any value of m must be a particular linear combination of these S(S + 1)/2 SSRs. The estimates of the break dates, the m-partition (S 1 , . . ., S m ), correspond to this linear combination with a minimal value. In the "strucchange" procedure, a dynamic programming algorithm was adopted to compare possible combinations corresponding to different m-partitions to achieve a minimum global SSR. In this study, we implemented the OLS-CUSUM method, as it does not employ a moving window, and does not perform an iterative regression that requires a training period. The breakpoints routine requires a minimum segment size (h), which makes it not fully automatic. The results reported in this article are produced with the default minimum segment size. Finally, the BIC was used to determine the optimal number of breakpoints for the time series. The temporal and spatial patterns of change points detected by the two methods were statistically analyzed, and evaluated by the 77 manually interpreted validation samples.

The processes of detecting change points
To better show the strategy differences between two methods, we showed the process of change detection through calculating the maximum CUSUM (random localization segmentation) and minimum CUSUM (global segmentation) and the corresponding breakpoints for an NDVI time series of eucalyptus with 306 clear observations ( Figure 4, Table 1 and Table 2). In the random localization segmentation method, the CUSUM was calculated for all possible intervals from 1986 to 2021. The higher CUSUM was mainly distributed in the upper left region (start time (s) was located in the 1st-150th observation and end (e) was in the 251st-306th observation), where the maximum CUSUM was located at the breakpointthe 224th observation in interval 2-306 (Table 1). Therefore, the whole time series was divided into two subsegment intervals, 1-224 and 225-306. For the two subsegments, the CUSUMs of all interval combinations were counted (Figure 4a,b) show the calculation results for all possible intervals, while in the actual procedure a random generation of M subsets was used), and the maximum CUSUM was located at the breakpoint -the 168th observation in interval 152-206. For the span of 1st-224th, its maximum CUSUM also occurred at the 168th observation, but the CUSUM value was smaller than that of interval 152-206, because the 151st observation was a potentially very strong change point. The remaining breakpoints can be detected successively, including the breakpoints of the 151st observation with the maximum CUSUM in interval 52-168 and the 211th observation in interval 170-215 (Table 1).
The OLS-CUSUM of global segmentation presents a pattern similar to maximum CUSUM of the random localization segmentation method (Figure 4c). The larger OLS-CUSUM was in the interval where the starting point was from the 1st to the 150th observation and the ending point was within the 210th-306th observation. There were obvious differences in reaching the minimum OLS CUSUM in different intervals (Figure 4d). However, for similar time periods, they all point to similar breakpoints. The preset numbers of breakpoints showed an important impact on the corresponding minimum OLS CUSUM and breakpoint locations, as the detected breakpoint location was obviously different when the number of breakpoints was different (Table 2). Based on RSS and BIC, the finally determined breakpoint locations were the 123rd, 168th, and 224th observations.

Change detection for different periods of Landsat NDVI
The results of change-point detection based on global segmentation had high variety, while those based on random localization segmentation had more robust results for different NDVI lengths of the same eucalyptus plots ( Figure 5). Random localization segmentation detected four breakpoints, and global segmentation detected three breakpoints in the NDVI time series of 1986-2021. The two methods detected the same two change points that had occurred when the eucalyptus returned to a stable state about one year after logging   Note: "start" and "end" denote start and end points of the intervals in which "cpt" (change-point candidates) were found; column "CUSUM" (cumulative sum) contains a corresponding value of the CUSUM statistic. The numbers in the brackets are marked in the order of 306 time-series observations for convenience, as in Figure 4b,d.
time spans. The detected breakpoints were consistent when using different lengths of NDVI time series, and even detected the same breakpoints in 2005-2021 as in 1986-2021 (Figure 5a). The Chow test was able to remove some spurious breaks. In the NDVI time series of 1986-2021, 1996-03-05 was a spurious break, and two more breakpoints were added (1992-07-16 and 2008-12-11) in the data of 1990-2021. These spurious breaks were successfully removed by the Chow test.

Temporal-spatial patterns of abrupt changes in eucalyptus plantations
Many cut and stability points were omitted in the global segmentation method, which led to very low accuracies for both cut and stability change points (Table 3). For example, the NDVI showed obvious decrease after the cut and recovered rapidly within one year in Figure 6a,b. The global segmentation missed the change points and located an irrelevant breakpoint. The random localization segmentation performed better to capture these land-cover conversions and change processes, and showed higher PA, UA, and OA (Table 3). When the changes were not so frequent or the cut or stability point was not obvious in the NDVI time series, the two methods could reach some degree of agreement (Figure 6a,b). But the global segmentation still missed some change points that were needed to describe the change process properly. Due to the change-detection result of the global segmentation not being robust, it was hard to identify the real generations and ages based on the detected change points. Therefore, the reconstruction of eucalyptus planting history was not further tested.
The spatial distribution of high-and lowfrequency changes of the two algorithms was relatively consistent (Figure 7). The high-frequency change points were mainly distributed in the central area and sides of the central road (a road that passes through the middle of the study area). The highest frequency was distributed in the northern part of the urban area, and the lower frequencies were mainly distributed in the western and northeastern areas. But the number of breakpoints identified by the random localization segmentation was significantly larger than that of the global segmentation in the eucalyptus planting area (Figure 8a), with averages of 3.5 and 2.0 change points, respectively. Random localization segmentation was mainly composed of 2-4 breakpoints, and global segmentation was dominated by 1-3 breakpoints. The proportion of random localization segmentation breakpoints greater than 4 was 26.6%, and global segmentation only accounted for 0.1% (Figure 8a). The area of detected abrupt changes in different times showed that random localization segmentation was able to detect the breakpoints at the beginning and end periods of the NDVI time series (Figure 8b), while the global segmentation did not detect any breakpoint before 1994 or after 2018. These omissions might lead to high uncertainty in eucalyptus disturbance detection, especially for monitoring recent forest changes.

Annual eucalyptus area from 1986 to 2021
The continuous land-cover classification had a satisfactory result with OA of 89.2%. The PA and UA for eucalyptus were 88.2% and 89.6%, respectively. Between 2007 and 2021, the eucalyptus area had increased considerably, accompanied by the emergence of a relatively stable area of logging slash (Figure 9). Whether the logging slash turned to eucalyptus can be determined based on the subsequent  (Table A1 in Appendix). The other forest area increased slightly during this period and decreased in the later stage, mainly due to the strict protection of natural forest and the fact that eucalyptus requires a specific growing environment (especially temperature). Farmers generally do not plant eucalyptus in high mountainous areas where it might be damaged by freezing.

Accuracy assessment of eucalyptus planting history
The OA of the eucalyptus's successive generation was 80.5% (Table 4). The errors were mainly induced by the omission of second-generation eucalyptus. The change-detection algorithm accurately identified the logging, but the corresponding segment was misclassified as non-eucalyptus and resulted in missing one generation. The misclassification of other forest as eucalyptus led to some first-generation samples being labeled as second generation. The eucalyptus age can be accurately estimated with the R 2 value of 0.91 and RMSE of 13.3 months. Misclassification and generation errors have important impacts on the stand age estimation. Three second-generation samples were estimated to be first generation (green circle in Figure 10) and two first-generation samples were mistakenly detected as second generation (red circle in Figure 10). They led to higher RMSE of the estimated age. When the five samples were deleted, R 2 and RMSE became 0.98 and 6.06 months, respectively.

Spatial patterns of eucalyptus planting history and age
The generations in 2021 were mainly first and second generations, with areas of 226.3 km 2 and 91.5 km 2 , accounting for 69.0% and 27.9% of the total area of eucalyptus, respectively. The first-generation eucalyptus was an average age of 6.6 years and was mostly in areas that were expanding outward from the secondgeneration forest area. That pattern shows the obvious expansion trend of eucalyptus in recent years. The second-generation eucalyptus was an average age of 3.3 years and mainly concentrated in the central area of the eucalyptus plantation ( Figure 11a). The rotation cycle of eucalyptus (only pixels that had complete rotation cycles) was mainly 5-8 years, accounting for 66.9% of the total area ( Figure 12a). The area of eucalyptus more than 10 years old accounted for only 2.2%, indicating the management was mainly short-term in the study area.   The vertical lines and their colors represent the order of results (red, green, blue, purple, and cyan represent the detected order -1st, 2nd, 3rd, 4th, 5th -of the breakpoints, respectively). R 1 , R 2 , and R 3 represent round 1, round 2, and round 3 of eucalyptus, respectively. The blue dots are observations, and the black curves are fitted based on harmonic function (Equation 4). A black line shows the change trend of each segment (c in Equation 4). The black line with two arrows indicates the land-cover types.
The spatial pattern of eucalyptus age was relatively fragmented (Figure 11b), with a few patches of younger forest having larger areas. A small number of eucalyptus stands located in the northwestern part of the study area were older (> 140 months), because they are ecological, public welfare forests and protected from logging. The management of eucalyptus in the study area was relatively fragmented, and the cutting boundary was irregular (Figure 11, c-f), which was mainly due to the fact that some of the forestland was contracted by enterprises for commercial management and other parts were managed by local farmers. Some small patches of eucalyptus forest were frequently cut down by individuals. In view of this, our detected result showed high confidence in reconstructing the eucalyptus history, and no postprocessing was performed. The 2021 distribution of eucalyptus age in the study area  presented two-stage characteristics (Figure 12b). The area aged 1-6 years had gradually decreased with a significant increase at 7 years, and then gradually decreased. One-year-old eucalyptus accounted for the highest percentage of 14.5%, followed by 7 years old (12.9%).

The advantages of using random localization segmentation for change detection
The number of breakpoints detected by random localization segmentation for eucalyptus was much greater than that of global segmentation, especially for the proportion of more than four breakpoints. The global segmentation showed high omission and commission errors in detection of both cut and stability change points in eucalyptus. The main reason was that global segmentation needs a certain length of interval between breakpoints to achieve the model fitting. The default parameter h = 0.15 adopted in the global segmentation led to the underestimation of the number of breakpoints in eucalyptus. If the h value were set to a smaller value, such as h = 0.03, the breakpoint monitoring would significantly increase and the detection accuracy of changed points likely would improve. In addition, when the maximum number of breakpoints is changed (the maximal number allowed by h = 0.15 is used by default), the breakpoint locations might also change because changing the parameters will change the interval used for fitting. For example, if the breaks value were set to 10 and h were in default (0.15), a maximum of 5 breakpoints could be detected for the pixel in Table 2. When h is set to 0.03, a maximum of 10 breakpoints could be detected. How to optimize the h value is an important issue of the global segmentation change-detection method. Watts and Laffan (2014)    in their case 0.2, provided some advantages over both smaller and larger values. This parameter is critical to defining the temporal extent of the period during which potential change points can be detected (Dutrieux et al. 2016). A smaller h value might detect some changes that have not been recorded or are unknown in current data. Therefore, a fixed h value is unreasonable due to the variant data density and largely different dynamics in different periods and regions. The flexible h value's definition, such as the localization mechanism, can be adopted in future change detection according to the change types, periods, and regions in order to get more robust temporal segmentations.
In order to avoid inclusion of possible change points in the model initialization, many algorithms strictly restrict the stage by searching for the duration that is as stable as possible. Thus, strict restrictions on the model fitting can lead to a higher omission error (Awty-Carroll et al. 2019;Cohen et al. 2017). That is why some  algorithms (BFAST,Composite2Change [C2C]) have a high UA for the change detection but low PA (DeVries et al. 2015;Hermosilla et al. 2016). Balancing omission and commission errors, and using less restrictive thresholds for finding change are urgently needed for identification of breakpoints. We overcame the issue of the global CUSUM being unsuitable for frequent disturbance and rapid growth of eucalyptus by localizing the CUSUM statistic in a randomized manner. This also dramatically reduced the permitted spacing between neighboring change points in comparison to traditional methods, such as standard BS (Fryzlewicz 2014) and model-fitting-based algorithms. The very fast growth of eucalyptus makes both cutting and stable times act as the change points in the NDVI time series. Traditional global segmentation is often affected by the interaction between the change points. Therefore, there are more uncertainties for eucalyptus forests that operate in short rotation, and Griffiths, Jakimow, and Hostert (2018) stated that approaches utilizing temporal segmentation (Hermosilla et al. 2016;Kennedy, Yang, and Cohen 2010) do not readily capture fast land-use changes, such as swift postdeforestation dynamics or circular slash-and-burn land-use practices. Our results show that the random localization calculation scheme was very robust and consistent at different lengths of NDVI time series and very short spacing between the breakpoints. Moreover, the problem of span or window size selection was avoided by randomly drawing intervals of different lengths. These advantages are very important to obtain consistent change-detection results among different periods and regions. While many forest disturbance methods' and change agents' attributes have been developed over the years (Hermosilla et al. 2016;Kennedy, Yang, and Cohen 2010;Zhu 2017) and some studies paid attention to the disturbance-recovery dynamics in temperate and tropical forests (DeVries et al. 2015;M. Liu et al. 2021), research about the continuous monitoring of generations and rotations in short-rotation forests (such as eucalyptus) is sparse. A study similar to ours was conducted by Dutrieux et al. (2016), who reconstructed land-use history for a swidden agriculture system in Brazil. That system has a very particular rotation cycle wherein the forest is cut and burned for cultivation, and then left to fallow until the next rotation. They used the BFAST method and obtained a number of cultivation cycles with a normalized residual mean squared error of 0.25. They set a small h value (0.07) for the change detection due to the rapid land-use dynamics. The difference in our study is that the rotation of eucalyptus is more diverse because of the different planting methods (coppice regeneration or seedling) and the land use was not changed in two adjacent eucalyptus rotations. We found the global segmentation method cannot accurately find the change points in eucalyptus rotations, mainly because the cut and stability change points are very close. The random localization segmentation proposed here is believed to be more robustness for these shortinterval changes.

Potential use of eucalyptus planting history information
Our study provides a combined framework to obtain planting history, which is a prerequisite for the study of eucalyptus forests' influence on ecosystem dynamics. Eucalyptus has the potential to deliver greater volumes of biomass from the same land area than alternative biomass crops (Cunningham and Tamang 2014). Economically optimal rotations are not frequently employed in plantation management despite the economic importance of eucalyptus to the pulp industry (Diaz-Balteiro and Rodriguez 2006). In our study area, the rotation cycles of eucalyptus are 5-8 years, which are similar to those in Brazil (Rodrigues et al. 2019), but shorter than those in Spain (Diaz-Balteiro and Rodriguez 2006). With the increase of successive generations, the productivity and carbon storage of eucalyptus generally show a downward trend (Wen et al. 2020). Providing wood production and sequestrating more carbon are not completely identical for the eucalyptus industry. Accurate identification of the temporal and spatial distribution characteristics of logging, generation, and rotation of eucalyptus has important applications in determining the optimal area and the optimal forest rotation times to provide carbon sequestration services for targeting carbon neutralization.
The management history was fragmented in the study area compared to the regions managed by forest companies (Deng et al. 2020;Le Maire et al. 2014). This was mainly caused by the complicated local policies and complex terrain. There are no strict restrictions on the cutting age of commercial forest, but the ecological forest is strictly regulated. Therefore, the cultivation mode of large-diameter wood or local small-area cutting might be adopted in ecological forests stands. Validation of stand age at the end of the study illustrated that the detection results of eucalyptus logging information yielded acceptable accuracy using random localization segmentation. However, some uncertainties in the detection of stand ages still exist because there was no cutting time, only stable time, in some pixels (stands). This is because eucalyptus regrows very rapidly in the first year, whether started with seedling or coppice regeneration, and these areas did not have enough data due to clouds. Filling and repairing the contaminated data to enhance the density of the time series will further improve the monitoring results. For example, Landsat-based temporal interpolation and MODIS-Landsat spatiotemporal fusion technology can provide time-series data on daily or other regular temporal intervals (Chen et al. 2021b;Roy and Yan 2020). The Sentinel-1/2 data have greatly improved time-series data frequency (Adrian, Sagan, and Maimaitijiang 2021;Drusch et al. 2012). Integration of Landsat and Sentinel data will improve temporal revisit frequency to every 2-4 days (Griffiths, Nender, and Hostert 2019;Wulder et al. 2015) and yield more timely detection of forest disturbance (Chen et al. 2021a;Shimizu, Ota, and Mizoue 2019). The integration of multiple data will increase the accuracy of eucalyptus planting history reconstruction, and will be our future research.

Limitations and potential improvement of random localization segmentation
Although the proposed random localization segmentation method showed very good performance in our local study, some limitations and potential improvement should be further considered. (1) The method was proposed to find the potential change points in a frequently changing eucalyptus plantation; its generality was not tested in ecosystems that do not change so frequently. Bullock, Woodcock, and Holden (2020) stated that having less restrictive thresholds is necessary to balance omission and commission errors. The more sensitive change-detection method of random localization segmentation with a further Chow test was proved to be a good combination in our study. However, forest changes are very complex, as the disturbance may be caused by storm damage, forest fires, drought stress, insect infestations, disease outbreaks, or harvesting operations. The change magnitude and recovery rate are very different. A higher-sensitivity change-detection method might be more accurate in finding these changes, but also bring more spurious change points in stable forest ecosystems. We believe the proposed method is able to catch these critical changes; particularly, it can provide the order of detected change points. However, it is necessary to further evaluate its generality and accuracy in different ecosystems and with different change agents. (2) The algorithm can be further improved by considering different functions in the Chow test to characterize the different types of land cover and change agents. The harmonic model was widely used in fitting vegetation seasonal changes, but it might overfit the vegetation change processes and fail to capture the rapid recovery trend in the early successional stage. In addition, the land use and cover in different time periods cannot be represented by a fixed harmonic model. (3) Appropriate weights associated with the time-series values can be defined to account for the unequally spaced data and change points that located away from the middle of the subsequence data. The mechanism of random localization segmentation is a CUSUM statistics-based method, and the "variance" of the process will be more comparable at each potential change point by applying weights to the standard CUSUM process. (4) The random localization mechanism also can be used to optimize the current change-detection methods; for example, the local minimum OLS-CUSUM can be calculated in BFAST and used to find the change points in local time-series data.

Conclusions
Time-series remote sensing has become an important means of long-term and continuous monitoring of forest resources. With different data availabilities and quantities in different regions, the development of more robust change-detection algorithms will be an important basis for producing more accurate and comparable forest monitoring data. Based on the idea of random localization segmentation combined with the Chow test, we proposed a more robust and flexible change-detection method and successfully applied it in a eucalyptus plantation to accurately reconstruct its planting history and stand age. Our research result shows that the proposed method is highly robust for different lengths of NDVI time series and can detect critical change points from a forest growth trajectory. The OA of generation identification reached 80.5%, and the eucalyptus age can be accurately estimated with the R 2 value of 0.91 and RMSE of 13.3 months. The temporal and spatial distribution of eucalyptus history information (including generation, rotation cycle, and forest age) can be clearly described. The rotation generations in 2021 were mainly first and second, accounting for 69.0% and 27.9% of the total eucalyptus area, respectively. The operational eucalyptus cycle has been mainly 5-8 years, accounting for 66.9% of the total area. In 2021, one-year-old eucalyptus accounted for the highest percentage of 14.5%, followed by 7-year-old trees (12.9%). These results can be an important input dataset for site quality assessment and forest ecosystem carbon models. They will improve the accuracy of carbon estimation and monitoring, and help to comprehensively assess the contribution to local carbon budget. The proposed change-detection method only needs to set the maximum numbers of breakpoints, avoiding the selection of window span. The results are consistent among different periods of data and also very effective for very short spacing between neighboring change-points. It can be easily applied to other ecosystems and regions with high comparability.