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

HAL
open science

# Unsupervised Monitoring Vegetation after the Closure of an Ore Processing Site with Multi-Temporal Optical Remote Sensing

Sophie Fabre, Rollin Gimenez, Arnaud Elger, Thomas Rivière

► To cite this version:

Sophie Fabre, Rollin Gimenez, Arnaud Elger, Thomas Rivière. Unsupervised Monitoring Vegetation after the Closure of an Ore Processing Site with Multi-Temporal Optical Remote Sensing. Sensors, 2020, 20 (17), pp.4800. 10.3390/s20174800. hal-02946841

HAL Id: hal-02946841

https://hal.science/hal-02946841

Submitted on 8 Aug 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Unsupervised Monitoring Vegetation after the Closure of an Ore Processing Site with Multi-Temporal Optical Remote Sensing

Sophie Fabre ¹,*, Rollin Gimenez ¹, Arnaud Elger ² and Thomas Rivière ¹

¹ Onera, Département Optique et Techniques Associées, 2 Avenue Edouard Belin, 31055 Toulouse, France; rollin.gimenez@onera.fr (R.G.); thomas.riviere@onera.fr (T.R.)

² EcoLab, Université de Toulouse, CNRS, INPT, UPS, 31055 Toulouse, France; arnaud.elger@univ-tlse3.fr

* Correspondence: sophie.fabre@onera.fr

Received: 31 July 2020; Accepted: 19 August 2020; Published: 25 August 2020

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

check for
updates

**Abstract:** Ore processing is a source of soil heavy metal pollution. Vegetation traits (structural characteristics such as spatial cover and repartition; biochemical parameters—pigment and water contents, growth rate, phenological cycle...) and plant species identity are indirect and powerful indicators of residual contamination detection in soil. Multi-temporal multispectral satellite imagery, such as the Sentinel-2 time series, is an operational environment monitoring system widely used to access vegetation traits and ensure vegetation surveillance across large areas. For this purpose, methodology based on a multi-temporal fusion method at the feature level is applied to vegetation monitoring for several years from the closure and revegetation of an ore processing site. Features are defined by 26 spectral indices from the literature and seasonal and annual change detection maps are inferred. Three indices—CI<sub>red-edge</sub> (CIREDEDGE), IRECI (Inverted Red-Edge Chlorophyll Index) and PSRI (Plant Senescence Reflectance Index)—are particularly suitable for detecting changes spatially and temporally across the study area. The analysis is conducted separately for phyto-stabilized vegetation zones and natural vegetation zones. Global and specific changes are emphasized and explained by information provided by the site operator or meteorological conditions.

**Keywords:** vegetation survey; change detection; fusion; former mining site; trace metal elements; multispectral; multitemporal; Sentinel-2 satellite

## 1. Introduction

The state of vegetation, its survey and mapping are of growing interest to local territory stakeholders. Vegetation maps are an essential tool in the land use planning, creation and management of protected areas and in the monitoring of natural or polluted sites. Changes in vegetation provide information on the climatic, edaphic, geological and physiographic characteristics of an area [1]. In the natural environment, plants react to a multitude of biotic and abiotic factors.

The presence of pollutants in the growth medium (soil or water) can prove stressful for the vegetation and have direct physiological effects on plants (visual symptoms, photosynthetic process alteration, impact on development...). These effects result from the contact of the pollutants with the plants' roots and their assimilation in the tissues [2–4]. Soil pollutants can also induce indirect effects on vegetation by changing soil properties such as the water regime and nitrogen availability to the plant [5]. Vegetation is then considered as an indicator of soil pollution which can be detected by means of optical imaging, with the changes in the characteristics of vegetation (pigment content, cell structure, water content, physiological state...) altering its optical properties [6–12].

Vegetation mapping and monitoring nearby mine sites are of interest in all phases of mining, from the planning (assessing biodiversity in the area, locating mining infrastructure, having a reference
map beforehand to assess the environmental impact afterwards), mining and processing of ores (detecting potential contamination through their impact on vegetation) to closure and reclamation (providing information on the results of rehabilitation or detecting and controlling residual contaminants through vegetation health) [13,14]. In rural regions, mining activities are well documented as one of the most significant sources of soil heavy metal pollution [15–17].

Several works have reviewed the effects of trace metal elements, so-called heavy metals (HM), on vegetation reflectance [4,9]. The increase in HM in the soil has the effect of reducing the photosynthetic activity and chlorophyll concentration. Chlorophyll a and b, main pigments in leaves, present two light absorption peaks at 0.44–0.45 µm (blue) and 0.65–0.67 µm (red) [18]. The *red edge*, usually centered around 0.72 µm and identified by a fast-growing reflectance between 0.67 µm and 0.76µm, is employed through many approaches, based on spectral indices or radiative transfer models, to characterize leaf chlorophyll content variations [9,18,19]. Some experiments have also highlighted the interest of the Near InfraRed domain (NIR) plateau and the water absorption bands in the Short-Wave InfraRed domain (SWIR) to detect metal-stressed vegetation [20]. Vegetation sensitivity to HM depends on the contamination level, the species, the phase within the growth cycle and the environment [21].

Vegetation indices are simple or normalized reflectance ratios involving two or more spectral bands. Indices are applied on reflectance values or on transformations such as reflectance derivate or Continuum Removal (CR) [22]. They provide information on vegetation growth and stress due to their sensitivity to specific traits such as leaf chlorophyll or moisture contents, vegetation cover or biomass ... Spectral vegetation indices have both spatial and temporal dimensions and can be used to identify changes in vegetation phenology [23]. Reflectance changes can then be related to the contaminants present in the growing environment through vegetation indices [9]. The most widely known and used index in optical remote sensing is the Normalized Difference Vegetation Index (NDVI) [24]. NDVI is an overall standardized way in which to measure healthy vegetation in a wide range of applications. Many other indices have been proposed to compensate for the disadvantages of NDVI related to its sensitivity to soil brightness and its saturation when used for dense canopy cover [25]. Some of these indices, such as Photochemical Reflectance Index (PRI) or Ratio Vegetation Index (RVI), have been applied on hyperspectral or multispectral reflectance data to detect vegetation stress related to heavy metals [9,26]. Recent indices have also been proposed for HM-stress detection in rice crops such as Heavy Metal Stress Sensitive Index (HMSSI), which is the ratio between two existing red-edge indices $CI_{\text{red-edge}}$ (named CIREDEDGE) and PSRI (Plant Senescence Reflectance Index) or Heavy metal Cd stress-sensitive Spectral Index (HCSI) [27,28]. It is difficult to select one or more indices from the large number existing in the literature, as the index performance depends on the context defined by the species, the pollutant, the pollution level and the environment.

Clustering analysis is an unsupervised method without the need for a priori information that is applied to vegetation monitoring and classification [29–31]. Pixels are grouped according to their similarities and a label is assigned to each group. These similarities concern spatial patterns, spectral signatures, texture and even temporal information [32]. The clustering algorithms largely used in remote sensing can be classified in three main classes: centroid-based methods such as K-means algorithm, neighborhood-based methods and hierarchical clustering [33–36]. The K-means approach is the best known clustering algorithm widely used in the last few decades [29–31,37], even if it can be unstable for large data sets. The two other classes limit this drawback but have higher computational costs. In order to enhance the performance and reduce processing time without losing variability, dimension reduction techniques such as Principal Component Analysis (PCA) or spectral indices are often applied before clustering [30,35,37]. Moreover, the exploitation of the vegetation spectral indices rather than spectral signatures allows the clustering algorithm to focus only on specific vegetation traits [33].

A high temporal resolution can provide an important advantage in capturing the dynamic information on vegetation states [38]. Time series of medium spatial resolution optical data (e.g., decametric spatial resolution) demonstrated a high capacity—for example, in providing
inter-annual variations in the phenological stages from green-up to senescence [39], illustrating
vegetation growth [40], or quantifying the temporal and spatial trends of vegetation cover [6,41].

Satellite-derived vegetation indices with a high temporal coverage provide an opportunity with which to obtain the biophysical parameters of vegetation over large areas (regional or global). The Sentinel-2 (S-2) (A and B) multispectral instruments, with a spatial resolution of 10 and 20 m depending on the spectral domains and a high temporal coverage (five days revisit time), are one of the best ways in which to detect inter-annual and intra-annual vegetation changes as well as stressed vegetation [7,33,42]. The presence of the red-edge bands provides the opportunity for HM stress monitoring [27]. The exploitation of multi-temporal data should make it possible to reduce false alarms that may be linked to other sources of stress, such as weather-related stress (flood, drought...) or pests and diseases, also known as abrupt stressors [28]. Abrupt stressors can affect vegetation over one or several growth cycles or one whole year [43]. Soil HM induces stable stress with a consistent long-term influence on vegetation growth.

Recent works introduce clustering to process multitemporal remote sensing and analysis land
cover evolution, essentially for agricultural monitoring [31,37] and land use change [29,30]. Such time
series data mining techniques allow us, by trend analysis, to examine shocks and unexpected variations
or to detect anomalies and stable similar patterns. In agriculture, this information permits us to
distinguish dissimilar levels of vegetation cover towards parcels as well as identify the planting pattern
and phenology anomalies across development stages [29,31]. Times series are exploited by selecting
sets of dates appropriate to the information to highlight [31,37]. The exploitation of multiple dates
provides then the evolution between the selected periods.

Change detection is the process that analyses multi-temporal remote sensing images acquired
over the same area in order to identify the changes occurring between the considered acquisitions
dates [44]. Many publications show a comprehensive review of the change detection techniques [45-48].
These techniques are fusion methods in the time domain and can be classified according to the
fusion level, whether feature-level or decision-level fusion [49]. Feature-level fusion methods are
unsupervised automatic detection methods avoiding the need to perform in situ measurements and
based on change indicators chosen to enhance multi-temporal information. Decision-level fusion
methods are essentially made of supervised or semi/partially supervised/unsupervised classification
methods. They include post-classification comparison, supervised direct multi-date classification
and compound classification [46]. Supervised classification approaches improve performance in
comparison to the unsupervised ones. They are, however, less usable in an operational context and
depend on the learning base quality (number of samples, distribution of samples per class...).

Selecting an appropriate change detection technique depends on many factors: the study objective and the context, the size of the study area, the multi-temporal image spatial resolution, the in situ information availability, etc. In this work, change detection is applied in order to automatically identify the long-term changes related to behavior modifications between time series, using the temporal signature of land cover, and also to monitor seasonal/annual changes without a priori information about the change (localization, impact date, etc.). The change detection method based on feature-level fusion is suited to our needs. Several mathematical operators, relying mainly on the difference operator, are provided in the literature to define suitable change indicators [46]. Spectral indices have been used as detection indicators in order to highlight changes in many applications (land or forest cover change, burned area extraction...) [30,50-52]. The fusion of multi-temporal data can then be performed by image comparison [46]. In order to provide simple and fast detection maps, different operations (such as ratio, difference) can be applied to the original images or to the detection index images [45,47].

The main purpose of this study is to propose an unsupervised methodology applied on multispectral satellite time series to monitor vegetation for several years from the closure and revegetation of an ore processing site. This methodology implements, on the one hand, a simple and fast multi-temporal fusion method at the feature level (features inferred from spectral index selection) to detect simultaneously abrupt and longer-term changes in the vegetation development and traits
(cover density, biochemical characteristics and phenology stages). On the other hand, a clustering
recognition of spatial patterns within areas of the study site is provided for retrieving temporally
stable (constant) patterns over the season and the year. The S-2 time series presents the geometric
details and the spectral information required for jointly operating spectral, spatial and temporal
information. The high temporal repeatability ensures both an analysis of the inter-annual variability of
habitat changes over a large temporal period and an analysis of the habitat phenology over one year.
Some areas of the study site have been phyto-stabilized, whilst in other areas natural (spontaneous)
vegetation has developed. Phyto-stabilization has been enhanced by using soil amendments that
immobilize HM combined with plant species that are tolerant of high levels of contaminants and
low-fertility soils [13]. This phyto-management process allows us to partially confine HM by limiting
its dissemination in the environment (through wind or water erosion, leaching).

The work is divided into several stages:

* Creation of the Sentinel-2 time series covering each season for many years;

* Construction of the spectral index base composed of numerous indices specified or used in various contexts (environment monitoring, biodiversity assessment, abiotic impact analysis, precision agriculture...);

* Development of an unsupervised diagnostic methodology based on the multi-temporal fusion methods for change detection and clustering recognition for pattern identification;

* Analysis of the intra-annual and inter-annual variabilities of the most pertinent spectral indices for the studied context in order to define the most appropriate phenological stage and to detect changes in vegetation over several years (temporal evolution comparison of phyto-stabilized vs. natural vegetation areas).

## 2. Materials and Methods

### 2.1. *Study Area*

The study area is a former ore processing site in France that was closed in 2004 and is managed
by a French public institute. In 2006, some areas of the site were phyto-stabilized. The site covers
approximately 120 hectares. Six zones of interest defined by the institute were analyzed in this study
(Figure 1). Z1, Z2 and Z5 are the three aided phyto-stabilized zones. The phyto-stabilization is
realized at the same period over most of the areas. Some differences in the management mode concern
zones Z1 and Z5. On the central part of Z1, natural uncontaminated soil was introduced. In Z5,
only the northeastern half was phyto-stabilized by a steel shot mixture. Spontaneous vegetation
developed on the other part of the site and, in addition, trees were planted in alignment. According
to the maps provided by the institute, Z2 and Z5 have the highest metal signatures in their soil.
In situ measurements cannot be provided because of the site sensitivity and for intellectual property
reasons, but a map with two metal signature classes is provided in the result analysis. Z3, Z4 and Z6,
identified as natural vegetation areas, have been vegetated with a mixture of local seeds on clean soil
brought and the natural vegetation has developed over the years. Species are known for zones 2, 3, 5
and 6 and are provided afterwards.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Former ore processing site and its six zones of interest, extracted from BD ORTHO® 50 cm, ortho-photography provided by IGN (Institut Géographique National) with a spatial resolution of 50 cm—scale: 1:10,000. Area of each zone (in ha): Z1 (6.5), Z2 (2.8), Z3 (9.4), Z4 (8.5), Z5 (19.4), Z6 (35.3).</figcaption></figure>

In order to provide information on the relief of the study site, the BD-ALTI® 25 m, Digital Terrain Model (DTM) provided by Institut Géographique National (IGN) and describing the relief of the French territory at a spatial resolution of 25 m, is processed by Geospatial Data Abstraction Library (GDAL) in order to extract the contour lines at 10 m elevation intervals [44]. The elevation in the study area ranges from approximately 180 to 280 m. Z6 has the greatest difference in height (about 75 m), and Z1 and Z4 present the lowest differences in height (about 15 m).

## 2.2. Sentinel-2 Satellite Imagery

The Sentinel-2A images covering the (0.4–2.5 µm) domain available on the study site are selected according to their quality (clear sky, processing level 2A surface reflectance products after geometric and atmospheric corrections) between May 2016 and May 2020—i.e., starting and ending in spring with, in theory, maximal green cover [53,54].

Some other considerations are taken into account when choosing the processed dates. For each year, the seasons must be the same in order to detect abrupt changes. Seasons must be properly sampled over the year to monitor the phenology evolutions. The choice is then made to select one image per season. The time series consists of 17 images listed in Table 1. No image is selected in February 2016 due to the presence of clouds.

<table><caption>**Table 1.** Time series of Sentinel-2 images.</caption><thead><tr><td>Image id Number</td><td>Acquisition Date <sup>1</sup></td></tr></thead><tbody><tr><td>SENTINEL2A_20160521-105553</td><td>2016/05/21</td></tr><tr><td>SENTINEL2A_20160826-104023</td><td>2016/08/26</td></tr><tr><td>SENTINEL2A_20161015-104513</td><td>2016/10/15</td></tr><tr><td>SENTINEL2A_20170225-105020</td><td>2017/02/25</td></tr><tr><td>SENTINEL2A_20170516-105322</td><td>2017/05/16</td></tr><tr><td>SENTINEL2A_20170821-104208</td><td>2017/08/21</td></tr><tr><td>SENTINEL2A_20171013-105315-</td><td>2017/10/13</td></tr><tr><td>SENTINEL2A_20180227-104236</td><td>2018/02/27</td></tr><tr><td>SENTINEL2A_20180511-105804</td><td>2018/05/11</td></tr><tr><td>SENTINEL2A_20180819-105124</td><td>2018/08/19</td></tr><tr><td>SENTINEL2A_20181025-104115</td><td>2018/10/25</td></tr><tr><td>SENTINEL2A_20190222-104853</td><td>2019/02/22</td></tr><tr><td>SENTINEL2A_20190513-104901</td><td>2019/05/13</td></tr><tr><td>SENTINEL2A_20190814-105857</td><td>2019/08/14</td></tr><tr><td>SENTINEL2A_20191030-104901</td><td>2019/10/30</td></tr><tr><td>SENTINEL2A_20200220-105848</td><td>2020/02/20</td></tr><tr><td>SENTINEL2A_20200520-105859</td><td>2020/05/20</td></tr></tbody></table>

<footnote n="1">¹ Date format: year/month/day.</footnote>
Ref. [7,55] provide the central wavelength positions and spatial resolutions for the Sentinel-2A spectral bands. A multispectral data cube (x, y positions and wavelengths) with a spatial resolution of 10 m is produced for each image referenced in Table 1. For this purpose, the spectral bands at 20 m spatial resolution are oversampled to a 10 m spatial resolution using a nearest neighborhood filter in order to preserve the spectral information [56]. For each image, a thumbnail image associated with the study site is extracted (Figure 2).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** RGB (Red-Green-Blue) representation of a thumbnail image extracted from a processed Sentinel-2 image (10 m spatial resolution).</figcaption></figure>

## 2.3. Methodology

As mentioned above, the most suitable method for our application case is a multi-temporal fusion method at the feature level. The fusion is processed in two stages:

1. Feature extraction: various spectral indices are used to highlight features;

2. Feature comparison: a mathematical operator is applied to highlight the changes that occurred in bi-temporal image pairs.

The processing chain includes four modules (Figure 3): a mask generation module to identify the vegetation pixels to be processed thereafter, a spectral indices database to produce feature maps, a spectral and temporal clustering module and an image comparison module to detect change. The inputs are spectral reflectance images and the outputs are change detection maps.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Methodology description.</figcaption></figure>

### 2.3.1. Mask Generation

For each S-2 image, a mask is generated so as to identify the vegetation pixels to be processed and to mask the pixels associated with other land use classes (bare soil, water, buildings...). The following combination of spectral indices is applied on the reflectance image to produce its corresponding mask:
* Normalized Built-up Area Index (NBAI) and Band Ratio for Built-up Area (BRBA) to identify road or building pixels [57,58];

* Normalized Difference Water Index (NDWI) and modified Normalized Difference Water Index (mNDWI) to highlight water areas [59];

* Soil Adjusted Vegetation Index (SAVI) and Dry Bareness Spectral Index (DBSI) to identify bare soil pixels [60,61];

* NDVI to provide vegetation cover pixels [24].

The combination of NDVI with indices related to other land use classes allows a reduction in the false alarm rate. A binary thresholding is applied to each index map, and a histogram is used to determine the appropriate threshold value. The index maps after thresholding are combined by an AND logical operator.

Figure 4 shows the mask obtained from the multispectral image acquired on 21/05/2016. The black pixels represent non-vegetated areas and the vegetation pixels are white. The vegetation masks are analyzed and compared in the next to monitor the development and evolution of the spatial distribution of vegetation between 2016 and 2020.

<figure><img src="image_8.png" /><figcaption>Figure 4. Vegetation mask from the multispectral image acquired on 21 May 2016.</figcaption></figure>

### 2.3.2. Spectral Index Database

Table 2 provides the description of the spectral indices database. Each index formulation is provided and the corresponding reference is cited. The S-2 spectral bands used by each index are provided ([7,55] defining the central wavelength and the spatial resolution for each band).

<table><caption>Table 2. Spectral indices database.</caption><thead><tr><th colspan="2">S2 Spectral Band</th><th>B1</th><th>B2</th><th>B3</th><th>B4</th><th>B5</th><th>B6</th><th>B7</th><th>B8</th><th>B11</th><th>B12</th></tr><tr><th colspan="2">Spectral Band Formulation <sup>1</sup></th><th>Ae.</th><th>B</th><th>G</th><th>R</th><th>RE1</th><th>RE2</th><th>RE3</th><th>NIR</th><th>SWIR</th><th>SWIR</th></tr><tr><th>Spectral Index</th><th>Description</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>ARVI</td><td>Atmospherically Resistant Vegetation Index [62]</td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C_1D</td><td>Maximum of first derivative reflectance at the red edge [63].</td><td></td><td></td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td><td></td></tr><tr><td>CIREDEDGE</td><td>Red-edge Chlorophyll Index [64]</td><td></td><td></td><td></td><td></td><td>x</td><td></td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>CTR</td><td>Carter index [65]</td><td></td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>EMEN2</td><td>Red and green ratio [66]</td><td></td><td></td><td>x</td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FIRSTDV_CHL</td><td>Ratio of the derivative at red-edge [67]</td><td></td><td></td><td></td><td></td><td>x</td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GREEN_NDVI</td><td>Green NDVI [68]</td><td></td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td></tr></tbody></table>
<table><caption>Table 2. *Cont.*</caption><thead><tr><th colspan="2">S2 Spectral Band</th><th>B1</th><th>B2</th><th>B3</th><th>B4</th><th>B5</th><th>B6</th><th>B7</th><th>B8</th><th>B11</th><th>B12</th></tr><tr><th colspan="2">Spectral Band Formulation <sup>1</sup></th><th>Ae.</th><th>B</th><th>G</th><th>R</th><th>RE1</th><th>RE2</th><th>RE3</th><th>NIR</th><th>SWIR</th><th>SWIR</th></tr><tr><th>Spectral Index</th><th>Description</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>HMSSI</td><td>Heavy Metal Stress Sensitive Index [27]</td><td></td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td></tr><tr><td>IRECI</td><td>Inverted Index Red-Edge Chlorophyll Index [7]</td><td></td><td></td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td></tr><tr><td>mARI</td><td>Modified Anthocyanin Reflectance Index [69]</td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>MCARI</td><td>Modified Chlorophyll Absorption Ratio Index [70]</td><td></td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>MCARI_NEW</td><td>New MCARI [71]</td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>MCARI_OSAVI</td><td>MCARI/Optimized SAVI [72]</td><td></td><td></td><td>x</td><td>x</td><td>x</td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>MCARI_OSAVI_NEW</td><td>[72]</td><td></td><td></td><td>x</td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td></tr><tr><td>MSAVI</td><td>Modified Soil Adjusted Vegetation Index [73]</td><td></td><td></td><td></td><td>x</td><td></td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>MTCI</td><td>MERIS Terrestrial Chlorophyll Index [74]</td><td></td><td></td><td></td><td></td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td></tr><tr><td>NBR</td><td>Normalized Burn Ratio [75]</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td><td>x</td></tr><tr><td>NDRE</td><td>Normalized Difference Red Edge [76]</td><td></td><td></td><td></td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td></td></tr><tr><td>NDVI</td><td>Normalized Difference Vegetation Index [77]</td><td></td><td></td><td></td><td>x</td><td></td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>PSRI</td><td>Plant Senescence Reflectance Index [78]</td><td></td><td>x</td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td><td></td></tr><tr><td>S2REP</td><td>Sentinel-2 red-edge position [7]</td><td></td><td></td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td></tr><tr><td>SAVI</td><td>Soil Adjusted Vegetation Index [60]</td><td></td><td></td><td></td><td>x</td><td></td><td></td><td></td><td>x</td><td></td><td></td></tr><tr><td>SIPI</td><td>Structure Insensitive Pigment Index [79]</td><td>x</td><td></td><td></td><td>x</td><td></td><td></td><td></td><td>x</td><td></td><td></td></tr><tr><td>TCARI_OSAVI</td><td>[72,79]</td><td></td><td></td><td>x</td><td>x</td><td>x</td><td></td><td>x</td><td></td><td></td><td></td></tr><tr><td>TCARI_OSAVI_NEW</td><td>[72]</td><td></td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td><td></td></tr><tr><td>VII</td><td>Vegetation Inferiority Index [80]</td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td></td><td></td></tr></tbody></table>

A total of 25 indices use spectral bands in the Visible-Near InfraRed (VNIR) spectral domain, and only one index, NBR, exploits the band B12 in the Short Wave InfraRed (SWIR) domain. There are many indices exploiting the red edge. This database includes:

* Indices specified for global vegetation monitoring such as NDVI [77] or MTCI (MERIS Terrestrial Chlorophyll Index) [74];

* Indices defined to characterize the impact of contaminants (heavy metals, hydrocarbons...) on vegetation such as EMEN2 (Red and green ratio) [66] or HSSMI (Heavy Metal Stress Sensitive Index) [27];

* Indices defined to estimate biochemical parameters (chlorophyll, anthocyanin...) such as ARVI [62] or SIPI [79].

* NDVI, widely used in many applications [24], is chosen as the reference index and the other indices are compared to this reference index. The best-performing indices are selected and the complementarity of the selected indices are analyzed.

<footnote n="1">¹ Spectral band formulation: Ae (Coastal aerosol), B (Blue), G (Green), R (Red), RE (Vegetation Red Edge), NIR (Near-InfraRed), SWIR (Short Wave InfraRed).</footnote>
### 2.3.3. Spectral and Temporal Clustering

#### Description of the Clustering Methods

Clustering is used to regroup and classify pixels with the same spectral and temporal behaviors [29–31,80]. After applying a vegetation spectral index on the time series, for each pixel a matrix regrouping spectral and temporal information is constructed:

$$M_{i,j} = [I_{i,j,t_1}, \dots, I_{i,j,t_n}] \quad (1)$$

where *I* is the considered spectral index or vector of several spectral indices, *i*, *j* are the coordinates of the pixel of interest and *t*<sub>*k*</sub> represents the date *k* of the time series composed of *n* dates *k* ∈ [1, *n*].

Several clustering methods used in the literature are compared in order to select the most adapted for our study (methods are provided by the scikit-learn library [81]). These methods are partitioned in three classes:

*   **Centroid-based methods, like the K-means algorithm [33,34]:** This consists of randomly initializing K centroids to assign each data point (in this case $M_{i,j}$) to its nearest centroid according to the Euclidean distance. Then, the centroids are recomputed as the mean (in the case of the K-means algorithm) of the data point of each group. These steps are finally repeated until attempting a stopping criterion. For its efficiency, K-means is selected and its instability is solved by running 100 times (value empirically fixed) the algorithm with different centroid seeds and selecting the best output in terms of inertia [33]. Still in this class, another chosen method consists of adding a principal component analysis (PCA) algorithm to reduce the dimension of the time series without losing its variance before applying K-means;

*   **Neighborhood-based methods, such as spectral clustering [35,36]:** The dataset is processed by creating a similarity graph by k-nearest neighbors (graph constructed by connecting each point to its K-nearest neighbors) or ε-neighborhood (graph obtained by associating to each point all points falling inside a circle of radius ε, where ε is a real value). Then, a corresponding Laplacian matrix is created to project the data onto a lower dimensional space, and K-means is applied on its rows. Both the ε-neighborhood and k-nearest neighbor methods are tested. The high computational cost obligates this class of method to provide for some adjustments such as image preprocessing [37,82]. In our case, the preprocessing consists of applying a spectral index.

*   **Hierarchical clustering algorithms:** partitions are obtained by the successive groupings of objects (bottom-up approach such as agglomerative clustering, the algorithm chosen in our study) or by successive division (top-down approach) [33,34].

#### Performance Comparison of the Clustering Methods

According to former works on clustering, evaluating the accuracy of a clustering classification is far from easy because clustering results do not necessary represent in situ measurable information [37]. Thus, to select an algorithm and fix its input parameters, three criteria defined according to published studies are considered [33,37]:

*   The computational cost of each algorithm.

*   The internal quality of clusters, defined by analyzing the compactness of each cluster and the separation between clusters. This information is evaluated by the silhouette score (or index), which indicates if a pixel belongs to the right cluster [83]. A score close to 1 indicates a good clustering. The average of these scores over all the data indicates if the method has appropriately clustered the data.

*   The external quality of clusters, defined by comparing the clusters to externally supplied information. To do this, the Moran's Index (Moran's I) is applied to spectral index maps. It is a measure of spatial autocorrelation which can be calculated over all the data (global Moran's I) or
in a specific zone (local Moran's I or LISA Local Indicators of Spatial Association) [84]. While the global index gives the spatial correlation of the entire dataset by returning a number between -1 and 1, the local index allows drawing a correlation map. Both are used to measure the efficiency of our clustering. The LISA map returns a map described by five classes (Figure 5). Two, usually called High-High or Low-Low, indicate regions which contribute significantly to a positive global spatial autocorrelation outcome; two others, called High-Low and Low-High, indicate those which contribute significantly to a negative autocorrelation outcome and one that indicate those with an insignificant contribution. Thus, the comparison of the classes obtained by the Moran index with the clusters allows us to verify the relevance of our clusters.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparison of the K-means clustering map with the LISA (Local indicators of spatial association) and NDVI (Normalized Difference Vegetation Index) maps. On LISA map: red clusters are pixels of high NDVI values (in comparison with the average NDVI value) surrounded by pixels of high NDVI values, orange cluster are pixels of high values surrounded by pixels of low values, pale blue those of low values surrounded by pixels of high values and blue those of low values in a group of pixel of blue values.</figcaption></figure>

This analysis is performed for a single date or combined dates and spectral indices in order to select the clustering algorithm among the following methods: K-means, PCA and K-means, spectral clustering with k-nearest neighbors, spectral clustering with ε-neighborhood, hierarchical clustering.

An example of the K-means performance evaluation is illustrated by Figure 5 for a specific date (21 May 2016). The silhouette score is equal to 0.53, which is a significant score. The global Moran's I is equal to 0.87, which indicates strongly grouped and correlated NDVI values across the map. The LISA map shows that the clusters returned by the LISA method are almost the same as those returned by the K-means algorithm.

## Clustering Method Selection

The simple K-means and principal component analysis followed by K-means are slightly quicker than the other methods. The silhouette scores are globally coherent according to the algorithms, but simple K-means presents barely better results than other algorithms. Finally, by comparing the maps returned by the different algorithms, all are consistent with Moran's I. Small differences observed only are the boundaries between clusters. Since the simple K-means algorithm is slightly better in terms of computational time and silhouette score, this algorithm is selected.

## Input Parameter

The cluster number, the input parameter of the K-means algorithm, has to be adapted to the application case. The silhouette score is initially used to fix the optimal cluster number. Two, three, four and five clusters are tested. The difference values of the silhouette score are not sufficient to decide the number of clusters on this criterion.

The cluster number selected has to represent the fluctuation range of the index values according to the season and the studied zones. Two or three clusters are not sufficient to represent the seasonal and inter-zones variabilities. Five clusters do not give more information than four clusters. Moreover, according to the literature [85], the NDVI value range is often divided into four classes: lower than 0.2 for bare soil, approximately 0.2 to 0.4 for sparse vegetation (like shrubs, grassland, or senescing crops),
approximately between 0.4 and 0.6 for moderate vegetation, and approximately higher than 0.6 for a high density of green vegetation, such as that found in temperate and tropical forests or crops at their peak growth stage. This leads to us retaining four clusters.

### 2.3.4. Image Comparison for Change Identification

Simple change detection techniques providing easy to interpret results are chosen [45]. A hypothesis of precise geometric registration between multi-temporal images to avoid false change areas being detected due to an image displacement is retained [47]. This assumption is realistic, as the multi-temporal images are acquired by a single instrument. The suggested change detection methods are based on the absolute, relative and normalized differences between two spectral index maps at two different dates. The difference is applied pixel by pixel. These methods are compared, and finally only the normalized difference is studied as normalization allows one to compare the difference values from one date to another. It is based on the following formulation:

$$ dif\_norm_{d1,d2} = \frac{I_{d1} - I_{d2}}{I_{d1} + I_{d2}} \qquad (2) $$

where $d_1$ and $d_2$ are two different dates, and $I_{di}$ is the map associated with a given spectral index (Section 2.3.2) and the date $d_i$ ($i = 1, 2$). $I_{di}$ can represent the spectral and temporal clustering map (Section 2.3.3) associated with a given time series $d_i$.

The vegetation masks associated with each date are combined by an OR operator in order to provide the vegetation mask of the detection map.

Figure 6 illustrates an example of change detection map obtained by applying the normalized difference to two NDVI maps. The decrease in the NDVI values in Z4 is quickly identified on the detection map.

<figure><img src="image_10.png" /><figcaption>Figure 6. Example of detection map displaying the evolution of NDVI between spring 2016 and spring 2019 (reference year: 2016).</figcaption></figure>

### 2.3.5. Data Analysis

A large volume of data has to be compared. A data analysis strategy is then put in place. Two types of result analysis are carried out:

* A global visual analysis and a visual analysis per zone of the various generated maps. The maps are vegetation masks (Section 2.3.1), spectral index maps (Section 2.3.2), change detection maps applied on spectral indices (Section 2.3.4), clustering maps (Section 2.3.3) and change detection applied on clustering maps (Section 2.3.4).

* A statistical analysis per zone. On one hand, the following statistics are calculated from the generated maps: mean, standard deviation (noticed STD), coefficient of variation (also known as RSD, Relative Standard Deviation) [86], median, first and third quartiles (Q1 and Q3, respectively).
On the other hand, the temporal evolution of the median and mean are analyzed and median
results are provided.

The statistics and maps obtained for the indices described in Table 2 are compared. At first, the indices characterized by temporal and spatial variabilities are chosen by a visual analysis of the index maps. Among these indices, the ones with similar temporal changes over seasons are grouped together by a statistical analysis of the index temporal evolution. For each group, the indices are compared in order to eliminate the indices leading to redundant information. Finally, the change detection maps and associated statistics obtained for the selected indices are then compared in order to define the best-suited indices for change detection.

Spectral and temporal clustering is applied to the selected indices (unique index or combination of indices) and a temporal series in order to highlight stable spatial patterns, monitor the annual evolution, and emphasize phenological dissimilarities over the year. Change detection is then applied to the clustering maps.

**3. Results**

Because of the significant amount of results obtained, only the fundamental outcomes are provided within this section. Results are completed by the Appendix A.

3.1. *Temporal Evolution of the Vegetation Cover*

Vegetation masks are used to monitor the temporal changes in the vegetation cover density and spatial distribution in each area.

Figure 7 represents the seasonal evolution of vegetation for the year 2017. In summer and autumn, most parts of the monitored areas have no vegetation. The two seasons with the most significant vegetation cover are winter and spring. This observation is valid for each year. For this studied site, these seasons are those to favor for vegetation monitoring.

<figure><img src="image_11.png" /><figcaption>Figure 7. Seasonal evolution of the vegetation mask for the year 2017.</figcaption></figure>

Figure 8 shows the temporal evolution of the percentage of vegetation pixels for each area. The same vegetation cycle is observed for phyto-stabilized and natural vegetation zones. Globally, this cycle corresponds to the reduction in vegetation cover in summer and growth in spring. Globally, there is an increase in vegetation cover over the years.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Temporal evolution of the percentage of vegetation pixels for each zone from 2016 to 2020, for phyto-stabilized (a) and natural vegetation (b) areas.</figcaption></figure>

Larger deviations are identified between summer and spring for two zones with natural vegetation (zones 3 and 4). The phyto-stabilized zones and zone 6 with natural vegetation maintain a greater vegetation cover in summer and autumn than zones 3 and 4. Zone 6 contains well-developed trees, mainly pines and poplars, and Spanish brooms (*Spartiumjunceum*). This explains the vegetation cover higher than 30% in summer and autumn. The difference in vegetation covers between zones 3 and 6 is explained by the presence of distinct species. The five dominant species in zone 3 are species present in the Mediterranean region: *Aphyllanthesmonspeliensis* (plant abundantly flowering in spring, forming rush-like clumps with blue flowers), *Bituminiariabituminosa* (a perennial plant that grows on the edges of fields or paths), *Dittrichiaviscosa* (an invasive perennial which flowers in late summer and early autumn, often growing in wasteland and rubble and along roadsides, forming abundant green tufts with yellow flower heads), *Pallenisspinosa* (an annual or biennial plant generally not more than 70 cm high with dark, stiff, straight stems with spaced hairs) and *Plantagolanceolata* (commonly known as Ribwort plantain, a perennial herbaceous plant of average size 15–50 cm, which exhibits variable forms depending on soil nutrient richness, sunlight exposure and soil moisture). Zone 3 also includes Spanish brooms and some developing trees (pines, poplars) on the edge of the zone.

The percentages of vegetation pixels in zone 1 are close to those in zone 5. Zones 1 and 5 have more vegetation pixels in summer-autumn than zone 2 and, conversely, fewer vegetation pixels in spring compared to zone 2 (differences in the order of 10% to 20% depending on the year considered until autumn 2019). After autumn 2019, the vegetation cover of zone 2 exceeds that of the other two areas. Zone 1 species have not been provided. Zones 2 and 5 are composed of the same species. The only dominant species common with those in zone 3 is *Dittrichiaviscosa*. The other three dominant species identified in zones 2 and 5 are *Blackstoniaperfoliata* (annual herbaceous plant that flowers from May to September with yellow flowers), *Securigeravaria* (perennial herbaceous plant that flowers from June to September), *Trisetumflavsecens* (perennial herbaceous plant of the *Poaceae* family, growing in thick clumps, up to 60 to 80 cm high). The presence of herbaceous vegetation in these areas explains the high pixel rates of vegetation cover in summer and autumn in comparison to the rates of zone 3.

After autumn 2018, an abrupt change in the vegetation cycle is observed. To find out the reason for this phenology alteration, the available meteorological data are analyzed. The change in vegetation cycle after autumn 2018 is explained by the heavy precipitations in October 2018, almost three times higher than the precipitations in October 2016 (Figure 9). Vegetation cover appears to increase overall across the years (Table 3). This increase is mostly highlighted after the autumn 2018 and can be correlated with the significant precipitation levels of autumn 2019 and spring 2020. A change in the phenological stages is observed and vegetation cover is increasing in autumn from 2018 (Table 3). Prior to this date, the vegetation cover in autumn was comparable to the one during summer. In spring,
the vegetation cover is the most important and remains relatively stable over the years (maximum increase of 5%). There are low apparent differences between the phyto-stabilized and natural vegetation areas, except for the highest rates of increase in autumn and the lowest increases in winter and spring over these 4 years for the natural vegetation areas.

<figure><img src="image_13.png" /><figcaption>**Figure 9.** (a) Temporal evolution of the percentage of vegetation pixels by zone (between May 2016 and May 2020); (b) seasonal evolution of the monthly rainfall during 5 consecutive years (in mm).</figcaption></figure>

<table><caption>**Table 3.** Rate of vegetation cover increase by season between 2016 (or 2017 in winter) and 2020. The summer results are not provided owing to the low vegetation cover.</caption><thead><tr><th rowspan="2">Zone</th><th colspan="3">Increase in Vegetation<br>Cover Pixels (%)</th></tr><tr><th>Spring</th><th>Autumn</th><th>Winter</th></tr></thead><tbody><tr><td>Z1 *</td><td>5</td><td>30</td><td>14</td></tr><tr><td>Z2 *</td><td>1</td><td>23</td><td>16</td></tr><tr><td>Z3</td><td>1</td><td>33</td><td>12</td></tr><tr><td>Z4</td><td>0</td><td>31</td><td>2</td></tr><tr><td>Z5 *</td><td>6</td><td>26</td><td>13</td></tr><tr><td>Z6</td><td>0</td><td>43</td><td>3</td></tr></tbody></table>

The evolution of the vegetation cover is analyzed as a function of the average monthly temperature. February 2018 is the coldest month within the analyzed period, and the increase in the number of vegetation pixels between autumn and winter 2018 is smaller than the increase observed in 2017 (Figure 9).

An analysis of the vegetation masks reveals areas of permanent bare soil with no vegetation cover regardless of the year and the season. These areas are identified on Figure 10a and are related to the HM field measurements in zones 2, 3 and 5 (Figure 10b). The bare soil of zone 5 corresponds to marked metal signatures in soil. The areas with the highest HM contents are not necessarily identified as areas of permanent bare soil in the image time series. Some vegetation has the capacity to grow on HM-contaminated soil.

<footnote>(*) Phyto-stabilized areas.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** (a) Identification of areas without vegetation cover; (b) superposition of the metal signatures classes provided by in situ measurements.</figcaption></figure>

## 3.2. Spectral Index Selection

The selection and evaluation of the spectral indices is carried out in several steps. Firstly, the NDVI reference index maps are visually analyzed. Figure 10 shows, for example, the seasonal evolution of the NDVI maps for the year 2017. In addition, the NDVI statistics (mean, standard deviation, RSD) are computed by the area and vegetation type (i.e., vegetation introduced for phyto-stabilisation, natural vegetation) (Section 2.3.4). Figure 11 provides the annual temporal evolution of NDVI statistics by zone.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Seasonal evolution of NDVI reference index maps for the year 2017 (masked pixels corresponding to non-vegetation pixels, area shapes representing the contours of the study zones).</figcaption></figure>

The NDVI analysis leads to the following common observations throughout the study area:

* The NDVI index presents the same phenology stages for all zones except zone 5. NDVI increases in spring and decreases in summer.

* A change is observed in winter 2019; the increase in the NDVI value is less marked than those of the other years (a decrease is even observed for zones 4 and 5). This can be explained by the heavy rainfalls in autumn 2018. This abrupt change is also identified during the analysis of the percentage of vegetation pixels (Figure 11).
* Starting in autumn 2019, the NDVI increases compared to the other years (especially for zone 6), and this can be related to the regular rainfall over the period November 2019 to May 2020.

* The highest interquartile range (IQR) values are observed for dates when the percentage of vegetation pixels per zone is less than 20%.

For natural vegetation zones (zones 3, 4 and 6), the NDVI is generally higher than for the phyto-stabilized zones, except for zone 3. Zone 3 has the lowest NDVI values until autumn 2017. Zone 6 has the highest NDVI values from autumn 2016 and the highest numbers of vegetation pixels in summer and autumn (Figure 12). Zone 6 is the zone with the most developed and tallest vegetation.

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Multi-annual evolution of NDVI statistics (median, Q1, Q3) by zone.</figcaption></figure>

For phyto-stabilized zones (zones 1, 2 and 5), the cycles are less pronounced and the differences in NDVI from one season to the next are smaller. There are areas for which the number of vegetation pixels remains greater than 20% whatever the season (Figure 9). Zone 5 has a different phenology compared to the other two zones. NDVI is high in autumn or winter (depending on the year) and decreases in spring. The in situ information indicates that zone 5 is the area with the highest species diversity.

In a second step, only indices with both spatial and temporal variabilities are chosen. For example, the SAVI index is selected (Figure 13a). The seasonal evolution of SAVI is different from that of NDVI owing to their formulations: in autumn SAVI decreases while NDVI increases (Figures 12 and 13b). The lowest SAVI values are then obtained in autumn, whereas those of NDVI are obtained in summer. NDVI allows the assessment of biomass importance and the monitoring of chlorophyll activity. It is often used for moderately dense canopies (no apparent soil but not too dense to avoid a saturation effect), without the mixing of standing dry matter with green matter. SAVI introduces an adjustment factor to minimize the ground effect. When this factor is equal to 0, SAVI and NDVI formulations are the same. SAVI is sensitive to soil color and brightness and is adapted in the case of a low chlorophyll vegetation cover [62,71,87]. NDVI and SAVI are then separated into two groups. This leads to the grouping of the indices with similar temporal evolution as follows:

* Group 1: NDVI (reference index), GREEN_NDVI, ARVI, CIREDEGE, MTCI, NDRE;

* Group 2: SAVI, MSAVI, IRECI, MCARI_NEW;

* Group 3: VII;

* Group 4: PSRI.

The SWIR index is not kept as the difference in spatial resolutions between VNIR and SWIR [7,54] should contribute to its disqualification. HSSMI (Table 2), specified for HM detection from Sentinel-2 images, is not used because of the very high IQR and RSD values (RSD between 50% and 100%).
<figure><img src="image_17.png" /><figcaption>**Figure 13.** (a) Seasonal evolution of SAVI (Soil Adjusted Vegetation Index) for the year 2017; (b) multi-annual evolution of SAVI statistics by zone: median, Q1 and Q3.</figcaption></figure>

These groups of indices are analyzed in detail below. Every index of group 1 is compared with the reference index belonging to this group. For Group 2, the SAVI index is compared to the reference index in order to check the complementarity of this index group. The other indices of Group 2 are then compared to SAVI. For groups 3 and 4 consisting of a single index, they are compared with the reference index.

Group 1 indices are compared to NDVI. Maps are visually compared as well as statistics. ARVI cannot be retained because its RSD values are twice as high as those of NDVI. The variation range of GREEN_NDVI is weaker than the one of NDVI, which leads to lower spatial variability and the temporal differences are less pronounced (Figure A1). GREEN_NDVI will not perform as well as NDVI in the detection of change. NDRE and NDVI provide the same information and the reference index is retained. CIREDEDGE presents a higher variation range, which highlights the temporal differences. This index, used to estimate chlorophyll, is retained in order to complete NDVI. MTCI supplies consistent information with CIREDEGE but provides higher RSD values. This analysis leads to the selection of CIREDEDGE in order to complete NDVI.

SAVI of Group 2 is compared to NDVI. Figure 4 compares the statistical values. In spring 2019 and the following seasons (exception of summer), for all zones, SAVI increases whilst NDVI varies differently depending on the considered zone. Zone 5 shows a temporal evolution of SAVI consistent with the zone 2 evolution, whereas the evolution of NDVI for zone 5 is not comparable to any other phyto-stabilized zone. The SAVI index is therefore complementary to the NDVI index.

The other indices of Group 2 are compared with SAVI. MSAVI is not retained because it provides the same information; the initial formulation of the SAVI index family is retained. The variation range of MCARI_NEW is lower and its RSD values are almost 1.7 times higher than those of NDVI. IRECI is of interest and is used to complete SAVI.

Group 3 consists of a single index PSRI specified for senescent vegetation. PSRI is sensitive to the ratio between carotenes and chlorophylls. It shows an inverse seasonal evolution in comparison to NDVI (Figures 12 and 14a): the lowest values are obtained in spring (in summer for NDVI) and the highest values in summer (in spring for NDVI). In addition, the areas with high NDVI values are the areas with low PSRI values (e.g., zone 6), which is explained by the role of each index. For zone 2, an increase in PSRI is observed for the date 2017/02/25, when the number of vegetation pixels is low. This difference is less pronounced when NDVI values are analyzed. The ranges variations of the IRQ and RSD values for PSRI and NDVI are similar. PSRI is complementary to NDVI.
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Multi-annual evolution of median indices by zone: (a) PSRI; (b) VII (Vegetation Inferiority Index).</figcaption></figure>

Group 4 consists of a single VII index developed for monitoring the state of vegetation growth in a mining area using hyperspectral data [80]. VII has been adapted because in the original formulation, it uses the green bands from 497 to 635 nm and the NIR and SWIR bands from 700 to 1200 nm. In our case, according to the 60 m spatial resolution of the two spectral bands B9 and B10 (Table 2), the second spectral band covers the domain from 704 nm to 865 nm. The VII seasonal evolution is different from the NDVI evolution. In autumn 2017 and for the following dates, zones 2, 3 and 4 present equivalent VII levels and global variations; the VII temporal evolution is comparable for zones 1, 5 and 6.

The high VII standard deviation values correspond to areas and dates with few pixels of vegetation. RSD values assessed from VII are lower for phyto-stabilized areas than those calculated from NDVI and, conversely, for areas with natural vegetation. This index is complementary to NDVI.

To conclude, the retained indices by group, each group corresponding to similar temporal behavior, are the following:

* Group 1: NDVI, CIREDEDGE;

* Group 2: SAVI, IRECI;

* Group 3: PSRI;

* Group 4: VII.

### 3.3. Change Detection Maps on Selected Indices

The inter-annual change detection maps (defined by Equation (1)) are provided for the selected indices and the two seasons with the highest percentage of vegetation cover (winter, spring) (Section 3.1). The reference year is 2016. The index maps of the years 2017, 2018, 2019 and 2020 are compared with those of 2016. The statistics (median, Q1, Q3, mean, standard deviation) of the change detection maps are calculated by zone.

Figure 15 shows the change detection maps for NDVI and IRECI indices. CIREDEDGE provides a better visibility for changes than NDVI does (Figure 15, Figure A3). Changes are more apparent in the IRECI detection map than in the SAVI map (Figure 15, Figure A3). PSRI highlights changes from one year to another (Figure A4). The changes are more difficult to see on the index VII maps (Figure A4). In spring, the vegetation cover remains stable over the years (Table 3), but significant variations for some indices depending on the area are observed. For examples, zones 4 and 6 have a significant diminution of IRECI values between 2016 and 2019, and all areas (except zone 4) present significant increase in the IRECI between 2016 and 2020.
Reference year: 2016

Group 1—NDVI

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

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

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

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

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

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

**Figure 15.** NDVI and IRECI (Inverted Red-Edge Chlorophyll Index) index change detection maps for Table 2017. 2019, 2020 compared to the year 2016. Maps of 2018 are not provided because the changes are not significant in comparison to the other years.

Statistics from the change detection maps are evaluated in spring and winter by zone. Figures 16 and 17 illustrate the mean for the three indices CIREDEDGE, IRECI and PSRI.

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

**Figure 16.** Median values of spring change detection maps on each zone for: (a) CIREDEDGE; (b) IRECI; (c) PSRI.

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

**Figure 17.** Median values of winter change detection maps on each zone for: (a) CIREDEDGE; (b) IRECI; (c) PSRI.

For every zone (except zone 4), the NDVI, CIREDEDGE, IRECI and SAVI values remain stable or even increase significantly from spring to winter between consecutive years until 2018 (positive values
< 0.1), decrease between spring 2018 and 2019 and increase until spring 2020 (deviation lower than 0.18). These variations are sensibly higher for zone 6. For zone 4, the same trend is observed but the deviation levels and temporal variations are higher. This point is confirmed in the field: this zone is undergoing regular changes. For zones 1 and 5 (phyto-stabilized), the indices present comparable temporal evolution and the evolution of the indices in phyto-stabilized zone 2 is the same as for the zone 3 with natural vegetation.

The PSRI values are roughly equal from spring 2016 to spring 2018. For zone 6, a slight decrease in index values is observed over the years until 2018. For zone 4, there is a gap between 2016 and 2017 and then stabilization until 2018. Between 2018 and 2019, PSRI values increase for every zone, which may reflect an increase in senescent vegetation during the spring season. Between 2019 and 2020, PSRI decreases and the median level values between spring 2018 and spring 2020 are comparable. This increase between 2018 and 2019 and decrease between 2019 and 2020 are twice as high for zones 4 and 6 with natural vegetation as well as for the other zones. This can be linked to the difference in terms of the dominant species in these areas. These statistics illustrate the abrupt change at the beginning of 2019, linked to the heavy rainfalls at the end of 2018.

The winter 2016 image is not available due to the presence of clouds. Only the comparison between winter 2017, winter 2019 and winter 2020 can then be performed (Figure 16).

For every area, there is either a decrease in the NDVI, CIREDEDGE, IRECI, SAVI index values or these values are stable between 2017 and 2018. Globally, the median values of NDVI, CIREDEDGE, IRECI and SAVI increase from 2018 to 2020. This increase rate depends on the area and the year. Zone Z4 (natural vegetation) has lower values of NDVI, CIREDEDGE and IRECI than the other zones (on the contrary, PSRI values are higher). With regards to the IRECI index, the phyto-stabilized zones (1, 2 and 5) show particular temporal behavior: the median values increase considerably until winter 2019 (deviation around 0.18) and then stabilizes. For natural vegetation zones, IRECI median values increase between 2019 and 2020. This can be seen on the 2019 and 2020 IRECI maps (Figure 18) as values around 0 in 2019 (pale beige color) are between 0.25 and 0.5 in 2020 (lime green color) for most pixels in these zones. There is a small variation in the PSRI values between 2017 and 2020. Senescent vegetation does not increase over the years in winter.

Reference year: 2017

Group 2—IRECI

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

<figure><img src="image_28.png" /><figcaption>**Figure 18.** IRECI index change detection maps for the winter of the years 2019 and 2020 compared to the year 2017.</figcaption></figure>

3.4. *Spectral and Temporal Clustering*

In a first time, the time series, called annual time series, is built with all the dates of the seasons winter, spring and autumn for each year. According to the previous results (see Section 3.1), the vegetation cover is very low in summer and the corresponding dates are not retained to construct the annual time series in order to avoid distorting vegetation pattern. The change detection is applied by year. The annual time series highlights major land use classes and provides multi-year monitoring. In a second time, the time series, called seasonal time series, is built by season (including the dates
of all the years for a given season) and the change detection is applied by season. The seasonal time series analysis allows highlighting phenological atypical behavior within zones.

Clustering is applied to each index of the previously selected indices and the annual time series. Change detection is then applied on the clustering maps. The results from Sections 3.1 and 3.2 are observed: bare soils are easily recognizable, vegetation areas are globally increasing and indices follow the trends observed Figure 16, slightly increasing from 2017 to 2018, stagnating or decreasing from 2018 to 2019 and sharply increasing, especially on zone 6, until 2020. Clustering in then applied on seasonal time series in order to distinguish phenological stages. The main results previously observed are: vegetation developing in spring and decreasing in summer, except in the south of zone 5 where vegetation decreases in spring and progresses in autumn, more apparent variations in natural zones than in phyto-stabilized zones and a greater vegetation cover in zones 1 and 5 than in zone 2 in summer and autumn. Clusters of sparse vegetation vary less than the others through the seasons. These results are consistent with those already obtained in Section 3.1 (Figure 8) and Section 3.3. IRECI, CIREDEDGE, NDVI and PSRI provide similar results. The main differences are (Figure A5):

* PSRI does not allow detecting change related to the senescence of the vegetation of the south of zone 5 in spring and autumn;

* IRECI and CIREDEDGE highlight greater progress in spring vegetation in the center of zone 6 than NDVI;

* The variability of the NDVI values is higher than those of other indices.

In order to take advantage of the sensibilities of spectral indices, clustering is then applied on combination of indices (Section 3.2) and the annual time series. The combination of NDVI and PSRI is retained to highlight the observed changes. Figure 19 illustrates the clustering map for each year. The size of bare soil clusters (in red) decreases while vegetation clusters (in green) are expending over the years. This progress is even more important over the entire zone 6, in the center of zone 3, in the south of zone 2 and in the south of zone 5 where areas of well-developed vegetation (dark green) spread out. Thus, both natural and phyto-stabilized vegetation seems to develop over the years. An area of bare soil is observed in the northeast of zone 1 in 2018 only. Those of zone 4 are constantly changing, a result which is consistent with those already observed (Section 3.3, Figure 15).

<figure><img src="image_29.png" /><figcaption>**Figure 19.** Clustering maps obtained with the combination of NDVI and PSRI for annual time series (built with dates in winter, spring and autumn of each year).</figcaption></figure>

The spectral and temporal clustering is applied on the combination of the selected indices (CIREDEDGE, IRECI and PSRI), the reference index NDVI and the seasonal time series. The change detection maps compare season by season (Figure 20). In the south of zone 5, a phenological stage change is observed (surrounded in blue) and can be explained by the presence of grass land drying in late spring with rising temperatures and returning in autumn with rainfalls (Figure 9). In the north of zone 6 (natural zone), surrounded in black, the cluster maps always group this area in the high vegetation cluster. The index values for this area vary only slightly between spring and autumn (brown color) and increase sharply in spring (dark green color). This area seems to be composed of persistent vegetation and field observations confirm this result: the vegetation is mainly composed of Spanish brooms and pines. The vegetation in the south of zone 2 weakens strongly in winter (surrounded in red), by falling from a cluster of high vegetation to one of sparse vegetation in winter.
The east of zone 5 (surrounded in orange) varies more strongly in winter and spring than the rest of the vegetation, passing from cluster of spare vegetation to high vegetation.

<figure><img src="image_30.png" /><figcaption>**Figure 20.** Clustering change detection maps using NDVI, CIREDEDGE, IRECI and PSRI for seasonal time series (combining all the dates of each season) (notation: season 1/season 2 means change detection map of the season 1 according to the season 2).</figcaption></figure>

**4. Discussion**

The main purpose of this study is to analyze the temporal evolution of vegetation growing in an ore processing site after closure and detect the abrupt change related to one-time events (for example related to meteorological conditions) and the longer-term change related to vegetation response to stress associated to the environmental conditions of development. The vegetation is divided into two classes in this study: vegetation introduced for phyto-stabilization and natural vegetation. For that purpose, the data of Sentinel-2 satellites, several spectral indices and change detection methods are used to provide an unsupervised diagnostic tool. The time series exploited, covering four seasons and five years, leads to large volumes of images and maps to be manipulated. The change detection tool proposed, based on multi-temporal fusion at feature level and involving the normalized differences of the index maps at two dates, allows a faster and more relevant analysis of the results. The clustering approach highlights some common behaviors across the time series and could be a good way to support change detection by emphasizing general trends and unexpected variations. The results of this investigation are carried out on the basis of the vegetation classes (phyto-stabilized and natural vegetation) by qualitative and quantitative analysis of the produced maps related to spectral index or representative of vegetation cover.

Among the 26 spectral indices tested from the literature covering different contexts (environmental monitoring, precision agriculture abiotic impact...), three indices, CIREDEDGE (based on spectral bands B5 and B7 and used for the estimation of chlorophyll) [64], IRECI (including the bands B4, B5, B6 and B7 in its formulation and specified for the estimation of chlorophyll content from Sentinel-2 data) [7] and PSRI (using the bands B2, B4 and B6 and defined for senescent vegetation, sensitive to the ratio of carotenoids, carotenes and chlorophylls) [78] are particularly suitable for detecting spatial and temporal changes over the studied area. These indices usually contain a wavelength sensitive to chlorophyll (often between 660 and 720 nm) and a reference wavelength (between 750 and 800 nm) to account for vegetation structure [18]. The common spectral bands used by these three indices are defined below with the purpose for each band: B4 (maximum chlorophyll absorption), B5 (REP—Red-Edge Position), B6 (REP) and B7 (LAI—Leaf Area Index, edge of the NIR plateau) [88]. PSRI introduces the B2 band, which is sensitive to vegetation senescing, carotenoid and browning. Many previous studies employed red-edge information in order to detect spectral changes due to the high soil metal contamination [9,18,89]: our index selection confirms the usefulness of this spectra region given an indication of the chlorophyll plant status. The signal in the Red Edge as a result of the increasing chlorophyll content is due to a broadening of the absorption feature centered at 680 nm, shifting the position of the red-edge to longer wavelengths [90]. The shift in the red-edge position towards shorter wavelengths around 680 nm is then associated with a decrease in chlorophyll content. CIREDEDGE and IRECI are related to chlorophyll contents in leaf indirectly impacted by the presence of HM in the
soil [4,9,18]. A high value of such indices indicates high chlorophyll content and healthy plant status. In addition to this, CIREDEDGE and PSRI are proposed as reference indices in the study conducted by Zhang and al for detecting heavy metal stress by using Sentinel-2 data [27]. IRECI has been shown to be correlated to canopy chlorophyll content and still performing well for LAI [7]. The vegetation indices selected allow capturing modifications of vegetation traits (mainly biochemical traits and LAI) even if minor vegetation cover changes are observed. In fact, during spring the vegetation cover remains relatively stable over the year, but specific spatial variations for each zone on the index map are detected and highlighted by the temporal statistical analysis.

Identifying and predicting phenological stages provide essential information for many applications. Some phenological stages are more sensitive to biotic or abiotic impacts. The Sentinel-2 time series offers a good opportunity to analyze phenological stages and transition dates [89,91,92]. The high variability of the vegetation cover cycle is analyzed in this study to detect change. Overall and abrupt changes of the vegetation cycle concerning every area, linked to heavy rainfalls, are highlighted. Vegetation development depends on meteorological conditions such as temperatures and rain fall. In special cases (floods, heatwaves, frost...), selected dates to construct the time series could be adjusted annually based on meteorological observations. A difference in the vegetation cycle across the area is observed (vegetation cover increases in autumn from 2017), likely related to the relative abundance of the different plant species. The most appropriate seasons to monitor vegetation on this study site are defined according to vegetation cover: winter and spring.

The spatial resolution offered by S2 allows to separate and analyze single areas, even when their size is relatively small (few hectares), with a high temporal resolution. Even if temporal and spectral clustering shows a reduction in bare soil surface over the years, some permanent bare soil areas on three zones without vegetation cover whatever the year and season since 2016 are highlighted. Permanent bare soils seem to be related to limestone soils that are poor in organic matter. Further in situ investigations (spectral signature measurements of bare soil area with very marked metal signature, soil texture analysis...) have to be conducted in order to identify the major factor preventing the development of vegetation (HM concentrations, high pH, low organic matter content...) and define if this factor can be detected or quantified via optical properties. Other specific changes are identified by the tool proposed—for example, the identification of the most impacted area or, conversely, the detection of the area with vigorous vegetation (most developed trees).

Only a subtle difference has been detected between phyto-stabilized and spontaneous vegetation areas (like differences noticed by IRECI in winter). No global trend can be identified for the phyto-stabilized areas and the zones with natural vegetation. This can be explained by the difference in terms of dominant species in these areas or by the differences in management methods for the aided phyto-stabilisation (partial phyto-stabilisation, tree planting, addition of natural soil, etc.). Seasonal and annual clustering gives information on land cover (bare soil, dense or sparse vegetation) and habitats (persistent, growing cycle-dependent). To remove certain ambiguities, major species and habitats on each zone have to be identified. Moreover, control zones spatially independent from the study site with similar species, soil mineralogy and exposition could be identified and the vegetation development of these control zones will be compared to the zones affected by the ore processing activity.

Most studied vegetation indices use spectral wavelengths alongside the red-edge region, sensitive to changes in chlorophyll and other pigments. This can lead to false alarms as the changes due to other biological and physiological stressors inducing chlorophyll alterations can be confused with the changes related to metal accumulation [89]. Excessive metal accumulation is known to affect leaf internal structure, cell structure and chloroplast ultrastructure [89,93], information available in the spectral domain (800–1300 nm). The combination of the red-edge domain with the SWIR domain can reduce false alarms and improve the detection of HM alteration. To further establish a direct link between HM in the soil and the health status of vegetation, in situ investigations are necessary and measurements of spectral signature of the major species with a field spectro-radiometer complemented along with measurements of leaf biochemical parameters (pigment and water contents...) and HM
concentrations in soil have to be acquired in order to define the most sensitive and useful spectral bands in the whole reflective domain (400–2500 nm) according to the context defined by the HM type and level as well as the plant species [11].

Previous studies have proven the interest of hyperspectral data (typically airborne or in-field measurements) for monitoring HM stress [9,26,94]. On the one hand, the approach, proposed to assess the potential of leaf optical properties for predicting Total Petroleum Hydrocarbon concentrations and based on the retrieval of leaf pigment using the PROSPECT model, could be applied in this context on a range of plant species over a full seasonal cycle in order to identify the most suitable conditions for detecting and predicting HM concentrations [10,19]. On the other hand, it would be interesting to process an airborne hyperspectral image acquired on the site in order to analyze the contribution of spectral and high spatial (metric) resolutions in comparison to temporal repeatability.

**5. Conclusions**

An unsupervised methodology for monitoring vegetation on an ore processing site closed and revegetated is proposed. It operates multispectral and multi-temporal Sentinel-2 images and generates change detection maps based on the exploitation of spectral indices and spectral and temporal clustering. Three indices, CIREDEDGE, IRECI and PSRI, related to the biochemical traits of vegetation are particularly suitable for detecting changes spatially and temporally across the study area.

Image time series reveal its capacity for monitoring vegetation in such a context, and this study proves that the multi-date analysis are important to emphasize difference in vegetation development and vegetation trait variation. Global and specific changes are emphasized by the exploitation of spectral index maps, vegetation masks related to vegetation cover, clustering maps and change detection maps. The major changes detected are explained by information provided by the site operator or meteorological conditions.

The results are specific to the study site owing to its history, the phyto-stabilization process and the species in development but this unsupervised diagnostic tool remains valid for other contexts such as phyto-management monitoring, the analysis of restored or polluted sites or precision agriculture.

Further study is needed in order to better understand the behavior difference between the vegetation introduced for phyto-stabilization and the natural vegetation, as well as the impact of trace metals on the vegetation health. In a first time, species or habitat mapping will be of great interest. Supervised classification methods based on machine learning algorithms could then be applied to the Sentinel-2 time series. In a second time, it seems necessary to carry out multi-temporal field measurements with a spectro-radiometer combined with measurements of biochemical traits of the major species and soil parameters (texture, HM, pH...) on contaminated and control zones in order to define the combination of wavelengths optimal to characterize the metal accumulation and/or metal stress in plants.

According to the relatively small size of the study areas, the spatial distribution of vegetation, the size of the plant patches and the major species, spatial information must be exploited through high spatial resolution imagery. Unmanned aerial vehicle (UAV) and airborne platforms taking hyperspectral camera on board could then be used over time to monitor the temporal dynamics of vegetation, map species and detect altered vegetation due to soil HM by combining high spatial, spectral and spatial resolution information.

**Author Contributions:** Conceptualization, S.F.; methodology, S.F.; software, S.F., R.G. and T.R.; validation, S.F., R.G. and A.E.; formal analysis, S.F. and A.E.; data curation, S.F.; writing—original draft preparation, S.F., A.E. and T.R.; supervision, S.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by CNES.

**Acknowledgments:** This work was performed in the frame of the COMPOST project between ONERA and EcoLab. The authors gratefully acknowledge the BRGM.

**Conflicts of Interest:** The authors declare no conflict of interest.
Appendix A

This part provides supplementary results to complete the results analysis and interpretation.

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

**Figure A1.** Multi-annual evolution of statistics (median, Q1, Q3) by zone: (a) NDVI; (b) GREEN_NDVI. The variation range of GREEN_NDVI is lower in comparison to NDVI's one, which leads to a lower spatial variability, and the temporal differences are less pronounced.

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

**Figure A2.** Multi-annual evolution of statistics (median, Q1, Q3) by zone: (a) MCARI_NEW; (b) IRECI. The variation range of MCARI_NEW is lower (mean values between 0.02 and 0.18) than those of NDVI. IRECI presents the temporal variability of interest for change detection.

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

**Figure A3.** SAVI and CIREDEDGE index change detection maps for the spring of the years 2017, 2019, 2020 compared to the year 2016. Maps of 2018 are not provided because the changes are not significant in comparison to the other years. CIREDEDGE provides a better visibility of changes than NDVI. Changes are more apparent in the IRECI detection map than in the SAVI map.
<figure><img src="image_34.png" /><figcaption>**Figure A4.** VII and PSRI index change detection maps for the spring of the years 2017, 2019, 2020 compared to the year 2016. Maps of 2018 are not provided because changes are not significant in comparison to the other years. The changes are more difficult to see in the index VII maps than in the other index maps. PSRI highlights changes by exploiting the full range of colors.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure A5.** Clustering change detection maps using for seasonal time series (spring according to winter). Main differences observed: PSRI does not allow detecting change related to vegetation senescence of the south of zone 5 in spring and autumn; IRECI and CIREDEDGE highlight greater progress in spring vegetation in the center of zone 6 than NDVI; the variability of the NDVI values is higher than that of other indices.</figcaption></figure>
## References

1. Verbesselt, J.; Hyndman, R.; Newnham, G.; Culvenor, D. Detecting trend and seasonal changes in satellite image time series. *Remote Sens. Environ.* **2010**, *114*, 106–115. [CrossRef]

2. Barceló, J.; Poschenrieder, C. Plant water relations as affected by heavy metal stress: A review. *J. Plant. Nutr.* **1990**, *13*, 1–37. [CrossRef]

3. Slonecker, T.; Fisher, G.B.; Aiello, D.P.; Haack, B. Visible and infrared remote imaging of hazardous waste: A review. *Remote Sens.* **2010**, *2*, 2474–2508. [CrossRef]

4. Ong, C.; Carrère, V.; Chabrillat, S.; Clark, R.; Hoefen, T.; Kokaly, R.; Marion, R.; Souza Filho, C.R.; Swayze, G.; Thompson, D.R. Imaging Spectroscopy for the Detection, Assessment and Monitoring of Natural and Anthropogenic Hazards. *Surv. Geophys.* **2019**, *40*, 431–470. [CrossRef]

5. Nie, M.; Lu, M.; Yang, Q.; Zhang, X.-D.; Xiao, M.; Jiang, L.-F.; Yang, J.; Fang, C.-M.; Chen, J.-K.; Li, B. Plants’ use of different nitrogen forms in response to crude oil contamination. *Environ. Pollut.* **2011**, *159*, 157–163. [CrossRef]

6. Barati, S.; Rayegani, B.; Saati, M.; Sharifi, A.; Nasri, M. Comparison the accuracies of different spectral indices for estimation of vegetation cover fraction in sparse vegetated areas. *Egypt. J. Remote Sens. Space Sci.* **2011**, *14*, 49–56. [CrossRef]

7. Frampton, W.J.; Dash, J.; Watmough, G.; Milton, E.J. Evaluating the capabilities of Sentinel-2 for quantification estimation of biophysical variables in vegetation. *ISPRS J. Photogramm. Remote Sens.* **2013**, *82*, 83–92. [CrossRef]

8. Zinnert, J.C.; Via, S.M.; Young, D.R. Distinguishing natural from anthropogenic stress in plants: Physiology, fluorescence and hyperspectral reflectance. *Plant. Soil* **2011**, *366*, 133–141. [CrossRef]

9. Slonecker, E.T. Analysis of the Effects of Heavy Metals on Vegetation Hyperspectral Reflectance Properties from Advanced Applications in Remote Sensing of Agricultural Crops and Natural Vegetation. In *Hyperspectral Remote Sensing of Vegetation*, 2nd ed.; Thenkabail, P.S., Lyon, J.G., Huete, A., Eds.; CRC Press: London, UK, 2018; Volume IV, pp. 49–66.

10. Lassalle, G.; Fabre, S.; Credoz, A.; Hédacq, R.; Borderies, P.; Bertoni, G.; Erudel, T.; Buffan-Dubaud, E.; Dubucq, D.; Elger, A. Detection and discrimination of various oil mixtures in soils using vegetation indices: A multi-scale approach. *Sci. Total Environ.* **2019**, *655*, 113–1124. [CrossRef]

11. Lassalle, G.; Fabre, S.; Credoz, A.; Hédacq, R.; Bertoni, G.; Dubucq, D.; Elger, A. Application of PROSPECT for estimating Total Petroleum Hydrocarbons in contaminated soils from leaf optical properties. *J. Hazard. Mater.* **2019**, *377*, 409–417. [CrossRef]

12. Gholizadeh, A.; Kopačková, V. Detecting vegetation stress as a soil contamination proxy: A review of optical proximal and remote sensing techniques. *Int. J. Environ. Sci. Technol.* **2019**, *16*, 2511–2524. [CrossRef]

13. Mendez, M.O.; Maier, R.M. Phytostabilization of Mine Tailings in Arid and Semiarid Environments—An Emerging Remediation Technology. *Environ. Health Perspect.* **2008**, *116*, 278–283. [CrossRef]

14. Davids, C.; Rouyet, L. *Remote Sensing for the Mining Industry*; Report, Project RESEM; Northern Research Institute: Tromsø, Norway, 2018.

15. Navarro, M.C.; Pérez-Sirvent, C.; Martínez-Sánchez, M.J.; Vidal, J.; Tovar, P.J.; Bech, J. Abandoned mine sites as a source of contamination by heavy metals: A case study in a semi-arid zone. *J. Geochem. Explor. Geol.* **2008**, *96*, 183–193. [CrossRef]

16. Escarré, J.; Lefèbvre, C.; Raboyeau, S.; Dossantos, A.; Gruber, W.; Marel, J.C.C.; Frérot, H.; Noret, N.; Mahieu, S.; Collin, C.; et al. Heavy Metal Concentration Survey in Soils and Plants of the Les Malines Mining District (Southern France): Implications for Soil Restoration. *Water Air Soil Pollut.* **2011**, *216*, 485–504. [CrossRef]

17. Sun, Z.; Xie, X.; Wang, P.; Hu, Y.; Cheng, H. Heavy metal pollution caused by small-scale metal ore mining activities: A case study from a polymetallic mine in South China. *Sci. Total Environ.* **2018**, *639*, 217–227. [CrossRef] [PubMed]

18. Sims, D.A.; Gamon, J.A. Relationships between Leaf Pigment Content and Spectral Reflectance across a Wide Range of Species, Leaf Structures and Developmental Stages. *Remote Sens. Environ.* **2002**, *81*, 337–354. [CrossRef]

19. Jacquemoud, S.; Baret, F. PROSPECT: A Model of Leaf Optical Properties Spectra. *Remote Sens. Environ.* **1990**, *34*, 75–91. [CrossRef]
20. Slonecker, T.; Haack, B.; Price, S. Spectroscopic analysis of arsenic uptake in in Pteris ferns. *Remote Sens.* **2008**, *1*, 644–675. [CrossRef]

21. Horler, D.; Barber, J.; Barringer, A.; Njoku, E. Effects of heavy metals on the absorbance and reflectance spectra of plants. *Int. J. Remote Sens.* **1980**, *1*, 121–136. [CrossRef]

22. Clark, R.N.; Roush, T.L. Reflectance spectroscopy: Quantitative analysis techniques for remote sensing applications. *J. Geophys. Res. Solid Earth* **1984**, *89*, 6329–6340. [CrossRef]

23. Beck, P.S.A.; Jönsson, P.; Høgda, K.-A.; Karlsen, S.R.; Eklundh, L.; Skidmore, A.K. A ground validated NDVI dataset for monitoring vegetation dynamics and mapping phenology in Fennoscandia and the Kola Peninsula. *Int. J. Remote Sens.* **2007**, *28*, 4311–4330. [CrossRef]

24. Pearson, R.L.; Miller, L.D. Remote mapping of standing crop biomass for estimation of the productivity of the shortgrass prairie, Pawnee national grasslands, Colorado. In Proceedings of the Eighth International Symposium on Remote Sensing of Environment, Ann Arbor, MI, USA, 2–6 October 1972; p. 1355.

25. Xue, J.; Su, B. Significant Remote Sensing Vegetation Indices: A Review of Developments and Applications. *J. Sens.* **2017**, *2017*, 1353691. [CrossRef]

26. Zhou, C.; Chen, S.; Zhang, Y.; Zhao, J.; Song, D.; Liu, D. Evaluating Metal Effects on the Reflectance Spectra of Plant Leaves during Different Seasons in Post-Mining Areas, China. *Remote Sens.* **2018**, *10*, 1211. [CrossRef]

27. Zhang, Z.; Liu, M.; Liu, X.; Zhou, G. A New Vegetation Index Based on Multitemporal Sentinel-2 Images for Discriminating Heavy Metal Stress Levels in Rice. *Sensors* **2018**, *18*, 2172. [CrossRef] [PubMed]

28. Liu, M.; Wang, T.; Skidmore, A.K.; Liu, X. Heavy metal-induced stress in rice crops detected using multi-temporal Sentinel-2 satellite images. *Sci. Total Environ.* **2018**, *637–638*, 18–29. [CrossRef] [PubMed]

29. Gonçalves, R.; Zullo, J., Jr.; Amaral, B.; Coltri, P.; Sousa, E.; Romani, L. Land use temporal analysis through clustering techniques on satellite image time series. In Proceedings of the 2014 IEEE Geoscience and Remote Sensing Symposium, Quebec City, QC, Canada, 13–18 July 2014; pp. 2173–2176. [CrossRef]

30. Espinoza-Molina, D.; Bahmanyar, R.; Bustamante, J.; Datcu, M.; Diaz-Delgado, R. Land-Cover Change Detection Using Local Feature Descriptors Extracted from Spectral Indices. In Proceedings of the IEEE IGARSS Conference, Fort Worth, TX, USA, 23–28 July 2017. [CrossRef]

31. Gonçalves, R.; Junior, J.; Amaral, B.; Sousa, E.; Romani, L. *Agricultural Monitoring in Regional Scale Using Clustering on Satellite Image Time Series*; IntechOpen: London, UK, 2018. [CrossRef]

32. Kaufman, L.; Rousseeuw, P. *Finding Groups in Data: An. Introduction to Cluster Analysis*; John Wiley & Sons: New York, NY, USA, 2009.

33. Xu, R.; Wunsch, D. Survey of clustering algorithms. *IEEE Trans. Neural Netw.* **2005**, *16*, 645–678. [CrossRef]

34. Han, J.; Kamber, M. *Data Mining—Concepts and Techniques*, 2nd ed.; Morgan Kaufmann Publishers: New York, NY, USA, 2006.

35. Zhao, Y.; Yuan, Y.; Wang, Q. Fast Spectral Clustering for Unsupervised Hyperspectral Image Classification. *Remote Sens.* **2019**, *11*, 399. [CrossRef]

36. Wang, Z.; Xia, G.; Xiong, C.; Zhang, L. Spectral active clustering of remote sensing images. In Proceedings of the IEEE Geoscience and Remote Sensing Symposium, Quebec City, QC, Canada, 13–18 July 2014; pp. 1737–1740. [CrossRef]

37. Pascucci, S.; Carfora, M.F.; Palombo, A.; Pignatti, S.; Casa, R.; Pepe, M.; Castaldi, F.A. Comparison between Standard and Functional Clustering Methodologies: Application to Agricultural Fields for Yield Pattern Assessment. *Remote Sens.* **2018**, *10*, 585. [CrossRef]

38. Sheeren, D.; Fauvel, M.; Josipovic, V.; Lopes, M.; Planque, C.; Willm, J.; Dejoux, J.F. Tree Species Classification in Temperate Forests Using Formosat-2 Satellite Image Time Series. *Remote Sens.* **2016**, *8*, 734. [CrossRef]

39. Huang, X.; Liu, J.; Zhu, W.; Atzberger, C.; Liu, Q. The Optimal Threshold and Vegetation Index Time Series for Retrieving Crop Phenology Based on a Modified Dynamic Threshold Method. *Remote Sens.* **2019**, *11*, 2725. [CrossRef]

40. Forkel, M.; Carvalhais, N.; Verbesselt, J.; Mahecha, M.D.; Neigh, C.S.R.; Reichstein, M. Trend Change Detection in NDVI Time Series: Effects of Inter-Annual Variability and Methodology. *Remote Sens.* **2013**, *5*, 2113–2144. [CrossRef]

41. Ohana-Levi, N.; Paz-Kagan, T.; Panov, N.; Peeters, A.; Tsoar, A.; Karnieli, A. Time series analysis of vegetation-cover response to environmental factors and residential development in a dryland region. *GI Sci. Remote Sens.* **2019**, *56*, 362–387. [CrossRef]
42. Drusch, M.; Del Bello, U.; Carlier, S.; Colin, O.; Fernandez, V.; Gascon, F.; Hoersch, B.; Isola, C.; Laberinti, P.; Martimort, P.; et al. Sentinel-2: ESA's optical high-resolution mission for GMES operational services. *Remote Sens. Environ.* **2012**, *120*, 25–36. [CrossRef]

43. Tian, L.; Liu, X.; Zhang, B.; Liu, M.; Wu, L. Extraction of rice heavy metal stress signal features based on long time series leaf area index data using ensemble empirical mode decomposition. *Int. J. Environ. Res. Public Health* **2017**, *14*, 1018. [CrossRef]

44. Bruzzone, L.; Bolovo, F. A Novel Framework for the Design of Change-Detection Systems for Very-High-Resolution Remote Sensing Images. *Proc. IEEE* **2012**, *101*, 609–630. [CrossRef]

45. Hussain, M.; Chen, D.; Cheng, A.; Wei, H.; Stanley, D. Change detection from remotely sensed images: From pixel-based to object-based approaches. *ISPRS J. Photogramm. Remote Sens.* **2013**, *80*, 91–106. [CrossRef]

46. Bolovo, F.; Bruzzone, L. The Time Variable in Data Fusion: A Change Detection Perspective. Data fusion in Remote Sensing. *IEEE Geosci. Remote Sens. Mag.* **2015**, *3*, 8–26. [CrossRef]

47. Ghamisi, P.; Rasti, B.; Yokoya, N.; Wang, Q.; Höfle, B.; Bruzzone, L.; Bovolo, F.; Chi, M.; Anders, K.; Gloaguen, R.; et al. Multisource and Multitemporal Data Fusion in Remote Sensing. A comprehensive review of the state of the art. *IEEE Geosci. Remote Sens. Mag.* **2019**, *7*, 6–39. [CrossRef]

48. Devi, R.N.; Jiji, G.W. Change detection techniques—A survey. *Int. J. Comput. Sci. Appl.* **2015**, *5*, 45–57.

49. Bhavani, M.; Sangeetha, V.H.; Kalaivani, K.; Ulagapriya, K.; Saritha, A. Change detection algorithm for multi-temporal satellite images: A review. *Int. J. Eng. Technol.* **2018**, *7*, 206–209. [CrossRef]

50. Melchori, A.E.; de Almeida Cândido, P.; Libonati, R.; Morelli, F.; Setzer, A.W.; de Jesus, S.C.; Garcia Fonseca, L.; Körting, T.S. Spectral indices and multi-temporal change image detection algorithms for burned area extraction in the Brazilian Cerrado. In Proceedings of the Anais XVII Simpósio Brasileiro de Sensoriamento Remoto—SBSR, JoãoPsso, Brasil, 25–29 April 2015.

51. Lu, D.; Mausel, P.; Brondizios, E.; Moran, E. Change detection techniques. *Int. J. Remote Sens.* **2004**, *25*, 2365–2407. [CrossRef]

52. Polykretis, C.; Grillakis, M.G.; Alexakis, D.D. Exploring the Impact of Various Spectral Indices on Land Cover Change Detection Using Change Vector Analysis: A Case Study of Crete Island, Greece. *Remote Sens.* **2020**, *12*, 319. [CrossRef]

53. Theia Data and Services Center for Continental Surfaces. Available online: https://www.theia-land.fr/pole-theia-2/ (accessed on 20 January 2020).

54. Louis, J.; Debaecker, V.; Pflug, B.; Main-Korn, M.; Bieniarz, J.; Mueller-Wilm, U.; Cadau, E.; Gascon, F. Sentinel-2 Sen2cor: L2A processor for users. In Proceedings of the Living Planet Symposium 2016, Prague, Czech Republic, 9–13 May 2016; Volume SP-740.

55. Martimort, P.; Fernandez, V.; Kirschner, V.; Isola, C.; Meygret, A. Sentinel-2 MultiSpectral imager (MSI) and calibration/validation. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, Munich, Germany, 22–27 July 2012. [CrossRef]

56. GDAL/OGR Contributors. GDAL/OGR Geospatial Data Abstraction Software Library. *Open Source Geospatial Foundation*. 2019. Available online: https://gdal.org (accessed on 20 January 2020).

57. Waqar, M.M.; Mirza, J.F.; Mumtaz, R.; Hussain, E. Development of new indices for extraction of built-up area & bare soil from Landsat data. *Open Access Sci. Rep.* **2012**, *1*, 2–4.

58. Valdiviezo-N, J.C.; Téllez-Quiñones, A.; Salazar-Garibay, A.; López-Caloca, A.A. Built-up index methods and their applications for urban extraction from Sentinel 2A satellite data: Discussion. *J. Opt. Soc. Am. A* **2018**, *35*, 35–44. [CrossRef] [PubMed]

59. Du, Y.; Zhang, Y.; Ling, F.; Wang, Q.; Li, W.; Li, X. Water Bodies' Mapping from Sentinel-2 Imagery with Modified Normalized Difference Water Index at 10-m Spatial Resolution Produced by Sharpening the SWIR Band. *Remote Sens.* **2016**, *8*, 354. [CrossRef]

60. Huete, A.R. A soil-adjusted vegetation index (SAVI). *Remote Sens. Environ.* **1988**, *25*, 295–309. [CrossRef]

61. Rasul, A.; Balzter, H.; Faqe Ibrahim, G.R.; Hameed, H.M.; Wheeler, J.; Adamu, B.; Ibrahim, S.; Najmaddin, P.M. Applying Built-Up and Bare-Soil Indices from Landsat 8 to Cities in Dry Climates. *Land* **2018**, *7*, 81. [CrossRef]

62. Kaufman, Y.J.; Tanré, D. Atmospherically Resistant Vegetation Index (ARVI) for EOS-MODIS. *IEEE Trans. Geosci. Remote Sens.* **1992**, *30*, 261–270. [CrossRef]

63. Filella, I.; Peñuelas, J. The red-edge position and shape as indicators of plant chlorophyll content, biomass and hydric status. *Int. J. Remote Sens.* **1994**, *15*, 1459–1470. [CrossRef]
64. Gitelson, A.A.; Gritz, Y.; Merzlyak, M.N. Relationships between leaf chlorophyll content and spectral reflectance and algorithms for non-destructive chlorophyll assessment in higher plant leaves. *J. Plant. Physiol.* **2003**, *160*, 271–282. [CrossRef]

65. Carter, G.A. Ratios of leaf reflectances in narrow wavebands as indicators of plant stress. *Int. J. Remote Sens.* **1994**, *15*, 697–703. [CrossRef]

66. Emengini, E.J.; Blackburn, G.A.; Theobald, J.C. Discrimination of plant stress caused by oil pollution and waterlogging using hyperspectral and thermal remote sensing. *J. Appl. Remote Sens.* **2013**, *7*, 073476. [CrossRef]

67. Smith, K.L.; Steven, M.D.; Colls, J.J. Use of hyperspectral derivative ratios in the red-edge region to identify plant stress responses to gas leaks. *Remote Sens. Environ.* **2004**, *92*, 207–217. [CrossRef]

68. Louhaichi, M.; Borman, M.M.; Johnson, D.E. Spatially located platform and aerial photography for documentation of grazing impacts on wheat. *Geocarto Int.* **2001**, *16*, 65–70. [CrossRef]

69. Gitelson, A.A.; Chivkunov, O.B.; Merzlyak, M.N. Nondestructive estimation of anthocyanins and chlorophylls in anthocyanic leaves. *Am. J. Bot.* **2009**, *96*, 1861–1868. [CrossRef] [PubMed]

70. Daughtry, C.; Walthall, C.L.; Kim, M.S.; Brown de Colstoun, E.C.; McMurtrey, J.E. Estimating Corn Leaf Chlorophyll Concentration from Leaf and Canopy Reflectance. *Remote Sens. Environ.* **2000**, *74*, 229–239. [CrossRef]

71. Haboudane, D.; Miller, J.R.; Pattey, E.; Zarco-Tejada, P.J.; Strachan, I.B. Hyperspectral Vegetation Indices and Novel Algorithms for Predicting Green LAI of Crop Canopies: Modeling and Validation in the Context of Precision Agriculture. *Remote Sens. Environ.* **2004**, *90*, 337–352. [CrossRef]

72. Wu, C.Y.; Niu, Z.; Tang, Q.; Huang, W.J. Estimating chlorophyll content from hyperspectral vegetation indices: Modeling and validation. *Agric. For. Meteorol.* **2008**, *148*, 1230–1241. [CrossRef]

73. Richardson, A.J.; Wiegand, C. Distinguishing vegetation from soil background information. *Photogramm. Eng. Remote Sens.* **1977**, *43*, 1541–1552.

74. Dash, J.; Curran, P.J. The MERIS terrestrial chlorophyll index. *Int. J. Remote Sens.* **2004**, *25*, 5403–5413. [CrossRef]

75. Keeley, J.E. Fire intensity, fire severity and burn severity: A brief review and suggested usage. *Int. J. Wildland Fire* **2009**, *18*, 116–126. [CrossRef]

76. Barnes, E.M.; Clarke, T.R.; Richards, E.; Colaizzi, D.; Haberland, J.; Kostrzewski, M.; Waller, P.; Choi, C.; Riley, E.; Thompson, T.; et al. Coincident detection of crop water stress, nitrogen status and canopy density using ground based multispectral data. In Proceedings of the Fifth International Conference on Precision Agriculture, Bloomington, MN, USA, 16–19 July 2000.

77. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, W.D. Monitoring vegetation systems in the Great Plains with ERTS. In Proceedings of the Third ERTS Symposium, Washington, DC, USA, 10–14 December 1973; Volume SP-351, pp. 309–317.

78. Merzlyak, M.N.; Gitelson, A.A.; Chivkunova, O.B.; Rakitin, V.Y. Non-destructive optical detection of pigment changes during leaf senescence and fruit ripening. *Physiol. Plant.* **1999**, *106*, 135–141. [CrossRef]

79. Haboudane, D.; Miller, J.R.; Tremblay, N.; Zarco-Tejada, P.J.; Dextraze, L. Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture. *Remote Sens. Environ.* **2002**, *81*, 416–426. [CrossRef]

80. Zhang, B.; Wu, L.; Zhang, L.; Jiao, Q.; Li, Q. Application of hyperspectral remote sensing for environment monitoring in mining areas. *Environ. Earth Sci.* **2012**, *65*, 649–658. [CrossRef]

81. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; et al. Scikit-learn: Machine Learning in Python. *JMLR* **2011**, *12*, 2825–2830.

82. Zhang, L.; You, J. A spectral clustering based method for hyperspectral urban image. In Proceedings of the 2017 Joint Urban Remote Sensing Event (JURSE), Dubai, UAE, 6–8 March 2017; pp. 1–3. [CrossRef]

83. Rousseeuw, P.J. Silhouettes A graphical aid to the interpretation and validation of cluster analysis. *J. Comput. Appl. Math.* **1987**, *20*, 53–65. [CrossRef]

84. Anselin, L. Local indicators of spatial association—LISA. *Geogr. Anal.* **1995**, *27*, 93–115. [CrossRef]

85. Al-doski, J.; Mansor, B. NDVI Differencing and Post-classification to Detect Vegetation Changes in Halabja City, Iraq. *IOSR JAGG* **2013**, *1*, 1–10. [CrossRef]

86. Everitt, B.; Skrondal, A. *The Cambridge Dictionary of Statistics*, 4th ed.; Cambridge University Press: Cambridge, UK, 2010.
87. Vani, V.; Mandla, V.R. Comparative study of NDVI and SAVI vegetation indices in Anantapur district semi-arid areas. *Int. J. Civ. Eng. Technol.* **2017**, *8*, 559-566.

88. Vanino, S.; Nino, P.; De Michele, C.; Bolognesi, S.F.; D'Urso, G.; Di Bene, C.; Pennelli, B.; Vuolo, F.; Farina, R.; Pulighe, G.; et al. Capability of Sentinel-2 data for estimating maximum evapotranspiration and irrigation requirements for tomato crop in Central Italy. *Remote Sens. Environ.* **2018**, *215*, 452-470. [CrossRef]

89. Rathod, P.H.; Rossiter, D.G.; Noomen, M.F.; van der Meer, F.D. Proximal Spectral Sensing to Monitor Phytoremediation of Metal-Contaminated Soils. *Int. J. Phytoremediation* **2013**, *15*, 405-426. [CrossRef] [PubMed]

90. Croft, H.; Chen, J.M. Leaf Pigment. *Compr. Remote Sens.* **2018**, *3*, 117-142.

91. Vrieling, A.; Meroni, M.; Darvishzadeh, R.; Skidmore, A.K.; Wang, T.; Zurita-Milla, R.; Oosterbeek, K.; O'Connore, B.; Paganini, M. Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier island. *Remote Sens. Environ.* **2018**, *215*, 517-529. [CrossRef]

92. Mercier, A.; Betbeder, J.; Baudry, J.; Le Roux, V.; Spicher, F.; Lacoux, J.; Roger, D.; Hubert-Moy, L. Evaluation of Sentinel-1 & 2 time series for predicting wheat and rapeseed phenological stages. *ISPRS J. Photogramm. Remote Sens.* **2020**, *163*, 231-256. [CrossRef]

93. Sridhar, B.B.; Han, F.X.; Diehl, S.V.; Monts, D.L.; Su, Y. Spectral reflectance and leaf internal structure changes of barley plants due to phytoextraction of zinc and cadmium. *Int. J. Remote Sens.* **2007**, *28*, 1041-1054. [CrossRef]

94. Wang, F.; Gao, J.; Zha, Y. Hyperspectral sensing of heavy metals in soil and vegetation: Feasibility and challenges. *ISPRS J. Photogramm. Remote Sens.* **2018**, *136*, 73-84. [CrossRef]

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

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).