Towards Automated Characterization of Canopy Layering in Mixed Temperate Forests Using Airborne Laser Scanning

Canopy layers form essential structural components, affecting stand productivity and wildlife habitats. Airborne laser scanning (ALS) provides horizontal and vertical information on canopy structure simultaneously. Existing approaches to assess canopy layering often require prior information about stand characteristics or rely on pre-defined height thresholds. We developed a multi-scale method using ALS data with point densities >10 pts/m to determine the number and vertical extent of canopy layers (canopylayer, canopylength), seasonal variations in the topmost canopy layer (canopytype), as well as small-scale heterogeneities in the canopy (canopyheterogeneity). We first tested and developed the method on a small forest patch (800 ha) and afterwards tested transferability and robustness of the method on a larger patch (180,000 ha). We validated the approach using an extensive set of ground data, achieving overall accuracies >77% for canopytype and canopyheterogeneity, and >62% for canopylayer and canopylength. We conclude that our method provides a robust characterization of canopy layering supporting automated canopy structure monitoring. OPEN ACCESS Forests 2015, 6 4147


Introduction
Forests provide a wide range of valuable ecosystem goods and services [1].Understanding forest ecosystems and their underlying processes help with forecasts using scenarios under global change conditions and supports the establishment of sustainable management strategies [2][3][4][5].The three-dimensional distribution of geometric objects and their topology within a forest canopy (termed "canopy structure" [6,7] is considered a particularly crucial constituent of forest ecosystems' functioning and processes.Canopy structure can, among other things, serve as an indicator of forest stand resistance to disturbances [8,9], enable the estimation of the conservation potential for biodiversity [10,11], or allow the identification of recruitment limitations [12,13].Canopy structure itself is not a measurable quantity, but the properties of canopy structure can be described by means of a wide variety of structural variables, such as tree height, tree diameter distribution, or canopy layering [14].The assessment of canopy layering, in particular, is of great importance.Canopy layers provide essential structural components of wildlife habitat [15], affect stand productivity [16], and are used to manage forest stands and to identify forest structural stages [17].From the various existing definitions of canopy layer (cf.[18]), we follow Franklin and Spies [19] and understand canopy layer as continuous vertical distribution of foliage within the canopy.
Traditionally, canopy layering is assessed by conventional fieldwork in relatively small areas, which is time-consuming and occasionally subjective [20,21].Advances in Earth observation systems and analysis techniques have greatly improved the ability to characterize canopy structure over large areas in not only the horizontal but also the vertical dimension [22][23][24].In particular, light detection and ranging (LiDAR) systems, especially airborne laser scanning (ALS), are suitable to provide very detailed vertical and horizontal information on canopy structure based on the physical measurement principle of active sensing and full-waveform digitization [25][26][27][28].This also means that definitions of canopy layering which, for example, focus more on composition of life forms or age classes cannot be taken into consideration, as the types of physical measurement possible with ALS do not meet these definitions.
ALS-based assessment of canopy layering can be either derived from LiDAR waveform data (mostly the case for large or medium-footprint LiDAR systems) [29][30][31], or based on ALS point clouds directly [32,33] or indirectly, e.g., by generating a kind of synthetic waveform using classified echo height distributions or height distribution probability functions [34][35][36].The assessment of canopy layering with ALS is mainly performed using area-based approaches (ABA).The spatial scale of ABA is currently often driven by user requirements and usually limited by the features present in the available data [37].However, the choice of a specific spatial scale affects the canopy layer analysis, as a decrease in horizontal resolution will result in more mixing of horizontal and vertical structure elements within each spatial unit [34,35,38,39].Most canopy-structure components, however, have inherent spatial scales and the choice of the spatial scale should be made considering the investigated structural component [40,41].Nevertheless, an acceptable compromise often needs to be found between the technical capabilities of ALS and user needs, such as the application domain or cost-benefit considerations [42,43].In this study, ABAs represented in regularly spaced grids were used, allowing flexible analysis in terms of spatial scales and cross-scale comparisons (e.g., testing reproducible up-and down-scaling methods).
ABA's were investigated by Wulder et al. [25,44], White et al. [45], and Naesset [46] who found the concept to be sufficiently robust for operational use.The assessment of canopy layering, however, is rarely considered.Existing approaches to assess canopy layers often require a substantial set of prior information about stand characteristics or rely on pre-defined height thresholds for an existing stratification of the canopy.Whitehurst et al. [30], for example, on one hand used the cumulative canopy cover and pre-defined height thresholds to distinguish between understory (0 to 5 m), midstory (5 to 15 m), and overstory (>15 m).On the other hand, the peaks of the LiDAR waveform were used to determine canopy layers, assuming that each peak represents the mean layer height.The peak definition parameters, however, were determined by visual analysis [30].Ferraz et al. [32] applied a mean shift segmentation algorithm to determine the modes of the distribution functions of echo heights and afterwards defined thresholds to divide the canopy into overstory, understory and ground vegetation.In this context, the definition of these thresholds remains critical as one component of the algorithm is linked to a function of canopy layer depth and thus depends on the tree species and/or forest biome [32].Wang et al. [35] also used the modes of echo height distribution functions to divide the canopy into single-layered and double-layered canopies.Furthermore, the length of the determined layers was estimated, whereby only changes in absolute point density (i.e., inflexion points in the height distribution function) were used as the criterion to define the borders of a layer, which are not necessarily coincident.Maltamo et al. [36] also distribution functions of echo heights, in a manner comparable to Whitehurst et al. [30], Ferraz et al. [32], and Wang et al. [35], to relate the modes of the function (i.e., unimodal or multimodal) to the amount of canopy layers.In contrast, however, a fundamental part of their algorithm was the prior definition of a threshold to separate the dominant tree layer and understory trees, which limits the determination of canopy layers to a maximum of two layers.Morsdorf et al. [33] applied a cluster analysis to echo heights and echo intensities to stratify the canopy into three layers, whereby each layer consists almost entirely of oak, pine or shrubs, respectively.Although the proposed approach results in high accuracy values, a successful transfer of the method to dense, mixed temperate forests is unlikely.Zimble et al. [43] used a LiDAR-based canopy height model for the identification of individual trees and to determine their respective heights.The variances in tree heights were subsequently analysed on a 30 × 30 m cell size and used to distinguish between single-layered and multi-layered canopies.Accordingly, only trees forming the top canopy layer were considered and the resulting canopy stratification is more related to the composition of age classes and less to the actual vertical foliage distribution.Possible canopy layers below the top canopy layer, however, could not be assessed.
Therefore, most of the existing approaches are limited in their transferability to other sites and they tend not to be directly comparable.The inherent spatial scales of canopy-structure components are barely considered and, to our knowledge, a multi-scale approach regarding canopy layering has so far only been presented by Wang et al. [35].To overcome these limitations, we propose using an automated and transferable ABA to provide quantitative assessments of canopy layering at different scales.

Aim and Research Objectives
The aim of this study is to develop a transferable, grid-based assessment of canopy layering.Requirements are defined using a stakeholder approach, originating from the forestry department of the Canton of Aargau, Switzerland, a regional authority responsible for forest management.The requirements were to assess canopy layering, canopy length, canopy type and small-scale heterogeneity of the canopy using an automatic and transferable approach with a target grid-cell size of 10 m × 10 m).Accordingly, the main objectives of this study were (i) to determine and classify the number and vertical extent of canopy layers; (ii) to distinguish between deciduous and evergreen forest stands; (iii) to characterize the small-scale heterogeneity of the canopy; and iv) to test the robustness and transferability of the method to a significant portion of forested land within the Canton of Aargau.The derived canopy-structure descriptors (canopylayer, canopylength, canopytype and canopyheterogeneity) as well as their spatial distribution were validated using an extensive set of field data.

Study Site
We developed and tested the method at the Laegern site (47°28′ N, 8°21′ E), a temperate forest located in Canton of Aargau, Switzerland (Figure 1).The site is approximately 800 ha and mainly composed of mixed deciduous forest, dominated by common beech (Fagus sylvatica L.), European ash (Fraxinus excelsior L.), and sycamore maple (Acer pseudoplatanus L.), with scattered silver fir (Abies alba Mill.) and Norway spruce (Picea abies L.) trees [47,48].To prove the transferability and robustness of the approach, the method was subsequently transferred to the entire area of Canton of Aargau (180,000 ha) without adapting the algorithm, including areas undergoing different forest management practices (ranging from natural forests to highly intensive regimes with silvicultural interventions) and a variety of forest types.

Reference Data
An extensive set of ground based reference measurements was collected for the Laegern site between 2012 and 2014.We applied a stratified sampling, where strata correspond to development stages (age classes) and canopy type (evergreen/deciduous), resulting in 200 sample plots with an approximate size of 100 m 2 per plot.Plot-wise measurements were taken for forest structure information on the individual tree level, such as tree position, social position, crown dimensions, tree vitality, tree height, tree species, and diameter at breast height for all trees with a diameter at breast height (DBH) larger than 20 cm (cf.[34]).The tree positions were precisely measured based on terrestrial land surveying using a total station and a differential GPS system, linked to fixed points of the Swiss national survey.Tree species were afterwards aggregated into the canopy type's "deciduous canopy" and "evergreen canopy".
Forest experts interpreted for an area within 10 m around each tree position's canopy layering and estimated the canopy length based on hypsometer measurements.In this regard, it must be noted that for canopies with two or more layers, it is assumed that measurements of the canopy length become increasingly uncertain.The small-scale heterogeneity was visually assessed in consideration of variations in DBH, crown base and the canopy length.Even if foresters' interpretation of canopy structure is often subjective and depends on the applied definitions (cf.[18]), it is the only non-destructive solution so far known to us to assess these types of canopy-structure descriptors.The information about canopy layering, canopy length, canopy type, and small-scale heterogeneity was transferred to a grid-based representation with a 10 m × 10 m grid-cell size.
For the entire forested area of Canton of Aargau, stand maps were available as polygon layers [49].The stand maps were produced by regional foresters, containing detailed descriptions of canopy type, canopy layering and canopy length as well as current forest management practices.To ensure the spatial and temporal comparability of the stand maps with the ALS-derived canopy structure information, we used only stand polygons updated in the years 2012-2014.Reference data thusly derived contain more than 16,000 stand polygons, distributed over the entire area of Canton of Aargau.

Airborne Laser Scanning Data
In the approach presented in this study, we used full-waveform ALS data.The ALS data for the Laegern site were acquired in 2010 under both defoliated (or leaf-off) and foliated (or leaf-on) conditions.The ALS data for Canton of Aargau were acquired in 2014.The used sensor specifications are summarized in Table 1 and described more in detail in [50,51].The benefit of full-waveform data is the approximation of the entire backscattered signal by digitization, which facilitates the extraction of additional features of each reflecting object within the ALS footprint.To detect and extract specific object reflections, Gaussian pulse estimation was applied using Riegl's software RiANALYZE 3.0.8(RIEGL, Horn, Austria) [52] in order to obtain representative echo descriptions (cf. Figure 2).This includes the derivation of a three-dimensional point cloud composed of planimetric coordinates (x and y), ellipsoidal heights (z) and echo type (i.e., first, intermediate, last echo).Additionally, it enables the physical description of each reflector with information such as amplitude, width and intensity of each specific echo.For the ALS data of the Laegern site, the terrain-corrected point cloud was derived by extracting potential ground returns from the point cloud while retaining single and last echoes, together with their geometrical characteristics and the corresponding echo-width information.Detection and filtering of the ground returns was carried out using an adaptive multi-scale filter [53].A 0.5 m × 0.5 m digital terrain model (DTM) was generated from the remaining points using an ordinary kriging interpolation.For the ALS data of Canton of Aargau, a DTM with a spatial resolution of 0.5 m × 0.5 m was processed by MILAN Geoservice GmbH using the Terrasolid software (Terrasoild, Helsinki, Finland).For both data sets, the height above ground was then calculated for each echo of the point cloud by subtracting the interpolated DTM value at the corresponding horizontal echo location.For the subsequent use of the point clouds, we used only the planimetric coordinates and terrain-corrected heights.

Stakeholder's Requirements
The thematic framework for our research was established by stakeholder requirements.According to these requirements, we defined canopy-structure descriptors and appropriate descriptor classes (Table 2).Small-scale variation of vertical foliage distribution within each 10 m × 10 m grid cell.Significance of variation should be defined using the p-value with a significance level of 0.05.Two classes should be distinguished: 1. homogeneous canopy (non-significant variation) 2. heterogeneous canopy (significant variation)

Calculation of the Relative Frequency Distribution
One of the essential feature of our method is the relative-frequency distribution (RFD) of echo heights within a given grid cell.This means that all echo returns within a grid-cell were binned according to height above ground [54], whereby the bin size was set to 1 m.Thus, the RFD shows the echo occurrences (echobin) for all height intervals/bins (γ) up to the maximum echo height (γN) as a percentage of the total number of echoes (echototal) (Equation ( 1)).
Those RFDs can be displayed as relative-frequency histograms (often called pseudo-waveforms in the literature) and should comprehensively represent the ALS-based characterization of the vertical canopy structure, although this is not necessarily a complete representation of the actual canopy structure (cf.[39,55]).Nevertheless, previous studies emphasized and proved the potential of RFDs for comprehensive canopy-structure characterization [35,54,56,57].
For our study, the definition of a grid-cell size of 10 m × 10 m to derive RFDs was according to the dimensions of the reference data and the user requirements, which is a common practice (cf.[25,45]).In a previous study [34], however, we determined based on a data-driven approach that the best feasible grid-cell size was 8 m × 8 m for canopy-structure characterization.Muss et al. [54] also concluded that a grid-cell size between 10 and 15 m was optimal to assess forest-canopy structure using RFDs.Therefore, the applied user-defined grid-cell size in our case was in line with "optimal" grid-cell sizes determined in recent research findings.
To assess the small-scale heterogeneity, we defined a small-scale grid-cell size of 1 m × 1 m and derived RFDs in the same way as described before.The definition of the fine-scale grid-cell size was deduced from the positional accuracy of the point cloud, the point density and the footprint diameter of the laser pulse (cf.[58]).As we showed in a previous study [34], the point density was not a limiting factor: Even for the small grid-cell size of 1 m × 1 m, the available point cloud has almost the maximum achievable information on the canopy structure that current ALS systems can provide.This means that a further increase in point density would not significantly change the derived RFDs.
Following Lefsky et al. [31], a threshold value of 1% was used to classify each RFD bin into either "filled" (relative frequency ≥ 1%) or "empty" (relative frequency < 1%).Filling single, isolated empty bins smoothed the derived RFD.Isolated empty bins can be either caused by smaller gaps in the canopy (below the stakeholder's definition of 3 m in vertical extent, cf.Table 2) or by occlusion effects and are thus not considered large enough to separate independent canopy layers.The same approach was used to remove isolated filled bins (below the stakeholder's definition of 3 m in vertical extent, cf.Table 2).Finally, the RFD was converted to a vector RFDvec with "1" for filled bins and "0" for empty bins.Figure 3 shows the processing from the terrain corrected ALS point cloud to the filtered RFDs, which were used for the subsequent calculation of canopy-structure descriptors.

Calculation of Canopy-Structure Descriptors
For the assessment of canopy layering (canopylayer), we first investigated the added value of the multi-seasonal ALS data.Kükenbrink et al. [59] investigated occlusion effects based on a 1 m voxel grid, showing that particularly for ALS acquisition under leaf-on conditions, large parts of the canopy volume were occluded.Therefore, we derived the RFDs of the leaf-on/leaf-off ALS data separately and correlated them to the RFD based on the combined point cloud (cf.[34]).It turned out that RFDs based on the leaf-on data show a high correlation of >95% to their respective histogram based on the combined point cloud, whereas the leaf-off derived RFDs mainly contain redundant information about the canopy layering (cf.[34]).This may indicate that the occluded volume in the leaf-on acquisition is often "empty" space without structural components.Thus, the underestimation of foliage and foliage distribution would be less significant than expected.Moreover, in particular for deciduous canopies, the leaf-off RFDs more likely represent stems and branches instead of actual foliage, which does not meet the definitions of canopy layering we used.We thus decided to use the leaf-on RFDs only to assess canopy layering, canopy length and the small-scale heterogeneity of the canopy.
Starting from the topmost filled bin, we applied a convolution kernel K to RFDvec, which determines vector positions with changes from "0" to "1" and vice versa (cf.Equation ( 2)).According to K, changes from "0" to "1" result in "1" and indicate the start of a layer, whereas changes from "1" to "0" result in "−1" and indicate the end of a layer.Accordingly, the sum of positive changes ("1") is related to the amount of canopy layers (canopynumber).
According to the requirements defined in Table 2, we classified canopynumber into canopylayer with the classes "1-layered", "2-layered" and "multi-layered" canopy.The length of the topmost canopy layer was subsequently determined in RFDvec by the sum of filled bins, starting from the topmost filled bin to the first change from "1" to "0".Following, we calculated the canopy ratio as the length of the topmost canopy layer divided by the total canopy height.To obtain canopylength, the canopy ratio was classified into "short canopies" (canopy ratio < 0.5) and "long canopies" (canopy ratio ≥ 0.5).
Using both the RFDs based on leaf-off and leaf-on data, we calculated the differences in relative frequencies in the topmost canopy layer as proxy for canopytype.We classified the difference values into significant differences and non-significant differences using the p-value with a significance level set to 0.05.To test the statistical hypothesis, we applied a Welch's test [60].Significant differences, which are mainly caused by changes in the leaf volume, indicate a "deciduous canopy", whereas non-significant differences indicate an "evergreen canopy".
The small-scale heterogeneity of the canopy (canopyheterogeneity) was derived using the RFDs derived on 1 m × 1 m (RFDsmall) and the corresponding RFD derived on 10 m × 10 m (RFDlarge).Each RFDsmall was compared to the specific RFDlarge using the empirical correlation coefficient rxy (e.g., [61]).The average rxy for the respective 10 m × 10 m grid cell was afterwards used to describe canopyheterogeneity.Comparable to the applied classification of canopytype, the classification into homogenous canopy ("non-significant variation") and heterogeneous canopy ("significant variation") was carried out based on the p-value with a significance level set to 0.05.The statistical hypothesis was tested using the t-test [60]."Significant variations" reveal a high variability in RFDsmall within a 10 m × 10 m grid cell and thus point to a worse representation of the actual vertical structure by RFDlarge in contrast to a low variability of RFDsmall, indicated by "non-significant variations".
The resulting canopy-structure descriptors canopylayer, canopylength, canopytype, and canopyheterogeneity were validated on the Laegern site with the forest inventory data for the canopy-structure descriptors based on both the 2010 and on the 2014 data.For the entire Canton of Aargau, we validated canopylayer, canopylength, and canopytype only, as until now canopyheterogeneity has not been a used/defined variable for canopy-structure characterization at the stand level.For the applied accuracy assessment, we calculated the widely used error matrix metrics "overall accuracy", "user's accuracy", "producer's accuracy", and the "Kappa coefficient" [62].Overall accuracies are computed by dividing the total correct classified grid cells by the total number of grid cells for each canopy-structure descriptor.Producer's accuracies indicate the probability of a "reference grid cell" being correctly classified and user's accuracies indicate the probability that the classified canopy-structure descriptor represents the "reference grid-cell".According to Rosenfield & Fitspatrick-Lins [63], the Kappa coefficient is the proportion of agreement after chance agreement is removed.
To estimate the comparability of the field measurements done for the Laegern site with the information contained in the stand polygons, we validated the canopy descriptors based on the 2010 and the 2014 data independently and each with both sources of reference data.

Results
For the Laegern site, the validation shows promising results for all canopy-structure descriptors.The validation results for the canopy-structure descriptors are listed in Table 3, including the differentiation between the ALS and reference data used.The different ALS and reference data used influence the accuracy values of the derived canopy-structure descriptors only to a small degree.Nevertheless, all canopy-structure descriptors derived from the ALS data of 2010 show higher accuracies except for canopytype.On average, the validation with the forest inventory data of the Laegern site results in slightly higher accuracies in comparison to the validation results based on the stand polygons of the Canton of Aargau.With 2.4% (cf.Table 3), however, the mean difference in overall accuracies caused by the used reference can be more or less neglected.The same applies for differences in the overall accuracies caused by the different ALS data used with a mean difference of 3.1%.Particularly in terms of the determination of canopytype and canopyheterogeneity, the developed algorithm shows a good performance.The classification of canopytype leads to high overall accuracies up to 89.5%.The highest overall accuracy for canopyheterogeneity is reached with 81.5%.Overall accuracies up to 66.0% are achieved for canopylength, whereby misclassifications mainly occur in 2-or multi-layered forest stands.Canopy layering (canopylayer) is derived with an overall accuracy of 67%, ranging from 80% (1-layered canopies) to 53.3% (multi-layered canopies).The highest overlap occurs between the 2-and multi-layered canopy classes.In general, most classification uncertainties are caused by continuous transition zones between forest stands differing in age and structure.For these regions an accurate separation (e.g., in terms of a definition of a borderline between a 1-layered and a 2-layered canopy) is complex, even more so in the field.Table 3. Accuracy assessment of the canopy-structure descriptors for the Laegern site.We use the accuracy metrics "overall accuracy" (OA), "user's accuracy" (UA), "producer's accuracy" (PA), and the "Kappa coefficient" (K).Resulting values in italics represent results obtained from stand polygons.The validation for the entire Canton of Aargau shows comparable results to the one of the Laegern site, with slightly higher overall accuracies.For canopytype we achieve an overall accuracy of 90.7%.Canopy layering (canopylayer) is derived with an overall accuracy of 69.2%.An overall accuracy of 70.3% is obtained for canopylength.The detailed validation results are given in Table 4. Table 4. Accuracy assessment of the canopy-structure descriptors for the whole of Canton of Aargau.We use the accuracy metrics "overall accuracy" (OA), "user's accuracy" (UA), "producer's accuracy" (PA), and the "Kappa coefficient" (K). Figure 5 shows a larger subset of canopylayer as an example of the derived canopy-structure descriptors for Canton of Aargau.In total for Canton of Aargau, 61% of the forested area is characterized by a 1-layered canopy, 38% by a 2-layered canopy, and only 1% by a multi-layered canopy.81% of the forested area is covered with "deciduous canopies" and thus 19% with "evergreen canopies".Long canopies occur in 88% of the forested area.

Discussion
The analysis of leaf-on/leaf-off conditions on the RFDs showed that the added value of leaf-off data is relatively low if leaf-on data are available, in particular in terms of characterizing vertical foliage distribution which does not include the wooden elements in the forest (e.g., stems and branches).This finding, however, largely depends on the combination of the underlying grid-cell size and the respective point density and only refers to the impact on the RFDs.It is important to stress this because in contrast to our findings, other studies reported an added value of leaf-off acquisitions over leaf-on acquisitions for area-based canopy-structure estimates (cf.[64][65][66]).Moreover, Kim et al. [67] pointed out that leaf-off data only is also suitable to discriminate between evergreen and deciduous canopies, which can be an important issue in terms of decreasing costs for ALS-based forest monitoring.Nevertheless, the ALS data used in these studies, had much lower point densities, and the direct comparability with their results is therefore limited.
In terms of the point density, we also need to point out that the structural information expressed in the RFDs depends on the point density in connection with underlying grid-cell sizes [68].Lim et al. [69] raised the issue of whether or not different point densities influence canopy-structure characterization, which is a contentious one.In general, RFDs based on lower point densities and calculated on smaller grid-cell sizes are difficult to interpret as a possible undersampling of the actual canopy structure increases.With increasing grid-cell sizes, however, the relevance of the point density on the RFDs decreases (cf.[34]).We showed in [34] that, for example, for a grid-cell size of 10 m × 10 m and a point density of 10 pts/m 2 , the derived RFDs match to ~95% with RFDs based on much higher point densities (i.e., >60 pts/m 2 ).A decrease in the point density to 5 pts/m 2 or 1 pts/m 2 decreases the match to 90% and 74%, respectively.For smaller grid-cell sizes, for example 5 m × 5 m, the match of the RFDs based on point densities of 5 pts/m 2 or 1 pts/m 2 becomes increasingly poor, at 85% and 59%, respectively.Whether this loss of information is significant or important still depends on the research question and/or the planned application.For example, Wilkes et al. [70], Hawbaker et al. [71], and Treitz et al. [42] pointed out that lower point densities (<5 pts/m 2 ) are sufficient for canopy-structure analysis regarding tactical forest management, but these findings are only valid for canopy-structure descriptors derived on larger scales than the one we used in this study.Hayashi et al. [72], in contrast, concluded that in complex forests, low-density ALS data was ineffective in estimating the defined canopy-structure descriptors with sufficient accuracy.This is in accordance with our results and highlights the need and advantage of higher point densities.We thus recommend point densities >10 pts/m 2 , particularly when the small-scale assessment of the vertical canopy structure is also of interest.Although ALS data with such high point densities are increasingly available for many Cantons in Switzerland, we are aware that this is not the standard.We are nevertheless convinced that in future nationwide ALS data with high point densities >5 pts/m 2 will be more and more available.For the second national height map of the Netherlands (AHN2), for example, LiDAR data with an approximate mean point density of 9 pts/m 2 were acquired already in 2011/2012 [73].
What should be also taken into account is a possible overestimation of upper canopy parts or an undersampling of lower canopy parts when using RFDs.This could be addressed by, for example, a weighting of echoes according to return number or intensity.However, as we only used the frequencies to filter the RFDs (i.e., determining filled and empty bins), the influence on the resulting canopy-structure descriptors is judged to be very small.Nevertheless, for canopy-structure descriptors which are more related to bio-physical variables, such as leaf area density, a proper weighting of the echoes should be taken into account [74].
The validation of the canopy descriptors proved good accuracies, particularly in terms of the seasonal and horizontal canopy variety (canopytype and canopyheterogeneity).The seasonal canopy variety was very well covered with canopytype and the reached overall accuracy is comparable to the highest obtained accuracies in existing studies (cf.[67,75,76]).Nevertheless, the validation of canopytype for Canton of Aargau could be influenced by misclassification of deciduous conifers, such as European Larch, as the reference data contains the distinction between coniferous and deciduous trees only.As European Larch covers less than 2.5% of the total area, the overall accuracies were likely only marginally influenced.The validation of the vertical canopy stratification (canopylayer and canopylength) showed that misclassifications mainly occurred for forest stands which are very heterogeneous, such as multi-layered canopies with a long topmost canopy layer.For these forest stands, occlusion by the dense vegetation layer hinders a full sampling of the lower canopy parts even for areas with point densities >40 pts/m 2 available.Moreover, for most applications in forest inventory/management or forest ecology, there is an important difference if a second layer below the topmost canopy layer belongs to the forest floor/ forest understory or still belongs to the tree canopy.Since, based on the proposed method, a distinction between these types of layers would be possible, it should be included in following investigations of canopy layering.
The plot-wise accuracy assessment of the canopy-structure descriptors, however, turned out to be very difficult.The subjective visual evaluation of canopylayer, canopyheterogeneity, and canopylength by the forestry experts includes a source of uncertainty, in particular for canopies with more than 2 layers (without a possible forest floor layer).Thus, this information cannot be regarded as an error-free reference: rather, it should be regarded as a data source for cross-comparison and feasibility checking.A possible alternative to create a more objective reference could be the modelling of canopy structure as presented by Scanlan et al. [77], using field measurements and allometric relationships.Furthermore, the canopy stratification approaches used in the forest inventory are more related to the composition of different tree development stages and less focused on the actual continuous vertical foliage distribution.For example, echoes from the vegetation in the lower canopy parts can be either from the forest floor or young growth or caused by low branches of larger trees.For the ALS-based canopy stratification we applied, it is difficult to distinguish between these different sources, even if a combination with methods at the individual tree level could help to address this issue, as, for example, understory and understory trees can be partly detected [32,78,79].
The assessment of the small-scale heterogeneity of the canopy (canopyheterogeneity) showed that specific forest plots are characterized by a higher horizontal variability of vertical canopy structure than others and thus can be used to further distinguish between forest types.The spatial unit defined in most forest inventories to assess the heterogeneity of canopies, however, is with respect to a larger spatial unit (1 ha) than the one we used in this study (10 m × 10 m).It turned out that this change in the spatial unit made it difficult for the forestry experts to apply the usual approach to assess heterogeneity.Although the consideration of the small-scale heterogeneity thus may be not appropriate for forest inventory, these small patches within the forest serve as important indicators for biodiversity estimation or habitat assessment [80].Therefore, if the available ALS data source facilitates a small-scale analysis it is useful to include LiDAR metrics derived on the smallest grid-cell size possible.Otherwise, a loss of spatial information regarding vertical canopy structure will occur.
The transferability of the method was proved by reaching comparably good results, at least for a similar type of temperate forest.The slightly better results for the canopy-structure descriptors derived for the whole of Canton of Aargau could be explained by the, on average, lower complexity of forest stands in comparison to the forest stands on the Laegern site.What still needs to be explored is the transferability of the method to forests with a highly complex canopy structure such as in the tropics or to forests with a very open/sparse canopy.In complex tropical forests it is likely that the density of the canopy hinders a full sampling of the vertical structure and the RFDs therefore cannot represent the actual canopy structure completely.In very open forests, in contrast, the grid-cell sizes need to be large enough to cover and describe the openness of the forest.In this context, more advanced self-calibration methods to derive grid-cell size in automated manner from the forest type and ALS data available need to be developed (cf.[34]).
For further applications, additional canopy-structure descriptors can be defined and/or adapted according to user's requirements or, to comply with existing definitions and spatial units, can be aggregated to spatial scales used in the specific forest inventory system, such as the stand level.

Conclusions
In this study we present a method to derive grid-based characterization of canopy layering using ALS data.The resulting canopy-structure descriptors meet stakeholders' requirements in terms of spatial resolution and information content.We demonstrated the transferability of the method from a small, well-characterized forest patch (800 ha) to a larger patch (180,000 ha), covering all forested area of the Canton of Aargau.Furthermore, we evaluated the advantages of using different grid-cell sizes and temporal information contained in multi-seasonal ALS data for further distinction and characterization of canopy layers.Our method to determine canopy-structure descriptors supports and improves operational, ALS-based canopy-structure monitoring in terms of inter-annual comparability.
Nevertheless, some limitations still need to be considered: (i) objective methods for field-based, non-destructive assessments of vertical foliage distribution are necessary in order to perform a more meaningful accuracy assessment; (ii) the transferability of the proposed method to very complex or open/sparse canopies remains unknown to date; and (iii) for the assessment of the small-scale structure, only ALS data with point densities >10 pts/m 2 are appropriate to obtain credible information.Additionally, proper interpretation of the canopy-structure descriptors still requires substantial forestry expert knowledge.Moreover, there is a need to link ground measurements and ALS-based canopy-structure descriptors considering differences in the used definitions and spatial scales.
The definitions of "significant/non-significant" in terms of temporal variations (canopytype) and small-scale heterogeneity (canopyheterogeneity) of the canopy still need to evolve, but can be used to represent stakeholder requirements, such as their use in forest inventory or forest ecology applications.We see the potential of our approach, however, in the field of forest ecology, for example wildlife habitat assessment (cf.[15]), rather than in the implementation in traditional forest inventories.The next steps will be to transfer the method to a broader variety of forests (including boreal and/or tropical forests) to ensure consistent application of the method.Once implemented, it can serve as a proxy for ecosystem-service-based assessments in forests, which includes the long term goal of monitoring forest change with high fidelity and accuracy.

Figure 1 .
Figure 1.Switzerland (left) and the location of Canton of Aargau (right, top panel).Laegern experimental site (800 ha forested area) plotted against ALS-derived canopy height (bottom panel).

Figure 2 .
Figure 2. Physical concept of airborne laser scanning: interaction of the emitted laser pulse with the vegetation canopy (left panel), the full and the decomposed waveform (middle panel), and the resulting echo descriptors.

Figure 3 .
Figure 3. Processing steps from the ALS point cloud to the filtered RFDvec for a selected 10 m × 10 grid-cell.

Figure 4
Figure4shows resulting maps of canopy-structure descriptors for a subset of the Laegern site.

Figure 4 .
Figure 4. Canopy layering (canopylayer) and small-scale heterogeneity of the canopy (canopyheterogeneity) (top panel).Length of the topmost canopy layer (canopylength) (middle panel).Evergreen and deciduous canopies (canopytype) (bottom panel).All panels are a subset of the Laegern site, with a hill-shaded digital terrain model as background.

Figure 5 .
Figure 5. Canopy layering (canopylayer) scaled to a larger part of Canton of Aargau.Background data are from a hill-shaded digital terrain model.

Table 1 .
Used sensor specifications of RIEGL's airborne laser scanner LMS-Q560 and LMS-Q680i (RIEGL Laser Measurements Systems GmbH, Austria) for the Laegern site and Canton of Aargau.

Table 2 .
Stakeholder's requirements for a grid-based assessment of canopy layering and according canopy-structure descriptors used in this study.