Assessing water vapor tomography in Hong Kong with improved vertical and horizontal constraints

In this study, we focused on the retrieval of atmospheric water vapor density by optimizing the tomography technique. First, we established a new atmospheric weighted average temperature model that considers the effects of temperature and height, assisted by Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) products. Next, we proposed a new method to determine the scale height of water vapor, which will improve the quality of vertical constraints. Finally, we determined the smoothing factor in the horizontal constraint based on Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim) products. To evaluate the advantages of the optimized technique over the traditional method, we used GPS datasets collected in Hong Kong in August 2016 to estimate the vertical distribution of water vapor density using both methods. We further validated the tomography results from the optimized technique using radiosonde products. The results show that the water vapor density quality obtained by the optimized technique is 13.8 % better below 3.8 km and 8.1 % better above 3.8 km than that obtained by the traditional technique. We computed the success rate of the tomography technique based on the Pearson product-moment correlation coefficient (PCC) and root mean square (RMS). The success rate of the optimized topography technique was approximately 10 % higher than that of the traditional tomography method.


Introduction
GPS technology has recently started being used to detect the Earth's atmosphere. Many studies have been carried out to retrieve the two-dimensional (2-D) or three-dimensional (3-D) distribution of atmospheric water vapor (Flores et al., 2000;Champollion et al., 2005;Nilsson and Gradinarsky, 2006;Esteban et al., 2013;Jiang et al., 2014;Chen and Liu, 2014). The obtained atmospheric water vapor product can be assimilated into a numerical weather prediction (NWP) model. By applying the NWP model to weather forecasting, we have discovered the usefulness of GPS tomography to estimate water vapor distribution (Jin et al., 2011;Esteban et al., 2013). Combined with the space-based GNSS (Global Navigation Satellite System) occultation technique, it can provide neutral-atmosphere products with high precision, high vertical resolution and low-cost near-real-time all-weather global coverage. In addition, it can contribute to scientific research on the ionosphere (Kursinski et al., 1997;Rocken et al., 1997;Hajj et al., 2002;Kuo et al., 2007).
In ground-based GPS meteorology, GPS signal propagation through the atmosphere is slowed, thus causing path delay on the GPS measurements, which is termed tropospheric delay (Kouba and Héroux, 2001). Zenith total delay (ZTD) is one of the most important error sources in GNSS navigation and positioning; however, it is a very reliable information source in GNSS meteorology (Jacob et al., 2007;Jin et al., 2007Falconer et al., 2009). ZTD consists of two parts: zenith wet delay (ZWD) and zenith hydrostatic delay (ZHD) (Davis et al., 1985). Usually, ZHD can be calculated with high accuracy from empirical models, and ZWD can then be easily derived from ZTD based on the formula ZWD = ZTD − ZHD. Afterward, slant wet delay (SWD) can be obtained from ZWD based on the wet Niell mapping function (Niell, 1996). Both ZWD and SWD are related to atmospheric water vapor, and thus precipitable water vapor (PWV) and slant water vapor (SWV) can be derived from ZWD and SWD using the humidity conversion coefficient (Song, 2004).
ZHD is usually estimated in GNSS meteorological research using the Saastamoinen model (Flores et al., 2000;Troller et al., 2006;Champollion et al., 2009;Perler et al., 2011;Jiang et al., 2014). The atmospheric weighted mean temperature T m is the key variable to obtain a high-precision humidity conversion coefficient (Mateus et al., 2014). T m will differ significantly as the season varies and the region changes (Jin et al., 2008). It can be determined by the surface temperature measurement, which is provided by a radiosonde product or other meteorological data analyses (Bevis et al., 1992;Wang et al., 2011).
In space-based GNSS meteorology, GNSS radio occultation (RO) is regarded as a valuable data source for atmospheric change studies (Rocken et al., 1997;Kursinski et al., 1997;Hajj et al., 2002;Beyerle et al., 2005). The Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) is housed within the University Corporation for Atmospheric Research (UCAR). The mission of the COSMIC RO is to develop the weather, climate, space weather and geodetic research (Liou et al., 2007). The University Corporation for Atmospheric Research COS-MIC Data Analysis and Archive Center (UCAR/CDAAC) supplies two different types of products from the COS-MIC mission: real-time data and postprocessed data products. Of these postprocessed products, wet atmospheric profiles (wetPrfs) offer water vapor pressure, temperature, etc. Shi and Gao (2009) compared the bias of PWV between wetPrf-derived and precise point positioning (PPP)-derived data and suggested that they have comparable accuracy levels. Kishore et al. (2011) discussed the difference in specific humidity between wetPrfs and radiosonde data. They concluded that both sources have good correlation (∼ 0.8) up to 8 km and that the humidity information of wetPrfs is reliable up to nearly 8 km. In addition, Wang et al. (2013) studied the accuracy of wetPrfs using the Radiosonde products as the reference and revealed that a global mean temperature deviation of −0.09 K and a global mean humidity deviation is less than −0.12 g kg −1 in the pressure range of 925 to 200 hPa.
To improve the accuracy of water vapor derived using the GNSS technique, we optimized several key techniques for GNSS tomography. First, we precisely derived the T m model using wetPrf profiles, and then determined the regional humidity conversion coefficient. Next, for vertical constraints, we used a new way to determine the scale height of water vapor in the exponential model. Finally, we derived the smooth-ing factors of the Gauss distance weighting function in the horizontal constraint using Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim) products. We used GPS datasets from Hong Kong in August 2016 to evaluate this new method. The results demonstrate better accuracy than those of the traditional method with radiosonde data.
The rest of this paper is organized as follows. Section 2 introduces the principles of GNSS tomography and the optimized technique for establishing the atmospheric weighted average temperature model and deriving the scale height of water vapor. Section 3 describes the data processing. Section 4 presents the validation of the optimized method, and the quality control process for the tomography results. The discussions and conclusions are given in Sect. 5.

GNSS tomographic formulation
In this section, we first introduce the GPS tomography model. We then illustrate the optimized techniques for the ZHD model and the humidity conversion coefficient determination. Finally, we present the constraint model.

Tomographic technique
To reconstruct 3-D images of water vapor density distributions, the SWV along ray paths traversing the imaged region should first be obtained from dual-frequency GNSS data. This is defined by the line integral of water vapor density along the ray path from satellite to receiver (Flores et al., 2000), as follows: where ρ w denotes the density of liquid water, s denotes the trajectory of GNSS signals in the troposphere, and ρ(s) indicates the water vapor density. Equation (1) reveals that the accuracy of water vapor density mainly depends on the quality of the SWV. Generally, ZTD can be precisely estimated using the double-difference or PPP method. ZWD can be obtained by removing ZHD from ZTD. After the humidity conversion coefficient is determined, the SWV will be computed, providing the SWD is known (MacMillan, 1995), as follows: L gradient = 1 sin (e) · tan (e) + C · (G N · cos (α) where STD and SHD are slant troposphere delay and slant hydrostatic delay, respectively; L gradient denotes the hor-izontal gradient; G N and G E are the north and east atmospheric horizontal gradients, respectively; e and α are the satellite elevation angle and the azimuth angle, respectively; C is a constant with as C = 0.003 (Chen and Herring, 1997); and denotes the humidity conversion coefficient. SWD and SHD can be projected to ZWD and ZHD based on the Niell mapping function (Niell, 1996). From Eqs.
(2) and (5), we know that the accuracy of the ZHD and the humidity conversion coefficient are the crucial aspects that affect SWV quality. Thus, it is essential to develop a high-precision ZHD model and humidity conversion coefficient.

Humidity conversion coefficient
The humidity conversion coefficient can be expressed as a function of T m . T m varies across seasons and areas and depends mainly on the surface atmospheric temperature (Bevis et al., 1994), as follows: where ρ w is the density of liquid water; k 1 , k 2 and k 3 are constants -k 1 = 77.6 K hPa −1 , k 2 = 70.4 K hPa −1 and k 3 = 3.739 × 10 5 K hPa −1 (Bevis et al., 1994); T m is the atmospheric weighted average temperature; m d and m w denote the molar masses of dry atmosphere and water vapor, respectively; R indicates the universal gas constant; P w indicates water vapor pressure in units of hPa; T is the atmospheric temperature; and h denotes the height.

Constraint model
Usually, the observation equation of the tomographic approach is rank deficient because the GPS signal cannot pass through all of the grids. Horizontal constraints, vertical constraints, priori information value constraints and boundary constraints must be added to avoid this deficiency. With these constraints, we can use an iterative reconstruction algorithm, or a noniterative reconstruction algorithm to resolve the tomography equation. The horizontal constraint is the Gauss distance weighting function (Song, 2004), as follows: where B is the horizontal smoothing; the subscripts i, j , k denote the index of voxel in 3-D space; nl and nn are the numbers of the grids in the east-west and north-south directions, respectively; d i,j,k indicates the distance between known and unknown water vapor grids; and δ denotes the smoothing factor, which will change at different levels. Section 3.3.1 explains how to estimate δ. The vertical distribution of water vapor does not follow the ideal-gas law, particularly in the lower levels. Currently, there is no accurate model function to fit the spatial distribution of water vapor. The vertical constraint of atmospheric tomography can be obtained using an exponential model (Jiang et al., 2014;Ye et al., 2016), as follows: where ρ(h) is the water vapor density at the height of h, ρ 0 is the water vapor density at the height of h 0 , and H we is the scale height of water vapor. The variables ρ 0 , h 0 and H we can usually be determined using radiosonde or COSMIC historical data. In this case, the estimated ρ(h) is only an experience value and will have a greater error than the true value. Therefore, we propose a new method to estimate ρ(h) and H we in near-real time. Based on Eq.
(2) and the Niell mapping function (Niell, 1996), ZWD can be estimated in real-time. PWV can then be obtained according to Eq. (4). The relationship between PWV and ρ(h) is established as follows: where ρ w is the density of liquid water, h 0 is the height of station, and h top is the height of tropopause. Combining Eqs. (11) and (12), we obtain the following: The parameter H we can be derived in real time using Eq. (13). Based on Eqs. (11) and (13), Eq. (14) can be utilized to establish the functional relationship in the vertical direction, as follows: where ρ i,j,k represents the water vapor value of datum voxel (i, j, k). The a priori humidity information can be used for the background field of troposphere tomography and will enhance the computing speed and tomography accuracy. The synoptic observation data include the atmospheric pressure, atmospheric temperature and relative humidity observed in the station, and the atmospheric temperature and relative humidity can be interpolated into all of the voxels using Eqs. (10) and (14). Thus, the water vapor density of every voxel can be calculated (Jiang et al., 2014).
3 Data processing 3.1 Data collection Data used to remotely sense atmospheric water vapor contain ground-based GNSS observations and meteorological data, as well as space-based COSMIC wet profiles. UCAR/CDAAC supplies two different types of products: real-time profiles and postprocessed profiles. The former can be available within a few hours and the latter can be available with a 6-week latency (http://www.cosmic.ucar.edu, last access: 20 October 2017). We selected postprocessed profiles in this study. Wet profiles (wetPrfs) are one type of COSMIC postprocessed products that are freely available for public access (http://cdaac-www.cosmic.ucar.edu/cdaac/, last access: 20 October 2017). The wetPrfs are interpolated products sampled at 100 m intervals and obtained using a nonstandard one-dimensional variation technique together with European Centre for Medium-Range Weather Forecasts (ECMWF) low-resolution analysis data from the altitude of the perigee point from the surface to a 40 km altitude (CDAAC, 2005). The average bias of temperature between wetPrfs and radiosonde is less than 0.1 K, and 70-90 % of the wetPrfs reach to within 1 km of the surface on a global basis (http://www.cosmic.ucar.edu/ro.html, last access: 20 October 2017).
We used ground-based GNSS observations and meteorological products from the Hong Kong SatRef network (https: //www.geodetic.gov.hk/, last access: 14 September 2016), from 12 continuously operating reference stations with an inter-station distance of 7 to 27 km, covering approximately 1100 km 2 . All 12 stations were equipped with LEICA GRX1200 + GNSS receivers and had data sampling rates of 5 s, as shown in Fig. 1. The meteorological data associated with each GPS station at 60 s intervals are freely available at https://www.geodetic.gov.hk/, last access: 14 Septem- horizontal grids and 17 vertical layers. A total of 8×5×17 = 680 voxels were divided in the 3-D space.
3.2 Regional weighted average temperature model Bevis et al. (1994) first put forward the global T m model using radiosonde products. Later, Wang et al. (2011) established the T m model in Hong Kong using radiosonde products. Ye et al. (2016) also assessed the relationship between T m and surface temperature based on radiosonde and COS-MIC products. However, these three models only consider the parameter of surface temperature. We propose considering the effects of temperature and height to establish a T m model using COSMIC products. The new model is given as follows (Yao et al., 2013): where a, b, c, e and f are constants that can be determined using COSMIC products; T h indicates the temperature at height h; h denotes the height; and TmN is the new model value of T m . The weighted average temperature T m is obtained using Eq. (9) with input wet pressure and temperature provided from wetPrfs. TmN can be derived using Eq. (15); its values are shown in Fig. 2. The wetPrfs described in Sect. 3.2 are used to derive the humidity conversion coefficient from Eq. (8).
As shown in Fig. 2   The statistical results comparing the model-derived and COSMIC-derived T m are given in Table 1. We provide a sum- mary of the T m deviation between radiosonde-derived and model-derived data in Table 2. As shown in Tables 1 and 2, the new T m model improves the accuracy of atmospheric weighted average temperature from the Bevis model and Wang model.

Estimating the smoothing factor
The smoothing factor δ in Eq. (10) is an uncertain parameter in the horizontal constraint. Usually, it is assigned a constant value of experience (Xia et al., 2013;Jiang et al., 2014). Because δ varies with regions and seasons and also changes with different levels of a tomography model, ERA-Interim data for Hong Kong from August 2009 to August 2015 were used to precisely estimate δ. ERA-Interim is a reanalysis of the global atmosphere covering the data period since 1989 and is continuing in real time (http://apps.ecmwf.int/datasets/ data/, last access: 20 October 2017). The specific humidity data with 60 levels of vertical spatial resolution and a minimum grid of 0.125 • × 0.125 • are publicly available. The main characteristics of the ERA-Interim system and many aspects of its performance are described in ECMWF newsletters 110, 115 and 119 (http://www.ecmwf.int/publications/ newsletters, last access: 20 October 2017). In addition, comprehensive documentation of ERA-Interim, including observation usage, is currently being prepared and will be made available at http://www.ecmwf.int/research/era, last access: 20 October 2017.
In each level, the humidity information of one grid point equals the weighted average of its neighbors (Rius et al., 1997), as follows: 0 = B 1 ρ 1 +B 2 ρ 2 +· · ·+B j −1 ρ j −1 ρ j +B j +1 ρ j +1 +· · · (14) According to the humidity information provided by ERA-Interim, Eq. (16) can be solved using the optimal parameter  Table 3 shows that the smoothing factors present a nonlinear change for increasing heights below 6 km, but do not change between 6 and 9 km. The horizontal constraint can be accurately determined based on the smoothing factor and the distance between known and unknown grids.

Vertical constraint
The purpose of GNSS tomography technique is to derive the 3-D distribution of water vapor. Thus, the accuracy of the vertical constraint will directly affect the quality of the tomography results. Because water vapor randomly varies in time and space, it is difficult to precisely probe the spatial distribution of water vapor. Traditionally, Eq. (14) was used as a vertical constraint and the parameter H we could be obtained using COSMIC or radiosonde historical data products (Ye et al., 2013(Ye et al., , 2016. Due to H we changes over time are obvious, so they need to be obtained once for each tomography epoch. In this paper, PWV was derived using Eq. (4), and H we was then derived in real time based on Eq. (13). To evaluate the accuracy of H we , the radiosonde-obtained water vapor is used as a reference to assess the water vapor calculated using Eq. (11). The statistical results are given in Table 4 using the 45 004th radiosonde station (latitude: 22.32 • N; longitude: 114.14 • E) and HKSC station (latitude: 22.32 • N; longitude: 114.14 • E) datasets from August 2016 under 10 km.
As shown in Table 4, the water vapor density derived from the H we obtained using the new technique and Eq. (11) is closer to the radiosonde-derived water vapor density. There-fore, it is more reasonable to use the H we obtained using the new technique and Eq. (14) as the vertical constraint.

Result validation and analysis
To evaluate our optimized method, we obtained ZTDs from the Hong Kong SatRef network in August 2016, based on Bernese 5.2 (nondifference) software. The ZHDs were estimated using the Saastamoinen mode. The SWV was then obtained using the Niell mapping function (Niell, 1996) and the calibrated humidity conversion coefficient. The WVLT was determined as 9.5 km from COSMIC historical data and Ye et al.'s (2016) method. Following the tomography model proposed by Flores et al. (2000), we estimated the 3-D water vapor distribution using the GPS tomography technique with the horizontal constraint from Eq. (10) and the vertical constraint from Eq. (14). The tomography equation was solved by adopting Kalman filtering. The tomography results were outputted once every 30 min. As we have limited space, Fig. 4 only shows the 3-D distribution of water vapor density on 1 and 2 August 2016. Figure 4 presents the 3-D tomographic water vapor distribution in Hong Kong for heights lower than 9.5 km. The results show that the water vapor changes significantly below 3.8 km, whereas it remains stable above 3.8 km. In addition, the water vapor is mainly concentrated below 2.6 km.

Compare the results between tomography-obtained and radiosonde-obtained data
Radiosonde products contain 3-D distribution of meteorological elements such as atmospheric temperature, atmospheric pressure, mixing ratio and relative humidity. The "wet" pressure can be obtained based on the pressure and mixing ratio and can be utilized to compute the water vapor density (Song, 2004). To verify the advantage of the optimized GPS tomography method, using radiosonde products as references, the tomography results were compared with those derived from the traditional tomography technique using the Saastamoinen dry model, a traditional humidity conversion coefficient (0.1538), a smoothing factor (Eq. 10), and H we obtained using the traditional technique and COSMIC historical products. Figure 5 compares the water vapor densities derived from radiosonde products and the traditional and optimized tomography techniques for 1 and 2 August 2016. It can be observed in Fig. 5 that the changing trends of water vapor with height across the tomography-obtained and radiosonde-obtained data have a good agreement. However,  . Water vapor densities obtained from tomography-derived and radiosonde-derived data. Rad is the water vapor density derived using radiosonde products, Trad is the water vapor density derived using the traditional tomography method, and Opti is the water vapor density derived using the optimized method. when the "inversion layer" occurs, GPS tomography cannot accurately reflect this situation. In Table 5, we present the deviation statistics for GNSS tomography-obtained and radiosonde-obtained water vapor density at heights above and below 3.8 km, using 31-day datasets from Hong Kong over the whole of August 2016. Table 5 provides the statistics values of the differences between GNSS tomography-obtained and radiosonde-obtained results. As seen from the statistical results, the root mean square (RMS) and mean values of troposphere tomography using the optimized technique is less than that based on the traditional method for altitudes below 3.8 km. In addition, compared with the radiosonde data, the test results show that the water vapor density quality obtained by the optimized technique is 13.8 % better below 3.8 km and 8.1 % better above 3.8 km than that obtained by the traditional technique. Table 5. Statistics for tomography-derived and radiosonde-derived water vapor density above and below 3.8 km (g m −3 ). Rad is the radiosonde-derived water vapor density, Opti is the optimized tomography-derived water vapor density, and Trad is the traditional tomography-derived water vapor density.

Height
Lower 3

Quality of GPS tomography technology
We also studied the differences in the entire humidity profile between the tomography-derived and radiosonde-derived results. We used the RMS and Pearson product-moment cor-  relation coefficient (PCC) as the evaluation index correlated between the two profiles. PCC is a commonly used measure of the degree of correlation of two sequences of parameters, and the mathematical model is as follows (Lee and Nicewander, 1988): (15) Figure 6 presents the PCC and RMS of tomography results (traditional and optimized) for August 2016. Here we set up a set of criteria to evaluate the tomography profile PCC > 0.90 and RMS < 2.0 g m −3 . When GPS tomography results meet these criteria, they are considered a success. According to the criteria, the success rate of the inversion is shown in Table 6. As shown in Table 6, the success rate of the optimized technique is nearly 10 % higher than that of the traditional technique, and the degree of improvement is evident. In fact, the principles of radiosonde and GPS tomography techniques are different. Radiosonde products reflect the state of the atmosphere at a certain time at the instrument's location, but GPS tomography techniques mirror the average water vapor state. Thus, it is difficult to determine an absolute standard to evaluate the success of GPS tomography results.

Conclusions
In this study, several key techniques in the GNSS tomography method were optimized to improve the accuracy of water vapor density. First, we re-established an atmospheric weighted average temperature model using COSMIC wet-Prfs. According to the spatial distributions of water vapor provided by COSMIC products, we used the exponential model to fit the vertical variation of water vapor. The exponential function is usually utilized as the vertical constraint, and we proposed a new method to compute the scale height of water vapor. We determined the smoothing factor of the Gauss distance weighting function using ERA-Interim products. Finally, we used GPS datasets from Hong Kong in August 2016 to compute the PWV and the vertical distribution of water vapor density.
To evaluate the quality of the optimized technique, we compared the optimized and traditional technique results with radiosonde-obtained water vapor. The statistical results show that the water vapor density quality obtained by the optimized technique is 13.8 % better below 3.8 km and 8.1 % better above 3.8 km than that obtained by the traditional technique. We then calculated the success rate of tomographic inversion according to PCC and RMS. The statistics show that the success rate of the optimized technique was approximately 10 % higher than that of the traditional technique.
Data availability. The Survey and Mapping Office of the Lands Department (2017), Hong Kong Special Administrative Region, provides GPS data and meteorological data. These datasets are from the Hong Kong Satellite Positioning Reference Station (SatRef) and can be made freely available for public access (https://www. geodetic.gov.hk/, last access: 14 September 2016). The 45 004 Kings Park radiosonde observations (University of Wyoming (UWYO), 2017) can be downloaded from the following website: http://weather.uwyo.edu/upperair/sounding.html, last access: 20 October 2017. In addition, the Constellation Observing System for Meteorology, Ionosphere and Climate-1 program (CDAAC, 2017) offers the COSMIC radio occultation data, which can be freely obtained from ftp://cosmic-io.cosmic.ucar.edu/cdaac/index. html, last access: 20 October 2017.
Author contributions. PX conceived and designed the work that led to the submission; SY acquired data and played an important role in interpreting the results. PJ and LP drafted and revised the manuscript, and MG approved the final version.