DETECTION AND QUANTIFICATION OF INDUSTRIAL METHANE PLUME WITH THE AIRBORNE HYSPEX-NEO CAMERA AND APPLICATIONS TO SATELLITE DATA

The control and monitoring of greenhouse gases is an important issue for the study of climate change, for the impact in terms of public health or for the risks related to industry. An algorithm has been developped dedicated to commercial airborne hyperspectral camera as HySpex-NEO for the detection of industrial methane plume and the quantification of the emission source. HySpex-NEO is an imager used in airborne campaign with a spatial resolution of 1.4 m at flight altitude of 2 km, a swath of 650 m and a spectral resolution of 6 nm. This algorithm has been validated over a controlled release less than 100 g/s during an airborne campaign over Lacq (France) industry. It has also been applied to the Aliso Canyon leakage data acquired with AVIRIS JPL (Airborne Visible Infrared Spectrometer) with a spatial resolution of 6.8 m at flight altitude of 6 km, a swath of 5.6 km and a spectral resolution of 10 nm. Application to satellite hyperspectral data is shown on artificial data derived from airborne hyperspectral acquisitions.


INTRODUCTION
The control and monitoring of greenhouse gases is a scientific, societal, public health and environmental issue. The reduction of anthropogenic gas emissions federates countries around the Paris Agreement (2015). Methane (CH4) influences climate as the second most important anthropogenic greenhouse gas and air quality factor (Forster et al., 2007). Carbon dioxide (CO2) has a higher concentration than CH4 in the atmosphere but CH4 has 21 times the radiative forcing power of CO2 (Lelieveld et al., 1998, Griggs, Noguer, 2002. The atmospheric residence time of CH4 about 7.9 years (IPCC et al., 2007) is shorter than that of CO2 (one hundred years), which explains the collective efforts to reduce CH4 emissions.
While anthropogenic sources of CH4 accounted for 4-34% in pre-industrial era (Houweling et al., 2000), they amounted to about 60-70% in 1998. Global atmospheric CH4 has increased significantly with about 0.65 ppm in the pre-industrial epoch for about 1.8 ppm today (Etheridge et al., 1998, Dlugokencky et al., 2009). In addition, sources and sinks of CH4 emissions are scattered randomly on Earth and show a high degree of interannual variability (Bousquet et al., 2006). If we focus on anthropogenic emissions, the main sources of CH4 emissions are energy sector (fossil fuel production, ...), agriculture (domestic ruminants, rice cultivation, ...), industry and waste treatment (Kirschke et al., 2013).
One of the powerful tools for the detection and quantification of gas emissions is hyperspectral remote sensing. It is based on the fact that the spectral signal received by a hyperspectral sensor is affected by the potential presence of a plume in the line of sight. There are two kinds of hyperspectral instruments. The first kind is satellite sounders as IASI-NG (Crevoisier et al., 2014) or TROPOMI (Veefkind et al., 2012) with a spatial resolution greater than 7 x 7 kilometers. This resolution allows studying very large scale plumes such as the volcanic eruptions or the evolution of background concentrations. However, it is not adapted to industrial plumes with subkilometric spread. To overcome this spatial resolution matter, hyperspectral imagers in the short-wavelength infrared (SWIR) domain are used either in airborne campaigns or more recently onboard satellites. The flight altitude of airborne imagers leads to finer spatial resolutions than satellite instruments. For altitudes under 2 km, the ground resolution is metric (Bradley et al., 2011. But such campaigns are expensive to implement and cover a small area. The new hyperspectral satellites would make data more accessible and remove the locks related to cost and spatial coverage. However, due to a larger ground pixel size (30 m for PRISMA (Labate et al., 2009) and EnMap (Guanter et al., 2015)) and higher atmospheric absorption, the ability of satellites to detect and quantify CH4 plumes has yet to be demonstrated. The section 2 shows the sets of data available. The section 3 explains the three steps required to obtain a source flow rate estimation. The section 5 is dedicated to the algorithm results in the case of available satellite data.

DATA SITES
The results presented in this paper come from two airborne data sets. The first have been acquired by the HySpex camera, during a campaign above Lacq industries (France) in June 2017. The hyperspectral camera was a SWIR imager with a spatial resolution of 1.4 m at flight altitude of 2 km with a swath of 650 m and a spectral resolution of 6 nm. The spectral range of acquisitions extends from 967 nm to 2501 nm with 256 bands. Of all the existing chimneys on the site, one platform was designed as a testing area to release gas in controlled way to simulate accidental gas leaks. In such a way that for each release, the in situ flow rate is recorded.
A second set of data were acquired by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) which provides 224 spectral bands with a spectral resolution of 10 nm. To reduce the size of the initial hyperspectral data, a selection of the 160 spectral bands between 966 and 2495 nm is applied. A spatial resolution of 6.8 m are obtained with the Twin Otter aircraft from the Aliso Canyon (USA). Aircraft flew over a CH4 superemitter during the explosion of natural gas wells in 2015 (Thompson et al., 2016).
The initial size of image is 320 x 2607 pixels for HySpex and 829 x 3966 pixels for AVIRIS. During the development of the algorithm, only a section around the chimney is taken, reducing the size to 170 x 285 pixels for HySpex and to 650 x 1000 pixels for AVIRIS. A representation of these two data is plotted in the figure 1.
A B Figure 1. A : Section of the HySpex data above the chimney of the Lacq industries. B : Section of the AVIRIS data. The two images are displayed at 2.3 µm.

METHODS
The principle of the algorithm is based on estimation the transmission of the observed plume to estimate its integrated concentration pixel by pixel. Therefore, a first detection step allows to separate plume-present pixels to plume-absent pixels, using Cluster-Tuned Matched Filter (CTMF) algorithm. Then, the quantification method developed in this paper can be applied to identified plume-present pixels.

Direct model of measured radiance
The total radiance measured by the instrument can be expressed as : where Ltot = total measured radiance (in W/m 2 /sr/µm) R = ground reflectance E0 = solar irradiance (in W/m 2 /µm) τ dir = direct transmission τ dif f = effective contribution of scattering effects S = spherical albedo Latm = atmospheric radiance (in W/m 2 /sr/µm) A sensitivity study of the impact of each term of equation 1 with COMANCHE radiative code (Poutier et al., 2002) leads to several simplifications. The (1 − RS) term tends to 1 on the range of reflectance and spherical albedo. At 2.3 µm, the scattering phenomena are negligible compared to the absorption. Therefore, the atmospheric radiance and the scattering transmission can be deleted of equation (1). The direct transmission corrresponds to the atmospheric transmission (τatm) multiplied by the CH4 transmission (τ ch 4 ). The CH4 transmission is composed of the downward transmission and the upward transmission (affected by the zenith angle θ of the sun). These transmissions can be expressed as the total CH4 transmission : To obtain the plume spectral signature from the total radiance, we use : where L pl is the measured radiance and L nopl is the reference radiance (no plume radiance).

Preprocessing step
A preprocessing is needed to define the wavelengths that will be used at different steps of the processing of hyperspectral data. This step consists of obtaining the transmission of the atmospheric gas. We consider that the transmission of all the gases present in the atmosphere can be approximated by the water transmission due to its high relative concentration. Using the monochromatic gas absorption (Sharpe et al., 2004), the target gas transmission is computed with the following relation : where λc = band center wavelength (in nm) τgas = transmission of the gas ρ = concentration in particle per million per meter (in ppm.m) Agas(λ) = monochromatic unitary absorption (λ in nm) S(λ) = spectral sensitivity of the instrument This equation leads to the gas transmission at the wavelengths of the chosen instrument by using typical value of water and CH4 concentration (respectively 10 7 ppm.m and 10 5 ppm.m). The algorithm calculates the transmission of water and CH4 by taking into account the spectral sensitivities of AVIRIS or HySpex. The SWIR spectrum is affected by the large water absorption bands around 1.4 and 1.9 µm. CH4 has also absorption bands at 1.65 and 2.3 µm.

Classification
The classification is applied by using the k-means function of Python 3.6 to cluster each hyperspectral pixel into a class number N (typically 15). The algorithm runs a first classification with N classes. All classes with few pixels are deleted and a second classification is running without taking into account these isolated pixels. The two previous classifications are running without the need to initialize class centers and with 300 iterations. An example of the classification result for HySpex data is shown on the figure 2.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition)

Detection
The detection step of the algorithm is based on the CTMF algorithm (Funk et al., 2001). Its aims to look for CH4 transmission spectral features in the hyperspectral data class by class. The CTMF algorithm and previous versions (Thorpe et al., 2014, Thompson et al., 2015, Hulley et al., 2016) take into consideration the theoretical transmission of a target gas to determine the level of similarity of each pixel to this theoretical spectrum. A high CTMF score reflects that this pixel is impacted by the target gas. This method assumes that the measured radiance can be expressed as a linear combination of the target gas signature and the signal of the rest of the scene as : where L = at sensor radiance b = spectral signature of the target gas α = strength of the plume gas L bkg = radiance of the scene without plume If we apply the CTMF method, the L of the equation (7) corresponds to the left term of the equation (5) and the b of the equation (7) is equal to τ ch 4 − 1. In addition, we define the optimum filter vector q of the equation (8) and we apply it on the hyperspectral data following the equation (9) : In order to obtain the CTMF score, we applied the q vector to the hyperspectral image using the equation (9). The CTMF score is shown in the figure (3) after applying a mask by threshold value.

Transmission estimation
The plume transmission is directly linked to the concentration as in the equation (2). But, the CH4 transmission requires to determine the reference radiance L nopl as shown by the equation (5). The goal of this step is to retrieve the spectrum of each pixel if the gas plume was not present. Using morphological transformations, a mask with the pixels that surround the plume is created. For each pixel of the plume, the reference is chosen among pixels of surrounding mask. To chose this specific pixel, 3 criteria are used : i) belonging to the surrounding mask, ii) belonging to the same class than the plume pixel, iii) having a spectrum similar to plume pixel at the wavelength not impacted by CH4 (by minimisation of the root-mean square error RMSE). The figure 4 represents an example of the chosen spectrum to rebuild the reference image. The blue and red curves are similar outside the CH4 range and we can observe small variations between 2.25 and 2.4 µm. Reference spectra become available for all detected pixels.  Reference spectra without plume are used to provide the observed transmission following the equation (5). Then, the comparison between this transmission and theoretical transmissions from a Look-Up Table (LUT) gives access to the concentration of the plume.

Detection and concentration map
A LUT of CH4 transmission has been created (from equation (6)) from 0 to 500 000 ppm.m for each instrument. For each plume pixel, the selected concentration corresponds to the minimum RMSE between the LUT transmissions and the measured transmission as given in the equation (5).
The figure 5 shows the concentration map of the HySpex (5A) and AVIRIS (5B) data. The emission point is identifiable in the both images and the algorithm detects in one the hand the main plume and on the other hand some areas around the plume that may be false alarms or gas accumulation locations.

Estimation of the flow rate
Determining the flow rate of a source is important to quantify industrial sources or natural leaks. In this paper, we chose to determine the flow rate using the central part of the plume, called "plume center" in opposition to the head (at the chimney localisation) and the plume tail (where the plume is largely dispersed). The figure 6A shows the entire concentration map of the HySpex image. The figure 6B represents a zoom on the plume after a rotation of 45 • . The next step is to reduce the plume to a line by summing the concentrations belong to the width of the plume perpendicular of the flow direction (a sum by columns in this case). The mass per unit length is plotted in the figure 6C.
The flow rate of the chimney can be obtained from the figure 6C by taking into account the mean of the mass per unit length. In this HySpex case, the mean of the mass per unit length is equal to 51 g/m.
To obtain an estimation of the flow rate, the last step is to multiply the mean mass per unit length by the wind speed. During the Lacq campaign, several sensors measured the speed of the wind at different altitudes. A lidar system ZEPHIR has recorded the wind with a frequency around 15 s at 11 m altitude (and other levels) and a sonic 3D instrument called METEK has measured the wind speed at 10 m of altitude every minute. During the Lacq campaign, the wind was very variable in strength and in direction as shown in the figure 7 for the wind speed. As we can see in the figure 7B, in only 2 minutes, the speed of the wind varies from 0.78 to 5.09 m/s for ZEPHIR data at 11 m of altitude and from 1.52 to 3.98 m/s for the METEK data at 10 m. With the mean of the mass per unit length of 51 g/m and a wind  speed ranging from 1 and 4 m/s, the algorithm estimates a flow rate between 51 and 204 g/s for the HySpex image.
The in situ measurement of the flow rate is to 75 g/s. The last value lies in the estimation range provided by our algorithm. Several points can explain the large range of the flow rate derived by the algorithm. The first point is related to the acquisition of measurements. In this respect, we have to note that the hyperspectral image acquisition is not a snapshot and during this acquisition flight, the wind can change significantly. Due to the temporal standard deviation, the measurement of the wind speed brings a lot of uncertainties to the flow rate estimation. In addition, the wind speed is measured at 10 m by the instruments while the emission source is close to the ground. This point can introduce a overestimation of flow rate because the wind is stronger as the altitude is increasing. However, even if the flow rate seems to be overestimated, the in situ value is in the uncertainty range. To conclude the airborne data part, we have developed an algorithm which estimate the flow rate of a source from SWIR hyperspectral data combining with wind data. The obtained results are in good agreement with in situ measurement despite the large uncertainty of the wind speed data.

APPLICATION TO SATELLITE LEVEL
The final purpose of this article is to demonstrate the ability of our algorithm to detect and quantify plumes with satellite data. The atmospheric layer above the airborne sensor is significant and the ground pixel size is larger from satellite than from airborne images. The following section of this paper will be dedicated to i) add the atmospheric layer above the airborne and ii) downgrade the spatial resolution. This part focuses on the AVIRIS data.

Top Of Atmosphere level
The initial data are recorded during airborne campaigns with 6 km flight altitudes by the AVIRIS imager. The satellite data will be acquired at more than 850 km. We have to take into consideration that the signal will pass through a significant additional atmospheric path. As in the previous part, we assume that the corresponding layer is only composed to water vapor. Using the radiative transfer code COMANCHE (Poutier et al., 2002), we run simulations to determine the direct upward transmission of the atmosphere between flight altitude campaign and the top of atmosphere (TOA). The initial radiance is then affected by this term following : where L (ima,T OA) = radiance in the TOA level (in W/m 2 /sr/µm) L (ima,f ly) = radiance of airborne campaign (in W/m 2 /sr/µm) τ (f ly−T OA) = direct transmission between the airborne and the top of atmosphere (TOA) As mentioned on the section 3.2.3, the scattering terms (in the SWIR) are negligible compared to the absorption terms. Then, only the direct transmission between the airborne and the TOA is considered. The algorithm is now applied to the TOA level radiance simulated from the airborne images and the results are shown in the figure 8. Only few differences are visible between the concentration map of the plume with the data from satellite altitude (figure 8A) and with airborne data ( figure 5B). This consistency is visible on the lineic masses shown in figure 8B showing the mass per unit length for airborne initial data (blue curve) and for the data translated to the TOA (dashed red curve). These two curves are similar and the differences can be due to the detection step : the mask of present/absent plume pixel can generate false alarms in the edge of plume if the gas occupies a small part of the pixel.

Spatial convolution
The spatial resolution of new satellite hyperspectral imager as PRISMA (Labate et al., 2009) or EnMap (Guanter et al., 2015) is of 30 m. As a reminder, the spatial resolution of AVIRIS image is of 6.8 m. A first order reduction of the spatial resolution of images can be done by averaging a pixels square of the initial matrix. This convolution does not take into consideration the noise introduced by the averaging process or other perturbations which can be present in future satellite data.
We applied this convolution on the initial (airborne) images and the algorithm is then run. The figure 9A presents the concentration map of the plume for a spatial resolution up to 27.2 m. The global shape of the plume is identifiable and stays similar to the 6.8 m spatial resolution image ( figure 8A). However, the mean mass per unit length presents more differences than the previous case. An important point is that the global shape of the curve at 27.2 m of spatial resolution is similar to the finer spatial resolution. As such, the plume detection and the estimation of its concentration seems still possible with satellite data. But, we can note that between 1000 and 1500 m, an overestimation of the concentration occurs where the plume has a larger width than the other part (the tail plume). This may be due to averaging of pixels square at the initial step that emphasizes the spread of the plume. This will be investigated in the future.
The satellite data allow to detect and to quantify the plume in the case of the AVIRIS campaign plume with kilometric extension. The uncertainties of (simulated) satellite results are larger as shown by the figure 8B. The present work has to be continued to study the limits of the algorithm to detect and to quantify small plume and low intensity plumes from satellite data. Mass per unit length (in g/m) along the plume (in m) for different initial images.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition)

CONCLUSION
In this paper, we demonstrate the ability of our algorithm to detect CH4 plume and to quantify the flow rate of the source, with an uncertainty mainly due to the wind information. This characterization is carried out from the commercial HySpex camera using CTMF method in the detection step and using the plume transmission computation and LUT comparison to quantify the source. The algorithm detects a 60 m wide plume and a flow rate between 50 and 200 g/s while the in situ measurement gives 75 g/s. The second part of this paper is dedicated to the algorithm results for simulated satellite images. For the both HySpex and AVIRIS, the translation to the TOA level has few impacts on the quantification of the plume. The study focuses on the absorption bands around 2.3 µm that leads to some advantage since scattering effects can be neglected. As the HySpex plume is 60 m wide, we focus to the AVIRIS plume for the degradation of the spatial resolution. From 6.8 m (airborne resolution) to 27.2 m (near PRISMA or EnMap resolution), the algorithm is able to detect the plume and to estimate its concentration but a 27.2 m spatial resolution leads to a overestimation of the plume mean mass per unit length in the tail part. In the future, several points have to be investigated. The algorithm can be upgraded by reducing false alarms. The detection and quantification sensitivities of the algorithm should be investigated. Of course they will depend on the spatial and spectral resolution. This study can be conducted using synthetic data. Testing and validating this algorithm with true satellite data will be our final goal.