An optimal tropospheric tomography approach with the support of an auxiliary area

Among most current tropospheric tomography studies, only the signals crossing out from the top boundary of the tomographic area are used for reconstructing the three-dimensional water vapour field, while signals penetrating from the side faces of the tomographic body are ignored as invalid information. Such a method wastes the valuable Global Navigation Satellite System (GNSS) observations and decreases the utilisation efficiency of GNSS rays. This is the focus of this paper, which tries to effectively use signals penetrating from the side faces of the tomographic body for water vapour reconstruction. An optimised tropospheric tomography method is proposed using an auxiliary area. The top height of the tomography body is determined based on the average water vapour distribution derived from the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) radio occultation (RO) products. In addition, the coefficients of a negative exponential function between the adjacent layers for vertical constraints are fitted using the COSMIC RO profiles. Thirteen GPS stations are selected in the CORS Network of Texas to perform the tomographic experiment and validate the performance of the proposed method at 00:00 and 12:00 UTC daily using the radiosonde data for a period of 15 days. Compared to the conventional method, the accuracy of the reconstructed water vapour information derived from the proposed method is increased by 14.37 % and 16.13 %, respectively, in terms of mean root mean square (rms) and mean absolute error (MAE). The tomographic results obtained from the proposed method are further validated with the slant water vapour (SWV) data derived using the GAMIT (GNSS processing software package). Results show that the rms and MAE accuracy of SWV values has been improved by 18.18 % and 27.62 %, respectively, when compared to the conventional method.


Introduction
Water vapour is one of the most important factors affecting atmospheric dynamics, evolution of weather systems, and global climate change.Therefore, knowledge of the water vapour distribution at a high spatio-temporal resolution is a vital prerequisite for the prediction of medium-to smallscale weather events as well as the monitoring of climate change (Zhang et al., 2015).Traditionally, radiosonde and microwave radiometers have been used for acquiring water vapour information at a low spatio-temporal resolution with high costs and the techniques do not work in all weather conditions (Wan, 2014).To overcome the drawbacks of the conventional methods, an improved method to detect atmospheric water vapour using Global Navigation Satellite System (GNSS) technology has emerged with a high spatiotemporal resolution and low cost.
Tropospheric tomography requires that the area of interest be discretised into many voxels and using those satellite rays crossing the voxels therein (Yao and Zhao, 2016); however, many voxels are not crossed by a signal ray due to the influence of the geometric distribution of the constellation of GNSS satellites as well as the ground-based GNSS receivers (Flores et al., 2000), which leads to a problem of inversion in the tomographic resolution (Bender et al., 2011a).Normally, some constraints, such as the horizontal constraints and vertical constraints, are imposed in tomographic modelling so as to overcome rank deficiency problems (Flores et al., 2000;Bi et al., 2006;Troller et al., 2006;Rohm and Bosy, 2011;Chen and Liu, 2014;Rohm et al., 2014).Ding et al. (2018) proposed a node parameterisation approach using a combination of three meshing techniques to dynamically adjust both the boundary of the tomographic region and the position of nodes.It must be noted that the satellite rays used in the conventional tomography method refer to signals crossing out from the top side of the tomographic body, while the signals penetrating from the side faces are ignored as "unusable" information (Yao et al., 2016).To maximise the usage of the signals penetrating from the side faces of the tomographic body, some studies have been carried out in recent years.A ray-tracing water vapour model has been applied by Rohm and Bosy (2011) to estimate the outer part of the GNSS signal rays using the UNB3m model.Notarpietro et al. (2011) proposed a method to obtain the values of GNSS signals inside the tomography area by extracting the values outside the area, where the values outside the area are estimated with the support of the CIRA-Q wet atmospheric climatologic model or the ECMWF reanalysis data.The water vapour content inside the tomography area is estimated by geometrical linear estimation using an exponential negative function (Benevides et al., 2014).Yao et al. (2016) proposed a method using signals crossing out from the side faces of the tomographic body by introducing the water vapour unit index which was further improved using a more sophisticated model and additional data resources (Zhao and Yao, 2017).Chen and Liu (2016) used numerical weather prediction (NWP) profile data to calculate the water vapour content outside the tomographic sides and then obtained the value of water vapour content inside the tomographic body.However, the methods mentioned above require the availability of specific external information, which restricts the application of those methods.
In this paper, a new tropospheric tomography method is proposed by introducing an auxiliary area into the tomographic process.The introduced auxiliary area can be used to incorporate the contribution of satellite rays crossing the side faces of a tomographic body into the tomographic results.In addition, the COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) radio occultation products are also used for determining the top height of the tomography body as well as the coefficients of a vertical constraint formula.The proposed method is capable of using all satellite rays with a given elevation angle for water vapour reconstruction.This paper is organised as follows: Sect. 2 presents a tropospheric tomography method with the support of an auxiliary area; data processing and the determination of the top height of the tomographic body are covered in Sect.3; in Sect.4, the tomographic results are analysed to validate the proposed method with radiosonde data, and the discussion and conclusions are presented in Sect. 5.
2 An optimal tropospheric tomography method with an auxiliary area

Tomographic principles
Generally, two classes of tropospheric parameters estimated from ground-based GNSS observations are used for the reconstruction of the water vapour field.The first kind of input information is slant wet delay (SWD), which can be used to obtain the atmospheric wet refractivity (Flores et al., 2000;Hirahara, 2000;Skone and Hoyle, 2005;Rohm and Bosy, 2009;Notarpietro et al., 2011).The second kind of input information is slant water vapour (SWV), which can be used to reconstruct the three-dimensional water vapour density field (Champollion et al., 2005;Bi et al., 2006;Chen and Liu, 2014).The corresponding output parameters from those two methods can be mutually interconverted with the help of atmospheric temperature information (Bender et al., 2011a).
In our study, the more sophisticated SWV parameters are selected as input information for a tropospheric tomography experiment.
The SWV can be described as a total water vapour content from satellite to receiver antenna along the propagation path, and can be expressed by where ρ(s) is the water vapour density (unit: g m −3 ) and ds is the distance travelled by a satellite ray.According to tropospheric tomography theory, the research area is discretised into many voxels under the assumption that the water vapour density value of each voxel is unchanged during a given period.Therefore, many satellite rays with different azimuth angles and elevation angles can be used to build the observation equation governing the tomographic modelling.The relationship between SWV and water vapour density can be expressed in the matrix form below: where y represents a column vector with many SWV values and A is a coefficient matrix of distances crossed by satellite rays, while x denotes the unknown parameter which refers to the water vapour density in our study.
Apart from the observation equation, some additional constraint equations are required to solve the rank deficiency issue arising from the inversion algorithm (Flores et al., 2000;Bender et al., 2011b;Rohm and Bosy, 2011).Usually, two different classes of constraints, the horizontal and vertical constraints, are employed in tropospheric tomography (Flores et al., 2000;Bi et al., 2006;Troller et al., 2006;Bender et al., 2001b;Rohm and Bosy, 2011;Chen and Liu, 2014).In our study, a horizontal constraint equation is established based on the assumption that the water vapour in a certain voxel is a weighted mean value of its horizontally near neighbours (Yao et al., 2017).As for the vertical constraint equation, the exponential negative function is introduced to reflect the relationship between two adjacent voxels vertically.The coefficients of the exponential negative function are fitted based on the average water vapour distribution at different altitudes derived from the COSMIC RO "wetPrf" profiles for the period of 2006 to 2016 (see Sect. 3.2).After the various constraints are combined with the observation equation, the conventional tomographic model can be obtained as where H and V represent the coefficient matrices of horizontal and vertical constraints, respectively.

An improved tropospheric tomography method
For most current tropospheric tomography studies, only those signals passing through the entire tomographic area are used, while the signals crossing out from the side faces of the tomographic body are ignored.Such behaviour not only neglects the valuable observations, but also decreases the numbers of satellite rays usable and voxels crossed by rays.As shown in Fig. 1, blue rays P1-P4 which cross the side faces of the tomographic body cannot be used for water vapour tomography based on the conventional methods.
To incorporate the value of rays crossing the side faces of a tomography body, a novel method is proposed by introducing an auxiliary area.The main idea is that, when the tomography area is extended in latitudinal and longitudinal directions over a certain distance, the satellite rays which penetrate from the side faces of the tomographic body are crossing out from the top boundary of the surrounding area (called the auxiliary area).Therefore, the rays crossing out from the side faces of the tomographic body can be used for the reconstruction of the water vapour field.For the original tomographic area, the reconstructed water vapour density using the auxiliary area can be used as the initial water vapour field.The specific steps are presented as follows.
1. Determine the horizontal distance of horizontal extension.As shown in Fig. 2, the blue rectangle marks the research area selected for tomography, while the red rectangle is the auxiliary area.Therefore, the distance between the research area and auxiliary area is first required to be determined.Assuming that one station is located on the border of the research area, as shown by in Fig. 3, the extended distance can be calculated using the following formula: where d is the extended distance, H is the height of the tomography area, and α is the satellite cut-off elevation angle.
2. Reconstruct the water vapour field of the auxiliary area: according to the tomographic model established in Eq. ( 3), the water vapour information for the auxiliary area can be obtained using GNSS measurements with elevation angles greater than α.
3. Select the water vapour density of voxels located in the tomographic area and consider those values initial constraints for the reconstruction of the water vapour field therein.Given this, a novel tomographic model for the proposed method can be established as where I is the unit matrix of the initial constraint equation, while ρ o is the initial water vapour density value of each voxel derived from the tomographic result in Step (2).
4. Singular value decomposition is used to solve Eq. ( 5), and the general inverse matrix can be obtained (Ran and Ge, 1997).In this way, the contribution of satellite rays which penetrate from the side faces of the tomographic body can be incorporated into the water vapour reconstruction process.
3 Data processing and determination of the top boundary

Data collection and processing
In the present study, the data from 13 ground-based stations ( in Fig. 2) derived from the CORS Network in Texas are selected for the period of 11 to 25 May 2015.One radiosonde station (• in Fig. 2) is located in this area, where radiosonde balloons are launched twice per day at 00:00 and 12:00 UTC, respectively.The space-based COSMIC postprocessed "wetPrf" profiles provided by the University Corporation for Atmospheric Research/COSMIC Data Analysis and Archival Center (UCAR/CDAAC) are also used for the period 2006 to 2016 (http://www.cosmic.ucar.edu/,last access: September 2017)."wetPrf" profiles include the latitude and longitude of the perigee point, refractivity, temperature, pressure, and mean sea level altitude from the surface to 40 km (Hudnut et al., 2007).The parameters provided by "wetPrf" profiles are interpolated at altitude intervals of 100 m and calculated by the non-standard 1DVar (onedimensional variational) technique (Hudnut et al., 2007).
The research area is divided as follows: latitude ranges from 32.1 to 33.3 • N in steps of 0.2 • and longitude ranges from 96.5 to 98.3 • W in steps of 0.3 • .The top height of the tomographic body is at an elevation of 10 km with steps of 1 km, which is determined based on the average water vapour variations with altitude from COSMIC "wetPrf" profiles.Consequently, there is a total of 6 × 6 × 10 voxels obtained in our study.The satellite cut-off elevation angle is selected as 10 • ; therefore, the auxiliary area can be determined based on the method proposed in Sect.2.2 with a latitude range from 31.7 to 33.7 • N, a longitude range from 96.2 • W to 98.6 • W, and a total of 10 × 8 × 10 voxels.
The ground-based GPS observations for the experiment period of 15 days are processed with GAMIT/GLOBK (v10.5)(Herring et al., 2010).The zenith total delay (ZTD) parameter and wet delay gradients in the east-west and north-south directions are estimated at intervals of 0.5 and 2 h, respectively.In the GPS data processing, four International GNSS Service (IGS) stations (INEG, NIST, NLIB, and PIE1) with a baseline length larger than 500 km are also used (see Fig. 2) to reduce the strong correlation of tropospheric parameters caused by the similar propagation paths of signals between satellites and stations in the local area (Rocken et al., 1995).The accurate zenith hydrostatic delay (ZHD) is calculated based on the Saastamoinen model (Saastamoinen, 1972) using the interpolated pressure parameters derived from the ECMWF ERA-Interim products with a resolution of 0.125 • × 0.125 • .The zenith wet delay (ZWD) can be obtained by extracting ZHD from ZTD.Therefore, the SWDs of satellite rays with various azimuth and elevation angles are obtained using the wet mapping function (Niell et al., 2001).The weighted mean temperature, which is highly related to the specific area and season, can be obtained by statistical analyses of a large number of radiosonde profiles (Bevis et al., 1992).Here, T m is calculated based on the empirical formula (T m = 70.20 + 0.72T s ) proposed by Bevis et al. (1992) using the interpolated surface temperature (T s ) from ERA-Interim products.Therefore, the SWV can be obtained using the calculated conversion factor with T m (Bevis et al., 1992).

Determination of the top boundary
As described in previous studies, a reasonable top height of the tomographic body is important for tropospheric tomography (Yao andZhao, 2016, 2017).If the height is too high, the estimated water vapour density values are possibly negative as more parameters near the altitude of the top boundary tend to zero.In contrast, a relatively large water vapour value is returned by the tomographic technique as some water vapour contents above the height of the top boundary are ignored.In our study, the top boundary of the research area is determined using the COSMIC RO product which is from the selected 178 "wetPrf" profiles from the period of 2006 to 2016.In addition, the water vapour density changes with altitude are also obtained from radiosonde recordings at station 72249 for the period from 2000 to 2015 to further verify the selected top boundary (Chen andLiu, 2014, 2016;Ye et al., 2016).Here, a principle is introduced to determine the final tomographic height.The principle is that the water vapour density is less than 0.2 g m −3 over the certain height, while the standard deviation (SD) is less than 0.05 g m −3 , which has been used by Yao andZhao (2016, 2017).
Figures 4 and 5 show the vertical distribution of water vapour density with height in a single event (left) and the water vapour distribution at different altitudes and exponential fitting curves for the selected period (right) derived from COSMIC RO and radio-sounding, respectively.It can be seen from Figs. 4 and 5 that the water vapour density generally decreases with the height increases, and the mean water vapour contents for the corresponding periods are both close to zero above a height of 10 km.The calculated water vapour density and SD for COSMIC RO and radiosonde aspects are 0.13/0.040and 0.17/0.044g m −3 at the height of 10 km, respectively, which both satisfy the principle mentioned above.Therefore, the top height of the tomographic body in our study was deemed to be 10 km.
In addition, a negative exponential function is introduced and the coefficients of the exponential model are estimated using the water vapour density from COSMIC "wetPrf" profiles as ρ = 10.5 • e (−0.437 h) , (6) where ρ represents the water vapour density while h refers to the geodetic height (unit of kilometres).This exponen-tial model is used for establishing the vertical constraints between the adjacent layers as described in Sect.2.1.

Result validation and analysis
Based on the method proposed in Sect.2, two schemes are designed to build the tomographic model and evaluate the performance of the various tomographic results.The two methods are presented as follows: -Method 1: only using the signals crossing out from the top side of the tomography area to reconstruct water vapour information based on the tomographic model (e.g.Eq. 3); -Method 2: both the signals penetrating from the top and side faces of a tomography body are used to reconstruct water vapour information based on the tomography model proposed in this paper (e.g.Eq. ( 5).

Utilisation rate of signals and number of voxels crossed by rays
If the satellite rays crossing the side faces of the tomography body are also used, the utilisation rate of signals and the number of voxels crossed by rays are expected to increase; therefore, the number of signals used and the number of voxels crossed by rays are first analysed using data from 13 stations derived from the CORS Network in Texas over a 15day period.Figure 6 shows the number of signals used and the number of voxels crossed by rays, from which it can be seen that the number of signals used and the number of voxels crossed by rays of the proposed method are both larger than those from the conventional method.Table 1 lists statistical results pertaining to the number of signals used as well as the number of voxels crossed by rays for various conditions.The numerical result shows that the utilisation rate of satellite rays is increased by 21.06 %, while the number of voxels crossed by rays is improved by 1.53 %, from 81.39 % to 82.78 %, respectively.

Water vapour comparison with radiosonde data
It has been proven that radiosonde data can provide fairly accurate vertical profiles of tropospheric water vapour; there-   fore, the water vapour profiles derived from radiosonde data are used as a reference to validate the tomographic results from various methods (Niell et al., 2001;Adeyemi and Joerg, 2012).The water vapour profiles for the location of radiosonde stations derived from various tomographic results are compared with that from radiosonde data at 00:00 and 12:00 UTC daily for the experimental period of 15 days.
Here, rms and MAE are introduced to quantify the quality of the tomographic result.Figures 7 and 8 present the average rms and MAE comparisons of tomographic results derived from two methods with radiosonde data during the experimental period, respectively.The rms and MAE comparisons clearly show that the accuracy of water vapour profiles from method 2 is superior to that from method 1. Statistical results for the experimental period (Table 2) show that the mean rms/MAE values of the proposed method are 1.49 and 1.04 g m −3 , respectively, while the values using the conventional method are 1.74 and 1.24 g m −3 , respectively.In addition, water vapour profile comparisons for the specific two epochs at 00:00 UTC, 14 May and 12:00 UTC, 19 May are shown in Fig. 9.Those two times are selected because they correspond to the minimum and maximum rms values for the tested period.It can be seen from Fig. 9 that the water vapour density values derived from various tomographic methods match the radiosonde data at most altitudes.The statistical results show that the rms/MAE values   of method 2 are 0.43 and 2.75 g m −3 , respectively, while the values of method 1 are 0.54 and 2.93 g m −3 , respectively.
To investigate the relationship between water vapour profile error and altitude, a relative error is introduced and expressed as where x ref is a reference value, which refers to the water vapour density derived from the radiosonde data, and x tomo is the reconstructed water vapour information based on different tomographic methods.For the selected period of 15 days, there are two epochs at 00:00 and 12:00 UTC daily.Therefore, a total of 30 sets of data are compared.Figure 10 presents the average water vapour profile comparison derived from different tomographic methods and radiosonde data (left), rms changes with altitude (middle) and relative error changes with altitude (right) for the experimental period.It can be seen that the average water vapour profiles at different altitudes derived from the proposed method better matched that from the radiosonde data, while the rms error and relative error of the proposed method are less than those arising from the use of the conventional method at all layers.Such results validate the superiority of the tomography method with an auxiliary area proposed in this paper.In addition, Fig. 10 shows that the rms error generally decreases with increasing altitude, while the relative error shows the opposite trend.This occurs because the water vapour is mainly concentrated in the lower layers and for the upper layers the water vapour density is very low; therefore, small differences between radiosonde and tomographic water vapour density values can lead to a large relative error.

SWV comparison
To validate the superiority of the proposed method, the data for 4 days are selected to compare the SWV values derived from the tomographic results and GAMIT 10.5.In this experiment, only 12 of 13 GPS stations in the research area participated in the reconstruction of the water vapour field, while the other station (TXSG) is regarded as a test station under the condition that other settings remained unchanged.The SWV values of the TXSG station calculated using the GPS observations are considered to be the true value.The comparison of SWV values derived from GPS-observed data and tomographic results of different methods is carried out with a tomography step of 0.5 h.Table 3 lists the average rms and MAE arising from the use of different tomographic methods and the statistical results for the selected 4-day period.It can be seen that the rms and MAE values of the proposed method (4.14 and 3.27 mm, respectively) are less than those of the traditional method (5.06 and 4.52 mm, respectively), which indicates the high accuracy of the proposed method.

Conclusion
An optimised tomographic method is proposed in conjunction with the concept of an auxiliary area for the inclusion of the signals crossing the side faces of the tomographic body.The proposed tomographic method is validated using the data of 13 GPS stations from the CORS Network in Texas over a 15-day period.The water vapour profile comparison with radiosonde data shows that the mean rms error and MAE of the proposed method are decreased by 14.37 % and 16.13 %, respectively, when compared to the traditional method.The water vapour profile error changes with altitude also show that the rms and relative errors of the proposed method are less than that from the conventional method at different heights, which indicates the superiority of the proposed method.The SWV comparison between the GAMITestimated and tomographic-derived values shows that the rms error and MAE of the proposed method have been improved by 18.18 % and 27.62 %, respectively, when compared to the conventional method, which further confirms its superiority.

Figure 1 .
Figure 1.Schematic diagram of satellite signals crossing the tomography area (a) and their front view (b), where signals P1-P4 are crossed from the side of the tomography area.

Figure 2 .
Figure2.Geographical distribution of CORS stations (black triangles) selected in the tomography area which corresponds to the blue rectangle and the radiosonde station (black circle).The blue rectangle is the original area, while the red rectangle is the auxiliary region.

Figure 3 .
Figure 3.The schematic diagram of the calculating distance of the auxiliary region.

Figure 4 .
Figure 4.The vertical distribution of water vapour density in a single COSMIC RO event (a) and the 3-D water vapour distribution of tomography area and exponential fitting curve from 2006 to 2016 (b); the x axis shows the mean sea level (MSL) altitude.

Figure 5 .
Figure 5.The vertical distribution of water vapour density in a single sounding event (a) and the 3-D water vapour distribution of the tomography area and exponential fitting curve from 2000 to 2015 (b).

Figure 6 .
Figure 6.The number of signals used and the number of voxels crossed by rays for the experimental period.

Figure 7 .
Figure 7.The rms comparison of the tomographic result derived from two methods with radiosonde data during the experimental period.

Figure 8 .
Figure 8.The MAE comparison of the tomographic result derived from two methods with radiosonde data during the experimental period.

Figure 10 .
Figure 10.Water vapour profile comparison derived from different tomographic methods and radiosonde (a), rms changes with altitude (b) and relative error changes with altitude (c) for the experimental period.
Both horizontal and vertical constraints are incorporated into the tomographic model.The top height of the tomography body is determined and the coefficients of the vertical constraint equation are fitted using the average water vapour density from the COSMIC "wetPrf" profiles for the period of 2006 to 2016.The proposed method improves the utilisation rate of signals used (21.06 %) as well as the number of voxels crossed by rays (1.53 %).

Table 1 .
Statistical result of the number of signals used and the number of voxels crossed by rays for the experiment period.

Table 2 .
Statistical information on the tomographic result compared with radiosonde data for the experimental period.

Table 3 .
Statistical results of SWV comparison between different tomographic methods and estimated SWV values for the selected 4 days.