Mireille Laigle
email: laigle@geoazur.unice.fr

Hans Agurto-Detzel

Anne Bécel

Milton Boucard

Caroline Chalumeau

Philippe Charvis

Jean-Xavier Dessa

Audrey Galve

Maria-José Hernandez

Sara Hussni

Frauke Klingelhoefer

Heidrun Kopp

Muriel Laurencin

Jean-Frédéric Lebrun

Boris Marcaillou

François Michaud

Michele Paulatto

Alessandra Ribodetti

Maria Sachpazi

Laure Schenini

Imaging the megathrust in subduction zones: lessons from Greece, Ecuador and the Lesser Antilles

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Understanding the physical parameters and processes that control the seismogenic behavior of subduction zones megathrust faults remains one of the outstanding challenges in Earth Sciences.

Here we present important results from several large seismic experiments aimed at addressing this question. These experiments focused on the three subduction zones off Greece, the Lesser Antilles islands, and Ecuador, with different convergence rates and seismic activities. Surveys included multibeam bathymetry, multichannel reflection seismic (MCS) and wide-angle seismic (WAS) acquisitions over the forearc domain, as well as teleseismic receiver-functions and local earthquakes monitoring with temporary deployments of seismological networks.

Our results demonstrate the needs of both dense and extensive geophysical investigations.

In the central Lesser Antilles subduction zone, the interplate has been imaged down to the backstop at 12-15 km depth over the 350-km-long Antigua to Martinique islands segment. The outer forearc crust is strongly faulted in response to the two subducting Tiburon and Barracuda ridges (SISMANTILLES1-and-2 surveys). Two WAS profiles constrained the deeper geometry of the interplate down to the forearc Moho located at 28 km depth (TRAIL survey). The OBS networks deployed over several months (OBSANTILLES and OBSISMER surveys) revealed mantle wedge supraslab earthquakes and M4-5 possible repeaters with flat-trust mechanisms. The joint active-source/local earthquake seismic tomography let us to unveil the Vp and Vp/Vs heterogeneity along the slab surface and derive unprecedented constraints on multi-stage fluid release from subducting slow-spread oceanic lithosphere. Farther northwest, where the convergence obliquity strongly increases, we constrained the geometry of the interplate down to the forearc Moho at 25 km depth. Strain partitioning localizes on inherited major structures within the forearc domain, like the left-lateral partitioning system of the Anegada Passage and the 850-km-long Bunce fault, located along the backstop (ANTITHESIS survey).

On the southwestern Hellenic subduction zone, MCS and WAS acquisitions highlight the existence of an outer forearc crust beneath the forearc Matapan Trough, but its highly complex structure prevented us to image the interplate (ULYSSE survey). Acquisition by the R/V Marcus Langseth with its 8-km-long streamer finally made it possible (SISMED survey). Dense receiver-function acquisition on a 300-km-long mobile seismic network constrained the 3D geometry of the slab top underneath central Greece. This imaging revealed that the subducting oceanic crust and backstop updip limit are segmented by 9 trench-normal subvertical faults, seismically active at intermediate depths and possibly of inherited origin (THALES WAS RIGHT survey).

South of the 1906 M8.8 Ecuador-Columbia rupture area, the April 2016 Mw7.8 Pedernales subduction earthquake and its ensuing postseismic phase revealed a combination of seismic/aseismic slip behavior. Fluid-enriched parts of the megathrust fault and structural margin segmentation are hypothesized to play a major role in controlling slip behavior but direct observations are still lacking. Previous MCS acquisitions revealed very locally a fluid-rich subduction channel along with severe damage effect of the forearc margin due to seamounts subduction (SISTEUR survey). Forthcoming 3D seismic acquisition along this segment will examine the impact of the along-strike and along-dip variations of the physical properties and fluid content on the slip mode (HIPER survey).

Powered by TCPDF (www.tcpdf.org)


Student Member, IEEE, Tales Imbiriba Ricardo Augusto Borsoi

Senior Member, IEEE José Carlos Moreira Bermudez

J C M Bermudez
email: j.bermudez@ieee.org.

C Richard
email: cedric.richard@unice.fr

Fast Unmixing and Change Detection in Multitemporal Hyperspectral Data

Keywords: Hyperspectral data, multitemporal, spectral unmixing, endmember variability, MESMA

Multitemporal spectral unmixing (SU) is a powerful tool to process hyperspectral image (HI) sequences due to its ability to reveal the evolution of materials over time and space in a scene. However, significant spectral variability is often observed between collection of images due to variations in acquisition or seasonal conditions. This characteristic has to be considered in the design of SU algorithms. Because of its good performance, the multiple endmember spectral mixture analysis algorithm (MESMA) has been recently used to perform SU in multitemporal scenarios arising in several practical applications. However, MESMA does not consider the relationship between the different HIs, and its computational complexity is extremely high for large spectral libraries. In this work, we propose an efficient multitemporal SU method that exploits the high temporal correlation between the abundances to provide more accurate results at a lower computational complexity. We propose to solve the multitemporal SU problem by separately addressing the endmember selection and the abundance estimation problems. This leads to a simpler solution without sacrificing the accuracy of the results. We also propose a strategy to detect and address abrupt abundance variations in time. Theoretical results demonstrate how the proposed method compares to MESMA in terms of quality, and how effective it is in detecting abundance changes. This analysis provides valuable insight into the conditions under which the algorithm succeeds. Simulation results show that the proposed method achieves state-of-the-art performance at a smaller computational cost.

I. INTRODUCTION

Hyperspectral images (HI) have become a central tool in an increasing number of applications due to their high spectral resolution, which offers important information about the materials in a scene [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF]. However, inherent limitations of imaging devices and large sensor-to-target distances typical of many applications, such as remote sensing, lead to HIs with low spatial resolution [START_REF] Shaw | Spectral imaging for remote sensing[END_REF]. Hence, each pixel in an HI is usually a mixture of spectral signatures of different pure materials, also called endmembers (EM) [START_REF] Keshava | Spectral unmixing[END_REF]. Spectral unmixing (SU) aims to decompose an HI into a collection of endmembers and their corresponding fractional abundances, thus revealing important information on the distribution of the materials in the scene [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF].

The Linear Mixing Model (LMM) is the simplest and most widely used to represent the interaction between light and the materials in a given scene [START_REF] Keshava | Spectral unmixing[END_REF]. The LMM represents the reflectance of each pixel as a convex combination of the endmembers in the scene. The combination coefficients can then be interpreted as the fractional abundances contributed to each pixel by the EMs. The simplest form of LMM models each material in the whole scene by a single endmember. Although allowing fast and simple SU strategies, such a model fails to account for the important phenomenon of endmember variability, observed in most practical scenes [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF]- [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF].

Spectral variability can be caused by many factors including, for instance, atmospheric, illumination and seasonal variations, and can propagate significant abundance and endmember estimation errors throughout the unmixing process [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF], [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF]. This motivated the use of more elaborated mixing models and algorithms that explicitly address endmember variability in SU [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF], [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF]. Several parametric and non-parametric models have been recently proposed which account for EM spectra variation within a single HI. These include Beta or mixture of Gaussian distributions [START_REF] Du | Spatial and spectral unmixing using the beta compositional model[END_REF], [START_REF] Zhou | A Gaussian mixture model representation of endmember variability in hyperspectral unmixing[END_REF], additive [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF] and multiplicative scaling factors [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF]- [START_REF]A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF], combinations of an uniform scaling and an additive variability dictionary [START_REF] Hong | An augmented linear mixing model to address spectral variability for hyperspectral unmixing[END_REF], reparameterizations using deep neural networks [START_REF] Borsoi | Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF], and low-rank tensor representations [START_REF] Imbiriba | Low-rank tensor modeling for hyperspectral unmixing accounting for spectral variability[END_REF]. Although these models were able to produce promising results, the most prominent strategy to deal with spectral variability in SU is to represent EMs as sets of spectral signatures, called spectral libraries or bundles [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF]. The spectral libraries are usually constructed a priori from laboratory or in situ measurements. They contain variants of the spectral signatures of each material, to better represent various acquisition conditions or physicochemical compositions. SU can then be formulated as a problem of selecting from the library the subset of spectral signatures that best represents the observed HI. This usually entails the use of either sparse unmixing [START_REF] Iordache | Sparse unmixing of hyperspectral data[END_REF], [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF] or Multiple Endmember Spectral Mixture Analysis (MESMA) [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF] algorithms. The MESMA algorithm is still the leading algorithm in practice due to its simplicity and interpretability, and because it provides good results if the spectral library adequately expresses the spectra contained in the scene [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF], [START_REF] Borsoi | Deep generative models for library augmentation in multiple endmember spectral mixture analysis[END_REF]. This allowed MESMA to be applied to many different environments and scenarios [7, p.1607].

A. Multitemporal SU and change detection

More recently, multitemporal SU has become a topic of great interest due to its ability to leverage temporal information in HI sequences for monitoring the evolution of the different materials and their distribution in the scene [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF]- [START_REF] Goenaga | Unmixing analysis of a time series of hyperion images over the guánica dry forest in puerto rico[END_REF]. Spectral variability becomes a more critical issue in multitemporal SU, when compared to EM variations within a single HI, because images are acquired at different instants [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF], [START_REF] Borsoi | Super-resolution for hyperspectral and multispectral image fusion accounting for seasonal spectral variability[END_REF]. Some works have proposed to extend parametric EM models developed for single-image SU to the multitemporal case [START_REF] Henrot | Dynamical spectral unmixing of multitemporal hyperspectral images[END_REF]- [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF]. These works considered either uniform [START_REF] Henrot | Dynamical spectral unmixing of multitemporal hyperspectral images[END_REF] or bandwise [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF] scaling variations of reference EMs, or additive temporal perturbations over a mean EM matrix [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]- [START_REF] Thouvenin | A hierarchical bayesian model accounting for endmember variability and abrupt spectral changes to unmix multitemporal hyperspectral images[END_REF]. Such works attempt to estimate the EMs from the data by employing batch (centralized [START_REF] Thouvenin | A hierarchical bayesian model accounting for endmember variability and abrupt spectral changes to unmix multitemporal hyperspectral images[END_REF], [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF] or distributed [START_REF] Sigurdsson | Sparse distributed multitemporal hyperspectral unmixing[END_REF]) or online [START_REF] Henrot | Dynamical spectral unmixing of multitemporal hyperspectral images[END_REF], [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF] processing strategies. Nonetheless, the multitemporal extension of MESMA remains the most frequently used solution in practical applications [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF]- [START_REF] Somers | Multi-temporal hyperspectral mixture analysis and feature selection for invasive species mapping in rainforests[END_REF], [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF], and naturally addresses both the spatial and temporal variations of EM spectra.

A closely related problem to SU consists of detecting and monitoring changes in material composition across an HI sequence [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF]. While many methods are focused on detecting only abrupt changes [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF], [START_REF] Borsoi | Online graph-based change point detection in multiband image sequences[END_REF], having detailed subpixel abundance evolution from which both subtle and abrupt variations can be discerned is important for many applications [START_REF] Goenaga | Unmixing analysis of a time series of hyperion images over the guánica dry forest in puerto rico[END_REF], [START_REF] Chakravortty | A technique for subpixel analysis of dynamic mangrove ecosystems with time-series hyperspectral image data[END_REF]. This motivated the development of unmixing-based change detection strategies based on, e.g., linear [START_REF] Ertürk | Informative change detection by unmixing for hyperspectral images[END_REF], sparse [START_REF] Ertürk | Sparse unmixing-based change detection for multitemporal hyperspectral images[END_REF] and locally adaptive [START_REF] Goenaga | Unmixing analysis of a time series of hyperion images over the guánica dry forest in puerto rico[END_REF] SU, where the estimated abundance results can be analyzed to detect abrupt changes [START_REF] Guo | Change detection for hyperspectral images via convolutional sparse analysis and temporal spectral unmixing[END_REF]. The ability of MESMA to produce good quality abundance estimates while addressing spatio-temporal EM variability makes it an important tool for this task. It is also important to note that learning-based approaches have recently shown good performance in classification [START_REF] Hong | Graph convolutional networks for hyperspectral image classification[END_REF], [START_REF] Hong | Invariant attribute profiles: A spatial-frequency joint feature extractor for hyperspectral image classification[END_REF] unmixing [START_REF] Hong | Endmember-guided unmixing network (EGU-Net): A general deep learning framework for self-supervised hyperspectral unmixing[END_REF], [START_REF] Li | Model-based deep autoencoder networks for nonlinear hyperspectral unmixing[END_REF] and related tasks in hyperspectral imaging [START_REF] Rasti | Feature extraction for hyperspectral imagery: The evolution from shallow to deep: Overview and toolbox[END_REF], [START_REF] Zhu | Deep learning in remote sensing: A comprehensive review and list of resources[END_REF]. Such approaches usually leverage the capability of deep neural networks to improve the quality of the results. However, the combination of simplicity, interpretability and robustness that underlies MESMA and related library-based approaches still makes them appealing, especially due to their potential for out-of-the-box solutions.

Nevertheless, the advantages offered by MESMA do not come without compromises. Besides depending on a spectral library, a significant drawback of MESMA lies in its high computational complexity, which increases very quickly with the size of the libraries and with the number of materials. This goes against the need for online processing of large amounts of hyperspectral data [START_REF] Ma | Remote sensing big data computing: Challenges and opportunities[END_REF], [START_REF] Chi | Big data for remote sensing: challenges and opportunities[END_REF], and motivates the search for new algorithms that are both efficient and accurate. An alternating angle minimization (AAM) approach has been recently proposed to provide an approximate but accurate solution to the MESMA problem with a significant reduction in computational complexity for large libraries [START_REF] Heylen | Hyperspectral unmixing with endmember variability via alternating angle minimization[END_REF]. Despite its merits, the AAM algorithm does not scale well with the number of materials in the scene, and thus may also lead to a large computational complexity in practical scenarios. Moreover, AAM and other existing multitemporal MESMA methods do not exploit the temporal correlation between the abundance maps, nor do they account for abrupt abundance variations.

B. Contributions and organization

In this paper, we propose a fast multitemporal SU algorithm, named FM-MESMA. Differently from previous works such as MESMA and AAM [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF], [START_REF] Heylen | Hyperspectral unmixing with endmember variability via alternating angle minimization[END_REF], the proposed method explores the abundance temporal information to deliver high-quality estimates at a lower complexity. Moreover, unlike previous change detection methods based on SU [START_REF] Goenaga | Unmixing analysis of a time series of hyperion images over the guánica dry forest in puerto rico[END_REF], [START_REF] Ertürk | Informative change detection by unmixing for hyperspectral images[END_REF], [START_REF] Guo | Change detection for hyperspectral images via convolutional sparse analysis and temporal spectral unmixing[END_REF], we integrate change detection in a specific EM selection stage of the algorithm, which allows us to account for EM variability while maintaining a small computational cost. In contrast to deep learning-based frameworks, the proposed method yields a simple, low-cost and robust solution which comes with rigorous theoretical guarantees. The main contributions of this paper are: a) We explicitly characterize slow and abrupt abundance variations in the multitemporal mixing model. This allows us to exploit abundance temporal correlation to propose an efficient and accurate SU algorithm. Differently from previous works, this is achieved by performing SU in two separate tasks: EM selection and abundance estimation, each of which can be solved more efficiently. This significantly reduces the complexity of SU compared to MESMA or AAM, with little impact on the results. b) We propose a methodology to detect pixels that undergo abrupt abundance changes based on the results of the EM selection task. Such pixels, which would otherwise degrade the performance of the proposed method, are then handled separately using a more sophisticated strategy. Unlike typical SU-based change detection, the changed pixels are identified by taking spectral variability into account but without solving the full SU problem, leading to a lower complexity. The resulting algorithm is robust and interpretable, and has only a single tuning parameter. c) We derive theoretical guarantees concerning the performance of FM-MESMA, both in the presence and in the absence of changes. Specifically, we first show under which conditions the proposed method is guaranteed to recover the correct EMs from the library. Then, we derive conditions under which the FM-MESMA correctly detects abrupt abundance changes. These theoretical results provide insight into which conditions are necessary for the proposed method to reach an accurate result. d) We provide an analysis of the computational complexity of FM-MESMA, MESMA and AAM, and show how they scale with the number of spectral bands, the size of the spectral libraries and the number of EMs. Simulation results with synthetic and real data illustrate the performance of FM-MESMA when compared to MESMA, AAM, and approaches that do not rely on spectral libraries.

The paper is organized as follows. Section II reviews the linear mixing model, MESMA and its multitemporal extensions. Section III presents the multitemporal mixing model and our algorithm. Section IV provides theoretical guarantees for the reconstruction accuracy and robustness of FM-MESMA, and an analysis of its computational complexity. Section V presents simulation results and comparisons. Conclusions are presented in Section VI.

II. SPECTRAL UNMIXING WITH MESMA

The basic idea behind the MESMA algorithm is to find the EMs and the fractional abundances that best represent each pixel with the LMM. The LMM represents each L-band pixel y n ∈ R L of an HI, for n = 1, . . . , N , as a convex combination of the spectral signatures of P endmembers:

y n = M a n + e n , s.t. 1 a n = 1 and a n ≥ 0 (1)
where matrix M ∈ R L×P contains the spectral signatures of the EMs as its columns m p , p = 1, . . . , P , a n is the abundance vector, and e n is an additive noise term.

Most SU strategies use a single spectral signature to represent each material in the scene. This can lead to significant abundance estimation errors in the presence of spectral variability. MESMA, on the other hand, considers P spectral libraries known a priori, one for each endmember, defined as:

M p = {m p,1 , . . . , m p,Cp }, m p,j ∈ R L , p = 1, . . . , P (2)
where C p is the number of available variations of spectral signatures of the p-th endmember. Then, for each pixel n, only one endmember is selected from each set M p to compose the endmember matrix M n for that pixel.

We define the set of all possible endmember matrices that can be composed this way as:

M = m 1 , . . . , m P : m p ∈ M p , p = 1, . . . , P (3)
Assuming M known, the MESMA SU problem corresponds to the search for the EM matrix (also called EM model) in M that best represents each pixel in the scene. This translates into the following optimization problem for the n-th pixel:

min M n ∈M min an y n -M n a n , s.t. a n ≥ 0, 1 a n = 1. (4)
Despite its widespread use and good performance in practical scenarios, the computational cost of MESMA is extremely high. Since solving (4) amounts to perform SU for every possible matrix M n extracted from M, its computational complexity scales with the product of the sizes of libraries M p as it consists of solving P p=1 |M p | FCLS (Fully Constrained Least Squares) problems [START_REF] Heylen | Hyperspectral unmixing with endmember variability via alternating angle minimization[END_REF].

Several works have attempted to circumvent this limitation by seeking approximate solutions to [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. The first approaches consisted of solving (4) for matrices M n randomly chosen from M until obtaining a reconstruction error below a threshold and well distributed across all spectral bands [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF]. Another approach ignores both constraints in (4) and performs unconstrained least squares for every possible M n ∈ M, and then selects the EM model resulting in the smallest reconstruction error without any negative abundances [START_REF] Combe | Analysis of OMEGA/Mars express data hyperspectral data using a multipleendmember linear spectral unmixing model (MELSUM): Methodology and first results[END_REF].

Although these approaches are simple, they are not guaranteed to achieve a good accuracy at a reduced processing cost, as only a relatively small subset of M could be tested. Also, ignoring the abundance constraints can make the results more sensitive to noise. Recent strategies attempt to provide lowcomplexity alternatives to MESMA with minimal impact on unmixing results. For instance, the approach in [START_REF] Heylen | Hyperspectral unmixing with endmember variability via alternating angle minimization[END_REF] employs an angle minimization strategy. A significant reduction of the computational complexity is obtained for P small and C p possibly very large, with unmixing accuracy similar to MESMA. Another approach formulates problem (4) as a mixed-integer optimization problem in order to benefit from advanced software packages [START_REF] Mhenni | Spectral unmixing with sparsity and structuring constraints[END_REF]. However, none of these works consider the multitemporal formulation of MESMA.

MESMA has recently been applied to multitemporal SU problems such as monitoring of rainforests [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF], [START_REF] Somers | Multi-temporal hyperspectral mixture analysis and feature selection for invasive species mapping in rainforests[END_REF], [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF] and shrublands [START_REF] Lippitt | Multidate MESMA for monitoring vegetation growth forms in southern california shrublands[END_REF]. Some of these methods improve MESMA performance by employing strategies such as band selection and weighting [START_REF] Somers | Invasive species mapping in hawaiian rainforests using multi-temporal hyperion spaceborne imaging spectroscopy[END_REF], [START_REF] Somers | Multi-temporal hyperspectral mixture analysis and feature selection for invasive species mapping in rainforests[END_REF], or library construction from multiple time instants [START_REF] Dudley | A multi-temporal spectral library approach for mapping vegetation species across spatial and temporal phenological gradients[END_REF]. However, they do not explicitly explore temporal correlation between the abundance maps at adjacent time instants. In the following, we propose a model for the evolution of the abundance maps over time, accounting for both small and large variations. This will allow us to devise an efficient algorithm to address multitemporal SU problems.

III. FAST MULTITEMPORAL MESMA

The multitemporal SU problem can be introduced with generic terms as follows: given a sequence of image pixels {y t,n }, for t = 1, . . . , T time instants and n = 1, . . . , N pixels, and a spectral library M as defined in (3), estimate the corresponding fractional abundances {a t,n } and EM models M t,n . For simplicity, we shall assume that all the images are spatially aligned, such that for each n ∈ {1, . . . , N } the pixels y t,n , t = 1, . . . , T refer to the same spatial location. A simple solution to this problem would be to directly apply the techniques discussed in Section II to each pixel individually. This, however, ignores important temporal information contained in the image sequence which can be used in order to devise an efficient algorithm.

We propose to model the evolution of the abundance maps by considering its changes to be composed of a small additive signal, and of large sparse changes. That is, the observation model represents the pixels at time t and t + 1 as follows:

y t,n = M t,n a t,n + e t,n , (5a)
y t+1,n = M t+1,n a t,n + δ t,n + s t,n + e t+1,n (5b) = M t+1,n a t+1,n + e t+1,n ,
where M t,n is the (true) endmember matrix for pixel n at time instant t, and e t,n is an additive noise vector. Changes taking place in the abundances between time t and time t + 1 are modeled as a combination of a small magnitude term δ t,n and a spatially sparse, high magnitude term s t,n , which represents abrupt variations taking place in a small number of image pixels. Note that, although model [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF] dictates the relationship between a single pair of images, the extension to multiple images is trivial and is thus omitted here for simplicity.

It turns out that the structure outlined in the model ( 5) can be explored in order to devise an efficient MESMA-based SU algorithm. We propose to use an online strategy to estimate the abundances and the EM matrices at time instant t+1 based on an estimate a t,n of the abundances at time instant t. The procedure is a two-step one:

1) Considering a t,n ≡ a t,n , estimate M t+1,n , δ t,n and s t,n that best represent pixel y t+1,n in the model (5); 2) Set a t+1,n = a t,n + δ t,n + s t,n and repeat for the next image.

Taking into consideration the prior information stated about the properties of δ t,n and s t,n , and the knowledge of the spectral library M, this can be translated into the following objectives and constraints for an optimization problem:

• Minimize the reconstruction error w.r.t. M ∈ M, δ t,n , and s t,n , given by y t+1,n -M a t,n + δ t,n + s t,n ; • Preserve the nonnegativity and sum-to-one constraints on a t+1,n = a t,n + δ t,n + s t,n ; • Consider that δ t,n should be small and s t,n spatially sparse (i.e., nonzero only in a small number of pixels).

Rather than accounting for all objectives outlined above at once by devising a single computationally demanding optimization problem, we adopt an alternative strategy to obtain an efficient solution. First, let us assume that s t,n = 0. If δ t,n is sufficiently small and a t,n is a good estimate of the true abundance a t,n , we have:

a t,n + δ t,n ≈ a t,n ≈ a t,n . (6)
In that case, we can isolate the problem of estimating M ∈ M from that of estimating δ t,n and s t,n in order to solve it separately in a much simpler manner. We formulate the optimization problem as follows:

RE t+1,n = min M ∈M y t+1,n -M a t,n . (7)
The endmember matrix M t+1,n obtained by solving problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] can then be used to compute abundance vector a t+1,n with a single run of the FCLS algorithm. However, this strategy relies on a strong hypothesis, namely, s t,n = 0 and a t,n ≈ a t,n , in order for the approximation in [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF] to hold. This hypothesis may not be satisfied for all pixels. Fortunately, it turns out that we can devise a simple strategy to address those cases without significantly compromising the performance of the algorithm. Specifically, by evaluating the magnitude of the reconstruction error RE t+1,n in (7), we can indirectly identify if there were any significant changes in the abundance vector by testing whether RE t+1,n is larger than a given threshold RE 0 and, if so, estimate the corresponding abundance vector from scratch using MESMA or, alternatively, an algorithm such as AAM.

Intuitively, the reason this works is that, if the spectral library M is not too large, and if a t+1,n ≈ a t,n is not satisfied, then we cannot accurately reconstruct y t+1,n from problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF]. In the next section, we shall formalize this intuition by providing a deeper theoretical analysis of this method. Setting the threshold RE 0 offers a trade-off between accuracy (which tends to MESMA's for small RE 0 ) and computational performance. We propose to choose it as:

RE 0 = K U y∈U min M ∈M min a≥0, 1 a=1 y -M a , (8)
where U = {y 1 , . . . , y U } is a set with U pixels and K ∈ R + .

The pixels in U should be similar to those contained in the image sequence {y t,n } to be unmixed, so that RE 0 approaches K times the average optimal reconstruction error of the data. The proportion K controls how much the RE t+1,n of ( 7) can deviate from the estimated optimal value before we decide an abrupt change occurred. The inner optimization problem in ( 8) can be solved using MESMA or AAM. The complete procedure is detailed in Algorithm 1.

Algorithm IV. THEORETICAL GUARANTEES FM-MESMA (Algorithm 1) relies on important assumptions (such as [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF]) in order to split SU into EM selection and abundance estimation, and to correctly detect abundance changes while maintaining a small complexity. This raises questions regarding how the accuracy of the method is affected by the underlying assumptions and by the different variables involved in the model, such as: the noise e t,n , the small and abrupt abundance changes δ t,n and s t,n , the library M, and the accuracy of the estimated abundances at the previous time instants. To investigate these questions, in this section we derive theoretical results in the form of two theorems, which provide conditions under which 1) the EM matrix can be correctly recovered from [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF], and 2) the abrupt abundance changes s t,n can be correctly identified based on the reconstruction error RE t+1,n . Informally, the main findings related to each theorem can be summarized as: 1) When s t,n = 0, FM-MESMA recovers the correct EMs from [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] if the pairwise difference between the signatures in each library M p is sufficiently large when compared to a measure related to the coherence between signatures at different libraries, to the noise e t,n and to the abundance temporal variations δ t,n (or to the errors in a t,n ). 2) When an abrupt abundance change occurs, the value of RE t+1,n will be significantly larger than when s t,n = 0 (thus making it easy to detect) as long as s t,n is sufficiently large compared to δ t,n , to e t,n , and to the maximum pairwise difference between the signatures of each EM in the library. To proceed further, we shall assume that the true EM matrix M t+1,n in ( 5) is an element of the library M1 .

Theorem 1. Let us assume that a t,n = a t,n , s t,n = 0 and that M t+1,n ∈ M. Also assume that e t,n < Ω e , that max M ∈M M δ t,n < Ω δ and that the libraries M p satisfy:

min p min m,m ∈Mp m =m m -m 2 > Ω M (9) max p,q max m,m ∈Mp m, m ∈Mq m -m, m -m ≤ µ (10) Then, if 2 √ P (Ω e + Ω δ ) < Ω M -(P -1)µ the solution of the optimization problem (7) is M = M t+1,n .
Proof. Under these assumptions, optimization problem ( 7) can be written equivalently as

min M ∈M y t+1,n -M a t,n = min M ∈M M t+1,n (a t,n + δ t,n ) + e t+1,n -M a t,n (11) = min M ∈M (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n .
In order for the solution M * to problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] to be the correct endmember model, we need the reconstruction error to be minimal only when M t+1,n -M * = 0. In that case, the reconstruction error is: RE * t+1,n = M t+1,n δ t,n + e t+1,n . Equivalently, this means that ∀ M ∈ M, M = M t+1,n , the following condition must be verified:

M t+1,n δ t,n + e t+1,n (12)
< (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n .
Using the Reverse Triangle Inequality with the r.h.s. of the above expression leads to:

(M t+1,n -M )a t,n -M t+1,n δ t,n + e t+1,n (13)
≤ (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n .
Now, we can lower bound the l.h.s. of the previous expression. First, note that by hypothesis we have:

M t+1,n δ t,n + e t+1,n ≤ M t+1,n δ t,n + e t+1,n
< Ω e + Ω δ .

Before proceeding to bound the first term in [START_REF] Borsoi | Improved hyperspectral unmixing with endmember variability parametrized using an interpolated scaling tensor[END_REF], let us first derive the following auxiliary result concerning the matrix M t+1,n -M . Since L > P , the Geršhgorin circle theorem [START_REF] Horn | Matrix Analysis[END_REF], [START_REF] Richard | Online prediction of time series data with kernels[END_REF] can be used to derive a lower bound on the smallest singular value of M t+1,n -M as follows:

σ 2 min (M t+1,n -M ) = λ min Ξ ≥ min p Ξ pp - j:j =p |Ξ pj | , (14)
for Ξ = M t+1,n -M M t+1,n -M with Ξ ij being its (i, j)-th element, where we used the fact that Ξ is diagonalizable. Functions σ min (•) and λ min (•) denote the smallest singular value and the smallest eigenvalue of a matrix, respectively. Using condition (10) (i.e., the differences between signatures belonging to distinct EM libraries have low coherence), the second term in [START_REF]A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF] satisfies

j:j =p |Ξ pj | = j:j =p m p t+1,n -(m ) p , m j t+1,n -(m ) j ≤ (P -1)µ , (15)
where m p t+1,n and (m ) p denote the p-th column of M t+1,n and M , respectively. This leads to the following lower bound of ( 14):

min p Ξ pp - j:j =p |Ξ pj | ≥ min p m p t+1,n -(m ) p 2 -(P -1)µ > Ω M -(P -1)µ , (16)
since min p m p t+1,n -(m ) p 2 > Ω M due to condition [START_REF] Zhou | A Gaussian mixture model representation of endmember variability in hyperspectral unmixing[END_REF]. Since Ω M -(P -1)µ > 0, by combining the above equation with [START_REF]A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF] and taking the square root we can bound σ min (M t+1,n -M ) as:

Ω M -(P -1)µ < σ min (M t+1,n -M ) . (17)
Due to a t,n being confined to the unit simplex, its L 2 norm satisfies (1/ √ P ) ≤ a t,n . When combined with [START_REF] Imbiriba | Low-rank tensor modeling for hyperspectral unmixing accounting for spectral variability[END_REF], this leads to:

(1/ √ P ) Ω M -(P -1)µ < (1/ √ P )σ min (M t+1,n -M ) ≤ a t,n σ min (M t+1,n -M ) ≤ (M t+1,n -M )a t,n . (18)
This last equation constitutes the desired bound on the first term of [START_REF] Borsoi | Improved hyperspectral unmixing with endmember variability parametrized using an interpolated scaling tensor[END_REF]. Note that by hypothesis 2(Ω e + Ω δ ) < (1/ √ P ) Ω M -(P -1)µ, which implies that:

2 M t+1,n δ t,n + e t+1,n < (M t+1,n -M )a t,n .
Finally, subtracting M t+1,n δ t,n + e t+1,n from both sides of the previous expression leads to:

M t+1,n δ t,n + e t+1,n < (M t+1,n -M )a t,n -M t+1,n δ t,n + e t+1,n ≤ (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n ,
for any M ∈ M, M = M t+1,n , which is precisely the necessary condition stated in equation [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF]. Therefore, any M = M t+1,n leads to a larger reconstruction error than using M = M t+1,n , which is the optimal solution to [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF].

The rationale behind Theorem 1 is that, when the changes between the abundances a t,n and a t+1,n are small, then the variations observed between the pixel spectra y t,n and y t+1,n are due to 1) changes between the EM signatures M t,n and M t+1,n , and 2) to the presence of additive noise. Thus, if the signal to noise ratio is high, which is common in HIs, and the EM matrices contained in the library M are sufficiently different from one another, then the EM matrix in M that best reconstructs the new image y t+1,n (with a t+1,n ≈ a t,n ) as the optimal solution to problem (7) will be the correct EM matrix M t+1,n , since the increase in reconstruction error due to the choice of any other matrix in M will be greater than that due to the noise and abundance variations.

Although Theorem 1 clarifies the circumstances under which the optimization problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] provides the correct EM matrix, another important question is whether the detection strategy in FM-MESMA can correctly identify the cases when it fails, so that the solution must then be processed using MESMA. Recall this mechanism allows the algorithm to deal with abrupt abundance variations. It turns out that under similar assumptions to those made in Theorem 1, if the EM matrices contained in M are not too different from one another, then abrupt abundance variations can be correctly identified by means of the behavior of the reconstruction error of problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF], as stated in the following theorem.

Theorem 2. Assume that a t,n = a t,n , that s t,n = 0 and satisfies min M ∈M M s t,n > Ω s and that M t+1,n ∈ M.

Assume that e t,n < Ω e , that max M ∈M M δ t,n < Ω δ and that the maximum distances between signatures in each EM library is bounded as

max p max m,m ∈Mp m -m < Ω M . ( 19
)
Then, if √ P Ω M + Ω δ + Ω e < Ω s /(F + 1) the reconstruction error will be at least F times larger than if s t,n = 0.

Proof. Under these assumptions, the optimization problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] can be written equivalently as:

min M ∈M y t+1,n -M a t,n = min M ∈M M t+1,n (a t,n + δ t,n + s t,n ) + e t+1,n -M a t,n
Using the Reverse Triangle Inequality, we have: [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF] for any solution M ∈ M.

M t+1,n s t,n -M t+1,n (a t,n + δ t,n ) + e t+1,n -M a t,n ≤ M t+1,n (a t,n + δ t,n + s t,n ) + e t+1,n -M a t,n
We can upper bound the second term in the l.h.s. of the above expression, that is, the reconstruction error without the effect of s t,n , as follows:

M t+1,n (a t,n + δ t,n ) + e t+1,n -M a t,n (21)
= (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n ≤ (M t+1,n -M )a t+1 + M t+1,n δ t,n + e t+1,n .
Since a t+1 is confined to the unit simplex, its L 2 norm satisfies a t+1 ≤ 1, which allows us to bound the first term above using the properties of matrix norms as:

(M t+1,n -M )a t+1 ≤ M t+1,n -M a t+1 ≤ M t+1,n -M F ≤ max p √ P m p t+1,n -m p ≤ max p max m,m ∈Mp √ P m -m < √ P Ω M (22)
where m p t+1,n and m p are the p-th columns of M t+1,n and M , respectively. Thus, (21) can be bounded as:

M t+1,n (a t,n + δ t,n ) + e t+1,n -M a t,n < √ P Ω M + Ω δ + Ω e < Ω s (F + 1) -1 < M t+1,n s t,n (F + 1) -1 (23)
Now, by multiplying both sides of ( 23) by F + 1, subtracting (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n from each side, and finally using the result in [START_REF] Roberts | Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[END_REF], we obtain:

F (M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n < M t+1,n s t,n -(M t+1,n -M )a t,n + M t+1,n δ t,n + e t+1,n ≤ RE * t+1,n ≤ M t+1,n (a t,n + δ t,n + s t,n ) + e t+1,n -M a t,n
for any solution M ∈ M, where RE * t+1,n is the optimal reconstruction error of problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF]. This result means that, if s t,n = 0, the reconstruction error RE * t+1,n will be at least F times larger than it would be if s t,n = 0.

The rationale behind Theorem 2 is that, when there are large abundance changes between two consecutive time instants, that is, s t,n is large, then as long as the maximum pairwise difference between the elements in M is not too large, it will not be possible to reconstruct the pixel y t+1,n with good accuracy with any EM matrix in M because a t,n is too distant from a t+1,n . This will lead to a large reconstruction error, what can be explored in order to identify such changes.

An important aspect of FM-MESMA is the assumption that an accurate estimate of the abundances is available at time t (i.e., a t,n ≈ a t,n ). Due to the sequential nature of FM-MESMA, this approximation affects the results at time t + 1. Since problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] actually relies on the approximation a t,n ≈ a t+1,n to perform EM selection, errors in a t,n will affect it in the same way as δ t,n (or, in the case of large errors, as s t,n ). Thus, while errors in a t,n can impact the performance of the method negatively, their effect can be controlled: if the errors are large and K is properly selected, they will have the same effect as sudden changes, and cause the algorithm to reprocess the given pixel with MESMA or AAM.

A. Computational complexity analysis

In this section we show how the computational complexity of the proposed algorithm compares to that of MESMA and AAM. Before proceeding, let us denote by κ ∈ [0, 1] the average proportion of pixels that undergo changes between consecutive time instants t and t + 1.

The operations in FM-MESMA at each time instant t with non-negligible computational complexity consist of 1) solving of optimization problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] for all N image pixels, and 2) running MESMA or AAM to estimate the abundances of pixels that did undergo significant changes. For the first step, a simple enumeration strategy can be employed to solve problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] by testing all possible EM models, which results in a complexity of about O( P p=1 C p LP ). For the second step, the complexity will be on average κ times the one of MESMA or of AAM, which will be discussed below.

Determining the computational complexity of MESMA is less direct, since the optimization problem (4) does not have a closed form solution and depends on iterative algorithms. For simplicity, we assume that the inner FCLS problem in MESMA ( 4) is solved by approximately translating it into a nonnegative least squares (NNLS) problem of dimensions (L + 1) × (P + 1) [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]. Different methods have been proposed to solve NNLS problems, including the active set methods, the interior point method, and other iterative approaches [START_REF] Chen | Nonnegativity constraints in numerical analysis[END_REF]. Although iterative approaches perform better in large scale problems, the interior point method works well for problems like MESMA and can give us an idea of its complexity [START_REF] Polyak | Projected gradient method for non-negative least square[END_REF].

The interior point method needs O(ln( -1 )) iterations, each with complexity (P + 1) 3 , in order to achieve a reconstruction error that is -close to the global optimum [55, p. 393], i.e.,

y n -M a n 2 -y n -M a * n 2 ≤ (24)
where a * n is abundance vector that minimizes the NNLS problem for a given M . Thus, since this problem is solved for each EM model, the computational complexity of MESMA is approximately O( P p=1 C p P 3 ln( -1 )). It can be seen that although the complexity of the proposed method and MESMA scale similarly as the size of the libraries C p increases, as long as κ is small the proposed algorithm can handle scenarios with moderately large P much more easily since problem (7) scales linearly with P . This characteristic needs to be emphasized because, while several methods can effectively address the problem of reducing the size C p of the libraries by removing redundant signatures [START_REF] Roth | Comparing endmember selection techniques for accurate mapping of plant species and land cover using imaging spectrometer data[END_REF]- [START_REF] Dennison | A comparison of error metrics and constraints for multiple endmember spectral mixture analysis and spectral angle mapper[END_REF], a larger number P of EM classes cannot be so easily circumvented. The AAM algorithm [START_REF] Heylen | Hyperspectral unmixing with endmember variability via alternating angle minimization[END_REF], for instance, requiring O(P 2 P (L 3 + P L max p C p )) operations at every iteration, has a larger base cost depending on terms such as L 3 and scales quickly with P , making it costly in scenarios where P is large and C p is small. Such cases in which the proposed algorithm is particularly faster (i.e., small C p ) are of special practical interest. This is because recently proposed state of the art approaches to spectral library reduction for MESMA have reported experimental results indicating that libraries could be reduced to between two and five (averaging three) signatures per EM without an appreciable drop in performance [START_REF] Meerdink | Developing spectral libraries using multiple target multiple instance adaptive cosine/coherence estimator[END_REF]. Finally, we note that the best choice, in terms of complexity, among using MESMA or AAM to unmix the changed pixels in the proposed algorithm may depend on each scenario. If P or C p is small, MESMA can be a good choice due to its smaller base cost compared to AAM, whereas for moderate C p AAM will perform faster.

V. EXPERIMENTAL RESULTS

We shall now evaluate the performance of the proposed FM-MESMA algorithm using simulations with synthetic, semireal, and real data. Our method is compared with MESMA and AAM [START_REF] Heylen | Hyperspectral unmixing with endmember variability via alternating angle minimization[END_REF], which are both library-based methods, with fully constrained least squares (FCLS) algorithm, and with the online unmixing algorithm (OU) [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]. The OU algorithm estimates both the abundances and one set of EMs for each time instant blindly from the HI using a two-stage stochastic optimization procedure. The EMs are modelled in OU as temporally smooth additive perturbation over a mean EM matrix. FM-MESMA is implemented using MESMA to unmix the significantly changed pixels in step 10. In all simulations, the endmembers for the FCLS were extracted from the HI using the VCA algorithm [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF]. For FM-MESMA, we computed RE 0 in (8) using the pixels in the HI at the initial time instant t = 1, i.e., U = {y 1,1 , . . . , y 1,N }. To evaluate the performance of the algorithms, different metrics were considered depending on the simulation setups. Performance metrics that are specific to the simulations with synthetic data will be defined in Section V-A. For the simulations with semi-real and real data, we considered as metrics the root mean squared error (RMSE) and the spectral angle mapper (SAM). The RMSE between two sequences of matrices X t and X * t , for t = 1, . . . , T , is defined as:

RMSE X = T t=1 1 T N X X t -X * t 2 F
where N X is the number of elements in X t . The SAM between the true and estimated endmembers is defined as:

SAM M = 1 T N P T t=1 N n=1 P p=1 arccos (m p t,n ) m p t,n m p t,n m p t,n
, where m p t,n and m p t,n are the p-th columns of the true and the estimated endmembers, respectively.

A. Synthetic Data

The simulations with synthetic data were designed to illustrate how FM-MESMA performs when compared to MESMA and AAM in four different ways, namely:

1) Its computational cost for different values of P and C p (Section V-A1); 2) Its accuracy when detecting abrupt abundance changes, and the effect of the proportion of changed pixels (i.e., κ) on the computational cost of FM-MESMA (Section V-A2); 3) Its accuracy when recovering the EMs from M for different amounts of abundance temporal variation δ t,n and different signal to noise ratios (SNRs) (Section V-A3); 4) Its accuracy when recovering the EMs from M for different amounts of library variance σ 2 M and different SNRs (Section V-A4).

These different experiments provide an empirical assessment of some of the theoretical results derived in Section IV. For the general simulation setup, we considered sequences of T = 11 images with L = 200 bands and N = 1000 pixels. Each pixel y t,n was generated according to the LMM in (5a), where the true endmember matrix M t,n was sampled uniformly from a library M, and e t,n was a white Gaussian noise. The EM library M was generated randomly. The mean µ p , p = 1, . . . , P of each material was first sampled from a uniform distribution over interval [0, 1] L . Then each EM signature in M p , p = 1, . . . , P was generated as a sample from an isotropic Gaussian distribution N (µ p , σ 2 M I) and truncated in the interval [0, 1] L . The abundances a 1,n were sampled from a Dirichlet distribution. Between each pair of images at instants t and t + 1, a proportion κ of the N pixels was changed with new samples, whereas the remaining ones were kept constant unless otherwise specified. The other parameters such as K, the SNR, the change ratio κ and the library variance σ 2 M will be specified in the following for each experiment. In general, the library variance is defined as:

σ 2 M = 1 LP P p=1 tr cov(m p , m p ) , (25)
where tr{•} is the matrix trace operator and cov(m p , m p ) denotes the covariance matrix of the signatures from the pth EM, which can be estimated using the samples in the library M p if unknown. 1) Computational complexity analysis: We evaluated the execution time of MESMA, AAM and FM-MESMA for different values of P ∈ {2, . . . , 9} and C p ∈ {2, . . . , 10}, p = 1, . . . , P . We considered an SNR of 40dB, a change ratio of κ = 0.01, a library variance of σ 2 M = 0.12 and K = 10 in Algorithm 1. The results are shown in Table I. It can be seen that FM-MESMA has a significantly smaller execution time when compared to MESMA. Moreover, there are also significant improvements over AAM when P is large and C p is small or moderate. Note that although AAM performs better when C p is large, this situation is avoided in practice without significant impact in the performance by removing some redundant signatures from the library [START_REF] Roth | Comparing endmember selection techniques for accurate mapping of plant species and land cover using imaging spectrometer data[END_REF]- [START_REF] Dennison | A comparison of error metrics and constraints for multiple endmember spectral mixture analysis and spectral angle mapper[END_REF], [START_REF] Roberts | Evaluation of the potential of hyperion for fire danger assessment by comparison to the airborne visible/infrared imaging spectrometer[END_REF]. However, this is not the case for larger P , in which case FM-MESMA leads to a performance improvement. We solved (7) using an exhaustive search procedure. More efficient solutions will be devised in the future.

2) Large change detection analysis: We evaluated the accuracy of FM-MESMA when detecting pixels containing large abundance changes to reprocess them using the MESMA or AAM algorithms. This experiment also allowed us to check numerically the theoretical results in Theorem 2. The experiment is divided in two parts. First, we evaluate the effect of K and κ on the computational burden of the proposed method by measuring the proportion of pixels marked as changes depending on these variables. Afterwards, we compare the change detection accuracy of FM-MESMA to that of the other algorithms (where change detection was performed using the strategy in [START_REF] Ertürk | Informative change detection by unmixing for hyperspectral images[END_REF]) by fixing κ = 0.2 and evaluating the probability of detection (PD) against the probability of false alarm (PFA), where PD and PFA are defined as:

PD A = T t=2 N n=1 χ 0 (a t,n -a t-1,n ) • ι s,t,n (T -1) N n=1 χ 0 (a t,n -a t-1,n ) PFA A = T t=2 N n=1 max ι s,t,n -χ 0 (a t,n -a t-1,n ), 0 (T -1) N - N n=1 χ 0 (a t,n -a t-1,n )
where χ 0 (•) is the indicator function of the set 0 = {0} (i.e., χ 0 (x) = 1 if x = 0 and 0 otherwise). For both cases, we considered an SNR of 30dB, a library variance of σ 2 M = 0.12, P = 4 EMs and C p = 3. The results are shown in Figure 1. It can be seen that when K was not too close to one, the number of times a pixel had to be reprocessed depended mostly on the actual amount of changed pixels in the scene κ. In terms of accuracy, the change detection performance was satisfactory for all methods, with MESMA showing the best performance followed closely by FM-MESMA, and OU and FCLS being slightly worse. This suggests that one can select a moderate value of K and obtain a good change detection accuracy without increasing the complexity unnecessarily. However, devising a strategy to select an optimal value for K is more complex and will be left as a subject for future work.

3) Library variance analysis: We evaluated how often the solution of problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] was the same as the true EM matrix M t,n when compared to that of MESMA, for different library variances. This experiment also allowed us to validate numerically one part of Theorem 1. We considered κ = 0.05, P = 4, C p = 3, and different values of K ∈ {1, . . . , 20} and σ 2 M = {0.02, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 0.7, 1, 1.5}. The performance was evaluated using the endmember positive predictive value (PPV)

PPV M = 1 T T t=1 N n=1 χ 0 (M t,n -M t,n ) N . ( 26
)
where M t,n ∈ M and M t,n ∈ M are the true and estimated EM matrices, respectively. The results are shown in Figure 2 for SNRs of 25, 35 and 45db, where the red dashed line depicts the MESMA results. The results show that the values of PPV M for FM-MESMA were very similar to those obtained by MESMA for the different values of σ 2 M and K. Moreover, the PPVs was consistently better for larger σ 2 M , what agrees with the conclusions of Theorem 1. While a smaller SNR negatively affected the overall accuracy of both algorithms, their relative behavior was not affected. Different values of K also did not have a significant effect on the results obtained by FM-MESMA for the selected proportion of large changes. This suggests that FM-MESMA is a computationally efficient and accurate alternative to MESMA in those circumstances.

4) Abundance variance analysis: We evaluated how often the solution to problem [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF] was the true EM matrix M t,n when compared to that of MESMA, for different amounts of small abundance variances, which are represented by δ n in (5b). This experiment also allowed us to validate numerically one part of Theorem 1. We considered κ = 0, P = 4, C p = 3, and different values of K ∈ {1, . . . , 20}. The variable abundances were generated by sampling each a t,n from Dirichlet distributions with mean µ a n , n = 1, . . . , N and standard deviations ranging from 10 -6 to 1. The abundance mean vectors µ a n were also sampled from a Dirichlet distribution. The performance were evaluated again using the endmember PPV as defined in [START_REF] Borsoi | Super-resolution for hyperspectral and multispectral image fusion accounting for seasonal spectral variability[END_REF].

The results are shown in Figure 3 for SNRs of 25, 35 and 45db, with the red dashed line depicting the MESMA results. It can be seen that for smaller abundance temporal variations (e.g., about ≤ 10 -2 ), FM-MESMA was able to obtain an endmember positive predictive value that was high and very similar to that of MESMA, even for larger values of K. This again supports the use of FM-MESMA as an alternative to MESMA under these circumstances. For larger temporal variations, a decrease in performance was observed for both MESMA and FM-MESMA. This is due to the fact that, when the variance of a t,n is large, its density tends to concentrate at the edges of the simplex, what causes many abundance fractions to be close to zero, making the identification of the correct M t,n through problems ( 4) or ( 7) more difficult. Although a smaller SNR negatively impacts the overall accuracy of the algorithms (like in the simulations of Section V-A3), the relative behavior between the algorithms remains approximately the same.

B. Abundance estimation performance in semi-real data

We evaluated FM-MESMA in terms of abundance and endmember estimation accuracy by comparing it to the FCLS, OU, MESMA, and AAM algorithms using a synthetic data set designed to closely emulate a practical scenario. To this end, first we manually extracted P = 3 sets of spectrally distinct pure pixels from different materials (tree, road and water) from the Jasper Ridge HI containing six signatures with L = 198 bands each. Then, we randomly sampled signatures from these sets in order to create one pair of disjoint spectral libraries M 1 p , M 2 p , i.e., M 1 p ∩ M 2 p = ∅ for each material p = 1, . . . , P , where both M 1 p and M 2 p have three signatures each. Afterwards, two sets of temporal abundances were generated. The first one was generated randomly to allow for a statistical evaluation, in which N = 1000 pixels in a sequence of T = 20 images were generated following the model in (5a), where the abundances a t,n were sampled from a Dirichlet distribution and a proportion of κ = 0.05 of the pixels undergoing significant changes. The second set of abundance maps consisted of a single sequence of T = 6 images containing N = 2500 pixels whose spatial compositions were adequate for a visual inspection of the results, and large changes were added in the form of random convex polygons. The latter sequence can be seen at the top row of Figure 4. For both cases, M t,n was sampled uniformly from the library M 1 p and the additive noise e t was selected as white Gaussian with an SNR of 30dB.

The library-based methods, namely, MESMA, AAM and FM-MESMA, were then used to unmix the images using the library M 2 p , with K = 10 selected for the proposed method. This generated a mismatch between the signatures in the HI and those used for SU, what is commonly observed in practice. To reliably measure the performance of the methods, we ran this simulation over 100 Monte Carlo realizations. The average performance of all algorithms are depicted in Table II. It can be observed that the library-based methods provided a considerable improvement in abundance and endmember estimation accuracy when compared to both FCLS and OU.

The results of AAM and MESMA were also very similar, with those of AAM being slightly worse. FM-MESMA provided an improvement of about 16% in abundance estimation when compared to MESMA. This confirms the benefits of exploring the temporal correlation between the abundances at adjacent time instants. In terms of endmember estimation accuracy, FM-MESMA performed similarly to MESMA, with a slightly better RMSE M but slightly worse SAM M . In terms of RMSE Y , the library-based methods achieved a similar reconstruction error (with MESMA's being slightly smaller as it solves the EM selection problem exactly), significantly smaller than FCLS but still larger than OU. Since the OU algorithm estimates the EMs from the scene, it has more degrees of freedom and is able to achieve a smaller RMSE Y . However, it is well-known that smaller reconstruction errors do not necessarily translate into better abundance map estimates.

For the second set of abundance maps, the quantitative and visual results are shown in Table III and in Figure 4. Due to space limitations 1) we do not show the AAM results since they were very similar to those of MESMA, and 2) only a composite image is shown, with red, green and blue corresponding to the abundances of the first, second, and third EMs, respectively. The quantitative metrics show that the algorithms behaved very similarly to the previous case, with the minor differences that MESMA and AAM now performed slightly better in terms of RMSE A and RMSE M , respectively, when compared to the remaining methods. Visually, it can be seen that the results of MESMA and of FM-MESMA are similar and approach the ground truth more closely when compared to FCLS and OU, whose results appear less accurate and more skewed towards blue for the FCLS, and towards red for OU when compared to the ground truth. The sudden changes in the ground truth sequence were also well captured by FM-MESMA. However, various pixels in the lower-left segment were also incorrectly marked as changes, indicating that the selection of the proportion K was slightly conservative.

C. Real Data

For the simulations with data, we considered the Lake Tahoe data set, which was originally presented in [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF]. This where three predominant materials (soil, water and vegetation) can be identified. The images were ordered according to the day/month in which they were acquired in order to make the seasonal changes more clear. The EM libraries for MESMA and for FM-MESMA were constructed as follows. First, reference EMs for each material were extracted by applying the VCA to the concatenation of the pixels of all six images. Then, a preliminary library was constructed by extracting pixels from each image that had small spectral angle to the reference endmembers. Finally, the libraries were obtained by removing/pruning the most redundant signatures (as measured according to their mutual Euclidean distance) so that the final libraries contained C p = 6 signatures for each EM. We set K = 10 for FM-MESMA, and the parameters for the OU algorithm were the same as those used in [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF].

The abundance maps estimated by the algorithms are shown in Figures 6, 7 and 8. Due to lack of space, we do not show the AAM results since they were very similar of MESMA. It can be seen that the FCLS performed poorly for some of the frames (e.g., the second and the fourth ones), with significant confusion between different materials. Although OU showed more consistent results, they were not as good as those by MESMA or FM-MESMA. Specifically, the OU abundances did not show a separation between vegetation and soil as clear as observed in Figure 5. Moreover, significant water abundances were predicted outside of the lake. The results of MESMA and FM-MESMA were very similar, and closely agreed with the distribution of the endmembers observed in Figure 5. However, some slight differences between the methods can be noticed, notably the abundances by the proposed algorithm being smoother in time as in the case of the water endmember. The changes detected by FM-MESMA, shown in Figure 9, clearly distinguish the seasonal variations at the borders of the lake and in the crop circles due to different water levels and in the growth states, respectively.

The reconstruction errors RMSE Y and execution times of the algorithms are provided in Table IV. The RMSE Y results behaved very similarly to the semi-real case discussed in Section V-B, with OU achieving the smallest and FCLS the largest reconstruction errors, and the library-based methods performing similarly to each other. The execution time of FM-MESMA was about half that of MESMA even though there were significant changes between some pairs of frames, what agrees with the results in Table I. Moreover, this scene contains a small number of materials and is thus not a good representative of the relative performance between the algorithms in other scenarios. For image sequences with larger values of P , the computational complexity gains will be much more significant.

VI. CONCLUSIONS

In this paper, we presented a new computationally efficient multitemporal unmixing algorithm (FM-MESMA) based on multiple endmembers spectral mixture analysis. The proposed strategy exploits the high temporal correlation of the abundance maps in order to improve both the accuracy and the computational complexity of the algorithm. Specifically, it approximates the solution to the multitemporal unmixing problem by separating it into two sub-problems, namely, endmember selection and abundance estimation, which are much easier to solve individually. A strategy was also proposed to detect abrupt abundance changes by analysing residuals of the endmember selection problem. Theoretical results demonstrated how FM-MESMA compares to MESMA in terms of quality and effectiveness in detecting abrupt abundance changes. Besides, these results also provide valuable insight into the conditions under which the approximate algorithm succeeds. Simulation results showed that the proposed method gives results with quality similar to, or better than, both MESMA and parametric models at a reduced computational complexity.

Figure 1 .

 1 Figure 1. Large changes analysis: Number of detected changes as a function of K and κ (top), and change detection performance for κ = 0.2 (bottom) (AAM results were omitted due to being similar to MESMA's).

Figure 2 .

 2 Figure 2. Library variance analysis: Correct EM recovery rate as a function of σ 2 M and K, for SNRs of 25dB (left), 35dB (middle) and 45dB (right).

Figure 3 .

 3 Figure 3. Abundance variance analysis: Correct EM recovery rate as a function of abundances standard deviation and K, for SNRs of 25dB (left), 35dB (middle) and 45dB (right).

Figure 4 .

 4 Figure 4. Composite color image of abundance maps (i.e., abundances for EMs 1, 2 and 3 correspond to red, green and blue colors) in the visual semireal dataset (top) and changes detected by the proposed method bottom).

Figure 5 .

 5 Figure 5. Lake Tahoe HI sequence (heading text means day/month).

Figure 6 .

 6 Figure 6. Multitemporal abundances of the Lake Tahoe HI for the water EM.

Figure 7 .

 7 Figure 7. Multitemporal abundances of the Lake Tahoe HI for the soil EM.

Figure 8 .

 8 Figure 8. Multitemporal abundances of the Lake Tahoe HI for the vegetation EM.

Figure 9 .

 9 Figure 9. Detected changes for the Lake Tahoe HIs.

  1: FM-MESMA Algorithm Input : Multitemporal HS images {y t,n }, endmember library M, threshold proportion parameter K.

	11	end
	12	end
	13 end

[START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF]

Perform SU for the first HI y 1,n , using the MESMA or AAM algorithms to obtain a1,n and M 1,n, for n = 1, . . . , N ; 2 Compute RE0 according to (8) ; 3 Initialize the change maps indicator function as ιs,t,n = 0, t = 1, . . . , T , n = 1, . . . , N ; 4 for t = 1, . . . , T -1 do 5 for n = 1, . . . , N do 6 Solve problem (7) to obtain the EM matrix M t+1,n and the reconstruction error REt+1,n ; 7 if REt+1,n ≤ RE0 then 8 Estimate at+1,n using the FCLS method with M t+1,n as the EM matrix; 9 else 10 Estimate at+1,n and M t+1,n using MESMA or AAM and set ιs,t,n = 1 ; 14 returnEstimated abundances at,n , detected change maps ιs,t,n , and endmembers M t,n ;

Table I EXECUTION

 I TIMES FOR DIFFERENT VALUES OF P AND Cp (BEST RESULTS MARKED IN BOLD, BEST RESULTS BY OUR METHOD MARKED IN RED).

		Method	Cp	2	3	4	5	6	7	8	9	10
	P = 2	MESMA		1.20	1.21	1.27	1.43	1.37	1.57	1.60	1.66	1.96
		AAM		6.61	6.48	6.67	7.11	6.79	6.88	6.98	7.06	7.17
		FM-MESMA	1.19	1.16	1.18	1.30	1.21	1.31	1.31	1.35	1.55
	P = 3	MESMA		1.55	1.78	2.32	2.88	4.52	6.26	6.79	9.46	15.90
		AAM		30.73	30.41	31.62	32.12	32.04	32.62	33.15	33.82	33.48
		FM-MESMA	1.44	1.52	1.60	1.75	2.08	2.45	3.05	4.03	4.93
	P = 4	MESMA		2.91	4.97	11.04	36.74	52.37	103.38	180.75	282.21	490.92
		AAM		98.53	98.12	101.76	103.84	103.99	106.57	107.61	108.82	109.05
		FM-MESMA	2.14	2.30	3.07	5.17	6.88	11.58	16.79	28.87	45.57
	P = 5	MESMA		5.04	24.34	74.94	272.21	680.85	1482.15	2780.66	4882.36	7783.21
		AAM		272.22	270.72	282.31	285.31	289.26	292.36	296.11	295.97	300.89
		FM-MESMA	2.68	4.21	8.80	25.38	50.61	110.60	234.83	383.27	645.80
	P = 6	MESMA		8.61	81.69	398.69	1600.86	4625.57	12153.74 24423.06	∞	∞
		AAM		678.18	691.35	715.88	728.04	745.54	751.10	761.53	760.70	760.47
		FM-MESMA	3.66	9.30	36.11	135.68	405.91	1069.17	2614.65	4707.32	8513.08
	P = 7	MESMA		19.58	257.71	1933.77	8554.18	31760.77	∞	∞	∞	∞
		AAM		1653.51 1671.48	1726.08	1753.13	1788.30	1808.63	1829.27	1843.13	1861.54
		FM-MESMA	5.07	23.60	158.68	736.29	2660.14	7728.61	21816.81	∞	∞
	P = 8	MESMA		36.29	748.21	8394.23	49221.23	∞	∞	∞	∞	∞
		AAM		3876.57 3932.96	4079.04	4112.40	4192.04	4368.23	4406.82	4453.85	4506.10
		FM-MESMA	6.50	68.31	656.97	3869.42	17330.46	∞	∞	∞	∞
	P = 9	MESMA		75.84	3167.33 38765.51	∞	∞	∞	∞	∞	∞
		AAM		8972.31 9126.39	9354.76	9505.15	9694.71	9861.76	9960.79	10041.00 10178.05
		FM-MESMA	9.78	211.21	2913.83	19922.95	∞	∞	∞	∞	∞

Table II AVERAGE

 II ABUNDANCE AND EM ESTIMATION RESULTS FOR THE SEMI-REAL SIMULATIONS (VALUES ARE ×10 2 ).

		FCLS	OU	MESMA AAM Proposed
	RMSE A	3.51	2.27	1.87	1.90	1.57
	RMSE M	-	45.6	37.3	38.5	36.2
	SAM M	-	24.3	10.7	10.7	11.2
	RMSE Y	1.89	0.18	0.85	0.88	0.96

Table III QUANTITATIVE

 III RESULTS FOR THE DATA IN FIGURE 4 (VALUES ARE ×10 2 ).

		FCLS	OU	MESMA AAM Proposed
	RMSE A	2.32	2.16	1.11	1.27	1.11
	RMSE M	-	35.7	33.8	36.4	34.0
	SAM M	-	16.4	11.2	10.9	11.6
	RMSE Y	5.07	0.26	0.80	0.83	0.96

Table IV EXECUTION

 IV TIMES AND RMSE Y FOR THE LAKE TAHOE HI.

		FCLS	OU	MESMA AAM Proposed
	Time [s]	12.3	85.5	112.0	228.8	62.8
	RMSE Y	0.412	0.013	0.064	0.068	0.072

Without loss of generality, one can always assume that M t+1,n ∈ M by incorporating any error ∆M t+1,n into the additive noise term e t+1,n by adding ∆M t+1,n a t+1,n .

This work has been supported by the National Council for Scientific and Technological Development (CNPq) under grants 304250/2017-1, 409044/2018-0, 141271/2017-5 and 204991/2018-8. R.A. Borsoi is with the Department of Electrical Engineering, Federal University of Santa Catarina (DEE-UFSC), , Nice, France.


Maha Bouhadida

Jean-Charles Beugnotb

Philippe Delayea

Kien Phan Huyb

Sylvie Lebruna

Conception and reproducibility study of efficient evanescent Raman converters based a nanofiber immersed in a liquid

Keywords: nanofiber, evanescent nonlinearity, Raman scattering

We demonstrate, build and optimize evanescent Raman converters at the sub-nanosecond regime based on a silica nanofiber immersed in ethanol. Two different standard silica fibers (SMF28, 460HP) are tested and compared. The converters are pumped at 532 nm and deliver pulses at 630 nm, which is the first Stokes order wavelength of ethanol. They present highly reproductible performances. A maximum output Stokes energy of 0.29 J is usually reached with an external conversion efficiency of 60%. Lowering the Raman threshold and pushing up the nanofiber breakdown allow a higher conversion operating range, that is the conception key of these converters.

INTRODUCTION

Optical nanofibers are fibers stretched until reaching diameters comparable or smaller than the light wavelength. For these diameters, nanofibers exhibit a strong confinement of light enabling the investigation of nonlinear effects [1] like the generation of supercontinuum [2]. They also have an intense evanescent field that can be exploited for different applications such as optical traps [3], optical sensing [4] and spectroscopy [5]. However, the generation of nonlinear effects in the evanescent field of the nanofiber remains limited. In this paper, we are interested in the Raman scattering in the evanescent field of a nanofiber immersed in a liquid as the "evanescent nonlinearity". The evanescent Raman effect was previously observed with a few-cm-length nanofiber immersed in ethanol and in toluene by using sub-nanosecond pump pulses at 532 nm. This first observation was investigated using nanofibers that allow a maximal modal Raman gain [6]. The output energy of the first Stokes order of the Raman liquid was limited to 0.11 J by the breakdown of the nanofiber and never exceeded the transmitted pump energy. Extracting more Stokes energy requires increasing the Raman conversion operating range. It is then mandatory to decrease the Raman threshold and to increase the nanofiber breakdown threshold. One alternative is to decrease the pump pulse laser duration until a few tens of ps. Another alternative is to keep the subnanosecond regime and to increase the diameter of the nanofiber. However, higher diameters would increase the Raman threshold because the intensity of the evanescent field would be decreased. A compromise has to be established by finding the nanofiber diameter that enables the extraction of the maximal output energy at the Stokes wavelength. In this context, we propose a design of a Raman converter based on a nanofiber immersed in a liquid that allows to extract the maximum of the Stokes energy [7]. Two different fibers are compared. The efficiency and the reproducibility of this evanescent Raman conversion are studied.

The paper is structured as follows: we dedicate the first section to model the evanescent Raman conversion. The structure of the converter is described in the second section. In the following section, we detail the nonlinear characterization of the proposed converter and we study the efficiency and the reproducibility of its conversion in the last section. We obtain an output Stokes energy of 0.29 J at 630 nm i.e. a value that is almost three times higher than the previous results. An external Raman conversion efficiency around 60% is reached with a high reproducibility opening the way to a new family of low-cost compact and efficient all fibered Raman sources that can be directly inserted in optical fibered networks.

EVANESCENT RAMAN CONVERSION MODELLING

The nanofiber is composed of a waist and two symmetrical tapers which are the transition sections between this waist and the unstreched parts of the fiber L1 and L2 as shown in Figure 1. The input pump beam is injected in L1 and propagates in the core of the fiber. It then enters the left taper where the fiber dimensions (core radius rco, cladding radius rclad) start to be progressively reduced. The taper is modeled as a three-layer structure combining the fiber core, the fiber cladding and the surrounding external medium. The pump beam enters the nanofiber waist which is assimilated to a two-layer structure composed by the fiber cladding and the external medium since the fiber core can be totally neglected in this part. The modes are considered to be scalar in the taper. For the nanofiber's radii equal or lower than the wavelength, the propagating modes are assumed to be vectorial because of the large index. The beam then propagates in the right taper and is collected at the end of the un-stretched part L2 as core mode. By immersing the nanofiber in a liquid, the incident pump photons of the evanescent field are scattered through Raman diffusion under the form of lower energy photons called Stokes photons.

That is what we call "evanescent Raman conversion" [7,8]. For example, if the pump pulses are emitted at p = 532 nm and the immersing liquid is pure ethanol, the first Stokes order photons are scattered at s = 630 nm. The evanescent part of the mode propagating in the external medium is significant. Figure 2 shows the fraction of the evanescent part of the fundamental mode field 𝑓 𝑒𝑣 as a function of the nanofiber radius for p = 532 nm and pure ethanol whose expression is the following:

𝑓 𝑒𝑣 = ∬ (𝑒 × ℎ). 𝑧ê 𝑣𝑎𝑛𝑒𝑠𝑐𝑒𝑛𝑡 𝑎𝑟𝑒𝑎 𝑑𝐴 ∬ (𝑒 × ℎ). 𝑧t 𝑜𝑡𝑎𝑙 𝑎𝑟𝑒𝑎

𝑑𝐴

Where 𝑒 and ℎ are the respectively electric and the magnetic fields. For the nanofiber's radii range we are considering (between 220 nm and 300 nm) and for the fundamental mode HE11, 𝑓 𝑒𝑣 is always higher than 0.38.

One key parameter behind the evanescent Raman conversion is the evanescent Raman threshold. We consider that the Raman threshold is achieved when the transmitted pump energy Ep,out equals the first Stokes order energy Es,out. We introduce the critical parameter  expressed as 𝛾 = 𝑔 𝑅𝑚 𝐿 𝑤 𝑃 𝑖𝑛 where 𝑔 𝑅𝑚 is the modal Raman gain (in m-1.W-1) and 𝑃 𝑖𝑛 is the peak pump power [9]. For sub-nanosecond pulses,   23 at the Raman threshold [6,7]. The expression of the modal Raman gain 𝑔 𝑅𝑚 is the following [10]:

𝑔 𝑅𝑚 = 𝜀 0 𝑐 2 ∬ 𝑔 𝐿 𝑛 𝐿 2 |𝑒 𝑝 . 𝑒 𝑠 | 𝑎𝑐𝑡𝑖𝑣𝑒 𝑎𝑟𝑒𝑎 2 𝑑𝐴 ∬ (𝑒 𝑝 × ℎ 𝑝 ). 𝑧t 𝑜𝑡𝑎𝑙 𝑎𝑟𝑒𝑎 𝑑𝐴 ∬ (𝑒 𝑠 × ℎ 𝑠 ). 𝑧t 𝑜𝑡𝑎𝑙 𝑎𝑟𝑒𝑎 𝑑𝐴 (2)
Figure 1. Evanescent Raman conversion scheme. Figure 2. Evanescent fraction of the pump field on HE11 at 532 nm versus the nanofiber radius.

(

) 1
where 𝑔 𝐿 (in m. W-1) is the Raman gain coefficient of the liquid and n 𝐿 is its refractive index. 𝜖 0 and 𝑐 are, respectively, the dielectric permittivity and the light celerity in vacuum; 𝑒 and h represent respectively the electric and the magnetic fields; p and s stand for, respectively, pump and Stokes. Since the modal Raman gain of the crossed interactions between modes is very weak, we assume that both the pump and the Stokes beams are propagating on the same mode.

The variations of 𝑔 𝑅𝑚 with the nanofiber radius for the modes HE11, TM01, TE01 and HE21 are shown in Figure 3. We observe that there is a radius optimizing 𝑔 𝑅𝑚 for each plot. When the radius is smaller, the beam structure approaches a plane wave spreading along important distances and so, the amplitude of the evanescent field decreases reducing the gain. However, larger radii make the mode more confined in the fiber core and so, the Raman gain decreases.

For nanofiber radii between 220 nm and 350 nm and an input peak power of 1kW, the evanescent Raman threshold can be achieved with waist lengths Lw between 5 cm and 8 cm which is possible with our pulling system. For the higher order modes, higher radii are needed to achieve the maximal modal Raman gain. This gain does not exceed 0.2 m-1. W-1 but it can be compensated by pulling longer nanofibers (Lw = 15 cm for TM01 and at R = 508 nm).

EVANESCENT RAMAN CONVERTER STRUCTURE

The proposed evanescent Raman converter is based on the component drawn in the section above and depicted in Figure 1.

The immersed parts in the Raman liquid are the nanofiber waist and the tapers. L 1 and L 2 are needed to make the manipulation easier and to facilitate the light injection and collection. L 1 and L 2 should be long enough for practical use but short enough to avoid parasitic nonlinear effects as it will be explained later.

The nanofiber is fabricated following the "pull and brush" technique with the pulling system described in [6]. A butane flame softens the central part of the fiber while two translation stages elongate it. The flame flow and intensity and the stages are fully computer-controlled, which allows to create reproducible nanofibers with high light transmission. Based on the target parameters of the nanofiber (nanofiber radius R, waist length Lw), the pulling process is divided into N cycles.

During each cycle, the radius of the fiber cladding is reduced by a constant ratio to reach the nanofiber waist radius at the end of the last cycle and the taper is elongated. At the n th cycle, the heat zone length is chosen as a power law of n, enabling to create few-cm length tapers. This fabrication process is performed in a class-5 cleanroom to avoid the presence on the dust. As soon as the nanofiber is immersed in the Raman liquid, the experiment can be carried on in laboratory conditions.

Based on the values of 𝑔 𝑅𝑚 represented in Figure 3, the converter's radii are chosen to be between 220 nm and 350 nm, which corresponds to a modal Raman gain between 0.35 m-1. W-1 and 0.98 m-1. W-1. Lw is fixed to 8 cm which is the highest length we can have with our pulling system.

We choose to pull the nanofiber from two different fibers : a single mode fiber at 532 nm (460HP from Thorlabs) and a single mode fiber at 1550 nm (SMF28 from Corning). Since the SMF28 fiber is highly multimode at p = 532 nm, high order modes should be filtered to allow only the propagation of the fundamental mode HE11. One way to filter higher order modes is to pull L1 over a few-cm length until reaching the cutoff frequency of the first higher order mode and to keep the fiber short and straight after the filter.

Figure 3. 𝑔 𝑅𝑚 versus nanofiber radius for HE11, TM01, TE01 and HE21.

To guarantee high light transmission, special attention should be paid to the tapers adiabaticity. The tapering angle has to be small enough to avoid coupling between the fundamental and higher order modes. We study the adiabaticity for the fundamental mode LP01. The adiabaticity criterion that is usually adopted is [11]: is the beating length. 𝛽 1 is the local propagation constant of the pump mode LP01 and 𝛽 2 is the local propagation constant of the mode having the highest chance to couple with it because of the azimuthal dependence, i.e. LP02.

Fixing L w to 8 cm, we consider three nanofiber samples : NF1 (R = 220 nm), NF2 (R = 300 nm) and NF3 (R = 350 nm). The tapers created for these three samples are the same for both fibers SMF28 and 460HP with lengths Lt respectively equal to 5 cm, 4.5 cm and 4 cm. The normalized slopes of the three tapers profiles are plotted as a function of the relative radius in Figure 4. The adiabaticity boundaries at 532 nm and 630 nm for both fibers are also represented.

Figure 4 shows that making the taper adiabatic with SMF28 is easier than with 460HP for the fabricated taper profiles. Indeed, the adiabaticity boundary for the SMF28 taper curve has a punctual minimum of 0.0035 m-1 for a relative radius of 0.15. This minimum is not as punctual for the 460HP taper and it is maintained over a relative radius range extending from 0.3 to 0.55 keeping a value of 0.0019 m-1. When the nanofiber is immersed in ethanol, the light transmission of the whole component at 532 nm defined as the ratio between the output energy 𝐸 𝑝,𝑜𝑢𝑡 and the input energy 𝐸 𝑖𝑛 is experimentally 90% for the SMF28 and 87% for the 460HP. These transmission values are almost the same for the three taper profiles. The adiabatic boundary curves are right-shifted when the wavelength increases which keeps the tapers adiabatic at the Stokes wavelength of ethanol S = 630 nm.

NONLINEAR CHARACTERIZATION OF THE CONVERTER

In order to obtain the maximal Raman conversion efficiency, controlling the undesired nonlinear effects is mandatory. The experimental setup presented in Figure 5. A pulsed pump laser (from HORUS) emits a beam at 532 nm with a maximum available peak pump power of 7 kW, a pulse duration of 900 ps and a frequency rate repetition of 4.7 kHz. A lens L1 with a focal distance of 10 cm collimates the pump light and the association of a half plate and a polarizing beam splitter helps us to control the incident pump power since the laser is linearly polarized. Part of the pump beam then passes through a glass plate. The converter is immersed in a tank filled with the Raman liquid. The pump beam is injected by a microscope objective (x20, NA 0.35) in the untapered input end of the fiber L1, propagates in the nanofiber and is collected at the output end of the fiber L2. The role of the glass plate is to reflect the backward beams into an optical spectrum analyzer OSA2 while OSA1 will analyze the forward beams. (

We depict in Figure 6 the spectra of the forward beam collected from the three considered nanofiber radii when the immersing liquid is ethanol. We observe that there is no forward silica Raman scattering for the whole converter at 545 nm. This is explained by the appropriate choice of L1 and L2 lengths (respectively 30 cm and 40 cm) that are short enough and the silica Raman nanofiber threshold energy that is higher than the evanescent Raman threshold energy in ethanol.

No significant spectral broadening of the pump line is observed showing the absence of Kerr effect in silica and in ethanol [12].

In the sub-nanosecond regime, as the effective length of interaction and the pump length are of the same magnitude ( ~ 10 𝑐𝑚), backward nonlinear effects can be favored [13]. On OSA2, we observed backward Raman scattering, neither in silica nor in ethanol. We measured the backward energy at 532 nm versus the incident pump energy. The plots are almost linear for the three nanofibers, showing that there is not Brillouin scattering. The backward energy is due to parasitical reflections.

CONVERSION EFFICIENCY AND REPRODUCIBILITY STUDIES

Conversion efficiency study of the converter

Ensuring high modal Raman gains 𝑔 𝑅𝑚 is not sufficient to guarantee a high conversion efficiency of the evanescent Raman converter. The two key parameters that influence the conversion efficiency of the evanescent Raman converter are the Raman threshold energy ER, th and the nanofiber optical breakdown energy Ebreak [7]. We define the gap between the two parameters by "the Raman conversion operating range" Δ𝐸. The higher Δ𝐸 is, the more efficient the converter becomes [7].

We estimate in Table 1 the optical breakdown energies of the nanofiber samples NF1, NF2 and NF3 for both fibers. We present the results for two different external mediums (air and ethanol). Almost the same values are found with for both fibers. In air, the optical breakdown energy increases with the radius. In ethanol, the optical breakdown energy becomes higher for the three radii. We focus on the three SMF28 nanofibers. We compare the experimental and the theoretical Raman threshold energies. The theoretical Raman threshold energies are deduced from the expression of . The experimental and the theoretical Raman thresholds for NF1 and NF2 are in a good agreement and confirm the value of the critical parameter   23 obtained in [6,7]. We were not able to reach the Raman threshold with NF3. A longer nanofiber would be needed. Theoretical and experimental Raman conversion operating ranges Δ𝐸 are represented in Table 2.

Based on the estimated Raman conversion operating ranges Δ𝐸, the evanescent Raman converter created with NF2 is expected to be more efficient than the one fabricated with NF1. We present in Figure 8 the transmitted pump energy 𝐸 𝑝,𝑜𝑢𝑡 and the output Stokes energy 𝐸 𝑠,𝑜𝑢𝑡 as a function of the incident pump energy 𝐸 𝑖𝑛 for the three considered nanofibers (SMF28-NF1, SMF28-NF2, SMF28-NF3).

We define the external efficiency of the Raman conversion ex as the fraction between the maximum output Stokes energy and the incident pump energy after reaching the Raman threshold. As expected, the highest Stokes energy was obtained with SMF28-NF2 and ex = 60%. An energy of 0.29 J at 630 nm is obtained just before the SMF28-NF2 broke. The Raman threshold is reached at 𝐸 𝑖𝑛 = 0.41J. For SMF28-NF1, the maximum extracted Stokes energy is only 0.13 J and ex is 51%. With SMF28-NF3, a Stokes energy of only 0.15 J for an incident energy 𝐸 𝑖𝑛 = 0.58 J is achieved and the Raman threshold is not reached.

Reproducibility study of the converter

In order to study the reproducibility of the converter, we draw 20 samples of each nanofiber (SMF28-NF1, SMF28-NF2, SMF28-NF3). Figure 9 represents ER,th, Ebreak and Esmax for the 20 samples of SMF28-NF2. The results of SMF28-NF1 and SMF28-NF3 show the same trends. Table 2. Raman conversion operating ranges in 𝜇J for the tested nanofibers.

We also compute the average and the standard deviation of ESmax, ER,th, Ebreak for the 20 samples of SMF28-NF2 just before the breakdown in Table 3.

CONCLUSION

CONCLUSION

We have proposed evanescent Raman converters in the sub-nanosecond regime based on a silica nanofiber immersed in ethanol. Two different standard silica fibers were tested and compared (SMF28 and 460HP). We obtained a maximum output Stokes energy of 0.29 J with an external conversion efficiency from the pump to the first Stokes order of ethanol of 60%. One key parameter to optimize the extracted Stokes energy is the Raman conversion operating range, which is the gap between the Raman threshold energy and the optical breakdown energy. Indeed, the optical breakdown energy is estimated for different nanofiber diameters and for different external mediums. The Raman conversion operating range is a guideline for the conception of other evanescent Raman converters based on other liquids and/or fibers. The performances of the evanescent Raman converters are highly reproducible and other perspectives can be studied like the deposition of nonlinear polymers on the nanofiber surface opening the way to a new family of robust all-fibered components that can be directly inserted in optical networks without additional losses.

Figure 4 .

 4 Figure 4. Tapers immersed in ethanol at 532 nm for the SMF28 (empty round) and 460HP (empty square). Same at 630 nm for the SMF28 (cross) and 460HP (star).

Figure 5 .

 5 Figure 5. Experimental setup scheme for the nonlinear characterization of evanescent Raman converter.

Figure 7 .

 7 Figure 7. Backward energy versus the incident pump energy for the three studied nanofibers.

Figure 6 .

 6 Figure 6. Collected forward spectra from OSA1 for nanofibers immersed in ethanol just before they break.

Figure 8 .

 8 Figure 8. Ep, out and ES, out versus the laser incident energy 𝐸 𝑖𝑛 for SMF28-NF1, SMF28-NF2 and SMF28-NF3.

Figure 9 .

 9 Figure 9. Experimental Raman threshold energies, incident nanofiber breakdown energies and maximal Stokes energies for the 20 samples of SMF28-NF2.

Table 1 .

 1 Optical breakdown energy in 𝜇J for different nanofibers diameters in different external mediums.

Rap. Publ. 13, 31 (2017).

[13] G. Fanjoux, J.