Unconstrained , robust Kalman filtering for GNSS troposphere tomography

Introduction Conclusions References


Introduction
The distribution and dynamics of water vapour (WV) and its closely associated meteorological phenomena, such as long persistent rainfalls, tropical cyclones, mid-latitude cyclonic storms and thunder storms (Ahrens and Samson, 2010), are on-going challenges for synoptic meteorology.These severe weather phenomena could potentially Figures in weather forecasting.Improving the understanding of WV distribution is important (Le Marshall et al., 2012), not only for meteorology, but also for an effective usage of the GNSS technology for precise positioning.For example, tropospheric effects are one of the important atmospheric errors that need to be removed or mitigated in both high accuracy differential positioning and precise point positioning (Wielgosz et al., 2012).
GNSS tomography is based on the inverse Radon transform theory and it has been intensively investigated by a number of research groups and universities across the globe (Bender et al., 2011;Perler et al., 2011;Brenot et al., 2012;Flores et al., 2000).A special working group on the troposphere tomography models' integration has been set up recently by the International Association of Geodesy (i.e.IAG WG4.3) to coordinate this IAG initiative (Rohm et al., 2012).In addition the GNSS tomography is an extended service of the Ground Based Augmentation System (GBAS).The standard approach to establish GNSS tomography models is to divide the troposphere into a 3-D voxel structure over the area of the GNSS CORS network coverage.The intercepted distance of a GPS ray passing through the voxel of concern is used in the design matrix for the calculation of the refractivity (Fig. 1).The design matrix is then inverted to estimate unknown wet refractivity values (Flores et al., 2000).
The major challenge of the tomography is to obtain a stable solution, in the presence of ill-conditioness and ill-posedness of the inverse process.However, these can be broadly divided into two categories: i.e. (1) increasing the number of pseudo observations (and decreasing the condition number) by adding horizontal and/or top layer constraints (Rohm and Bosy 2009;Perler et al., 2011;Hirahara, 2000;Flores et al., 2000;Bender et al., 2011); (2) extending the observation system with additional observations in zenith direction, point observations, and radiosonde and radiometer profiles (Perler et al., 2011).As a consequence all solutions are similar in the sense that refractivity values are given for all voxels in the model even though not all are intercepted by GNSS signals (Perler et al., 2011).
In addition to the studies of tomography observation system, the 3-D model structure was investigated (e.g. by ETH Zurich tomography research group).Perler et al. (2011) Introduction

Conclusions References
Tables Figures

Back Close
Full recently showed that it is feasible to indirectly calculate the coefficients of a wet refractivity trilinear spline function instead of the wet refractivity inside each voxel.The most common inversion technique applied in GNSS tomography is based on a Singular Value Decomposition (SVD), technique that allows only for complete A matrix inversion (Rohm and Bosy, 2009;Perler et al., 2011).A slightly different approach was proposed by Bender et al. (2011) where an algebraic reconstruction technique is applied to estimate refractivity in an iterative way.Each iteration step updates wet refractivity only in voxels that are actually intercepted by the GNSS signals.
In this study, to overcome the ill-conditioness of the inverse tomography problem without applying explicit constraints, the SVD method proposed by Xu (1998) and Lynch (2005) is used.The novelty in this approach is that the information provided in the observation matrix is used exclusively and singular values in the design matrix are sensibly selected.This paper aims to extend the previous research (e.g.Rohm and Bosy, 2011;Bosy et al., 2012;Rohm, 2012Rohm, , 2013) ) by a new concept of a robust Kalman filter.Unlike all predecessors, the solution of the tomography model presented in this paper is not affected by implicit constraints (i.e.no horizontal and no vertical constraints are applied), and it may not rely on additional observations (i.e.no need for NWP observations).It delivers solutions only for voxels that are intercepted by GNSS signals with a full variance-covariance matrix.The robust Kalman filter allows for variations of the refractivity field in time and prevents the noise propagation in outputs.This technique is discussed in Sect. 4.Moreover, GNSS signal delay and the tomography model structure is presented in Sects. 2 and 3, respectively.Case study investigation using real and simulated data is performed in order to assess the quality and effectiveness of the new approach proposed.Conclusions and summary is given in the closing section of the paper.Introduction

Conclusions References
Tables Figures

Back Close
Full

GNSS signal delay
GNSS carrier frequencies reside in a microwave spectrum reserved for navigation (L-band, 1-2 GHz) (Hofmann-Wellenhof et al., 2008).Such spectrum's allocation is designed to minimise signal attenuation in the atmosphere, and hence allows for all weather operation of the system.The microwave signal's refractivity in the neutral atmosphere is given as (Thayer, 1974): where p is the pressure (hPa), e is the water vapour partial pressure (hPa), T is the temperature (K), k i (i = 1, 2, 3) are known empirical coefficients (Kleijer, 2004), Z −1 d is an inverse compressibility factor for dry air and Z −1 v is an inverse compressibility factor for water vapour, respectively (both values are assumed to be 1 in this study).The analysis in this paper is focused only on the phase speed changes of the GNSS signals' propagation (delay) in the neutral atmosphere, hence signal bending is neglected, and no effects on the signal's energy are considered.
Tomography processing assumes that signal is sufficiently modelled as a straight line between the satellite and the receiver and the Slant Total Delay (STD) in neutral atmosphere is given by the following equation (Kleijer, 2004): where, SHD is the slant dry delay, and SWD is the slant wet delay.Usually tomography models utilise the SWD to reconstruct the water vapour distribution.

Model structure
The tomography technique is founded on the theory of the Radon transform and its inverse (Kak and Slaney, 2001).In principle any function's integral along the path line, Introduction

Conclusions References
Tables Figures

Back Close
Full executed along an infinite number of lines could be converted into the distribution of the medium affecting the signal path.According to the Radon principle (Kak and Slaney, 2001), a single scanning ray SWD n from a satellite to a receiver, is given as: where, N vm is the water vapour refractivity in the voxel m and d mn is the intercepted distance in the voxel m of signal n.
The full functional model of the tomography solution in a matrix form is given as The observations in Eq. (4) (SWDs) are appended with an a priori value of refractivity N vapriori .The design matrix consists of three blocks; A inner is the matrix containing distances in the inner model and A outer is the equivalent matrix for the outer model, A a priori is the matrix containing value 1 when there are external observations and value 0 when there are no external observations in the voxel.The unknowns N v and No v are the wet refractivities in the inner model and the water vapour refractivity in the outer model, respectively.The general equation for tomography using relation (4) in the matrix notation is given as: Equation ( 5) is ill-conditioned (no constraints, with a limited number of observations) and ill-posed (a limited number of observations).Therefore, an effective inversion of matrix A is a central problem of all GNSS tomography applications.The unconstrained tomography solution studies (Rohm, 2012) show that it is feasible to obtain precise wet refractivity exclusively from the GNSS SWD observations using 9138 Introduction

Conclusions References
Tables Figures

Back Close
Full stacked observations from at least ten 1 h consecutive epochs.However, ten hours is to long time period to be represented by the single value of refractivity, it is therefore convenient to use the robust Kalman filter as a dynamic model of troposphere.This is the focus of the next section.

Kalman filter application
A classic Kalman filter formulation follows notation given in Grewal et al. (2001) where observations and process are separated.In this study, the process is set to be a wet refractivity field N v k,k+1 with the time evolution given as a linear dynamic system (Yang, 2010): Where Φ k,k+1 is a state transition matrix (in this study it is an identity matrix Φ k,k+1 = I).
The w k parameter is the noise with the characteristics of mean E (w k ) = 0 and covariance E (w k w T k ) = Q k , which is called the dynamic disturbance noise matrix.The observation linear model for epoch k is given in Eq. ( 6): where SWD k are uncorrelated normally distributed observations, and measurement noise ν k has mean However, in a robust Kalman filter observations are assumed to be of a normal distribution contaminated with outliers (Yang, 2010), therefore The prediction step of Kalman filtering is given as a set of equations:

Conclusions References
Tables Figures

Back Close
Full where Nv k (−) and Nv k−1 (+) is the predicted and the corrected estimates of wet refractivity in the voxels of GNSS tomography model.The matrices P k (−) and P k−1 (+) are the prediction and the correction P k−1 (+) of the covariance matrix of the estimated state.The Kalman gain matrix K is: the covariance matrix R R k of the Robust Kalman Filter is calculated using the following equation: where The parameter e i is a posteriori residual value calculated via: Where, c = 1.5 is a scaling factor and σ is a reference variance (usually 1 mm).
Usually, a robust Kalman filter is applied to observations contaminated with outliers, to minimise or remove their impact by increasing the selected observation's covariance in the estimation process.Therefore the process of estimating R R k is iterative and might need to be repeated several times.In the paper by Koch and Yang (1998) downweighting is applied on the parameters, however in this study to be consistent with the previously developed SVD methodology (Rohm and Bosy, 2011;Rohm, 2013Rohm, , 2012) ) Introduction

Conclusions References
Tables Figures

Back Close
Full the downweighting of the parameters is not used.In this paper the structure of a design matrix A is evaluated to reveal and remove linearly dependent observations (in a numeric sense).This technique reduces matrix A condition number and improves inversion stability.The design matrix A filtering process follows the methodology developed by Rohm (2013), based on the work of (Xu, 1998;Lynch, 2005).In essence, matrix A is decomposed into three vectors (Strang and Borre, 1997): Where, U, V are a set of orthonormal bases and Σ is a set of singular values (s k ).
A condition number of any matrix (including matrix A) is calculated as (Anderson et al., 1999): where k is the rank of matrix A.
Previous investigation by Rohm (2013) shows that the uncertainty expressed as a covariance of wet refractivity is linked with the singular values on the diagonal of matrix Σ.The smaller the singular value considered in the design matrix (A) inversion (Eq.10) the stronger the amplification of observation uncertainty R R k .The last, useful singular value (k) considered in the processing is found using functional analysis of singular values function (Xu, 1998;Lynch, 2005;Hansen and O'Leary, 1993), this method to improve the matrix condition number is named Truncated Singular Values Decomposition (TSVD).The new filtered Ã matrix is obtained by composing the A matrix back Introduction

Conclusions References
Tables Figures

Back Close
Full according to the equation: To reflect the changes in matrix A observation matrix SWD has to be converted to the matrix S WD to eliminate the linearly dependent observations.The identification of linearly dependent rows (f ) is based on comparing rows from matrices A and Ã according to the following equation: Consequently, the Kalman filter sequence as shown in Koch and Yang (1998) for filtering observations will be transformed to the following sequence: S WD = TSVD(SWD) (20) The robust estimation process of wet refractivity starts with the filtering of matrix A to produce Ã (Eq.19), as well as truncating observations SWD to obtain S WD (Eq.20).9142 Introduction

Conclusions References
Tables Figures

Back Close
Full Initial estimates of parameters Nv k are calculated via Eq.( 21).Residuals e are derived then, which form the base for downweighting of the outliers (Eqs.22, 12, and 13), and calculation of R R k matrix.The following step comprises of a Kalman gain K derivation (Eq.24).Equations ( 21)-( 24) are repeated several times to remove outliers from observations.This operation is followed by an update step (Eq.25), and the propagation of covariance of parameters to the next epoch (Eq.27).
The matrix Q k adds a noise to each voxel in a covariance matrix P k−1 .The amount of noise in each element of matrix Q k (q m,i for the inner domain and, q o,i for the outer domain) is driven by three factors: (1) location in the model h; (2) time since last update t; (3) location in inner or outer model, according to the following formula: where In this study, uncertainty parameters σ

NWP model
The meteorological data covers ACCESS-R model outputs (analysis run) with the time resolution of 6 h and spatial outline covering Australia and a ∼ 35 • buffer zone.The model is based on the UK Meteorological Office Unified Model, and a number of data sources is used to produce forecasts (e.g.COSMIC, AIRIS, SYNOP) (Le Marshall et al., 2012).The model in the horizontal plane contains 320 × 220 nodes, with the grid spacing of 0.375 • (∼ 37.5 km) and the model utilises in the vertical direction terrain following hybrid (pressure/height) coordinates with 50 levels.This study, from all possible NWP model parameters, considers only pressure, temperature and WV partial pressure (given as a mixing ratio) (Fig. 3).

Common data
Two ZTD datasets, with the same time, terrain and identical receiver network settings are prepared, one is simulated and the other one is real data.cuboids (in this study called voxels), from ground level to 10 km with increasing layer thickness (varying between 500 to 1700 m), whereas the last voxel spans from 10 to 15 km and has 5 km thickness.The horizontal plane of the tomography model consists of 6 voxels in the north direction and 12 voxels in the east direction (to accommodate the model outline to the GPSnet shape).As a consequence of the vertical and horizontal settings, the size of each voxel is approximately 75 km × 45 km × 0.5 ∼ 5 km, which roughly represents the average inter station distance.
Pressure, temperature and water vapour produced by ACCESS-R model ( 1) are utilised to separate wet and dry delays, and (2) act as a reference value for tomography model reliability investigation.To separate wet and dry delays, pressure values from 8 NWP model nodes (4 below, 4 above) surrounding GNSS stations are interpolated to the antenna reference point.The second use of NWP model outputs is realised interpolating water vapour partial pressure and temperature values from NWP model nodes to the tomography's voxel centre points (Bosy et al., 2012(Bosy et al., , 2010)).Alternatively the pressure and temperature values are obtained from the Global Pressure and Temperature (GPT) model (Boehm et al., 2007) and water vapour pressure is calculated from UNB3m (Leandro et al., 2008).

Simulation data
The first dataset comprises wet refractivities N R derived from NWP temperature and WV interpolated to the centre of each voxel of the tomography model.Then using simple analytical raytracing (Rohm, 2013;Rohm and Bosy, 2009) the signal intercepting distances in each voxel of the inner A inner and the outer models A outer (Eq.4) are calculated, along with observations SWD S (Eq.7).Uncertainty of the raytracing through the model of unknown precision (there is no impartial measure of weather model parameters precision) remains unidentified.Therefore all elements of the covariance matrix R S are equal in weighting matrix.The first observations dataset is a simulation of real observations based on the NWP model data; it also constitutes the reference data collection.

Conclusions References
Tables Figures

Back Close
Full

Real data
The second dataset consists of the real observations from a GNSS network processed with the Bernese GPS Software 5.0 (Dach et al., 2007).The wide/narrow lane (L5/L3) GNSS processing procedure are applied (Dach et al., 2007).The ambiguities are solved with the wide-lane L5 (98 % success rate)/narrow-lane L3 (90 % success rate) strategy.Final coordinates are estimated with the minimum constraint conditions imposed on the translation parameters of coordinates and velocities of IGS reference stations (MOBS, HOB2, STR1 and CEDU).The mean accuracy of the solutions in the horizontal directions and in the vertical direction, based on repeatability score are 1.5 mm and 4 mm, respectively.
The troposphere estimates in 30 min resolution are obtained in the next step by fixing the translation parameters of the network and pre-eliminating the velocities as well as the coordinates from weekly solution.In Bernese GPS software the standard approach to estimate the ZTD (Dach et al., 2007) is used.The ZTD parameters are estimated as corrections to an a priori standard atmosphere model using piecewise linear functions.
The Niell (1996) mapping function is used to map the observed delays to the vertical direction.The atmospheric gradients are estimated.This is not an optimal configuration (state of the art mapping functions are not used) to estimate troposphere parameters.However it is a common approach used in Bernese 5.0 GPS Software.Output TRO and TRP files comprise the most important part of the second dataset.This set also includes: pressure parameters interpolated to the antenna heights from the NWP model, final station coordinates as well as precise orbits from IGS.The dry part is subtracted from the total delay based on the Saastamoinen model of dry delay (Saastamoinen, 1972) with pressure values from the NWP model.
In this study SWD G are calculated using either zenith part of the delay or the zenith delay and horizontal gradients, in either case the DD residuals are not considered.
The observations in the slant direction for the wet part of the delay SWD G (Eq. 7) are determined by applying the wet Niell mapping function (Niell, 1996).Therefore the Introduction

Conclusions References
Tables Figures

Back Close
Full SWD G s are not uncorrelated and the mapping function used to map the delay from zenith to slant direction contains implicit information on the vertical distribution of WV.Using simple analytical raytracing (Rohm, 2013;Rohm and Bosy, 2009) the signal interception distance in each voxel of the inner and outer models (elements of matrix A G (Eq. 7)) are calculated.The uncertainty measures R G are based on the estimated error of particular ZTD value by applying the law of error propagation (Rohm, 2012).

Tomography processing results discussion
A number of test runs of the tomography model are performed to precisely assess the impact of particular methodological improvements.The following major groups of settings are adopted in regards to observations: (1) simulated observations (M) (with different levels of noise); (2) real observations with gradients (R); and finally (3) real observations without gradients (Z).Alternatively experiments may be grouped together in relation to the a priori models adopted, the following settings are considered: ( 1 3) NWP derived outer model values for all epochs and UNB3m and GPT derived inner model for all epochs (G2).To assess the impact of the innovative robust Kalman filter processing procedure, the following three levels validation is adopted: (1) firstly all equations related to the observation selection criteria are applied including: SWD observation removal (Eq.20) (S); reconditioning of matrix A (Eq. 18) (A) and downweighting of selected observations (Eqs.11, 23, and 24) (D); (2) secondly the downweighting scheme is not applied (O) but observation removal (S) and reconditioning (A) is, (3) thirdly no robust improvements of the Kalman filter are considered (OOO), so the filter runs like a classic Kalman filter.
The experiment setup is shown in Table 1.In total 18 different settings are investigated, and the most significant results are presented in Table 2.The solutions for all 325 epochs are depicted on Figs.4-7.Figures 4 and 6 present the mean bias of the 9147 Introduction

Conclusions References
Tables Figures

Back Close
Full tested solution against reference data, whereas Figs. 5 and 7 show the mean standard deviation of the tested solution against reference data.The major outcomes of this experiment are summarised below.
The most important conclusion drawn from the set of experiments is that the a priori value N vapriori for inner model (N, G1, and G2) is the main factor in all processing schemes.Figures 5-7 show that whenever the a priori value for all epochs and all voxels are set, even with a simple deterministic model such as UNB3m and GPT the quality of the reconstruction is much higher than those in all other cases.The values in Table 2 show also higher accuracy of the G2 solution.Introducing initial values into the tomography system (Eq.4) can effectively stabilise the tomography solution.In this study the initial wet refractivity field is a function of the day of the year, latitude, longitude and altitude.The variance of the a priori observations in Eq. ( 17) is a set to rather large values (i.e. 30 mm km −1 ).Hence, the results show that the quality of the a priori observations is not an issue for tomography reconstruction.
The second most important outcome of this research is that the robust filtering helps to reduce noise in outputs.Clearly, the dashed red line showing the standard deviation of the real data solution on Fig. 5 (RNSAD) is much lower than the one with partial robust algorithms (RNSAO) and no robust procedures in place (RNOOO).However the difference between the last two is not significant which in turn means that the most significant improvement in real time data processing is due to the downweighting not reconditioning.The same effect is visible when the processing covers the real observations with a large number of a priori data (RG2SAD -black dashed line on Figs. 4  and 5).However the effect is not strong, at least in the solution scatter.In terms of systematic errors, the mean difference is effectively removed by the robust algorithm (RG2SAD) as depicted on Fig. 4.
The third conclusion is that the best achievable performance using this tomography model and simulated observations (without noise) is 4.2 mm km −1 (Table 2) (M1G2SAD shown on the Fig. 7).However, introduction of the realistic noise and bias to the observations (0.025 m of the random noise and 0.007 m of the bias) results in the tomography Introduction

Conclusions References
Tables Figures

Back Close
Full quality degradation to 6.8 mm km −1 (Table 2, M2NSAD).In comparison, the quality of the tomography retrieval based on real data (Table 2, RG2SAD) is 6.7 mm km −1 , (black dashed line in Fig. 7).This suggests that both solutions (real-RG2SAD and simulated-M2NSAD) are equally accurate and hence all real data outliers are effectively filtered out.
Another important issue clearly visible on Figs.4-7 is that there is not much difference between tomography solution fed by the observations with and without gradients (Table 2 RG2SAD 6.5 mm km −1 standard deviation, ZG2SAD 6.7 mm km −1 standard deviation).The same level of bias has been also observed for both type of measurements (Table 2 RG2SAD and ZG2SAD 0.4 mm km −1 bias).Therefore using gradients in the signal delay modelling does not improve solution.
Many authors (Bender et al., 2011;Manning et al., 2012;Perler et al., 2011;Rohm, 2012) report that the tomography quality varies between 4 to 10 mm km −1 and is lower for the bottom level of troposphere and increase with high until the amount of water vapour is lower than the sensitivity of the method.This suggests that the effectiveness of the tomography method in resolving vertical structure of troposphere needs further investigation.In this study to validate whenever this method has some advantage over a deterministic model, we simply subtract GPT and UNB3m wet refractivity from NWP derived refractivities and calculate the statistics such as standard deviation and bias.The results in Table 2 (UNB3m, GPT) shows that the standard deviation is slightly higher than that of tomography model (7.2 mm km −1 ), but the bias is much higher (−3.6 mm km −1 ).Hence, the tomography processing has advantage over the deterministic models; the question remains if the level of the obtained quality is satisfactory for meteorological and positioning applications.

Conclusion
In this paper, the new GNSS tomography model TOMO2 is presented.This model employs a robust Kalman filter to solve the unconstrained (in the implicit sense) to-9149 Introduction

Conclusions References
Tables Figures

Back Close
Full mography problem.This study demonstrates that the real slant wet delay dataset is affected by noise and outliers and the estimated zenith delay uncertainties are overly optimistic.Therefore the real GNSS data requires advanced processing beyond the ordinary Kalman filter.In this study, both the robust Kalman filter and a truncation of the design matrix (with TSVD method) are investigated to limit the noise impact on the model updates.The estimations of wet refractivities and their associated uncertainties in the troposphere above a network of GNSS receivers for selected voxels can be determined through these methods.Results show that the STD discrepancy between the reference wet refractivity and the tomography model outputs is in the order of 6.5 mm km −1 (or ppm), which is the equivalent to 1 hPa of WV (Rohm, 2012).The results are in a good agreement with GNSS tomography simulation studies with an intermediate level of noise 5.2 mm km −1 (or ppm) (Bosy et al., 2012).The most important contribution of this paper is an effective GNSS tomography reconstruction without using implicit constraints, which allows for a quicker tomography model to be established in a response to the changing environment conditions.The constraint free approach investigated in this research produces more realistic wet refractivity uncertainties that are unbiased by inner constraints.The results presented in this paper shows current level of the quality achievable for tomography reconstruction.Further discussion with meteorological community is needed to investigate an efficient way to assimilate the GNSS tomography outputs into NWP models.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2h,i and σ 2 o
are calculated from NWP model outputs, rescaled to the tomography model voxels.The performance analysis takes into account vertical variability of the wet refractivity as well as time autocorrelations of these parameters.In a more general case (without access to NWP data) the uncertainty parameters may be derived from climatological data.5 Case study To demonstrate the capability of this new GNSS tomography model TOMO2, a case study based on simulated and real data is performed and the results are validated against NWP model outputs.Discussion Paper | Discussion Paper | Discussion Paper | The same tomographic model setup is used.The time span covered by this case study is limited to 325 epochs of ZTD estimation between 3 March 2010 and 9 March 2010, whereby each ZTD epoch covers an interval of 30 min.GNSS observations from 45 stations were taken during the development, transition and dissipation of a heavy hail storm(Choy et al., 2011).The GNSS network employed in this study (GPSnet) (Fig.3) is owned and operated by the Victorian Government Department of Sustainability and Environment (Victoria, Australia).Standard GPSnet stations are equipped with Trimble NetR5 receivers and high quality antennas (mostly TRM55971.00).A few International GNSS Service (IGS) stations were also processed, where their receivers/antennas are from other manufacturers (e.g.Leica, Ashtech).The inter receiver distance is roughly 50 km in the investigated case.Terrain undulation, especially in the east mountainous part of Victoria, is in favour of tomography because of large height differences, some receivers are located on the mountains' peaks (e.g.MTBU -Mount Buller 1600 m).The troposphere in the vertical direction, above a network of the GNSS receivers, is divided into a number of Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) NWP derived outer model values for all epochs and NWP derived inner model values for the first value (as an initial value) (N); (2) NWP derived outer model values for all epochs and UNB3m and GPT derived inner model for the first epoch (G1); and finally ( Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Strang, G. and Borre, K.: Linear Algebra, Geodesy, and GPS, Wellesley Cambridge Pr., 1997.Wielgosz, P., Paziewski, J., Krankowski, A., Kroszczyñski, K., and Figurski, M.: Results of the application of tropospheric corrections from different troposphere models for precise GPS rapid static positioning, Acta Geophys., 60, 1236-1257, 2012.Xu, P.: Truncated SVD methods for discrete linear ill-posed problems, Geophys.J. Int., 135Discussion Paper | Discussion Paper | Discussion Paper | ) and inner model (initial epoch) (N), in Kalman filter A matrix reconditioning (A), SWD removal and downweighting (D) has been used.OBSERVATIONS SWD type A priori type R -Real observations with gradients N -NWP outer (all epochs) and inner Z -Real observations without gradients (first epoch) M1 -Simulated observations without noise G1 -GPT + UNB3m outer(all epochs) M2 -Simulated observations with and inner (first epoch) realistic noise G2 -NWP outer (all epochs) and GPT + UNB3m outer inner (all epochs) The signal from satellite (modelled as a straight line BA) intersects with the ontal plane given by three points (0, 1 and 2) in the pierce point P. The distance)) m between the pierce point P and the GNSS station (A) is an element of design matrix e tomography processing

Fig. 1 .
Fig. 1.The signal from satellite (modelled as a straight line BA) intersects with the horizontal plane given by three points (0, 1 and 2) in the pierce point P. The distance (d R S (m)) between the pierce point P and the GNSS station (A) is an element of design matrix in the tomography processing

Fig. 6 .Figure 6 .Figure 7 .
Fig. 6.Total bias of estimated parameters over 325 epochs, for number of a priori modes (N -NWP inner and outer, G2 -UNB3m and GPT inner NWP outer) and types of observations (R -real with gradients, Z -real without gradients, M1 -simulated without noise, M2 -simulated with realistic noise).

Fig. 7 .
Fig. 7. Standard deviation of estimated parameters over 325 epochs, for number of a priori modes (N -NWP inner and outer, G2 -UNB3m and GPT inner, NWP outer) and types of observations (R -real with gradients, Z -real without gradients, M1 -simulated without noise, M2 -simulated with realistic noise).

Table 1 .
The list of tested tomographic solutions.The naming conventions explained in two bottom rows.