Forest disturbance decreased in China from 1986 to 2020 despite regional variations
Disturbance detectionWe used a well-established spectral-temporal segmentation method, Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr), to detect disturbances within the Google Earth Engine (GEE) cloud-computing platform57,58. The core of the LandTrendr is to extract a set of disturbance-related metrics by breaking pixel-level annual time-series spectral trajectories into linear features using Landsat observations. The LandTrendr has been widely used for change detection in various forest settings, and details about the algorithms can be found in previous publications57. Here we briefly described the key steps in generating the year and type of disturbances in China’s forests using the LandTrendr within the GEE platform. The overall analytic flow can be found in Supplementary Fig. 10.First, we generated annual spectrally consistent time-series data by using all available, good quality (cloud cover ≤ 20) Tier 1 Landsat 5 (Thematic Mapper), Landsat 7 (Enhanced Thematic Mapper Plus), and Landsat 8 (Operational Land Imager) images acquired during the peak growing seasons (June 1—September 30) from 1986 to 2020. The peak growing seasons were selected to exclude compounding influences from ice, snow, and soil, and to maximize the spectral changes after forest disturbances. To tackle the spectral inconsistency among Landsat sensors, we harmonized spectral values via linear transformations according to band-respective coefficients presented in59. Clouds, cloud shadows, snow, and water were masked out using the Fmask algorithm60. The annual band composites at 30-meter spatial resolution during 1986–2020 were computed using the Medoid method61.Secondly, we ran the LandTrendr using five spectral indices, including two spectral bands (shortwave infrared I and II that were B5 and B7), tasseled cap wetness (TCW), normalized burn ratio (NBR), and normalized difference vegetation index. These five indices were effective indictors to represent vegetation greenness and structures, and were commonly used for detecting changes in forest disturbance and recovery62. For each spectral index, the LandTrendr produced a set of parameters to describe a possible disturbance event at the pixel level, including spectral values at pre-disturbance level (preval), magnitude of change (mag), duration (dur) and rate of change (rate), and the signal-to-noise ratio (dsnr) (n = 5). Using these five spectral indices, we generated a stack of disturbance-related parameter layers (n = 25, 5 spectral indices × 5 parameters), which were later used to detect and classify disturbances using machine learning models derived from reference data (described below).Disturbance classificationReference dataHigh-quality consistent reference data is key to train and classify disturbance types. To do so, we generated a total of 31225 reference points using a hierarchical approach. We first generated a large number of potential disturbance points using forest loss data from 2001 to 20203. Then we separated fire disturbances from non-fire disturbances by overlaying MODIS burned area (BA) with potential disturbance points following the procedure used by63. Specifically, fire disturbances were determined if the MODIS BA data coincided with the Landsat-derived forest loss for the fire year and 2 years postfire (i.e., t + 0, t + 1, t + 2) to account for delayed post-fire tree mortality. Following this step, we derived points as potential disturbances that consisted of fires and non-fire disturbances (including forest conversion to other land use types and silvicultural practices at various intensities). We also generated roughly the same number of points that experienced no disturbances (e.g., persistent forests), which were determined by selecting pixels with very few changes in spectral indices. These reference points, including fire, non-fire disturbances, and persistent forests, were then used to sample the time-series spectral data from 1986 to 2020. Finally, time-series spectral data from each reference point were visually checked to make sure they accurately represented disturbance events. This process resulted into a total of 31225 reference data points, including 2356 fire disturbance points, 13,242 non-fire disturbance points, and 15,627 no disturbance points (persistent forests) (Supplementary Fig. 2).Random forest classificationWe used machine learning modeling to classify each pixel into fire disturbance, non-fire disturbance, or no disturbance. The reference data points were used to sample the LandTrendr-derived disturbance-related parameter layers described above, which resulted into a dataset consisting of disturbance types. We divided the dataset into 70% of training data, and 30% as validation data. Using the training data, a Random Forest (RF) model was trained to classify each reference point into fire, non-fire disturbance, or no disturbance. Our RF approach showed that short-wave infrared (SWIR)-based moisture indices (e.g., B7, TCW) were strong predictors for detecting forest disturbances (Supplementary Fig. 11) likely because of their sensitivity to vegetation water content and canopy structure64. Finally, we applied the trained RF model to the full classification stack to consistently map the disturbance types from 1986 to 2020 across China’s forests, assuming that the spectral trajectories derived from reference data period 2001–2020 can be extrapolated to the whole mapping period 1986–2020. However, note that our approach was meant to detect relatively acuate and discrete disturbances that caused canopy opening, rather than subtle changes of forest structure or composition resulted from low intensive silvicultural practices and chronic disturbances.Year of disturbanceWe used the LandTrendr to determine the year of disturbance as the onset of magnitude of spectral change. Since we ran LandTrendr on five spectral indices, there were five possible years of disturbance for each pixel. Thus, we determined the year of disturbance using the median value from at least three different indices (i.e., NDVI, NBR, TCW, B5, B7). In this way, we only kept pixels that were detected as disturbances using at least three indices, thus reducing commission errors. The year with the greatest spectral changes generated by the LandTrendr often had an accuracy within 3 years11. A confidence level was also assigned to each disturbed pixel based on numbers of indices which showed possible disturbance events. Specially, low, medium, and high confidence were assigned if the disturbance was detected by three, four, or five spectral indices, respectively.ValidationsWe validated the disturbance map at the pixel and national levels. At the pixel level, we validated the final map using the validation sub-sample described in the previous section. We derived a confusion matrix to report user’s and producer’s accuracy (Supplementary Table 1) as the main accuracy assessment metrics. At the national level, we compared forest disturbance detected in this study to available existing dataset. Specifically, we compared the area of forest fire disturbance between our study and the national fire records during 2003–2009 (Supplementary Fig. 5). We compared the disturbance rates between our study and Landsat-derived global forest cover changes from 2001 to 20193 (Supplementary Fig. 4).Post-processingWe applied a series of spatial filters to minimize the unrealistic outliers from two potential sources of uncertainty, including speckle in time-series spectral trajectories or misregistration among images. This may lead to individual pixel or small patches including only a few pixels, which were (a) detected as disturbances, thus increasing the commission errors, or (b) not detected as disturbances, while their surrounding pixels were mostly disturbed, thereby increasing the omission errors. To address the issue (a), we removed all single-pixel disturbance patches through setting the minimum mapping unit as two 30 × 30 m2 pixels (0.18 ha). To address the issue (b), we applied a 3 by 3 moving window to fill holes through assigning the year of disturbance based on the years in the surrounding pixels. Finally, we smoothed the year of disturbance by assigning the center pixel using majority rules from surrounding pixels within the 3 by 3 windows, thus accounting for artefacts associated with uncertainties in the correct identification of the disturbance year.Characterizing disturbance regimes and their trendsWe characterized the disturbance regime using five indicators within each 0.5° grid cell (n = 1946) across China’s forests based on annual forest disturbance maps generated from the previous step. Within each grid cell, we calculated (1) total annually disturbed forest area (km2 yr−1), (2) percentage of forest disturbed annually (% yr−1), as annual disturbed forest area divided by the total forested area, (3) disturbance size (ha), as the number of disturbed pixels for each individual patch using an eight-neighbor rule, (4) disturbance frequency (# of patches per 1000 km2 forested area each year), as the number of disturbance patches per year divided by the total forested area, (5) disturbance severity (ΔNDVI = NDVIt−1 − NDVIt+1), as magnitude of NDVI change 1 year before and 1 year after disturbance, obtained from the LandTrendr analysis. We used (1) and (2) to characterize the disturbance rate, and (3)–(5) to describe the patch characteristics. The (2) and (4) were normalized by forest area within each grid cell, thus making them comparable among grid cells. For (3)–(5), we only calculated the patch size >0.45 ha (five 30 × 30-m2 pixels), because patches TC2000), and the expansion of forested area from 1986 to 2000 (e.g., TC1986 20% following Liu et al., (2019). We should note that our study area did not include the newly afforested area after 2000. All analyses were performed within the forest mask, thus excluding the potential confounding factors from other land cover types. The description of TC1986 and TC2000 can be found in3,32. More
