LEO Constellation-Augmented Multi-GNSS for 3D Water Vapor Tomography

: Global navigation satellite systems (GNSS) water vapor tomography is an important technique to obtain the three-dimensional distribution of atmospheric water vapor. The rapid development of low Earth orbit (LEO) constellations has led to a richer set of observations, which brings new expectations for water vapor tomography. This paper analyzes the inﬂuence of LEO constellation-augmented multi-GNSS(LCAMG)on the tomography, in terms of ray distribution, tomography accuracy, and horizontal resolution, by simulating LEO constellation data. The results show that after adding 288 LEO satellites to GNSS, the 30-min ray distribution effect of GNSS can be achieved in 10 min, which can effectively shorten the observation time by 66.7%. In the 10-min observation time, the non-repetitive effective observation value of LCAMG is 2.38 times that of GNSS, and the accuracy is 1.27% higher than that of GNSS. Compared with GNSS and the global positioning system (GPS), at a horizontal resolution of 13 × 14, the proportion of empty voxels in LCAMG reduces by 5.22% and 22.53%, respectively.


Introduction
Water vapor is the main driving force for weather and climate change [1]. The capability and feasibility of atmospheric water vapor monitoring, by using a global positioning system (GPS), was first proposed in [2]. Since then, the technique used to obtain atmospheric precipitable water vapor (PWV) has been well studied. However, the obtained PWV only reflects one-dimensional water vapor variation over a GPS ground-based station. To provide the distribution of water vapor at different altitudes, and in the vertical direction, the technology of obtaining the three-dimensional (3D) spatial distribution of water vapor, by using GPS observation data, was first proposed in [3]. This discretizes the study area into finite voxels, and assumes that water vapor density values within voxels are stable at a certain time. Then, the slant wet delay (SWD) that is obtained from the GPS observations can be used to invert the water vapor density values for each voxel.
However, constrained by the spatial geometry relationship between global navigation satellite system (GNSS) satellites and ground-based stations, there are still many voxels that cannot be observed, leading to the singular problem in the design matrix [4]. Generally, it can be solved by adding constraint equations [3][4][5]. In addition, optimizing the voxel distribution [6], using rays passing from the side [7], and fusing GNSS station observations outside the tomographic region [8], were all proven to be effective in reducing the proportion of empty voxels in a tomographic region. A new algorithm for discrete tomography regions into inverted cones also proved to be very effective [9]. Introducing moderateresolution imaging spectroradiometer (MODIS) [10], interferometric synthetic aperture radar (InSAR) [11], and ground-based meteorological observation [12] data into the tomographic model, also helps to solve the above problem. Introducing wet refractivity profiles that are derived from radio occultation data can also considerably improve tomographic solutions [13][14][15]. In addition to the effective methods mentioned above, increasing the number of observations is also effective. The multi-GNSS can produce more observations, which can achieve better accuracy and space-time resolution [16,17]. The potential of multi-GNSS tomography was verified by the triple system tomography experiments of GPS, BDS, and GLONASS [18]. Zhao et al. [19] studied the accuracy and stability of GNSS tomography, and noted that the accuracy improvement did not meet the expectations, in addition to the fact that 50% of the voxels were still unobserved under certain conditions.
The current rapid development of low Earth orbit (LEO) navigation satellites results in a large number of observations that could reduce the number of empty tomographic voxels. Some companies and institutes, such as OneWeb, SpaceX, and Boeing, have announced ambitions for their own LEO constellations. Once all the constellations are completed, the number of available LEO satellites would be in the tens of thousands [20]. Some of these LEO constellations can provide navigation augmentation services, which can contribute to precision point position (PPP) ambiguity solving, shorten convergence time [21,22], and improve the accuracy of the ionospheric model [23,24]. Compared with GNSS, LEO satellites have some of the following characteristics: lower orbits, stronger signals, less interference, longer trajectories across the sky in the same time, and faster geometry changes. These also facilitate the rapid modelling of the atmosphere and the 3D tomography of water vapor. Apart from the aforementioned advantages, LEO satellites require a large number of satellites to cover the entirety of Earth, and have a short lifetime, which are also challenges for the development of LEO satellites. In this paper, a tomographic model was constructed on the basis of the Hong Kong satellite positioning reference station network (SatRef), and two different numbers of LEO constellations were simulated. The performance of different combination schemes was analyzed, in terms of observation ray distribution, tomographic time resolution, accuracy, and spatial resolution. This paper is organized as follows: Section 2 describes the LEO data simulation and the tomography approach. The LEO-augmented GNSS water vapor tomographic experiments and results are described in detail in Section 3. Sections 4 and 5 present the discussion and conclusion, respectively.

LEO Constellation Simulation
In navigation, a number of selected LEO satellites can be equipped with a transmitter to transmit similar navigation signals to ground users, so that they can serve as GNSS satellites, but with much faster geometric change to enhance GNSS capability; this is called LEO constellation-augmented GNSS [25]. Figure 1 shows the comparison of the trajectory of GPS and LEO satellites across the sky within 10 min. It can be seen that LEO has a longer trajectory in the same amount of time. The ability of LEO satellites to broadcast navigation signals brought about the possibility of 3D water vapor tomography. Figure 2 shows that LEO satellites can receive GNSS navigation signals, broadcast GNSS navigation-augmentation signals, and broadcast navigation signals, similarly to GNSS satellites. In addition, LEO satellites orbit at a higher altitude than that of the tropospheric region, so the signal can cross the troposphere. LEO satellite's orbit and clock error accuracy are the key factors to determine the accuracy of LEO water vapor inversion. For orbit accuracy, GNSS-based orbits of satellites in LEO have reached generally very high qualities due to advances in dynamical modeling and carrier phase ambiguity fixing. At present, satellite laser ranging (SLR) verification shows that the orbit determination accuracy of LEO satellites is even better than that of GNSS orbit determination [26,27]. Moreover, the function of inter-satellite link between LEO satellites will also help to improve satellite orbit accuracy [28]. Meng et al. [29] also proposed a method to obtain high medium-long stability LEO clock. Besides, the real observation results from the Luojia-1A satellite also show that the accuracy of the pseudorange and carrier phase measurement at high elevation angles is better than 1.5 m and 1.7 mm, respectively [30]. With the support of these studies, we believe that LEO satellites are expected to be introduced into water vapor tomography in the future.
Precise satellite orbit information is the basis for the tomographic experiment. In this study, the real orbit of GNSS satellites was applied, while the orbit of LEO satellites was generated by Satellite Tool Kit (STK) software [31] since most of the LEO constellation is still under construction. For the number of LEO satellites, 60 and 288 satellites were selected for analysis. The 60-satellite scheme ensured that at least 1 satellite was visible from any location on the globe, while the 288-satellite scheme ensured that at least 4 satellites were globally visible. Their planes were 6 and 12 equidistant orbital planes, respectively, and the orbit inclination was set as 90 • . Orbit altitudes were set as 1000 km, and orbit types were set as polar orbit. In this way, LEO satellite constellation could travel north-south at both poles of Earth, which can cover all of Earth's surface. The specific track design scheme is shown in Table 1. More details of the LEO orbit simulation can be found in [32].

Tomographic Theory and Approach
Since the orbital altitude of LEO satellites is 1000 km, the signal of LEO and GNSS satellites passes through the troposphere when it reaches the ground-based receiver. Hence, tropospheric water vapor information can be derived from LEO and GNSS signals. Slant water vapor (SWV), as an input observation for the tomographic model, could be obtained by the method proposed in [3]. The observations of the tropospheric water vapor tomography are the integration of water vapor in the direction of the ray-path SWV which can be expressed as: where ρ is water vapor density, s denotes the distance of the line of sight in the troposphere, and ds represents the distance of a satellite ray path element.
For the tomographic observation model, the studied region can be divided into several voxels in the horizontal and vertical directions, as shown in Figure 3. The observation rays can be determined by the positions of station and satellites that pass through the discrete voxels. The water vapor density value parameter is considered to be a constant for the selected time period. The slant path is replaced by a straight line. Bender et al. [17] pointed out that when the elevation angle is less than 30 degrees, there will be an error in treating the observation ray as a straight line, but this error is relatively small. Therefore, the observation ray will be replaced by a straight line in this paper. Thus, the total SWV value could be expressed as the sum of the delayed fractions in each voxel that were discrete along the satellite ray path, as follows: where m, n, p are the total number of voxels in the longitude, latitude, and vertical directions, respectively; a ijk is the distance of satellite rays; and x ijk is the unknown water vapor density in voxel (i, j, k). The linear equation of SWV can be described as follows: where y is the derived SWV; x is the unknown state vector that is to be estimated in the voxel; and A is the kernel matrix containing the sub paths of each slant path in each voxel, which could be obtained by calculating the relative positions from the satellite, ground receiver, and tomographic voxels. The satellite position includes the positions of the GNSS and LEO satellites, of which the latter could be obtained from simulated LEO satellite orbital information. However, the observed values of y in the course of this experiment were simulated from reanalysis datasets. Fifth-generation accurate global atmospheric reanalysis (ERA5) data released by the European Centre for Medium-Range Weather Forecasts (ECMWF) were employed for SWV simulation. In this study, ERA5 hourly data of pressure levels from 1979 to the present were employed [33]. The parameters of ERA5 data are shown in Table 2. Water vapor density is calculated for each ERA5, which is 3D grid data, and water vapor density values for each tomographic voxel centroid could be obtained by inverse distance weighted interpolation. Pressure-level data include relative humidity and temperature, which can be used to calculate water vapor density x ERA5 as follows [34]: where RH is relative humidity expressed as percentage; T is temperature (K); and e denotes the saturated water vapor pressure (Pa), which can be calculated by Wexler method [35,36], as follows: where c 0 = −6.0436117 × 10 3 , c 1 = 1.89318833 × 10 1 , c 2 = −2.8238594 × 10 −2 , c 3 = 1.7241129 × 10 −5 , and c 4 = 2.858487 in Equation (5), and c 0 = −5.8653696 × 10 3 , c 1 = 2.22410300 × 10 1 , c 2 = 1.3749042 × 10 −2 , c 3 = −3.4031775 × 10 −5 , c 4 = 2.6967687 × 10 −8 , and c 5 = 6.918651 × 10 −1 in Equation (6). Input value y simu for the tomographic simulation can be calculated by kernel matrix A and water vapor density values x ERA5 , as follows: In order to use uniform data for research, the SWV observations of all GNSS inversions were also obtained by ERA5 data simulation as described above. The value of x ERA5 was regarded as the true water vapor density in voxels without errors. Therefore, y simu did not contain any errors. All observation errors and those in the functional model are reflected in the SWV residuals. If a suitable residual value is constructed, the simulated observations are likely to be close to the true values. In this study, the following method is employed to calculate the residual values of SWV.
Ware et al. [37] investigated the observation with GPS receivers and radiometers at Platteville and Table Mountain, Colorado on 21-23 May, 1996. Bi et al. [38] studied the SWV difference between the water vapor radiometer and GNSS in Xi'an and Nanjing, China, and concluded that the root mean square (RMS) of the difference is less than 4 mm. Bender et al. [39] also compared the data of water vapor radiometer and GNSS in Germany. The SWV differences in different regions are all of the same order of magnitude. Moreover, Ware et al. [37] provided the SWV differential scattering maps between GPS and water vapor radiometer at different elevation angles. Therefore, we quoted the data in the scatter plot by Ware et al. [37]. By vectorizing this scatter plot, we extracted the maximal-difference data corresponding to different elevation angles. Extracted data are shown in Figure 4; the horizontal coordinate is the elevation angle, and the vertical coordinate is the maximal difference corresponding to each elevation angle. On the basis of the above extracted data, a fitting equation was constructed using linear interpolation. The function was employed to fit the relationship between elevation angle and the maximal difference as follows: where ∆ l,max is the maximal difference in the lth ray; ele l is the elevation angle of the lth ray; then the random error ε l of the lth ray could be obtained, and ε l ∼ N(0, ∆ 2 l,max ), ε l will be used as the residual term of the lth SWV.
Lastly, we gave a distribution map of the SWV residual scatter of GNSS and 288 LEO satellite constellations in 30 min. In Figure 5, the abscissa is the satellite elevation angle, and the ordinate is the residual value of SWV. Each point represents an observation ray. Once the residual values were obtained, SWV observations with an error term could be described by the following equation: where ε is the residual term of SWV. After obtaining observed data y simu , water vapor density in the tomographic region could be solved by using Equation (3). Since matrix A was singular, the simultaneous iterations reconstruction technique (SIRT) was used to solve the inverse of Equation (3) [40], as follows: where k denotes the kth iteration; x k j denotes the result of the jth voxel in the kth iteration; r 0 is the relaxation parameter; a i,j is the value of the element for the ith row and jth column of matrix A; A i denotes the ith row of the matrix; x k=0 is the initial density value for water vapor density, which is obtained by averaging the sounding data of the last three days. Relaxation parameter r 0 gives the correction weight for each voxel with respect to the initial voxel value [41]. The number of iterations determines whether the algebraic reconstruction algorithm can achieve a stable result in the end. We set r 0 to 1.2094, and the maximal number of iterations to 150.
The data-processing flow is shown in Figure 6. The GNSS and LEO satellites in the line-of-sight direction, and ground-base stations can build a design matrix for tomographic analysis. ERA5 data are mainly used to provide water vapor density values. GPS and radiometer comparison data provide SWV residuals. Radiosonde data are used to provide initial water vapor density. Combining the above data, water vapor density valuation in voxels can be calculated by using SIRT. Lastly, the ray distribution, and tomographic time resolution, accuracy, and spatial resolution are analyzed in the paper. Then, the ratio of errors in water vapor density before and after tomography wa selected as an accuracy index for assessment, described as follows: Then, the ratio of errors in water vapor density before and after tomography was selected as an accuracy index for assessment, described as follows: where η is the improvement in water vapor density before and after tomography; x ERA5 is water vapor density that was considered to be a true value; x cal is the voxel density of water vapor obtained from tomography; and x k=0 is the initial density value of water vapor density.

Experimental Area and Data Classification
Eighteen GNSS stations of the Hong Kong Satellite Positioning Reference Station Network and a radiosonde station (ID: 45004) in Hong Kong, were selected for the experiment (as shown in Figure 7). Sounding data were used to calculate the initial values and determine the maximal height for tomography. The study area ranged from 22.190   The maximal height of tomography could influence the choice of observations. The sounding station data for 2 years were used to analyze the appropriate maximal height of tomography. The relationship between water vapor density and height, calculated from the sounding stations in Hong Kong from 2016 to 2017, is plotted in Figure 8. The maximal height was about 9600 m; therefore, the maximal altitude of the tomography voxels was set to be 9600 m. The height of each voxel layer was 800 m, which was converted to orthometric height. Therefore, the total number of voxels was 7 × 8 × 12 (7 latitudinal, 8 longitudinal, and 12 vertical).
In order to study the performance of tropospheric tomography combined with the LEO constellation, observational data combinations were classified into five schemes, namely, G, L288, GREC, GRECL60, and GRECL288, where G, R, E, C, L60, and L288 denote the GPS, GLONASS, Galileo, BDS, 60, and 288 satellite LEO constellations, respectively.

Tomographic Observation Distribution
The observations that are collected by a single navigation system are distributed unevenly, and vary with time [42]. Figure 9 shows the distribution of observations received by ground-based GNSS stations, within one epoch. The center point in Figure 9 denotes the zenith direction of the station; 0~360 • clockwise represents the azimuth angle of the satellites, and the length of the arrow represents the cosine of the elevation angle. In total, 47 satellites were observed for the station, including 8 GPS satellites, 6 GLONASS satellites, 7 Galileo satellites, 18 BDS satellites, 7 LEO constellations with 288 satellites, and 1 LEO constellation with 60 satellites. The BeiDou navigation system has more visible satellites in the Asia-Pacific region. The LEO constellation, with 60 satellites only, has one visible satellite, which is related to the purpose for which the constellation was designed. Sixty satellites were primarily designed from a communications perspective, ensuring that they can be seen from anywhere in the world. L288 is comparable to GPS in terms of the number of visible satellites on the same station, although it has a greater total number of satellites. The multisystem combination is richer in the number of observations and more uniform in its distribution. Therefore, the LEO constellation can enrich the data sources for tomography. In traditional tomography, in order to obtain reliable results and solve the ill-posed problem of the tomography equation, it is necessary to assume that the spatial distribution of water vapor does not change within a certain period (usually 30 min). Thus, enough observations can be obtained by extending the observation time. Figure 10 shows the distribution of GPS and L288 visible satellites during the 30-min observation period at HKSC station. From Figure 10a, it is shown that at HKSC station and at 12 o'clock UTC time, six GPS satellites can be seen during the entire 30-min observation period. Figure 10b is the visible satellite distribution of L288. In the L288, 17 LEO satellites can be successively observed during the entire observation period. It can also be seen, from the figure, that some LEO satellites only take about 7 min from entering to leaving the tomographic area. This shows that LEO satellites have the characteristics of rapid geometric movement. Then, GPS satellites have high orbits, slow angular velocity, and slow geometric changes. Even if the observation time is 30 min, the observation rays of the six satellites are still within the tomographic range. This situation will cause GPS satellite rays to always pass through the same voxel repeatedly, and cause many identical equations in the tomographic design equation. This not only increases the ill-posedness of the inversion of the tomographic design equation, but also increases the time to build the tomographic model, and even increases the computational load. This shows that simply increasing the observation time, although increasing the number of collected observation rays, will also bring many invalid repeated observations. The distribution of all the observed GPS rays at a station is shown in Figure 11a. The position of the rays varied with the motion of the satellite, within 30 min. The length of the trajectory of a GPS satellite across the sky within 30 min was finite, leading to the rays only being within a finite number of voxels. In addition, the distribution of the L288 rays is given in Figure 11b. The distribution of these rays was more extensive than that of the GPS, reaching the edge of the voxel due to the rapid motion of the LEO satellites. Furthermore, 30 min was a sufficient amount of time for the LEO satellites to enter and leave the study area. The figure also reflects a denser distribution of rays along the north-south direction than that in the east-west direction, because LEO satellites were designed to be polar satellites.
We have counted all the observed rays, and defined an effective observation ray. If the observation rays are observed by the same receiver, the numbers of passing voxels are exactly the same, and the distances of the rays passing through each voxel are also similar; we consider these observation rays as repeated observation rays. The repeated observation rays are deleted and only one of them is kept, then the remaining observation rays are defined as effective observation rays. Figure 12 shows the heat map of the number of times that the top voxels of five different satellite combinations were passed by the number of effective rays under six observation times. The sampling rate of G, R, E, C is set to 30 s, and the sampling rate of the LEO satellites is set to 15 s. All the sub-pictures in the first row in Figure 12 are the first epoch data. The observation time of each row of the sub-picture is the same, and each column is the same scheme. The figure reflects that increasing the observation time is beneficial to increase the number of effective observation rays. As the observation time increases, the number of times that the top voxel of GRECL288 is passed by effective rays varies greatly, followed by L288. In the whole scheme, the best result appears in (d7). In addition, the effect of the GRECL288 observation for 10 min (d3) is comparable to the effect of GNSS for 30 min (b7). This shows that when GRECL288 achieves the 30-min ray distribution effect of GNSS, the observation time can be shortened by 66.7%. Shortening the observation time is conducive to the rapid modeling of tomography.   Table 3 shows the number of effective rays (effective rays), and the total number of times that the top voxel observations were passed by effective rays (voxel passed). It can be seen from the table that when the observation time is 10 min, the number of effective rays of GRECL288 is 2.38 times that of GREC, and the number of voxel passed is 3.69 times that of GREC. When the observation time is 30 min, the number of effective rays of GRECL288 is 2.77 times that of GREC, and the number of voxel passed is 2.83 times that of GREC. The above experiment shows that appropriate sampling intervals should be used for different constellations. For constellations with slow angular velocity, the sampling interval should be appropriately increased. If the sampling interval is too small, there will be more repeated observations. For fast-moving constellations, the sampling interval should be reduced. Sampling intervals that are too large will cause the effective observation value to be missed.
To test the appropriate sampling rate, two indicators were employed. The first indicator was the number of voxels without observation rays crossing (empty voxels). The other indicator was the number of voxels that were observed by more than two ground stations (two-station voxels). There were relatively few cases where the rays that were observed by different ground stations were parallel, and they mostly crossed in space. The two-station voxel index indicates that voxels are traversed by nonparallel rays. Figure 13 shows the number of empty and two-station voxels for different constellations, such as G, L60, and L288, with sampling rates of 1, 5, 10, 15, 30, 60, 120, and 300 s, respectively, within 30 min. The number of empty voxels for the three satellite systems with different sampling rates is plotted in Figure 13a. L60 had the largest number of empty voxels, followed by G and L288. The sampling rate had little effect on the number of empty voxels for the G scheme because the G scheme had a small trajectory across the sky in a short period of time, and the ray movement was limited in a short period of time, which made it difficult to enter other new voxels. Compared to the G scheme, the variation in empty L60 and L288 voxels with the sampling rate was larger, because LEO satellites can traverse longer sky tracks in a short period of time, and some valid LEO observations are not effectively captured after increasing the length of the sampling interval. Figure 13b shows the two-station voxels for three systems with different sampling intervals. The value of L288 was the largest, followed by those of G and L60. Similarly, the value of G was nearly the same with different sampling rates. The two-station voxels of the L60 and L288 also varied more with the sampling rate than the G scheme. The two-station voxels of L288 were twice as large as those of G. The above indicators show that L60 ray distribution was the weakest when viewed from a single system. The L288 scheme could increase the observed ray distribution and decrease the number of empty voxels. More importantly, the L288 scheme could greatly increase the number of nonparallel observed rays. Therefore, the sampling rates of GNSS and LEO, for all the following experiments, were set to be 300 and 15 s, respectively.
Moreover, the distribution of the five schemes' rays, namely, G, GREC, GRECL60, GRECL288, and L288, at different times within 1 day, is shown in Figure 14; the observation time was 30 min. Figure 14a shows that, for the empty-voxel indicator, G had the largest value and varied the most throughout the day. The variation was smoother for the other combined systems. The value of GRECL288 was the smallest at any time of the day. The performance of L288 was very promising, and although the number of visible satellites was about the same as that for G, the number of empty voxels was smaller than that of the integrated model of GREC. This suggests that, with enough LEO satellites, LEO satellites can obtain the same ray distribution as that of the GNSS constellation. Figure 14b shows the two-station voxels of the five schemes at different times. GRECL288 was the largest at any time, followed by L288, and the smallest was still G. The two-station voxels of L288 were, on average, 37.8 more than those of GREC, indicating that L288 was better than GREC in terms of distribution. Compared with GREC and G, the average empty voxels of GRECL288 were improved by 2.12% and 10.68%, respectively; the average two-stations voxels of GRECL288 were improved by 7.67% and 33.29%, respectively. This suggests that LEO could generate more useful redundant observations than GRECE and G, even in relatively small discrete regions. The increase in two-station voxels was greater than that in empty voxels, indicating that LEO brought more useful cross-rays with the same station distribution.

Tomographic Accuracy Evaluation
Three observation durations were chosen to analyze the precision and time resolution of the tomography, namely, 10, 20, and 30 min, respectively. The accuracy index was calculated according to Equation (11). Figure 15 shows that the accuracy of all five schemes improved with observation duration. The performance of the GRECL288 scheme is the best, while the G scheme is the worst. In addition, the accuracy of all the schemes improves with increasing the observation time, and the most obvious improvement is the L288 scheme. When the observation time was 10 min, the accuracy of the L288 scheme was worse than that of the GREC scheme, while the accuracy of L288 was better than that of GREC when the observation time was 30 min. Meanwhile, LEO was also better than GNSS for tomography accuracy improvement. A side-by-side comparison showed that GRECL288 took less time than GREC did to achieve the same level of accuracy. GRECL288 had an accuracy factor of 9.82% at 10 min of observation time, while GREC had an accuracy factor of 9.64% at 30 min of observation time. When the observation time is 10 min, the accuracy of GRECL288 is 1.27% and 5.24% higher than that of GREC and G, respectively. When the observation duration was 30 min, the accuracy of GRECL288 was improved by 0.75% compared to that of GREC, and 5.83% compared to that of G. However, a certain amount of observation time is beneficial for accumulating valid observation information. However, a shorter observation time is conducive to rapid tomographic modeling. The water vapor density profiles that were extracted from ERA5, the sounding station, and the five schemes within 30 min, are presented in Figure 16, respectively. Data from the sounding station were used as the initial background field for the tomography, and ERA5 data were used for the tomography without errors and could be considered to be true values. The closer the calculated value profile line was to the value of ERA5, the better the calculation was. There were some errors between the initial and true values. All five schemes showed a tendency to approach the true value, and showed good agreement with the true value at the middle and high levels, but a large gap with the true value at the bottom level.
More information on the accuracy at different altitude levels can be derived from Figures 17 and 18. Figure 17 shows the distribution of the difference between the five schemes, and the ERA5 data at different altitudes. The maximal difference values in all five schemes occurred at an altitude of 0.4 km, followed by 1.2 km. This indicates that, despite the increase in the number of observations, the worst accuracy at the bottom still existed. Comparing the five schemes, the G scheme showed the sharpest fluctuations in horizontal distribution, while the GRECL288 and L288 schemes showed significantly flatter fluctuations. This indicates that the addition of LEO satellites made the tomography results more stable, especially those of the L288 scheme, which showed better stability, even without GNSS data.  Similarly, the RMS and the mean absolute error (MAE) values of the difference between the five schemes and ERA5 data, are given in Figure 18. Among the five schemes, GRECL288 presented the best accuracy. Similarly, the maximal error occurred at the bottom layer, where the distribution of observations entirely depended on the distribution of ground stations. Increasing the number of ground stations, or introducing additional observations, can fundamentally solve the problem of the poor accuracy of the tomographic bottom.

Horizontal Tomography Resolution
Nilsson and Gradinarsky [43], and Chen and Liu [44], analyzed the influence on the size of the voxel for tomography. Manning et al. [45] also indicated that one of the main limitations to solve the horizontal domain of the storm system is horizontal resolution, which is constrained by the GPS network resolution. To investigate the effect of LEO constellation on the horizontal resolution, eight different horizontal resolution strategies were proposed. As the initial grid resolution, 7 × 8 was selected. Then, the latitudinal direction was fixed, and the number of longitude grids was changed, followed by fixing the longitudinal direction and changing the number of latitudinal grids. The detailed settings are shown in Table 4. The percentage of overall empty and two-station voxels, and precision factor η, were also used for analysis.  Figure 19 shows the calculation for the five schemes of the percentage of empty to total voxels, the percentage of two-station to total voxels, and accuracy factor values. Figure 19a shows the empty-voxel values of all five constellation combination schemes increasing with the increase in horizontal resolution. However, GRECL288 had the smallest empty voxels, followed by L288, while the largest were in G. The comparison shows that the slope of the data change of the G scheme is also the largest. This shows that when the observation data is insufficient, the proportion of empty voxels to total voxels will be accelerated when the horizontal grid division becomes dense. Meanwhile, this trend was much flatter with a larger number of satellites. Figure 19b shows two-station voxels for different constellations at different horizontal resolutions. The results still showed that GRECL288 performed the best. Comparing the data of empty voxels and two-station voxels, the gap between the empty-voxel data is relatively close. The two-station voxels data gap is relatively sparse. This shows that LEO has a greater impact on two-station voxels indicators, and the growth of two-station voxels indicates the growth of non-parallel observations. Figure 19c shows the accuracy-factor values for different constellations at different resolutions. Compared with the empty and two-station voxels, the accuracy factor also responded to similar findings as described above. Moreover, GRECL288 had the smallest percentage of empty voxels and the largest percentage of two-station voxels at different horizontal resolutions. The addition of multiple systems and LEO constellations facilitated the improvement of the horizontal resolution of the layer analysis, without the need for additional ground stations. Compared with the G scheme, the multisystem had a higher stability in ray distribution with increased horizontal resolution, and a smoother transition between empty voxels and two-station voxels. L288 performed better than GREC and worse than GRECL288. Even without GNSS, L288 presented the best performance. We calculated the difference between GRECL288 and GREC and G under all grid divisions, and counted the differences of the maximum, minimum, and average values. The data are shown in Table 5. At a horizontal resolution of 13 × 14, compared with GREC and G, GRECL288 improves the empty voxels by 5.22% and 22.53%, two-station voxels by 11.03% and 49.36%, and accuracy factor η by 2.06% and 6.32%, respectively. Interestingly, the largest difference appears in the 13 × 14 grid division scheme. This shows that when the horizontal grid is more densely divided, the GRECL288 can effectively control the increase in empty voxels, and stabilize the decrease in the two-station voxels, so as to effectively restrain the influence of the increased horizontal resolution on the observation data.

Discussion
By comparing five different constellation combinations, LEO has improved the time resolution, tomographic accuracy, and horizontal resolution of water vapor tomography. Especially after excluding the repeated observation rays, the tomographic scheme with LEO has great potential in shortening the time of tomographic modeling. Within the same observation period, the number of non-repetitive observations is also doubled. This is an effect that cannot be achieved by increasing the number of GNSS satellites. In addition, the addition of the LEO satellite tomography scheme can also reduce the proportion of voxels without observations in the tomography, and increase the proportion of voxels that are penetrated by non-parallel crossing rays. In terms of tomography accuracy, the scheme with LEO has improved the tomography accuracy, but the improvement of water vapor tomography accuracy is not as large as expected. This conclusion also echoes some previous studies [18,19]. The reason for this is because of the 30-min observation, so there are already many observation rays in the tomographic area. Another reason that affects the accuracy of tomography is that even if observation satellites are added, there are still voxels without observations in the bottom layer of the tomography. The loss of accuracy in this area affects the overall tomographic accuracy. However, it can also be seen from our experiments that when the observation time is shortened to 10 min, the accuracy of the scheme with added LEO satellites is higher. From the perspective of horizontal resolution, when the horizontal grid is more densely divided, the scheme of adding LEO satellites can effectively suppress the increase in the proportion of empty voxels in the tomographic area. Of course, the biggest key to LEO adding water vapor tomography is the accuracy of the LEO calculation of SWV. Since there are no real data at present, many studies still use simulation data. However, the current studies have given more positive results. The main purpose of this paper is also to quantitatively analyze the changes brought about by the addition of LEO to water vapor tomography. In the future, obtaining the SWV value of LEO inversion with high precision is still a problem that needs to be paid attention to and solved in the future.

Conclusions
In this paper, LEO constellation-augmented multi-GNSS for 3D water vapor tomography was investigated by simulated LEO constellation data. Three aspects of the study were proposed, namely, observed ray distribution, tomographic water vapor accuracy, and horizontal resolution. Then, 60 and 128 LEO constellations were simulated for tomographic water vapor experiments.
From the tomographic observation value distribution experiment, the scheme with LEO satellites showed a more uniform distribution of observation rays. The GRECL288 scheme could shorten the observation time by 66.7% when it reaches the 30-min ray distribution effect of the GREC scheme. In the 10-min and 30-min experiments, the number of non-repetitive effective observation rays of GRECL288 was 2.38 times and 2.77 times that of GREC, respectively.
From the experiment with different sampling rates, it could be concluded that a large amount of data both increases the difficulty of designing the matrix calculations and makes tomographic results unstable. However, for the fast angular rate of LEO, a low sampling rate would miss some valid observations. Hence, appropriate sampling rates of multi-GNSS and LEO were analyzed, and set to be 300 and 15 s, respectively. In addition, L60 was less effective than L288 for the smaller number of satellites, so a higher number of LEO satellites could bring better advantages.
From the results of 1 day, compared to GREC and G, the average empty voxels of GRECL288 were improved by 2.12% and 10.68%, respectively; the average two-station voxels of GRECL288 were improved by 7.67% and 33.29%, respectively.
From the experiment of accuracy comparison, GRECL288 had the highest accuracy, with an improvement of 0.75% and 5.83% compared with that of GREC and G, respectively. Otherwise, the integrated strategy with adding LEO constellation requires less observation time to achieve the same level of accuracy. In this respect, introducing the LEO constellation enhances the time resolution of tomography.
The horizontal resolution experiment shows that when the horizontal resolution increases, adding LEO constellation can effectively suppress the increase in the number of empty voxels, the decrease in the number of two-station voxels, and the decrease in accuracy. At a horizontal resolution of 13 × 14, compared to GREC and G, the empty-voxel ratio of GRECL288 increased by 5.22% and 22.53%, respectively; the two-station voxel ratio of GRECL288 increased by 11.03% and 49.36%, respectively; and the accuracy factor η increased by 2.06% and 6.32%, respectively.
LEO constellations can bring obvious benefits to water vapor tomography. This result provided a new solution for the improvement of spatial and temporal resolution, and accuracy, in water vapor tomography. However, the above improvements are still affected by the distribution of ground stations. As the density of ground stations increases, integrated with LEO constellation data, the water vapor stratification yields even more promising results.
Author Contributions: This paper was conceptualized by X.Z. and S.X.; LEO constellations were designed by F.M.; software was written by S.X.; validation was done by X.R.; data were collated by J.C.; original draft was written by S.X. All authors have read and agreed to the published version of the manuscript.