Sun–induced fluorescence heterogeneity as a measure of functional diversity

Plant functional diversity, defined as the range of plant chemical, physiological and structural properties within plants, is a key component of biodiversity which controls the ecosystem functioning and stability. Monitoring its variations across space and over time is critical in ecological studies. So far, several reflectance-based metrics have been tested to achieve this objective, yielding different degrees of success. Our work aimed at exploring the potential of a novel metric based on far-red sun-induced chlorophyll fluorescence (F760) to map the functional diversity of terrestrial ecosystems. This was achieved exploiting high-resolution images collected over a mixed forest ecosystem with the HyPlant sensor, deployed as an airborne demonstrator of the forthcoming ESA-FLEX satellite. A reference functional diversity map was obtained applying the Rao's Q entropy metric on principal components calculated on key plant functional trait maps retrieved from the hyperspectral reflectance cube. Based on the spectral variation hypothesis, which states that the biodiversity signal is encoded in the spectral heterogeneity, two moving window-based approaches were tested to estimate the functional diversity from continuous spectral data: i) the Rao's Q entropy metric calculated on the normalized difference vegetation index (NDVI) and ii) the coefficient of variation (CV) calculated on hyperspectral reflectance. Finally, a third moving window approach was used to estimate the functional diversity based on F760 heterogeneity quantified through the calculation of the Rao's Q entropy metric. Results showed a strong underestimation of the functional diversity using the Rao's Q index based on NDVI and the CV of reflectance. In both cases, a weak correlation was found against the reference functional diversity map (r2 = 0.05, p < .001 and r2 = 0.04, p < .001, respectively). Conversely, the Rao's Q index calculated on F760 revealed similar patterns as the ones observed in the reference map and a better correlation (r2 = 0.5, p < .001). This corroborates the potential of far-red F for assessing the functional diversity of terrestrial ecosystems, opening unprecedented perspectives for biodiversity monitoring across different spatial and temporal


Introduction
Understanding the ecological consequences of biodiversity loss recently became a central issue in ecological research (Oliver et al., 2015;Cardinale et al., 2012). The reason for that lies in the well-established relationship existing between biodiversity and provision of ecosystem services (e.g., Isbell et al., 2015;Musavi et al., 2017;Schwalm et al., 2017;Midgley, 2012;Wang and Gamon, 2019) and in particular in the critical role of biodiversity in controlling the ecosystem productivity and its temporal stability (Isbell et al., 2009(Isbell et al., , 2015Tilman et al., 1997).
Climate change and anthropogenic activities have already affected biodiversity severely, leading to a significant loss of ecosystem functions (Hooper et al., 2012;Oliver, 2016). For these reasons, the international community is engaging urgent actions to avert further biodiversity loss and degradation of ecosystem services, as framed e.g. in the European Biodiversity Strategy for 2020 (European Commission, 2011).
Biodiversity is a multidimensional concept that includes taxonomic diversity (the number and relative abundance of species), functional diversity (reflecting species' ecological roles), phylogenetic diversity (evolutionary differences of species), and genetic diversity (referring to genetic variation of species). In particular, the plant functional diversity, defined as the range of plant chemical, physiological and structural properties within plants, takes into account the variability of plant traits (PTs), which varies both inter-and intra-species. The unevenness of PTs constitutes a determinant of the ecosystem functioning (Diaz et al., 2007;Musavi et al., 2015;Ruiz-Benito et al., 2014;Tilman et al., 1997), thus quantifying the degree of functional variation in plant communities constitutes a key challenge in biodiversity research .
Although measuring biodiversity consistently across space and over time is critical to engage conservation actions, this objective is still far to be reached (Pettorelli et al., 2017). As a matter of fact, traditional field sampling of biodiversity is generally considerably demanding in terms of time, costs and human resources Gamon, 2019, Rocchini et al., 2019). Due to its inherent capacity of providing spatially and temporally continuous information about the Earth's surface, remote sensing has a great potential for detecting, quantifying, assessing and forecasting biodiversity at different spatio-temporal scales in a reliable, consistent and repeatable way Skidmore and Pettorelli, 2015;Kuenzer et al., 2014;Pettorelli et al., 2016).
So far, the studies exploring biodiversity using remote sensing have exploited three main approaches: i) mapping the distribution of habitat or species, ii) estimating biodiversity through spectral diversity and iii) estimating the functional diversity through plant functional traits (Wang and Gamon, 2019).
The first approach consists in the exploitation of remote sensing data to produce maps of habitat, functional types or species depending on the spectral and spatial resolution of the data (Fassnacht et al., 2016;Ghosh et al., 2014). These maps can also be used to assess the taxonomic diversity using traditional biodiversity metrics based on the species presence and abundance (e.g., species richness, Shannon's index).
The second approach directly exploits the spectral dimensionality of remote sensing data to infer the functional diversity from variations in the spectral patterns across space (Asner et al., 2017;Gholizadeh et al., 2019, Gholizadeh et al., 2018aNagendra, 2001). This approach is based on the spectral variation hypothesis (SVH) (Palmer et al., 2000(Palmer et al., , 2002, which states that the biodiversity signal is encoded in the spectral variability, therefore the larger the spectral heterogeneity is, the higher is biodiversity. An advantage of the SVH approach is that the spectral variability can be readily computed from remotely sensed data without the need to map the species or to retrieve plant functional traits (Rocchini et al., 2013). The spectral heterogeneity (and biodiversity in turn) can be quantified through a number of statistical metrics applied to full-spectra, vegetation indices or transformed spectra synthesizing the spectral information in a principal component space (Wang and Gamon, 2019). These metrics can be simple dispersion measures, e.g. the coefficient of variation (CV) of spectral reflectance , or more advanced measures based on the information theory, e.g. the Rao's Q index (Rao, 1982;Botta-Dukát, 2005). Among the diversity metrics, the latter index was widely used as a measure of functional diversity (e.g., Botta-Dukát, 2005;Scherer-Lorenzen et al., 2007;Weigelt et al., 2008) due to its capacity to include both the relative abundance of species and the functional differences within species. Its strong potential as remote sensing-based measure to quantify the functional diversity was demonstrated by recent studies (Rocchini et al., 2019. Compared to simple metrics such as the Shannon's index-which only accounts for the relative abundance of each pixel value-the Rao's Q index considers the sum of all the pixel pairwise distances in a spectral space, each of which is multiplied by the relative abundance of each pair of pixels within the moving window (Rao, 1982).
The third approach is based on the quantification of the functional diversity through plant traits, which reflect phylogenetic components, resource limitations, environmental conditions and photoprotection strategies (Schneider et al., 2017;Schweiger et al., 2017). Plant traits influence the absorption, scattering and transmission of the incoming solar radiation, inducing variations in the optical properties of plants that can be detected by optical remote sensing . Based on this link, plant traits such as foliar pigment content, leaf water content and leaf area index can be retrieved from remotely sensed data using physically based approaches, statistical models trained with ground-measured plant traits or hybrid approaches (refer to Verrelst et al., 2018 for a review). Recent developments in imaging spectroscopy systems have greatly enriched the dimensionality of remote sensing data (Asner et al., 2012;Thompson et al., 2017) and expanded the range of detectable plant biochemical, physiological and structural properties that can contribute to assess diversity (Hill et al., 2019;Ustin and Gamon, 2010;Asner et al., 2012). Among them, recent studies demonstrated the feasibility of retrieving sun-induced chlorophyll fluorescence (F) from ultra-fine spectral resolution airborne imagery (e.g., Rascher et al., 2015;Rossini et al., 2015;Middleton et al., 2017;Tagliabue et al., 2019). F is a dynamic plant functional trait  that synthesizes the variability of multiple plant traits often considered relatively stable and integrates information about the plant physiological status. F is in fact a physical signal originating from the core of the photosynthetic machinery. Together with heat dissipation, F is a pathway that plants use to dissipate the excess of absorbed photosynthetically active radiation (APAR) not used for photochemistry (Mohammed et al., 2019). Since the photosynthetic efficiency affects the F emission, this link constitutes the rationale behind the use of F to track variations in the plant functional state.
In this framework, the aim of this study is to test the potential of a novel approach based on the use of F heterogeneity-defined as the spatial variability of F across neighbouring pixels-as a measure of the functional diversity of a mixed forest ecosystem. We hypothesize that F might be a promising signal for mapping the functional diversity due to its capacity to integrate information about plant trait variability and plant physiology. The F signal is in fact related to the variability of key leaf and canopy traits, such as nitrogen and chlorophyll concentration, which strongly influence photosynthesis Porcar-Castell et al., 2014). Moreover, F can track heterogeneous physiological behaviours related to divergent downregulation strategies (Ac et al., 2015;Paul-Limoges et al., 2018). While previous studies conducted at the leaf level using active fluorescence measurements laid the foundations for these hypotheses (Gamon et al., 2005;Gamon, 2015;Rascher et al., 2004;van der Tol et al., 2014), we aim at testing the potential of passive sun-induced fluorescence in the detection of functional diversity patterns. As far as we know, this is the first attempt trying to exploit sun-induced fluorescence in the context of the functional diversity estimation. We aim at answering the following research questions: • does far-red sun-induced chlorophyll fluorescence provide valuable information to improve the diversity mapping compared to state-ofthe-art measures based on reflectance?
• does far-red sun-induced chlorophyll fluorescence provide valuable information related to functional rather than taxonomic/species diversity?
To achieve these objectives, high-resolution airborne images acquired with the HyPlant imaging sensor  over a temperate mixed forest ecosystem were used to map the forest diversity using different approaches. Firstly, the species richness and the Shannon's index calculated on a species classification map were used to evaluate the taxonomic diversity across the study area. Secondly, three different approaches were evaluated to estimate the functional diversity: i) the calculation of the Rao's Q index on sun-induced chlorophyll fluorescence measured in the far-red (F 760 ); ii) the calculation of the Rao's Q index on the normalized difference vegetation index (NDVI) (Tucker, 1979), which is widely used in ecological applications (e.g., Gamon et al., 1995;Pettorelli et al., 2005) and iii) the calculation of the CV of hyperspectral reflectance, which is a benchmark approach for the diversity estimation using remotely sensed data Somers et al., 2015;Sakowska et al., 2019, Gholizadeh et al., 2019. The heterogeneity maps obtained were compared against a reference functional diversity map to evaluate the performances of the different metrics as functional diversity measures.

Study site and field data acquisition
The site selected for this study is the Hardt Forest, a temperate forest located in France (47°48′29" N, 7°26'53" E; Mulhouse; Alsace). The analysis was limited to a subset of the forest (i.e.,~90 ha) covered by the airborne overpasses.
The Hardt Forest is dominated by the presence of broadleaved species (~90%), with a sparse presence of deciduous and evergreen coniferous ones (~10%) (Tagliabue et al., 2016). The most common species in the main canopy layer are: European hornbeam (Carpinus betulus L.), pedunculate and sessile oak (Quercus robur L., Quercus petraea (Matt.) Liebl.), field maple (Acer campestre L.), small-leaved linden (Tilia cordata Mill.), Scots pine (Pinus sylvestris L.) and European larch (Larix decidua Mill.). The forest is structured in stands of at least 500 m size, which are characterised by a relative variability in terms of management stage. The regeneration areas, which are characterised by the presence of young trees (mainly hornbeams and oaks) planted with high density, are clearly distinguishable from the mature areas of the forest, where the species variability is higher, and the stem density is lower. The average diameter of the tree crowns is 10-15 m in the mature forest and < 10 m in the regeneration stands.
The region where the forest is located is temperate, with an average temperature of 22°C in summer and of 4°C in winter. The mean annual rainfall is 680 mm, with a maximum typically occurring between May and August.
A field campaign was conducted in the study site in the summer of 2013 to collect ground-based measurements to validate the airborne products. Plant traits were sampled within sampling sites measuring 20 × 20 m 2 . The position of the sampling sites, hereafter referred to as elementary sampling units (ESUs), was selected by forest experts along the forest tracks and recorded using a high precision Geo-XT GPS (Trimble, Sunnyvale, USA). The leaf chlorophyll content (LCC), leaf mass per area (LMA) and leaf water content (LWC) were estimated from destructive measurements on leaves collected in 12 ESU. In each ESU, ten leaves were collected for each dominant species in the main canopy layer (n ≃ 250). This sampling scheme only applies to broadleaf species, since the coniferous species were only sparse in this area of the forest. The samples were collected from sunlit leaves that were shotgun-sampled from the top branches of the canopy and were immediately placed in plastic bags that were then stored at −80°C until the laboratory analysis. The LCC (i.e., sum of chlorophyll a and b contents) was quantified measuring the absorbance at 645, 662 and 710 nm using a UVIKON XL spectrophotometer (BioTek Instruments, Winooski, USA) after extracting the pigments with hydroxide carbonate magnesium buffered with acetone. The extinction coefficients to quantify chlorophyll a and b contents were taken from Lichtenthaler and Buschmann (2001). The LMA and LWC were calculated after determining the fresh weight (Wf) and dry weight (Wd) (i.e., after oven-drying for 24 h at 80°C) of leaves as LMA = Wd/Area and LWC = (Wf-Wd)/Area, respectively. The area was measured using a LI-3000 planimeter (LI-COR Biosciences, Lincoln, USA). The LAI was estimated from digital hemispherical photos using the CAN-EYE software (https://www6.paca.inra. fr/can-eye/). The images -seven looking upward per ESU -were acquired in correspondence of 14 ESUs using a Sigma Camera (Sigma Corporation, Ronkonkoma, USA) equipped with a fisheye lens.

Airborne data acquisition and pre-processing
Hyperspectral images were acquired on June 16, 2013 using the airborne imaging spectrometer HyPlant . The sensor was developed within the framework of the FLuorescence EXplorer (FLEX) mission preliminary activities and was specifically designed as airborne demonstrator of the FLEX satellite. HyPlant is a lineimaging push-broom scanner consisting of two modules: i) the Dual Channel Imager (DUAL) covering the visible (VIS), near-infrared (NIR) and shortwave infrared (SWIR) spectral regions (370-2500 nm) with a FWHM of 4.0 nm (VIS-NIR) -13.3 nm (SWIR) for the reflectance (ρ) and vegetation indices (VIs) calculation; and ii) the Fluorescence Imager (FLUO) covering the red and far-red spectral regions (670-780 nm) at ultra-fine spectral resolution (FWHM≃0.25 nm) for the fluorescence retrieval. The two modules operate synchronously, ensuring the spatial matching of the DUAL and FLUO hyperspectral cubes.
The airborne acquisition started at 11:55 Central European Summer Time (CEST) under clear sky conditions and was completed in approximately 157 s. The flight was set to 610 m above the ground level (1 × 1 m 2 ground sampling distance) with a heading of 195°.
HyPlant raw data were dark current subtracted, radiometrically calibrated and geo-referenced with the CaliGeo software (Specim Ltd., Oulu, Finland) using the data recorded by the position and attitude sensor during the acquisition. The DUAL data were atmospherically corrected using the ATCOR-4 software (ReSe Applications GmbH, Langeggweg, Switzerland) to obtain top-of-canopy reflectance (Siegmann et al., 2019). The FLUO data were processed using a dedicated processing-chain specifically developed to retrieve F. A shadow mask obtained with a threshold based on red and near-infrared reflectance was applied to both the DUAL and FLUO images to exclude the inter-crown shadows from the analysis.

Forest species mapping
The reflectance cube obtained from the HyPlant DUAL module was used to map the forest species distribution over the study site. The classification was performed through a supervised machine learning method, targeting seven different tree species (i.e., Carpinus betulus L., Quercus robur L., Quercus petraea (Matt.) Liebl., Acer campestre L., Tilia cordata Mill., Pinus sylvestris L., Larix decidua Mill.) grouped into six classes: "Hornbeam", "Oak", "Maple", "Linden", "Pine" and "Larch". An additional "Defoliated trees" class was defined in order to detect trees, mainly hornbeams and maples, affected by crown defoliation due to a caterpillar infestation at the time of the campaign. Mature and regeneration patches of the forest were classified separately. This choice was motivated by the high variability of the spectral signal among trees of the same species at different development stage, caused by differences in the leaf biochemical properties as well as in the canopy architecture. The classification scheme described hereafter was repeated likewise for both, masking out alternatively either the mature or the regeneration stands. The training set was prepared selecting pure spectral endmembers for each class by integrating the knowledge derived from the field surveys and the visual interpretation of high resolution orthophotos. These endmembers consisted of polygons targeting pure tree crowns over a wide patch of the study area in order to capture the intra-specific spatial variability. Depending on the occurrence of the species and the crown size, 15 to 30 polygons, each consisting of 6 to 10 pixels, were collected for each class, constituting a training set of~700 pixels. The training set was used to train a support vector machine (SVM) algorithm with a radial basis function kernel (Vapnik, 1995) at clustering the pixels into the defined classes. The SVM is a non-parametric classifier based on the statistical learning theory proposed by Vapnik and Chervonenkis (1971). The rationale behind the classifier is the definition of an optimal n-dimensional hyperplane that maximises the separation among the classes while minimising the misclassification errors (Vapnik, 1995). The SVM algorithms have been used increasingly in RS in the past few years because of their effectiveness and suitability in handling high-dimensional data (Melgani and Bruzzone, 2004;Mountrakis et al., 2011;Pal and Mather, 2005). In particular, they have been extensively applied for tree species classification in different biomes using hyperspectral data, showing outperforming accuracies compared to other classification methods (Dalponte et al., 2009(Dalponte et al., , 2013Feret and Asner, 2013;Melgani and Bruzzone, 2004;Pal and Mather, 2004). A relevant property of the SVMs, making them suitable for hyperspectral data classification, is their low sensitivity to the so-called Hughes phenomenon (Hughes, 1968). This effect consists in a decrease of the classification accuracy when the ratio between the number of input features and the training samples exceeds a certain threshold, due to the fact that the estimation of the classifier parameters (e.g., the estimation of the covariance matrices in the case of the Maximum Likelihood classifier) becomes complicated. This phenomenon is frequent in RS because the availability of training samples is usually limited and is critical when using hyperspectral data since it requires using feature reduction techniques at the cost of a loss of information and time. The SVMs proved to be unaffected by this issue (Camps-valls and Bruzzone, 2005;Melgani and Bruzzone, 2004;Pal and Mather, 2004;Waske and Benediktsson, 2010), hence, all the features were used as input for the classifier.
For the accuracy assessment, a testing set composed of~400 pixels was selected on the image using a randomly stratified sampling scheme. In order to ensure the representativeness of the testing sample for all the classes, the numerosity of the testing sample was defined for each class according to the abundance of the corresponding species in the study area. Each testing sample was labelled through visual interpretation and used as ground truth in the validation process. The standard accuracy metrics (i.e., overall accuracy (OA), producer's accuracy (PA) and user's accuracy (UA)) were calculated from the confusion matrix generated by crossing the classification result with the ground truths. The PA was calculated for each class as the ratio between the number of correctly classified pixels and the total number of testing samples. The UA was calculated for each class as the ratio between the number of correctly classified pixels and the total number of pixels assigned to that class. The OA was calculated as the ratio between the number of correctly classified pixels and the total number of testing samples.

Plant functional trait retrieval
Four key PTs were mapped exploiting the airborne hyperspectral DUAL data cube. Based on their ecological relevance, we selected leaf chlorophyll content (LCC), leaf area index (LAI), leaf water content (LWC) and leaf mass per area (LMA). The selected traits were retrieved from the reflectance data through a Lookup Table (LUT) based inversion of the coupled leaf and canopy PROSPECT-4-INFORM radiative transfer model (RTM) (Atzberger, 2000;Jacquemoud and Baret, 1990). The RTM was parameterised based on previous knowledge about the variability of the model input parameters in the forest and on the results of a global sensitivity analysis. The model was run in forward to generate a LUT of 30,000 simulated spectra. Hence, the traits of interest were retrieved by inverting the model with a LUT-based approach. Regularisation options were used to minimise the drawbacks of ill-posedness: LCC was retrieved using a logarithmic minimum contrast cost function (Leonenko et al., 2013) and averaging the first ten solutions of the inversion; LAI was retrieved with a divergence measure cost function formalised by Kullback and Leibler (1951) and averaging the ten best solutions; LWC and LMA were obtained using the Root Mean Square Error (RMSE) as cost function and averaging the ten best solutions. The airborne retrievals were validated against the field data collected close ( ± 1 day) to the airborne overpass. This comparison showed consistency between the ground-based and airborne retrievals for all traits (RMSE LCC = 5.66 μg cm −2 , RMSE LAI = 0.51 m 2 m −2 , RMSE LWC = 0.0019 g cm −2 , RMSE LMA = 0.0035 g cm −2 ). Further details about the validation can be found in Tagliabue et al. (2019).

Far-red fluorescence retrieval
The far-red fluorescence map used in this study was retrieved from HyPlant ultra-fine resolution data at 760 nm (F 760 ) based on the spectral fitting method (SFM) (Cogliati et al., 2015) adapted to HyPlant observations (Cogliati et al., 2018). The description of the retrieval approach can be found in Tagliabue et al. (2019). The map was validated against ground-based F 760 retrievals measured in correspondence of seven different vegetated targets close to the airborne overpasses. This comparison showed a good agreement between the ground-based and airborne retrievals (RMSE = 0.41 mW m −2 sr −1 nm −1 ) (Tagliabue et al., 2019).

Diversity measures from remotely sensed data
The biodiversity across the study area was measured using three different approaches: i) estimating the taxonomic diversity through the forest species map; ii) estimating the functional diversity through the spectral diversity and iii) estimating the functional diversity through plant functional trait and sun-induced fluorescence variability.
A moving window approach was used to derive the above-mentioned diversity measures from the airborne maps. For each diversity measure, different moving window sizes (i.e., 3 × 3, 5 × 5 and 9 × 9 pixels) were tested to evaluate the dependency of the results on the size of the moving window. Results are shown for the 9 × 9 pixels size which allows capturing both the intra-and inter-crown variability.
The taxonomic diversity was estimated using traditional metrics calculated on the forest species map: the species richness and the Shannon's index (H) (Shannon, 1948). The species richness was calculated for each pixel as the number of different tree species within the considered moving window. The Shannon's index was calculated for each pixel of the species classification map as: where n is the number of species and p is the proportion of the area occupied by the species i respect to the moving window size.
The functional diversity across the study area was estimated using three approaches based on spectral heterogeneity, functional traits and sun-induced fluorescence heterogeneity.
Firstly, two different moving window-based approaches were applied to HyPlant-derived products to estimate the functional diversity based on the spectral variability hypothesis: the calculation of the Rao's quadratic entropy (Q) (Rao, 1982;Botta-Dukát, 2005;Ricotta, 2005) on NDVI and the calculation of the CV on hyperspectral reflectance. The Rao's Q index is a popular multidimensional functional diversity index, which is traditionally computed on functional attributes measured in the field. The index is calculated as: where d ij is the dissimilarity between two values i and j based on the Euclidean distance, and p i and p j are the relative proportions of the values i and j, respectively, within the considered moving window. Based on its definition, the Rao's Q index can be described as the expected spectral difference in the numerical value between two pixels drawn randomly with replacement from the pixels within a moving window. Respect to the Shannon's index that only accounts for the relative abundance of each value, the Rao's Q index also considers the overall spectral difference among the pixels within a moving window. This allows us to know if two pixels drawn randomly from the moving window are different and to quantify how much different they are. Among vegetation indices, NDVI was used because it is a benchmark index in most studies aimed at the remote quantification of biodiversity G. Tagliabue, et al. Remote Sensing of Environment 247 (2020) 111934 using index-based methods (see Wang and Gamon, 2019 for a review). Nevertheless, other vegetation indices were tested, which yielded similar results to those obtained with NDVI. The CV of spectral reflectance was calculated as the average CV for all the wavelengths between 400 and 2500 nm, as: where sd(ρ λ ) and mean(ρ λ ) are the standard deviation and the mean of reflectance at the wavelength λ, respectively, and n is the number of spectral bands.
Secondly, the Rao's Q index was calculated on plant functional traits to derive a reference functional diversity map. Here, the Rao's Q index was applied to the functional trait maps retrieved from the DUAL hyperspectral reflectance through inversion of the coupled PROSPECT-4-INFORM radiative transfer model. A Principal Component Analysis (PCA) was used to combine LCC, LAI, LWC and LMA while removing the statistical redundancy in the multiple selected traits. Before calculating the PCA, all traits were re-scaled between 0 and 1 to homogenise the units. The Rao's Q diversity was calculated on the first three PCs, which together explained 96% of the total variance of the four selected PTs. The calculation was performed for the entire image within the R software (R Core Team, 2019) using the spectralrao function stored in the GitHub repository https://github.com/mattmar/spectralrao . The mode multidimension was set to use multiple matrices as input for the Rao's Q diversity calculation.
Finally, the Rao's Q index was calculated on sun-induced chlorophyll fluorescence measured in the far-red (F 760 ).
Both NDVI and F 760 were converted to 8 bits and re-scaled between 0 and 1 before calculating the Rao's Q index to harmonize the range and variability of the two variables. The calculation of the Rao's Q index was performed in R (R Core Team, 2019) using the spectralrao function as for the reference map based on plant functional traits.
The performances of the F 760 -based Rao's Q index proposed in this study and of two state-of-the-art approaches based on spectral diversity (i.e., CV of reflectance and NDVI-based Rao's Q index) in estimating the functional diversity were finally evaluated by fitting regression models against the reference PC-based Rao's Q index map. The statistical analysis was performed in R (R Core Team, 2019). To avoid spatial autocorrelation, the models were fitted on 1000 pixels extracted fully randomly from the images with a minimum distance of 15 m among them.

Mapping of species and species diversity
The classification process applied on HyPlant high-resolution imagery enabled us to produce a thematic map of the distribution of the main tree species in the Hardt forest (Fig. 2a). Table 1 shows the confusion matrix obtained crossing the ground truth with the classification results. The overall accuracy (OA) obtained was 75.1%, while the producer's (PA) and user's (UA) accuracies ranged between 59 and 94% and between 63 and 97%, respectively, depending on the considered class. Scots pine was mapped with the highest accuracy due to its spectral dissimilarity compared to the broadleaves, which were characterised by a higher misclassification. The defoliated tree class showed a high UA (89.2%), but the low PA (58.6%) highlighted that the defoliated trees were sometimes not identified, probably depending on the degree of defoliation.
The thematic product obtained highlighted an uneven abundance and distribution of the different forest species. Overall, the most common species were oak and hornbeam, but the two species were characterised by a different spatial distribution. While oaks were found to prevail over hornbeams in the mature forest, the regeneration areas were overwhelmingly dominated by hornbeam. Pine was mainly limited to confined areas of the forest (e.g., the mature forest patch in the left-central part of the image). Maple, linden and larch presented a scattered distribution throughout the entire forest. The defoliated trees were widespread across the mature forest, but were almost absent in the The species richness and the Shannon's index calculated on the species map are shown in Fig. 2b and Fig. 2c, respectively. Both the diversity metrics based on the spatial distribution of the tree species across the site showed a low heterogeneity and range of variation across space. The patterns are consistent in the two images and reflect a low biodiversity in the study site in terms of species composition.

Estimation of the reference functional diversity map based on retrieved plant functional traits
An RGB composite of the principal components (PCs) obtained performing a PCA on key forest traits (i.e., LCC, LAI, LWC and LMA) is shown in Fig. 3. Different patterns in the forest are clearly distinguishable in the map (Fig. 3a). These patterns are related to the diversity in terms of the biochemical and structural traits mapped and reflect the differences in the ecological functioning. Magenta-orange areas correspond to high values in the first PC (PC1), which is strongly related to LAI and moderately correlated with LCC and LWC. These patches correspond to regeneration areas of the forest, where the trees are younger, and the canopy density is higher. Green-bluish colours, indicating a high influence of the second (PC2) and partly of the third PC (PC3), dominate the mature part of the forest. The sub-regions marked in Fig. 3 with "A", "B" and "C" (i.e., white squares) are representative of three different conditions encountered across the area. "A" marks an area at the edge between the regeneration and mature forest, where the managed plantation is characterised by high values in PC1 and PC2. "B" is a 100% mature forest area where the canopy is complexly structured. "C" marks a mixed regeneration/mature forest area where the regeneration stand presents high values in PC1 and PC3 and the mature stand has a different species composition compared to "A". PC1, which explained 77% of the total variance, was strongly dominated by LAI with some influence of LCC and LWC, as shown by the loadings and explained variance bar plots displayed in Fig. 3b and e. LCC showed the highest contribution in PC2 (Fig. 3c and f), which explained 11% of the variance. LMA, LWC and LAI also drove some of the variance in PC2, but with a negative contribution respect to LCC. Species richness calculated from the species map using a moving window of 9 × 9 pixels. The species richness is expressed as number of species; c) Shannon's index calculated from the species map using a moving window of 9 × 9 pixels. High Shannon's index values correspond to high taxonomic diversity. The map projection is UTM zone 32 N with datum WGS84.

Table 1
Confusion matrix obtained from the validation scheme. Columns represent the true classes while rows represent the classification results. The % of correctly classified pixels is reported along the upper-left to lower-right diagonal. The omission errors (% of pixels that belong to a class but were classified into another) are reported down the columns, while the commission errors (% of pixels that were classified into a class but do not belong to that class) are reported across the rows. The user's (UA) and producer's (PA) accuracies are reported in the last column and row, respectively. The overall accuracy (OA) is also displayed in the lower-right corner. G. Tagliabue, et al. Remote Sensing of Environment 247 (2020) 111934 PC3, which explained 8% of the variance, was strongly driven by LMA, with a minor contribution of LWC, LAI and LCC ( Fig. 3d and g). Overall, PC1, PC2 and PC3 explained 96% of the variance in the selected traits. The functional diversity map obtained applying the Rao's Q diversity metric on the PCs calculated on plant functional traits is shown in Fig. 4a. The patterns in the map revealed a high diversity in the forest, highlighting the presence of patches characterised by different entropy. The highest Rao's Q values were observed in the regeneration stands of the forest, which were also characterised by a considerable spatial variability. In contrast, the mature forest appeared in general more homogeneous, but high diversity spots were also observed. This map constitutes the reference against which the functional diversity maps based on the spectral and sun-induced fluorescence heterogeneity will be discussed.
The use of moving windows of different sizes (from 3 × 3 pixels up to 9 × 9 pixels) to calculate the Rao's Q diversity did not significantly affect the observed entropy patterns (r 2 = 0.88, p < .001 and r 2 = 0.71, p < .001 comparing the Rao's Q index calculated using a moving window of 9 × 9 against moving windows of 5 × 5 and 3 × 3 pixels, respectively).

Estimation of functional diversity based on spectral and sun-induced fluorescence heterogeneity
The functional diversity map obtained as CV of spectral reflectance between 400 and 2500 nm is shown in Fig. 4b. The map showed a clear separation between mature and regeneration forest, with generally low diversity values in the young stands, which are located in the middle and close to the top of the image. The lowest diversity was observed in correspondence of the young patches of the forest, while the mature forest was characterised by a considerable variability.
The functional diversity maps obtained calculating the Rao's Q index on NDVI and F 760 are shown in Fig. 4c and d, respectively. The NDVI-based Rao's Q was characterised by generally low values across the entire image, with patterns similar to those observed in the CV of reflectance. In particular, very low entropy values were observed in the regeneration stands, with Rao's Q values around 0.008-0.01. The highest diversity was observed in correspondence of the mature areas of the forest characterised by a significant presence of coniferous species, e.g., the area on the left in respect to the regeneration stand in the centre of the image. In these areas, the Rao's Q diversity reached values up to~0.05. The F 760 -based Rao's Q showed a considerably larger heterogeneity across the image. The regeneration stands showed the highest entropy (~0.2), but high entropy values were also observed in the mature forest. Overall, the F 760 -based Rao's Q was substantially higher than the NDVI-based Rao's Q.

Comparison between the reference and the spectral and sun-induced fluorescence heterogeneity maps
The results of the regression analysis performed between the functional diversity maps based on the spectral variability (i.e., CV of reflectance, NDVI-based Rao's Q index), the one based on the Rao's Q index calculated on F 760 and the reference PC-based functional diversity map are shown in Fig. 5. A weak correlation was observed between the reference PC-based Rao's Q index and the CV of reflectance (r 2 = 0.04, Fig. 3. a) RGB composite of the first three Principal Components (PCs) derived from the Principal Component Analysis performed on leaf chlorophyll content (LCC), leaf area index (LAI), leaf water content (LWC) and leaf mass per area (LMA) retrieved from the airborne hyperspectral images. The colours are the result of the combination of the first (PC1), second (PC2) and third (PC3) principal components in the R, G and B channels, respectively. The sub-regions marked with "A", "B" and "C" (white squares) are representative of three different conditions encountered across the area. "A" marks an area at the edge between the regeneration and mature forest, where the managed plantation is characterised by high values in PC1 and PC2. "B" is a 100% mature forest area where the canopy is complexly structured. "C" marks a mixed regeneration/mature forest area where the regeneration stand presents high values in PC1 and PC3 and the mature stand has a different species composition compared to "A". On the right, bar plots of the loadings of b) PC1, c) PC2 and d) PC3 and of the explained variance of e) PC1, f) PC2 and g) PC3. The map projection is UTM zone 32 N with datum WGS84. Fig. 4. Functional diversity maps obtained using different metrics calculated with a 9 × 9 pixels moving window: a) Rao's Q entropy calculated on principal components (PCs); b) average coefficient of variation (CV) calculated for all wavelengths from 400 to 2500 nm; c) Rao's Q entropy calculated on NDVI and d) Rao's Q entropy calculated on F 760 . The map projection is UTM zone 32 N with datum WGS84. G. Tagliabue, et al. Remote Sensing of Environment 247 (2020) 111934 p < .001) and between the PC-based Rao's Q index and the NDVIbased Rao's Q index (r 2 = 0.05, p < .001). Conversely, the Rao's Q index calculated on F 760 showed a strong positive correlation with the reference PC-based map (r 2 = 0.5, p < .001). The relationship between the latter two variables was also positive in terms of absolute Rao's Q index values. A relative bias of 10.8% was observed between the F 760 -based and the PC-based Rao's Q index, while the Rao's Q index based on NDVI showed a bias of −83.5% respect to the reference PCbased Rao's Q index.

Discussion
Measuring biodiversity across space and over time is a pivotal objective in ecology. In this respect, remote sensing can be a powerful tool as it provides the means to overcome several drawbacks related to the traditional biodiversity estimation approaches based on field data (Wang and Gamon, 2019;Pettorelli et al., 2016;Nagendra, 2001). However, there is no agreement yet on the remotely sensed data and metrics to be used to infer biodiversity, as well as the pros and cons of different products (i.e., reflectance, VIs and F) to map biodiversity (Skidmore and Pettorelli, 2015;Ma et al., 2019;Turner, 2014).
In this study, we tested different approaches to estimate the taxonomic and functional diversity of a mid-latitude forest ecosystem. The taxonomic diversity across the study area was inferred calculating the species richness and the Shannon's index on a species map obtained using a semi-automatic classification approach. This traditional approach for assessing biodiversity provides information about the species distribution based on their taxonomic identity and can be relevant for e.g. monitoring the biodiversity patterns and detecting species changes or losses over time. However, it cannot provide any information about the intra-specific variability in functional traits, which constitutes a determinant of the functional diversity and presumably it is more related to the ecosystem stability rather than to the taxonomic diversity alone (Díaz and Cabido, 2001). As a matter of fact, the functional traits might vary within the same species as much as between different species; vice versa, different species might be characterised by similar functional traits thus not contributing to the functional diversity (Ma et al., 2019;Schneider et al., 2017). For this reason, disregarding the intra-specific variability might strongly risk either masking or overestimating the functional diversity. Schneider et al. (2017) further suggested that in relatively species-poor temperate forests, such as the one analysed in this study, ignoring the functional diversity typically leads to a strong underestimation of biodiversity. This statement was confirmed by our findings, which showed that the use of diversity metrics based on the species composition led to diametrically opposed results compared to their application to remotely sensed, spatially continuous variables capable of grasping the intra-specific variability of PTs.
However, applying the entropy metrics to continuous RS data under the spectral variation hypothesis (Palmer et al., 2000) is not the only requisite to properly map the functional diversity. In order to capture all the intra-specific variability, two other critical factors need to be addressed: the choice of the entropy metric to be used to infer the functional diversity and the spectral variable to be used as input for the calculation. In this study, we used the Rao's Q index that was found to be a good performing indicator of the spatial heterogeneity compared to simpler metrics that do not consider the pairwise numerical difference between the values encountered within the moving window Rocchini et al., 2018Rocchini et al., , 2018bRocchini et al., 2019). Regarding the choice of the spectral variable to be used to infer the functional diversity, it must be considered that RS offers a variety of different tools. So far, different approaches have been proposed to synthesize the intra-specific variability of PTs, e.g., the use of VIs related to the vegetation structure and greenness Rocchini et al., 2018Rocchini et al., , 2018b and the use of spatialised retrievals of PTs (Schneider et al., 2017). Regardless of the methodology used to quantify biodiversity, both approaches are essentially based on reflectance, i.e., on the amount of radiation reflected by vegetation as a function of its biochemical and structural properties. In this study, we tested a completely novel approach for assessing the functional diversity based on the exploitation of the F emission. The strength of the F signal relies in its inherent capacity of providing insights into the ecosystem functioning (refer to Mohammed et al., 2019 for a review). Since F is linked to photochemistry being one of the two mechanisms of dissipating the excess of absorbed light, F carries information on the light use efficiency, which constitutes the unknown in the relationship linking the APAR to the gross primary productivity (Porcar-Castell et al., 2014).
Our results demonstrated the potential of using F to map the functional diversity over reflectance or traditional VIs used in ecology (e.g., Wang et al., 2016;Féret and Asner, 2014). Compared to the reflectancebased approaches tested, the use of the F 760 -based Rao's Q metric showed the highest correlation against the PC-based Rao's Q index used as reference (r 2 = 0.5, p < .001). Conversely, the average CV of reflectance and the NDVI-based Rao's Q index only explained a small part of the variance in the functional diversity patterns (r 2 = 0.04, p < .001 and r 2 = 0.05, p < .001, respectively).
In terms of patterns, both the reflectance-and NDVI-based functional maps clearly showed a strong underestimation of the entropy across the study area compared to the PC-and F-based maps (Fig. 4b  and c). This was particularly evident in the regeneration stands, where the extremely low entropy obtained using the approaches based on reflectance highlighted that reflectance cannot grasp most of the heterogeneity in the functional diversity. The fact that this effect is emphasized in the regeneration patches is probably related to the saturating effect of near-infrared reflectance in these areas, which is related to the high biomass and compact structure of the canopy in the managed stands. This suggests that F 760 might be a powerful tool for mapping the functional diversity especially in environments characterised by dense vegetation with high biomass and greenness (e.g., tropical and boreal forest ecosystems). This hypothesis is underpinned by previous studies conducted in such environments using active fluorescence measurements. Gamon et al. (2005) used optical indices, active fluorescence measurements and gas exchange measurements to investigate the photosynthetic performance of tropical forest canopies. While optical indices such as NDVI and simple ratio were closely related to the canopy structure and APAR of different species, they were not able to capture their heterogeneous physiological behaviour. Conversely, chlorophyll fluorescence, measured using active measurements, was able to detect divergent downregulation strategies among evergreen species. Similarly, Rascher et al. (2004) investigated the functional diversity of photosynthesis in tropical rainforest mesocosms under drought stress, observing that chlorophyll fluorescence was able to capture the heterogeneous response of individual rainforest trees expressing different downregulation strategies. While these studies showed the possibility to detect variations related to divergent downregulation mechanisms using active fluorescence, the potential of suninduced chlorophyll fluorescence in this respect needs to be further tested in dedicated studies performed in stress conditions.
In this study, the better performance of F 760 in predicting the functional diversity compared to NDVI is largely explained by the strong influence of the variability of PTs on the F signal. While NDVI is completely driven by the APAR, which is a function of the green biomass of vegetation (i.e., LAI), F 760 carries information on APAR as well as on other critical drivers of the ecosystem functioning. In particular, the F 760 signal is strongly influenced by LCC and LMA, as shown in the sensitivity analysis performed by Verrelst et al. (2015). LCC has been shown to be related to the maximum carboxylation rates (Vcmax) (e.g. Houborg et al., 2013). Vcmax has been shown to influence the F emission (Frankenberg and Berry, 2018, Martini et al., 2019. Moreover, LCC also scales positively with the leaf nitrogen (Niinemets et al., 1999;Croft et al., 2017), while LMA scales negatively with the leaf nitrogen, as demonstrated by the leaf economics spectrum (LES) (Wright et al., 2004;Osnas et al., 2013) and here observed in the opposed contribution of LMA and LCC on PC2 (Fig. 3c). PC2 and PC3 are strongly driven by LCC and LMA, respectively, meaning that a considerable portion of the variability in the reference functional diversity map is determined by the variation of LCC and LMA. Given the considerations above, our results show that F 760 is a promising indicator of functional diversity, given the capability of F of capturing the spatial variability of photosynthesis, which is ultimately influenced by functional traits such as LCC, and LMA.
It is worth noting that the spatial scale of this study might foster the spectral diversity-functional diversity relationships observed. In fact, the moving window of 9 × 9 pixels applied to the high spatial resolution HyPlant data (i.e., 1 × 1 m 2 ) allowed detecting both intraand inter-crown variations in the functional traits, enhancing the detection of biodiversity patterns. In our study, the functional diversity map estimated using a moving window of 9 × 9 pixels correlates quite well with the one estimated using a moving window of 3 × 3 pixels (r 2 = 0.71, p < .001). Such a small window size (3 × 3) allows tracking the intra-crown variability, so the good correlation between the above-mentioned maps suggests that the 9 × 9 maps are capturing a lot of intra-crown variability plus some inter-crown variability. The possibility to detect the intra-crown variability of functional traits may not be that strong when the spatial resolution is of the same order or larger than the tree crown size, since the intra-crown variability would be integrated in the pixel signal (Féret and de Boissieu, 2020). The biodiversity estimation based on the spectral heterogeneity is also typically hampered by the influence of confounding factors such as the presence of shadows or soil background (Gholizadeh et al., 2018;Wang et al., 2018). These elements increase the spectral heterogeneity, potentially leading to misleading conclusions. In this study, the availability of imaging spectroscopy data at high spatial resolution minimised this issue, since these confounding elements were excluded from the analysis. On the other hand, the exploitation of airborne data intrinsically limits the possibility to monitor biodiversity at large spatiotemporal scales. This recalls the need of integrating remote sensing data from multiple sources as well as in situ observations of functional traits and species distribution in order to detect and explain the biodiversity patterns at the global scale (Jetz et al., 2016). The rapid increase of the availability of Earth Observation satellites is already pushing forward in this direction, but future missions able to monitor the plant traits and functions at the global scale are still considered a priority of the scientific community (Jetz et al., 2016;Ma et al., 2019;Somers et al., 2015). In this respect, the ESA-FLEX satellite mission (Drusch et al., 2017), which will be specifically dedicated to observing F globally at unprecedented high spatial resolution (i.e., 300 × 300 m 2 ), might be a promising candidate for advancing the understanding and monitoring of the functional diversity across the Earth. However, previous studies indicate that the power of detecting biodiversity through spectral variability varies with spatial scales (Hakkenberg et al., 2018) and pixel size (Rocchini, 2007). Thus, we might expect that the ESA-FLEX pixels are too coarse for directly assessing biodiversity. The increasing availability of hyperspectral data at high spatial resolution provided by new spaceborne sensors (e.g., ASI-PRISMA, DLR-DESIS) might provide the means to overcome this issue. Hyperspectral observations might in fact be exploited in the future to either downscale the F signal using machine learning techniques (Duveiller and Cescatti, 2016;Wen et al., 2020) or unmix the contribution of shadows and soil background (Asner and Lobell, 2000;Asner and Heidebrecht, 2002).

Conclusions
In this study, multiple approaches for mapping the diversity of a mixed forest ecosystem were tested exploiting remotely sensed highresolution HyPlant imagery. Firstly, we used a classic approach for estimating the taxonomic diversity across the forest. The hyperspectral reflectance cube was used to classify the forest species with a semiautomatic approach, and two traditional biodiversity metrics (i.e., species richness and Shannon's index) were calculated based on the species map obtained. Secondly, we tested two approaches for assessing the functional diversity directly through the spectral diversity (the calculation of the Rao's Q entropy metric on NDVI and the calculation of the CV of hyperspectral reflectance). Finally, we explored for the first time the potential of using F 760 heterogeneity, quantified through the calculation of the Rao's Q entropy metric on F 760 , to estimate the functional diversity.
Our results strongly support the use of F 760 heterogeneity as a synthetic measure of the diversity of terrestrial ecosystems. The exploitation of the F 760 signal in fact outperformed reflectance and reflectance-based indices for estimating the functional diversity across the study area. Furthermore, evidence was shown that the F 760 heterogeneity carries information related to the functional diversity rather than the taxonomic diversity.
This opens unprecedented perspectives for assessing biodiversity from remote, paving the way of future studies aimed at exploring the F 760 diversity-biodiversity links at different spatial and temporal scales and across biomes. G. Tagliabue, et al. Remote Sensing of Environment 247 (2020) 111934

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.