Estimating Forest Inventory Information for the Talladega National Forest Using Airborne Laser Scanning Systems

1. Introduction

The ability to adequately characterize and assess forest structure with a high level of accuracy is not only important for the development of a reliable forest plan but is also informative for assessments that demonstrate the sustainable management of natural resources [1,2]. Within this context, an understanding of forest conditions, particularly growing stock (tree and stand volume), aboveground biomass, and basal area, is crucial for planners, forest managers, and landowners. These attributes directly influence the potential revenue and the potential habitat a forest can provide and facilitate opportunities for addressing other management objectives [2,3]. Furthermore, a regularly updated forest inventory is essential for monitoring the spatiotemporal dynamics of forest ecosystems over the length of a planning horizon. For example, an estimate of aboveground biomass can provide insights into the capacity of a forest to sequester carbon, which is considered a critical factor in addressing climate change. This issue may become more important in the future as the forestry sector faces increasing pressure to assess the ability and rate of forests to sequester carbon [4,5,6,7]. Additionally, estimates of biomass can provide valuable information for assessing forest health related to the outbreak of Southern pine beetle [8] and for assessing fire risk associated with fuel management [9].

Traditionally, an estimate of a forest inventory has heavily relied on labor-intensive and time-consuming field measurements. Timely field measurements may be limited in spatiotemporal coverage, may include sampling error, and may not be representative of large forested areas [10]. While measuring the diameter at breast height (dbh) of trees, identifying the tree species, and counting trees within a sample unit (e.g., plot or prism point) are relatively straightforward methods, estimating the height of each tree may require more effort and include more uncertainty [11], especially in natural, mountainous forest landscapes [7]. Endeavors meant to obtain a sufficient number of well-distributed sampling plots that properly represent an entire forest area remain challenging due to limited resources and accessibility. To account for these limitations, traditional field measurements have been complemented by products derived from remote sensing systems, which may help address spatial and temporal challenges in developing a forest inventory [3,12]. This approach is sometimes referred to as an enhanced forest inventory [13,14]. Within this framework various remote sensing systems have been demonstrated to provide relatively accurate and cost-effective forest information including the development of forest metrics [2,3,10,15,16] and greenhouse gas inventories [17].

Data acquired from Light Detection and Ranging (LiDAR) devices can assist in the development of forestry information. Unlike other remote sensing processes that only provide two-dimensional information, LiDAR can characterize three-dimensional forest structure to a certain extent based on the point density of the LiDAR data [5,16,18]. LiDAR can help with deriving height estimates of an object through the time interval between the emission of a pulse of energy by the LiDAR sensor and the moment that the reflected signal has been returned to the device [19]. This type of active remote sensing technology has advanced significantly in the last ten years, resulting in a diverse range of LiDAR systems that can be installed in satellites, mounted on airplanes and unmanned aerial vehicles (UAVs), or used as hand-held devices. Each type of system has advantages and disadvantages. For instance, one advantage of using satellite-based laser scanning platforms such as NASA’s Global Ecosystem Dynamics Investigation (GEDI) and Ice Cloud and Land Elevation Satellite 2 (ICESat-2) is that they can provide multi-temporal data for the entire Earth. For instance, Potapov et al. [20] was able to create a global canopy height model using GEDI data along with Landsat imagery. Additionally, Dubayah et al. [21] presented the estimate of mean biomass densities for every country covered by GEDI with 1 km resolution.

Terrestrial-based LiDAR platforms can be categorized as either static or mobile. Terrestrial laser scanners (TLS) capture point clouds by looking upward from ground level, so they are advantageous in capturing details of forest structure from under the canopy. Especially, static TLS platforms (consisting of a sensor, a tripod, and a GNSS receiver) produce the highest-quality point clouds among LiDAR systems [5,22]. Appropriately employed TLS can facilitate visualization of tree branches and leaves; therefore, these systems have enormous potential for assisting in the development of new allometric biomass and wood quality relationships [23]. Arseniou et al. [24], for example, was able to estimate woody aboveground biomass for urban and rural settings using a TLS platform to identify various tree parts other than the main stem. Nevertheless, TLS platforms have several limitations when employed for forest inventory and operational forest management purposes. One, an occlusion effect, is caused by features hidden or obstructed behind larger diameter trees within the point cloud data, and it poses a challenge when a fixed-position, single-scan approach is used [25]. In addition, the weight and size of a TLS platform can make the effort of moving between field measurement plots challenging. These limitations, along with the cost of the platform, may influence opinions of whether TLS is a practical alternative for collecting forest inventory [26,27].

Mobile laser scanning (MLS) platforms, which is mobile TLS platforms, capture point clouds by looking horizontally from a height near ground level, where they are held, making them advantageous for capturing details of forest structure from a perspective very similar to human-collected field measurements. MLS platforms can collect information while a person traverses sampling plots with the instrument held in the person's hand or carried in a backpack. Moreover, the development of point clouds that are georeferenced to a local coordinate reference system, using Simultaneous Localization and Mapping (SLAM) algorithms, reduces the need for GNSS. While numerous researchers have illustrated the applicability of MLS for diverse forest inventory tasks [7,10,28,29], these studies often focus on MLS data collected at the plot-level. Among them, Vatandaşlar et al. [7] estimated several forest attributes (tree counts, dominant height, basal area, dbh, stand volume, and relative density) within plots located in a near-natural forest landscape. Employing a handheld MLS platform, Vatandaşlar et al. [7] mapped every stem within each field measurement plot and estimated these attributes with RMSEs ranging between 4.5% and 16.4%. However, as with TLS platforms, the number of sampling plots employed will likely influence the accuracy of forest or stand estimates [23]. As with any measurement system, the accuracy of forest attribute estimations can be positively correlated with the number of plots measured [30]. For instance, the diverse forests of the Talladega Division of the Talladega National Forest in the southeastern US cover almost 93,694 ha [31], and thus the number of TLS or MLS plots needed to describe forest character relatively accurately may be substantial. Therefore, an alternative solution needs to be sought to effectively and efficiently characterize the forest inventory of large, diverse areas such as this.

In this context, airborne laser scanning (ALS) systems have been considered a suitable choice for helping to describe the forest character of broad areas. Notably, the ALS data acquisition process is not constrained by the accessibility restrictions related to TLS including static and mobile platforms [5]. The versatility of an ALS system allows the collection of information across diverse temporal and spatial scales [32]. Recent advancements in sensor technologies have further encouraged the adoption of ALS systems, allowing the development of regularly updated and increasingly dense point clouds. For instance, a LiDAR-based forest inventory effort conducted in Ontario two decades ago resulted in a point cloud dataset with about 0.5 points per m², yet today the development of a point cloud dataset above 40 points per m² can be obtained [33].

While forest attributes such as tree heights and canopy coverage can be directly estimated from ALS LiDAR data, other forest characteristics such as aboveground biomass, growing stock (tree volume), and basal area can be inferenced from LiDAR-derived metrics [5]. The development of these estimates relies on modeling methods, which can range in complexity from regression to random forest models and other machine learning techniques [34]. Distinct models for estimating characteristics of different forest types (conifer, broadleaved, and mixed) might also be developed, rather than a general model that is applicable to an entire forested area. Further research in this domain is necessary, with a specific emphasis on investigating the nature by which additional spectral data can enhance the predictive capability of LiDAR point clouds. Additionally, as mathematical techniques and remotely sensed data evolve, the most effective combination of methods and data sources needs to be assessed.

A map of forest characteristics for an extensive forest area is an ideal outcome of remote sensing methods, yet this outcome is complicated by two underlying factors: the multitude of potential independent variables that can be derived from remotely sensed data, and the potential correlation amongst these which can induce a multicollinearity problem [35,36]. Furthermore, a large number of independent variables within predictive models can challenge the application of these models for developing broad scale GIS databases. Consequently, the selection of independent variables during the model development process is important. Tibshirani [37] suggested a method for developing linear models that estimate forest conditions, while enhancing prediction accuracy by reducing the number of independent variables. Adhikari et al. [35] recommended the use of ALASSO, as it effectively eliminated highly correlated independent variables from prediction models.

The main goal of this study is to develop models to estimate forest conditions across a broad area from information provided by ALS and aerial imagery. This research effort seeks to: (i) evaluate the quality of predictive models developed using the ALASSO method, (ii) evaluate whether additional remotely sensed data (multispectral aerial imagery) can enhance the quality of predictive models that rely on LiDAR point cloud data, and (iii) determine the suitability of general versus species group-specific models for characterizing mixed coniferous and deciduous forests located in the southern United States. We use the Talladega Division of the Talladega National Forest as our case study area because it represents typical characteristics of natural, pine-dominated forests of southeastern USA. National forest managers do not have comprehensive inventories of the extent under management limiting global decision space for at risk resources be it wildland fire, forest health, endangered species habitat condition to help assess where management is needed or has achieved desired future conditions. Therefore models, maps, and other outcomes derived from predictive models that are based on LiDAR (and other) data may provide forest managers, researchers, and policymakers with valuable insight to monitor and manage forests throughout the southeastern USA.

2. Materials and Methods

2.1. Study Area

This study was conducted within the boundaries of the Talladega Division of the Talladega National Forest, located in northwestern Alabama, USA. This 93,694-ha part of the national forest lies within the Piedmont and Ridge and Valley ecoregions [31]. The climate in this area generally consists of mild winters and hot summers, which are characteristic of a humid, subtropical climate. Elevation of the lands in the Talladega Division varies between approximately 160 m to 735 m above sea level, and annual precipitation is around 1260 mm. Historically, the forests in the study area have been composed of coniferous species such as longleaf pine (Pinus palustris Mill.), shortleaf pine (P. echinata Mill.), and loblolly pine (P. taeda), particularly in the uplands and on south-facing slopes. Deciduous tree species are often found in the riparian areas and on north-facing slopes. Oak (Quercus spp.), hickory (Carya spp.), maple (Acer spp.), and yellow-poplar (Liriodendron tulipifera) are among the more prevalent deciduous tree species in the study area.

2.2. Data Collection

2.2.1. Field Data Collection

A set of 254 fixed area, circular plots (402.6 m² with a 11.32 m radius) were measured by U.S. Forest Service crews between February and April 2022. These plots were pseudo-randomly located within the operable (upland) lands of the study area. An equal number of plots were measured within eight forest structure classes, defined by canopy density and tree height. Plot locations were limited to the interior (rather than edges) of management units (stands), and only one plot per management unit was allowed. Plot centers were mapped using a Trimble R1 GNSS receiver (Trimble, Colorado, USA). The field data collection procedure followed methods described in Laes et al. [38]. Within each plot, several characteristics of each live tree (height, dbh, tree species, tree crown status (dominant and co-dominant), azimuth, and distance from the plot center) were recorded. A dbh larger than 7.6 cm (3 inches) and a height greater than 0.6 m (2 feet) represented the minimum sizes of measured live trees.

2.2.2. Remote Data Collection

True color and near-infrared images of five counties that cover the study area were captured during the leaf-off season between 2020 (December) and 2021 (February). These images were used to create county mosaic orthophotographs with a spatial resolution of 0.3 m. The county mosaic images were acquired from the U.S. Department of Agriculture National Agricultural Imagery Program (NAIP) [39]. The ALS data was collected using a Leica Terrain Mapper. Within the study area, the ALS data contained 13.92 billion points. The average point density achieved was 9.7 points per m² and ranged from 3.5 points per m² to 20.3 points per m². The horizontal and vertical accuracy of the LiDAR data (0.71 m and 0.051 m, respectively) was estimated using 441 survey points spread across the State of Alabama. These characteristics satisfy the requirements for at least topographic quality level 2 (QL2) [40].

2.3. Data Processing

2.3.1. Field Data Processing

From the field measurement plots we estimated a variety of forest attributes related to tree heights (minimum, maximum, dominant, and mean merchantable height), tree dbh (minimum, maximum, quadratic mean of dbh, arithmetic mean of dbh, and coefficient of variation of dbh), canopy conditions (mean crown ratio), volume per unit area, and aboveground biomass per unit area. Among these, we focus on basal area (m² ha^-1), tree volume (m³ ha^-1), and aboveground biomass (Mg ha^-1) as dependent variables for the modeling effort. Basal area was estimated based on dbh, and tree volume and aboveground biomass were estimated using species-specific allometric equations from the U.S. Forest Service (National Biomass Estimator Library (NBEL) and National Volume Estimator Library (NVEL) which are Excel Add-ins developed by the Forest Management Service Center, U.S. Forest Service). Additionally, we classified sampling plots by the dominant tree species present, based on the total basal area for a species in a given plot exceeding 70% of the total basal area. As a result, there were only 14 oak dominant sampling plots and 149 pine dominated sampling plots (i.e., pine plots). All other plots indicated no dominance towards a specific tree species (i.e., mixed plots). Accordingly, we separately developed models based on two datasets (a) using all plots (n=254), and (b) pine plots (n=149).

2.3.2. Remote Data Processing

From NAIP imagery we created vegetation indices (i.e., greenness, Normalized Difference Vegetation Index (NDVI), and Enhanced Vegetation Index (EVI)) for the study area using ArcGIS Pro (Esri, Redlands, CA, USA). A total of 24 vegetation indices were developed and utilized as NAIP-derived independent variables (Table 1).

Due to the large file size of the raw LiDAR point cloud data, processing could not be completed with a desktop personal computer. Thus, the LiDAR point cloud data was processed using RStudio (R version 2023.06.1 Build 524) and the University of Georgia’s Georgia Advanced Computing Resource Center’s Sapelo2 Linux (64-bit CentOS 7.9) high-performance computing cluster. To reduce the processing time, parallel computing was used with 21 cores utilized simultaneously with a total of 900 GB RAM. The raw LiDAR data were first utilized to create a digital terrain model using only those points in the point cloud classified as bare ground. In creating the digital terrain model, selecting an appropriate algorithm and spatial resolution are important because it directly affects the result of estimated vegetation metrics [18]. Because of this, the spatial resolution was chosen based on point spacing and the number of points in the point cloud as suggested by McCullagh [41] and Hengl [42]. The triangular irregular network algorithm, a vector terrain model using Delaunay triangles, was applied to create a 1 m spatial resolution digital terrain model. With the digital terrain model representing the ground, we then generated a normalized LiDAR point cloud including only the above ground points. The first returns in the normalized point cloud were then used to create a canopy height model. Noise, resulting from points within the point cloud with negative values or values greater than 95% of the height, were removed to increase the quality of regression models. LiDAR-derived metrics were then calculated from point clouds that were clipped using the boundary of each sampling plot using R package (‘lidR’, version 4.0.3). We extracted 56 LiDAR-derived metrics from the normalized LiDAR point clouds for 254 of sampling plots. Among LiDAR-derived metrics, we excluded certain metrics with absolute values (such as maximum height and intensity, mean intensity, area, point counts, and total intensity), assuming that the LiDAR sensor was not calibrated for light conditions before data collection, and to avoid developing a model with variables that were return-density dependent. Therefore, a total of 74 independent variables were developed including 50 LiDAR-derived metrics (Table 2) and 24 NAIP-derived metrics which are used as independent variables in the modeling effort.

2.4. Modeling

We created 12 models to estimate three forest attributes (basal area, volume, aboveground biomass) with two different sets of field measurement sampling plots (all plots, pine plots) and two types of data sources (LiDAR only, LiDAR + NAIP). Depending on which field measurement sampling plots were employed, two sets of regression models were created: General models and pine models. Due to the high number of independent variables derived from LiDAR and NAIP data sources, inter-correlation between independent variables and the complexity of models was inevitable [3]. To reduce the complexity of models and to avoid multicollinearity issues, we applied the ALASSO regression. As ALASSO regression is a regularization technique, it addresses multicollinearity issues by shrinking regression coefficients to zero based on the lambda value which minimizes the sum of the squared differences between predicted values and observed values. The best lambda was obtained by applying 10-fold cross-validation. Prior to applying ALASSO regression, the normality of data was assessed using the Box-Cox method. Given the results of the Box-Cox method, a natural logarithmic transformation was applied to the dependent variables, improving the linearity of the data, the homogeneity of residual variances, and the normality of residuals. Additionally, potential outliers were investigated using Cook’s distance and studentized residuals. The observed outliers without leverage were eliminated. The best models were selected based on various parameters including the R²_adj., the number of independent variables, root mean square of error (RMSE), Collin Mallow’s C_p, Akaike’s information criteria (AIC), and Bayesian Information Criterion (BIC). The models were back transformed by exponentiating the dependent variable and were applied to estimate the forest attributes for the broader study area. Lastly, 10-fold cross-validation was applied to evaluate the quality of the models. While the modeling procedures were conducted using RStudio (R version 2023.06.1 Build 524) on a desktop computer, broad scale estimation maps were processed using the Sapelo2 Linux based high-performance computing cluster.

3. Results

3.1. Field-Based Forest Inventory

The average basal area, tree volume and aboveground biomass of the 254 field measurement sampling plots was about 23.4 m² ha^-1, 180.9 m³ ha^-1, and 40.1 Mg ha^-1, respectively (Table 3). For the pine dominated field measurement plots, the average basal area, tree volume, and aboveground biomass were about 22.3 m² ha^-1, 169.0 m³ ha^-1, and 34.6 Mg ha^-1, respectively. The range of values for the sampling plots was large since the plots were meant to be representative of all forest conditions in the study area, from early successional to mature forest stages.

3.2. Regression Models

3.2.1. Estimation of Forest Attributes Based on General Models

In nearly every case, the LiDAR + NAIP regression models included a larger number of independent variables than the LiDAR-only models. In the basal area model based on LiDAR + NAIP, five vegetation indices (G_MIN, NDVI_MIN, NDVI_MEDIAN, EVI_MAX, and EVI_PCT90) were additionally selected instead of the zq30 variable (30th percentile height distribution from the ground) found in the basal area model based on LiDAR-only (Table 4). The R²_adj. values for the basal area models were 0.72 (cross-validation: 0.69) for LiDAR + NAIP and 0.71 (cross-validation: 0.71) LiDAR-only, while RMSE values were 5.6 m² ha^-1 (cross-validation: 5.90 m² ha^-1) for LiDAR + NAIP and 5.7 m² ha^-1 (cross-validation: 5.91 m² ha^-1) LiDAR-only (Table 5). Other quality metrics except AIC also indicated that the model based on LiDAR + NAIP performed slightly better than the model based on LiDAR-only. The predicted basal area tended to be underestimated as its basal area increased regardless of data sources (Figure 2).

Figure 1. Locations of the field measurement sample plots within the study area (Sources: Esri. “World Topographic Map” [basemap]. January 31, 2024).

Figure 1. Locations of the field measurement sample plots within the study area (Sources: Esri. “World Topographic Map” [basemap]. January 31, 2024).

In the volume models, a higher R²_adj. value (0.77) was observed (Table 5). Although the overall quality metrics for LiDAR + NAIP and LiDAR-only models represented a similar level of performance, the increment in the number of independent variables in LiDAR + NAIP models was noticeable as it increased from 12 to 21 (Table 4). Specifically, eight independent variables derived from NAIP imagery were selected. Additionally, a LiDAR-derived metric, zpcum5 (cumulative percentile of return in the 5th layer), was also selected as an independent variable in LiDAR + NAIP model yet was not used in the LiDAR-only model. As observed in the basal area models, the difference between observed volume and estimated volume tended to increase with higher volume regardless of data sources (Figure 2).

The aboveground biomass models had R²_adj. values that were 0.73 and 0.72 (cross-validation: 0.64 and 0.65) for LiDAR + NAIP and LiDAR-only models, respectively (Table 5). Similar to the basal area models, the LiDAR + NAIP model yielded slightly more accurate RMSE values than the LiDAR-only model (LiDAR + NAIP model: 11.7 Mg ha^-1 (cross-validation:13.05 Mg ha^-1); LiDAR-only model: 11.8 Mg ha^-1(cross-validation:13.09 Mg ha^-1)). Nevertheless, the difference between aboveground biomass models depending on data sources was slight. Five vegetation indices were added into the LiDAR + NAIP model as independent variables (Table 4). Both aboveground biomass models underestimated biomass for the sample plots having higher aboveground biomass, suggesting estimation errors increase as forest stands mature (Figure 2).

The average R²_adj. and average R² values of 10-fold cross-validation results ranged from 0.64 to 0.73 and from 0.68 to 0.75, respectively (Table 6). The difference between statistics of cross-validation and developed models were generally smaller in LiDAR based models compared to the LiDAR + NAIP models. Regarding the selection of independent variables, several metrics were commonly selected regardless of data source, including: pzabove2, zq95, zpcum6, isd, iskew, ikurt, ipcumzq90 and p2th (all from LiDAR data); and G_MIN, NDVI_MIN and EVI_PCT90 (all from NAIP data) (Table 7).

3.2.2. Estimation of Forest Attributes Based on Pine Models

The R²_adj. and RMSE values of basal area models were 0.81 and 0.80 (R²_adj. of cross-validation: 0.79 and 0.75), and 4.8 m² ha^-1 and 5.1 m² ha^-1 (RMSE of cross-validation: 5.21 m² ha^-1 and 5.71 m² ha^-1) for LiDAR + NAIP and LiDAR, respectively (Table 8). Noticeable improvements in quality metrics of LiDAR + NAIP models were not observed by supplementing independent variables using NAIP derived metrics. Interestingly, there was no increase in the number of independent variables, but three NAIP derived metrics (NDVI_MIN, NDVI_MEDIAN, and EVI_PCT90) were selected instead of LiDAR derived metrics (zpcum6, ipcumzq30, ipcumzq90) (Table 4). The visual interpretation of the scatter plots showed that basal area tended to be underestimated as basal area values increased regardless of the data sources (Figure 3). Similar trends in scatter plots were also observed in general models, but the distribution of points was closer to the 1:1 line.

Regarding the volume models, R²_adj. values of 0.84 (cross-validation: 0.78) for LiDAR + NAIP and 0.82 (cross-validation: 0.79) for LiDAR-only models achieved (Table 8). The RMSE values ranged from 37.9 m³ ha^-1 to 43.5 m³ ha^-1 (cross-validation: from 47.44 m³ ha^-1 to 48.94 m³ ha^-1) which were lower than comparable values of the general regression models. However, AIC and BIC of general volume regression models were lower than those of the pine volume regression models. The increment in number of independent variables in LiDAR + NAIP models was notable as it increased from 11 to 30 in volume models. Specifically, eight independent variables derived from NAIP images were added to the LiDAR + NAIP regression model (Table 4). Additionally, 12 independent variables derived from LiDAR point clouds were added on the LiDAR + NAIP regression model.

The aboveground biomass models had R²_adj. values of 0.83 (cross-validation: 0.80) for LiDAR + NAIP and 0.82 (cross-validation: 0.78) for LiDAR-only models (Table 8). The RMSE values of aboveground biomass models were 7.9 Mg ha^-1 and 8.9 Mg ha^-1 (cross-validation: 9.4 Mg ha^-1 and 10.2 Mg ha^-1) for LiDAR + NAIP and LiDAR-only, respectively. Regarding AIC and BIC values, our results indicated that the quality of LiDAR-only models was better than the LiDAR + NAIP models. The number of independent variables also noticeably increased in LiDAR + NAIP models (LiDAR + NAIP: 34; LiDAR-only: 10). Ten NAIP derived metrics were selected in addition to 16 LiDAR derived metrics in LiDAR + NAIP model (Table 4). Similar to the general regression models, both models (LiDAR + NAIP and LiDAR-only) underestimated aboveground biomass for the sample plots having higher biomass estimates, suggesting estimation errors increased as stands mature (Figure 3).

The average R²_adj. values of 10-fold cross-validation results of pine models ranged from 0.75 to 0.80 regardless of forest attributes and data sources (Table 9). The average R² values of 10-fold cross-validation results for the pine model ranged from 0.78 to 0.84. Developed models based on LiDAR + NAIP for basal area and aboveground biomass were more robust than developed models based on LiDAR-only. Otherwise, the developed model based on LiDAR + NAIP for the volume less robust that developed models based on LiDAR-only. Similar to the general models, there were some LiDAR and NAIP derived metrics which were commonly selected as independent variables in every pine model regardless of forest attributes: pzabove2, zq5, zq95, zpcum5, iskew, ikurt, and p2th (all from LiDAR data); and NDVI_MIN, EVI_PCT90, and NDVI_MEDIAN (all from NAIP data) (Table 7).

4. Discussion

Estimating forest attributes such as basal area, volume, and aboveground biomass and updating this information regularly are critical activities for forest management and planning efforts. This study attempted to improve the performance of forest attribute estimations for large areas using LiDAR point clouds and high-resolution, multispectral remotely sensed data. We investigated the effect of different combinations of remotely sensed data (LiDAR-only or LiDAR + NAIP) on the quality of regression models. Also, we evaluated the quality of regression models depending on the classification of sampling plots according to species mixture. To avoid the overfitting issue resulting from multicollinearity between independent variables, the ALASSO method was employed during the modeling process as this method performance was suggested by earlier studies [35,43]. Eventually, a total of 12 models were developed for three forest attributes (basal area, volume, aboveground biomass) based on the species mixture (all species and pine), as well as data source (LiDAR-only, LiDAR + NAIP). When the general models are applied to the case study landscape, broad scale maps of these resource conditions can be visualized (Figure 4, Figure 5 and Figure 6).

The R²_adj. values of developed models ranged from 0.71 to 0.84 which were comparable to other ALS study results [3,34,44]. While many researchers have developed regression models for estimating forest attributes over relatively small forested areas [10,15,19,34,36], there have been few attempts to work in large areas such as U.S. National Forests. For instance, Leboeuf et al. [44] mapped a merchantable wood volume of very large area (440,000 km²) using airborne LiDAR. However, this study had a few limitations, including significant temporal differences between the field measurements (2003-2018) and LiDAR data (2011-2020) collection periods, limited spatial resolution, and poor representativeness of sampling plots. In our study, however, field and remote sensing datasets were collected during a relatively similar period of time (2020-2022). The distribution of sample plots is also important in developing robust regression models for forest attribute estimation. To enhance the quality of regression models, we used a stratified pseudo-random sampling design by classifying the entire operable study area with consideration of recent management activities to improve the balance and representation of the heterogeneity of the forest in the sampling design. Additionally, we surveyed all tree species (both merchantable and non-merchantable > 7.62 cm in dbh) within sampling plots to enhance the accuracy of regression models, as suggested by Brown et al. [45].

The R²_adj. values and R² values of pine models were higher than those of the general models regardless of forest attributes. The overall quality metrics also indicated pine models were higher compared to the general models. Bouvier et al. [3] also confirmed that separate models may result in higher accuracy when compared to general models. Regarding forest attributes, the highest R²_adj. values were observed in tree volume models and the lowest R²_adj. values were observed in basal area models regardless of data sources and sampling plots (all plots or pine plots). This trend has also been observed in other studies using ALS systems [3,34,46,47]. Sumnall et al. [48] suggested that basal area models have relatively lower R² values than models estimating tree height and biomass. This is likely because, unlike tree height, dbh cannot be directly measured with ALS systems. Dbh is used to estimate basal area, yet dbh and tree height are often needed to estimate volume and aboveground biomass [3], therefore it may be logical to observe that the basal area estimation models would underperform the volume and biomass models.

There have been many attempts recently to enhance the performance of regression models for estimating forest conditions. They generally involve the development of LiDAR metrics [3,49], and perhaps supplementary remotely sensed data [34,45]. In this study, the vegetation indices derived from NAIP imagery were included as supplemental data. As the additional NAIP imagery had higher spatial resolution (0.3 m), we expected the vegetation indices might improve the quality of regression models to a considerable extent. However, it was observed that the addition of NAIP-derived vegetation metrics was not very influential in improving the quality of prediction models. Although LiDAR-only derived models had slightly larger RMSE values and slightly lower R² values compared to the LiDAR + NAIP derived models, the increase in the number of independent variables with the addition of NAIP metrics led us to conclude that LiDAR-only models were more appropriate for broad-scale mapping efforts.

Of the 74 metrics derived from LiDAR and NAIP data sources, 45 were selected as independent variables in at least one regression model. Further, five LiDAR-based metrics and two NAIP-based metrics were selected for every regression model (Table 7). These were, specifically, the LiDAR-based metrics pzabove2, zq95, iskew, ikurt, and p2th, and the NAIP-based metrics NDVI_MIN and EVI_PCT90. LiDAR-based metrics pzabove2 and p2th metrics help eliminate the effect of understory vegetation that is not measured in a typical forest inventory survey. The LiDAR-based metric zq95 is widely used to represent forest canopy height. In addition to height related LiDAR metrics, intensity related metrics, such as iskew and ikurt, were also crucial in developing regression models as they provide information related to stand density. The inclusion of NAIP-based NDVI_MIN and EVI_PCT90 can be explained by the concentration of active chlorophyll in pine tree crowns. Ozkan et al. [34] also confirmed that metrics such as these may be significantly correlated with forest attributes.

Although inclusion of NAIP-based metrics improved the accuracy of the regression models we developed, we suggest that LiDAR-derived metrics may be sufficient for developing robust regression models used for operational forest management purposes. Interestingly, the general trend in the spatial pattern of basal area, volume, and aboveground biomass is due to the high correlation among these forest attributes (Figure 4, Figure 5 and Figure 6). Since aboveground biomass of a forest is typically a derivation of growing stock, volume estimates often rely on dbh measurements, and as basal area is directly related to dbh, these relationships may help in developing simple ratios between forest conditions without the need for additional, elaborate regression models. While this may potentially decrease the reliability of some of the outcomes, the practicality of the overall effort may be worth investigating further.

5. Conclusions

The objective of this study was to estimate and map key forest attributes across a wide area of interest, utilizing ALS data and aerial imagery. To achieve this objective, we employed the ALASSO modeling method, relying on 254 field measurement plots from a national forest in Alabama (USA), a mostly natural ecosystem composed of coniferous and deciduous tree species. One of the main conclusions drawn from our findings was that the LiDAR data collected by the ALS system with topographic quality level QL2 character seems to be a sufficient input for the development of regression models that can estimate basal area, volume, and aboveground biomass at accuracy levels that are acceptable for operational forest inventories. A second conclusion was that the added value of using optical data (aerial imagery) and associated vegetation indices as inputs for the development of regression models for estimating basal area, volume, and aboveground biomass was negligible, considering the increased model complexity and extra time required to process and analyze the additional model inputs. A third conclusion was that the models exclusively developed for areas dominated by pine species seem to outperform (with respect to the coefficient of determination from cross-validation) general models that were developed for all tree species within the study area. While we affirmed these conclusions for natural, pine-dominated forests in the study area, further research was needed to assess the scalability of the results to other forest ecosystems of the southern USA.