Observing geometry effects on a 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 1 Least Squares (LSQ) and a Compressive Sensing (CS) approach for water vapor tomography is compared based on synthetic 2 Global Navigation Satellite System (GNSS) Slant Wet Delay (SWD) estimates. In this context, the term observing geometry 3 mainly refers to the number of GNSS sites situated within a specific study area subdivided into a certain number of volumetric 4 pixels (voxels) and to the number of signal directions available at each GNSS site. The novelties of this research are 1) the 5 comparison of the observing geometry’s effects on the tomographic reconstruction accuracy when using LSQ resp. CS for the 6 solution of the tomographic system and 2) the investigation of the effect of the signal directions’ variability on the tomographic 7 reconstruction. The tomographic reconstruction is performed based on synthetic SWD data sets generated, for many samples 8 of various observing geometry settings, based on wet refractivity information from the Weather Research and Forecasting 9 (WRF) model. The validation of the achieved results focuses on a comparison of the refractivity estimates with the input WRF 10 refractivities. The results show that the recommendation of Champollion et al. (2004) to discretize the analyzed study area into 11 voxels with horizontal sizes comparable to the mean GNSS inter site distance represents a good rule of thumb for both LSQ 12 and CS based tomography solutions. In addition, this research shows that CS needs a variety of at least 15 signal directions 13 per site in order to estimate the refractivity field more accurately and more precisely than LSQ. Therefore, the use of CS is 14 particularly recommended for water vapor tomography applications for which a high number of multi-GNSS SWD estimates 15 are available. 16

In addition to slant wet delay estimates from GNSS, Hurter and Maier (2013) introduce wet refractivity profiles from radio 30 occultation and radiosonde observations into a combined least squares collocation. Rather than using slant wet delay estimates 31 as input observations, Nilsson and Elgered (2007) apply a solution that relies directly on GPS phase observations. 32 Independently of the reconstruction strategy, due to the point-wise GNSS observing geometry, the tomographic system of 33 tomography, to the determination of SWD estimates -fulfills certain prerequisites. For general applications of CS, Rauhut 23 (2010) e.g. states that randomness in the acquisition step helps to utilizing the minimum number of measurements. When 24 reconstructing images based on frequency data, e.g. Candès et al. (2006) alternatively recommend to randomly measure fre-25 quency coefficients such that sparse objects are sensed by taking as few measurements as possible. For CS-based water vapor 26 tomography approaches, no explicit requirements for the SWD acquisition or for designing advantageous observing geometry 27 settings have been established so far. 28 For LSQ, Champollion et al. (2004) state that the optimal horizontal size of a voxel should correspond to the mean inter-site 29 distance between the used GNSS sites. Given a certain cutoff elevation angle, the height layers' thicknesses in their approach 30 should be defined such that signals received at a GNSS site situated within a voxel's center are able to cross neighboring voxels. 31 Due to the small wet refractivity values in the upper layers and in order to make the tomographic solutions less sensitive to 32 errors in the input data, Rohm (2012) recommends to increase the height layer thicknesses with increasing altitude. Bender and 33 Raabe (2007) resp. Bender et al. (2009) estimate the spatial distribution of the geometric intersection points between different 34 ray paths in order to compute the information density contained in a given set of GPS signals. They then use this information 35 as a precondition to an optimal tomographic reconstruction resp. in order to identify regions that are well covered by GPS 1 slant paths. Although Bender et al. (2011b) and Zhao et al. (2019) state that changing the observing geometry by combining 2 multi-GNSS observations instead of GPS only observations does not substantially improve the reconstruction quality, Rohm 3 (2012) realizes that the uncertainty of the tomographic solution is largely influenced by the mathematical properties of the 4 design matrix, depending itself on the observing geometry. With the aim of giving advice for the installation of new permanent 5 sites and for the solution of future water vapor tomographies, this work therefore investigates the observing geometry's effect 6 on the quality of both a LSQ and a CS solution to the tomographic system. In order to analyze the observing geometry's effect on the quality of the LSQ and CS solution to water vapor tomography, 9 different observing geometry settings are defined. Based on synthetic SWD estimates derived from WRF, 3D water vapor 10 distributions are reconstructed for each of the defined observing geometry settings using both LSQ and CS. The quality of the 11 LSQ and CS solutions to water vapor tomography is then compared w.r.t. the respective observing geometry settings.

15
where SWD i, cont stands for the integrated slant wet delay observations between a certain satellite and a certain GNSS site and 16 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 17 ray path, i.e. the slant ray path of the radiowave signal between a certain satellite and a certain receiver. The variable N wet 18 contains the wet refractivity along this path. The index i attains the values 20 where N corresponds to the number of observations available between any receiver and any satellite. When tomographically 21 reconstructing the wet refractivity, however, the continuous functional model from Equation 1 is usually replaced by a discrete 22 functional model That is, the 3D water vapor distribution is discretized into L = P × Q × K voxels in longitude, latitude, and height, assuming a 25 constant refractivity value for each voxel. As in Heublein et al. (2018) and Heublein (2019), in this work, a uniform voxel dis-26 cretization is selected in the horizontal directions, while the voxel sizes in the vertical direction increase with increasing height.

27
Horizontally, the voxels are limited by constant longitudes and latitudes. In the vertical direction, the voxels are separated by 28 layers of constant ellipsoidal height.
in a parameter vector x ∈ R L×1 , and all distances d ij in a design matrix Φ data ∈ R N ×L , the discrete tomographic system from 4 Equation 3 can be rewritten as 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 9 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 2 large study 10 area subdivided into 5×5×5 voxels, disposing of seven GNSS sites and ten signals per site). For voxels that are not crossed by 11 any signals, Φ data has a zero column. Therefore, the tomographic model and the mathematical properties of the design matrix 12 largely depend on the observing geometry settings described in Section 3.3. 13 3.2 Solution of the inverse tomographic model using LSQ resp. CS

14
The LSQ solution to Equation 5 is derived by solving the minimization problem regularization constraints and prior knowledge 16 regularized by means of t = 3 regularization terms, namely by horizontal and vertical smoothing constraints as well as by prior knowledge from surface meteorology. Namely, the wet refractivity at the surface is computed from in-situ observations 18 of pressure, temperature, and dew point temperature at a single weather site within study area. As described in Heublein

23
The weights can e.g. be derived using inverse distance weighting 2 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.

3
Moreover, Davis et al. (1993) state that an average refractivity profile can be approximated assuming the refractivity to expo-4 nentially decrease with height: The variable h k is the height of the kth layer, h 0 stands for some reference height at which the refractivity equals N wet (h 0 ), and 7 H scale represents the scale height of the local troposphere. As H scale is essential for defining an exponential decay with height,     The iii) number of signal directions per site is motivated by the GPS resp. by a multi-GNSS orbit geometry. According to   The horizontal distribution of the synthetic GNSS sites within the URG study area is shown in Figure 2. The signal direc-3 tions result from selecting at random the defined number of signal directions from a synthetic multi-GNSS orbit constellation 4 composed of GPS, GLONASS, and Galileo. Both signals entering the study area on its top and on its side are included. 5 From WRF, simulations of water vapor mixing ratio, temperature, pressure, and geopotential height are available at a 900 m   Figure 2. Distribution of the seven GNSS permanent sites (blue squares) as well as of the five 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. the CS refractivity estimates and the WRF refractivity of the considered voxel approaches zero for most samples. However, 1 e.g. for 27 sites and 20 signal directions per site, there are some samples in which the CS based refractivity estimate differs 2 from the WRF refractivity by up to 3.3 ppm. I.e. for many signal directions, CS is able to accurately and precisely reconstruct 3 the voxel's refractivity, but for some signal directions, the voxel's refractivity estimate does not match well with the voxel's 4 validation refractivity from WRF. In contrast, in the case of few sites and few signal directions per site (e.g. twelve sites and 5 ten signal directions per site), LSQ yields refractivity estimates differing from −5.9 ppm to −0.7 ppm from the WRF refractiv-6 ity, while the CS refractivity estimates differ much more from the WRF refractivity (differences of −42.9 ppm to 26.9 ppm). 7 Consequently, when investigating the observing geometry's effect on the quality of the tomographic reconstruction, the chosen 8 solution strategy as well as the effect of varying signal directions absolutely need to be taken into account. Therefore, in this 9 research, a representative set of 48 half-hourly samples of synthetic GNSS orbits is considered in order to analyze the observing 10 geometry's effect on the tomographic reconstruction quality. As expected, an increased number of sites and an increased number of signal directions per site, in general, decrease the 1 mean of the absolute difference (called mean difference in the following) and the standard deviation of the difference between 2 estimated refractivities and WRF refractivities. Yet, as shown in Figure 4 averaged over all voxels, in the case of a LSQ solution 3 to the tomographic system, the mean difference decrease by means of introducing more SWD estimates into the tomographic 4 reconstruction is much smaller than that in the case of a CS solution. When averaged over 48 samples per observing geometry, 5 introducing more SWD estimates improves the mean difference by up to 1.3 ppm resp. 1.9 ppm (maximum improvement 6 observed for 20 resp. 15 signal directions per site in the case of LSQ resp. CS).  The legend in the upper left subplot holds for all subplots. In each subplot, the improvement by introducing 32 sites instead of seven sites is given in red resp. blue. Degradations are given with a minus sign.

8
When investigating the standard deviation of the differences between estimated refractivities and WRF refractivities for the CS   Moreover, as the presented approach only relies on a synthetic data set deduced from WRF, the synthetic SWDs introduced 30 within the tomographic system in this research are too optimistic, when compared to real GNSS SWD estimates. Therefore, the 31 conclusions drawn in Section 4 cannot necessarily be transferred to tomographic applications involving real SWD estimates. 32 1 of adding different types of noise to the synthetic SWD estimates should be investigated (e.g. measurement and sensor noise 2 and uncertainties resulting from the observing geometry). In the presented approach, instead of mapping ZWDs to the slant 3 signal directions as in the case of a real GNSS processing, the synthetic SWD data set is computed based on a direct raytracing 4 within the same voxels in which the tomographic reconstruction is thereafter performed. Yet, Heublein (2019) shows that this 5 involves neglecting both a voxel discretization error and a mapping error committed in the case of real data.
6 Furthermore, Section 4 shows that the standard deviation of the difference between LSQ refractivity estimates and WRF re-7 fractivities is 6 % to 65 % smaller than that computed based on the CS refractivity estimates, if at most ten different signal 8 directions per site are available. In contrast, in the case of a high number of sites and a high number of signal directions per 9 site, the LSQ reconstruction is not able to yield as accurate and as precise estimates of the water vapor distribution as CS.   Rauhut, H.: Compressive sensing and structured random matrices, Theoretical foundations and numerical methods for sparse recovery, 9,