Effect of Leaf Occlusion on Leaf Area Index Inversion of Maize Using UAV–LiDAR Data

The leaf area index (LAI) is a key parameter for describing crop canopy structure, and is of great importance for early nutrition diagnosis and breeding research. Light detection and ranging (LiDAR) is an active remote sensing technology that can detect the vertical distribution of a crop canopy. To quantitatively analyze the influence of the occlusion effect, three flights of multi-route high-density LiDAR dataset were acquired at two time points, using an Unmanned Aerial Vehicle (UAV)-mounted RIEGL VUX-1 laser scanner at an altitude of 15 m, to evaluate the validity of LAI estimation, in different layers, under different planting densities. The result revealed that normalized root-mean-square error (NRMSE) for the upper, middle, and lower layers were 10.8%, 12.4%, 42.8%, for 27,495 plants/ha, respectively. The relationship between the route direction and ridge direction was compared, and found that the direction of flight perpendicular to the maize planting ridge was better than that parallel to the maize planting ridge. The voxel-based method was used to invert the LAI, and we concluded that the optimal voxel size were concentrated on 0.040 m to 0.055 m, which was approximately 1.7 to 2.3 times of the average ground point distance. The detection of the occlusion effect in different layers under different planting densities, the relationship between the route and ridge directions, and the optimal voxel size could provide a guideline for UAV–LiDAR application in the crop canopy structure analysis.


Introduction
The leaf area index (LAI) is an important crop canopy structure parameter usually defined as half of the total leaf surface area of a unit of surface area [1,2].The LAI can be used to reflect the growth status of crops, which is affected by factors such as size, age, and plant spacing.LAI is closely related to processes, such as respiration, transpiration, precipitation interception of crops, and the carbon cycle [3].The study of the vertical structure of maize is of great significance to the analysis of maize photosynthesis, pollen transmission, and stress resistance.Quantitative analysis of vertical distribution and dynamic changes of LAI for crops is of great importance for early nutrition diagnosis and breeding research.Therefore, it is necessary to accurately estimate the LAI of crops.
Currently, LAI measurement uses mainly direct and indirect methods [4,5].Direct measurement methods include destructive sampling, allometric growth equations, and the litter collection method [6].The litter collection method is more widely used because it is not destructive.Direct measurement methods are labor-intensive, time-consuming, and costly, and obtains data at only a limited sampling point, which is impractical in a real-life situation.Indirect measurement methods include a parameter inversion method and a remote sensing optical measurement method.The parameter inversion method uses some measurement parameters, such as height, to obtain the LAI [7].This method not only avoids the damage caused by the direct method, but are also convenient and fast.However, the measurement accuracy is lower than that of the direct measurement methods.The optical measurement of remote sensing uses optical instruments, such as a Tracing Radiation and Architecture of Canopies (TRAC), an LAI-2000 plant canopy analyzer (LI-COR, USA), and hemispherical photography to obtain the LAI [8].However, it is difficult for optical instruments to penetrate the maize canopy and to provide precise vertical structural information.Its spectral response is affected by only the upper canopy of the vegetation, and it cannot express an LAI that increases due to the inner leaves of the canopy [9].
Light Detection and Ranging (LiDAR) is an active remote sensing technology that provides detailed three dimensional information of canopy structure, which is particularly useful for studying the vertical distribution of LAI [10][11][12][13][14]. Unmanned Aerial Vehicle (UAV)-LiDAR is mounted on UAV to obtain a wider range of three-dimensional information.Currently, LiDAR has been applied to some crop parameter extractions [15,16], and many experts and scholars have discussed it.According to the difference of the format, LiDAR data can be divided into two types-point cloud and waveform [17,18], and in this study, the point cloud format was used.
For the estimation of LAI using LiDAR data, the gap fraction method and the voxel-based method are commonly used [19].Chen et al. [19] presented a method to reduce the clumping effect and estimate LAI using Terrestrial Laser Scanning (TLS) data, based on the gap fraction theory.Zhu et al. [20] obtained clumping index by using the gap size distribution method, while the woody contribution was evaluated, based on an improved point classification between woody and foliar materials.However, the spatial heterogeneity caused by various path lengths within canopies has not been considered using a simple clumping index.It was found to be another important and non-negligible cause for LAI underestimation [21].The voxel-based method is based on Wilson's contact frequency theory, and the voxel size is sensitive for LAI inversion.Su et al. [22] separated stems and leaves of maize, using the Difference of Normal algorithm (DON), and inverted the LAI using TLS, based on Wilson's contact frequency theory, and the selection of voxel size was related to the width of maize leaves.Zhang et al. [23] used voxel-based method to extract height-related and density-related metrics over Unmanned Aerial Vehicle (UAV)-LiDAR, to predict forest structure parameters, and the factors affecting the optimal voxel size included plot size, LiDAR, and forest structure characteristics.Many other unknown factors affect the optimal voxel size; these need to be considered in the determination of the optimal voxel size.
The occlusion effect occurring in UAV-LiDAR data is one of the main limitations of the acquisition of crop parameters.Many scholars have discussed the influences of the occlusion effect and the procedures for reducing its influence.Holmgren et al. [24] concluded that lower height percentiles were affected more than upper height percentiles, due to obscuration and the longer path of laser pulse penetrating into the canopy.Bauwens et al. [25] assessed and compared a hand-held mobile laser scanner (HMLS), with two TLS approaches (single scan (SS) and multi scan (MS)) for the estimation of several forest parameters in a wide range of forest types and structures, and found that MS reduced the occlusion effect.Hamraz et al. [26] modeled the occlusion effect in terms of point density, estimated the fractions of points representing the different canopy layers (one overstory and multiple understory) and also pinpointed the required density for reasonable tree segmentation (where the accuracy plateaus).However, none of these studies mentioned a quantitative analysis of the occlusion effect.Su et al. [27] used a terrestrial laser scanner to calculate plant height, plant area index, and projected leaf area at the individual plant level, and each individual plant was divided into five height strata to validate the occlusion effect.Qin et al. [28] divided the maize into five 0.5 m high layers, they used LiDAR waveform data to invert the Fraction of Photosynthetically Active Radiation (FPAR) and studied the vertical distribution of the FPAR of the maize.The above studies analyzed the occlusion effect as a way of layering, and the occlusion effect at different planting densities was not considered.
Incidence angle refers to the angle between the spot direction and the vertical direction of the laser beam emitted by the LiDAR.Most studies estimated LAI from spaceborne or airborne LiDAR data, with a large distance from the ground, so the incidence angle was often overlooked or a nadir view was assumed [10,17,29].However, with the development of UAV technology, the flight altitude of UAV got increasingly lower, due to its short action distance and large field of view, the optimal incidence angle of UAV-LiDAR could better reflect the structure of maize [30].Ridge direction is the planting direction of maize.The occlusion effect of route direction, parallel to the planting ridge direction and perpendicular to the planting ridge direction might be different, and few people have studied it.
The specific objectives of this study were-(1) different layers and different planting densities were used to quantitatively analyze the occlusion effect and further analyze the effectiveness of LAI estimation; (2) relationship between route direction and ridge direction was compared to acquire a better LiDAR data; (3) optimal voxel size was determined and the relationship between optimal voxel size and point cloud density was analyzed; and (4) factors affecting the selection of optimal voxel size, and the effect of different incidence angle size and the size of incidence angle range were discussed.

Experimental Design
This study was carried out at the National Experimental Station for Precision Agriculture, which is in the Changping district of Beijing (40 • 10.6 N, 116 • 26.3 E).This site has a mean annual rainfall of 508 mm and a mean annual temperature of 13 • C [31].The Jinghua 38 maize cultivar was used.The experiment was carried out at the Silking and Blister stages (R1 and R2) [32], the maize at two time points, all belonged to the reproductive growth stage, and there was no significant difference in plant height, no leaf was added, and only the bottom two leaves showed drying and yellowing.Figure 1a1,b2 show the route design diagrams for the silking and blister stages, respectively, including the distance between the routes, and the location relations between the experimental area and each route; the route distance in a1 and b2 were 25 m and 15 m, respectively.The size of the experimental area was 15 m × 15 m.As shown in Figure 1a2,b1, 30 plots-each 2.5 m × 3 m-were investigated.In the north-south direction, there were five different maize planting densities, and in the east-west direction there were six replicated experiments of the same density.Three layers of de-leaf treatment was conducted in the experiment.According to the agronomic function of maize leaves, all the leaves of each maize in the experimental area were divided into upper, middle, and lower layers, as shown in Figure 1c.The lower layer included all leaves that ranged from the second leaf below the ear to the root, the middle layer was four leaves above the lower layer, and the top layer was the rest of the leaves.Three UAV flights with LiDAR sensor were carried out to acquire the point cloud data with the same sensor configuration and route planning.The first flight obtained all plant point cloud data, the second obtained the data after removing the lower layer manually, and the third obtained the data after removing the middle layer manually.Figure 2 shows the three flights of point cloud data obtained in this study.In the silking and blister stages, each flight had eight or seven routes and each route corresponded to a different incidence angle.Figure 1d is a schematic diagram of the incidence angles of the R 11 , R 12 , R 13 , and R 14 routes and the settings of different densities in the experimental area.

LiDAR Data
The UAV-LiDAR data were acquired on 28 August and 14 September 2018, using the Riegl VUX-1 (RIEGL Co., Austria) sensor, flown at 15 m above ground level, with a flight speed of 3 m•s -1 .The scanning frequency was 550 kHz, the UAV-LiDAR obtained the data using pendulum scanning; the laser foot point distribution of that scanning method is shown in Figure 3, the spot diameter was 0.0075 m, with an average ground point distance of approximately 0.0239 m, the format of LiDAR

LiDAR Data
The UAV-LiDAR data were acquired on 28 August and 14 September 2018, using the Riegl VUX-1 (RIEGL Co., Austria) sensor, flown at 15 m above ground level, with a flight speed of 3 m•s -1 .The scanning frequency was 550 kHz, the UAV-LiDAR obtained the data using pendulum scanning; the laser foot point distribution of that scanning method is shown in Figure 3, the spot diameter was 0.0075 m, with an average ground point distance of approximately 0.0239 m, the format of LiDAR

LiDAR Data
The UAV-LiDAR data were acquired on 28 August and 14 September 2018, using the Riegl VUX-1 (RIEGL Co., Austria) sensor, flown at 15 m above ground level, with a flight speed of 3 m•s −1 .The scanning frequency was 550 kHz, the UAV-LiDAR obtained the data using pendulum scanning; the laser foot point distribution of that scanning method is shown in Figure 3, the spot diameter was 0.0075 m, with an average ground point distance of approximately 0.0239 m, the format of LiDAR data was LAS, and the flight parameters of the two time points were same.The incidence angle range of the different routes obtained and the point cloud density of each route [point cloud number in the experimental area / (length × width)] are shown in Table 1, and we can see from the table, the incidence angles of the east-west direction were similar to the north-south direction, such as R 11 (69 • to 77   1, and we can see from the table, the incidence angles of the east-west direction were similar to the north-south direction, such as R11 (69° to 77°) and R18 (68°-76°), R12 (−64° to −32°) and R17 (-60° to −25°), R13 (−34° to 29°) and R16 (−54° to 0°), R14 (60° to 74°) and R15 (63° to 73°), R22 (−12° to −5°) and R26 (3° to 18°), and R23 (8° to 32°) and R25 (−30° to 0°), where "to" refers to the UAV flight direction; from west to east and north to south was indicated as positive, and the opposite was indicated as negative.

Field Data
To provide a reference value for comparison with the point cloud data, the leaf area (LA) was manually measured.Six maize plants were randomly selected from each plot for marking, and the lower, middle, and upper layers of the six maize plants in each plot were collected, in order.The length and width of the leaves were measured by a ruler indoors, the leaf area of each leaf was

Field Data
To provide a reference value for comparison with the point cloud data, the leaf area (LA) was manually measured.Six maize plants were randomly selected from each plot for marking, and the lower, middle, and upper layers of the six maize plants in each plot were collected, in order.The length and width of the leaves were measured by a ruler indoors, the leaf area of each leaf was calculated using Equation (1) [33].Six ground control points were distributed evenly on both sides of the experimental area, and real-time kinematic GPS was used to measure the control points.
where L is the length, W is the width of the widest part of the leaf, and a is a correction coefficient.Many measurement results have shown that a is 0.75 for maize [34].

Data Preprocessing
To obtain a real maize point cloud, the data acquired by a laser LiDAR must be preprocessed.Preprocessing includes two parts-a trajectory solution comprising the raw laser data processing and point cloud data preprocessing, and a point cloud data layering.The trajectory solution, the raw laser data processing, and the point cloud data preprocessing are shown in Figure 4. First, a difference method and a GPS level arm inverse calculation were used to calculate the accurate trajectory data for the GPS base station data and the Position and Orientation System (POS) raw data.Then, the RiPROCESS software (version 1.7.2,RIEGL Co., Austria) was used to process the raw laser data, including a waveform solution, track data matching, and three-dimensional point cloud data visualization.Finally, point cloud data in the LAS file format was exported.Those data were de-noised, registered, filtered to remove the ground points, and segmented using the LiDAR360 Suite software (version 3.0, GreenValley Co., Beijing, China) to complete the preprocessing.Finally, the point cloud data were stratified.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 16 calculated using Equation (1) [33].Six ground control points were distributed evenly on both sides of the experimental area, and real-time kinematic GPS was used to measure the control points.
where L is the length, W is the width of the widest part of the leaf, and a is a correction coefficient.Many measurement results have shown that a is 0.75 for maize [34].

Data Preprocessing
To obtain a real maize point cloud, the data acquired by a laser LiDAR must be preprocessed.Preprocessing includes two parts-a trajectory solution comprising the raw laser data processing and point cloud data preprocessing, and a point cloud data layering.The trajectory solution, the raw laser data processing, and the point cloud data preprocessing are shown in Figure 4. First, a difference method and a GPS level arm inverse calculation were used to calculate the accurate trajectory data for the GPS base station data and the Position and Orientation System (POS) raw data.Then, the RiPROCESS software (version 1.7.2,RIEGL Co., Austria) was used to process the raw laser data, including a waveform solution, track data matching, and three-dimensional point cloud data visualization.Finally, point cloud data in the LAS file format was exported.Those data were de-noised, registered, filtered to remove the ground points, and segmented using the LiDAR360 Suite software (version 3.0, GreenValley Co., Beijing, China) to complete the preprocessing.Finally, the point cloud data were stratified.As shown in Figure 5, an agronomic stratification method was adopted to process the data of the three flights, to obtain the point cloud data of the upper, middle, and lower layers of maize.First, the stalk data were removed from the third flights to obtain the upper layer point cloud data.Then, the upper layer and the stalk data were removed from the second flights, to obtain the middle layer point cloud data.Finally, after the upper and middle layers were removed from the data of the first flights, the point cloud data by the lower layer was obtained, and the layering processing of the point cloud data was completed in this step.As shown in Figure 5, an agronomic stratification method was adopted to process the data of the three flights, to obtain the point cloud data of the upper, middle, and lower layers of maize.First, the stalk data were removed from the third flights to obtain the upper layer point cloud data.Then, the upper layer and the stalk data were removed from the second flights, to obtain the middle layer point cloud data.Finally, after the upper and middle layers were removed from the data of the first flights, the point cloud data by the lower layer was obtained, and the layering processing of the point cloud data was completed in this step.

Estimation of LAI
As a method of LiDAR point cloud data visualization, the voxel-based method has been widely used in LiDAR point cloud data processing [16,33,35].The voxel size plays an important role in the estimation of maize structure parameters.In this study, the voxel-based method was adopted for LAI extraction of point cloud data.By determining the boundaries of the point cloud in the experimental area, the point cloud was divided into equal-sized grids, as shown in Figure 6, and the point cloud was judged to exist in the grid.As more factors affect the size of the voxel, the voxel-based method is more sensitive, so the optimal voxel size is different in different situations.In this study, the LAI of each route was inverted with a voxel size of 0.01 m to 0.1 m, at intervals of 0.005 m, and then the optimal voxel size of each route was determined.LAI was given by Equation (2).
where  () is the number of grids with point clouds in the kth layer and  () is the total number of grids in the kth layer.

Estimation of LAI
As a method of LiDAR point cloud data visualization, the voxel-based method has been widely used in LiDAR point cloud data processing [16,33,35].The voxel size plays an important role in the estimation of maize structure parameters.In this study, the voxel-based method was adopted for LAI extraction of point cloud data.By determining the boundaries of the point cloud in the experimental area, the point cloud was divided into equal-sized grids, as shown in Figure 6, and the point cloud was judged to exist in the grid.As more factors affect the size of the voxel, the voxel-based method is more sensitive, so the optimal voxel size is different in different situations.In this study, the LAI of each route was inverted with a voxel size of 0.01 m to 0.1 m, at intervals of 0.005 m, and then the optimal voxel size of each route was determined.LAI was given by Equation (2).
where n I (k) is the number of grids with point clouds in the kth layer and n T (k) is the total number of grids in the kth layer.

Estimation of LAI
As a method of LiDAR point cloud data visualization, the voxel-based method has been widely used in LiDAR point cloud data processing [16,33,35].The voxel size plays an important role in the estimation of maize structure parameters.In this study, the voxel-based method was adopted for LAI extraction of point cloud data.By determining the boundaries of the point cloud in the experimental area, the point cloud was divided into equal-sized grids, as shown in Figure 6, and the point cloud was judged to exist in the grid.As more factors affect the size of the voxel, the voxel-based method is more sensitive, so the optimal voxel size is different in different situations.In this study, the LAI of each route was inverted with a voxel size of 0.01 m to 0.1 m, at intervals of 0.005 m, and then the optimal voxel size of each route was determined.LAI was given by Equation (2).
where  () is the number of grids with point clouds in the kth layer and  () is the total number of grids in the kth layer.

Analysis and Validation
In order to quantitatively analyze the occlusion effect, the following points were analyzed in this study.(1) As the upper layer was less affected by the occlusion effect, we used the upper layer to obtain the optimal voxel size for each route, and through comparison and analysis, the relationship between the optimal voxel size and the average ground point distance was obtained.(2) To determine the influence of the ridge direction on the accuracy of the point cloud data acquisition of UAV-LiDAR, the whole plant maize was considered.Therefore, the route direction of the UAV, parallel to the ridge direction and perpendicular to the ridge direction of maize planting was compared.(3) R 12 and R 22 with the highest inversion accuracy in the two time points were used to quantitatively analyze the occlusion effect on different densities and layers.
In this study, we compared the LAI of the UAV-LiDAR point cloud data inversion with the measured LAI.By evaluating the accuracy of the R 2 and the normalized root-mean-square error (NRMSE), the optimal incidence angle, effective detecting depth, and point cloud data inversion of each layer were determined.The coefficient of determination (R 2 ), the root-mean-square error (RMSE), and the normalized root-mean-square error (NRMSE) were calculated by Equations (3)-( 4).
where x i is the LAI measured value of the plot i; x is the average of the measured values of LAI in each plot; f i is the i-plot LAI value obtained by the point cloud; n is the number of plots.In this study, n is 15; i is the plot number.

Determination of the Optimal Voxel Size
The measured LAI and the LAI of different densities are shown in Table 2, and compared with the upper and lower layers.We can see from the table that the LAI of the lower layer was the least concentrated.Considering the occlusion effect of UAV-LiDAR, we compared the LAI inverted from the upper layer to the LAI measured in the field of 15 routes, R 2 and NRMSE were obtained to determine the optimal voxel size.As shown in Figure 7, the x coordinate was the routes, the y coordinate was the voxel sizes, and the range of voxel size was 0.01 m to 0.1 m.As can be seen from Figure 7, in each route, with the increase of voxel size, the R 2 showed a little increase, and the NRMSE showed a decreasing and then increasing trend.Therefore, the voxel size corresponding to the minimum value of the NRMSE of each route was the optimal voxel size.3 shows the optimal voxel size for each route and the corresponding R 2 and NRMSE.Figure 8 shows the relationship between the point cloud density and the optimal voxel size.As can be seen from Figure 8, the optimal voxel size became smaller with an increase in point cloud Table 3 shows the optimal voxel size for each route and the corresponding R 2 and NRMSE.Figure 8 shows the relationship between the point cloud density and the optimal voxel size.As can be seen from Figure 8, the optimal voxel size became smaller with an increase in point cloud densities, and the optimal voxel sizes were concentrated on 0.040 m to 0.055 m, considering the average ground points distance of approximately 0.0239 m.Therefore, the optimal voxel size is approximately 1.7 to 2.3 times of the average ground point distance.Among them, the highest R 2 was 0.94 and the lowest was 0.70, and the highest NRMSE was 1.9% and the lowest was 39.6%.densities, and the optimal voxel sizes were concentrated on 0.040 m to 0.055 m, considering the average ground points distance of approximately 0.0239 m.Therefore, the optimal voxel size is approximately 1.7 to 2.3 times of the average ground point distance.Among them, the highest R 2 was 0.94 and the lowest was 0.70, and the highest NRMSE was 1.9% and the lowest was 39.6%.

Relationship between Route Direction and Ridge Direction
Figure 9 shows the accuracy comparison and analysis of the measured LAI and the point cloud data inversion LAI, of the whole maize plant in each plot.Based on the optimal voxel size of each route in 3.1, the inversion accuracy of the whole plant was obtained.As shown in Table 4, the similar incidence angles of east-west direction and north-south direction were selected, and the LAI inversion accuracy was compared, R 2 and NRMSE of the six east-west routes were better than the six north-south routes, respectively.In the north-south direction, the LAI inversion accuracy of R22 was the highest (NRMSE and R 2 were 4.6% and 0.96, respectively), in the east-west direction, the LAI inversion accuracy of R26 was the highest (NRMSE and R 2 were 17.7% and 0.84, respectively).
Considering the relation between the route direction and the maize planting ridge direction, in the experimental area, we concluded that when the UAV-LiDAR was collected data, the direction of flight perpendicular to the maize planting ridge was better than that parallel to the maize planting ridge.

Relationship between Route Direction and Ridge Direction
Figure 9 shows the accuracy comparison and analysis of the measured LAI and the point cloud data inversion LAI, of the whole maize plant in each plot.Based on the optimal voxel size of each route in 3.1, the inversion accuracy of the whole plant was obtained.As shown in Table 4, the similar incidence angles of east-west direction and north-south direction were selected, and the LAI inversion accuracy was compared, R 2 and NRMSE of the six east-west routes were better than the six north-south routes, respectively.In the north-south direction, the LAI inversion accuracy of R 22 was the highest (NRMSE and R 2 were 4.6% and 0.96, respectively), in the east-west direction, the LAI inversion accuracy of R 26 was the highest (NRMSE and R 2 were 17.7% and 0.84, respectively).densities, and the optimal voxel sizes were concentrated on 0.040 m to 0.055 m, considering the average ground points distance of approximately 0.0239 m.Therefore, the optimal voxel size is approximately 1.7 to 2.3 times of the average ground point distance.Among them, the highest R 2 was 0.94 and the lowest was 0.70, and the highest NRMSE was 1.9% and the lowest was 39.6%.

Relationship between Route Direction and Ridge Direction
Figure 9 shows the accuracy comparison and analysis of the measured LAI and the point cloud data inversion LAI, of the whole maize plant in each plot.Based on the optimal voxel size of each route in 3.1, the inversion accuracy of the whole plant was obtained.As shown in Table 4, the similar incidence angles of east-west direction and north-south direction were selected, and the LAI inversion accuracy was compared, R 2 and NRMSE of the six east-west routes were better than the six north-south routes, respectively.In the north-south direction, the LAI inversion accuracy of R22 was the highest (NRMSE and R 2 were 4.6% and 0.96, respectively), in the east-west direction, the LAI inversion accuracy of R26 was the highest (NRMSE and R 2 were 17.7% and 0.84, respectively).
Considering the relation between the route direction and the maize planting ridge direction, in the experimental area, we concluded that when the UAV-LiDAR was collected data, the direction of flight perpendicular to the maize planting ridge was better than that parallel to the maize planting ridge.Considering the relation between the route direction and the maize planting ridge direction, in the experimental area, we concluded that when the UAV-LiDAR was collected data, the direction of flight perpendicular to the maize planting ridge was better than that parallel to the maize planting ridge.

Quantitative Analysis of the Occlusion Effect
As stated in Section 3.1, the LAI inversion accuracy of the R 12 and R 22 was the highest in the two time points, respectively, so the R 12 and R 22 routes were selected to analyze the influence of the occlusion effect.Figure 10 shows the LAI inversion accuracy comparison for the different planting densities and layers.As shown in Figure 10, when the planting densities were D1, D2, D3, and D4, as the densities decreased, the NRMSE showed a downward trend.When the planting density was D5, an abnormal trend occurred, because there were few maize plants in the plots, and the sample size collected was not large enough.In addition, the inverted accuracy of the upper layer was the highest (and good) in different planting densities (21.4%, 17.7%, 15.7%, 10.8%, and 49.5%), less in the middle layer (45.1%, 34.8%, 32.9%, 12.4%, and 53.8%), and worst in the lower layer (66.7%, 47.7%, 38.3%, 42.8%, and 63.3%).When the planting density was 27,495 plants/ha, the inversion accuracy of LAI was the highest (the upper, middle, and lower layers were 10.8%, 12.4%, and 42.8%).NRMSE of the lower layer was an abnormal value, this might be due to the effect of the wheat stubble on the ground.When the planting density was 88,050 plants/ha, the inversion accuracy of LAI was the lowest (the upper, middle, and lower layers were 21.4%, 45.1%, and 66.7%).

Quantitative Analysis of the Occlusion Effect
As stated in Section 3.1, the LAI inversion accuracy of the R12 and R22 was the highest in the two time points, respectively, so the R12 and R22 routes were selected to analyze the influence of the occlusion effect.Figure 10 shows the LAI inversion accuracy comparison for the different planting densities and layers.As shown in Figure 10, when the planting densities were D1, D2, D3, and D4, as the densities decreased, the NRMSE showed a downward trend.When the planting density was D5, an abnormal trend occurred, because there were few maize plants in the plots, and the sample size collected was not large enough.In addition, the inverted accuracy of the upper layer was the highest (and good) in different planting densities (21.4%, 17.7%, 15.7%, 10.8%, and 49.5%), less in the middle layer (45.1%, 34.8%, 32.9%, 12.4%, and 53.8%), and worst in the lower layer (66.7%, 47.7%, 38.3%, 42.8%, and 63.3%).When the planting density was 27,495 plants/ha, the inversion accuracy of LAI was the highest (the upper, middle, and lower layers were 10.8%, 12.4%, and 42.8%).NRMSE of the lower layer was an abnormal value, this might be due to the effect of the wheat stubble on the ground.When the planting density was 88,050 plants/ha, the inversion accuracy of LAI was the lowest (the upper, middle, and lower layers were 21.4%, 45.1%, and 66.7%).

Analysis of the Optimal Incidence Angle
Many studies have discussed the incidence angles.Korhonen et al. [36] have discussed the effect of incidence angle range on the accuracy of the estimated Gap fraction (Pgap) from airborne

Analysis of the Optimal Incidence Angle
Many studies have discussed the incidence angles.Korhonen et al. [36] have discussed the effect of incidence angle range on the accuracy of the estimated Gap fraction (P gap ) from airborne LiDAR, including a small range of scan angles (0 • -15 • ) and a large range of scan angles (0 • -75 • ), which proved the relationship between the accuracy of the estimated P gap and the incidence angle range.However, this comparison might not be appropriate because the incidence angle of 0 • -75 • was too large.Liu et al. [37] acquired airborne LiDAR data from nadir (0 • -7 • ), small off-nadir (7 • -23 • ), and large off-nadir (23 • -38 • ) directions, to calculate both P gap and vertical P gap profile, and Digital Hemispherical Photographs (DHP) were used as references, for validation.The results showed that angular P gap from airborne LiDAR correlated well with the angular P gap from DHP (R 2 = 0.74, 0.87, and 0.67 for nadir, small off-nadir, and large off-nadir directions).From the results of this research, the scan angle was found to have a large impact on the measured P gap from airborne LiDAR.Cao et al. [11] divided all of the correctly detected trees into two categories, to assess the effect of scan angle on species discrimination, one group included tree crowns with small scan angles (i.e., the mean of the absolute scan angles ≤15 • ), whereas the other group included tree crowns with large scan angles (i.e., the mean of the absolute scan angles between 15 • and 30 • ).Tree crowns with small scan angles had the highest classification accuracies (overall accuracy = 69.4%,Kappa coefficient = 0.633), which was slightly higher than trees crowns with large scan angles (overall accuracy = 67.3%,Kappa coefficient = 0.608).Compared with the two studies, the incidence angle of the study was wider and covered more cases.Table 5 shows the comparison of the different incidence angles.The incidence angles in angle 1 were all greater than 60 • , and the inversion accuracy was the lowest; angle 2 and angle 4 were orthographic data, angle 3 and angle 5 were non-orthogonal data, and compared with angle 4 and angle 5, angle 2 and angle 3 had a wider range of incidence angles.As can be seen from Table 5, the inversion accuracy of angle 4 and angle 5 were all better than angle 2 and angle 3 (NRMSE and R 2 were 11.7% and 0.79 for angle 4, NRMSE and R 2 were 4.6% and 0.96 for angle 5, NRMSE and R 2 were 43.8% and 0.77 for angle 2, NRMSE and R 2 were 42.5% and 0.93 for angle 3), which showed a smaller range of incidence angle and a higher accuracy of LAI inversion.In angle 4 and angle 5, comparing R 21 , R 27 with R 26 , R 22 , the inversion accuracy of R 26 , R 22 were found to be higher than R 21 , R 27 (NRMSE and R 2 were 17.7% and 0.84, 4.6% and 0.96, 31.6% and 0.69, 11.7% and 0.79, for R 26 , R 22 , R 21 , and R 27 , respectively), which proved that the non-orthogonal data had a higher inversion accuracy.To sum up, the incidence angle was too large or orthonormal; the LAI inversion accuracy was the highest when the incident angle was small-smaller the range of the incidence angle, higher the accuracy of LAI inversion.); angle 2 (large incident angle range for the orthographic data); angle 3 (large incident angle range for the non-orthogonal data); angle 4 (small incident angle range for the orthographic data); and angle 5 (small incident angle range for the non-orthogonal data).

Factors Influencing the Selection of the Optimal Voxel Size
The voxel-based method representing the canopy elements were abstracted by a volume grid and placed in a 3D grid [23].Therefore, the selection of voxel size was very important [38,39].If voxels were too small, it might produce redundant unfilled voxels of empty volume, containing few canopy elements, which might lead to an underestimation of the canopy structural parameters.However, too large voxels might lead to too few voxels and result in statistically insignificant descriptions of canopy features [23,40].In this study, several factors that affected the selection of the optimal voxel size in estimating LAI, using LiDAR data based on the voxel-based method, were explored.Voxel size is a key parameter pertaining to the scale of the crop's structural parameter estimates for the physical dimensions of the canopy components [41].Su et al. [22] determined the optimal voxel size for the two phases of data (10 cm in 2011 and 5 cm in 2014), respectively, and concluded that the difference between the two optimal voxel size probably resulted from the different morphologies of the maize, in these two years, so the optimal voxel size depended on the maize stalk and the leaf dimensions, in addition to dependence on the laser beam's diameter and the range and resolution of the laser scanner.After comparing and analyzing the point cloud densities of the routes in Table 1, we concluded that the factors affecting the selection of the optimal voxel size should include the point cloud density, the flight height, flight speed, the laser beam's diameter, the range and resolution of the laser scanner, and the structure of the maize canopy, which were all related to the point cloud density.To sum up, we concluded that the optimal voxel size depended on the characteristics of the LiDAR instrument used (such as incidence angle, flight height, flight speed, laser beam's diameter, the range and resolution of the laser scanner, and so forth) and experimental area characteristics (such as canopy structure and plant height).In addition, there were many approximate assumptions in the formula of the voxel-based method, which reduced the accuracy of the LAI estimation.Further studies, which take these variables into account, needs to be undertaken.

Conclusions
In this study, three flights of multi-route and five planting densities LiDAR dataset were acquired, using a UAV-mounted RIEGL VUX-1 laser scanner, at the altitude of 15 m.The novelty of our results lies in the occlusion effect of the maize LAI inversion, which was quantitatively analyzed from the aspects of the different layers and planting densities, analysis of the relationship between route direction and ridge direction, and the determination of the optimal voxel size.The results showed that: (1) The optimal voxel size was approximately 1.7 to 2.3 times of the average ground point distance.In addition, the optimal voxel size was related to the density of the point cloud; as the density of the point cloud increased, the optimal voxel size decreased, within a certain range.
(2) The direction of flight, perpendicular to the maize planting ridge, was better than that parallel to the maize planting ridge (R 22 -perpendicular to the maize planting ridge: NRMSE and R 2 were 4.6% and 0.96, respectively; R 26 -parallel to the maize planting ridge: NRMSE and R 2 were 17.7% and 0.84, respectively).
(3) The occlusion effect decreased as the planting density decreased, along with the NRMSE (21.4%, 17.7%, 15.7%, 10.8%, and 49.5% for densities D1D5 of the upper layer).The effect of occlusion on the upper layer was the lowest (10.8%), followed by the middle layer (12.4%), and the lower layer had the greatest influence (42.8%, which was an abnormal value due to the effect of wheat stubble).When the planting density was 27,495 plants/ha, the occlusion effect had the least impact (the upper, middle, and lower layers were 10.8%, 12.4%, and 42.8%), but with the planting density increasing to the highest densities, the NRMSE increased about 10%-30% for the upper, middle, and lower layers.
This study quantitatively analyzed the leaf occlusion effect and provided a guideline for obtaining the inner structural information of crop plants and their canopy, in agricultural remote sensing applications.The method developed in this study could be useful for early nutrition diagnosis and breeding research.Further research is needed to study the factors that affect the accuracy of the inversion of maize parameters, to expand its applicability.

Figure 2 .
Figure 2. Three flights of point cloud data.The figure shows the R13 route data, which is the point cloud data of 15 meters in the north-south direction and 0.625 m in the east-west direction of the experimental area, that is, a column of corn data planted on each flight.

Figure 1 . 16 Figure 1 .
Figure 1.The illustration of the experimental design.(a1,b2) Positional relationship between the experimental areas and each air route of the silking and blister stages.(a2,b1) Design of the different planting densities and repeat groups at two time points.(c) Diagram of maize stratification.(d) Schematic diagram of the incidence angle of each route and the setting of different densities.Note: R 11 , R 12 , R 13 , and R 14 represent the incidence angles of 69 • to 77 • , -64 • to -32 • , -34 • to -29 • , and 60 • to 74 • ; D1, D2, D3, D4, and D5 represent planting densities of 88,050 plants/ha, 64,005 plants/ha, 43,995 plants/ha, 27,495 plants/ha, and 7995 plants/ha.The same as below.

Figure 2 .
Figure 2. Three flights of point cloud data.The figure shows the R13 route data, which is the point cloud data of 15 meters in the north-south direction and 0.625 m in the east-west direction of the experimental area, that is, a column of corn data planted on each flight.

Figure 2 .
Figure 2. Three flights of point cloud data.The figure shows the R 13 route data, which is the point cloud data of 15 meters in the north-south direction and 0.625 m in the east-west direction of the experimental area, that is, a column of corn data planted on each flight. • ) and R 18 (68 • -76 • ), R 12 (−64 • to −32 • ) and R 17 (-60 • to −25 • ), R 13 (−34 • to 29 • ) and R 16 (−54 • to 0 • ), R 14 (60 • to 74 • ) and R 15 (63 • to 73 • ), R 22 (−12 • to −5 • ) and R 26 (3 • to 18 • ), and R 23 (8 • to 32 • ) and R 25 (−30 • to 0 • ), where "to" refers to the UAV flight direction; from west to east and north to south was indicated as positive, and the opposite was indicated as negative.Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 16 data was LAS, and the flight parameters of the two time points were same.The incidence angle range of the different routes obtained and the point cloud density of each route [point cloud number in the experimental area / (length × width)] are shown in Table

Figure 3 .
Figure 3.The Unmanned Aerial Vehicle-Light detection and ranging (UAV-LiDAR) obtained the data by the way of pendulum scanning and the distribution of laser foot point.

Figure 3 .
Figure 3.The Unmanned Aerial Vehicle-Light detection and ranging (UAV-LiDAR) obtained the data by the way of pendulum scanning and the distribution of laser foot point.

Figure 4 .
Figure 4. Trajectory calculation, raw laser data processing, and the preprocessing of the point cloud data.

Figure 4 .
Figure 4. Trajectory calculation, raw laser data processing, and the preprocessing of the point cloud data.

16 Figure 5 .
Figure 5.The point cloud data of the upper, middle, and lower layers of maize were obtained from the data of three flights.

Figure 5 .
Figure 5.The point cloud data of the upper, middle, and lower layers of maize were obtained from the data of three flights.

16 Figure 5 .
Figure 5.The point cloud data of the upper, middle, and lower layers of maize were obtained from the data of three flights.

Figure 7 .
Figure 7. LAI inversion accuracy comparison of 15 routes with different voxel sizes in the upper layer.R 2 and normalized root-mean-square error (NRMSE) are the accuracy comparison and analysis of the measured LAI and the point cloud data of the 15 plots.

Figure 7 .
Figure 7. LAI inversion accuracy comparison of 15 routes with different voxel sizes in the upper layer.R 2 and normalized root-mean-square error (NRMSE) are the accuracy comparison and analysis of the measured LAI and the point cloud data of the 15 plots.

Figure 8 .
Figure 8. Optimal voxel size corresponding to the cloud densities of different routes.

Figure 9 .
Figure 9. LAI inversion accuracy comparison of different routes with different voxel size of the whole plant.R 2 and NRMSE are the accuracy comparison and analysis of the measured LAI and the point cloud data of the 15 plots. .

Figure 8 .
Figure 8. Optimal voxel size corresponding to the cloud densities of different routes.

Figure 8 .
Figure 8. Optimal voxel size corresponding to the cloud densities of different routes.

Figure 9 .
Figure 9. LAI inversion accuracy comparison of different routes with different voxel size of the whole plant.R 2 and NRMSE are the accuracy comparison and analysis of the measured LAI and the point cloud data of the 15 plots. .

Figure 9 .
Figure 9. LAI inversion accuracy comparison of different routes with different voxel size of the whole plant.R 2 and NRMSE are the accuracy comparison and analysis of the measured LAI and the point cloud data of the 15 plots.

Figure 10 .
Figure 10.LAI inversion accuracy comparison for the different planting densities and layers.NRMSE is the accuracy comparison and analysis of the measured LAI and the point cloud data of 6 plots with 5 different planting densities.Note: D1, D2, D3, D4, and D5 represent the planting densities of 88,050 plants/ha, 64,005 plants/ha, 43,995 plants/ha, 27,495 plants/ha, and 7,995 plants/ha.

Figure 10 .
Figure 10.LAI inversion accuracy comparison for the different planting densities and layers.NRMSE is the accuracy comparison and analysis of the measured LAI and the point cloud data of 6 plots with 5 different planting densities.Note: D1, D2, D3, D4, and D5 represent the planting densities of 88,050 plants/ha, 64,005 plants/ha, 43,995 plants/ha, 27,495 plants/ha, and 7995 plants/ha.

Table 1 .
Incidence angle data and point cloud densities of different routes.Direction Route Angle (°) Point cloud density (pts/m 2 )

Table 2 .
The measured leaf area index (LAI) and LAI of different layers and planting densities.

Table 3 .
The optimal voxel size of each route.

Table 3 .
The optimal voxel size of each route.

Table 4 .
Comparison of LAI inversion accuracy of the different routes in the optimal voxel sizes.

Table 5 .
Comparison of different incidence angles.