Tropospheric Wet Refractivity Tomography Based on the BeiDou Satellite System

Water vapor plays an important role in atmospheric transmission, the evolution of weather systems, the radiation budget of the climate system, and global climate change (Liu et al., 2005). Global Navigation Satellite System (GNSS) technology provides a new means for mapping the dynamics of atmospheric water vapor. With its features of real-time continuity, all-weather detection and high precision, it has become the least costly and most convenient means for water vapor detection, and is a key part of the broader global atmospheric observation system (Gutman et al., 2004).

Although the precision of precipitable water vapor (PWV) sensed by GNSS receivers can reach 1-2 mm (Bevis et al., 1992, 1995; Rocken et al., 1993, 1997), unfortunately it is hard to achieve tropospheric water vapor tomography with high precision using only slant wet delays (SWD) (Hirahara, 2000, Flores et al., 2000; Nilsson and Gradinarsky, 2006; Bender et al., 2009; Rohm and Bosy, 2009; Rohm and Bosy, 2011). The primary cause is that the number of in-view satellites is limited, the GNSS signal spatial density between satellites and receivers is usually insufficient, and the effective SWDs are not enough to obtain a stable tropospheric water vapor field. Experts have tried many improvements to solve the above deficiency of ground-based GNSS technology; for example, by introducing horizontal and vertical constraint functions to make tomography equations solvable, or by using iteration methods such as algebraic reconstruction techniques (ART) instead of the least square method (LSM) to avoid the problem of matrix inversion (Bender et al., 2011a).

Up until the end of October 2012, the Chinese BeiDou navigation satellite system (BDS) had launched 16 satellites, and then started providing formal operation services for China and surrounding regions from the beginning of 2013. The additional SWD signals from the BeiDou system may improve the weak spatial structure between satellites and stations based purely on the existing GPS system. Previous research has focused on tropospheric water vapor tomography based on GPS, GLONASS and GALILEO systems (Bender et al., 2011b), but so far there have been no published studies on atmospheric water vapor tomography based on the BeiDou system alone.

Traditional water vapor tomography uses LSM to solve the tomography equations. The disadvantages of LSM (Flores et al., 2000; Nilsson and Gradinarsky, 2006; Perler et al., 2011) include the difficulty of data organization, the inverse problem of a huge matrix, and the oscillation of the results because of the ill-posed nature of tomography equations. ART, based on an iterative method for water vapor tomography, can overcome the above shortcomings of LSM. In this context, the present paper examines tropospheric wet refractivity tomography based on the BeiDou system for the first time, using multiplicative algebraic reconstruction techniques (MART) and simulation data from the Shenzhen and Hongkong (SH) GNSS network. The intention is to answer the following two questions: (2) Can the precision of the tropospheric wet refractivity field be improved when introducing additional SWDs from the BeiDou system? (3) Can the tropospheric wet refractivity field be obtained with high precision based purely on the BeiDou system when GPS becomes unavailable due to certain factors?

Station data were selected from the SH GNSS network with a latitude range of 22^°12'N to 22^°36'N and a longitude range of 113^°53'E to 114^°17'E, as shown in Fig. 1. In total, there are 15 stations in this network, including four (NANS, SZYL, JIAN and SHTJ) in Shenzhen and 11 in Hong Kong, covering an area of approximately 40× 40 km.

The simulation data are based on precise ephemeris, and with the SH GNSS network under an uneventful atmospheric condition (to be termed "standard atmosphere"). The simulation time is 0000-1700 UTC 14 August 2012. The precise ephemeris of GPS satellites was downloaded from the IGS public service network (http://igscb.jpl.nasa.gov; retrieved 19 December 2011), and then interpolated to gain the satellite coordinates with a sampling rate of 30 seconds according to a Lagrangian interpolation method (Ma et al., 2011). The precise ephemeris of the BeiDou satellites was obtained from Wuhan University (Shi et al., 2012). For BeiDou satellite system, there are five satellites on geosynchronous earth orbit (GEO), 5 satellites on inclined geosynchronous satellite orbit (IGSO), and 4 satellites on medium earth orbit (MEO). Satellites are named after the respective orbit they are on and the sequence they were launched (MEO1, GEO1, GEO3, IGSO1, GEO4, IGSO2, IGSO3, IGSO4, IGSO5). Signals from above nine satellites are used in the paper.

Figure 1.  Geographic distribution of the SH GNSS network (approximate size: 40× 40 km).

The grid model for tomography is as follows: grid center coordinates (B=22.4, L=114.1, H=18.105; B, L and H stand for latitude, longitude, and height respectively) act as the coordinate origin of a unified station-origin coordinate system. There are four voxels in the east-west direction, with 10 km per voxel. The same is the case in the north-south direction. There are 16 voxels in each layer. The height of the tropopause is set to 10 km and 500 m per layer. Thus, in total there are 20 layers in the vertical direction and 320 voxels altogether.

The simulation process takes place as follows: (2) the true values of wet refractivity for all voxels and SWDs are calculated from standard atmosphere; (3) 50% system error is uniformly added to the true value of the wet refractivity field as an initial value for iteration; (4) to be more close to the real situation, certain error with multipath characteristics can be added to the SWD true value. The actual situation is that the lower elevation angle, the greater the SWD error. During data processing, σ _ZWD means the error of the ZWD and is set to 3% of the ZWD, and the error of a SWD with elevation angle α is set to σ _ZWD sinα^-1.

Figure 2.  Convergences of (a) δ, (b) σ, (c) Bias and (d) RMS with iterations using MART based purely on BDS, purely on GPS and the BDS+GPS mode. The simulation time is 0800-0830 UTC 14 August 2012; 100 iterations and "standard atmosphere" are used.

We propose two improvements for MART (Bender et al., 2011a): one is that the Gaussian constraint is realized to update the value of the voxels not penetrated by any SWD; the other is that two new parameters of Bias and RMS, calculated from the result of all the voxels, are adopted to evaluate the precision of the tomography result of the entire gridded space, whereas the existing δ and σ parameters, shown in Eqs. (2) and (3), can only evaluate the precision of the voxels penetrated by SWDs. In the above equations, m denotes slant wet delay; I denotes the total number of slant wet delay; and k denotes the k th iteration.

The Bias and RMS parameters are calculated as shown in Eqs. (4) and (5): In the above equations, N_w denotes voxel wet refractivity; n denotes the total number of voxels; and k denotes the k th iteration.

The simulation results, from 0800 to 0830 UTC 14 August 2012, respectively based on BDS, GPS and BDS+GPS, are shown in Fig. 2. As can be seen, the δ, σ and Bias parameters converge rapidly with iterations and stay steady under three modes, yet the RMS parameter keeps a weak divergent trend after initial convergence. Figure 2a shows that there is little difference between δ values under three modes; Fig. 2b shows the σ value under the BDS mode is the minimum in the whole iteration procedure; and Figs. 2c and d show the Bias and RMS values under the BDS mode is the minimum after several initial iterations. Figure 2 shows the number of necessary iterations for MART can be about 20, which is inconsistent with the conclusion of 100-200 iterations given by (Bender et al., 2011a). Too many iterations may result in numerical artifacts and can completely destroy former good results. A small variance of σ does not always indicate a good reconstruction. In several cases the result became worse with an increasing number of iterations, even though σ decreased, as indicated in Fig. 2.

Information from the simulation such as the number of valid SWDs penetrating out of the top of the grid and the number of voxels penetrated by SWDs in the three modes is shown in Table 1. In Table 1, Bias (0-10 km) means the Bias calculated from the results of the total voxels, and Bias (0-5 km) means that of the partial voxels with a layer height from 0 to 5 km; the meanings of RMS (0-10 km) and RMS (0-5 km) are equivalent.

Table 1 shows that evaluating the tomography results based only on δ and σ parameters (Bender et al., 2011a; Notarpietro et al., 2011; Wang and Wang, 2011) is imperfect. The δ and σ decrease with an increasing number of iterations, the Bias changes little after several initial iterations, while RMS (0-10 km) or RMS (0-5 km) diverge slightly with an increasing number of iterations; for example, RMS (0-10 km) increases from 3.52 mm km^-1 at the 20 th iteration to 3.92 mm km^-1 at the 100 th iteration in the BDS mode. So, the number of appropriate iterations should be chosen by not only the δ and σ parameters, but also the Bias and RMS parameters, and not be more than 20 times. It is exciting that the precision of the tomography results of the lower troposphere (0-5 km) is also good in all three modes.

The above analysis is based on a simulation time of 0800-8300 UTC 14 August 2012 and standard atmosphere. In order to make the analysis convincing, Fig.3 gives comparisons of the δ, σ, Bias and RMS parameters at the 100 th iteration between the three modes, based on nearly a whole day's simulation time of 0000-1700 UTC 14 August 2012. A sample is taken every 30 min, and there are 34 samples in Fig. 3. Figure 3 shows that no matter whether based on δ and σ or Bias and RMS, in most cases the BDS mode can achieve the best tomography result. Among them, the σ value in the BDS mode is significantly less than that in the GPS and BDS +GPS mode.

In Figs.2 and 3, standard atmosphere is used to calculate the true wet refractivity field. In order to study the ability of MART for tropospheric water vapor tomography in bad weather, a dynamic atmosphere is also simulated, as shown in Fig.4. Based on standard atmosphere, wet refractivity values for six voxels are changed with 50% amplitude to simulate the dynamic atmosphere, and they are: (-20000, -20000)-(-10000, -10000) from the 5 th to 7 th layer, and (10000, 0)-(20000, 10000) from the 7 th to 9 th layer. Figure 4 gives the true value, initial value and iteration value for profiles above (-20000, -20000)-(-10000, -10000) and (10000, 0)-(20000, 10000) in the BDS mode, and it shows that MART has detected the voxels with abrupt change at the time of 0800-0830 UTC 14 August 2012. The case for GPS and BDS+GPS is similar.

Figure 5 provides comparisons of the δ , σ , Bias and RMS parameters at the 100 th iteration between standard and dynamic atmosphere in the three modes. The simulation time is 0000-1700 UTC 14 August 2012. As shown, MART has the ability to achieve tropospheric water vapor tomography, in bad or stable weather conditions.

The results of the above analysis show that, for the SH GNSS network, and for the simulation time of 0000-1700 UTC 14 August 2012, the tropospheric wet refractivity tomography based purely on the BDS mode (only nine satellites used) can achieve results superior to that based purely on the GPS mode. Furthermore, the BDS+GPS mode shows no obvious improvement in terms of the spatial structure of the existing SWDs from the GPS system.

Figure 3.  The change of (a) δ, (b) σ, (c) Bias and (d) RMS with time under the different modes. The simulation time is 0000-1700 UTC 14 August 2012; 100 iterations and "standard atmosphere" are used.

Figure 4.  The simulated dynamic atmosphere and the profiles of wet refractivity for voxels with abrupt change based purely on BDS mode using MART. There are six voxels with abrupt change: (-20 000, -20 000)-(-10 000, -10 000) from the 5 th to 7 th layer, and (10 000, 0)-(20 000, 10 000) from the 7 th to 9 th layer. The simulation time is 0800-0830 UTC 14 August 2012; 100 iterations are used.

Figure 5.  Comparisons of (a) δ, (b) σ, (c) Bias and (d) RMS parameters between the standard and dynamic atmosphere in the three modes. The simulation time is 0000-1700 UTC 14 August 2012; 100 iterations are used. D.A means dynamic atmosphere, S.A means standard atmosphere.

Figure 6.  Distribution of effective SWD signals of BDS and GPS with elevation angle.

Following the above analysis, we ask why the BDS mode achieves the best results for the SH GNSS network compared with the GPS and BDS+GPS modes. Figure 6 shows the distribution of effective SWD signals with the elevation angle, demonstrating that the SWD signals with high elevation angle in BDS are considerably more than that in GPS during the simulation time of 0800-0830 UTC 14 August 2012. The higher the elevation angle, the smaller the error of the SWD, the δ and σ parameters are calculated from SWDs, which may explain why the σ parameter in the BDS mode is significantly less than that in the GPS and BDS+GPS modes.

Although the BDS+GPS mode includes SWDs nearly double those in the GPS or BDS modes (Table 1), there is no doubt that the BDS+GPS mode has not obviously improved the spatial structure of the existing SWDs from the GPS system for the SH GNSS network. This may be because there are too many ground stations for the SH area covering only about 40× 40 km; the existing GPS mode provides enough SWDs from the above direction and is lacking SWDs from the horizontal direction, whereas the BDS mode can only provide redundant SWDs from the above direction and has no contribution in providing SWDs from the horizontal direction. This also explains why the precision of tomography based on the BDS mode is comparable to that based on the GPS mode and their SWDs are similar rather than complementary.

Figure 7.  Convergences of (a) δ, (b) σ, (c) Bias and (d) RMS with iterations using MART based purely on BDS, purely on GPS, and the BDS+GPS mode. The simulation time is 0800-0830 UTC 14 August 2012; 100 iterations and the SH GNSS network with only four stations (NANS, HKST, HKSL and HKOH) and "standard atmosphere" are used.

The used nine BeiDou satellites include one MEO, three GEOs, and five IGSOs. The relative positions of GEOs to the GNSS stations remain unchanged, and the high utilization of GEOs and IGSOs allows the BeiDou satellites to keep high elevation angles and be visible nearly all day long, which is favorable for Chinese ground-based GNSS technology in tropospheric water vapor tomography. Although there are 24 GPS satellites in service, their utilization is far below the BeiDou satellites for China.

An experiment using the SH GNSS network with inadequate stations was conducted. Four stations (NANS, HKST, HKSL and HKOH) were selected. The simulation results, from 0800 to 0830 UTC 14 August 2012, respectively based on the BDS, GPS and BDS+GPS systems, are shown in Fig. 7. Statistical results for the sparse SH GNSS network in the three modes are shown in Table 2. Table 2 shows that the number of voxels penetrated by SWDs in the BDS mode, with relative concentration of SWDs with high elevation angle, is far smaller than that in the GPS mode. Figures 7a and b show that the δ and σ values in the BDS mode are minimum, but Fig. 7d shows that the SWD structure in the BDS mode is highly ill-posed. Figure 7 and Table 2 show that the BDS+GPS mode also cannot improve the precision of tropospheric water vapor tomography for the GNSS network with very inadequate stations. Good precision of water vapor tomography is expected as long as there are adequate and independent observations, when the spatial distribution of SWDs is decent, and regardless of whether GPS or BDS is used, or both are used together.

The present study used BeiDou satellite signals for the first time to reconstruct the tropospheric wet refractivity field with a simulation method. Multiplicative algebraic reconstruction techniques were used to solve the tomography equations. Compared with previous research, we added two parameters (Bias and RMS), computed from the tomography results of the total voxels, and they were able to evaluate the precision of the wet refractivity field more comprehensively and directly than mean deviation (δ) and variance (σ) parameters computed from SWD. The main conclusions are as follows:

(2) Based on a simulation time of the period 0000-1700 UTC 14 August 2012, the precision of tomography using the BeiDou system (only 9 satellites used) is superior to that based on the GPS system in most cases for the SH GNSS network. This may have great significance for achieving independent tropospheric water vapor tomography operations in China based on ground-based GNSS technology. With the completion of spatial network for the BeiDou satellite system in 2020, the spatial structure of SWD signals purely from the BeiDou system may be good enough to achieve a 4D tropospheric water vapor field.

(3) For the SH GNSS network, either with adequate stations or not, the BDS+GPS mode does not significantly improve the spatial structure between satellites and stations, and the precision of tropospheric wet refractivity tomography does not improve. The station density of the network has a strong effect on the precision of water vapor tomography; it is necessary to maintain adequate stations in a network in order to obtain the water vapor field with high precision. What the most reasonable station density of the network is will be investigated in future work.