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

Mailys Lopes, Mathieu Fauvel, Annie Ouin, Stéphane Girard

► To cite this version:

Mailys Lopes, Mathieu Fauvel, Annie Ouin, Stéphane Girard. Spectro-Temporal Heterogeneity Measures from Dense High Spatial Resolution Satellite Image Time Series: Application to Grassland Species Diversity Estimation. Remote Sensing, 2017, 9 (10), pp.993:1-23. 10.3390/rs9100993. hal-01613722

HAL Id: hal-01613722

https://hal.science/hal-01613722

Submitted on 10 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spectro-Temporal Heterogeneity Measures from Dense High Spatial Resolution Satellite Image Time Series: Application to Grassland Species Diversity Estimation

Published in *Remote Sensing* 2017, 9(10), 993; doi:10.3390/rs9100993.

Mailys Lopes¹, Mathieu Fauvel¹, Annie Ouin¹ and Stéphane Girard²

¹ Dynafor, University of Toulouse, INRA, INPT, INPT-EI PURPAN, Chemin de Bordeaux, 31326 Castanet Tolosan, France

² Team Mistis, Inria Grenoble Rhône-Alpes, LJK, 655 Avenue de l'Europe, 38334 Montbonnot, France

## Abstract

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.

**Keywords:** Spectral Variation Hypothesis; Spectral heterogeneity; Dense satellite image time series; Alpha-diversity; Grasslands.

## 1 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.
# 2 Materials

## 2.1 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°17′N, 0°54′E, Figure 1). This hilly area of around 900 km² 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].

<figure><img src="image_2.png" /><figcaption>Figure 1: Location of the study area in southwest France and of the grasslands within the study area. The background is an aerial photograph issued from the French orthophoto database "BD ORTHO®" (©IGN).</figcaption></figure>
## 2.2 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).

<table><caption>Table 1: Characteristics of SPOT5 (Take5) imagery used in this study.</caption><tr><td>Pixel Size</td><td>10 m</td></tr><tr><td>Spectral bands</td><td>B1 "Green" (500–590 nm)<br>B2 "Red" (610–680 nm)<br>B3 "Near-Infrared" (780–890 nm)<br>B4 "Short Wave Infrared" (1580–1750 nm)</td></tr><tr><td>Acquisition dates</td><td>20-04-2015, 25-04-2015, 30-04-2015, 10-05-2015, 20-05-2015, 04-06-2015,<br>24-06-2015, 29-06-2015, 04-07-2015, 09-07-2015, 14-07-2015, 19-07-2015,<br>24-07-2015, 13-08-2015, 18-08-2015, 28-08-2015, 02-09-2015, 07-09-2015</td></tr></table>

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 Whitaker 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 $\mathbf{x}_k$ per pixel $k$, such as

$$ \mathbf{x}_k = [\mathbf{x}_{kB_1}(t_1), \dots, \mathbf{x}_{kB_1}(t_T), \mathbf{x}_{kB_2}(t_1), \dots, \mathbf{x}_{kB_2}(t_T), \mathbf{x}_{kB_3}(t_1), \dots, \mathbf{x}_{kB_3}(t_T), \mathbf{x}_{kB_4}(t_1), \dots, \mathbf{x}_{kB_4}(t_T)]^\top, $$

where $\mathbf{x}_{kB_1}(t_j)$ is the value of pixel $k$ in band B1 at the $j$-th acquisition, and $T = 18$ is the number of acquisitions.

## 2.3 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®" French database of orthophotos, ©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², 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):

$$H = - \sum_{s=1}^{R} p_s \ln p_s \qquad (1)$$

where $p_s$ is the proportion of the s-th species with $\sum_{s=1}^{R} 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.

## 3 Method

The objective of this part is to describe the proposed method to measure the spectro-temporal 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
<figure><img src="image_3.png" /><figcaption>Figure 2: Histogram of (a) Shannon index H and (b) grasslands' size in number of pixels $n_i$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: SPOT5 NDVI temporal profiles of all the pixels belonging to three grasslands along the H gradient: (a) grassland with a low level of biodiversity, (b) grassland with a medium level of biodiversity and (c) grassland with a high level of biodiversity. The floristic record of these three grasslands can be found in Appendix ??, 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.</figcaption></figure>

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.

## 3.1 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:

$$MDC_i = \frac{1}{n_i} \sum_{k=1}^{n_i} \|\mathbf{x}_{ik} - \mathbf{\mu}_i\|^2 \quad (2)$$

where $n_i$ is the number of pixels in the plot $i$, $\mathbf{x}_{ik}$ is the spectral vector associated with pixel $k$, $\mathbf{\mu}_i$ is the plot's spectral centroid and $\| \cdot \|^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:

$$V_i = \frac{1}{n_i} \sum_{k=1}^{n_i} (\mathbf{x}_{ik} - \mathbf{\mu}_i)(\mathbf{x}_{ik} - \mathbf{\mu}_i)^{\top} \quad (3)$$

and:

$$\text{trace}(V_i) = \frac{1}{n_i} \sum_{k=1}^{n_i} \|\mathbf{x}_{ik} - \mathbf{\mu}_i\|^2. \quad (4)$$

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.

<figure><img src="image_5.png" /><figcaption>Figure 4: Simulated pixels' distributions for two different plots (a) and (b). Pixels are displayed in blue, and the centroids of the plots are displayed in red. The estimated MDC are very close while the spectral distributions of the plots are clearly different.</figcaption></figure>

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.

<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Simulated distributions of three spectral species (blue, yellow, red) in a 3-dimensional space. (b) Clustering of the three distributions with PCA and k-means. (c) Clustering of the three distributions with Gaussian mixture models.</figcaption></figure>

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.

## 3.2 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]:

$$p(\mathbf{x}) = \sum_{c=1}^{C} \pi_c f_c(\mathbf{x} | \mu_c, \Sigma_c) \qquad (5)$$
where C is the number of clusters, $\pi_c$ is the proportion of cluster c and $f_c(x|\mu_c, \Sigma_c)$ is a Gaussian distribution with parameters $\mu_c$ and $\Sigma_c$, i.e.,

$$f_c(x|\mu_c, \Sigma_c) = \frac{1}{(2\pi)^{\frac{d}{2}}|\Sigma_c|^{\frac{1}{2}}} \exp \left( -\frac{1}{2}(x - \mu_c)^{\top}\Sigma_c^{-1}(x - \mu_c) \right). \quad (6)$$

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:

$$Q_c(x) = \frac{1}{2} \times \left( -(x - \mu_c)^{\top} \Sigma_c^{-1} (x - \mu_c) - \log(|\Sigma_c|) + 2 \log(\pi_c) - d \log(2\pi) \right). \quad (7)$$

The computation of Equation (7) requires the inversion of the covariance matrix $\Sigma_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 \in \mathbb{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 \in \{1, ..., n_i\}$; $n_i$ is the number of pixels in $g_i$, $C_i \in \{1, ..., C\}$ and $C_i \ll C$.

For each cluster c in the grassland, the mean vector $\mu_{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} = \frac{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.

### 3.2.1 Between- and Within-Class Variabilities

The covariance matrix of $g_i$ can be decomposed as [37]:

$$V_i = B_i + W_i \quad (8)$$

with:

* $B_i = \sum_{c=1}^{C_i} p_{ic}(\mu_{ic} - \mu_i)(\mu_{ic} - \mu_i)^{\top}$ is the between-class covariance matrix,

* $\mu_{ic}$ is the spectro-temporal mean of pixels in $g_i$ assigned to cluster c,

* $\mu_i$ is the mean spectro-temporal value computed from all the pixels of $g_i$,

* $W_i = \frac{1}{n_i} \sum_{c=1}^{C_i} \sum_{k \in c} (x_{ik} - \mu_{ic})(x_{ik} - \mu_{ic})^{\top} = \sum_{c=1}^{C_i} p_{ic} V_{ic}$ is the within-class covariance matrix,

* $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:

$$ MDC_i = \operatorname{trace}(\mathbf{V}_i) = \operatorname{trace}(\mathbf{B}_i) + \operatorname{trace}(\mathbf{W}_i). \quad (9) $$

From Equation (9), two measures of spectral heterogeneity can be extracted: the between-class variability and the within-class variability. The trace of $\mathbf{B}_i$ quantifies the between-class variability:

$$ \operatorname{trace}(\mathbf{B}_i) = \sum_{c=1}^{C_i} p_{ic} \|\boldsymbol{\mu}_{ic} - \boldsymbol{\mu}_i\|^2. \quad (10) $$

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 $\operatorname{trace}(\mathbf{B}_i) = 0$.

The within-class variability is quantified by the trace of $\mathbf{W}_i$:

$$ \operatorname{trace}(\mathbf{W}_i) = \frac{1}{n_i} \sum_{c=1}^{C_i} \sum_{k \in c} \|\mathbf{x}_{ik} - \boldsymbol{\mu}_{ic}\|^2. \quad (11) $$

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, $\operatorname{trace}(\mathbf{W}_i)$ will be low. Contrary to the between-class variability, if there is only one cluster, the within-class variability still provides information on the heterogeneity in the grassland.

### 3.2.2 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:

$$ E_i = \sum_{c=1}^{C_i} -p_{ic} \log(p_{ic}). \quad (12) $$

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, $\pi_{ic}$:

$$ \pi_{ic} = \frac{1}{n_i} \sum_{k \in g_i} \pi_{ick} \quad (13) $$

with $\pi_{ick}$ the assignment probability of pixel k from $g_i$ to cluster c provided by the algorithm, $\sum_{c=1}^{C} \pi_{ick} = 1$ (Figure 6). Therefore, a finer measure of entropy can be written by replacing $p_{ic}$
by $\pi_{ic}$ in Equation (12) and by summing it to the total number of clusters C:

$$E_{s_i} = \sum_{c=1}^{C} -\pi_{ic} \log(\pi_{ic}). \qquad (14)$$

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.

<figure><img src="image_7.png" /><figcaption>Figure 6: Grassland clustered with an initial clustering of the landscape into 8 clusters. (a) Hard assignment of the pixels. One color corresponds to one cluster (orange, yellow, blue). (b–d) Soft assignment of the pixels. The grey-scaled color corresponds to the assignment probability $\pi_{ick}$ to cluster (b) orange, (c) yellow and (d) blue.</figcaption></figure>

The entropy reflects the grassland's clusters evenness: whether it is dominated by one cluster or numerous equally-distributed clusters.

### 3.3 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 determination $\bar{R}^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] and
<figure><img src="image_8.png" /><figcaption>Figure 7: Overview of the method to compare the Spectral Heterogeneity (SH) measures (explanatory variables) to the Shannon index (response variable). Square rectangles correspond to data, and rounded rectangles correspond to a process.</figcaption></figure>

pandas (http://pandas.pydata.org) libraries.

The initial number of clusters for the clustering has an influence on the clusters found in each grassland (Figure 8). Hence, the correlations are studied for different numbers of initial clusters, from 2–150 clusters (every 2 clusters in the range [2, 60] and every 25 clusters in the range [75, 150]). For each number of clusters, 10 runs of the algorithm with different random initializations are performed, and the best result in terms of the Integrated Classification Likelihood (ICL) [41] is kept.

<figure><img src="image_9.png" /><figcaption>Figure 8: (a) False color image of a grassland acquired on 30 April 2015. The same grassland clustered using HDDC on multitemporal data with an initial clustering into (b) 8 clusters and (c) 150 clusters. Each cluster is represented by one color.</figcaption></figure>

# 4 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.

<figure><img src="image_10.png" /><figcaption>Figure 9: Adjusted coefficient of determination in the multivariate linear regression between different combinations of SH measures (V: log-transformed global variability or MDC, W: log-transformed within-class variability, B: log-transformed between-class variability, E: entropy) computed from multitemporal data and the Shannon index (response variable) depending on the number of clusters.</figcaption></figure>

## 4.1 Univariate Correlation with Multitemporal Data

The global variability, or MDC computed at the grassland scale, is significantly correlated with the Shannon index ($\bar{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 is $\bar{R}^2 = 0.099$ and $\bar{R}^2 = 0.105$, respectively (*p*-value <0.001). Its minimum is $\bar{R}^2 = -0.005$ and not significant for *C* = 2. The within-class variability reaches maximum correlation values of $\bar{R}^2 = 0.097$ for *C* = 2 and $\bar{R}^2 = 0.074$ for *C* = 6 (*p*-value <0.001). Its minimum value is $\bar{R}^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 ($\bar{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.

## 4.2 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
<figure><img src="image_11.png" /></figure>

Figure 10: Shannon index (H) best univariate linear correlations with different SH measures (E: entropy, V: log-transformed global variability or MDC, W: log-transformed within-class variability, B: log-transformed between-class variability) computed from multitemporal data. C is the corresponding number of clusters, $\bar{R}^2$ is the adjusted coefficient of determination and ** signifies *p*-value <0.001. The black line is the linear regression line.

entropy (E) and the within-class variability (W) and the model combining the four proposed SH measures (E, W, B and V) always has higher coefficients of determination than SH measures alone, regardless of the number of clusters. Moreover, for the number of clusters reaching the maximum adjusted coefficient of determination (C = 8), the variables contributing the most to explain H are W and E (*p*-value <0.05), whereas V and B do not contribute much (*p*-value >0.05, Table 2). Therefore, the combination of W and E ($\bar{R}^2$ = 0.131) explains H better than the combination of W, B, V and E ($\bar{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.

## 4.3 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 and $\bar{R}^2$ = 0.169 (*p*-value <0.001) with the model combining W, B, V and E for C = 150 (for higher numbers of clusters, $\bar{R}^2$ was lower, data not shown). Using the image of 29 June, $\bar{R}^2$ = 0.137 with the model combining W and E, and $\bar{R}^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:
<table><caption>Table 2: Multivariate linear models for C = 8 to explain the Shannon index (H) from the SH measures (V: log-transformed global variability or MDC, W: log-transformed within-class variability, B: log-transformed between-class variability, E: entropy) computed from multitemporal data. Reg. Coeff. is the regression coefficient; Std Err. is the standard error; F stands for the F-value with degrees of freedom in brackets; R² is the coefficient of determination; and R̄² is the adjusted coefficient of determination.</caption><thead><tr><th>Response Variable</th><th>Explanatory Variables</th><th>Reg. Coeff.</th><th>Std Err.</th><th>p-Value</th></tr></thead><tbody><tr><td rowspan="5">H</td><td>W</td><td>0.29</td><td>0.14</td><td>0.04</td></tr><tr><td>B</td><td>0.01</td><td>0.02</td><td>0.61</td></tr><tr><td>V</td><td>-0.15</td><td>0.14</td><td>0.30</td></tr><tr><td>E</td><td>0.40</td><td>0.13</td><td>0.003</td></tr><tr><td>intercept</td><td>0.73</td><td>0.51</td><td>0.16</td></tr><tr><td colspan="5">Model summary: F<sub>(4, 187)</sub> = 8.0, p-value &lt;0.001, R<sup>2</sup> = 0.145, R̄<sup>2</sup> = 0.127</td></tr><tr><td rowspan="3">H</td><td>W</td><td>0.16</td><td>0.06</td><td>0.005</td></tr><tr><td>E</td><td>0.37</td><td>0.09</td><td>&lt;0.001</td></tr><tr><td>intercept</td><td>0.65</td><td>0.51</td><td>0.20</td></tr><tr><td colspan="5">Model summary: F<sub>(2, 189)</sub> = 15.4, p-value &lt;0.001, R<sup>2</sup> = 0.140, R̄<sup>2</sup> = 0.131</td></tr></tbody></table>

W and E; 30 April image: C = 150, explanatory variables: W, B, V and E; and 29 June image: C = 20, explanatory variables: W, B, V and E) by conducting a Wilcoxon signed-rank test. 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.

# 5 Discussion

## 5.1 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.
<figure><img src="image_12.png" /><figcaption>Figure 11: Adjusted coefficient of determination in the multivariate linear regression using one image acquired on (a) 30 April and (b) 29 June between different combinations of SH measures (V: log-transformed global variability or MDC, W: log-transformed within-class variability, B: log-transformed between-class variability, E: entropy) and the Shannon index (response variable) depending on the number of clusters.</figcaption></figure>

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 cor- related 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 in-
formation 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 complemen-
tary. 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.

<figure><img src="image_13.png" /><figcaption>Figure 12: Maps of spectral heterogeneity inside three grasslands (a–c). The first row shows the grasslands’ polygon limits in yellow on the SPOT5 false color image acquired on 10 May 2015. The second row shows the clusters after an HDDC clustering into eight clusters using multitemporal data. The color scale corresponds to the log-transformed variability of each cluster *c* in the grassland *g*<sub>i</sub>. (a) H = 0.10, E = 0, V = 10.13, W = 10.13, B = 0; (b) H = 1.57, E = 0.68, V = 10.06, W = 9.41, B = 9.33; (c) H = 2.89, E = 1.06, V = 9.58, W = 9.22, B = 8.42. The floristic record of these three grasslands can be found in the Appendix ??, Table 3.</figcaption></figure>

**5.2 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><img src="image_14.png" /><figcaption>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.</figcaption></figure>

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 best $\bar{R}^2$. From an operational viewpoint, this strategy can be time consuming, but it can adapt to any spatial configurations (size and location).
5.3 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.

<figure><img src="image_15.png" /><figcaption>Figure 14: Clustering of the same grassland (false color image of (a) 30 April and (d) 29 June) with an initial clustering into 150 clusters, using (b) the image of 30 April and (c) the full SITS, and into 20 clusters, using (e) the image of 29 June and (f) the full SITS.</figcaption></figure>

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.

5.4 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 re-
lated 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² 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 spec-
tral 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 cor-
relations between the species richness and the spectral heterogeneity measured with the NIR
band (R² = 0.08) and between the species turnover and the spectral heterogeneity measured
in the red band (R² = 0.10), the NIR band (R² = 0.19) and the NDVI (R² = 0.14). Better
correlations were found with multivariate models, but only the model predicting the species
turnover included the spectral heterogeneity (NIR, R̄² = 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–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 be-
cause it is known that the topography influences the reflectance.
## 5.5 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:

$$BC_{ij_{\text{spectral}}} = \frac{\sum_{c=1}^{C} |\pi_{ic} - \pi_{jc}|}{\sum_{c=1}^{C} \pi_{ic} + \pi_{jc}} \quad (15)$$

where $\pi_{ic}$ and $\pi_{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–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.

## 6 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 commonly-used 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.

# Appendix
<table><caption>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.</caption><thead><tr><th>Species</th><th>a</th><th>b</th><th>c</th><th>Species</th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td>Agrimonia eupatoria</td><td></td><td></td><td>+</td><td>Medicago spp.</td><td>+</td><td></td><td></td></tr><tr><td>Agrostis capillaris</td><td></td><td>1</td><td></td><td>Muscari comosum</td><td></td><td></td><td>+</td></tr><tr><td>Anthoxanthum odoratum</td><td></td><td>1</td><td></td><td>Orchis purpurea</td><td></td><td></td><td>+</td></tr><tr><td>Arrhenatherum elatius</td><td></td><td>1</td><td></td><td>Plantago lanceolata</td><td></td><td></td><td>1</td></tr><tr><td>Bellis perennis</td><td></td><td></td><td>1</td><td>Poa pratensis</td><td></td><td></td><td>2</td></tr><tr><td>Bromus erectus</td><td></td><td></td><td>1</td><td>Poa trivialis</td><td>+</td><td>5</td><td></td></tr><tr><td>Carex divulsa</td><td></td><td></td><td>+</td><td>Potentilla reptans</td><td></td><td>1</td><td>1</td></tr><tr><td>Carex flacca</td><td></td><td></td><td>1</td><td>Prunus spinosa</td><td></td><td></td><td>1</td></tr><tr><td>Centaurea nigra</td><td></td><td></td><td>+</td><td>Rafanus spp.</td><td>+</td><td></td><td></td></tr><tr><td>Cirsium arvense</td><td></td><td>1</td><td></td><td>Ranunculus acris</td><td></td><td></td><td>1</td></tr><tr><td>Cirsium dissectum</td><td></td><td>1</td><td></td><td>Ranunculus bulbosus</td><td></td><td></td><td>1</td></tr><tr><td>Cirsium vulgare</td><td></td><td>+</td><td></td><td>Ranunculus repens</td><td></td><td>2</td><td></td></tr><tr><td>Convolvulus arvensis</td><td></td><td>1</td><td>1</td><td>Rasica oleacera</td><td>+</td><td></td><td></td></tr><tr><td>Crepis capillaris</td><td></td><td></td><td>1</td><td>Rhinanthus minor</td><td></td><td></td><td>+</td></tr><tr><td>Crepis spp.</td><td></td><td>+</td><td></td><td>Rubus spp.</td><td></td><td></td><td>+</td></tr><tr><td>Dactylis glomerata</td><td></td><td>1</td><td>3</td><td>Rumex acetosa</td><td></td><td></td><td>1</td></tr><tr><td>Daucus carota</td><td></td><td></td><td>1</td><td>Rumex crispus</td><td></td><td>1</td><td>+</td></tr><tr><td>Festuca arundinacea</td><td></td><td>2</td><td>3</td><td>Senecio jacobaea</td><td></td><td></td><td>1</td></tr><tr><td>Festuca rubra</td><td></td><td></td><td>1</td><td>Sonchus asper</td><td>+</td><td></td><td></td></tr><tr><td>Galium mollugo</td><td></td><td></td><td>1</td><td>Stachys officinalis</td><td></td><td></td><td>+</td></tr><tr><td>Gaudinia fragilis</td><td></td><td></td><td>1</td><td>Taraxacum officinalis</td><td></td><td>1</td><td>1</td></tr><tr><td>Holcus lanatus</td><td></td><td>1</td><td></td><td>Tragopogon pratensis</td><td></td><td>++</td><td></td></tr><tr><td>Hypericum perforatum</td><td></td><td></td><td>+</td><td>Trifolium dubium</td><td></td><td></td><td>2</td></tr><tr><td>Hypochaeris radicata</td><td></td><td>1</td><td>1</td><td>Trifolium pratense</td><td></td><td>1</td><td>2</td></tr><tr><td>Lathyrus pratensis</td><td></td><td></td><td>2</td><td>Trifolium repens</td><td></td><td></td><td>1</td></tr><tr><td>Leucanthemum vulgare</td><td></td><td></td><td>1</td><td>Verbena officinalis</td><td></td><td></td><td>1</td></tr><tr><td>Linum usitatissimum</td><td></td><td></td><td>1</td><td>Veronica arvensis</td><td></td><td></td><td>+</td></tr><tr><td>Lolium perenne</td><td></td><td>5</td><td></td><td>Vicia sativa</td><td>+</td><td></td><td>1</td></tr><tr><td>Lotus corniculatus</td><td></td><td></td><td>1</td><td></td><td></td><td></td><td></td></tr></tbody></table>
Acknowledgements

This work was partially supported by the French National Institute for Agricultural Research (INRA) and the French National Institute for Research in Computer Science and Automation (INRIA) Young Scientist Contract (CJS INRA-INRIA) and by the project SEBIOREF ("Promouvoir les Services Ecosystémiques rendus par la Biodiversité à l'agriculture", Région Occitanie, INRA, IRSTEA (French National Institute for Research in Sciences and Technologies for Environment and Agriculture)). The authors would like to thank the botanists who made possible the constitution of the field dataset used in this study: Philippe Caniot, Jérôme Willm, Emilie Andrieu and Gérard Balent. We would also like to thank the people accompanying the botanists: François Calatayud, Romain Carrié, Jean-Philippe Choisis, Donatien Dallery, Bruno Dumora, Camille Gouwy, Wilfried Heintz, Marc Lang and Anne-Sophie Mould. Special thanks to Clélia Sirami, Romain Carrié, Rémi Duflot and Nicolas Gross for sharing their knowledge about grasslands' species diversity. The authors would like to thank CNES for providing the pre-processed SPOT5 (Take5) data. We thank the anonymous reviewers for their valuable comments that greatly improved the manuscript.
References

[1] A. Eriksson, O. Eriksson, and H. Berglund, "Species abundance patterns of plants in Swedish semi-natural pastures," *Ecography*, vol. 18, no. 3, pp. 310-317, 1995.

[2] C. Critchley, M. Burke, and D. Stevens, "Conservation of lowland semi-natural grasslands in the UK: a review of botanical monitoring results from agri-environment schemes," *Biological Conservation*, vol. 115, no. 2, pp. 263–278, 2004.

[3] F. P. O'Mara, "The role of grasslands in food security and climate change," *Annals of Botany*, vol. 110, no. 6, p. 1263, 2012.

[4] B. P. Werling, T. L. Dickson, R. Isaacs, H. Gaines, C. Gratton, K. L. Gross, H. Liere, C. M. Malmstrom, T. D. Meehan, L. Ruan, B. A. Robertson, G. P. Robertson, T. M. Schmidt, A. C. Schrotenboer, T. K. Teal, J. K. Wilson, and D. A. Landis, "Perennial grasslands enhance biodiversity and multiple ecosystem services in bioenergy landscapes," *Proceedings of the National Academy of Sciences*, vol. 111, no. 4, pp. 1652-1657, 2014.

[5] A. Magurran, *Ecological Diversity and Its Measurement*. Croom Helm, 1988.

[6] C. E. Shannon, "A mathematical theory of communication," *Bell System Technical Journal*, vol. 27, no. 3, pp. 379-423, 1948.

[7] E. H. Simpson, "Measurement of diversity," *Nature*, 1949.

[8] D. Rocchini, N. Balkenhol, G. A. Carter, G. M. Foody, T. W. Gillespie, K. S. He, S. Kark, N. Levin, K. Lucas, M. Luoto, H. Nagendra, J. Oldeland, C. Ricotta, J. Southworth, and M. Neteler, "Remotely sensed spectral heterogeneity as a proxy of species diversity: Recent advances and open challenges," *Ecological Informatics*, vol. 5, no. 5, pp. 318–329, 2010. Special Issue on Advances of Ecological Remote Sensing Under Global Change.

[9] A. K. Skidmore, N. Pettorelli, N. C. Coops, G. N. Geller, M. Hansen, R. Lucas, C. A. Mücher, B. O'Connor, M. Paganini, H. M. Pereira, M. E. Schaepman, W. Turner, T. Wang, and M. Wegmann, "Environmental science: Agree on biodiversity metrics to track from space," *Nature*, vol. 523, p. 403, July 2015.

[10] D. Rocchini, D. S. Boyd, J.-B. Féret, G. M. Foody, K. S. He, A. Lausch, H. Nagendra, M. Wegmann, and N. Pettorelli, "Satellite remote sensing to monitor species diversity: potential and pitfalls," *Remote Sensing in Ecology and Conservation*, vol. 2, no. 1, pp. 25-36, 2016.

[11] N. Pettorelli, W. F. Laurance, T. G. O'Brien, M. Wegmann, H. Nagendra, and W. Turner, "Satellite remote sensing for applied ecologists: opportunities and challenges," *Journal of Applied Ecology*, vol. 51, no. 4, pp. 839-848, 2014.

[12] A. F. Cord, K. A. Brauman, R. Chaplin-Kramer, A. Huth, G. Ziv, and R. Seppelt, "Priorities to advance monitoring of ecosystem services using earth observation," *Trends in Ecology & Evolution*, vol. 32, no. 6, pp. 416-428, 2017.
[13] J. T. Kerr and M. Ostrovsky, "From space to species: ecological applications for remote sensing," *Trends in Ecology & Evolution*, vol. 18, no. 6, pp. 299-305, 2003.

[14] W. Gould, "Remote sensing of vegetation, plant species richness, and regional biodiversity hotspots," *Ecological Applications*, vol. 10, no. 6, pp. 1861-1870, 2000.

[15] M. W. Palmer, P. G. Earls, B. W. Hoagland, P. S. White, and T. Wohlgemuth, "Quantitative tools for perfecting species lists," *Environmetrics*, vol. 13, no. 2, pp. 121-137, 2002.

[16] J. Wilson, S. J. Fuller, and P. B. Mather, "Formation and maintenance of discrete wild rabbit (Oryctolagus cuniculus) population systems in arid Australia: Habitat heterogeneity and management implications," *Austral Ecology*, vol. 27, no. 2, pp. 183-191, 2002.

[17] J. Tews, U. Brose, V. Grimm, K. Tielbörger, M. C. Wichmann, M. Schwager, and F. Jeltsch, "Animal species diversity driven by habitat heterogeneity/diversity: the importance of keystone structures," *Journal of Biogeography*, vol. 31, no. 1, pp. 79-92, 2004.

[18] J. Oldeland, D. Wesuls, D. Rocchini, M. Schmidt, and N. Jürgens, "Does using species abundance data improve estimates of species diversity from remotely sensed spectral heterogeneity?," *Ecological Indicators*, vol. 10, no. 2, pp. 390-396, 2010.

[19] T. Möckel, J. Dalmayne, B. C. Schmid, H. C. Prentice, and K. Hall, "Airborne hyperspectral data predict fine-scale plant species diversity in grazed dry grasslands," *Remote Sensing*, vol. 8, no. 2, 2016.

[20] B. O. Oindo and A. K. Skidmore, "Interannual variability of NDVI and species richness in Kenya," *International Journal of Remote Sensing*, vol. 23, no. 2, pp. 285-298, 2002.

[21] D. H. K. Fairbanks and K. C. McGwire, "Patterns of floristic richness in vegetation communities of California: Regional scale analysis with multi-temporal NDVI," *Global Ecology and Biogeography*, vol. 13, no. 3, pp. 221-235, 2004.

[22] D. Rocchini, "Effects of spatial and spectral resolution in estimating ecosystem α-diversity by satellite imagery," *Remote Sensing of Environment*, vol. 111, no. 4, pp. 423-434, 2007.

[23] J.-B. Féret and G. P. Asner, "Mapping tropical forest canopy diversity using high-fidelity imaging spectroscopy," *Ecological Applications*, vol. 24, no. 6, pp. 1289-1296, 2014.

[24] M. Drusch, U. D. Bello, S. Carlier, O. Colin, V. Fernandez, F. Gascon, B. Hoersch, C. Isola, P. Laberinti, P. Martimort, A. Meygret, F. Spoto, O. Sy, F. Marchese, and P. Bargellini, "Sentinel-2: ESA's optical high-resolution mission for GMES operational services," *Remote Sensing of Environment*, vol. 120, pp. 25-36, 2012. The Sentinel Missions - New Opportunities for Science.

[25] D. U. Hooper, "The role of complementarity and competition in ecosystem responses to variation in plant diversity," *Ecology*, vol. 79, no. 2, pp. 704-719, 1998.

[26] S. Sakai, "Phenological diversity in tropical forests," *Population Ecology*, vol. 43, pp. 77-86, May 2001.
[27] J. Ryschawy, N. Choisis, J. P. Choisis, A. Joannon, and A. Gibon, "Mixed crop-livestock systems: an economic and environmental-friendly way of farming?," *Animal*, vol. 6, pp. 1722-1730, Oct 2012.

[28] R. Carrié, E. Andrieu, S. A. Cunningham, P. E. Lentini, M. Loreau, and A. Ouin, "Relationships among ecological traits of wild bee communities along gradients of habitat amount and fragmentation," *Ecography*, vol. 40, no. 1, pp. 85-97, 2017.

[29] O. Hagolle, M. Huc, D. Villa Pascual, and G. Dedieu, "A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENUS, LANDSAT and SENTINEL-2 images," *Remote Sensing of Environment*, vol. 114, pp. 1747-1755, Aug. 2010.

[30] P. H. C. Eilers, "A perfect smoother," *Analytical Chemistry*, vol. 75, no. 14, pp. 3631-3636, 2003. PMID: 14570219.

[31] C. Atzberger and P. H. C. Eilers, "Evaluating the effectiveness of smoothing algorithms in the absence of ground reference measurements," *International Journal of Remote Sensing*, vol. 32, no. 13, pp. 3689-3709, 2011.

[32] J. Braun-Blanquet, G. Fuller, and H. Conard, *Plant Sociology: The Study of Plant Communities: Authorized English Translation of Pflanzensoziologie*. McGraw-Hill publications in the Botanical Sciences, McGraw-Hill, 1932.

[33] L. Parsons, E. Haque, and H. Liu, "Subspace clustering for high dimensional data: A review," *SIGKDD Explor. Newsl.*, vol. 6, pp. 90-105, June 2004.

[34] C. Bouveyron, S. Girard, and C. Schmid, "High-dimensional data clustering," *Computational Statistics & Data Analysis*, vol. 52, no. 1, pp. 502-519, 2007.

[35] A. Lagrange, M. Fauvel, and M. Grizonnet, "Large-scale feature selection with gaussian mixture models for the classification of high dimensional remote sensing images," *IEEE Transactions on Computational Imaging*, vol. 3, pp. 230-242, June 2017.

[36] C. Bouveyron, S. Girard, and C. Schmid, "High-dimensional discriminant analysis," *Communications in Statistics - Theory and Methods*, vol. 36, no. 14, pp. 2607-2623, 2007.

[37] S. Girard and J. Saracco, "Supervised and Unsupervised Classification Using Mixture Models," in *EAS Publications Series*, vol. 77 of *EAS Publications Series*, pp. 69-90, May 2016.

[38] M. C. Neyrinck, I. Szapudi, and A. S. Szalay, "Rejuvenating the matter power spectrum: Restoring information with a logarithmic density mapping," *The Astrophysical Journal Letters*, vol. 698, no. 2, p. L90, 2009.

[39] K. Hall, T. Reitalu, M. T. Sykes, and H. C. Prentice, "Spectral heterogeneity of Quick-Bird satellite data is related to fine-scale plant species spatial turnover in semi-natural grasslands," *Applied Vegetation Science*, vol. 15, no. 1, pp. 145-157, 2012.

[40] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, "Scikit-learn: Machine learning in Python," *Journal of Machine Learning Research*, vol. 12, pp. 2825-2830, 2011.
[41] C. Biernacki, G. Celeux, and G. Govaert, "Assessing a mixture model for clustering with the integrated completed likelihood," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 22, pp. 719-725, jul 2000.

[42] J. C. Price, "How unique are spectral signatures?," *Remote Sensing of Environment*, vol. 49, no. 3, pp. 181-186, 1994.

[43] H. Nagendra, "Using remote sensing to assess biodiversity," *International Journal of Remote Sensing*, vol. 22, no. 12, pp. 2377-2400, 2001.

[44] D. Moog, P. Poschlod, S. Kahmen, and K.-F. Schreiber, "Comparison of species composition between different grassland management treatments after 25 years," *Applied Vegetation Science*, vol. 5, no. 1, pp. 99-106, 2002.

[45] E. Schäfer, J. Heiskanen, V. Heikinheimo, and P. Pellikka, "Mapping tree species diversity of a tropical montane forest by unsupervised clustering of airborne imaging spectroscopy data," *Ecological Indicators*, vol. 64, pp. 49-58, 2016.

[46] E. E. Maeda, J. Heiskanen, K. W. Thijs, and P. K. Pellikka, "Season-dependence of remote sensing indicators of tree species diversity," *Remote Sensing Letters*, vol. 5, no. 5, pp. 404-412, 2014.

[47] M. Pärtel and M. Zobel, "Small-scale plant species richness in calcareous grasslands determined by the species pool, community age and shoot density," *Ecography*, vol. 22, no. 2, pp. 153-159, 1999.

[48] H. H. Bruun, "Patterns of species richness in dry grassland patches in an agricultural landscape," *Ecography*, vol. 23, no. 6, pp. 641-650, 2000.

[49] D. Franzén and O. Eriksson, "Small-scale patterns of species richness in swedish semi-natural grasslands: the effects of community species pools," *Ecography*, vol. 24, no. 5, pp. 505-510, 2001.

[50] J. R. Bray and J. T. Curtis, "An ordination of the upland forest communities of Southern Wisconsin," *Ecological Monographs*, vol. 27, no. 4, pp. 325-349, 1957.

[51] S. L. Ustin and J. A. Gamon, "Remote sensing of plant functional types," *New Phytologist*, vol. 186, no. 4, pp. 795-816, 2010.

[52] L. Homolová, Z. Malenovský, J. G. Clevers, G. García-Santos, and M. E. Schaepman, "Review of optical-based remote sensing for plant trait mapping," *Ecological Complexity*, vol. 15, pp. 1-16, 2013.

[53] W. Jetz, J. Cavender-Bares, R. Pavlick, D. Schimel, F. W. Davis, G. P. Asner, R. Guralnick, J. Kattge, A. M. Latimer, P. Moorcroft, M. E. Schaepman, M. P. Schildhauer, F. D. Schneider, F. Schrodt, U. Stahl, and S. L. Ustin, "Monitoring plant functional diversity from space," *Nature plants*, vol. 2, no. 3, 2016.
[54] O. J. Abelleira Martínez, A. K. Premier, S. Günter, Z. Ramos Bendaña, L. Vierling, S. M. Galbraith, N. A. Bosque-Pérez, and J. C. Ordoñez, "Scaling up functional traits for ecosystem services with remote sensing: concepts and methods," *Ecology and Evolution*, vol. 6, no. 13, pp. 4359-4371, 2016.