Atmospheric bending effects in GNSS tomography

In Global Navigation Satellite System (GNSS) tomography, precise information about the tropospheric water vapor distribution is derived from integral measurements like ground-based 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 2-D piecewise linear ray-tracing 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 ray-tracing 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 ray-tracing approach allows not only for the correct treatment of low elevation observations but is also fast and applicable for near-real-time applications.

1 Introduction

For the conversion of precise integral measurements into 2- or 3-D 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 build-up of densified ground-based 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), straight-line 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 ray-tracing 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 3-D 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 ray-tracing techniques. Here, the modified piecewise linear ray-tracing approach is described – including its ability for reconstruction of the GNSS signal geometry. In Sect. 5, the defined ray-tracing 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.

2 Atmospheric bending effects in GNSS signal processing

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.

$$\begin{array}{}\text{(1)}& n={\displaystyle \frac{c}{\mathit{\nu }}}=\sqrt{{\displaystyle \frac{\mathit{ϵ}\cdot \mathit{\mu }}{{\mathit{ϵ}}_{\mathrm{0}}\cdot {\mathit{\mu }}_{\mathrm{0}}}}}\end{array}$$

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⁻¹ (ppm).

$$\begin{array}{}\text{(2)}& N={\mathrm{10}}^{\mathrm{6}}\cdot (n-\mathrm{1})\end{array}$$

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):

$$\begin{array}{}\text{(3)}& \text{STD}={\mathrm{10}}^{-\mathrm{6}}\cdot \underset{R}{\int }N\cdot \mathrm{d}s+\left[\underset{R}{\int }\mathrm{d}s-\underset{S}{\int }\mathrm{d}s\right].\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }\text{STD}& {\displaystyle }=\text{SHD}+\text{SWD}={\mathrm{10}}^{-\mathrm{6}}\cdot \underset{R}{\int }{N}_{\mathrm{h}}\cdot \mathrm{d}s+{\mathrm{10}}^{-\mathrm{6}}\cdot \underset{R}{\int }{N}_{\mathrm{w}}\cdot \mathrm{d}s\\ \text{(4)}& {\displaystyle }& {\displaystyle }+\left[\underset{R}{\int }\mathrm{d}s-\underset{S}{\int }\mathrm{d}s\right].\end{array}$$

The slant wet delay (SWD) depends on the wet refractivity along the true ray path R.

$$\begin{array}{}\text{(5)}& \text{SWD}={\mathrm{10}}^{-\mathrm{6}}\cdot \underset{R}{\int }{N}_{\mathrm{w}}\cdot \mathrm{d}s\end{array}$$

The slant hydrostatic delay (SHD) results from the hydrostatic refractivity along R and, by definition, from the additional path length due to atmospheric bending.

$$\begin{array}{}\text{(6)}& \text{SHD}={\mathrm{10}}^{-\mathrm{6}}\cdot \underset{R}{\int }{N}_{\mathrm{h}}\cdot \mathrm{d}s+\left[\underset{R}{\int }\mathrm{d}s-\underset{S}{\int }\mathrm{d}s\right]\end{array}$$

While the signal path S follows from the straight-line 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:

$$\begin{array}{ll}{\displaystyle }\text{STD}(\mathit{\epsilon },\mathit{\alpha })& {\displaystyle }=\text{SHD}+\text{SWD}\\ \text{(7)}& {\displaystyle }& {\displaystyle }=\text{ZHD}\cdot m{f}_{\mathrm{h}}\left(\mathit{\epsilon }\right)+\text{ZWD}\cdot m{f}_{\mathrm{w}}\left(\mathit{\epsilon }\right)+G(\mathit{\epsilon },\mathit{\alpha }),\end{array}$$

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 elevation-dependent and azimuth (α)-dependent first-order 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 ray-traced 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).

3 The principles of GNSS tomography

According to Iyer and Hirahara (1993), the general principle of tomography is described as follows:

$$\begin{array}{}\text{(8)}& f_s=\underset{R}{\int }{g}_s\cdot \mathrm{d}s,\end{array}$$

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⁶ vanishes if ds is provided in kilometers and SWD in millimeters).

$$\begin{array}{}\text{(9)}& \text{SWD}=\underset{R}{\int }{N}_{\mathrm{w}}\cdot \mathrm{d}s\end{array}$$

A full nonlinear solution of Eq. (9) for wet refractivity is not of practical relevance since according to Fermat's principle, first-order changes of the ray path lead to second-order 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

$$\begin{array}{}\text{(10)}& \text{SWD}=\sum _{k=\mathrm{1}}^m{N}_{\mathrm{w},k}\cdot d_k,\end{array}$$

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

$$\begin{array}{}\text{(11)}& \mathbf{SWD}=\mathbf{A}\cdot {\mathit{N}}_{\mathrm{w}},\end{array}$$

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.

$$\begin{array}{}\text{(12)}& \mathbf{A}=\left[\begin{array}{ccc}\frac{\mathit{\delta }{\text{SWD}}_{\mathrm{1}}}{\mathit{\delta }{N}_{\mathrm{w},\mathrm{1}}}& \mathrm{\cdots }& \frac{\mathit{\delta }{\text{SWD}}_{\mathrm{1}}}{\mathit{\delta }{N}_{\mathrm{w},m}}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \frac{\mathit{\delta }{\text{SWD}}_l}{\mathit{\delta }{N}_{\mathrm{w},\mathrm{1}}}& \mathrm{\cdots }& \frac{\mathit{\delta }{\text{SWD}}_l}{\mathit{\delta }{N}_{\mathrm{w},m}}\end{array}\right]\end{array}$$

Solving Eq. (11) for N_w requires the inversion of matrix A.

$$\begin{array}{}\text{(13)}& {\mathit{N}}_{\mathrm{w}}={\mathbf{A}}^{-\mathrm{1}}\cdot \mathbf{SWD}\end{array}$$

The inverse A⁻¹ 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 ill-posed, 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:

$$\begin{array}{}\text{(14)}& {\mathbf{A}}^{+}=\mathbf{V}\cdot {\mathbf{S}}^{-\mathrm{1}}\cdot {\mathbf{U}}^T,\end{array}$$

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:

$$\begin{array}{}\text{(15)}& {\mathit{N}}_{\mathrm{w}}={\mathit{N}}_{\mathrm{w}\mathrm{0}}+\mathbf{V}\cdot {\mathbf{S}}^{-\mathrm{1}}\cdot {\mathbf{U}}^T\cdot {\mathbf{A}}^T\cdot \mathbf{P}\cdot (\mathbf{SWD}-\mathbf{A}\cdot {\mathit{N}}_{\mathrm{w}\mathrm{0}}),\end{array}$$

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 non-diagonal elements of C and C_c are zero and the diagonal elements are defined as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{C}}^{\mathrm{2}}={\sin }^{\mathrm{2}}\mathit{\epsilon }\cdot {\mathit{\sigma }}_{\mathrm{ZWD}}^{\mathrm{2}}\text{(17)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{{\mathrm{C}}_c}^{\mathrm{2}}={\left({\displaystyle \frac{\partial N_{\mathrm{w}}}{\partial T}}\cdot {\mathit{\sigma }}_T\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial N_{\mathrm{w}}}{\partial q}}\cdot {\mathit{\sigma }}_q\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial N_{\mathrm{w}}}{\partial p}}\cdot {\mathit{\sigma }}_p\right)}^{\mathrm{2}},\end{array}$$

whereby σ_ZWD=2.5 mm reflects the uncertainty of the ZWD. The values for σ_T, σ_q and σ_p were taken from height-dependent error curves for pressure (p), temperature (T) and specific humidity (q) as provided by Steiner et al. (2006) for the ECMWF (European Centre for Medium-Range Weather Forecasts) analysis data. For further details, the reader is referred to Möller (2017).

4 Reconstruction of GNSS signal paths

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 ray-tracing shooting techniques (Hofmeister, 2016; Nievinski, 2009). Thereby, the basic equation for ray tracing, the so-called eikonal equation, has to be solved for obtaining optical path length L.

$$\begin{array}{}\text{(18)}& \left|\right|\mathit{▿}L|{|}^{\mathrm{2}}=n(\mathit{r}{)}^{\mathrm{2}}\end{array}$$

From Eq. (18), a number of 3-D and 2-D ray-tracing approaches have been derived for the reconstruction of ground-based and space-based 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 space-based GNSS observations derived from limb sounding, the bending angle is usually described as a function of impact parameter a, for ground-based observations, elevation and azimuth angles are used for characterizing the signal geometry. In consequence, the optimal ray-tracing approach will be significantly different for various observation geometries.

In order to find an optimal approach for the operational analysis of ground-based measurements, Hofmeister (2016) carried out a number of exploratory comparisons. Based on the outcome, the 2-D 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 3-D ray tracer. However, for use in GNSS tomography, the ray-tracing approach had to be further modified. In the following, the developed ray-tracing 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 2-D 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 straight-line geometry.

Figure 1 Geometry of the ray-tracing approach with the geocentric coordinates (y,z), the geocentric angles (η,θ), elevation angle ε and d as the distance between two consecutive ray points.

Therefore, the initial parameters for ray tracing (see Fig. 1), i.e., the geocentric coordinates (y₁,z₁) and the corresponding geocentric angles (η₁,θ₁), read

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{y}_{\mathrm{1}}=\mathrm{0}\text{(20)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{1}}=R_{\mathrm{G}}+h_{\mathrm{1}}\text{(21)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathrm{1}}=\mathrm{0}\text{(22)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{1}}={\mathit{\epsilon }}_k,\end{array}$$

where R_G is the Gaussian radius, an adequate approximation of the Earth radius:

$$\begin{array}{}\text{(23)}& R_{\mathrm{G}}={\displaystyle \frac{a^{\mathrm{2}}\cdot b}{(a\cdot \cos {\mathit{\phi }}_{\mathrm{1}}{)}^{\mathrm{2}}+(b\cdot \sin {\mathit{\phi }}_{\mathrm{1}}{)}^{\mathrm{2}}}},\end{array}$$

with a and b as the semi-axes 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 ray-tracing 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 ray-tracing 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:

$$\begin{array}{ll}\text{(24)}& {\displaystyle }& {\displaystyle }{y}_{i+\mathrm{1}}=y_i+d_i\cdot \cos {\mathit{\epsilon }}_i\text{(25)}& {\displaystyle }& {\displaystyle }{z}_{i+\mathrm{1}}=z_i+d_i\cdot \sin {\mathit{\epsilon }}_i\text{(26)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{i+\mathrm{1}}=\arctan {\displaystyle \frac{y_{i+\mathrm{1}}}{z_{i+\mathrm{1}}}}\text{(27)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{i+\mathrm{1}}=\arccos \left({\displaystyle \frac{n_i}{n_{i+\mathrm{1}}}}\cdot \cos ({\mathit{\theta }}_i+{\mathit{\eta }}_{i+\mathrm{1}}-{\mathit{\eta }}_i)\right){\displaystyle }& {\displaystyle }{d}_i=-(R_G+h_i)\cdot \sin {\mathit{\theta }}_i\\ \text{(28)}& {\displaystyle }& {\displaystyle }+\sqrt{(R_G+h_{i+\mathrm{1}}{)}^{\mathrm{2}}-(R_G+h_i{)}^{\mathrm{2}}\cdot {\cos }^{\mathrm{2}}{\mathit{\theta }}_i}\text{(29)}& {\displaystyle }& {\displaystyle }{\mathit{\epsilon }}_{i+\mathrm{1}}={\mathit{\theta }}_{i+\mathrm{1}}-{\mathit{\eta }}_{i+\mathrm{1}},\end{array}$$

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}/\sin \mathrm{3}{}^{\circ }$).

The ray-tracing 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:

$$\begin{array}{ll}{\displaystyle }{\mathit{\phi }}_{i+\mathrm{1}}& {\displaystyle }=\arcsin \left(\sin {\mathit{\phi }}_{\mathrm{1}}\cdot \cos ({\mathit{\eta }}_{i+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{1}})\right.\\ \text{(30)}& {\displaystyle }& {\displaystyle }\left.+\cos {\mathit{\phi }}_{\mathrm{1}}\cdot \sin ({\mathit{\eta }}_{i+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{1}})\cdot \cos \mathit{\alpha }\right)\text{(31)}& {\displaystyle }{\mathit{\lambda }}_{i+\mathrm{1}}& {\displaystyle }={\mathit{\lambda }}_{\mathrm{1}}+\arctan \left({\displaystyle \frac{\sin \mathit{\alpha }}{\cot ({\mathit{\eta }}_{i+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{1}})\cdot \cos {\mathit{\phi }}_{\mathrm{1}}-\sin {\mathit{\phi }}_{\mathrm{1}}\cdot \cos \mathit{\alpha }}}\right),\end{array}$$

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 ray-tracing 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):

$$\begin{array}{}\text{(32)}& g_{\mathrm{bend}}\left[{}^{\circ }\right]={\displaystyle \frac{\mathrm{0.02}\cdot {\exp }^{\frac{-h}{\mathrm{6000}}}}{\tan {\mathit{\epsilon }}_k}},\end{array}$$

where h is replaced by h_t, the height of the voxel top layer. After convergence of the ray-tracing 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 ray-traced 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 Ray-traced signal path differences (b) caused by differences in the a priori refractivity field (a).

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 3 Convergence behavior (a) and ray-traced signal path differences after convergence (b). All iterations are based on the same first guess (standard atmosphere or ALARO numerical weather model data) but differ with respect to the refractivity field used for reconstruction of the bended signal paths.

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.

Figure 4 Point error at the voxel model top (h=13.6 km) caused by the bending model of Hobiger et al. (2008) – computed on a global $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$ grid over the period of one year, 2014, by comparison with ray-traced bending angles based on ECMWF analysis data.

4.2.2 The empirical ray-bending 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 ray-traced 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.

$$\begin{array}{}\text{(33)}& N\left(f\right)={\mathrm{10}}^{\mathrm{6}}\cdot {\displaystyle \frac{-\mathrm{40.2993}\cdot N_{\mathrm{e}}}{f^{\mathrm{2}}}}\end{array}$$

By making use of Eq. (33), the obtained electron density profiles were converted into profiles of refractivity (N), assuming signal frequency f₁=1575.42 MHz (GPS L1) and f₂=1227.60 MHz (GPS L2). Figure 5 shows the obtained vertical profiles of ionospheric refractivity as an example for frequency f₁. The higher the signal frequency f, the lower the phase velocity through the ionosphere and the less its refraction.

Figure 5 Profiles of ionospheric refractivity N(f) assuming signal frequency f=1575.42 MHz.

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 straight-line 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.

5 Impact of atmospheric bending on the tomography solution

In the following, the differences between straight-line 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 straight-line and bended ray tracing. The only difference is that in the case of straight-line 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.

Figure 6 Additional ray path caused by the straight-line assumption (a), the resulting drying effect due to the additional ray path (b) and the resulting drying effect caused by the fact that the straight-line ray travels through lower atmospheric levels than the true bended ray (c).

The additional ray path decreases rapidly with increasing elevation angle. Thus, a mixed ray-tracing 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 straight-line 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 ray-tracing approaches is negligible (be aware that only the entire distance through all voxels is comparable for both ray-tracing 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 straight-line geometry. Here, it is distinguished between the drying effect caused by the additional ray path (dNw₁; panel b) and the drying effect caused by the fact that the straight line travels through lower layers of the voxel model (dNw₂; panel c). Both effects were assessed as follows:

$$\begin{array}{}\text{(34)}& {\displaystyle }& {\displaystyle }\mathrm{d}N{w}_{\mathrm{1}}={\text{SWD}}_{\mathrm{b}}\cdot (d_{k,\mathrm{s}}-d_{k,\mathrm{b}})\text{(35)}& {\displaystyle }& {\displaystyle }\mathrm{d}N{w}_{\mathrm{2}}=({\text{SWD}}_{\mathrm{b}}-{\text{SWD}}_{\mathrm{s}})\cdot d_{k,\mathrm{s}},\end{array}$$

where SWD_b and SWD_s are the slant wet delays obtained by ray tracing through ALARO model data along the bended and the straight-line 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.

Table 1 Summary of GNSS tomography test case settings.

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 straight-line solution are traversed due to bending.

Figure 7b shows the differences in wet refractivity between the first two iterations of the mixed ray-tracing 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.

Figure 7 Error in wet refractivity caused by the straight-line assumption (a) and the differences in wet refractivity between first and second iterations (b). Here, voxel number 1 is dedicated to the southwest corner and number 20 to the northeast corner of the voxel model. For visualization, a bilinear interpolation method was applied between the grid points. Analyzed period: 4 May 2013 at 15:00 UTC.

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 ray-tracing approach against straight-line ray tracing, the tomography-derived 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 tomography-derived 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 ray-tracing 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.

Figure 8 Differences in wet refractivity between radiosonde, ALARO and the two tomography solutions based on straight-line (blue) and bended ray tracing (red) as an example for 1 May 2013 at 03:00 UTC (a) and 31 May at 03:00 UTC (b).

6 Conclusions

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 2-D piecewise linear ray-tracing 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 ray-tracing approach. Due to significant optimization, the mixed ray-tracing 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 ray-tracing approach is also applicable in near-real time and therefore well suited for operational purposes.

Appendix A: Unmodeled bending effects in the Vienna hydrostatic mapping function

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:

$$\begin{array}{}\text{(A1)}& a_{\mathrm{h}}=-{\displaystyle \frac{m{f}_{\mathrm{h}}\left(\mathit{\epsilon }\right)\cdot \sin \mathit{\epsilon }-\mathrm{1}}{\frac{m{f}_{\mathrm{h}}\left(\mathit{\epsilon }\right)}{\sin \mathit{\epsilon }+\frac{b_{\mathrm{h}}}{\sin \mathit{\epsilon }+c_{\mathrm{h}}}}-\frac{\mathrm{1}}{\mathrm{1}+\frac{b_{\mathrm{h}}}{\mathrm{1}+c_{\mathrm{h}}}}}},\end{array}$$

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(ε,α), ray-traced slant hydrostatic delays were compared with mapped slant hydrostatic delays as follows:

$$\begin{array}{ll}{\displaystyle }\mathrm{d}{g}_{\mathrm{bend}}(\mathit{\epsilon },\mathit{\alpha })\left[m\right]& {\displaystyle }=\text{ZHD}\left[m\right]\cdot m{f}_{\mathrm{h}}\left(\mathit{\epsilon }\right)-\text{ZHD}\left[m\right]\cdot m{f}_{\mathrm{h}\mathrm{0}}\left(\mathit{\epsilon }\right)\\ \text{(A2)}& {\displaystyle }& {\displaystyle }-g_{\mathrm{bend}}(\mathit{\epsilon },\mathit{\alpha })\left[m\right],\end{array}$$

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₀(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).

Figure A1 The unmodeled geometric bending effect in VMF1 hydrostatic mapping function (dg_bend), as an example for VLBI sites in Fortaleza, Brazil, and Wettzell, Germany. Analyzed period: January–February 2014.

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.