Prediction of Forest Structural Parameters Using Airborne Full-Waveform LiDAR and Hyperspectral Data in Subtropical Forests

Abstract

Accurate acquisition of forest structural parameters, which is essential for the parameterization of forest growth models and understanding forest ecosystems, is also crucial for forest inventories and sustainable forest management. In this study, simultaneously acquired airborne full-waveform (FWF) LiDAR and hyperspectral data were used to predict forest structural parameters in subtropical forests of southeast China. The pulse amplitude and waveform shape of airborne FWF LiDAR data were calibrated using a physical process-driven and a voxel-based approach, respectively. Different suites of FWF LiDAR and hyperspectral metrics, i.e., point cloud (derived from LiDAR-waveforms) metrics (DPC), full-waveform (geometric and radiometric features) metrics (FW) and hyperspectral (original reflectance bands, vegetation indices and statistical indices) metrics (HS), were extracted and assessed using correlation analysis and principal component analysis (PCA). The selected metrics of DPC, FW and HS were used to fit regression models individually and in combination to predict diameter at breast height (DBH), Lorey’s mean height (H_L), stem number (N), basal area (G), volume (V) and above ground biomass (AGB), and the capability of the predictive models and synergetic effects of metrics were assessed using leave-one-out cross validation. The results showed that: among the metrics selected from three groups divided by the PCA analysis, twelve DPC, eight FW and ten HS were highly correlated with the first and second principal component (r > 0.7); most of the metrics selected from DPC, FW and HS had weak relationships between each other (r < 0.7); the prediction of H_L had a relatively higher accuracy (Adjusted-R² = 0.88, relative RMSE = 10.68%), followed by the prediction of AGB (Adjusted-R² = 0.84, relative RMSE = 15.14%), and the prediction of V had a relatively lower accuracy (Adjusted-R² = 0.81, relative RMSE = 16.37%); and the models including only DPC had the capability to predict forest structural parameters with relatively high accuracies (Adjusted-R² = 0.52–0.81, relative RMSE = 15.70–40.87%) whereas the usage of DPC and FW resulted in higher accuracies (Adjusted-R² = 0.62–0.87, relative RMSE = 11.01–31.30%). Moreover, the integration of DPC, FW and HS can further improve the accuracies of forest structural parameters prediction (Adjusted-R² = 0.68–0.88, relative RMSE = 10.68–28.67%).

1. Introduction

As the dominant terrestrial ecosystem on earth, forests occupy approximately 30% of the land surface area and contribute to 75% of land gross primary production [1,2]. Subtropical forests have high diversity, dense carbon and complex structure, and cover approximately one quarter of China’s total area. They provide valuable ecosystem goods and services to humanity and play a key role in the mitigation of climate change [3,4]. Forest structures, shaped by silvicultural practices and natural events, provide considerable information about ecosystem values such as biodiversity, water conservation and erosion control [5]. Forest structural parameters (e.g., diameter at breast height (DBH), tree height (H), tree density (N), basal area (G), etc.) are essential for the parameterization of forest growth models and understanding forest ecosystems [6,7]. Timely, accurate and reliable acquisition of forest structural parameters across large areas is crucial for sustainable and multifunctional forest management [8]. Traditionally, forest structural parameters were collected using field inventory, which is labor intensive and time-consuming [9,10]. Remote sensing technology is able to provide detailed continuous-spatial, multi-dimensional and massive-spectral information, allowing for precise forest structural parameter prediction based on their structural and spectral signatures [11,12,13]. Remote sensing data have advantages such as spatial information quantification, high geometric precision and vast geographic coverage [14,15], and have been used in the prediction of forest structural parameters over a range of forest types [16,17,18].

Light detection and ranging (LiDAR) has been applied as a promising technology for predicting forest structural parameters due to its capability to provide three-dimensional information about forest structures with high accuracy [19,20,21]. Airborne discrete-return LiDAR systems record multiple return signals, which contain the three-dimensional position and intensity of reflected light from each transmitted pulse. The metrics extracted from discrete point clouds represent the vertical structural characteristics of canopy such as height measures, Weibull distribution parameters, and classes of crown volume zones, etc., which can be used to predict forest structural parameters [22,23,24]. Forest structural parameters, i.e., H, G, V and AGB, have been predicted using discrete point cloud data in previous studies [21,22,24,25]. However, discrete-return systems record limited information in each returned signal, and can only detect the surfaces which are separated sufficiently in space [26]. Airborne full-waveform (FWF) LiDAR systems record the whole backscattered returns, thus they can record the geometric and biophysical attributes of forests [27]. The point cloud can be derived from FWF LiDAR data, and the number of returns extracted from FWF data is much higher than that from discrete point cloud data [28]. The waveforms obtained by full-waveform systems depend on many factors such as target backscattering characteristics, LiDAR sensor types, scan geometry [29], etc. In order to make pulse amplitude possible to be used as its true value, the pulse amplitude needs to be calibrated and corrected in full-waveform data processing [30,31]. Owing to the waveform recorded by full-waveform system being stretched by the increases of off-nadir angle and the waveforms from different trajectories often being non-vertical, the waveform processing approach needs to be developed to synthesize the raw waveforms in multiple directions into composite vertical waveforms. The full-waveform metrics (FW) extracted from FWF data describe the canopy response using parameters of waveform shape and can be applied in the prediction of forest structural parameters [13,32]. Lindberg et al. [33] extracted the amplitude of waveforms in height intervals to predict the volume in hemi-boreal forest in the southwest of Sweden, and found that the predictive model of total volume had high accuracy for waveform data (relative RMSE = 31.9%). In the prediction of forest structural parameters (e.g., DBH, H, N and AGB, etc.), the accuracies of predictive models were improved when using full-waveform metrics [34,35].

Hyperspectral data offer large amounts of continuous-narrow bands which contain detailed spectral signatures associated with forest biophysical properties, and which can be applied to predict forest structural parameters [16,36,37]. Airborne hyperspectral data usually have finer spatial resolution than space-borne data because airborne platforms commonly having lower altitudes than space-borne platforms [38,39]. Previous studies have demonstrated that airborne hyperspectral data perform well in forest tree species classification [40,41], and structural parameters prediction, e.g., DBH [42], tree height [43], basal area [44], stem density [23] and biomass [44,45,46]. The visible (VIS) and near-infrared (NIR) regions of hyperspectral narrow bands are usually considered to be correlated with forest structure properties [39,47]. Latifi et al. [23] found that a number of atmospheric window bands in the domains of VIS and NIR such as 540 nm, 680–730 nm and 970 nm were the most important and stable predictors of forest stem density and biomass. The hyperspectral narrow band metrics, formulated using the bands in VIS and NIR domains, rely on the pigments (e.g., chlorophyll, carotene and anthocyanin, etc.), structure and physiology of tree canopy and have great potential in the prediction of forest structural parameters [48,49]. Particularly, using hyperspectral narrow band metrics can relatively weaken the influences of soil background reflectance, illumination and atmospheric absorption [17,50]. The vegetation indices are the most commonly used narrow band metrics in the prediction of forest structural parameters [37,51,52]. The vegetation indices such as the soil-adjusted vegetation index (SAVI), atmospherically resistant vegetation index (ARVI), and normalized difference vegetation index (NDVI) are considered to be correlated with N, G, V and AGB [23,53]. Zhang et al. [54] used HJ-1 hyperspectral data to predict forest AGB in a subtropical forest and found that SAVI is strongly related to the AGB (r = 0.91). Nevertheless, in a densely forest area, the hyperspectral data metrics are prone to asymptotically reach a saturation level [55]. Moreover, since the hyperspectral data usually provide horizontal information, it has certain limitations in quantifying the vertical structure of forests [56]. These limitations can influence the accuracy of predictions of forest structural parameters.

An integration of airborne LiDAR and hyperspectral data is expected to provide more information on the prediction of forest structural parameters. However, few studies have attempted to improve forest structural parameters predictions by integrating airborne LiDAR and hyperspectral data. Dalponte et al. [42] integrated airborne discrete point cloud LiDAR and hyperspectral data to predict stem diameter and volume in a temperate forest and found that the improvement in accuracies of stem diameter and volume were 0.4% and 0.5%, respectively. Luo et al. [57] integrated airborne discrete point cloud LiDAR and hyperspectral data to estimate AGB in a northern temperate deciduous forest. The results indicated that by using integrated discrete point cloud LiDAR and hyperspectral data, 2.2% more of the variability in AGB was explained. In previous studies, most have only focused on the prediction of forest biomass by integrating discrete-return LiDAR and hyperspectral data in temperate and boreal forests. Since the acquisition modalities and structures of FWF LiDAR and hyperspectral data are disparate, it is difficult to integrate FWF LiDAR and hyperspectral data at the raw data level [58]. The integration commonly includes transforming FWF LiDAR data into two-dimensional images and adding suites of hyperspectral metrics or the information of tree-species derived from hyperspectral data into FWF data [23,32,53].

However, the predictions of forest structural parameters have been mainly implemented in temperate and boreal forests, and there are few published studies from subtropical forest. Moreover, few studies have used radiometrically and geometrically calibrated FWF LiDAR data in analysis, therefore, the synergetic effects of FWF LiDAR and hyperspectral data could be influenced by the target backscattering characteristics and scan geometry. In addition, previous studies did not comprehensively extract and optimize suites of FWF LiDAR and hyperspectral metrics, and the relativities and synergetic effects of point cloud, full-waveform and hyperspectral metrics have not been fully explored in the prediction of forest structural parameters. To the best of our knowledge, no previous study has integrated simultaneously acquired airborne FWF LiDAR and hyperspectral data to predict forest structural parameters in subtropical forests. The objectives of this paper are: (1) to calibrate the airborne FWF LiDAR data by the physical process-driven and voxel-based models; (2) to integrate and assess the synergetic effects of FWF LiDAR and hyperspectral data-derived metrics for predicting forest structural parameters in subtropical forests; and (3) to validate the predictive models fitted by DPC, FW and HS individually, and in combination using field measured data and to analyze the residuals of the prediction.

2. Materials and Methods

Figure 1 shows the workflow for forest structural parameters prediction. First, the FWF LiDAR data and hyperspectral data were preprocessed to extract discrete point cloud data, calibrate the pulse amplitude and waveform shape, and reduce the influences of atmospheric interference and terrain distortion. Second, different suites of FWF LiDAR and hyperspectral metrics, i.e., point cloud (derived from full-waveform) (DPC) metrics and full-waveform (geometric and radiometric features) (FW) metrics and hyperspectral (original reflectance bands, vegetation indices and statistical indices) (HS) metrics, were extracted and selected using correlation analysis and principal component analysis (PCA), and the relativities of selected metrics were assessed using Pearson’s correlation analysis. Finally, the selected DPC, FW and HS were used to fit regression models individually, and in combination to predict DBH, H_L, N, G, V and AGB, and the capability of the predictive models and synergetic effects of metrics were assessed.

2.1. Study Area

The nearly 1103 ha study area is located in the Yushan forest (120°42′9.4″ E, 31°40′4.1″ N), situated in the southern Jiangsu provinces, southeast China (Figure 2). The annual mean temperature and precipitation are 15.6 °C and 1062.5 mm, respectively. The elevation of the Yushan forest ranges from 20 to 261 m above sea-level. The Yushan forest is north subtropical secondary forest and has three types of forests: coniferous tree species dominated, broadleaved tree species dominated and mixed tree species forests [59]. Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.) and Masson pine (Pinus massoniana Lamb.) are the main coniferous tree species. Sweet gum (Liquidambar formosana Hance) and Sawtooth oak (Quercus acutissima Carruth.) are the major broadleaved tree species in the study area.

2.2. Field Data

Field surveys were conducted under leaf-on condition in June and August 2012 and August 2013. Guided by the pre-stratified stand inventory data in 2012, a total of 67 square (30 × 30 m²) field plots were established. These plots covered multiple site indices, age classes and tree species, which can be classified into three types according to the composition of the tree species: (i) coniferous tree species forest (n = 15); (ii) broadleaved tree species forest (n = 18); and (iii) mixed tree species forest (n = 34).

The coordinates of the plot corners were acquired using Trimble GPS measurements with the result of sub-meter accuracy. All the live trees within each plot, which have a DBH > 5 cm were measured. The measurement of individual tree parameters can be seen in [41]. The dead wood and small trees which has a DBH < 5 cm within the plot were also recorded, but excluded in the calculations of biomass and volume. The six plot-level forest structural parameters, including DBH, H_L, N, G, V and AGB, were calculated using the measured individual tree data. Species-specific allometric equations and general volume equations of local or nearby provinces were used to calculate AGB and V, respectively (

Table A1. Allometric equations used in this study for aboveground biomass calculation.

Tree Species	Component	a	b	R2	References
Masson pine	Stem wood (Ws)	0.141	1.092	0.9970	[102]
	Live branches (Wb)	0.065	0.991	0.9871
	Foliage (Wf)	0.132	0.745	0.9827
Chinese fir	Stem wood (Ws)	0.124	0.680	0.9704	[103]
	Live branches (Wb)	0.203	0.385	0.7223
	Foliage (Wf)	0.850	0.189	0.6567
Slash pine	Stem wood (Ws)	0.235	0.900	0.9523	[104]
	Live branches (Wb)	0.080	1.064	0.8520
	Foliage (Wf)	0.456	0.610	0.8802
Sawtooth oak	Stem wood (Ws)	0.018	1.034	0.9864	[105]
	Live branches (Wb)	0.00008	1.468	0.9745
	Foliage (Wf)	0.004	0.769	0.8662
Sweet gum	Stem wood (Ws)	0.093	0.801	0.9310	[106]
	Live branches (Wb)	0.083	0.649	0.9890
	Foliage (Wf)	1.084	0.217	0.6940
Other broadleaves	Stem wood (Ws)	0.023	0.985	0.9903	[107]
	Live branches (Wb)	0.00004	3.785	0.9623
	Foliage (Wf)	0.00003	1.378	0.9456

Note: Each biomass component was calculated using the equation of W = a(D²H)^b, where H is the tree height (m), D is the DBH (cm) and a, b are the parameters.

and

Table A2. Equations used in this study for volume calculation.

Tree Species	a	b	c	R2	References
Masson pine	0.000056	1.9980	0.8470	0.9649	[108]
Chinese fir	0.000062	1.8520	1.0080	0.9956	[109]
Slash pine	0.000067	1.9952	0.7908	0.9523	[110]
Sawtooth oak	0.000033	2.0000	1.0000	0.9930	[111]
Sweet gum	0.000053	1.8822	1.0093	0.9810	[112]
Other broadleaves	0.000067	1.9600	0.8144	0.9900	[111]

Note: Tree volume was calculated using the equation of V = aD^bH^c, where H is the tree height (m), D is the DBH (cm) and a, b, c are the parameters.

). Within each plot, the AGB and V of each individual tree were calculated according to the DBH and H measured in the field, and then summed to the plot-level AGB and V.

Table 1. Summary statistics of field-measured forest structural parameters in the study area (n = 67, size = 900 m²).

Parameters	Coniferous Forests (n = 15)		Broadleaved Forests (n = 18)		Mixed Forests (n = 34)
	Mean	Std. Dev.	Mean	Std. Dev.	Mean	Std. Dev.
DBH	13.92	3.09	16.90	4.01	14.33	3.16
HL	10.41	1.84	12.25	3.00	10.60	1.57
N	1692.73	619.81	986.28	441.67	1467.59	413.82
G	25.96	6.01	24.44	6.56	25.23	6.24
V	124.80	24.51	142.05	34.31	126.71	26.91
AGB	77.87	23.27	98.69	47.45	88.81	30.04

Note: DBH = diameter at breast height (cm); H_L = Lorey’s mean height (m); N = stem number (ha⁻¹); G = basal area (m²·ha⁻¹); V = volume (m³·ha⁻¹); AGB = above ground biomass (Mg·ha⁻¹). Std. Dev. = standard deviation.

provides a summary of the six forest structural parameters at plot-level.

2.3. Remote Sensing Data

In August 2013, the airborne full-waveform LiDAR and hyperspectral data were simultaneously obtained using the LiCHy System [60]. The platform was flown at the height of 900 m above ground and the datasets covered the whole Yushan Forest. Full-waveform LiDAR data were obtained using the Riegl LMS-Q680i scanner. The scanning angle was ±15° from nadir, and the pulse repetition frequency was 360 kHz. The temporal sample spacing for recording returned waveforms was 1 ns (15 cm in distance approximately), and the size of the footprint at nadir was 0.45 m in diameter. In the overlapping regions, the pulse density was three times higher than a single strip. Hyperspectral data were acquired using an AISA Eagle II sensor with 64 bands and the spectral resolution was 3.3 nm. The sensor obtained hyperspectral images in the pattern of push-broom imaging and the spectrum ranges covered from 400 nm to 970 nm. The spatial and radiative resolution of the hyperspectral data were 0.6 m and 12 bit, respectively. The geometric accuracy of each pixel was less than one meter with an inertial measurement unit (IMU), which utilized real-time differential corrections by a 12-channel GPS receiver.

Table 2. The summary of flight parameters and properties of the full-waveform LiDAR and hyperspectral sensors.

Parameters	Riegl LMS-Q680i	AISA Eagle II
Data of acquisition	17 August 2013	17 August 2013
Flight height	900 m	900 m
Flight speed	55 m/s	55 m/s
Beam divergence	0.5 mrad	–
IFOV	-	0.65 mrad
Spatial resolution	0.45 m	0.6 m
Wavelength	1550 nm	400–970 nm
Bands	1	64
Swath width	1040 m	612 m
Bit depth	–	12 bits
Average pulse distance	0.49	–
Average point density	8.37	–

summarizes the characteristics of the full-waveform LiDAR and hyperspectral sensors.

2.4. Remote Sensing Data Pre-Processing

2.4.1. Full-Waveform LiDAR Data Pre-Processing

First, a de-noising algorithm and a Gaussian filter were applied to suppress and smooth the background noise of each returned waveform. The full width at half-maximum (FWHM) was used to calculate the kernel size of Gaussian filter [61]. Then, the locations and amplitudes of each peak within the waveform were extracted using a local maxima peak detection filter [62]. Finally, the returned waveform was decomposed using the Gaussian decomposed algorithm.

The LiDAR point clouds can be derived from the FWF LiDAR data using the Gaussian decomposed algorithm. Generally, the following equation can be used to decompose the backscattered waveform into Gaussian components:

$$f\left(x\right)=b+{{\displaystyle \sum }}_{i=1}^n{a}_i{e}^{-{\left(x-t_i\right)}^2/2{\sigma }_i^2}$$

(1)

where f(x) stands for the returned waveform, b represents the background noise, and n is the number of decomposed Gaussian components. The a_i, t_i and σ_i are the parameters corresponding to pulse amplitude, time of round trip, and the pulse width, respectively [31]. Then, the Levenberg-Marquardt algorithm and a nonlinear least squares method were applied to fit multiple Gaussian components into the backscattered waveform. The LiDAR point clouds extracted from FWF LiDAR data were stored as the format of LAS 1.3 and used for analysis.

In this study, the points in the ground and upper surface of the forest canopy were applied to create the digital terrain model (DTM) and digital surface model (DSM), respectively. The cell size of DTM and DSM was 0.6 m, the same as the resolution of the hyperspectral data. The value in each cell was calculated as the mean elevation of these points, and the cells which had no points were interpolated using neighboring cells by a linear interpolation approach. The value of the DTM was subtracted from each point elevation to calculate the normalized point cloud of whole study area.

The returned pulse width (W_i^Г, the standard deviation of pulse) and amplitude (I_i^Г, the integral of returned waveform, which represent the pulse energy) were derived from Gaussian components [63]. In this study, the pulse width (W_i^Г) and amplitude (I_i^Г) were calibrated using a physical process-driven approach [31]. The values of the pulse width (W^e) and amplitude (I^e) of the scanner emitted pulses were used to calibrate W_i^Г and I_i^Г, and the I_i^Г was corrected for the loss of signal using the distance between the sensor and the object (D_i) and the normalization distance of D_o [30]:

$$W_i^c=\frac{W_i^{\mathsf{\Gamma }}}{W^e}$$

(2)

$$I_i^c=\frac{I_i^{\mathsf{\Gamma }}}{I^e}\times {\left(\frac{D_i}{D_o}\right)}^k$$

(3)

where W_i^c and I_i^c is the calibrated pulse width and amplitude, respectively. The value of k which depended on the attenuation of signal occur in the atmosphere was set to 2 [30], and the D_o was set to 900 m (the mean height of the platform).

It has been demonstrated in previous studies that the waveform is stretched by the increase in the off-nadir angle [64,65]. Moreover, due to the obtained airborne FWF LiDAR data are normally comprised of multiple overlapping strips, the waveforms in a specific location may come from several strips [66,67]. In this study, a voxel-based approach to composite waveforms was used to correct FWF data to avoid the influences of off-nadir angle in the waveform shape and to integrate non-vertical waveforms from multiple strips into composited vertical waveforms. This approach first decomposed the forest canopies into voxels by vertical space partition (0.6 × 0.6 × 0.3 m³), and then synthesized raw waveforms from multiple strips into composite vertical waveforms using the maximum amplitude value in each voxel (Figure 3). Each composite vertical waveform was normalized using the digital terrain model (DTM).

2.4.2. Hyperspectral Data Pre-Processing

The radiance hyperspectral images covering the whole Yushan forest were geometrically rectified with Global Navigation Satellite System (GNSS) and Inertial Navigation System (INS) data. Then, the geometric-rectified images were mosaicked into a single scene. Atmospheric correction was applied using an empirical line model, combined with field-measured reflectance spectra of different target objects obtained by ASD FieldSpec spectrometer, to get the surface reflectance of covers. In this study, the FWF LiDAR and hyperspectral data were integrated at feature level based on a common coordinate frame. The framework of integration usage of these two datasets can be seen in Figure 3. In order to have the best possible geographical matches between the FWF LiDAR and hyperspectral data, the hyperspectral data were co-registered to the digital surface model (DSM) which calculated from FWF LiDAR data. In the area of each plot, more than 30 ground control points (GCPs) were used on the hyperspectral image (30 × 30 m²). The root mean square error of co-registration was lower than 0.3 m (half of one pixel).

2.5. Full-Waveform LiDAR Metrics

2.5.1. Point Cloud Metrics

The metrics derived from the height normalized LiDAR point cloud were applied to describe the canopy structure of the plots. In this study, the calculated point cloud metrics (DPC) were: (i) the selected height measures (n = 11); (ii) the Weibull parameters fitted to the profile of apparent foliage density (n = 2); and (iii) the crown volume zones (n = 4). A summary of the point cloud metrics and their descriptions is given in

Table 3. The summary of the FWF LiDAR metrics (i.e., DPC and FW) extracted from FWF LiDAR data (the code and description of each metric are listed).

FWF LiDAR Metrics	Description
a Point cloud metrics (DPC)
Percentile heights (h25, h50, h75 and h95)	The percentiles of the canopy height distributions (25th, 50th, 75th and 95th).
Canopy return density (d1, d3, d5, d7 and d9)	The proportion of points above the height percentiles (10th, 30th, 50th, 70th and 90th).
Coefficient of variation of heights (hcv)	Coefficient of variation of heights of all points.
Canopy cover above 2 m (Cover)	Percentage of all points above 2 m.
α and β parameter of Weibull distribution (i.e., Wα and Wβ)	The α and β parameters of the Weibull distribution fitted to foliage density profile.
Open and Closed gap zones of canopy volume models (CVM) (i.e., Og and Cg)	The empty voxels located above and below the canopy respectively.
Euphotic and Oligophotic zones of CVM (i.e., E and O)	The voxels located within an uppermost percentile (65%) of all filled grid cells of that column, and voxels located below the point in the profile.
b Full-waveform metrics (FW)
Height of median energy (HOME)	The distance from waveform centroid to the ground.
Waveform distance (WD)	The distances from waveform beginning to the ground.
Vertical distribution ratio (VDR)	The differences between the WD and the HOME, divided by WD.
Number of peaks (NP)	The number of detected peaks within each normalized composite waveform.
Roughness of outermost canopy (ROUGH)	The distance from the waveform beginning to the first peak.
Front slope angle (FS)	The vertical angles from waveform beginning to the first peak of canopy return energies.
Return waveform energy (RWE)	The total received energy, i.e., the area below the waveform between beginning and end.
Intensity of Gaussian component (Int)	Mean of the intensity of Gaussian components within one waveform.
Full width at half maximum (FWHM)	Full width at half maximum of one waveform.

Note: ^a These DPC were calculated using point cloud derived from FWF LiDAR data; ^b The FW were extracted from the waveforms which were pre-processed and calibrated. See text for details.

.

To exclude the influences of below-canopy and non-canopy returns, the point cloud metrics such as percentile heights and canopy return densities were calculated using the points that were two meters above ground [68]. The parameters α and β of the Weibull curve were extracted from the profile of apparent foliage density as follows [22]:

$$L\left(z\right)=1-\left[e^{-{\left(\frac{1-z/H}{\alpha }\right)}^{\beta }}\right]$$

(4)

where α and β are the parameters of Weibull, z is the height and H is the maximum canopy height in a plot.

The zones of crown volume model (i.e., O_g, C_g, E and O) were used to characterize the forest crown volume and spatial arrangement of the canopy materials in three-dimensions [69]. First, the forest canopy was decomposed into a matrix of voxels, the size of each voxel was 0.6 × 0.6 × 0.3 m³. Second, the voxels within the matrix were classified into “filled” if there was energy returned from the voxel and classified into “empty” if there was no energy returned from the voxel. Third, the “filled” voxels were classified as “euphotic” and “oligophotic” depending on whether the voxel was above or below the threshold height of the uppermost 65% for all “filled” voxels. Finally, the “empty” voxels were classified into “open” and “closed” gap zones depending on whether they were located above or below the filled voxels.

2.5.2. Full-Waveform Metrics

Full-waveform metrics (FW) provide three-dimensional forest structure information by extracting radiometric and geometric properties of recorded backscattered waveforms. In this study, 18 full-waveform metrics including the mean (μ) and standard deviation (σ) within each plot were extracted from composite waveforms (Figure 3). First, the full-waveform metrics (Table 3) of each composite waveform were calculated; second, the mean and standard deviation of all the full-waveform metrics in each plot were calculated as the full-waveform metrics at plot-level. Table 3 gives the summary of these full-waveform metrics and descriptions.

2.6. Hyperspectral Metrics

The hyperspectral metrics (HS) are good indices in the prediction of forest structural properties, due to their ability to describe crown structures, which are related to vegetation pigments, physiology and stress, directly or indirectly. In this study, 112 hyperspectral metrics were derived from the preprocessed hyperspectral image, including: (i) reflectance values from AISA Eagle II channels; (ii) vegetation indices; (iii) first 10 components of the principal component transformation (PCT), independent components transformation (ICT) and minimum noise fraction transformation (MNF).

The spectral reflectance was strongly correlated with the structural properties (e.g., leaf area index, the amount of biomass and spatial arrangement of structures) of forests [70]. In this study, all channels in the domains of VIS, RE, and NIR were used. The mean values of 50 × 50 pixels within the plots were calculated from the reflectance of the entire 64 channels. The same procedure for calculation was followed to extract the other hyperspectral metrics (vegetation indices and first 10 components of PCT, ICT and MNF).

Hyperspectral vegetation indices, which rely on specific absorption features, are the most commonly used narrow band metrics. The vegetation indices calculated from the hyperspectral image have great advantages in predicting forest structural parameters [23,49]. In this study, 18 vegetation indices were extracted and summarized in

Table 4. A summary of the vegetation indices with respective equations and references. Wavelengths chosen were the closest AISA wavelengths to the equations in the cited literature.

Vegetation Index	Equation	Reference
Simple ratio (SR)	$\rho 801/\rho 676$	[71]
Normalized difference vegetation index (NDVI)	$\left(\rho 801-\rho 676\right)/\left(\rho 801+\rho 676\right)$	[72]
Enhanced vegetation index (EVI)	$2.5\times \left[\left(\rho 801-\rho 676\right)/\left(1+\rho 801+6\times \rho 676-7.5\times \rho 479\right)\right]$	[73]
Green normalized difference vegetation index (GNDVI)	$\left(\rho 801-\rho 553\right)/\left(\rho 801+\rho 553\right)$	[74]
Soil adjusted vegetation index (SAVI)	$1.5\times \left(\rho 801-\rho 676\right)/\left(\rho 801+\rho 676+0.5\right)$	[75]
Atmospherically resistant vegetation index (ARVI)	$\left[\rho 801-\left(2\times \rho 676-\rho 479\right)\right]/\left[\rho 801+\left(2\times \rho 676-\rho 479\right)\right]$	[76]
Red-edge vegetation stress index (RVSI)	$\left[\left(\rho 724+\rho 753\right)/2\right]-\rho 733$	[77]
Plant senescence reflectance index (PSRI)	$\left(\rho 676-\rho 498\right)/\rho 753$	[78]
Vogelmann red edge index 1 (VOG1)	$\rho 743/\rho 724$	[79]
Vogelmann red edge index 2 (VOG2)	$\left(\rho 733-\rho 743\right)/\left(\rho 714-\rho 724\right)$	[79]
Red green ratio index (RGRI)	$RED/GREEN$	[80]
Photochemical reflectance index (PRI)	$\left(\rho 534-\rho 572\right)/\left(\rho 534+\rho 572\right)$	[81]
Photochemical reflectance ratio (PRR)	$\rho 534/\rho 572$	[82]
Water band index (WBI)	$\rho 898/\rho 975$	[83]
Carotenoid reflectance index 1 (CRI1)	$\left(1/\rho 507\right)-\left(1/\rho 553\right)$	[84]
Carotenoid reflectance index 2 (CRI2)	$\left(1/\rho 507\right)-\left(1/\rho 705\right)$	[84]
Anthocyanin reflectance index 1 (ARI1)	$\left(1/\rho 553\right)-\left(1/\rho 705\right)$	[85]
Anthocyanin reflectance index 2 (ARI2)	$\rho 801\times \left[\left(1/\rho 553\right)-\left(1/\rho 705\right)\right]$	[85]

.

The principal component analysis (PCA), minimum noise fraction analysis (MNF) and independent components analysis (ICA) are three algorithms which are commonly used to de-noise and extract primary information from hyperspectral images. We used these three approaches to calculate 192 components (64 components for each approach), out of which we used the first 10 of each approach for analysis, to explore whether there was a component that summarized the forest structural parameters-related channels to one value, and therefore, to ensure the models’ conciseness [48,86].

2.7. Metrics Optimization and Regression Analysis

Previous studies have demonstrated that optimization of the candidate metrics can reduce irrelevant and redundant information and help create highly efficient, transferable and robust productive models. In this study, all of the FWF LiDAR and hyperspectral metrics were first optimized using correlation analysis. The 15 metrics, which had relatively high correlation with the forest structural parameters, were correspondingly selected from point cloud metrics (DPC), full-waveform metrics (FW) and hyperspectral metrics (HS). Then, the 45 metrics (15 point cloud, 15 full-waveform and 15 hyperspectral metrics) were analyzed using the biplot of PCA, which can be used to select the important metrics in the clusters [87]. The 12 point cloud, 8 full-waveform and 10 hyperspectral metrics which highly correlated with the first and second principal component (r > 0.7) and were selected from the three groups divided by the PCA. Finally, the five metrics that had the highest correlations with the first and second principal component in each group were selected as the best metrics to fit the combo models.

The backward stepwise regression approach was applied to relate FWF LiDAR and hyperspectral metrics to field-measured forest structural parameters. In the models, three predictor variables at the 5% significance level were selected. To ensure the metrics in the models had no serious collinearity, the models which had the condition number (k) < 30 were selected. Finally, according to the value of the Akaike information criterion (AIC), the best fitting models with the lowest AIC were selected.

In the study, three types of predictive models of DBH, H_L, N, G, V and AGB were developed using DPC, FW, HS and an integration of two or three of these for the combo models. First, the DPC models (DPC based models) were fitted using 12 DPC alone to predict the six forest structural parameters; second, the FW models (DPC and FW based models) were fitted using the integration of 12 DPC and 8 FW to predict the six forest structural parameters; third, the combo models were fitted using the integration of the best metrics (each of the five metrics selected from PCA groups) to predict the six forest structural parameters. All of the models were assessed by adjusted coefficient of determination (Adj-R²), Root-Mean-Square-Error (RMSE), and relative RMSE (rRMSE). The leave-one-out (LOO) cross validation was applied to assess the accuracy of prediction models and assess the synergetic effects of FWF LiDAR and hyperspectral metrics.

3. Results

3.1. Full-Waveform LiDAR and Hyperspectral Metrics Extraction

All of the metrics, including the point cloud metrics (n = 17), full-waveform metrics (n = 18) and hyperspectral metrics (n = 112) were extracted from the LiDAR point cloud (with height normalized), composite waveform and preprocessed hyperspectral image, respectively. Figure 4 shows the profiles of point cloud and apparent foliage (I), profiles of intensity of energy and composite waveform (II), and the spectral reflectance from 400 nm to 1000 nm (III). The profiles of apparent foliage and Weibull distribution appropriately describe the vertical distribution of point cloud. The peak of the Weibull distribution curve in broadleaved forest plot is relatively higher (height = 8.55 m), followed by the mixed forest plot (height = 6.34 m), and the peak of the Weibull distribution curve in coniferous plots is relatively lower (height = 5.41 m). The profiles of composite waveform truly describe the space distribution of energy, and the height of peaks of composite waveform are same as the height of energy concentration. The profiles of intensity of energy and composite waveform are similar to the profiles of point cloud and Weibull distribution, respectively. The spectral reflectance describes the biophysical and biochemical properties of the canopy. The spectral reflectance in the broadleaved forest plot is highest (mean = 0.06–0.52), followed by spectral reflectance in the coniferous forest plot (mean = 0.06–0.41), and spectral reflectance in the mixed forest plot is lowest (mean = 0.05–0.33).

3.2. Metrics Selection and Optimization

Figure 5 shows the projection of the first two PCA scores from the selected point cloud (n = 15), full-waveform (n = 15) and hyperspectral (n = 15) metrics. The first (PCA1) and second (PCA2) components of PCA account for 57.61% and 31.02% of the total variance, respectively. In the four groups divided by PCA, 37 metrics in three groups were highly correlated with the PCA1 and PCA2. The 12 point cloud (h₂₅, h₅₀, h₇₅, h₉₅, d₁, d₅, d₉, h_cv, Cover, W_α, E, O), 8 full-waveform (HOME_μ, WD_μ, VDR_μ, NP_μ, RWE_μ, Int_μ, HOME_σ, WD_σ) and 10 hyperspectral (B45, B63, VOG1, NDVI, RVSI, EVI, ARVI, CRI1, PCA1, PCA2) metrics were correlated with the PCA1 and PCA2 higher than 0.7. Most of the point cloud and full-waveform metrics were in the opposite direction. Moreover, most of hyperspectral metrics and FWF LiDAR metrics (point cloud and full-waveform metrics) were orthonormal. Therefore, the FWF LiDAR metrics and hyperspectral metrics could be used as complementary metrics in forest structural parameters prediction.

The five point cloud metrics (h₅₀, h₇₅, d₁, E, O), five full-waveform metrics (HOME_μ, VDR_μ, NP_μ, RWE_μ, and WD_σ) and five hyperspectral metrics (NDVI, RVSI, ARVI, CRI1 and PCA2), which had the highest correlation with PCA1 and PCA2 in each group are shown in Figure 6. According to the result of the correlation analysis, most of the metrics had weak relationships with each other (r < 0.7) (Figure 6). Therefore, the metrics in the combo models were non-collinear and the combo models were robust. The hyperspectral metrics of CRI1 and PCA2 had the weakest relationships with the other metrics. However, the metrics related to crown height were relatively strongly correlated with each other (h₅₀, h₇₅ and HOME_μ) and the waveform metric HOME_μ had a relatively strong positive correlation with the hyperspectral metrics of RVSI.

3.3. Predictive Models with Point Cloud Derived Metrics

The DPC models obtained using point cloud derived metrics are summarized in

Table 5. Summary of the forest structural parameters prediction models with point cloud data based metrics and plot-level accuracy assessment results.

Dependent	Final Models	a Adj-R2	RMSE	rRMSE (%)	k
DBH	18.527 + 0.644h25 − 1.402Wα − 33.340E	0.54 **	5.53	37.02	17.00
HL	1.441 + 0.901h50 + 8.703d9 + 3.206Wα	0.81 **	1.73	15.70	11.26
N	−4699.282 − 71.688h50 + 5796.245Cover + 4437.070E	0.63 **	425.41	30.63	15.77
G	−44.816 + 68.045d1 − 4.286d5 + 12.927O	0.52 **	10.29	40.87	13.90
V	−105.911 + 10.211h50 + 166.787d1 − 229.444d9	0.69 **	33.59	25.76	10.52
AGB	−181.605 + 8.059h95 + 237.778d1 − 158.626hcv	0.68 **	22.45	25.22	9.53

Note: See Table 1 and Table 3 for code of forest structural parameters and point cloud metrics, respectively. ^a Level of significance: ** p < 0.01.

. The point cloud metrics performed well with models significant at p < 0.01. All of the forest structural parameters were well predicted, and the accuracies of the predictive models ranged from 0.52 (Adj-R², rRMSE = 40.87%) to 0.81 (Adj-R², rRMSE = 15.70%). The Lorey’s mean height model had the highest accuracy (Adj-R² = 0.81, rRMSE = 15.70%), followed by the volume model (Adj-R² = 0.69, rRMSE = 25.76%), the above ground biomass model (Adj-R² = 0.68, rRMSE = 25.22%), the stem number model (Adj-R² = 0.63, rRMSE = 30.63%) and the DBH model (Adj-R² = 0.54, rRMSE = 37.02%). The basal area model had the lowest accuracy (Adj-R² = 0.52, rRMSE = 40.87%).

3.4. Predictive Models with Full-Waveform LiDAR Metrics

The FW models obtained using FWF LiDAR metrics are summarized in

Table 6. Summary of the forest structural parameters prediction models with FWF LiDAR metrics (integration of point cloud based metrics and full-waveform based metrics) and plot-level accuracy assessment results.

Dependent	Final Models	a Adj-R2	RMSE	rRMSE (%)	k
DBH	10.925 + 1.051h25 − 31.923E + 20.846VDRμ	0.62 **	3.77	25.24	13.18
HL	2.549 + 0.408h95 + 0.431HOMEμ − 0.007Intμ	0.87 **	1.21	11.01	19.58
N	−5313.482 + 6506.593Cover + 4449.668E − 65.711WDμ	0.70 **	304.40	21.92	9.60
G	−51.975 + 68.144d1 + 11.438Cover − 0.080HOMEσ	0.64 **	7.88	31.30	14.84
V	−28.826 + 44.895d1 + 7.627HOMEμ + 0.177RWEμ	0.80 **	22.58	17.32	14.64
AGB	−446.300 + 12.749h50 + 398.166d1 + 262.543VDRμ	0.82 **	13.72	15.41	15.81

Note: See Table 1 and Table 3 for code of forest structural parameters and FWF LiDAR metrics, respectively. ^a Level of significance: ** p < 0.01.

. The FWF LiDAR metrics performed well with models significant at p < 0.01. All of the forest structural parameters were well predicted, and the accuracies of predictive models ranged from 0.62 (Adj-R², rRMSE = 25.24%) to 0.87 (Adj-R², rRMSE = 11.01%). The Lorey’s mean height model had the highest accuracy (Adj-R² = 0.87, rRMSE = 11.01%), followed by the above ground biomass model (Adj-R² = 0.82, rRMSE = 15.41%), volume model (Adj-R² = 0.80, rRMSE = 17.32%), stem number model (Adj-R² = 0.70, rRMSE = 21.92%) and basal area model (Adj-R² = 0.64, rRMSE = 31.30%); the DBH model had the lowest accuracy (Adj-R² = 0.62, rRMSE = 25.24%). Compared with the models only using DPC, the improvement in models using DPC and FW were significant. The improvements in accuracy ranged from 0.06 to 0.14 (ΔAdj-R², ΔrRMSE = 4.69–11.78%).

3.5. Predictive Models with Combined Full-Waveform LiDAR and Hyperspectral Metrics

The combo models obtained using integrated point cloud, full-waveform and hyperspectral derived metrics are summarized in

Table 7. Summary of the forest structural parameters prediction models with integrated point cloud, full-waveform and hyperspectral metrics and plot-level accuracy assessment results.

Dependent	Final Models	a Adj-R2	RMSE	rRMSE (%)	k
DBH	32.009 − 41.480VDRμ − 0.786NPμ − 1.454CRI1	0.69 **	3.45	23.11	6.93
HL	1.549 + 0.916h75 − 0.047WDσ + 1.533NDVI	0.88 **	1.18	10.68	6.34
N	−66.329 + 2940.122E + 1879.922O − 90635.301RVSI	0.72 **	279.98	20.16	5.18
G	16.252 + 9.265NPμ + 0.050RWEμ − 36.932ARVI	0.68 **	7.22	28.67	4.56
V	−168.279 + 253.243d1 + 6.371HOMEμ − 24.130PCA2	0.81 **	21.34	16.37	5.66
AGB	−113.778 + 11.946h50 + 58.117VDRμ + 124.417NDVI	0.84 **	13.48	15.14	5.59

Note: See Table 1, Table 3 and Table 4 for code of forest structural parameters, FWF LiDAR metrics and hyperspectral metrics, respectively. ^a Level of significance: ** p < 0.01.

. The integration of point cloud, full-waveform and hyperspectral metrics performed well with models significant at p < 0.01. All of the forest structural parameters were well predicted, and the accuracies of the predictive models ranged from 0.68 (Adj-R², rRMSE = 28.67%) to 0.88 (Adj-R², rRMSE = 10.68%). The Lorey’s mean height model had the highest accuracy (Adj-R² = 0.88, rRMSE = 10.68%), followed by the above ground biomass model (Adj-R² = 0.84, rRMSE = 15.14%), volume model (Adj-R² = 0.81, rRMSE = 16.37%), stem number model (Adj-R² = 0.72, rRMSE = 20.16%) and DBH model (Adj-R² = 0.69, rRMSE = 23.11%); the basal area model had the lowest accuracy (Adj-R² = 0.68, rRMSE = 28.67%). Compared with the models only using DPC, the models using integrated DPC, FW and HS significantly improved the accuracies of prediction. The improvements in accuracy ranged from 0.07 to 0.16 (ΔAdj-R², ΔrRMSE = 5.02–13.91%). Moreover, compared with the models using FWF LiDAR metrics, the models using integrated DPC, FW and HS slightly improved the accuracies of prediction. The improvements in accuracy ranged from 0.01 to 0.07 (ΔAdj-R², ΔrRMSE = 0.27–2.63%).

3.6. Assessment of the Predictive Models

Cross validation of the combo models demonstrated that the relationship between field surveyed and predicted forest structural parameters were close to the 1:1 line (Figure 7). The mean differences between field surveyed and cross-validated models predicted forest structural parameters were not significant statistically (

Table 8. The results of cross validation for the forest structural parameters predictive combo models.

Parameters	CV-R2	CV-rRMSE (%)	Differences in Cross-Validation
			a Mean	Std. dev.	Range
DBH	0.65	23.80	0.33 NS	3.56	−6.08–6.22
HL	0.85	11.50	−0.18 NS	1.26	−2.66–2.63
N	0.70	22.90	44.87 NS	317.19	−504.61–554.17
G	0.63	29.04	0.37 NS	7.36	−12.21–12.48
V	0.77	17.95	4.68 NS	23.10	−38.39–45.96
AGB	0.80	17.82	2.05 NS	15.85	−22.92–26.91

Note: See Table 1 for code of forest structural parameters. CV-R²: R² of the result of cross-validation; CV-rRMSE: relative RMSE of the result of cross-validation; ^a Level of significance: NS represents the mean of differences is not significant (p > 0.05).

). The Lorey’s mean height model had the highest accuracy (CV-R² = 0.85, CV-rRMSE = 11.50%), followed by the above ground biomass (CV-R² = 0.80, CV-rRMSE = 17.82%), volume (CV-R² = 0.77, CV-rRMSE = 17.95%), stem number (CV-R² = 0.70, CV-rRMSE = 22.90%) and DBH (CV-R² = 0.65, CV-rRMSE = 23.80%); the basal area model had the lowest accuracy (CV-R² = 0.63, CV-rRMSE = 29.04%). This means that the point cloud, full-waveform and hyperspectral metrics had desirable synergetic effects in the prediction of forest structural parameters.

Figure 8 shows the residuals of six forest structural parameters predicted by models fitted using point cloud metrics, full-waveform LiDAR metrics, and combined full-waveform LiDAR and hyperspectral metrics. It shows that the combo models, which fitted using full-waveform LiDAR and hyperspectral metrics have the smallest variation, followed by the models fitted using full-waveform LiDAR metrics, while the models fitted using point cloud metrics have the largest variation.

4. Discussion

4.1. Synergetic Effects of the Full-Waveform LiDAR and Hyperspectral Data

In this study, the forest structural parameters in subtropical forests were predicted using simultaneously acquired airborne FWF LiDAR and hyperspectral data. The airborne LiDAR data turned out to be the most applicable remote sensing data in forest structural parameters prediction due to its ability to measure the three dimensional structures of forest canopy [16,88]. The hyperspectral data recorded large amounts of continuous narrow bands from VIS to NIR which can provide continuous and detailed spectral signatures of forest biophysical attributes [89,90]. Previous studies have demonstrated that the integration of these two complementary datasets can improve the prediction of forest structural parameters [52,57]. It was also proved in this study, although the FWF LiDAR metrics (only DPC or the integration of DPC and FW) can predict forest structural parameters with relatively high accuracy (Adj-R² = 0.52–0.81, rRMSE = 15.70–40.87%; Adj-R² = 0.62–0.87, rRMSE = 11.01–31.30%), the integration of FWF LiDAR metrics and hyperspectral metrics performs better (Adj-R² = 0.68–0.88, rRMSE = 10.68–28.67%). Zhang et al. [91] predicted six forest structural parameters using only small footprint point cloud LiDAR data in 51 square plots (30 × 30 m²) in the same study site of subtropical forest. Compared with the results in this study, less variabilities in forest structural parameters were explained (45–77%). Dalponte et al. [42] predicted stem diameter and volume using integrated small footprint point cloud LiDAR and hyperspectral data in 52 plots in a temperate forest. The results indicated that the improvements in the prediction of stem diameter and volume using integrated LiDAR and hyperspectral data were 0.4% and 0.5% (rRMSE), respectively. Luo et al. [57] predicted AGB using integrated airborne small footprint point cloud LiDAR and hyperspectral data in 33 plots in a northern temperate deciduous forest. It reported that 2.2% more of the variability in AGB was explained, and 7.9% less predicted error occurred by using integrated point cloud LiDAR and hyperspectral data. The improvement from using integrated point cloud LiDAR and hyperspectral data was slightly higher than this study. This may be caused by the complex structure of subtropical forest, which is multilayered and has greater variability in DBH and tree height, which reduce the ability of the hyperspectral data to predict forest structural parameters.

In the combo models, the five metrics selected from point cloud, full-waveform and hyperspectral metrics had improved synergetic effects for forest structural parameters prediction. Compared with the DPC based models and FWF LiDAR metrics models, the combo models of six forest structural parameters all performed better, explaining large amounts of variability in forest structural parameters (Table 7). Moreover, in Figure 8 the residuals of the combo models were smaller than FWF LiDAR metrics models, and the residuals of DPC-based models were maximal. The distribution of residuals in the combo models were more convergent and the mean of residuals were closer to zero. Consequently, the combo models were more accurate and robust, which indicated that the integration of FWF LiDAR and hyperspectral data has great potential in the prediction of forest structural parameters in subtropical forest.

4.2. Implications of Predictive Models

The parametric regression modeling approach (i.e., backward stepwise regression), which is easy to transplant and can help to understand the relationships between metrics and forest structural parameters, was used to fit predictive models in this study. In Table 5, Table 6 and Table 7, the final models fitted by backward stepwise regression approach are shown. In these tables, the forest structural parameter of H_L had the highest accuracy, and this indicated that full-waveform LiDAR data has unique advantages in the prediction of tree height. Previous studies have demonstrated that the tree height predicted by LiDAR data was similar to field-based measurements [92], and the forest structural parameters related to vegetation height were predicted accurately [93,94]. In this study, it was also confirmed that the outcomes for volume and AGB had a relatively high accuracy.

The full-waveform LiDAR data can provide more information compared with point cloud LiDAR data. The full-waveform metrics, such as HOME and WD, were found to be sensitive to the vertical arrangement of canopy elements and canopy height [95,96]. With the addition of full-waveform metrics, the improvement in the accuracy of predictive models was significant (Table 5 and Table 6). In particular, for G, V and AGB, the improvement of rRMSE is over 10%. Although the hyperspectral data has limitations in quantifying the vertical structure of forest, the hyperspectral metrics (e.g., NDVI, CRI1 and RVSI, etc.) still have the capability to improve the accuracy of forest structural parameters prediction (ΔrRMSE = 0.27–2.63%). This indicates that hyperspectral metrics have great potential in the prediction of forest structural parameters.

4.3. Availability of Full-Waveform LiDAR and Hyperspectral Metrics

Point cloud metrics are significantly related to the three-dimensional structural properties of forest canopy, which allow the prediction of the forest structural parameters. Five point cloud metrics, i.e., the 50th percentile height (h₅₀), 75th percentile height (h₇₅), 10th canopy return density (d₁), euphotic zones of CVM (E), and the oligophotic zones of CVM (O) were selected to fit the combo models (Table 7). h₅₀ and h₇₅ were respectively defined as the heights above ground that are 50% and 75% of the LiDAR returns below the height. h₅₀ likely provides the height of the over-story and h₇₅ presents the spatial variability of the canopy height [25]. Thomas et al. [97] reported that the point cloud metrics of h₅₀ and h₇₅ were strongly related to mean dominated height, basal area, crown closure and AGB. Tsui et al. [25] found that h₅₀ and h₇₅ were the main metrics to explain the majority of the variation in AGB. In the combo models, the h₅₀ and h₇₅ were selected in the predictive models of AGB and H_L, respectively, and the rRMSE of AGB and H_L were 15.14% and 10.68%, respectively. The d₁ was the percentage of points above the 10th quantile to total quantity of points in the plot. This can be used to represent the structure of tree canopy. Previous studies have shown that the d₁ performed well in the prediction of mean dominated tree height and volume [18,53]. In this study, d₁ was selected in the prediction of volume and the rRMSE was 16.37%. E and O were the metrics related to the light extinction in the tree canopy, and are thought to be a three-dimensional analog of cover [69]. They perform well in the prediction of forest structural parameters (e.g., DBH, height, basal area and biomass etc.) [98]. Here, the point cloud metrics of E and O were selected in the prediction of stem number and the predictive model of stem number performed well (rRMSE = 20.16%).

The returned waveforms of the laser system record the specific reflections within the footprint, and are influenced by the materials and vertical structures of canopy. The shape of the returned waveforms can be used to predict forest structural parameters appropriately [7,33]. Sumnall et al. [34] extracted full-waveform metrics (amplitude and echo width variables) from full-waveform data and reported that the full-waveform metrics performed well in the prediction of forest structural parameters (the NRMSE values for the best fit models were in the range of 16–48%). In this study, five waveform metrics, i.e., mean of height of median energy (HOME_μ), mean of vertical distribution ratio (VDR_μ), mean of number of peaks (NP_μ), mean of returned waveform energy (RWE_μ), and the standard deviation of waveform distance (WD_σ), were selected to fit the combo models (Table 7). HOME_μ describes the arrangement of canopy materials in the vertical direction and VDR_μ depicts the change of HOME_μ relative to the height of canopy. These two metrics help to enhance the interpretation of intermediate tree canopies in the mid-story and suppress the trees in the understory. Therefore, the two full-waveform metrics potentially improved the prediction of forest structural parameters. In the combo models, HOME_μ and VDR_μ were selected to predict V, DBH and AGB, and the rRMSE were 16.37%, 23.11% and 15.14%, respectively. NP_μ and RWE_μ were the full-waveform metrics used to describe the number of peak and energy of composited waveform, respectively. NP_μ and RWE_μ were likely related to the arrangement and density of forest materials, which can be used to represent the forest structure. Thus, NP_μ and RWE_μ were selected in the combo predictive models of DBH and G, and the rRMSE were 23.11% and 28.67%, respectively. WD_σ was the full-waveform metric defined as the standard deviation of distances of waveform, which is related to the variation of crown height to some extent. It potentially helped in the prediction of forest structure variables related to tree height. In this study, WD_σ was selected in the combo prediction of H_L and it performed well (rRMSE = 10.68%).

Previous studies have demonstrated that hyperspectral metrics have unique effects in the prediction of forest structural parameters. In this study, the integrated usage of FWF LiDAR and hyperspectral metrics improved the accuracies of forest structural parameters. Compared with the prediction using point cloud metrics only, the improvements in prediction were significant (ΔAdj-R² = 0.07–0.16, ΔrRMSE = 5.02–13.91%). Compared with the prediction using point cloud and full-waveform metrics, the improvements in prediction were distinct (ΔAdj-R²=0.01–0.07, ΔrRMSE = 0.27–2.63%). In this study, five hyperspectral metrics, i.e., normalized difference vegetation index (NDVI), red-edge vegetation stress index (RVSI), atmospherically resistant vegetation index (ARVI), carotenoid reflectance index 1 (CRI1) and the second component of PCA transformation (PCA2), were selected to fit the combo models (Table 7). Many studies have indicated that the NDVI and ARVI are correlated with forest structures such as the leaf area index (LAI) and forest AGB [58,99,100]. Latifi et al. [23] reported that the ARVI was also correlated with mean height. Wang et al. [101] extracted NDVI and ARVI from hyperspectral data and reported that the biomass predictive models performed well using the two metrics (rRMSE = 13.9% and 15.8%). In the combo models, NDVI and ARVI were selected to predict H_L, AGB and G, and the rRMSE was 10.68%, 15.14% and 28.67%, respectively. RVSI was the hyperspectral metric used to identify vegetation stress trends based on the spectral shifting of red-edge. CRI1 is a metric representing reflectance of carotenoid, which was directly correlated with the content of carotenoid. The RVSI and CRI1 potentially improved the prediction of forest structural parameters due to the profound influences of stress and photosynthetic pigments in the forest structure formation. The RVSI and CRI1 were selected in the combo models of N and DBH, and the rRMSE were 20.16% and 23.11%, respectively. PCA analysis was widely used to reduce the dimensions of the dataset. In general, the first several components contain the most of variance in the dataset. In this study, PCA2 contained much more variances of hyperspectral reflectance channels and might d perform well in the prediction of forest structural parameters. In the combo models, PCA2 was selected to predict V and the rRMSE of prediction was 16.37%.

The full-waveform LiDAR metrics (including point cloud and full-waveform metrics) have the ability to record the forest structure in detail. In Figure 4, the vertical distribution of point cloud and pulse energy indicate the distribution of forest materials and the transmission of energy in the canopy. The hyperspectral metrics are associated with the biophysical and biochemical attributes of forest, which can be used to reflect the health conditions and growth stage of forest vegetation (e.g., Figure 4 (III)). Therefore, the full-waveform LiDAR metrics, hyperspectral metrics and the forest structural parameters predicted in this study can be used to assess the structural, compositional and functional properties of forest ecosystems. Moreover, the forest nutritional status, wildlife habitat, and ecosystem service value can also be evaluated using the full-waveform LiDAR and hyperspectral data.

4.4. Potential Improvements and Indications for Future Research

This study combined full-waveform LiDAR and hyperspectral data in the feature level, and the addition of hyperspectral metrics in full-waveform LiDAR metrics improved the accuracy of forest structural parameters. In the future, more advanced approaches associated with physical theory or the mechanisms of sensors will be used to explore the better integration of full-waveform LiDAR and hyperspectral data. Moreover, the non-parametric modeling approach can be used to enhance the capability of prediction of forest structural parameters in subtropical forests.

5. Conclusions

Simultaneously acquired airborne full-waveform (FWF) LiDAR and hyperspectral data were used to predict forest structural parameters in subtropical forests of southeast China. The pulse amplitude and waveform shape of airborne full-waveform LiDAR data (FWF) were calibrated using a physical process-driven and a voxel-based approach, respectively. Different suites of LiDAR and hyperspectral metrics were calculated and selected using correlation analysis and principal component analysis (PCA). The selected point cloud, full-waveform and hyperspectral metrics were used to fit regression models individually and in integration to predict six forest structural parameters, and the capability of predictive models and synergetic effects of metrics were assessed using leave-one-out cross validation. The results showed that: most of the metrics selected from three groups divided by the PCA analysis, which were highly correlated with the first and second principal component, had a weak relationship with each other (r < 0.7); the prediction of H_L using the metrics of h₇₅, WD_σ and NDVI had a relatively higher accuracy (Adjusted-R² = 0.88, relative RMSE = 10.68%) than the other forest structural parameters; and the usage of DPC and FW resulted in higher accuracies (Adjusted-R² = 0.62–0.87, relative RMSE = 11.01–31.30%) than the models only including DPC (Adjusted-R² = 0.52–0.81, relative RMSE = 15.70–40.87%). Moreover, the integration of DPC, FW and HS had a positive synergetic effect for forest structural parameters prediction (Adjusted-R² = 0.68–0.88, relative RMSE = 10.68–28.67%).