Measurement of fine-spatial-resolution 3D vegetation structure with airborne waveform lidar: Calibration and validation with voxelised terrestrial lidar

Vegetation structure controls habitat availability, ecosystem services, weather, climate and microclimate, but current landscape scale vegetation maps have lacked details of understorey vegetation and within-canopy structure at resolutions finer than a few tens of metres. In this paper, a novel signal processing method is used to correctly measure 3D voxelised vegetation cover from full-waveform ALS data at 1.5 m horizontal and 50 cm vertical resolution, including understorey vegetation and within-canopy structure. A new method for calibrating and validating the instrument specific ALS processing using high resolution TLS data is also presented and used to calibrate and validate the ALS derived data products over a wide range of land cover types within a heterogeneous urban area, including woodland, gardens and streets. This showed the method to accurately retrieve voxelised canopy cover maps with less than 0.4% of voxels containing false negatives, 10% of voxels containing false positives and a canopy cover accuracy within voxels of 24%. The method was applied across 100 km2 and the resulting structure maps were compared to the more widely used discrete return ALS and Gaussian decomposed waveform ALS data products. These products were found to give little information on the within-canopy structure and so are only capable of deriving coarse resolution, plot-scale structure metrics. The detailed 3D canopy maps derived from the new method allow landscape scale ecosystem processes to be examined in more detail than has previously been possible, and the new method reveals details about the canopy understorey, creating opportunities for ecological investigations. The calibration method can be applied to any waveform ALS instrument and processing method. All code used in this paper is freely available online through bitbucket (https://bitbucket.org/StevenHancock/voxel_lidar).


Introduction
The 3D structure of vegetation canopies is a key determinant of ecological function and processes, providing an indicator of habitat (Ashcroft et al., 2014), biomass (Calders et al., 2015), impacting on weather and climate (Ni-Meister and Gao, 2011) and modulating microclimate (Clinton, 2003). For example, in urban systems the pattern and distribution of greenspace mitigates the "heat island" effect (Myint et al., 2015), with implications for human health.
The distribution and quality of greenspace affect mental wellbeing, directly and by providing corridors for wildlife (Vaz et al., 2015;Shanahan et al., 2017). Understanding and quantifying how characterising 3D vegetation structure over large areas but can be used to calibrate and validate larger scale measurements (Hopkinson et al., 2013).
Airborne laser scanning (ALS) measures the location and radiometric properties of reflected laser light over landscape scales, allowing the characterisation of 3D structure. They operate in two different modes, "discrete return" and "waveform". Discrete return uses proprietary algorithms to produce a point cloud (Disney et al., 2010). This allows measurement of canopy height  and has been used to estimate canopy density from the ratio of points returned from the canopy and ground (Stark et al., 2012). However, these algorithms have been developed for measuring hard targets and can be biased over vegetation (Disney et al., 2010), requiring ground based calibration to correct . In addition the return strength may not be related to target reflectance , complicating its use in canopy characterisation. These discrete return instruments methods only return a few (around 4) points per laser shot with no way of knowing what is not being measured (Gaveau and Hill, 2003;Disney et al., 2010), potentially preventing the measurement of within-canopy and understorey structure.
Full-waveform lidar measures the reflected laser intensity as a function of range (Baltsavias, 1999). This gives information on all objects visible to the ALS but requires processing to extract target properties from the signal . Fig. 1 illustrates how an ALS waveform is made up of the vertical distribution of objects that are to be measured, referred to as the "target profile", ( Fig. 1 (a)), attenuation as laser light is blocked by targets (black line in Fig. 1 (b)), blurring by the lidar system pulse (black line in Fig. 1 (c)) and noise to give the measured signal (red line in Fig. 1 (d)). The effects of noise, system pulse and attenuation must be removed in order to measure high resolution (<2 m) vegetation structure. The extra information available to waveform lidar has been used to measure leaf area index (Hopkinson et al., 2013), gap fraction (a) Lidar beam and true target profile  (Musselman et al., 2013), biomass (Drake et al., 2002), land cover (Mallet et al., 2011) and plot scale vertical foliage profiles (Harding et al., 2001). Past waveform studies have either not corrected for the system pulse, limiting vertical resolution to around 2 m (Harding et al., 2001;Fieber et al., 2015), which is too coarse a resolution to detect short understorey, or have aggregated to tens of metres horizontal resolution (Hopkinson et al., 2013), which is too coarse to detect shrubs and under-canopy paths, both of which are ecologically important. Some previous studies have used discrete return ALS to investigate within-canopy and understorey structure. They have attempted to overcome the sampling issues by either processing at coarse resolution (>20 m) (Martinuzzi et al., 2009;Hilker et al., 2010;Miura and Jones, 2010;Korpela et al., 2012), using data collected during leafoff periods (Hill and Broughton, 2009) or used statistical approaches in sparse canopies (Wing et al., 2012). Plot scale metrics do not capture the fine scale detail needed to understand the effects of structure, particularly in heterogeneous areas such as towns. Leafoff measurements are not possible in many situations (e.g. evergreen vegetation) and their use assumes that the understorey structure is constant as the overstorey changes, which is not the case for deciduous understorey. The study of Wing et al. (2012) estimated understorey density at 2 m resolution, but in a sparse canopy (from Fig. 9 in Wing et al. (2012) the maximum canopy cover was 80% with a mean around 25%), so it is not known whether their approach will work in a closed canopy. Their work also used a statistical approach, deriving metrics such as number of lidar returns from different height bands and canopy top roughness, then empirically relating those to ground measurements of understorey density. A similar statistical approach was used by Miura and Jones (2010) and Latifi et al. (2015). This assumes that understorey density and overstorey structure are correlated and requires local calibration. It is questionable if this approach can be implemented over large or heterogeneous areas, where management and disturbance regimes may vary.
Full-waveform lidar is a direct measure of light reflected from throughout the canopy and so, after correcting for instrument noise and attenuation, offers the potential to directly measure the complete vertical structure at the resolution of the ALS footprint density (typically 1 m or better). This paper aims to fully characterise vegetation structure over large, heterogeneous areas using physically based waveform ALS inversion. Detailed, plot scale TLS measurements were used to calibrate the ALS instrument specific parameters and validate the ALS processing. The sensitivity of the results to the calibrated parameters was tested. The ALS processing method was then applied to a dataset covering 100 km 2 . Table 1 Plot descriptions along with number of TLS scans.

Plot
Description Scans   1  Woodland with dense understorey  3  2  Woodland with medium understorey  3  3  Well manicured garden  2  4 Road and front garden with shrubs 2 5 Roadside with isolated trees 2 6 Grass field transitioning to dense bushes 2 7 Car park, front lawn and isolated tree 2 8 Lawn divided by hedge with large oak tree, ringed by scrubby bushes 2 2. Materials and methods

Field site and lidar data
Luton in the UK, a predominantly Victorian terraced town with areas of woodland, scrubland and parkland, was surveyed by the NERC-ARSF Dornier 228 aircraft in September 2012. This carried a Leica ALS50-II ALS with the waveform WDM65 add-on, recording data as separate discrete return and waveform streams. The discrete returns used Leica's proprietary algorithm to record up to four returns per laser shot with a reported accuracy of 20 cm (Kukko and Hyyppä, 2007). The waveform stream recorded the returned laser intensity every 1 nm (15 cm) for an extent of 38.4 m after the first return (with a buffer to capture the leading edge) per laser shot. The system pulse (the combination of laser pulse shape and detector response) was measured from returns from Luton airport runway (low reflectance) and a grass football field (high reflectance) and was found to have a width (one standard deviation, s) of 53 cm (shown in Fig. 3 Hancock et al. (2015)), causing a blurring of 2 m. The laser wavelength was 1064 nm. The aircraft flew at an altitude of 2.6 km and airspeed of 140 knots, giving a footprint of around 33 cm in diameter with a density of between 0.25 and 4 footprints per m 2 depending on ground elevation, scan angle and flightline overlap.
Aerial photos and lidar waveform shapes were examined in a 1 km 2 area around a site of intensive ecological surveys (Cox and Gaston, 2016). Eight ground plots were selected to cover the range of observed land covers across the urban extent (which included woodland, isolated trees, flower garden, scrubland, road, building, grass field and hedges) and ALS waveform shapes (in terms of maximum return extent and number of separate returns per beam). These plots are shown in Fig. 2 and described in Table 1. Ground data were collected in August 2014 with a Riegl VZ400 TLS. An angular resolution of 0.35 mrad was used with up to four discrete returns per laser shot recorded. The laser had a wavelength of 1545 nm (different to the ALS). The reported range accuracy was 5 mm. Between two and three separate hemispherical scan locations were collected per plot in order to avoid errors from attenuation, with two for sparse vegetation and three for denser vegetation. Scan locations were manually selected to give a clear view to each side of every tree crown within a plot (Hancock et al., 2014). The TLS was geolocated by manually aligning roofs of buildings to the discrete return ALS. This gave an accuracy of 10 cm vertically and 30 cm horizontally.

Lidar processing
The ALS processing is described in Section 2.3, the TLS processing in Section 2.4 and Section 2.4.2 describes how the TLS data was used to calibrate and validate the ALS processing. A flowchart of the ALS and TLS processing methods and how they were calibrated and validated is shown in Fig. 5.

ALS processing
To characterise vegetation at high spatial resolution (<2 m) from waveform ALS, the target profile (red line in Fig. 1 (d)) must be extracted from the measured lidar waveform (black line in Fig. 1 (d)). Noise must be removed (extract 1(c) from Fig. 1 (d)), the system pulse blurring removed ( Fig. 1 (c) to (b)), then attenuation must be corrected ( Fig. 1 (b) to (a)) to retrieve the true foliage profile (defined as the canopy cover within each vertical bin).
The most common waveform processing method is to denoise the signal by subtracting a noise threshold, then decompose the signal by fitting Gaussians (Hofton et al., 2000;Wagner et al., 2006), also referred to as Gaussian decomposition. This gives a discrete point cloud where each point has an associated width and amplitude, giving details on canopy structure (Mallet et al., 2011), however it assumes that the target profile can be represented as a sum of Gaussians. Previous work carried out by the authors suggests that this assumption is valid for large footprint lidar (30 m diameter), for which Gaussian decomposition was developed, but may not hold for small footprint (33 cm) ALS due to increased heterogeneity . Alternative methods have been proposed to retrieve the target profile by deconvolving waveforms. Jutzi and Stilla (2006) applied the Wiener filter, Parrish and Nowak (2009)  This paper tests the ability of Gold's method (referred to as "Gold" from this point onwards) to characterise vegetation canopies and benchmarked against the more common Gaussian decomposition method and discrete return ALS data. Richardson-Lucy deconvolution was also tested, but the results for Gold were found to be slightly more accurate (1% difference) and a factor of ten less computationally expensive than Richardson-Lucy deconvolution due to faster convergence, so only the results from Gold are given in this paper. A comparison with the other deconvolution methods was beyond the scope of this paper. Sections 2.3.1 to 2.3.4 describe the denoising, system pulse removal and attenuation correction steps needed for a direct, physically based measurement of understorey and overstorey density.

Denoising
Noise can create false returns and distort true ALS returns. Fig. 5 Hancock et al. (2015) shows that the Leica ALS50-II has stable background noise, therefore it can easily be removed by ignoring all signal beneath a threshold (Hofton et al., 2000). There are multiple combinations of denoising methods, each of which has potential advantages and disadvantages. As these have not yet been assessed in as much detail as the combination of ALS and TLS data allows, six combinations were assessed in order to achieve the highest accuracy possible. The threshold can either be fixed or set by the statistics of the signal (Hofton et al., 2000), referred to as "fixed" and "variable" thresholds. The variable noise threshold was set as the modal value of each waveform plus a multiple of the mode of the deviation about the mode, "threshScale" in Table 2. The mode was used, rather than the mean of a section of blank waveform, due to the short extent of the waveforms (38.4 m). To avoid occasional spikes in background noise, only features above a minimum width were accepted, minWidth. If this width was set to one, all signal above the noise threshold would be accepted.
Distortions in the signal could be removed by convolution with a Gaussian, either before or after the removal of background noise, referred to as "pre-smoothing" and "post-smoothing". Smoothing could be removed by setting the width, "smooWidth" in Table 2, to 0. Thus there were four denoising parameters (Table 2).
Using a hard noise threshold will truncate the leading and trailing edge of the signal, potentially preventing the accurate deconvolution of weak signals (such as ground returns beneath dense canopies). This can be avoided by tracking from the point at which the signal crosses the threshold to the mean noise level (Hancock et al., 2011). Therefore the threshold can either be "noise tracked" or "hard". When using noise tracking with fixed noise thresholds, the modal noise level for the Leica ALS50-II was found to be DN = 13, as reported in Hancock et al. (2015). These three denoising options were combined to give six ALS processing methods, listed in Table 3.

System pulse
Two methods for removing the blurring effect of the system pulse were tested, Gaussian decomposition and Gold's method. Gaussian decomposition assumes that the denoised waveform (A vis ) is made up of a sum of n Gaussians (Hofton et al., 2000), given by: where; r is range from the lidar l i the range of the centre of the ith Gaussian and s i is the width of the ith Gaussian.
The MINPACK implementation of the Levenberg-Marquardt optimisation method (Garbow et al., 1980) was used to decompose denoised waveforms into Gaussians. Inflection points of a smoothed waveform were used to decide upon the number of Gaussians and estimates of the initial parameters (Hofton et al., 2000). This does not extract the target profile from the measured waveform, instead extracting target properties as a sum of Gaussians.
Gold's method attempts to deconvolve the system pulse from the measured waveform to retrieve the target profile but deconvolution operators are notoriously sensitive to noise. To overcome this Gold (1964) proposed iteratively reblurring and deconvolving the signal. This reduces the effect of noise whilst ensuring that no deconvolved signal appears outside the bounds of the original, raw waveform. Gold's method is given by (Jansson, 1997): where; i is the raw, measured waveform s is the ALS system pulse o (k) is the kth estimate of the visible target profile o (0) = i, the initial estimate and ⊗ is the convolution operator  Deconvolution is mathematically ill-posed; there are multiple possible solutions. Wu et al. (2011) proposed iterating until the root mean square difference between subsequent estimates falls below a tolerance (deconTol in Table 2).

Hard targets
Testing revealed that deconvolution was unable to retrieve the Dirac-delta functions of returns from single, hard objects, such as bare earth. Instead, the retrieved profiles were spread over several bins, which would be interpreted as short vegetation where none exists. To remove this error, waveforms with returns from hard targets were not deconvolved and instead replaced by a single point located at the centre of gravity (Jutzi and Stilla, 2005). Waveforms were identified as coming from a hard target if there was a single feature above noise and either correlated strongly with the system pulse shape (within 4.6% RMSE) or was narrower than the system pulse width. For a 33 cm diameter footprint (s=8.25 cm), slopes will not significantly broaden the ground return. Even a 60 o slope would only increase the ground return width from the system pulse of 53 cm to 55 cm. Large footprint systems (30 m diameter) may need an alternative approach to detect hard targets in the presence of slope. Some waveforms contain diffuse targets and finish with a hard target (such as the ground) and so are not identified as hard targets by the above method. It was observed that some processed waveforms contained subterranean signal, which would have been due to a combination of multiple scattered light and electronic noise after true returns. Whilst it may be possible to filter this out using knowledge of the expected increase in the mean noise level after true hits, this was beyond the scope of this paper. Instead the false returns were removed by ignoring all returns beneath a likely ground surface. This was defined by fitting a polynomial plane (5th order polynomial in x and y fitted over 50 m by 50 m patches) through the discrete return lidar points that were identified as ground by lastools (Isenburg, 2011). This solved for the case of diffuse targets over ground (a common occurrence), but it may be the case that some hard targets with overhanging diffuse objects were not at ground level, such as roofs or cars, and in these cases there could be false positives between the hard target and ground elevations. There were no examples of this within the field sites to allow this to be investigated.

Attenuation correction
Correct denoising and deconvolution will give the target profile visible to the lidar instrument, A vis (r). The true target profile is related to this by attenuation, described by: where; A p (r) is the target area within a lidar beam projected in that direction at range r and gap(r) is the gap fraction up to range r.
This assumes that the visible target area at a given range is representative of the obscured target. The true area can be found by dividing the deconvolved signal by gap fraction; equal to one minus the cumulative visible area up to that point: This is a more direct, physically based method than the statistical approach used in previous studies (Wing et al., 2012). As it is known that the total cumulative gap fraction is equal to one when looking down from above (the ground will block all light), only the relative effective reflectance of target elements (e.g. ground and canopy) is needed. It can either be assumed that this is homogeneous throughout the canopy, in which case the A vis (r) waveform's integral can be normalised to unity and the area directly calculated, or an attempt to calculate the canopy and ground reflectance using the method of Armston et al. (2013) can be made. Whilst this method has been shown to perform well over forests, it has not been tested over the more heterogeneous urban areas. Therefore this paper assumes homogeneous effective reflectance.

TLS processing
To calculate canopy cover and generate waveforms from the TLS data, the point cloud was converted into a set of spheres with the radius representing the projected area of the target. Rays were traced through the spheres to calculate the gap fraction in a given direction, following the method presented in Hancock et al. (2014). Sphere radii were set from the lidar beam radius, accounting for partial hits and attenuation. This is slightly different to the voxel gap fraction methods used by Hosoi and Omasa (2006) and Béland et al. (2011) as the explicit geometric information is retained and will account for clumping and angular distribution effects. The setting of TLS point radii, r p , is illustrated in Fig. 3 and described by Eq. (5) (note that the equation in Hancock et al. (2014) left out the square root in error).
where; r b is the beam radius at that point q app is the target's apparent reflectance (described in Eq. (6)) I r is the TLS return intensity, corrected for range and scaled to lie between 0 and 1 and P gap is the gap fraction up to that point.
The apparent target reflectance, q app , is the product of the target's reflectance, q t , and the phase function at that view angle, X(h), as darker targets or targets at high angles of incidence will give weaker returns whilst still filling the field of view. Current TLS data is not of sufficient resolution to determine the angle of incidence throughout the canopy and so h cannot be calculated, so no attempt was made to separate the elements of q app and it was used as an effective parameter. The q app and I r terms account for partial hits, following Hancock et al. (2014), as shown in Fig. 3 (b).
To correct for attenuation of the TLS beams through the canopy, the radius was divided by the gap fraction up to the voxel containing that point, P gap . Fig. 3 (a) illustrates the attenuation correction (P gap ); for the voxel outlined in green, the lidar beams in blue reach the voxel (either pass through or hit targets within) whilst those in red would have passed through but were blocked by earlier elements. The gap fraction to that voxel, P gap , is then the ratio of the number of beams that reach the voxel (blue) to the number that could have passed through (blue plus red). In order to avoid errors by dividing by small numbers, a minimum gap fraction that could be trusted, minGap, was specified and P gap was not allowed to be smaller than minGap.
To reduce computational expense, TLS points were read into memory and mapped into voxel space. Only TLS points within voxels of interest were searched through in subsequent analysis, using an efficient voxel intersection test (Amanatides et al., 1987). TLS waveforms were generated by tracing rays from each TLS point within an ALS beam in the direction of that ALS beam (Hancock et al., 2014), rejecting rays that were obscured by nearer points. This is illustrated in Fig. 4. To calculate the vertically projected cover, the process in Fig. 4 was repeated using vertical square columns.

TLS parameters
Thus there were two TLS parameters that must be specified, q app and minGap (Table 2). These were found by optimising the TLS generated waveform, convolved with the ALS system pulse, against the raw ALS waveform, less the known modal noise level (step 1 in Fig. 5). The optimisation used the Levenberg-Marquardt method (Garbow et al., 1980). With these parameters, TLS generated waveforms without the ALS system pulse but including attenuation (TLS visible target profile, black line in Fig. 1 (b)) and without attenuation (TLS true target profile, red line in Fig. 1 (a)) were generated. Note that these will give the area projected towards the ALS rather than the surface area. Correction for angular projection and clumping are needed to get surface area (Chen et al., 1997). As this paper is primarily interested in estimating vegetation cover rather than calculating surface area (such as LAI), this was not calculated. The method presented is a necessary first step in the calculation of surface area.

Calibrating ALS with TLS
In this study the TLS generated waveforms were used to calibrate ALS processing to extract target profiles from individual ALS beams; which is the maximum amount of information available to waveform ALS. The ALS processing can be optimised to match the retrieved target profile (Fig. 1 (a)) to the TLS generated target profile (Fig. 4) through step 2 in Fig. 5. As the generation of the TLS waveforms is It should be noted that the ALS parameters in Table 2 are instrument specific rather than site specific. So calculating them for the Leica ALS50-II over one area allows the method to be applied to any Leica ALS50-II datasets. Other ALS instruments may require different values.

Tests performed
The analysis was conducted in three stages, testing each step of processing described above to ensure that each was suitable for measuring accurate 3D vegetation structure. The first assessed the suitability of the TLS data to be used in calibration and validation of ALS data (method described in Section 2.4 and step 1 in Fig. 5). The second examined the calibration procedure, identifying any limitations (method in Section 2.4.2 and step 2 in Fig. 5). The accuracy of the ALS derived voxel product was then assessed and a best combination chosen (method in Section 2.3 and step 3 in Fig. 5). Finally the results were compared to the traditional discrete return and Gaussian decomposed ALS data products.

Step 1: TLS voxelisation
The TLS parameters were found for each individual ALS beam to see whether a global set of parameters could be found. Sensitivity to voxel size was also tested. Whilst Calders et al. (2015) showed that TLS reaches the tree tops and Hancock et al. (2014) showed that TLS point clouds give accurate gap fractions, Hancock et al. (2014)  Optimal denoising parameters along with mean RMSEs between ALS and TLS derived target profiles. Method labels are described in Table 2 and parameter meanings in Table 3 (Béland et al., 2011) or else covered a canopy with much higher scan densities than are possible in a wide area field study (Hosoi and Omasa, 2006). To test the ability of the attenuation correction to overcome any errors at the canopy top, TLS generated waveforms with the system pulse were compared to raw ALS signal in areas of tall, dense foliage.

Step 2: ALS beam level optimisation
A calibration dataset was made up of 317 ALS waveforms from three plots (1, 3 and 4). These were manually selected to cover the range of observed waveform shapes (trees, trees with understorey, short bushes and hard ground). To avoid areas with changes in vegetation structure in the two years between data acquisition (whether growth or human management such as hedge trimming or shrubbery planting), only waveforms with accurate fits of the raw ALS to the TLS generated waveforms with system pulse were used for calibration. An analysis of the impact of changes during the time between measurements in such a heterogeneous and managed area was beyond the scope of this paper. ALS processing parameters were found by optimising these 317 waveforms against the TLS generated target profiles using the Levenberg-Marquardt method (Garbow et al., 1980). The Levenberg-Marquardt method searches for local minima of the error surface. In order to find the global minimum the optimisation was performed with a range of initial parameter values covering the expected bounds, giving 625 separate optimisations for each method. The parameter set with the lowest root mean square error (RMSE) of the ALS to TLS derived foliage profiles was selected.
Errors due to the attenuation correction were assessed. The accuracy was examined for the different ALS waveform shapes to assess how the method performed over different land cover types.

Step 3: Plot scale validation
A voxel map of vertically projected fractional cover of objects in voxels in a 35 m sided cube around each plot centre was calculated separately from the TLS and ALS data for each plot. The TLS tunable parameters (minGap and q app Section 2.4) were calculated for each ALS beam and the mean of all ALS beams intersecting a voxel used. The ALS parameters calculated in Section 2.5.2 were used. The ALS data were limited to 1.5 m horizontal (controlled by footprint density) and 15 cm vertical (controlled by digitisation rate) resolution, therefore analysis was carried out at 1.5 m horizontal and 50 cm vertical resolution. The omission-commission errors were assessed, as given in Eqs. (7) and 8.
where; E om is the omission error E com is the commission error N (tls>0)∧(als=0) is the number of voxels with positive TLS derived cover and zero ALS derived cover N (als>0)∧(tls=0) is the number of voxels with positive ALS derived cover and zero TLS derived cover N (tls>0) is the number of voxels with positive TLS derived cover and N (tls=0) is the number of voxels with zero TLS derived cover The mean TLS derived cover within the N (tls>0)∧(als=0) voxels and the ALS derived cover within the N (als>0)∧(tls=0) voxels were calculated to determine the magnitude of the errors. To test the sensitivity of the final product to the ALS processing parameters, the noise threshold and deconvolution convergence tolerance were varied about the optimum values and the omission-commission errors examined.
The number of TLS, N tlsbeams , and ALS beams, N alsbeams , intersecting each voxel was used as a measure of uncertainty within that voxel, as the greater the number of beams intersecting a voxel, the greater the confidence in the derived canopy cover. In addition the mean RMSE of the TLS generated waveforms, including system pulse, compared to the raw ALS waveforms, D tls−als was recorded per voxel as a measure of uncertainty due to any change in structure between ALS and TLS measurements or violations of the assumptions made in the TLS voxelisation. This is given by: where; N waves is the number of ALS waveforms intersecting a voxel N bins is the number of waveform bins (256 in this case) I tls,n,i is intensity within the ith bin of the nth TLS derived waveform and I als,n,i is intensity within the ith bin of the nth ALS derived waveform

Benchmarking
The agreements between the TLS generated voxel map (used as truth) and the ALS generated voxels, the discrete return ALS and the Gaussian decomposed ALS data were visually compared for plots 1 (woodland with dense understorey) and 8 (oak tree with hedge and lawn) in order to assess the relative information content of the different methods. Vegetation cover maps for different heights above ground were created over 100 km 2 and visually assessed. Fig. 6 shows examples of raw ALS waveforms along with waves generated from TLS point clouds using the method in Section 2.4. This shows that, once attenuation and pulse shape are taken into account, the TLS data are capable of near perfectly recreating the observed ALS waveforms. This was observed in the majority of ALS waveforms close to scan centres, unless there were an obvious change in structure in the two years between measurements. It was found that no single set of TLS parameters was appropriate for all ALS beams across all areas, which is not surprising due to target reflectance heterogeneity. Therefore a separate pair of TLS parameters was calculated for each ALS waveform. For the canopy scale voxel map, the mean TLS parameters for all ALS beams intersecting that voxel were used. No dependence of accuracy on voxel size was found and so unless otherwise stated, 50 cm cubic voxels were used for the TLS data. Fig. 7 shows a raw ALS waveform and TLS generated waveform with system pulse and ALS attenuation (as in Fig. 1 (c)), with and without correcting for TLS attenuation. This shows that without correcting for TLS attenuation (green line), the TLS derived vegetation density at the tree top was underestimated compared to ALS (blue line). Scaling TLS points by the gap fraction up to those points (purple line) corrected this. This was observed for all ALS  waveforms with significant TLS attenuation. This demonstrates that the TLS processing described in Section 2.4 can be used to generate the necessary synthetic waveforms to calibrate and validate ALS data.

ALS optimisation results
Fig . 8 shows the method successfully retrieving the canopy structure from measured waveforms using the optimum parameters across the 317 calibration waveforms, given in Table 4. The RMSE is low and the layers of vegetation have been correctly identified. Fig. 9 shows the difference between attenuation corrected projected area from processed ALS data and "true" projected area from TLS data. This suggests that the attenuation correction will produce estimates of vegetation area with 90% accuracy until 95% of the signal has been attenuated, at which point errors increased to 21%. Therefore understorey is correctly identified under a dense canopy. Fig. 10 shows two examples of hard targets and one diffuse. The hard target identification was successful in Fig. 10 (a) whilst that in Fig. 10 (b) was not correctly identified, leading to a diffuse processed waveform that could be interpreted as short vegetation. This suggests that the correlation threshold used to match returns to the system pulse may have been too stringent. Increasing the correlation threshold would correct the examples in Fig. 10, but with the danger of incorrectly treating diffuse returns as hard targets, such as Fig. 10 (c), where a strong canopy return has been incorrectly identified as a hard target. The next section will test whether this error had a significant impact at the plot scale when voxelising. Table 7 Omission-commission errors for voxelisation using GFps on plot 2 for a range of signal processing parameters. Values are in % of voxels with false positives and negatives, and mean cover of those voxels.

Plot scale analysis
No dependence of errors on the three uncertainty metrics described in Section 2.5.3 (N tlsbeams , N alsbeams and D tls−als ) was found, suggesting that all the data in the 35 m cubes could be used. The six denoising methods in Table 3 were compared to choose the most accurate. The omission-commission errors are shown in Table 5 for plot 2 (woodland) in terms of the percentage of voxels that show false positive or negative cover and the mean cover of those false voxels. The relative differences between the methods at this plot were representative of all plots.
All methods had false positives in 4% to 11% of the voxels, suggesting that either not all the noise was removed or that deconvolution was not removing all the system pulse blurring. The false positive canopy covers were low (12% to 21%) and so would have given small deviations during the optimisation. These low cover false positives could be filtered out after the signal processing or the optimisation could be weighted to more strongly penalise weak false positives.
There were far fewer false negatives (0.33% to 1.7%), with covers of 8% to 20%, except for GVh, which was an outlier with 7% of voxels having false negatives with a cover of 25%. For this reason GVh was rejected. Note that the false negatives are 0.3% of voxels that contain 9% vegetation cover to give a 0.03% total canopy cover error.
A hard noise threshold gave half the number of voxels with false positives as noise tracking, but those voxels had twice the canopy cover, giving the same final error. Noise tracking gave a much smaller fraction of voxels with false negatives with a similar mean cover of those voxels. Examining profiles revealed that the main difference was in the understorey, where noise tracking's ability to detect weak signals was more apparent in plots with dense understorey (plot 1) than those with sparser understorey (plot 2).
From these results, Gold's method with a fixed noise threshold, pre-smoothing and noise tracking (GFps) gave the lowest voxel error, with the second lowest fraction of voxels with false negatives and only a 0.07% greater false positive error than GVnt (the lowest false positive error other than the rejected GVh). The results for GFps over Fig. 11. Transect through ALS and TLS derived voxels (1.5m horizontal resolution by 50cm vertical resolution). For ALS GFps and TLS colour represents vertically projected cover (blue, low, to red, high). ALS colouring has been rescaled to remove bias. For discrete return and Gaussian decomposed ALS colour represent return strength in order to illustrate the amount of information available. all plots are shown in Table 6. The fraction of voxels with false negatives was very low in all plots, suggesting that nearly all understorey is detected. Changes in vegetation structure in the two years between ALS (2012) and TLS (2014) measurements may have contributed to the errors, although the dominance of false positives suggests that this was not the case as we would expect vegetation to grow when no tree felling was observed by visual comparison of ALS and TLS voxel maps. The overall accuracy of voxel covers was 24% (RMSE) with a mean bias of 5.6%, but with no trend with height within the canopy.

Parameter sensitivity
In order to test how sensitive the results were to the denoising parameters, the plot scale analysis was repeated, varying the threshold (thresh) and deconvolution tolerance (deconTol) to see how sensitive the voxel maps were to these. The noise threshold was varied from 14 (1 DN above mean noise, the lowest possible threshold) to 24 (10 standard deviations above mean noise, a very high noise threshold). The deconvolution tolerance was varied by a factor of 10,000 either side of the optimum value. Table 7 shows the omissioncommission errors using the GFps method over plot 2. This shows that there was a negligible difference in errors with varying denoising parameters, and so the method is robust to the parameter choice. All plots showed a similar tolerance to changing parameters, but note that these are instrument rather than site specific parameters and so we would expect all sites to behave similarly. The robustness to threshold is not surprising given the use of threshold insensitive noise tracking (Hancock et al., 2011) and minimum feature widths (minWidth). That varying the deconvolution tolerance by a factor of 10,000 did not significantly affect the omission-commission errors suggests that the solution had already converged to an acceptable limit at the highest tolerance threshold (10.00 × 10 −8 ), but that the RMSE relative to TLS did continue to decrease as the threshold was lowered, and so a lower tolerance was selected as the optimum in Table 3. This suggests that computational expense could be lowered by raising the deconvolution tolerance without significantly affecting the results.

Benchmarking
Transects: Fig. 11 shows a side view of a 2 m wide, 35 m long transect through plots 1 (forest with understorey) and 8 (large tree over a hedge and lawn) for the TLS and ALS derived voxels. Point clouds for discrete return ALS and Gaussian decomposed waveform ALS are shown to illustrate the amount of information available to these methods. The ALS voxel colour scales were adjusted to remove the bias shown in Fig. 9.
Whilst Gaussian decomposition has been shown to give accurate estimates of canopy cover over tens of metres (Armston et al., 2013) and be useful for land cover classification (Mallet et al., 2011), Fig. 11 shows that it gives limited information on structure within the canopy and none on the understorey. The discrete return ALS data gave similarly limited structural information but had the added disadvantage of giving low intensity returns from the tree tops due to their diffuse nature . Thus no useful radiometric information was contained to help calculate canopy cover or land cover class. Given the low amount of information available to these methods, no attempt was made to estimate voxelised canopy cover at high-resolution from them.
In both plots, GFps has retrieved the correct pattern of canopy cover, with gaps within the tree crown and details of understorey evident and areas of flat ground correctly identified. There are some differences, with TLS giving higher covers at the base of the trees and lower estimates at the top of tall shrubs. ALS will not see trunks and so a difference is not surprising. The shrubs in plot 1 were very dense so that disagreement at the top may have been due to high TLS attenuation. Fig. 12 shows vegetation cover maps at different heights above ground over an area with woodland, shrubs and buildings. The ALS voxel map has successfully detected the path through the trees and captured the understorey vegetation elsewhere. Variations in understorey vegetation, not apparent from the 2D aerial photos, are visible in the voxel map. Of particular note is that the vegetation at 1 m above ground is denser along the edge of the woodland than in the centre, which we would expect due to increased light availability. Figs. 11 and 13 show that neither the discrete return or Gaussian decomposed waveform data would be able to identify the difference between paths and understorey under tree canopies or accurately determine understorey density at resolutions finer than a few tens of metres. The lack of information in the centre of Fig. 13 (a) is particularly notable.

Conclusions
A method to produce high resolution voxel maps of vegetation cover (1.5 m by 1.5 m horizontally by 50 cm vertically), including understorey, has been presented. Comparison to the commonly used discrete return and Gaussian decomposition of waveform ALS shows that the new method captures far more detail on within-canopy and understorey structure, even through dense woodland canopies, and can be applied at the landscape scale.
A new method to use terrestrial laser scanner (TLS) data to calibrate and validate the ALS, instrument specific signal processing was presented. Some blurring of the voxel map was apparent,  leading to 10% of voxels containing low density false positives, but less than 0.4% of voxels containing vegetation were not detected, suggesting that nearly all understorey is detected. These false positives could easily be filtered by removing voxels with less than a threshold (e.g. 15%) cover. The overall canopy cover accuracy was 24% (RMSE). The results were insensitive to perturbing the processing parameters, showing that the method is robust and can be applied to large areas. Voxel maps of vegetation cover derived by this method can be used in a wide range of applications such as monitoring understorey health and habitat availability, quantifying viewsheds in complex environments, mapping easy walking routes in forests with closed canopies, as a direct measure of sunlit and shaded leaf area for radiative transfer (Kotchenova et al., 2004), investigating ecosystem services (Grafius et al., 2016) and calculation of explicit clumping factors for calibrating coarser resolution spaceborne sensors.
The use of TLS data to optimise and validate ALS processing can be applied to any instrument, biome or ALS processing method. This optimisation method has been implemented as a C library and is available online (https://bitbucket.org/StevenHancock/voxel_lidar). An inter-comparison of all ALS denoising methods using this tool would allow a thorough comparison and give an insight into their performance.