Generalized models for subtropical forest inventory attribute estimations using a rule-based exhaustive combination approach with airborne LiDAR-derived metrics

ABSTRACT

Airborne LiDAR has been widely used to map forest inventory attributes at various scales. However, most of the developed models on airborne LiDAR-based forest attribute estimations are specific to a study site and forest type (or species), so it is essential to develop predictive models with excellent generalization capabilities across study sites and forest types for the consistent estimation of forest attributes. In this study, 13 LiDAR-derived metrics, which depicted the three-dimensional structural aspects of stand canopy and had clear forest mensuration and ecology significance, were categorized into three groups (height, density, and vertical structure). A rule-based exhaustive combination was then used to construct 86 multiplicative power formulations consisting of 2–5 predictors for estimating the stand volume and basal area. By calibrating and validating these formulations using data from four forest types in the three study regions, we obtained the 24 best local models. Based on these models we proposed a set of accuracy criteria to determine generalized formulations and models. By applying two selection methods (the mean and mixed data methods), we finally archived the eight best region-generalized models, which could be used for estimating the stand volume and basal area of four forest types across study sites on a province scale. This study highlights the accuracy criteria and procedures for developing generalized formulations and models for consistent estimations of forest inventory attributes using airborne LiDAR data.

KEYWORDS:

• Forest structure
• model generalization
• accuracy criteria
• multiplicative power model
• remote sensing

1. Introduction

Accurate, detailed, and reliable information on forest resources at various scales not only helps forest managers optimize forest management activities; but also helps government departments make strategic decisions. Airborne laser scanning (ALS; also referred to as light detection and ranging (airborne LiDAR)) has the ability to range and penetrate tree canopies accurately (Watt and Watt 2013) and capture precise three-dimensional (3D) structural information of the canopy (Coops et al. 2021). Since it became commercially applicable in the mid-1990s (Vauhkonen et al. 2014), airborne LiDAR has gradually developed as an advanced technological tool for forest inventory. In Scandinavian countries and Canada, airborne LiDAR has been applied to operational forest resource inventories since 2002 (Næsset et al. 2004; White et al. 2017). In some other countries and regions, airborne LiDAR has been applied in large-scale ecological monitoring and planted forest inventories (Watt and Watt 2013; Görgens et al. 2015; Silva et al. 2017; Matasci et al. 2018).

Regional airborne LiDAR forest attribute estimation was performed using an area-based approach (ABA) based on the statistical relationship between the predictor variables derived from airborne LiDAR data and the response variables measured from field plots (Næsset et al. 2004). In the past 20 years, numerous studies have been conducted on various forest types and forest attributes (mean stand height, H; mean stand diameter at breast height, DBH; basal area, BA; stand volume, VOL; aboveground biomass, AGB; etc.) and have produced numerous estimation models (Zolkos, Goetz, and Dubayah 2013; Latifi et al. 2015), including parametric regression and nonparametric models. Most of these studies have been devoted to finding the best prediction model, i.e. maximizing explained variability (e.g. R²), minimizing prediction error (e.g. RMSE), and minimizing systematic bias (Zolkos, Goetz, and Dubayah 2013) for specific forest attribute, forest type, and study site (Næsset et al. 2005; Hudak et al. 2008; Penner, Pitt, and Woods 2013; White et al. 2017). Parametric models or multiple linear regression (MLR) are the main tools for airborne LiDAR forest parameter estimation (da Silva et al. 2020; Novo-Fernández et al. 2019). Their main advantages are the simplicity and clarity of the resulting models, but they need to satisfy the assumptions of normality, homoscedasticity, independence, and linearity (García-Gutiérrez et al. 2015). Non-parametric models or machine learning models are data-driven models with the ability to capture the implied, potentially nonlinear, and complex relationships between the dependent and independent variables (García-Gutiérrez et al. 2014; Gajardo, García, and Riaño 2014) and do not depend on any priori assumptions about the data (Millie et al. 2012). The main problems with non-parametric models are their higher intrinsic dependence on training data (usually field measurement data) (Zaffalon 2005), the high risk of overfitting (Hamedianfar et al. 2022; Hawkins 2004), and lack of transferability (Zhao et al. 2017; García-Gutiérrez et al. 2014). Many studies have demonstrated that the estimation accuracies of machine learning models are slightly better than that of multiple linear regression models (van Ewijk et al. 2019; García-Gutiérrez et al. 2015). However, some studies have also demonstrated the superiority of linear regression models over machine learning models (Xu, Manley, and Morgenroth 2018; Gao et al. 2018; Watt et al. 2015; da Silva et al. 2020). This study is dedicated to the development of transferability models among study sites and forest types, so we focused on the parametric models. Given the flexibility of multiplicative models (Hollaus et al. 2009), the multivariate multiplicative power models were used to build generalization models in this study.

Multiple linear (generally requires a logarithmic transformation, it is a nonlinear model) and nonlinear regressions are the most commonly used models for estimating forest attributes (Görgens et al. 2015; Giannico et al. 2016; Montealegre et al. 2016; Silva et al. 2017). There are various approaches to selecting model variables, such as stepwise regression techniques (Means et al. 2000; Treitz et al. 2012; Ene et al. 2012; Montealegre et al. 2016), the random forest approach (Silva et al. 2017; Luo et al. 2018; van Ewijk et al. 2019), the Akaike information criterion (AIC) (Silva et al. 2017; Sheridan et al. 2015), the Pearson correlation coefficient, and empirical analysis (Zhang, Cao, and She 2017; Cao et al. 2019), all of which have advantages and limitations. No matter which method is used, the selected model variables varied greatly across study areas, forest types, and stand conditions (Xu, Manley, and Morgenroth 2018; Görgens et al. 2015; Giannico et al. 2016; Montealegre et al. 2016; Silva et al. 2017; Sheridan et al. 2015; Treitz et al. 2012; Ene et al. 2012; Zhang, Cao, and She 2017; Kim et al. 2016). Thus, most models are not generalizable and transferable. During the pre-experiment of this study, we performed ten tests using a rule-based combination of 13 LiDAR variables, each with 90% of the plots randomly selected from the same dataset for model calibration. The final results obtained six different model formulations. This experiment indicated that, even with minor changes in the data, the model variables were different and the model structures were unstable. This situation is common to both MLR and machine learning; even though the model, forest type, and forest parameters are the same, the model variables are different in different study areas. (Singh et al. 2015; Gleason and Im 2012; Liu et al. 2021). Therefore, constructing formulations and models applicable to various regions, forest types, and stand structures—i.e. generalized or transferable models – the models built with data from one region, one forest type, or one year are transferable to other study areas, other forest types, or another year – is a current challenge (Knapp et al. 2020; Bouvier et al. 2015; Fekety, Falkowski, and Hudak 2015). Another major problem with existing LiDAR-based models for estimating forest attributes is that some LiDAR-derived variables, such as the (all_return_above_3 m)/(total_1st_return)×100 (Xu, Manley, and Morgenroth 2018), 95% density percentile (Montealegre et al. 2016), and cubic mean height (Görgens et al. 2015), have no clear forest mensuration or ecological significance. This problem is also present in machine learning models (Silva et al. 2017; Pearse et al. 2017).

Based on the relationship between AGB and DBH in field plots, Asner et al. (2012) proposed a generalized model for estimating forest aboveground carbon density (ACD). This model was developed using data from four tropical forest study areas and has been applied successfully to different tropical regions (Asner et al. 2012; Coomes et al. 2017; Getzin et al. 2017). However, this model had only one LiDAR variable among its three variables; therefore, it was not a genuine LiDAR-based model. Knapp et al. (2020) categorized LiDAR-derived metrics depicting forest structure into four groups: height, density, maximum height, and canopy heterogeneity variables, and obtained optimal prediction models for AGB and BA containing 4–5 variables selected by different combinations of variables. Similar to the model of Asner et al. (2012), non-LiDAR variables, such as maximum stand density in the study area obtained from field plot measurements or expert knowledge, were used in the model of Knapp et al. (2020). However, the idea of model construction was scientific and imaginative. Bouvier et al. (2015) proposed an empirical model with the same model shape and the same four LiDAR metrics depicting the 3D structure of the stand canopy (mean point cloud height and its standard deviation, canopy cover, and coefficient of variation of leaf area density) for forest inventory attribute estimation. The model achieved satisfactory results in estimating different forest attributes for different forest types in France. Thus, it was a generalized model for forest attribute estimation constructed entirely by LiDAR metrics. Yet, questions remain on how to evaluate the generalization performance of a model, how to develop a generalized model, and how to construct a generalized model with high accuracy while following the principles of forest mensuration and ecology. All these issues require additional investigation.

The overall goal of this study was to investigate the procedures for developing generalized predictive models for forest inventory attributes using airborne LiDAR data that can be applied across study sites and forest types. To fulfill this goal, airborne LiDAR data and field plot data were collected from a large subtropical region. Three specific objectives were pursued: (1) to establish the accuracy criteria for determining the generalized models; (2) to develop the best region-generalized formulations and models for estimating the stand volume and basal area across study sites without causing prediction bias at any individual site; and (3) to develop the best type-generalized formulations and models for estimating the stand volume and basal area across study sites and forest types without causing prediction bias at any individual site and forest type.

2. Materials and Methods

2.1. Study area

The study area were the same as the another study of Li et al. (2022), it covered the entire Guangxi Zhuang Autonomous Region (104°28′-112°04′E, 20°54′-26°24′N) in southern China, with an area of 237.6 × 10³ km². The forest area covers 117.21 × 10³ km², which represents 49.3% of the total area of the study area. In this study, airborne LiDAR data acquisition and field plot measurements were carried out over three years in three regions based on financial allocations, namely, the Nanning region (with an area of 22.1 × 10³ km²), the eastern region (128.4 × 10³ km²), and the western region (87.1 × 10³ km²) (Figure 1).

Figure 1. Location of the study area. (a) Geographic location of study area in China; (b) distribution of field plots in three regions in study area; (c) locations of field plots.

The Tropic of Cancer traverses the central part of the study area, which is in a subtropical monsoon climate zone. The mean annual temperature is 17.6–23.8°C, and the mean annual rainfall is 1080–2760 mm, with 70%–85% of the rainfall concentrated from April to September.

Due to the extensive extent of the study area, the climate, topography, and forest composition of the three regions are different. The elevations in the Nanning region range from 50 to 1760 m above sea level. It has a mixture of mountains, hills, basins, plains, and karst landforms, with forests dominated by eucalyptus plantations (mainly Eucalyptus urophylla and E. grandis × E. urohpylla), natural pine forests (mainly Pinus massoniana), and natural broad-leaved forests, with large karst rocky mountain shrub forests and small planted Chinese fir (Cunninghamia lanceolata) forests. The eastern region has an elevation ranging from 5 m to 2145 m. The area is interspersed with mountains, hills, and plains. The north-central part is dominated by natural broad-leaved forests (typical evergreen broad-leaved forests), pine forests (approximately 80% are natural forests and 20% are planted forests), and planted fir forests, with some eucalyptus plantations. In the central, southern, and southwestern parts of the region, eucalyptus plantations and natural pine forests are dominant, with a few planted fir forests. The elevations in the western region range from 60 m to 2062 m. The northern and northwestern parts of the area are mountainous landforms with natural broad-leaved forests (mainly Quercus variabilis, Q. acutissima, Q. fabri, and Q. aliena), natural pine forests, and planted fir forests, with small eucalyptus plantations. The eastern and south-central parts of the region have a karst landscape with natural broad-leaved and shrub forests in the karst mountains, along with some natural pine forests and eucalyptus plantations. The southern part of the region has a mountainous landscape with natural pine forests and eucalyptus plantations, some natural broad-leaved forests, and a few planted fir forests. According to the 5^th Forest Management Inventory of Guangxi (GX-FMI, 2017–2020), the areas of the fir, pine, eucalyptus, and broad-leaved forests accounted for 16.5%, 17.5%, 24.8%, and 41.2% of the total forest area, respectively. Approximately 80% of the fir forests are pure forests, and the remaining 20% are mixed forests. Approximately 70% of the pine forests are pure forests, and 30% are mixed forests. Almost all the eucalyptus forests are pure forests, and almost all the broad-leaved forests are natural mixed forests.

2.2. Data collection

2.2.1. Field plot data

The forests in the study area were categorized into four types: fir, pine, eucalyptus, and broad-leaved. The cluster-distributed field plots were installed across the study area using the 2015 forest resource database and the ArcGIS 10.0 software based on the representativeness of the forest attributes, such as stand age, mean height, and stand volume. The smallest spatial interval between the plots was 500 m (200 m in the Nanning region). There were 107 clusters of plots, each with 5–30 field plots (Figure 1). A hand-held dual-frequency differential GNSS receiver was utilized to navigate to the stand at the plot location. The rectangular plot was 30 × 20 m in size, set in a north-south orientation, and was subdivided into four subplots with an area of 15 × 10 m. The boundaries of the plot and subplots were measured and set using a compass and a handheld laser rangefinder and marked by nylon ropes. The northwest and southwest corners of the plot were positioned using a Trimble Global Navigation Satellite System (GNSS) receiver (Trimble Navigation Ltd., Sunnyvale, CA, USA) based on a real-time kinematic (RTK) positioning method. Two RTK-GNSS instruments were used as the base stations, which were located in a nearby open area. Using the post-correction approach, the positioning accuracy was better than 1 m. Within each subplot, we measured the DBH of all live trees with a DBH (1.3 m) greater than or equal to 5 cm using a diameter tape and recorded the tree species. Tree height was measured for three average-height trees and the tallest tree in each subplot using a Haglöf Vertex IV hypsometer. Subplot-level forest attributes included DBH, H, BA, VOL, and stem density. The VOL was calculated using the provincial species-specific allometric equations (Liao and Huang 1986), which were based on BA and H as predictors. The following methods were employed to calculate the plot-level forest attributes: the BA and VOL were summed by the subplot-level BA and VOL, respectively, the DBH was calculated for all trees in the four subplots; and the H was calculated using the BA-weighted-average method. Moreover, the dominant species, forest origin, average age, and canopy coverage were measured and recorded at the plot level.

The field plots in the Nanning, eastern, and western regions were measured from October 2016 to January 20 November 201718 to May 2019, and August 2019 to January 2020, respectively. There was a long time interval between the measurement of some plots and the LiDAR data acquisition in the eastern region (over 70% of plots had an interval of over 90 d). In order to characterize the forest structure states during the laser-scanning process, a set of stand growth models was established based on the reviewed field plot data from the Continuous Forest Inventory of Guangxi (CFI-GX) from 2005 to 2015, and the DBH, H, BA, and VOL of the field plots were adjusted via these stand growth models and tree growth rhythm tables (Huang et al. 1993, 1996). A summary of the field data based on forest type is presented in Table 1. For the same forest type, there were differences in the measured forest attributes of the plots in three regions, as shown by the great differences in both the mean and coefficient of variation.

Table 1. Summary statistics of the field plot data. CV is the coefficient of variation.

				DBH		H		BA		VOL
Region	Forest type	Sample size	Stem density (stem ha^−1)	Mean (cm)	CV (%)	Mean (m)	CV (%)	Mean (m^2ha^−1)	CV (%)	Mean (m^3ha^−1)	CV (%)
Eastern	Fir	61	1233–5567	12.7	21.05	10.65	27.57	32.04	29.49	186.82	50.85
	Pine	81	417–3917	18.0	30.47	12.56	26.58	25.07	34.67	159.93	49.04
	Eucalyptus	83	667–3117	10.6	20.53	14.95	19.10	14.56	40.35	114.03	50.31
	Broad-leaved	52	333–4050	13.5	35.54	11.38	29.82	22.78	36.55	141.99	56.77
Western	Fir	77	900–6883	12.7	27.10	11.70	24.81	34.58	28.81	227.51	44.06
	Pine	82	467–3967	20.1	27.04	15.89	22.19	32.32	24.07	245.44	35.10
	Eucalyptus	79	667–3167	11.8	24.00	17.04	22.31	18.08	30.71	154.49	27.81
	Broad-leaved	105	250–3067	13.9	34.52	10.94	23.55	19.74	35.48	114.12	50.01
Nanning	Fir	84	683–6650	11.7	31.12	10.47	26.96	30.16	24.93	173.67	36.16
	Pine	97	350–3667	19.0	27.72	13.05	26.56	26.66	32.03	173.67	47.79
	Eucalyptus	107	517–3350	11.4	19.43	16.28	16.64	18.29	31.29	150.12	38.74
	Broad-leaved	95	233–4800	13.6	40.28	9.58	26.25	17.83	47.92	93.10	63.24

2.2.2. LiDAR data

The LiDAR data were collected in three separate acquisitions in the Nanning, eastern, and western regions between October 2016 and April 20 October 201718 and October 2019, and August 2019 and January 2020, respectively. The Riegl VQ-1560 and Riegl VQ-1560i laser-scanning systems (Riegl, Austria) were applied with the same standards to collect LiDAR data in all three regions. The LiDAR survey flight and sensor parameters are shown in Table 2. The point clouds were georeferenced to a projection system of the China Geodetic Coordinate System 2000 (CGCS 2000). The mean square error of the laser point cloud height was smaller than 0.15 m.

Table 2. Summary of LiDAR survey flight and sensor parameters.

Parameter	Value
Laser wavelength	1064 nm
Pulse repetition frequency	2×1000 kHz
Laser beam divergence	0.5 mrad
Maximum Scan angle	30°
Return per pulse	5
Flight altitude above ground	2500 m
Average flying speed	200–240 km h^−1
Footprint diameter	60 cm
Side overlaps	21%-25%
Average ground point density + stdev.	5.54±2.14 pts·m^−2

Before being delivered to the end users, the raw data was preprocessed by the LiDAR data provider. The TerraScan software (TerraSolid Ltd., Helsinki, Finland) was used to label the point clouds as ground or non-ground returns using the adaptive triangulation network (TIN) filter algorithm. Utilizing ground returns, a digital terrain model (DTM) with a grid cell size of 2 m was finally generated. We removed the influence of topography and obtained DTM-normalized LiDAR point clouds based on the DTM. For the LiDAR height metrics and vertical structure metric, a 2.0-m and a 0.0-m height cutoff were employed in the computation of all features, respectively. All point clouds were used for density metric extraction. To extract LiDAR-derived variables, we used the package version 3.8 of Python (Python Software Foundation)

Height metrics are the significant variables in estimation models of forest attributes (Bouvier et al. 2015). The mean height of the point cloud height distribution (Hmean) is one of the most important variables (Asner et al. 2012; Lefsky, Cohen, and Spies 2001). Hmax, obtained from point clouds or the canopy height model (CHM), is an essential variable characterizing canopy height distribution and can be used to estimate forest attributes (Knapp et al. 2020). However, some studies suggested that the variation in Hmax extracted from laser point clouds was extremely unstable and was unsuitable for forest attribute estimation (Gobakken and Næsset 2008). Therefore, in this study, the 95^th height percentile (hp95) was used instead of Hmax. Additionally, the standard deviation (Hstdev), and coefficient of variation (Hcv) of the point cloud height distribution, characterized the heterogeneity of canopy height distribution, were also critical variables (Bouvier et al. 2015; Silva et al. 2017). Stand density depicts the horizontal distribution of trees within a field plot and is an essential variable for describing stand structure. Among various LiDAR-derived density variables, canopy cover (CC), which characterizes the cover of tree stems, branches, and leaves, effectively correlates the LiDAR point cloud with the horizontal distribution of trees and is one of the most commonly used LiDAR density metrics. The 50^th and 75^th canopy density percentiles (dp50 and dp75), which characterize foliage distribution in the middle and upper parts of the canopy, are also commonly used as density metrics (Næsset 2014). As subdominant trees and the understory exist under the canopy, LiDAR height and density metrics alone are not sufficient for accurately characterizing the complex canopy structure. In contrast, the leaf area density (LAD) profiles and vertical foliage profiles (VFPs) contain information on the vertical structure of the stand, providing information on canopy vertical heterogeneity and understory vegetation (Hopkinson et al. 2013; Lovell et al. 2003). In this paper, we extracted the canopy LAD-related metrics, including the mean (LADmean), standard deviation (LADstdev), and coefficient of variation (LADcv), following the methods described by Bouvier et al. (2015). In addition, we calculated the VFP-related variables, including the mean (VFPmean), standard deviation (VFPstdev), and coefficient of variation (VFPcv) following the methods described by Knapp et al. (2020).

2.3. Methods for developing the generalized model

TIn this study, a region-generalized model is defined as one that estimates a specific forest attribute of a forest type across the study sites with an accuracy that meets the given criteria. Similarly, a type-generalized model is defined as one that estimates a specific forest attribute across the forest types and study sites with an accuracy that meets the given criteria. The generalized models in this study involve four classes of models: region-generalized formulation, region-generalized model, type-generalized formulation, and type-generalized model. The difference between the generalized formulation and the generalized model is that the regression parameters of the latter have been determined, while those of the former have not.

The development of the generalized model involves a series of methods and procedures, as shown in Figure 2. This framework comprises the following steps: (1) the construction, calibration, and validation of the model. Thirteen LiDAR variables were used to construct model formulations using a rule-based exhaustive combination method, and 86 formulations were archived for each forest attribute; (2) after calibrating and validating all 86 models using datasets from each region, the best local model with the smallest rRMSE and the largest R² for a forest inventory attribute of a region and forest type was selected; (3) establishing the accuracy criteria for the generalized models based on the rRMSE and R² of the best local models; (4) using two methods to select the candidates for the best region- or type-generalized formulations; (5) calculating the differences of the rRMSE and R² between the candidates for the best region- or type-generalized formulations and the best local models, and then determining the best region- or type-generalized formulation; (6) using mixed data to calibrate and validate the best generalized formulation, and then determining the best generalized model based on the accuracy criteria of the generalized models. The detailed methods for each step will be described in the following sections.

Figure 2. A flowchart depicting a brief description of the procedures of determining the best region- or type-generalized model for a specific attribute. Note: a) While determining the region-generalized formulation, there are 86 models for each region; while determining the type-generalized formulation, there are four forest types in each region, so there are 4×86 models. b) While determining the region-generalized formulation, there is one best local model for each region; while determining the type-generalized formulation, there are four best local models for each region. c) for the mean method, while determining the region-generalized formulation, average the rRMSE and R² of the models in all regions (three models); while determining the type-generalized formulation, average the rRMSE and R² of the models in all regions and all forest types (12 models). d) for the mean method, while determining the region-generalized formulation, calculate the differences in the rRMSE and R² between the candidates for the best region-generalized models and the three best local models in three regions; while determining the type-generalized formulation, calculate the differences in the rRMSE and R² between the candidates for the best type-generalized models and the 12 best local models in three regions and four forest types. e) While determining the best region-generalized formulation, the mixed data consist of data from three regions of the same forest type; while determining the best type-generalized formulation, the mixed data consist of data from four forest types in three regions.

2.3.1. Model construction

In regional forest resource inventories and monitoring, the VOL is calculated using the following allometric equation:

(1) $VOL=a_{0}B{A}^{a_1}{\bar{H}}^{a_2}+\epsilon$(1)

where $VOL$ is the stand volume per hectare (m³ ha⁻¹); $BA$ is the stand basal area per hectare (m² ha⁻¹); $\bar{H}$ is the mean stand height (m); $a_0$, $a_1$, and $a_2$ are the regression model parameters; and $\epsilon$ is the error term. According to a previous study (Li and Li 2021), Eq. (1) can be written as follows when using LiDAR variables for estimating forest attributes.

(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2)

where $\hat{y}$ is the forest attribute, which can be VOL, BA, DBH, or H, etc.; P, H, and S are the density variable, height variable, and vertical structure variable extracted from LiDAR data, respectively, and they can be regarded as three groups of LiDAR variables, each group including 1–6 variables.

The 13 above-mentioned LiDAR-derived metrics could be categorized into three groups: the height metric (Ph), the density metric (Pd), and the vertical structure metric group (Pv) (Table 3), each group was supposed to capture a different structural aspect of the forest canopy. Based on the forest mensuration and ecological significance of LiDAR variables, and their correlations with the target variables, we specified that Hmean and hp95 were the primary height variables, CC was the primary density variable, and the remaining height and density variables were the secondary variables. The six vertical structure variables were not divided into primary and secondary variables.

Table 3. List of LiDAR-derived metrics used for establishing the prediction model.

Acronym	Explanation of metric	Structural aspect	Predictor (Px)
VOL	Stand volume (m^3 ha^−1)	Target variable	-
BA	Basal area (m^2 ha^−1)	Target variable	-
Hmean	Mean height of point clouds (m)	Canopy height	Phm
hp95	95^th height percentile	Canopy height	Phm
Hstdev	Standard deviation of point height distribution (m)	Canopy height	Ph
Hcv	Coefficient of variation of point height distribution	Canopy height	Ph
CC	Canopy cover	Canopy density	Pdm
dp50	50^th density percentiles	Canopy density	Pd
dp75	75^th density percentiles	Canopy density	Pd
LADmean	Mean of vertical leaf area density (LAD) profile	Vertical heterogeneity	Pv
LADstdev	Standard deviation of vertical LAD profile	Vertical heterogeneity	Pv
LADcv	Coefficient of variation of vertical LAD profile	Vertical heterogeneity	Pv
VFPmean	Mean of vertical foliage profile (VFP)	Vertical heterogeneity	Pv
VFPstdev	Standard deviation of VFP	Vertical heterogeneity	Pv
VFPcv	Coefficient of variation of VFP	Vertical heterogeneity	Pv

To explore the optimal combinations of variables to estimate the stand volume and basal area based on Equation (2)(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2) , we chose 1–2 variables in each variable group and added them to Equation (2)(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2) . The selection and combination of the variables are carried out according to the following rules (Li et al. 2022):

1. The variable combination scheme is as follows: 1–2 height variables+1–2 density variables+0–1 vertical structure variable.

2. The model must comprise at least a primary height variable (Hmean or hp95) and a primary density variable (CC). When two height variables are chosen, only a primary height variable and a secondary height variable can be chosen. When two density variables are chosen, a secondary density variable is chosen next to the CC.

3. Each variable in the group of vertical structure variables can appear in the model separately.

4. When the model is composed of two variables, Equation (2)(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2) consists of a primary height variable and a primary density variable. When the model comprises three variables, Equation (2)(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2) consists of a primary height variable, a primary density variable, and a vertical structure variable. When the model comprises four variables, Equation (2)(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2) consists of a primary height variable, a secondary height variable, a primary density variable, and a vertical structure variable. When the model is composed of five variables, Equation (2)(2) $\hat{y}=a_0{P}^{a_1}{H}^{a_2}{S}^{a_3}+\epsilon$(2) consists of a primary height variable, a secondary height variable, a primary density variable, a secondary density variable, and a vertical structure variable.

An exhaustive combination approach was utilized to combine the variables. Finally, 86 formulations consisting of 2–5 variables were obtained, as shown in Table A1.

2.3.2. Determining the best local model

Calibration and validation of the models were performed separately for the forest attributes, forest types, and study areas. Each model underwent 50 iterations of calibration and validation. The Gauss-Newton algorithm was used to estimate the model parameters using Python software. For each iteration, 80% of the randomly selected samples from the dataset were used as the calibration dataset, and the remaining 20% were used as the validation dataset. The evaluation statistics for the model include the relative root mean square error (rRMSE) and coefficient of determination (R²). After all 50 iterations, the mean rRMSE and the mean R² were calculated for each model. Repeating the above procedure, all 86 models were calibrated and validated for the forest attribute of a forest type in a study area. The best local model was the one that had the smallest mean rRMSE and the largest mean R² of the validation datasets among all 86 modes. There were 24 best local models for estimating two stand attributes (VOL and BA) of four forest types (Chinese fir, pine, eucalyptus, and broad-leaved forests) in three regions (the eastern, western, and Nanning regions).

2.3.3. Determining the best region-generalized formulation and model

The best region-generalized formulation or model is the formulation or model that has the smallest rRMSE and the largest R² among the region-generalized formulations or models. We proposed two methods: the mean method and the mixed data method, to determine the best region-generalized formulation and model. The difference between the two methods was the determination of the candidate for the best region-generalized formulation, and the remaining methods and steps were the same (Figure 2). The detailed methods and procedures are as follows.

1) Selection of the candidates for the best region-generalized formulations.

1. The mean method: For a specific attribute of a specific forest type in three study sites, average the rRMSE and the R² of three models achieved in Section 2.3.2, respectively. The formulation with the smallest mean of rRMSE and the largest mean of R² among all 86 models was the candidate for the best region-generalized formulation.

2. The mixed data method: The data from three regions of a specific forest type were mixed to calibrate and validate all 86 formulations using a similar method as in Section 2.3.2. The formulation with the smallest rRMSE and the largest R² among all 86 models was the candidate for the best region-generalized formulation.

2) Determining the best region-generalized formulation. Differences ($\mathrm{\Delta }rRMSE=\frac{rRMS{E}_g-rRMS{E}_l}{rRMS{E}_l}\times 100$ and $\mathrm{\Delta }{R}^2=\frac{R_g^2-R_l^2}{R_l^2}\times 100$) between the rRMSE and R² of the candidate for the best region-generalized formulation ($rRMS{E}_g$ and $R_g^2$) and the rRMSE and R² of the best local model ($rRMS{E}_l$ and $R_l^2$) in three regions were calculated. If the following accuracy criteria are satisfied, the candidate for the best regional-generalized formulation is regarded as the best region-generalized formulation: (a) the differences in rRMSE in all three regions are smaller than 10%, and the differences in R² in all three regions are larger than −10%; and (b) the average difference in rRMSE is smaller than 5%, and the difference in R² is larger than −5% (Table 4). Otherwise, the candidate formulation is not the best regional-generalized formulation. Thus, there was no best region-generalized model.

1. Synthesis selection: If both selection methods yielded a candidate for the best region-generalized formulation, the one with the smallest average difference in rRMSE was the best region-generalized formulation.

2. Determining the best region-generalized model. Field plots from three regions of a specific forest type were mixed, and 50 repeats of model calibration and validation were performed for the best region-generalized formulation. Given a large number of field plots after data mixing, for each repetition, 70% of randomly selected samples from the mixed dataset were utilized to calibrate the model, and the remaining 30% were used to validate the model. The mean rRMSE and the mean R² of the validation datasets in each region were calculated separately based on the validation data source. The model was confirmed as the best region-generalized model if the differences between the means rRMSE (and mean R²) of the models in three regions and the rRMSE (and mean R²) of the best local model met the accuracy criteria specified in Table 4. Otherwise, there was no best region-generalized model.

Table 4. The accuracy criteria for the generalizing formulation and model based on the difference in rRMSE and R² between the candidates for the best region- or type-generalized formulations or models and the best local model.

Generalized formulation or model	Max. difference (%) of all regions		Average difference (%)
	ΔrRMSE	ΔR^2	ΔrRMSE	ΔR^2
Region-generalized formulation	≤10	≥-10	≤5	≥-5
Region-generalized model	≤15	≥-15	≤7.5	≥-7.5
Type-generalized formulation	≤15	≥-15	≤7.5	≥-7.5
Type-generalized model	≤20	≥-10	≤10	≥-10

2.3.4. Determining the best type-generalized formulation and model

The best type-generalized formulation or model is the formulation or model that has the smallest rRMSE and the largest R² for estimating among the type-generalized formulation or model. The method and procedure for determining the best type-generalized formulation and model are similar to those described in Section 2.3.3 for determining the best region-generalized formulation and model. The main difference is that the former involves not only the study area but also the forest type, and the amount of dataset is much larger than that of the latter, so the methods used to determine the candidates for the best type-generalized formulation are different.

1. The mean method: For a specific attribute of all four forest types at all three study sites, 12 models were achieved in the above procedures (Section 2.3.2) for each of the 86 formulations. Averaging the rRMSE and the R² of these 12 models, we could thus determine the formulation with the smallest mean of rRMSE and the largest mean of R² as the candidate for the best type-generalized formulation.

2. The mixed data method: The data from three regions and four forest types were mixed and used to calibrate and validate all 86 formulations (50 repeats). The candidate for the best type-generalized formulation was the one with the smallest mean rRMSE and the largest mean R².

If the differences in the rRMSE (and in the R²) between the candidate for the best type-generalized formulation and the best local models of three regions of four forest types were all smaller than 15% (and larger than −15%) and the average difference in the rRMSE (and in the R²) was smaller than 7.5% (and larger than −7.5%) (as specific in Table 4), the candidate formulation was identified as the best type-generalized formulation. The best type-generalized model was further determined using a procedure similar to that described above.

3. Results

3.1. The best local model for stand attribute prediction

The performances of the 86 VOL estimation formulations varied considerably across the three study sites and the four forest types. For example, in the fir forests, formulations 8, 3, and 21 performed best in the eastern region, while in the western and Nanning regions, the best-performing formulations were formulations 3, 4, and 21, and formulations 27, 28, and 33, respectively. There were seven formulations among the nine best-performing models in the three regions, but none of them were present in all three regions simultaneously. In the pine forests, the three best-performing models were completely different in the three regions. In the eucalyptus and broad-leaved forests, they had seven best-performing formulations in all three regions, but none of them in all three regions simultaneously (Table 5). The above results indicated that the performances of all formulations varied widely across regions, even if the forest types were the same.

Table 5. Model validation statistics for the best-performing models for the VOL and BA predictions in the three regions of the four forest types.

			East			West			Nanning
Forest type	Attribute	Order	No. of model	rRMSE (%)	R^2	No. of model	rRMSE (%)	R^2	No. of model	rRMSE (%)	R^2
Fir	VOL	1	8	18.52	0.837	3	21.56	0.727	27	19.05	0.684
		2	3	18.82	0.832	4	21.97	0.717	28	19.14	0.681
		3	21	18.90	0.831	21	22.05	0.715	33	19.16	0.680
	BA	1	21	18.12	0.582	21	19.22	0.505	79	16.24	0.550
		2	15	18.19	0.579	15	19.24	0.504	78	16.25	0.549
		3	22	18.20	0.578	3	19.28	0.501	43	16.29	0.547
Pine	VOL	1	81	16.57	0.874	67	16.59	0.755	61	23.30	0.709
		2	86	16.59	0.874	79	16.71	0.751	49	23.31	0.709
		3	85	16.60	0.873	63	16.85	0.746	55	23.35	0.710
	BA	1	85	16.81	0.750	67	16.41	0.500	7	19.68	0.560
		2	86	16.83	0.751	79	16.52	0.494	55	19.93	0.548
		3	81	16.86	0.749	68	16.63	0.487	43	19.95	0.547
Eucalyptus	VOL	1	81	24.67	0.746	69	17.70	0.770	85	17.84	0.767
		2	86	24.67	0.746	81	17.73	0.769	73	17.91	0.765
		3	69	24.68	0.746	73	17.75	0.769	49	18.11	0.761
	BA	1	81	22.85	0.660	74	16.81	0.618	73	16.44	0.689
		2	82	22.87	0.660	86	16.84	0.616	85	16.45	0.689
		3	69	22.88	0.659	72	16.97	0.612	72	16.60	0.684
Broadleaf	VOL	1	7	34.94	0.575	7	31.16	0.588	6	39.79	0.574
		2	8	34.99	0.571	8	31.30	0.584	26	40.66	0.552
		3	5	35.32	0.562	3	31.40	0.581	20	40.72	0.550
	BA	1	8	27.02	0.393	8	26.12	0.419	6	31.80	0.531
		2	7	27.25	0.382	3	26.20	0.415	62	32.28	0.515
		3	6	27.36	0.378	5	26.21	0.414	50	32.35	0.514

Similar to VOL estimation, the performances of all 86 formulations for BA estimation varied greatly across the regions and forest types, with none ranking among the three best models in all three regions and all forest types simultaneously (Table 5). Among the nine best-performing formulations in three regions, there were seven, nine, eight, and seven formulations in the fir, pine, eucalyptus, and broad-leaved forests, respectively.

Further analysis revealed that the difference between the largest rRMSE and the smallest rRMSE did not exceed 10% among the 20 best-performing VOL estimation models in regions and forest types. The differences in R² were similar to those in rRMSE (except for R² of the VOL model for the broad-leaved forests in the eastern region), indicating that the performances of the 20 best-performing VOL estimation models in three regions and four forest types were analogous to each other, although their formulations differed considerably. The performances of the 20 best-performing models for BA estimation in three regions of four forest types were comparable to the performances of the 20 best-performing models for VOL estimation. The best local models for VOL and BA estimation in three study areas and four forest types are shown in Tables 6 and 7, respectively.

Table 6. The best local models for VOL estimation of four forest types in three regions and the best region-generalized model for VOL estimation of four forest types and their goodness-of-fit and validation statistics.

			Parameter estimates														Fit static		Validation static
Forest type	Study site	No. of model	a0	Hmean	hp95	Hstdev	Hcv	CC	dp50	dp75	LADmean	LADstdev	LADcv	VFPmean	VFPstdev	VFPcv	R^2	rRMSE (%)	R^2	rRMSE (%)
Fir	East	8	10.6383	1.6153				0.9414								−0.3709	0.901	15.98	0.837	18.52
	West	3	9.1707	1.3638				0.4880			0.1034						0.810	18.71	0.727	21.56
	Nanning	27	1.9457		2.0135	−0.7331		1.2081			0.08881						0.772	17.10	0.684	19.05
	Region- generalizing	21	9.2016	1.1881			−0.2487	1.3008			0.08589						0.783	21.19	0.735	22.65
Pine	East	81	5.9902		1.3105		−0.1889	0.4012		0.2219	−0.02103						0.915	14.32	0.874	16.57
	West	67	4.0108		1.7129	−0.3182		0.0229	0.5847						−0.1701		0.825	14.58	0.755	16.59
	Nanning	61	10.9203	1.1902			0.0878	0.1686		0.1604					0.1490		0.780	19.75	0.709	23.30
	Region- generalizing	49	6.4800	1.3281		0.01814		0.1886		0.09635					0.1425		0.841	18.70	0.812	19.88
Eucalyptus	East	81	1.9994		1.5900		−0.0667	0.1036		0.4001	0.0139						0.827	20.95	0.746	24.67
	West	69	2.6458		1.4860	−0.03930		0.2495		0.2788	0.0446						0.874	14.46	0.770	17.70
	Nanning	85	13.0395		1.1879		−0.0850	0.3714		0.3623					−0.3905		0.810	16.73	0.767	17.84
	Region- generalizing	85	5.0809		1.3428		−0.07944	0.3134		0.3250					−0.1492		0.829	18.13	0.801	19.16
Broadleaf	East	7	1.5380	1.5468				−0.2860							0.3532		0.725	29.41	0.575	34.94
	West	7	3.7393	1.2614				0.6556							0.2668		0.661	29.09	0.588	31.16
	Nanning	6	0.7975	2.1910				0.3678						1.1365			0.695	33.92	0.574	39.79
	Region- generalizing	6	4.5387	1.3929				0.2621						0.1601			0.647	34.55	0.596	36.22

Hmean, hp95, … , VFPcv are the estimates of the model variables Hmean, hp95, … , VFPcv, respectively; Region generalized is the best forest type-specific-region-generalized model.

Table 7. The best local models of BA estimation of four forest types in three study sites and the best region-generalized model of BA estimation of four forest types and their goodness-of-fit and validation statistics.

			Parameter estimates														Fit static		Validation static
Forest type	Study site	No. of model	a0	Hmean	hp95	Hstdev	Hcv	CC	dp50	dp75	LADmean	LADstdev	LADcv	VFPmean	VFPstdev	VFPcv	R^2	rRMSE (%)	R^2	rRMSE (%)
Fir	East	21	7.0246	0.6049			−0.1420	0.8668			0.06619						0.709	15.63	0.582	18.12
	West	21	7.1746	0.5548			−0.1876	0.5861			0.08123						0.656	16.66	0.505	19.22
	Nanning	79	1.4097		0.6468		−0.7607	1.0852	−0.2836						0.1591		0.662	14.37	0.550	16.24
	Region- generalizing	39	5.6103	0.8950		−0.3141		1.0094	−0.09915		0.05261						0.622	17.39	0.564	18.44
Pine	East	85	3.0794		0.7482		−0.1996	0.7311		0.1463					0.1318		0.819	14.71	0.750	16.81
	West	67	4.6217		0.9335	−0.2664		0.1434	0.4678						−0.1600		0.620	14.63	0.500	16.41
	Nanning	7	1.4333	0.9636				0.1488							0.3590		0.668	18.29	0.560	19.68
	Region- generalizing	55	5.5222	0.6716			0.02978	0.09258	0.3509						0.1093		0.701	17.27	0.657	18.21
Eucalyptus	East	81	0.8547		1.1372		−0.07231	0.2911		0.2879	0.01208						0.759	19.75	0.660	22.85
	West	74	1.8442		1.0257	−0.07322		0.3473		0.1747						−0.1999	0.758	14.42	0.618	16.81
	Nanning	73	4.3178		0.8585	−0.0622		0.3365		0.2963					−0.3315		0.752	15.20	0.689	16.44
	Region- generalizing	73	2.3821		0.9392	−0.07953		0.3657		0.2530					−0.1493		0.747	17.19	0.704	18.22
Broadleaf	East	8	2.8256	0.6586				−0.2044								0.2286	0.564	22.99	0.393	27.02
	West	8	3.9090	0.6088				0.7449								0.1542	0.515	24.31	0.419	26.12
	Nanning	6	0.6663	1.5649				0.6187						0.8739			0.629	29.00	0.531	31.80
	Region- generalizing	6	3.3333	0.8105				0.5059						0.1404			0.517	28.07	0.458	29.52

3.2. The best region-generalized formulation and model

After averaging the rRMSE and R² of all 86 models for a specific attribute of a specific forest type in three regions, it was found that the performances of all 86 formulations differed greatly from their performances in all three regions. For example, in the VOL estimation of the fir forests, the three best-performing formulations selected by the mean method were 1, 19, and 13, respectively, and model formulation 13 was not among the three best-performing models in all three regions (Table 5). Similarly, the performances of all 86 formulations obtained by the mixed data method differed distinctly from their performances in three regions, e.g. the three best-performing formulations in the VOL estimation of the fir forest were formulations 19, 20, and 55, and formulations 20 and 55 were not among the three best-performing models in all three regions. The sixteen candidates for the best region-generalized formulations obtained by two selection methods are shown in Table 8. The candidates for the best region-generalized formulations for the VOL and BA estimations obtained by the two selection methods were different, except for the BA estimation in the eucalyptus forests.

Table 8. The number of the candidates for the best region-generalized formulations for the VOL and BA estimations of four forest types selected by the mean method and mixed data method.

		Fir			Pine			Eucalyptus			Broadleaf
Selection method	Attribute	No. of model	rRMSE (%)	R^2	no. of model	rRMSE (%)	R^2	No. of model	rRMSE (%)	R^2	No. of model	rRMSE (%)	R^2
Mean	VOL	3	18.82	0.832	49	19.55	0.759	85	20.17	0.758	6	35.82	0.563
	BA	39	16.24	0.526	55	18.19	0.580	73	18.83	0.651	6	28.52	0.438
Mixed data	VOL	21	22.23	0.745	51	19.69	0.816	86	19.14	0.801	7	36.08	0.600
	BA	81	18.38	0.567	60	10.33	0.848	73	18.22	0.704	7	29.40	0.462

Model formulation 3 was the candidate for the best region-generalized formulation for VOL estimation of the fir forests selected by the mean method. In the eastern region, the rRMSE and R² of this formulation were 18.82% and 0.832, respectively. In the eastern region, the rRMSE and R² of the best local model (model (8)) for the VOL estimation of the fir forests were 18.52% and 0.837, respectively (Table 5), and the differences in rRMSE and R² between the two models were 1.62% and −0.60%, respectively. Using the above-mentioned approach, we obtained the differences in rRMSE and R² between the sixteen candidates for the best region-generalized formulations for the VOL and BA estimations of four forest types and the best local models for VOL and BA estimations, as shown in Table 9.

Table 9. The difference in rRMSE (δrrmse) and R² (ΔR²) between the sixteen candidates for the best region-generalized formulations for the VOL and BA estimations obtained by two selection methods and the best local models for four forest types in three study sites.

				East		West		Nanning		Average
Selection method	Forest type	Attribute	No. of model	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)
Mean	Fir	VOL	3	1.62	−0.60	0.00	0.00	8.94	−8.62	3.52	−3.07
		BA	39	2.59	−3.68	2.07	−3.94	1.75	−2.78	2.14	−3.47
	Pine	VOL	49	2.83	−0.83	10.26	−6.83	0.06	−0.02	4.39	−2.56
		BA	55	5.38	−3.52	3.09	−6.22	1.29	−2.13	3.25	−3.96
	Eucalyptus	VOL	85	0.27	−0.19	1.35	−0.95	0.00	0.00	0.54	−0.38
		BA	73	0.55	−0.55	1.52	−1.99	0.00	0.00	0.69	−0.84
	Broadleaf	VOL	6	3.60	−7.05	1.01	−1.37	0.00	0.00	1.54	−2.81
		BA	6	1.27	−3.85	1.05	−3.22	0.00	0.00	0.77	−2.36
Mixed data	Fir	VOL	21	2.06	−0.64	2.26	−1.64	6.10	−5.79	3.47	−2.69
		BA	81	3.87	−5.77	3.37	−6.42	6.48	−9.66	4.57	−7.28
	Pine	VOL	51	17.20	−5.43	6.00	−3.96	2.15	−1.36	8.45	−3.58
		BA	60	4.52	−3.01	7.41	−15.34	7.66	−12.90	6.53	−10.42
	Eucalyptus	VOL	86	0.01	−0.03	0.45	−0.22	4.42	−2.72	1.63	−0.99
		BA	73	0.55	−0.55	1.52	−1.99	0.00	0.00	0.69	−0.84
	Broadleaf	VOL	7	0.00	0.00	0.00	0.00	7.11	−10.75	2.37	−3.58
		BA	7	0.84	−2.75	0.52	−1.46	8.18	−14.93	3.18	−6.38

Table 9 shows that: (1) the differences in rRMSE between the eight candidates for the best region-generalized formulations of VOL and BA estimations selected by the mean method and the best local models in three regions were all smaller than 10%, except for the pine forests in the western region, for which the difference in the VOL model was 10.26%; (2) the differences in R² between the eight candidates for the best region-generalized formulations of VOL and BA estimations selected by the mean method and the best local models in three regions were all larger than −10%; and (3) the averaged differences in rRMSE in the three regions were smaller than 5% for all four VOL and BA formulations selected by the mean method, and their average differences in R² were larger than −5%. The above result indicated that seven of the eight candidates for the best region-generalized formulation selected by the mean method met the criteria for the region-generalized formulation. The only exception was the VOL estimation model of the pine forests. The performances of the eight candidates for the best region-generalized formulations of VOL and BA estimation selected by the mixed data method were worse than those selected by the mean method, except for the VOL estimation in the fir forests.

Comprehensively considering the performances of all candidate formulations and the fact that the candidate VOL model (model 49) of the pine forest selected by the mean method slightly exceeded the criteria only in the western region, the eight best formulations for VOL and BA estimation were obtained based on two selection approaches, as shown in Tables 6 and 7. The best models for VOL and BA estimations for broad-leaved forests were the same; therefore, there were seven generalized formulations for the VOL and BA estimations in four forest types.

The mixed data from all three regions of the same forest type was used to calibrate and validate the best region-generalized formulations for each forest attribute. The parameter estimates and the goodness-of-fit statistics of the eight best region-generalized models are shown in Tables 6 and 7. The models were referred to as the pending best region-generalized models because their generalizability had not been evaluated.

When calibrating and validating the pending models, the validation datasets were divided into three groups (regions) according to their sources (regions), and the rRMSE and R² of the validation datasets in each region were calculated. The results showed that the variation in rRMSE for the VOL and BA estimation models of the fir, pine, and eucalyptus forests ranged from 16.52% to 23.15%, and the rRMSE of the VOL and BA estimation models of the broad-leaved forest was smaller than 35% (except for the VOL model in the Nanning region) (Table 10).

Table 10. The rRMSE and R² of the eight pending best-region-generalized models in three regions of four forest types.

			East region		West region		Nanning region
Forest type	Attribute	No. of model	rRMSE(%)	R^2	rRMSE(%)	R^2	rRMSE(%)	R^2
Fir	VOL	21	18.18	0.838	22.87	0.682	21.84	0.587
	BA	39	17.66	0.641	19.40	0.542	17.94	0.483
Pine	VOL	49	17.45	0.870	18.60	0.715	22.58	0.775
	BA	55	16.91	0.757	16.52	0.521	20.39	0.596
Eucalyptus	VOL	85	23.15	0.786	17.69	0.807	18.14	0.778
	BA	73	21.72	0.705	17.65	0.636	16.64	0.706
Broadleaf	VOL	6	34.67	0.610	33.41	0.545	41.27	0.559
	BA	6	26.99	0.427	27.01	0.412	34.63	0.461

From Tables 10 and 5, we obtained the differences in rRMSE and R² between the eight pending best region-generalized models and the best local model of each forest type in all regions, as shown in Table 11.

Table 11. The differences in rRMSE (δrrmse) and R² (ΔR²) between the eight pending best region-generalized models and the best local models for the stand volume and basal area estimations of four forest types in three regions.

			East		West		Nanning		Average
Forest type	Attribute	No. of model	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)
Fir	VOL	19	0.11	−1.83	−6.23	6.07	−14.18	14.65	−6.77	6.29
	BA	37	10.18	−2.55	7.42	0.93	−12.26	10.48	1.78	2.95
Pine	VOL	47	−0.44	5.32	−5.35	12.11	9.36	−3.10	1.19	4.77
	BA	53	0.97	0.57	4.21	0.66	6.39	3.60	3.85	1.61
Eucalyptus	VOL	83	5.33	−6.18	4.84	−0.05	1.46	1.68	3.88	−1.52
	BA	71	6.89	−4.93	2.84	5.01	2.54	1.21	4.09	0.43
Broadleaf	VOL	4	6.02	−0.78	−7.29	7.22	−2.64	3.72	−1.30	3.39
	BA	4	8.73	−0.11	−1.63	3.39	−13.12	8.89	−2.01	4.06

All eight pending region-generalized models met the criteria for region-generalized models specified in Table 4, indicating that these models are generalizable. The region-generalized models shown in Tables 6 and 7 could be identified as the region-generalized models for the VOL and BA estimations in three regions. Comparing the model estimates with the field measurements, we found that the model estimates were closer to the 1:1 line (Figure 3), demonstrating that the region-generalized models provided excellent predictions.

Figure 3. Scatterplots of survey VOL versus predicted VOL (a, c, e, g) and survey BA versus predicted BA (b, d, f, h) of the best region-generalized model of validation data for four forest types. The black, blue, and red dots the indicate validation data from the Eastern, Western, and Nanning regions, respectively.

In addition, Figure 3 shows that the predictions of the two forest attributes of the eucalyptus forest, which had a simple stand structure, were slightly better than those for other forest types. Since the study area was too large to demonstrate the differences in forest inventory attribute estimates, we selected a county for mapping the forest attributes, and the results showed that the estimates of the region-generalized models are very close to those of the best local models (Figure 4).

Figure 4. Wall-to-wall maps of the stand volume and basal area generated by the models of the pine forest in a selected area of 4550 km²: (a) and (b) are the estimates of the VOL and BA predicted by the region-generalized models, respectively; (c) and (d) are the estimates of the VOL and BA predicted with the best locale models in the Western region, respectively; (e) is the geographic location of the selected area in the study area.

After performing a logarithmic transformation, eight linear models were built using the same variables as the eight region-generalized multiplicative power models, and the standardized coefficients of the models are shown in Table 12.

Table 12. Standardized coefficients of the linear regression (logarithmic conversion) of the best region-generalized model of VOL and BA estimations of four forest types.

Forest type	Attribute	Hmean	hp95	Hstdev	Hcv	CC	dp50	dp75	LADmean	VFPmean	VFPstdev
Fir	VOL	0.828			−0.108	0.187			0.080
	BA	1.011		−0.386		0.350	−0.134		0.096
Pine	VOL	0.820		−0.010		0.141		0.119			0.022
	BA	0.667			0.032	0.145	0.246				0.071
Eucalyptus	VOL		0.663		−0.130	0.104		0.256			−0.075
	BA		0.660	−0.159		0.134		0.298			−0.086
Broadleaf	VOL	0.874				0.219				0.163
	BA	0.716				0.362				0.175

Overall, the primary height variables (Hmean and hp95) contributed the most to the target variables, followed by the density variables (sometimes, the secondary density variables contributed more than the primary density variables), and the vertical structure variables and secondary height variables contributed slightly less.

3.3. The performances of the best type-generalized formulations

The determination of the best type-generalized formulations of a specific attribute involved 12 datasets from three study sites and four forest types. For each attribute, 12 models of each of the 86 formulations were obtained in Section 2.3.1. After averaging and sorting the rRMSE and R² of these 12 models, we determined that formulations 57 and 61 were the candidates for the best type-generalized formulations for VOL and BA estimations selected by the mean method. For the maxing data methods, formulations 52 and 67 were the candidates for the best type-generalized formulations for VOL and BA estimations. The differences in rRMSE and R² between these candidates for the best type-generalized formulations for VOL and BA estimation and the best local models in the three regions and four forest types are shown in Table 13.

Table 13. The differences in rRMSE (δrrmse) and R² (ΔR²) between the eight candidates for the best type-generalized formulation of the VOL and BA estimations selected by two selection methods and the best local models in three regions of four forest types.

				East		West		Nanning
Selection method	Attribute	Forest type	No. of model	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)	ΔrRMSE(%)	ΔR^2(%)
Mean	VOL	Fir	57	4.45	−1.66	4.03	−3.10	7.77	−7.24
		Pine	57	4.02	−1.12	9.82	−6.60	2.71	−2.23
		Eucalyptus	57	1.00	−0.82	6.16	3.77	1.32	−0.16
		Broadleaf	57	15.12	−25.22	2.23	−3.16	4.28	−7.83
	BA	Fir	61	3.53	−5.24	5.28	−10.33	7.28	−11.04
		Pine	61	1.24	−0.79	6.89	−14.14	4.27	−7.36
		Eucalyptus	61	0.34	−0.52	3.90	−4.98	2.12	−1.72
		Broadleaf	61	5.57	−18.71	2.16	−6.01	1.96	−3.57
Mixed data	VOL	Fir	52	2.06	−0.64	2.26	−1.64	6.10	−5.79
		Pine	52	17.20	−5.43	6.00	−3.96	2.15	−1.36
		Eucalyptus	52	0.01	−0.03	0.45	−0.22	4.42	−2.72
		Broadleaf	52	0.00	0.00	0.00	0.00	7.11	−10.75
	BA	Fir	67	3.87	−5.77	3.37	−6.42	6.48	−9.66
		Pine	67	4.52	−3.01	7.41	−15.34	7.66	−12.90
		Eucalyptus	67	0.55	−0.55	1.52	−1.99	0.00	0.00
		Broadleaf	67	0.84	−2.75	0.52	−1.46	8.18	−14.93

Calibrating models 57 and 61 using all the data from four forest types and three regions, the results were that, for the VOL estimation, the rRMSE and R² were 25.56% and 0.562, respectively; for the BA estimation, the rRMSE and R² were 28.78% and 0.489, respectively.

The differences in rRMSE and R² between the candidates for the best type-generalized formulations and the best local models were significantly larger than those of the best region-generalized formulations and models (Tables 9, 11, and 13). For both the VOL and BA predictive models, there were always cases in which some forest types had obvious increases in rRMSE or obvious decreases in R² in some regions. The largest differences in rRMSE between the two candidate VOL formulations and the best local VOL models were 15.12% and 17.20%, respectively. The largest differences (absolute value) in R² between the two candidate BA formulations and the best local models were 18.71% and 15.34%, respectively (Table 13). Since the increase in errors (rRMSE) of the two candidate formulations and the decrease in the explanatory power (R²) of the model variables were considerably larger than those of the best local model, the criteria for the best type-generalized formulations specified in Table 4 were not satisfied. Hence, there were no best type-generalized formulations, and thus the best type-generalized models were not found.

4. Discussion

A generalized model or model formulation would facilitate airborne LiDAR-based estimation of forest attributes across regions and forest types (or species), which would not only be beneficial for the consistent estimation of forest inventory attributes but also help to reduce the number of field plots and the cost of field measurements. The novelties of this study are the accuracy criteria to determine the region- and type-generalized formulations and models for estimating the stand volume and basal area in subtropical forests, and the procedures for determining the generalized formulations and models.

4.1. Accuracy criteria for determining the generalized formulations and models

Consistent with the findings of most studies, the accuracy of the best local models of the sample forest type differed among the three regions (Tables 6 and 7); we suppose the possible reason is that the difference in the 3D structure of the forest canopy and the variation of the target variables vary in different regions (Table 1). Therefore, the determination of the generalized model needs to take into account the differences in the structure of forest stands in different forest types and regions, and it is necessary to maintain relatively high accuracy in all regions and all forest types instead of just considering the overall accuracy of the generalizing model. As the study area increases, the accuracy of the generalizing model is generally lower than that of the best local model (Tables 6 and 7); hence, the generalization of the predictive model requires us to make a trade-off between the estimation accuracy and the consistency of the estimate across study sites and forest types, as shown in the study by Bouvier et al. (2015). Some scholars have worked on developing generalized models for airborne LiDAR-based forest attribute estimation to achieve consistent estimates. However, to the best of our knowledge, no accuracy criteria have been proposed for the generalized formulation and model. Since there are four classes of generalized models, we have established four accuracy criteria. These accuracy criteria were based on the rRMSE and R² of the best local model, which took into account the differences in the characteristics of different forest types in different regions. We anticipated that the resulting generalized formulations or models would be highly accurate across all study sites or forest types. Otherwise, the generalized formulations or models would be meaningless.

Model accuracy is affected by numerous factors, such as the extent of the study area, heterogeneity of the forest, sampling method, number and area of the plots, modeling method, and so on (Zolkos, Goetz, and Dubayah 2013; Coops et al. 2021), among which the number of sample plots (sample size) is an essential factor. Normally, model accuracy gradually increases as the sample sizes increases (Gobakken and Næsset 2008; da Silva et al. 2020). However, the sample size required for forest attribute estimation varies by forest type, the extent of the study site, sampling, and modeling approach (Næsset 2014). In this study, the sample sizes for all forest types in most regions were less than 100; therefore, whether the accuracy of the best local model reached stability or not, additional study is needed.

Our study showed that it is difficult to obtain the best type-generalized formulations and models for the VOL and BA estimations because the accuracy of the candidates for the best type-generalized formulations was significantly lower than that of the best local models in some regions, and did not met the accuracy criteria of the type-generalized formulation. This finding differed from those of some studies (Bouvier et al. 2015; Knapp et al. 2020). The major reason is that none of the previous studies proposed an accuracy criterion for the generalized model, but only determined its generalizability with the rRMSE and R². In our study, the rRMSE and R² of the pending type-generalized formulation for the VOL estimation were 25.56% and 0.562, respectively. Although the statistics were not too low, the rRMSE was significantly larger than that of the best local models for some forest types in some regions, and the opposite was true for the R². Therefore, we believed that these two pending models had poor generalization capabilities and limited practical value. From the perspective of forest mensuration and ecology, this may be attributed to the various biological and ecological characteristics of different tree species, they lead to the 3D structures of the canopy varying considerably among forest types, e.g. most eucalyptus forests are single-storied pure forests, most broad-leaved forests are multi-storied mixed forests, and the fir and pine forests consist of both single-storied pure forests and multi-storied mixed forests. In addition, since the stem- and stand-form factors vary with tree species, the allometric equations for calculating the stem and stand volumes also vary considerably. Therefore, a small set of LiDAR variables can not accurately characterize the diverse canopy structures of the different forest types. As a result, there is no best type-generalized formulation. One potential solution is to stratify the vertical structure of the forest canopy and then estimate the forest attributes separately by stratum (Zhou et al. 2022; Yu et al. 2022). Perhaps adding 1–2 LiDAR variables to the model is also an alternative solution.

4.2. Accuracy of generalized model

Only a few studies on airborne LiDAR-based forest attribute estimation in subtropical forests have been reported, and all of them have small study areas. In two studies conducted in planted forests in Guangxi Province (with an area of 5000 ha) and Jiangsu Province (1260 ha), China, the rRMSEs for estimating the stand volume using multiple linear regression models were 21.34% and 16.47%, respectively (Liu et al. 2021; Zhang, Cao, and She 2017). The results of these two studies are not comparable to this study due to the small extent and homogeneity of the forests. In some studies conducted in temperate forests larger than 5,000 km², the estimation accuracy of forest inventory attributes using parametric models varied considerably, with rRMSE ranging from 11.04% to 46.7% for the stand volume and 13.8% to 37.1% for the basal area (Hauglin et al. 2021; Hill, Buddenbaum, and Mandallaz 2018; Nilsson et al. 2017; Watt and Watt 2013; Treitz et al. 2012; Nord-Larsen and Schumacher 2012; Dalponte et al. 2011; Woods, Lim, and Treitz 2008). In 34 studies worldwide, the mean residual standard error (RSE) of discrete airborne LiDAR-based aboveground biomass estimation was 27% (Zolkos, Goetz, and Dubayah 2013). The comparability of model accuracy between different studies is poor due to the influence of various factors mentioned earlier. In this study, the rRMSEs of the regional-generalized models used for stand volume estimation were 22.65%, 19.88%, 19.16%, and 36.22% for the fir, pine, eucalyptus, and broad-leaved forests, respectively, although they were slightly lower than the best local models. Given the large extent of the study area and the heterogeneity of forest structures in different regions, we believe that the accuracies of the regional-generalized models in this study are achieved as desired.

4.3. Interpretability of the model variables

Numerous LiDAR-derived metrics can be extracted from LiDAR point clouds and CHM. All these metrics have explicit statistical meaning and characterize aspects of the 3D structure of the forest canopy. However, some of them, such as the 5% and 10% percentile heights, the 0.15–1.37 m echo ratio (Fekety, Falkowski, and Hudak 2015), and the 10% and 99% density percentiles (Xu, Manley, and Morgenroth 2018), have no clear meanings in forest mensuration or ecology. Although these metrics are closely correlated with some specific attributes for certain forest types in certain regions and can be used for forest attribute estimation, these models are poorly interpretable. The 13 LiDAR-derived metrics selected for this study have clear forest mensuration and ecological significance, and accurately depict the 3D structure of the forest canopy in terms of height, density, and vertical structure; thus, the constructed models had excellent interpretability. In addition, the correlation coefficients between these 13 variables are small, effectively reducing the multicollinearity of the independent variables in the model. There are various approaches to selecting model variables (Montealegre et al. 2016; Luo et al. 2018; van Ewijk et al. 2019; Silva et al. 2017; Sheridan et al. 2015; Zhang, Cao, and She 2017; Cao et al. 2019), all of which have advantages and limitations. In this study, the rule-based exhaustive combination of 13 LiDAR variables not only enabled the high accuracy of the resulting model but also made the model interpretable in terms of forest mensuration and ecology. Although the calibration and validation of 86 models were computationally intensive (this is no longer an issue today due to the powerful computing power), it had the advantage of achieving the best combination of LiDAR-derived variables suitable for different forest attribute estimations for different forest types (i.e. different canopy structures), and the model variables had excellent explanatory power.

The vertical structure variables are calculated in this paper in terms of gap fractions using the Beer-Lambert law. The clumping effect of leaves might increase the gap probability and thus lead to an estimate of the leaf area index (LAI) being smaller than the actual LAI (Hu et al. 2018). Therefore, it is necessary to investigate the physical models based on stochastic radiative transfer theory, such as the model of Detto et al. (2015), to extract the vertical structure information of forest stands more accurately.

The predictors of the best local models varied significantly across study sites and forest types. The major reason may be the differences in the 3D structure of the forest canopy among diverse forest types and diverse study sites. Therefore, future work will focus on classifying forest canopy structures using LiDAR data (Zhou et al. 2022) and then developing generalized formulations and models for different canopy structures (single-, double-, and multi-storied) (Yu et al. 2022), to further improve model generalization and estimation accuracy while maintaining consistent estimation of forest attributes. Furthermore, despite this study covering a whole province, the extent of the distribution of major tree species is broad regardless of the biogeographic region; e.g. the distribution of Cunninghamia lanceolata and Pinus massoniana cover approximately half of China. Therefore, the adaptability of the region-generalized model needs to be validated in a larger region than a province.

5. Conclusion

We obtained 86 formulations in this study for estimating the stand volume and basal area using a rule-based exhaustive combination of 13 LiDAR-derived metrics that characterized the 3D structure of the forest canopy and had explicit forest mensuration and ecological significance. Then, we obtained the 24 best local models by calibrating and validating these formulations using data from four forest types in the three study sites. We further proposed a set of accuracy criteria for determing the generalized formulations and models based on these models. Finally, by applying two selection methods based on the rRMSE and R2, we achieved the eight best region-generalized formulations and models. However, since the 3D structure of the canopy varies considerably among forest types, it was difficult to achieve type-generalized formulations and models. The accuracy criteria and procedures presented in this paper would be helpful in the development of generalized models for the consistent estimation of forest inventory attributes using airborne LiDAR data. However, the accuracy criteria still need to be further validated in a larger region than a province. Furthermore, additional stratification of forest types (or species) based on the vertical structure is required to obtain stratum-generalized formulations and models for forest attribute estimations that are consistent across forest types and study sites.

CRediT authorship contribution statement

Chungan Li: Conceptualization, Methodology, Calculation, Data collection and analysis, Writing – original draft and – Review and editing, Project administration, Funding acquisition. Zhongchao Chen:and Xiangbei Zhou Calculation, Writing – original draft. Mei Zhou and Zhen Li: Data collection and analysis, Calculation.

Acknowledgement

This project is a part of the Fifth Forest Management Inventory Project of Guangxi Zhuang Autonomous Region (5^th FMI-GX, 2017-2020), China. The airborne LiDAR data acquisition and prepossessing and field data measurement were founded by the finance department of the Guangxi Zhuang Autonomous Region. The authors would like to express their sincere gratitude to Chengling Yang and Yao Liang from the Guangxi Forest Inventory and Planning Institute (FIPI-GX) and 120 field crews who worked on the field measurements. The authors also thank Guangxi 3D Remote Sensing Engineering Technology Co., Ltd., Feiyan Aero Remote Sensing Technology Co., Ltd., Zhongke Remote Sensing Technology Group Co., Ltd., and Guangzhou Jiantong Surveying and Mapping Geographic Information Technology Co., Ltd., which were responsible for the acquisition and preprocessing of the airborne LiDAR data utilized herein. In addition, the authors are grateful to Professor Guofan Shao from the Department of Forestry and Natural Resources of Purdue University in West Laffayette, IN, USA, and Dr. Marc Bouvier from UMR TETIS Irstea-Cirad-AgroParisTech/ENGREF, Maison de la Télédétection en Languedoc-Roussillon, 500, rue J.F. Breton BP 5095, 34196 Montpellier Cedex 05, France, for their advice on the name of the generalization model. We acknowledge the anonymous reviewers and the editor for their insightful suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available at https://doi.org/10.57760/sciencedb.j00001.00467

Additional information

Funding

This work received financial support from the Forest Department of Guangxi Zhuang Autonomous Region, China [GXLYKJ2016-001].

Appendix

Table A1. 86 formulations constructed with a rule-based exhaustive combination of 13 LiDAR variables.

No. of model	hmean	hp95	hstdev	hcv	CC	dp50	dp75	LADmean	LADstdev	LADcv	VFPmean	VFPstdev	VFPcv
1	√				√
2		√			√
3	√				√			√
4	√				√				√
5	√				√					√
6	√				√						√
7	√				√							√
8	√				√								√
9		√			√			√
10		√			√				√
11		√			√					√
12		√			√						√
13		√			√							√
14		√			√								√
15	√		√		√			√
16	√		√		√				√
17	√		√		√					√
18	√		√		√						√
19	√		√		√							√
20	√		√		√								√
21	√			√	√			√
22	√			√	√				√
23	√			√	√					√
24	√			√	√						√
25	√			√	√							√
26	√			√	√								√
27		√	√		√			√
28		√	√		√				√
29		√	√		√					√
30		√	√		√						√
31		√	√		√							√
32		√	√		√								√
33		√		√	√			√
34		√		√	√				√
35		√		√	√					√
36		√		√	√						√
37		√		√	√							√
38		√		√	√								√
39	√		√		√	√		√
40	√		√		√	√			√
41	√		√		√	√				√
42	√		√		√	√					√
43	√		√		√	√						√
44	√		√		√	√							√
45	√		√		√		√	√
46	√		√		√		√		√
47	√		√		√		√			√
48	√		√		√		√				√
49	√		√		√		√					√
50	√		√		√		√						√
51	√			√	√	√		√
52	√			√	√	√			√
53	√			√	√	√				√
54	√			√	√	√					√
55	√			√	√	√						√
56	√			√	√	√							√
57	√			√	√		√	√
58	√			√	√		√		√
59	√			√	√		√			√
60	√			√	√		√				√
61	√			√	√		√					√
62	√			√	√		√						√
63		√	√		√	√		√
64		√	√		√	√			√
65		√	√		√	√				√
66		√	√		√	√					√
67		√	√		√	√						√
68		√	√		√	√							√
69		√	√		√		√	√
70		√	√		√		√		√
71		√	√		√		√			√
72		√	√		√		√				√
73		√	√		√		√					√
74		√	√		√		√						√
75		√		√	√	√		√
76		√		√	√	√			√
77		√		√	√	√				√
78		√		√	√	√					√
79		√		√	√	√						√
80		√		√	√	√							√
81		√		√	√		√	√
82		√		√	√		√		√
83		√		√	√		√			√
84		√		√	√		√				√
85		√		√	√		√					√
86		√		√	√		√						√