Methane retrievals from airborne HySpex observations in the shortwave infrared

Monitoring anthropogenic emissions is a crucial aspect in understanding the methane budget. Moreover, a reduction of methane emissions could help to mitigate global warming on a short timescale. This study compares various retrieval schemes for estimating localized methane enhancements around ventilation shafts in the Upper Silesian Coal Basin in Poland using nadir observations in the shortwave infrared acquired from the airborne imaging spectrometer HySpex. Linear and nonlinear solvers are examined and compared, with special emphasis put on strategies that tackle degeneracies between the surface reflectivity and broad-band molecular absorption features – a challenge arising from the instrument's low spectral resolution. Results reveal that the generalized nonlinear least squares fit, employed within the Beer InfraRed Retrieval Algorithm (BIRRA), can measure enhanced methane levels with notable accuracy and precision. This is accomplished by allowing the scene's background covariance structure to account for surface reflectivity statistics. Linear estimators such as matched filter (MF) and singular value decomposition (SVD) are able to detect and, under favorable conditions, quantify enhanced levels of methane quickly. Using k-means clustering as a preprocessing step can further enhance the performance of the two linear solvers. The linearized BIRRA fit (LLS) underestimates methane but agrees on the enhancement pattern. The non-quantitative spectral signature detection (SSD) method does not require any forward modeling and can be useful in the detection of relevant scenes. In conclusion, the BIRRA code, originally designed for the retrieval of atmospheric constituents from spaceborne high-resolution spectra, turned out to be applicable to hyperspectral airborne imaging data for the quantification of methane plumes from point-like sources. Moreover, it is able to outperform well-established linear schemes such as the MF or SVD at the expense of high(er) computing time.

1 Introduction

Methane (CH₄) is the second most important anthropogenic greenhouse gas next to carbon dioxide (CO₂), according to the IPCC (Intergovernmental Panel on Climate Change; Masson-Delmotte et al., 2021) report. Due to its comparatively short lifetime of approximately 9 years, a reduction of methane emissions could help to mitigate global warming on a relatively short timescale. Despite improvements in monitoring regional and global CH₄ emissions in recent years, the IPCC report points out that fundamental uncertainties pertaining to the methane budget remain (Intergovernmental Panel on Climate Change, 2014).

The vast majority of anthropogenic CH₄ emissions are caused by small-scale processes such as agriculture (enteric fermentation and manure), waste management (landfills), and fossil fuel exploitation, where the last is responsible for 20 %–30 % of all anthropogenic CH₄ emissions. Consequently, there is a need for continuous long-term methane observations on local to global scales in order to foster understanding on the global methane cycle, devise future reduction measures, and monitor their effectiveness. The monitoring of anthropogenic emissions of CH₄ and CO₂ is also part of the United Nations Framework Convention on Climate (2015), as nationally determined contributions should be assessed via global stock takes on a 5-year basis from 2023 (Article 13 and 14 of the Paris Agreement).

Satellite observations are typically used for continuous and global long-term monitoring of atmospheric composition, but also ground-based networks such as the Global Atmosphere Watch (GAW) Programme of the World Meteorological Organization (WMO) or the European Integrated Carbon Observation System (ICOS) are crucial assets. Spaceborne spectrometers measuring shortwave infrared (SWIR) solar radiation reflected at the Earth's surface are especially well suited to observe atmospheric CH₄ in the lower atmosphere by measuring its absorption around 1.6 and 2.3 µm. In contrast, the thermal infrared is less sensitive to variations in CH₄ concentration close to the surface, while mid-infrared sensors often have lower spatial resolution making them less favorable for emission monitoring (Richter, 2010).

Operational CH₄ products from contemporary atmospheric composition missions such as TROPOMI (TROPOspheric Monitoring Instrument; Veefkind et al., 2012) and GOSAT/GOSAT-2 (Greenhouse gases Observing SATellite; Kuze et al., 2009, 2016) measure trace gas concentrations with very high accuracy. Nevertheless, they are not optimally suited to measure emissions of point sources. This limitation is due to their focus on rapid global coverage, which entails a comparatively coarse spatial resolution of several square kilometers per pixel. Since the emission of a single point source inside a pixel is averaged over the entire resolution cell, even large sources rarely elevate the mean CH₄ concentration within one pixel by more than 1 % compared to the undisturbed background (Lauvaux et al., 2022). A way to increase the contrast of enhancements is to operate typical atmospheric remote sensing spectrometers at lower altitudes (e.g., on aircraft), thus increasing the spatial resolution while leaving the overall optical design untouched. This strategy is followed by instruments such as MAMAP/MAMAP-2D (Gerilowski et al., 2011) or GHOST (Humpage et al., 2018), which are well suited for the calibration and validation of their spaceborne counterparts.

Another way to increase the sensitivity towards smaller sources is to increase the instrument's spatial resolution. This in turn necessitates a trade-off in spectral resolution because the loss of photons caused by the smaller ground pixels reduces the signal-to-noise ratio (SNR) of the image, which has to be compensated for by broadening the width of the spectral channels. Imaging spectrometers for land surface remote sensing (often referred to as hyperspectral cameras) are typical examples of instruments optimized for spatial resolution this way. Their technology matured over the last 30 years, and a variety of airborne instruments and several spaceborne versions are either in orbit (Cogliati et al., 2021, PRISMA), (Guanter et al., 2015; Chabrillat et al., 2020, ENMAP) or going to be launched in the future (Rast et al., 2021, CHIME). Yet other sensors dedicated for the detection of methane or carbon dioxide, e.g., GHGSat (Jervis et al., 2021), CO2Image (Hochstaffl et al., 2023), or MethaneSat, have slightly higher spectral resolution than their hyperspectral counterparts but still offer a much higher spatial resolution than classical atmospheric composition missions.

Thorpe et al. (2013) were the first to demonstrate that localized CH₄ emissions over land can be detected from hyperspectral cameras with the Airborne Visible/Infrared Imaging Spectrometer (Green et al., 1998, AVIRIS) and that a limited quantitative analysis is possible (Thorpe et al., 2014). Similar studies were repeated with airborne instruments (AVIRIS-NG, Frankenberg et al., 2016; Duren et al., 2019; Borchardt et al., 2021; HySpex, Nesme et al., 2020) and spaceborne instruments (Thompson et al., 2016; Guanter et al., 2021). Varon et al. (2019) and Jervis et al. (2021) demonstrated that CH₄ sources can even be detected with the Multispectral Instrument (MSI) aboard the Sentinel-2 satellites, but these measurements are restricted to “favorable conditions” (i.e., strong sources and high surface albedo).

One of the core challenges when retrieving methane from measurements with high spatial (⪅100 m) and moderate spectral resolution (⪆1 nm) is the separation of spectral variations caused by molecular absorption and surface reflectivity (Ayasse et al., 2018). Classical methods for trace gas retrievals from high-spectral resolution instruments such as RemoTeC (Lorente et al., 2021), Weighting Function Modified Differential Optical Absorption Spectroscopy (Buchwitz et al., 2005, WFM-DOAS), or the Beer InfraRed Retrieval Algorithm (Gimeno García et al., 2011, BIRRA) exploit the high-frequency characteristics of gaseous absorption and attribute the smooth varying part to the surface albedo (and scattering). Instruments with coarse spectral resolution, however, are unable to sufficiently resolve those molecular signatures which causes ambiguities that often lead to surface-type-related biases in the classical retrieval schemes (e.g., Borchardt et al., 2021, Sect. 3.3, or Thorpe et al., 2014, Sect. 9.2). Alternative more data-driven retrieval schemes such as matched filter (MF) or singular value decomposition (SVD) employ methods from linear algebra and statistics that deal with spectral correlations and often yield results of sufficient accuracy (Thorpe et al., 2013, 2014; Thompson et al., 2015).

This study aims to compare concentration enhancements from different retrieval methods using measurements of the German Aerospace Center's (DLR) HySpex sensor system. The objective is to evaluate the retrievals' performance in terms of accuracy, precision, and speed and show advantages and drawbacks for each method. Another goal is to assess the latest BIRRA updates and its applicability to moderately resolved spectra from airborne sensors. Therefore, the paper is structured as follows.

First, the experimental setup is briefly described, followed by a quick review of atmospheric radiation and an introduction to the various BIRRA setups examined in this study. Afterwards, other simpler but faster retrieval schemes employed in this work are briefly discussed. The result section presents the CH₄ retrievals from HySpex observations over the Pniówek V ventilation shafts and compares the inferred concentrations and errors from the different methods. In the last section, results are summarized and put into perspective.

2 Methodology

Both linear and nonlinear methane retrieval schemes are examined. While the former are very fast but often lack sufficient accuracy, nonlinear iterative solvers require more computing power and time to come up with a best estimate. The study utilizes measurements collected by the DLR airborne HySpex sensor system (see Table 1) within the scope of the COMET (Carbon Dioxide and Methane) campaign on 7 June 2018 which focused on the detection and characterization of CO₂ and CH₄ sources in the Upper Silesian Coal Basin (USCB) in southern Poland.

Table 1 Summary of some important HySpex properties. The sensor is described in detail in Köhler (2016) and references therein.

To compare the performance of various retrieval methods, we limit our analysis to the two flight lines shown in Fig. 1, namely flight line 9 (called scene 09) and flight line 11 (called scene 11). The weather during the survey was well suited for remote sensing measurements. Apart from very few occasional patches of thin cirrus clouds, no further low- or mid-level clouds were near. Actual wind data for the USCB area are presented in Luther et al. (2022, Figs. 4 and 6).

Figure 1 (Left, top) Flight line 9 (dashed red line) was obtained around 09:55 UTC, while flight line 11 (solid red line) was acquired around 10:10 UTC (© OpenStreetMap contributors, 2022; distributed under the Open Data Commons Open Database License (ODbL) v1.0). The aircraft flew at an altitude of approximately ≈1200 and ≈2600 m above ground level, respectively, while heading eastward at 115^∘. (Left, bottom) Photograph of the ventilation shafts from the Pniówek V site. Photo credit: Leon Scheidweiler (Heidelberg University). (Right) False-color image from the SWIR-320m-e camera around the three Pniówek V shafts in scene 09.

Figure 2 (top) displays an ensemble of along-track-averaged HySpex observations. The spectral coverage of the HySpex SWIR-320m-e camera ranges from 967–2496 nm, with the exact number depending on the across-track detector ($\approx \pm \mathrm{1}$ cm⁻¹). The full width at half maximum (FWHM) of the SWIR-320m-e camera in the 1500–2500 nm (4000–6500 cm⁻¹) region ranges from 6.0–9.5 nm and is provided with the level 1b data for a sampling distance of 1.2 nm. The figure shows that the radiative intensity in the interval around 1.6 µm (≈6000 cm⁻¹) is significantly larger than the one around 2.3 µm (≈4300 cm⁻¹) mostly due to H₂O absorption (see Fig. 3). A possible bad pixel (number 104) is shown in the right plot around 1.65 nm (a descending cyan line). Also the surface reflectivity, depicted in Fig. 1, causes spectral variations in the observed radiance.

Figure 2 (a) HySpex average spectrum with the span (minimum to maximum) depicted in gray for measurements across the 320 across-track detectors for scene 09. (b) A bunch of individual spectra around 1.67 µm with the black lines indicating the pixel positions and sampling distance. The radiance values of pixel 104 (cyan) at ≈1.677 µm (5960 cm⁻¹), which is relevant for the CH₄ retrieval, appear to be problematic. (c) Reference reflectances for different surface types (measured at the Johns Hopkins University, Baldridge et al., 2009; Meerdink et al., 2019).

Figure 3 (a) Monochromatic transmissions of CH₄, CO₂, and H₂O for the SWIR spectral range and a nadir-looking observer at 1.5 km at a solar zenith angle (SZA) of 30^∘. (b) Total monochromatic transmission (black) vs. degraded to HySpex resolution (magenta). The spectral intervals used for the CH₄ fit are indicated by the yellow background.

2.1 Radiative transfer

In the SWIR spectral range the radiative transfer through the atmosphere under clear-sky conditions (cloud and scattering free) is well described by Beer's law (Zdunkowski et al., 2007) with the monochromatic transmission in wavenumbers ν given by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathcal{T}}_m(\mathit{\nu },s)& =\exp \left(-\sum _m{\mathit{\tau }}_m(\mathit{\nu },s)\right)\\ & =\exp \left(-\underset{\text{path}}{\int }\mathrm{d}s\sum _m{n}_m\left(s\right)k_m(\mathit{\nu },p(s),T(s\left)\right)\right).\end{array}\end{array}$$

The model assumes a pure gas atmosphere of molecules m, i.e., CH₄, CO₂, and H₂O. Optical depth τ_m is calculated by the path integral along s over the molecular number densities n_m and the pressure p and temperature T dependent absorption cross section k_m. The study utilizes the 2020 spectroscopic line data from GEISA (Gestion et Etude des Informations Spectroscopiques Atmosphériques; Delahaye et al., 2021) for molecular absorption calculations.

The decision to exclude aerosol modeling for HySpex observations was encouraged by findings from Borchardt et al. (2021), who concluded that different aerosol scenarios in the SWIR do not induce errors greater than 0.2 %. Moreover, since the spectra were observed at low flight altitudes on a rather clear day, retrieval errors induced by aerosol scattering should be negligible in our scenario as well (also see Fig. 3 and Thorpe et al., 2013; Thompson et al., 2015).

In Fig. 3, the top panel shows the individual components of the monochromatic total transmission for the US Standard Atmosphere, including methane's first overtone of the fundamental vibrational transition 2ν₃ (with its P and R branches) around 6000 cm⁻¹ (1560–1660 nm, tetradecad band), as well as additional strong absorption lines ranging from 4200–4600 cm⁻¹ (2090–2290 nm, octad). The bottom panel illustrates how the observer's coarse spectral resolution smooths the total monochromatic transmission (shown in black). There are 67 and 28 HySpex detectors used by the retrievals (see Fig. 3) within the range of 4100–4900 cm⁻¹ (4K) and 5700–6300 cm⁻¹ (6K), respectively.

2.2 Model atmosphere setup

The model atmosphere's vertical extent ranges from 0–80 km with 39 levels in total. The atmosphere is composed of pure gaseous layers. The highest vertical resolution is found in those layers below z_pl=2 km where the enhancement is expected to take place. The CH₄ optical depth is divided in two components, i.e., a climatological background τ_bg and a low-level (Gaussian) plume τ. The vertical profile of the initial guess plume is not crucial since nadir spectra in the SWIR do not contain sufficient information on the vertical distribution of trace gases (see Hochstaffl et al., 2020, Fig. 7).

The CH₄ and CO₂ initial guess background profiles are modeled according to the Air Force Geophysics Laboratory (AFGL) but adjusted to current concentrations in the atmosphere (1875 ppb and 400 ppm, respectively; cf. Sect. 2.3.2). H₂O and temperature and pressure are taken from reanalysis data provided by the National Centers for Environmental Prediction (Kalnay et al., 1996, NCEP).

2.3 Beer InfraRed Retrieval Algorithm (BIRRA)

The classical BIRRA level 2 processor, developed at DLR, uses the Generic Atmospheric Radiation Line-by-line Infrared Code (Schreier et al., 2014, GARLIC) as a forward model and a separate (SLS) or nonlinear least squares solver (NLS) for trace gas retrieval in the SWIR spectral region. It has been successfully applied to SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY; Gimeno García et al., 2011; Hochstaffl et al., 2018) and TROPOMI (TROPOspheric Monitoring Instrument, Hochstaffl et al., 2020) observations. In this study, however, the new Python version of BIRRA which is based on Py4CAtS (Python for Computational Atmospheric Spectroscopy, Schreier et al., 2019), a Python reimplementation of the validated Fortran code GARLIC (Schreier et al., 2013) is used.

The mathematical forward model Φ(x,ν) describes the measured intensity spectrum I(ν) for a nadir-looking observer according to

$$\begin{array}{}\text{(2)}& \mathrm{\Phi }(\mathit{x},\mathit{\nu })={\displaystyle \frac{r\left(\mathit{\nu }\right)}{\mathit{\pi }}}\cos \left(\mathit{\theta }\right)I_{\text{sun}}\left(\mathit{\nu }\right){\mathcal{T}}_m^{\downarrow }\left(\mathit{\nu }\right){\mathcal{T}}_m^{\uparrow }\left(\mathit{\nu }\right)\otimes S(\mathit{\nu },\mathit{\gamma }(\mathit{\nu }\left)\right),\end{array}$$

where r refers to the surface reflectivity and θ represents the solar zenith angle. The terms ${\mathcal{T}}_m^{\downarrow }$ and ${\mathcal{T}}_m^{\uparrow }$ denote the total transmission between the Sun and the reflection point (e.g., the Earth) and between the reflection point and the observer, respectively (see Eq. 1).

The transmission is described by

$$\begin{array}{}\text{(3)}& {\mathcal{T}}_m(\mathit{\nu },s)=\exp \left(-\sum _m{\mathit{\alpha }}_m{\mathit{\tau }}_m\left(\mathit{\nu }\right)\right),\end{array}$$

where the molecular scaling factors α_m adjust initial guess profiles. The simple scaling approach recognizes the significantly under-determined vertical profile information in the observed spectrum and enables an unconstrained least squares fit. All unknown (to be estimated) parameters are collected in the state vector x, which includes α_m and the polynomial coefficients for the surface reflectivity r_j (with j∈ℕ₀) (Gimeno García et al., 2011, Fig. 1). Finally, the instrument’s spectral response function (ISRF) is described by S. Its parameters such as the half width γ or a spectral shift can (optionally) be part of the state vector (also see Thorpe et al., 2014, Sect. 5.2).

2.3.1 Nonlinear solvers

This study examines various nonlinear retrieval schemes that are implemented in the BIRRA level 2 processor and are briefly introduced below. Nonlinear least squares methods are iterative and require calculating derivatives for each of the nonlinear state vector elements, represented by a Jacobian matrix J.

Nonlinear least squares (NLS) and separable least squares (SLS)

The nonlinear least squares fit minimizes the residual norm ($\Vert \cdot \Vert$ represents the 2-norm) for given measurements y when the model function Φ is nonlinear in one or more parameters of x according to

$$\begin{array}{}\text{(4)}& \underset{\mathit{x}}{min}{∥\mathit{y}-\mathbf{\Phi }\left(\mathit{x}\right)∥}^{\mathrm{2}}.\end{array}$$

The SLS splits (separates) the state vector x into nonlinear and linear parameters $\mathit{x}=(\mathit{\eta },\mathit{\zeta })$, where the elements in ζ enter the forward model Φ linearly (see Sect. 2.4.1). The minimization problem is hence given by

$$\begin{array}{}\text{(5)}& \underset{\mathit{\eta },\mathit{\zeta }}{min}{∥\mathit{y}-\mathbf{\Phi }\left(\mathit{\eta }\right)\mathit{\zeta }\left(\mathit{\eta }\right)∥}^{\mathrm{2}}.\end{array}$$

This setup is also known as the variable projection (VarPro, Golub and Pereyra, 2003) method, where η is independent of ζ in the matrix product Φ(η)ζ(η) (for details see Bärligea et al., 2023). The parameters in η can hence be fitted in the usual way by means of Gauss–Newton or Levenberg–Marquardt algorithms (see Hansen et al., 2013).

Generalized least squares (GLS)

A generalized least squares fit is used to account for correlated errors. The covariance matrix C encompasses spectral variations of the scene's background and the sensor's measurement noise. The motivation is that the matrix compensates for background variations that could be mistakenly attributed to methane band absorption. It is computed from background pixels that are assumed to not be affected by the CH₄ enhancements. This of course requires some information on the point source's location and prevailing wind direction.

The error covariance matrix C is a symmetric positive semi-definite matrix computed for each flight track. Figure 4 shows the square root matrix ${\mathbf{S}}^{-\mathrm{1}}={\mathbf{C}}^{-\mathrm{1}/\mathrm{2}}$ for the two methane retrieval intervals. In order to reduce fitting errors caused by degeneracies, S⁻¹ is included according to

$$\begin{array}{}\text{(6)}& \underset{\mathit{x}}{min}\Vert {\mathbf{S}}^{-\mathrm{1}}(\mathit{y}-\mathbf{\Phi }(\mathit{x}\left)\right){\Vert }^{\mathrm{2}}.\end{array}$$

Figure 4 Scene 09 inverse square root matrix of C: (a) 4100–4900 cm⁻¹ (4K) and (b) 5700–6300 cm⁻¹ (6K) spectral range. The background area was defined outside of the pixels: $\text{along-track}=(\mathrm{6300},\mathrm{6670})$; $\text{across-track}=(\mathrm{180},\mathrm{285})$. Note that beside the bad HySpex pixel mentioned in Fig. 2 at 5992.74 cm⁻¹ there appears to be another suspect pixel at 4691.04 cm⁻¹.

2.3.2 Enhancement estimates for the nonlinear solvers

A scene-averaged background spectrum, excluding ground pixels around the suspected CH₄ sources, is employed to estimate actual H₂O, CO₂, and CH₄ background concentrations. The CO₂ background level of the scene is inferred from the 1.6 and 2 µm bands via a multi-interval (4K and 6K spectral windows) fit. For scene 09 and scene 11, a scaling factor of ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{0.96}$ (≈385 ppm) and ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{0.93}$ (≈375 ppm) was determined for the initial guess, respectively. Due to the degeneracy between H₂O and the reflectivity polynomial at HySpex's spectral resolution, the scene-averaged H₂O scaling factor constitutes an effective parameter partly capturing low-frequency components in the spectrum. The scene-averaged CH₄ background profile was found to be within 5 % of the initial guess of 1875 ppbv; hence it is not (pre)scaled.

The state vector x for the CH₄ enhancement fit comprises the CH₄ scaling factor and the coefficients for a second-order reflectivity polynomial per spectral interval $\mathit{x}=(\mathit{\alpha },r_{\mathrm{0}},r_{\mathrm{1}},r_{\mathrm{2}})$. In this setup the parameter α only applies to the plume component (up to 2.0 km) of the CH₄ optical depth

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}}={\mathit{\tau }}_{\text{bg}}+\mathit{\alpha }\mathit{\tau }.\end{array}$$

This setup was found to be robust toward lower SNR values and less susceptible to correlations among state variables, which in turn enhances the condition number of the Jacobian matrix.

The actual CH₄ total column is then given by the background concentration plus the retrieved enhancement and includes corrections for light path modifications via the prefitted scene-averaged background CO₂ given by

$$\begin{array}{}\text{(8)}& N_{{\mathrm{CH}}_{\mathrm{4}}}=N_{\text{bg}}+{\displaystyle \frac{\mathit{\alpha }}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{{\mathrm{CO}}_{\mathrm{2}}}}}\widehat{N}\left(z_{\mathrm{0}}\right),\end{array}$$

with

$$\begin{array}{}\text{(9)}& \widehat{N}\left(z_{\mathrm{0}}\right)=\underset{z_{\mathrm{0}}}{\overset{z_{\text{pl}}}{\int }}{n}_{\text{pl}}\left(z\right)\mathrm{d}z\end{array}$$

and z₀ representing the bottom of the atmosphere and n_pl the plume's number density. This approach assumes that the CO₂ profile upon which ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{{\mathrm{CO}}_{\mathrm{2}}}$ is estimated corresponds to the true profile and that ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{{\mathrm{CO}}_{\mathrm{2}}}$ is 1 in absence of scattering.

2.4 Linear solvers

In contrast to nonlinear fitting schemes, linear solvers for x can only be used when equations can be expressed as a linear combination of the variables in x. To utilize such methods, it is usually required to linearize the forward model with respect to the variables of interest.

2.4.1 Linear least squares (LLS)

Assuming that the increase in optical depth caused by the plume, τ, is relatively small, the BIRRA forward model from Sect. 2.3 is linearized with respect to α by approximating the transmission spectrum of the plume by Taylor expansion according to

$$\begin{array}{}\text{(10)}& \exp \left(-\mathit{\tau }\right)\approx \left(\mathrm{1}-\mathit{\alpha }\mathit{\tau }\right).\end{array}$$

The linear least squares problem of M measurements can then be formulated according to

$$\begin{array}{}\text{(11)}& \underset{\mathit{x}}{min}\Vert \mathit{y}-\mathbf{\Phi }\mathit{x}\Vert ,\end{array}$$

where the model functions in Φ for the linear parameters of the state vector $\mathit{x}=(r_{\mathrm{0}},r_{\mathrm{0}}\mathit{\alpha })$ are given by

$$\begin{array}{}\text{(12)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}={\displaystyle \frac{\cos \left(\mathit{\theta }\right)}{\mathit{\pi }}}{I}_{\text{sun}}{\mathcal{T}}^{\downarrow }{\mathcal{T}}^{\uparrow }\otimes S,\text{(13)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2}}=-{\displaystyle \frac{\cos \left(\mathit{\theta }\right)}{\mathit{\pi }}}{I}_{\text{sun}}{\mathcal{T}}^{\downarrow }{\mathcal{T}}^{\uparrow }\mathit{\tau }\otimes S.\end{array}$$

It is important to note that in this setup the reflectivity coefficient r₀ is present in two elements of the state vector. In order to avoid this degeneracy and allow for higher-order reflectivity polynomials in the fit, which are required for large spectral intervals, the retrieval is performed in two steps. First, only the reflectivity coefficients are fitted, while in a second step only α is estimated with the prefitted reflectivity coefficients provided as input. The setup can be complemented by de-weighting individual pixels in the albedo fit that are impacted by methane. This approach minimizes interference between the two fits, preventing the reflectivity polynomial from capturing absorption of CH₄.

Another aspect that should be kept in mind is that since $\mathrm{1}-\mathit{\alpha }\mathit{\tau }\le \exp (-\mathit{\alpha }\mathit{\tau })$ for α≥0 the linearized model underestimates the CH₄ enhancement for a given optical depth τ compared to the nonlinear setup.

2.4.2 Matched filter (MF)

The MF is a well-established method for estimating molecular concentration enhancements from hyperspectral sensors, with numerous studies supporting its effectiveness (Villeneuve et al., 1999; Funk et al., 2001; Thorpe et al., 2013; Thompson et al., 2015). The linear enhancement factor is inferred by perturbing an average (background) radiance spectrum μ with a known target spectrum t. The approach is analogous to that used by Thompson et al. (2015), where CH₄ enhancements are estimated by linearly scaling a target signature that perturbs the mean radiance:

$$\begin{array}{}\text{(14)}& {\mathit{\alpha }}_i\left(\mathit{y}\right)={\displaystyle \frac{\left(\mathbf{t}\right(\mathit{\mu }\left){)}^{\mathrm{T}}{\mathbf{C}}^{-\mathrm{1}}\right({\mathit{y}}_i-\mathit{\mu })}{\sqrt{\left(\mathbf{t}\right(\mathit{\mu }\left){)}^{\mathrm{T}}{\mathbf{C}}^{-\mathrm{1}}\right(\mathbf{t}\left(\mathit{\mu }\right))}}}.\end{array}$$

This equation constitutes the linear minimizer that solves the Gaussian log likelihood:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\underset{{\mathit{\alpha }}_i}{min}\Vert {\mathbf{C}}^{-\mathrm{1}/\mathrm{2}}\mathit{d}{\Vert }^{\mathrm{2}}& \text{with}\mathit{d}={\mathit{y}}_i-(\mathit{\mu }+{\mathit{\alpha }}_i\mathbf{t}(\mathit{\mu }\left)\right)\\ & \text{and}\mathbf{t}\left(\mathit{\mu }\right)=-\mathit{\mu }\mathit{\tau }.\end{array}\end{array}$$

The method assumes that the measured spectrum can be represented as a linear superposition of the CH₄ plume's optical depth and the mean unperturbed radiance μ and tests an observed vector y_i against a base vector while accounting for the background covariance C. The mean background spectrum μ and C are computed per scene; C⁻¹ is approximated by decomposing C into eigenvalues and eigenvectors (Thompson et al., 2015, Eqs. 6–8).

In order to improve accuracy, a per-measurement target spectrum is computed, which accounts for the pixel's albedo (Foote et al., 2020, II. Methods, C.). This normalized matched filter includes an albedo factor r_i for each measurement spectrum according to

$$\begin{array}{}\text{(16)}& {\mathit{d}}_{\mathit{r}}={\mathit{y}}_i-(\mathit{\mu }+r_i{\mathit{\alpha }}_i\mathbf{t}(\mathit{\mu }\left)\right)\text{with}{r}_i={\displaystyle \frac{{\mathit{y}}_i^{\mathrm{T}}\mathit{\mu }}{{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }}}.\end{array}$$

However, the MF method has its limitations; for example, it suffers from a heterogeneous background and correlation between the plume and the background which limits the detection quality even for strong plumes (Theiler and Foy, 2006). According to Guanter et al. (2021), a way to mitigate the effect is by k-means clustering of the scene. This approach reduces within-class variance, which in turn should minimize the albedo sensitivity of α. In the so-called cluster-tuned matched filter, instead of a single background covariance statistic, a per-cluster background statistic C_i is computed (Thorpe et al., 2013; Nesme et al., 2020).

2.4.3 Singular value decomposition (SVD)

The retrieval of methane enhancements from hyperspectral AVIRIS data using singular vectors of the observed spectrum plus a target signature was first demonstrated by Thorpe et al. (2014). The SVD method is well suited for parameter estimation from moderately resolved spectral data because it allows us to consider only the most significant components of the spectrum while preserving the main spectral information.

The orthogonal singular vectors are obtained from HySpex spectra that are not impacted by the plume. The matrix containing the scene's log-space background spectra is decomposed into USV^T, where $\mathbf{U}\in {\mathbb{R}}^{m\times m}$ and $\mathbf{V}\in {\mathbb{R}}^{n\times n}$ are unitary matrices, and $\mathbf{S}\in {\mathbb{R}}^{m\times n}$ is a diagonal matrix. The target signature (spectrum) is represented by the CH₄ plume's optical depth τ, which is computed with Py4CAtS.

The basic idea is analogous to the MF, i.e., to represent the general variability in spectral radiance by a linear combination of singular vectors and a target signal. The minimization problem is then given by

$$\begin{array}{}\text{(17)}& \underset{\mathit{w}}{min}\Vert \mathit{y}-\mathbf{A}\mathit{w}{\Vert }^{\mathrm{2}}\text{with}\mathbf{A}\mathit{w}=\sum _k^N{\mathit{u}}_k{\mathit{w}}_k+\mathit{t}{\mathit{w}}_{{\mathrm{CH}}_{\mathrm{4}}},\end{array}$$

where A represents the concatenated matrix of the first N columns of the unitary matrix U (see Fig. 5). The vector w contains the corresponding weights, with $\mathit{\alpha }={\mathit{w}}_{{\mathrm{CH}}_{\mathrm{4}}}$ scaling the contribution of enhanced methane in the lowest atmospheric layers t=τ. In the cluster-tuned variant, the background spectra are clustered by k-means clustering, and the SVD is performed for each cluster separately. The respective base vectors per cluster are then used in the linear fit.

Figure 5 Standardized singular vectors u and the methane plume's target signature t in 4K (a) and 6K (b) spectral intervals, respectively. Standardization removes the mean and scales to unit variance. The u vectors are defined by the SVD and the t vector by the radiative transfer model Py4CAtS. Modeling the plume's optical depth with the same tools and for an equivalent setup (<2 km) is crucial for comparability with the nonlinear BIRRA setups.

2.4.4 Spectral signature detection (SSD)

A straightforward approach to identify methane absorption is the SSD, which compares the ratio of spectral residual norms to produce a score. Unlike other methods, this approach does not require any radiative transfer calculations, lookup tables, or initial guess information – only calibrated sensor data for a specific interval.

The algorithm is based on a simple polynomial fit of spectral pixels and the calculation of spectral residuals. The idea behind this method is similar to the continuum interpolated band ratio (CIBR) from Green et al. (1989) and Thompson et al. (2015, Eq. 2), which also measures absorption depths (Pandya et al., 2021). The method splits the spectral interval into pixels where CH₄ is absorbed and where it is not (or only weakly). A polynomial of degree P is fitted to the M out-of-band pixels:

$$\begin{array}{}\text{(18)}& \underset{x}{min}\left|\right|\mathit{y}-\mathbf{P}\mathit{x}|{|}^{\mathrm{2}}\text{with}P\in {\mathbb{R}}^{m\times (P+\mathrm{1})}.\end{array}$$

Next the residual norms for the in-band and out-of-band pixels are computed. The ratio of the residual norms yields an absorption band depth score for each observation which indicates variations in the CH₄ absorption given that the in-band and out-of-band pixels were properly chosen.

The algorithm constitutes a fast scheme which can also be applied for real-time detection of enhancements, e.g., determine whether or not a CH₄ ventilation shaft is active at the time of instrument overpass. When a zero-order polynomial is used for the out-of-band fit, the method is comparable to the CIBR algorithm. However, by using higher-order polynomials, the method can model the surface reflectivity and other interfering species more precisely, especially over larger spectral intervals.

3 Results

This section presents the results for the CH₄ estimates over the Pniówek V shaft(s). Except as otherwise stated, the retrievals were performed on 3×3 pixels averaged spectra in order to increase the signal-to-noise ratio and thereby reduce clutter of the CH₄ fits across pixels.

3.1 NLS and SLS fits

Figure 6 shows the results of the classical BIRRA NLS fit. The position of the source is indicated by the intersection of the dashed line. The fits reveals a significant enhancement of CH₄ in both spectral intervals downwind of the ventilation shaft. However, both BIRRA configurations exhibit biases, with the SLS fit displaying a somewhat more pronounced sensitivity to surface variations (therefore not shown). As depicted in Fig. 6, the combination of multiple spectral intervals can alleviate these adverse effects to a considerable extent, and the downwind shape of the plume is captured better (see Table 2).

Figure 6 Methane enhancements for 3×3 spatially averaged HySpex observations in the (a) 4150–4900 cm⁻¹ (4K) interval and the (b) 5700–6300 cm⁻¹ (6K) range. (c) Multi-interval fit, i.e., combining the 4K and 6K ranges (4K6K). Note that the former two NLS fits suffer from albedo correlations with methane in opposite direction.

Table 2 Mean and standard deviation (SD) for the background pixels of the t test. Relating the standard deviation to the mean is a good indication of the accuracy and precision of a method.

3.2 GLS fits

Figure 7 displays the retrieved columns using the generalized least squares (GLS) fit from averaged spectra for scene 09 in the 4K and 6K intervals. Compared to other methods, it reduces the correlation between the methane enhancement and surface reflectivity significantly, resulting in a more distinct plume signal and less background clutter.

Figure 7 Methane plume depicted for the single-window covariance-weighted fits for scene 09. The background pixel concentration is rather stable in the 4K interval depicted in panel (a), while there is still some overestimation of CH₄ in the 6K range in panel (b) which might partly be caused by the bad pixel close to the methane absorption band (see Figs. 2 and 4).

Figure 8 shows the multi-window covariance-weighted GLS fits for scene 09 and 11. In both cases the retrieval yields a distinct plume that separates well from background clutter. The figure depicts the impact of decreasing ground pixel resolution (from higher altitudes) on the inferred concentrations as enhancements are less pronounced for scene 11. However, this could also partly be attributed to a decreased amount of emissions since the observation was taken at another point in time. Furthermore, winds could have changed as the plume's shape is different compared to scene 09.

Figure 8 (a) Multi-window (4K6K) retrieval output for scene 09 and (b) enhancements for scene 11. The stripe pattern in the along-track direction is a multi-window retrieval artifact.

Figure 9 depicts the fits from individual (non-averaged) HySpex spectra for scene 09 and 11 for the GLS multi-window retrieval setup. The single-pixel total columns are more affected by retrieval noise caused by the lower signal-to-noise ratio (SNR), which varies significantly over different surface types. However, the method still identifies elevated methane concentrations.

Figure 9 In the single-pixel spectra depicted for (a) scene 09 and (b) scene 11, retrieval noise is significantly dependent on the underlying surface.

3.3 MF fits

The albedo-normalized, cluster-tuned, and classical matched filters are examined for scene 09. Figure 10 shows that all three variants are able to identify the methane plume, although absolute CH₄ concentrations differ in certain parts of the scene. The cluster-tuning MF variant in the middle panel yields more homogeneous enhancements downwind across various surface types, but pixels along class boundaries such as streets show some artifacts.

Figure 10 (a) Albedo-normalized MF, (b) the cluster-tuned variant, and (c) the classical MF fit shown for the 4K interval (4100–4900 cm⁻¹).

3.4 SVD fits

The SVD-based retrieval method illustrated in Fig. 11 is able to identify elevated levels of CH₄ in the HySpex spectrum. The method yields consistent results for both spectral intervals, employing four base vectors and the CH₄ Jacobian for the lowest 2 km (see Fig. 5). Including more than four base vectors significantly increases the condition number of A as column five interferes with the methane signal. The plume is also identified for the purely data-driven approach, where the base vector mimicking the CH₄ absorption (the fifth column in U) is used instead of the CH₄ spectrum. Thus, this approach does not require any forward model and is hence purely data-driven. Cluster tuning in general improves the fit due to a reduction in variance within each cluster; however, the results become more sensitive to the selected number of base vectors. It was found that within a cluster the number of base vectors required to resemble A should be reduced.

Figure 11 (a) Standard SVD fit and (b) background cluster-tuned SVD, both for the 6K spectral range. Three clusters reduce background clutter but suppress some enhancements close to the source. However, also false positives like the spot around the coordinate (200,6350) are diminished.

3.5 LLS fits

The linear least squares fit is able to identify CH₄ enhancements, although it differs significantly in the absolute values in the 4K and 6K spectral range (see Fig. 12). As pointed out in Sect. 2.4.1, the method is prone to underestimating enhancements. Moreover, the selected weights for the reflectivity coefficient fit were found to impact the CH₄ result. However, for the sake of simplicity and since the optimal selection of weights was not clear initially, no weighting was applied. Similar to its nonlinear counterpart (NLS) the fit is also affected by albedo-related offsets in opposite directions in the two intervals. However, relative enhancements between plume and background values are rather similar.

Figure 12 CH₄ enhancements for scene 09 estimated with the LLS setup. The results in panel (a) show the results for the 4K spectral window, while panel (b) shows the 6K outcome. In the latter method, the methane enhancements are less pronounced, but the reflectivity-related bias is also smaller.

3.6 SSD fits

In Fig. 13, results for the SSD method are shown. The results show that relative variations are more pronounced in the zero-order fit, while the higher-order fit better captures the downwind plume by suppressing background clutter.

Figure 13 The ratio of the spectral residuals in the 6K range for the in-band and out-of-band pixel is depicted. In panel (a) the in-band residuals were computed with respect to a quadratic polynomial, while in panel (b) a constant was used.

It is important to note that the method yields better results for the 6K absorption since the 4K absorption features are distributed over a larger spectral range which causes more uncertainty in the out-of-band polynomial fit since many pixels need to be omitted.

3.7 Statistical significance of results

In order to provide a more quantitative measure on the quality and confidence of the fits, a Student t test was applied to the results (Varon et al., 2018). The test helps to measure how well the plume is represented with respect to the background for a given retrieval setup. This is accomplished by testing for pixels that contradict the null hypothesis, which assumes that all pixels belong to the background (methane concentration). Moreover, samples need to be independent and identically normally distributed (Bruce et al., 2020).

The null hypothesis was rejected at the 1 % significance level, which can be considered a strong evidence. Although some fit results may ask for a tighter significance level in the t test to isolate the plume and get rid of most outliers, for the sake of comparison 1 % is used throughout this study.

Figure 14 depicts the result of the t test applied to the retrieval output for scene 09 and 11 from the covariance-weighted nonlinear solver (GLS) in the 4K range. The plume is well pronounced, and the test is able to isolate enhanced CH₄ values from the background. In particular, the higher-ground-resolution scene 09 shows almost no outliers at the selected significance level, indicating that the depicted values occur only in ≤1 % of the cases, assuming the null hypothesis (background methane concentrations) holds.

Figure 14 Plume pixels according to the Student t test for the nonlinear multi-window GLS fit. Panel (a) shows scene 09, while panel (b) shows scene 11.

The Student t test was also applied to the linear solvers, with results reported in Fig. 15. The test was performed with the same significance level set for the previous cases. Each of the linear methods provides enough pixels within the confident range to isolate the plume pixels. While MF provides the most accurate enhancement values compared to the GLS (see Fig. 14), the SVD better captures the downwind plume; however, peak enhancements are ≈30 % lower. The LLS method does capture the downwind plume but is much less sensitive to enhancement as it significantly underestimates these.

Figure 15 Plume pixels identified by the t test in scene 09 for the different linear schemes. Note that the color scale was adapted. (a) Classical MF output from the 4K range, (b) plume pixels according to the SVD method in the 6K interval, and (c) the LLS fit in 6K.

3.8 Errors and correlations

In general the retrieval's fit quality is assessed with respect to the discrepancy between the measurement y and the modeled spectrum I according to $\mathit{\sigma }=\Vert \mathit{y}-\mathit{I}\left(\widehat{\mathit{x}}\right)\Vert$, also known as the residual norm. In order to get the uncertainties in the estimates of the model parameter $\widehat{{\mathit{x}}_i}$ for a particular fit, the residual norm is multiplied by the least squares covariance matrix:

$$\begin{array}{}\text{(19)}& \mathit{V}={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{M-N}}{\left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}\right)}^{-\mathrm{1}}\text{with}\widehat{{\mathit{x}}_i}\pm \sqrt{{\mathit{V}}_{\mathit{i}\mathit{i}}},\end{array}$$

with J representing the Jacobian matrix for the state vector $\widehat{\mathit{x}}$.

The errors of the individual state vector parameters are represented in the square root of the diagonal elements of V. The standard error $\sqrt{{\mathit{V}}_{\mathbf{11}}}$ for the CH₄ scaling factor is shown in Fig. 16. The uncertainty varies with different surface types according to Eq. (19). A different way to evaluate the quality of the retrieval for a scene is to estimate the fit error from the variability of pixels identified as background by the t test. This method calculates a score by comparing the means of pixels from the target area and the background area and dividing this by the standard deviation of the background. These values are also obtainable for all the linear fit variants.

Figure 16 Uncertainties in the estimated CH₄ according to Eq. (19) for the covariance-weighted fit in the (a) 4K and (b) 6K spectral windows. The 6K range shows larger errors as it contains less than half the number of pixels than the 4K window. Beside the higher variances, the bad pixel close to the methane lines around 1.6 µm also increases the spectral residual norm.

Tables 2 and 3 present the findings of this analysis for the nonlinear and linear solvers, respectively. The analysis shows that the GLS fit performs best and that SLS and NLS yield similar results, while the MF scores highest amongst the linear solvers. In accordance with Fig. 16, fits in the 4K window score higher compared to the 6K. The less sensitive the retrieval is to CH₄ enhancements, the less variations will be observed in the background. Therefore, the standard deviation in the last column should not be overemphasized in the evaluation of the setups.

Table 3 Same as Table 2 but for the linear retrieval setups.

Figure 17 shows the correlation matrix of the retrieval outputs for the different solvers and spectral intervals. It reveals that most solvers have rather good correlations with the GLS solver (sort of benchmark), particularly in the 4K and multi-window 4K6K spectral ranges. Moreover, the GLS, MF, and SVD show blocks of high correlation. Blue colors indicate that inferred concentrations tend to move in opposite directions, which is the case for example in the single-window NLS fits shown in Fig. 6.

Figure 17 Pearson correlation coefficients for inferred methane from scene 09 for the nonlinear and linear solvers in the two examined spectral intervals.

4 Discussion

The study found that nonlinear setups which utilize background pixel covariance statistics (GLS) are suited to quantify CH₄ concentrations with good accuracy and precision and should also allow us to quantify emissions. The NLS and SLS fits encounter challenges due to degeneracies between the surface reflectivity and the broad-band molecular absorption signal at the HySpex resolution. In accordance with Borchardt et al. (2021) and Guanter et al. (2021), surface brightness and homogeneity were found to be important factors in detecting and quantifying methane plumes. The retrieval noise can vary significantly depending on surface characteristics. A given type of surface can lead to a positive bias in one spectral window, while the opposite may be true in another window (see Fig. 9).

In order to scan for potential CH₄ leakages on large datasets with millions of pixels, linear solvers such as the SVD, MF, or LLS are more appropriate due to their significantly higher speed (Thompson et al., 2015). While the iterative setups require roughly 1 s per fit, the linear methods are up to 3 orders faster. In particular, the SVD and MF solvers yield enhancements that often agree well with the more sophisticated nonlinear methods, although their sensitivity can be significantly hampered by the lack of uniformity in background reflectance (Thorpe et al., 2014; Foote et al., 2020).

The LLS fit ignores background statistics, and hence the inversion suffers from albedo correlations similar to its nonlinear counterparts. The fit also significantly underestimates enhancements, although it is able to capture parts of the pattern.

Polynomials up to the second order were able to capture the enhanced methane signal in the rather simple SSD method. The selection of an adequate polynomial needs to consider the width of the spectral interval and its surface reflectivity. Moreover, the method is not designed to quantify methane but only allows for the detection of anomalies in the spectral residuum.

Cluster tuning of linear retrieval setups can help to mitigate background clutter and surface-reflectivity-induced biases (Nesme et al., 2020), especially if large background areas are selected. Predicting the right cluster for pixels impacted by the methane plume is crucial for this method in order to improve fit results. Cluster tuning can moreover be a beneficial preprocessing step as it allows us to reduce the base vectors per cluster in the SVD method so that fewer base vectors are sufficient to model the background spectrum. However, in this case a separate model matrix A needs to be compiled for each cluster.

While linear methods are well suited to survey vast datasets and pinpoint potential sources, iterative solvers such as BIRRA are adequate to quantify enhanced concentrations at known locations as the slower speed is not of much concern for some thousands of observations.

5 Conclusions

The study examines the feasibility of methane retrievals from hyperspectral imaging observations for various retrieval methods. It was found that localized CH₄ enhancements close to the ground can be quantified from HySpex airborne observations. The generalized BIRRA retrieval is well suited to investigate potential methane emissions. The statement is underpinned by the relatively low background variations and distinct CH₄ enhancement pattern in the surface-albedo covariance-weighted BIRRA fits in (see Figs. 7 and 8 and Table 2).

The BIRRA NLS and SLS fits were found to be sensitive to spectral variations in the albedo, leading to surface-type-dependent biases known from previous studies utilizing data from hyperspectral sensors. This effect is more pronounced for single spectral intervals but less evident when multiple intervals are combined.

The linear estimators proved to be highly efficient and effective for many cases, making them suitable in the survey of large hyperspectral datasets. The well-established SVD and MF method produced results that often agree well with the BIRRA inferred enhancements; however, the sensitivity is lower. The LLS method turned out to be the least sensitive one. For detection purposes the SSD was found to be a useful tool.

In conclusion, covariance-weighted methods are able to quantify methane enhancements from hyperspectral SWIR observations at high spatial resolution with good accuracy. In particular, the GLS solver is suited to capture enhancements with an accuracy that should allow for emission estimation. Considering the significant speedup and reasonable accuracy of the linear methods MF and SVD, both constitute a valuable tool in examining plumes on vast datasets.

The methods are also applicable to spaceborne sensors, which will be considered in a next step. Overall, the new Python version of the BIRRA code used in this study turned out to be a flexible toolbox for prototyping.