DUAL-WAVELENGTH TERRESTRIAL LASER SCANNING AS A CALIBRATION TOOL FOR SATELLITE ESTIMATION OF FOREST CANOPY MOISTURE CONTENT

Satellites can estimate forest canopy moisture content at the landscape level, and thus have been widely utilized in forest health monitoring. However, the calibration and validation of the estimation models can be challenging. Collecting a sufficient number of leaf samples from the canopy top layers in sampling plots that match the pixel size of the sensor is needed, which is a cost, time and effort consuming process. Dual-wavelength Terrestrial Laser Scanning (TLS) has been successfully used in estimating canopy moisture content of individual trees in three dimensions (3D) in several recent studies. Such 3D estimates, if produced at the plot level, can be used in the calibration and validation of satellite forest canopy moisture content estimation models. In this study, forest canopy moisture content, quantified as the leaf Equivalent Water Thickness (EWT), was estimated in 3D at the plot level in a mixedspecies deciduous broadleaf forest plot using dual-wavelength TLS intensity data (808 nm near infrared and 1550 nm shortwave infrared wavelengths). The relative error in the EWT estimation was 6%, and the EWT point cloud revealed vertical heterogeneity in the EWT distribution. EWT was 37% higher in the canopy top layers than in the canopy bottom layers. The results obtained in this study showed that dual-wavelength TLS has the potential to be used in operational landscape-scale EWT estimation, and can be a useful tool for the calibration and validation of satellite EWT estimation models. * Corresponding author


INTRODUCTION
Satellite optical sensors have been widely used in measuring forest canopy moisture content, which serves as an early indicator of tree stress and risk of wildfires (Colombo et al., 2008). Satellites can provide measurements at the landscape level and produce time series of the change in forest canopy moisture content (Foley et al., 1998). In addition, free access is available to data from a number of the earth observation satellite sensors, including Landsat, MODIS (Moderate-Resolution Imaging Spectroradiometer), Sentinel-2, and Hyperion, reducing the cost needed for continuous monitoring of forest health.
Typically, canopy moisture content at the landscape level is quantified as the product of leaf Equivalent Water Thickness (EWT) and Leaf Area Index (LAI) (Clevers et al., 2010). EWT is the amount of liquid water in a given leaf area (Danson et al., 1992), while LAI is the total one-sided green leaf area in canopy per unit ground surface area (Jonckheere et al., 2004). EWT can be estimated using vegetation indices that combine the reflectance measured in near and shortwave infrared bands, or by the inversion of Radiative Transfer Models (RTMs) (Ceccato et al., 2001;Yebra et al., 2013). RTMs are more widely used because they are not site-dependent, as they simulate the vegetation spectra using the radiative transfer equation, taking into consideration the leaf and canopy biophysical and biochemical characteristics (Yebra et al., 2008). Similarly, LAI can be estimated using vegetation indices and RTMs, utilizing the reflectance measured in visible and near infrared spectral bands (Eriksson et al., 2006).
One challenge associated with the use of satellites in mapping canopy EWT is collecting large number of leaf samples from the upper canopy layers in sampling plots that match the pixel size of the sensor to calibrate and validate the EWT estimation models. For instance, Trombetti et al. (2008) considered the use of destructive sampling to directly calibrate and validate MODIS estimates of canopy EWT in woodlands, shrublands, and grasslands sites in continental USA to be extremely difficult and costly. It was not possible to collect a sufficient number of leaf samples in a large number of sampling plots to directly scale the EWT measurements to 500 m MODIS pixels, and Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data was used for this purpose. Jurdao et al. (2013) reported a 30% error in the estimation of canopy EWT using MODIS data in deciduous broadleaf and coniferous woodlands in Spain, highlighting the challenge of collecting a sufficient number of leaf samples for calibration and validation. Such difficulty was also highlighted by Quan et al. (2017) when Landsat 8 data were used to estimate canopy EWT in evergreen broadleaf, deciduous broadleaf, and coniferous woodlands in Sichuan province, China. Canopy EWT was overestimated, and this was explained as a result of the number of leaf samples collected to calibrate and validate the models being insufficient to represent all canopy properties in the sampled plots.
There is a need for a fast, non-destructive EWT estimation method that can retrieve canopy EWT in large sampling plots. Terrestrial Laser Scanning (TLS) can serve as such a tool because of the high speed and long range of the commercial, modern instruments, making them suitable for collecting data at the plot level. Also, TLS can easily provide information about the canopy top layers, even if these layers are not accessible from the ground. Recently, two studies showed that the intensity data from near and shortwave infrared TLS instruments, after calibrating the intensity to apparent reflectance, can be combined into a Normalized Difference Index (NDI) and used to estimate canopy EWT in woodlands (Elsherif et al., 2019b;Junttila et al., 2019). However, these studies estimated EWT at the canopy level only for individual trees.
The main aim of this study was to use NDI of 808 nm near infrared and 1550 nm shortwave infrared wavelengths, as utilized in the Leica P20 and P40 TLS instruments respectively, in estimating forest canopy EWT in three dimensions (3D) at the plot level. Canopy LAI was not the focus of this study. However, a number of recent studies have shown the potential of using TLS in measuring LAI at the canopy and plot levels (Antonarakis et al., 2010;Meng et al., 2019;Moorthy et al., 2008;Vincent et al., 2015;Zheng et al., 2013).

TLS instrumentation and calibration
Two commercial, time-of-flight, pulsed TLS systems were used in this study: the Leica P20 (808 nm near infrared wavelength) and the Leica P40 (1550 nm shortwave infrared wavelength). The same instruments have been previously utilized in measuring canopy EWT of individual trees in Wytham Woods, Oxfordshire, UK, and calibration of their intensity data to retrieve apparent reflectance is fully described in Elsherif et al. (2019b). The point clouds from the two instruments, when acquired from the same surveying station, can be accurately aligned, and their intensity data can be combined to calculate NDI on a point-bypoint basis (Elsherif et al., 2018;Elsherif et al., 2019b). This is mainly a result of the similarities between the instruments' chassis, scanning mechanism, and laser beam exit location. However, the slight differences between their beam diameter and divergence (Table 1) required further filtering and point matching for the calculation of NDI, as a P20 point cloud always had more points than a corresponding P40 point cloud.
The point matching and filtering are described in Section 2.2.

Study area and TLS data processing
A mixed-species forest plot (40 m × 40 m) in Wytham Woods, Oxfordshire, UK (51.78 °N, 1.31 °W) was scanned with the two TLS instruments from five scanning positions ( Figure 1). The plot was dominated by Acer pseudoplatanus (sycamore) and Fagus sylvatica (beech) trees. Four Leica black and white registration targets were placed in the plot, mounted on tripods with different heights, to be used in aligning the point clouds.
At each scanning position, full-hemisphere scans (360° × 270°) were conducted with a resolution (point spacing) of 3 mm at 10 m, with a duration of fifteen minutes for each instrument. Scanning the plot began at 11 am, and the total duration of the scanning was approximately two and a half hours. This was followed by the leaf sampling described in Section 2.4.
Leica Cyclone version 9.1 (Leica Geosystems HDS) was used for the point cloud registration, using the registration targets placed in the scans. Then, a 3D nearest neighbour function was applied to match each point in a P40 point cloud to its nearest neighbour in a corresponding P20 point cloud. Afterwards, only the nearest neighbour points were retained in the P20 point clouds, and remaining points were disregarded. The point matching and filtering was carried out in Matlab (The MathWorks Inc., USA, 2016). The intensity data were then calibrated to apparent reflectance, and NDI was calculated on a point-by-point basis, following Equation (1), and the NDI point cloud of the plot was generated.
where P20R = reflectance from the P20 instrument P40R = reflectance from the P40 instrument

EWT estimation model
Typically, leaf samples are needed to calibrate the EWT estimation model. However, the NDI used in this study has previously been reported to be species-independent to an extent, meaning that a general EWT estimation model can be transferred to different sites with no need for a recalibration (Elsherif et al., 2018;Elsherif et al., 2019a;Elsherif et al., 2019b). This was examined in this study by developing a general EWT model, using PROSPECT radiative transfer model simulations and the canopy NDI and EWT values reported in the aforementioned studies. This included nineteen trees from six deciduous broadleaf species, including sycamore, beech, Sorbus intermedia (Swedish whitebeam), Fraxinus excelsior (ash), Quercus robur (oak), and Acer davidii (snake-bark maple). No additional leaf samples for the EWT model calibration were collected from the site, aiming at examining the transferability of the general EWT model.
The PROSPECT-5 radiative transfer model (Feret et al., 2008;Jacquemoud and Baret, 1990) was then used to develop a general EWT model. The model was trained using the EWT and Leaf Mass per Area (LMA) values retrieved from the leaf sampling conducted in the aforementioned studies, calculated following Equations (2) and (3), respectively.
The additional parameters of the model, chlorophyll a and b content, carotenoid content, and brown pigment content, were kept constant at the model defaults (47.7 µg cm -2 , 4.4 µg cm -2 and 0, respectively). Such parameters do not affect the reflectance in the near and shortwave infrared wavelengths, and using generic, fixed values for them in the simulations is sufficient (Ceccato et al., 2001;Gaulton et al., 2013;Zarco-Tejada et al., 2003). The leaf structural coefficient parameter (N), which represents the number and thickness of the leaf mesophyll cell layers, was changed between its minimum and maximum values for dicot leaves (1.5 and 2.5, respectively) (Jacquemoud and Baret, 1990), and simulation were carried out with an interval of 0.1.

Leaf sampling for validation
Leaf samples for validation of the EWT estimates were collected immediately after scanning the plot. The total number of leaf samples collected was 53. Table 2 shows the number of leaf samples collected from each species. The canopy top layers were not accessible from the ground, and thus all leaf samples were collected from the canopy bottom using a tree pruner. As there were more sycamore trees in the plot than beech trees, more sycamore leaf samples were collected, in order to better represent the average EWT. The leaf samples were collected from multiple trees from a specific canopy layer with a thickness of approximately 2 m (layer extending between 4 m and 6 m above ground). The leaf samples were collected from the lower branches of the main canopies, and no leaf samples were collected from the understory vegetation. For each leaf, the fresh weight was measured immediately on collection using a precise scale (0.001g division). the surface area of each leaf sample was measured using Image-J 1.50i software (Schneider et al., 2012) after scanning them with a CanoScan LiDE 110 photo scanner (600 dpi). Next, the leaf samples were dried in an oven for 72 hours at 60 °C and then their dry weight was measured. EWT of each leaf sample was calculated following Equation (2). Figure 2 shows the EWT -NDI relationship at the canopy level for the calibration data, in addition to the EWT models derived from the PROSPECT simulations. It was found that all trees had leaves with N values ranging between approximately 1.5 and 2, meaning that a general EWT model with N value of 1.7 can be sufficient to represent all trees from the different species involved. Equation (4) describes the model. Figure 2. EWT -NDI relationship at the canopy level for the calibration data, in addition to the EWT models derived from PROSPECT simulations (black corresponds to N = 1.5, blue corresponds to N = 1.7, and green corresponds to N = 2).

General EWT estimation model
EWT (g cm -2 ) = 0.0720 × NDI -0.0076 (4) However, Figure 2 also shows that for some individual trees, the general EWT model can introduce more errors in the EWT estimation than species and site-specific models. For instance, the snake-bark maple tree, one beech tree, and two oak trees were better represented by an EWT model corresponding to N = 1.5, while two sycamore trees were better represented by a model corresponding to N = 2. According to Figure 2, applying the EWT model corresponding to N = 1.7 to these trees can cause errors in the EWT estimates at the canopy level ranging between 10% and 14%. At the plot and landscape levels, the error in the EWT estimates is expected to be lower, as the errors of the individual trees will be averaged. However, the species present in the site must have N values between 1.5 and 2.
In general, dicot leaves have N values between 1.5 and 2.5 (Jacquemoud and Baret, 1990), and Elsherif et al. (2019a) showed that only species with very thick leaves (high LMA) had N value > 2, for instance, Ilex aquifolium (holly), which is an evergreen species. This, in addition to the number of leaves and species used in training the model, suggested that the EWT model developed in this study can be considered a general model for deciduous broadleaf species, as such species were found to have N values between 1.5 and 2. The model can be further improved by adding more species and leaf samples to the calibration data.

Canopy level
The general EWT model, described in Equation (4), was applied to the NDI point cloud of the plot. Afterwards, all the points that had EWT less than zero were disregarded from the point cloud, as these points corresponded to woody materials or noise, as discussed in Elsherif et al. (2019b). The resulted 3D EWT point cloud of the plot is shown in Figure 3. As shown in the figure, there remained some points in the point cloud that corresponded to woody materials but had EWT slightly above zero. The majority of these points, especially the ones that were clearly part of the trunks, were visually identified and manually removed.
The canopy layer from which leaf samples for validation were collected was extracted from the point cloud. EWT values of all points in the layer were average, and the average EWT was considered to be the layer's EWT. The estimated EWT was compared to the actual EWT measured from destructive sampling, revealing a relative error of 6% in the EWT estimation. The observed error was considered low, as previous studies that utilized TLS in the estimation of EWT have reported errors in the EWT estimation ranging between 3% and 14% (Elsherif et al., 2018;Elsherif et al., 2019a;Elsherif et al., 2019b;Zhu et al., 2015;Zhu et al., 2017). The low error observed indicated that the general EWT model represented the NDI -EWT relationship of the trees present in the plot, although no leaf samples were collected from them to calibrate the model. This further confirmed the transferability of the general EWT model. However, the species present in the site, sycamore and beech, were both included in the model's training data. The general EWT model still needs to be tested in forest plots with broadleaf deciduous species that were not included in the training data to further examine its transferability. Furthermore, the transferability of the general EWT model to other types of forests (boreal forests and tropical rainforests) still needs to be examined. General EWT models also still need to be developed for coniferous species, as the model introduced in this study was developed for broadleaf species only.

EWT vertical profile
The 3D EWT point cloud of the plot was split into 10 horizontal layers plus the ground layer ( Figure 4). The first layer, extending between the ground level and 4 m above ground, represented the understory vegetation. The remaining nine layers, 2 m in depth each, represented the canopy layers. Vertical heterogeneity in EWT was observed in the plot. This agreed with the findings reported in the few previous studies found in the literature that investigated the vertical distribution of EWT at the canopy and plot levels (Arellano et al., 2017;Elsherif et al., 2018;Elsherif et al., 2019a;Elsherif et al., 2019b;Gara et al., 2018;Liu et al., 2015;Zhu et al., 2017).
Average EWT in the top three canopy layers was 57% higher than EWT of the ground, 42% higher than EWT of the understory, and 37% higher than EWT of the canopy bottom layer. Assuming that the plot corresponded to a pixel or more in a satellite image, the operator can choose which canopy layers to use in calibration and validation of the satellite EWT estimation model. EWT of the canopy top layers was produced, although such layers were not accessible from the ground. Furthermore, as the reflectance obtained from satellite sensors is expected to be dominated by the canopy top layers (Liu et al., 2015), the 3D EWT estimates can be useful in providing information about EWT in the canopy bottom, the understory, and the ground, and also in studying the vertical variation of EWT within the plot.

The potential of dual-wavelength TLS in operational landscape-scale EWT estimation
The use of general EWT models, being independent of species and transferable to different sites, can make dual-wavelength TLS a useful tool for operational, landscape-scale EWT estimation in mixed-species woodlands, with no need for a prior knowledge of the species present in the site. In case of the presence of trees from species that are known to have thick leaves (N > 2), such trees can be auto-detected in the point cloud using their NDI values, as such trees will have significantly higher NDI values than the other species, making them appear as statistical outliers. For instance, NDI values of the trees shown in Figure 2 ranged between 0.21 and 0.32, with mean NDI of 0.26 and standard deviation of 0.044. On the other hand, average NDI of holly trees, as an example of species with thick leaves and high EWT and LMA, was reported to be 0.42 (Elsherif et al., 2019a). Using the mean and standard deviation values, holly trees can be auto-detected in the point cloud by using a threshold of the mean + one standard deviation (the threshold = 0.3). Afterwards, the general EWT model can be applied to all trees in the point cloud, whilst a species-specific model can be applied to the auto-detected holly trees.
TLS has the potential to produce time series of the change in EWT at much higher temporal resolution than spaceborne sensors, as TLS is independent of the sun illumination and cloud coverage. For instance, during drought conditions, the daily change in EWT can be monitored, with the potential to even detect the changes in EWT throughout the day. As the EWT estimates are provided in 3D, the change in EWT vertical profiles can be observed, and trees reaction to drought conditions, and how they redistribute resources and water, can be studied at very high spatial and temporal resolutions.
However, one challenge associated with the use of TLS to estimate EWT at the landscape level is the scan time required to cover large areas of woodlands. For instance, applying this method to generate 3D EWT estimates in a one-hectare forest plot, using 20 m × 20 m scan grid and 3 mm point spacing, would require at least 25 hours of scanning (5 days of scanning, 5 hours per day). This would reduce the temporal resolution of the time series of the change in EWT produced for such large areas, which may also be further reduced because of the wind effects, as scans cannot be conducted in windy conditions. The scanning time can be reduced by increasing the point spacing (reducing the scan resolution) or increasing the scan grid size. However, at the landscape level, TLS may only be suitable for detecting the weekly or monthly changes in EWT, depending on the size of the area to be scanned. Scanning only a set of plots that represent the different species in the site can be more useful, as these plots can be scanned daily and the rapid changes in EWT can be monitored. This can also reduce the data collection campaign's cost, as the personnel and their instruments need to move between limited number of scanning locations. In general, cost-wise, TLS systems are becoming more and more affordable, faster, and easier to use, reducing their operational costs, and making them a potential tool for forest health monitoring.
Another factor that needs to be taken in consideration while applying this approach at the landscape level is the effect of the variation of the canopy LAI in heterogeneous sites. Although the use of the NDI can be sufficient to reduce the effects of incidence angle and leaf internal structure, canopy LAI can influence the accuracy of the EWT estimation using vegetation indices. Thus, LAI measurements are needed, and canopy EWT should be quantified as the product of EWT and LAI. Furthermore, the effects of occlusion and forest density on the proposed approach still need to be further investigated. There existed canopy gaps and areas of low-density canopy cover in the study area involved in this research. Setting the scanning positions in locations corresponding to such areas reduced the occlusion effects, thus large number of laser beam returns were acquired from the canopy top layers. In denser forests, the understory and lower canopy layers can reduce the number of laser beams reaching the canopy top layers, which may reduce the accuracy of this approach in quantifying EWT in such layers. In case of inaccessible forest sites, this method cannot be used. One proposed solution for very dense or inaccessible forest sites is applying similar approach using airborne LiDAR data. Recently, the Optech Titan (Teledyne Optech) multispectral airborne LiDAR system was launched, being the first multispectral airborne LiDAR system commercially available. The system utilizes three laser wavelengths: 1550 nm shortwave infrared, 1064 nm near infrared, and 532 nm visible. The NDI of the 1550 nm and 1064 nm wavelengths can be linked to EWT, thus the system has the potential to provide 3D EWT estimates at the landscape level. However, this process will mainly depend on the accuracy of aligning the point clouds from the two different wavelengths, which still needs to be examined.

CONCLUSIONS
In this study, the NDI of 808 nm near infrared and 1550 nm shortwave infrared wavelengths, utilized in the Leica P20 and P40 commercial TLS instruments respectively, was used to estimate forest canopy EWT at the plot level. A mixed-species (sycamore and beech) deciduous broadleaf forest plot was scanned, and the PROSPECT RTM was used to develop a general EWT estimation model. The model was used to estimate EWT in the plot, and validating the estimates using leaf destructive sampling revealed a relative error of 6% in the EWT estimation. The low error suggested that the developed general EWT model was transferable to different sites, but more study areas are still needed to further test its transferability. The plot 3D EWT point cloud revealed vertical heterogeneity in the canopy EWT distribution. EWT in the canopy top layers was higher than that in the canopy bottom layers, understory, and ground by 37%, 42%, and 57%, respectively. EWT of the canopy top layers can be easily extracted from the point cloud and used in calibrating and validating satellite EWT estimation models.
This study showed that dual-wavelength TLS can estimate canopy EWT at the plot level with low error. Combining the proposed approach with an automatic leaf-wood separation algorithm can further improve its applicability in larger areas. Furthermore, combining the EWT estimates with canopy LAI measurements can lead to the use of dual-wavelength TLS in canopy EWT estimation at the landscape level in heterogeneous sites. It is worth mentioning that this was the first study to use dual-wavelength TLS in estimating canopy EWT at the plot level, and until the method is transferred to and tested in more sites with different species, the method should be considered a proof-of-concept.