Tracking vegetation degradation and recovery in multiple mining areas in Beijing, China, based on time-series Landsat imagery

ABSTRACT Monitoring of the vegetation change trajectory in mining areas is crucial to understand ecosystem degradation induced by mining activities and to evaluate the effectiveness of ecological restoration. Beijing, the capital of China, is rich in mineral resources and has had many small-scale mine sites. In this study, we aimed to propose a new method to automatically characterize vegetation degradation and recovery trajectory for multiple mine sites. First, the method constructs temporal composites of the Normalized Differential Vegetation Index (NDVI) on a pixel-by-pixel basis using Landsat satellite observations on the Google Earth Engine platform; second, the Pettitt test and Sen+Mann–Kendall analysis are used to identify the year of change and pre- and post-change trends; then, the time-series NDVI is classified into five vegetation trajectory types [recovery ( R ), degradation ( D ), degradation–recovery ( D–R ), recovery–degradation ( R–D ), and no change ( NC )] and 13 subtypes; and finally, the recovery status of the mine sites is analyzed. The method was applied to track vegetation change in 500 mine sites in the Beijing mountainous area and compared to widely used the LandTrendr algorithm. The results showed that our method achieved satisfactory accuracies in trajectory type classification with an overall accuracy of 91.10%. The year of change detected by the Pettitt test was generally consistent with historical Google Earth high resolution imagery. Compared to the LandTrendr algorithm, our method yielded much less omission and commission errors of R, D, R–D, and D–R types and had better capability in identifying gradual recovery or degradation. As the method required few parameters, it is suitable for automatic trajectory monitoring of multiple mine sites. In the study area, 1469.07 ha out of 3746.25 ha of the mine sites have been recovered from 2000 to 2019. The recovery mainly occurred during 2009–2013, around the same time when the Green Mine Construction Campaign was launched in Beijing. 622.62 ha of the mine sites have experienced degradation probably due to illegal mining activities. Our results are expected to support evaluation of mine restoration effects and detection of illegal mine sites.


Introduction
Mining has been and continues to be an essential economic component to satisfy people's demand and industrial development. However, mining activities have significant impact on the surrounding environment and ecosystem health as the land cover change by mining may cause land degradation, vegetation destruction, and biodiversity loss, posing great challenge to sustainable mining development. In recent years, many countries and regions have implemented a series of reclamation and ecological restoration measures to achieve low-impact or sustainable mining goals (Moomen et al. 2019). However, there are also some places where illegal mining activities still exist, causing continuous ecosystem degradation. In forested area, both mining activities and vegetation restoration result in abrupt or continuous changes in the vegetation cover. Timely monitoring of vegetation change is particularly important to understand the impacts of mining activities on the surrounding ecosystem and to evaluate the effects of ecosystem recovery measures, which further supports regional sustainable development (UNDP and UN Environment 2018).
With the advancement of satellite observation technology, remote sensing has become an effective tool providing timely and accurate data for land surface change monitoring over large area. Previous research has utilized satellite imageries acquired using a Moderate Resolution Imaging Spectroradiometer (MODIS), Landsat-series satellite sensors, and commercial high-spatial-resolution satellites to monitor vegetation change in mining areas (Erener 2011;Yu et al. 2018;Yang et al. 2018;Wang et al. 2021). Although MODIS provides time series imagery with high revisiting frequency (1 day), its low spatial resolution (≥ 250 m) hinders the application in vegetation change detection in small mine sites (Ma et al. 2011;Lei, Ren, and Bian 2016). High-spatial resolution satellite images (< 5 m), such as SPOT or WorldView-2/3, are more beneficial for small and medium scale change detection (Wu, Ma, and Liu 2009;Zhang et al. 2020), but the high cost of these commercial satellite images limit their wide applications in long-term change analysis. As Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper (ETM+), and Landsat 8 Operational Land Imager (OLI) provide freely available time-series imagery with 30 m spatial resolution acquired over the past 30 years, Landsat images have been increasingly used to identify vegetation changes in mining areas (Pericak et al. 2018;Mi et al. 2019;Pouliot and Latifovic .2016).
The existing methods for monitoring vegetation disturbance and restoration in mining areas can be categorized into two types: (1) classification method and (2) time-series analysis method. The former methods classify multi-temporal remote sensing images into land cover classes, and the temporal changes of vegetation loss and vegetation gain are analyzed based on the land cover maps (Grinand et al. 2013;Petropoulos et al., 2014;Vasuki et al. 2019;Chen et al. 2021). The second type of methods aims to detect change and identifies change types based on timeseries spectral indices such as normalized differential vegetation indices (NDVI). Typical approaches include threshold analysis (Li, Feng, and Qin 2020), linear or non-linear regression (Liu, Zhou, and Bai 2016;Vidal-Macua et al. 2020), trend analysis (Liu et al. 2019;Wang et al. 2021), and temporal segmentation algorithm such as Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) Xiao et al. 2020). Li, Feng, and Qin (2020) identified six types of land damage and restoration in rare earth mining areas by setting NDVI thresholds based on Landsat images acquired in growing season from 1995 to 2017. Vidal-Macua et al. (2020) utilized the linear regression model to evaluate the trend of MODIS NDVI after reclamation of mining sites. Wang et al. (2021) applied Mann-Kendall trend analysis to monitor the monotonous trend of vegetation restoration. Generally, these research studies did not either characterize complex vegetation change trajectory types such as transition from vegetation degradation to restoration or from restoration to degradation or automatically identify the time stamp when vegetation degradation or recovery started.
The LandTrendr algorithm was developed by Kennedy, Yang, and Cohen (2010). It implements a temporal segmentation approach to divide timeseries NDVI into segments. The change events and the associated attributes such as year of the event and duration are further labeled. The LandTrendr algorithm has been widely used to characterize forest disturbance Xiao et al. 2020). Yang et al. (2018) adopted the LandTrendr algorithm to assess vegetation disturbance and recovery in Curragh coal mine in Australia. As LandTrendr requires parameters setting specific to a single mining site, its capability in automatic change monitoring over multiple mining sites with complex change trajectories remains unclear.
In China, the mining industry has played a significant role in social and economic development over the past few decades. Beijing, the Capital of China, is rich in coal, metallic minerals (such as iron ore), and nonmetallic minerals (such as limestone, dolomite, and marble). Mineral mining in Beijing was mainly conducted through surface mining. Years of mining activities have led to serious vegetation destruction, soil erosion, ecological damage, and environmental pollution around mining sites (Liu 2014). In order to protect the ecological environment, the Beijing Municipal Government has taken a series of measures to restrict the continued exploitation of minerals and reduce the total number of active mining sites since early 2000s. Reclamation and ecological restoration measures have also been implemented gradually (Gan 2007;Li and Zhang 2009). However, reports showed that illegal mining activities still existed (Jiang, Tian, and Zhan 2017). Some abandoned mine sites were experiencing continued degradation due to soil erosion. It is in urgent need to characterize the vegetation change trajectories in these mine sites to support the assessment of ecological restoration effects and provide insight into improving future reclamation actions.
In this study, we aim to present a new method to track vegetation degradation and recovery trajectory in multiple mine sites. The new method is designed that few parameters are required and has the capability to detect vegetation change trajectory types automatically for a large number of mining sites. The method is implemented and evaluated in the 500 mine sites in Beijing mountainous area, aiming to provide a scientific baseline for the formulation of ecological restoration and management policies and to support sustainable development of the region.

Study area
The study area is located in the northern mountainous region of Beijing (115°25′18"-117°30′23.41 "E, 39°49′12.77"-41°3′24.84 "N) with a total area of 14,803.40 km 2 ( Figure 1). The study area has a typical warm, temperate, semi-humid continental monsoon climate with an average annual temperature of 10°C-12°C and annual precipitation ranging from 400 mm to 1000 mm. The elevation ranges from 43 m to 1191 m. Vegetation is dominated by shrubs. Coniferous trees also existed. The distribution of mineral resources is relatively concentrated, such as limestone and metal minerals in Yanqing District, iron ore in Miyun District, and coal mines in Mentougou District. The study area contains a total of 500 mine sites, with an average area of 7.5 ha per mine site. The types of minerals include coal, metallic, and nonmetallic minerals. The General Planning of Beijing Mineral Resources (2008Resources ( -2015 reported that most of these mine sites were developed before 2000, and the intensive mining activities have aggravated the degradation of regional ecological environment and the occurrence of geological disasters (Jiang, Tian, and Zhan 2017). In response, Beijing started ecological restoration and treatment over the mine sites (Liu 2012). Natural restoration and artificial reclamation were the main methods for ecological restoration of the abandoned mines. While the former process was exceptionally slow and did not match the required effect of vegetation recovery, the latter one greatly accelerated the speed of ecological restoration (Li and Zhang 2009).

Datasets and pre-processing
In this study, we collected all available Landsat TM/ ETM+/OLI Level 2 surface reflectance (SR) image collections acquired during the growing season (May to September) in 2000-2019 on the Google Earth Engine (GEE) platform ( Figure 2). The total number of images was 898 covering 3-4 Landsat scenes (Path122-124, Row 42-43, Figure 1). For each image, the quality assurance (QA) bands of the SR products were used to mask out the pixels covered by cloud, cirrus, shadow, and snow. Figure 2a shows that the number of Landsat images varied from 2000 to 2019 ranging from 30 to 60 scenes. The exact numbers of Landsat images are listed in Table S1. Figure 2b shows the percentage of pixels with valid observations, i.e. observations not contaminated by clouds, shadows, snows, ice, etc. Over 50% pixels had more than six valid observations in the growing season in 2000-2012, and over 50% pixels had more than 11 valid observations in 2013-2019.
Google Earth's high-spatial resolution images are commonly used as auxiliary data for visual interpretation. Google Earth provides images taken in different periods over many years with high spatial resolution (< 5 m), and detailed land cover within the mine sites such as vegetation and mine dumps can be visually observed in these images. In this paper, the historical Google Earth high-spatial resolution images were used as reference data for accuracy assessment of vegetation change trajectory characterization.

Methods
The workflow of this study is outlined in Figure 3. First, temporal statistical composites of NDVI were calculated using time-series Landsat imagery on a pixel-by-pixel basis (Section 3.1) based on which the unchanged pixels were detected, and the changed pixels were selected for further processing (Section 3.2). Then, Savitzky-Golay filter was applied to annual mean NDVI time series to remove the spikes and smoothen the NDVI trajectory (Section 3.3). Third, the Pettitt test and Sen +Mann-Kendall trend analysis were implemented, and the year of change and the NDVI trajectory types were identified to characterize vegetation degradation and recovery within the 500 mine sites (Section 3.4). Finally, the results were evaluated using reference samples visually interpreted based on historical Google Earth high-spatial resolution imagery, and our method was compared to the LandTrendr approach (Section 3.5).

Developing temporal composites of NDVI
For each image, NDVI was calculated using the red band (Band 3 for Landsat 5/7 and Band 4 for Landsat 8) and near-infrared band (Band 4 for Landsat 5/7 and Band 5 for Landsat 8). Previous research has reported good consistencies of vegetation indices between Landsat TM/ETM+ and Landsat 8 OLI regardless of slightly different NIR waveband designation (Ke et al. 2015;Roy et al. 2016). For each pixel, the mean NDVI during May-September in each year was calculated (NDVI y mean ) and the time-series NDVI y mean were sorted in chronological order. Thus, a total of 20 scenes of yearly NDVI y mean from 2000-2019 were obtained. Based on the 20 NDVI y mean layers, the maximum NDVI (NDVI max ) and the minimum NDVI (NDVI min ) were calculated.

Detecting unchanged pixels
Our field surveys and careful examination of the historical Google Earth high-spatial resolution imagery showed that some area in the mine sites were consistently vegetated or non-vegetated land during the past 20 years. The stable non-vegetated lands were covered by either bare soil, buildings or mineral deposits and were considered as unrecovered land. In this study, we extracted the stable unchanged pixels using NDVI thresholds determined by examining NDVI y mean values of each type of pixels, i.e. stable vegetated (SV) pixels, stable non-vegetated (SNV) pixels, and pixels with land cover change. Here, we manually selected 100 pixels from each of the three types that were visually interpreted from Google Earth images. Figure 4 illustrates that the SNV pixels had consistently low NDVI y mean , and the maximum value of NDVI y mean (NDVI max ) was less than 0.2. Likewise, the SV pixels had consistently high NDVI y mean , and the minimum value of NDVI y mean (NDVI min ) was greater than 0.35. Therefore, we used the thresholding rules NDVI max < 0.2 to identify SNV pixels and used NDVI min > 0.35 to identify SV pixels. The remaining pixels were considered as changing pixels during 2000-2019 and selected for further processing.

Smoothing NDVI time series with Savitzky-Golay filter
In order to suppress the influence of noise and outliers in the time series of NDVI y mean , we adopted the Savitzky-Golay (S-G) filter to smoothen the NDVI y mean trajectory on a pixel-bypixel basis. The S-G filter is a weighted moving average convolution with the convolutional coefficients determined by a polynomial least-square-fit within the sliding window (Savitzky and Golay 1964;Chen et al. 2004). Based on the S-G filtering principle, the S-G filtering process for the NDVI time series dataset can be described by the following equation: where Y is the original NDVI value, Y � is the resultant NDVI value, C i is the coefficient for the ith NDVI value of the filter (smoothing window), and N is the number of convoluting integers and is equal to the smoothing window size (2 m + 1). The index j is the running index of the original ordinate data table. The smoothing array (filter size) consists of 2 m + 1 points, where m is the half-width of the smoothing window. In this study, the window size was set to 5 for timeseries NDVI y mean , and a third order polynomial function was used.

Change point detection using Pettitt test
Time-series NDVI y mean shows decreasing trend after disturbance by mining and increasing trend after vegetation recovery. The rate and magnitude of NDVI y mean decline are affected by the manner of mining and the degree of disturbance; the rate and magnitude of NDVI y mean increase depend on the recovery status and effectiveness of restoration. In this study, some area in the mine sites have experienced continued NDVI y mean increase (or decrease) during 2000-2019, indicating consistent vegetation recovery (or degradation). However, some areas have experienced opposite trends before and after a certain year (as illustrated in Figure 5), indicating changing status from disturbance (or recovery) to recovery (or disturbance).
In this study, we applied the Pettitt test to identify the year of change in the time-series NDVI y mean . The Pettitt test is the Mann-Whitney nonparametric test. It determines the change point of the sequence by examining the longterm sequence changes (Pettitt 1979). The Pettitt test has been used widely in hydro-climatological studies to detect changes in long-term precipitation, runoff, and temperature (Duhan and Pandey 2013;Pandey and Khare 2018;Shawul and Chakma 2020). The calculation formulas are as follows: (3) where U t;n is the Mann-Whitney nonparametric statistics. x i is the ith value of the point in timeseries NDVI y mean . n is the length of the time-series NDVI y mean , which is 20 in this study. sgn x t À x i ð Þ is the sign function to determine the parameter x t À x i . K is the change point derived from U t;n , and K represents the change point, i.e. the year of change. P is the result of the significance test for the change point. In this study, the significance parameter P is set to 0.1. When P ≤ 0.1, it is considered that change point exists in time-series NDVI y mean . If multiple change points are detected in the time series, the first one is considered. In this study, less than 1% of the pixels have two change points.

Trend analysis of NDVI time series
After performing the Pettitt test on the NDVI y mean time series, the pixels can be classified as those with change point (pass Pettitt test) and those without change point (fail to pass Pettitt test). For the pixels where the change point was detected, Sen + Mann-Kendall trend analysis and NDVI threshold were applied to the NDVI y mean time series before and after the change point, respectively ( Figure 5); For the pixels where the change point was not detected, Sen + Mann-Kendall trend analysis and NDVI threshold were performed on the entire NDVI y mean time series.
In this study, two non-parametric tests, Mann-Kendall significance test and Sen's slope trend test, were used to detect trends in NDVI y mean time series (Mann 1945;Sen 1968;Kendall 1970). The Mann-Kendall significance test is a non-parametric test, which is usually used to test the significance of the trend (Khan et al. 2021;Valipour et al. 2020). The Sen's slope trend test was used to detect trends of time series. The formulas of the two non-parametric tests are as follows: where Z C is the standardized test statistic; S is the test statistic of the time-series NDVI y mean . x k and x i are the kth and ith values in time-series NDVI y mean . n is the length of the time-series NDVI y mean , which is 20 in this study. sgn x k À x i ð Þ is the sign function to determine the parameter x k À x i . var S ð Þ is the variance of S, and β is the trend degree of the time-series NDVI y mean , which is used to identify the upward or downward trend of the NDVI change in time series. In this study, the significance level α ¼ 0:05 and Z=1.96. |Z|> 1.96 indicates a significant change in NDVI time series; |Z|< 1.96 indicates an insignificant change in the NDVI time series. A positive value of β for the Sen's slope test indicates a rising trend and a negative value of β of the Sen's slope test indicates a downward trend. Based on the β and |Z| values derived from Sen + Mann-Kendall analysis, the trend of NDVI time series can be divided into one of the following three categories (Table 1): (1) Increasing trend, denoted as "+." A NDVI trajectory shows an increasing trend when β>0 and |Z| >1.96 or when β>0 and |Z|≤1.96, but the difference between NDVI at the end of the time series (NDVI end ) and that at the start of the time series (NDVI start ) is greater than 0.3.
(3) No change, denoted as "0". A NDVI trajectory shows a flat trend when |Z|≤1.96, indicating that NDVI does not change significantly.

Classification of trajectory types
The entire NDVI trajectories were further classified into five categories based on the Pettitt test and trend analysis as follows: "Recovery," "Degradation-Recovery," "Degradation," "Recovery-Degradation," and "No Change" (Table 2).
Type 1: Recovery (denoted as "R", Figure 6a). The NDVI time series of the pixel shows an overall upward trend during the study period. This type mainly includes (1) continuous increase of NDVI when no change point is detected (denoted as '+'), (2) flat pattern of NDVI (no change) before the change point and upward trend after the change point ("0+"), (3) increase before the change point and no change after the change point ("+0"), and (4) different magnitudes of increase before and after the change point ("++").
Type 2: Degradation-Recovery (denoted as "D-R", Figure 6b). The change point is detected during the study period, and the NDVI time series shows a decreasing trend before the change point and increasing trend after the change point ("-+").
Type 3: Degradation (denoted as "D", Figure 6c). The NDVI time series shows an overall downward trend. This type mainly includes (1) continuous decrease of NDVI when no change point is detected ('-'), (2) no change before the change point and decrease after the change point ("0-"), (3) decrease before the change point and no change after the change point ("-0"), and (4) different magnitudes of decrease before and after the change point ("-").
Type 4: Recovery-Degradation (denoted as "R-D", Figure 6d). The change point is detected during the study period, and the NDVI time series shows an upward trend before change and downward trend after change point ("+-").
Type 5: No Change (denoted as "NC," Figure 6e). The NDVI time series basically do not change significantly throughout the studied time period, maintaining a stable status. The SV and SNV pixels are also included in this type.
In order to evaluate the restoration status of the mine sites, we calculated the percentage of pixels with R, D-R types, and SV subtype within each mine site (Perc Rec ). We defined the mine site as not recovered, partially recovered, and fully recovered if Perc Rec < 25%, 25% � Perc Rec < 75%, and Perc Rec � 75%, respectively.

Validation of the results and comparison with LandTrendr algorithm
A total of 1000 sample pixels were randomly selected within the mine sites using the stratified sampling strategy (locations shown in Figure S1). The vegetation degradation and recovery type of each pixel were were discarded or replaced with new samples with high confidence. The visually interpreted trajectory types were compared to our classification results to derive confusion matrix. Overall accuracy, producer's accuracy, user's accuracy, and Kappa coefficients were calculated for accuracy assessment. Note that the year of change derived from the method was not evaluated quantitatively because many sample pixels did not have full collection of yearly high-resolution imagery on Google Earth. For comparison purpose, we implemented LandTrendr over the 500 mine sites. LandTrendr utilized linear regression and point-to-point fitting of time-series spectral indices to identify features such as segments, homogeneous periods of change or stability, and occurrence of events (Kennedy, Yang, and Cohen 2010). It was originally developed for forest disturbance monitoring. In this study, we used GEE-version LandTrendr (Kennedy et al. 2018). For consistency, the input datasets of LandTrendr were Landsat imagery during May to September. NDVI gain and NDVI loss in time series were identified separately using "runLT" and "getChangeMap" functions. If more than one gains (or losses) were detected, the first event was identified. The five vegetation trajectory types (R, D, D-R, R-D, and NC) were determined using the rule sets in Table S2. LandTrendr requires a set of control parameters to ensure the quality of change detection. The optimal parameters were selected by testing different combinations of parameter values (Table 3, Table S3). The classification results of the LandTrendr algorithm were also evaluated and compared to those from our method. Note that the definition of the occurrence of an event is different from that of the change point in our method. If both disturbance and recovery existed in one time series, the occurrence of disturbance (or recovery) was defined as the start year of disturbance (or recovery) in LandTrendr (as illustrated in Figure S2). In our method, the year of change was defined as the breakpoint year where NDVI started to show the opposite trend.

Accuracy assessment and comparison with LandTrendr algorithm
Tables 4 and 5 list the confusion matrices of our method and LandTrendr algorithm, respectively. It is obvious that our method yielded significantly higher The minimum number of valid observations required to perform output fitting.
Despike 0.9 A threshold used to dampen the spikes.
Note: the full meaning of the parameter above can be found in (Kennedy, Yang, and Cohen 2010) and https://emapr.github.io/LT-GEE/lt-gee-requirements.html (Kennedy et al., 2018).   Figure 7 demonstrates the resultant maps of vegetation trajectory types in an example mine site using our method and LandTrendr algorithm. Overall, both methods show similar spatial patterns of recovery distribution (R type), while our method detected more pixels with D-R types. From Google Earth imagery, it can be seen that a forest patch in the northwest of the mine site (Figure 7a1, Point 1) turned to bare land in 2009 ( Figure 7a2) and remained un-recovered by 2015. The NDVI trajectory type was determined as D (subtype "-0") ( Figure 7b2-b4), and the change point detected by our method is around the year 2006, which is consistent with Google Earth imagery. In comparison, LandTrendr failed to detect vegetation disturbance in this area. Point 2 had relatively low vegetation coverage in 2002 (Figure 7a1). The vegetation was completely removed by 2005 and recovered by 2015. This area was identified as D-R type (Figure 7b4) and started to recover around 2010 (Figures 7b1 and 7d2). The recovery year detected by LandTrendr was similar to those detected by our method; however, LandTrendr again failed to detect the decreasing trend of NDVI before the recovery year, thus identifying as R types (Figure 7c3 was identified as R type (subtype "0+") and started to recover around 2011 (Figure 7b1-b4) by our method. For this point, LandTrendr generated the same trajectory type (Figure 7c3) and similar recovery year (2012, Figure c2) as our method.
The area estimates of vegetation trajectory types show large discrepancies between our method and LandTrendr (Table 6). Our method estimated a much larger area with recovery and degradation than LandTrendr (1345.41 ha vs. 456.93 ha for R type and 509.76 ha vs. 405.72 ha for D type). It also detected a considerable number of pixels with changing trend from degradation to recovery (D-R type, 123.66 ha) and from recovery to degradation (R-D type, 112.86 ha). Pixels with stable status accounted for 44.17% of the study area. In comparison, pixels with no change event are dominant type from the LandTrendr algorithm (76.25%). Figure 8 demonstrates spatial distribution of vegetation trajectory types and the year of change in two example mine sites. Mining activities started before 2000 in both sites, while they exhibited different patterns of trajectory. Site 1 (Figure 8a1-a10) mainly showed degradation pattern. In particular, the north part of site 1 experienced vegetation destruction after recovery. As illustrated in the Google Earth imagery, green vegetation coverage increased from 2011 to 2015, while the vegetation patches were destroyed during 2015 to 2017. The year of degradation detected by our method was around 2015, and the type was R-D, which is consistent as shown in Google Earth imagery. Site 2 (Figure 8b1-b10) is interesting as  both recovery and degradation existed. The Google Earth imagery showed obvious increasing tree coverage since 2010 in the east part of site 2 and destruction of vegetation in the west and south. Trajectory types identified by our method generally show R type in the east and D type in the west. These examples illustrate that the disturbance and restoration patterns were complex in the study area, and the trajectory types show considerable spatial variations among the mine sites and within the mine sites. Table 7 lists the area of each type and sub-types for the 500 mine sites. Over 55% of the area experienced vegetation change during the study period. The remaining area with NC type includes stable vegetation (406.17 ha), stable non-vegetation (93.33 ha), and those failed to pass Pettitt test (1155.06 ha). For the R type (1345.41 ha), over 50% area had time series NDVI transitioned from flat to upward trend (subtype "0+"), indicating vegetation recovery started sometime during the study period. There was a substantial area with "+ +" subtype (269.82 ha), indicating different magnitude of NDVI increase between pre-change and postchange period. This might be associated with different stages of ecological restoration in mine sites. Around 40% of the area with D type had "0-" subtype (216.27 ha). In addition, there was a considerable area with status transiting from recovery to degradation (112.86 ha with R-D type). The "0-" subtype and R-D type maybe associated with vegetation disturbance due to emerging mining activities during the study period, and some of the restored area was at the risk of degradation. Overall, vegetation restoration and disturbance coexisted in the study area from 2000 to 2019, while the area with recovery trend (both R type and D-R type) was significantly larger than that with degradation trend (both D type and R-D type). Figure 9 showed the distribution of change point for "0+" subtype of R, D-R type, "0-" subtype of D, and R-D type. The change point of the former two types indicates the year when recovery started, while that of the latter two types indicates the year when degradation started. For the "0+" subtype of R type, the recovery start years concentrated between 2009 and 2012. The area is 421.38 ha, accounting for 58.65% of all pixels with "0+" subtype. Similarly, for D-R type, the recovery start years concentrated between 2011 and 2013, with an area of 85.41 ha (69.07% of D-R type). The coal mine sites recovered slightly earlier than metallic and nonmetallic sites. Interestingly, degradation start years seems to coincide with recovery start years, mainly distributed from 2009 to 2011 ( Figure 9c). The area with R-D type was dominated by nonmetallic mine sites, and they are at the risk of degradation mainly in 2014, later than the change year of other types.

Recovery status of mine sites
According to the percentage of the recovered area within the mine sites, we classified the mine sites into un-recovered, partially recovered, and fully recovered. Figure 10 and Figure 11 show that 461 out of the 500 mine sites have been partially or fully recovered, accounting for 92% of the total number of mine sites. The remaining 39 mine sites (23 metallic mine sites and 16 nonmetallic mine sites, Figure 11) have less than 25% pixels covered by vegetation, indicating that they should be paid close attention and should be the key mine sites for restoration in the future. The spatial distribution of the recovery types ( Figure 10) presented obvious clustering pattern in that fully recovered mine sites are mainly located in Mentougou District, while un-recovered mine sites are mainly located in Miyun and Changping Districts (13 in Miyun District and 9 in Huairou District). Around 78% of the coal mine sites have been fully recovered and 22% of the coal mine sites have been partially recovered. As almost all coal mine sites were within Mentougou District    ( Figure 1), this indicates that the ecological restoration measures taken by the local district have been generally effective.

Advantages of the method in tracking vegetation changes in mining areas
The above results showed better accuracies of our method in identifying different types of vegetation change trajectory compared to LandTrendr. With respect to the classification results, our method yielded low omission and commission errors for each of the vegetation trajectory types. In comparison, LandTrendr resulted in greater underestimation of the area with R, R-D, D, and D-R types and greater overestimation of NC area (Table 5 and  (Rodman et al. 2021) examined the performances of LandTrendr in detecting forest disturbances induced by wildfire and spruce beetle and reported that LandTrendr significantly underestimated the area affected by gradual spruce beetleinduced mortality and was difficult to detect lowseverity disturbance. This is consistent with our findings.
In addition, vegetation change in different stages of ecological restoration can be detected by our method. Figure 12 shows an example of a nonmetallic mine site from mine development in 2002 (Figure 12a) to land surface reshape, limestone stope removal, and gradual vegetation recovery since 2010 (Figure 12b), and to replanting of trees in 2017 and 2019 (Figures 12c and 12d). Correspondingly, our method detected two-stage increase of NDVI ("++" subtype) with the change point occurring in 2010 (Figure 12e). The area witnessed more rapid increment of NDVI after 2010 than that before 2010.
Another advantage of the method is that it requires few parameters. Thus, it is possible to perform fully automatic identification of trajectory types for multiple mine sites. Compared to previous research utilizing cloud-free image scenes to detect inter-annual NDVI change Mi et al. 2019;Li, Feng, and Qin 2020), our method made use of the GEE platform and adopted a pixel-based strategy based on all Landsat images during the growing season. In Beijing mountainous area, over 70% of the annual precipitation concentrates in summer. The pixel-based strategy guarantees sufficient valid observations over each pixel and minimizes the inter-annual fluctuation of NDVI induced by different temporal distribution of images. It is expected that our method will perform well in other mine sites with similar trajectory types and climate conditions. In different climate regions, growing seasons need to be examined to guide selection of appropriate months.

Effects of restoration policies on mining sites recovery
Our results showed that most mine sites have been partially or fully recovered by 2020, confirming the effectiveness of mine restoration measures taken in Beijing. According to the Beijing Municipal Bureau of land and resources (2006) Our results reported that some mine areas demonstrated slightly increasing NDVI trend before 2010 ("+," "++" and "+0" subtypes), indicating natural and slow vegetation succession, and most area with recovery type witnessed sharp increase in NDVI during 2011 and 2012 after the GMCC was launched (Figures 9a and 9b).
The effectiveness of the mine restoration measures varied among the mineral resource types and in different districts. All coal mines have been partially or fully recovered, while some metallic and nonmetallic mineral mine sites have not been recovered yet. As most coal mine sites were in Mentougou District, this district witnessed better recovery status than other districts. The variations in the recovery status can be partially explained by topographic factor. Nonmetallic and metallic mine sites are mostly distributed in the mountains with high elevation and steep slope, where vegetation is relatively difficult to recover. Figure 13 clearly illustrated that the area with flat terrain or gentle slope is associated with vegetation recovery (R or D-R type), while the area with steep slope showed a degrading trend of vegetation. For the mine sites with gentle slope, reclamation was implemented by stabilizing the slope, creating a surface soil environment, and then planting suitable vegetation. Steep slope increases the soil erosion risk and is unfavorable for soil stability and water supply to plants (Vidal-Macua et al. 2020).
Another reason that explained the variations in recovery status is the strength of the restoration measures in different districts. Mentougou District was the pilot district that first carried out ecological restoration for abandoned mines in 2006. Active and strong restoration measures were implemented from 2008 to 2010. However, in other districts, illegal and private mining activities continued in these years despite repeated prohibitions, mainly in the form of cross- border mining or continued mining after the expiration of mining rights. Previous research reported that by the end of 2015, there were still 35 illegal private mining areas in Miyun District, Mentougou District, and Fangshan District (note that Fangshan District is not in our study area). Especially, the mining areas in Miyun District showed severe vegetation degradation and environment pollution (Jiang, Tian, and Zhan 2017). This is consistent with our findings, which showed that over 22.81% of the mine sites in Miyun District have experienced vegetation degradation ( Figure 10). Variations in vegetation recovery status in the study area indicate that mine sites restoration and management are still challenging in Beijing. Supervision of the activities in the mine sites need to be strengthened. As the automatic method developed in our study has the ability for timely monitoring of vegetation change over large area, it has great potential to help environmental agencies detect and identify illegal mining activities.

Uncertainties and implications for future work
Although our method yielded positive performance in identifying vegetation change trajectory types and detecting the year of change in the 500 mine sites in Beijing mountainous area, we note several limitations in our research. First, the current method is suitable for surface mining assessment, while underground mining may not have a great impact on surface vegetation (Bao et al. 2014). Second, the current method only considers five trajectory types due to the short time span (20 years) because the reclamation measures were not taken until 2006 in the study area. More complex vegetation trajectory types (e.g. multiple times of vegetation disturbance and recovery events) may exist during a longer time period such as 30 or 40 years. In this case, trajectory classification system considering multiple change points (e.g. recoverydegradation-recovery, degradation-recoverydegradation) need to be considered. As the Pettitt test has the ability to detect multiple breakpoints in time series, our future work will evaluate the performance of this method in identifying complex vegetation trajectory types. We will also evaluate the method in study areas in different climate regions. When implementing the method in other areas, we recommend that the vegetation trajectory classification system should be determined with the guidance of prior knowledge of the mining/restoration activities as vegetation change types may differ from one study area to another Vidal-Macua et al. 2020).

Conclusions
In this study, we presented a new approach for characterizing vegetation degradation and recovery trajectory for multiple mine sites. Based on time-series Landsat satellite imagery, the method integrates Pettitt test and Sen+Mann-Kendall trend analysis to classify five vegetation trajectory types including 13 subtypes and to identify the change point during the trajectory. Using 500 mining areas in Beijing mountainous area as an example, the results show that the overall classification accuracy reached 91.10%, and all five types achieved high user's accuracies and producer's accuracies (>85%). The year of change detected by the method is generally consistent with historical Google Earth high resolution imagery. Compared to the LandTrendr algorithm, our method yielded much less omission and commission errors and has better capability in detecting gradual change. In addition, our method requires few parameters so that it is suitable to monitor multiple mine sites.
In the study area, an area of 1469.07 ha out of 3746.25 ha of the mine sites have been recovered from 2000 to 2019, including 1345.41 ha with R type and 123.66 ha with D-R type. The recovery mainly occurred during 2009-2013, around the same time when the Green Mine Construction campaign was launched in Beijing (2010), indicating the effectiveness of the mine restoration measures. 622.62 ha of the mine sites have experienced degradation (D type and R-D type), suggesting existence of illegal mining activities and risk of vegetation degradation. The recovery status of mine sites showed spatial clustering pattern, with fully recovered mine sites mainly located in Mentougou District and un-recovered mine sites mainly located in Miyun and Changping Districts. Overall, the approach presented in this study provided reliable results for recovery and degradation process monitoring in the mine sites in Beijing mountainous area. It is expected that this study can support regional mine restoration and sustainable mining development planning. Future work will include implementation and evaluation of the method in other study area or over longer periods with more complex vegetation trajectory types.