Spectro-Temporal Heterogeneity Measures from Dense High Spatial Resolution Satellite Image Time Series: Application to Grassland Species Diversity Estimation

Grasslands represent a significant source of biodiversity that is important to monitor over large extents. The Spectral Variation Hypothesis (SVH) assumes that the Spectral Heterogeneity (SH) measured from remote sensing data can be used as a proxy for species diversity. Here, we argue the hypothesis that the grassland’s species differ in their phenology and, hence, that the temporal variations can be used in addition to the spectral variations. The purpose of this study is to attempt verifying the SVH in grasslands using the temporal information provided by dense Satellite Image Time Series (SITS) with a high spatial resolution. Our method to assess the spectro-temporal heterogeneity is based on a clustering of grasslands using a robust technique for high dimensional data. We propose new SH measures derived from this clustering and computed at the grassland level. We compare them to the Mean Distance to Centroid (MDC). The method is experimented on 192 grasslands from southwest France using an intra-annual multispectral SPOT5 SITS comprising 18 images and using single images from this SITS. The combination of two of the proposed SH measures—the within-class variability and the entropy—in a multivariate linear model explained the variance of the grasslands’ Shannon index more than the MDC. However, there were no significant differences between the predicted values issued from the best models using multitemporal and monotemporal imagery. We conclude that multitemporal data at a spatial resolution of 10 m do not contribute to estimating the species diversity. The temporal variations may be more related to the effect of management practices.


Introduction
Grasslands represent one of the largest land covers on Earth. They are an important source of biodiversity in farmed landscapes, thanks to their plant and animal composition [1,2]. This biodiversity supports many ecosystem services such as carbon regulation, erosion regulation, food production, biological control of pests and crop pollination [3,4]. However, global grassland surface area is decreasing, and grassland diversity is declining because of agriculture intensification, abandonment and urbanization [3], leading to a loss of biodiversity and associated services. To understand these effects, it is of utmost importance to determine and monitor grassland species diversity and composition over large extents.
Biodiversity can be characterized by alpha-diversity [5], which is related to the diversity in species of a community. Alpha-diversity is commonly measured by the species richness (number of species in the sampling area). However, this diversity can also be quantified with heterogeneity measures, such as the Shannon index [6], which combines richness and evenness (even abundance between species), and the Simpson index [7], which measures the dominance of species over the others.
Usually, ecologists measure and monitor biodiversity during field surveys. However, these surveys are time consuming, and they require important human and material resources, making them costly and limited in time and space [8,9]. Moreover, they tend to be influenced by the assessor [10], which can make the comparison between study areas difficult. Ecological field surveys are thus limited to a local scale, whereas there is an important need to monitor biodiversity over larger extents (national to international scales). To circumvent this issue, remote sensing appears to be an appropriate tool. Indeed, thanks to the broad spatial coverage and regular revisit frequency of satellite sensors, remote sensing provides continuous, regular and repeatable observations over large extents [11,12]. It has already proven its ability for habitat mapping [11,13], and it can be seen as an indirect approach for biodiversity estimation [8,9].
Considerable progress has been made in the remote sensing of biodiversity during the last few decades [10]. Many works are based on the Spectral Variation Hypothesis (SVH) [14,15], which assumes that the spectral heterogeneity in the image is correlated with the heterogeneity of the habitat. The diversity of species being related to the heterogeneity of the habitat [16,17], the spectral heterogeneity can be used as a proxy for species diversity [10].
The measures of grasslands' species diversity in the context of SVH have been discussed in the work of Oldeland et al. [18]. Most of the studies are based on species richness. However, this measure gives equal weight to every species, regardless of their proportion in the community. The contribution of rare individuals in the spectral heterogeneity can be doubtful. Abundance-based measures of species diversity, such as the Shannon index, give more weight to species with higher proportions. Therefore, these measures should be preferred to species richness in the context of SVH [18,19].
Many ways to quantify the Spectral Heterogeneity (SH) and to relate it to alpha-diversity have been developed in the remote sensing community [8]. SH has been quantified with the standard variation or the coefficient of variation of the Normalized Difference Vegetation Index (NDVI) [14,20], using Principal Components Analysis (PCA) [18,21] and with the mean Euclidean distance to the spectral centroid [15,18,19,22]. However, these measures do not describe well the variability in the spectral space. Recently, Féret and Asner [23] developed an original approach to account for both the spatial and the spectral information of imaging spectroscopy. Their approach is based on the hypothesis that species can be categorized according to their spectral reflectance. They performed an unsupervised clustering using the k-means algorithm, assigning each pixel of the image to a cluster called a "spectral species". Then, they computed the "spectral species distribution", which is the entropy (Shannon index) of the spectral species and which was found highly correlated with the ground Shannon index of a tropical forest. However, since they were using hyperspectral data with very high spatial resolution (2 m), simplifying steps were necessary prior to the clustering because k-means is not suitable for high dimensional data, resulting in a significant loss of spectral information and possibly of the heterogeneity.
Most of the aforementioned works were performed with hyperspectral data issued from a field spectroradiometer or an airborne sensor, thus with a very high spatial resolution. Although these works showed good results, they were limited to a very local scale, because of the costs involved by such a mission. Conversely, new satellite missions for continuous vegetation monitoring, such as Sentinel-2 [24], provide freely multispectral time series with high spatial and high temporal resolutions. Therefore, a tradeoff could be considered by using time series of satellite images to monitor grasslands biodiversity over large extents. Indeed, species communities differ in their temporal and seasonal behaviors, i.e., their phenology, making the phenological diversity related to the species diversity [25,26]. Therefore, in this study, we argue the hypothesis that the spectro-temporal heterogeneity of a community can be related to its species diversity, such as suggested by Rocchini et al. [10]: "Multispectral satellite sensors with high to very high spatial resolution and short revisit period, such as Sentinel-2, Venµs, and other high spatial resolution multispectral sensors may be good candidates for biodiversity mapping based on spectro-temporal variations". We could name this hypothesis the "Spectro-Temporal Variation Hypothesis" (STVH) in reference to the SVH.
However, the use of both the spectral and the temporal information in dense time series involves big data issues. Indeed, we have to deal with a high number of spectro-temporal variables, but with a small number of samples, because grasslands in Europe are relatively small objects in the landscape (around one hectare). Even with high spatial resolution sensors (around 10 m), on average, only a hundred pixels compose these grasslands, while there is about the same number of spectro-temporal variables during a year of acquisitions. Clustering algorithms and measures of spectral heterogeneity that are suitable for high dimensional data are thus required.
The objective of this study is to verify if the spectro-temporal variations in grasslands are related to their species diversity, using dense Satellite Image Time Series (SITS) with a high spatial resolution (10 m). To verify this hypothesis, we propose to link spectro-temporal and spectral-only heterogeneity measures derived from the unsupervised clustering of grasslands to their species diversity through linear regression models. To address the high dimensional issue, we propose to use a robust clustering algorithm that does not require dimension reduction prior to the clustering. We introduce new SH measures derived from the clustering and computed at the grassland level, to be consistent with ecological studies that usually estimate the biodiversity at the grassland level. The SH measures make possible the comparison of grasslands of varying sizes.
The proposed method is experimented on the Shannon index measured in 192 grasslands from southwest France with a dense intra-annual SPOT5 (Take5) multispectral time series and with single images extracted from this SITS. Note that contrary to [23], the spatial units are determined by the grasslands' spatial limits that are defined in a GIS, such as a land cover database.
In the next section, we present the materials used in this study. Then, the method proposed to measure the spectro-temporal heterogeneity is detailed in Section 3. Finally, the results are given in Section 4 and discussed in Section 5. Conclusions are given in Section 6.

Study Area
The study area is part of the Long-Term Ecological Research site "Coteaux et Vallées de Gascogne" (LTER_EU_FR_003), located in Gascony, in southwest France near the city of Toulouse (43 • Figure 1). This hilly area of around 900 km 2 is characterized by a mosaic of crops, small woods and grasslands. It is dominated by mixed crop-livestock farming. Grasslands provide food for cattle by grazing and/or producing hay or silage. They range from monospecific grasslands sown with ryegrass (improved with mineral fertilizing and mown up to three times a year) to semi-natural grasslands composed of spontaneous plant species (not fertilized and mown once a year). Grasslands are mainly located on steep slopes, whereas annual crops are in the valleys on the most productive lands. The climate is sub-Atlantic with sub-Mediterranean and mountain influences (mean annual temperature, 12.5 • C; mean annual precipitation, 750 mm) [27,28].

Satellite Image Time Series
The time series issued from the SPOT5 (Take5) mission (https://www.spot-take5.org) was used in this study. Eighteen images at a spatial resolution of 10 m were available from April-September 2015 ( Table 1).  The images were orthorectified, radiometrically and atmospherically corrected by the French Spatial Agency (CNES). They were provided in reflectance with a mask of clouds and shadows issued from the MACCS (Multi-sensor Atmospheric Correction and Cloud Screening) processor [29].
To reconstruct the time series due to missing data (clouds and their shadows), the Whittaker filter [30,31] was applied pixel-by-pixel on the reflectances in each spectral band. The smoothing parameter was the same for all the pixels and all the spectral bands. It was fixed to 10 4 after an ordinary cross-validation done on a subset of the pixels.
The smoothed time series associated with each of the spectral bands were concatenated to get a unique spectro-temporal vector x k per pixel k, such as where x kB1 (t j ) is the value of pixel k in band B1 at the j-th acquisition, and T = 18 is the number of acquisitions.

Field Data
Grasslands composing the dataset have been monitored for several years in the frame of different research projects. They represent more than 200 managed grasslands. The management practices and their intensity (i.e., number of mowings, intensity of grazing) are known for some grasslands.
The grasslands were digitalized in a GIS from aerial photographs ("BD ORTHO R " French database of orthophotos, c IGN). For this study, an inner buffer of 10 m was removed from all the grasslands' polygons to avoid edge effects due to mixed pixels at the parcel edges. After rasterizing the polygons, only the grasslands composed of at least 10 pixels of 10-m resolution, i.e., having an area higher than 1000 m 2 , were kept to ensure a minimum number of pixels per grassland. After this treatment, the dataset is composed of 192 grasslands. Their location can be seen in Figure 1.
A botanical survey was conducted in the spring of 2015 and 2016 after the flowering and before the mowing (April-May), to record the botanical composition of these grasslands. The grassland composition is supposed to remain stable from one year to the following year. The survey consisted of an exhaustive visual recording of all the species present in the grassland. The recording was processed while walking on a "W-shaped" transect, and the percentage of cover of each species was estimated at the grassland scale. The cover was estimated with the Braun-Blanquet abundance-dominance coefficients [32] for each present species (*: one individual, +: cover <1%, 1: 1-5%, 2: 5-25%, 3: 25-50%, 4: 50-75%, 5: 75-100%). An average abundance was kept for each coefficient (*: 0.1%, +: 0.2%, 1: 2.5%, 2: 15%, 3: 37.5%, 4: 62.5%, 5: 87.5%). From these absolute abundance-dominance covers, relative proportions of cover in the grassland can be retrieved for each species.
In this study, the species richness (number of species in the plot) was not considered because it accounts for rare species in the grassland, which might impact its functional diversity, but which do not have an impact on the spatio-spectral heterogeneity [18]. Therefore, an abundance-based biodiversity index was preferred to measure the alpha-diversity, the Shannon index (H): where p s is the proportion of the s-th species with ∑ R s=1 p s = 1 and R is the total number of species in the grassland (species richness). H values usually range between 0 and 5, increasing as the diversity increases. The Shannon index is a measure of the entropy in the grassland. It reflects the evenness of a population: a community with one or two dominating species is considered less diverse than a community that has different species with a similar number of individuals [18].
Most of the grasslands are semi-natural grasslands with a medium to high level of biodiversity (H > 2) (Figure 2a). Only a few are monospecific grasslands (sown with one species, H < 0.5). Three examples of grasslands' temporal profiles along the H axis, from a low to a high level of biodiversity, are shown in Figure 3.
The average grassland size in pixels is 135 pixels, and the median is 94 pixels (Figure 2b). In total, there are 25,903 pixels in the dataset.

Method
The objective of this part is to describe the proposed method to measure the spectrotemporal heterogeneity in the grasslands from SITS and to link it with the Shannon index measured in the field in order to verify the STVH.
Each pixel k is represented by a vector x k of dimension d, d being the number of variables. For instance, in the case of hyperspectral data, d is the number of spectral bands, which is a few hundred. In the case of multitemporal data, usually a vegetation index is used, and d   Table 3. The x-axis corresponds to the day of the year of 2015, and the y-axis corresponds to the NDVI. Grasslands have been voluntarily chosen by their high number of pixels for better visualization.
corresponds to the number of temporal measurements. In this study, where we use both the spectral and the temporal information, d = n B n T , with n B the number of spectral bands and n T the number of temporal acquisitions, as presented in Section 2.2.

Measures of Spectral Heterogeneity in the Literature
In the literature, the measure of spectral heterogeneity is based on measures of dispersion [8] such as the standard deviation [14] or the coefficient of variation [20]. However, these measures require selecting single bands or performing band reduction, such as using a vegetation index or using ordination methods like Principal Components Analysis (PCA), and thus, they lose some relevant information. To enable the use of all the spectral information, Rocchini [22] proposed the mean of the pairwise Euclidean distances from the spectral centroid (MDC) [15] for all the pixels covering the sampling plot: where n i is the number of pixels in the plot i, x ik is the spectral vector associated with pixel k, µ i is the plot's spectral centroid and · 2 stands for the Euclidean distance. In our case, the plot is the grassland. Hence, the centroid is the grassland's pixels' centroid, i.e., the mean spectro-temporal value of the grassland's pixels.
To reduce the dimensional-space, some studies compute the MDC on the first few components of PCA performed on the spectral variables [18,19,22]. Theoretically, it is almost equivalent to the original MDC.
MDC is in fact the trace of the pixels' empirical covariance matrix V i , which measures the spectral variability in the plot: and: However, several drawbacks can be raised from the MDC. This measure is flawed because it assumes homoscedasticity of the variables (same variance), and it does not differentiate between monomodal and multimodal distributions in the spectral space. Different spectral configurations can have the same distance to the centroid [23] as illustrated in a two-dimensional space in Figure 4.  To address this issue, Féret and Asner [23] introduced a clustering step to account for the global distribution in the spectral space. The unsupervised clustering is used to obtain "spectral species" (clusters), which are related to one or several species sharing similar spectral signatures. The clusters were estimated through a PCA and a k-means procedure. Although this clustering algorithm is often used, it is not robust in a high dimensional space [33]. Moreover, it assumes homoscedasticity of the clusters (same variance for each cluster), and it does not reflect well the distribution in the spectral space. An illustration can be seen in Figure 5 where (a) a simulated plot made of three different spectral species was clusterized by (b) the pipeline PCA + k-means and (c) by Gaussian mixture models. With PCA + k-means, the spectral species (clusters) are not well found contrary to the clustering using Gaussian mixture models. Thus, in the following, we propose a clustering technique that is robust to high dimensional data. Hence, no dimension reduction such as PCA is necessary. Moreover, it assumes heteroscedasticity of the clusters, i.e., each cluster could have a different variability. This clustering algorithm enables the computation of other measures of spectral heterogeneity.

Spectral Clustering Algorithm for High Dimensional Data and Derived Measures of Spectral Heterogeneity
We suggest to use a robust clustering algorithm that encompasses k-means, but that is suitable in a context of a small sample size with a large number of variables: the High Dimensional Data Clustering (HDDC) [34] (available in the R package at https://cran.r-project. org/web/packages/HDclassif/index.html). HDDC is based on a mixture model where each mixture component follows a Gaussian distribution. Under this model, a given sample (i.e, pixel) x is the realization of a random vector, for which distribution p is such that [35]: where C is the number of clusters, π c is the proportion of cluster c and f c (x|µ c , Σ c ) is a Gaussian distribution with parameters µ c and Σ c , i.e., The parameters of the mixture model are estimated using a conventional expectation-maximization algorithm. Once the optimal parameters are found, each sample x is assigned to the cluster c for which the log-probability Q c (x) is maximal. Q c (x) is computed as: The computation of Equation (7) requires the inversion of the covariance matrix Σ c and the computation of the logarithm of its determinant, which can be numerically unstable in a high dimensional context. To circumvent these issues, HDDC assumes that the last (lowest) eigenvalues of the covariance matrix are equal. It results that the inverse of the covariance matrix and its determinant can be computed explicitly while the numerical stability is controlled [36].
In this study, the clustering is applied to all the grasslands' pixels x k ∈ R d , regardless of the grassland to which they belong. The clustering splits all the pixels into C clusters. Then, for each grassland g i , its corresponding pixels x ik are assigned to C i clusters with k ∈ {1, ..., n i }; n i is the number of pixels in g i , C i ∈ {1, . . . , C} and C i C. For each cluster c in the grassland, the mean vector µ ic and the covariance matrix associated with the pixels belonging to this cluster are updated from the initial clustering on all the pixels. The proportion of this cluster p ic is also updated, p ic = n ic n i with n ic the number of pixels of g i associated with c. Hence, by considering several clusters inside a given grassland, it is possible to assess the between-class variability, the within-class variability and the entropy within the grassland. They provide additional information with respect to MDC to assess the spectral heterogeneity, and they are defined in the next subsections.

Between-and Within-Class Variabilities
The covariance matrix of g i can be decomposed as [37]: with: • µ ic is the spectro-temporal mean of pixels in g i assigned to cluster c, • µ i is the mean spectro-temporal value computed from all the pixels of g i , • V ic is the empirical covariance matrix of pixels of g i assigned to cluster c.
Therefore, using Equation (4), the MDC associated with g i can be written as: From Equation (9), two measures of spectral heterogeneity can be extracted: the betweenclass variability and the within-class variability. The trace of B i quantifies the between-class variability: This measure reflects the inter-classes variance: how the means of clusters in the grassland are different or similar. The more the clusters are different, the higher the between-class variability. If there is only one class in the grassland, then trace(B i ) = 0.
The within-class variability is quantified by the trace of W i : This measure represents the mean of the clusters variances. The more the grassland has heterogeneous clusters with high variance, the higher the within-class variability. However, if the grassland has many homogeneous clusters, trace(W i ) will be low. Contrary to the betweenclass variability, if there is only one cluster, the within-class variability still provides information on the heterogeneity in the grassland.

Entropy
Another measure of spectral heterogeneity that can be derived from the spectral clustering of grasslands is the entropy. The entropy is linked to the proportions of clusters within the grassland g i and is quantified by: The more the dataset has equally-balanced clusters, the higher E i . The least balanced is the dataset, the closer E i is to 0. If there is only one cluster in the grassland, C i = 1 and p i1 = 1, then E i = 0.
The HDDC algorithm provides for each pixel the probability that it belongs to each cluster, which can be understood as a soft assignment with respect to the hard assignment of k-means ( Figure 6). Therefore, a "finer" measure of entropy can be computed using the probability of assignment of the grassland's pixels to each cluster c of the global clustering, π ic : with π ick the assignment probability of pixel k from g i to cluster c provided by the algorithm, ∑ C c=1 π ick = 1 ( Figure 6). Therefore, a finer measure of entropy can be written by replacing p ic by π ic in Equation (12) and by summing it to the total number of clusters C: This measure of entropy hardly ever reaches null values, unless all the pixels of the grasslands are assigned to the same cluster with a probability of 1. The entropy reflects the grassland's clusters evenness: whether it is dominated by one cluster or numerous equally-distributed clusters.

Methodology
To link the proposed SH measures issued from SITS to the Shannon index measured from the species, univariate and multivariate (combining several SH measures), linear regressions are performed. The response variable is the Shannon index, and the explanatory variables are the global variability or MDC (Equation (9)), the between-class variability (Equation (10)), the within-class variability (Equation (11)) and the entropy with soft assignment (Equation (14)).
Since the linear regressions assume normality of the distributions, the global variability, the between-class variability and the within-class variability are log-transformed to Gaussianize them [38], as done in [19,39]. In the following, the entropy with soft assignment is denoted by E, and the log-transformed global (or MDC), within-class and between-class variabilities are denoted by V, W and B, respectively.
The adjusted coefficient of determinationR 2 is used to measure the goodness of fit of the regressions. It is defined as the proportion of variance explained by the regression model adjusted for the number of explanatory variables.
The proposed methodology including the clustering is synthesized in Figure 7. To assess the contribution of temporal variations to the SVH through the use of multitemporal data, we also applied the same methodology using only one acquisition issued from the SITS. We compared the results obtained with the SITS and obtained with a single image by computing a Wilcoxon signed-rank test between the two distributions of predicted values issued from their best models.
The clustering algorithm, the computation of the SH measures and the statistical analysis were performed in Python through the SciPy (https://www.scipy.org), scikit-learn [40

Results
The proposed SH measures were computed using the spectro-temporal data for all the grasslands for different numbers of clusters. Then, the adjusted coefficient of determination of the linear regression between the Shannon index measured in the field (H) and the individual or combined SH measures was calculated (Figure 9). The entropy computed with soft assignment (Equation (14)) was slightly better correlated with H than the "simple" entropy (Equation (12)); therefore, the usual entropy is not shown.

Univariate Correlation with Multitemporal Data
The global variability, or MDC computed at the grassland scale, is significantly correlated with the Shannon index (R 2 = 0.071, p-value <0.001). Depending on the number of clusters, the entropy and the within-class variability reach higher correlation coefficients than the global variability. For instance, for the entropy, with C = 8 and C = 75, the adjusted coefficient of determination isR 2 = 0.099 andR 2 = 0.105, respectively (p-value <0.001). Its minimum is R 2 = −0.005 and not significant for C = 2. The within-class variability reaches maximum correlation values ofR 2 = 0.097 for C = 2 andR 2 = 0.074 for C = 6 (p-value <0.001). Its minimum value isR 2 = 0.034 for C = 75 (p-value <0.05). The between-class variability never reaches as high values as MDC except for C = 14 (R 2 = 0.072, p-value < 0.001). The problem with this measure lies with its null values, which make it not continuous.
The relationships between H and the proposed SH measures for each grassland that reach the highest coefficients of determination are shown in Figure 10. The entropy is the SH measure showing the highest coefficient of determination with H. All the SH measures tend to increase with the Shannon index suggesting that the spectro-temporal heterogeneity is linked to the species diversity.

Multivariate Correlation with Multitemporal Data
Multivariate linear regressions were run with different combinations of SH measures to assess which combinations of variables are the most related to H.
The models combining several variables explain better the Shannon index than the univariate models (Figure 9, blue and red lines). Indeed, the multivariate model combining the  Table 2). Therefore, the combination of W and E (R 2 = 0.131) explains H better than the combination of W, B, V and E (R 2 = 0.127, Figure 9 and Table 2). This is the case for most of the numbers of clusters. Additionally, the combination of V and E (data not shown) is worse to explain H than the combination of W and E.

Univariate and Multivariate Correlation with Monotemporal Data
To evaluate the contribution of multitemporal data to the SVH, we compared the above results to results obtained from monotemporal data. We chose two acquisitions dates from the time series: 30 April (near the growth peak, before the occurrence of the management practices such as mowing and grazing) and 29 June (after most of the management practices occurred).
Higher coefficients of correlations are obtained with only one acquisition ( Figure 11). Using the image of 30 April, the maximum adjusted coefficient of determination is 0.139 (p-value <0.001) with the model combining W and E for C = 150 andR 2 = 0.169 (p-value <0.001) with the model combining W, B, V and E for C = 150 (for higher numbers of clusters,R 2 was lower, data not shown). Using the image of 29 June,R 2 = 0.137 with the model combining W and E, andR 2 = 0.140 (p-value <0.001) with the model combining W, B, V and E, both for C = 20.
For both images, the combination of the four proposed SH and the combination of W and E are better at explaining the Shannon index than MDC. However, the contributions of each SH measure in the model are not the same as for the model using multitemporal data.
We compared the distributions of the predicted Shannon index values issued from the best models using the three types of data (i.e., multitemporal data: C = 8, explanatory variables: There were no significant differences between the values predicted by the best models for each data type (p-value >0.05 for each pair of distributions). Therefore, the models are equivalent in terms of predicted H values.

Spectral Heterogeneity Measures
Globally, the variance of the species diversity (represented by the Shannon index) explained by the SH measures derived from the clusterings of grasslands using SITS is weak.
Indeed, low species diversity grasslands are not associated with low SH measures, except for E (Figure 10). High species diversity grasslands can be associated with any SH measures: low to medium E values, medium V values, high W values and medium B values. Medium species diversity grasslands, that represent a great percentage of the dataset, can be associated with any SH measures. As a result, the unexplained variance of H is very high (87%) with the model combining E and W (C = 8). Indeed, it predicts a much smaller range of values than the actual ones (Figure 2a): the maximum predicted H is 2.8, while the minimum is 1.8 with a mean of 2.3 and a standard deviation 0.2.
The spatial, but also the spectral resolution of the sensor may have limited the analysis. Indeed, individual grassland's species are particularly small and mixed. A spatial resolution of 10 m is too coarse to detect individual species. In a pixel of 10 m, there can be a large number of mixed grassland plant species. Moreover, some species are spectrally too similar to be discriminated with low spectral resolution [42,43]. Consequently, if a grassland has a high level of biodiversity, with a large number of species, and if these species are homogeneously mixed within the grasslands, the pixels will be spectrally similar even if they contain a mix of species. Thus, there would be one cluster in the grassland. This would result in a low spectral entropy, although this grassland has a high Shannon index. This could explain the wide range of SH values associated with grasslands with high species diversity. These SH measures seem to reflect more the variability of homogeneous species assemblages (i.e., the clusters) in the grasslands than the diversity of species, explaining the low, but significant relationship with the ground Shannon index. Despite the weak relationship with the Shannon index, we proposed SH measures that provide supplementary information on the grassland's heterogeneity with regard to MDC. Indeed, the combination of the entropy and the within-class variability was always more correlated with the species diversity than the MDC alone, regardless of the type of imagery used. These two SH measures contributed the most to explain H.
The within-class variability is an interesting measure since it provides a quantitative information on the grassland's heterogeneity even if there is only one cluster, which cannot be provided by the entropy. Thus, the within-class variability and the entropy are complementary. Their combination was better than the combination of the four proposed SH measures together, or the univariate models, regardless of the number of clusters (Figures 9 and 11, blue line).
Beyond the relationship with species diversity of the SH measures, the proposed method makes possible the detection of assemblages of species within the grasslands, which share common spectro-temporal properties. These assemblages can indirectly give information about the heterogeneity within the grasslands. The heterogeneity of these groups of species can be quantified with their spectral variability ( Figure 12). Such a spectral map at the grassland scale would not be possible using MDC.  Table 3.

Clustering
To understand the meaning of the clusters found using multitemporal data (C = 8), the mean vectors corresponding to each cluster were extracted, and their NDVI temporal profiles were computed (Figure 13). The profiles seem consistent with typical profiles of grasslands. The pixels associated with cluster C7 in Figure 13 belong to grasslands varying in species diversity, but all intensively-used. These pixels usually represent the whole grasslands, while grasslands less intensively used can be associated with several other clusters.  Figure 13: Mean NDVI temporal profiles of each cluster from the clustering into C = 8 clusters using multitemporal data. The x-axis is the month of year 2015, and the y-axis is the NDVI.
Hence, the clusters seem to be more related to phenological profiles linked to management practices than to phenological profiles of species. Indeed, the management practices have an influence on the species distribution and composition [44], but they may have a stronger impact on the spatial, temporal and spectral profiles of the grasslands because they induce abrupt changes in the grassland (mowing, grazing, fertilizing). In particular, due to the use of acquisitions from April-September, the effect of management practices that usually occur within this period may be very significant.
More precisely, we suspect that clusters are related to the intensity of practices. Indeed, an intensive use with constant defoliation does not allow for the expression of the phenology of species. However, when the grassland is extensively used, species can express different phenologies (during the regrowth after the mowing for instance), and different clusters related to these phenologies can be detected. This could explain the multiple clusters found in extensively-used grasslands (for instance, in Figure 12c), while intensively-used grasslands are represented by one cluster, which has a typical signature of intensively-used grassland ( Figure 12a).
Hence, at a spatial resolution of 10 m, the clusters found using multitemporal data seem to reflect more the intensity of practices in the grasslands than the species composition. This could explain the weak correlations with the Shannon index.
Regarding the number of clusters, the proposed clustering algorithm (HDDC) provides model selection criteria (ICL and BIC) that were not efficient in our experiment. Indeed, the theoretical optimal number of spectral clusters may not correspond to the number of expected clusters of species [45]. Therefore, our strategy was to test a wide range of numbers of clusters and to keep the one that gives the bestR 2 . From an operational viewpoint, this strategy can be time consuming, but it can adapt to any spatial configurations (size and location).

Contribution of Multitemporal Imagery
In this study, we made the assumption that the spectro-temporal variations of a grassland could be related to its species diversity. The results obtained with the monotemporal imagery showed that the multitemporal data do not improve the relationship with the Shannon index. Indeed, higher coefficients of determination were reached with the two dates proposed than with the full SITS. Hence, the SVH with temporal variations, the so-called STVH, is not verified in this work, at a spatial resolution of 10 m.
However, the clusters found in the grasslands using the full SITS create homogeneous patterns within the grasslands, contrary to the clusters found with one image, which are quite "pixelized" and do not seem spatially consistent ( Figure 14). Considering that the predicted values with models issued from these three datasets were not significantly different, we can also doubt the relationship of the clusters found using one image with the species diversity. However, this would require verification in the field, but this was not possible in the frame of this work. As previously suggested, the temporal variations measured by the sensor seem to be more related to the management practices than to the species diversity. Indeed, we used a time series covering the period from mid-April-September. Most of the management practices such as mowing and/or grazing usually occur within this period. To circumvent this effect, the time series could be limited to a period or a combination of periods when there is no management practices, such as the beginning of the growing season, before the growth peak. Moreover, previous studies have shown that the relationships between species diversity and remote sensing metrics can be season-dependent [46]. Therefore, specific care should be considered regarding the dates of the imagery selected.

Limitations
The weak relationships found in this work can be due to the unbalanced H values present in our dataset. Indeed, it was mostly composed of grasslands with a medium level of biodiversity (H between two and 2.5, Figure 2a). The H gradient was not very well sampled, with very few grasslands low in species diversity (H < 1) and rich in species diversity (H > 3). Hence, the regression models were more calibrated on medium level biodiversity and lacked generality. This may be why the models have an average predicted H value of around 2.3 (Section 5.2).
Moreover, we obtained lower coefficients of determination than in other studies that related the species diversity in grasslands with the SH using monotemporal imagery at a very high spatial resolution. For instance, Oldeland et al. [18] used airborne hyperspectral data at a spatial resolution of 5 m and found significant correlations between the Shannon index of savannah plots and the MDC computed from PCA. They reached significant R 2 ranging from 0.31-0.62 depending on the 20 m × 50 m plots. Möckel et al. [19] investigated the prediction of grassland species diversity (species richness and inverse Simpson's diversity) in Sweden from airborne hyperspectral data with a spectral response approach and a spectral heterogeneity approach. However, they failed to detect a significant relationship between species diversity and spectral heterogeneity (PCA + MDC) at the plot scale, contrary to the spectral response approach. In the same study area, Hall et al. [39] related the species richness (alpha-diversity) and the species turnover (beta-diversity) measured in three plots per grassland with the spectral heterogeneity in the four bands of the QuickBird sensor (2.4 m resolution) and other field variables. The spectral heterogeneity was measured as the mean difference between the mean of each individual 3 × 3 pixel window (corresponding to each plot) and the mean of all three pixels windows within each grassland site. It can be assimilated to the between-class variability, but with three plots of the same size. They found low, but significant linear correlations between the species richness and the spectral heterogeneity measured with the NIR band (R 2 = 0.08) and between the species turnover and the spectral heterogeneity measured in the red band (R 2 = 0.10), the NIR band (R 2 = 0.19) and the NDVI (R 2 = 0.14). Better correlations were found with multivariate models, but only the model predicting the species turnover included the spectral heterogeneity (NIR,R 2 = 0.33).
However, these studies were conducted at the plot scale, both for the floristic record and the associated spectral information. They used the pixels corresponding only to the sampling unit. Our protocol was different, since the botanical survey was conducted at the grassland scale by a random walk strategy, and only one biodiversity index was computed from it. Yet, grasslands are characterized by patterns of small scale species composition and spatial distribution [47][48][49]. Hence, this estimation of the biodiversity at the grassland level may be difficult to relate to remote sensing data and might have limited our analysis.
Furthermore, the influence of the topography should be considered for future studies because it is known that the topography influences the reflectance.

Outlooks
In terms of methodology, the proposed method could be used to assess the beta-diversity among grasslands. Indeed, in light of the Bray-Curtis dissimilarity [50], a spectral dissimilarity could easily be computed from the proportions of clusters in each grassland [23]. It could be improved by using the probability of belonging to each cluster, similarly to the way done with the entropy. The pairwise spectral Bray-Curtis dissimilarity between grasslands g i and g j would be defined as: where π ic and π jc are the mean assignment probabilities to cluster c of pixels of grasslands g i and g j , respectively, defined in Equation (13).
In terms of application, the proposed method is not specific to grasslands, and it could be used to assess the species diversity of other habitats. For instance, it could be used on forest, since the method is not required to work at a specific object scale: it can be applied to a plot of fixed size.
In addition, this work could be extended to the relationship of the spectral heterogeneity with the functional diversity of the habitat. Indeed, some functional traits are related to the way plants reflect light and thus to the signal measured by the sensor [9,[51][52][53] and may be related to the spectral heterogeneity. However to our knowledge, functional diversity has not yet been related to remotely-sensed measures [53,54] and has not been discussed in the context of SVH [8,10]. The stakes would be to determine which traits and which measures of functional diversity are the most consistent with SVH. Using SITS, we would suggest to select functional traits that are linked to the phenology of the species such as the flowering date, the flowering length and the leaf life span.

Conclusions
The aim of this work was to attempt to verify the Spectral Variation Hypothesis (SVH) in grasslands under the assumption that the temporal variations could be used in addition to the spectral variations of the habitat as a proxy of its species diversity: the Spectro-Temporal Variation Hypothesis (STVH). To do so, we proposed a method based on an unsupervised clustering of the grasslands using multitemporal and multispectral data, allowing for the derivation of spectro-temporal heterogeneity measures computed at the grassland level: the within-class variability, the between-class variability and the entropy. We compared them to the commonlyused mean distance to the centroid. The method was applied on 192 grasslands from southwest France using an inter-annual multispectral time series of SPOT5 images. Univariate and multivariate regression models combining several spectro-temporal heterogeneity measures were run with different numbers of clusters to assess their correlation with the Shannon index measured from field data.
The tested spectral heterogeneity measures were found significantly, but weakly correlated with the Shannon index. The combination of the within-class variability and the entropy was found always better correlated with the Shannon index than the mean distance to the centroid, regardless of the number of clusters. The best regression model explained 13.1% of the variance of the ground Shannon index while the mean distance to the centroid explained 7.1% of the variance. Hence, the clustering makes possible the extraction of spectral heterogeneity measures that give supplementary information to the mean distance to the centroid. However, equivalent results were obtained using monotemporal imagery.
Therefore, the spectro-temporal variation hypothesis was not verified using multispectral multitemporal imagery at a spatial resolution of 10 m. The proposed spectro-temporal heterogeneity measures seemed to be more related to the management practices performed in the grasslands than to the species diversity. The use of a whole time series covering the growing season or the season when the management practices occur does not seem to be suitable to detect the diversity in species. A period when no practice occurs should be more appropriate.
More research should be conducted on the extension of the SVH to the functional diversity. The STVH might be more related to functional traits linked to the phenology of species. Table 3: Braun-Blanquet abundance-dominance coefficients associated with each plant species recorded in three grasslands a, b and c having a Shannon index of 0.10, 1.57 and 2.89, respectively. "spp." means that the species from the given genus was not identified.