UAV Based Estimation of Forest Leaf Area Index (LAI) through Oblique Photogrammetry

: As a key canopy structure parameter, the estimation method of the Leaf Area Index (LAI) has always attracted attention. To explore a potential method to estimate forest LAI from 3D point cloud at low cost, we took photos from different angles of the drone and set ﬁve schemes (O (0 ◦ ), T15 (15 ◦ ), T30 (30 ◦ ), OT15 (0 ◦ and 15 ◦ ) and OT30 (0 ◦ and 30 ◦ )), which were used to reconstruct 3D point cloud of forest canopy based on photogrammetry. Subsequently, the LAI values and the leaf area distribution in the vertical direction derived from ﬁve schemes were calculated based on the voxelized model. Our results show that the serious lack of leaf area in the middle and lower layers determines that the LAI estimate of O is inaccurate. For oblique photogrammetry, schemes with 30 ◦ photos always provided better LAI estimates than schemes with 15 ◦ photos (T30 better than T15, OT30 better than OT15), mainly reﬂected in the lower part of the canopy, which is particularly obvious in low-LAI areas. The overall structure of the single-tilt angle scheme (T15, T30) was relatively complete, but the rough point cloud details could not reﬂect the actual situation of LAI well. Multi-angle schemes (OT15, OT30) provided excellent leaf area estimation (OT15: R 2 = 0.8225, RMSE = 0.3334 m 2 /m 2 ; OT30: R 2 = 0.9119, RMSE = 0.1790 m 2 /m 2 ). OT30 provided the best LAI estimation accuracy at a sub-voxel size of 0.09 m and the best checkpoint accuracy (OT30: RMSE [H] = 0.2917 m, RMSE [V] = 0.1797 m). The results highlight that coupling oblique photography and nadiral photography can be an effective solution to estimate forest LAI.


Introduction
Leaf area index (LAI) is defined as half of the total leaf area per unit ground surface area [1][2][3] and is a critical parameter in many forest process models. As a key canopy structure feature, LAI controls the biophysical processes of the forest (photosynthesis, respiration, transpiration, carbon cycle and precipitation interception, etc.) [4]. Therefore, it is important to estimate LAI quickly and accurately.
Many field-based methods have been developed for LAI distribution estimation, but they are still time-consuming and labor-intensive work, cause permanent damage to the forest [5,6] or are affected by many factors (such as leaf distribution, canopy shape and sampling scale) [7]. Multispectral data from satellite or unmanned aerial vehicle (UAV) have been widely used in the inversion of forest LAI in the past. However, two-dimensional images based on optical passive remote sensing often ignore the leaf area of the middle and lower canopy, which will cause inaccurate estimation of LAI [8,9]. Point clouds data can directly reflect the three-dimensional structure of the canopy, which is widely used in

UAV Parameters and Flight-Scheme Design
DJI Phantom 4 drone equipped with multispectral sensors (Phantom 4, DJ-Innovations, Shenzhen, China) was used to collect aerial images. The DJI Phantom 4 multispectral version is a vertical-lift quadrotor drone with real-time kinematics (RTK). The network RTK produced an average accuracy of 1.5 cm horizontally and 1 cm vertically. The maximum photo resolution was 1600 × 1300 px (4:3.25); the ground sampling distance was (H/18.9) cm/pixel (H is flying height). The effective lens field of vision was 62.7°, and focal length was 5.74 mm. The rotation angle of the PTZ was pitch: −90° to + 30°.
To obtain more detailed forest canopy information, we used a visible light camera on a UAV to take photos in the same route with different lens angles. Considering that an excessively large lens tilt angle would result in a lower homology point-matching relationship, we acquired two oblique image sets (15°, 30°) and one nadir photo set (0°). The overlap rate of the UAV flight course heading and side direction was 85% at a flying height of 120 m. Taking into account the particularity of the view of oblique photography, we used an orthogonal flight to obtain oblique photos (Figure 2), which was implemented in two phases (N-S and E-W). Despite the field of view of the nadir photo being the same in all directions, orthogonal flight plans are also used to obtain nadir photos to ensure that the results are comparable. RTK was used to accurately locate the pre-arranged ground control points (GCP) around the 72 sample plots to determine the coordinate accuracy of the derived point cloud and the point cloud coordinate registration ( Figure 1).

UAV Parameters and Flight-Scheme Design
DJI Phantom 4 drone equipped with multispectral sensors (Phantom 4, DJ-Innovations, Shenzhen, China) was used to collect aerial images. The DJI Phantom 4 multispectral version is a vertical-lift quadrotor drone with real-time kinematics (RTK). The network RTK produced an average accuracy of 1.5 cm horizontally and 1 cm vertically. The maximum photo resolution was 1600 × 1300 px (4:3.25); the ground sampling distance was (H/18.9) cm/pixel (H is flying height). The effective lens field of vision was 62.7 • , and focal length was 5.74 mm. The rotation angle of the PTZ was pitch: −90 • to + 30 • .
To obtain more detailed forest canopy information, we used a visible light camera on a UAV to take photos in the same route with different lens angles. Considering that an excessively large lens tilt angle would result in a lower homology point-matching relationship, we acquired two oblique image sets (15 • , 30 • ) and one nadir photo set (0 • ). The overlap rate of the UAV flight course heading and side direction was 85% at a flying height of 120 m. Taking into account the particularity of the view of oblique photography, we used an orthogonal flight to obtain oblique photos (Figure 2), which was implemented in two phases (N-S and E-W). Despite the field of view of the nadir photo being the same in all directions, orthogonal flight plans are also used to obtain nadir photos to ensure that the results are comparable. RTK was used to accurately locate the pre-arranged ground control points (GCP) around the 72 sample plots to determine the coordinate accuracy of the derived point cloud and the point cloud coordinate registration ( Figure 1).

Methods
The specific methods followed in this study are summarized in Figure 3, which is mainly divided into three parts: (1) 3D-reconstruction process; (2) point-cloud data processing; (3) LAI estimation. The code for voxelization model and LAI fast extraction is detailed in the attachment.

Methods
The specific methods followed in this study are summarized in Figure 3, which is mainly divided into three parts: (1) 3D-reconstruction process; (2) point-cloud data processing; (3) LAI estimation. The code for voxelization model and LAI fast extraction is detailed in the attachment.

Methods
The specific methods followed in this study are summarized in Figure 3, which is mainly divided into three parts: (1) 3D-reconstruction process; (2) point-cloud data processing; (3) LAI estimation. The code for voxelization model and LAI fast extraction is detailed in the attachment.  SFM-MVS workflow can couple photos well from different angles to restore 3D structure and consists of two steps: (1) structure-from-motion (SFM) and (2) multi-view stereo (MVS). The SFM algorithm is a motion photogrammetric technique used to reconstruct 3D point clouds from high-overlapping images [26], the derived point clouds of which are sparse. In the second step, multi-view stereo (MVS) reconstructs dense point clouds based on sparse point clouds and image-feature-matching results, camera position, attitude and internal parameters calculated by SFM [27]. To explore the influence of UAV images taken from different angles on the 3D reconstruction of Masson pine forest, we set up five solutions for point-cloud reconstruction through the forms of single angle and multiple angles: The images were processed on Pix4 D (version 4.3, Pix4 D, Prilly, Switzerland), which can automatically construct 3D structure based on the SFM-MVS algorithm. The SFM process of Pix4 D starts with image-feature detection. In this critical step, homology-matching points are obtained, which are later used for generating image correspondence. Next, a sparse point cloud is generated based on iterative bundle adjustment (BA) [28]. The image matching result is the decisive factor to decide whether the surface model is complete or not. After using different parameters for image matching, we found that when the number of adjacent images is set to 4 or more and the self-calibration scheme is set to "precise geolocation and direction", an effective homologous point matching result is derived ( Figure 4). Meanwhile, geometric verification matching can provide accurate and effective image matching. In the phase of generating point cloud and encrypting point cloud, we set the density of point cloud and image scale to the highest to ensure the integrity of canopy structure. SFM-MVS workflow can couple photos well from different angles to restore 3D structure and consists of two steps: (1) structure-from-motion (SFM) and (2) multi-view stereo (MVS). The SFM algorithm is a motion photogrammetric technique used to reconstruct 3D point clouds from high-overlapping images [26], the derived point clouds of which are sparse. In the second step, multi-view stereo (MVS) reconstructs dense point clouds based on sparse point clouds and image-feature-matching results, camera position, attitude and internal parameters calculated by SFM [27]. To explore the influence of UAV images taken from different angles on the 3D reconstruction of Masson pine forest, we set up five solutions for point-cloud reconstruction through the forms of single angle and multiple angles: (1) O (0°), (2) T15 (15°), (3) T30 (30°), (4) OT15 (0° and 15°) and (5) OT30 (0° and 30°). The images were processed on Pix4 D (version 4.3, Pix4 D, Prilly, Switzerland), which can automatically construct 3D structure based on the SFM-MVS algorithm. The SFM process of Pix4 D starts with image-feature detection. In this critical step, homology-matching points are obtained, which are later used for generating image correspondence. Next, a sparse point cloud is generated based on iterative bundle adjustment (BA) [28]. The image matching result is the decisive factor to decide whether the surface model is complete or not. After using different parameters for image matching, we found that when the number of adjacent images is set to 4 or more and the self-calibration scheme is set to "precise geolocation and direction", an effective homologous point matching result is derived ( Figure  4). Meanwhile, geometric verification matching can provide accurate and effective image matching. In the phase of generating point cloud and encrypting point cloud, we set the density of point cloud and image scale to the highest to ensure the integrity of canopy structure.

Point-Cloud Data Processing
To compare the LAI estimates of different schemes (O, OT15, OT30, T15, T30) with comparability, we used GCP coordinates pre-arranged around the sample site [29] for point-cloud registration, which was performed in CloudCompare (http://cloudcompare.org/ (accessed on 17 July 2020)). Redundant ground point clouds overestimate the final LAI estimation result, so all ground point clouds and noise should be eliminated before constructing the voxel model. Even though the automatic point cloud classification function of the Pix4 d software is very convenient, the classification results often have wrong estimates and cause the separated ground points to be too sparse. Therefore, we chose another method (Cloth Simulation Filter) to extract ground point clusters. Cloth

Point-Cloud Data Processing
To compare the LAI estimates of different schemes (O, OT15, OT30, T15, T30) with comparability, we used GCP coordinates pre-arranged around the sample site [29] for pointcloud registration, which was performed in CloudCompare (http://cloudcompare.org/ (accessed on 17 July 2020)). Redundant ground point clouds overestimate the final LAI estimation result, so all ground point clouds and noise should be eliminated before constructing the voxel model. Even though the automatic point cloud classification function of the Pix4 d software is very convenient, the classification results often have wrong estimates and cause the separated ground points to be too sparse. Therefore, we chose another method (Cloth Simulation Filter) to extract ground point clusters. Cloth Simulation Filter (CSF) is a tool that can extract ground points by simply setting a few parameters [30,31]. The cloth resolution refers to the size of cloth used to cover the terrain. Excessive cloth resolution will result in rough ground point clouds. The CSF process was executed in Cloudcompare using the following parameters: the general parameter was set to a steep slope, the cloth resolution was 0.5 m, the maximum iteration was 500 and the classification threshold was 0.1 m. Then, we used the simple kriging method in the geostatistical module of ArcGIS (version 10.2, ESRI, Redlands, CA, USA) to interpolate the ground point clouds to obtain the DEM data. Finally, the Z-value of each point coordinate was subtracted from the height value of the corresponding position in the DEM to obtain the height value of the point clouds relative to the ground. The derived point clouds collection also included forest canopy, ground vegetation, soil, etc. Notably, the collection of non-canopy point clouds can lead to the overestimation of leaf area in the lower part of the stand. Therefore, it was necessary to eliminate the non-canopy point-cloud collections. The height of the Masson pine branches in the study area was generally about 2 m, so we used 2 m as a dividing limit to delete the non-canopy point-cloud collection. The elimination of noise was performed in MATLAB software (version 2018 a, MathWorks, Natick, MA, USA) according to the following specific steps: 1. Use KNNSearch to retrieve the 30 nearest points of each point. 2. Calculate the average distance between each point and the 30 nearest neighbors and the mean and standard deviation of these average distances. If the difference between the average distance from a point to the 30 nearest neighbors and the mean is not within 1 times the standard deviation, it is considered a noise point and deleted.

LAI Calculation Method
Point-cloud data derived from UAV images do not reflect details at a single-tree level well. Thus, we followed the method developed by Hosoi and Omasa [32] to construct a voxel model. Using the voxel model to calculate LAI can avoid recalculating the unilateral leaf area of the same leaf. The point-cloud coordinates was converted in MATLAB into voxelized array, which contained both a 0.5 × 0.5 × 0.5 m voxelized model and a subvoxelized model. To compare the leaf area estimation capabilities of point clouds derived from different schemes at different heights of the canopy, we extracted the leaf area density (LAD) distribution, defined as the sum of the unilateral leaf area per unit volume [33], at different heights of the canopy from the ground from the voxel model. Then, the derived LAI was calculated based on the principle that the integral of LAD in the vertical direction is equal to the LAI. Clearings in the forest that did not contain point clouds could not be retrieved with the voxelized array, which would cause leaf-area overestimation. The subvoxel size was the most critical factor in determining the estimated value of leaf area. To explore the appropriate sub-voxel size of the derived point clouds for leaf area estimation, we set up ten 0.06 to 0.15 m side-length sub-voxels in steps of 0.01 m, according to the point-cloud density in the report from Pix4 D. For these areas, we added the coordinates of the blank voxels in the voxelized array and set the LAD value of these blank voxels to 0. Then, the forest stands were stratified according to the height, and the forest LAD was calculated using a stratification interval of 0.5 m. We used the average LAD values of all voxels in each 0.5 m layer within each plot range as the LAD value of the layer.

LAI Field Measurement
To obtain the optimal sub-voxel size interval for leaf area calculation, we adopted the leaf area index LAI of each plot. The LAI-2200 canopy analyzer-equipped with 90 • covering cap-was used to measure the LAI values of 72 plots on a cloudy day without rain. For each sample plot, we took the center and four corners as measurement sites ( Figure 5). When measuring at the corner of the sample plot, the probe of LAI-2200 points to the center, and when measuring at the center of the sample plot, the probe of LAI-2200 points to the corner ( Figure 5), and finally, the average value of eight measurement results is taken as the LAI value of sample plot. The measurement was repeated three times on each sample plot according to the same method. The outliers were removed, and the average of the remaining data were taken as the final LAI value of each plot. to the corner ( Figure 5), and finally, the average value of eight measurement results is taken as the LAI value of sample plot. The measurement was repeated three times on each sample plot according to the same method. The outliers were removed, and the average of the remaining data were taken as the final LAI value of each plot. Although the LAI obtained by the LAI-2200 is actually the effective leaf area index (LAIe), many studies use it to replace LAI because it is difficult to remove the influence of non-leaf components [34]. The LAI values derived from the voxelized models of different sub-voxel sizes and the LAI values measured on-site using the LAI 2200 canopy analyzer were used for linear fitting. Then, the interval of the sub-voxel size most suitable for LAI estimation for each scheme was judged by comparing the slope (a), intercept (b), correlation coefficient (R 2 ) of the fitted linearity and the root mean square error (RMSE) between the two sets of data. Generally speaking, the derived LAI is closer to LAIe when the R 2 and slope (a) are closer to 1 and the intercept (b) and RMSE are closer to 0.

Point Clouds Coordinates Accuracy Assessments
All the point clouds collections of the five schemes were generated from photos combining of different angles without using GCP. The coordinates of GCP collected from derived point clouds were compared with the coordinates collected by RTK ( Figure 6). On the premise that GCP coordinates are not incorporated in the bundle adjustment (BA),  Although the LAI obtained by the LAI-2200 is actually the effective leaf area index (LAIe), many studies use it to replace LAI because it is difficult to remove the influence of non-leaf components [34]. The LAI values derived from the voxelized models of different sub-voxel sizes and the LAI values measured on-site using the LAI 2200 canopy analyzer were used for linear fitting. Then, the interval of the sub-voxel size most suitable for LAI estimation for each scheme was judged by comparing the slope (a), intercept (b), correlation coefficient (R 2 ) of the fitted linearity and the root mean square error (RMSE) between the two sets of data. Generally speaking, the derived LAI is closer to LAIe when the R 2 and slope (a) are closer to 1 and the intercept (b) and RMSE are closer to 0.

Point Clouds Coordinates Accuracy Assessments
All the point clouds collections of the five schemes were generated from photos combining of different angles without using GCP. The coordinates of GCP collected from derived point clouds were compared with the coordinates collected by RTK ( Figure 6). On the premise that GCP coordinates are not incorporated in the bundle adjustment (BA),

LAI Retrieval from Five Schemes
The derived LAI values were obtained from the 10-voxel models generated by each of the five schemes ( Figure 7). Overall, the LAI derived from the five schemes increased with the increase in the sub-voxel size. For different voxelization models, OT30 always provided LAI estimates higher than O, OT15, T15 and T30. The LAI estimates provided by T15 and T30 were very close, and O always provided the lowest LAI estimate among the five schemes. For the five schemes, the sub-voxel sizes corresponding to the slope a, intercept b, the highest R 2 and the lowest RMSE of the linear fitting were explored. O, T15 and T30 give the best LAI estimates at the sub-voxel size of 0.12 m, while OT30 gives the best LAI estimates at the sub-voxel size of 0.09 m and OT15 gives the best LAI estimates at the subvoxel size of 0.10 m (Table 1). This result is consistent with the RMSE result (sub-voxel

LAI Retrieval from Five Schemes
The derived LAI values were obtained from the 10-voxel models generated by each of the five schemes ( Figure 7). Overall, the LAI derived from the five schemes increased with the increase in the sub-voxel size. For different voxelization models, OT30 always provided LAI estimates higher than O, OT15, T15 and T30. The LAI estimates provided by T15 and T30 were very close, and O always provided the lowest LAI estimate among the five schemes.

LAI Retrieval from Five Schemes
The derived LAI values were obtained from the 10-voxel models generated by each of the five schemes ( Figure 7). Overall, the LAI derived from the five schemes increased with the increase in the sub-voxel size. For different voxelization models, OT30 always provided LAI estimates higher than O, OT15, T15 and T30. The LAI estimates provided by T15 and T30 were very close, and O always provided the lowest LAI estimate among the five schemes. For the five schemes, the sub-voxel sizes corresponding to the slope a, intercept b, the highest R 2 and the lowest RMSE of the linear fitting were explored. O, T15 and T30 give the best LAI estimates at the sub-voxel size of 0.12 m, while OT30 gives the best LAI estimates at the sub-voxel size of 0.09 m and OT15 gives the best LAI estimates at the subvoxel size of 0.10 m (Table 1). This result is consistent with the RMSE result (sub-voxel For the five schemes, the sub-voxel sizes corresponding to the slope a, intercept b, the highest R 2 and the lowest RMSE of the linear fitting were explored. O, T15 and T30 give the best LAI estimates at the sub-voxel size of 0.12 m, while OT30 gives the best LAI estimates at the sub-voxel size of 0.09 m and OT15 gives the best LAI estimates at the sub-voxel size of 0.10 m (Table 1). This result is consistent with the RMSE result (sub-voxel size is 0.  (Figure 8 and Table 2). Although the three single-angle schemes also provide good LAI estimation accuracy, the most suitable sub-voxel size is larger than that of OT15 and OT30. This shows that the completeness of O, T15 and T30 canopy information is lower than that of OT15 and OT30 (Figure 8).  Each of the five schemes had a good linear fitting relationship with the LAI value measured on the field under the sub-voxel size that is most suitable for estimating LAI. The results of OT30 were stably distributed around Y = X, and its accuracy verification results are the best among all the schemes (Figure 9). The LAI values derived from the other four schemes had a relatively high degree of dispersion with the LAI values measured on field. Generally speaking, between the same kind of schemes, the LAI estimation provided by the multi-angle scheme was better than with the single-angle scheme (OT30 better than T30), and the LAI derived from adding the 30° photo was better than the one with the 15° photo (OT30 better than OT15; T30 better than T15) (Figure 9). However, the best LAI estimation accuracy of OT15 was worse than the other four schemes, which may be related to the poor image matching between the oblique photos of OT15. Each of the five schemes had a good linear fitting relationship with the LAI value measured on the field under the sub-voxel size that is most suitable for estimating LAI. The results of OT30 were stably distributed around Y = X, and its accuracy verification results are the best among all the schemes (Figure 9). The LAI values derived from the other four schemes had a relatively high degree of dispersion with the LAI values measured on field. Generally speaking, between the same kind of schemes, the LAI estimation provided by the multi-angle scheme was better than with the single-angle scheme (OT30 better than T30), and the LAI derived from adding the 30 • photo was better than the one with the 15 • photo (OT30 better than OT15; T30 better than T15) (Figure 9). However, the best LAI estimation accuracy of OT15 was worse than the other four schemes, which may be related to the poor image matching between the oblique photos of OT15.  Each of the five schemes had a good linear fitting relationship with the LAI value measured on the field under the sub-voxel size that is most suitable for estimating LAI. The results of OT30 were stably distributed around Y = X, and its accuracy verification results are the best among all the schemes (Figure 9). The LAI values derived from the other four schemes had a relatively high degree of dispersion with the LAI values measured on field. Generally speaking, between the same kind of schemes, the LAI estimation provided by the multi-angle scheme was better than with the single-angle scheme (OT30 better than T30), and the LAI derived from adding the 30° photo was better than the one with the 15° photo (OT30 better than OT15; T30 better than T15) (Figure 9). However, the best LAI estimation accuracy of OT15 was worse than the other four schemes, which may be related to the poor image matching between the oblique photos of OT15.

Leaf Area Distribution in Vertical Directions under the Same Sub-Voxel Size
To explore the structural integrity of the point cloud of five schemes at different heights and its influence on the overall LAI estimation of the canopy, the leaf area distribution at different heights of the canopy, namely the leaf area density (LAD), was compared under a same sub-voxel size. Then, scheme O was taken as a reference to further compare the LAD increments of the remaining four schemes at different heights. The vertical LAD distribution trends of the five schemes were very consistent with the morphological characteristics of the masson pine forest, but there was a large difference between the LAD values of the five schemes ( Figure 10). The LAD estimate provided by the single tilt angle schemes (T15, T30) was always lower than the multiple angle schemes (OT15, OT30). In the middle and lower part of the canopy, T15, T30, OT15 and OT30 all gave leaf area estimates higher than O, but the LAD increments of these four schemes relative to O were not consistent in stands with different LAI values. Here, we divide the sample plots into four groups according to the value of LAI and calculate the LAD values and LAD increments of these four plans at different heights of the canopy relative to the O schemes ( Figure 10). For the middle and lower part of the canopy, the LAD increment of the four schemes relative to O in the area where the LAI of 0.5-1.5 was more significant than that in the area with a LAI of 1.5-2.5. In the area where the LAI was 0.5-1.5, the LAD increment of the scheme with 30 • photos was higher than that of the scheme with 15 • photos, was not significant in the area where the LAI interval was 1.5-2. This shows that the acquisition for the lower part of the canopy from UAV gradually weakened with the increase of LAI and the decrease in the UAV lens angle. It is worth noting that T15 and T30 all showed a large decrease in LAD relative to O at the upper part of the canopy. This suggests that the point cloud density of O in the upper part of the canopy was higher than T15 and T30 and close to OT15 and OT30 ( Figure 10). logical characteristics of the masson pine forest, but there was a large difference between the LAD values of the five schemes ( Figure 10). The LAD estimate provided by the single tilt angle schemes (T15, T30) was always lower than the multiple angle schemes (OT15, OT30). In the middle and lower part of the canopy, T15, T30, OT15 and OT30 all gave leaf area estimates higher than O, but the LAD increments of these four schemes relative to O were not consistent in stands with different LAI values. Here, we divide the sample plots into four groups according to the value of LAI and calculate the LAD values and LAD increments of these four plans at different heights of the canopy relative to the O schemes ( Figure 10). For the middle and lower part of the canopy, the LAD increment of the four schemes relative to O in the area where the LAI of 0.5-1.5 was more significant than that in the area with a LAI of 1.5-2.5. In the area where the LAI was 0.5-1.5, the LAD increment of the scheme with 30° photos was higher than that of the scheme with 15° photos, was not significant in the area where the LAI interval was 1.5-2. This shows that the acquisition for the lower part of the canopy from UAV gradually weakened with the increase of LAI and the decrease in the UAV lens angle. It is worth noting that T15 and T30 all showed a large decrease in LAD relative to O at the upper part of the canopy. This suggests that the point cloud density of O in the upper part of the canopy was higher than T15 and T30 and close to OT15 and OT30 ( Figure 10). Figure 10. In areas with different LAI levels, the leaf area density (LAD) distribution and the LAD increment relative to O of five schemes at different heights from ground (sub-voxel size is 0.09 m). a is the area with LAI interval 0.5-1, b is the area with LAI interval 1-1.5, c is the area with LAI interval 1.5-2, d is the area with LAI interval 2-2.5. Figure 10. In areas with different LAI levels, the leaf area density (LAD) distribution and the LAD increment relative to O of five schemes at different heights from ground (sub-voxel size is 0.09 m). a is the area with LAI interval 0.5-1, b is the area with LAI interval 1-1.5, c is the area with LAI interval 1.5-2, d is the area with LAI interval 2-2.5.

Three-Dimensional Reconstruction
Our research result was that coupling oblique photos and nadir photos can effectively improve the vertical coordinate accuracy of the point cloud. However, it is worth noting that under the premise of not using GCP, the height coordinate error provided by the scheme with nadir-photos-only reached 10.8 m (Figure 11), which is quite different from the results of other researchers [35][36][37]. The specific reason for this result is not clear, but it may be related to the camera calibration result or ground type, which can be further explored in future research. Even if there is still a coordinate error higher than the image resolution after point clouds registration, this degree of checkpoint error will not have much impact on the calculation of the canopy leaf area index on the plot scale. A large number of GCP can significantly improve the accuracy of point cloud coordinates [38]. However, setting a large number of GCPs in a complex forest scene is very difficult, and the low-cost advantages of the SFM-MVS process decrease as the GCP deployment points increase. Many studies have pointed out that coupling oblique photographs and nadiral photographs with RTK (without GCP) can also restore the surface model well in some complex scenes, which is consistent with the results of this study [39,40]. Restricted by the terrain and the coverage of the base station, RTK will produce large errors when applied to the coordinate acquisition of GCP under a forest that is heavily shaded by trees. Hence, based on our research results, in forest scenes, we recommend combining the nadir photo and the 30 • tilted photo of the drone to obtain the canopy 3D point cloud.
but it may be related to the camera calibration result or ground type, which can be further explored in future research. Even if there is still a coordinate error higher than the image resolution after point clouds registration, this degree of checkpoint error will not have much impact on the calculation of the canopy leaf area index on the plot scale. A large number of GCP can significantly improve the accuracy of point cloud coordinates [38]. However, setting a large number of GCPs in a complex forest scene is very difficult, and the low-cost advantages of the SFM-MVS process decrease as the GCP deployment points increase. Many studies have pointed out that coupling oblique photographs and nadiral photographs with RTK (without GCP) can also restore the surface model well in some complex scenes, which is consistent with the results of this study [39,40]. Restricted by the terrain and the coverage of the base station, RTK will produce large errors when applied to the coordinate acquisition of GCP under a forest that is heavily shaded by trees. Hence, based on our research results, in forest scenes, we recommend combining the nadir photo and the 30° tilted photo of the drone to obtain the canopy 3D point cloud. In the SFM process, homology point-matching between images is the key step that directly determines the canopy 3D point-cloud density and the utilization of image information [41]. The result of homologous point matching between photos of different viewing angles changes and becomes unstable. Thus, a large number of oblique photos were used in this research to ensure the integrity of the canopy structure, which led to the point cloud derived from the multi-angle scheme being more refined and not having a large range of noise compared with the single-oblique angle schemes. In addition, lens resolution and flight height of UAV are also important factors that determine the density of derived point clouds [42]. Even if the high-resolution photos of drones flying at low altitudes are very effective in capturing canopy details [43], low content of photo content from a low flight altitude requires a high image overlap rate to get an excellent image matching result, while numerous photo sets generated in this way will greatly increase the cost of photogrammetry and computer hardware requirements. Thus, adjusting the UAV flight height according to lens resolution is very important for the reconstruction accuracy of the 3D model.

Leaf Area Estimation
Numerous remote sensing estimation methods on forest LAI have been reported in the past [44,45], while the instability of multispectral data for model fitting and the high In the SFM process, homology point-matching between images is the key step that directly determines the canopy 3D point-cloud density and the utilization of image information [41]. The result of homologous point matching between photos of different viewing angles changes and becomes unstable. Thus, a large number of oblique photos were used in this research to ensure the integrity of the canopy structure, which led to the point cloud derived from the multi-angle scheme being more refined and not having a large range of noise compared with the single-oblique angle schemes. In addition, lens resolution and flight height of UAV are also important factors that determine the density of derived point clouds [42]. Even if the high-resolution photos of drones flying at low altitudes are very effective in capturing canopy details [43], low content of photo content from a low flight altitude requires a high image overlap rate to get an excellent image matching result, while numerous photo sets generated in this way will greatly increase the cost of photogrammetry and computer hardware requirements. Thus, adjusting the UAV flight height according to lens resolution is very important for the reconstruction accuracy of the 3D model.

Leaf Area Estimation
Numerous remote sensing estimation methods on forest LAI have been reported in the past [44,45], while the instability of multispectral data for model fitting and the high cost of LiDAR data have restricted forest LAI estimation. In this study, based on the SFM-MVS workflow, a 3D point cloud of the forest canopy was generated from a collection of photos from different angles of the drone, and the LAI was calculated based on the voxelized model. Our results show that the 3D point cloud derived from photogrammetry using nadir photos and 30 • oblique photos provides the best LAI estimation. Although Lisein et al. [46] suggested that the canopy point cloud generated based on photogrammetry is incomplete because the optical image is not penetrating, they did not discuss the complement ability of oblique photography for canopy structure. Combining the LAD incremental comparison results of different oblique photography schemes relative to the nadir-photo-only scheme, we believe that coupling oblique photography and nadiral photography can compensate to a certain extent for the weak acquisition ability of the nadir-photo-only scheme for the canopy structure. Moreover, the UAV that executes the orthogonal flight plan can obtain canopy information from four directions, which further improves the canopy information obtained by oblique photography. Although detecting the complete canopy internal information is still difficult, oblique photography can obtain a very fine leaf point cloud for tree species with exposed leaves, which could avoid the overestimation of LAI caused by the existence of branches inside the canopy. The canopy structure of Masson pine in our study area is simple, and the leaves are relatively sparse due to poor soil [47], so the leaf information can be obtained very clearly and completely from the oblique image ( Figure 12). metry is incomplete because the optical image is not penetrating, they did not discuss the complement ability of oblique photography for canopy structure. Combining the LAD incremental comparison results of different oblique photography schemes relative to the nadir-photo-only scheme, we believe that coupling oblique photography and nadiral photography can compensate to a certain extent for the weak acquisition ability of the nadirphoto-only scheme for the canopy structure. Moreover, the UAV that executes the orthogonal flight plan can obtain canopy information from four directions, which further improves the canopy information obtained by oblique photography. Although detecting the complete canopy internal information is still difficult, oblique photography can obtain a very fine leaf point cloud for tree species with exposed leaves, which could avoid the overestimation of LAI caused by the existence of branches inside the canopy. The canopy structure of Masson pine in our study area is simple, and the leaves are relatively sparse due to poor soil [47], so the leaf information can be obtained very clearly and completely from the oblique image ( Figure 12). Sub-voxels of different sizes have different abilities to fill gaps in the point clouds [48]. The LAI estimation derived from the voxelization model increase with the increase of the sub-voxel size [49]. As long as the relatively complete canopy structure is retained, even point clouds with low density can also improve leaf area estimation by increasing the sub-voxel size (O, T15, T30). Therefore, the determination of the optimal sub-voxel Sub-voxels of different sizes have different abilities to fill gaps in the point clouds [48]. The LAI estimation derived from the voxelization model increase with the increase of the sub-voxel size [49]. As long as the relatively complete canopy structure is retained, even point clouds with low density can also improve leaf area estimation by increasing the sub-voxel size (O, T15, T30). Therefore, the determination of the optimal sub-voxel size is extremely critical for forest leaf area estimation. However, the density of the point cloud derived from the nadir-photo only scheme (O) has a big difference in vertical space [50]. The point cloud derived from O is very fine due to excellent image details at the top of canopy but is sparse due to the limitation of nadir photo viewing angle in the middle and lower parts of the canopy. Thus, the good LAI estimation given by O is based on the fact that the voxels in the lower part of the canopy contain more blank areas in exchange, which determines that the nadir-photo-only scheme is not suitable for LAI estimation. The lower and middle structure of the two single-oblique angle schemes (T15, T30) is more complete than that of O. Nevertheless, the point cloud derived from oblique-only photos schemes is sparse and has many noise points because of the poor homology points matching results, which fail to reflect the details of canopy well. Thus, even if the sub-voxel size is increased to accommodate more blanks, the oblique-only photo schemes alone still can not give an accurate LAI estimate.
Multi-angle schemes combine the excellent image details of nadir photos for the top part of the canopy and the strong acquisition ability of oblique photos for the middle and lower canopy, which led to good performance in LAI estimation. The top and middle parts of the canopy are less affected by the canopy occlusion, so the multi-angle schemes in different LAI areas always provide a fairly accurate leaf area estimation. For the lower part of the canopy, the point cloud details become more complete as the angle of the oblique photo increases in the low-LAI area, which is not very significant in high-LAI areas because of the overlapping occlusions between canopies. In low-LAI areas, overlapping occlusions are rare because the lower part of the canopy is largely exposed due to the scarcity of trees, so the addition of oblique photos can significantly improve the LAD estimation of the middle and lower canopies. In the high-LAI area, the canopy overlap between the single wood blocks the lower canopy layer, so the UAV cannot get fine point clouds in this part of the canopy, which also leads to the underestimation of the leaf area under the canopy in the high LAI area.

Limitation of this Study
Since there are no real leaf area data for the lower part of the canopy, we cannot evaluate the LAD estimates for different heights of the canopy generated by the multi-angle scheme. Even though the LAD increment in this study is relative to the nadir-photo-only scheme, it can also reflect the improvement in leaf area estimation by point clouds derived from multi-angle schemes to a certain extent. In this study, we did not set up different altitude UAV flight plans to test their ability for forest LAI estimation, and the impact of increasing the resolution of the UAV camera on the estimation of LAI using point clouds generated from multi-angle images is also still unknown. Even if they are not the focus of this study, they may provide the potential for insufficient leaf area estimation in the lower part of the canopy in areas with canopy occlusion. Finally, obtaining sub-centimeter-level three-dimensional structures of the forest surface without using a large amount of GCP is one of the most important issues for future SFM-MVS applications in forest surveys.

Conclusions
Under the background that the SFM-MVS workflow technology is very mature, we explored the feasibility of using point cloud data derived from photogrammetry to estimate forest LAI. The results show that the scheme with oblique photos (OT15, OT30, T15, T30) always provides better checkpoint accuracy than schemes with nadir-photos only (O). For the four oblique photography schemes, the surface model accuracy of the multi-angle scheme (OT15 and OT30) is better than that of the single oblique angle scheme (T15 and T30). The point clouds data derived from the coupled nadir photos and 30 • oblique photos provide the best checkpoint accuracy (RMSE (H) = 0.2917 m, RMSE (V) = 0.1797 m) and the best LAI estimation by sub-voxel size of 0.09 m (RMSE = 0.1790 m 2 /m 2 ), and its improvement in the estimation of leaf area in the middle and lower part of the canopy is better than the single angle scheme (O, T15, T30) and the multi-angle solution with 15 • oblique photos (OT15). Moreover, the LAI estimation provided by the scheme with 30 • photos is always better than that of the scheme with 15 • photos (T30 better than T15, OT30 better than OT15).

Data Availability Statement:
The source codes developed in this study were donated to GitHub (https:// github.com/TOTOROLLC/Forest-Stand-LAI-Remote-Sensing-Retrieval-Based-on-Photogrammetry (accessed on 7 January 2021)). And the data used to support the findings of this study are available from the corresponding author upon request.