Estimation of alpine forest structural variables from imaging spectrometer data

: Spatial information of forest structural variables is crucial for sustainable forest management planning, forest monitoring, and the assessment of forest ecosystem productivity. We investigate a complex alpine forest ecosystem located in the Swiss National Park (SNP) and apply empirical models to retrieve the structural variables canopy closure, basal area, and timber volume at plot scale. We used imaging spectrometer (IS) data from the Airborne Prism EXperiment (APEX) in combination with in-situ measurements of forest structural variables to develop empirical models. These models are based on simple and stepwise multiple regressions, while all potential two narrow-band combinations of the Simple Ratio (SR), the Normalized Difference Vegetation Index (NDVI), the perpendicular vegetation index (PVI), the second soil-adjusted vegetation index (SAVI2), and band depth indices were tested. The accuracy of the estimated structural attributes was evaluated using a leave-one-out cross-validation technique. Using stepwise multiple regression models, we obtained a moderate to good accuracy when estimating canopy closure (R² = 0.81, rRMSE = 10%), basal area (R² = 0.68, rRMSE = 20%), and timber volume (R² = 0.73, rRMSE = 22%). We discuss the reliability of empirical approaches for estimates of canopy structural parameters considering the causality of light interaction and surface information. Abstract: Spatial information of forest structural variables is crucial for sustainable forest management planning, forest monitoring, and the assessment of forest ecosystem productivity. We investigate a complex alpine forest ecosystem located in the Swiss National Park (SNP) and apply empirical models to retrieve the structural variables canopy closure , basal area , and timber volume at plot scale. We used imaging spectrometer (IS) data from the Airborne Prism EXperiment (APEX) in combination with in-situ measurements of forest structural variables to develop empirical models. These models are based on simple and stepwise multiple regressions, while all potential two narrow-band combinations of the Simple Ratio (SR), the Normalized Difference Vegetation Index (NDVI), the perpendicular vegetation index (PVI), the second soil-adjusted vegetation index (SAVI2), and band depth indices were tested. The accuracy of the estimated structural attributes was evaluated using a leave-one-out cross-validation technique. Using stepwise multiple regression models, we obtained a moderate to good accuracy when estimating canopy closure (R 2 = 0.81, rRMSE = 10%), basal area (R 2 = 0.68, rRMSE = 20%), and timber volume (R 2 = 0.73, rRMSE = 22%). We discuss the reliability of empirical approaches for estimates of canopy structural parameters considering the causality of light interaction and surface information.


Introduction
Forest inventories seek to enumerate trees over a defined area and to gain forest structural attributes, such as volume, basal area, canopy cover, stem density, diameter at breast height (DBH), and maximum height at individual tree, plot, stand, regional, national, and global scale [1]. Structural attributes are crucial to determine information facilitating the management of forest ecosystems [2,3], to modify silvicultural practices [4,5], to quantify deforestation impacts on the carbon cycle and related impacts for global warming [6], and to assess productivity of forested areas in relation to forest composition [7] as well as in view of specific forest management practices [4,5]. Further, accurate estimates of structural variables are practical indicators for determining forest evolutionary history [8], monitoring forest sustainability [9], assessing insect infestation susceptibility [10,11], studying wildlife management and biodiversity [12], as well as forest water flux modeling [13]. Accordingly, there is an increasing need to generate accurate information regarding forest structural dynamics [14].
Traditionally, forest structure data have been collected by means of field surveys [15]. Field surveys are the most accurate way to collect forest structural data [16,17], but require in-situ measurements that are generally limited to a small area (plot area). Consequently, they do not band depth analysis [62], narrow-band normalized difference water indices (NDWI) [31], or from MS data by mostly either employing original spectral bands or computing broad-band vegetation indices (VI) [63][64][65]. As an example, Schlerf et al. [61] used HyMap data of highly managed and relatively homogenous Norway spruce (Picea abies) stands and reported significant linear relationships between forest leaf area index (LAI) and crown volume with a narrow-band perpendicular vegetation index (PVI). Cho et al. [39] tested the relationship between reflectance factors derived from HyMap data and forest structural attributes including mean DBH, mean tree height and tree density in a homogenous (i.e., homogeneous in DBH, height, tree density, and species) and closed canopy beech (Fagus sylvatica) stand. They conclude that their new vegetation indices, generated from contrasting spectral regions around the red-edge shoulder (756-820 nm) and the water absorption feature centered at 1200 nm (1172-1320 nm), showed higher correlation with structural variables compared with the standard vegetation indices derived from near-infrared and visible reflectance. Latifi et al. [66] have noted that the original bands of HyMap data located in the green and near-infrared parts of the spectrum contain useful information correlated to structural variables in a temperate forest site in Germany. All these studies suggest that the usage of narrow-band vegetation indices generated from the whole spectrum (350-2500 nm) offers a strategy to overcome the saturation problem at high canopy cover or LAI [67].
Relationships between vegetation indices and important stand characteristics are rarely documented [19]. Furthermore, since empirical relationships are at the cost of causality, their transferability is reduced and the usage of obtained results requires specific attention. We consequently aim to demonstrate the capability of IS data to empirically map structural attributes and to discuss the reliability of such approaches. We selected a heterogeneous alpine ecosystem comprising a natural forest including different tree species at varying proportions to empirically retrieve the forest structural attributes canopy closure, basal area, and volume from IS data. We tested various data analysis methods such as continuum removal (CR) and band depth analysis introduced by Kokaly and Clark [68] that were proposed to enhance spectral absorption features, to compensate disturbing background reflectance, and to enable improved linkage of spectral features to vegetation structural properties [69][70][71].

Study Area
The study was conducted in the Trupchun valley (Val Trupchun) within the Swiss National Park (SNP), located in the southeastern part of Switzerland ( Figure 1). Val Trupchun covers an area of 22 km 2 and the elevation ranges from 1775´3145 meters above sea level (a.s.l.) [72]. Based on climate data recorded by MeteoSwiss at the SNP's weather station, Buffalora at 1977 m a.sl., the climate of the SNP is dry and harsh with an average annual precipitation of 744˘160 mm (mean˘standard deviation (SD)) and an average annual temperature of 0.9˘0.5˝C [73]. The study site is a heterogeneous alpine landscape with a rough topography and steep slopes. The forest is classified as boreal type forests with European larch (Larix decidua L.) as a dominant tree. The Norway spruce (Picea abies (L.) Karst) and Swiss stone pine (Pinus cembra L.) are associated species [40]. The Norway spruce is present at a lower altitude, whereas the Swiss stone pine occurs at a high altitude. The trees are approximately 250 years old and have been strictly protected and unmanaged since 1914 [74].

Field Data
Field reference data of forest structural variables were collected in June 2012 in the Trupchun valley of the SNP. 35 Sample plots with a size of 30 m × 30 m were identified by applying a stratified random sampling, and excluding protected areas that were not accessible. Identified plots represent most of the present tree species and possible ranges of their properties, i.e., canopy cover, density, and species composition. A differential global positioning system (DGPS) providing an accuracy better than 1 m was used to record the geographic position of the four corners of the field plots. Within each plot, tree species, DBH, and tree height were measured, while only trees having a DBH (DBH measured at 1.3 m above the ground surface) greater than 5 cm were considered [75]. Basal area (BA), i.e., the cross-sectional areas of a tree stem measured at 1.3 m above ground [76], was estimated according to Equation (1): where BA is the estimated basal area (m 2 /ha), DBH is the diameter of a tree stem (m), and Plot is the area of the field plot (m 2 ). Height was measured for five randomly selected trees within each plot. A Haga hypsometer (Haga GmbH and Co. KG, Nuremberg, Germany) [39] was used to measure the height of the largest tree in terms of DBH and the trees located at the four corners of each plot. It should be noted that we had no a-priori information in term of DBH, height, and species type about the distribution of trees within each randomly selected plot, meaning that each tree has an equal probability of being selected from the respective population. To increase the selection chance of large trees (in terms of DBH), which are highly contributable to various forest structural properties, i.e., stand height, stand timber volume, stand biomass, and canopy cover, the height of a tree with the largest DBH in each sample plot was also measured. In total, the heights of 175 trees were measured in the study area. A typespecific height-DBH function was developed while considering tree species, DBH, and the measured height [77]. According to tree species, the corresponding type-specific height-DBH function was then applied to estimate the height of the unmeasured trees. The volume of the trees was calculated as a function of the measured DBH and estimated tree height. Tree volume was summed up to estimate the volume of the respective plot.

Field Data
Field reference data of forest structural variables were collected in June 2012 in the Trupchun valley of the SNP. 35 Sample plots with a size of 30 mˆ30 m were identified by applying a stratified random sampling, and excluding protected areas that were not accessible. Identified plots represent most of the present tree species and possible ranges of their properties, i.e., canopy cover, density, and species composition. A differential global positioning system (DGPS) providing an accuracy better than 1 m was used to record the geographic position of the four corners of the field plots. Within each plot, tree species, DBH, and tree height were measured, while only trees having a DBH (DBH measured at 1.3 m above the ground surface) greater than 5 cm were considered [75]. Basal area (BA), i.e., the cross-sectional areas of a tree stem measured at 1.3 m above ground [76], was estimated according to Equation (1): where BA is the estimated basal area (m 2 /ha), DBH is the diameter of a tree stem (m), and Plot is the area of the field plot (m 2 ). Height was measured for five randomly selected trees within each plot. A Haga hypsometer (Haga GmbH and Co. KG, Nuremberg, Germany) [39] was used to measure the height of the largest tree in terms of DBH and the trees located at the four corners of each plot. It should be noted that we had no a-priori information in term of DBH, height, and species type about the distribution of trees within each randomly selected plot, meaning that each tree has an equal probability of being selected from the respective population. To increase the selection chance of large trees (in terms of DBH), which are highly contributable to various forest structural properties, i.e., stand height, stand timber volume, stand biomass, and canopy cover, the height of a tree with the largest DBH in each sample plot was also measured. In total, the heights of 175 trees were measured in the study area. A type-specific height-DBH function was developed while considering tree species, DBH, and the measured height [77]. According to tree species, the corresponding type-specific height-DBH function was then applied to estimate the height of the unmeasured trees. The volume of the trees was calculated as a function of the measured DBH and estimated tree height. Tree volume was summed up to estimate the volume of the respective plot.
Remote Sens. 2015, 7, Canopy closure (CC) is defined as the percentage of ground within 30 mˆ30 m, covered by the vertical projection of the overlying trees' crown [21]. The CC was visually estimated using a digital aerial photograph acquired in 2006 with 25 cm spatial resolution.
It is important to mention that there was a two-year time lag between field reference data collection and imaging spectrometer data acquisition and also a four-year time lag between aerial photograph and imaging spectrometer data acquisition. However, the newest reference available providing a sub-meter resolution were aerial photographs from 2006. Since the average tree ring width increase is 3 mm/year [78], a stem growth of less than 1 cm/year can be expected and canopy cover change from 2006-2010 is assumed to be negligible. Further, the SNP is a protected area without any human activity and no occurrence of significant natural disturbance (i.e., fire, windstorm, or bark beetle attack) has been reported in Val Trupchun after 2006 or even before that time, thus, we assumed that the forest structure attributes remained relatively constant during this period [13,49].

APEX Imaging Spectrometer Data
Airborne Prism EXperiment (APEX) imaging spectrometer data were used in this study. APEX is an airborne pushbroom device with 1000 imaging pixels across track and covers a swath width of 1.5-5 km (depending on flight altitude) with a field of view (FOV) of 28˝ [43]. The spectral wavelength range of 372-2500 nm is recorded by a Visible/Near Infrared (VNIR), and a Shortwave Infrared (SWIR) detector. The spectral resolution (FWHM) ranges between 0.7 and 9.7 nm for the VNIR and between 6.2 and 12 nm for the SWIR detector. After removal of noisy bands, an APEX data cube contains 285 spectral bands.
APEX data over the SNP were acquired on 26 June 2010 under cloud free conditions and with a 2 m pixel size. The acquisition time was close to solar noon. The sun zenith and azimuth angle ranged between 31.8˝and 28.1˝and between 127.5˝and 139.1˝. A parametric geocoding approach (PARGE) using 15 ground control points with the Swiss national coordinate system (CH 1903 LV03) was applied to the data [79], resulting in a geometric accuracy of 3.2 m˘1.4 m. This corresponds to a total pixel shift of 1.3-1.9 pixels˘0.5-0.8 pixel. Measured APEX radiance data were atmospherically corrected using an atmospheric and topographic correction approach for rugged terrain (ATCOR-4) [80]. ATCOR-4 comprises the atmospheric radiative transfer code MODATRAN-5 [81] to pre-calculated look-up tables (LUT) of atmospheric variables. We set the atmosphere type and the aerosol model to mid-latitude summer and a rural aerosol model. These settings in addition to image based estimates of atmospheric water vapor and visibility provided by the pixel-wise selection of LUT entries allowed to compensate for atmospheric effects and to determine top-of-canopy hemispherical conical reflectance factor (HCRF) data (for terminology see [82]).
Previous studies have shown that the narrow band indices can provide significant improvement as compared to broadband indices for estimating vegetation structural variables [39,61,67]. Consequently, we only focused on narrow band indices.

Band Depth Indices
Band depth indices generated from continuum-removed spectra were used to estimate structural properties. This approach has shown its capability to improve the performance of structural variable estimation [62,71]. Based on previous research, the continuum removal was applied to four regions of the spectrum, ranging from the VIS to the SWIR part ( Table 1). The selected region in VIS is affected by chlorophyll absorption, whereas the NIR [62] and SWIR regions are mainly affected by water, lignin, and cellulose absorption [68]. The inversion of the reflectance factors was applied on SWIR1 and SWIR2 using Equation (6). This procedure was necessary to reshape the concave form of the spectrum to a convex form in order to calculate banddepth indices [71].
where ρ inverse is the inverse reflectance factor at band i and ρ i is the original reflectance factor at band i. Table 1. Start, end, and center points of selected spectral regions for continuum removal and subsequent calculation of band depth indices (VIS = visible domain, NIR = near-infrared, SWIR = shortwave infrared).

Spectral Domain Location of the Selected Wavelength Start of Spectral Region
Analysis (nm)

End of Spectral Region
Analysis ( Band depth indices including band depth (BD), normalized band depth ratio (BDR), normalized band depth index (NBDI), the area under curve (AUC) of a continuum removed spectrum, as well as the area under continuum-removed curve normalized to the chlorophyll absorption band depth (ANCB) [90], were implemented. They were calculated according to the following Equations: NBDI " BD´BD c BD`BD c (10) ANCB " AUC BD c (12) where R 1 is the continuum removed spectrum, R is the original reflectance value, R c is the reflectance value of the continuum line [68], BDC is the band depth at the band center, λ i`1 and λ i are wavelengths of the i and i + 1 bands, BD i + 1 and BD i are band depths at the i and i + 1 bands, n is the number of used spectral bands.

Correlograms
In order to identify sensitive wavelength regions to predict forest structural attributes a 2D correlation plot (correlogram) was calculated for all possible band combinations of SR, NDVI, SAVI2, and PVI type. These correlograms represent the coefficients of determination (R 2 ) of the individual linear regression models relating APEX derived VIs and in-situ measurements of each target variable (i.e., canopy closure, basal area, and forest volume) [39]. The spectral bands with strong water absorption features at 1335-1490 nm and at 1780-1990 nm were excluded from further analysis. Therefore, 193 APEX spectral bands were used to generate correlograms, allowing to compute 37,249 different narrow-band combinations. The optimal band combination was selected based on R 2 (i.e., we used an R 2 threshold of 0.5 with a confidence level of 95%), considering multi-collinearity issues, causality explanations, and focusing on "hot-spot areas" (spectral range) [40]. The selected band combinations were then used in regression analysis.

Regression Analysis
Regression analysis is an empirical approach to determine the relationships among variables, i.e., dependent and independent variables (or predictors). In remote sensing literature regression analysis has seen a widespread use due to its ease of implementation and efficiency [91]. We used simple and stepwise multiple regression to predict forest structural attributes from APEX IS metrics. An implemented standard deviation outlier labeling method, i.e., the 3SD method [92], identified one field plot as an outlier that was removed from further analysis. Normality distribution of dependent variables is an important assumption associated with regression [39]. The Shapiro-Wilk test indicated that the BA (W = 0.975, p = 0.579), volume (W = 0.974, p = 0.566), and CC (W = 0.965, p = 0.328) followed a normal distribution (p > 0.05). A 15ˆ15 pixels window (30 mˆ30 m) was used to extract average reflectance values corresponding to each field plot from APEX IS data.
Stepwise multiple regression was run in two phases. First, we performed a stepwise multiple regression with optimal narrow-band indices to predict their performance to estimate structural variables. Afterwards, we ran the stepwise multiple regression with optimal narrow-band and band depth indices together to assess the influence of the band depth indices on the performance of the model [65]. The SPSS 20.0 package (IBM Crop. 2011) was used for statistical analysis. The probability of F was employed as a criterion to include (F < 0.05) or to remove (F > 0.10), an independent variable in the regression analysis [70]. The best model was selected by comparing R 2 , RMSE, and variance inflation factor (VIF). The acceptable VIF value was set to smaller than 10 to avoid multi-collinearity among variables and potential over-fitting of the model [93].

Validation
The results were tested by comparing the estimated forest structural attributes with the actual measurements obtained from field inventory. Because of the limited number of field plots, a eave-one-out cross-validation technique (LOOCV) was used to evaluate the accuracy of the estimated structural attribute [49]. The root mean square error (RMSE) and mean relative RMSE (rRMSE (%)) according to [1,66] were employed to quantify the reliability of the predictions.
rRMSE " RMSE yˆ1 00 (14) whereŷ i is the estimated basal area or canopy closure or volume, y i is the measured basal area or canopy closure or volume, n is the number of observations and y is the mean of the measured variable (basal area or canopy closure or volume). The LOOCV has the advantage of providing an unbiased estimation of the prediction error [23]. In addition, we evaluated the performance of selected models using a bootstrapping approach with 5000 replications. The bootstrapping approach enables to increase the robustness of a model where the number of field plots is limited [94,95]. Three metrics, i.e., the mean of R 2˘o ne SD, the mean of RMSE˘one SD, and the mean rRMSE˘one SD were calculated to assess the model fit and predictive performance. The uncertainty was quantified at the 95% confidence interval.

In-Situ Forest Structural Measurements
Descriptive statistics of the in-situ data are shown in Table 2. The measured volume ranged between 77.6 m 3 /ha and 807 m 3 /ha. The average volume was 421.7˘181.7 m 3 /ha, roughly representing the full range of volume values reported elsewhere for alpine areas [96,97]. The mean basal area was 42.5˘16 m 2 / ha and the canopy closure ranged between 22% and 85% with an average of 61%˘15.3%. Our findings are in line with results reported in other studies performed in the SNP [74,98]. Table 2.
Forest structural characteristics derived from reference data in the study area (DBH = diameter at breast height, SD = standard deviation, CV = coefficient of variation). We calculated correlograms representing the coefficients of determination (R 2 ) of individual linear regression models relating APEX derived VIs and in-situ measured CC for all possible two-band combinations of SR, NDVI, PVI, and SAVI2 type. Based on these correlograms, we identified potential band combinations and wavelength regions showing high sensitivity to variations of forest canopy closure ( Figure S1).

Parameter
In general, similar patterns in terms of sensitive wavelengths combinations are observed for SR and NDVI ( Figure S1a,b)). Sensitive combinations are mainly located in the shortwave infrared (SWIR). PVI type VIs shows a different pattern and higher R 2 compared to SR, NDVI, and SAVI2. Potential band combinations for PVI type VIs are located in (i) the VIS part in combination with NIR and SWIR, (ii) the red-edge region in combination with red wavelengths, and (iii) the SWIR part in combination with the NIR or the SWIR. Table S1 provides the results of the selection of candidate VIs derived from 2D correlograms based on a set of criteria introduced in the method section. The generated VIs obtained different results. The maximum R 2 values were 0.61, 0.62, 0.65, and 0.76 using NDVI, SR, SAVI2, and PVI, respectively. This indicates that a PVI type index generated from red-edge and NIR regions (883 and 763 nm) provided more information to estimate canopy cover than other categories of VIs (SR, NDVI, and SAVI2 type). The NIR and red-edge regions (895 and 754 nm) again were found suitable to estimate canopy closure using a SAVI2 type index (R 2 = 0.65). For SR and NDVI type VIs, spectral wavelengths with highest R 2 are located at 2204 and 2253 nm.
In a second step of the regression analysis, a stepwise multiple regression was applied considering all narrowband SR, NDVI, SAVI2, PVI type VIs that performed best. These are listed in Table S1 in addition to selected band depth indices (i.e., BD, BDR, NBDI, AUC, and ANCB) as independent variables to test whether band depth indices improve the prediction performance. As shown in Table 3, using the stepwise multiple regression increased the R 2 to 0.90. The best performing stepwise multiple regression model has fulfilled the statistical collinearity requirement with the VIF values < 10 [93]. This model uses band depth indices generated from the NIR and SWIR2 regions as well as an SR narrow-band VI employing wavelengths from the SWIR (2204 nm, 2253 nm) region. The corresponding wavelengths for each regression variable are summarized in Table 3. The two models, i.e., the simple and stepwise multiple regression models, were tested for significance using analysis of variance (ANOVA) at 5% significance level. The obtained p value (p = 0.0005) indicates that there is a statistically significant difference between the two models. In conclusion, the combination of band depth indices and an SR index increases the reliability of the prediction model for canopy closure compared to a model based on a single PVI type index. Consequently, stepwise multiple regression improved the model performance. The final stepwise multiple regression model presented in Table 3 was applied to the entire APEX data set to generate the canopy closure map (Figure 2). The estimated values range from 0% to 90% with an average of 56%˘15%.
The relationship between measured and estimated canopy closure obtained by applying a LOOCV approach is depicted in Figure 3. The model employing a single variable (i.e., a PVI type index) gained reasonable results (R 2 = 0.60, RMSE = 9.34, and rRMSE = 15%). Predicted values showed a good agreement with high canopy closure values (CC > 60%). However, this relationship is weak for low and medium canopy closure values of 20%-60% (Figure 3a). The results indicate that PVI generated from red-edge wavelengths shows a high sensitivity to denser canopy closure, and thus low sensitivity to open canopies.   Table 3. Coefficient of determination (R 2 ± SD), root mean square error (RMSE ± SD), and relative root mean square error of prediction (rRMSE ± SD) are shown.   Table 3. Coefficient of determination (R 2 ± SD), root mean square error (RMSE ± SD), and relative root mean square error of prediction (rRMSE ± SD) are shown.  Table 3. Coefficient of determination (R 2˘S D), root mean square error (RMSE˘SD), and relative root mean square error of prediction (rRMSE˘SD) are shown.
The stepwise multiple regression model yielded higher accuracy with an R 2 of 0.81, and a rRMSE of 10% (Figure 3b). This model significantly increased the R 2 and reduced the relative RMSE error. Table S2 and Figures S2 and S3 illustrate the obtained mean values and distribution of R 2 , RMSE, and rRMSE generated from the bootstrapping approach, yielding slightly better results compared to the cross validation approach (i.e., R 2 of 0.82 and an rRMSE of 10%). As it can be seen in Figure 3b, the combination of band depth indices and an SR type index substantially improved prediction power, especially for low and medium canopy closure (20%-60%). This indicates that the selected VI and band depth indices are highly sensitive to canopy closure and can minimize disturbing effects such as soil background, atmospheric, and sun-target-sensor geometry effects [99].

Basal Area
The maximum R 2 obtained from the correlation plots for basal area considering SR, NDVI, and PVI is almost identical (R 2 > 0.60). The pattern of variation in R 2 is similar for SR and NDNI type VIs ( Figure S4a,b)). In addition, these figures show that the narrow-band combinations generated mainly from narrow portions of the SWIR can explain variability in basal area most effectively. Contrary to SR and NDVI type VIs, the PVI and SAVI2 type VIs ( Figure S4c,d) resulted in different correlograms in terms of sensitive wavelengths. Figure S4 shows that a large portion of the spectrum could potentially be used to predict basal area. Sensitive combinations are mainly generated from (i) the VIS part in combination with NIR and SWIR parts, (ii) the red-edge region in combination with NIR wavelengths, and (iii) a combination of NIR wavelengths themselves or together with the SWIR part. Interestingly, the highest variability explained by SAVI2 type indices (R 2 = 0.44) is significantly smaller than the highest R 2 for SR, NDVI, and PVI (0.69, 0.69, 0.64).
The optimal narrowband VIs derived from 2D correlograms to predict basal area are presented in Table S3. An SR type index that can be calculated from the SWIR (2385 nm, 1993 nm), and a NDVI type index located in the SWIR (1993 nm, 2391 nm) contained most information to estimate the basal area. The stepwise multiple regression using all optimal VIs did not improve model performance (data not shown). Several studies have reported a strong correlation between the SWIR domain and basal area [10,11,65,100]. The SWIR region is highly affected by canopy water content [50], canopy biochemical properties like cellulose, lignin, and protein [101], and shadowing [102]. Forest structural attributes such as basal area are related to the addressed variables. Our results show that the selected VIs from the SWIR region are best correlated to basal area (Table S3).
Stepwise multiple regression was subsequently performed to assess the joint usage of optimal narrow-band VIs and band depth indices to predict basal area. As is evident from Table 4, the best performing stepwise multiple regression model (including an SR type VI and a normalized band depth index (NBDI) obtained a higher R 2 compared to the best performing simple linear model. The result of an ANOVA test (p = 0.017) indicated a statistically significant difference between the two models at 5% significance level. The significant wavelengths in both the simple and the stepwise multiple regression modesl are located in the SWIR part of the spectrum. For forest basal area estimation NDVI, PVI, and SAVI2 type indices did not lead to improved model fits. The basal area thematic map was generated by employing the stepwise multiple regression model and the APEX IS data set (Figure 4). The estimated values ranged from 0 to 75 m 2 /ha with an average of 42˘17 m 2 /ha.   Table 4. Coefficient of determination (R 2 ± SD), root mean square error (RMSE ± SD), and relative root mean square error of prediction (rRMSE ± SD) are shown.    Table 4. Coefficient of determination (R 2 ± SD), root mean square error (RMSE ± SD), and relative root mean square error of prediction (rRMSE ± SD) are shown.  Table 4. Coefficient of determination (R 2˘S D), root mean square error (RMSE˘SD), and relative root mean square error of prediction (rRMSE˘SD) are shown.
The relationship between measured and estimated basal area obtained by applying a LOOCV approach is illustrated in Figure 5. The simple linear model yielded an R 2 value of 0.60, a RMSE value of 9 m 2 /ha, and an rRMSE value of 23% (Figure 5a). The stepwise multiple regression model produced a higher R 2 (R 2 = 0.68), and decreased the RMSE and the rRMSE values to 8 m 2 /ha and 20%. This result is in line with the findings of Mutanga & Skidmore [62], who concluded that the use of stepwise multiple regression using band depth information improved structural parameter estimation (biomass). The suitability of band depth indices to assess vegetation structural attributes has been indicated by previous studies [52,69,90,103,104].
As can be observed from Figure 5a,b, the linear relationship between estimated and measured basal area is stronger for the stepwise multiple regression model than for the simple model. In fact, for the simple model, the trend line is closer to power form (non-linear) than to a linear relationship in a way that applying an exponential model would increase the model performance (R 2 = 0.63). The bootstrapping approach (Table S4 and Figures S5 and S6) confirms the statistics of the leave-one-out cross validation approach and yielded a R 2 of 0.68 and a rRMSE of 20%.

Volume
Considering all potential band combinations for SR, NDVI, PVI, and SAVI2 type indices for forest volume, the strength of relationship increased from the VIS to the SWIR for all generated VIs ( Figure S7). Maximum R 2 was higher for SR, NDVI, and PVI type indices compared to SAVI2. The similar sensitive spectral region to predict volume using SR and NDVI type indices observed and mainly located in the SWIR (Figure S7a,b). These two types of VIs are functionally similar and thus contain the same information [105]. In addition to the SWIR region, parts of the VIS and NIR regions showed sensitivity to predict volume using PVI and SAVI2 type indices ( Figure S7c,d). Table S5 lists the best narrow-band VIs derived from 2D correlograms based on a set of criteria introduced in the method section to estimate volume. The maximum R 2 value varied from 0.61, 0.61, 0.64-0.38 using the SR, NDVI, PVI, and SAVI2. A PVI type index located in the SWIR (1993 nm and 2091 nm) provided more information to estimate forest volume than the other categories of VIs (i.e., SR, NDVI, and SAVI2 type). The predictive performance of an SR and an NDVI type index (1993 nm and 2385 nm), both with R 2 = 0.61, was slightly lower than in the case of PVI. Table S5 indicates the suitability of the SWIR region for forest volume estimation. Our findings correspond with Schlerf et al. [61], who reported highest correlations of narrow-band indices with forest crown volume in wavelengths related to water absorption features in the SWIR. In contrast, Cho et al. [39] found that VIs based on the contrast between reflectance values in the red-edge shoulder of the spectrum and in water absorption features are related to forest structural attributes. The relationship between the SWIR spectral information and forest volume has been reported elsewhere [7,85]. However, some studiesreported red and near-infrared spectral bands as good predictors to estimate forest volume [17]. A stepwise multiple regression analysis on all optimal VIs listed in Table S5 did not improve the model performance (data not shown).
The combined use of optimal narrow-band VIs and band depth indices was assessed to predict forest volume using a stepwise multiple regression. It is evident from Table 5 that the stepwise multiple regression model (including a PVI type VI and an NBDI index) obtained a higher R 2 compared to the best-performing simple linear model. The result of an ANOVA test (p = 0.024) indicated a statistically significant difference between the two models at 5% significance level. The significant wavelengths in both the simple and the stepwise multiple regression model are located in the SWIR part of the spectrum. From our results, it appears that the band depth indices have a high potential for an improved estimation of forest volume and forest structure attributes in general. This might be due to the capability of such indices to minimize effects of external influences, such as atmospheric water vapor, soil background [71], and BRDF effects [68].
The volume thematic map generated from APEX IS data is shown in Figure 6. The average estimated forest volume was 431˘149 m 3 /ha. Figure 7 illustrates the agreement between measured and estimated forest volume for the simple and the stepwise multiple regression model. The simple model explained 68% of variability in volume estimation (R 2 = 0.68, RMSE = 101 m 3 /ha, and rRMSE = 24) based on a PVI type index involving wavelengths at 1993 nm and 2091 nm.     Table 5. Coefficient of determination (R 2 ± SD), root mean square error (RMSE ± SD), and relative root mean square error of prediction (RMSE% ± SD) are shown.     Table 5. Coefficient of determination (R 2 ± SD), root mean square error (RMSE ± SD), and relative root mean square error of prediction (RMSE% ± SD) are shown.  Table 5. Coefficient of determination (R 2˘S D), root mean square error (RMSE˘SD), and relative root mean square error of prediction (RMSE%˘SD) are shown.
The stepwise multiple regression model showed an improvement in volume prediction (R 2 = 0.73, RMSE = 90 m 3 /ha, and rRMSE = 22%) including a combination of a PVI type index and a band depth index (NBDI). The bootstrapping (Table S6 and Figures S8 and S9) yielded similar results compared to cross validation approach, with a R 2 of 0.72 and a rRMSE of 22%.

Canopy Closure
Sensitive band combinations to predict canopy closure were selected from the red-edge, NIR, and SWIR part of the spectrum. These results are in line with previous studies suggesting optimal combinations for estimation of forest structural attributes in the red-edge region [62,67,106] and in the NIR and SWIR spectral regions [55,58,107]. Our canopy closure estimates are comparable or even higher than results reported elsewhere. For example, Schlerf et al. [61] employed a PVI index based on wavelengths at 885 nm and 954 nm with an R 2 of 0.79 to map crown volume in a relatively homogenous Norway spruce (Picea abies) stand. Although we assessed the capability of APEX IS data in a heterogeneous ecosystem with complicated stand structure and tree composition, our results are slightly higher than those reported in Schlerf et al. [61]. We conclude that heterogeneity helps to increase the spectral sensitivity required to infer forest structural variables due to increased contrast between canopy and background reflectance [108]. In other words, a highly dynamic spectral range, as found in heterogeneous areas, contributes to improved forest canopy closure estimations [109]. Reported accuracies for canopy closure estimates in previous studies utilizing multispectral [7,[110][111][112], hyperspectral [22], or a fusion of active and passive data [21], ranged between 0.50-0.88 (R 2 ) and 15%-30% (relative RMSE). The performances of our simple and stepwise multiple regression model lie in the reported accuracy ranges. However, the relative RMSE error of 10% indicates a remarkably higher accuracy of canopy closure estimation. It is worth noting that an acceptable level of error for operational forest inventory is 10% [113].

Basal Area
The accuracy of the stepwise multiple regression to retrieve basal area is better than previously reported results that employed IS data to estimate basal area. Hyyppä et al. [15], for example, achieved an R 2 = 0.58 and a relative RMSE of 35% using data from an AISA sensor. Defibaughy Chávez & Tullis [20] estimated the basal area with an R 2 less than 0.30 and a relative RMSE of 50%. Anderson et al. [114] reported an R 2 value of 0.40 for the basal area prediction using Airborne Visible/Infrared Imaging Sensor (AVIRIS) data. Lefsky et al. [25] predicted basal area with an R 2 value of 0.36 using AVIRIS data. However, our R 2 value was lower than the R 2 value of Welter et al. [110]. They used SPOT5 data to model the basal area. Literature reports R 2 and RMSE% for estimating basal area between 0.20-0.75 and 20%-50%, respectively [20,110,114]. The performance of our simple and stepwise multiple regression models lies in the accuracy range reported in previous studies.

Volume
Our two models tend to overestimate low values and underestimate high volume values (Figure 7), an effect known as "local bias" [39]. Deviations between estimated and measured forest volume are largest for values around 300-500 m 3 /ha, strongly affecting model accuracies. A possible explanation of this effect may be related to changes in composition and vertical structure of the respective forest stands. In general, they consist of a dominant tree (larch), which forms the upper tree story, and two associated species (i.e., Swiss stone pine and Norway spruce), which dominate the middle and lower parts of the canopy. As the upper part of the forest story is clearly seen by an imaging device, the underlying tree stories may be hidden. Therefore, the true spectral variation cannot be accurately modeled by a simple relationship [57], resulting in a poor relationship between measured and estimated forest volume. This problem might be partly solved, if a type specific calibration is applied or additional forest attributes (e.g., stand age, site index) are considered [17]. Nevertheless, the performance of our models is more accurate than what was achieved by previous studies that employed IS data to estimate forest volume. Hyyppä et al. [15] predicted forest volume with an R 2 of 0.55 and a relative RMSE of 45% using AISA IS data with 30 spectral bands (466-870 nm). In the case of APEX, the availability of a larger spectral range incorporating the SWIR region might explain why our models yielded a higher accuracy. Recently, Hill et al. [96] assessed the capability of LIDAR Canopy Height Models (CHM) and forest inventory data to map timber volume in a mountainous study site in Eastern Switzerland using a multiple linear regression model approach. They obtained an R 2 value of 0.64 and a RMSE of 123 m 3 /ha (31% of the mean). Interestingly, both our simple and stepwise multiple regression model performed slightly better [96].
The relative error of volume estimation using optical remote sensing has been reported between 20%-75% [63,115,116]. It should be mentioned that for operational forest management the acceptable level of estimation error is 15% at stand level [115]. Although APEX derived forest volume showed a significant improvement in estimation performance, the method is not yet accurate enough for operational forest management from an accuracy requirements point of view.

Reliability of Structural Parameter Retrieval from APEX Data
Our results are in line with many other studies [39,46,47,61] and clearly demonstrate that empirical-statistical methods allow relating spectroscopy measurements with forest structural variables and, thus, quantifying them. However, it is important to interpret such empirical results in the context of causality. Measurements of reflected radiance in the optical spectral domain are per se less or not sensitive for forest structural variables, particularly for the variables basal area or timber volume. In fact, optical measurements are mostly representative for the upper canopy part, and derived information tends, for example, to saturate with increasing canopy volume [40,117,118]. This indicates that results obtained in our study are likely based on secondary relationships. The explanation why we found significant correlations between IS data and forest structural variables is that certain forest structural variables (i.e., basal area, volume) are related to structural variables (i.e., canopy cover) directly affecting the radiative transfer of optical measurements. In fact, based on our field data we observed significant correlation between canopy cover and both volume (r(35) = 0.62, p < 0.01) and basal area (r(35) = 0.74, p < 0.01) (Table S7). In addition, the combination of forest structural variables and sun illumination causes specific patterns of cast shadow. These illumination effects often persist after atmospheric compensation, particularly for heterogeneous surfaces measured in high spatial resolution [119]. These processing artifacts eventually translate into optical indices derived and allow relating them with canopy structural variables as shown with our results.
The contradiction between good empirical findings obtained and missing causality implies two things. First, it is possible to apply relatively simple methods to retrieve forest structural information from spectroscopy data, if extensive reference data are available and the empirical models are only applied to a specific canopy at particular time, observed with a specific sensor [41]. Second, obtained results must be interpreted very carefully: extreme values or observed anomalies are likely caused by the limited representativeness of underlying empirical models.

Strategies for Improved Simultaneous Retrieval of Biochemical and Structural Plant Traits
The great advantage of optical RS data is that they allow retrieving various plant functional traits and these data are frequently used to estimate forest structural attributes. According to our discussion above, several problems are, however, related to the use of optical RS data in combination with empirical approaches to retrieve structural variables. Particularly in closed dense forests with a complex and multi-layered canopy [25,120], spectral signals and eventually derived empirical models saturate. Further, calculated surface reflectances are typically affected by shadowing and BRDF effects, contributions of the soil background, and atmospheric conditions. Active RS systems such as LiDAR are more suitable for capturing structural forest information, but also suffer from certain limitations (i.e., data are expensive, their availability is limited in coverage, they may be prone to low point density, and algorithms to generate digital surface and terrain model suffer from shortcomings) and, importantly, do not allow obtaining functional or biochemical plant traits.
With increasing evidence it can be concluded that there is no single RS technology available to consistently retrieve forest structural attributes [20,21,121,122]. The combination of both technologies is considered the most promising strategy to provide highly accurate estimates of forest structural attributes [66,114,123] and other plant traits [124,125]. Future research attempts are needed to properly evaluate the capability of combined approaches and to tackle various technical issues such as the temporal and geometric integration of these complementary observational approaches.

Conclusions
We conclude that APEX data in combination with empirical-statistical approaches allow estimating forest structural attributes including canopy closure, basal area, and volume in heterogeneous alpine ecosystems. The two tested simple and stepwise multiple regression models, applied to narrow-band indices and band depth indices, provide a relatively simple analytical framework for forest structural parameter retrievals but at the cost of causality between surface information and light interaction. Among the structural attributes, canopy closure was predicted more accurately than basal area and volume. The SWIR region of the reflectance spectrum was found to be particularly sensitive to forest structural attributes and could be used to predict relevant parameters. The availability of a high quality imaging spectrometer with high signal to noise ratio (SNR), especially in the SWIR region, is, thus, crucial for forest research. Once a simple model was applied, the use of a PVI type narrow-band index showed a strong correlation with canopy closure and volume. However, more variation was explained by using a SR type index to estimate basal area. Stepwise multiple regression models based on a combined use of narrow-band VIs and band depth indices provided improved estimates of forest structural attributes compared to simple models that apply only narrow-band type indices. We therefore conclude that band depth indices generated from continuum removed spectra in combination with simple narrow-band indices can be used to predict forest structural variables, in particular for canopy closure assessment. This methodology, however, needs further investigation in different forest ecosystems with different structural and species compositions.
Considerable effort has been made by numerous researchers to successfully estimate forest structural attributes using empirical approaches [91]. However, results obtained from empirical approaches are limited in their representativeness and are only valid for specific vegetation types, phenological state, and sensors. Further, underlying empirical models rely on extensive in-situ data. Both aspects question the operational applicability of empirical approaches to retrieve complex surface information such as forest structural variables. Physically-based approaches may partly overcome these limitations but will still face the problem of a limited sensitivity of optical data for structural vegetation variables.
We suggest combining active and passive RS, two complementary techniques that provide optimal capability and sensitivity to characterize complex forest ecosystems. Active systems are the preferred solution to obtain various structural information, while optical systems can provide complementary information about other important plant functional traits.