Assessment of Light Environment Variability in Broadleaved Forest Canopies Using Terrestrial Laser Scanning

Light availability inside a forest canopy is of key importance to many ecosystem processes, such as photosynthesis and transpiration. Assessment of light availability and within-canopy light variability enables a more detailed understanding of these biophysical processes. The changing light-vegetation interaction in a homogeneous oak (Quercus robur L.) stand was studied at different moments during the growth season using terrestrial laser scanning datasets and ray tracing technology. Three field campaigns were organized at regular time intervals (24 April 2008; 07 May 2008; 23 May 2008) to monitor the increase of foliage material. The laser scanning data was used to generate 3D representations of the forest stands, enabling structure feature extraction and light interception modeling, using the Voxel-Based Light Interception Model (VLIM). The VLIM is capable of estimating the relative light intensity or Percentage of Above Canopy Light (PACL) at any arbitrary point in the modeled crown space. This resulted in a detailed description of the dynamic light environments inside the canopy. Mean vertical light extinction profiles were calculated for the three time frames, showing significant differences in light attenuation by the canopy between April 24 on the one hand, and May 7 and May 23 on the other hand. The proposed methodology created the opportunity to link these within-canopy light distributions to the increasing amount of photosynthetically active leaf material and its distribution in the considered 3D space.


Introduction
The forest canopy is a central part of the forest ecosystem and acts as an interface, fulfilling the key role of cycling material and energy between vegetation and atmosphere.Light interception is a driving variable for many ecosystem processes, such as photosynthesis and transpiration, and was proven to be the factor on which canopy structure has the most pronounced impact [1].The interaction between stand structure and light is fundamental to understanding biophysical and ecological processes, such as forest growth, forest regeneration, and understory biodiversity.In addition, a detailed analysis of the within-season dynamics of these processes requires the measurement and modeling of a time series of light-vegetation interactions.Light interception models for plant canopies differ widely, depending on their detail in canopy structure representation and intended use.Model complexity ranges from simple 1D turbid medium models to complex geometrically explicit 3D ray tracing models that rely on a detailed structural description of the forest canopy [2].Light models with an explicit crown architectural representation use ray tracing technology to track simulated solar rays.These rays originate at a given point in the upper hemisphere and then travel through a geometric representation of the forest canopy [3].When interacting with a canopy element, a fraction of the radiation energy is absorbed, transmitted, or reflected.The path of the reflected radiation can be tracked through the canopy, resulting in second-order and higher-order estimations [1,2].Ray tracing technology is limited by the technical difficulties of acquiring realistic 3D descriptions of the object of interest.A 3D structural description of vegetation is not only tedious and time consuming, but also difficult to obtain in a non-destructive manner for large trees or entire forest stands.This makes time series of detailed structural assessments of full grown forest stands during a growth season very difficult.
Several techniques have been documented describing the acquisition of detailed 3D structural information of small to medium sized canopies (i.e., <1.5 m).Stereo-photography uses a stereo imaging system to create partial 3D models by digitizing the upper layer soybean canopy leaves.Two photographs from the canopy are simultaneously taken from known positions, and corresponding points on each image are matched, allowing the determination the points' 3D coordinates [4,5].This allows the acquisition of structural information on the level of single leaves (e.g., 3D leaf position, leaf inclination angles).Frasson et al. [6] presented a similar methodology using stereoscopy to develop highly detailed 3D digital models of maize canopies.Still, applying the method to full grown forest stands with dense canopies is difficult because of occlusion effects and the difficulty of matching corresponding points.Alternatively, Sinoquet et al. [7] described the use of an electromagnetic digitizer to map the position and shape of leaves and other vegetative elements.The digitizer measures the spatial co-ordinates of any point within a given volume by creating three perpendicular magnetic fields, in which the location of a pointer can be determined.While this 3D structure information was used to characterize the light environment in small canopies [8], the technique is tedious, time consuming, and often not precise enough to accurately capture detailed leaf geometry.For example, it may not capture measure the leaf edges [9].Also, the technique is limited to small canopies because of its maximal range (<4 m).
In the case of full grown forest stands, terrestrial laser scanning systems are capable of producing the necessary structural information for geometrically explicit models with a level of detail that depends solely on the resolution of the laser scanning dataset.Omasa et al. [10] performed pioneering work in the fusion between airborne and terrestrial laser scanning datasets to obtain the highest level of detail possible with current technologies.The requirements of geometrically explicit light models are best met when the actual 3D distribution of the above-ground density of phyto-elements (leaf, stem, twig, etc.) is available.This 3D structure is often represented by the leaf area density (LAD) per volume unit, where LAD is defined as the total one-sided leaf area per unit volume [11].LAD estimates can be computed from multi-angular laser datasets using the inclined point quadrate method by [12,13] or by gap fraction inversion, which is similar to the LAI determination process (for forest stands) using hemispherical photography.Voxel-based methods for LAD estimation have been presented by [14].With a sufficient level of detail, the LiDAR voxel size will approximate the leaf size, resulting in an accurate description of the actual leaf distribution.The detailed structural description can then be used as the input for a ray tracing program to determine the 3D light environment.
The main objective of this paper focuses on the description of the changing light-vegetation interaction in a homogeneous oak (Quercus robur L.) stand at different instants in the growth season using LiDAR datasets and ray tracing technology.With this approach, the potential of the innovative 'Voxel-based Light Interception Model' (VLIM), a light interception model designed to estimate the Percentage of Above Canopy Light (PACL) at any given point in the forest scene, is tested for multi-temporal purposes.

Study Site
An oak (Quercus robur L.) stand near Antwerp (Belgium) was selected as a test site.The stand was characterized by a uniform upper-story vegetation of evenly aged trees with a planting distance of 9 m by 9 m.The average tree height was approximately 17 m.No understory vegetation was present.

Structure Data Acquisition
The terrestrial laser scanning device used in this case study consisted of a commercially available Laser Measurement System (LMS200, Sick AG, Germany).The LMS200 device has a non-contact optical active sensor, which was mounted on a dynamic measurement platform, enabling measurements in a 3D hemispherical pattern.A full description of the characteristics of the laser device and the hemispherical platform is given in [15].In case of a first return laser device like the LMS200, the spatial information related to the position of the vegetation elements located behind the target object is not available due to the shadow effect.To circumvent this problem, these background elements have to be measured from different angles in order to obtain comprehensive laser coverage [10,16,17].The scanner was positioned in a diamond shaped setup.This consisted of a single hemispherical laser measurement in the center of each plot and one in each of the four cardinal points (N, E, S, W), with a displacement of 0.5 m from the center.Four measurement plots were systematically distributed throughout the forest stands at 20 m intervals.Additionally, digital hemispherical photos, based on the approach of [11], were taken in each measurement plot to estimate the Leaf Area Index (LAI).A total of four photos were collected in a regular pattern at 10 m intervals.Three field campaigns were organized in the growth season at regular intervals (24 April 2008;07 May 2008;23 May 2008) to capture the increase of foliage material.

Data Pre-Processing
The laser scanning datasets were used to create a detailed structural representation of the measured forest canopies.The scale at which this 3D representation can be determined strongly depends on the laser beam distribution inside the forest canopy.The laser scanning datasets, which were gathered according to the predefined sampling design, required preprocessing before the 3D structure of the measured forest plot could be extracted.The primary processing chain consisted of four different steps: (i) Transformation of the LiDAR data from a polar to a Cartesian coordinate system.(ii) Registration of the separate datasets per measurement plot to one central coordinate system, enabling the merging of the datasets into a comprehensive scan of the 3D scene with minimal shadowing effect [16].(iii) Transformation of vector data to raster data with cubic voxels.The 3D space, considered with a ground surface of 30 m × 30 m and a height of 25 m, was subdivided into cubic voxels of 0.1 m × 0.1 m × 0.1 m (i.e., 'small' voxels).This resulted in arrays of 300 × 300 × 250 voxels.The voxel size was set at 0.1 m to obtain sufficient detail and to avoid the use of a statistical description of the leaf distribution, as the voxel size resembles the actual leaf size.All laser beams emitted from the different measurement positions were traced through the voxel array representing the forest stand.Following the methodology of [14], the voxels were characterized depending on beam/voxel interaction.Voxels with at least one intercepted laser beam were given the attribute F (filled).Attribute E (empty) was assigned to voxels that were intersected by laser beams without interception.Attribute X (no data) was granted to voxels that were not touched by any laser beam.While the diamond setup was used to minimize the amount of voxels with attribute X, shadowing is always present when dealing with a first return laser scanner.To obtain an accurate and complete structure description of the canopy from LiDAR datasets, these 'X' voxels needed to be given an attribute F or E. 1,000 small voxels were grouped into large voxels of 1 m × 1 m × 1 m.For each of these large voxels, the contact frequency (CF) was calculated.
The CF can be defined as the proportion of filled voxels (F) to the measured voxels (F + E).For each large voxel of 1,000 small voxels, a proportion of 'X' voxels equal to the calculated CF was given attribute F. These 'X' voxels were chosen randomly in the large voxel.(iv) Differentiating voxels representing stem material from voxels filled with leaves by using a specially developed 'stem locator' algorithm.The algorithm is based on the density distributions of vertically projected leaf/beam intersections.As stems are characterized by a strong vertical structure, one can assume that the vertical organization of the laser hits reflected by stems will form dense clusters in the horizontal projection, thus revealing the stem positions.A cone with its base, height, and taper equal to the mean characteristics of the stems of the considered forest plots is positioned at each of the determined stem locations.Voxels contained in this cone are identified as 'stem' voxels and the remaining voxels as 'leaf' voxels.
By preprocessing the laser scans, a comprehensive 3D raster image of each measurement plot was obtained, where each single voxel (side of 10 cm) was characterized by an ID dependent on the contents of that particular voxel (i.e., leaf material, stem material).
The hemispherical photos were analyzed using an IDL image-processing software package for Windows (RSI, Version 4.0, CO, USA).As absorption of leafy material is maximal and light scattering tends to be highest in the blue channel, the blue channel was chosen for subsequent analysis.Image binarization was performed using manual thresholding.An experienced operator selected for each image, individually, the best visual threshold, in order to discriminate foliage from the sky background.This image was then used to calculate the LAI using the gap fraction inversion methodology.The assumption of randomly distributed foliage allowed the use of the Poisson model.The LAI values for each plot were averaged to enable comparison with the LAI values obtained with the terrestrial laser scanning system.

Leaf Area Distribution (LAD)
To visualize the increase in leaf biomass throughout the growing season, a detailed LAD distribution of each of the measured plots was determined for the different time frames.This allowed a description of the increase in photosynthetically active material during the growth season.To obtain accurate LAD estimates from the voxel-based LiDAR image, voxels containing non-photosynthetic material had to be separated from voxels containing leaves.LAD values for volumes of 1 m 3 (big voxels = 1,000 small voxels) were reproduced using the calculated contact frequency and the following equation [14]:

LAD (voxel) = cos (θ beam ) / G (θ beam ) CF (voxel)
where θ beam represents the mean zenith angle for all laser beams penetrating the voxel.G(θ beam ) is the mean projection of a unit leaf area on a plane perpendicular to the direction of the laser beam at θ beam [11,[18][19][20].The term cos(θ beam )[G(θ beam )] −1 describes the probability of a laser beam that enters a filled voxel of actually hitting the leaf or leaves inside that voxel as a function of laser beam and leaf inclination angles.This term corrects for the situation where voxels were wrongfully given the attribute E while they were actually filled.

Voxel-based Light Interception Model (VLIM)
The VLIM is a light interception model designed to estimate the Percentage of Above Canopy Light (PACL), or the Relative Light Intensity (RLI), at any given point in a forest scene, using a voxel-based representation of the trees derived from ground-based LiDAR data.The term, light, will be used for photosynthetically active radiation (PAR).The model is created by simulating virtual rays of light in the 3D scene and registering the interaction of those beams with the vegetative elements represented as voxels (Figure 1).Voxels representing leaf material were abstracted by disks with a fixed area of 0.01 m 2 , a random azimuth angle, and a fixed zenith angle corresponding to the average zenith angle of the trees (Figure 1c).All leaves were assigned a constant reflectance of 6.5% and a transmittance of 6.5%.These data were derived by averaging measured reflectance spectra of healthy leaves and trunks over the 400-700 nm range (own measurements, 2008, ASD spectroradiometer with leaf probe, 10 leaves/set, not published) and are supported by the findings of [21].All objects were assumed to be Lambertian scatterers.Voxels representing stems were abstracted by a participating media (volume grid) and placed in a 3D grid, of which each stem voxel was assigned an infinite optical thickness, while other empty voxels were assigned zero optical thickness.The forest scene illumination was composed of a combination of direct and diffuse light that was modeled to closely agree with the average circadian illumination for a time period of 14 days per measurement date (Figure 1d).For each time period, the average hourly direct and diffuse irradiance and solar positions for Brussels, Belgium (lat 50.85, lon 4.37) were obtained from Meteosat data recorded from 1996 until 2000 (www.satellight.com).
Per light simulation, 40,000 voxels (i.e., leaves) in the measurement plots were randomly sampled for light interception.This sampling procedure allowed the light conditions in the 3D space to be defined in great detail.This allowed these light conditions to be up-scaled to any larger scale.The resulting light distribution can then be directly used in dynamic models describing biophysical processes, like photosynthesis and transpiration.

3D Canopy Representation
An accurate geometrical description of the measured forest plots was imperative in obtaining high quality light regimes using VLIM [16].Figure 2a visualizes the voxel-based representation of a measured forest plot.This figure illustrates the capacity of the laser scanner to capture the actual structure built-up of the forest plots.The comprehensive structural dataset was used to describe the 3D Leaf Area Density (LAD) distribution inside the forest canopy.The LAD was calculated for volumes of 1 m 3 and the vertical profile of the LAD distribution is presented in Figure 2b.The extracted LAD distributions show a significant increase in the amount of leaf material between the first two time frames.From dates two to three, only a small increase between a height of 9 m and 13 m, is observable, which indicates that the forest plots were reaching their maximum foliage capacity by May 23, 2008.Integrating these LAD values (Table 1) for the complete plots resulted in mean LAI values of 3.13, 4.94, and 5.38 for dates 1, 2, and 3 respectively.Due to the lack of detailed reference data, an overall validation of the structural description was obtained by LAI comparison.The LAI values yielded by hemispherical photography are presented in Table 1.In average, LAI values of 1.34, 3.42, and 3.45 were obtained for dates 1, 2, and 3 respectively.It has to be stated that the use of hemispherical photographs is not ideal for this type of validation, as the laser scanner provides information of a volume of 30m × 30m × 25m while the hemispherical photographs consider a much bigger volume LAI value calculation.This makes a direct comparison of the values difficult.Still, the relative LAI increase between dates 1 and 2 and stagnation between dates 2 and 3 can be compared with the results from the laser scanner data, and show similar behavior.Table 1.LAI values calculated for the four different plots using terrestrial laser scanning (TLS) and hemispherical photography (HEM).

VLIM
Using the extracted structure data, in combination with the VLIM, light interception datasets of the three time frames were calculated for each measurement plot.The sampled leaves for each light interception dataset were classified into five classes, according to their average incoming global radiation, compared to the initial incident radiation at the top of the canopy: 0-20%, 20-40%, 40-60%, 60-80%, 80-100%.This enabled a direct insight into the light use efficiency of the upper-story oak canopies during the considered time frames.The level of detail of the model light environment inside the canopy is illustrated in Figure 3a, which shows a measured plot of the first date.No light extinction models, like Beer-Lambert's law [22] were necessary for determining the PACL at the sampled location.Abstract leaves were used instead; this was possible because actual leaves coincided with the voxel size.The proposed methodology created the potential to study the impact of actual crown geometry, leaf distributions, leaf angle distributions, etc. on the light use efficiency of the canopy.When looking at the position of the light interception classes inside the canopy, it is clear that light penetrates differently from what could be expected from layer-based models, where the distribution of the vegetative elements is assumed to be homogeneous.The distribution of the light classes was not only dependent on the solar zenith and azimuth angle, but more importantly, on the structural buildup of the canopy.
Shading by neighboring crowns is a major constraint on the available light of trees in the forest stand.The shading was also simulated as realistically as possible.The considered voxel volume (30 m × 30 m × 25m) was divided into two sub volumes: (i) the central and (ii) the buffer zone surrounding the central zone.The volume of the central zone measured 15 m × 15 m × 25 m.The width of the buffer around this central area was 7.5 m in all directions.Per light interception calculation, a total of 40,000 randomly chosen leaves in the central zone of the forest stand were sampled to obtain a detailed overview.
Mean vertical light extinction profiles were calculated for the three time frames (Figure 3b); the profiles show significant differences in light attenuation within the canopy between April 24 on the one hand, and May 7 and May 23 on the other hand.These results are to be expected.The vertical LAD profiles (Figure 2b) show the presence of a larger amount of photosynthetically active leaf material for these last two time frames in comparison to the first one.When evaluating time frames 1 and 2, it can be observed that the difference in light attenuation is largest between 9 m and 19 m and decreases below 9 m.This is caused by the fact that the increase in leaf surface in the upper-story oak canopy between April 24 and May 7 is larger than the biomass increase in the understory.This higher light attenuation in the upper-story of the canopy causes the vertical light extinction profile to initially shift to the left, resulting in a lower amount of light reaching the understory.The direct link between the presence of leaf material and light attenuation is also observable between the mean light extinction profiles for May 7 and May 23.While the differences in LAD and PACL distributions are minimal, the small increases in LAD between 9 m and 13 m on the one hand, and 2 m and 6 m on the other hand, cause a proportional decrease in light penetration for the last date.
As expected from the LAI values for each stand, the PACL at the bottom of the stand is directly linked to the amount of leaves in the canopy.The increase in PACL at lower heights is caused, in combination with a low local leaf density, by an increased light penetration at low solar positions.

Conclusions
The proposed methodology creates the possibility of studying photon-vegetation interactions at a high level of detail without expensive equipment or labor-intensive direct measurements.The speed and ease of structural data collection with the LiDAR device enables the acquisition of detailed structural descriptions of forest stands, in a non-destructive manner.This allows one to map the seasonal structural changes.Low resolution ground-based LiDAR technology was used to obtain a detailed structural representation of the forest scenes during three different times in the growing season for light interception modeling purposes.A processed LiDAR image consisted of an array of voxels with following attributes: (i) filled (leaf or stem) and (ii) empty.This LiDAR image was transformed into an abstract forest scene of discs (i.e., leaves) and filled voxels (i.e., stems and branches), which were imported in a ray tracing algorithm to finally render the light interception through the canopy.Light interception of the forest canopy at near leaf level scale can be directly implemented in the canopy research of biophysical processes, like photosynthesis and transpiration, by linking the light availability to leaf level models [23].As general illumination depends on geographical location and climate conditions, the VLIM was designed so that different initial light conditions could be used.This flexibility enabled instantaneous light interception estimates, as well as average light distributions for entire seasons.Future work will focus on method optimizations, especially in increasing the accuracy of the voxel-based representation of the forest canopies derived from LiDAR measurements.Work performed by [24] can provide additional options to optimize the VLIM.The described technique is able to accurately measure the true shape of the leaves, and how the leaves change depending on incoming light quantities.This is done by using a double digitizing method, combining a hand-held electromagnetic digitizer and a non-contact 3D laser scanner.This approach enables the construction of 3D virtual trees with non-planar leaves.Such information could be incorporated into, and improve, the abstract representation of the voxel space (as it is used in the VLIM).Also, as shadowing will be an ever present factor, methods to correct this factor will be further investigated.

Figure 1 .
Figure 1.Schematic overview of the Voxel-based Light Interception Model (VLIM): The forest stand consisted of several trees (a), the leaf distribution of the canopy (b) was represented by a collection of voxels based on the acquired LiDAR data (c), and each of the filled leaf voxels was abstracted by a disc for ray tracing purposes (d).

Figure 2 .
Figure 2. (a) Voxel-based representation of a studied forest plot from 24 April 2008 after preprocessing: 'leaf' and 'stem' voxels are presented in grey and black, respectively.(b) Vertical Leaf Area Density (LAD) profiles for three different time frames.

Figure 3 .
Figure 3. (a) 3D visualization of light interception calculated with the VLIM for one measurement plot on April 24, 2008.(b) Mean vertical light extinction in terms of Percentage Above Canopy Light (%) for the three subsequent time frames.