Effect of layer thickness and voxel size inversion on leaf area density based on the voxel-based canopy profiling method

Yan Chen·Jian Liu·Xiong Yao·Yangbo Deng·Zhenbang Hao·Lingchen Lin·Nankun Wu·Kunyong Yu

Abstract Voxel-based canopy profliing is commonly used to determine small-scale leaf area.Layer thickness and voxel size impact accuracy when using this method.Here，we determined the optimal combination of layer thickness and voxel size to estimate leaf area density accurately.Terrestrial LiDAR Stonex X300 was used to generate point cloud data for Masson pines (Pinus massoniana).The canopy layer was stratified into 0.10-1.00-m-thick layers，while voxel size was 0.01-0.10 m.The leaf area density of individual trees was estimated using leaf area indices for the upper，middle，and lower canopy and the overall canopy.The true leaf area index，obtained by layered harvesting，was used to verify the inversion results.Leaf area density was inverted by nine combinations of layer thickness and voxel size.The average relative accuracy and mean estimated accuracy of these combined inversion results exceeded 80%.When layer thickness was 1.00 m and voxel size 0.05 m，inversion was closest to the true value.The average relative accuracy was 92.58%，mean estimated accuracy 98.00%，and root mean square error 0.17.The combination of leaf area density and index was accurately retrieved.In conclusion，nondestructive voxel-based canopy profiling proved suitable for inverting the leaf area density of Masson pine in Hetian Town，F(xiàn)ujian Province.

Keywords Terrestrial LiDAR·Leaf area density·Pinus massoniana ·Voxel-based canopy profiling method·Layer thickness·Voxel size

Introduction

Carbon storage，water-use efficiency，and total photosynthetic contribution rate of leaves differ at different heights above the ground.The vertical distribution of leaf area provides a theoretical basis on carbon sequestration and water consumption in forest ecosystems，as well as technical support for forest ecosystem management and maintenance(Zhang et al.2009).Leaf area density reflects the distribution of leaf area at different heights above the ground，and depends on certain factors，including species，growth stage，and environmental disturbance.In particular，leaf area density is an important input parameter for ecological process models that simulate the role of terrestrial vegetation in the global carbon cycle (Treuhaft et al.2 002).Therefore，the rapid and accurate determination of leaf area density is important for research on ecosystems and carbon cycles.

Leaf area density is measured either directly or indirectly.For direct measurements，layered sampling is mainly used(Li and Wang 1995).Compared to indirect measurements，data derived from direct measurements of leaf area density are more accurate.However，direct measurements are timeconsuming，laborious，destructive，and，although robust，are not suitable on a large scale.Indirect measurements have the advantage of being highly efficient and far less damaging and can also be used to estimate leaf area density at large scales.Many methods exist for indirectly measuring leaf area density (Ryu et al.2010;Yao et al.2011;Mac-Farlane et al.2019).Most of these methods rely on various instruments and software，including hand scanners，laser optic apparatus，portable scanning planimeters，and image analysis software.However，traditional optical instruments that determine leaf area density cannot distinguish between leaves and branches (Neumann et al.1989;Whitford et al.1995;Chen 1996) and are affected by various factors such as tree species (Deblonde et al.1994)，zenith angle (Jie et al.2009)，azimuth angle (Jie et al.2009)，and season.The correction factor does not adequately correct the measurement results.With the development of active remote sensing technology，such as LiDAR，the 3D structural characteristics of forest canopies can now be obtained by establishing 3D coordinate system data and digital images (MacFarlane et al.2007;Hu et al.2014，2018)，which are widely used(Jupp et al.2009;Tang et al.2014;Hu et al.2018).

Leaf area density measured using LiDAR is commonly derived using regression，gap probability，and voxel methods.For the regression method，parameters for a single tree(such as tree height，breast diameter，and crown width) are obtained by LiDAR to derive a regression equation for leaf area density.This method efficiently and accurately estimates leaf area density (Roberts et al.2005).However，the regression equation can only be applied to the examined tree species (Richardson et al.2009).The gap probability method estimates leaf area density based on gap probability and the angle distribution law (Bao et al.2018).This method effectively distinguishes non-intercepted and intercepted laser pulses.The multiple reflections of the laser pulse and the relatively large laser pulse size allow the canopy to be recognized.However，issues exist with small and mediumsized gaps (Lovell et al.2003;Seidel et al.2012).The voxel method estimates leaf area density based on contact frequency.The target tree is scanned from multiple sites to reduce errors caused by uneven laser penetration (Van der Zande et al.2006;Zhang et al.2016).The obtained point cloud data is voxelized and，then，the contact frequency of different layers is calculated.The leaf area index of different layers is subsequently obtained and accumulated to determine the leaf area density of the canopy.Occlusion within the canopy generates results with low estimates (Zhao et al.2015)，with estimated voxel size also having an impact.

The Masson pine (Pinus massonianaLamb.;Pinaceae)is a pioneer tree species during forest succession.It is characterized by its resistance to drought and barren conditions，ease of planting，high survival rate，and wide use for pine wood.It is one of the most representative tree species of southern China.Leaf area density generally represents forest quality (McFadyen et al.2009;Ghoreishi et al.2012;Lin et al.2021).It is hard to estimate the leaf area density of Masson pine using LiDAR because its leaves are not flat.As a result，the effect of this type of leaf morphology on leafinterception remains poorly understood.

The voxel-based canopy analysis method was first proposed by Hosoi and Omasa based on the frequency of contact.In brief，point cloud data obtained from multiple sites are registered in the same coordinate system.Using these data，a cubic element model is established，the contact frequency of each layer of canopy needles is estimated using the laser，and a leafinclination correction factor is constructed to obtain the leaf area density.Using the voxel canopy analysis method in the present study，we divided the canopy of single trees into layers of different horizontal thickness and combined different voxel sizes to estimate the leaf area density of Masson pine.The projected area of the horizontal layer represented the largest projected area on the z-axis.Thus，our main objectives were to (1) estimate the area density of coniferous needles with ground-based LiDAR and (2) determine the optimal voxel size and layer thickness in the estimation process.The results of the present study are expected to facilitate nondestructive and accurate estimations of leaf area density of pine trees.

Materials and methods

Study area

The study area was located in Changting County，F(xiàn)ujian Province，China (25°35′-25°46′ N，116°16′-116°30′ E)(Fig.1).The study area is bordered by Nanshan in the east，Tufang and Putian in the south，Xinqiao in the north，and Cewu in the west.The study area has a subtropical monsoon humid climate，with a mean annual temperature of 18.3°C.The mean annual rainfall is 1700.0 mm.The terrain is broken with soil in the middle of valley.The soil is mainly red loam and sandy loam and has low resistance to erosion.The forest is mostly composed of Masson pine andDicranopteris dichotoma(Thunb.) Berhn.

Data collection

Specific leaf area acquisition

Twenty-six Masson pine trees were selected from the study site as the research objects (Fig.1).True leaf area density of the canopy was obtained using the stratified sampling method.The target trees were felled，and whole leaves were harvested from the bottom up.All leaves were removed and weighed to obtain fresh leaf mass.Leaf area was determined by taking 1000 g leaves using the crossing method (Ren and Peng 1997).Specific leaf mass was calculated by dividing the leaf area with 1000 g.We calculated the leaf area index according to the following formula (Ren and Peng 1997;Tan and Zhao 2008):

whereWis the total leaf mass of the target tree，Kis the specific leaf mass，andSis the vertical projection area of the target tree.

Point cloud data acquisition

We obtained point cloud data for the 26 pine samples using a Stonex X300 laser scanner，which is a pulsed-static 3D laser scanner.Combined with a high-resolution digital camera，this scanner provides 3D point cloud data，and color images of a target scene quickly and nondestructively.The technical parameters are shown in Table 1.

Table 1 Technical parameters of Stonex X300

Table 2 Measured model parameters of the target tree

To reduce the occlusion of foliage and obtain complete information about the canopy，three observation sites were set up per plot.The distribution of the stations is shown in Fig.2.To unify the point cloud data merged from the three stations independently under a unified coordinate system，we deployed two black and white targets of different heights.One of the three stations was selected for a panoramic scan of the target tree.The horizontal field of view of the scanning area of the other two stations was subject to partial overlap between stations.Laser scanning was carried out when there was no cloud cover or wind.

Fig.2 Stations from where individual tree data were acquired.·station location

Data pre-processing

Data pre-processing was mainly completed for point cloud registration and denoising.Point cloud registration involved registering cloud data from different measurement sites based on the common target method.In brief，the target was accurately scanned from all stations during the scanning process.The target center was the common point between two stations.Target registration involved identifying the common point between stations (**Lin 2019).In general，good registration results can be obtained based on the common target registration method.The registration error was 0-1 mm，which met the accuracy requirements of the study.Point cloud denoising was achieved by removing nontarget point clouds that may have occurred due to errors caused by noise signals from the instruments.Erroneous point clouds can be removed by visual interpretation of the target using Cyclone7.0 software (Leica Geosystems AG，Heerbrugg，Switzerland).Point clouds generated by instrument noise can be deleted using an algorithm.We used the method of Rusu et al.(2008) to remove error points.

Canopy leaf area density inversion method

Removal of nonphotosynthetic tissue

Eliminating nonphotosynthetic tissue point clouds (such as branches) is important to extract leaf point clouds and accurately estimate leaf area density.Nonphotosynthetic tissue was manually removed by visual interpretation (Fig.3).

Fig.3 Point cloud data after removing part of nonphotosynthetic tissue

Constructing a voxel model

Avoxel in a 3D image corresponds to a pixel in a 2D image.Based on the range of point cloud data after pre-processing，the minimum values (Xmin，Ymin，andZmin) of coordinatesX，Y，andZwere taken as the starting point.Voxel size was used as the step size to determine the corresponding units in the point element coordinate system.Voxel size was determined by the lengthl，widthw，and heighthof the voxel.Voxels with equall，w，andhwere cubes.The point cloud was divided intoNl，NwandNhindividual elements，where(Nilsson 1996).The coordinate values after the point cloud was normalized were calculated as

where Int is a function to round offthe coordinates to one decimal place to the nearest integer，andi，j，kis a voxel coordinate of point cloudX，Y，Z.In the present study，the voxel sizelwhwas consistent with the dot spacing used for scanning.

The voxel value was determined by judging the number of laser points contained in the voxel.If the number of laser points was greater than or equal to 1，then the laser beam was intercepted，and a voxel value of 1 was assigned;otherwise，a voxel value of 0 was assigned.

Leaf area density computational model

The voxel-based canopy analysis method was used to calculate the leaf area density (LAD) with the following formula:

where LAD is leaf area density，ΔHis the horizontal thickness layer，θis the average value of the zenith angle of all incident laser beams for the horizontal thickness layer ΔHrange，nI(k) is the number of voxels in thekth layer intercepting the laser，nP(k) is the number of voxels penetrated by the laser on the kth layer，is used to indicate the frequency at which the laser contacts the canopy，is the correction factor of the needle inclination angle and the direction of the laser beam，andG(θ) is an average projection with a uniform azimuth distribution of the needles and a unit leaf area in a plane perpendicular to the direction of the laser beam.As defined by Hosoi and Omasa (2006):

wherex=cos-1(cosθcosθL)，θis the eccentric angle of the incident laser beam，θLis the leaf-inclination angle，φLis the normal azimuth angle of the leaf surface，andnBandnLare the laser beam direction and the expected inclination of the leaf surface，respectively.The zenith angle and the leaf dip angle were obtained from the laser point cloud data.

As described by Nilsson (1996)，the actual measured distribution of leaf-inclination angle was expressed as

whereqis the category of different leaf dip angles，Nqis the total number of categories，if the spacing is 5°，Nqis 18，θL(q) is theq-type intermediate leaf dip angle value，andg(q)is theq-type leaf dip angle.The distribution is the ratio of the leaf area of theqclass to the total leaf area.

Contact frequency calculation

The laser-canopy contact frequency is calculated by counting the number of intercepting laser voxelsnI(k) and the number of light penetrating voxelsnP(k) in each horizontal thickness layer.If a 3D voxel model is constructed according to the maximum and minimum range of the regional point cloud，due to the irregularity of the canopy structure，invalid voxels outside the canopy occur，which cannot be considered as voxels penetrated by the laser beam.Therefore，it is important to validate the canopy boundary and eliminate invalid voxels before calculating the contact frequency.

Graham’s scan，a simple and efficient 2D convex hull algorithm，was used to determine the outer contour of each canopy (Graham 1972).The point set convex hull refers to a minimum convex polygon.The points satisfying the point set are inside the polygon or on the sides of the polygon.The Graham scan convex hull algorithm included the following criteria:

(1) The point with the smallest coordinate on theY-axis was the base point，and the starting pointp0of the convex hull was used.If there were multiple coordinate points with the same minimumY-value，then the point with the smallestX-coordinate was selected as the reference point.

(2) The point at which the vector formed wasp iand was rotated in a counterclockwise direction to becomep i+1.If the newly added point satisfiedp i+1on the left side of the vector＜p i-1，p i＞，then scanning was continued.Ifp i+1was not satisfied on the left side of the vector＜pi-1，p i＞，thenp iwas not the convex cullet vertex and was deleted until a point was found that meets the stated conditions.

(3) Whenp i=p0，the graph was closed，and the convex hull was completed.

A top projection was performed on the acquired point cloud data to obtain a canopy projection map.We used Graham’s scan algorithm to obtain the horizontal thickness layer of the outer contour of the canopy.Figure 4 shows the needles intercepting laser beams in a horizontal thickness layer，reflecting the horizontal thickness of needle coverage.The convex hull vertices were connected to obtain a convex polygon representing the canopy range.

Fig.4 Canopy boundary contour

Correction for leafinclination

The angle of inclination of the needles is the angle between the normal direction of the needle surface and the direction of the zenith.The distribution of the inclination angle of a needle directly determines the interception of the incident laser beam by the vegetation canopy.The correction for needle inclinationsis the correction factor of the needle inclination angle and laser beam direction.Equations (7) and(8) were used to calculate the angle of needle inclination and the zenith angle of the laser beam.

The ground-based 3D laser Stonex X300 has a scanning accuracy of＜6 mm，while the width of Masson pine needles is＜6 mm;consequently，a single pine needle point cloud could not be obtained (Fig.5).However，a clustered pine needle point cloud could be obtained，and the inclination angle of the same cluster-shaped pine needles was similar.The average pine needle plane was estimated by randomly selecting pine needles at different heights，then using the eigenvalue method to fit the plane.The angle between the average direction and zenith direction was the leafinclination angle (Abbas et al.2012).

Fig.5 Leaf point cloud of the Masson pine (Pinus massoniana)

wherex，y，zis the Cartesian coordinate of the point，r，α，βis the polar coordinate of the point，βis the angle between the scanning direction of the instrument and the horizontal line.The elevation angle was positive，the depression angle was negative，and the Cartesian coordinates could be converted to each other.The apex angle wasθ.

Precision verification

Direct and indirect methods for measuring leaf area density are available.An indirect measurement was used to compare and verify the leaf area index，which was obtained by accumulating the calculated leaf area density with the canopy leaf area index using a leaf area index measuring instrument.The method was relatively easy to implement;however，it failed to express the leaf area density inversion for horizontal thickness layers accurately.The direct method involved collecting the leaves layer by layer and measuring leaf area.This method reflected the leaf area density of the horizontal thickness layers accurately;however，this approach is timeconsuming and destructive.

The canopy was divided into three layers:upper，middle，and lower.All leaves were obtained from the bottom to the top.The cross-section method was used to obtain 1000 g of leaves from each layer to measure leaf area.The leaf area index was subsequently calculated.Root mean square error(RMSE) (10)，average relative accuracy (RA) (11)，and mean estimated accuracy (RM) (12) were used to evaluate the accuracy of the inversion.

wherenis the number of samples，yiis the actual measured value，andis the estimated value.

Results

Leaf area density inversion

Contact frequency

The overall contact frequency changed with layer thickness in a consistent manner (Fig.6).The smaller the layer thickness，the more detailed is the description of canopy contact.Figure 7 shows the contact frequency distribution under different voxel sizes at different height levels.When the layer thickness was the same，the contact frequency increased with increasing voxel size because the pore area in the canopy was gradually replaced by the leaf area (Fig.6).At different voxel sizes and different layer thickness，the contact frequency of the canopy flat layer differed (Figs.6 and 7).

Leafinclination estimation

One hundred leaf point clouds were randomly selected at different heights in the canopy.Based on these data，the plane was estimated by the eigenvalue method，and the leafinclination angle was calculated.The leafinclination angle was calculated as 0-90°.The distribution range of the leafinclination angle was divided into nine intervals covering 10°each.We calculated the angular probability of each interval to obtain a leaf probability distribution map.Figure 8 shows that the leafinclination angles of Masson pine were distributed at a range of 10°-90° and were concentrated at a range of 30°-70°.The average leafinclination angle was 49.90°.

Fig.7 Contact frequency for different voxel sizes (0.5 m layer thickness as an example)

Fig.8 Distribution of leafinclination angle probabilities

Zenith angle and correction factor

The average zenith angle was obtained by converting the rectangular coordinates to polar coordinates using Eq.9.The distribution range of the zenith angle was relatively large (33.45°-79.10°).The large range of zenith angles was attributed to the frequently short distance between the terrestrial LiDAR and target tree and to the high forest canopy during scanning.The correction factors for different levels were calculated from Eqs.7 and 8.Whenθthe correction factor was 1.07;however，whenθthe correction factor changed with the zenith angle.

Leaf area density inversion

The contact frequency and the correction factor were used to calculate the leaf area density with Eq.5.Leaf area density at the same height differed for the thickness value of each layer (Fig.9).The height of the maximum leaf area index differed with different layer thickness.The correction factors calculated from the zenith angle for each layer separately，along with the correction factors calculated from the average zenith angle of the canopy，also impacted estimates of leaf area density.Compared with the correction factors for layers，the average correction factors for canopy had a greater overall leaf area density，and the distribution at different heights was more uniform.Voxel size had a certain effect on estimates of leaf area density (Fig.10).The contact frequency of each horizontal layer increased with increasing voxel size.The pattern of leaf area density changed with height.The leaf area density of each layer was similar for different correction factors.

Fig.9 Effect of correction factors and average correction factors on inversion results with different layer thickness (0.05 m voxel as an example).Black line shows leaf area density (LAD) for each horizontal layer correction factor.Red line shows LAD for the average correction factor

Fig.10 Effect of correction factors and average correction factors on inversion results with different voxel sizes (0.50 m layer thickness as an example).Black line shows leaf area density (LAD) for each horizontal layer correction factor.Red line shows LAD for average correction factor

Precision verification

The accumulated real leaf surface index from the layered harvesting method was 0.65-3.32 m2/m2(Table 2).Root mean square error (RMSE)，average relative accuracy (RA)，and mean estimated accuracy (RM) were calculated using Eqs.10，11，and 12，respectively.There were nine combinations with an average relative accuracy higher than 80.00%(Table 3).When the stratified thickness was 1 m and voxels 0.05 m，the inversion accuracy of the stratified correction factors and the average canopy correction factor were highest.

Discussion

Method applicability

The leaf area index (LAI) was measured by the LAI-2200C plant canopy analyzer，and can be impacted by many factors (Kuusk 2016;Pearse et al.2016).Obtaining the LAI with this instrument provides an effective value of LAI.If the effective LAI is used to verify the true LAI，then the verification accuracy differs significantly from the accuracy determined using true LAI (Wang 2015;Wang et al.2016).Therefore，estimated leaf area density was verified using a destructive layered harvesting method that obtained the true LAI.

The overall LAI was used to verify the results.The optimal combination of voxel size and layer thickness for estimating Masson pine was found.The leaf area density estimated from the nine combinations was the average canopy LAI，and the estimation accuracy was high (Table 2).The LAI of Masson pine ranged from 0.65 to 3.32 m2/m2，values that were relatively low compared with those reported in other studies (Norman and Jarvis 1974;Deblonde et al.1994;Chen et al.1997;Yan et al.2019a，b).The current study area was located in Hetian Town，Changting County，Longyan City，F(xiàn)ujian Province.Masson pine was planted in this region to combat a serious soil erosion problem;however，it exhibits poor growth and development，with a low LAI.

Table 3 Accuracy verification table of overall leaf area density under nine different combinations

To verify whether the leaf surface density at the same height above the ground was consistent with the measured value，we divided the estimated leaf area density into upper，middle，and lower parts for verif ciation.Unfortunately，layerwise accuracy was below 80.00%.The error in leaf area density calculated by the average correction factor was much larger than that of leaf area density calculated by the layered correction factor.This phenomenon arose because the average correction factor cannot accurately correct the leafinclination angle and laser incident zenith angle at different canopy heights.However，because the leaf area index is obtained by the accumulation of leaf area density，it is mutually “compensated” at different heights，generating higher accuracy in the obtained leaf area index.The layered correction factor accurately corrects the leafinclination angle and the laser incident zenith angle at different heights;however，this accuracy is still not high enough for layered verification.This result was attributed to canopy layer thickness being 0.40 m，0.50 m，0.60 m，0.70 m，0.90 m，and 1.00 m when verifying accuracy.This layer thickness is large.When the canopy is divided into upper，middle，and lower layers for verification，the cumulative height of the three layers differed from that of the actual canopy.Consequently，the accuracy of the layered estimates was lower，whereas the accuracy of the combined estimate was higher.Thus，the voxel-based canopy profiling method was suitable for estimating the leaf area of Masson pine.Future studies verifying leaf area density at different heights should optimize layer thickness to maximize the accuracy of estimates.

Effects of different voxel sizes and layer thickness

Leaf area density，the sum of the areas of photosynthetic unilateral leaves in a unit volume at a certain canopy height，is used to characterize differences in the leaf area of the different layers and was estimated at various heights.When the leaf area at different heights is estimated based on the voxel canopy analysis method，the layer thickness and voxel size impact the contact frequency and，therefore，the estimation results.When voxel size or layer thickness is large，the small gap between the canopy layers is overlooked，and the estimated results are larger than the actual value.When the layer thickness or voxel size is small，the estimated result is closer to the true leaf area density (Wang 2015)，with more calculations.When the voxel size is very small，the conditions assumed by Poisson’s law cannot be met (Béland et al.2015).Therefore，appropriate voxel size and layer thickness are needed to estimate leaf area density accurately.Béland et al.(2011) used 0.05-2.00 m voxel sizes to estimate the area of broad leaves，whereas Wang et al.(2016，2018) used 0.50 m and 0.20 m，as layer thicknesses to estimate leaf area density.On the basis of these studies，here we delineated the best combination of voxel size and layer thickness to estimate leaf area density for Masson pine by evaluating a range of voxel sizes from 0.01 to 0.10 m and layer thicknesses from 0.10 m to 1.00 m.Previous studies reported that small gaps in voxel size in the canopy are excluded when voxel or layer thickness is large，generating a more general canopy expression;however，the frequency of exposure is overestimated.Conversely，the laser can provide detailed information on the internal structure of the canopy.We found that the leaf area density was accurately estimated when layer thickness was 1.00 m and voxel size was 0.05 m.

The combination of estimated leaf area density of Masson pine in our study included layer thicknesses of 0.40，0.50，0.60，0.90，and 1.00 m，and voxel sizes of 0.01，0.02，0.03，0.04，0.05，0.06，and 0.07 m.Our study area was the same as that investigated by Deng et al.(2017).However，the authors used regression analysis to estimate leaf area and reported a layer thickness of 0.10 m and voxel size of 0.05 m，which contrasted with our results.This difference might be attributed to the effects of leafinclination and zenith angle used in the present study.Furthermore，our results were corrected by correction factors to obtain accurate leaf area density.By examining nine combinations of voxel and layer size，we showed that voxel size increased with layer thickness.Future studies should explore whether there is a relationship between voxel size and layer thickness and expand the research scope of the voxel-based canopy profiling method.

Effects of the correction factor

Based on the voxel-based canopy profiling method，different scholars have used different correction factors.Wang(2015) and Wang et al.(2016) documented fir a large error the correction coefficient calculated from the average zenith angle of the entire canopy.The canopy should be layered horizontally，and the correction factor is obtained based on the average zenith angle of each horizontal layer.Hosoi and Omasa (2006) and Dai (2019) found that，when the apex angle is 57.5°，the correction factor is 1.1 and that it is minimally affected by needle inclination.The zenith angle in the present study was considerably different from 57.5°.The average leafinclination angle obtained here was similar to that calculated by Deng et al.(2017).To explore the difference between the average correction factor and layer correction factor，the same combinations of voxel and layer thickness were used in the present study.Two different correction factors were used to estimate leaf area density and verify its accuracy.The current study showed that the layered correction factor estimates leaf area density more accurately than the average correction factor.Thus，when the estimation accuracy requirement is not high，the average correction factor can be used for calculation.If the accuracy requirement is high，then the layered correction factor should be used.However，calculating stratified correction factors is extremely time-consuming，so future studies should explore how to optimize that process.

Conclusions

Based on the voxel-based canopy profiling method，volume density and layer thickness of 26 Masson pine trees were evaluated as the key factors for estimating leaf area density.The true leaf area index obtained in the field was used to verify the estimation accuracy.This study verified the feasibility of estimating leaf area density of Masson pine based on the voxel-based canopy profiling method.The most accurate estimates of leaf area density were obtained when the voxel size was 0.05 m and layer thickness was 1.00 m.This method used terrestrial LiDAR to scan three sites and obtain target tree point cloud data.Visual interpretation was used to delete nonphotosynthetic tissues.The 2D convex hull method was used to calculate contact frequency and obtain leaf area density.This method eliminates the effect of nonphotosynthetic tissue on the results，and it eliminates the need to assume the spatial distribution，size，and shape of branches and leaves in the canopy.However，the optimal voxel size and layer thickness identified here was specific for Masson trees and must be separately calculated for tree species with different leaf types.

Author contributionsYC:methodology，data analysis，verification，writing-original draft;JL:review，editing;XY:concept，methodology;YD:survey，data analysis;LL:survey，data processing;ZH:survey，verification methodology;NW:data analysis;KY:grant writing，review，modification.