Reconstruction of Wet Refractivity Field Using an Improved Parameterized Tropospheric Tomographic Technique

Abstract

In most previous studies of tropospheric tomography, water vapor is assumed to have a homogeneous distribution within each voxel. The parameterization of voxels can mitigate the negative effects of the improper assumption to the tomographic solution. An improved parameterized algorithm is proposed for determining the water vapor distribution by Global Navigation Satellite System (GNSS) tomography. Within a voxel, a generic point is determined via horizontal inverse distance weighted (IDW) interpolation and vertical exponential interpolation from the wet refractivities at the eight surrounding voxel nodes. The parameters involved in exponential and IDW interpolation are dynamically estimated for each tomography by using the refractivity field of the last process. By considering the quasi-exponential behavior of the wet refractivity profile, an optimal algorithm is proposed to discretize the vertical layers of the tomographic model. The improved parameterization algorithm is validated with the observational data collected over a 1-month period from 124 Global Positioning System (GPS) stations of Hunan Province, China. Assessments by GPS, radiosonde, and European Centre for Medium-Range Weather Forecasts (ECMWF) ReAnalysis 5 (ERA5) data, demonstrate that the improved model outperforms the traditional nonparametric model and the parameterized model using trilinear interpolation. In the assessment by GPS data, the improved model performs better than the traditional model and the trilinear parameterized model by 54% and 10%, respectively. Such improvements are 31% and 10% in the validation by radiosonde profiles. In comparison with the ERA5 reanalysis, the improved model yields a minimum overall root mean square (RMS) error of 8.94 mm/km, while those of the traditional and trilinear parametrized models are 10.79 and 9.73 mm/km, respectively. The RMS errors vertically decrease from ~20 mm/km at the bottom to ~5 mm/km at the top layer.

1. Introduction

Water vapor in the troposphere represents a mere fraction of the total atmospheric volume but is strongly associated with climate change, atmospheric radiation, weather pattern, and hydrologic cycle [1,2,3,4]. Accurate information on water vapor not only leads to a better understanding of the aforementioned fields but also to an enhanced natural hazard mitigation (e.g., floods and landslides) because water vapor observations are crucial for initializing the numerical weather prediction (NWP) models [5,6,7]. Nevertheless, atmospheric water vapor is one of the poorly described components in the atmosphere because of its highly spatiotemporal variability [8,9,10].

The global navigation satellite system (GNSS)-based tropospheric tomography has become a powerful technique for retrieving the water vapor fields with both high spatial and temporal resolutions owing to the rapid development of the GNSS [11,12,13,14,15,16]. The first research work was carried out by Flores et al.; they reconstructed the 3D wet refractivity fields with the tomographic method by using rays from a global positioning system (GPS) network in Hawaii, USA [11]. After this successful trial, a number of significant studies have been performed in terms of the theoretical models and experimental analysis for GNSS-based tropospheric tomography [5,16,17,18,19,20,21,22]. The vital significance of tomographic water vapor products for scientific research (e.g., heavy precipitation monitoring [23,24,25], precise point positioning (PPP) augmentation [26], and assimilation into NWP models [27]) has justified the various efforts in tomographic modeling.

The tomographic space is normally partitioned into many 3D closed voxels assuming that the water vapor of each voxel is constant and evenly distributed during the modeled time period. The wet refractivity field can be retrieved from a large number of slant wet delays (SWDs) penetrating the probed space from various directions via the tomographic technique. The number of crossing SWDs per voxel is dependent on the geometry defined by the constellation of GNSS satellites, the geographic distribution of ground-based receivers, and the integration time and on the model resolution [28]. The tomographic equation is often ill-conditioned, and some voxels are underdetermined because having enough GNSS satellites and ground sites to allocate sufficient rays for each voxel is impossible. The following are the four ways to mitigate the ill-posed problem in the tomographic equation: (1) Addition of intervoxel constraints (e.g., horizontal and vertical constraints) to tomographic equations [11,16,29]; (2) assimilation of non-GNSS measurements (e.g., radiosonde, NWP, and radio occultation) [8,17,30]; (3) optimization of voxel discretization [16,20,31]; and (4) usage of advanced algorithms, such as singular value decomposition, Kalman filter approach, and algebraic reconstruction technique, to solve the tomographic equations [11,32,33].

As previously stated, considerable progress has been achieved in the tropospheric tomography in the past decades. In most previous studies, a critical deficiency in the tomographic model is to assume that the water vapor content in each voxel follows a homogeneous distribution. Water vapor greatly varies with space in the voxel, particularly in the vertical direction. The negative effects caused by the improper assumption can be mitigated by applying a high resolution. This approach will increase the computational complexity and the effect of intervoxel constraints on the solutions. In the field of 2D image reconstruction, Andersen and Kak [34] applied the discrete approximation to the ray integrals of a continuous image by using bilinear interpolation. Their study proved the superior performance of the continuous image representation with bilinear elements over the traditional pixel-based method. Ding et al. [31] reported a method to determine the water vapor density at a certain point via inverse distance weighted (IDW) interpolation in the horizontal direction for the troposphere tomography. However, water vapor is assumed to have no vertical variations within a layer, which is unreasonable in cases of the large thickness of the voxel layer or strong vertical changes in water vapor. Perler et al. [33] proposed a new voxel parameterization method by modeling the wet refractivity at a certain point by utilizing trilinear/spline functions from its eight adjacent voxel nodes. The new parameterized tomographic model is shown to be a valid means to reduce the impacts of discretization while negligibly increasing the computational costs. Nevertheless, bilinear/spline interpolations adopted in the parameterization do not consider the physical characteristic of the water vapor changes. Chen et al. [35] applied the method of Perler et al. [33] in ionospheric tomography with modified interpolations, showing significantly better performance than the traditional nonparametric method. Compared with the troposphere, the spatiotemporal distribution of the ionospheric electron density is more stable, thus constant parameters were used in the interpolations. On the basis of the studies of Perler et al. [33] and Chen et al. [35], we developed an improved parameterized algorithm to refine the tropospheric tomographic model to enhance the performance of the wet refractivity reconstruction. Horizontal IDW interpolation and vertical exponential interpolation are introduced to our improved model, and their parameters are dynamically estimated for every tomographic process. In addition, an optimal algorithm is proposed to determine the vertical discretization of the tomographic model.

Section 2 describes the methodology of the improved parameterized water vapor tomography. Section 3 presents the voxel discretization and the datasets exploited to carry out the tomographic experiments. The assessments of the parameterized tomographic model by GPS, radiosonde, and European Centre for Medium-Range Weather Forecasts (ECMWF) ReAnalysis 5 (ERA5) data, are also shown in this section. Finally, Section 4 provides the conclusions and outlook of this study.

2. Retrieval of the Wet Refractivity Field with Improved Parameterized Tomography

GNSS signals will be significantly delayed due to the refraction of the neutral atmosphere as they travel through the troposphere. The tropospheric delay is normally divided into 2 components: A hydrostatic part caused by the neutral hydrostatic atmosphere and a wet part induced by the water vapor. At present, the zenith tropospheric delay (ZTD) can be estimated with a high accuracy of several millimeters. High-accuracy zenith hydrostatic delay (ZHD) can be attained using empirical models with surface meteorological observations; thus, the zenith wet delay (ZWD) can be extracted from ZTD minus ZHD. The SWD along the ray path from a receiver to a satellite can then be derived as follows:

$$\mathrm{SWD}=\left(\mathrm{ZTD}-\mathrm{ZHD}\right)f\left(z\right)+\frac{\partial f}{\partial z}\left[G_{NW}·\cos \left(\varphi \right)+G_{EW}·\sin \left(\varphi \right)\right]+R$$

(1)

where $z$ and $\varphi$ are the satellite zenith distance and azimuth angle, respectively; $f$ refers to the wet mapping function (global mapping function is used here); $G_{NW}$ and $G_{EW}$ are the components of the wet delay gradient in the north–south and east–west directions, respectively; and $R$ refers to the post-fit residuals. The exploitation of the post-fit residuals can mitigate the adverse effects of only using the gradient terms for accounting for the anisotropy of the local troposphere [17].

The relationship between SWD and wet refractivity along the signal from a satellite to a receiver can be expressed by:

$$\mathrm{SWD}={{\displaystyle \int }}_l{N}_w{\mathrm{d}}_l$$

(2)

where $N_w$ is the wet refractivity, and $l$ is the propagation path of the signal through the troposphere. Given that the effect of a ray bending to SWD can be neglected for elevations greater than 15° [36], $l$ is usually taken as a straight line in the tomography. The model space is discretized into many voxels to reconstruct the wet refractivity field from the massive SWDs interweaving in the troposphere across different directions (Figure 1). The water vapor distribution is generally assumed to be homogeneous for each voxel over the reconstruction period. In this case, each SWD is approximately equal to the sum of the product of wet refractivity and the length of the ray path crossing each voxel. Therefore, Equation (2) can be approximated by:

$$\mathrm{SWD}={\displaystyle \sum }_{i=1}^n{N}_w^i{d}_i$$

(3)

where $n$ is the number of voxels along the SWD ray path, $N_w^i$ is the wet refractivity in voxel $i$, and $d_i$ is the intercept of ray by voxel $i$.

In the parameterized tomographic model, the wet refractivity of each voxel is no more regarded as invariable but varies with the position. The wet refractivity of a generic point is expressed by a weighted mean of the wet refractivity values at the 8 voxel nodes, where the point is located. For example, Figure 1 demonstrates that the wet refractivity of any point along P1–P5 can be determined from the 8 nodes (i.e., ${\mathrm{N}}_1,{\mathrm{N}}_2,\cdots ,{\mathrm{N}}_8$) of voxel 4. Accordingly, the SWD can be expressed as an integral of the wet refractivities at the voxel nodes. The integral can hardly be analytically expressed; thus, the Newton–Cotes quadrature is used to solve the integral [33]. Figure 1 displays that the integral of wet refractivity along P1–P5 can be discretized as follows:

$${{\displaystyle \int }}_{\mathrm{P}1}^{\mathrm{P}5}{N}_w\mathrm{d}l=\frac{D_{\mathrm{P}1\mathrm{P}5}}{90}\left\{7\left[N_w\left(\mathrm{P}1\right)+N_w\left(\mathrm{P}5\right)\right]+32\left[N_w\left(\mathrm{P}2\right)+N_w\left(\mathrm{P}4\right)\right]+12N_w\left(\mathrm{P}3\right)\right\}$$

(4)

where $D_{\mathrm{P}1\mathrm{P}5}$ is the intercept of ray $l$ by voxel 4; $\mathrm{P}1$, $\mathrm{P}2$, $\mathrm{P}3$, $\mathrm{P}4$, and $\mathrm{P}5$ are 5 equidistant points on the intercept; and the 4 constant values (i.e., 90, 7, 32, and 12) are coefficients for the 4-order Newton–Cotes quadrature formula. Perler et al. [33] adopted trilinear and bilinear/spline functions to interpolate the wet refractivity of point $\mathrm{P}i$. In this work, an improved parameterization method was developed by considering the characteristic of the water vapor changes. The wet refractivity vertically follows the exponential distribution [8], thus taking $\mathrm{P}3$ as an instance, and its wet refractivity can be vertically interpolated by using points ${\mathrm{V}}_1$ and ${\mathrm{V}}_2$:

$$N_w\left(\mathrm{P}3\right)=\frac{h_{\mathrm{P}3}-h_{\mathrm{V}1}}{\left(h_{\mathrm{V}2}-h_{\mathrm{V}1}\right)·e^{\alpha \left(h_{\mathrm{V}2}-h_{\mathrm{P}3}\right)}}{N}_w\left(\mathrm{V}2\right)+\frac{h_{\mathrm{V}2}-h_{\mathrm{P}3}}{\left(h_{\mathrm{V}2}-h_{\mathrm{V}1}\right)·e^{\alpha \left(h_{\mathrm{V}1}-h_{\mathrm{P}3}\right)}}{N}_w\left(\mathrm{V}1\right)$$

(5)

where $h_{\mathrm{V}1}$, $h_{\mathrm{V}2}$, and $h_{\mathrm{P}3}$ refer to the altitudes of points $\mathrm{V}1$, $\mathrm{V}2$, and $\mathrm{P}3$, respectively; and $\alpha$ is a parameter describing the exponential variation of the wet refractivity, and it can be estimated from the following expression:

$$N_w\left(h_i\right)=N_w\left(h_0\right)e^{\alpha \left(h_{\mathrm{i}}-h_0\right)}$$

(6)

where $h_0$ represents the elevation of the lower surface of the vertical layer, and $h_i$ represents the elevation of a generic point within the layer. Variable $\alpha$ is estimated for each voxel by using the wet refractivity profiles of the last tomographic period to improve the modeling performance. The initial profiles are derived from the National Centers for Environmental Prediction (NCEP) FNL Analysis products.

The wet refractivities of ${\mathrm{V}}_1$ and ${\mathrm{V}}_2$ are determined by the IDW interpolation:

$$N_w\left({\mathrm{V}}_i\right)=\frac{{{\displaystyle \sum }}_{j=1}^4\frac{1}{{\left(d_j\right)}^u}·N_w\left({\mathrm{N}}_j\right)}{{{\displaystyle \sum }}_{j=1}^4\frac{1}{{\left(d_j\right)}^u}}$$

(7)

where $d_j$ is the distance between ${\mathrm{V}}_i$ and ${\mathrm{N}}_j$; $u$ is the power of the distance; and $N_w\left({\mathrm{N}}_j\right) \left(j=1,2,3,4\right)$ represent the 4 neighboring voxel nodes of the grid surface, where point ${\mathrm{V}}_i$ is located. For example, Figure 1 shows that the 4 surrounding nodes of point ${\mathrm{V}}_1$ are ${\mathrm{N}}_1$, ${\mathrm{N}}_2$, ${\mathrm{N}}_3$, and ${\mathrm{N}}_4$. Here, we propose to estimate $u$ for each tomographic process to refine the modeling. In each voxel layer, $u$ is estimated from Equation (7) by using the wet refractivity field of the last tomographic period.

A large sparse system of linear equations is established by collecting all the SWD measurements over the tomographic period (30 min in this study):

$$y=Ax$$

(8)

where $y$ is a column vector with a set of SWDs, $x$ is the unknown parameter vector that consists of the wet refractivities of all voxel nodes, and $A$ is a matrix with elements denoting the contributions of $x$ on the SWDs. An inversion algorithm should be carried out to solve the unknowns. However, not all the voxels have ray crossings in most cases; thus, design matrix $A$ in Equation (8) is a large sparse matrix. To overcome the problem of ill-posedness, the horizontal and vertical constraints were imposed to regularize the linear least-square inversion. These constraints were added according to Equations (4) and (6). The a priori wet refractivity field provided by the National Centers for Environmental Prediction Final (NCEP FNL) analysis dataset was used to constrain the solution. The Helmert variance component estimation method was adopted to determine the weight of the a priori information for the tomographic solution [8].

However, the tomographic solution obtained from Equation (8) was just an approximate wet refractivity field. We thus further implemented the multiplicative algebraic reconstruction technique (MART) to improve the least square solution from Equation (8) due to its advantage of converging to a good result within a short processing time [16,32]. In the $k$th MART iteration, the ratio between the observed y and reconstructed $\langle A,x^{k-1}\rangle$ is computed to produce corrections for involved voxel nodes. Given the generic $i$th ray and the generic $j$th voxel node, the $x_j^k$ wet refractivity after the $k$th iteration is calculated as follows:

$$x_j^k=x_j^{k-1}·{\left(\frac{y_i}{\langle A_i,x^{k-1}\rangle }\right)}^{\frac{\lambda {\mathrm{A}}_{ij}{x}_j^{k-1}}{{{\displaystyle \sum }}_{j=1}^n{\mathrm{A}}_{ij}{x}_j^{k-1}}}$$

(9)

where $\lambda$ is the relaxation parameter (an empirical value of 0.9 used here). The wet refractivity field solved by Equation (8) was used as the initial for the iteration. The iteration was terminated when the standard deviation of the differences between GNSS-estimated and tomographically reconstructed SWDs was less than 0.5 mm. For cases the MART solution was unable to converge to 0.5 mm, the maximum iterations were set to 50. An accurate wet refractivity field will be obtained after performing the combined reconstruction algorithm [30,36,37,38]. Tomographic results solved from the parameterized method were wet refractivities of the voxel nodes. The vertical interpolation in Equation (5) and horizontal interpolation in Equation (7) must be implemented using the wet refractivities of the 8 voxel nodes of that voxel to obtain the wet refractivity of a generic point within a voxel. In the traditional tomography, the wet refractivity of a generic point is equal to that of its located voxel.

3. Validation of the Improved Parameterized Tomography

3.1. Experiment Description and Voxel Discretization

Various tests have been conducted to validate the performance of the proposed improved parameterized tomographic model. The tomographic experiment is carried out based on GPS observations collected from 124 stations with an average separation distance of 41 km (Figure 2) from the CORS network of Hunan Province, China. The time span of the GPS data used in the tests was from 1 to 30 June 2018, which is the most humid month in that year in Hunan. Tomography was consecutively implemented with an interval of 30 min. The SWDs from HNRC (113.34°E, 25.54°N, 499.061 m), SYDK (110.61°E, 27.03°N, 321.878 m), and XTXX (112.51°E, 27.75°N, 70.008 m) stations (black triangles shown in Figure 2) were excluded in the input dataset; they were used for self-consistency validation purposes. Most GPS stations were equipped with Trimble or Leica receivers and had a typical sampling rate of 30 s. In this work, Bernese 5.2 software was exploited to estimate the ZTDs with the PPP technique [39]. The ZTDs and horizontal gradients were estimated every 30 min and 12 h, respectively, while the global mapping function (GMF) was adopted [40] in the estimation. The comparison with radiosonde measurements revealed an accuracy of ~9 mm for our estimated ZTDs [41]. The quality-assured atmospheric profiles observed at 3 radiosonde stations (blue diamonds in Figure 2) from the Integrated Global Radiosonde Archive (IGRA) [42] will be used to validate the tomographic solutions.

The tomographic region covered from 108.85°E–114.05°E in longitude and 24.85°N–30.05°N in latitude. Radiosonde profiles revealed that wet refractivities approached zero at the altitude of 11 km. Therefore, the selection of 11 km as the top boundary of the tomographic region in Hunan and regarding the atmosphere above 11 km as dry air was reasonable. The water vapor variations in the latitude and longitude directions were comparable; thus, a uniform resolution of 0.4° in the 2 horizontal directions was adopted. The water vapor content rapidly decreased with altitude and was negligible in the upper troposphere. Considering the quasi-exponential behavior of the wet refractivity profile, we derived the following expression to determine the vertical layer distribution:

$$h_i=\{\begin{array}{c}{h}_{min}+ln\left(\frac{i}{n}{e}^{\alpha \left(h_{max}-h_{min}\right)}+\frac{n-i}{n}\right)/\alpha i=1,\cdots ,n-1 \\ h_{max} i=n \end{array}$$

(10)

where $h_i$ denotes the altitude of the upper boundary of layer number $i$, $n$ is the total number of vertical layers, $h_{min}$ is the minimal surface altitude within the target area, $h_{max}$ is the top height of the tomographic region, and $\alpha$ is the exponential variation parameter that can be determined using (6) with history radiosonde data. Within the bottom and uppermost layer, this expression was established to distribute the intermediate layers for ensuring that the integral of the wet refractivity (i.e., ZWD) in each layer was approximately constant.

In this study, $h_{min}$ and $h_{max}$ were set as 0 and 11 km, respectively. A value of −0.28 was estimated for $\alpha$ by using the historical radiosonde profiles collected over the whole month of June 2017. Flores et al. [11] suggested that the thickness of a vertical layer should be more than 350 m; otherwise, the noise will affect the tomographic solutions. Accordingly, the total number of layers was determined as 10 to ensure that the thickness of the lowest layer was larger than 350 m. Finally, 10 nonuniform layers were discretized with their thicknesses of 358, 398, 448, 513, 598, 719, 902, 1209, 1842, and 4013 m. The SWDs with elevation angles <10° were rejected in the tomography to minimize the multipath effects. Three schemes were used in the tomographic modeling to assess the performance of the improved method.

Tomo-I: Using the traditional nonparametric method that water vapor was assumed to have a homogeneous distribution within each voxel.

Tomo-II: Using the trilinear parameterization method developed in Perler et al. [33]. The bilinear/spline approach was not included here since it has a worse performance than the trilinear one in the assessment with real data. The wet refractivity at a generic point within a voxel was determined by trilinear interpolation from the wet refractivities at the 8 nodes of that voxel.

Tomo-III: Using the improved parameterized method developed in this study. As previously stated, the wet refractivity of any point within a voxel was expressed via vertical exponential interpolation and horizontal IDW interpolation by using the 8 refractivity values of the voxel corners.

3.2. Self-Consistency Validation by GPS Data

As previously mentioned, 3 GPS stations (i.e., HNRC, SYDK, and XTXX) were excluded from the tomographic experiment; however, they were adopted for self-consistency validation. The SWD along a generic ray path was calculated by an integral of wet refractivities with respect to its propagation path by using the tomographic results. The tomographic SWDs were then directly compared with those estimated from GPS measurements. Figure 3 exhibits the overall statistical results derived from the 3 stations during the period of 1–30 June 2018. The SWD derived from Tomo-I performed badly because its root mean square (RMS) error of 24.68 mm was approximately twice those of schemes Tomo-II and Tomo-III. With regard to the mean bias, the 3 schemes consistently yielded positive values in the range of 1.50 to 3.50 mm. This phenomenon was likely due to the neglect of the troposphere above altitude 11 km in the tomographic model as GPS-derived SWDs contained a small portion from the space above 11 km. RMS errors obtained from both the parameterized schemes were much smaller than the nonparametric scheme. The improvements attained by the parameterized method were 49% and 54% for schemes Tomo-II and Tomo-III, respectively. Tomo-III achieved the best performance with an RMS error of 11.40 mm, which corresponds to an improvement of approximately 10% with respect to Tomo-II.

Figure 4a shows the RMS errors of the SWD comparison at 8 different elevation intervals. The RMS errors quickly decreased with the increase in elevation in all comparisons. The RMS errors of the SWD differences for elevations between 10° and 20° were 2.7, 4.0, and 4.5 times those for elevations 80°–90° for Tomo-I, Tomo-II, and Tomo-III, respectively. The significant increase of the RMS error with elevation decrease occurred because the GPS rays will cost a longer time to travel through the troposphere at a low elevation, thereby leading to a larger wet delay. For this reason, Figure 4b further displays the change of relative RMS with elevations. The relative RMS was defined as the averaged GPS-estimated SWD divided by the corresponding RMS error. Tomo-I obtained relative RMS varying from 3% to 5%, whilst much smaller relative RMS values in the range of 1.5% to 2% were yielded for parameterized methods Tomo-II and Tomo-III. The relative RMS of Tomo-I, in general, increased with the elevation increase; however, slight decreases were found for Tomo-II and Tomo-III. This finding shows that the parameterized method was more robust than the traditional one because no evident changes in performance were observed at different elevations.

Table 1. Biases and RMS errors of the differences between GPS-estimated SWDs and tomographic SWDs using schemes Tomo-I, Tomo-II, and Tomo-III.

Station	Tomo-I		Tomo-II		Tomo-III
	Bias (mm)	RMS (mm)	Bias (mm)	RMS (mm)	Bias (mm)	RMS (mm)
HNRC	5.20	24.37	0.82	12.33	2.87	11.24
SYDK	0.54	24.12	6.12	16.40	0.51	11.37
XTXX	4.50	25.49	0.68	9.16	0.86	9.06

further illustrates the respective statistics of the reconstructed SWDs by the 3 schemes at the 3 GPS stations. Consistent with the overall statistics, Tomo-III performed best at all the stations, followed by Tomo-II.

3.3. Validation of the Tomographic Solutions by Radiosonde Profiles

Although the tomographic SWDs were in agreement with the GPS-estimated ones, the vertical profiles of wet refractivity were not necessarily correctly modeled. In this section, we further exploit radiosonde data to assess the tomographic wet refractivity profiles. Figure 2 shows that 3 radiosonde stations were located in Hunan Province. However, only 2 stations (i.e., RSCZ and RSHH) can provide observations for the period of the tomographic experiment. The measured radiosonde and the reconstructed tomographic wet refractivity profiles were resampled to heights with an interval of 200 m to conduct a straightforward comparison. In the traditional method, the wet refractivity of an arbitrary point was equal to that of the voxel where the point was located. The matchup wet refractivities of the parameterized methods were obtained by two steps: (1) Searching the voxel where the point was located; and (2) interpolating the wet refractivities of the 8 voxel nodes to the point by trilinear interpolation for Tomo-II or exponential/IDW interpolation for Tomo-III. The comparison between radiosonde and Tomo-I yielded a bias of 0.69 mm/km and an RMS error of 10.17 mm/km, respectively. The bias and RMS error of the wet refractivity profiles for Tomo-II were 0.27 and 7.81 mm/km, respectively. In scheme Tomo-III, the obtained bias and RMS error were −0.33 and 7.00 mm/km, respectively. The overall statistics showed that the tomographic profiles by Tomo-III have a great agreement with those observed by the radiosonde.

Figure 5 exhibits the change of RMS errors and relative RMS at various altitudes. Here, the relative RMS was defined as the averaged radiosonde-observed wet refractivity divided by the corresponding RMS of the layer. Tomo-III consistently showed the optimal performance with the RMS error decrease from ~10 mm/km at the bottom layers (0–1 km) to ~3 mm/km at the upper layers (9–11 km). The RMS errors of Tomo-I were larger than those of Tomo-II at various altitudes, particularly exceeding 30 mm/km at the bottommost layer. The relative RMS values for Tomo-III increased from 8% at the lowest layer to 443% at the uppermost layer with altitude. Figure 5b demonstrated that the relative RMS exceeded 100% at an altitude above 8 km for all the schemes. In the worst scheme Tomo-I, the relative RMS approached 1100% at the uppermost layer. This finding indicates the difficulty of tomography in reconstructing water vapor profiles of high-altitude layers. The water vapor content in the upper layers is small, and a minor error in the tomographic modeling would cause relatively large discrepancies in wet refractivity for the top layers.

3.4. Comparison of the Wet Refractivity Fields between Tomography and ERA5 Reanalysis

The limited spatial coverage of the benchmark datasets of the above 2 assessments hampered the comprehensive understanding of the tomographic solutions. ERA5 was the latest (5th generation) European Centre for Medium-Range Weather Forecasts atmospheric reanalyses of the global climate, which will replace its predecessor ERA-Interim within several years [43]. ERA5 reanalysis has been greatly upgraded in the spatiotemporal resolution and assimilation method compared with ERA-Interim [44]. ERA5 can provide hourly atmospheric parameters at 37 pressure levels from 1000 to 0.1 hPa at horizontal grids of 0.25° × 0.25°. ERA5 reanalysis offers us a chance to assess our tomographic solutions from the perspective of high spatial and temporal resolutions. Tomo-I, Tomo-II, and Tomo-III yielded biases of 3.27, 3.67, and 2.79 mm/km, respectively, by stacking 1-month comparison data of all voxels. The RMS errors of 10.79, 9.73, and 8.94 mm/km were obtained for Tomo-I, Tomo-II, and Tomo-III, respectively. Tomo-III had the optimal overall agreement with ERA5 data.

Figure 6a–c present the spatial distribution of the RMS errors of the wet refractivity differences between ERA5 and tomography over the study region. At each horizontal grid, the RMS error was calculated considering all the vertical voxels over this grid. The RMS errors of Tomo-I, Tomo-II, and Tomo-III varied from 7.0 mm/km to 16.8 mm/km, 5.9 mm/km to 15.8 mm/km, and 6.0 mm/km to 11.0 mm/km, respectively, depending on the location. In Tomo-I, large RMS errors exceeding 15 mm/km occurred in the southeast portion of the study area. In Tomo-II, the north part of the study area achieved a worse performance with RMS errors greater than 12 mm/km, with reasons unknown. In Tomo-III, a majority of the area was populated with RMS errors less than 10 mm/km, thereby showing an evidently enhanced consistency with ERA5. Figure 6d–f display the RMS differences between each of the two schemes. Figure 6d exhibits that the RMS errors of Tomo-I are larger than those of Tomo-II in most regions, except for the north above 29°N where the RMS differences of −1 to −3 mm/km are found. Tomo-III significantly performed better than Tomo-I because positive RMS differences were observed in the vast majority of the study area (Figure 6e). In Figure 6f, positive values could be observed everywhere, thereby demonstrating the improvements brought by Tomo-III versus Tomo-II for the parameterized method. The large RMS errors occurred in the boundary regions as less GPS rays interweaved in the troposphere. The southwest consistently obtained relatively good performance in all the 3 schemes, while no GPS sites were located there. The tomographic solutions of voxels over the southwest were highly dependent on the initial information due to the lack of GPS ray crossings. In this study, the NCEP FNL analysis was used to provide the a priori water vapor fields for the tomography. The NCEP FNL analysis had a good consistency with the ERA5 reanalysis, thereby leading to a low RMS error in the southwest. Figure 7 further shows the RMS error and relative RMS at 10 vertical layers. The RMS errors basically decreased with the increase in altitude from ~20 mm/km at the bottom to ~5 mm/km at the top. The relative RMS values decreased from approximately 20% at the bottom to over 90% at the top. The Tomo-III again outperformed the other 2 schemes at various vertical layers.

4. Conclusion and Outlook

The water vapor within each voxel is assumed to have homogeneous distribution in the tomographic modeling, which is unreasonable for cases with coarse voxel discretization and highly variable water vapor changes in the space. The parameterization of voxels can reduce the effects of discretization. In this study, we presented an improved parameterized algorithm for tropospheric tomography and validated its superiority in several tests. In the improved algorithm, the wet refractivity of a generic point is expressed via vertical exponential interpolation and horizontal IDW interpolation by using the eight refractivity values at the voxel nodes in which the point is located. The parameters involved in exponential and IDW interpolation are dynamically estimated for each tomography by using the wet refractivity field of the last process. In addition, an optimal expression is derived to discretize the vertical layers of the tomographic model, considering the quasi-exponential behavior of the wet refractivity profile. Various tomographic tests were carried out using SWD measurements estimated from 124 GPS sites of Hunan, China, over the whole month of June of 2018 to examine the performance of the improved parameterized method. Tomographic tests using the traditional nonparametric model and parameterized model with trilinear interpolation were also performed for a straightforward comparison with the improved model.

The tomographic water vapor results were fully evaluated with independent datasets derived from GPS, radiosonde, and ERA5 reanalysis. All assessments demonstrated the better performance of the improved model over the nonparametric model and the trilinear parameterized model. In the assessment by GPS-inferred SWD measurements, the improved model outperformed the traditional model and the trilinear parameterized model by 54% and 10%, respectively. In the evaluation of the wet refractivity profiles by radiosonde, the improved model yielded an RMS error of 7.00 mm/km with respect to 10.17 and 7.81 mm/km for the traditional model and the trilinear parameterized model, respectively. The RMS error vertically decreases from ~10 mm/km at the lowest layers (0–1 km) to ~3 mm/km at the uppermost layers (9–11 km). The relative RMS values increase from 8% (from the bottom) to 443%. The improved model achieved an optimal consistency with ERA5 reanalysis data with an overall RMS error of 8.94 mm/km. The RMS errors of the refractivity differences between ERA5 and the improved model vary from 6.0 mm/km to 11.0 mm/km throughout the study area. In the vertical profiles, the relative RMS increases from ~20% at the bottom to ~90% at the altitude of 9 km. Both assessments of the vertical profiles by radiosonde and ERA5 reanalysis reveal the difficulty of tomography in the reconstructing wet refractivity of altitudes above 8 km because the relative RMS may reach up to 1000% in the uppermost layer.

The high-quality water vapor fields retrieved by the tomography have many application potentials (e.g., atmospheric research, rainfall monitoring and forecasting, and GNSS positioning). In our study, the improved voxel parameterization methods have been developed to refine the spatial modeling. Future work will focus on the parameterized tomographic modeling with high temporal resolution (e.g., 5 min). The improvement in the standard and precise point positioning brought by the tomographic SWDs will be examined. The assimilation of the tomographic refractivity fields into a numerical prediction model to enhance the capability of precipitation forecasting should be considered in the future.