A Spectra Classification Methodology of Hyperspectral Infrared Images for Near Real-Time Estimation of the SO2 Emission Flux from Mount Etna with LARA Radiative Transfer Retrieval Model

: Fast and accurate quantification of gas fluxes emitted by volcanoes is essential for the risk mitigation of explosive eruption, and for the fundamental understanding of shallow eruptive processes. Sulphur dioxide (SO 2 ), in particular, is a reliable indicator to predict upcoming eruptions, and its systemic characterization allows the rapid assessment of sudden changes in eruptive dynamics. In this regard, infrared (IR) hyperspectral imaging is a promising new technology for accurately measure SO 2 fluxes day and night at a frame rate down to 1 image per second. The thermal infrared region is not very sensitive to particle scattering, which is an asset for the study of volcanic plume. A ground based infrared hyperspectral imager was deployed during the IMAGETNA campaign in 2015 and provided high spectral resolution images of the Mount Etna (Sicily, Italy) plume from the North East Crater (NEC), mainly. The LongWave InfraRed (LWIR) hyperspectral imager, hereafter name Hyper-Cam, ranges between 850–1300 cm − 1 (7.7–11.8 µ m). The LATMOS (Laboratoire Atmosph è res Milieux Observations Spatiales) Atmospheric Retrieval Algorithm (LARA), which is used to retrieve the slant column densities (SCD) of SO 2 , is a robust and a complete radiative transfer model, well adapted to the inversion of ground-based remote measurements. However, the calculation time to process the raw data and retrieve the infrared spectra, which is about seven days for the retrieval of one image of SO 2 SCD, remains too high to infer near real-time (NRT) SO 2 emission ﬂuxes. A spectral image classiﬁcation methodology based on two parameters extracting spectral features in the O 3 and SO 2 emission bands was developed to create a library. The relevance is evaluated in detail through tests. From data acquisition to the generation of SO 2 SCD images, this method requires only ~40 s per image, which opens the possibility to infer NRT estimation of SO 2 emission ﬂuxes from IR hyperspectral imager measurements.


Introduction
More than 500 million people live within the potential exposure range of a volcano [1]. Monitoring volcanoes is then essential and involves different types of measurements, such as volcanic degassing, seismicity, and ground deformation detection [2]. Volatiles, in particular, are a crucial component of volcanic systems. The explosivity of an eruption depends, in large part, on the amount and composition of gas contained in the erupted magma. For that reason, measurements of volcanic degassing have been one of the most widely used methods in volcanic monitoring networks for more than 40 years. The shallow exsolution depth of SO 2 , compared to carbon dioxide (CO 2 ), makes it a good indicator of the presence of a magmatic body near the surface, and, therefore, provides a tool to forecast eruptions. Fluctuations in degassing levels may reflect changes in the magma supply rate and help inform a short-term forecast of on-going eruptions [3][4][5][6][7][8]. Moreover, the composition of volcanic gases offers insight into physical processes occurring at depth [9,10]. In addition to the risks induced by volcanic eruptions, monitoring the volcanoes degassing emissions in the atmosphere is important for an environmental impact and a hazardous effect on human health. After H 2 O and CO 2 , sulfur dioxide is the main volcanic gas emitted by volcanoes. It is widely monitored due to its low background atmospheric abundance. The impact of SO 2 emission on atmospheric chemistry and the estimation of the global budget of SO 2 in the atmosphere induced by volcanoes are part of climate and environmental monitoring as well as public health prevention [11]. The oxidation of SO 2 leads to the formation of sulfur aerosols responsible for acid rain [12], problems on vegetation growth close to volcanoes, and cause asthma or respiratory problems to humans [13,14].
Mount Etna mean SO 2 flux is equivalent, for example, to the total anthropogenic SO 2 emissions from France [12]. It is the largest active volcano in Europe and the world's strongest permanent contributor of volcanic volatiles to the atmosphere [15], which makes it one of the most monitored volcanoes of the planet [2]. The accurate characterization of its gas emission fluxes is, hence, the purpose of many studies.
Volcanic gas emissions are mostly monitored by satellites as well as by ground-based instruments. The satellite-based instruments remotely measuring SO 2 use different channels and absorption bands such as OMI (Ozone Monitoring Instrument) [16,17] and TROPOMI (Tropospheric Monitoring Instrument) [18,19] in the ultraviolet (UV), MODIS (Moderate-Resolution Imaging Spectroradiometer), SEVIRI (Spinning Enhanced Visible and Infrared Imager) [20,21], IASI (Infrared Atmospheric Sounding Interferometer), and AIRS (Atmospheric InfraRed Sounder) [22,23] in the infrared (IR). Ground-based monitoring is performed remotely using ultraviolet (UV) cameras [24], mini-DOAS (Differential Optical Absorption Spectroscopy) instruments [25,26], and IR cameras [27]. The NOVAC (Network for Observation of Volcanic and Atmospheric Change) project with the deployment of fully automated mini-DOAS systems used in a scanning mode has permit to get the first real-time determinations of SO 2 flux in a monitoring system [28]. This system allows the monitoring of nearly 20 volcanoes around the world [29]. In 2004, the FLAME (FLux Automatic MEasurements) network was developed and installed by the Istituto Nazionale Geofisica e Vulcanologia (INGV) (Catania, Sicily) on Mount Etna flanks with nine UV scanning spectrometers to automatically get SO 2 fluxes in real-time and with a high frequency [30]. The network collects data daily during the daytime with a complete scan everỹ 5 min [31]. The SO 2 flux emissions of the different vents of Mount Etna from the FLAME network reveals shifts in the activity of each vent approaching an eruptive episode [5,32].
Ground-based hyperspectral infrared imaging systems of volcanic degassing offer high spatial, spectral, and temporal resolution measurements [33,34]. Such imagers can provide continuous day and night measurements, which is likely the most important advantage over UV imagers, but also essential to monitor volcanoes located in high latitude regions during the wintertime. Commercial hyperspectral IR imagers are still an expensive technology, which may be more attractive in the near future with the increase of possible applications and the development of studies and algorithms aiming at exploiting all the capacities of such instruments. Nevertheless, some recent advances on the development of uncooled detector technology allow the use of cheap thermal infrared imagers to quantify hazardous gases. Wright et al. (2013) presented the characteristics of a thermal hyperspectral imager with~30 spectral bands in the 8-14 µm wavelength region [35].
Hyperspectral imaging is an interesting and powerful tool in a wide range of applications such as agriculture, forestry, and environmental management [36], geological exploitation and mineralogy, urban planning, and disaster prediction, military and defense [37], but also remote sensing of volcanic plume for gas emissions monitoring. The vast amount of data comprised in the spectra of the images obtained with this kind of instrument makes it heavy to process and induce to set up a methodology for classifying the pixels of an image to get the most information of the spectral signature in near real-time. Different deep-learning methods exist, including supervised or unsupervised, parametric or non-parametric ones. Paoletti et al. (2019) reviewed the different popular models and architectures used for the classification of remotely sensed hyperspectral spectral images [38]. To achieve a successful classification of the image pixels, many factors and steps are involved, starting with the determination of a suitable classification system, selection of training samples, image preprocessing, feature extraction, selection of suitable classification approaches, post-classification processing, and accuracy assessment [39].
When pixels of images are spectra like in this study, the goal of image pixels classification is to identify different patterns in spectra with an automated methodology in order to extract information from data and create a set of reference spectra. Many classification approaches have been developed depending on the domain of application. For example, Van Damme et al. (2017) present a version of a neural-network-based algorithm developed to retrieve in near real-time the ammonia columns from the satellite observations from the hyperspectral imager IASI [40].
A methodology to analyze the large amount of data produced by such technology is necessary to respond to the need of NRT knowledge of volcanic gas emission fluxes.
In this paper, we present a spectral image classification methodology for infrared hyperspectral images from Mount Etna volcanic plume. The objective is to achieve a fast and accurate retrieval of volcanic plume in order to get the SO 2 emission flux in NRT, which is essential for early warning and mitigation of volcanic risk.
In the next section, we describe the dataset from the IMAGETNA field campaign and the LARA model (LATMOS (Laboratoire Atmosphères Milieux Observations Spatiales) Atmospheric Retrieval Algorithm) used to retrieve the data, and we present the method. More detailed information about the IMAGETNA and LARA model are presented in Huret et al. 2019 [34]. Then, we bring step-by-step the pixels classification methodology that we applied to our dataset. The results of the classification and the flux estimations are presented in Section 3 and discussed in Section 4. The last section is dedicated to the conclusions and perspective of this study.

Overview of the Data
The IMAGETNA campaign was held from 21-26 June 2015 from the Pizzi De Neri Volcano observatory on the North side of Mount Etna at 2850 m altitude and at~2 km from the plume. One of the purposes of the field campaign was to make remote observations of the plume of Mount Etna using the infrared imager Hyper-Cam developed by Telops company, during a quiescent stage of its activity. The first days of the campaign were dedicated to instrumentation deployment and basic tests. The weather during the campaign was very good, with mainly cloudless skies, allowing us to carry out a large number of measurements in good conditions. The Hyper-Cam is a passive infrared hyperspectral imaging system with a spectral range of 840-1300 cm −1 (7.7-11 µm). This camera combines a high spatial resolution with images containing up to 320 × 256 pixels, a maximum field of view of 6.4 • × 5.1 • , and spectral resolutions from 0.25 to 150 cm −1 . Considering the geographical arrangement of the campaign and the parameterization of the hyperspectral imager in relation to the Etna plume, the surface of each pixel is 7.84 m 2 (i.e., 2.8 × 2.8 m). The Hyper-Cam is radiometrically calibrated using two blackbodies with different temperatures with a spectral radiometric accuracy of 0.5 K at a blackbody temperature of 30 • C [41]. Each measured image has three dimensions: two spatial dimensions and one spectral dimension.
As shown in Figure 1a, Mount Etna has several active craters pointed out with red arrows: the two central craters (Bocca Nuova (BN) and Voragine (V)), the North East crater (NEC), and the South East craters (SEC), which is the one created most recently [42]. To get the maximum thermal contrast between the plume and the background, the measurements were done early in the morning starting at around 8:00 UTC. Considering the direction of the plume propagation and the location of the observatory where the measurements have been done, the plume captured by the Hyper-Cam imager during the field campaign is mainly related to the North East Crater (NEC) of Mount Etna. An example of projection of the contours of the acquired images is superimposed on Figure 1b. This projected image corresponds to an image of [320 × 64] pixels.
Remote Sens. 2020, 12, x FOR PEER REVIEW  4 of 21 with a spectral radiometric accuracy of 0.5 K at a blackbody temperature of 30 °C [41]. Each measured image has three dimensions: two spatial dimensions and one spectral dimension. As shown in Figure 1a, Mount Etna has several active craters pointed out with red arrows: the two central craters (Bocca Nuova (BN) and Voragine (V)), the North East crater (NEC), and the South East craters (SEC), which is the one created most recently [42]. To get the maximum thermal contrast between the plume and the background, the measurements were done early in the morning starting at around 8:00 UTC. Considering the direction of the plume propagation and the location of the observatory where the measurements have been done, the plume captured by the Hyper-Cam imager during the field campaign is mainly related to the North East Crater (NEC) of Mount Etna. An example of projection of the contours of the acquired images is superimposed on Figure 1b. This projected image corresponds to an image of [320 × 64] pixels. A trade-off is necessary between the high spectral and spatial resolutions and the temporal resolution to capture the dynamism of the measured scene. Table 1 presents the main characteristics of the three sequences reported in this study such as spectral resolution, from 2 to 4 cm −1 , and the acquisition time, from ~1.3 to ~4.6 s, resulting from the trade-off made for each sequence. All the images have 20,480 pixels (i.e., 320 × 64 pixels). The field of view includes part of the ground, the plume, and the clear sky. The dataset is composed of ~900 images representing over 18 million of pixels and ~35 min of measurements.
The three sequences, to investigate SO2 emissions fluxes, have been selected among the other sequences because of their cloud-free conditions but also because they have a common measurement geometry with an unchanged location and distance of the camera with respect to the plume. Moreover, the meteorological conditions are close in term of temperature and wind direction, which is roughly horizontal southeasterly. With such similarities, we can extract information from the data by developing data processing to the three cases. The acquisition times are of the order of a few seconds, allowing us to finely capture the dynamics of the plume and, hence, the SO2 emission flux. A trade-off is necessary between the high spectral and spatial resolutions and the temporal resolution to capture the dynamism of the measured scene. Table 1 presents the main characteristics of the three sequences reported in this study such as spectral resolution, from 2 to 4 cm −1 , and the acquisition time, from~1.3 to~4.6 s, resulting from the trade-off made for each sequence. All the images have 20,480 pixels (i.e., 320 × 64 pixels). The field of view includes part of the ground, the plume, and the clear sky. The dataset is composed of~900 images representing over 18 million of pixels and 35 min of measurements.   The three sequences, to investigate SO 2 emissions fluxes, have been selected among the other sequences because of their cloud-free conditions but also because they have a common measurement geometry with an unchanged location and distance of the camera with respect to the plume. Moreover, the meteorological conditions are close in term of temperature and wind direction, which is roughly horizontal southeasterly. With such similarities, we can extract information from the data by developing data processing to the three cases. The acquisition times are of the order of a few seconds, allowing us to finely capture the dynamics of the plume and, hence, the SO 2 emission flux.

Retrieval Model: LARA
The data retrieval of SO 2 slant column densities has been performed with the LATMOS Atmospheric Retrieval Algorithm (LARA). LARA is a radiative transfer model associated with a minimization algorithm of the Levenberg-Marquardt type [43], which has been developed and adapted, in the last decade [44,45], to both nadir and limb geometries, and to balloon, satellite, or ground-based experiments [46][47][48]. It takes into account the geometry of the plume, the thickness and the altitude, and different chemical species, H 2 O, CO 2 , O 3 , CH 4 , NO 2 , and SO 2 . The radiative transfer is calculated in each pixel with a profile along the line of sight divided in more than 40 layers from the instrument altitude of 2.85 km up to the top of the atmosphere at 100 km. The model also considers the contribution of particles in the IR spectra as well as the water vapor continuum [49]. Huret et al. 2019 [34] present the different tests to configure the model and results of the SO 2 slant column densities (SCD) retrievals. The sensitivity of the retrieval to the parameters of the geometry of the plume, such as the plume altitude, its thickness, and the additional temperature with respect to ambient temperature, highlighted a low influence of those parameters on the SO 2 SCD in the dilute part of the plume. That study demonstrated that LARA is a robust and well adapted model to ground-based remote sensing of volcanic plume measurements.
The main downside of LARA is the calculation time to retrieve each image pixel-by-pixel. For example, the retrieval of an image of [320 × 64] pixels with 60 spectral points in the SO 2 spectral band requires one week of calculation on the computers of AERIS data center (www.aeris-data.fr) to be able to image the SO 2 slant column densities. This method is very expensive in term of calculation time and cannot be applied to process all other sequences from the IMAGETNA campaign. To make the most of the capabilities of this kind of IR hyperspectral imager and to reach quasi real-time estimation of the SO 2 flux, developing a methodology to significantly reduce the calculation time of SO 2 SCD retrieval is crucial.
Five images of one of the sequences measured on 26 June 2015 have been previously processed pixel-by-pixel using the LARA model for the study presented in Huret et al. 2019 [34]. Already retrieved images of SO 2 SCD (ppm m) can easily be converted into mass of SO 2 per surface unit (g m −2 ) using Equation (1).

Massive Retrieval Methodology
The objective here is to define the architecture of a network, based on characteristic features in the brightness temperature (BT) spectra, from which information leading to fast and accurate determination of the SO2 SCD in each pixel can be extracted.

Hyper-Cam Spectral Band Analysis
The spectral range of the Hyper-Cam encompasses in particular the O3 and the SO2 emission bands at (1000-1100) cm −1 and (1100-1200) cm −1 , respectively, as shown in Figure 3a. The emission line intensities of O3, SO2, and H2O extracted from the GEISA (Gestion et Etude des Informations Spectroscopiques Atmosphériques) spectroscopic database and simulated with a GEISA graphical tool (https://geisa.aeris-data.fr/), with an average spectral resolution greater than 0.2 cm −1 , are presented in Figure 3a,b. Figure 3a shows well-defined spectral emission regions for O3 and SO2, which are independent from each other and can, therefore, provide specific information from either spectral window. The signal due to O3 emission lines can be considered as two continuous well identified branches. The SO2 emission line intensities are less intense, with two flatted wings. H2O has several emission lines in the SO2 spectral range [1100-1200] cm −1 (see Figure 3b), which contribute to the signal. Since the water vapor column is retrieved simultaneously to SO2, H2O emission lines are then well fitted. Aerosols and droplets are parameterized and considered as a first-order modelling with a spectral dependency of the plume optical thickness (see Huret et al. 2019 for more details) [34]. Figure 3c presents different examples of BT spectra extracted from the images of the third sequence of 26 June 2015 (sequence B in Table 1). The presence of O3, SO2, and H2O is recognizable with more or less intensity in those example spectra depending on the part of the image they are related to. The spectral signatures in those images denote the ground in green, the sky without plume in red, the plume with different levels of dilution (in orange and blue), with more or less water vapor, droplets, and aerosols. The brightness temperature is the temperature a black body would have if it has the same radiance as the one measured by Hyper-Cam.

Massive Retrieval Methodology
The objective here is to define the architecture of a network, based on characteristic features in the brightness temperature (BT) spectra, from which information leading to fast and accurate determination of the SO 2 SCD in each pixel can be extracted.

Hyper-Cam Spectral Band Analysis
The spectral range of the Hyper-Cam encompasses in particular the O 3 and the SO 2 emission bands at (1000-1100) cm −1 and (1100-1200) cm −1 , respectively, as shown in Figure 3a. The emission line intensities of O 3 , SO 2 , and H 2 O extracted from the GEISA (Gestion et Etude des Informations Spectroscopiques Atmosphériques) spectroscopic database and simulated with a GEISA graphical tool (https://geisa.aeris-data.fr/), with an average spectral resolution greater than 0.2 cm −1 , are presented in Figure 3a,b. Figure 3a shows well-defined spectral emission regions for O 3 and SO 2 , which are independent from each other and can, therefore, provide specific information from either spectral window. The signal due to O 3 emission lines can be considered as two continuous well identified branches. The SO 2 emission line intensities are less intense, with two flatted wings. H 2 O has several emission lines in the SO 2 spectral range [1100-1200] cm −1 (see Figure 3b), which contribute to the signal. Since the water vapor column is retrieved simultaneously to SO 2 , H 2 O emission lines are then well fitted. Aerosols and droplets are parameterized and considered as a first-order modelling with a spectral dependency of the plume optical thickness (see Huret et al. 2019 for more details) [34]. Figure 3c presents different examples of BT spectra extracted from the images of the third sequence of 26 June 2015 (sequence B in Table 1). The presence of O 3 , SO 2 , and H 2 O is recognizable with more or less intensity in those example spectra depending on the part of the image they are related to. The spectral signatures in those images denote the ground in green, the sky without plume in red, the plume with different levels of dilution (in orange and blue), with more or less water vapor, droplets, and aerosols. The brightness temperature is the temperature a black body would have if it has the same radiance as the one measured by Hyper-Cam.

O 3 Emission Region: Spectral "Band 1"
The first spectral band considered here is situated in the range of 1000-1100 cm −1 characteristic of O 3 . Different examples of BT spectra from the measurements on 26 June 2015 with a spectral resolution of 2 cm −1 are given on Figure 3c. We can easily identify the ozone emission region, defined as Band 1, thanks to its typical rovibrational branches in the clear sky, the diluted plume, and the less diluted plume pixels, corresponding respectively to red, orange, and blue spectra. The ozone band in the clear sky is strong and characteristic of stratospheric ozone. The less diluted and opaque the plume is, the flatter the ozone branches are. As for the green spectra corresponding to ground pixels, the ozone is not relevant. Hence, the spectra are flat with an average BT value of 300 K. By calculating the integrated BT values in Band 1, we extract a spectral feature, identifying if we have a ground, a clear sky, or a plume pixel spectrum. The lowest values correspond to a clean atmosphere (i.e., plume free atmosphere),~22,950 K cm −1 for red spectra in Figure 3c, and the highest values to a very opaque plume or even ground pixel,~26,450 and~29,950 K cm −1 for blue and green spectra in Figure 3c, respectively.

SO 2 Emission Region: Spectral "Band 2"
The second spectral band considered here is in the range of 1100-1200 cm −1 , which is the SO 2 spectral band. Inside Band 2 in Figure 3c, we can see that the clear sky BT is the lowest with an average value of 200 K. A bump shape of the spectra in the diluted and less diluted pixels (orange and blue spectra, respectively) highlights the presence of SO 2 in the plume and the few more intense peaks are significant of the water vapor emission lines. The BT in the SO 2 spectral region increases along with the increase of SO 2 concentration inside the plume. By calculating the mean BT values in Band 2, we extract another spectral feature from the data, which is characteristic of the SO 2 concentration in the spectra. For the images of the IMAGETNA campaign previously presented in Table 1, the plume mean BT varies between 200 and 300 K in the clear sky and ground, respectively.
The architecture of the classification network is composed of two layers. The first layer categorizes the pixels with regard to their main characteristic, which is either a ground, plume, or clear sky pixel using Band 1. An index I O 3 corresponding to the integrated value of the BT in the O 3 spectral region is attributed to each spectrum. Then, the second layer gives information about the concentration in SO 2 . The second index, I SO 2 , attributed to each spectrum corresponds to the mean BT in the SO 2 spectral region. Each spectral feature is extracted and classified following the mapping function Y = f I O 3 , I SO 2 , which assigns a vector label to each pixel.

SO 2 Emission Flux Estimation
In this section, we present the chosen method to determine the Etna plume transport speed by analyzing the images and then the most adapted method used to estimate the SO 2 emission flux with regard to our campaign configuration.

Plume Transport Speed
A key parameter to estimate the SO 2 emission flux is the plume transport speed. Such data has not been measured during the remote measurement campaign. Each sequence contains a few hundred images, corresponding to a few minutes of measurement. Following the method explained in Aiuppa et al. (2015) [50], it is possible to determine an average transport speed after determining its vertical and horizontal components using the reconstructed images of SO 2 SCD (or mass of SO 2 per surface unit). Two cross sections are fixed vertically, separated by a distance d v and two horizontally separated by a distance d h , corresponding to the number of pixels between the cross sections multiplied by the length of a pixel. By calculating the mean SO 2 SCD along the cross sections, we obtain two time-series, one vertically and one horizontally tracking the evolution over time of SO 2 SCD in the plume. The cross correlation of both time-series let us determine the time shift t s , which represents the time necessary for the plume passing through the first vertical (or horizontal) cross section to reach the second cross section. The appropriate time shift is the one giving the best cross-correlation coefficient between the two time-series. A representation of the cross sections is presented in Figure 4a. correlation coefficient between the two time-series. A representation of the cross sections is presented in Figure 4a.    Figure 4b,c present the SO 2 mean SCD values evolutions along 200 images of sequence B, for the vertical cross sections separated by a distance d v corresponding to 5 pixels (i.e.,~14 m) and the horizontal cross sections separated by a distance d h corresponding to 20 pixels (i.e.,~56 m), to characterize the plume transport speed. The plume transport speed v in m s −1 is then calculated by following Equation (2).

Box Method for the Emission Flux Estimation
Different methods to derive the SO 2 emission flux from the SO 2 slant column densities (or mass of SO 2 ) exist such as a traverse method, delta-M method, box method, or an inverse modelling method [51]. Most methods are adapted to satellite or airborne measurements of volcanic plumes. Instantaneous SO 2 mass loading measurements from remote instruments represent the budget between emitted and lost SO 2 . Therefore, when calculating gas fluxes, the SO 2 loss has to be taken into account. The loss is mainly due to SO 2 oxidation into sulfuric acid droplets and SO 2 volume dilution, making its detection below a given threshold impossible. The rate of gas loss is related to the meteorological conditions and its inclusion may depend on the measurement settings. In our case, measurements were carried very close to the source vent by using high temporal and spatial resolutions, and located to altitudes larger than 3 km above sea level, close to the free troposphere. These conditions of acquisition allow us to neglect the SO 2 loss factor [52]. One of the simplest and most adapted method to derive the first estimation of the emission flux of SO 2 is the box method. It consists in defining a box, as presented with the yellow box in Figure 4a, calculating the total mass of SO 2 in the box and then multiplying it by the plume age. The plume age is the time the plume needs to cross the box. It is obtained by multiplying the box length with the average plume transport speed v. In order to lower the uncertainties of the SO 2 emission flux estimation, the fixed box is placed as close as possible to the emission source.

Training Dataset
The spectral image classification methodology was conducted using five images of the third sequence recorded on 26 June 2015. Those images have been retrieved pixel-by-pixel with the LARA model for a previous study presented in Huret et al. (2019) [34]. They are used here as a training dataset composed of 102,400 spectra. Different ranges of values assigned to the two indexes were tested to determine the best compromise between the number of classes obtained, which means the number of different patterns of spectra, and the range size limit to keep as much information as possible by classifying the images.

Interval Width for I O 3 Index
Four ranges of values for the integrated BT values have been tested: 10, 50, 100, and 500 K cm −1 . The main results obtained for each tested value with a fixed value of 1 K for the range of mean BT in the spectral Band 2 is presented in Table 2. The only I O 3 range value showing a significant difference is 500 K cm −1 with only~26% of the pixels with less than 10% of relative standard deviation for plume pixels against~48% for ranges of 10, 50, and 100 K cm −1 . Three ranges of values for the mean BT have been tested: 1, 5, and 10 K. The main results obtained for each value with a fixed range value of 100 K cm −1 for the integrated values in the spectral Band 1 are presented in Table 3. The tested range values for the mean BT, attributed to I SO 2 , only gives satisfying results for the range value of 1 K. The ranges of 5 and 10 K induce an important loss of pixels in classes with less than 10% of a relative difference in the SO 2 standard deviation in classes. The chosen spectral image classification configuration with the best compromise is a range of 100 K cm −1 , between two values of integrated BT for I O 3 , and a range of 1 K, between two values of mean BT for I SO 2 . Figure 5a,b represent the weight attributed to each index I O 3 and I SO 2 , respectively. The different bumps in the distributions highlight three main features from the images, which we can call clusters. The first cluster (1) corresponds to the clear sky pixels with I O 3 index < 25 and I SO 2 index < 30, which is equivalent to an integrated value of the BT below 23,500 K cm −1 and a mean BT below 230 K. The second cluster (2) encompasses the plume pixels with I O 3 indexes in the range of 25 to 75 and I SO 2 indexes in the range of 30 to 90, equivalent to integrated BT values in the range of [23,500:28,500] K cm −1 and mean BT in the range of [230:290] K. Then, the third cluster (3) collects the ground pixels with I O 3 index > 75 and I SO 2 index > 90, equivalent to an integrated BT greater than 28,500 K cm −1 and a mean BT above 290 K. The distribution of the number of pixels for each I SO 2 is wider in the plume cluster than the distribution of pixels for I O 3 . The pixels have a thinner distribution along I SO 2 indexes with a total number of pixels around 3000 for the higher peaks while the higher peak of the total number of pixels along I O 3 reaches 4000. The diversity of I SO 2 values for one I O 3 is greater when considering the plume part of the image. This highlights the added value of the I SO 2 indexes to get information about SO 2 SCD from data using the two layers classification methodology.  [34], the retrieval of SO2 SCD must respect two threshold values to be considered of good quality: a normalized χ² < 10 and an error on the retrieved scaling factor on the SO2 SCD calculated by the LARA model that must be lower than 10%. Taking into account those threshold values, the SO2 SCD corresponding to each pair of ( , ) is calculated by taking the average value in each class of the SO2 SCD obtained with the pixel-by-pixel method for the training dataset. Around 50% of the classes corresponding to the plume part have a relative standard deviation of <10% and ~90% have a relative standard deviation of <15%. Some classes only contain a couple of pixels from the five images. To ensure that a class corresponds to a sufficiently characteristic spectral pattern of the plume, we consider that each class must contain at least 0.01% of the total number of plume pixels from the dataset. Hence, the classes from the training dataset must contain at least seven pixels. After applying this criterion, we have 195 remaining classes corresponding to ~68,000 pixels (99.7% of the plume pixels). In the end, 66% of the pixels from the training dataset are used to create a library to reconstruct images of SO2 SCD (ppm m) or SO2 mass per surface unit (g m −2 ) for the other sequences of Table 1. The other 34% corresponds to clear sky pixels, earth surface pixels, and plume pixels with LARA retrieval parameters above threshold values.

Analysis of the Classification Accuracy
A way to evaluate the accuracy of the classification using the five images of the training dataset is to see how correlated are the SO2 SCD values from the pixel-by-pixel retrieval and the SO2 SCD from the library created thanks to the classification, as shown in Figure 6. As presented in Huret et al. 2019 [34], the parameters used to define the Etna plume captured with Hyper-Cam (altitude, thickness, and an additional temperature with regard to ambient temperature) can have a significant impact on the SO2 SCD retrieval close to the crater. Hence, for the interpretation of the results, we applied the same procedure and separated the five images of the training dataset in two parts, "diluted plume part" and "dense plume part", along the same vertical line. The first order regression of the correlation points of the diluted plume part, in red, gives a slope of 0.97 following Equation (3). The 2D distribution of the number of pixels in each class (see Figure 5c) has a very linear shape. However, the third dimension corresponding to the weight of each class, let us find the three main features from the images with lighter blue contours. Most of the pixels are situated in the cluster related to the plume.
As presented in Huret et al. (2019) [34], the retrieval of SO 2 SCD must respect two threshold values to be considered of good quality: a normalized χ 2 < 10 and an error on the retrieved scaling factor on the SO 2 SCD calculated by the LARA model that must be lower than 10%. Taking into account those threshold values, the SO 2 SCD corresponding to each pair of (I O 3 , I SO 2 ) is calculated by taking the average value in each class of the SO 2 SCD obtained with the pixel-by-pixel method for the training dataset. Around 50% of the classes corresponding to the plume part have a relative standard deviation of <10% and~90% have a relative standard deviation of <15%. Some classes only contain a couple of pixels from the five images. To ensure that a class corresponds to a sufficiently characteristic spectral pattern of the plume, we consider that each class must contain at least 0.01% of the total number of plume pixels from the dataset. Hence, the classes from the training dataset must contain at least seven pixels. After applying this criterion, we have 195 remaining classes corresponding to~68,000 pixels (99.7% of the plume pixels). In the end, 66% of the pixels from the training dataset are used to create a library to reconstruct images of SO 2 SCD (ppm m) or SO 2 mass per surface unit (g m −2 ) for the other sequences of Table 1. The other 34% corresponds to clear sky pixels, earth surface pixels, and plume pixels with LARA retrieval parameters above threshold values.

Analysis of the Classification Accuracy
A way to evaluate the accuracy of the classification using the five images of the training dataset is to see how correlated are the SO 2 SCD values from the pixel-by-pixel retrieval and the SO 2 SCD from the library created thanks to the classification, as shown in Figure 6. As presented in Huret et al. 2019 [34], the parameters used to define the Etna plume captured with Hyper-Cam (altitude, thickness, and an additional temperature with regard to ambient temperature) can have a significant impact on the SO 2 SCD retrieval close to the crater. Hence, for the interpretation of the results, we applied the same procedure and separated the five images of the training dataset in two parts, "diluted plume part" and "dense plume part", along the same vertical line. The first order regression of the correlation points of the diluted plume part, in red, gives a slope of 0.97 following Equation (3).
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 21 = 0.97 × + 88.9 The correlation of the diluted plume part points is compact. When SO2 SCD values exceed 7000 ppm m, we observe a wider range of SO2 SCD values from the pixel-by-pixel retrieval attributed to one class. The correlation points over 7000 ppm m extending over a wider range of values represent only 2.7% of the plume pixels of the training dataset. The results obtained with the training dataset composed of five images are, thus, very satisfying. The library created with the 195 plume classes was applied to retrieve the images of the other sequences.

Tested Dataset
We apply the method to the 900 images of the three sequences previously presented in Table 1. Those sequences have been captured with similar instrument orientation as the plume of Mount Etna was transported in the same direction over the days of the campaign. Each image has ground pixels situated at the bottom right part of the images, clear sky pixels at the upper right part of the images, and the rest being plume pixels. All the sequences are classified independently, and each pixel is labeled with its corresponding and indexes. As mentioned in Section 2.2, the classes with less than 0.01% of the total number of plume pixels from the sequences are excluded. Moreover, since our interest is in plume pixels, we only consider classes corresponding to indexes with BT integral values in the range of [23,500:28,500] K cm −1 . Table 4 sums up the characteristics of the classification of the three tested datasets such as the time to classify an image, the proportion of patterns necessary to reconstruct the SO2 SCD images, and those who are not listed in the library. The time necessary to classify a sequence depends on three elements: the size of the image (i.e., number of pixels), the number of images in the dataset, and the spectral resolution. The retrieval is in the order of ~15 s for an image with a spectral resolution of 4 cm −1 and in the order of ~35 s for an image with a spectral resolution of 2 cm −1 . In comparison with the pixel-by-pixel method, which demands a week of calculation for an image of sequence B, the classification and retrieval using the library only takes ~34 s to retrieve an image from that same sequence. Figure 6. Example of correlation between the SO 2 SCD of the plume pixels of the training dataset of 26 June 2015 from the pixel-by-pixel method with the reconstructed values using the classification library. The black points are from the "diluted part of the plume" and the blue point are from the "dense part of the plume". The first order regression of the "diluted part of the plume" is in red with a slope of 0.97 and a determination coefficient of 0.94.
The correlation of the diluted plume part points is compact. When SO 2 SCD values exceed 7000 ppm m, we observe a wider range of SO 2 SCD values from the pixel-by-pixel retrieval attributed to one class. The correlation points over 7000 ppm m extending over a wider range of values represent only 2.7% of the plume pixels of the training dataset.
The results obtained with the training dataset composed of five images are, thus, very satisfying. The library created with the 195 plume classes was applied to retrieve the images of the other sequences.

Tested Dataset
We apply the method to the 900 images of the three sequences previously presented in Table 1. Those sequences have been captured with similar instrument orientation as the plume of Mount Etna was transported in the same direction over the days of the campaign. Each image has ground pixels situated at the bottom right part of the images, clear sky pixels at the upper right part of the images, and the rest being plume pixels. All the sequences are classified independently, and each pixel is labeled with its corresponding I O 3 and I SO 2 indexes. As mentioned in Section 2.2, the classes with less than 0.01% of the total number of plume pixels from the sequences are excluded. Moreover, since our interest is in plume pixels, we only consider classes corresponding to I O 3 indexes with BT integral values in the range of [23,500:28,500] K cm −1 . Table 4 sums up the characteristics of the classification of the three tested datasets such as the time to classify an image, the proportion of patterns necessary to reconstruct the SO 2 SCD images, and those who are not listed in the library. The time necessary to classify a sequence depends on three elements: the size of the image (i.e., number of pixels), the number of images in the dataset, and the spectral resolution. The retrieval is in the order of~15 s for an image with a spectral resolution of 4 cm −1 and in the order of~35 s for an image with a spectral resolution of 2 cm −1 . In comparison with the pixel-by-pixel method, which demands a week of calculation for an image of sequence B, the classification and retrieval using the library only takes~34 s to retrieve an image from that same sequence. Even though Mount Etna was in a quiescent stage of activity (i.e., not in eruption) and the plume was transported in the same direction, the training dataset represents only a part of the variety of the SO 2 spectral patterns of the Etna NEC plume. The library created allows the massive retrieval of three sequences with only 15% to~22% of unlisted plume pixels patterns, distributed in the whole sequence. One image example from each of the three datasets is presented in Figure 7 with the missing pixels colored in blue. As can be seen in Figure 7a, the sequence A has higher values of mass of SO 2 per surface unit in the range of~2.5-35 g m −2 in comparison with sequence B and C with values ranging from~1-20 g m −2 . The example images from sequences A and B suggest that most of the missing pixels are close to the emission source. Nevertheless, the unlisted pixels of each sequence are not only limited to the denser part of the plume. If we look at those pixels in more detail, 70% of the missing pixels of sequence A, 92% from the one of sequence B, and 29% from the unlisted ones of sequence C, are in the denser part of the plume (i.e., I O 3 indexes > 50).  Even though Mount Etna was in a quiescent stage of activity (i.e., not in eruption) and the plume was transported in the same direction, the training dataset represents only a part of the variety of the SO2 spectral patterns of the Etna NEC plume. The library created allows the massive retrieval of three sequences with only 15% to ~22% of unlisted plume pixels patterns, distributed in the whole sequence. One image example from each of the three datasets is presented in Figure 7 with the missing pixels colored in blue. As can be seen in Figure 7a, the sequence A has higher values of mass of SO2 per surface unit in the range of ~2.5-35 g m − ² in comparison with sequence B and C with values ranging from ~1-20 g m − ². The example images from sequences A and B suggest that most of the missing pixels are close to the emission source. Nevertheless, the unlisted pixels of each sequence are not only limited to the denser part of the plume. If we look at those pixels in more detail, 70% of the missing pixels of sequence A, 92% from the one of sequence B, and 29% from the unlisted ones of sequence C, are in the denser part of the plume (i.e., indexes >50).

SO 2 Emission Flux
The three tested datasets have more than 80% of their plume pixels included in the library, which is enough to estimate the SO 2 emission flux.
Over the few minutes of measurement, the vertical contribution of the plume transport speed was null or negligible with cross correlation coefficients of determination from 0.83 to 0.89. The vertical cross correlations indicate that the Etna plume propagated horizontally with a plume transport speed going from~5.8 m s −1 (sequence A) to~6.6 m s −1 (sequences B and C). Mount Etna was in a passive degassing period and the sequences recorded lasted a few minutes. Based on the cross-correlation calculation results, the plume transport speed remained uniform and was considered constant.
The boxes are fixed with a horizontal length corresponding to 50 pixels (~140 m) and a vertical length corresponding to the 64 vertical pixels (~179 m) of each of the three sequences and are placed close to the ground part of the image. The time evolution of the SO 2 emission flux of the three sequences is presented in Figure 8 and the average results from the three datasets are presented in Table 5. The position and the small size of the boxes,~25.2 km 2 , allow us to capture the heterogeneity of the concentration of SO 2 coming out of the crater. The variability of the SO 2 emission flux is important, with values going from single to double within seconds, especially for the sequences A and B. Short videos showing the dynamism of the plume for the three sequences are available in the additional material of this paper (externally hosted Supplementary Materials Files 1-3).
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 21 The three tested datasets have more than 80% of their plume pixels included in the library, which is enough to estimate the SO2 emission flux.
Over the few minutes of measurement, the vertical contribution of the plume transport speed was null or negligible with cross correlation coefficients of determination from 0.83 to 0.89. The vertical cross correlations indicate that the Etna plume propagated horizontally with a plume transport speed going from ~5.8 m s −1 (sequence A) to ~6.6 m s −1 (sequences B and C). Mount Etna was in a passive degassing period and the sequences recorded lasted a few minutes. Based on the cross-correlation calculation results, the plume transport speed remained uniform and was considered constant.
The boxes are fixed with a horizontal length corresponding to 50 pixels (~140 m) and a vertical length corresponding to the 64 vertical pixels (~179 m) of each of the three sequences and are placed close to the ground part of the image. The time evolution of the SO2 emission flux of the three sequences is presented in Figure 8 and the average results from the three datasets are presented in Table 5. The position and the small size of the boxes, ~25.2 km 2 , allow us to capture the heterogeneity of the concentration of SO2 coming out of the crater. The variability of the SO2 emission flux is important, with values going from single to double within seconds, especially for the sequences A and B. Short videos showing the dynamism of the plume for the three sequences are available in the additional material of this paper (externally hosted Supplementary Materials files 1-3).

Discussion
The number of patterns of the training dataset is small in comparison to the diversity of Mount Etna behaviors [2]. Nevertheless, we managed to retrieve the SO 2 SCD of three sequences from the IMAGETNA campaign.  [54] obtained SO 2 flux from 270 to 1000 t day −1 with a standard error of 30% using a ground-based lightweight open-path Fourier transform infrared spectrometry in active mode with a portable infrared lamp. Gliß et al. (2018) [55] presented an average estimation, from 15 min of ground-based UV camera measurements on 16 September 2015, of the NEC SO 2 emission flux of~8 kg s −1 (equivalent to~690 t day −1 ). In Prata and Bernardo (2014), the obtained SO 2 emission flux is in the range of 10-20 kg s −1 (equivalent to~860-1730 t day −1 ) along a period of 7 h of measurement in September 2003 with data acquired every 4-6 min using the multifilter thermal infrared camera system Cyclops [27]. Our results are in the range of values obtained in the literature. Nevertheless, our estimation of the North East Crater SO 2 emission flux are only partial estimations since the plume of Mount Etna was trimmed by the field of view of the Hyper-Cam imager, capturing only a part of the total plume of the NEC, which, moreover, would be the gas emission of only one of the four main craters of Mount Etna. Oppenheimer et al. (2006) presented SO 2 column amounts from ground-based UV spectrometer measurements with a scattered skylight as a UV source, ranging from 0.5 to 3 g m −2 for measurements within La Voragine summit crater, which was the strongest gas sources at the moment of the measurements, and values up to 2 g m −2 for measurements of the plume of Etna from a location at 2 km from the craters [56]. The SO 2 emission flux also measured by UV spectroscopy on that day was~11 kg s −1 (equivalent to~950 t day −1 ). The SO 2 column amount from our three datasets are a bit higher with average values ranging from 4.5 to 10.5 g m −2 .
One of the main sources of uncertainties is related to the estimation of the plume transport speed. The error on the retrieval of SO 2 SCD also impact the estimation of the SO 2 emission flux. For the training dataset, the average error is~5.8%. The uncertainties of the SO 2 emission flux estimation are mostly related to the plume transport speed determination and the averaging of the SO 2 slant column densities calculated in each class outcoming from the training dataset classification. The error on the plume transport speed is directly linked to the estimation error on the distance from the instrument to the plume that was determined visually with the relief. We consider this error as being~10%, so that the average error on the estimate of the SO 2 emission flux is~16% for the presented datasets. It is in the range or even lower than the uncertainties on SO 2 emission flux of similar studies like the one from Gabrieli et al. on Kïlauea volcano in Hawaii with less than 25% of error mainly induced by the determination of wind velocities [57].

Conclusions
In this paper, we have presented a spectral classification methodology of the Hyper-Cam images pixels from the IMAGETNA campaign. The network architecture is composed of two layers corresponding to two features extracted from the spectral bands of O 3 and SO 2 . One feature characterizing the opacity of the plume and the second one characterizing the concentration in SO 2 . The accuracy of the classification results was evaluated in correlation with the results of the time-consuming pixel-by-pixel retrieval, giving a first order regression slope of 0.97 and a very good determination coefficient of 0.94 in the diluted part of the training dataset. This method allows a fast and accurate reconstruction of SO 2 slant column density images of three sequences from the IMAGETNA campaign. Each image is processed in less than 40 s instead of 7 days with a pixel-by-pixel retrieval methodology.
The plume transport speed was determined by calculating the cross correlation between time series of mean SO 2 slant column densities along two horizontal and two vertical cross sections. The SO 2 emission flux relative of the North East Crater from Mount Etna was estimated using the "box method" and the results ranged from~5 to 11 kg s −1 (~430-940 t day −1 ) with an error of~16%.
A way to improve the developed methodology would be to include neighborhood pixels information. The results presented here were performed on sequences of a few minutes with similar weather conditions, plume level of activity, and measurement characteristics. In the future, it would be interesting to test this image pixels classification method on datasets with other characteristics to enrich the training dataset and develop a fast procedure to enlarge the library, and, thus, retrieve the missing features. A wider field of view or a different geometrical configuration of measurement, making it possible to capture the entire Etna plume, would also allow comparisons to be made with other methods of spectra retrieval and estimation of the SO 2 emission flux.
The use of a hyperspectral infrared imager, with an acquisition time of images of a few s, associated with the classification and retrieval time of the spectra of an image lower than a min, would allow the monitoring of the evolution of SO 2 emission flux in near-real-time at a high frequency. All of this would contribute to help prevent any hazards coming from active volcanoes but also strengthen the attractiveness and interest of this technology for volcanic gas emissions monitoring.