Topographic effects amplify forest disturbances detected by yearly wide-time-window Landsat time series

ABSTRACT

Widespread topographic effects in remote sensing images of mountainous regions with rugged terrain severely hinder satellite-based global land surface monitoring and change detection. However, topographic effects in Landsat time series (LTSs) disturbance detection remain highly debated and unquantified. In this study, we proposed a novel postprocessing approach to quantify the impact of topographic effects on long time-series (1989 − 2020) forest disturbance detection, taking the Hengduan Mountains Region (HDMR) of southwest China as an example. This approach applied a pixel-by-pixel simulation based on a semiempirical sun-canopy-sensor with a C-correction (SCS+C) model to identify and remove topography-induced disturbances from the LandTrendr-detected forest disturbances. Filtering of the detected forest disturbances with different patch sizes was conducted, and its effectiveness in removing topography-induced disturbances was quantitatively evaluated. The results showed that 10.43% of the total area of LandTrendr-detected forest disturbances was topography-induced. Removing topography-induced disturbances increased the traditional and area-adjusted overall accuracies (OAs) of the forest disturbance map by 3.50% and 0.62%. Topography-induced disturbances occurred mainly on northwest-facing (300°−360°) slopes with a gradient of 30°−55° and higher latitude and between the day of year (DOY) pairs (such as DOYs 90, 330, and 360) with a considerable disparity in solar azimuth and zenith angle. The pixel-by-pixel simulation approach outperformed spatial noise filtering with a minimum mapping unit (MMU), which considerably decreased the area percentage of topography-induced disturbances at the expense of decreasing the producer’s accuracy (PA) in detecting disturbed forests. This study highlighted the importance of removing topographic effects from the detected forest disturbances and shed new light on the promising potential of the postprocessing topographic correction method to replace the preprocessing topographic correction of the LTSs to map forest disturbances in mountainous regions globally.

KEYWORDS:

• Topographic effects
• Landsat time series change detection
• Forest disturbance
• SCS+C topographic correction model

1. Introduction

Mountain regions cover only approximately 25% of the Earth’s land surface (Meybeck, Green, and Vorosmarty 2001) but harbor more than 85% of the terrestrial global biodiversity (Rahbek et al. 2019). However, mountainous forest disturbances that cause permanent or temporary loss of forest threaten the habitat and resource availability (Newbold et al. 2015), species richness (Daskalova et al. 2020; Duraes et al. 2013), and long-term persistence of species (Blandon et al. 2016). Recently, mountainous forests have become hotspots and key areas for global land surface change detection (Adamczyk and Osberger 2015). With the explosive growth of open-access earth observation data, forest disturbances have been increasingly detected at finer spatial and temporal scales from a long-term perspective (Hansen and Loveland 2012; Woodcock et al. 2020; Zhu et al. 2020). Their mapping accuracies relied on the degree of consistency of images (Qiu et al. 2019). However, widespread topographic effects in images caused by terrain orientation and solar geometry (Fan, Koukal, and Weisberg 2014; Van den Hoek et al. 2021) often ruin the consistency of images and potentially affect the accuracy of forest disturbance detection in mountainous areas (Vanonckelen, Lhermitte, and Van Rompaey 2015).

Numerous topographic correction models have been proposed and applied to eliminate topographic effects. These models can be grouped into three broad categories: empirical models (Ekstrand 1996; Gao and Zhang 2009), physical models (Gu and Gillespie 1998; Yin et al. 2018), and semiempirical models (Teillet, Guindon, and Goodenough 1982). Each topographic correction model has a suitable application (Lu et al. 2008). A semiempirical sun-canopy-sensor with C-correction (SCS+C) model preserves the vertical growth of trees (Soenen, Peddle, and Coburn 2005) and blends the merits of empirical and physical models. It has been demonstrated to perform better than other topographic correction models in mountainous areas of southwest China (Wu, Jin, and Fan 2016). Hence, the SCS+C model is especially suitable for mountainous forested areas.

Since Landsat images are freely available (Woodcock et al. 2008), numerous change detection algorithms have been developed, reviewed and widely applied over local to global scales (DeVries et al. 2015; Galiatsatos et al. 2020; Huang et al. 2010; Kennedy, Yang, and Cohen 2010; Zhu and Woodcock 2014). Some of them used yearly Landsat time series (LTSs) composited by a pixel-based image compositing (PBIC) approach (Huang et al. 2010; Kennedy, Yang, and Cohen 2010; Ochtyra, Marcinkowska-Ochtyra, and Raczko 2020), and the other employed LTSs included all available observations (Bullock, Woodcock, and Holden 2020; DeVries et al. 2015; Zhu and Woodcock 2014). However, topographically corrected Landsat image products remain unavailable (Chance et al. 2016; Yin et al. 2022). Whether topographic correction should be considered for LTSs-based change detection is user-dependent (Adhikari et al. 2016; Chance et al. 2016; Vanonckelen, Lhermitte, and Van Rompaey 2015). Several studies indicated that mapping accuracies of forest change using topographically corrected LTSs were improved by 1.15%−10% (Shimizu et al. 2016; Tan et al. 2013; Van den Hoek et al. 2021). Conversely, topographic correction exerted very limited or even negative impacts on LTSs-based forest disturbance detection (Adhikari et al. 2016; Chance et al. 2016; Qiu et al. 2019). Contradictory conclusions were drawn on the necessity and effectiveness of topographic correction when conducting yearly LTSs-based change detection. Unfortunately, little quantitative simulation work of topographic effects has been reported until now.

In regions of rain-fed forests, the growing seasons of forests maintain a step with rainy seasons characterized by frequent and heavy cloud cover; therefore, very limited cloud-free Landsat images are available so that intra-year change detection algorithms are unsuitable in these regions (Pasquarella et al. 2022; Zhang et al. 2021; Zhu 2017). Yearly LTSs composed of wide-time-window images with diverse illumination conditions have been used for forest disturbance detection (Grogan et al. 2015; Shimizu et al. 2019; Yang et al. 2019), which could cause erroneous identification of forest disturbance (topography-induced disturbance) (Banskota et al. 2014). However, to what extent the interannual variations at yearly wide-time-window LTSs related to topographic effects affect forest disturbance detection remains to be quantified. Furthermore, some studies have argued that filtering small-sized patches can effectively eliminate topography-related spatial noise in the results of pixel-based forest disturbance detection (Chance et al. 2016; Kennedy et al. 2015). However, the performance of the filtering postprocessing operation in removing topography-induced disturbances was also unclear.

Therefore, taking the Hengduan Mountains Region (HDMR) in southwest China as the study area, we combined an SCS+C model with a LandTrendr algorithm (Kennedy, Yang, and Cohen 2010) to quantify topographic effects in forest disturbance detection based on yearly wide-time-window LTSs. The main objectives of this study are to (1) quantify topography-induced disturbances when mapping forest disturbances using yearly wide-time-window LTSs through a pixel-by-pixel simulation, (2) calculate the area contribution of topography-induced disturbances to LandTrendr-detected forest disturbances, and (3) compare performances between filtering postprocessing and removing the simulated topography-induced disturbances.

2. Data and methodology

2.1 Study area

The HDMR, situated in southwestern China (Figure 1), is an important geographical unit of the Tibetan Plateau (Li 1987), featuring a series of north-south-oriented mountain ranges (i.e. Gaoligong Mountains) and rivers (i.e. Nu River and Jinsha River). This region exhibits extremely rugged terrain, with elevations ranging from approximately 850 m in the Nu River valley to 7556 m at the peak of Gongga Mountain. Areas with a slope above 25° cover approximately 53.03% of the whole HDMR, and forested areas steeper than 30° occupied approximately half of the initial forested area in the 1990s (Figure 1a). Major vegetation types in the HDMR include succulent thorny shrubs, montane mixed needle- and broad-leaved forests, montane dark coniferous forests, alpine shrub meadows, and alpine sparse vegetation. The growing season of vegetation in the HDMR is in step with the rainy (May to October) seasons (He et al. 2005; Li et al. 2010). The dry season presents little precipitation but is foggy, whereas the rainy season is rainy and foggy, in which 60% to 90% of the annual precipitation falls (Li et al. 2010).

Figure 1. Geographical location and topography of the study area. The insets (a) and (b) show the area percentage of forested regions with different slope levels in the HDMR and footprints of the 36 Landsat scenes covering the study area, respectively. (c) high-resolution Google Earth images and the LandTrendr-derived NBR pixel-based trajectories, which are used for visual validation of the detected forest disturbances, at three representative sites (A, B and C). DEM data used here are sourced from the Shuttle Radar Topography Mission V3 (SRTM−3) in GEE.

2.2 Data and preprocessing

The major steps of this study include (Figure 2): (1) simulating the normalized burn ratio (NBR) on different days of the year (DOYs 30, 60, 90, 300, 330, and 360) that correspond to acquisition months of composite LTSs in forested areas across the HDMR using the SCS+C model to determine topography-induced disturbances; (2) detecting forest disturbances using LandTrendr based on yearly LTSs and then excluding topography-induced disturbances from the resultant forest disturbances; (3) filtering the LandTrendr-detected forest disturbances with different patch sizes; (4) evaluating the accuracy and analyzing spatiotemporal characteristics of the topography-induced disturbances; and (5) comparing the areas and accuracies of forest disturbances when excluding simulated topography-induced disturbances and filtering forest disturbances with different patch sizes.

Figure 2. A flowchart of detecting and filtering forest disturbances and simulating topography-induced disturbances. α, β, θ, Ω and i denote slope, aspect, solar zenith angle, solar azimuth angle and solar incidence angle, respectively.

2.2.1 Landsat time series images and preprocessing

Landsat Collection−1 Surface Reflectance Tier 1 products consisting of images from Landsat sensors (TM, ETM+ and OLI) that were atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al. 2006) were adopted in this study. The normalization method was used to reduce the differences in reflectance among sensors (Roy et al. 2016). It was impossible to use clear observations of Landsat images acquired in growing seasons to composite yearly cloud-free LTSs in the HDMR due to the heavy cloud cover. Therefore, we selected non-growing season (Jan.-Mar. and Oct.-Dec.) images within the 36 Landsat footprints (Figure 1b) from 1989 to 2020 to composite yearly LTSs on the Google Earth Engine (GEE) cloud platform (Gorelick et al. 2017; Kumar and Mutanga 2018). Clouds and cloud shadows in images were removed to obtain all clear observations using the CFmask algorithm (Foga et al. 2017). Subsequently, the resulting clear observations were used to composite yearly cloud-free continuous LTSs through a PBIC method, which calculates the best available pixel value from all available cloud-free pixels for the location of each pixel over the time window of interest following the compositing criteria (Roy et al. 2010). In this case, we adopted the pixel-based median of clear observations per year as the best available pixel value.

A total of 216 Landsat images with cloud cover below 20% near DOYs 30, 60, 90, 300, and 360 were selected to produce monthly cloud-free composite images in January, February, March, October, November, and December, respectively, with a pixel-by-pixel median algorithm. Subsequently, these composite images were further used to calculate the monthly values of the C parameter of the SCS+C model.

2.2.2 Initial forest map and topographic data

The forest baseline map is required to exclude disturbances occurring in non-forested areas (DeVries et al. 2015). China Land-Use/Cover datasets (CLUDs) with a 1 km resolution in the late 1980s (Liu et al. 2003) were adopted to extract the raw forest baseline map, which was then resampled to a 30 m resolution. Moreover, a composite image of the maximum value of the normalized difference vegetation index (NDVI) in 1990 was produced from Landsat thematic mapper (TM) images. Removing the pixels with negative NDVI values from the raw forest baseline map yielded the final forest baseline map (Figure 2).

A 30-m-resolution DEM from the Shuttle Radar Topography Mission V3 (SRTM−3) provided by NASA’s Jet Propulsion Laboratory was used to extract topographic measures, including elevation, slope and aspect, in the GEE platform.

2.2.3 Reference data

Two groups of reference samples (Figure 1b) were separately collected in this study. One group consists of 566 blocks with more than 10 samples each year, which were generated in the overlapping areas of the topography-induced disturbances simulated by the SCS+C model and LandTrendr-detected forest disturbances using a stratified random sampling method. These blocks were used to compute the simulating accuracy of topography-induced disturbances. The other group comprised 804 disturbed forest sampling patches and 1087 stable forest sampling patches with patch sizes ranging from 3 to 1154 pixels. Different from the study of Li et al. (2021), where the size and allocation of samples were designed according to the existing experts’ experience, we determined the sample size and sample allocation according to Olofsson et al. (2014) to balance among user’s, producer’s and overall accuracies in order to produce a more uniform distribution of samples in space and time (Figure 1b). For each year of the study period, at least 10 disturbed sampling patches were separately collected to calculate the mapping accuracies of forest disturbances. To obtain these sampling patches, we conducted three field campaigns in August of 2019 and September of 2020 (Figure 1). During the field campaigns, occurrences of forest disturbances were recorded with a digital camera with a global positioning system (GPS) receiver device and subsequently the years of forest disturbances were obtained by interview investigations. Combining with visual interpretations of high spatial resolution images from Google Earth^TM and/or the LandTrendr-derived NBR pixel-based trajectories, disturbed forest samples were generated (Figure 1c). Specifically, patches featuring considerably decreasing spectral segments were identified as disturbed forest, and the year of occurrence was recorded, whereas patches with constant spectral trajectories were labeled as stable forest. For the years with fewer than 10 disturbed samples, a purposive sampling method along with visual interpretation was employed to collect supplementary samples.

2.3 Detecting forest disturbances

The LandTrendr method, a widely used algorithm based on the segmentation analysis of temporal trajectories (Kennedy, Yang, and Cohen 2010), on the GEE, was employed to identify forest disturbances (Kennedy et al. 2018). This method can find vertices and fit the temporal trajectory of each pixel into a series of straight-line segments with varying lengths by the point-to-point or regression method. From these resulting segments, information on the stability, disturbance, and regrowth of forests can be extracted based on the properties of individual segments (such as the magnitude and duration of change) (Kennedy, Yang, and Cohen 2010). Due to its high sensitivity to disturbance events (Grogan et al. 2015; Kennedy, Yang, and Cohen 2010) and better detection performance compared to other spectral indices (Cohen et al. 2018; Li et al. 2022), NBR (Eq. 1) was chosen as an indicator for detecting forest disturbance in the LandTrendr algorithm. As suggested by Grogan et al. (2015) and Li et al. (2021), the segments representing forest disturbance were extracted by a slide threshold value (0.277) of NBR change magnitude. Consequently, the segments with a change magnitude above 0.277 and duration below 3 years were labeled forest disturbances.

(1) $NBR=\frac{{\rho }_{nir}-{\rho }_{swir2}}{{\rho }_{nir}+{\rho }_{swir2}}$(1)

where ρ_nir and ρ_swir2 denotes near-infrared (NIR) and shortwave infrared 2 (SWIR2) reflectance, respectively.

2.4 Simulating topography-induced disturbances and filtering the LandTrendr-detected disturbances with different patch sizes

The process of simulating topography-induced disturbances was implemented in GEE (Brovelli, Sun, and Yordanov 2020; Pratico et al. 2021). The SCS+C topographic correction model was used to remove topographic effects in the surface reflectance of NIR and SWIR2 wavelengths when calculating the NBR. The formulation of the SCS+C model is given as follows.

(2) ${\rho }_n=\rho \frac{cos\alpha cos\theta +\mathrm{C}}{cosi+\mathrm{C}}$(2)

(3) $cosi=cos\theta cos\alpha +sin\theta sin\alpha cos(\mathrm{\Omega }-\beta )$(3)

(4) $cos\mathrm{\Omega }=\frac{sinhsinlatsin\delta }{\mathrm{c}\mathrm{o}\mathrm{s}hcoslat}$(4)

(5) $sinh=sinlatsin\delta +coslatcos\delta cost$(5)

(6) $\theta =90{ }^{\circ }-h$(6)

(7) $\rho =a+bcosi$(7)

(8) $\mathrm{C}=\frac{a}{b}$(8)

where ρ_n is the normalized reflectance, ρ is the uncorrected reflectance, α is the slope, θ is the solar zenith angle, i is the solar incidence angle calculated by Eq. 3, C represents a semiempirical parameter, Ω represents the solar azimuth angle, and β is the aspect, h is the solar elevation angle, lat is the latitude, δ is the solar declination angle, and t is the hour angle. θ and Ω were calculated for every pixel on different DOYs according to Eqs. 4 − 6, and C was calculated by linear regression between ρ and cos i in Eqs. 7 − 8. According to Eqs. 3 − 6, we deduced that i can be described by the following function (Eq. 9)

(9) $i=f(DOY,T,\beta ,\alpha ,lat,lon)$(9)

where lon is the longitude, and T is the passing time of Landsat satellites. For brevity, only six representative groups of passing-time parameters on DOYs 30, 60, 90, 300, 330 and 360 were used to compute solar angles for simulating topographical effects.

We randomly generated 50,000 sampling points in forested areas of the HDMR. The values of cosi and ρ at these points were calculated for monthly cloud-free composite images generated in section 2.2.1. Then, they were used to determine the C values of the NIR and SWIR2 bands on DOYs 30, 60, 90, 300, 330 and 360 through Eqs. 7 − 8. In our study, the median NIR and SWIR2 reflectance (0.1771 for NIR, 0.0521 for SWIR2) of stable forest samples occurring in the areas with slopes below 5° were assumed to be flat surface reflectance (ρ_n). The i for each DOY was simulated according to Eq. 3. The slope surface reflectance for different DOYs was computed by Eq. 10.

(10) $\rho ={\rho }_n\frac{cosi+\mathrm{C}}{cos\alpha cos\theta +\mathrm{C}}$(10)

Subsequently, the simulated NIR and SWIR2 reflectance were further used to calculate the NBR values for different DOYs. The abnormal NBR values that were > 1 or < −1 were replaced with the surrounding maximum or minimum pixel values within a window of 3 × 3 pixels, respectively. The NBR difference images between each pair of different DOYs (i.e. DOYs 30, 60, 90, 300, 330 and 360) were then generated; as a result, a total of 15 difference images were obtained. For each difference image, the pixels with an absolute value greater than 0.277 ($\left|\mathrm{\Delta }\mathrm{N}\mathrm{B}\mathrm{R}\right|\ge$0.277) were identified as topography-induced disturbances.

Four regions covering 1°×1° (Figure 1) were selected to examine the relations of topography-induced disturbances with latitude. For each region, 150,000 forest sampling points were randomly generated. Then, the propensity score matching method (Rosenbaum and Rubin 1983) was applied to select points with similar slopes and aspects from the four different regions, which was performed using the MatchIt package in RStudio (Ho et al. 2011; Racine 2012). Five matching results, each of which included 3000 similar points of every region, were produced. Subsequently, one-way analysis of variance (ANOVA) was applied to test whether the differences in topography-induced disturbances in different latitude zones were statistically significant.

As reported in existing studies (Chance et al. 2016; Kennedy et al. 2015; Vogeler et al. 2020), filtering can effectively remove the spatial noise of disturbance results. We filtered the forest disturbance patches mapped in this study as less than different patch sizes varying from one to ten pixels considering 8-connectivity. The process of filtering was coded and performed in GEE by calling the function “connectPixelCount” (Tamiminia et al. 2020). Assuming that the minimum mapping unit (MMU) of the prefiltered disturbance map was one, the MMU of the resultant disturbance map filtered with a patch size of one was changed to two.

2.5 Accuracy assessment

We used the percentages of topography-induced disturbances that were consistent or inconsistent with actual forest disturbances to evaluate the accuracy of topography-induced disturbance simulation in this study. The traditional confusion matrix method was used to assess the accuracy of forest disturbance detection. Moreover, considering the small area percentages of topography-induced and actual disturbances in the study area (Li et al. 2021), the area-adjusted confusion matrix method (Equation 11(11) ${\hat{p}}_{ij}=W_i\frac{n_{ij}}{n_{i\bullet }}$(11) ) advised by Olofsson et al. (2013) was also used to compute the accuracies of the forest disturbance map following Equations 12(12) $\hat{O}=\sum _{j=1}^q{\hat{p}}_{jj}$(12)  − 14(14) ${\hat{U}}_i=\frac{{\hat{p}}_{ii}}{{\hat{p}}_{i\bullet }}$(14) . The accuracy measures were presented with 95% confidence intervals, which were obtained from the area variances (Equations 15(15) $\hat{V}(\hat{O})=\sum _{i=1}^q{W}_i^2{\hat{U}}_i\left(1-{\hat{U}}_i\right)/\left(n_{i\bullet }-1\right)$(15)  − 17(17) $\hat{V}(\hat{U})={\hat{U}}_i\left(1-{\hat{U}}_i\right)/\left(n_{i\bullet }-1\right)$(17) ) in the disturbed and stable forests.

(11) ${\hat{p}}_{ij}=W_i\frac{n_{ij}}{n_{i\bullet }}$(11)

where W_i is the area proportion of class i, n_ij are cell entries in the common confusion matrix, and n_i· is the total number of sample units of class i in the evaluated map.

(12) $\hat{O}=\sum _{j=1}^q{\hat{p}}_{jj}$(12)

(13) ${\hat{P}}_j=\frac{{\hat{p}}_{jj}}{{\hat{p}}_{\bullet j}}$(13)

(14) ${\hat{U}}_i=\frac{{\hat{p}}_{ii}}{{\hat{p}}_{i\bullet }}$(14)

where $\hat{O}$, $\hat{P}$, and $\hat{U}$ represent the overall accuracy, producer’s accuracy, and user’s accuracy from the area-adjusted confusion matrix of q classes, respectively.

(15) $\hat{V}(\hat{O})=\sum _{i=1}^q{W}_i^2{\hat{U}}_i\left(1-{\hat{U}}_i\right)/\left(n_{i\bullet }-1\right)$(15)

(16) $\hat{V}\left(\hat{P}\right)=\frac{1}{{\hat{N}}_j}\left[\frac{N_j^2{\left(1-{\hat{P}}_j\right)}^2{\hat{U}}_j\left(1-{\hat{U}}_j\right)}{n_j-1}+{\hat{P}}_j^2\sum _{i\ne j}^q{N}_i^2\frac{n_{ij}}{n_i}\left(1-\frac{n_{ij}}{n_i}\right)/\left(n_i-1\right)\right]$(16)

(17) $\hat{V}(\hat{U})={\hat{U}}_i\left(1-{\hat{U}}_i\right)/\left(n_{i\bullet }-1\right)$(17)

where ${\hat{N}}_{\bullet j}=\sum _{i=1}^q\frac{N_{i\bullet }}{n_{i\bullet }}{n}_{ij}$ is the estimated number of reference class j, n_j· is the total number of sample units of class j in the evaluated map, and N_i· and N_j· are the number of class i and j in the evaluated map, respectively.

3. Results

3.1 Accuracies and spatiotemporal patterns of topography-induced disturbances simulated by the SCS+C model

Among all sampling blocks of LandTrendr-detected forest disturbances overlapping with the topography-affected disturbances obtained from simulation, 86.39% of blocks were stable forests but wrongly identified as disturbed forests due to topographic effects (Table 1), and only 12.37% of blocks were verified to be actual forest disturbances. However, approximately 1.24% of the sampled blocks failed to be used to determine whether forest disturbances actually occurred due to the lack of reliable reference data.

Table 1. The percentages of actually stable and disturbed forested areas in the overlapping areas between topography-induced disturbances simulated by the SCS+C model and LandTrendr-detected forest disturbances.

	Actual stable forests (correctly simulated but wrongly detected)	Actual disturbed forests (correctly detected but wrongly simulated as topography-induced disturbances)	Uncertain due to the lack of high spatial resolution images
Number of blocks	489	70	7
Percentage (%)	86.39	12.37	1.24

Topography-induced disturbances widely occurred in the HDMR, especially in the northeastern, eastern, and western regions (Figure 3). Larger areas of topography-induced disturbances were observed on DOYs 60 − 360 (3.72%), 90 − 360 (5.02%), 300 − 360 (4.65%), and 330 − 360 (3.40%). The difference in areas of topography-induced disturbances in the same DOY interval lengths was considerable, such as on DOYs 30 − 60 (0.98%) and 330 − 360 (3.40%). The area percentages of topography-induced disturbances were lower than 1% for the DOY pairs among DOY30, DOY60, and DOY300 but substantially increased to 2.75%−5.02% for the DOY pairs among DOY90, DOY330, and DOY360.

Figure 3. Spatiotemporal pattern of all topography-induced disturbances simulated by the SCS+C model in the forested regions of the HDMR and zoomed-in images of two typical sites (D and E). DOYs i−j indicate the potential disturbances caused by the differences in topographic effects between DOY i and DOY j.

Topography-induced disturbances mainly occurred on shady slopes (300°−360°) with a slope of 30°−55°, and their areas also varied considerably with different DOY pairs (Figure 4). In summary, the areas of topography-induced disturbances peaked at a slope aspect of approximately 330°. Variations in topography-induced disturbances with slopes were relevant to the examined DOY pairs. Larger areas of topography-induced disturbances occurred on slopes of 30°−45° on DOYs 60 − 360, 90 − 330, 90 − 360, 300 − 330, 300 − 360, and 330 − 360 and on slopes of 40°−55° on DOYs 30 − 60, 30 − 90, 30 − 300, 30 − 330 and 60 − 330. Moreover, no topography-induced disturbances were observed on slopes of less than 30°.

Figure 4. Variations in the areas of all topography-induced disturbances simulated by the SCS+C model with slopes and aspects in the forested regions of the HDMR. DOYs i−j indicate the potential disturbances caused by the differences in topographic effects between DOY i and DOY j.

As illustrated in Figure 5, the implementation of the propensity score matching method substantially lowered the absolute values of the standardized mean differences (SMDs) of slopes and aspects between the examined regions to approximately 0 (Figure 5a). Accordingly, no statistically significant (p>0.05) differences in slope and aspect existed for the points selected from different regions, except for one matching result (Figures 5b,c). After removing the effects of the differences in slope and aspect, the number of topography-induced disturbances varied significantly (p<0.001) with latitude, and more topography-induced disturbances occurred in higher latitude regions (Figure 5d).

Figure 5. Results of propensity score matching and one-way ANOVA. The standardized mean differences (SMDs) of slopes and aspects between the examined regions and the t-statistics and p value before and after propensity score matching are presented in (a), (b) and c), respectively; (d) illustrates the number of topography-induced disturbances and the significance tests using one-way ANOVA of the five matching results. Regions I−IV are presented in Figure 1. *** means p < 0.001.

3.2 Forest disturbances before and after excluding topography-induced disturbances

Table 2 shows that the traditional and area-adjusted OAs of forest-disturbed mapping increased from 88.82% and 93.67%±0.12% to 92.32% and 94%.29 ± 0.12%, respectively, after excluding topography-induced disturbances. Following the removal of topography-induced disturbances, the traditional and area-adjusted UAs of disturbed forest mapping increased from 79.08% and 79.08%±0.30% to 87.65% and 87.65%±0.26%, respectively, whereas the traditional and area-adjusted PAs of disturbed forest mapping decreased weakly from 90.13% and 52.43%±0.61% to 89.45% and 51.12%±0.60%, respectively. Accordingly, the traditional and area-adjusted PAs of stable forest mapping increased considerably from 88.17% and 98.41%±0.02% to 93.75% and 99.18%±0.02%, respectively, with a slight decrease (0.03%) in the UAs of stable forest mapping.

Table 2. Accuracy results of forest disturbances in the HDMR from 1989 to 2020 before and after excluding topography-induced disturbances.

Accuracy assessment method		Before excluding topography-induced disturbances		After excluding topography-induced disturbances
		Disturbed	Stable	Disturbed	Stable
Traditional	PA (%)	90.13	88.17	89.45	93.75
	UA (%)	79.08	94.74	87.65	94.71
	OA (%)	88.82		92.32
Area-adjusted	PA (%)	52.43±0.61	98.41±0.02	51.12±0.60	99.18±0.02
	UA (%)	79.08±0.30	94.74±0.13	87.65±0.26	94.71±0.12
	OA (%)	93.67±0.12		94.29±0.12

Note: OA, PA, and UA represent overall accuracy, producer’s accuracy, and user’s accuracy, respectively.

After removing topography-induced disturbances, the spatial patterns of forest disturbances changed greatly in the HDMR (Figure 6). Zoomed-in images further show that topography-induced disturbances were clumped or scattered (Figure 6d). Most topography-induced disturbances occurred in stable forests (F−H in Figures 6d,e). However, a few actual forest disturbances were incorrectly removed (I in Figures 6d,e).

Figure 6. Forest disturbances in the HDMR from 1989 to 2020 before and after excluding topography-induced disturbances. (a) LandTrendr-detected forest disturbances, (b) topography-induced disturbances generated from the SCS+C simulation, (c) the resultant forest disturbances after excluding topography-induced disturbances, (d) magnified images of four typical sites (F, G, H, and I) in (a), and (e) LandTrendr temporal segment trajectories of representative points (+) of four typical sites in (d). Topography-induced disturbances in F, G, and H were successfully removed from the LandTrendr-detected forest disturbances, whereas actual forest disturbances in I were incorrectly removed.

3.3 Forest disturbances filtered by different patch sizes

A sharp decrease in areas of forest disturbances occurred with the increase in filtered patch size from 0 to 3 (Figures 7a,b). Without filtering, approximately 10.43% (11.75$\times$10² km²) of the total area (112.66$\times$10² km²) of LandTrendr-detected forest disturbances were identified as topography-induced disturbances (Figures 7a,b); following spatial filtering, the area percentage of disturbed forests contributed by topography-induced disturbances decreased with increasing filtered patch sizes. For a filtered patch size of 10, the area percentage dropped to only 2.38% (Figures 7a,b). Additionally, with the removal of small-sized, scattered disturbance patches by the filtering process, the OAs of the disturbed forest decreased from 93.67%±0.12% to 91.94%±0.14% with increasing filtered patch sizes (Figure 7a). A similar trend was observed in PAs (from 52.43%±0.61% to 22.02%±0.37%) of disturbed forests (Figure 7a). Conversely, the UAs of disturbed forests increased from 79.08%±0.30% to 94.41%±0.20% with increasing filtered patch sizes (Figure 7a). Figure 8 further shows that the filtering process successfully removed the small-sized, scattered disturbance patches but retained the large-sized, lumped disturbance patches (Figures 8F-H), which were most relevant to topography-induced disturbances (Figures 6cF-H). Accordingly, some small-sized, scattered actual disturbance patches were also removed by the filtering process (Figure 8l).

Figure 7. Variations in areas and accuracies of forest disturbances generated by filtering LandTrendr-detected forest disturbances with different patch sizes before and after excluding topography-induced disturbances. (a) LandTrendr-detected forest disturbances, (b) forest disturbances after excluding topography-induced disturbances, and (c) forest disturbances occurring in the regions with topography-induced disturbances. Filtered patch size 0 indicates no filtering; filtered patch size i (i.e. 1–10) indicates that disturbance patches containing the number of pixels with an 8-connected rule greater than i were retained. OA, PA, and UA represent the overall accuracy, producer’s accuracy, and user’s accuracy based on the area-adjusted confusion matrix method, respectively.

Figure 8. Zoomed-in images of forest disturbances generated by filtering with different patch sizes. Four typical sites (F, G, H, and I) are marked in Figure 6a. Filtered patch size 0 indicates no filtering; filtered patch size i (i.e. 1–10) indicates that disturbance patches containing the number of pixels with an 8-connected rule greater than i were retained.

4. Discussion

4.1 Necessity of topographic correction for forest disturbance detection based on yearly wide-time-window LTSs

Our results showed that approximately 10% of the total area of LandTrendr-detected forest disturbances was contributed by topography-induced disturbances simulated by the SCS+C model, and the traditional and area-adjusted OAs of LandTrendr-based forest disturbance detection respectively increased by 3.50% and 0.62% following the exclusion of topography-induced disturbances (Table 2). The SCS+C model exhibited high performance in determining topography-induced disturbances, with an accuracy of approximately 86% (Table 1). These results demonstrated the substantial impacts of topographic effects on yearly wide-time-window LTSs-based forest disturbance detection in complex terrain mountainous regions. The improvement in the accuracy of forest change detection after excluding topographic effects was supported by the results of previous studies (Shimizu et al. 2016; Tan et al. 2013; Van den Hoek et al. 2021). Topographically corrected bitemporal Landsat image pairs considerably improved the OA of forest change detection by over 10% in mountainous areas with mean slopes of 11.2°−13.0° (Tan et al. 2013); topographic correction led to a good OA improvement (3.6%) in forest change detection from trajectory-based analysis of LTSs in mountainous areas with a mean slope of 11.4° and approximately 40% of slopes steeper than 30° (Shimizu et al. 2016); the inclusion of a topographic correction into the classification improved the classifier accuracy by an average of 1.15%, producing more conservative estimates of forest change in a predominantly mountainous country (Nepal) (Van den Hoek et al. 2021), which had important implications for improving the accuracies of postclassification change detection results. Although the OA of the resultant map excluding topography-induced disturbances is similar to our recent study (Li et al. 2021), the traditional UA and PA of disturbed forests were more balanced in this study. The area-adjusted UAs and PAs of the results with and without topography-induced disturbances were individually increased by 2.26% and 0.95%, and 1.69% and 10.26% compared with the previous study (the UA and PA were 77.39%±0.32 and 50.17%±0.88). The UAs and PAs of stable forests were more balanced in this study. Therefore, more high-quality forest disturbance results were achieved in the present study by excluding topographical effects. However, other studies argued that topographic correction contributed little improvement or even adversely affected the accuracies of forest disturbance detection, the consistency of LTSs, and fractional tree cover (Fcover) prediction (Adhikari et al. 2016; Chance et al. 2016; Qiu et al. 2019). Chance et al. (2016) reported that there was a weak difference in the forest change detection results before and after topographic correction in mountainous regions with a mean slope of 20°; Qiu et al. (2019) found that topographic correction had very limited or even negative impacts on the consistency of LTSs at five sites with a remarkable elevation gradient; and Adhikari et al. (2016) showed that topographic correction failed to improve the prediction of Fcover using ratio-based vegetation indices in tropical mountains. Moreover, numerous studies asserted that topographic correction required little deliberation for LTSs-based change detection using normalized difference spectral indices, such as NDVI and NBR, because these spectral indices can effectively remove or mitigate topographic effects in the LTSs (Banskota et al. 2014; Liang et al. 2014; Shimizu et al. 2016). However, our results demonstrated that topographic effects still exist when using NBR to detect disturbances based on LTSs. Considering that non-normalized difference spectral indices are more sensitive to topographical effects, the postprocessing topographical correction method proposed in this study could more considerably improve the performance of LTSs change detection based on these indices. The widespread uncertainties of topographic overcorrections and undercorrections were also widely used reasons to reject the implementation of topographic correction in LTSs-based change detection (Hantson and Chuvieco 2011; Yin et al. 2022). Rather than directly applying a topographic correction algorithm in the LTSs, we quantified topography-induced disturbances in the wide-time-window PBIC-generated NBR LTSs using a pixel-by-pixel simulation with the SCS+C model and then overlaid these disturbances with the LandTrendr-detected disturbances to produce the final forest disturbances. Thus, the effects of topographic overcorrections and undercorrections were avoided, thereby substantially increasing the accuracy of forest disturbance detection in this study area. Our findings shed light on the necessity of excluding topographic effects from the yearly wide-time-window LTSs for mapping long-term forest disturbances in areas with rugged terrain and steep slopes.

The results also indicated that filtering considerably improved the UA of disturbed forest mapping by removing small-sized, scattered topography-induced disturbances. After filtering, the area percentage of disturbed forests contributed by topography-induced disturbances decreased greatly from 10.43% (without filtering) to 2.38% (filtering with a patch size of 10) (Figures 7a,c). This effectiveness of filtering spatial noises was demonstrated by previous studies, which individually applied MMUs with patch sizes of 4 (Grinand et al. 2013), 9 (Masek et al. 2013), and 11 (Kennedy et al. 2015; Vogeler et al. 2020) to filter the disturbance detection results. However, filtering failed to remove the topography-induced disturbances with large-sized, lumped patches (Figure 8F-H). Moreover, incorrect removals of real forest disturbances (Figure 8l) induced by filtering led to a substantial decrease in the PA of disturbed forest mapping by approximately 30% (Figure 7a), and the OA of forest disturbance detection substantial decreased by approximately 2.7% (Figure 7a). Therefore, filtering postprocessing cannot substitute the inclusion of topographic correction into forest disturbance mapping based on yearly wide-time-window LTSs.

4.2 Spatiotemporal patterns of topographic effects

The yearly LTSs consisting of images acquired on DOYs 30, 60, and 300 produced the smallest topography-induced disturbances (all less than 1%), followed by those composed of the images acquired on DOYs 90 and 330 (Figure 3); however, the yearly LTSs containing the images acquired on DOY 360 yielded relatively more topography-induced disturbances (Figure 3). These illumination effects resulted from different image acquisition time points and led to large discrepancies in the detection accuracies of forest disturbances in mountain areas, as reported by previous studies based on yearly wide-time-window LTSs (Griffiths et al. 2014; Langner et al. 2018; Shimizu et al. 2016; Shimizu, Ota, and Mizoue 2019; Tang et al. 2019). The results of Tang et al. (2019) and Griffiths et al. (2014) produced by the composited yearly LTSs excluding the images acquired in December presented higher OAs of mapping forest disturbance, with values of 89.72%±0.67% and 86%, respectively; conversely, the results of Shimizu et al. (2016), Shimizu, Ota, and Mizoue (2019) and Langner et al. (2018) derived from the yearly wide-time-window LTSs containing the images acquired in December exhibited relatively low OAs of detecting forest disturbances, with values of 73.3%, 78.3%, and 78%, respectively. The simulated area percentages of topography-induced disturbances between DOYs 330 − 360 (higher solar zenith and azimuth angle) and between DOYs 30 − 60 (lower solar zenith and azimuth angle) were 3.40% and 0.98%, respectively (Figure 3); however, the values on wider intervals, such as on DOYs 30 − 300 (0.90%), remained low (Figure 3), which suggested that the areas of topography-induced disturbances were unrelated to the interval length of DOY pairs but depended on the differences in solar geometry between DOY pairs. Moreover, lower solar elevation angles at higher latitude regions led to more topography-induced disturbances (Figure 5d). These results agreed with the theory of illumination effects (Adhikari et al. 2016; Teillet, Guindon, and Goodenough 1982; Yin et al. 2022). Our findings highlighted the importance of selecting images with low illumination effects to composite the wide-time-window LTSs for detecting forest disturbances in mountainous regions with frequent cloud cover, such as Southeast Asia (Grogan et al. 2015).

The overwhelming majority of topography-induced disturbances occurred on shady slopes (300°−360°) with sloping angles above 30° (Figure 4), where topographic effects were considerable due to highly varying illumination (Soenen, Peddle, and Coburn 2005) and needed to be removed to detect forest disturbances based on yearly wide-time-window LTSs. Similar results of the topographic effects of steep slopes contributing to more disturbance detection errors were reported by previous studies (Shimizu et al. 2016; Tan et al. 2013). The critical sloping angle related to topography-induced disturbances was far higher than 10° (Hall-Könyves 1987) and approximately 20° (Kumar, Skidmore, and Knowles 1997), where topography-induced reflectance variations were substantial, as proposed by previous studies. This result showed that the NBR used for detecting forest disturbances in this study can partly offset the multiplication factors of topographic effects, which led to negligible topography-induced disturbances on slopes below 30°. Therefore, topographic effects can be neglected when using normalized difference spectral indices, such as NBR, for detecting forest disturbances in mountainous areas with a maximum slope of 30°.

Although the SCS+C model’s good performance was demonstrated in forested regions (Soenen, Peddle, and Coburn 2005), it failed to correct the effects of cast shadows (Li et al. 2016), which exist widely in dramatically undulating mountainous regions, such as the HDMR (Figure 1). Following topographic corrections using the SCS+C model, overcorrections or undercorrections still exist because of widespread cast shadows (Li et al. 2016). Moreover, the quality of the DEM affects the performance of DEM-dependent topographic correction methods. Although it reportedly performed better than other global DEMs for removing topographic effects from Landsat images (Adhikari et al. 2016; Pimple et al. 2017; Wu, Jin, and Fan 2016), the SRTM DEM has a far coarser spatial resolution than the theoretically required value, which should be finer than one-third of the pixel size of the image to be corrected (Hantson and Chuvieco 2011; Wu, Jin, and Fan 2016; Yin et al. 2022). The low vertical precision of the SRTM DEM also negatively affects the effectiveness of topographic correction (Elkhrachy 2018). Unfortunately, high-quality DEMs with a fine spatial resolution and a high vertical precision have poor availability, especially in forested mountainous regions (Adhikari et al. 2016; Pimple et al. 2017). Recently, commercial global DEMs with higher spatial resolution, such as AW3D (5 m), NEXTMAP 5™ (5 m), and WorldDEM™ (12 m), have become available. With the further development of earth observation platforms and technologies, the gap in vertical accuracy and spatial resolution between global DEMs and local DEMs is being closed (Li et al. 2018; Tang et al. 2021), which will boost the high-precision removal of topographic effects in forest disturbance detection over mountainous areas.

5. Conclusion

Accurately mapping forest disturbances in mountainous regions has implications for the assessment of carbon targets, biodiversity loss and water security. However, uncertainties caused by topographic effects inherent in yearly wide-time-window LTSs negatively affect the performance of forest disturbance detection. In this study, we conducted a quantitative and pixel-by-pixel simulation of topography-induced forest disturbances by integrating the SCS+C topographic correction model with a LandTrendr time-series change detection algorithm, taking the HDMR of southwest China as a case study. Approximately one-tenth of the forest disturbance areas detected by the LandTrendr algorithm were topography-induced. The removal of these false disturbances resulting from topographic effects led to a considerable increase in the traditional (3.50%) and area-adjusted (0.62%) OAs of the forest disturbance map. This result indicated that the inclusion of the simulation of topographic effects in the LandTrendr-based change detection procedure performed well and was necessary. Topography-induced disturbances mainly occurred on northwest-facing (300°−360°) slopes with gradients of 30°−55° and higher latitudes. Their occurrences relate to larger differences in the solar zenith and azimuth angle between image pairs of the LTSs to be detected. Our study confirmed that filtering spatial noise had limited impacts on improving the accuracy of forest disturbance mapping because it failed to remove topography-induced disturbance patches with large areas.

In this paper, a postprocessing approach with a pixel-by-pixel simulation of topography-induced disturbance was proposed to improve LandTrendr-based forest disturbance detection. Its good performance implies a promising potential for replacing the preprocessing method of correcting the topographical effects of images to improve the reflectance consistency of the LTSs before performing a change detection algorithm. Although high accuracies in forest disturbance mapping have been obtained, the postprocessing approach proposed here was DEM-dependent. With the free availability of global DEM data with higher spatial resolution and more precise vertical accuracy, high-precision removal of topographical effects from forest disturbances detected by yearly wide-time-window LTSs in mountainous areas using our proposed approach can be expected.

Acknowledgments

This research was jointly supported by the Second Tibetan Plateau Scientific Expedition and Research Program (grant number: 2019QZKK0402), the National Natural Science Foundation of China (grant number: 41971239), and the Scientific Research Foundation of Education Department of Yunnan Province (grant number: 2022Y061).

Disclosure statement

The authors declare that no known competing financial interests or personal relationships could have affected the work reported in this paper.

Data availability statement

The China Land-Use/Cover datasets are available in the Resource and Environment Science and Data Centre at https://www.resdc.cn/data.aspx?DATAID=197. All the Landsat images and SRTM DEM that support the findings of this study are available in the Google Earth Engine platform at https://developers.google.com/earth-engine/datasets/catalog/landsat and at https://developers.google.com/earth-engine/datasets/catalog/USGS_SRTMGL1_003?hl=en, respectively.