High-dimensional detection of Landscape Dynamics: a Landsat time series-based algorithm for forest disturbance mapping and beyond

ABSTRACT

Time series analysis of medium-resolution multispectral satellite imagery is critical to investigate forest disturbance dynamics at the landscape scale. In particular, the spatial, temporal, and radiometric consistency of Landsat time series data provides unprecedented insight into past disturbances that occurred over the last four decades. Several Landsat time series-based algorithms have been developed to automate the detection of forest disturbances. However, automated detection of non-stand-replacing disturbances based on Landsat time series remains a challenging task due to the difficulty of effectively separating them from spectral noise. Here, we present the High-dimensional detection of Landscape Dynamics (HILANDYN) algorithm, which exploits spatial and spectral information provided by Landsat time series to detect forest disturbance dynamics retrospectively. A novel and unsupervised procedure for changepoint detection in high-dimensional time series allows HILANDYN to perform the temporal segmentation of inter-annual time series into linear trends. The algorithm embeds a noise filter to remove spurious changepoints caused by residual spectral noise in the time series. We tested HILANDYN to detect disturbances that occurred in the forests of the European Alps over a period of 39 years, i.e. between 1984 and 2022, and evaluated its accuracy using a validation dataset of 3000 plots randomly located inside and outside the disturbed patches. We compared HILANDYN with the Bayesian Estimator of Abrupt change, Seasonality, and Trend (BEAST), which is a well-established and high-performing time series-based algorithm for changepoint detection. The quantitative results highlighted that the number of bands, i.e. original Landsat bands and spectral indices, included in the high-dimensional time series and the threshold controlling the significance of changepoints strongly influenced the user’s accuracy (UA). Conversely, changes in the combinations of bands primarily affected the producer’s accuracy (PA). HILANDYN achieved an F₁ score of 0.801, which increased to 0.833 when we activated the noise filter, allowing the algorithm to balance UA (83.1%) and PA (83.5%). The qualitative results showed that disturbed forest patches detected by HILANDYN were characterized by a high spatio-temporal consistency, regardless of the disturbance severity. Furthermore, our algorithm was able to detect forest patches associated with secondary disturbances, such as salvage logging, that occur in close succession with respect to the primary event. The comparison with BEAST evidenced a similar sensitivity of the algorithms to non-stand-replacing events, as both achieved comparable PA. However, BEAST struggled to balance UA and PA when using a single parameter set, achieving a maximum F₁ score of 0.717. Moreover, the computational efficiency of BEAST in processing high-dimensional time series was very limited due to its univariate nature based on the Bayesian approach. The adaptability of HILANDYN to detect a wide range of disturbance severities using a single parameter set and its computational efficiency in handling high-dimensional time series promotes its scalability to large study areas characterized by heterogeneous ecological conditions.

KEYWORDS:

• European Alps
• Forest disturbance
• high-dimensional analysis
• changepoint analysis
• temporal segmentation
• spatial context
• Natural disturbances

1. Introduction

Disturbances such as wildfires, windthrows, and insect or pathogen outbreaks are key components of forest ecosystems (Turner 2010). They alter the state and trajectories of ecosystems and generate heterogeneity in space and time. The intensification of climate-driven changes in disturbance regimes has been observed globally, with negative impacts on the structure, functions, and composition of forests (Forzieri et al. 2021; Seidl et al. 2017). Shifts in disturbance patterns toward more frequent, larger and more severe events can drastically impair the resilience of forests and trigger transitions to a non-forest state (Johnstone et al. 2016; Seidl et al. 2016). Besides stand-replacing disturbances, long-term studies based either on field (Andrus et al. 2021; Camarero et al. 2015) or remote sensing data (Cohen et al. 2016; Senf et al. 2018) highlighted increasing trends in tree mortality rates, associated with gradual, non-stand-replacing disturbances. Warming temperatures and increases in severe and prolonged droughts driven by climate change resulted in extensive forest decline globally (Allen, Breshears, and McDowell 2010, 2015). Apart from the direct effects of climate change, modifications of the interactions between disturbance agents (Seidl et al. 2017) and land-use change patterns, such as land abandonment (Mantero et al. 2020), further drove shifts in disturbance regimes. Identifying changes in disturbance regimes that occurred over the last decades has crucial implications for current and future management policies (Leverkus et al. 2021).

Satellite-based remote sensing is a fundamental data source for forest ecologists as it enables a comprehensive understanding of ecological processes both in the spatial and temporal dimensions (Kennedy et al. 2014; Senf 2022). In particular, continuously acquired data by Landsat satellites over the last four decades provide a unique opportunity to reconstruct long-term forest disturbance dynamics at the landscape scale (Senf 2022; Wulder et al. 2019). The opening of the U.S. Geological Survey (USGS) Landsat archive in 2008 prompted a rapid increase in the use of time series for analyzing forest ecosystem dynamics from regional to global scales (Banskota et al. 2014; Wulder et al. 2012).

Several Landsat time series-based change detection algorithms have been developed so far, either to explicitly map forest dynamics or target a broader range of land cover changes (Zhu 2017). Some of these algorithms, such as LandTrendr (Landsat-based Detection of Trends in Disturbance and Recovery; Kennedy, Yang, and Cohen 2010) and CCDC (Continuous Change Detection and Classification; Zhu and Woodcock 2014), gained considerable popularity, thanks also to their implementation in the Google Earth Engine platform (Pasquarella et al. 2022). Both LandTrendr and CCDC perform a temporal segmentation, i.e. they partition time series into discrete intervals based on the similarity between consecutive observations, at the pixel level (Pasquarella et al. 2022). Whilst, based on a similar theoretical framework, these algorithms present fundamental differences in the characteristics of the input time series and the mathematical approaches they rely on. LandTrendr is an offline algorithm, i.e. it operates on completely available time series, that segments inter-annual data into linear trends. Conversely, CCDC is an online algorithm, i.e. it processes continuously acquired data, that fits harmonic functions for modeling intra-annual phenological cycles. LandTrendr is univariate as it can segment one band, i.e. original reflectance bands and spectral indices, sensu Cohen et al. (2020), at each run, while CCDC is multivariate as it simultaneously analyzes multiple bands (Zhu 2017). C2C (Composite2Change; Hermosilla et al. 2015) and VeRDET (Vegetation Regeneration and Disturbance Estimates through Time; Hughes, Kaylor, and Hayes 2017) are among the algorithms that share characteristics with LandTrendr, as they are offline, univariate, and segment inter-annual time series into linear trends. On the contrary, EWMACD (Exponentially Weighted Moving Average Change Detection; Brooks et al. 2014) and COLD (COntinuous monitoring of Land Disturbance; Zhu et al. 2019) are similar to CCDC, as they model phenological cycles using continuously acquired data.

The sensitivity of Landsat time series-based algorithms to non-stand-replacing disturbances, i.e. those causing partial canopy loss, has been reported to be somewhat limited compared, for example, to visual interpretation of spectral-temporal profiles (Cohen et al. 2017) or bi-temporal burn severity indices (Rodman et al. 2021). Indeed, extracting the weak signals associated with non-stand-replacing disturbances is hindered by noise in Landsat time series, which is generated by uncertainties in pre-processing operations like georeferencing and atmospheric correction (Cohen et al. 2017; Rodman et al. 2021; Vogelmann et al. 2016). A key issue when targeting non-stand-replacing disturbances is the proper optimization of the parameters of an algorithm in order to balance omission and commission errors (Cohen et al. 2017; Ye et al. 2021).

Various ensemble approaches have been proposed in the context of Landsat time series-based algorithms as they improved the accuracy of disturbance maps, primarily by increasing the sensitivity to low-severity events (Cohen et al. 2018; Healey et al. 2018; Schultz et al. 2016). Algorithm ensembles are built by aggregating results from multiple change detection algorithms that are executed in parallel (Cohen et al. 2020; Healey et al. 2018; Hislop et al. 2019). Conversely, multispectral ensembles (Cohen et al. 2018; Marzo et al. 2021; Schultz et al. 2016) involve the analysis of several bands through parallel runs of a univariate algorithm, which outputs are then dissolved into a single result using a supervised classifier, e.g. random forest. Considering individual algorithms, BEAST (Bayesian Estimator of Abrupt change, Seasonality, and Trend; Zhao et al. 2019) employs an ensemble learning technique termed Bayesian model averaging (BMA). BMA involves combining the estimates of many candidate models by quantifying their usefulness according to a Bayesian posterior probability. This approach allows BEAST to explicitly quantify model uncertainty and deal with complex nonlinear signals in the seasonal and trend components. BEAST exhibited solid performance in detecting both stand-replacing and non-stand-replacing disturbances using time series of medium-resolution satellite data, such as Landsat (Hu et al. 2021; Moreno-Fernández et al. 2021) and Sentinel-2 (Giannetti et al. 2021; Mulverhill, Coops, and Achim 2023).

Spectral information from the spatial context has been integrated into Landsat time series-based algorithms, either before or after temporal segmentation (Hermosilla et al. 2015; Hughes, Kaylor, and Hayes 2017; Meng et al. 2021). Moving window-based approaches have been employed before temporal segmentation to increase the accuracy of forest disturbance (Hamunyela, Verbesselt, and Herold 2016) and recovery (Meng et al. 2021) maps. Object-based image analysis (OBIA) is another effective approach for analyzing change processes in forest landscapes at different spatial scales (Gómez, White, and Wulder 2011). In this sense, the VeRDET algorithm (Hughes, Kaylor, and Hayes 2017) spatially segments individual bands into patches of constant values to reduce the heterogeneity within the same land cover class prior to perform temporal segmentation. Applying spatial context approaches after temporal segmentation is also effective for enhancing forest disturbance maps. The C2C algorithm (Hermosilla et al. 2015), for example, performs postprocessing operations based on moving windows to improve the spatial-temporal consistency of disturbance patches. In a similar fashion, spatial filtering operations allowed Senf and Seidl (2020) and Ye et al. (2021) to reduce omission and commission errors associated with small patches in the disturbance maps. Recently, Ye et al., (2023) have integrated OBIA into a COLD-like algorithm for identifying change objects, i.e. regions of homogeneous change in space and time, at different spatial scales.

A marked increase in the availability of statistical methods for changepoint detection in multivariate and high-dimensional time series, i.e. time series with dimensions of the same order of magnitude of the observations or even larger, has been observed in recent years (Cho et al. 2021; Truong, Oudre, and Vayatis 2020). To deal with the complexity associated with multivariate time series, one common approach is to aggregate a test statistic (Groen, Kapetanios, and Price 2013) or project the time series to a single dimension and then process the data with a univariate method (Wang and Samworth 2018). Assumptions on the proportion of time series that undergo a change determine the most appropriate aggregation method (Groen, Kapetanios, and Price 2013). For example, the widely used CUSUM statistic for detecting structural changes (Brown, Durbin, and Evans 1975) has been aggregated across time series using either the average or the maximum (Groen, Kapetanios, and Price 2013; Jirak 2015). Generally, the average of a test statistic performs better when a change occurs in the majority of the time series, i.e. a dense change, while the maximum is more appropriate for changes occurring only in a fraction of them, i.e. a sparse change (Groen, Kapetanios, and Price 2013; Jirak 2015). Changepoint methods have primarily focused on detecting changes in the mean associated with constant signals, though some approaches use piecewise polynomial models to segment time series (Cho et al. 2021).

Here, we aimed to develop a Landsat time series-based algorithm for retrospectively mapping forest disturbance and recovery dynamics that overcomes some of the major limitations posed by existing approaches. These include the ability to simultaneously analyze all the dimensions of Landsat data, i.e. spatial, temporal, and spectral, the need for extensive parametrization tuning, and the balancing of omission and commission errors while targeting non-stand-replacing events. The algorithm, named High-dimensional detection of Landscape Dynamics (HILANDYN), harnesses an offline and unsupervised temporal segmentation procedure for detecting changes in linear trends derived from high-dimensional time series containing multispectral and spatial context information. Though HILANDYN can detect changes associated with forest disturbances and recovery processes, here, we focused only on disturbance detection. In particular, we aimed to (i) assess the sensitivity of HILANDYN to stand-replacing and non-stand-replacing disturbances both quantitatively and qualitatively, and (ii) to compare the proposed approach with the BEAST algorithm.

The remainder of the article is organized as follows. In Section 2, we provide details on the study area and a technical description of the algorithm. In Section 3, we present quantitative and qualitative results on forest disturbance detection obtained with HILANDYN, including the comparison with a recent time series-based algorithm. In Section 4, we discuss the strengths and weaknesses of the proposed algorithm and how it compares to state-of-the-art change detection algorithms.

2. Materials and methods

2.1. Study area

Our study area encompasses the core mountainous area of the European Alps as defined by the Alpine Convention (Plangger 2020); its area is 190,706 km², of which 51% is covered by forests according to the Copernicus Dominant Leaf Type 2018 map (European Environmental Agency 2020; Figure 1). Elevation ranges between 0 and 4807 m a.s.l. (Mont Blanc), with an average altitude of 1293 m a.s.l. According to the official European biogeographical regions (European Environment Agency 2016), the study area includes Alpine, Continental, and Mediterranean regions. Climate ranges from oceanic to dry, with mean temperatures ranging from less than 0° to 10° (Nigrelli and Chiarle 2023) and annual precipitation from 400 to 3000 mm (Isotta et al. 2014). Forests dominated by broadleaved species occupy 45,088 km², while those dominated by conifers amount to 52,694 km², according to the Dominant Leaf Type 2018 map (European Environmental Agency 2020). Oaks (Quercus spp.) dominate the forests at lower elevations, while silver fir (Abies alba Mill.), beech (Fagus sylvatica L.), and Scots pine (Pinus sylvestris L.) dominate in the montane belt. The most common species in the subalpine belt are European larch (Larix decidua Mill.), Norway spruce (Picea abies (L.) H. Karst.), and Swiss stone pine (Pinus cembra L.).

Figure 1. Geographic location of the study area and validation plots with information from the Dominant Leaf Type 2018 forest cover layer (European Environmental Agency 2020) and Landsat Path/Row coverage (World Reference System 2). The base layer corresponds to the hill shade derived from the 25 m EU-DEM digital surface model.

Natural disturbances in the European Alps mainly include windthrows, wildfires, insect outbreaks, and snow avalanches (Bebi et al. 2017). Windthrows are by far the most important disturbances in terms of size and severity (Kulakowski et al. 2016) and are usually associated with winter storms. Wildfires, which occur primarily in the southern and central Alps, are largely driven by the land abandonment of rural areas, which has led to fuel accumulation and increased connectivity (Bebi et al. 2017). Bark beetle outbreaks are the most important biotic disturbances in the European Alps and typically occur shortly after a windthrow. A strong increasing trend in damage to timber volume caused by bark beetles has been observed in the Alps over the last decades due to favorable climatic conditions, such as warming temperatures and prolonged drought periods (Patacca et al. 2023).

2.2. Landsat data

Our Landsat dataset included Collection 2 imagery acquired in 30 scenes (Figure 1) between 1984 and 2022 by the Thematic Mapper (TM), Enhanced Thematic Mapper + (ETM+), and Operational Land Imager (OLI) sensors. We retrieved all Level 2 Tier 1 imagery with <80% cloud cover acquired between May and October from the U.S. Geological Survey EarthExplorer tool. We removed pixels contaminated by clouds, cloud shadows, and snow using the Quality Assessment (QA) band bundled with each Landsat image and produced through the C Function of Mask (CFmask) algorithm (Foga et al. 2017).

We created a total of 39 pixel-based annual reflectance composites using the weighted geometric median (Morresi et al. 2022; Roberts, Wilford, and Ghattas 2019) based on the six reflectance bands (Blue, Green, Red, NIR, SWIR1, and SWIR2) of the TM/ETM+ sensors and the corresponding OLI bands relative to the growing season, i.e. from June 1 to September 30. Following Morresi et al. (2022), the compositing time window iteratively widened up to 20 days on both of its sides until at least three clear observations were found in order to minimize data gaps. To reduce the presence of artifacts generated by residual clouds and cloud shadows in reflectance composites, we computed pixel-wise weights relative to the NDVI value ($w_{\mathit{NDVI}}$) and the Euclidean distance to clouds and cloud shadows ($w_{\mathit{DIST}}$) for each valid observation within the compositing time window. In particular, the compositing algorithm prioritized observations with a high NDVI, as these were typically not contaminated by clouds or cloud shadows (Qiu et al. 2023). Furthermore, $w_{\mathit{NDVI}}$ increased the importance of those observations acquired close to the peak of the growing season, when NDVI typically reaches its maximum. The inclusion of $w_{\mathit{DIST}}$ in the compositing algorithm limited the influence of undetected cloud shadows, which can exhibit high NDVI values due to low reflectance of both the Red and NIR. To compute $w_{\mathit{DIST}}$, we used the sigmoid function proposed by Frantz et al. (2017), which can assume values in the interval [0, 1] depending on the requested distance ($d_{\mathit{req}}$) parameter. We set $d_{\mathit{req}}$ equal to 1500 m, i.e. 50 Landsat pixels, as proposed by White et al. (2014).

We tested HILANDYN using the six Landsat bands included in the annual reflectance composites and eight spectral indices (Table 1).

Table 1. List of the spectral indices used in this study.

Index	Formulation	Direction of change caused by a disturbance	Reference
Normalized Difference Vegetation Index (NDVI)	(NIR – Red)/(NIR + Red)	decrease	(Rouse et al. 1973)
Normalized Difference Moisture Index (NDMI)	(NIR – SWIR1)/(NIR + SWIR1)	decrease	(Wilson and Sader 2002)
Normalized Burn Ratio (NBR)	(NIR – SWIR2)/(NIR + SWIR2)	decrease	(García and Caselles 1991)
Moisture Stress Index (MSI)	NIR/SWIR1	increase	(Hunt and Rock 1989)
Tasseled Cap Brightness (TCB)	0.3037*Blue+0.2793*Green+0.4743*Red+0.5585*NIR+0.5082*SWIR1+0.1863*SWIR2	increase	(Crist 1985)
Tasseled Cap Greenness (TCG)	−0.2848*Blue−0.2435*Green−0.4743*Red+0.7243*NIR+0.0840*SWIR1−0.1800*SWIR2	decrease	(Crist 1985)
Tasseled Cap Wetness (TCW)	0.1509*Blue+0.1973*Green+0.3279*Red+0.34065*NIR−0.7112*SWIR1−0.4572*SWIR2	decrease	(Crist 1985)
Tasseled Cap Angle (TCA)	arctan(TCB/TCG)	decrease	(Powell et al. 2010)

2.3. Overview of the algorithm

HILANDYN aims to map forest patches characterized by disturbance and recovery dynamics through a changepoint analysis based on linear trends derived from high-dimensional Landsat time series. Specifically, a procedure called High-dimensional Trend Segmentation (HiTS, Maeng 2019, publicly available at https://github.com/hmaeng/HiTS) is at the core of our algorithm. The HiTS procedure is a generalization to higher dimensions of the TrendSegment procedure (Maeng and Fryzlewicz 2023), which was developed to detect multiple changepoints in linear trends for univariate data sequences.

HILANDYN builds high-dimensional time series by extracting data from time-ordered sequences of raster cubes using a three-dimensional spatial kernel (Figure 2). The number of variables in the high-dimensional time series corresponds to the product between the number of pixels in the spatial kernel and the number of bands. Here, we employed a three-by-three kernel, i.e. nine connected Landsat pixels. Our preliminary tests showed that this kernel dimension was adequate to provide spatial context information when detecting forest dynamics with Landsat data. Other studies (Hermosilla et al. 2015; Meng et al. 2021) used the same kernel dimension with Landsat time series-based change detection algorithms. Although not specified in the HiTS procedure, we recommend that high-dimensional time series contain at least six valid observations for the algorithm to work properly.

Figure 2. Flowchart of the main operations performed by HILANDYN. Input data consist of four-dimensional raster cubes where bands and time form the third and fourth dimensions, respectively. A three-dimensional spatial kernel is employed on a pixel basis to build $\mathit{n}$-dimensional time series with a matrix form of the dimension $\mathit{n}\times \mathit{T}$. The core operations performed by HILANDYN on $\mathit{n}$-dimensional time series are described in the orange box.

Initially, HILANDYN removes gaps with one-year length from high-dimensional time series while discards those with longer gaps. It then employs an iterative procedure to discern between changepoints associated with forest dynamics and those caused by impulsive noise, i.e. outliers in the spectral signal that are caused, for example, by undetected clouds, cloud shadows, or haze (Hermosilla et al. 2015; Kennedy, Yang, and Cohen 2010, Figure 2). If changepoints and gaps co-occur in a multispectral time series, HILANDYN updates the corresponding year of change depending on the number of gaps preceding each changepoint. It imputes data gaps using values of linear segments previously estimated by the HiTS procedure. Specifically, gaps located at the vertices of segments of at least 3 years in length are imputed using linear extrapolation based on the two preceding values. Otherwise, gaps are filled using the preceding value. Gaps located in the middle of a segment are filled using linear interpolation.

HILANDYN classifies changepoints as disturbances or growth events if at least half of the bands (Table 1) in the focal pixel follow the pattern of the spectral change associated with that event. In particular, the algorithm evaluates only those bands whose original and estimated values feature the same direction of spectral change. The spectral change magnitude of each event corresponds to the median of the relative change magnitudes among bands, using all the pixels in the spatial kernel. Following Cohen et al. (2016), change magnitudes are relativized using pre-disturbance values.

The implementation of HILANDYN in the C++ language through the Application Programming Interface (API) provided by the Armadillo open-source linear algebra library (Sanderson and Curtin 2016) ensured computational efficiency even when the number of variables in high-dimensional time series was high, e.g. more than 100.

2.4. HiTS procedure

The changepoint model proposed in the HiTS procedure (Maeng 2019) considers a high-dimensional time series data containing $\mathit{n}$ variables of length $\mathit{T}$ as follows:

(1) $X_{\mathit{i},\mathit{t}}=f_{\mathit{i},\mathit{t}}+{\epsilon }_{\mathit{i},\mathit{t}},$(1)

where $f_i={\left({\mathit{f}}_{\mathit{i},1},\dots ,f_{\mathit{i},\mathit{T}}\right)}^{\mathrm{\top }}$ is the piecewise linear signal of the time series $X_i={\left({\mathit{X}}_{\mathit{i},1},\dots ,X_{\mathit{i},\mathit{T}}\right)}^{\mathrm{\top }}$and ${\epsilon }_i={\left({\mathit{\epsilon }}_{\mathit{i},1},\dots ,{\epsilon }_{\mathit{i},\mathit{T}}\right)}^{\mathrm{\top }}$ is the independent Gaussian random error with mean zero and variance ${\sigma }^2$. The model assumes that the signal vectors ${\left\{f_i\right\}}_{\mathit{i}=1}^{\mathit{n}}$ have the form of a piecewise linear function and share $\mathit{N}$ distinct changepoints at unknown locations ${\eta }_1,\dots ,{\eta }_N$ in the sense that at each changepoint, at least one signal vector undergoes a change in its linear trend, whether in the intercept or the slope or both. HiTS was designed to work well in detecting multiple changepoints corresponding to linear trend changes or point anomalies, which are large deviations of the signal from its neighboring segments. It consists of four main steps: (1) High-dimensional Tail-Greedy Unbalanced Wavelet (HiTGUW) transform, (2) thresholding, (3) inverse HiTGUW transformation, and (4) post-processing.

The standard deviation of the error, σ, can vary across the $\mathit{n}$ time series, and it is estimated using the median absolute deviation (Hampel 1974), which is adjusted for achieving asymptotic normal consistency (Equation 2(2) ${\hat{\sigma }}_{\mathit{i}}=\left(\mathit{MAD}\left(X_{\mathit{i},1}-2X_{\mathit{i},2}+X_{\mathit{i},3},\dots ,X_{\mathit{i},\mathit{T}-2}-2X_{\mathit{i},\mathit{T}-1}+X_{\mathit{i},\mathit{T}}\right)\text{break}\times 1.4826\right)/\sqrt{6}$(2) ).

(2) ${\hat{\sigma }}_{\mathit{i}}=\left(\mathit{MAD}\left(X_{\mathit{i},1}-2X_{\mathit{i},2}+X_{\mathit{i},3},\dots ,X_{\mathit{i},\mathit{T}-2}-2X_{\mathit{i},\mathit{T}-1}+X_{\mathit{i},\mathit{T}}\right)\text{break}\times 1.4826\right)/\sqrt{6}$(2)

At the beginning of the processing, the HiTS procedure normalizes each time series $X_i$ by its estimated standard deviation ${\hat{\sigma }}_{\mathit{i}}$. Values are transformed back to their original scales at the end of the processing.

The core ingredient of the HiTS procedure is the HiTGUW transform, which employs a bottom-up approach to construct a data-driven wavelet basis common to all the variables in the high-dimensional time series. Starting from the finest scale, i.e. using raw input data of dimension $\mathit{n}$ × $\mathit{T}$, the HiTGUW transform recursively merges neighboring regions of the data from bottom to top and is completed after $\mathit{T}$ − 2 local orthonormal transformations that result in a multiscale decomposition of the input matrix with $\mathit{n}$ × ($\mathit{T}$ − 2) detail-type coefficients and $\mathit{n}$ × 2 smooth-type coefficients. The detail coefficient plays an important role in deciding which region should be merged first, as its size indicates the strength of local linearity; the detail coefficient becomes zero only when the raw observations in the corresponding merged region have a perfect linear trend. The merges are performed by giving priority to the region whose aggregated detail coefficient has the smallest size, as they unlikely contain changepoints. The aggregation strategy of detail coefficients at a given time step is computed as the column-wise maximum value over $\mathit{n}$ variables. Thanks to its bottom-up approach, the HiTGUW transform encodes the bulk of the variance of the input data in only a few detail coefficients obtained at later merges, thus enabling the sparse representation of the data. This latter justifies thresholding as a way of deciding the significance of each aggregated detail coefficient, i.e. the strength of the local deviation from linearity.

In the thresholding step, the $\mathit{n}$ × ($\mathit{T}$ − 2) detail coefficient matrix produced by the HiTGUW transform is aggregated column-wise using the maximum value over the $\mathit{n}$ time series. A changepoint is detected when, at a given time step, the size of the aggregated detail coefficient is greater than a pre-specified threshold. The threshold used to detect significant deviations from linearity is computed as follows:

(3) $\mathit{\lambda }=\mathit{C}\sqrt{2\mathit{log}\left(\mathit{nT}\right)},$(3)

where $\mathit{n}$ is the number of time series, $\mathit{T}$ is their length and $\mathit{C}$ is a constant value.

Next, the inverse HiTGUW transformation performs inverted, i.e. transposed, orthonormal transformations in reverse order to that in which they were initially performed. This step uses the thresholded detail coefficients to produce the estimated piecewise linear signal composed of best linear regression fits (i.e. minimizing the sum of squared errors) for each estimated segment. Lastly, in the post-processing step, non-significant changepoints are removed by performing the first three steps of the HiTS procedure using the estimated functions computed by the inverse HiTGUW transformation as input data.

The outputs of the procedure include $\mathit{n}$ × $\mathit{T}$ estimated values of the high-dimensional time series, the position of changepoints in time, and a matrix of dimension $\mathit{n}$ × $\mathit{T}$ with the position of changepoints among $\mathit{n}$ variables. This latter is obtained by applying the threshold (Equation 3(3) $\mathit{\lambda }=\mathit{C}\sqrt{2\mathit{log}\left(\mathit{nT}\right)},$(3) ) to the detail coefficients of each variable in the high-dimensional time series. The reader is referred to Maeng (2019) for more details on the HiTS procedure.

2.5. Modifications of the HiTS procedure

The original HiTS procedure prioritizes the detection of high-magnitude and possibly ephemeral changes that are highly sparse across the time series (Maeng 2019). This characteristic is mainly related to the aggregation strategy of the detail coefficients in the HiTGUW transform and thresholding steps (Section 2.4). To detect forest patches characterized by disturbance and recovery dynamics, HILANDYN accounts for the similarity in spectral-temporal information from adjacent Landsat pixels. Indeed, within the spatial kernel, changes may occur in a subset of pixels, for example, when the focal pixel lies at the edges of a disturbed forest patch (Figure 3a) or belongs to a linear feature (Figure 3b–c). We have therefore modified the HiTS procedure to fully exploit the spectral information in the spatial neighborhood included in high-dimensional Landsat time series.

Figure 3. Examples of the spatial arrangement of disturbed pixels (orange) within the spatial kernel used by HILANDYN where the focal pixel is surrounded by non-disturbed pixels (grey).

First, in the HiTGUW transformation, we modified the aggregation strategy of the detail coefficients so that the merging order of the time intervals depends on the sum of their maximum and average. This procedure delays the merging of those time intervals that contain both large and dense deviations from linearity among the $\mathit{n}$ variables of the high-dimensional time series, thus favoring the detection of true changepoints in the thresholding step.

Second, in order to preserve the spatial pattern of disturbed patches, we embedded the spectral-temporal information provided by the neighboring pixels in the HiTGUW transformation. Specifically, the wavelet basis is constructed by prioritizing the detail coefficients of those neighboring pixels that are more similar to the focal pixel, taking into account the spectral and temporal characteristics. In particular, for each band, we computed the spectral angle mapper (SAM; Equation 4(4) $\mathit{SAM}\left(\mathit{p},\mathit{r}\right)=\mathit{co}{s}^{-1}\left(\frac{{\sum }_{i=1}^T{p}_i{r}_i}{\sqrt{{\sum }_{i=1}^T{p}_i^2}\sqrt{{\sum }_{i=1}^T{r}_i^2}}\right)$(4) ; Kruse et al. 1993) between vectors $\mathit{p}$ and $\mathit{r}$, corresponding to one neighboring and focal Landsat pixel, respectively, and having length $\mathit{T}$, i.e. the length of time series.

(4) $\mathit{SAM}\left(\mathit{p},\mathit{r}\right)=\mathit{co}{s}^{-1}\left(\frac{{\sum }_{i=1}^T{p}_i{r}_i}{\sqrt{{\sum }_{i=1}^T{p}_i^2}\sqrt{{\sum }_{i=1}^T{r}_i^2}}\right)$(4)

The weight assigned to each neighboring pixel, $w_{\mathit{ngb}}$, falls in the interval [0, 1] and is computed as the sum of the SAM values obtained for each band (Equation 5(5) $w_{\mathit{ngb}}=1-\frac{{\sum }_{j=1}^N\mathit{SA}{M}_j}{{\sum }_{i=1}^C{\sum }_{j=1}^N\mathit{SA}{M}_{\mathit{i},\mathit{j}}}$(5) ).

(5) $w_{\mathit{ngb}}=1-\frac{{\sum }_{j=1}^N\mathit{SA}{M}_j}{{\sum }_{i=1}^C{\sum }_{j=1}^N\mathit{SA}{M}_{\mathit{i},\mathit{j}}}$(5)

Where $\mathit{N}$ is the number of bands and $\mathit{C}$ is the number of neighboring pixels, i.e. eight.

Third, the column-wise aggregation of the $\mathit{n}\times \left(\mathit{T}-2\right)$ detail coefficient matrix in the thresholding step of the HiTS procedure is performed using two alternative strategies, depending on which one produces the highest aggregated detail coefficient. One approach is based on the mean of the detail coefficients relative to all the $\mathit{n}$ variables in the high-dimensional time series, where each band is weighted using the pixel-wise weights, computed as in Equation 5.(5) $w_{\mathit{ngb}}=1-\frac{{\sum }_{j=1}^N\mathit{SA}{M}_j}{{\sum }_{i=1}^C{\sum }_{j=1}^N\mathit{SA}{M}_{\mathit{i},\mathit{j}}}$(5) The alternative approach involves computing the mean of the detail coefficients considering only the bands in the focal pixel. HILANDYN can use these two approaches interchangeably, as the threshold for evaluating the significance of the deviation from linearity depends solely on the number of bands included in the high-dimensional time series (Equation 3(3) $\mathit{\lambda }=\mathit{C}\sqrt{2\mathit{log}\left(\mathit{nT}\right)},$(3) ). Finally, we applied the finite sample bias correction factor proposed by Park et al. (2019) to the MAD (Equation 2(2) ${\hat{\sigma }}_{\mathit{i}}=\left(\mathit{MAD}\left(X_{\mathit{i},1}-2X_{\mathit{i},2}+X_{\mathit{i},3},\dots ,X_{\mathit{i},\mathit{T}-2}-2X_{\mathit{i},\mathit{T}-1}+X_{\mathit{i},\mathit{T}}\right)\text{break}\times 1.4826\right)/\sqrt{6}$(2) ) to improve the estimation of the standard deviation.

2.6. Noise filter

HILANDYN typically ignores impulsive noise caused by artifacts in reflectance composites when it is sparse among the time series (Figure 4a) thanks to the modifications to the HiTS procedure. Nevertheless, changepoints caused by dense impulsive noise (Figure 4b) can be confounded with those associated with real changes such as forest disturbances (Figure 4c). Therefore, we designed a two-stage noise filter that is iteratively applied to remove spurious changepoints (Figure 2). At each iteration, the noise filter first identifies candidate changepoints, i.e. changepoints potentially caused by impulsive noise, and then processes them. The activation of the noise filter and the maximum number of iterations allowed for processing impulsive noise are controlled through the $\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$ parameter.

Figure 4. Example of time series of three focal pixels containing sparse noise (a), dense noise (b) or spectral changes caused by a wildfire (c). HILANDYN ignored sparse noise (a) while detected the dense ones in 1986 and 2002 (b). The wildfire occurred in 1991 (c) caused a changepoint characterised by a spike similar to impulsive noise though some bands, i.e. SWIR1, SWIR2 and TCW, did not exhibit a rapid spectral recovery.

During the initial stage (Figure 5), the noise filter searches for spurious changepoints by analyzing individual intervals, $t_1,\dots ,t_n$, which include one or multiple consecutive changepoints and the preceding time point. When $t_1$ is at the beginning of the time series, HILANDYN evaluates the reliability of changepoints in the interval through condition $C_1$(Equation 6(6) $t_1=1\wedge \left(\mathit{no}{b}_1<\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}\vee \mathit{no}{b}_2<\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}\right)$(6) , Figure 5).

Figure 5. Diagram of the operations performed by the noise filter during the screening stage for each interval. The number of observations per year and the parameter $\mathit{nob}\_{init}_{min}$ are required for evaluating condition $C_1$. Original and estimated spectral values are required for evaluating conditions $C_2$ and $C_3$. Time points at the beginning of the time series are first processed through condition $C_1$. Conditions $C_2$ and $C_3$ are checked on time points that exhibit the maximal spectral deviation from a reference value. Blue boxes indicate input data, black boxes indicate processing and green boxes indicate outputs.

(6) $t_1=1\wedge \left(\mathit{no}{b}_1<\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}\vee \mathit{no}{b}_2<\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}\right)$(6)

Where $\mathit{no}{b}_1$ and $\mathit{no}{b}_2$ are the number of observations available at $t_1$ and $t_2$, and $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ is a constant corresponding to the minimum number of clear observations required for producing reflectance composites. The number of clear observations available for each time point is computed as the median of the per-pixel clear observations within the spatial kernel. The value assigned to $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ should maximize the reliability of reflectance composites while not exceeding the maximum number of observations theoretically available for a given location and time period. Therefore, $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ primarily depends on the compositing approach, on the temporal frequency of Landsat data at the beginning of the time series, i.e. 8 or 16 days, and whether the study area is located where Landsat tiles overlap.

When condition $C_1$ is not satisfied, the noise filter finds a time point, $\mathit{ma}{x}_i$, in each band, whose values in the spatial kernel differ the most from the reference ones. For a single band, given a set of $\mathit{n}$ vectors, $\mathit{A}$ = $\left\{v_1,\dots ,v_n\right\}$, where $v_i\in {\mathbb{R}}^P$, with $\mathit{P}$ corresponding to the number of pixels in the spatial kernel, the noise filter finds the vector $v_{\mathit{max}}$ that maximizes the Euclidean distance to the reference vector $v_{\mathit{ref}}$ (Equation 7(7) $v_{\mathit{max}}=arg\mathit{ma}{x}_{\mathit{v}\in \mathit{A}}∥v_i-v_{\mathit{ref}}∥$(7) ).

(7) $v_{\mathit{max}}=arg\mathit{ma}{x}_{\mathit{v}\in \mathit{A}}∥v_i-v_{\mathit{ref}}∥$(7)

Where $∥.∥$ denotes the Euclidean norm and $v_{\mathit{ref}}$ contains the spectral values estimated by the HiTS procedure at the time point preceding $t_1$ or at that following $t_n$, depending on the position of $t_1$. The vectors in the set $\mathit{A}$ contain values fitted through the HiTS procedure except for $v_1$, which contains original input values.

A time point $\mathit{ma}{x}_i$ becomes a candidate time point, $\mathit{ch}{k}_i$, if at least one band at that location satisfies both conditions $C_2$ (Equation 8(8) $∥v_{\mathit{max}-1}-v_{\mathit{max}+1}∥<∥v_{\mathit{max}}-v_{\mathit{max}+1}∥$(8) ) and $C_3$ (Equation 9(9) $\sum _{\mathit{i}=1}^P\mathit{sgn}\left(v_{\mathit{i},\mathit{max}-1}-v_{\mathit{i},\mathit{max}}\right)\ne \mathit{sgn}\left(v_{\mathit{i},\mathit{max}}-v_{\mathit{i},\mathit{max}+1}\right)>0$(9) ).

(8) $∥v_{\mathit{max}-1}-v_{\mathit{max}+1}∥<∥v_{\mathit{max}}-v_{\mathit{max}+1}∥$(8)

(9) $\sum _{\mathit{i}=1}^P\mathit{sgn}\left(v_{\mathit{i},\mathit{max}-1}-v_{\mathit{i},\mathit{max}}\right)\ne \mathit{sgn}\left(v_{\mathit{i},\mathit{max}}-v_{\mathit{i},\mathit{max}+1}\right)>0$(9)

where $\mathit{v}\in {\mathbb{R}}^P$, with $\mathit{P}$ corresponding to the number of pixels within the spatial kernel, and $\mathit{sgn}$ is the sign function.

Condition $C_2$ (Equation 8(8) $∥v_{\mathit{max}-1}-v_{\mathit{max}+1}∥<∥v_{\mathit{max}}-v_{\mathit{max}+1}∥$(8) ) requires that the spectral distance between the time points preceding and following $\mathit{ma}{x}_i$ is lower than that between $\mathit{ma}{x}_i$ and the following time point. Condition $C_3$ (Equation 9(9) $\sum _{\mathit{i}=1}^P\mathit{sgn}\left(v_{\mathit{i},\mathit{max}-1}-v_{\mathit{i},\mathit{max}}\right)\ne \mathit{sgn}\left(v_{\mathit{i},\mathit{max}}-v_{\mathit{i},\mathit{max}+1}\right)>0$(9) ) requires that at least one pixel in the spatial kernel has to form a spike at $\mathit{ma}{x}_i$.

During the processing stage, the noise filter analyzes candidate time points individually, by iteratively reprocessing a subset of the high-dimensional time series (Figure 5). The algorithm removes all the data at $\mathit{ch}{k}_i$ and selects the variables relative to those bands that exhibited a significant change. Specifically, we deemed a changepoint in a specific band significant if the number of changed pixels within the spatial kernel is equal to or greater than the median number of changed pixels among all the bands. A changepoint is labeled as an actual change if a new changepoint is detected at the beginning of the interval. Spectral values within unreliable intervals are replaced with the average of the two observations following the interval, while time points associated with spurious changepoints are replaced through linear interpolation.

2.7. Accuracy assessment

Our definition of forest disturbance was as broad as possible and included changes associated with either stand-replacing, e.g. high-severity crown fires and windthrows, or non-stand-replacing events, e.g. drought-induced mortality and defoliating insect outbreaks. We considered individual events in space, i.e. at the scale of a Landsat pixel, and time, i.e. per year, following other studies (Bullock, Woodcock, and Holden 2020; Cohen et al. 2017). Therefore, we treated consecutive disturbances, e.g. windthrows followed by salvage logging, as separate events. Our validation dataset included 3000 plots (Figure 1, Table 2), each corresponding to a Landsat pixel, that were forested at any time during the analysis period. To locate validation plots, we used a stratified random sampling based on two strata delineated for each disturbed patch detected using an initial version of HILANDYN. Each stratum corresponded to either detected disturbed patches, i.e. the inner stratum, or a two-pixel positive buffer around these patches, i.e. the outer stratum. We classified disturbance events as stand-replacing or non-stand-replacing using the relative TCW change magnitude and a threshold of 50%, as proposed by Cohen et al. (2018). The number of stand-replacing disturbance events in the inner stratum was noticeably higher than that in the outer stratum (Table 2, Figure S2).

Table 2. Information on the validation dataset based on the validation strata. We discriminated between stand-replacing and non-stand-replacing events using a threshold of 50% of the relative TCW change magnitude.

Validation stratum	Total plots	Disturbed plots	Stand-replacing events	Non-stand-replacing events
Inner	1500	1439	1268	632
Outer	1500	706	408	494
Inner and outer	3000	2145	1676	1126

Using a buffer stratum was recommended for applications involving classifying rare phenomena such as forest disturbances (Olofsson et al. 2020). To discriminate between noise and disturbances, we considered only those events that occurred until the second-last year of the time series, i.e. until 2021.

Operationally, one operator determined the occurrence of a disturbance by simultaneously visualizing temporal profiles of annual Landsat time series at the plot-level and raster data in QGIS software, using the RasterDataPlotting and RasterTimeseriesManager plugins (https://raster-data-plotting.readthedocs.io/en/latest/index.html). Raster data comprised RGB false-color composites derived from Landsat yearly reflectance data (R = SWIR2, G = NIR, B = Red), nationwide aerial orthophotos, and high-resolution satellite imagery provided free-of-charge by Google, Bing, and Esri. By evaluating individual years for each plot, we obtained 111,000 sample units, of which 2802 were disturbance events as some plots were disturbed multiple times. Using a confusion matrix, we computed three accuracy metrics relative to the disturbed class: user’s accuracy (UA; Equation 10(10) $\mathit{UA}=\frac{\mathit{True}\mathit{Positive}}{\left(\mathit{True}\mathit{Positive}+\mathit{False}\mathit{Positive}\right)}\times 100$(10) ), i.e. precision, producer’s accuracy (PA; Equation 11(11) $\mathit{PA}=\frac{\mathit{True}\mathit{Positive}}{\left(\mathit{True}\mathit{Positive}+\mathit{False}\mathit{Negative}\right)}\times 100$(11) ), i.e. recall, and F₁ score (Equation 12(12) $F_1\mathit{score}=\frac{\mathit{UA}\times \mathit{PA}}{\mathit{UA}+\mathit{PA}}/100$(12) ).

(10) $\mathit{UA}=\frac{\mathit{True}\mathit{Positive}}{\left(\mathit{True}\mathit{Positive}+\mathit{False}\mathit{Positive}\right)}\times 100$(10)

(11) $\mathit{PA}=\frac{\mathit{True}\mathit{Positive}}{\left(\mathit{True}\mathit{Positive}+\mathit{False}\mathit{Negative}\right)}\times 100$(11)

(12) $F_1\mathit{score}=\frac{\mathit{UA}\times \mathit{PA}}{\mathit{UA}+\mathit{PA}}/100$(12)

2.8. Algorithm parametrization

HILANDYN requires a few parameters for processing high-dimensional Landsat time series. The primary ones are the bands to include in the time series and the constant factor $\mathit{C}$ (Equation 3(3) $\mathit{\lambda }=\mathit{C}\sqrt{2\mathit{log}\left(\mathit{nT}\right)},$(3) ). The constant factor $\mathit{C}$ and the number of bands influence the threshold λ (Figure S1; Section 2.5), which controls the sensitivity toward deviations from linearity in the HiTS procedure. Furthermore, the selection of bands is crucial as they are characterized by different sensitivities to forest disturbances (Cohen et al. 2018; DeVries et al. 2016; Schultz et al. 2016). We assessed how the parameters required by HILANDYN (Table 3) affected the accuracy metrics across the two validation strata (Section 2.7). Specifically, we first searched for the best performing combination of input bands (Table 1) and constant $\mathit{C}$ in terms of F₁ score and then tested the parameters for configuring the noise filter. The latter is controlled by the parameters $\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$ and $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ (Section 2.6).

Table 3. List of the main parameters required by HILANDYN and values tested for the sensitivity analysis.

Parameters	Tested alternatives or values
Input bands	All the unique combinations with length in the interval [1, 14], obtained using six original Landsat bands (Blue, Green, Red, NIR, SWIR1, SWIR2) and spectral indices listed in Table 1.
Constant $\mathit{C}$	Values in the interval [0.6, 1.2] increased by 0.2.
$\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$	Values in the interval [0, 4] increased by 1.
$\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$	Values in the interval [2, 10] increased by 1.

2.9. Comparison with BEAST

We selected BEAST (Zhao et al. 2019) to perform a comparison with HILANDYN based on our validation dataset. BEAST builds additive models by considering time series as linear combinations of seasonality, trend, abrupt changes, and noise, similar to BFAST (Breaks For Additive Seasonal and Trend; Verbesselt et al. 2010). Here, we employed the beast123 function provided in the Rbeast (version 0.9.9) R package to analyze the annual Landsat time series (Section 2.2) that included data from the spatial and spectral dimensions. Specifically, we tested time series that included data relative to one band and one pixel, one band and nine pixels, i.e. within a three-by-three spatial kernel, multiple bands and one pixel, and multiple bands and nine pixels. Input time series included either NBR due to its effectiveness in disturbance detection, e.g. Cohen et al. (2020), or those bands that achieved the best F₁ score in the sensitivity analysis of HILANDYN (Section 2.8).

Since our time series have annual frequency, we excluded the seasonal component, forcing BEAST to estimate linear trends, i.e. a polynomial order of one. As pointed out by Zhao et al. (2019), BEAST is sensitive to two hyperparameters: the maximum number of seasonal and trend changepoints and the minimum distance in time between two neighboring changepoints in a model. We therefore tested the maximum number of trend changepoints (trendMaxKnotNum) for values in the interval [1, 10], while we set 1 year as the minimum distance in time between two neighboring changepoints. This latter parameter allowed us to make a direct comparison with HILANDYN. We used the pre-defined values of the algorithm for all the other parameters.

Though BEAST can handle high-dimensional time series, being a univariate algorithm, it computes the probability distributions associated with the number of trend changepoints (ncpPr) and their position in time (cpOccPr) separately for each variable. Therefore, we used the mean to summarize these probability distributions among the variables in the high-dimensional time series. To find the position of the most probable changepoints, we followed the approach proposed by the authors of the algorithm. That is, we first obtained the mean number of changepoints, ncp, using the corresponding mean probability distribution. We then selected a set of changepoints, whose number was equal to ncp, among those with the highest mean occurrence probability. Finally, we used the same approach employed by HILANDYN (Section 2.3) for identifying trend changepoints associated with disturbances.

We compared BEAST and HILANDYN in terms of computation time by including from 1 to 14 bands in the high-dimensional time series within a three-by-three spatial kernel. Each combination of bands included those that performed best with HILANDYN in terms of producer’s accuracy when the constant factor $\mathit{C}$ was equal to one (Table S1).

3. Results

3.1. Quantitative assessment: detection capabilities and sensitivity analysis of the parameters

The constant factor $\mathit{C}$ and the number of bands included in the high-dimensional Landsat time series influenced the accuracy metrics noticeably (Figure 6, Table S1). Increasing the value of $\mathit{C}$ and the number of bands markedly gained UA while partly reduced PA, though this latter remained stable until the number of bands was equal or lower than seven. The highest increase in UA occurred when more than one band was included in the high-dimensional time series (Figure 6, Table S1). Without the noise filter, HILANDYN achieved the highest PA (93.3%) with the constant $\mathit{C}$ equal to 0.6 and seven bands in the high-dimensional time series, while it achieved the highest UA (85.5%) with $\mathit{C}$ equal to 1.2 and 12 band. The algorithm achieved the highest F₁ score (0.802) when we set $\mathit{C}$ equal to 1 and included seven bands (Figure 6, Table S1).

Figure 6. Maximum values of user’s accuracy (UA), producer’s accuracy (PA) and F₁ score based on the validation dataset and computed from all the unique combinations of bands with lengths from one to 14. The value of the constant $\mathit{C}$ ranged from 0.6 to 1.2. The noise filter was deactivated.

The selection of bands to include in the high-dimensional time series influenced primarily the PA, while the UA was mainly determined by their number (Figure 7). Considering the results obtained with the constant factor $\mathit{C}$ in the interval [0.6, 1.2] and all the unique combinations of bands, seven of them were included most often in the best-performing combinations in terms of PA (Figure 8). Specifically, MSI, NBR, SWIR2, TCW, SWIR1, TCA, and NDMI emerged as the most effective bands for forest disturbance detection (Figure 8). Differences with the remaining bands in terms of selection frequency were noticeable.

Figure 7. Distribution of values of user’s accuracy (UA), producer’s accuracy (PA) and F₁ score based on the validation dataset and computed from all the unique combinations of bands with lengths from one to 14. The value of the constant factor $\mathit{C}$ was equal to one. The noise filter was deactivated.

Figure 8. Frequency count relative to the number of times each band was in the best-performing combination in terms of PA, using from one to 14 bands in the high-dimensional time series. Results obtained with values of the constant factor $\mathit{C}$ ranging from 0.6 to 1.2 were pooled together. The noise filter was deactivated.

Using the configuration that produced the best performance in terms of F₁ score, i.e. the combination of SWIR1, SWIR2, NDMI, NBR, MSI, TCW, and TCA, and $\mathit{C}$ equal to one (Figure 6, Table S1), the noise filter significantly increased UA while marginally decreased PA (Figure S3, Table S2). Specifically, UA increased up to 9.7%, i.e. from 73.7% to 83.4%, while PA decreased up to 4.6%, i.e. from 87.8% to 83.2%. The noise filter increased the F₁ score slightly, from 0.801 to 0.833 when $\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$ was equal to four and $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ was equal to five (Figure S3, Table S2). HILANDYN achieved a substantial balance between UA and PA when we activated the noise filter, and this was consistent across the different validation strata (Table 4). The increase in UA caused by the noise filter was greater in the outer stratum than in the inner one, with a variation of 13% and 6.6%, respectively (Table 4). Similarly, the reduction in PA caused by the noise filter was greater in the outer stratum than in the inner one, with a reduction of 6.5% and 3.3%, respectively (Table 4). HILANDYN achieved comparable accuracies in validation plots dominated by either broadleaf or conifer tree species (Table S3).

Table 4. User’s accuracy (UA), producer’s accuracy (PA), and F₁ score achieved by HILANDYN with the noise filter activated or not. Accuracy metrics were computed using either both strata of the validation dataset or a single one. High-dimensional Landsat time series included data from seven bands (SWIR1, SWIR2, NDMI, NBR, MSI, TCW, and TCA) within a three-by-three spatial kernel. We parametrized the algorithm using $\mathit{C}$ = 1, $\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$ = 4, and $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ = 5.

Validation stratum	Noise unfiltered			Noise filtered
	UA (%)	PA (%)	F1 score	UA (%)	PA (%)	F1 score
Inner	80.9	92	0.861	87.5	88.7	0.881
Outer	60.5	78.9	0.685	73.5	72.4	0.730
Inner and outer	73.7	87.8	0.801	83.1	83.5	0.833

Overall, the noise filter allowed HILANDYN to reduce the influence of spurious changepoints caused by impulsive noise, as highlighted by the decrease from 23.5% to 14.9% in the median relative TCW change magnitude of commission errors (Figure 9, Table S4). For omission errors, the noise filter did not change the median relative TCW change magnitude, as this latter slightly decreased from 46.6% to 45.3% when we filtered noise (Figure 9, Table S4).

Figure 9. Distribution of the relative TCW change magnitude of commission (100 – UA) errors and omission (100 – PA) errors for HILANDYN and BEAST. High-dimensional Landsat time series included data of seven bands (SWIR1, SWIR2, NDMI, NBR, MSI, TCW and TCA) within a three-by-three spatial kernel. We parametrised each algorithm using the settings that provided the best results in terms of F₁ score. We constrained the relative TCW change magnitude between −50% and 250% by assigning these values to those outside that range.

3.2. Comparison with BEAST

The accuracy of BEAST, especially the PA and F₁ score, increased when we added both spectral and spatial data in the input time series (Table 5).

Table 5. User’s accuracy (UA), producer’s accuracy (PA), and F₁ score achieved by BEAST with the validation dataset as a function of the maximum number of trend changepoints (trendMaxKnotNum) and the data in the time series. Input time series included either data of one band (NBR), or seven bands (SWIR1, SWIR2, NDMI, NBR, MSI, TCW, and TCA) of one or nine pixels, i.e. within a three-by-three spatial kernel. The highest value of each accuracy metric for every input data type is highlighted in bold.

Maximum number of trend changepoints (trendMaxKnotNum)	One band (one pixel)			One band (nine pixels)			Seven bands (one pixel)			Seven bands (nine pixels)
	UA (%)	PA (%)	F1 score	UA (%)	PA (%)	F1 score	UA (%)	PA (%)	F1 score	UA (%)	PA (%)	F1 score
1	65.1	53.0	0.584	70.7	57.4	0.634	69.3	58.8	0.636	73.0	62.5	0.674
2	65.0	52.9	0.584	70.9	57.6	0.635	69.0	58.4	0.633	73.1	62.7	0.675
3	58.9	64.1	0.614	64.3	73.2	0.684	62.7	71.9	0.669	67.2	76.9	0.717
4	50.2	68.4	0.579	54.8	78.6	0.646	53.6	77.5	0.633	55.3	82.5	0.662
5	46.0	70.1	0.555	51.1	80.4	0.625	49.4	79.0	0.608	52.7	84.4	0.648
6	44.0	70.7	0.542	48.9	81.2	0.611	47.6	79.4	0.595	50.7	85.3	0.636
7	43.2	70.9	0.537	47.9	81.6	0.604	46.8	79.4	0.589	49.9	85.5	0.630
8	43.1	71.0	0.536	48.1	82.1	0.607	46.8	80.0	0.590	49.6	85.6	0.628
9	43.0	71.2	0.536	48.0	81.7	0.605	46.7	79.7	0.589	49.5	85.8	0.628
10	42.8	70.8	0.534	47.6	81.7	0.601	46.5	79.7	0.587	49.7	85.7	0.629

In particular, the PA peaked at 85.8% and the F₁ score at 0.717, when we used high-dimensional time series that included data of the seven best-performing bands in terms of PA (SWIR1, SWIR2, NDMI, NBR, MSI, TCW, and TCA; Figure 8) extracted using a three-by-three spatial kernel. The maximum number of trend changepoints markedly influenced the accuracy of BEAST as changes in its values between 1 and 10 caused significant variations of UA and PA up to approximately 23% (Table 5). In particular, the PA increased when we raised the maximum number of trend changepoints, while the UA followed the opposite trend. Therefore, we achieved the highest F₁ score by setting the maximum number of trend changepoints equal to three (Table 5).

Overall, the best F₁ score provided by BEAST (0.717) was lower than that provided by HILANDYN (0.801) when the noise filter was deactivated. With regard to the best F₁ score, the UA achieved by BEAST (67.2%, Table 5) was lower than that of HILANDYN (73.7%, Table S2). Similarly, differences in PA between the algorithms were noticeable, with BEAST and HILANDYN achieving 76.9% (Table 5) and 87.8% (Table S2), respectively.

The distribution of commission, i.e. 100 – UA, and omission, i.e. 100 – PA, errors in terms of relative TCW change magnitude highlighted slight differences between HILANDYN and BEAST (Figure 9, Table S4). The median relative TCW change magnitude of commission errors was 23.5% for HILANDYN (noise filter deactivated) and 25.1% for BEAST. Conversely, for omission errors, these values were 46.6% for HILANDYN (noise filter deactivated) and 37.6% for BEAST.

BEAST required considerably more time to process Landsat time series of the validation dataset compared to HILANDYN (Figure 10, Table S5).

Figure 10. Execution time required by HILANDYN and BEAST to process the 39-year high-dimensional time series of the validation dataset. The total number of variables in each time series was equal to the product of the number of bands and nine, i.e. the number of pixels within a three-by-three spatial kernel. The bands selected in each combination were those that gave the best results in terms of PA with HILANDYN when the constant $\mathit{C}$ was equal to one (Table S1). We configured the noise filter of HILANDYN by setting the $\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$ to four and $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ to five. For BEAST, we used the predefined prior parameters except for trendMaxknotnum and trendMinsepdist that we set to eight and one, respectively. Each process ran on a single core of an AMD Ryzen 9 5950X processor.

In particular, differences between the two algorithms increased with the number of bands in high-dimensional time series, as BEAST required proportionally higher computational time compared to HILANDYN. For instance, processing high-dimensional time series with 63 variables, i.e. the values of seven bands extracted from nine connected pixels, required 40 and 1117 seconds for HILANDYN and BEAST, respectively (Figure 10, Table S5). When we activated the noise filter, the computation time required by HILANDYN increased to 64 seconds.

3.3. Qualitative assessment of the detection capabilities

HILANDYN was sensitive to a wide range of disturbance severities as it detected forest patches disturbed by stand-replacing events, such as wildfires, windthrow, and sleet damage, and non-stand-replacing events, such as tree crown dieback caused by defoliating insects and severe drought and heat waves (Figure 11). The spatio-temporal pattern of disturbed forest patches detected by HILANDYN exhibited a high degree of contiguity among pixels, which was consistent across different disturbance agents and severities (Figure 12). Moreover, HILANDYN effectively detected forest patches that were repeatedly disturbed, either a few years apart or in close succession, i.e. one year apart (Figure 13). Secondary disturbances occurring shortly after the main event were often associated with salvage logging in our study area. In particular, the temporal profiles of NBR and TCW highlighted different sensitivities to the complete removal of tree cover associated with salvage logging when it occurred shortly after the natural disturbance (Figure 13h–p).

Figure 11. Examples of stand-replacing and non-stand-replacing forest disturbances detected by HILANDYN in the study area between 1985 and 2022. Panels A-F depict the following disturbance events: (a) dieback induced by defoliating insect outbreak (Asian chestnut gall wasp, Dryocosmus kuriphilus) in 2011; (b) dieback caused by a severe drought and heatwave in 2003; (c) wildfires in 1990, 1991, 2003, and 2017; (d) windthrow in 2007; (e) ice storm in 2014; (f) windthrow in 2018. The grey background in the central panel indicates the presence of forest cover either in 1990 or 2018 according to the Corine Land Cover within the European Alps borders (blue line).

Figure 12. Examples of disturbed forest patches detected by HILANDYN and characterised by different disturbance agents and severities. Please refer to the legend in Figure 11 for the colour corresponding to each year in panels a-e. Clearcuts (a, f; 5°45‘18,102“E; 44°3’9,791“N), debris flows (b, g; 8°1‘54,389“E; 46°46’4,057“N), defoliating insect outbreak (c, h; 7°35‘10,709“E; 45°22’18,178“N), wildfire (d, i; 13°16‘54,707“E; 46°28’25,202“N), windthrow (e, j; 7°21‘17,61“E; 44°39’31,641“N). The grey background indicates the presence of forest cover either in 1990 or 2018 according to the Corine Land Cover.

Figure 13. Examples of forest patches, depicted in yellow, that were disturbed multiple times during the analysis period. Each row shows results from HILANDYN, recent high-resolution satellite imagery, and temporal profiles of NBR and TCW. We parametrised the algorithm according to the configuration that achieved the highest F₁ score (Section 3.1). Recurrent wildfires (a-d; 7° 9’1.02“E; 45° 9’35.57“N). Wildfire followed by salvage logging (e-h; 9°39’23.26“E; 46°11’59.51“N). Bark beetle attack followed by salvage logging (i-l; 14°42’0.63“E; 46°16’1.43“N). Windthrow followed by salvage logging (m-p; 11°21’13.93“E; 46°17’41.66“N).

4. Discussion

4.1. Key features of HILANDYN

Here, we presented a novel, unsupervised algorithm for retrospectively mapping forest disturbance dynamics, although its capabilities are not limited to any specific land cover. In the context of Landsat time series-based algorithms, the novelty introduced by HILANDYN lies in its ability to temporally segment high-dimensional Landsat time series containing multispectral information from adjacent pixels. By integrating information from the spectral and spatial domains into the time series, HILANDYN was able to achieve satisfactory accuracy (Table 4) while providing spatio-temporal consistency of disturbed forest patches (Figures 12, 13). These results are made possible by the modified HiTS procedure, which provides a computationally sound approach to reduce the dimensionality of high-dimensional time series in the early stages of the analysis. The construction of a unique, data-driven wavelet basis from all the variables in the time series (Section 2.4) and the aggregation strategy of the detail coefficients (Section 2.5) are key elements for detecting changepoints by exploiting the spectral and spatial domains. Widely used algorithms that can harness information from multiple bands, such as CCDC (Zhu and Woodcock 2014) and COLD (Zhu et al. 2019), primarily work in parallel, processing data in a band-wise manner and combining univariate tests to detect changepoints in later stages of the analysis. This results in high computational requirements and the need to run the algorithms on geospatial cloud computing platforms, such as Google Earth Engine, to analyze large areas (Bullock, Woodcock, and Holden 2020; Pasquarella et al. 2022).

Algorithms designed to exploit the spatial correlation between adjacent pixels via OBIA, e.g. VeRDET (Hughes, Kaylor, and Hayes 2017) and OB-COLD (Ye, Zhu, and Cao 2023), typically perform temporal segmentation either after or before the spatial analysis. In addition to increasing computational and algorithmic complexity, sequential approaches require appropriate parametrization for each stage of the analysis, such as setting the optimal scale for spatial segmentation. Conversely, the inclusion of spectral information from the spatial context in high-dimensional time series enhanced the capabilities of HILANDYN by extending the spatial scale of the analysis without increasing the number of parameters. Indeed, the pixel-wise weighting scheme and the thresholding approach in the modified HiTS procedure (Section 2.5) allowed the algorithm to take advantage of the spatial context while maintaining its sensitivity to changes occurring at the individual Landsat pixel.

The bottom-up, i.e. agglomerative, approach employed by the HiTS procedure to construct the adaptive wavelet basis has proved to be effective for the temporal segmentation of high-dimensional time series, as demonstrated by the accuracy metrics achieved by HILANDYN (Table 4). When applied to the segmentation of different types of signals, bottom-up approaches outperform those based on sliding windows and top-down algorithms (Keogh et al. 2004). They have also been successfully used for temporal segmentation in other Landsat time series-based algorithms, such as in C2C (Hermosilla et al. 2015). The HiTGUW transformation, i.e. the first stage of the HiTS procedure, focuses on local features in its early stages before identifying global features, which enables HILANDYN to perform well in detecting abrupt local changes, including point anomalies, e.g. Figure 4, and consecutive disturbances occurring in close succession, e.g. Figure 13. Bottom-up temporal segmentation also proved to be useful in the first stage of the noise filter built into HILANDYN, i.e. to find changepoints potentially caused by spectral noise. Indeed, when the noise filter was activated, the UA increased significantly (Table 4) and the change magnitude associated with commission errors decreased (Figure 9).

4.2. Algorithm parametrization

HILANDYN proved to scale well to a considerably large study area when configured with a single set of parameters that maximized the F₁ score (Table 4, Figures 11, 12). Indeed, we achieved an overall balance between UA and PA while detecting forest patches disturbed by different agents and a wide range of severity. Landsat time series-based algorithms typically require to fine-tune the main parameters to adapt to different study areas (Pasquarella et al. 2022) and specific disturbance agents (Ye et al. 2021). HILANDYN mostly missed non-stand-replacing events when using the best performing set of parameters, as shown by the relative TCW change magnitude for omission errors (Figure 9, Table S4). This behavior is common among Landsat time series-based automated algorithms, as the detection of low-severity disturbances is typically hindered by residual spectral noise in Landsat time series (Cohen et al. 2017; Rodman et al. 2021). Nevertheless, the high PA achieved by HILANDYN in the outer validation stratum, i.e. 78.9% with the noise filter deactivated (Table 4), evidenced its overall sensitivity to non-stand-replacing events.

Since the threshold $\mathit{\lambda }$ used by the modified HiTS procedure depends directly on the number of bands and the constant $\mathit{C}$ (Equation 3(3) $\mathit{\lambda }=\mathit{C}\sqrt{2\mathit{log}\left(\mathit{nT}\right)},$(3) ), the optimization of these parameters is critical to discriminate between low-severity disturbances and spurious changepoints. In fact, increasing the value of $\mathit{\lambda }$ reduces the sensitivity of HILANDYN to low-severity disturbances, which are typically characterized by small and transient deviations from linearity in the linear trends. It is worth noting that the optimization of $\mathit{\lambda }$ in HILANDYN allows the user to maximize the balance between UA and PA (Figure 6, Table S1), thus reducing the need to post-process the disturbance metrics generated by the algorithm. In contrast, other temporal segmentation algorithms such as LandTrendr, C2C, and BEAST require setting a predefined maximum number of changepoints or linear segments.

When we included the best performing bands in the high-dimensional time series, their number and the constant $\mathit{C}$ had a significant effect on UA, whereas PA was less sensitive to their variations (Figure 6, Table S1). As observed with other Landsat time series-based algorithms, such as LandTrendr (Cohen et al. 2018, 2020; Pasquarella et al. 2022) and COLD (Cohen et al. 2020), the use of a single band introduces a significant amount of uncertainty in temporal segmentation, primarily due to increased commission errors. The heterogeneity of spectral changes across wavelengths associated with disturbance events (Pasquarella et al. 2022), and the difficulty of effectively distinguishing between low-severity disturbances and noise are the main sources of uncertainty. Secondary classification models that combine the results of multispectral ensembles obtained from univariate algorithms, e.g. LandTrendr, were primarily aimed to increase UA of disturbance maps (Cohen et al. 2018, 2020; Senf and Seidl 2020). HILANDYN benefited from the increase in the number of bands in terms of UA (Figure 6; Figure 7) because this likely reduced the influence of band-dependent noise, as shown in Figure 4a. The latter can be caused by several factors, such as the presence of cirrus clouds, which alter the surface reflectance unevenly across wavelengths, with a smaller effect on the SWIR bands than on the others (Qiu, Zhu, and Woodcock 2020). Conversely, forest disturbances and band-independent noise typically cause dense changes in high-dimensional time series (Figure 4b–c) and further processing is required to correctly separate them (Section 2.6).

Regarding the selection of bands, their effectiveness in detecting forest disturbances largely depends on the residual information they carry after the normalization by the estimated standard deviation of the random error ($\hat{\sigma };$Equation 3(3) $\mathit{\lambda }=\mathit{C}\sqrt{2\mathit{log}\left(\mathit{nT}\right)},$(3) ). In this study, the best-performing bands in terms of PA (Figure 8) were seven and, except for TCA, they were based on the SWIR wavelengths, either the original reflectance bands, i.e. the SWIR2 and SWIR1 bands of the Landsat sensors, or the spectral indices contrasting the NIR and the SWIR bands (MSI, NBR, TCW, and NDMI). The effectiveness of the SWIR bands for forest disturbance detection has been reported for different forest ecosystems and different change detection algorithms (Cohen et al. 2018; DeVries et al. 2016; Hislop et al. 2019). In particular, algorithms that use one image per year appear to rely mostly on SWIR bands to detect forest disturbance (Cohen et al. 2018, 2020). Conversely, algorithms that fit harmonic functions to all the available Landsat images, such as COLD, rely primarily on information in the NIR band due to its effectiveness in reconstructing forest phenology (Cohen et al. 2020). Results in Figure 8 also highlighted the usefulness of the TCA band, as it complemented the information provided by the SWIR-based bands. Consistent with previous studies (Cohen et al. 2018, 2020), our results showed that TCA outperformed the widely used NDVI in detecting forest disturbances when using one image per year due to its greater sensitivity at high levels of vegetation cover (Cohen et al. 2016).

Optimizing the two parameters required by the integrated noise filter, i.e. $\mathit{noise}\_\mathit{ite}{r}_{\mathit{max}}$ and $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$, further improved the balance between UA and PA (Table S2, Figure S4), as HILANDYN effectively discriminated between noise and trustworthy changes (Table 4, Figure 9). The $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ parameter is aimed to reduce over-segmentation at the beginning of the time series, which is caused by residual noise in the reflectance composites due to data scarcity. Indeed, $\mathit{nob}\_\mathit{ini}{t}_{\mathit{min}}$ strictly depends on the compositing algorithm, as this latter can significantly affect the quality of the reflectance composites (Qiu et al. 2023).

4.3. Comparison with BEAST

Despite being a univariate algorithm, BEAST benefited from the use of multispectral information and from the inclusion of spatial context information in the time series, achieving the best results when we combined this information (Table 5). In fact, increasing the number of variables in high-dimensional time series increased both UA and PA to values similar to those achieved by HILANDYN (Section 3.1). Nevertheless, BEAST struggled to achieve a balance between UA and PA, as evidenced by the relatively low F₁ score (Table 5). This was determined by the strong dependence of UA and PA on the maximum number of trend changepoints, which are inversely correlated with this parameter. In particular, the relatively low UA achieved by BEAST, i.e. around 55% or less, when we set the maximum number of trend changepoints higher than three (Table 5), seems to confirm the larger proportion of commission over omission errors initially reported by Zhao et al. (2019). The analysis of the relative TCW change magnitude associated with commission and omission errors showed that BEAST and HILANDYN had comparable sensitivity to disturbances in terms of severity (Figure 9, Table S4). Nevertheless, the relative TCW change magnitude of omission errors evidenced that BEAST was slightly more sensitive to non-stand-replacing events compared to HILANDYN (Figure 9, Table S4). This highlighted the effectiveness of the Bayesian approach in detecting subtle trend changepoints in Landsat time series. Overall, the main limitations associated with this approach stem primarily from the significant computing power required, which increases linearly with the number of variables included in the time series (Figure 10). This can be a significant challenge for applications over large areas or when using high spatial resolution data, unless high-performance computing platforms are used.

4.4. Limitations and future developments

Compared to other Landsat time series-based algorithms, such as BFAST Monitor (DeVries et al. 2015) and COLD (Zhu et al. 2019), we developed HILANDYN to analyze forest disturbance dynamics retrospectively. Therefore, the algorithm lacks crucial features dedicated to real-time monitoring . First, the HiTS procedure is a complete offline approach; thus, a reanalysis of the entire time series is required when new data is added. A possible solution could come from online approaches such as SWAB (Keogh et al. 2001), which combine sliding windows and bottom-up methods by applying the latter within sliding windows of constant size $\mathit{w}$. Notably, a sliding window approach has been adopted by Mulverhill et al. (2023) to add monitoring capabilities to BEAST. We note that online algorithms, such as CCDC and COLD, need to reinitialize the time series model after a disturbance event, which may result in the omission of consecutive disturbances occurring in close succession. For this reason, backward fitting of the data has been proposed after a disturbance event for the OB-COLD algorithm (Ye, Zhu, and Cao 2023).

Second, the annual frequency of the time series employed by HILANDYN hinders the timely detection of forest disturbances, which is a fundamental requirement for monitoring purposes. While using time series with evenly spaced observations is required by the HiTS procedure, seasonal reflectance composites could increase the number of observations available during the year, as proposed by some studies based on LandTrendr (Meng et al. 2021). Similarly, approaches for normalizing phenological variations offer the opportunity to model intra-annual time series through linear trends (Yang et al. 2023). An increase in the temporal frequency of high-dimensional time series from inter- to intra-annual would also limit the reduction in PA caused by the noise filter (Table 4). A higher temporal frequency would reduce the omission of low-severity disturbances caused by the masking effect of post-disturbance vegetation recovery (Yang et al. 2023). However, Giannetti et al. (2021) reported that the optimal time interval required to achieve the highest accuracy of disturbance maps, regardless of the algorithm, was between 7 and 12 months.

Third, linear trends may be suboptimal for modeling satellite time series as nonlinear trends may be better suited to represent forest dynamics (Moisen et al. 2016; Senf, Müller, and Seidl 2019). A possible improvement to HILANDYN could be the extension to piecewise-quadratic signals, as proposed by Maeng and Fryzlewicz (2023) for the TrendSegment procedure, i.e. the univariate version of the HiTS procedure. Nevertheless, changes in linear trends have been found to be useful for detecting non-stand-replacing disturbances and forest dieback in general (Moreno-Fernández et al. 2021; Mulverhill, Coops, and Achim 2023).

5. Conclusions

As changes in disturbance regimes are critical to understand the effects of climate change on forest ecosystems, automated and accurate methods are needed to reconstruct regional-scale forest disturbance dynamics from long-term satellite time series. HILANDYN offers new opportunities to study forest disturbance dynamics retrospectively by exploiting the temporal, spatial, and spectral domains provided by Landsat data. In fact, temporal segmentation of high-dimensional time series through the modified HiTS procedure is the key element of our algorithm. It ensures an accurate detection of disturbance events and guarantees computational efficiency by reducing the dimensionality of the time series in the early stages of the processing by locating deviations from linearity that are common between variables. The comparison with BEAST confirms the robustness of HILANDYN for the detection of stand-replacing and non-stand-replacing disturbances. Its scalability to large study areas is supported by the high computational efficiency and flexibility to heterogeneous conditions when using a single set of parameters. The methods implemented in HILANDYN are general enough to analyze time series of optical satellite imagery other than Landsat, such as that acquired by the Sentinel-2 satellites.

Acknowledgments

We acknowledge the U.S. Geological Survey for distributing Landsat imagery and Dr Andreas Rabe for providing the Raster Timeseries Manager plugin for QGIS.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The code of the algorithm is implemented in the R package “hilandyn” publicly available via the GitHub repository of the first author at https://github.com/donatomorresi/hilandyn.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2024.2365001

Additional information

Funding

This work was supported by the European Union through the MOSAIC Interreg Alpine Space project (Managing prOtective foreSt fAcIng clImate Change compound events) under Grant [ASP0100014].