Improving LandTrendr Forest Disturbance Mapping in China Using Multi-Season Observations and Multispectral Indices

Abstract

Forest disturbance detection is of great significance for understanding forest dynamics. The Landsat-based detection of the Trends in Disturbance and Recovery (LandTrendr) algorithm is widely used for forest disturbance mapping. However, there are still two limitations in LandTrendr: first, it only used for summer-composited observations, which may delay the detection of forest disturbances that occurred in autumn and winter by one year, and second, it detected all disturbance types simultaneously using a single spectral index, which may reduce the mapping accuracy for certain forest disturbance types. Here, we modified LandTrendr (mLandTrendr) for forest disturbance mapping in China by using multi-season observations and multispectral indices. Validations using the randomly selected 1957 reference forest disturbance samples across China showed that the overall accuracy (F1 score) of forest disturbance detection in China was improved by 21% with these two modifications. The mLandTrendr can quickly and accurately detect forest disturbance and can be extended to national and global forest disturbance mapping for various forest types.

1. Introduction

Forest disturbance, such as fire, harvest, and insect infestation, is a main driving force in forest ecosystem dynamics [1]. It can even cause the transformation of forests from carbon sinks to carbon sources through carbon release by forest fire and the loss of carbon stocks caused by forest harvest and insects [2,3]. Information on the specific timing, location, and extent of forest disturbance is critical for forest management and decision support for achieving carbon neutrality [4]. This information requires accurate forest disturbance mapping based on consistent and continuous remote sensing data such as the Landsat time series [5,6].

Many algorithms have been developed for forest disturbance detection based on the Landsat time series [7,8,9,10]. The day-scale algorithms, such as CCDC [11], COLD [8], and EWMACD [12] can detect the specific date of forest disturbance using all available Landsat images but they require a large amount of data storage and computing resources. The year-scale algorithms, which detect the occurring year of forest disturbance using the yearly-composited Landsat observations, require less data storage and computing resources and are more suitable for large-area forest disturbance mapping [13], particularly in the case of local computing applications where researchers need to consider both storage and computational efficiency. The Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) [9] is one of the most widely-used year-scale algorithms in forest disturbance mapping [14,15,16,17] due to its efficient operation and suitability for parallel computing on the Google Earth Engine (GEE) cloud platform [18,19]. However, there are still two limitations in LandTrendr which may hinder its applicability and introduce uncertainties in downstream forest carbon modeling [20].

First, the LandTrendr only used yearly summer-composited (usually from June to September) observations to monitor forest disturbance [21,22,23] which could delay detected forest disturbances that occurred in autumn and winter by one year. Figure 1 shows an example of delayed forest disturbance detection. A forest fire occurred between 19 October 2009 and 4 November 2009 but the occurring year was detected as 2010 by LandTrendr; while using the winter-composited observations, the occurring year of the forest fire could be correctly detected as 2009. This difference would cause an overestimation of forest carbon consequence before the disturbance [24]. Therefore, it required the integration of observations in other seasons to improve LandTrendr for forest disturbance mapping.

Second, the LandTrendr used a single input spectral index such as NBR, NDVI, TCW, or TCA [25,26,27,28] to detect all forest disturbance types simultaneously which may reduce the accuracy of detecting certain forest disturbance types. Figure 2 shows an example of using a single spectral index of NBR in forest disturbance detection. A forest harvest occurred between 27 July 2000 and 13 September 2000, but it was omitted by LandTrendr when only using NBR as the detector; while using TCW, the occurring year of forest harvest could be correctly detected as 2000. The importance of using multiple spectral indices has been demonstrated in many studies [1,11,29,30,31] and the disturbance detection from different input indices has been merged together to obtain a better result of the detected forest disturbance [7,8,32]. However, it was the ensemble of the detected results using different indices, not the improvement of the LandTrendr algorithm. Therefore, incorporating multiple spectral indices was required to improve LandTrendr for forest disturbance mapping.

In this study, we aimed to modify the LandTrendr algorithm (hereinafter referred to as mLandTrendr) by using multi-season observations and multispectral indices for forest disturbance mapping. There were three objectives: (1) to incorporate multi-seasonal data into the annual composition instead of only using summer observations; (2) to incorporate multispectral indices into the forest disturbance mapping instead of only using a single spectral index; and (3) to evaluate the performance of the mLandTrendr for forest disturbance mapping in China not only at the pixel level but also at three Military Grid Reference System (MGRS) tiles. The modified algorithm can be useful for mapping forest disturbance at the national and global scales.

2. Data and Methods

2.1. Study Area

A total of 3 MGRS tiles and 3914 reference samples across China were selected as the study area (Figure 3). The three MGRS tiles, including Daxinganling, Heilongjiang Province (tile ID: T52UCU), Shennongjia, Hubei Province (tile ID: T49RDQ), and Sanming, Fujian Province (tile ID: T50RNP) were selected for their different locations, forest cover types, and forest disturbance types. T52UCU was selected for its deciduous forest in North China with forest disturbance types of fire and harvest. T49RDQ was selected for its central location and rich forest cover types and it mainly includes the forest disturbance types of harvest and fire. T50RNP was selected for its evergreen forests in South China with forest disturbance types of harvest, insects, and typhoons. In view of the wide similarity of natural conditions and anthropogenic pressures in the study area, the forest disturbance types in the study are representative across China, meeting the need to verify the applicability of the algorithm nationwide.

2.2. Data

2.2.1. Reference Disturbance Samples

A total of 3914 randomly selected reference forest disturbance samples from 1986 to 2021 in China (Figure 3) based on a 30-m spatial resolution forest map were used to calibrate and validate the forest disturbance detection algorithm. Each reference forest disturbance sample was interpreted with the help of all available Landsat 5/7/8 time series images, Sentinel-2 images, PlanetScope images (https://www.planet.com/explorer/), and Google Earth high-resolution images to obtain the accurate occurring date of forest disturbance. Figure 4 shows an example of the interpretation of each forest disturbance event confirmed by at least two clear images before and after the disturbance. There were 645 reference samples with at least one disturbance event between 1986 and 2021 and 21% of these samples had more than one disturbance event. The disturbance types in the reference sample were classified as harvest, fire, insects, and others.

2.2.2. Landsat Time Series Data

All available Landsat 5/7/8 TM, ETM +, and OLI Collection 2 Tier 1 surface reflectance from 1986 to 2021 were used to detect the forest disturbance. The long record of Landsat data and 30-m spatial resolution [33,34] made it sufficient to monitor forest disturbances [35,36]. The Time-series-based Reflectance Adjustment (TRA) algorithm [37] was used to harmonize the reflectance difference between data from different Landsat sensors. Clouds and cloud shadows were screened using the quality control (QA) values [38] and the remaining clouds and cloud shadows were screened using the time series filter proposed by Shang et al. [7]. Except for the summer observations, observations in other seasons were also needed for improving forest disturbance detection. A set of commonly-used spectral indices such as Normalized Burn Ratio (NBR) [39], Normalized Difference Moisture Index (NDMI) [40], Normalized Difference Vegetation Index (NDVI) [41], Tasseled Cap Wetness (TCW) [42], and Tasseled Cap Angle (TCA) [43] (

Table 1. The candidate spectral indices that would be selected and evaluated for improving forest disturbance detection in China.

Index	Equation	Reference
NBR	(NIR − SWIR2)/(NIR + SWIR2)	[39]
NDMI	(NIR − SWIR1)/(NIR + SWIR1)	[40]
NDVI	(NIR − Red)/(NIR + Red)	[41]
TCW	0.0315 × Blue + 0.2021 × Green + 0.3102 × Red + 0.1594 × NIR − 0.6806 × SWIR1 − 0.6109 × SWIR2	[42]
TCA	Arctan[(0.2043 × Blue + 0.4158 × Green + 0.5524 × Red + 0.5741 × NIR − 0.3124 × SWIR1 + 0.2303 × SWIR2)/(−0.1603 × Blue − 0.2819 × Green − 0.4934 × Red + 0.7940 × NIR − 0.0002 × SWIR1 − 0.1446 × SWIR2)]	[43]

) were calculated as the candidate spectral indices to modify Landtrendr for forest disturbances detection in China.

2.3. Methods

We modified LandTrendr (mLandTrendr) for forest disturbance mapping by using multi-season observations and multispectral indices. A total set of 3914 reference forest disturbance samples were used to calibrate the optimal threshold for each parameter in mLandTrendr and to validate its performance of forest disturbance detection. Figure 5 shows the flowchart of the algorithm modifications.

The first modification was to incorporate multi-season observations into the yearly composition instead of only using the summer observations. The mLandTrendr would monitor the summer composition and an additional composition twice before confirming a forest disturbance. Sensitivity analysis showed that the optimal compositing period was from 1 October to 31 December (

Table A1. Sensitivity analysis using different time periods for the yearly composition.

Compositing Periods	Commission Error (%)	Omission Error (%)	F1 Score (%)
03/21–06/23	55.59	64.66	39.36
06/20–08/31	30.96	64.40	46.98
06/21–09/23	34.26	62.83	47.49
10/01–12/31	40.93	59.95	47.74
12/21–03/23	48.04	61.78	44.04

). It should be noted that a direct combination of the two detections would cause a commission by the summer detection. Therefore, a confirming procedure was proposed to remove the duplicates if the occurring year of a forest disturbance detected by the summer detection is one year later than that by the winter detection. As there were two detections, the output of mLandTrendr can accurately pinpoint the timing of detected disturbances to specific seasons, such as summer and winter.

The second modification was to incorporate multispectral indices into the forest disturbance detection instead of only using a single index. All possible combinations with 1–5 indices from the five candidate spectral indices (Table 1) were evaluated, and the combination of NBR, TCW, NDMI, and TCA was determined as the optimal combination by the sensitivity analysis (

Table A2. Sensitivity analysis using different combinations of spectral indices. The number from 1 to 5 in the name of the combination indicates using the index of NBR, NDMI, TCW, NDVI, and TCA. “Tc” is the optimal change threshold used in disturbance confirmation.

Index Combinations	Commission Error (%)	Omission Error (%)	F1 Score (%)	Tc
C12345	34.51	31.69	66.87	0.999
C1234	32.43	30.77	68.39	0.999
C1235	30.58	30.15	69.63	0.999
C1245	35.57	38.15	63.11	0.999
C1345	34.13	32.31	66.77	0.999
C2345	34.45	33.84	65.85	0.999
C123	31.6	31.38	68.51	0.999
C124	33.77	39.07	63.46	0.999
C125	31.97	38.46	64.62	0.999
C134	36.06	34.55	64.68	0.99
C135	33.42	35.86	65.33	0.99
C145	40.1	39.79	60.05	0.95
C234	32.17	33.84	66.98	0.999
C235	30.19	33.85	67.93	0.999
C245	39.71	35.6	62.28	0.99
C345	41.82	41.36	58.41	0.99
C12	33.04	39.68	63.47	0.99
C13	29.85	38.87	65.33	0.99
C14	29.05	43.7	62.78	0.99
C15	24.18	44.5	64.09	0.99
C23	32.24	39.14	64.12	0.99
C24	42.28	42.9	57.41	0.95
C25	27.74	46.92	61.21	0.99
C34	37.89	46.38	57.55	0.95
C35	34.53	46.11	59.12	0.95
C45	37.12	61.39	47.84	0.9
C1	43.04	39.68	58.59	0.95
C2	40.22	42.63	58.55	0.9
C3	32.75	48.26	58.48	0.9
C4	30.37	64.34	47.16	0.9
C5	29.67	60.59	50.52	0.9

). In mLandTrendr, a normalized change matrix (Equation (1)) derived from the above spectral indices was calculated as the indicator of forest disturbance detection, and change in any one of the spectral indices would result in a large value in the normalized change vector [8]. The chi-squared distribution with a change threshold (Tc) was used to measure the normalized change matrix (Equation (2)) and a forest disturbance would be confirmed if there were several (Tn) consecutive normalized change matrices exceeding the change magnitude calculated from the chi-squared distribution. Indeed, 1 was determined as the optimal threshold of Tn by the sensitivity analysis (

Table A3. Sensitivity analysis using different thresholds for the parameter Tn.

Tn	Commission Error (%)	Omission Error (%)	F1 Score (%)
1	30.58	30.15	69.63
2	18.87	44.84	65.67
3	13.30	51.85	61.91

) and 0.999 was determined as the optimal threshold of Tc by the sensitivity analysis using different thresholds of 0.9, 0.95, 0.99, 0.999, and 0.9999 (Figure 6). Other parameters in the mLandTrendr were also calibrated and

Table 2. The optimal thresholds for the parameters in the mLandTrendr.

Parameter	Meaning	Optimal Threshold	Sensitivity Analysis
Compositing periods	One observation value is selected to represent the annual spectral value	1 October to 31 December	Table A1
Indices	Spectral indices combinations are used to monitor forest disturbance	NBR, NDMI, TCW, TCA	Table A2
Tn	Several consecutive normalized change matrices	1	Table A3
Tc	Change threshold used in the chi-squared distribution to measure the normalized change matrix	0.999	Figure 6
maxSegments	Maximum number of segments to be fitted on the time series	6	Table A4
pvalThreshold	If the p-value of the fitted model exceeds this threshold, then the current model is discarded and another one is fitted using the Levenberg–Marquardt optimizer	0.01	Table A5
spikeThreshold	Threshold for dampening the spikes (1.0 means no dampening)	0.85	Table A6
recoveryThreshold	If a segment has a recovery rate faster than 1/recoveryThreshold (in years), then the segment is disallowed	0.5	Table A7
bestModelProportion	Takes the model with the most vertices that have a p-value that is at most this proportion away from the model with the lowest p-value	0.5	Table A8

shows the summary of the optimal values for these parameters.

$$CM\left(i\right)={\sum }_{j=1}^k{\left(\frac{p_j\left(i\right)-{\widehat{p}}_j\left(i\right)}{\mathrm{m}\mathrm{a}\mathrm{x}(rmse,mini\_rmse)}\right)}^2$$

(1)

$$\mathrm{m}\mathrm{i}\mathrm{n}(CM\left(i\right),CM\left(i+1\right),\dots ,CM(i+Tn-1\left)\right)>X^2\left(Tc,k\right)$$

(2)

where i is the i-th observation in the yearly composited time series, j is the j-th spectral index, k is the number of indices used for forest disturbance detection, $p_j\left(i\right)$ is the observation for the i-th spectral index, ${\widehat{p}}_j\left(i\right)$ is the prediction for the ith spectral index by the linear regression, CM(i) is the i-th normalized change matrix, Tc is the change threshold used in the chi-squared distribution, Tn is the number of consecutive normalized change matrices exceeding the change magnitude derived from chi-squared distribution, $rmse$ is the root mean square error between observations and predictions, and $mini\_rmse$ is a minimum RMSE used to prevent the model from overfitting.

2.4. Evaluation Methods

The commission error (Equation (3)), omission error (Equation (4)), and F1 score (Equation (5)) were used to indicate the accuracy of forest disturbance detection. Commission refers to the forest disturbances that are only detected by the LandTrendr algorithm but are not detected by forest disturbances in the reference samples. Omission refers to the forest disturbances that existed in the reference samples but are not monitored by the LandTrendr algorithm. The F1 score is an indicator to quantify whether commission and omission are balanced [8]. If the detected forest disturbance is consistent with the reference sample in both time and space, the correctness of a disturbance event can be determined. Sensitivity analysis was conducted using different thresholds to calibrate the optimal threshold for each parameter in the algorithm. In total, 50% of reference samples (shortened as calibrating samples) were randomly selected to calibrate the parameters and the remaining 50% (shortened as validating samples) was used to evaluate the algorithm accuracy.

Except for the point evaluation using the reference forest disturbance samples, three MGRS tiles were selected to evaluate the performance of wall-to-wall forest disturbance detection by LandTrendr and mLandTrendr. The cloud-free RGB images before and after the forest disturbance were used to indicate whether the entire forest disturbance events were captured correctly.

$$commission=\frac{a}{b}\times 100\%$$

(3)

$$omission=\frac{c}{d}\times 100\%$$

(4)

$$F1 score=\frac{\left(1-commission\right)\times (1-omission)}{(2-commission-omission)}\times 200\%$$

(5)

where a is the number of detected forest disturbances that disagree with the reference; b is the total number of detected forest disturbances; c is the number of reference forest disturbances that disagree with the detection; and d is the total number of forest disturbances in the reference.

3. Results

3.1. Forest Disturbance Detection at the Reference Samples

The LandTrendr is modified (mLandTrendr) by integrating multi-season observations and multispectral indices. The accuracy of the mLandTrendr using the optimal thresholds was evaluated based on the validating samples and it achieved an F1 score of 70% with a commission error of 29% and an omission error of 32%.

Table 3. The confuse matrix calculated from the mLandTrendr using the optimal thresholds based on the 1957 validating samples.

Detection	Reference		Accuracy
	No Disturbance	Disturbance
No disturbance	1516	122	92.6%
Disturbance	105	260	71.2%
Accuracy	93.5%	68.1%	88.7%

shows the calculated confusion matrix using the validating samples and it shows that the overall accuracy of forest disturbance detection by mLandTrendr was 89%.

To indicate the improvement, the original LandTrendr using its default parameter settings was also evaluated based on the same validating samples and it got an F1 score of 47%. The default parameter settings might not be optimal for forest disturbance detection in China, and for a fair comparison, each parameter in the original LandTrendr was calibrated using the calibrating samples. Different input spectral indices of NBR, NDMI, TCW, or TCA were compared and the yearly composition of the winter period was also added as a comparison. The highest F1 score of 49% was achieved by using NBR with the summer-composited observations (

Table A9. Accuracy comparison using different single input spectral indices with the summer-composited observations in the original LandTrendr algorithm.

Index	Commission Error (%)	Omission Error (%)	F1 Score (%)
NBR	51.12	51.57	48.62
NDMI	60.70	57.64	40.77
TCW	55.81	53.08	45.51
TCA	61.29	62.30	38.20

), while for using the winter-composited observations, NBR also achieved the highest F1 score of 47% among the four spectral indices (

Table A10. Accuracy comparison using different single input spectral indices with the winter-composited observations in the original LandTrendr algorithm.

Index	Commission Error (%)	Omission Error (%)	F1 Score (%)
NBR	53.21	53.52	46.63
NDMI	41.30	63.81	44.78
TCW	41.67	66.22	42.78
TCA	55.66	63.09	40.29

).

These results suggested a large accuracy improvement of 21% in forest disturbance detection across China by the mLandTrendr, demonstrating the necessity of integrating multi-season observations and multispectral indices into forest disturbance detection across China.

3.2. Forest Disturbance Detection at the MGRS Tiles

The applicability of the mLandTrendr forest disturbance detection in China was also evaluated at the three MGRS tiles. Figure 7, Figure 8 and Figure 9 show the comparisons of forest disturbance detection by the original LandTrendr and mLandTrendr at the MGRS tiles of T50RNP (Figure 7), T49RDQ (Figure 8), and T52UCU (Figure 9). The RGB composited images were used as the reference to indicate the performance of forest disturbance detection. For the original LandTrendr, NBR was used as the single input spectral index and the summer months were used as the yearly composition period according to their highest performance across China (Table A9). Compared to the original LandTrendr, the forest disturbance detected by the mLandTrendr was significantly improved, not only in the area of forest disturbance but also in the integrity of each forest disturbance event. The spatial pattern of the forest disturbance monitored by the mLandTrendr algorithm was close to the disturbance area in the Landsat-8 RGB-composited (R: SWIR1, G: NIR, B: R) images and only a little difference could be observed in some pixels.

4. Discussions

Forest disturbance has significant effects on the forest carbon budget and ecosystem services [44]. Detecting forest disturbance is essential for sustainable forest management and climate change mitigation [45,46]. LandTrendr is one of the most widely-used algorithms in forest disturbance detection [1,29,31,47] but it has two limitations: first, only used summer-composited observations, which may delay the detection of winter forest disturbances by one year, and second, detecting all forest disturbance types simultaneously using a single variable, which may reduce the detection accuracy for certain forest disturbance types. This study focused on improving LandTrendr for forest disturbance detection in the above two aspects using multi-season observations and multispectral indices. With the improvement, the F1 score of forest disturbance detection in China was improved from 49% to 70%.

The use of multispectral indices made it possible to monitor all forest disturbance types simultaneously [7,32,48]. This study selected five commonly used spectral indices of NBR, NDMI, TCW, NDVI, and TCA as the potential detectors and tested all possible combinations of multispectral indices from 1 to 5. The results showed that the combination of NBR, TCW, NDMI, and TCA performed the best (Table A2). It should be noted that the integration of multispectral indices was not simply merging the detection results of each index, but a direct modification of the algorithm on GEE. Similarly to the Continuous Monitoring of Land Disturbance (COLD) [8] and Near-real-time Monitoring of Land Disturbance (NRT-MONITOR) [7], we used a chi-square distribution to integrate the indices and then confirmed forest disturbance with a change probability of 0.999. This integration can make it possible to monitor all forest disturbance types because the change caused by forest disturbance in one index could be expanded in the final change matrix calculated from the chi-square distribution [7,8].

Using multi-season observations could avoid the delayed detection of forest disturbance that occurred in autumn and winter. The mLandTrendr model would monitor the summer composition and the winter composition twice before confirming a forest disturbance. Sensitivity analysis showed that the optimal winter compositing period was 1 October to 31 December (Table A1). It should be noted that the simple combination of the two detections would cause duplicates for a summer or winter disturbance event in the occurring year. To avoid this kind of commission error, we proposed a confirming procedure to remove the duplicates if the occurring year of a forest disturbance found by summer detection was to be one year later than that by winter detection.

The detected forest disturbance by mLandTrendr was compared to a forest loss and gain product (shortened as Hansen’s product) at the 30-m spatial resolution [49]. In Hansen’s product, forest loss was defined as a stand-replacement disturbance [49] and, to keep consistency, the forest loss disturbance from the year 2001–2021 was compared. Based on the same validating samples, Hansen’s product achieved an F1 score of 53% with a commission error of 12% and an omission error of 62%, which was 17% lower than the F1 score of mLandTrendr.

The time efficiencies of LandTrendr and mLandTrendr were compared to CCDC, a year-scale algorithm [11] implemented in GEE. Running 1000 pixels with these 3 algorithms took 10.01 s, 55.45 s, and 888.87 s, respectively, suggesting that the year-scale algorithms are more time efficient than the day-scale algorithms. This also explained why we focused on improving the year-scale algorithm rather than the day-scale algorithm in this study. There were also some limitations in the original LandTrendr and mLandTrendr. First, the calibrated optimal thresholds in China may not be optimal for other countries and new calibrations are needed if using the mLandTrendr to monitor forest disturbance in other areas [50]. Second, the mLandTrendr can only determine the occurring year of the monitored forest disturbance [9], not the exact date with a month and a day [51,52]. Third, just like the original LandTrendr implemented on GEE, mLandTrendr also required a minimum size of 11 × 11 pixels (330m × 330m) for detecting a forest disturbance event, and forest disturbance events smaller than this size would be screened and excluded. This size could be reduced to one Landsat pixel (30 m) if a smaller disturbance event was also needed. Fourth, mLandTrendr was also a pixel-level algorithm and the integration of spatial information from its nearby pixels was suggested for future modifications to achieve higher accuracy of forest disturbance detection. Last, for boreal forests, if the snow cover lasts for the entire period from 1 October to 31 December, the forest disturbance that occurred in winter could be hardly detected by both LandTrendr and mLandTrendr. However, even in the north of China, only a small proportion of boreal forests had a long winter starting from 1 October and mLandTrendr was applicable if there had been one available snow-free and cloud-free observation during the composition period.

5. Conclusions

In this study, we modified LandTrendr for forest disturbance mapping by using multi-season observations and multispectral indices. The modified algorithm was calibrated and validated by a set of 3914 reference forest disturbance samples in China. It used multi-season remote sensing data rather than just summer data and incorporated multiple vegetation indices rather than just one index. With the two modifications, the accuracy of forest disturbance monitoring across China was improved by 21%. The modified algorithm can quickly and accurately detect forest disturbance and can be extended to global forest disturbance mapping with various forest types.