COMPARISON OF ESTIMATION METHODS TO QUANTIFY METHANE PLUME CONCENTRATION AT HIGH SPATIAL RESOLUTION FROM HYPERSPECTRAL IMAGES

The detection and the quantification of greenhouse gases is essential to climate change studies, avoid leakages in industrial site, prevent accidental explosions. Because of these properties, methane (CH4) is an important target gas in remote sensing quantification. We focused on industrial plume quantification with data obtained during airborne campaign with HySpex-Neo hyperspectral camera. The 1.4 m of spatial resolution allows comparing quantification methods on real data combined with full or semi-synthetic plume case. A linear method largely used in the literature is compared with a quantification method based on the rebuilt background image and the estimation of plume transmission (PTE method). We have developed a hybrid approach using intermediate results of two previous methods. The hybrid method is based on the optimal estimation (OE) formalism and is providing uncertainty estimates. We show that the linear method underestimates the concentration of plumes for concentration above 5000 ppm.m. For low reflectance pixels, the hybrid method is more robust than the PTE method. The uncertainty of the hybrid method is about 30% for pixels with concentrations above 5000 ppm.m. For a HySpex-Neo image, the total mass of the plume is underestimated by 30% with the linear method compared to the hybrid method.


INTRODUCTION
Methane (CH 4 ) is the second most important anthropogenic greenhouse gas (GHG) due to its strong radiative forcing. While carbon dioxide has a higher concentration in the atmosphere, CH 4 has a global warming potential about 28 times larger than that of CO 2 (IPCC, 2007;Borchardt, 2020;Krings et al. 2013). Due to the residence time of CH 4 in the atmosphere of less than 10 years (IPCC, 2007) in comparison to several hundred years lifetime of CO 2 , countries aim to reduce GHG emissions and those of CH 4 in particular (Paris Agreement, 2015). This goal is linked to the exponential growth rate of methane concentration: from about 0.65 ppm in pre-industrial times to about 1.8 ppm today (Dlugokencky, 2009;Etheridge, 1998).
As all other species, CH 4 has sources and sinks that affect its atmospheric concentration. On one hand, CH 4 sinks are natural and rather constant over time: oxidation by the OH radical in the troposphere, absorption by soils and oxidation in the stratosphere (Wetch, 2014). On the other hand, the sources are spread between natural, as wetlands termites and oceans (Bousquet, 2006), and anthropogenic emissions. The latter ones followed the industrial revolution: whereas they accounted for 4-34 % of methane emissions in the pre-industrial epoch (Houweling, 2000), they amount to about 60-70 % in 1998. Anthropogenic methane comes from the energy sector (fossil fuel production...), agriculture, industry, and waste treatment (Kirschke, 2013). The study presented in this paper is focusing on anthropogenic CH 4 emissions. In addition to the overall environmental issue, CH 4 is of particular interest for safety. Methane has well known explosive properties and inhaling high concentrations of methane can lead to respiratory complications (Thorpe, 2020). A permanent methane leakage is a profit loss for industries. It then clearly appears that the detection and quantification of CH 4 emissions are essential today.
Corresponding author: nicolas.nesme@onera.fr Hyperspectral instruments provide valuable data to monitor CH 4 . Three kinds of data are now available. The first one is provided by satellite sounders such as TROPOMI (Veefkind, 2012) with a kilometric spatial resolution. However, industrial plumes have commonly smaller extension (except large accidents as Aliso Canyon (Thompson, 2016)) and cannot be quantified with kilometric resolution. To overcome this spatial resolution problem, hyperspectral imagers in the shortwavelength infrared (SWIR) spectral domain are used either in airborne campaigns or more recently on-board satellites. Indeed, PRISMA, a satellite imager with 30 m spatial resolution (Labate, 2009), provides hyperspectral data and others instrument will increase the amount of available satellite imagers data like EnMap (Guanter, 2015). Even if recent papers show detection and quantification with satellite data, we focus here on hyperspectral airborne images to compare different methods used to quantify CH 4 plumes.
A linear method commonly used in the literature is compared to the plume transmission estimation (PTE) method. The intermediate results of these two methods provide values used in a third method. In this article, a new method based on the optimal estimation has been developed and permits an a posteriori study of the results with uncertainty values. The data available for comparison is described in the section 2. Section 3 presents the three methods and different study cases are reported in section 4. A summary of the results obtained and concluding remarks appear in the section 5.

Direct model of radiance
The radiance measured by an imaging instrument without any plume in the line of sight can be expressed as:
In the SWIR domain, scattering effects can be neglected. Indeed, a sensibility study with a radiative transfer code allows quantifying the impact of each term and shows that a simplification of this equation under reasonable hypotheses is leading to the following expression (Nesme, 2020): Where = atmospheric transmission.

Plume observation
As shown by Nesme, 2020, when the line of sight is crossing the plume, the transmission of the gas (for single gas plume) is to be compounded with the atmospheric transmission: with: Where = plume concentration (in ppm.m) = monochromatic absorption (in ppm -1 .m -1 ) at the instrument wavelengths by considering its spectral sensibility ( Figure 1). = solar zenith angle Figure 1. Methane monochromatic absorption in ppm -1 .m -1 of the database (Sharpe, 2004) in blue line and the methane absorption in ppm -1 .m -1 convolved by the spectral sensitivity of the hyperspectral instrument in red line.
Using the information on the location and time of the simulation, on the spectral sensitivity of the instrument and on the instrument altitude, a radiative transfer code can estimate the atmospheric quantities of Equation 2. The radiance in the absence of plume case is while in the presence of plume the radiance is :

Data
In this paper, three kinds of data are used: full-synthetic data, semi-synthetic data and real case data from an airborne campaign. A semi-synthetic image is composed of a real background and a synthetic gas plume.

Full-synthetic cases:
The ECOSTRESS JPL Spectral Library (https://speclib.jpl.nasa.gov/) provides a wide range of reflectance spectra. Nine of them were selected because of their amplitude and variability: water, black paint on aluminium rod, wood, asphalt, brick, green paint on wood, 10-year-old copper metal, paving stones and shingles. These spectra are plotted in Figure 2 as a function of wavelength. A full-synthetic with a plume present in the image has been created by using the atmospheric quantities obtained with the radiative transfer code: on one side, Equation 2 allows calculating radiance from database reflectance spectrum and on the other side, Equations 3 and 5 allow introducing a synthetic plume on the radiance image. The atmospheric quantities used for simulating a scene in summer at 13:00 UTC are from in situ atmospheric profiles.

Figure 2.
Reflectance spectra for nine surface types chosen for their different spectral variabilities.

Real hyperspectral image case:
Two hyperspectral radiance images were acquired during an airborne campaign above Lacq industries (France) in June 2017 by ONERA in collaboration with TOTAL. The Hyspex hyperspectral camera from NEO was used to acquire SWIR images with 1.4 m spatial resolution at flight altitude of 2 km with a swath of 650 m. This commercial camera has a spectral resolution of 6 nm and covered the domain from 967 nm to 2501 nm with 256 channels. The initial images have about 320 x 2607 pixels and a section is presented in the Figure 3. A test-stack started emitting methane in a controlled way before the airborne flight above the stack to record plume images. The emission point of this scene are represented by the red circle in the Figure 3. For the two hyperspectral radiance images, a CH 4 plume is present in the data. Results will be focus on one of the two images.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition) Figure 3. Hyperspectral radiance image acquired with HySpex-NEO camera and plotted at 967 nm. The circle represents the emission point.

Semi-synthetic case:
An atmospheric correction was applied on a section without real plume of the two HySpex radiance images to obtain reflectance hyperspectral maps on the industrial site. The Figure 4 represents the synthetic CH 4 plume used on this paper (the same intensity and variation are used for full-synthetic cases). The background displayed in Figure 4 is the reflectance map at 1538 nm. The semi and full-synthetic cases have been created with a radiance affected by the plume in a unique pass through the plume, in the ascending path on the nadir direction. In this way, A function to add HySpex instrument noise is applied to the full and semi-synthetic data. Concerning the full-synthetic case, the brick reflectance map reaches a SNR about 150 for reflectance 54% at 2200 nm. For the water case, the SNR is about 10 for a reflectance of 2% at 2200 nm. In comparison, new satellite imagers as PRISMA specified signal noise ratio (SNR) about 180 at 2200 nm for a reflectance of 30%.

METHODS
A pre-processing of the data is applied. First a classification is running to cluster hyperspectral pixels outside the wavelengths impacted by CH 4 . A detection method is then applied to create a mask of the pixels affected by the CH 4 plume. Nesme, 2020 detail the pre-processing step. In this paper, the detection step is considered to be robust and to provide a correct identification of plume-affected pixels. This step allows focusing the comparison of quantification methods on a reduced number of pixels.

Linear method
The first method presented in this paper is a linear method (LM) commonly used in the literature. The LM is derived from the Clustering-Tuned Matched-Filter method (CTMF) that uses a probability score for the difference between observed and theoretical spectra (Funk, 2001). For remote sensing of CH 4 , the signal is the transmission of the gas. The CTMF algorithm (Thorpe, 2014;Hulley, 2016) takes into consideration the CH 4 theoretical transmission to determine the similarity between the spectrum of each pixel and the theoretical spectrum. The higher the CH 4 column, the higher the CTMF score. This method assumes that the measured radiance can be expressed as a linear combination of the CH 4 signature and the signal of the background. Dividing the CTMF score by the product of the theoretical signal and the optimal filter vector (Thompson, 2015) leads to the CH 4 concentration of the pixel. Nevertheless, LM assumes that the transmission through the plume can be linearized as: 4 = 1 − (for a single pass). This assumption is valid for low concentrations and leads to an underestimation for highly concentrated plumes.

Plume Transmission Estimation (PTE) method
This method is based on Equation 3 which relates the plume concentration to the transmission. If the CH 4 transmission is known, the concentration is obtained by matrix calculation. The Equations 4 and 5 allow calculating the observed CH 4 transmission as: The value of is unknown and must be estimated. The estimation method selects for each pixel under study a similar pixel in the image. To be selected, the pixel must i) be outside the plume, ii) belong to the same class as the studied pixel, iii) have a spectrum similar to plume pixel spectrum at the wavelength not impacted by CH 4 (by minimization of the rootmean square error RMSE). Once this pixel is identified, its spectrum is used to rebuild the radiance of the plume pixel for the wavelength impacted by CH 4 . The Equation 6 gives then the transmission of observed plume and a concentration is calculated for each plume-present pixel.

Hybrid method
The hybrid approach is based on the optimal estimation (OE) formalism often used for solving atmospheric retrieval problems (Rodgers, 2000). The OEM compares the measured spectrum with spectrum obtained by a direct forward model. Iterations of the direct model are run to decrease the difference between observation and model, through minimization of a cost function. The direct model is commonly based on a radiative transfer model (RTM) that is quite time-consuming for such iterative scheme. Here, the direct model is simply the product of the rebuilt background ( = ) and the CH 4 transmission, as derived from Equation 6. The hybrid concentration is given by:

Method comparison metrics
The concentration of the synthetic plume is known and each line of the plume has the same concentration. There are therefore 10 pixels with the same real value. The first metric is the difference between the mean concentrations found for this line according to the synthetic concentration. For each concentration, the standard deviation (stdv) of the retrieved concentration is then calculated over the 10 pixels of the line. The stdv can be expressed as a function of the synthetic concentration.
For real cases, the concentration is not available. A concentration map for each method will therefore be displayed. A second type of graph will consist of positioning the concentrations found by one method in relation to those of another method. The plume total mass will also be a point of comparison between the methods.

Ideal full-synthetic case: brick reflectance
The Figure 5 shows the retrieved concentrations according to the synthetic concentration for the different methods. As expected, the linear method, in red, underestimates the concentration above 10000 ppm.m for this case as shown in the insert. The underestimation with the linear method is increasing with increasing concentrations while the PTE and hybrid methods, respectively in blue and green, retrieve properly the synthetic (i.e. true) concentration. The insert of Figure 5 shows similar values of the stdv for the PTE and the hybrid methods about 800 ppm.m. The stdv of the linear method is close to 300 ppm.m. The differences between the PTE and hybrid methods are really small for this case.

Unfavourable full-synthetic case: water reflectance
For this unfavourable case due to low SNR, the hybrid method is more robust than PTE as shown in Figure 6. This figure presents zooms of the synthetic concentration map (A panel), the retrieved concentrations with the PTE (B panel) and the hybrid (C panel) methods. The PTE is not able to produce significant results over most of the plume while the hybrid method is providing a concentration map similar to the synthetic plume. The previous result is confirmed by Figure 7 where the retrieved concentrations are plotted for the three methods: the linear method underestimates the concentration, whereas the PTE method is not robust on this case. Indeed, the PTE method is based on the calculation of the plume transmission. In the case of very small radiances (SNR of about 10), the method does not succeed to isolate the plume signature in the spectrum and the estimation of the concentration fails for the most of pixels. The mean and standard deviation do not have a robust statistical significance in this case. The mean value of the retrieved concentration for the PTE method cannot be exploited even if it is closer to the synthetic concentration than linear method. The hybrid method provides the best results in terms of retrieved concentration level and stdv.
The Figure 8 represents the stdv calculated for the water case. The PTE curve values are not relevant due to the small number of estimated points. The mean stdv of the linear method stays the lowest but is 15 times larger for the water case than for the brick case (related to the SNR differences: 10 vs 150). For several runs of the algorithm, both stdv values remain close together but higher (factor 2) than the stdv of the linear method.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B3-2021XXIV ISPRS Congress (2021 The hybrid method based on OEM takes into consideration the image covariance matrix and the uncertainty on the background radiance rebuilding which are ignored by the PTE method. The same result is seen for the mean hybrid stdv (by a factor 11). From the OE formalism, the a posteriori uncertainty and the stdv are very similar which makes it possible the generation of a posteriori uncertainties in real case. The other surface types plotted in Figure 2 lead to the same conclusion: the higher the SNR, the better the results in terms of retrieved concentrations.

Semi-synthetic case
When a synthetic plume is added to a real reflectance image, the inter and intra-class variability is important. Using the HySpex reflectance map, a linear plume is added during the pass from reflectance to radiance image. The rebuilding uncertainties increase with the real reflectance map and the hybrid method must take this variation into consideration. Figure 9 shows the retrieved concentrations as for previous cases. The differences between the PTE and the hybrid methods are small. The hybrid method seems to retrieve concentrations slightly lower than with the PTE method. That might be a small bias effect of the underestimated a priori concentration used in the hybrid method, but this effect is not significant according to the stdv level. In that case, both methods are in good agreement with synthetic concentration. The Figure 10 presents the stdv of the three methods. For this case, the PTE has the lowest stdv of the two major methods. The differences are small in front of the level between mean PTE/hybrid methods and the linear method. The hybrid mean stdv was about 800 ppm.m for brick case and about 9000 ppm.m for water case. The mean hybrid stdv has the same value as for the water case. In the semi-synthetic case, the large stdv level is due to the reflectance heterogeneity of the scene whereas the low SNR level was responsible of the larger stdv in the water case.
Based on the OE analysis, the ratio of the a posteriori uncertainty and the retrieved concentration gives the relative uncertainty of the hybrid method. The Figure 11 presents this ratio for the brick (in orange) and the water (in blue) fullsynthetic cases and for the semi-synthetic case (in green), in percentage. The uncertainty is calculated for each concentration, which allows determining a threshold value for the minimum concentration to achieve such an uncertainty. The three cases show an exponential decrease with higher concentrations. For 2000 ppm.m, the uncertainty for the water case, the semi-synthetic case and the brick case are respectively about 68 %, 47 % and 38%. As expected, due to the good SNR, the brick full-synthetic case reaches the smallest uncertainty about 3%. The water full-synthetic and the semi-synthetic cases converge to 5%. Figure 11. Uncertainty of the hybrid method for the two fullsynthetic cases and the semi-synthetic case, in percentage.
These values show the robustness of the hybrid method to quantify a CH 4 plume even for rather dark pixels and low SNR. The validation of the method allows using it on real data. The results of this application are presented on the following section.

Results on a real methane plume
The airborne campaign data with the HySpex-NEO camera are acquired above a stack that emits methane with a flow rate of 75 g/s. During this campaign, the wind speed is extremely variable (Nesme, 2020). A mean wind speed of 2 m/s can be assumed. With the 1.4 m spatial resolution of the camera and this wind, the concentration expected for the plume is about 30000 ppm.m for the pixels close to the source. This value is in good agreement with the retrieved concentration of the PTE and the hybrid methods shown in the Figure 12. As previous, the linear method underestimates the source point concentration. Previously, the PTE method has shown limitations for unfavourable cases in terms of uncertainty. In addition, the linear method underestimates the retrieved concentrations for values over 5000 ppm.m. In order to compare the latter method with the hybrid method, Figure 13 shows the concentration from the hybrid method versus the concentration from the linear method. Above 5000 ppm.m, most of the points are above the first bisector which confirms the underestimation of the linear method observed throughout the paper. In terms of integrated mass enhancement (IME), the underestimation of the linear method is significant. The IME with the linear method is 16.7 kg compared to 21.6 kg for the hybrid, i.e. a mass difference of 30%.
Finally, the OE formalism allows the estimation of the a posteriori errors and uncertainties. For example, the Figure 14 represents the mean a posteriori uncertainty over the detected map as a function of a threshold value used to cut out the plume. For example, the mean uncertainty is about 31% when only pixels with more than 5000 ppm.m are considered. This figure shows a rapid decrease of the uncertainty: it is almost divided by two by keeping pixel above 1000 ppm.m and is reaching 20% for pixels above 15000 ppm.m Nevertheless, using only pixels with a concentration higher than 15000 ppm.m significantly reduces the number of pixels assigned to the plume itself. Figure 14. Uncertainty of the hybrid method for the real data as a function of the threshold concentration value use to calculate the mean uncertainty, in percentage.

CONCLUSION
The study of this paper compares three methods of CH 4 plume quantification (of the industry type). To do this, we benefited from the spectral sensitivity of the HySpex-NEO hyperspectral camera with a spectral resolution of 10 nm and 256 spectral channels. This camera was deployed during an airborne campaign above an industrial site that was emitting a CH 4 plume during the flight. We have tested the methods with semisynthetic data and full-synthetic data. The semi-synthetic data corresponds to a reflectance map of this industrial site (obtained by atmospheric correction) where a synthetic plume was added and converted in a radiance image. Synthetic reflectance map with various surface types have been used to test the methods. A brick reflectance with noise applied represents conditions with high SNR while a water surface with noise applied is representative of an unfavourable case with low SNR.
The linear method is based on the CTMF approach and linearizes the CH 4 transmission model. It leads to a bias on the retrieved concentration (underestimation). This expected result is commonly known but the robustness ans the short computation time of this algorithm makes it a very commonly used method in the literature. The PTE method begins with a rebuilding of the background radiance. The rebuilding step allows estimating the transmission of the plume by a differential calculation between an image with plume and a rebuilt image where the plume is absent. This method is much longer to run than the linear method, but the retrieved concentrations are in good agreement with the synthetic plume concentrations. Nevertheless, the success of this method is directly linked to the proper estimation of the transmission. In the case of a synthetic water surface, the low SNR does not provide robust results and a large part of the synthetic plume is not quantifiable. The hybrid method succeeds in quantifying the plume concentration where the PTE method fails. Moreover, the corresponding retrieved concentrations are generated with an estimation of the a posteriori uncertainty.
The results obtained with hyperspectral images of an airborne campaign are in good agreement with the flow rate of the source and the wind speed during the acquisition period. The total mass estimations present a difference about 30% between the LM and the hybrid method. The performance of the hybrid method can be estimated by the a posteriori uncertainty and the retrieved concentration values. If we consider only pixels with retrieved concentration higher than 5000 ppm.m, the uncertainty of the method is about 30 % and decreases to about 25% for pixels with columns higher than 10000 ppm.m. Future work will focus on the flow rate uncertainties estimation and the application of the hybrid method to satellite observations.