Estimating Aboveground Carbon Stock at the Scale of Individual Trees in Subtropical Forests Using UAV LiDAR and Hyperspectral Data

: Accurate estimation of aboveground carbon stock for individual trees is important for evaluating forest carbon sequestration potential and maintaining ecosystem carbon balance. Airborne light detection and ranging (LiDAR) data has been widely used to estimate tree-level carbon stock. However, few studies have explored the potential of combining LiDAR and hyperspectral data to estimate tree-level carbon stock. The objective of this study is to explore the potential of integrating unmanned aerial vehicle (UAV) LiDAR with hyperspectral data for tree-level aboveground carbon stock estimation. To achieve this goal, we ﬁrst delineated individual trees by a CHM-based watershed segmentation algorithm. We then extracted structural and spectral features from UAV LiDAR and hyperspectral data respectively. Then, Pearson correlation analysis was conducted to assess the correlation between LiDAR features, hyperspectral features, and tree-level carbon stock, based on which, features were selected for model development. Finally, we developed tree-level carbon stock estimation models based on the Schumacher–Hall formula and stepwise multiple regression. Results showed that both LiDAR and hyperspectral features were strongly correlated to tree-level carbon stock. Both tree height ( H , r = 0.75) and Green index ( GI , r = 0.83) had the highest correlation coefﬁcients with tree-level carbon stock in LiDAR and hyperspectral features, respectively. The best model using LiDAR features alone includes the metrics of H , the 10th height percentile of points ( PH10 ), and mean height of points ( H mean ), and can explain 74% of the variations in tree-level carbon stock. Similarly, the best model using hyperspectral data includes GI and modiﬁed normalized differential vegetation index ( mNDVI ), and has similar explanatory power (r 2 = 0.75). The model that integrates predictors, namely, GI and the 95th height percentile of points ( PH95 ) from hyperspectral and LiDAR data, substantially improves the explanatory power (r 2 = 0.89). These results indicated that while either LiDAR data or hyperspectral data alone can estimate tree-level carbon stock with reasonable accuracy, combining LiDAR and hyperspectral features can substantially improve the explanatory power of the model. Such results suggested that tree-level carbon stock estimation can greatly beneﬁt from the complementary nature of LiDAR-detected structural characteristics and hyperspectral-captured spectral information of vegetation.


Introduction
With the increase in global warming and Greenhouse Effect, carbon cycle has become a hot topic in global climate change research [1][2][3][4][5][6][7]. Forests have the ability to absorb and fix CO 2 from the atmosphere and are the largest carbon pool in terrestrial ecosystems [8,9]. The United Nations Intergovernmental Panel on Climate Change (IPCC) has The overall objective of this study is to investigate whether the combination of UAV LiDAR and hyperspectral data can improve the estimation of aboveground carbon stocks at the scale of individual trees. Specifically, we first evaluate the ability of LiDAR structural features and hyperspectral features in estimating aboveground carbon stocks of individual trees separately. We then examine whether combining UAV LiDAR and hyperspectral data can improve the accuracy of tree-level aboveground carbon stock estimation.

Study Area
The study area is located on Julongshan Park, a subtropical broadleaf forest situated in Shenzhen in southern China (Figure 1). The center coordinates of the study region are 114°23′28″E and 22°43′50″N. The topography of the study region is plain, with a mean altitude of 64 m. The climate is subtropical maritime, with an annual average temperature of 22.4 °C and an annual rainfall of 1933.3 mm. There are several tree species in the study area, such as Terminalia mantaly, Elaeocarpus grandifloras, Litchi chinensis, Lagerstroemia speciosa, Acacia mangium, and Cinnamomum camphora. Due to the limitations of field accessibility and sample number, this research was only conducted on Terminalia mantaly which accounted for about 50% of total individuals.

LiDAR Data
The LiDAR data were acquired on August 2019, by a Velodyne LiDAR PUCK-16 laser mounted on the DJI Matrice 600 Pro six-rotor aircraft (DJI, Shenzhen, China). The laser emitted pulse at a frequency of 5 kHz with a wavelength of 903 nm. The flight altitude

Remote Sensing Data 2.2.1. LiDAR Data
The LiDAR data were acquired on August 2019, by a Velodyne LiDAR PUCK-16 laser mounted on the DJI Matrice 600 Pro six-rotor aircraft (DJI, Shenzhen, China). The laser emitted pulse at a frequency of 5 kHz with a wavelength of 903 nm. The flight altitude was about 70 m above ground and the average flying speed was 3.6 m/s. Two returns were recorded for each pulse, and the average point cloud density was larger than 100 points/m 2 . The point clouds were first filtered into ground points and nonground points by the TerraSolid software (Terrasolid, Helsinki, Finland). The density of the final point cloud containing only ground-level observations was 8 points/m 2 . TerraSolid software uses an adaptive TIN algorithm for filtering, which was first proposed by Axelsson et al. [39] Next, digital terrain model (DTM) and digital surface model (DSM) were produced based on ground points and all points, respectively. Then, a canopy height model (CHM), representing the height of the vegetation above ground surface, was calculated by subtracting DTM values from DSM values. Finally, z-values of points were normalized by subtracting DTM values from the initial, absolute, z-values [40].

Hyperspectral Data
UAV hyperspectral imagery were acquired on August 2019 by the DJI Matrice 600 Pro UAV with a Resonon Pika-L sensor, which is a push-broom hyperspectral sensor. The hyperspectral images had 150 spectral bands between 400 nm and 1000 nm, with a spectral resolution of 2.2 nm. The flight altitude was about 200 m above ground, obtaining the hyperspectral images with a spatial resolution of 0.3 m. The hyperspectral images were first radiometrically corrected according to reference data, then geometrical corrections were carried out based on LiDAR-derived DTM and ground-based GPS data. Finally, atmospheric corrections were applied using the FLAASH module of ENVI software.

Field Data
A fieldwork campaign was carried out in August 2019. A total number of 60 Terminalia mantaly trees were measured. The tree height (H) and DBH were measured using a laser hypsometer and a tape measure respectively, and the position of each tree was acquired using a real-time kinematic GPS (RTK-GPS). Based on the carbon coefficient (CC) and aboveground biomass (AGB), tree-level carbon stock (CS) was calculated to serve as control data (Equation (1)) [15]. After an in-depth literature analysis, no available carbon coefficient and allometric growth equation for Terminalia mantaly have yet been discovered. Thus, since this species is a subtropical deciduous broad-leaved tree species in China, we decided instead to use the general carbon coefficient and allometric growth equation of subtropical deciduous broad-leaved forest from China. The allometric growth equation and carbon coefficient were obtained by cutting down trees, then drying, burning, and weighing them, and finally, we performed a regression analysis. In this study, the carbon coefficient of Terminalia mantaly was set to 0.4956 according to Li et al. [41] The aboveground biomass of each individual tree was calculated as the sum of different forest components, including stem biomass (stemAGB), branch biomass (branchAGB), and leaf biomass (leafAGB) (Equation (2)) [42]. The biomass of stem, branch, and leaf were calculated based on their DBH and height according to the allometric growth equation as shown in Equations (3)-(5), whose reported r 2 values are 0.98, 0.97, and 0.96, respectively [42].

Methods
Tree-level carbon stocks were estimated by UAV LiDAR and hyperspectral data. The procedure consisted of three main steps: (1) individual tree segmentation; (2) treelevel structure and spectral features extraction; and (3) tree-level carbon stocks estimation ( Figure 2).

Methods
Tree-level carbon stocks were estimated by UAV LiDAR and hyperspectral data. The procedure consisted of three main steps: (1) individual tree segmentation; (2) tree-level structure and spectral features extraction; and (3) tree-level carbon stocks estimation (Figure 2).

Individual Tree Segmentation
In this study, a CHM-based watershed segmentation algorithm was used to delineate tree crowns, which had been successfully used for individual tree segmentation in many studies [43][44][45][46][47][48]. The procedure of CHM-based watershed segmentation ( Figure 3) is presented in detail by Chen et al. [44]. Initially, we subtracted CHM from its maximum value to obtain the inverted CHM. Next, a local minimum filter algorithm was used to find local minimum points from the inverted CHM. Then, these local minimum points were used as markers to simulate the immersion process. The influence area of each local minimum point expanded outward gradually, and "dams" were built at the confluence of multiple "basins", which formed a "watershed". Finally, the CHM was segmented to obtain individual trees.

Individual Tree Segmentation
In this study, a CHM-based watershed segmentation algorithm was used to delineate tree crowns, which had been successfully used for individual tree segmentation in many studies [43][44][45][46][47][48]. The procedure of CHM-based watershed segmentation ( Figure 3) is presented in detail by Chen et al. [44]. Initially, we subtracted CHM from its maximum value to obtain the inverted CHM. Next, a local minimum filter algorithm was used to find local minimum points from the inverted CHM. Then, these local minimum points were used as markers to simulate the immersion process. The influence area of each local minimum point expanded outward gradually, and "dams" were built at the confluence of multiple "basins", which formed a "watershed". Finally, the CHM was segmented to obtain individual trees.

LiDAR Features Extraction
To estimate tree-level carbon stocks, we extracted 10 structural features from UAV LiDAR data, including tree height, crown width, and the distribution characteristics of point clouds, which have been commonly used in forest aboveground biomass and carbon stocks estimation [26,[49][50][51]. LiDAR features and their formulas are shown in Table 1. LiDAR features extraction was implemented by programming with Matlab software.

LiDAR Features Extraction
To estimate tree-level carbon stocks, we extracted 10 structural features from UAV LiDAR data, including tree height, crown width, and the distribution characteristics of Remote Sens. 2021, 13, 4969 6 of 15 point clouds, which have been commonly used in forest aboveground biomass and carbon stocks estimation [26,[49][50][51]. LiDAR features and their formulas are shown in Table 1. LiDAR features extraction was implemented by programming with Matlab software.

Features Formula/Description
Tree height The 10th height percentile of point clouds PH25 The 25th height percentile of point clouds PH50 The 50th height percentile of point clouds PH75 The 75th height percentile of point clouds PH90 The 90th height percentile of point clouds PH95 The 95th height percentile of point clouds Where CD S_N and CD E_W are the crown diameters in south-north direction and eastwest direction, respectively. h i is the height of the ith point, and n is the number of LiDAR points.

Hyperspectral Features Extraction
Nineteen tree-level narrowband spectral features, which are commonly used in forest carbon stocks or biomass estimation, were extracted from hyperspectral imagery, including principal components, band reflectance, and vegetation indices, to estimate carbon stock. The first three principal components were extracted by principal component analysis (PCA) algorithm. The mean vegetation indices and band reflectance of each tree were calculated based on hyperspectral data. Hyperspectral features and their formulas are shown in Table  2. Hyperspectral features extraction was implemented using ENVI software. Where ρ wavelength means the reflectance in the wavelength. For example, ρ 798 means the reflectance in the wavelength of 798 nm. n means the number of bands.

Tree-Level Carbon Stock Estimation
Pearson correlation analysis was carried out to explore the relationship between LiDAR features, hyperspectral features, and control data. Pearson correlation analysis was implemented using IBM SPSS Statistics 22 software. The LiDAR and hyperspectral features having a significant correlation with control data (with an absolute value of correlation coefficient (r) greater than 0.5) were selected for subsequent carbon stocks estimation of individual trees. The tree-level carbon stock estimation model has the form of the Schumacher-Hall formula as expressed in Equation (6), where LI i refers to the ith LiDARderived structural feature, SI j refers to the jth narrowband hyperspectral feature, c, a i , and b j are the regression coefficients. In order to facilitate modeling, the natural logarithm transformation of model (6) was carried out to obtain the model form (Equation (7)).
To investigate which combination of LiDAR and hyperspectral features can best explain the variance of tree-level carbon stocks, a stepwise multiple regression analysis was carried out. Stepwise multiple regression was implemented using IBM SPSS Statistics 22 software. Stepwise multiple regression is a commonly-used and effective method to eliminate the multicollinearity issue and select the optimal regression equation, which has been widely used to estimate forest carbon stocks or biomass [65,66]. To estimate carbon stocks of individual trees accurately, stepwise multiple regression models, which included LiDAR metrics and hyperspectral metrics as predictors, were developed independently and in combination.

Accuracy Assessment
To assess the predictive performance of tree-level carbon stock estimation models, five accuracy assessment measures were applied. The coefficient of determination (r 2 ) (Equation (8)), root mean square error (RMSE) (Equation (9)), and mean absolute error (MAE) (Equation (10)) have been widely used to evaluate the prediction accuracy [26]. Percentage RMSE (PRMSE) is RMSE expressed as a percentage of the mean value (Equation (11)). Root mean square percentage error (RMSPE) is an average deviation from a true value (Equation (12)). PRMSE and RMSPE are all scale-independent metrics to assess the predictive performance of a model. In Equations (8)- (12), n is the number of samples, y i and y i are the true and predicted values of the ith sample, and y is the average true value of all samples. In this study, MAE, RMSE, PRMSE, and RMSPE were used to evaluate the prediction accuracy of tree-level carbon stock estimation models, and the model with a smaller metric value was considered to have a higher tree-level carbon stock accuracy. The overall prediction performance (OPP) of tree-level carbon stock estimation model can be calculated according to Equation (13). The model's uncertainty can be calculated according to Equation (14). Lower values for Uncertainty increase the reliability of predictions. These accuracy metrics were calculated by programming with Matlab software.

Results
Pearson correlation analysis was carried out to explore the relationship between each remote sensing feature and control data. The correlation coefficients between LiDAR features, hyperspectral features, and control data are plotted in Figure 4. Nine of the ten LiDAR features, including crown width, tree height, H mean , the height percentiles of point clouds, had positive correlations with control data. All of the above had correlation coefficients larger than 0.48. Only one LiDAR feature, i.e., H std , had no correlation with control data. For all the LiDAR features, tree height had the highest correlation with control data (r = 0.75), followed by PH95 (r = 0.70) and PH90 (r = 0.68). Among the 19 hyperspectral features, six features (PCA3, Datt, GI, mNDVI, VOG, and MRESRI) had positive correlations greater than 0.5 with control data, and four features (PCA2, ACI, PSRI, and SIPI) had negative correlations lower than -0.3 with control data. For all hyperspectral features, GI had the highest correlation with control data (r = 0.83), followed by VOG (r = 0.75) and MRESRI (r = 0.75).
As shown in Table 3, the tree-level carbon stock models derived from LiDAR metrics alone, hyperspectral metrics alone, and the combination of them were coded with "L", "H", and "LH", respectively. These models were derived by multiple stepwise regression. Among the ten LiDAR metrics, tree height, 10th height percentile of point clouds, and mean point height were selected by stepwise regression to develop the tree-level carbon stock estimation model. From all hyperspectral metrics, GI and mNDVI were selected to establish the tree-level carbon stock estimation model. When combining LiDAR and hyperspectral metrics, the tree-level carbon stock estimation model was established based on GI and 10th height percentile of point clouds. data. For all the LiDAR features, tree height had the highest correlation with control data (r = 0.75), followed by PH95 (r = 0.70) and PH90 (r = 0.68). Among the 19 hyperspectral features, six features (PCA3, Datt, GI, mNDVI, VOG, and MRESRI) had positive correlations greater than 0.5 with control data, and four features (PCA2, ACI, PSRI, and SIPI) had negative correlations lower than -0.3 with control data. For all hyperspectral features, GI had the highest correlation with control data (r = 0.83), followed by VOG (r = 0.75) and MRESRI (r = 0.75).  The coefficient of determination (r 2 ) of those tree-level carbon stock estimation models ranged from 0.66 to 0.89. Among the three kinds of tree-level carbon stock models, the LiDAR-derived model had the smallest r 2 , where the r 2 values of model L1, L2, and L3 were 0.66, 0.70, and 0.74, respectively. The hyperspectral-derived model had a higher r 2 than the LiDAR-derived model, in which the r 2 values of model H1 and H2 were 0.72 and 0.75, respectively. The combined model (GI and PH95 are the predictors) had the highest r 2 (0.89).
Beyond the coefficient of determination, four accuracy statistical indices, MAE, RMSE, PRMSE, and RMSPE, were calculated to evaluate the bias of tree-level carbon stock estimation models. As shown in Figure 5, the MAE of the four models ranged from 2.36 kg to 5.86 kg with an average value of 4.55 kg. The RMSE ranged from 3.98 kg to 8.21 kg with an average of 6.56 kg. The PRMSE ranged from 14% to 30% with an average of 23%. The RMSPE ranged from 15% to 30% with an average of 24%. Among the three kinds of tree-level optimal carbon stock estimation models, the combined model (LH2)   As shown in Figure 6, the control data ranged from 5.98 kg to 56.28 kg with an average of 27.71 ± 11.30 kg. LiDAR-derived models had slightly smaller mean and standard deviation than control data, and the maximum value and the minimum value were all higher than control data. Hyperspectral-derived models also had slightly smaller mean and standard deviation, and their ranges of maximum and minimum values were significantly smaller than control data. The minimum value, maximum value, range, and average value of the combined model were all closest to control data. Therefore, the combined model had the optimal tree-level carbon stock estimation results. The overall prediction performance and uncertainty of all models are shown in Table  4. The values of OPP of all models ranged from 70 to 85. The uncertainties of all models ranged from −5 to −1. The OPP and uncertainty of these models are similar, so these models can account for variations in tree-level carbon stocks. As shown in Figure 6, the control data ranged from 5.98 kg to 56.28 kg with an average of 27.71 ± 11.30 kg. LiDAR-derived models had slightly smaller mean and standard deviation than control data, and the maximum value and the minimum value were all higher than control data. Hyperspectral-derived models also had slightly smaller mean and standard deviation, and their ranges of maximum and minimum values were significantly smaller than control data. The minimum value, maximum value, range, and average value of the combined model were all closest to control data. Therefore, the combined model had the optimal tree-level carbon stock estimation results. As shown in Figure 6, the control data ranged from 5.98 kg to 56.28 kg with an average of 27.71 ± 11.30 kg. LiDAR-derived models had slightly smaller mean and standard deviation than control data, and the maximum value and the minimum value were all higher than control data. Hyperspectral-derived models also had slightly smaller mean and standard deviation, and their ranges of maximum and minimum values were significantly smaller than control data. The minimum value, maximum value, range, and average value of the combined model were all closest to control data. Therefore, the combined model had the optimal tree-level carbon stock estimation results. The overall prediction performance and uncertainty of all models are shown in Table  4. The values of OPP of all models ranged from 70 to 85. The uncertainties of all models ranged from −5 to −1. The OPP and uncertainty of these models are similar, so these models can account for variations in tree-level carbon stocks. The overall prediction performance and uncertainty of all models are shown in Table 4. The values of OPP of all models ranged from 70 to 85. The uncertainties of all models ranged from −5 to −1. The OPP and uncertainty of these models are similar, so these models can account for variations in tree-level carbon stocks.

Discussion
Among all LiDAR-derived structural features, tree height showed the highest correlation with tree-level carbon stocks. This is an indication that tree height is the dominant factor of carbon stock estimation, which is consistent with previous studies [67,68]. The point height distribution features, including PH95, PH90, H mean , PH75, PH50, PH25, and PH10, are also correlated well with carbon stock. This is because these features can describe the vertical structure of forest, which have an influence on carbon stock [68]. Crown diameter had weak correlation with carbon stock (r = 0.48), and no correlation was found between H std and carbon stock. Additionally, stepwise multiple regression results indicated that the combination of tree height, PH10, and H mean could provide accurate carbon stock estimation of individual trees. This may be owing to the fact that carbon stock is mainly determined by tree height and vertical structure [68,69].
Previous studies have found that hyperspectral metrics can be used to estimate forest aboveground biomass at plot level [36,70,71]. Luo et al. estimated plot-level forest biomass using hyperspectral-derived narrow-band vegetation indices, and the estimated biomass showed a good correlation with field biomass [72]. Gong et al. used airborne hyperspectral imagery to estimate forest biomass, and they found vegetation indices and red edge positions were effective for forest biomass retrieval [73]. This can be explained by the fact that hyperspectral data could describe the forest growth status, which is related to biomass [36,70]. However, our approach, in contrast to previously cited research, was to explore the ability of hyperspectral features in estimating tree-level carbon stock. For hyperspectral metrics, GI was most relevant to tree-level carbon stock as it explained 72% of the variance of carbon stocks. The absolute correlation coefficient between calculated carbon stock and PCA2, PCA3, Datt, PSRI, VOG, and MRESRI were all greater than 0.5. These results indicated that hyperspectral features had potential in estimating tree-level carbon stock.
Stepwise multiple regression showed that the combined model provided the highest tree-level carbon stock estimation accuracy. In addition, the predictions of the combined model were closest to control data. Therefore, combining UAV LiDAR and hyperspectral features could improve the accuracy of tree-level carbon stock estimation than using LiDAR and hyperspectral data alone. The improvement might be due to the complementarity between LiDAR-detected structural characteristics and hyperspectral-captured vegetation spectral information, and the combined model took full advantage of the predictive ability of structural features and hyperspectral features, which is consistent with previous studies [72,74,75]. Several studies have combined hyperspectral and LiDAR data to estimate carbon stocks at individual tree level [34,76,77]. They used hyperspectral information to classify tree species, and then developed carbon stock estimation model based on LiDARderived structural features. Kandare et al. [38] classified tree species using hyperspectral data and estimated broadleaf tree-level volume using LiDAR-derived tree height and DBH, resulting in an estimation accuracy of 58%. Alonzo et al. [74] estimated tree-level carbon stocks of a broadleaf tree species (e.g., Quercus agrifolia) by 28 LiDAR-derived structural metrics, with an accuracy of 83%. In this study, hyperspectral features and LiDAR structural features were combined to develop tree-level carbon stock estimation model, and they explained 89% of the variance in tree-level carbon stock. The higher estimation accuracy of tree-level carbon stocks reported in our study owes to the introduction of hyperspectral features. This can be explained by the fact that carbon stock is not only related to tree-level structural features that can be extracted from LiDAR data, but also to biomass conversion factor and carbon coefficient which can be reflected in hyperspectral information [37].

Conclusions
In this study, we explored the ability of combining UAV LiDAR and hyperspectral data in estimating tree-level carbon stocks in subtropical forests. We extracted tree-level structural and spectral features from UAV LiDAR and hyperspectral data, respectively, and then evaluated the capacities of LiDAR and hyperspectral features alone and in combination to predict aboveground carbon stock of individual trees. Results indicated that LiDAR or hyperspectral features alone can yield reasonable carbon stock (r 2 = 0.74 and 0.75, respectively). As expected, the combination of LiDAR and hyperspectral data could improve the accuracy of tree-level carbon stock estimation to some degree (r 2 = 0.89), and LiDAR-derived PH95 and hyperspectral-derived GI were included in the combination model. The improvement can be attributed to the fact that carbon stock is not only related to tree-level structural features that can be extracted from LiDAR data, but also to biomass conversion factor and carbon coefficient which can be reflected in hyperspectral information. Therefore, if both LiDAR and hyperspectral data are available, the fusion of LiDAR and hyperspectral data is the best method to accurately estimate tree-level aboveground carbon stock. This study could be a valuable resource for researchers and forest managers to obtain more accurate tree-level carbon stocks. Further efforts should focus on improving this approach or applying it to more forest types that have specific allometric equations, as well as in larger areas. In addition, based on the method of tree-level carbon stock estimation, multi-temporal remote sensing data should be obtained and then used to estimate forest carbon sequestration to ameliorate the unavoidable impacts of climate change.