A troposphere tomography method considering the weighting of input information

Troposphere tomography measurement using a global navigation satellite system (GNSS) generally consists of several types of input information including the observation equation, horizontal constraint equation, vertical constraint equation, and a priori constraint equation. The reasonable weightings of input information are a prerequisite for ensuring the reliability of the adjustment of the parameters. This forms the focus of this research, which tries to determine the weightings, including the observations, for the same type of equation and the optimal weightings for different types of equations. The optimal weightings of the proposed method are realized on the basis of the stable equilibrium relationship between different types of a posteriori unit weight variances, which are capable of adaptively adjusting the weightings for different types of equations and enables the ratio between the two arbitrary a posteriori unit weight variances to tend to unity. A troposphere tomography experiment, which was used to consider these weightings, was implemented using global positioning system (GPS) data from the Hong Kong Satellite Positioning Reference Station Network (SatRef). Numerical results show the applicability and stability of the proposed method for GPS troposphere tomography assessment under different weather conditions. In addition, the root mean square (RMS) error in the water vapor density differences between tomography-radiosonde and tomography-ECMWF (European Centre for Medium-Range Weather Forecasts) are 0.91 and 1.63 g m−3, respectively, over a 21-day test.


Introduction
Tomography is one of the most cost-effective means of obtaining the spatio-temporal three-dimensional distribution of integral measurements within the study region.Since the use of tomography was first publicized (Bramlet, 1978), it has developed as a powerful and practical tool to retrieve information of interest in many application areas, such as geology (Bourjot and Romanowicz, 1992), earthquakes (Kissling et al., 1994) and ionosphere and troposphere modeling (Hajj et al., 1994;Rius et al., 1997;Braun et al., 1999;Flores et al., 2000), wind (Gao et al., 1999).
Global navigation satellite system (GNSS) troposphere tomography usually refers to the three-dimensional reconstruction of water vapor information based on the GNSS signal rays that cross the tomographic modeling area (Flores et al., 2000).The tomographic area is usually discretized into many voxels in the latitudinal, longitudinal, and vertical directions, by which integral observations can be used to calculate the wet refractivity value and water vapor density value of each voxel (Hirahara, 2000;Champollion et al., 2005;Rohm, 2013;Rohm et al., 2014); however, for voxels without signals crossing within the area of interest, which is mainly caused by the unfavorable geometry of ground-based GNSS receivers and the constellation of GNSS satellites, the parameter information value is unavailable.This drawback also leads to a rank deficiency in the tomographic observation modeling.To overcome this ill-posed problem, some constraint conditions were mentioned in previous studies for superposition onto the tomographic modeling (Flores et al., 2000;Hirahara, 2000;Troller, et al., 2002;Skone and Hoyle, 2005;Nilsson and Gradinarsky, 2006;Bender and Raabe, 2007;Perler et al., 2011;Brenot et al., 2014).The most com-Q.Zhao et al.: A troposphere tomography method monly used constraints include a horizontal constraint, vertical constraint, and a priori constraint.
Regardless of the type of constraint condition imposed on the tomographic modeling, these different types of constraint information need to be considered as the independent input information, which means that (1) the reasonable weightings among the observations for the same type of input information should be considered, (2) and the reasonable unit weight variances for different types of input equations should also be given before the tomographic modeling resolution.Flores et al. (2000) first used the weightings of three constraints (horizontal, vertical, and boundary constraints) to a constant k and increased the value of k, which made the minimum eigenvalue of the coefficient matrix form of the normal equation larger than the threshold of 8.1 in their tomographic modeling.Song (2004) proposed a posterior variance component estimation method that also considers the robust estimation of parameters and variance components, and the termination condition was considered to have been met when the values of different posterior unit weight variances are equal.Bosy et al. (2012) and Xia et al., (2013) both introduced the same weight value for various equations.Chen and Liu (2014) tuned the weight on the horizontal constraint following the method of Flores et al. (2000), while the weight on the vertical constraint and observation equation are the same.Guo et al. (2016) proposed an optimal weighting method to determine the reasonable weights for various types of equations on the basis of being statistically equal, while the weightings among observations for the same type of input information are the same.
The aforementioned methods of selecting the weights of different kinds of input equations usually include a determination of the weights by using an iterative method until the final unit weight variances are equal on the basis of statistics or mathematics; however, some points still need to be discussed as follows: (1) as for the same type of equation, how to determine the weightings among observations was seldom mentioned in the results of previous published studies, and this is important when trying to obtain an accurate tomographic result; (2) the scheme used to select equal weights was unable to embody a reasonable correlation among different types of equations; (3) even although iterative methods were applied to calculate the weights, it may be difficult to obtain a reasonable result due to the various practical conditions in regions subject to tomography; (4) and it may be unreasonable and time-consuming merely from the perspective of mathematics or statistics when determining the optimal weights for various equations.Consequently, in this work, (1) the determination of weights for the same type of input information is presented based on the characteristics of their practical relationship, (2) and for the different types of equations, a weighting algorithm is proposed based on the concept of a stable equilibrium relationship for their a posteriori unit weight variances.This paper is organized as follows: the basic theory of GNSS troposphere tomography modeling is introduced in Sect. 2. The determination of the weightings is described in Sect.3. In Sect.4, the test results from tomographic experiments under various weather conditions (21-day duration) are validated by comparison with radiosonde and ECMWF data.Conclusions and recommendations for future research are provided in Sect. 5.

Theory of GNSS troposphere tomography modeling
Previous studies have shown that two kinds of input information derived from GNSS observations are usually considered to build the tropospheric observation equation.On the one hand, the slant wet delay (SWD) was introduced to obtain the troposphere wet refractivity information (Flores et al., 2000;Skone and Hoyle, 2005;Troller et al., 2006;Rohm and Bosy, 2009;Notarpietro et al., 2011;Guo et al., 2016).On the other hand, the slant-integrated water vapor (SIWV) was applied to reconstruct the three-dimensional atmospheric water vapor field (Champollion et al., 2005;Bi et al., 2006;Xia et al., 2013;Jiang et al., 2014;Yao et al., 2016).In our project, we are trying to assimilate the three-dimensional water vapor data into a weather research and forecasting model (WRF), and this paper thus focuses on GNSS water vapor field reconstruction using SIWV data.

Estimation of SIWV
Generally, two types of processing modes are adopted for GNSS observations when retrieving tropospheric parameters: the first is the un-differencing mode, which only needs one GNSS receiver to estimate the absolute troposphere parameter, also known as precise point positioning (PPP; Zumbergre et al., 1997).Secondly, GNSS observations can also be processed by using a double-differencing approach, which requires at least two receivers and only the relative troposphere parameter can be obtained for a regional network if additional receivers are not considered (Duan et al., 1996).
In our study, the latter method is applied and the global positioning system (GPS) observations are processed using GAMIT/GLOBK (v.10.5; Herring et al., 2010) with the following configuration: the sampling rate is 30 s, and a cutoff elevation angle of 10 • is selected for satellite rays so as to ignore the influence of ray bending (Mendes, 1999).A global mapping function (GMF; Böhm et al., 2006) is adopted and the troposphere parameters are estimated; the intervals of estimated parameters for zenith troposphere delay (ZTD) and delay gradients are 0.5 and 2 h, respectively.To obtain absolute troposphere estimation, three International GNSS Service (IGS) stations (SHAO, BJFS, and SHAO) with a baseline longer than 500 km (Rocken et al., 1995) are also considered during final data processing, and the zenith wet delay (ZWD) can therefore be readily separated from ZTD based on the formula where ZHD is zenith hydrostatic delay, which can be accurately calculated using the Saastamoinen model (Saastamoinen, 1972) based on the measured surface pressure from the closest meteorological station.The SWD can thus be acquired as (Chen and Liu, 2014) where mf and gf refer to the wet mapping function (wet GMF) and gradient mapping function, respectively, ele and azi are the satellite cutoff elevation angle and azimuth angle, respectively, G w NS and G w WE are the wet horizontal gradients in the north-south and east-west directions, and R e refers to the zero difference (ZD) post-fit residuals, which is converted from the double difference (DD) residuals without an elimination of systematic effects using the approach proposed by Alber et al. (2000).
Therefore, the final SIWV can be obtained based on the following formula (Bevis et al., 1992): where ρ is the density of the water vapor (unit: g m −3 ), R v is the specific gas constant for water vapor (461.495J kg −1 K −1 ), k 3 = 3.776 × 10 5 K 2 hPa −1 , k 2 = 16.48K hPa −1 , and T m is the weighted mean tropospheric temperature, which is obtained from the empirical formula proposed by Chen et al. (2007) using measured surface temperatures.

GNSS troposphere tomography modeling
Following the principle mentioned in previous studies, the tomography area is discretized into a number of voxels, in which the water vapor density is assumed to be the same during a given period.Subsequently, the tomographic observation equation of a single signal ray can be established: where m, n, p represent the number of tomographic voxels in the latitudinal, longitudinal, and vertical directions, d i,j,k is the distance in the voxel (i, j, k) as shown in Fig. 1, which can be calculated based on the two intersections of two sides of the voxel and the signal ray, and x i,j,k is the water vapor density in voxel (i, j, k).Subject to the constellation of GNSS satellites and the specified distribution of GNSS receivers in a regional network, it is possible that many voxels are not crossed by any rays during the given period.Therefore, the observation equation established above is a large sparse matrix that will lead to a rank deficiency upon inversion of the equation.To overcome this problem, methods of imposing constraints or using the singular value decomposition (SVD) can be adopted.In this paper, the former method is used.A Gauss-weighted functional method (Song, 2004) is applied for the horizontal inner voxels, while the negative exponential function (Flores et al., 2000) is established to describe the relationship between vertical voxels.In addition, the initial value of water vapor density derived from radiosonde data over the first 3 days in the tomography area is also considered as a prior constraint on the location of the radiosonde station and the voxels vertically above it.Therefore, four kinds of input information are used in our tomographic modeling: where y and ρ rs are the column vectors of a series of SIWV observations and the initial water vapor density value derived from radiosonde data, respectively.A SIWV , A H , A V , and A P are the coefficient matrices for the observation equation and the horizontal, vertical, and prior constraints, respectively; m1, m2, m3, and m4 are the numbers of different types of input information.
the column vector of the water vapor density value to be estimated; n1 is the total number of voxels in the area of interest, which is equal to m × n × p; and ε SIWV , ε H , ε V , and ε P are the noise in the observation equation and the horizontal, vertical, and prior constraints, respectively.The covariance matrices for tomographic modeling can be expressed as follows: , where D 0 q (q = SIWV, H, V , P ), σ 2 0 q (q = SIWV, H, V , P ), and P 0 q (q = SIWV, H, V , P ) refer to the variance matrices, the unit weight variances, and the weight matrices for the four types of equations, respectively.3 Troposphere tomographic method considering the weights on input information As mentioned above, it is necessary to apply the appropriate weightings to guarantee the reliability and accuracy of the estimated parameters.There are two kinds of weightings to be considered for troposphere tomographic modeling.On the one hand, the weightings for the same type of input information should be determined, which will affect the accuracy of tomographic modeling if the improper weightings are applied.On the other hand, the weightings for different types of equations are also important in the final tomographic result.In the present study, a method considering both of the weightings mentioned above is proposed.

Determination of weightings for the same type of input information
For various types of equations, the methods for determining the weightings among the same type of input are different.
For the tomographic observation equation, the SIWV value is closely related to satellite elevation angle, and therefore the weighted function of the cutoff elevation angle should be considered: In addition, the SIWV data used when constructing the tomographic observation equation are derived from multi-epoch measurements during a given period, and thus the time span between the observed epoch and the tomographic epoch should also be considered.Here, the weighted function considering different time spans is expressed as where Tcorr = |Obs time −Tom time | , Obs time refers to the epoch of the measured signal ray, and Tom time is the tomographic epoch.T interval is the given period in which the water vapor density is assumed to be constant.In our study, 30 min was selected as a suitable value for T interval .Hence, the final weighting for each signal ray used was P i = P ele(i) • P Tcorr(i) , and the final weight matrix of observation P 0 SIWV can be expressed as For the horizontal weighted matrix P 0 H , a unit weight matrix is selected as For the vertical weighted matrix P 0 v and a priori weighted matrix P 0 p , radiosonde data for the first 3 days of the tomographic epoch were used to obtain the water vapor density for a specified height, and then the SD of water vapor density at different heights can be calculated.Therefore, the weighting for each voxel can be expressed as where SD k represents the SD calculated based on radiosonde data.The final weighted matrices of P 0 V and P 0 P are deter-mined as

Determination of weightings for the different types of equations
Usually, constraint equations are selected empirically, and the accurate unit weight variances of different constraint conditions are unavailable for most cases.Consequently, the appropriate unit weight variances should be given before the tomographic modeling resolution.Here, a method is proposed with which to determine the weightings for different types of equations.The flowchart of the proposed method is shown in Fig. 2 and the specified steps are as follows.
1. Initialize the unit weight variances for various types of equations as 1, and set initial weight matrix to P = D −1 .
2. Calculate the a posteriori residuals for all types of input equations v SIWV , v H , v V , and v P .In our study, singular value decomposition (SVD; Ran and Ge, 1997) is used to obtain the estimated water vapor density values X, as many elements in the design matrix are zero.Hence, the posterior residuals v SIWV , v H , v V , and v P can be calculated as where L = y m1×1 0 m2×1 0 m3×1 ρ rs m4×1 T .
3. Update the unit weight variances of all types of equations σ 2 0 SIWV , σ 2 0 H , σ 2 0 V , and σ 2 0 P .Here, simplified covariance component estimation (Mikhail and Ackerman, 1976) was used and the posterior unit weight variance σ 2 0 q (q = SIWV, H, V , P ) was as follows: where N = A T PA, N q = B T q P q B q , and B = [A SIWV A H A V A P ] T , n q (q = SIWV, H, V , P ) are the numbers of different types of equations, and tr is the rank of the related matrix.
In addition, the updated mean square error of the unit weight is also exploited to remove outliers from the tomography observation equation if v SIWV i > λ σ0 SIWV , λ is an empirical value, and λ = 3 is selected here after preliminary testing.
4. The co-integration test (Engle and Granger, 1987) is introduced to judge whether or not the estimated posterior unit weight variances are acceptable.This test is based primarily on whether or not the linear combinations of those unit weight variances are in a stable equilibrium relationship.The corresponding co-integration test procedure is as follows.
(a) Establish the relationship for the estimated unit weight variances.Here, the first-order autoregression variable sequence is selected for those variances, and the relationship can be established as where the values of ϕ are set to 1 for the firstorder auto-regression variable sequence used in this work.
(b) Calculate the test statistic t (ϕ) based on the Dickey-Fuller test (Dickey andFuller, 1979, 1981): where φ is the estimated value of ϕ by using a leastsquares method.
where N is the number of types of input equations.
(c) Give the null hypothesis and the alternative; due to the accurate weightings used for all types of equations being unavailable on the first time of use, the hypothesis is determined as follows.
Original hypothesis H 0 : the variable sequence is unstable, which means that where c is an arbitrary constant, and c = σ 2 0 SIWV is selected in our test; ite is the number of iterations.
6.The final unit weight variances for all types of equations are determined if the hypothesis H 1 is accepted.
4 Tomographic experiment and analysis

Tomography strategy
As mentioned in Sect.2.2, a tomographic experiment was carried out based on the data from 12 SatRef stations (Fig. 3) over 21 days (DoY 84-104, 2014).Different weather conditions (sunny days and rainy days) are included during the selected period.The time span for each step in the tomographic model solution is 0.5 h.The tomography area spans 113.87 • to 114.45 • longitude and 22.19 • to 22.54 • latitude (Fig. 1).The tomography vertical boundary is selected as 10.4 km a.s.l., assuming no water vapor above this altitude.The horizontal resolutions are 0.06 • (about 6.7 km) and 0.05 • (about 5.6 km) in the longitudinal and latitudinal directions, respectively, while the vertical resolution is 0.8 km.Therefore, the total number of voxels is 7 × 8 × 13: 8 longitudinal, 7 latitudinal, and 13 vertical.
To analyze the method proposed above, four schemes are designed.Only one step of the tomography experiment is presented at 00:00-30:00 UTC for DoY 87 in 2014 (Scheme 1); the whole day tomography experiment is presented for DoY 87 in 2014 (Scheme 2), the whole day tomography experiment is presented for DoY 89 in 2014 (Scheme 3), and a 21-day analysis of the tomography experiment is presented from DoY 84-104 in 2014 (Scheme 4).Two of the days mentioned above (DoY 87 and 89) are selected as they correspond to the two different weather conditions assessed; DoY 87 was a sunny day and DoY 89 was a rainy day with a total precipitation of 115.6 mm according to the records of the Hong Kong Observatory.

Analysis of one-step tomography modeling resolution
The a posteriori unit weight variances for all types of equations are first analyzed for Scheme 1. Table 1 lists the numerical results of troposphere tomographic modeling resolved by using the proposed method, which includes the number of iterations, the posterior unit weight variances for four kinds of input information, and the statistics from the cointegration test.It can be seen from Table 1 that the corresponding statistic is 1.777 after the first iteration, which means that the unit weight variances for four types of inputs do not show a stable equilibrium relationship.Therefore, the weightings of horizontal, vertical, and a priori equations are tuned, while c = σ 2 0 SIWV is selected for each iteration mentioned above.After five iterations, the unit weight variances for different inputs are 7. 286, 7.283, 7.273, and 7.256 mm 2 .The test statistic of 0.064 is also less than the given threshold of 0.1, and therefore the alternative hypothesis H 1 (Sect.3) should be accepted after five iterations and the estimated unit weight variances are assumed to be equal.It can be clearly noted that large discrepancies between the different unit weight variances were present after the first iteration, and the discrepancies decreased with an increasing number of iterations.

Validation of the proposed method under different weather conditions
Tomography experiments for Schemes 2 and 3 are performed to evaluate the proposed method on a sunny day (DoY 87) and a rainy day (DoY 89).According to the statistical results from the two days, the maximum number of iterations is 11 and the minimum is 4 when the alternative hypothesis H 1 is accepted.In addition, the residuals of the SIWV change with the cutoff elevation angle and the SD at each hour.Schemes 2 and 3 are also analyzed under different weather conditions (Figs. 5 and 6).Here, the residuals of the SIWV refer to the differences between the calculated SIWV values using the reconstructed water vapor density and the SIWV values obtained based on Sect.2.1.Generally, the  residual decreased as the cutoff elevation angle increased.
From Fig. 4 we can see that, compared to the residuals of SIWV on a rainy day, the residuals of SIWV are relatively small on a sunny day. Figure 5 also shows that the standard deviation (SD) of Scheme 2 is less than that of Scheme 3; however, the range of residuals is [−0.6, 0.6 mm] and the SD is less than 0.2 mm even on a rainy day, which means that the tomographic modeling has a high internal precision based on the proposed method under different weather conditions.

A comparison with independent techniques
Some studies have shown that accurate vertical water vapor information can be derived from radiosonde data (Niell, 2001;Adeyemi and Joerg, 2012;Liu et al., 2013); therefore, radiosonde data are selected as a reference to validate the water vapor density derived from other methods.Fortunately, there is a radiosonde station in the research area (red circle, Fig. 3).In addition, ECMWF also provides meteorological data that can be used to calculate water vapor density (Böhm et al., 2015;Wang et al., 2016).In our study, the spatial resolution of ECMWF data (ERA-Interim reanalysis data) for our study is 0.125 • × 0.125 • , and a total of 12 grid points were located in the tomography region.The geographic location  of the grid points is given in Table 2. Here, the water vapor information derived from troposphere tomography modeling under different weather conditions is compared with that calculated using radiosonde and ECMWF data.Figures 6  and 7 show the water vapor density change with height at 00:00-00:30 and 12:00-12:30 UTC for DoY 87 and 89, respectively.It can be seen that the tomographic water vapor density has good consistency with that calculated from radiosonde and ECMWF data for different heights.For a sunny day (Figs.6a and 7a), the water vapor density decreased rapidly with increasing height, while the rate of change was slower on a rainy day (Figs.6b and 7b); this was mainly because the water vapor content is large and the water vapor is more energetic on rainy days.Table 2 also compares the water vapor density derived under different weather conditions, which includes the bias, RMS error, and standard deviation (SD).It can be concluded from Table 2 that, compared to radiosonde results, the RMS error of the tomographic re-sults (0.888 and 1.621 g m −3 ) on a sunny day was not always better than that (1.738 and 1.359 g m −3 ) on a rainy day.For ECMWF data comparison, the RMS error of the tomographic results (1.133 and 1.009 g m −3 ) on a sunny day was even worse than that (0.938 and 0.720 g m −3 ) on a rainy day.This may have been because, on the one hand, the radiosonde and ECMWF data also have certain errors; on the other hand, the water vapor distribution is more complicated on a sunny day (Zhang, 2016).

Validation of the proposed method for many days
The integrated water vapor (IWV) profiles derived from tomographic modeling, radiosonde, and ECMWF data were first compared for the period DoY 84 to 104 in 2014.The IWV value for this radiosonde station was calculated by vertical integration (Brenot et al., 2006) using the water vapor information derived from the tomographic result at a specified epoch: 00:00-00:30 and 12:00-12:30 UTC, at which the  sounding balloon is launched daily.Figure 8 shows a direct IWV comparison between tomography and radiosonde data.For ECMWF data, only the data at the epoch of 00:00 and 12:00 UTC are selected to unify the temporal resolution with radiosonde data.Consequently, there are two epochs compared on each day, and the total number of compared epochs is 42.The IWV derived from the tomographic modeling using the proposed method provided good agreement with that from the radiosonde and ECMWF.Numerical results over 21 days (Table 3) were analyzed: the bias, RMS, and SD are −2.9, 5.8, and 5.0 mm and −2.1, 6.7, and 6.4 mm, respec-tively, when compared to those derived from radiosonde and ECMWF data.
Although the experimental comparison above indicates that the IWV time series derived from troposphere tomography conforms to that from radiosonde data, more tests are needed to validate the proposed method (as such, the compared result may not guarantee that the three-dimensional water vapor distribution is correctly modeled; Chen and Liu, 2014).Here, the radiosonde and ECMWF data are both considered to validate the water vapor distribution at different altitudes.For ECMWF data, the water vapor information at  different heights for the location of this radiosonde station is interpolated using the calculated water vapor density across the 12 grid points.The relative error is defined as follows to evaluate the relationship between the tomographic result error and height: where re refers to the relative error, x tomo is the water vapor density derived from tomographic modeling at different altitudes, and x 0 is the reference water vapor value calculated using radiosonde or ECMWF data.
Figure 9 shows the water vapor density profiles and the RMS and relative error change with height between tomography, radiosonde, and ECMWF data.Figure 9 shows that the mean water vapor density decreases with altitude, and the same trend can also be seen with regard to the RMS error.In contrast, the relative error increases with height (especially above 2 km) and the maximum relative errors reached approximately 200 and 180 % in the upper layers for radiosonde tomography and ECMWF tomography, respectively.The main reason is that the water vapor content of the upper layers is close to zero; therefore, a large relative error could be obtained even for a small discrepancy between radiosonde, ECMWF, and tomographic water vapor fields.In addition, Table 4 also lists the statistical results from the period of comparison.

Comparison of tomographic resolution with and without VCA
For the comparisons above, the proposed method with variance component analysis (VCA) for GNSS tomography has been validated and the tomographic result is good.To further validate the performance of VCA proposed in this paper, another comparison is performed between the results of a GNSS tomography solution based on the VCA and one without VCA methods for the period of 21 days.The two schemes are defined as the VCA scheme and the N-VCA scheme (without considering the weightings of different input information).The reconstructed water vapor density of voxels over the radiosonde stations using the two schemes are compared with that from the radiosonde data for the selected period.Figure 10 presents the RMS comparison between the VCA and N-VCA schemes during the selected 21 days: for most cases, the RMS error derived from the VCA method proposed in this paper is smaller than that from the N-VCA method because the VCA can be applied to find a realistic scenario for different observation types if the ac-  curate information is unavailable; however, it should be handled with caution since it tends to increase data scatter in the tomographic solution, which probably leads to some RMS errors in the VCA scheme being greater than those of the N-VCA scheme in our experiment (Möller, 2017).According to the statistical results, the bias, RMS, and SD of water vapor density differences between the VCA scheme and radiosonde data are 1.28, 1.33, and 0.38 g m −3 , while the values between the N-VCA scheme and radiosonde are 1.55, 1.59, and 0.36 g m −3 , respectively.

Conclusions
For GNSS troposphere tomography, some external constraints are required for the reasons mentioned in Sect. 2. Therefore, reasonable weightings among the observations for the same type of input information and weightings for different types of equations are a prerequisite for obtaining an accurate tomographic result.Consequently, this paper presents a method with which to determine the two types of weightings mentioned above.For the determination of the first type of weightings, some a priori knowledge was used, e.g., the functional relationship between the cutoff elevation angle and time span are implied for the observation equation; the SD of the different vertical heights derived from radiosonde data was introduced to determine the weightings for vertical and a priori constraints.For the weightings for different types of equations in the tomographic modeling, the method ensures that the posterior unit weight variances are in a stable equilibrium based on an iterative process.
A tomographic experiment was carried out to validate the proposed method using data from 12 stations in the SatRef network and data from radiosonde and ECMWF.Troposphere parameters were estimated using a double-differenced method, and the input SIWV for tomography experiments was calculated as described in Sect. 2. Different weather conditions (a sunny day and a rainy day) were both included during tomographic testing and the numerical results indicated that the proposed method could tune the unit weight variances of various types of inputs.A comparison of water vapor densities derived from tomography, radiosonde, and ECMWF showed that the applicability of the proposed method and the bias, RMS, and SD were 0.090, 1.401, and 1.398 g m −3 and 0.428, 0.950, and 0.848 g m −3 , respectively, under sunny and rainy conditions.A long-period comparison was performed using radiosonde and ECMWF data, and the Q.Zhao et al.: A troposphere tomography method bias, RMS, and SD of IWV differences were −2.9, 5.8, and 5.0 mm and −2.1, 6.7, and 6.4 mm, respectively.In addition, water vapor density profiles over 21 days were also compared with those from radiosonde and ECMWF: the water vapor density profiles and their RMS errors decreased with increasing height, while the relative error increased with height.Statistical analysis showed that the bias, RMS, and SD were 0.30, 0.91, and 0.86 g m −3 and 0.36, 1.63, and 1.59 g m −3 , respectively.
GNSS troposphere tomography considering two kinds of weightings is proposed before tomographic modeling resolution.In this study, although only the GPS-derived SIWVs were used to validate the effectiveness and applicability of the proposed method, it also can be used for other systems (Galileo, GLONASS, and BDS).In addition, more observations can also be used for water vapor tomography, such as data from a radio occultation apparatus, solar spectrometers, a water vapor radiometer (WVR), an interferometric synthetic aperture radar (InSAR), and the Constellation Observing System for Meteorology, Ionosphere, and Climate (COS-MIC) radio occultation (Liu et al., 2013;Xia et al., 2013;Alshawaf, 2013;Heublein et al., 2015;Benevides et al., 2015), which can also be regarded as a new type of input information in the proposed method.Hence, more voxels are expected to be crossed and the geometric structure of the observation equation will be enhanced.
Competing interests.The authors declare that they have no conflict of interest.
Figure 1.Geometric illustration of a signal ray crossing the discretized voxels.

Figure 2 .
Figure 2. Flowchart of the proposed method considering the weightings for different types of equations.

Figure 3 .
Figure 3. Geographic distribution of the Hong Kong Satellite Positioning Reference Station Network (SatRef) and radiosonde station.

Figure 4 .
Figure 4.The change in SIWV residuals with cutoff elevation angle for Schemes 2 and 3.

Figure 5 .
Figure 5.The comparison of SD for Schemes 2 and 3.

Figure 9 .
Figure 9. Water vapor density, RMS, and relative error change with height for Scheme 4.

Figure 10 .
Figure 10.A comparison of RMS for the VCA and N-VCA schemes during the period of 21 days.

Table 1 .
Posterior unit variances for different input equations and test statistics calculated by using the proposed method for Scheme 1.

Table 2 .
Statistical results of water vapor density between radiosonde, ECMWF, and tomography for Schemes 2 and 3 (unit: g m −3 ).Comparison of IWV time series derived from radiosonde, ECMWF, and tomography modeling for Scheme 4.

Table 4 .
Statistical results of water vapor density between radiosonde, ECMWF, and tomography for Scheme 4 (unit: g m −3 ).