Exploring Digital Surface Models from Nine Different Sensors for Forest Monitoring and Change Detection

Digital surface models (DSMs) derived from spaceborne and airborne sensors enable the monitoring of the vertical structures for forests in large areas. Nevertheless, due to the lack of an objective performance assessment for this task, it is difficult to select the most appropriate data source for DSM generation. In order to fill this gap, this paper performs change detection analysis including forest decrease and tree growth. The accuracy of the DSMs is evaluated by comparison with measured tree heights from inventory plots (field data). In addition, the DSMs are compared with LiDAR data to perform a pixel-wise quality assessment. DSMs from four different satellite stereo sensors (ALOS/PRISM, Cartosat-1, RapidEye and WorldView-2), one satellite InSAR sensor (TanDEM-X), two aerial stereo camera systems (HRSC and UltraCam) and two airborne laser scanning datasets with different point densities are adopted for the comparison. The case study is a complex central European temperate forest close to Traunstein in Bavaria, Germany. As a major experimental result, the quality of the DSM is found to be robust to variations in image resolution, especially when the forest density is high. The forest decrease results confirm that besides aerial photogrammetry data, very high resolution satellite data, such as WorldView-2, can deliver results with comparable quality as the ones derived from LiDAR, followed by TanDEM-X and Cartosat DSMs. The quality of the DSMs derived from ALOS and Rapid-Eye data is lower, but the main changes are still correctly highlighted. Moreover, the vertical tree growth and their relationship with tree height are analyzed. The major tree height in the study site is between 15 and 30 m and the periodic annual increments (PAIs) are in the range of 0.30–0.50 m.


Introduction
Height is one of the most important attributes to characterize forest stands.Besides being indispensable for estimating forest timber volume, it is very helpful for forest structure analyses, classification, mapping and change detection [1].New satellite sensors and advanced image processing techniques provide data and tools for an improved forest monitoring.The question is: can forest height and timber volume estimation be carried out efficiently using remote sensing data?To answer this question, an extensive understanding about the quality and potential applicability of height information derived from the latest sensors and processing technologies is necessary.
Although forest height is currently monitored mainly through forest inventory plots, Digital Surface Models (DSMs) and Canopy Height Models (CHMs) from remote sensing data have recently gained more attention.In general, four main remote sensing techniques are available to acquire three-dimensional data in forests: LiDAR (Light Detection And Ranging), aerial photogrammetry, satellite photogrammetry and Interferometric Synthetic Aperture Radar (InSAR) [2][3][4].
In the last 20 years, different airborne LiDAR (also known as airborne laser scanning, ALS) systems have been used to measure canopy height [1,[5][6][7][8][9][10].Full-waveform LiDAR systems provide an insight into forest canopies, delivering even detailed information on the vertical structure of forests [11,12].Recent research presented solutions for single tree identification even if crowns are interlaced [13,14].Moreover, LiDAR systems are able to penetrate the forest canopy, capturing additional information about the vegetation structure below the top forest layer.
Besides LiDAR, photogrammetry uses overlapping images taken from different viewing-angles to derive height information [15].Aerial photogrammetry is not a new technology, but improvement in forest monitoring applications has been relatively slow due to limits in image quality and image processing capabilities.Height information from LiDAR is more direct and more accurate than data derived from photogrammetry [16], and it is presently considered as reference for forest height determination.However, the fast developments in computer vision methodology in recent years, especially the development of dense image matching algorithms, have largely improved the quality of DSMs generated from photogrammetry, thus boosting the research in aerial photogrammetry for forest structure monitoring [3,15,[17][18][19][20][21].The costs for acquiring aerial imagery are much lower than the costs for LiDAR data [22] and image data acquisitions are performed regularly for all German federal states.Additionally, multispectral aerial photographs allow identifying tree species groups [15,21].
Satellite photogrammetry has also been assessed for forest attribute estimation [3,4,23], since it can be performed more efficiently and at a lower cost with respect to aerial photogrammetry.In addition to a larger field of view, satellite images often provide spectral information in the visible and Near Infrared (NIR) ranges (e.g., WorldView-2 satellite data).In general, two main acquisition types for satellite sensors are available to derive a DSM: across-track and along-track.Across-track stereoscopy captures data from two or more adjacent orbits, which means that changes in surface reflectance between the two orbits can impair the quality of the generated DSM.Therefore, since its introduction by the experimental MOMS-02 system in the early 1990s [24], this technology is now used by the majority of recent satellites, for instance SPOT-5, IKONOS, GeoEye-1/2, WorldView-2/3, QuickBird, ALOS/PRISM, Cartosat-1, and Ziyuan-3 [25].These satellites obtain stereo imagery from the same orbit from two or more views, either by using multiple cameras or by adjusting the viewing angles [25].The spatial resolution of the satellite stereo imagery has also improved steadily in recent years, with the recently launched WorldView-3 and WorldView-4 satellites able to capture stereo imagery with up to 0.31 m resolution.With the availability of these data, research on space-borne photogrammetry is gaining more attention, including for forest monitoring [4].
Elevation data can also be provided through SAR techniques.Using radargrammetry, height information is obtained from two or more multi-view SAR images, and Interferometric SAR (InSAR) derives height information through phase differences.Recent advances in very high resolution SAR imagery and the launch of TerraSAR-X and TanDEM-X [26] have enabled the use of this technique for forest biomass monitoring [27,28].
Tree height, which is closely related to tree age (young stands) and volume (especially old stands), is an important parameter in forest management.At regional to enterprise level, remote sensing data are capable of quickly locating and estimating damages after a disturbance event, especially in the case of a sudden storm.Forest monitoring at national/global level or after a sudden event usually employs data from various sensors.Therefore an evaluation of the data quality for the different sensors, their derived DSMs and behavior in change detection is of high importance.
The primary objectives of this paper are to introduce the most important sensors which can provide DSMs, including optical stereo data, InSAR systems, and airborne LiDAR, and assess the quality of the extracted height information through investigating their ability for forest change monitoring.As a case study, we chose a highly-structured forest area near Traunstein in Bavaria, Germany, since sufficient data from various sensors including reference data are available for this region.
The remainder of this paper is organized as follows.In Section 2, the test region and the characteristics of the imagery presented and DSM data are introduced.The change detection approaches, including forest damage detection and forest annual growth detection are presented in Section 3. As the accuracy of the generated DSMs directly influence the height-oriented change detection, a quality evaluation of these DSMs is necessary before performing the change detection experiments.The DSM quality is investigated at pixel-and inventory plot-level.The evaluation approaches and experimental results are described in Section 4 as Experiment I.The forest change detection results are grouped as Experiment II and presented in Section 5.In Section 6, a discussion on the datasets and experimental results is reported.Finally, the conclusions are presented in Section 7.

Study Site
The study site "Stadtwald Traunstein" is located in southeastern Bavaria, close to Lake Chiemsee and the Austrian border (Figure 1).It is a long-term research site of the Technical University of Munich (TUM) and is currently managed by the Chair of Forest Growth and Yield Science at TUM.The forest covers an area of about 242 ha.
Remote Sens. 2017, 9, 287 3 of 26 The primary objectives of this paper are to introduce the most important sensors which can provide DSMs, including optical stereo data, InSAR systems, and airborne LiDAR, and assess the quality of the extracted height information through investigating their ability for forest change monitoring.As a case study, we chose a highly-structured forest area near Traunstein in Bavaria, Germany, since sufficient data from various sensors including reference data are available for this region.
The remainder of this paper is organized as follows.In Section 2, the test region and the characteristics of the imagery presented and DSM data are introduced.The change detection approaches, including forest damage detection and forest annual growth detection are presented in Section 3. As the accuracy of the generated DSMs directly influence the height-oriented change detection, a quality evaluation of these DSMs is necessary before performing the change detection experiments.The DSM quality is investigated at pixel-and inventory plot-level.The evaluation approaches and experimental results are described in Section 4 as Experiment I.The forest change detection results are grouped as Experiment II and presented in Section 5.In Section 6, a discussion on the datasets and experimental results is reported.Finally, the conclusions are presented in Section 7.

Study Site
The study site "Stadtwald Traunstein" is located in southeastern Bavaria, close to Lake Chiemsee and the Austrian border (Figure 1).It is a long-term research site of the Technical University of Munich (TUM) and is currently managed by the Chair of Forest Growth and Yield Science at TUM.The forest covers an area of about 242 ha.For more than 30 years this study site has been managed with the aim of increasing the forest stability against biotic and abiotic calamities.This results in highly structured, mixed forest stands.Presently, the forest is composed of a mixture of spruce (Picea abies), including aged pure spruce stands, fir (Abies alba), Douglas fir (Pseudotsuga menziesii), beech (Fagus sylvatica), sycamore maple (Acer pseudoplatanus), ash (Fraxinus excelsior), and others.As Traunstein forest is a managed forest, all For more than 30 years this study site has been managed with the aim of increasing the forest stability against biotic and abiotic calamities.This results in highly structured, mixed forest stands.Presently, the forest is composed of a mixture of spruce (Picea abies), including aged pure spruce stands, fir (Abies alba), Douglas fir (Pseudotsuga menziesii), beech (Fagus sylvatica), sycamore maple (Acer pseudoplatanus), ash (Fraxinus excelsior), and others.As Traunstein forest is a managed forest, all changes that took place in the last 100 years are well documented.The main sources of changes are small local logging activities, tree growth and exceptional events such as storm damages.On 18 January 2007, storm Kyrill caused widespread damages in the area, especially in the aged pure spruce stands [29].In Traunstein forest, trees are harvested applying a selective logging system, based on local single tree cutting.

Remote Sensing Data
The remote sensing datasets used in this paper were collected in different years for different purposes, and are summarized in Table 1.The main characteristics and technical parameters of each sensor are reported.The data are captured in a time span ranging from 2001 to 2013 and include four satellite stereo sensors-(ALOS/PRISM, Cartosat-1, RapidEye and WorldView-2), one satellite InSAR sensor (TanDEM-X) and two aerial stereo camera systems (HRSC and UltraCam), as well as two LiDAR data sets.Most of the nine DSMs displayed in Figure 2 have been derived using data acquired during leaf-on season.
• Airborne LiDAR data LiDAR data acquired in 2001 and 2010 with different point density were provided by the Bavarian Administration for Surverying (Bayerische Landesamt für Digitalierung, Breitband und Vermessung, LDBV).The LiDAR data from 2001 were acquired in February with a point density of around 0.5 points/m 2 .The 2010 dataset was acquired from March to April 2010, and exhibits a much higher density of 5.65 points/m 2 .The DSM and DTM used in this paper were prepared by using the filtering method implemented in the software TreesVis [30], and are also described in [4].In the following text, the generated LiDAR DSMs are referred to as LiDAR2001 and LiDAR2010, respectively (see Figure 2a,b).

• Satellite optical data
Data from four optical satellite systems are involved in this research, including ALOS/PRISM, Cartosat-1, RapidEye and WorldView-2.The Panchromatic Remote-Sensing Instrument for Stereo Mapping (PRISM) is one of the three instruments on board of the Japanese satellite ALOS (Advanced Land Observation Satellite).ALOS was launched on 24 January 2006 and stopped acquiring data on 21 April 2011.PRISM provided 2.5 m resolution stereo imagery with up to 70 km swath width and a revisit time of 2 days [25,31].The ALOS data used in this research was captured in June 2008.Cartosat-1 was launched on 5 May 2005 by the Indian Space Research Organization (ISRO).It carries two panchromatic cameras that obtain stereo imagery with a forward view of +26 • (Band F) and a nadir view of −5 • (Band A).The panchromatic images obtained have a spatial resolution of 2.5 m and cover a swath of 30 km.The satellite has a maximum revisit time of five days [25].Three pairs of Cartosat data captured in February, July and November 2008 were available for this study.The RapidEye satellite system, launched at the end of 2008, provides data in five bands, namely blue, green, red, red edge and NIR, with a spatial resolution of 6.5 m.The RapidEye sensor is not designed to provide stereo data but has a very frequent revisit interval of 2-3 days, allowing across-track stereo matching [32].The RapidEye data captured in 17 and 23 May 2009 are used as stereo pair in the DSM generation procedure.The acquisition times of these two datasets are very close, which leads to similar spectral information of the tree leaves.WorldView-2 is a commercial Earth observation satellite that provides very high resolution space borne images.Being Digital Globe's third satellite, it was launched on 8 October 2009.It provides 0.5 m panchromatic images and 2 m resolution eight band multi-spectral images, adding to the usual four bands (red, green, blue, and NIR1) a set of four new bands (coastal, yellow, red edge and NIR2).WorldView-2 is a pushbroom sensor able to collect large areas of stereo-view or multi-view imagery in a single pass by adjusting the viewing angle of the camera by steering the satellite.The average revisit time of the satellite can be up to 1.1 days by using also off-nadir viewing [25].
For the satellite stereo sensors, DSMs were generated from ALOS/PRISM (Figure 2d), Cartosat-1 (Figure 2e), RapidEye (Figure 2f) and WorldView-2 (Figure 2g) stereo imagery using the Semi Global Matching (SGM) algorithm implemented at DLR [33].The DSM generation and gap filling procedures are described in detail in [34] and [4].It has to be mentioned that the DSM from Cartosat-1 is a fused result from the three available pairs of stereo imagery [4].Among these sensors ALOS/PRISM, Cartosat-1 and WorldView-2 can provide along-track stereoscopy data, resulting in a less complex image matching, as the two stereo images are captured with the same weather and illumination conditions.The DSMs from ALOS/PRISM and Cartosat-1 are sampled to a resolution of 5 m, while the DSM generated from WorldView-2 is characterized by a resolution of 1 m.Limited by the lower image resolution and across-track stereoscopy (acquired at different dates), the lower density of original matched pixels for RapidEye data only enables a 10 m resolution DSM.

• Satellite SAR data
TanDEM-X is a German satellite mission and consists of the two X-band SAR satellites TerraSAR-X and TanDEM-X, which are nearly identical in construction and fly in a close formation (see also [26]).The TanDEM-X mission provided for the first time single pass interferometric SAR data from space.The main objective of the mission is the acquisition of data for a global DSM with a spatial resolution of 12 m and an absolute height accuracy of 2 m and 4 m, respectively.The nominal revisit time of TanDEM-X is eleven days.However, by changing the incidence angle revisit times up to three days are possible.TanDEM-X provides a wide range of different acquisition modes, while usually it operates in strip map mode with VV (vertical transmit and vertical receive) polarization, for experimental issues a dual-pol mode is provided and in a special experimental mode even quad-pol acquisitions are feasible.Additionally, the TanDEM-X mission provides a ScanSAR, a spotlight, and a sliding spotlight imaging mode which can be operated at different combination of polarizations [26].The TanDEM-X2013 DSM (Figure 2i), was calculated with an experimental software package developed by the DLR-Microwaves and Radar Institute (HR) and has an original resolution of 6 m.Calculation steps for the processing of interferometric data are well described in [35], including a description of the parameters affecting the expected accuracy.The TanDEM-X specific interferometric configuration and its performance issues are extensively described in [26].As radar penetrates into the vegetation layer (for instance in forest areas), the DSM obtained represents a height somewhere between the top of the vegetation layer and the underlying ground surface [36].

•
Aerial image data One of the airborne stereo datasets was captured by an UltraCam Eagle camera on 16 June 2012 as part of the regularly scheduled aerial survey of Bavaria [37].The original color images have a ground sampling distance (GSD) of around 0.20 m, while the aerial data from 2003 was captured by the High Resolution Stereo Camera (HRSC) with GSD between 0.15 and 0.25 m [38].The HRSC camera was originally designed and developed for the Mars Missions [39], but one camera version has been used for land cover monitoring applications [38].The DSM used in this study (displayed in Figure 2c) was provided by the Institute of Space Sensor Technology and Planetary Exploration, German Aerospace Center (DLR).The data set from 2012 was prepared using the SGM method as in the case of the satellite images, but was generated with the remote sensing software package Graz (RSG) [37].The DSMs delivered from RSG included gaps, i.e., unmatched pixels, which were filled with the same method used for the satellite images [34].We refer to the DSM generated from these data as Aerial2012 (Figure 2h).
All data sets were converted and co-registered to the geographic reference system of the Bavarian State Forest Administration, which is the German Gauss-Krüger coordinate system, zone 4. The LiDAR2010 DSM was chosen as reference for the co-registration.The co-registration procedure is described in [34] and leads to accuracies of half a pixel.As the shifts among these data sets are very small, only linear shifting is considered.We use a simple version of least square matching.The shift distances in three dimensions are estimated via iterative 3D shifts adjustments based on the • Aerial image data One of the airborne stereo datasets was captured by an UltraCam Eagle camera on 16 June 2012 as part of the regularly scheduled aerial survey of Bavaria [37].The original color images have a ground sampling distance (GSD) of around 0.20 m, while the aerial data from 2003 was captured by the High Resolution Stereo Camera (HRSC) with GSD between 0.15 and 0.25 m [38].The HRSC camera was originally designed and developed for the Mars Missions [39], but one camera version has been used for land cover monitoring applications [38].The DSM used in this study (displayed in Figure 2c) was provided by the Institute of Space Sensor Technology and Planetary Exploration, German Aerospace Center (DLR).The data set from 2012 was prepared using the SGM method as in the case of the satellite images, but was generated with the remote sensing software package Graz (RSG) [37].The DSMs delivered from RSG included gaps, i.e., unmatched pixels, which were filled with the same method used for the satellite images [34].We refer to the DSM generated from these data as Aerial2012 (Figure 2h).
All data sets were converted and co-registered to the geographic reference system of the Bavarian State Forest Administration, which is the German Gauss-Krüger coordinate system, zone 4. The LiDAR2010 DSM was chosen as reference for the co-registration.The co-registration procedure is described in [34] and leads to accuracies of half a pixel.As the shifts among these data sets are very small, only linear shifting is considered.We use a simple version of least square matching.
The shift distances in three dimensions are estimated via iterative 3D shifts adjustments based on the minimization of the summarized squared height difference.Instead of using the whole image, we have manually selected some ground regions for the calculation.The calculated shifts in each of the three directions are documented in Table 2. To compare the datasets at pixel-level, all DSMs displayed in Figure 2 were resampled to a resolution of 2.5 m.
As mentioned above, some parts of the Traunstein forest were affected by a windthrow in 2007.The forest regions that were damaged are clearly visible by comparing the two LiDAR DSMs.The main affected area is marked by an ellipse in Figure 2a,b.This change in the forest structure can also be seen in all other DSMs captured after 2007.

Terrestrial Forest Inventory
Two terrestrial forest inventories were conducted at the study site by staff of the Chair of Forest Growth and Yield Science at TUM in summer 2008 and 2013, respectively, based on the standard design of the Bavaria State Forest Enterprise [40,41].
For this inventory, 221 field plots were placed in the forest area using a regular grid of 100 m × 100 m, with each one assumed to represent a forest area of 1 ha.The plot centers are permanently marked and each field plot consists of three concentric circles: an inner one with radius r 1 = 3.15 m (31.15 m 2 ), a middle one with radius r 2 = 6.31 m (125 m 2 ) and an outer one with radius r 3 = 12.62 m (500 m 2 ).In this study, the last (500 m 2 ) were used as field reference.The trees measured within the concentric circles depended on their diameter at breast height (DBH): all trees with a DBH < 10 cm, 10 cm ≤ DBH < 30 cm and DBH > 30 cm were callipered inside the smallest, middle, and largest circle, respectively.During field work, the species and age class were determined for each measured tree, which was then assigned to one of the following vertical layers: regeneration, understory, overstory or emergent trees (i.e., trees with crowns that emerge above the rest of the canopy).Then, at least one tree height was measured for every existing combination of tree species, age class and vertical layer affiliation occurring at a plot.If available, the height of at least two of the highest trees was measured for each of the following species: Picea abies, Pinus sylvestris, Abies alba, Larix decidua, Pseudotsuga menziesii, Fagus sylvatica, and Quercus spp.
In order to take only the inventory plots of the unchanged regions (after damage), a forest mask is manually extracted from the multispectral channels of WorldView-2 data.Thus, 184 plots are selected for the evaluation procedure.The distribution of the inventory plots is shown in Figure 3. Using these data, the top height H100 (which is defined as the average height of the 100 stems with the largest diameter at breast heights (DBHs) per hectare) was modeled for each plot using the algorithm described in [37].

Terrestrial Forest Inventory
Two terrestrial forest inventories were conducted at the study site by staff of the Chair of Forest Growth and Yield Science at TUM in summer 2008 and 2013, respectively, based on the standard design of the Bavaria State Forest Enterprise [40,41].For this inventory, 221 field plots were placed in the forest area using a regular grid of 100 m × 100 m, with each one assumed to represent a forest area of 1 ha.The plot centers are permanently marked and each field plot consists of three concentric circles: an inner one with radius r1 = 3.15 m (31.15 m 2 ), a middle one with radius r2 = 6.31 m (125 m 2 ) and an outer one with radius r3 = 12.62 m (500 m 2 ).In this study, the last (500 m 2 ) were used as field reference.The trees measured within the concentric circles depended on their diameter at breast height (DBH): all trees with a DBH < 10 cm, 10 cm ≤ DBH < 30 cm and DBH > 30 cm were callipered inside the smallest, middle, and largest circle, respectively.During field work, the species and age class were determined for

Methodology for Change Detection
The task of forest change detection using the introduced DSMs is applied to two scenarios.The first one is identifying decreases in forest height, typically caused by natural disturbance or forest management, and seldom due to illegal logging activities.The second one is estimating the increase in forest height due to natural tree growth.

Forest Decrease Detection
The forest decrease change detection techniques applied in this study rely on the computation of height differences.In the selected test region, trees are harvested applying a selective logging system, the windthrow in 2007 is therefore the main reason for height decrease in large regions.The DSMs taken before and after 2007 can be used as pre-and post-event data, respectively.Tree change monitoring from pixel based DSM data is a difficult task, as the tree crown representation in the DSMs may differ in different acquisitions due to local variations in reflection properties.Moreover, DSMs generated using stereoscopic or similar imaging data may contain imprecise height values due to matching errors.Finally, the resolution differences among the DSMs as well as tree growths might influence local change detection results.Resampling all DSMs to the same resolution and median filtering by a small kernel does not significantly change the original quality of the DSMs.Thus, in order to allow a reliable comparison, a median filter with window size of 5 × 5 pixels was applied to all DSMs before the post-change DSMs were subtracted from the pre-change DSMs.The height change ∆H(i.j) of pixel H(i, j) at coordinates i, j is defined as, where H pre (i, j) and H a f ter (i, j) represent the height values in the pre-and post-event DSMs after applying a median filtering, respectively.The windowsize of 5 × 5 pixels is set due to the resolution differences among the DSMs.Based on the generated height difference map, a change mask ChangeMask can be derived by setting a threshold T as: Remote Sens. 2017, 9, 287 9 of 26

Tree Growth Monitoring
Tree height and tree growth are both of great importance in forest monitoring, as they can indicate and predict forest biomass productivity [42].In this paper forest growth is assessed using periodic annual increment (PAI) values [37].
According to our investigation, the maximum height from the CHMs is the best predictor for H100 (Experiment I).The PAI is computed using the difference between the maximum heights of two CHMs within each inventory plot region (500 m 2 ), divided by the number of years between the data acquisition of both datasets: The LiDAR2001CHM was used as the initial height for each inventory plot, and represented as InvenHeight Pre in Equation (3).
Due to the storm event, many inventory plots have negative changes instead of positive tree growth change.Therefore, in a preprocessing step, all inventory plots with negative height changes were removed, which reduced the number of inventory plots for tree growth evaluation to 137.The PAIs calculated using LiDAR2010 as InvenHeight a f ter are recorded as reference.
We assume that tree heights have a positive correlation with tree ages, which influences the PAIs [31].Thus the obtained tree growth can be further analyzed by calculating the relationship between tree height and PAIs.

Experiment I-DSM Quality Evaluation
Forest change detection based on 3D data may be able to deliver good results, but it highly depends on the quality of the 3D data.For a meaningful interpretation, having deep understanding of the quality of the height information is necessary.For a continuous change monitoring it is almost impossible to have sufficient reference data for each dataset.In addition to that the storm event in 2007 has significantly changed the forest structure in Traunstein.Therefore, we have separated the DSMs into two groups for evaluation: (a) DSMs derived from the data acquired before the storm event in 2007; and (b) DSMs derived from the data acquired after the storm event in 2007.
The quality analysis was carried out separately using LiDAR data from 2001 for group (a) and LiDAR data from 2010 for group (b).Profiles of the DSMs were first analyzed visually; subsequently DSM heights were plotted against the LiDAR data and the top height H100 as derived for the field plots of the forest inventory.Finally, for a quantitative assessment, the DSMs were compared at pixel level for different land cover types using several statistical accuracy features.

Profile Comparison
To visually evaluate the quality of the DSMs, pixel-wise profiles were generated for two groups (Figure 4).The first group relates to the moderately-high resolution, space-borne data sets, which include ALOS/PRISM, Cartosat-1, RapidEye and TanDEM-X.The LiDAR data from 2010 is plotted as the height reference.The second group comprises the two very-high resolution airborne DSMs and the WorldView-2 DSM.For comparison, the two LiDAR DSMs are plotted in the same figure.
In Figure 4, the upper white line passes through three types of land cover, namely high-dense forest region, low-dense forest region and ground area (agricultural land or meadows).Except for RapidEye2009, all DSMs can clearly identify the forest boundary (between low-dense forest and ground), although the boundaries from ALOS2008 and TanDEM-X2013 are less sharp compared to the other investigated DSMs.
In the ground and high-dense forest areas, the difference between the DSMs derived from the high and low density LiDAR point clouds (LiDAR2010 and LiDAR2001, respectively) is not obvious.Using the LiDAR2010 as reference, the first profile plot displays the quality difference among the moderately high resolution DSMs.The profiles show that Cartosat2008 and TanDEM-X2013 match the reference data better than ALOS2008 and RapidEye 2009, because the canopy heights are smooth.Nevertheless, none of the mentioned DSMs is able to capture small gaps in the high-dense forest region, as the reference LiDAR DSM does.On the other hand, Cartosat2008 displays more details than TanDEM-X2013 showing some of the gaps in the low-dense forest region correctly, while RapidEye2009 can only show a rough height change along the profile line.
For both ground and high-dense forest areas, there is no clear difference between the LiDAR DSM and the very high resolution photogrammetrically generated DSMs.However, in the low-dense forest region, the DSMs vary in detail.As shown in the second profile plot in Figure 4, in the low-dense forest area the LiDAR DSMs, even the lower quality LiDAR2001 DSM, show many gaps in the canopy, which can only be partly seen in the HRSC2003, Aerial2012 and WorldView2012 DSMs.However, the upper layers of the canopy from HRSC2003, Aerial2012 and WorldView2012 are well matched with both LiDAR DSMs.A mismatch among the LiDAR2001, HRSC2003 and other DSMs can be identified in the second low-dense forest region (profile pixels 170-230).This can be explained by the storm event in 2007, which caused this region to change from high-dense to low-dense forest.

Grid-Based Comparison
In order to get a pixel-wise quality assessment in this section a grid-based evaluation was performed using the LiDAR 2001 and LiDAR2010 as reference data [20].The HRSC2003 has been compared against LiDAR2001, while all other DSMs have been compared against LiDAR2010.Differing from the workflow in [20], we use two different grid cell sizes: one of 500 m 2 , which corresponds to the size of the largest concentric circle of the forest inventory plots in Bavaria Germany (see Section 2.3), and another of 1 ha (10,000 m 2 ), which is an important unit for forest management.Focusing on quality evaluation for the forest regions, we use only the grids that are full or partly covered by the forest mask (Figure 3) for the comparison.Figure 5 shows scatterplots for the 500 m 2 grid cells; for each cell, the mean height was derived from the DSM and compared to the corresponding mean height of the LiDAR2001 and LiDAR2010 data.In total, 6931 mean height

Grid-Based Comparison
In order to get a pixel-wise quality assessment in this section a grid-based evaluation was performed using the LiDAR 2001 and LiDAR2010 as reference data [20].The HRSC2003 has been compared against LiDAR2001, while all other DSMs have been compared against LiDAR2010.Differing from the workflow in [20], we use two different grid cell sizes: one of 500 m 2 , which corresponds to the size of the largest concentric circle of the forest inventory plots in Bavaria Germany (see Section 2.3), and another of 1 ha (10,000 m 2 ), which is an important unit for forest management.Focusing on quality evaluation for the forest regions, we use only the grids that are full or partly covered by the forest mask (Figure 3) for the comparison.Figure 5 shows scatterplots for the 500 m 2 grid cells; for each cell, the mean height was derived from the DSM and compared to the corresponding mean height of the LiDAR2001 and LiDAR2010 data.In total, 6931 mean height values were compared.The names of the DSMs are reported below the y-axis.
Remote Sens. 2017, 9, 287 11 of 26 than in the profile comparison, and it yields similar results as the ones obtained by Cartosat2008 and TanDEM-X2013.Pearson's correlation coefficients were also derived for the mean, maximum and 68.3% and 95% quantiles of each DSM height metrics as derived from the grid cells [20,43,44].Table 3 shows the correlation coefficients for the 500 m 2 and the 1 ha grid cells, in which mean height values yield the highest correlation coefficients for all DSMs.Aerial2012 achieved the highest correlation coefficients, while WorldView2012 showed satisfactory results, and RapidEye2009 had the lowest correlation coefficients.For the larger (1 ha) grid size the correlation coefficients are overall slightly higher than for the 500 m 2 grid.The correction coefficients values are generally higher than those in the next section, as only grids intersecting the forest mask have been taken into account, as described above.In this comparison, the Aerial 2012 data match the LiDAR data best, while WorldView2012 and HRSC2003 achieve similar results.In this grid-based comparison, ALOS2008 performs much better than in the profile comparison, and it yields similar results as the ones obtained by Cartosat2008 and TanDEM-X2013.
Pearson's correlation coefficients were also derived for the mean, maximum and 68.3% and 95% quantiles of each DSM height metrics as derived from the grid cells [20,43,44].Table 3 shows the correlation coefficients for the 500 m 2 and the 1 ha grid cells, in which mean height values yield the highest correlation coefficients for all DSMs.Aerial2012 achieved the highest correlation coefficients, while WorldView2012 showed satisfactory results, and RapidEye2009 had the lowest correlation coefficients.For the larger (1 ha) grid size the correlation coefficients are overall slightly higher than for the 500 m 2 grid.The correction coefficients values are generally higher than those in the next section, as only grids intersecting the forest mask have been taken into account, as described above.In this section, the normalized canopy heights extracted from the DSMs are evaluated with the aid of Forest inventory plots introduced in Section 2.3.The top height H100 is used as reference height for each inventory plot.
Nine CHMs were generated by subtracting the terrain height from the DSMs.For this purpose, a digital terrain model (DTM) was generated from the last returns of the 2010 airborne laser scanning data.The H100 height from 2008 and 2013 are used as reference heights for each inventory plot.The locations of these inventory plots are indicated by black dots in Figure 3. Thus, for each inventory plot, the corresponding pixels in the CHMs represented by the inventory plots are used for the evaluation.The DSMs acquired before 2010 were compared against the H100 heights from 2008.The DSMs acquired after 2010 were compared against the H100 heights from 2013.
The maximum heights from the CHMs yielded the highest correlation coefficients with the H100 heights in [37].Thus, following the same workflow, we plot the maximum CHM heights separately for all datasets against the H100 height of the field plots in Figure 6.The CHMs from LiDAR2010, WorldView2012, and Aerial2012 match much better the H100 heights than the other datasets.The main scatter points from LiDAR2001 and HRSC2003 cluster lower values than the center 1:1 line, which means that, due to the time difference the H100 heights are significantly higher than the height from these two DSMs.ALOS2008, Cartosat2008, RapidEye2009 and TanDEM-X2013 show many inventory plots with lower height estimation and therefore less correlation to the H100 heights.
Besides the maximum CHM values, the mean, 68.3% and 95% quantiles within the forest inventory plot regions are derived.The Pearson's correlation coefficients of these DSM metrics with the H100 heights are listed in Table 4.By taking the H100 heights (h) as reference, the relative root mean square error (RRMSE) is calculated for each CHM as: where n is the total number of inventory plots, ĥ represents the canopy height calculated from the CHMs by taking the maximum, mean, 68.3% or 95% quantiles of the pixels covered by the inventory plot, and h is the H100 height of the corresponding plot.In   By using the maximum heights of each inventory plot (as shown graphically in Figure 6), Aerial2012 has the highest correlation with the inventory H100 heights and the lowest RMSE,  By using the maximum heights of each inventory plot (as shown graphically in Figure 6), Aerial2012 has the highest correlation with the inventory H100 heights and the lowest RMSE, followed by the LiDAR2010 and WorldView2012 DSMs, while Cartosat2008 yields slightly worse results than WorldView2012.Correlation coefficients based on the 95% quantile are similar to the results derived from the maximum height value.

Statistical Evaluation for Different Land Cover Types
In addition to the visual comparison and correlation coefficients, the DSMs have been evaluated using more statistical features.Over the 13 years in which the data were acquired, some changes may have influenced the forest structure.Therefore, we have manually selected some high-dense and low-dense forest areas unrelated to visible change.
As the inventory plots only cover the forest region, LiDAR2001 and LiDAR2010 are again used as reference for the two defined groups.As shown in Table 5, the accuracy assessment is based on the absolute vertical height differences (∆h) between the DSMs and the reference height model.The mean error ( μ) represents the mean absolute shift value for all pixels, while the standard deviation ( σ) is derived from the difference of the two DSMs, which also indicates the surface smoothness of the generated DSM.Another robust, statistically-based accuracy measurement is the Normalized Median Absolute Deviation (NMAD) [44].NMAD can be used to evaluate the quality of the DSMs in the presence of outliers and non-normal distributions, and it is proportional to the median of the absolute differences between errors and the median error.It can be considered as an estimate for the standard deviation which is more resilient to outliers [44].In Table 6 the calculated "ME", "STD" and "NMAD" values illustrate the characteristics of all DSMs.More than 10,000 pixels are selected for each land cover type.The best results in each comparison group are highlighted in bold with black color and the worst results are highlighted in grey-bold.However, by comparing the accuracies of these DSMs for each land cover type, it becomes clear that DSMs in high-density forest regions generally match better the reference DSM than in the low-density ones.In high-density forest, the mean errors range from 3.62 to 5.72 m: the range is very small, considering the resolution of the original images.The differences are also larger for low-density forest regions.All DSMs are instead less accurate in low-density forest areas.One reason could be that in low-density forests, LiDAR data tend to measure terrain height rather than canopy height.Matching the optical image is also more challenging in low-density forest regions, as explained in [4], which influences the quality of the photogrammetrically derived DSMs.

Reference Map Generation
In this experiment, we test the capability of height difference maps as derived from the DSMs to detect windthrow areas, caused by a storm event in 2007 in the described test site.Due to the lack of field data, we use the height difference map derived from both available LiDAR data sets as reference (Figure 7a).
Remote Sens. 2017, 9, 287 15 of 26 In Table 6 the calculated "ME", "STD" and "NMAD" values illustrate the characteristics of all DSMs.More than 10,000 pixels are selected for each land cover type.The best results in each comparison group are highlighted in bold with black color and the worst results are highlighted in grey-bold.However, by comparing the accuracies of these DSMs for each land cover type, it becomes clear that DSMs in high-density forest regions generally match better the reference DSM than in the low-density ones.In high-density forest, the mean errors range from 3.62 to 5.72 m: the range is very small, considering the resolution of the original images.The differences are also larger for low-density forest regions.All DSMs are instead less accurate in low-density forest areas.One reason could be that in low-density forests, LiDAR data tend to measure terrain height rather than canopy height.Matching the optical image is also more challenging in low-density forest regions, as explained in [4], which influences the quality of the photogrammetrically derived DSMs.

Reference Map Generation
In this experiment, we test the capability of height difference maps as derived from the DSMs to detect windthrow areas, caused by a storm event in 2007 in the described test site.Due to the lack of field data, we use the height difference map derived from both available LiDAR data sets as reference (Figure 7a).As discussed in the previous section, although the two LiDAR data sets exhibit similar characteristics, the elevation values in the forest regions do not perfectly match, even in the regions with no changes.Therefore, the same median filter as mentioned in Equation ( 1) is applied separately for the two LiDAR DSMs.The filtered height difference map is shown in Figure 7b.Here, As discussed in the previous section, although the two LiDAR data sets exhibit similar characteristics, the elevation values in the forest regions do not perfectly match, even in the regions with no changes.Therefore, the same median filter as mentioned in Equation ( 1) is applied separately for the two LiDAR DSMs.The filtered height difference map is shown in Figure 7b.Here, the windthrow areas are well highlighted in red, and the noise in the central area is largely removed.The light green color in Figure 7 indicates a height increase from natural tree growth, in which Figure 7a assumes a value of about 5-10 m for the nine year time difference (around 0.6 m per year), which is slightly higher than the 0.3-0.4m annual growth speed obtained in [37].However, in the filtered height change map (Figure 7b), the height increase is only 3-6 m.Some regions exhibit a height change of about 9 m, which is assumed to be related to young, fast-growing trees.
A major task after storms in forest areas is to determine and delineate the damaged regions.Thus, we generated a reference mask for the windthrow areas based on the refined height difference map (Figure 7b).First, a height threshold of 12 m (value selected according to the recommendation given by Knoke [45]) was applied to obtain an initial change mask, which shows the areas where the filtered difference map indicates a reduction in height of 12 m or more.Next, morphological erosion and dilation operators were applied to clean the generated change reference mask.In the morphological filter procedure, a squared structuring element of size of 3 × 3 was used, and regions smaller than 100 pixels were removed.

DSM Based Storm Damage Detection
To generate the height change maps we used LiDAR2001 and HRSC2003 as pre-change (i.e., before storm event) DSMs, and all DSMs acquired after 2007 as post-change (after storm) DSMs. Figure 8 shows results from using LiDAR2001 as the pre-change DSM, while results based on HRSC2003 are displayed in Figure 9.To allow comparison, both figures use the same color bar as in Figure 7a,b: red represents changes in the forest structure related to the storm event, while green shows positive changes due to tree growth.
Focusing first on Figure 8, the DSMs show promising results when compared to the change reference mask (Figure 7c), with each map highlighting the forest damages with varying degrees of accuracy.Benefitting from the higher quality of WorldView2012 and Aerial2012 data, the difference maps in Figure 8d,e match better the reference data.Although in Figure 8b,f the residential regions exhibit more noisy structures, the changes highlighted in the forest region show similar results as in Figure 8d,e the windthrow areas are well highlighted in red, and the noise in the central area is largely removed.The light green color in Figure 7 indicates a height increase from natural tree growth, in which Figure 7a assumes a value of about 5-10 m for the nine year time difference (around 0.6 m per year), which is slightly higher than the 0.3-0.4m annual growth speed obtained in [37].However, in the filtered height change map (Figure 7b), the height increase is only 3-6 m.Some regions exhibit a height change of about 9 m, which is assumed to be related to young, fast-growing trees.
A major task after storms in forest areas is to determine and delineate the damaged regions.Thus, we generated a reference mask for the windthrow areas based on the refined height difference map (Figure 7b).First, a height threshold of 12 m (value selected according to the recommendation given by Knoke [45]) was applied to obtain an initial change mask, which shows the areas where the filtered difference map indicates a reduction in height of 12 m or more.Next, morphological erosion and dilation operators were applied to clean the generated change reference mask.In the morphological filter procedure, a squared structuring element of size of 3 × 3 was used, and regions smaller than 100 pixels were removed.

DSM Based Storm Damage Detection
To generate the height change maps we used LiDAR2001 and HRSC2003 as pre-change (i.e., before storm event) DSMs, and all DSMs acquired after 2007 as post-change (after storm) DSMs. Figure 8 shows results from using LiDAR2001 as the pre-change DSM, while results based on HRSC2003 are displayed in Figure 9.To allow comparison, both figures use the same color bar as in Figure 7a,b: red represents changes in the forest structure related to the storm event, while green shows positive changes due to tree growth.
Focusing first on Figure 8, the DSMs show promising results when compared to the change reference mask (Figure 7c), with each map highlighting the forest damages with varying degrees of accuracy.Benefitting from the higher quality of WorldView2012 and Aerial2012 data, the difference maps in Figure 8d,e match better the reference data.Although in Figure 8b,f the residential regions exhibit more noisy structures, the changes highlighted in the forest region show similar results as in Figure 8d,e.As the HRSC2003 data were captured two years later than the LiDAR2001 DSM, in Figure 9 the green color, representing positive height differences, is not as dark as in Figure 8.In addition to tree growth, the lighter green may also reflect the lower quality of the HRSC DSM, and the fact that HRSC is smoother than LiDAR2001.In Figure 9, most of the regions which suffered damages from the storm (in red) match the ones in Figure 8, but the size of the regions is different in some cases.
To further analyze the height difference maps, the forest damage changes are evaluated.For forest damage, change masks are generated by taking the pixels with a height decrease of more than 12 m.The resulting change masks from Figure 8 are displayed in Figure 10a-f.The results from Figure 9 are shown in Figure 10g-l.In each subfigure, the forest mask is shown in light grey and used as background.The extracted change mask is overlaid with the reference change mask.True positives are represented in gray, while red and blue represent false positives and false negatives, respectively.No further refinements were carried out to keep unaltered the original detail of the datasets.
A visual comparison shows that all change masks in Figure 10 correctly highlight the regions destroyed by the storm as in the reference change mask (Figure 9c), but with different accuracies.By using the LiDAR2001 as pre-change DSM, the results from Cartosat2008 (Figure 9b), WordView2012 (Figure 9d) and Aerial2012 (Figure 9e) are similar and exhibit less false positives than others.Despite of a few false alarms, ALOS2008 (Figure 9a) and even RapidEye2009 (Figure 9c) correctly highlight the changed regions.Using HRSC2003 as pre-change DSM leads instead to more false positives and false negatives than using LiDAR2001.The changed regions are correctly located but with incorrect boundaries and size.
To quantitatively evaluate the accuracy of these change masks, the user and producer accuracy, together with the overall and Kappa statistic accuracy [46] are reported in Table 7 for the 12 change masks.When the LiDAR2001 data are adopted, all change masks achieve Kappa accuracies higher than 0.5.Based on the HRSC2003, two of the Kappa accuracies are below 0.5.Within both comparisons, WorldView2012 achieves the best producer, overall and Kappa accuracies.Aerial2012 exhibits the best user accuracy across both comparisons.As the HRSC2003 data were captured two years later than the LiDAR2001 DSM, in Figure 9 the green color, representing positive height differences, is not as dark as in Figure 8.In addition to tree growth, the lighter green may also reflect the lower quality of the HRSC DSM, and the fact that HRSC is smoother than LiDAR2001.In Figure 9, most of the regions which suffered damages from the storm (in red) match the ones in Figure 8, but the size of the regions is different in some cases.
To further analyze the height difference maps, the forest damage changes are evaluated.For forest damage, change masks are generated by taking the pixels with a height decrease of more than 12 m.The resulting change masks from Figure 8 are displayed in Figure 10a-f.The results from Figure 9 are shown in Figure 10g-l.In each subfigure, the forest mask is shown in light grey and used as background.The extracted change mask is overlaid with the reference change mask.True positives are represented in gray, while red and blue represent false positives and false negatives, respectively.No further refinements were carried out to keep unaltered the original detail of the datasets.
A visual comparison shows that all change masks in Figure 10 correctly highlight the regions destroyed by the storm as in the reference change mask (Figure 9c), but with different accuracies.By using the LiDAR2001 as pre-change DSM, the results from Cartosat2008 (Figure 9b), WordView2012 (Figure 9d) and Aerial2012 (Figure 9e) are similar and exhibit less false positives than others.Despite of a few false alarms, ALOS2008 (Figure 9a) and even RapidEye2009 (Figure 9c) correctly highlight the changed regions.Using HRSC2003 as pre-change DSM leads instead to more false positives and false negatives than using LiDAR2001.The changed regions are correctly located but with incorrect boundaries and size.
To quantitatively evaluate the accuracy of these change masks, the user and producer accuracy, together with the overall and Kappa statistic accuracy [46] are reported in Table 7 for the 12 change masks.When the LiDAR2001 data are adopted, all change masks achieve Kappa accuracies higher than 0.5.Based on the HRSC2003, two of the Kappa accuracies are below 0.5.Within both comparisons, WorldView2012 achieves the best producer, overall and Kappa accuracies.Aerial2012 exhibits the best user accuracy across both comparisons.

DSM Based Tree Growth Monitering
We calculate the PAIs with the described method in each inventory plot, taking the LiDAR2001 as reference pre-change data.The PAIs sets calculated from satellite data and aerial imagery are compared to the PAIs calculated by using two LiDAR datasets.Figure 11 shows the results.and RapidEye2009 are obviously larger than those from LiDAR2010.
The five best performing DSMs visually observed from Figure 11 are further used to analyze the relationship between tree height and PAIs (Figure 12).Most plots with lower (younger) trees have larger PAIs with respect to plots with higher trees.In the test region, most trees are 15-30 m high, with PAIs of 0.3-0.5 m per year.A linear least square fitting regression is applied to get the relationship between the tree heights and the PAIs.The obtained regression lines are plotted in Figure 12.All datasets follow similar trends and exhibit decreasing PAIs as the height increase.The scatter plots in Figure 11 show that most of the PAIs are in the range between 0 and 1 m.The results from Aerial2012 and WorldView2012 match well the results derived from LiDAR2010.PAIs detected from Cartosat2008 and TanDEM-X2013 show some differences in comparison to LiDAR2010, but most points are still quite close to the diagonal.The PAIs detected from ALOS2008 and RapidEye2009 are obviously larger than those from LiDAR2010.

Discussion
In this paper, DSMs from nine different sensors were analyzed and adopted for change monitoring in a complex forest region.Height, as a fundamental feature for forest monitoring, is related to management activities and natural disasters.Sensors have a defined life time, and limited acquisition capacity.A continuous monitoring with the same sensor is in many cases impossible to The five best performing DSMs visually observed from Figure 11 are further used to analyze the relationship between tree height and PAIs (Figure 12).Most plots with lower (younger) trees have larger PAIs with respect to plots with higher trees.In the test region, most trees are 15-30 m high, with PAIs of 0.3-0.5 m per year.A linear least square fitting regression is applied to get the relationship between the tree heights and the PAIs.The obtained regression lines are plotted in Figure 12.All datasets follow similar trends and exhibit decreasing PAIs as the height increase.

Discussion
In this paper, DSMs from nine different sensors were analyzed and adopted for change monitoring in a complex forest region.Height, as a fundamental feature for forest monitoring, is related to management activities and natural disasters.Sensors have a defined life time, and limited acquisition capacity.A continuous monitoring with the same sensor is in many cases impossible to obtain.Thus an extensive understanding of the characteristics, especially the quality of height information from various sensors is necessary.Earlier studies have often assessed one remote sensing data set [6,9,14,19] or compared two to four types of data sets [3,4,16,19,20,23,27,37,47].In this paper, besides forest change detection, in the experimental part DSMs from both satellite and aerial images were compared to reference LiDAR DSMs and to forest inventory plots.Moreover, within the analysis of satellite photogrammetry results, the applicability of stereo imagery with resolutions ranging from 0.5 m to 6.5 m has been assessed.
The evaluated datasets available for the "Stadtwald Traunstein" represent a variety of sensors with different spatial resolution and stereo capability.These sensors have diverse advantages and disadvantages.LiDAR is usually considered as the most suitable remote sensing technique to acquire forest height data.However, compared to photogrammetric surveys, LiDAR data are acquired with lower flying heights and slower flying speeds which results in higher costs [48].Thus, it can be challenging to monitor large forest regions within a short time period with LiDAR.
Nowadays, DSMs from airborne stereo sensors benefit from the new matching techniques and can achieve accuracies similar to LiDAR [37].In previous literature [20], stereo images were matched with different window sizes, which influenced the DSM quality.However, currently available dense matching approaches have greatly improved accuracy and are not anymore influenced by this parameter.Therefore, it may be possible to adopt space-borne sensors with relatively low resolution to provide sufficient accurate DSMs [49].Space-borne sensors exhibit larger coverage in comparison to airborne sensors, which is very helpful for monitoring larger regions.When considering the cost, lower resolution images like Cartosat-1 are considerably cheaper and more easily available for monitoring larger regions than very high resolution WordView-2 data.Among active sensors, TanDEM-X is a new active satellite constellation operation in the X-band frequency, which provides a worldwide consistent DSM.The X-band has less penetration compared to other longer wavelengths such as L-or P-band, and allows estimating the forest canopy surface.TanDEM-X can continually operate from a space platform and is not influenced by weather.In general, it is difficult to quantify the applicability of the mentioned sensors.A detailed understanding of the quality of the specific DSMs and their behavior in forest change detection can help to select some of them for practical application.

Accuracy Evaluation of the DSMs
Forest change detection based on DSMs delivers promising results, but highly depends on the quality of the DSM, making a quality evaluation of the DSMs is necessary.Using field data, a inventory plot area of 500 m 2 is statistically extrapolated to represent areas of 1 ha.In order to get pixel-wise evaluation results, LiDAR data are assumed to provide the highest level of detail and are therefore used as additional reference.The Traunstein forest is one of the best investigated long term forest test sites in Germany.The LiDAR2010 dataset has been used and evaluated already by other publications for height and timber volume estimations [4,15,18].Up to now, LiDAR is considered as the most accurate remote sensing technology for forest height determination.In addition, the LiDAR DSM covers the complete area of the investigated forest canopy and is better suited with respect to field inventory plots for comparing DSMs derived from other sensors.
DSMs derived from the HRSC data (HRSC2003) as well as from aerial photographs (Aerial2012) show fewer details than LiDAR and display a weakness in the gap assessment, most probably due to shade effects.The DSM from very high resolution satellite data (WorldView-2) can almost reach similar accuracies as the aerial data, and shows sharp boundaries of forest regions.Although the DSMs from high resolution optical ALOS/PRISM, Cartosat-1 and RapidEye (~2.5-5 m pixel size), and InSAR (TanDEM-X) satellite data are less accurate than WorldView2012, they can still distinguish high from low level objects (e.g., trees from ground), and show quite accurately the average canopy height of densely forested regions.
The quality of the generated DSMs is related to the resolution of the sensors, even though the DSM quality in high density forest region is less sensitive to image resolution than in low-density ones.As already shown, Aerial2012 has the highest accuracy, followed by HRSC2003 and WorldView2012.Based on the latest matching and post-processing techniques, the quality of the DSMs derived from stereo data has greatly improved.Previous studies have compared DSMs to investigate correlations at forest inventory plot level [20,43].Because they analyzed different test sites with different forest types, results from these studies cannot be directly compared, but they can be selectively used as references.The correlation coefficients between LiDAR and Aerial2012 calculated within this paper can reach 0.99, much higher than those reported in [20], in which the highest correlation between LiDAR and Photogrammetry-LiDAR(to calculate CHM) was 0.95.The lowest correlation we obtained (from RapidEye2009) is about 0.94.These results may reflect the latest dense matching DSM generation technique.Moreover, different forest types with different structure might also be the reason for other correlation results.In [43], the correlation at the plot level (size 500 m 2 ) was analyzed for four different land cover types.As our results do not differentiate tree species, we compare them with the species mixture class, which was 0.64 in [43].In our study, referring to H100 heights, the DSMs from Cartosat-1 reached a similar correlation value of 0.64.Aerial2012 reaches the highest correlation coefficients when the maximum height was calculated, which was around 0.83.The correlation coefficient between LiDAR2010 and 2008-H100 was 0.79.This result could be explained by the time difference between the DSMs and inventory H100, as the LiDAR2010 and 2008-H100 have two years distance and Aerial2012 and 2013-H100 have only one year time difference.
The DSM quality evaluation in forest regions is not as easy as in urban regions, because the measured objects "trees" do not have as distinct heights and boundaries as buildings.In addition, taking LiDAR as reference is controversial.Firstly, LiDAR data are also remote sensing data, and thus they still show some discrepancies to field measurements.Secondly, the CHMs derived from remote sensing data give a mixture of canopy objects and gaps.Therefore not all of the pixels in the CHM are indicating tree height.However, typical "heights" associated to management relevant development stages are 0-10 (−12) m (stand established), 10-12 m to 20-24 m (qualification, stabilizing measures) and above 20-24 m (dimensioning, up to harvesting, establishing of a new stand by successive extraction of single trees, etc.) [45].Such categories and additionally gaps and disturbances are interesting for the practitioner.When using a region-based height analyses by taking the average or maximum values of each region, the mentioned problems regarding distinct heights and boundaries can be ignored.Thirdly, the data used in this paper have been captured for various purposes over a 13-year timespan.Besides tree growth, the shape of the tree canopies can change in different years and seasons.Two LiDAR data from 2001 and 2010 were used as reference data.Thus, the longest distance between LiDAR measurement and another remote sensing data acquisition is three years.In this study, almost all photogrammetric data were captured in leaf-on season, which enables a better accuracy for the image matching.Unfortunately the LiDAR data were captured in early spring, which is leaf-off season.However, the main tree species in the test site are coniferous trees, which are less influenced by seasonal differences.More detailed research is required to analyze the influences of leaf-on/leaf-off status to the forest heights obtained from remote sensing sensors.Furthermore, the statistical accuracy can certainly be improved if the data are captured in the same year and month.However, a time difference does allow change monitoring.

Change Detection
The change detection accuracy is another important aspect of DSM accuracy.Forest decrease/growth analysis is still an on-going research topic for both environmental research and forest management.In this paper, height change maps and detected forest change masks were generated and evaluated by a comparison with reference change masks.Change detection in forest regions using optical 2D remote sensing data can be a challenging task, because of the similar spectral and textural features between trees and other vegetation.However, changes in height can provide an efficient and more accurate indication for forest change.Focusing on forest changes, our previous results [34,50] prove the good performance of the DSMs from satellite images.The quality of these DSMs directly influences the change detection accuracy [34].Therefore, comparing the change detection results from various sensor combinations helps to select proper data sources and enhance the forest change algorithms.
By using LiDAR2001 and HRSC2003 as the pre-change DSMs separately, two groups of results are generated.The height change accuracies in Table 7 are much higher when LiDAR2001 data are used as the pre-change DSM, possibly because the change reference is also produced using LiDAR2001 data as the pre-change DSM.When using LiDAR2001 as the pre-change DSM, all results show the correct location and approximate size of the damages.More false detections were introduced when using less accurate DSM, for example the ones derived from RapidEye2009.
Although the results from using HRSC2003 as the pre-change DSM show some false positives, the locations of the main changes are still highlighted correctly.The comparison confirms that satellite data are well suited to support a fast response for detecting damage from sudden events, such as the investigated storm event Kyrill, because all main changes are detected through the automatic change detection procedure.According to the accuracy requirements and available funding, proper sensor combinations can be selected.Forest changes due to management activities can also be detected through comparing these DSMs [29].However, these changes are not analyzed in this paper.In addition, more research is required to analyze the relationship between forest density and storm damage.
In addition to forest damage monitoring, the tree growth was extracted and evaluated.Our results confirm the findings from [37], whose analysis of the same test site was mainly based on aerial images.Forest growth is barely detectable when only spectral information is available.By dividing the height difference with the time distance (measured in years), PAIs were derived as an indicator for the average forest growth per year.In our study, more data sources, especially DSMs from satellites, were tested; these data were acquired over a relatively long time interval.Therefore, it is possible to provide reasonable average tree growth PAIs which are less affected by annual variations of growth.The relationship between tree height and PAIs also matches the results from [37], in which shorter trees, which are assumed to be younger, grow faster than taller trees.Besides tree age, tree species, environment [37,51] and health conditions might have influenced the tree growth.For mixed temperate forests an approximate mean annual growth rate of 0.30-0.40m can be expected over a 100-year period [37].With WorldView-2 data and aerial datasets, as shown in Figure 12, the main tree height in this forest is 25-30 m and the results for the PAIs are in the range of 0.30-0.50m.The change detection is performed based on the assumption that the pre-event LiDAR data are available, which is usually true for most parts of Europe.

Conclusions
Detecting and predicting stand productivity is fundamental to forest monitoring.Besides, LiDAR high spectral and spatial resolution stereo imaging systems from both space-borne and airborne sensors are capable of providing height information.Therefore, adopting remote sensing data in forest monitoring can be a good choice, especially for frequent forest monitoring purposes.As stated in [20], the DSM quality could be further improved with upgraded image acquisition sensors and stereo-matching techniques, which are now available through new camera systems and computer vision approaches [17,25,33,34,48].Therefore, a systematic investigation of DSMs derived from the different currently available sensors is performed in this research to show their potential usability in change detection in detail.
A novel approach of this paper was not only forest decrease change and forest growth monitoring, but also a quality evaluation of DSMs from nine different sensors.As an overall conclusion, the data can be categorized in three sets according to their accuracy.The first set includes the DSMs from LiDAR and stereo imagery with less than 1 m resolution derived from airborne and WorldView-2 data.These are the most suitable images for a detailed monitoring of forest structure and height properties.With these data, even medium size gaps (5-10 m diameter) can be analyzed and small changes in the mean stand to single tree height detected.These data might be a good replacement of LiDAR for regular forest surveying tasks.The second set is the 2.5 m resolution Cartosat-1 and TanDEM-X data, where the quality is degraded but still sufficient for estimating the main changes, also in smaller areas.TanDEM-X has its strength in its large scale coverage, its operation independent of cloud cover, and its capability of delivering worldwide monitoring results.Data in the third set are ALOS/PRISM and RapidEye.The resolution and stereo capability influence the quality of DSMs derived from Rapid Eye, which is significantly lower than DSMs derived from other sensors.This is also caused by the cross-track instead of along-track capturing of stereo pairs, meaning that the time difference between stereo pair acquisition is at least several days.The cross-track RapidEye data are not the first choice for DSM generation, but based on the latest DSM generation approach they are able to represent roughly the forest canopy height and locate the main change regions.The ALOS/PRISM sensor has stopped acquiring data since 21 April 2011, and its valuable historical data would be used as pre-change event DSM.With the available DSMs from the pre-and post-change event, forest change can be highlighted.Problems may arise if forests are mixed with other land covers, like buildings.In this case spectral information will be needed to highlight the tree covered areas as a preprocessing step.
We conclude that satellite data have a strong potential for future forest height and damage (change) monitoring, especially since more and more sensors providing DSMs are and will be in orbit within the next years, making the rapid response to forest calamities through automatic change detection using satellite data feasible.In combination with information from very high to high resolution multi-seasonal multispectral data, an annual update of forest data bases, capable to support regular management and create the basis for fast response needs seems possible [41].After large windthrow events, space-borne stereo data could provide, in a short time, a first order evaluation of the affected forest area and even allow a rough estimate of the damaged timber volume.

Figure 1 .
Figure 1.Geographical location of study site Stadtwald Traunstein.The left part of the figure shows the location of the study site in Bavaria, Germany.The right part of the figure shows a color-infrared orthophoto of the study site (coordinates in Gauss-Krüger system).

Figure 1 .
Figure 1.Geographical location of study site Stadtwald Traunstein.The left part of the figure shows the location of the study site in Bavaria, Germany.The right part of the figure shows a color-infrared orthophoto of the study site (coordinates in Gauss-Krüger system).

Figure 3 .
Figure 3. Distribution of 184 selected inventory plots overlaid on the forest mask which was manually digitized from WorldView-2 data, where black pixels represent the center of each plot and the forest mask is shown in light grey color.

Figure 3 .
Figure 3. Distribution of 184 selected inventory plots overlaid on the forest mask which was manually digitized from WorldView-2 data, where black pixels represent the center of each plot and the forest mask is shown in light grey color.
Remote Sens. 2017, 9, 287 10 of 26 the low-dense forest area the LiDAR DSMs, even the lower quality LiDAR2001 DSM, show many gaps in the canopy, which can only be partly seen in the HRSC2003, Aerial2012 and WorldView2012 DSMs.However, the upper layers of the canopy from HRSC2003, Aerial2012 and WorldView2012 are well matched with both LiDAR DSMs.A mismatch among the LiDAR2001, HRSC2003 and other DSMs can be identified in the second low-dense forest region (profile pixels 170-230).This can be explained by the storm event in 2007, which caused this region to change from high-dense to low-dense forest.

Figure 4 .
Figure 4. Profiles of the DSMs covering low-dense forest, high-dense forest and ground areas.(The DSMs are separated into two groups-above: moderately-high resolution DSMsl and below: very high resolution DSMs.)

Figure 4 .
Figure 4. Profiles of the DSMs covering low-dense forest, high-dense forest and ground areas.(The DSMs are separated into two groups-above: moderately-high resolution DSMsl and below: very high resolution DSMs.)

Figure 7 .
Figure 7. Reference change mask generation procedure: (a) Original height difference map; (b) refined height difference map (in (a,b), the red areas represent a decrease in height, while green ones represent the increase in height and white represents regions with no height difference and no change); and (c) generated change reference mask (green: a decrease in height of 12 m or more; grey: forest mask extracted from WorldView-2 data).

Figure 7 .
Figure 7. Reference change mask generation procedure: (a) Original height difference map; (b) refined height difference map (in (a,b), the red areas represent a decrease in height, while green ones represent the increase in height and white represents regions with no height difference and no change); and (c) generated change reference mask (green: a decrease in height of 12 m or more; grey: forest mask extracted from WorldView-2 data).

Table 1 .
Summary of the remote sensing datasets.

Table 3 .
Correlation coefficients between height metrics of image-and LiDAR-derived DSMs for 500 m 2 and 1 ha grid cells (highest and lowest values are marked in bold).
4.3.Comparison of DSMs to the Top Height H100 at Forest Inventory Plots

Table 4
, the best results with highest correlation and lowest RRMSE are marked in bold black, while the worst result with lowest correlation and highest RRMSE are marked in bold grey.Remote Sens. 2017, 9, 287 13 of 26 correlation and lowest RRMSE are marked in bold black, while the worst result with lowest correlation and highest RRMSE are marked in bold grey.

Table 4 .
Correlation coefficients between DSM metrics and with the H100 height of forest inventory plots and RRMSE by taking H100 from 2008 and from 2013 as reference (the best and worst results are marked in bold).

Table 4 .
Correlation coefficients between DSM metrics and with the H100 height of forest inventory plots and RRMSE by taking H100 from 2008 and from 2013 as reference (the best and worst results are marked in bold).

Table 5 .
Accuracy measures for DSM quality assessment.

Table 6 .
Quality assessment for low-density and high-density forest based on the quality measures in Table5by using LiDAR data from 2001 as reference for HRSC2003 and from 2010 as references for the other DSMs (the highest and lowest values are marked in bold). .

Table 7 .
Comparison of the change masks generated from different combination of DSMs (highest values are marked in bold).

Table 7 .
Comparison of the change masks generated from different combination of DSMs (highest values are marked in bold).