A New Method of GPS Water Vapor Tomography for Maximizing the Use of Signal Rays

: The spatio-temporal distribution of atmospheric water vapor information can be obtained by global positioning system (GPS) water vapor tomography. GPS signal rays pass through the tomographic area from different boundaries because the scope of the research region (latitude, longitude, and altitude) is designated in the process of tomographic modeling, the inﬂuence of the geographic distribution of receivers, and the geometric location of satellite constellations. Traditionally, only signal rays penetrating the entire tomographic area are considered in the computation of water vapor information, whereas those passing through the sides are neglected. Therefore, the accuracy of the tomographic result, especially at the bottom of the area, does not reach its full potential. To solve this problem, this paper proposes a new method that simultaneously considers the discretized tomographic voxels and the troposphere outside the research area as unknown parameters. This method can effectively improve the utilization of existing GPS observations and increase the number of voxels crossed by satellite signals, especially by increasing the proportion of voxels penetrated. A tomographic experiment is implemented using GPS data from the Hong Kong Satellite Positioning Reference Station Network. Compared to the traditional method, the proposed method increases the number of voxels crossed by signal rays and the utilization of the observed data by 15.14% and 19.68% on average, respectively. Numerical results, including comparisons of slant water vapor (SWV), precipitable water vapor (PWV), and water vapor density proﬁle, show that the proposed method is better than traditional methods. In comparison to the water vapor density proﬁle, the root-mean-square error (RMS), mean absolute error (MAE), standard deviation (SD), and bias of the proposed method are 1.39, 1.07, 1.30, and − 0.21 gm − 3 , respectively. For the SWV and PWV comparison, the RMS/MAE of the proposed method are 10.46/8.17 mm and 4.00/3.39 mm, respectively.


Introduction
Although atmospheric water vapor comprises only a small percentage of the atmosphere, it plays a key role in a series of weather phenomena [1,2]. A good understanding of the spatiotemporal changes in water vapor is helpful in improving weather prediction, water resource management,

GPS Water Vapor Estimation
In this study, the SWV derived from GPS-based receivers is exploited to obtain the water vapor distribution information. The SWV can be obtained by the following formula [24,25]: where ρ water is the density of liquid water; R ω = 461J/(kg · K) is the specific gas constants for water vapor; k 2 = 16.48KhPa −1 and k 3 = 3.776 × 10 5 K 2 hPa −1 are constants; and T m is the weighted mean tropospheric temperature, which is calculated by an empirical formula [26]. In order to obtain the SWD, the ZWD and the wet delay gradient parameters should be mapped into the elevation direction and the unmodeled residuals should be taken into account [27,28]. The specific calculation formula is as follows: SWD = m w (e) * [ZWD + cot(e)(G w NS * cos ϕ + G w WE * sin ϕ)] + R where e and ϕ are the satellite elevation and azimuth, respectively; m w is the wet mapping function; R denotes the unmolded residuals; G w NS and G w WE are wet delay gradient parameters in the north-south and east-west directions, respectively; and ZWD is the zenith wet delay, which can be calculated by subtracting the zenith hydrostatic delay from the zenith total delay.

Traditional Tomographic Equation
In water vapor tomography, the tomographic equation consists of two types of equations, one is the observation equation and the other is the constraint equation. The observation equation is the core of water vapor tomography, which is based on the water vapor integral along the signal ray path. The constraint equation contains horizontal and vertical constraints.

Observation Equation
The whole research area and the water vapor density in the discretized tomographic voxel are shown in Figure 1, and the observation equation can be established on the basis of the GPS signal crossing as follows: where superscript q denotes the satellite index; n is the total number of discretized tomographic voxels; s q i represents the length of the satellite signals within the voxel, which can be calculated by the coordinates of two intersections; and x i is the water vapor density of the voxel, i.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 15 In water vapor tomography, the tomographic equation consists of two types of equations, one is the observation equation and the other is the constraint equation. The observation equation is the core of water vapor tomography, which is based on the water vapor integral along the signal ray path. The constraint equation contains horizontal and vertical constraints.

Observation Equation
The whole research area and the water vapor density in the discretized tomographic voxel are shown in Figure 1, and the observation equation can be established on the basis of the GPS signal crossing as follows: where superscript q denotes the satellite index; n is the total number of discretized tomographic voxels; q i s represents the length of the satellite signals within the voxel, which can be calculated by the coordinates of two intersections; and i x is the water vapor density of the voxel, i.

Constraint Equations
The tomographic result only utilizing the observation equation is typically undetermined in the inversion process. The source of this undetermined nature is related to SWV paths not penetrating every voxel and that the SWV observations have similar geometry (i.e., from a GPS satellite to a station). To overcome this inversion problem, the mostly widely-used method is to impose constraint information [29].
In the research area, the water vapor density, i x , within a certain voxel can be represented by the weighted average of its neighbors because of the relatively stable horizontal distribution [30,31]. The specific horizontal constraint equation is as follows: where m is the total number of voxels in the same layer and w denotes the horizontal weighting coefficient, which can be computed by the Gaussian weighting function as follows: where d is the distance between the centers of the two corresponding voxels.

Constraint Equations
The tomographic result only utilizing the observation equation is typically undetermined in the inversion process. The source of this undetermined nature is related to SWV paths not penetrating every voxel and that the SWV observations have similar geometry (i.e., from a GPS satellite to a station). To overcome this inversion problem, the mostly widely-used method is to impose constraint information [29].
In the research area, the water vapor density, x i , within a certain voxel can be represented by the weighted average of its neighbors because of the relatively stable horizontal distribution [30,31]. The specific horizontal constraint equation is as follows: where m is the total number of voxels in the same layer and w denotes the horizontal weighting coefficient, which can be computed by the Gaussian weighting function as follows: where d is the distance between the centers of the two corresponding voxels. The vertical constraint equation is established using the empirical formula of the water vapor density between two adjacent layers: where x j and x j+m denote the water vapor density in the voxel, j and j + m, respectively, and p j+m is the vertical weighting coefficient, which can be obtained as follows: where e is the Euler's number and is approximately equal to 2.718, and H is the water vapor scale height, which is set as 1-2 km.

New Tomographic Equation
As shown in Figure 1, the research area is divided into n voxels and the water vapor density of each voxel is x i (i = 1, 2 · · · n) in the traditional method. To improve the utilization of GPS data, the water vapor density in the four directions outside the research region are added to the observation equation as unknown parameters. Figure 2 displays the four new parameters (x n+1 , x n+2 , x n+3 , and x n+4 for the north, east, south, and west, respectively) in the proposed method from a top view (left) and a side view (right). The boundaries of these four regions are selected based on the location of GPS receivers and the specifics of the tomographic area (see the process in the following experiment). The observation equation of the proposed method is as follows: The vertical constraint equation is established using the empirical formula of the water vapor density between two adjacent layers: where j x and jm x + denote the water vapor density in the voxel, j and jm + , respectively, and jm p + is the vertical weighting coefficient, which can be obtained as follows: where e is the Euler's number and is approximately equal to 2.718, and H is the water vapor scale height, which is set as 1-2 km. As shown in Figure 1, the research area is divided into n voxels and the water vapor density of each voxel is

New Tomographic Equation
This equation is similar to that in the traditional method; only the four parameters in the where A from Equation (8), H from Equation (4), and V from Equation (6) represent the coefficient matrices of the new observation equation, the horizontal constraint equation, and vertical constraint equation, respectively; X = [x 1 , x 2 , · · · , x n , · · · , x n+4 ] is a vector of the water vapor density containing all voxels to be estimated; ∆ 1 , ∆ 2 , ∆ 3 are the noises of the new observation equation, the horizontal constraint equation, and vertical constraint equation, respectively.

Experiment Description
The tomographic experiment is implemented using data from 12 stations (HKLT, HKMW, HKNP, HKOH, HKPC, HKSC, HKSL, HKSL, HKSS, HKST, HKWS, T430 for tomographic modeling and HKKT for SWV comparison) of the Hong Kong Satellite Positioning Reference Station Network (SatRet) for a week (June 1-7, 2017). Figure 3 shows the GPS station distribution, the division of the grids, and radiosonde station 45004. GAMIT (v.10.61) is used to process the GPS data with a sampling interval of 30 s. The wet mapping function and the empirical formula of the weighted mean tropospheric temperature are used to project SWV on the basis of Equations (1) and (2). The tomographic period is 30 min and 48 tomographic experimental solutions can be processed every day.

Experiment Description
The tomographic experiment is implemented using data from 12 stations (HKLT, HKMW, HKNP, HKOH, HKPC, HKSC, HKSL, HKSL, HKSS, HKST, HKWS, T430 for tomographic modeling and HKKT for SWV comparison) of the Hong Kong Satellite Positioning Reference Station Network (SatRet) for a week (June 1-7, 2017). Figure 3 shows the GPS station distribution, the division of the grids, and radiosonde station 45004. GAMIT (v.10.61) is used to process the GPS data with a sampling interval of 30 s. The wet mapping function and the empirical formula of the weighted mean tropospheric temperature are used to project SWV on the basis of Equations (1) and (2). The tomographic period is 30 min and 48 tomographic experimental solutions can be processed every day.   The scope of the tomographic region is (113.87 • , 114.35 • ) for longitude, (22.19 • , 22.54 • ) for latitude, and (0, 8 km) for altitude. The horizontal resolutions in the east-west and north-south directions are 0.06 • and 0.05 • , respectively, and the vertical resolution is 800 m. It should be noted that water vapor density of the outer region (the north, east, south, and west) can be expressed by an unknown parameter (x n+1 , x n+2 , x n+3 , and x n+4 ), provided that the water vapor density in that region changes are relatively stable and as many as possible signal rays from the side face can be included in the region. To find the boundaries of the outer region that meet the above two requirements, the law of water vapor change, the size of the research area, and the distribution of the stations need to be carefully studied. The specific steps are as follows: 1.
To analyze the law of water vapor change. Figure 4 shows the relationship between water vapor density and altitude during the study period, with the red line being the fitted curve (upper for 45005, middle for 59,280, and lower for 59,316, respectively). The stations, 59,280 and 59,316, are located outside the tomographic area and their coordinates are (23.66 • , 113.05 • ) and (23.35 • , 116.66 • ), respectively. We can see that the water vapor concentrates in the lower troposphere, where it also changes more rapidly than above. In order to obtain the function between water vapor density (x) and altitude (h), the polynomial fitting method was used. Equation (10) was estimated by taking the altitude as the independent variable, while the water vapor density was regarded as the independent variable in Equation (11): 2.
To find the suitable altitudinal boundary of the outer region. It can be computed by Equation (10) that the value of the water vapor density near the surface is 26.34 g/m 3 . The range of the water vapor density in this research area is considered as around 0-26.34 g/m 3 . For 10 layers in the altitude direction, if the change of the water vapor density with altitude in the outer region is within 0 to 2.634 g/m 3 , it can be defined as a stable outer region. From Equation (11), the water vapor density of 2.634 g/m 3 corresponds to the height of 5586 m. Considering the convenience of rectification and the error caused by fitting, 5500 m can be selected as the alternative to the suitable altitudinal boundary of the outer region.

3.
To decide the boundaries of the outer region in the east-west and north-south directions. The signal rays of all experimental periods acquired by satellite receivers in the research area are collected and their paths in the tomographic grid are calculated. Then, the signal rays passing through from the side face are extracted and divided to two parts, one of which are signals penetrating from the side face at a height greater than 5500 m. Additionally, the intersections of this part of the signal rays with the top layer, namely the layer of 8 km, are counted. The boundaries of the outer region in the east-west and north-south directions are decided by the distribution of these intersections.

4.
To check the boundaries of the outer region. If the horizontal range of the outer region is within twice the length of a voxel, the reasonable boundaries of the outer regions are achieved.
Otherwise, go back to Step 2, increase the value of the altitudinal boundary, and re-proceed to the next steps.

Experiment Description
The tomographic experiment is implemented using data from 12 stations (HKLT, HKMW, HKNP, HKOH, HKPC, HKSC, HKSL, HKSL, HKSS, HKST, HKWS, T430 for tomographic modeling and HKKT for SWV comparison) of the Hong Kong Satellite Positioning Reference Station Network (SatRet) for a week (June 1-7, 2017). Figure 3 shows the GPS station distribution, the division of the grids, and radiosonde station 45004. GAMIT (v.10.61) is used to process the GPS data with a sampling interval of 30 s. The wet mapping function and the empirical formula of the weighted mean tropospheric temperature are used to project SWV on the basis of Equations (1) and (2). The tomographic period is 30 min and 48 tomographic experimental solutions can be processed every day.   To evaluate the performance of the proposed method, we made tomographic solutions using the traditional method with 8 × 7 × 10 voxels, identified as Scheme #1, and the new method above-mentioned, identified as Scheme #2. In addition, Scheme #3 is to extend the size of the entire tomographic area, that is, the scope is changed to (113.81 • , 114.41 • ) for longitude, (22.14 • , 22.59 • ) for latitude. The new scope, which results in 10 × 9 × 10 voxels in the tomographic modeling, is nearly to the outer boundaries of Scheme #2. For Scheme #3, the horizontal resolution in the east-west and north-south direction are 0.06 and 0.05, respectively, and the vertical resolution is 800 m, which is the same as that in Scheme #1. Since the outer regions in Scheme #2 are included, Scheme #3 can also utilize the increased signal rays of Scheme #2 in the tomographic solutions.

Utilization of Voxels and Signal Rays
To compare the influences of different methods on the number of signal rays used, the experimental results were analyzed. Figure 5 shows the number of signal rays used in different schemes at day of year (DOY) 152, 2017, the upper one for the case of per 30 min, and the lower one for each tomographic solution. It can be seen that only a small part of the signal rays is used in Scheme #1 compared with the total number of signal rays observed. Scheme #2 improves the utilization of signal rays in every epoch and solution, and Scheme #3 has the same effect in improving the usage of signal rays, which illustrates that Scheme #3 can be a good comparison to Scheme #2. Moreover, the average number of signal rays observed, used in different schemes, can be seen in Figure 6. The figure illustrates that Scheme #2 improves the average utilization of signal rays used by 19 The new scope, which results in 10 × 9 × 10 voxels in the tomographic modeling, is nearly to the outer boundaries of Scheme #2. For Scheme #3, the horizontal resolution in the east-west and north-south direction are 0.06 and 0.05, respectively, and the vertical resolution is 800 m, which is the same as that in Scheme #1. Since the outer regions in Scheme #2 are included, Scheme #3 can also utilize the increased signal rays of Scheme #2 in the tomographic solutions.

Utilization of Voxels and Signal Rays
To compare the influences of different methods on the number of signal rays used, the experimental results were analyzed. Figure 5 shows the number of signal rays used in different schemes at day of year (DOY) 152, 2017, the upper one for the case of per 30 min, and the lower one for each tomographic solution. It can be seen that only a small part of the signal rays is used in Scheme #1 compared with the total number of signal rays observed. Scheme #2 improves the utilization of signal rays in every epoch and solution, and Scheme #3 has the same effect in improving the usage of signal rays, which illustrates that Scheme #3 can be a good comparison to Scheme #2. Moreover, the average number of signal rays observed, used in different schemes, can be seen in Figure 6. The figure illustrates that Scheme #2 improves the average utilization of signal rays used by 19.68% compared with Scheme #1. For Scheme #3, the percentage of increment is 22.2%.  On the other hand, the number of voxels passing through by the signals is analyzed. The upper image in Figure 7 shows the number of voxels crossed by the signal rays once per 30 min between Schemes #1, #2, and #3, whereas the lower image shows the situation of every tomographic solution in that day. Evidently, the number of voxels crossed by the signal rays in Scheme #2 is more than that of Scheme #1, both for an epoch or a tomographic solution. For Scheme #3, it has the largest number of voxels crossed by signal rays in most cases, but it does not represent a good performance in penetrating voxels for Scheme #3. This is because the total number of voxels in Scheme #3 is 900, which is much larger than 560 in the other schemes. It may make more sense to focus on the percentage of the voxels crossed by signal rays. On the other hand, the number of voxels passing through by the signals is analyzed. The upper image in Figure 7 shows the number of voxels crossed by the signal rays once per 30 min between Schemes #1, #2, and #3, whereas the lower image shows the situation of every tomographic solution in that day. Evidently, the number of voxels crossed by the signal rays in Scheme #2 is more than that of Scheme #1, both for an epoch or a tomographic solution. For Scheme #3, it has the largest number of voxels crossed by signal rays in most cases, but it does not represent a good performance in penetrating voxels for Scheme #3. This is because the total number of voxels in Scheme #3 is 900, which is much larger than 560 in the other schemes. It may make more sense to focus on the percentage of the voxels crossed by signal rays.   On the other hand, the number of voxels passing through by the signals is analyzed. The upper image in Figure 7 shows the number of voxels crossed by the signal rays once per 30 min between Schemes #1, #2, and #3, whereas the lower image shows the situation of every tomographic solution in that day. Evidently, the number of voxels crossed by the signal rays in Scheme #2 is more than that of Scheme #1, both for an epoch or a tomographic solution. For Scheme #3, it has the largest number of voxels crossed by signal rays in most cases, but it does not represent a good performance in penetrating voxels for Scheme #3. This is because the total number of voxels in Scheme #3 is 900, which is much larger than 560 in the other schemes. It may make more sense to focus on the percentage of the voxels crossed by signal rays.  Moreover, the statistical results, which include the average number and percentage of voxels crossed by signal rays once per day in different schemes for the period of DOY 152-158, 2017, can be seen in Figure 8. The finding shows that the average number of Scheme #1 increases by 15.14% from 65.36% to 80.50% under Scheme #2. Additionally, Scheme #3 has the largest value in the average number of voxels, but its percentage (54.74%) of voxels crossed by signal rays is the smallest. For models with different total voxels, the percentage of voxels passing though by signal rays is more important than the number of that. Moreover, the statistical results, which include the average number and percentage of voxels crossed by signal rays once per day in different schemes for the period of DOY 152-158, 2017, can be seen in Figure 8. The finding shows that the average number of Scheme #1 increases by 15.14% from 65.36% to 80.50% under Scheme #2. Additionally, Scheme #3 has the largest value in the average number of voxels, but its percentage (54.74%) of voxels crossed by signal rays is the smallest. For models with different total voxels, the percentage of voxels passing though by signal rays is more important than the number of that.

PWV Comparison
To assess the accuracy of the tomographic result, the PWV value for the location of the radiosonde station is calculated using the water vapor density of the voxels acquired from the three schemes and compared with PWV derived from the radiosonde data at UTC 0:00 and 12:00. In this experiment, eight tomographic solutions are selected daily to calculate PWV. As shown in Figure 9, the trends of the PWV time series are basically consistent and the results of Scheme #3 agree worst with that from the radiosonde. Compared with Scheme #1, it is clear that Scheme #2 is more consistent with the radiosonde data, indicating that the proposed method can improve the tomographic accuracy. The statistical results are listed in Table 1, which show an advantage of Scheme #2 compared with other schemes in terms of RMS, MAE, and bias.

PWV Comparison
To assess the accuracy of the tomographic result, the PWV value for the location of the radiosonde station is calculated using the water vapor density of the voxels acquired from the three schemes and compared with PWV derived from the radiosonde data at UTC 0:00 and 12:00. In this experiment, eight tomographic solutions are selected daily to calculate PWV. As shown in Figure 9, the trends of the PWV time series are basically consistent and the results of Scheme #3 agree worst with that from the radiosonde. Compared with Scheme #1, it is clear that Scheme #2 is more consistent with the radiosonde data, indicating that the proposed method can improve the tomographic accuracy. The statistical results are listed in Table 1, which show an advantage of Scheme #2 compared with other schemes in terms of RMS, MAE, and bias. Moreover, the statistical results, which include the average number and percentage of voxels crossed by signal rays once per day in different schemes for the period of DOY 152-158, 2017, can be seen in Figure 8. The finding shows that the average number of Scheme #1 increases by 15.14% from 65.36% to 80.50% under Scheme #2. Additionally, Scheme #3 has the largest value in the average number of voxels, but its percentage (54.74%) of voxels crossed by signal rays is the smallest. For models with different total voxels, the percentage of voxels passing though by signal rays is more important than the number of that.

PWV Comparison
To assess the accuracy of the tomographic result, the PWV value for the location of the radiosonde station is calculated using the water vapor density of the voxels acquired from the three schemes and compared with PWV derived from the radiosonde data at UTC 0:00 and 12:00. In this experiment, eight tomographic solutions are selected daily to calculate PWV. As shown in Figure 9, the trends of the PWV time series are basically consistent and the results of Scheme #3 agree worst with that from the radiosonde. Compared with Scheme #1, it is clear that Scheme #2 is more consistent with the radiosonde data, indicating that the proposed method can improve the tomographic accuracy. The statistical results are listed in Table 1, which show an advantage of Scheme #2 compared with other schemes in terms of RMS, MAE, and bias.

Water Vapor Profile Comparison
In general, if two vertical layers are exchanged arbitrarily, then the PWV remains the same, but the vertical distribution of the water vapor changes. Therefore, the PWV comparison mentioned above may not denote a correct tomographic result obtained, although the PWV derived from the different schemes are in good agreement with that from the radiosonde. To validate the accuracy of the vertical water vapor density, the water vapor profiles between the radiosonde and the two schemes are compared, as shown in Figure 10. Two dates are selected because they correspond to the maximum and minimum RMS during the week-long experimental period.

Water Vapor Profile Comparison
In general, if two vertical layers are exchanged arbitrarily, then the PWV remains the same, but the vertical distribution of the water vapor changes. Therefore, the PWV comparison mentioned above may not denote a correct tomographic result obtained, although the PWV derived from the different schemes are in good agreement with that from the radiosonde. To validate the accuracy of the vertical water vapor density, the water vapor profiles between the radiosonde and the two schemes are compared, as shown in Figure 10. Two dates are selected because they correspond to the maximum and minimum RMS during the week-long experimental period.  Figure 10 shows that the tomographic water vapor profiles of the three schemes agree well with those from the radiosonde data. The water vapor density profile of Scheme #2 better matches that from radiosonde than the other two schemes, especially in the low layers. Scheme #3 has the worst performance. Compared with radiosonde, the numerical results of the tomography from Scheme #2 for the two selected epochs (RMS are 1.90 and 0.78 g/m 3 ) are superior to those of Scheme #1 (RMS are 2.98 and 1.32 g/m 3 ) and Scheme #3 (RMS are 3.37 and 1.67 g/m 3 ).   Figure 10 shows that the tomographic water vapor profiles of the three schemes agree well with those from the radiosonde data. The water vapor density profile of Scheme #2 better matches that from radiosonde than the other two schemes, especially in the low layers. Scheme #3 has the worst performance. Compared with radiosonde, the numerical results of the tomography from Scheme #2 for the two selected epochs (RMS are 1.90 and 0.78 g/m 3 ) are superior to those of Scheme #1 (RMS are 2.98 and 1.32 g/m 3 ) and Scheme #3 (RMS are 3.37 and 1.67 g/m 3 ). Table 2 lists the statistical results of the three schemes in all water vapor profile comparisons, including RMS, MAE, bias, and SD. The numerical results reveal that the average RMS/MAE/bias/SD values are 1.97/1.51/−0.33/1.87, 1.39/1.07/−0.21/1.30, and 2.23/1.68/-0.41/2.11 gm−3 for Schemes #1, #2, and #3, respectively. The maximum and minimum values of each statistic are shown in the table. These findings show that the proposed method is more consistent with the radiosonde data than the other two traditional methods. It should be noted that Scheme #3 does worse than Scheme #1, which indicates the total number of voxels and the selection of the entire tomographic area size are important to the tomographic result. To make a further comparison of the relationship between the altitude and the errors in different schemes, the tomographic results (DOY 152-158, 2017) are analyzed. The average RMS error and relative error at different layers for the three schemes are calculated in this period to show the pattern of the vertical water vapor distribution and altitude. Figure 11 gives the RMS and relative error changes with height throughout the experimental period. It can be seen that the RMS, generally, decreases with altitude, whereas the relative error shows an opposite trend. A large relative error appears in the upper layers because the water vapor density is very low in those layers and a small discrepancy between the tomographic result and radiosonde data will result in a large value. In addition, the RMS and relative error of the proposed method (Scheme #2) are generally less than those of the traditional methods (Scheme #1 and #3) for the different layers, which validates the improved nature of the proposed method over the traditional methods. These findings show that the proposed method is more consistent with the radiosonde data than the other two traditional methods. It should be noted that Scheme #3 does worse than Scheme #1, which indicates the total number of voxels and the selection of the entire tomographic area size are important to the tomographic result. To make a further comparison of the relationship between the altitude and the errors in different schemes, the tomographic results (DOY 152-158, 2017) are analyzed. The average RMS error and relative error at different layers for the three schemes are calculated in this period to show the pattern of the vertical water vapor distribution and altitude. Figure 11 gives the RMS and relative error changes with height throughout the experimental period. It can be seen that the RMS, generally, decreases with altitude, whereas the relative error shows an opposite trend. A large relative error appears in the upper layers because the water vapor density is very low in those layers and a small discrepancy between the tomographic result and radiosonde data will result in a large value. In

SWV Comparison
To further evaluate the performance of the proposed method, the SWV of station HKKT is computed using the three schemes and their differences against the GAMIT-estimated SWV are identified. It should be noted that only the SWV penetrated from the top boundary in Scheme #1 were used in this comparison. To better show the results, SWV residuals are grouped into individual elevation bins of 5 • , such that all residuals with an elevation angle between 15 • and 20 • are evaluated as a single unit. The RMS of each elevation bin are calculated and are shown in Figure 12, in which the upper and lower one represent the rainy and rainless scenario, respectively. It is clearly visible that the RMS of SWV residuals reduced as the elevation angle decreased. Colors in the figure indicate that better RMS results can be achieved by Scheme #2 in different weather conditions. Scheme #3 also is worse than Scheme #1 in this comparison. In summary, the residuals of SWV are reduced by the proposed method compared with the other two schemes.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 15 used in this comparison. To better show the results, SWV residuals are grouped into individual elevation bins of 5°, such that all residuals with an elevation angle between 15° and 20° are evaluated as a single unit. The RMS of each elevation bin are calculated and are shown in Figure 12, in which the upper and lower one represent the rainy and rainless scenario, respectively. It is clearly visible that the RMS of SWV residuals reduced as the elevation angle decreased. Colors in the figure indicate that better RMS results can be achieved by Scheme #2 in different weather conditions. Scheme #3 also is worse than Scheme #1 in this comparison. In summary, the residuals of SWV are reduced by the proposed method compared with the other two schemes. To further show the superiority of the proposed method in the SWV comparison, the RMS and MAE of every tomographic solution in DOY 152, 2017 are compared in Figure 13. The red marks acquired by Scheme #2 are smaller than the blue markers obtained by Scheme #1, and also smaller than the green markers obtained by Scheme #3 in the solutions of every period of time, which illustrates the superiority of the proposed method over the traditional methods. Moreover, Table 3 shows the average RMS and MAE calculated from the three schemes for 7 days. A similar situation as that in Figure 13 shows that the RMS and MAE of SWV reduced by the proposed method in the experimental period from 11.74/9.15 mm (Scheme #1) and 12.51/9.87 mm (Scheme #3) to 10.46/8.17 mm, respectively. In summary, it is clear evidence that the proposed method improves the accuracy of the tomographic result. To further show the superiority of the proposed method in the SWV comparison, the RMS and MAE of every tomographic solution in DOY 152, 2017 are compared in Figure 13. The red marks acquired by Scheme #2 are smaller than the blue markers obtained by Scheme #1, and also smaller than the green markers obtained by Scheme #3 in the solutions of every period of time, which illustrates the superiority of the proposed method over the traditional methods. Moreover, Table 3 shows the average RMS and MAE calculated from the three schemes for 7 days. A similar situation as that in Figure 13 shows that the RMS and MAE of SWV reduced by the proposed method in the experimental period from 11.74/9.15 mm (Scheme #1) and 12.51/9.87 mm (Scheme #3) to 10.46/8.17 mm, respectively. In summary, it is clear evidence that the proposed method improves the accuracy of the tomographic result.

Conclusions
The observation equations of GPS water vapor tomography are constructed by signal rays, and an increase in the quantity of signal rays added to the tomographic model contributes to an improved tomographic result. This paper proposed a new method that utilizes the signal rays coming from the top boundary and the side face of the tomographic area. Four outer regions were added to the tomographic equation as unknown parameters, and their boundaries were reasonably determined in the experiment. The traditional method using 8 × 7 × 10 voxels with the selected research area (Scheme #1) and another method using 10 × 9 × 10 voxels with an extended research area (Scheme #3) were utilized as comparisons.
The tropospheric parameters were estimated using GAMIT software based on the Hong Kong Satellite Positioning Reference Station Network. The tomographic experiment was subsequently conducted using the three schemes. Numerical results from DOY 152 -158, 2017 show that the proposed method increases the number of voxels crossed by signal rays (15.14% on average) and improves the utilization of the observed data (19.68% on average) in comparison with the traditional method (Scheme #1). For Scheme #3, the cost of increasing the number of signal rays and the quantity of voxels crossed by rays was to add voxels in total from 560 to 900, which makes the proportion of crossed voxels drop.
Radiosonde data were used as a reference to validate the tomographic result. In a comparison of the PWV time series, the RMS, MAE, and bias values of the proposed method were 4.00, 3.39, and 2.76 mm, respectively, which were smaller than those of the two traditional methods. A comparison

Conclusions
The observation equations of GPS water vapor tomography are constructed by signal rays, and an increase in the quantity of signal rays added to the tomographic model contributes to an improved tomographic result. This paper proposed a new method that utilizes the signal rays coming from the top boundary and the side face of the tomographic area. Four outer regions were added to the tomographic equation as unknown parameters, and their boundaries were reasonably determined in the experiment. The traditional method using 8 × 7 × 10 voxels with the selected research area (Scheme #1) and another method using 10 × 9 × 10 voxels with an extended research area (Scheme #3) were utilized as comparisons.
The tropospheric parameters were estimated using GAMIT software based on the Hong Kong Satellite Positioning Reference Station Network. The tomographic experiment was subsequently conducted using the three schemes. Numerical results from DOY 152 -158, 2017 show that the proposed method increases the number of voxels crossed by signal rays (15.14% on average) and improves the utilization of the observed data (19.68% on average) in comparison with the traditional method (Scheme #1). For Scheme #3, the cost of increasing the number of signal rays and the quantity of voxels crossed by rays was to add voxels in total from 560 to 900, which makes the proportion of crossed voxels drop.
Radiosonde data were used as a reference to validate the tomographic result. In a comparison of the PWV time series, the RMS, MAE, and bias values of the proposed method were 4.00, 3.39, and 2.76 mm, respectively, which were smaller than those of the two traditional methods. A comparison of the water vapor density profiles at different altitudes showed that the profile of the proposed method matches the radiosonde data well. This finding was supported by the statistical results for RMS/MAE/bias/SD, which changed from 1.97/1.51/−0.33/1.87 gm −3 of Scheme #1 and 2.23/1.68/-0.41/2.11 gm −3 of Scheme #3 to 1.39/1.07/−0.21/1.30 gm −3 of Scheme #2. In addition, an SWV comparison using data from station HKKT was conducted using different schemes; the RMS and MAE values of the SWV changed from 11.74/9.15 mm of Scheme #1 to 10.46/8.17 mm of Scheme #2 and 12.51/9.87 mm of Scheme #3. In all comparisons mentioned above, the proposed method achieved the best performance, and Scheme #3 with the extended research area performed worst. This suggests that increasing the percentage of voxels crossed by rays is key to improving tomographic accuracy, rather than simply increasing the number of signal rays.
The experiments were only carried out in Hong Kong at a specific time period, in which the boundaries of the outer regions were suitable. When changes occur in the tomographic region, GPS receivers' location, or resolution of voxels, the boundaries of outer regions should be appropriately adjusted. So, the feasibility and superiority remain to be further verified for various periods, regions, and weather conditions.