Recovery of sparse urban greenhouse gas emissions

To localize and quantify greenhouse gas emissions from cities, gas concentrations are typically measured at a small number of sites and then linked to emission fluxes using atmospheric transport models. Solving this inverse problem is challenging because the system of equations often has no unique solution and the solution can be sensitive to noise. A common top–down approach for solving this problem is Bayesian inversion with the assumption of a multivariate Gaussian distribution as the prior emission field. However, such an assumption has drawbacks when the assumed spatial emissions are incorrect or not Gaussian distributed. In our work, we investigate sparse reconstruction (SR), an alternative reconstruction method that can achieve reasonable estimations without using a prior emission field by making the assumption that the emission field is sparse. We show that this assumption is generally true for the cities we investigated and that the use of the discrete wavelet transform helps to make the urban emission field even more sparse. To evaluate the performance of SR, we created concentration data by applying an atmospheric forward transport model to CO₂ emission inventories of several major European cities. We used SR to locate and quantify the emission sources by applying compressed sensing theory and compared the results to regularized least squares (LSs) methods. Our results show that SR requires fewer measurements than LS methods and that SR is better at localizing and quantifying unknown emitters.

1 Introduction

Understanding anthropogenic greenhouse gas (GHG) emissions is important for scientists and decision makers fighting climate change. Based on a growing amount of atmospheric observations, studies estimating emission fields of GHG sources and sinks from these observations have been performed on local (Chen et al., 2016; Viatte et al., 2017; Toja-Silva et al., 2017), metropolitan (Jones et al., 2021; Turner et al., 2020; Hase et al., 2015), country (Miller et al., 2013; Shekhar et al., 2020), and global (Hirsch et al., 2006; Mueller et al., 2008; Turner et al., 2015; Jacob et al., 2016) scales. One of the main reasons for such studies is to verify and improve GHG emission inventories created by bottom–up methods. Verification and improvements include but are not limited to

• determining the difference between the real emissions and the emissions captured by inventories

• determining differences between the real and bottom–up estimated emissions for individual emitters

• finding emitters which are not captured by inventories (unknown emitters).

Atmospheric inverse modeling methods use column or in situ GHG concentration measurements to estimate emission fields. Due to a lack of measurements and high modeling and measurement uncertainties, estimating each grid cell of an emission field independently is not possible. Instead, sectors (Jones et al., 2021), spatial correlations (Wesloh et al., 2020), and/or temporal correlations (Jones et al., 2021; Wesloh et al., 2020) are used to construct alternative parameterizations of emission fields to prevent overfitting.

An alternative to overcoming these issues are sparse reconstruction (SR) methods (Ray et al., 2015). SR methods can use concentration measurements to estimate sparse emission fields, meaning that only a small number of large emitters contribute significantly to the total emissions. These methods determine the critical emission grid cells and adjust the emissions of only those cells until the model best matches the observations. All other grid cells are set to zero. Once conditions of compressed sensing (CS) are fulfilled, SR methods are guaranteed to determine the best possible emission fields and provide a good estimation of their emissions.

SR for the recovery of GHGs has been proposed by several recent studies. Ray et al. (2014, 2015) used stagewise orthogonal matching pursuit (StOMP), a reconstruction method known from compressed sensing, and modified it to enforce positive emission estimates. To overcome the restriction of sparse reconstruction of only being able to reconstruct sparse fields, they used a multi-scale resolution field based on wavelets. In Ray et al. (2015), StOMP has been modified so that prior information can be included. Both studies reconstructed fossil fuel CO₂ emissions in an idealized scenario with synthetic measurements and very low measurement noise (SNR >40 dB). Hase et al. (2017) demonstrated sparse reconstruction with enforced positive emission estimates of anthropogenic CH₄ emissions from synthetic observations in the US. To increase the sparsity of the emission field, Hase et al. (2017) used a redundant dictionary representation, where the representation of the emission field is not unique. Yan et al. (2012) proposed compressed sensing, i.e., sparse reconstruction with guaranteed best feature selection, for environmental monitoring with the focus on undersampling.

In this paper, we apply SR to assessing urban GHG emissions. Current GHG emission monitoring systems in cities, such as Dietrich et al. (2021), Shusterman et al. (2016) and Sargent et al. (2018), acquire GHG concentration in the atmosphere as column or in situ concentration measurements. These measurements are then related to city emissions and background concentrations, where the city emissions are the unknowns of interest while the background is (partially) known. Göckede et al. (2010) have shown that in smaller domains, uncertainties in the background have a high influence on the estimation of the city emissions. Therefore, modern approaches, such as Jones et al. (2021) and Klappenbach et al. (2021), use the measurements acquired to additionally improve the certainty of background concentrations using a Bayesian approach. In this work, we ignore the background and make the assumption that it is known in full detail. Extending our approach to include background concentrations is straightforward.

As urban emissions, we are using anthropogenic emission inventories from multiple European cities. To overcome the sparsity constraints of SR for non-sparse emissions, we use a wavelet transformation.

We are the first to apply SR to the estimation of urban GHG emissions. The findings of our work are the following:

• Urban emissions are mostly sparse and a third-level wavelet transform performs well in sparsifying urban emissions further.

• SR needs fewer measurements than Gaussian prior methods to achieve a similar performance if the emissions are sparse enough.

• SR performs well in localizing and quantifying large emitters, leading to the application of finding unknown emitters not captured by emission inventories.

The paper is structured as follows. Section 2 gives a formulation of inverse problems, introduces the reader to sparse reconstruction methods, compressed sensing, and compressible emissions, and provides a description of the algorithms used in this paper. The compressibility of the anthropogenic emissions in European cities is discussed in Sect. 3. Section 4 shows selected scenarios of our reconstruction method, highlighting beneficial conditions and use cases of our reconstruction method as well as discussing measurement noise. In Sect. 5 we create case studies for different European cities in an idealized and noisy case.

2 Methodology

This section gives the problem statement of atmospheric inverse problems (Sect. 2.1), provides an introduction to the theory of sparse reconstruction (Sect. 2.2, Sect. 2.3, and Sect. 2.4), and introduces measures (Sect. 2.5) and algorithms (Sect. 2.6) for sparse reconstruction.

2.1 Inverse problems

An inverse problem is a problem in which input parameters should be determined from the observation of a process. For the problem in this paper, those input parameters x∈ℝⁿ are the GHG emission fluxes for each grid cell in an emission field and the measurements y∈ℝ^m are in situ or column measurements of GHG concentrations in the atmosphere. These quantities are connected by an atmospheric process, referred to as forward model F: y=F(x). In practice, such a forward model can be a linear, non-linear, or even a stochastic process. For this analysis, we limit the forward model to linear cases. Therefore, we can write y=Ax, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ is called the sensing matrix. A least squares estimation of the GHG emission fluxes x is given by

$$\begin{array}{}\text{(1)}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}{∥\mathbf{A}\mathit{x}-\mathit{y}∥}_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

Often, however, such inverse problems are ill-posed, as in the cases we deal with in this paper. In such cases, no or no unique solution exists or the solution does not depend smoothly on the data, therefore, being sensitive to noise. For ill-posed inverse problems, the least squares estimation without regularization does not provide a useful reconstruction technique. For a more detailed discussion of ill-posed problems we refer to chap. 3 of Nakamura and Potthast (2015).

2.2 Bayesian inversion

A typical approach in atmospheric sciences to solve inverse problems is Bayesian inversion. In such a setup, the unknown emissions are assumed to follow a known probability distribution. This probability distribution is referred to as a priori. Measurements are used to update the a priori, which results in a new probability distribution referred to as a posteriori. From this posteriori distribution a parameter estimation can be made using a maximum likelihood (ML) detector on the a posteriori. Since the ML detector acts on the a posteriori, this is commonly referred to as Maximum a posteriori (MAP) detector. Let us call the probability distribution of the a priori p_X(x). Furthermore, probability distributions are assigned both to the observations and the model to account for uncertainties. The measurement distribution is written as p_Y(y) and the distribution which maps x to y is written as $p_{Y|X}\left(\mathit{y}\right|\mathit{x})$. Using Bayes' theorem, a posteriori distribution for x under the condition of the observations y can be derived:

$$\begin{array}{}\text{(2)}& p_{X|Y}\left(\mathit{x}\right|\mathit{y})={\displaystyle \frac{p_{Y|X}\left(\mathit{y}\right|\mathit{x}\left)p_X\right(\mathit{x})}{p_Y\left(\mathit{y}\right)}}.\end{array}$$

On this derived distribution the MAP detector can be applied. Since p_Y(y) is only a constant factor for some specific measurements, the MAP detector can be written as

$$\begin{array}{}\text{(3)}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmax}}{p}_{Y|X}\left(\mathit{y}\right|\mathit{x}\left)p_X\right(\mathit{x}),\end{array}$$

where $\widehat{\mathit{x}}$ are the estimated physical quantities. Typically the distributions are assumed to be Gaussian, where the a priori distribution is assumed to be centered around an initial guess, x_A. This allows for easier error analysis and makes the problem computationally feasible.

2.3 From Bayesian inversion to regularization

To show the relation between Bayesian inversion and regularization methods, we show how a Bayesian inversion problem, using Gaussian priors, can be converted to a regularization problem. Assume that $p_{Y|X}\left(\mathit{y}\right|\mathit{x})$ is Gaussian distributed with the covariance matrix S_o,

$$\begin{array}{}\text{(4)}& p_{Y|X}\left(\mathit{y}\right|\mathit{x})={\displaystyle \frac{\mathrm{1}}{\sqrt{\left(\mathrm{2}\mathit{\pi }{)}^m\text{det}\right({\mathbf{S}}_o)}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥{\mathbf{S}}_o^{-\mathrm{1}/\mathrm{2}}(\mathbf{A}\mathit{x}-\mathit{y})∥}_{\mathrm{2}}^{\mathrm{2}}\right),\end{array}$$

where m are the number of measurements, and p_X(x) to be a Gaussian prior of the form

$$\begin{array}{}\text{(5)}& p_X\left(\mathit{x}\right)=C\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{∥{\mathbf{S}}_A^{-\mathrm{1}/\mathrm{2}}(\mathit{x}-{\mathit{x}}_A)∥}_{\mathrm{2}}^{\mathrm{2}}\right),\end{array}$$

where C is a normalization constant. Applying the MAP detector from Eq. (3) gives an estimation of

$$\begin{array}{}\text{(6)}& {\displaystyle }\widehat{\mathit{x}}& {\displaystyle }=\underset{\mathit{x}}{\mathrm{argmax}}{p}_{Y|X}\left(\mathit{y}\right|\mathit{x}\left)p_X\right(\mathit{x})\text{(7)}& {\displaystyle }& {\displaystyle }=\underset{\mathit{x}}{\mathrm{argmin}}\left[{∥{\mathbf{S}}_o^{-\mathrm{1}/\mathrm{2}}(\mathbf{A}\mathit{x}-\mathit{y})∥}_{\mathrm{2}}^{\mathrm{2}}+{∥{\mathbf{S}}_A^{-\mathrm{1}/\mathrm{2}}(\mathit{x}-{\mathit{x}}_A)∥}_{\mathrm{2}}^{\mathrm{2}}\right].\end{array}$$

The idea of regularization methods on the other hand is to add a penalty term to the least squares problem from Eq. (1) to prefer solutions of a certain kind,

$$\begin{array}{}\text{(8)}& \underset{\mathit{x}}{\mathrm{argmin}}\left[{∥{\mathbf{C}}_{\mathrm{1}}\left(\mathbf{A}\mathit{x}-\mathit{y}\right)∥}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\lambda }R\left({\mathbf{C}}_{\mathrm{2}}\mathit{x}\right)\right],\end{array}$$

where $R:{\mathbb{R}}^n\to \mathbb{R}$ is the regularization function and C₁,C₂ are correlation matrices. This equation is equivalent to Eq. (7) with λ=1, ${\mathbf{C}}_{\mathrm{1}}={\mathbf{S}}_o^{-\mathrm{1}/\mathrm{2}}$, $R={∥{\mathbf{S}}_A^{-\mathrm{1}/\mathrm{2}}(\mathit{x}-{\mathit{x}}_A)∥}_{\mathrm{2}}^{\mathrm{2}}$, and C₂=I. For such a regularization function, the regularization scheme is known as Tikhonov regularization or also ridge regression (Golub et al., 1999). While in Bayesian inversion it is most often assumed that a prior value x_A is known, in regularization x_A is often unknown and assumed to be 0.

In this paper, we investigate sparse reconstruction (SR) methods. To achieve SR, the regularization term has to be changed so that sparse solutions are preferred over non-sparse solutions. In statistics, such a regularization function is the Lasso regularization function, presented by Tibshirani (1996). The Lasso is given by

$$\begin{array}{}\text{(9)}& R\left(\mathit{x}\right)=\sum _j\left|{\mathit{x}}_j\right|={∥\mathit{x}∥}_{\mathrm{1}}\end{array}$$

and is used especially when x is approximately sparse. The equivalent of the Lasso in Bayesian inversion is the assumption of a Laplacian distributed prior. The Lasso is expected to select those elements in x which are important and meaningful, while the irrelevant features are estimated to be zero. Su et al. (2017) showed that this is not necessarily the case, as coefficients which are zero in x are sometimes estimated to be important (which is referred to as false discovery). In the next section, we introduce compressed sensing (CS), which provides sufficient conditions to prevent false discoveries in x using Lasso regularization.

2.4 Compressed sensing

Compressed sensing (CS) is a theory which provides sufficient conditions to guarantee best possible reconstruction using the Lasso regularizer, therefore, preventing false discoveries. The conditions of CS apply to the forward model A and are hard to examine. For this reason, the conditions are normally already considered in the design process. In the following, we provide the very basics of CS needed to understand our work. For a more comprehensive introduction to CS, we refer to Boche et al. (2015).

CS states that an s-sparse signal $\mathit{x}\in {\mathrm{\Sigma }}_s^n$, where s-sparse is defined as $s\ge \left|\left\{j|{\mathit{x}}_j\ne \mathrm{0}\right\}\right|$ and ${\mathrm{\Sigma }}_s^n$ is the set of all s-sparse signals which are in ℝⁿ, can be uniquely reconstructed by m measurements y∈ℝ^m, defined by y=Ax, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, if certain conditions are satisfied for A. An overview of the symbols, norms, and measures used is given in Table 1. The unique solution is found by solving the l₀ minimization problem given by

$$\begin{array}{}\text{(P0)}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}{∥\mathit{x}∥}_{\mathrm{0}}\text{ s.t. }\mathbf{A}\mathit{x}=\mathit{y}.\end{array}$$

Solving this minimization problem is NP-hard and not applicable to real-world applications. Candès et al. (2006) showed that for additional conditions in A, one can solve the l₁ regularization problem instead, which is a convex problem:

$$\begin{array}{}\text{(P1)}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}{∥\mathit{x}∥}_{\mathrm{1}}\text{ s.t. }\mathbf{A}\mathit{x}=\mathit{y}.\end{array}$$

Then, solving Eq. (P1) leads to the same solution as Eq. (P0). A sufficient condition for recovering Eq. (P0) with Eq. (P1) is the restricted isometry property (RIP) introduced in Candès and Tao (2005). This property determines a restricted isometry constant (RIC) δ_s for a certain sparsity level s, which is calculated by

$$\begin{array}{}\text{(10)}& (\mathrm{1}-{\mathit{\delta }}_s){∥\mathit{x}∥}_{\mathrm{2}}^{\mathrm{2}}\le {∥\mathbf{A}\mathit{x}∥}_{\mathrm{2}}^{\mathrm{2}}\le (\mathrm{1}+{\mathit{\delta }}_s){∥\mathit{x}∥}_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

where $\mathit{x}\in {\mathrm{\Sigma }}_s^n$. The RIC tells us how close the singular values of m×s submatrices of A are to 1. For δ_2s<1, any s-sparse solution can be uniquely determined by solving the Eq. (P1) problem, and for ${\mathit{\delta }}_{\mathrm{2}s}<\left(\sqrt{\mathrm{2}}-\mathrm{1}\right)$ this is even possible if the signal is superimposed by noise (noisy case) (Candès, 2008). In practice, calculating the RIC is NP-hard (see Tillmann and Pfetsch, 2014) and it is not applicable to calculate this constant for a given matrix. However, the RIP might be used within a design process, since there are known random distributions of matrices, which satisfy the RIP for large s considering sufficiently large n and m (Baraniuk et al., 2008). Another property, which very loosely upper-bounds the RIC, is the incoherence property (Wang et al., 2015). The coherence μ of a matrix is defined by

$$\begin{array}{}\text{(11)}& \mathit{\mu }=max\left|<{\mathit{a}}_i|{\mathit{a}}_j>\right|i\ne j,\end{array}$$

where a_i and a_j are distinct column vectors of $\stackrel{\mathrm{̃}}{\mathbf{A}}$, where $\stackrel{\mathrm{̃}}{\mathbf{A}}$ is the column-normalized matrix of A. The coherence bounds the RIC by

$$\begin{array}{}\text{(12)}& {\mathit{\delta }}_{\mathrm{2}s}\le (\mathrm{2}s-\mathrm{1})\mathit{\mu }.\end{array}$$

Therefore, s-sparse solutions can be uniquely recovered if $\mathit{\mu }<\frac{\mathrm{1}}{\mathrm{2}s-\mathrm{1}}$ holds in the noiseless case or $\mathit{\mu }<\frac{\sqrt{\mathrm{2}}-\mathrm{1}}{\mathrm{2}s-\mathrm{1}}$ in the noisy case.

In real-world scenarios, coefficients are rarely sparse but often compressible. This means that the coefficients can be well approximated by sparse coefficients. A more detailed explanation for compressible coefficients is found in Sect. 2.5. CS guarantees good estimates of compressible solutions if the RIP is satisfied for the noisy case. Then the reconstruction error can also be bounded (see Candès, 2008, for the exact definitions and bounds).

2.5 Sparsifying emissions

Sparsity is one of the key elements for SR. However, emissions are not always sparse. In order to make a non-sparse emission field sparse, a transformation into a different domain can be used. Such transformations include the Fourier transform, wavelet transform, transformations tailored to specific data sets, e.g., by SVD truncation (see Hong et al., 2011), or over-complete dictionaries, where the representation does not have to be orthogonal and multiple representations for the same emission field exist (see Candès et al., 2011). In this paper, we only deal with the discrete wavelet transform (DWT) to sparsify emissions. This transform is used for image compression (Lewis and Knowles, 1992) and was also used in Ray et al. (2014) and Ray et al. (2015) to parameterize fossil fuel CO₂ emissions for the US.

There are several ways to quantify the sparsity of an emission field. One possibility is to measure the error, using any l_p norm, of an emission field x to its best s-sparse approximation:

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_s(\mathit{x}{)}_p=inf\left\{{∥\mathit{x}-\mathit{z}∥}_p,\mathit{z}\in {\mathrm{\Sigma }}_s^n\right\}.\end{array}$$

Independent of the norm used, the best approximation is given by an emission field z, which contains the same s highest values of x and is zero otherwise. Often, it is more intuitive to give the relative sparsity s_rel, which is the fraction of non-zero entries as a percentage of all entries, instead of the sparsity level s. In this paper, both notations are used, e.g., σ_10%(x)₂ is the l₂ error of the signal which best approximates x and maximally possesses 10 % non-zero elements while σ₁₀(x)₂ is the l₂ error of the signal approximating x with maximally 10 non-zero elements. To show the distribution of values in x, a plot showing σ_s(x) over all possible s can be used.

Figure 1 Visualization of the compressibility of emission fields by plotting σ_s(x)₁ over all s. The more hyperbolic the curve is, the more compressible the emission field.

An example is depicted in Fig. 1. The more hyperbolic the plot is, the more compressible the signal is. In order to get a measure for the compressibility of the distribution, the Gini index is used, which was also proposed as a sparsity measure for CS by Zonoobi et al. (2011). This index measures how unequally x is distributed, which is a similar measure to the hyperbolicity of a curve in the depicted figure. For the case of an equal distribution (yellow curve) the Gini index becomes 0, and for a highly compressible distribution (orange curve) the Gini index gets close to 1. A Gini index of exactly 1 is equal to a 1-sparse signal.

In this paper, a DWT is used to sparsify emissions, using Haar wavelets. Throughout the paper, we use the third-level wavelet transform and refer to its matrix as W and its inverse transformation as W⁻¹. For an introduction to wavelets, we refer to Graps (1995) and Mallat (1999). For the visualization of a third-level DWT, see Fig. 1.1 in Mallat (1999).

The emission map x is transformed into its wavelet coefficients c by Wx=c and vice versa by $\mathit{x}={\mathbf{W}}^{-\mathrm{1}}\mathit{c}$. To transform the forward model Ax=y into the wavelet domain, we write

$$\begin{array}{}\text{(14)}& \mathit{y}=\mathbf{A}\mathit{x}={\mathbf{AW}}^{-\mathrm{1}}\mathit{c}=\stackrel{\mathrm{‾}}{\mathbf{A}}\mathit{c},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathbf{A}}={\mathbf{AW}}^{-\mathrm{1}}$ is the new forward model. The wavelet coefficients are then determined by the modified Eq. (P1) problem.

$$\begin{array}{}\text{(15)}& \widehat{\mathit{c}}=\underset{\mathit{c}}{\mathrm{argmin}}{∥\mathit{c}∥}_{\mathrm{1}}\text{ s.t. }\stackrel{\mathrm{‾}}{\mathbf{A}}\mathit{c}=\mathit{y}.\end{array}$$

Note that the CS conditions for $\stackrel{\mathrm{‾}}{\mathbf{A}}$ are not the same as for A. Once the wavelet coefficients $\widehat{\mathit{c}}$ are determined, the emissions are estimated by $\widehat{\mathit{x}}={\mathbf{W}}^{-\mathrm{1}}\widehat{\mathit{c}}$.

2.6 Reconstruction algorithms

There have been many algorithms proposed for the task of CS and SR. However, our goal in this paper is not to compare these algorithms, instead, we want to demonstrate the applicability of SR in general for urban GHG emission assessments. We, therefore, only solve the initial SR problem, given in Eq. (P1). To find the minimum, we use the cvx library, a matlab package for specifying and solving convex programs (Grant and Boyd, 2014, 2008), using the gurobi optimizer as a backend (Gurobi Optimization, LLC, 2021).

For the noisy case, we solve a modified form of Eq. (P1), given by

$$\begin{array}{}\text{(P1e)}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}{∥\mathit{x}∥}_{\mathrm{1}}\text{ s.t. }{∥\mathbf{A}\mathit{x}-\mathit{y}∥}_{\mathrm{2}}^{\mathrm{2}}\le {∥\mathit{ϵ}∥}_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

where ϵ is the noise vector. This is equivalent to the Lasso with the right choice of λ. Since we generate the noise by design, the optimization process simplifies by choosing ϵ without having to find the right λ value.

We compare our results to regularized least squares, which we refer to as least squares (LSs) hereafter. The equation of the LS is given by

$$\begin{array}{}\text{(P2)}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}{∥\mathit{x}∥}_{\mathrm{2}}\text{ s.t. }\mathbf{A}\mathit{x}=\mathit{y},\end{array}$$

in the noiseless case. There, the solution is given by the pseudo-inverse. In the noisy case, the equation is given by

$$\begin{array}{}\text{(P2e)}& {\displaystyle }\widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}{∥\mathit{x}∥}_{\mathrm{2}}\text{ s.t. }{∥\mathbf{A}\mathit{x}-\mathit{y}∥}_{\mathrm{2}}^{\mathrm{2}}\le {∥\mathit{ϵ}∥}_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

Table 1 Table of symbols and measures used in the paper.

2.7 Data evaluation

Table 1 gives an overview of the most important symbols and measures used in this paper. All of the measures we are using for evaluating the reconstruction results compare the estimations to the true emissions, which are the city emission inventories in this paper. If the l_p error is not specified, the l₂ error is used. The primary measure we use for the error evaluation of our estimates is the relative l₂ error, since it provides a measure for the highest spatial resolution of the emissions (1 km × 1 km), while the relative l₂ smoothed error shows an error for a lower spatial resolution (5 km × 5 km) and eliminates errors which are due to spatial errors. Because of the high spatial resolution of the emission fields reconstructed in this paper, relative l₂ error values greater than 1 are possible. Those are because of spatial errors, where an emission is assigned to the wrong spatial location. The relative total error disregards these spatial errors and only evaluates the difference between the total estimated emission and the true total emission of a city.

The mean relative error is used if the relative error of emission grid points independent of their contribution to the total amount of emissions is of interest. This metric is useful to evaluate how well a reconstruction method estimates emitters of a certain kind.

Table 2 Sparsity of the reported emission fields in different European cities. In the wavelet domain, in all of the cases a better or at least as good a sparse approximation of the emission map exists.

3 The sparsity/compressibility of emissions in European cities

In the following, we determine the sparsity and compressibility of real-world emissions. To do so, we study the CO₂ emissions of major European cities, such as Berlin, Hamburg, Munich, London, Paris, and Vienna using data from the TNO_GHGco_v1.1 emission inventory (van der Gon et al., 2019). This database provides annual gridded anthropogenic emissions with a resolution of about 1 km × 1 km for the areas and species. Furthermore, the emissions are classified by a proxy, such as public power, industry, and stationary combustion. Emissions of all proxies of the data set are included for the emission data used, making the assumption that this provides a good representation of typical emission fields for the cities investigated. We measure the sparsity of the cities using the Gini index and σ_10 %(x)₂ in the spatial and wavelet domain. The results are given in Table 2. For all the cities we consider, the wavelet domain of the emission fields achieves a higher Gini index compared to the spatial domain. Furthermore, the approximation error of s_rel=10 % of the signal is lower in the wavelet domain, except for Vienna, where this error is identical for both domains. From these data, we conclude that a representation of the city emission fields in the wavelet domain should generally be better suited for SR. We also split the cities into two groups: cities whose emission fields are highly compressible (Berlin, Hamburg, Vienna, and Munich; Gini index >0.7) and cities where the fields are significantly less compressible (London and Paris; Gini index <0.7).

4 Evaluating sparse reconstruction of GHG emissions

In the following, we define the estimation problem (Sect. 4.1) and apply SR to several European cities (Sect. 4.2).

4.1 Formulating the estimation problem

In the following, we assume that the influence of the background GHG concentration on the measurements is fully known. Subtracting this influence from the measurements yields the enhancement produced by the GHG emissions in the domain of interest. Therefore, the background can be ignored, and the model is simplified to a pure transport emission system. Let the sensitivity of the GHG concentration measurements y to the city emissions x be given by the sensing matrix A:

$$\begin{array}{}\text{(16)}& \mathit{y}=\mathbf{A}\mathit{x}.\end{array}$$

The rows in A contain vectorized footprints, which determine the sensitivity of the measurements to the GHG fluxes of different emission grid cells in the domain. These footprints are determined by backward transport models, such as the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al., 2003; Gerbig et al., 2003). In this paper, a simplified linear model, the Gaussian plume model, is used as the transport model. This approach allows varying parameters within the transport model without the computationally costly calculations of the STILT model. In Appendix B, we explain how the Gaussian plume footprints are created and show two footprints: one calculated by STILT, the other by the Gaussian plume model (Fig. B1). In our work, we use artificial wind data. These wind fields used, except in Sect. 4.2.1, possess a wind coverage of ∼143^∘. All of the footprints generated by those wind fields are available at https://doi.org/10.5281/zenodo.5901298 (Zanger et al., 2022a).

Figure 2 Visualization of a sensing matrix, where the values in the matrix are color encoded. The entry A_ji in the matrix gives the sensitivity of the jth measurement to the ith grid cell.

Figure 2 visualizes the sensing matrix A, which represents the sensitivity of all measurements to each grid cell. Depending on the footprints in A, there might be emissions in x which are not strongly sensed. The total sensitivity to a single grid cell x_i is determined by the sum of the corresponding column in the sensing matrix. If this sensitivity is below a certain threshold, ${\sum }_j{\mathbf{A}}_{ji}\le \mathit{\gamma }$, where $\mathit{\gamma }={\mathrm{10}}^{-\mathrm{9}}\mathrm{ppm}\cdot {\mathrm{km}}^{\mathrm{2}}\cdot \mathrm{h}\cdot {\mathrm{mol}}^{-\mathrm{1}}$ is the threshold, we remove the ith column from the sensing matrix and do not reconstruct this emission grid cell. For Fig. 2 this would be the first column of the sensing matrix, where all values are 0. This introduces some small perturbation between y and Ax, so that Ax≈y. However, since we choose our threshold small enough, this is not reflected in the solution but improves the conditioning of A and reduces the computational effort.

From a physical point of view, these removed emission grid cells are not situated upwind of the measurement stations and, therefore, cannot be well reconstructed using the measurements. For the wavelet domain, we perform this step on the wavelet transformed sensing matrix. Therefore, wavelet coefficients which are weakly sensed and below the threshold are removed.

Figure 3 Map of exemplary measurement locations for Munich using seven measurement stations, which represents a ratio $\frac{m}{n}$ of approximately 75 % (background: OpenStreetMap contributors, 2020).

In our study, we use both artificially created emission fields and real emission fields (van der Gon et al., 2019). Both field types have a spatial resolution of 1 km × 1 km. For the artificial emissions, a domain size of 32 km × 32 km is used, resulting in n=1024 unknown emission grids. For the emission inventory, the number of emission grids depends on the size of the city (see Table 2). In order to compare emission estimates of the different domains, we vary the number of measurement stations used for each domain, so that the number of measurements per unknown emission grid ($\frac{m}{n}$) is roughly constant. During the measurement period, each station takes 50 samples. The measurement stations were randomly distributed over a normalized area of the emission map and scaled to the domain size of each city. An exemplary distribution of the stations for Munich is depicted in Fig. 3, using $\frac{m}{n}\approx \mathrm{75}$ %.

The measurements y are then created using the sensing matrix A and the emissions x: y=Ax. In the noisy case, we add a Gaussian distributed noise vector to the measurements: $\mathit{y}=\mathbf{A}\mathit{x}+\mathit{ϵ}$, $\mathit{ϵ}\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}})$. σ_ϵ is chosen based on the SNR used.

Figure 4 Visualization of the (a) emission inventory of Munich and estimated emission fields of the inventory by using (b) LS, (c) SR, and (d) SR in the wavelet domain. The color scale is logarithmic and given in the unit of µmol (m² s)⁻¹. Negative emissions are colored dark blue and cannot be visually distinguished from low or 0 emitters. The emission inventory does not include negative emitters.

4.2 Applying sparse reconstruction to European cities

Reconstruction results for Munich in the noiseless case using LS, SR, and SR in the wavelet domain are depicted in Fig. 4, and the reconstruction errors for the figure can be found in Table 3. A visual inspection indicates that SR in the wavelet domain does resemble the real emissions best, while SR in the spatial domain is second best. Negative emitters are not visible because of the color scale. LS estimates a total of 171 such negative emitters with a total of 2139 µmol (m² s)⁻¹ negative emissions, while SR estimates only 33 negative emitters with a total of 9 µmol (m² s)⁻¹ negative emissions. SR in the wavelet domain estimates 50 negative emitters with a total of 59 µmol (m² s)⁻¹ negative emissions.

In the following, we identify properties of SR and compare them to the LS using the European city emission inventories. The properties we identify are interconnection between SR and CS (Sect. 4.2.1), the wavelet domain to make emissions sparser (Sect. 4.2.2), the ability of identifying unknown emitters (Sect. 4.2.3), and the number of measurements needed (Sect. 4.2.4). We first analyze these properties in the noiseless case (Sect. 4.2.1 to 4.2.4) and then extend those to the noisy case (Sect. 4.2.5).

Additionally, in Appendix A, we present an example using artificial, sparse emissions to better establish further connections between SR and CS.

4.2.1 Influence of wind coverage

We examine the effect of wind coverage on the effectiveness of the sensing matrix for SR. The term wind coverage is used to measure the range of wind directions during the measurement period. A wind coverage of 0^∘ corresponds to no changes in the wind direction during the observations, while a wind coverage of 360^∘ means that the wind blew from all directions during the time of observation.

Figure 5 Estimation errors for SR in the spatial and wavelet domain and LS reconstruction of (a) Munich's and (b) Paris's emissions with footprints covering different amount of wind directions. With a higher wind coverage, the estimation result using SR improves. The LS has no clear improvement for higher wind coverages, but its performance decreases with high coverages for Munich, while for Paris the performance first increases until about a coverage of 150^∘ before decreasing again.

The wind coverage is changed in an interval from 96 to 360^∘, with a step size of 24^∘, using a Gaussian plume model. As emission field we use the CO₂ city emission inventories of Munich and Paris from Sect. 3. Figure 5 depicts the relative errors of the SR in the spatial and wavelet domain and LS estimates of (a) Munich's and (b) Paris's emission fields. Errors of the estimates for both SR in the spatial and wavelet domain decrease with increasing wind coverage for Munich and Paris. In contrast, the LS estimates reveal no clear improvement for higher wind coverages, but the relative error of LS increases with high coverages for Munich, while for Paris the error initially decreases, until about a coverage of 150^∘, before increasing again. A sensitivity analysis, given in Appendix F, reveals that error variations given by the different wind coverages are due to the shift in sensitivity to different parts of the domain, while the sum of the sensitivities does not change.

This major improvement in SR for higher wind coverage can be explained using the incoherence property from CS. The coherence parameter, given in Eq. (11), calculates the maximum similarity of how two emission grid cells are measured. For the spatial domain, by increasing the wind coverage, there is greater variety in the measurements and thus between the measurements of two different grid cells. In the wavelet domain, we determined the coherences at different wind coverages and found a similar behavior.

Even though for both domains the coherence for a wind coverage of 360^∘ is too high to offer reconstruction guarantees, our example shows that including CS parameters in the design process can improve reconstruction results.

4.2.2 Sparse reconstruction in the wavelet domain

In the following, we compare reconstruction results of SR in the wavelet and spatial domains. Accordingly, s-sparse city emission fields are used to compare these domains. More precisely, the s highest emitters of the CO₂ inventory of London are used to generate an s-sparse emission field. This allows us to determine at which trade-off in sparsity the wavelet domain can achieve better results compared to the spatial domain. We chose the inventory of London, since it has the emission field with the lowest compressibility (see Table 2).

Figure 6 Comparison of the LS (yellow) and SR in the spatial (blue) and wavelet (red) domains in terms of relative sparsity s_rel. An emission map artificially generated from London's emission inventory with a total of 2205 unknown emitters was used. The x axis represents the portion of these emitters that were used (the emitters are sorted from large to small). In (a) 22 measurement stations with a total of 1100 measurements and (b) 33 measurement stations with a total of 1650 measurements are used.

Figure 6 shows the relative reconstruction error in the spatial and third-level wavelet domain for different sparsity levels using $\frac{m}{n}\approx \mathrm{50}$ % and $\frac{m}{n}\approx \mathrm{75}$ %. As expected, SR in the spatial domain performs better for sparse emission maps (low s). For less sparse emission maps (higher s), however, the wavelet domain is superior.

With fewer measurements, the error in the wavelet domain already plateaus for s_rel≈15 %, while for $\frac{m}{n}\approx \mathrm{75}$ % the error plateaus at s_rel≈40 %. The results also indicate that for the entire emission map of London (s_rel=1), the performance of SR in the wavelet domain does not increase significantly when the measurements are changed from $\frac{m}{n}\approx \mathrm{50}\mathit{\%}$ to $\frac{m}{n}\approx \mathrm{75}$ %, while the performance in the spatial domain benefits considerably. This improved performance is caused by the higher compressibility of the emission map in the wavelet domain, where the highest wavelet coefficients already provide a good representation of the entire emission map, while the spatial domain needs more coefficients for a representation with the same relative error. Therefore, the wavelet domain requires fewer measurements compared to the spatial domain, to achieve the same performance for cities with emission fields that have a low compressibility (for similar results see Sect. 4.3.3).

Our result demonstrates that the wavelet domain can help to improve SR in those instances for which the spatial domain is not well suited. However, changing the domain also alters the conditions of CS; therefore, conclusions made for CS in the spatial domain cannot be directly transferred to the wavelet domain.

Figure 7 Qualitative and quantitative measure of how well SR and LS estimate the highest emissions for (a) Munich and (b) Paris.

4.2.3 Discovering unknown emitters

In the following, we evaluate how well SR can identify unknown emitters. Assuming a good prior estimation of the emission field, where the relative error for the real emissions for each emitter is approximately constant and small, unknown emitters can be assumed to make a huge contribution to the difference between the prior expected measurements and real measurements. Therefore, we evaluate how well the highest emitters are reconstructed to assess the performance of finding unknown emitters.

To do this, we consider the emission inventory data for Munich and Paris. While Munich's emissions inventory is highly compressible, the emission inventory of Paris is not. For the setup, seven measurement stations are used in Munich and 39 in Paris. This results in $\frac{m}{n}\approx \mathrm{75}$ % for both cities, making the estimation results comparable.

We employ a qualitative and a quantitative measure, both of which are depicted in Fig. 7. The x axes show bins of emission grids ranked by their emission strength. For example, the first bin (0 %–0.4 %) contains the first 0.4 % highest-emission grids, while the last bin in the lower panels of Fig. 7 contains emission grids which are among the 36.2 %–53.2 % highest-emission grids. Since the emissions have a high compressibility (the first highest emissions make a huge contribution to the total emissions), we use a logarithmic scale for the bins.

The qualitative measure (see upper panels in Fig. 7) compares the fraction of the highest emitters in the inventory matching those of the estimation (the higher the ratio, the better). In Munich, all three methods reconstruct the same largest 0.4 % emitters. For Paris, SR in the spatial domain performs best followed by the LS in the case of the largest 0.4 % emitters. For lower emitters, SR in the wavelet domain performs best for both cities. While the performance gain of SR and SR DWT over the LS method is quite large in Munich, which has a highly compressible emission inventory, the gain is smaller in Paris, which has a slightly compressible inventory.

The quantitative measure (see lower panels in Fig. 7) compares the relative error in the reconstructed emissions (the lower, the better). For the largest 0.4 % emitters, SR in the spatial domain performs best and can reconstruct these emissions with an error of less than 0.7 % in Munich and 3 % in Paris. When analyzing slightly lower emissions, SR in the spatial domain performance drops significantly and the wavelet domain performs best in both cities.

The results indicate that SR in the spatial domain works particularly well for the highest emitters (largest 0.4 % in our examples), while SR in the wavelet domain performs well for a broader range of high emitters. The LS performance is much less sensitive to the emission strength of individual strong emitters, which is not beneficial for estimating large unknown emitters.

These results for the same scenario but with fewer measurement stations ($\frac{m}{n}\approx \mathrm{50}$ %) can be found in Appendix D.

4.2.4 Number of measurements needed

Previous scenarios have presented comparisons between SR and LS using $\frac{m}{n}\approx \mathrm{50}$ % or $\frac{m}{n}\approx \mathrm{75}$ %. However, the question of how many measurements m are needed to produce results of a certain accuracy remains unanswered. To determine this, we vary the number of measurements per station, which leads to a variation in $\frac{m}{n}$.

Figure 8 Relative error for SR, SR with DWT, and LS at different degrees of undersampling, for (a) Munich and (b) Paris.

Figure 8 depicts the relative estimation error for Munich and Paris for varying $\frac{m}{n}$. In Munich, where the emissions are highly compressible, the reconstruction error for SR in the spatial domain already drops sharply for $\frac{m}{n}\gtrsim \mathrm{10}$ %. SR in the wavelet domain needs more measurements but performs similar to the spatial domain using $\frac{m}{n}\gtrsim \mathrm{23}$ % and even performs slightly better than the spatial domain for $\frac{m}{n}>\mathrm{50}$ %. The LS also sees a steep improvement with $\frac{m}{n}\approx \mathrm{30}$ % but does not produce results as good as SR.

For Paris, where the emissions have a low compressibility, SR in the spatial domain performs worst for $\frac{m}{n}<\mathrm{20}$ %, while the wavelet domain performs about the same as the LS for this number of measurements. For a higher number of measurements ($\frac{m}{n}>\mathrm{0.2}$), both SR methods show increased performance and both perform better than the LS. Overall, SR in the wavelet domain performs best.

These trends are also confirmed for the emissions of other cities (see Supplement). These results demonstrate that when using SR instead of the LS, a higher undersampling with fewer errors is possible, especially for highly compressible emissions (as in Munich).

4.2.5 Measurement noise

To make the results applicable to real-world scenarios, noise also has to be taken into consideration. There are different types of noise, including measurement noise, transport error, and representation error. In the following, we assume measurement errors with an SNR typical of column measurements. Chen et al. (2016) report measurement errors for column measurements for an average time of 10 min for CO₂ as 0.04 to 0.05 ppm and for CH₄ as 0.2 ppb. Jones et al. (2021) show that CH₄ concentration enhancements measured for the city of Indianapolis are in the order of several parts per billion. Below, we assume an SNR of 20 dB, which corresponds to enhancements of 2 ppb for methane or 0.4 to 0.5 ppm for CO₂.

Figure 9 Reconstruction errors of SR in the spatial and wavelet domain as well as LS for Munich when varying the SNR and $\frac{m}{n}\approx \mathrm{1.5}$.

As emission fields, we use the inventory data for Munich. Compared to the noiseless case, we increase the number of measurement stations to $\frac{m}{n}\approx \mathrm{1.5}$, which is 15 measurement stations with a total of 750 measurements. Figure 9 depicts how the relative l₂ error changes for different SNRs in the Munich simulation. The figure shows that the results we obtain for an SNR of 20 dB should not differ significantly for slightly lower or higher SNRs between 15 and 25 dB. Therefore, we believe that despite our noise assumptions, which are optimistic compared to real-world scenarios, our results are qualitatively valid. We provide a variance map for SR with an SNR of 20 dB in Appendix E.

Figure 10 Reconstruction errors use SR in the spatial and wavelet domain as well as the LS for Munich with an SNR of 20 dB, where (a) shows the average relative error for the reconstruction of 250 different noise vectors at different levels of resolutions using $m/n\approx \mathrm{1.5}$ and (b) shows the relative error for the highest resolution (1 km × 1 km) when varying the number of measurements.

Here, we show the relative error of the estimated emission fields for the noisy case at different spatial resolutions of the emissions, both for SR in the spatial domain and wavelet domain and for the LS (see Fig. 10a). For this case we have run the reconstruction of 250 different noise vectors and show the mean estimation result. Furthermore, we assess the undersampling capability in the noisy case (see Fig. 10b). For SR in the spatial domain, the relative error is only slightly higher for high resolutions and remains fairly constant across other spatial resolutions, suggesting that almost all errors occur in the reconstruction of emission strength rather than in the localization of emitters. Both SR in the wavelet domain and the LS have a high relative error for high resolutions (below 7 km × 7 km), while at lower resolutions both perform better than SR in the spatial domain, indicating that most of the error is due to incorrect localization of emitters.

As a result, they do not provide accurate localization of the emitters. To overcome this, more measurements are needed for SR in the wavelet domain and LS. To determine the number of measurements required, we show the relative errors for the highest resolution (1 km × 1 km) when varying the number of measurements in Fig. 10b. For SR in the spatial domain, the relative error drops significantly for $\frac{m}{n}<\mathrm{50}$ % and performs better than the other reconstruction methods for more measurements. The wavelet domain performs similarly to the LS for $\frac{m}{n}<\mathrm{1}$ and better for $\frac{m}{n}>\mathrm{1}$. These results indicate that SR in the spatial domain, for Munich, is more robust against noise, while the wavelet domain is not as robust (many more measurements are needed than in the noiseless case for similar results).

Figure 11 Qualitative and quantitative measure of how well SR and LS estimate the highest emissions with a noise of 20 dB for Munich. The qualitative plot shows how many percent of the highest emissions in the inventory are also contained in the same highest amount of emission of the reconstruction. The quantitative plot shows the mean relative errors for different emission strengths.

Next, we analyze the performance of finding unknown emitters in the noisy case. For this purpose, we use the same setup as in Sect. 4.2.3 but solely look at the emissions in Munich. We again use $\frac{m}{n}\approx \mathrm{1.5}$ and add the noise with an SNR of 20 dB. We employ a qualitative and a quantitative measure, both of which are depicted in Fig. 11 (for the description of the measures, see Sect. 4.2.3). SR is the only method to identify the 0.4 % highest emitters (qualitative measure) and estimates their emissions with a significantly lower relative error of 2 % (quantitative measure). Slightly smaller emitters are better estimated by the SR in the wavelet domain. For even lower emitters, there is no clear winner of which method estimates them best. Compared to the noiseless case, the estimation errors are in general larger, and SR in the wavelet domain does not achieve as promising results as in the noiseless case. SR in the spatial domain remains sensitive to the largest emitters in the noisy case. These results for Paris can be found in the Supplement.

5 Broader comparison of sparse reconstruction for European cities

In the following, case studies on emission fields of all cities considered in Sect. 3 are introduced to more broadly compare the performance of SR. The measurement configuration corresponds to the specifications in Sect. 4.1. The emissions are reconstructed using SR in the spatial and wavelet domains as well as the LS, first with $\frac{m}{n}\approx \mathrm{75}$ % in the noiseless case and then with $\frac{m}{n}\approx \mathrm{1.5}$ in the noisy case.

Table 3 Reconstruction performances for SR in the spatial and wavelet domain as well as for LS by measuring the relative error (1 km × 1 km), relative smoothed error (5 km × 5 km), and the relative total error for different European cities. The best approach for each city and type of error is highlighted in bold.

Table 4 Reconstruction performances for SR in the spatial and wavelet domain as well as for LS by using $\frac{m}{n}\approx \mathrm{1.5}$ with an SNR of 20.0 dB for different European cities. The results are evaluated using the relative error (1 km × 1 km), relative smoothed error (5 km × 5 km), and the relative total error. The best approach for each city and type of error is highlighted in bold.

To compare the performance, we measure relative errors, relative smoothed errors (for a spatial resolution of 5 km × 5 km), and relative total errors. The results for the noiseless case are in Table 3. At the highest spatial resolution (1 km × 1 km, right column in Table 3), SR performs much better than LS for the highly compressible emission fields in both domains. For the cities with emission fields that are slightly compressible, the SR performance is significantly worse than for the other cities but still slightly better than LS. For a lower spatial resolution of 5 km × 5 km (middle column in Table 3), the error of all methods decreases. Among the highly compressible emissions, SR still performs significantly better than LS. Both for London and Paris (slightly compressible), all methods produce a similar error, with the LS giving the smallest error in London and SR in the spatial domain giving the smallest error in Paris. For the total emissions, all reconstruction methods perform similarly.

These results support our previous findings that SR performs well at high resolutions for cities with highly compressible emissions.

Next, we consider the noisy case with an SNR of 20.0 dB and $\frac{m}{n}\approx \mathrm{1.5}$ (see the results in Table 4).

For the cities that are slightly compressible, SR in the spatial domain performs the worst, while SR in the wavelet domain and the LS perform similarly.

For the highly compressible cities, SR in the spatial domain produces the best results for the 1 km × 1 km and 5 km × 5 km resolution and is worse than the other reconstruction methods for the total emissions (except for Vienna). Furthermore, SR in the wavelet domain performs better than the LS. In contrast to the noiseless case, SR in the wavelet domain performs worse than SR in the spatial domain for high resolutions (≲5 km × 5 km). Both results support our findings from Sect. 4.2.5 that SR in the spatial domain accurately localizes emissions and that SR in the wavelet domain is not as robust to noise as the spatial domain for high resolutions.

6 Conclusions

In this paper, we introduced sparse reconstruction (SR) as a novel method for the inversion of urban GHG emissions and further provided key examples to identify the advantages of SR to inversely model emissions.

SR can be easily integrated into existing top–down frameworks for estimating emission. We examined the applicability of this method by evaluating the sparseness of emissions from several European cities. Our results indicate that the emissions from most of these cities are sparse and that SR is applicable. We also showed that a wavelet transform increases the sparsity of urban emissions, making SR applicable to cities with less sparse emissions. SR is known to have reconstruction guarantees if conditions of compressed sensing (CS) are satisfied. We tested different CS conditions for various wind fields and showed that wind fields with better CS conditions do significantly increase the performance of SR.

Compared to state-of-the-art inversion methods using Gaussian priors, our method requires fewer measurements and provides better localization and quantification of unknown emitters.

SR works best if the underlying representation of the emissions is sparse. In this paper, we employed the wavelet transform to increase the sparsity; however, other transformations, such as a curvelet transform or more general dictionary representations, might be even more suited for specific spatial domains. Finding such transformations which also work well with CS conditions is challenging, and future studies should be devoted to them.

Appendix A: Boundaries of the emission inventories

Table A1 shows the longitude and latitude boundaries used to create the CO₂ emission fields for the European cities from the TNO_GHGco_v1.1 emission inventory.

Table A1 Longitude and latitude boundaries used to create the CO₂ emission inventories for the European cities from the TNO_GHGco_v1.1 emission inventory.

Appendix B: Wind fields and generation of the sensing matrix

The sensing matrices we use consist of vectorized footprints generated by the Gaussian plume model with artificial wind data. The Gaussian plume model provides a good enough approximation for Lagrangian particle dispersion models for the purpose of this paper. The footprints are generated for a high-resolution grid using a dynamic Gaussian plume with varying wind speed, direction, and diffusion. The footprints are then scaled to the domain size of the city. Because of the scaling, the wind speed and diffusion changes for each city. Nevertheless, we think that this approach makes the reconstruction performances for different cities more comparable than using different footprints for every city.

For Sect. 4.2.1, we use a static Gaussian plume model to generate footprints for different wind coverages. This allows us to ignore additional effects from the dynamical case, for example, in the dynamic case the change in the wind direction produces spiral-looking footprints. These footprints have a different shape compared to the footprints generated using a static Gaussian plume model. By using a static model, the footprints are more predictable and systematic, which makes the results of different wind coverages more comparable to each other.

Figure B1 shows a comparison between a STILT footprint for Munich and a Gaussian plume footprint used in our work, where the two footprints are generated using different wind fields.

Figure B1 Comparison of a STILT footprint to a Gaussian plume footprint, where different wind fields have been used.

Appendix C: Recovery of sparse emissions

Consider s emitters producing an s-sparse emission field in a city. We distribute these emitters over the city using a Gaussian distribution, such that the probability for an emitter in the center of the city is higher compared to outside of the city. Furthermore, we use a Gaussian distribution to model the emission strength of each emitter with a mean value of 10 µmol (m² s)⁻¹ and a standard deviation of 1 µmol (m² s)⁻¹, so that all emitters have roughly a similar contribution to the total emissions. In Fig. C1, the reconstruction errors for the spatial domain using SR and LS are depicted. SR allows exact recovery for s_rel≲13 % non-zero emissions using $\frac{m}{n}\approx \mathrm{50}\mathit{\%}$ and for s_rel>35 % non-zero emissions using $\frac{m}{n}\approx \mathrm{75}\mathit{\%}$. In contrast, the relative error of LS is generally much higher but decreases with a larger number of measurements.

Figure C1 Reconstruction errors for sparse cities with (a) 50 % and (b) 75 % measurements m per emission fields n. Using SR (blue line), it is possible to reconstruct sparse emissions. By increasing the number of measurement stations, and therefore observations, more emitters can be reconstructed using SR (a up to 13 %; b up to 35 % non-zero emission cells). For the LS method (red line), the relative error decreases both for less sparse signals and with a larger number of observations.

Figure C2 Visualization of falsely discovered (dark blue), undiscovered (yellow), and discovered (cyan) emissions in x by (a) SR and (b) LS.

These results demonstrate the power of SR for sparse emission fields and illustrate the link between SR and CS. As shown in Sect. 2.4, the sensing matrix A must satisfy some sufficient conditions to ensure that CS yields the best possible reconstruction. However, in real-world scenarios, it is challenging to verify these conditions. By using a Monte Carlo algorithm, we found counterexamples that showed δ_2s>1 for already small s (both for Gaussian plume footprints and STILT footprints). Thus, the sensing matrix we use here does not guarantee CS for all cases. Nevertheless, A can be useful if it satisfies the CS conditions for some specific emission distributions of interest.

Next, we consider the case with $\frac{m}{n}\approx \mathrm{75}\mathit{\%}$ and choose a number of emitters s, such that the relative error produced by SR and the LS is similar.

Therefore, both methods produce the same relative error, which allows us to compare the spatial distribution of the error on an equal basis. This is the case for s_rel=46.9 %, where the relative error for SR is 28.1 %, while the relative error for LS is 27.3 %. We focus on the following questions:

• which emitters are found correctly (discovered)?

• which emitters are not found (undiscovered)?

• which emitters are discovered, even though there is no emitter there (falsely discovered)?

Figure C2 depicts the discovered, not discovered, and falsely discovered emitters for SR (panel a) and LS (panel b). SR reconstructs most emitters correctly, only 1.1 % of the emitters are undiscovered (yellow), and 8.0 % are falsely discovered (dark blue). In contrast, the LS has no emitters which have not been found but 45.6 % falsely discovered emitters. This is because the LS prefers smooth estimates, which is equivalent to estimating emitters everywhere. This is also seen by the fraction between the emissions estimated for falsely discovered emitters and the total estimated emissions. While this fractions is only 3.0 % for SR, for LS this fraction is 12.2 %.

These results demonstrates the advantage of SR for the localization of emitters. This property is especially useful to find unknown emitters, since the spatial certainty of emitters for SR is much higher.

Appendix D: Discovering unknown emitters using fewer measurements

In Sect. 4.2.3 we showed that SR is better at finding high emitters compared to LS and can reconstruct them with smaller errors. This is visualized by Fig. 7, which has been created using $\frac{m}{n}\approx \mathrm{75}\mathit{\%}$. In the following, we visualize this result using only $\frac{m}{n}\approx \mathrm{50}\mathit{\%}$.

Figure D1 Qualitative and quantitative measure of how well SR and LS estimate the highest emitters for Munich (a) and Paris (b) using $\frac{m}{n}\approx \mathrm{50}\mathit{\%}$.

The results are depicted in Fig. D1. Compared to the case with more measurements, for Munich, SR in the spatial domain performs better in the qualitative measure compared to the wavelet domain. Furthermore, for Paris, the margin between the performance of SR in the wavelet domain and the other methods has increased.

Appendix E: Uncertainty quantification

Uncertainty quantification for sparse reconstruction is distinctly different compared to uncertainty quantification for LS reconstruction. First, while for the LS fit with Gaussian noise a closed-form solution exists, there does not exist any closed-form solution for sparse reconstruction.

In the following we use the uncertainty quantification approach from Hase et al. (2017) in order to derive a variance plot of the reconstruction result by applying bootstrapping. This gives us the covariance matrices for SR in the spatial and wavelet domain. We only use the variances of the matrix and plot the standard deviations for every single reconstructed emission, which is depicted in Fig. E1. For the spatial domain, the standard deviation is localized around certain areas in the map, while for the wavelet domain, the standard deviation covers much larger areas of the map. This lets one believe that a reconstructed emission in the spatial domain has a higher certainty of location compared to SR in the wavelet domain. This assumption is verified by Fig. 10.

Figure E1 Uncertainty of the emission estimates for Munich with $\frac{m}{n}\approx \mathrm{1.5}$ and an SNR of 20 dB for SR in the (a) spatial domain and (b) wavelet domain.

Appendix F: Least square sensitivity for different wind coverages

In the following, we show that the changes in relative error for LS in Sect. 4.2.1 are due to the fact that the sensitivity to different areas of the map changes, while the total sensitivity does not change. Figure F1 depicts the sensitivity map of Munich for a wind coverage of (a) 72^∘ and (b) 360^∘. For 72^∘, the southeast has good sensitivity to the estimation, while other parts, especially the west of the city, are weakly sensitive. For 360^∘, there is no area in the city which is not sensitive to the estimation. However, some areas are less sensitive than in the 72^∘ case, especially the southeast. If we compare this sensitivity map to the emission map of Munich in Fig. 4, we see that the strongest emitter is positioned in the southeast. This explains why, for the case of 72^∘, the relative error of the LS is lower than for the 360^∘ case. If we look at the mean sensitivity of an emission grid cell for both cases, we find that those are exactly the same (0.6629), and a change in wind coverage does not improve the total sensitivity.

Figure F1 Sensitivity map of LS in Munich for a wind coverage of (a) 72^∘ and (b) 360^∘. In (a), the southeast has good sensitivity to the estimation, while the west and north of the city are only weakly sensitive. In (b), there is no area which is not sensitive to the estimation, but the areas in the southeast are less sensitive.