Modelling Three-Dimensional Spatiotemporal Distributions of Forest Photosynthetically Active Radiation Using UAV-Based Lidar Data

The three dimensional (3-D) spatiotemporal variations of forest photosynthetically active radiation (PAR) dictate the exchange rates of matter and energy in the carbon and water cycle processes between the plant-soil system and the atmosphere. It is still challenging to explicitly simulate spatial PAR values at any specific position within or under a discontinuous forest canopy. In this study, we propose a novel lidar-based approach to estimate both direct and diffuse forest PAR components from a 3-D perspective. An improved path length-based direct PAR estimation method was developed by incorporating the point density along a light transmission path, and we also obtained the diffuse PAR components using a point-based sky view analysis by assuming the anisotropic sky diffuse distribution. We compared the total PAR modelled using three light path length-based parameters with reference data measured by radiometers on a five-minute time scale during a daily solar course. Our results show that, in a discontinuous forest canopy, the effective path length is a feasible and powerful (R2 = 0.92, p < 0.01) parameter to capture the spatiotemporal variations of total PAR along a light transmission path with a mean bias of −53.04 μmol·m−2·s−1(−6.8%). Furthermore, incorporating point density and spatial distribution factors will further improve the final estimation accuracy (R2 = 0.97, p < 0.01). In the meantime, diffuse PAR tends to be overestimated by 17% at noon and underestimated by about 10% at sunrise and sunset periods by assuming the isotropic sky diffuse distribution. The proposed lidar-based 3-D PAR model will provide a solid foundation to various process-based eco-hydrological models for simulating plant physiological processes such as photosynthesis and evapotranspiration, intra-species competition and succession, and snowmelt dynamics purposes.


Introduction
The forest radiation regime is the basic driving factor for most physiological processes such as photosynthesis and respiration [1,2].The spatiotemporal distributions of forest photosynthetically active radiation (PAR, 400~700 nm) in the three dimensional (3-D) space has great effects on sub-canopy snowmelt dynamics [3,4], understory evapotranspiration [5,6], and the growth and succession of tree seedlings [7,8].It is determined by absorbing, scattering, and transmitting processes between solar PAR and foliage elements as it penetrates through a forest canopy [9,10].These processes are dominated by factors including leaf orientation and morphology, the distribution pattern of foliage elements, solar position, and topographic conditions [11,12].Extremely high heterogeneity could usually be observed for spatiotemporal distributions of forest PAR in both horizontal and vertical dimensions, especially at a sub-daily scale [13,14].
The spatiotemporal variations of PAR could be conventionally and directly measured by pyranometers [15,16] and indirectly observed by digital hemispherical photography (DHP) [17][18][19] at different locations within a forest plot.However, it is difficult to set up a pyranometer observation network or to obtain DHPs in large or inaccessible forested areas, and the high spatial heterogeneity of PAR in complex forest areas cannot be effectively represented by data at specific sampling points.Therefore, their applications can't be extended to larger spatial scales.To solve this issue, canopy light models became a powerful tool to investigate the 3-D spatiotemporal distributions of forest PAR [20].The spatially explicit 3-D radiative transfer models have been developed in the past decades for non-randomly distributed heterogeneous forest canopies.They included the Beer law-based radiative transfer model [21][22][23], and the ray tracing model using Monte-Carlo simulation [24][25][26].The Beer-Lambert law describes the exponential attenuation of solar radiation, and it is only applicable to situations where randomly distributed foliage elements can be assumed [27].The detailed 3-D structure information of forest canopies is a prerequisite input for these computer-based models.For example, some models represented the tree crowns using different geometric objects such as a cone or cylinder, and the direct light could be traced during the daily solar course [28][29][30].However, it is usually extremely time-consuming and labor-intensive to obtain forest canopy structural parameters such as tree height, diameter at breast height (DBH), crown size, and crown base height [14].Moreover, compared with the real 3-D forest scene, the artificially reconstructed forest scene usually fails to capture the real gap distribution patterns, which makes the modelled radiation at the specific location deviate from the actual PAR values [31,32].In addition, previous studies have shown that the estimation accuracy of aerial laser scanning (ALS)-derived forest canopy structural parameters will decrease slightly as data point density decreases [33,34].However, how the point cloud density affects the accuracy of ALS-based 3-D canopy light models remains unclear.
An accurate description of the tree structure is the key step to simulating the radiation regime of a forest canopy.With the development of light detection and ranging (lidar) technology, a detailed 3-D structure represented by lidar-based high-density point cloud data can be acquired efficiently.ALS has shown great potential for forest ecological applications at the landscape level [35,36].To investigate the radiation regime of forest canopies, several ALS-based 3-D canopy light models have been successfully developed over the past decades [32,37,38].Most of these models characterize the probability of direct or diffuse solar radiation penetrating through forest canopies based on various lidar-derived metrics such as the laser penetration index [39], leaf area density [40], point projection area on the ground surface [41], number of points around a solar ray [42], and synthetic hemispherical images derived from projected point cloud data [43,44].More detailed 3-D canopy light models have been developed based on terrestrial laser scanning (TLS) to simulate radiation transmittance by voxel-based canopy reconstruction [45,46].As a bridge linking ALS and TLS, unmanned aerial vehicle (UAV)-based lidar systems provide great flexibility in selecting data point density, as well as better delineation to the upper part of a canopy with top-down scanning [47,48].However, a 3-D forest canopy radiation model driven by high-density UAV-based lidar data remains a necessary development.
The PAR arriving at one point in the forest canopy mainly consists of direct components originating from the transmission and attenuation of beams and diffuse components derived from single and multiple scattering.Naturally, most canopy light models simulate them separately [12,39].In terms of direct solar radiation, the length of the transmission path, defined as the distance between the first contacted point and the target point [13,30], is a widely used parameter to determine the extent of its potential attenuation [14,29].It was found that a path-length based radiation simulation model is more flexible in terms of reproducing sub-daily and seasonal PAR variations [49].Moreover, existing research found that a spatially explicit 3-D ray trace model incorporating path length could simulate the spatiotemporal distribution patterns of light intensity in the forest floor and along the vertical gradient [14].However, it is usually not good enough to only consider the total path length to accurately characterize the light attenuation effect along its transmission path due to the non-random distribution of foliage elements within a forest canopy [50].The gaps between clumped foliage elements along a light transmission path will alter their attenuation ability for direct solar radiation.In addition, the point density and distribution pattern along a light transmission path will also affect the ability of a direct beam penetrating through a forest canopy [33,51].However, few studies have incorporated the gap size information into the path length-based method to characterize the light attenuation, and it is still unclear whether it is possible to improve the simulation accuracy of the light attenuation effect by incorporating the point density and distribution information.In terms of diffuse components, most research treats the diffuse light above the hemispherical sky as an isotropic distribution.The diffuse solar radiation reaching a specific location in a forest canopy is usually estimated based on the sky view factor or canopy closure [15,19,49].To assume the isotropic diffuse solar radiation will affect the estimated radiation regime in 3-D space [52,53], while the difference remains unexplored.
Therefore, in this study, we aimed to develop a spatially explicit forest canopy radiation model to investigate the spatiotemporal variations of forest PAR within and under a forest canopy.The specific objectives of this study are to: 1.
Develop and validate a spatially explicit lidar-based 3-D forest canopy PAR simulation model treating direct and diffuse solar radiation separately; 2.
Characterize the attenuation effects of direct solar PAR along its transmission path penetrating through a forest canopy, and 3.
Explore the spatiotemporal variations of forest PAR in 3-D space and the effects of point density on the estimation accuracy of forest PAR simulation.

Study Sites
Our study area is located in the Baima (BM) experimental forest site in the city of Nanjing (119 • 11 8"E, 31 • 36 51"N), China (Figure 1).The site is a planted forest site with a relatively flat area and different tree species, densities, and ages.We set up a medium density (0.15 trees / m 2 ) deciduous broadleaf circular forest plot with a radius of 15 m.The average tree height and crown size of this forest plot are 7 m and about 2~3 m, respectively.The effective leaf area index (LAIe) is about 1.5.The tree species of this forest plot is wheel wingnut (Cyclocarya paliurus) with an average spacing of 2 m between tree stems.diffuse solar radiation reaching a specific location in a forest canopy is usually estimated based on the sky view factor or canopy closure [15,19,49].To assume the isotropic diffuse solar radiation will affect the estimated radiation regime in 3-D space [52,53], while the difference remains unexplored.Therefore, in this study, we aimed to develop a spatially explicit forest canopy radiation model to investigate the spatiotemporal variations of forest PAR within and under a forest canopy.The specific objectives of this study are to: 1. Develop and validate a spatially explicit lidar-based 3-D forest canopy PAR simulation model treating direct and diffuse solar radiation separately; 2. Characterize the attenuation effects of direct solar PAR along its transmission path penetrating through a forest canopy, and 3. Explore the spatiotemporal variations of forest PAR in 3-D space and the effects of point density on the estimation accuracy of forest PAR simulation.

Study Sites
Our study area is located in the Baima (BM) experimental forest site in the city of Nanjing (119°11′8″E, 31°36′51″N), China (Figure 1).The site is a planted forest site with a relatively flat area and different tree species, densities, and ages.We set up a medium density (0.15 trees / m 2 ) deciduous broadleaf circular forest plot with a radius of 15 m.The average tree height and crown size of this forest plot are 7 m and about 2~3 m, respectively.The effective leaf area index (LAIe) is about 1.5.The tree species of this forest plot is wheel wingnut (Cyclocarya paliurus) with an average spacing of 2 m between tree stems.

UAV Lidar Data
We acquired the UAV-based lidar data in the BM site on July 26, 2017 on a clear day with the DJI M600 pro (DJI Technology Co. Ltd, Shenzhen, China) drone mounted with a Velodyne VLP-16 (Velodyne LiDAR Inc.San Jose, USA) lidar sensor coupled with a high precision inertial measurement unit (IMU).The scanning field of view (FOV) was set as ±15 • with the laser wavelength of 903 nm and beam divergence of 0.18 • (3.0 mrad) during the flight.The maximum scan rate of the lidar sensor was 300,000 points/s.Two cross flight paths (i.e., East-West and North-South) were designed to acquire relatively comprehensive forest point cloud data, which resulted in high-density (about 500 points/sqm) 3-D forest point cloud data (Figure 1b).

TLS Data
In the meantime, we also collected terrestrial lidar data using the Leica Scan Station 2 (Leica Geosystem AG, St. Gallen, Switzerland) at five different scanning locations with one center position and four corner positions, respectively.At each station, we scanned the plot using hemispherical mode (i.e., horizontal: 0~360 • ; vertical: −45~90 • ) by setting sampling spacing as 1 cm @ 10 m (1 cm at 10 m away from TLS station) with the laser wavelength of 532 nm.All TLS data from five different stations were registered into a comprehensive forest plot TLS point cloud.

Lidar Data Pre-Processing
In this study, we filtered ground points for both UAV and TLS lidar data using the tool 'CSF' in the 'CloudCompare' software package (version 2.6.2) (GPL software and available from http: //cloudcompare.org).Then, the digital elevation model (DEM) and canopy height model (CHM) with the spatial resolution of 0.1 m were interpolated in ArcGIS software using ground points from UAV lidar data.In addition, after removing the noise points of TLS data, we removed all points lower than each PAR observation location to estimate the sky view factor for diffuse solar radiation calculation purposes.Accurately locating the observation points in 3-D both UAV-and TLS-based lidar data is a key step to assure the effectiveness of comparison between the lidar-and field-based PAR values.All nine observation locations were precisely captured by TLS data.By visually selecting two control points in both UAV-and TLS-lidar data space, we first calculated distances from each observation location to these two control points.Then, by using the distance intersection method we successfully located nine observation locations in the UAV-lidar data space (Table 1).

Field-Based PAR Measurements
Nine different light quantum sensors recorded photosynthetic photon flux density (PPFD) (µmol•m −2 •s −1 ) with the average logging time interval of 5 minutes and heights of 1~2 m above the ground at different locations at the BM site from July 27 to July 28, 2017.In the meantime, we also measured the total (PPFD total ) and diffuse PAR (PPFD 0,di f ) above forest canopy using the BF-5 sunshine sensor (Delta-T Devices Ltd., Cambridge, UK) placed on a roof 20 meters above the ground.By doing so, we obtained the direct PAR (PPFD 0,dir ) by subtracting PPFD 0,di f from PPFD 0,total .

PAR Measurements Normalization
To compare the PAR measurements from nine different light sensors, we conducted the data calibration process using a standard reference light quantum sensor (LI-190R, LI-COR, Inc., Lincoln, Nebraska) under the same light condition and indoor environment.By sampling the PAR readings from the to-be-calibrated sensor with the reference one at a step of 50 µmol•m −2 •s −1 within the range from 0 to 2500 µmol•m −2 •s −1 , we calculated the difference between two readings.The quadratic polynomial statistical model best explains the variations of reference PAR measurements with the to-be-calibrated measurements with significant correlations (R 2 = 0.99 and P < 0.01) and is expressed as: where ∆PAR i is the difference of PPFD between the i-th pyranometer and LI-190R, PAR LI−COR is the corresponding PPFD of LI-190R, and a i , b i and c i are calibration coefficients for the i-th pyranometer as shown in Table 1.

Model Development
We developed a lidar-based radiation model based on the following assumptions: (1) The total solar radiation (PPFD total ) reaching a specific location within or under a forest canopy in 3-D space consists of both the direct PAR (PPFD dir ) and diffuse PAR (PPFD di f ) components.(2) The solar radiation reaching a specific location is mainly determined by the non-random distributions of foliage elements within the 3-D space of a forest canopy which is represented by UAV-lidar data in this study.The solar radiation absorbed by leaves was ignored in this study.We calculated the direct PAR (PPFD dir ), diffuse PAR (PPFD di f ) and scattering PAR (PPFD dir,scat ) components separately.

Direct PAR Estimation Model
In the current study, we developed a lidar-based approach to characterize the attenuation of direct solar radiation along their transmission path as penetrating through a forest canopy.As shown in the Figure 2, for a given direct light beam l with the zenith (Z) and azimuthal (A) angles reaching a given target point P (x i , y i , z i ), we could build a 3-D cylindrical buffer (S) region around the light beam l with the radius of r as it penetrates through a forest canopy.The radius r was determined by the sensitivity analysis in Section 2.4.The 3-D cylindrical buffer could be either occupied by points (i.e., "vegetation cluster") or empty (i.e., "gaps") resulting from non-random distribution of foliage elements.We projected each point within the buffer S onto the light beam l and computed their one-dimensional distances between each projected point and point P along l (d p ).The point with the longest d p was considered as the first contact point between light beams and forest canopy (P 0 ).We then applied the K-Means++ algorithm [54] to group the vegetation points and identify the gaps between foliage elements within the buffer S. The initial number of clusters (k) was set as half of the total number of points (N) within S, and the buffer was considered as an empty one if N was less than 2. Two clusters would be merged if the distance between their boundaries was less than 0.5 m determined by measuring the mean gap size manually.By iteratively clustering all points within the buffer S, we obtained five "vegetation clusters" (G 1 ~G5 ) and five "gaps" with varied sizes in the buffer S, respectively.The direct PAR values at the point P' could be obtained if the attenuation effects of direct light beams could be characterized as they penetrate through a forest canopy.By removing the gaps from the buffer S whose total path length is total PL , we computed the effective total path length ( e PL ) as: where Gi L is the length for the i-th vegetation cluster cylinder; total PL is the total path length of direct light beam starting from the first contacted point within forest canopy ( 0 P ) to the target point P' .As for the value of Gi L , it could be computed as: where ′ p ,i d is the distance of each point within the i-th vegetation cluster between its projected point and the target point P' along l, and  The direct PAR values at the point P could be obtained if the attenuation effects of direct light beams could be characterized as they penetrate through a forest canopy.By removing the gaps from the buffer S whose total path length is PL total , we computed the effective total path length (PL e ) as: where L Gi is the length for the i-th vegetation cluster cylinder; PL total is the total path length of direct light beam starting from the first contacted point within forest canopy (P 0 ) to the target point P .As for the value of L Gi , it could be computed as: where d p ,i is the distance of each point within the i-th vegetation cluster between its projected point and the target point P along l, and max(d p ,i ) and min(d p ,i ) are the maximum and minimum distances, respectively.To better characterize the non-distribution pattern of foliage elements along the light transmission path, besides the path length, we also proposed another parameter named the point density-based length (PL k ) to characterize the light attenuation effect by incorporating the point density and spatial distribution pattern within vegetation clusters.The relative volume point density (RD i ) of the i-th separated vegetation cluster G i is computed as: where AD is the average point density of tree crowns in the study area with the unit of points/m 3 , n i is the total number of points contained in G i , r is the radius of the cross-sectional circle of the cylindrical buffer S, and L Gi is the length of the vegetation cluster G i .We determined the value of AD by averaging the point volume density of 50 randomly selected spheres with a radius of 0.5 m within tree crowns.For the vegetation clusters with extremely high or low densities along all possible direct beams reaching the target point P , we plotted the histogram of all AD i values during a daily solar course and replaced abnormal (i.e., extremely large or small) AD i values using the ones at the point of 95% confidence interval.
In addition, the 3-D spatial distribution patterns of all points within each G i also affect the attenuation effects of direct solar beams.By assuming that the ability to intercept direct solar radiation was inversely related with the distances of each point to the axis of 3-D cylindrical buffer S, we proposed a parameter named "relative point distribution factor (RF i )" to characterize the spatial distribution pattern of points within vegetation clusters.In each vegetation cluster G i , the value of RF i was computed as: where d i,j is the distance of the j-th point within the i-th vegetation cluster G i to the axis of cylindrical buffer S and n i is the total number of points within G i .We obtained the average value DF based on all possible values of DF i for all direct beams reaching nine observation points during a daily solar course.AD i and DF i are the absolute volume point density and the absolute point distribution factor for the i-th vegetation cluster, respectively.Statistically, they approximately show a normal distribution with AD or DF as the mean value, so the values of RD i and RF i exhibit a positively skewed distribution with the mode of 1.When both RD i and RF i are equal to 1, we call it the standard case.A RD i value greater than 1 indicates that more foliage elements exist in the cluster G i , enhancing its attenuation effects compared to the standard case.Meanwhile, a RF i value greater than 1 indicates that more foliage elements are located closer to the light path in the cluster G i and has the same effects as RD i .Finally, we computed the value of point density-based path length (PL k ) weighted by the point volume density and point spatial distribution pattern as: In terms of the direct PAR reaching the target point P , we applied an exponent-based direct light attenuation model as: where a, b, c and d are unknown coefficients for the exponent statistical model and PL i is one of the three light path length-based parameters (PL total , PL e , PL k ).To determine the unknown parameters, we conducted the non-linear regression analysis based on the No. 5 and No. 6 pyranometer measurements in the field to fit the exponential model between PL i and direct light transmittance (PPFD dir /PPFD 0,dir ).We obtained the observed PPFD dir by subtracting the PPFD di f and PPFD dir,scat from the pyranometer readings during a daily solar course.

Diffuse PAR Estimation Model
The sky diffuse distributes ununiformly in the upper hemisphere under a clear sky condition and it was controlled by the solar position.In this study, we incorporated our diffuse sky model with an anisotropic sky diffuse distribution model [52].The model describes the variation of relative contributions of sky diffuse irradiance in each zenith-angle zone to horizontal irradiance with the various solar positions under clear sky conditions.We divided the upper hemisphere into nine zenith-angle zones with a 10 • zenith angle interval, and their relative contributions varying with the solar zenith angle are shown in Figure 3.

Diffuse PAR Estimation Model
The sky diffuse distributes ununiformly in the upper hemisphere under a clear sky condition and it was controlled by the solar position.In this study, we incorporated our diffuse sky model with an anisotropic sky diffuse distribution model [52].The model describes the variation of relative contributions of sky diffuse irradiance in each zenith-angle zone to horizontal irradiance with the various solar positions under clear sky conditions.We divided the upper hemisphere into nine zenith-angle zones with a 10° zenith angle interval, and their relative contributions varying with the solar zenith angle are shown in Figure 3.The angular gap fraction is defined as the proportion of partial or whole hemispherical sky not obscured by foliage elements of a forest canopy when viewed upwards from a specific location [55].To calculate dif PPFD , the gap fraction of each annulus zone needs to be obtained, which can be estimated as follows.Firstly, all points above the target location and within a radius of R were cut out of the original point cloud and then projected onto a plane using a stereographic projection [56].All projected points are distributed in a circle with a unit radius as shown in Figure 4. Secondly, the unit circle was divided into nine concentric annulus rings with 10° zenith angle (Figure 4).A specific zenith angle interval ( Δθ ) and azimuth angle interval ( Δφ ) were selected to divide each annulus ring into several trapezoidal facets.Due to the geometric distortion of stereographic projection, the point density of rings with higher zenith angles is much larger than those with lower zenith angles.Therefore, different values of Δθ and Δφ should be applied to each annulus ring.The Δ k θ and Δ k φ for the k-th (k=1 ~ 9) annulus ring were determined by the sensitivity analysis in Section 2.4.Thirdly, we could obtain the index ( k i , k j ) of the trapezoid that contained point ( θ, φ ) in the k-th annulus ring as: The angular gap fraction is defined as the proportion of partial or whole hemispherical sky not obscured by foliage elements of a forest canopy when viewed upwards from a specific location [55].To calculate PPFD di f , the gap fraction of each annulus zone needs to be obtained, which can be estimated as follows.Firstly, all points above the target location and within a radius of R were cut out of the original point cloud and then projected onto a plane using a stereographic projection [56].All projected points are distributed in a circle with a unit radius as shown in Figure 4. Secondly, the unit circle was divided into nine concentric annulus rings with 10 • zenith angle (Figure 4).A specific zenith angle interval (∆θ) and azimuth angle interval (∆ϕ) were selected to divide each annulus ring into several trapezoidal facets.Due to the geometric distortion of stereographic projection, the point density of rings with higher zenith angles is much larger than those with lower zenith angles.Therefore, different values of ∆θ and ∆ϕ should be applied to each annulus ring.The ∆θ k and ∆ϕ k for the k-th (k = 1 ~9) annulus ring were determined by the sensitivity analysis in Section 2.4.Thirdly, we could obtain the index (i k , j k ) of the trapezoid that contained point (θ, ϕ) in the k-th annulus ring as: where θ and ϕ are zenith angle and azimuth angle of each projected point.Considering that the area of facets near zenith is smaller than that adjacent to the horizontal plane, it is necessary to perform sinusoidal weighting on each facet [57].The weight value assigned to the trapezoid with index (i, j) in the k-th annulus ring (ω k i,j ) could be calculated as: Then, we defined l k i, j as a flag indicating whether the trapezoid with index (i, j) in the k-th annulus ring contained no projected points (empty) or not (non-empty), and the empty one was denoted as 1 while the non-empty one as 0. Finally, the angular gap fraction of the k-th annulus ring (P gap (θ k )) was calculated as the weighted ratio of the number of empty trapezoids to the total number of trapezoids in that annulus ring: = int( ) where θ and φ are zenith angle and azimuth angle of each projected point.Considering that the area of facets near zenith is smaller than that adjacent to the horizontal plane, it is necessary to perform sinusoidal weighting on each facet [57].The weight value assigned to the trapezoid with index ( i, j ) in the k-th annulus ring ( k i,j ω ) could be calculated as: Then, we defined k i,j l as a flag indicating whether the trapezoid with index ( i, j ) in the k-th annulus ring contained no projected points (empty) or not (non-empty), and the empty one was denoted as 1 while the non-empty one as 0. Finally, the angular gap fraction of the k-th annulus ring (   Therefore, PPFD di f under anisotropic clear sky diffuse distribution (PPFD di f ,ani ) was estimated as: where PPFD 0,di f is the sky diffuse PAR above the canopy, P gap (θ k ) and CR k are the gap fraction and sky diffuse proportional contribution of the k-th annulus ring, respectively.Furthermore, PPFD di f under an isotropic sky diffuse distribution (PPFD di f ,iso ) was estimated through the sky view factor (SVF): This part of PAR originates from multiple scattering of direct PAR inside the canopy, and it can be calculated as a function of the elemental scatter magnitude of PPFD 0,dir and two canopy structural parameters [58,59]: where α is the elemental scattering coefficient of broad-leaved forest leaves in the PAR band and specified as 0.15, Ω θ is the clumping index, PAI e is the plant area index of forest stands, and θ s is solar zenith angle.Ω θ describes the non-random distribution of canopy elements, which is closely related to canopy structure and tree species [10].In this study we adopted values retrieved using TLS data for the same study area [60].PAI e can be calculated using Beer's law [27,61]: where G(θ) is the mean projection of unit leaf area on the plane perpendicular to beam direction.If we divide the upper hemisphere into n annulus rings, then PAI e will be retrieved by: where θ k is the central zenith angle of the k-th annulus ring.The calculation of G(θ) requires leaf angle distribution, in this study a simple method was used [62]: where θ L is the mean leaf inclination angle.According to the result based on TLS data for the same study area [33], 66 • was used.

Model Parameters Determination
There are four parameters required to be set manually in the model.They are the radius of the spatial ray path buffer (r), the radius around the target location for cutting the point cloud (R), the zenith angle interval (∆θ k ), and the azimuth interval (∆ϕ k ) for each annulus ring.Sensitivity analysis was conducted to determine their optimal values.In terms of r, equidistant values with 2.5 cm interval ranging from 5~25cm were selected to calculate the values of PL total , PL e and PL k at two calibration observation locations during a daily solar course.Then, we fitted Equation ( 9) with these values.The buffer radius with which the model fitted best was taken as the optimal value of r.Values of 2, 3, 5 m, and every 5 m thereafter up to 70 m were selected for R. Values of 0.1 • , 0.2 • , and every 0.1 • thereafter up to 1 • along with 2 • , 3 • and 5 • were considered for ∆θ k .In addition to the above, 10  , and 30 • were also considered for ∆ϕ k .By using all combinations of R, ∆θ k and ∆ϕ k , we calculated P gap (θ k ), SVF and PAI e for each annulus ring at all observation locations [63].TLS data of nine plots were projected using stereographic projection and binarized to produce synthetic images (Figure 4b).The resolution was determined by measuring the mean distance between adjacent projected points.We used these synthetic images to calculate the angular gap fraction of each annulus ring at the corresponding observation point with in Digital Hemispherical Photography software (Natural Resources Canada, Canada Centre for Remote Sensing, Canada) as the reference value for P gap (θ k ).The combinations of R, ∆θ k and ∆ϕ k with which the root mean square error (RMSE) between calculated P gap (θ k ) and reference values reaches minimum were taken as the optimal values.
In order to further analyze how the location of annulus rings influence the optimal angular intervals, we enlarged the unit projection circle and assumed a radius of 90m.The mean area of facets formed by ∆θ k and ∆ϕ k in the k-th annulus ring was approximated by the following equation: where the unit of ∆θ k is • and the unit of ∆ϕ k is rad.

UAV Lidar Data Thinning
To investigate the effects of point densities on parameter optimization and further the estimation accuracy of the forest PAR model, we conducted a sensitivity analysis by changing the neighboring point distance (NPD) from 2.5 cm to 25 cm with an interval of 2.5 cm.The NPD was defined as the regular and minimum distance between two neighboring points of a forest lidar data [33].By doing so, we can also investigate the applicability of the proposed forest PAR model to the aerial lidar data acquisition system which usually has a higher (> 1,000 m) flight altitude and low cloud point density.

Spatiotemporal Variations of Forest PAR
The model was evaluated qualitatively by comparing modelled total PAR with the measurement at each observation location.The solar position at each logging time of measured PAR was calculated [64] and the total PAR was modelled using optimal values for the direct light attenuation model and the diffuse PAR model.Then the mean absolute error (MAE) between modelled and measured total PAR at 7 verification locations during a daily solar course was calculated.The MAE is less influenced by large outliers than RMSE and is a robust estimate of the error.Furthermore, taking the central part of the plot as an example, we mapped the spatial patterns of PAR in the forest floor at 14:00 and 16:00 on July 27, 2017, as well as 8:00, 10:00 and 12:00 on July 28, 2017 at local time.We also mapped the spatial patterns of PAR at different vertical gradients (1, 3, 5, 6, and 7 m above the forest floor) at 12:00 on July 28, 2017.In order to visually understand the PAR regime inside the canopy, the spatial patterns of total PAR in the longitudinal profiles along the solar principal plane (SPP) and its perpendicular cross-plane (PCP) were also mapped.

Determined Parameters of Direct Light Attenuation Model
In the direct PAR model, we adopted the value of 294.68 points/m 3 and 0.3342 as the average point density AD in Equation ( 5) and the average point distribution factor DF in Equation (7), respectively.Figure 5a shows the fitting results of the direct light attenuation model (Equation (9)) using direct PAR transmittance and three light path length-based parameters (PL total , PL e and PL k ) calculated from original data with various buffer radius (r).PL k showed the best result and PL e was a little worse, while PL total showed a much worse performance.All three parameters gave poor performances when a small r was used.As r increased, their performances improved significantly.After r reached 12.5 cm, the fitting residual sum of squares (RSS) and adjusted R 2 remained stable.When r was assigned 15.0 cm, the fitting RSS had a minimum of 92.47, 113.1 and 343.85 (µmol•m −2 •s −2 ) 2 while adjusted R 2 had a maximum of 0.88, 0.84 and 0.51 for PL total , PL e and PL k , respectively.Therefore, 15 cm was chosen as the optimal value of r for the original data.Figure 5b shows the scatter plot between PL k calculated with optimal r and direct PAR transmittance at two calibration locations during a daily solar course, in which direct PAR attenuates to zero rapidly with the increase of PL k .Transmittance and PL k were exponentially correlated (R 2 = 0.88, p < 0.01).Therefore, 15 cm was chosen as the optimal value of r for the original data.Figure 5b shows the scatter plot between

Determined Parameters of Sky Diffuse PAR Model
The analysis on the sensitivity of estimated SVF and e PAI to varied R shows that 35 m is the optimal plot radius for cutting the UAV-lidar data.Table 2 shows the optimal values of Δ k θ and Δ k φ for each annulus ring to calculate angular gap fraction using data with different NPDs.In general, the optimal angular intervals become larger as point density decreases and become smaller as the zenith angle of the annulus ring increases, despite some exceptions when Δ k θ and Δ k φ are considered individually.
Figure 6 shows the comparison of angular gap fraction at each observation location estimated using data with several NPDs and corresponding optimal parameters (Table 2) with reference values retrieved from TLSderived synthetic images.Obviously, nine observation locations can be divided into two categories according to the change of angular gap fraction with varied zenith angle.In terms of the first seven annular rings, with an increase of zenith angle, the angular gap fraction firstly increased and then decreased at locations 1, 3, 5, and 8, and consistently decreased at other locations.The figure manifests that data with a large range of NPDs resulting in a similar angular gap fraction, which is consistent with the reference value.∆θ1 ∆φ1 ∆θ2 ∆φ2 ∆θ3 ∆φ3 ∆θ4 ∆φ4 ∆θ5 ∆φ5 ∆θ6 ∆φ6 ∆θ7 ∆φ7 ∆θ8 ∆φ8

Determined Parameters of Sky Diffuse PAR Model
The analysis on the sensitivity of estimated SVF and PAI e to varied R shows that 35 m is the optimal plot radius for cutting the UAV-lidar data.Table 2 shows the optimal values of ∆θ k and ∆ϕ k for each annulus ring to calculate angular gap fraction using data with different NPDs.In general, the optimal angular intervals become larger as point density decreases and become smaller as the zenith angle of the annulus ring increases, despite some exceptions when ∆θ k and ∆ϕ k are considered individually.
Table 2.The optimal value of zenith angle interval (∆θ k ) and azimuth angle interval (∆ϕ k ) for calculating angular gap fraction (P gap (θ k )) of first eight annular rings using data with different NPDs.  Figure 6 shows the comparison of angular gap fraction at each observation location estimated using data with several NPDs and corresponding optimal parameters (Table 2) with reference values retrieved from TLS-derived synthetic images.Obviously, nine observation locations can be divided into two categories according to the change of angular gap fraction with varied zenith angle.In terms of the first seven annular rings, with an increase of zenith angle, the angular gap fraction firstly increased and then decreased at locations 1, 3, 5, and 8, and consistently decreased at other locations.The figure manifests that data with a large range of NPDs resulting in a similar angular gap fraction, which is consistent with the reference value.

Comparisons between the UAV-and Field-Based PARs
. Angular gap fraction at nine observation locations estimated with NPDs varying from origin to 25 cm and corresponding optimal parameters (Table 2).Reference values based on TLS-derived synthetic images are also given.Reference values for SVF and PAI e are shown in the upper left and upper right corner of each panel respectively.

Comparisons between the UAV-and Field-Based PARs
Based on the optimal value of R, ∆θ k , and ∆ϕ k , we calculated PPFD total arriving at seven verification locations at each logging moment.The model was evaluated by the mean value of MAE and bias of modelled PPFD total against measurements with various buffer radius (Figure 7a).PL k produced the highest accuracy, followed by PL e .PL total produced a much higher deviation of modelled PPFD total from the measured values.After r reached 15 cm, the means of MAE and bias approached minimums of 136.85 and −40.99 µmol•m −2 •s −1 for PL k , 164.97 and −53.044 µmol•m −2 •s −1 for PL e , 62.41 and −203.76 µmol•m −2 •s −1 for PL total , respectively.Figure 7b shows the scatter plot between modelled PPFD total using PL k with r of 15 cm and observed total PAR at seven verification locations.Modelled and observed PAR are positively correlated (R 2 = 0.97, p < 0.01).The modelled values are slightly lower than the observations.The models using PL e and PL total resulted in R 2 of 0.92 and 0.54 (scatter plots are not shown).Based on PL k and all optimal values of r, R, ∆θ k and ∆ϕ k , we estimated the dynamic trend of direct, diffuse, and total PAR at each observation location during a daily solar course and showed them in Figure 8.The model can capture the dominating variation tendency of total PAR over most of the time period but may cause some errors when the PAR fluctuates dramatically in a short period of time.

Spatiotemporal Distribution Patterns of Forest PAR
Figure 9 shows the spatial patterns of PAR in the forest floor at 14:00 and 16:00 on July 27, 2017 as well as 8:00, 10:00 and 12:00 on July 28, 2017 at local time.It is shown that not only the intensity of PAR but also the distribution patterns of shadowed and lighted ground changed obviously.In addition, Figure 10 shows the spatial patterns of PAR at different vertical gradients (1, 3, 5, 6, and 7 m above the forest floor) at 12:00 on July 28, 2017.Obviously, the shadowed area gradually decreased until it finally disappeared with the increase of the height above the ground surface, indicating the variation of the foliage distribution in the vertical direction.

Spatiotemporal Distribution Patterns of Forest PAR
Figure 9 shows the spatial patterns of PAR in the forest floor at 14:00 and 16:00 on July 27, 2017 as well as 8:00, 10:00 and 12:00 on July 28, 2017 at local time.It is shown that not only the intensity of PAR but also the distribution patterns of shadowed and lighted ground changed obviously.In addition, Figure 10 shows the spatial patterns of PAR at different vertical gradients (1, 3, 5, 6, and 7 m above the forest floor) at 12:00 on July 28, 2017.Obviously, the shadowed area gradually decreased until it finally disappeared with the increase of the height above the ground surface, indicating the variation of the foliage distribution in the vertical direction.Figure 11 shows the PAR variations in the vertical profiles along SPP and PCP at 14:00 and 16:00 on July 27, 2017 as well as 8:00, 10:00 and 12:00 on July 28, 2017.At times with relatively lower solar altitude, PAR in SPP decreased from canopy surface to crown interior along the direction of solar zenith.However, PAR in PCP showed a gradual decline from the surface of the canopy to the interior and an inconspicuous relationship with the solar zenith angle.
Figure 11 shows the PAR variations in the vertical profiles along SPP and PCP at 14:00 and 16:00 on July 27, 2017 as well as 8:00, 10:00 and 12:00 on July 28, 2017.At times with relatively lower solar altitude, PAR in SPP decreased from canopy surface to crown interior along the direction of solar zenith.However, PAR in PCP showed a gradual decline from the surface of the canopy to the interior and an inconspicuous relationship with the solar zenith angle.values for all observation locations during a daily solar course (Figure 12a) show that both the mean value and standard deviation of keep rising as r becomes larger.This indicates that a larger value of r not only enlarges parameters themselves but also introduces more severe temporal fluctuations to daytime trends of their values.The probability of each parameter in their different value ranges (Figure 12b-d) change greatly as r increases.The spatial buffer built with a small r fails to contain all points that direct light may encounter, bringing large errors to the direct light attenuation model (Figure 5).These effects can be explained in two aspects.Firstly, if r is not big enough, the two endpoints that make up each discrete vegetation clusters will probably not be the true points that separate foliage elements from gaps and will likely be located inside the canopy.This may cause insufficient clustering of foliage element points, leading to a lower estimated path length.Secondly, if r is close to or smaller than the NPD, points representing foliage elements may be misjudged as isolated noise points and removed by the model.This may cause a deficiency of vegetation clusters, resulting in underestimated parameters.However, the model accuracy will be slightly lower when r exceeds 15 cm (Figure 11.The total PAR variations in the vertical profiles along the solar principal plane (SPP) and its perpendicular cross-plane (PCP) at certain times during July 27 ~28, 2017 in the central part of the plot.

Sensitivity Analysis of 3D Ray Trace Model
Although PL k might provide the highest accuracy for the direct light attenuation model, the performance of each light path length-based parameter was determined by the spatial buffer of 3D ray trace model.How the spatial buffer radius influences each parameter can partially explain how these parameters diverge.The statistics of PL k , PL e and PL total values for all observation locations during a daily solar course (Figure 12a) show that both the mean value and standard deviation of keep rising as r becomes larger.This indicates that a larger value of r not only enlarges parameters themselves but also introduces more severe temporal fluctuations to daytime trends of their values.The probability of each parameter in their different value ranges (Figure 12b-d) change greatly as r increases.The spatial buffer built with a small r fails to contain all points that direct light may encounter, bringing large errors to the direct light attenuation model (Figure 5a).These effects can be explained in two aspects.Firstly, if r is not big enough, the two endpoints that make up each discrete vegetation clusters will probably not be the true points that separate foliage elements from gaps and will likely be located inside the canopy.This may cause insufficient clustering of foliage element points, leading to a lower estimated path length.Secondly, if r is close to or smaller than the NPD, points representing foliage elements may be misjudged as isolated noise points and removed by the model.This may cause a deficiency of vegetation clusters, resulting in underestimated parameters.However, the model accuracy will be slightly lower when r exceeds 15 cm (Figure 7a), since the spatial buffer contains many points that are far away from the solar ray and have no impact on the light attenuation.These redundant "outlying points" due to oversampling of the canopy, may overestimate the length of vegetation clusters.In models using cubic or pyramidal voxels to represent canopy structure, the level of canopy simplification and the fineness of capturing structure features depend on voxel size.Therefore, it should be properly determined firstly [2,34].Similarly, it is a primary task to select the width of the "long tube" to appropriately sample the vegetation structure and extract the part of the canopy where light attenuation mainly occurs.
Remote Sens. 2019, 11, x FOR PEER REVIEW 18 of 27 7a), since the spatial buffer contains many points that are far away from the solar ray and have no impact on the light attenuation.These redundant "outlying points" due to oversampling of the canopy, may overestimate the length vegetation clusters.In models using cubic or pyramidal voxels to represent canopy structure, the level of canopy simplification and the fineness of capturing structure features depend on voxel size.Therefore, it should be properly determined firstly [2,34].Similarly, it is a primary task to select the width of the "long tube'' to appropriately sample the vegetation structure and extract the part of the canopy where light attenuation mainly occurs.

Comparison among Three Different Path Length-Based PAR Models
In this study, both periods, and especially higher in the early morning or at dusk (Figure 13c).In some time periods when sunflecks and shadows alternate, total PL changes more sharply.This indicates that the existence of gaps between and inside tree crowns results in much instability.Canopy gaps along the ray path originate from leaves clumped around branches.If structural information for continuous canopy cannot be described at the branch level with low density ALS data, total path length and effective path length may not be discriminated clearly.Therefore, variations of direct PAR on a short time scale will be smoothed out.This could explain the result in the previous study that the total path length could capture 92% of PAR variation on a 10-min average time scale [14].However, The process

Comparison among Three Different Path Length-Based PAR Models
In this study, both PL e and PL k are feasible and powerful parameters for estimating direct light transmittance despite little difference.However, PL total fails to capture the variation of direct PAR on a five-min time scale.This result needs to be further explored.Therefore, we compared the statistics of differences among three path length-based parameters calculated with various spatial buffer radius at the same moment and showed them with each PL k value ranges (Figure 13a,b).Moreover, the true values of the three parameters during a daily solar course are also shown in Figure 13c.Overall, PL total is more than two times of PL e or PL k in most daytime periods, and especially higher in the early morning or at dusk (Figure 13c).In some time periods when sunflecks and shadows alternate, PL total changes more sharply.This indicates that the existence of gaps between and inside tree crowns results in much instability.Canopy gaps along the ray path originate from leaves clumped around branches.If structural information for continuous canopy cannot be described at the branch level with low density ALS data, total path length and effective path length may not be discriminated clearly.Therefore, variations of direct PAR on a short time scale will be smoothed out.This could explain the result in the previous study that the total path length could capture 92% of PAR variation on a 10-min average time scale [14].However, The process of light attenuation occurred almost entirely in the low value range of PL k between 0 and 2 m (Figure 5b), where the largest difference between PL k and PL e was less than 0.5m (Figure 13a).Then, as PL k became larger the absolute difference between PL k and PL e increased.The direct light completely attenuated to 0 when the absolute difference reached 1m, having no impact on the fitting of parameters in Equation (9).In contrast, PL total is higher than PL k by 2~3 m at its range of 0~2 m and the difference between PL total and PL k is 10 times higher than that between PL total and PL e (Figure 13a,b).Given the general delineation of discontinuous canopy structure on the branch level using UAV-based lidar data, we can trace the frequent movement of tiny sunflecks below canopy during the daytime.Although it is unpractical to acquire leaf number and clumping status using UAV-based lidar data, they can be implicitly expressed by the point density and distribution information.Therefore, incorporating point density and distribution factors might enhance the ability of effective path length to estimate direct light transmittance.

Sensitivity Analysis of Sky Diffuse PAR Model
In this study, the appropriate radius for cutting the point cloud (R) in sky diffuse model depends on tree species, distribution patterns, and point density.Furthermore, different optimal values of angle intervals for different annulus rings can guarantee the accuracy and reliability of angular gap fraction estimation and diffuse sky simulation for any lidar point cloud dataset.On the one hand, gap fraction analysis depends on the canopy structure within a certain radius around target locations in the canopy.The changes of SVF and e PAI with R show opposite "S" type trends (Figure 14a,b), and both of them tend to be steady after R reaches 35 m.It is insufficient for a small R to build a plot containing a substantially complete canopy structure, while a plot built with large R covers redundant points far away from the plot center.These points distribute at the edge of the unit circle.Thus, they have little impact on the distribution pattern of foliage and gaps in the hemisphere space.Most ALS-based models for retrieving canopy structural parameters need to clip the original point cloud using a

Sensitivity Analysis of Sky Diffuse PAR Model
In this study, the appropriate radius for cutting the point cloud (R) in sky diffuse model depends on tree species, distribution patterns, and point density.Furthermore, different optimal values of angle intervals for different annulus rings can guarantee the accuracy and reliability of angular gap fraction estimation and diffuse sky simulation for any lidar point cloud dataset.On the one hand, gap fraction analysis depends on the canopy structure within a certain radius around target locations in the canopy.The changes of SVF and PAI e with R show opposite "S" type trends (Figure 14a,b), and both of them tend to be steady after R reaches 35 m.It is insufficient for a small R to build a plot containing a substantially complete canopy structure, while a plot built with large R covers redundant points far away from the plot center.These points distribute at the edge of the unit circle.Thus, they have little impact on the distribution pattern of foliage and gaps in the hemisphere space.Most ALS-based models for retrieving canopy structural parameters need to clip the original point cloud using a specific radius ranging from 25 m to 60 m [65,66].Our diffuse PAR model adopts a projection strategy, but the prerequisite for a clipping radius is consistent with traditional methods based on statistical regression [67,68].On the other hand, angle intervals are also important for angular gap fraction estimation, since how the unit circle is divided dominates the probability of whether a trapezoid contains projected points.It is acceptable to ignore the uneven distribution of points in the unit projection circle derived from TLS data due to its extremely high point density and precise characterization of foliage structure.However, for UAV-based lidar data, this effect should be considered as a critical factor.Existing studies solve this problem by assigning different sizes to the corresponding pixels in synthetic images according to the position of the projected points in the unit circle [43,44], but this method cannot be applied to the case of no image generation.In this study, we divided different regions of the unit circle with different angular intervals.Figure 14c,d shows the estimation error of P gap (θ k ) for each annulus ring varying with one angle interval when the other interval is optimized and R is 35m.The differences in the trends of various annulus rings demonstrate variations of density distribution of projected points in different regions of a unit circle.The lowest value of RMSE for the outermost ring is about 2~3 times larger than that for the other rings.This can be explained by the deviation of the gap distribution between the unit circle and TLS-based synthetic image.UAV lidar fails to capture trunks and other wooden structures at the bottom of the canopy which is located in the outermost ring (Figure 4).In addition, it should be highlighted that the uncertainty may arise when this method is applied to a non-managed system where trees are not uniformly spaced.The difference between P gap (θ k ) based on ALS and TLS may increase due to more complex geometric distribution and the mutual occlusion of foliage elements.Therefore, more tests should be done when this method is applied to randomly distributed natural forests. .Our diffuse PAR model adopts a projection strategy, but the prerequisite for a clipping radius is consistent with traditional methods based on statistical regression [67,68].On the other hand, angle intervals are also important for angular gap fraction estimation, since how the unit circle is divided dominates the probability of whether a trapezoid contains projected points.It is acceptable to ignore the uneven distribution of points in the unit projection circle derived from TLS data due to its extremely high point density and precise characterization of foliage structure.However, for UAV-based lidar data, this effect should be considered as a critical factor.Existing studies solve this problem by assigning different sizes to the corresponding pixels in synthetic images according to the position of the projected points in the unit circle [43,44], but this method cannot be applied to the case of no image generation.In this study, we divided different regions of the unit circle with different angular intervals.Figure 14c,d shows the estimation error of ( ) for each annulus ring varying with one angle interval when the other interval is optimized and R is 35m.The differences in the trends of various annulus rings demonstrate variations of density distribution of projected points in different regions of a unit circle.The lowest value of RMSE for the outermost ring is about 2~3 times larger than that for the other rings.This can be explained by the deviation of the gap distribution between the unit circle and TLS-based synthetic image.UAV lidar fails to capture trunks and other wooden structures at the bottom of the canopy which is located in the outermost ring (Figure 4).In addition, it should be highlighted that the uncertainty may arise when this method is applied to a non-managed system where trees are not uniformly spaced.The difference between ( ) based on ALS and TLS may increase due to more complex geometric distribution and the mutual occlusion of foliage elements.Therefore, more tests should be done when this method is applied to randomly distributed natural forests.

Effects of Isotropic and Anisotropic Sky Diffuse Distributions
Our diffuse PAR model shows that the subcanopy sky diffuse PAR under the anisotropic sky diffuse distribution is up to 17% lower at noon and 10% higher in early morning or evening than that under the isotropic sky diffuse distribution.Figure 15a shows the normalized difference of estimated sky diffuse PAR at nine observation locations assuming anisotropic sky diffuse distribution (PPFD di f ,ani ) and isotropic sky diffuse distribution (PPFD di f ,iso ) during a daily solar course.The variation trend of the difference during a day shows symmetry at noontime with the highest solar elevation.The explanation for this phenomenon is that the region with the highest proportion of diffuse PAR in the hemispherical space switches following the solar zenith angle under anisotropic sky diffuse distribution (Figure 2), whereas PPFD di f ,iso depends only on the temporal variation of sky diffuse above canopy rather than solar position.Therefore, significant differences can be found at different moments between the distribution of diffuse PAR in each annulus ring.Taking observation location No. 5 as an example, despite low PPFD di f ,ani in the ring of 0~10 • due to obstruction by vegetation elements, PPFD di f ,ani in the top three rings still account for 50% and 67% of the total diffuse PAR at 10:00 AM and 12:00 PM, respectively (Figure 15b).It should be noted that at certain locations the normalized difference narrows near noon.This can be explained by the fact that some observation locations (such as No.5) are covered by a large proportion of foliage elements in the vicinity of the zenith direction, obstructing the penetration of a large percent of diffuse sky PAR.Considering that the sky diffuse PAR above canopy reaches its maximum value at noon, the gap distribution and canopy structure within the range of 0~30 • zenith angle above one location will dominate the magnitude of sky diffuse PAR reaching that location.

Effects of Isotropic and Anisotropic Sky Diffuse Distributions
Our diffuse PAR model shows that the subcanopy sky diffuse PAR under the anisotropic sky diffuse distribution is up to 17% lower at noon and 10% higher in early morning or evening than that under the isotropic sky diffuse distribution.Figure 15a shows the normalized difference of estimated sky diffuse PAR at nine observation locations assuming anisotropic sky diffuse distribution ( dif,ani PPFD ) and isotropic sky diffuse distribution ( dif,iso PPFD ) during a daily solar course.The variation trend of the difference during a day shows symmetry at noontime with the highest solar elevation.The explanation for this phenomenon is that the region with the highest proportion of diffuse PAR in the hemispherical space switches following the solar zenith angle under anisotropic sky diffuse distribution (Figure 2), whereas dif,iso PPFD depends only on the temporal variation of sky diffuse above canopy rather than solar position.Therefore, significant differences can be found at different moments between the distribution of diffuse PAR in each annulus ring.Taking observation location No. 5 as an example, despite low dif,ani PPFD in the ring of 0~10° due to obstruction by vegetation elements, dif,ani PPFD in the top three rings still account for 50% and 67% of the total diffuse PAR at 10:00 AM and 12:00 PM, respectively (Figure 15b).It should be noted that at certain locations the normalized difference narrows near noon.This can be explained by the fact that some observation locations (such as No.5) are covered by a large proportion of foliage elements in the vicinity of the zenith direction, obstructing the penetration of a large percent of diffuse sky PAR.Considering that the sky diffuse PAR above canopy reaches its maximum value at noon, the gap distribution and canopy structure within the range of 0~30° zenith angle above one location will dominate the magnitude of sky diffuse PAR reaching that location.

Effects of Point Density
Our lidar-based canopy PAR model showed little sensitivity to data point density when parameters were optimized.The model accuracy will only slightly deteriorate as point density decreases.Although point cloud with NPD greater than 5 cm fails to capture the detailed structure of twigs and to effectively separate the vegetation clusters inside the canopy from small gaps, it is still possible to have a comprehensive delineation of large branches and locally clumped foliage.Therefore, large gaps between adjacent canopies that dominate the radiation regime can still be effectively detected by the model.For the direct PAR model, point density determines the optimal radius of spatial ray path buffer.The correlation coefficients of fitting the direct light attenuation model with k PL illustrates that, for data with any point density from origin to NPD of 25 cm, it is possible to find an optimal r to ensure a fitting correlation coefficient above 0.85 (Figure 16a).This result indicates that the developed 3-D canopy ray trace model has great stability and good applicability to the direct light attenuation model in a large range of point density.For the sky diffuse PAR model, point density determines the optimal values of zenith and azimuth angle intervals to segment the unit circle.Assuming a 90 m radius of the unit projection circle, the variation range of k A based on data with different NPD levels is evidently larger in annulus rings with a low zenith angle (Figure 16b).This phenomenon can be explained by the uneven distribution

Effects of Point Density
Our lidar-based canopy PAR model showed little sensitivity to data point density when parameters were optimized.The model accuracy will only slightly deteriorate as point density decreases.Although point cloud with NPD greater than 5 cm fails to capture the detailed structure of twigs and to effectively separate the vegetation clusters inside the canopy from small gaps, it is still possible to have a comprehensive delineation of large branches and locally clumped foliage.Therefore, large gaps between adjacent canopies that dominate the radiation regime can still be effectively detected by the model.For the direct PAR model, point density determines the optimal radius of spatial ray path buffer.The correlation coefficients of fitting the direct light attenuation model with PL k illustrates that, for data with any point density from origin to NPD of 25 cm, it is possible to find an optimal r to ensure a fitting correlation coefficient above 0.85 (Figure 16a).This result indicates that the developed 3-D canopy ray trace model has great stability and good applicability to the direct light attenuation model in a large range of point density.For the sky diffuse PAR model, point density determines the optimal values of zenith and azimuth angle intervals to segment the unit circle.Assuming a 90 m radius of the unit projection circle, the variation range of A k based on data with different NPD levels is evidently larger in annulus rings with a low zenith angle (Figure 16b).This phenomenon can be explained by the uneven distribution characteristics of projected points in different annulus rings.The compression(expansion) effect in regions with a low(high) zenith angle enhances(weakens) the impact of point density variation.Moreover, the largest difference between A k for all annulus rings is confined to 0.2 m 2 with point density ranging from origin to NPD of 10cm.This indicates that only when point density is reduced to a certain extent, will the size of optimal angular interval alter significantly.To summarize, we recommend selecting an appropriate point density on the basis of different tree species, canopy morphology, and spatial distribution patterns in order to achieve a balance between model accuracy and efficiency.This indicates that only when point density is reduced to a certain extent, will the size of optimal angular interval alter significantly.To summarize, we recommend selecting an appropriate point density on the basis of different tree species, canopy morphology, and spatial distribution patterns in order to achieve a balance between model accuracy and efficiency.

Model Evaluation
The newly developed PAR model is superior in the following aspects.(1) The 3D ray trace model can intuitively calculate the total ray path length or effective path length at a particular moment with great accuracy.To estimate direct light transmittance, only the original lidar point cloud is needed and other raster or vector data is not required.Many spatially explicit canopy light models based on lidar data adopt a voxelization procedure to the entire canopy to track the geometric relationship between a light ray and voxels with return points contained inside [2,14,42].However, through these methods, the selection of voxel size and the voxel itself will reduce the accuracy of path length calculation to some extent.This problem can be avoided in our ray trace model.
(2) Accurate estimation of angular gap fraction for each annulus ring can be obtained, and thus can be perfectly combined with anisotropic sky diffuse distribution in the sky diffuse PAR model.TLS derived synthetic images were used as a substitution for hemispherical photographs as reference data (Figure 4).Although these two methods observe canopy structure in a similar way [69], TLS data obtained from a 3-D perspective is more stable than the fisheye photograph in extracting canopy structural parameters [70].Meanwhile, several sources of error related to acquiring and processing hemispherical photographs can be avoided.These sources include light environment and exposure settings during shooting, threshold selection for separating canopy and sky pixels, accurate georeference of observation locations, and height deviation between the camera and observation locations [69,71].(3) Our model has expanded traditional methods based on hemispherical photographs to largescale applications with great flexibility.The temporal variation of PAR at any location is easily available, which can support ecological models focusing on soil temperature, water content and snowmelt dynamics [19,72].Since the canopy PAR model simulates PAR pointwise, the sampling density can be flexibly selected according to the size of study area and the degree of horizontal heterogeneity to improve processing efficiency.
There are several sources of error in simulating PAR.The primary one is derived from the occlusion effect of laser beams.Although UAV can easily densify data points by overlapping flight paths in different directions, the occlusion effect in the dense canopy is still inevitable [48,73].It can be effectively eliminated by changing the flight path direction or increasing the scan angle, but the influence of scan angles on the inversion of canopy

Model Evaluation
The newly developed PAR model is superior in the following aspects.(1) The 3D ray trace model can intuitively calculate the total ray path length or effective path length at a particular moment with great accuracy.To estimate direct light transmittance, only the original lidar point cloud is needed and other raster or vector data is not required.Many spatially explicit canopy light models based on lidar data adopt a voxelization procedure to the entire canopy to track the geometric relationship between a light ray and voxels with return points contained inside [2,14,42].However, through these methods, the selection of voxel size and the voxel itself will reduce the accuracy of path length calculation to some extent.This problem can be avoided in our ray trace model.(2) Accurate estimation of angular gap fraction for each annulus ring can be obtained, and thus can be perfectly combined with anisotropic sky diffuse distribution in the sky diffuse PAR model.TLS derived synthetic images were used as a substitution for hemispherical photographs as reference data (Figure 4).Although these two methods observe canopy structure in a similar way [69], TLS data obtained from a 3-D perspective is more stable than the fisheye photograph in extracting canopy structural parameters [70].Meanwhile, several sources of error related to acquiring and processing hemispherical photographs can be avoided.These sources include light environment and exposure settings during shooting, threshold selection for separating canopy and sky pixels, accurate georeference of observation locations, and height deviation between the camera and observation locations [69,71].(3) Our model has expanded traditional methods based on hemispherical photographs to large-scale applications with great flexibility.The temporal variation of PAR at any location is easily available, which can support ecological models focusing on soil temperature, water content and snowmelt dynamics [19,72].Since the canopy PAR model simulates PAR pointwise, the sampling density can be flexibly selected according to the size of study area and the degree of horizontal heterogeneity to improve processing efficiency.
There are several sources of error in simulating PAR.The primary one is derived from the occlusion effect of laser beams.Although UAV can easily densify data points by overlapping flight paths in different directions, the occlusion effect in the dense canopy is still inevitable [48,73].It can be effectively eliminated by changing the flight path direction or increasing the scan angle, but the influence of scan angles on the inversion of canopy structural parameters should also be considered [34].UAV should maintain stable speed and uniformity of flight path during scanning to avoid uneven distribution of point density.The difference in data acquisition methods is the posterior source of model error.UAV-based lidar provides a relatively complete depiction of the upper part of the canopy structure, while it fails to obtain a comprehensive description to the bottom part, especially for the woody components.On the contrary, TLS or a fisheye camera scans the upper hemisphere space based on a location below the canopy.They can capture the precise understory structural information but lack the morphological description for the upper part of the canopy [44,48].When using hemispherical photographs or TLS-based synthetic images as reference, many woody components such as trunks, low branches, and shrubs are included, which will influence the estimation of canopy structural parameters [73].As shown in Figure 4, after point projection, this difference is relatively limited in the region with a low zenith angle and more pronounced at regions with a high zenith angle.This error can be limited by removing the outermost ring or increasing point density.In addition, transforming the 3-D point cloud to a 2-D plane by projection will result in some loss of partial canopy structure information.Subjectivity remains when selecting appropriate binary image resolution after TLS point cloud projection and may bring some error to sky diffuse PAR estimation.
Hemispherical photographs or TLS derived synthetic images is necessary in the sky diffuse PAR model as reference data for angular gap fraction estimation.TLS data used in our study has several advantages mentioned above.Nevertheless, hemispherical photography is faster, cheaper, and easier to process, especially when the spatial coverage is huge.We recommend that a reasonable choice should be made according to the coverage of study area and cost.The direct PAR and sky diffuse PAR reaching an arbitrary location in the canopy are dominated by return points around that location and the ray path respectively.Under the premise that the point cloud can describe canopy structure substantially and completely, the modelled subcanopy PAR is determined by the optimal model parameters selected manually.Values of these parameter are not universal or constant but need adjusting according to tree species, canopy distribution patterns, and data point density.Therefore, this model is supposed to be further tested in more comprehensive conditions.Moreover, topography should be considered in this model for further direction.

Conclusions
We proposed a novel spatially-explicit canopy PAR model according to UAV-based lidar data to simulate spatiotemporal distributions of PAR within or under a discontinuous broad-leaved forest canopy.The study compared direct PAR simulated by three light path length-based parameters and sky diffuse PAR simulated under conditions of isotropic and anisotropic sky diffuse distribution.The effects of data point density on model parameterization and accuracy were also discussed.For a discontinuous canopy, the effective path length of light passing through foliage elements is a feasible and powerful parameter to capture the spatiotemporal variations of total PAR along a light transmission path.In addition, incorporating point density and spatial distribution factors will further improve the accuracy of final estimation.However, the total path length fails to capture direct light variation on a five-minute time scale.In the meantime, the estimated sky diffuse PAR is up to 17% higher at noon and 10% lower at sunrise or sunset by assuming the anisotropic sky diffuse distribution rather than the isotropic sky diffuse distribution.The point density has a limited impact on model accuracy as long as parameters are determined as their optimal values.The developed canopy PAR model

Figure 1 .
Figure 1.The illustration for the study sites.(a) Location of the study site.(b) 3-D point cloud of the plot.(c) Distribution of calibration and verification locations with canopy height model (CHM) shown below.

2. 2 .Figure 1 .
Figure 1.The illustration for the study sites.(a) Location of the study site.(b) 3-D point cloud of the plot.(c) Distribution of calibration and verification locations with canopy height model (CHM) shown below.

Figure 2 . 1 GP 5 G 5 G ( 5 RF
Figure 2. The 3D canopy ray-tracing model.(a) The straight dashed line l indicates a direct sunbeam reaches the target point P' within a 3-D cylindrical buffer S after penetrating through a forest canopy.(b) Details in the spatial ray path buffer S. According to the distance of each point in S along l to the point P' ( ′ p d ), S can be separated into maximum and minimum distances, respectively.To better characterize the non-distribution pattern of foliage elements along the light transmission path, besides the path length, we also proposed another parameter named the point density-based length ( k PL ) to characterize the light attenuation effect by incorporating the point density and spatial distribution pattern within vegetation clusters.The relative volume point density ( i RD ) of the i-th separated vegetation cluster i G is computed as:

Figure 2 .
Figure 2. The 3D canopy ray-tracing model.(a) The straight dashed line l indicates a direct sunbeam reaches the target point P within a 3-D cylindrical buffer S after penetrating through a forest canopy.(b) Details in the spatial ray path buffer S. According to the distance of each point in S along l to the point P (d p ), S can be separated into five vegetation clusters (G 1 ~G5 ) and five gaps.The effective path length of light passing through the canopy (PL e ) is the sum of L G1 ~LG5 , and the total path length (PL total ) is the distance from P to P 0 along l.The cross-sectional view of point distribution in vegetation clusters G 5 along the ray direction is displayed in the right.The relative point distribution coefficient of G 5 (RF 5 ) is calculated by the point-to-ray distance (d 5, j ) of each point within G 5 .

i
where a , b, c and d are unknown coefficients for the exponent statistical model and i PL is one of the three light path length-based parameters ( total PL , To determine the unknown parameters, we conducted the non-linear regression analysis based on the No. 5 and No. 6 pyranometer measurements in the field to fit the exponential model between i PL and direct light transmittance ( readings during a daily solar course.

Figure 3 .
Figure 3.The proportional contributions of sky radiance in each zenith-angle zone to the diffuse irradiance on a horizontal surface varied with solar zenith angle under a clear sky condition.

Figure 3 .
Figure 3.The proportional contributions of sky radiance in each zenith-angle zone to the diffuse irradiance on a horizontal surface varied with solar zenith angle under a clear sky condition.
θ ) was calculated as the weighted ratio of the number of empty trapezoids to the total number of trapezoids in that annulus ring:

Figure 4 .
Figure 4.The stereographic projection for lidar data.(a) The projected points of UAV lidar point cloud centered at the No. 5 observation location with a radius of 30 m.(b) The synthetic image derived from the TLS point cloud at the same location.The circle with a unit radius and the synthetic image were divided into 9 annulus rings for comparison.

Figure 4 .
Figure 4.The stereographic projection for lidar data.(a) The projected points of UAV lidar point cloud centered at the No. 5 observation location with a radius of 30 m.(b) The synthetic image derived from the TLS point cloud at the same location.The circle with a unit radius and the synthetic image were divided into 9 annulus rings for comparison.

2 R
Sens. 2019, 11, x FOR PEER REVIEW 12 of 27 improved significantly.After r reached 12.5 cm, the fitting residual sum of squares (RSS) and adjusted remained stable.When r was assigned 15.0 cm, the fitting RSS had a minimum of 92.47, 113.1 and 343.85 (μmol•m −2 •s −2 ) 2 while adjusted R 2 had a maximum of 0.88, 0.84 and 0.51 for total PL , e PL and k PL , respectively.

k
PL calculated with optimal r and direct PAR transmittance at two calibration locations during a daily solar course, in which direct PAR attenuates to zero rapidly with the increase of k PL .Transmittance and k PL were exponentially correlated (R 2 = 0.88, p < 0.01).

Figure 5 . 2 R
Figure 5. Fitting results of direct light attenuation model with three path length-based parameters calculated at two calibration locations.(a) The residual sum of squares and adjusted 2 R of fitting the direct light attenuation model

Figure 5 .
Figure 5. Fitting results of direct light attenuation model with three path length-based parameters calculated at two calibration locations.(a) The residual sum of squares and adjusted R 2 of fitting the direct light attenuation model with PL k , PL e and PL total calculated with various r.(b) Scatter plot of direct PAR transmittance and PL k calculated with optimal r.The red line indicates the fitting equation.

Figure 6 .
Figure 6.Angular gap fraction at nine observation locations estimated with NPDs varying from origin to 25 cm and corresponding optimal parameters (Table 2).Reference values based on TLS-derived synthetic images are also given.Reference values for SVF and PAI e are shown in the upper left and upper right corner of each panel respectively.

2 R
verification locations.Modelled and observed PAR are positively correlated ( R = 0.97, p < 0.01).The modelled values are slightly lower than the observations.The models using of 0.92 and 0.54 (scatter plots are not shown).Based on k PL and all optimal values of r, R, Δ k θ and Δ k φ , we estimated the dynamic trend of direct, diffuse, and total PAR at each observation location during a daily solar course and showed them in Figure8.The model can capture the dominating variation tendency of total PAR over most of the time period but may cause some errors when the PAR fluctuates dramatically in a short period of time.

Figure 7 .
Figure 7. Evaluation results of the canopy radiation model using three light path length-based parameters at seven verification locations.(a) Variations of the mean value of MAE and bias between modelled and measured PAR with varied spatial buffer radius.(b) Scatter plot of modelled (y) total PAR based on

Figure 7 . 27 Figure 8 .
Figure 7. Evaluation results of the canopy radiation model using three light path length-based parameters at seven verification locations.(a) Variations of the mean value of MAE and bias between modelled and measured PAR with varied spatial buffer radius.(b) Scatter plot of modelled (y) total PAR based on PL k with 15cm of r and observed (x) total PAR.The red line indicates the fitting equation and thin-dotted line indicates the 1:1 line.Remote Sens. 2019, 11, x FOR PEER REVIEW 15 of 27

Figure 8 .
Figure 8. Variation trends of estimated direct, diffuse, total PAR and observed total PAR at each observation locations during one daily solar course.

Figure 9 .
Figure 9. Spatial patterns of total PAR in the forest floor at certain times during July 27~28, 2017 in the central part of the plot.

Figure 10 .
Figure 10.Spatial patterns of total PAR along the vertical gradients (1, 3, 5, 6, and 7 m above the forest floor) at 12:00 on July 28, 2017 in the central part of the plot.

Figure 9 . 27 Figure 9 .
Figure 9. Spatial patterns of total PAR in the forest floor at certain times during July 27~28, 2017 in the central part of the plot.

Figure 10 .
Figure 10.Spatial patterns of total PAR along the vertical gradients (1, 3, 5, 6, and 7 m above the forest floor) at 12:00 on July 28, 2017 in the central part of the plot.

Figure 10 .
Figure 10.Spatial patterns of total PAR along the vertical gradients (1, 3, 5, 6, and 7 m above the forest floor) at 12:00 on July 28, 2017 in the central part of the plot.

Figure 11 .
Figure 11.The total PAR variations in the vertical profiles along the solar principal plane (SPP) and its perpendicular cross-plane (PCP) at certain times during July 27 ~ 28, 2017 in the central part of the plot.
Although k PL might provide the highest accuracy for the direct light attenuation model, the performance of each light path length-based parameter was determined by the spatial buffer of 3D ray trace model.How the spatial buffer radius influences each parameter can partially explain how these parameters diverge.The statistics of

Figure 12 .
Figure 12.Sensitivity analysis of three attenuation parameters for all observation locations at all recording times to the radius of spatial ray path buffer (r, in the 3D ray trace model).(a) Mean value and standard deviation of k PL , and powerful parameters for estimating direct light transmittance despite little difference.However, total PL fails to capture the variation of direct PAR on a five-min time scale.This result needs to be further explored.Therefore, we compared the statistics of differences among three path length-based parameters calculated with various spatial buffer radius at the same moment and showed them with each k PL value ranges (Figure 13a,b).Moreover, the true values of the three parameters during a daily solar course are also shown in Figure 13c.Overall,

Figure 12 .
Figure 12.Sensitivity analysis of three attenuation parameters for all observation locations at all recording times to the radius of spatial ray path buffer (r, in the 3D ray trace model).(a) Mean value and standard deviation of PL k , PL e and PL total with varied r. (b)-(d) Percentage of PL k , PL e and PL total values within their different value ranges with varied r.The solid square represents the zero value and the solid circle, triangle and diamond represent each value range.
Remote Sens. 2019, 11, x FOR PEER REVIEW 19 of 27 of light attenuation occurred almost entirely in the low value range of k PL between 0 and 2 m (Figure 6b), where the largest difference between k PL and e PL was less than 0.5m (Figure 13a).Then, as k PL became larger the absolute difference between k PL and e PL increased.The direct light completely attenuated to 0 when the absolute difference reached 1m, having no impact on the fitting of parameters in Equation (9).In contrast, total PL is higher than k PL by 2~3 m at its range of 0~2 m and the difference between 13a,b).Given the general delineation of discontinuous canopy structure on the branch level using UAV-based lidar data, we can trace the frequent movement of tiny sunflecks below canopy during the daytime.Although it is unpractical to acquire leaf number and clumping status using UAV-based lidar data, they can be implicitly expressed by the point density and distribution information.Therefore, incorporating point density and distribution factors might enhance the ability of effective path length to estimate direct light transmittance.

Figure 13 .
Figure 13.Path length-based parameters comparison: (a) and (b) changes of means and standard deviations of the difference between

Figure 13 .
Figure 13.Path length-based parameters comparison: (a) and (b) changes of means and standard deviations of the difference between PL k and PL e or PL total calculated with various spatial buffer radius (represented by different solid symbols) at exact same moments for all observation locations with PL k .(c) variations of three path length-based parameters during a daily solar course at the No. 5 observation location with a spatial buffer radius of 15cm.

Figure 14 .
Figure 14.Sensitivity analysis of parameters in the sky diffuse PAR model: (a) and (b) Changes of means and standard deviations of SVF and

Figure 14 .
Figure 14.Sensitivity analysis of parameters in the sky diffuse PAR model: (a) and (b) Changes of means and standard deviations of SVF and PAI e at 9 observation locations with R (original data with the optimal angle intervals given in Table 2 were used).(c) and (d) show changes of RMSE of estimated P gap (θ k ) against reference values for each annulus ring with zenith and azimuth angle intervals, respectively.For a given zenith/azimuth angle interval, azimuth/zenith angle interval was optimized with R = 35 m.

Figure 15 .
Figure 15.Comparison between isotropic and anisotropic sky diffuse distributions: (a) The normalized difference of

Figure 15 .
Figure 15.Comparison between isotropic and anisotropic sky diffuse distributions: (a) The normalized difference of PPFD di f ,ani and PPFD di f ,iso during a daily solar course.The red line and shading area indicate the mean value and the range of normalized difference at nine observation locations respectively.(b) The distribution of sky diffuse PAR arriving at the No. 5 observation location in each annulus ring at different moments under anisotropic sky diffuse distribution.
Remote Sens. 2019, 11, x FOR PEER REVIEW 22 of 27 characteristics of projected points in different annulus rings.The compression(expansion) effect in regions with a low(high) zenith angle enhances(weakens) the impact of point density variation.Moreover, the largest difference between k A for all annulus rings is confined to 0.2 m 2 with point density ranging from origin to NPD of 10cm.

Figure 16 .
Figure 16.The relationship between point density and model parameters.(a) The correlation coefficients of fitting the direct light attenuation model with

Figure 16 .
Figure 16.The relationship between point density and model parameters.(a) The correlation coefficients of fitting the direct light attenuation model with PL k calculated from data with various NPD levels varied with spatial buffer radius.(b) The distribution of A k in each annulus ring based on data with different NPD levels.

Table 1 .
Locations of nine observation points in the experimental forest plot and calibration coefficients for each pyranometer at the corresponding location.

Table 2 .
The optimal value of zenith angle interval ( Δ