Comparison of Coniferous Plantation Heights Using Unmanned Aerial Vehicle (UAV) Laser Scanning and Stereo Photogrammetry

: Plantation forests play a critical role in forest products and ecosystems. Unmanned aerial vehicle (UAV) remote sensing has become a promising technology in forest related applications. The stand heights will reﬂect the growth and competition of individual trees in plantation. UAV laser scanning (ULS) and UAV stereo photogrammetry (USP) can both be used to estimate stand heights using different algorithms. Thus, this study aimed to deeply explore the variations of four kinds of stand heights including mean height, Lorey’s height, dominated height, and median height of coniferous plantations using different models based on ULS and USP data. In addition, the impacts of thinned point density of 30 pts to 10 pts, 5 pts, 1 pts, and 0.8 pts/m 2 were also analyzed. Forest stand heights were estimated from ULS and USP data metrics by linear regression and the prediction accuracy was assessed by 10-fold cross validation. The results showed that the prediction accuracy of the stand heights using metrics from USP was basically as good as that of ULS. Lorey’s height had the highest prediction accuracy, followed by dominated height, mean height, and median height. The correlation between height percentiles metrics from ULS and USP increased with the increased height. Different stand heights had their corresponding best height percentiles as variables based on stand height characteristics. Furthermore, canopy height model (CHM)-based metrics performed slightly better than normalized point cloud (NPC)-based metrics. The USP was not able to extract exact terrain information in a continuous coniferous plantation for forest canopy cover (CC) over 0.49. The combination of USP and terrain from ULS can be used to estimate forest stand heights with high accuracy. In addition, the estimation accuracy of each forest stand height was slightly affected by point density, which can also be ignored.


Introduction
Plantation forests, accounting for 7.3% (290 million ha) of the world's forest cover (3999 million ha) [1], play a critical role in the restoration and reconstruction of forest ecosystems, the provision of forest products, the increase in forest carbon sink, the improvement of ecological environment, meeting demands for wood and other forest products [2,3]. The covering area of planted forests in China is the largest worldwide and is about approximately 36% of the total national forest coverage [4]. As is widely-known, regardless of the purpose of plantation forest management, ecological benefits, or economic benefits, plantation forests are inseparable from forest management decision-making. Forest structural attributes and biophysical properties reflect the growth and health of the forests to a certain extent. Therefore, efficient and accurate estimations of the properties of plantation forests are critically important for forest managers to quantitatively understand the forests' condition [5].
Forest height is an important parameter required for forest management, inventory, and monitoring [6]. Additionally, it is also an ecological measure that provides essential information related to ecological, hydrological, biophysical, and micro-meteorological properties [7]. As such, the estimation accuracy and variations of forest height is sensitive to the estimation accuracy of above ground biomass (AGB) [8], which plays a critical role in programs aimed at reducing the global emissions of carbon from deforestation and the degradation of forests [9]. Forest canopy height is defined as the height of the upper envelope of tree crown above terrain, while forest stand height is defined as the statistic height of tree height in a forest stand [10].
Light Detection and Ranging (LiDAR) has the ability to penetrate forest canopy to acquire 3-dimentional information of the observation target based on the time difference between the emission and return of the laser pulses [11]. In recent years, with the development of photogrammetry algorithms, computer technology and stereo photogrammetry images also have the ability to generate dense point clouds with 3-dimentional characteristics. Subsequently, the exploitation of unmanned aerial vehicle (UAV) systems promotes the flexibility of forest attribute estimation at small-scale [12] or large-scale forest inventory [13], which mainly include the two technologies, UAV laser scanning (ULS) and UAV stereo photogrammetry (USP). The common characteristic of the two technologies is its corresponding capability of generating 3D point clouds data, which have great suitability for forest attribute estimation [14,15]. Some comparisons considering the capacity of USP and ALS on forests have been done to analyze the suitability of USP to be an alternative of ULS [16,17]. A variety of comparable studies on the performance of LiDAR and USP point clouds for forest attribute parameters have been conducted [16] including diverse forest environments such as complex boreal forest types, complex coastal forest [18], and subtropical planted forests [19]. Most studies have concluded that LiDAR provides higher accuracy for modeling and predicting forest attributes than stereo photogrammetry data [20][21][22]. However, few studies have analyzed the performance of USP on temperate coniferous plantations in mountainous area.
In contrast to LiDAR, dense point clouds generated from photogrammetry images are limited to the canopy envelope of tree crown above terrain. The vertical distributions of NPC or CHM from ULS and USP were more similar with the increasing heights [23]. Metrics mainly including height metrics and density metrics representing the vertical distribution of canopy having strong correlation with forest structural attributes [24,25], especially with forest height. In operational forest inventories, an area-based approach (ABA) is more widely used to estimate forest stand heights and other forest attributes using metrics calculated from point clouds or CHM by linear regression [12], non-linear regressions [16,20], and machine learning methods [26,27]. Noordermeer et al. [20] estimated Lorey's height and dominated height in five districts dominated by coniferous forests in southeastern Norway using height metrics and density metrics from ALS point clouds and aerial photogrammetry by nonlinear prediction models. Mean height and median height were successfully estimated by UAV photogrammetry using simple linear regression with high accuracy [28]. Järnstedt et al. [27] compared the estimation accuracy of mean and dominant height using aerial imagery and ALS data. Puliti et al. [12] compared the performance of USP and ULS data to estimate Lorey's height and dominated height. Tompalski et al. [29] assessed the capability of ALS and aerial photogrammetry data to estimate Lorey's height. What we can observe from previous publications is that the commonly used observed forest stand heights are mean height, Lorey's height, dominated height, and median height, but they have been rarely systematically and comprehensively analyzed in one study using the metrics from ULS or USP, respectively. More study needs to be conducted to deeply understand the questions related with forest heights estimated using ULS and USP. For example, are the estimation accuracies of forest heights using a low-cost consumer-grade camera similar to those using ULS in temperate coniferous Remote Sens. 2021, 13, 2885 3 of 23 plantations taking into account the studied forest type in the current literature? What kind of forest height such as Lorey's height, dominated height, arithmetic mean height, and median height has a better performance using ULS and USP? Additionally, the difference of flight parameters between LiDAR and stereo photogrammetry flight task and the point cloud density of ULS and USP point clouds may be different. While point density is one of important factors affecting the estimating accuracy of forest structure attributes [30][31][32][33], it may be more persuasive for comparisons between ULS and USP point clouds with similar point cloud density. However, previous publications on the comparison between the ULS and USP have rarely considered the effect of point cloud density difference on the performance of the ULS and USP on forest structural estimation.
In this paper, we selected a typical coniferous plantation in a mountainous area of northern China to testify the potential of ULS and USP point clouds to estimate different forest heights. The objectives of this study were to (1) compare the capability of ULS and USP data to estimate forest stand heights; (2) identify the optimal forest stand height model among the four commonly used forest stand heights; and (3) evaluate the effects of point density on the estimation of forest heights. These will be very useful for forest inventory and management.

Study Site
The study site is located at the Wangyedian forest farm (118 • 9 ~118 • 30 E, 41 • 21 ~41 • 39 N) in Inner Mongolia Autonomous Province, China, and covers approximately 5.52 × 104 hm 2 ( Figure 1). The elevation ranges from 800 m to 1890 m above sea level and it belongs to a semi-arid warm temperate continental climate. The annual temperature changes from −31 • to 36 • C with an annual average temperature of 4.2 • C. The annual precipitation is about 400~600 mm, which mainly occurs in July and August. The dominate tree species of the plantation are L. principis-rupprechtii (Larix principis-rupprechtii) and P. tabuliformis (Pinus tabuliformis) while other areas are covered by Betula platyphylla and Betula dahurica Pallas. The plantation stage ranges from young to near-mature forest.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 23 low-cost consumer-grade camera similar to those using ULS in temperate coniferous plantations taking into account the studied forest type in the current literature? What kind of forest height such as Lorey's height, dominated height, arithmetic mean height, and median height has a better performance using ULS and USP? Additionally, the difference of flight parameters between LiDAR and stereo photogrammetry flight task and the point cloud density of ULS and USP point clouds may be different. While point density is one of important factors affecting the estimating accuracy of forest structure attributes [30][31][32][33], it may be more persuasive for comparisons between ULS and USP point clouds with similar point cloud density. However, previous publications on the comparison between the ULS and USP have rarely considered the effect of point cloud density difference on the performance of the ULS and USP on forest structural estimation.
In this paper, we selected a typical coniferous plantation in a mountainous area of northern China to testify the potential of ULS and USP point clouds to estimate different forest heights. The objectives of this study were to (1) compare the capability of ULS and USP data to estimate forest stand heights; (2) identify the optimal forest stand height model among the four commonly used forest stand heights; and (3) evaluate the effects of point density on the estimation of forest heights. These will be very useful for forest inventory and management.

Study Site
The study site is located at the Wangyedian forest farm (118°9′~118°30′ E, 41°21′~41°39′ N) in Inner Mongolia Autonomous Province, China, and covers approximately 5.52 × 104 hm 2 ( Figure 1). The elevation ranges from 800 m to 1890 m above sea level and it belongs to a semi-arid warm temperate continental climate. The annual temperature changes from −31° to 36 °C with an annual average temperature of 4.2 °C. The annual precipitation is about 400 mm~600 mm, which mainly occurs in July and August. The dominate tree species of the plantation are L. principis-rupprechtii (Larix principis-rupprechtii) and P. tabuliformis (Pinus tabuliformis) while other areas are covered by Betula platyphylla and Betula dahurica Pallas. The plantation stage ranges from young to near-mature forest.

Field Data
There are 19 blocks regularly distributed in the study area, and four plots of 25 m × 25 m were designed in each block (Figure 1b). The diameter at breast heights (DBHs), species, tree heights (m), and crown widths (m) at perpendicular directions of north-south and east-west of individual trees with DBH ≥ 5 cm were measured in September 2019. The number of plots dominated with L. principis-rupprechtii (denoted as LYS) and P. tabuliformis (denoted as YS) were 30 and 41, respectively, and used as references to assess the accuracy of the estimated stand heights.
The heights of individual trees of all plots varied from 2 m to 24.1 m with a mean of 11.23 m. The mean of tree heights of the LYS and YS plots were 11.67 m and 11.05 m, respectively. The histogram of tree heights of the LYS and YS plots has a unmoral and bimodal distribution, respectively, as shown in Figure 2. The stand densities varied from 435 to 4097 stems/hm 2 . The descriptive statistics of plots are summarized in Table 1.

Field Data
There are 19 blocks regularly distributed in the study area, and four plots of 25 m × 25 m were designed in each block (Figure 1b). The diameter at breast heights (DBHs), species, tree heights (m), and crown widths (m) at perpendicular directions of north-south and east-west of individual trees with DBH ≥ 5 cm were measured in September 2019. The number of plots dominated with L. principis-rupprechtii (denoted as LYS) and P. tabuliformis (denoted as YS) were 30 and 41, respectively, and used as references to assess the accuracy of the estimated stand heights.
The heights of individual trees of all plots varied from 2 m to 24.1 m with a mean of 11.23 m. The mean of tree heights of the LYS and YS plots were 11.67 m and 11.05 m, respectively. The histogram of tree heights of the LYS and YS plots has a unmoral and bimodal distribution, respectively, as shown in Figure 2. The stand densities varied from 435 to 4097 stems/hm 2 . The descriptive statistics of plots are summarized in Table 1.    [34]. H Med is the median height of each plot.
where h i is the height of i th individual tree; n is the number of sample plots; G i is basal area at 1.3 m height of i th individual tree; and D i is diameter of i th individual tree. The boxplots of the four forest stand heights within different forest cases are shown in Figure 3, which indicated that the interquartile range of ALL plots ( Figure 3a) and LYS plots ( Figure 3b) were similar, while the medians of the four forest heights from LYS plots were significantly greater than that of ALL plots. For YS plots (Figure 3c), the interquartile range of the three stand heights including Lorey's height, mean height, and median height were shorter than those of the LYS and ALL plots. The Lorey's height (HL), mean height (HA), dominated height (HDom), and median height (HMed) were calculated as forest stand height. HL, calculated using DBH and height, the two critical parameters of each tree, could better represent the stand height. HA is the average height of all individual trees with DBH ≥5 cm in each ground plot. HDom is the average height of 20% top trees [34]. HMed is the median height of each plot.
where hi is the height of i th individual tree; n is the number of sample plots; Gi is basal area at 1.3 m height of i th individual tree; and Di is diameter of i th individual tree. The boxplots of the four forest stand heights within different forest cases are shown in Figure 3, which indicated that the interquartile range of ALL plots ( Figure 3a) and LYS plots ( Figure 3b) were similar, while the medians of the four forest heights from LYS plots were significantly greater than that of ALL plots. For YS plots (Figure 3c), the interquartile range of the three stand heights including Lorey's height, mean height, and median height were shorter than those of the LYS and ALL plots.

ULS and USP Datasets
The ULS datasets were acquired in September and October 2019 over the 19 blocks. A Riegl VUX-1 laser scanner mounted on a RC6-2000 UAV aircraft with eight-rotors was used to acquire LiDAR data at an altitude ranging from 933 m to 1322 m above ground level. The flying height was determined by considering the variation of terrain and safety of flight. The differences in height between the lowest and highest points of 19 blocks varied from approximately 162 m to 263 m with a mean of 203 m. The actual flying heights above terrain varied from 58 m to 340 m with a mean of 202 m. The original point clouds were georeferenced using the high-precision position and orientation (POS) data, and

ULS and USP Datasets
The ULS datasets were acquired in September and October 2019 over the 19 blocks. The USP data were acquired in September 2019 over the 19 blocks on a sunny day using a DJI Phantom 4 RTK, a light-weight and low-cost aircraft carrying a consumer-grade camera. The designed flying heights were approximately 190 m above the ground level. The images were obtained with along-and across-track overlap of 80% and 70% in the 19 blocks, which had a ground sampling resolution of 0.1 m. A total of 14,000 images were acquired during 48 flights across the 19 blocks. The detailed specifications of the UAV laser scanning and photogrammetry system are listed in Table 2. The dense three-dimensional point clouds for the 19 blocks were generated with high overlapped images using Agisoft Metashape professional (Version 1.6.2 build 10247). The processing mainly includes images checking, feature identification, matching, and bundle adjustment [35]. The specific data processing instructions of this study can be found in [36]. The point densities of 19 blocks varied from 40 to 75 pts/m 2 at 1.0 m resolution.

Point Cloud Normalization and CHM Generation
To make the comparison more accurate and convincing, the registration of LiDAR point clouds and photogrammetry point clouds was primarily conducted by visual and manual editing based on the artificial objects such as buildings and roads [37]. The point clouds were classified into four classes including noise points, ground points, building points, and vegetation points using the Terrasolid software. Noise points were isolated points with obviously higher or lower height value compared with other points, and could be identified using the isolation algorithm or visual evaluation. The isolation algorithm was identified based on the straight distance between two points. As for one point A, if there is only point A in the spherical space with a radius of 5 m and point A as the center, then point A is identified as one noise point. If there are two or more points in the spherical space, point A is not a noise point. The radius in this study was set to 5 m based on the empirical value. Building points such as house, electric power, and lines could easily be extracted by visual and manual editing based on their regular shape. The triangulated irregular network (TIN) algorithm was applied to extract ground point clouds. Subsequently, the digital elevation model (DEM) with a resolution of 0.2 m was generated from ground points using the (TIN) interpolation method [38]. Then, the classified point clouds were normalized by subtracting the DEM in which the noise and building point clouds were excluded to improve the accuracy of the normalized point clouds. The height value of each normalized point equaled the difference value between the height value of the original point and the corresponding terrain elevation. DSM with a resolution of 0.2 m was generated from the ground and vegetation points using the maximum interpolation method [39]. CHM was directly generated by subtracting DEM from DSM. The cell value of the canopy height model (CHM) equaled the difference value between the height value of the highest point within that cell and the corresponding terrain elevation.
To make the results from ULS and USP more comparable, the NPCs of UAV photogrammetry (P-NPC) were created by subtracting the ULS DEM (L-DEM), and then the CHM of photogrammetry (P-CHM) was created by subtracting L-DEM from P-DSM.
In addition, to analyze the impact of point density on the estimated stand heights, we thinned the original density datasets to 0.8, 1, 5, 10, and 30 pts/m 2 using the ThinData model in Fusion software, respectively [40]. The algorithm first figures out the cover area of each input point cloud data, then divides the data into grids with a certain cell size, and identifies the thinned ratio based on the target point cloud density and the pulse density of each grid cell. Finally, the thinning is performed by randomly removing until achieving the desired pulse density. In this study, the certain grid cell size was set to the default value 5 m × 5 m. To evaluate the uncertainty in the thinning process and reduce the random effects from the thinned point clouds, we executed 30 random repetitions for original point cloud datasets to the target density. An example of thinned point clouds can be found in Figures 4 and 5. It can be found from Figure 4 that the thinning algorithm evenly thinned the point clouds. Correspondingly, the vertical distribution profiles did not have significant differences except for lower point densities (0.8 pts/m 2 , 1 pts/m 2 ), as seen in Figure 5. were normalized by subtracting the DEM in which the noise and building point clouds were excluded to improve the accuracy of the normalized point clouds. The height value of each normalized point equaled the difference value between the height value of the original point and the corresponding terrain elevation. DSM with a resolution of 0.2 m was generated from the ground and vegetation points using the maximum interpolation method [39]. CHM was directly generated by subtracting DEM from DSM. The cell value of the canopy height model (CHM) equaled the difference value between the height value of the highest point within that cell and the corresponding terrain elevation.
To make the results from ULS and USP more comparable, the NPCs of UAV photogrammetry (P-NPC) were created by subtracting the ULS DEM (L-DEM), and then the CHM of photogrammetry (P-CHM) was created by subtracting L-DEM from P-DSM.
In addition, to analyze the impact of point density on the estimated stand heights, we thinned the original density datasets to 0.8, 1, 5, 10, and 30 pts/m 2 using the ThinData model in Fusion software, respectively [40]. The algorithm first figures out the cover area of each input point cloud data, then divides the data into grids with a certain cell size, and identifies the thinned ratio based on the target point cloud density and the pulse density of each grid cell. Finally, the thinning is performed by randomly removing until achieving the desired pulse density. In this study, the certain grid cell size was set to the default value 5 m × 5 m. To evaluate the uncertainty in the thinning process and reduce the random effects from the thinned point clouds, we executed 30 random repetitions for original point cloud datasets to the target density. An example of thinned point clouds can be found in Figures 4 and 5. It can be found from Figure 4 that the thinning algorithm evenly thinned the point clouds. Correspondingly, the vertical distribution profiles did not have significant differences except for lower point densities (0.8 pts/m 2 , 1 pts/m 2 ), as seen in Figure 5.

Feature Metrics Generation
The normalized point clouds and CHM of 71 plots from ULS and USP were all clipped out according to their corresponding boundaries. The point clouds or the CHM cells with height <2 m were excluded, which were considered as echoes from shrubs, grass, or ground surface rather than from the canopy. For comparison of the performance of ULS and USP in describing vertical distribution and height estimation, the commonly used height metrics and density metrics were both derived from NPC or CHM. In this study, a total of 20 metrics (Table 3) were extracted from L-NPC, P-NPC, L-CHM, and P-CHM respectively, including 10 height percentiles (h10, h20, h30, ….., h90, and h95), 6 height statistical metrics (hmax, hcv, hsd, hmean, hmed, and hIQ), and 4 density statistical metrics [41]. The hcv and hsd were the coefficients of variation and the standard deviation of heights representing the variation and heterogeneity of forest canopy height. In summary, four datasets including normalized point clouds from ULS and USP (L-NPC, P-NPC), CHM from ULS, and USP (L-CHM and P-CHM) were used in this study. Percentage of point cloud above 2 m CCmean Percentage of point cloud above mean CCmode Percentage of point cloud above mode

Feature Metrics Generation
The normalized point clouds and CHM of 71 plots from ULS and USP were all clipped out according to their corresponding boundaries. The point clouds or the CHM cells with height <2 m were excluded, which were considered as echoes from shrubs, grass, or ground surface rather than from the canopy. For comparison of the performance of ULS and USP in describing vertical distribution and height estimation, the commonly used height metrics and density metrics were both derived from NPC or CHM. In this study, a total of 20 metrics (Table 3) were extracted from L-NPC, P-NPC, L-CHM, and P-CHM respectively, including 10 height percentiles (h 10 , h 20 , h 30 , . . . , h 90 , and h 95 ), 6 height statistical metrics (h max , h cv , h sd , h mean , h med , and h IQ ), and 4 density statistical metrics [41]. The h cv and h sd were the coefficients of variation and the standard deviation of heights representing the variation and heterogeneity of forest canopy height. In summary, four datasets including normalized point clouds from ULS and USP (L-NPC, P-NPC), CHM from ULS, and USP (L-CHM and P-CHM) were used in this study.

Data Analysis
To have a systematic and comprehensive analysis of ULS and USP, we quantified the spatial distribution of the two data source based on their metrics and compared their performance on the forest common stand height estimation. In detail, two methods were Remote Sens. 2021, 13, 2885 9 of 23 used to compare and analyze the differences between the ULS and USP. The first was by comparing the correlation and difference between each pair of metrics from the four datasets such as h 10 from L-NPC and h 10 from P-NPC using the Pearson's coefficients (r) and mean difference (MD).
where l i is the estimated stand height of i th from ULS data and p i is the estimated stand height of i th from USP data.
The second method was to analyze the capability of the four datasets to estimate the forest stand height. In this paper, the simple linear regression was employed to model the forest stand height using metrics. The best metrics with the highest R 2 from one dataset were extracted as the explanatory variable for modeling simple linear regression. Then, 10-fold cross validation was applied for estimating the accuracy of each linear regression model using R 2 , RMSE, and rRMSE based on the observed values and the corresponding estimated values. In this study, three forest cases were considered including ALL plots, LYS plots alone, and YS plots alone, respectively. In total, 48 models for original point cloud density and 40 models for five thinned point cloud densities were analyzed in this study.
where n is the number of plots; y i is the measured stand height of i th plot;ŷ i is the predicted stand height of i th plot; and y is the mean of observed stand height.
In addition, to explore the effect of point densities on the estimation accuracy of forest stand heights and analyze the estimation accuracy of ULS and USP data at the same point density level, four stand heights were estimated using the thinned data and the same explanatory variable as the original point density of the ULS and USP data. The predication accuracy was also assessed using R 2 , RMSE, and 10-fold cross validation.

Visual Comparison of ULS Point Clouds and USP Point Clouds
Based on the CC of all plots, three plots (Table 4) representing three different CC levels (high CC, median CC, low CC), respectively, were selected to demo and visually compare the differences between ULS and USP point clouds. The height ranges of NPC from USP were mostly consistent with those from ULS. The minimum heights of the NPC were all equal to 0 except for USP NPC within the high CC plot, the maximum heights of P-NPC (8.59 m~17.02 m) were all less than the corresponding L-NPC (8.65 m~17.36 m). From the 3D distribution of point clouds within different CC plots from ULS (Figure 6a1-a3), it can be seen that there were more and more point clouds hitting the ground with the decreasing CC for ULS. In contrast, there were no corresponding changes and there were obviously less point clouds from the understory and ground for USP point clouds (Figure 6b1-b3), even for the plot with low CC = 0.49. In addition, the extracted ground point clouds from USP mostly gathered together where there was a big gap between crowns. It can be seen from the profiles of the original point clouds (Figure 6c1-c3) and normalized point cloud (Figure 6d1-d3) that the point clouds generated from USP were mostly located around the canopy surface, and had high coincidence with the point clouds extracted from ULS, except for some point clouds from tree tops. The different distribution of USP dense point clouds and ULS under the canopy was subjected to the different principles of the two data sources. The USP did not have the ability to penetrate the canopy, which is the difference from ULS.  From the 3D distribution of point clouds within different CC plots from ULS ( Figure  6a1-a3), it can be seen that there were more and more point clouds hitting the ground with the decreasing CC for ULS. In contrast, there were no corresponding changes and there were obviously less point clouds from the understory and ground for USP point clouds (Figure 6b1-b3), even for the plot with low CC = 0.49. In addition, the extracted ground point clouds from USP mostly gathered together where there was a big gap between crowns. It can be seen from the profiles of the original point clouds (Figure 6c1-c3) and normalized point cloud (Figure 6d1-d3) that the point clouds generated from USP were mostly located around the canopy surface, and had high coincidence with the point clouds extracted from ULS, except for some point clouds from tree tops. The different distribution of USP dense point clouds and ULS under the canopy was subjected to the different principles of the two data sources. The USP did not have the ability to penetrate the canopy, which is the difference from ULS.

Forest Stand Heights Modeling Using ULS and USP Metrics
The correlations between forest stand heights and height metrics from L-NPC, L-CHM, P-NPC, and P-CHM are shown in Figure 7a1,a2,b1,b2 and Figure 8a1,a2,b1,b2, respectively, from which we could find that the correlations between NPC or CHM metrics from USP ( Figure 8) and forest stand heights were consistent with the corresponding metrics from ULS ( Figure 7). Overall, each forest height had strong correlation with height percentiles (h 10~h95 ), h med and h mean , and less correlation with height statistical metrics h cv , h sd , and h IQ with R 2 ranging from 0.006 to 0.460 and RMSE varying from 10.015 to 13.06 m.
The R 2 and RMSE of height percentiles (h 10~h95 ), h med , h mean , and forest stand heights ranged from 0.734 to 0.944, 0.861 m to 7.681 m for L-NPC, respectively, and from 0.735 to 0.952, 0.822 m to 7.518 m for L-CHM, respectively. Like with the USP data, the R 2 and RMSE of the height metrics ranged from 0.800 to 0.942, 0.890 m to 5.517 m for the NPC metrics, respectively, and from 0.795 to 0.942, 0.868 m to 5

Forest Stand Heights Modeling Using ULS and USP Metrics
The correlations between forest stand heights and height metrics from L-NPC, L-CHM, P-NPC, and P-CHM are shown in Figures 7a1,a2,b1,b2 and 8a1,a2,b1,b2, respectively, from which we could find that the correlations between NPC or CHM metrics from USP ( Figure 8) and forest stand heights were consistent with the corresponding metrics from ULS (Figure 7). Overall, each forest height had strong correlation with height percentiles (h10~h95), hmed and hmean, and less correlation with height statistical metrics hcv, hsd, and hIQ with R 2 ranging from 0.006 to 0.460 and RMSE varying from 10.015 to 13.06 m.
The R 2 and RMSE of height percentiles (h10~h95), hmed, hmean, and forest stand heights ranged from 0.734 to 0.944, 0.861 m to 7.681 m for L-NPC, respectively, and from 0.735 to 0.952, 0.822 m to 7.518 m for L-CHM, respectively. Like with the USP data, the R 2 and RMSE of the height metrics ranged from 0.800 to 0.942, 0.890 m to 5.517 m for the NPC metrics, respectively, and from 0.795 to 0.942, 0.868 m to 5.463 m for the CHM metrics, respectively.
HA and HMed had a higher correlation with the lower height percentiles (h30, h40), hmed, and hmean, while HL and HDom had a higher correlation with the upper height percentiles such as h80 and h90.   Based on the strong correlations between height metrics and each forest stand height, simple linear regression for HA, HL, HDom, and HMed using height metrics generated from the ULS and USP datasets was applied to fit each forest stand height. The best explanatory variables for different models, and corresponding modeling accuracy of simple linear regression models using NPC and CHM metrics alone extracted from ULS and USP, respectively, are summarized in Figure 9 and Table 5.
Compared with the best explanatory variables obtained from L-NPC, the explanatory variables obtained from P-NPC mostly took the lower height percentiles. For instance, h80, h90, and h80 of L-NPC were the best explanatory variables for HL in the corresponding forest types (ALL plots, LYS plots, and YS plots, respectively), while h70, h80, and h60 obtained from P-NPC were the corresponding best explanatory variables for HL. This tendency was also obvious for CHM metrics.
From Figure 9, it can also be seen that the best explanatory variables for the four forest stand heights modeled from P-NPC and P-CHM metrics were basically the same, except for the HDom of LYS plots, while the best explanatory variables from L-NPC metrics were higher or the same as that from the L-CHM metrics.
The best explanatory variables varied with different stand heights from lower height percentiles to upper height percentiles. In an item of HL and HDom, their best explanatory variables were mostly the upper height percentiles such as h80, h90, and h95, while the lower and median height percentiles were the best independent variables for HA and HMed such as h30, h40, h50, and h60. For different forest types, the best explanatory variables for the same stand height varied within a certain range. The best independent variables for forest stand heights from the LYS plots were higher than those of the YS plots. For example, for HA, HL, HDom, and HMed modeled from the L-NPC metrics, their corresponding best variables h60, h90, h95, and h70 from the LYS plots were higher than h40, h80, h90, and h30 from the YS plots, respectively. Based on the strong correlations between height metrics and each forest stand height, simple linear regression for H A , H L , H Dom , and H Med using height metrics generated from the ULS and USP datasets was applied to fit each forest stand height. The best explanatory variables for different models, and corresponding modeling accuracy of simple linear regression models using NPC and CHM metrics alone extracted from ULS and USP, respectively, are summarized in Figure 9 and Table 5.
Compared with the best explanatory variables obtained from L-NPC, the explanatory variables obtained from P-NPC mostly took the lower height percentiles. For instance, h 80 , h 90 , and h 80 of L-NPC were the best explanatory variables for H L in the corresponding forest types (ALL plots, LYS plots, and YS plots, respectively), while h 70 , h 80 , and h 60 obtained from P-NPC were the corresponding best explanatory variables for H L . This tendency was also obvious for CHM metrics.
From Figure 9, it can also be seen that the best explanatory variables for the four forest stand heights modeled from P-NPC and P-CHM metrics were basically the same, except for the H Dom of LYS plots, while the best explanatory variables from L-NPC metrics were higher or the same as that from the L-CHM metrics.
The best explanatory variables varied with different stand heights from lower height percentiles to upper height percentiles. In an item of H L and H Dom , their best explanatory variables were mostly the upper height percentiles such as h 80 , h 90 , and h 95 , while the lower and median height percentiles were the best independent variables for H A and H Med such as h 30 , h 40 , h 50 , and h 60 . For different forest types, the best explanatory variables for the same stand height varied within a certain range. The best independent variables for forest stand heights from the LYS plots were higher than those of the YS plots.     Figure 10 shows the scatterplots of the estimated stand heights using 10-fold cross validation with L-NPC and L-CHM metrics alone. Overall, the L-NPC and L-CHM metrics both performed best for estimating the forest stand heights with R 2 = 0.88~0.96, RMSE = 0.77~1.23 m, and rRMSE = 5.69~9.72%. The performances of the L-CHM metrics models were slightly higher than the results of L-NPC (∆R 2 = 0~0.01, ∆RMSE = −0.08~0 m, ∆rRMSE = −0.61~0.05%).

Forest Stand Height Estimation Accuracy
Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 23 Figure 10 shows the scatterplots of the estimated stand heights using 10-fold cross validation with L-NPC and L-CHM metrics alone. Overall, the L-NPC and L-CHM metrics both performed best for estimating the forest stand heights with R 2 = 0.88~0.96, RMSE = 0.77~1.23 m, and rRMSE = 5.69~9.72%. The performances of the L-CHM metrics models were slightly higher than the results of L-NPC (ΔR 2 = 0~0.01, ΔRMSE = −0.08~0 m, ΔrRMSE = −0.61~0.05%).  The scatterplots of the estimated stand heights using the best explanatory variables from the P-NPC and P-CHM metrics, respectively, are exhibited in Figure 11, from which it can be found that the P-NPC and P-CHM metrics both performed best for estimating forest stand heights with R 2 = 0.

Estimation Results of Forest Stand Heights with Different Point Density
The influence of point density on the estimation accuracy of the four forest stand heights using L-NPC and P-NPC metrics alone are shown in Figure 12. It can be found that the estimation accuracy of each forest stand height was little affected by point density, slightly increased with the point density changing from 0.8 pts to 30 pts/m 2 , and became more stable, which could also be found from the corresponding RMSE results. As such, mean R 2 varied from 0.933 (±0.002) to 0.935 (±0) with increasing point density for H A using the L-NPC data, and from 0.923 (±0.002) to 0.931 (±0) with increasing point density for H A using the P-NPC data. Mean RMSE varied from 1.040 (±0.015) to 1.036 (±0.003) with increasing point density for H A using the L-NPC data, and from 1.056 (±0.014) to 1.037 (±0.001) with increasing point density for H A using the P-NPC data.
In addition, the estimation accuracy of the four forest stand heights modeled by L-NPC metrics was slightly better than those of USP, except for H Med (Figure 12a4,b4). The H Med estimated using USP was higher than ULS. H Dom exhibited the least variation with increasing point density (Figure 12a3,b3).

Estimation Results of Forest Stand Heights with Different Point Density
The influence of point density on the estimation accuracy of the four forest stand heights using L-NPC and P-NPC metrics alone are shown in Figure 12. It can be found that the estimation accuracy of each forest stand height was little affected by point density, slightly increased with the point density changing from 0.8 pts to 30 pts/m 2 , and became more stable, which could also be found from the corresponding RMSE results. As such, mean R 2 varied from 0.933 (±0.002) to 0.935 (±0) with increasing point density for HA using the L-NPC data, and from 0.923 (±0.002) to 0.931 (±0) with increasing point density for HA Figure 11. The scatterplot of estimated forest stand height and field measured height using USP metrics. (a1-a4,b1-b4) for P-NPC. (c1-c4,d1-d4) for P-CHM.
using the P-NPC data. Mean RMSE varied from 1.040 (±0.015) to 1.036 (±0.003) with increasing point density for HA using the L-NPC data, and from 1.056 (±0.014) to 1.037 (±0.001) with increasing point density for HA using the P-NPC data.
In addition, the estimation accuracy of the four forest stand heights modeled by L-NPC metrics was slightly better than those of USP, except for HMed (Figure 12a4 and b4). The HMed estimated using USP was higher than ULS. HDom exhibited the least variation with increasing point density (Figure 12a3 and b3).

Discussion
In this study, a low-cost and consumer-grade UAV was used to acquire photogrammetry images. Meanwhile, LiDAR point clouds were acquired using a Riegl VUX-1 laser

Discussion
In this study, a low-cost and consumer-grade UAV was used to acquire photogrammetry images. Meanwhile, LiDAR point clouds were acquired using a Riegl VUX-1 laser scanner mounted on a RC6-2 UAV aircraft with eight-rotors to better assess the threedimensional mapping of photogrammetry images, mainly including the capacity of ULS and USP data for forest stand height estimation, and the effect of point cloud density on the estimation accuracy of the four stand heights. The results demonstrated that the USP performed basically as well as ULS for the continuous coniferous plantation, except for the ability of generating terrain information, which were basically in line with the outcomes from previous studies [28].

The Capability of USP and ULS to Estimate Forest Stand Heights
For the four forest stand heights, our results showed that USP performed basically as well as ULS in estimating the forest stand heights, indicating that USP with LiDAR DEM could also achieve the same level of accuracy for stand height estimation at the plot level as well as ULS. Meanwhile, the estimated results from CHM were also similar with those of the point clouds.
Among the four forest stand heights, Lorey's height was the most accurately estimated, followed by dominated height, mean height, and median height, indicating that Lorey's height was the best parameter representing the forest stand height, which can also be found in [4,42]. The four stand heights were calculated using different methods. Lorey's height is the weighted height of DBH and tree height, and has better representativeness and stability for a stand with different ages. The dominate height was calculated using the 20% tallest trees, which mostly came from the upper trees rather than depressed trees and could be easily detected from CHM or point clouds as a result of having a strong correlation with the upper height percentiles. Furthermore, with the increasing height, the characterization of the canopy became more and more accurate, especially for CHM and point clouds from photogrammetry. However, mean height and median height took the average and median of all heights within each plot, respectively, which were directly related to the height distribution of each plot including the pressed tree in the understory, and also easily affected by the measurement error, causing a lower estimation accuracy than Lorey's height and dominated height.
Additionally, compared with the results using other parametric [43] or non-parametric methods [27] conducted in other studies, the simple linear regression utilized in this study was confirmed with high accuracy and having a better performance, which may be attributed to the characteristics of coniferous plantations.

The Best Explanatory Variables for Forest Stand Heights
The selected best explanatory variables had direct theoretical significance with the four stand heights and were mostly consistent with other studies. The difference between the best variables for ULS and USP were mainly caused by the vertical height distribution characteristics of its height metrics (Tables A1 and A2). The phenomenon that the best explanatory variable from P-NPC and P-CHM were basically consistent was attributed to the lower variation in height for the USP data, causing the same height percentile to take a similar height value. On the other hand, the variety of the best explanatory variables from L-NPC and L-CHM revealed that ULS was more capable of penetrating the forest canopy.
The best explanatory variables for all stand heights were all the height percentiles, which demonstrated the strong correlation between stand heights of the coniferous plantation with the height percentiles. The best explanatory variables for Lorey's height were mostly upper height percentiles such as h 90 and h 95 , which was consistent with the results in [4,18]. Mean height was strongly correlated with median height percentiles (h 40 , h 50 , and h 60 ), which were affected by the whole distribution of all heights of each plot. Jensen and Mathews [28] estimated the H A and H Med using h mean and h med metrics by simple linear regression with R 2 = 0.90 and R 2 = 0.89, respectively. The high height percentiles were the best variables for dominated height because the dominated height was calculated using the 20% tallest heights of each plot, representing the characteristics of upper trees. Puliti et al. [12] estimated the H L and H Dom using multiple linear regression by taking the height percentiles and density metrics as explanatory variables with R 2 = 0.71, RMSE = 1.4 m and R 2 = 0.97, RMSE = 0.7 m, respectively.

The Influence of Point Density on the Estimation Accuracy
The varied flying height above the terrain will surely affect the spatial distribution of point density. To examine the sensitivity of forest height estimation to point density, five point densities varying from 0.8 pts to 1 pts, 5 pts, 10 pts, and 30 pts/m 2 for ULS point cloud and USP point cloud data were analyzed in this study. The predication accuracy of forest stand heights were slightly affected by thinned point clouds, which was mostly attributed to the height percentile metrics from different point densities showing little variability [4]. The results were also in keeping with other studies [44]. Consequently, it may be possible to take a higher flight height or lower overlap to acquire lower point cloud density, which could collect the data over the same area using a shorter time, thus potentially reducing the cost.

The Performances of Point Clouds and CHM
The performances were similar between the point clouds and CHM for the estimation accuracy of forest stand heights. P-CHM metrics models were mostly consistent with the results of P-NPC, rather than the correlation of the L-CHM and L-NPC metrics. As is widely known to all, USP is passive remote sensing, where the dense point clouds are extracted depending on the features existing in the images and SFM algorithms. The dense point clouds of USP were mostly gathered at the mid-to-upper canopy with lower penetration through the canopy to the ground, resulting in less height difference for the point clouds and CHM while the 3D point clouds from ULS were distributed throughout the forest canopy and ground surface [23], thus having the ability to depict more detailed characterization of the forest canopy. The distribution of point clouds resulted in differences between the corresponding height metrics.

The Limitation of This Study
Stereo photogrammetry has the ability to accurately estimate the large solitary crowns by individual tree detection approaches such as canopy local maxima filtering [45], watershed segmentation [45,46], and region growing algorithm [47,48]. However, in our study site, the tree crowns are severely occluded by each other so it is difficult to accurately isolate individual trees. Additionally, the growth of the plantation in our study site is limited by crucial environmental conditions such as precipitation and energy of photosynthesis. Tall trees will usually fall down and saplings will be planted in the open areas in a rotation period. It is almost impossible to find trees with a height >35 m. While the study might be more representative to consider mature forest with >35 m heights, we would like to explore tall trees in different study sites in future study.
The point cloud generated from USP could be normalized using the DEM based on LiDAR point cloud or other datasets. LiDAR has the ability to accurately acquire flat or complex terrain. Additionally, a good DEM can also be acquired by other methods such as ground surveying, which will be rather adequate for relatively flat terrain.
The difference between different forest types may be attributed to the distribution of heights and the morphological characteristics of Larixprincipis-rupprechtii and Pinus tabuliformis [49]. It can be seen from Figure 2 that the height distributions between LYS plots and YS plots were different, which might be caused by the ratio of different composition of stand ages for the LYS plots and YS plots, as the LYS plots had many young trees. Additionally, the growth rates and canopy shapes for the LYS and YS species were different and the canopy shape of LYS was an ellipsoid shape, while it was an inverted triangle distribution for YS within the growing age, which can also be found from the best explanatories of forest stand heights in Figure 9 for LYS and YS. The best explanatories for LYS were higher than those for YS. Racine et al. [50] also found that the species, crown closure, and age of balsam fir (Abies balsamea (L.) Mill.), paper birch (Betula papyrifera Marsh), black spruce (Picea mariana (Mill.) BSP), white spruce (Picea glauca Moench), and aspen (Populus tremuloides Michx.) within even-aged stands had an influence on their vertical distribution based on the airborne LiDAR data.

Conclusions
UAV stereo photogrammetry UAV laser scanning can be used to extract plantation height metrics in north China. For the four forest stand heights, the predication accuracy of Lorey's height was the highest, followed by dominated height, mean height, and median height, which meets the requirements of the plantation forest resource survey in this study. These heights can be estimated with high accuracy by simple linear regression based on the plot-level metrics. The USP cannot obtain accurate terrain information in median and high forest CC in the leaf-on season, but can be used to estimate forest heights with high accuracy combined with terrain from ULS, even with low point density. There were strong correlations between the metrics from the two UAV techniques. The forest stand heights estimated using the metrics extracted from NPC or CHM of USP were in agreement with those of ULS. The USP data including both NPC metrics and CHM metrics performed well in estimating forest stand heights, similar to the prediction accuracy of the ULS data. Additionally, the point cloud density that reduced from 30 pts to 0.8 pts had little influence on the prediction accuracy. Furthermore, point cloud and stereo images should be acquired in more circumstances such as broadleaf forest and leaf-off season to investigate the forest structure parameters in future studies.