Airborne laser scanning proxies of canopy light transmission in forests

Reliable estimates of canopy light transmission are critical to understanding the structure and function of vegetation communities but are difficult and costly to attain by traditional field inventory methods. Airborne laser scanning (ALS) data uniquely provide multi-angular vertically resolved representation of canopy geometry across large geographic areas. While previous studies have proposed ALS indices of canopy light transmission, new algorithms based on theoretical advancements may improve existing models. Herein, we propose two new models of canopy light transmission (i.e., gap fraction, or Po, the inverse of angular canopy closure). We demonstrate the models against a suite of existing models and ancillary metrics, validated against convex spherical densiometer measurements for 950 field plots in the foothills of Alberta, Canada. We also tested the effects of synthetic hemispherical lens models on the performance of the proposed hemispherical Voronoi gap fraction (Phv) index. While vertical canopy cover metrics showed the best overall fit to field measurements, one new metric, point-densitynormalized gap fraction (Ppdn), outperformed all other gap fraction metrics by two-fold. We provide suggestions for further algorithm enhancements based on validation data improvements. We argue that traditional field measurements are no longer appropriate for ‘ground-truthing’ modern LiDAR or SfM point cloud models, as the latter provide orders of magnitude greater sampling and coverage. We discuss the implications of this finding for LiDAR applications in forestry. 2 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40


Introduction
The light environment is a critical factor for the structure and function of vegetation communities (Dengel and Grace, 2010;Gamon, 2014;Gamon and Bond, 2013;Saeki, 2005, 1953). In northern forests, tree crown geometries are well suited to a low solar elevation, occluding less light from neighboring trees (Aakala et al., 2016). Understory light is an important factor in the successional trajectory of forests through vegetation establishment and growth, making it a critical parameter required to forecast forest ecosystems . Although understory light is a function of quantifiable variation in local stand geometry, topographic position, atmospheric conditions, and solar position, it remains difficult and costly to measure. While the importance of understory light has long been understood (Monsi and Saeki, 1953), it is notoriously difficult to measure with remote sensing methods. The advent of multiangular remote sensing technologies such as airborne laser scanning LiDAR (ALS) and photogrammetric computer vision have made it possible to map canopy light transmission as a proxy for understory light by assuming beam canopy penetration equivalent to a Poisson process. Monsi & Saeki (1953) were the first to represent contact frequency as a Poisson process, equivalent to the Beer-Lambert law (Hancock, 2010).
Due to limitations in spaceborne sensor resolution and coverage, given the large footprint and limited coverage (Coyle et al., 2015;Dubayah et al., 2014). Data assimilation or imputation techniques are required to generate wall-to-wall maps from sparse spaceborne LiDAR data, increasing the uncertainty of estimates through the inclusion of an additional model. Despite recent advances in deriving forest canopy geometry from commercial passive optical spaceborne sensors (Shean et al., 2016), active optical airborne LiDAR systems remain ideal instruments for estimating understory light conditions at the landscape scale. This is due to the high precision (i.e., point density and geolocation error), broad spatial coverage, and availability of data in many countries, allowing direct measurement of canopy light transmission with multi-angular pulses of near-infrared photons and multi-return or waveform detectors (i.e., photodiodes).
Airborne laser scanning (ALS) is used throughout boreal forests and contains detailed information on forest geometry at scales ranging from stands to landscapes .
Many of these ALS metrics may be readily applied as indices of canopy light transmission, individually or in combination. Some of the earliest, simplest, and most effective metrics of ACC and thus P o are based on the ratio of ground-to-canopy returns (Korhonen et al., 2011;Morsdorf et al., 2006;Riaño et al., 2004;Solberg et al., 2009). The metric of Solberg et al. (2009) differs in that it corrects for pulses that have returns from both the canopy and ground, assigning a partial cover value to these. A pulse intensity-based approach was designed to correct for two-way transmission loss (Hopkinson and Chasmer, 2007), also novel for utilizing target reflectance information. More recent approaches provide hemispherically projected LiDAR metrics comparable to traditional ground measurements (Parent and Volin, 2014;Varhola et al., 2012), while others further utilize geometric operations to improve the estimation of cover (Alexander et al., 2013). An opportunity exists to improve simple transmission metrics and advanced representations of forest geometry to estimate cover, as the theory surrounding both continue to improve. While future studies should apply deep neural networks designed for scattered, unordered point data, such as using models based on the PointNet++ architecture (Qi 5 et al., 2017), we focus on simple geometric operations for their diminished need for labeled data and speed/ease of computation in large-area mapping applications.
Calculations of forest structural parameters from ALS are often distinct from those of traditional ground methods, due to differences in sampling bias (top-vs. bottom-of-canopy), lending to variation in terminology and methodology. Canopy light attenuation calculations based on ALS often assume canopy light transmission (T) equivalent to canopy gap fraction (P o ), each inverses of vertical canopy cover (VCC) and angular canopy closure (ACC), as provided in the following equation (Gonsamo et al., 2013;Hopkinson and Chasmer, 2009;Morsdorf et al., 2006): Traditionally, VCC quantifies the 2-D areal canopy coverage, while T is a function of incident photosynthetically active radiation (PAR), fraction of absorbed PAR (fPAR) by leaf absorptance, leaf transmissivity, and scattering, incorporating leaf chemistry, geometry, position, and orientation effects on the bidirectional reflectance distribution function, or BRDF (Gastellu-Etchegorry et al., 1996). While the equivalence of T and P o holds in the absence of detailed information, the two metrics remain distinct, providing different -though complementaryinformation (Gonsamo et al., 2013).
Although ALS pulses are typically emitted at narrow zenith angles less than 20 degrees from nadir, they provide an empirical test of angular light penetration through the canopy, making ALS suitable for estimating P o . Meanwhile, VCC is often calculated from ALS for each 6 cell using narrow incoming zenith angles between 0 and 10, opposite to scan and beam divergence source angle (Morsdorf et al., 2006;Weiss et al., 2004). Hence, the measurement of VCC with ALS is often a field-of-view, or scope, function (Lee et al., 2008), rather than a true measure of 2-D areal coverage (although simple grid-based methods exist), making it sensitive to neighborhood effects. Here, as with leaf area index (L), gridded ALS-derived metrics (e.g., the ratio of canopy first-returns to ground first-returns) are more compatible with the classical definition of VCC. Similar challenges of sampling bias have been reported for gap fraction (P o ) estimates derived from terrestrial laser scanning (TLS) LiDAR (Vaccari et al., 2013).
The objective of this study was to develop new ALS metrics and regression models of T that can be extended to forest landscape models to simulate understory irradiation across large areas. Four new ALS metrics for retrieving T are presented, including hemispherical Voronoi gap fraction (P hv ), point-density normalized gap fraction (P pdn ), and their inverses, hemispherical Voronoi angular canopy closure (ACC hv ) and point-density normalized angular canopy closure (ACC pdn ). While P hv and ACC hv are intended to improve estimates of canopy light interception from LiDAR with varying sensor properties, P pdn and ACC pdn are intended to reduce sensor effects by normalizing hemispherical sectors by their surface area and the overall point density.
The four new hemispherical canopy metrics (P hv , P pdn , ACC hv , and ACC pdn ), nine vertical canopy cover (VCC) metrics, twelve stem and crown metrics, and five other metrics, for a total of 30 metrics (Table 1), were validated against traditional coarse-resolution convex spherical densiometer ground measurements of angular canopy closure (ACC), representing the inverse of T. The P hv metric was applied using four different hemispherical lens geometries at canopy height thresholds varying from one meter to five meters in 0.25 m steps, for a total of 68 different P hv configurations for each plot. In doing so, we provide key innovations that are readily deployable across a range of forested systems with available ALS data, as the P pdn method is designed to overcome common shortcomings related to changes in LiDAR sensor design over time. Thus, it is anticipated that the P pdn method may be highly valued by the forestry industry for operational use. Furthermore, we provide a future direction for research along both detailed geometric and generalized coefficient approaches. Finally, we make all of our innovations openly available for use in the gapfraction package for R

Materials and methods
Vegetation ground plot measurements were collected in the Hinton Forest Management Area in the early 2000s during summer (leaf-on) conditions. While details of the area have been documented in previous research (Nielsen, 2005;Nielsen et al., 2006Nielsen et al., , 2004, the foothills region is generally characterized by monospecific stands of lodgepole pine (Pinus contorta Douglas ex Louden), well-drained post-glacial soils, moderate temperatures and precipitation, and extensive forest management (Natural Regions Committee, 2006). Angular canopy closure (ACC), and thus canopy gap fraction (P o = 1-ACC), was measured from breast-height using a convex spherical densiometer. Densiometer measurements were recorded for each of the four cardinal directions and averaged for each plot (Lemon, 1956;Nielsen, 2005). the Optech ALTM 3100 is one of the first commercial ALS systems capable of full-waveform digitization, the system used in this study is a discrete-return system, recording up to four returns for every laser pulse, each with 12-bit dynamic range intensity information (Hilker et al., 2013).
Ground and non-ground returns were classified using Terrasolid TerraScan version 0.6 consumer-off-the-shelf (COTS) software, which applies previously demonstrated methods (Kraus and Pfeifer, 1998

Data pre-processing
Using LAStools (Isenburg, 2015), the ALS tiles were height-normalized before extracting circular field plots with a 50 m radius, based on previous research exhibiting a saturation of edge effects below this radius threshold (Alexander et al., 2013;Zhao and Popescu, 2009). Delaunay triangle-position elevation from each return's z value. LAStools implements an optimized variant of the best available ground plane extraction algorithm (Axelsson, 1999;Maguya et al., 2014), modified to include Delaunay streaming or triangulated irregular network (TIN) streaming (Isenburg et al., 2006b(Isenburg et al., , 2006a(Isenburg et al., , 2006c for improved computational efficiency on large datasets. Maximum point height was filtered at 40 m, based on local tree species ground measurements. The ALS plots were processed with a series of point cloud metrics implemented in custom R scripts (R Core Team, 2015), described below. Finally, the top performing ALS metric (VCC fci ) was applied to an expanded set of ALS plots to analyze variation related to species composition and age class.

Spike-free canopy height model algorithm
The first step required the generation of continuous canopy height models (CHMs) without smoothing-or sampling-related artifacts. This was due to pitting in the simple gridded maxima CHMs given a mean point density below 2 points m -2 , known to affect the accuracy of tree detection. In order to improve CHM inputs for individual tree crown (ITC) detection, a layered 2-D adaptation of the spike-free CHM algorithm (Khosravipour et al., 2016(Khosravipour et al., , 2014  Equivalent in output to the original, our modified implementation of the spike-free CHM algorithm vertically stratifies all returns into user-defined windows or slices to constrain 12 Delaunay triangulations, which can be absolute distances or height percentiles. A 2 m height threshold was used with steps at 5, 10, and 15 m, as in the pit-free CHM work (Khosravipour et al., 2014). Delaunay triangles with edge lengths exceeding a user-defined threshold are filtered to limit smoothing, set to the default value (Khosravipour et al., 2014). The final CHM consists of continuous height maxima along raster grid points. This adaptation takes advantage of vertical stratification to generate non-overlapping points necessary for 2-D Delaunay triangulation. The theoretical advantage over the 3-D Constrained Delaunay approach (Khosravipour et al., 2016) is chiefly computational for the sake of speed and simplicity. These and other functions are provided in the gapfraction package for R (https://adamerickson.xyz/gapfraction/).

Hemispherical Voronoi gap fraction
The hemispherical Voronoi gap fraction (P hv ) index represents P o as the areal coverage of Voronoi tessellation cells above a given canopy height threshold from the perspective of standing at the plot center and looking toward the zenith, identical to a traditional hemispherical photograph. The plot center at 3-D local Cartesian coordinate (x=0, y=0, z=0) is set equal to the hemispherical camera model principal point, or intersection of the optical axis and image plane.
The ground plane is set equal to the image plane, with the optical axis pointing skyward at the zenith. Once the LiDAR data is pre-processed into normalized heights and local Cartesian coordinates, the first step is to re-project the LiDAR points into image coordinates based on a model of a fisheye (hemispherical) lens.

The projection of a 3-D point
requires a mathematical model of a fisheye lens, consisting of a series of transformations with extrinsic and intrinsic camera parameters (Abraham and Förstner, 2005;Ray, 2002). The extrinsic parameters map the real-world coordinates into camera coordinates, while the intrinsic parameters map the camera coordinates onto the image plane. The image coordinate calculations take the following form (Abraham and Förstner, 2005): Here, c x and c y are the principal distances (this allows for non-square pixels), φ and θ are the azimuthal and polar angles, respectively, r*(θ) is the radial projection function, or mapping function, and, x' H and y' H are the coordinates of the principal point, or the intersection of the optical axis and the image plane. The distortion model parameters used for real-world lenses, Δx' and Δy', typically added to the end of their corresponding equations, are omitted. To change to a different hemispherical camera model, the radial projection function can be simply modified.
The classical pinhole camera is described by the perspective projection function of the form r' = c tan(θ), where r' is the radial distance from the principal point on the image plane and c is the principal distance, a function of the focal length and focal distance (Fourcade, 1928). Here, all four projections are implemented with a circular design in the gapfraction package for R. The radial projection function, or mapping function, for each projection is as follows (Abraham and Förstner, 2005;Ray, 2002): To transform the real-world coordinates to camera coordinates, the normalized point clouds were projected into 3-D local Cartesian coordinates with an (x, y, z) tuple centroid of (0, 0, 0). A function was developed that allows this calculation without plot center geolocation information to ease LiDAR plot processing. The function sets the midpoint of the vector of X and Y values to half of the range, as shown below: To transform the camera coordinates into image plane coordinates, the 3-D local Cartesian coordinates are projected into 2-D polar coordinates (azimuth angle and radial distance, or φ and r) before projecting the 2-D polar coordinates into 2-D Cartesian space with standard trigonometric equations, where x' = r cos(φ) and y' = r sin(φ). The calculations were implemented in their normalized image plane form (Abraham and Förstner, 2005), as the 3-D local Cartesian coordinates were normalized to their true distance values in meters, rather than the typical unit sphere. This was done to preserve 3-D Cartesian distances for calculations that do not require hemispherical or image plane projections.
Once the LiDAR data were projected onto the 2-D hemispherical image plane, the 2-D Delaunay triangulation and Voronoi tessellation were computed for the planar point sets using the deldir package for R (Turner, 2015), filtering points below a user-defined canopy threshold.
The summed area of filtered cells, or gaps, was calculated as a percentage of the overall plot area, providing the hemispherical Voronoi gap fraction (P hv ). This assumes 100% light occlusion by non-filtered cells. The implication of this simplification is that light attenuation is The φ values were rescaled from (-π, π) to the interval (0, 2π) by adding 2π to φ values where φ is less than zero. Based on previous research (Zhao & Popescu, 2009) coordinates were sectioned at polar and azimuthal increments of 5° and 45°, respectively, producing 18 x 8 sky sectors for a total of 144 sectors. A polar resolution of 15° is also commonly used in LiDAR studies (Korhonen and Morsdorf, 2014), but is likely coarser than necessary for modern sensors. The number of first returns per hemispherical sector was calculated using the following equation: is the number of elements contained in a set defined by the intersection of polar and azimuthal angle subsets, θ returns i and φ returns j , at hemisphere sector intervals defined by steps i and j, respectively. A matrix is produced containing the frequency of returns within each sector of the hemisphere. In order to account for varying sector sizes, the values are adjusted by the hemispherical surface area of each sector. To do so, the surface area of each hemispherical sector is first calculated, as follows: This produces a second matrix of equal dimensions, i x j. Here, A i,j is the area of a sector for polar angle ϴ i and azimuth angle φ j at intervals defined by steps i and j, while R is the radius This mitigates issues related to sensor effects (e.g., point density). The filtering of non-firstreturns is necessary to also reduce sensor effects along the z-axis, as vertical resolution can vary due to a number of factors. Point-density normalized canopy gap fraction (P pdn ) was calculated with the following equation:

Comparison with other ALS metrics
A set of standard metrics were also implemented to assess their performance against new methods and ground measurements. The method comparison framework includes estimates of canopy gap fraction, angular canopy closure, vertical canopy cover, individual tree detection, crown area, distance to crown and canopy, leaf area index, and clumping. First, these methods are described in the following paragraphs.
Based on previous research on the estimation of leaf area index (Lang and Yueqin, 1986;Miller, 1967;Ryu et al., 2010;Zhao and Popescu, 2009), the effective leaf area index (L e ) was calculated using the following equation (Korhonen and Morsdorf, 2014): The apparent clumping index (Ω app ) was calculated based on a ratio of two L e estimation methods (Ryu et al., 2010). The previous approach was modified by approximating the integral as a summation, with each L e method weighted by the sine of the given polar angle, θ (Korhonen and Morsdorf, 2014): Next, the L e vector is used for n polar angles θ to calculate the canopy gap fraction per the Beer-Lambert Law Saeki, 2005, 1953): Other metrics include the following vertical canopy cover (VCC) metrics: canopy-tototal-return ratio (VCC r ) (Morsdorf et al., 2006), canopy-to-total-first-return ratio (VCC fr ) (Morsdorf et al., 2006), intensity-return ratio (VCC ir ) (Hopkinson and Chasmer, 2009), Beer's Law-modified-intensity-return ratio (VCC bl ) (Hopkinson and Chasmer, 2009) Zhao and Popescu, 2007), crown area (G) using detected tree heights with an empirical height-to-crown-radius function, distances and directions to nearest crown(C dist , C dir ) and canopy pixels (Cr dist , Cr dir ) from the plot center (Moeser et al., 2015), effective leaf area index (L e ) based on the Beer-Lambert Law (Korhonen and Morsdorf, 2014;Monsi and Saeki, 1953), L e based on the ground-to-total-return ratio (Richardson et al., 2009), and L e based on contact frequency (Morsdorf et al., 2006), apparent clumping index (Ω app ) (Ryu et al., 2010), and Beer-Lambert Law canopy gap fraction (P bl ) Saeki, 2005, 1953;Ryu et al., 2010). While the ITC results may not be physically meaningful in this case, as they were not locally validated, we analyze these values for correlation with T in a classical feature engineering approach.
Correlations with convex spherical densiometer measurements were calculated before testing univariate and multivariate linear models with stepwise-AIC and -BIC model selection.

Tree and crown metrics
In order to perform individual tree crown (ITC) detection and crown area estimation, empirical data from recent research in the study area (Cortini et al., 2011) was applied to model the height-to-crown-area relationship for deciduous and conifer species, as well as all species as one group. The ground data consist of aggregated minima, means, and maxima for major regional tree species height-to-crown-area, with standard deviations provided. Models for heightto-crown-area were developed for aggregated native species in the study area from these statistical moments. 23

Correcting for temporal mismatch
The effect of filtering sites likely disturbed between spherical densiometer and ALS sampling campaigns was tested, in order to correct for a half-decade mismatch in data collection.
This filtering process was also used to correct for discontinuity between ground and remote sensing observations due to seasonal changes in leaf area index, as ground observations were generally collected during summer leaf-on conditions while ALS sorties were conducted in fall leaf-off conditions. The error contribution of leaf state is likely minimal, as evergreen forest is dominant in the study area (Nielsen, 2005). Observations with ground-based angular canopy closure (ACC) values below 0.10 were filtered or removed, where disturbances or leaf condition discontinuities were apparent in ground-to-ALS ACC plots. Observations were filtered if the ground ACC value, collected at a later date (i.e., potentially subject to disturbance), was less than 0.1 and showed a reduction of 0.1 or more.

Results
Estimation of ACC and P o as a proxy for T using ALS data showed good performance.
Regression models using multiple metrics substantially outperformed any single ALS metric, yet individual metrics have utility for their simplicity and physical basis, facilitating interpretation.
Of the individual metrics, VCC fci , showed the best performance.

ALS Estimates of ACC and P o
To test for correlations, given the perfectly inverse relationship between gap fraction ( The mean R 2 improvement attributable to filtering out disturbances was ΔR 2 = +0.05. The largest gains were shown by VCC cv (ΔR 2 = +0.20), VCC ir (ΔR 2 = +0.18), VCC fr (ΔR 2 = +0.16), VCC p (ΔR 2 = +0.15), and VCC r (ΔR 2 = +0.15), while the largest loss was shown by the stereographic and equidistant fisheye lens model P hv metrics at a minimum canopy height of five meters (ΔR 2 = -0.01). Overall, VCC metrics, ITC metrics, and the equisolid angle P hv metrics showed the greatest model improvement, indicating sensitivity to disturbance-or leaf arearelated noise. Figure 3 shows the full P hv calculation process conducted for each site tested.  Applying the VCC fci calculation to the full dataset of 950 ALS and ground plots, model fit improvement is again exhibited by filtering out disturbances. Both second-order polynomial (R 2 = 0.39) and exponential (R 2 = 0.35) models show reasonable model fit before filtering disturbed sites, followed by a simple linear model (R 2 = 0.32). After filtering out disturbed sites, model fit improved for the second-order polynomial model (R 2 = 0.43), exponential model (R 2 = 0.42), and linear model (R 2 = 0.40). Thus, linear and exponential models showed the greatest improvement in model fit, which is logical given their relatively inflexible behavior compared to polynomials.
Meanwhile, P pdn showed strong linearity with ACC and thus P o ( Figure A2.5). Errors were higher at lower values of ACC, with the presence of a few strong outliers. The application of exponential and polynomial linear models were tested in terms of their impact on model performance (Table 3).

Point-density normalized canopy gap fraction
The P pdn algorithm produced reasonable results, showing agreement with other P o estimates and measurements. A visualization of point-density-normalized gap fraction (P pdn ), Beer-Lambert Law gap fraction (P bl ), and Beer-Lambert Law effective leaf area index (Le bl ), and apparent clumping index (Ω app ) are provided for an example ALS field plot ( Figure 5). Variants of the ITC detection algorithms implemented here underwent validation in a number of previous studies (Kaartinen et al., 2012;Popescu et al., 2002). The algorithms were applied to generate predictor variables to test for variable importance in estimating canopy gap fraction (P o ), and its inverse, angular canopy closure (ACC). Herein, ITC results are treated as features for estimating T, rather than tree crown counts, as the purpose was to extract additional information from ALS data. Hence, the accuracy of their results is not a consideration in this work. From a visual analysis of ITC estimates, reasonable algorithm performance is assumed. The ITC algorithms implemented include standard and hierarchical watershed segmentation, as well as standard and hierarchical variable-size moving window methods.
Standard and hierarchical variable-size moving window ITC detection counts of tree crowns performed the best in predicting ACC of the ITC methods, each with an R 2 above 0.4, despite not undergoing calibration. While ITC methods were not inferred to be able to predict ACC on their own, as ITC counts and ACC are considered dependent variables (Falkowski et al., 2008;Kaartinen et al., 2012;Wang et al., 2016), they are complimentary to other metrics as an additional feature of forest geometry, as is the apparent clumping index (Ω app ).

Discussion
While solar position, topography, and atmospheric conditions are known to effect the quantity and quality of understory light (Dengel et al., 2015), in this paper, we focus on canopy light transmission (T) indices best captured by LiDAR. This follows longstanding hemispherical photography research on canopy light transmission indices, including the gap light index or GLI/C (Canham, 1995(Canham, , 1988 and the related Gap Light Analyzer or GLA (Frazer et al., 1999), as well as recent LiDAR methods aimed at characterizing broad areas at reduced time and cost (Korhonen and Morsdorf, 2014). Our proposed LiDAR canopy light transmission indices are intended for later application with statistical (e.g., machine learning) models to capture nonlinear effects between canopy geometry, solar position, topography and atmospheric conditions on understory solar irradiation levels in large-area mapping efforts. This obviates the need for computationally expensive physical simulations at every grid cell. Although previous studies show strong agreement with ground measurements for a number of ALS metrics of forest structure (Korhonen and Morsdorf, 2014), notable challenges remain. Models of canopy light transmission often utilize physically-based ray-tracing (Disney et al., 2000), which can be thought of as a synthetic LiDAR system, or are derived from simple canopy metrics such as Lorey's canopy height or leaf-area index (Niinemets and Anten, 2009).
While the latter method lacks physical-geometric realism readily visible in existing point cloud datasets, the former also has its challenges. While radiative transfer models using ray-tracing may improve landscape-scale understory light estimates (Gastellu-Etchegorry et al., 2015;Moeser et al., 2014;Reich et al., 2012), ray-tracing requires high-point-density data (> 10 returns/m 2 ) from ALS or terrestrial laser scanning (TLS) LiDAR systems along with ancillary information beyond standard (x, y, z, intensity) information. Ray-tracing methods are also computationally demanding, making them slower and more expensive to apply. While deep reinforcement learning methods designed to accelerate ray-tracing algorithms through improved importance sampling may partially alleviate these challenges (Dahm and Keller, 2017), as demonstrated by Nvidia's latest RTX GPUs, ray-tracing remains computationally expensive.
In contrast, simple return-ratio approaches of quantifying canopy radiation attenuation may offer improved functioning with low-point-density data, simple, accelerated wall-to-wall mapping, and improved compatibility with historical ground-based methods needed to validate models with existing datasets or to analyze historical changes in forest strcuture. Furthermore, canopy attenuation-based ALS metrics may be comparable to methods used in the synthetic aperture RADAR community to estimate aboveground volume, such as the semi-empirical Water Cloud Model (Attema and Ulaby, 1978;Graham and Harris, 2003 showed that ITC detection methods provided one of the best proxies for ACC, which was unexpected and thus noteworthy. Our new P hv metric showed a low ACC R 2 saturation near 0.2 for all lens geometries even after filtering disturbed sites. Maximum R 2 values for the P hv index were consistently shown for a canopy height threshold of 5 m. Of the lens geometries tested, the equisolid angle (equiangle) projection was shown to be the most sensitive to disturbances present in the observational record.
Meanwhile, after filtering disturbed sites, differences in accuracy were more attributable to canopy height threshold than to lens model, with each lens model showing a similar R 2 pattern across tested threshold values. Meanwhile, the P pdn metric may be considered a step toward the harmonization of ground-based and airborne estimates of P o , which remains an outstanding challenge due to the different nature of ground and LiDAR measurement techniques. Finally, the excellent result for P pdn and poor results for P hv begs the question: why do simple ratio-based models continue to outperform detailed geometric models? We believe this is due to the sensitivity of highly detailed models to discrepancies in the validation data, which brings us to our study limitations.

Limitations
A fundamental limitation of this work was the half-decade difference in time between ground and ALS data acquisition, yielding strong disagreement between ground and ALS metrics of ACC for some sites. From ALS and field data scatterplots, it was apparent that disagreement arose either from disturbance or regrowth on previously disturbed sites. This temporal mismatch diminished the utility of ground ACC data for use in model validation, as shown by model performance after filtering disturbed sites. This gave rise to a second question present throughout the duration of this study: why do we still use spherical densiometers for remote sensing model validation in 2019? Although the ALS data had a low mean point density of 1.64 points/m 2 , these active data are of greater geolocation accuracy, precision, and sampling density than passive coarse spherical densiometer measurements. Presently, it would not be unreasonable to treat LiDAR itself as ground-truth data, given its superior characteristics by most metrics.
Thus, we question the use of coarse ground measurements of ACC (e.g., spherical densiometers), instead arguing for modern LiDAR systems, structure-from-motion (SfM), realtime simultaneous localization and mapping (SLAM), 360-degree spherical imagers (e.g., FLIR Ladybug), or digital hemispherical imagers. Today, the average smartphone imager provides greater information about canopy geometry than spherical densiometers, including an ability to produce 3-D SfM or SLAM point clouds and to display the produced 3-D models using built-in augmented reality (AR) interfaces running on onboard graphics accelerators (e.g., ARM Mali, Apple A12 Bionic, Qualcomm Adreno 640). The use of full-waveform data may further add state-of-the-art vertical canopy sampling and canopy penetration essential for modeling canopy light transmission. Yet, historical spherical densiometer data was essential for the completion of this study and methods will continue to be in demand that are able to cope with densiometer data for global change studies. For such applications, we provide the new P pdn metric and for detailed geometric datasets, we provide the new P hv metric.
As a result of the aforementioned data limitations, none of the P hv methods tested show strong performance, requiring further validation against hemispherical photography measurements closer to the time of ALS acquisition. This is indicated by the strong agreement between multiple LiDAR-derived predictors of ACC showing only moderate agreement with convex spherical densiometer measurements. While step-wise AIC and BIC linear regression models included high numbers of coefficients without substantial performance gains, univariate linear models showed equivalent performance. We infer that this temporal mismatch poses a fundamental limitation on algorithm performance in this study, as top-performing metrics saturate near the same accuracy level.

Conclusion
This work demonstrated two important new algorithms for modeling of forest structure applicable to multiple types of point cloud data (e.g., ALS, TLS, SfM), as well as a method for filtering disturbed sites. While our study was limited by the quality and acquisition timing of field data, we found that the P pdn metric in particular showed strong performance. In addition, we showed that filtering sites and canopy threshold height have a greater effect on P hv performance than the synthetic lens model. Meanwhile, traditional VCC metrics still showed the best overall corrrespondances to ACC measurements, despite being fundamentally different in principle.
From these results, we concluded that the new ALS-based models of T are promising, yet require further development with higher point densities closer to the time of ground data acquisition. Those with high point density LiDAR datasets may nonetheless benefit from the methods presented above, necessary for pursuing similar studies in regions where there is limited ground sampling coverage, as is often the case in boreal forests. These new metrics in turn are 41 likely to be overcome by unsupervised feature learning (i.e., deep learning) applied to highpoint-density datasets.
As point densities increase with technological advances, and spectral data are embedded to points (e.g., SfM or multi-spectral ALS systems), traditional ground measurement techniques may be less relevant for model validation. We argue that point cloud models are sufficient in their own right for the estimation of canopy geometric properties, such as coverage or closure.
Future studies should move beyond historical ground measurement techniques of canopy light transmission to explore the use of synthetic data under idealized conditions. By generating idealized point clouds of forests (evenly spaced 1 point/mm 2 ) using a latest generation 3-D simulation framework, and iterating over random samples from these, robust physical features may be engineered that function across a variety of forest conditions. Such physically-based rendering tools are also ideal for the generation of large labeled datasets needed to train state-ofthe-art supervised learning models, overcoming the central factor limiting the application of deep learning in LiDAR remote sensing of forests.