Detection of Spatio-Temporal Changes of Norway Spruce Forest Stands in Ore Mountains Using Landsat Time Series and Airborne Hyperspectral Imagery

The study focuses on spatio-temporal changes in the physiological status of the Norway spruce forests located at the central and western parts of the Ore Mountains (northwestern part of the Czech Republic), which suffered from severe environmental pollution from the 1970s to the 1990s. The situation started improving after the pollution loads decreased significantly at the end of the 1990s. The general trends in forest recovery were studied using the tasseled cap transformation and disturbance index (DI) extracted from the 1985–2015 time series of Landsat data. In addition, 16 vegetation indices (VIs) extracted from airborne hyperspectral (HS) data acquired in 1998 using the Advanced Solid-State Array Spectroradiometer (ASAS) and in 2013 using the Airborne Prism Experiment (APEX) were used to study changes in forest health. The forest health status analysis of HS image data was performed at two levels of spatial resolution; at a tree level (original 2.0 m spatial resolution), as well as at a forest stand level (generalized to 6.0 m spatial resolution). The temporal changes were studied primarily using the VOG1 vegetation index (VI) as it was showing high and stable sensitivity to forest damage for both spatial resolutions considered. In 1998, significant differences between the moderately to heavily damaged (central Ore Mountains) and initially damaged (western Ore Mountains) stands were detected for all the VIs tested. In 2013, the stands in the central Ore Mountains exhibited VI values much closer to the global mean, indicating an improvement in their health status. This result fully confirms the finding of the Landsat time series analysis. The greatest difference in Disturbance Index (DI) values between the central (1998: 0.37) and western Ore Mountains stands (1998:  ́1.21) could be seen at the end of the 1990s. Nonetheless, levelling of the physiological status of Norway spruce was observed for the central and western parts of the Ore Mountains in 2013 (mean DI values  ́1.04 (western) and  ́0.66 (central)). Although the differences between originally moderately-to-heavily damaged, and initially damaged stands generally levelled out by 2013, it is still possible to detect signs of the previous damage in some cases.


Introduction
Forests play an important role in regulating global climate; moreover, they provide human beings with a whole range of ecosystem services.In Europe forest health and functioning have been primarily influenced by anthropogenic activities (e.g., air pollution, surface mining, contamination).The situation of deteriorated forest health status has led to numerous remote sensing applications focusing on monitoring forest damage in the last 40 years, especially after the launch of the ERTS/Landsat program [1,2].Remote sensing vegetation studies are based on canopy-level vegetation spectral properties, which are determined by a combination of factors, including canopy density and structure, as well as foliar pigments, water, and other constituent contents [3].The structure of a forest stand is also directly linked with the stand's disturbance history [4].
The Landsat program offers an opportunity to study the long-term development of forest ecosystems through the availability of a long-term data series (since 1972).Numerous forest monitoring methods have been developed, including simple techniques, such as classification or vegetation indices analysis [5,6] as well as more sophisticated approaches, such as change vector analysis [7,8], trajectory based segmentation [9,10] or multi-temporal spectral mixture analysis [11].Although the Landsat data have been used for several decades, their use for forest monitoring is still relevant [12][13][14][15].
The Landsat data can be widely used in forest monitoring; however, their applicability for detailed vegetation studies is limited, especially by their spatial and spectral resolution.Classification of forest status using the Landsat data is, therefore, too general, usually allowing only very basic forest damage classes to be distinguished, such as no damage, intermediate damage, and heavy damage/dead forest [1].Significant advances have been achieved with the use of imaging spectrometers enabling identification of not only significant forest damage, but also the very initial phases of vegetation stress due to detailed records (narrowband) of vegetation spectral signatures [16][17][18].The availability of airborne hyperspectral imagery has, thus, led to the development of techniques for detailed studies and modelling of selected biochemical, biophysical, and structural vegetation parameters, thus enabling not only a qualitative comparison, but also a quantitative estimation [19][20][21][22].Many of these techniques are based on regression between hyperspectral vegetation indices and vegetation properties [23][24][25][26].A large part of vegetation indices are usually affected simultaneously by several vegetation parameters [25,27], so these indices can be used as indicators of the overall vegetation status rather than as predictors of particular vegetation parameters.Nevertheless, this is an advantage when looking for non-specific vegetation stress markers.Therefore, hyperspectral data have been successfully used for detecting vegetation stress and damage of various levels and causes.Huesca et al. [28] used a combination of AHS airborne hyperspectral data and MODIS superspectral satellite data to study recovery of forest vegetation after fire by the means of spectral unmixing analysis and the Normalized Difference Vegetation Index (NDVI) change detection technique.Wu et al. [29] studied the effects of late frost stress using a wide range of spectral parameters extracted from the red-edge spectral domain.Fassnacht et al. [30] applied a combination of a genetics algorithm (a heuristic search approach inspired by natural evolution generating solutions for optimization problems) and a supervised support vector machine on airborne HyMap hyperspectral data to classify different tree mortality classes caused by a bark-beetle outbreak.Mahlein et al. [31] developed a set of specific spectral disease indices optimized for various crops.Sanches et al. [32] developed a Plant Stress Detection Index based on the parameters of the continuum removed chlorophyll absorption feature which was then applied at both leaf and canopy level to reveal stressed plants.Behmann et al. [33] used a Support Vector Machine and a set of relevant vegetation indices to quantify distribution of progressive senescence stages and to distinguish well-watered and drought stressed plants.They pointed out that drought stress can be detected much earlier in such a case than in the case of using NDVI.Nevertheless, one of the disadvantages of using hyperspectral data is the fact that due to high costs they are usually not acquired repeatedly over a long period, so as a stand-alone they cannot be used for monitoring the change in forest health status.
The main aim of this paper is to assess spatio-temporal changes of the physiological status of homogeneous (i.e., single species and even aged) Norway spruce (Picea abies L. Karst) forests located at two sites of the Ore Mountains located in the northwestern part of the Czech Republic differing by forest disturbance history and pollution load.The study was based on a long-term (1985-2015) series of Landsat imagery in combination with analysing the physiological status of the forest stands under study using the airborne hyperspectral data acquired during the peak of the disturbance level (1998) and after 15 years (2013).The image information thus derived is then coupled with the results of an in situ field survey organized together with both hyperspectral data acquisitions.The focus was on answering whether there is currently any difference in tree physiological status between the originally damaged stands (central Ore Mountains) and those stands located in the western Ore Mountains region, which had not been found to be seriously damaged in 1998.
The present study is a direct continuation of the work conducted in these sites by Campbell et al. [17] who used ASAS airborne hyperspectral imagery for detection of initial stress in Norway spruce stands and for establishing the links between various forest parameters (describing the particular damage levels) and the image-extracted spectral signatures.

Study Area
The study was performed at two sites situated in the Ore Mountains in Western Bohemia (Figure 1).The central part of the Ore Mountains region was suffering high loads of SO 2 , NO X , and PM 10 [34] originating mainly from coal mining and heavy industry.The high levels of air pollution led to massive dieback of the local Norway spruce forests between 1970 and 2000.The position of the pollution sources and the prevailing eastward wind direction contributed to a sharp west-to-east increasing pollution gradient and a corresponding gradient in tree damage levels.In the early 1990s, the west-east gradient of imissions in the Ore Mountains was SO 2 81.4 µg¨m ´3 (central) and 24.0 µg¨m ´3 (western); NO X 34.8 µg¨m ´3 (central) and 16.5 µg¨m ´3 (western) [34].Although the pollution loads had fallen by 1998, the difference between the central and western Ore Mountains still remained (SO 2 17.0 and 6.0 µg¨m ´3; NO X 18.0 and 14.0 µg¨m ´3) [34].There was an observed relationship between loss of forest and the prevailing wind directions.As such, the central and western parts exhibited 70% and 50% forest cover loss between 1972 and 1989, respectively [35].The forest cover loss was due to high mortality and rapid "emergency cuttings" of heavily damaged trees.The emission loads dropped significantly from the end the of 1990s; however, the west-east gradient was still observable in 2012: SO 2 7.4 µg¨m ´3 (central), 2.5 µg¨m ´3 (western); NO X 11.5 µg¨m ´3 (central), 7.3 µg¨m ´3 (western) [34].
Remote Sens. 2016, 8, 92 3 and after 15 years (2013).The image information thus derived is then coupled with the results of an in situ field survey organized together with both hyperspectral data acquisitions.The focus was on answering whether there is currently any difference in tree physiological status between the originally damaged stands (central Ore Mountains) and those stands located in the western Ore Mountains region, which had not been found to be seriously damaged in 1998.
The present study is a direct continuation of the work conducted in these sites by Campbell et al. [17] who used ASAS airborne hyperspectral imagery for detection of initial stress in Norway spruce stands and for establishing the links between various forest parameters (describing the particular damage levels) and the image-extracted spectral signatures.

Study Area
The study was performed at two sites situated in the Ore Mountains in Western Bohemia (Figure 1).The central part of the Ore Mountains region was suffering high loads of SO2, NOX, and PM10 [34] originating mainly from coal mining and heavy industry.The high levels of air pollution led to massive dieback of the local Norway spruce forests between 1970 and 2000.The position of the pollution sources and the prevailing eastward wind direction contributed to a sharp west-to-east increasing pollution gradient and a corresponding gradient in tree damage levels.In the early 1990s, the west-east gradient of imissions in the Ore Mountains was SO2 81.4 µg•m −3 (central) and 24.0 µg•m −3 (western); NOX 34.8 µg•m −3 (central) and 16.5 µg•m −3 (western) [34].Although the pollution loads had fallen by 1998, the difference between the central and western Ore Mountains still remained (SO2 17.0 and 6.0 µg•m −3 ; NOX 18.0 and 14.0 µg•m −3 ) [34].There was an observed relationship between loss of forest and the prevailing wind directions.As such, the central and western parts exhibited 70% and 50% forest cover loss between 1972 and 1989, respectively [35].The forest cover loss was due to high mortality and rapid "emergency cuttings" of heavily damaged trees.The emission loads dropped significantly from the end the of 1990s; however, the west-east gradient was still observable in 2012:   The central Ore Mountains region is represented by Norway spruce stands located in sites within the neighbouring the town of Kovářská and the Přísečnice water dam.The local forests showed a moderate to high level of damage at the end of the 1990s.In 1998 the trees had an average of 6 to 8 needle age classes for moderately to heavily damaged trees, respectively [17].Nevertheless, in extreme cases of heavily damaged trees only 4 needle age classes remained.Since the reduction in air pollution the process of forest recovery has stared and currently the Norway spruce stands show only initial to moderate damage symptoms.
The second study site is located in the western part of the Ore Mountains region in the neighbourhood of Přebuz town.In contrast with the central Ore Mountains this area was outside of the region of the highest pollution load and the local Norway spruce forests therefore remained almost intact (healthy to initial damage).In 1998 the trees had 9 to 10 needle age classes for healthy and initially damaged sites, respectively.In 2013 all the trees had an average of eight needle age classes.

Landsat Time Series
The general trends of forest disturbance/recovery (i.e., defining the periods of improvement or worsening of the forests physiological status) were studied using the time series of Landsat TM/ETM+/OLI data covering the time period from 1985 until 2015.The Landsat data were acquired from the USGS data archive via the Earth Explorer interface.Note that only the scenes with no clouds over the compared areas obtained during the highest vegetation season (June to September) were selected for further processing.The Landsat scenes used are listed in Table 1.

Airborne Hyperspectral Data
ASAS airborne hyperspectral data, consisting of 512 ˆ662 pixels, 62 spectral bands (with FWHM ca. 10 nm) covering the spectral range from 410 to 1013 nm, were acquired between 20 August and 1 September 1998.However, only the bands between approx.550-850 nm were calibrated to reflectance and used for further analyses due to poor atmospheric conditions, affecting the spectral calibration of the sensor in the VIS-B and far NIR regions.The original ground spatial resolution of this dataset was 1.5 ˆ2.0 m.ASAS radiance data was calibrated to at-surface reflectance using a vicarious hybrid calibration approach utilizing the 6S atmospheric correction code, as well as in situ ground reflectance measurements from dark and bright calibration targets (asphalt and lime) collected by a GER-2500 field spectroradiometer [17].
APEX (Airborne Prism Experiment) data, consisting of 286 spectral bands (with FWHM ca.7 nm) covering the spectral range from 430 to 2450 nm, were acquired on 6 September 2013 consisting of 286 spectral bands (with FWHM ca.7 nm) covering the spectral range between 430 and 2450 nm.The original ground spatial resolution of this dataset was 2.0 m.A supportive calibration/validation ground campaign was organized simultaneously with the acquisition of hyperspectral imagery.Surface level spectra of several calibration and validation targets of an adequate spatial extent and homogeneity were taken by an ASD Fieldspec-4 spectroradiometer.The targets included water bodies (2ˆ), asphalt (3ˆ), concrete tiles (4ˆ), bare soil/sand (1ˆ), fresh grass (2ˆ), and dry grass (1ˆ).APEX data pre-processing, including radiometric calibration, geometric, and atmospheric corrections was performed by VITO (Flemish Institute for Technological Research).Radiometric calibration was conducted using APEX Calibration Home Base data [36].Geometric correction was performed using the data from a GPS/IMU unit, boresight corrections, and a 5 m grid LiDar-based digital elevation model.Atmospheric correction and other necessary pre-processing (e.g., a smile-aware atmospheric correction and retrieval of hemispherical-conical reflectance factors-HCRF) was conducted within the Central Data Processing Centre (CDPC) as well as at VITO using a module based on MODTRAN-5 radiative transfer code.A detailed description of the whole APEX data pre-processing provided by CDPC was given in [37].
The challenge when comparing the images acquired with different optical sensors is arriving at a common set of spectral and spatial parameters.Identifying such a common set of parameters for comparison is a key step, which may not be possible when using multispectral sensors, providing discrete spectral bands.However, imaging spectrometers provide spectrally continuous calibrated reflectance data (e.g., for ASAS and APEX ~10 nm resolution, 550-850 nm common region) to compute comparable spectral parameters (VIs).To address the difference in the calculated VIs due to the sensitivity of the two instruments (e.g., signal strength and noise characteristics), and the different correction procedures the computed VIs were normalized as described in Section 3.2.

Field Reference Data
A supportive ground campaign was organized simultaneously with the ASAS data acquisition in late August 1998 as described in [17].In total, 38 even aged mature Norway spruce stands (21 at the western and 17 at the central Ore Mountains) were selected for evaluation of tree physiological status using the five-grade damage class scale as described in Table 2. Forest physiological condition was assessed using a methodology for evaluating the state of individual tree crowns, modified for use with high spectral resolution digital imagery and field surveys [17].The classification scheme [38] is based on the evaluation of the degree of defoliation and the presence of chlorosis in individual tree crowns.The crown characteristics of a stand were then summarized to classify whole stands into damage levels.For more details see [17].In addition, five trees were selected at each stand for needle sampling to determine the content of photosynthetic pigments (chlorophylls and carotenoids).The needle samples were collected from the sunlit production portion of tree crowns using pruning poles.The collected samples were placed into the dark conditions of a portable freezer and transported to a laboratory for further processing.Photosynthetic pigments were extracted in dimethylformamide (DMF) in cold dark conditions following the procedure of [39].The foliar pigments content was then determined spectrophotometrically using the equations of [40] and normalized by dry sample weight (mg¨g ´1).All sampled trees fell within the 30 m ˆ30 m buffer area at each particular stand.The coordinates of the centroids of these buffer areas were measured using differentially-corrected GPS [17].The Norway spruce stands defined during the 1998 field campaign that were considered essentially as initially damaged (including the undamaged or initially damaged stands of DC 0 -DC 1 ; located at western Ore Mountains) or moderately to heavily damaged DC 2 -DC 3 ; located at the central Ore Mountains, were clearly distinguished using a visual assessment of the defolitation level and the presence of chlorosis as described in Table 2.Note that no stand was classified into the DC 4 class either at the western or at the central Ore Mountains.In 2013, a new ground campaign was organized simultaneously with APEX data acquisition.Eleven Norway spruce stands were selected for detailed study; of these, eight coincided with the stands selected during the 1998 campaign: P 9 , P 10 , P 11 , P 40 at the western Ore Mountains and K 25 , K 37 , K 46 , and K 47 at the central Ore Mountains A similar sampling design as in 1998 [17] was applied; five trees were sampled per stand to assess foliar pigments content.The defoliation level and presence of chlorosis were assessed for twelve trees per stand.All stands in both sites were considered as initially to moderately damaged (DC 1 and DC 2 ) since the trees showed a mean defoliation level of 30% in the central Ore Mountains and 33% in the western Ore Mountains in 2013.
To validate the results of the hyperspectral (HS) data analysis of ASAS and APEX image data, only those stands that could have been found in both ASAS and APEX datasets were selected (e.g., some of the stands were covered by clouds in case of the ASAS dataset or were harvested in the 1998-2013 period).This included a total of 16 stands: 14 at the western Ore Mountains (P 1 , P 2 , P 3 , P 4 , P 5 , P 6 , P 9 , P 10 , P 11 , P 15 , P 17 , P 43 , P 44 , and P 45 ) and two at the central Ore Mountains (K 33 and K 37 ) (Figure 1).

Methods
The focus was on Norway spruce forests and, thus, it was necessary to mask the other surfaces.Therefore prior to the analysis Norway spruce masks were derived for individual Landsat images as well as for the HS image data sets.These Norway spruce masks were then used for all following calculations and transformations applied to Landsat scenes as well as to both HS datasets.The long-term trend was studied using a Landsat time series , whereas a detailed analysis on the forest health status was conducted using the available multi-date hyperspectral (HS) image data acquired in 1998 and 2015.Due to technological differences between ASAS and APEX sensors (e.g., spectral resolution and signal-to-noise (SNR) ratio) and differences in atmospheric corrections applied to the ASAS and APEX data, the original VI values had to be normalized to allow their mutual comparability (see Section 3.2 Equation ( 6) and Section 3.3).Only after that was it possible to further assess the spatio-temporal trends.The forest health status analysis of multi-date HS image data was performed at two levels of spatial resolution of the hyperspectral imagery, at a tree level (original 2.0 m spatial resolution), as well as at a forest stand level (generalized to 6.0 m spatial resolution).A statistical assessment for the Landsat time series, as well as the HS data sets and field reference were done at a stand level.

Landsat Disturbance Index Time Series Analysis
The Landsat data were firstly transformed from the raw DNs into top-of-atmosphere (TOA) reflectance using the calibration tool included in the ENVI 5.2 software and the appropriate calibration coefficients included in the Landsat metadata files (MTL).Note that regarding [41,42] TOA, reflectance is recommended for the Tasseled Cap transformation rather than at-surface reflectance.
The Tasselled Cap transformation-TCT [43] was then performed to obtain three components related to vegetation characteristics: Brightness (B), Greenness (G) and Wetness (W).The TCT was performed using the appropriate tool included in the ENVI 5.2 software for the scenes acquired by Landsat 5 TM and Landsat 7 ETM+.Unfortunately, this version of the ENVI software does not allow the TCT transformation for the data acquired by Landsat 8 OLI.Thus, the TCT coefficients (describing the relationship of the original image bands and the TCT components) were used as recommended by [42].The coefficients optimized for the growing season were used for more details see [42].
In the next step, the Disturbance Index (DI) was calculated from the TCT components following the definition of [44] as: DI " B r ´pG r `Wr q (4) where: B r , G r , and W r are rescaled Brightness, Greenness, and Wetness components; B, G, and W are the original Brightness, Greenness, and Wetness values; B µ , G µ , and W µ are the means of B, G, and W calculated from the forest land areas; B σ , G σ , W σ are standard deviations of B, G, and W calculated from the forest land areas defined by the appropriate binary mask.The concept of DI assumes that high B r and low G r and W r values are typical for disturbed stands.The DI value is high (positive) in such a case.On the other hand, undisturbed, fully-regenerated and young vital stands are typical for low B r and high G r and W r values resulting in low (negative) DI values [44].
The Landsat 7 scenes were acquired in SLC-off mode after 2003 due to a malfunction of the Scan Line Corrector of the ETM+ sensor (compensating the orbital motion of the satellite).This results in stripes of zero-value pixels in the imagery.Thus, these areas were masked out and only the valid values were used for the all analysis.
A raster layer of the DI values calculated for the forest land areas was constructed for each Landsat scene.The DI values were then classified into three classes using mean (µ) and standard deviation (σ) values.Class 1 includes low DI values (less than µ ´0.5σ) representing undisturbed and most vital forest, whereas Class 3 includes high DI values (more than µ + 0.5σ) representing the most disturbed stands.Class 2 represents the most common DI values between µ ´0.5σ and µ + 0.5σ.

Vegetation Indices Derived from Hyperspectral Image Data and Stand Separability
The forest health status analysis was performed at two levels of spatial resolution of the hyperspectral imagery: original (2.0 m) considering individual tree crowns and generalized (6.0 m) considering whole forest stands.Therefore, the stands had to be defined differently for each spatial resolution level.The stands were defined by 30 m circular buffer zones whose centroids were localized by GPS measurements during the 1998 campaign (see Section 2.4).The pixels representing sunlit parts of Norway spruce crowns defined by the Mahalanobis distance classificator were taken into account in the case of the original resolution data (Figure 2), whereas all the pixels within the defined buffers were used in the case of the generalized level data.
( ) where: Br, Gr, and Wr are rescaled Brightness, Greenness, and Wetness components; B, G, and W are the original Brightness, Greenness, and Wetness values; Bµ, Gµ, and Wµ are the means of B, G, and W calculated from the forest land areas; Bσ, Gσ, Wσ are standard deviations of B, G, and W calculated from the forest land areas defined by the appropriate binary mask.
The concept of DI assumes that high Br and low Gr and Wr values are typical for disturbed stands.The DI value is high (positive) in such a case.On the other hand, undisturbed, fully-regenerated and young vital stands are typical for low Br and high Gr and Wr values resulting in low (negative) DI values [44].
The Landsat 7 scenes were acquired in SLC-off mode after 2003 due to a malfunction of the Scan Line Corrector of the ETM+ sensor (compensating the orbital motion of the satellite).This results in stripes of zero-value pixels in the imagery.Thus, these areas were masked out and only the valid values were used for the all analysis.
A raster layer of the DI values calculated for the forest land areas was constructed for each Landsat scene.The DI values were then classified into three classes using mean (µ) and standard deviation (σ) values.Class 1 includes low DI values (less than µ − 0.5σ) representing undisturbed and most vital forest, whereas Class 3 includes high DI values (more than µ + 0.5σ) representing the most disturbed stands.Class 2 represents the most common DI values between µ − 0.5σ and µ + 0.5σ.

Vegetation Indices Derived from Hyperspectral Image Data and Stand Separability
The forest health status analysis was performed at two levels of spatial resolution of the hyperspectral imagery: original (2.0 m) considering individual tree crowns and generalized (6.0 m) considering whole forest stands.Therefore, the stands had to be defined differently for each spatial resolution level.The stands were defined by 30 m circular buffer zones whose centroids were localized by GPS measurements during the 1998 campaign (see Section 2.4).The pixels representing sunlit parts of Norway spruce crowns defined by the Mahalanobis distance classificator were taken into account in the case of the original resolution data (Figure 2), whereas all the pixels within the defined buffers were used in the case of the generalized level data.In the case of the ASAS dataset, only the bands between 550-850 nm were usable for further processing.This limited the number of vegetation indices (VIs) that could be calculated from these data (e.g., all the indices requiring bands in the VIS-B region had to be excluded).The APEX data were resampled to the spectral resolution of the ASAS data prior to calculating the VIs to ensure their mutual comparability.The vegetation indices used within this study are listed in Table 3 [17,[46][47][48][49][50][51][52][53][54][55][56].In the case of the ASAS dataset, only the bands between 550-850 nm were usable for further processing.This limited the number of vegetation indices (VIs) that could be calculated from these data (e.g., all the indices requiring bands in the VIS-B region had to be excluded).The APEX data were resampled to the spectral resolution of the ASAS data prior to calculating the VIs to ensure their mutual comparability.The vegetation indices used within this study are listed in Table 3 [17,[46][47][48][49][50][51][52][53][54][55][56].
To check their sensitivity to forest damage, a simple separability index (SI) was calculated between the initially damaged and moderately to heavily damaged forest stands for all the VIs.The SI was calculated using the following formula used in [45]: where: SI is the separability index, µ h and σ h are the mean and standard deviation calculated from the pixels within initially damaged stands, µ d and σ d are the mean and standard deviation calculated from pixels within the moderately to heavily damaged forest stands.The original values of the given VIs obtained from the ASAS dataset are not directly comparable with the ones extracted from the APEX dataset.This is mainly due to technological differences between ASAS/APEX sensors and differences in atmospheric correction applied to the ASAS and APEX data etc.The original VI values were, therefore, normalized to allow their mutual comparability using the following formula: where: VI 1 x is the normalized value of vegetation index VI for pixel x, VI x is the original value of vegetation index VI for pixel x, VI µ and VI σ are the mean and standard deviation values of the vegetation index VI calculated from the all pixels representing sunlit crowns within the buffer areas representing the selected forest stands.

Detection of Spatio-Temporal Differences
The values of the VIs used (or VI 1 s) were extracted from both ASAS and APEX data using the buffer zones (defining the particular stands) and tree crowns mask (in case of the original resolution data).The Student's two-sample test was used to analyse spatial differences as well as temporal changes in the used VI 1 s values (considered as the markers of forest physiological status).In the case of the spatial differences , the VI values were grouped regarding the classification of the stands as moderately to heavily damaged or initially damaged (1998), and formerly moderately to heavily damaged or still initially damaged (2013).The test was applied to analyze whether there are any significant differences in VI values between these groups (separately for 1998 and 2013).
The normalized VI values describe the relative distance of the given pixel from the global mean (calculated from all Norway spruce pixels within the defined stands).Thus, the analysis was whether any significant change of the relative distance to the global mean occurred at each given stand.Both tests were performed at the 95% significance level (α = 0.05).

Statistical Assessment Using Field Reference Data
Two non-specific stress markers based on photosynthetic pigments ratios were used as the indicators of tree physiological status: chlorophyll-a to chlorophyll-b ratio (C a /C b ) and total carotenoids to total chlorophylls ratio (C x /C ab ).Worsening of tree physiological status is usually exhibited by an increase in the C a /C b ratio as chlorophyll-b is regarded to be more sensitive to vegetation stress than chlorophyll-a [57].Similarly, stress conditions lead to an increase in the C x /C ab ratio [58][59][60].
Temporal differences in photosynthetic pigments contents and their ratios were evaluated for the stands sampled in both years (1998 and 2013) and covered by both hyperspectral image datasets.A two-way ANOVA (performed at the 95% significance level α = 0.05) was used for evaluating the differences between the western and central Ore Mountains sites and temporal changes.The differences between the sites in particular years were determined by the Tukey-Kramer multiple comparison test.

General Trends of Forest Recovery
Temporal changes of DI values were studied across the time series of the Landsat data for both western and central Ore Mountains sites.Figure 3 shows examples of the DI raster layers.The initial situation in the second half of the 1980s is represented by the scene acquired in 1985.The scene from 1998 represents the phase of the highest difference between the central and western Ore Mountains, whereas the scene acquired in 2015 shows the current situation.Figure 4 shows the relative frequencies of the three defined DI classes and changes of the DI values at the forest stands where a health status assessment and needle sampling were performed during the 1998 and 2013 campaigns.
It can be seen that there is significant difference between the two study areas.The central Ore Mountains exhibited higher DI values indicating a higher disturbance level than the forests in the western part during the 1980s and 1990s with peak in 1998, while in the western Ore Mountains the forests disturbance level did not change much during the same period.It can also be seen that the relative frequency of Class 3 has been declining slightly since around the year 2000, representing regeneration processes in the central Ore Mountains area.
In addition, mean DI values were calculated for the Norway spruce stands sampled during both field campaigns in 1998 and 2013 (i.e., P 9 , P 10 , P 11 , P 40 , K 25 , K 37 , K 46 , and K 47 )-see Figure 4. Again, the higher DI of the forest stands in the central Ore Mts indicate a worse health status compared to the ones at the western Ore Mountains.The highest difference between the central and western Ore Mountains sites can be seen at the end of the 1990s.After that a levelling trend of the two sites can be observed.

Assessment of Forest Health Status Change Using ASAS and APEX Datasets
The sensitivity analysis proved that the majority of the VIs used is highly sensitive to vegetation health status and forest damage.Two sets of five indices were defined (one set for the original 2 m spatial resolution data and one for the generalized 6 m spatial resolution) providing the best results to distinguish initially and moderately to heavily damaged regions.The fist set designed for the original

Assessment of Forest Health Status Change Using ASAS and APEX Datasets
The sensitivity analysis proved that the majority of the VIs used is highly sensitive to vegetation health status and forest damage.Two sets of five indices were defined (one set for the original 2 m spatial resolution data and one for the generalized 6 m spatial resolution) providing the best results to distinguish initially and moderately to heavily damaged regions.The fist set designed for the original (2 m) resolution data included the following indices: VOG 1 , NDVI, TCARI/OSAVI, MSR, and NDVI 705 , whereas the set designed for the generalized (6 m) resolution included the following indices: N 714 , VOG 1 , TCARI/OSAVI, MSR and NDVI 705 .The values of these most sensitive vegetation indices together with the calculated separability scores are shown in Figure 5.
Significant differences were detected for the VI values between the moderately to heavily damaged stands (located at the central Ore Mountains) and initially damaged ones (located at the western Ore Mountains) in both the original (2 m) and generalized (6 m) spatial resolution of the ASAS (1998) data.The p-values of the tests carried out were far less than 0.01 in all cases.Moreover, these differences were also detectable in the case of the APEX (2013) data, although they were not as noticeable as in the case of the ASAS dataset.
Remote Sens. 2016, 8, 92 12 these differences were also detectable in the case of the APEX (2013) data, although they were not as noticeable as in the case of the ASAS dataset.The temporal changes were studied primarily using the VOG1 index, as it was showing high and stable sensitivity to forest damage for both spatial resolutions considered.Significant changes to the normalized VI values were proven for a major part of the stands studied except the P5 (p-values 0.28) and P17 (p-value 0.46).Examples of the normalized VOG1 index values for two selected stands can be seen in Figure 6.
A simple visualization was developed to facilitate interpretation of the results obtained (Figure 7).The charts show the relative distance of the particular stands (defined by the local mean of the given VI') from the baselines (defined by the global means of the given VI') in both time horizons.Four possible situations can be defined from this point of view: • positive stagnation (+/+): local mean of the normalized VI' values was positive in both time horizons.The given stand was above the global mean in both years.
• negative stagnation (−/−): local mean of the normalized VI' values was negative in both years.
The stand was below the global mean in both years.The temporal changes were studied primarily using the VOG 1 index, as it was showing high and stable sensitivity to forest damage for both spatial resolutions considered.Significant changes to the normalized VI values were proven for a major part of the stands studied except the P 5 (p-values 0.28) and P 17 (p-value 0.46).Examples of the normalized VOG 1 index values for two selected stands can be seen in Figure 6.
A simple visualization was developed to facilitate interpretation of the results obtained (Figure 7).The charts show the relative distance of the particular stands (defined by the local mean of the given VI') from the baselines (defined by the global means of the given VI') in both time horizons.Four possible situations can be defined from this point of view:  The position of the sites in Figure 7 describes the temporal change at the particular stands using the normalized VOG1 vegetation index.For the stands in the central Ore Mountains, the site K33 exhibited negative stagnation while site K37 tended to recover.Unfortunately, only two sites could have been included in the analyses at the central Ore Mountains because of the limited spatial coverage of the ASAS dataset in this site.The sites in the western Ore Mountains exhibited mainly positive stagnation and some of them negative stagnation or worsening.For this site, it was possible to include 14 stands in the analysis.The position of the sites in Figure 7 describes the temporal change at the particular stands using the normalized VOG 1 vegetation index.For the stands in the central Ore Mountains, the site K 33 exhibited negative stagnation while site K 37 tended to recover.Unfortunately, only two sites could have been included in the analyses at the central Ore Mountains because of the limited spatial coverage of the ASAS dataset in this site.The sites in the western Ore Mountains exhibited mainly positive stagnation and some of them negative stagnation or worsening.For this site, it was possible to include 14 stands in the analysis.The position of the sites in Figure 7 describes the temporal change at the particular stands using the normalized VOG1 vegetation index.For the stands in the central Ore Mountains, the site K33 exhibited negative stagnation while site K37 tended to recover.Unfortunately, only two sites could have been included in the analyses at the central Ore Mountains because of the limited spatial coverage of the ASAS dataset in this site.The sites in the western Ore Mountains exhibited mainly positive stagnation and some of them negative stagnation or worsening.For this site, it was possible to include 14 stands in the analysis.

Differences in Photosynthetic Pigments Content
The descriptive statistics of the studied Norway spruce needles parameters are listed in the Table 4.The statistical analysis did not prove any significant differences in biochemically determined total chlorophyll content in Norway spruce needles samples at either site in 1998 and 2013, if the site (p = 0.66) and year (p = 0.49) were used as fixed factors in the two-way ANOVA (Table 5).Nevertheless, differences were detected in the ratios of photosynthetic pigments (C a /C b and C x /C ab ).The C a /C b ratio was significantly affected by year (p less than 0.001), but not by site (p = 0.32).Similarly the C x /C ab ratio was high at the central Ore Mountains in 1998 in comparison to the western Ore Mountains in 1998 and both localities in 2013 (year effect: p = 0.04; locality effect: p = 0.01).According to the Tukey-Kramer multiple comparison test, the C a /C b ratio at central Ore Mountains in 1998 was higher in comparison to western Ore Mountains in 1998.Both pigments ratios (C a /C b and C x /C ab ) for both sites in 2013 decreased in comparison to values shown by central Ore Mountains in 1998.This implies that the trees at the central Ore Mountains recovered their photosynthetical apparatus, whilst trees in the western Ore Mountains did not show a significant temporal change in the ratio of photosynthetical pigments.This implies that the physiological status of needles of the first three needle age classes remained unchanged in the western Ore Mountains, while in the central Ore Mountains it had significantly improved by 2013 compared to 1998.

Discussion
The ASAS airborne hyperspectral data were acquired in 1998, the year of the maximum forest disturbance level.Significant differences between the moderately to heavily damaged (central Ore Mountains) and initially damaged (western Ore Mountains) stands were detected for the all tested VIs calculated from the ASAS dataset.The normalized VI' values were far below the global mean in the case of the central Ore Mountains stands, whereas for the western Ore Mountains they were closely distributed around the global mean value.This result fully confirms the finding of the Landsat time series analysis.The highest difference in DI values between the central and western Ore Mountains stands could be seen at the end of the 1990s.Afterwards, however, a levelling trend of the two sites can be observed.
In 2013, the stands at the central Ore Mountains exhibited VI's values much closer to the global mean (compared to 1998) indicating an improvement in the health status of the local Norway spruce forests.The VI's values for the stands in the western Ore Mountains indicate stagnation or, in two cases (P 1 and P 6 ), a slight worsening of tree physiological status.Nevertheless, there are still minor differences between the stands in the western and central Ore Mountains, which are detectable by the majority of the VIs used.This demonstrates the ability of the chosen VIs to detect not only a high level of forest damage, but also slight differences in the physiological status of trees, which did not exhibit a high level of defoliation (on average 30% in the central Ore Mountains and 33% in the western Ore Mountains).
The concentrations of photosynthetic pigments themselves did not prove to be a sensitive indicator of the change in tree physiological status.Ratios of photosynthetic pigments are known to be sensitive to air pollution [58,61].The differences in C a /C b and C x /C ab ratios generally corresponded with the spatial gradient of the forest physiological status observed by remote sensing tools in both comparison years 1998 and 2013.In 2013, the analysis confirmed an improvement in photosynthetic pigments ratios in the case of the central Ore Mountains, whereas the values for the western Ore Mountains site remained the same.
The results obtained from non-specific stress markers based on the ratios of photosynthetic pigments and VIs analyses led to the assumption that the VIs values are a very powerful tool for determining of slight changes in the vegetation's physiological status.VIs seem to be the same or even more sensitive than biochemically-derived ratios of photosynthetical pigments.It is assumed that VIs can be also influenced by other biophysical factors which are not directly linked to foliar biochemistry such as needle orientation on a shoot, needle retention, primary or secondary shoots, etc. Crown architecture changes including formation of secondary shoots can provide information of retrospective tree response to the impact of stress factors [62,63].As such, this may have an influence on the forest stand spectral signature recorded by space or airborne sensors.
Noticeable peaks of DI values derived from the Landsat time series can be observed in the case of the central Ore Mountains in 1998 as well as in 2000 in the case of the western Ore Mountains area.Nevertheless, it is unclear whether they are caused by a delayed vegetation stress reaction, climate change events or by any issue related to the data used and their processing.However, the maximum extent of the area covered by highly disturbed forest was reached between 1998 and 2000.This can even be seen by just a visual analysis of the Landsat imagery.
One problematic aspect seems to be the use of TOA reflectance for the TCT time series analysis.On the one hand, TOA reflectance is recommended for the TCT rather than at-surface reflectance by [41,42].On the other hand, this means that atmospheric effects were not removed from the data, which gives rise to the question of comparability of the TCT components values obtained from the particular scenes.Another questionable issue is the comparability of the TCT components values extracted from Landsat 5/7 TM/ETM+ sensors with the ones obtained from Landsat 8 OLI imagery, as the band settings are similar but not the same.
The improvement in forest physiological status after 1998, which was detected by analyzing Landsat time series, as well as ASAS and APEX datasets, is in good agreement with the desulfurization of all the coal power stations in Northern Bohemia conducted between 1996 and 1999 [64].As such, the levels of fly ash decreased by 95% and the SO 2 and NO x emissions decreased by 92%, respectively; 50% compared to the situation at the beginning of 1990s [64].The area covered by dead forests reached its maximal extent in the 1998-2000 period [65].
Despite the dramatic drop in emissions, soil acidification and nitrification still remains a serious problem.The soil pH was reported generally less than 4.0 with an extreme value of 2.2 in the Ore Mountains region in 1989 [66].However, the critical values were also exceeded at a much later time [67].For this reason it seems that the full regeneration of forest stands in the formerly-disturbed areas of the Ore Mountains region will be a long-term process.

Conclusions
Regarding spatio-temporal changes in forest status, in 1998 significant differences were observed between moderately to heavily damaged (central Ore Mountains) and initially damaged (western Ore Mountains) stands; whilst in 2013 a levelling of the physiological status of Norway spruce stands has been observed for both sites.The stands in the western Ore Mountains exhibited stagnation or in two cases a slight worsening of the physiological status, while the Norway spruce forests located in the central Ore Mountains have significantly improved during the period under study.This trend was confirmed by ground reference data and comparison of the VIs extracted from ASAS and APEX airborne hyperspectral imagery.The most sensitive hyperspectral vegetation indices used in this study (VOG 1 , TCARI/OSAVI, MSR, NDVI, and NDVI 705 ) demonstrated their ability to detect not only severe forest damage (as in the case of the ASAS dataset), but also slight differences in the physiological status of trees, which did not exhibit high defoliation levels (in average 30% in the central Ore Mountains and 33% in the western Ore Mountains).
The improvement of the forest physiological status observed at the central Ore Mountains corresponds well with the desulfurization of the coal power stations in the Northern Bohemia region conducted between 1996 and 1999, leading to a tremendous reduction in sulfur air pollution load.
Due to the limited spatial extent of the hyperspectral imagery acquired in 1998, the study is a methodical demonstration showing the necessary steps regarding the diverse sensor data normalization and consequent spatio-temporal analysis.Nevertheless, the results obtained seem to be promising and demonstrate a great potential for further use of the proposed workflow.A higher temporal resolution of hyperspectral data (e.g., one acquisition per each vegetation season) would positively improve knowledge about forest recovery processes.The main environmental drivers affecting forest productivity and services are currently regarded to be climate change together with anthropogenic pressure including air pollution [68,69].Since forests represent a key component in carbon cycling, the effective large-scale methods of remote sensing to detect their physiological status and the capacity of carbon sequestration are becoming crucially important.This study demonstrates the capability of highly sensitive VIs derived from hyperspectral data combined with the long-term Landsat series analysis to facilitate this demand for remote sensing-based monitoring techniques.
It was shown that the processing used is only slightly affected by the spatial scale as the results of the spatio-temporal analysis performed at crown (2 m spatial resolution) and stand (6 m spatial resolution) are very similar.Moreover, the vast majority of the vegetation indices which were found to be the most appropriate for studying forest physiological status changes could be calculated not only from hyperspectral data (e.g., airborne HS imagery or up-coming EnMap satellite HS sensor), but also from superspectral satellite imagery as providing the appropriate spectral bands.It is, therefore, highly probable that the same data processing approach might be applied to the data provided by new satellite sensors like Sentinel-2 MSI (with 10-20 m spatial resolution) with a comparable performance.In this case, the relatively narrow bands 5 and 6, localized in the red-edge region (705 and 740 nm), seems to be the most suitable.The use of such satellite data will provide opportunity to study changes in forest physiological status on a far larger spatial extent with much higher temporal resolution, which are the two most limiting parameters when using airborne imagery.

Figure 1 .
Figure 1.(A) Main sources of air pollution (coal power stations, heavy industrial facilities, etc.) in the Ore Mountains region; (B) Detailed maps of the western (WE) and central (CE) Ore Mountains sites with positions of the Norway spruce stands sampled during the field campaigns in 1998 and 2013.Background map provided by ČUZK.

Figure 1 .
Figure 1.(A) Main sources of air pollution (coal power stations, heavy industrial facilities, etc.) in the Ore Mountains region; (B) Detailed maps of the western (WE) and central (CE) Ore Mountains sites with positions of the Norway spruce stands sampled during the field campaigns in 1998 and 2013.Background map provided by ČUZK.

Figure 2 .
Figure 2. Definition of the Norway spruce sampling stands in the original APEX data (A) and binary mask of the sunlit tree crowns (B).

Figure 2 .
Figure 2. Definition of the Norway spruce sampling stands in the original APEX data (A) and binary mask of the sunlit tree crowns (B).

R
λ stands for reflectance at wavelength λ D λ stands for 1st derivation of reflectance at wavelength λ

Figure 3 .
Figure 3. Disturbance index (DI) values for coniferous forests calculated from Tasseled Cap-transformed Landsat data for western (WE) and central (CE) Ore Mountains localities for 1985, 1998, and 2015.The lower DI values (shown in green) correspond to relatively less disturbed stands, whereas high DI values (shown in red) correspond to relatively more disturbed stands.

Figure 3 .
Figure 3. Disturbance index (DI) values for coniferous forests calculated from Tasseled Cap-transformed Landsat data for western (WE) and central (CE) Ore Mountains localities for 1985, 1998, and 2015.The lower DI values (shown in green) correspond to relatively less disturbed stands, whereas high DI values (shown in red) correspond to relatively more disturbed stands.

Figure 4 .
Figure 4. Relative frequency of the Disturbance Index (DI) classes calculated for the forest areas at the western (WE) and central (CE) Ore Mountains Class 1-less disturbed areas, Class 2-most common disturbance level, and Class 3-more disturbed areas (a and b).Mean DI values calculated for the Norway spruce stands sampled during the field campaigns in 1998 and 2013 (c).

Figure 4 .
Figure 4. Relative frequency of the Disturbance Index (DI) classes calculated for the forest areas at the western (WE) and central (CE) Ore Mountains Class 1-less disturbed areas, Class 2-most common disturbance level, and Class 3-more disturbed areas (A and B).Mean DI values calculated for the Norway spruce stands sampled during the field campaigns in 1998 and 2013 (C).

Figure 5 .
Figure 5. Values of the most sensitive vegetation indices calculated from original (2 m) and generalized (6 m) spatial resolution ASAS imagery from 1998.ID stands for initially damaged stands (in total 14 sites), whereas HD stands for moderately to heavily damaged forest stands (in total two sites).SI represents separability scores.

Figure 5 .
Figure 5. Values of the most sensitive vegetation indices calculated from original (2 m) and generalized (6 m) spatial resolution ASAS imagery from 1998.ID stands for initially damaged stands (in total 14 sites), whereas HD stands for moderately to heavily damaged forest stands (in total two sites).SI represents separability scores.

‚
positive stagnation (+/+): local mean of the normalized VI' values was positive in both time horizons.The given stand was above the global mean in both years.‚ negative stagnation (´/´): local mean of the normalized VI' values was negative in both years.The stand was below the global mean in both years.‚ recovery (´/+): local mean of the normalized VI' values was negative in 1998, but positive in 2013.The stand was below the global mean in 1998, but above the global mean in 2013.‚ worsening (+/´): the given stand was above the global mean in 1998, but below the global mean in 2013.• recovery (−/+): local mean of the normalized VI' values was negative in 1998, but positive in 2013.The stand was below the global mean in 1998, but above the global mean in 2013.• worsening (+/−): the given stand was above the global mean in 1998, but below the global mean in 2013.

Figure 6 .
Figure 6.Normalized values of the VOG1 vegetation index calculated from the ASAS (1998) and APEX (2013) hyperspectral imagery datasets in the original 2 m spatial resolution for P11 (western Ore Mountains-less damaged) and K37 (central Ore Mountains-more damaged) stands.The lower values (shown in red colour) correspond to relatively worse physiological status, while high values (shown in green) correspond to a relatively good physiological status.

Figure 6 .
Figure 6.Normalized values of the VOG 1 vegetation index calculated from the ASAS (1998) and APEX (2013) hyperspectral imagery datasets in the original 2 m spatial resolution for P 11 (western Ore Mountains-less damaged) and K 37 (central Ore Mountains-more damaged) stands.The lower values (shown in red colour) correspond to relatively worse physiological status, while high values (shown in green) correspond to a relatively good physiological status.
+): local mean of the normalized VI' values was negative in 1998, but positive in 2013.The stand was below the global mean in 1998, but above the global mean in 2013.• worsening (+/−): the given stand was above the global mean in 1998, but below the global mean in 2013.

Figure 6 .
Figure 6.Normalized values of the VOG1 vegetation index calculated from the ASAS (1998) and APEX (2013) hyperspectral imagery datasets in the original 2 m spatial resolution for P11 (western Ore Mountains-less damaged) and K37 (central Ore Mountains-more damaged) stands.The lower values (shown in red colour) correspond to relatively worse physiological status, while high values (shown in green) correspond to a relatively good physiological status.

Figure 7 .
Figure 7. Visualization of Norway spruce health status temporal change using normalized VOG 1 vegetation index constructed for original (A) and generalized (B) spatial resolution of the used image data.P-western Ore Mountains (Přebuz) sites, and K-central Ore Mountains (Kovářská) sites.

Table 2 .
[17,38]health status classification scheme developed during the sampling campaign in 1998.Adopted from[17,38].The same scheme was used during the campaign in 2013.

Table 3 .
List of vegetation indices used in this study.

Table 4 .
Mean ˘standard deviation of the chlorophyll a + b (mg¨g ´1), pigments ratios, and relative water content (%) in the Norway spruce needles collected at the western and central Ore Mountains sites in 1998 and 2013.Different letters (a, b) indicate significant difference between year and site combinations according to the performed Tukey-Kramer test (α = 0.05).

Table 5 .
Effect of year sampling, site and their interaction on content of photosynthetic pigments, their ratios and water content.Two-way ANOVA, n.s.= not significant, * significant with p < 0.05.