the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Atmospheric bending effects in GNSS tomography
Gregor Möller
Daniel Landskron
In Global Navigation Satellite System (GNSS) tomography, precise information about the tropospheric water vapor distribution is derived from integral measurements like groundbased GNSS slant wet delays (SWDs). Therefore, the functional relation between observations and unknowns, i.e., the signal paths through the atmosphere, have to be accurately known for each station–satellite pair involved. For GNSS signals observed above a 15^{∘} elevation angle, the signal path is well approximated by a straight line. However, since electromagnetic waves are prone to atmospheric bending effects, this assumption is not sufficient anymore for lower elevation angles. Thus, in the following, a mixed 2D piecewise linear raytracing approach is introduced and possible error sources in the reconstruction of the bended signal paths are analyzed in more detail. Especially if low elevation observations are considered, unmodeled bending effects can introduce a systematic error of up to 10–20 ppm, on average 1–2 ppm, into the tomography solution. Thereby, not only the raytracing method but also the quality of the a priori field can have a significant impact on the reconstructed signal paths, if not reduced by iterative processing. In order to keep the processing time within acceptable limits, a bending model is applied for the upper part of the neutral atmosphere. It helps to reduce the number of processing steps by up to 85 % without significant degradation in accuracy. Therefore, the developed mixed raytracing approach allows not only for the correct treatment of low elevation observations but is also fast and applicable for nearrealtime applications.
Please read the corrigendum first before continuing.

Notice on corrigendum
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Article
(1370 KB)

The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(1370 KB)  Corrigendum
 BibTeX
 EndNote
For the conversion of precise integral measurements into 2 or 3D structures, a technique called tomography has been invented. In the field of GNSS meteorology, the principle of tomography became applicable with the increasing number of Global Navigation Satellite System (GNSS) satellites and the buildup of densified groundbased GNSS networks in the 1990s (Flores Jimenez, 1999; Raymond et al., 1994). Since then, a variety of tomography approaches based on raw GNSS phase measurements (Nilsson, 2005), double difference residuals (Kruse, 2001), slant delays (Flores Jimenez, 1999; Hirahara, 2000) or slant integrated water vapor (Champollion et al., 2005) have been developed for the accurate reconstruction of the water vapor distribution in the lower atmosphere. An overview about the major developments within this field of research since Flores Jimenez (1999) is provided by Manning (2013).
While in most tomography approaches, observations gathered at low elevation angles are discarded (Bender et al., 2011; Champollion et al., 2005; Hirahara, 2000), straightline signal path reconstruction is sufficient for the determination of the path lengths. However, Bender and Raabe (2007) showed that especially low elevation observations can be a very useful source of information in GNSS tomography. In addition to their information content about the lower troposphere, the additional observations strengthen the observation geometry and therewith contribute to a more reliable tomography solution. However, the correct treatment of low elevation observations requires more advanced raytracing algorithms. The first paper which deals with bended ray path reconstruction in GNSS tomography was published by Zus et al. (2015), with a main focus on the reconstruction of the signal paths for delay estimation but also for the assimilation of GNSS slant delays into numerical weather prediction systems. Most recently, Aghajany and Amerian (2017) published their results about 3D ray tracing in water vapor tomography and briefly analyzed its impact on the tomography solution.
Based on the existing studies, in the following, a more detailed discussion of possible error sources in signal path reconstruction is provided. Therefore, Sect. 2 describes the effect of atmospheric bending and its handling in GNSS signal processing. Section 3 describes the principles of GNSS tomography and how the basic equation of tomography is solved for wet refractivity. Section 4 introduces the concept of the reconstruction of signal paths using raytracing techniques. Here, the modified piecewise linear raytracing approach is described – including its ability for reconstruction of the GNSS signal geometry. In Sect. 5, the defined raytracing approach is applied to real slant wet delays (SWDs) and its impact on the tomography solution is assessed and validated against radiosonde data. Section 6 concludes the major findings.
The effect of atmospheric bending on GNSS signals is related to the propagation properties of electromagnetic waves. In a vacuum, GNSS signals travel at the velocity of light. When entering into the atmosphere, the electromagnetic wave velocity changes, dependent on the electric permittivity (ϵ) and magnetic permeability (μ) of the atmospheric constituents and the frequency of the electromagnetic wave. The ratio between the velocity of light c in a vacuum and the velocity ν in a medium defines the refractive index n.
For signals in the microwave frequency band, n ranges from 0.9996 to 1.0004. Thus, n is usually replaced by refractivity N, expressed in mm km^{−1} (ppm).
The GNSS signal delay in the lower atmosphere, also known as slant total delay (STD), is related to refractivity by the following equation (Bevis et al., 1992):
The first term of Eq. (3) describes the change in travel time due to velocity changes along the true ray path R. The second term (about 3 orders of magnitude smaller than the first term) is related to the difference in geometrical path length between the true (R) and the chord signal path (S). According to Dalton's law, the refractivity of air can be split up into a hydrostatic and a wet component: $N={N}_{\mathrm{h}}+{N}_{\mathrm{w}}$. Therefore, the GNSS signal delay reads
The slant wet delay (SWD) depends on the wet refractivity along the true ray path R.
The slant hydrostatic delay (SHD) results from the hydrostatic refractivity along R and, by definition, from the additional path length due to atmospheric bending.
While the signal path S follows from the straightline geometry between the satellite and the receiver, the true signal path R depends in addition on the hydrostatic and the wet refractivity distribution along the signal path (see Sect. 3 for more details).
In GNSS signal processing, the integral along the signal path is usually replaced by the zenith delay and a mapping function. Therefore, Eq. (4) is rewritten as follows:
where ZHD is the zenith hydrostatic delay, ZWD is the zenith wet delay and mf_{h} and mf_{w} are the corresponding mapping functions, which describe the elevation (ε) dependency of the signal delay. The elevationdependent and azimuth (α)dependent firstorder horizontally asymmetric term G(ε,α) reflects local variations in the atmospheric conditions; see MacMillan (1995), Chen and Herring (1997) or Landskron and Böhm (2018). In practice, e.g., when using the VMF1 mapping function (Böhm et al., 2006) or similar mapping concepts, the tropospheric delay due to atmospheric bending is absorbed by the hydrostatic mapping function term mf_{h}. Comparisons between raytraced SHD(ε) and “mapped” SHD$\left(\mathit{\epsilon}\right)=\text{ZHD}\cdot m{f}_{\mathrm{h}}\left(\mathit{\epsilon}\right)$ slant hydrostatic delays reveal that about 97 % of the atmospheric bending effect is compensated by the VMF1 hydrostatic mapping function (see Appendix A for further details).
According to Iyer and Hirahara (1993), the general principle of tomography is described as follows:
where f_{s} is the projection function, g_{s} is the object property function and ds is a small element of the ray path R along which the integration takes place. In GNSS tomography, g_{s} is usually replaced by wet refractivity N_{w}, and integral measure f_{s} by SWD (the prefactor of 10^{6} vanishes if ds is provided in kilometers and SWD in millimeters).
A full nonlinear solution of Eq. (9) for wet refractivity is not of practical relevance since according to Fermat's principle, firstorder changes of the ray path lead to secondorder changes in travel time. In consequence, by ignoring the path dependency in the inversion of N_{w} along ds and by assuming the ray path as a straight line, a linear tomography approach can be defined which is well applicable to SWDs above a 15^{∘} elevation angle (Möller, 2017). However, with decreasing elevation angle, the true signal path deviates significantly from a straight line. In consequence, by ignoring atmospheric bending, a systematic error is introduced in the tomography solution. In order to overcome this limitation, in the following, an iterative tomography approach is defined in which the bended signal path is approximated by small line segments. It is similar to the linear tomography approach in that the neutral atmosphere or parts of it are discretized in volume elements (voxels) in which the refractivity N_{w,k} in each voxel k is assumed as constant. Consequently, Eq. (9) can be replaced by
where d_{k} is the traveled distance in each voxel. Assuming l observations and m voxels, a linear equation system can be set up. In matrix notation it reads
where SWD is the observation vector of size (l,1), N_{w} is the vector of unknowns of size (m,1) and A is a matrix of size (l,m) which contains the partial derivatives of the slant wet delays with respect to the unknowns, i.e., the traveled distances d_{k} in each voxel.
Solving Eq. (11) for N_{w} requires the inversion of matrix A.
The inverse A^{−1} exists if A is squared and if the determinant of A is nonzero, otherwise matrix A is called singular. Unfortunately, singularity appears in GNSS tomography in most cases since the observation data are “incomplete” and matrix A is not of full rank. Therefore, Eq. (13) becomes illposed, i.e., not uniquely solvable. In order to find a solution which preserves most properties of an inverse, in the following, matrix A is replaced by the pseudo inverse A^{+}. According to Hansen (2000) the pseudo inverse is defined as follows:
where U and V are orthogonal, normalized left and right singular matrices of A and matrix S is a diagonal matrix, which contains the singular values in descending order. In case a priori information (N_{w0}) can be made available, it enters the tomography solution as first guess as follows:
where matrices U, V and S are obtained by singular value decomposition of matrix ${\mathbf{A}}^{T}\cdot \mathbf{P}\cdot \mathbf{A}+{\mathbf{P}}_{c}$. The weighting matrices P and P_{c} are defined as the inverse of the variance–covariance matrix C for the observations and C_{c} for the first guess, respectively. Assuming that the observations are uncorrelated, the nondiagonal elements of C and C_{c} are zero and the diagonal elements are defined as follows:
whereby σ_{ZWD}=2.5 mm reflects the uncertainty of the ZWD. The values for σ_{T}, σ_{q} and σ_{p} were taken from heightdependent error curves for pressure (p), temperature (T) and specific humidity (q) as provided by Steiner et al. (2006) for the ECMWF (European Centre for MediumRange Weather Forecasts) analysis data. For further details, the reader is referred to Möller (2017).
Assuming that the geometrical optics approximation is valid and that the atmospheric conditions change only inappreciably within one wavelength, the signal path is well reconstructible by means of raytracing shooting techniques (Hofmeister, 2016; Nievinski, 2009). Thereby, the basic equation for ray tracing, the socalled eikonal equation, has to be solved for obtaining optical path length L.
From Eq. (18), a number of 3D and 2D raytracing approaches have been derived for the reconstruction of groundbased and spacebased GNSS measurements and of their signal paths through the atmosphere (Hobiger et al., 2008; Zou et al., 1999).
The main difference between both observation types is related to the observation geometry. While for spacebased GNSS observations derived from limb sounding, the bending angle is usually described as a function of impact parameter a, for groundbased observations, elevation and azimuth angles are used for characterizing the signal geometry. In consequence, the optimal raytracing approach will be significantly different for various observation geometries.
In order to find an optimal approach for the operational analysis of groundbased measurements, Hofmeister (2016) carried out a number of exploratory comparisons. Based on the outcome, the 2D piecewise linear ray tracer was defined as the optimal reconstruction tool for the iterative reconstruction of the atmospheric signal delays including atmospheric bending. It is limited to positive elevation angles but it is fast and almost as accurate as the 3D ray tracer. However, for use in GNSS tomography, the raytracing approach had to be further modified. In the following, the developed raytracing approach but also its impact on the GNSS tomography solution are discussed in more detail.
4.1 Piecewise linear ray tracer
The starting point for the 2D piecewise linear ray tracer is the receiver position in ellipsoidal coordinates $({\mathit{\phi}}_{\mathrm{1}},{\mathit{\lambda}}_{\mathrm{1}},{h}_{\mathrm{1}})$, the “vacuum” elevation angle ε_{k} (see Fig. 1) and the azimuth angle α under which the satellite is observed. In the case of GNSS tomography, these parameters can be determined with sufficient accuracy from satellite ephemerides and the receiver position – assuming straightline geometry.
Therefore, the initial parameters for ray tracing (see Fig. 1), i.e., the geocentric coordinates (y_{1},z_{1}) and the corresponding geocentric angles (η_{1},θ_{1}), read
where R_{G} is the Gaussian radius, an adequate approximation of the Earth radius:
with a and b as the semiaxes of the reference ellipsoid (e.g., GRS80). The z axis connects the geocenter with the starting point; the y axis is defined perpendicular to the z axis in direction (azimuth angle) of the GNSS satellite in view. After setting the initial parameters, the “true” ray path is reconstructed iteratively by making use of raytracing shooting techniques. Therefore, total refractivity derived from an a priori field is read in and preprocessed for ray tracing. Here, the input data are interpolated vertically and horizontally to the vertical plane, spanned by the y and z axis.
In the raytracing loop, for each height layer h_{i+1} with $i=\mathrm{1}:(t\mathrm{1})$, whereby t defines the top layer of the voxel model, the geocentric coordinates and the corresponding angles are computed as follows:
where d_{i} is the reconstructed path length between the height layer h_{i} and h_{i+1} (${h}_{i+\mathrm{1}}>{h}_{i}$). It depends on the observation geometry but also on the atmospheric conditions (refractive indices n_{i} and n_{i+1}). By default, for our analysis, the spacing between two height layers h_{i} and h_{i+1} was set to 5 m, which corresponds to a maximum path length d_{i} of 100 m – assuming an elevation angle of 3^{∘} ($\mathrm{5}\mathrm{m}/\mathrm{sin}\mathrm{3}{}^{\circ}$).
The raytracing loop stops when the ray reaches the top layer t of the voxel model. Assuming spherical trigonometry, the spherical coordinates $({\mathit{\phi}}_{i+\mathrm{1}},{\mathit{\lambda}}_{i+\mathrm{1}})$ of the ray path segments are defined as follows:
where φ_{i+1} and λ_{i+1} are defined in the range $[\mathit{\pi}/\mathrm{2},\mathit{\pi}/\mathrm{2}]$ and $[\mathit{\pi},\mathit{\pi}]$, respectively. The ray coordinates are necessary for interpolation of the refractive indices n_{i} and n_{i+1} for the next processing step i but also for computation of the intersection points with the voxel model boundaries.
The raytracing loop is repeated until ${\mathit{\epsilon}}_{t}{\mathit{\epsilon}}_{k}+{g}_{\mathrm{bend}}$ is smaller than a predefined threshold (e.g., 10^{−6}^{∘}). While the elevation angle ε_{t} is obtained by Eq. (29) for $i=t\mathrm{1}$, the correction term g_{bend} accounts for the additional bending above the voxel model. Since the atmosphere is almost in a state of hydrostatic equilibrium, g_{bend} can be well approximated by a bending model, like the one of Hobiger et al. (2008):
where h is replaced by h_{t}, the height of the voxel top layer. After convergence of the raytracing loop, the path length in each voxel is obtained by summing up the distances d_{i} in each voxel. Thereby, allocation of the ray parts is carried out by comparison of the ray coordinates $({\mathit{\phi}}_{i},{\mathit{\lambda}}_{i},{h}_{i})$ with the coordinates of the voxel model. The obtained ray paths in each voxel – for each station and each satellite in view  are used for setting up design matrix A (see Eq. 12).
4.2 Quality of reconstructed ray paths
4.2.1 The refractivity field
The quality of the raytraced signal paths depends primarily on the quality of the refractivity field. Especially if no good a priori data can be made available, e.g., if standard atmosphere (StdAtm) is used instead of numerical weather model data (ALARO), the reconstructed signal path might deviate significantly from the true signal path.
Figure 2 shows the impact of the refractivity field on the signal geometry as an example for a GNSS signal observed at Jenbach station, Austria $(\mathit{\phi}=\mathrm{47.4}{}^{\circ}$, $\mathit{\lambda}=\mathrm{11.8}{}^{\circ}$, h=545 m), with ε=5^{∘} and $\mathit{\alpha}=\mathrm{230}{}^{\circ}$. At this particular epoch (4 May 2013 at 15:00 UTC), standard atmosphere deviates by about 30 ppm from the ALARO model data. Assuming ALARO as the reference, ray tracing through the standard atmosphere causes a ray deviation of 100–200 m (see Fig. 2b).
In order to reduce the impact of possible refractivity errors on the reconstructed ray paths and in further consequence on the tomography solution, ray tracing was carried out iteratively. Therefore, the refractivity field obtained from the first tomography solution replaces the initial refractivity field for ray tracing for the next iteration and so on. The processing is repeated until N_{w} converges.
Figure 3a shows the convergence behavior assuming standard atmosphere (StdAtm) and ALARO model data as input. In both cases, the standard deviations of the differences in path length between two consecutive epochs (${d}_{k,i+\mathrm{1}}{d}_{k,i}$) were selected as convergence criteria. Both solutions converge after two iterations, thereby “improving” the path lengths within each voxel by about 22 m in the case of the standard atmosphere and by 11 m in the case of ALARO data. This result was expected, since ALARO data are closer to the true atmospheric conditions. By comparison of Fig. 3b with Fig. 2b, it is clearly visible that the two additional iterations help to reduce the ray offset caused by errors in the standard atmosphere from 100–200 to 30–40 m. In Sect. 5 the resulting effect on the tomography solution is assessed.
4.2.2 The empirical raybending model
In addition to the refractivity field, the quality of the reconstructed ray paths might also be affected by errors in the bending model as defined by Eq. (32). Comparisons of the bending model with raytraced bending angles on a global $\mathrm{10}{}^{\circ}\times \mathrm{10}{}^{\circ}$ grid over the period of 1 year reveal that the error in bending is usually kept below 0.8 arcsec. Assuming a GNSS site near sea level and an elevation angle of 5^{∘}, an error in bending angle of ±0.8 arcsec causes an error in path length of up to ±10 m; i.e., the reconstructed GNSS signal enters the voxel model slightly earlier or later than the observed GNSS signal. In Fig. 4, the bending error is visualized as a pointing error at the voxel model top. However, for the tomography solution, this effect is too small to be significant. Thus, it can be concluded that the bending model of Hobiger et al. (2008) is well applicable for the reconstruction of the bending angle above the voxel model, in particular if the voxel model height h_{t} is set to 12 km or higher.
4.2.3 Ionospheric bending effects
Beyond, the ionosphere also influences GNSS signal propagation. In order to assess the impact of free electrons in the ionosphere (above 80 km altitude) on the signal path, the electron density model by Anderson et al. (1987) was executed in three scenarios, assuming a vertical total electron content $(\text{VTEC}=\int {N}_{\mathrm{e}}\cdot \mathrm{d}h)$ of 34 TECU (average daytime), 120 TECU (solar maximum) and 455 TECU (maximum possible; see Wijaya, 2010), respectively.
By making use of Eq. (33), the obtained electron density profiles were converted into profiles of refractivity (N), assuming signal frequency f_{1}=1575.42 MHz (GPS L1) and f_{2}=1227.60 MHz (GPS L2). Figure 5 shows the obtained vertical profiles of ionospheric refractivity as an example for frequency f_{1}. The higher the signal frequency f, the lower the phase velocity through the ionosphere and the less its refraction.
Following the approach by Wijaya (2010), the ray paths in the ionosphere were reconstructed separately for GPS L1 and L2. The analysis revealed significant path differences between the true ray path and its chord line but also between the two signal frequencies. Assuming a VTEC of 455 TECU and an elevation angle of 3^{∘}, the maximum deviation from the straightline signal path is 800 m for L1 and 550 m for L2, respectively, at h=400 km, slightly below the layer of peak electron density. Fortunately, ray path deviation decreases significantly with decreasing VTEC and altitude to a few tens of meters at h=13.6 km (the upper rim of the troposphere at which the top of the voxel model was defined). In consequence, the impact of free electrons on the signal path in the lower atmosphere is negligible under moderate and low ionospheric conditions.
In the following, the differences between straightline and bended ray tracing are further analyzed. For a high degree of consistency, the ray tracing approach defined in Sect. 4 was used for both straightline and bended ray tracing. The only difference is that in the case of straightline ray tracing the ratio ${n}_{i}/{n}_{i+\mathrm{1}}$ in Eq. (27) was set to 1, thereby guaranteeing that only the impact of atmospheric refraction is assessed.
5.1 Expected drying effect
In the beginning, the ray position is equal for both methods but diverges with increasing height. Thereby, in most cases, the bended ray is traveling “above” the straight ray; i.e., the straight ray enters the voxel model top “earlier” than the bended ray. This leads to the effect that the straight ray remains in the voxel model longer than the bended ray; i.e., the straight ray path within the voxel model (h<13.6 km) is systematically longer than the bended ray path. The differences between both ray paths are plotted in Fig. 6a as a function of elevation angle. Therefore, ALARO model data were selected as input for the bended ray tracer.
The additional ray path decreases rapidly with increasing elevation angle. Thus, a mixed raytracing approach can be defined, which considers ray bending only for $\mathit{\epsilon}\le \mathrm{15}{}^{\circ}$. At higher elevation angles, the additional ray path is below 0.1 km, and straightline ray tracing is sufficient for ray path reconstruction.
Figure 6a also shows that in some cases, even for low elevation angles, the difference in path length is small (below 0.1 km). This appears when the ray enters the voxel model not through the top layer but through the lateral surface of the voxel model. In this particular case, the difference in path length between both raytracing approaches is negligible (be aware that only the entire distance through all voxels is comparable for both raytracing approaches, not the individual distances in each voxel). Figure 6b and c show the expected drying effects in the tomography solution caused by errors in the reconstructed signal paths assuming straightline geometry. Here, it is distinguished between the drying effect caused by the additional ray path (dNw_{1}; panel b) and the drying effect caused by the fact that the straight line travels through lower layers of the voxel model (dNw_{2}; panel c). Both effects were assessed as follows:
where SWD_{b} and SWD_{s} are the slant wet delays obtained by ray tracing through ALARO model data along the bended and the straightline ray path, respectively. The variables d_{k,b} and d_{k,s} are the corresponding path lengths within the voxel model. The sum of their differences along the ray paths is identical to the additional ray paths plotted in Fig. 6a.
Both drying effects have to be considered as additive and are strongly connected to the current atmospheric conditions as well as to the parametrization applied for interpolation of the refractivity field. In our analysis we assumed an exponential decrease of refractivity between the vertical layers of the voxel model and applied a bilinear interpolation method for horizontal interpolation between the grid points.
5.2 Results from the Austrian GNSS tomography test case
In order to study the impact of bended ray tracing on the tomography solution, a GNSS tomography test case was defined. The corresponding settings are summarized in Table 1.
Figure 7a shows the differences in wet refractivity between Sol1 and Sol2 (as defined in Table 1). Even though on average over all voxels no bias in wet refractivity is observed, specific voxels show differences in wet refractivity of up to 10 ppm, particularly if different voxels than in the straightline solution are traversed due to bending.
Figure 7b shows the differences in wet refractivity between the first two iterations of the mixed raytracing approach (Sol2). In this particular case, refractivity differences are smaller than 0.05 ppm, which implies that the a priori model used for ray tracing is already close to the true atmospheric conditions; i.e., in this particular case no further iteration was necessary.
From all differences in wet refractivity over 248 epochs in May 2013, a maximum of 14.2 ppm, a bias of 0.12 ppm and a standard deviation of 0.24 ppm were obtained. Although the bias and the standard deviation over all voxels are small, differences of about 1 ppm were observed on average at each epoch, especially when observations below a 10^{∘} elevation angle enter the tomography solution.
5.3 Validation with radiosonde data
For validation of the mixed raytracing approach against straightline ray tracing, the tomographyderived wet refractivity fields were compared with radiosonde data at the airport of Innsbruck (${\mathit{\phi}}_{i}=\mathrm{47.3}{}^{\circ}$, $\mathit{\lambda}=\mathrm{11.4}{}^{\circ}$, h=579 m). First, the radiosonde data obtained once a day between 02:00 and 03:00 UTC were preprocessed; i.e., outliers in temperature were removed and dew point temperature was converted to water vapor pressure and further to wet refractivity. Finally, the radiosonde profiles were vertically interpolated to the height layers of the voxel model and the tomographyderived wet refractivity fields were horizontally interpolated to the ground position of the radiosonde launching site, respectively. As an example, Fig. 8 shows the differences in wet refractivity as a function of height above surface for two epochs in May 2013. In both cases, the bended raytracing approach helps to reduce the tomography error by about 1–2 ppm, especially in the lower 4 km of the atmosphere. Largest differences are visible when the bended ray traverses other voxels than its chord line. This appeared in about 2 % of the test cases, especially if observations below a 10^{∘} elevation angle enter the tomography solution.
GNSS signals which enter the neutral atmosphere at low elevation angles (ε<15^{∘}) are significantly affected by atmospheric bending. In the case that the bending is neglected when setting up design matrix A, a systematic error of up to 10–20 ppm, on average 1–2 ppm, is introduced into the GNSS tomography solution. This error can be widely reduced if atmospheric bending is considered in the reconstruction of the signal paths. Therefore, a 2D piecewise linear raytracing approach was defined, which describes the bended GNSS signal path by small line segments. By limiting the length of the line segments to 100 m in the case of $\mathit{\epsilon}=\mathrm{3}{}^{\circ}$ or even shorter for higher elevation angles, the true signal path can be widely reconstructed. However, the quality of the reconstructed signal paths depends primarily on the quality of the a priori refractivity field. Comparisons between refractivity fields derived from standard atmosphere and ALARO weather model data reveal that a refractivity error of 30 ppm can cause a ray deviation of up to several hundred meters; i.e., the distance traveled in each voxel but also the number of traversed voxels are prone to misallocations. In consequence, reliable a priori data, e.g., derived from numerical weather model data, are recommended for GNSS tomography.
Nevertheless, if reliable a priori data are not available or if the quality is unknown, iterative ray tracing helps to reduce the impact of wet refractivity errors on the tomography solution. Therefore, the wet refractivity field obtained from an initial tomography solution is used for reconstruction of the signal paths for the next iteration. The processing is repeated until the tomography solution converges. This ensues usually after two iterations. Further, a bending model, like the one provided by Hobiger et al. (2008), helps to significantly reduce computational cost by describing the remaining bending in the higher atmosphere (above the voxel model). In consequence, the ray tracer can be stopped right after the reconstructed signal leaves the voxel model. In the case of h_{t}=13.6 km, the number of processing steps is reduced by 85 %, which is a tremendous reduction in processing time without a significant loss of accuracy.
In contrast, ionospheric bending effects have less impact on the GNSS tomography solution. Even during periods of solar maximum, ray path deviation caused by ionospheric bending is negligible for signals in the L band (1–2 GHz). However, even if ionospheric bending has no impact on the tomography solution, first and higher order ionospheric effects should be taken into account when processing GNSS phase observations.
In addition, comparisons with radiosonde data revealed that if atmospheric bending effects are considered in GNSS tomography, the quality of the tomography solution can be improved by 1–2 ppm. Within the defined test case, especially voxels in the lower 4 km of the atmosphere benefitted from the applied mixed raytracing approach. Due to significant optimization, the mixed raytracing approach ensures processing of large tomography test cases in adequate time. A test case with 72 GNSS sites and $\mathrm{7}\times \mathrm{9}\times \mathrm{15}$ voxels can be processed in less than 2 min. Thus, the developed mixed raytracing approach is also applicable in nearreal time and therefore well suited for operational purposes.
The 2D piecewise linear ray tracer for GNSS tomography as well as the RADIATE ray tracer are part of the Vienna VLBI and Satellite Software (VieVS). The code of the RADIATE ray tracer is available at https://github.com/TUWVieVS/RADIATE (last access: 19 December 2018). For more details on VieVS, the reader is referred to http://vievswiki.geo.tuwien.ac.at (last access: 19 December 2018).
In the case of VMF1 (Böhm et al., 2006) or similar mapping concepts, azimuthal asymmetry is not considered, and for convenience, only a single hydrostatic mapping coefficient per site (a_{h}) is determined as follows:
where b_{h} is 0.0029, c_{h} depends on the day of year and latitude and mf_{h}(ε) is defined as the ratio between SHD(3^{∘}) and ZHD, obtained by ray tracing through numerical weather model data. For assessing the remaining unmodeled geometric bending dg_{bend}(ε,α), raytraced slant hydrostatic delays were compared with mapped slant hydrostatic delays as follows:
where ZHD is the zenith hydrostatic delay obtained by vertical integration, g_{bend}(ε,α) is the geometric bending effect as obtained by ray tracing, mf_{h}(ε) is the VMF1 hydrostatic mapping function determined by $\text{SHD}\left(\mathrm{3}{}^{\circ}\right)/\text{ZHD}$ and mf_{h0}(ε) is the hydrostatic mapping function determined by ${\text{SHD}}_{\mathrm{0}}\left(\mathrm{3}{}^{\circ}\right)/\text{ZHD}$, where SHD(3^{∘}) and SHD_{0}(3^{∘}) are the slant hydrostatic delays obtained by ray tracing for a vacuum elevation angle ${\mathit{\epsilon}}_{k}=\mathrm{3}{}^{\circ}$ with and without geometric bending, respectively. Figure A1 shows the remaining unmodeled geometric bending as obtained for six elevation angles (and 16 equidistant azimuth angles) as an example for the two VLBI sites in Fortaleza, Brazil ($\mathit{\phi}=\mathrm{3.9}{}^{\circ}$, $\mathit{\lambda}=\mathrm{321.6}{}^{\circ}$; h=23 m), and Wettzell, Germany ($\mathit{\phi}=\mathrm{49.1}{}^{\circ}$, $\mathit{\lambda}=\mathrm{12.9}{}^{\circ}$; h=669 m).
In the case of $\mathit{\epsilon}=\mathrm{3}{}^{\circ}$, almost no bending error is visible since mf_{h}(ε) was tuned for this elevation angle. However, for other elevation angles, the unmodeled geometric bending is about 3 % of the slant hydrostatic delay, e.g., up to ±5 mm at a 5^{∘} elevation angle. In the case of Wettzell, dg_{bend}(ε,α) is mostly negative; i.e., the mapped SHD is smaller than the observed SHD, and vice versa for Fortaleza. So far, these small variations have been neglected when using the VMF1 hydrostatic mapping function in GNSS signal processing.
GM, as main author, did most of the analysis, drafted the manuscript and designed the figures. DL provided the raytraced delays based on the RADIATE ray tracer, contributed to the definition of the raytracing approach for GNSS tomography and gave feedback on the draft.
The authors declare that they have no conflict of interest.
This article is part of the special issue “Advanced Global Navigation Satellite Systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC) (AMT/ACP/ANGEO interjournal SI)”. It is not associated with a conference.
Open access funding was provided by the Austrian Science Fund (FWF). The authors
would like to thank the Austrian Science Fund (FWF) for financial support of
this study within the project RADIATE ORD (ORD 86) and the Austrian Research
Promotion Agency (FFG) for financial support within the project GNSSATom
(840098).
Edited by: Jonathan Jones
Reviewed by: two anonymous referees
Aghajany, S. H. and Amerian, Y.: Three dimensional ray tracing technique for tropospheric water vapor tomography using GPS measurements, J. Atmos. Sol.Terr. Phy., 164, 81–88, https://doi.org/10.1016/j.jastp.2017.08.003, 2017. a
Anderson, D. N., Mendillo, M., and Herniter, B.: A semiempirical lowlatitude ionospheric model, Radio Sci., 22, 292–306, https://doi.org/10.1029/RS022i002p00292, 1987. a
Bender, M. and Raabe, A.: Preconditions to ground based GPS water vapour tomography, Ann. Geophys., 25, 1727–1734, https://doi.org/10.5194/angeo2517272007, 2007. a
Bender, M., Stosius, R., Zus, F., Dick, G., Wickert, J., and Raabe, A.: GNSS water vapour tomography  Expected improvements by combining GPS, GLONASS and Galileo observations, Adv. Space Res., 47, 886–897, https://doi.org/10.1016/j.asr.2010.09.011, 2011. a
Bevis, M., Businger, S., Herring, T. A., Rocken, C., Anthes, R. A., and Ware, R. H.: GPS meteorology: Remote sensing of atmospheric water vapor using the Global Positioning System, J. Geophys. Res., 97, 15787–15801, https://doi.org/10.1029/92JD01517, 1992. a
Böhm, J., Werl, B., and Schuh, H.: Troposphere mapping functions for GPS and Very Long Baseline Interferometry from European Centre for MediumRange Weather Forecasts operational analysis data, J. Geophys. Res., 111, 1–9, https://doi.org/10.1029/2005JB003629, 2006. a, b
Champollion, C., Masson, F., Bouin, M.N., Walpersdorf, A., Doerflinger, E., Bock, O., and Bealen, J. V.: GPS water vapour tomography: preliminary results from the ESCOMPTE field experiment, Atmos. Res., 74, 253–274, https://doi.org/10.1016/j.atmosres.2004.04.003, 2005. a, b
Chen, G. and Herring, T. A.: Effects of atmospheric azimuthal asymmetry on the analysis of space geodetic data, J. Geophys. Res., 102, 20489–20502, https://doi.org/10.1029/97JB01739, 1997. a
Flores Jimenez, A.: Atmospheric tomography using satellite radio signals, PhD thesis, Universitat Politecnica de Catalunya, Departament de Teoria del Senyal i Comunicacions, available at: http://www.tdx.cat/handle/10803/6881 (last access: 19 December 2018), 1999. a, b, c
Hansen, P. C.: The Lcurve and its use in the numerical treatment of inverse problems, in: Computational Inverse Problems in Electrocardiology, Advances in Computational Bioengineering, ISBN: 1853126144, 119–142, WIT Press, Southampton, UK, 2000. a
Hirahara, K.: Local GPS tropospheric tomography, Earth Planets Space, 52, 935–939, https://doi.org/10.1186/BF03352308, 2000. a, b
Hobiger, T., Ichikawa, R., Koyama, Y., and Kondo, T.: Fast and accurate raytracing algorithms for realtime space geodetic applications using numerical weather models, J. Geophys. Res., 113, 1–14, https://doi.org/10.1029/2008JD010503, 2008. a, b, c, d, e
Hofmeister, A.: Determination of path delays in the atmosphere for geodetic VLBI by means of raytracing, PhD thesis, TU Wien, Department of Geodesy and Geoinformation, available at: http://repositum.tuwien.ac.at/obvutwhs/content/titleinfo/1370367 (last access: 19 December 2018), 2016. a, b
Iyer, H. M. and Hirahara, K.: Seismic tomography: Theory and practice, Springer Netherlands, iSBN: 0412371901, 1993. a
Kruse, L. P.: Spatial and Temporal Distribution of Atmospheric Water Vapor using Space Geodetic Techniques, Geodätischgeophysikalische Arbeiten in der Schweiz, 61, 2001. a
Landskron, D. and Böhm, J.: Refined discrete and empirical horizontal gradients in VLBI analysis, J. Geod., 92, 1387–1399, https://doi.org/10.1007/s0019001811271, 2018. a
MacMillan, D. S.: Atmospheric gradients from Very Long Baseline Interferometry observations, Geophys. Res. Lett., 22, 1041–1044, https://doi.org/10.1029/95GL00887, 1995. a
Manning, T.: Sensing the dynamics of severe weather using 4D GPS tomography in the Australien region, PhD thesis, Royal Melbourne Institute of Technology, School of Mathematical and Geospatial Sciences, available at: http://researchbank.rmit.edu.au/eserv/rmit:160752/Manning.pdf (last access: 19 December 2018), 2013. a
Möller, G.: Reconstruction of 3D wet refractivity fields in the lower atmosphere along bended GNSS signal paths, PhD thesis, TU Wien, Department of Geodesy and Geoinformation, available at: http://repositum.tuwien.ac.at/obvutwhs/content/titleinfo/2268559 (last access: 19 December 2018), 2017. a, b
Nievinski, F. G.: Raytracing to mitigate the neutral atmosphere delay in GPS, Master's thesis, University of New Brunswick, Department of Geodesy and Geomatics Engineering, available at: http://www2.unb.ca/gge/Pubs/TR262.pdf (last access: 19 December 2018), 2009. a
Nilsson, T.: Assessment of tomographic methods for estimation of atmospheric water vapor using groundbased GPS, Ph.D. thesis, Chalmers University of Technology, Göteborg, Department of Radio and Space Science, Sweden, 2005. a
Raymond, T. D., Franke, S. J., and Yeh, K. C.: Ionospheric tomography: its limitations and reconstruction methods, J. Atmos. Terr. Phys., 56, 637–657, https://doi.org/10.1016/00219169(94)90104X, 1994. a
Steiner, A. K., Löscher, A., and Kirchengast, G.: Error characteristics of Refractivity Profiles Retrieved from CHAMP Radio Occultation Data, Springer Berlin Heidelberg, Atmosphere and Climate: Studies by Occultation Methods, 27–36, https://doi.org/10.1007/3540341218_3, 2006. a
Wijaya, D. D.: Atmospheric correction formulae for space geodetic techniques, PhD thesis, Graz University of Technology, Faculty of Technical Mathematics and Technical Physics, available at: http://www.shaker.de/OnlineGesamtkatalogDownload/2018.12.1920.30.09207.151.220.142rad8778A.tmp/3832289941_INH.PDF (last access: 19 December 2018), 2010. a, b
Zou, X., Vandenberghe, F., Wang, B., Gorbunov, M. E., Kuo, Y., Sokolovskiy, S., Chang, J. C., Sela, J. G., and Anthes, R. A.: A raytracing operator and its adjoint for the use of GPS/MET refraction angle measurements, J. Geophys. Res.Atmos., 104, 22301–22318, https://doi.org/10.1029/1999JD900450, 1999. a
Zus, F., Dick, G., Heise, S., and Wickert, J.: A forward operator and its adjoint for GPS slant total delays, Radio Sci., 50, 393–405, https://doi.org/10.1002/2014RS005584, 2015. a
 Abstract
 Introduction
 Atmospheric bending effects in GNSS signal processing
 The principles of GNSS tomography
 Reconstruction of GNSS signal paths
 Impact of atmospheric bending on the tomography solution
 Conclusions
 Code availability
 Appendix A: Unmodeled bending effects in the Vienna hydrostatic mapping function
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 References
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(1370 KB)  Fulltext XML
 Abstract
 Introduction
 Atmospheric bending effects in GNSS signal processing
 The principles of GNSS tomography
 Reconstruction of GNSS signal paths
 Impact of atmospheric bending on the tomography solution
 Conclusions
 Code availability
 Appendix A: Unmodeled bending effects in the Vienna hydrostatic mapping function
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 References