Assessing the Robustness of Multispectral Satellite Imagery with LiDAR Topographic Attributes and Ancillary Data to Predict Vertical Structure in a Wet Eucalypt Forest

1. Introduction

Vertical forest structure (VFS) is an essential component of ecological processes [1], habitat quality and biodiversity [2,3]. Knowledge about VFS is helpful in relation to timber production and carbon sequestration [4,5,6], the habitat requirements of different organisms [7], and the distribution of fuels and fire behaviour [8,9,10]. Detailed knowledge of VFS is required to manage forests for estimating forest water storage capacity [11], wildlife habitat management, wildfire mitigation [4,10,12], and understanding habitat restoration plans [13]. However, measuring VFS directly in the field is an ongoing challenge for forest managers [14] and is expensive, time-consuming, and even logistically difficult [15,16]. Thus, ecologically meaningful and robust measures of VFS for remote, inaccessible, and large geographical areas have been lacking [3,17] and require further research, in particular for sub-canopy metrics [4].

Remote Sensing satellite data have been used for modelling and mapping forest structural attributes in large geographical areas [16,18,19]). Mapping forest structure attributes using remotely sensed data has been a topic of research studies for over four decades [18,20]. Although aerial photographs were previously used for mapping forest stands, those techniques were highly subjective, time-consuming, and depended on the experience of the interpreter [20,21]. With the advancements in technologies, automated image analysis techniques were explored to retrieve forest stand structural attributes [20], and satellite imagery was found to meet the needs of spatially contiguous information about forests over large geographical areas [22,23]. Although traditional passive remote sensing techniques have demonstrated the potential to provide useful information on forest structure attributes based on spectral reflectance and derived vegetation indices [24], these techniques are not currently sufficient to predict and map three-dimensional (3D) VFS for forestry and ecological applications [25]. Also, the successful application of remote sensing data depends on identifying the optimum spatial resolution [26] for monitoring multiscale biodiversity attributes [27]. To date, the assessment of optimum spatial resolution remains an area of active research, especially for heterogeneous ecosystems [28].

As an alternative to passive remote sensing techniques, Light Detection and Ranging (LiDAR), an active remote sensing technology, can directly measure the 3D vegetation structure characteristics by penetrating through canopy gaps [4,29,30,31], and provide valuable information for the ecological and forestry applications [16]. This technology has been widely used in ecological studies [32] and provides a valuable tool in forest inventory [8,9,33]. In addition, LiDAR has been used for the assessment of structure response to a range of bushfire events [12,34], wildlife-habitat relationships [35], and bird species richness [36,37]. LiDAR is the most popular approach to survey VFS and can bridge the gap between satellite and field data [38,39]. The increasing capabilities of LiDAR technology allow researchers to characterise forest structure [16,40] and could also be used to quantify structural diversity at multiple scales [41].

LiDAR can capture VFS even in high biomass forests, where optical multispectral satellite data tend to be saturated [18,42,43] and cannot provide information about the vegetation strata underneath the canopy cover [44]. Although LiDAR point clouds can capture VFS, their application to large geographical areas is challenging due to logistical, cost, and data volume constraints [38]. High-resolution LiDAR is likely to be cost-prohibitive for broadscale applications. So, there is considerable interest in integrating operational resolution LiDAR and multispectral satellite data by natural resource managers and the forest industry to retrieve forest structural attributes [8,14,45]. Considering the different aspects of VFS, no single sensor can exhibit all the required information relevant to forest managers, and the integration of multispectral satellite imagery and LiDAR data could advance the prediction and mapping of forest structure characteristics [42].

Topography plays a significant role in the development of VFS, affecting abiotic conditions, i.e., microclimate and edaphic factors [46,47] that can influence the structure, function, and dynamics of ecological communities [46]. Topographic attributes provide a basis for understanding ecological relationships with vegetation, landforms, and soils [48].

There is relatively limited research into predicting VFS from LiDAR-derived topographic attributes with other sources of remote sensing satellite data, especially in mixed-species forests [8,49]. Due to the complex forest structure as well as topographic variations, stratifying a complex forest on a large geographical area is still challenging [50]. This study aims to evaluate the robustness of predicting vertical forest strata in a wet eucalypt forest (Tasmania, Australia) based on multispectral satellite imagery combined with topographic attributes derived from a digital terrain model (DTM). Moreover, this study assesses whether selected vegetation indices and image texture features improve the prediction accuracy of the VFS attributes. To provide insights into the different sources of datasets owing to different scales, this study examines the importance of predictor variables and investigates how three different pixel sizes influence model outcomes. Thus, this study addresses data complexities, including multidimensionality and nonlinearity in multisource data, and provides a robust approach for assessing forest structure.

2. Materials and Methods

2.1. Study Site

We conducted this research within a 5 km by 5 km area centred at 43.104° S and 146.656° E, located at the Warra Supersite within the Tasmanian World Heritage Area, Australia. The selected site (Figure 1) is topographically and geologically complex, comprising harvested sites, rivers, and gullies. The ground elevation ranges from 51 m to 648 m above mean sea level according to the 30 m resolution DTM. The Warra site is a member of the Terrestrial Ecosystem Research Network (TERN) Australian Supersite Network [51], which covers the cool, temperate wet forest biome. Mean annual rainfall is 1,707 mm and mean daily temperature ranges are 8.3°C to 19.3°C in January (summer) and 2.5°C to 8.6°C in July (winter) [52].

The Warra Supersite is one of the most productive terrestrial ecosystems in the world, and its management generates a high level of political and social interest. This Supersite has a long history of fire, with several fires since 1850, and ranges from intensively managed to protected forests [53]. There is a general pattern of ecological succession from wet sclerophyll-dominated younger stands to older mixed forests consisting of eucalypts with a rainforest understorey [54]. The forest at the study site is mostly dominated by Eucalyptus obliqua, which forms a tall (~50 m) overstorey overtopping secondary stratum of trees, for example, Acacia dealbata and Acacia melanoxylon, or rainforest trees such as Nothofagus cunninghamii, Atherosperma moschatum and Eucryphia lucida. Understorey tree and shrub species vary from dense Melaleuca squarrosa and Leptospermum lanigerum on soils with impeded drainage to Nematolepis squamea and Pomaderris apetala on well-drained soils [55]. The understorey includes tree ferns (Dicksonia antarctica), Bauera rubioides, the tall sedge Gahnia grandis, and other understorey shrubs and trees of up to ~10 m in height.

2.2. Remote Sensing Data

We used airborne LiDAR point clouds in combination with WorldView-3 and Landsat-8 (OLI) satellite multispectral remote sensing datasets and ancillary data, i.e., a geology map layer [56]. Airborne LiDAR data used in this study was collected by Airborne Research Australia using a Riegl LMS-Q560 laser sensor on Diamond Aircraft HK36TTC ECO-Dimonas. LiDAR data with a spatial resolution of 28.66 points m^-2 (horizonal spacing = 0.19 m) and a 1064 nm laser wavelength was collected from a flying height of 500 m above terrain on May 30, 2014. The scan angle rank ranged from -360 to 440. Full-waveform data was collected and discretised up to 7 returns per pulse. The scale factor of the data was 0.01 m.

In this study, the discrete return airborne LiDAR data of 28.66 points m^-2 were thinned to reduce the point density to 4.94 points m^-2 and spacing of 0.45 m (all returns) to simulate operational low-density LiDAR data typically available for forest surveys. This study did not consider field plot data and used an approach solely depending on remote sensing data, simulating operational LiDAR data from high-density LiDAR data for predicting VFS. Potapov et al. [57] used airborne LiDAR data as a proxy for field surveys and the reference for model calibration and validation. Mäyrä et al. [58] also suggested high-resolution airborne data for training for landscape-level species detection utilising spaceborne data. WorldView-3 satellite imagery acquired on October 05, 2015, was resampled from 1.6 m for the visible near-infrared (VNIR) 8 bands to 7.5 m using pixel aggregation [59] to match the spatial resolution of the shortwave infrared (SWIR) 8 bands of the same WorldView-3 imagery. We then used all 16 bands for further analysis. Similarly, all 16 bands of WorldView-3 imagery were again resampled to 30 m pixels to match the spatial resolution of the Landsat 30 m (OLI) data acquired on October 21, 2014 (See Table 1 for the detailed specifications of multispectral remote sensing data). Landsat-8 (OLI) Collection 2 Level 2 surface reflectance imagery was acquired from the Earth Explorer of the United States Geological Survey (USGS) portal. WorldView-3 imagery was resampled to examine the impact of pixel size and to scale up the VFS models to larger pixel sizes. The underlying rock types of the study site were extracted from a geology data layer based on 1:25,000 scale mapping [56]. This study focused on the fusion of multispectral satellite imagery and DTM derivatives to predict VFS following the workflow presented in Figure 2.

2.3. Data Pre-Processing

The LiDAR data pre-processing included noise removal, derivation of three DTMs, normalisation, and extraction of densities from three vegetation layers for each of the 1.6 m, 7.5 m, and 30 m grid cell sizes matching the spatial resolution of the visible and NIR bands of Worldview-3 (1.6 m), shortwave infrared bands of Worldview-3 (7.5 m) and Landsat OLI imagery (30 m) respectively. WorldView-3 satellite imagery required radiometric calibration, atmospheric correction, and extraction of spectral indices and texture features from each band. Similarly, the spectral indices and texture features were extracted from each band of Landsat-8 (OLI) surface reflectance imagery. To examine the effects of scales on the prediction accuracies, we conducted this study on the three spatial resolutions, combining the derivatives of LiDAR and multispectral satellite imagery, i.e., WorldView-3 and Landsat-8 datasets.

First, we clipped LiDAR point clouds for the study site, classified them into the ground and non-ground classes using lasground, and then normalised them using lasheight of LAStools software. The normalised data was filtered to provide smooth data for further analysis. Finally, we derived height density rasters for three vegetation strata using lascanopy to use them as response variables. For this study, we classified the vertical structure of forest areas into three vegetation layers: understorey or a lower vegetation layer (≥ 2 m to ≤ 10 m), mid-storey or a middle vegetation layer (> 10 m to ≤ 30 m), and overstorey or an upper vegetation layer (> 30 m to ≤ 50 m), based on expert knowledge of the vegetation communities [60,61]. Next, the original LiDAR point clouds were thinned to simulate operational LiDAR data (approximately 5 points m-2), and the three DTMs (1.6 m, 7.50 m, and 30 m) were generated using blast2dem. LiDAR returns below 2 m above the ground layer were discarded, avoiding low ferns, small shrubs, and coarse woody debris [9,62,63].

Like other multispectral remote sensing data, we processed the raw WorldView-3 satellite image to convert from digital number values to top-of-atmosphere spectral radiance and then to surface reflectance before performing further analysis [64]. We performed radiometric calibration using the information given in Table 2. We then performed atmospheric correction on the imagery to achieve the top of atmospheric reflectance by applying a QUick Atmospheric Correction (QUAC) model. No corrections were performed on the Landsat-8 (OLI) surface reflectance imagery, which could be directly used as input to the biophysical models [65]. All the datasets were registered and reprojected to GDA94/MGA zone 55.

2.4. Response and Predictor Variables

We deployed three density grid layers derived from high-density LiDAR data for the response variables, representing the LiDAR point density in the lower vegetation layer (≥ 2 m to ≤ 10 m), a middle vegetation layer (> 10 m to ≤ 30 m), and an upper vegetation layer (> 30 m to ≤ 50 m). The point density values are taken as a proxy for the presence of a forest understorey, mid-storey, and overstorey, respectively. This definition of VFS was motivated by prior research conducted by White [67] using the same high-resolution LiDAR dataset for directly mapping the understorey in our study area. White [67] validated the use of canopy structure as a proxy for understorey type based on 40 field plots (20 m radius) that were located using a stratified random approach to capture the variation in forest structure in the region. White [67] found a strong correlation between LiDAR-derived canopy density values and understorey within the field plots he sampled. White [67] also checked for the potential blocking effect of the overstorey canopy on the understorey data and reported that the canopy had a minimal impact on the density of the understorey data. Our study was conducted on the same LiDAR data and the same study site as White's (2017) study. We, therefore, were able to build on these earlier findings and use point density as a proxy for understorey presence. To limit the potential bias caused by the overstorey, this study used the seven highest-intensity returns discretised from high-density full-waveform LiDAR collected from a low flying height (500 m above ground level), achieving a point density of >28 points m^-2. The discretised high-density data ensured that understorey vegetation was captured in sufficient detail.

We tested 13 schemes of different combinations of datasets for the predictor variables combining topographic attributes, texture features, spectral indices, and a geology data layer as listed below. We used twelve topographic attributes, i.e., slope, aspect, catchment area, solar radiation, profile curvature, plan curvature, convergence index, terrain ruggedness index, slope length and steepness (LS) factor, SAGA wetness index, topographic position index, and stream power index, extracted from a DTM. We derived eight texture features (Appendix, Table A1) using the grey-level co-occurrence matrix for each band of multispectral satellite imagery based on their ability to characterise vegetation structure [68,69]. Thus, we derived 128 raster layers of texture features (16 bands × 8) from WorldView-3 of 7.5 m, 64 layers (8 bands × 8) from Worldview-3 of 1.6 m, and 56 layers (7 bands × 8) from Landsat-8 (OLI) 30 m data. We also extracted fifteen spectral indices from both multispectral satellite imagery (Appendix, Table A2). Variable selection is mostly not required in non-parametric modelling [70]; however, as this study dealt with a large number of variables, we applied a variance inflation factor (VIF) to remove the collinear variables using a threshold of 5.0 [71,72]. We selected all the variables based on their appropriateness and their established applications in forestry-related studies [73,74]. The following 13 classification schemes were assessed:

• Spectral bands (B);
• Topographic attributes and geology (A+G);
• Spectral indices (I);
• Spectral bands, topographic attributes, and geology (B+A+G);
• Spectral bands and spectral indices (B+I);
• Topographic attributes, geology, and spectral indices (A+G+I);
• Spectral bands, topographic attributes, geology, and spectral indices (B+A+G+I);
• Texture features (T);
• Spectral bands and texture features (B+T);
• Topographic attributes, geology, and texture features (A+G+T);
• Spectral indices and texture features (I+T);
• Topographic attributes, geology, spectral indices, and texture features (A+G+I+T);
• Spectral bands, topographic attributes, geology, spectral indices, and texture features (B+A+G+I+T).

In this study, we converted raster datasets into points and then sampled randomly from all the datasets without replacement. All the data covering the harvested sites, rivers, roads, and no data were removed. Out of the remaining point values, 10,000-point locations were randomly selected [75,76] and then divided into a model training dataset (50%) and an independent validation dataset (50%). From the training dataset, 70% of point locations were randomly drawn for model training and 30% for cross-validation [4,77,78]. To achieve robust and stable results for the models, the cross-validation was repeated 100 times [79,80] for predicting VFS. The predicted results were finally validated using the independent validation datasets.

2.5. Random Forest Modelling

We used random forest (RF) regression modelling to examine the predictive power of different datasets at different spatial resolutions for VFS. The RF model is capable of efficiently incorporating a large number of continuous and categorical variables [14]. The RF regression algorithm is a bagging technique which employs recursive partitioning to divide the input data into many homogenous subsets called regression trees (ntree) and then averages the results of all trees. Each tree is independently grown to its maximum size based on bootstrap samples from the training dataset (approximately 67%) without pruning. In each tree, RF selects a random subset of variables (mtry) to determine the split at each node [81]. The ensemble predicts the data that is not in the tree (OOB (out of bag data), approximately 33%), and by calculating the difference in the mean square errors between the OOB data and the dataset used to grow the regression trees. The RF algorithm allows us to identify important predictor variables to predict vegetation strata [74]. Variable importance is evaluated based on how much worse the prediction would be if the dataset for that variable were permuted randomly [82], and this can be used in feature selection by determining the importance of each variable in the regression process [83]. In this study, model performance was evaluated using the coefficient of determination (R²) e.g., [84]. Because the R² from RF modelling is derived from cross-validation, it provides useful comparisons between models with different numbers of predictors, unlike R² from standard linear models. The RF modelling analysis processes were conducted in the R programming language [85] using libraries randomForest [86], caret, and caTools [87].

We examined the importance of predictor variables to understand the relative importance of those used in this study. We utilised a percentage increase in mean square error (%IncMSE) as a metric for variable importance, which is one of the most widely used scores of importance [88]. The higher the %IncMSE value, the more important the variable [89].

3. Results

3.1. Model Accuracy Assessment

The RF model prediction accuracy for the 30 m pixels of Landsat-8 (OLI) data varied with the inclusion or exclusion of predictor variables, with mean R² (hereafter R²) values ranging from 0.15 to 0.65. The R² values increased with the number of predictor variables in most cases (Figure 3). Here, the RF model had the greatest predictive power for the overstorey (R² = 0.65) followed by the mid-storey (R² = 0.46), and the understorey vegetation density had the least predictive power (R² = 0.44). In this study, for the overstorey vegetation density, the combined scheme of B+A+G+I+T (scheme #13) produced the highest R² (0.65) followed by the scheme A+G+I+T (#12) dataset (R² = 0.63). Spectral indices (R² = 0.48) and texture features (R² = 0.47) were also useful for predicting the density of the overstorey vegetation of wet eucalypt forests. If the textures and spectral indices are combined, this dataset produced a higher model accuracy (I+T; #11; R² = 0.56) than either dataset alone.

In contrast to the overstorey vegetation density, the R² values for the mid-storey did not increase with the increase in the number of predictor variables when comparing the scheme B+A+G+I+T (#13; R² = 0.45) with A+G+I+T (#12; R² = 0.46) which yielded the highest R² for the mid-storey vegetation density.

For the understorey vegetation density model, the B+A+G+I+T scheme had equivalent explanatory power to A+G+I+T (both schemes' R² = 0.44). Overall, the R² values of combined datasets were mostly greater than for a single dataset. For example, R² values of spectral indices (I; #1; R² = 0.17) and texture features (T; #8; R² = 0.25) were lower compared to the scenarios with combined datasets (I+T; #11; R² = 0.27), although the gain from combining them was less than for the overstorey vegetation density.

For the models in the schemes with the resampled 30 m WorldView-3 datasets, the R² values ranged from 0.19 to 0.55 (Figure 4). Unlike the 30 m Landsat-8 data, the mid-storey vegetation density was the least predicted in most scenarios. The scheme of B+A+G+I (#7) explained the greatest percentage of variance (R² = 0.55) in the overstorey vegetation density, followed by that of B+A+G+I+T (#13; R² = 0.53), and the model using spectral indices (I) had the lowest explanatory value (R² = 0.19). For the mid-storey vegetation density, the combined spectral indices and texture features (I+T; #11) explained the highest percentage of variation (R² = 0.38), whereas the B+A+G+I+T explained the highest variance (R² = 0.39) in the understory vegetation density, closely followed by B+A+G (#4; R² = 0.38).

Comparing the schemes of the 30 m Landsat-8 (OLI) with the resampled 30 m WorldView-3 dataset, the WorldView-3 imagery could predict the understorey vegetation density better than the mid-storey vegetation density. In contrast, the schemes with the Landsat-8 (OLI) dataset could better predict the mid-storey vegetation density than the understorey vegetation density. Overall, at a 30 m scale, the Landsat-8 data had better predictive power than the WorldView-3 imagery for all three vegetation strata.

For 7.5 m pixel sizes with WorldView-3 imagery, combining the resampled VNIR bands and original SWIR bands, the predicted R² values ranged from 0.02 to 0.32 (Figure 5). The overstorey vegetation density had the highest predictive power (R² = 0.32), followed by the understorey vegetation density (R² = 0.21), and the mid-storey vegetation density (R² = 0.14) had the least predictive capacity. Combinations of two or more datasets provided better predictions than any of the individual datasets in this case.

For the 1.6 m pixel size of the original VNIR eight band WorldView-3 imagery and LiDAR-derived topographic attributes, the overall result demonstrated poor predictive power, with R² values ranging from 0.0 to 0.08 (Figure 6). The overstorey vegetation density had a stronger predictive capacity than the mid-storey and understorey vegetation strata. Thus, we inferred that the WorldView-3 imagery with a pixel size of 1.6 m and their combinations were inappropriate for predicting the vertical strata of wet eucalypt forests.

3.2. Model Validation

We evaluated the prediction model accuracies using an independent validation dataset to confirm the robustness and transferability of the acquired models (Figure 7 for Landsat-8 and Figure 8 for WorldView-3). In this study, the differences between the predicted and validated R² values (validated–predicted) were small (<5.67%) across all models. We considered two criteria for robustness and transferability: (1) the model with higher prediction and validation accuracies, and (2) the smaller the difference between validated and predicted R² values, the more robust the model. This study indicates that the developed models are robust and, thus, appropriate for transfer to other areas in the forest landscape.

3.3. Importance of Predictor Variables

We present the top 20 important predictor variables using Landsat-8 (OLI) and WorldView-3 datasets (1.6 m, 7.5 m, and 30 m) combined with topographic attributes, spectral indices, and texture features. For predicting the overstorey vegetation density, the Normalised Difference Water Index (NDWI) (%IncMSE = 104.0) was the most important variable, followed by geology (%IncMSE = 96.3) using the Landsat-8 (OLI) dataset (Figure 9). Surface Water Capacity Index (SWCI) (%IncMSE = 66.4) ranked as the most important predictor variable, followed by Normalised Burn Ratio (NBR) (%IncMSE = 41.5) using the resampled 30 m WorldView dataset (Figure 10). Saga Wetness Index, GNDVI, and LS factor were important predictors for these models. Other variables for the overstorey vegetation density, for example, SWI, CI, SIMI, and SPI using the Landsat-8 (OLI) dataset, and GNDVI, RDVI, TRI, and SIMI using the 30 m WorldView dataset produced lower variable importance scores.

For the mid-storey vegetation density, a variance of band#4 (%IncMSE = 90.2) was the most important predictor, followed by NDWI (%IncMSE = 65.5) and then geology (%IncMSE = 45.0) using Landsat-8 (OLI) dataset. Solar radiation (%IncMSE = 65.9) was the most important predictor variable followed by SATVI (%IncMSE = 28.4) and then SWCI (%IncMSE = 27.9) using the Worldview 30 m dataset. Similarly, for the understorey vegetation density, geology (%IncMSE = 58.0) was the most important predictor, followed by LS factor (%IncMSE = 40.6) and then solar radiation (%IncMSE = 36.7) using Landsat-8 (OLI) dataset. In contrast, SWCI (%IncMSE = 55.5) was the most important predictor variable followed by GARI (%IncMSE = 46.1) and then band#1 (%IncMSE = 29.1) using the Worldview 30 m dataset. For the WorldView dataset of 7.5 m resolution, band#16 (SWIR band#8) (%IncMSE = 46.7) and WorldView band#1 (%IncMSE = 34.7) produced the highest variable importance scores for the overstorey and understorey vegetation strata respectively. The 1.6 m resolution WorldView dataset provided the mean of band#1 (%IncMSE = 27.3) for the understorey, geology (%IncMSE = 15.2) for the overstorey, and again geology (%IncMSE = 8.9) for the mid-storey vegetation density as the most important variables.

4. Discussion

Comparing all the thirteen schemes, the spatial resolution of the freely available 30 m Landsat-8 (OLI) dataset and its combinations with the topographic attributes derived from the resampled LiDAR and geology data produced the best result (R² = 0.65) for the overstorey vegetation density followed by the mid-storey (R² = 0.46) and then the understorey vegetation density (R² = 0.44) (Figure 3). Models based on Landsat-8 (OLI) data combinations were consistently better than those found on the WorldView-3 dataset. The highest performing WorldView-3 model was produced by the 30 m resampled WorldView-3 dataset with topographic attributes and geology (R² = 0.55) for the overstorey vegetation density (Figure 4). Across datasets, predictive accuracy was generally increased with an increasing number of predictor variables. The models that included pixel-values of spectral bands from the multispectral datasets, topographic attributes derived from the resampled LiDAR data, spectral indices, texture features, and ancillary geology vector data produced the best overall prediction accuracy (R² = 0.65 for the overstorey vegetation density), providing higher model performance than the combinations with the spectral indices and texture features in isolation. Although the relationships between image texture metrics and forest structure attributes can be used to characterise complex forest structures and enhance vegetation properties, the uses of vegetation indices, particularly in the case of closed canopies, are difficult and challenging [90].

In this study, the spatial resolution of the 1.6 m WorldView-3 dataset was not considered appropriate for modelling vertical forest structure, as the R² values ranged from 0.0 to 0.08 (Figure 6). The results demonstrate that the resolution of remote sensing datasets was inversely proportional to the prediction accuracy for three vegetation density layers. The poor performance of the 1.6 m resolution WorldView-3 dataset is likely related to the mismatch between the high spatial resolution of the satellite pixels and the spatial scale of the response variables. Since the LDAR points were divided into three canopy layers and the ground, the 1.6 m spatial resolution dataset might have produced the lower prediction accuracy due to lack of LiDAR points. Overstorey, mid-storey, and understorey vegetation layers are difficult to capture and quantify at 1.6 m spatial resolution.

Despite focusing on a different forest type, our results were in broad accordance with those of Latifi, et al. [91], who used RF modelling of Landsat-5 Thematic Mapper (TM) imagery with LiDAR data to predict forest structural attributes. They also reported that accuracy was increased with the number of predictor variables and emphasised the utility of topographic variables derived from LiDAR data over the multispectral satellite data.

The optimal pixel size observed here (30 m) is likely to be related to tree crown size, with Azaele, et al. [92] and Cohen and Spies [45] finding that the spatial resolution of the dataset should match the object to be predicted. Cohen and Spies [45] postulated that the pixel size of Landsat TM data is roughly equivalent in size to the tree crowns. In contrast, a large tree may appear in many pixels if pixel sizes are smaller, for example, 7.5 m or 1.6 m, as tested in this study, and this consequently produced lower model accuracy. Munsamy, et al. [93] tested resolutions from 1 m to up to 9 m and found that coarser spatial resolutions produced the best results for predicting dominant and mean tree height in mature stands in South Africa. Therefore, a study requiring information on individual tree canopies may be expected to require a larger pixel size than studies requiring information on individual branches or smaller plants in lower canopy layers. We infer to do more research on why the 30 m pixel size performed best in our study, presumably due to using indirect predictors, and some of the key predictors were at a coarser scale regardless of the vegetation strata, with individual plants in the mid-storey and understorey layers having canopies far smaller than 30 m in diameter. Kamal, Phinn and Johansen [26] applied original and resampled WorldView-3 data and demonstrated that a pixel size ≤ 2 m was appropriate for mapping inter-canopy features, such as canopy gaps, and a pixel size ≥ 4 m was more suitable for mapping vegetation formation and communities in mangrove forests. In contrast to our results, their prediction accuracy improved as the pixel size decreased. If information is required on forest species composition as well as vertical structure, our previous research has demonstrated the utility of combining LiDAR with hyperspectral data to identify the canopies of overstorey and mid-storey trees [94]. It is important to note that we could not compare direct canopy density estimates for operational LiDAR with those of high-resolution LiDAR. Operational LiDAR was unavailable for much of the study area, and it was deemed inappropriate to compare direct density estimates from the down-sampled high-resolution LiDAR since this subset would have included the same points as for the comparison dataset, artificially inflating their explanatory value. Therefore, this study focused on combining multispectral satellite data and topographic attributes derived from the resampled LiDAR data.

The Landsat-8 (OLI) dataset produced the best overall results for each vegetation layer. However, the 30 m Landsat-8 (OLI) and the 30 m resampled WorldView-3 datasets produced contrasting outcomes for relative abilities to predict the mid-storey and understory vegetation densities. The original and resampled WorldView-3 data could predict understory vegetation density better than the mid-storey vegetation density, whereas Landsat-8 (OLI) data predicted mid-storey vegetation density better than understory vegetation density. This study recommends that when high-density LiDAR data are unavailable, using freely available Landsat-8 imagery combined with geology and DTM derivatives (topographic attributes), texture measures and vegetation indices can be adequate for forest managers and planners to assess the VFS of wet eucalypt forests, to contribute to sustainable management, planning, and monitoring of forests.

In general, the predicted accuracy of the models was compared favourably with those of previous research in other forest systems and confirmed the validity of the approach of this study. A study conducted by Zald, Wulder, White, Hilker, Hermosilla, Hobart and Coops [22] advocated that prediction accuracy depended on structural variable type using derivatives from LiDAR and Landsat data and presented R² ranging from 0.0 to 0.77, with the highest value for structural variables, including live trees and the lowest for downed wood. They found that metrics derived from Landsat data improved the prediction accuracies of all the structural variables, although those improvements were less than the inclusion of LiDAR data. Similarly, Wallner, et al. [95] stratified forest plots based on forest types and showed R² values ranging from 0.37 to 0.63 for modelling stand structural attributes using RapidEye data, particularly R² values of 0.4 for stand density. Kayitakire, Hamel and Defourny [20] reported the model prediction with an R² value of 0.38 for stand density. Therefore, this research can act as a benchmark for including a broad range of explanatory variables; however, further research is required to translate this approach to other types of forest ecosystems and satellite imagery (e.g., Sentinel-2) with different topography and understorey conditions.

Our study validated the predicted accuracies with promising goodness of fit statistics (differences between the predicted and validated R² values were <5.67%), suggesting our best models could be applied to larger geographical areas and indicating the capability of LiDAR-derived metrics to predict vertical structural variables with reasonable accuracy. Bolton, White, Wulder, Coops, Hermosilla and Yuan [21] achieved negative model bias (-1.2 to -2.1%) and validated their research in other study systems to examine transferability. They suggested that the structural variability in the new systems might not have been captured in the training plots used in the model development, and the bias varied from -3.9% to –8.0%. Therefore, further research would be required to test the robustness and transferability of our models to other wet forest landscapes, with due consideration given to the forest age structure, management history (our models excluded previously harvested and wildfire-impacted forests), underlying geology, and rainfall.

The best predictor for overstorey density derived from the Landsat-8 (OLI) dataset with topographic attributes, indices, and texture features (scheme #13) was the Normalised Difference Water Index (NDWI) (%IncMSE = 104.0). The second-ranked variable was geology (%IncMSE = 96.3) for the overstorey vegetation density (Figure 9). For mid-storey vegetation density, the variance of band#4 (one of the texture features) ranked first, followed by NDWI. Similarly, SWCI (%IncMSE = 55.5) ranked first, followed by GARI (%IncMSE = 46.1) for the understorey vegetation density.

Li, Li, Li and Liu [72] reported that texture features from Landsat-8 showed a significant correlation with the aboveground biomass, and the variance of Landsat band#4 performed best in predicting biomass. Using spectral indices and textural variables in combination with Landsat bands, Halperin, et al. [96] compared the Landsat-8 (OLI) dataset with RapidEye for canopy cover estimation in Zambia. They found Landsat-8 (OLI) consistently performed better than RapidEye and soil data improved the model accuracy. They inferred that vegetation indices were the most important variables. Nasiri, et al. [97] evaluated predictor variables and reported that vegetation indices were dominant predictors in most of the models. Forest health, foliage density, chlorophyll content, water stress, and tree biomass influenced forest canopy cover, resulting in the determination of important variables [96].

The topographic attributes derived from a DTM could be deployed in combination with Landsat-8 (OLI) datasets, as their contributions were each high in the models, demonstrating the advantage of combining the multispectral data rather than relying solely on topographic attributes derived from a DTM [89]. The freely available Landsat dataset could be used for modelling and mapping larger geographical areas and is easily adopted by forestry professionals, bushfire managers, and ecologists who require stand-level details of VFS.

5. Conclusions

The aim of this study was to predict the density of three vegetation layers in a wet eucalypt forest in Tasmania, Australia, based on multispectral satellite imagery and topographic attributes derived from airborne LiDAR data. We applied a random forest machine learning approach to model the density of the forest understorey, mid-storey, and overstorey based on a range of spectral (15) and texture (8) features per spectral band derived from multispectral satellite data and 12 terrain attributes derived from LiDAR data. We tested the impact of spatial resolution by applying the random forest algorithm at three resolutions: 1.6 m, 7.5 m, and 30 m, corresponding to the visible bands of Worldview-3, the SWIR bands of Worldview-3, and Landsat 8 OLI, respectively.

The fusion of the derivatives from the 30 m Landsat-8 (OLI) satellite imagery and DTM derivatives (topographic attributes) and geology produced the best overall results (R² = 0.65 for the overstorey layer). The classification schemes utilising Landsat-8 (OLI) datasets outperformed models based on high-resolution WorldView-3 imagery.

Spectral indices and texture features were useful for predicting the overstorey vegetation density of wet eucalypt forests. If the texture features and spectral indices were combined with topographic attributes, the combined dataset produced higher model accuracy than either dataset alone.

The analysed output demonstrated that the resampled 30 m WorldView-3 data showed better model performance than at original pixel sizes of 7.5 m and 1.6 m. The schemes with WorldView-3 datasets with a pixel size of 1.6 m were inappropriate for predicting the density of three vegetation layers of a wet eucalypt forest. Landsat-8 outperformed WorldView-3 imagery for all three-pixel sizes and all three vegetation layers.

The results showed an increase in model accuracies with an increasing number of predictor variables in most schemes, illustrating the merit of combining topographic attributes, geology, texture features, and spectral indices.

The differences between the predicted and validated accuracies were less than 5.7% for all models. This indicates that the developed models are robust and could be transferred to similar forests outside the focal research landscape. Our findings also show that texture features and spectral indices are important variables that can be used for vertical forest structure modelling.

Further work would be required to validate this approach and extend the results to different forest types and age classes since this study was confined to mature wet eucalypt forests. It would be highly valuable to test this approach on the topographic attributes derived from the spaceborne Global Ecosystem Dynamics Investigation (GEDI) LiDAR data (footprints averaging 25 m in diameter) [98] combined with the Landsat-8 (OLI) and Sentinel-2 imagery for predicting and mapping VFS over larger geographic areas. Our study suggests that freely available spatial data, including Landsat, DEM, and geology, could be used as predictor variables, enabling forest managers and planners to adopt the approach presented in this study even where LiDAR data are unavailable.