Assessment of Individual Tree Detection and Canopy Cover Estimation using Unmanned Aerial Vehicle based Light Detection and Ranging (UAV-LiDAR) Data in Planted Forests

: Canopy cover is a key forest structural parameter that is commonly used in forest inventory, sustainable forest management and maintaining ecosystem services. Recently, much attention has been paid to the use of unmanned aerial vehicle (UAV)-based light detection and ranging (LiDAR) due to the ﬂexibility, convenience, and high point density advantages of this method. In this study, we used UAV-based LiDAR data with individual tree segmentation-based method (ITSM), canopy height model-based method (CHMM), and a statistical model method (SMM) with LiDAR metrics to estimate the canopy cover of a pure ginkgo ( Ginkgo biloba L.) planted forest in China. First, each individual tree within the plot was segmented using watershed, polynomial ﬁtting, individual tree crown segmentation (ITCS) and point cloud segmentation (PCS) algorithms, and the canopy cover was calculated using the segmented individual tree crown (ITSM). Second, the CHM-based method, which was based on the CHM height threshold, was used to estimate the canopy cover in each plot. Third, the canopy cover was estimated using the multiple linear regression (MLR) model and assessed by leave-one-out cross validation. Finally, the performance of three canopy cover estimation methods was evaluated and compared by the canopy cover from the ﬁeld data. The results demonstrated that, the PCS algorithm had the highest accuracy ( F = 0.83), followed by the ITCS ( F = 0.82) and watershed ( F = 0.79) algorithms; the polynomial ﬁtting algorithm had the lowest accuracy ( F = 0.77). In the sensitivity analysis, the three CHM-based algorithms (i.e., watershed, polynomial ﬁtting and ITCS) had the highest accuracy when the CHM resolution was 0.5 m, and the PCS algorithm had the highest accuracy when the distance threshold was 2 m. In addition, the ITSM had the highest accuracy in estimation of canopy cover ( R 2 = 0.92, rRMSE = 3.5%), followed by the CHMM ( R 2 = 0.94, rRMSE = 5.4%), and the SMM had a relative low accuracy ( R 2 = 0.80, rRMSE = 5.9%). The UAV-based LiDAR data can be e ﬀ ectively used in individual tree crown segmentation and canopy cover estimation at plot-level, and CC estimation methods can provide references for forest inventory, sustainable management and ecosystem assessment.


Introduction
Planted forests are a key component of global forests and play vital roles in the restoration and reconstruction of the ecological environment, mitigating global climate change and promoting sociological and economic development [1][2][3]. China has the largest planted forest area, which is nearly 6.93 × 10 7 ha and accounts for 26.2% of the world's total planted forest area [4,5]. Ginkgo (Ginkgo biloba L.) is a characteristic planted forest species in China [6]. According to historical statistics, there are approximately 280,000 mature ginkgo trees in China, and the annual fruit yield has reached 16,000 tons [7]. The areas and stocks of ginkgo planted forests in China, which not only help to supplement national timber and forest products but also play a crucial role in the sustainable management and development of forestry and the maintenance of the carbon balance, have increased rapidly in recent years [7][8][9].
Forest canopy cover is defined as the ratio of the vertical projection area of tree crowns to the forest area [10,11]. Forest canopy cover is one of the main indicators used to determine whether land cover is a forest [12] and is usually utilized to assess surface microclimate and light conditions that affect ecosystem energy transfer and carbon sink contents [13][14][15]. The accurate acquisition of canopy cover information for ginkgo planted forests is essential for the characterization of forest ecosystems and estimation of the harvest of ginkgo products [10,16]. The traditional approach for canopy cover measurements was field measurements with artificial, most of which were collected by looking upward from sample points, which was time consuming and labor intensive. Notably, this approach can only obtain canopy cover data at a small scale [13,17,18].
Remote sensing can acquire forest canopy cover data rapidly and accurately at a large scale [19][20][21]. Light detection and ranging (LiDAR) is an active remote sensing technology that transmits laser pulses to the surface of objects and records return signals. The three-dimensional information of surface entities, such as terrain and vegetation, can be obtained by LiDAR with high precision in real time. In theory, a vertical accuracy of 15-30 cm can be achieved [22][23][24][25]. Compared with traditional optical remote sensing technologies, LiDAR has an increased penetration and is not easily affected by weather conditions; therefore, this technology has unique advantages in obtaining forest structure information. LiDAR data may even be used to assess the quality of field measurements [26,27]. Canopy cover field measurements are obtained by determining the ratio of the total crown area to the sample plot area [13]. Uncertainty factors such as artificial measurement errors and inconsistent standards usually affect ground measurements of crown amplitudes, resulting in low-accuracy canopy cover (CC) estimations [18]. Previous studies have shown that high-density LiDAR point clouds can separate individual trees and extract crown width data accurately [28][29][30][31].
In recent years, the estimation of CC from LiDAR data has been mainly classified into two categories: the number of canopy pixels in a threshold of rasterized canopy height model (CHM) and the ratio between the number of canopy pixels and the total number of pixels, which is known as the CHM-based method [20,[32][33][34]. For example, Korhonen et al. (2011) extracted the CHM from airborne discrete LiDAR data and smoothed the CHM surface, thus estimating the coniferous forest CC in southern Finland; the root mean square error (RMSE) of the predictive model was 6.6-16.1%. Ma et al. (2017) compared estimations of the CC of broad-leaved mixed forest using airborne LiDAR, aerial images and satellite images in the western United States. The results showed that the CC was estimated with the highest accuracy by airborne LiDAR data, and the estimation accuracy of the CHM-based method was 0.58 (R 2 , RMSE = 19%). The other methods constructed mathematical models by combining LiDAR metrics with field-measured CC, which is known as the statistical model method [33,34]. Smith et al. (2009) compared estimations of the CC of mixed-conifer forest forests using multispectral data and LiDAR data in the northeastern United States and found that the canopy density estimation with LiDAR metrics was the most accurate (R 2 = 0.78, RMSE = 16.1%) [27]. Melin et al. (2017) used several different remote sensing approaches to estimate the CC of coniferous forests in eastern Finland and found that an inversion model constructed by extracting characteristic variables and the measured canopy density on the ground had the highest accuracy when using LiDAR (R 2 = 0.65, rRMSE = 10.3%) [35].
Unmanned aerial vehicle (UAV)-based remote sensing technology has developed rapidly in recent years. Previous researchers found that UAV-based LiDAR technology can accurately and efficiently estimate forest structure parameters [36,37]. UAV-based LiDAR has a lower cost, more convenient operation and more flexible flight route design than airborne LiDAR as well as unique advantages in the high-precision acquisition of high-density point clouds. Nagai et al. (2009) equipped an unmanned helicopter with a low-frequency laser scanner (SICK LMS-291) and a single-reflection camera (Canon EOS 10D) to extract canopy digital surface model (DSM) and texture features with an average mapping error of approximately 3-10 cm [38]. Jaakkola et al. (2010) integrated a laser scanner (Sick LMS 151), a Charge Coupled Device (CCD) camera (AVT Pike F-421C), and a thermal infrared camera (Specim V10H) in a lightweight multirotor UAV-LiDAR system (FGI Sensei) to estimate the height of coniferous forests in Finland (standard deviation was 13 cm, R 2 = 0.92) [39]. Wallace et al. (2014Wallace et al. ( , height 2016 developed a low-cost multirotor UAV-LiDAR system (Terra Luma) equipped with a low-frequency scanner (Ibeo LUX) and a high-definition video recorder for individual tree detection (the detection accuracy was 98%) and tree height estimation (R 2 = 0.69 -0.84, RMSE = 0.92 -1.3 m) in southeastern Australia [40,41]. Obtaining high-precision measured CC is the key to verifying the accuracy of CC estimates derived from LiDAR data.
However, most of the previous studies used only a CHM-based method or a statistical model method to estimate CC (few studies using individual tree segmentation method to estimation CC); most of them were applied in boreal and temperate natural forests of Europe and America (few studies on planted forests in subtropics), and most of the data were high-cost airborne LiDAR data obtained by manned aircraft (few studies using UAV-based LiDAR data). In this study, high-point-density UAV-LiDAR data were used to detect individual trees and estimate canopy cover in planted forests. The specific aims of the study were as follows: (i) To compare the accuracy of four LiDAR individual tree algorithms in detecting individual ginkgo trees within plots with three different stem densities. (ii) To compare the accuracy of three CC-estimated methods, these methods were employed to predict the canopy cover and the accuracies of prediction were validated using field data. (iii) Sensitivity analyses of individual tree detection and canopy cover estimation (ITSM and CHMM) were carried out to optimize the parameters of the algorithms, and the accuracy improvement achieved by the statistic model method was evaluated by adding LiDAR metrics.

Materials and Methods
A general overview of the technical route for estimating CC is shown in Figure 1. First, a typical and representative ginkgo planted forest on the northern plain of Jiangsu Province was taken as the study area. The high-density point cloud data obtained by the multirotor UAV-LiDAR system were combined with individual tree coordinate information from ground measurements. In this study, 23 typical plots (three different densities, 892 individual ginkgo trees) of ginkgo planted forests were selected on the plains of China (Tiefu Town, Pizhou City, Jiangsu Province). High-density point cloud data (average point density 159 pts·m −2 ) were obtained by multirotor UAV-LiDAR. Then, using a watershed algorithm, polynomial fitting, individual tree crown segmentation (ITCS) and point cloud segmentation (PCS) to segment individual trees, the accuracy of the four methods was assessed and analyzed (individual tree segmentation). Next, the CC at the sample plot level was estimated by extracting the canopy amplitude of the sample plot with the highest segmentation precision (method 2), the different height-based CHM (method 1) and the LiDAR metrics statistical model (method 3). Finally, the accuracy of three CC estimation methods were evaluated. In order to explore the different results of individual tree detection in different CHM resolutions (for watershed, polynomial fitting and ITCS algorithm) and different distance thresholds (for PCS algorithms), we generated the CHM at 0.3 m, 0.5 m and 0.7 m resolution, and the distance threshold was set as 1 m, 2 m and 3 m.

Study Area
The study area is located in Tiefu Town, Pizhou City, northern Jiangsu Province (34°33'49''-34°34'23''N, 118°05'1''-118°06'06''E). The region has a semi-humid temperate monsoon climate with an annual rainfall of 903 mm. The maximum rainfall is concentrated in the plum rain season of July and August. The annual average temperature is approximately 13.9 °C, and the frost-free period is 211 days. The main soil type is acidic black clayey soil. During 1993 to 1996, ginkgo trees were planted in Tiefu Town with ultra-high-density (approximately 67,500 n/ha). With the growth of ginkgo trees, some forest areas were extensively harvested for improving forest growth conditions (e.g., Light, nutrition area and growing space, etc.). Then, the local farmers managed the ginkgo trees under different cultivation (e.g., seedlings planting, weed clearing and agroforestry management, etc.) and fertilizer treatments from 1996 to the present. Thus, currently, the ginkgo trees represent various ages and tree height in the study area. The low-stem-density plots are mainly mature forests, and the highstem-density plots are mainly young forests.

Field Data
According to historical forest resource survey data, we selected five 1×1 km square UAV sampling areas in the core distribution area of the ginkgo planted forest with the typical sampling method and collected data by using UAV-LiDAR. The ground survey of this study was conducted in late October 2016, and 45 circular plots with a radius of 15 m were established according to the typical sampling method in five square areas. In the process of sample plot investigation, we measured the DBH (Diameter at Breast Height (DBH measurement)), tree height and branch height (Vertex V®

Study Area
The study area is located in Tiefu Town, Pizhou City, northern Jiangsu Province (34 • 33 49"-34 • 34 23" N, 118 • 05 1"-118 • 06 06" E). The region has a semi-humid temperate monsoon climate with an annual rainfall of 903 mm. The maximum rainfall is concentrated in the plum rain season of July and August. The annual average temperature is approximately 13.9 • C, and the frost-free period is 211 days. The main soil type is acidic black clayey soil. During 1993 to 1996, ginkgo trees were planted in Tiefu Town with ultra-high-density (approximately 67,500 n/ha). With the growth of ginkgo trees, some forest areas were extensively harvested for improving forest growth conditions (e.g., Light, nutrition area and growing space, etc.). Then, the local farmers managed the ginkgo trees under different cultivation (e.g., seedlings planting, weed clearing and agroforestry management, etc.) and fertilizer treatments from 1996 to the present. Thus, currently, the ginkgo trees represent various ages and tree height in the study area. The low-stem-density plots are mainly mature forests, and the high-stem-density plots are mainly young forests.

Field Data
According to historical forest resource survey data, we selected five 1×1 km square UAV sampling areas in the core distribution area of the ginkgo planted forest with the typical sampling method and collected data by using UAV-LiDAR. The ground survey of this study was conducted in late October 2016, and 45 circular plots with a radius of 15 m were established according to the typical sampling method in five square areas. In the process of sample plot investigation, we measured the DBH (Diameter at Breast Height (DBH measurement)), tree height and branch height (Vertex V ® hypsometer (Langsele, Sweden)) individually. DBH (cm), Lorey's tree height (m), and stem Remote Sens. 2019, 11, 908 5 of 21 density (n/ha) were included in the plot-level stand characteristic information. The average Lorey's tree height is the average height obtained by weighting the area of the breast height and section of each tree. Moreover, the position of each individual tree in 23 sample plots was acquired using an ultrasound-based Haglöf PosTex positioning instrument (Långsele, Sweden). In this study, only the 23 plots with individual tree position data were used. The canopy cover of these 23 sample plots was measured by the line-intercept method using ultra-high resolution ortho-imagery (0.05 m) in this study. In each sample plot, three 15 m slope-correct lines transects radiating out from the center of each subplot were set (azimuth 0 • , 120 • , and 240 • ) [42]. It includes the entire length within the outline of a crown as cover, each individual tree in a maximum of three vertical canopy layers was added, and the canopy cover was calculated as the ratio of total cover length to the overall three transects lines length (45 m). According to the plot stem density calculated by ground survey data, we divided the plots into three groups (first group: 0-500 n·ha −1 ; second group: 501-600 n·ha −1 ; and third group: 601-810 n·ha −1 ). The study area and plot distribution are shown in Figure 2. The results of the ginkgo planted forest plot groups and information summary are shown in Table 1. hypsometer (Langsele, Sweden)) individually. DBH (cm), Lorey's tree height (m), and stem density (n/ha) were included in the plot-level stand characteristic information. The average Lorey's tree height is the average height obtained by weighting the area of the breast height and section of each tree. Moreover, the position of each individual tree in 23 sample plots was acquired using an ultrasound-based Haglöf PosTex positioning instrument (Långsele, Sweden). In this study, only the 23 plots with individual tree position data were used. The canopy cover of these 23 sample plots was measured by the line-intercept method using ultra-high resolution ortho-imagery (0.05 m) in this study. In each sample plot, three 15m slope-correct lines transects radiating out from the center of each subplot were set (azimuth 0°, 120°, and 240°) [42]. It includes the entire length within the outline of a crown as cover, each individual tree in a maximum of three vertical canopy layers was added, and the canopy cover was calculated as the ratio of total cover length to the overall three transects lines length (45m). According to the plot stem density calculated by ground survey data, we divided the plots into three groups (first group: 0-500 n·ha -1 ; second group: 501-600 n·ha -1 ; and third group: 601-810 n·ha -1 ).The study area and plot distribution are shown in Figure 2. The results of the ginkgo planted forest plot groups and information summary are shown in Table 1.

LiDAR Data Acquisition and Preprocessing
The LiDAR data in our study were obtained from a Velodyne LiDAR carried by an eight-rotor UAV during October 13-17, 2016. The data acquisition parameters of the UAV-based LiDAR are as follows: flight altitude 60 m, flight speed 4.8 m·s -1 , flight interval 50 m, laser wavelength 903 nm, laser scanner scanning angle ±30°, divergence 3 mrad, spot diameter 18 cm, pulse transmitting frequency 21.7 kHz, scanning frequency 16 lines·s -1 . Therefore, we obtained a heading distance of 0.3 m, a distance between side points of 0.05 m and an average point density of approximately 159 pts·m -2 .

LiDAR Data Acquisition and Preprocessing
The LiDAR data in our study were obtained from a Velodyne LiDAR carried by an eight-rotor UAV during noise points of the original data via the height threshold method. Then, the point clouds were classified and the ground points were extracted. Next, the ground points were interpolated by the inverse distance weighted (IDW) method to generate the digital elevation model (DEM). Finally, we used the DEM to normalize the processed point cloud height (z) to generate the CHM.

LiDAR Metrics Calculation
Using normalized discrete point cloud data, the LiDAR metrics at the sample plot scale were calculated. This study selected three suites of metrics, including distributional metrics (n = 14), Weibull-fitting metrics (n = 2) and canopy volume metrics (n = 4). The UAV-LiDAR metrics are summarized in Table 2. Table 2. Description of metrics derived from UAV-LiDAR data.

Metrics Description
Distributional metrics The percentiles of the canopy height distributions by first echo (25th, 50th, 75th and 95th). h mean The mean height of all points after normalized.
h cv The coefficient of variation of height of all points after normalized (the ratio of the standard deviation to the mean). h skewness /h kurtosis The skewness and kurtosis of the heights of all points by first echo The proportion of points above the quantiles(10th, 30th, 50th,70th, and 90th) to total number of points CC 1m The first return points above 1m accounts for the percentage of all return points Weibull-fitting metrics α/β The α and β parameter of the Weibull distribution fitted to foliage density profile. Canopy volume metrics Open/Closed The empty voxels located above and below the canopy respectively.

Euphotic/Oligophotic
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 The Weibull-fitting metrics α and β are calculated by following steps. The foliage profile (FP) is the spatial vertical distribution of tree branches and leaves. FP describes the variation in leaf area index (LAI) with height in the canopy [43]. The FP value is the sum of the branch and leaf area per unit horizontal area within the specified height per unit volume of the surface canopy [44]. The FP is closely related to the vertical distribution of the leaf area index (LAI) [45].
In the above formula, L(z) is the cumulative effective LAI at height Z, and Z 1 and Z 2 are the height of the canopy. In this study, the Weibull function is used to fit FP. The Weibull parameter in this study is not a parameter (model-dependent variable) representing the diameter order of the sample plot but, rather, is a characteristic variable (model-independent variable) calculated by Weibull fitting of the LiDAR vertical branch and leaf section. The formula is as follows: In formula 2, H max is the greatest tree height, α is the scale parameter, and β is the shape parameter [46]. The parameters α and β in the Weibull function are usually solved by maximum likelihood estimation. Previous studies have shown that maximum likelihood estimation fits well in predicting forest Weibull distribution functions [47]. For this purpose, the formula is as follows: In the above formulas, n is the sample number, x i is the canopy volume of all individual trees, and α and β are calculated by iterative operation.

Watershed Algorithm
The watershed algorithm is based on the CHM for single tree segmentation. First, the canopy boundary is determined by the watershed algorithm, and then the local maximum is found as the tops of the trees [48][49][50][51]. In this study, we used fixed windows with a size of 3 × 3 pixels to detect the local maximum, i.e., the tree tops [52][53][54]. The canopy boundary was obtained by dividing the CHM into the watershed algorithm, which is a mathematical morphology algorithm similar to the simulation of the immersion process. The image was considered to be a 3D topographic map and was composed of the corresponding image in 2D space and the corresponding image gray values as the heights. There were local minimum points in the image. Assuming that a hole forms at each minimum and water is injected into the image, the water should flow through the holes and soak the surface, forming basins. Starting from the lowest minima in the image, the water would gradually flood the basin of the image. As the water level in the basin rises, a dam will form between two areas when the water surfaces from different minimum points converge. When the entire process is finished, the area of each minimum point will be surrounded by a dike, and all the dikes will assemble to form a watershed. Each dike represents the boundary of image segmentation, and the basin is the dividing area [53,55,56]. Finally, the highest value in each basin was chosen as the tree top for each area.

Polynomial Fitting Method
Individual tree segmentation using polynomial fitting was achieved by detecting local maxima in the CHM with a movable dynamic circular filter as the tree tops. Then, each probed tree top was regarded as the polynomial vertex, and the crown was fitted by the polynomial curve according to the morphological growth characteristics of the trees. Because the canopy of a single tree resembles a circle from immediately above the tree, this study used a dynamic circular window (which can be changed according to the crown width of individual trees of different sizes) to detect each tree top [57]; the diameter of the window was the average value of the north-south direction of the single tree and the east-west crown length. When observing adjacent individual trees from the side, the best polynomial fitting condition can be obtained when the edge of the top of the tree and the crown have convex and concave forms in the middle. We can obtain the maximum tree top value and two minimum values at the edge of the tree crown. The horizontal distance between the minimum points of the two crown edges is the crown size [58,59]. In this algorithm, the crown radius was calculated using the following formulas: where X in formula (5) is the horizontal distance and a, b, c, d and e in formula (6) are the coefficients of the fitting curve (a-e are constants). With formula (6), the polynomial of the fitting curve can be derived. A first derivative of 0 represents the extreme point. In formula (7), the two derivatives can be obtained at the extreme value point; when the extreme value is negative, the curve shape is concave. If the extremum is positive, the shape of the curve is convex, and the values of x 1 and x 2 , which correspond to formulas (6) and (7), can finally be calculated. The difference between these two values is the crown width estimated by the algorithm [60]. The ITC delineation approach first detects the local maximum value and marks the CHM through a fixed circular detection window. In this way, the over segmentation phenomenon is reduced; then, the crown tree can be created by the decision tree method at the top of the detected tree [28]. The approach has the following steps: (i) A low-pass filter is applied to the rasterized CHM to smooth the surface and reduce the number of local maxima. (ii) The local maximum value is detected with the circular search window, and the local maximum value is marked. A pixel of the CHM is labeled as the local maximum if its value is greater than all other values in the window. (iii) Each local maximum is labeled as an "initial region" around which a tree crown can grow; the heights of the four neighboring pixels are extracted from the CHM. When these heights differ from the vertical height marked with the local maximum less than the user-defined value (a percentage of the local maximum) and the horizontal distance less than the user-defined horizontal distance, these pixels fall within the growth range from the "initial region"; this step is repeated until the surrounding pixels do not satisfy the above conditions and is then stopped. (iv) The point cloud of the first echo of the LiDAR is added to the previously labeled "initial region." (v) A 2D convex hull is applied to these points, and the resulting polygons become the final ITCS [61].

Point Cloud Segmentation (PCS)
PCS is a method of single tree segmentation based on normalized point cloud data, and its principle is regional growth combined with threshold determination. This method takes full advantage of the phenomenon that a certain distance separates trees, especially tree tops. First, assuming that the tree vertex is the highest point in the normalized point cloud data, a tree is segmented by region growing from the vertex point. Finally, the tree is segmented by iterations of the algorithm until all individual trees in the stand are segmented. The specific process is as follows: (i) Set the horizontal distance as threshold d. (ii) Search the top point of the LiDAR point cloud and choose the horizontal distance between two points larger than d as different tree tops. (iii) Measure the horizontal distance between the target point and the vertex of a tree from top to bottom. If the distance is less than d, it belongs to the given single tree; if the distance is greater than d, it belongs to another single tree. (iv) If the horizontal distance between a target point and the adjacent tree top exceeds d, the target point is segmented as the tree vertex. (v) If the horizontal distance between a target point and the vertices of at least two nearby trees is less than d, the horizontal distances between the target point and the tops of these nearby trees are compared, and the target point belongs to the nearest tree. (vi) Continue this process downward until the bottom points are segmented [62].

Indicators of Individual Tree Segmentation Accuracy
In this study, for each detected tree, the tree position, tree height and crown area were estimated and compared to the corresponding tree in the field. It was considered correct when a detected tree was located within the crown of the field inventory tree. We used the confusion matrix F-score as the accuracy evaluation index, and the calculation methods were as follows [63,64]: where r represents the tree detection rate, p represents the tree detection accuracy, and F represents the overall accuracy taking both omission and commission into consideration. Nt is the number of trees detected by the algorithm and corresponds to the measured trees on the ground. No is the number of trees that cannot be detected by the algorithm but exist on the ground. Nc is the number of trees detected by the algorithm but that do not exist on the ground.

Estimation of Forest Canopy Cover by the Individual Tree Segmentation-Based Method
Firstly, each individual tree crown was extracted from the individual tree segmentation algorithms with highest accuracy. Then, each individual tree crown in a sample plot was added and minus the overlap area to calculated the total area. Finally, the radio of the area of tree crown to the area of the sample plot was calculated as canopy cover.

Estimation of Forest Canopy Cover by the CHM-based Method
In this study, remote sensing image analysis used the filtered and smoothed point cloud data of UAV-LiDAR to interpolate and generate a DEM with a grid resolution of 0.5 m. Then, we generated CHM from height normalized point cloud. According to the traditional method of CC calculation and the tree height of ginkgo in the study area, the predicted forest CC value can be obtained by calculating the ratio of the area of CC pixels in the CHM images above 2 m to the total area of the sample plot. Each covered pixel is recorded as 1, whereas each non-covered pixel is recorded as 0 [20,34]. The total CC area is calculated by counting the number of CC pixels. Finally, the forest CC prediction is obtained by dividing the total canopy cover area by the sample area (this method is shown in Figure 3). The formula is as follows: Canopy Cover = CHM canopy CHM total (11) Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 22 were CHM with grid resolution of 0.5m in three plots; c, f, i were the distribution of canopy cover pixels above 2 m in three plots (black was canopy cover pixels).

Estimation of Forest Canopy Cover by the LiDAR Metrics Statistics Model
In this study, the multiple linear regression (MLR) method was used to construct the linear equation to predict the CC (statistical model method). The dependent variable was the CC calculated from line-intercept method, and the independent variable was a suite of metrics extracted by LiDAR. The independent variables were introduced into the model one by one, and the F test (p value <0.05) was applied to each independent variable. When the previously introduced independent variable was no longer significant due to the subsequent introduction of an independent variable, it was were orthophoto images and individual tree crown widths of three different density plots; (b,e,h) were CHM with grid resolution of 0.5m in three plots; (c,f,i) were the distribution of canopy cover pixels above 2 m in three plots (black was canopy cover pixels).

Estimation of Forest Canopy Cover by the LiDAR Metrics Statistics Model
In this study, the multiple linear regression (MLR) method was used to construct the linear equation to predict the CC (statistical model method). The dependent variable was the CC calculated from line-intercept method, and the independent variable was a suite of metrics extracted by LiDAR. The independent variables were introduced into the model one by one, and the F test (p value <0.05) was applied to each independent variable. When the previously introduced independent variable was no longer significant due to the subsequent introduction of an independent variable, it was deleted, and this process was repeated all the independent variables were filtered out. Finally, the three independent variables with the highest correlation to the CC (the number of independent variables does not increase) were selected for modeling [65,66].

Model Accuracy Evaluation
In this study, the determinant coefficient (R 2 ), root mean square error (RMSE) and relative root mean square error (rRMSE) were used to evaluate the accuracy of the regression model. R 2 reflects the correlation between independent variable and dependent variable. The higher the R 2 value, the higher the correlation between the independent variable and the dependent variable. The RMSE reflects the standard error between measured and predicted values and has a good relationship with the evaluation object value. The smaller the RMSE value is, the better the prediction effect of the model is. The rRMSE is the ratio of the RMSE to the mean predicted value. The rRMSE has nothing to do with the value of the evaluation object itself and reflects the overall prediction accuracy of the model. The smaller the rRMSE is, the higher the prediction accuracy of the model is. The formula is as follows: where x is the ground-measured forest canopy, x is the average of the ground-measured forest CC,x is the forest CC estimated by the model, and n is the number of plots.

Results
Three typical plots were segmented by four algorithms (Figure 4). Figure 4 shows the point clouds of the three plots (a-c), individual trees delineated using the watershed algorithm (d-f), and individual trees delineated using polynomial fitting (g-i); the red frames are individual tree crown boundaries, and the black lines are individual tree stems. The crown boundaries and stems were created by these two algorithms. Figure 4 also shows individual trees delineated using individual tree crown segmentation (j-l) and using PCS (m-o). The different colors represent different individual trees.
A total of 892 individual trees in the study area were detected by the four segment methods. In all plots, PCS had the highest accuracy (F = 0.83), followed by ITCs (F = 0.82), and the watershed algorithm had a higher accuracy (F = 0.79) than polynomial fitting (F = 0.77). The four methods all had the highest segment accuracy in low-stem-density plots (F = 0.87, 0.87, 0.9 and 0.91), followed by medium-stem-density plots (F = 0.79, 0.78, 0.81 and 0.83); the lowest segment accuracy occurred in high-stem-density plots (F = 0.75, 0.72, 0.81 and 0.83) ( Table 3). clouds of the three plots (a-c), individual trees delineated using the watershed algorithm (d-f), and individual trees delineated using polynomial fitting (g-i); the red frames are individual tree crown boundaries, and the black lines are individual tree stems. The crown boundaries and stems were created by these two algorithms. Figure 4 also shows individual trees delineated using individual tree crown segmentation (j-l) and using PCS (m-o). The different colors represent different individual trees.    (8), (9) and (10). Figure 5 shows the individual trees segmented by four segmentation algorithms (i.e., watershed, polynomial fitting, individual tree segmentation and point cloud segmentation) in three stem density plots. The accuracy assessment of three CHM-based segmentation algorithms is shown in Table 4. In all different stem density plots, the F-score of all three algorithms was highest at 0.5 m resolution CHM. As the CHM resolution increased, the counts of individual trees detected by the algorithms increased. Compared with 0.3 m-resolution CHM, the r-values of 0.5 m-and 0.7 m-resolution CHM were lower because the number of individual trees segmented was less than that in the 0.3 m-resolution CHM. The p-values of 0.3 m-resolution CHM were higher than those of 0.5 m-and 0.7 m-resolution CHM.
A sensitivity analysis of PCS is shown in Figure 6. We set three different distance thresholds (i.e., 1 m, 2 m and 3 m), and as the distance threshold increased, the number of individual trees decreased and the tree volume increased ( Figure 6A). In addition, when the distance threshold was 1 m, over-segmentation occurred in all the plots ( Figure 6B). The accuracy assessment of PCS algorithms is shown in Table 5. In all distance thresholds, low-stem-density plots had the highest accuracy, and in all plots, when the distance threshold was 2 m, the PCS had the highest accuracy.    Table 6 presents a comparison and accuracy assessment of forest CC prediction models with different combinations of LiDAR metrics as prediction metrics. As shown in the table, the prediction accuracy was the lowest for only height metrics (R 2 = 0.59, rRMSE = 9.9%) and was improved by adding other metrics. The prediction accuracy based on the height and density metrics was the highest (R 2 = 0.83, rRMSE = 5.6%), followed by that based on the height and canopy volume model (R 2 = 0.78, Remote Sens. 2019, 11, 908 13 of 21 rRMSE = 6.0%). Table 7 presents the variance inflation factor (VIF) values of LiDAR metrics used in the CC models. The result indicated that the independent variables in predictive models were low inter-correlated. Note: the mean of r, p, F can be seen in equation (8), (9) and (10). Figure 5 shows the individual trees segmented by four segmentation algorithms (i.e., watershed, polynomial fitting, individual tree segmentation and point cloud segmentation) in three stem density plots. The accuracy assessment of three CHM-based segmentation algorithms is shown in Table 4. In all different stem density plots, the F-score of all three algorithms was highest at 0.5 m resolution CHM. As the CHM resolution increased, the counts of individual trees detected by the algorithms increased.     Note: the mean of r, p, F can be seen in Equations (8), (9) and (10).  Note: the code of metrics can be seen in Table 2. VIF: variance inflation factor.
The results of the cross validation of the CC values estimated by the metric model method, CHM-based method CHM-based method and individual tree segmentation-based method are shown in Figure 7. The accuracy of the model with only the height metric and field-measured CC was R 2 = 0.48 (rRMSE = 8.3%) (Figure 7a). The accuracy of the model that combined the height metric with the density metric and field-measured CC was R 2 = 0.69, rRMSE = 6.5% (Figure 7b). The accuracy of the model that combined the height metric with the canopy volume metric and field-measured CC was R 2 =0.80, rRMSE = 5.9% (Figure 7c). Figure 7a-c show that the addition of different metrics improved the prediction effect of the model established with only the height metric to varying degrees. CC estimated by the combination of the height metric and density metric had the highest accuracy; CC estimated by the combination of the height metric and canopy volume metric had the second highest accuracy. Three canopy pixels with different heights were extracted from a CHM image with a grid resolution of 0.5 m × 0.5 m to estimate CC. The accuracy of CC estimation via the extraction of canopy pixels above 1.3 m was R 2 = 0.92 (rRMSE = 5.9%) (Figure 7d). The accuracy of CC estimation via the extraction of canopy pixels above 2 m was R 2 = 0.94 (rRMSE = 5.4%) (Figure 7e). The accuracy of CC estimation via the extraction of canopy pixels above 5 m was R 2 = 0.87 (rRMSE = 8.1%) (Figure 7f).
The canopy crown extracted from PCS algorithm in three different distance thresholds to estimate CC (ITSM). The accuracy of CC estimation by the distance threshold at 1m was R 2 = 0.54 (rRMSE = 11.2%) (Figure 7g). The accuracy of CC estimation by the distance threshold at 2m was R 2 = 0.92 (rRMSE = 3.5%) (Figure 7h). The accuracy of CC estimation by the distance threshold at 1m was R 2 = 0.74 (rRMSE = 9.3%) (Figure 7i). The accuracy of CC estimation via the ITSM was higher than that via the CHMM and the SMM, and the highest accuracy was obtained by distance threshold at 2 m (R 2 = 0.92, rRMSE = 6.4%). The results of the cross validation of the CC values estimated by the metric model method, CHM-based method CHM-based method and individual tree segmentation-based method are shown in Figure 7. The accuracy of the model with only the height metric and field-measured CC was R 2 = 0.48 (rRMSE = 8.3%) (Figure 7a). The accuracy of the model that combined the height metric with the density metric and field-measured CC was R 2 = 0.69, rRMSE = 6.5% (Figure 7b). The accuracy of the model that combined the height metric with the canopy volume metric and field-measured CC was R 2 =0.80, rRMSE = 5.9% (Figure 7c). Figure 7a, b and c show that the addition of different metrics improved the prediction effect of the model established with only the height metric to varying

Discussion
In this study, the CC of ginkgo-planted forests in northern Jiangsu Province was estimated using UAV-LiDAR data. We selected 23 plots of ginkgo (892 individual trees with three different densities) and separated individual trees by the watershed algorithm, polynomial fitting algorithm, individual tree crown segmentation method and PCS method. Calculating the detection rate, accuracy rate and overall accuracy index values of individual trees allowed the segmentation accuracies of the four algorithms in sample plots with different densities were compared. We also calculated the CC at the sample plot scale using the method with the highest segmentation accuracy as validation and then estimated the CC using the CHM extraction approach (CHM-based method) and metric model method (statistic model method). The results showed that the four segmentation methods had the highest segmentation accuracy ( Figure 5 and Table 3) when the plot density was low, with F values of 0.87, 0.87, 0.90 and 0.91. With the increase in plot density, the segmentation accuracy of the four methods showed a decreasing trend. In low-density plots, the spacing between individual tree crowns is relatively large and there is little overlap between trees. Therefore, when using the CHM-based method, the local maximum is usually only the top of a tree, which reduces over segmentation and missed segmentation [52,67]. In a comparison of the segmentation accuracy of the PCS method and the CHM individual tree crown segmentation method, the PCS method had the highest accuracy, and the difference between the two methods increased with increasing sample density. This result is similar to the results of Tao et al. (2014) [68]. Taking a coniferous forest in California, USA, as the research object, Tao et al. compared the segmentation accuracy of PCS with that of a watershed algorithm. The results showed that the PCS accuracy (average F value of 0.9) was higher than the watershed algorithm accuracy (average F value of 0.81). Especially for small trees beneath the canopy of large trees, the PCS method was more powerful [28,61,62], as this method is directly based on point cloud data, without interpolation and other processing, thus avoiding information loss and other conditions impacting accuracy [69].
In a comparison of the three types of CHM-based segmentation methods, the watershed algorithm had the highest error rate (lowest p value) because the watershed algorithm first divides the crown by the CHM, and the local maximum value of the crown is detected as the tree tops. This process often divides a single tree into several trees, resulting in multiple phenomena [49,70]. The detection rate r obtained by the polynomial fitting method was the lowest. The polynomial fitting method is based on a dynamic circular detection window, which considers the local maximum on the CHM as the tree tops. Then, based on the extreme points for crown segmentation, the detection window radius is determined according to the linear relationship between tree height and crown width. The window radius used in this study was based on Popescu et al. (2004) [53]. The detection rate was low because the target of the study by Popescu et al. was a Pinus massoniana forest in North America, which has a different structure than the ginkgo planted forests in this study. Popescu et al. (2004) and Huo et al. (2015) [53,60] reported a higher precision in coniferous forests than in this study because their research objects were coniferous tree species. Gingko trees have an intermediate height and relatively narrow crown, but in some cases, the crown width does not exhibit this shape and the extremum points are not at the top or at both sides of the crown, decreasing the detection of the extremum points and affecting the r-value. The CHM individual tree crown segmentation method had the highest segmentation accuracy. When using the CHM individual tree crown segmentation method to segment a single tree, the over segmentation phenomenon is reduced by determining the detected local maximum. At the same time, the decision tree growth principle is used to classify the points in each pixel near the local maximum to form crown width, which effectively increases the accuracy of individual tree segmentation. Therefore, in this study, the individual tree crown segmentation method had the highest overall accuracy among the CHM-based methods [61].
Different CHM raster resolutions affect the detection of the local maximum. Therefore, we changed the raster resolution of CHM to analyze the accuracy of three types of CHM individual tree segmentation algorithms. As shown in Table 4, varying the resolution of the CHM grid had a great influence on segmentation accuracy. The results show that the overall segmentation accuracy of the three algorithms was the highest when the CHM raster resolution was 0.5 m × 0.5 m. Similar to the results of Koch et al. (2006) and Chen et al. (2006) [49,71], when the CHM resolution was too low, local maxima were divided into several parts, so the number of segmentations increased and the p value decreased; when the CHM resolution was too high, leakage occurred and the number of local maxima was reduced, affecting the r-value [64]. Sensitivity analysis of key parameters of PCS was carried out by changing the distance threshold because the basic principle of this method is to determine the segmentation based on the distance between point clouds. The results ( Figure 6 and Table 5) show that when the distance threshold was 2 m, the segmentation accuracy was the highest. Li et al. (2012) set the default value of the distance threshold to 2 m when PCS was proposed and yielded results similar to those of the present study [62]. When the selected distance threshold was less than 2 m, such as 1 m, the over-segmentation phenomenon increased and the p value was low; however, when the selected distance threshold was greater than 2 m, the single tree crowns increased, the missing segmentation phenomenon increased, and the r-value was low.
In this study, we used the statistical model method to estimate canopy density based on the prediction of forest canopy density only through the height metric. The prediction accuracy of CC was improved by introducing the density metric and canopy volume metric step by step (Figure 7), similar to previous research results. Naesset et al. used airborne LiDAR to retrieve the forest structure parameters of different forest ages in Norway. By extracting canopy height and canopy density metrics from LiDAR point clouds, the results showed that the estimation accuracy resulting from a combination of height with density metrics was better than that resulting from the use of only height or density metrics [22]. The model established by the combination of height and density metrics had the highest precision (R 2 = 0.80, rRMSE = 5.9%). Compared with the forest canopy accuracy (R 2 = 0.48, rRMSE = 8.3%) estimated by using only the height metric, the prediction accuracy was improved by adding the canopy density metric; the canopy height metric was larger than the height of vegetation in a proportion of all laser spots, and the density metric selected in this study had 50% and 70% quantile variables (as shown in Tables 2 and 6), which can better represent upper and middle canopy structures. Considering the growth of ginkgo planted forests, the branches of ginkgo were concentrated in the middle and upper parts of the trunk, so the point cloud were mainly concentrated in the middle and upper parts of the canopy.
When the canopy density was predicted by the CHM-based method, the number of canopy pixels was counted by pixel decomposition based on the CHM generated by LiDAR point cloud interpolation. In the process of generating the CHM from point clouds by grid processing, the ground point clouds usually appeared to be closed after interpolation. Therefore, the predicted value of the CHM was often higher than the measured value [34], which is similar to the results of Lee and Lucas (2007) and Korhonen (2011) [20,72], as shown in Figure 7. The fitting precision of the canopy density estimated by the CHM method was higher (R 2 = 0.94, rRMSE = 5.4%) than that estimated by the metrics statistical model method. The fitting precision of the canopy density was the highest (R 2 = 0.94, rRMSE = 5.4%) based on the CHM-based method. This result occurred because the density of the UAV-LiDAR used in this study was approximately 159 pts·m −2 , so the point density was high, which can describe the forest vertical structure in relatively high detail and is conducive to improving the prediction effect of forest CC.

Conclusions
The accurate acquisition of canopy cover information for ginkgo-planted forests is essential for the characterization of forest ecosystems. In this study, we used UAV-based high-density LiDAR data to detect individual tree and estimate the CC of ginkgo planted forests in China. The results demonstrate that, for individual tree segmentation in the total plot, PCS had the highest accuracy (F = 0.83), followed by ITCS (F = 0.82). The watershed algorithm had a higher accuracy (F = 0.79) than the polynomial fitting algorithm (F = 0.77). With the increase in stem densities, the segmentation accuracy of the four algorithms was reduced. For the three CHM-based methods, the segmentation accuracy was the highest in the CHM with a 0.5-m resolution, and for PCS, the segmentation accuracy was the highest with the 2-m distance threshold. All three of the CC estimation methods had a relatively high accuracy. For the ITSM, the canopy cover estimation accuracy was the highest with PCS algorithms when the distance threshold was 2 m (R 2 = 0.92, rRMSE = 3.5%). For CHM-based method estimations, the accuracy was the highest with a CHM value above 2 m (R 2 = 0.94, rRMSE = 5.4%). For statistic model method estimations, the combo models (R 2 = 0.80, rRMSE = 5.9%) with both height and density metrics had a higher performance than the models fitted by height metrics only (R 2 = 0.48, rRMSE = 8.3%) and the models fitted by height metrics and canopy volume metrics (R 2 = 0.69, rRMSE = 6.5%).
Author Contributions: X.W. and X.S. analyzed the data and wrote the paper. L.C. helped in project and study design, paper writing and analysis. G.W. and F.C. helped with field work and data analysis.