Deriving Stand Structural Complexity from Airborne Laser Scanning Data—What Does It Tell Us about a Forest?

: The three-dimensional forest structure is an important driver of several ecosystem functions and services. Recent advancements in laser scanning technologies have set the path to measuring structural complexity directly from 3D point clouds. Here, we show that the box-dimension (D b ) from fractal analysis, a measure of structural complexity, can be obtained from airborne laser scanning data. Based on 66 plots across di ﬀ erent forest types in Germany, each 1 ha in size, we tested the performance of the D b by evaluating it against conventional ground-based measures of forest structure and commonly used stand characteristics. We found that the D b was related (0.34 < R < 0.51) to stand age, management intensity, microclimatic stability, and several measures characterizing the overall stand structural complexity. For the basal area, we could not ﬁnd a signiﬁcant relationship, indicating that structural complexity is not tied to the basal area of a forest. We also showed that D b derived from airborne data holds the potential to distinguish forest types, management types, and the developmental phases of forests. We conclude that the box-dimension is a promising measure to describe the structural complexity of forests in an ecologically meaningful way.


Introduction
The spatial structure of forests is of great interest to forest scientists as it is related to many ecosystem functions and services [1][2][3][4][5]. Great advances have been made in addressing the three-dimensional (3D) character of forest structure by approaching it with airborne remote sensing as well as ground-based (close-range) remote sensing techniques. Among these techniques, radio detection and ranging (RaDAR; e.g., [6][7][8]), airborne light detection and ranging (LiDAR; e.g., [9,10]), spaceborne LiDAR [11], terrestrial laser scanning (TLS; e.g., [12][13][14]), aerial and satellite imagery (photogrammetry; e.g., [15,16]), or structure-from-motion (SfM, e.g., [17]) can be named as prominent examples. All of them aim at capturing the real forest and representing it based on 3D data in a virtual model space, probably most commonly in the form of 3D point clouds or digital elevation models. Quite recently, the unique opportunity to comprehensively address the detailed 3D structure of something as complex in structure as a tree or forest, gave rise to research focusing on the relationship between the 3D-structural complexity and ecosystem functions and services. Microclimate regulation (e.g., [18]), carbon storage (e.g., [19]), productivity (e.g., [20]), and habitat provisioning (e.g., [21]) are some examples. In the past, such investigations had to rely on surrogates of 3D structure, such as the diameter distribution of trees (e.g., [22]) or indices of structural complexity that are based on tree heights and tree positions, e.g., the structural complexity index [23], the enhanced structural complexity index [24], or the Clark and Evans index of aggregation [25]. For a comprehensive review of structural indices, the interested reader is referred to [26] and [4].
With technological advancements, the latest efforts have aimed at using the full potential of 3D data from high resolution, ground-based point cloud representations across scales, from single trees [27][28][29][30], forest layers [31,32], or complete stands [12,[33][34][35]. Addressing "structural complexity" based on 3D models of the real forest offers new research opportunities. Therefore, we need to be clear about how we define structural complexity.
The structural complexity of a tree can be defined as a summarizing term describing all dimensional, architectural, and distributional patterns of a tree's organs at a given point in time [14]. For a forest, structural complexity may be defined as all dimensional, architectural, and distributional patterns of plant individuals and their organs in a given forest space at a given point in time. While high-resolution 3D data from ground-based approaches proved suitable for deriving different measures of structural complexity on the tree and stand level, the potential of airborne data has received lesser attention [9,[36][37][38]. Airborne data would be particularly interesting though, since large-scale measurements from the ground are still tedious, even if ground-based LiDAR is used in the efficient single-scan sampling mode [35]. Some of the more recently introduced measures, like the stand structural complexity index [12], are specifically designed to make use of the ground-based perspective and are hence not suitable for application to data derived from airborne perspectives. However, a promising approach lies in the use of fractal analysis, more precisely the box-dimension (abbrev: D b ; cf. [29]). This measure can be determined for point clouds from all sources and platforms, including mobile, airborne, or even spaceborne. The D b can theoretically be calculated for any kind of object (cf. [39]) and for single trees as well as stands [13]. D b is considered a measure that integrates the distribution and density of material in space or, in other words, a measure that characterizes the way in which plants physically occupy space [40]. To calculate the box-dimension of forest stands, the laser scanning point cloud is to be transferred into voxel models with voxels of varying size. Usually, one starts with the minimum bounding cube (box) enclosing the entire point cloud as one large voxel representing the stand. Of course, this representation is so simple that it does not contain any other information beyond the extent of the object investigated. However, in subsequent steps the voxel size is reduced and structures are resolved in more and more detail. With each step, both the voxel size and the corresponding number of voxels are recorded and later contrasted with each other (see Methods chapter for details). In this way, in contrast to earlier measures of structural complexity, the D b summarizes the mathematical complexity of an entire point cloud by expressing it as a single number [39]. The potential of the application of the box-dimension approach to airborne LiDAR data has, to the best of our knowledge, not been evaluated so far. In our study, we evaluated the suitability of the box-dimension approach as a measure of structural complexity derived from airborne laser scanning data. To test whether the D b obtained this way is actually meaningful, we first tested whether it could successfully discriminate deciduous from coniferous stands and whether it was able to distinguish the different developmental phases of the stands. We then hypothesized that (i) D b would be related to ground-based measures of structural complexity, namely the stand structural complexity index (SSCI [12]) and the structural complexity index (SCI [23]), as well as measures of structural heterogeneity, namely the coefficient of variation of diameters and the effective number of layers (ENL [41]). We do not assume very close relationships between the different indices as two measures of structural complexity hardly refer to the exact same aspect of complexity. In addition, there are different perspectives on the forest plot when using airborne scanning vs. terrestrial scanning or inventory data from the ground. Finally, terrestrial scans are based on the single-scan perspective but contain very high resolutions (mm) while airborne data are generated from hundreds of perspectives during the fly-over but with much coarser resolutions. We do, however, expect them to correlate significantly. We also hypothesized that (ii) there would be a relationship between the easy to observe stand characteristic basal area and D b . Furthermore, we hypothesized that (iii) the airborne box-dimension would be related to stand age and management intensity as well as (iv) to the microclimatic variability in the stands, which is addressed here using the diurnal temperature range.

Study Sites
In our study, we used 66 forest plots, each 100 by 100 m in extent (see Table 1), that were part of the Biodiversity Exploratories project (see [42] or https://www.biodiversity-exploratories.de/ for further detail). Table 1.
Overview of some key characteristics of the study sites with HAI = Hainich; SCH = Schorfheide-C.; ALB = Swabian Alb. For these 66 plots, both airborne laser scanning data and a full all ground-based enumeration of all trees with diameter at breast height (DBH) > 7 cm was available. The plots were in the three Exploratories (see Figure 1), with 16 plots in the Schorfheide-Chorin region (Northeast Germany), 24 in the Hainich-Dün region (Central Germany), and 26 in the Swabian Alb region (Southwest Germany).

Airborne Laser Scanning Data and Point Cloud Processing
Airborne laser scanning data were acquired in July, 2015 using a Q780 Riegl Sensor at an operating frequency of 400 kHz from a flight height of approximately 950 m above ground. This resulted in an average pulse density of 13 pls × m −2 with a mean pulse spacing of 29 cm, as observed for the areas of the 66 plots. As multiple returns were recorded per pulse, the average point density was 36.24 pts × m −2 . The pre-processing was performed using LAStools [43]. This included the removal of isolated returns, retiling into 500 × 500 m tiles, classification of returns into ground and vegetation, and normalization of the elevation values to above ground level (AGL). Finally, the point clouds were clipped using the plot boundaries and the vegetation returns were exported as CSV files with the xyz-coordinates of each return. A visualization of the airborne LiDAR data of two exemplary plots from different perspectives is provided in Figure 2.  Table 1. Colors indicate return intensity were not of importance in our analysis. We kept colors here only for better visualization.
Based on these CSV files, we calculated the box-dimension of each plot for the full 100 × 100 m with Mathematica software (Wolfram Research, Champaign, USA) and as described in the following. Based on all points in the point clouds, the algorithm determined the number of boxes of varying size that were needed to encapsulate all points. Starting with the minimum bounding cube (only one needed to enclose all points), the box sizes were consecutively reduced in steps that always cut the box edge length in half. Consequently, starting with 100 × 100 × 100 m, the next size was 50 × 50 × 50 m, 25 × 25 × 25 m, and so on, until the lower cut-off was reached. The lower cut-off was set as 50 cm as a highly conservative estimate of the point cloud resolution.
The box-dimension of a tree or forest can then be considered the slope of the fitted straight line through a graph of log(N) over log(1/r) [29,30,44]. Here, log() is the natural logarithm, N is the number of boxes (here: boxes) of size r (edge length in relation to the initial edge length used for the first box) needed to enclose all points of a tree or forest. Figure 3 provides a visual outline of the box-dimension calculation.

Terrestrial Laser Scanning and Point Cloud Processing
Terrestrial laser scans were conducted in Summer 2014 on all 66 plots, using a Faro Focus 3D 120 terrestrial laser scanner (Faro Technologies Inc., Lake Mary, USA). On a regular grid with 33 m spacing, we made nine scans per 1-ha plot; for details see [12,41]. All scans were set to scan with an angular resolution of 0.035 • in horizontal and vertical directions. The total 594 single scans were then filtered for erroneous points using standard filters in Faro Scene (Faro Technologies Inc., Lake Mary, USA) and exported as xyz-files. Since the scans were used each on their own (single-scan approach), no merging procedure (co-registration) was necessary. Based on these xyz-files, we calculated the effective number of layers (ENL2D; [41]) as a measure of vertical structure, as well as the stand structural complexity index (SSCI; [12]), using the same algorithms as in the original work and Mathematica software. In short, ENL2D is a measure that first quantifies the number of vertical layers (1 m in thickness) that actually contain vegetation and secondly weighs the numbers based on the filling of the layers using the inverse Simpson index (more plant material in a layer = greater weight). The SSCI index uses vertical cross-sections through the point cloud to quantify the complexity of the surrounding forest scene based on the relationship between the perimeter of the polygon, connecting all points in a cross-section and the area enclosed by it. By using a 1280 cross-section for a 360 • field of view in the azimuthal direction, a detailed description of the structure of the forest scene was provided. The interested reader is referred to [12,41] for more explanation and mathematical details on the two indices.
Furthermore, we used the mean diurnal temperature range (DTR; see [18]; based on Dataset ID: 19007) as a measure of the "stability" of the forest microclimate.

Statistical Analysis
The open source software R (Vers.3.5, R Development Core Team) was used for all statistical analyses. The general relationship between the different variables and the box-dimension (D b ) was described through the Spearman's rank correlation coefficient (rho) as a non-parametric correlation measure, because a linear relationship between the variables could not be assumed. The significance of the correlation was assumed if the p-value of the correlation coefficient was <0.05, indicating the true rho to be not equal to 0.
We used Welch's T-test to test for significant differences in structural complexity between coniferous and deciduous forest plots. The regression analysis between explanatory and response variables was conducted with generalized linear modeling techniques by applying generalized additive models to the data [48]. This method requires no predefinition of the relationship between response and explanatory variable, which allows the unbiased detection of trends in the data themselves. The exponential family distribution of the response variable was specified as normal along with an identity link function. The effective degrees of freedom were limited to a maximum of 3 to avoid over-fitting the data and allow the detection of general trends within the data [49]. The amount of smoothing was chosen automatically through generalized cross-validation [50].

Results
Despite some overlap, we found a significantly (p < 0.05) higher structural complexity (D b ) for the investigated deciduous stands when compared to the coniferous stands ( Figure 4).
We observed significant differences between the management types "age class forests", pooling together all phases (thickets, polewoods, immature and mature stands), and "selection forests" (Figure 5a), considering all 66 plots. The structural complexity of unmanaged forests did not differ significantly from age class forests or selection forests.  Considering only plots dominated by European beech (largest sample), we found an increasing trend in the box-dimension from thickets (mean: 2.04) over that of polewood (mean: 2.07) to immature (mean: 2.11) and mature forests (mean: 2.14) (Figure 5b). We did not test for significance when comparing the age classes due to the limited sample sizes of polewoods (n = 3) and thickets (n = 2).
The box-dimension (D b ) derived from the airborne laser scanning data was significantly positively correlated with the stand structural complexity measures obtained from the ground-based measurements, namely the SSCI and SCI (see Figure 6A+B). It was also positively related to the coefficient of variation of diameters (Figure 6 C) and the effective number of layers (E), both measures describing forest structure. Furthermore, we found a positive relationship with stand age (F). Significant negative relationships were discovered for D b and management intensity (G) as well as the diurnal temperature range (H), the latter indicating an increased microclimatic stability (lower DTR) with increasing complexity. Basal area ( Figure 6C) was not significantly related to the D b from airborne LiDAR. An overview of some numerical characteristics of the generalized additive models that have been visualized in Figure 6 are provided in Table 2.

Discussion
We first tested the performance of D b in differentiating between forest types (coniferous vs. deciduous). D b differed significantly for coniferous vs. deciduous stands even though we observed some overlap in the range of values. We argue that this is to be expected since, aside from forest type, the contrasted stands also differed in management regime, developmental phase, and geographical location. Furthermore, earlier research [12] already revealed that the coniferous stands pooled here have considerable differences in structural complexity depending on the main tree species (pine vs. spruce).
Secondly, the D b from airborne laser scanning (ALS) was evaluated for its performance in separating different management types. The unmanaged forests investigated here showed only intermediate structural complexity and did not differ from the age class stands forest or selection forests, which is a result of the short period since management ceased (~20 yrs). Furthermore, management ceased when the stands were in the optimum phase, which is naturally a phase of low structural complexity, particularly in the case of European beech [35]. The significant differences between age class forests and selection forests supports the hypothesis that the D b from ALS data is indeed capable of distinguishing structural patterns as a consequence of different management. While age class forests have a high variability in structural complexity if all age classes are considered together, they also tend to have the greatest complexity in the younger developmental phase, particularly in thickets [35]. The age class forests investigated here also showed great variability in structural complexity ( Figure 6A) but may seem to contrast when it comes to the complexity of thickets. Instead of the highest values for thickets, we found a trend from low to high complexity for the successive developmental phases. We argue that this is due to the very limited number of shelter trees in our thickets. Usually, the remaining overstory trees add great complexity to thickets [35], but, in the case of the thickets investigated here, these trees have already been removed to a great extent. However, our findings must be interpreted with caution due to the small sample sizes of thickets (n = 2) and polewood (n = 3).
We further compared D b with various ground-based measures, including two measures of structural complexity, namely the stand structural complexity index (SSCI) and the structural complexity index (SCI). While SSCI is a holistic measure of structural complexity that is derived from terrestrial laser scanning (cf. [12]), SCI can be calculated from conventional inventory data. Both have proven useful in the past as surrogates of microclimatic stability (SSCI; [41]), classifiers of forest types (SCI: [23]; SSCI: [12] and [35]), or predictors of biodiversity, e.g., for ant species (SSCI: [51]) or tree species diversity (SSCI: [12]). Here, the birds-eye perspective of airborne laser scanning (ALS) yielded intermediate but significant correlations (R = 0.35 − 0.51) between D b and the two ground-based measures of structural complexity. It is hence not surprising that D b was also correlated with the coefficient of variation of the diameters at breast height (CV_dbh) as the latter is known to be strongly correlated to stand structural heterogeneity [46]. We showed that D b from airborne data, despite the possible effects of the occlusion of lower forest strata, was correlated with several ground-based measures that are related to stand structural complexity. The occlusion effects may have reduced the strength of the correlations but did not impede them. The observed correlation between D b and the effective number of layers (ENL2D), a measure of vertical structural evenness, supports this argument. Morsdorf et al. [52] also highlighted that ALS data can be used to discriminate forest strata even in multi-layered forests.
A relationship between structural complexity and basal area may be expected, since D b is dependent on stand density and material distribution in space [39]. Earlier studies found that basal area was related to Zenner and Hibbs' SCI [53]. In addition, previous research showed that structural complexity measured using SSCI was dependent on stand density as well, even though space-filling was used instead of basal area as a measure of density [12]. However, space-filling considers the actual canopy space exploration and goes far beyond basal area (sensu [34,54]). Here, we argue that basal area was not related to D b since it is not a holistic measure of structural complexity. Basal area provides no information on the crown dimensions (those parts that are responsible for most of the complexity in a forest), but also trees below a certain threshold, usually 7 cm in diameter at breast height, are not included in its calculation based on standard inventory data. In contrast, in the point cloud, all trees are included and accounted for, regardless of their size. The absence of a significant correlation between D b and basal area may be a result of these fundamental discrepancies and leads to a rejection of hypothesis (ii).
We further hypothesized that D b would be related to the age of a stand. In our data, we found clear indications to support this hypothesis. With increasing stand age, structural complexity as measured by D b increased (Figure 6f). We argue that this is due to the increased structural complexity of the individual trees with age [30]. The effect was still evident when the forests were grouped according to the main tree species (data not shown), meaning that this relationship did not result from the fact that the deciduous stands in our study were on average older than the coniferous ones (see Table 1).
The observed negative relationship between the management intensity index (SMI) and D b confirmed the findings of Stiers et al. [35], who studied beech-dominated forests across a management gradient from intensively managed to unmanaged old-growth forests, clearly highlighting that the greatest complexity is to be found in old-growth stands (low SMI). ALS data has already shown a potential for addressing forest structure-management relationships in earlier studies. For example, Valbuena et al. [55] showed that the Gini coefficient derived from ALS data was useful to identify changes in forest structure due to management. Our findings on the relationship between D b and stand age and SMI are in line with these results and support our third hypothesis (iii).
Finally, hypothesis (iv) was also supported by our identification of a significant positive relationship between the microclimatic stability, here addressed using the diurnal temperature range, and the structural complexity (D b ) of the stands. Earlier pilot studies already showed that structural complexity, when addressed from the ground using the TLS-based SSCI, was positively related to the microclimatic stability, reducing daily temperature fluctuations and the vapor pressure deficit [12].
Past research showed that a large number of forest metrics can be derived from ALS point clouds, including species identities [56,57], a plethora of structural measures [9,58], and even forest type classifications [10,59,60]. As a measure of structural complexity, D b was related to all measures tested in our study except basal area. The strength of the relationships may often only be weak or intermediate (R = 0.34 − 0.51), but we argue that this is to be expected if a holistic measure like D b, that uses basically all data in the point cloud to combine them into a single number, is related to the individual characteristics of a stand.
In earlier studies based on digitizers, measurements of D b were restricted to tree compartments, e.g., roots [61]. Terrestrial laser scanning (TLS) data expanded the application of fractal analysis to single trees or tree groups [13,29]. Here, we show that D b from ALS can be a useful descriptor enabling an assessment of forest structural complexity for 1-ha plots and theoretically much larger areas and even in difficult terrain. This supports the idea that applying fractal techniques at the stand level can yield new insights into forest structure [62].

Conclusions
Quantifying the structural complexity of forest stands is important to support a management for complexity, thereby strengthening ecosystem functioning and ecosystem services. We showed that the box-dimension (D b ), as a measure of structural complexity that already showed potential in ground-based LiDAR applications, can also be derived from airborne laser scanning data. Using different ground-based measures for reference, we were able to show that ALS-derived D b has the potential to identify differences in forest type, management type, and developmental phase. Furthermore, it was sensitive to stand age, management intensity, and structural changes. Finally, it was related to an indicator of microclimatic stability in terms of temperature fluctuations in the stands.