Hongliang

Pierre Pascal

Alexandre Constantin
email: alexan-dre.constantin@inria.fr

Mathieu Fauvel
email: mathieu.fauvel@inrae.fr

Stéphane Girard
email: stephane.girard@inria.fr

Serge Iovleff

Keywords: Satellite Image Time Series (SITS), Sentinel-2, classification, reconstruction, irregular sampling, Gaussian processes, Earth Observation (EO), remote sensing. Chapter 5 Multivariate Gaussian processes, Classification, Multivariate imputation of missing data, Irregular sampling, Satellite Image Time Series (SITS), remote sensing

Je tiens également à remercier l'équipe Mistis (maintenant Statify) pour son accueil, en commençant par Florence la chef, JB, Sophie, Julyan et Pedro. Merci aussi aux plus anciens doctorants: Alexis, Brice, Thibault, Clément, Karina, Veronica et Fabien. Une mention spéciale aux doctorants avec qui j'ai partagé ces trois années: Meryem pour les randonnées, Masha pour les anniversaires, toutes les petites attentions et d'avoir proposé le café à 10h.. Enfin Benoit, merci pour les discussions toujours intéressantes ainsi que ton aide sur la rigueur mathématique ! Je souhaite maintenant bon courage aux doctorants qui continuent: Dasha (et ses animaux), Théo, Lucrezia, Hana, MinhTri, Louise et Yuchen. Je salue aussi les post-doctorants actuels et anciens de l'équipe: Antoine, Fei,

Résumé

Cette thèse s'inscrit dans le contexte de l'exploitation des données issues de la mission Sentinel-2. Cette mission, initiée par l'Agence spatiale européenne et lancée en 2015, produit massivement des séries temporelles d'images satellite (SITS). Parmi les axes d'analyse de ces images, cette thèse se concentre plus particulièrement sur la classification, c'est-à-dire la production de cartes d'occupation ou d'utilisation des sols en utilisant l'aspect spectrotemporel des SITS issues de Sentinel-2.

Les deux principales difficultés auxquelles ces travaux de thèse se confrontent sont les suivantes. Tout d'abord, la quantité sans précédent de données nécessite à la fois la mise en oeuvre de classifieurs capables de passer à l'échelle et l'utilisation de techniques d'optimisation de code (telles que le traitement parallèle). Deuxièmement, le bruit d'acquisition (nuages, ombres) combiné à l'aspect temporel des données résulte en un échantillonnage irrégulier des séries temporelles. Les approches conventionnelles ré-échantillonnent les séries temporelles sur une grille unique, puis elles utilisent des méthodes d'apprentissage vectorielles pour classer à grande échelle (échelle nationale). Cette démarche en deux étapes a pour inconvénient principal d'alourdir le nombre de traitements appliqués aux SITS, rendant les traitements plus complexes. Dans une moindre mesure, l'étape de ré-échantillonnage peut légèrement altérer les caractéristiques temporelles de la donnée.

Les contributions présentées dans cette thèse sont les suivantes. Nous introduisons une nouvelle approche statistique ayant la capacité de classer des séries temporelles avec un échantillonnage irrégulier basée sur un mélange de processus gaussiens multivariés. Une approche en deux étapes a été proposée, en définissant dans un premier temps un modèle sur des séries temporelles uni-dimensionnelles et indépendantes au sens de l'indépendance spectrale, puis en considérant dans un second temps conjointement les informations spectrales et temporelles des SITS. Ces modèles permettent, de surcroît, une reconstruction de données à des instants non observés ou bruités. L'estimation des deux modèles repose sur un code python parallélisé afin d'être exploitable sur les données de très grande taille. Les deux modèles sont évalués numériquement sur les SITS issues de Sentinel-2 en terme de classification et de reconstruction et sont comparés aux approches conventionnelles. L'analyse des résultats illustre la pertinence des deux modèles et le bénéfice de disposer de modèles paramétriques interprétables.

Mots clés : Processus gaussiens, Analyse de séries temporelles, Optimisation, Données massives, Observation de la Terre.

i

Abstract

This thesis takes place in the context of the processing of the data from Sentinel-2 mission. This mission, initiated by the European Space Agency and launched in 2015, produces an unprecedented amount of Satellite Image Time-Series (SITS). Among the key analyses of these images, this thesis focuses on the classification task, i.e. land use or land cover maps that can be produced using spectro-temporal aspect of the Sentinel-2 SITS.

Two main difficulties are identified in this thesis for the process of Sentinel-2 SITS. First, the unprecedented amount of data requires both scalable classifiers and code optimization techniques (such as parallel processing). Second, the acquisition noise (clouds, shadows) combined with the temporal aspect results in irregular and unevenly sampled time-series. Conventional approaches re-sample time-series to a set of time stamps, then they use machine learning techniques to classify vectors at a large-scale (national scale). The main disadvantage of this two-step processing approach is that it increases the number of operations applied to the SITS, implying a more difficult transition to massive amount of data. To a lower extent, the re-sampling step may slightly alter the temporal characteristics of the data. This thesis contributions are the following. We introduce a novel model-based approach with the ability to classify irregularly sampled time-series based on a mixture of multivariate Gaussian processes. A two-step approach has been used, by defining on one hand a model of uni-variate time-series, independent from the spectral wavelength point of view, then by considering on the second hand both spectral and temporal information from SITS. These models allow jointly a reconstruction of unobserved or noisy data. Estimation of both models has been implemented using a parallelized python code to be scalable to large-scale data-sets. The two models are evaluated numerically on Sentinel-2 SITS in terms of classification and reconstruction accuracy and are compared with conventional approaches. Analyses of the results illustrate the relevance of the two models and the benefit of using interpretable parametric models.

Index terms: Gaussian Processes, Time Series analysis, Optimization, Big Data, Earth Observation.

iii

General Introduction

The context of this PhD thesis is part of a multi-disciplinary research project between remote sensing and applied statistics. This PhD has been done in the Statify team, a joint team between LJK (Jean-Kuntzmann laboratory) and Inria Grenoble and in collaboration with the CESBIO lab, a joint unit with the French national space center (CNES), the French National Centre for Scientific Research (CNRS), the French National Research Institute for Sustainable Development (IRD), Toulouse University (UT2) and the National Research Institute for Agriculture, Food and Environment (INRAe), Toulouse. It has been done under a grant from the CNES and the Grenoble institute of technology (Grenoble INP) in the context of an IRS project (ANR-15-IDEX-02). CNES has also provided data and computational resources.

Finally, in parallel to this thesis, a so-called label has been also obtained with focuses on teaching. It consists of a reflection on the teaching practices and applications on concrete cases in the classroom (experimentation of formats or modes of evaluation for example).

Earth observation using satellite image time-series

With the increasing number of successful launches of orbital satellites, earth observations are more and more used for environmental and climatic monitoring such as land use or land cover. The context of this work focuses on Sentinel-2 constellation, initiated by the European Space Agency with the Copernicus mission. Launched in 2015, it produces an unprecedented amount of time-series thanks to its revisit cycle. These data are used for military and civilian applications. The very high ground resolution allows ecosystems mapping, forests and agricultural plots monitoring. Furthermore, the Sentinel-2 mission covers the complete European area which makes the data usable at a country scale for a country like France. At such spatial and temporal scale, the data contains numerous clouds, shadows, snows (without taking into account perpetual snow) or any other interference with the ground information.

In order to produce a clean map, the time-series at large scale is conventionally done in two steps:

Passive remote sensing satellite images are images taken by an optical sensor known as Multi-or Hyper-Spectral Instrument (MSI/HSI) sensor, onboarded on the satellite's payload. They differ from digital photographies taken by a camera, as aerial ortho-photos, from both mechanical and data point of views. More specifically, an optical sensor measures the reflected light at different wavelengths from the sun by a surface, see Figure 1.1, and divides it by the solar illumination. This quantity is known as reflectance, the reflectance is a function of the wavelength and is a property of the surface from which the light has been reflected, it takes its values between 0 and 1 (see [START_REF]Radiometric quantities : Irradiance, radiance, reflectance[END_REF]). A wavelength, often denoted by λ 0 , is the characteristic length of a wave (from visible spectrum to infra-red in the order of magnitude of a micrometer for optical satellites, see Figure 1.2). The wavelength is inversely proportional to the wave's frequency f 0 : λ 0 = c/ f 0 where c is the speed of light in vacuum.

Physically, the image sensor (MSI/HSI) of a satellite is different from a camera on an electronic point of view [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]Section 3.1.2] as the materials are adapted to space constraints. Indeed, a satellite has a higher orbit than a plane (or helicopter) so an embedded telescope can be found to complete the MSI and additional protections to space constraints (electro-magnetic waves or extreme variations of temperature for example) are also embedded. An other aspect is the reception of images which is more complex for a satellite. A satellite is not able to come back to earth to copy the data. All satellites are completed with a complete array of stations on earth (see [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]Fig. 6] for Sentinel-2 mission) to download images from the satellite. [START_REF] Keshava | Spectral unmixing[END_REF]). Surface reflectance measured thanks to an embedded multi-spectral instrument. 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -9 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -8 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -7 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -6 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -5 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -2 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1 10 -1

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Gamma rays X-rays Infrared 0.4µm 0.8µm

Ultraviolet radio waves

Blue Green Red Satellite images are quantized in three dimensions: spatial, as any image has a finite number of pixels; spectral, as MSI/HSI have a finite number of spectral bands; and temporal, satellite images are taken at different times. The spatial characteristic is denoted by the ground resolution. It is the shape of the pixel at ground level, for the considered satellites, the pixel is usually a square and ground resolution corresponds to its side-length. Then, the main difference between multi-and hyper spectral instruments is the number of spectral bands acquired. Indeed MSI captures between three and a dozen spectral bands whereas HSI captures a hundred of them (higher spectral resolution). The choice of the ground resolution will differ with respect to the targeted application, for example high resolution satellites (as Pléiades family with 50cm ground resolution) are used for civil or military applications [START_REF] Pu | Mapping urban tree species by integrating multi-seasonal high resolution pléiades satellite imagery with airborne LiDAR data[END_REF] and coarser resolution satellites (as Sentinel-3 with 300m ground resolution) with applications to oceanography [START_REF] Donlon | The Global Monitoring for Environment and Security (GMES) Sentinel-3 mission[END_REF]. Finally, the temporal quantization is characterized by the revisit cycle. It is the time between two satellite images acquisitions of the same spatial areas. A complete set of images acquired in a given time window is called a Satellite Image Time-Series (or SITS).

Among the Sentinel missions 1 , Sentinel-2 satellites carry a MSI and are used for environmental monitoring, Land Use or Land Cover (LULC) applications among others. Sentinel-2 SITS are a set of images with the same views (same incidence angle, same altitude, ...). We refer to [13, Table 1] for an overview of the Sentinel family including optical and non-optical (Synthetic Aperture Radar) images. Sentinel-2 satellites have a 10 to 60m ground Figure 1.3: Sentinel-2X tiles, represented by red squares, over South-east Europe (image downloaded from [START_REF]The sentinel-2 tiles, how they work ?[END_REF]). The doubled-lines represent a 10km overlap between two consecutive tiles. resolution sensor including 13 spectral bands and have a revisit cycle of 5 days per satellite.

Sentinel-2 multispectral satellite image time-series are detailed in Section 1.1, then a review of Satellite image time-series classification is presented in Section 1.2. Section 1.3 concludes this chapter.

Sentinel-2 multi-spectral satellite image time-series

This Section presents the SITS produced by Sentinel-2 constellation. It illustrates the time-series aspect and some basic processing done to these data for Land Use or Land Cover (LULC) applications. Section 1.1.1 presents the Sentinel-2 mission in depth. Section 1.1.2 presents the production of the reference data and the problem of mislabeled data in large scale.

Sentinel-2 data

Sentinel-2 [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF], also called Sentinel-2X (or simply S2), is a constellation. Two optical satellites, Sentinel-2A and Sentinel-2B, are currently in orbit (polar orbit, altitude of 786 kilo-meters) and take high resolution images from earth. As an optical satellite, it mainly captures visible light bands and near infra-red spectral bands (from 400 to 2400nm approximatively). The data providers (European Space Agency, ESA and Centre National d'Études Spatiales, CNES for the French territory) provide Sentinel-2 data by tiles. A tile is a 100km side-length square [START_REF]The sentinel-2 tiles, how they work ?[END_REF]. One tile is revisited every five days thanks to low orbital path of Sentinel-2. The tiles over south-east of Europe are presented in Figure 1.3.

Sentinel-2 ground resolution and spectral wavelengths

To provide sufficient information for classification purpose, Sentinel-2 has a wide spectral covering from the visible range to short wave infrared wavelengths with a total of thirteen spectral bands. The choice of a 400 to 2400nm spectral range is the ability to focus on agricultural particularities for different crops [START_REF] Whitcraft | Meeting Earth Observation Requirements for Global Agricultural Monitoring: An Evaluation of the Revisit Capabilities of Current and Planned Moderate Resolution Optical Earth Observing Missions[END_REF]. All spectral bands, their central wavelength and the ground resolution are reported in 1.1: Sentinel-2 bands specifications from [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]. Reported wavelengths are central wavelengths only; bandwidths and Signal to noise ratios are detailed in [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]. Figure 1.4 highlights the wavelengths positions on the spectrum (400 to 2400 nm).

Among the thirteen spectral bands, four wavelengths, including blue, green and red visible light and one infrared band, have a ground resolution of 10m. These bands are sufficient to detect vegetations or non-vegetation areas [START_REF] Tucker | Red and photographic infrared linear combinations for monitoring vegetation[END_REF] as they allow the computation of vegetation indexes as NDVI, the Normalized Difference Vegetation Index. This index is defined as the normalized difference between the infrared (IR) and the red (R) spectral band: NDVI = (IR -R)/(IR + R). Figure 1.4 shows the reflectance w.r.t. to the wavelength for grass and dried grass repeated three times, one for each ground resolution and the spectral filters for all bands are highlighted within the coloured areas. In the third graph of Figure 1.4 (10m ground resolution), it is easy to see that the difference between the reflectance in the grey zone (IR) and the red zone (R) is larger for grass than dried grass. [START_REF] Beck | Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI[END_REF] provides an application to vegetation's monitoring.

Six wavelengths in the near infrared and short-wave infrared have a ground resolution of 20m. Together with the 10m ground resolution spectral bands, they cover the complete range of S2 optical abilities. The infrared spectral bands around 800nm may be used to focus on vegetation and higher values of wavelengths focus on other behaviours.

Finally three wavelengths have a ground resolution of 60m. In general they are used to calibrate atmospheric correction.

Sentinel-2 revisit cycle

As the vegetation is changing with the seasons (for example leafs are falling from deciduous trees in Autumn), the previous spectra (Grass or Dried Grass) depend on the acquisition day in a year: the spectral signature may also change with time [START_REF] Boori | Monitoring crop phenology using ndvi time series from sentinel 2 satellite data[END_REF]. The following focuses on the temporal aspect thanks to the high revisit cycle of S2.

Despite the intrisinc noise of electronic components which is not of interest here, the noise considered in the following concerns environmental effects which makes the acquisition unusable. As an optical satellite, the noise sources we have is the presence of clouds and shadows. Presented in [START_REF] Whitcraft | A Framework for Defining Spatially Explicit Earth Observation Requirements for a Global Agricultural Monitoring Initiative (GEOGLAM)[END_REF] and in [START_REF] Whitcraft | Meeting Earth Observation Requirements for Global Agricultural Monitoring: An Evaluation of the Revisit Capabilities of Current and Planned Moderate Resolution Optical Earth Observing Missions[END_REF], an average revisit cycle with less than 8 days allows a 70% of noise-free data, also called clear data. Thus, about 30% of the data are corrupted by clouds or shadows. The Figure 1.5 represents an extract of one SITS (true color composition) from France for a complete year of satellite images acquisition. We can clearly see variations of colors (which means a variation of the reflectance with respect to the time) with seasonal effects, clouds or shadows at some locations. An other tile from the city center of Toulouse, France, is represented in Figure 5.1. In the following, T = [t min , t max ] denotes the time window for a SITS: t min is close to January the 1st and t max around December the 31st for a one-year SITS.

Sentinel-2 level 2A SITS and noise detection

Now that the Sentinel-2 SITS are described, what follows describes the final product, i.e. level 2A SITS reflectance, and the construction of the mask of data (noise detection). 1.1 and with associated bandwidths [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]. Visible light starts from 380nm and ends around 780nm: blue, green and red bands are highlighted in the 10m ground resolution bands. The spectral reflectance signature of healthy grass is reported on each grid (full line) and the signature of dried grass (dashed line). [START_REF] Hagolle | The product level names, how they work ?[END_REF]). Sentinel-2 product levels: level 1C (with atmospheric interference and clouds) and level 2A (reflectance and clouds removal). Level 3A shows reconstructed satellite image.

The SITS are proposed by the data provider at different level of processing ranging from the raw data to clean an reconstructed surface reflectance. These levels are called product levels. Product levels available for Sentinel-2 are the ones derived from the NASA (National Aeronautics and Space Administration) in 1996, some of the product levels are presented in Figure 1.6. The multi-spectral instrument from the satellite captures a raw reflectance (level 0) which is completed with metadata (level 1A: acquisition time, geo-references, ...) and then ortho-rectified (level 1C). At levels 0 and 1, reflectances are still interfered by the atmosphere and have to be corrected to be used as surface reflectances (see Figure 1.1).

Image products are downloaded2 at level 2A where the data are composed by surface reflectance and a mask [START_REF] Hagolle | A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENµS, LANDSAT and SENTINEL-2 images[END_REF]. At time t ∈ T , we denote by Y(t) ∈ E the surface reflectance at multiple spectral bands and M(t) ∈ R + the associated mask. M(t) > 0 means that clouds or shadows are detected. Figure 1.7 shows a true color satellite image (left image) and the same image where noisy data are masked (right image) thanks to the associated mask. All product's levels for Sentinel-2 are presented and discussed in [START_REF] Hagolle | The product level names, how they work ?[END_REF]. These pre-processing steps are based on the MACCS (Multi-sensor Atmospheric Correction and Cloud Screening) / MAJA (MACCS-ATCOR Joint Algorithm) a processing chain [START_REF]MACCS/MAJA, how it works[END_REF] from CNES (French space agency), CESBIO (National center for space and biosphere studies) and the German Aerospace Center (DLR, [START_REF] Baetens | Validation of Copernicus Sentinel-2 Cloud Masks Obtained from MAJA, Sen2Cor, and FMask processors using reference cloud masks generated with a supervised active learning procedure[END_REF]).

Finally, some data are in the edge of a tile or, for some reason, the data are not exploitable. In that case there is no data. The data provider indicates it using "reflectance" value -10000. No data "values" are automatically removed from the time-series.

Land use or land cover maps

The construction of LULC maps, especially at large scale (country scale or larger), is often computationally expensive and the reference data are partially mislabelled. Mainly two reasons are identified, firstly the labels change over the years (constructions of urban areas, clear-cutting, etc) and yields to a heavy and expensive update at large scale. Secondly the change of ground resolution between two satellites missions imply a re-labelling of reference data. For example Formosat SITS are well-known for LULC maps and have now corrected the reference data throughout their use but has a different ground resolution to be used directly with Sentinel-2 SITS. The following focuses on the construction of the sets of Sentinel-2 reference data.

SITS reference data

The classes in this manuscript are extracted from [93, Section 4.2] and grouped by polygons. Each polygon refers to a unique class designated by CORINE Land Cover (CLC) dataset [START_REF] Bossard | CORINE land cover technical guide: Addendum 2000[END_REF], French National Geographic Institute (IGN) maps [START_REF] Maugeais | Dérivation cartographique multi échelles de la BDTopo de l'IGN france: Mise en oeuvre du processus de production de la nouvelle carte de base[END_REF] and other more specialized providers (grasslands and glaciers). The classification problem aims to separate around twenty classes. In a tile-based approach, the number of classes to predict may change from one tile to another, for example Glaciers polygons only appear in high altitudes. An example of resulting polygons is reported in Figure 1.8 (Other classes are presented in Figure 4.7).

As the authors mentioned in [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF], they took into account multiple reference datasets to obtain large-scale maps but the proportion of mislabeled data is still significant. To overcome this issue, the authors used a Random Forests classifier (discussed later) known as robust for mislabeled data.

The diversity of land use or land cover maps

The production of Land-Use or Land-Cover (LULC) maps in the literature is wide. Depending on the context, the data source, the geographical area and the land cover may vary.

From [START_REF] Chen | Global land cover mapping at 30m resolution: A POK-based operational approach[END_REF][START_REF] Gong | Finer resolution observation and monitoring of global land cover: First mapping results with Landsat TM and ETM+ data[END_REF][START_REF] Yu | Improving 30 m global land-cover map FROM-GLC with time series MODIS and auxiliary data sets: A segmentation-based approach[END_REF], global land cover maps of high resolution (Landsat 30m ground resolution images) are discussed and reached classification scores from 64% to 70% for a total of 9 to 10 classes including agricultural areas or artificial areas among others. Only recently, [START_REF] Karra | Global land use / land cover with sentinel 2 and deep learning[END_REF] presented a LULC map at 10m using Sentinel-2 images with 85% classification score for 11 classes. The LULC has also been studied for smaller maps as continents or countries. [START_REF] Giri | Land Cover Characterization and Mapping of South America for the Year 2010 Using Landsat 30 m Satellite Data[END_REF] presented the case of South America with 5 classes with an overall accuracy of 89%; [START_REF] Homer | Conterminous United States land cover change patterns 2001-2016 from the 2016 National Land Cover Database[END_REF] for United State of America from 2001 to 2016 with 16 classes and 80% to 83% of estimated accuracy. France Metropolitan area has also been mapped [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF] with 17 classes and about 90% accuracy score using Sentinel-2 10m ground resolution. China was mapped in [START_REF] Li | Using Land Long-Term Data Records to Map Land Cover Changes in China Over 1981-2010[END_REF] with 19 classes and about 75% accuracy score. [START_REF] Constantin | Joint supervised classification and reconstruction of irregularly sampled satellite image times series[END_REF]). Polygons classes over a true colors tile in the center of the metropolitan France. Red squares are zoomed areas presented in [START_REF] Constantin | Joint supervised classification and reconstruction of irregularly sampled satellite image times series[END_REF].

The definition of the LULC problem will not be discussed further. Next Section only focuses on the classification techniques and comparison for a given LULC problem. This part highlights that a direct comparison between different studies is difficult.

Review of satellite image time-series classification

In the following E denotes the space where the SITS objects y live, i.e. y ∈ E, and C is the total number of classes in the classification problem. In the following we consider non-noisy data, y can be represented by a vector of fixed length, and we consider supervised classification problems, each vector y is associated with a class membership. A common example of a vector input space is E = R pq where pq is the dimension of the vector. q represents the number of temporal features within SITS and p the number of spectral bands. Other cases of E are considered in the next Chapter when reviewing state of the art in statistical modelling. The supervised framework implies that the models' parameters are learned from a training set including the class membership. The problem of supervised classification is recalled in Section 2.1.

Production of land use or land cover maps

As presented previously, SITS are complex for multiple reasons: they have spatio-spectro temporal information, they are noisy (optical satellite images) and they are produced in very large scale. The processing and classification tasks are challenging as the use of some classifiers is not possible at large scale or some additional assumptions must be done such as spatial or temporal independence, etc. A recent review can be found in [START_REF] Talukdar | Land-Use Land-Cover Classification by Machine Learning Classifiers for Satellite Observations-a Review[END_REF] for different LULC purpose and scales. The "conventional" classifiers for LULC applications are presented, they were intensively used in the past decade and were issued from the Machine Learning community. Then more complex classifiers are presented, they are gaining in interest and are part of the deep neural network family.

Conventional classifiers

Among the conventional classifiers, this chapter will not review model-based classifiers as Gaussian Mixture models [START_REF] Lagrange | Large-Scale Feature Selection With Gaussian Mixture Models for the Classification of High Dimensional Remote Sensing Images[END_REF]. These models are described in-depth in the Chapter 2. Two major classifiers are reviewed: Random The class membership associated with y is given by the red path (c = C -1) in this example.

Forests (RF) and Support Vector Machines (SVM). Random forests (RF, [START_REF] Breiman | Random Forests[END_REF]) are probably the most known classifier in the remote sensing literature. RF is a tree-based methods, it is a collection of numerous decision trees. A decision tree is a binary tree which can take as input a finite-dimensional vector, y ∈ E. At each node, a decision rule gives a direction to the decision (left or right) and results into a set of labelled nodes by c ∈ {1, . . . , C}, a schematic view is reported in Figure 1.9. The data is classified thanks to its final node, for example the red path in Figure 1.9 leads to class c = C -1. We refer to [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 9] for additional properties. Each decision tree within the RF classifier returns a class, the class membership returned by RF is the majority voted class from decision trees.

Considering a time-series observed at q time-stamps with p spectral bands and a spatial independence (as done in the most cases [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF]), then the input space E is often a vector of size pq and E = R pq where all the p time-series are stacked. In [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF], the authors presented promising results on Landsat 8 spectra-temporal SITS before the launch of Sentinel-2. They trained multiple RF classifiers, one per tile over the France territory and one per climatic area. They used linear interpolation to re-sample time-series. Spatio-temporal information can also be exploited, in [START_REF] Gbodjo | Toward Spatio-Spectral Analysis of Sentinel-2 Time Series Data for Land Cover Mapping[END_REF] they have shown good results on Sentinel-2 at regions level. Finally [START_REF] Verhegghen | Accuracy assessment of the first EU-wide crop type map with LUCAS data[END_REF] used RF to produce crops maps over Europe on Sentinel-1 time-series.

In addition to the good classification results, [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF] highlighted that RF is robust when some target classes within the training set are wrong. RF is known to scale well w.r.t. the number of samples.

Support Vector Machines (SVM, [START_REF] Cortes | Support-vector networks[END_REF]) were investigated for LULC applications before RF [START_REF] Sheykhmousa | Support Vector Machine Versus Random Forest for Remote Sensing Image Classification: A Meta-Analysis and Systematic Review[END_REF]. In binary classification (C = 2), it aims to find an optimal hyperplane to separate the two classes, see Figure 1.10. The hyperplane is optimal if it maximizes the margin between the two classes and minimizes the classification error in the training set. From there, two versions of SVM are found in the literature, linear SVM with an optimum hyperplane within the data space and non-linear SVM which transforms the data space into a feature space where the data are more likely linearly separable.

In linear SVM, the hyperplane is linear w.r.t. the input y ∈ E, it corresponds to:

ω, y E + b = 0, (1.1)
where {ω, b} are parameters of the model. Non-linear SVM can be defined by replacing the dot product in (1.1) by a kernel function (kernel trick), we refer to Section 2.2.1 for more details on kernels or to [START_REF] Gómez-Chova | A Review of Kernel Methods in Remote Sensing Data Analysis[END_REF] for kernel methods in remote sensing. However, despite higher abilities to separate classes (polynomial or RBF hyperplanes among others, see Section 2.2.1), it is harder to scale to large data-sets. Indeed, only linear SVM are considered as "light".

Finally the problem is extended to multi-class (C > 2) using the so-called one-vs-all or one-vs-one strategies. The latter one is often preferred [START_REF] Hsu | A comparison of methods for multiclass support vector machines[END_REF]. Multiple comparisons between SVM and RF have been done on LULC applications. [START_REF] Huang | An assessment of support vector machines for land cover classification[END_REF] compared SVM with decision trees and two other classifiers (model-based and Neural network based approaches) on a 6 classe problem. [START_REF] Noi | Comparison of Random Forest, k-Nearest Neighbor, and Support Vector Machine Classifiers for Land Cover Classification Using Sentinel-2 Imagery[END_REF] provided a more recent comparison between SVM, RF and k-Nearest Neighbors classifiers on Sentinel-2 images. However, despite good results, the choice of non-linear SVM with more complex kernels Figure 1.10: (Image from [START_REF] Sheykhmousa | Support Vector Machine Versus Random Forest for Remote Sensing Image Classification: A Meta-Analysis and Systematic Review[END_REF]). Two classes (red and green) are separated thanks to the optimal hyperplane (blue line), i.e. the hyperplane which maximizes the margin between the two groups. In this example the data y are real-input vectors.

scales quadraticly or cubicly with respect to the number of samples within the training set an large data-set remains a problem.

Applications in LULC are also wide, in addition to the above comparisons: [START_REF] Mathur | Crop classification by support vector machine with intelligently selected training data for an operational application[END_REF] used SVM on spectral features for crop classification and [START_REF] Camps-Valls | Kernel-Based Framework for Multitemporal and Multisource Remote Sensing Data Classification and Change Detection[END_REF] used time-series for data change detection. The use of both SVM and RF techniques is still very important. As presented in [START_REF] Sheykhmousa | Support Vector Machine Versus Random Forest for Remote Sensing Image Classification: A Meta-Analysis and Systematic Review[END_REF], RF is becoming more and more important than SVM in the remote sensing literature over the last four years. Both methods are still of interest as the number of parameters is reasonable to reach high accuracy on the produced map and the interpretability of some of the model parameters can be done. The following focuses on more complex techniques with a larger number of trainable parameters.

Artificial neural network classifiers

Artificial deep neural networks are widespread in LULC applications on smaller area as they have millions of parameters to learn. The number of parameters induces naturally a huge amount of computation time and resources.

A Neural Network (NN) is a structure based on the multi-layer perceptron (MLP) from Rosenblatt [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF]. A perceptron is a parametric function with parameters {ω, b}, ω ∈ E and b ∈ R that has been transformed through the Heaviside step function h : R → {0, 1} as in (1.2). A perceptron is also known as an artificial neuron as the output is activated (Output of h equal to 1) or inactivated otherwise:

z = h ( ω, y E + b) .
(1.

2)

The MLP has also been studied in remote sensing as in [START_REF] Huang | An assessment of support vector machines for land cover classification[END_REF] which compared RF/SVM with the MLP structure. However, the advent of new kind of layers within neural networks (derived from MLP) and the big data era brought the intensive study of Deep Neural Network (DNN) with a large number of intermediate layers, also called hidden layers, over the past ten years. Every layer, including hidden layers, are composed by weights with different shapes and each layer is terminated by a non-linear function (more general than heaviside step function, like relu, elu, etc). The dimensions of input, output and hidden layers and the parameters (or weights) assigned to each layer with the associated activation function, together, is called a network architecture. A DNN architecture with a large number of hidden layers easily reaches millions of parameters. Inference of DNN is complex [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Section 11.5] and done by back propagation of the gradient, issues related to a large number of parameters arise but also overfitting problems which have to be taken into account within the training of such network (dropout techniques, data augmentation, etc).

In the context of remote sensing, mainly two architecture's families of DNN are studied in the literature: the first family is composed by convolutional layers and is known as Convolutional NN (CNN), they are specialized in Figure 1.11: (Image from [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF]). U-net architecture for pixel-wise classification. image analysis (or spatial features). The second family is composed by recurrent layers and is known as Recurrent NN (RNN). They are specialized in time-series (temporal features).

A Convolution layer is a convolution product between a so-called convolutional filter and a tensor input. From there two structures are derived from the CNN, the first architecture takes as input an image (real-input matrix, or higher order tensor) and has C outputs (real-valued vector) where C is the number of classes in the classification problem. The number of parameters is correlated with the size of the convolutional filter and the complexity decreases significantly compared to a fully-connected network (for all hidden layers, each output is a linear combination of all inputs). The convolutions are mostly used for application to images. In remote sensing they are used to exploit spatial information of satellite images. For example [START_REF] Syrris | Evaluation of the Potential of Convolutional Neural Networks and Random Forests for Multi-Class Segmentation of Sentinel-2 Imagery[END_REF] uses CNN to produce maps from Sentinel-2 images by upsampling the final layer to the size of the image. An other architecture derived from CNN is the autoencoder architecture. The auto-encoder shares some dimensions of the output with the input of the network. The most known architecture is the U-net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF], it has the spatial dimensions shared between inputs and outputs. U-net is used for pixel-wise classification as illustrated in Figure 1.11 (outputs are tensors of order 3, the last dimension can be increased for different classification tasks. An application to aerial images is presented in [START_REF] Kim | Objects Segmentation From High-Resolution Aerial Images Using u-Net With Pyramid Pooling Layers[END_REF].

The second family, Recurrent Neural Network (RNN), uses the temporal features of SITS. This family contains LSTM (Long-Short Term Memory, [START_REF] Hochreiter | Long Short-Term Memory[END_REF]) networks which are interesting for multiple views of the same scene at different times to predict, in this context, pixel-wise class membership. It has been used in [START_REF] Qiu | Local climate zone-based urban land cover classification from multi-seasonal Sentinel-2 images with a recurrent residual network[END_REF][START_REF] Rußwurm | Multi-temporal land cover classification with long short-term memory neural networks[END_REF] on Sentinel-2 images and shows promising results.

Some recent works also focus on combining recurrent networks with convolutions for images [START_REF] Garnot | Satellite image time series classification with pixel-set encoders and temporal self-attention[END_REF][START_REF] Pelletier | Temporal Convolutional Neural Network for the Classification of Satellite Image Time Series[END_REF] on Sentinel-2 images. For two in-depth reviews, we refer to [START_REF] Zang | Land-Use Mapping for High-Spatial Resolution Remote Sensing Image Via Deep Learning: A Review[END_REF] for deep NN models applied to high spatial resolution satellite images and to [5] for application to hyper-spectral images. At this time, [START_REF] Karra | Global land use / land cover with sentinel 2 and deep learning[END_REF] presents a large scale map with 10m ground resolution but without temporal information from SITS. When dealing with large data set at a country scale with large time domain (one year of images) at high resolution, the conventional machine learning techniques are preferred when DNN architecture reached millions and millions of parameters to train.

Conventional and neural network classifiers, take as input a vector of fixed length, a situation that can be hardly reached with SITS. The next Section discusses on problems induced by missing data.

Classification with missing data

When dealing with missing data, the previous classifiers are not able to perform without re-sampling to a finite and fixed grid of time stamps. This grid must be the same within the training and the validation sets of data. For Sentinel-2 spatio-spectro-temporal images, the reconstruction of missing data to a fixed temporal grid may differ for different studies. We refer to [START_REF] Shen | Missing Information Reconstruction of Remote Sensing Data: A Technical Review[END_REF] for a review of these techniques and how to combine them to reconstruct missing data. Before going further, let us recall that S2 SITS are processed at level 2A (Section 1.1.1) and are brought with a mask with noise detection. As the 60m ground resolution are used to process the data from 1A to 2A, the 60m spectral bands will not be used in this manuscript.

Spatial re-sampling

Let us recall that the spatial location of each pixel has not been taken into account in this manuscript, the timeseries will be classified thanks to their spectral and temporal signature (using non-noisy time stamps). We refer to [START_REF] Shen | Missing Information Reconstruction of Remote Sensing Data: A Technical Review[END_REF][START_REF] Zhang | Missing Data Reconstruction in Remote Sensing Image With a Unified Spatial-Temporal-Spectral Deep Convolutional Neural Network[END_REF] for missing data reconstruction using spatial information.

The process will only concern 20m ground resolution spectral bands (Table 1.1). A simple up-sampling from 20m to 10m ground resolution will be done. It uses the "bicubic" algorithm from orfeo toolbox [171], in particular superimpose; this toolbox is an open-source project for remote-sensing applications.

Temporal re-sampling

Sentinel-2 SITS have different temporal samplings for each time series. The reasons for this two-fold: firstly the Sentinel-2 tiles have different time stamps from one tile to another (provided by the data provider) thanks to the orbital path of the satellites, an example is provided in Figure 1.12 for three different tiles by the blue and red dots. The red dots correspond to masked data (noisy acquisitions). Secondly the noise sources as clouds or shadows appear locally, i.e. the clear data (or noise-free data) from one time-series to another appear at different times in a year. To overcome these issues, [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF] defined a unique temporal grid over a complete territory (metropolitan France in their study). The authors removed noisy acquisition thanks to the mask values (red dots in Figure 1.12) and linearly re-sampled time-series to the new fixed grid. This technique, simply mentioned as gap-filling, is implemented in orfeo toolbox [171]. It consists in taking the linear interpolation between two clear dates, an illustration is provided in Figure 1.12 where dashed lines represent the gap-filled time-series.

Other techniques are available but are not used to process Sentinel-2 SITS. In Chapter 5 we review statistical methods for missing data reconstruction.

Challenges of satellite image time-series classification

This chapter presented the Satellite image time-series and the Sentinel-2 constellation which operates to produce high resolution images. We finished the presentation of Sentinel-2 SITS with the emphasis on mislabelled data sets, particularly in large scale areas. These data will be used throughout this manuscript in order to assess the performances of the proposed classifiers thanks to the spectro-temporal aspect of the data (a spatial independence assumption will be done). We also presented the state-of-the-art classifiers with finite-sized input vectors, with promising classification scores and working on national scale to produce Land Use or Land Cover (LULC) maps. Application to Sentinel-2 SITS requires an additional pre-processing step to re-sample time-series to a fixed set of time stamps. We illustrated this processing step on one time-series using the linear gap-filling technique which is usually done to process national LULC map [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF].

The main challenges are the classification of these SITS using classifiers that scales to large data-sets and to remove the processing step of re-sampling the time-series to a known set of time-stamps.

Recent statistical approaches are of interest thanks to their interpretability (unlike more complex classifiers as Neural Networks) and the wide variety of applications. An example is the use of the Gaussian processes [START_REF] Camps-Valls | A Survey on Gaussian Processes for Earth-Observation Data Analysis: A Comprehensive Investigation[END_REF] for the analysis of the three dimensions of SITS

French introduction

Ce chapitre présente la classification de séries temporelles par modélisation statistique. La modélisation statistique peut être séparée en deux grande catégories : l'apprentissage non supervisé et l'apprentissage supervisé. Ce chapitre se concentre plus particulièrement sur ce dernier.

La première partie présente la classification supervisée. Pour cela la règle de décision (assignation d'un élément à une classe) est présentée pour différents types d'objets à classer : des vecteurs, des matrices ou bien des objets plus complexes (données textuelles, etc). Différents modèles génératifs sont alors présentés, en particulier les modèles supposant une distribution gaussienne. Plus particulièrement la loi gaussienne multivariée est étudiée puis étendue au cas matriciel.

La seconde partie présente la statistique fonctionnelle en détaillant les processus gaussiens (Gaussian Processes, GP). Deux cas sont abordés : l'utilisation des processus gaussiens dans le cadre de problèmes de régression et de classification. La régression repose sur le principe que les échantillons sont des observations d'une fonction unique puis, par la loi des probabilités conditionnelles pour des vecteurs gaussiens, permet la prédiction de valeurs à des instants non observés. La classification repose sur l'hypothèse que la probabilité d'appartenance à une classe est composée d'une fonction à valeurs dans [0, 1] et d'une fonction latente modélisée par un GP. Le cas de la classification binaire est présenté puis étendu au cas multi-classes. Enfin, cette partie se termine par un bref état de l'art sur les GP multi-dimensionnels, i.e. des fonctions aléatoires à valeurs dans R p avec p > 1. Ces méthodes de classification ont, pour l'essentiel, une complexité cubique par rapport au nombre d'échantillons. Ces notions sont nécessaires pour la compréhension des contributions.

Ces méthodes ne permettent pas, à notre connaissance, de classer les séries temporelles à échantillonnage irrégulier, notamment les SITS Sentinel-2. Ces méthodes reposent sur des objets de dimensions finie et fixe, elles impliquent un ré-échantillonnage temporel des SITS, ce qui n'est pas compatible avec notre problématique.

Let E be an arbitrary space containing all objects to classify (vectors of dimension p, time-series, etc) and let {1, . . . , C} be a set of C classes. Let S = {(y i , z i )} n i=1 be a set of n independent and identically distributed (i.i.d.) samples, or realizations, from a random pair (Y, z) ∈ E × {1, . . . , C}. For any given class c ∈ {1, . . . , C}, the model-based classification problem is to assume that the conditional distribution of y|z = c, y ∈ E, belongs to some parametric family. For any sample i, y i |z i = c means that y i belongs to the class c with associated density p(y i |z i = c). This approach differs from the previous approaches where the data mechanism is not modelled by a statistical distribution [START_REF]Statistical Modeling: The Two Cultures[END_REF] as presented in Chapter 1. The statistical approaches are interesting for their interpretability and, often, reduced numerical complexity thanks to a reduced number of model parameters. In this chapter, we focus on statistical modelling by describing supervised model-based classification at first in Section 2.1 with a focus on Gaussian distributions. Afterwards an overview of Gaussian processes for regression and classification is given in Section 2.2. Finally a brief conclusion is provided in Section 2.3.

Supervised model-based classification

Model-based classification is based on the definition of a decision rule δ which returns a class c considered as the best candidate from an observation y ∈ E. We say that ĉ is the class membership of y:

δ : E → {1, . . . , C}, y → ĉ.
In a supervised framework, all classes are known and represented in the set of data. The probabilistic model is inferred on the complete set {(y i , z i )} n i=1 . On the opposite, an unsupervised framework supposes that the classes are unknown and are inferred from the data inputs {y i } n i=1 . Firstly the decision rule is presented Section 2.1.1, then the discrimination problem with a Gaussian assumption is presented in Section 2.1.2 and finally extensions to matrix-variate Gaussian distributions are presented in Section 2.1.3.

Decision theory and supervised framework

Maximum a posteriori rule

In the Bayesian decision theory, the optimal decision rule is called the maximum a posteriori (MAP) rule. The Bayes rule gives us the a posteriori probability by:

P(z = c|y) = π c p(y|z = c) p(y) ∝ π c p(y|z = c),
where p(y|z = c) is the likelihood density and π c = P(z = c) the a priori probability on class c. p(y) is the marginal distribution of y which can be written as a finite mixture:

p(y) = C c=1 π c p(y|Z = c).
The MAP rule assigns to a non-labeled object the class ĉ such as:

δ(y) = ĉ = arg max c P(z = c|y) = arg max c π c p(y|z = c). (2.1)
From (2.1), two approaches are found in the literature: the discriminative approach models the a posteriori probability and the generative approach which models the a priori probability and the likelihood. The latter approach is the one considered in the following.

In a supervised framework, π c is usually estimated by its empirical counterpart and the generative approach is to model the likelihood p(y|z = c). A model is a distribution which can be parametrized using a finite number of parameters denoted by θ. If we suppose that the problem is decoupled w.r.t. each class, we define C models where each model in class c ∈ {1, . . . , C} is parametrized by its own parameters θ c . We write the parametrized likelihood density in class c by p(y|z = c; θ c ). In the following we introduce the subset S c of S such that S c = {(y i , z i )} i|z i =c . We write n c = |S c |, then n c ≤ n and C c=1 n c = n. Remark. Let us highlight that the optimization problem (2.1) is equivalent to:

ĉ = arg min c k log (π c p(y|z = c; θ c )) ,
where k is a non-zero negative constant (k = -2 for Gaussian distributions as log is strictly increasing on R + .

Maximum likelihood estimation

Maximum likelihood estimator (or MLE) is an estimator which maximizes the likelihood function w.r.t. parameters θ. Optimal parameters are denoted by θMLE and defined by: θMLE = arg max θ L(θ, y).

(2.2)

L(θ, y) = C c=1 p(y|z = c, θ c
) is the complete likelihood function of the model with parameters θ = {θ 1 , . . . , θ C } on data y ∈ E. As all samples from S are independent and identically distributed, the complete likelihood is the product of density p(y i |z i = c; θ c ), also called marginals. Additionally, when the problem is decoupled w.r.t. the class, then (2.2) is computed for each class, i.e. for each sample in subset S c : θc,MLE = arg max

θ c i|z i =c p(y = y i |z i = c; θ c ).
(2.3)

For the sake of clarity, the realization y = y i given z i = c will be directly written

y i |z i = c. Generally,

Probabilistic models and distributions

As mentioned previously, the data y live in an arbitrary space E. In the following we discuss different examples of space E and associated distributions found in the literature.

Let q be the number of features of y. If E = R q then the samples are real-valued vectors, for example regularly sampled time-series. In this case, the multivariate (q-variate) Gaussian distribution is often adopted, leading to the well-known linear and quadratic discriminant analysis classifiers (see Section 2.1.2). The probability density function (pdf) of the multivariate normal distribution is given in (2.4).

Recent works focus on non-Gaussian distributions. One example is the skew-normal distribution which has been introduced in the past century [START_REF] Azzalini | A class of distributions which includes the normal ones[END_REF]. Considering the one-dimensional case (E = R), the pdf of the skew-normal distribution, denoted by p s-n , is defined with a parameter γ which controls the skewness (γ = 0 yields the standard Gaussian distribution):

p s-n (y) = 2p(y)Φ(γy), where p(.) denotes the standard normal pdf and Φ(.) is the standard Gaussian cumulative density function (Φ(y) = y -∞ p(t)dt). This distribution is used to deal with asymmetric data [START_REF] Chamroukhi | Skew t mixture of experts[END_REF][START_REF] Theodossiou | Financial Data and the Skewed Generalized t Distribution[END_REF], the multivariate skew-normal distribution is described in [START_REF] Azzalini | The Multivariate Skew-Normal Distribution[END_REF].

An other example is the t-distribution with heavier tails. The distribution depends on a shape parameter ζ (ζ → ∞ yields the standard Gaussian distribution), the multivariate t-distribution pdf is given by:

p(y) = Γ ((ζ + q)/2) Γ(ζ/2)(ζπ) q/2 |Σ| 1/2 1 + 1 ζ (y -µ) Σ -1 (y -µ) -(ζ+q)/2 , y ∈ R q ,
where Γ is the gamma function and |Σ| is the determinant of Σ. This distribution is used to deal with outliers [4,[START_REF] Murray | A mixture of SDB skew-t factor analyzers[END_REF]. The distributions are illustrated in Figure 2.1. It presents a right skewed pdf (γ > 0) and a t-distribution with heavy tails (ζ = 1) with comparison to the standard Gaussian pdf. More details about these two parameters are provided in [START_REF] Bouveyron | Model-Based Clustering and Classification for Data Science: With Applications in R[END_REF]Chapter 9]. When y are categorical data, for example nominal data as species analysis [START_REF] Elith | Species distribution models: Ecological explanation and prediction across space and time[END_REF] (presence or absence of species in an environment), the space E is discrete. In that context, some studies focus on multinomial distributions [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF] or Dirichlet distributions [START_REF] Bouguila | Novel Mixtures Based on the Dirichlet Distribution: Application to Data and Image Classification[END_REF]. If E is ordered, then dedicated distributions are proposed as in [START_REF] Biernacki | Model-based clustering of multivariate ordinal data relying on a stochastic binary search algorithm[END_REF] with ordinal data. Other data-sets, or spaces E, are presented in [START_REF] Bouveyron | Model-Based Clustering and Classification for Data Science: With Applications in R[END_REF]Chapter 6]. Lastly, when E is finite and high dimensional or infinite dimensional (i.e. when E = R q and q → ∞), specific parametric and non-parametric models are preferably used. Section 2.1.2 presents the High Dimensional Discriminant Analysis (HDDA) and Section 2.2.1 discusses the kernel functions that can be used in non-linear classifiers for infinite dimensional data.

General discriminative problem

As introduced before, when E = R q we can use linear and quadratic classifiers. The following presents the Linear and Quadratic Discriminant Analysis (known as LDA and QDA) with Gaussian assumption on the likelihood density. Linear and Quadratic Discriminant Analysis (LDA/QDA) Let y|z = c ∼ N q (µ c , Σ c ) be a multivariate Gaussian model with mean µ c and covariance matrix Σ c . The model parameters are θ c = {µ c , Σ c }. Then the associated q-variate Gaussian likelihood density (or pdf) for the i-th sample, or realization from the Gaussian, y i |z i = c is given by:

p(y i |z i = c; θ c ) = (2π) -q/2 |Σ c | -1/2 exp - 1 2 (y i -µ c ) Σ -1 c (y i -µ c ) . (2.4)
From (2.3) in supervised framework (recall that a-priori probabilities are estimated by their empirical counterpart) with Gaussian density (2.4), the minimization of the negative log-likelihood is:

ĉ = arg min c i|z i =c (y i ; θ c ).
with:

(y i ; θ c ) = -log(p(y i ; θ c )) ∝ D 2 Σ c (y i , µ c ) + log |Σ c | + q log(2π). D Σ (y, µ) = (y -µ) Σ -1 (y -µ)
stands for the Mahalanobis distance with covariance Σ. The optimization problem minimizes both the quadratic term (D Σ (y, µ)) and the complexity (|Σ c |). Then the MAP rule δ(y i ) = arg max c P(z i = c|y i ) returns the class with the highest probability, this method is known as the Quadratic Discriminant Analysis (or QDA).

If all covariance matrices Σ c , c ∈ {1, . . . , C}, are assumed to be equal to Σ, i.e. θ c = {µ c , Σ}, then the covariance matrix is shared for all classes. This particular case is known as the Linear Discriminant Analysis (LDA) as explained below.

Decision boundaries

A decision boundary describes a delimitation between groups where the probabilities to belong to all classes are equal. Let us consider the case C = 2 with E = R q , then the decision boundary delineates two groups. Let s(y) be the quantity defined by: If the probabilities are equal, then s(y) = 0, ∀ y ∈ E = R q . In particular, the space {y ∈ R q |s(y) = 0} defines the decision boundary.

s(y) = log P(z = 2|y) P(z = 1|y) = log(p(y|z = 2)) -log(p(y|z = 1)), ∀ y ∈ E. (2.5)
If the distribution is Gaussian, s(y) = 0 is equivalent to:

D 2 Σ 2 (y, µ 2 ) + log |Σ 2 | -D 2 Σ 1 (y, µ 1 ) + log |Σ 1 | = 0, ⇔ (y -µ 2 ) Σ -1 2 (y -µ 2 ) -(y -µ 1 ) Σ -1 1 (y -µ 1 ) + log |Σ 2 | |Σ 1 | = 0, ⇔ y Σ -1 2 y -y Σ -1 1 y -2µ 2 Σ -1 2 y + 2µ 1 Σ -1 1 y + µ 2 Σ -1 2 µ 2 -µ 1 Σ -1 1 µ 1 + log |Σ 2 | |Σ 1 | = 0.
(2.6)

The boundary equation in (2.6) is quadratic w.r.t. y, explaining why it refers to the Quadratic discriminant analysis.

Taking

Σ 1 = Σ 2 = Σ in (2.6) yields: 2(µ 1 -µ 2 ) Σ -1 y + µ 2 Σ -1 µ 2 -µ 1 Σ -1 µ 1 = 0,
which is linear w.r.t. y, referring to the linear discriminant analysis. Figure 2.2 illustrates the linear and quadratic decision boundaries between two classes (blue and red dots) for samples y i ∈ E = R 2 , ∀i ∈ {1, . . . , n}.

Estimation of parameters

Recall that estimation of parameters is done by maximum likelihood in each class c independently. Inference in class c is done on the set S c with cardinal n c , the number of samples in class c. θc,MLE is solution of (2.3) and πc = n c /n where n = c n c . In QDA, the optimal parameters are well known [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 4] and are given by: μc

= i|z i =c y i /n c , Σc = i|z i =c (y i -μc ) (y i -μc )/n c .

Penalized and regularized discriminant analysis

In the Gaussian discriminative analysis, some additional hypothesis are usually made to give the model more robustness to real world data-sets especially in high dimension. Most of them consist of modifying the estimator of the covariance matrix.

The penalized approach is done by adding a-priori information on the covariance matrix Σ as:

Σ = Σ + ηP.
P is the penalizer which contains a-priori behavior and η is the strength of the penalization. It also increases the stability of convergence (in case of singularity of matrix Σ).

The regularized discriminant analysis (RDA) is more complex, it combines both LDA and QDA behavior. We refer to [START_REF] Friedman | Regularized Discriminant Analysis[END_REF] for a complete description of the technique.

High dimensional data analysis

The problem of high-dimensional data arises when the cardinal of the set n is small compared to the dimension of the samples q. Although regularized discriminant analysis increases the stability, other methods are found in the literature as dimension reduction and parsimonious methods.

We just mention here the so-called High Dimensional Discriminant Analysis (HDDA) [START_REF] Bouveyron | High-Dimensional Discriminant Analysis[END_REF] which is part of the parsimonious models. It assumes that the covariance matrix Σ is approximated by: Σ = Q∆Q , where Q is the matrix of eigenvectors from Σ and more specifically, ∆ is the diagonal matrix which contains eigenvalues with some simplification. An example is a diagonal matrix which contains only two eigenvalues (λ 0 and λ 1 ): ∆ = diag(λ 0 , . . . , λ 0 , λ 1 , . . . , λ 1 ) with the same size as Σ. We refer to [START_REF] Bouveyron | Model-Based Clustering and Classification for Data Science: With Applications in R[END_REF]Chapter 8] for a wide overview of the dimensionality problem and applications.

The dimensionality issue in this manuscript does not refer to high dimensional data but to large data-sets (i.e. n large enough to encounter numerical issues as the processing of Terabytes of data in satellite images).

Matrix-variate Gaussian distribution

Previously, we illustrated distributions and examples on vectors (E = R q ). In the context of satellite image timeseries, both spectral, temporal and spatial information are provided from the satellites and it is interesting to combine them in order to increase the precision of the analysis.

In the following we decide to illustrate the case where two dimensions are used jointly. Let p be the number of spectral bands and q the number of observations (time-stamps), then the observation Y is a p × q real matrix, we write Y ∈ E with E = M p,q (R). M p,q (R) denotes the set of real matrices sized p × q. One may stack all vectors (see Definition 2.1) in order to retrieve the vector case (with E = R pq ) but the dimensions could become too large and numerical issues may appear.

The matrix-variate normal distribution has been used in different contexts, for example in electro-encephalography [START_REF] Gonzalez-Navarro | Spatio-Temporal EEG Models for Brain Interfaces[END_REF][START_REF] Spinnato | Finding EEG space-time-scale localized features using matrix-based penalized discriminant analysis[END_REF] where the study combines time samples and sensors (or spatial and spectral parts [START_REF] Zhu | EEG Signal Classification Using Manifold Learning and Matrix-Variate Gaussian Model[END_REF]), remote sensing [START_REF] Glanz | An expectation-maximization algorithm for the matrix normal distribution with an application in remote sensing[END_REF], where the authors combine time samples and spectral information, or dimension reduction [START_REF] Greenewald | Robust kronecker product PCA for spatio-temporal covariance estimation[END_REF] with time samples and spatial locations. [START_REF] Manceur | Maximum likelihood estimation for the tensor normal distribution: Algorithm, minimum sample size, and empirical bias and dispersion[END_REF] shows that it may be extended to higher order tensor distributions.

Definitions and properties

For the sake of self-containedness, the vec operator and the Kronecker product are defined (Definitions 2.1 and 2.2). Then the matrix-variate Gaussian distribution is formerly introduced in Definition 2.3.

C ⊗ D =             c 11 D . . . c 1n D . . . . . . . . . c m1 D . . . c mn D             . (2.7)
Remark. Let us highlight that, from (2.7), the matrix C ⊗ D does not allow us to retrieve C and D uniquely. Indeed, let η > 0, then C ⊗ D = (η)C ⊗ (1/η)D. With higher-order tensor distribution as in [START_REF] Manceur | Maximum likelihood estimation for the tensor normal distribution: Algorithm, minimum sample size, and empirical bias and dispersion[END_REF], the identifiability issue is increasing.

Besides, we refer to [START_REF] Schott | Matrix Analysis for Statistics[END_REF] for several properties of these two operators. From [START_REF] Dawid | Some matrix-variate distribution theory: Notational considerations and a Bayesian application[END_REF], a popular approach to reduce the dimension of the problem is to consider uncorrelated lines and columns of Y ∈ M p,q (R). This is the hypothesis of separability between lines and columns, illustrated in Figure 2.3. With a Gaussian assumption [START_REF] Dawid | Some matrix-variate distribution theory: Notational considerations and a Bayesian application[END_REF][START_REF] Srivastava | Models with a kronecker product covariance structure: Estimation and testing[END_REF], it implies the definition of two covariance matrices: Σ of size q × q which defines the covariance between two columns of Y and Λ of size p × p between two lines. This illustration does not show any mean, however it is easy to shift the matrix using a mean matrix M. Definition 2.3 (Matrix-variate Gaussian distribution, [START_REF] Dawid | Some matrix-variate distribution theory: Notational considerations and a Bayesian application[END_REF]). Let Y be a random matrix of size p × q. Y follows a matrix-variate Gaussian distribution, denoted by MN p,q (M, Σ, Λ) with mean M and covariance matrices Λ (sized p × p) between lines and Σ (sized q × q) between columns if and only if y = vec(Y), the vector of size pq, follows a pq-variate Gaussian distribution with mean µ = vec(M) and covariance matrix Σ ⊗ Λ. We write:

Y 1,1 Y 1,2 . . . Y 1,l . . . Y 1,q Y 2,1 Y 2,2 . . . Y 2,l . . . Y 2,q . . . . . . . . . . . . Y k,1 Y k,2 . . . Y k,l . . . Y k,q . . . . . .
Y ∼ MN p,q (M, Σ, Λ) if and only if y ∼ N pq (µ, Σ ⊗ Λ).
Following this definition, since y has its associated multivariate Gaussian probability density function (pdf), it is possible to define a pdf on matrix-variate Gaussian random variable Y. Proposition 2.1. Let Y ∼ MN p,q (M, Σ, Λ), the probability density function of Y is:

p(Y) = (2π) -pq/2 |Σ| -p/2 |Λ| -q/2 etr -1/2(Y -M)Σ -1 (Y -M) Λ -1 ,
where etr(.) denotes the exponential of the trace.

Proof. By Definition 2.3, y = vec(Y) ∼ N pq (µ = vec(M), Σ ⊗ Λ), then the pdf of y is: p(y) = (2π) -pq/2 |Σ ⊗ Λ| -pq/2 exp - 1 2 (y -µ) {Σ ⊗ Λ} -1 (y -µ) .
On one hand: |Σ ⊗ Λ| -pq/2 = |Σ| -p/2 |Λ| -q/2 thanks to the Kronecker product property with determinant [155, Chapter 8]. On the other hand, thanks to the properties between vec operator, the Kronecker product and the Trace (Tr) in [START_REF] Schott | Matrix Analysis for Statistics[END_REF]Theorem 8.12]:

- 1 2 (y -µ) {Σ ⊗ Λ} -1 (y -µ) = - 1 2 vec(Y -M) {Σ ⊗ Λ} -1 vec(Y -M), = Tr - 1 2 (Y -M)Σ -1 (Y -M) Λ -1 .
Let us add that the separability is not restricted to the Gaussian distribution. [START_REF] Thompson | Classification With the Matrix-Variate-t Distribution[END_REF] proposed a classification method with matrices following a t-distribution.

Maximum likelihood estimation

Let Y ∈ E a random matrix and z ∈ {1, . . . , C} a discrete random variable. The generative model is defined by:

Y|z = c ∼ MN p,q (M c , Σ c , Λ c ), (2.8)
with parameters θ c = {M c , Σ c , Λ c } and c ∈ {1, . . . , C}.

For class c we have the set S c = {Y i , z i = c} n c i=1 of n c i.i.d. samples, distributed according to a matrix-variate Gaussian distribution with mean M c and covariance matrices Σ c and Λ c (2.8). From Proposition 2.1 and the linearity of the trace, the likelihood is:

(2π) -nc pq 2 |Σ c | -n c p/2 |Λ c | -n c q/2 etr         - 1 2 i|z i =c (Y i -M c )Σ -1 c (Y i -M) Λ -1 c         .
(2.9)

The optimal parameters θc,MLE = { Mc , Σc , Λc } of (2.9) are (see [START_REF] Spinnato | Modèles de covariance pour l'analyse et la classification de signaux électroencéphalogrammes[END_REF]):

Mc = 1 n c i|z i =c Y i ,
and each covariance matrix is estimated given the other covariance matrix as:

Σc = 1 n c p i|z i =c (Y i -M) Λ-1 c (Y i -M), Λc = 1 n c q i|z i =c (Y i -M) Σ-1 c (Y i -M) .
However, this solution is, firstly, not unique because of identifiability issues in the Kronecker product. Secondly, convergence is obtained by iteratively estimating both matrices.

Flip-flop estimation

The two covariance matrices are estimated iteratively using the so-called flip-flop algorithm (see Algorithm 1, steps 3) [START_REF] Dutilleul | The MLE algorithm for the matrix normal distribution[END_REF][START_REF] Mardia | Spatial-temporal analysis of multivariate environmental monitoring data[END_REF][START_REF] Spinnato | Modèles de covariance pour l'analyse et la classification de signaux électroencéphalogrammes[END_REF]. In the following we consider a normalizing constant on Σ c written r(Σ c ) for class c ∈ {1, . . . , C}, i.e. Σ has a norm equal to 1. One may also consider Λ c which is normalized. The normalizing constant is usually a norm of a matrix, it appears in steps 3.b and 3.c of Algorithm 1.

Algorithm 1: Flip-flop algorithm

Input : Sample {(Y i , z i ) ∈ R p×q × {1, . . . , C}, i = 1, . . . , n} and initialization (Λ (0) 0 , . . . , Λ (0) C ). Output: θc,MLE , c = 1, . . . , C. 1 for c = 1 to C do 2 (1) Initialize k ← 1; 3 (2) Compute Mc ← 1 n c i|z i =c Y i ; 4 repeat 5 (3.a) Update Σ(k) c ← 1 n c p i|z i =c (Y i -M) {Λ (k-1) c } -1 (Y i -M); 6 (3.b) Compute η ← r Σ(k) c ; 7 (3.c) Update Σ(k) c ← Σ(k) c /η; 8 (3.d) Update Λ(k) c ← 1 n c q i|z i =c (Y i -M){ Σ(k) c } -1 (Y i -M) ; 9
until until (2.9) as converged;;

The normalizing constant r(.) may be, for example, the Frobenius norm (η

= Σ c F = Tr(Σ c Σ c ))
. An other solution is to remove the normalization (η = 1) and add a special initialization of Σ to ensure that (Σ) q,q = 1 as a constraint. The latter one has been presented and demonstrated in [START_REF] Srivastava | Models with a kronecker product covariance structure: Estimation and testing[END_REF], they have also proved that the associated algorithm converges to a unique extrema.

In our context, the use of this distribution will be extended using Gaussian processes in Chapter 5 with normalizing constant as in flip-flop.

Gaussian processes

This section focuses on Gaussian processes for both regression and classification tasks. Gaussian process extends the multivariate Gaussian distribution to a continuous set of inputs (also known as function-space view [START_REF] Williams | Gaussian processes for machine learning[END_REF]). Recall that S = {(y i , z i )|(y i , z i ) ∈ E × {1, . . . , C}} n i=1 denotes a set of n i.i.d. samples. In Section 2.2.1, kernel functions are introduced as they can be used within Gaussian processes, then the Gaussian processes are defined and presented from the function-space view in Section 2.2.2, Gaussian processes for classification are described in Section 2.2.3. Finally Section 2.2.4 introduces multi-outputs Gaussian processes.

Kernel functions

Before introducing Gaussian processes, this Section presents a key element which measures the similarity between two samples in space E, denoted by y i and y j , (i, j) ∈ {1, . . . , n} 2 . Kernels functions, or operators, are symmetric positive-definite functions.

Definition 2.4 (positive semi-definite functions). Let E be an arbitrary finite or infinite space. Let K : E × E → R be a binilear function. K is said to be positive semi-definite (psd) if and only if, for all y ∈ E:

K(y, y) ≥ 0.
Positive-definiteness holds if and only if K(y, y) = 0 ⇔ y = 0.

In a so-called "feature space" (y finite or infinite) denoted by H, kernel functions are equivalent to an inner product.

Definition 2.5 (Kernel). Let H be a Hilbert space with the associated inner product ., . H , let φ : E → H be a linear or non-linear mapping function. Let y i and y j two samples in E, the kernel function

K : E × E → R is defined by: K(y i , y j ) = φ(y i ), φ(y j ) H .
Let K the matrix of size n × n defined by (K) i, j = K(y i , y j ), (i, j) ∈ {1, . . . , n} 2 , the matrix K is called the Gram matrix.

Theorem 2.1. K is a positive (semi-)definite function if and only if the Gram matrix is symmetric positive (semi-)definite.

Theorem 2.1 is a weaker result from the Mercer's theorem [START_REF] Mercer | Functions of positive and negative type, and their connection the theory of integral equations[END_REF]. The family of positive (semi-)definite kernels is wide. Firstly, kernels can be separated into two categories, stationary and non-stationary kernels. In both cases, we define ., . E an inner product on space E and ||.|| E the associated norm. Definition 2.6 (Stationary kernels). Let y i and y j two samples from E. A stationary kernel is a kernel function K which only depends on a distance d E in space E between the two samples:

K(y i , y j ) = K d E (y i -y j ) .
Among the family of stationary kernels, some examples are presented in the following. Firstly the Matérn kernel is, with a positive parameter ν, defined by:

K(y i , y j ) = 2 1-ν Γ(ν)        √ 2ν y i -y j E        ν K ν        √ 2ν y i -y j E        , y i , y j ∈ E, (2.10)
where Γ is the Gamma function and K ν a modified Bessel function. is a length scale which measures, for a covariance operator, how strong two samples are correlated. Taking ν → ∞ in (2.10) gives us the well-known Squared-Exponential kernel, or Radial Basis Function (RBF):

K(y i , y j ) = exp       - y i -y j 2 E 2 2       , y i , y j ∈ E. (2.11)
Many other examples of stationary kernels can be found in [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 4]

) RBF kernel ( = 0.25) Matérn ( = 0.25, ν = 1.5) 1/2 DP (d = 1, σ 0 = 0.5)
Figure 2.4: Different kernels K(y, y 0 ), evaluated for y ∈ [0, 2] and y 0 = 1: (full line) the RBF kernel with a lengthscale = 0.25, (dashed line) the Matérn kernel with ν = 1.5 and same length-scale as RBF, and (dotted line) the dot product (DP) kernel at order 1 and a constant σ 0 = 0.5, this kernel's outputs are divided by 2 to fit in window.

The non-stationary family is wider, we just mention here the dot product kernel defined at order d and up to an additive constant σ 2 0 by:

K(y i , y j ) = y i , y j E + σ 2 0 d , y i , y j ∈ E.
If d > 1, we also described them as polynomial kernels are defined [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 4]. Other non-stationary kernels can be found in the literature such as neural network kernels [START_REF] Williams | Computation with infinite neural networks[END_REF] or spectral kernels [START_REF] Remes | Non-Stationary Spectral Kernels[END_REF]. We also mention here that colored noise (opposed to white noise) can be represented by a non-stationary kernel. Figure 2.4 illustrates some kernels for a one-dimensional space (E = R). Kernels can also be used in non-parametric techniques for classification [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 6]. If the data are linearly not separable (for example using LDA or QDA classifiers), the mapping function, introduced in Definition 2.5, can be used to transform the data into a new space, called feature space, where the data are separable. The mapping function can be unknown, this technique is also known as kernel trick and based on the Mercer's theorem [START_REF] Mercer | Functions of positive and negative type, and their connection the theory of integral equations[END_REF].. Support Vector Machines [START_REF] Williams | Gaussian processes for machine learning[END_REF]Section 6.4] from the preview in Chapter 1 are using this principle.

Kernels methods are popular among non-parametric techniques in machine learning [START_REF] Hofmann | Kernel methods in machine learning[END_REF] and can be used on vectors (E = R d ) but also with more complex data (for example when d → ∞). Thus, kernels can be defined on many different type of data such as strings [START_REF] Lodhi | Text classification using string kernels[END_REF], graphs [START_REF] Kriege | A survey on graph kernels[END_REF] or vector-valued functions [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF][START_REF] Flaxman | Scalable high-resolution forecasting of sparse spatiotemporal events with kernel methods: A winning solution to the NIJ "Real-Time Crime Forecasting Challenge[END_REF]. Different types of benchmarks where kernels are applied are presented in [START_REF] Bouveyron | Kernel discriminant analysis and clustering with parsimonious gaussian process models[END_REF].

Gaussian processes for regression

Gaussian processes are part of stochastic processes. In a function-space view [START_REF] Williams | Gaussian processes for machine learning[END_REF], its main advantage is to generalize the distribution on vectors to a continuous process. Definition 2.7 (Gaussian process). A (uni-dimensional) Gaussian process (GP) is a stochastic process y such that any finite-dimensional marginal follows a multivariate Gaussian distribution.

We write:

Y ∼ GP(m, K),
where m(t) = E(Y t ) is the mean function and of size q is distributed according to a Gaussian distribution as:

K(t, t ) = E [(Y t -m(t))(X t -m(t ))]
y =                  Y t 1 Y t 2 . . . Y t q                  ∼ N q                                      m(t 1 ) m(t 2 ) . . . m(t q )                  ,                     K(t 1 , t 1 ) K(t 1 , t 2 ) . . . K(t 1 , t q ) K(t 2 , t 1 ) K(t 2 , t 2 ) . . . . . . . . . . . . . . . . . . K(t q , t 1 ) . . . . . . K(t q , t q )                                         .
In the following we write the mean µ = [m(t 1 ), m(t 2 ), . . . , m(t q )] and the covariance matrix Σ such that (Σ) i, j = K(t i , t j ), i, j ∈ {1, . . . , q}. Figure 2.5 illustrates two samples Y 1 and Y 2 from a GP with zero mean and radial basis function (RBF) covariance operator. y 1 is a marginal of Y 1 and y 2 a marginal of Y 2 on the same set of indices. The covariance function can be defined by any kernel functions as the symmetric positive semi-definite properties ensure the definition of covariance, see [START_REF] Williams | Gaussian processes for machine learning[END_REF].

Bayesian regression problem

Let us now consider a vector y ∈ R q that is a marginal of size q from a Gaussian process Y:

Y ∼ GP(m, K).
The vector of inputs is denoted by t ∈ T q . The regression problem is to find a mean function m and a covariance operator K, from which the input-output data {t, y} are sampled, to predict new outputs from an unobserved input t ∈ T .

The Bayesian framework considers the GP distribution on the function Y as a prior and a new sample y is an output from the input t by the process Y [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 2]. As a consequence, as y and y are marginals from the same Gaussian process, we have:

y y ∼ N q+1 µ m(t ) , Σ k k K(y , y ) , (2.12)
with µ = [m(t 1 ), . . . , m(t q )] , (Σ) i, j = K(t i , t j ) and k = [K(t , t 1 ), . . . , K(t , t q )]. If we consider a vector of q outputs y from inputs t , from (2.12) it is straightforward that:

y y ∼ N q+q µ µ , Σ K K Σ , 0
µ = [m(t 1, ), . . . , m(t q , )], (K) k, j = K(t k, , t j ) and (Σ ) k,l = K(t k, , t l, )
where (k, l) ∈ {1, . . . , q } 2 and j ∈ {1, . . . , q}.

By the properties of Gaussian distributions, y |t , y, t is a q -variate Gaussian vector with mean and covariance matrix:

E(y |t , y, t) = µ + KΣ -1 (y -µ) cov(y |t , y, t) = Σ -KΣ -1 K .
(2.13) Figure 2.6 illustrates the mean (E(y |t , y, t)) and 95% confidence interval (±1.96 V(y |t , y, t)) on T = [0, 1], computed thanks to (2.13) with m(t) = 0 and K a RBF kernel function. We refer to [190, Chapter 2, Algorithm 2.1] for a practical implementation.

Gaussian processes are well-known in regression studies [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF][START_REF] O'hagan | Curve fitting and optimal design for prediction[END_REF][START_REF] Roberts | Gaussian processes for timeseries modelling[END_REF]. They are used for example to smooth observations [START_REF] Tronicke | Steering kernel regression: An adaptive denoising tool to process GPR data[END_REF] or in the context of earth-observation to analyse data [START_REF] Camps-Valls | A Survey on Gaussian Processes for Earth-Observation Data Analysis: A Comprehensive Investigation[END_REF]. Extensions in the regression context are proposed as Manifold GP [START_REF] Calandra | Manifold Gaussian Processes for regression[END_REF], or GP with constraints [START_REF] Koyejo | Constrained Gaussian Process Regression for Gene-Disease Association[END_REF][START_REF] Pensoneault | Nonnegativity-enforced Gaussian process regression[END_REF].

Finally, non-Gaussian processes are also studied in the literature. One of the most promising distribution is to use Student-t processes [START_REF] Shah | Student-t Processes as Alternatives to Gaussian Processes[END_REF]. A student-t distribution allows heavier tails as presented in Figure 2.1, they are robust w.r.t. outliers in the latent function. In [START_REF] Chen | Multivariate Gaussian and Student-t process regression for multioutput prediction[END_REF] the authors presented multivariate Student-t processes and compared them to Gaussian distributions (see Section 2.2.4).

This manuscript does not go any further in the field of regression as the context of our work is oriented towards classification. However, it is important to introduce GPs in a function-space view and their conditional distribution on a regression framework as they are used in Chapter 4 for temporal reconstruction of satellite images with Gaussian processes.

Gaussian processes for classification

Gaussian processes are flexible and can be adapted to different situations. In this section, our aim is to transform real-valued inputs into C discrete outputs (or classes). To this end, Gaussian processes are used in a discriminative approach, i.e., given c ∈ {1, . . . , C}, one GP models the posterior probability P(z = c|y).

Gaussian processes as latent functions

As introduced, the use of Gaussian processes for classification is based on a set of centered GP priors on a functional f , called the latent function, from an input space to the space of real numbers ( f : E → R).

The main idea of GP for classification is to define a conditional Bernoulli distribution thanks to a transformation of the latent function f through a logit function σ : R → [0, 1]. Let c be a given class, the conditional Bernoulli distribution for class c is denoted by P(Z = c|y) (with the previous notations) and is defined by:

P(Z = c|y) = σ( f c (y)), (2.14)
with f c ∼ GP(0, K c ), i.e. f c is completely described by the covariance operator K c .

Let us now consider C classes. Recall that S is the set defined by S = {Y, z} ∈ R n×q × {1, . . . , C} n . We define f as the vector of latent functions of size nC defined by:

f = [ f 1 1 , . . . , f 1 n , f 2 1 , . . . , . . . , f C 1 , . . . , f C n ] ,
where f c i is the latent function for sample (Y) i of size q for class c. Then f ∼ GP(0, K) where K = diag(K 1 , . . . , K C ) as the C processes are assumed to be independent and K c is the covariance operator in class c.

Inference and approximations for the a-posteriori probabilities

Inference is complex, it aims to return the a-posteriori probabilities for each class c ∈ {1, . . . , C}. It is done in two steps [190, Chapter 3]:

1. Compute the likelihood density distribution over the latent function f = [ f 1 , . . . , f C ] for a test case y :

p(f |Y, z, y ) = p(f |Y, y , f)p(f|Y, z)df, (2.15)
where p(f|Y, z) = p(z|f, Y)p(f|Y)/p(z|Y) is the posterior over the latent variables for known data.

2. Produce the conditional Bernoulli distribution on y for all classes c ∈ {1, . . . , C}:

P(Z = c|y ) ∝ σ( f c )p( f c |Y, z, y )d f , (2.16)
with the likelihood defined in (2.15). The notation P refers to the use of f = f, i.e. the latent function which maximizes (2.15).

In this situation, the integral (2.16) is intractable for some choice of sigmoid functions. Additionally the likelihood (2.16) is non-Gaussian [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 3]. To overcome this problem, a solution is to approximate the posterior with a Gaussian distribution, two solutions are found in the literature: the Laplace approximation [START_REF] Williams | Computation with infinite neural networks[END_REF], we refer to [START_REF] Williams | Gaussian processes for machine learning[END_REF]Algorithms 3.3 and 3.4] for practical implementation, and the expectation propagation (EP) method [START_REF] Minka | A family of algorithms for approximate Bayesian inference[END_REF][START_REF] Nickisch | Approximations for binary gaussian process classification[END_REF], see [START_REF] Williams | Gaussian processes for machine learning[END_REF]Algorithms 3.5 and 3.6]. The following describes the two approximations.

On one hand, the Laplace approximation approaches the posterior density p(f|Y, z) over f by a Gaussian approximation, denoted by q(f|Y, z) = N(f| f, A -1 ) [START_REF] Williams | Gaussian processes for machine learning[END_REF]Section 3.4] where f is the latent function f which maximizes the likelihood (2.15) and A = ∇∇ f= f log p(f|Y, z) is the Hessian of the negative log posterior evaluated at f. The new likelihood is Gaussian and denoted by q(f |Y, z, y ). In [START_REF] Williams | Gaussian processes for machine learning[END_REF]Section 3.5], they showed that the predictive Gaussian likelihood mean and covariance are:

E(f |Y, z, y ) = Q (z c -P(Z = c|y 1 ), . . . , P(Z = c|y n ) ), cov(f |Y, z, y ) = Σ + Q K -1 (K -1 + W) -1 K -1 Q , (2.17)
where Q is a block diagonal matrix defined by

Q = diag(k 1 (y ), . . . , k C (y )) of size Cn × C with k c (y ) = [K c (y 1 , y ), . . . , K c (y n , y )] , K = diag(K 1 , . . . , K C
) is the previous diagonal matrix where each K c is the Gram matrix in class c and W is a diagonal matrix [START_REF] Williams | Gaussian processes for machine learning[END_REF]Equation (3.38)]. Finally Σ is a diagonal matrix of size C × C where (Σ) c,c = K c (y , y ) -k c (y )K -1 c k c (y ). On the other hand, the Expectation propagation algorithm (EP) has been presented by [START_REF] Minka | A family of algorithms for approximate Bayesian inference[END_REF]. A main difference with the Laplace approximation is that the EP algorithm is used in the binary classification case (C = 2), corresponding to (2.14) where σ can be changed easily. It means that only one posterior probability is necessary, for example we consider the first class P(Z = 1|y ) (as P(Z = 2|y ) = 1 -P(Z = 1|y )), i.e. one latent function is necessary and f = [ f 1 , . . . , f n ] . The non-Gaussian posterior is now given by the Bayes' rule as:

p(f|Y, z) = p(f|Y) p(z|Y) n i=1 p(y i | f i ),
where the likelihood p(y i | f i ) remains non-Gaussian. However, with the use of the probit likelihood (cumulative density function of a standard Gaussian distribution), the integral of the posterior can be computed. EP framework approximates the likelihood thanks to a local likelihood [START_REF] Williams | Gaussian processes for machine learning[END_REF]Section 3.6]. We also refer to [190, Section 3.6, Algorithms 3.5 and 3.6] for a complete and practical implementation.

To understand the limitation of Gaussian processes for classification when dealing with large data-sets (as land use of land cover applications), it is important to highlight that these algorithms scale in O((C + 1)n 3 ), i.e. cubic w.r.t. the number of samples. Additionally, within the Laplace approximation algorithm, a Monte-Carlo loop has to be done to obtain the posterior probabilities (2.14).

Large data-sets approximations and scaling

The following presents how, in the literature, the problem of large data-sets with GP for classification has been tackled.

Recall that, for regression and classification, the computation complexity is O(n 3 ) as the Gram matrix K has to be inverted and depends on the number of samples n in the set of data S. In order to reduce complexity, a first case to consider is to work on a well chosen subset of S. Also called sub-sampling, it takes a smaller number m of data where m n. The algorithm scales in O(m 3 ) and [START_REF] Hayashi | On Random Subsampling of Gaussian Process Regression: A Graphon-Based Analysis[END_REF] presented some theoretical results on the error made by this strategy. A second case is to sparsify the matrix K to make the inversion easier. One example is the Nyström technique, it replaces the matrix K by K defined in [START_REF] Williams | Using the nyström method to speed up kernel machines[END_REF] by:

K = K n,m K -1 m,m K m,n , (2.18)
where K p,q are blocks of size p × q from K constructed by the evaluation of the covariance operator K between the first p samples and the first q samples. In practice, m n and the inversion cost of K in (2.18) is O(m 2 n). We refer to [START_REF] Gittens | Revisiting the nyström method for improved large-scale machine learning[END_REF][START_REF] Williams | Using the nyström method to speed up kernel machines[END_REF] and [START_REF] Williams | Gaussian processes for machine learning[END_REF]Section 8.3] for discussions on the limits of this method and a criteria to select m. An other approach is to approximate the marginal likelihood. A typical one is the use of the variational inducing point framework [START_REF] Hensman | Scalable Variational Gaussian Process Classification[END_REF] where the likelihood is lower bounded by a quantity with lower complexity. It results in a complexity in O(m 3 ).

A complete review, considering these approximations and many more can be found in [START_REF] Liu | When Gaussian process meets big data: A review of scalable GPs[END_REF]. They also reviewed extensions as deep GP [START_REF] Damianou | Deep Gaussian Processes[END_REF] or more complex structures as multi-task GP [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF] which is reviewed in the next Section. Finally, let us highlight that similar problems are found in the use of GP for regression as in Section 2.2.2. A review of scalable approximations can be found in [START_REF] Quiñonero-Candela | A unifying view of sparse approximate gaussian process regression[END_REF].

The GP classifiers will not be used as a latent function to model the a-posteriori probability. The use of GPs in our work will be based on a generative approach. We will show that the inference algorithm scales linearly w.r.t. n.

Multi-output Gaussian processes

This section presents a brief review of GP techniques to deal with multi-output processes. Let Y(t) = [Y 1 (t), . . . , Y p (t)] , t ∈ T a process with p outputs. This section aims to present GP models that take into account correlation between outputs Y b , b ∈ {1, . . . , p}.

The models of multi-output GP were originally introduced in the geostatistics literature [START_REF] Cressie | Statistics For Spatio-Temporal Data[END_REF][START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF]. For any given P ∈ N (let us highlight that the number of components (P) and the number of outputs (p) are not necessarily equals), the p models are composed by P latent functions, or processes, such that the b-th function, b ∈ {1, . . . , p}, is a linear combination of the latent functions. ∀(t, t ) ∈ T 2 :

Y b (t) = P =1 a b, u (t), (2.19)
where the latent functions {u (t)} P =1 are centered GP with covariance operator given by cov(u (t), u (t )) = k (t, t ) if = , zero otherwise. This model is known as the Linear Model of Coregionalization [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF] (or LMC). The LMC covariance operators are known to have a sum of separated kernels [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF], i.e. the covariance operator in time does not model the covariance between (outputs) and the linear combination of outputs does not depend on the time.

Let P ≤ P be the total number of functions {u k (t)} P k=1 sharing the same covariance [3, Equation ( 18)], the general covariance operator matrix, sized p × p, associated with (2.19) is:

K(t, t ) = P =1 C k (t, t ), (2.20)
where C , a p × p kernel matrix, is known as the coregionalization matrix, C = A A where (A ) b,k = a k b, is a matrix of size p × P . This model is used mostly in regression context with data of fixed length, [START_REF] Cheng | Sparse multi-output Gaussian processes for online medical time series prediction[END_REF] uses LMC on time-series to monitor patient state using clinical covariates as features. [START_REF] Lark | Fitting a linear model of coregionalization for soil properties using simulated annealing[END_REF] uses spatial information on one hand and soil properties on the other hand to model soil's variograms.

For a fiven ∈ {1, . . . , P}, let B be a p × p matrix which does not depend on . A particular case of LMC is to consider C equal to B up to a multiplicative constant c : C = Bc . Then

K(t, t ) = B P =1 c k (t, t ) = B k(t, t ). (2.21)
This model is known as Intrinsic Coregionalization Model [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF] (ICM) where the covariance operator matrix (2.20) is given in (2.21). ICM has been derived under multiple names in the machine learning literature where, most of the time, the introduction of the model is different. An example, from [START_REF] Bonilla | Multi-task Gaussian Process Prediction[END_REF], is the multi-task Gaussian processes model. It assumes that the family of latent functions {Y b } p b=1 in one direction has a GP prior to induce correlation between "tasks". For (b, b ) ∈ {1, . . . , p} 2 , (t, t ) ∈ T 2 , the model is defined by the correlation:

Y b (t), Y b (t ) = (B Y ) b,b k(t, t ), Y b (t) ∼ N( f b (t), σ 2 b ),
where Y b (t) is the b-th output at input t ∈ T , B Y is the covariance between outputs and k is the covariance operator on inputs. Other works derived from ICM are presented in [START_REF] Bilionis | Multi-output local Gaussian process regression: Applications to uncertainty quantification[END_REF][START_REF] Moreno-Muñoz | Heterogeneous Multi-output Gaussian Process Prediction[END_REF]. These models are mostly applied in the context of regression with inputs of same length, we refer to [START_REF] Chen | Multivariate Gaussian and Student-t process regression for multioutput prediction[END_REF][START_REF] Liu | Remarks on multi-output Gaussian process regression[END_REF] for recent reviews with various applications as species distribution models [START_REF] Ingram | Multi-output Gaussian processes for species distribution modelling[END_REF] or remote sensing [START_REF] Pipia | Fusing optical and SAR time series for LAI gap filling with multioutput Gaussian processes[END_REF]. Some strategies with large data-sets are also proposed as in [START_REF] Li | Hierarchical Gaussian Processes model for multi-task learning[END_REF] but most of them remain quadratic or cubic w.r.t. n the number of samples. An other example is the Semi-parametric Latent Factor Model [START_REF] Teh | Semiparametric latent factor models[END_REF] for multiple response variables using GP as a prior, it assumes the same structure as (2.21) but the complexity of the conditional dependency between responses imposes to use sparse approximations of the posterior. Nowadays, the multioutputs GP are used also for classification [START_REF] Dezfouli | Scalable Inference for Gaussian Process Models with Black-Box Likelihoods[END_REF] where the principle is similar to GP for classification (Section 2.2.3) or [START_REF] Skolidis | Bayesian Multitask Classification With Gaussian Process Priors[END_REF] with a variational approach.

Recently, extensions of LMC are proposed as the convolutional Gaussian processes. The convolutional Gaussian processes have been defined in [2]. Starting from (2.19) and adding a convolution between covariance kernels on outputs {a b (t)} p b=1 (smoothing kernels in [2]) and P latent functions {u } P =1 . Let (t, s) ∈ T 2 , one of the most general form of the model is:

f b (t) = P =1 T a b, (t -s)u (s)ds, (2.22)
with b ∈ {1, . . . , p} and ∈ {1, . . . , P} and {u (t)} P =1 are centered GP. It results in a dependency between outputs which is varying with inputs.

Other extensions from convolutional Gaussian processes are presented in [START_REF] Boyle | Dependent Gaussian Processes[END_REF] where outputs share similar latent functions. In [START_REF] Cardona | Convolved Multi-output Gaussian Processes for Semi-Supervised Learning[END_REF] they added a graph regularization and proposed an Expectation Maximization (EM) algorithm for inference. These extensions are often applied in the context of regression.

Statistical modelling for Sentinel-2 satellite image time-series

This chapter has presented a state-of-the-art in supervised model-based classification with an in-depth study of Gaussian distributions. These techniques are presented to understand both discriminative analysis and matrixvariate Gaussian distributions as they are keys to understand the contributions of this thesis. To the best of our knowledge, no technique is well suited for application to SITS with missing time stamps at large scale with the challenges arise from Sentinel-2 mission.

This chapter has also presented Gaussian processes for regression and classification. Learning Gaussian processes for classification has been shown to be cubic with respect to the number of samples in the provided training set. The multi-output extensions of Gaussian processes have been reviewed, this thesis will contribute specifically to this field by formerly defining a mixture of multivariate GP which scales to large datasets with missing values. Recent extensions with convolutional GP are presented for a discussion on these models.

The following introduces the contributions. We propose two generative models based on Gaussian processes with application to classification of irregularly and unevenly sampled time-series.

French introduction

Ce chapitre a pour objectif d'introduire les deux contributions de cette thèse. Ces deux contributions peuvent être spécifiées par un modèle de processus gaussiens. Ces deux modèles, nommés MIMGP et M2GP, font l'hypothèse, pour le premier, d'indépendance entre les bandes spectrales des SITS de Sentinel-2, et pour le second, d'une combinaison linéaire de processus gaussiens latents ayant le même opérateur de covariance temporel.

Ces deux modèles de classification sont introduits pour modéliser d'un point de vue fonctionnel les séries temporelles. Considérer des séries temporelles observées à des instants différents revient alors à considérer une discrétisation temporelle de la fonction modélisée. Ensuite la reconstruction de données manquantes (nuages ou ombre) est détaillée.

Ce chapitre se conclut sur la présentation des annexes consacrées à la description du code, pour lequel le choix du langage python a été fait et optimisé avec le langage cython (langage compilé).

Let Y (the notation will be argued later) be a random matrix of size p × q, it contains the SITS reflectance at p spectral bands and q temporal acquisitions after noise removal. Let Z be a discrete random variable taking its values in {1, . . . , C}. The decision rule, from (2.1), is given by:

ĉ = arg max c P(Z = c|Y ) = arg max c π c p(Y |Z = c), c ∈ {1, . . . , C}.
Our aim is to build a generative model, i.e. to assume a distribution on the conditional random variable Y |Z = c where Y ∈ E is the reflectance, and is distributed conditionally to its class membership Z = c. This chapter describes how Gaussian processes can model irregularly and unevenly sampled time-series and links our contributions to state of the art. Section 3.1 explains how Gaussian processes overcome the issue of irregular samples and presents their use in the classification step. Section 3.2 shows how the models can be used to perform the continuous reconstruction of missing time-stamps using the complete mixture. Finally Section 3.3 introduces the implementation of two models.

34

Chapter 3 -Model-based classification for irregularly sampled time-series

Classification of irregularly sampled time-series

Section 3.1.1 presents the general model and explains why such model is able to handle irregularly and unevenly sampled time-series. Section 3.1.2 presents a one-dimensional irregularly sampled time-series classifier and presents the classifier for multiple spectral bands with independence assumption. Finally Section 3.1.3 presents the model with linear dependency between wavelengths.

Continuous representation of SITS

The following defines the general model behind the two contributions.

Preliminary definitions

In the following, E denotes the data space for re-sampled time-series data (E = R q ). In the following we assume to observe the time-series continuously on a given time window T = [t min , t max ]. The multi-dimensional functions with p square integrable processes is denoted by Y(t) and defined as:

Y : T → R p , t → [Y 1 (t), . . . , Y p (t)] . (3.1)
A discretized observation from a continuous process Y is called marginal and is denoted by Y . As presented previously, it is a random matrix of size p × q.

In practice, the observations are Sentinel-2 SITS. The reflectance are included in the set S with the associated class membership. Considering the i-th sample yields to the marginal matrix Y i, of size p × q i . If, ∀i ∈ {1, . . . , n}, q i = q and, moreover, the time stamps are equal, then the data space E = M p,q (R) is the space of matrices with same size. The classifiers are explored in Section 2.1.

In our context, clouds and shadows noises occur at random spatial locations and time stamps. The resulting effect is different sample sizes, i.e. q i q j , and different time stamps. An example of two simulated time-series is given in Figure 3.1: in this case q 1 = 18 and q 2 = 14. By taking into account level 2A products (surface reflectance Ỹi, on p spectral bands, the associated mask M i ) and removing noisy data (mask values above '0') yields the set S = {Y i, , z i } n i=1 of n irregularly sampled SITS where each i-th sample has its own size: p × q i as the number of clear dates may change; z i is the associated class membership to sample i, i ∈ {1, . . . , n}.

General model

The general model assumes that, conditionally to z = c: Equivalently, conditionally to Z = c, the b-th spectral band is a linear combination of the p latent processes:

Y(t) = A c W(t) + m c (t), ∀t ∈ T , (3.2
Y b (t) = m b,c (t) + p =1 (A) b, W k (t). (3.3)
where W k ∼ GP(0, K k,c ).

One-dimensional process and the Mixture of Independent multivariate Gaussian processes

The first contribution assumes that the p spectral bands are idependent, i.e. we consider p one-dimensional processes. For one spectral bands, the model (3.2) becomes: where W W 1 ∼ GP(0, Kc ). Then it is straightforward that:

Y(t) Y 1 (t) = a c W 1 (t) + m c (t), ∀t ∈ T , 3
Y ∼ GP(m c , a 2 c Kc ).
Then the problem is to consider, for a fixed spectral band, the i-th sample which is assumed to follow a multivariate normal distribution (see Definition 2.7).

Mixture of Independent multivariate Gaussian processes

This section introduces our first contribution in Chapter 4, with associated supplementary materials in Appendix B. From (3.2), the model assumes that, conditionally to Z = c, the p processes (Y 1 , . . . , Y p ) are independent, i.e. A c = I, it comes the following decision rule:

ĉ = arg max c π c p(Y |Z = c) = arg max c π c p b=1 p(Y b |Z = c),
where Y b ∼ GP(m b,c , K b,c ). It consists of p independent processes, models inference is detailed in Section 4.2.3. The classification accuracy results, computed on a separate validation set, are presented in Section 4.6 and are compared with the state of the art methods introduced in Chapter 1. State-of-the-art methods take as input linearly re-sampled time-series as presented in Figure 1.12.

Inference of model parameters (mean coefficients and kernel parameters) on the set S for is provided in Section 4.2.3. It is inferred using a maximum likelihood estimator (MLE). This estimation is possible as every sample is a marginal from a continuous process and the complete process is inferred from MLE, despite irregular time stamps for each sample i.

To the best of our knowledge, this contribution introduces the first model, namely Mixture of Independent Multivariate Gaussian Processes (MIMGP), which is able to classify irregularly sampled data. We complete the presentation of MIMGP with an application to Sentinel-2 SITS using non-noisy data, p = 10 spectral bands and time stamps distributed as in Figure 4.6. Additionally, this model is able to reconstruct time-series at any time t † ∈ T and is described in Section 3.2.

When the class memberships are unknown (unsupervised framework), the definition of an unsupervised classifier is straightforward and inference with EM (Expectation-Maximization) algorithm is simple thanks to the Gaussian assumption.

Mixture of multivariate Gaussian processes

Highlighted in the first contribution, one of the weaknesses of the MIMGP model is the independence assumption.

To this end we propose a new model which deals linearly with the p outputs of the vector function Y (see (3.1)): the mixture of multivariate Gaussian processes (M2GP, Chapter 5). Let us describe the construction of M2GP.

Our aim in Chapter 5 is to obtain the mean and the variance of the Gaussian distribution Y i, |z i = c (see Section 5.8) for the i-th matrix sample of size p × q i or, similarly, y i, = vec(Y i, ). The mean is straightforward to compute, and the covariance is given by:

cov Y b (t i j ), Y b (t i j ) = cov a c,b W b (t i j ), a b ,c W b (t i j )
where a c,b is the b-th line of A c , then, as

W b W b for b b , then: cov Y b (t i j ), Y b (t i j ) = a c,b Σ b,c,i j, j a c,b , if b = b , = 0, otherwise.
Applied to all time stamps t = [t i 1 , . . . , t i q i ] and all p spectral wavelengths, it comes:

cov(y i, |z i = c) =                          A c diag [Σ 1,c,i 1,1 , . . . , Σ p,c,i 1,1 ] A c A c diag [Σ 1,c,i 1,2 , . . . , Σ p,c,i 1,2 ] A c . . . A c diag [Σ 1,c,i 2,1 , . . . , Σ p,c,i 2,1 ] A c A c diag [Σ 1,c,i 2,2 , . . . , Σ p,c,i 2,2 ] A c . . . . . . . . . . . .                          .
The covariance matrix is of size pq i × pq i which makes the inference numerically difficult. One solution, from (3.2), is to assume that K b,c = K c , ∀b ∈ {1, . . . , p}, i.e. ∀b, W b ∼ GP(0, K c ). This defines the M2GP model (5.1). The hypothesis of the same time-covariance operator yields a Kronecker product structure which allows a significant decrease of numerical complexity. All marginals are now distributed according to a matrix-variate Gaussian distribution (see Definition 2.3).

Parameters estimation is presented in Section 5.6.2. Figure 5.6 shows significant performance improvements when adding some dependence compared with MIMGP. Finally comparisons on Sentinel-2 SITS are reported in Table 5.2 with LULC maps in Figure 5.12.

Reconstruction of noisy time stamps

An additional and promising application is the reconstruction of noisy reflectance values. An example is the use of Sentinel-2 SITS for phenological research, a review of applications can be found in [START_REF] Misra | Status of Phenological Research Using Sentinel-2 Data: A Review[END_REF]. As presented in Chapter 1, reconstruction of missing values is also an active area of study in remote sensing [START_REF] Shen | Missing Information Reconstruction of Remote Sensing Data: A Technical Review[END_REF].

Imputation with independence assumption

The imputation of missing values is possible thanks to the Gaussian marginal at a new time input t † ∈ T from the continuous process modeled by a GP. The joint distribution (2.12) allows the imputation of missing values when the class membership is known using conditional expectation and conditional covariance [17, p.63] defined in (2.13). The reconstruction of the p spectral bands is done independently using the conditional expectation (4.11).

An advantage of this reconstruction is the confidence interval computed thanks to (4.12).

Imputation with M2GP

Imputation of missing values using the M2GP model differs. Indeed, Proposition 5.3 presents the expectation and variance of the conditional distribution. The expectation does not depend on the linear dependency between spectral bands but all p reflectances are imputed at once. The covariance retrieves the Kronecker structure induced by the model. These observations are directly linked to the uncorrelated lines and columns in the matrix-variate Gaussian distribution (Definition 2.3). The reconstruction quality is the same using M2GP or MIMGP. Figure 5.6 compares the two models on toy data-sets by computing the normalized mean square error (5.21).

Imputation using the complete mixture

When the class membership is unknown, the two models allow to compute the reconstruction of missing values using the mixture of GPs. (4.13) (or (5.17) for M2GP) allows the reconstruction using the posterior probabilities to belong to all classes.

Section4.7 presents reconstructions and a comparison with state of the art non parametric smoothing techniques. As in [START_REF] Vuolo | Smoothing and gap-filling of high resolution multi-spectral time series: Example of Landsat data[END_REF], we compared the reconstruction using the mean of the conditional distribution (4.13) with the Whittaker smoother [START_REF] Eilers | A Perfect Smoother[END_REF].

Implementation -python code

The code has been written using python and the optimization has been done using cython (see Appendix C). The choice of a second language which is compiled has been done to avoid heavy "for" loops. Indeed, these loops can not be rewritten on a vectorial form where python performs well (the processing of each sample within the training or validation set induces a heavy loop). This work has been done for MIMGP but it can be adapted similarly for M2GP. However, both models are using parallel implementations on the classes.

During the tests of the code, some numerical issues arose for the estimation of the mean coefficients (Section 4.2.2) and, consequently, have interfered with the convergence of the gradient descent algorithm. The latter one uses fortran "f_min_l_bfgs_b" function [START_REF] Zhu | Algorithm 778: L-BFGS-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF] on kernel's parameters. These numerical issues are discussed in Appendix A with illustrations on a toy data-set.

Introduction

In the last decade, the successful launching of two satellites Sentinel-2 A and B offers a unique opportunity to record, analyze and predict the evolution of the Earth's land surface. Sentinel-2 mission provides a high resolution multispectral (13 spectral bands at 10m, 20m or 60m per pixel) acquisition with a 5 day revisit cycle [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]. It is planed for long-term operational observations (more than 15 years) and optical Sentinel-2 satellite image timeseries (SITS) are available to users under a free and open data policy. Such a mission provides several terabytes of worldwide data per day [START_REF] Soille | A versatile data-intensive computing platform for information retrieval from big geospatial data[END_REF].

Thanks to their spectral content and frequent update, SITS have found applications in various scientific field: water management [START_REF] Batur | Assessment of Surface Water Quality by Using Satellite Images Fusion Based on PCA Method in the Lake Gala, Turkey[END_REF][START_REF] Pereira-Sandoval | Calibration and validation of algorithms for the estimation of chlorophyll-a in inland waters with sentinel-2[END_REF], agricultural systems mapping [START_REF] Feng | Crop Type Identification and Mapping Using Machine Learning Algorithms and Sentinel-2 Time Series Data[END_REF][START_REF] Moeini Rad | Developing an Automatic Phenology-Based Algorithm for Rice Detection Using Sentinel-2 Time-Series Data[END_REF][START_REF] Useya | Comparative Performance Evaluation of Pixel-Level and Decision-Level Data Fusion of Landsat 8 OLI, Landsat 7 ETM+ and Sentinel-2 MSI for Crop Ensemble Classification[END_REF], urban area analysis [START_REF] Haas | Urban Land Cover and Ecosystem Service Changes based on Sentinel-2A MSI and Landsat TM Data[END_REF][START_REF] Iannelli | Urban Extent Extraction Combining Sentinel Data in the Optical and Microwave Range[END_REF] or ecological monitoring [START_REF] Fauvel | Prediction of plant diversity in grasslands using Sentinel-1 and -2 satellite image time series[END_REF][START_REF] Li | Estimating Plant Traits of Alpine Grasslands on the Qinghai-Tibetan Plateau Using Remote Sensing[END_REF]. However, such abundance of images raises new challenges in terms of large scale multi-resolution SITS processing [START_REF] Chi | Big Data for Remote Sensing: Challenges and Opportunities[END_REF]. Issues related to Big Earth Data were obviously explored in the remote sensing community in order to calibrate and distribute the images seamlessly [START_REF] Baetens | Validation of Copernicus Sentinel-2 Cloud Masks Obtained from MAJA, Sen2Cor, and FMask processors using reference cloud masks generated with a supervised active learning procedure[END_REF][START_REF] Kempeneers | Optimizing Sentinel-2 Image Selection in a big data context[END_REF][START_REF] Revel | Sentinel-2A and 2B absolute calibration monitoring[END_REF][START_REF] Skakun | Multispectral Misregistration of sentinel-2A Images: Analysis and Implications for Potential Applications[END_REF][START_REF] Sudmanns | Big Earth data: Disruptive changes in earth observation data management and analysis?[END_REF].

Image and signal processing issues were also investigated specifically to SITS properties. Spatial differences in terms of resolution were tackled using super resolution strategies to resample all the spectral bands at the same finest spatial resolution (10m) [START_REF] Paris | A Novel Sharpening Approach for Superresolving Multiresolution optical images[END_REF][START_REF] Ulfarsson | Sentinel-2 Sharpening Using a Reduced-Rank Method[END_REF][START_REF] Wang | Sentinel-2A Image Fusion Using a Machine Learning Approach[END_REF]. Fusion with other sensors, such as Landsat [START_REF] Shao | Deep learning-based fusion of Landsat-8 and Sentinel-2 images for a harmonized surface reflectance product[END_REF] or Sentinel-1 [START_REF] Fernandez-Beltran | Multimodal Probabilistic Latent Semantic Analysis for Sentinel-1 and Sentinel-2 Image Fusion[END_REF], was also considered to complement Sentinel-2 satellite images.

Classification of land cover or land use using machine learning techniques has emerged as a major application of SITS. Deep learning was intensively studied either using only one or few temporal acquisitions with convolutional neural networks (CNN) [START_REF] Shendryk | Deep learning for multi-modal classification of cloud, shadow and land cover scenes in planetscope and Sentinel-2 imagery[END_REF][START_REF] Syrris | Evaluation of the Potential of Convolutional Neural Networks and Random Forests for Multi-Class Segmentation of Sentinel-2 Imagery[END_REF], or using the full temporal stack of acquisitions using recurrent neural networks (RNN) [START_REF] Qiu | Local climate zone-based urban land cover classification from multi-seasonal Sentinel-2 images with a recurrent residual network[END_REF][START_REF] Rußwurm | Multi-temporal land cover classification with long short-term memory neural networks[END_REF]. Combination of RNN applied on SITS with CNN applied on SPOT-6 image was also investigated for land cover mapping over Reunion Island [START_REF] Benedetti | M 3 Fusion: A Deep Learning Architecture for Multiscale Multimodal Multitemporal Satellite Data Fusion[END_REF]. Temporal convolution was investigated in [START_REF] Pelletier | Temporal Convolutional Neural Network for the Classification of Satellite Image Time Series[END_REF]. Yet, due to the huge amount of data to be processed when large geographical area or large temporal domain (e.g. 1 year of acquisition) are considered, lighter machine learning techniques, such as Random Forest (RF), were proposed as well, with very good results in terms of classification accuracy [START_REF] Gbodjo | Toward Spatio-Spectral Analysis of Sentinel-2 Time Series Data for Land Cover Mapping[END_REF][START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. Similarly, Dersken et al. have shown that, when the number of spectro-temporal measurements per pixel is high, conventional hand-crafted features and RF perform as well as CNN in terms of classification accuracy but with a drastically reduced processing time [START_REF] Derksen | Geometry Aware Evaluation of Handcrafted Superpixel-Based Features and Convolutional Neural Networks for Land Cover Mapping Using Satellite Imagery[END_REF].

One major issue arising when dealing with large SITS, i.e. using several tracks or tiles 1 , is related to the irregular temporal sampling. Indeed, each Sentinel-2 track is acquired on different dates, and using images from different tracks results in pixels with different acquisition times. For instance, Figure 4.1 shows the acquisition dates for three Sentinel-2 tiles located over three parts of France. Even though some dates are similar, the total number of acquisitions is different from one tile to another. Furthermore, clouds and shadows, occurring at random, are another cause of irregular temporal sampling: Dates with shadow/cloud are usually considered as missing values and are not taken into account to build the feature vector of the corresponding pixels. Hence, even pixels belonging to the same tile can be represented by feature vectors of different sizes [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF].

Classical machine learning algorithms take as input a time-series of finite and constant size. Random Forest, Support Vector Machine or RNN use a vector representation of input features, i.e., a time-series, or pixel, is a collection of features stacked into a p-dimensional vector and each pixel has the same number of features. Similarly, CNN use a constant size patch representation of the image. Therefore, a pre-processing step, sometimes called missing information reconstruction, is usually applied on the SITS to recover pixels of same size [START_REF] Shen | Missing Information Reconstruction of Remote Sensing Data: A Technical Review[END_REF]. Among the various existing techniques, temporal filters, such as the Savitzky-Golay filter [START_REF] Luo | Savitzky-Golay smoothing and differentiation filter for even number data[END_REF] or the best index slope extraction (BISE) [START_REF] Viovy | The Best Index Slope Extraction (BISE): A method for reducing noise in NDVI time-series[END_REF], and parametric or non parametric curve fitting [START_REF] Beck | Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI[END_REF][START_REF] Eilers | A Perfect Smoother[END_REF][START_REF] Jonsson | Seasonality extraction by function fitting to time-series of satellite sensor data[END_REF][START_REF] Lopes | Object-based classification of grasslands from high resolution satellite image time series using gaussian mean map kernels[END_REF] are popular due to their simplicity and efficiency. At large scale (entire metropolitan France territory), linear temporal interpolation has shown to perform very well compared to spline interpolation [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. Linear interpolation may also be combined with non parametric smoothers [START_REF] Vuolo | Smoothing and gap-filling of high resolution multi-spectral time series: Example of Landsat data[END_REF].

Methods taking account both spatial and temporal information have provided interesting results, see for instance non-local filters [START_REF] Buades | A non-local algorithm for image denoising[END_REF] or deep learning based approaches [START_REF] Cresson | Optical image gap filling using deep convolutional autoencoder from optical and radar images[END_REF][START_REF] Zhang | Thick cloud and cloud shadow removal in multitemporal imagery using progressively spatio-temporal patch group deep learning[END_REF][START_REF] Zhang | Missing Data Reconstruction in Remote Sensing Image With a Unified Spatial-Temporal-Spectral Deep Convolutional Neural Network[END_REF]. However, due to their high computational cost, they are hardly applicable in national/continental scale settings [START_REF] Derksen | Geometry Aware Evaluation of Handcrafted Superpixel-Based Features and Convolutional Neural Networks for Land Cover Mapping Using Satellite Imagery[END_REF]. Furthermore, they are also known as "black-box" and hardly interpretable.

Recently, Gaussian Processes (GP) have gained attention in the remote sensing community [START_REF] Camps-Valls | A Survey on Gaussian Processes for Earth-Observation Data Analysis: A Comprehensive Investigation[END_REF][START_REF] Fauvel | Parsimonious Gaussian Process Models for the Classification of Hyperspectral Remote Sensing Images[END_REF]. GP mixed Bayesian and kernel methods to build statistical learning machines, for classification or regression problems. They have proved to be accurate and interpretable through their hyperparameters (mean and covariance functions) [START_REF] Gewali | Gaussian Processes for Vegetation Parameter Estimation from Hyperspectral Data with Limited Ground Truth[END_REF][START_REF] Morales-Alvarez | Remote sensing image classification With Large-scale Gaussian Processes[END_REF][START_REF] Pipia | Fusing optical and SAR time series for LAI gap filling with multioutput Gaussian processes[END_REF][START_REF] Verrelst | Retrieval of Vegetation Biophysical Parameters Using Gaussian Process Techniques[END_REF]. In this work, a Gaussian Process model is introduced, relying on a linear projection of the mean function on a well-chosen basis, conditionally to the classes of interest. It is able to deal with irregularly sampled pixels in the learning and classification steps as well as reconstructing missing data.

The contribution of this work is three-fold:

1. The definition of a GP model that handles irregular temporal sampling (Section 4.2), without resampling all the data on a common temporal grid.

2. Jointly to the classification, the proposed model allows for the reconstruction of pixels on any temporal grid, together with confidence intervals (Section 4.3). Furthermore, the mean and covariance functions associated with each class are available for analysis.

3. The reconstruction used the GP model with optimal parameters (from the likelihood point of view) and takes into account the class membership. Rather than relying only on the reconstruction error, as with conventional missing data reconstruction techniques, the proposed method used the classification model likelihood to control the reconstruction and its subsequent smoothing level.

The remainder of this paper is organized as follows. The proposed GP model and its estimation are described in Section 4.

Irregularly sampled Gaussian processes model

In the following S = (y i , z i ) n i=1 denotes a set of n independent and identically distributed (i.i.d.) random multivariate and irregularly sampled pixels from a SITS, with their associated class labels, see Table 4.1 for an overview of the main notations. A pixel y is modeled by a random vector Y containing p random square integrable processes T → R p , where T = [t min , t max ] is the time interval where the SITS are observed:

Y : T → R p t → Y 1 (t), . . . , Y b (t), . . . , Y p (t) .
This property is denoted by Y ∈ L p 2 (T ). The associated class z is modeled by a discrete random variable Z with possible values in {1, . . . , C}. In the context of this work, p is the number of wavelengths and/or spectral indices, t min = 0 and t max = 365 (in day-of-year unit).

A (univariate) Gaussian process (GP) f is a stochastic process such that any finite-dimensional marginal follows a multivariate Gaussian distribution [START_REF] Williams | Gaussian processes for machine learning[END_REF]. It is specified by its mean function m and covariance function K:

m(t 1 ) = E[ f (t 1 )], K(t 1 , t 2 ) = E[( f (t 1 ) -m(t 1 ))( f (t 2 ) -m(t 2 ))],
with t 1 ∈ T , t 2 ∈ T two given times where the SITS are observed, and we note f ∼ GP(m, K). The covariance function K usually encodes a priori knowledge about the phenomena to be modeled. Figure 4.2 shows several realizations of a GP with a squared exponential covariance function (see Section 4.2.2).

In the following, a multivariate Gaussian process model is proposed to cope with the multivariate nature of SITS and their irregular temporal sampling. where b ∈ {1, . . . , p} is the spectral band and c ∈ {1, . . . , C} is the class. The first property states that each band of a pixel from a given class is a realization of a univariate GP with class and band specific mean and covariance functions. The second property is introduced for computational reasons, and is discussed in Section 4.6.2. Let y i be an observed pixel i at times {t i 1 , . . . , t i T i }. Its bth coordinate is represented by a vector of dimension T i :

Mixture of Independent Multivariate Gaussian Processes Model

y i,b = [Y i b (t i 1 ), . . . , Y i b (t i T i )] .
By definition, conditionally to Z = c, this vector follows a T i -variate Gaussian distribution

y i,b |Z i = c ∼ N T i µ i,b,c , Σ i b,c , (4.1)
where N T i is the Gaussian distribution on R T i with vector mean µ i,b,c = [m b,c (t i 1 ), . . . , m b,c (t i ), . . . , m b,c (t i T i )] and covariance matrix (Σ i b,c ) , = K b,c (t i , t i ). Finally, in view of the independence assumption on the p wavelengths, y i follows a product of p Gaussian densities

y i |Z i = c ∼ p b=1 N T i µ i,b,c , Σ i b,c . (4.2)

Mean and covariance functions

The unknown parameters of the model are the mean functions m b,c and the covariance operators K b,c , for each band b ∈ {1, . . . , p} and each class c ∈ {1, . . . , C}. To deal with irregularly sampled observations, these parameters should be defined (and estimated) for any time t in the interval T . To this end, projection techniques and parametric models are adopted respectively for the mean and the covariance functions. Parametric modeling of the latter ones are usual with GP, see for instance [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 4]. However, in our specific problem, using a projection method for the mean function allows to cope nicely with the irregular temporal sampling.

Mean function

Let {ϕ j } J j=1 be a subset of basis functions of L 2 (T ). Then, m b,c can be written as

m b,c (t) = J j=1 α b,c, j ϕ j (t) + ε b,c (t), t ∈ T ,
with α b,c, j the projection coefficient of m b,c on ϕ j and where ε b,c is an approximation error term. Any functional basis can be used, e.g., Fourier, exponential, splines, . . . , see for instance [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 5]. In case of the Fourier basis, α b,c, j represents the amplitude of the corresponding frequency in the Fourier expansion. In the following, the basis is assumed to be fixed while J, the number of basis functions, is an hyperparameter.

Introducing α b,c = [α b,c,1 , . . . , α b,c,J ] ∈ R J and B i the T i × J design matrix associated with pixel y i defined by (B i ) , j = ϕ j (t i ) with ∈ {1, . . . , T i }, the mean vector in (4.1) can be written as

µ i,b,c = B i α b,c . (4.3)

Covariance function

The covariance function is modeled using functions issued from the GP literature [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 4]. A typical one would be a squared exponential covariance function with an additional colored noise [START_REF] Le | Heteroscedastic Gaussian process regression[END_REF] covariance function:

K b,c (t, s) = γ 2 b,c exp        - (t -s) 2 2h 2 b,c        + σ 2 b,c (t)δ t,s . (4.4)
Introducing θ b,c = {γ 2 b,c , h 2 b,c , σ 2 b,c }, the covariance matrix in (4.1) is denoted in the following by

Σ i b,c = Σ i (θ b,c ). (4.5)
Parameters α b,c and θ b,c are estimated by maximizing the marginal log-likelihood, as explained in the next section.

Estimation

By plugging (4.3) and (4.5) in (4.2), it follows that

y i |Z i = c ∼ p b=1 N T i B i α b,c , Σ i (θ b,c ) . (4.6)
Algorithm 2: Estimation of model parameters.

Input : S,α 0 ,θ 0 Output: α, θ

i|Z i =c log Σ i (θ b,c ) + (y i,b -B i α b,c ) Σ i (θ b,c ) -1 (y i,b -B i α b,c ) + κ,
where κ is a constant independent of the parameters. Since independence is assumed between each spectral component, the optimization ends up with p × C independent optimization problems:

( αb,c , θb,c ) = arg min α b,c ,θ b,c b,c (α b,c , θ b,c ). (4.7)
Marginal likelihood optimization is common practice in GP regression, while for binary classification, Laplace approximation is usually employed (and more computationally demanding) [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 2 and 3]. Therefore, to the best of our knowledge, the use of marginal likelihood for multiclass classification is novel is this context. Each sub-problem (4.7) is solved by gradient descent as usually done in likelihood optimization [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 5]. The algorithm is based on an alternate optimization with respect to α and θ (see Algorithm 2). At iteration (k), the update rule for α is given by

α (k+1) b,c =         i|Z i =c B i Σ i (θ (k) b,c ) -1 B i         -1         i|Z i =c B i Σ i (θ (k) b,c ) -1 y i,b         . (4.8)
If the design matrix has no redundancies, it can be shown that the matrix in the left hand side of (4.8) is indeed non singular provided that the number of basis functions J is smaller or equal to the total number of unique observations in the training set. A proof is given in Section I of the supplementary materials. For θ, there is no close-form expression and a gradient step is required. The gradient is computed using the following partial derivative w.r.t. each coordinate θ m of θ:

∂ ∂θ m b,c = i|Z i =c tr         Σ i (θ (k) b,c ) -1 -β i β i ∂Σ i (θ (k) b,c ) ∂θ m         (4.9) with β i = Σ i (θ (k) b,c ) -1 (y i -B i α (k+1) b,c
). Such an optimization procedure applied on simulated data from As a final remark, the p × C optimization problems can be solved in parallel since parameters are not shared between classes and spectral components.

Numerical Complexity

Usual GPs have complexity that scales in O(n 3 ), making them unsuitable for very large scale problems. In contrast, the proposed method has a reduced complexity, that scales in O(n(T 3 M + J 3 )) where T M is the maximal length of observed time-series (T M ≥ T i , ∀i ∈ 1, . . . , n). The first term comes from the inversion of Σ i (θ b,c ) and the second term comes from the update rule (4.8). In our application, T M and J are much smaller than n, typically by several orders of magnitude. Furthermore, MIMGP allows for the classification of any irregularly new time-series as well as the reconstruction of the observed time-series using reconstruction techniques, as described in the following.

Classification and Reconstruction of Missing Values

Once the parameters are estimated, it is possible to classify a new pixel without any temporal resampling, as for training. In addition, MIMGP is also able to reconstruct missing values on any temporal scheme.

In the following, y j denotes a new SITS observed at T j times denoted by {t 1 , . . . , t T j } which may not have been observed in the training set.

Classification of a new time-series

The a posteriori probability P(Z j = c|y j ) to belong to a class c given y j is computed using the product of Gaussian densities given in (4.6) and Bayes' rule:

P(Z j = c|y j ) ∝ π c P(y j |Z j = c) ∝ π c p b=1 f T j y j,b ; B j α b,c , Σ(θ b,c ) , (4.10)
where f d (•; µ, Σ) is the d-variate Gaussian density with mean µ and covariance matrix Σ and π c is the prior probability P(Z = c) of the class c ∈ {1, . . . , C}.

In practice, π c is estimated by its empirical counterpart πc = n c /n where n c is the number of pixels assigned to class c in the set S. Parameters αb,c and θb,c are estimated thanks to Algorithm 2. The new time-series y j is then assigned to the class of maximum posterior probability (MAP rule):

ẑ j = max c P(Z j = c|y j ).

Time-series reconstruction

Time-series reconstruction is achieved using conditional expectations of Gaussian distributions. Two cases are considered: Either the class membership of the considered pixel is known, or the class membership is estimated using the posterior probability.

Reconstruction when class membership is known

Let us write Y j b (t * ) the unobserved value at wavelength b and time t * . The following reconstruction rule, based on the conditional expectation, is considered:

Ŷ j b,c (t * ) := E[Y j b (t * )|y j , Z j = c].
Using properties of Gaussian distributions [17, p.63] and replacing the unknown quantities by their estimated counterparts yield (see Appendix 4.9 for a proof):

Ŷ j b,c (t * ) = b * αb,c + k(t * , t j 1:T j | θb,c ) Σ j ( θb,c ) -1 (y j,b -B j αb,c ) (4.11)
and

V Ŷ j b,c (t * ) = K(t * , t * | θb,c ) -k(t * , t j 1:T j | θb,c ) Σ j ( θb,c ) -1 k(t * , t j 1:T j | θb,c ), (4.12)
where k(t * , t 1:

T j | θb,c ) = [K(t * , t 1 | θb,c ), . . . , K(t * , t T j | θb,c )] , b * = [ϕ 1 (t * ), .
. . , ϕ J (t * )] contains the evaluation of the basis functions at time t * and B j is the T j × J design matrix associated with y j . Equations (4.11) and (4.12) provide respectively the reconstructed value and the variance of the reconstruction. They are given here for a single time t * but similar formulas can be derived for multiple times t * 1 , . . . , t * T * . Interestingly, Equation (4.11) shows that the reconstruction is given by b * αb,c , the estimated mean of the GP, corrected by a value proportional to the error made at the acquisition time of pixel j. A similar remark holds for the variance: It can be interpreted as the estimated variance of the GP corrected by the variance of the process observed at the acquisition times of pixel j.

Reconstruction when class membership is unknown

The reconstruction of Ŷ j b (t * ) is done as the average of the previous reconstructions Ŷ j b,c (t * ) in each class (see (4.11)) weighted by the posterior probabilities estimated with (4.10):

Ŷ j b (t * ) = C c=1 P(Z j = c|y j ) Ŷ j b,c (t * ). (4.13)
A similar formula holds for the variance with an additional between-classes variance term:

V[ Ŷ j b (t * )] = C c=1 P(Z j = c|y j )V[ Ŷ j b,c (t * )] + C c=1 P(Z j = c|y j ) Ŷ j b,c (t * ) 2 -Ŷ j b (t * ) 2 ,
where Ŷ j b,c (t * ) and Ŷ j b (t * ) are the reconstructions at t * respectively when the class is known (4.11) and when the class is unknown (4.13). See Appendix 4.9 for a proof.

Sentinel-2 Satellite Image Time-Series Datasets

Three Sentinel-2 tiles of level 2A over the French metropolitan territory were downloaded from the Theia Land Data Center2 . All available acquisitions between January 2018 and December 2018 for the two orbits of satellites Sentinel-2A and 2B were used. Figure 4.5 shows the location of the three tiles. They correspond to different climatic regions [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF], with varying meteorological and topographical conditions.

Surface reflectance time-series were produced using the MAJA (Multi-sensor atmospheric correction and cloud screening-ATCOR Joint Algorithm) processing chain developed by the CNES-CESBIO and DLR [START_REF] Baetens | Validation of Copernicus Sentinel-2 Cloud Masks Obtained from MAJA, Sen2Cor, and FMask processors using reference cloud masks generated with a supervised active learning procedure[END_REF]. It involves orthorectification, atmospheric correction, clouds and shadows detection. Spectral bands at 10m/pixel and 20m/pixel were used, for a total of 10 spectral bands. Bands at 20m/pixel were spatially up-sampled to 10m/pixel using the Orfeo ToolBox3 [171]. Furthermore, the data from each tile were re-sampled and gap-filled (missing information due to clouds or shadows were reconstructed using linear interpolation) onto the same set of dates (every 5 days, starting from 2018-01-01 and ending 2018-12-27) as in [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. Figure 4.11 shows the raw data and the re-sampled data for one random pixel from each tile in the same class. Hence, for each pixel location, three temporal informations are available after the pre-processing:

1. The raw multispectral time-series with irregular acquisition dates; 2. The mask time-series, indicating for each acquisition date the presence/absence of clouds/shadows. In our model, the presence of clouds/shadows implies that the corresponding raw spectral values are considered as missing values;

3. The re-sampled multispectral time-series with regular acquisition dates.

The reference data arise from the work of Inglada et al. [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. They were extracted from freely available data source. Seventeen land cover classes were defined, ranging from artificial areas to vegetation and water bodies. Table 4.2 shows the exhaustive list of classes.

The reference data are provided as a set of spatial polygons overlapping the 3 tiles. Figure 4.7 shows an extract of these polygons. The training and validation set were constructed by stratifying pixels according to the polygons membership information: Pixels from one polygon fully belong to either the training or the validation set. Depending on the number of available referenced pixels per class, 10,000 (or less) pixels were extracted for the training and validation set, except for winter and summer crops, for which 30,000 and 40,000 pixels were extracted, respectively. Ten independent train/validation sets were generated for statistical validation. Table 4.2 shows the average number of pixels per class and per tile4 .

Experimental set-up

The model parameters are α b,c and θ b,c for each class c and each wavelength b. They are estimated using the training set of pixels. MIMGP has also some hyperparameters that are set before the learning step: The design matrix B with the associated number of basis functions J and the family of parametric kernels for the covariance operator. These settings are common to each class and each wavelength.

Functional bases

Four functional bases have been investigated: Two of them are local bases (non-zero only on subsets of [0, T ]), the Gaussian and BSplines bases, while the remaining ones are global bases, the well-known Fourier and Polynomial bases. An user-defined hyperparameter, denoted by J * , is used to select the number of basis functions, as explained in the next paragraph. The choice of a basis setting corresponds to an assumption on the temporal behavior of the time-series: For instance, a Fourier basis can represent a periodic signal in the time domain.

• Fourier: The number of basis functions is J = 2J * + 1 and

m(t) = α 0 + J * j=1 α j cos(2π j t T ) + α j+J * sin(2π( j + J) t T ) . • Polynomial: m(t) = α 0 + J * j=1 α j t j , with J = J * + 1. • Gaussian: m(t) = J * j=1 α j exp - (t-t j ) 2 d 2 j
, with J = J * . The t j , j ∈ {1, . . . , J * } can be equidistant in T or chosen as quantiles of the distribution of the clean dates. The hyperparameter d j is set such that d 2 j = 8|t j+1 -t j | to ensure a sufficient overlap between two consecutive local exponential functions.

• Cubic splines: m(t) = J * j=0 α j S j (t), with J = J * + 1. S j is the j th bicubic spline on [0, T ]. The knots t j are chosen as in the Gaussian case.

The invertibility condition associated with (4.8) implies that J should be smaller than 92, 100 and 103 for tiles T31TCJ, T31TDN and T31TGK, respectively.

Covariance function

Any convex combination of positive semi-definite kernels is a valid covariance function [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 4]. In this work, the squared exponential covariance function added with a colored noise covariance function, as given in (4.4), is used.

The parameters

θ b,c = {γ 2 b,c , h 2 b,c , σ 2 b,c }, ∀{b, c} ∈ [1, p] × [1,
C] are learned as described in Algorithm 2. The estimated covariance parameters provide some insights about the observed processes. In particular, the lengthscale h b,c is related to the temporal behavior of the reflectance. For a given band b and conditionally to class c, the longer h b,c grows, the more correlated two distant dates are.

Supervised classification

The classification accuracy is assessed by the F 1 score, computed as the harmonic mean of the precision and recall for each class. We report the "mean F 1 score", which is the mean of F 1 scores computed on all the classes of the dataset.

First, the influence of the basis functions on the classification score is investigated. Then, comparison with other classifiers is reported and discussed. Convergence and model parameter analyses are provided in Section II of the supplementary materials. The results are reported for the tile T31TCJ.

Influence of the basis functions

The influence of the basis functions and its dimension are investigated. Figure 4.8 represents the averaged "mean F 1 scores" obtained on the 10 independent train/validation data sets for the tile T31TCJ. From the figure, it appears that the dimension of the basis should be larger than 9 for this tile to allow for a reasonable classification score, independently of the basis itself. Furthermore, MIMGP is robust w.r.t. basis selection since similar accuracy are obtained for all bases. However, numerical instability has been observed for very large number of basis functions, in particular for the exponential basis. For instance, when the number of functions is greater than 50, consecutive local exponential bases overlap too much and the invertibility conditions are violated (see Section I of the supplementary materials). In practice, the Fourier basis is the most stable and was used to compare our method to other classifiers.

Comparison with other classifiers

The performances of MIMGP are compared to three other methods: Quadratic Discriminant Analysis (QDA) which involves a similar Gaussian assumption on regularly sampled data, a linear Support Vector Machine (SVM) fitted with a Stochastic Gradient Descent (SGD) [START_REF] Zhang | Solving Large Scale Linear Prediction Problems Using Stochastic Gradient Descent Algorithms[END_REF], and Random Forests (RF) [START_REF] Breiman | Random Forests[END_REF] which have shown state of the art results in large scale pixel-wise classification of SITS [START_REF] Derksen | Geometry Aware Evaluation of Handcrafted Superpixel-Based Features and Convolutional Neural Networks for Land Cover Mapping Using Satellite Imagery[END_REF]. QDA, SVM and RF methods take the re-sampled multispectral time-series as input. The 10 spectral bands are stacked together to let these three methods account for dependence between spectral bands.

RF is used with 100 trees and a maximum depth of tree set to 25. MIMGP model is implemented with the Fourier basis and a RBF kernel with an additive colored noise which represents a total of 32 parameters for θ. Table 4.3 summarizes the results for the four methods 5 . We can see that some significant improvements are achieved compared to the standard QDA classifier, even though MIMGP imposes independence between wavelengths. The QDA classifier suffers from numerical instability due to the high number of spectro-temporal features. However, MIMGP is not as accurate as the RF and the SVM classifiers applied on the reconstructed time-series.

Figure 4.9 shows three extracts of the classification map for tile T31TCJ 6 . Visually, significant differences can be observed both in terms of homogeneity and in terms of pixel-wise prediction. For instance, for the airport area, SVM predicts wrongly a lot of natural grasslands while MIMGP and RF predict correctly meadow for vegetation between runways. Overall, MIMGP and SVM seem to predict more homogeneous maps than RF. This can be clearly seen in the second row, where the city center classification maps obtained with RF contain salt and pepper noise.

Results provided in this section show that MIMGP is not at the level of state-of-the art methods such as RF or SVM. One critical point is the spectral independence assumption that is not true in practice and used here for computational purposes. Removing this constraint would help to improve the classification accuracy.

Yet, the computation of MIMGP does not require any temporal resampling processes and directly handles the raw time-series in the training and prediction steps. In comparison with standard GPs for classification, its complexity is linear w.r.t to the number of samples and cubic w.r.t to the number of temporal acquisitions, thus allowing large scale processing. Furthermore, it can also be used to reconstruct pixel time-series, as discussed in the next section.

Time-series reconstruction

Following Section 4.3.2, MIMGP is able to reconstruct missing values conditionally to a known, or predicted, class. Class mean functions and covariance operators are also obtained after the training step (for clarity, mean functions and covariance operators are provided in Section III of the supplementary materials). In the following, the quality of the reconstruction is compared with a state-of-the-art method: the Whittaker smoother [START_REF] Eilers | A Perfect Smoother[END_REF][START_REF] Vuolo | Smoothing and gap-filling of high resolution multi-spectral time series: Example of Landsat data[END_REF].

To assess the quality of the reconstruction, the validation set was used to ensure that reconstructed pixels were not seen during the training phase. We randomly remove one clear date from each pixel and reconstruct it with the trained MIMGP model and with the Whittaker smoother for which the regularization parameter is set using the V-curve technique for each individual pixel [START_REF] Frasso | L-and v-curves for optimal smoothing[END_REF][START_REF] Vuolo | Smoothing and gap-filling of high resolution multi-spectral time series: Example of Landsat data[END_REF]. The quality of the reconstruction is estimated using the normalized Mean Absolute Error (nMAE) computed as:

nMAE = n v i=1 |Y i b (t * ) -Ŷi b (t * )| n v i=1 |Y i b (t * ) -Y b | with Y b = 1 n v n v i=1 Y i b (t * )
, n v is the number of validation pixels. Table 4.4 summarizes the results obtained on the three tiles and for the 10 independent data sets. The scatter plots associated with nMAE results are provided in the supplementary materials.

From the table, MIGMP outperforms the Whittaker smoother in terms of nMAE for each band for tiles T31TDN and T31TGK. For tile T31TCJ, Whittaker is better for 5 bands. It is important to note that MIMGP provides better nMAE and does not need any additional fitting, apart the training step, in contrast to the Whittaker smoother that requires to be fitted for each pixel, as well its regularization parameter.

Figure 4.10 shows a part of the tile T31TGK reconstructed for the band 4 (red) on August 6, 2018. For this simulation, all the clear dates are used to reconstruct the time series for both methods. Interestingly, for the Whittaker smoother, the clouds mask is visible indicating a non continuous reconstruction of the dynamic. This artifact is not visible for the proposed method. Figure 4.11 illustrates the red wavelength (R) for two pixels of the same class (summer crops). They were selected to enlighten the robustness of the proposed method to an inaccurate cloud mask file. The first pixel has an undetected cloud for the third temporal acquisition (blue point in the upper Figure 4.11). A clear drop of the reflectance can be seen in the gap-filled infra-red (dashed red line). The reconstructed reflectance for MIMGP (black line) and for Whittaker smoother (black dotted line) do not exhibit such a drop in the reflectance. Similar comments can be done for the second pixel of the lower Figure 4.11 which has a undetected saturation (pixel value equal to zero). Other wavelengths are reported in the supplementary materials.

Conclusion

This work introduces a novel approach to jointly classify and reconstruct irregularly sampled multidimensional time-series. The proposed model, namely MIMGP, only involves a small number of parameters and is scalable to large datasets. The performances of the method were illustrated on a full-year of Sentinel-2 satellite image time-series dataset from 2018 involving a high number of temporal acquisitions. A by-product of MIMGP is also to infer a confidence interval on the reconstruction. The reconstruction has shown a good behavior on noisy pixels, while the model does not compete with state of the art classifiers such as Random Forest or SVM. One limitation comes from the independence hypothesis between spectral wavelengths, which is not true in practice. Therefore, our future work will be dedicated to the 4.9. Appendix -Time-series reconstruction 57 definition of multivariate Gaussian Processes able to take into account the correlation between wavelengths.

Another work in progress concerns the estimation of the mean function. MIMGP is able to handle the irregular sampling caused by multiple tiles, but the mean function is assumed to be spatially stationary. Rather than using a linear combination of fixed basis functions, a more complex approximating function depending on the geographical location could be used.

Finally, this work can be extended in a number of directions that could be of interest in the remote sensing field. A first extension would be the use of the MIMGP model to perform classification in the unsupervised case thanks to an EM-like algorithm. A second direction of future research would be the joint use of both Sentinel-2 and Sentinel-1 time-series thanks to the MIMGP model which allows for arbitrary temporal sampling in each spectral band.

Appendix -Time-series reconstruction

Let us first consider the case where the class of the missing value is known to be c. In such a situation, the reconstruction can be achieved using the conditional expectation:

Ŷ j b,c (t * ) := E[Y j b (t * )|y j , Z j = c]. The independence assumption A2 then yields E[Y j b (t * )|y j , Z j = c] = E[Y j b (t * )|y j,b , Z j = c]. Besides, conditionally on Z i = c, one has        Y j b (t * ) y j,b        ∼ N T j +1 b * α b,c B j α b,c ,        K(t * , t * |θ b,c ) k(t * , t j 1:T j |θ b,c ) k(t * , t j 1:T j |θ b,c ) Σ j (θ b,c )               ,
where k(t * , t 1:

T j |θ b,c ) = [K(t * , t 1 |θ b,c ), . . . , K(t * , t T j |θ b,c )] , b * = [ϕ 1 (t * ), . . . , ϕ J (t * )
] and B j is the T j × J design matrix associated with the time-series y j such that (B j ) ,k = ϕ k (t j ) for all ( , k) ∈ {1, . . . , T j } × {1, . . . , J}. From classical properties of conditional Gaussian distributions (see for instance [17, p.63]), it follows that

E(Y j b (t * )|y j,b ; Z j = c) = b * α b,c , +k(t * , t j 1:T j |θ b,c )Σ j (θ b,c ) -1 (y i,b -B j α b,c ), V(Y j b (t * )|y j,b ; Z j = c) = K(t * , t * |θ b,c ) -k(θ b,c |t * , t j 1:T j )Σ j (θ b,c ) -1 k(t * , t j 1:T j |tθ b,c
) . Finally, replacing the unknown quantities by their estimated counterparts yields the desired results.

Second, when the class is unknown, the reconstruction rule is given by

Ŷ j b (t * ) = E[Y j b (t * )|y j ]. Since E[Y j b (t * )|y j ] = C c=1 P(Z j = c|y j )E[Y j b (t * )|y j , Z j = c],
it straightforwardly follows that

Ŷ j b (t * ) = C c=1 P(Z j = c|y j ) Ŷ j b,c (t * ),
where the reconstructed value Ŷ j b,c (t * ) is provided in (4.11) while the posterior probabilities are given by (4.10). The variance is obtained by a similar calculation:

V(Y j b (t * )) = E[Y j b (t * ) 2 |y j ] -E[Y j b (t * )|y j ] 2 , = C c=1 P(Z j = c|y j )E[Y j b (t * ) 2 |y j , Z j = c] -Ŷ j b (t * ) 2 .
Besides, remarking that

E[Y j b (t * )|y j , Z j = c] = Ŷ j b,c (t * ), E[Y j b (t * ) 2 |y j , Z j = c] = V(Y j b,c (t * )|y j , Z j = c) + E[Y j b,c (t * )|y j , Z j = c] 2 , it follows that V( Ŷ j b (t * )) = C c=1 P(Z j = c|y j )V( Ŷ j b,c (t * )|y j , Z j = c) + C c=1 P(Z j = c|y j ) Ŷ j b,c (t * ) 2 -Ŷ j b (t * ) 2 ,
and the result is proved. Let us highlight that these derivations were conducted when a single time t * is considered. Similar calculations can be achieved when reconstructed simultaneously several values. In such a case, this estimation procedure provides the covariance matrix of the reconstructed values at each time.

Introduction

Satellite images availability has exponentially grown in the last decade. Thanks to free data access policy, optical satellite image time-series (SITS) such as Landsat or Sentinel-2, offer an unique opportunity to monitor the state and evolution of our living planet. Therefore, SITS have found many applications in ecological monitoring [START_REF] Li | Estimating Plant Traits of Alpine Grasslands on the Qinghai-Tibetan Plateau Using Remote Sensing[END_REF][START_REF] Fauvel | Prediction of plant diversity in grasslands using Sentinel-1 and -2 satellite image time series[END_REF], meteorology [START_REF] Liu | A spatio-temporal modeling framework for weather radar image data in tropical Southeast Asia[END_REF][START_REF] Bertolacci | Climate inference on daily rainfall across the Australian continent, 1876-2015[END_REF] or agricultural system mapping [START_REF] Useya | Comparative Performance Evaluation of Pixel-Level and Decision-Level Data Fusion of Landsat 8 OLI, Landsat 7 ETM+ and Sentinel-2 MSI for Crop Ensemble Classification[END_REF][START_REF] Moeini Rad | Developing an Automatic Phenology-Based Algorithm for Rice Detection Using Sentinel-2 Time-Series Data[END_REF][START_REF] Feng | Crop Type Identification and Mapping Using Machine Learning Algorithms and Sentinel-2 Time Series Data[END_REF], among others. SITS are characterized by their spatial and spectral resolutions, and their revisit cycle. The spatial resolution corresponds to the size of a pixel on the ground, e.g., a square of 10 meters while the spectral resolution is related to the number of wavelengths collected by the sensor, ranging typically in the visible and near infra-red part of the spectrum [START_REF] Manolakis | Hyperspectral Imaging Remote Sensing: Physics, Sensors, and Algorithms[END_REF]. The revisit cycle stands for the time between two acquisitions over the same location: SITS have constant and short (e.g. few days) revisit time. Hence, for a given temporal period, a pixel is the collection of spectral measurements made at different times over the same location.

These properties lead to an unprecedented amount of numerical data, for which statistical methods are used to extract meaningful information such as land cover, crops yields . . . For a pixel-wise based analysis, the predictor variables are multivariate time-series and the output variables represent the information to be extracted. While spatial independence is usually assumed [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF], temporal and spectral correlations are commonly taken into account in statistical models [START_REF] Lopes | Spectro-temporal heterogeneity measures from dense high spatial resolution satellite image time series: Application to grassland species diversity estimation[END_REF].

However, external random meteorological factors interfere with the availability of the acquired data at the pixel scale. Indeed, as displayed in Fig. 5.1, shadows and clouds result in missing data in the time-series. Furthermore, orbital trajectory generates an irregular temporal sampling: Even though the acquisition scheme is regular, acquisition days are different for pixels located at different places [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. As such, each pixel of the SITS has its own size in the temporal domain: Fig. 5.2 illustrates the irregular temporal sampling on the data under consideration in this paper.

Specific models are thus required to properly analyze such time-series, as described in Section 5.2. Conventional approaches usually start by resampling the data onto a common temporal grid. In this work, we aim at analyzing irregularly sampled multidimensional SITS without any temporal resampling. In particular, the supervised pixel classification task is considered, i.e. the assignment of each pixel of the time-series to a predefined class.

To this end, a mixture of multivariate Gaussian Processes is proposed. A linear dependence model is assumed between the spectral variables leading to a separable covariance function in time and spectral domains. The resulting model provides statistical information on the underlying process for each class (mean and covariance functions) and scales linearly w.r.t. the number of samples. The remainder of this paper is organized as follows. Section 5.2 reviews the state-of-the-art on classification with missing data and Gaussian processes. The statistical model is introduced in Section 5.3 while inference aspects are discussed in Section 5.4 including the estimation of the model parameters, the supervised classification, and the imputation of missing values. These statistical procedures are validated on simulated data in Section 5.5. Section 5.6 is dedicated to the application of our methodology to the classification of SITS from a Sentinel-2 data-set. Section 5.7 concludes with a discussion on possible extensions of this work.

Related Work

This section briefly reviews state-of-the-art methods for model-based classification, classification dealing with missing values and classification with Gaussian processes.

Supervised model-based classification

Supervised model-based classification (also referred to as model-based discriminant analysis) starts from a training set of n independent realizations from a random pair (Y, Z) ∈ E × {1, . . . , C} and assumes that the conditional distribution of y|Z = c belongs to some parametric family: p(y|Z = c) = p c (y; θ c ), for all c ∈ {1, . . . , C} and y ∈ E, where E is an arbitrary space. Letting π c = P(Z = c), the marginal distribution of Y can be written as a finite mixture

p(y) = C c=1 π c p c (y; θ c ),
whose parameters can be estimated by the maximum likelihood principle. A non-labeled observation can then be classified thanks to the Maximum a posteriori (MAP) criteria:

ĉ = arg max c∈{1,...,C} p(Z = c|y) = arg max c∈{1,...,C} π c p c (y; θ c ),
thanks to Bayes' rule. When E = R q , the multivariate Gaussian distribution is often adopted for p c (y; θ c ) and gives rise to the well-known Quadratic discriminant analysis (QDA) method. We refer to [84, Section 4.3] for a discussion on the advantages and drawbacks of QDA and for possible extensions. Recent studies extend the model-based classification framework to non-Gaussian distributions such as the skew-normal distributions [START_REF] Theodossiou | Financial Data and the Skewed Generalized t Distribution[END_REF][START_REF] Chamroukhi | Skew t mixture of experts[END_REF] to deal with asymmetric data, or t-distributions [4,[START_REF] Murray | A mixture of SDB skew-t factor analyzers[END_REF] to deal with outliers. We refer to [START_REF] Bouveyron | Model-Based Clustering and Classification for Data Science: With Applications in R[END_REF]Chapter 9] for an in-depth review. The case E = R q also encompasses the situation of discretized time-series on a common grid. Specific models can be then defined, as in [START_REF] Povinelli | Time series classification using gaussian mixture models of reconstructed phase spaces[END_REF] for temporal signatures.

If E is discrete, including for example the case of categorical data, extensions focus on the multinomial [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF] or the Dirichlet [START_REF] Bouguila | Novel Mixtures Based on the Dirichlet Distribution: Application to Data and Image Classification[END_REF] distributions. In the case of ordinal data, other extensions are proposed using a dedicated model of the process generating data [START_REF] Biernacki | Model-based clustering of multivariate ordinal data relying on a stochastic binary search algorithm[END_REF]. Finally, when E is more complex, e.g. infinite dimensional, nonparametric techniques are used. Kernel methods are probably the most popular non-parametric techniques in this situation [START_REF] Hofmann | Kernel methods in machine learning[END_REF]. Recall that a kernel is a positive definite function that corresponds to a dot product in a feature space. It allows for the construction of non-linear and non-parametric classifiers on E without computing explicitly the feature space. Kernels can be defined, for instance, on strings [START_REF] Lodhi | Text classification using string kernels[END_REF], graphs [START_REF] Kriege | A survey on graph kernels[END_REF], vector-valued functions [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF][START_REF] Flaxman | Scalable high-resolution forecasting of sparse spatiotemporal events with kernel methods: A winning solution to the NIJ "Real-Time Crime Forecasting Challenge[END_REF], or combinations of several data types [START_REF] Bouveyron | Kernel discriminant analysis and clustering with parsimonious gaussian process models[END_REF].

Classification with missing data

When dealing with remote sensing data, i.e. spatial-spectro temporal data such as SITS, handling missing values [START_REF] Shen | Missing Information Reconstruction of Remote Sensing Data: A Technical Review[END_REF] is a recurrent problem. Classification dealing with missing data occurs when some inputs in the training set are incomplete, i.e. the number of available coordinates in y can be different from one sample to another, see [START_REF] Schafer | Analysis of Incomplete Multivariate Data[END_REF][START_REF] García-Laencina | Pattern classification with missing data: A review[END_REF][START_REF] Bagnall | The great time series classification bake off: A review and experimental evaluation of recent algorithmic advances[END_REF] for reviews.

Three main approaches can be found in the literature. A first solution is to impute missing values before the classification itself. The pre-processing gives rise to a training set with observations re-sampled on a common grid that can be considered as vectors in a finite space E = R q , opening the door to classical model-based classification methods. We refer to [START_REF] Lin | Missing value imputation: A review and analysis of the literature (2006-2017)[END_REF] for a review on imputation techniques. For SITS, [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF] used such two-stage approaches on Sentinel-2 SITS where a linear interpolation was applied before performing the classification with a Random Forests classifier. Yet, by applying imputation techniques without any connection to the actual processing, propagated errors from the interpolation may degrade the results.

Alternative solutions are based on functional data analysis [START_REF] Ramsay | Functional Data Analysis[END_REF]. Each observation is interpreted as a sample from a random function. As such, it can be approximated by an expansion on some basis functions. The statistical analysis is then performed on the random vectors of coefficients, see [START_REF] Schmutz | Clustering multivariate functional data in group-specific functional subspaces[END_REF] for an application to clustering. Nonparametric smoothing techniques may also adopted, see [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF]Chapter 8] for an overview.

Finally, purely non-parametric methods can also be implemented by defining an appropriate dissimilarity measure between samples of varying size. In the context of time-series, Dynamic time warping (DTW) [START_REF] Chouakria | Adaptive dissimilarity index for measuring time series proximity[END_REF] is one of the most popular algorithms. It computes an optimal match between two vectors with different lengths. This map defines a dissimilarity that can be used for comparison in order to cluster samples into multiple groups.

Classification with Gaussian processes

A recent approach for supervised classification is based on the use of Gaussian processes (GPs) in a Bayesian framework. More specifically, Gaussian processes are used as prior distributions on the regression function linking the label Z to the explanatory variable X. In the binary classification case, the conditional Bernoulli distribution of Z is defined through a logit transformation: logit(p(Z = 1|X = x)) := f (x) where f (x) is a centered Gaussian process. The considered prior Gaussian process is, most of the time, one-dimensional. Extensions to the multi-dimensional case include the so-called multi-tasks or multi-outputs GP models, see [START_REF] Bonilla | Multi-task Gaussian Process Prediction[END_REF][START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF]. Finally, some recent works focus on non Gaussian processes such as Student-t processes which have gain attention over the past years [START_REF] Shah | Student-t Processes as Alternatives to Gaussian Processes[END_REF][START_REF] Chen | Multivariate Gaussian and Student-t process regression for multioutput prediction[END_REF].

The discrete nature of Z makes the exact inference of model parameters infeasible. To overcome this difficulty, several techniques have been proposed, including the Laplace approximation, or through the expectationpropagation algorithm [START_REF] Nickisch | Approximations for binary gaussian process classification[END_REF]. Such approaches rely on the inversion of a n × n covariance matrix and thus scale in O(n 3 ) which makes the inference computationally demanding for large data sets. Scalable GPs were proposed to overcome this vexing effect, using for instance variational inference as in [START_REF] Hensman | Scalable Variational Gaussian Process Classification[END_REF]. We refer to [START_REF] Liu | When Gaussian process meets big data: A review of scalable GPs[END_REF] for a review on this topic. In the next Section, we define a mixture of multivariate Gaussian processes which can be used for classification or imputation tasks without resort to approximate inference techniques.

Mixture of Multivariate Gaussian processes

The mixture of multivariate Gaussian processes model is introduced in Paragraph 5.3.1 and some associated properties are derived in Paragraph 5.3.2.

Model

Let T be a compact subset of R, throughout this document, we denote by GP 1 (0, K) a continuous univariate centered Gaussian process on T with covariance function K : T 2 → R. Recall that, by definition, W ∼ GP 1 (0, K) implies that, for all (t 1 , . . . , t q ) ∈ T q , the random vector (W(t 1 ), . . . , W(t q )) follows a multivariate centered Gaussian distribution N q (0, Σ) such that Σ i, j = K(t i , t j ), see for instance [START_REF] Williams | Gaussian processes for machine learning[END_REF].

For all p > 0, let us similarly denote by IGP p (0, K) a p-dimensional, independent, centered Gaussian process defined as

W = (W 1 , . . . , W p ) ∼ IGP p (0, K) if and only if        W b ∼ GP 1 (0, K), ∀b ∈ {1, . . . , p}, W b W b , ∀b b ∈ {1, . . . , p} 2 ,
where stands for independence. The above defined multivariate Gaussian processes are the building blocks to define more general multivariate Gaussian processes denoted by MGP p (m, K, A) where m : T → R p is the mean function, K : T 2 → R is the covariance operator and A a non-singular p × p matrix:

Y ∼ MGP p (m, K, A) if and only if Y = AW + m with W ∼ IGP p (0, K). (5.1)
Let us remark that model (5.1) is not identifiable without additional constraints. Indeed, MGP p (m, K, A) and MGP p (m, λK, A/ √ λ) yield the same process for all λ > 0. This issue is discussed in further details in Section 5.4, see also the next paragraph for some basic properties of multivariate Gaussian processes defined in (5.1).

Introduce Z a discrete random variable taking its values in c ∈ {1, . . . , C} with π c = P(Z = c). The mixture of multivariate Gaussian processes (M2GP) is defined by:

Conditionally to Z = c, Y ∼ MGP p (m c , K c , A c ), (5.2)
where m c : T → R p , K c : T 2 → R and A c is a non-singular p × p matrix, for all c ∈ {1, . . . , C}. In the context of SITS classification, Y represents the (unobserved) multidimensional process and p denotes the number of spectral bands. The particular case A c = I p yields a mixture of independent Gaussian processes (MIGP) whose applications to classification have been investigated in [START_REF] Constantin | Joint supervised classification and reconstruction of irregularly sampled satellite image times series[END_REF]. Let us also note that multivariate Gaussian processes have already been used in the machine learning community, without formal definition though, see for instance the so-called multi-task Gaussian process [START_REF] Bonilla | Multi-task Gaussian Process Prediction[END_REF] or the multivariate Gaussian process regression [START_REF] Chen | Multivariate Gaussian and Student-t process regression for multioutput prediction[END_REF].

First properties

Let C and D be two matrices of size m × n and p × q respectively. Recall that the Kronecker product C ⊗ D is the mp × nq matrix such that Keeping these definitions in mind, the matrix-variate normal distribution MN p,q [START_REF] Dawid | Some matrix-variate distribution theory: Notational considerations and a Bayesian application[END_REF][START_REF] Srivastava | Models with a kronecker product covariance structure: Estimation and testing[END_REF] is defined for all p × q random matrix Y as:

C ⊗ D =             c 11 D . . . c n1 D . . . . . . . . . c m1 D . . . c mn D             and vec(C) ∈ R mn
Y ∼ MN p,q (M, Σ, Λ) if and only if vec(Y ) ∼ N pq (vec(M), Σ ⊗ Λ), (5.3)
where M is a p × q matrix, Σ and Λ are symmetric positive definite matrices of size q × q and p × p respectively. We refer to [START_REF] Dawid | Some matrix-variate distribution theory: Notational considerations and a Bayesian application[END_REF] for an early definition of the matrix-variate normal distribution (as well as some of its derivatives), to [START_REF] Gupta | Matrix Variate Distributions[END_REF] for a general account on matrix-variate distributions and to [1] for an application to missing data imputation. The associated density function is defined for all p × q matrix y by

p(y) = (2π) -pq/2 det(Σ) -p/2 det(Λ) -q/2 exp - 1 2 tr Λ -1 (y -M)Σ -1 (y -M) , (5.4)
where tr(•) denotes the trace operator. The next Proposition establishes that the finite sized marginals of the multivariate Gaussian process (5.1) can be interpreted as random matrices from a matrix-variate normal distribution.

Proposition 5.1. Let Y ∼ MGP p (m, K, A) and introduce Y the p×q random matrix defined as Y = (Y(t 1 ), . . . , Y(t q )) where (t 1 , . . . , t q ) ∈ T q . Then, Y ∼ MN p,q (M, Σ, AA ), (

where M = (m(t 1 ), . . . , m(t q )) and Σ is the covariance matrix defined by Σ k, = K(t k , t ) for all (k, ) ∈ {1, . . . , q} 2 . Equivalently, vec(Y ) ∼ N pq (µ, Σ ⊗ AA ), with µ = vec(M).

In the SITS framework, Y represents the observed q-dimensional SITS which is a discretized version of Y at q timestamps. An illustration is provided in Fig. 5.3 where T = [0, 1] and p = q = 10. Only the first two coordinates are represented for lack of space reasons. Let •, • denote the Euclidean scalar product on R p and • be the associated norm. For all non zero vectors (u, v) ∈ R p × R p , we also introduce cos(u, v) = u, v /( u v ). As a direct consequence of the covariance structure in (5.5), the correlation ρ between the elements of the random matrix Y can be derived: (with a b the bth line of A) and is thus independent of j ∈ {1, . . . , q}.

(ii) For all ( j, j ) ∈ {1, . . . , q} 2 and b ∈ {1, . . . , p}, one has

ρ(Y b, j , Y b, j ) = Σ j, j Σ j, j Σ j , j , (5.6)
and is thus independent of b ∈ {1, . . . , p}.

Inference

65

It appears that A tunes the dependence between the lines of Y (i.e. the spectral bands in the SITS context) while Σ drives the dependence between the columns (i.e. the acquisition times of the SITS).

A likelihood ratio test is introduced in [START_REF] Lu | The likelihood ratio test for a separable covariance matrix[END_REF] to check whether the separability of the covariance (5.5) is adapted to the data in hand. However, this test has not been extended to irregularly sampled time-series. Let us also mention that, in [START_REF] Mahanta | Regularized LDA based on separable scatter matrices for classification of spatio-spectral EEG patterns[END_REF], the same Kronecker product model is used to regularize the estimation of the covariance matrix in high dimension.

Inference

This section addresses several inference aspects associated with the M2GP model. Consider {(Y 1 , Z 1 ), . . . , (Y n , Z n )} a set of n random pairs identically distributed from the M2GP model. Clearly, π c can be estimated by its empirical counterpart πc = n c /n where n c = n i=1 I{Z i = c} is the number of samples in class c (and I{•} is the indicator function). Besides, from (5.2), Y i ∼ MGP p (m c , K c , A c ) conditionally to Z i = c, for all i ∈ {1, . . . , n}. The unknown quantities to be estimated are m c : T → R p , K c : T 2 → R and the matrix A c . The use of parametric models for mean and covariance functions is discussed in Subsection 5.4.1 and the Maximum likelihood estimation (MLE) of all resulting parameters is presented in Subsection 5.4.2. The associated classification method based on the Maximum a posteriori (MAP) rule and the imputation of missing values are described in Subsection 5.4.3 and Subsection 5.4.4 respectively.

Parametric mean and covariance functions

Let J > 0 and introduce {ϕ 1 , . . . , ϕ J } a subset of J basis functions of L 2 (T ). For all b ∈ {1, . . . , p}, the bth coordinate

(m c (t)) b of m c (t) is expanded as (m c (t)) b = J j=1 α c,b, j ϕ j (t), (5.7)
with t ∈ T , and where a c,b, j is the projection coefficient of (m c (•)) b on ϕ j (•). Denoting by α c the p × J matrix defined by:

α c =                     α c,1,1 α c,1,2 . . . α c,1,J α c,2,1 . . . . . . α c,2,J . . . . . . . . . . . . α c,p,1 . . . . . . α c,p,J                    
and letting b : t ∈ T → (ϕ 1 (t), . . . , ϕ J (t)) ∈ R J , then (5.7) can rewritten matricially as m c (t) = α c b(t).

The covariance operator K c is assumed to belong to a family of symmetric positive-definite kernels, [START_REF] Williams | Gaussian processes for machine learning[END_REF]Chapter 4]. A typical kernel is the squared exponential kernel (also known as Gaussian or RBF kernel) with an additive white noise covariance function:

K c (t, t |θ c ) = γ 2 c exp - (t -t ) 2 2h 2 c + σ 2 c I{t = t }, (5.8)
where (t, t ) ∈ T 2 .The parameters are collected in θ c with, in this case, θ c = {γ c , h c , σ c }.

Maximum likelihood estimation

Assume each multivariate Gaussian process Y i is observed on its own finite grid of distinct q i timestamps (t i 1 , . . . , t i q i ) ∈ R q i and note Y i, = (Y i (t i 1 ), . . . , Y i (t i q i )) the associated p × q i random matrix. Let us stress that this formalism naturally allows to deal with irregularly sampled SITS since the size of Y i, may depend on i. From Proposition 5.1, one has

Conditionally to

Z i = c, Y i, ∼ MN p,q i (α c B i , Σ c,i (θ c ), A c A c ), (5.9)
where the covariance matrix Σ c,i (θ c ) is defined for all ( j, j ) ∈ {1, . . . , q i } 2 by Σ c,i (θ c ) j, j = K c (t i j , t i j |θ c ) and B i = (b(t i 1 ), . . . , b(t i q i )) is a J × q i design matrix. Parameters {α c , θ c , A c } are estimated by minimizing the negative loglikelihood given hereafter.

Lemma 5.1. The negative log-likelihood associated with (5.9) can be expanded as

L = 1 2 C c=1 c (α c , θ c , A c A c ),
(up to an additive constant) where, for all c ∈ {1, . . . , C},

c (α c , θ c , A c A c ) = Q c log det(A c A c ) + p i|Z i =c log det Σ c,i (θ c ) + tr         i|Z i =c (Y i, -α c B i ){Σ c,i (θ c )} -1 (Y i, -α c B i ) {A c A c } -1         , (5.10
)

with Q c = i|Z i =c q i .
It appears that the likelihood only involves the product of matrices A c A c and not the matrix A c itself. This is a direct consequence of (5.5): The matrix-variate normal distribution of the sampled process Y only depends on the above product. The parameters of interest are thus α c , θ c and A c A c and the MLE is obtained by solving C independent optimization problems:

( αc , θc , A c A c ) = arg min α c ,θ c ,A c A c c (α c , θ c , A c A c ), (5.11)
for all c ∈ {1, . . . , C}. The solution is partially explicit as explained in the next Proposition.

Proposition 5.2. Let c ∈ {1, . . . , C}.

(i) Solutions of (5.11) satisfy the following two properties. Given θc , one has:

αc =         i|Z i =c Y i, {Σ c,i ( θc )} -1 (B i )                 i|Z i =c B i {Σ c,i ( θc )} -1 (B i )         -1 , (5.12
)

A c A c = 1 Q c i|Z i =c (Y i, -αc B i ){Σ c,i ( θc )} -1 (Y i, -αc B i ) . (5.13)
(ii) The partial derivative of (5.10), w.r.t. the kth coordinate of θ c is given by:

∂ c (α c , θ c , A c A c ) ∂(θ c ) k = i|Z i =c tr p{Σ c,i (θ c )} -1 -β c,i (α c , θ c ) {A c A c } -1 β c,i (α c , θ c ) ∂Σ c,i ∂(θ c ) k (θ c ) , (5.14)
where

β c,i (α c , θ c ) = (Y i, -α c B i ){Σ c,i (θ c )} -1 .
In practice, the computation of the MLE is achieved thanks to an iterative procedure based on (5.12)-(5.14), described in Algorithm 3 and discussed in Paragraph 5.4.5.

Supervised classification

Starting from a training set from the M2GP model {(Y 1, , Z 1 ), . . . , (Y n, , Z n )} with Y i, = (Y i (t i 1 ), . . . , Y i (t i q i )) , our goal is to assign a label c ∈ {1, . . . , C} to a new p × q random matrix Ỹ = (Y( t1 ), . . . , Y( tq )) . We focus on the MAP rule which consists in maximizing w.r.t. c the posterior probability

P(Z = c| Ỹ ) ∝ π c p( Ỹ |Z = c),
where p( Ỹ |Z = c) is matrix-variate normal density defined as

p( Ỹ |Z = c) = (2π) -pq/2 det Σc (θ c ) -p/2 det A c A c -q/2 × exp - 1 2 tr {A c A c } -1 ( Ỹ -α c B) Σc (θ c ) -1 ( Ỹ -α c B) , (5.15)
see (5.22) in the proof of Lemma 5.1. Here, the covariance matrix Σc (θ c ) is defined for all ( j, j ) ∈ {1, . . . , q} 2 by Σc (θ c ) j, j = K c ( t j , t j |θ c ) and B = (b( t1 ), . . . , b( tq )) is a J × q design matrix. In practice, all parameters are replaced using their MLE counterparts and c is selected by minimizing the negative log posterior probability, that is:

c = arg min c p log det Σc ( θc ) + q log det A c A c -2 log(n c /n) + tr A c A c -1 ( Ỹ -αc B) Σc ( θc ) -1 ( Ỹ -αc B) .
In the SITS framework, the above formula provides a natural way to classify a new multivariate time-series even though it is not observed at the same timestamps as the examples from the training set.

Imputation of missing values

The next result provides the distribution of the MGP process at time t † conditionally to its label and to observations at times t 1 , . . . , t q .

Proposition 5.3. Assume that, conditionally to Z = c, Y ∼ MGP p (α c b, K c , A c ) and introduce Y the p×q random matrix defined as Y = (Y(t 1 ), . . . , Y(t q )) where (t 1 , . . . , t q ) ∈ T q . Let t † ∈ T such that t † t k for all k ∈ {1, . . . , q}. Then, conditionally to Z = c and

Y = y , Y(t † ) ∼ N p µ c (t † , y ), Λ c (t † ) , with µ c (t † , y ) = α c b(t † ) + (y -α c B){Σ c (θ c )} -1 k c (t † ), Λ c (t † ) = K c (t † , t † |θ c ) -k c (t † ) Σ c (θ c ) -1 k c (t † ) ⊗ A c A c ,
and where

k c (t † ) = (K c (t † , t 1 |θ c ), . . . , K c (t † , t q |θ c )) . Recall that the covariance matrix Σ c (θ c ) is defined for all ( j, j ) ∈ {1, . . . , q} 2 by Σ c (θ c ) j, j = K c (t j , t j |θ c ) and B = (b(t 1 ), . . . , b(t q )) is a J × q design matrix.
As a consequence, when Y(t † ) is not observed (but its label is known to be c), this missing value can be imputed by the conditional expectation given in Proposition 5.3, where the unknown parameters are replaced by their associated MLE:

Ŷc (t † ) = αc b(t † ) + (Y -αc B){Σ c ( θc )} -1 kc (t † ). (5.16)
This allows for the reconstruction of SITS values at unobserved times. If the label of Y is unknown, the distribution of the MGP process at time t † conditionally to observations at times t 1 , . . . , t q can still be derived from Proposition 5.3:

conditionally to Y = y , Y(t † ) ∼ C c=1 P(Z = c|Y = y )N p µ c (t † , y ), Λ c (t † ) , leading to µ(t † , y ) = C c=1 P(Z = c|Y = y )µ c (t † , y ), Λ(t † ) = C c=1 P(Z = c|Y = y ) Λ c (t † ) + µ c (t † , y ) µ c (t † , y ) -µ(t † , y ) µ(t † , y ).
Thus, when both Y(t † ) and its label are not observed, Y(t † ) can be imputed by

Ŷ(t † ) = C c=1 P(Z = c|Y ) Ŷc (t † ), (5.17)
where Ŷc (t † ) is given in (5.16), Y = (Y(t 1 ), . . . , Y(t q )) and

P(Z = c|Y ) = πc p(Y |Z = c) C k=1 πk p(Y |Z = k) ,
with p(Y |Z = k) the estimated matrix-variate density defined similarly to (5.15) by

p(Y |Z = k) = (2π) -pq/2 det Σ k ( θk ) -p/2 det A k A k -q/2 × exp - 1 2 tr A k A k } -1 (Y -αk B)Σ c ( θk ) -1 (Y -αk B) .
distribution, c ∈ {1, 2}, b ∈ {1, . . . , 10} and j ∈ {1, . . . , 11}. The covariance operator is identical for both classes:

K 1 (•, •) = K 2 (•, •).
It is defined following (5.8) as the sum of a RBF kernel and a white noise covariance function.

The associated parameters are θ 1 = {γ 1 , h 1 , σ 1 } = {1.5, 150, 0.05} = θ 2 . We also set A 1 = A 2 with

A 1 A 1 =                     1 β • • • β β 1 • • • . . . . . . . . . . . . β β • • • β 1                     , (5.18)
so that β tunes the pairwise correlation between the 10 coordinates of the Gaussian processes. In the following, we shall consider β ∈ {0, 1/4, 1/2}. In practice, M2GP processes are simulated on random grids of varying size q ∈ {10, 20, . . . , 100}, see Fig. 5.3 for an illustration in the case q = 10 and β = 0.

Estimation results

All estimation procedures are evaluated on 100 replications of the above described simulation model. First, for all c ∈ {1, 2}, the quality of the reconstructed mean Mc = αc B is measured by the normalized Mean Squared Error (nMSE) defined as:

nMSE( Mc , M c ) = M c -Mc 2 F M c -M c 2 F , (5.19)
where M c is the empirical mean of the processes in class c. The lower this score is, the better the estimation. An example of reconstructed mean is presented on Fig. 5.4, for one replication. Second, the quality of the estimation of the covariance structure A c A c (see 5.18) by A c A c is assessed by the cosine score defined as:

C( A c A c , A c A c ) = 1 - A c A c , A c A c F A c A c F A c A c F . (5.20) Let us note that C( A c A c , A c A c ) ∈ [0, 2] with C( A c A c , A c A c ) = 0 when A c A c and A c A c are proportional.
Finally, turning to the estimation of the kernel part (5.8) of the dependence structure, we focus on the estimation accuracy of the length-scale by computing the absolute difference between the true length-scale h 1 = h 2 = 150 and its estimated counterpart. The results are averaged over the 100 independent replications and are reported on Fig. 5.5 for the first class. Similar results are obtained for the second one. It appears that, unsurprisingly, the quality of the estimates increases with the number q of discretization times. At the opposite, the dependence parameter β does not seem to influence much the accuracy of the estimation. One can nevertheless note that, as expected, the variability of the estimators increases with β, as the information carried by correlated coordinates decreases. Besides, the estimated length-scales do not depend on β, this may be explained by the separability property exhibited in Corollary 5.0.1.

Classification and imputation results

Here, we focus on the comparison between results associated with M2GP and MIGP models. To assess the classification and imputation performances, 4, 000 samples are generated following the model described in Paragraph 5.5.1 and then split into two disjoint balanced sets. The first one is used as a training set (of size n = 2, 000) to estimate model parameters. The second one is used as a test set where the accuracy of the classification and imputation steps associated with the two above methods are compared. The classification performance is assessed thanks to the Overall Accuracy (OA), that is the ratio of the number correctly classified test observations and the total number of test observations, while the nMSE is used for the imputation task. Similarly to (5.19), we let

nMSE( Ŷ , Y ) = Ŷ -Y 2 F Y -Y 2 F , (5.21)
where Ŷ is the imputed discretized process when the class is unknown thanks to (5.17), given the observed discretized process on q points. Y is the empirical mean of discretized processes in the test set. The above Frobenius norms are computed on a fixed regular grid of T defined as {t = /100, = 1, . . . , 100}. The results are reported in Fig. 5.6. It appears that the classification scores associated with M2GP increase with the dependence coefficient β and the number q of discretization times. On the opposite, MIGP scores are decreasing with β, due to the independence assumption. When there is no dependence between coordinates (β = 0), both methods provide similar classification scores. Unsurprisingly, M2GP outperforms MIGP as soon as a dependence occurs.

In terms of reconstruction, both methods feature similar performances, increasing with q. The dependence strength only impacts the variance of the reconstructed processes: The larger β is, the larger the variability.

Time-series classification: Application to satellite data

This section is devoted to multivariate SITS classification using the M2GP model. The data were acquired by the Sentinel-2 satellite, and are presented in Paragraph 5.6.1, with a focus on the irregular temporal sampling. The estimated M2GP parameters are interpreted and discussed in Paragraph 5.6.2. Finally Paragraph 5.6.3 concludes this section with classification results and comparisons to state-of-the-art methods.

Sentinel-2 satellite image time-series

Since 2016, the Sentinel-2 mission [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF] produces massive multispectral images, 2 around 1.6TBytes a day, with a spatial resolution of 10 m/pixel and 13 spectral bands (only 10 bands are used for the analysis). The frequency of revisit is 5 days and clouds as well as shadows are present in the data, at random locations. Most of the clouds and shadows positions are automatically extracted by the data provider. Yet, thin clouds may remain in the data. The selected images cover the area of Toulouse, France (Fig. 5.7) and all available acquisitions for the year 2018 were used. The image is of spatial size 10,000×10,000 pixels (10,000 km 2 ). Each extracted time-series i has a dimension of p = 10 channels (or bands) and its own number of timestamps q i . The distribution of the q i s is represented in Fig. 5.8 for this area in 2018. The supervised classification task consists in assigning a pre-defined label to every pixel of the image. Fourteen classes were extracted from national data-bases and 10 pairs of training and validation data-sets are generated independently for the experiments by randomly selecting samples for the training and testing sets. Training and testing sets were carefully constructed to avoid spatial dependence between pixels.

Table 5.1 shows the number of extracted samples for each training and validation set. The number of samples per class is unbalanced but represents the actual proportion of land cover classes in the region.

Parameters estimation

M2GP is fitted to the satellite image time-series using the estimators described in Section 5.4. A Fourier basis is adopted for estimating the means using J = 19 functions while the time dependence structure is modeled by a RBF kernel combined with an additive white noise. The choice of the basis and the selection of the dimension J are Estimated mean functions are reported in Fig. 5.9 for four selected channels: blue, green, red and near infrared (nIR) and four selected classes: continuous urban fabric, summer crops, broad-leaved forest and water bodies. In the context of remote sensing data, nIR is often correlated with the presence or absence of vegetation: Large values of nIR associated with small values of red, indicate that the vegetation is abundant. This behavior is observed in agricultural classes such as summer crops or broad-leaved forest during spring and summer.

The estimated covariance matrices between all 10 channels A c A c are reported in Fig. 5.10 for the same classes. Similar covariance matrices have already been observed on mono-temporal Sentinel-2 data, we refer to [START_REF] Wang | Assessment of Sentinel-2 MSI Spectral Band Reflectances for Estimating Fractional Vegetation Cover[END_REF]Fig. 8] for similar results on crops classes.

Finally, the time covariance structure is illustrated on Fig. 5.11. The estimated RBF kernel on the same four classes is drawn when centered at day 180. The temporal correlation associated with natural elements, such as summer crops or broad-leaved forest, is short since their reflectance evolves along the year (e.g. because of the vegetation cycle, or anthropic events). In contrast, man-made materials, such as continuous urban fabrics, exhibit longer temporal correlation because their reflectance does not evolve along the time. 3

Classification results

In this section, the classification performances of M2GP are compared to state-of-the-art methods. Four competitors are considered: Random forests (RF) [START_REF] Breiman | Random Forests[END_REF], Quadratic discriminant analysis (QDA) which is based on a finite-dimensional Gaussian model, linear Support vector machine (SVM) classifier fitted with a Stochastic Gradient Descent [START_REF] Zhang | Solving Large Scale Linear Prediction Problems Using Stochastic Gradient Descent Algorithms[END_REF], and, finally, Mixture of independent Gaussian processes (MIGP) [START_REF] Constantin | Joint supervised classification and reconstruction of irregularly sampled satellite image times series[END_REF]).

The time-series have been resampled on a common temporal grid of size 73 (every 5 days of year 2018) using a linear interpolation for RF, QDA and SVM methods since they require a fix vectorial representation of the sample. All the spectral bands have been stacked together to obtain a vector of dimension 73 dates × 10 spectral bands = 730 features. RF is trained with 100 trees of depth 25, and QDA is used with a regularized version of the estimated covariance matrix [START_REF] Friedman | Regularized Discriminant Analysis[END_REF], Σ = (1 -) Σ + I, with = 10 -2 .

The F 1 -score is computed to assess numerically the classification accuracy. The F 1 is defined as the harmonic mean of the precision and recall scores [START_REF] Tharwat | Classification assessment methods[END_REF]. Classification maps are also presented in order to qualitatively evaluate the spatial coherency of the results (despite a spatial pixel-wise independence assumption made by all considered methods).

Means F 1 scores and their standard deviations computed on 10 independent sets are reported in Table 5.2 for each class as well as the "average F 1 score" computed on all classes. Nnon-parametric methods (RF and SVM) provide the best classification results in terms of F1-score. The uni-modal assumption induced by Gaussian models may thus be ill-adapted to this data-set. M2GP and QDA provide lower and similar accuracy, even though M2GP is based on stronger assumptions on the covariance structure than QDA.

The obtained classification maps are reported in Fig. 5.12 for 3 different sites. Large differences are observed in these scenes. For the first column, corresponding to the airport zone, most of the inner vegetations are wrongly classified to natural grasslands with QDA, while RF, SVM and M2GP classify correctly them as meadow. Runway are mostly confused with industrial/commercial units using RF while runways are almost recovered by M2GP. Overall, strong differences between thematic maps are observed, but visual assessment from a mono-date color image is difficult. Yet, without taking into account the spatial dependence, M2GP recovers most of the spatial structure of the image, and the salt and pepper classification noise is limited, as for RF and SVM.

Discussion

A multivariate Gaussian process model has been introduced for the classification of irregularly sampling satellite image time-series. The multivariate model involves a specific structure of the covariance operator that exploits the data features and also reduces the number of parameters to estimate. Furthermore, the proposed formulation scales linearly w.r.t. the number of samples. Experimental results on simulated and real data sets show the importance of modeling the dependence between coordinates of the process, in particular for classification accuracy.

Current development concerns the use of two satellite sources. Sentinel-2 satellites are complemented with Sentinel-1 ones (which are not affected by clouds) which acquire radar data (with a different physical content): An extension of the proposed model will consist in combining these two time-series with irregular temporal and spectral sampling. Another possible extension would be to consider a non spatially stationary mean function, as in [START_REF] Cressie | Statistics for Spatial Data[END_REF].

Finally, these models can be extended to non-Gaussian processes, e.g. Student-t as in [START_REF] Shah | Student-t Processes as Alternatives to Gaussian Processes[END_REF][START_REF] Chen | Multivariate Gaussian and Student-t process regression for multioutput prediction[END_REF] and applied to the unsupervised classification problem. A) and introduce Y the p × q random matrix defined as Y = (Y(t 1 ), . . . , Y(t q )) where (t 1 , . . . , t q ) ∈ T q . From (5.1), we have Y = AW + m with W ∼ IGP p (0, K).

Let W = (W(t 1 ), . . . , W(t q )) be the associated p × q random matrix. Our first goal is to prove that W ∼ MN p,q (0, Σ, I p ) or, equivalently, from (5.3), to prove that vec(W ) ∼ N pq (0, Σ ⊗ I p ). To this end, let us consider the random variable

S = p b=1 q j=1
λ b, j W b, j , and let us prove that S is a Gaussian random variable pour all λ b, j ∈ R pq . Clear, one also has

S = p b=1 S b , with S b := q j=1 λ b, j W b, j = q j=1 λ b, j W b (t j ),
where S 1 , . . . , S p are independent centered Gaussian random variables with variance var(S b ) = q j=1 q j =1 λ b, j λ b, j Σ j, j .

As a consequence, S is a centered Gaussian random variable with variance

var(S ) = p b=1 var(S b ) = p b=1 p b =1 q j=1 q j =1 λ b, j λ b, j Σ j, j × (I p ) b,b .
As a conclusion, vec(W ) ∼ N pq (0, Σ ⊗ I p ) and thus W ∼ MN p,q (0, Σ, I p ). Finally, Y = AW + M ∼ MN p,q (M, Σ, AA ), see [START_REF] Dawid | Some matrix-variate distribution theory: Notational considerations and a Bayesian application[END_REF]Example 1].

Proof of Lemma 5.1 Combining (5.4) and (5.9) yields that the density of Y i, conditionally to Z i = c is given for all i = 1, . . . , n by

p i,c (y) = (2π) -pq i /2 det(Σ c,i (θ c )) -p/2 det(A c A c ) -q i /2 × exp - 1 2 tr (A c A c ) -1 (y -α c B i )(Σ c,i (θ c )) -1 (y -α c B i ) . (5.22)
The likelihood is thus defined as

C c=1 i|Z i =c p i,c (Y i, ),
and the negative log-likelihood can be written as

L = - C c=1 i|Z i =c log p i,c (Y i, ) := 1 2 C c=1 c (α c , θ c , A c A c ) + p log(2π) 2 C c=1 i|Z i =c q i ,
with, for all c = 1, . . . , C,

c (α c , θ c , A c A c ) = p i|Z i =c log det(Σ c,i (θ c )) + i|Z i =c q i log det(A c A c ) + i|Z i =c tr (A c A c ) -1 (Y i, -α c B i )(Σ c,i (θ c )) -1 (Y i, -α c B i ) .
The conclusion follows.

Proof of Proposition 5.

2 (i) Let us first consider the differential of c (α c , θ c , A c A c ) w.r.t. α c : d c (α c ) := i|Z i =c tr {A c A c } -1 d(Y i, -α c B i ){Σ c,i (θ c )} -1 (Y i, -α c B i ) = - i|Z i =c tr {A c A c } -1 (dα c )B i {Σ c,i (θ c )} -1 (Y i, -α c B i ) - i|Z i =c tr {A c A c } -1 (Y i, -α c B i ){Σ c,i (θ c )} -1 (B i ) (dα c ) = -2 i|Z i =c tr {A c A c } -1 (Y i, -α c B i ){Σ c,i (θ c )} -1 (B i ) (dα c ) ,
by remarking that both terms are equal in view of the properties of the trace operator. Moreover, from Kronecker product properties [START_REF] Schott | Matrix Analysis for Statistics[END_REF]Theorem 8.12], one has

d c (α c ) = -2 i|Z i =c vec(dα c ) B i {Σ c,i (θ c )} -1 ⊗ {A c A c } -1 vec(Y i, -α c B i ) = -2(dvec(α c )) vec         {A c A c } -1 i|Z i =c (Y i, -α c B i ){Σ c,i (θ c )} -1 (B i )         .
Interpreting the above result as a scalar product and using the "broad" definition of matrix derivative defined in [START_REF] Magnus | On the concept of matrix derivative[END_REF], if follows:

∂ c (α c , θ c , A c A c ) ∂α c = -2vec         {A c A c } -1 i|Z i =c (Y i, -α c B i ){Σ c,i (θ c )} -1 (B i )         .
Setting this partial derivative to zero yields

i|Z i =c (Y i, -α c B i ){Σ c,i (θ c )} -1 (B i ) = 0,
or equivalently,

α c =         i|Z i =c Y i, {Σ c,i (θ c )} -1 (B i )                 i|Z i =c B i {Σ c,i (θ c )} -1 (B i )         -1
, which is the desired result. Second, let us consider the differential of

c (α c , θ c , A c A c ) w.r.t. A c A c : d c (A c A c ) = Q c d log det(A c A c ) + dtr N(θ c ){A c A c } -1 , 5.8. Appendix -Proofs 81 where N(θ c ) = i|Z i =c (Y i, -α c B i ){Σ c,i (θ c )} -1 (Y i, -α c B i ) .
From [START_REF] Schott | Matrix Analysis for Statistics[END_REF]Example 9.6], the associated partial derivative vanishes for

A c A c = N(θ c ) Q c = 1 Q c i|Z i =c (Y i, -α c B i ){Σ c,i (θ c )} -1 (Y i, -α c B i ) ,
and the result is proved.

(ii) Introduce β c,i (α c , θ c ) = (Y i, -α c B i ){Σ c,i (θ c )} -1 and consider the kth coordinate of the gradient of c (α c , θ c , A c A c ) w.r.t. θ: ∂ c (α c , θ c , A c A c ) ∂θ k = p i|Z i =c ∂ ∂θ k log det(Σ c,i (θ c )) + ∂ ∂θ k tr N(θ c ){A c A c } -1 = p i|Z i =c tr {Σ c,i (θ c )} -1 ∂Σ c,i (θ c ) ∂θ k - i|Z i =c tr (Y i, -α c B i ){Σ c,i (θ c )} -1 ∂Σ c,i (θ c ) ∂θ k {Σ c,i (θ c )} -1 (Y i, -α c B i ) {A c A c } -1 = p i|Z i =c tr {Σ c,i (θ c )} -1 ∂Σ c,i (θ c ) ∂θ k - i|Z i =c tr β c,i (α c , θ c ) {A c A c } -1 β c,i (α c , θ c ) ∂Σ c,i (θ) ∂θ k = i|Z i =c tr p{Σ c,i (θ c )} -1 -β c,i (α c , θ c ) {A c A c } -1 β c,i (α c , θ c ) ∂Σ c,i ∂θ k (θ c ) .
The result is proved. c ) is defined for all ( j, j ) ∈ {1, . . . , q} 2 by Σ c (θ c ) j, j = K c (t j , t j |θ c ) and B = (b(t 1 ), . . . , b(t q )) is a J × q design matrix. Let t † ∈ T be an unobserved time, i.e. t † t k , for all k ∈ {1, . . . , q}, and k c (t † ) = (K c (t † , t 1 |θ c ), . . . , K c (t † , t q |θ c )) . Then, classical properties on conditional Gaussian random vectors (see for instance [17, p. 63]) entail that, conditionally to Z = c and vec(Y ) = vec(y ), Y(t † ) follows the p-variate Gaussian distribution N p (µ c (t † , y ), Λ c (t † )) with, on the one hand

µ c (t † , y ) = α c b(t † ) + [k c (t † ) ⊗ A c A c ]{Σ c (θ c ) ⊗ A c A c } -1 vec(y -α c B) = α c b(t † ) + [k c (t † ) ⊗ A c A c ]{Σ c (θ c ) -1 ⊗ (A c A c ) -1 }vec(y -α c B) = α c b(t † ) + {k c (t † ) Σ c (θ c ) -1 } ⊗ {(A c A c )(A c A c ) -1 } vec(y -α c B) = α c b(t † ) + {k c (t † ) Σ c (θ c ) -1 } ⊗ I p vec(y -α c B) = α c b(t † ) + vec I p (y -αc B){k c (t † ) Σ c (θ c ) -1 } = α c b(t † ) + (y -α c B)Σ c (θ c ) -1 k c (t † ),
and on the other hand,

Λ c (t † ) = K c (t † , t † |θ c ) ⊗ A c A c -[k c (t † ) ⊗ A c A c ] Σ c ⊗ A c A c (θ c ) -1 [k c (t † ) ⊗ A c A c ] = K c (t † , t † |θ c ) ⊗ A c A c -(k c (t † ) Σ c (θ c ) -1 ) ⊗ I p [k c (t † ) ⊗ A c A c ] = K c (t † , t † |θ c ) ⊗ A c A c -k c (t † ) Σ c (θ c ) -1 k c (t † ) ⊗ (I p A c A c ) = K c (t † , t † |θ c ) -k c (t † ) Σ c (θ c ) -1 k c (t † ) ⊗ A c A c .
The result is proved.

Conclusion and perspectives

Summary

This thesis defines two statistical models for remote sensing data-sets, i.e. multi-spectral satellite images timeseries (SITS) from Sentinel-2 constellation, to classify irregularly sampled time-series. To this end we defined a Gaussian process (GP) over functions. Any irregularly sampled time-series is a marginal from the process and is distributed according to a multi-variate Gaussian distribution. We presented two models to overcome the classification problem or irregularly sampled SITS. The first model introduced in Chapter 4 assumes independence between spectral bands. It learns as many one-dimensional GP as the number of spectral bands. This model shows interesting behaviour, particularly in time-series reconstructions where the reconstructed reflectance value is as good as with non-parametric techniques, or sometimes better. The reconstruction also provides an uncertainty measure which is not always provided by other techniques. It has been applied to Sentinel-2 SITS and the provided information may be useful for the study of some phenomena (agricultural mowing period, etc). The classification scores are small when compared to state-of-the-art methods. The choice of various kernels and the combination of them have not shown significant improvements, the presented scores are the best we could obtain. Among others, we identified the independence assumption as too strong for the representation of SITS. The second model has been introduced in Chapter 5 to tackle the issue of independence between spectral bands. The model takes into account a linear dependency between the bands. The reconstruction quality has not been changed when compared to the previous model. However, the classification scores have been significantly improved. This model is slightly bellow state-of-the-art performances but the land use or land cover maps from irregularly sampled SITS is more reliable than the first approach.

To conclude, the issue of irregularly sampled time-series has been tackled. We were able to demonstrate the possibility to classify time-series without temporal re-sampling on large real-world remote sensing data-sets. The second model provides sufficient classification accuracy to rely on the produced maps. There are still some issues, in my opinion, some of them are linked to the Gaussian assumptions. Indeed, mislabelled data combined with variations within the same class linked to external factors (for example altitudes which modify the reflectance) yields an increasing noise learned by GP. This has been observed on areas with mountains where altitude changes the seasonality.

Future works

We identified three axes to extend the model and improve the performances of classification of irregularly sampled time-series. The first axis discusses non-Gaussian models, particularly Student-t processes. The second one discusses the three dimensional aspect of SITS by incorporating spatial variability. Finally the last axis concludes on more general covariance models with application to multi-sensor data fusion.

Non-Gaussian modelling

The Gaussian assumption done on the processes is, for SITS at least, limiting. The noise within SITS and the mislabelled data induce undesired outliers and it is known that the Gaussian distribution is very sensitive to outliers.

One usual solution within the past ten years is to assume a Student-t distribution. The Student-t distribution, or t-distribution, is more robust to outliers thanks to a heavier tail (see Figure 2.1). In [START_REF] Shah | Student-t Processes as Alternatives to Gaussian Processes[END_REF] they proposed a definition of the Student-t process and showed interesting properties to model noisy functions. It may provide interesting reconstruction of SITS.

Additionally, the multi-variate Student-t processes has also been studied in [START_REF] Chen | Multivariate Gaussian and Student-t process regression for multioutput prediction[END_REF] and is interesting to take into account spectral dependency of SITS.

Incorporating spatial variability

An other interesting work is to take into account the spatial information from SITS. Spatio-temporal [START_REF] Cressie | Statistics For Spatio-Temporal Data[END_REF] studies of SITS are of interest in Remote Sensing and already applied to Sentinel-2 data-set [START_REF] Gbodjo | Toward Spatio-Spectral Analysis of Sentinel-2 Time Series Data for Land Cover Mapping[END_REF].

A perspective is to learn the spatial information within the mean of the Gaussian process as a function of time and space. It can by extended to a three-dimensional space by taking into account the altitude of the pixel location. It should scale to large data-sets as the pixels are seen independently and because the spatial information only changes the mean.

An other perspective is to use this information within the covariance operator (following the idea of three-order tensor distribution) but results in very complex estimation (M2GP is defined up to a multiplicative constant, a three-order tensor is defined up to two multiplicative constants). The use of large spatial areas (up to a complete Sentinel-2 tile) would result to a numerically non-invertible covariance matrix but restrictions to small areas may increase the classification accuracy.

More general covariance models

As introduced in Section 2.2.4, the LMC model (2.19) has been introduced in the context of Gaussian distributions over multi-outputs processes. LMC and its derived models issued from the geostatistics literature [START_REF] Goovaerts | Geostatistics for Natural Resources Evaluation[END_REF] and are used in a regression context, for example on time-series [START_REF] Cheng | Sparse multi-output Gaussian processes for online medical time series prediction[END_REF] (in this work, it is interesting to highlight that authors considered irregular time stamps but removed this specificity by considering one patient). LMC is a covariance model which assumes independent Gaussian distributions over latent processes resulting to a Gaussian posterior distribution [2,[START_REF] Teh | Semiparametric latent factor models[END_REF] over a linear combination over these latent processes. The number of latent processes may be different from the number of outputs. However, if all the latent processes share the same covariance operator, it yields the ICM model (2.20) which structure, based on the Kronecker product, generalizes M2GP. Inference on LMC, or ICM, is then also done using the Maximum Likelihood Estimator [2, Section 6.2] and is cubic w.r.t. the number of samples whereas M2GP scales linearly. LMC and ICM could be adapted to our context to scale to large data-sets.

An other idea is to consider a correlation between spectral bands which evolves with time input. The distribution is known as convolutional Gaussian processes [START_REF] Àlvarez | Kernels for vector-valued functions: A review[END_REF][START_REF] Cardona | Convolved Multi-output Gaussian Processes for Semi-Supervised Learning[END_REF]. It may be seen as a two-dimensional Gaussian process with different applications. As an example, the M2GP model assumes that, despite an irregular sampling, the spectral bands are observed at the same time stamps. However, within the remote sensing literature, one key application is to use jointly multi-sensors to classify SITS, as an example one may consider the Sentinel-1 (SAR images) and Sentinel-2 (optical images) data fusion [START_REF] Pastorino | Multisensor and Multiresolution Remote Sensing Image Classification through a Causal Hierarchical Markov Framework and Decision Tree Ensembles[END_REF]. Indeed, Sentinel-1 provides complementary information as the information is not modified by the presence of clouds or shadows. The joint use implies an irregular sampling between spectral bands at the same time stamps. In [START_REF] Boyle | Dependent Gaussian Processes[END_REF] they studied the Dependent GP where the idea is to build latent features to retrieve the same assumptions as M2GP model or LMC.

In spite of that, estimation of convolutional GP is challenging. In [2, Section 3] the spectral covariance has been simplified to a Dirac delta function, resulting to the LMC. In [START_REF] Cardona | Convolved Multi-output Gaussian Processes for Semi-Supervised Learning[END_REF] the authors considered a Gaussian kernel between spectral features and between temporal features and simplified the resulting kernel. In both cases, the numerical complexity is as complex as in LMC, adaptation to our context may require additional assumptions particularly for large data-sets. the upper and lower bounds of the Rayleigh quotient, we have:

λ B MB max = sup x x B MBx x x = sup x x B MBx x x x B Bx x B Bx , because B is injective, = sup x x B MBx x B Bx x B Bx x x = sup x (Bx) MBx (Bx) Bx x B Bx x x ≤ sup y y My y y sup x x B Bx x x ≤ λ M max λ B B max . Finally, λ B MB max ≤ λ M max λ B B max .
With the lower bound of the Rayleigh quotient, it is straightforward that:

λ B MB min ≥ λ M min λ B B
min , and yields to the expected upper bound.

In (4.8), a sufficient condition to compute α c numerically is to consider both a kernel and the design matrix such that B c B c that are simple to invert.

A.2.2 Illustration on a toy data-set

This toy example aims to present the estimation of the mean for an ill-conditioned design matrix. To this end we define three means, one for each class (p = 1 spectral band) for the Gaussian processes with different behaviors: one constant line, one generated using Fourier functions and one generated using polynomial functions. The three means are illustrated in Figure A.1. Besides, we define a RBF kernel (with a length scale of 63 days and a multiplicative constant of 1) with an additive white noise (set to 0.5). The noise may is not representative of the Sentinel-2 SITS, neither the range of the mean values, but our aim is to observe the estimation of the mean.

Afterwards, a total of n c = 350 GP samples are generated from the distribution for each class, and are marginalized such as a temporal window remains unobserved. The size of each marginal is randomly selected according to a discrete uniform distribution and ranges between 3 and 13. The simulation is replicated 100 times by generating new samples/ Finally, five simulations settings are computed with an increasing size of the unobserved temporal window, starting from day 100 and ending at 100 + N day with N such as: N 1 = 10 days, N 2 = 20 days, N 3 = 50 days, N 4 = 100 days and N 5 = 170 days.

II. NUMERICAL IMPLEMENTATION AND CONVERGENCE

The model has been implemented with Python and inherits from Scikit Learn the Gaussian Process and Kernel classes [1]. The optimization problem is solved using L-BGFS-B subroutine [2] that allows for additional constraints or bounds (i.e. positivity or minimal/maximal values) on the covariance function parameters.

Minimization of [START_REF] Azzalini | The Multivariate Skew-Normal Distribution[END_REF] in the main document is known to converge to a local minimum, which depends on initial values of the hyperparameters. Figure 2 illustrates the evolution of the parameter h of the RBF kernel as well as the marginal log-likelihood with respect to the number of iterations. Here, a Fourier basis of dimension 19 was used to estimate the mean function. Figure 1 represents the estimated mean function in the infra-red for class "broad-leaved forest" for different dimensions of the Fourier basis. For J > 25, some oscillations appear in the reconstruction, due to a possible overfit. J = 19 seems to be the most appropriate dimension, even though, in terms of classification, J = 25 yields a slightly better accuracy score. In practice, a compromise should be found between a good reconstruction and a better classification score (see Figure 8 in the main document), depending on the final objective.

Figure 3 presents the RBF kernel computed as exp -(t-180) 2 2h 2

, t ∈ {1, ...365} on the infrared wavelength from three different classes. It appears that the continuous urban fabric has a higher temporal correlation than Broad-Leaved forests and Water bodies. Indeed, natural elements reflectance such as vegetation evolve along the year (e.g. phenology) while man-made material's reflectance does not evolve along the year and thus exhibits longer temporal correlation.

IV. CONSEQUENCES OF THE INDEPENDENCE ASSUMPTION

To illustrate the impact of the independence assumption on the classification accuracy, two versions of Quadratic Discriminant Analysis (QDA) have been implemented and compared. The first version, referred as u-QDA is based on the assumption that all the wavelengths are uncorrelated. More specifically, a QDA model is fitted independently to each wavelength. The likelihood is computed as the product of the likelihoods associated with each wavelength, and the posterior probability is obtained using Bayes' rule as usual. The second version is the classical QDA method presented in the main document. The results are reported in Table I. It appears that the classification accuracy of u-QDA is significantly smaller than the one of the usual QDA method. Therefore, one can conclude that the independence assumption has a significant impact on the performance of the classifiers. Fig. 4: Reference data extraction where each color represents an extract of land cover classes. The color code is described in Figure 7. The axes correspond to the geographical coordinates.The red continuous squares correspond to zoomed areas of the classification thematic, and reconstruction, maps provided in Figure 7 and Figure 15 respectively. The background image is an extract from Sentinel-2 optical images.

V. SUPERVISED CLASSIFICATION

This section provides an extract of the thematic classification maps obtained with RF, MIMGP and SVM on tiles T31TCJ (Figure 6), T31DN (Figure 7) and T31TGK (Figure 8). The original tiles are too large to be correctly displayed, we thus choose to extract the center of each tile. All extracts have 1000 × 1000 pixels. We also provide a zoom (200 × 200 pixels) on several extracts.

The first extract from T31TCJ is presented in Figure 7 of the main document. The second and third extracts from, respectively, T31TDN and T31TGK, are presented in Figures 4 and5.

The F 1 scores for each class and for each tiles, averaged over the 10 independent runs, are provided in Tables II, III and IV.

VI. TIME-SERIES RECONSTRUCTION

This Section displays additional temporal reconstructions. Figures 9 and 10 present the reconstructed time-series for all days in year 2018 for two pixels on four wavelengths of different classes, summer crops and broad-leaved forests.

Figure 11 compares reconstructions obtained with MIMGP model and Whittaker smoother on the 10 wavelengths. The more intense the colors are (yellow) on the line of equation y = x, the better the reconstruction is. Both methods have higher density Fig. 5: Reference data extraction where each color represents an extract of land cover classes. The color code is described in Figure 8. The axes correspond to the geographical coordinates. The red continuous squares correspond to zoomed areas of the classification thematic maps provided in Figure 8 and the dashed square corresponds to the reconstruction map given in Figure 16. The background image is an extract from Sentinel-2 optical images. on the y = x line, confirming that there is no bias in the reconstructions, but Whittaker exhibits higher dispersion. This is especially true for tiles T31TDN and T31TGK. Finally, MIMGP and Whittaker image reconstructions are compared visually on the three tiles. The four wavelengths natively at a resolution of 10 meters per pixel (Blue, Green, Red and Infrared) are reported in Figures 14, 15 and 16 for, respectively, tile T31TCJ, T31TDN and T31TGK. For the visualization sake, each image of a given tile has been displayed using the same histogram stretching function and identical parameters. The mask is a binary image, with white color for valid pixels and black color for invalid pixels (clouds, saturation, ....).

Figure 1 . 1 :

 11 Figure 1.1: (Image from[START_REF] Keshava | Spectral unmixing[END_REF]). Surface reflectance measured thanks to an embedded multi-spectral instrument.

Figure 1 . 2 :

 12 Figure 1.2: Electromagnetic spectrum (wavelengths in m). Optical satellites measure reflectance in the visible spectrum and in the infrared. The concerned wavelengths values range from 400 to 2400nm.

1. 1 .Figure 1 . 4 :

 114 Figure 1.4: Sentinel-2 spectral band filters (coloured areas) centered around the wavelength provided in Table1.1 and with associated bandwidths[START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]. Visible light starts from 380nm and ends around 780nm: blue, green and red bands are highlighted in the 10m ground resolution bands. The spectral reflectance signature of healthy grass is reported on each grid (full line) and the signature of dried grass (dashed line).

Figure 1 . 5 :

 15 Figure 1.5: Satellite Image Time-Series true colors. The data come from year 2018 over the center of French territory. Images were downloaded from Theia Land Data Center (https://www.theia-land.fr/en/products/).

Figure 1 . 6 :

 16 Figure 1.6: (Image from [80]). Sentinel-2 product levels: level 1C (with atmospheric interference and clouds) and level 2A (reflectance and clouds removal). Level 3A shows reconstructed satellite image.

Figure 1 . 7 :

 17 Figure 1.7: True color noisy image (left side) with clouds, shadows and saturations. The left image presents the data detected as noisy by[START_REF] Hagolle | A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENµS, LANDSAT and SENTINEL-2 images[END_REF]. Most of the noise has been captured, however some clouds on the top have not been detected. Images were downloaded from Theia Land Data Center (https://www.theia-land.fr/en/ products/).

Figure 1 . 8 :

 18 Figure 1.8: (Image from the supplementary materials of [41]). Polygons classes over a true colors tile in the center of the metropolitan France. Red squares are zoomed areas presented in [41].

Figure 1 . 9 :

 19 Figure1.9: Schematic view of a decision tree which returns as output a class membership for a given input y ∈ E. The class membership associated with y is given by the red path (c = C -1) in this example.

14 Chapter 1 -Figure 1 . 12 :

 141112 Figure 1.12: Sentinel-2 pixel time series for the near infra-red band (B8) for three different tiles over the French territory. The clear dates are marked as blue whereas red marks indicate a detection of noise from the mask. Dashed lines represent the gap-filled data.

18 Chapter 2 -

 182 Statistical modelling for time-series classification

Chapter 2 -Figure 2 . 1 :

 221 Figure 2.1: Probability density functions (p) for 1-dimensional Gaussian distribution (solid line), skew-normal distribution (dashed line) with γ = 2 and the t-distribution (dotted line) with ζ = 1.

Figure 2 . 2 :

 22 Figure 2.2: LDA and QDA decision boundaries between two classes (c = 1 in red and c = 2 in blue).

Definition 2 . 1 (

 21 vec operator). Let C a m × n matrix. vec(C) is a vector of size mn such that: vec(C) = (c 11 , . . . , c m1 , c 12 , . . . , c m2 , . . . , . . . , c 1n , . . . , c mn ).

Definition 2 . 2 (

 22 Kronecker product). Let C a m 1 × n 1 matrix and D an other matrix of size m 2 × n 2 . The left Kronecker product, or Kronecker product, C ⊗ D is a matrix of size m 1 m 2 × n 1 n 2 such that:

Figure 2 . 3 :

 23 Figure 2.3: Separability between lines and columns for a centered random matrix Y. The k-th line of Y does not depend on covariance matrix Λ and the l-th line does not depend on Σ.

Chapter 2 -

 2 Statistical modelling for time-series classification

Figure 2 . 5 :

 25 Figure 2.5: Two samples from the same centered GP (∀t ∈ T = [0, 1], m(t) = 0) with a RBF kernel, (2.11) with length-scale = 0.25. The dots represent the marginals from each sample.

Chapter 3

 3 Model-based classification for irregularly sampled time-series Outline French introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1 Classification of irregularly sampled time-series . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.1.1 Continuous representation of SITS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.1.2 One-dimensional process and the Mixture of Independent multivariate Gaussian processes 34 3.1.3 Mixture of multivariate Gaussian processes . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Reconstruction of noisy time stamps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.1 Imputation with independence assumption . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.2 Imputation with M2GP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2.3 Imputation using the complete mixture . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3 Implementation -python code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

  ) where (b, b ) = {1, . . . , p} 2 , W b ∼ GP(0, K b,c ) are independent latent processes. For b b , the independence property is denoted by W b W b . m c : T → R p is a vector mean function and is presented in Section 5.3 and K b,c is the b-th covariance operator for class c. Model (3.2) results in the function-space view [190, Section 2.2] framework as it defines a distribution over functions.

Figure 3 . 1 :

 31 Figure 3.1: Simulated satellite image time series at bands 1 and p. The dots below the temporal axis correspond to the temporal sampling.

2 .

 2 The associated classification and reconstruction methods are derived in Section 4.3. Section 4.4 details the SITS datasets while the experimental set-up is presented in Section 4.5. Experimental classification and reconstruction results are provided in Sections 4.6 and 4.7 respectively. Finally, Section 4.8 concludes this paper and opens discussion on future work.

Figure 4 . 1 :

 41 Figure 4.1: Acquisition dates for tile T31TDN (red), T31TCJ (black) and T31TGK (blue) in 2018. The location of the tiles is given in Fig. 4.5. Large temporal gaps between two dates correspond to very cloudy acquisitions: Images tagged as composed of more than 90% of cloudy pixels are not processed by the French data provider.

Figure 4 . 2 :

 42 Figure 4.2: Simulated univariate Gaussian processes with squared exponential covariance function (length scale parameter set to 50). The continuous black line is the mean function and the dashed lines are the corresponding 10 realizations.

Definition 4 . 1 (

 41 Mixture of Independent Multivariate Gaussian Processes). The proposed model, namely Mixture of Independent Multivariate Gaussian Process (MIMGP), relies on two main assumptions: A1 Each process Y b , conditionally to Z = c, follows a Gaussian Process: Y b |Z = c ∼ GP(m b,c , K b,c ); A2 Conditionally to Z = c, all components Y b of Y are independent,

1 for c=1 to C do 2 for b=1 to p do 3 repeat 4 5

 2345 Update α b,c using (4.8); Do a gradient step w.r.t. θ b,c using (4.9); 6 until b,c (α b,c , θ b,c ) has converged; Such expression is a consequence of the use of non-zero mean GPs [190, Section 2.7]. The associated negative marginal log-likelihood is given by (α, θ) = p,C b,c=1 b,c (α b,c , θ b,c ), with b,c (α b,c , θ b,c ) =

Figure 4 .

 4 2 (i.e., C = 1 and p = 1) leads to the estimation of the mean function displayed on Figure 4.3.

Figure 4 . 3 :

 43 Figure 4.3: The continuous line is the mean function associated with the GP displayed in Figure 4.2 and the dashed line is the estimated one using Algorithm 2 and J = 20 Gaussian functions (see Section 4.5.1 for details on bases).

Figure 4 .

 4 6 displays the distribution of the number of clear dates (dates not tagged clouds, shadows or no-data in the raw time-series): It appears that the number of clear acquisitions per pixel significantly varies depending on the tile.

Figure 4 . 4 :

 44 Figure 4.4: Time-series reconstruction. The black continuous line is one (continuous) realization of the GP from Figure 4.2. The black dots are the observed noisy acquisitions at several (discrete) times. The dashed line is the reconstructed time-series. The gray region displays a pointwise confidence region for the reconstruction: Ŷ j b (t * ) ± √ V[ Ŷ j b (t * )].

Figure 4 . 5 :

 45 Figure 4.5: Location of the Sentinel-2 tiles. The label inside the square is the tile name.

Figure 4 . 6 :

 46 Figure 4.6: Proportion of pixels as functions of the number of clear dates.

Figure 4 . 7 :

 47 Figure 4.7: Reference data extraction where each color represents an extract of land cover classes from Table4.2. The color code is described in Figure4.9. The axes correspond to the geographical coordinates. The white continuous squares correspond to zoomed areas of the classification thematic maps provided in Section 4.6. The background image is an extract from Sentinel-2 optical images.

4. 5 .

 5 Experimental set-up 51

Figure 4 . 8 :

 48 Figure 4.8: Averaged "mean F 1 scores" of the MIMGP model with the 4 bases as a function of the basis dimension. The results are reported for the tile T31TCJ.

Figure 4 . 9 :

 49 Figure 4.9: Thematic maps for 3 sites from tile T31TCJ. These sites are located as illustrated on Figure 4.7.

54 Chapter 4 -

 544 Joint supervised classification and reconstruction of irregularly sampled satellite image times series

Figure 4 . 10 :

 410 Figure 4.10: Image reconstruction from a site in tile T31TGK. This site is located as illustrated on Figure 5 of the supplementary materials. The top left panel is the original spectral band 4 (red) on August 6, 2018. The top right panel corresponds to the mask: white color for a clear pixel and black color for pixel with clouds/saturation as detected by the pre-processing chain. Bottom left and right are the reconstruction with the Whittaker smoother and the MIMGP, respectively.

Figure 4 . 11 :

 411 Figure 4.11: Time-series reconstruction where top and bottom figures represent the red reflectance for two different pixels. Dot points are the original values from the irregularly sampled time-series. Blue points are non-detected clouds and saturation by the satellite data preprocessing chain, for the top and bottom figure, respectively. The black continuous line is the conditional reconstruction of the signal and the gray region corresponds to the 95% confidence interval. Each day of the year was reconstructed with MIMGP. The red dashed line represents the linear interpolation taking into account all the dates tagged "clear", the black dotted line represents the smoothed signal using Whittaker smoother.

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: True color Sentinel-2 satellite image time-series. Data were acquired in 2018 at different time steps over the area of Toulouse, France (images were downloaded from Theia Land Data Center: http://www.theialand.fr/en/presentation/products).

  is the vector obtained by stacking the n columns of C: vec(C) = (c 11 , . . . , c m1 , c 12 , . . . , c m2 , . . . , c 1n , . . . , c mn ) .

Corollary 5 . 0 . 1 .

 501 Suppose the assumptions of Proposition 5.1 hold. (i) For all (b, b ) ∈ {1, . . . , p} 2 and j ∈ {1, . . . , q}, one has ρ(Y b, j , Y b , j ) = cos (a b , a b ),

5 •Figure 5 . 3 :

 553 Figure 5.3: Two simulated M2GP processes (transparent lines) in dimension p = 10 observed at q = 10 timestamps (dots), from two classes (c = 1: blue, c = 2: red). The mean functions are depicted as continuous opaque lines. Top panel: first coordinates, bottom panel: second coordinates (only the first two coordinates p 1 and p 2 are represented).

5. 6 . 71 - 1 - 1 -Figure 5 . 4 :

 6711154 Figure 5.4: Estimation of mean functions by M2GP on simulated data for all coordinates b ∈ {1, . . . , 10}, classes c ∈ {1, 2} and β = 0 on one replication. The dashed line is the true mean, the red line is the estimated GP mean from a discretization on a grid of size q = 10.

3 CFigure 5 . 5 :

 355 Figure 5.5: Estimation of M2GP parameters on simulated data as a function of the number q of discretization times on class c = 1. From left to right: normalized mean squared error (5.19), cosine score (5.20) and absolute difference of length-scales. From top to bottom: β = 0, β = 1/4 and β = 1/2.

5. 6 .Figure 5 . 6 :

 656 Figure 5.6: Classification overall accuracy (OA, left panel) and reconstruction normalized mean-squared error (nMSE, right panel in log scale) boxplots computed on simulated data. Comparison between M2GP (blue) and MIGP (red) results as functions of the number q of discretization times. From top to bottom: β = 0, β = 1/4 and β = 1/2.

Figure 5 . 7 :

 57 Figure 5.7: The study area is located in the south of France (right bottom image). The left bottom image corresponds to the entire area (100 km×100 km) and the upper image is a zoom over the red rectangle (11 km×5 km).

Figure 5 . 8 :

 58 Figure 5.8: Normalized histogram of the q i s within the SITS data-set.

5. 6 .

 6 Time-series classification: Application to satellite data 75

76 Chapter 5 -Figure 5 . 9 :

 76559 Figure 5.9: Estimated means for four channels and four classes (continuous urban fabric, summer crops, broadleaved forest and water bodies). The horizontal axis represents the days of the year and the vertical axis represents the reflectance value.

Figure 5 . 10 :Figure 5 . 11 :

 510511 Figure 5.10: Estimated covariance matrices A c A c on four land-cover classes.

Proof of Proposition 5 . 3

 53 Let Y = (Y(t 1 ), . . . , Y(t q )) be a p×q random matrix where, conditionally to Z = c, Y ∼ MGP c (α c b, K c , A c ). Recall that Proposition 5.1 yields vec(Y ) ∼ N pq (vec(α c B), Σ c (θ c ) ⊗ A c A c ), where Σ c (θ

3 Figure A. 1 :

 31 Figure A.1: Means of the 3 classes.

Figure A. 2 :

 2 Figure A.2: Estimated mean function for the 1st class when the average number of time stamps and the number of basis functions J are identical but the size of the unobserved window is increasing (N 4 > N 3 > N 2 > N 1 ), delimited by the dotted vertical lines.

1

 1

IIIFig. 1 :

 1 Fig.1: Reconstructed mean near-infrared (IR) time-series from broad-leaved forest class with parameters α (continuous line) estimated on the raw data, for different dimensions of basis. The dashed line corresponds to the gap-filled data mean from QDA on the resampled grid.

Fig. 2 :Fig. 3 :

 23 Fig.2: Evolution of the length-scale h ∈ θ, given in (4) in the main document, and associated log-likelihood for the IR wavelength and three different areas: artificial (continuous urban fabric), semi-natural (broad-leaved forests) and water (water bodies) areas.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Thematic maps for 3 sites from tile T31TCJ. The first line is the whole 10km long side square while the other lines are enlargements (2km length) of some extracts. The associated coordinates are presented in Figure 7 of the main document.

Fig. 8 :

 8 Fig. 8: Thematic maps for 3 sites from tile T31TGK. The first line is the whole 10km long side square while the other lines are enlargements (2km length) of some extracts. The associated coordinates are presented in Figure 5.

Fig. 9 :

 9 Fig. 9: Blue, Green, Red and IR wavelengths reconstruction for the Summer crops class. The red dashed line represents the linear interpolation, the black dotted line the Whittaker reconstruction and the black continuous line the MIMGP reconstruction.

Fig. 10 :

 10 Fig. 10: Blue, Green, Red and IR wavelengths reconstruction for the broad-leaved forests class. The red dashed line represents the linear interpolation, the black dotted line the Whittaker reconstruction and the black continuous line the MIMGP reconstruction.

Fig. 14 :Fig. 15 :Fig. 16 :Figure C. 1 :

 1415161 Fig. 14: Comparison of reconstructions on a 2km long side square from T31TCJ tile on July 12, 2018. The associated site is presented in Figure 7 of the main document.

Figure C. 2 : 4 : 1 #C. 3 . Cython code examples 113 Algorithm 5 : 112 X [ . . . ] = B 114 #

 24131135112114 Figure C.2: Computation time of the computation of posterior probabilites on all classes averaged on 10 runs.

Table 1 .

 1 1.

	Band name	Ground resolution (m) Wavelength (nm)
	B2 (Blue) B3 (Green) B4 (Red) B5 B6 B7 B8 (near infra-red) B8A B11 B12	10 10 10 20 20 20 10 20 20 20	490 560 665 705 740 783 842 865 1610 2190
	B1 B9 B10	60 60 60	443 945 1380
	Table

  . Supervised model-based classification is, more generally, already used in time-series classification. The following Chapter reviews statistical modelling for classification of multidimensional time-series. French introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1 Supervised model-based classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.1 Decision theory and supervised framework . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.2 General discriminative problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1.3 Matrix-variate Gaussian distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Gaussian processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.1 Kernel functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.2 Gaussian processes for regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.3 Gaussian processes for classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2.4 Multi-output Gaussian processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Statistical modelling for Sentinel-2 satellite image time-series . . . . . . . . . . . . . . . . . . 31

	Chapter 2
	Statistical modelling for time-series
	classification
	Outline

  π c is a parameter of the model and is encompassed in θ c , for example in unsupervised classification they have to be estimated through the EM algorithm (Expectation Maximization) [22, Eq. (2.10), p.25]. In our supervised framework, π c are estimated by their empirical counterpart: πc = n c /n and are removed from parameters θ c . MLE has interesting properties like convergence (when n → ∞) to a normal distribution centered on the true value of θ [84, Chapter 8, p.266].

  .1. Classification of irregularly sampled time-series

																						35
											Band 1
	Reflectance	0.5 1																	y 1,1, y 2,1, Y 1
		0
		0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20
												. . .
											Band p
	Reflectance	0.5 1																	y 1,p, y 2,p, Y p
		0
		0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20
												Times

Table 4 .

 4

	f d µ,Σ y|x p(y|x) S n p b T Y y i,b T i C c z i π c α, θ J |K|	Density associated with N d (µ, Σ). Random variable y conditionally on x, Probability distribution of y|x. A set of SITS or pixels. The total number of SITS available in S. Number of wavelenghts of a SITS. Band number, in {1, . . . , p}. Temporal window where SITS are observed. p-dimensional process: T → R p . ith pixel from wavelength b. Number of temporal acquisitions associated with pixel i. Number of classes. Class value, in {1, . . . , C}. Discrete random variable associated with pixel i and representing its class membership. Prior probability of class c. Model parameters. Number of basis functions i.e. dimension of α. Determinant of the matrix K.

1: Symbols and Notations Symbol Meaning ∼ Distributed according to (a probability distribution). GP Gaussian Process: f ∼ GP(m, K), the function f is distributed as a GP with mean function m and covariance operator K. The function is indexed in time by t. N d (µ, Σ) d-variate Gaussian (Normal) distribution with mean vector µ and covariance matrix Σ .

Table 4 .

 4 2: Land cover classes and number of extracted SITS in each tile

	Class	T31TCJ T31TDN T31TGK
	Artificial areas
	Continuous urban fabric Discontinuous urban fabric Industrial or commercial units Road surfaces	10,000 10,000 10,000 10,000	8,292 10,000 10,000 9,906	959 10,000 10,000 3,664
	Agricultural areas
	Winter crops Summer crops	30,000 40,000	30,000 40,000	15,975 24,912
	Forest and semi-natural areas
	Meadow Orchards Vines Broad-leaved forest Coniferous forest Natural grasslands Woody moorlands	10,000 10,000 10,000 10,000 9,957 9,939 9,972	10,000 2,775 8,719 10,000 10,000 3,022 10,000	10,000 10,000 153 10,000 10,000 10,000 10,000
	Open spaces with little or no vegetation
	Bare rock Beaches, dunes and sand plains Glaciers and perpetual snow Water bodies	0 0 0 10,000	0 5,355 0 10,000	10,000 10,000 10,000 10,000
	Total	189,868 178,069	165,663

Table 4 .

 4 3: Averaged "mean F 1 score" (mean(%) ± standard deviation) computed on 3 tiles. MIMGP model is parametrized by a Fourier basis with 19 parameters for T31TCJ, and 41 for T31TDN and T31TGK tiles.

	QDA	SVM	RF	MIMGP
	T31TCJ 36.

2 ± 2.5 70.1 ± 2.1 72.2 ± 1.6 57.6 ± 2.0 T31TDN 30.5 ± 2.8 74.8 ± 1.9 77.6 ± 1.5 65.1 ± 1.0 T31TGK 38.9 ± 2.1 62.1 ± 1.6 63.9 ± 1.6 45.5 ± 2.7

Table 4 . 4 :

 44 Averaged normalized MAE (mean (%) ± Standard deviation) for each wavelength for the different tiles with the a Fourier basis of dimension J = 19.

	S2 Band	Method	T31TCJ	T31TDN	T31TGK
	Band 2	Whittaker 47.1 ± 0.5 123 ± 1.5 40.1 ± 1.2 MIMGP 46.9 ± 2.4 46.5 ± 1.3 36.3 ± 1.4
	Band 3	Whittaker 45.2 ± 0.5 114 ± 1.2 42.3 ± 1.3 MIMGP 43.9 ± 3.7 44.0 ± 1.6 39.6 ± 1.4
	Band 4	Whittaker 37.7 ± 0.4 94.2 ± 1.0 43.5 ± 1.3 MIMGP 42.7 ± 4.4 45.2 ± 2.2 41.3 ± 1.5
	Band 5	Whittaker 40.0 ± 0.5 99.9 ± 1.0 44.7 ± 1.3 MIMGP 45.8 ± 0.4 40.9 ± 2.0 41.3 ± 2.2
	Band 6	Whittaker 32.6 ± 0.4 83.2 ± 1.0 51.1 ± 1.3 MIMGP 38.6 ± 1.8 46.4 ± 0.9 46.9 ± 2.6
	Band 7	Whittaker 30.0 ± 0.3 73.5 ± 0.9 51.1 ± 1.2 MIMGP 24.7 ± 2.7 45.1 ± 0.9 41.8 ± 1.1
	Band 8	Whittaker 32.6 ± 0.4 71.8 ± 0.9 52.9 ± 1.2 MIMGP 33.1 ± 1.9 44.6 ± 2.3 42.9 ± 1.1
	Band 8A Whittaker 29.4 ± 0.3 69.3 ± 0.9 51.9 ± 1.2 MIMGP 33.2 ± 3.4 43.8 ± 1.5 42.6 ± 1.2

Band 11 Whittaker 36.1 ± 0.5 53.1 ± 0.7 46.9 ± 0.6 MIMGP 37.0 ± 2.0 38.3 ± 1.7 43.0 ± 1.3 Band 12 Whittaker 34.1 ± 0.4 46.0 ± 0.6 51.3 ± 0.7 MIMGP 39.5 ± 2.6 38.7 ± 3.3 46.6 ± 2.6

Table 5 .

 5 1: Land cover classes and number of extracted samples n c per class for each training and validation set.

	Class	n c
	Summer crops Winter crops Broad-leaved forest Continuous urban fabric Discontinuous urban fabric Industrial or commercial units 10,000 40,000 30,000 10,000 10,000 10,000 Meadow 10,000 Orchards 10,000 Road surfaces 10,000 Vines 10,000 Water bodies 10,000 Woody moorlands 9,972 Coniferous forest 9,957 Natural grasslands 9,939
	Total	189,868
	discussed in the MIGP framework by [41, Fig. 8, and Fig. 1 in the supp. mat.].

Table 5 .

 5 2: Mean F 1 score (mean(%) ± standard deviation) on the 10 independent data-sets. ± 12.88 30.7 ± 16.90 19.4 ± 14.68 15.4 ± 7.86 20.6 ± 8.46

	QDA	RF	SVM	MIGP	M2GP
	Summer crops Winter crops Broad-leaved forest Cont. urban fabric Discont. urban fabric Ind. or commercial units 31.3 ± 2.14 96.5 ± 0.27 91.6 ± 0.48 77.4 ± 3.91 39.8 ± 6.18 58.5 ± 1.39 Meadow 58.3 ± 4.14 Orchards 72.9 ± 4.05 Road surfaces 73.1 ± 1.92 Vines 71.1 ± 4.35 Water bodies 98.7 ± 0.35 Woody moorlands 23.9 ± 7.70 Coniferous forest 76.6 ± 7.24 Natural grasslands 29.8 Average F 1 score 70.5 ± 0.75	96.8 ± 0.45 94.0 ± 0.77 86.2 ± 2.35 58.0 ± 1.55 57.3 ± 3.44 40.2 ± 12.61 42.5 ± 3.17 54.5 ± 0.80 95.6 ± 0.81 90.0 ± 0.83 95.9 ± 0.44 93.9 ± 0.66 80.2 ± 0.83 92.2 ± 0.64 85.3 ± 2.63 75.7 ± 5.03 81.5 ± 3.10 55.9 ± 2.49 21.4 ± 3.49 30.9 ± 5.51 60.3 ± 1.35 48.3 ± 4.05 27.4 ± 0.92 38.4 ± 2.34 64.8 ± 2.94 63.0 ± 3.17 43.3 ± 3.80 55.0 ± 4.19 81.0 ± 2.64 76.4 ± 3.11 51.9 ± 5.46 77.6 ± 3.58 87.1 ± 1.87 78.7 ± 2.79 54.2 ± 5.79 75.0 ± 2.06 78.9 ± 6.86 78.5 ± 6.57 60.9 ± 7.61 71.7 ± 5.18 99.4 ± 0.08 99.3 ± 0.10 84.9 ± 5.38 96.8 ± 0.84 56.6 ± 3.50 56.1 ± 3.85 14.1 ± 5.52 10.6 ± 12.00 86.9 ± 2.76 87.0 ± 2.56 61.2 ± 5.41 82.4 ± 6.61 78.2 ± 1.17 75.2 ± 1.11 57.4 ± 1.04 70.1 ± 0.43
	5.8 Appendix -Proofs

Proof of Proposition 5.1 Let Y ∼ MGP p (m, K,

TABLE I :

 I Averaged mean F 1 score (mean(%) ± standard deviation) associated with u-QDA and QDA computed on three tiles.

	u-QDA	QDA
	T31TCJ T31TDN 21.3 ± 3.6 26.2 ± 3.1 T31TGK 29.2 ± 2.1	36.2 ± 2.5 30.5 ± 2.8 38.9 ± 2.1

TABLE II :

 II F 1 scores for tile T31TCJ. ± 11.8 53.6 ± 01.1 56.9 ± 01.2 27.8 ± 03.6 Discontinuous urban fabric 24.2 ± 05.6 53.7 ± 02.7 38.9 ± 10.1 51.2 ± 01.1 Industrial or commercial units 13.1 ± 03.7 55.5 ± 00.9 47.8 ± 03.8 27.7 ± 01.6 Road surfaces 36.7 ± 06.1 85.4 ± 01.7 78.1 ± 04.0 65.7 ± 03.2 Winter crops 40.7 ± 07.9 93.2 ± 00.7 93.4 ± 01.1 86.4 ± 02.0 Summer crops 66.5 ± 08.4 96.6 ± 00.3 95.5 ± 00.9 92.1 ± 01.0 Meadow 25.4 ± 07.7 63.4 ± 02.9 63.2 ± 03.7 53.7 ± 06.2 Orchards 48.8 ± 05.9 79.2 ± 02.7 76.4 ± 03.3 53.1 ± 04.4 Vines 21.6 ± 13.9 74.4 ± 06.6 78.7 ± 07.0 64.6 ± 09.7 Broad-leaved forest 58.6 ± 11.0 85.3 ± 02.4 84.6 ± 02.9 71.7 ± 03.6 Coniferous forest 28.4 ± 16.5 86.1 ± 02.5 86.4 ± 02.4 73.4 ± 03.6 Natural grasslands 06.0 ± 07.1 29.3 ± 15.3 22.5 ± 12.2 21.8 ± 13.0 Woody moorlands 25.9 ± 03.6 54.6 ± 03.5 55.4 ± 04.0 20.6 ± 05.9 Water bodies 87.7 ± 07.7 99.3 ± 00.0 99.2 ± 00.1 95.1 ± 01.9

	Class Continuous urban fabric 23.0 Sentinel-2 QDA MIMGP	RF	RF	SVM	MIMGP	SVM
	Continuous urban fabric Discontinuous urban fabric Industrial/commercial units Road surfaces Winter crops Summer crops Meadow Orchards Broad-leaved forest Coniferous forest Natural grasslands Woody moorlands Water bodies Vines

TABLE III :

 III F 1 scores for tile T31TDN. ± 11.6 43.1 ± 11.3 36.7 ± 14.3 39.3 ± 10.0 Woody moorlands 4.9 ± 4.7 69.2 ± 2.2 68.3 ± 3.2 49.8 ± 2.5 Beaches, dunes and sand plains 43.1 ± 16.5

	Sentinel-2	MIMGP		RF	SVM
	Continuous urban fabric Discontinuous urban fabric Winter crops Summer crops Broad-leaved forest Coniferous forest Bare rock	Industrial/commercial units Meadow Natural grasslands	Road surfaces Orchards Woody moorlands
	Class		QDA	RF	SVM	MIMGP
	Continuous urban fabric Discontinuous urban fabric Industrial or commercial units Road surfaces Winter crops Summer crops Meadow Orchards Vines Broad-leaved forest Coniferous forest Natural grasslands Water bodies	20.3 ± 7.7 14.0 ± 6.5 13.4 ± 2.2 34.9 ± 6.4 23.0 ± 5.0 26.5 ± 9.9 32.1 ± 3.0 24.4 ± 12.9 42.8 ± 13.9 7.6 ± 10.1 61.6 ± 7.6 24.3 80.6 ± 5.6 67.7 ± 2.3 64.7 ± 3.6 55.1 ± 2.6 83.3 ± 1.5 93.4 ± 0.7 95.5 ± 0.7 85.0 ± 2.1 64.1 ± 9.6 89.6 ± 2.1 86.0 ± 1.2 91.9 ± 0.3 84.1 ± 3.3 94.0 ± 1.6	66.6 ± 4.0 54.8 ± 11.0 48.3 ± 4.3 75.6 ± 2.9 91.4 ± 1.2 94.2 ± 1.0 85.4 ± 2.2 66.0 ± 20.1 88.4 ± 3.0 83.5 ± 2.0 91.2 ± 0.5 78.7 ± 4.4 92.3 ± 2.6	58.5 ± 2.8 51.3 ± 3.2 38.3 ± 1.6 65.2 ± 2.6 85.3 ± 1.6 92.9 ± 0.8 71.2 ± 1.8 41.8 ± 6.8 85.8 ± 2.0 70.9 ± 2.8 80.3 ± 1.5 59.7 ± 6.8 86.4 ± 2.6

Beaches, dunes and sand plains Glaciers and perpetual snow Water bodies

TABLE IV :

 IV F 1 scores for tile T31TGK. ± 18.3 54.6 ± 4.5 45.6 ± 11.0 27.7 ± 15.5 Beaches, dunes and sand plains 47.3 ± 14.3 56.8 ± 9.3 61.3 ± 10.7 33.8 ± 7.5 Glaciers and perpetual snow 76.6 ± 10.2 89.3 ± 7.3 90.9 ± 2.7 78.7 ± 10.8 Water bodies 14.5 ± 11.1 86.3 ± 2.4 83.3 ± 4.1 43.2 ± 13.1

	Class	QDA	RF	SVM	MIMGP
	Continuous urban fabric Discontinuous urban fabric 23.2 Industrial or commercial units 0.3 ± 0.9 27.0 ± 6.1 Road surfaces4 33.3 ± 12.2 63.2 ± 8.5 9.1 ± 3.6 55.5 ± 3.5 Winter crops 61.6 ± 9.2 91.3 ± 1.9 Summer crops 80.8 ± 4.7 94.3 ± 1.4 Meadow 25.2 ± 4.7 48.1 ± 2.6 Orchards 84.6 ± 3.6 91.5 ± 1.2 Broad-leaved forest 20.9 ± 15.8 79.9 ± 5.2 Coniferous forest 43.8 ± 4.6 79.4 ± 2.5 Natural grasslands 35.2 ± 5.6 49.3 ± 3.2 Woody moorlands 46.7 ± 5.3 50.8 ± 2.9 Bare rock 40.7 ± 7.0 72.3 ± 4.9	1.7 ± 1.0 50.3 ± 4.1 61.3 ± 9.0 88.0 ± 3.7 92.8 ± 2.2 53.1 ± 2.2 88.9 ± 2.0 78.5 ± 5.2 79.2 ± 2.6 45.5 ± 9.3 48.4 ± 5.4 73.2 ± 4.4	8.7 ± 4.2 28.2 ± 4.7 42.1 ± 5.8 72.4 ± 3.2 81.9 ± 3.2 19.2 ± 3.2 73.5 ± 2.3 56.0 ± 17.0 45.4 ± 4.6 37.9 ± 4.4 33.2 ± 2.3 46.8 ± 4.2

European Copernicus program, https://www.copernicus.eu/en

All products are available in the Theia Land Data Center: http://www.theia-land.fr/en/presentation/products.

Chapter 2 -Statistical modelling for time-series classification

Sentinel-2 products are available as a collection of elementary tiles of size 100×100 km, see https://sentinel.esa.int/web/ sentinel/missions/sentinel-2/data-products.

http://www.theia-land.fr/en/presentation/products.

Using "SuperImpose" application, see https://www.orfeo-toolbox.org/CookBook/Applications/app_Superimpose.html.

Since a stratification w.r.t polygons is done, and polygons size can vary in terms of area, the number of pixels may vary slightly from one experiment to another. The model was trained on a total average of 178,000 pixels and was evaluated on 178,000 pixels.

Per class F 1 scores are provided in the supplementary materials.

Full classification map and extracts for the three tiles are provided in the supplementary materials.

https://sentinel.esa.int/web/sentinel/missions/sentinel-2/data-products.

This is true when the period of observation is not too long, few years, otherwise the material property might be altered and its reflectance could vary.

# F r e e memory 120 f r e e ( < v o i d * > i p i v _ ) ; f r e e ( < v o i d * > i w o r k _ ) f r e e ( < v o i d * > r _ ) ; f r e e ( < v o i d * > c_ ) 122 f r e e ( < v o i d * > work_ ) ; f r e e ( < v o i d * > Af )

Acknowledgements / Remerciements

To start in English, I would like to thank all the members of the jury. Thank you Julien and Gabriele for accepting to review my manuscript, I really enjoyed reading your comments and feedback on my work. Thanks also to Marie and Lionel, it was a pleasure to answer your questions during my oral defence. Finally, thanks to Serge and Yannick, I enjoyed talking to you both during the three years of my thesis. The following is written in French.

Je tiens tout d'abord à remercier Stéphane et Mathieu pour m'avoir offert l'opportunité de faire cette thèse. C'était un plaisir d'avoir pu collaborer avec vous pendant trois ans et je pense que toutes vos exigences pour le moindre détail (chaque graphe, chaque mot, jusqu'au contenu du code), m'ont vraiment permis de clôturer sur ce manuscrit que j'espère être clair, concis et bien illustré. J'esp'ere que nous continuerons à l'avenir d'échanger sur les problématiques soulevées dans cette thèse.

Acknowledgements v

General Introduction

Acknowledgment

The authors thank the reviewers for their many constructive and helpful comments. The authors would like to thank S. Iovleff for his support and advices during the construction of the model. The authors would also like to thank Y. Tanguy for his help when using the CNES computational resources to run the experiments presented in this paper.

Acknowledgment

The authors would like to thank S. Iovleff for his support and advices during the design of the model. The authors would also like to thank Y. Tanguy for his help when using the CNES computational resources to run the experiments presented in this paper.

This work is supported by the French National Research Agency in the framework of the Investissements dAvenir program (ANR-15-IDEX-02) and by the Centre National dEtudes Spatiales (CNES).

Chapter 4

Joint supervised classification and reconstruction of irregularly sampled satellite image times series

The following content has been published: A. Constantin, M. Fauvel, and S. Girard, "Joint supervised classification and reconstruction of irregularly sampled satellite image times series", IEEE Transactions on Geoscience and Remote Sensing, pp. 1-13, 2021, To appear. doi: 10.1109/TGRS.2021.3076667 Appendices are reported in Section 4.9 and supplementary materials in Appendix B.

Numerical implementation

The computation of the MLE is implemented as detailed in Algorithm 3 using the results of Proposition 5.2. To deal with the idenfiability issue mentioned in Paragraph 5.3.1, A c A c is normalised by η c such that A c A c F = 1 (where • F denotes the Frobenius norm) and each covariance matrix Σ c,i (θ c ) is modified accordingly so that the likelihood remains unaffected (step (d) of Algorithm 3). The gradient step (e) is performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, see [START_REF] Zhu | Algorithm 778: L-BFGS-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF]. More specifically, the L-BFGS-B version is used which allows for box and positivity constraints. As described in [START_REF] Zhu | Algorithm 778: L-BFGS-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF], the gradient step is obtained by line search and the algorithm stops when: the objective function (i.e. the likelihood) does not change significantly, the (infinite) norm of the projected gradient is sufficiently small or when the maximum number of iterations is reached. Since the objective function is not convex, the optimization process is sensitive to the initialization. In practice, multiple random restarts are used and the best solution is retained. Let us highlight that, in practice, steps (a)-(e) are computed for all classes in parallel since the model parameters are decoupled w.r.t. the classes.

Algorithm 3: Computation of MLE of model parameters.

Input : Sample (Y i, , Z i ) ∈ R p×q i × {1, . . . , C}, i = 1, . . . , n and initialization (θ 1 , . . . , θ C ).

Output: MLE αc , A c A c , θc , c = 1, . . . , C. 7 (e) Update θ c with a gradient step using (5.14);

The numerical complexity of one iteration for all classes of Algorithm 3 is O(n(q 3 ∞ + p 3 + J 3 )) where n is the sample size and q ∞ = max{q i , i = 1, . . . , n}. The computation of the MLE thus scales linearly w.r.t. n. In constrast, the cost associated with standard classfication methods based on Gaussian processes is O((C +1)n 3 ) [START_REF] Williams | Gaussian processes for machine learning[END_REF]Algorithm 3.3]. Here, the computation of the MLE only relies on the inversion of p × p and q i × q i matrices whose sizes do not depend on the sample size.

Let us note that Algorithm 3 can be interpreted as an extension of the so-called Flip-flop method introduced independently by [START_REF] Mardia | Spatial-temporal analysis of multivariate environmental monitoring data[END_REF][START_REF] Dutilleul | The MLE algorithm for the matrix normal distribution[END_REF]. This latter method in an iterative way to compute the MLE associated with the matrixvariate normal distribution. As such, it is limited to the situation where q 1 = q 2 = • • • = q n which only occurs when all Gaussian processes are observed on a common grid. Identifiability issues are discussed in [START_REF] Srivastava | Models with a kronecker product covariance structure: Estimation and testing[END_REF] and the method is extended to higher order tensor distributions in [START_REF] Manceur | Maximum likelihood estimation for the tensor normal distribution: Algorithm, minimum sample size, and empirical bias and dispersion[END_REF]. Applications of matrix-variate normal distribution are found in different contexts such as electro-encephalography [START_REF] Spinnato | Finding EEG space-time-scale localized features using matrix-based penalized discriminant analysis[END_REF] or remote sensing [START_REF] Glanz | An expectation-maximization algorithm for the matrix normal distribution with an application in remote sensing[END_REF].

Finally, all the above estimations procedures have been implemented in Python using the Scikit-Lean API, see [START_REF] Buitinck | API design for machine learning software: Experiences from the scikit-learn project[END_REF]. The Fourier basis {ϕ 1 , . . . , ϕ J } was chosen to estimate the mean function (see [START_REF] Constantin | Joint supervised classification and reconstruction of irregularly sampled satellite image times series[END_REF] for other bases), while the family of symmetric positive-definite kernels was selected among the Kernels class in the Scikit-Learn library.

Validation on simulated data

The performance of the inference procedure associated with the M2GP model is illustrated on simulated data. 1 The simulated model is described in Paragraph 5.5.1. First, the influence of the dependence between coordinates as well as the influence of the number of observation times are investigated in Paragraph 5.5.2. Second; the consequences on the classification and imputation accuracy are discussed in Paragraph 5.5.3.

Experimental design

A binary classification problem is considered. Two classes are simulated from a 10-dimensional M2GP model on T = [0, 1] with 1, 000 samples per class leading to n = 2, 000 and p = 10. Mean functions are generated following (5.7) with a Fourier basis of size J = 11. Coefficients α c,b, j are simulated independently from a N 1 (0, 0.02) 1 The code and a notebook are available at https://gitlab.inria.fr/aconstan/mixture-of-multivariate-gaussianprocesses-for-classification-of-irregularly-sampled-satellite-image-time-series.

About the estimation of the mean

This appendix describes the estimation of the mean. Section A.1 presents the study of the two matrix derivatives to assess the optimum found in Chapter 5. Section A.2 highlights computational issues encountered with random signals (irregularly and unevenly sampled).

A.1 On the broad and the narrow definitions of matrix derivative

This section discusses the two matrix derivatives defined in [START_REF] Magnus | On the concept of matrix derivative[END_REF] for the computation of optimal mean in Chapter 5. Let f : M p,q (R) → R be a function which returns a scalar from a p×q matrix X ∈ M p,q (R), then the broad definition is defined as:

and the narrow definition as:

The first equation defines the derivative as a matrix of size p × q and the second definition returns a pq-dimensional vector.

It has been shown in [START_REF] Magnus | On the concept of matrix derivative[END_REF] that the so-called narrow definition of the matrix derivative generalizes better the derivative on vectors. However, what follows proves that the two definitions give us the same result for computing the partial derivative with respect to the mean coefficient α c ∈ M p,J (R) in Proposition 5.2.

Proof. Firstly, recall that the differential w.r.t. α c is given by the Proposition 5.2, on one hand, as:

And, on the second hand, using Kronecker product properties [START_REF] Schott | Matrix Analysis for Statistics[END_REF]Theorem 8.12] gives:

Then:

1. Using the broad definition (A.1), and the matrix inner product (

Finally,

2. Using the narrow definition (A.2), the derivative w.r.t. vec(α c ) from (A.4) is:

Then we have

and yields to the same result as (A.5).

As a conclusion, considering one definition or the other returns the same optima in the maximum likelihood estimation.

A.2 On the conditioning of the design matrix

In what follow, we discuss the problem of time-series reconstruction. Whereas the reconstruction of regularly sampled time-series is well-known (for example using the Fourier decomposition, it is known that the frequency of the associated basis must not be too high w.r.t. the Shannon criteria), the theoretical study of irregularly sampled time-series involves fewer theoretical criteria to use properly the two models.

We defined in Appendix B (Section I) a necessary condition based on the injectivity of the design matrix B to compute the optimal mean. The matrix B{Σ -1 }B is invertible if B is injective. This section goes further by defining a sufficient criteria to obtain a numerically stable inversion in Section A.2.1. This criteria involves the matrix B B. Indeed if matrix B is injective, numerical issues occur when a temporal window is not observed by any sample, this is ilustrated in Section A.2.2.

A.2.1 Theoretical study

i=1 a set of n i.i.d. irregularly sampled time-series from a Gaussian process within the same class c ∈ {1, . . . , C}. For any given i, following Appendix B, E(y

where B c is the design matrix of size T c × J, with T c the total number of unique time stamps within S c , and N i of size T i × T c a matrix which selects the observed time stamps in y i .

From (4.8) (or (A.5)), αc is obtained by inverting the matrix

, with the previous notations and M c = i|z i =c N i Σ i -1 N i a non-singular matrix as presented in Appendix B, equation (2) (the number corresponds to the one of the included supplementary materials file). What follows aim to study the conditioning of this matrix to be numerically stable. Proof of Proposition A.1. The proof is in two steps, firstly the upper bound of the largest eigenvalue is presented, then the lower bound for the smallest eigenvalue. Recall that the largest and the smallest eigenvalues correpond to

Model settings and estimation

The MIMGP model is learned using a Fourier base with 19 basis functions, a Matérn kernel (smoother than RBF [START_REF] Genton | Classes of kernels for machine learning: A statistics perspective[END_REF] with parameter ν = 2.5) and an additive white noise. The estimated means for one run are reported in Figures A.2 The first class and the second class are well estimated for a size of N sufficiently small. Indeed, these two classes are generated according to the model. The third class is more complex, the reconstruction is poor as the basis functions are not chosen correctly. However, the three classes present a huge mis-estimation of mean for large unobserved windows which can sometimes shift the estimated mean where the data are observed (see class 3 on Figure A.4).

The same experience holds for a fixed unobserved temporal window and an increasing number of basis functions.

Condition number within class 1

The condition numbers for the quantities in Proposition A. 1 The inequality in Proposition A.1 is verified. No lower bound has been found but the upper bound becomes broad with the condition number of matrix B c B c , where B c is the design matrix within class c.

Despite no improvements on the classification accuracy, a regularized approach can improve the conditioning of the matrix to be inverted (see [START_REF] Tabeart | The conditioning of least-squares problems in variational data assimilation: The conditioning of least squares problems in variational data assimilation[END_REF]Section 3.1]). Nevertheless, a regularization on the mean has to take into account the complete temporal window of interest (T ) with an increasing number of unobserved time-stamps as the number of basis functions increases to avoid mis-estimation.

Supplementary Materials of Chapter 4

This appendix describes the supplementary materials of [START_REF] Constantin | Joint supervised classification and reconstruction of irregularly sampled satellite image times series[END_REF] (compiled before inclusion, pagination does not follow the one of this manuscript neither the references). They contain additional details related to the main paper. The identifiability of the model depends on the ability to compute α b,c . In (1) below, we remind that the optimal parameter is computed using the inverse of the matrix G c , defined as the sum of terms involving the design and precision matrices:

Supplementary Materials

The goal of this Section is to establish a necessary and sufficient condition for G c to be non-singular. Let T c be the total number of unique temporal acquisitions for a given class c. Let us also introduce the T i × T c matrix N i (with T i ≤ T c for all i such as z i = c) and the global design T c × J matrix B c such that B i = N i B c . The matrix N i is composed of ones and zeros to select the observed samples in signal i from the global design matrix. Then, the matrix G c can be rewritten as

where we set

From the definition of M c , and since Σ i (θ b,c ) -1 is a positive definite matrix, v * ⊤ M c v * is a sum of non-negative terms. Consequently, this entails that Section C.1 presents additional details on implementation and identify key elements to optimize and Section C.2 presents the optimization done.

C.1 Implementation of the code

Firstly, the code has been designed to inherit scikit-learn API [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF]. To do so, some classical functions are provided as "score" or "fit" and "predict". It allows the use of scikit-learn's functions as GridSearch.

The model is a python class. The use of the model is standard. However the main problem with python is the use of "for" loops. The language is well optimized in a vectorized approach of numerical objects. In out context two loops are identified as heavy from a computational point of view:

1. A loop on each sample for each class (and each spectral band for MIMGP) is hardly vectorized.

2. A loop on each sample when computing the log-marginal likelihood and the prediction of posterior probabilities functions.

Additionally the code has been implemented in the CNES computational resources in order to use both the HPC clusters and the simplified Sentinel-2 data requests.

C.2 Optimization of the code

An important part of this work is to optimize the code to compute the large number of Sentinel-2 SITS, the two optimizations are described briefly.

The first loop has been easily parallelized using multi-processing techniques. It consists of p × C independent loops for MIMGP and C independent loops for M2GP, where C is the number of classes and p the number of spectral bands.

The second loop does not allow parallelization. Indeed the loop computes the sum of the likelihood and each element of the sum is heavy to compute: computation of kernel operators, inversion of matrices, etc. The latter part also has to be computed to predict posterior probabilities. The latter one uses cython, a compiled language, to process the functions. Algorithms 4 and 5 in Section C.3 are an example of the cython code to compute α. This code is then compiled using the cython compiler.