Landsat Imagery Spectral Trajectories—important Variables for Spatially Predicting the Risks of Bark Beetle Disturbance

Tree mortality caused by bark beetle infestation has significant effects on the ecology and value of both natural and commercial forests. Therefore, prediction of bark beetle infestations is critical in forest management. Existing predictive models, however, rarely consider the influence of long-term stressors on forest susceptibility to bark beetle infestation. In this study we introduce pre-disturbance spectral trajectories from Landsat Thematic Mapper (TM) imagery as an indicator of long-term stress into models of bark beetle infestation together with commonly used environmental predictors. Observations for this study come from forests in the central part of the Šumava Mountains, in the border region between the Czech Republic and Germany, Central Europe. The areas of bark beetle-infested forest were delineated from aerial photographs taken in 1991 and in every year from 1994 to 2000. The environmental predictors represent forest stand attributes (e.g., tree density and distance to the infested forest from previous year) and common abiotic factors, such as topography, climate, geology, and soil. Pre-disturbance spectral trajectories were defined by the linear regression slope of Tasseled Cap components (Wetness, Brightness and Greenness) calculated from a time series of 16 Landsat TM images across years from 1984 until one year before the bark beetle infestation. Using logistic regression and multimodel inference, we calculated predictive models separately for each single year from 1994 to 2000 to account for a possible shift in importance of individual predictors during disturbance. Inclusion of two pre-disturbance spectral trajectories (Wetness slope and Brightness slope) significantly improved predictive ability of bark beetle infestation models. Wetness slope had the greatest predictive power, even relative to environmental predictors, and was relatively stable in its power over the years. Brightness slope improved the model only in the middle of the disturbance period (1996). Importantly, these pre-disturbance predictors were not correlated with other predictors, and therefore bring additional explanatory power to the model. Generally, the predictive power of most fitted model decreases as time progresses and models describing the initial phase of bark beetle outbreaks appear more reliable for conducting near-future predictions. The pre-disturbance spectral trajectories are valuable not only for assessing the risk of bark beetle infestation, but also for detection of long-term gradual changes even in non-forest ecosystems.


Introduction
Forest disturbances caused by wind, fire and insects are prominent forces controlling forest ecosystem dynamics across the world [1].Insect forest disturbances represent slower events relative to wind and fire, but have strong effects because of irreversible defoliation and subsequent dieback of mature trees [2].Bark beetles are among the most common insect disturbance agents in temperate coniferous forests, affecting tens of millions of hectares in recent decades [3,4].Moreover, the affected forest area is predicted to double in the next two decades because of climate change and cultivation of large areas of coniferous plantations outside their natural distribution [5].Even though bark beetle disturbances are a part of the natural dynamics of coniferous forest [6][7][8] they also represent a significant threat for timber production and certain forest-provided ecosystem services [9].Therefore, the need for both detection and prediction of bark beetle infestation is increasing.
Bark beetle outbreaks are commonly detected from satellite imagery using spectral indices, which are primarily based on differences and/or ratios between multiple spectral bands [10].While the first studies of bark beetle infestation used changes in red and near-infrared (NIR) bands (i.e., normalized difference vegetation index (NDVI)), indices based on NIR and shortwave infrared (SWIR) have become the standard, mostly because of their greater sensitivity to forest canopy defoliation caused by bark beetle infestation [11].The widely used Tasseled Cap (TC) linear transformation has three components: Brightness, Greenness, and Wetness [12,13].Among these three components, Wetness is often used to identify bark beetle forest disturbances.It is a combination of all Landsat TM or Enhanced Thematic Mapper Plus (ETM+) spectral bands except the thermal band, with the greatest weight given to the SWIR channel (e.g., [14]).
While the use of remote sensing for detection of bark beetle infestation is well established, its application for prediction or risk assessment is rather rare [15][16][17].Empirical bark beetle risk-assessment models commonly use environmental predictors including topography, climate, or soil types [18][19][20].While these environmental predictors are important for bark beetle population dynamics and tree vitality, the forest stand attributes like structure, composition and spatial configuration have at least the same importance; however, they differ at various phases of infestation [16,21,22].
Information on tree vigor, which is the crucial factor shaping the risk of infestation, is usually missing in the standard forest inventory data.Earth observation imagery, specifically from the Landsat sensor family, can fill this lack of information by providing data on tree vitality and forest health status [23].When analyzed as a time series, Landsat can provide information regarding pre-disturbance change events (e.g., windstorm) or gradual changes in forest vitality (e.g., soil depletion or air pollution), which increase the stand's vulnerability to bark beetle infestation [24][25][26][27].Besides the tree to stand-level factors, the risk of bark beetle infestation is increased by large-scale drivers, which include climatic variables and/or landscape connectivity [28].Most of these latter changes are well reflected in spectral indices and can be documented using multi-temporal remote sensing [29].Free access to the Landsat data archive by the US Geological Survey (USGS) in 2008 accelerated the use of these data for multi-temporal change detection [30].Subsequently, temporal segmentation algorithms and spectral trajectory processing tools like LandTrendr [31] and TimeSync [32] were developed to extract the main features (e.g., greatest disturbance segment, segment duration, and magnitude) from spectral time series trajectories.These approaches have enabled a wide use of Landsat data for broad scale ecosystem studies [33].The use of spectral trajectories based on Landsat imagery helped to describe forest disturbance history since 1972 [34], encouraging forest recovery [35], and detecting bark beetle-affected forest stands more precisely [36].The following studies distinguished different defoliator and bark beetle effects [37], duration and severity of insect disturbances [38], and large-scale drivers of primary forest disturbance events [39].A similar trajectory-based approach could be used to detect less pronounced pre-disturbance forest state changes which predispose the forest to disturbances [40].
We hypothesize that gradual changes in forest vitality significantly increase forest vulnerability to bark beetle infestation.Therefore our primary objective is to determine if pre-disturbance spectral trajectories are useful predictors of bark beetle infestation.In particular, we examine how explanatory and predictive power of models that include pre-disturbance spectral trajectories change relative to models in which these variables are excluded.For this purpose, we combine pre-disturbance predictors with environmental predictors in a logistic regression model combined with multimodel inference using several sets of models to assess the importance of pre-disturbance spectral trajectories.

Study Area
The study area (1376 ha, centroid N 48 • 57.78435 , E 13 • 28.33847 ) is located in the central part of the Šumava Mountains, in the border region between the Czech Republic and Germany, Central Europe (Figure 1).It is part of the Šumava National Park (Czech Republic) and is adjacent to the Bavarian Forest National Park (Germany).The dominant natural vegetation unit is Norway spruce forest (Picea abies L. Karst.)[41,42].Where the area is not covered by forest, mostly peats and mountain meadows occur.The relief of the study area is created by a peneplain with numerous local hills.The elevation ranges from 700 m to 1453 m with an average elevation of 1093 m. distinguished different defoliator and bark beetle effects [37], duration and severity of insect disturbances [38], and large-scale drivers of primary forest disturbance events [39].A similar trajectory-based approach could be used to detect less pronounced pre-disturbance forest state changes which predispose the forest to disturbances [40].
We hypothesize that gradual changes in forest vitality significantly increase forest vulnerability to bark beetle infestation.Therefore our primary objective is to determine if pre-disturbance spectral trajectories are useful predictors of bark beetle infestat ion.In particular, we examine how explanatory and predictive power of models that include pre-disturbance spectral trajectories change relative to models in which these variables are excluded.For this purpose, we combine pre-disturbance predictors with environmental predictors in a logistic regression model combined with multimodel inference using several sets of models to assess the importance of pre-disturbance spectral trajectories.

Study Area
The study area (1376 ha, centroid N 48°57.78435′,E 13°28.33847′) is located in the central part of the Šumava Mountains, in the border region between the Czech Republic and Germany, Central Europe (Figure 1).It is part of the Šumava National Park (Czech Republic) and is adjacent to the Bavarian Forest National Park (Germany).The dominant natural vegetation unit is Norway spruce forest (Picea abies L. Karst.)[41,42].Where the area is not covered by forest, mostly peats and mountain meadows occur.The relief of the study area is created by a peneplain with numerous local hills.The elevation ranges from 700 m to 1453 m with an average elevation of 1093 m.Over the last 25 years, several outbreak waves of the European spruce bark beetle (Ips typographus L.) have caused extensive spruce forest dieback in this area.The outbreak originated in wind-fallen trees (173 ha) in 1983 and 1984, following two windstorms in the western part of the Bavarian Forest National Park.These fallen and heavily damaged trees provided suitable conditions for the start of the outbreak [43,44].At the beginning of the 1990s the first infested trees occurred in the Czech region.Because of the subsequent bark beetle dispersal, the Šumava NP administration applied two management approaches to the infested forests in 1995: (1) a small portion of the stands in the core zone was left unmanaged (more than 2.000 ha); and (2) the stands in the surrounding buffer zone were cut to prevent the spread of bark beetles.We analyzed the area in the unmanaged zone where we expect a dominant role of natural predictors of tree infestation.The bark beetle infestation accelerated between 1995 and 1997 when the largest area in the central part of Šumava Mountains was affected [45].By the end of the year 2000, most of the spruce stands lost all canopy trees and only small patches of living trees persisted, mostly near peat bogs or at forest edges.After the main dieback in the study area the forest regenerated successfully towards a natural mountain spruce forest [46,47].

Aerial and Satellite Data
Groups of trees infested by bark beetle (grey-attack stage) were delineated from ortho-rectified aerial photographs, which were acquired in 1991 and every year for the period 1994-2000.The smallest polygon was defined as a group of five adjacent dead trees (0.03 ha).Simultaneously, living spruce forest, clear-cuts and non-forest areas (meadows, peat-bogs) were also delineated.To avoid their influence on values of mixed pixels, we excluded all pixels up to 100 m to clear-cuts and non-forest areas.The resulting land cover polygon layers were transformed to UTM Zone 33 North coordinate system and overlaid on a 30 m × 30 m grid to make them compatible with Landsat imagery.Pixels were defined as declined when the dieback was higher than 50% of its area.We used 16 Landsat TM images from 1984 to 1999 (Table 1).For some years (1988,1994), multiple images were used to cover no data areas because of clouds and their shadows.To avoid phenological differences among the images, only data acquired during the high summer season were used (nominally July and August with an exception in 1999).The imagery was obtained from the USGS data archive (glovis.usgs.gov).All images were radiometrically normalized to a common image (acquired on 4 August 1990) using the MADCAL approach of Canty et al. [48], following Schroeder et al. [49].

Spectral Trajectories
Based on the study of Hais et al. [40], which shows the potential of spectral trajectories to describe forest vigor decline before bark beetle infestation, we used this approach to test the significance of these predictors for modelling bark beetle infestation risk.The abovementioned approach uses linear regression to fit multi-temporal spectral data.In our study, we consider all three TC components (Wetness, Brightness and Greenness), though mainly Wetness and Brightness have been shown as the appropriate multi-spectral vegetation indices to describe forest structure [50] and forest health in relation to insect disturbance [14,[51][52][53].To obtain the pre-disturbance spectral trajectories, we calculated the linear regression slope of each TC component from 1984 to one year before the year for which the model was calibrated.The slope was calculated for each 30 m × 30 m pixel except those already infested in previous years.Slopes for non-affected stands were calculated from the same years as the affected ones in the one-year model.An example of Wetness slope is shown in Figure 2.These slopes, which represent historical stand condition (e.g., stress, stability, growth), were included as a model predictor together with commonly used environmental predictors in models of bark beetle infestation risk.
Based on the study of Hais et al. [40], which shows the potential of spectral trajectories to describe forest vigor decline before bark beetle infestation, we used this approach to test the significance of these predictors for modelling bark beetle infestation ris k.The abovementioned approach uses linear regression to fit multi-temporal spectral data.In our study, we consider all three TC components (Wetness, Brightness and Greenness), though mainly Wetness and Brightness have been shown as the appropriate multi-spectral vegetation indices to describe forest structure [50] and forest health in relation to insect disturbance [14,[51][52][53].To obtain the pre-disturbance spectral trajectories, we calculated the linear regression slope of each TC component from 1984 to one year before the year for which the model was calibrated.The slope was calculated for each 30 m × 30 m pixel except those already infested in previous years.Slopes for non -affected stands were calculated from the same years as the affected ones in the one-year model.An example of Wetness slope is shown in Figure 2.These slopes, which represent historical stand condition (e.g., stress, stability, growth), were included as a model predictor together with commonly used environmental predictors in models of bark beetle infestation risk.

Environmental Predictors of Bark Beetle Infestation
Despite a large amount of environmental predictors tested for the risk of bark beetle infestation [19,21], we selected representatives of each group (see Table 2) to cover all types of possible influences.Forest stand susceptibility to bark beetle infestation is determined by forest condition itself or abiotic environmental factors, which in turn influence forest vigor.One of the most significant forest attributes for bark beetle dispersal is distance to the nearest infested forest [21].Despite just 100 m buffer zones of forest in each year, it is still possible to test the distance to the nearest declined forest inside these buffer zones.In addition, we tested the direction to the nearest declined forest in the previous year.The direction azimuth value was transformed to eight 45° wide classes (N, NE, E, ES, S, SW, W and NW).Distance to the nearest clear -cut area and distance to the nearest natural forest edge were taken as a static variable from aerial photograph acquired in 1991.We hypothesize that the clear-cut edges might ease the entry of bark beetles into the stand while the natural forest edges represent a physical barrier for the beetles.Moreover, the trees in forest edge adjacent to clear -cut may suffer from water stress due to higher amount of solar radiation [20,54].The degree of forest "naturalness" was calculated as a difference between current forest inventory and potential natural tree species composition [55], derived from modified Czech typological maps (natural forest types [56]), which followed a similar concept as potential natural vegetation [57].Stand age represents an important attribute and we used age of dominant (>50%) tree layer because bark beetles preferentially attack trees older than 60 years [44].If any tree layer does not reach 50%, we calculated age as a mean value from the first and second dominant tree layer age.Similarly, the dominant tree layer was used

Environmental Predictors of Bark Beetle Infestation
Despite a large amount of environmental predictors tested for the risk of bark beetle infestation [19,21], we selected representatives of each group (see Table 2) to cover all types of possible influences.Forest stand susceptibility to bark beetle infestation is determined by forest condition itself or abiotic environmental factors, which in turn influence forest vigor.One of the most significant forest attributes for bark beetle dispersal is distance to the nearest infested forest [21].Despite just 100 m buffer zones of forest in each year, it is still possible to test the distance to the nearest declined forest inside these buffer zones.In addition, we tested the direction to the nearest declined forest in the previous year.The direction azimuth value was transformed to eight 45 • wide classes (N, NE, E, ES, S, SW, W and NW).Distance to the nearest clear-cut area and distance to the nearest natural forest edge were taken as a static variable from aerial photograph acquired in 1991.We hypothesize that the clear-cut edges might ease the entry of bark beetles into the stand while the natural forest edges represent a physical barrier for the beetles.Moreover, the trees in forest edge adjacent to clear-cut may suffer from water stress due to higher amount of solar radiation [20,54].The degree of forest "naturalness" was calculated as a difference between current forest inventory and potential natural tree species composition [55], derived from modified Czech typological maps (natural forest types [56]), which followed a similar concept as potential natural vegetation [57].Stand age represents an important attribute and we used age of dominant (>50%) tree layer because bark beetles preferentially attack trees older than 60 years [44].If any tree layer does not reach 50%, we calculated age as a mean value from the first and second dominant tree layer age.Similarly, the dominant tree layer was used for determining the tree density.In fact, there was mostly one dominant tree layer and one small layer of young trees.Soil is characterized by edaphic category, using the Czech typological system.It combines soil texture, nutrient richness and soil moisture [58].The age, tree density and edaphic category are from the Forest Management Institute (FMI) database (www.uhul.cz).The geology parameters are from Czech Geological Survey maps (www.geology.cz).We calculated the distance to two merged geological classes (7th moor, peat and 8th sediment), because we anticipated that the distance to these geology classes reflect wet or aquatic conditions or water supply.Altitude and other topographic variables were derived from a digital elevation model (DEM) with high spatial resolution (5 m) and vertical error of 1 m (DMR 4G from ČÚZK, www.cuzk.cz).We have used the heat load index (HLI) which describes the variation in potential solar radiation and increasing temperature on SW slopes [59], and the area solar radiation (ASR12), which is a cumulative value of potential solar radiation for entire year (watt hours per square meter-Wh/m 2 ) with maximum value on south slopes.Topographical wetness index (TWI), reflecting potential wetness was calculated using SAGA GIS 2.1.0.We expect that water availability may positively influence the tree resistance against the bark beetle infestation.

Statistical Models
Bark beetle dispersal results in gradual infestation of adjacent stands from year to year, which carries the possibility of changes in predictor significance during the outbreak period [21].To assess predictor significance changes, we built an individual statistical model of bark beetle infestation for every year within the studied period.The model for a given year included declined and living forest limited by the 100 m buffer zone around declined forest in the previous year because more than 90% of affected trees for every year were located in this zone.Because of gaps in data on bark beetle dispersal for 1992-1993, we used a 250 m buffer zone between the years 1991 and 1994.
To assess if pre-disturbance spectral trajectories improve the ability of a model to predict bark beetle infestation, we first composed a maximum model that contains all environmental variables plus all spectral trajectories (slopes).We used the square root transformation of the ASR12 variable to get values comparable in order to those of the other variables.To assess potential collinearity between any of the candidate predictors, we looked at the correlations between all combinations of variable pairs, as well as variance inflation factors (VIFs), calculated from the corvif function, which quantifies the severity of multicollinearity [61].The highest correlation we found was between Brightness and Greenness slopes, ranging between −0.89 and −0.97 for various years.Also, Greenness slope achieved the highest VIF value every year, ranging between 80 and 130, but reaching as high as 430 in one year.A strategy to avoid multicollinearity is to remove the predictor with the highest VIF value and recalculate the VIFs, repeating this process until all VIFs are smaller than a threshold value commonly set to 3-10 [61].Both approaches suggest removing Greenness slope from our subsequent analyses.Doing this, maximum observed correlations between any two remaining variables reach 0.6-0.7,and all VIFs are mostly below 3, with one or two VIFs reaching 3.1-3.5 in some years.Thus, except for Greenness slope, we keep all explanatory variables for the subsequent statistical modeling.
Since we fix three explanatory variables as mandatory parts of any model (DEM, Dist-source, ASR12), 12 other explanatory variables are left as candidate predictors.From the previous analyses we found that these three variables have relatively stable contribution to the models across years.While the DEM has relatively low importance, the latter two variables played dominant roles for models.Because these variables do not show any trend between the years, we decided to fix them.Using the all-subset approach, we investigate the effect of the pre-disturbance variables by running a total of M = 2 12 = 4096 candidate models.
Since our response variable is binary (30 m × 30 m forest on pixel is declined or not), we use binary logistic regression for fitting a regression curve.Our data for years 1995 and 2000 suffer from the problem of complete separation, meaning that a predictor or a subset of predictors is associated with only one outcome value when the predictor is greater than a threshold.This occurs commonly in small data sets when an event is rare, and this is also the case here.To cope with this, we use the Bayesian generalized linear model approach (bayesglm function) implemented in the library arm in the statistical software R version 3.1.1[62,63].We used the default procedure settings, i.e., an independent Cauchy prior distribution for the model coefficients.

Multimodel Inference
Because of high uncertainty in model selection (several candidate models or even several dozens of candidate models may describe data similarly well), we used the multimodel inference method as a basis for investigating explanatory and predictive performance of the considered model variables.In this approach, inference is based not on one, but several models in the candidate set [64][65][66].
In particular, the Akaike information criterion (AIC) is calculated for each model in the candidate set (we actually calculate AICc, since our sample sizes are relatively small) and these models are ranked according to their Akaike weights, defined for a model i as: In this formula, ∆AIC is the difference between AIC of the model with the minimum AIC value and AIC of each respective model.The Akaike weight w i expresses the likelihood that the model i is the best approximating model [66].We note that AIC (and any of its variants) provides only a relative measure of how good a model is, given the data and the candidate model set.
The Akaike weights allow estimating relative importance of every predictor variable in the full model by summing the Akaike weights for each model in which an individual variable occurs [64,66].We calculated this for all predictor variables, including with and without the pre-disturbance Wetness and Brightness slopes, allowing assessment of importance of the remote sensing predictors relative to the environmental variables.Similar to the Akaike weight, this relative importance (or predictor weight) can be interpreted as equivalent to a likelihood that the respective variable is a part of the best model [64,66].
We conducted the above-described procedure for four nested model groups.The first group consists of all M candidate models (full model set).The second group consists of all candidate models in which the predictor variable Brightness slope is absent and the only candidate pre-disturbance predictor is the Wetness slope (Wetness slope model set).Similarly, the third group consists of all candidate models in which the Wetness slope is absent and the only candidate pre-disturbance predictor is Brightness slope (Brightness slope model set).Our final set of models did not contain any pre-disturbance variable (plain model set).We used these four model sets to assess predictive importance of the Wetness slope and Brightness slope variables.All calculations involving multimodel inference were performed using the R package MuMIn [67].

Predictions
From each of the four model sets (the full, Wetness slope, Brightness slope and plain model sets), candidate models with the highest Akaike weights (with summed Akaike weight of at least 0.95) were selected.These were then used to predict probability of bark beetle attack in the years different from those for which the models were fitted.This was done to assess whether models fitted for a specific year in the outbreak cycle could be used to predict dynamics of the bark beetle attack in its other part of the cycle or how robust the fitted models are across the outbreak.These "best" model sets contained tens of candidate models, as opposed to thousands of models in the original model sets.We used models fitted for the year F (referred to as F-models) to calculate probability of bark beetle attack in the year P (referred to as P-predictions) as follows.Any F-model was used individually to calculate P-predictions, using data from the year P as its "validation" data.These individual P-predictions were then averaged, weighting specific P-predictions by the respective F-model Akaike weights [64,66].
Binary logistic regression models allow modeling binary responses (presence/absence of infestation) by linear regression via a logistic link function.As such, they can predict the probability of infestation, which can then be classified into two discrete classes choosing a probability threshold [15].Receiver Operating Characteristics (ROC) curves are a way of evaluating such binary models, since they are threshold-independent [68,69].The area under the ROC curve (AUC) is often taken as a typical performance measure for a binary response variable since it provides a single measure of classification accuracy [68,69].We quantified a predictive power of our model-averaged predictions using the ROC curves and AUC values.In particular, to assess an impact of considering the pre-disturbance spectral trajectories for making predictions, we statistically tested the null hypothesis that relevant pairs of the ROC curves (e.g., when both pre-disturbance variables are used vs.only one such variable is used) are statistically insignificant.We used the R package pROC [70] for calculating the ROC curves, AUC values, and statistical differences between any two ROC curves.

Pre-Disturbance Forest Condition
Both pre-disturbance predictors Wetness slope and Brightness slope show a similar spatial pattern (Figure 3).The spatial variability of both of these predictors is quite high and the slope values vary from negative to positive, with a clear spatial autocorrelation.Except for the highest Wetness slope and Brightness slope in affected stands at the beginning of the outbreak relative to latter years of infestation [40], the most obvious observation here is that the slopes of both indices vary most near non-forested areas (Figures 3 and 4).While values of the Wetness slope decrease near clear-cuts and locally increase near natural non-forested areas, the Brightness slope values demonstrate quite opposite behavior (Figures 3 and 4).Its values increase near clear-cuts and decrease near natural non-forested areas.Still, the correlation between Wetness and Brightness slopes is not that high (r = −0.54).The linear regression of Wetness from 1984 to 1994 has mean values of coefficient of determination for pixels with affected forest higher (R 2 = 0.35; SD 0.32) than pixels with unaffected forest until 1994 (R 2 = 0.17; SD 0.33).Similar results show the linear regression of Brightness from 1984 to 1994 for affected forest, which has higher coefficient of determination (R 2 = 27; SD 0.37) than unaffected (R 2 = 0.14; SD 0.38).More detailed information about the linear fitting of pre-disturbance trajectories in our study area is given in [40].
Remote Sens. 2016, 8, x FOR PEER 9 of 21 locally increase near natural non-forested areas, the Brightness slope values demonstrate quite opposite behavior (Figures 3 and 4).Its values increase near clear-cuts and decrease near natural nonforested areas.Still, the correlation between Wetness and Brightness slopes is not that high (r = −0.54).
The linear regression of Wetness from 1984 to 1994 has mean values of coefficient of determination for pixels with affected forest higher (R 2 = 0.35; SD 0.32) than pixels with unaffected forest until 1994 (R 2 = 0.17; SD 0.33).Similar results show the linear regression of Brightness from 1984 to 1994 for affected forest, which has higher coefficient of determination (R 2 = 27; SD 0.37) than unaffected (R 2 = 0.14; SD 0.38).More detailed information about the linear fitting of pre-disturbance trajectories in our study area is given in [40].locally increase near natural non-forested areas, the Brightness slope values demonstrate quite opposite behavior (Figures 3 and 4).Its values increase near clear-cuts and decrease near natural nonforested areas.Still, the correlation between Wetness and Brightness slopes is not that high (r = −0.54).
The linear regression of Wetness from 1984 to 1994 has mean values of coefficient of determination for pixels with affected forest higher (R 2 = 0.35; SD 0.32) than pixels with unaffected forest until 1994 (R 2 = 0.17; SD 0.33).Similar results show the linear regression of Brightness from 1984 to 1994 for affected forest, which has higher coefficient of determination (R 2 = 27; SD 0.37) than unaffected (R 2 = 0.14; SD 0.38).More detailed information about the linear fitting of pre-disturbance trajectories in our study area is given in [40].

Predictors of Bark Beetle Infestation
Statistical modelling based on the all-subsets approach and multimodel inference suggests that when both pre-disturbance predictors are taken into account, Wetness slope is an important variable to be considered for understanding the bark beetle infestation patterns (Table 3).Conversely, the likelihood that Brightness slope is a part of the best model is often as low as 0.3 (Table 3).In addition, whereas the high importance weights of Wetness slope are relatively stable over the years (with the exception of 2000), the Brightness slope weights show an increase in the years 1994-1996 and then a steady decline to lower values (≥0.28, see Table 3).The weights of other environmental predictors show relatively high variability over the years.However, similar to the pre-disturbance predictors, we recognize two basic trend-groups (Table 3).One group of predictors has a relatively high weight in the initial years of bark beetle infestation, while gradually losing weight towards the last modelled year.This group includes the edaphic category, which has high weight values until 1998 (with the exception of 1995); heat load index (HLI) and distance to wet or aquatic conditions (Geol-wet) show similar trends.The second group of predictors (TWI, AGE, and Dg-natur) reach the highest importance weights in the middle of bark beetle outbreak (1996)(1997) when the largest forest areas are declined in our study area [45].A special case of the opposite trend is represented by distance to clear-cut (D-clear-cut), which has the highest weights at the beginning and the end of the outbreak.The remaining predictors show somewhat fluctuating weights over the years.
To examine how explanatory power of models that include our two considered remote sensing variable changes relative to models in which one of these variables is excluded, we repeated the above analysis considering all environmental predictors but only with either Wetness or Brightness slope.With just Wetness slope (i.e., without Brightness slope), its importance weight and those of all environmental variables remained virtually unchanged (compare Table 3 and Table A1).On the contrary, with just Brightness slope (i.e., without Wetness slope), its importance weight increased in several years (1994,1996), but decreased in other ones (1995,1997; compare Table 3 and Table A2).These changes in the Brightness slope importance weights were accompanied by changes in the HLI importance weights.In particular, the HLI importance weight generally increased, with the highest increase in the years 1996 to 1998 being in the range of 10% to 20%.Similar importance weight changes were observed when the results based on the Wetness model set were compared with those based on the plain model set (high importance of Wetness slope; Table A3), and when the results based on the Brightness model set were compared with those based on the plain model set (low to medium importance of Brightness slope; Table A3).

Model Predictions
We assessed the predictive power of our fitted models by calculating the ROC curves and the corresponding AUC values for all years except the training year of the model.Regardless of the modelling scenario used, the trend is quite similar.Table 4 shows the results of models including both Wetness and Brightness slope.Generally, the predictive power of most fitted models decreases as time progresses.Further, models fitted at the beginning of bark beetle outbreak have higher predictive power than models fitted after outbreak culmination.The models describing an initial phase of bark beetle outbreaks thus appear more reliable for conducting near-future predictions.The higher predictive power of a few models for the year 2000 breaks this rule, but at the end of outbreak, reliability of models somewhat decreases because of low sample size.Table 4. Ability of fitted full models for each individual year (rows) to predict bark beetle forest infestation in the years 1994 to 2000 (columns), except the same year for which the model was fitted, evaluated using AUC.Higher AUC values mean higher predictive power of the respective models; value 0.5 represents random prediction.Inclusion of the pre-disturbance indices generally increases model predictive power (Table 5).Although the improvement may seem marginal, it is statistically significant in many cases (Table A4).Interestingly, an increase is generally higher with backward prediction (Table 5, values below the diagonal) than forward one (Table 5, values above the diagonal).Differences in Table 5 would not change much if we replace the plain model set with the Brightness model set.On the other hand, if we consider differences between AUC values corresponding to the full model set and the Wetness model set, they are all virtually zero and statistically insignificant.This also suggests that with respect to the model predictive power, the Wetness slope plays a more important role than the Brightness slope.

Spectral Indices Reflecting the Forest Health Change
Our results clearly demonstrate that long-term changes of TC Wetness and Brightness are related to the probability of bark beetle infestation during outbreak.This corroborates our hypothesis that bark beetles preferably affect the weakened forest stands, and that these changes in forest health or tree vigor are detectable in spectral signal change.In general the nature of spectral signal change is connected with the tree defoliation [37], chlorophyll fluorescence change and/or water content decrease in needles due to water stress [71].This spectral change is well detected in NIR and SWIR spectral region [11].In particular the SWIR-based indices (NDMI, Wetness, Brightness, NBR) are often used to assess change in forest health due to acid rain [72], defoliation following insect attack (e.g., [51]), or forest mortality (e.g., [73]).For this study, we selected Wetness and Brightness indices.These two indices show almost opposite behavior in reaction on forest decline [45].In contrast, some authors have reported success using NIR-based indices.For example, Senf et al. [37] found stronger association of spruce budworm disturbance with Greenness than two other TC indices.However, the same study describes the highest sensitivity of Wetness and Brightness to mountain pine beetle disturbance.We removed Greenness early during our statistical analysis due to correlation and collinearity.Moreover, previous results from the same study area [45] support this removal, because they showed the inconsistent response of Greenness to the forest decline, when the magnitude of change was lower and the values changed with delay after forest decline.A possible reason could be the variability in understory composition; a dying canopy may uncover either soil and litter, or a green shrub/herbaceous layer.The sensitivity in spectral response to these contrasting cover types may be greater in Greenness than Wetness or Brightness and lead to inconsistent spectral-temporal trends associated with infestation-related tree mortality.
Further understanding of ecological processes underlying changes in Wetness and Brightness can be inferred from index trends we observed near forest edges.The decrease of Wetness and increase of Brightness and/or high variability of both predictors near clear-cuts may be related to water stress due to higher amount of solar radiation, particularly in southern and western forest edges [20,54].Similar results were reported by Dantas de Paula et al. [74] in a Landsat-based study, where the tree cover variation increases towards the forest edge.The authors suggested the underlying causes for this variability to be degradation debt, hyperdynamism and environmental factors (topography, soil conditions etc.).In contrast to clear-cuts, the natural forest edges do not show obvious change of variability in the pre-disturbance indices.Nevertheless, some of these edges have higher Wetness slope values than do areas inside the forest.This effect is pronounced at the margin between forest and meadows (e.g., Lusen valley), where the high canopy closure protects the soil against drying, or denser canopies may serve as a physical defense to bark beetle infestation.

Pre-Disturbance Forest Conditions
If the trend in both indices indicates changes in forest health status or stand vigor, we should be able to describe their main causes.One of the strongest factors responsible for forest decline in Central Europe and North America in the second half of the last century was anthropogenic air pollution causing acid deposition.Even though the Šumava Mountains do not belong to the most heavily polluted areas in Czech Republic (e.g., the "black triangle" in northern Czech Republic with large historical forest decline), there is evidence from ground observations of forest decline due to air pollution in the last century [75,76].Other retrospective methods suggest the consequences of acid deposition on forest health can be still recognized.Šantr ůčková et al. [76] found a proportional relation between ∆ 13 C in tree rings and soil acidification (pH decrease).Polák et al. [77] describe tree defoliation followed by crown regeneration using secondary shots that were detectable even after 30 years.It is therefore likely that multi-temporal imagery captures this forest health decline due to acid deposition, which is in concordance with results reported by Rock et al. [72].Another possible long-term stress factor in the Šumava Mountains is stand damage due to repeated windstorms or snow [7].The repeated windstorms may damage the tree roots and cause tree vigor decrease.Drought is another possible long-term factor [28], and its role probably increases during bark beetle infestation.
All the aforementioned factors may reflect the spatial variability resulting from climate and topography differences from local to regional scales.Even though some studies mention no significant relation between forest health and bark beetle attack (e.g., [78]), Bytnerowicz et al. [79] suggested that long-term elevated levels of atmospheric N and S depositions and elevated O 3 predispose trees to insect attacks and other stresses.

Spectral Trajectories
The use of linear regression for Wetness and Brightness fit impose simple trends over time.This assumption is consistent with results of Hais et al. [40], who reported lower Wetness slope values in the first years of the bark beetle outbreak in comparison to latter years, which were more stable.This suggests that the initial bark beetle attack occurred in the most vulnerable stands that show long-term Wetness/Brightness change before the disturbance.Nevertheless, the resulting spectral trajectory may be influenced by other factors, which can act antagonistically.While long-term stress factors decrease Wetness values (or increase Brightness) over time, canopy growth may support an opposite trend, in particular for fast-growing tree species.For example, partial forest defoliation may result in small spectral change, which is quickly followed by spectral recovery after disturbance [37,38].Similarly, forest thinning is reflected in Wetness change [16] when the values first decrease due to lower canopy closure (and higher role of ground cover in the spectral signal), but after several years, the values increase due to faster tree growth [16].Users should therefore be careful when using these indices in managed forests.Additionally, recent forest history often includes high magnitude changes that are not always linear in their spectral-temporal response.Because of this complexity, the non-linear fit of spectral trajectory is often needed to describe forest history for a given pixel [35].Nevertheless, much easier interpretable linear fits are valuable for description of many long-term processes with unaltered trend like gradual systematic change in natural vegetation communities [80].Similarly we used the linear fit to describe pre-disturbance forest condition assuming it mostly follows simple trend.However, to avoid complexity in the course of spectral-temporal response we use just forest pixels without any type of disturbance during assessed pre-disturbance time.

Risk of Bark Beetle Infestation and Model Predictions
Among all predictors used in our study, Wetness slope provided the greatest improvement of bark beetle infestation risk modelling, while Brightness slope was a significant predictor only for a few years.Moreover, Wetness slope and Brightness slope are not correlated with the environmental predictors, and thus are able to explain additional observed variability.The Landsat spatial resolution (30 m) provides a good dataset to describe the spatial variability on a stand level.Relatively low weights of importance for some environmental variables may be caused by the small dataset, especially when the spatial resolution of the data is relatively coarse and/or the predictor is categorical.This may be the case for very low weights of broad edaphic categories in models for years 1995 and 2000 (lowest number of cases), which are high in other years.This limitation has lower importance for continuous variables like pre-disturbance indices.This supports the potential of the Landsat data for local to global-level modelling in general.In agreement with Lausch et al. [21], we conclude that the importance of individual predictors changes during the bark beetle infestation.The predictors with the highest model weights (edaphic category, Geol-wet, HLI) in the first years of the bark beetle outbreak document the selective infestation of the most susceptible forest stands.The selective infestation is critical until bark beetle population size reaches the threshold where tree defense is overwhelmed [78].After successful selective infestation, bark beetle propagation produces population sizes so great that selective behavior is an unnecessary strategy to seize trees, which shifts increased importance to other predictors.This was demonstrated in 1997, when the largest forest area was affected, and tree-related predictors had the highest weights (age, density, DG-natur).The trees were probably more susceptible individually, reflecting stand structure and not the environment predictors as in the beginning of the bark beetle outbreak.The decrease in almost all predictor weights further into the future, shows that the high bark beetle density led to random-like infestation, where the beetles were able to overcome any tree's defense.The predictor D-clear-cut has high importance in the first and last years of bark beetle outbreak, which might be connected with low bark beetle numbers in both periods and need for selective behavior.
Predictive power of our fitted models (including all predictors) calculated using ROC curves and corresponding AUC values was quite similar.The high AUC values when data were used for fitting and also for prediction (0.76-0.91) documented reasonably good fit to the data, when assessed using AUCs.The decrease of predictive power in fitted models over the years documents not only the general decrease of fits to the distant years, but it also supports the explanation of predictor importance changes during a bark beetle outbreak [21].Moreover, the slower decrease of predictive power (higher predictive power for distant years) for models fitted at the beginning of the bark beetle outbreak nicely corresponds with the hypothesis that at the beginning of the bark beetle outbreak beetles are more selective, while later, when beetle densities rise and/or the most susceptible stands are already destroyed, beetles need to be more opportunistic.

Conclusions
Using pre-disturbance slopes of Tasseled Cap components Wetness and Brightness, spectral trajectories improved our understanding and prediction of bark beetle infestation.Linear regression slopes of these two variables calculated for the time period before the disturbance significantly increase the predictive power of bark beetle infestation models from year to year and even for longer predictions.Moreover, these predictors were not correlated with environmental predictors and therefore are able to explain additional observed variability.This also suggests that long-term gradual decrease of forest vitality significantly increases forest vulnerability to bark beetle infestation.When comparing the importance weights of both pre-disturbance predictors, the Wetness slope weights are high and relatively stable over the years (with the exception of 2000), while the weights of Brightness slopes show an increase only in the middle of the bark beetle outbreak (1996) and steadily decline to low values.We distinguished two components of relationship between these two indices and tree susceptibility to bark beetle infestation.First, they reflect long-term stress factors (acid deposition, drought, snow, windstorms).The second component is their change near forest clear-cuts edges.This probably relates to higher solar radiation at open edges or other disturbances caused during the clear-cut.
We have also identified two groups of environmental factors important for prediction of bark beetle infestation risk.The first group has high importance weights in the initial years of the outbreak (edaphic category, distance to aquatic conditions, heat load index), which might be connected with high pest selectivity for the most susceptible stands.The second group of tree-related predictors has the highest weights in the middle of the outbreak, when the largest forest area was affected, and the pest numbers were highest (age of dominant tree layer, tree density of dominant tree layer, Distance to the natural non-forested area).The first can be used to identify potentially vulnerable areas, and the second for predicting bark beetle spread after the initial phase of infestation.
But forest ecosystems also face long-term, gradual changes due to hidden stress factors that contribute to their vulnerability and susceptibility to bark beetle infestation.And this is where non-correlated pre-disturbance spectral trajectories, based on Landsat data with its vast and accessible archive, can really help.

Figure 1 .
Figure 1.Study are a locate d in the Šumava Mountains, Cze ch Re public -Ge rmany borde r re gion.

Figure 1 .
Figure 1.Study area located in the Šumava Mountains, Czech Republic-Germany border region.

Figure 2 .
Figure 2. Spe ctral traje ctory of We tne ss TC compone nt in a give n pixe l calculate d as the slope of linear re gre ssion of We tne ss for conse cutive ye ars from 1984 to one ye ar be fore the year for which the model was calibrate d.Slope was calculate d for both infe ste d and non-infe ste d pixe ls; the infe ste d one is shown as an e xample .

Figure 2 .
Figure 2. Spectral trajectory of Wetness TC component in a given pixel calculated as the slope of linear regression of Wetness for consecutive years from 1984 to one year before the year for which the model was calibrated.Slope was calculated for both infested and non-infested pixels; the infested one is shown as an example.

Figure 3 .
Figure 3. Spatial distribution of We tne ss slope and Brightne ss slope value s.Similar to Figure 2, the pixe l value s we re calculate d as the slope of line ar re gre ssion of We tne ss/Brightne ss for consecutive ye ars from 1984 to one ye ar be fore the ye ar for which the mode l was calibrate d.The maps are compose d from pixe ls combine d among all e valuate d ye ars (1994-2000) according to the year of bark be e tle infe station.The lowe st We tne ss slope value s and highe st Brightne ss slope value s are mostly adjace nt to cle ar-cuts, while the slope value s in the vicinity of natural non-fore ste d are as show opposite be havior, e specially in the Luse n valle y.

Figure 3 .
Figure 3. Spatial distribution of Wetness slope and Brightness slope values.Similar to Figure 2, the pixel values were calculated as the slope of linear regression of Wetness/Brightness for consecutive years from 1984 to one year before the year for which the model was calibrated.The maps are composed from pixels combined among all evaluated years (1994-2000) according to the year of bark beetle infestation.The lowest Wetness slope values and highest Brightness slope values are mostly adjacent to clear-cuts, while the slope values in the vicinity of natural non-forested areas show opposite behavior, especially in the Lusen valley.

Figure 3 .
Figure 3. Spatial distribution of We tne ss slope and Brightne ss slope value s.Similar to Figure 2, the pixe l value s we re calculate d as the slope of line ar re gre ssion of We tne ss/Brightne ss for consecutive ye ars from 1984 to one ye ar be fore the ye ar for which the mode l was calibrate d.The maps are compose d from pixe ls combine d among all e valuate d ye ars (1994-2000) according to the year of bark be e tle infe station.The lowe st We tne ss slope value s and highe st Brightne ss slope value s are mostly adjace nt to cle ar-cuts, while the slope value s in the vicinity of natural non-fore ste d are as show opposite be havior, e specially in the Luse n valle y.

Figure 4 .
Figure 4. Wetness slope and Brightness slope values as they vary with distance to the nearest clear-cuts (A,B) and natural non-forested areas (C,D).The trend (blue line) of slope among distances is visualized using local polynomial regression (loess) ±1 standard error mean (grey area).

Table 2 .
List of predictors for bark beetle infestation including pre-disturbance indices describing forest susceptibility and predictors of bark beetle dispersal.

Table 3 .
Predictor importance weights based on multimodel inference when both pre-disturbance indices are used as predictors.The darkest/lightest green show the highest/lowers predictor importance weights.

Table 5 .
Difference between AUC values corresponding to the full model set (with both Wetness and Brightness slopes) and the plain model set (with no remote sensing variable).