Tropical Forests of Réunion Island Classified from Airborne Full-Waveform LiDAR Measurements

From an unprecedented experiment using airborne measurements performed over the rich forests of Réunion Island, this paper aims to present a methodology for the classification of diverse tropical forest biomes as retrieved from vertical profiles measured using a full-waveform LiDAR. This objective is met through the retrieval of both the canopy height and the Leaf Area Index (LAI), obtained as an integral of the foliage profile. The campaign involved sites ranging from coastal to rain forest, including tropical montane cloud forest, as found on the Bélouve plateau. The mean values of estimated LAI retrieved from the apparent foliage profile are between ~5 and 8 m2/m2, and the mean canopy height values are ~15 m for both tropical montane cloud and rain forests. Good agreement is found between LiDARand MODIS-derived LAI for moderate LAI (~5 m2/m2), but the LAI retrieved from LiDAR is larger than MODIS on thick rain forest sites (~8 against ~6 m2/m2 from MODIS). Regarding the characterization of tropical forest biomes, we show that the rain and montane tropical forests can be well distinguished from planted forests by the use of the parameters directly retrieved from LiDAR measurements.


Introduction
Tropical forest areas are difficult to monitor/to classify using either remote sensing or in situ approaches, because of their tremendous heterogeneity and complex structure.The fundamental challenge is thus to acquire information about the forest vegetation structure given the fact that forest vegetation limits the ability to acquire information.Forest horizontal patterns are accessible using passive multispectral sensors [1,2] and hyper-spectral sensors [3][4][5], but these sensors are not adequate to penetrate beyond the upper canopy layer [6].Active remote sensing instruments, including LiDAR and radar, have more of a chance to peer through the forest canopy down to the ground level [7].Radar yields volumetric scattering information in addition to surface scattering observations, but the retrieval of the vegetation vertical structure is not direct, unlike with LiDAR.
Ground-based LiDAR systems, either terrestrial or portable, can accurately estimate canopy structural parameters [8][9][10][11][12]; however, covering large areas with such systems is impractical.Airborne/spaceborne LiDAR technology has been used to rapidly describe forest structure over large areas; whereas the observations of optical remote sensing is often limited by cloud in tropical areas.Several airborne discrete return LiDAR datasets have been acquired over tropical forests and have been successfully used to derive structural characteristics, such as canopy height, canopy cover and aboveground biomass [13][14][15].A full description of the forest vertical structure (including canopy top, tree crown base height and understory structures) has also been obtained by airborne full-waveform LiDAR, both with infrared wavelengths [16,17] and ultraviolet wavelengths [18,19].Recently, the airborne demonstration instrument called the Laser Vegetation Imaging Sensor (LVIS) [20] has shown that a full-waveform infrared LiDAR with a large footprint can reliably extract the vertical structure and Leaf Area Index (LAI) of a tropical rainforest (Costa Rica [21]) as well as a mid-latitude forest (California [22]), even with a dense canopy cover.The LVIS team however acknowledges the need of a broader dataset on multiple tropical biomes to confirm these findings and compare the extracted features.
The overarching goal of this paper is to report on a methodology of classification using a full-waveform ultraviolet airborne LiDAR with a large footprint, of varied tropical forest types on Réunion Island, which is a rich diversity of tropical ecosystems listed as World Heritage by UNESCO.The classification, obtained using LiDAR-derived canopy height and LAI, has distinguished native forest from plantations/exotic forests.The study sites and data collection methods will be described in Section 2, where the main steps of LiDAR processing for forest studies and the sampling strategy will also be presented.The retrieved forest structural and ecological properties will be analyzed in Section 3, along with comparisons to ground-based census and spaceborne observations.The classification of tropical forest sites will be also presented in this section.

Study Sites
Réunion Island is a French overseas department located in the Indian Ocean (20 ˝06'52"S, 55 ˝31'57"E; Figure 1).It is a small (2512 km 2 ) tropical volcanic island, which reaches 3070 m in altitude at its highest point (Piton des Neiges).In spite of the transformation of its habitats [23], the island still shelters 100,000 ha of native ecosystems (included in a national park) and is home to the last remnants of intact tropical forests in the Mascarenes archipelago (Réunion, Mauritius, Rodrigues).
Remote Sens.2016, 8, 43 2 Ground-based LiDAR systems, either terrestrial or portable, can accurately estimate canopy structural parameters [8][9][10][11][12]; however, covering large areas with such systems is impractical.Airborne/spaceborne LiDAR technology has been used to rapidly describe forest structure over large areas; whereas the observations of optical remote sensing is often limited by cloud in tropical areas.Several airborne discrete return LiDAR datasets have been acquired over tropical forests and have been successfully used to derive structural characteristics, such as canopy height, canopy cover and aboveground biomass [13][14][15].A full description of the forest vertical structure (including canopy top, tree crown base height and understory structures) has also been obtained by airborne full-waveform LiDAR, both with infrared wavelengths [16,17] and ultraviolet wavelengths [18,19].Recently, the airborne demonstration instrument called the Laser Vegetation Imaging Sensor (LVIS) [20] has shown that a full-waveform infrared LiDAR with a large footprint can reliably extract the vertical structure and Leaf Area Index (LAI) of a tropical rainforest (Costa Rica [21]) as well as a mid-latitude forest (California [22]), even with a dense canopy cover.The LVIS team however acknowledges the need of a broader dataset on multiple tropical biomes to confirm these findings and compare the extracted features.
The overarching goal of this paper is to report on a methodology of classification using a fullwaveform ultraviolet airborne LiDAR with a large footprint, of varied tropical forest types on Réunion Island, which is a rich diversity of tropical ecosystems listed as World Heritage by UNESCO.The classification, obtained using LiDAR-derived canopy height and LAI, has distinguished native forest from plantations/exotic forests.The study sites and data collection methods will be described in Section 2, where the main steps of LiDAR processing for forest studies and the sampling strategy will also be presented.The retrieved forest structural and ecological properties will be analyzed in Section 3, along with comparisons to ground-based census and spaceborne observations.The classification of tropical forest sites will be also presented in this section

Study Sites
Réunion Island is a French overseas department located in the Indian Ocean (20°06'52"S, 55°31'57"E; Figure 1).It is a small (2512 km 2 ) tropical volcanic island, which reaches 3070 m in altitude at its highest point (Piton des Neiges).In spite of the transformation of its habitats [23], the island still shelters 100,000 ha of native ecosystems (included in a national park) and is home to the last remnants of intact tropical forests in the Mascarenes archipelago (Réunion, Mauritius, Rodrigues).Seven plots on Réunion Island were used for forest sampling in our study (Figure 1; Table 1).The coastal test site (CT) has only exotic vegetation.The Cryptomeria (CM) and Tamarind (TM) plots, located on Mount Maïdo, and the Bélouve (BF) site, located on a central plateau, are tropical montane cloud forests, which in particular still cover large areas of Réunion Island (60,000 ha), extending from 800 to 1900 m above mean-sea-level (amsl) on the windward side and from 1100 to 2000 m amsl on the leeward side of the island.This dense forest within cultivated forests on moderate slopes is very similar to Acacia koa forests in Hawaii.The CM is a monoculture of Cryptomeria japonica, an introduced species in the Taxodiaceae family.These trees commonly reach 20 to 25 m in height on the site; they produce a dense canopy, under which light is very scarce.Most acacia stands, like the TM site, are composed of secondary forest and display a monospecific acacia (highland tamarind) canopy with shrubby vegetation in the understory, of which the structure can vary with the intensity of human activities (stock farming in particular).The Mare-Longue sites (three plots dubbed ML-150, -250 and -550, according to their altitude) are located in the National Park of Réunion Island in the former Mare-Longue nature reserve, which shelters the last remnant of lowland tropical rainforest in the Mascarene Islands with around 4000 mm of yearly rainfall.This lowland forest grows on a non-altered basaltic pahoehoe lava flow dated between four and six centuries old [24].This forest displays the greatest tree species diversity on Réunion Island with an average richness of 40 tree species per hectare [25].Whereas average tree height remains very low (15 to 20 m), the stem density exceeds 1000 trees/ha (diameter at breast height >10 cm).The most abundant tree species in the sampled plots is Labourdonnaisia calophylloides (Sapotaceae), endemic from the Mascarene Islands.

Data Collection
The following estimation of forest parameters on Réunion Island was performed in May 2014, combining mainly airborne LiDAR measurements with in situ approaches.

Airborne LiDAR and Instrumentation
The LiDAR system used during the campaign is the Ultraviolet LiDAR for Canopy Experiment (ULICE; [19]) developed at Laboratoire des Sciences du Climat et de l'Environnement (LSCE) with the support of CNES (Centre National d'Etudes Spatiales).It was integrated into an autonomous payload flown on an ultra-light aircraft (ULA) shown in Figure 2. The ULICE system characteristics and the airborne payload are given in Table 2.The ultraviolet domain is well suited both for eye safety at low carrier altitude and for precise retrievals over dense forests with little distortion due to multiple scattering (Shang and Chazette, 2014).As recommended by several authors [20,[26][27][28][29], a large LiDAR footprint is preferred so that the laser can better penetrate the dense tropical forests.Shang and Chazette [26] estimated an optimal laser footprint diameter around 20 m for dense forests, whereas Riaño et al. [30] found that the Leaf Area Index (LAI) was better estimated using laser footprints between 7.5 and 12.5 m.The ULICE system was thus modified to obtain a large and controllable sounding area on the ground (approximately a 1-m to 10-m footprint diameter for a flight altitude of ~350 m above the ground level (agl)).The effect of LiDAR footprint size will be studied in Section 2.4.
Remote Sens.2016, 8, 43 4 between 7.5 and 12.5 m.The ULICE system was thus modified to obtain a large and controllable sounding area on the ground (approximately a 1-m to 10-m footprint diameter for a flight altitude of ~350 m above the ground level (agl)).The effect of LiDAR footprint size will be studied in Section 2.4.An ancillary positioning instrument (inclinometer and GPS), an MTi-G system by XSense, is also onboard the ULA.It provides the horizontal geolocation of the ULA with 5-m accuracy and the direction of the laser beam with 0.7° accuracy (i.e., 3.6 m at the ground for a flight altitude of 300 m•agl).With such uncertainties, the study performed at Réunion Island should be statistical, because we cannot distinguish one tree from another.A Tetracam ADC (Agricultural Digital Camera) air camera is also onboard to map the photosynthesis activity index (Normalized Difference Vegetation Index (NDVI)) over the forest canopy to check the scene heterogeneity.However, its   An ancillary positioning instrument (inclinometer and GPS), an MTi-G system by XSense, is also onboard the ULA.It provides the horizontal geolocation of the ULA with 5-m accuracy and the direction of the laser beam with 0.7 ˝accuracy (i.e., 3.6 m at the ground for a flight altitude of 300 m¨agl).With such uncertainties, the study performed at Réunion Island should be statistical, because we cannot distinguish one tree from another.A Tetracam ADC (Agricultural Digital Camera) air camera is also onboard to map the photosynthesis activity index (Normalized Difference Vegetation Index (NDVI)) over the forest canopy to check the scene heterogeneity.However, its images showed high NDVI (>0.65) over all of the observed sites, making it rather irrelevant for the validation of other ecological parameters, such as LAI [31].

Field Data Collection
During two months around the airborne measurements, four representative sub-plots of ~0.2 ha were set up in four forest sites (CM, TM, ML-150 and ML-250 sites, as shown in Table 1), where in situ measurements were performed.Within each sub-plot, all trees with diameter at breast height (DBH) >7 cm were identified.The tree top height (TTH) and the DBH were measured using a dendrometer and forestry measuring tapes, respectively.For trees with multiple stems, each significant stem was recorded individually.The uncertainties on the TTH from in situ measurements in a dense tropical forest have been evaluated during the experiment in the order of ˘4 m (several measurements on the same tree with different operators with a dendrometer).This is due to the difficulty in identifying the tree top among other branch extremities.

Other Data Collections
Digital terrain models (DTM) of 500-m or 5-m resolution for the whole Réunion Island were provided by the Parc National de la Réunion (J.-C.Notter, personal communication).The topography of Réunion Island is given in Figure 1, using the DTM-500 m.The slope of each sampled site was evaluated using the DTM-5 m.
MODIS (Moderate Resolution Imaging Spectroradiometer) Level 3 land products were compared with LiDAR observations.The 8-day LAI products derived from MODIS onboard Terra and Aqua are considered [31].On 8-day syntheses from May to August 2014, LAI values retrieved on 1-km pixels were averaged after screening for cloud contamination (i.e., values below the median were removed).

LiDAR Data Processing
From airborne LiDAR measurements, the forest structural and optical parameters are estimated; they are in turn used to evaluate ecological parameters, such as LAI.In this study, three key parameters are estimated from the LiDAR backscatter profiles to characterize the sampled forest sites: the canopy height (CH), the vertical profile of apparent foliage (F app ), which informs both the canopy density and the vertical distribution of leaf biomass along the profile, and the LAI, which is linked to the integral of the latter parameter.

Forest Structural Parameter: Canopy Height
The canopy height (CH) parameter is assessed from the full-waveform LiDAR profile using the threshold approach documented in Chazette et al. [32] and applied to forest detection by Cuesta et al. [18] and Shang and Chazette [19].CH is estimated as the distance between the first return, at the upper surface of the vegetation, and the last return, which is normally the ground echo.For dense forests, the laser beam cannot always penetrate the leaves and reach the ground, so the last return of the backscattered LiDAR signal is not necessarily the ground echo.Nevertheless, as frequent measurements were performed (33 Hz), allowing some overlap, the ground level can be correctly located thanks to time-integrating signal processing in almost all cases and be the reference to estimate the CH.An example is given in Figure 3.Note that a parasitic echo (undershot) can be observed beneath the ground echo, which is due to the rebounding, non-linear response of the detector to the strong pulse returned by the ground.The standard deviation of LiDAR-derived CH was assessed to be ~1.5 m when only considering measurement noise and signal processing errors.Shang and Chazette [26] assessed that LiDAR signal distortion due to the surface slope can lead to a relative CH uncertainty of ~5% for a slope of 30 ˝(see Table 1) and a 10-m footprint, as is the case in our present study.As a result, the standard deviation of our LiDAR-derived CH should be of the order of ~2 m.We do not consider geolocation errors in this statistical study.) for a section of flight over the tropical forest of Bélouve (BF).The difference between the ranges of these two points yields the canopy height (CH).Note that the y-axis is not the ground elevation, but the distance from the emitter.ULA, ultra-light aircraft.

Forest Optical Parameters
The range-corrected backscattered airborne LiDAR signal Sv [33,34], taken at a height above ground level (agl) h inside the forest cover, can be expressed by the LiDAR equation [35]: where K is the instrumental constant and Ta is the atmospheric transmission.The backscatter to extinction ratio BER is a classical parameter used in LiDAR analyses [26,36], which characterizes the probability that an intercepted photon would be backscattered by a scattering layer; it is assumed to be constant for all canopy levels in this study.The canopy extinction coefficient αcanopy (h) is defined as the sum of the absorption and scattering coefficients in the canopy.The forest optical thickness (FOT) is defined only in the forest layer between the considered height (h) and the canopy height (CH) and is given by: Following the method proposed by Ni-Meister et al. [37], we define the transmittance height profile (THP) by taking a ratio of the energy from canopy returns to the total energy, which characterizes the amount of skylight intercepted by vegetation at a given level [38]: with: where ρv and ρg are the canopy and ground reflectance, respectively.The integrated range-corrected canopy return Rv(h) (respectively ground return Rg) is defined as the integral of the LiDAR signal from the canopy top CH to height level h (respectively in the equivalent width of the ground echo ∆hGE, ΔhGE ~4 m): The difference between the ranges of these two points yields the canopy height (CH).Note that the y-axis is not the ground elevation, but the distance from the emitter.ULA, ultra-light aircraft.

Forest Optical Parameters
The range-corrected backscattered airborne LiDAR signal S v [33,34], taken at a height above ground level (agl) h inside the forest cover, can be expressed by the LiDAR equation [35]: where K is the instrumental constant and T a is the atmospheric transmission.The backscatter to extinction ratio BER is a classical parameter used in LiDAR analyses [26,36], which characterizes the probability that an intercepted photon would be backscattered by a scattering layer; it is assumed to be constant for all canopy levels in this study.The canopy extinction coefficient α canopy (h) is defined as the sum of the absorption and scattering coefficients in the canopy.The forest optical thickness (FOT) is defined only in the forest layer between the considered height (h) and the canopy height (CH) and is given by: Following the method proposed by Ni-Meister et al. [37], we define the transmittance height profile (THP) by taking a ratio of the energy from canopy returns to the total energy, which characterizes the amount of skylight intercepted by vegetation at a given level [38]: with: where ρ v and ρ g are the canopy and ground reflectance, respectively.The integrated range-corrected canopy return R v (h) (respectively ground return R g ) is defined as the integral of the LiDAR signal from the canopy top CH to height level h (respectively in the equivalent width of the ground echo ∆h GE , ∆h GE ~4 m): where the normalized ground echo g(h) is modelled as a Gaussian function [39] and can be calibrated by using the returned laser pulse at nadir over a flat surface.Thus, as BER is equal to ρ v π , THP can be expressed as a function of FOT: The ε parameter is usually estimated using a known ratio of canopy and ground reflectance, which was estimated around 2.5 [21] or 2 [38,40] at a 1064 nm wavelength.However, the reflectance values are not available in our study area.Nevertheless, the reflectance ratio can be estimated using only LiDAR measurements as follows.Assuming FOT(0) >> 1 (i.e., ε « 1), which is realistic for thick tropical forests, as highlighted by Shang and Chazette [26], an initial FOT estimator can be evaluated from Equations ( 3) to (6), and is given for h > 0 by: This leads to a second assessment: which is made after a correction of the first order using the ε parameter explained as: with h 0 chosen inside the undergrowth layer of the forest (2 to 3 m above the ground level).
Meanwhile, the reflectance ratio (ρ v /ρ g ) can be determined using the LiDAR signal by: A final estimate is obtained after correcting for bias towards high values, which can be assessed very reliably using a simulator of the LiDAR measurements taking inversed profiles for each sampling site as an input.Such an algorithm converges within a relative uncertainty of ~20% after corrections of the bias.
This iterative approach was chosen as an alternative to the one based on ground echo normalization of forest transmittance as proposed by Ni-Meister et al. [37].The reliance of the latter on an accurate estimate of the ground to vegetation reflectance ratio is incompatible with the very diverse and variable forest grounds found at Réunion Island (leaves and debris, soil, lava) and leads to important errors on the retrieved ecological parameters.This inversion process has been applied to each suitable LiDAR profile acquired during the flights above the tropical forests of Réunion Island, in order to characterize the various tropical forest sites.

Forest Ecological Parameters
Several studies have shown that LiDAR is a powerful instrument to retrieve the LAI [41,42].The LAI can be derived from LiDAR measurements by (e.g., [21]): with the "apparent foliage profile" (F app ), which can be identified as the vertical profile of vertical projections of foliage elements, defined by the following equation as in Ni-Meister et al. [37]: Comparing Equations ( 2) and ( 8), we find that the F app is actually the canopy extinction coefficient α canopy in the classical LiDAR equation.We will use F app in the following.
All of this mathematical development can be used whatever the wavelength, but for visible or infrared wavelengths, it may be necessary to consider the multiple scattering effect due to leaves and branches [26].The multiple scattering enhances the backscatter LiDAR signal and makes the LiDAR signal distorted.It can be taken into account by using a multiple scattering parameter as in Platt [43], Berthier et al. [44] and Shang and Chazette [26].
The LAI is an integration of F app .The random orientation of foliage [37] is introduced as G = 0.5.Clumping coefficient C has been assessed by Chen et al. [45] around 1.58 using the bidirectional reflectances derived from the Polarization and Directionality of the Earth Reflectance (POLDER, [46]) instrument onboard the Advanced Earth Observing Satellite (ADEOS), as well as by He et al. [47] using the MODIS BRDF (Bidirectional Reflectance Distribution Function) product, yielding C = 1.54 ˘0.05 over most tropical forests and on Réunion Island.Note that this clumping coefficient is the reciprocal of the clumping index defined as the ratio of the effective LAI and the true LAI in some literature [45,47].The LAI here calculated is a crude estimate of the true LAI, because it takes into account the contributions of branches and trunks.Tang et al. [21] considered that the majority of the backscattered energy measured (93%) was due to leaves, whereas only 7% came from the rest of the tree.Nevertheless, such a value is not justified in their article, and we do not have the capacity to verify it in the current study.

Sampling Strategy
Spatial sampling is a key parameter when using airborne LiDAR to characterize forest plots.Figure 4 gives an overview of the study sites and examples of LiDAR measurements performed continuously from the ULA.Note that the availability of ground echoes with a good signal-to-noise ratio (SNR) is highly variable with the forest site, depending on the vegetation density or the existence of gaps due to dead trees.We therefore had to adapt the sampling approach or find a representative sample for all of the sites.
Weather conditions with almost daily cloud formation over the tropical forest sites, coupled with trade winds or recirculation currents often exceeding 10 m¨s ´1, forced us to revise our initial sampling strategy of forest sites.It was not realistic to expect ground traces to be sufficiently numerous and close to each other to reproduce a 3D vision of forest structures.Nevertheless, we had to check that our samples remained representative.
The horizontal sampling grid during our airborne LiDAR experiment is defined by three independent parameters: the diameter of the laser footprint (d), the sampling along the ground track of the ULA (∆X) and the sampling along the perpendicular to the ground track (∆Y).The last one must be more specifically defined, as it is based on successive passes of the ULA above the same forest site.
Laser footprint: In order to evaluate the influence of the laser footprint size (d) on the retrieval of the ground echo for a flat surface, three specific flights have been conducted at the same flight altitude over the Tamarind site (TM) with laser footprints at the ground level of 4, 10 and 20 m, respectively.We note no significant difference in the statistical distributions of tree structures between these experiments.Indeed, the treefall gaps help to identify the ground echo when using a sufficiently large laser footprint.A laser footprint of ~10 m associated with a ~350-m flight altitude was therefore considered for the entire sampling campaign.This footprint size is comparable to the overall span of dominant trees, and the ground echoes could be perceived from the optically thinner areas between the trees at each laser shot.Such a value is adequate for a correct assessment of the LAI, as shown by Riaño et al. [30], who found that LAI was better estimated using laser footprints between 7.5 and 12.5 m.Along-track and cross-track samplings: The influence of the along-track and cross-track sampling distances (ΔX and ΔY) on the horizontal sampling was evaluated, in order to ensure that the sampling was sufficient to accurately retrieve the canopy structure, i.e., the correct CH distribution.After the accumulation of a sufficient number of samples, ΔX and ΔY were found to be log-normally distributed.On the long flight (~40 km) performed over the Bélouve site (BF), we assessed the mean value and standard deviation of ΔX to be 0.9 ± 0.5 m.Such a value is fully suitable for sampling dense tropical forest from an airborne LiDAR.The sampling distance ΔY between each ground track was around 100-times (respectively 10-times) larger than ΔX in the BF site (respectively other sites), because a flight pattern including too many overpasses over the forest plot is not feasible.To assess the effect of reducing the horizontal sampling frequency, the CH distribution is computed with artificially increased ΔX values ranging from 1 to 100 m, i.e., using only a fraction of the LiDAR Along-track and cross-track samplings: The influence of the along-track and cross-track sampling distances (∆X and ∆Y) on the horizontal sampling was evaluated, in order to ensure that the sampling was sufficient to accurately retrieve the canopy structure, i.e., the correct CH distribution.After the accumulation of a sufficient number of samples, ∆X and ∆Y were found to be log-normally distributed.On the long flight (~40 km) performed over the Bélouve site (BF), we assessed the mean value and standard deviation of ∆X to be 0.9 ˘0.5 m.Such a value is fully suitable for sampling dense tropical forest from an airborne LiDAR.The sampling distance ∆Y between each ground track was around 100-times (respectively 10-times) larger than ∆X in the BF site (respectively other sites), because a flight pattern including too many overpasses over the forest plot is not feasible.To assess the effect of reducing the horizontal sampling frequency, the CH distribution is computed with artificially increased ∆X values ranging from 1 to 100 m, i.e., using only a fraction of the LiDAR shots.The results, presented in Figure 5, show that the CH distribution at ∆X = 10 m is similar to the reference CH distribution, which is the one at the native resolution of ∆X « 1 m.Moreover, in Figure 5d, CH histograms appear similar for values of ∆X between 1 and 100 m.Hence, the distance ∆Y, which is most difficult to keep during the flights, does not significantly affect the statistical studies performed on the different tropical forest sites.These results on the horizontal sampling also demonstrate the strong homogeneity of the dense tropical cloud forest of Bélouve.Note that similar results are obtained for the tropical rain forest of Mare-Longue.

Results and Discussion on Retrieved Tropical Forest Parameters
Each parameter estimated from the airborne LiDAR measurements over the tropical forests of Réunion Island will be analyzed and discussed in this section.The results derived from LiDAR measurements on the various sites will also be compared.

Results and Discussion on Retrieved Tropical Forest Parameters
Each parameter estimated from the airborne LiDAR measurements over the tropical forests of Réunion Island will be analyzed and discussed in this section.The results derived from LiDAR measurements on the various sites will also be compared.

LiDAR-Derived Canopy Height
The number of available laser profiles on each sampled site is closely related to weather conditions during the experiment (Table 3).The two last sites of the Mare-Longue area are likely insufficiently characterized because the number of samplings is not enough for a reliable statistic.Nevertheless they are also taken into account in Table 3 where the mean, median and maximum values of retrieved CH are given, together with the standard deviation around the mean value, which represents so-called canopy rugosity.The statistic has been established on ~90% of the LiDAR profiles, when the ground echoes can be well located.The largest CH values are in the same range on all sites (~30 m), except for ML-250 and -550.The mean and median values are similar, which indicates that there is no significant bias due to outliers in the statistics.The standard deviation is larger than 5 m for the coastal site (CT), i.e., large canopy rugosity, pointing out larger differences in terms of tree maturity on these sites.On the Cryptomeria site (CM), LiDAR observations also show large canopy rugosity.That is because the sub-plot of interest, which has uniform tree height (~22 m), is surrounded with low vegetation.As expected, for the Mare-Longue sites, we observe that CH decreases when altitude increases.It is less noticeable elsewhere, because tree species vary between plots.Table 3. Statistics (mean, median and standard deviation (SD)) for both the canopy height (CH) and the assessment of the LAI on the 7 sites: coastal (CT), Tamarind (TM), Cryptomeria (CM), Bélouve (BF) and Mare-Longue (ML-150, -250 and -550) sites.Profiles with CH < 5 m are not considered.Bold characters highlight the 4 forest sites where in situ measurements are available (see Table 1).The maximal CH derived from the LiDAR is also indicated.

Comparison with in Situ Measurements
In the four sampled sub-plots, the LiDAR-derived canopy heights (CH) were compared to the tree top heights (TTH) from in situ measurements.Their statistical moments and normalized distributions are given in Table 4 and in Figure 7, respectively.Understandably, for in situ measurements, the sampling areas were reduced (~0.2ha), but the sub-plots were chosen so as to be as representative as possible of the extended sites.The agreement is quite good between the two distributions for the Cryptomeria site (CM).This is not the case for the others.With a LiDAR footprint diameter of 10 m, numerous smaller trees are hidden by the higher trees; the LiDAR-derived CH distribution is then biased toward the higher trees, whereas ground measurements can underestimate actual TTH as the tree top may not always be well identified because of the complex canopy.It is thus necessary to consider the apparent foliage profile to better identify the underlying trees.
Table 4. Statistics (mean, median, maximal values and standard deviation (SD)) for both the tree top height (TTH) as measured from the in situ census and the canopy height (CH) retrieved by LiDAR measurements in 4 sub-plots of ~0.2 ha: Tamarind (TM), Cryptomeria (CM) and Mare-Longue (ML-150, -250) sites.The values are given for sub-plots well identified in the main sites (Table 1).All The locations of the sub-plots where in situ measurements were performed are highlighted using thick black lines.

Comparison with in Situ Measurements
In the four sampled sub-plots, the LiDAR-derived canopy heights (CH) were compared to the tree top heights (TTH) from in situ measurements.Their statistical moments and normalized distributions are given in Table 4 and in Figure 7, respectively.Understandably, for in situ measurements, the sampling areas were reduced (~0.2ha), but the sub-plots were chosen so as to be as representative as possible of the extended sites.The agreement is quite good between the two distributions for the Cryptomeria site (CM).This is not the case for the others.With a LiDAR footprint diameter of 10 m, numerous smaller trees are hidden by the higher trees; the LiDAR-derived CH distribution is then biased toward the higher trees, whereas ground measurements can underestimate actual TTH as the tree top may not always be well identified because of the complex canopy.It is thus necessary to consider the apparent foliage profile to better identify the underlying trees.Table 4. Statistics (mean, median, maximal values and standard deviation (SD)) for both the tree top height (TTH) as measured from the in situ census and the canopy height (CH) retrieved by LiDAR measurements in 4 sub-plots of ~0.2 ha: Tamarind (TM), Cryptomeria (CM) and Mare-Longue (ML-150, -250) sites.The values are given for sub-plots well identified in the main sites (Table 1).All trees with diameters at breast height higher than 7 cm are considered in in situ measurements.

TM CM ML-150 ML-250
Location 21 ˝41 25"S, 55 ˝21 1 44"E 21 ˝41 3"S, 55 ˝20 1 13"E 21 ˝21 1 29"S, 55 ˝44 1 43"E  In Figure 6, the Tamarind (TM) and Cryptomeria (CM) sub-plots are highlighted.We notice that the TM sub-plot has scarcer vegetation than its surroundings, whereas the CM sub-plot is denser than In Figure 6, the Tamarind (TM) and Cryptomeria (CM) sub-plots are highlighted.We notice that the TM sub-plot has scarcer vegetation than its surroundings, whereas the CM sub-plot is denser than its surroundings.Images of the ADC-air vegetation camera obtained over Mare-Longue (ML) and Bélouve (BF) present on the contrary a good homogeneity, as well as high NDVI (>0.8), typical of primary tropical forest growing on regular slopes, which is coherent with the results previously discussed in Figure 5.We also found out that LiDAR-derived CHs are comparable between Tables 3  and 4 for the TM and ML-150 sites.Note that for BF and ML-500, it was difficult access the site and almost impossible to identify the top of a tree from neighbors.Consequently, in situ measurements have not been considered as valid for these two sites.For the CM site where trees have about the same maturity, the differences between CH derived from LiDAR and TTH derived from in situ measurements are ~3 m, which is included in the standard deviation of CH and TTH as discussed previously (~2 m and ~4 m, respectively).For the others, the discrepancy is larger (>6 m) as for the TM site.

Understanding Results with Apparent Foliage Profiles
The canopy height (CH) does not provide enough resolved information on the vertical structure of the forest systems, which can be very complex because of the presence of multiple layers of saplings, as well as undergrowth (i.e., tree ferns, ferns, bushes).Hence, it is preferable to consider the vertical profile of the apparent foliage (F app ).Such a profile is corrected from the extinction of the upper canopy.Parker et al. [48] have shown that F app is an important constraint for energy, water and nutrient flows through forest cover.This is due to the contrasted contributions of the different canopy levels to both photosynthesis and carbon storage [48].
Remote Sens.2016, 8, 43 14 nutrient flows through forest cover.This is due to the contrasted contributions of the different canopy levels to both photosynthesis and carbon storage [48].
As an example, we focus here on two different flight segments obtained on the Bélouve (BF) and Tamarind (TM) sites.The first site is very dense with continuous vegetation from ground to the canopy (as seen in the field), whereas the second one is composed of several distinct internal structures.Figure 8 shows the evolution of the apparent foliage as a function of distance along the transect in the BF.Two typical vertical profiles are also given.Certain profiles can show pronounced peaks, which identify the precise position of the tree crown or, on the contrary, smoother shapes due to the likely contribution of branches of nearby trees, lianas and, in the lower part of the profile, important undergrowth.Overall, because of the high density of trees in the tropical montane forest of BF, the LiDAR profiles mainly highlight only one vertical structure with one peak.This is also the case for the Cryptomeria (CM) site, but not for the same reason, because it is an exploited plot and there are generally no overlapping trees.On the contrary, for the TM site, which lies on the slope of Piton Maïdo, the convoluted tree trunks and complex development in response to storm winds lead to the existence of two superimposed layers, as can be seen in Figure 9.The profiles show there is generally an area with a lower density of vegetation between the two layers (between 4 and 8 m•agl), leading to less backscattered signal.This complex structure may be the source of discrepancies between airborne measurements and the in situ census made from the ground level.Another interesting and concrete conclusion is that the energy and water vapor fluxes between the forest and the atmosphere are mainly at the crown level of the trees for Tamarinds, even if undergrowth also contributes below ~8 m•agl.In contrast, for the BF site, these fluxes are distributed over the whole vertical forest structure.As an example, we focus here on two different flight segments obtained on the Bélouve (BF) and Tamarind (TM) sites.The first site is very dense with continuous vegetation from ground to the canopy (as seen in the field), whereas the second one is composed of several distinct internal structures.Figure 8 shows the evolution of the apparent foliage as a function of distance along the transect in the BF.Two typical vertical profiles are also given.Certain profiles can show pronounced peaks, which identify the precise position of the tree crown or, on the contrary, smoother shapes due to the likely contribution of branches of nearby trees, lianas and, in the lower part of the profile, important undergrowth.Overall, because of the high density of trees in the tropical montane forest of BF, the LiDAR profiles mainly highlight only one vertical structure with one peak.This is also the case for the Cryptomeria (CM) site, but not for the same reason, because it is an exploited plot and there are generally no overlapping trees.On the contrary, for the TM site, which lies on the slope of Piton Maïdo, the convoluted tree trunks and complex development in response to storm winds lead to the existence of two superimposed layers, as can be seen in Figure 9.The profiles show there is generally an area with a lower density of vegetation between the two layers (between 4 and 8 m¨agl), leading to less backscattered signal.This complex structure may be the source of discrepancies between airborne measurements and the in situ census made from the ground level.Another interesting and concrete conclusion is that the energy and water vapor fluxes between the forest and the atmosphere are mainly at the crown level of the trees for Tamarinds, even if undergrowth also contributes below ~8 m¨agl.In contrast, for the BF site, these are distributed over the whole vertical forest structure.

Leaf Area Index
LAI is also a key parameter linked to the plant respiration and photosynthesis, as explained by Gower and Norman [49].It is important for vegetation growth (carbon sequestration) estimation [50].Such a parameter is also a strong constraint for forest ecosystem modelling.It characterizes the forest interaction surface and exchange efficiency with the atmosphere.

LiDAR-Derived LAI
The LAI has been retrieved from each individual LiDAR profile to complement the characterization of the sampled forest sites.The LAI values are found to be log-normal distributed for all sampled forest sites, with the mean LAI (LAI) ranging from about 3.5 to 7.8 m 2 /m 2 (Table 3).The standard deviation is between 2.5 and 3.9 m 2 /m 2 , not necessarily correlated with the one of CH (Table 3).The tropical rain forests of Réunion Island (ML-150) have been shown to be associated with the higher mean LAI of 7.8 m 2 /m 2 .Such a value contrasts with the one retrieved for the tropical montane cloud forest of Bélouve, which is shown to be ~5 m 2 /m 2 on average.This difference may be explained in terms of plant nutrient supply between these two forest categories [51].As an example, the histogram of LAI derived from the LiDAR for the Bélouve site is also presented in Figure 10.The LAI value varies much for the same site from one point to another, likely due to very important difference in terms of nutrient availability in the ground.
15 explained in terms of plant supply between these two forest categories [51].As an example, the histogram of LAI derived from the LiDAR for the Bélouve site is also presented in Figure 10.The LAI value varies much for the same site from one point to another, likely due to very important difference in terms of nutrient availability in the ground.

Inter-Comparison and Discussion
In this part, we will compare the LiDAR-and MODIS-derived LAI and discuss previous results published in the scientific literature.
The comparison with our observations was possible only on the large sampled areas due to the low spatial resolution of MODIS.Even though the LAI algorithm from satellite observations usually assumes that the ground is flat and does not handle mutual shadows due to the terrain, the agreement is very good on Bélouve (ground almost flat; Table 1) and sites situated on the foothills of the Piton Maïdo (TM and CM, ~5 m 2 /m 2 ; ground with a certain slope; Table 1).LiDAR-derived mean LAIs for ML-150, -250 and -550 (7.8, 7.5 and 6.7 m 2 /m 2 , respectively) present higher differences with MODIS-derived LAI (5.9 m 2 /m 2 ).The small slopes in these sites (Table 1) are not related to a better agreement, as the sampled surface is too small.We believe that LiDAR measurements offer a better assessment of LAI than MODIS, as this LiDAR metric has already been validated by Tang et al. [21] Figure 10.Histogram of the LAI derived from the airborne LiDAR for the tropical mountain rain forests of Bélouve (BF).

Inter-Comparison and Discussion
In this part, we will compare the LiDAR-and MODIS-derived LAI and discuss previous results published in the scientific literature.
The comparison with our observations was possible only on the large sampled areas due to the low spatial resolution of MODIS.Even though the LAI algorithm from satellite observations usually assumes that the ground is flat and does not handle mutual shadows due to the terrain, the agreement is very good on Bélouve (ground almost flat; Table 1) and sites situated on the foothills of the Piton Maïdo (TM and CM, ~5 m 2 /m 2 ; ground with a certain slope; Table 1).LiDAR-derived mean LAIs for ML-150, -250 and -550 (7.8, 7.5 and 6.7 m 2 /m 2 , respectively) present higher differences with MODIS-derived LAI (5.9 m 2 /m 2 ).The small slopes in these sites (Table 1) are not related to a better agreement, as the sampled surface is too small.We believe that LiDAR measurements offer a better assessment of LAI than MODIS, as this LiDAR metric has already been validated by Tang et al. [21] and MODIS-derived LAI saturates at certain value levels because of the nature of the tool [31].In addition, the clumping of vegetation structure is less considered by the MODIS LAI algorithm; so it is normal that the LAI is underestimated.As previously explained, it is difficult to conclude for both ML-250 and ML-550, because the number of samples obtained over these sites is not significant.
Mature evergreen tropical forests usually have large LAI, more than 4 m 2 /m 2 , as shown for example in Doughty and Goulden [52] using monthly MODIS observations over tropical forests of Brazil.Cristiano et al. [53] also found mean LAI to be larger than 7 m 2 /m 2 for native subtropical forests of Argentina, Brazil and Paraguay.Moreover, the mean LAI derived from our LiDAR measurements is within the range of values deduced from the LiDAR-or tower-derived cumulative LAI of Tang et al. [21], which give LAI between 5 and 9 m 2 /m 2 for secondary and old-growth forests.Asner et al. [54] performed a global synthesis of LAI from various ecological and remote sensing studies.The mountain tropical forest of Bélouve, the TM and the CM sites are associated with mean LAI very close to the one compiled by these authors, which found for tropical deciduous and evergreen broadleaves LAI between 3.9 and 4.8 m 2 /m 2 on 78 samplings, with a maximal value close to 9 m 2 /m 2 .Our standard deviations are in the same range as obtained by Asner et al. [54], which give values from 0.7 to 4.3 m 2 /m 2 .new technical insights on the of LiDAR to penetrate through dense forest, on the choice of the laser footprint and spatial sampling, taking into account the forest heterogeneity, and on the retrieval of the LAI, in complementarity with the method of Tang at al. [21].
Finally, this campaign was an opportunity to compose an original and diverse LiDAR database that will help further works on remote sensing of tropical forests, mainly for the inter-annual evolution of the forest cover of Réunion Island, which is a rich United-European (French) diverse tropical ecosystem listed as World Heritage by UNESCO.

Figure 1 .
Figure 1.Location of the study sites and topography of Réunion Island.

Figure 1 .
Figure 1.Location of the study sites and topography of Réunion Island.

Figure 2 .
Figure 2. Autonomous payload (~80 kg) implemented on an ultra-light aircraft, including (a) the Ultraviolet LiDAR for Canopy Experiment (ULICE).The other instruments are also onboard: (b) a Tetracam ADC (Agricultural Digital Camera) air camera is used to get the photosynthesis activity index (Normalized Difference Vegetation Index (NDVI)) images over the forest canopy; the ancillary positioning instrument, called the MTi-G system, consists of a Global Positioning System (5-m accuracy) and an inclinometer (0.7° accuracy); (c) a Vaisala PTU-300 pressure/temperature/relative humidity probe is used for altitude correction and control of the tropical high humidity conditions that could affect the transmittance of LiDAR optics.

Figure 2 .
Figure 2. Autonomous payload (~80 kg) implemented on an ultra-light aircraft, including (a) the Ultraviolet LiDAR for Canopy Experiment (ULICE).The other instruments are also onboard: (b) a Tetracam ADC (Agricultural Digital Camera) air camera is used to get the photosynthesis activity index (Normalized Difference Vegetation Index (NDVI)) images over the forest canopy; the ancillary positioning instrument, called the MTi-G system, consists of a Global Positioning System (5-m accuracy) and an inclinometer (0.7 ˝accuracy); (c) a Vaisala PTU-300 pressure/temperature/relative humidity probe is used for altitude correction and control of the tropical high humidity conditions that could affect the transmittance of LiDAR optics.

Figure 3 .
Figure 3. Example of ground echo and canopy top detection from a range-corrected LiDAR signal explained in arbitrary units (a.u.) for a section of flight over the tropical forest of Bélouve (BF).The difference between the ranges of these two points yields the canopy height (CH).Note that the y-axis is not the ground elevation, but the distance from the emitter.ULA, ultra-light aircraft.

Figure 3 .
Figure 3. Example of ground echo and canopy top detection from a range-corrected LiDAR signal explained in arbitrary units (a.u.) for a section of flight over the tropical forest of Bélouve (BF).The difference between the ranges of these two points yields the canopy height (CH).Note that the y-axis is not the ground elevation, but the distance from the emitter.ULA, ultra-light aircraft.

Figure 5 .
Figure 5. Distribution of the canopy height (CH) computed using a varying fraction of the LiDAR footprints in order to simulate an effective sampling distance between two consecutive footprints along the ULA ground-track of 1 m (a); 10 m (b) and 100 m (c).The two-dimensional representation is given in (d).

Figure 5 .
Figure 5. Distribution of the canopy height (CH) computed using a varying fraction of the LiDAR footprints in order to simulate an effective sampling distance between two consecutive footprints along the ULA ground-track of 1 m (a); 10 m (b) and 100 m (c).The two-dimensional representation is given in (d).

Figure 6
Figure 6 gives examples of samplings performed over the TM and CM sites with sampling distances ∆X « 1 m and ∆Y « 10 m.The presence of both valleys and clear areas observed on the ADC-air vegetation camera image can explain the inhomogeneity in the LiDAR CH measurements.

Figure 6 .
Figure 6.Horizontal samplings performed on the (a) Tamarind (TM) and (b) Cryptomeria (CM) sites with a laser footprint of 10 m.The horizontal sampling distances along and perpendicular to the ground track are ~1 and ~10 m, respectively.The brown color corresponds to the ground numerical model at a 5-m resolution provided by the Parc National de la Réunion (J.-C.Notter, personal communication).The locations of the sub-plots where in situ measurements were performed are highlighted using thick black lines.

Figure 7 .
Figure 7. Distributions of the tree top height (TTH) from in situ (red) and canopy height (CH) from LiDAR (blue) measurements for the sub-plots of Tamarind site (TM), the Cryptomeria site (CM) and Mare-Longue tropical rain forest sites at 150 m (ML-150) and 250 m (ML-250) amsl.The overlapping parts are in purple.Uncertainty on retrieved heights is 4 m (2 m) for in situ (LiDAR) measurements.

Figure 7 .
Figure 7. Distributions of the tree top height (TTH) from in situ (red) and canopy height (CH) from LiDAR (blue) measurements for the sub-plots of Tamarind site (TM), the Cryptomeria site (CM) and Mare-Longue tropical rain forest sites at 150 m (ML-150) and 250 m (ML-250) amsl.The overlapping parts are in purple.Uncertainty on retrieved heights is 4 m (2 m) for in situ (LiDAR) measurements.

Figure 8 .
Figure 8. Evolution of the forest vertical profile of the LiDAR-derived apparent foliage (Fapp) along a transect of the Bélouve site (BF, left panel).Two specific vertical profiles are shown in the right panel, and their locations are highlighted in transparent gray in the left panel.

Figure 9 .
Figure 9. Evolution of the forest vertical profile of the LiDAR-derived apparent foliage (Fapp) along a transect of the Tamarind site (TM, left panel).Two specific vertical profiles are given in the right panel, and their locations are highlighted in transparent gray in the left panel.The gap between the distances

Figure 8 .
Figure 8. Evolution of the forest vertical profile of the LiDAR-derived apparent foliage (F app ) along a transect of the Bélouve site (BF, left panel).Two specific vertical profiles are shown in the right panel, and their locations are highlighted in transparent gray in the left panel.

Figure 8 .
Figure 8. Evolution of the forest vertical profile of the LiDAR-derived apparent foliage (Fapp) along a transect of the Bélouve site (BF, left panel).Two specific vertical profiles are shown in the right panel, and their locations are highlighted in transparent gray in the left panel.

Figure 9 .
Figure 9. Evolution of the forest vertical profile of the LiDAR-derived apparent foliage (Fapp) along a transect of the Tamarind site (TM, left panel).Two specific vertical profiles are given in the right panel, and their locations are highlighted in transparent gray in the left panel.The gap between the distances

Figure 9 .
Figure 9. Evolution of the forest vertical profile of the LiDAR-derived apparent foliage (F app ) along a transect of the Tamarind site (TM, left panel).Two specific vertical profiles are given in the right panel, and their locations are highlighted in transparent gray in the left panel.The gap between the distances of 70 to 90 m corresponds to LiDAR shots with big pointing angles (the angle between the actual LiDAR line of sight and the nadir direction is larger than 20 ˝), which were not considered in our study.

Figure 10 .
Figure 10.Histogram of the LAI derived from the airborne LiDAR for the tropical mountain rain forests of Bélouve (BF).

Table 1 .
Main characteristics of the study sites and associated available LiDAR profiles.

Table 2 .
Summary of the ULICE characteristics.

Table 2 .
Summary of the ULICE characteristics.