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.

1 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; Jin and Luo, 2009; 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., 2007, 2009; Falconer 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 COSMIC Data Analysis and Archive Center (UCAR/CDAAC) supplies two different types of products from the COSMIC 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⁻¹ 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 smoothing 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.

2 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.

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

$$\begin{array}{}\text{(1)}& \text{SWV}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot \left(\underset{s}{\int }\mathit{\rho }\left(s\right)\mathrm{d}s\right),\end{array}$$

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:

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }\text{SWD}=\text{STD}-\text{SHD}-\mathrm{\Delta }{L}_{\mathrm{gradient}},{\displaystyle }& {\displaystyle }\mathrm{\Delta }{L}_{\mathrm{gradient}}={\displaystyle \frac{\mathrm{1}}{\sin \left(e\right)\cdot \tan \left(e\right)+C}}\cdot \left(G_{\mathrm{N}}\cdot \cos \left(\mathit{\alpha }\right)\right.\\ \text{(3)}& {\displaystyle }& {\displaystyle }\left.+G_{\mathrm{E}}\cdot \sin \left(\mathit{\alpha }\right)\right),\text{(4)}& {\displaystyle }& {\displaystyle }\text{PWV}=\mathrm{\Pi }\cdot \text{ZWD},\text{(5)}& {\displaystyle }& {\displaystyle }\text{SWV}=\mathrm{\Pi }\cdot \text{SWD},\end{array}$$

where STD and SHD are slant troposphere delay and slant hydrostatic delay, respectively; ΔL_gradient denotes the horizontal 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.

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

$$\begin{array}{}\text{(6)}& \mathrm{\Pi }={\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{{\mathit{\rho }}_{\mathrm{w}}\cdot \frac{R}{m_{\mathrm{w}}}\cdot \left[\frac{k_{\mathrm{3}}}{T_{\mathrm{m}}}+k_{\mathrm{2}}-\frac{m_{\mathrm{w}}}{m_{\mathrm{d}}}\cdot k_{\mathrm{1}}\right]}},\end{array}$$

$$\begin{array}{}\text{(7)}& T_{\mathrm{m}}={\displaystyle \frac{{\int }_{h_{\mathrm{0}}}^{\mathrm{\infty }}\frac{P_{\mathrm{w}}}{T}\mathrm{d}h}{{\int }_{h_{\mathrm{0}}}^{\mathrm{\infty }}\frac{P_{\mathrm{w}}}{T^{\mathrm{2}}}\mathrm{d}h}}={\displaystyle \frac{\sum _{i=\mathrm{0}}^{i=n-\mathrm{1}}\left(\frac{P_{\mathrm{w}}^i}{T_i}+\frac{P_{\mathrm{w}}^{i+\mathrm{1}}}{T_{i+\mathrm{1}}}\right)\left(h_{i+\mathrm{1}}-h_i\right)}{\sum _{i=\mathrm{0}}^{i=n-\mathrm{1}}\left(\frac{P_{\mathrm{w}}^i}{T_i^{\mathrm{2}}}+\frac{P_{\mathrm{w}}^{i+\mathrm{1}}}{T_{i+\mathrm{1}}^{\mathrm{2}}}\right)\left(h_{i+\mathrm{1}}-h_i\right)}},\end{array}$$

where ρ_w is the density of liquid water; k₁, k₂ and k₃ are constants – k₁=77.6 K hPa⁻¹, k₂=70.4 K hPa⁻¹ and $k_{\mathrm{3}}=\mathrm{3.739}\times {\mathrm{10}}^{\mathrm{5}}$ K hPa⁻¹ (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.

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

$$\begin{array}{}\text{(8)}& B={\displaystyle \frac{{\exp }^{\frac{-d_{i,j,k}^{\mathrm{2}}}{\mathrm{2}{\mathit{\delta }}^{\mathrm{2}}}}}{\sum _{i=\mathrm{1}}^{nl}\sum _{j=\mathrm{1}}^{nn}{\exp }^{\frac{-d_{i,j,k}^{\mathrm{2}}}{\mathrm{2}{\mathit{\delta }}^{\mathrm{2}}}}}},\end{array}$$

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:

$$\begin{array}{}\text{(9)}& \mathit{\rho }\left(h\right)={\mathit{\rho }}_{\mathrm{0}}\cdot \exp \left(-{\displaystyle \frac{h-h_{\mathrm{0}}}{H_{\mathrm{we}}}}\right),\end{array}$$

where ρ(h) is the water vapor density at the height of h, ρ₀ is the water vapor density at the height of h₀, and H_we is the scale height of water vapor. The variables ρ₀, h₀ 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:

$$\begin{array}{}\text{(10)}& \text{PWV}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot \underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{top}}}{\int }}\mathit{\rho }\left(h\right)\mathrm{d}h,\end{array}$$

where ρ_w is the density of liquid water, h₀ is the height of station, and h_top is the height of tropopause. Combining Eqs. (11) and (12), we obtain the following:

$$\begin{array}{ll}{\displaystyle }\text{PWV}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot \underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{top}}}{\int }}\mathit{\rho }\left(h\right)\mathrm{d}h={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot \underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{top}}}{\int }}{\mathit{\rho }}_{\mathrm{0}}\cdot \exp \left(-{\displaystyle \frac{h-h_{\mathrm{0}}}{H_{\mathrm{we}}}}\right)\mathrm{d}h\\ \text{(11)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}\cdot H_{\mathrm{we}}}{{\mathit{\rho }}_{\mathrm{w}}}}\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{h_{\mathrm{top}}-h_{\mathrm{0}}}{H_{\mathrm{we}}}}\right)\right]\cong {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}\cdot H_{\mathrm{we}}}{{\mathit{\rho }}_{\mathrm{w}}}}.\end{array}$$

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:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{{\mathit{\rho }}_{i,j,k+\mathrm{1}}}{{\mathit{\rho }}_{i,j,k}}}={\exp }^{-\left(\frac{h_{k+\mathrm{1}-h_k}}{h_{\mathrm{we}}}\right)},\end{array}$$

where ${\mathit{\rho }}_{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². 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 September 2016. GPS datasets from 1 to 31 August 2016 were collected daily in Hong Kong. The wetPrfs in or near Hong Kong in August of 2009–2015 were downloaded.

Figure 1 Distribution of the Hong Kong SatRef sites (blue triangles) inside the tomography horizontal grid (black dotted lines) and the King's Park radiosonde station (red star). The region was discretized into an $\mathrm{8}\times \mathrm{5}\times \mathrm{17}$ cell grid for the GPS water vapor tomography. The layer heights are 0, 400, 800, 1400, 2000, 2600, …, 8600 m from ground to water vapor layer top.

The reconstruction region covered an area ranging from latitude 22.22 to 22.52^∘ N, longitude 113.85 to 114.35^∘ E, and from ground to water vapor layer top (WVLT) in height. Thus, the entire area of Hong Kong was divided into 5×8 horizontal grids and 17 vertical layers. A total of $\mathrm{8}\times \mathrm{5}\times \mathrm{17}=\mathrm{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 COSMIC 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):

$$\begin{array}{}\text{(13)}& \text{TmN}=a+b\cdot T_h+c\cdot T_h^{\mathrm{2}}+e\cdot h+f\cdot h^{\mathrm{2}},\end{array}$$

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

Figure 2 Considering height and surface temperature to establish the T_m model using wetPrf products. T₀ is the surface temperature, and TmC is the fitted atmospheric weighted average temperature obtained from COSMIC products for 2009 to 2015 in Hong Kong.

As shown in Fig. 2, the new model's T_m values agree well with the true values. To evaluate the new T_m model, its values are compared with those obtained from radiosonde and COSMIC products. Figure 3 shows the results for Hong Kong in August 2016.

Figure 3 New T_m model values are compared with those derived from COSMIC and radiosonde products in Hong Kong in August 2016. TmC is the T_m derived from COSMIC products, TmN is the T_m derived from the new model, TmB is the T_m derived from the Bevis model, TmW is the T_m derived from the Wang model, and TmR is the T_m derived from radiosonde products.

The statistical results comparing the model-derived and COSMIC-derived T_m are given in Table 1. We provide a summary of the T_m deviation between radiosonde-derived and model-derived data in Table 2.

Table 1 Summary of the T_m deviation between COSMIC-derived and model-derived (K) data.

Table 2 Summary of the T_m deviation between radiosonde-derived and model-derived (K) data.

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.

3.3 Tomography constraint condition

3.3.1 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 $\mathrm{0.125}{}^{\circ }\times \mathrm{0.125}{}^{\circ }$ 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.

Table 3 Smoothing factor derived by ERA-Interim products at different heights.

In each level, the humidity information of one grid point equals the weighted average of its neighbors (Rius et al., 1997), as follows:

$$\begin{array}{}\text{(14)}& \mathrm{0}=B_{\mathrm{1}}{\mathit{\rho }}_{\mathrm{1}}+B_{\mathrm{2}}{\mathit{\rho }}_{\mathrm{2}}+\mathrm{\cdots }+B_{j-\mathrm{1}}{\mathit{\rho }}_{j-\mathrm{1}}{\mathit{\rho }}_j+B_{j+\mathrm{1}}{\mathit{\rho }}_{j+\mathrm{1}}+\mathrm{\cdots }\end{array}$$

According to the humidity information provided by ERA-Interim, Eq. (16) can be solved using the optimal parameter search method. The search step is set to 1 and the search range is [0, 20]. The value of δ is exactly equal to the number of grid points in each level, and we defined the mean of δ as the smoothing factor of the level. Table 3 lists the δ values at different heights using ERA-Interim data for Hong Kong from August 2009 to August 2015.

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.

3.3.2 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, 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.

Table 4 Statistical results from Eq. (11)-derived and radiosonde-derived PWV (g m⁻³). RWV is the water vapor density obtained from the radiosonde product; TWV is the water vapor density derived from the H_we obtained by the traditional method using Eq. (11); NWV is the water vapor density derived from the H_we obtained by the new method using Eq. (11).

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. Therefore, it is more reasonable to use the H_we obtained using the new technique and Eq. (14) as the vertical constraint.

4 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 3-D tomographic water vapor distribution in Hong Kong 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.

4.1 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.

Figure 5 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.

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, 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 Statistics for tomography-derived and radiosonde-derived water vapor density above and below 3.8 km (g m⁻³). 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.

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.

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

$$\begin{array}{}\text{(15)}& \text{PCC}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(X_i-\overline{X_i}\right)\left(Y_i-\overline{Y_i}\right)}{\sqrt{\sum _{i=\mathrm{1}}^N{\left(X_i-\overline{X_i}\right)}^{\mathrm{2}}}\cdot \sqrt{\sum _{i=\mathrm{1}}^N{\left(Y_i-\overline{Y_i}\right)}^{\mathrm{2}}}}}.\end{array}$$

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⁻³. 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.

Figure 6 Time series of PCC and RMS for August 2016. Opti is the optimized tomography-derived water vapor density, and Trad is the traditional tomography-derived water vapor density.

Table 6 Statistical results of PCC and RMS for August 2016 (%).

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.

5 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 wetPrfs. 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.