AEROSOL PLUMES CHARACTERIZATION BY HYPERSPECTRAL IMAGES COUPLED WITH SENTINEL-2 PRODUCTS

In this paper, we focus on the retrieval of microphysical and optical properties of industrial aerosol plumes through a process using airborne hyperspectral and Sentinel-2 multi-spectral images. The process allows first to perform atmospheric correction and second to determine background aerosols thanks to a comparison between hyperspectral and Sentinel-2 reflectances. Hyperspectral methodologies use the radiance differential between the measurement in the plume and the corresponding measurements out of the plume to estimate plume properties. To retrieve the surface reflectance under the plume, a principal component analysis coupling hyperspectral and multispectral data class by class is achieved. The developed method aims to compare the difference between measured and estimated reflectance with a radiative transfer model accounting for plume properties (particle radius and aerosol optical thickness of the plume). We have applied the method to a steel plant in the south of France. The retrieved plume show an aerosol mean radius between 0.05 and 0.2μm with a mean aerosol optical thickness about 0.05 along the plume.


INTRODUCTION
Human activities are responsible for large emissions of particulate matter. Atmospheric particulate matter (also called aerosols) is responsible for air quality degradation and the exposure to ambient fine particulate matter (PM2.5) is a global health concern (World Health Organization, 2006; Burnett et al., 2018).
The industrial sector is the second largest source of primary PM10 emissions (Guerreiro et al., 2014) representing 20 to 25% of the global emissions in Europe between 2003 to 2012. The industrial inventory for particle matter emissions is updated each year according to specific guidelines (European Environment Agency, 2019; Andre et al., 2019). Emissions are reported by industries and depends on manufacturing process, fuel consumption, emission factors and abatement technology. As emission are reported yearly, the instantaneous emission fluxes are subject to large uncertainties. The monitoring of atmospheric concentrations in the vicinity of industrial sources can provide a constrain on the estimation of emission fluxes.
Ground-based remote sensing techniques are commonly used for fine particles monitoring (Amaral et al., 2015) and provide a precise measurement of the physical and optical particle properties. For large areas without a dense network of instruments, airborne and satellite remote sensing techniques may be considered as an alternative for ground-based systems. To this end , hyperspectral imagers can be used to analyse particulate matter. A large number of spectral channels at a high spatial resolution can provide information on the particle properties like size distribution on the one hand and absorbing and scattering ratio on the other (Alakian et al., 2008;Deschamps, 2012). Hyperspectral methodologies use the differential between the radiance or surface reflectance modified by the plume, respectively L plume sensor and ρ plume soil and the corresponding radiance or reflectance without the plume, Lsensor and ρ soil , to estimate plumes properties. Commonly L plume sensor is the measured hyperspectral data and Lsensor is estimated from pixels out of the plume (Philippets et al., 2018) or multi-temporal data (Foucher et al., 2019).
This paper focus on airborne hyperspectral images at 2m resolution over a steel factory in southern France. A three steps process has been applied to a selected image over the factory in 2016. First, we simulate a database of atmospheres containing different plume types. Second, we retrieve the surface reflectance under the plume ρ corrected soil thanks to a coupling between the Sentinel-2 reflectance product and the hyperspectral reflectance image. At the end, we retrieve plume properties thanks to a differential model. Differential model consists in comparing the simulated reflectance differential with the measured reflectance differential. The best matching between the simulation and the measurement is performed with the Cluster-Tuned Matched Filter (CTMF) criterion.

Airborne hyperspectral data
The study focuses on a single plume on a radiance W/(m 2 .sr.µm) image taken the 17 February 2016 over the ArcelorMittal steel plant. The Fos-sur-Mer ArcelorMittal site is the second biggest steel plant in France to product steel, coils and tubes. According to the IREP 1 database, this site is respectively the first and the second transmitter of PM10 and CO2 in France in 2018 with respectively emissions of 1.230 t/year and 7.460.000 t/year. The studied plume (see figure 1) comes from the stack sinter plant. According to Philippets et al. (2018), the main releasing materials belong to the scattering aerosol family defined by Quinn et al. (1995). The plume image was measured with the HYSPEX hyperspectral camera aboard the SAFIRE ATR-42 research aircraft. HYSPEX camera has a spectral range between 0.41 and 2.5 µm with by 160 spectral channels in the VNIR (410 -996 nm) and 262 in the SWIR (970 -2500 nm). HYSPEX spatial resolution is 1 meter in the VNIR and 2 meters in the SWIR. The georeferencing of the hyperspectral image in UTM zone 31N is performed with QGIS 2 by using a 50cm spatial resolution ortho-rectified picture of the scene from the IGN website 3 . A change in scale is realised from 1 meter resolution to 10 meters resolution on the geo-referenced hyperspectral radiance image to match the Sentinel-2 spatial resolution.

Sentinel-2 product
The method is based on the coupling with a Sentinel-2 reflectance product from THEAI website 4 . The Sentinel-2 image is a level-2A surface reflectance product produced by the MAJA algorithm Hagolle et al. (2017). The Sentinel-2 image was taken the 5 February 2016, 15 days before the hyperspectral measurement. This product is the closest picture in the time which is operable.

Atmospheric correction and spatial matching
An atmospheric correction is performed in order to obtain a surface reflectance image from the hyperspectral image. First, the method uses an atmospheric correction algorithm named COCHISE. COCHISE is an iterative algorithm allowing to apply an atmospheric correction to a radiance image in W/(m 2 .sr.µm) in order to retrieve a surface reflectance map. Second, the coupling with Sentinel-2 surface reflectance product allows to help to choose the "best" background aerosol to compute the atmospheric correction. At the end, a spatial matching between the HYSPEX reflectance map and the Sentinel-2 is then applied thanks to the GEFOLKI algorithm (Brigot et al., 2016;Plyer et al., 2015). GEFOLKI algorithm allows to realise a pixel matching through a flow optic calculation and a resampling steps if needed between two images. At the end, we have a hyperspectral picture which spatially matches pixel-by-pixel with the Sentinel-2 picture (see figure 1).

METHODS
The retrieval method is divided into three parts : (i) the definition of the direct model, (ii) the estimation of ρ corrected soil and (iii) the plume characterisation by the differential method.

Direct model : Plume simulations
The first step of the process is the simulation of the radiative terms of an atmospheric column with the presence of a plume in the lower part. We compare this simulation with the measurements in the plume (see section 3.3)). The simulated atmospheres described by different radiative terms are: the atmospheric radiance Latm, the total solar irradiance E surf , the atmospheric transmittance Tatm and the atmospheric spherical albedo Satm. All these radiative terms are computed with COMANCHE (Poutier et al., 2002), a radiative transfer model which is a frontend of MODTRAN. The computation of the radiative terms depends on the flight altitude, the sensor spatial and spectral specifications and the background aerosols.
The method used a Look-up Table (LUT) computed with COCHISE algorithm of reflectance images computed with the same hyperspectral radiance image resulting for different atmospheric corrections. Reflectance images are computed for visibilities ranging from 10 to 100km −1 and several background aerosol types. Visibility is a MODTRAN parameter associated to the aerosol concentration in the atmosphere.
Background aerosols used are defined in MODTRAN as rural, urban and marine aerosols. The best background aerosol is obtained thanks to a comparison between the mean spectra out of the plume in surface reflectance for each class for the hyperspectral and the Sentinel-2 images. The hyperspectral surface reflectance map associated to the best background aerosol is chosen as a first guess for the reflectance estimation under the plume (see section 3.2). In the hyperspectral image, the best background aerosol is a marine aerosol (as defined in MODTRAN) with a 18km −1 visibility. For rural aerosols the surface reflectances computed are too low, contrary to the urban aerosols which generate a signal too strong.
Aerosol optical properties of the plume are simulated using Mie theory for scattering and absorbing aerosols. The plume is represented as a local layer with an 100 meters thickness located 10 meters above the ground. The aerosol optical thickness (AOT) of the plume, τ 550 ref , is defined at 0.1 at 550nm. Philippets et al. (2018) and Alakian et al. (2008) have considered the plume as an infinite layer. In this case all the radiative terms of the radiative equations are modified. When the plume is considered as punctual, the direct and diffuse solar irradiance are not modified. The figure 2 shows the plume signature under different hypothesis: an infinite plume, a semipunctual and a punctual plume. In this paper we use the punctual plume hypothesis. Figure 2. Signature of plume scattering residual aerosols according to several plume modellings compared to the plume measured signature over water. The hypotheses 1, 2 and 3 correspond respectively to an infinite plume (yellow), a semi-punctual plume (red) and a punctual plume(green).

Reflectance Estimation
In order to apply the differential model described in the section 3.3, the surface reflectance under the plume, ρ corrected soil is required. The existing methods used by Alakian et al. (2009) or Bojinski et al. (2002) suppose that the plume is transparent beyond 1.5 µm. For each pixel it is possible to compare ρ plume soil beyond 1.5 µm with ρ corrected soil in the same class. The nearest value is assigned to the pixel. This method was adapted in the VNIR (0.4 to 1.0 µm) by Deschamps (2012) for the CASI (Compact Airborn Spectrography Imager) images. Philippets et al. (2018) have adapted once again the Deschamps' approach by assigning to each plume pixel the mean spectrum of the class with the highest correlation.
For our study, we decide to associate the hyperspectral data with the Sentinel-2 reflectance product to compute ρ corrected soil . The use of the Sentinel-2 product allows first to already have a classification of the ground under the plume. This method is based on the description of the intraclass variability by eigenvectors and their associated weights.
The Sentinel-2 reflectance image is classified with a Random-Forest method according to user-defined classes: water, sparse vegetation, dense vegetation , concrete soil, dark soil and bright soil (see figure 3).
Second, Sentinel-2 reflectance allows to compute the weights associated to hyperspectral eigenvectors for each pixel to rebuild the hyperspectral signal.
We assume that the intraclass variation can be described by a few principal components. For a given N-dimensional spectrum dataset Si(λ), it is possible to define the mean spectrum by class SA as follow.
The Principal Components Analysis (PCA) is applied on a standardised spectrum defined as S * i = Si -SA. The PCA returns a set of eigenvectors Vj(λ) describing the spectrum variability.
In fact, this spectrum variability can be expressed as a weighted sum of the eigenvectors: where M is the number of eigenvectors and αj is the eigenvector associated weights.
This method is linear and can be applied to pixels containing "pure" or "mixed" surface types. Next, the aim is to choose a vector with weights minimizing the distance between the measured and the estimated spectra. The optimum weight vector is obtained by the following equation: where D is the matrix of the M eigenvectors sampled on the Sentinel-2 wavelengths, FA is the mean spectrum SA sampled at the same wavelengths and F is the Sentinel-2 spectrum of the studied pixel. The optimum weights vector formulation allows us to suppose that the Sentinel-2 intraclass variability is the same as the hyperspectral intraclass variability.
For each pixel in each class, a weight vector is computed and the hyperspectral spectrum can be reconstructed as : The use of 4 eigenvectors is sufficient to describe from 90 to 99% of the class variance. The PCA is applied on a dataset without the plume and after removing the first and the last percentiles of spectra sorted by norms. This removal criteria is applied empirically in order to delete spectra which were The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 4. Estimated reflectance spectrum (blue) compare to the hyperspectral reflectance spectrum (orange) and the Sentinel-2 reflectance spectrum (green). The estimated reflectance, the measurements and the Sentinel-2 spectra are respectively associated to the left, the middle and the right panels.
not properly classified. The figure 4 shows an example of a below plume estimated reflectance and the hyperspectral image corrected by the plume impact at 550nm. We are able to delete the impact of the plume on the reflectance image. The figure 5 shows the approach to estimate ρ corrected soil out of the plume to compute the estimated reflectance error δ∆Rmes. (2018) have used the differential radiance or surface reflectance defined as the products of the aerosol optical depth (AOT) and the aerosol signature to retrieve the plume properties, In this study we use the differential surface reflectance described by the radiative transfer equation. It consists in comparing ρ soil and ρ plume .

Retrieval of the plume aerosol properties : Differential model
Assuming a flat, homogeneous and Lambertian surface and for a monochromatic radiance acquired by a hyperspectral sensor, the surface reflectance can be expressed by the following equation in the reflective domain Chandrasekhar (1960): where L sensor = at-sensor radiance in (W/(m 2 .sr.µm)) L atm = atmospheric radiance, without interactions with the ground E surf = total solar irradiance T atm = total atmospheric transmittance All these radiative terms are computed with COMANCHE (Poutier et al., 2002) algorithm.
In the presence of a plume, all of these terms are modified. Assuming an infinite plume layer, the surface reflectance equation can be written as follow : with ∆L sensor plume , ∆L atm plume , ∆E surf plume , ∆T atm plume are respectively the at-sensor radiance differential, the upwelling atmospheric radiance differential, the downwelling irradiance differential and the upward transmittance differential. ∆ρ plume soil is the differential signal containing the desired aerosol signature.
The aerosols signature can be expressed as: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 6. Surface reflectance theoretical error for a 0.15µm modal radius plume for different types of ground.
The reflectance differential is computed for measured spectra inside the plume and plume simulations. The best simulated aerosol signature can be retrieved thanks to the higher Cluster-Tuned Mached Filter (CTMF) score of the simulation with the measurements. The CTMF was developed and used by Funk et al. (2001) to retrieve gas thermal signatures (CO2, SO2, N2O and O3) in the thermal infrared. The CTMF was adapted by Thorpe et al. (2012Thorpe et al. ( , 2016 and Dennison et al. (2013) for other gases in the reflective domain. Most recently, Philippets et al. (2018) used this method on aerosols.As opposite to gases that are characterised by well-defined absorption bands, aerosols have an impact on the entire reflective signal (reflective and infrared domains). The objective of this section is to determine the best correlation between the hyperspectral reflectance signal with the aerosol signature ∆ρ plume soil .
For a given hyperspectral image of N pixels with m spectral bands, the CTMF model can be described as the combination of ρ soil with the spectral reflectance signature of the plume aerosol ∆ρ plume soil . Both of these parameters are m dimensional vectors. After a classification of the image, the mean reflectance spectrum ρ soil j and the inverse of the correlation matrix C −1 j are computed for each class.For each class, an optimal filter qj can be computed as : With the optimal filter qj, it is possible to obtain the detection score fi,j for each pixel of the hyperspectral image : This approach by classes allows to not be impacted by the modification of the aerosol magnitude signature according to the ground type. The figure 6 shows the different aerosol signature magnitudes for a scattering plume with an AOT of 0.1 over different ground types.
Additionally, with a computation of the Root Mean Square Error (RMSE) of the reflectance differentials, it is possible to characterise at which wavelengths the plume signature will be detected (see figure 5)). For water, dark soils, sparse vegetation and dense vegetation, the aerosol signal can be studied on the whole wavelength band (0.4 to 0.9 µm). For the bright soil class, the reflectance surface estimation error δ∆Rerr is too high compared to the aerosol signals. Finally, for the concrete soil class, just the spectral domain from 0.4 to 0.7µm can be used.

CTMF application for plume properties retrieval
The retrieval process is applied to the HYSPEX hyperspectral image over Fos-sur-Mer and provides a set of CTMF maps. Each map represents the correlation of the measurement with an aerosol signal computed thanks to the atmospheric database.
On the different CTMF maps, the plume is clearly visible for aerosol signatures with a modal radius from 0.05 and 0.2 µm. For higher modal radius the plume is undetectable. The figure 7 shows the CTMF map corresponding to a 0.1µm aerosol radius. We can easily detect the plume above water due to the low surface reflectance. Bright soils data are flagged (see figure 5) due to an higher error on the surface reflectance estimation compared to the aerosol signal. A part of the plume is visible above sparse dense vegetation classes. Over concrete soils, its footprint detection is more complex. The major part of pixel signal is lower than δ∆Rmes.
By selecting the best CTMF score for each pixel, we can build a modal-radius aerosol map. This map represents the contribution of the major mode inside the plume. For the Fos-sur-Mer image, the plume major aerosol mode is a 0.1 to 0.15µm modal radius scattering aerosol.

The AOT map estimation
We retrieve the AOT by computing the optical thickness ratio. The optical thickness ratio is obtained by comparing the reflectance differential of the measurement with the reflectance differential of the simulation. The optical thickness ratio can be expressed by the following expression : The figure 8 shows the AOT map associated to the aerosol radius corresponding to the best CTMF score. For the studied plume, the total plume AOT varies between 0.02 to 0.08.

Conclusion and perspectives
We have developed a method based on hyperspectral measurements coupled with a multi-spectral data (Sentinel-2) at 10 meters spatial resolution. The proposed method allows first to perform the image atmospheric correction thanks to the retrieving of background aerosols with a comparison between hyperspectral and Sentinel-2 reflectance. Second, the rebuilding of the hyperspectral signal under the plume ρ corrected soil is possible by computing the hyperspectral eigenvectors and their associated weights obtained with the Sentinel-2 reflectance.
We compute first for each pixel a measured differential spectrum defined as the difference between the measured surface reflectance and the surface reflectance estimation. This differential represents the measured plume signature. We compare this signature to a differential model to retrieve plume parameters.
To retrieve the best aerosol modal radius and to estimate the AOT, we use the CTMF algorithm with the differential method. The CTMF compares the simulated aerosol plume signature with the measured plume signature from the reflectance image.
However, the Sentinel-2 product seems to underestimate the surface reflectance, especially over dark surface (ex : water). Using another dataset with Sentinel-2 products to retrieve the background aerosols should be a possibility in order to improve the atmospheric correction step. More, for classes with a high spectral variability as the "concrete soils" class, the reflectance reconstruction method is hard to apply. The reconstruction errors are too large compare to the aerosol plume signature.
As a future work, in order to improve the plume characterisation, we should develop first a multi-modal retrieval method for the radius. Second, a 3-dimensional representation of the plume is planned in order to better model the radiative impact of a plume according to the studied pixel location. At the end, we intend to use this process to compute a priori of radiative parameters for the use of an optimal estimation approach in order to characterise plumes.