Observing geometry effects on a Global Navigation Satellite System (GNSS)-based water vapor tomography solved by least squares and by compressive sensing

In this work, the effect of the observing geometry on the tomographic reconstruction quality of both a regularized least squares (LSQ) approach and a compressive sensing (CS) approach for water vapor tomography is compared based on synthetic Global Navigation Satellite System (GNSS) slant wet delay (SWD) estimates. In this context, the term “observing geometry” mainly refers to the number of GNSS sites situated within a specific study area subdivided into a certain number of volumetric pixels (voxels) and to the number of signal directions available at each GNSS site. The novelties of this research are (1) the comparison of the observing geometry's effects on the tomographic reconstruction accuracy when using LSQ or CS for the solution of the tomographic system and (2) the investigation of the effect of the signal directions' variability on the tomographic reconstruction. The tomographic reconstruction is performed based on synthetic SWD data sets generated, for many samples of various observing geometry settings, based on wet refractivity information from the Weather Research and Forecasting (WRF) model. The validation of the achieved results focuses on a comparison of the refractivity estimates with the input WRF refractivities. The results show that the recommendation of Champollion et al. (2004) to discretize the analyzed study area into voxels with horizontal sizes comparable to the mean GNSS intersite distance represents a good rule of thumb for both LSQ- and CS-based tomography solutions. In addition, this research shows that CS needs a variety of at least 15 signal directions per site in order to estimate the refractivity field more accurately and more precisely than LSQ. Therefore, the use of CS is particularly recommended for water vapor tomography applications for which a high number of multi-GNSS SWD estimates are available.

1 Introduction

In this paper, we intend to determine the three-dimensional (3-D) atmospheric water vapor distribution for each point in time. This adds further essential information to the spatio-temporal analyses of two-dimensional (2-D) water vapor fields commonly used in weather forecasting and climate research. In addition, atmospheric water vapor delays the microwave signal propagation within the atmosphere and thus represents an error source in, e.g., Global Navigation Satellite System (GNSS) and Interferometric Synthetic Aperture Radar (InSAR) observations. Therefore, a precise knowledge of the water vapor field, for example, is required for accurate deformation monitoring using InSAR (Hansen and Yu, 2001). However, the atmospheric water vapor distribution is difficult to model because it is highly variable in time and space. Several approaches exist for reconstructing the 3-D tomographic water vapor reconstruction using one-dimensional (1-D) GNSS slant wet delays (SWDs); see Sect. 2.

One of the main limiting factors in water vapor tomographies consists of the point-wise GNSS observing geometry, which causes an ill-conditioned inverse tomographic model that needs to be regularized. Yet, even after regularization, the observing geometry composed, e.g., of the number and the geographic distribution of the GNSS sites, the SWD signal directions, and the voxel discretization still affect the quality of the tomographic solution. This work therefore meets the challenge of comparing the observing geometry's effect on a GNSS-based water vapor tomography solved by means of least squares (LSQ) or by means of compressive sensing (CS). By investigating the observing geometry's effect on the LSQ and CS solution strategies, the differences between the LSQ solution and a CS solution approach benefiting from the signal's sparsity in an appropriate transform domain for regularization are better understood, and recommendations can be given for future water vapor tomography campaigns and the processing of their measurements. Based on synthetic data sets deduced from the Weather Research and Forecasting (WRF) model described in Skamarock et al. (2005), the presented work answers the research question of (1) to what extent the rule of thumb of Champollion et al. (2004), derived for LSQ and recommending a voxel size corresponding to the mean GNSS intersite distance, can be transferred from an LSQ solution to a CS solution. In addition, this research investigates (2) in which settings CS is able to more accurately and more precisely reconstruct the tomographic water vapor field than LSQ and (3) to what extent multi-GNSS SWD observations improve the tomographic solution obtained by means of LSQ or CS when compared to solutions obtained from SWDs originating from the Global Positioning System (GPS) only.

2 Related work

Current water vapor tomographies can be distinguished, e.g., based on the methodology and the data sets applied for solving the tomographic model. The tomographic model is commonly established based on the directions along which space-geodetic SWD estimates are acquired and based on a discretization of the investigated atmospheric volume into volumetric pixels (voxels), e.g., of constant refractivity. The existing tomography solution approaches applied to such a discretized atmosphere are subdivided into iterative and non-iterative techniques. Bender et al. (2011a) propose different iterative algebraic reconstruction techniques (ARTs), while Hirahara (2000), Flores et al. (2000), Champollion et al. (2004), Troller (2004), Song et al. (2006), Notarpietro et al. (2008), and Rohm (2013) apply different non-iterative methods for solving the tomographic system using an LSQ adjustment. Thanks to its good capability to estimate dynamically changing parameters, Flores et al. (2000), Gradinarsky and Jarlemark (2004), and Rohm et al. (2014) choose a Kalman filter approach. Hirahara (2000) proposes a damped least squares solution known from seismic tomography to solve the tomographic problem. Xia et al. (2013) combine iterative and non-iterative techniques. Instead of using voxels in water vapor tomographies with a small number of GNSS sites, Ding et al. (2018) discretize the tomographic field based on the perimeter of the tomographic boundary on the plane and based on meshing techniques. They then determine tomographic fields by means of fitting the real distribution of GNSS signals on different tomographic planes at different tomographic epochs.

In addition to slant wet delay estimates from the GNSS, Hurter and Maier (2013) introduce wet refractivity profiles from radio occultation and radiosonde observations into a combined least squares collocation. Rather than using slant wet delay estimates as input observations, Nilsson and Elgered (2007) apply a solution that relies directly on GPS phase observations.

Independently of the reconstruction strategy, due to the point-wise GNSS observing geometry, the tomographic system of equations is usually ill posed and needs to be regularized, e.g., (i) by constraining the tomographic system by means of pseudo observations, (ii) by introducing additional observations from models, from simulations, or from other sensors, or (iii) by decreasing the number of voxels crossed by no rays at all.

Both Flores et al. (2000) and Gradinarsky and Jarlemark (2004) regularize the solution by means of adding horizontal and vertical smoothing constraints to the tomographic system and by means of introducing a boundary constraint assuming the refractivity to approach zero above a certain height. Alternatively, as proposed in Elosegui et al. (1998), the refractivity field can be assumed to decrease exponentially with increasing height. Yet, while regularizing the solution significantly, both geometric constraints and the exponential decay usually are not able to accurately model the real atmospheric state. As an alternative to using horizontal and vertical constraints relying on physical approximations to the atmospheric behavior, the work in Heublein et al. (2018) and Heublein (2019) exploits the signal's sparsity in a particular, predefined transform domain as prior knowledge for regularizing the tomographic system in order to then reconstruct the signal using an L₁-norm minimization. Similarly, CS and sparse reconstruction are applied here for the tomographic reconstruction of the 3-D water vapor field, and the CS solution to water vapor tomography is compared to a solution obtained using a classical LSQ approach. Initially presented by Candès (2006), Donoho (2006), Baraniuk (2007), and Candès and Wakin (2008) for the image or signal recovery from a number of samples below the desired resolution or the Nyquist rate, CS has been, since then, applied to many remote sensing problems in which sparse signals occur. For example, Potter et al. (2010) and Alonso et al. (2010) describe the use of CS for synthetic aperture radar (SAR) imaging, Pruente (2010) applies CS for ground-moving target identification, Zhu and Bamler (2010), Budillon et al. (2011), Aguilera et al. (2013), and Zhu and Bamler (2014) apply CS to SAR tomography, and Li and Yang (2011), Zhu and Bamler (2013), Grohnfeldt et al. (2013), Jiang et al. (2014), and Zhu et al. (2016) use CS for pan-sharpening and hyperspectral image enhancement. When compared to classical LSQ adjustments usually applying L₂-norm regularizations, compressive sensing and sparse reconstruction based on a small number of measurements led to promising results.

However, compressive sensing only yields encouraging results if the input data acquisition – corresponding, in water vapor tomography, to the determination of SWD estimates – fulfills certain prerequisites. For general applications of CS, Rauhut (2010), e.g., states that randomness in the acquisition step helps in utilizing the minimum number of measurements. When reconstructing images based on frequency data, Candès (2006), e.g., alternatively recommend randomly measuring frequency coefficients such that sparse objects are sensed by taking as few measurements as possible. For CS-based water vapor tomography approaches, no explicit requirements for acquiring the SWD or for designing advantageous observing geometry settings have been established so far.

For LSQ, Champollion et al. (2004) state that the optimal horizontal size of a voxel should correspond to the mean intersite distance between the used GNSS sites. Given a certain cutoff elevation angle, the height layers' thicknesses in their approach should be defined such that signals received at a GNSS site situated within a voxel's center are able to cross neighboring voxels. Due to the small wet refractivity values in the upper layers and in order to make the tomographic solutions less sensitive to errors in the input data, Rohm (2012) recommends increasing the height layer thicknesses with increasing altitude. Bender and Raabe (2007) and Bender et al. (2009) estimate the spatial distribution of the geometric intersection points between different ray paths in order to compute the information density contained in a given set of GPS signals. They then use this information as a precondition for an optimal tomographic reconstruction and in order to identify regions that are well covered by GPS slant paths. Although Bender et al. (2011b) and Zhao et al. (2019) state that changing the observing geometry by combining multi-GNSS observations instead of GPS-only observations does not substantially improve the reconstruction quality, Rohm (2012) realizes that the uncertainty of the tomographic solution is largely influenced by the mathematical properties of the design matrix, depending itself on the observing geometry. With the aim of giving advice for the installation of new permanent sites and for the solution of future water vapor tomographies, this work therefore investigates the observing geometry's effect on the quality of both an LSQ solution and a CS solution to the tomographic system.

3 Methodology

In order to analyze the observing geometry's effect on the quality of the LSQ and CS solutions to water vapor tomography, different observing geometry settings are defined. Based on synthetic SWD estimates derived from WRF, 3-D water vapor distributions are reconstructed for each of the defined observing geometry settings using both LSQ and CS. The quality of the LSQ and CS solutions to water vapor tomography is then compared with respect to the respective observing geometry settings.

3.1 Tomographic model

For tomography using GNSS SWDs, Flores et al. (2000) introduce the functional model

$$\begin{array}{}\text{(1)}& {\mathrm{SWD}}_{i,\mathrm{cont}}={\mathrm{10}}^{-\mathrm{6}}\cdot \underset{{\mathrm{sp}}_i}{\int }{N}_{\mathrm{wet}}\text{d}l,\end{array}$$

where SWD_i, cont stands for the integrated slant wet delay observations between a certain satellite and a certain GNSS site and dl is a differential along the slant ray path. As in Heublein et al. (2018) and Heublein (2019), the variable sp_i is the ith slant ray path, i.e., the slant ray path of the radio-wave signal between a certain satellite and a certain receiver. The variable N_wet contains the wet refractivity along this path. The index i attains the values

$$\begin{array}{}\text{(2)}& i\in \mathbb{N}\mathrm{with}\mathrm{1}\le i\le N,\end{array}$$

where N corresponds to the number of observations available between any receiver and any satellite. When tomographically reconstructing the wet refractivity, however, the continuous functional model from Eq. (1) is usually replaced by a discrete functional model:

$$\begin{array}{}\text{(3)}& {\mathrm{SWD}}_{i,\mathrm{disc}}={\mathrm{10}}^{-\mathrm{6}}\cdot \sum _{j=\mathrm{1}}^L{N}_{\mathrm{wet}j}\cdot d_{ij}.\end{array}$$

That is, the 3-D water vapor distribution is discretized into $L=\mathcal{P}\times \mathcal{Q}\times \mathcal{K}$ voxels in longitude, latitude, and height, assuming a constant refractivity value for each voxel. As in Heublein et al. (2018) and Heublein (2019), in this work, a uniform voxel discretization is selected in the horizontal directions, while the voxel sizes in the vertical direction increase with increasing height. Horizontally, the voxels are limited by constant longitudes and latitudes. In the vertical direction, the voxels are separated by layers of constant ellipsoidal height.

As in Heublein et al. (2018) and Heublein (2019), summarizing all observations SWD_i, disc in an observation vector ${\mathit{y}}_{\mathrm{data}}\in {\mathbb{R}}^{N\times \mathrm{1}}$, all unknown refractivities N_wet j for

$$\begin{array}{}\text{(4)}& j\in \mathbb{N}\mathrm{with}\mathrm{1}\le j\le L\end{array}$$

in a parameter vector $\mathit{x}\in {\mathbb{R}}^{L\times \mathrm{1}}$, and all distances d_ij in a design matrix ${\mathbf{\Phi }}_{\mathrm{data}}\in {\mathbb{R}}^{N\times L}$, the discrete tomographic system from Eq. (3) can be rewritten as

$$\begin{array}{}\text{(5)}& {\mathit{y}}_{\mathrm{data}}={\mathbf{\Phi }}_{\mathrm{data}}\cdot \mathit{x},\end{array}$$

with

$$\begin{array}{}\text{(6)}& {\mathrm{\Phi }}_{\mathrm{data}}(i,j)=\mathit{\{}\begin{array}{cl}{d}_{ij}& \mathrm{if}\mathrm{signal}i\mathrm{crosses}\mathrm{voxel}j,\\ \mathrm{0}& \mathrm{otherwise}.\end{array}\end{array}$$

As each signal only passes a small subsection of the study area, most entries of the matrix Φ_data are zero and only a few matrix elements are non-zero (e.g., only about 5 % of the entries of Φ_data are non-zero in the case of an about 95×99 km² large study area subdivided into $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ voxels, consisting of seven GNSS sites and 10 signals per site). For voxels that are not crossed by any signals, Φ_data has a zero column. Therefore, the tomographic model and the mathematical properties of the design matrix largely depend on the observing geometry settings described in Sect. 3.3.

3.2 Solution of the inverse tomographic model using LSQ or CS

The LSQ solution to Eq. (5) is derived by solving the minimization problem

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \widehat{\mathit{x}}=\underset{\mathit{x}}{\mathrm{argmin}}\mathit{\{}\underset{\mathrm{data}\mathrm{fidelity}\mathrm{term}}{\underbrace{{∥{\mathit{y}}_{\mathrm{data}}-{\mathbf{\Phi }}_{\mathrm{data}}\cdot \mathit{x}∥}_{\mathrm{2}}^{\mathrm{2}}}}\\ & +\underset{\mathrm{regularization}\mathrm{constraints}\mathrm{and}\mathrm{prior}\mathrm{knowledge}}{\underbrace{\sum _{t=\mathrm{1}}^{\mathrm{3}}{\mathrm{\Gamma }}_{{\mathrm{constraints}}_t}^{\mathrm{2}}\cdot {∥{\mathit{y}}_{{\mathrm{constraints}}_t}-{\mathbf{\Phi }}_{{\mathrm{constraints}}_t}\cdot \mathit{x}∥}_{\mathrm{2}}^{\mathrm{2}}}}\mathit{\}},\end{array}\end{array}$$

regularized by means of t=3 regularization terms, namely by horizontal and vertical smoothing constraints as well as by prior knowledge from surface meteorology. That is, the wet refractivity at the surface is computed from in situ observations of pressure, temperature, and dew point temperature at a single weather site within the study area. As described in Heublein et al. (2018), the horizontal smoothing constraints assuming the refractivity of a voxel $(a,b,k)$ to equal the weighted mean refractivity of the surrounding voxels $(p,q,k)$ with voxel indices p≠a and q≠b within the same height layer k are defined by

$$\begin{array}{}\text{(8)}& N_{{\mathrm{wet}}_{a,b,k}}=\sum _{p,q}{w}_{p-a,q-b}\cdot N_{{\mathrm{wet}}_{p,q,k}}.\end{array}$$

The weights can, e.g., be derived using inverse distance weighting:

$$\begin{array}{}\text{(9)}& w_{p-a,q-b}=\left\{\begin{array}{ll}\frac{\frac{\mathrm{1}}{d_{p-a,q-b}}}{{\sum }_{p,q}\frac{\mathrm{1}}{d_{p-a,q-b}}}& \mathrm{if}(a,b)\ne (p,q),\\ -\mathrm{1}& \mathrm{if}(a,b)=(p,q),\end{array}\right.\end{array}$$

with distances $d_{p-a,q-b}$ between the center of voxel (p,q) and the center of voxel (a,b) of the considered kth height layer.

Moreover, Davis et al. (1993) state that an average refractivity profile can be approximated assuming the refractivity to exponentially decrease with height:

$$\begin{array}{}\text{(10)}& N_{\mathrm{wet}}\left(h_k\right)=N_{\mathrm{wet}}\left(h_{\mathrm{0}}\right)\cdot \exp \left(-{\displaystyle \frac{h_k-h_{\mathrm{0}}}{H_{\mathrm{scale}}}}\right).\end{array}$$

The variable h_k is the height of the kth layer, h₀ stands for some reference height at which the refractivity equals N_wet(h₀), and H_scale represents the scale height of the local troposphere. As H_scale is essential for defining an exponential decay with height, its value is determined within the solution of the tomographic system from a set of realistic values for H_scale between 1000 and 2000 m. Both the weights for the horizontal and vertical smoothing constraints and for the prior knowledge from surface meteorology are determined with respect to the data fidelity term using the place holder trade-off parameter ${\mathrm{\Gamma }}_{\mathrm{constraint}}^{\mathrm{2}}$ in Eq. (7). The selection of the trade-off parameters from a certain number of logarithmically scaled possible trade-off parameters and the selection of H_scale are described in Heublein et al. (2018) and Heublein (2019).

When aiming at a tomographic reconstruction of atmospheric water vapor by means of compressive sensing, the parameters x are sparsely represented in some transform domain

$$\begin{array}{}\text{(11)}& \mathit{x}=\mathbf{\Psi }\cdot \mathit{s}\end{array}$$

as sparse parameters s. Estimates $\widehat{\mathit{s}}$ for these sparse parameters are obtained by

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \widehat{\mathit{s}}=\underset{\mathit{s}}{\mathrm{argmin}}\mathit{\{}\underset{\mathrm{data}\mathrm{fidelity}\mathrm{term}}{\underbrace{{∥\mathit{y}-\mathbf{\Phi }\cdot \mathbf{\Psi }\cdot \mathit{s}∥}_{\mathrm{2}}^{\mathrm{2}}}}+\underset{{\mathrm{L}}_{\mathrm{1}}\text{-norm regularization}}{\underbrace{{\mathrm{\Gamma }}_{\mathrm{CS}}^{\mathrm{2}}\cdot {∥\mathit{s}∥}_{\mathrm{1}}}}\\ & +\underset{\mathrm{prior}\mathrm{knowledge}\mathrm{from}\mathrm{surface}\mathrm{meteorology}}{\underbrace{{\mathrm{\Gamma }}_{\mathrm{constraints}}^{\mathrm{2}}\cdot {∥{\mathit{y}}_{\mathrm{constraints}}-{\mathbf{\Phi }}_{\mathrm{constraints}}\cdot \mathbf{\Psi }\cdot \mathit{s}∥}_{\mathrm{2}}^{\mathrm{2}}}}\mathit{\}},\end{array}\end{array}$$

as described in Heublein et al. (2018) and Heublein (2019). Instead of adding horizontal and vertical constraints to the data fidelity term as in Eq. (7), an L₁-norm regularization term is introduced into the CS solution to promote sparse solutions for s, as described in Heublein et al. (2018) and Heublein (2019). The L₁-norm minimization of the sparse parameters reduces the solution space. The wet refractivity estimates $\widehat{\mathit{x}}$ are then reconstructed using

$$\begin{array}{}\text{(13)}& \widehat{\mathit{x}}=\mathbf{\Psi }\cdot \widehat{\mathit{s}},\end{array}$$

with a dictionary $\mathbf{\Psi }\in {\mathbb{R}}^{L\times M}$. As in Heublein et al. (2018) and Heublein (2019), the dimension M of the parameters $\mathit{s}\in {\mathbb{R}}^{M\times \mathrm{1}}$ in the transform domain varies with the number of base functions or atoms in Ψ. Similarly to Heublein et al. (2018) and Heublein (2019), we assert that a sparse representation of the refractivity distribution can be obtained by means of, e.g., a dictionary composed of Kronecker products of inverse discrete cosine transform (iDCT) letters in longitudinal and latitudinal directions and of Euler letters and Dirac letters in the height direction. When thinking of languages, an atom would stand for a word within a dictionary. Comparably to a word composed of different letters within a language dictionary, each atom results from the Kronecker product of smaller items, namely letters, within the dictionary for sparse representation. That is, for each column, i.e., line, of the 3-D wet refractivity signal, a Kronecker product of the 1-D letters along the longitudinal, latitudinal, and height directions is computed. In the longitudinal and latitudinal directions, iDCT letters will represent horizontal refractivity variations. The Euler letters model an exponential refractivity decay with height, and the Dirac letters describe deviations from a decay described by a linear combination of Euler letters.

3.3 Characteristics of the study area and observing geometry settings

In Sect. 4, tomographic solutions obtained based on a high number of different observing geometry settings are compared. The observing geometry settings result from (i) a fixed voxel discretization, (ii) 7 to 32 sites, (iii) 5 to 20 signal directions per site, and (iv) 48 signal direction samples per number of sites and signals. Champollion et al. (2004) recommend (i) horizontal voxel sizes for an LSQ solution to water vapor tomography greater than or equal to the mean intersite distance between the available GNSS sites, i.e., voxel sizes greater than or equal to about 37×37 km² or to about 17×17 km² in the case of 7 or 32 uniformly distributed GNSS sites within the investigated study area of about 95×99 km² in size. In this work, the study area is discretized into $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ voxels of horizontal sizes of about 19×20 km². In the vertical direction, five height layers are distinguished. With increasing height, the height layer thicknesses increase from 1300 up to 2900 m. The lowest layer's thickness is set to 1300 m in order to ensure at least for signals with very low elevation angles that a signal arriving at the center of a voxel is able to pass the horizontally neighboring voxel within the same height layer. This is only possible if the minimum thickness Δdl_min of the height layers is related to the horizontal voxel size Δhz=20 km and to the cutoff elevation angle ϵ_cut=7^∘ by means of

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }\mathrm{d}{l}_{\min }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \mathrm{\Delta }hz\cdot \tan {\mathit{ϵ}}_{\mathrm{cut}}.\end{array}$$

The (ii) minimum number of seven sites originates from the real GNSS Upper Rhine Graben (URG) network site distribution within the analyzed study area. The maximum number of sites is chosen such that the rule of thumb of Champollion et al. (2004), introduced in Sect. 2, is clearly fulfilled. The horizontal position of the synthetic GNSS sites corresponds, for seven sites, to the position of real GNSS sites within the analyzed study area. The horizontal position of the additionally defined synthetic GNSS sites is chosen such that they are uniformly distributed within the study area. The vertical position of the synthetic GNSS sites corresponds to the height of the WRF digital elevation model at the horizontal position of the sites. The (iii) number of signal directions per site is motivated by the GPS or by a multi-GNSS orbit geometry. According to Feairheller and Clark (2006), the Global'naya Navigatsionnaya Sputnikova Sistema (GLONASS) constellation consists of 21 active plus 3 spare satellites on three orbital planes inclined by 64.8^∘ with respect to the Equator. In contrast, Hofmann-Wellenhof et al. (2008) describe the GPS orbit constellation as consisting of 21 satellites plus 3 spares on six orbital planes inclined by 55^∘ with respect to the Equator and the Galileo orbit constellation as consisting of 27 operational plus 3 spare satellites on three planes inclined by 56^∘ with respect to the Equator. Therefore, five signal directions and five visible satellites correspond to a pessimistic GPS setting, e.g., with site-specific shadowing. Eight signal directions per site are here considered to be a typical GPS setting at a GNSS permanent site (there may even be up to 12 observations for GPS-only), at which about 30 % of the total number of GPS satellites is visible at a time. Similarly, a total number of 20 signal directions per site corresponds to a visibility of about 30 % of the total number of satellites at a time within a multi-GNSS constellation composed of GPS, GLONASS, and Galileo. Similarly, observation of the BeiDou Navigation Satellite System could also be used in order to increase the number of signal directions in tomography applications. For each of the mentioned numbers of sites and numbers of visible satellites, 48 signal direction samples are defined. Given the GPS repeat cycle of about 1 d, the number of 48 signal direction samples is chosen in order to emulate about half-hourly orbit samples. Using the orbit characteristics mentioned above, such synthetic half-hourly satellite positions are approximated by means of circular orbits.

For each of the described observing geometry settings, synthetic SWD observations as input for the tomographic system are deduced from one single WRF simulation covering an about 200 km × 200 km large area centered around the longitude and latitude (λ, φ) = (8.15, 49.15^∘) at a 900 m horizontal resolution. The vertical resolution increases with height, ranging from about 50 to about 500 m. As schematically illustrated in Fig. 1, for each synthetic GNSS site, this is done by means of averaging the refractivity information of all WRF cells situated within the defined tomographic voxels, a direct ray tracing within these tomographic voxels, and adding together the SWD along each signal direction within the tomographic voxels using Eq. (3). As the voxels are limited by ellipsoidal upper and lower borders, the ray tracing is performed according to Perler (2011).

Figure 1 Schematic illustration of the generation of synthetic GNSS SWDs according to Heublein (2019): within each tomographic voxel (grey), the WRF refractivities of all those WRF cells (blue) situated within that voxel are averaged. A direct ray tracing along the considered signal direction then yields the SWD introduced into the tomographic system.

The horizontal distribution of the synthetic GNSS sites within the URG study area is shown in Fig. 2. The signal directions result from selecting at random the defined number of signal directions from a synthetic multi-GNSS orbit constellation composed of GPS, GLONASS, and Galileo. Both signals entering the study area at its top and on its side are included.

Figure 2 Distribution of the 7 GNSS permanent sites (blue squares) as well as of the 5 to 25 additional, synthetic sites (black symbols) within the URG study area. The additional, synthetic sites are distributed within a grid that uniformly covers the study area. Triangles, pentagons, hexagons, diamonds, and circles represent the first, second, third, fourth, and fifth group of five additional sites each.

From WRF, simulations of water vapor mixing ratio, temperature, pressure, and geopotential height are available at a 900 m horizontal resolution for generating the synthetic GNSS SWDs within the 95×99 km² large study area situated in the URG as shown in Fig. 2. The topography within the Rhine Valley is flat. Height differences mainly occur at the foot of the Black Forest mountain range. The height difference between the highest and lowest synthetic GNSS sites used for this study is about 494 m.

4 Results

For the most humid acquisition date (27 June 2005) for which WRF simulations were provided for this research and for an exemplary voxel in the lower middle of the lowest voxel layer, Fig. 3 shows that variations in the SWD signal directions available within the tomographic system cause variations in the estimated refractivities. The magnitudes of the absolute wet refractivity values during the analyzed weather conditions range from 0 to 74 ppm. As variations in the signal directions imply a change in the observed atmospheric volume, these variations in the estimated refractivities seem obvious. Yet, Fig. 3 illustrates that the variations in the refractivity estimates vary with the selected solution strategy. When considering many sites and many signal directions per site (e.g., at least 27 sites and at least 15 signal directions), the difference between the CS refractivity estimates and the WRF refractivity of the considered voxel approaches zero for most samples. However, e.g., for 27 sites and 20 signal directions per site, there are some samples in which the CS-based refractivity estimate differs from the WRF refractivity by up to 3.3 ppm. That is, for many signal directions, CS is able to accurately and precisely reconstruct the voxel's refractivity, but for some signal directions, the voxel's refractivity estimate does not match well with the voxel's validation refractivity from WRF. In contrast, in the case of few sites and few signal directions per site (e.g., 12 sites and 10 signal directions per site), LSQ yields refractivity estimates differing from −5.9 to −0.7 ppm from the WRF refractivity, while the CS refractivity estimates differ much more from the WRF refractivity (differences of −42.9 to 26.9 ppm). Consequently, when investigating the observing geometry's effect on the quality of the tomographic reconstruction, the chosen solution strategy as well as the effect of varying signal directions should be taken into account. Therefore, in this research, a representative set of 48 half-hourly samples of synthetic GNSS orbits is considered in order to analyze the observing geometry's effect on the tomographic reconstruction quality.

Figure 3 Absolute differences between estimated refractivity and the WRF refractivity in ppm for an exemplary voxel in the lower middle of the lowest voxel layer for the 48 samples of each investigated observing geometry setting. The two left columns use an ordinate ranging from −25 to 25 ppm, the third column plots the differences within the range −10 to 10 ppm, and the right column plots the differences within the range −5 to 5 ppm. The legend in the upper left subplot holds for all the subplots: red circles stand for the LSQ results, while blue squares represent the CS results. In each subplot, the minimum and maximum absolute differences in ppm of the LSQ or CS refractivity estimate with respect to the WRF input refractivities and the mean and the standard deviation over all samples are given in red and blue. Moreover, the mean and the standard deviation over all the samples are indicated for LSQ by a red dashed line, for CS by a blue dash-dotted line, and in the corresponding colors by error bars.

As expected, an increased number of sites and an increased number of signal directions per site, in general, decrease the mean of the absolute difference (called mean difference in the following) and the standard deviation of the difference between estimated refractivities and WRF refractivities. Yet, as shown in Fig. 4 averaged over all voxels, in the case of an LSQ solution to the tomographic system, the mean difference decrease by means of introducing more SWD estimates into the tomographic reconstruction is much smaller than that in the case of a CS solution. When averaged over 48 samples per observing geometry, introducing more SWD estimates improves the mean difference by up to 1.3 ppm or 1.9 ppm (maximum improvement observed for 20 or 15 signal directions per site in the case of LSQ or CS).

Figure 4 Averaged over all voxels, this figure shows the mean of the absolute difference and the standard deviation (SD) of the difference between estimated refractivities and WRF refractivities in ppm, deduced from 48 samples of each investigated observing geometry setting composed of a certain number of synthetic GNSS sites and various numbers of signal directions per site. The dashed and dotted lines serve for better following the variation of the represented quantities with the number of sites, but only the discrete values indicated by the markers should be evaluated. The legend in the upper left subplot holds for all the subplots. In each subplot, the improvement by introducing 32 sites instead of 7 sites is given in red and blue. Degradations are given with a minus sign.

A more in-depth analysis does not show any significant differences among the individual voxels. Moreover, no significant systematics between the LSQ and CS solutions appeared, neither at the boundary layer nor at the top of the atmosphere. For most of the tested scenarios, LSQ provides better results than CS. However, if there are many sites and many different signal directions available, CS yields more accurate and more precise results than LSQ.

When investigating the standard deviation of the differences between estimated refractivities and WRF refractivities for the CS case, considering an increased number of synthetic GNSS sites only while keeping a constant number of five signal directions per site is not advantageous. However, as of 15 different signal directions per site, a clear improvement in standard deviation is visible when increasing the number of sites in the tomographic setting solved by means of CS. Independently of the number of sites, for realistic GPS-like observing geometry settings with 5 to 10 signal directions per site, the LSQ refractivity estimates are more precise than the CS refractivity estimates. In contrast, as of 15 signal directions per site, the CS solution yields more accurate and more precise refractivity estimates than the LSQ solution if at least 22 sites are available. That is, this study shows that LSQ is less sensitive to the number of signal directions than CS. Therefore, we recommend using LSQ for water vapor tomography with GPS-only observations and CS for water vapor tomography with multi-GNSS SWD estimates.

In the case of the maximum number of sites and the maximum number of signal directions per site (32 sites and 20 signal directions per site), when averaged over the 48 considered samples per observing geometry, the mean difference and the standard deviation of the LSQ or CS reconstruction attain values of about 0.3 ppm or 0.0 ppm. Therefore, the number of sites and the number of signal directions per site are of particular interest when aiming at a very accurate and very precise tomographic reconstruction using CS.

For the given $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ voxel discretization with horizontal voxel sizes of 19 km and 20 km, the rule of thumb of Champollion et al. (2004) requires the mean intersite distance to correspond to no more than 19 to 20 km. The results show that, using the investigated synthetic data set with at least 15 signal directions per site, the CS solution is able to more accurately and more precisely reconstruct the atmospheric water vapor distribution than LSQ in the case of 22 sites within the 95×99 km² large study area, i.e., at a site density of about 1 site per 20.7×20.7 km², which is a bit lower than that required by the rule of thumb of Champollion et al. (2004). That is, if 15 signal directions are available per site, the rule of thumb can be transferred from the LSQ solution to water vapor tomography to CS solutions.

Consequently, the following three main results are summarized from this study.

1. The rule of thumb of Champollion et al. (2004) can be transferred from LSQ to CS.

2. Based on site distributions obeying the rule of thumb of Champollion et al. (2004), CS needs a variety of at least 15 signal directions per site in order to estimate the 3-D refractivity field more accurately and precisely than LSQ.

3. LSQ seems to be less sensitive to the number of signal directions than CS. Therefore, CS should only be used in the case of multi-GNSS SWD estimates yielding a variety of at least 15 signal directions per site.

5 Discussion and outlook

Section 4 states not only that the rule of thumb of Champollion et al. (2004) holds for a tomographic solution based on LSQ, but also that it ensures a good tomographic reconstruction in the case of a CS solution. Although this finding is based on many different observing geometry settings, it only refers to a single voxel discretization and to a single study area with a single topography and a single site distribution within that study area. As a consequence, this research mainly investigates the validity of the rule of thumb of Champollion et al. (2004) for CS for the given study area, weather condition, and voxel discretization. For generalization, further tests should be performed that repeat the described methodology for other study areas, weather conditions, and voxel discretizations and for site distributions varying not only in the number but also the position of the sites. In particular, using only five vertical layers potentially represents limits to the validity of the results of the proposed approach.

Moreover, as the presented approach only relies on a synthetic data set deduced from WRF, the synthetic SWDs introduced within the tomographic system in this research are too optimistic when compared to real GNSS SWD estimates. Therefore, the conclusions drawn in Sect. 4 cannot necessarily be transferred to tomographic applications involving real SWD estimates. In order to get a better idea of the transferability of the results, the analysis should be repeated based on real data or the effect of adding different types of noise to the synthetic SWD estimates should be investigated (e.g., measurement and sensor noise and uncertainties resulting from the observing geometry). In the presented approach, instead of mapping ZWDs (zenith wet delays, related to SWDs by means of mapping as described, e.g., in Niell, 1996) to the slant signal directions as in the case of a real GNSS processing, the synthetic SWD data set is computed based on a direct ray tracing within the same voxels in which the tomographic reconstruction is thereafter performed. Yet, Heublein (2019) shows that this involves neglecting both a voxel discretization error and a mapping error committed in the case of real data.

Furthermore, Sect. 4 shows that the standard deviation of the difference between LSQ refractivity estimates and WRF refractivities is 6 % to 65 % smaller than that computed based on the CS refractivity estimates if at most 10 different signal directions per site are available. In contrast, in the case of a high number of sites and a high number of signal directions per site, the LSQ reconstruction is not able to yield as accurate and as precise estimates of the water vapor distribution as CS. That is, when solving the tomographic system by means of LSQ, increasing the number of SWD signal directions improves the tomographic reconstruction quality less than when using CS. This may be due to the geometric smoothing constraints forming the basis of the LSQ solution. In the case of a small number of observations, the smoothing constraints ensure a smooth solution free of outliers that does not necessarily correspond to the prevailing atmospheric conditions. In the case of a high variety of observations, the smoothing constraints become less important with respect to the data fidelity term within the LSQ solution to the tomographic system, but they still affect the tomographic solution. Even in the case of a very high number of observations, the tomographic system cannot be solved in a purely data-driven way. Instead, the tomographic solution always takes into account the chosen model assumptions; i.e., the LSQ solution always applies a certain amount of smoothing.

In addition, a low or high number of signal directions chosen from a synthetic multi-GNSS constellation for the recommendation of LSQ or CS for GPS-only or multi-GNSS water vapor tomography applications should not be set equal to considering a real GPS-only or real multi-GNSS setting. Choosing a small number of signal directions from a multi-GNSS constellation yields a higher variability in the signal directions than choosing the same small number of signal directions from a GPS-only constellation. Since a high number of signal directions proved to be of particular importance in the case of a CS solution, the quality of the refractivity estimates deduced using CS may decrease if real GPS-only signal directions are chosen.

Finally, future research should analyze in more detail which signal directions are necessary in an LSQ- or CS-based water vapor tomography in order to reconstruct well the refractivities of as many voxels as possible. A two-step CS LSQ may then first yield accurate refractivity estimates for most voxels by means of CS and then use the geometric smoothing constraints applied in LSQ to improve the refractivity estimates of those voxels in which CS yields inaccurate refractivity estimates even if a high number of sites and a high number of signal directions per site are available.