Near real-time water vapor tomography using ground-based GPS and meteorological data : long-term experiment in Hong Kong

Water vapor tomography is a promising technique for reconstructing the 4-D moisture field, which is important to the weather forecasting and nowcasting as well as to the numerical weather prediction. A near real-time 4-D water vapor tomographic system is developed in this study. GPS slant water vapor (SWV) observations are derived by a sliding time window strategy using double-difference model and predicted orbits. Besides GPS SWV, surface water vapor measurements are also assimilated as real time observations into the tomographic system in order to improve the distribution of observations in the lowest layers of tomographic grid. A 1-year term experiment in Hong Kong was carried out. The feasibility of the GPS data processing strategy is demonstrated by the good agreement between the time series of GPS-derived Precipitable Water Vapor (PWV) and radio-sounding-derived PWV with a bias of 0.04 mm and a root-mean-square error (RMSE) of 1.75 mm. Using surface humidity observations in the tomographic system, the bias and RMSE between tomography and radiosonde data are decreased by half in the ground level, but such improved effects weaken gradually with the rise of altitude until becoming adverse above 4000 m. The overall bias is decreased from 0.17 to 0.13 g m−3 and RMSE is reduced from 1.43 to 1.28 g m −3. By taking the correlation coefficient and RMSE between tomography and radiosonde individual profile as the statistical measures, quality of individual profile is also improved as the success rate of tomographic solution is increased from 44.44 to 63.82 % while the failure rate is reduced from 55.56 to 36.18 %.


Introduction
In the last decades, GPS water vapor tomography has shown its great potential to retrieve the 3-D/4-D distribution of atmospheric water vapor (e.g., Bastin et al., 2005;Lutz, 2008;Champollion et al., 2009;Rohm and Bosy, 2009).Wet refractivity or humidity fields at high spatiotemporal resolution can be reconstructed by a tomographic technique with combining numerous tropospheric slant observations from a dense ground-based GPS network.Information provided by these fields is very important for improving the accuracy of modern numerical weather prediction models, especially for the quantitative precipitation forecasting (Emanuel et al., 1995;Bauer et al., 2011).Furthermore, monitoring water vapor by GPS has the advantages of low-cost and all-weather measurement in comparison to other means, e.g., the radio sounding.
The classic GPS tomography model, which was first realized by Flores et al. (2000), is defined by discretizing the atmosphere space using a spatial grid consisting of several voxels.Water vapor distribution inside each voxel is assumed to be homogeneous.The input data of the tomographic system, the GPS tropospheric slant wet delay (SWD), the or slant water vapor (SWV) observation, can then be represented by a summation of separated SWDs/SWVs distributed in each cell which the ray passes through.The estimated parameters are the wet refractivity or absolute humidity values in all of voxels, and the reference point of each parameter is approximated to be the central point of corresponding voxel.Perler et al. (2011) presented new voxel parameterized approaches using trilinear and bilinear/spline interpolations.Instead of one constant wet refractivity parameter in each voxel, estimated point wet refractivity at grid nodes are the unknown Published by Copernicus Publications on behalf of the European Geosciences Union.

P. Jiang et al.: Long-term experiment in Hong Kong
parameters in tomographic equation.Adverse effects of the classical discretization method on tomographic solution can be reduced by the parameterized approaches especially when the voxels are not small in size (Nilsson, 2005).Rohm (2012) studied the error sources and their propagation in tomography solution, and the topology of cofactor matrix shows the most significant influence on the final solution uncertainty.Meanwhile, research about the distributions of signal slant paths and their intersection points (Bender et al., 2009) demonstrated that the density of slants is maximum in the boundary layer while the intersection points distribute the most sparsely, which leads to poor topology of coefficient matrix for the estimated parameters in corresponding tomography level.Therefore it is necessary to introduce the surface water vapor measurements to overcome aforementioned disadvantage.The classic method is to use the synoptic observations to provide prior information about the water vapor field, e.g., modeling the prior 3-D water field with surface observations or constraining directly the estimated parameters at the lowest grid level to equal the measured point moisture value (Troller et al., 2006;Bender et al., 2011).However, these methods will produce additional errors because of the spatial distance between the voxel central points and the observation stations.Furthermore, in a classic tomographic processing batch, GPS slant observations observed in a time period are combined to form the equation, but at the same time, only one prior field is used.The time variation information of water vapor contained in the surface observations are lost.
In this study, we develop a near real-time (NRT) water vapor tomographic system by applying a bilinear/spline parameterized method.In the proposed approach, not only GPS slant observations but also surface meteorological measurements are used as real time observations; A Kalman filter is applied to reconstruct the 4-D humidity fields.An experiment over a whole year of 2010 in Hong Kong was carried out to investigate the performance of our tomographic algorithm.Due to its location on the edge of the South China Sea and in the path of the wet Asian monsoon, Hong Kong experiences more rainstorms than most other cities (Liu et al., 2013).The influence of the surface humidity on the tomographic solution is also studied.
In the following sections, we present a near real-time water vapor tomography approach.In Sect.2, the tomographic algorithm is presented.Section 3 describes our data processing strategy in detail.Tomographic results are compared with other data in Sect. 4. Conclusions are given is Sect. 5.

Tomographic method
The SWV, one of the important products from GPS remote sensing atmosphere, is physically the integrated water vapor amount along the curved ray path from GPS satellite to ground-based antenna.It can be expressed by the formula (Miimla et al., 2007) where ρ represents the density of liquid water, H s denotes the absolute humidity (we will refer the absolute humidity as the humidity in the following paper) and s is the slant path.
Strictly the signal path is bent when passing through the atmosphere.However, the contribution of this bending term to Slant Total Delay (STD) is less than 1 cm, which means about 1.5 mm to SWV for elevations greater than 15 • (Bevis et al., 1992) therefore it is neglected for simplicity.The slant path then can be regarded as a line of sight from GPS antenna to satellite.
Intensive observations of regional atmosphere are available with a large number of SWVs obtained from a ground based GPS network.In order to use this type of integral observation for reconstructing the humidity field, space above the network is discretized into voxels both in horizontal and vertical directions.We made the hypothesis that all atmospheric water vapor is concentrated in the altitude lower than the top height of grid model (10 km in our study).Only SWVs that penetrate into the grid from the top boundary are used in tomography.Each SWV is then partitioned into several segments in the volumes.Individual SWV part in each voxel expressed as Eq. ( 1) can be solved approximately by Newton-Cotes quadrature rules (Forsythe et al., 1977).The closed Newton-Cotes formula of degree 4 requires known water vapor density at the 5 equally spaced points.It is a simple way that interpolating the estimated humidity parameters of grid nodes to the interesting points.Bilinear and spline interpolation formulas, respectively, in the horizontal and vertical directions are adopted in our solution with the bilinear one performed prior to the other.Positions of all points have to be expressed in ellipsoidal coordinates to take the Earth's curvature into account.Then adding the coefficient vectors of every discrete SWV part, the sum is the final coefficients vector for a whole individual SWV.Combining all SWVs' coefficient vectors together, the data set of SWV are correlated to the humidity field according to the equation: where SW V is the SWV observations vector, x is the state vector of the humidity at all voxel vertexes, and A is the design matrix containing the coefficients obtained from the bilinear/spline parameterized approach.Besides SWVs, surface synoptic observations can also be assimilated into tomographic system.It will be helpful to improve the accuracy of humidity results near the lower atmospheric boundary layer especially when the SWV observations distributed inhomogeneously in the groundlevel voxels.Furthermore, specific surface observing network can provide even more surface observations.Sensors included in such system generally are barometric pressure, temperature and relative humidity ones.It means that the atmospheric state around the GPS receivers can be monitored accurately in high time resolution as well as GPS observations.With data provided by these systems, humidity at the considered location is calculated with the following formula as (Sheng et al., 2003) where H is humidity in g m −3 , R v is the specific gas constant for water vapor, T is temperature in K and e denotes vapor pressure in hPa, which is related to temperature and relative humidity.Bilinear/spline interpolation of voxel vertexes' parameters as described previously can also be performed to construct the observed equation of surface measurements: where H sfc is the surface humidity measurements vector and B is the coefficient matrix.Due to the limited distribution of SWVs and surface observations, in the observation matrix there are usually some parameters of which coefficients are always zero or very small, resulting in the ill-posed equations combined by real observations.The humidity fields retrieved from such solution are unstable because of the probable existence of unreasonable sharp change.To solve the issue, horizontal and vertical inter voxel constraints are added as pseudo observations to our tomographic system.These constraints are based on the assumptions that the humidity at a voxel corner is determined by the horizontally and vertically neighbors in different forms.Horizontally Gaussian smoothing with a controllable width of one voxel's side length and vertically exponentially decreasing assumptions with statistical scale height of 2.53 km (Song, 2004) is implemented.The constraints are introduced as the formula Combining all the observations obtained above, the final form of tomographic observational equation of each epoch is Finally, assuming the time evolution of humidity values is a random walk stochastic process, we apply a standard Kalman filtering (Herring et al., 1990;Gradinarsky and Jarlemark, 2004) to identify the time variation of the studied humidity field.The Kalman filtering method has the advantage of obtaining a 4-D humidity field.In detail, it updates the reconstructed field once new data are generated, so that the time variation information observed by real observations can be represented better in the final tomography solution.In the filter predict step, the state transition matrix is set to be a unit matrix.The process noise matrix is designed as a diagonal matrix of which the diagonal elements are obtained by statistics of surface measurements with assuming an exponential law in the vertical humidity distribution.The details will be described in Sect.3. In the filter update step, uncertainty of SWV in the zenith direction is considered to be 1.7 mm (Deng et al., 2011).It is reported that accuracies of surface pressure and temperature are better than ±0.08 hPa and ±0.2 • C, respectively, and surface relative humidity performance is better than ±2 % at 25 • C (More information in http://www.paroscientific.com/pdf/met4DataSheet.pdf).By taking all the measurement errors into account, a noise level lower than 3 % should be appropriate for surface humidity observation.Weights of both constraints are designed as small as 10 −2 to reduce the smoothing effect of constraints on the variability of humidity fields in real case.

Processing strategy
12 continuously operating reference stations are mounted in Hong Kong Satellite Positioning Reference Station Network (SatRef).The distance between two neighboring reference stations is roughly 10 to 15 km.Horizontal distribution of these stations is approximately uniform (see Fig. 1).All stations are equipped with Leica GPS receivers and Paroscientific Met3/Met3A/Met4A meteorological sensors.Chokering antennas (Leica AT504 + LEIS) are used to minimize the influence caused by the multipath errors in observations.GPS and surface meteorological data in 5 and 60 s intervals are collected in real time.These data can be freely downloaded afterwards from website (http://www.geodetic.gov.hk/smo/gsi/programs/en/GSS/satref/satref.htm).In our project, a complete data set of the year 2010 was obtained.
For GPS data, to simulate near real-time tomographic experiment, a sliding time window strategy (Foster et al., 2005) was applied with a 6-hour-wide time window and was stepped forward 1 h at a time.We employed GAMIT software (Herring et al., 2010) to process the GPS data.The atmospheric parameters were estimated by a piecewise linear (PWL) interpolation model, together with the site coordinates and orbit parameters in each time window.Observations with elevation lower than 10  (Kačmařík and Skřivánková, 2011).Interval of the atmospheric ZTD solution was 1 h when the gradient parameters were 2 h.However, in the case of GAMIT software processing data with a double-difference model, it is necessary to introduce three or more distant reference GPS stations to form long baselines for estimating absolute ZTD (Duan et al., 1996).IGS sites include the SHAO site in Shanghai, the KUNM site in Kunming, and the TWTF site in Taiwan were added to GPS network.
Using double-difference network solution mode can eliminate the common satellite and receiver clock errors contained in the phase observations and greatly reduce the impact of orbit errors, especially for a small dense network.Consequently the observations are more sensitive to the small but important variations in atmospheric delays.On the other hand, the precise point positioning (PPP) solution needs to be given accurate satellite clock and orbit estimates before processing.In consideration of the unavailability of realtime high-precision orbit and clock products from IGS at the present, the double difference solution is favored for NRT practical application in spite of the increasing computational costs.
The SWV, main input data of the tomographic equation, is not estimated directly in GPS data processing.It can be computed by formula (Braun et al., 2002): where PWV is the perceptible water vapor in zenith direction; M w is the wet mapping function, we used global mapping function (Boehm et al., 2006) as GAMIT software did; is the scaling factor for mapping wet delay to water vapor.It is dependent on the atmospheric weighted mean temperature T m ; ε and φ denote the slant path elevation and azimuth; the horizontal gradient L g and the post-fit phase residual R represent the anisotropic component.L g formula varies with software.In GAMIT software, it is expressed as where C is a constant equal to 0.003, G N and G E are the north and east gradient parameters.
PWV is transformed from zenith wet delay (ZWD) also by the factor.ZWD can be separated using formula ZWD = ZTD − ZHD, since the zenith hydrostatic delay ZHD (Zenith Hydrostatic Delay) can be calculated very precisely using surface pressure and the Saastamoinen model (Saastamoinen, 1972).The model is a function of latitude θ in radian and height above the geoid H in km: where P 0 represents surface pressure in hPa, and The horizontal gradient estimates also contain hydrostatic and wet parts.The hydrostatic part is determined primarily by the pressure and temperature horizontal gradients at the GPS station.Considering that the experimental region is a small moist area with highly variable humidity, we are confident that the dry gradients account for a rather small proportion of total gradients which even may lower than the noise.The horizontal gradient estimates therefore are added to the mapping values of PWVs in the path directions without subtracting the dry parts.The post-fit residuals are included to compensate for the higher-order terms of asymmetric slant delay not modeled by gradient estimation scheme.Doubledifferenced residuals have to be transformed to line-of-sight residuals using the method proposed by Alber et al. (2000).
The relational expression between the conversion constant and T m is (Bevis et al., 1994) = where ρ is the density of liquid water, R v is the specific gas constant for water vapor, k 2 = (17 ± 10) K mbar −1 , k 3 = (3.776± 0.014) × 10 5 K 2 mbar −1 .T m can be estimated using the surface temperature measurement T s .Bevis et al. (1992) firstly studied the linear regression between T m and T s based on radiosonde profiles from 9000 sites over the United States.
Varying relationships were derived in different areas to better fit the local atmospheric features.Empirical linear formula of T m and T s in Hong Kong was given by Wang et al. (2011) as which was used in our study.The benefit of using a local specific formula is shown in Sect. 4. The SWVs were generated corresponding to the sampling interval of surface meteorological measurements.At the tomographic epoch, the ZWDs and horizontal gradients were obtained by linear interpolating ZWD and gradient estimates at the nearby knots of time window.It was quite necessary since the atmospheric delay model used by GAMIT software was PWL.The SWVs which intersect the discrete grid model at the lateral boundaries instead of the top should not be included in tomographic input observations, so a pre-judgment of the intersection was carried out before each SWV computation.
Humidity at the GPS sites was calculated using Eq.(3) with surface measurements.A simple data quality control method was performed to the observations.It is reasonable to regard the atmospheric state around the sensor to be steady during a very short period of 5 minutes.Before every epoch, the mean values of previous five humidity observations were stored at all stations.Measurements that differed from the corresponding averages by greater than 5 % were marked as outliers and rejected in tomographic computation.In our study, a very limited quantity of deleted data was found at each station, so the stability of surface observations was proved.
The reconstruction area covered from longitude 113.797 to 114.432 • E and latitude 22.157 to 22.586 • N, and the height is from 0 to 10 000 m.The tomography grid is shown in Fig. 1.The horizontal discretization is about 10 km, not much better than the mean GPS interstation distance.16 layers are distributed in vertical direction, and the thickness of these layers increase with height which is max to 1000 m in the high layers.

Characteristic model of temporal variations in atmospheric humidity
In the Kalman filter predicted step of tomographic processing, as described in Sect.2, the diagonal elements of process noise matrix need to be determined.We applied structure functions D to characterize the temporal variability in humidity (Jarlemark and Elgered, 1998).D is defined as where t is the time epoch of the measurement, τ is the time lag, and the angle brackets indicate the expected value which was replaced by the average value in our calculations.Time series of humidity measurements at the three lowest sites in the year 2010 were analyzed.Although the measurements were available at intervals of 60 s, structure functions at small time lags are affected by the short term instrumental white noise badly.This effect can be ignored at larger time lags.So we selected the time lags from 0.5 to 4 h every 0.5 h. Figure 2 shows the calculated structure functions at different time lags.The differences between these structure functions are very small.The structure function generally follows a power law, where the c represents a constant and α is the power-law index, with values varying between 1 (Jarlemark and Elgered, 1998) and 2/3 (Treuhaft and Lanyi, 1987).Assuming a value for α of 1, we fit a constant c to the mean value of the three structure functions.A variance rate of 0.005 (g m −3 ) 2 min −1 is derived.It is a small value and varies little over our study region so we believe this value indicates the structure function of the lowest atmospheric layer.The upper layers' structure functions are derived based on the assumption of water vapor vertical exponential decrease distribution, so the structure function at altitude h can be computed using equation where h 0 is the altitude of the lowest layer and h sc is the wet humidity scale height of 2.53 km.Structure function at each tomographic voxel vertex can then be determined depending on its height, which can be used in the Kalman filter predicted step to form the diagonal of process noise matrix.

Result and discussion
In this section, the data processing strategy and tomographic algorithm are assessed.in the year 2010 were used.The data can be freely downloaded from the internet (http://weather.uwyo.edu/upperair/sounding.html).The sounding balloons were launched automatically every 12 h at UTC time 00:00 and 12:00 every day and the observation heights were up to higher than 15 km.However, radiosonde data missing are found at 00:00 of Day Of Year (doy) 95, 12:00 of doy 103 and 12:00 of doy 210.The available radiosonde data are obtained to make comparisons with the atmosphere products derived from our solutions.Except for the missing epochs, our GPS meteorological products derived from the data processing as described in Sect. 3 at all sounding balloon launched time are extracted for comparison purpose.

Representativity of surface data
The representativity of surface measurements should be taken into account because they may be greatly influenced by the surrounding of surface meteorological sensors.Therefore radiosonde data were interpolated (or extrapolated) to the altitude of each surface site and compared with the surface humidity at balloon launched time.The mean and standard deviation of the difference at each site were given in Table 1.Nearly zero bias and the smallest standard deviation were found at GPS site named HKSC which is the nearest one to the radio sounding site.The relative larger differences between the other stations and radiosonde are attributed to the spatial humidity variations between the sites, so the farthest GPS site named HKNP had the largest difference.We think that the reason for the good representativity is that all surface meteorological stations are located in very open area.Another probable factor that could impact the surface data is the weather situation.We selected the HKSC site to study such issues.The data observed at UTC time 00:00 and 12:00 were divided into two sets and compared with corresponding radiosonde data, respectively (see Table 2).The mean differences were both almost zero, and the standard deviation varies little between two data sets.Moreover, the data from opposite seasons, winter (month January, February and December) and summer (month June, July and August), were also investigated (see Table 2).The biases were both rather small.The larger standard derivation of difference in summer is expected because of the increasing water evaporation and more air flow in the boundary layer in this season.The accuracy, however, can still be accepted to improve our tomographic solution in consideration of the high humidity value in summer.

Comparison of PWV time series
In order to generate SWV, as Eq. ( 7) shows, PWV need to be mapped to the signal path by mapping function, so that the errors of PWV can be amplified in the SWV according to propagation of errors.It is important to confirm that the accuracies of PWVs derived by our approach can meet the accuracy requirement of meteorological application (Shoji, 2009).We compared all of the GPS-derived PWVs with Radiosonde (RS)-derived PWVs at radiosonde launched time.The GPS site named HKSC was selected for comparison.Its horizontal distance from the RS site was around 2.4 km so the effects of horizontal variations in atmospheric properties were neglected (see Table 1).GPS-derived PWV time series was generated by extracting the PWVs from the end time knot of each time window.Two time series of GPS-derived PWV in the year 2010 were generated using different T m computational equations: Bevis T m equation and local specific T m equation as Eq. ( 12).However, the altitude of GPS site is lower than that of RS site, which the former is about 20 m and the latter is about 66 m.There would be a PWV bias between them duo to the water vapor distributed in the altitude offset part.Therefore corrections to all RS-derived PWVs were implemented by regarding the surface measurement at HKSC site as the bottom measurement for the radiosonde profile.The three time series of PWV and their differences are shown in Fig. 3a.Table 3 shows the statistics of these data.The mathematical expressions of the statistics are given in Appendix A. Scatter plot (Fig. 3b) illustrates that the  GPS-derived PWV, see Fig. 3 and RSPWVs are highly consistent.The linear expressions between GPS-derived PWVs and RSPWV are given as Eq. ( 16).
It can be seen that there is a bias between the RSPWV and GPSPWV based on Bevis T m equation.Using a local specific T m equation can remove the bias largely from 0.41 to 0.04 mm.The RMAD and RMSE decrease slightly.In our following work, the Eq. ( 12) was selected for all GPS sites.The good agreement between GPSPWV and RSPWV also indicates the feasibility of the time window processing strategy that we used to retrieve the water vapor content and makes sure the accuracy of our tomographic solution.

Comparison of humidity results
Except for the missing data, a total of 741 radio sounding humidity profiles with nearly 30 000 point observations in the year 2010 were obtained.As shown in Fig. 1, the RS station doesn't locate on the intersection of our tomographic grid lines.Therefore we interpolated the four neighboring tomographic profiles to retrieve the profile over the RS site.Horizontally bilinear interpolation and vertically spline interpolation were performed.The interpolated tomographic profiles could be compared with RS profiles directly.

Quality of point humidity
With assimilation of surface humidity observations or not, two different sets of reconstructed humidity fields could be achieved as -TSSH: Tomographic Solution with Surface Humidity observations assimilated -TSnoSH: Tomographic Solution with NO Surface Humidity observations assimilated.
Comparing both sets of tomographic reconstruction with radiosonde data could investigate the impact of such assimilation to the tomographic results.Histogram of the differences between both tomographic solutions and measured values retrieved from radiosonde data is shown in Fig. 4. It indicates that the mean value and standard deviation of TSSH-RS differences are both smaller than those of TSnoSH-RS differences.However, the vertical structure of atmospheric humidity is highly related to the altitude and observation time.For example, a higher tomographic humidity result at low altitude has better quality than the lower one at high altitude even if their differences with corresponding RS observations were equal.So quality of tomographic reconstruction cannot be fully assessed by such statistics of the whole data set without regarding the relevant factors.Therefore, we classified all of the nearly 30 000 points' humidity results into 16 layers which were defined in the tomography grid.In each layer, mean humidity of radio sounding was computed, statistics of the differences between tomography and radio sounding were performed.The results are shown in Fig. 5.It illustrates that the mean humidity shows an exponential decreasing gradient.For both tomographic solutions, there are negative biases in the layers that are lower than 3500 m except for the bottom layer, while positive biases in the layers that higher than 3500 m.The RMSE of TSSH is much smaller than, nearly a half of, the RMSE of TSnoSH in the lowest layer, and the difference decreases with altitude to nearly zero above 2000 m.In the upper air higher than 4000 m, the RMSEs of TSSH are slightly larger than TSnoSH.The relative RMSE, ratio between RMSE and corresponding mean humidity, reduces from 0.16 for TSnoSH to 0.08 for TSSH in the lowest level.However, relative RMSEs are small in the low levels while  large in high altitudes.It indicates that the humidity in high atmosphere is difficult to retrieve precisely because of the small humidity value and the limited accuracy of tomography solution (Nilsson et al., 2004).For example, in the two layers between 3000 and 4000 m, although the RMSE of TSSH decreased from 1.41 to 1.35 g m −3 , the humidity also declined from 5.00 to 3.64 g m −3 leading to the increase of the relative RMSE from 0.28 to 0.37.
Furthermore, in order to reveal the effect of surface data inclusion more clearly, we divided the atmosphere into two layers: a lower layer (altitude < 4 km) and an upper one (altitude > 4 km).We repeated the comparisons between tomographic solutions and radiosonde in the same way as previously described.The results are shown in Table 4.It indicates that the accuracy of tomographic solution is better in the lower layer than that of the upper layer.TSSH has smaller RMSE than TSnoSH at lower level while slightly larger in upper level.The results demonstrate that the assimilation of surface humidity observations can improve the quality of tomographic solution near the surface of network.These positive impacts, however, are bound in the lower atmospheric layer and weaken with the rise of height because all of the ground meteorological sites located in low altitude.The corrections to the only-GPS tomographic profile in low parts should be compensated and the compensations are "spread" in all of the higher layers, which leads to a slight increase of differences.Considering that the water vapor is concentrated in the lower atmosphere, which has great impact on weather change, we believe that this influence produced by adding surface observations to the tomographic solution is positive for the application of tomography in real time weather services.

Quality of entire humidity profile
The consistency of the individual humidity profile between tomography and RS data was also investigated.For this purpose, Pearson's product moment correlation coefficient (PCC, see Appendix B) and RMSE of each individual profile were selected as the statistical measures.Larger PCC value indicates better consistency between two profiles.Figure 6 shows the time series of the PCCs and RMSEs for all obtained humidity profiles.Visually most of PCCs of TSnoSH are lower, while RMSEs are larger, than those of TSSH.For each profile, the PCCs of both solutions with their difference below 0.001 were considered equal as well as the RMSEs with difference smaller than 0.01 g m −3 .Using this criterion, 56.23 % of TSSH's PCCs and 72.90 % of TSSH's RMSEs, which take majority of corresponding statistic's amount, are better than TSnoSH's (see Table 5).Only 21.41 % TSSH's PCCs are smaller than TSnoSH's and even no TSSH's RMSE is larger.
It is difficult to define general criteria which can be used to automatically check the success and failure of a tomographic There are several problems in evaluating the tomographic solution by using radiosonde data.The measure principles of tomography and radio sounding are significantly different.Radiosonde measures the atmosphere close to the sensors point by point when the sounding balloon ascends, while our tomographic solutions can merely provide point humidity of the grid nodes, which only represent the mean state of each atmospheric layer.Radiosonde data can be influenced by complicated factors which are difficult to estimate using statistics.For example, when radiosonde penetrates the clouds, the observations maybe suffer from perturbations (Wang and Rossow, 1995;Chernykh et al., 2001).However, these perturbation factors usually do not occur in the normal atmosphere.Furthermore, although the accuracy of radiosonde data is good overall, there are still unexpected errors in RS profiles and sometimes the errors are rather big (Wang et al., 2002;Ciesielski et al., 2010).In that case, it is not appropriate to compare the tomographic solution with radio sounding.However, empirical criteria were still established in our study to provide an overview of the success rate of tomographic solutions.As shown in Table 6, only a tomographic profile with both a high PCC (> 0.97) and a small RMSE (< 1.3 g m 3 ) was regarded as a success.The remaining profiles were categorized into failure set.Apparent increase in the total amount of success was gained by using surface humidity measurements.All of the comparisons listed above obviously demonstrate that surface humidity measurements inclusion can also improve the quality of entire tomographic profile.
Some profiles successfully retrieved from TSSH are shown in Fig. 7.These profiles have different characteristics.Panel a and b show profiles with regular exponential decrease in altitude, but the profile in panel b is dryer than panel a and the decline rate is slower.The profile a also shows smaller vertical gradients in the high levels than that at low altitudes, which is similar to the radiosonde profile.The Profile in panel c, differently, shows a nearly linear vertical gradient.In panel d, there is an irregular increase of humidity below 1000 m shown by radiosonde data, and tomographic solutions also present the similar features.It demonstrates that the inversions in low altitude can still be reconstructed.Panel e shows a profile in which the humidity decrease very rapidly from 19 g m −3 to nearly zero between 0 and 3000 m.The resolvability of different humidity profiles using tomography technique is demonstrated by our NRT experiments.
This study also indicates the inability of our tomography approach to reconstruct the spike layers in upper levels.As shown in Fig. 7f, there is an apparent humidity inversion structure between 2200 and 5000 m represented by radiosonde profile.However, the tomographic solution failed to retrieve it while it was "spread" into the neighbor layers.The humidity inversion, which means water vapor concentration at higher altitudes, appeared more frequently in the colder dryer weather during our experiment period.This fact induced more failures of tomographic solution in colder seasons as illustrated by Fig. 6.Studies from different researchers (e.g., Flores et al., 2001b;Champollion et al., 2005;Nilsson and Gradinarsky, 2006;Perler et al., 2011) also revealed the same disadvantage.These studies indicated that the geometry of ground-based GPS network, especially in vertical direction, has great influence on the tomographic ability to retrieve the spike layers above the network.Unfortunately, the Hong Kong network is a relative flat network that has deficient geometry in vertical direction.Another reason for this incompetency may be that we added vertical constraints as pseudo observations, as Eq. ( 5) expresses, to the tomographic equation.However, the weights of the constraints are much lower than the real observations, so the smoothing effects should not be great.Our studies about the tomographic solutions without vertical constraints showed that in several cases the spikes could be obtained but the reconstructed magnitudes of perturbations were always weaker than the radiosonde.In most cases, the accuracy of tomographic solution declined without humidity inversions being resolved.

Conclusion and outlook
A near real-time water vapor tomography system is developed in our study.Sliding time window strategy and double-difference network solution are implemented to retrieve the GPS water vapor observations.A tomographic bilinear/spline parameterized algorithm is applied.With the specific advantage of such algorithm, not only GPS SWVs but also surface humidity measurements with high temporal resolution are used as observations to form the tomographic equation.
A long-term tomography experiment in Hong Kong over the whole year of 2010 is carried out.4-D humidity structure is reconstructed by using Kalman filter.Good agreement between tomographic solution and radiosonde data demonstrates the superior performance of our tomography algorithm.Compared results using different statistical methods all indicate that the significant improvement in tomographic solutions can be gained by introducing surface humidity observations.The distribution of errors is improved by reductions in the mean value and deviation.Accuracies in low layers increase substantially, especially in the bottom layer the RMSE between tomographic results and radiosonde data is reduced by near 50 %.However, the positive impacts reduced quickly with altitude increasing.An assessment of individual profile quality shows that the percentage of good tomographic profiles rises from 44.44 to 63.82 % and the proportion of poor solutions drops from 55.56 to 36.18 %, which attributed to the inclusion of surface observations.We are confident that assimilation of surface observations is conducive to the practical application of tomography technique in real time atmospheric monitoring.
There are still some issues needed to solve.The tomography method we proposed in this paper is unable to reconstruct the spike layers in the upper air.One possible solution is introducing COSMIC data into tomography (Xia et al., 2013).COSMIC can provide the humidity profile with high accuracy in high altitudes.Besides this, automatic quality control without comparison to other independent data should be developed.Robustness of the real time tomography system can be greater by independent quality control methods.

Figure 2 .
Figure 2. Humidity structure functions for the three lowest Hong Kong SatRef stations and their mean value.Altitude of GPS site HKPC is 18 m, HKSC is 20 m and HKKT is 34 m.Average structure functions are calculated for the time lags from 30 min to 4 h using 1 year of data.

Figure 3 .
Figure 3.Comparison between PWV derived from radiosonde and GPS.RSPWV represents the PWV derived from radiosonde data.GPSPWV B and GPSPWV S respectively represent the GPS-derived PWV using Bevis T m equation and local specific T m equation.

Figure 4 .
Figure 4. Histogram of the differences between different tomographic solution and radiosonde data.Mean value of the differences is 0.1294 g m −3 for TSSH and 0.1659 g m −3 for TSnoSH, with standard deviation 1.2797 g m −3 for TSSH and 1.4295 g m −3 for TSnoSH.

Figure 5 .
Figure 5. Statistics of the humidity differences between tomographic solutions and RS data in 16 layers.Statistical measures are mean humidity (upper left), BIAS (upper right), RMSE (down left) and relative RMSE (down right).

Figure 6 .
Figure 6.Time series of PCCs (upper) and RMSEs (down) between different tomographic humidity profiles (TSSH blue line and TSnoSH red line) and the ones derived from radiosonde data.

P.
Figure 7. Representative cases of tomographic solution and radiosonde data during our long-term study.

Jiang et al.: Long-term experiment in Hong KongTable 1 .
The differences between surface humidity measurements and radiosonde interpolation (or extrapolation) data at all groundbased sites.

Table 2 .
Divided surface humidity measurements at HKSC site using different basis of classification and their comparison with radiosonde data.

Table 4 .
Mean humidity derived from radiosonde in lower and upper layer of atmosphere.The BIAS, RMSE and relative RMSE of tomographic solutions compared with radiosonde are also computed in two layers.

Table 5 .
Comparison of the PCC and RMSE between individual profiles derived from two tomographic solutions and radiosonde data.

Table 6 .
Three sets of profiles divided by the quality.